青藏高原多尔索洞错水深反演与水量变化估算

崔剑 乔宝晋 郭恒亮 杨洪

摘要：利用由实测水深空间插值获得的水下地形数据估算湖泊水量及其变化往往存在较大误差。以Landsat系列影像数据及多尔索洞错12 m以下实测水深数据为基础，建立多波段组合模型进行水深模拟，结合多时相Landsat影像获取的多尔索洞错边界数据获取水深变化，估算1996～2016年多尔索洞错水量变化。实验结果表明：多因子反演模型的相关系数（R2）均在0.90以上，平均绝对误差低至0.48 m，相比空间插值方法能更精确地模拟多尔索洞错水深在12 m以下区域的水深分布情况；近20 a来多尔索洞错不断扩张，水位升高约 0.40 m/a，水量增加约0.18 km3/a。

关键词：水深反演； 水量变化； 归一化水体指数法（NDWI）； 多波段组合模型； 多尔索洞错

中图法分类号： P332

文献标志码： A

DOI：10.16232/j.cnki.1001-4179.2024.04.019

0引 言

湖泊作为岩石圈、大气圈、生物圈、陆地水圈相互作用的重要中间组成部分［1］，其参与自然界的水分循环，是揭示全球气候变化与区域响应的重要信息载体，对区域及全球气候变化显得极为敏感［2］。因此，对湖泊面积、水位、水量等进行研究十分重要。近年来，随着科学技术的发展与经济条件的提高，中国在青藏高原湖泊研究中已有大量利用卫星测高数据及数字高程模型直接获取水位变化的研究［3］。李龙等［4］以2003～2009年ICESat/GLAS衛星测高数据为基础研究可可西里地区主要湖泊水位变化，并分析了该地区湖泊变化对气候的响应。戴玉凤等［5］利用2003～2011年Landsat ETM+数据和2003～2009年ICESat激光测高数据，获取青藏高原佩枯错的面积和水位变化，估算其水量变化。Wang等［6］利用ICESat/GLAS数据研究发现中国各地区湖泊水位变化情况各不相同。Long等［7］通过整合陆地卫星专题制图影像、冰、云及地面高程卫星来提高水位测量的精确度（ICESat/GLAS），使得水位误差在1 m以内。Yang等［8］利用SRTM数据估算了青藏高原的湖泊水量变化情况。但结合遥感反演技术反演湖泊水深并估算水量变化的研究还很少。

遥感水深反演是一种间接估算水深的方法。反演模型主要有统计相关模型、理论模型、半理论半经验模型［9］。然而，水体内部的许多光学参数受条件限制无法确定，有条件可获取的参数也与卫星过境时间不能同步，使得相关统计模型在遥感水深反演中应用广泛［10］。田庆久等［11］利用多光谱遥感信息传输方程推导出水深信息对数反演模式，根据TM影像的第3和第4波段，建立了江苏近海辐射沙脊群海域水深反演模型，平均绝对误差达到1.857 m，但15 m以上反演误差较大；张靓等［12］提出了适用于海岸带地区的非线性优化水深遥感探测方法，模型决定系数R2高于0.81；李畅游等［13］以LandsatTM/ETM+影像数据及呼伦湖分布的实测水深点为基础，应用太阳光波段和热红外辐射波段组合模型进行多个时间段的水深模拟，模型绝对误差大多在0.8 m以内，相对误差大多在30%以内；Guo等［14］将统计模型用于反演青藏高原湖泊邦达错的浅水区水深，平均相对误差低于20%；Shen等［15］对Lyzenga方法进行了改进，引入广义相加线性模型来描述水深与图像回归方程之间的非线性关系，精度提高了20%左右；林征等［16］利用Stumpf水深反演模型对北极3 187个湖泊水深进行反演，避免了底质对反演的影响，说明了遥感反演在极地湖泊反演中的适用性，验证其中一个湖泊的平均绝对误差为 0.37 m；Yunus等［17］利用Sentinel-2和Landsat-8影像建立经验模型与随机森林模型，后者平均绝对误差均低于1 m，说明了两种影像在海岸和湖泊环境反演制图中的适用性。相关研究主要分析模型精度、水深分布情况，很少进一步研究水量及其变化。

