王志剛 顧雪峰
(1.91550部隊 大連 116023)(2.海軍工程大學(xué)兵器工程系 武漢 430033)(3.海軍工程大學(xué)兵器科研部 武漢 430033)
?
基于二維高斯樣條函數(shù)的水下地磁被動定位
王志剛1,2顧雪峰3
(1.91550部隊 大連 116023)(2.海軍工程大學(xué)兵器工程系 武漢 430033)(3.海軍工程大學(xué)兵器科研部 武漢 430033)
目前利用水下地磁信息對慣導(dǎo)誤差進(jìn)行校正多以各類匹配算法為核心,對此論文提出一種基于二維高斯樣條函數(shù)的水下地磁被動定位新模式。文章首先介紹了一種利用二維高斯樣條函數(shù)逼近的連續(xù)局部地磁場模型,隨后在該模型的基礎(chǔ)上將實測地磁表示為連續(xù)的解析形式,最后以實測地磁作為包含目標(biāo)位置信息的量測值,并結(jié)合擴(kuò)展卡爾曼濾波算法對目標(biāo)位置進(jìn)行最優(yōu)估計。這種方法無需使用常用的匹配算法,因而也就擺脫了匹配算法的諸多限制。以分辨率為2′×2′的某區(qū)域地磁異常數(shù)據(jù)為背景場進(jìn)行仿真,最終的仿真結(jié)果表明:經(jīng)高斯樣條函數(shù)逼近的局部地磁場模型平均誤差小于0.26nT,水下平均經(jīng)、緯定位誤差分別小于0.96和0.33海里。
被動定位; 高斯樣條函數(shù); 局部地磁場; 匹配算法
Class Number P207
目前,對于水下潛器的被動定位是導(dǎo)航、定位領(lǐng)域研究的一個熱點和難點問題。水下潛器的慣導(dǎo)系統(tǒng)(Inertial Navigation System,INS)誤差隨時間積累,因此必須要通過其它導(dǎo)航方式實時或定期修正INS。出于隱蔽性要求,水下潛器又很難利用衛(wèi)星或無線電信息,此時利用水下地理特征信息的輔助導(dǎo)航定位技術(shù)就是很好的選擇。常用的輔助導(dǎo)航定位系統(tǒng)有地形、重力、地磁輔助導(dǎo)航系統(tǒng)[1~4]等,目前這些輔助導(dǎo)航系統(tǒng)都是以各種匹配算法為核心來獲得潛器的最佳匹配位置(真實位置)。傳統(tǒng)的匹配算法主要分為單點迭代和序列迭代兩類,分別以SITAN[5~6](Sandia Inertial Terrain-Aided Navigation)和ICCP[7~8](Interval Closest Contour Point)算法為典型代表。這些匹配算法分別在實時性和穩(wěn)定性方面有其優(yōu)點,但是為了保證良好的匹配效果,算法對潛器的運動規(guī)律和初始誤差環(huán)境有很多嚴(yán)格的限制[9~11],因而也限制了其應(yīng)用范圍。針對這一問題,本文提出一種直接利用觀測地磁進(jìn)行水下被動定位的新模式,應(yīng)用該方法的前提就是確定地磁異常與潛器位置(大地經(jīng)、緯度坐標(biāo))的精確解析表達(dá)式,之后與航空領(lǐng)域的目標(biāo)被動定位原理相似,由于觀測地磁中包含了目標(biāo)的位置坐標(biāo)信息,直接通過非線性濾波算法,就可對潛器的航行位置進(jìn)行最優(yōu)估計。
直接利用地磁數(shù)據(jù)進(jìn)行定位,在很大程度上取決于觀測地磁與目標(biāo)所在位置關(guān)系的確定,也就是全球地磁場模型的確定。目前常用的各類全球地磁場模型[12]無論在階次和精度方面都已經(jīng)達(dá)到了很高的水平,它本質(zhì)上屬于球坐標(biāo)系下的調(diào)和分析,因而應(yīng)用過程中需要計算大量的Legendre系數(shù),使得模型計算量巨大,而這對于定位的實時性是一大考驗,另外現(xiàn)有的大部分全球地磁場模型非線性比較嚴(yán)重,也不太利于濾波估計時的線性化處理。
近年來,高斯樣條函數(shù)在醫(yī)學(xué)影像、數(shù)字地形構(gòu)建、計算機(jī)圖形學(xué)及地球物理等領(lǐng)域有較為廣泛的應(yīng)用,文獻(xiàn)[13]討論了高斯樣條在協(xié)方差函數(shù)代數(shù)確定中的應(yīng)用,文獻(xiàn)[14]采用高斯函數(shù)作為樣條基函數(shù)對區(qū)域重力進(jìn)行二維整體逼近,文獻(xiàn)[15]討論了高斯樣條函數(shù)逼近局部重力的相關(guān)影響因素。高斯樣條函數(shù)運算簡單且最終解析式是統(tǒng)一的,文章將其引入用于地磁建模,采用高斯樣條函數(shù)逼近局部地磁場,由某區(qū)域的試算結(jié)果可以看出函數(shù)逼近的平均誤差只有0.26nT。鑒于二維高斯樣條函數(shù)的解析形式以及高精度,使得該模型非常適合在估計潛器位置時作為量測方程使用。
2.1 二維高斯樣條函數(shù)逼近數(shù)學(xué)原理
=span{Gx((x-xi)/Δx)Gy((y-yj)/Δy)}
(1)
即有
×Gy((y-yi)/Δy)
(2)
不妨設(shè)
(3)
則上式可以簡化為矩陣形式:
(4)
則根據(jù)插值條件{L(xi,yj)=zi,j|i=1,…,m;j=1,…,n},有如下線性方程:
(5)
易知X、Y均為非奇異矩陣,即式(5)有唯一解,解方程即得到系數(shù)矩陣C,將系數(shù)矩陣代入式(2)便得到該局部地磁異常基準(zhǔn)圖二維高斯樣條函數(shù)逼近解析式。
由式(5)可以看出系數(shù)矩陣C由X、Y與Z唯一確定,而X、Y由矩陣的階數(shù)和ax、ay決定。文獻(xiàn)[14]曾對求逆計算誤差與a和計算階數(shù)的關(guān)系進(jìn)行過詳盡的探討。經(jīng)討論可知:系數(shù)矩陣C的求解過程中如X、Y矩陣階數(shù)過高且ax、ay過小或階數(shù)過低且ax、ay過大時,X、Y矩陣的求逆存在較大誤差。而當(dāng)X、Y矩陣階數(shù)小于5或ax、ay值小于5時,矩陣求逆運算誤差較小。考慮到求逆計算誤差以及實際應(yīng)用中用到的格網(wǎng)區(qū)域遠(yuǎn)大于5×5,即可采用文獻(xiàn)[15]給出的斐波那契數(shù)列方法對ax∈[1,5],ay∈[1,5]區(qū)間內(nèi)進(jìn)行尋優(yōu)。
2.2 二維高斯樣條函數(shù)逼近精度分析
取范圍i為117°E-121°E,j為21°N-25°N,分辨率為2′×2′的地磁異常場數(shù)據(jù)作為已知地磁基準(zhǔn)圖。文章在該已知2′×2′地磁異常基準(zhǔn)圖的基礎(chǔ)上,取4′×4′格網(wǎng)點處地磁異常數(shù)據(jù)作為基礎(chǔ)數(shù)據(jù)點,進(jìn)行二維高斯樣條函數(shù)逼近,獲得該范圍局部地磁異常基準(zhǔn)圖解析表達(dá)式,并據(jù)此計算分辨率為2′×2′的地磁異常逼近值,然后與已知的2′×2′格網(wǎng)處地磁異常值進(jìn)行比對分析。




