靳宇航，王海涌，劉 濤，賈平會(huì)，王永海

(1.北京航空航天大學(xué)宇航學(xué)院，北京100191；2.北京航天萬(wàn)源科技有限公司，北京100176；3.北京航天長(zhǎng)征飛行器研究所，北京100076)

0 引言

在導(dǎo)彈的導(dǎo)航系統(tǒng)中,慣性/天文/GPS三者或兩兩之間的組合導(dǎo)航是目前組合導(dǎo)航的主要研究范疇。其中,慣性/天文組合導(dǎo)航系統(tǒng)、慣性/GPS組合導(dǎo)航系統(tǒng)已經(jīng)得到應(yīng)用,但仍有一定的局限性。星敏感器成本高、體積大、價(jià)格昂貴,應(yīng)用時(shí)需要主動(dòng)控制以保障姿態(tài)穩(wěn)定,避免成像拖尾。GPS及北斗導(dǎo)航信號(hào)易受干擾,在實(shí)戰(zhàn)中難以使用。地磁導(dǎo)航具有自主性,長(zhǎng)期工作無(wú)積累誤差,隨著器件和算法的發(fā)展,地磁導(dǎo)航越來(lái)越引起了人們的關(guān)注。因此,某些型號(hào)應(yīng)用自主導(dǎo)航、轉(zhuǎn)向慣性和地磁的組合尋求技術(shù)解決方案。地磁場(chǎng)的磁強(qiáng)范圍在50000nT～60000nT之間,目前已有的三軸磁強(qiáng)計(jì)分辨率可達(dá)0.1nT,MEMS磁強(qiáng)計(jì)的分辨率可以達(dá)到10nT左右,可以用于地磁輔助導(dǎo)航。慣性/地磁組合導(dǎo)航主要有基于批處理算法的地磁匹配輔助導(dǎo)航和地磁序貫濾波輔助導(dǎo)航兩種方式。前者是將地磁匹配定位算法得到的載體位置估計(jì)作為量測(cè)輸入,再利用濾波技術(shù)估計(jì)慣導(dǎo)系統(tǒng)的誤差,該算法需要存儲(chǔ)大量的地磁數(shù)據(jù),并且累積一段路徑才能完成定位,實(shí)時(shí)性較差；后者直接利用磁強(qiáng)計(jì)的測(cè)量值作為濾波器的量測(cè)輸入,需要的存儲(chǔ)空間較小,同時(shí)結(jié)合遞推濾波技術(shù),直接把地磁測(cè)量值作為觀測(cè)量,每獲得一個(gè)測(cè)量值便對(duì)系統(tǒng)進(jìn)行實(shí)時(shí)更新,具有連續(xù)修正能力,導(dǎo)彈和衛(wèi)星通常采用后者。

雷達(dá)高度表(Radar Altimeter,RA)作為一種主動(dòng)式距離傳感器,與氣壓高度表和激光高度表相比,雷達(dá)高度表具有測(cè)高范圍大、精度高、功耗低、體積小等優(yōu)點(diǎn)。雷達(dá)測(cè)距技術(shù)已成功應(yīng)用于軌道測(cè)量、彈道修正、衛(wèi)星定位、交會(huì)對(duì)接、巡航制導(dǎo)、航天著陸、大地測(cè)量、微波遙感等方面。本文利用雷達(dá)高度表的測(cè)量信息,建立了一種估計(jì)導(dǎo)彈位置的方法,結(jié)合磁強(qiáng)計(jì)的輸出進(jìn)行姿態(tài)矯正,用一個(gè)觀測(cè)式就可以包含姿態(tài)和位置信息,實(shí)現(xiàn)對(duì)導(dǎo)彈的組合導(dǎo)航。

1 地磁場(chǎng)模型

