适用于Ku波段雷达高度计海冰干舷高度反演的积雪校正方法

张 晰，赵全芳,孟俊敏，刘眉洁

（1.山东科技大学 测绘科学与工程学院，山东 青岛 266590；2.自然资源部第一海洋研究所，山东 青岛 266061;3.青岛大学 物理学院，山东 青岛 266071）

海冰直接影响着全球的气候变化，并通过与大气的相互作用反馈于全球环境系统。海冰厚度是重要的气候环境变化表现因子。准确估算海冰厚度，对于气候变化研究、极区航行保障具有重要意义。利用遥感手段准确估算大范围海冰厚度始终是当前的国际难题。相比于人工实测、电磁感应和仰视声呐等手段，近年来发展的雷达高度计测高技术使得获取连续大范围的海冰厚度成为可能[1]。到目前为止，ERS-1/2、ENVISAT、HY-2、AltiKa、CryoSat-2和Sentinel-3A等雷达高度计相继发射，为海冰厚度估算提供了多源的卫星观测手段。在众多雷达高度计中，工作在Ku波段的CryoSat-2（CS-2）和Sentinel-3A（S3）合成孔径雷达高度计相较于传统雷达高度计，具有更高的空间分辨率（沿轨约0.3 km，CS-2/S3交轨约1.5/1.64 km），能够提供较高分辨率的海冰厚度信息，是当前最先进的卫星雷达高度计。

雷达高度计并不直接计算海冰厚度，是通过计算海冰的冰面出水高度（海冰干舷高度），并结合浮体定律实现海冰厚度估算[2-3]。因此，准确反演海冰干舷高度对于估算海冰厚度至关重要。由于卫星高度计信号回波的实际时间跟踪点与预设时间跟踪点存在位置上的偏差，需要对回波波形进行重跟踪处理，以计算实际时间跟踪点，进而校正得到卫星高度计质心到地面点的真实距离。Ricker等[4]使用波形重跟踪方法指出重跟踪阈值的选择对海冰干舷和厚度的估计值有显著影响。在极区，海冰表面普遍覆盖着积雪，当雷达高度计探测海冰表面高度时，Ku波段电磁波须穿过覆盖于海冰之上的积雪层，因此在计算海冰干舷高度时，需考虑电磁波在积雪中穿透速度较慢的距离校正问题。在实际情况下，积雪层受海冰析盐过程的影响，常会在积雪的上下层形成强盐度梯度变化[5-7]；另外，积雪层内粒径也会与电磁波发生体散射[8]。这些因素使得微波雷达高度计的电磁波很难完全穿透积雪层直达海冰表面。对于短波长的Ku波段雷达，该现象更为明显，这将导致雷达信号主散射面向上偏离雪-冰交界面，从而影响海冰干舷高度的估计。但大多数研究者在使用雷达高度计反演海冰干舷高度时，常假设Ku波段电磁波可以完全穿透积雪层，直接将雪-冰交界面当作雷达信号的主散射面[2-4,9-12]，利用海冰上的整体积雪深度进行距离校正。Armitage和Ridout[13]分析了CS-2雷达高度计电磁波在积雪中的雷达穿透系数，发现Ku波段电磁波在垂直入射时无法完全穿透积雪层，并指出若不考虑积雪层对电磁波的影响，海冰干舷高度将被高估，从而大大影响海冰厚度的估算。Ricker等[14]发现，未知雷达信号穿透到积雪层中的不确定性可能导致大约0.06～0.12 m的海冰干舷高度偏差。

因此，选择合适的冰间水道与海冰的重跟踪阈值以及确定雷达信号在积雪中的穿透深度，对精确地反演海冰干舷高度起决定性作用。基于此，本文分别针对CS-2与S3两种高度计卫星开展了确定最优波形重跟踪阈值的分析，并详细分析了积雪层对Ku波段CS-2和S3合成孔径雷达高度计主散射面位置的影响，计算了Ku波段电磁波的积雪穿透系数，针对海冰干舷高度精确反演的需求，改进了积雪校正算法。

1 数据

1.1 CryoSat-2雷达高度计数据

