1 引言

卫星雷达测高是现代大地测量学的重要技术手段之一，已成功用于确定全球海面高和海洋重力异常等应用^{[1, 2]}。但是在近海，由于卫星雷达测高技术自身的局限和该区域复杂的地理环境，卫星测高的观测数据质量较差，限制了卫星测高在近海的应用。文献^[3]分析澳大利亚近海的ERS-2和Poseidon数据后发现：距海岸22 km以内的海面高数据均受到不同程度的污染，而且距离海岸小于4 km的数据几乎不能使用。若按照通常的数据处理方法和编辑标准，近岸大部分测高数据会在预处理过程中被删除^[4]，导致近海区域可用的数据很少甚至空白。

改善近海卫星测高数据质量需要从两方面进行，其一是改进近海区域海潮、大气延迟等地球物理和环境改正的精度；另一方面是提高卫星至海面距离观测值的确定精度。后者主要通过对测高回波波形进行重跟踪实现。为提高测高观测值精度，研究人员提出了多种波形重跟踪方法，应用较多的有Beta-5参数法^[5]、重心偏差法^[6]和阈值法^[7]等，基于这些方法又提出了许多改进的算法^{[8, 9, 10, 11, 12]}。虽然波形重跟踪技术已被证明能够有效改善测高数据的精度，但是没有一种方法能够使海洋环境极为复杂的近海地区各种测高波形的重跟踪效果均达到最优。这是因为各种波形重跟踪方法都针对特定地区的应用提出的，只对某些特定的波形有较好的重跟踪效果，例如Beta-5参数法适用于海洋波形、阈值法是为了处理冰面波形提出的。地理环境和海况复杂的近海存在多种多样的波形，仅用一种算法难以满足各种波形的重跟踪，也就不能最大限度地提高测高数据的质量。若能根据波形形状自动选取最佳重跟踪算法，则可以极大提高近海可用测高数据的质量和数量。文献^[13]曾开发了一个专家系统用于陆地测高波形数据的处理，据称该系统集成了十几种重跟踪算法，能够对不同波形进行自动处理。基于此思想，本文以Topex/Poseidon(T/P)卫星数据为例，在对近海测高波形的进行有效分类的基础上，提出了一种根据波形的形状选择最佳重跟踪算法的自适应方法。 2 波形分类 2.1 聚类方法与有效性指标

测高波形是雷达回波功率随时间变化的曲线，由测高仪的自动增益装置在设定的时间窗内对回波功率进行连续采样得到，其形状与反射面的性质直接相关。对测高波形进行准确的分类是合理选择重跟踪方法的关键，也可以用于研究反射面的性质，例如用于研究南极海区的海冰密集度^[14]。测高波形可以根据回波功率分布特征(如脉冲峰值、波形振幅、波宽等)^{[14, 15, 16, 17]}或者后向散射系数和有效波高等参数^[18]来进行分类，这些方法一般将测高波形分为两类：漫射波形和镜面波形，或者是海洋波形和非海洋波形。在测高波形分类数未知的情况下，聚类分析方法可以对波形进行有效分类^{[19, 20]}。本文采用文献^[19]提出的以移动最小欧氏距离作为相似性统计量、以最小方差法进行聚类的分类方法，这里对距离计算公式做了改进

点击浏览公式

利用聚类分析进行分类最关键的问题是确定最佳聚类数，即聚类有效性分析。聚类结果的质量可以从类内紧密度和类间分离度两个方面来衡量，好的聚类结果应使类内距离尽可能小，而类间距离尽可能大。为了使类内紧密度和类间分离度两者达到一定平衡，聚类有效性分析通常采用这2个量的组合作为有效性指标，用来评价聚类的效果和确定最佳聚类数。本文在确定波形的最佳分类数时，根据文献^[21]设计了一个不依赖于聚类算法、只与聚类的划分有关的有效性指标。

假设给定样本集的一个划分C^k={C₁,C₂,…,C_k}，k即为聚类数，将类内紧密度Scat(C^k)和类间分离度Sep(C^k)分别定义为

点击浏览公式

点击浏览公式

点击浏览公式

要使聚类质量最优，类内紧密度和类间离散度应达到平衡，即有效性指标Q(C^k)取得最小值，由此可确定最佳的聚类数k。 2.2 T/P卫星近海测高波形分析

以穿过台湾海峡的T/P卫星第6周第164号轨迹为例分析近海测高波形的分类。该轨迹为降轨道，在台湾海峡西岸卫星轨迹从福建长乐市南部入海，穿过海坛海峡和海坛岛，进入台湾海峡，最后在海峡西岸的台湾彰化县登陆，如图 1所示。利用GMT软件的gmtselect命令去除陆地轨迹点后该段轨迹共有308个10 Hz点，即有308个测高波形。由于各个测高波形的振幅不同，先采用极值归一化方法将所有波形规格化，以规格化后的波形为样本集进行聚类分析^[19]。本文使用的T/P卫星波形数据为压缩后的波形，每个波形有64个功率采样(N=64)，在利用公式(1)计算距离值时根据经验将k₀取8。为了确定最佳聚类数k，分别计算了将样本集分为1类至10类时的有效性指标Q值，结果如图 2所示。当聚类数k=3时，有效性指标Q取得最小值0.083 5，说明将波形分为3类最合理。

