基于Kalman 濾波的船舶磁化干擾系數測量算法

閆輝，周國華

1 信陽學院 大數據與人工智能學院，河南 信陽 464000

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

0 引 言

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

為了消除船舶磁化磁場對地磁測量傳感器的干擾，需預先測量船舶的感應干擾系數矩陣和固定干擾系數向量[3]。根據以船舶為載體的三分量地磁場測量系統數學模型，肖昌漢[4]提出了“四航向法”，但該方法只能測量感應干擾系數矩陣的前2 列元素，且需在不同的兩地進行補充測量。肖昌漢等[5]提出了“8”字航向法，即要求船舶在已知地磁場的測量場地進行高速、大舵角的變航向機動，進而造成較大幅度的船舶姿態變化，并需要大量的有效磁場測量數據作為支撐。閆輝等[6]對船載三分量地磁場的測量過程進行了實驗模擬，發現即使在實驗室條件下，船舶磁化干擾系數的測量也存在一定困難。因此，船舶磁化磁場干擾系數的測量問題是限制船載三分量地磁場測量系統大規模工程應用的難點。

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

1 三分量地磁場測量的數學模型

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

式 中：K為3×3 感 應 磁 化 干 擾 系 數 矩 陣[7]；A為3×3 地理坐標系與船舶坐標系之間的轉換矩陣[8-9]，該矩陣元素為測量船姿態角的函數；Be為3×1 三分量地磁場向量，nT；Bp為3×1 固定磁化干擾系數向量[10]，nT。其中，

式中，φ，θ，ψ 分別為船舶的航向角、橫搖角和縱傾角，可以由船舶的導航系統進行實時測量。

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

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

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

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

2 Kalman 濾波算法

2.1 算法原理

Kalman 濾波算法需以系統的時間離散化狀態空間為基礎[14]，在船舶結構及磁傳感器安裝位置不變的條件下，船舶磁化系數K和Bp均為時域不變量，當以4×3 的擴展矩陣x=[KBp]T（T 表示矩陣的轉置）為狀態量時，則測量過程的計算方程為

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

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

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

由式(5)和式(6)構成狀態空間，則自適應Kalman 濾波迭代算法的步驟如下：

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

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

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

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

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

2.2 算法仿真

為了驗證2.1 節算法的正確性及其對隨機誤差的抑制能力，本節將開展仿真試算。設定三分量地磁場向量為Be=[34 000, 2 000, 35 000] nT，感應磁化干擾系數矩陣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]，固定磁化干擾系數向量Bp=[1 500, -800, 1 800] nT。船舶姿態角φ，θ，ψ 分別在區間[0°，360°]，[-20°，20°]，[-15°，15°]取值，并生成轉換矩陣A。按照式(1)生成12 組傳感器輸出向量B，并施加最大值為100 nT 的隨機誤差。記錄每一組向量B所對應的φ，θ，ψ，并施加最大值為0.05°的隨機誤差。將數據進行延長處理，即可重復使用多次，從而保證算法具有足夠數量的迭代步數。此時，以Be，A，B為已知量，即可按照2.1 節的Kalman 濾波算法解算K和Bp，算法收斂情況如圖1 所示。

由圖1 可知，該Kalman 濾波算法經過約300 次迭代即可實現收斂，且具有良好的穩定性。由于本文在仿真過程施加了隨機誤差，其解算結果為：

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

與設定值相比，K的誤差出現在萬分位，Bp的誤差出現在個位。鑒于地磁場為104nT 量級，按式（3）的傳遞關系估算誤差，可知因船舶磁化干擾系數測量誤差所造成的地磁場解算誤差約為10 nT 量級，該誤差水平與文獻[3]中僅考慮隨機測量誤差所得的地磁場解算誤差水平相當。

3 船模驗證實驗

3.1 實驗模型及測量過程

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

圖2 實驗裝置示意圖Fig. 2 Experiment device

為了簡化測量過程，船模分別在磁東、磁北、磁西、磁南4 個主航向進行了3 種水平姿態下的磁干擾測量，共計獲得12 組測量數據，并同步記錄了三分量磁傳感器和水平角度傳感器的測量數據。

3.2 測量結果及計算分析

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

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

固定干擾系數向量為

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

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

圖3 干擾系數的實驗收斂曲線Fig. 3 Convergence curves of magnetization factors in experiment

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

圖4 三分量地磁場計算收斂曲線Fig. 4 The convergence curves in computing three-component geomagnetic field

表1 船模磁場測量數據表Table 1 Measured magnetic field data of the mockup

4 結 語

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

在實驗室模擬測量過程中，影響測量精度的主要原因為：三分量磁傳感器沒有經過固有誤差校正；測量干擾系數時，環境地磁場存在測量誤差；傳感器存在安裝誤差，尤其是傳感器x軸向與模型艏線、艉線方向的一致性誤差；角度傳感器僅能測量橫搖角和縱傾角，航向角僅采用了實驗室軌道方向作為近似值。因此，在實際工程應用中，還需進一步處理這些誤差，以提高船舶感應干擾系數的測量精度。