局部連續場下重力輔助導航模型構建算法

文超斌，王躍鋼田 琦郭志斌左朝陽楊家勝

（1.第二炮兵工程大學 304教研室，西安 710025；2.中國人民解放軍96124部隊，通化 134000）

局部連續場下重力輔助導航模型構建算法

文超斌1，2，王躍鋼1，田 琦1，郭志斌1，左朝陽1，楊家勝1

（1.第二炮兵工程大學 304教研室，西安 710025；2.中國人民解放軍96124部隊，通化 134000）

為了研究高精度的重力輔助導航模型以克服傳統模型的局限，必須建立精度高且具有良好解析性質的局部重力異常場解析模型，同時考慮模型誤差方程中的重力精準補償問題。針對二維高斯樣條函數逼近局部重力異常場中的局部支撐參數選擇問題，通過對所涉及的系數矩陣、解誤差、插值模型精度評估等問題進行分析，提出了一種新的最優局部支撐參數計算方法；基于此提出了一種高精度的基于高斯插值的重力輔助導航模型構建算法，該精準模型補償了重力擾動矢量、標準重力值誤差、厄特弗斯修正計算值對導航模型的影響。實驗結果表明利用新型重力輔助導航模型構建算法，可使輔助導航系統位置精度提高1倍左右，姿態、速度精度提高1～2倍，定位誤差保持在100～200 m。

重力輔助導航；高斯樣條插值；重力誤差補償；導航模型

重力匹配輔助慣性導航作為一種利用地球重力場特征信息結合匹配算法對對慣性導航系統位置信息進行初次估計，而后采用導航數學模型方程結合一定的濾波算法對慣性導航系統位置信息進行進一步精準修正的新技術，近年來得到了快速的發展。目前關于重力匹配算法的研究文獻較多，按照匹配算法設計原理可分為序列相關匹配方法和遞推濾波方法兩種，序列相關匹配方法主要包括最近等值線迭代算法（ICCP），相關極值分析算法兩大類[1-5]；獲取高精度、高分辨率重力異常圖是進行準確匹配的關鍵所在。高斯樣條函數具有局部支撐的良好特性，用其作為樣條基函數對計算區域重力異常進行二維整體逼近，能在保證重力異常局部特性不失真的前提下獲取計算區域重力異常的統一解析式等優點，近年來有關學者進行了相關應用研究。但針對該問題中涉及到的局部支撐參數計算問題分析并不多，局部支撐參數過大或過小都將導致很大的逼近誤差，太大或太小均失去插值意義，有必要對其進行尋優，研究相應的尋優算法[6-7]。

目前重力輔助導航系統狀態方程建模中通過地球正常重力模型得到的正常重力矢量進行重力補償，未考慮重力異常擾動矢量、標準重力值誤差（標準重力值誤差由測量位置維度誤差解算引起）和厄特弗斯修正計算值三類誤差影響，這些誤差都將引起帶來導航誤差，隨著慣性元件精度的不斷提高，慣性元件帶來的系統誤差與上述三類誤差帶來的系統誤差處在同一個量級，系統誤差不能被忽略，有必要采用精度更高的重力輔助導航系統狀態方程模型對正常重力矢量進行補償[8-12]。

所以，本文首先針對二維高斯樣條函數逼近局部重力異常場中的局部支撐參數選擇問題，通過對所涉及的插值系數矩陣可逆性進行研究給出了局部支撐參數的上限值，而后通過對插值系數矩陣解誤差研究，給出了矩陣條件數和局部支撐參數定性關系，最后結合插值模型精度評估理論，提出了一種新的高斯基插值函數最優局部支撐參數計算方法，接著局部連續重力異常場提出了一種考慮重力擾動矢量、標準重力值誤差、厄特弗斯修正計算值的精準重力輔助導航模型構建算法，最終為重力輔助慣性導航理論的高精度實踐化應用打基礎。

1 高斯基插值函數最優局部支撐參數計算

1.1 二維高斯樣條優化插值數學模型

式中，Δx、Δy為格網點分辨率[7]。這樣符合上述條件的xy平面二維高斯基函數可寫成矩陣形式：

插值條件：

將矩陣組合則有如下線性方程成立：

