基于EGM2008地球重力场模型-Griddata函数组合模型的似大地水准面精化探讨

付来迎，周世健，许亚男，王奉伟，周 青

(1.东华理工大学测绘工程学院，330013，南昌；2.江西省科学院，330096，南昌)

基于EGM2008地球重力场模型-Griddata函数组合模型的似大地水准面精化探讨

付来迎1，周世健2，许亚男1，王奉伟1，周 青1

(1.东华理工大学测绘工程学院，330013，南昌；2.江西省科学院，330096，南昌)

介绍了EGM2008地球重力场模型以及MATLAB提供的Griddata函数，结合工程实例，探讨基于EGM2008地球重力场模型-Griddata函数组合模型对似大地水准面进行精化，在Griddata函数提供的4种不同插值方法下进行分析比较，并根据4种不同插值的特点进行改进，最后得出基于EGM2008地球重力场模型-Griddata函数法中利用四格点样条插值法替代3次多项式插值法中凸包点法的拟合效果最好，具有很好的实际推广意义。

EGM2008；Griddata函数；移去-恢复法；似大地水准面精化

0 引言

由GPS测量定位得到的基线向量，经平差后可得到高精度的大地高程。若网中有一点或多点具有精确的WGS-84大地坐标系的大地高程，则在GPS网平差后，可得各个GPS点的WGS-84大地坐标系的大地高程。GPS相对定位高程方面的相对精度一般可达(2×3)×10-6，在绝对精度方面，对于10 km以下的基线边长，可达几个厘米，如果在观测和计算时采用一些消除误差的措施，其精度将优于±1 cm。但在实际应用中，地面点一般采用正常高程系统[1]。因此，为了尽可能减少外业工作量、广泛应用新技术，充分发挥GPS测量方便、快捷、精度高、成本低等优点，就必须知道WGS-84参考椭球面与似大地水准面之间的差距，即高程异常。知道了每个GPS点的高程异常ζ，就能通过式(1)，计算出相应的正常高，实现用GPS测量代替水准测量。

H大地高=H正常高+ζ

(1)

确定高程异常方法大致可分为几何曲面拟合法和重力法两大类。几何曲面拟合法就是根据区域内若干个公共点上的高程异常值，构造某一种曲面逼近似大地水准面，由于所构造的曲面不同，计算方法也不同，其主要方法有：平面拟合法、曲面拟合法、多面函数拟合法、样条函数法等[2]。重力法是指利用局部或者全球重力数据求解高程异常，将GPS大地高转换为正常高。重力法充分利用了地球物理场，但由于重力数据缺乏和重力场模型精度对于工程应用而言较低，所以工程建设一般都不采用重力法直接进行GPS高程转换[3]。在区域似大地水准面精化中，最严密、有效的方法是利用由地球重力场模型、地面重力数据和DEM数据获得的该区域剩余重力异常，采用移去-恢复法确定重力似大地水准面，再用GPS水准数据对重力似大地水准面进行拟合，求得与国家或地方高程系统定义一致的似大地水准面[4]。由于EGM2008模型的精度明显好于其他模型[5]，而且它的相对精度比绝对精度高[6]，所以将应用EGM2008模型计算的高程异常差和GPS水准数据确定局部似大地水准面的方法，并通过算例对其应用的可行性进行分析。

1 EGM2008全球重力场模型简介

EGM2008是美国国家地理空间情报局(National Geospatia1-Intelligence Agency)，简称NGA，经过多年的研究和总结，在以往构建地球重力场模型的经验和理论基础上，采用最先进的建模技术与算法，以PGM2007B(PGM2007A的变种模型)为参考模型，利用GRACE卫星采集的重力数据和全球5′×5′的重力异常数据，TOPEX卫星测高数据以及现势性好、分辨率高的地形数据，结合精度高，覆盖面广的地面重力数据完成的最新一代全球重力模型[7]。该模型的阶次完全至2 159(另外球谐系数的阶扩展至2 190次)，相当于模型的空间分辨率约为5′(约9 km)。地球重力场球谐函数模型EGM2008由一系列完全正常化的球谐函数系数组成，根据球谐函数连续叠加可计算扰动位T，利用Bruns公式可求得地面上任意点的高程异常[8]。

(2)

2 EGM2008-Griddata组合模型原理

由于EGM2008模型定义的全球高程基准与我国采用的国家高程基准不同，两者之间存在着一个系统差，因而，首先将高程异常可如式(3)可分成两部分

ζ=ζGM+△ζ

(3)

式中：ζGM为EGM2008地球重力场模型求得的高程异常；△ζ为实测高程异常与EGM2008地球重力场高程异常的残差。从而，当利用EGM2008模型计算高程异常进而推算正常高时，式(1)相应改为：

H大地高=H正常高+ζGM+△ζ

(4)

利用“移去-恢复”法，首先，将基于已有的m个GPS水准联测点，可通过式(1)的变形式：ζi=H大地高i-H正常高i,(i=1，2，3，…，k)，计算出这k个点的高程异常值ζi；再利用地球重力场模型EGM2008计算这些点的高程异常的近似值ζGMi，然后根据式(3)计算出这k个点的高程异常残差△ζi。第2，将求出的k个高程异常残差值△ζi作为已知数据，用MATLAB所提供的Griddata函数插值方法进行拟合计算，最后再内插出未知点的高程异常残差△ζ。最后，利用地球重力场模型EGM2008求出未联测水准的点上的高程异常近似值ζGM，再与拟合出的该点的高程异常残差值△ζ相加，得出最终的高程异常值ζ，再通过利用GPS测得的大地高数据由式(1)计算出所有待求点的正常高。

