一、地磁传感器介绍

1、地磁传感器概念：

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

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

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

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

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

暂不考虑ωb\omega ^bωb，对上式移项取逆后可得地磁传感器输出校正模型： B^ee=CbeAM−1(B^b−bb)(2)\hat{B}^{e}_e=C^e_bA_M^{-1}(\hat{B}^{b}-b^b)\tag{2} B^ee=CbeAM−1(B^b−bb)(2) 其中B^ee\hat{B}^{e}_eB^ee为地磁场矢量的估计值，CbeC^e_bCbe为从传感器坐标系到地理坐标系的方向余弦矩阵，且有Cbe=(Ceb)−1C^e_b=(C^b_e)^{-1}Cbe=(Ceb)−1。

二、地磁传感器标定方法

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

1.自标定

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

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

对于式（4），有 R02=∣∣RA(B^ib−bb)∣∣2,i=1,2,...,m(5)R^2_0= ||R_A( \hat{B}^{b} _i-b^b)||^2,i=1,2,...,m\tag{5} R02=∣∣RA(B^ib−bb)∣∣2,i=1,2,...,m(5) 其中，R02R^2_0R02为当地地磁场矢量的模值；B^ib\hat{B}^{b}_iB^ib为传感器输出矢量；iii为测量序列次数；m为采样个数。将上式展开，可得， a1xi2+a2xiyi+a3xizi+a4yi2+a5yizi+a6zi2+a7xi+a8yi+a9zi+a10=0,i=1,2,...,m,(6){a_1}{x^2_i}+{a_2}{x_i}{y_i}+{a_3}{x_i}{z_i}+{a_4}{y^2_i}+{a_5}{y_i}{z_i}+ {a_6}{z^2_i}+{a_7}{x_i}+{a_8}{y_i}+{a_9}{z_i}+{a_{10}}=0, i=1,2,...,m,\tag{6} a1xi2+a2xiyi+a3xizi+a4yi2+a5yizi+a6zi2+a7xi+a8yi+a9zi+a10=0,i=1,2,...,m,(6) 式中：xi,yi,zix_i,y_i,z_ixi,yi,zi分别为B^ib\hat{B}^{b}_iB^ib的3个分量，a1−a10a_1-a_10a1−a10是由参数RAR_ARA和bbb^bbb构成的联合变量。将(6)式写成矩阵形式，有 HX=W,(7)HX=W,\tag{7} HX=W,(7) 式中： X=[x12x1y1x1z1y12y1z1x1y1z11x22x2y2x2z2y22y2z2x2y2z21......xm2xmymxmzmym2ymzmxmymzm1]X=\begin{bmatrix} x_1^2 & x_1y_1 & x_1z_1 &y_1^2 &y_1z_1 &x_1 &y_1 &z_1 &1\\ x_2^2 & x_2y_2 & x_2z_2 &y_2^2 &y_2z_2 &x_2 &y_2 &z_2 &1\\ ... &...\\ x_m^2 & x_my_m & x_mz_m &y_m^2 &y_mz_m &x_m &y_m &z_m &1\\ \end{bmatrix} X=⎣⎢⎢⎡x12x22...xm2x1y1x2y2...xmymx1z1x2z2xmzmy12y22ym2y1z1y2z2ymzmx1x2xmy1y2ymz1z2zm111⎦⎥⎥⎤ X=[a1/a6a2/a6a3/a6a4/a6a5/a6a7/a6a8/a6a9/a6a10/a6]TX=\begin{bmatrix} 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{bmatrix}^T X=[a1/a6a2/a6a3/a6a4/a6a5/a6a7/a6a8/a6a9/a6a10/a6]T W=[−z12−z22...−zm2]TW=\begin{bmatrix} -z^2_1 &-z^2_2 &...&-z^2_m \end{bmatrix}^T W=[−z12−z22...−zm2]T H∈ℜm×n,X∈ℜn×1,W∈ℜm×1H\in \Re^{m\times n},X\in \Re^{n\times 1},W\in \Re^{m\times 1}H∈ℜm×n,X∈ℜn×1,W∈ℜm×1;一般来说，m>nm>nm>n,因此(7)式为超定方程¹。以上理论推导参考文献¹。

