卫星测高具有覆盖范围广、观测精度高等优势，在冰冻圈研究中发挥着重要作用。随着测高理论的不断完善，以及平台、控制、传感器等技术的进步，新的观测模式(如CryoSat-2的SAR、SARIn模式)逐渐被引入，大大拓展了雷达测高卫星的应用领域^[1-5]。多位学者利用CryoSat-2开展内陆湖泊变化等研究^[4-6]，但CryoSat-2数据产品仍处于不断更新完善中，特别是在联合其他测高卫星(如ICESat)研究长周期地球物理现象方面仍有待加强。在联合利用ICESat(2003~2009)和CryoSat-2(2010~2017)的SARIn模式二级产品(基线C)提取青藏高原及天山地区湖泊高程时间序列时发现，2套卫星数据存在几十米的偏差。为此，本文从CryoSat-2一级产品入手，采用重跟踪方法重新获取湖面高程，分析导致偏差的原因，并结合实测数据对重跟踪结果及二级产品作对比分析。

1 CryoSat-2卫星及测高产品

CryoSat-2的主要科学目标是观测地球上永久性冰体(包括南北极冰盖、山地冰川和多年海冰)的厚度变化，并预测气候和海平面变化^[7]。其主要载荷是合成孔径干涉雷达测高仪(SAR interferometric radar altimeter, SIRAL)，它具有3种测量模式：低分辨率模式(low resolution mode，LRM)、合成孔径雷达测量模式(synthetic aperture radar, SAR)和合成孔径雷达干涉测量模式(SAR interferometric, SARIn)，3种模式的关键设计参数见表 1。LRM模式主要用于相对平坦的两极冰盖内陆区域和海洋观测；SAR模式能提高沿轨空间分辨率，主要用于海冰测量；SARIn模式通过增加一副接收天线，对接收信号进行相位差分处理，提高交轨方向的分辨率，适合地势较为复杂的冰盖边缘和山地冰川^[7]。

图 1给出SARIn模式的相位差分示意图。在倾斜地面上，卫星测量的点并不是星下点，但SARIn模式可以精确测量其点位^[8]。

目前ESA官网上发布的CryoSat-2数据产品版本为基线C，包括L1b数据(一级产品)和L2I数据(二级产品)。其中，L1b数据主要包括时间、未经改正的经纬度、回波强度和地球物理改正等参数，L2I数据提供时间、改正后的经纬度以及高程等参数。本文要用到的数据是SARIn模式下的L1b数据(SIR_SIN_L1)和L2I数据(SIR_SIN_L2I)。为了获得更长时间的湖泊高程观测，我们联合ICESat陆地测高产品(GLAH14)，其包含时间、经纬度、高程及大地水准面差异等参数。

2 L2I数据偏差

在研究纳木错水位变化时，由于ICESat与CryoSat-2两套卫星数据采用的参考椭球及高程基准不同，为进行联合分析，本文将2套数据统一到WGS84坐标系统下。图 2给出ICESat与CryoSat-2在纳木错的几条轨迹。

结合纳木错边界数据提取出落在边界内的光斑。在提取湖泊高程时，考虑到测高数据在湖边易受岸上地物反射信号的干扰，将湖泊边界向中心缩进200 m，另采取以下3个步骤去除粗差、减少误差的影响：1)剔除所有数据中超过中值一定范围(设定为100 m)的数据；2)对某一时间段(如同一轨迹)的高程数据，剔除超过中值2 m的数据；3)在进行以上两步处理后，对某一时间段(如同一轨迹)的高程数据用3倍中误差迭代剔除误差。经过上述处理得到的纳木错高程时间序列中发现，2套数据存在几十米的偏差(图 3(a))。

为了调查两套数据是否在所有湖泊都存在一个几十米的偏差，本文提取青藏高原及天山地区所有湖泊的高程，从中挑选出满足要求的数据来查验。挑选标准是2套数据同时满足以下3个条件：1)湖泊中既有ICESat轨迹又有CryoSat-2轨迹；2)ICESat在2008年以后、CryoSat-2在2011年以前要有光斑点落在湖泊内；3)落在湖泊中的光斑点要足够多，在时间序列图中至少有6个点，以增加统计的可靠性。最终选取158个湖泊并得到ICESat与CryoSat-2的湖水高程偏差。经统计，所有湖泊都存在几十米的偏差，其中153个湖泊(约占96.8%)偏差在57~62 m之间。这表明，在内陆水域普遍存在一个系统性偏差。各湖泊高程偏差值的散点分布见图 3(b)。

多位学者对ICESat数据的精度及可靠性进行验证^[9-10]。以往研究使用的ICESat数据多为旧版本(如V31、V33)，本文数据为最新版本(V34)，数据质量有所提高。王莉等^[11]对ICESat的最新数据(V34)与V33数据进行高程差统计表明，均值和标准差均在精度范围内。综合得出，ICESat V34版本数据在湖泊高程监测上是可靠的, 几十米的系统性偏差很可能是由CryoSat-2数据产品引起的。

