基于稳健估计的无人车磁干扰补偿方法

A Robust-Estimation-Based Magnetic Interference Compensation Method for Unmanned Ground Vehicles

摘要

地磁导航因其无源、 稳定和抗干扰等优点， 已广泛应用于无人载体自主定位与导航中。高质量地磁信息获取是实现精确导航的关键。然而， 无人车受机动特征限制和复杂磁环境影响， 测量数据中常常同时存在病态、 噪声及粗差等问题， 降低了磁补偿的精度与鲁棒性。现有载体磁干扰补偿模型中的噪声一般假设服从高斯白噪声或近似高斯白噪声， 这并不全面。为此， 提出融合t型抗差估计的截断总体最小二乘（Truncated Total Least Squares， TTLS）补偿方法（TTLS-t）， TTLS能够解决磁测模型两侧都具有误差的严重病态方程问题， 在此基础上， 引入重尾t分布可增强对多类型误差的鲁棒性， 通过t自由度参数自适应检测来灵活调节尾部厚度， 利用迭代重加权最小二乘（Iteratively Reweighted Least Squares， IRLS）并基于残差构建反映观测点粗差置信度的权重矩阵， 削弱误差影响， 实现抗差估计， 最后， 基于EM（Expectation Maximization）算法实现磁干扰的稳健补偿。实测结果表明， TTLS-t算法在数据质量不理想条件下， 相较传统方法， 补偿精度提升大于75%， 具有更高的精度和鲁棒性。

Abstract

Geomagnetic navigation， characterized by its passive nature， stability， and strong anti-interference capability， has been widely applied in autonomous positioning and navigation of unmanned vehicles. High-quality geomagnetic information acquisition is a prerequisite for achieving accurate navigation performance. However， constrained maneuverability and complex magnetic environments often lead to the coexistence of ill-conditioned structures， noise， and gross errors in magnetic measurements， which significantly degrade the accuracy and robustness of magnetic compensation. Most existing magnetic disturbance compensation models assume Gaussian or near-Gaussian noise， which is insufficient to characterize real measurement errors. To address these limitations， this paper proposes a truncated total least squares compensation method with t-type robust estimation （TTLS-t）. The TTLS framework effectively handles severely ill-conditioned magnetic models with errors on both sides， while the incorporation of a heavy-tailed t-distribution enhances robustness against multiple error sources. The tail thickness is adaptively adjusted through the degrees-of-freedom parameter， and an iteratively reweighted least squares （IRLS） scheme is employed to construct a residual-based weighting matrix that reflects the confidence in individual observations， thereby suppressing the influence of gross errors. The overall magnetic interference compensation is realized within the expectation-maximization （EM） framework. Experimental results based on measured data demonstrate that， under degraded data quality conditions， the proposed TTLS-t method achieves more than 75% improvement in compensation accuracy compared with conventional approaches， exhibiting superior accuracy and robustness.

Graphical abstract

关键词

地磁导航 / 磁干扰补偿 / 截断总体最小二乘 / t型估计 / 迭代重加权最小二乘 / EM算法

Key words

geomagnetic navigation / magnetic interference compensation / truncated total least squares / t-type estimation / iteratively reweighted least squares / expectation maximization algorithm

0 引　言

地磁导航是一种利用地磁场空间分布， 根据指定区域内地磁基准图与载体运动中获取的区域磁场测量序列来确定载体位置信息的导航方案。随着全球导航卫星系统（Global Navigation Satellite System， GNSS）受限与复杂环境下自主导航需求的提升， 地磁导航技术以其无源、 隐蔽性强与误差不随时间累积等优势， 已在多类载体与应用场景中得到广泛研究与验证^［1‐2］。例如， 在航空与无人机平台， 地磁异常匹配被用于GNSS拒止条件下的辅助或备份导航^［3‐4］； 在地面载体与机器人领域， 地磁特征被用于复杂环境下的长期稳定定位与多传感器融合导航^［5］； 水下环境中， 地磁导航作为水下航行器的重要自主定位手段， 常与惯性与声学导航相结合^［6］。相关研究表明， 地磁导航在空中、 地面和水下等多种军事与民用场景中均具备明确的应用背景与工程潜力。

高质量、 高精度地磁场信息的获取是实现精确地磁导航的基础。无人车以其测量地磁异常信息丰富、 成本低、 机动强等特点在磁测领域具备独特优势^［7‐8］。

目前， 无人车及其他小型平台的地磁导航/磁补偿常采用椭球拟合（Ellipsoid Fitting， EF）^［9］、 最小二乘（Least Square， LS）^［10］、 总体最小二乘（Total Least Square， TLS）及岭回归^{［11‐12］}等模型驱动校正方法。这类方法基于最小残差准则， 在噪声较小、 数据稳定时能取得良好效果， 但其高斯噪声假设难以应对多源干扰与病态问题^［13］。为此， 研究者提出截断总体最小二乘法（Truncated Total Least Squares， TTLS）， 通过截断奇异值缓解病态矩阵影响， 在噪声与病态共存条件下提升了补偿精度^［14］。然而， TTLS仍对异常观测高度敏感， 难以适应地面复杂磁环境中噪声与粗差并存的情况， 参数估计稳定性欠佳。

针对最小二乘估计对异常观测点敏感、 鲁棒性不足的问题， 统计学领域研究者引入重尾分布模型与加权估计框架以提高估计对粗差与非高斯噪声的适应能力^{［15‐16］}。t分布作为高斯分布的重尾扩展形式， 可通过自由度参数灵活调节尾部厚度， 以在不同噪声条件下兼顾估计精度与抗干扰能力^［17］。

