联合Landsat影像和ICESat测高数据估计青海湖湖泊水量变化

吴红波，陈艺多

(1. 陕西理工大学 地理科学系，陕西 汉中 723000； 2.西北大学 陕西省地表系统与环境承载力重点实验室，陕西 西安 710127； 3.中国科学院青藏高原地球科学卓越创新中心，北京 100101)

1 研究背景

湖泊作为陆地水体的一部分，也是陆地水循环的重要储蓄场所。在气候变暖背景下，青藏高原的高山湖泊对气候变化的信息较为敏感，是指示气候变化的重要载体和指示器，能够敏感地指示某一区域的降水、蒸发、温度和湿度气象要素变化趋势[1]，因而深受国内外学者的广泛关注[2-4]。青海湖是我国湖泊面积最大的内陆湖泊，位于青藏高原东北边缘，处于青藏高寒区、黄土高原和西北干旱区的过渡地带，是研究气候响应、水量平衡及水文过程的理想场所[5]。内陆湖泊水量波动不仅受到气候变化影响，而且也受到人类活动干扰。为了揭示青海湖湖泊水量变化规律和特征[6]，急需借助星载遥感技术对青海湖面积、水位、水量的变化进行监测，揭示高山内陆湖泊水量变化规律及气候响应，为青海湖流域生态环境保护提供科学依据。

随着星载测高技术的发展与应用，星载测高雷达和激光雷达测高技术具有了全球覆盖能力[7]，能够较大范围、周期性探测和记录地表水体和海平面的水位变化，在地表高程和地物垂直结构观测方面具有较大潜力和优势[8]。关于利用星载ICESat、Jason、EnviSat、Cryosat等估计全球典型湖泊水位变化的研究[9-11]，国内外已有报道[12-13]。李生生等[14]利用Landsat-8 OLI数据的绿光和近红外波段提取青海湖湖体边界，在水体与湖滩地的边界准确性和连续性方面具有明显的改进。赵云等[15]利用主波峰重心偏移法、主波峰阈值法、主波峰5-β参数法、传统重心偏移法、传统阈值法和传统5-β参数法6种算法对Cryosat-2/SIRAL LRM 1级数据进行波形重跟踪，提取青海湖水域日均水位，并结合EnviSat/RA-2 GDR数据，获得2002-2015年的青海湖水位时变序列。此外，青海湖湖泊水域面积的遥感监测研究，主要以季节性和年际变化特征为主[16-18]。骆成凤等[19]利用长时间序列的MODIS影像，通过人工提取湖岸线信息，重建青海湖水域面积年变化序列。张洪源等[20]用MODIS影像和LEGOS(Laboratoire d’Etudes en Géophysique et Oc-éanographie Spatiale)的水位数据，基于湖泊水位-面积关系，揭示了青海湖2001-2016年期间湖泊水量年内、年际变化规律。卢善龙等[21]基于空间分辨率250 m的MOD09数据，采用阈值分割法提取湖泊水域信息，构建了2000-2012年青藏高原地区湖泊面积逐年变化数据集。虽然基于Landsat、MODIS多光谱遥感数据提取湖泊水域面积较为成熟[22]，但是联合星载测高数据和多光谱遥感影像估计青海湖湖泊水量变化的相关研究还较少，亟待通过较长时间序列的湖泊水位与水域面积数据估计青海湖湖泊水量变化，探讨湖泊水量变化估计的可行性。

青藏高原的高山湖泊地处偏僻，人类活动较少，湖泊的气象水文观测资料匮乏，且观测时段缺乏连续性，对湖泊水文特征、水量平衡和水储量变化的认知不足。为此，本文首先借助1988-2018年Landsat TM/ETM/OLI影像数据和归一化差异水体指数，提取青海湖湖泊的当日水域面积；其次，利用2003-2009年ICESat测高资料提取湖泊水域的日均水位；再次，利用湖泊水位-面积的相关关系式，重构1988-2018年青海湖面积-水位-水量变化序列，并结合水文观测和实地观测数据对湖泊水位估计精度进行验证；最后，分析了青海湖湖泊水位年际和年内变化特征和趋势，其结果可为青海湖流域水资源利用、生态环境保护和气候变迁研究提供理论参考和技术支持。

2 数据来源与研究方法

2.1 数据来源