图 1 T/P卫星第6周第164号轨迹及波形分类 Fig. 1 Sample track of T/P waveform data(cycle 6, pass 164) and waveform classifications

图选项

图 2 有效性指标变化曲线 Fig. 2 Variation curve between validity index and cluster number

图选项

对样本集划分为3类的聚类结果为C={C1(240),C2(39),C3(29)}，圆括号内数值是类中的波形样本个数。图 1给出了3类波形的空间分布，其中轨迹线上方浅灰色表示的是C1类，轨迹线上用黑色圆点表示的属于C2类，轨迹线下方为黑色则属于C3类。图 3显示的是每一类的聚类中心，由于测高波形的前缘位置可能存在平移，所以不适合采用平均值作为聚类中心，本文采用平均距离最小原则确定聚类中心，即把类中到其他样本的平均距离最小的样本作为聚类中心。聚类中心反映了各类波形的主要特征，C1类的中心表现为典型的海洋波形，该类波形主要分布在台湾海峡中间的开阔海域。C3类的波形呈尖锥状，类似镜面反射，回波功率上升和下降的速度都很快，大多出现在海-陆转换区域。上述2类中各波形的形状比较规则单一，而C2类中的波形则较为复杂，一些不规则的波形均被划分到此类。由图 1可以看出C2类波形主要出现在卫星轨迹从陆地到海洋的转换带上，可能因为存在多个不同高度的反射面所致，这一类的聚类中心也表现出具有多个前缘的特征。

图 3 聚类数k＝3时的聚类中心 Fig. 3 Cluster centers for cluster number k=3

图选项

为了进一步验证分类的效果以及聚类中心的代表性，利用最小距离判别法对样本集重新分类。计算每个样本波形分别到3个聚类中心的距离，将样本归为距离最小的一类。重新分类的结果为C={C1(242),C2(33),C3(33)}，与聚类分析的结果差别很小，各类波形的空间分布也与图 1类似。可见，利用聚类分析方法得到的波形类型划分是可靠的，特别是海洋波形与非海洋波形之间能够准确的区别开。同时也表明得到的聚类中心对于该轨迹上的测高波形具有很好的代表性，利用聚类中心可以简单根据最小距离判别法就能有效进行波形类型的识别和分类。 3 测高波形的自适应重跟踪方法 3.1 各类波形的最佳重跟踪方法

现有的波形重跟踪方法都是针对不同的测高环境提出的，它们分别有各自的适用情况。重心偏差法和阈值法属于统计方法，没有顾及波形形状及其物理意义，所以这两种方法几乎能够成功解算所有类型的波形，但是精度得不到保证。Beta-5参数法是严格的数学方法，利用一个含有5个参数的标准数学模型去拟合实际波形，通过迭代计算解算出这5个参数。相对于前两种方法，Beta-5参数法具有较高的精度，而且这些参数都具有物理意义，但是它只对接近于标准波形的海洋波形有效。对于近海浅海区域，回波波形与标准波形相差较大，Beta-5参数法在迭代计算时难以收敛，不能成功解算出所有的5个参数，导致波形重跟踪失败，因此，这种方法的成功率较低。为了改善近海测高数据，针对多前缘的测高波形提出了一类改进的阈值算法，取得了较好的效果^{[9, 10, 11, 12]}。根据前面波形分类的结果，这里对各类波形分别采用重心偏差法、50%阈值法、改进阈值法和Beta-5参数法进行波形重跟踪计算，以确定适合各类波形的最佳重跟踪方法。

以EGM08重力场模型计算的大地水准面高为参考值评价重跟踪后的海面高精度。计算海面高时进行了各项环境、地球物理改正以及高程基准的统一，其中湿对流层延迟改正和海潮改正分别采用ECMWF和NAO99b模型计算值。海面高中的海面地形利用MDT_CNES-CLS09平均海面地形模型^[22]进行了扣除。表 1给出了各类波形重跟踪前以及不同重跟踪方法改正后的海面高与EGM08大地水准面高差值的标准差，圆括号内数字为参加统计的点数。这里目的只是评价波形重跟踪方法的效果，故海面高数据没有按GDR数据处理的编辑准则进行预处理，而是保留了全部海面高数据，所以有些数据存在一些比较大的误差导致统计出标准差比较大。例如C3类波形经过重跟踪后的标准差扔接近5 m，在确定平均海面高的应用时，这些数据若不进行特殊处理通常都是当作粗差值被剔除。产生这些大误差的原因除了测距误差外，很重要的一个因素是近岸海潮模型误差。表 1的结果显示，各类波形经过重跟踪后的海面高与重跟踪前的海面高相比精度都有较大的提高，但是对于不同类型的波形每种重跟踪方法的改善程度不同。对于C1类波形，Beta-5参数法的效果明显优于其他方法，其次是改进阈值法。对于C2、C3两类波形，Beta-5参数法重跟踪的成功率很小，显然不适合这两种波形。C2类波形的重跟踪结果改进阈值法精度最高，而C3类波形各种方法的效果相差不大，其中重心偏差法略微占优。