在数据驱动与智能算法方面， 神经网络和深度学习已广泛应用于无人车、 无人机的地磁导航与磁干扰补偿， 此类方法利用非线性逼近与特征自学习机制， 克服了传统模型在复杂电磁环境下的建模不足。典型研究包括基于径向基函数（Radial Basis Function， RBF）或反向传播（Back Propagation， BP）网络的非线性补偿^{［18‐19］}、 基于卷积神经网络（Convolutional Neural Network， CNN）和长短期记忆网络（Long Short-Term Memory， LSTM）的时空特征提取^［20］以及结合物理模型与神经网络的混合建模框架^{［21‐22］}。这些方法具备强大的拟合与自适应能力， 可实现高精度补偿， 但对训练数据分布和质量依赖较高， 在存在粗差、 重尾噪声或病态问题时， 稳定性与鲁棒性无法保证。

综上所述， 模型驱动类方法虽具备良好的可解释性与数学严谨性， 但现有方法在多源干扰与异常观测条件下稳定性不足； 而数据驱动类方法虽能刻画复杂非线性关系， 却过度依赖训练样本分布， 对噪声与粗差敏感且缺乏数值约束。并且在实际工程应用中， 无人车磁补偿算法通常需部署于资源受限的嵌入式平台， 对计算复杂度、 存储开销及在线实时性具有严格要求， 而数据驱动类方法往往依赖大规模离线训练与较高的在线推理成本， 难以在保证实时性与稳定性的前提下实现可靠部署。

针对无人车实测环境中噪声、 粗差与观测病态并存的问题， 亟需一种无需训练数据且兼具稳健性与数值稳定性的通用方案。本文在模型驱动框架下， 融合TTLS的病态抑制能力与t型估计的抗差特性， 提出TTLS-t抗差估计方法。该方法基于TTLS， 在EM算法框架下以IRLS形式实现迭代加权， 通过E步计算观测权重、 M步执行加权TTLS拟合， 有效抑制异常观测影响。TTLS-t在保持TTLS数值稳定性的同时， 具备对重尾噪声与多类型误差的鲁棒性， 为解决复杂磁环境下无人车磁干扰的精确建模难题， 提出了一种高精度、 强稳健的补偿方案。本文算法整体实现框架如图 1 所示。

1 磁干扰分析与模型构建

无人车搭载的磁通门传感器在地磁导航与磁测任务中， 受到车体结构、 电子设备及环境磁场等多源干扰影响， 导致三轴磁场数据出现偏移、 拉伸和畸变， 影响磁干扰补偿精度^{［23‐24］}。本节围绕磁干扰机理、 建模方法及抗差补偿算法展开分析。

1.1 磁干扰源分析

无人车磁测数据会受车体软硬磁及外部环境干扰影响： 软磁效应源于车体金属结构和电子元件， 使磁力计响应产生非正交和尺度不均的偏移； 硬磁效应源于永久磁体、 通电线缆和电机部件， 引入方向无关的恒定偏置^{［25‐26］}。此外， 随机噪声和粗差常叠加于观测数据中， 噪声近似高斯分布且方差随环境变化， 粗差则可能因突发电磁干扰或传感器瞬态饱和而呈重尾分布。

在地磁矢量恒定的背景磁场下测量时， 传感器的观测矢量必然是其经过软硬磁效应变换， 并叠加观测误差后的结果。设真实地磁矢量为

{\mathit{b}}_{\mathrm{e}}\in {\mathbb{R}}^{\mathrm{3}},‖{\mathit{b}}_{\mathrm{e}}‖=F,（1）

式中： F为当地地磁场强度常数。传感器测量值{\mathit{b}}_{\mathrm{T}}可建模为

{\mathit{b}}_{\mathrm{T}}=\mathit{S}{\mathit{b}}_{\mathrm{e}}+\mathit{H}+\mathit{\epsilon },（2）

式中： \mathit{S}\in {\mathbb{R}}^{\mathrm{3}\times \mathrm{3}}为软磁效应矩阵， 表征载体结构材料对外部磁场的线性畸变效应， 反映传感器三轴之间的比例因子误差、 轴间非正交性以及坐标轴方向耦合等特性； \mathit{H}\in {\mathbb{R}}^{\mathrm{3}}为硬磁效应向量， 反映来源于载体自身的恒定磁场； \epsilon为观测噪声以及粗差。

由式（2）可知， 磁通门在存在软、 硬磁干扰时， 其测量输出分布将符合一个椭球二次曲面方程。为恢复传感器输出与真实地磁场之间的线性关系， 构造补偿矩阵\mathit{M}\in {\mathbb{R}}^{\mathrm{3}\times \mathrm{3}}与补偿向量\mathit{H}\in {\mathbb{R}}^{\mathrm{3}}， 使

{\mathit{b}}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}=\mathit{M}\left(\begin{array}{c}{\mathit{b}}_{\mathrm{T}}-\mathit{H}\end{array}\right),（3）

式中： {\mathit{b}}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}为补偿后的磁场三分量， 理想状况下， 其模长满足‖{\mathit{b}}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}‖\approx F， 即所有补偿后的样本点应尽可能分布于以原点为球心， 半径为F的圆球面上。

软磁效应矩阵与硬磁效应向量的准确求解是精确补偿的关键， 故针对原始磁通门观测三分量进行误差建模。

设{\mathit{b}}_i=(m_{i,x},m_{i,y},m_{i,z}{)}^{\mathrm{T}}为第i次测量值， 其满足如下二次曲面关系，

\begin{array}{c}{q}_{\mathrm{1}}{m}_{i,x}^{\mathrm{2}}+q_{\mathrm{2}}{m}_{i,y}^{\mathrm{2}}+q_{\mathrm{3}}{m}_{i,z}^{\mathrm{2}}+q_{\mathrm{4}}{m}_{i,x}{m}_{i,y}+q_{\mathrm{5}}{m}_{i,x}{m}_{i,z}+\\ q_{\mathrm{6}}{m}_{i,y}{m}_{i,z}+q_{\mathrm{7}}{m}_{i,x}+q_{\mathrm{8}}{m}_{i,y}+q_{\mathrm{9}}{m}_{i,z}+q_{\mathrm{10}}=\mathrm{0},\end{array}（4）