可見，高斯基插值函數存在兩個重要參數需要確定，第一：插值系數矩陣λ；第二：局部支撐參數a。如果X、Y均為非奇異矩陣，則式(2)有唯一解，將插值系數矩陣λ代入式(1)便得到該局部重力異常基準圖二維高斯基函數逼近解析式；局部支撐參數參數a決定了X、Y矩陣求逆的具體求逆誤差，應賦予局部支撐參數a值，使矩陣X、Y求逆誤差較小；當a過小時，由于局部支撐區間很小，使得插值曲線剛度小，容易出現鋸齒狀；當a過大時，由于局部支撐區間很大，插值曲線剛度大，插值曲線過于平滑甚至使插值解析式在已知點的插值誤差很大，所以a太大或太小均失去插值意義，因此有必要對參數a進行尋優，結合插值模型精度評估方法進行解算獲取參數a的最優值。

1.2 插值系數矩陣可逆性和求逆誤差

1.2.1 插值系數矩陣可逆性研究

設上述式(2)中高斯對稱矩陣階數K，如下示：

上述系數矩陣記為Ψ，考慮Ψ在什么條件下滿足非奇異，由數值計算方法理論知只要該矩陣是嚴格對

角占優即可。也就是對于第i行，：

嚴格對角占優：

成立，上述式(5)左邊可放大為：

所以只要確保式(7)成立則可以保證系數矩陣Ψ非奇異，對曲線的高斯插值總是可以保證的。可見a值最大限度為。

1.2.2 求逆矩陣誤差研究

對于矩陣Ψ，設為定向量空間nR上的向量范數，矩陣范數由式(8)給出：

可逆矩陣Ψ的條件數由式(9)給出：

由于計算機計算系數矩陣Ψ時，每個元素存在計算誤差，于是設擾動矩陣為E，則存在如下引理：

引理1：如果對任意滿足的擾動矩陣E和一個與E無關的正常數δ，矩陣Ψ的條件數由：給出，對于方程Ψx=b和解x和y，式(10)成立，則。

又由式(11)知：2/(1-β)＞1，所以0＜β＜1。于是當β單調增時，1-β單調減小，2(1-β)單調增，單調增。所以當時，此時對應的值就是方程的解x和y的最小的相對誤差上界。于是有以下定理結論：

定理1：如果對任意滿足的擾動矩陣E，矩陣Ψ的條件數由：給出，方程的解x和y，于是關于矩陣Ψ、E、x、y存在如下不等式：

定理1說明在誤差估計式(13)中，用求逆條件數可以得到最小的相對誤差上界。該等式給出了矩陣擾動相對誤差和由擾動引起的矩陣方程解相對誤差值間關系。

圖1 Ψ矩陣的條件數分析圖Fig.1 The condition number analysis of Ψ

1.3 插值模型精度評估

對局部支撐參數a進行尋優，為了保證插值函數與已知數據點保形性良好，同時考慮更多的函數非線性變化因素，研究利用分段二次拋物線代替所求插值函數對應的準確值，尋求分段二次拋物線與插值函數距離最小時對應的局部支撐參數amin，同時將二維局部重力異常基準圖問題解耦為x方向一維高斯基函數和y方向一維高斯基函數分別構造進行尋優建模。求解其對應的x方向最優局部支撐參數axmin，y方向最優局部支撐參數aymin，由此選擇可使插值函數與已知點具有最好的保形性，差值誤差較小。

具體建模方法為：對n行格網點，每行分別進行上述二維高斯插值局部支撐參數的尋優建模，而后利用所有尋優建模準則函數之和作為二維高斯樣條函數插值x方向的尋優模型，局部支撐優化函數(Local Support Optimize Function)LSOF尋優模型如式(15)示；對m列格網點，同x方向的尋優處理方法，LSOF如式(17)示，于是將局部支撐參數a限定在[0.2，2.5]區間內，分別按照式(15)(17)尋優準則進行計算，即可得到最優化局部支撐參數。

式中，ai，j、bi，j、ci，j為方程組(16)的解：

其中，ai，j、bi，j、ci，j為方程組(18)的解，

2 精準重力輔助導航模型構建算法

2.1 速度誤差精準模型方程

實際系統中總存在各種誤差，所以實際系統中考慮各種誤差情況下，實際速度計算值應由下述方程確定：

用式(20)減去式(19)，略去二階小量可得線性近似的捷聯慣導誤差方程為：

2.2 狀態方程

除速度誤差方程外導航系統其他誤差方程利用傳統思路方法參考相關文獻建立[13]，于是可以定義SINS的狀態變量如下，共15維狀態矢量：

于是得到局部連續重力輔助導航精準模型方程:

式中，矩陣wSINS為6維系統干擾白噪聲，矩陣GSINS是15×6噪聲分配矩陣，矩陣FSINS是15×15狀態矩陣。

2.3 量測方程

本文系統量測值取為重力儀量測重力與慣導指示位置對應重力圖上重力之差作為觀測方程。設定此時系統的量測方程Zg可由式(24)右邊計算得出：