3 实例分析

MATLAB是MathWorks公司推出的一套高性能的数值计算和可视化软件，集数值分析、矩阵运算与图形显示于一体，可方便地应用于数学计算、数据分析和工程绘图，所采用的数值计算方法均采用公认的先进、可靠的算法，其程序均由世界一流专家编制并经高度优化[9]。MATLAB中的Griddata函数可以将位于同一空间坐标系下的散点插值为规则格网，提供了4种插值方法：线性插值(linear)、3次多项式插值(cubic)、最近点插值(nearest)、四格点样条插值(v4)，可以方便地实现结合邻近离散点分布特征的光滑曲面拟合。其中线性插值和最近点插值结果构成的曲面不光滑不连续，3次多项式插值和四格点样条插值结果构成的曲面较光滑[10]。

为对Griddata函数提供的4种插值方法进行合理的选择，本文以某地区16个GPS水准数据为例，该测区南北跨度约20 km，东西跨度约15 km，地形较为复杂，测区起伏较大，最大高差约为170 m之多，其点位分布情况如图1所示。

图1 点位分布情况

图2 4种方法4个学习样本拟合残差对照图

图3 4种方法6个学习样本拟合残差对照图

图4 4种方法8个学习样本拟合残差对照图

通过以上3个图对比可以看出：较于最近点插值法替代，四格点样条插值法替代线性插值法和3次多项式插值法中出现的“凸包”点的效果更好。同时，从以上3幅图可以看出“凸包”点替代后的3次多项式插值法插值残差值最小，替代后的线性插值法效果次之，四格点样条插值法效果稍差，最近点插值法残差波动性较大且残差值较大；并且4种方法都会随着学习样本数量的增加，插值的残差减小且波动性会降低。

通过4种方法不同数量学习样本得到的拟合成果还可以分别通过内符合精度和外符合精度进行评定，其中内符合精度可按式(5)计算求得，外符合精度可按式(6)求得

(5)

(6)

其中,vi为各点的拟合残差，m为学习集样本个数，n为总体样本数量。其结果如表1所示。

由表1可以看出，随着样本数量的增多，4种不同的插值方法插值效果都有所改善，最大偏差值都在减小，外符合精度不断提高；同时，在线性插值法与3次多项式插值法插值中利用四格点样条插值法进行凸包点替代的外符合精度高于最近点插值法进行替代的外符合精度。

表1 不同插值方法不同数量学习样本精度统计表

4 总结

通过工程实例采用基于EGM2008重力场模型的Griddata函数的4种不同方法进行插值比较，利用四格点样条插值法进行替代的线性插值法与3次多项式插值法有着较好的插值效果，基于插值效果与拟和曲面光滑性的综合考虑，利用四格点样条插值法进行凸包点替代的3次多项式插值法在采用基于EGM2008重力场模型的Griddata函数模型的似大地水准面拟合中有更好的适用性，并在实际生产中具有较好的实际推广意义。

[1] 张福荣,田倩.GPS测量技术与应用[M].成都:西南交通大学出版社,2013:177-178.

[2]张勤,李家权,等.GPS测量原理及应用[M].北京:科学出版社,2005:214-219.

[3]谷延超,范东明.顾及EGM2008和残差地形模型的GPS高程转换方法研究[J].测绘工程,2013,22(2):26-29.

[4]冯林刚,郅军义,宝因乌力吉.应用EGM2008模型和GPS/水准数据确定局部似大地水准面[J].测绘通报,2011(1):18-20.

[5]章传银,郭春喜,陈俊勇,等.EGM2008地球重力场模型在中国大陆适用性分析[J].测绘学报,2009,38(4):283-289.

[6]张兴福,刘成,刘红新.利用GPS冰准数据检核EGM2008重力场模型的精度[J].测绘通报,2009(2):7-9.

[7]Pavlis N K,Holmes S A,Kenyon S C,etal.An Earth Gravitational Model to Degree 2160:EGM2008[C].Vienna:Presented at the 2008 General Assembly of the European Geosciences Union,2008.

[8]游为,范东明,付淑娟,等.GPS高程转换的新方法研究[J].工程勘察,2009,37(3):60-62.

[9]求是科技.Matlab7.0从入门到精通[M].北京:人民邮电出版社,2006.

[10]吴守亮.基于Matlab的三维数字地形模拟及空间分析[J].测绘工程,2011,20(3):54-57.

DiscussionofRefiningtheGeoidBasedonEGM2008-GriddataFunctionCombinationModel

FU Laiying1,ZHOU Shijian2,XU Yanan1,WANG Fengwei1,ZHOU Qing1

(1.East China Institute of Technology.Faculty of Geomatics,330013,Nanchang,PRC;2.Jiangxi Academy of Science,330096,Nanchang,PRC)

The paper introduced the EGM2008 gravity field model and the Griddata function that was provided by MATLAB.At the same time,the paper combining with the engineering example presented a new method which was the EGM2008 gravity field model-Griddata function model to refine the geoid.Meanwhile,Analysis and comparison was made with four kinds of interpolation methods provided by Griddata function.At the last,we got a conclusion that the new method had a good fitting effect on the geoid refinement by analyzing the result of the date,and it had an important significance in the actual production.

EGM2008;Griddata function;remove-restore method;refine the geoid

2014-10-12;

2014-11-21

付来迎(1990-)，男，山东莱芜人，硕士研究生，主要研究方向为现代测量数据处理。

国家自然科学基金资助项目(41374007)；2014年江西省研究生创新专项资金项目(YC2014-S317)。

10.13990/j.issn1001-3679.2014.06.012

P223

A

1001-3679(2014)06-0794-04