杨兴超， 郭斌*， 邢文雪， 项颖颖， 郭瑞莲， 徐莹， 季民
1 山东科技大学测绘与空间信息学院， 山东青岛 266590 2 泰安市林业保护发展中心， 山东泰安 271000 3 潍坊市峡山水库管理服务中心， 山东潍坊 261325
湖泊水库是国土资源的重要组成部分，具有调节河川径流、减轻洪涝灾害、维系流域生态平衡、提供生产生活用水等诸多功能(杨桂山等，2010；
Shu et al.，2020).水位是湖泊水库的重要特征，及时、精确地获取湖库水位对于政府决策、水资源管理、气候变化研究等具有重要的意义(袁翠，2017；
文京川等，2018).
传统获取湖库水位的方法主要是依靠水文站进行实地监测，该方法精度高，可以连续获取水位，但存在成本高、时效性差、数据难共享等问题，给水资源管理和相关科学研究带来了一定障碍(田山川等，2018).卫星测高技术的出现和迅猛发展为高精度、低成本、大尺度、周期性地监测湖库水位提供了新思路(付波霖等，2021).近年来，利用该技术监测湖库水位的研究日益增多(Zhang et al.，2011；
Huang et al.，2018；
Li et al.，2019；
Wang et al.，2019；
Jiang et al.，2020；
王文种等，2020；
Xu et al.，2022).同时，测高卫星搭载的高度计也逐渐从传统雷达高度计过渡到了合成孔径雷达(Synthetic Aperture Radar，SAR)高度计.Cryosat-2是第一颗搭载SAR高度计的卫星，而Sentinel-3则是第一颗以SAR模式覆盖全球的卫星(Jiang et al.，2020).
卫星测高技术最初是针对海洋问题设计的，将其应用于内陆水体不可避免地会遇到许多新问题.内陆湖库环境复杂，回波波形常被陆地信号污染(赵云等，2017；
刘琪等，2021)，限制了卫星测高的精度.波形重跟踪(也称波形重定)能有效地降低回波被污染带来的影响(Shu et al.，2020).经典的重跟踪算法主要有两类，建立模型对回波进行拟合的物理算法和基于统计规律的经验算法(Villadsen et al.，2016).国内外学者在此基础上又提出了许多改进的算法.Guo等(2006)提出的改进阈值重跟踪算法(Improved Threshold Retracker，ITR)、Jain 等(2015)提出的狭窄主波峰重跟踪算法(Narrow Primary Peak Retracker，NPPR)、Villadsen等(2016)提出的多回波波峰一致重跟踪算法(Multiple Waveform Persistent Peak retracker，MWaPP)提高了对于被污染回波的处理能力，得到了广泛的应用.Huang等(2018)提出的TIC(Threshold and Ice-1 Combined)算法提高了卫星测高提取高山区狭窄河流水位的能力.Shu等(2020)提出的双峰校正算法(Bimodal Correction Algorithm，BCA)在一定程度上解决了卫星测高对高纬度地区冰封湖泊水位明显低估的问题.此外，对回波进行波形分类(彭福凯和沈云中，2015；
Villadsen et al.，2016；
袁翠，2017；
廖静娟等，2020)和波形净化(黄征凯等，2021)可有效提高重跟踪算法的精度.
现有重跟踪算法中，ITR具有计算简便、精度高等优点(刘琪等，2021).对于多子波回波，该算法计算回波中每一子波对应的水位，其中只有一个是正确的水位，需要利用外部信息进行判断(Guo et al.，2006).对于湖泊而言，可以利用数字地形模型(Digital Terrain Model，DTM)、数字高程模型(Digital Elevation Model，DEM)、重力场模型计算的大地水准面起伏，湖泊水量模型计算的湖泊水位，上一周期卫星测高反演的水位等参考信息判断正确水位(高永刚等，2008；
Huang et al.，2018).然而，有些湖泊难以获取参考信息(Schwatke et al.，2015)，且当多子波回波中子波之间的距离较小时，无法通过外部参考准确判断正确水位(Huang et al.，2018).高贤文等(2021)以同一轨道单峰回波对应的水位为参考确定正确水位，该方法精度高，但并非所有轨道都有单峰回波.Huang等(2018)以当前周期标准差最小为准则确定正确水位，而对于周期内足迹较多的湖泊水库，该方法计算量巨大，因此，难以适用于湖泊水库.正确水位的选择决定了ITR的可靠性(Huang et al.，2018)，但目前确定正确水位的方法并不完善.为此，本文对ITR进行优化，提出了一种确定正确水位的新方法，该方法不依赖外部信息，通过寻找正确水位所在区间，自动确定正确水位.选择东平湖和峡山水库为试验区，评估和对比优化后的ITR、MWaPP和OCOG提取湖库水位的能力.
1.1 研究区概况
本文以东平湖和峡山水库为研究区，评估各波形重跟踪算法的性能.东平湖和峡山水库具有充足的实测日水位数据，可用于验证卫星测高提取水位的精度.
东平湖(图1a)位于山东省泰安市东平县西部(116°06′E—116°18′E，35°43′N—36°07′N).湖区属暖温带大陆性半湿润季风气候，年均气温13.6 ℃，年均降水640.5 mm，降水集中在夏秋两季.东平湖水面面积约120 km2，平均水深2.5 m，蓄水总量3亿m3，是山东省第二大湖泊.它与黄河、大汶河相连接，起到蓄纳洪水的作用，是南水北调东线工程的重要调蓄湖泊，也是山东省西水东送工程的水源地(韩非等，2020).
图1 东平湖(a)、峡山水库(b)地理位置及Sentinel-3卫星的地面轨迹Fig.1 Geographical location of Dongping Lake (a) and Xiashan Reservoir (b) and the ground track of Sentinel-3 satellite
峡山水库(图1b)位于山东省中东部(118°42′E—119°42′E，35°38′N—36°30′N)，隶属于潍坊市.库区气候与东平湖相似，年均气温12.3 ℃，年均降水627.3 mm，降水集中在夏季.峡山水库水面面积约97 km2，总库容14亿m3，是山东省第一大水库.水库除防洪功能外，还为昌邑、安丘、高密、寒亭等灌区供给农业用水，灌区总面积26.2万ha，为周边市区提供工业及生活用水，是一座功能齐全的大型水利枢纽(许传志，2015).
1.2 数据来源
(1)Sentinel-3卫星测高数据
Sentinel-3A/B卫星分别发射于2016年和2018年，轨道倾角98.56°，运行周期27天.卫星搭载的SRAL高度计，能发射Ku和C两种波段的微波，其中Ku波段用于测距，C波段用于确定电离层误差，拥有低分辨率(Low Resolution Mode，LRM)和合成孔径雷达(Synthetic Aperture Radar，SAR)两种模式(Jiang et al.，2020).本文采用SAR模式下采集的2016—2020年Level-2陆地非实时关键(Non-Time Critical，NTC)产品计算湖库水位.产品获取地址为https:∥scihub.copernicus.eu/dhus/#/home.产品沿轨分辨率约300 m，在东平湖和峡山水库的平均轨道长度分别为7.6 km、12.9 km，包含简化测量、标准测量和增强测量三个数据文件，本文使用的是增强测量数据.
(2)水体数据
本文所用水体数据来自欧盟联合研究中心于 2016 年发布的JRC(Joint Research Centre)全球水体数据集.该数据集中的月水体历史数据为1984年3月—2020年12月期间逐月全球水体数据，空间分辨率为30 m(Pekel et al.，2016).本文在Google Earth Engine云平台结合Landsat影像对月水体历史数据检查编辑，使用检查后的数据生成水体掩膜.
(3)实测数据
用于评估算法性能的实测数据为东平湖(2016年3月—2020年11月)和峡山水库(2018年12月—2020年11月)逐日实测水位数据.东平湖布设有遥测水位计和水尺(图1)，同时监测多个位置的水位，遥测水位计与邻近水尺监测结果互相核对，以保证数据的准确性.峡山水库以水尺为监测水位的主要手段.
2.1 水位提取的基本方法
当卫星经过湖库上方时，雷达高度计向水面发射微波脉冲，同时接收来自水面的回波，通过测量发射与接收的时间间隔来确定卫星到水面的距离，再结合卫星高度便可得到湖库水位(郭金运等，2013；
赵云等，2017).本文利用全球水体数据集中的月水体历史数据计算像元在研究期内被识别为水体的频率(水体频率)，生成东平湖、峡山水库永久性水体掩膜(水体频率为100%)以筛选湖库范围内的测高数据(Yang et al.，2020)，之后根据式(1)提取湖库水位：
Hwater=Halt-R-V-G,
(1)
式中，Hwater是湖库水位；
Halt是以WGS-84为参考椭球的卫星高度；
V是大气及地球物理改正，包括：电离层校正、湿对流层校正、干对流层校正、固体潮校正和极潮校正；
G是大地水准面EGM 2008到参考椭球的改正；
R是卫星到水面的距离，一般需经过波形重跟踪校正后才能得到精准的卫星-水面距离.
2.2 波形重跟踪
2.2.1 基本原理
利用回波前缘中点(图2)可以确定卫星-水面距离.理想情况下，前缘中点(Gr)应当位于预设跟踪门(Gnominal)上，但现实往往有偏移(郭金运等，2013).利用波形重跟踪可以计算出前缘中点与预设跟踪门之间的偏移量，从而对卫星-水面距离进行改正.如式(2)、(3)所示：
图2 Sentinel-3单峰回波波形及波形重跟踪原理示意图Fig.2 Sentinel-3 single-peak echo waveform and principle of waveform retracking
(2)
(3)
式中，R是波形重跟踪后的卫星-水面距离；
Runretracked是未重跟踪的距离；
VR是距离改正；
c为光速；
t为两相邻采样门间隔的时间.对于Sentinel-3测高数据，Gnominal=43(从0开始数起)，t=3.125 ns(Shu et al.，2020).确定Gr常用的方法包括重心偏移法、阈值法(Threshold Retracker，TR)、多回波波峰一致法和改进阈值法等重跟踪算法.
2.2.2 重心偏移重跟踪算法和阈值重跟踪算法
OCOG与TR确定前缘中点(Gr)的方法如下：
OCOG算法基于统计规律计算回波振幅、宽度以及重心的位置，进而确定出一个矩形，矩形与回波的第一个交点就是前缘中点(褚永海等，2005).Sentinel-3数据包含该算法改正后的卫星-水面距离，字段名称为range_ocog_20_ku，本文直接使用了该距离.
TR算法首先确定一个阈值T，然后在该阈值的几个邻近采样门之间进行线性内插，确定前缘中点(郭金运等，2010).本文参考已有研究(高贤文等，2021；
Villadsen et al.，2016)，设阈值为回波振幅的80%，回波振幅可由OCOG计算得到.
2.2.3 多回波波峰一致重跟踪算法和改进阈值重跟踪算法
对于内陆湖库，高度计接收到的回波通常是被污染的，即回波中包含从水面和从陆地返回的信号，形成具有多个子波的回波(图3).回波中每一子波都对应着一种反射面的高程，但只有一个子波对应水面高程，在这种情况下，OCOG和TR很难获取准确的水位.MWaPP和ITR通过寻找水面信号形成的子波，进而确定水位(Guo et al.，2006；
Villadsen et al.，2016).
图3 Sentinel-3多子波回波及子波提取示意图Fig.3 Sentinel-3 multi-subwaveform echo and subwaveform extraction
MWaPP实现步骤：
(1)将所有回波波形的横坐标——采样门，转换为高程，并将其插值到1 cm以下.
(2)对于转换后的回波，计算每一回波与前后两个回波的平均波形.
(3)寻找回波波形中第一个超过平均波形最高功率20%的点F，取回波波形中F点前后各三个采样门作为子波，认为该子波即水面信号形成的子波，最后利用TR算法确定该子波的前缘中点(Villadsen et al.，2016).
ITR实现步骤：
(1)利用回波差分的标准差提取回波中所有子波(图3).
(2)利用阈值法确定出每个子波的前缘中点(图3)，进而每个子波都可以确定出一个与之对应的水位，但只有一个水位是正确的(图4).
(3)利用外部信息，如模型值、已有测量值来确定正确水位(田山川等，2018).对于中小水体，确定正确水位较好的方法是以同一轨道的单峰回波(可以视为没有受污染)对应的水位为参考，与参考水位距离最近的即为正确水位(高贤文等，2021；
袁翠，2017)(图4a).位于峡山水库的所有轨道均有参考水位，而对于东平湖而言，43%的轨道没有参考水位，无法通过该方法来确定正确水位.
2.2.4 优化后的ITR
优化后的ITR根据湖库水面近似水平的特点，不依赖外部信息可以自动从多个子波对应的水位中确定正确水位.
正常情况下，因湖库水面近似水平，一条轨道上的正确水位应大致相同，集中分布于某个较小的高度区间内(田山川等，2018；
文京川等，2018).而湖岸由于地形起伏，各部反射率不同，且卫星不在其正上方，所以由湖岸反射信号计算的水位分布散乱.如图4所示，正确水位集中在一个小区间内，理想情况下，其他区间所包含水位的个数要远小于该区间，基于该特性，便可找到正确水位所在区间，进而筛选出正确水位.为了便于表述，本文约定，一条轨道上的所有回波按纬度从低到高分别为w1,w2,…,wn，若wi中有k个子波，则这些子波对应的水位分别为hi1,hi2,…,hik，记为Hi=[hi1,hi2,…,hik].此外，正确水位所在的区间记为正确区间.确定正确水位的步骤如下：
(1)设w1中有m个子波，分别以h11,h12, …,h1m(H1)为中心，生成m个宽度为d的区间.若H2中有水位不在已生成的区间之中，则分别以这些水位为中心生成新的宽度为d的区间.之后，同样以H3中不在已生成区间内的水位为中心，生成新区间，以此类推，直到Hn为止.图4b展示了生成的部分区间.
图4 回波子波对应水位及生成的区间与正确区间(a) 2019年3月16日，东平湖轨道上有参考水位；
(b) 2018年10月5日，东平湖轨道上没有参考水位.Fig.4 Corresponding water level of subwaveform and generated interval and correct interval(a) Dongping Lake on March 16, 2019, with reference water level on the track; (b) Dongping Lake on October 5, 2018, without reference water level on the track.
(2)计算各区间内水位的平均值，若Hi中有多个水位位于同一区间，则只保留最接近平均值的水位.
(3)计算各区间包含水位的个数，包含水位个数最多的区间记为q，若有区间与q相交，则将其与q合并为Q，Q为正确区间.若没有区间与q相交，则q为正确区间.筛选位于正确区间内的水位.若某回波wi中的所有子波对应的水位Hi均不在正确区间中，则该点处的正确水位为空值，这样就使得提取的水位均在正确区间内，从而避免产生异常水位.
该算法涉及区间宽度的选取.一般来说，同一轨道正确水位围绕水位真值上下波动，区间宽度可视为波动的合理范围.本文对单峰回波占比超过80%的轨道去除多子波回波后，利用TR算法反演这些轨道的水位并计算每一轨道最高与最低水位的差值，以差值的平均值作为区间宽度.通过计算，在东平湖和峡山水库，区间宽度分别为0.43 m、0.42 m.该算法可根据不同输入数据生成不同的区间，同时可利用步骤(3)的区间合并自动调整正确区间的宽度，因此具有一定程度的自适应性.
2.3 异常水位剔除与评价指标
若提取结果中存在异常值，则利用稳健统计分析将其剔除.对于同一轨道上的水位，根据式(4)、(5)计算该轨道所有水位的统计分数(ZMAD)，剔除统计分数大于3的水位(Shu et al.，2020).
MAD=1.4826×med{|Hall-med{Hall}|},
(4)
(5)
式中：Hall为一条轨道上的所有水位；
med{Hall}为水位中值；
MAD为中值绝对偏差.剔除异常水位后取中值，作为该湖泊的日水位.由于水文站使用的高程基准未知，本文在反演的水位序列和实测水位序列(东平湖实测水位为多位置水位均值)中减去各自的均值得到水位距平时间序列.计算反演水位距平和对应实测水位距平的决定系数(R2)与均方根误差(Root Mean Squared Error，RMSE)，评估不同重跟踪算法的精度.
3.1 精度验证
图5为不同算法提取的水位距平与实测水位距平的结果对比.图中,P<0.01表示在0.01显著性水平上显著．从图中可以看出，三种重跟踪算法中，优化后的ITR精度最高，在研究区内平均RMSE为6.5 cm，与MWaPP和OCOG相比，精度分别提高了2.7 cm和24.2 cm.三种算法在峡山水库的RMSE均达到厘米级，精度较高.优化后的ITR和MWaPP在东平湖表现较好，但OCOG精度远低于前两者，RMSE为53.7 cm，且有多个严重偏离实测值的点(图5).主要原因是东平湖回波被污染的情况更为严重，导致OCOG提取的水位异常值占比过高.如图6所示，即使多次利用稳健统计分析，仍有较多的异常值残留.MWaPP提取的水位异常值较少，容易剔除，优化后的ITR则很少产生异常值.由于东平湖回波质量差，且单条轨道上足迹较少(平均23个，峡山水库平均39个)，三种算法提取水位的精度在东平湖比在峡山水库有所降低，而优化后的ITR精度下降程度最小，RMSE在东平湖比在峡山水库升高了2.0倍，MWaPP与OCOG则分别升高了2.9倍、7.0倍.
图5 反演的水位距平与实测水位距平对比Fig.5 Comparison of retrieved water level anomaly and measured water level anomaly
图6 东平湖异常值剔除前后的所有水位Fig.6 All water levels before and after removing outliers in Dongping Lake
3.2 不同重跟踪算法的成功率
为更好地评价不同算法的重跟踪效果，本文统计了各周期重跟踪成功波形个数与回波个数之比，作为各周期的重跟踪成功率(高贤文等，2021).将明显偏离正常区间的重跟踪结果认定为重跟踪失败，主要利用稳健统计分析进行识别.当个别周期某重跟踪算法提取的结果异常值占比过高，稳健统计分析失效时，参考其他算法提取的结果，人工计算重跟踪成功率.图7为不同重跟踪算法成功率的频率分布直方图.图中，μ表示期望，σ表示标准差.优化后的ITR、MWaPP、OCOG在研究区内平均重跟踪成功率分别为0.92、0.83、0.82.由于回波质量较差，使用三种算法在东平湖进行重跟踪时，分别有2%、14%和18%的周期成功率在0.60以下.三种算法在峡山水库进行重跟踪时，每一周期的成功率均在0.60以上.受水库形状影响，峡山水库范围内的足迹可分为南、北两部分(图1)，质量较差的回波主要集中在南半部分.三种算法在北半部分的平均重跟踪成功率分别为0.96、0.91、0.89，而在南半部分分别为0.83、0.45、0.70.MWaPP在南半部分效果较差，导致该算法在峡山水库平均成功率最低.优化后的ITR平均成功率最高，且在东平湖和峡山水库均相差不大，说明该算法较另外两种算法更为稳定.
图7 不同重跟踪算法在东平湖(a—c)、峡山水库(d—f)的成功率频率分布直方图Fig.7 Histograms of frequency distribution of success rate of different retracking algorithms in Dongping Lake (a—c) and Xiashan Reservoir (d—f)
3.3 悬挂效应对不同重跟踪算法的影响
本文提取的大部分水位表现正常，而利用MWaPP与OCOG提取的2018年6月19日东平湖水位沿卫星轨迹呈V型分布(图8)，该现象是由悬挂效应引起的(Frappart et al.，2006；
Da Silva et al.，2010；
田山川等，2018).MWaPP与OCOG基于回波功率识别水面信号，有跟踪强信号的倾向(Jiang et al.，2020)，因此当回波中湖岸信号过强，而湖面信号较弱时，易将湖面信号误判为噪声，从而跟踪湖岸信号，这时计算的卫星-湖面距离实际上是卫星与湖岸的斜距.卫星距离湖岸越远，斜距越大，提取的水位越低，最终整个沿轨水位呈现出先下降后上升的V型.优化后的ITR提取的结果比较正常(图8).优化后的ITR、MWaPP、OCOG提取的该日湖面点水位的标准差分别为9.0 cm、109.4 cm、114.2 cm.
为更好地分析悬挂效应对不同算法的影响，本文将提取结果分为三部分(区域a—c).如图8所示，优化后的ITR在区域c表现最好，而在区域a、b无明显优势.当卫星位于区域a(靠近湖岸)时，湖岸信号与水面信号基本重叠，卫星-湖岸的斜距与卫星-水面距离相差不大，利用湖岸信号也能近似得到湖库水位，三种算法的结果基本一致.当卫星位于区域c(远离湖岸)时，水面信号与湖岸信号完全分离，MWaPP和OCOG跟踪湖岸信号得到的卫星-湖岸斜距远大于卫星与水面的距离，因此计算的水位偏低.优化后的ITR根据同一轨道水位大致相同的特点，识别出水面信号，得到了较为理想的结果.当卫星位于区域b时，两种信号处于半分离状态，此时MWaPP的效果优于OCOG，但仍不理想.优化后的ITR只提取位于正确区间内的水位，因此在这些位置得到的水位为空值，而不会产生异常值.
图8 2018年6月19日东平湖沿卫星轨迹水位及不同区域的典型波形Fig.8 Water level along satellite track in Dongping Lake on June 19, 2018, and typical waveforms of different regions
卫星测高技术是监测湖库水位的重要手段，波形重跟踪可以有效提高其精度.本文提出了一种针对ITR重跟踪算法的优化方案，同时，利用Sentinel-3A/B合成孔径雷达高度计观测数据，基于优化后的ITR、MWaPP和OCOG重跟踪算法提取了东平湖2016年3月—2020年11月、峡山水库2018年12月—2020年11月的水位序列.主要结论如下：
(1)三种重跟踪算法中，优化后的ITR表现最佳，在研究区内平均RMSE为6.5 cm，MWaPP精度稍差，OCOG精度最差，且在东平湖由于提取结果异常值占比过高，因而精度远低于另外两种算法，MWaPP提取结果中的异常值较容易剔除，优化后的ITR提取结果几乎不含异常值.
(2)优化后的ITR、MWaPP、OCOG在研究区进行波形重跟踪时平均成功率分别为0.92、0.83、0.82.优化后的ITR平均成功率最高，在东平湖和峡山水库成功率均相差不大，表明该算法较为稳定.
(3)受悬挂效应影响，MWaPP与OCOG提取的2018年6月19日东平湖水位沿卫星轨迹呈V型分布，优化后的ITR提取的该日东平湖水位比较正常.优化后的ITR、MWaPP、OCOG提取的该日湖面点水位的标准差分别为9.0 cm、109.4 cm、114.2 cm.
正确水位的选择决定了ITR的可靠性.但当同一轨道没有单峰回波时，现有方法难以准确确定正确水位，本文提出的方法有效地缓解了这一问题.目前没有一种通用的确定正确水位的方法，特别是对于环境复杂地区，需要综合利用多种方法才能确定正确水位.这些地区更易产生同一轨道没有单峰回波的情况，因此将本文方法与已有方法结合，能有效地提高ITR算法在环境复杂地区的反演精度.本文仅使用了Sentinel-3数据监测研究区内水位变化，未来将基于多种测高数据，对比不同重跟踪算法在各类型水体上的表现，同时从气候变化、人类活动、水源补给等多角度分析水位变化原因.
猜你喜欢东平湖湖库子波一类非线性动力系统的孤立子波解数学物理学报(2021年6期)2021-12-21——东平湖增殖放流活动实施">修复渔业资源，改善水域环境——东平湖增殖放流活动实施科学养鱼(2021年5期)2021-11-30东平湖水环境质量评价及趋势研究皮革制作与环保科技(2021年14期)2021-11-12湖库富营养化形成原因和处理策略山西农经(2020年4期)2020-04-07再去看看东平湖（外二首）核桃源(2019年3期)2019-11-14卫星测高数据筛选方法研究南水北调与水利科技(2018年3期)2018-11-13秀山县湖库水质特征分析及富营养化评价绿色科技(2017年14期)2017-08-22基于Landsat数据的近30年东平湖湿地植被覆盖演变研究山东林业科技(2017年1期)2017-06-29地震反演子波选择策略研究采矿与岩层控制工程学报(2015年3期)2015-12-16加强生态清洁小流域建设 推进湖库型水源地水土保持工作中国水利(2015年6期)2015-02-28