0 引言

近年来，随着无人机的应用发展，无人机航空磁力探测在航空物探领域发挥着重要的作用^[1-2]。在航磁测量中，飞机上磁性元件产生的磁场叠加在地磁场之上会对磁力仪产生很大干扰，此叠加干扰磁场随着飞机飞行方向和姿态变化。干扰补偿模型是1950年Tolles首先根据磁干扰产生的原因及性质将磁干扰分解为恒定磁场、感应磁场和涡流磁场，并建立了磁干扰场的数学模型Tolles-Lawson方程(T-L方程)。1961年，Leliak^[3]基于T-L方程设计了一套北—东—南—西航向的参考机动，干扰场是关于飞机机动的函数，从而建立了16项系数方程，根据系数求出干扰，最后将干扰从航磁数据中剔除；但该方法通过姿态角解算方向余弦的过程复杂。Leach^[4]给出了利用三轴磁通门数据获取与飞机机动有关的方向余弦。Zhang等^[5]认为不能直接获得目标磁场的方向余弦值，只能间接使用总场的强度来表示它。Zhou等^[6]将磁通门传感器的进行数据滤波后用于表示方向余弦。王婕等^[7]在固定翼无人机上安装磁通门开展补偿试验，但模型中忽略了涡流场对应的干扰项，模型中方向余弦的准确性对系数的稳定性有很大影响。吴佩林等^[8]在无人直升机上以三轴磁通门求解方向余弦，并通过前馈网络实现了数据补偿。

本文以Leach磁通门数据求解方向余弦的方法进行航磁补偿，建立三轴磁通门数据理论模型，并以此验证补偿方法效果。最后，将此方法应用云南某无磁异常地区试飞磁通门数据进行干扰补偿。

1 飞机磁场模型建立及求解 1.1 坐标系与飞机姿态

飞机坐标系^[9-10](图 1)坐标原点选在磁力仪探头的中心O；x_b轴与飞机横轴平行，左机翼方向为正；y_b轴与飞机纵轴平行，机首方向为正；z_b轴与x_b、y_b轴所在平面正交，垂直机身，向下为正。在飞行试验平台上，磁力仪的三轴分别平行于飞机坐标系的x_b、y_b、z_b轴。当地水平坐标系^[9-10]原点与飞机坐标系原点重合，x_g轴指向当地北子午线，y_g轴垂直于x_g轴指向东，z_g轴与x_g、y_g所在平面正交垂直向下。

飞机飞行姿态^[3]可以分解为侧滚角(ψ)、俯仰角(λ)和侧滑角(θ)。其中：ψ为x_b轴与x_gOy_g所在平面之间的夹角(x_h为x_b轴在x_gOy_g平面上的投影)，自左翼向下为正；λ为y_b轴与x_gOy_g所在平面之间的夹角，自机头向下为正；θ为y_b轴在x_gOy_g平面上投影y_h与x_g轴的夹角，由北逆时针为正。

1.2 磁补偿模型建立

飞机磁干扰场主要由恒定场、感应场和涡流场组成^[11]。恒定场(H_p)为由飞机内铁磁性物体的剩余磁化产生的磁场，感应场(H_i)为由飞机内铁磁性物体受地磁场磁化产生的磁场，涡流场(H_ec)为由飞机机身、机翼等软磁性金属材料切割地磁场产生的磁场，分别表示为：

式中：T表示地磁场；c₁，c₂，…，c₁₈为变量系数；α、β、γ分别表示飞机x_b、y_b、z_b轴与T之间的夹角；cos′α、cos′β、cos′γ表示方向余弦cosα、cosβ、cosγ关于时间t的微分。则飞机干扰总磁场H_t表示为

根据余弦函数关系式

并结合式(1)—(4)整理得到

式中, d₁，d₂，…，d₁₆为16项补偿系数。

1.3 方向余弦获取

三轴磁力仪测地磁场方向余弦存在测量误差^[5](图 1)，可能会对方向余弦产生一定影响。图 1中H_g表示H_t与T的矢量和。在磁通门求解方向余弦模型中，光泵磁力仪为原点O，飞机产生的干扰磁场强度大多在100~1 000 nT之间，并且传感器与飞机之间的距离越远，检测到的干扰越小。在本次试验飞行平台上，可将正北方向平飞时的三轴磁通门数据与地磁场在地面坐标系下的三分量值进行对比。经分析，干扰对地磁场矢量和磁通门矢量夹角余弦值产生的误差最大为1%；因光泵磁力仪固有特点，数据可用于16项补偿系数方程中干扰值H_t的求解。

假设T_bx，T_by，T_bz分别为地磁总场T_t在飞机坐标系下x_b、y_b、z_b轴方向的投影大小，在忽略飞机干扰产生误差的条件下，可认为

则关于时间连续的方向余弦可表示为

1.4 模型参数求解

最小二乘法是求解多元线性回归的常用方法。此外，还有Praga-Alejo等^[12]提出的岭回归估计改进算法等。假设光泵磁力仪和三轴磁通门有效采样点数为N，则根据式(6)可以列出含有n个方程的超定方程^[13]，线性方程可表示为

其中：

式中:

结合Leliak给出的参考机动，采用加拿大国立研究院航空处(National Aeronautical Establish-ment, NAE)提出的“1/2 FOM机动飞行”^[14]方法，即在北—东—南—西上做±10°侧滚、±5°俯仰、±5°侧滑的机动(图 2)，每种姿态时长5s，并间隔设置5s的水平飞行。

2 理论模型验证

T-L方程中假设补偿飞行的场地地磁场强度恒定，因此当地水平坐标系下的地磁场三分量与飞机坐标系下的磁场三分量之间存在着坐标系旋转的关系^[9-10]，即地磁场三分量T_gx、T_gy、T_gz依次转动θ、ψ、λ，可获得飞机坐标系下的磁场三分量T_bx、T_by、T_bz：

式中：

b_{1, 1}=cosθcosψ；b_{1, 2}=sinθcosψ；b_{1, 3}=-sinψ；b_{2, 1}=cosθsinλsinψ-sinθcosλ；b_{2, 2}=sinθsinλ· sinψ+cosθcosλ；b_{2, 3}=sinλcosψ；b_{3, 1}=cosθ· cosλsinψ+sinθsinλ；b_{3, 2}=sinθcosλsinψ- cosθ·sinλ；b_{3, 3}=cosλcosψ; Q取值在北航向为±5°, 在东航向为90°±5°，在南航向为180°±5°，在西航向为270°±5°。

以云南某地试验场磁场为例，由WMM2015^[15]模型计算得到T=46 861.6 nT, μ=-1.297 6°, φ=36.662 7°, (μ为磁偏角，φ为磁倾角)，则可得到T_gx、T_gy、T_gz，进而得到三轴与飞机坐标系平行的磁通门T_bx、T_by、T_bz数据模型。此模型要求飞机严格按照参考机动飞行(北—东—南—西依次进行)，否则难以验证模型误差。同时在光泵数据中设置20 nT的异常信号，用于验证补偿效果。

依照图 2所示参考机动设计姿态角理论模型，通过式(9)计算得到磁通门、光泵磁力仪理论模型分别如图 3、图 4所示。

本次研究基于Matlab软件编写了补偿程序。只需输入三轴磁通门磁力仪三分量数据、光泵磁力仪总场数据以及飞行姿态所在时间段数据，即可完成补偿。

通过理论模型验证补偿程序，图 5表示自行设计的补偿程序能检测出磁异常信号。

3 应用实例 3.1 飞行平台整合与补偿场地

试验飞行平台采用吉林大学自研重载荷智能化物探专用无人直升机^[16](图 6)，采用自研发靶向溅射技术喷射低磁防静电材料作为机身涂层，导电性能大幅提高，易于放电，由机动产生的干扰较小。平台搭载TFM100-G2型三轴磁通门磁力仪和加拿大GEM公司的GSMP-35A型高精度钾光泵磁力仪，配套航磁收录系统，该系统磁场范围为20 000~ 100 000 nT，采样率10 Hz。

依托前期在云南某试验场地开展的航磁系统性能测试，该试验场地地磁场变化小，地形较为平坦。并在此进行了补偿飞行试验，飞行路线如图 7所示，飞行高度平均155 m(基于大地水准面)。

3.2 数据滤波

测量数据包含许多不同类型的噪音，为得到只受飞机磁干扰的总场数据，在补偿过程中须对数据进行滤波^[17-18]。为了叙述方便，本文采用巴特沃斯滤波器以三轴磁通门数据为列进行滤波处理(图 8)，滤波后可更精确地估算16项补偿系数。

3.3 补偿计算

将滤波后的数据以最小二乘法解出16项补偿系数，再从总场值中剔除干扰场，从而得到补偿后的地磁场(图 9)。从图 9中补偿曲线可以看出补偿效果明显，自行设计的程序可以完成航磁补偿任务。

4 补偿质量评价

目前实际探测工作中主要使用的量化指标为NAE提出的补偿后残留磁干扰标准差σ^[19]，并以“改善比”(R_i)辅助评估补偿效果。σ是根据参考机动补偿后剩余磁干扰的标准差，小于标准值0.08视为合格；R_i为补偿前后标准差的比值，其值越高，说明补偿效果越好。经计算，σ=0.144，R_i=10.611；在高精度航磁测量中，磁日变将在航磁数据中引入较大的误差，因此，进行磁日变矫正是十分必要的。去除日变后，σ=0.069，R_i=19.364，对比去除日变前效果有明显提升，达到补偿标准值的要求。

5 结论

1) 机载磁测系统需要高精度磁干扰补偿算法，基于T-L方程建立了飞机磁干扰补偿模型，并计算了16项补偿系数。以坐标旋转关系构建的磁通门理论模型，可在飞行试验前用于验证补偿模型的可信度。三轴磁通门理论模型数据经补偿后，补偿系统能消除干扰场，并保留异常信号。

2) 利用磁通门三轴数据解算方向余弦可免去方向余弦繁琐的计算过程。将本文设计的补偿算法用于云南某试验区的实测航磁数据，采用标准差能准确地评价磁补偿效果，实际补偿效果较好；基于磁通门估计的航磁补偿方法是有效的。