本篇地球物理論文研究地球重力場模型，介紹了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函數模型的似大地水準面擬合中有更好的適用性，并在實際生產中具有較好的實際推廣意義。