圖1 局部地磁異常基準(zhǔn)圖二維高斯逼近仿真實驗
最終的仿真結(jié)果如圖1所示,其中圖1(a)為已知2′×2′地磁異常基準(zhǔn)圖,圖1(b)為由4′×4′地磁異常基準(zhǔn)圖逼近的2′×2′地磁異常基準(zhǔn)圖,由此可以看出二維高斯樣條函數(shù)逼近的局部地磁異常基準(zhǔn)圖解析式能較好地描述這一范圍的地磁異常基準(zhǔn)圖。
圖1(c)、圖1(d)分別為逼近誤差圖及相應(yīng)的誤差統(tǒng)計圖,表1為誤差統(tǒng)計表,由以上圖、表可以得出,試算區(qū)域4′×4′格網(wǎng)點處二維高斯逼近計算值絕對誤差均值為0.26nT且98.21%計算點絕對誤差小于1nT,逼近精度可以滿足匹配導(dǎo)航要求。鑒于地磁二維高斯樣條逼近函數(shù)的連續(xù)形式,與經(jīng)、緯度坐標(biāo)具有明確的解析關(guān)系以及自身的高精度,使得這一逼近函數(shù)十分適合用于地磁被動定位時所需的量測方程。

表1 計算值和真值比較(單位:nT)
水下被動定位較水面與航空領(lǐng)域的被動定位難度更大,究其原因就是在強(qiáng)調(diào)隱蔽性的前提下水下可利用的觀測量更少,即便有可用的觀測量其解析形式也十分復(fù)雜,因而當(dāng)今潛器的水下定位、導(dǎo)航多是基于特定的地理信息數(shù)字圖,通過采用各種的匹配算法和觀測信息在數(shù)字圖上估算出潛器的位置坐標(biāo),這類方法的核心就是匹配算法,但是匹配算法在應(yīng)用時為了保證具有好的定位效果,對潛器的初始誤差、運動模式具有諸多限制,且有的匹配算法還不能保證定位的實時性(如ICCP算法),能保證實時性的又容易出現(xiàn)誤匹配(如Sitan算法)。水面和航空領(lǐng)域的被動定位相對簡單,它只需要建立合適的目標(biāo)運動狀態(tài)方程再由具體的觀測方程,就能將真實觀測值中所包含的我們感興趣的目標(biāo)狀態(tài)(如經(jīng)、緯度坐標(biāo))通過非線性濾波技術(shù)[14~17]估算出來。
建立潛器狀態(tài)方程如式(6):
x(k+1)=F(k)x(k)+Γ(k)v(k)
(6)
其中

