基于二維高斯樣條函數的水下重力被動定位?

王志剛 顧雪峰

（1.91550部隊 大連 116023）（2.海軍工程大學兵器工程系 武漢 430033）（3.海軍工程大學兵器科研部 武漢 430033）

1 引言

目前，對于水下潛器的被動定位是導航、定位領域研究的一個熱點和難點問題。水下潛器的慣導系統（Inertial Navigation System，INS）誤差隨時間積累，因此必須要通過其它導航方式實時或定期修正INS。出于隱蔽性要求，水下潛器又很難利用衛星或無線電信息，此時利用水下地理特征信息的輔助導航定位就是很好的選擇。常用的輔助導航定位系統有地形、重力輔助導航系統［1～2］等，目前這些輔助導航系統都是以各種匹配算法為核心來獲得潛器的最佳匹配位置（真實位置），以實現對慣導系統的修正。傳統的匹配算法主要分為單點迭代和序列迭代兩類，分別以SITAN［3～4］（Sandia Inertial Terrain－Aided Navigation）和ICCP［5～6］（Interval Closest Contour Point）算法為典型代表。這些匹配算法分別在實時性和穩定性方面有其優點，但是為了保證良好的匹配效果對潛器的運動規律和初始誤差環境有很多嚴格的限制［7］，因而也限制了其應用范圍。針對這一問題，本文提出一種直接利用觀測重力進行水下被動定位的新模式，應用該方法的前提就是確定重力異常與潛器位置（大地經、緯度坐標）的精確解析表達式，之后與航空領域的目標被動定位原理相似，由于觀測重力中包含了目標的位置坐標信息，直接通過非線性濾波算法，就可對潛器的航行位置進行最優估計。

2 二維高斯樣條函數逼近

直接利用地球重力進行定位，在很大程度上取決于觀測重力與目標所在位置關系的確定，也就是地球重力場模型的確定。目前常用的各類全球重力場模型［8］無論在階次和精度方面都已經達到了很高的水平，它本質上屬于球坐標系下的調和分析，因而應用過程中需要計算大量的Legendre系數，使得模型計算量巨大，而這對于定位的實時性是一大考驗，另外現有的大部分全球重力場模型非線性比較嚴重，也不太利于濾波估計時的線性化處理。

近年來，高斯樣條函數在醫學影像、數字地形構建、計算機圖形學及地球物理等領域有較為廣泛的應用，文獻［8］討論了一維高斯樣條在協方差函數代數確定中的應用，文獻［9］采用高斯函數作為樣條基函數對計算區域重力異常進行二維整體逼近，文獻［10］討論了高斯樣條函數逼近局部重力異常的相關影響因素。采用高斯樣條函數逼近局部重力異常運算簡單且最終解析式是統一的，文章在該樣條逼近函數的基礎上建立了重力量測方程并結合具體的目標運動狀態方程，成功運用擴展Kalman濾波實現了對潛器位置的最優估計。

2.1 二維高斯樣條函數逼近數學原理

即有

不妨設

則上式可以簡化為矩陣形式：

則根據插值條件｛L（xi，yj）＝zi，j｜i＝1，…，m；j＝1，…，n｝，有如下線性方程：

易知X、Y 均為非奇異矩陣，即式（5）有唯一解，解方程即得到系數矩陣C，將系數矩陣代入式（2）便得到該局部重力異常基準圖二維高斯樣條函數逼近解析式。

由式（5）可以看出系數矩陣C由X、Y與Z唯一確定，而X、Y由矩陣的階數和ax，ay決定。文獻［9］曾對求逆計算誤差與和計算階數的關系進行過詳盡的探討。經討論可知：系數矩陣的求解過程中如X、Y矩陣階數過高且ax、ay過小或階數過低且ax、ay過大時，X、Y矩陣的求逆存在較大誤差。而當X、Y矩陣階數小于5或ax、ay值小于5時，矩陣求逆運算誤差較小。考慮到求逆計算誤差以及實際應用中用到的格網區域遠大于5×5，我們即可采用文獻［10］給出的斐波那契數列方法對ax∈［1，s］，ay∈［1，s］區間內進行尋優。

2.2 二維高斯樣條函數逼近精度分析