式中： q_{\mathrm{1}}~q_{\mathrm{6}}对应二次项系数， 描述空间拉伸与耦合（软磁效应）； q_{\mathrm{7}}~q_{\mathrm{9}}对应一次项系数， 描述整体偏移（硬磁效应）； q_{\mathrm{10}}为常数项。

为便于矩阵化处理， 将式（4）重写为向量内积形式

{\mathit{d}}_i^{\mathrm{T}}\mathit{q}=\mathrm{0},（5）

其中，

\begin{array}{c}{\mathit{d}}_i=[m_{i,x}^{\mathrm{2}},m_{i,y}^{\mathrm{2}},m_{i,z}^{\mathrm{2}},m_{i,x}{m}_{i,y},m_{i,x}{m}_{i,z},\\ m_{i,y}{m}_{i,z},m_{i,x},m_{i,y},m_{i,z},\mathrm{1}],\hspace{1em}\\ \mathit{q}=[q_{\mathrm{1}},q_{\mathrm{2}},q_{\mathrm{3}},q_{\mathrm{4}},q_{\mathrm{5}},q_{\mathrm{6}},q_{\mathrm{7}},q_{\mathrm{8}},q_{\mathrm{9}},q_{\mathrm{10}}]\mathrm{。}\end{array}（6）

因此， 每一条观测满足约束{\mathit{d}}_i^{\mathrm{T}}\mathit{q}=\mathrm{0}， 假设有N条观测数据， 便可得到约束矩阵\mathit{D}\in {\mathbb{R}}^{N\times \mathrm{10}}为

\mathit{D}={\left[{\mathit{d}}_{\mathrm{1}}^{\mathrm{T}},{\mathit{d}}_{\mathrm{2}}^{\mathrm{T}},\cdots ,{\mathit{d}}_N^{\mathrm{T}}\right]}^{\mathrm{T}}\mathrm{。}（7）

从而， 总体约束为

\mathit{D}\mathit{q}\approx \mathrm{0}\mathrm{。}（8）

由于实测数据常含有各种形式的噪声与粗差， 式（8）在实际中仅能近似满足。若\mathit{D}\mathit{q}=\mathrm{0}成立， 则磁测三分量分布完全落在某一椭球面上； 若利用数值方法求得\mathit{q}， 则其参数将包含软磁效应矩阵、 硬磁效应向量等物理补偿量的全部信息。

1.2 截断总体最小二乘法

在实际应用中， 将式（8）等价为未知量与特定分量分离的形式， 令其中一类变量移至等号右边， 从而构造矩阵\mathit{A}\in {\mathbb{R}}^{N\times \mathrm{9}}与向量\mathit{b}\in {\mathbb{R}}^N。具体而言， 以m_z^{\mathrm{2}}作为等号右端， 便于在三轴磁测方程中形成标准的增广矩阵结构， 同时可在保持物理约束关系的前提下， 改善矩阵条件数， 提高TTLS求解的数值稳定性， 为后续的加权迭代与抗差估计提供可靠基础。

\mathit{A}\mathit{\varphi }=\mathit{b},（9）

式中： \mathit{\varphi }为9维待估参数向量； \mathit{b}=[m_{\mathrm{1},z}^{\mathrm{2}},\dots ,m_{N,z}^{\mathrm{2}}{]}^{\mathrm{T}}。

于是， 得到超定方程组

\mathit{A}\varphi \approx \mathit{b}\mathrm{。}（10）

为了避免普通最小二乘在高噪声和病态矩阵下的不稳定性， 首先采用TLS。其核心思想是： 允许\mathit{A}与\mathit{b}同时存在扰动， 寻找最小扰动使系统严格可解。构建增广矩阵为

\mathit{G}=[\mathit{A},\mathit{b}]\in {\mathbb{R}}^{N\times \mathrm{10}}\mathrm{。}（11）

对式中\mathit{G}作奇异值分解（Singular Value Decomposition， SVD）为

\mathit{G}=\mathit{U}\mathit{\Sigma }{\mathit{V}}^{\mathrm{T}},（12）

式中： \mathit{U}， \mathit{V}为正交矩阵； \mathit{\Sigma }为奇异值对角矩阵。

在理想TLS中， 解由\mathit{V}的最小奇异值对应列（最后一列）给出。然而， 在强噪声条件下， 最小奇异值可能放大误差。为此， 引入TTLS， 通过舍弃若干小奇异值（即\mathit{\Sigma }对角线的值）以平衡拟合精度与稳定性。

利用L曲线法来确定截断参数k。L曲线即以参数向量的范数‖\mathit{\varphi }‖为纵坐标， 残差的范数‖\mathit{A}\mathit{\varphi }-\mathit{b}‖为横坐标来绘制图线， 其结果会拟合为一条“L”形曲线。曲线的垂直与水平部分分别代表解与残差的范数随截断参数变化的趋势， 曲线的拐点处恰好对应残差与解范数的折中点， 此时既可保证较好的拟合精度， 又能保持解的稳定性， 因此该位置对应的k认为是最优截断参数。

将正交矩阵\mathit{V}分块为

\mathit{V}=\left[\begin{array}{cc}{\mathit{V}}_{\mathrm{11}}& {\mathit{V}}_{\mathrm{12}}\\ {\mathit{V}}_{\mathrm{21}}& {\mathit{V}}_{\mathrm{22}}\end{array}\right],（13）

式中： {\mathit{V}}_{\mathrm{12}}\in {\mathbb{R}}^{\mathrm{9}\times (\mathrm{10}-k)}， {\mathit{V}}_{\mathrm{22}}\in {\mathbb{R}}^{\mathrm{1}\times (\mathrm{10}-k)}。则有TTLS解为

\widehat{\varphi }=-{\mathit{V}}_{\mathrm{12}}{\mathit{V}}_{\mathrm{22}}^{†}\mathrm{。}（14）

