王立伟 金涛勇 张胜军 李大炜

1 武汉大学测绘学院，武汉市珞瑜路129号，430079

CryoSat－2卫星是继ERS－1、ERS－2、Envisat及ICESat后首个针对极地进行观测的卫星，在高纬区域具有30d的子周期［1］。卫星搭载Ku波段雷达高度计SIRAL，有3种测量模式：LRM（低分辨率测量模式）、SAR（合成孔径雷达测量模式）、InSAR（干涉合成孔径雷达测量模式）［2］。考虑到数据源及数据处理方式的一致性，本文选择基于SAR 模式的观测数据。为顾及数据的现势性，使计算出的海冰干舷高能更好地反映某一时刻某一区域的特征，计算过程中采用沿轨内插方法，而没有采用常用的区域格网化方法。将两种方法的计算结果进行比较发现，二者在多年冰区域存在略微差异，但整体趋势较为一致。

1 CryoSat－2卫星波形重跟踪

在近岸区域或海冰区域，卫星测高波形通常会受到较为严重的污染，在数据处理之前，首先需要进行波形重跟踪，计算出星载预设门与重跟踪门之间的距离，进而对测量得到的卫星到星下点之间的距离进行改正。常见的波形重跟踪方法包括重心偏移法（OCOG 法）、阈值法、改进阈值法、神经网络算法、Beta－5参数算法、Beta－9参数算法以及海洋算法等［3－9］。

欧空局所发布的SAR 数据分为3 个级别：L1B数据，即波形数据；L2数据，直接给出相对于WGS84椭球的表面高；L2I数据，在L2数据基础上增加海冰区域的表面识别信息，以便计算海冰干舷高（海冰露出水面的高度）。由L1B 数据得到的L2 数据和L2I数据的处理方式统一采用OCOG 波形重跟踪算法，而OCOG 算法最初是针对海洋波形提出的，并不适合海冰的返回波形，而阈值法更适合于海冰区域的波形重跟踪［1，8－11］。

阈值法的思路是：首先计算或者直接利用OCOG 方法获得波形的振幅，然后对振幅乘以一个阈值因子，获得重跟踪点的幅值，最后线性内插获得重跟踪点的采样位置（图1）。计算公式如下：

图1 阈值法重跟踪示意图Fig.1 Sketch map of threshold retracking method

式中，p（i）为第i个采样点所对应的回波功率；Amax为回波功率的最大幅值，即波形振幅；DC 为利用前c个采样点所计算的噪声水平；TL为重跟踪点所对应的幅值；t为幅值大于TL值点所对应的最近采样门；LEG 为通过线性插值所获得的波形前缘中点，即波形重跟踪点所对应的采样门［12］。需要说明的是，很多学者在计算波形的噪声水平时通常采用波形的前几个采样值，一般默认从第一个采样值开始。有的测高卫星由于第一个采样值或前几个采样值存在较大偏差，实际计算时是利用后面的一段采样值对噪声水平进行估计［8］。本文在实际计算时利用每个波形的前5个采样值来估计噪声水平，即a＝0，c＝5。

本文采用阈值法对CryoSat－2卫星在海冰区域的SAR 波形进行重跟踪。由于雪密度和雪湿度等的影响，实际中雷达波形并不能完全穿透雪深。一些研究者将阈值系数直接选为50%，是假设雷达波形能够完全穿透雪深，使得计算出的冰层厚度存在较大误差。本文将阈值系数选为40%，即假设雷达波形只能穿透一定深度的雪深［1］。

2 海冰干舷高反演

2.1 海冰干舷高计算原理

计算海冰干舷高的关键是要获得高精度海冰表面高和可靠的瞬时海面高，通常采用格网化的方法，计算格网内的平均海冰高度和瞬时海面高度［8］。高精度海冰表面高可通过重跟踪方法获得，平均瞬时海面高一般取为格网内浮冰之间的海面高度。因此，对于海冰干舷高反演来说，最重要的就是进行浮冰之间的海水面（Lead）识别，准确测量出浮冰之间的局部瞬时海面高，然后内插计算海冰处的局部瞬时海面高。