下面介绍在求得X之后，如何得到RA,bbR_A,b^bRA,bb。
由向量2范数的定义∣∣X∣∣2=XTX||X||^2=X^TX∣∣X∣∣2=XTX，(5)式可展开为：
R02=[RA(B^ib−bb)]T[RA(B^ib−bb)]=(B^ibT−bbT)RATRA(B^ib−bb),i=1,2,...,m.(9)R^2_0= [R_A( \hat{B}^{b} _i-b^b)]^T[R_A( \hat{B}^{b} _i-b^b)]=( {\hat{B}^{b}_i}^T -{b^b}^T){R_A}^TR_A( \hat{B}^{b} _i-b^b),i=1,2,...,m.\tag{9} R02=[RA(B^ib−bb)]T[RA(B^ib−bb)]=(B^ibT−bbT)RATRA(B^ib−bb),i=1,2,...,m.(9)
将上式继续展开，可得
R02=B^ibTRATRAB^ib−B^ibTRATRAbb−bbTRATRAB^ib+bbTRATRAbb,i=1,2,...,m.(10)R^2_0= {\hat{B}^{b}_i}^T{R_A}^TR_A {\hat{B}^{b}_i}-{\hat{B}^{b}_i}^T{R_A}^TR_Ab^b-{b^b}^T{R_A}^TR_A\hat{B}^{b} _i+ {b^b}^T{R_A}^TR_Ab^b,i=1,2,...,m.\tag{10} R02=B^ibTRATRAB^ib−B^ibTRATRAbb−bbTRATRAB^ib+bbTRATRAbb,i=1,2,...,m.(10)
由于B^ib∈ℜ3×1,RA∈ℜ3×3,bb∈ℜ3×1\hat{B}^{b}_i \in \Re^{3\times 1},R_A \in \Re^{3\times 3},b^b \in \Re^{3\times 1}B^ib∈ℜ3×1,RA∈ℜ3×3,bb∈ℜ3×1上式中B^ibTRATRAbb{\hat{B}^{b}_i}^T{R_A}^TR_Ab^bB^ibTRATRAbb为标量，对其取转置值不变即
B^ibTRATRAbb=(B^ibTRATRAbb)T=bbTRATRAB^ib，i=1,2,...,m.(11){\hat{B}^{b}_i}^T{R_A}^TR_Ab^b=({{\hat{B}^{b}_i}^T{R_A}^TR_Ab^b})^T={b^b}^T{R_A}^TR_A\hat{B}^{b} _i，i=1,2,...,m.\tag{11} B^ibTRATRAbb=(B^ibTRATRAbb)T=bbTRATRAB^ib，i=1,2,...,m.(11)
于是，（10）式可简化为：
R02=B^ibTRATRAB^ib−2bbTRATRAB^ib+bbTRATRAbb,i=1,2,...,m.(12)R^2_0= {\hat{B}^{b}_i}^T{R_A}^TR_A {\hat{B}^{b}_i}-2{b^b}^T{R_A}^TR_A\hat{B}^{b} _i+ {b^b}^T{R_A}^TR_Ab^b,i=1,2,...,m.\tag{12} R02=B^ibTRATRAB^ib−2bbTRATRAB^ib+bbTRATRAbb,i=1,2,...,m.(12)
其中，第一项为二次型，第二项为一次项，最后一项为常数项。
显然RATRA{R_A}^TR_ARATRA为对称阵，令RATRA=M=[m1m2m3m2m4m5m3m51]，B^ib=[xyz]T{R_A}^TR_A=M=\begin{bmatrix} m_1 &m_2 &m_3\\ m_2 &m_4 &m_5\\ m_3 &m_5 &1\\ \end{bmatrix}，{\hat{B}^{b}_i}={\begin{bmatrix} x &y &z \end{bmatrix}}^TRATRA=M=⎣⎡m1m2m3m2m4m5m3m51⎦⎤，B^ib=[xyz]T,则
B^ibTRATRA=B^ibTMB^ib=[xyz][m1m2m3m2m4m5m3m51][xyz]=m1x2+m4y2+z2+2m2xy+2m3xz+2m5yz.(13){\hat{B}^{b}_i}^T{R_A}^TR_A= {\hat{B}^{b}_i}^TM{\hat{B}^{b}_i}=\begin{bmatrix} x &y &z \end{bmatrix}\begin{bmatrix} m_1 &m_2 &m_3\\ m_2 &m_4 &m_5\\ m_3 &m_5 &1\\ \end{bmatrix}\begin{bmatrix} x\\y\\z \end{bmatrix}\\=m_1x^2+m_4y^2+z^2+2m_2xy+2m_3xz+2m_5yz.\tag{13} B^ibTRATRA=B^ibTMB^ib=[xyz]⎣⎡m1m2m3m2m4m5m3m51⎦⎤⎣⎡xyz⎦⎤=m1x2+m4y2+z2+2m2xy+2m3xz+2m5yz.(13)
对比上式与（7）式，可得
M=[m1m2m3m2m4m5m3m51]=[X(1)X(2)/2X(3)/2X(2)/2X(4)/2X(5)/2X(3)/2X(5)/21](14)M=\begin{bmatrix} m_1 &m_2 &m_3\\ m_2 &m_4 &m_5\\ m_3 &m_5 &1\\ \end{bmatrix}= \begin{bmatrix} 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{bmatrix}\tag{14} M=⎣⎡m1m2m3m2m4m5m3m51⎦⎤=⎣⎡X(1)X(2)/2X(3)/2X(2)/2X(4)/2X(5)/2X(3)/2X(5)/21⎦⎤(14)
对MMM取乔伊斯基变换可得RA=[cholesky(M)]TR_A=[cholesky(M)]^TRA=[cholesky(M)]T
另外，对于一次项，可采用类似比较系数的方法，得
−2bbTRATRA=[X(6)X(7)X(8)](15)-2{b^b}^T{R_A}^TR_A=\begin{bmatrix} X(6) &X(7) &X(8) \end{bmatrix}\tag{15} −2bbTRATRA=[X(6)X(7)X(8)](15)
又RATRA=M{R_A}^TR_A=MRATRA=M，所以有
bbT=−12[X(6)X(7)X(8)]M−1(16){b^b}^T=-\frac{1}{2}\begin{bmatrix} X(6) &X(7) &X(8) \end{bmatrix}M^{-1}\tag{16} bbT=−21[X(6)X(7)X(8)]M−1(16)