取范圍i為117°E～121°E，j為21°N～25°N，分辨率為2′×2′的衛星測高反演重力異常數據作為已知重力異常基準圖。文章在該已知2′×2′重力異常基準圖的基礎上，取4′×4′格網點處重力異常數據作為基礎數據點，進行二維高斯樣條函數逼近，獲得該范圍局部重力異常基準圖解析表達式，并據此計算分辨率為2′×2′的重力異常逼近值，然后與已知的2′×2′格網處重力異常值進行比對分析。

最終的仿真結果如圖1所示，其中圖1（a）為已知2′×2′重力異常基準圖，圖1（b）為由4′×4′重力異常基準圖逼近的2′×2′重力異常基準圖，由此可以看出二維高斯樣條函數逼近的局部重力異常基準圖解析式能較好地描述這一范圍的重力異常基準圖。

圖1（c）、圖1（d）分別為逼近誤差圖及相應的誤差統計圖，表1為誤差統計表，由圖1、表1可以得出，試算區域4′×4′格網點處二維高斯逼近計算值絕對誤差均值為0.19mGal且99.5%計算點絕對誤差小于1mGal，逼近精度可以滿足導航定位要求。鑒于重力二維高斯樣條逼近函數的連續形式、與經緯度坐標具有明確的解析關系以及自身的高精度，使得這一逼近函數十分適合用于重力被動定位時所需的量測方程。

圖1 局部重力異常基準圖二維高斯逼近仿真實驗

表1 計算值和真值比較（單位：mGal）

3 水下被動定位模型的建立

水下被動定位較水面與航空領域的被動定位難度更大，究其原因就是在強調隱蔽性的前提下水下可利用的觀測量更少，即便有可用的觀測量其解析形式也十分復雜，因而當今潛器的水下定位、導航多是基于特定的地理信息數字圖，通過采用各種匹配算法和觀測信息在數字圖上估算出潛器的位置坐標，這類方法的核心就是匹配算法，但是匹配算法在應用時為了保證具有好的定位效果，對潛器的初始誤差、運動模式具有諸多限制，且有的匹配算法還不能保證定位的實時性（如ICCP算法），能保證實時性的又容易出現誤匹配（如Sitan算法）。水面和航空領域的被動定位相對簡單，它只需要建立合適的目標運動狀態方程再由具體的觀測方程，將真實觀測量中所包含的我們感興趣的目標狀態（如經、緯度坐標）通過非線性濾波技術［11～14］估算出來。

建立潛器狀態方程如下：

其中

其中F（k）、Γ（k）為k時刻的狀態和擾動矩陣，φ（k）和λ（k）為k時刻的潛器經緯度坐標，Vφ、Vλ為潛器在緯度和經度方向上的運動速度，T為采樣間隔，v（k）為零均值過程白噪聲。在目標定位領域這一建模方式包含了目標所有的勻速和勻加速運動情況。

對于水下潛器定位的觀測值取為重力儀觀測重力

其中γ＝9.780325（1＋0.00530244sin2（φ）－0.00000587sin2（2φ））為 WGS84橢球正常重力公式。

設

其中

z（φ，λ）為水下定位時將要采用的量測值，式（7）就是本文最終建立的量測方程。

式（7）與經、緯度坐標是一種非線性關系，本文在對潛器位置誤差進行最優估計時采用的是經典的擴展卡爾曼濾波（EKF），因而需要對式（7）進行泰勒展開的線性變換。

將式（7）分別對φ、λ求導可得

其中

對式（7）進行泰勒展開并舍去泰勒展開高階項，并直接給出線性化后的量測方程如下

其中H（k）為線性化后的量測矩陣，W（k）為觀測噪聲。

由狀態方程（6）及線性化后的量測方程（9）通過非線性濾波算法就可對目標的水下位置進行最優估計。

4 仿真分析

以2.2節給出的區域范圍，i為117°E～121°E，j為21°N～25°N，分辨率為2′×2′的衛星測高重力異常數據作為已知重力基準圖進行仿真分析。水下潛器初始位置為（21.5°，121.7°），北向和東向航速分別為 9nmile／h 和 －5nmile／h，加速度都為1nmile／h2，整個仿真階段取為100個采樣點，采樣間隔為6min，總航行時間為10h；重力儀測量數據的仿真采用重力異常數據加入測量噪聲的方法。現今海洋重力儀的動態精度已經達到1mGal量級（如L＆R的S型海洋重力儀），因此重力測量誤差一般取方差為1mGal2的白噪聲，另外，考慮到各種數據濾波及誤差綜合影響（如交叉耦合、厄特沃什改正），取一項白噪聲作為濾波估計的誤差噪聲，這個綜合值根據經驗取為9mGal2。采用50次的蒙特卡羅仿真，結果如圖2～圖5所示。

