CN108919156B - 基于噪声补偿的三轴磁强计线下校正方法 - Google Patents

Abstract

本发明属于三轴磁强计校正领域，针对大噪声未知环境下三轴磁强计校正，提出一种基于无偏估计噪声幅值的测量噪声补偿方法。首先，在对固有误差和干扰磁场充分分析的基础上，合理构建载体坐标系下三轴磁强计误差校正模型，以使磁强计能够适应载体周围复杂的电磁环境；其次，针对经典代数法求解模型参数时忽略噪声项所带来的缺陷，建立无偏代价函数求解测量数据矩阵的噪声补偿矩阵，并通过未知环境下噪声方差幅值的合理求解，来减小测量噪声对参数估计的影响；最后，利用最小二乘法对校正模型的各个参数做最优估计，并加入椭球约束，使得到的结果更符合实际物理意义。

Description

基于噪声补偿的三轴磁强计线下校正方法

技术领域

本发明涉及三轴磁强计校正领域，用于解决大噪声未知环境下机载磁强计的线下校正问题。

背景技术

三轴磁强计是测量磁场强度的传感器，可以提供精确和可靠的近地磁场向量信息，并且不会随时间漂移，可和陀螺仪等其他传感器一起提供位置信息，因此广泛应用于导航系统等领域。然而，由于受制造工艺的限制，三轴磁强计的输出有着固有误差(如非正交误差、尺度因素误差、零偏)，并且经常被各种干扰磁场(如软磁、硬磁、测量噪声)影响。因此，在应用之前进行校正具有重要意义。

磁强计校正经典方法称为摇摆算法。在已知航向角的情况下，水平仪绕着唯一垂直方向的轴旋转。这种方法虽然简单易操作，但是有需要外部独立航向测量信息的缺点。与此同时，该方法是航向校正方法，不适合除航向确定以外的其他校正。而椭球校正方法不需要外部独立航向信息，并且能提供多种误差源的校正。椭球校正都是基于一个事实，即在补偿地点，无论传感器是什么方向，地磁场幅值应相同。在这种情况下，通过旋转设备的三个轴得到的测量值应该在一个球面上，球面半径为地磁场强度。

然而由于固有误差和干扰磁场的存在，测量值生成的是椭球而不是标准的球。因此，在建立误差校正模型之后，三轴磁强计的校正问题可转化为三维椭球参数辨识问题，待得到椭球参数之后再进行误差补偿。经典的普通最小二乘法以代数拟合的方式极小化测量值与椭球体之间的误差，计算简便，不需要迭代，是应用最广泛的方法，目前大多数椭球参数辨识方法都是该方法的延伸。但该方法是一种有偏估计，鉴于噪声的存在将无法收敛到真值。目前大多数研究在椭球参数辨识的过程中均直接或者间接地忽略掉了测量噪声的影响。而事实上，这样做在某种程度上是不合理的。由于安装磁强计的载体不同，噪声水平有很大差别，如果精度要求低或噪声水平低，误差可忽略不计，但在某些高精度或者大噪声应用背景下忽略噪声引起的误差将导致结果不可用。有研究单独针对椭球拟合问题设计了一种自适应最小二乘算法，得到了无偏估计，但是需要已知噪声方差。实际情况中，由于载体的结构及环境不同，噪声方差难以准确预知。可见，如何构建大噪声未知环境下的机载磁强计线下校正方法是当前需要解决的关键问题。

针对现有线下校正方法在处理噪声方面的不足，本发明提出了带有噪声补偿的三轴磁强计线下校正方法。该方法在建立了较完整的误差校正模型基础上，提出三轴磁强计校正过程中的噪声补偿问题，在噪声方差幅值未知的情况下设计了测量数据噪声补偿方法；然后，采用最小二乘方法对校正模型参数进行最优估计，并加入了更符合实际物理意义的椭球约束，使得结果不是任意二次曲面。通过这种方法实现大噪声未知环境下机载三轴磁强计的线下高精度校正。

发明内容

本发明的目的在于突破传统最小二乘法的限制，提升三轴磁强计校正方法的精度，扩大应用范围，解决大噪声未知环境下机载机载磁强计的线下校正难题，提出一种带有噪声补偿的磁强计线下校正新方法。其中主要解决的问题包括：