1.1.两步法

第一步：对于（7）式，假定矩阵H是精确的，W含有噪声，根据Gauss-markov定理，X的最小二乘估计（LS）为 XLS=(HTH)−1HTW,(17)X_{LS}=(H^TH)^{-1}H^TW,\tag{17} XLS=(HTH)−1HTW,(17) 第二步：根据算法X中间量的定义，求出RAR_ARA和bbb^bbb。

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

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

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

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

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

递推最小二乘法主要将批量处理地磁传感器数据改为递推形式，更适合在线执行³。
文献³所述的参数方程与上面推导的（6）式略有不同，但也只是系数归一化的所参照的参数项不同而已。
对（6）式稍做变换，写成矩阵形式，有
h(k)=[xi2xiyixiziyi2yizixiyizi1]h(k)=\begin{bmatrix} x_i^2 & x_iy_i & x_iz_i &y_i^2 &y_iz_i &x_i &y_i &z_i &1 \end{bmatrix} h(k)=[xi2xiyixiziyi2yizixiyizi1]
z(k)=[−zi2]\\z(k)=\begin{bmatrix} -z_i^2 \end{bmatrix}\\ z(k)=[−zi2]
θ(k)=[a1/a6a2/a6a3/a6a4/a6a5/a6a7/a6a8/a6a9/a6a10/a6]\theta(k)=\begin{bmatrix} 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{bmatrix} θ(k)=[a1/a6a2/a6a3/a6a4/a6a5/a6a7/a6a8/a6a9/a6a10/a6]