3 处理方法 3.1 波形重跟踪原理

测高卫星通过测量脉冲在卫星和地表之间的往返时间来确定卫星与地面间的距离^[12]，波形前缘半功率点(前缘中点)对应地面点。高度计在设计时往往设定某个波形采样点为波形预设中点^[12]。卫星在运行过程中，记录的是脉冲从发射到预设中点的往返时间，波形重跟踪的目的是找到波形的前缘中点并与波形的预设中点比较，得到重跟踪改正值。

利用重跟踪获得波形前缘中点采样值，可以计算出重跟踪改正及改正后的星地距离。根据卫星测高基本原理，结合卫星到参考椭球面的高度进行一些地球物理改正(如大气、潮汐)，就可以得到地表高程。

3.2 波形重跟踪方法

准确地确定波形前缘中点是精确获取地面高程的关键。由于地物的多样性，回波的波形特征也不尽相同，目前还没有一种重跟踪方法适用于所有的地物类型。许多学者针对不同的地形及波形特征提出不同的波形重跟踪方法^[13-17]，每一种算法在其适用的地物类型和波形特征范围内对测高精度都有一定的改善。

本文采用经典的Threshold算法进行波形重跟踪，主要采用插值法, 以计算得到的振幅(A)为基础，选定一个阈值因子(T_h)，如25%、50%、75%、85%，振幅与阈值因子的乘积为阈值(T_l)，该阈值即认为是波形前缘中点对应的回波功率值，通过相邻采样点内插出阈值处的采样值(G_r)。

Threshold计算公式^{[12, 14]}如下：

$ G_{r}=K-1+\frac{T_{l}-P_{k-1}}{P_{k}-P_{k-1}} $

(1)

其中，

$ \begin{array}{c}{T_{l}=\left(A-P_{N}\right) \times T_{h}+P_{N}} \\ {A=\max \left(P_{i}\right)} \\ {P_{N}=\sum\limits_{i=1}^{5} P_{i} / 5}\end{array} $

式中，P_i和P_k分别为第i和第k个采样点的功率值，P_N为前5个采样点功率值的平均，K为振幅大于T_l的点所对应的最近采样点。

4 重跟踪结果及分析 4.1 重跟踪结果

在纳木错选取2条(ICESat和CryoSat-2各1条)非常靠近的轨迹，如图 4(a)所示。利用Threshold方法，分别取阈值因子为25%、50%、75%和85%(记为Th25、Th50、Th75和Th85)对CryoSat-2的SARIn模式L1b数据进行重跟踪来获取轨迹上各光斑点的高程结果。图 4(b)为L1b数据重跟踪后得到的高程及ICESat得到的高程。

由图 4(b)可知，利用Threshold算法进行波形重跟踪后，ICESat与CryoSat-2之间几十米的偏差被消除。重跟踪后的CryoSat-2轨迹高程没有与ICESat完全重合，主要原因如下：1)两条轨迹的时间不同，ICESat观测时间是2008年，CryoSat-2为2012年，在此期间湖水水位可能发生年际变化；2)季节性差异，ICESat轨迹时间为10月，CryoSat-2则在3月，季节不同导致的冰雪消融、降雨等因素会影响湖泊的水位；3)重跟踪过程中，阈值的选取也影响高程结果；4)CryoSat-2与ICESat卫星间的偏差(如雷达穿透性、仪器偏差等)。对比结果发现，波形重跟踪能够去除几十米偏差，说明系统性偏差主要是由CryoSat-2二级产品引起的。

4.2 原因分析

针对基线B产品存在的问题，如与姿态有关的异常、测距和时间偏差等^[18]，ESA进行部分修正，发布了新的产品(基线C)。在新的产品中，L1b波形采样数由原来的512变为1 024，因此窗口延迟(window delay)^[7]参数也应更新，使其在L1b波形中的参考位置为更新后波形的中点，即参考位置由256变为512。但由于参数的更新会影响SARIn二级产品其他参数的质量^[19]，因此在SARIn模式由一级产品生成二级产品时并未使用更新后的窗口延迟，由此产生偏差。由于波形采样间距约为0.234 2 m(表 1)，使用未更新的窗口延迟参数产生的偏差约为59.955 2 m(256个采样间距)。而在基线C版本的L1b数据中记录的是更新后的窗口延迟参数，因此可以对L1b数据进行波形重跟踪来消除偏差，或者直接在二级产品上加一个偏差改正^[19]。

4.3 跟踪方法对比分析

为了评估哪个阈值更加合适，将不同阈值得到的纳木错高程与实测数据比较。评估参数主要是2套数据的差值平均值、标准偏差σ和相关系数r，标准偏差和相关系数的表达式为：