格网法的缺点是，一个格网内的数据可能会含有多条轨迹，几条轨迹间的时间跨度可能很长，通过取中值或平均值所得结果并不一定代表某一时刻的瞬时海面高，从而影响海冰干舷高的计算精度。

本文采用单条轨迹的沿轨内插方法计算海冰位置处的瞬时海面高。先将DTU10平均海平面高内插到卫星的每一个星下点，然后通过波形特征参数识别出海冰的位置。以每一个海冰位置为节点，向轨迹的两端同时搜索，识别出与该海冰节点距离最近的Lead。如果两端均搜索到Lead，则同时计算两个Lead处海平面异常，并采用距离加权方法将其内插到海冰位置处，再加上海冰处的平均海平面高得到海冰处瞬时海平面高；如果仅一端搜索到Lead，则直接将该Lead处的海平面异常视为海冰位置处的海平面异常。事实上，仅一端能搜索到Lead的情况非常少见，可以忽略。

2.2 波形识别方法

为了进行Lead和海冰的识别，将OSI SAF（海洋与海冰监测）高纬处理中心发布的海冰浓度格网数据内插到CryoSat－2数据中，删除那些海冰浓度小于70%的数据（主要是海洋数据），然后再利用脉冲峰值、波束标准差、波峰陡峭度这3个波形特征参数对Lead和海冰进行区分。各参数的门限值设置情况如表1。

表1 波形识别参数设置Tab.1 Set the waveform parameter identification

脉冲峰值是波形回波功率的函数［1］：

pp值越大，说明波形越尖锐。波束标准差用来描述波形后向反射角的变化［2］，而波峰陡峭度用来描述波形顶端的陡峭程度［13］。

2.3 实验结果及比较

德国AWI利用40%阈值法分析2011－01～2013－03共27个月的数据，利用格网化方法计算得到海冰干舷高。选取其最后3个月的数据与本文结果进行比较，如图2。可以看出，二者得到的北冰洋海冰干舷高整体分布趋势较为一致，AWI在多年冰区域的海冰干舷高要更高一些。图3给出2013－01～03 二者的差值直方图。可以看出，这3个月二者的差值直方图非常接近于正态分布。最后利用OSI SAF发布的每天多年冰和一年冰格网数据，分别计算AWI和本文数据的多年冰区域和一年冰区域的平均海冰干舷高（图4）。

从图4看出，本文计算的2011～2013年北冰洋多年冰和一年冰区域的平均海冰干舷高结果与AWI较为一致，海冰干舷高呈现明显的周期性变化。由于夏季海冰融化，使得波形识别受到影响，利用返回波形特征识别Lead和海冰的方法不再可用，所以没有对6～8月的海冰进行处理。AWI利用波形识别方法得到的6～8月北冰洋海冰干舷高很不理想。对于夏季海冰的识别，还需采用其他更为有效的方法。

图5为2011～2013年北冰洋各年同期数据对比图，无论是本文结果还是AWI结果，各年同期海冰干舷高结果均较为接近。但同时也要注意到，在2～5月，本文得到的2013年数据无论是在多年冰区域还是在一年冰区域都较同期其他数据偏低，而在9～12月略微偏高，尤其是2013－09的数据 较2011－09和2012－09呈现明显的增长趋势。由于AWI数据只公布到2013－03，所以为了验证本文结果，还需要采用实地测量数据或其他卫星测量数据进行更好的解释和说明。

图2 本文海冰干舷高计算结果与AWI计算结果的比较Fig.2 Comparison between the results of sea ice freeboard calculation in this paper and the results of AWI

图3 2013－01～03AWI与本文海冰干舷高差值直方图Fig.3 Histogram of the sea ice freeboard height difference（2013－01～03）

3 结 语

本文采用沿轨内插方法计算北冰洋区域的平均海冰干舷高，与常用的格网化方法计算结果整体趋势较为一致，说明沿轨内插方法完全可以用来进行海冰干舷高的反演。在计算中发现，2013－09以后北冰洋区域无论是多年冰区域还是一年冰区域，平均海冰干舷高均呈现明显上升趋势，这还需要通过其他数据进行验证。为了对本文结果和AWI结果进行更好的比较，下一步将采用实测数据和机载数据，对二者的计算结果进行更详细的比较，进而评价格网化方法和沿轨内插方法在计算海冰干舷高中的优劣，最终将两种方法合理组合，更好地反演海冰干舷高。