CryoSat-2是欧洲空间局发射的合成孔径雷达高度计，其工作在Ku波段，中心频率为13.575 GHz，带宽约为320 MHz，空间覆盖范围达到南北纬88°，运行重复子周期为30 d，其沿轨道分辨率可达0.3 km，交轨道分辨率优于1.5 km，相比于足迹为10 km的传统雷达高度计[15]，其测量精度显著提升。该卫星有LRM（低分辨率）、SAR和SARin 3种工作模式，本文采用的是SAR模式L1b级数据。在该模式下，CS-2的每一个波形采样窗口大小为60 m，包含有256个距离门[16]。为了进行波形重跟踪校正，得到海冰/水的表面高程，还需对L1b数据进行波形噪声去除和地球物理校正（地球物理校正包含：干/湿对流层、电离层、逆气压、海洋潮汐、长周期平衡潮、海洋负荷潮、固体潮和地心极潮校正等处理）。本文使用的是北极区域、时间为2017年3月-4月和2018年4月的雷达数据。

1.2 Sentinel-3A雷达高度计数据

Sentinel-3A也是欧空局发射的合成孔径雷达高度计，它是第一个100%以SAR模式覆盖海洋、冰区和内陆水域的雷达高度计。其工作的Ku波段中心频率与CS-2相同，带宽约为350 MHz，空间覆盖范围为南北纬81.35°，运行重复子周期为4 d，沿轨道分辨率约为0.3 km，交轨道分辨率大约为1.64 km[17]。本文采用的是SAR模式的L2 WAT Enhanced（Water Enhanced）数据。在该模式下波形采样窗口大小为60 m，包含有128个距离门。与CS-2相同，使用的也是北极区域2017年3月-4月和2018年4月的数据，并对这些数据进行了与CS-2相同的地球物理校正处理。

1.3 Operation IceBridge数据

本文使用美国国家航空航天局（NASA）发布的OIB机载数据作为海冰干舷高度的验证数据，其每年均会在北极、格陵兰岛和南极上空开展机载探测飞行。该数据由机载地形扫描测绘仪（Airborne Topographic Mapper,ATM）、数字测绘相机（Digital Mapping System Camera,DMS） 和雪雷达（Snow Radar）3种装备采集得到。ATM用于测量表面高度，空间分辨率为0.4 m，垂直分辨率为0.03 m；DMS用于识别海冰、海水和冰间水道等，空间分辨率为0.1 m；雪雷达用于测量积雪厚度，空间分辨率为40 m，垂直分辨率为0.06 m[18]。综合上述3种传感器，机载OIB可提供海冰的激光干舷高度和积雪厚度等数据。图1展示了2017-2018年春季机载OIB北极飞行路线及激光干舷高度分布图。

对于NASA发布的机载OIB数据，其海冰激光干舷高度、积雪厚度和产品空间分辨率都统一为40 m[18]。在此需要指出的是：机载OIB数据携带的测高仪器为激光雷达，因此，得到的海冰激光干舷高度并非海冰干舷高度，而是海冰干舷高度与积雪厚度之和。为得到海冰干舷高度，需将OIB得到的激光干舷高度减去雪雷达得到的积雪厚度，才能得到真实的海冰干舷高度。根据文献记载，机载OIB的海冰激光干舷高度探测精度为0.015±0.06 m[19]，积雪厚度的探测精度为0.01±0.05 m[20]。本文所用的OIB数据均与CS-2和S3雷达高度计数据时间、地点范围相同（北极，2017年3月-4月和2018年4月）。由于CS-2、S3雷达高度计与OIB数据的空间分辨率不同，为进行逐一比较，将雷达高度计数据与OIB数据在25 km的网格中取平均，实现数据一一对应。

图1 春季机载OIB北极飞行路线及激光干舷高度分布图

2 Ku波段雷达高度计波形重跟踪阈值确定

计算海冰干舷高度可通过反演海冰面到雷达质心的高程和冰间水道面到雷达质心的高程，并计算二者的差值得到。因此，首先需要区分雷达高度计的回波波形，以识别海冰、开阔水域和冰间水道等不同地物类型，然后再反演各种地物类型表面的高程。

