地磁传感器标定技术综述

一、地磁传感器介绍

1、 地磁传感器概念：

地磁传感器，又称磁强计、磁力计、电子罗盘，是一种通过测量地磁场的磁感应强度来确定方位的传感器。
在MEMS技术大力发展下，地磁传感器广泛应用于与MEMS 陀螺仪、加速度计进行数据融合，组成AHRS（航姿参考系统 Attitude and heading reference system），从而以较少的成本获得准确的姿态。应用领域有手机、无人机、路基惯性/GPS组合导航、人体动作捕捉等。

2、地磁传感器误差分析：

理论上，地磁传感器的输出在空间上的分布应该是一个球心位于原点的圆球面。由于地磁传感器自身以及所处环境的影响，地磁传感器的实际输出通常为一个偏离原点的椭球面。

地磁传感器的误差主要分为两类： 1）传感器结构、材料和电路引起的误差，如非正交误差、刻度因子误差和零偏误差等； 2）外部干扰磁场的叠加影响，如固定磁场、感应磁场等。

对地磁传感器输出进行建模，并考虑安装误差、非正交误差、刻度因子误差、硬磁误差和软磁误差后，将输出模型进行简化，可以得到输出模型为：¹ $${\hat{B}}^b=A_M{C}_e^b{B}_e^e+b^b+{\omega }^b\text{\hspace{1em}(1)}$$ $A_M$为包含安装误差、非正交误差、刻度因子误差和软磁系数的误差阵；$b^b$为包含误差阵$A_M$、硬磁误差和零偏误差的误差矢量，${\omega }^b$是测量随机误差和载体电路系统产生的杂散磁场造成的噪声和，认为是0均值白噪声。

暂不考虑${\omega }^b$，对上式移项取逆后可得地磁传感器输出校正模型： $${\hat{B}}_e^e=C_b^e{A}_M^{-1}({\hat{B}}^b-b^b)\text{\hspace{1em}(2)}$$ 其中${\hat{B}}_e^e$为地磁场矢量的估计值，$C_b^e$为从传感器坐标系到地理坐标系的方向余弦矩阵，且有$C_b^e=(C_e^b{)}^{-1}$。

二、地磁传感器标定方法

地磁传感器标定方法主要分为三类，一类为自标定，仅仅利用此传感器自身采集的数据进行分析，得到各项误差系数；第二类为利用陀螺仪/加速度计等其他传感器辅助标定；第三类为利用无磁转台进行基准比较²。本文主要介绍前两类方法。

1.自标定

地磁传感器自标定的方法主要依赖于地磁矢量的模值大小在当地保持基本不变。

在公式(2)中，单独使用地磁传感器无法得到$C_b^e$，故对公式等号两侧取模可得：
$$||{\hat{B}}_e^e||=||C_b^e{A}_M^{-1}({\hat{B}}^b-b^b)||=||A_M^{-1}({\hat{B}}^b-b^b)||\text{\hspace{1em}(3)}$$ 如此一来，便摆脱了对$C_b^e$的依赖，在无法获得地磁传感器姿态的情况下，可利用式（3）进行标定补偿。 对矩阵$A_M^{-1}$进行QR分解，可得 $$A_M^{-1}=Q_M{R}_A\text{\hspace{1em}(4)}$$ 其中，显然$A_M^{-1}$非奇异，故分解唯一，$Q_M$为正交矩阵，$R_A$为上三角阵，所有对角线元素为正。 将（4）式带入（3）式，可得 $$||{\hat{B}}_e^e||=||R_A({\hat{B}}^b-b^b)||\text{\hspace{1em}(4)}$$ 于是，参数空间的维数变为$R_A$的六维加$b^b$的3维共9维。由于QR分解唯一性，保证了参数的唯一性，从而避免了（3）式中由于丢失姿态信息以及求模操作而导致的多值问题。但由于缺失正交阵$Q_M$，按照（4）进行标定的地磁传感器输出相比理想标定输出相差一个姿态变换，这可以通过引入外部传感器例如陀螺仪/加速度计解决¹。

