基于移动二次曲面拟合的GNSS点高程异常计算

吴新辉

(漳州市测绘设计研究院，福建 漳州 363000)

包括北斗卫星导航系统在内的全球导航卫星系统(GNSS)及连续运行参考站(CORS)已逐步成为测绘等相关领域所采用的主要测量手段[1]。目前我国使用的高程是正常高，GNSS获得的高程为大地高，大地高数据对于高程控制测量只有几何意义，只有获得大地水准面精化模型成果，将测量所得的大地高转换为正常高，才有实际应用意义[2]。

GNSS直接获取的高程为大地高H，而我国采用正常高h，所以大地高数据对于高程控制测量仅具有几何意义。因此，在实际生产应用过程中，应采用正常高才能起到高程控制作用，所以如何将获取的大地高快速、准确地转换为正常高已成为测绘工作中的一项重要内容。

当前测绘生产中，一般采用静态测量或者精密单点定位技术(PPP)获得一些数量且分布合理的点的大地高，同时采用传统水准测量方式进行高程联测，然后拟合出该“空白区”的似大地水准面，即高程异常值，最后得到该区域其他待求点的正常高[3]。该方法被广泛应用于地形平坦地区，且精度可达到四等水准或更高[4]，但效率较低，人力、物力投入较多，生产周期长，对于地形起伏较大、时间紧迫的任务不适用。随着重力场模型的精度及分辨率越来越高，采用少量GNSS水准点进行高程异常拟合快速获取正常高成为可能。

EGM2008重力场模型在中国大陆区域具有一定的适应性，且精度优于其他模型[5]。本研究在EGM2008重力场模型的基础上，采用漳州市区实测C级控制网数据，比较了移动曲面拟合模型中平面法、二次曲面法、加权平均法这3种方法在残余高程异常计算中的应用，结果表明在EGM2008重力场模型的基础上，对残余高程异常进行拟合获得的高程异常值，可以满足一定范围内的测绘工作需求。

1 原理与方法

1.1 技术流程

正常高是以似大地水准面为基准面的，两者之间的差值为高程异常ζ，三者关系如下：

h=H-ζ。

(1)

当前建立一个区域的似大地水准面模型，流行的方法是移去-恢复法，即对测区范围内离散的控制点进行GNSS水准联测后，把控制点的高程异常ζ分解成中长波分量(ζGM)、短波分量(ζT)和残余部分(ζres)，其关系如下：

ζ=ζGM+ζT+ζres。

(2)

其中，用重力场模型高程异常来表示中长波部分高程异常，剩余地形模型计算短波分量高程异常。

在计算时先把中长波和短波部分移去，再对残余部分用数学或统计的方法进行拟合，然后将拟合结果与移去的中长波和短波部分相叠加，得到最终的高程异常，最后拟合出该区域的似大地水准面模型。

在似大地水准面精化研究中，刘瑛等[6]对多个重力场模型进行了对比，研究表明EGM2008重力场模型精度要优于其他模型，并肯定了其在中国大陆区域的适用性。本研究亦选用EGM2008重力场模型。

在不考虑短波分量部分时，可以将短波分量部分归入残余高程异常，即

ζ=ζGM+ζres。

(3)

在对残余高程异常部分进行拟合时，移动曲面拟合模型是日常测绘中较为常用的一种方法。

1.2 计算高程异常长波分量

EGM2008重力场模型是美国国家地理空间情报局(NGA)在充分利用GRACE卫星重力、TOPEX卫星测高、5′×5′地面重力等多种数据的基础上，研制发布的2 190次球谐系数的全球超高阶地球重力场模型。该模型的基本格网分辨率为5′，目前格网分辨率已经达到2.5′×2.5′和更加精细的1′×1′，是当前计算重力高程异常的理想模型。目前可以通过NGA网站下载格网分辨率1′×1′的数据，利用Alltrans EGM2008 Calculator软件可以直接计算高程异常ζGM。

1.3 残余高程异常项处理

残余高程异常项拟合方法较多，主要有二次曲面函数法、多面函数法、神经网络法、移动曲面法，每种方法都有其优缺点。比如：二次曲面法仅能满足地形起伏小且控制区域小的一般工程需求；多面函数由于其参数的选取而对拟合精度影响较大且工作量也大；神经网络法的算法复杂，在工程应用中使用神经网络法尚有诸多问题有待进一步研究[7]。考虑到模型复杂性、精度和可靠性，本研究的残余高程异常拟合采用距离定权的移动曲面模型。该模型是在每一个待求点的邻域半径R范围内拟合出一个曲面，同时引入距离定权的权函数，在该区域内，已知点与待求点距离越近权越大，反之越小[8-11]。移动曲面模型又包含平面法、二次曲面法和加权平均法。

1.3.1平面法

平面法模型如下：

ζresi=a0+a1xi+a2yi，i=0，1，2，…，n。

(4)

平面法要求已知点大于3个，利用最小二乘法的参数平差原理求解未知数a0、a1、a2，则式(4)的矩阵形式如下：

V=AX+L，

(5)

根据最小二乘原理VTPV可得