对于地物类型识别，主要是利用Ocean and Sea Ice SAF（OSI-SAF）发布的海冰密集度数据以及CS-2或S3的波形区分海冰、开阔水和冰间水道3种类型。常用后向散射系数（Radar Backscatter Coefficient,Sigma0）、波形前缘宽度（Leading Edge Width,LEW）和脉冲峰值（Pulse Peakiness,PP）3种回波波形特征进行类型识别[21]。Sigma0为雷达高度计接收到的地物的表面后向反射系数；LEW是最大峰值功率5%和95%点位之间的波形前缘宽度；PP是雷达波形最大峰值功率与同一采样波形里的所有波形总功率的比率[4]，计算方法如下：

式中：NWF代表一个波形内的距离门数；WF为雷达波形在第i个距离门处的功率。通常冰间水道的表面较平缓，雷达信号多为镜面反射，PP值较高；对于开阔水域和海冰则以发生漫反射为主，PP值相对较低，LEW较宽（开阔水域的LEW相比海冰更大）。常用的分类参数值设置[22-23]见表1。

表1 CryoSat-2和Sentinel-3A的雷达波形分类参数设置

雷达高度计是通过计算信号发射与返回的时间差和光速来计算雷达到星下地物点之间的距离，因此常需要预先设置信号发射与返回时的时间跟踪点。但在实际情况中，雷达接收信号的实际时间跟踪点与预设时间跟踪点之间存在偏差，需要根据实际情况进行波形重跟踪校正[24]。本质上，波形重跟踪是将发生主散射的波形位置定为重跟踪点，然后计算重跟踪点和预设跟踪点之间的偏移，以校正雷达高度计质心到星下地物点之间的真实距离[25]。

TFMRA方法（Threshold First Maximum Retracker Algorithm）[26]是海冰干舷高度反演中常用的波形重跟踪方法。该方法针对冰间水道和海冰等不同的地物类型，经验性地设置重跟踪点阈值位置。表2总结了当前最为主要的一些重跟踪阈值组合。为了更全面地比较不同重跟踪阈值组合对干舷反演的影响，并得到最优阈值组合，本文在表2的3种阈值组合的基础上重新构建了8种阈值组合，分别对应为(40%，40%)，(50%，40%)，(50%，50%)，(60%，40%)，(60%，50%)，(70%，40%)，(70%，50%)，(70%，60%)（括号内第一个阈值对应于冰间水道，第二个阈值对应于海冰）。以S3数据为例，这些阈值在雷达回波波形中的位置示意图见图2。同时，需要指出的是，冰间水道的阈值通常大于海冰的阈值，这是因为对于雷达高度计冰间水道的后向散射常高于海冰[27]。

表2 文献中常用的重跟踪阈值组合

图2 Sentinel-3A雷达高度计海冰回波波形中第30～60距离门处归一化回波功率阈值位置示意图

在完成重跟踪阈值设置后，利用波形重跟踪校正后的海冰表面高程减去相应的平均海表面高度与海表面高程异常值，得到海冰雷达干舷高度FR（Radar Freeboard）。由于电磁波在穿过积雪层时的速度与真空中不同，所以还需要利用积雪厚度数据对雷达海冰干舷高度FR进行距离校正，以得到真正的海冰干舷高度F。因此，当雷达信号完全穿透积雪层时，进行积雪校正的海冰干舷高度计算公式[4]如下：

式中：cs为积雪层中雷达信号的传播速度；c为雷达信号在空气中的传播速度；hs(1-cs/c)≈0.22hs为积雪校正值。需要注意的是，公式（2）的使用条件为：电磁波能够完全穿透雪层直达雪-冰交界面。若电磁波无法穿透雪层，应用上式计算海冰干舷高度将会不可避免地引入误差。

图3 2017年3月、4月与2018年4月北极地区不同阈值组合下海冰干舷高度与OIB海冰干舷高度散点图

表3 CryoSat-2卫星数据不同阈值组合计算的海冰干舷高度与OIB海冰干舷高度对比

将2017年3月、4月与2018年4月CS-2与S3在8种阈值组合方案下得到的海冰干舷高度分别与同期的OIB海冰干舷实测数据进行对比，从而确定两种Ku波段合成孔径雷达高度计针对海冰干舷高度反演的最优波形重跟踪阈值组合。为保证精度，在本文的处理中，积雪厚度数据采用的机载OIB数据提供的积雪厚度测量值。图3和表3给出了不同阈值组合得到的海冰干舷高度与机载OIB得到的海冰干舷高度对比结果。

