基于双调和样条内插和高斯曲率极值的多面函数拟合高程异常方法

孙佳龙，焦明连，郭淑艳

(1. 淮海工学院 测绘工程学院，江苏 连云港 222001；2. 江苏省海洋资源开发研究院，江苏 连云港 222005)

一、引 言

随着GPS测量技术的不断完善和发展，求取高精度的高程异常在实际工程应用中越来越重要[1-3]。在拟合区域较小的情况下，一般的拟合方法可以达到较高的精度，但由于这些方法未顾及似大地水准面的物理性质，拟合出的曲面只能是高程异常的趋势面，与高程异常的实际数值存在一定差异[4-7]。而多面函数在理论上可以以任意精度逼近任意复杂曲面，因此可以更加逼近实际的高程异常曲面[8]。而当GPS/水准点数量较多时，如果随机选取拟合点，则选择的方案数过多，可能超出普通计算机的计算能力。鉴于高程异常曲面具有均衡性，且在数学上可用调和函数来描述，本文提出了基于双调和样条及高斯曲率极值的多面函数拟合方法，通过计算高程异常曲面网格数据的高斯曲率及极值，选取距离高斯曲率极值点最近的GPS/水准数据点作为多面函数的结点，参与高程异常曲面的拟合，而剩余点作为检核点，对拟合的曲面进行外部检核。

二、双调和样条函数

对于一组不规则的空间数据点，通常可以利用三次样条函数找出最平滑的曲线或曲面。从力学的角度看，该过程就是强迫一个弹性杆或弹性板去匹配各个数据点[9]。

样条的格林函数的点压力满足双调和方程[10]

(1)

式中，x为空间位置；δ(x)是Delta函数；格林函数Φ(x)就是方程式(1)的一个解，表示内插曲线或曲面的位移，是空间位置x的函数，可以写成

(2)

当该格林函数对N个在xi的数据点ωi进行插值时，式(1)就成为

(3)

ω(xi)=ωi

(4)

式中，ωi是内插曲线或曲面在空间位置xi处的位移；αj为加权系数。

三、基于高斯曲率极值的结点判定方法

在曲面论中，Weingarten公式刻画了曲面单位法向的运动状况，可以进一步了解曲面的弯曲状况。Weingarten公式的矩阵形式为[11]

(5)

式中，矩阵ω为曲面在参数(u1,u2)下的Weingarten矩阵。ω可以表示为

曲面上一点的法曲率关于切线方向的两个最值分别为曲面在该点处的主曲率，而使法曲率达到最值的两个切线方向分别为曲面在该点处的主方向。曲面的主曲率和主方向分别为Weingarten变换的特征值和特征方向。

曲面在一点处的两个主曲率的乘积K称为在该点处的高斯曲率。高斯曲率用Weingarten矩阵表示为[12]

对于曲面上某分片区域的一点，如果该点的2个主曲率(K1和K2)中任何一个沿着对应的主方向上为极值，则该点称为曲率极值点。由于高斯曲率极值可以很好地刻画曲面的弯曲情况，因此，本文选取高程异常曲面上距离高斯曲率极值点最近的GSP/水准数据点作为最后参与拟合的数据点，即多面函数中核函数的结点，参与多面函数的拟合。

四、算例分析

某市E级GPS控制网，进行了三等水准测量，共包含71个GPS/水准点，本文以此数据为例，利用双调和样条函数插值出了该市的高程异常曲面，并计算了高程异常曲面网格数据的高斯曲率及极值点。采用距离高斯曲率极值点最近的原则对所有GPS/水准数据点进行筛选，利用剩余数据点对拟合的曲面进行检核。将筛选出的数据点作为结点进行多面函数拟合，而其余点作为检核点对拟合的高程异常曲面进行外部精度检核，各种数据点如图1所示。

在图1中，星号表示高斯曲率极值点，三角符号为所有的GPS/水准数据点，而三角符号加四边形符号为从所有GPS/水准数据中选择的参与拟合的数据点，即GPS/水准数据点中距离高斯曲率极值点最近的点。图中颜色条表示高程异常值，单位为m。从图1中可以看出，整个观测区域的高程异常曲面呈现北高南低的趋势，但在中部却又有一个低凹的异常区域出现，因此，观测区域的高程异常变化趋势比较复杂。所有的GPS/水准数据又呈现出南多北少的分布特点，特别是在东北区域，没有一个数据点。通过计算高斯曲率极值，即该点的高斯曲率均大于(极大值)或小于(极小值)周围4个网格点的高斯曲率，如图1中各种颜色的最深处和两种颜色之间的过渡处，表示曲面弯曲变化较大较快的情况，因此，利用高斯曲率极值可以很好地反映曲面的变化情况。由于计算出的高斯曲率均位于插值出的网格点上，而GPS/水准数据点一般并不与这些网格点重合，从图1中可以看出，高斯曲率极值点有时与GPS/水准点相距较近，而有时却很远，但是，距离这些高斯曲率极值网格点最近的GPS/数据点都可以在一定程度上反映曲面的弯曲情况，在图1中也可以看出，被选择的数据点也比较均匀地分布在观测区域里，因此，选择这些数据点参与曲面的拟合可以在一定程度上综合反映高程异常曲面的变化趋势。

图1 高斯曲率极值点、GPS/水准数据点及选择的拟合点