ATPAX+ATPL=0。

(6)

设ATPA=N、ATPL=U，则式(6)可改写成如式(7)的法方程：

NX+U=0，

(7)

式中：N为系数阵；U为常数项。当N为非奇异阵时，可以求出未知数X：

X=-N-1U。

(8)

将已知数据代入函数，解出未知数a0、a1、a2，这样就确定了一个拟合平面。把待求点的平面坐标代入式(4)，可求出待求点的残余高程异常ζresi，其与重力场模型高程异常ζGMi之和即为待求点的高程异常ζi。最后，将ζi代入式(1)可得到正常高hi。

1.3.2二次曲面法

二次曲面法模型如下：

(9)

与平面法模型类似，由大于6个点的已知点可以求解出未知数a0、a1、a2、a3、a4、a5，拟合得到二次曲面，进而求解出残余高程异常ζresi，再恢复重力场模型高程异常ζGMi，最后得到正常高hi。

1.3.3加权平均法

加权平均法模型待求点的残余高程异常ζres可由周围已知点计算，公式如下：

(10)

式中：n为设定半径区域内的已知点个数；Pi为各待求点的权。求出ζres后，叠加待求点的重力高程异常项，即可求出待求点的高程异常值。

权P是与待求点、已知点的距离有关的值，具体关系如下：

(11)

式中：S为拟合半径，是待求点与已知点的水平距离；e为定权常数，是任意常数，根据移动曲面模型的特点，本研究的定权常数e取0.01；n为距离定权次方，一般取大于零的数。

另外，为防止S的距离过小导致P值无限大，规定当S

1.4 精度评定

采用最小二乘法中的单位权中误差作为高程拟合精度分析的标准，如式(12)所示：

(12)

式中：v表示已知高程异常(H-h)与拟合所得高程异常(ζGM+ζres)的差值；n为检查点个数。

2 实例与精度评定

2.1 实验区情况

漳州市位于福建省最南部，东临台湾海峡，东北部与厦门市和泉州市接壤，西北部与龙岩市相连，西南部与广东省梅州市和潮州市毗邻。漳州市主要由山地、丘陵、平原组成，地势从西北向东南倾斜。

图1 已知点和检查点分布图 Fig.1 Distribution map of given points and check points

图2 实验区高程异常晕渲图

利用漳州市C级GNSS控制网303个控制点作为实验数据，所有控制点均包含CGCS2000坐标系下平面坐标、大地高和1985国家高程基准下的正常高成果，其中274个控制点作为已知点、29个控制点作为检查点，已知点和检查点的分布如图1所示。根据303个控制点数据求得的高程异常，绘制漳州市高程异常晕渲图，如图2所示。

2.2 高程异常长波项计算

在高程异常长波项计算过程中，将303个已知点用Alltrans EGM2008 Calculator软件计算长波项高程异常ζGM，表1列举了部分已知GNSS水准点的重力场模型ζGM、大地高H、正常高h等信息。残余高程异常晕渲图见图3，残余高程异常值最小为0.092 m，最大为1.08 m。

表1 部分已知GNSS水准点信息Tab.1 Partially given GNSS bench marking point information m

图3 残余高程异常晕渲图Fig.3 Dizziness map of residual elevation abnormal

2.3 残余高程异常拟合

依次对274个已知点用移动曲面模型法中的平面法、二次曲面法和加权平均法进行残余高程异常值拟合，其中待求点视为已知点的距离d取100 m，定权常数e取0.01，距离定权次方n依次取1、2、3，移动曲面拟合搜索半径依次取25 km、30 km、35 km。

2.4 检查点精度评价

为检验上述3种方法及不同的拟合参数求出的待定点高程异常值的精度，需要对待定点拟合高程异常值进行精度评定。本研究对前面选取的29个检查点进行检验计算，具体实测高程异常值与拟合高程异常值较差的最小值(取较差绝对值)、最大值(取较差绝对值)和中误差如表2至表4所示。

表2 检查点移动曲面之平面法精度统计Tab.2 Accuracy statistics of plane method of check point moving surface

表3 检查点移动曲面之二次曲面法精度统计Tab.3 Accuracy statistics of quadric surface method for checkpoint moving surface

表4 检查点移动曲面之加权平均法精度统计Tab.4 Accuracy statistics of weighted average method of checkpoint moving surface

从表2至表4可以看出，29个检查点的实测高程异常值与拟合高程异常值较差，采用二次曲面法，距离定权次方n=2和n=3时效果较好，最大值均小于10 cm，中误差均小于4 cm，完全满足CJJ/T 8—2011《城市测量规范》[13]中图根点的精度要求。

3 结语

采用移去-恢复法的高程拟合思想，用EGM2008模型计算高程异常长波分量，利用在地形起伏较大且控制面积较大的区域效果较好的移动曲面模型进行残余高程异常值的拟合。实践表明，该模型是一种有效的GNSS高程拟合方法，特别是利用移动曲面之二次曲面法，根据实验区内的检查点，实测高程异常值与拟合高程异常值较差最大值小于10 cm、中误差小于4 cm。