从表3中可以看出，对于CS-2卫星数据，最优的阈值组合为70%（冰间水道）、60%（海冰）。在该组合下得到的海冰干舷高度与机载OIB数据的海冰干舷高度非常接近。观测表3可知，其平均绝对值差值和均方根误差均为最小值，分别为0.060 7 m和0.077 6 m。

对于S3卫星数据，由表4可知，最优的阈值组合为50%（冰间水道）、50%（海冰）。在该组合下的平均绝对值差值和均方根误差分别为0.075 0 m和0.098 2 m。

同时，由表3、表4中的数据可以发现，CS-2卫星数据在最优阈值组合下得到的海冰干舷高度，无论是反演精度还是相关性，均高于S3卫星数据。

表4 Sentinel-3A卫星数据不同阈值组合计算的海冰干舷高度与OIB海冰干舷高度对比

进一步分析图3，可以发现，无论是哪种阈值组合（即便采用本文得到的最优阈值组合），均普遍存在雷达高度计探测的海冰干舷高度大于机载OIB探测的海冰干舷高度的情况。这说明对于Ku波段雷达，积雪层不能视为透明，需要考虑其无法穿透积雪层的情况。当积雪覆盖在一年冰上时，海冰在析盐过程中形成的薄盐水层向上渗入到积雪层中，从而改变了积雪的介电特性，降低了雷达信号的积雪穿透能力[31]。对于覆盖在多年冰上的积雪，不仅积雪表面的粗糙度更高[32]，而且经过多年的消融冻结循环过程，会形成包含多个不同密度层的积雪层[33]。这些因素综合在一起，一方面会增加积雪层的后向散射强度，另一方面还会使雷达信号的主散射面由雪-冰交界面上移至积雪层中。这将导致海冰干舷高度被高估，从而引起海冰厚度的高估[14]。

3 积雪层对海冰干舷高度反演的影响分析

为详细评估积雪层对雷达信号的影响，计算了CS-2与S3的雷达穿透系数f[13]：

式中：d为积雪厚度hs与式（2）中的积雪校正值0.22hs之和；FR_sat为雷达高度计测的雷达海冰干舷高度；FR_OIB为OIB得到的雷达海冰干舷高度；OIB雷达海冰干舷高度由OIB的海冰干舷高度减去积雪校正值0.22hs得到。式（3）中，当雷达穿透系数f=0时，表示雷达信号没有穿入雪层，雷达的主散射面位于空气-雪交界面；f=1表示雷达信号完全穿透雪层，雷达信号的主散射面位于雪-冰交界面。0＜f＜1表示雷达信号未完全穿透雪层，雷达信号的主散射面位于积雪层内部。

图4和表5给出了2017年3月-4月和2018年4月北极区域CS-2与S3各自的雷达穿透系数分布与统计结果。在图4中还利用高斯分布对雷达穿透系数进行了拟合，并将雷达穿透系数的平均值作为拟合曲线的中心值。

图4 雷达穿透系数分布直方图

表5 CryoSat-2和Sentinel-3A的雷达穿透系数平均值

需要注意的是，图4存在雷达穿透系数f＜0或者f＞1的情况，对于该问题Armitage和Ridout[13]已给出过分析和评估。这主要是由于：（1）在计算海冰干舷时，为了有足够的覆盖范围与机载OIB数据进行比较，卫星数据的平均时间是远长于机载OIB数据的；（2）卫星高度计和机载OIB系统的足印大小存在较大差异。由于卫星海冰干舷与机载OIB海冰干舷的两种测量差异，使得卫星高度计的海冰干舷分布范围大于机载OIB的海冰干舷分布范围，从而出现雷达穿透系数f不总在[0,1]范围区间的情况。

从图4和表5可知，Ku波段雷达高度计对一年冰上积雪的穿透能力普遍强于多年冰。例如：CS-2对于一年冰上积雪的雷达穿透系数f为0.950，对于多年冰为0.889；S3对于一年冰上积雪的雷达穿透系数f为0.873，对于多年冰为0.856；对于全部海冰类型CS-2的雷达穿透系数f为0.912，S3的雷达穿透系数f为0.867。总体而言，CS-2的雷达穿透系数f大于 S3。

