黃 炎,李?yuàn)檴?李新星,宋星光,范 雕,萬(wàn)宏發(fā)

1.軍事科學(xué)院國(guó)防科技創(chuàng)新研究院,北京 100071; 2.信息工程大學(xué)地理空間信息學(xué)院,河南 鄭州 450001; 3.北京航天飛行控制中心,北京 100094

水下潛航器搭載的慣性導(dǎo)航系統(tǒng)(inertial navigation system,INS)因慣性元器件等誤差導(dǎo)致的定位誤差隨時(shí)間積累,定位精度逐漸發(fā)散。為了抑制慣性導(dǎo)航系統(tǒng)誤差發(fā)散,同時(shí)不造成潛航器隱蔽性能損失,國(guó)內(nèi)外學(xué)者對(duì)各類校準(zhǔn)慣性導(dǎo)航的方法進(jìn)行了研究,其中利用海洋先驗(yàn)重力場(chǎng)信息進(jìn)行匹配導(dǎo)航的思想得到了廣大學(xué)者的認(rèn)同。各類重力匹配算法中,基于遞推估計(jì)最優(yōu)濾波原理的單點(diǎn)迭代匹配算法(尤其是Kalman濾波匹配算法)以其良好的實(shí)時(shí)性和在低信噪比條件下的高精度特性被廣泛研究和應(yīng)用[1-3]。但是,在采用Kalman濾波方法匹配導(dǎo)航過(guò)程中隨著重力觀測(cè)量的增加,若僅采用傳統(tǒng)Kalman濾波進(jìn)行計(jì)算,則各類觀測(cè)量的量測(cè)誤差會(huì)相互影響,降低濾波精度,同時(shí)也不方便對(duì)量測(cè)故障進(jìn)行診斷。因此,本文采用聯(lián)邦Kalman濾波(federated Kalman filter,FKF)對(duì)多觀測(cè)量進(jìn)行分散降權(quán)處理,將每一個(gè)獨(dú)立的重力觀測(cè)量作為一個(gè)子系統(tǒng)的觀測(cè)值,從而避免某個(gè)觀測(cè)量出現(xiàn)粗差時(shí)影響整個(gè)系統(tǒng)的濾波精度,同時(shí)也便于對(duì)各個(gè)子系統(tǒng)的故障問(wèn)題進(jìn)行隔離與診斷。

為保證濾波精度,還需要建立較為準(zhǔn)確的子濾波器量測(cè)方程,該過(guò)程需要對(duì)重力量測(cè)值及其變化率進(jìn)行局部建模,而構(gòu)建模型的精確與否又直接影響濾波結(jié)果精度。目前,應(yīng)用于重力匹配導(dǎo)航的局部重力場(chǎng)構(gòu)建方法有平面擬合法、雙二次曲面法和重力場(chǎng)球諧模型法等[4-6]。其中平面擬合法效率較高,但是僅是線性展開(kāi)精度難以保證;在此基礎(chǔ)上文獻(xiàn)[4]使用雙二次曲面對(duì)局部重力場(chǎng)進(jìn)行擬合,計(jì)算較為高效,但由于展開(kāi)階次的限制,恢復(fù)重力場(chǎng)的信息仍然不夠精細(xì);文獻(xiàn)[5]使用傅里葉階數(shù)表示局部重力場(chǎng)的連續(xù)解析形式,但該方法會(huì)由于添零導(dǎo)致部分區(qū)域精度較差,且計(jì)算效率不高;文獻(xiàn)[6]使用重力場(chǎng)球諧模型計(jì)算局部重力場(chǎng)信息,使用并行技術(shù)提高計(jì)算效率,盡管在階次方面達(dá)到了應(yīng)用需求,但是對(duì)于匹配導(dǎo)航局域性、實(shí)時(shí)性特點(diǎn)而言,重力場(chǎng)球諧模型計(jì)算復(fù)雜、耗時(shí)長(zhǎng)且在局部重力場(chǎng)的精細(xì)結(jié)構(gòu)和高頻信息的表述上則顯得精度不足,盡管通過(guò)并行方法可提高其計(jì)算效率,但并行計(jì)算設(shè)備增加了成本和潛航器計(jì)算機(jī)功耗。