综上，湖泊水位与水量变化研究主要利用测高数据及数字高程模型，而遥感水深反演局限于研究一时一湖的反演模型精度、水深分布情况、水下地形等，对反演的最大水深也有一定的限制。同时，遥感水深反演在湿润地区的湖泊及浅海海域研究较多［18-19］，而对于寒旱区湖泊水深反演的研究很少，青藏高原的大部分湖泊地理位置偏远，水深测量需要花费大量的人力物力，目前只有少数已发表的文章中涉及到个别湖泊的水深研究［20-21］。因此，利用遥感水深反演结合湖泊边界变化估算水位及水量变化在可行性与精度方面都有待进一步研究。本文通过选取2016年9月多尔索洞错实测水深在12 m以下的数据与Landsat-8影像各波段反射率及其组合建立水深反演模型模拟湖泊水深，得到12 m以下区域的水深分布情况，利用1996～2016年的湖泊边界提取反演深度并求其均值，进而估算近20 a的水深变化及水量变化情况。

1研究区与数据

1.1研究区概况

多尔索洞错位于中国西藏自治区那曲地区，唐古拉山腹地，双湖特别区东部，海拔4 920 m左右，地理坐标为89°37′E～90°0′E，33°13′N～33°32′N。根据相关气象站的数据，受南部季风和北部干旱大陆风的影响，该地区的年平均气温约为-4°C，年降水量小于500 mm，其中5～9月降水量占全年降水量的90%左右［22］。湖面分布形状呈不规则斜长方形，轴为东北至西南方向，湖面水域多年平均面积在400 km2左右，平均水深在30 m左右，是长江源区大型湖泊之一。与赤布张错通过一条河道相连接，且水位略低于赤布张错。因此，多尔索洞错水位变化除受降水、冰川融水等影响外，赤布张错通过河道补给多尔索洞错，两个湖泊形成了一个内陆单一湖泊体系，在过去的10 a里，多尔索洞错的水位迅速上升［23］。

1.2遥感影像数据来源与处理

多尔索洞错地处青藏高原，夏季多风，湖面受风浪起伏影响，约在每年的9月底天气开始转冷，水温迅速下降，这一时期风浪较小，湖面较为平静，入湖水量也因进入枯水期而大量减少，湖水深度波动小，水气温度相差大，湖水温度下降快，而11月中旬开始进入结冰期，直到来年5月初［24］。因此用于反演的最佳影像时间为9～11月。此外，影像的选取优先考虑多尔索洞错上无云覆盖时拍摄的。研究中使用的Landsat系列遥感影像从美国地质勘探局（https：∥earthexplorer.usgs.gov/）获取，详细信息见表1。利用ENVI 5.5对影像进行辐射定标、大气校正预处理后获得各波段反射率并进行水陆分离。

1.3实测水深数据来源与处理

本文使用的实测水深数据来自Qiao等［25］已公开发表的数据（https：∥www.scien- cedirect.com/Science/article/pii/S0022169419307796#s0105），于2016年用Lowrance HDS5测深仪测定，该仪器每秒记录一次水深，垂直精度达0.01 m，实测点分布如图1所示。受船只吃水深度及湖面环境条件影响，本次测量的最浅水深为0.93 m，每一个数据都记录了测深点的深度及相应的经纬度坐标。因使用多波段组合模型来反演水深，实测水深样本点的多少在很大程度上决定了模型精度，以及光在水中的衰减性，利用光学影像进行水深反演的方法一般只在浅水区有效。Qiao等［25］利用ICESat和ICESat-2对水位的研究结果表明，在2011～2014年间多尔索洞错出现最大平均水位变化速率为0.53 m/a，为保证本次反演的深度包含多尔索洞错近20 a的总体水位变化，首先选取0～12 m部分的实测水深数据，利用ArcGIS 10.2软件“创建渔网”工具生成与Landsat 8影像空间分辨率相同的渔网，然后求取每一个30 m×30 m网格里的均值作为对应影像像元水深的表征值（见图2）。

2研究方法

2.1反演因子選取与模型建立

由于水体类型与底质空间变异性较大，单一波段建立的反演模型容易缺失水体的部分光学特性，导致模型的精度较低［26］。因此本文拟建立多波段组合模型，首先对水深表征值和Landsat-8影像各波段的反射率及任意两波段比值进行相关性分析，以选取建模因子，水深与各因子的相关系数（R2）见表2。其中OLI1、OLI2、OLI3、OLI4波段分别为蓝绿、绿、红、近红外波段；OLI6为热红外波段；OLI5、OLI7为短红外波段。