所以普遍来说，对于Ku波段雷达高度计，其主散射面均位于雪-冰交界面上方的积雪层中，不能再利用公式（2）进行积雪校正，必须发展新的校正方法。

4 积雪校正新方法与海冰干舷高度反演

由雷达穿透系数f的定义可知，当0＜f＜1时，雷达高度计测得的海冰干舷高度会比实际的海冰干舷高度高(1-f)hs，所以需要补偿这一距离差。同时，还需对雷达信号穿透的那部分积雪厚度进行距离校正0.22fhs。因此，改正的积雪校正模型可表述如下[13]：

式中：雷达穿透系数f可取为表5中的统计值，所以式（4）可根据雷达高度计和海冰类型的不同，分别针对CS-2与S3，得到本文新改进的积雪校正公式。

对于CS-2一年冰：

对于CS-2多年冰：

对于S3一年冰：

对于S3多年冰：

为了评估新改进的积雪校正方法的海冰干舷高度反演性能，本文以机载OIB实测海冰干舷高度数据为基础，将其与改进的积雪校正法的海冰干舷反演结果、通用积雪校正方法（公式（2））的反演结果和ESA发布的海冰干舷高度产品进行了对比分析。ESA海冰干舷高度产品是欧空局利用CS-2和S3数据制作的海冰干舷高度产品（L2级），该数据可公开下载（ftp://science-pds.cryosat.esa.int）。

表6 通用积雪校正方法、提出的积雪校正方法、ESA海冰干舷高度产品与机载OIB数据的海冰干舷高度对比

表6和图5给出的是3种海冰干舷反演方法的结果与机载OIB海冰干舷高度的差异统计表和分布图。根据图与表的统计结果可知，对于ESA发布的海冰干舷高度产品，CS-2发布的海冰干舷高度产品值明显偏大，而S3发布的海冰干舷高度产品值明显偏小，但总体上S3海冰干舷高度产品的误差低于CS-2。对于通用积雪校正方法，CS-2和S3反演的海冰干舷高度均偏高。而本文改进的积雪校正方法显著校正了通用积雪校正方法存在的高估问题，CS-2和S3海冰干舷高度的反演精度无论是平均绝对值差值还是均方根误差均小于前两种方法。另外需要指出的是，本文所改进的积雪校正方法对于一年冰和多年冰，其精度均有所提高，呈现了良好的稳定性。

图5 海冰干舷高度差异分布图

5 结论与讨论

工作在Ku波段的CryoSat-2和Sentinel-3A合成孔径雷达高度计是当前最先进的高度计。本文针对雷达回波重跟踪阈值的位置，以及积雪对雷达信号穿透能力的干扰影响，发展了改进的积雪校正算法，提高了海冰干舷高度反演精度。主要结论有：

（1）对于CS-2雷达高度计最优的波形重跟踪阈值组合为70%（冰间水道）、60%（海冰）；对于S3雷达高度计，其最优的阈值组合为50%（冰间水道），50%（海冰）。

（2）Ku波段雷达信号未完全穿透一年冰与多年冰上的积雪层，并且CS-2的雷达穿透系数大于S3。CS-2对于一年冰上积雪层的雷达穿透系数为0.95，对于多年冰为0.889；S3对于一年冰上积雪层的雷达穿透系数为0.873，对于多年冰为0.856；对于全部海冰类型而言，CS-2的雷达穿透系数为0.912，S3的雷达穿透系数为0.867。

（3）相比于通用积雪校正法和ESA产品法，在本文的改进积雪校正法中得到的海冰干舷高度与OIB数据之间的误差最小，精度最高。

本文的研究仅使用了北极春季的数据，为提高算法的普适性，未来将会利用北极其他季节的数据开展雷达穿透系数分析，从而进一步优化算法。在下一步的研究工作中，准备利用本文改进的方法开展南极地区的海冰干舷高度反演与校正。

致谢：感谢欧洲空间局提供的CryoSat-2和Sentinel-3A雷达高度计数据及海冰干舷高度数据；感谢美国冰雪中心提供Operation IceBridge数据。