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]，该矩阵元素为测量船姿态角的函数；B_e为3×1三分量地磁场向量，nT；B_p为3×1固定磁化干扰系数向量^[10]，nT。其中，

式中，φ，θ，ψ分别为船舶的航向角、横摇角和纵倾角，可以由船舶的导航系统进行实时测量。

根据式（1），如果需从传感器输出B的测量值中分离出地磁场B_e，则必须先测量船舶磁化干扰系数K和B_p。将式（1）展开并整理，可知在不同的舰船姿态下，传感器的任意输出分量为

式中：B_i为磁传感器第i轴向的输出，其中i=1，2，3，分别对应x轴、y轴、z轴；k_i1，k_i2，k_i3分别为矩阵K中第i行的第1，2，3列元素；f₁，f₂，f₃分别为式（1）中的函数映射关系；B_pi为3×1固定磁化干扰系数向量B_p在第i轴向的分量。

在已知地磁场和船舶姿态的条件下，需至少记录4组不同姿态下的传感器输出，才能通过求解一个线性方程组来得到船舶的磁化干扰系数。将该系数代入式(1)，则未知量地磁场B_e的解为

上述地磁场的解算过程虽然在理论上成立，但在实际工程运用中却难以保证解算精度，主要原因如下：1）三分量磁通门传感器普遍存在固有的零位误差、标度误差和正交误差^[11]；2）船舶姿态角φ，θ，ψ数据主要来源于导航系统，存在一定的漂移误差，以及信号传输过程中的延迟误差；3）由于背景环境的磁干扰及地磁场静日变化，将导致一定的磁传感器随机误差^[12]。由于船舶磁化干扰系数对以上各种误差源较敏感，所以这些误差将在测量过程中累积，进而通过式（4）向地磁场解算值传递。目前，除了第1种误差之外，其他2类误差暂无较好的处理方法。因此，为了提高三分量地磁场的解算精度，需采用相应措施来消除船舶磁化干扰系数的测量误差，从而提高式（1）数学模型的准确性。广泛应用于动态测量系统的Kalman滤波方法在原理上可以较好地抑制系统测量误差^[3,13]，且可以通过信息的引入对变化状态进行及时更新。将Kalman滤波应用于磁化干扰系数测量时，必须以磁化干扰系数作为状态量，但该状态具有相对不变性^[4]，因此磁化干扰系数测量具有静态测量系统的特点。如果将静态系统视为动态系统的一种特例，则Kalman滤波可为解决船舶磁化干扰系数测量问题提供一条新途径。

2.   Kalman滤波算法

2.1   算法原理

Kalman滤波算法需以系统的时间离散化状态空间为基础^[14]，在船舶结构及磁传感器安装位置不变的条件下，船舶磁化系数K和B_p均为时域不变量，当以4×3的扩展矩阵x=[K B_p]^T（T表示矩阵的转置）为状态量时，则测量过程的计算方程为

式中：x(n)为状态量，测量点数n=0，1，2，…，N，其中N为最大测量点数；4×3向量v₁(n)为过程噪声向量，用于描述状态转移误差。

根据式(1)，即可建立观测方程：

式中：3×1向量y(n)= B^T ，1×4向量$\boldsymbol{C}(n) = [ {{{\boldsymbol{B}}_{\text{e}}}^{\text{T}}{{\boldsymbol{A}}^{\text{T}}}}\;\;1 ]$，3×1向量v₂(n)为观测噪声向量，用于描述随机噪声。

由式(5)和式(6)构成状态空间，则自适应Kalman滤波迭代算法的步骤如下：

1）设置初始条件：状态x(0)，中间变量P(0)（4×4对角矩阵）。

2）已知参数：4×4过程噪声方差阵Q，一维观测噪声方差σ²。

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

4）计算中间量及状态量，其中：R(n)，P(n)， m(n)均为中间变量；E为单位矩阵。

（1）${\boldsymbol{R}}(n) = {\boldsymbol{P}}(n - 1) + {\boldsymbol{Q}}$

（2）${\boldsymbol{m}}(n) = {\boldsymbol{R}}(n){\boldsymbol{C}}{(n)^{\text{T}}}[{\boldsymbol{C}}(n){\boldsymbol{R}}(n){\boldsymbol{C}}{(n)^{\text{T}}} + {\boldsymbol{\sigma }^2}]$

（3）${\boldsymbol{P}}(n) = [{\boldsymbol{E}} - {\boldsymbol{m}}(n){\boldsymbol{C}}(n)]{\boldsymbol{R}}(n)$

