闫辉，周国华
1 信阳学院 大数据与人工智能学院，河南 信阳 464000
2 海军工程大学 电气工程学院，湖北 武汉 430033
置于地球磁场中的铁磁物体，例如舰艇、坦克、车辆、飞机等，都会因地磁场作用而产生磁化磁场[1]，所以安装在船舶上的磁传感器会受到船舶磁化磁场的干扰。如何从磁传感器的测量数据中分离出三分量地磁场，是以船舶为载体的地磁场测量系统所面临的核心问题。舰载消磁系统是一种以三分量地磁场作为系统控制信号[2]来补偿船舶磁化磁场的船用设备，该系统采用安装于舰艇桅杆上经过抗干扰调整的三分量磁传感器来实时监测地磁场，进而控制消磁绕组电流以补偿舰船磁场，即“测地磁消舰磁”，其中抗干扰调整的目的是消除船舶磁化磁场对地磁测量传感器的干扰。
为了消除船舶磁化磁场对地磁测量传感器的干扰，需预先测量船舶的感应干扰系数矩阵和固定干扰系数向量[3]。根据以船舶为载体的三分量地磁场测量系统数学模型，肖昌汉[4]提出了“四航向法”，但该方法只能测量感应干扰系数矩阵的前2 列元素，且需在不同的两地进行补充测量。肖昌汉等[5]提出了“8”字航向法，即要求船舶在已知地磁场的测量场地进行高速、大舵角的变航向机动，进而造成较大幅度的船舶姿态变化，并需要大量的有效磁场测量数据作为支撑。闫辉等[6]对船载三分量地磁场的测量过程进行了实验模拟，发现即使在实验室条件下，船舶磁化干扰系数的测量也存在一定困难。因此，船舶磁化磁场干扰系数的测量问题是限制船载三分量地磁场测量系统大规模工程应用的难点。
为了降低测量船舶磁化磁场干扰系数的工程难度和成本，本文拟基于船舶为载体的三分量地磁场测量系统数学模型，分析解决船舶磁化磁场干扰的难点问题。首先，将根据三分量地磁场解算的Kalman 滤波方法[3]，重新建立适用于船舶磁化干扰系数测量的状态空间方程，提出采用Kalman滤波方法求解干扰系数的自适应算法，并通过计算机仿真以验证该算法的理论正确性；
然后，将开展实验室模拟测量，在有限的测量数据样本条件下，应用解算所得的干扰系数来有效排除船舶磁化磁场的干扰，从而分离出较高精度的三分量地磁场，用以为船舶磁化干扰系数测量提供一条高效率、低成本的可行路径。
三分量磁传感器作为磁场的敏感元件，可以采集其所处位置的各种矢量合成磁场，其中既包括地磁场，也包括来自测量船自身的磁化磁场干扰。安装于船舶桅杆的三分量磁传感器的动态输出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.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.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 滤波算法，即可得到：
固定干扰系数向量为
由于测量值的数据量过小（仅12 组），导致滤波器的迭代步数过少，因此无法收敛。为了解决这一问题，本文将这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
对于以船舶为载体的三分量地磁场测量系统，其核心问题是船舶磁化干扰系数的测量。本文提出了基于Kalman 滤波的船舶磁化干扰系测量算法，根据安装于船舶的三分量磁场测量系统的传感器测量输出数学模型，建立了观测方程；
根据船舶干扰系数的相对不变性，建立了过程方程。由观测方程和过程方程所构成的状态空间表达式较好地契合了Kalman 滤波对于模型准确性的严苛要求，并通过计算机仿真验证了该算法理论的正确性。在较少实验数据样本的情况下，仅通过数据的简单重复利用，即可保证Kalman 滤波算法的良好收敛性。相较于“四航向法”，Kalman滤波方法可以测量所有磁化干扰系数，且无需严格按照东、西、南、北4 个主航向实施测量；
相较于“8”字航向法，该方法的数据量小，且无需按某种特定航迹实施测量。因此，从干扰系数测量的完备性，所需数据量的大小，以及测量实施条件的难易程度而言，Kalman 滤波方法可以显著提高船舶磁化干扰系数的测量效率，并降低工程实施的难度。
在实验室模拟测量过程中，影响测量精度的主要原因为：三分量磁传感器没有经过固有误差校正；
测量干扰系数时，环境地磁场存在测量误差；
传感器存在安装误差，尤其是传感器x轴向与模型艏线、艉线方向的一致性误差；
角度传感器仅能测量横摇角和纵倾角，航向角仅采用了实验室轨道方向作为近似值。因此，在实际工程应用中，还需进一步处理这些误差，以提高船舶感应干扰系数的测量精度。