由于繪制成的重力圖值，實質上等于實際重力測量值經過厄特弗斯修正而后減去標準重力值。所以對于實際位置點(L，)λ有：

代入方程(24)經線性化處理后，可表示為方程(26)所示系統的量測方程為：

VGravity為重力傳感器、重力異常圖制作噪聲和線性化誤差總和，在實際計算中，重力異常梯度，可以采用九點擬合、全平面擬合等方法計算。

所以，根據式(32)確定其測量矩陣為

式中，Zgline為式(24)中計算得到的重力測量誤差值，為1×15維觀測矩陣。為重力傳感器的噪聲、重力異常圖制作噪聲和線形化誤差總和。

3 仿真試驗與分析

本文采用第二炮兵工程大學慣性導航實驗室基于Matlab軟件自行開發研制的慣性導航軌跡生成軟件生成仿真需要的真實航跡和慣導指示航跡。重力異常數據使用渤海地區經度范圍 118°～119°，緯度范圍 37°～38°的真實數據，該重力異常場的三維圖如圖2所示。重力圖分辨率為 1′×1′，重力儀的實時量測數據是利用真實航跡在重力數據庫中的采樣值疊加0.4 mGal的隨機量測噪聲獲取，陀螺常值漂移取 0.01°/h，隨機噪聲漂移為0.01/h，加速度計隨機常值零偏1×10-4g，隨機噪聲漂移為 1×10-5g；初始狀態：速度 10 m/s、經度118.02°、緯度37.24°、慣導測量經維度誤差0.02°、對準水平姿態誤差角2′、對準航向誤差角10′。運載體保持俯仰、橫滾角 10°、航向角 30°的姿態，0.005 m/s2加速度每隔10 min加速一次，3 min采樣匹配濾波一次，行進了70 min時間，得到24個測量點[12]。捷聯慣導的情況下進行的仿真，對比分析了ICCP算法初次估計解算，局部連續場情況下傳統重力輔助導航模型和本文敘述的基于高斯插值的精準重力輔助導航模型導航解算法導航效果（分別稱其為tradition-ICCP和gausscompensation-ICCP算法，簡記為TRAI和GACOI算法）。

圖2 重力異常場三維圖Fig.2 3D map of gravity anomaly field

圖3(a)～圖3(b)為TRAI和本文GACOI新算法（最優局部支撐參數）導航的經、緯度位置誤差，航向、俯仰、橫滾姿態誤差角，東、北、天向速度誤差分布圖。如圖3中所示：TRAI算法導航位置誤差在 125～330 m、平均 210 m，水平姿態誤差在0.2′～0.75′、平均 0.50′，航向姿態誤差在 4′～9′、平均誤差7.5′，速度誤差在0.027～0.085 m/s、平均誤差0.045 m/s；本文 GACOI新算法導航能夠保持在100～200 m、平均110 m的位置誤差，0.01′～0.25′、平均 0.18′的水平姿態誤差，2′～6′、平均 2.5′的航向姿態誤差，0.002～0.05 m/s、平均0.016 m/s的速度誤差。可見，本文 GACOI導航新算法將傳統重力輔助慣性導航算法的位置導航精度提高了1倍左右，姿態角導航精度、速度導航精度提高了1～2倍。可見新算法引入了先進的局部支撐模型參數確定方法，同時建了考慮重力擾動矢量、標準重力值誤差、厄特弗斯修正計算值的精準重力輔助導航模型進行導航解算提高了導航精度，達到了接近于高精度導航的水平。

圖3 導航仿真誤差分析圖Fig.3 The analysis of navigation simulation errors

4 結 論

本文從分析重力輔助導航所需的地球重力場數學建模開始，通過對高斯樣條優化插值所涉及到的系數矩陣可逆性進行研究，得出了局部支撐參數的最大限度值，緊接著通過分析與之相關的求逆誤差進一步印證了結論的準確性，給出了局部支撐參數的數值區間，接著對插值模型精度評估理論進行研究，設計了給定區間上的最優局部支撐參數求解方法，而后建了考慮重力擾動矢量、標準重力值誤差、厄特弗斯修正計算值的精準重力輔助導航模型。最終實現了重力輔助導航算法的精確導航定位計算。本文研究內容不但可以為重力輔助慣性導航算法的實踐化做一定的理論鋪墊，而且相關理論研究對推進數學高斯插值相關研究也有一定的意義，當慣性導航系統導航觀測模型為非線性情況下，如何在高斯樣條逼近模型的基礎上展開相應的濾波算法有待進一步研究。

References)

[1]Mcphail S.Autosub6000:A deep diving long range AUV[J].Journal of Bionic Engineering，2009(6):55-62.