（4）${\boldsymbol{x}}(n) = [{\boldsymbol{E}} - {\boldsymbol{m}}(n){\boldsymbol{C}}(n)]{\boldsymbol{x}}(n - 1) + {\boldsymbol{m}}(n){\boldsymbol{y}}(n)$

通过以上算法，即可得到状态量x(n)^T，前3列可以构成感应干扰系数矩阵K，第4列则可以构成固定干扰系数向量B_p。在上述迭代过程中，初始条件的设置具有一定的任意性，对于状态先验信息不足的条件下，初始状态x(0)取零值，中间变量初始值P(0)取较小值将有利于迭代过程的收敛。作为已知参数，过程噪声方差阵Q和观测噪声σ²的取值对迭代过程的收敛速度和收敛稳定性具有决定性作用：当这2个参数取较大值时，收敛速度快，但稳定性较差；反之，则收敛速度慢，而稳定性较好，因此，本文将通过模拟仿真试算的方法来平衡确定参数值。

2.2   算法仿真

为了验证2.1节算法的正确性及其对随机误差的抑制能力，本节将开展仿真试算。设定三分量地磁场向量为B_e=[34 000, 2 000, 35 000] nT，感应磁化干扰系数矩阵K=[ k₁₁，k₁₂，k₁₃; k₂₁，k₂₂，k₂₃; k₃₁，k₃₂，k₃₃]=[0.95, −0.003, 0.003; 0.003, 0.95, −0.03; 0.03, −0.03, 1.05]，固定磁化干扰系数向量B_p=[1 500, −800, 1 800] nT。船舶姿态角φ，θ，ψ分别在区间[0°，360°]，[−20°，20°]，[−15°，15°]取值，并生成转换矩阵A。按照式(1)生成12组传感器输出向量B，并施加最大值为100 nT的随机误差。记录每一组向量B所对应的φ，θ，ψ，并施加最大值为0.05°的随机误差。将数据进行延长处理，即可重复使用多次，从而保证算法具有足够数量的迭代步数。此时，以B_e，A，B为已知量，即可按照2.1节的Kalman滤波算法解算K和B_p，算法收敛情况如图1所示。

图  1  磁化干扰系数的仿真收敛曲线

Figure  1.  Convergence curves of magnetization factors in simulation

下载: 全尺寸图片 幻灯片

由图1可知，该Kalman滤波算法经过约300次迭代即可实现收敛，且具有良好的稳定性。由于本文在仿真过程施加了随机误差，其解算结果为：

与设定值相比，K的误差出现在万分位，B_p 的误差出现在个位。鉴于地磁场为10⁴ nT量级，按式（3）的传递关系估算误差，可知因船舶磁化干扰系数测量误差所造成的地磁场解算误差约为10 nT量级，该误差水平与文献[3]中仅考虑随机测量误差所得的地磁场解算误差水平相当。

3.   船模验证实验

3.1   实验模型及测量过程

为了检验船舶磁化干扰系数测量时Kalman滤波方法的有效性，本节将在无磁实验室开展某型船模的验证实验。测量装置由三分量磁通门传感器、水平角度传感器及硬质无磁不锈钢支架组成，其中磁传感器与角度传感器通过支架进行刚性连接，需注意调整相应测量轴。将该测量装置安装于某型磁性船模上，测量轴与舰船坐标系保持一致，如图2所示。将磁性船模放置于无磁平台，该平台可以满足船模4个主航向旋转、横摇和纵倾要求，据此即可构建船载地磁三分量测量平台。

图  2  实验装置示意图

Figure  2.  Experiment device

下载: 全尺寸图片 幻灯片

为了简化测量过程，船模分别在磁东、磁北、磁西、磁南4个主航向进行了3种水平姿态下的磁干扰测量，共计获得12组测量数据，并同步记录了三分量磁传感器和水平角度传感器的测量数据。

3.2   测量结果及计算分析

测量数据如表1所示，其中磁传感器三分量、横摇角、纵倾角分别为传感器实测值，磁东、磁北、磁西、磁南4个航向角则直接取值90°，0°，270°，180°，并转换为弧度。

表  1  船模磁场测量数据表

Table  1.  Measured magnetic field data of the mockup

参数序号	磁传感器三分量			航向角ψ/rad	横摇角θ/rad	纵倾角φ/rad
	B1 /nT	B2 /nT	B3 /nT