表2 定位誤差統計表（單位：海里）

圖2、圖3為基于重力觀測值的水下被動定位示意圖，由兩圖可以看出估計航跡在大部分的航行階段都能較好地反映真實航跡，只是由于運動模式轉換在航行轉彎時有一定程度的波動，由以上分析可以得出結論：

圖2 重力被動定位示意圖

圖3 定位放大示意圖

圖4 緯度誤差示意圖

圖5 經度誤差示意圖

1）以觀測重力作為量測值是可以用來提取潛器位置坐標的，由圖2、圖3的被動定位示意圖可以看出，估計航跡在大部分的航行階段都能較好地反映真實航跡，這說明文章建立的量測方程是有效的，通過非線性濾波算法對量測值中的目標位置信息進行最優估計，可實現潛器的被動導航與定位。

2）圖4、圖5、表2給出的是本文重力被動定位方法的定位誤差及統計結果，在10h的航行階段最大緯度誤差不超過4.02海里，平均緯度誤差只有0.74海里，最大經度誤差2.28海里，平均經度誤差0.59海里，從以上數據及誤差示意圖可以得出結論，雖然經緯度誤差存在一定程度上的波動，但總體的定位誤差已經不隨時間積累，這使得該定位模式在水下可長期使用。

5 結語

被動定位在軍事領域具有極高的應用價值，文章將該原理引入潛器的水下被動定位，首先基于二維高斯樣條函數逼近的重力場模型建立所需的量測方程，隨后通過非線性濾波技術對水下潛器的位置坐標進行實時估計。該定位方法使用簡單，無需使用傳統的匹配算法，因而避免了匹配算法對潛器的初始環境以及慣導誤差等的種種限制，最終的仿真結果證明，直接由實測重力數據對水下潛器進行定位是可行的，它的定位誤差不隨時間積累，另外隨著海洋重力數據測量精度的不斷提高，最終的定位精度也必將得到進一步的提高。

［1］Rice H，Mendelsohn L，Robert Aarons R，et al.Next generation marine precision navigation system［C］／／IEEE：Position Location and Navigation Symposium，2000：200－205.

［2］Moryl，J，Rice，H，Shinners，S.The Universal Gravity Module for Enhanced Submarine Navigation［C］／／IEEE Position Location and Navigation Symposium（PLANS'98）Record，20－23April 1998.

［3］Hollowell J.Heli／SITAN：a terrain referenced navigation algorithm for helicopters［C］／／IEEE Position，Location，and Navigation Symposium 1990（PLANS'90），Las Vegas，NV，USA，March 20－23，1990，616：625.

［4］許大欣.利用重力異常匹配技術實現潛艇導航［J］.地球物理學報，2005，48（4）：812－816.

［5］Behzad K.P，Behrooz K.P.Vehicle location on gravity maps［J］.Proceedings of SPIE－The International Society for Optical Engineering，1999：191.

［6］劉繁明，孫楓，成怡.基于ICCP算法及其推廣的重力定位［J］.中國慣性技術學報，2004，12（5）：36－39.

［7］Bishop G.C.Gravitational field maps and navigational errors［J］.IEEE Journal of Oceanic Engineering，2002，27（3）：726－737.

［8］BIAN Shaofeng，JOACHIM M.determining the parameter of a covariance function by analytical rules［J］.ZfV，1999（7）：212－216.

［9］童余德，邊少鋒，蔣東方，等.基于高斯樣條函數的局部重力異常場解析重構［J］.測繪學報，2012，41（5）：754－760.

［10］Overholt K J.Efficiency of the Fibonacci search method［J］.BIT，1973，13（1）：92－96.

［11］楊元喜.自適應動態導航定位［M］.北京：測繪出版社，2006.

［12］王丹力，張洪鉞.慣導系統初始對準的非線性濾波算法［J］.中國慣性技術學報，1999，7（3）：17－21.

［13］吳俊偉，曾啟明，聶莉娟.慣性導航系統的誤差估計［J］.中國慣性技術學報，2002，10（6）：1－5.

［14］秦永元，張洪鉞，汪叔華.卡爾曼濾波與組合導航原理［M］.西安：西北工業大學出版社，1998.