式（14）给出了在截断参数k下的稳健解估计， 该结果即为后续椭球参数恢复和补偿矩阵构造的核心输入。根据截断后所得的特征子空间， 可将降维后的稳健解重新映射至原始参数空间， 从而恢复完整的原始二次曲面参数q。

基于TTLS解， 构建描述矢量磁传感器的零偏中心在三维坐标中位置的中心向量为

{\mathit{x}}_{\mathrm{0}}=-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{A}}_{\mathrm{e}}^{-\mathrm{1}}[q_{\mathrm{7}},\begin{array}{c}{q}_{\mathrm{8}},\end{array}{q}_{\mathrm{9}}{]}^{\mathrm{T}},（15）

{\mathit{A}}_{\mathrm{e}}=\left[\begin{array}{ccc}{q}_{\mathrm{1}}& \frac{q_{\mathrm{4}}}{\mathrm{2}}& \frac{q_{\mathrm{5}}}{\mathrm{2}}\\ \frac{q_{\mathrm{4}}}{\mathrm{2}}& q_{\mathrm{2}}& \frac{q_{\mathrm{6}}}{\mathrm{2}}\\ \frac{q_{\mathrm{5}}}{\mathrm{2}}& \frac{q_{\mathrm{6}}}{\mathrm{2}}& q_{\mathrm{3}}\end{array}\right]\mathrm{。}（16）

未知参数q_{\mathrm{3}}的求解可通过构造非线性残差方程， 将TTLS解与椭球一致性条件结合， 形成关于q_{\mathrm{3}}的标量方程。再利用数值迭代求解直至残差收敛， 从而得到稳定的q_{\mathrm{3}}估计。

经由TTLS估计得到二次项系数后， 便可构造完整的软磁效应矩阵\mathit{S}与硬磁效应向量\mathit{H}， 进而获取补偿矩阵\mathit{M}， 其具体过程如下。

首先， 对矩阵{\mathit{A}}_{\mathrm{e}}进行奇异值分解

{\mathit{A}}_{\mathrm{e}}=\mathit{J}\mathit{K}{\mathit{J}}^{\mathrm{T}},（17）

式中： \mathit{J}为正交矩阵； \mathit{K}为奇异值对角阵。利用其主轴方向与尺度信息， 可重构补偿矩阵

\widehat{\mathit{M}}=F\mathit{U}\mathit{J}\sqrt[]{\mathit{K}}{\mathit{J}}^{\mathrm{T}},（18）

式中： F为当地地磁场强度， 可通过IGRF获取或在无人车上加装伸杆并安装传感器测量。进而， 软磁效应矩阵与硬磁效应向量可表示为

\widehat{\mathit{S}}={\mathit{M}}^{-\mathrm{1}}-{\mathit{I}}_{\mathrm{3}},\widehat{\mathit{H}}=\mathit{M}{\mathit{x}}_{\mathrm{0}}，（19）

式中： {\mathit{I}}_{\mathrm{3}}为\mathrm{3}\times \mathrm{3}单位矩阵。上述参数共同描述了磁传感器的畸变特性： \widehat{\mathit{S}}表征比例因子差异与非正交性引起的软磁效应； \widehat{\mathit{H}}反映由零偏和结构不对称导致的硬磁效应； \widehat{\mathit{M}}则表征软、 硬磁畸变对观测磁场的整体影响， 是实现磁补偿的核心参数。

对于任意一次原始测量磁矢量

{\mathit{b}}_{\mathrm{T}}=\left[\begin{array}{c}{m}_x\\ m_y\\ m_z\end{array}\right],（20）

其补偿后的三分量可统一表示为

{\mathit{b}}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}=\widehat{\mathit{M}}\left(\begin{array}{c}{\mathit{b}}_{\mathrm{T}}-\widehat{\mathit{H}}\end{array}\right)\mathrm{。}（21）

即对磁测原始数据{\mathit{b}}_{\mathrm{T}}首先减去硬磁效应向量\widehat{\mathit{H}}进行平移补偿， 再通过矩阵\widehat{\mathit{M}}矫正各类尺度与方向畸变， 进而获取补偿磁场。

1.3 t型抗差估计融合

1.3.1 基于t分布的误差建模

TTLS具有良好的数值稳定性， 然而， 其仍依赖于残差分布较接近高斯的假设。针对无人车磁测数据中普遍存在的非高斯噪声、 重尾残差及局部粗差问题， 传统TTLS方法在高噪声、 强异常观测条件下仍易出现参数估计失真。为提高对各类噪声/粗差的鲁棒性， 本文引入具有长尾特性的t分布对观测误差进行建模。t分布能够在统计意义上降低大残差样本的影响权重， 将其视为低置信度观测， 为后续稳健估计提供合理的概率基础。

t型估计与TTLS融合的核心在于利用t模型导出的权重重构加权增广矩阵， 并执行加权TTLS估计， 从而实现截断总体最小二乘与重尾分布建模结合的稳健估计。

记第i个残差为

r_i={\mathit{a}}_i^{\mathrm{T}}\mathit{\varphi }-b_{\mathrm{T},i},（22）

式中： {\mathit{a}}_i^{\mathrm{T}}为\mathit{A}的第i行； b_{\mathrm{T},i}为{\mathit{b}}_{\mathrm{T}}的第i个元素。若假设残差服从t分布r_i\sim t_{\nu }(\mathrm{0},{\sigma }^{\mathrm{2}})， 其参数化形式的概率密度为

f(r_i\left|\nu ,\sigma \right.)=\frac{\Gamma \left(\frac{\nu +\mathrm{1}}{\mathrm{2}}\right)}{\Gamma \left(\frac{\nu }{\mathrm{2}}\right)\sqrt[]{\nu \mathrm{\pi }}\sigma }{\left(\mathrm{1}+\frac{r_i^{\mathrm{2}}}{\nu {\sigma }^{\mathrm{2}}}\right)}^{-\frac{\nu +\mathrm{1}}{\mathrm{2}}}\mathrm{。}（23）

