基于EGM2008重力場球諧模型的水平重力擾動計算方法

鐵俊波，吳美平，蔡劭琨，張開東，練軍想

（國防科技大學 機電工程與自動化學院，長沙 410073）

基于EGM2008重力場球諧模型的水平重力擾動計算方法

鐵俊波1，吳美平1，蔡劭琨1，張開東1，練軍想1

（國防科技大學 機電工程與自動化學院，長沙 410073）

垂線偏差數據對于提高慣性導航精度具有重要意義，闡述了使用重力場球諧模型計算水平重力擾動的方法，使用EGM2008重力場球諧模型計算水平重力擾動，并將其與美國國家地理空間情報局提供的全球2.5′2.5′× 垂線偏差網格數據進行對比，驗證了所述水平重力擾動計算方法的正確性。將水平重力擾動計算結果用于航空、車載慣性導航數據離線處理，結果表明EGM2008模型計算的水平重力擾動可用于補償慣性導航；將水平重力擾動計算結果用于長航時船載慣性導航數據離線處理，結果表明慣性導航精度最大提升1.5 n mile，平均提升0.8 n mile。

慣性導航；重力擾動；重力場球諧模型；垂線偏差補償

為進一步提高慣性導航精度，補償水平重力擾動對慣性導航的影響是必須要解決的問題[4-7,15-18]，而獲取高精度的水平重力擾動數據就是首先要面對的問題。水平重力擾動數據的測量是困難的，使用天文方法可以靜態、單點測量，但效率很低，難以獲取區域數據。捷聯式重力矢量測量方法有望在未來成為高效的水平重力擾動測量手段，通過直接求差法[2-3]或姿態求差法[19-21]可以獲得水平重力擾動數據，但其對水平姿態測量精度要求非常高，為提高重力矢量測量精度需進一步研究姿態誤差分離、傳感器漂移補償等問題[3]。目前獲取大面積水平重力擾動的可行方法是利用全球重力場球諧模型計算水平重力擾動，近年來國內學者研究了使用重力場模型計算水平重力擾動用于提高慣性導航精度的方法。文獻[4-5]以美國國家大地數據測量局（National Geospatial-Intellgence Agency, NGA）提供的垂線偏差網格數據對慣性導航系統進行了補償，但使用網格數據需要插值及延拓處理。文獻[6-7]使用重力場球諧模型直接計算水平重力擾動，但其給出的水平重力擾動計算方法不完整。為此本文闡明了使用球諧模型計算重力擾動的方法，將計算結果與NGA官方提供的垂線偏差網格數據進行比較，并結合航空、車載及船載激光陀螺捷聯慣性導航數據對所計算的水平重力擾動數據能否提高慣性導航精度進行了驗證。

1 擾動重力位與水平重力擾動

其中：δgN、δgE為北向、東向重力擾動；θ為計算點的地心余緯度，又稱為極距（Polar Distance）；λ為計算點經度；ρ為計算點到參考橢球球心的距離。擾動重力位滿足拉普拉斯方程，可用球諧函數展開：

圖1 垂線偏差示意圖Fig.1 The definition of vertical deflection

由重力擾動位與水平重力擾動關系(1)(2)及重力擾動位的球諧函數表達式(3)，可以得到水平重力擾動計算公式[1,10-12]：

完全正常化勒讓德函數通過地推的方式進行求解，需要注意的是，計算以階n為順序遞推，即n為變量、m為常量，其遞推公式如下（根據n與m的取值關系，遞推時分為3種情況）[7,9]：

完全正常化勒讓德函數遞推初值為

由式(4)(5)看到，使用球諧模型計算東向重力擾動僅需要完全正常化勒讓德函數，但求解北向重力擾動還需要完全正常化勒讓德函數的導數。文獻[7]給出的求解完全正常化勒讓德函導數的公式是相對于cosθ的導數，而由公式(4)看到，計算北向水平重力擾動的導數應是相對于θ的導數，使用如下遞推公式計算[8]：

2 擾動重力位球諧函數計算水平重力擾動方法

完全正常化勒讓德函數及其導數的計算是使用球諧模型計算水平重力擾動的關鍵，由公式(6)～(11)可以看到：m=n的完全正常化勒讓德函數及其導數可以直接遞推計算，將其計算過程稱之為“等階次完全正常化勒讓德函數及其導數計算”；而m＜n的完全正常化勒讓德函數及其導數的計算有賴于等階次完全正常化勒讓德函數及其導數的計算結果，將此計算過程稱之為“非等階次完全正常化勒讓德函數及其導數計算”。水平重力擾動算法流程如圖2所示。

圖2 球諧模型計算水平重力擾動算法流程圖Fig.2 The calculation scheme of vertical deflection based on gravity harmonic model

第一步：參考橢球參數初始化