1)现有的基于代数意义距离的最小二乘椭球校正方法是一种有偏算法，计算时忽略了噪声项。在噪声较大，或者噪声水平难以准确预估的情况下，校正精度明显下降，甚至有可能难以满足要求；

2)现有的最小二乘椭球校正方法是一种代数距离意义下的校正方法，如果没有额外添加椭球约束，求解得任意二次曲面，而不是椭球，不符合物理意义，将导致校正误差的发散。

本发明所述的基于噪声补偿的三轴磁强计线下校正方法，其特征在于包括以下技术措施：

步骤一：利用无偏代价函数建立卷积等式，通过卷积等式求解噪声矩阵，从而将测量数据矩阵中的噪声项分离出来；

步骤二：利用行列式为零的特性解得噪声矩阵中的噪声方差幅值，从而完整得到去掉噪声项后的地磁场测量数据；

步骤三：将完整去掉噪声项后的地磁场测量数据矩阵带入无偏代价函数，联合椭球约束建立约束最小二乘问题，利用广义特征值问题求解方法求解该约束最小二乘问题，得到椭球具体参数；

步骤四：利用奇异值分解的方法分解椭球参数矩阵以得到三轴磁强计测量误差参数，将该参数带入磁强计测量模型，补偿测量误差，从而校正磁强计的测量值。

对比现有技术，本技术方案所述的基于噪声补偿的三轴磁强计线下校正方法，有益效果在于：

(1)本发明提出的噪声补偿方法不需要已知噪声方差幅值，可以滤除测量噪声对磁强计校正精度的影响，使得校正方法能够适应大噪声未知环境；

(2)本发明在数据拟合的过程中考虑了椭球约束，使得拟合的结果不仅仅是代价距离的极小化，不会获得任意曲面，而是严格符合椭球方程，与实际物理意义更吻合，并且精度更高。

附图说明

附图1是基于噪声补偿的三轴磁强计线下校正方法实施流程图。

附图2是三轴磁强计测量误差对测量数据的影响原理图。

附图3是测量地磁场数据噪声补偿基本原理图。

具体实施方式

针对大噪声未知环境下三轴磁强计校正问题，提出一种基于无偏估计的测量噪声补偿方法。首先，在对固有误差和干扰磁场充分分析的基础上，合理构建载体坐标系下三轴磁强计误差校正模型，以使磁强计能够适应载体周围复杂的电磁环境；其次，针对经典代数法求解模型参数时忽略噪声项所带来的缺陷，建立无偏代价函数求解测量数据矩阵的噪声矩阵，并通过未知环境下噪声方差幅值的合理求解，来减小测量噪声对参数估计的影响；最后，利用最小二乘法对校正模型的各个参数做最优估计，并加入椭球约束，使得到的结果更符合实际物理意义。

以下结合说明书附图1对本发明做进一步的详细描述。参照说明书附图1，本发明的处理流程分以下步骤：

1)噪声矩阵求解

磁强计的测量误差包括：传感器零位偏差、尺度因素误差、非正交误差、测量噪声、硬磁干扰误差和软磁干扰误差，下面分别加以分析。

受制造工艺限制，各个传感器特性不一致，不同的TAM具有不同的零位偏差c_hi和不同的尺度因子。尺度因素误差矩阵D_sf可以表示为

其中，s_x、s_y、s_z为三个轴的尺度因子。

TAM的敏感轴不能完全正交，从而导致非正交误差，误差矩阵D_no为

其中，ρ为实际y_s轴与理想y'_s轴的夹角；φ为实际z_s轴与理想z'_s轴的夹角在x_sO_sz_s平面的投影；λ为实际z_s轴和理想z'_s轴的夹角在y_sO_sz_s平面的投影。

受热特性和电路测量因素的影响，磁强计输出不可避免地具有测量噪声，可建模为高斯白噪声ε，服从均值为0，方差为σ²的正态分布。

由于飞行器的结构件或电子设备中包含铁磁性材料，所以TAM的测量值会受到除地球磁场外的其他磁场干扰。如果干扰磁场为永久时不变磁场，则称为硬磁干扰。硬磁干扰在磁强计的输出中表现为固定零位偏移c₀。如果干扰磁场是随载体姿态变化的，则为软磁干扰。软磁干扰在TAM的输出中表现为变化的零位偏差，软磁干扰误差矩阵D_si为