在71个GPS/水准数据中，共有44个高斯曲率极值点，利用高斯曲率极值点作为多面函数的结点对高程异常曲面进行拟合。当高斯曲率极值K>0时，多面函数采用椭圆抛物面函数进行拟合；而当高斯曲率极值K<0时，采用双曲抛物面函数进行拟合。以剩余的27个数据点作为检核点，二次多项式、移动二次曲面、基于双曲抛物面的多面函数(hyperbolic paraboloid of polyhedral function，HPPF)、椭圆抛物面的多面函数(elliptic paraboloid of polyhedral function，EPPF)及本文提出的基于高斯曲率极值的多面函数(Gaussian curvature extremum of polyhedral function，GCEPF)方法的拟合残差如图2所示。

从图2中可以看出，二次多项式拟合的残差序列波动性较大，在多个检核点的残差都较大，最大绝对值约为10 cm；而移动二次曲面和EPPF的残差序列的趋势虽较为平缓，但最大残差绝对值也接近6 cm；HPPF和本文提出的GCEPF方法得到的残差序列的整体趋势趋于平缓，最大绝对值约为4 cm。为了更细致地比较几种方法的优劣性，各种拟合方法拟合得到残差的最大值、最小值、平均值和标准差见表1。

图2 5种方法的拟合残差

表1 残差序列的统计特征 cm

从表1中可以看出，GCEPF方法残差最大值虽然比其他几种方法(EPPF除外)大1～2 cm，但最小值却比其他几种方法增加4～9 cm，使极值的误差控制在5 cm以内(其他4种方法均超过5 cm)，满足了一般测量工作对高程精度的要求。

HPPF和EPPF残差序列的标准差比二次多项式分别减小了27.2 mm和25.7 mm，但却比移动二次曲面增加了2.5 mm和4 mm，说明采用单一某种核函数的多面函数拟合的精度优于二次多项式，但却不及移动二次曲面。可见，由于移动二次曲面拟合方法顾及了拟合点周边的高程异常起伏的影响，因此拟合的曲面更加逼近真值，拟合精度较高。

而以高斯曲率极值为参考的GCEPF方法能更细致地识别高程异常曲面的突变情况，其残差序列的标准差比二次多项式减小了33.5 mm，比其他采用单一某种曲面的多面函数HPPF和EPPF减少了6.3 mm和7.8 mm，比移动二次曲面也减小了3.8 mm，因此，GCEPF残差序列更趋于平缓，更逼近真值，优于以上4种方法。

五、结束语

对于区域高程异常曲面，基于某种单一核函数的多面函数很难较好地进行拟合，而从某种意义上说，基于单一核函数的多面函数已经失去了多面函数的优势，而在实际情况中，想要找到适合不同区域的核函数，实现多个核函数的组合拟合，从而真正体现其“多面性”是比较困难的。

本文利用双调和样条函数对高程异常曲面进行插值，利用高斯曲率极值的条件选择核函数及其结点，既顾及了高程异常曲面的均衡性，同时又能有效地对核函数及其结点加以选定，从而实现了利用“多面”函数进行拟合的目的。最终的拟合结果也说明，本文提出的基于双调和样条函数及高斯曲率极值的高程异常拟合方法的拟合精度比以上4种方法得到的精度更高，更逼近真值。

参考文献：

[1] 陈传法, 岳天祥. 基于预处理共轭梯度法的高精度曲面模型解算[J]. 中国矿业学报, 2010, 39(2): 290-294.

[2] 程芦颖, 许厚泽. 顾及测站点上重力场信息的大地水准面高的拟合方法[J]. 武汉大学学报：信息科学版，2008, 33(7): 701-705.

[3] 车延国. GPS高程水准拟合模型与精度分析[J]. 辽宁工程技术大学学报, 2008, 27(1): 35-38.

[4] 刘大杰，姚连璧，周全基. GPS水准的拟合基准面高程系统[J]. 测绘学报, 2000, 29(S1):11-15.

[5] 胡伍生, 华锡生, 张志伟. 平坦地区转换GPS高程的混合转换方法[J]. 测绘学报, 2002,31(2): 128-133.

[6] 郭东美, 许厚泽. 应用GPS水准与重力数据联合解算大地水准面[J]. 武汉大学学报：信息科学版, 2011, 36(5): 621-624.

[7] 马洪滨, 董仲宇. 多面函数GPS水准高程拟合中光滑因子求定方法[J]. 东北大学学报：自然科学版，2008, 29(8):1176-1178,1191.

[8] 张菊清, 刘平芝. 抗差趋势面与正交多面函数结合拟合DEM数据[J]. 测绘学报, 2008, 37(4): 526-530.

[9] SANDWELL D. Bihamonic Spline Interpolation of GEOS-3 and SEASAT Altimeter Data[J]. Geophysical Research Letter, 1987, 14(2): 139-142.

[10] JUN Yongtae, CHOI Kuiwon. Automated Feature-based Registration for Reverse Engineering of Human Model[J]. Journal of Mechanical Science and Technology, 2005, 19(12): 2213-2223.

[11] 马骊溟, 徐毅, 李泽湘. 基于高斯曲率极值点的散乱点云数据特征点提取[J]. 系统仿真学报, 2008, 20(9): 2341-2344.

[12] RABBARI T, DIJKMAN S. An Intergrated Approach for Modeling and Global Registration of Point Clouds[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2007,61(2): 355-370.

[13] YANG M, LEE E. Segmentation of Measured Data Using a Parametric Quadric Surface Approximation [J]. Computer Aided Design, 1999, 31(7): 449-457.