[2]李姍姍，吳曉平，馬彪.水下重力異常相關極值匹配算法[J].測繪學報，2011，40(4)：464-470.LI Shan-shan，WU Xiao-ping，MA Biao.Correlative extremum matching algorithm using underwater gravity anomalies[J].Acta Geodaetica et Cartographica Sinica，2011，40(4):464-470.

[3]袁贛南，張紅偉，袁克非，吳簡彤.基于重力梯度輔助定位的概率神經網絡改進方法[J].中國慣性技術學報，2013，21(3)：369-374.YUAN Gan-nan，ZHANG Hong-wei，YUAN Ke-fei，WU Jian-tong.Improved probabilistic neural network method based on gravity gradient aided location[J].Journal of Chinese Inertial Technology，2013，21(3):369-374.

[4]Jiang F，Wu Y M，Zhang Z S，et al.Combinational seabed terrain matching algorithm basing on probability data associate filtering and iterative closest contour point[C]//2009 Second International Conference on Intelligent Computation Technology and Automation.IEEE Computer Society，Piscataway，2009:245-249.

[5]Yuan G N，Zhang H W，Yuan K F.A Combinational underwater aided navigation algorithm based on TERCOM/ICCP and Kalman filter[J].IEEE Journal of Sciences and Optimization，2011，23(1):952-955.

[6]Wang Z G，Bian S F，Xiao S H.A local geopotential model for implementation of underwater passive navigation[J].Progress in Natural Science，2008，18(9):1139-1145.

[7]童余德，邊少鋒，將東方.一種新的基于局部重力圖逼近的組合匹配算法[J].地球物理學報，2012，55(9)：2917-2925.TONG Yu-de，BIAN Shao-feng，JIANG Dong-fang.A new integrated gravity matching algorithm based on approximated local gravity map[J].Chinese Journal of Geophysics，2012，55(9):2917-2925.

[8]Dransfield M，ZENG Y.Airbone Gravity gradiometry:Terrian corrections and elevation error[J].Geophysics，2009，74(5):137-142.

[9]李勝全，歐陽永忠，常國賓.慣性導航系統重力擾動量補償技術[J].中國慣性技術學報，2012，20(4)：410-414.LI Sheng-quan，OUYANG Yong-zhong，CHANG Guo-bin.Compensation technology of gravity disturbance vector in inertial navigation system[J].Journal of Chinese Inertial Technology，2012，20(4):410-414.

[10]DeGregoria A.Gravity gradiometry and map matching an aid to aircraft inertial navigation systems[D].Ohio:Air Force Institute of Technology，2010.

[11]Senobari M.New results in airborne vector gravimetry using strapdown INS/DGPS[J].Journal of Geodesy，2010，84:277-291.

[12]Neumeyer J，Sch?ferb U，Kremerc J，et al.Derivation of gravity anomalies from airborne gravimeter and IMU recordings-Validation with regional analytic models using ground and satellite gravity data[J].Journal of Geodynamics，2009，47:191-200.

Constructive model algorithm for gravity aided navigation in local anomaly field

WEN Chao-bin1，2，WANG Yue-gang1，TIAN Qi1，GUO Zhi-bin1，ZUO Zhao-yang1，YANG Jia-sheng1
(1.304 Unit，The Second Artillery Engineering University，Xi’an 710025，China;2.Unit 96124 of the Chinese people’s liberation army，Tonghua 134000，China)

In studying high-precision gravity aided navigation，it is indispensable to build a local continuous gravity anomaly field model with high precision and good analytic property，and accurately compensate the gravity in error equation.To solve the choice problem of local support optimal parameters in local continuous gravity anomaly field model with 2D gauss spline interpolation，an local support optimal calculation method is proposed by studying the related problems such as，coefficient of matrix，solution errors，and interpolation model evaluation，etc.Then a high-precision construct algorithm based on the gravity aided navigation model is put forward，which has compensated the gravity disturbance vector，standard gravity value error and eotvos effect.The experiment results show that the positioning accuracy is increased by about 1 times，the attitude and the velocity are both increased by 1～2 times.Its latitude and longitude errors are kept within 100～300 m.

gravity aided navigation;gauss spline interpolation;gravity error compensation;navigation model

U666.1

A

1005-6734(2014)03-0368-06

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

2013-12-25；

2014-04-11

國防預研（103030203）

文超斌（1986—），男，博士研究生，從事導航、制導與控制研究。E-mail：weijing123@126.com

聯 系 人：王躍鋼（1958—），男，教授，博士生導師。E-mail：wangyueg@163.com