由表2可以看出，2016年实测水深数据与各波段相关性较好的有波段OLI1、OLI2、OLI3、OLI4，其中OLI3与水深的相关系数最高，两波段比值与水深的相关性有较大提高，由于两波段之间比值组合较多，表中只列出OLI1、OLI2、OLI3、OLI4四个波段的相互组合中与水深相关性较高的组合。拟选择波段OLI1、OLI2、OLI3、OLI4及OLI3/OLI1、OLI3/OLI2分别取对数和非对数进行组合后建立反演模型，并对建模与检验结果进行误差分析，选出8种模型中的最优模型。使得模拟的水深分布尽可能地接近真实情况，对比模型如下。

2.2模型的计算与检验

对8种模型进行回归分析，用相关系数（R2）与均方根误差（RMSE）、平均绝对误差（MAE）和平均相对误差（MRE）评估模型精度［27］，从而选出最优模型。其中：

2.3多尔索洞错面积提取与水量计算

（1） 面积提取。

水体指数法是通过分析水体的光谱特性， 得到水体最弱反射波段与最强反射波段，将最强波段与最弱波段通过比值运算得到一幅比值增强图像，增大相邻像元的差别，削弱外界环境的影响。在这个过程中其他地物均受到抑制，从而达到提取研究需要的水体的目的。本文采用归一化水体指数法（NDWI）提取研究区多尔索洞错，公式为［29］

（2） 水量计算。

本文利用经验公式计算多尔索洞错水储量，计算公式为［30］

3结果与分析

3.1水深反演结果及精度验证

3.1.1水深反演模型与精度验证

对8种模型进行精度分析，得到结果见表3。

由表3可知，8种模型相关系数R2均在0.80以上，采用对数关系拟合反射率与实测水深值相关性较高的有对数反射波段、对数反射波段+比值波段、对数反射波段+对数比值波段3种组合模型，R2分别达到了0.936 7、0.938 7和 0.937 6；均方根误差分别为0.62，0.61 m和0.62 m；平均绝对误差分别为0.49，0.48，0.49 m；平均相对误差分别为8.29%，8.22% 和8.29%。检验相关系数均略小于建模相关系数，说明拟合数据的整体相关性较好。计算这3种模型的实测值与模拟值的绝对误差与相对误差，以更好说明模型拟合效果的优劣，绝对误差主要在0.2～0.8 m之间，均为检验大于建模；相对误差值主要在5%～14%之间，也均为检验大于建模。

综合对比8种模型，对数反射波段+比值波段模型优于其他模型，该模型为本次研究中反演多尔索洞错的最优反演模型，能更精确地模拟多尔索洞错12 m以下水深分布情况，其表达式为

3.1.2最优模型反演结果

运用式（15）对2016年的多尔索洞错进行反演，用参与建模与检验的实测水深数据提取最优模型反演的水深并进行相关性分析，结果见图3。利用空间分析工具划分不同水深区域分布情况，再生成12 m等深线，结果见图4（a）。

由于本次研究为了提高水深反演的精度用于计算水量变化，只选取了12 m以内的实测水深数据作为建模数据与检验数据，超过12 m深度并未参与模型建立，因此模型反演的多尔索洞错最大深度只有18.56 m，与实测最大深度68.71 m有较大误差。但从1996～2016年多尔索洞错边界线及12 m等深线比较发现，12 m等深线在12 m实测点附近且1996年多尔索洞错边界在12 m等深线以外，即1996～2016年多尔索洞错水深变化没有超过12 m，实测水深与反演水深相关性R2=0.932 7，也进一步说明建立的多波段组合模型用于反演多尔索洞错0～12 m之间的水深取得了较好的结果。为了进一步说明在没有实测值区域遥感水深反演的准确性与优越性，对所有实测值经 ArcGIS 10.2交集制表后利用地形转栅格工具进行插值，插值结果如图4（b），分析见3.3节。

3.2多尔索洞错面积变化

