基于Kalman 濾波的船舶磁化干擾系數(shù)測(cè)量算法

閆輝，周國(guó)華

1 信陽(yáng)學(xué)院 大數(shù)據(jù)與人工智能學(xué)院，河南 信陽(yáng) 464000

2 海軍工程大學(xué) 電氣工程學(xué)院，湖北 武漢 430033

0 引 言

置于地球磁場(chǎng)中的鐵磁物體，例如艦艇、坦克、車輛、飛機(jī)等，都會(huì)因地磁場(chǎng)作用而產(chǎn)生磁化磁場(chǎng)[1]，所以安裝在船舶上的磁傳感器會(huì)受到船舶磁化磁場(chǎng)的干擾。如何從磁傳感器的測(cè)量數(shù)據(jù)中分離出三分量地磁場(chǎng)，是以船舶為載體的地磁場(chǎng)測(cè)量系統(tǒng)所面臨的核心問題。艦載消磁系統(tǒng)是一種以三分量地磁場(chǎng)作為系統(tǒng)控制信號(hào)[2]來補(bǔ)償船舶磁化磁場(chǎng)的船用設(shè)備，該系統(tǒng)采用安裝于艦艇桅桿上經(jīng)過抗干擾調(diào)整的三分量磁傳感器來實(shí)時(shí)監(jiān)測(cè)地磁場(chǎng)，進(jìn)而控制消磁繞組電流以補(bǔ)償艦船磁場(chǎng)，即“測(cè)地磁消艦磁”，其中抗干擾調(diào)整的目的是消除船舶磁化磁場(chǎng)對(duì)地磁測(cè)量傳感器的干擾。

為了消除船舶磁化磁場(chǎng)對(duì)地磁測(cè)量傳感器的干擾，需預(yù)先測(cè)量船舶的感應(yīng)干擾系數(shù)矩陣和固定干擾系數(shù)向量[3]。根據(jù)以船舶為載體的三分量地磁場(chǎng)測(cè)量系統(tǒng)數(shù)學(xué)模型，肖昌漢[4]提出了“四航向法”，但該方法只能測(cè)量感應(yīng)干擾系數(shù)矩陣的前2 列元素，且需在不同的兩地進(jìn)行補(bǔ)充測(cè)量。肖昌漢等[5]提出了“8”字航向法，即要求船舶在已知地磁場(chǎng)的測(cè)量場(chǎng)地進(jìn)行高速、大舵角的變航向機(jī)動(dòng)，進(jìn)而造成較大幅度的船舶姿態(tài)變化，并需要大量的有效磁場(chǎng)測(cè)量數(shù)據(jù)作為支撐。閆輝等[6]對(duì)船載三分量地磁場(chǎng)的測(cè)量過程進(jìn)行了實(shí)驗(yàn)?zāi)M，發(fā)現(xiàn)即使在實(shí)驗(yàn)室條件下，船舶磁化干擾系數(shù)的測(cè)量也存在一定困難。因此，船舶磁化磁場(chǎng)干擾系數(shù)的測(cè)量問題是限制船載三分量地磁場(chǎng)測(cè)量系統(tǒng)大規(guī)模工程應(yīng)用的難點(diǎn)。

為了降低測(cè)量船舶磁化磁場(chǎng)干擾系數(shù)的工程難度和成本，本文擬基于船舶為載體的三分量地磁場(chǎng)測(cè)量系統(tǒng)數(shù)學(xué)模型，分析解決船舶磁化磁場(chǎng)干擾的難點(diǎn)問題。首先，將根據(jù)三分量地磁場(chǎng)解算的Kalman 濾波方法[3]，重新建立適用于船舶磁化干擾系數(shù)測(cè)量的狀態(tài)空間方程，提出采用Kalman濾波方法求解干擾系數(shù)的自適應(yīng)算法，并通過計(jì)算機(jī)仿真以驗(yàn)證該算法的理論正確性；然后，將開展實(shí)驗(yàn)室模擬測(cè)量，在有限的測(cè)量數(shù)據(jù)樣本條件下，應(yīng)用解算所得的干擾系數(shù)來有效排除船舶磁化磁場(chǎng)的干擾，從而分離出較高精度的三分量地磁場(chǎng)，用以為船舶磁化干擾系數(shù)測(cè)量提供一條高效率、低成本的可行路徑。