2.1.1 GLAS数据 搭载在ICESat卫星上的GLAS(geoscience laser altimeter system)激光测高仪采用绿光543 nm和近红外1 064 nm处激光脉冲获取地表信息特征[23]，在平地上形成一个直径约70 m的圆形光斑，沿升/降轨道方向上的GLAS光斑脚点间隔约170 m。GLAS测高仪每秒发射40个激光脉冲，脉冲宽度为4 ns，1 ns的脉冲对应于0.15 m的高度[24]。本文采用Level 1B级GLA14产品提取GLAS光斑脚点经度、纬度、高程、过境日期和记录时间、大地水准高度等参数，用于陆地高程和湖泊水位信息估计。GLA01产品中提取光电信号延迟校正、轨道高度、波形记录参数等，用于进行湖泊水位估计的修正。GLA01产品版本为33，GLA14产品版本为34。青海湖流域及湖泊水域内有效的GLAS光斑分布见图1，共计10 678个；青海湖湖面的GLAS光斑汇总见表1，可从美国国家冰雪数据中心NSIDC(National Snow and Ice Data Center，URL：http://nsidc.org/data/icesat/)获取。

图1 青海湖流域及ICESat测高有效GLAS光斑分布图

2.1.2 Landsat TM/ETM/OLI数据 选取1988-2018年期间青海湖湖泊区无积雪、云量小于15%且湖泊水域无封冻的Landsat TM/ETM/OLI影像，用于提取湖泊水域面积。Level 1B级的Landsat影像可通过美国地质勘探局USGS(United States Geological Survey，URL：https://glovis.usgs.gov/)下载。行号Path=133、列号Row=34和行号Path=133、列号Row=35的两景Landsat TM/ETM/OLI影像可覆盖青海湖水域范围，共计175景。

2.1.3 观测资料 青海湖年初、年末、年均水位和水量变化数据来源于《青海省水资源公报》和青海水利信息网(URL：http://slt.qinghai.gov.cn/)，其中水域面积包括沙岛湖、海晏湾，不含尕海、错果。下社水文站的水位观测时段为当年5-10月份，当日水位观测值与ICESat的水位估计值进行对比验证。1995-2018年青海湖湖泊水位资料从法国地球物理学和海洋学太空观测研究中心LEGOS(URL：http://www.legos.obs-mip.fr/observations)获取。SRTM3数字高程模型产品，空间分辨率为90 m×90 m，从中国科学院地理空间数据云(URL：http://www.gscloud.cn/)下载。2018年9-10月期间，在青海湖湖泊水域南缘(S1区和S2区)和北缘(N1区和N2区)(见图1)，用分米级差分GPS采集仪对湖泊水域边界和平地高程进行测量，并随机选取平地上32个地面控制点用于GLAS光斑脚点高程的精度检验。考虑到GLAS光斑脚点在投影时，会有一定水平位移，利用二阶多项式模型将SRTM3像元和GLAS光斑脚点的地理位置精准配准，再做GLA14产品中高程和SRTM3的高程的误差分析。

表1 青海湖湖面的GLAS光斑汇总

2.2 技术路线

首先，借助1988-2018年Landsat TM/ETM/OLI影像和归一化差异水体指数(normalized difference water index，NDWI)提取湖泊水域面积；其次，联合2003-2009年期间ICESat-GLAS数据中的GLA01和GLA14产品提取湖泊水域瞬时水位高程和地表高程，并结合青海湖下社水文站水位数据和地表高程测量值，对GLAS光斑脚点高程做配准误差和标准误差分析；再次，通过湖泊水位-面积的关系建立线性回归关系式，用于重建非监测时段内的湖泊水位；最后，根据湖泊面积-水位-水量的绳套关系，构建1988-2018年青海湖湖泊水量时变序列。技术路线详见图2。

图2 技术路线

2.3 湖泊水量估计方法

2.3.1 湖泊水位估计 ICESat卫星测高数据估计湖泊瞬时水位Hlake的公式为：

Hlake=Hsat-Crange+Cdelay+Cpressure+Cwet+Cst+

Cpt+e

(1)

式中：Hsat为ICESat卫星飞行高度，km；Crange为ICESat卫星到地表距离，km；Cdelay为电离层传播延迟修正，m；Cpressure为大气气压变化所所引起的信号延迟修正，m；Cwet为大气湿度变化所引起的信号延迟修正，m；Cst为地壳运动所引起的垂直高度修正，m；Cpt为潮汐变化所引起的高程修正，m；e为未考虑的不确定性误差，m。

(2)

式中：halt,lake为Earth Gravitational Model 2008(EGM2008)重力位模型下的地表高程[8]，m；a为斜率；b为高度偏移常数[25]，m。

(3)

2.3.2 湖泊水域面积提取步骤

(1)对1988-2018年期间青海湖湖泊水域清晰且无云、少云、无雪覆盖的Landsat TM/ETM/OLI影像绿色和近红外波段做地理配准、辐射校正和大气校正；

(2)选取1988年1月22日的Landsat TM影像的湖泊水域边界，作为湖泊水域面积变化的参考边界；

(3)由于裸地与水体在近红外和绿色波段反射率的差异，归一化差异水体指数RNDWI可突出水体信息[28]，抑制地表土壤、植被信息。利用RNDWI阈值分割进行陆地和水体分类，其计算公式为：

(4)

式中：rgreen为Landsat TM/ETM/OLI影像的绿色波段反射率；rnir为近红外波段的反射率。RNDWI取值范围为-1～1；当像元的RNDWI≤0时，地物类型为陆地、植被、建筑等；当像元的0

(4)根据预设的RNDWI≥0.15，统计出第ti时期Landsat TM/ETM/OLI影像中青海湖水域水体像元的数量，并计算青海湖湖泊的水域面积Alake,ti：

(5)

式中：l为第ti时期水体像元数量，i=1，2，…，n；alake,ti为任一水体像元的面积，km2。

(5)当确定青海湖水域界限和面积后，湖泊水域面积在提取过程中的误差se为[29]：

(6)

式中：se为任一时期湖泊水域面积误差，km2；λ为Landsat TM/ETM/OLI影像中绿色和近红外波段空间分辨率，像元大小为30 m×30 m；εgeo为配准误差，m。GPS手持仪获取的湖泊边界线和地面控制点，通过均匀布设尽可能多的控制点，使水域边界的配准误差绝对值小于5 m。

2.3.3 湖泊水量变化估计 湖泊水量变化量ΔVlake(t)为湖泊水域面积和水位随着时间t变化的关系式[30]，其估计式为：

(7)

式中：ΔVlake(t)为湖泊水量随时间t变化的关系式，km3；S(t)为湖泊水域面积随时间t变化的关系式，km2；H(t)为在t时刻的湖泊水位，m；H(t-1)为在t-1时刻的湖泊水位，m。

Landsat4-8卫星重访间隔16 d，但由于云层、薄雾限制，有的月份无影像可用，本文利用湖泊面积和水位的关系式插补缺失的观测时段湖泊水域面积。如果观测时段内湖泊水位波动较小，那么湖泊水量变化可以用一系列连续的截锥体之和近似代替。则从ti+1时刻到ti时刻的湖泊水量变化ΔVlake(t)可用定积分求导：

(8)

式中：ti和ti+1为观测时段跨度，d；f(t)为湖泊水位和面积的关系式； dh为水位变化量。

假设湖泊水域范围随水位变化呈线性增加或者减少，在未考虑湖岸线地形、湖盆变形等因素的影响下，任一时段Tti+1-ti=ti+1-ti内，湖泊水量变化值ΔVlake由湖泊水位和湖泊水域面积共同决定[31]，ΔVlake可简化为：

(9)

2.4 误差分析

GLAS光斑脚点、SRTM3和GPS测量点的高程偏差，主要源于配准误差、地形误差和系统误差[32]。每个GPS测量点为35 m×35 m样方，共计32个，样方内中心点的高程测量5次，取5次高程测量值平均值作为控制点的高程真值；当样地平均坡度大于5°时，需考虑地形坡度对GLAS光斑脚点高程的偏差[33]。绝对误差Ea可用于描述GLAS脚点与GPS测量点、SRTM3高程之间的误差大小，即：

Ea=|zsat,i-zg,i|

(10)

式中：zsat,i为ICESat卫星下轨迹上第i个GLAS脚点高程值，m；zg,i为第i个GPS测量点或者SRTM3像元的高程值，m。

作为评价地理空间数据质量的重要指标之一，均方根误差(root mean square error，RMSE)亦称标准误差，被用于数字地图和地理空间数据点的位置精度评价[34]，湖泊日均水位估计的均方根误差RMSEz为：

(11)

式中：n为湖泊水域内GLAS光斑的数量。

标准化均方根误差(normalized root mean square error，NRMSE)用于分析湖泊水域GLAS光斑脚点高程与SRTM3像元高程、GPS测量点高程的偏差，标准化均方根误差NRMSEn为：

(12)

式中：zgps,max为GPS测量点(样方)的高程最大值，m；zgps,min为GPS测量点(样方)的高程最小值，m。

3 结果与分析

3.1 误差分析与验证

为了描述GLAS光斑脚点高程与SRTM3像元高程的误差大小以及空间异质性，本文随机从青海湖湖区周边陆地上N1区选取49个、N2区选取94个、S1区选取71个、S2区选取202个GLAS光斑，做GLAS光斑脚点高程与SRTM3高程的标准误差分析，结果见图3。由图3可知，N1、N2、S1和S2区中GLAS光斑脚点与SRTM3像元高程的标准误差平均值分别为0.19、0.12、0.17和0.14 m。

为了评估GLAS光斑脚点、SRTM3像元和GPS测量点高程之间的绝对误差，从湖泊水域边界向外2 km的缓冲区内，选择32个陆地GLAS光斑脚点与SRTM3像元、GPS测量点高程做绝对误差分析和关系式拟合，结果见图4。由图4可知，在陆地上，GLAS脚点高程与SRTM3高程的绝对误差均值为0.26 m，最大值为0.78 m，复相关系数为0.997 6；GLAS脚点高程与GPS测量点高程的绝对误差均值为0.14 m，最大值为0.46 m，复相关系数为0.967 6，且存在线性相关。GLAS光斑脚点与地表GPS测量点或者SRTM3高程的绝对误差，主要受到地形坡度、粗糙度、地表覆被、卫星姿态、透射深度等因素影响[35-36]，由于湖区周边地势较为平坦，植被覆盖度较低或者为裸地，对GLAS光斑脚点的高程误差影响较小，可用经验关系式修正。

为评估LEGOS水位和ICESat估计水位的差异，选取2003-2009年青海湖LEGOS日均水位值与ICESat-GLAS的日均水位估计值进行相关性分析，结果见图5。由图5可知，青海湖LEGOS日均水位值与ICESat-GLAS的日均水位估计值呈线性正相关，复相关系数R2为0.879 9，通过显著性水平p<0.05检验。

下社水文站水位观测时段为每年5-10月份，为了检验ICESat的当日估计水位平均值与下社站日均水位观测值的差异，选取2003-2009年期间5-10月份的湖泊水域GLAS光斑，求出湖泊的日均水位。ICESat估计的当日水位平均值与下社站水位观测值的Pearson相关系数为0.865，通过显著性水平p<0.01检验，二者存在系统性偏差，偏差均值为(1.08±0.22) m，见图6。

Landsat影像提取的青海湖水域面积与ICESat-GLAS的水位估计值相关性分析见图7，下社水文站年均水位与当年湖泊水域面积相关性分析见图8。由图7可知，两者呈线性正相关，复相关系数R2为0.789 1，通过显著性水平p<0.05检验。当青海湖湖泊水位持续上升较快时，青海湖大湖区与尕海、错果融为一体，青海湖湖泊水域面积增加较快，导致水位和水域面积时间匹配不一致，出现拟合异常点，偏离拟合线(见图8)。从湖泊水域面积和ICESat的年均水位估计值、年均水位观测值的拟合效果来看，水位实测值和湖泊水域面积的拟合曲线的相关性优于ICESat的估计水位与水域面积的拟合曲线，GLAS回波信号在水体表面的透射深度、波浪形状、重力位变化的影响使GLAS脚点高程产生一定随机误差。

图3 各GLAS光斑脚点高程与SRTM3高程的标准误差 图4 GLAS光斑脚点高程与SRTM3、GPS测量点高程

图5 2003-2009年青海湖ICESat估计的日均水位与LEGOS日均水位相关性分析 图6 2003-2009年青海湖ICESat估计的日均水位与下社水文站日均水位的差异

图7 Landsat影像提取的青海湖水域面积与ICESat-GLAS的水位估计值相关性分析 图8 下社水文站年均水位与当年青海湖水域面积相关性分析

3.2 湖泊水位、面积和水量的时变曲线

利用公式(3)、(5)和(9)分别求出1988-2018年青海湖的日均水位、水域面积和水量变化，结果见图9。借助2003-2009年期间ICESat卫星在轨运行期间的GLA14数据，建立湖泊水位-面积和水位-面积-水量绳套关系，推算出非观测时段内湖泊日均水位、水量变化。由图9可见，ICESat估计的1988-2018年青海湖水位与LEGOS水位、水域面积和水量变化均呈增加趋势，曲线峰谷节点较为一致，可计算出其Pearson相关系数分别为0.879 9、0.891 0、0.861 6。

ICESat卫星仅有18个工作周期，且工作周期间隔12～55 d，观测时段连续性不足，监测时段内湖泊水位的偏差和湖泊水域面积估计误差，会导致湖泊水位、水量变化值与LEGOS的水位和水域面积值存在系统性偏差和时间错位。因此，应用星载测高卫星资料估计较大面积的湖泊水量变化时，需要与水文过程、湖泊水动力模型相结合，减少湖泊水位估计中的随机误差。

图9 1988-2018年青海湖的日均水位、水域面积和水量变化曲线

3.3 湖泊水位、面积、水量变化的年际与年内变化特征

作为青藏高原地区的内陆湖，青海湖湖泊水量主要由降水、冰雪融水和冻土冰融水等补给[37]，山区降水量和冰雪融水量的变化会直接影响湖泊水域面积和水位。1988-2018年青海湖湖泊年均水位、湖泊水域面积和湖泊水量变化曲线见图10，图10中LEGOS的水位、面积和水量变化资料的时段为1995-2018年，与Landsat提取的水域面积在2004年左右存在偏差，是由于文中未统计海晏湾、沙岛湖等水域面积所致。由图10可知，1988-2018年青海湖湖泊年均水位、面积和水量变化，总体上呈增加趋势。青海湖年均水位由1988年的(3 194.19±0.15) m上升至2018年的(3 196.13±0.18) m，增幅为(1.93±0.22) m，见图10(a)；青海湖年均面积由1988年的(4 282.25±8.5) km2变为2018年的(4 480.0±12) km2，增加了(197.75±6.3) km2，见图10(b)；与1988年的湖泊水量相比，2018年青海湖湖泊水量增加了(8.93±0.12) km3，见图10(c)。

1988-2018年间青海湖的年均水位有两个显著变化：(1)1988-2004年期间，青海湖年均水位由1988年的(3 194.19±0.15) m下降到2004年的(3 193.0±0.16) m，年平均水位降低了(1.19±0.15) m；(2)2004-2018年期间，青海湖年均水位以0.19 m/a的速率上升，2018年年均水位比2004年增加了(2.69±0.25) m，湖泊水域面积增加了(233.0±7.3) km2。2018年下社水文站和LEGOS的年均水位比2004年的年均水位分别增加了2.55 m、(2.61±0.14) m，湖泊水域面积分别增加了(250.2±8.1) km2、(215.7±7.0) km2。

选取观测时段连续性较好的2010年青海湖湖泊水位、面积和水量变化，描述青海湖湖泊水位、面积和水量年内变化特征，结果见图11。由图11可看出，2010年1-4月期间，湖泊水量补给较少，湖泊水位降低，湖泊水量减少；同时，湖泊水面会有封冻、风吹雪或者浮冰[37]，湖冰或湖面积雪以升华为主，引起湖泊水面高程降低，湖泊水域面积萎缩。在5-6月期间，山区积雪和冻土层的地下冰开始融化，补给湖泊水量，但湖泊水域内的湖冰也开始融化，水位会略有下降。7-10月期间，受到山区降水、冰雪融水的补给使湖泊水量增加，湖泊水位逐渐回升。11-12月份，山区开始降雪，降雨和冻土层的地下冰融水量逐渐减少，湖泊水位逐渐下降，湖泊水面开始封冻。

图10 1988-2018年青海湖湖泊年均水位、湖泊水域面积和湖泊水量变化曲线 图11 2010年青海湖湖泊水位、面积和水量变化曲线

4 结 论

(1)在平地上，GLAS光斑脚点高程与SRTM3高程的标准误差呈线性相关，绝对误差最大值为0.78 m，SRTM3高程产品的C波段透射深度会使GLAS脚点高程的绝对误差增大；GLAS脚点高程与GPS测量点高程的绝对误差均值为0.14 m，绝对误差最大值为0.46 m，主要受地形坡度、粗糙度、地表覆被、卫星姿态等影响。ICESat卫星观测时段的不连续性和GLAS光斑空间分布的异质性，使得湖泊水位估计值与下社站观测值存在一定系统性偏差。

(2)1988-2018年的青海湖湖泊年均水位、水域面积、湖泊水量变化呈增加趋势，但在2004年湖泊年均水位达到最低值(3 193.0±0.16) m。2004-2018年湖泊年均水位持续上升，2018年年均水位已达到(3 196.13±0.18) m。青海湖水域面积也由1988年的(4 282.25±8.5) km2增加到2018年的(4 480.0±12) km2。

(3)湖泊水量作为湖泊水文观测的重要参数，当缺少湖泊面积或者水位观测资料时，ICESat测高数据可较好地估计湖泊水位和水量的变化，为青藏高原高山湖泊变化数据集重建提供理论参考和技术支持。虽然青藏高原高山湖泊的水量波动受到人为活动干扰较少，但是湖面蒸发、冰雪融水、山区降水、冻土退化等对湖泊水量波动的作用不可忽视。