图4 2011～2013年北冰洋平均海冰干舷高变化趋势Fig.4 Variation trend of the average sea ice freeboard height from 2011to 2013year

图5 2011～2013年北冰洋平均海冰干舷高同期数据比较Fig.5 Arctic sea ice freeboard period average data comparison from 2011to 2013

［1］Wingham D J，Francis C R，Baker S，et al.CryoSat：A Mission to Determine the Fluctuations in Earth’s Land and Marine Ice Fields［J］.Adv Space Res，2006，37（4）：841－871

［2］Ricker R，Hendricks S，Helm V，et al.Sensitivity of CryoSat－2Arctic Sea－Ice Volume Trends on Radar－Waveform，Interpretation［J］.The Cryosphere Discuss，2014，8（2）：1 831－1 871

［3］高永刚，郭金运，黄金维，等.湖泊范围内TOPEX／Poseidon卫星测高波形分析［J］.大地测量与地球动力学，2008，28（1）：45－49（Gao Yonggang，Guo Jinyun，Huang Jinwei，et al.Waveform Analysis of Topex／Poseidon Satellite Altimeter over Lake［J］.Journal of Geodesy and Geodynamics，2008，28（1）：45－49）

［4］褚永海，李建成，张燕，等.ENVISAT 测高数据波形重跟踪分析研究［J］.大地测量与地球动力学，2005，25（1）：76－80（Chu Yonghai，Li Jiancheng，Zhang Yan，et al.Analysis and Investigation of Waveform Retracking Data of Envisat［J］.Journal of Geodesy and Geodynamics，2005，25（1）：76－80）

［5］褚永海.卫星测高波形处理理论研究及应用［D］.武汉：武汉大学，2007（Chu Yonghai.The Research of Satellite Altimetry Waveform Theory and Application［D］.Wuhan：Wuhan University，2007）

［6］常晓涛，李建成，郭金运，等.一种多前缘多阈值的波形重构算法［J］.地球物理学报，2006，49（6）：1 629－1 634（Chang Xiaotao，Li Jiancheng，Guo Jinyun，et al.A Multileading Edge and Multi－threshold Waveform Retrackef［J］.Chinese J Geophys，2006，49（6）：1 629－1 634）

［7］European Space Research Institute（ESRI N）.European Space Agency and Mullard Space Science Laboratory－University College London［EB／OL］.http：／／earth.esa.int／web／guest／missions／cryosat／ipf－baseline，2013

［8］Rose S K.Measurements of Sea Ice by Satellite and Airborne Altimetry［D］.Denmark：DTU Space，National Space Institute，2013

［9］Davis C H.A Robust Threshold Retracking Algorithm for Measuring Ice－Sheet Surface Elevation Change from Satellite Radar Altimeters［J］.IEEE Trans Geosci Remote Sens，1997，35（4）：974－979

［10］Helm V，Humbert A，Miller H.Elevation Change of Greenland and Antarctic Derived from CryoSat－2［J］.The Cryosphere Discuss，2014，（8）：1 673－1 721

［11］Stenseng L，Andersen O B.Preliminary Gravity Recovery from CryoSat－2Data in the Baffin Bay［J］.Adv Space Res，2012，50（8）：1 158－1 163

［12］阮海林.卫星测高波形重跟踪算法在台湾周边海域的应用研究［D］.厦门：国家海洋局第三海洋研究所，2010（Ruan Hailin.Research and Application of Algorithm in the Waters Surrounding Taiwan Satellite Altimetry Waveform Retracking［D］.Xiamen：The Third Oceanography Research Institute，2010）

［13］Brenner A C，Zwally H J，Bentley C R，et al.Derivation of Range and Range Distributions from Laser Pulse Waveform Analysis for Surface Elevations，Roughness，Slope，and Vegetation Heights［EB／OL］.http／／glas.wff.nasa.govl，2000