ICESat-2/ATLAS测高数据在青海湖湖泊水位估计中的应用

吴红波， 王宁练， 郭忠明

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

1 研究背景

随着先进对地观测技术的发展和应用，星载遥感技术已成为获取地球各圈层及交互过程信息的重要手段之一[1-2]。其中，星载激光雷达(light detection and ranging，LiDAR)测高技术将激光雷达系统、全球定位系统(global positioning system，GPS)、卫星通信-测控系统、惯性导航系统融为一体[3]，具有覆盖范围广、测量精度高、时效性较好、多时相等优势，为准确获取地球表面各圈层物质质量变化及迁移以及海平面、地表水体水位和地物三维空间信息提供了有效的技术支撑[4]。2003年1月13日全球首颗冰、云和陆地高程卫星ICESat-1(ice，cloud and land elevation satellite-1)发射成功，其搭载的地球科学激光测高系统GLAS(geoscience laser altimeter system)能够在全球或者较大空间范围内全天候、不间断地观测地球表面空间和结构[5]，主要用于测量冰川和冰盖质量变化、气溶胶厚度、海冰厚度、地表水体水位、海平面变化、植被冠层高度和生物量以及三维地形信息[6]。ICESat-1/GLAS在2003至2009年运行期间，为林学、冰冻圈科学、全球气候与环境变化、大地测量学、海洋科学、水文模型与循环等研究提供了连续、周期性、高精度的历史记录数据。

为了克服ICESat-1卫星航迹间隔较大、数据覆盖和空间代表性不足等问题，ICESat-2卫星搭载了先进地形激光测高仪系统ATLAS(advanced topographic laser altimeter system)，其中，光子计数激光测高仪(photon-counting laser altimeter)每秒发射10 000个激光脉冲来进行地表高度测量，测量精度约为5 mm，将继续观测山地冰川和极地冰盖高度变化、陆地和植被高度、内陆水体水位、海平面、海冰厚度、云层厚度等空间信息[7]。ICESat-2卫星提供的测量数据填补了ICESat-1与“冰桥行动”[8]、欧空局的CryoSat-2之间的空隙[9]。国内外研究人员利用ICESat-2产品估计湖泊水位变化[10-11]、山地冰川高程变化[12]、植被垂直结构[13]、海冰厚度等[14-15]的研究成果已有报道，研究结果侧重于ATLAS产品的基础理论、应用和估测精度验证。Yuan等[16]选取了中国境内30座水库和水域面积大于10 km2的湖泊，利用ICESat-2测高资料提取水位信息并与观测值进行了精度检验，结果表明，ATLAS测高仪的水位估计值的相对误差为0.06 m，每条轨道上光斑高程的标准偏差小于0.17 m；Kwok等[17]利用ICESat-2和CryoSat-2资料对2018年10月14日至2019年4月底北极海冰上的积雪厚度进行了估算，讨论了积雪深度估计的不确定性；Neuenschwander等[18]利用ATL08产品对ICESat-2卫星下轨道内的三维地形重建和地表植被冠层进行了定量评估，表明光斑脚点位置存在5 m水平偏移，地表和植被冠层高度估计的均方根误差分别为0.85和3.20 m；Zhang等[19]利用多光束高度计实验激光雷达MABEL(multiple altimeter beam experimental lidar)数据，估算出美国北卡罗莱纳州东海岸区域的海平面的背景噪声率，不仅填补了水面噪声光子理论的空白，也提供了一种快速有效的地表分类方法；Wang等[20]评估了ICESat-2卫星数据在地表高程变化中的估测性能，得出ICESat-2/ATLAS与机载LiDAR系统反演的地面标高的均方根误差分别为-0.61和1.96 m，在森林、冻土和裸地上的地面标高的RMSE(root mean square error)值分别为-0.64、-0.61和-0.59 m；Brunt等[21]利用全球导航卫星系统的实时动态高程数据、ICESat-2卫星的ATL03和ATL06产品提取了南极冰盖高度变化值，表明GNSS(global navigation satellite system)高程数据的绝对误差为5.6 cm，ATL03产品的绝对误差小于5 cm，而ATL06产品的绝对误差小于3 cm。虽然研究人员对ATLAS产品在理论方法、测量精度、实地验证方面开展了一些尝试和应用[22]，但是基于ATL13产品重建湖泊水位时变序列的相关研究较少，亟待借助连续实地水位观测，分析湖泊水位估计精度的不确定因素(如湖岸线、波浪、湖流、浮冰等)，探讨ATL13产品重构日、月、季节和年际湖泊水位时变序列的可行性。