局部重力場(chǎng)的譜方法是地球重力場(chǎng)理論的研究重點(diǎn),如何構(gòu)造合適的譜函數(shù)并對(duì)其進(jìn)行快速有效計(jì)算是核心問(wèn)題。本文基于局部重力場(chǎng)球冠諧分析(spherical cap harmonic analysis,SCHA)方法,采用海洋重力場(chǎng)先驗(yàn)信息,使用改進(jìn)球諧分析技術(shù)(adjusted spherical harmonic analysis,ASHA)快速構(gòu)建局部重力場(chǎng)球冠諧模型,并在此基礎(chǔ)上組建FKF子濾波器非線性量測(cè)方程,實(shí)現(xiàn)對(duì)潛航器慣性導(dǎo)航系統(tǒng)的誤差估計(jì)。

1 FKF濾波器構(gòu)建

聯(lián)邦Kalman濾波器是一種典型的分散式濾波器[7],由一個(gè)參考系統(tǒng)(主濾波器)和若干個(gè)子系統(tǒng)(局部濾波器)組成,屬于兩步級(jí)聯(lián)的分散式濾波系統(tǒng)[8]。根據(jù)慣性/重力/重力梯度組合導(dǎo)航原理,其FKF濾波器結(jié)構(gòu)如圖1所示。

圖1 重力匹配FKF濾波結(jié)構(gòu)

1.1 濾波觀測(cè)量選取

利用重力和重力梯度信息作為子濾波器觀測(cè)量進(jìn)行濾波匹配時(shí),重力異常和擾動(dòng)重力梯度共有6個(gè)獨(dú)立分量(Δg、Tzz、Txx(或Tyy)、Tzx、Tzy和Txy)。根據(jù)適定邊值問(wèn)題理論,Δg和Tzz分量能夠反映完整的重力場(chǎng)信息,可以直接作為觀測(cè)量進(jìn)行濾波匹配;而Txx(或Tyy)、Tzx、Tzy和Txy則無(wú)法單獨(dú)使用。為此引入復(fù)數(shù)坐標(biāo),根據(jù)適定邊值問(wèn)題理論,獲得包含完整重力場(chǎng)信息的重力矢量與擾動(dòng)重力梯度的復(fù)數(shù)組合形式。

以局部北東天坐標(biāo)系為基準(zhǔn),定義復(fù)數(shù)坐標(biāo)系(i-,i0,i+),則其可用北東天坐標(biāo)系坐標(biāo)(X,Y,Z)表示為[9-10]

(1)

因此在復(fù)數(shù)坐標(biāo)系下,擾動(dòng)重力矢量?T和擾動(dòng)重力梯度張量?·?T可表示為

(2)

式中,T表示擾動(dòng)位。

1.2 主濾波器狀態(tài)方程

重力/重力梯度信息對(duì)速度誤差的反應(yīng)不夠敏感,故在進(jìn)行重力匹配導(dǎo)航時(shí),僅以位置誤差作為狀態(tài)變量。又由于慣性系統(tǒng)的垂直通道是不穩(wěn)定的,因此高程信息常使用單獨(dú)通道進(jìn)行測(cè)量,一般通過(guò)測(cè)深測(cè)潛儀等來(lái)改進(jìn)垂直通道的不穩(wěn)定性。因此本文選用X=(δφ,δλ)T作為濾波狀態(tài)變量。由慣性導(dǎo)航誤差方程和濾波狀態(tài)變量可知,慣性導(dǎo)航系統(tǒng)平面定位誤差微分方程表達(dá)式為[6]

(3)

式中,δφ和δλ分別表示緯向和經(jīng)向誤差;h表示高度;ve表示東向速度;RN表示指示位置點(diǎn)卯酉圈曲率半徑;φ表示地心緯度;Wφ與Wλ表示由系統(tǒng)噪聲產(chǎn)生的緯向和經(jīng)向誤差。

Xk=Φk/k-1Xk-1+Wk-1

(4)

式中,Xk表示濾波狀態(tài)向量;Φk/k-1表示狀態(tài)轉(zhuǎn)移矩陣,可由系統(tǒng)狀態(tài)矩陣F通過(guò)Laplace變換得到Φk/k-1=I+tF,t表示濾波周期。

1.3 子濾波器量測(cè)方程

由于重力量測(cè)值與狀態(tài)參數(shù)之間呈非線性關(guān)系,對(duì)于非線性系統(tǒng)常通過(guò)構(gòu)建局部重力場(chǎng)解析函數(shù),獲得匹配點(diǎn)重力場(chǎng)值在x和y方向上的變化率,而后利用慣導(dǎo)指示位置在重力基準(zhǔn)圖上的比對(duì)值、潛航器搭載的重力儀測(cè)量數(shù)據(jù)經(jīng)過(guò)計(jì)算得到的實(shí)際重力值之差及其變化率構(gòu)建濾波量測(cè)方程,如圖2所示。