對于水下潛器定位的觀測值取為地磁儀觀測地磁異常
z(φ,λ)=Δg(φ,λ)
(7)
其中
z(φ,λ)為水下定位時將要采用的量測值,式(7)就是本文最終建立的量測方程。
式(7)與經(jīng)、緯度坐標(biāo)是一種非線性關(guān)系,本文在對潛器位置誤差進(jìn)行最優(yōu)估計時采用的是經(jīng)典的擴(kuò)展卡爾曼濾波[18~20](EKF),因而需要對式(7)進(jìn)行泰勒展開的線性變換。
將式(7)分別對φ,λ求導(dǎo)可得
(8)
其中
對式(7)進(jìn)行泰勒展開并舍去泰勒展開高階項,直接給出線性化后的量測方程如下:
z(φ(k),λ(k)) ≈H(k)x(k)+W(k)
(9)其中H(k)為線性化后的量測矩陣,W(k)為觀測噪聲。
由狀態(tài)方程(6)及線性化后的量測方程(9)通過非線性濾波算法就可對目標(biāo)的水下位置進(jìn)行最優(yōu)估計。
以2.2節(jié)給出的區(qū)域范圍i為117°E-121°E,j為21°N-25°N,分辨率為2′×2′的地磁異常數(shù)據(jù)作為已知地磁基準(zhǔn)圖進(jìn)行仿真分析。水下潛器初始位置為(21.5°,121.7°),北向和東向航速分別為9nmile/h和-5nmile/h,加速度都為1nmile/h2,整個仿真階段取為150個采樣點,采樣間隔為3min,總航行時間為450min;地磁儀測量數(shù)據(jù)的仿真采用地磁異常數(shù)據(jù)加入測量噪聲的方法,考慮到各種測量及地圖誤差,取一項白噪聲作為濾波估計的誤差噪聲,這個綜合值根據(jù)經(jīng)驗取為10nT2。采用50次的蒙特卡羅仿真,結(jié)果如圖2~圖5所示。