1	4 206	25 290	36 913	4.712 39	0.025 22	−0.004 26
2	−30 595	−2 785	40 620	3.141 59	0.024 73	−0.003 02
3	2 478	−31 825	43 890	1.570 80	0.023 21	−0.002 62
4	37 118	−3 690	40 110	0.000 00	0.023 09	−0.003 91
5	37 485	−9 719	39 480	0.000 00	−0.174 01	−0.014 15
6	2 608	−36 985	36 580	1.570 80	−0.173 47	−0.012 90
7	−30 197	−8 205	41 110	3.14159	−0.171 83	−0.013 28
8	4 761	19 020	43 985	4.712 39	−0.171 27	−0.014 54
9	37 811	2 408	37 845	0.000 00	0.222 70	−0.021 92
10	4 514	30 860	28 158	4.712 39	0.224 57	−0.022 46
11	−29 836	3 038	39 796	3.141 59	0.223 66	−0.022 48
12	3 615	−25 322	49 387	1.570 80	0.222 16	−0.022 11

下载: 导出CSV 

| 显示表格

已知环境地磁场为B_e=[34 425, 1961, 35 898] nT，采用上文Kalman滤波算法，即可得到：

固定干扰系数向量为

B_p=[2 011.1, −2 886.8, 2 054.6] nT。

由于测量值的数据量过小（仅12组），导致滤波器的迭代步数过少，因此无法收敛。为了解决这一问题，本文将这12组数据进行了延长处理，即进行12组数据的多次重复使用，并以此作为实验数据集。从实验数据集中抽取一定数量的数据，即可解算磁化干扰系数，图3所示为感应干扰系数和固定干扰系数的迭代过程曲线。

图  3  干扰系数的实验收敛曲线

Figure  3.  Convergence curves of magnetization factors in experiment

下载: 全尺寸图片 幻灯片

由图3可知，对于测量数据较少的情况，通过对数据进行重复性延长处理，Kalman滤波算法仍然具有较好的收敛性；但与仿真结果相比，其收敛之后的稳定性略有下降，存在微幅震荡。由于仿真所采用的数据均为无误差的理想值，而实验中传感器所测量的磁场、航向角、横摇角、纵倾角均存在误差，因此，造成稳定性下降的主要原因是实验数据的测量误差。

为了验证表1船舶磁化干扰系数的正确性，将其代入文献[3]的三分量地磁场解算算法，并从实验数据集中抽取一定样本的数据用于解算地磁场，其收敛过程如图4所示，其中B_ei为3×1三分量地磁场向量B_e在第i轴向的分量（i=1，2，3）。根据计算结果，收敛稳定后的地磁场解算数据的最大相对误差为1.9%（此处定义为：误差最大值与地磁场总量的比值），均方根误差为266 nT。由此可见，在仅有12组测量数据的条件下，Kalman滤波算法具有较高的计算精度，其对测量误差具有较好的抑制能力。

图  4  三分量地磁场计算收敛曲线

Figure  4.  The convergence curves in computing three-component geomagnetic field

下载: 全尺寸图片 幻灯片

4.   结　语

对于以船舶为载体的三分量地磁场测量系统，其核心问题是船舶磁化干扰系数的测量。本文提出了基于Kalman滤波的船舶磁化干扰系测量算法，根据安装于船舶的三分量磁场测量系统的传感器测量输出数学模型，建立了观测方程；根据船舶干扰系数的相对不变性，建立了过程方程。由观测方程和过程方程所构成的状态空间表达式较好地契合了Kalman滤波对于模型准确性的严苛要求，并通过计算机仿真验证了该算法理论的正确性。在较少实验数据样本的情况下，仅通过数据的简单重复利用，即可保证Kalman滤波算法的良好收敛性。相较于“四航向法”，Kalman滤波方法可以测量所有磁化干扰系数，且无需严格按照东、西、南、北4个主航向实施测量；相较于“8”字航向法，该方法的数据量小，且无需按某种特定航迹实施测量。因此，从干扰系数测量的完备性，所需数据量的大小，以及测量实施条件的难易程度而言，Kalman滤波方法可以显著提高船舶磁化干扰系数的测量效率，并降低工程实施的难度。

在实验室模拟测量过程中，影响测量精度的主要原因为：三分量磁传感器没有经过固有误差校正；测量干扰系数时，环境地磁场存在测量误差；传感器存在安装误差，尤其是传感器x轴向与模型艏线、艉线方向的一致性误差；角度传感器仅能测量横摇角和纵倾角，航向角仅采用了实验室轨道方向作为近似值。因此，在实际工程应用中，还需进一步处理这些误差，以提高船舶感应干扰系数的测量精度。