研究区统计分析得到1996～2016年多尔索洞错面积变化情况如图5～6所示，结果表明多尔索洞错在20 a间，面积持续增加，年均变化速率为4.25 km2/a，具有显著扩张趋势。从图6曲线可以看出，近几年多尔索洞错面积扩张趋势有所减缓。湖泊快速扩张阶段处于2001～2007年间，年均变化速率为6.50 km2/a。从2008年开始湖泊变化处于缓慢增长状态，2008～2012年年平均变化速率仅为2.12 km2/a；2013～2016年湖泊面积基本保持稳定。此外多尔索洞错的西岸、北岸、南岸面积扩张较大，东南岸与西北岸也有明显扩张，西南岸扩张相对较小，经过近邻分析后统计得到多尔索洞错各方向每年平均扩张距离见表4。经计算，各方位平均每年水平扩张距离约为29.08 m，接近影像像元大小，说明同一方向相邻两年湖泊边界能够从不同像元提取水深值，用边界提取反演水深作为水深变化具有合理性。

3.3多尔索洞错水位与水量变化

利用多尔索洞错边界分别对最优模型反演结果和空间插值结果进行掩膜提取，得到每一年多尔索洞错边界处所有像元的均值即为相应的水位变化。通过分析，对于遥感反演结果，1996～2016年湖泊水位变化也呈现出显著的升高趋势。1996～2016年的水位总共上升了约8.05 m，湖泊水位在 1996～2002年稳步上升，6 a水位上升约为1.20 m，上升幅度为0.20 m/a，而在2003～2006年间却呈现出剧烈升高的趋势，年均升高幅度高达0.80 m/a，接下来的2008～2012年湖泊水位又恢复了平稳升高的趋势，年均升幅约为0.46 m/a（图7）。而对于空间插值，由于实测点分布不均匀，在远离实测点的区域插值结果明显偏小，如图8所示，分别对两种结果生成3，6，9，12 m等深线。从图8的A、B、C能清楚地发现在A（图8（a））、C（图8（c））处有实测点，两种方法生成的等深线则基本一致，趋近重合；B（图8（b））处没有实测点，相应的等深线偏离较远，即没有实测点的区域通过空间插值生成等深线均在通过反演生成相应的等深线内侧，且趋势随着远离实测点的距离逐渐增加。通过计算空间插值的水位变化发现1996～2016年水位共上升 3.89 m，年均增加仅有0.19 m，与实际水位变化相差较大，而Zhang等［31］利用测高卫星数据估算了2003～2018年多尔索洞错水位的年均增加速率为0.34 m/a。因此，本文通过反演得到2003年到2016年平均水位增加速率为0.44 m/a与真实结果相近，进一步证实了通过遥感水深反演方法计算水位及水量变化的可行性。

利用反演结果得到的水位变化计算多尔索洞错的水量，1996～2016年增加水量约3.66 km3，年均增加速率为0.18 km3/a。其中1996～2002年这6 a间增加速率相对较慢，共增加约0.48 km3，年均以0.08 km3/a的速率增加；2003～2007年增加速率加快，4 a共增加水量约1.35 km3，年均增加速率0.30 km3/a，2007年以后增加速率变缓（见图8）。年际变化与水位变化呈现出一致的变化趋势。

讨 论

4.1反演因子选取

本文使用Landsat系列遥感影像的分辨率为30 m，因此影像中一个像元在实际空间中代表的面积为900 m2。当底质起伏变化很大、性状不均匀时，像元点所反映出的水深信息并不能夠很好地代表该范围内的水底地形，对应点的波段反射率与相应水深应有的反射率不能很好地对应，从而降低了实测水深与反射率之间的相关性，增加了误差。此外，在建模因子选取中，7个波段的组合很多，为了减小工作量，本文只考虑了OLI1～OLI4波段，对这4个波段进行部分组合得到其中的最优模型，但实际上还有很多组合没有尝试，得到的模型不一定是所有模型中的最佳模型，模型的选择对水深反演精度的影响，有待进一步的研究。

4.2水深提取准确性分析

本文研究通过多尔索洞错边界对其2016年反演水深掩膜提取得到水深变化，尽管使用的是30 m空间分辨率的Landsat系列影像，但一个象元大小为30×30 m，得到的水深是边界所在象元的平均水深，并不能很好地代表边界处的深度，增大了提取水深的误差。但边界经过的象元数量较多，均在8 000个左右，最终水深变化由边界经过的所有象元求均值得到，避免了通过湖泊局部的水深变化代表整体的水深变化，从而很大程度上减小了误差，使得整体水深变化接近真实值。因此，本文通过水深反演进而提取水深变化的方法避免了对测高数据进行处理得到水位变化的繁杂过程，对获取水深变化进而研究水量变化具有一定的参考价值。