1 三分量地磁場(chǎng)測(cè)量的數(shù)學(xué)模型

三分量磁傳感器作為磁場(chǎng)的敏感元件，可以采集其所處位置的各種矢量合成磁場(chǎng)，其中既包括地磁場(chǎng)，也包括來自測(cè)量船自身的磁化磁場(chǎng)干擾。安裝于船舶桅桿的三分量磁傳感器的動(dòng)態(tài)輸出3×1 向量B為[4]

式 中：K為3×3 感 應(yīng) 磁 化 干 擾 系 數(shù) 矩 陣[7]；A為3×3 地理坐標(biāo)系與船舶坐標(biāo)系之間的轉(zhuǎn)換矩陣[8-9]，該矩陣元素為測(cè)量船姿態(tài)角的函數(shù)；Be為3×1 三分量地磁場(chǎng)向量，nT；Bp為3×1 固定磁化干擾系數(shù)向量[10]，nT。其中，

式中，φ，θ，ψ 分別為船舶的航向角、橫搖角和縱傾角，可以由船舶的導(dǎo)航系統(tǒng)進(jìn)行實(shí)時(shí)測(cè)量。

根據(jù)式（1），如果需從傳感器輸出B的測(cè)量值中分離出地磁場(chǎng)Be，則必須先測(cè)量船舶磁化干擾系數(shù)K和Bp。將式（1）展開并整理，可知在不同的艦船姿態(tài)下，傳感器的任意輸出分量為

式中：Bi為磁傳感器第i軸向的輸出，其中i=1，2，3，分別對(duì)應(yīng)x軸、y軸、z軸；ki1，ki2，ki3分別為矩陣K中第i行的第1，2，3 列元素；f1，f2，f3分別為式（1）中的函數(shù)映射關(guān)系；Bpi為3×1 固定磁化干擾系數(shù)向量Bp在第i軸向的分量。

在已知地磁場(chǎng)和船舶姿態(tài)的條件下，需至少記錄4 組不同姿態(tài)下的傳感器輸出，才能通過求解一個(gè)線性方程組來得到船舶的磁化干擾系數(shù)。將該系數(shù)代入式(1)，則未知量地磁場(chǎng)Be的解為

上述地磁場(chǎng)的解算過程雖然在理論上成立，但在實(shí)際工程運(yùn)用中卻難以保證解算精度，主要原因如下：1）三分量磁通門傳感器普遍存在固有的零位誤差、標(biāo)度誤差和正交誤差[11]；2）船舶姿態(tài)角φ，θ，ψ 數(shù)據(jù)主要來源于導(dǎo)航系統(tǒng)，存在一定的漂移誤差，以及信號(hào)傳輸過程中的延遲誤差；3）由于背景環(huán)境的磁干擾及地磁場(chǎng)靜日變化，將導(dǎo)致一定的磁傳感器隨機(jī)誤差[12]。由于船舶磁化干擾系數(shù)對(duì)以上各種誤差源較敏感，所以這些誤差將在測(cè)量過程中累積，進(jìn)而通過式（4）向地磁場(chǎng)解算值傳遞。目前，除了第1 種誤差之外，其他2 類誤差暫無較好的處理方法。因此，為了提高三分量地磁場(chǎng)的解算精度，需采用相應(yīng)措施來消除船舶磁化干擾系數(shù)的測(cè)量誤差，從而提高式（1）數(shù)學(xué)模型的準(zhǔn)確性。廣泛應(yīng)用于動(dòng)態(tài)測(cè)量系統(tǒng)的Kalman 濾波方法在原理上可以較好地抑制系統(tǒng)測(cè)量誤差[3,13]，且可以通過信息的引入對(duì)變化狀態(tài)進(jìn)行及時(shí)更新。將Kalman 濾波應(yīng)用于磁化干擾系數(shù)測(cè)量時(shí)，必須以磁化干擾系數(shù)作為狀態(tài)量，但該狀態(tài)具有相對(duì)不變性[4]，因此磁化干擾系數(shù)測(cè)量具有靜態(tài)測(cè)量系統(tǒng)的特點(diǎn)。如果將靜態(tài)系統(tǒng)視為動(dòng)態(tài)系統(tǒng)的一種特例，則Kalman 濾波可為解決船舶磁化干擾系數(shù)測(cè)量問題提供一條新途徑。