地磁場(chǎng)是一個(gè)矢量場(chǎng),它是時(shí)間和位置的函數(shù),但是地磁場(chǎng)隨著時(shí)間變化緩慢,目前普遍采用的全球地磁模型有國(guó)際地磁參考場(chǎng)模型(IGRF)和世界地磁模型(WMM)。本文取用IGRF作為研究的背景場(chǎng),地磁場(chǎng)勢(shì)函數(shù)V表示為:

式中,Re表示地球半徑,r表示載體的地心距,λ表示地理經(jīng)度,θ代表了地理緯度,gmn(t)和hmn(t)是Gauss系數(shù),IGRF模型提供相應(yīng)的Gauss系數(shù)表,Pmn(cosθ)為Schmidt準(zhǔn)歸一化n次m階Legendre函數(shù)。

地磁場(chǎng)強(qiáng)度B在地磁球面坐標(biāo)系f與地球坐標(biāo)系e中的轉(zhuǎn)換關(guān)系如圖1所示。

從而得到地磁場(chǎng)矢量B在地球坐標(biāo)系e下的矢量Be,表示為:

2 捷聯(lián)慣導(dǎo)/地磁/高度表組合系統(tǒng)模型

本文涉及的坐標(biāo)系有:導(dǎo)彈發(fā)射點(diǎn)慣性坐標(biāo)系l系、地球坐標(biāo)系e系、地心慣性系i系、本體坐標(biāo)系b系。 文中以發(fā)射點(diǎn)慣性系l系為導(dǎo)航坐標(biāo)系,建立捷聯(lián)慣導(dǎo)/地磁/雷達(dá)高度表組合導(dǎo)航系統(tǒng)的狀態(tài)方程和量測(cè)方程。

2.1 組合導(dǎo)航系統(tǒng)狀態(tài)方程

系統(tǒng)狀態(tài)模型為慣性導(dǎo)航系統(tǒng)的誤差模型,慣性器件加速度計(jì)、陀螺建模零偏和白噪聲的疊加,建立15維狀態(tài)方程:

其中,狀態(tài)變量為:

其中,?x、?y、?z為捷聯(lián)慣導(dǎo)數(shù)學(xué)平臺(tái)失準(zhǔn)角,δvx、δvy、δvz為發(fā)射點(diǎn)慣性系3個(gè)軸上的速度誤差, δx、 δy、 δz為位置誤差,εx、εy、εz為陀螺儀常值漂移,為加速度計(jì)常值偏置。系統(tǒng)噪聲為:

式中,各項(xiàng)分別為三軸陀螺儀誤差和三軸加速度計(jì)誤差,均建模為Guass白噪聲。

2.2 雷達(dá)高度表量測(cè)模型

導(dǎo)彈在飛行過(guò)程中,通過(guò)向地球反射表面(海平面或地表平面)發(fā)射一系列雷達(dá)脈沖,由彈上雷達(dá)信號(hào)接收機(jī)測(cè)得回波信號(hào)后,可根據(jù)時(shí)間間隔Δt得到雷達(dá)高度表測(cè)量高度H′為:

其中,c為電磁波傳播速度。為了獲得雷達(dá)高度表測(cè)得的海拔觀測(cè)高度Ho,需在實(shí)測(cè)值H′的基礎(chǔ)上進(jìn)行高度補(bǔ)償,即:

補(bǔ)償值Δh由當(dāng)?shù)氐乇砗０蔚纫蛩貨Q定。

接下來(lái)推導(dǎo)高度表的輸出和捷聯(lián)慣導(dǎo)解算高度的關(guān)系,導(dǎo)彈的高度示意圖如圖2所示。其中,Re為當(dāng)前所在位置的地球半徑,Re0為發(fā)射點(diǎn)至地心距離。

在l系Ol－XlYlZl下,位置矢量OlP與發(fā)射點(diǎn)水平面夾角θ0滿(mǎn)足:

式中,x、y、z分別為慣導(dǎo)位置解算結(jié)果在l系下的坐標(biāo),則該時(shí)刻慣導(dǎo)測(cè)得的海拔計(jì)算高度Hc為:

Hc與Ho之差ΔH滿(mǎn)足如下關(guān)系:

其中,

wh為誤差項(xiàng),受導(dǎo)航誤差及量測(cè)噪聲影響。δx、δy、δz分別為捷聯(lián)解算與標(biāo)稱(chēng)彈道的位置誤差。

2.3 地磁量測(cè)模型

(1)姿態(tài)誤差角與平臺(tái)失準(zhǔn)角的關(guān)系

在捷聯(lián)慣性導(dǎo)航系統(tǒng)中,載體姿態(tài)角是通過(guò)姿態(tài)矩陣計(jì)算出來(lái)的。理想情況下,導(dǎo)航計(jì)算機(jī)計(jì)算的導(dǎo)航坐標(biāo)系(?系)應(yīng)和理想的導(dǎo)航坐標(biāo)系(l系)一致,然而,實(shí)際計(jì)算的導(dǎo)航系與理想的導(dǎo)航系之間將產(chǎn)生偏差,對(duì)應(yīng)的誤差角為平臺(tái)失準(zhǔn)角,用角矢量表示。 在小角度假設(shè)下,l系與?系之間滿(mǎn)足如下關(guān)系:

式中,[φ×]為平臺(tái)失準(zhǔn)角φ的反對(duì)稱(chēng)矩陣。定義姿態(tài)誤差角矢量η為:

式中,Cb?為捷聯(lián)解算的姿態(tài)矩陣。忽略2階小量,可得l系下姿態(tài)誤差角與平臺(tái)失準(zhǔn)角關(guān)系:

式中,M為誤差角轉(zhuǎn)換矩陣:

(2)地磁觀測(cè)模型

三軸磁強(qiáng)計(jì)測(cè)量值Bm可表示為:

磁強(qiáng)計(jì)在載體預(yù)估狀態(tài)下,根據(jù)地磁場(chǎng)模型計(jì)算的地磁場(chǎng)強(qiáng)度值Bc可表示為:

定義上述磁強(qiáng)計(jì)測(cè)量值Bm與根據(jù)地磁場(chǎng)模型計(jì)算的地磁場(chǎng)強(qiáng)度值Bc之差為量測(cè)Z1:

地磁場(chǎng)矢量Bb與其估計(jì)值之間滿(mǎn)足如下數(shù)學(xué)關(guān)系:

將式(22)和式(23)代入式(21)中,可得:

由式(25)可知,磁強(qiáng)計(jì)測(cè)量誤差主要由兩部分組成:姿態(tài)誤差引起的測(cè)量誤差和位置誤差引起的測(cè)量誤差,本文重點(diǎn)研究姿態(tài)部分引起的誤差。

式(25)右邊第一項(xiàng)可以簡(jiǎn)寫(xiě)成:

根據(jù)微分運(yùn)算法則:

式中,Bl為地磁場(chǎng)矢量在l系下的表示。定義為:

得到:

由于組合系統(tǒng)狀態(tài)變量中代表姿態(tài)信息的是平臺(tái)失準(zhǔn)角φ,而式(31)中的量測(cè)是用姿態(tài)誤差角η表示的,將前面推導(dǎo)的它們的關(guān)系代入式(31)中,得到:

因此,地磁量測(cè)模型為:

其中,量測(cè)矩陣H1滿(mǎn)足:

2.4 系統(tǒng)觀測(cè)方程

在推導(dǎo)了雷達(dá)高度表和磁強(qiáng)計(jì)的觀測(cè)模型之后,得到如下的系統(tǒng)觀測(cè)方程:

式中,[η×]為姿態(tài)誤差角η的元素組成的反對(duì)稱(chēng)矩陣。將式(28)和式(29)代入式(26),可得:

Z(t)包含地磁誤差Z1和高度誤差Z2兩部分,。 其中,H1的 表達(dá)如式(34)所示,

V為觀測(cè)噪聲,仿真設(shè)定參考實(shí)際器件誤差。

2.5 濾波器設(shè)計(jì)

