利用ICESat、Cryosat-2、ICESat-2三颗测高卫星数据提取2003-2022年青藏高原湖泊水位变化序列，并以青海湖为例对比分析Cryosat-2 LRM三种重跟踪算法(Ocean-CFI、UCL Land-ice、OCOG)，结果表明，OCOG算法最优。利用青海湖、纳木错、美玛错、鲁玛江冬错的水文站数据，对通过卫星测高数据提取的湖泊水位变化时间序列精度进行评价，结果发现，本文提取的湖泊水位变化值精度较高，青海湖水位变化时间序列与下社站实测结果差异的均方根仅为0.092m，纳木错湖水位变化时间序列与实测结果的相关系数达0.989。相同时段内不同测量技术的精度分析表明，Cryosat-2 SARIn测量模式的精度略高于ICESat-2。提取青藏高原106个湖泊2003-2022年水位序列，并分析其变化趋势，结果表明，106个湖泊中19个湖泊水位呈下降趋势，83个湖泊水位呈上升趋势。
关键词：多源遥感数据;卫星测高;青藏高原;湖泊水位变化
论文《联合多源卫星测高数据获取青藏高原湖泊水位变化》发表在《大地测量与地球动力学》，版权归《大地测量与地球动力学》所有。本文来自网络平台，仅供参考。
引言
青藏高原被誉为“亚洲水塔”，湖泊密布，面积广阔，冰雪储量巨大，滋养了亚洲多条河流。该区域人口密度低，受人类活动影响小，为探究湖泊演变机制提供了理想环境。随着卫星遥感技术的发展，学者们将其广泛运用于湖泊动态监测，发现多源测高卫星数据能显著增加监测湖泊的数量并延长时间序列的长度[2]。这些卫星凭借良好的轨道连续性和更替重叠期，促进了数据融合，实现了对湖泊的长期监测[7-10]。
传统雷达高度计测高卫星虽时间跨度长、重访周期短，但地面足印大，易受陆地干扰，难以精确监测小湖泊。而ICESat及ICESat-2卫星，尤其是后者，以更小的地面足印(约17m)，有效解决了小于100km²湖泊的水位监测难题[11-12]。但激光测高卫星普遍面临任务不连续、监测时间跨度短的挑战，限制了长期动态监测。相比之下，Cryosat-2的SARIn模式在青藏高原复杂地形中表现优异，更适用于该地区[13]。
鉴于传统卫星高时间分辨率、低空间分辨率与新一代卫星低时间分辨率、高空间分辨率的互补性，本文结合ICESat、Cryosat-2等多源卫星数据，利用其高精度、连续观测的优势，成功提取2003-2022年青藏高原湖泊的水位信息，并分析其变化趋势。
1 数据介绍
本文使用ICESat、Cryosat-2与ICESat-2卫星测高数据提取水位信息，选择具有明显变化特征的湖泊，包括内流区羌塘盆地的色林错、纳木错、多尔索洞错、赤步张错，黄河源流域的青海湖及雅鲁藏布江流域的羊卓雍错。
使用的卫星测高数据主要包括ICESat卫星的GLAH14产品、Cryosat-2的L2 GDR产品及ICESat-2的ATL13产品，其中GLAH14产品提供全球地表测高数据，包括采样点的经度、纬度、时间、大地高(参考椭球为T/P椭球，采用EGM2008大地水准面)及各项地球物理改正等信息，高程数据采取Hall和O'Loughlin使用的方法，利用高程使用标志(elev_use_flg)和饱和度指数标志(sat_corr_flg)去除或校正观测值;Cryosat-2卫星采用L2 GDR数据，该产品包括采样点经度、纬度、时间、大地高等信息，高程数据已经过仪器校正、传输延迟改正、几何改正和地球物理改正(主要包括海洋负荷潮汐、固体潮、极潮改正，干、湿对流层改正及电离层改正);ICESat-2 ATL13产品提供内陆水体的高度信息。研究区轨迹分布见图1，其中图1(a)和(b)分别为ICESat和Cryosat-2的轨道分布，因区域内ICESat-2数据较多，展示效果不好，图1(c)仅显示纳木错地区ICESat-2卫星轨道。
图1 研究区域内卫星轨道分布 Fig.1 Distribution of satellite orbits in the studied area
2 数据处理方法
2.1 参考椭球的统一
联合多源测高卫星可充分发挥各测高卫星的优势，增加观测湖泊的数量与时间序列的长度，但ICESat卫星采用T/P参考椭球，Cryosat-2、ICESat-2采用的是WGS84椭球，在数据处理之前首先需要统一基准。本文将所有卫星参考椭球统一至WGS84椭球，转换公式为：
[
egin{aligned}
Delta B&=frac{N}{(M+H)^{2}} e^{2} sin B cos Bleft(alpha_{TP}-alpha_{84} ight)+ \
&quad frac{Mleft(2-e^{2} sin ^{2} B ight)}{(M+H)(1-f)} sin B cos Bleft(f_{TP}-f_{84} ight) quad(1) \
Delta H&=-frac{N}{a}left(1-e^{2} sin ^{2} B ight)left(alpha_{TP}-alpha_{84} ight)+ \
&quad frac{M}{1-f}left(1-e^{2} sin ^{2} B ight) sin ^{2} Bleft(f_{TP}-f_{84} ight) quad(2)
end{aligned}
]
[
Delta L=0 quad(3)
]
[
B=B_{0}+Delta B quad(4)
]
[
H=H_{0}+Delta H quad(5)
]
[
L=L_{0} quad(6)
]
式中，(N) 和 (M) 分别为卯酉圈和子午圈曲率半径;(e) 为椭球第一偏心率;(a) 为椭球长半轴;(f) 为相应椭球的扁率;(B_{0})、(L_{0})、(H_{0}) 和 (B)、(L)、(H) 分别为转换前后测高点的大地坐标。
2.2 Cryosat-2重跟踪算法的选择
2010-2016年Cryosat-2卫星LRM测量模式的L2 GDR产品提供了Ocean-CFI、UCL Land-ice、OCOG三种重跟踪算法处理得到的湖泊水位高度数据，为选择适用于内陆湖泊水位变化研究的重跟踪算法，本文以青海湖为例，以单日轨道观测值的剔除率作为不同重跟踪算法数据质量检验的标准进行对比分析。
以天为单位，剔除单日所有观测值中超过3倍MAD(median absolute deviation)的数据。单日平均水位是利用卫星沿轨观测的高程统计出准确的湖面高程，由于沿轨的雷达回波质量不一，如靠近湖岸的回波可能被严重污染，由其得到的水位不准确，而湖中心的水位往往比较精确。若利用平均值或3倍标准差准则剔除异常观测值，容易受到离群值的干扰导致剔除结果不准确。因此，利用单日轨道观测值的MAD作为剔除异常值准则更为合理(式(7))，以单日所有观测值的剔除率作为数据质量检验标准(式(8))：
式中，(X_{i}) 为单日中第(i)个测量值;(X) 为单日所有观测值的一元序列;(R_{r}) 为单日数据剔除率;(N_{r}) 为单日数据剔除个数;(N) 为单日数据总数。
2010-2016年Cryosat-2卫星L2 GDR产品中3种重跟踪算法对应的湖泊水位高度观测值的数据质量见表1。OCOG算法的基本思想是通过经验算法找到每个雷达回波波形的重心位置，计算由波形值确定的矩形重心和面积，进而精确标识出波形前缘的中点位置[14]。根据表1所示，OCOG表现最佳，3种重跟踪算法剔除前的标准差分别为7.58m、24.60m、6.64m，(R_{r})的最大值、最小值、平均值分别为0.58、0、0.11。在数据离散度上，OCOG重跟踪算法同样表现最佳。因此，LRM测量模式下，本文选用OCOG重跟踪算法得到的湖泊水位数据。
2010-2016年Cryosat-2在青海湖区域的主要测量模式为LRM，而2016-2022年主要测量模式为SARIn。SARIn模式的数据产品采用UCL margin重跟踪算法，对该数据也采用类似的剔除粗差策略，最大单日数据剔除率为0.32，平均单日数据剔除率为0.07，均优于LRM测量模式下3种重跟踪算法数据产品，证明了SARIn测量模式的优势。
表1 LRM与SARIn测量模式的(R_{r})值统计 Tab. 1 (R_{r}) statistics of LRM and SARIn measurement modes
| 测量模式 | 重跟踪算法 | 最大值 | 最小值 | 平均值 |
| LRM | Ocean-CFI | 0.72 | 0 | 0.17 |
| | UCL Land-ice | 0.82 | 0 | 0.26 |
| | OCOG | 0.58 | 0 | 0.11 |
| SARIn | UCL margin | 0.32 | 0 | 0.07 |
2.3 湖泊水位时间序列的批量提取
利用3颗测高卫星数据提取2003-2022年青藏高原湖泊的水位序列，并分析水位变化趋势。湖泊边界沿用青藏高原大于1km²湖泊水量变化(1976-2020年)V2.0数据集，在理想情况下，将水面看作一个水位变化连续的静态平面，但偏离水位中心或陆地污染均会造成异常值，因此，生成日平均水位时间序列之前需剔除预处理后测高数据的异常值。
本文在剔除异常值时分为两步，分别为单日异常值剔除与时间序列(经过单日异常值剔除的联合多源测高卫星得到的水位时间序列)平滑。对单日观测水位剔除异常值，并求平均得到日水位，同时计算单日观测值的总数量与剔除率。因为根据较少的沿轨数据统计出来的湖泊水位不确定性较大，且本文均使用测量频率为20Hz的数据，总数小于20的单日测量值具有很大的不确定性，不能代表湖泊水位，因此不作考虑。
以青海湖为例，展示Cryosat-2卫星时间序列的生成过程。图2(a)为卫星观测点序列，存在较多异常值且分布离散;图2(b)为剔除异常值后的日平均水位;图2(c)统计了异常值剔除率;图2(d)为每日平均水位剔除不合格数据后的时间序列。2010-2016年青海湖观测数据共剔除29d不合格数据，而2016-2022年观测数据仅剔除2d不合格数据。
图2 青海湖Cryosat-2卫星单日异常值剔除过程 Fig.2 Single-day outlier removal process of Cryosat-2 satellite of Qinghai lake
2.4 时间序列小波去噪
由于地表情况及地球物理等因素的影响，直接利用测高卫星数据获得的水位时间序列含有噪声，需要对其进行滤波处理。同样以青海湖为例，图3显示了2003-2022年青海湖水位时间序列，ICESat、Cryosat-2与ICESat-2数据覆盖的时间段分别为2003-2009年、2010-2022年和2018-2022年，其中2018-2022年Cryosat-2与ICESat-2两类卫星数据同时存在。
图3 多源卫星数据散点图 Fig.3 Scatter plot of multi-source satellite data
本文采用基于sym小波基的小波方法去除噪声[15]，去噪后效果见图4。可以看出，2003-2022年青海湖水位持续上涨，2013-2016年水位增速减缓，趋于稳定;但2016年6月开始，水位上涨速度较2003-2013年快。根据采用的卫星数据的种类将时间序列分为2个时段，第1个时段为2003-2010年，数据来源于ICESat卫星，数据量较少;第2个时段为2010-2022年，数据来源于Cryosat-2与ICESat-2。2个时段观测数量不一致，无法对青海湖水位的周期变化进行分析，因此仅考虑长期趋势。利用式(9)对测高反演的湖泊水位变化时间序列进行拟合分析：
式中，(y) 为湖泊水位;(t) 为时间;(A) 为常数项;(B) 为长期趋势项。图5为青海湖时间序列分析，可以看出，2003-2022年青海湖水位变化速率为23.10±0.1cm/a，其中2003-2010年为10.84±0.1cm/a，2010-2022年为26.63±0.82cm/a。
图4 小波去噪前后的青海湖水位变化时间序列 Fig.4 Time series of water level changes in Qinghai lake before and after wavelet denoising
图5 青海湖水位时间序列分析 Fig.5 Time series analysis of water level in Qinghai lake
3 数据分析
3.1 青藏高原湖泊水位变化趋势
由于不同测高卫星的重访周期、轨道覆盖密度及运行时间不同，对2003-2022年青藏高原湖泊的年均观测频率进行统计，将年观测频率大于1的湖泊视为可监测动态变化的湖泊。结果发现，2003-2022年具有长期观测的湖泊数量为106个，大部分湖泊位于内流区的羌塘盆地，其中5个湖泊面积小于10km²，26个湖泊面积为10~50km²，21个湖泊面积为50~100km²，面积大于100km²的有54个。湖泊水位变化速率见图6。
图6 青藏高原湖泊水位变化速率 Fig.6 Change rate of water level in lakes in Qinghai-Xizang plateau
就湖泊水位变化趋势而言，19个湖泊水位呈下降趋势，下降速率最快的是位于雅鲁藏布江流域的羊卓雍错，水位变化速率为-0.23m/a，水位最大值与最小值之差约为5.4m。4个湖泊水位保持平稳趋势，83个湖泊水位呈上升趋势，上升速率最快的湖泊位于羌塘盆地库塞湖东侧约30km处，湖泊面积约为150km²，其上升速率达到0.8m/a，水位最大值与最小值之差约为12.75m。
从流域尺度分析，雅鲁藏布江流域长期观测的3个湖泊中，仅普莫雍错的水位呈缓慢上升趋势，上升速率为0.02m/a，其余2个湖泊均呈显著下降趋势。恒河流域仅佩估错存在长期监测，呈缓慢下降趋势，水位变化速率为-0.06m/a。印度河流域监测到的3个湖泊中，拉昂错水位呈显著下降趋势，速率达到-0.16m/a，其他2个湖泊水位呈缓慢上升趋势。羌塘盆地长期观测的90个湖泊中，14个湖泊水位呈下降趋势，其中下降速率最快的为纳木错北部面积约为50km²的湖泊，水位变化速率达到-0.16m/a;73个湖泊的水位呈上升趋势，上升速率最快的湖泊位于35.52°N、93.43°E，水位变化速率达到0.8m/a。柴达木盆地长期观测的4个湖泊中，面积小于50km²的2个湖泊水位保持稳定，面积大于100km²的2个湖泊水位均呈现显著上升趋势;位于37.62°N、93.68°E的湖泊面积从2005年的8km²增至2020年的384km²，是近20a来的典型新增湖泊，但湖泊水位仅升高约2.5m，说明湖泊较浅且向四周扩张。黄河流域与长江流域具有长期观测的湖泊水位均呈上升趋势，其中青海湖水位上升速率最快，约为0.23m/a。河西流域观测到1个湖泊呈上升趋势，位于38.29°N、97.59°E，变化速率为0.27m/a。
根据统计，水位变化剧烈的湖泊多位于内流区羌塘盆地，原因是羌塘盆地的湖泊数量巨大，面积超过1km²的湖泊有近500个，湖泊数量和面积分别达到内流区湖泊的86%和70%。相对于外流湖水位变化(变化速率约0.12m/a)更加剧烈，其中内流区的冰川补给湖数量为60个，占湖泊总数和总面积的80%和67%，平均变化速率约为0.24m/a;非冰川补给湖的平均变化速率约为0.20m/a。雅鲁藏布江流域、恒河流域湖泊整体呈下降趋势，其他流域湖泊水位均呈上升趋势。
3.2 利用地面测量数据验证卫星测高探测水位精度
为验证青藏高原湖泊水位反演精度，收集青海湖、纳木错、美玛错、鲁玛江东错的水位站实测数据，选取皮尔逊相关系数(( ho_{xy}))、均方根(RMS)来评估水位探测精度。利用青海湖下社站2003-2022年平均水位时间序列，对联合ICESat、Cryosat-2、ICESat-2卫星测高数据得到的青海湖水位年平均值进行验证，如图7所示。
图7 青海湖水位验证 Fig.7 Verification of Qinghai lake water level
卫星测高提取的水位精度较高，两者差异的最大值、最小值、平均值、均方根误差分别为0.109m、-0.170m、-0.037m、0.092m;最大差异出现在2009年，分析原因发现，2009年青海湖范围内ICESat卫星共有5个观测点，其中4个观测点采样时间为1-3月，以此计算年度水位平均值并不准确，因此与下社站年平均水位差异较大。另外，与下社站水位的比较也可反映不同测高卫星、不同测量模式测量水位的精度差异。由图7可见，2003-2010年，ICESat卫星在该区域的总观测值数量仅为43个，平均每年观测值数量约为6个，水位差异的均方根误差较大，约为0.127m;2010-2016年，Cryosat-2 LRM模式在该区域的总观测值数量为108个，年均观测值数量为18个，均方根误差约为0.078m;2016-2022年为Cryosat-2 SARIn模式与ICESat-2共同观测，总观测值数量为274个，年均观测值数量为45个，均方根误差为0.054m，可见Cryosat-2 SARIn模式和ICESat-2卫星探测内陆湖泊方面的精度优于Cryosat-2 LRM模式，ICESat卫星的精度最低。
为进一步对比Cryosat-2 SARIn模式与ICESat-2卫星的精度，分别利用下社站数据对2018-2022年两类观测数据精度进行检核，结果表明，ICESat-2数据的精度较低，均方根误差为4.15cm;Cryosat-2 SARIn数据的精度较高，均方根误差为3.47cm。
利用纳木错(2019-04-12)、美玛错(2017-2020年)、鲁玛江东错(2016-2022年)的水位数据，对3个湖泊的卫星测高提取的水位精度进行评价(图8)，并对2类数据的差异进行统计(表2)。由于无法得知验证数据集的高程基准，根据对应时间序列距平值，验证本文反演的水位精度。
表2 反演水位与验证数据差异统计 Tab. 2 Statistics of differences between inversion water level and validation data
| 湖泊 | 面积/km² | 最大值/m | 最小值/m | 平均值/m | 均方根误差/m | 相关系数 |
| 纳木错 | 2020.98 | 0.104 | -0.049 | 0.016 | 0.047 | 0.989 |
| 鲁玛江东错 | 409.04 | 0.178 | -0.132 | 0.063 | 0.094 | 0.986 |
| 美玛错 | 184.99 | 0.170 | -0.268 | -0.013 | 0.112 | 0.934 |
图8 利用地面测量验证反演水位精度 Fig.8 Verifying the accuracy of inverted water level using ground measurements
本文利用卫星测高提取的湖泊水位精度较高，与水位站数据的相关系数均达到0.93以上。由于面积较大的湖泊受陆地信号污染较少，使得面积大的湖泊较面积小的湖泊具有更高的水位反演精度。而纳木错湖泊面积最大，为2020.98km²，其对应的相关系数最高，为0.989，且均方根误差最低，为0.047m，说明提取的纳木错水位精度更高。
3.3 与Hydroweb水位数据集差异对比
因实测水位站数量和数据量有限，本文同时采用Hydroweb水位数据集对卫星测高提取的湖泊水位进行检核，Hydroweb水位数据集未采用ICESat、Cryosat-2、ICESat-2卫星数据，因此可作为独立数据集与本文反演结果进行对比分析。选取Hydroweb水位数据集中具有长时间序列的赤布张错、多格仁强错、鄂陵湖、拉昂错4个湖泊的水位值，与卫星测高提取的水位进行对比，差异统计信息见表3。
表3 与Hydroweb数据集的差异统计 Tab.3 Difference statistics with Hydroweb dataset
| 湖泊 | 2020年面积/km² | 均方根误差/m | 相关系数 |
| 赤布张错 | 342.56 | 0.25 | 0.99 |
| 多格仁强错 | 453.97 | 0.18 | 0.98 |
| 鄂陵湖 | 673.10 | 0.25 | 0.75 |
| 拉昂错 | 252.64 | 0.30 | 0.78 |
以上4个湖泊反演水位与Hydroweb水位数据集具有较高的吻合度，皮尔逊相关系数均大于0.75，其中赤布张错、多格仁强错与Hydroweb的相关性分别达0.99、0.98。观察差异值分布(图9)可知，差异较大值主要集中在2010年之前，ICESat观测值与Hydroweb存在明显的系统差，系统差可能来自于ICESat观测值与其他测高卫星重跟踪算法的不同。已有研究表明，对于Envisat和Jason-2来说，ICE-1和OCEAN重跟踪算法不同导致的高程偏差分别为0.24m和0.23m[6]。2016之后，Cryosat-2观测模型转换为SARIn模式，且融入了ICESat-2观测值，这也说明Cryosat-2 SARIn测量模式与ICESat-2数据的使用提高了反演水位时间序列的精度，因此卫星提取的湖泊水位与Hydroweb数据集吻合度更高。
图9 反演结果与Hydroweb水位数据集差异统计 Fig. 9 Statistical difference between inversion results and Hydroweb water level dataset
4 结语
本文利用卫星测高数据提取青藏高原湖泊水位变化时间序列的相关技术，选用ICESat、Cryosat-2与ICESat-2卫星测高数据获取青藏高原106个湖泊的水位变化信息。主要结论如下：
1) 对Cryosat LRM测量模式下Ocean-CFI、UCL Land-ice和OCOG三种重跟踪算法的水位变化信息提取效果进行比对分析，就精度和数据离散度而言，OCOG算法最优，但低于SARIn测量模式的精度。
2) 采用实测水文站数据对利用卫星测高数据提取的湖泊水位变化时间序列进行外部检核，结果表明，本文提取的湖泊水位变化信息精度较高，且相对于ICESat-2，Cryosat-2 SARIn精度略高。与Hydroweb独立数据集对比，结果表明，两者具有较高的拟合度，但Hydroweb独立数据集陈列的湖泊数量少，且均为中大型湖泊，缺少对小型湖泊的监测，而本文提取的湖泊水位序列可以弥补这一缺点。
3) 提取青藏高原106个湖泊水位序列并分析其变化趋势，结果表明，19个湖泊水位呈下降趋势，83个湖泊水位呈上升趋势。
参考文献
[1] Ma R H, Yang G S, Duan H T, et al. China's Lakes at Present; Number, Area and Spatial Distribution [J]. Science China Earth Sciences, 2011, 54(2): 283-289.
[2] 褚永海, 李建成, 姜卫平, 等. 利用Jason-1数据监测呼伦湖水位变化[J]. 大地测量与地球动力学, 2005, 25(4): 11-16.
[3] Cheng K C, Kuo C Y, Tseng H Z, et al. Lake Surface Height Calibration of Jason-1 and Jason-2 over the Great Lakes [J]. Marine Geodesy, 2010, 33(S1): 186-203.
[4] Tong X H, Pan H Y, Xie H, et al. Estimating Water Volume Variations in Lake Victoria over the Past 22 Years Using Multi-Mission Altimetry and Remotely Sensed Images [J]. Remote Sensing of Environment, 2016, 187: 400-413.
[5] Yuan C, Gong P, Liu C X, et al. Water-Volume Variations of Lake Hulun Estimated from Serial Jason Altimeters and Landsat TM/ETM+ Images from 2002 to 2017 [J]. International Journal of Remote Sensing, 2019, 40(2): 670-692.
[6] 孙明智, 刘新, 汪海洪, 等. 多源卫星测高数据监测拉昂错1992年-2020年水位变化[J]. 遥感学报, 2022, 26(1): 126-137.
[7] Zhang G Q, Xie H J, Kang S C, et al. Monitoring Lake Level Changes on the Qinghai-Xizang Plateau Using ICESat Altimetry Data(2003-2009)[J]. Remote Sensing of Environment, 2011, 115(7): 1733-1742.
[8] Jiang L G, Nielsen K, Andersen O B, et al. Monitoring Recent Lake Level Variations on the Qinghai-Xizang Plateau Using CryoSat-2 SARIn Mode Data[J]. Journal of Hydrology, 2017, 544: 109-124.
[9] 常翔宇, 蔡宇, 夏文韬, 等. Cryosat-2测高数据的太湖水位监测与变化分析[J]. 测绘科学, 2020, 45(10): 85-91.
[10] 廖静娟, 薛辉, 陈嘉明. 卫星测高数据监测青藏高原湖泊2010年-2018年水位变化[J]. 遥感学报, 2020, 24(12): 1534-1547.
[11] Schutz B E, Zwally H J, Shuman C A, et al. Overview of the ICESat Mission[J]. Geophysical Research Letters, 2005, 32(21).
[12] Cooley S W, Ryan J C, Smith L C. Human Alteration of Global Surface Water Storage Variability [J]. Nature, 2021, 591(7848): 78-81.
[13] Meloni M, Bouffard J, Parrinello T, et al. CryoSat Ice Baseline-D Validation and Evolutions[J]. The Cryosphere, 2020, 14(6): 1889-1907.
[14] 褚永海, 李建成, 张燕, 等. ENVISAT测高数据波形重跟踪分析研究[J]. 大地测量与地球动力学, 2005, 25(1): 76-80.
[15] Torrence C, Compo G P. A Practical Guide to Wavelet Analysis[J]. Bulletin of the American Meteorological Society, 1998, 79(1): 61-78.
[16] 袁翠. 基于多源卫星测高数据的全球湖泊动态监测(1992-2020)[D]. 北京: 清华大学, 2021.
转载请注明来自：http://www.lunwencheng.com/lunwen/lig/22722.html