2 Kalman 濾波算法

2.1 算法原理

Kalman 濾波算法需以系統(tǒng)的時(shí)間離散化狀態(tài)空間為基礎(chǔ)[14]，在船舶結(jié)構(gòu)及磁傳感器安裝位置不變的條件下，船舶磁化系數(shù)K和Bp均為時(shí)域不變量，當(dāng)以4×3 的擴(kuò)展矩陣x=[KBp]T（T 表示矩陣的轉(zhuǎn)置）為狀態(tài)量時(shí)，則測(cè)量過程的計(jì)算方程為

式中：x(n)為狀態(tài)量，測(cè)量點(diǎn)數(shù)n=0，1，2，···，N，其中N為最大測(cè)量點(diǎn)數(shù)；4×3 向量v1(n)為過程噪聲向量，用于描述狀態(tài)轉(zhuǎn)移誤差。

根據(jù)式(1)，即可建立觀測(cè)方程：

式中：3×1 向量y(n)=BT，1×4 向量C(n)=[AT1]，3×1 向量v2(n)為觀測(cè)噪聲向量，用于描述隨機(jī)噪聲。

由式(5)和式(6)構(gòu)成狀態(tài)空間，則自適應(yīng)Kalman 濾波迭代算法的步驟如下：

1）設(shè)置初始條件：狀態(tài)x(0)，中間變量P(0)（4×4 對(duì)角矩陣）。

2）已知參數(shù)：4×4 過程噪聲方差陣Q，一維觀測(cè)噪聲方差σ2。

3）輸入：y(n)，C(n)。

4）計(jì)算中間量及狀態(tài)量，其中：R(n)，P(n)，m(n)均為中間變量；E為單位矩陣。

通過以上算法，即可得到狀態(tài)量x(n)T，前3 列可以構(gòu)成感應(yīng)干擾系數(shù)矩陣K，第4 列則可以構(gòu)成固定干擾系數(shù)向量Bp。在上述迭代過程中，初始條件的設(shè)置具有一定的任意性，對(duì)于狀態(tài)先驗(yàn)信息不足的條件下，初始狀態(tài)x(0)取零值，中間變量初始值P(0)取較小值將有利于迭代過程的收斂。作為已知參數(shù)，過程噪聲方差陣Q和觀測(cè)噪聲σ2的取值對(duì)迭代過程的收斂速度和收斂穩(wěn)定性具有決定性作用：當(dāng)這2 個(gè)參數(shù)取較大值時(shí)，收斂速度快，但穩(wěn)定性較差；反之，則收斂速度慢，而穩(wěn)定性較好，因此，本文將通過模擬仿真試算的方法來平衡確定參數(shù)值。

2.2 算法仿真

為了驗(yàn)證2.1 節(jié)算法的正確性及其對(duì)隨機(jī)誤差的抑制能力，本節(jié)將開展仿真試算。設(shè)定三分量地磁場(chǎng)向量為Be=[34 000, 2 000, 35 000] nT，感應(yīng)磁化干擾系數(shù)矩陣K=[k11，k12，k13;k21，k22，k23;k31，k32，k33]=[0.95, -0.003, 0.003; 0.003, 0.95, -0.03;0.03, -0.03, 1.05]，固定磁化干擾系數(shù)向量Bp=[1 500, -800, 1 800] nT。船舶姿態(tài)角φ，θ，ψ 分別在區(qū)間[0°，360°]，[-20°，20°]，[-15°，15°]取值，并生成轉(zhuǎn)換矩陣A。按照式(1)生成12 組傳感器輸出向量B，并施加最大值為100 nT 的隨機(jī)誤差。記錄每一組向量B所對(duì)應(yīng)的φ，θ，ψ，并施加最大值為0.05°的隨機(jī)誤差。將數(shù)據(jù)進(jìn)行延長(zhǎng)處理，即可重復(fù)使用多次，從而保證算法具有足夠數(shù)量的迭代步數(shù)。此時(shí)，以Be，A，B為已知量，即可按照2.1 節(jié)的Kalman 濾波算法解算K和Bp，算法收斂情況如圖1 所示。