式（23）所描述的分布中， 其尾部概率随离群值的增加呈幂次级下降， 相较于高斯分布的指数级衰减， 这种缓变的尾部结构显著增强了模型对大残差的鲁棒性， 体现出其分布的重尾特性。

t分布可表达为高斯分布与随机精度变量{\lambda }_i的混合， 其条件分布为

r_i\left|{\lambda }_i\sim \mathrm{𝒩}\right.\left(\mathrm{0},\frac{{\sigma }^{\mathrm{2}}}{{\lambda }_i}\right)\mathrm{。}（24）

即残差在给定{\lambda }_i下是一个方差被缩放的高斯分布， 其潜变量分布为

{\lambda }_i\sim Gamma\left(\frac{\nu }{\mathrm{2}},\frac{\nu }{\mathrm{2}}\right),（25）

式中： Gamma(\alpha ,\beta )的密度为

p(\lambda )=\frac{{\beta }^{\alpha }}{\Gamma (\alpha )}{\lambda }^{\alpha -\mathrm{1}}{\mathrm{e}}^{-\beta \lambda },\hspace{1em}\lambda >\mathrm{0}\mathrm{。}（26）

联合式（24）与式（25）可得

p(r_i)={\int }_{\mathrm{0}}^{\mathrm{\infty }}p(r_i\left|{\lambda }_i)p({\lambda }_i)\mathrm{d}{\lambda }_i\right.\mathrm{。}（27）

式（27）的积分结果为式（23）的t分布密度。由上述推导可知， t分布可被解释为一种带有随机缩放因子的高斯层次模型， 每个观测的方差由随机变量调节， 从而产生长尾特性。

在此表示下， t分布等价于加权模型， 每个样本对应一个隐含权重， 反映其置信度， 通过后验期望计算可获得“有效权重”， 迭代中自动削弱异常观测的影响。

1.3.2 基于EM框架的鲁棒求解策略

在引入t分布误差模型后， 其关键在于如何将概率建模结果有效地融入TTLS框架中。根据1.3.1节结论， 每个观测样本可对应一个反映其置信度的隐含权重， 本节通过对该权重的后验期望进行估计， 将t型抗差机制转化为加权形式， 从而构建与TTLS结构兼容的加权增广矩阵， 实现鲁棒建模与病态抑制的统一。

在含大量观测误差的磁测数据中， 直接对t分布似然极大化计算复杂且依赖初值。EM算法利用t分布的随机权重特性， 将优化过程拆为E步更新观测权重与M步加权TTLS拟合， 在保证数值稳健的同时， 具备对TTLS矩阵求解的兼容性。

在EM框架中， 将潜变量{\lambda }_i视作未观测的精度因子。

E步： 根据当前残差， 计算每个观测点的潜变量{\lambda }_i的后验期望， 从而量化各点的置信度。根据式（22）， 记第i个残差为r_i， 其后验分布为

{\lambda }_i\left|r_i\right.\sim Gamma\left(\frac{\nu +\mathrm{1}}{\mathrm{2}},\frac{\nu +r_i^{\mathrm{2}}/{\sigma }^{\mathrm{2}}}{\mathrm{2}}\right)\mathrm{。}（28）

其期望为

w_i\equiv \mathbb{E}[{\lambda }_i\left|r_i]\right.=\frac{\nu +\mathrm{1}}{\nu +\frac{r_i^{\mathrm{2}}}{{\sigma }^{\mathrm{2}}}}\mathrm{。}（29）

故残差越大， 式（29）中的权重w_i越小， 即对应粗差点被抑制。

1.3.3 基于IRLS的迭代加权TTLS

在EM框架下， M步的参数更新可等价表示为一个带权最小二乘问题。通过将E步得到的观测权重引入TTLS增广矩阵， 本文将鲁棒估计过程转化为基于IRLS的迭代加权TTLS求解。该方法在不改变TTLS原有奇异值分解结构的前提下， 实现了权重自适应更新与参数稳健估计的统一。

M步： 在给定权重的条件下， 执行加权TTLS拟合以更新参数向量\mathit{\varphi }， 同时依据残差的加权平方和更新尺度参数\sigma。

固定权重w_i与尺度\sigma， 更新参数向量\mathit{\varphi }。目标函数为

Q(\mathit{\varphi },\sigma )={\displaystyle \sum _{i=\mathrm{1}}^N}\frac{{\omega }_i}{\mathrm{2}{\sigma }^{\mathrm{2}}}{\left({\mathit{a}}_i^{\mathrm{T}}\mathit{\varphi }-b_{\mathrm{T},i}\right)}^{\mathrm{2}}+N\mathrm{l}\mathrm{n}\sigma \mathrm{。}（30）

当\sigma固定时， 式（30）等价为一个加权最小二乘问题

\underset{\mathit{\Phi }}{\mathrm{m}\mathrm{i}\mathrm{n}}{\displaystyle \sum _{i=\mathrm{1}}^N}{w}_i{\left({\mathit{a}}_i^{\mathrm{T}}\mathit{\varphi }-b_{\mathrm{T},i}\right)}^{\mathrm{2}}\mathrm{。}（31）

式（31）的w_i即为E步中式（29）所得， 其反映第i个观测点的“置信度”， 残差大的点对应较小的w_i， 因此， 在后续拟合中对其影响被削弱， 实现抑制粗差的效果。

为实现与TTLS方法的融合， 将上述权重以矩阵形式进行表征参与后续求解。

\mathit{W}=diag(\sqrt[]{w_{\mathrm{1}}},\sqrt[]{w_{\mathrm{2}}},\dots ,\sqrt[]{w_N})\mathrm{。}（32）

定义加权后的设计矩阵与观测向量为\tilde{\mathit{A}}=\mathit{W}\mathit{A},\tilde{\mathit{b}}=\mathit{W}{\mathit{b}}_{}， 则式（31）可写作