式中，α_ij(i,j＝x,y,z)为由j轴外部磁场诱导产生的磁场在i轴上的投影。软磁干扰误差可等效为尺度因素误差、非正交误差和非对准误差的综合效应。

综合这几种误差，TAM校正数学模型可以表示为

B_m＝D_sfD_noD_si(B_b+c_hi)+c₀+ε

其中，B_m和B_b为载体坐标系下TAM的实际输出值与无误差源干扰的理论真实值，令D_E＝D_sfD_noD_si，c＝D_Ec_hi+c₀，则TAM误差校正数学模型可表示为

B_m＝D_EB_b+c+ε

忽略ε，设

B'_m＝D_EB_b+c

令x表示实际测量值B_m，

为理想测量值B'_m，则

其中，x＝[xyz]。磁强计静止并且只改变方向时，磁场向量的幅值B₀为常数。给定一系列磁强计理想测量点

此时可以发现

严格分布在一个三维椭球体上(校正误差对测量值的影响见附图2)：

其中，

为了求解椭球参数，用代数法表示理想测量点与待求椭球之间的距离，则代价函数可以写为

上式是目前常用的校正模型参数求解方法，普通方法在代价函数(1)中用实际测量值x直接代替理想测量值

认为噪声项ε较小，可直接将其忽略掉。但事实上，对于机载磁强计来说，由于周围环境的差异，ε有时是相当大的。给出式(1)中用x代替

时忽略的噪声项：

ε'＝ε^TD_mε+(x-c)^TD_mε+ε^TD_m(x-c)

从上式可以看出，虽然Eε＝0，但是Eε'＝tr{D_mσ²}>0。因此，1)当测量数据噪声σ²较大时，忽略ε'会引起较大估计偏差，求解结果不能收敛到真值。同时，2)ε'和待求的参数D_m，也即尺度因素误差、非正交误差和软磁干扰的大小有直接关系。当误差参数较大时，省略ε'项也会引起较大误差。综合以上分析，目前的方法并不能准确估计磁强计校正模型参数。

为了消除实际测量数据中的白噪声，定义了无偏代价函数，根据无偏代价函数求得了测量数据噪声补偿矩阵，参照附图3，具体过程如下：

在式(1)中，用x代替

然后将其展开得

其中，b＝-2D_mc，

令

则子代价函数q(·)描述了单个测量点与椭球面F的偏离程度。式(2)即可写为

为了得到式(2)结果，引入β表示待求参数向量。令

为将对称矩阵D_m的上三角部分中的各个元素放入向量里的操作。那么待求参数向量为

β:＝[vec_s(D_m)^T b^T d]^T

它是空间

的一个子空间，n_β:＝(n+1)n/2+n+1。

定义Kronecker乘积

为

这样子代价函数q(·)也可表示为

因此代价函数(2)可以表示为

其中，

式(2)描述的最小二乘法虽然简单，但直接最小化代数误差不是无偏的。当测量噪声较大时，这种偏差不可以忽略。因此，为得到无偏估计，本发明通过采用无偏代价函数消除数据矩阵Y中噪声的影响，无偏子代价函数q₀(·)可以通过下面的式子进行定义

则无偏代价函数Q₀(·)可以表示为

令

从而

可以看出，随着采样点数n的增大，Q₀将逐步逼近

如果

有全局最优解

那么Q₀就是绝对无偏的。下面详细描述无偏估计方法：

因此

q₀(β；x)＝β^Tψ₀(x)β

其中

假设数据噪声方差

μ表示噪声方差幅值，

σ＝diag(σ_x,σ_y,σ_z)。其中，

已知，噪声方差幅值μ未知。由于ε_l独立同分布，因此

这样就可以得到如下卷积公式

通过上式可以得到多项式ψ₀：

ψ₀＝ψ+Δψ₀

那么，消除白噪声影响后的数据矩阵Y₀即为

其中，矩阵C即为噪声补偿矩阵，

目前，虽然得到了噪声补偿矩阵C，但是矩阵C中的噪声方差幅值μ仍然是未知的。