由圖1 可知，該Kalman 濾波算法經(jīng)過約300 次迭代即可實(shí)現(xiàn)收斂，且具有良好的穩(wěn)定性。由于本文在仿真過程施加了隨機(jī)誤差，其解算結(jié)果為：

圖1 磁化干擾系數(shù)的仿真收斂曲線Fig. 1 Convergence curves of magnetization factors in simulation

與設(shè)定值相比，K的誤差出現(xiàn)在萬分位，Bp的誤差出現(xiàn)在個(gè)位。鑒于地磁場(chǎng)為104nT 量級(jí)，按式（3）的傳遞關(guān)系估算誤差，可知因船舶磁化干擾系數(shù)測(cè)量誤差所造成的地磁場(chǎng)解算誤差約為10 nT 量級(jí)，該誤差水平與文獻(xiàn)[3]中僅考慮隨機(jī)測(cè)量誤差所得的地磁場(chǎng)解算誤差水平相當(dāng)。

3 船模驗(yàn)證實(shí)驗(yàn)

3.1 實(shí)驗(yàn)?zāi)Ｐ图皽y(cè)量過程

為了檢驗(yàn)船舶磁化干擾系數(shù)測(cè)量時(shí)Kalman濾波方法的有效性，本節(jié)將在無磁實(shí)驗(yàn)室開展某型船模的驗(yàn)證實(shí)驗(yàn)。測(cè)量裝置由三分量磁通門傳感器、水平角度傳感器及硬質(zhì)無磁不銹鋼支架組成，其中磁傳感器與角度傳感器通過支架進(jìn)行剛性連接，需注意調(diào)整相應(yīng)測(cè)量軸。將該測(cè)量裝置安裝于某型磁性船模上，測(cè)量軸與艦船坐標(biāo)系保持一致，如圖2 所示。將磁性船模放置于無磁平臺(tái)，該平臺(tái)可以滿足船模4 個(gè)主航向旋轉(zhuǎn)、橫搖和縱傾要求，據(jù)此即可構(gòu)建船載地磁三分量測(cè)量平臺(tái)。

圖2 實(shí)驗(yàn)裝置示意圖Fig. 2 Experiment device

為了簡(jiǎn)化測(cè)量過程，船模分別在磁東、磁北、磁西、磁南4 個(gè)主航向進(jìn)行了3 種水平姿態(tài)下的磁干擾測(cè)量，共計(jì)獲得12 組測(cè)量數(shù)據(jù)，并同步記錄了三分量磁傳感器和水平角度傳感器的測(cè)量數(shù)據(jù)。

3.2 測(cè)量結(jié)果及計(jì)算分析

測(cè)量數(shù)據(jù)如表1 所示，其中磁傳感器三分量、橫搖角、縱傾角分別為傳感器實(shí)測(cè)值，磁東、磁北、磁西、磁南4 個(gè)航向角則直接取值90°，0°，270°，180°，并轉(zhuǎn)換為弧度。

已知環(huán)境地磁場(chǎng)為Be=[34 425, 1961, 35 898] nT，采用上文Kalman 濾波算法，即可得到：

固定干擾系數(shù)向量為

由于測(cè)量值的數(shù)據(jù)量過?。▋H12 組），導(dǎo)致濾波器的迭代步數(shù)過少，因此無法收斂。為了解決這一問題，本文將這12 組數(shù)據(jù)進(jìn)行了延長(zhǎng)處理，即進(jìn)行12 組數(shù)據(jù)的多次重復(fù)使用，并以此作為實(shí)驗(yàn)數(shù)據(jù)集。從實(shí)驗(yàn)數(shù)據(jù)集中抽取一定數(shù)量的數(shù)據(jù)，即可解算磁化干擾系數(shù)，圖3 所示為感應(yīng)干擾系數(shù)和固定干擾系數(shù)的迭代過程曲線。

由圖3 可知，對(duì)于測(cè)量數(shù)據(jù)較少的情況，通過對(duì)數(shù)據(jù)進(jìn)行重復(fù)性延長(zhǎng)處理，Kalman 濾波算法仍然具有較好的收斂性；但與仿真結(jié)果相比，其收斂之后的穩(wěn)定性略有下降，存在微幅震蕩。由于仿真所采用的數(shù)據(jù)均為無誤差的理想值，而實(shí)驗(yàn)中傳感器所測(cè)量的磁場(chǎng)、航向角、橫搖角、縱傾角均存在誤差，因此，造成穩(wěn)定性下降的主要原因是實(shí)驗(yàn)數(shù)據(jù)的測(cè)量誤差。