\underset{\mathit{\Phi }}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖\tilde{\mathit{A}}\varphi -\tilde{\mathit{b}}‖}_{\mathrm{2}}^{\mathrm{2}}\mathrm{。}（33）

后续， 保持TTLS的求解逻辑， 对增广矩阵\tilde{\mathit{G}}=[\tilde{\mathit{A}},\tilde{\mathit{b}}]进行奇异值分解， 并根据截断参数k对右奇异向量分块， 得到更新后的参数向量为

{\mathit{\varphi }}^{(t+\mathrm{1})}=-{\mathit{V}}_{\mathrm{12}}{\mathit{V}}_{\mathrm{22}}^{†},（34）

式中： {\mathit{V}}_{\mathrm{12}}与{\mathit{V}}_{\mathrm{22}}为\tilde{\mathit{G}}的SVD子块， 其分块形式可参考式（13）。在此基础上利用残差的加权平方和更新尺度，

{\left({\sigma }^{(t+\mathrm{1})}\right)}^{\mathrm{2}}=\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}{w}_i^{(t+\mathrm{1})}{\left(r_i^{(t+\mathrm{1})}\right)}^{\mathrm{2}}\mathrm{。}（35）

上述迭代过程将持续进行， 直到参数收敛或相对变化量小于阈值， 最终得到稳定的\varphi、 \sigma与权重分布。

同时， 为保证算法的数值稳定性与计算收敛， 设置收敛判据： 当连续两次迭代间的参数相对变化率满足式（36）时， 认为参数收敛。其中， \eta是为避免除零错误而设置的最小阈值。

（36）

后续将收敛的参数向量还原为椭球矩阵与中心向量， 基于此构建软磁效应矩阵\mathit{S}与硬磁效应向量\mathit{H}， 进而获得补偿矩阵\mathit{M}\in {\mathbf{R}}^{\mathrm{3}\times \mathrm{3}}， 实现对三轴磁通门测量的稳健补偿。本文算法流程如图 2 所示。

1.3.4 自由度参数的自适应估计

在t分布鲁棒建模中， 自由度v决定了分布的尾部厚度， 是控制抗差能力的关键因素。为避免经验设定自由度v所带来的主观性， 本文引入一种基于初始残差统计特性的自由度自适应估计策略， 使模型能够根据实际数据分布特征自动调节鲁棒强度。

首先， 利用标准TTLS方法获取初始参数估计， 并构造对应残差向量\mathit{r}=\{r_i{\}}_{i=\mathrm{1}}^N。假设残差服从t分布， 其概率密度如式（23）。通过最大似然准则对自由度参数v进行估计， 从而使所得t分布在统计意义上最优地刻画初始残差的重尾特性。该过程不依赖先验经验参数， 实现了对自由度的自动估计。在后续迭代过程中， 通过残差驱动权重更新以反映异常观测的影响程度。该策略使算法能够依据残差分布特性自适应调节抗差强度， 在抑制异常观测影响的同时保持对有效数据的充分利用， 实现鲁棒性与适应性的兼顾。

相比传统最小二乘或TLS方法， TTLS-t具备两大优势： 通过L曲线自适应截断奇异值， 抑制病态矩阵与噪声放大； 引入t型误差模型， 将残差映射为动态权重， 削弱噪声与粗差影响。通过自由度的调节， 算法可兼顾高斯噪声下的效率与重尾分布下的稳健性， 实现对复杂电磁环境中磁测数据的可靠、 适应性补偿。

2 实验验证

实验选用无人车进行载体磁干扰补偿算法性能验证。无人车安装有数据采集模组、 三轴磁通门、 光泵磁力仪以及GNSS， 如图 3 所示。数据采样率为100 Hz。

测试地点选定为一空旷步道路段， 实验区域占地面积约100 m×10 m， 当地磁场环境稳定。

实验采用匀速圆周运动作为标定工况。通过在二维场景限制下进行各自由度机动的遍历， 以充分激励补偿参数。无人车在程控模式下平稳运行， 测试区域远离干扰源， 磁场环境稳定， 数据理想。该组数据记为A， 其轨迹与三轴磁通门数据如图 4、 图 5 所示。

为分析无人车磁测数据的统计特性与误差分布， 首先对原始观测进行分布拟合。受姿态变化、 车体磁扰及环境扰动影响， 数据含非高斯噪声与粗差， 传统高斯假设难以准确表征。通过绘制概率密度直方图并对比正态与t分布拟合， 验证采用t分布建模的合理性。

以无人车匀速直线行驶为例验证分布拟合情况。实验在磁环境稳定区域进行， 无人车以0.5 m/s速度匀速行驶765 s， 沿100 m道路折返两次， 轨迹如图 6 所示。受多源干扰影响， 测得磁场总场数据含大量观测误差， 如图 7 所示。

对该组磁测数据的残差序列分别进行正态分布与t分布拟合。经验证， t分布能够更准确地描述干扰特性， 如图 8 所示。由此说明实测噪声更符合t分布， 为TTLS-t基于t型误差模型的鲁棒补偿提供数据支撑。

2.1 扰动仿真验证

为验证TTLS-t在无人车磁补偿中的鲁棒性与数值稳定性， 本节在实测纯净数据上注入不同自由度的t型粗差， 构造多类型干扰样本， 模拟随机噪声、 突发粗差及系统扰动场景。通过统一基准对比试验， 评估算法在复杂观测下的稳定性与异常抑制能力。

在数据A的三分量磁测数据基础上仿真添加自由度v=\mathrm{5}的t型粗差， 将其作为标定数据， 记为A_t5。数据A_t5在注入粗差前后的磁场强度对比如图 9 所示。

采用数据A_t5通过TTLS-t计算补偿参数验证可用性。因标定区磁场稳定， 补偿磁场强度应接近常值， 故以标准差评估补偿效果， 更小的标准差代表更优的补偿效果。