2)噪声幅值求解

由于μ未知，此时C是关于μ的二次型，将C中的一次、二次项分开表示，Y₀可写为

Y₀(μ)＝Y-C(μ)＝Y-μC₁-μ²C₂

其中，C₁是由σ_x、σ_y、σ_z以及测值的期望确定的(例如：Ex²、Ey²、Ez²、Exy、Exz等)，C₂是由σ_x、σ_y、σ_z确定的。校正后的矩阵Y₀(μ)则用于求解带约束的二次曲面拟合问题。由于我们寻找的是符合方程y^Tβ＝0的线性拟合，因此Y₀是奇异矩阵，可以通过下式求解μ：

det(Y-μC₁-μC₂)＝0

事实上，β和Y₀的最后一个参数都是冗余的，因此可以将Y₀经行分块，得到

其中，Y_*表示非常数项(二次和线性项)，y_*表示常数和非常数的混合项，y₀＝1。按照同样的方法也可以将参数向量分开

其中，

用分解Y₀和β的方法分解Y和C可得到

由于常数项的噪声应为0，因此矩阵C右下角的元素为0。写出矩阵Y₀右下角元素的Schur补，然后可以得到一阶更新：

这是一个关于μ的二次多项式，接下来即需求解未知数μ：

det(Y_*(μ))＝det(M₀-μM₁-μ²M₂)＝0

线性化上式：

其中，W为任意满秩矩阵，令W＝-M₀即可得到一个广义特征值问题。至此，根据本文的方法就可以得到μ，也就可以完全得到噪声补偿后的数据矩阵Y₀。

3)椭球参数的求解

如果没有约束的存在式(2)将求解得任意二次曲面，而在实际情况中，可能只关心某种具体类型的二次函数。由于磁强计校正模型最终形式为椭球体，因此本文加入椭球约束，使得拟合结果为严格的椭球体。这样做虽然有可能牺牲了精度，但结果符合实际物理意义，更有价值。根据上述分析，磁强计校正问题就可以转化成下面带椭球约束的最小二乘问题

s.t.β^TTβ＝1

这就是一个广义特征值问题

Y₀β＝λTβ (3)

矩阵T代表约束。在约束矩阵T中只涉及x²、y²、z²、xy、yz、xz这几个非线性项，而和线性项x、y、z无关。也可以说T中的大部分元素为0。式(3)中的向量和矩阵可以分解为

其中，T₁代表可以确定二次曲面为椭球的约束矩阵。那么，式(3)可以写为

令

Ypol＝Y₁+Y₂Y'

其中，除测量点都在一个平面的情况，Y₃ ^-1总是存在。如果Y₃ ^-1是奇异的，则测量点

分布在同一个平面上，此时拟合椭球面也就没有意义。因此

通过上述分析，现可以通过求解Ypol满足约束

的特征向量β₁来估计磁强计误差校正参数，求解如下广义特征值问题：

Ypol·β₁＝λT₁β₁

对于椭球拟合的约束，可选择为

其中，k≥4，Ones为全1矩阵，I为单位阵，0为零阵。令

β＝[a₁ a₂ a₃ a₄ a₅ a₆ a₇ a₈ a₉ a₁₀]^T

β_I＝a₁+a₂+a₃

当

时，参数β代表椭球体。同时，椭球体短半径长度是长半径一半以上时，可以得到

即k＝4。

4)三轴磁强计误差补偿

在得到椭球体参数之后，就可以对磁场测量值进行误差补偿。对D_m进行奇异值分解可得

D_m＝UΣU^T

式中，U为与对称矩阵

特征值对应的特征矢量矩阵。

假设

那么

由于

所以可以利用上式求取D_E ^-1。从而TAM的补偿方程可以写为

Claims (1)

1.基于噪声补偿的三轴磁强计线下校正方法，其特征在于，包括以下步骤：

步骤一：利用无偏代价函数建立卷积等式，通过卷积等式求解噪声矩阵，从而将测量数据矩阵中的噪声项分离出来；

步骤二：利用行列式为零的特性解得噪声矩阵中的噪声方差幅值，从而完整得到去掉噪声项后的地磁场测量数据；