$ \sigma=\sqrt{\frac{\sum\limits_{i=1}^{n}\left(d_{i}-\overline{d}\right)^{2}}{n-1}} $

(2)

$ r=\frac{\sum\limits_{i=1}^{n}\left(X_{i}-\overline{X}\right)\left(Y_{i}-\overline{Y}\right)}{\sqrt{\sum\limits_{i=1}^{n}\left(X_{i}-\overline{X}\right)^{2} \sum\limits_{i=1}^{n}\left(Y_{i}-\overline{Y}\right)^{2}}} $

(3)

式中，d_i为CryoSat-2数据与实测数据的差值，d为差值平均值，X_i、Y_i为2套数据的高程值，X、Y为平均值。标准偏差σ反映偏离平均值的程度，值越小代表数据精度越高；相关系数反映2套数据的相关性，绝对值越大表示相关性越强。

本文使用了纳木错水位实测数据，相关数据由国家科技基础条件平台建设项目“青藏高原科学数据中心”提供(www.tpedatabase.cn)。由于实测数据是水位相对变化量，而ICESat的高程基准为WGS84，根据水位观测点与ICESat公共时间段(2007~2009年)的水位数据计算二者差值的平均值，然后将实测水位归算到ICESat高程基准。在时间序列图中，水位观测点采用的数据只是一点的水位，而卫星测高采用的是一条轨迹或几条轨迹上所有点高程的平均值，且纳木错东西向跨度较大(最宽处约73 km)，2种水位结果可能存在较大的差异。因此，在测高数据与实测数据联合使用时，只选取离观测站较近的轨迹(经度大于90.58°)。同理，通过CryoSat-2二级产品与实测数据间的公共时间段(2010~2016年)计算二者之间的差值平均值(-59.553 m)，将二级产品归算到同一基准下。图 5为将ICESat、实测数据以及CryoSat-2二级产品归算到同一基准下的时间序列。ICESat与实测数据的相关系数高达0.983，二级产品与实测数据的相关系数达0.672，ICESat和CryoSat-2二级产品与实测数据吻合较好。

利用波形重跟踪得到Th25、Th50、Th75和Th85的测高结果，CryoSat-2二级产品与实测数据在公共时间段内所有的重叠点共41对，选用差值平均值、标准偏差及相关系数等参数对重跟踪方法质量进行评定。表 2为几种方法的比较，其中差值为测高数据与实测数据的高程差，二级产品高程加上了因窗口延迟产生的偏差改正(59.955 2 m)。

由表 2可知，使用重跟踪方法后，二级产品与Th25、Th50的标准偏差及相关系数都优于Th75和Th85；二级产品与Th50的标准偏差基本相同，优于Th25；Th25和Th50的相关系数比二级产品好；综合来看Th50更有优势。另外，实测数据与CryoSat-2之间存在一个小偏差，主要由几部分组成：CryoSat-2与ICESat卫星间的偏差^[20-21](如雷达穿透性^[22-23]、不同仪器间偏差^[23]、接收器灵敏度等)及不同重跟踪方法产生的偏差等。此外，实测数据是经过改正归算到ICESat高程基准上的，这一过程也可能存在误差。本文只对纳木错实测数据与CryoSat-2数据比较，资料有限，因此还不能确定CryoSat-2与ICESat卫星间偏差的具体值。

图 6为Th50与CryoSat-2二级产品的比较，都归算到了实测数据上(二级数据加了59.553 m的偏差改正，Th50减去0.651 m)。可以看出，二者整体趋势基本一致，二级产品跳动性较大的点较多，Th50时间序列相对稳定。相对于二级产品，重跟踪结果质量有一定的改善。

因此，如果考虑用Threshold法进行重跟踪，50%的阈值因子比较合适。如果直接在二级产品上加一个偏差改正，由于偏差改正值的大小受实测数据精度、重跟踪方法等因素的影响，本文仅根据纳木错实测数据得到2套数据间的偏差改正值为59.553 m，因此，对于其他的湖泊并不能给出统一的偏差改正值。

5 结语

本文在联合CryoSat-2和ICESat提取青藏高原及天山地区湖泊高程时发现2套数据存在几十米的偏差，证实偏差主要由CryoSat-2产品质量引起。利用重跟踪方法对SARIn模式L1b数据(基线C)进行处理后，能够很好地消除偏差，并利用纳木错实测数据对几种阈值因子的Threshold法及二级产品进行对比分析。通过调查发现，出现偏差是由于现有的二级产品在生成过程中使用了错误的窗口延迟参数，因此CryoSat-2的SARIn模式二级产品(L2I)不能直接用于多源数据的联合分析，需要对一级产品重跟踪或者直接在二级产品上加一个偏差改正。

致谢: CryoSat-2数据来源于欧洲航天局(ESA)，ICESat数据来源于美国国家冰雪数据中心(NSIDC)，水位实测数据由中国科学院青藏高原研究所纳木错站提供，在此表示衷心感谢！