5结 论

本文利用多尔索洞错2016年实测水深数据以及与测深时间相近的Landsat-8遥感影像的波段反射率进行相关分析，选取反演因子，建立了最优多波段水深反演模型。采用归一化水体指数法（NDWI）提取1996～2016年多尔索洞错边界，进而提取对应边界的模型反演深度作为水深变化用于计算整体水量变化并与空间插值结果对比。主要结论有：

（1） 与单一波段相比，多波段组合反演模型结合了多个波段的水体反射率特性，能在很大程度上提高波段反射率与实测水深之间的相关性，从而提高模型反演精度。结果表明，实验建模与检验误差在0.50 m左右，能较为精确地模拟多尔索洞错水深12 m以下区域的水深分布情况，进而估算水量变化。

（2） 多尔索洞错的面积从1996～2016年持续增大，扩张趋势明显，20 a间面积增加近85 km2，年均增加约4.25 km2。在多尔索洞错的北岸、东南岸、南岸、西岸均有明显的扩张，西南岸的扩张相对较小，平均每年水平扩张距离29.08 m，与影像像元分辨率相近，用边界提取反演水深作为水深变化具有合理性。

（3） 结合反演的水深以及湖泊面积数据，1996～2016年间多尔索洞错水位上升约8.05 m，年均升高 0.40 m，该结果与Zhang等［31］利用ICESat和ICESat-2得到的结果一致；水量增加约3.66 km3，年均增加水量0.18 km3。该研究结果表明，结合水深数据及遥感反演不仅能够有效估算湖泊水量，还能够进一步提升湖泊水位及水量变化的估算精度，为快速扩张类型湖泊水量变化估算提供新思路。

参考文献：

［1］鲁安新，姚檀栋，王丽红，等.青藏高原典型冰川和湖泊变化遥感研究［J］.冰川冻土，2005，27（6）：783-785.

［2］CARROLLL M，TOWNSHEND J，DIMICELI C M，et al.A new global raster water mask at 250m resolution［J］.International Journal of Digital Earth，2009，2（4）：291-308.

［3］王志杰，李畅游，张生，等.基于水平衡模型的呼伦湖湖泊水量变化［J］.湖泊科学，2012，24（5）：667-674.

［4］李龙，姚晓军，李风贤，等.基于ICESat/GLAS数据的可可西里地区湖泊水位变化研究［J］.中国农业资源与规划，2019，40（3）：45-52，60.

［5］戴玉凤，高杨，张国庆，等.2003～2011年青藏高原佩枯错相对水量变化及其对气候变化的响应［J］.冰川冻土，2013，35 （3）：723-732.

［6］WANG X W，GONG P，ZHOU Y Y.Water-level changes in China′s lakes determined from ICESat/GLAS data［J］.Remote Sensing of Environment，2013，132：131-144.

［7］LONG Y N，YAN S X，JIANG C B，et al.Inversion of lake bathymetry through integrating multi-temporal Landsat and ICESat imagery［J］.Sensors，2019，19（13）：2896.

［8］YANG R M，ZHU L P，WANG J B，et al.Spatiotemporal variations in volume of closed lakes on the Tibetan Plateau and their climatic responses from 1976 to 2013［J］.Climatic Change，2018，140：621-633.

［9］王艷姣，董文杰，张培群，等.水深可见光遥感方法研究进展［J］.海洋通报，2007，26（5）：92-101.

［10］ LIU Z.Bathymetry and bottom albedo retrieval using hyperion：a case study of Thitu Island and Reef［J］.Chinese Journal of Oceanology and Limnology，2013，31（6）：1350-1355.

［11］田庆久，王晶晶，杜心栋.江苏近海岸水深遥感研究［J］.遥感学报，2007（3）：373-379.

［12］张靓，滕惠忠，孟婵媛，等.基于半分析模型的高光谱遥感水深探测方法［J］.海洋测绘，2011，31（4）：17-21.