本文推導(dǎo)得到的測(cè)量方程是非線(xiàn)性方程,采用工程上常用的擴(kuò)展Kalman濾波方法,并采用反饋矯正濾波器。圖3為組合導(dǎo)航系統(tǒng)原理圖。H2表達(dá)式如下所示:

擴(kuò)展Kalman濾波方程為以下形式:

3 仿真分析

3.1 仿真條件

1)導(dǎo)彈發(fā)射點(diǎn)位于北緯 39.98°,東經(jīng)116.34°,采用垂直向東發(fā)射,主動(dòng)段飛行時(shí)間為160s,全程飛行時(shí)間為1000s。

2)初始俯仰角為90°,偏航角為0°,橫滾角為0°；俯仰角誤差為3′,偏航角誤差為6′,橫滾角誤差為 3′。

3)采樣周期:陀螺及加速度計(jì)為0.01s,磁強(qiáng)計(jì)和高度表為0.1s,Kalman濾波周期為0.1s。

4)元器件參數(shù):陀螺常值漂移為0.003(°)/s,陀螺測(cè)量白噪聲標(biāo)準(zhǔn)差為0.0003(°)/s；加速度計(jì)常值偏置為1mg,加速度計(jì)白噪聲標(biāo)準(zhǔn)差為0.1mg；磁強(qiáng)計(jì)測(cè)量誤差為100nT；雷達(dá)高度表測(cè)量誤差為50m。

3.2 仿真結(jié)果及分析

捷聯(lián)慣性導(dǎo)航及與磁強(qiáng)計(jì)高度表組合導(dǎo)航相應(yīng)的位置誤差、速度誤差、姿態(tài)角誤差如圖4～圖8所示。經(jīng)過(guò)1000s的飛行,純捷聯(lián)解算的位置誤差為44.58km,速度誤差為52.40m/s,俯仰角誤差為3.21°,偏航角誤差為3.07°,橫滾角誤差為2.94°。利用磁強(qiáng)計(jì)和高度表進(jìn)行組合導(dǎo)航之后,上述誤差的發(fā)散均得到有效抑制,位置誤差為2.68km,速度誤差為8.4m/s,俯仰角誤差為0.31°,偏航角誤差為0.25°,橫滾角誤差為0.19°。

由仿真結(jié)果可知,加上磁強(qiáng)計(jì)和高度表之后進(jìn)行組合,導(dǎo)航精度得到明顯的提高。因?yàn)闇y(cè)量方程中既包含了姿態(tài)誤差,又包含了位置誤差,引入了磁強(qiáng)計(jì)和高度表的測(cè)量值,采用濾波估計(jì)的方法能夠較好地估計(jì)系統(tǒng)的誤差。

4 結(jié)論

本文基于雷達(dá)高度表和磁強(qiáng)計(jì)的測(cè)量信息,提出一種導(dǎo)彈捷聯(lián)慣性/地磁/高度表組合導(dǎo)航方法。以捷聯(lián)慣導(dǎo)誤差方程為基礎(chǔ)建立系統(tǒng)的狀態(tài)模型,以磁強(qiáng)計(jì)測(cè)量值、地磁場(chǎng)模型計(jì)算的地磁場(chǎng)強(qiáng)度之差和高度表測(cè)量與捷聯(lián)慣導(dǎo)解算高度之差共同作為量測(cè),只用一個(gè)觀測(cè)表達(dá)式即同時(shí)包含載體的位置和姿態(tài)信息。引入狀態(tài)反饋,通過(guò)混合校正的Kalman濾波進(jìn)行信息融合,得到系統(tǒng)導(dǎo)航信息的最優(yōu)估計(jì)。仿真結(jié)果表明,該算法能有效抑制捷聯(lián)解算誤差的發(fā)散。值得一提的是,慣性器件和磁強(qiáng)計(jì)可以選用低成本、中等精度的MEMS器件,經(jīng)過(guò)組合導(dǎo)航之后可以到達(dá)較高的精度,該方法完全自主,具有一定工程應(yīng)用價(jià)值。