由于地球重力場(chǎng)是一個(gè)與位置有關(guān)的連續(xù)物理場(chǎng),因此,假設(shè)水下潛航器真實(shí)位置的重力異常Δg和擾動(dòng)重力梯度組合在慣導(dǎo)指示位置點(diǎn)(φj,λj)的鄰域內(nèi)存在一階導(dǎo)數(shù),利用泰勒級(jí)數(shù)將其展開(kāi)可得[4]

(5)

式中,vgj、vTzzj、vTzxj與vTzyj表示泰勒級(jí)數(shù)展開(kāi)過(guò)程中的截?cái)嗾`差;ΔgM(φj,λj)、TzzM(φj,λj)、TzxM(φj,λj)與TzyM(φj,λj)表示利用慣導(dǎo)指示位置坐標(biāo)在重力異常和擾動(dòng)重力梯度基準(zhǔn)圖上提取出的重力異常值和擾動(dòng)重力梯度值;Δg(φ,λ)、Tzz(φ,λ)、Tzx(φ,λ)與Tzy(φ,λ)則表示水下潛航器真實(shí)位置處的重力異常值和擾動(dòng)重力梯度值。

又因?yàn)樗聺摵狡髯陨泶钶d有重力儀和重力梯度儀,因此,Δg(φ,λ)與Tzz(φ,λ)又可表示為[4]

(6)

式中,vgs、vTzzs、vTzxs與vTzys表示由于重力儀和重力梯度儀測(cè)量產(chǎn)生的觀測(cè)誤差;Δgs(φj,λj)、Tzzs(φj,λj)、Tzxs(φj,λj)與Tzys(φj,λj)表示由重力儀和重力梯度儀測(cè)量并計(jì)算得到的重力異常和擾動(dòng)重力梯度分量觀測(cè)值。

聯(lián)合式(5)、式(6)可得到重力異常和擾動(dòng)重力梯度組合的匹配導(dǎo)航的濾波量測(cè)方程分別為

(7)

(8)

(9)

式中,Pk.k-1表示一步預(yù)測(cè)誤差方陣;Rk*表示量測(cè)噪聲方差陣;δ為閾值常數(shù),本文取δ=2.5。

(10)

(11)

2 基于SCHA建模的局部重力場(chǎng)模型

由于采用了多個(gè)重力觀測(cè)量作為不同子濾波器的量測(cè)值,因此若直接使用離散格網(wǎng)數(shù)據(jù)對(duì)觀測(cè)值進(jìn)行局部建模,則各個(gè)觀測(cè)量需要通過(guò)各自數(shù)據(jù)分別進(jìn)行建模,計(jì)算效率受到損失。當(dāng)研究區(qū)近似一個(gè)球冠時(shí),球冠諧函數(shù)可以作為該區(qū)域的譜函數(shù)[17-18]。而使用SCHA建模方法,可以有效利用重力場(chǎng)各個(gè)擾動(dòng)場(chǎng)元量與擾動(dòng)位之間的泛函關(guān)系,僅需要進(jìn)行一次建模即可,減少整體計(jì)算量,提高匹配算法實(shí)時(shí)性。球冠諧函數(shù)相較于球諧函數(shù),球冠諧函數(shù)使用非整階勒讓德函數(shù)替代整階勒讓德函數(shù),余緯θ的取值范圍由[0,π]變化為[0,θ0](θ0表示球冠半角),同時(shí)將地球自轉(zhuǎn)軸由北極點(diǎn)沿子午線方向旋轉(zhuǎn)至球冠區(qū)域中心點(diǎn)(圖3)。由于球冠諧分析建模區(qū)域范圍較小,因此球冠諧函數(shù)可以用較少的位系數(shù)反映較高分辨率的重力場(chǎng)信息。例如,為了達(dá)到某一分辨率,利用球諧分析(spherical harmonic analysis,SHA)建立全球重力場(chǎng)球諧模型需要使用NOSHA個(gè)位系數(shù),而在半角為θ0的球冠區(qū)域內(nèi)使用球冠諧分析建立球冠諧模型僅需使用NSCHA個(gè)位系數(shù),NOSHA與NSCHA的關(guān)系為[18]

圖3 球冠諧坐標(biāo)系

(12)

式中,Scap表示球冠的表面積;Searth表示地球的表面積。

(13)

2.1 廣義勒讓德函數(shù)

擾動(dòng)位T在球冠坐標(biāo)系下的級(jí)數(shù)表達(dá)式為

(14)

廣義勒讓德函數(shù)是指階數(shù)為實(shí)數(shù)的勒讓德函數(shù),包括整階、非整階勒讓德函數(shù),直接給出規(guī)格化的非整階勒讓的函數(shù)計(jì)算公式為[18]

(15)

式中,Jmax表示級(jí)數(shù)的最高階。Aj(lk,m)可通過(guò)遞推計(jì)算得到,遞推公式為[18]

(16)

式中,Klkm表示規(guī)格化因子。Klkm近似計(jì)算公式為[18]

(17)

2.2 基于移動(dòng)窗口的球冠諧建模

相對(duì)于球諧分析,球冠諧分析建模就是在所選中的某一地球球冠上進(jìn)行局部球諧分析。根據(jù)Helmholtz方程及其邊界條件,擾動(dòng)位在區(qū)域內(nèi)必須要滿足一定的邊界條件,在球冠極點(diǎn)(即余緯θc=0°)時(shí),球冠諧分析與球諧分析相同;而在球冠邊界處(即θc=θ0)時(shí),擾動(dòng)位需滿足以下邊界條件[18-19]

(18)

文獻(xiàn)[18—19]給出了式(18)的等價(jià)形式

(19)

不同的半角θ0與模型級(jí)數(shù)m對(duì)應(yīng)式(18)的解不同,且lk一般為非整數(shù),因而通過(guò)Helmholtz邊界條件反解非整階勒讓德函數(shù)階數(shù)的計(jì)算量較大、耗時(shí)較長(zhǎng),不利于重力匹配導(dǎo)航實(shí)時(shí)性需求。為提高計(jì)算效率,引入ASHA技術(shù)[20],該技術(shù)能夠在保證區(qū)域模型構(gòu)建精度的前提下減少計(jì)算時(shí)間,提高匹配算法實(shí)時(shí)性(圖4)。ASHA技術(shù)的核心是將余緯θc由原球冠范圍(0,θ0)映射到(0,π/2)上,即映射到半球上,達(dá)到將非整階勒讓德函數(shù)用整階勒讓德函數(shù)進(jìn)行替換的目的,從而化簡(jiǎn)勒讓德函數(shù)的階數(shù)求解過(guò)程,減少計(jì)算耗時(shí),提高計(jì)算效率。

圖4 ASHA技術(shù)

使用ASHA技術(shù),首先需要對(duì)原球冠坐標(biāo)系的坐標(biāo)進(jìn)行轉(zhuǎn)換,將其轉(zhuǎn)換至半球坐標(biāo)系[20]。具體公式為

(20)

式中,ρ′、λ′、θ′為半球坐標(biāo)系中的坐標(biāo);s表示坐標(biāo)轉(zhuǎn)換系數(shù)。

(21)

根據(jù)式(21)即可獲得球冠諧模型系數(shù)與重力異常觀測(cè)值間的誤差方程為

(22)

式中,g表示由DTU18全球重力異常模型獲得的觀測(cè)數(shù)據(jù)。

依據(jù)式(22)的誤差方程即可通過(guò)最小二乘法構(gòu)建局部重力場(chǎng)球冠諧模型。當(dāng)球冠半角θ0≤20°時(shí),ASHA技術(shù)均可以相當(dāng)?shù)木劝亚蚬谟蛴成涞桨肭蛏稀?/p>

在貧困心理學(xué)中，研究者們普遍認(rèn)為導(dǎo)致短視行為的原因是物質(zhì)資源缺乏對(duì)個(gè)體的認(rèn)知和情感層面產(chǎn)生的消極影響。在認(rèn)知層面，個(gè)體會(huì)因?yàn)榻?jīng)濟(jì)需要得不到滿足而產(chǎn)生壓力感和不安全感，產(chǎn)生無(wú)法從身邊的重要他人獲得支持的信念。在情感層面，個(gè)體會(huì)產(chǎn)生一系列的消極情感，包括焦慮、恐慌、抑郁、自我貶損以及自主性喪失感等等(Fabio & Maree, 2016)。這些因素也會(huì)消耗人們大量的心理資源，使個(gè)體沒(méi)有足夠的心理資源去處理其他任務(wù)，因而失去了為長(zhǎng)遠(yuǎn)打算的能力，表現(xiàn)為短視(Haushofer & Fehr, 2014)。

2.2.2 移動(dòng)窗口實(shí)時(shí)球冠諧模型

盡管球冠諧模型可以用較少的模型系數(shù)反映區(qū)域較為精細(xì)的重力場(chǎng)結(jié)構(gòu),但是隨著球冠半角和數(shù)據(jù)分辨率的增加,對(duì)應(yīng)球冠諧模型的階數(shù)也隨之增加。然而,高階球冠諧模型存在不穩(wěn)定性,當(dāng)階數(shù)lk序號(hào)k達(dá)到110以上時(shí)穩(wěn)定性急速下降[21]。

目前,海洋重力異常和擾動(dòng)重力梯度數(shù)據(jù)的分辨率可以達(dá)到1′[6],球冠諧模型階數(shù)應(yīng)與實(shí)際數(shù)據(jù)分辨率相契合,若數(shù)據(jù)分辨率為fs,則階數(shù)lk的最大序數(shù)Kmax可近似計(jì)算得到

(23)

通過(guò)式(23)易得分辨率為1′,Kmax=110時(shí),球冠諧模型能達(dá)到的最大球冠半角θmax約為55′。但是,水下潛器在進(jìn)行長(zhǎng)時(shí)間航行時(shí),航行范圍往往大于55′[22-24]。因此直接對(duì)整個(gè)航行區(qū)域進(jìn)行建模不僅建模范圍不足,而且隨著匹配點(diǎn)遠(yuǎn)離建模中心,計(jì)算精度也會(huì)有所損失。為有效解決該問(wèn)題,本文選擇以每一個(gè)匹配點(diǎn)為球冠中點(diǎn),以匹配點(diǎn)附近一定范圍區(qū)域重力異常作為觀測(cè)數(shù)據(jù),構(gòu)建以匹配點(diǎn)為中心的小范圍高分辨率實(shí)時(shí)移動(dòng)窗口球冠諧模型。

(24)

計(jì)算得到所有時(shí)刻置信概率為P的圓概率誤差半徑RPi后,整個(gè)慣性導(dǎo)航系統(tǒng)的圓概率誤差半徑RP取RPi中的最大值,即RP=max(RPi)。據(jù)此,將兩倍的RP作為構(gòu)建球冠諧模型時(shí)的移動(dòng)窗口最小半徑。

3 試驗(yàn)驗(yàn)證與分析

利用慣性/重力/重力梯度組合導(dǎo)航仿真程序進(jìn)行水下潛航器航跡與測(cè)量數(shù)據(jù)仿真,慣性導(dǎo)航采樣率為100 Hz。仿真數(shù)據(jù)包含慣性元器件輸出(角增量與速度增量)、重力異常觀測(cè)數(shù)據(jù)、擾動(dòng)重力梯度觀測(cè)數(shù)據(jù)及水下潛航器真實(shí)航跡。慣性元器件誤差設(shè)置見(jiàn)表1。當(dāng)前,海洋重力儀動(dòng)態(tài)測(cè)量精度已達(dá)到1 mGal(如德國(guó)的KSS系列和美國(guó)的L&R系列重力儀),因此重力儀測(cè)量誤差一般取標(biāo)準(zhǔn)差為1 mGal的白噪聲,此外還需要增加一個(gè)誤差項(xiàng)用以表示由于厄特弗斯改正等因素造成的影響,根據(jù)經(jīng)驗(yàn)取為3 mGal[11]。重力梯度儀的動(dòng)態(tài)測(cè)量精度可達(dá)到7 E(1 E=10-9/s2)(如美國(guó)的HD-AGG和英國(guó)的EGG系列重力梯度儀),因此重力梯度儀測(cè)量誤差取標(biāo)準(zhǔn)差為7 E的白噪聲[26]。重力異?；鶞?zhǔn)圖選用DTU利用衛(wèi)星測(cè)高數(shù)據(jù)、實(shí)測(cè)重力異常數(shù)據(jù)等反演計(jì)算發(fā)布的DTU18全球重力異常模型,數(shù)據(jù)分辨率1′,中誤差為2～3 mGal。擾動(dòng)重力梯度基準(zhǔn)圖則在EIGEN-6C4地球重力場(chǎng)模型與DTU18全球重力異常模型基礎(chǔ)上,基于Stokes理論利用移去-恢復(fù)技術(shù)計(jì)算得到,數(shù)據(jù)分辨率計(jì)算至1′,基于Txx+Tyy+Tzz=0的約束條件,計(jì)算結(jié)果Txx+Tyy+Tzz絕對(duì)值最大值為0.035 E,標(biāo)準(zhǔn)差為0.002 1 E,能夠一定程度上反映所構(gòu)建擾動(dòng)重力梯度基準(zhǔn)圖的正確性。同時(shí),由于重力觀測(cè)量的基準(zhǔn)圖數(shù)據(jù)是以格網(wǎng)形式離散存儲(chǔ),導(dǎo)致重力觀測(cè)量在基準(zhǔn)圖上的值存在半個(gè)格網(wǎng)的等值區(qū)域。若量測(cè)更新時(shí)間較短,則觀測(cè)值存在強(qiáng)相關(guān)性,導(dǎo)致量測(cè)方程不可觀測(cè);若量測(cè)更新時(shí)間較長(zhǎng),則慣導(dǎo)累積誤差較大難以校正。因此結(jié)合本文所用的基準(zhǔn)圖分辨率及水下潛航器的航行速度,FKF量測(cè)方程更新時(shí)間間隔選取180 s(即每隔180 s進(jìn)行一次組合匹配導(dǎo)航校正)。慣導(dǎo)解算采用雙子樣算法,每次慣導(dǎo)更新均進(jìn)行一次狀態(tài)方程更新,因此FKF狀態(tài)方程更新時(shí)間間隔為0.02 s。

表1 慣性元器件誤差設(shè)置

仿真試驗(yàn)區(qū)域大小為2.5°×4°(范圍為7°N—9.5°N,110.5°E—114.5°E),試驗(yàn)區(qū)重力異常與擾動(dòng)重力梯度張量基準(zhǔn)圖如圖5所示。

圖5 試驗(yàn)區(qū)基準(zhǔn)

選取5個(gè)重力場(chǎng)特征參數(shù),分別為標(biāo)準(zhǔn)差、緯向粗糙度、經(jīng)向粗糙度、緯向信息熵和經(jīng)向信息熵[27],用于表達(dá)區(qū)域重力場(chǎng)分布、描述區(qū)域代表性,并作為評(píng)價(jià)區(qū)域適配性的指標(biāo)。各個(gè)重力場(chǎng)分量的特征參數(shù)見(jiàn)表2。

表2 區(qū)域重力場(chǎng)特征參數(shù)

由表2可知,區(qū)域內(nèi)Δg、Tzz和Tzx的標(biāo)準(zhǔn)差數(shù)值較大,而Tzy則由于其自身數(shù)值較小,標(biāo)準(zhǔn)差相對(duì)較小;各個(gè)重力觀測(cè)量的緯向和經(jīng)向粗糙度量級(jí)相當(dāng),表明該區(qū)域內(nèi)各分量緯向和經(jīng)向的起伏變化相當(dāng);各個(gè)重力觀測(cè)量的經(jīng)向和緯向信息熵均接近0.5(即信息熵極大值),表明該區(qū)域各個(gè)重力觀測(cè)量特征變化明顯,符合重力匹配導(dǎo)航的應(yīng)用要求。

試驗(yàn)中,用于匹配導(dǎo)航的航跡(M航跡)航行時(shí)長(zhǎng)24 h,初始速度(0 m/s,0 m/s,0 m/s),初始姿態(tài)為(0°,0°,-10°),起始坐標(biāo)為(7°N,110°E,-100 m),慣導(dǎo)初始誤差設(shè)置為1.1 n mile。潛航器首先沿著北偏東45°方向先以0.008 6 m/s2的加速度勻加速至10 n mile/h,然后勻速直線航行12 h;而后右轉(zhuǎn)45°,再勻速航行12 h。

3.1 移動(dòng)窗口范圍選取

表3 不同置信概率圓概率誤差半徑與

由表3可知,慣性導(dǎo)航系統(tǒng)置信概率為95%的圓概率誤差半徑為4.852 n mile,約為5′,因此移動(dòng)窗口范圍最小應(yīng)選取半角為10′的球冠區(qū)域,并以此構(gòu)建移動(dòng)窗口球冠諧模型。

表4 誤差統(tǒng)計(jì)

由表4可知,使用不同移動(dòng)窗口半徑的濾波匹配結(jié)果精度近似,但隨著球冠半角的增加,相同分辨率對(duì)應(yīng)的階數(shù)序號(hào)最大值呈線性增長(zhǎng),球冠諧模型構(gòu)建耗時(shí)則呈指數(shù)增加,當(dāng)球冠半角為20′時(shí),模型構(gòu)建耗時(shí)超過(guò)10 s,不利于濾波的實(shí)時(shí)性計(jì)算。因此,在保證精度的前提下,顧及匹配導(dǎo)航實(shí)時(shí)性需求,移動(dòng)窗口應(yīng)選擇半徑為10′的球冠區(qū)域。

3.2 FKF-SCHA重力匹配算法匹配效果分析

表5 計(jì)算耗時(shí)統(tǒng)計(jì)

表6 24 h航行誤差統(tǒng)計(jì)

由以上仿真試驗(yàn)結(jié)果可知:

(1) 由表5可知,使用SHA法構(gòu)建量測(cè)方程耗時(shí)超過(guò)1 s,在使用基于GPU的并行計(jì)算方法進(jìn)行加速后,計(jì)算效率得到提高,PSHA法計(jì)算耗時(shí)為0.090 s,而本文所提移動(dòng)窗口SCHA法計(jì)算耗時(shí)與PSHA法計(jì)算耗時(shí)相當(dāng),為0.086 s。但是PSHA法需要使用GPU設(shè)備輔助提高計(jì)算效率,考慮到潛航器搭載的計(jì)算機(jī)硬件計(jì)算能力和成本,移動(dòng)窗口SCHA法則更具實(shí)用價(jià)值。

(2) 結(jié)合圖6(a)與圖6(d)可知,本文所提FKF-SCHA重力匹配算法經(jīng)過(guò)5次量測(cè)更新(15 min)即可對(duì)初始誤差進(jìn)行修正,且修正后的航跡相較于INS航跡更貼近于參考航跡,說(shuō)明本文算法是有效的,能夠?qū)T導(dǎo)誤差進(jìn)行一定程度的修正。

(3) 結(jié)合表6與圖6(b)、(c)可知,慣性導(dǎo)航緯向誤差呈周期性變化,經(jīng)向誤差則隨航行時(shí)間發(fā)散;在使用FKF-SCHA重力匹配算法后,緯向誤差和經(jīng)向誤差絕對(duì)值均值分別為0.298和0.158 n mile。盡管誤差曲線存在一定程度的波動(dòng),但均是以接近0的均值進(jìn)行震蕩,且誤差總體上并未隨時(shí)間積累,證明重力FKF-SCHA重力匹配算法能夠以相當(dāng)?shù)木葘?duì)慣導(dǎo)位置誤差進(jìn)行實(shí)時(shí)估計(jì)并修正。

(4) 結(jié)合試驗(yàn)區(qū)域重力特征參數(shù)與誤差曲線,當(dāng)局部粗糙度較大時(shí),匹配導(dǎo)航效果較優(yōu);反之,局部粗糙度較小時(shí),則匹配效果較差。當(dāng)移動(dòng)窗口范圍內(nèi)重力觀測(cè)量信息熵?cái)?shù)值在(0.4,0.6)范圍內(nèi)時(shí),匹配效果顯著;當(dāng)移動(dòng)窗口范圍內(nèi)重力觀測(cè)量信息熵?cái)?shù)值在(0,0.2]∪[0.8,1)時(shí),則匹配效果較差。

(5) 結(jié)合表6與圖6(d)可知,在24 h航行過(guò)程中,INS導(dǎo)航定位誤差最大值達(dá)到5.015 n mile,平均誤差為2.744 n mile。在使用FKF-SCHA重力匹配算法后導(dǎo)航定位誤差最大值為1.092 n mile,減少了78.2%;平均誤差為0.402 n mile,減少了85.2%;導(dǎo)航定位精度提高了85.3%。相較于FKF-SHA重力匹配算法,FKF-SCHA算法導(dǎo)航定位誤差最大值減少了14.4%,平均誤差減少了21.9%,導(dǎo)航定位精度提高了22.3%。統(tǒng)計(jì)結(jié)果表明,重力FKF-SCHA重力匹配算法可以有效減小INS累積誤差,提高匹配精度。

3.3 長(zhǎng)時(shí)間航行仿真試驗(yàn)

為驗(yàn)證FKF-SCHA重力匹配算法在潛航器長(zhǎng)時(shí)間水下航行時(shí)的有效性,在3.2節(jié)的仿真試驗(yàn)基礎(chǔ)上,進(jìn)行潛航器長(zhǎng)時(shí)間航行試驗(yàn)。試驗(yàn)航跡航行時(shí)長(zhǎng)240 h(10 d),初始速度與姿態(tài)分別為(0 m/s,0 m/s,0 m/s)和(0°,0°,0°),起始坐標(biāo)為(46°S,76°E,-100 m)。試驗(yàn)區(qū)域內(nèi)各個(gè)重力觀測(cè)量緯向和經(jīng)向信息熵均在(0.45,0.55)區(qū)間,表明該區(qū)域滿足重力匹配導(dǎo)航適配性需求。潛航器首先沿著正北方向以0.008 6 m/s2的加速度,勻加速至速度達(dá)到10 n mile/h,隨后勻速直線航行50 h;然后左轉(zhuǎn)90°,再勻速航行80 h;最后右轉(zhuǎn)90°,再勻速航行110 h。航跡與誤差曲線如圖7所示,匹配導(dǎo)航誤差統(tǒng)計(jì)結(jié)果見(jiàn)表7。

表7 10 d航行誤差統(tǒng)計(jì)

圖7 10 d航行航跡與誤差曲線

由表7與圖7可知,10 d航行過(guò)程中,慣性導(dǎo)航系統(tǒng)導(dǎo)航定位誤差最大達(dá)到46.074 n mile,平均值為16.739 n mile。由于地球自轉(zhuǎn)等周期性誤差影響,慣導(dǎo)經(jīng)向和緯向誤差均存在周期性變化,且經(jīng)向方向存在隨時(shí)間游走誤差。使用FKF-SCHA重力匹配算法后,導(dǎo)航定位誤差最大值降低為6.643 n mile,平均值為1.596 n mile,定位精度提高了88.7%。試驗(yàn)結(jié)果表明,FKF-SCHA重力匹配算法在潛航器長(zhǎng)時(shí)間無(wú)源導(dǎo)航時(shí)能夠有效抑制慣性導(dǎo)航系統(tǒng)隨時(shí)間發(fā)散的缺陷,提高系統(tǒng)整體導(dǎo)航定位精度。

4 結(jié) 論

本文針對(duì)傳統(tǒng)重力匹配導(dǎo)航算法各類觀測(cè)量的量測(cè)誤差會(huì)相互影響而導(dǎo)致濾波精度受損和局部重力場(chǎng)模型構(gòu)建不準(zhǔn)確引起濾波發(fā)散的問(wèn)題,提出了一種基于聯(lián)邦Kalman濾波和局部重力場(chǎng)球冠諧建模的水下重力匹配導(dǎo)航算法。首先通過(guò)FKF濾波器對(duì)多觀測(cè)量進(jìn)行分散降權(quán)處理,同時(shí)基于適定邊值問(wèn)題構(gòu)建FKF子系統(tǒng)濾波器的重力梯度復(fù)數(shù)組合觀測(cè)量;然后基于SCHA理論,采用ASHA技術(shù),快速建立以匹配點(diǎn)為中心的移動(dòng)窗口球冠諧模型,依據(jù)該模型將重力量測(cè)值表示為連續(xù)的解析形式并組建子濾波器非線性量測(cè)方程;最后利用預(yù)測(cè)殘差向量設(shè)計(jì)自適應(yīng)信息分配因子將各子濾波器狀態(tài)估值及協(xié)方差進(jìn)行融合得到慣導(dǎo)位置誤差最優(yōu)估計(jì)量。試驗(yàn)表明:采用SCHA建模的重力FKF匹配算法在24 h航行過(guò)程中導(dǎo)航定位誤差保持在1.1 n mile以內(nèi),導(dǎo)航定位精度提高了85%以上,并根據(jù)水下潛航器長(zhǎng)期水下作業(yè)航行特性,設(shè)計(jì)了一組長(zhǎng)時(shí)間水下無(wú)源定位導(dǎo)航試驗(yàn),試驗(yàn)航行時(shí)長(zhǎng)10 d,重力FKF-SCHA匹配算法導(dǎo)航定位精度提高了88.7%。試驗(yàn)結(jié)果分析表明,重力FKF-SCHA匹配算法能夠在一定程度上克服慣性導(dǎo)航系統(tǒng)由于時(shí)間推移誤差積累的缺陷,提高系統(tǒng)導(dǎo)航定位精度,增加匹配算法的穩(wěn)健性。