根據給定的參考橢球的4個基本參數，計算其它參數，4個基本參數包括：長半軸a、橢球的扁率f或2階帶狀球諧系數J2、地心引力常數GM及地球自轉角速率Ω。這4個基本不僅決定了參考橢球的形狀，還確定了與此參考橢球對應的正常重力場。NGA官方提供的球諧系數是在WGS-84橢球模型下計算得到，慣性導航一般也是以WGS-84為參考橢球，當使用其它橢球模型時，NGA提供的球諧系數模型需要進行一定處理，WGS-84橢球坐標系的4個基本參數取值如表1所示。

表1 WGS-84橢球模型基本參數Tab.1 Fundamental parameter of WGS-84

第二步：讀取球諧模型系數并處理

EGM2008官方網站[10]給出了無潮汐假設（Tide-Free）下EGM2008重力場球諧模型系數文件：系數文件的第1列為階n；第2列為次m；第3、4列為與完全正常化勒讓德函數對應的球諧系數該系數是在WGS-84坐標系下計算得到，當使用的參考橢球非WGS-84時，需要對其進行改正。水平重力擾動是與擾動重力位相關，NGA提供的EGM2008球諧系數是重力位球諧系數而非擾動重力位球諧系數，因此需要將正常重力場從其中扣去后得到所需的擾動重力位球諧系數。令正常重力位球諧系數中與cosmλ相關項為使用如下公式(14)將其從EGM2008球諧系數中扣除，得到擾動重力位的球諧系數

e為橢球的第一離心率，可由扁率f計算得到：

第三步：讀取計算點坐標及非勒讓德函數項計算

重力場擾動位球諧模型公式中的θ為地心余緯度，當輸入計算點坐標為地心緯度時可以直接用π/2減去地心緯度后得到地心余緯度θ。而慣性導航、衛星導航中一般使用大地緯度，因此需要將大地緯度轉換為地心緯度。擾動位球諧函數中的項與完全正常化勒讓德函數無關，可以單獨計算。需要注意的是，在計算時，NGA的官方程序中給項乘以了系數乘以此系數的目的是為了提高完全正常化勒讓德函數的數值計算精度[13]。

第四步：等階次完全正常化勒讓德函數及其導數計算

第五步：非等階次完全正常化締合勒讓德多項式及其導數計算

第六步：水平重力擾動計算將第二步得到的球諧系數，第三步得到的計算點地心余緯度、非勒讓德函數項，第四、五步得到的完全正常化勒讓德函數及其導數代入公式(4)(5)得到計算點處的水平重力擾動。

3 與NGA提供垂線偏差網格數據比較

將水平重力擾動計算結果與 NGA官方提供的全球垂線偏差網格數據比較，以驗證所述算法計算結果的正確性。NGA官方提供了全球分辨率為2.5′2.5′× 的垂線偏差網格數據，需要注意的是，該網格數據是橢球面上的垂線偏差網格數據，將其乘以正常重力γ后即可得到水平重力擾動，正常重力計算公式參見文獻[13]，水平重力擾動與垂線偏差關系為

取 NGA 提供的網格數據中范圍為 10°N～11°N（大地緯度），110°E～112°E 的網格垂線偏差數據，利用公式(12)(13)將其轉換為水平重力擾動。獲取其坐標后使用所述方法計算水平重力擾動，兩者差異如表 2、圖 3～4所示，存在差異的主要原因在于NGA給出的網格數據是該網格區域的平均值，而使用球諧模型計算的水平重力擾動是計算點處的單點值，因此存在微小差別。

表2 水平重力擾動計算結果差異統計Tab.2 Differences between grid datum and computation

圖3 水平重力擾動南北向分量差異Fig.3 Differences of north-south component

圖4 水平重力擾動東西向分量差異圖Fig.4 Differences of east-west component

4 實測數據驗證

計算水平重力擾動的目的是提高慣性導航精度，使用航空及車載試驗數據對所述水平重力擾動計算方法進行驗證。車載數據為湖南某地90型激光捷聯慣導系統車載試驗數據，試驗車向西北方向行駛，時間約4 h；航空數據為山東某沿海地區的90型激光捷聯慣導系統飛行試驗數據，載機沿東西向往返飛行，時間約4 h。根據采樣率為2 Hz的差分GPS給出的大地緯度、經度及高度，以所述算法計算載體出發點及運動路徑上的水平重力擾動，如圖5～6所示。兩次試驗路徑上的水平重力擾動統計信息見表3所示。

圖5 飛行試驗路徑上的水平重力擾動Fig.5 Gravity disturbance on the airborne trajectory

圖6 車載試驗路徑上的水平重力擾動Fig.6 Gravity disturbance on the terrestrial trajectory

表3 飛行、車載試驗 水平重力擾動統計值Tab.3 Statistical value of gravity disturbance