对于式（4），有 $$R_0^2=||R_A({\hat{B}}_i^b-b^b)|{|}^2,i=1,2,...,m\text{\hspace{1em}(5)}$$ 其中，$R_0^2$为当地地磁场矢量的模值；${\hat{B}}_i^b$为传感器输出矢量；$i$为测量序列次数；m为采样个数。 将上式展开，可得， $$a_1{x}_i^2+a_2{x}_i{y}_i+a_3{x}_i{z}_i+a_4{y}_i^2+a_5{y}_i{z}_i+a_6{z}_i^2+a_7{x}_i+a_8{y}_i+a_9{z}_i+a_{10}=0,i=1,2,...,m,\text{\hspace{1em}(6)}$$ 式中：$x_i,y_i,z_i$分别为${\hat{B}}_i^b$的3个分量，$a_1-a_{1}0$是由参数$R_A$和$b^b$构成的联合变量。 将(6)式写成矩阵形式，有 $$HX=W,\text{\hspace{1em}(7)}$$ 式中： $$X=\left[\begin{array}{ccccccccc}{x}_1^2& x_1{y}_1& x_1{z}_1& y_1^2& y_1{z}_1& x_1& y_1& z_1& 1\\ x_2^2& x_2{y}_2& x_2{z}_2& y_2^2& y_2{z}_2& x_2& y_2& z_2& 1\\ ...& ...& & & & & & & \\ x_m^2& x_m{y}_m& x_m{z}_m& y_m^2& y_m{z}_m& x_m& y_m& z_m& 1\end{array}\right]$$ $$X={\left[\begin{array}{ccccccccc}{a}_1/a_6& a_2/a_6& a_3/a_6& a_4/a_6& a_5/a_6& a_7/a_6& a_8/a_6& a_9/a_6& a_{10}/a_6\end{array}\right]}^T$$ $$W={\left[\begin{array}{cccc}-z_1^2& -z_2^2& ...& -z_m^2\end{array}\right]}^T$$ $H\in {\Re }^{m\times n},X\in {\Re }^{n\times 1},W\in {\Re }^{m\times 1}$;一般来说，$m>n$,因此(7)式为超定方程¹。 以上理论推导参考文献¹。

下面介绍在求得X之后，如何得到$R_A,b^b$。
由向量2范数的定义$||X|{|}^2=X^{T}X$，(5)式可展开为：
$$R_0^2=[R_A({\hat{B}}_i^b-b^b){]}^T[R_A({\hat{B}}_i^b-b^b)]=({{\hat{B}}_i^b}^T-{b^b}^T){R_A}^T{R}_A({\hat{B}}_i^b-b^b),i=1,2,...,m.\text{\hspace{1em}(9)}$$
将上式继续展开，可得
$$R_0^2={{\hat{B}}_i^b}^T{R_A}^T{R}_A{\hat{B}}_i^b-{{\hat{B}}_i^b}^T{R_A}^T{R}_A{b}^b-{b^b}^T{R_A}^T{R}_A{\hat{B}}_i^b+{b^b}^T{R_A}^T{R}_A{b}^b,i=1,2,...,m.\text{\hspace{1em}(10)}$$
由于${\hat{B}}_i^b\in {\Re }^{3\times 1},R_A\in {\Re }^{3\times 3},b^b\in {\Re }^{3\times 1}$上式中${{\hat{B}}_i^b}^T{R_A}^T{R}_A{b}^b$为标量，对其取转置值不变即
$${{\hat{B}}_i^b}^T{R_A}^T{R}_A{b}^b=({{\hat{B}}_i^b}^T{R_A}^T{R}_A{b}^b{)}^T={b^b}^T{R_A}^T{R}_A{\hat{B}}_i^b，i=1,2,...,m.\text{\hspace{1em}(11)}$$
于是，（10）式可简化为：
$$R_0^2={{\hat{B}}_i^b}^T{R_A}^T{R}_A{\hat{B}}_i^b-2{b^b}^T{R_A}^T{R}_A{\hat{B}}_i^b+{b^b}^T{R_A}^T{R}_A{b}^b,i=1,2,...,m.\text{\hspace{1em}(12)}$$
其中，第一项为二次型，第二项为一次项，最后一项为常数项。
显然${R_A}^T{R}_A$为对称阵，令${R_A}^T{R}_A=M=\left[\begin{array}{ccc}{m}_1& m_2& m_3\\ m_2& m_4& m_5\\ m_3& m_5& 1\end{array}\right]，{\hat{B}}_i^b={\left[\begin{array}{ccc}x& y& z\end{array}\right]}^T$,则
$$\begin{gathered} {{\hat{B}}_i^b}^T{R_A}^T{R}_A={{\hat{B}}_i^b}^{T}M{\hat{B}}_i^b=\left[\begin{array}{ccc}x& y& z\end{array}\right]\left[\begin{array}{ccc}{m}_1& m_2& m_3\\ m_2& m_4& m_5\\ m_3& m_5& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right] \\ =m_1{x}^2+m_4{y}^2+z^2+2m_{2}xy+2m_{3}xz+2m_{5}yz.\text{\hspace{1em}(13)} \end{gathered}$$
对比上式与（7）式，可得
$$M=\left[\begin{array}{ccc}{m}_1& m_2& m_3\\ m_2& m_4& m_5\\ m_3& m_5& 1\end{array}\right]=\left[\begin{array}{ccc}X(1)& X(2)/2& X(3)/2\\ X(2)/2& X(4)/2& X(5)/2\\ X(3)/2& X(5)/2& 1\end{array}\right]\text{\hspace{1em}(14)}$$
对$M$取乔伊斯基变换可得$R_A=[cholesky(M){]}^T$
另外，对于一次项，可采用类似比较系数的方法，得
$$-2{b^b}^T{R_A}^T{R}_A=\left[\begin{array}{ccc}X(6)& X(7)& X(8)\end{array}\right]\text{\hspace{1em}(15)}$$
又${R_A}^T{R}_A=M$，所以有
$${b^b}^T=-\frac{1}{2}\left[\begin{array}{ccc}X(6)& X(7)& X(8)\end{array}\right]M^{-1}\text{\hspace{1em}(16)}$$