为了检验ATL13产品估计湖泊水位的估计精度和时空差异，文中选择青藏高原地区的青海湖作为研究对象，基于2018年10月31日至2019年11月8日期间的ATL13测高数据提取青海湖水域的水面高程，借助风浪、水位、气象观测等资料，分析ATL13产品中6束激光脉冲的湖泊水位估计精度。考虑到大风天气下的湖泊水域会有较大波浪，选择刚察站风力等级0～2级和3～5级天气条件下的ATL13产品中的水面高程记录，分析水面波浪对湖泊水位估计精度的影响。由于青海湖湖泊水域范围较大，提取同一时期的6束脉冲所记录的湖泊水位，分析湖泊水面高程剖面上的空间差异，并重构了湖泊水位日均和月均变化特征，其结果可为内陆湖泊、水库、河流水位变化研究提供理论参考和技术借鉴。

2 数据来源与研究方法

2.1 数据来源

2.1.1 ATLAS测高数据 为了继续执行ICESat-1卫星的测量任务，2018年9月15日ICESat-2卫星在美国加利福尼亚范登堡空军基地成功升空，ICESat-2卫星轨道运行高度为500 km，对地飞行速度为6.76 km/s，回归周期为91 d，轨道倾角为92°，共预设1 387个轨道。搭载在ICESat-2卫星上的ATLAS激光测高仪重约155 kg，采用绿色波段532 nm处激光脉冲和单光子敏感探测器来测量地表空间信息[23]。ATLAS测高仪使用3对光束(共6束)，每对光束之间间隔约3 km，间距约90 m，激光脉冲需要3.3 ms时长就能到达地球并返回。单束脉冲在地表形成的光斑直径为17 m，脉冲频率为10 kHz，沿升/降轨道方向光斑采样间隔为0.7 m，单束脉冲半峰全宽(full width at half maximum，FWHM)为1.6 ns。由于ATL13产品可用于湖泊、河流、水库、人工运河等水位估计[24]， 故采用2018年10月31日至2019年11月8日期间ATL13中的6束测高脉冲提取青海湖水域范围内水面高程(见图1)。ICESat-2卫星的6束激光脉冲在地表形成6 km宽幅的条带状区域，激光脉冲的地表覆盖范围和光斑数量的增加，可降低地表坡度和粗糙度所引起的高程偏差，缩小升/降轨道之间的空隙。

图1 青海湖流域及ATL13数据分布

Level 3A级ATL13产品用于青海湖湖泊水域表面瞬时水位提取，数据版本为V002，数据类型有shape矢量、csv、hdf、asccii、二进制等格式，可从美国国家冰雪数据中心(National Snow and Ice Data Center，NSIDC)获取，数据下载网址URL：https://snsidc.org/data/atl13/versions/2。每个数据包提取光斑脚点的经度、纬度、高程、过境日期、大地水准高度、水面高度、标准偏差等参数，青海湖水域内有效的光斑数量共计213 351个，见表1。

表1 2018年10月-2019年11月青海湖湖区水域ATL13产品光斑汇总