经验证， 标定数据自身补偿后模值标准差为19.03 nT， 相较原始数据整体改善约95.6%， 补偿效果明显。相较传统代表性方法EF与TTLS， 本文提出的TTLS-t方法在观测误差扰动下的磁补偿精度提升大于89%， 优化效果显著， 表现出良好的鲁棒性与数值稳定性。具体结果如表 1 与图 10 所示。

EF在无人车磁测中效果有限， 因三维姿态变化受限、 Z轴自由度不足， 观测点分布接近扁平椭球而非理想球面， 导致拟合偏差大。此外， EF基于最小二乘， 对异常点与粗差敏感， 观测误差易扭曲拟合结果。

TTLS通过奇异值分解与截断抑制部分噪声， 并改善矩阵病态问题。但其基于高斯噪声假设， 对长尾误差和异常观测缺乏抑制， 观测误差过多时会导致参数估计偏移， 补偿精度下降。

TTLS-t结合TTLS的数值稳定性与t分布的抗差特性， 可在噪声、 观测误差及病态数据下稳健发挥。通过截断奇异值抑制病态， 并以迭代加权调整观测权重， 有效抑制异常点， 实现多类型误差的稳健补偿， 性能优于EF与TTLS。

为验证算法对各类观测误差的处理能力， 基于数据A分别加入自由度3~9的t型粗差， 生成A_t3~A_t9的7组数据。利用t分布自由度模拟不同干扰， 评估算法鲁棒性与数值稳定性。以A_t5计算补偿参数， 并对7组数据分别用EF、 TTLS和TTLS-t补偿， 结果见表 1。

经多组数据验证， EF和TTLS在异常观测下易偏移或发散， 补偿效果有限； TTLS-t能有效抑制异常点影响， 补偿标准差均小于20 nT， 相较现有方法精度提升均大于89%， 稳健性与精度显著优越。

2.2 实测数据验证

为进一步评估算法应对地面环境下复杂、 不可控误差的能力， 采用无人车实测数据， 以“8”字型与“S”型轨迹为例， 验证算法在真实条件下的适用性、 稳定性及补偿精度。

通过人工操控， 使无人车在设备附近磁干扰显著区域测试， 分别进行“8”字型轨迹与“S”型轨迹机动， 行驶轨迹如图 11 所示。

测试过程中受场地、 设备及人工操作等多重干扰影响， 磁测数据出现大量观测误差， 导致磁场测量不稳定， 其模值如图 12 所示。

同扰动仿真验证中所用方法， 分别以EF、 TTLS和TTLS-t对实测数据进行磁干扰补偿， 并对3种结果进行对比分析， 具体结果见表 2。

图 13 与图 14 分别为实测验证“S”型轨迹与“8”字型轨迹磁场补偿效果对比。实测表明， 当数据存在大量异常观测点或粗差时， EF和TTLS因最小二乘敏感性高， 补偿参数易偏移， 导致补偿磁场分布离散或补偿失效（如“8”字型轨迹补偿结果）。相比之下， TTLS-t在TTLS基础上引入t型抗差估计， 通过迭代权重抑制异常点， 有效减弱了重尾残差影响。在复杂磁环境下， TTLS-t保持数值稳定的同时， 补偿后磁场模值更集中、 标准差更小， 补偿精度较传统方法提升约75%~89%， 整体性能显著优越。

3 结　论

本文针对无人车实测磁场数据同时存在各类观测误差、 粗差及观测病态的问题， 提出TTLS-t抗差估计补偿算法， 实现模型驱动下的稳健磁补偿。该方法将t型抗差估计融入TTLS， 首先以IRLS形式实现迭代加权， 将残差映射为动态权重， 以抑制各类观测误差、 粗差及重尾观测， 然后， 基于EM算法框架实现稳健求解， 与此同时， 保持TTLS自身的数值稳定性和病态抑制能力。 算法融合二者优势在多源干扰与异常观测条件下能够显著削弱异常数据的不良影响。多类实验验证表明， TTLS-t在数据质量不理想、 含大量观测误差与病态问题的情况下仍能实现高精度、 强鲁棒的磁补偿， 整体性能优于TTLS和EF， 为无人车地磁导航及组合导航提供了一种稳健的建模策略。

后续研究中， 该方法可扩展至空中与水下地磁导航场景。但因其主要针对地面观测误差与病态问题进行优化， 在数据质量理想、 载体机动充分的海空测试条件下， 本文提出方法的改进幅度可能有限。

参考文献

[1] 林沂， 孙晶京， 闫旭. 地磁导航定位技术原理与方法综述［J］. 全球定位系统， 2023， 48（6）： 32-41.

[2] Lin Yi， Sun Jingjing， Yan Xu. A review of the principles and methods of geomagnetic navigation and positioning technology［J］. GNSS World of China， 2023， 48（6）： 32-41. （in Chinese）

[3] 张樱子， 刘改改， 申雅丽， 等. 基于归一化磁源强度的磁目标实时定位方法［J］. 测试技术学报， 2024， 38（5）： 535-542.

[4] Zhang Yingzi， Liu Gaigai， Shen Yali， et al. Real-time magnetic target localization method based on normalized source strength［J］. Journal of Test and Measurement Technology， 2024， 38（5）： 535-542. （in Chinese）

[5] Canciani A， Raquet J. Airborne magnetic anomaly navigation［J］. IEEE Transactions on Aerospace and Electronic Systems， 2017， 53（1）： 67-80.

[6] Moradi M， Zhai Z， Nielsen A， et al. Random forests for detecting weak signals and extracting physical information： a case study of magnetic navigation［J］. APL Machine Learning， 2024， 2： 016118.

[7] Chun Y， Chen G W， Si Y B， et al. Improved strong tracking Kalman filter algorithm based SINS/GNSS/ODO integrated navigation［J］.Journal of Measurement Science and Instrumentation， 2026，17（1）： 61-71.