圖2 地磁被動定位示意圖

圖3 定位放大示意圖

圖4 緯度誤差示意圖

圖5 經(jīng)度誤差示意圖
圖2、圖3為基于地磁觀測值的水下被動定位示意圖,由兩圖可以看出估計航跡在大部分的航行階段都能較好地反映真實航跡,只是由于運動模式轉(zhuǎn)換在航行轉(zhuǎn)彎時有一定程度的波動,由以上分析我們可以得出結(jié)論:

表2 定位誤差統(tǒng)計表(單位:海里)
1) 以觀測地磁作為量測值是可以用來提取潛器位置坐標(biāo)的,由圖2、圖3圖的被動定位示意圖可以看出,估計航跡在大部分的航行階段都能較好地反映真實航跡,這說明文章建立的量測方程是有效的,通過非線性濾波算法對量測值中的目標(biāo)位置信息進(jìn)行最優(yōu)估計,可實現(xiàn)潛器的被動導(dǎo)航與定位。
2) 圖5、圖6、表2給出的是地磁被動定位誤差及統(tǒng)計結(jié)果,在450min的航行階段最大緯度誤差不超過6.68海里,平均緯度誤差只有0.96海里,最大經(jīng)度誤差4.07海里,平均經(jīng)度誤差0.33海里,從以上數(shù)據(jù)及誤差示意圖可以得出結(jié)論,雖然經(jīng)、緯度誤差存在一定程度上的波動,但總體的定位誤差已經(jīng)不隨時間積累,這使得該定位模式在水下可長期使用。
被動定位在軍事領(lǐng)域具有極高的應(yīng)用價值,本文將該原理引入潛器的水下被動定位,首先基于二維高斯樣條函數(shù)逼近的地磁場模型建立所需的量測方程,隨后通過非線性濾波技術(shù)對水下潛器的位置坐標(biāo)進(jìn)行實時估計。該定位方法使用簡單,無需使用傳統(tǒng)的匹配算法,因而避免了匹配算法對潛器的初始環(huán)境以及慣導(dǎo)誤差等的種種限制,最終的仿真結(jié)果證明,直接由實測地磁數(shù)據(jù)對水下潛器進(jìn)行定位是可行的,且定位誤差不隨時間積累,另外隨著海洋地磁測量精度的不斷提高,最終的定位精度也必將得到進(jìn)一步的提高。
[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,1998:20-23.
[3] 蔡兆云,魏海平,任志新.水下地磁導(dǎo)航技術(shù)研究綜述[J].尖端科技,2007,1(3):28-30.
[4] 趙建虎,王勝平,王愛學(xué).一種改進(jìn)型TERCOM水下地磁匹配導(dǎo)航算法[J].武漢大學(xué)學(xué)報·信息科學(xué)版,2009,34(11):1320-1323.
[5] 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.
[6] 許大欣.利用重力異常匹配技術(shù)實現(xiàn)潛艇導(dǎo)航[J].地球物理學(xué)報,2005,48(4):812-816.
[7] Behzad K. P, Behrooz K. P. Vehicle location on gravity maps[C]//Proceedings of SPIE-The International Society for Optical Engineering,1999,3693,182:191.
[8] 劉繁明,孫楓,成怡.基于ICCP算法及其推廣的重力定位[J].中國慣性技術(shù)學(xué)報,2004,12(5):36-39
[9] Bishop G. C. Gravitational field maps and navigational errors[J]. IEEE Journal of Oceanic Engineering,2002,27(3):726-737.
[10] 王向磊.地磁匹配導(dǎo)航中幾項關(guān)鍵技術(shù)研究[J].測繪工程,2011,20(1):1-5.
[11] 喬玉坤,王仕成,張琪.地磁匹配制導(dǎo)技術(shù)應(yīng)用于導(dǎo)彈武器系統(tǒng)的制約因素分析[J].控制與制導(dǎo),2006,5(8):39-41.
[12] 彭富清.地磁模型與地磁導(dǎo)航[J].海洋測繪,2006,26(2):73-75.
[13] BIAN Shaofeng, JOACHIM M. determining the parameter of a covariance function by analytical rules[J]. ZfV,1999(7):212-216.
[14] 童余德,邊少鋒,蔣東方,等.基于高斯樣條函數(shù)的局部重力異常場解析重構(gòu)[J].測繪學(xué)報,2012,41(5):754-760.
[15] Overholt K J. Efficiency of the Fibonacci search method[J]. BIT,1973,13(1):92-96.
[16] 楊元喜.自適應(yīng)動態(tài)導(dǎo)航定位[M].北京:測繪出版社,2006.
[17] 王丹力,張洪鉞.慣導(dǎo)系統(tǒng)初始對準(zhǔn)的非線性濾波算法[J].中國慣性技術(shù)學(xué)報,1999,7(3):17-21.
[18] 吳俊偉,曾啟明,聶莉娟.慣性導(dǎo)航系統(tǒng)的誤差估計[J].中國慣性技術(shù)學(xué)報,2002,10(6):1-5.
[19] 秦永元,張洪鉞,汪叔華.卡爾曼濾波與組合導(dǎo)航原理[M].西安:西北工業(yè)大學(xué)出版社,1998.
[20] 晏登洋,任建新,宋永軍.慣性/地磁組合導(dǎo)航技術(shù)研究[J].機(jī)械與電子,2007,6(1):19-22.
Underwater Geomagnetic Passive Location Based on the 2-D Gauss Spline Function
WANG Zhigang1,2GU Xuefeng3
(1. No. 91550 Troop of PLA, Dalian 116023) (2. Department of Weapon, Naval University of Engineering, Wuhan 430033) (3. Department of Science, Naval University of Engineering, Wuhan 430033)
As varies of matching algorithms have been used as the major methods for correcting the errors of INS on underwater geographic information, a new pattern of geomagnetic passive location which is based on 2-D Gauss spline function is given in this paper. Firstly, the algorithm to reconstruct the local grid geomagnetic base-map with the 2-D Gauss spline function is introduced, then by which the measured geomagnetic is expressed as continual analytical equation. Finally, the classical extended Kalman filter is introduced to estimate the target’s position with measured geomagnetic as measurement containing the information of target’s position. This method can be applied without using any matching algorithms, as a result it is free from the influence of the restriction by using these matching algorithms. Simulation is done on 0.2′×0.2′ geomagnetic anomaly database, and the simulation result proves that the mean error of analytic equation of local geomagnetic anomalies is less than 0.26nT, and the mean location errors in longitude and latitude are 0.96 and 0.33 sea mile respectively.
passive location, Gauss spline function, local geomagnetic field, matching algorithm
2015年1月13日,
2015年2月20日 基金項目:國家自然科學(xué)基金項目(編號:41201478,41374018);中國博士后科學(xué)基金項目(編號:2013M542533)資助。作者簡介:王志剛,男,工程師,研究方向:水下匹配輔助導(dǎo)航與非線性濾波技術(shù)。顧雪鋒,男,工程師,研究方向:非線性濾波技術(shù)與重力,地磁場匹配輔助導(dǎo)航。
P207
10.3969/j.issn1672-9730.2015.07.020