圖3 干擾系數(shù)的實(shí)驗(yàn)收斂曲線Fig. 3 Convergence curves of magnetization factors in experiment

為了驗(yàn)證表1 船舶磁化干擾系數(shù)的正確性，將其代入文獻(xiàn)[3]的三分量地磁場(chǎng)解算算法，并從實(shí)驗(yàn)數(shù)據(jù)集中抽取一定樣本的數(shù)據(jù)用于解算地磁場(chǎng)，其收斂過程如圖4 所示，其中Bei為3×1 三分量地磁場(chǎng)向量Be在第i軸向的分量（i=1，2，3）。根據(jù)計(jì)算結(jié)果，收斂穩(wěn)定后的地磁場(chǎng)解算數(shù)據(jù)的最大相對(duì)誤差為1.9%（此處定義為：誤差最大值與地磁場(chǎng)總量的比值），均方根誤差為266 nT。由此可見，在僅有12 組測(cè)量數(shù)據(jù)的條件下，Kalman濾波算法具有較高的計(jì)算精度，其對(duì)測(cè)量誤差具有較好的抑制能力。

圖4 三分量地磁場(chǎng)計(jì)算收斂曲線Fig. 4 The convergence curves in computing three-component geomagnetic field

表1 船模磁場(chǎng)測(cè)量數(shù)據(jù)表Table 1 Measured magnetic field data of the mockup

4 結(jié) 語(yǔ)

對(duì)于以船舶為載體的三分量地磁場(chǎng)測(cè)量系統(tǒng)，其核心問題是船舶磁化干擾系數(shù)的測(cè)量。本文提出了基于Kalman 濾波的船舶磁化干擾系測(cè)量算法，根據(jù)安裝于船舶的三分量磁場(chǎng)測(cè)量系統(tǒng)的傳感器測(cè)量輸出數(shù)學(xué)模型，建立了觀測(cè)方程；根據(jù)船舶干擾系數(shù)的相對(duì)不變性，建立了過程方程。由觀測(cè)方程和過程方程所構(gòu)成的狀態(tài)空間表達(dá)式較好地契合了Kalman 濾波對(duì)于模型準(zhǔn)確性的嚴(yán)苛要求，并通過計(jì)算機(jī)仿真驗(yàn)證了該算法理論的正確性。在較少實(shí)驗(yàn)數(shù)據(jù)樣本的情況下，僅通過數(shù)據(jù)的簡(jiǎn)單重復(fù)利用，即可保證Kalman 濾波算法的良好收斂性。相較于“四航向法”，Kalman濾波方法可以測(cè)量所有磁化干擾系數(shù)，且無需嚴(yán)格按照東、西、南、北4 個(gè)主航向?qū)嵤y(cè)量；相較于“8”字航向法，該方法的數(shù)據(jù)量小，且無需按某種特定航跡實(shí)施測(cè)量。因此，從干擾系數(shù)測(cè)量的完備性，所需數(shù)據(jù)量的大小，以及測(cè)量實(shí)施條件的難易程度而言，Kalman 濾波方法可以顯著提高船舶磁化干擾系數(shù)的測(cè)量效率，并降低工程實(shí)施的難度。

在實(shí)驗(yàn)室模擬測(cè)量過程中，影響測(cè)量精度的主要原因?yàn)椋喝至看艂鞲衅鳑]有經(jīng)過固有誤差校正；測(cè)量干擾系數(shù)時(shí)，環(huán)境地磁場(chǎng)存在測(cè)量誤差；傳感器存在安裝誤差，尤其是傳感器x軸向與模型艏線、艉線方向的一致性誤差；角度傳感器僅能測(cè)量橫搖角和縱傾角，航向角僅采用了實(shí)驗(yàn)室軌道方向作為近似值。因此，在實(shí)際工程應(yīng)用中，還需進(jìn)一步處理這些誤差，以提高船舶感應(yīng)干擾系數(shù)的測(cè)量精度。