[8] Li H， Liu M， Zhang F. Geomagnetic navigation of autonomous underwater vehicle based on multi-objective evolutionary algorithm［J］. Frontiers in Neurorobotics， 2017， 11： 34.

[9] Chen K， Liang W， Liu M， et al. Comparison of geomagnetic aided navigation algorithms for hypersonic vehicles［J］. Journal of Zhejiang University： Science A， 2020， 21（8）： 673-683.

[10] 王劲楠， 张晓明， 刘俊， 等. 基于矢量磁传感器模型的扩充T-L航磁补偿方法［J］. 导航定位与授时， 2025， 12（2）： 78-85.

[11] Wang Jinnan， Zhang Xiaoming， Liu Jun， et al. An extended T-L aeromagnetic compensation method based on vector magnetic sensor model［J］. Navigation Positioning and Timing， 2025， 12（2）： 78-85. （in Chinese）

[12] 李志敏， 赖际舟， 吕品， 等. 基于椭球拟合的三轴磁传感器快速标定实验设计［J］. 实验技术与管理， 2023， 40（10）： 42-47.

[13] Li Zhimin， Lai Jizhou， Pin Lü， et al. Experimental design for rapid calibration of three-axis magnetic sensors based on ellipsoidal fitting［J］. Experimental Technology and Management， 2023， 40（10）： 42-47. （in Chinese）

[14] Shakerian A， Bilel S， Jr Landry R. Calibration of 3-axis low-cost magnetometer using the least square ellipsoid fitting algorithm［J］. International Journal of Sensors and Sensor Networks， 2023， 11（1）： 18-24

[15] Chi C， Wang D， Tao R， et al. Error calibration of cross magnetic gradient tensor system with total least-squares method［J］. Mathematical Problems in Engineering， 2023（1）： 6974834.

[16] 孙同贺. 病态与稳健总体最小二乘平差方法若干关键问题研究［J］. 测绘学报， 2020， 49（3）： 398.

[17] Sun Tonghe. Research on several key problems of ill-posed and robust total least squares method［J］. Acta Geodaetica et Cartographica Sinica， 2020， 49（3）： 398. （in Chinese）

[18] Ji C， Song C， Li S， et al. A compensation method of carrier magnetic interference under pathological conditions for geomagnetic navigation［J］. Journal of Applied Physics， 2024， 135（5）： 054901.

[19] Ning Z， Jiang P， Pan X， et al. Compensation technology of vehicle magnetometer based on TTLS-Tikhonov regularization［J］. IEEE Sensors Journal， 2024， 24（22）： 36594-36603.

[20] Adomaityte U， Defilippis L， Loureiro B， et al. High-dimensional robust regression under heavy-tailed data： asymptotics and universality［J］. Journal of Statistical Mechanics： Theory and Experiment， 2024（11）： 114002.

[21] Wu J， Jiang J， Tang Y， et al. Gaussian-student’s t mixture distribution-based robust Kalman filter for global navigation satellite system/inertial navigation system/odometer data fusion［J］. Remote Sensing， 2024， 16（24）： 4716.

[22] Liu B， Liu Q， Wu Z. A novel robust Student’s t scale mixture distribution based Kalman filter［J］. Signal Processing， 2026， 239： 110296.

[23] Zhou S， Yang C， Su Z， et al. An aeromagnetic compensation algorithm based on radial basis function artificial neural network［J］. Applied Sciences， 2022， 13（1）： 136.

[24] Chen Z， Yu Z， Liu C， et al. Intelligent dynamic-enhanced compensation for UAV magnetic interference［J］. Sensors， 2025， 25（16）： 5059.

[25] 刘强， 王涛， 张晓明， 等. Tolles-Lawson模型—一维卷积神经网络磁补偿方法［J］. 导航定位与授时， 2025， 12（2）： 69-77.

[26] Liu Qiang， Wang Tao， Zhang Xiaoming， et al. A one-dimentional convolutional neural network aeromagnetic compensation method based on Tolles-Lawson model［J］. Navigation Positioning and Timing， 2025， 12（2）： 69-77. （in Chinese）

[27] 王劲楠， 张晓明， 张鸽， 等. 基于GAN-LSTM改进T-L磁补偿方法［J］. 电子测量技术， 2025， 48（5）： 50-56.

[28] Wang Jinnan， Zhang Xiaoming， Zhang Ge， et al. Improved T-L magnetic compensation method based on GAN-LSTM［J］. Electronic Measurement Technology， 2025， 48（5）： 50-56. （in Chinese）

[29] Xu Y， Liu Z， Zhang Q， et al. Interference model guided neural network for aeromagnetic compensation［J］. IEEE Sensors Journal， 2024， 24（8）： 12266-12275.

[30] Zhang G， Zhang X， Gao L， et al. Real-time attitude estimation for spinning projectiles by magnetometer based on an adaptive extended Kalman filter［J］. Micromachines， 2023， 14（11）： 2000.

[31] 杜江波， 张晓明， 薛晓华. 稳健估计下三轴磁传感器误差补偿方法研究［J］. 传感器与微系统， 2025， 44（1）： 29-33.

[32] Du Jiangbo， Zhang Xiaoming， Xue Xiaohua. Research of error compensation method of three-axis magnetic sensor by robust estimation［J］. Transducer and Microsystem Technologies， 2025， 44（1）： 29-33. （in Chinese）

[33] Huang C， Hendeby G， Skog I. Joint magnetometer-IMU calibration via maximum a posteriori estimation［EB/OL］. ［2005-10-23］.

[34] RodríGuez-MartíNez S， Troni G. Towards a factor graph-based method using angular rates for full magnetometer calibration and gyroscope bias estimation［C］//2024 IEEE/RSJ International Conference on Intelligent Robots and Systems （IROS）， 2024： 1199-1205.

基金资助

国家自然科学基金资助项目(61873247)