Sentinel-2和ICESat-2密集时序数据反演浅平型内陆湖泊水量变化

吴浩儒,李均力,包安明,张久丹,马英莲

1. 中国科学院新疆生态与地理研究所荒漠与绿洲生态国家重点实验室,新疆 乌鲁木齐 830011; 2. 中国科学院大学,北京 100049; 3. 新疆遥感与地理信息系统应用重点实验室,新疆 乌鲁木齐 830011; 4. 新疆维吾尔自治区测绘科学研究院,新疆 乌鲁木齐 830002

湖泊水量是湖泊水情的关键因素,相对于面积和水位,水量更能客观地指示湖泊的时空变化[1]。干旱区内陆尾闾湖是封闭型流域水资源的最终汇集点,其水量变化能够客观反映流域的水量平衡过程[2]。然而,大部分干旱区内陆湖泊地处偏远,水文基础资料稀缺,尚未有湖盆地形及水量的公开信息[3],重建无资料地区内陆湖泊准确的水量变化可为分析内陆湖泊变化的机理、科学配置干旱区水资源及制定河湖生态治理策略提供科学依据[4]。

无资料或缺资料地区湖泊水量变化提取方法目前主要有基于地形测量和基于遥感的方法。其中,湖盆地形测量多采用船载回声测深法[5]和RTK地形测量法[6],这种方法测量精度高,但作业成本高、测量效率低,对于水深变动大、季节性变化显著的尾闾湖,开展相关测量工作的难度更大。遥感技术具有大尺度、周期性及全天候的观测特性,目前已成为获取内陆水体水文信息的重要手段[7],相关文献对光谱特征法、重力卫星法及水位-面积关系法等基于遥感的水量反演方法进行了归纳和总结[8-9]。其中,光谱特征法通过建立湖泊水深与水体光谱特征的经验关系进行水深反演[10-12],适用于水质良好的浅水湖;重力卫星法通过测量研究区域的重力变化来估算水量变化[13-15],适用于大尺度、大区域内陆流域的水量变化,如塔里木河流域或咸海流域;湖泊水位-面积关系法主要利用测高卫星和光学卫星分别获取水位、面积信息,建立水位、面积与湖泊水量变化之间的关系,并根据面积或水位推算湖泊水量或水量变化[16-19],适用于有时间一致的水位和水面积监测数据的湖泊。

内陆干旱区尾闾湖泊是河流流至广袤剥蚀盆地中心潴积而成,在长期湖泊沉积和强烈的风生湖流作用下,通常湖水较浅且湖底地形平坦,湖盆的最大深度往往只有4～5 m,湖岸线极其复杂[20]。这种湖泊的变化主要表现为水面变化迅速而水位变化微小,现有方法很难根据水位-水面关系来确定水量变化;而且在多风条件下,湖水中泥沙和悬浮物较多导致水体光谱特征不稳定,光谱特征法反演水深的精度变动较大[21]。近年来,文献[22]利用湖泊水面的等高特性,获取不同时相湖泊水面的激光雷达测高信息,形成若干等高线,并构建湖盆地形,在边界规则、地形渐进变化的自然湖泊上取得了较好的效果[23-24]。

随着遥感技术的不断发展,卫星遥感分辨率、覆盖周期和测高精度都有了长足的进步,其中Sentinel-2光学遥感卫星可实现5 d的重访周期[25],ICESat-2的激光测高仪发射的激光脉冲可获取沿轨道方向间隔0.7 m、光斑直径为17 m的光子信号,测高精度达到了2～3 cm[26-27],可以满足微小地形变化的测量要求。本文结合以上两种卫星的优势,以塔里木河干流尾闾湖台特玛湖为例,提出一种基于湖盆地形重建的水量变化估算方法。该方法通过湖盆内多时相ICESat-2非水面高程点和多时相Sentinel-2水面线等高点构建高精度湖盆地形信息,并在此基础上利用时序面积信息重建水量变化时序。

1 研究方法