由水平重力擾動計算結果及正常重力構成重力矢量用于慣性導航，以差分GPS定位結果作為位置基準，得到未補償及補償后的慣性導航位置誤差，用未補償情況的位置誤差減去補償后的位置誤差，如圖7～8所示。

文獻[4]中定量分析了在4 h慣性導航條件下，在一般區域及異常較大區域，水平重力擾動對慣性導航精度的影響，飛行及車載試驗的慣性導航時間也恰好是4 h，因此文獻[4]中的仿真分析結論對此處的實際數據計算結果分析有較大參考意義。將表3統計結果與文獻[4]進行比較，飛行及車載試驗區域應屬于文獻[4]所定義的水平重力擾動一般區域，在水平重力擾動一般區域，水平重力擾動造成的慣性導航位置誤差在百米量級。如圖6～7所示，4 h飛行及車載試驗中，在使用計算的水平重力擾動數據對慣性導航進行補償后，慣性導航定位精度最大提高120 m及110 m，與文獻[4]中的仿真分析結論一致。為進一步驗證長航時情況下水平重力擾動數據對慣性導航精度的提升效果，使用南海某海域24 h船載激光捷聯慣導系統試驗數據進行驗證，使用GPS精密單點定位結果作為位置基準，使用位置基準信息及所述方法計算航路上的水平重力擾動數據，如圖9、表4所示。

圖7 飛行試驗 補償前后定位誤差比較Fig.7 Position error comparison of airborne test

圖8 車載試驗 補償前后定位誤差比較Fig.8 Position error comparison of terrestrial test

圖9 船載試驗路徑上的水平重力擾動Fig.9 Gravity disturbance on the shipborne trajectory

表4 船載試驗中水平重力擾動統計值Tab.4 Statistical value of gravity disturbance

圖10 船載試驗 北向誤差對比Fig.10 Latitude error comparison of shipborne test

圖11 船載試驗 東向誤差對比Fig.11 Longitude error comparison of shipborne test

使用水平重力擾動數據構建重力矢量用于慣性導航，將補償后位置誤差與未補償時的位置誤差對比，緯度誤差、經度誤差及位置誤差對比如圖10～12所示。由誤差對比圖可以看出，在使用了水平重力擾動數據補償后，慣性導航精度高于使用正常重力模型情況，因此認為所計算的水平重力擾動數據可用于提高長航時慣性導航精度。對比表3、表4可以看到，相比于航空、車載試驗區域，船載試驗區域的水平重力擾動變化更大，因此使用水平重力擾動數據補償后，慣性導航系統精度提升更加明顯。

為定量評估垂線偏差數據對長航時船載慣性導航系統精度補償效果，將圖12中補償前后在時間上對應的最大位置誤差相減作為評估的依據，慣性導航精度最大提升1.5 n mile，平均提升0.8 n mile。

圖12 船載試驗位置誤差對比Fig.12 Position error comparison of shipborne test

5 結 論

為解決高精度慣性導航對水平重力擾動數據需求，闡明了使用 EGM2008重力場球諧模型計算水平重力擾動的方法，水平重力擾動計算結果與 NGA官方提供的網格數據基本一致。

將水平重力擾動計算結果用于航空、車載及船載激光捷聯慣性導航系統數據離線處理，計算結果表明，使用所述的水平重力擾動數據計算方法得到的水平重力擾動數據可用于提升慣性導航精度。

（References）：

[1]Hofmann-Wellenhof B, Moritz H.Physical geodesy[M].Austria: Springer-Verlag Wien, 2006.

[2]張開東.基于 SINS/DGPS的航空重力測量方法研究[D].長沙: 國防科學技術大學, 2007.Zhang K D.Research on the methods of airborne gravimetry based on SINS/DGPS[D].Changsha: National University of Defense Technology, 2007.

[3]蔡劭琨.航空重力矢量測量及誤差分離方法研究[D].長沙: 國防科學技術大學, 2014.Cai S K.The research about airborne vector gravimeter and methods of errors separation[D].Changsha: National University of Defense Technology, 2014.

[4]趙忠, 王鵬.高精度慣性導航系統垂線偏差影響與補償[J].中國慣性技術學報, 2013, 23(6): 701-705.Zhao Z, Wang P.Analysis and compensation of vertical deflection effect on high accuracy inertial navigation system[J].Journal of Chinese Inertial Technology, 2013,23(6): 701-705.

[5]周瀟, 楊功流, 蔡慶中.基于小波神經網絡的高精度慣導重力擾動補償方法[J].中國慣性技術學報, 2016,24(5): 571-576.Zhou X, Yang G L, Cai Q Z.Compensation on gravity disturbance for high-precision INS based on wavelet neural network[J].Journal of Chinese Inertial Technology,2016, 24(5): 571-576.

