地球物理论文地球重力场模型

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

　　推荐期刊：《地球物理学报》创刊于1948年，是中国地球物理学会、中国科学院地质与地球物理研究所联合主办的有关地球物理科学的综合性学术刊物。主要刊载固体地球物理、应用地球物理、地磁和空间物理、大气和海洋地球物理，以及与地球物理密切相关的交叉学科研究成果的高质量论文。作者和读者对象主要为从事地球物理学、地球科学及其他相关学科的国内外科技工作者和大专院校师生。

　　关键词：EGM2008;Griddata函数;移去-恢复法;似大地水准面精化

　　由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水准数据确定局部似大地水准面的方法，并通过算例对其应用的可行性进行分析。

　　1EGM2008全球重力场模型简介

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

　　式中：(r，θ，λ)为以地球质心为坐标原点，Z轴与地球自转轴重合的坐标系中的球坐标(其中r为GPS 水准点的地心向径，θ为地心余纬，λ为地心经度)，γ-为该点的正常重力值，P-mn(cos θ)为完全正常化Legendre函数，C-mn，S-mn为完全正常化扰动位系数，其中偶次带谐系数代表实际引力位与正常引力位的系数差，m为实际地球引力位减去正常引力位的系数，α=6.378 136 46×106 m为参考椭球长半轴，nmax为系列地球重力场模型展开的球谐函数最高阶数，fM=3.986 004 415×1014 m3/s2为地心引力常数。

　　2EGM2008-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所示。

　　本实例将基于“移去-恢复法”模型分别采用griddata函数提供了4种插值方法对该测区不同数量的学习集样本进行插值比较。方法1：四格点样条插值 (v4);方法2：最近点插值(nearest);方法3：线性插值(linear);方法4：3次多项式插值(cubic)。在用 Griddata函数的4种方法进行数据内插时，线性插值和3次多项式插值法有时会因为计算舍入很难确定是否靠近临近边界而会出现“凸包”点问题，而导致结果输出为“NaN”的情况，本文采用利用四格点样条插值法或最近点插值法对“凸包”点进行代替。最后根据参与拟合计算已知点的高程异常残差值△ζi与拟合高程异常残差值△ζ′i，用vi=△ζ′i-△ζi求得拟合残差vi，得到不同学习样本数情况下的残差图，如图2~图4所示。

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

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

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

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

　　4总结

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