湖盆地形重建主要技术流程如图1所示,主要包括多时序水面边界线的提取、多时序湖盆激光雷达高程点提取,水面边界线的高程内插、地形构建和水量时序重建等。其中构建高精度的湖盆地形是获取浅平型内陆湖泊水文参数的关键。单光子计数激光雷达(photon-counting LiDAR,PCL)的光子高程点可获取平坦地面高精度的高程信息,但由于轨道间距较大,无法全面覆盖湖盆,单靠湖盆内的激光雷达高程点无法构建完整的湖盆地形信息。浅平型内陆湖泊在上游河流来水时面积迅速扩张,河流断流时水面迅速消退乃至干涸。特别是在水面消退期间,湖泊分裂成众多小水面,且水体流动少、水面相对静止,可以近似将每个小水面看作一个等高面。为此,本文提出改进的湖盆反演方法,利用湖盆内无水区域的激光光子高程点对处于退水期的多时相水面边界线进行高程插值,不同水面边界线的高程点加密了湖盆内的高程点。在此基础上,利用加密的高程点构建TIN和数字高程模型,并采用精确的湖盆地形和多时相水面边界重建湖泊水量变化。

图1 浅平型湖盆的地形构建技术流程Fig.1 Flowchart of basin topographic retrieval for shallow-flat lakes

1.1 湖泊水面边界线提取

干旱区内陆湖泊的水陆边界明显,水体指数(normal different water index,NDWI)[28]是区分陆地和水面的一种主要方法。为此,本文利用光学遥感时序影像生产的NDWI图,采用“全局-局部”自适应阈值的方法[29]提取湖泊边界的时序信息。该方法提取的水面时序信息的生产者精度在96.42%～99.73%之间。

1.2 PCL光子高程点提取

PCL激光脉冲信号从发射到回波接收的过程中会受到太阳辐射、大气散射等影响,形成较多的噪声,如何去除噪声光子会直接影响高程信息的反演精度[30]。为此,本文提出一种改进的密度自适应的信号光子检测算法来区分信号与噪声(图2)。其主要原理是激光光子到达反射面并回波时,信号光子沿反射面的分布密度远高于噪声光子[31]。本文的信号与噪声光子检测方法如下。

图2 PCL光子高程信息提取Fig.2 PCL photons elevation information extraction

(1) 建立椭圆邻域计算光子密度。对每个激光光子m,根据光子信号与噪声的分布规律[32],采用长轴a=10 m、短轴b=1 m、倾角θ=0°构建椭圆邻域范围(图2(b)),根据以下公式判断光子点n是否在中心光子点m的椭圆邻域内

dx=(xm-xn)cosθ+(ym-yn)sinθ

(1)

dy=(ym-yn)cosθ-(xm-xn)sinθ

(2)

(3)

式中,a、b和θ分别表示椭圆邻域的长轴、短轴和倾角;xm和ym分别为中心光子点m的纬度和海拔;xn和yn分别为光子点n的纬度和海拔。如果dist(m,n)<1则说明光子n为中心光子点m邻域范围内的光子点。

根据邻域内的所有光子计算中心光子点m的光子密度如下

(4)

(5)

(6)

式中,Pm为中心光子点m的邻域点集;N为Pm中光子点数量;Wn是采用高斯核计算的椭圆邻域中每个光子点的权重[33];δ是高斯核的标准差;Dm为中心光子点m的点密度,由邻域内所有点的高斯权重求和得到。

(2) 计算光子最大点密度。顺时针每旋转15°计算一次中心光子密度,最后取得到的最大密度作为中心光子的点密度。通过自适应调节椭圆邻域的倾角θ,可以使信号光子的椭圆邻域减小地形变化的影响,捕捉到更多密集光子,更容易与噪声光子进行区分(图2(a)和图2(c))。

(3) 确定信号与噪声光子的阈值。当所有光子的最大点密度计算完毕,采用OTSU阈值法[34]计算信号光子和噪声光子密度之间的阈值,从而提取信号光子。

(4) 基于信号点生成地形剖面线。根据提取单激光束所有的信号光子点,采用最小二乘回归方法拟合出该轨迹下的地形剖面线,对信号光子较少的区域进行加密。

在计算PCL的地面光子高程信息后,利用与之时间对应的湖泊图层,提取湖盆范围内、湖面边界外的高程点数据,并将其高程基准面转换为大地水准面,作为湖盆内无水区的高程点。对不同时相的所有高程点进行提取,最后合并为覆盖湖盆的PCL光子高程点。