步骤三：将完整去掉噪声项后的地磁场测量数据矩阵带入无偏代价函数，联合椭球约束建立约束最小二乘问题，利用广义特征值问题求解方法求解该约束最小二乘问题，得到椭球具体参数；

步骤四：利用奇异值分解的方法分解椭球参数矩阵以得到三轴磁强计测量误差参数，将该参数带入磁强计测量模型，补偿测量误差，从而校正磁强计的测量值；

步骤一中测量数据矩阵中的噪声项分离方法为：

利用无偏代价函数，设计噪声矩阵求解方法，从而分离测量矩阵中的白噪声，

测量噪声模为高斯白噪声ε，服从均值为0，方差为σ²的正态分布，ε_l表示第l个测量值的噪声，l＝1，...,n，由于ε_l独立同分布，根据卷积公式得到测量数据矩阵多项式ψ₀，那么，消除白噪声影响后的数据矩阵Y₀即为

其中，矩阵C即为噪声矩阵，Y表示分离出噪声矩阵后的测量数据矩阵，x_l表示磁强计测量点，

μ表示噪声方差幅值，

表示单位方差，

σ＝diag(σ_x,σ_y,σ_z)，σ_x、σ_y、σ_z分别为x、y、z轴的噪声方差，其中，

已知，噪声方差幅值μ未知；

步骤二中噪声方差幅值求解方法为：

在分离出测量数据矩阵中的噪声项的基础上，再将测量数据矩阵中关于噪声幅值的一次项和二次项分开，利用行列式为零的特性，解得噪声幅值，从而确定噪声矩阵的具体值，

将C中的一次、二次项分开表示，Y₀写为

Y₀(μ)＝Y-C(μ)＝Y-μC₁-μ²C₂

其中，C₁是由σ_x、σ_y、σ_z以及测值的期望确定的，C₂是由σ_x、σ_y、σ_z确定的，Y₀是奇异矩阵，通过下式求解μ：

det(Y-μC₁-μC₂)＝0

将Y₀经行分块，得到

其中，Y_*表示非常数项，即二次和线性项，y_*表示常数和非常数的混合项，y₀＝1，令β表示待求参数向量，按照同样的方法将参数向量分开

其中，

用分解Y₀和β的方法分解Y和C得到

写出矩阵Y₀右下角元素的Schur补，然后得到一阶更新：

上式中的未知数μ求解方法为

det(Y_*(μ))＝det(M₀-μM₁-μ²M₂)＝0

线性化式(2)：

其中，W为任意满秩矩阵，令W＝-M₀即得到一个广义特征值问题，至此就得到μ，也就完全得到噪声补偿后的数据矩阵Y₀；

步骤三中带椭球约束的椭球参数求解方法为：

磁强计校正问题转化成下面带椭球约束的最小二乘问题

s.t.β^TTβ＝1

这就是一个广义特征值问题

Y₀β＝λTβ (4)

矩阵T代表约束，式(3)中的向量和矩阵分解为

其中，T₁代表确定二次曲面为椭球的约束矩阵，那么，式(4)写为

令

Ypol＝Y₁+Y₂Y'

其中，除测量点都在一个平面的情况，

总是存在，

现通过求解Ypol满足约束

的特征向量β₁来估计磁强计误差校正参数，求解如下广义特征值问题：

Ypol·β₁＝λT₁β₁

对于椭球拟合的约束，选择为

其中，k≥4，Ones为全1矩阵，I为单位阵，0为零阵，令

β＝[a₁ a₂ a₃ a₄ a₅ a₆ a₇ a₈ a₉ a₁₀]^T

β_I＝a₁+a₂+a₃

当

时，参数β代表椭球体，同时，椭球体短半径长度是长半径一半以上时，得到

步骤四中三轴磁强计误差补偿方法为：

在得到椭球体参数之后，就进一步设计磁场测量值误差补偿方法，对D_m进行奇异值分解得

D_m＝UΣU^T

式中，U为与对称矩阵

特征值对应的特征矢量矩阵，令

那么

由于

所以利用上式求取D_E ^-1，从而三轴磁强计的补偿方程写为

式中，B_m和B_b为载体坐标系下TAM的实际输出值与无误差源干扰的理论真实值。

Images