篇幅有限，以下内容开篇分述。

1.4.非线性寻优法

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

1.5.椭球拟合法

1.6.形态滤波与HHT法

2、陀螺仪/磁强计辅助法

三、总结

参考文献

吴志添，武元新，胡小平，吴美平. 基于总体最小二乘的捷联三轴磁力仪标定与地磁场测量误差补偿[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. ↩︎ ↩︎ ↩︎ ↩︎ ↩︎ ↩︎
高可,林新华,储志伟,武靖博.无磁转台的电子罗盘误差分离标定方法[J].传感器与微系统,2017,36(02):21-24. ↩︎
李保国, 陈克川, 李淑影. 一种基于递推最小二乘法的磁强计现场快速标定方法:. ↩︎ ↩︎

地磁传感器标定技术综述（1）相关推荐

地磁传感器还可以用于室内定位？
我们都知道无线地磁传感器普遍应用于路边停车系统以及路外停车系统的车位检测上,其具有检测精度高,不易受环境因素影响的特点,获得可用户的好评,是目前停车位检测系统应用比较广泛的方案.地磁传感器不经可以应用 ...
【图像融合】多传感器图像融合技术综述
多传感器图像融合多传感器图像融合属于多传感器信息融合的范畴,是指将不同传感器获得的同一景物的图像或同一传感器在不同时刻获得的同一景物的图像,经过去噪.时间配准.空间配准和重采样后,再运用某种融合 ...
自动驾驶系统的传感器标定方法
作者 | 黄浴编辑 | 汽车人原文链接:https://zhuanlan.zhihu.com/p/57028341 点击下方卡片,关注"自动驾驶之心"公众号 ADAS巨卷干货 ...
ICRA2022论文汇总 | 跟踪/检测/分割/传感器标定/SLAM/轨迹预测/导航规划等
点击下方卡片,关注"自动驾驶之心"公众号 ADAS巨卷干货,即可获取点击进入→自动驾驶之心技术交流群后台回复[ECCV2022]获取ECCV2022所有自动驾驶方向论文! Ro ...
【AI有识境】如何掌握好图像分割算法？值得你看的技术综述
大家好,这是专栏<AI有识境>的第二篇文章,讲述如何掌握好图像分割算法. 进入到有识境界,可以大胆地说自己是一个非常合格的深度学习算法工程师了,能够敏锐地把握自己研究的领域,跟踪前沿和能落 ...
【技术综述】你真的了解图像分类吗？
言有三毕业于中国科学院,计算机视觉方向从业者,有三工作室等创始人作者 | 言有三编辑 | 言有三图像分类是计算机视觉中最基础的任务,基本上深度学习模型的发展史就是图像分类任务提升的发展历史,不 ...
matlab 光平面标定代码,线结构光视觉传感器/线激光深度传感器标定工具（matlab）...
线结构光视觉传感器/线激光深度传感器标定工具(matlab) 发布时间:2018-03-13 17:32, 浏览次数:655 , 标签: matlab 线结构光视觉系统有着结构简单.使用灵活.受周围光 ...
领域最全！多传感器融合方法综述！(Camera/Lidar/Radar等多源异构数据)
点击下方卡片,关注"自动驾驶之心"公众号 ADAS巨卷干货,即可获取点击进入→自动驾驶之心技术交流群后台回复[ECCV2022]获取ECCV2022所有自动驾驶方向论文! 自动 ...
360度全景标定方法_全景摄像机标定方法综述
全景摄像机标定方法综述张静 1 , 柴兴华 2 , 李小英 1 , 裴春琴 1 [摘要] 简要阐述了折反射系统的几何构成和全景摄像机系统的两种常用模型, 详细分析了近年来全景摄像机的标定方法 . ...

地磁传感器标定技术综述（1）

一、地磁传感器介绍

1、地磁传感器概念：

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

二、地磁传感器标定方法