1.3 湖泊水面边界线高程内插

在获取湖盆内无水区的PCL光子高程信息点后,还需要确定不同时相湖面边界线的高程值,从而加密湖盆内的PCL光子高程信息点。本文采用反距离加权(inverse distance weighted,IDW)[35]的内插方法,用PCL光子点的高程值内插湖面边界线上点的高程值,进而获取湖边界线高程值,具体方法如图3所示。对于每条湖面边界线都会有若干条PCL光子轨迹与其相交(图3(a)),考虑到台特玛湖湖盆较为平坦,且PCL光子轨迹上点较为密集,可以在轨迹上选取交点两侧各4个距离最近的湖盆高程点进行反距离加权插值得到交点高程(图3(b))。对于每条湖面边界线,可以求出边界交点的高程值,选定远离平均值3倍标准差的值为异常值并将其剔除[36],然后取平均作为水面边界的高程。

图3 内插湖泊边界高程Fig.3 Altitude interpolation of lake boundaries

1.4 湖盆地形构建

多时相的湖泊水面边界高程点既能加密湖盆内的PCL光子高程点,水面等深线也能更好模拟湖盆的地形地貌特征。为此,利用加密后的湖盆高程点构建不规则三角网,并在此基础上利用Kriging插值生成湖盆的数字高程模型(DEM),DEM的格网大小与提取水面边界线的遥感影像分辨率一致。

1.5 湖泊水量时序重建

浅平型湖泊的季节性变化显著,湖泊往往分裂为众多的小湖,且每个小湖泊的水面高程存在一定差异,常用的“水位-面积-水量”方法不能适用这种情况,需要分别对每个小湖分别计算水量,并累加为湖泊的总水量。为此,每个时相的水面对应的湖泊水量计算方法为

(7)

式中,Nt为时相t湖泊的小水面数量;Ai为第i个小水面的像素数量;Hi为第i个小水面的高程;Hi,j为第i个小水面、第j个像素的高程;Si,j为第i个小水面、第j个像素的面积。对每个时相的湖泊水面线,根据式(7)分别计算湖泊的水量,最后形成湖泊的水量时序信息。

2 研究区与数据

2.1 研究区概况

本文选择塔里木河干流下游的尾闾湖——台特玛湖作为试验区构建湖盆并重建湖泊水量,并且以最大湖面为基准,根据地势起伏程度向外扩张3～5 km作为湖盆范围,如图4所示。台特玛湖位于新疆维吾尔自治区若羌县的北部,是塔里木河干流和车尔臣河交汇形成的尾闾湖。台特玛湖的西部是车尔臣河下游的康拉克湖群,南部的若羌河与米兰河是间歇河流,仅在丰水期有地表径流汇入湖区。湖泊东部为库木塔格沙漠,历曾有河道与罗布泊相连,后因河水量减少及风沙阻塞河道,台特玛湖开始在车尔臣河与塔里木河交汇处潴水成湖,湖泊周边地形平坦,地势呈西高东低的特征。

2.2 数据源选择

湖泊水面线信息提取采用2016—2022年共132期无云和质量较好的Sentinel-2卫星遥感数据,覆盖了不同季节从湖盆无水到最大水面的各种情况。高程信息数据从冰、云和陆地高程卫星2号(the ice,cloud,and land elevation satellite-2,ICESat-2)的地形激光测高仪系统(advanced topographic laser altimeter system,ATLAS)数据提取,数据源为美国冰雪数据中心发布的ICESat-2/ATLAS全球地理定位光子数据产品ATL03[37]。ATL03产品是获取陆地表面高程的基础数据,记录了包括大量噪声光子在内的传感器接收的所有光子经纬度及椭球体高度信息[38]。本文挑选了28对时相相近的ICESat-2测高点数据和Sentinel-2影像用于湖盆高程点提取,详细信息见表1。

表1 选择的ICESat-2和Sentinel-2数据对应列表

3 试验与结果分析

3.1 试验过程