根据上面的分析结果，提出一种近海波形重跟踪的自适应算法，根据波形自动选择最合适的重跟踪方法，图 4给出了算法流程图。将T/P卫星近海测高波形分为3类，首先计算每个波形与类中心波形的距离，采用最小距离判别法判断该波形属于哪一类，再选择相应的方法进行重跟踪。若波形属于C1类，选择Beta-5参数法，如果Beta-5参数法重跟踪失败则选择改进阈值法。若波形属于C2或C3类，则分别选择改进阈值法、重心偏差法。

图 4 自适应重跟踪算法流程图 Fig. 4 Flow chart of adaptive retracking algorithm

图选项

利用该算法对整条轨迹308个波形进行重跟踪计算，图 5给出了重跟踪改正后的海面高(点折线)。作为比较，重跟踪改正前的海面高(带叉折线)和EGM08大地水准面(黑线)也绘于图 5中。尽管轨迹两端近岸的海面高仍偏离大地水准面很大，重跟踪后海面高的精度无论在海峡中央还是在两岸附近都有明显改善。从图中可以看出，近岸的海面高在大地水准面高上下振荡，表现出波的特性，这说明近岸复杂的海况引起的海面高误差是主要误差源。表 2给出了各种重跟踪方法的精度对比，Beta-5参数法不能应用于所有波形，有效点数比其他方法少很多，故未对其精度进行统计。由于海面高数据之前未经预处理编辑，在进行精度统计之前先将与大地水准面之差大于11 m的数据作为粗差剔除，这个剔除标准是以T/P卫星GDR数据手册中有效波高的编辑标准为参考选取的。表中给出的精度是用粗差剔除后的海面高与相应大地水准面高的差值计算的标准差，剔除粗差后的有效点数列于标准差后的括号内。除了第6周期的数据外，还随机选取了其他5个周期的数据进行实验，结果同样列在表 2中。统计结果表明改进阈值法是仅采用一种重跟踪方法时整体精度最好的，明显优于重心偏差法和阈值法。从第193周期的结果还显示，如果仅采用某一种重跟踪算法，有时虽然能增加有效数据，但是海面高的整体精度却可能比未作重跟踪改正时还差。而本文提出的自适应重跟踪算法综合了各种重跟踪方法，使每一类波形的海面高获得最优的重跟踪改正，海面高的精度比重跟踪改正前有分米级的提高，也显著优于其他单一方法，且有效点数与其他方法相当。

图 5 重跟踪前后的海面高比较 Fig. 5 Comparison of sea surface heights before and after retracking

图选项

本文在讨论测高波形分类的基础上，提出了一种根据波形自动选择最佳重跟踪算法的自适应重跟踪方法，可以更加有效地提高近海卫星测高数据的精度。本文虽然只对T/P卫星近海测高数据进行分析，但是该方法对其他卫星数据同样适用，而且可以推广应用于内陆湖泊测高数据的处理。该方法的基本思路是：先对样本数据进行聚类分析，确定最佳的聚类数和相应的聚类中心；然后针对每类波形比较各种重跟踪算法的精度，确定适合该类波形的最优重跟踪算法；根据前面两步的分析结果制定自适应重跟踪算法，即利用聚类中心判断波形属于哪一类，然后选择适用该类波形的最佳重跟踪算法计算海面高改正量。

卫星测高近海波形极为复杂，细分可达10余种之多^[23]，但是对台湾海峡T/P卫星测高波形数据的聚类分析表明，T/P卫星近海测高波形分为3类时聚类效果最佳。除海洋波形和尖锥状波形外，其他复杂波形都归为多前缘波形这一类。海洋波形和多前缘波形分别采用Beta-5参数法和改进阈值法进行重跟踪效果最优，尖锥状波形除了不能用Beta-5参数法以外其他方法的效果相差不大。后两类波形的海面高虽然经重跟踪改正后精度有明显提高，但仍然存在较大误差，主要是由于海潮改正等环境改正项不准确所致，因此，这些数据若不经进一步的特殊处理应谨慎使用。由于算例在波形分类时只用了1条轨迹的数据，数据样本量较小，所以分类的结果可能不能完全代表近海所有的测高波形，这一问题需要加大样本量作进一步的分析。

致 谢：宁津生院士是我国著名大地测量学家，为我国测绘事业发展和人才培养作出了突出贡献。值此宁津生院士80华诞之际，特撰写此文以表达对他的崇高敬意，并感谢他对笔者多年的培养！