[6]何昆鵬, 王曉雪.重力擾動對高精度慣性導航初始對準的影響與補償[J].導航與控制, 2015, 14(5): 84-88.He K P, Wang X X.The effect of gravity disturbance on high-precision inertial navigation system[J].Navigation and Control, 2015, 14(5): 84-88.

[7]Wang J, Yang G L, Li X Y, et al.Application of the spherical harmonic gravity model in high precision inertial navigation systems[J].Measurement Science and Technology, 2016, 27: 095103(1-10).

[8]魏子卿.完全正常化締合勒讓德函數及其導數與積分的遞推關系[J].武漢大學學報(信息科學版), 2016,41(1): 27-36.Wei Z Q.Recurrence relations for fully normalized associated legendre functions and their derivatives and integrals[J].Geomatics and Information Science of Wuhan University, 2016, 41(1): 27-36.

[9]Jekeli C, Lee J K, Kwon J H.On the computation and approximation of ultra-high-degree spherical harmonic series[J].Journal of Geodesy, 2007, 81(9): 603-615.

[10]http://earth-inf.nga.mil/GandG/wgs84/gravitymod/egm2008.

[11]Rapp R H.A FORTRAN program for the computation of gravimetric quantities from high degree spherical harmonic expansions[R].Department of Geodetic Science and Surveying, The Ohio State University, 1982.

[12]Colombo O L.Numerical methods for harmonic analysis on the sphere[R].Department of Geodetic Science and Surveying, The Ohio State University, 1981.

[13]Singh A.On numerical evaluation of normalized associated legendre functions[R].Department of Geodetic Science and Surveying, The Ohio State University, 1982.

[14]Titterton D H.Strapdown inertial navigation technology[M].2nd Edition.American Institute of Aeronautics and Astronautics, Inc., 2004.

[15]Wu R N, Wu Q P, Han F T, et al.Gravity compensation using EGM2008 for high-precision long-term inertial navigation systems[J].Sensors, 2016, 16: 2177.

[16]Zhou X, Yang G L, Cai Q Z, et al.A novel gravity compensation method for high precision free-INS based on extreme learning machine[J].Sensors, 2016, 16: 2019.

[17]Zhou, X Yang G L, Wang J, et al.An improved gravity compensation method for high-precision free-INS based on MEC-BP-AdaBoost[J].Measurement Science and Technology, 2016, 27: 125007.

[18]Fang J C, Chen L Z, Yao J F.An accurate gravity compensation method for high-precision airborne POS[J].IEEE Transactions on Geoscience and Remote Sensing,2014, 52(8): 4564-4573.

[19]Dai D, Wang X, Zhan D, et al.An improved method for dynamic measurement of deflections of the vertical based on the maintenance of attitude reference[J].Sensors, 2014,14: 16322-16342.

[20]Dai D, Wang X, Zhan D, et al.Dynamic measurement of high-frequency deflections of the vertical based on the observation of INS/GNSS integration attitude error[J].J.Appl.Geophys, 2015, 119: 89-98.

[21]Zhu J, Zhou Z B, Li Y, et al.Further development of the attitude difference method for estimating deflections of the vertical in real time[J].Measurement Science and Technology, 2016, 27: 075004.

Gravity disturbance calculation method based on Earth Gravitational Model 2008

TIE Jun-bo1, WU Mei-ping1, CAI Shao-kun1, ZHANG Kai-dong1, LIAN Jun-xiang1
(College of Mechatronics Engineering and Automation, National University of Defense Technology,Changsha 410073, China)

Vertical deflection can be used to improve the performance of inertial navigation system.The calculation method of vertical deflection based on harmonic model is described.EGM2008 is used to calculate vertical deflection.Comparison between the calculation result and 2.5′ 2.5′× DOV grid data from NGA verifies the validity of the described calculation method.Through the airborne and terrestrial strapdown inertial navigation calculation, it is shown that the gravity disturbance from EGM2008 can be used to improve the position accuracy of inertial navigation system, and the maximum and average location accuracies of inertial navigation system are increased by 1.5 n mile and 0.8 n mile, respectively, in the longtime shipborne test.

inertial navigation; gravity disturbance; gravity harmonics model; compensation for vertical deflection

U666.1

A

1005-6734(2017)05-0624-06

10.13695/j.cnki.12-1222/o3.2017.05.012

2017-05-27；

2017-09-11

國家自然科學基金項目（61603401）；國土資源部航空地球物理與遙感地質重點實驗室 青年創新基金課題（2016YFL06）

鐵俊波（1989—），男，博士研究生，從事慣性導航系統研究。E-mail: tiejunbo11@nudt.edu.cn

聯 系 人：吳美平（1970—），男，教授，博士生導師。E-mail: meipingwu@263.net