［13］LI C Y，SUN B，JIA K L，et al.Multi-band remote sensing based retrieval model and 3D analysis of water depth in Hulun Lake，China［J］.Mathematical and Computer Modelling，2013，58（3/4）：765-775.

［14］GUO H L，YANG H，QIAO B J，et al.Multi-resolution satellite images bathymetry inversion of Bangda Co in the western Tibetan Plateau［J］.International Journal of Remote Sensing，2021（21）：8077-8098.

［15］SHEN X，JIANG L M，LI Q Q.Retrieval of near-shore bathymetry from multispectral satellite images using generalized additive models［J］.IEEE Geoscience and Remote Sensing Letters，2019，16（6）：922-926.

［16］林征，黎夏，乔纪纲.阿拉斯加北极滨海平原极地湖泊的水深遥感反演［J］.中山大学学报（自然科学版），2012，51（3）：128-134.

［17］YUNUS A P，DOU J，SONG X，et al.Improved bathymetric mapping of coastal and lake environments using sentinel-2 and landsat-8 Images［J］.Sensors，2019，19（12）：2788.

［18］TRAGANOS D，POURSANIDIS D，AGGARWAL B，et al.Estimating satellite-derived bathymetry （SDB） with the Google Earth Engine and Sentinel-2［J］.Remote Sensing，2018，10（6）：859-878.

［19］MA S，TAO Z，YANG X，et al.Bathymetry retrieval from hyper-spectral remote sensing data in optical-shallow water［J］.IEEE Transactions on Geoscience and Remote Sensing，2014，52（2）：1205-1212.

［20］ZHANG B，WU Y H，ZHU L P，et al.Estimation and trend detection of water storage at Nam Co lake，central Tibetan Plateau［J］.Journal of Hydrology，2011，405（1-2）：161-170.

［21］JIANG L G，NIELSEN K，AADERSEN B，et al.Montoring recent lake level variations on the Tibetan-Plateau using cryoSat-2 SARIn mode data［J］.Journal of Hydrology，2017，544：109-124.

［22］陈军，汪永丰，郑佳佳，等.中国阿牙克库木湖水量变化及其驱动机制［J］.自然资源学报，2019，34（6）：1345- 1356.

［23］QIAO B J，ZHU L P，WANG J B，et al.Estimation of lake water storage and changes based on bathymetric data and altimetry data and the association with climate change in the central Tibetan Plateau［J］.Journal of Hydrology，2017，578：56-67.

［24］KLEINHERENBRINK M，LINDENBERGH R C，DITMAR P G.Monitoring of lake level changes on the Tibetan Plateau and Tian Shan by retracking cryosat SARIn waveforms［J］.Journal of Hydrology，2015，521：119-131.

［25］QIAO B J，ZHU L P，WANG，J B，et al.Estimation of lakes water storage and their changes on the northwestern Tibetan Plateau based on bathymetric and Landsat data and driving force analyses［J］.Quaternary International，2017，454：56-67.

［26］曹斌，邱振戈，曹彬才.四種遥感浅海水深反演算法的比较［J］.测绘科学技术学报，2016，33（4）：388-393.

［27］MOKNATIAN M，PIASECKI M，MOSHARY F，et al.Development of digital bathymetry maps for lakes azuei and enriquillo using sonar and remote sensing techniques［J］.Transactions in GIS，2019，23（4）：841-859.

［28］薛峭，梁开，甘华阳，等.基于Landsat8 OLI数据的后海水深反演［J］.海南热带海洋学院学报，2017，24（5）：63-68.

［29］ASIAN A，RAHMAN A F，WARREN M W，et al.Mapping spatial distribution and biomass of coastal wetland vegetation in indonesian papua by combining active and passive remotely sensed data［J］.Remote Sensing of Environment，2016，183：65-81.

［30］PACHECO A，HORTA J，LOUREIRO C，et al.Retrieval of nearshore bathymetry from Landsat 8 Images：a tool for coastal monitoring in shallow waters ［J］.Remote Sensing of Environment，2015，159：102-116.

［31］ZHANG G Q，CHEN W F，XIE H J.Tibetan Plateau′s Lake level and volume changes from NASA′s ICESat/ICESat-2 and Landsat Missions［J］.Geophysical Research Letters，2019，46（22）：13107-13118.

（编辑：黄文晋）