注意：在上述推导中，将Z轴上的分量的大小作为标定后地磁矢量的模值大小，显得不那么严谨。好在磁强计在使用时是用比值进行计算角度，因此问题不大。更为准确的做法我认为是两边同除以模值,让（7）式中a10不必移项，直接变为1即可。标定完毕后，根据位置和时间查询地磁矢量的模值大小（MatLab 中函数igrfmagm和wrldmagm，二者所用模型不同，暂未查询到二者的区别）。

1.1.两步法

第一步： 对于（7）式，假定矩阵H是精确的，W含有噪声，根据Gauss-markov定理，X的最小二乘估计（LS）为 $$X_{LS}=(H^{T}H{)}^{-1}{H}^{T}W,\text{\hspace{1em}(17)}$$ 第二步： 根据算法X中间量的定义，求出$R_A$和$b^b$。

两步法中，H和W均由地磁传感器量测值构造而成，二者都含有噪声，当测量值误差较大值，两步法标定较差。

1.2.基于总体最小二乘（TLS）的算法

基于TLS的方法利用扰动矩阵和扰动向量同时干扰H和W，以便矫正二者内在的扰动，由于考虑到了测量误差，理论上其估计精度比LS更高。 假设H和W含有测量误差，分别记为E和e，（7）式可改写为： $$(H+W)X=W+e\text{\hspace{1em}(18)}$$ TLS将其转化为约束最优化问题进行求解¹，

最优解可通过对扩展矩阵$G=\left[\begin{array}{cc}-W& H\end{array}\right]$进行奇异值分解得到¹

1.3.递推最小二乘法（RLS)

递推最小二乘法主要将批量处理地磁传感器数据改为递推形式，更适合在线执行³。
文献³所述的参数方程与上面推导的（6）式略有不同，但也只是系数归一化的所参照的参数项不同而已。
对（6）式稍做变换，写成 矩阵形式，有
$$h(k)=\left[\begin{array}{ccccccccc}{x}_i^2& x_i{y}_i& x_i{z}_i& y_i^2& y_i{z}_i& x_i& y_i& z_i& 1\end{array}\right]$$
$$z(k)=\left[\begin{array}{c}-z_i^2\end{array}\right]$$
$$\theta (k)=\left[\begin{array}{ccccccccc}{a}_1/a_6& a_2/a_6& a_3/a_6& a_4/a_6& a_5/a_6& a_7/a_6& a_8/a_6& a_9/a_6& a_{10}/a_6\end{array}\right]$$

1.4.非线性寻优法

利用两步法的结果作为初值，然后应用非线性寻优方法例如牛顿法、Levenberg-Marquardt法等进行迭代求解。

1.5.椭球拟合法

此方法是直接拟合公式（6）中所表示的椭球，可以参考各种球面拟合方法。matlab中有相关函数可供调用。

1.6.信号预处理—形态滤波与HHT法

此方法主要是利用形态学滤波去除传感器信号中的脉冲噪声，用希尔伯特-黄变换HHT识别高频噪声，以消除随机噪声对信号的干扰，后续的标定过程还是参考上面的方法⁴。

2、陀螺仪/加速度计辅助法

可以参考武元新老师发表的这篇文章⁵，文章还公布了文章的数据，方便复现研究。该文章中还考虑到了磁强计与IMU间微小姿态转换问题。
不过当时复现时候，文章中（20）中Hm调整第一项$SCm\times$的正负号才收敛。自己推一下就看出来了，应该是笔误。

三、总结

参考文献

---

1. 吴志添， 武元新， 胡小平， 吴美平. 基于总体最小二乘的捷联三轴磁力仪标定与地磁场测量误差补偿[J]. 兵工学报, 2012, 33(10): 1202-1209. WU Zhi-tian, WU Yuan-xin, HU Xiao-ping, WU Mei-ping. Calibration of Strapdown Three-axis Magnetometer and Measurement Error Compensation of GeomagneticField Based on Total Least Squares. Acta Armamentarii, 2012, 33(10): 1202-1209. ↩︎ ↩︎ ↩︎ ↩︎ ↩︎ ↩︎

2. 高可,林新华,储志伟,武靖博.无磁转台的电子罗盘误差分离标定方法[J].传感器与微系统,2017,36(02):21-24. ↩︎

3. 李保国, 陈克川, 李淑影. 一种基于递推最小二乘法的磁强计现场快速标定方法. ↩︎ ↩︎

4. 刁云云,高金耀,吴国超,等.基于形态学-HHT算法的船载地磁三分量信号分析与预处理[J].海洋学研究, 2021, 39(3):9.DOI:10.3969/j.issn.1001-909X.2021.03.005. ↩︎

5. Y. Wu, D. Zou, P. Liu and W. Yu, “Dynamic Magnetometer Calibration and Alignment to Inertial Sensors by Kalman Filtering,” in IEEE Transactions on Control Systems Technology, vol. 26, no. 2, pp. 716-723, March 2018, doi: 10.1109/TCST.2017.2670527. ↩︎