本文基于ICESat-2测高卫星数据和Sentinel-2遥感影像数据,提取了密集的湖盆高程点,构建了较为精确的湖盆DEM,具体过程如下:①分别提取多时相Sentinel-2水面信息和ICESat-2高程点信息,28个时相的激光光子经过去噪后的有效光子数见表1,共10 4653个光子高程点;②根据与ICESat-2时相对应湖泊水面,剔除湖泊水面内的ICESat-2高程点,从而得到不同时相湖盆内无水区的高程点共82 050个;③对所有不同时相Sentinel-2提取的水面边界线,根据水面等高特性[39],采用ICESat-2高程点内插出每个单独水面线的高程,得到967 689个高程点;④得到湖盆内的高程点共计1 049 739个高程点。经过水面线高程点的加密,不仅湖盆内的高程点数增加了11.8倍,而且这些水面线还可以作为地形构建的高程约束线。经过加密后的湖盆高程点构建的湖盆DEM如图5(b)所示。

图5 构建台特玛湖DEMFig.5 Construction of Taitema lake DEM

在获取精确的湖盆DEM后,结合2016—2022年台特玛湖的水面边界,根据式(7)可以计算台特玛湖的水量时序,如图6所示。由图6可知,湖泊的水量变化与面积变化是高度一致的,湖泊在的面积在2017—2020年处于面积和水量较大的阶段,这个阶段湖泊在1—2月水量达到最大,其中2019年1月19号水量达到最大值,为0.894 km3。而在2016年和2021年以后,水量在3—4月达到最大,最大面积在0.189～0.234 km3之间,仅为水量高峰期的1/4。湖泊水量最小的时段均在7—8月,在2016—2022年几乎没有变化。

图6 台特玛湖水量时序重建Fig.6 Time-series reconstruction of water volumes of the Taitema lake

3.2 精度评价与误差分析

由于台特玛湖尚未有公开的湖盆地形或水量数据,很难直接验证湖泊水量时序的精度。由于湖泊在不同水深下面积变化极大,且分裂成数量众多的小水面,同时由图5可知,台特玛湖的湖盆地形较为复杂, 这种情况下湖泊水量的计算很难用面积-水位关系方法计算,而湖泊水量的精度主要取决于湖盆地形的精度。为此,本文采用较为常用的检查点法和等高线回放方法[40]评估湖盆DEM的精度。

3.2.1 检查点法精度验证

检查点法是在DEM上设置一些检查点作为高程真值,并从DEM上提取检查点的内插高程,比较各个检查点的误差,然后算出中误差用贝塞尔公式[41]计算高程中误差,如果中误差的离散度小,则表明DEM的精度高。这里选择Sentinel-2影像2018年5月26日提取的湖泊水面线,分别从湖盆DEM上内插得到每条水面线的高程点信息,共得到185 493个检查点。由于湖泊的每个小水面可以看作是等高的,分别计算每个水面线的平均高程作为真值,进而计算中误差。由于中误差只有在服从正态分布时才能反映数据的精度,本文先采用正态QQ图对误差分布进行了正态检验(图7(a))。

图7 检查点的误差分布Fig.7 Error distribution of the check points

由图7可知,正态QQ图围绕对角线呈近似直线分布(R2=0.984),检查点误差服从正态分布统计特征[42]。计算这些水面检查点高程的平均中误差RMSE=0.103 m,从误差分布图来看,99.82%的误差在3倍中误差内。这表明同一水面线高程点的中误差离散度很小,残差的实际值与理论值越接近,模型的拟合效果越好。因此,从检查点法定精度验证表明DEM的误差满足极限误差的要求。

3.2.2 等高线回放法精度验证

等高线回放法[43]指的是对一条湖泊水面线,计算水面线点的平均高程,并根据该高程从湖盆DEM内插出等高线,并与实际的湖泊水面线的点一一对应,通过计算两者的距离差来反映DEM的误差。本文以2018年5月26日Sentinel-2影像提取的水面边界线为例,首先计算每条水面线的平均高程值,并以此值内插出等高线,这里选取6个大小不一的水面及对应的等高线(图8),分别计算水面与等高面的面积、对应点的最大和平均距离,来反映湖盆DEM与水面线的相似度(表2)。

表2 湖盆等高线与水面边界的误差对比

图8 用于等高线回放法验证的水面Fig.8 The water samples for accuracy assessment by the contour playback method