2.1.2 观测数据 青海湖年初、年末、年均水位和水量变化数据来源于《青海省水资源公报》和青海水利信息网(URL：http://slt.qinghai.gov.cn/)，其中水域面积包括青海湖大湖区、沙岛湖和海晏湾，不含尕海、果错。下社水文站的水位观测时段为当年5月份至10月份，用当日水位观测值与ICESat-2卫星的水位估计值进行对比验证。2018年10月至2019年11月青海湖湖泊水位资料从法国地球空间物理和海洋学研究实验室(Laboratoire d’Etudes en Géophysique etceanographie Spatiales，LEGOS)获取，下载网址URL：http://www.legos.obs-mip.fr/observations。刚察气象站的风速日值、时值记录资料来源于国家气象科学数据中心，数据下载网址URL：http://data.cma.cn/，资料内容包括最大风速、最小风速和平均风速等。2018年9-10月和2019年4-5月期间，在青海湖湖泊水域南缘和北缘(见图1)，使用厘米级差分GPS采集仪对湖泊水域边界和平地高程进行测量，并随机选取53个地面控制点用于ATL13产品光斑脚点高程的精度分析。湖泊水面波浪高度采用SBF3-1型波浪浮标遥测系统测量，测量范围为0.1～2.0 m，采样周期为15 s。

2.2 研究方法

首先，利用水体、流域边界等基础地理空间数据对ATL13产品中光斑中心经度、纬度信息做掩膜处理，提取青海湖湖区内ICESat-2卫星的星下航迹，并做光子信号延迟、地球物理修正等；其次，根据ATLAS测高仪的6束激光脉冲在水面的高程初始值，对原始信号做背景噪声处理并提取每个数据包内有效信号的相关参数(波峰、波谷、起止时间、能量、频率等)；再次，对6束脉冲信号所对应的光斑高程进行修正(重力位异常、水准面高度、投影变形)，计算出修正后的光斑高程相关参数(最大值、最小值、标准差、高程差范围等)和湖泊水位瞬时水位；最后，剔除波浪、漂浮物、湖岸线、湿地植物等引起的水位异常值，利用湖泊水域内6束激光脉冲对应的有效光斑高程进行线性加权，求出湖泊水位平均值，技术路线见图2。

图2 技术路线

2.3 湖泊水位估计及误差分析

2.3.1 湖泊水位估计 ICESat-2卫星测高数据估计湖泊、水库、河流等地表水体的瞬时水位Hlake计算式如下：

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

Cst+Cpt+e

(1)

式中：Hsat为卫星在近极地轨道面上的飞行高度，km；Crange为ATLAS测高仪到地球表面或者物体表面的欧氏距离，km；Cdelay为电离层折射或者延迟误差，m；Cpressure为激光脉冲在大气层传输过程中，大气气压变化所引起的光子信号延迟，m；Cwet为激光脉冲穿过对流层时，大气湿度变化引起的延迟校正，m；Cst为地壳运动或者形变所引起的垂直高度异常值，m；Cpt为潮汐变化所引起的湖区水面高程误差，m；e为未考虑的其他不确定性偏差，m，包括激光脉冲的后效应[25]、背景噪声干扰等。

(2)

式中：halt,lake为2 160阶重力位模型Earth Gravitational Model 2008(EGM2008)模拟出的湖泊水面高程，m，格点大小为5′×5′；a为斜率；b为高程偏移常数[26]。

(3)

在任一时段ti+1-ti内，湖泊水域日均水位变化量ΔH为:

(4)

2.3.2 误差分析 目前，光子散射、高背景噪声和大气损耗将星载激光测高仪技术的扫描宽幅限制在几十公里以内[29]。光子散射引起的高程误差源于抽样误差、背景噪声干扰、复杂地表变化、第一光子偏压、大气前向散射和地表散射6个方面[30]。在湖泊水域的湖心区和边缘区，水面可能被视为一个粗糙的平面，ATLAS测高仪发射的光子在表面产生早期光子和晚期光子，并返回到系统接收器，返回的光子束在水面光斑内的高度具有近似高斯分布，那么，湖区水面的均方根误差σR等于脉冲传输、表面坡度和粗糙度引起的高度偏差的平方根：

(5)

式中，φ为激光脉冲与地表法线的夹角；c为单位m/s的光速值，取值为c=299 792 458；σtx为激光脉冲在时间尺度上的标准差，约为0.68 ns，相当于1.6 ns的半波宽；σbeam为激光脉冲在空间尺度上的标准差，约为4.25 m；R为激光脉冲的均方根偏差，m，一般指由地表粗糙度和地形坡度共同引起激光脉冲的高程偏移。

湖泊水面的有效波高Hs估计为：

Hs=4σh

(6)

式中：σh为湖泊水面的估计标准偏差，设定有效波高阈值为0.10 m。

ATL13产品光斑高程和GPS测量点的高程偏差主要源于配准误差、地形误差和系统误差[31]。每个GPS测量点为35 m×35 m矩形样方，共计53个，样方内中心点的高程测量5次，取5次高程测量值的平均值作为控制点的高程真值；当样地平均坡度大于20°时，需考虑ATLAS激光脉冲受坡度所引起光斑脚点高程的偏差[32]。ATL13产品中光斑脚点高程与GPS测量高程的绝对误差Ea为：

Ea=|zsat, i-zg, i|

(7)

式中：zsat, i为ICESat-2卫星下轨道上第i个光斑脚点的高程值，m；zg, i为第i个GPS测量点的高程值，m。

均方根误差RMSE可用于数字地图或者地理空间数据中点的二维、三维空间位置精度评价[33]。因此，湖泊水面光斑中心点的水位估计值的均方根误差RMSEz为：

(8)

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

采用标准化均方根误差(normalized root mean square error，NRMSE)分析湖泊水域ATLAS光斑脚点高程与厘米级差分GPS测量点高程的偏差，每个测量点重复采集15次。标准化均方根误差NRMSEn的计算式为：

(9)

式中：zgps,max和zgps,min分别为厘米级差分GPS采集仪所采集控制点的高程最大值和最小值，m。

3 结果与分析

3.1 不确定性分析与精度检验

为了分析ATL13产品的6束激光脉冲gt1L、gt1R、gt2L、gt2R、gt3L、gt3R对应湖泊水域的高程变化，沿着ICESat-2卫星上升轨道，由南向北依次提取湖泊水域内6束光斑高程值，见图3(a)、3(b)和3(c)。地处高海拔地区的青海湖流域常年多风，湖泊水面波浪会给湖泊水位估计带来不确定性误差，需要根据刚察气象站当日观测资料，选取风力等级0～2级、3～5级的天气，分析波浪对湖泊水面高程轮廓线的影响。在风力等级小于2级或者无风条件下，ATL13产品估计的湖泊水域轮廓线近似呈U形，靠近湖岸线的水域湖泊水位高于湖心区的水位，水位高差绝对值为0.3 m，见图3(a)；当风力等级较大时，青海湖湖区水面波浪较大，迎风岸的水面高程轮廓线的水位线较高，青海湖南缘水域的水位高于湖心区和北缘水域的水位，迎风岸和背风岸水位高差绝对值为0.5 m，见图3(b)。由图3(a)和3(b)可知，虽然青海湖湖泊水位轮廓线空间形状受湖底地形、湖流速率、湖岸线形状等影响，但当风力等级较大时，迎风岸和背风湖岸线处湖泊水位的绝对高差也会随着风力等级的增大而增加。

为了检验ATL13产品中在6 km幅宽区域内的湖泊水位变化的空间差异，选择风力等级为0～2级时的湖泊水面，分别提取6束激光脉冲的湖泊瞬时水位轮廓线，见图3(c)。由图3(c)可知，6束激光脉冲估计的水位值存在水位偏差，绝对偏差为0.08 m。为了降低ATLAS测高仪激光光斑的空间分布异质性引起的水位估计偏差，用ATL13产品中6束激光脉冲对应的瞬时水位值进行线性加权，估计出青海湖湖区水位的平均值。

综合图3(a)、3(b)和3(c)可知，ATL13产品中的6束激光脉冲估计水位的标准偏差均小于0.10 m；湖心区水位的估计偏差较小，标准偏差均小于0.06 m；因而湖心区的水位估计值的精度优于湖岸线区域的水位估计值的精度，湖心区的水位估计值具有良好的代表性，更能说明湖泊整个水域的水位变化和趋势。为了评估湖岸线周边水域瞬时水位受水面波浪的影响程度，文中提取距湖岸线0～1 km，1～2 km、2～3 km、3～4 km、4～5 km、5～6 km、6～7 km、7～8 km、8～9 km、9～10 km区间内的平均水位，见图3(d)。

由图3(d)可知，青海湖南、北缘水域的平均水位在9～10 km处趋于3 195.85 m，瞬时水位变化幅度小于0.05 m。受波浪和迎风湖岸线影响，湖泊水域的水位变化存在空间差异，青海湖南缘水域的平均水位随着距湖岸线距离的增加而下降速率较快，北缘水域的平均水位随着距湖岸线距离的增加而上升速率较慢。因此，为了降低波浪和湖岸线对湖泊水域平均水位估计的干扰，文中剔除距湖岸线0～9 km区间内的光斑后再估算湖泊瞬时水位的平均值更具有合理性。

图3 青海湖水域ATL13数据轮廓线和实测高程值对比

2018年9-10月和2019年4-5月期间，在青海湖南缘、北缘区域的水域内连续获取水位和陆地表面高程信息，选取2018年10月31日和2019年5月1日的ATL13产品中ATLAS光斑的水位估计值和水位观测值进行拟合分析，由于ICESat-2卫星过境和实地观测日期存在一定的时间差，采用湖岸线附近浮标水位实测值与ATLAS光斑的水位估计值进行线性拟合，结果见图4。由图4可知，2018年10月31日和2019年5月1日的ATLAS光斑的水位估计值与水位观测值均存在相关性，复相关系数分别为0.850 7和0.644 1，绝对误差分别为0.13和0.18 m，浮标观测的系统误差和记录时间滞后也会引起水位观测值与ATLAS光斑水位估计值之间绝对误差的增大。

图4 ATLAS光斑水位估计值与浮标水位观测值的拟合关系

3.2 湖泊水位日变化

为了描述LEGOS水位、ATL13产品的水位估计值和当日水位实测值，分别提取2018年10月至2019年11月初ATL13产品的6束激光脉冲gt1L、gt1R、gt2L、gt2R、gt3L、gt3R对应的瞬时水位，根据公式(2)求出当日估计水位平均值，并与LEGOS水位和下社水文站当日观测水位进行对比分析，见图5。分析图5可知，在研究时段内ATL13产品中6束激光光束gt1L、gt1R、gt2L、gt2R、gt3L、gt3R估计的日均水位分别上升了0.47±0.05 m、0.36±0.04 m、0.48±0.06 m、0.45±0.08 m、0.47±0.03 m、0.44±0.07 m；与2018年10月31日的日均水位相比，2019年11月8日的青海湖湖区的日均水位总体上升了0.45±0.06 m。ATL13产品中的6束激光光束估计同一区域的水面高程时，日均水位的绝对误差均小于0.07 m。青海湖的LEGOS水位从2018年10月26日的3 195.93±0.20 m上升到2019年11月27日的3 196.22±0.03 m，水位增加了0.29±0.20 m。青海湖下社水文站附近的日均水位观测值从2018年10月的3 195.71±0.12 m上升到2019年11月初的3 196.29±0.08 m，水位增加了0.58±0.10 m。通过ATL13产品的日均水位估计值、LEGOS日均水位和下社水文站水位观测值对比，发现ATL13产品的日均水位估计值与日均水位观测值、LEGOS日均水位的绝对误差分别为0.07、0.26 m，因此，ATL13产品的日均水位估计值与下社水文站水位观测值较为相关，具有良好一致性，而LEGOS水位的记录有数据质量、仪器测量精度、时间滞后等误差源，水位数据与实际水位记录偏差较大。

图5 2018年10月—2019年11月ATL13产品中6束激光光束估计的青海湖日均水位与观测水位和LEGOS水位对比

3.3 湖泊月均水位变化

青海湖水域面积较大、湖岸线较长，文中通过对ATL13产品中6束脉冲对应的瞬时水位进行线性加权，拟合出的青海湖湖区水域的月均水位变化曲线见图6。由图6可知，2018年10月青海湖的平均水位估计值为3 195.75 m，2019年10月的平均水位估计值为3 196.21 m，年内湖泊月均水位上升了0.46 m，与下社水文站月均水位观测值相差-0.06 m。ICESat-2卫星航迹之间的间隔和光斑空间分布的异质性导致水位高差和变化速率的差异，虽然已剔除波浪、湖岸线等影响，但在ATL13产品估计的青海湖湖区月均水位与水文站水位实测值仍存在一定偏差，这可能与河流补水、湖泊底盆、水深、湖流等动力学因素有关，需要借助水动力学模型分析湖泊水位动态变化过程中的水位误差。

图6 2018年10月-2019年10月ATL13产品估计青海湖月均水位变化

4 结 论

(1)湖岸线、波浪、ICESat-2卫星航迹间隙、湖泊水面异质性等因素会直接影响到星载测高数据ATL13产品中湖泊瞬时水位的估计精度。与ICESat-1卫星相比，ICESat-2卫星测高资料在空间分辨率、时间分辨率和数据覆盖范围方面已明显改进。6束激光脉冲同步观测约6 km宽幅的条带区域，地表光斑采样密度较大，消除了ICESat-2卫星轨道的间隙和地表坡度引起的垂直误差。与水位观测值相比，ATL13产品的湖泊水位估计的误差绝对值小于0.07 m，标准误差为0.18 m。

(2)综合2018年10月至2019年11月期间水位实测值、水文站观测值、LEGOS水位值和ATL13产品的青海湖水位估计值进行分析，ATL13产品估计的2018年10月青海湖月均水位为3 195.75 m，2019年11月的月均水位估计值为3 196.21 m，年内湖泊月均水位上升了0.46 m。此外，青海湖下社水文站附近的月均水位增加了0.58±0.10 m，而LEGOS的月均水位增加了0.29±0.20 m，这也与LEGOS原始数据源、数据质量和水位计算方法等有关。

(3)青海湖湖泊水域范围较大，湖岸线较长，为了减少湖面波浪、湖岸线对ATLAS测高资料估计湖区平均水位的不确定性，应尽量选择靠近湖心区的ATLAS光斑。ICESat-2卫星航迹之间的空间异质性会引起水位变化速率和水位高度差的变化，可借助线性加权法对6束激光脉冲的瞬时水位进行拟合，以降低ATL13产品估计湖泊水位的偏差。总之，随着ICESat-2卫星测高数据的历史累积和处理方法的完善，ATLAS测高产品在空间分布和测量精度方面，将为湖泊、水库、沼泽、河流水位变化提供可靠的数据源，也会为流域水量平衡、湖泊水动力学模型和水库库容变化等方面研究提供数据支持。