结果显示图8中的6个湖面线与生成的等高线空间吻合度很高。等高面与水面的面积误差最大仅为1.564%,模拟的等高面的误差已经低于湖泊水面提取的误差,这表明湖盆DEM精度能很好地描述湖泊水面的精细变化。而等高线与水面边界线的最大距离和最大平均间距分别为39.525及9.794 m,这表明每条模拟的等高线的平均间距在1个像素(10 m)内,最大误差不超过4个像素。结合图9可知,等高线与水面边界线距离小于10、20及30 m的像素点比例分别大于71.5%、94.4%及99.8%,这表明利用DEM生成都等高线与水面线具有非常高的一致性,本文方法构建的DEM所提取的等高线能够很好地吻合实际的湖面边界线。

图9 等高线与湖面边界间距的累积百分比Fig.9 Cumulative frequency of distance between contour and water boundary

综上所述,本文构建的湖盆地形中误差仅为0.103 m,这表明湖盆的单点高程精度很高,经过水面线加密后的ICESat-2高程点能更好地体现。而且湖盆DEM提取的等高线能与实际的湖面线吻合度高,说明湖盆DEM已经能模拟浅平型湖盆地形的精细变化。

3.3 讨 论

干旱区内陆尾闾湖泊经过长时间的沉积和演化,湖盆地形平坦,湖泊水深较浅,虽湖泊总水量较少,但常年保有最小生态水面是维持湖泊周边生态系统健康的关键所在。本文就是为了精确计算缺资料浅平型内陆湖泊的水量变化,进而为湖泊生态补水提供科学依据,保障水资源的高效利用。由图5的台特玛湖湖盆地形可以看出,台特玛湖湖盆地形整体上西高东低,整个湖盆地区东西横跨约96 km,而海拔仅变化了约15 m,地形起伏非常细微。这类浅平型的湖泊的水面积变化通常较为明显,而水量和水位变化相对非常微小。由图6构建的台特玛湖水量时序可以看出,台特玛湖在2016—2019年间水面积和水量均呈动态增加的趋势,年均变化率分别为144.64 km2/a和0.28 km3/a,在2019—2022年间水面积和水量均呈动态减少的趋势,年均变化率分别为150.48 km2/a和0.29 km3/a。

传统的基于遥感的“面积-高程-库容”方法在应用于这类湖盆水量反演时,由于SRTM(shuttle radar topography mission)高程精度在2～10 m范围内[44-45],难以获得精确的水面-水位关系。利用ICESat-2高精度的激光光子高程点构建缺资料地区的湖盆地形,并利用水面等高的特性作为地形约束线加密湖盆内的高程点,湖盆地形精度可达0.1 m左右,能够描述湖盆内地形的微小起伏,对于水面周期性、季节性变化显著的湖泊,可大大提高水量时序反演的精度。然而,该方法的前提条件是,ICESat-2过境的激光光子束能够探测到湖泊从水面最大到干涸的全过程,对于横向距离小于3 km(ICESat-2轨道间距),以及水面变化不大的湖泊该方法难以适用。

4 结 论

针对干旱区浅平型尾闾湖水量难以估算的问题,本文基于湖泊水面边界等高的特性,提出利用湖泊水面线高程加密ICESat-2地面高程点,构建的湖盆地形结合水面时序数据模拟湖盆水量变化时序,主要思路是通过密集的湖盆内无水区高程点构建湖盆地形,并结合时序面积信息重建湖泊水量时序,重点解决缺资料地区湖泊水文参数的信息反演。以台特玛湖为主要试验区,重建2016—2022年台特玛湖湖泊水量变化,验证了湖盆地形的反演精度。研究得出以下结论。

(1) 本文提出了基于改进湖盆反演技术的湖泊水量时序重建方法,能够提取精度要求更高的浅平型湖泊地形构建,适宜于水位变化不敏感的湖泊水量估算。

(2) 改进湖盆反演技术构建的湖盆地形信息的高程中误差仅为0.103 m,而且DEM提取的等高线与实际水面线高度吻合,这表明改进湖盆反演技术满足浅平尾闾湖的高精度水量信息的反演。

(3) 基于改进湖盆反演技术的水量时序重建方法能够适用于干旱区浅平型的内陆尾闾湖,也可以推广到全球范围内水面变化显著且有卫星测高数据过境的湖泊,为无资料浅平型湖泊水量信息的时序重建提供一种思路。