CN113281824B - 一种考虑飞机非刚性以及极化电流因素的航空磁补偿方法 - Google Patents

Abstract

本发明公开了一种考虑飞机非刚性以及极化电流因素的航空磁补偿方法。本发明的方法包括：建立经典的磁补偿T‑L模型，建立飞机平台干扰磁场模型，建立三分量磁通门磁力仪误差模型的磁补偿模型，得到极化电流产生的干扰磁场，对地磁场进行补偿，建立光泵磁力仪方向误差，建立最终的飞机平台磁干扰模型，使用总体最小二乘方法求解出磁补偿系数向量，将求解出的磁补偿系数向量补偿系数，用于飞机平台的实时补偿；本发明将原先的18项T‑L模型扩展到50项甚至更多，解决三了轴磁通门磁力仪误差问题和飞机在非刚性情况下和极化电流产生的干扰磁场的数学模型问题；本发明大大地提高了航磁补偿的精度，减小了飞机平台本身产生的磁场对探测目标信号磁场的干扰。

Description

一种考虑飞机非刚性以及极化电流因素的航空磁补偿方法

技术领域

本发明涉及航空磁补偿技术，具体涉及一种考虑飞机非刚性以及极化电流因素的航空磁补偿方法。

背景技术

航空磁补偿技术是为了补偿掉飞机平台产生的干扰磁场，有利于目标信号的检测。经典的T-L模型由三部分组成：恒定磁场、感应磁场和涡流磁场，通过分析每种磁场产生的原因来建立数学模型，即托尔斯·劳森(T-L)方程。1961年，Leach在文章《Aeromagneticcompensation as a linear regression problem》里表示补偿系数的大小与飞机本身的结构和磁力仪的安放位置有关系，所以如果使用同一架飞机并且磁力仪放置位置不变化，补偿系数也不会发生变化。同时Leach在论文《Automatic aeromagnetic compensation》中写出了系数的求解方案，通过飞机在规定航向上做俯仰、偏航、侧滚机动性工作来得到测量数据，由此求解出补偿系数。2016年，韩琦团队在《An aeromagnetic compensationcoefficient-estimating method robust to geomagnetic gradient》一文中建立了地磁梯度模型，补偿掉了环境磁场的一部分干扰。另外，他们团队发表的文章《A modifiedTolley-Lawson model robust to the errors of the three-axis strapdownmagnetometer》中指出三轴磁通门由于存在零偏性、非正交性、测量噪声等误差会导致数据存在误差，并对磁通门进行了修正。此外，郭弘老师团队发表了关于多光系光泵磁力仪的方向误差校正的模型的专利；在胡凯的硕士论文里提到了关于单光系光泵方向误差的处理方法。目前仍然使用广泛的是T-L模型，但是T-L模型建立的条件是飞机是刚性结构，在飞行过程中不会发生形变，且三轴磁通门磁力仪是正交的，没有零飘，飞机飞行区域的环境磁场是固定的，而且没有考虑到环境电场产生的极化电流因素。但是实际上这些条件在真实环境中是发生变化的，因此建立一个更精确的飞机平台磁干扰补偿模型是有重要意义的。

发明内容

本发明提出了一种考虑飞机非刚性以及极化电流因素的航空磁补偿方法，用以解决三轴磁通门磁力仪误差问题和飞机在非刚性情况下和极化电流产生的干扰磁场的数学问题。

本发明的考虑飞机非刚性以及极化电流因素的航空磁补偿方法，包括以下步骤：

1)建立经典的磁补偿T-L模型：

以飞机平台建立坐标系，建立经典的磁补偿T-L模型，得到飞机平台产生的干扰磁场B_TL(t)：

其中，a_i为飞机永久干扰磁场的补偿系数，b_i，j为飞机感应干扰磁场的补偿系数，c_i，j为飞机涡流干扰磁场的补偿系数，u_i(t)和u_j(t)为飞机平台坐标系与地磁场的方向余弦u(t)的分量，i＝1，2，3，j＝1，2，3，在飞机平台坐标系下分别为x、y和z轴，把它们展开写为u₁(t)＝cos X，u₂(t)＝cos Y，u₃(t)＝cos Z，其中X、Y和Z分别为飞机平台坐标系x、y和z轴与地磁场的夹角，cos X、cos Y和cos Z表达式分别如下：

其中，B₁(t)、B₂(t)和B₃(t)分别为在t时刻三分量磁通门磁力仪测量的沿飞机平台坐标系x、y和z轴的地磁场测量值；

2)建立飞机平台干扰磁场模型：

经典的T-L模型是在飞机平台为刚体这个基本假设条件建立的，但是飞机在转弯和做机动动作的时候，方向舵、升降舵和副翼这些活动部件发生了状态变化，机翼发生了弯曲，因此飞机在做机动动作的时候不能看成是一个理想刚体，当飞机在做侧滚、俯仰和偏航机动动作时，飞机的结构状况采用运动状态参数来描述，飞机平台坐标系下地磁场的方向余弦u(t)能够描述整体的飞机姿态，但是不能描述具体的各个活动部件的状态；

a)以飞机平台坐标系下地磁场的方向余弦的微分

反映活动部件的运动状态，把式(1)中的补偿系数写成u(t)和

的泛函：

b)将补偿系数以近似的线性化表示为：

c)建立飞机平台干扰磁场模型，得到修改补偿系数之后飞机平台产生的干扰磁场B_coe(t)：

其中，d_i，j，k为飞机主体感应磁场与活动部件永久磁场的耦合的补偿系数，g_i，j，k为感应磁场与活动部件涡流磁场的耦合的补偿系数；

3)建立三分量磁通门磁力仪误差模型的磁补偿模型：

三分量磁通门磁力仪存在非正交性和零漂问题，地磁场的真实值为：

其中，

为磁通门磁力仪测量的沿飞机平台坐标系x、y和z轴的地磁场真实值，

为修正矩阵，该矩阵是未知的，用来修正磁通门存在的非正交性问题，

为误差向量，也是未知的向量，用来修正三分量磁通门磁力仪存在的零漂问题；

将式(4)写成简洁形式为：

因此在飞机平台坐标系下，修正后的飞机平台坐标系下地磁场的方向余弦

表示为：

其中，σ(t)为时间缓变误差函数，

将

代入式(1)，且

整理后得到包括修正磁通门数据建立的磁补偿模型，得到修正磁通门数据后飞机平台产生的干扰磁场 B_cal(t)：

其中，a′_i为修正后的永久磁场的补偿系数，b′_i，j为修正后的感应磁场的补偿系数，c′_i，j为修正后的涡流磁场的补偿系数，h_i为由磁通门零漂和非正交性导致的磁干扰的补偿系数，β(t)是随时间变化的一次函数；

由于飞机的运动姿态也是随时间变化的，因此对式(5)中的补偿系数a′_i(t)、b′_i，j(t) 和c′_i，j(t)进行与式(2)形式一致的线性展开，再与式(5)合并之后得到通过校准磁通门数据和改变补偿系数形式而建立的飞机干扰磁场模型，得到通过校准磁通门数据和改变补偿系数形式之后飞机平台产生的干扰磁场B_interference(t)：

其中，h_i为由磁通门零漂和非正交性导致的磁干扰的补偿系数；

4)得到极化电流产生的干扰磁场：

在环境场中，不仅有磁场还有电场，经典的磁补偿T-L模型里并未考虑到环境电场产生的极化电流因素，因此忽略了极化电流产生的一部分干扰磁场；

海洋上空的环境电场E_env(r，t)为垂直指向海面的，则

E₀为环境电场的大小，e(t)为环境电场方向；

为由整体的飞机姿态组合成的矩阵，整体的飞机姿态以飞机平台坐标系下地磁场的方向余弦u(t)表示，因此把

写成关于u(t)的函数，于是把e(t)用u(t)来表示： e(t)＝e(u(t))；

环境电场E_env(r，t)使飞机的金属机壳产生正比于

的极化电流，这个极化电流产生的干扰磁场B_ele(t)表示为：

其中，

为环境电场方向e(t)的微分，

为环境电场方向的微分

的第j分量(j＝1，2，3)，

表达为：

其中，α_j为常数项系数，β_j，k为一次项系数，γ_j，k，l为二次项系数，M_i，j为极化电流产生的干扰磁场系数，将该表达式与极化电流产生的干扰磁场结合，得到极化电流产生的干扰磁场B_ele(t)：

其中，σ_i，k为极化电流磁场干扰的补偿系数，τ_i，k，l为飞机感应磁场与极化电流磁场的耦合的干扰磁场的补偿系数；

5)对地磁场进行补偿：

采用国际地磁参考场(IGRF)，用IGRF模型算出地磁场值后，再用光泵磁力仪测量的磁场值减去地磁场值，得到补偿掉地磁场后的剩余干扰磁场值；

6)建立光泵磁力仪的方向误差：

在实际实验中，会用到光泵磁力仪测量的总磁场，光泵磁力仪测量的总磁场都存在方向误差磁场，因此分别建立单光系或多光系光泵磁力仪的方向误差B_h(t)；

7)建立最终的飞机平台磁干扰模型：

将光泵磁力仪的方向误差加入到通过校准磁通门数据、改变补偿系数形式和增加极化电流磁场而建立的飞机干扰磁场模型中，建立最终的飞机平台磁干扰模型，得到最终的飞机干扰磁场B_f(t)：

B_f(t)＝B_interference(t)+B_ele(t)+B_h(t)

上式简化为：

B_f(t)＝A(t)*N (7)

其中，A(t)为飞机平台坐标系下地磁场的方向余弦组成的基函数矩阵，是u_i(t)和

的组合，N为待求的磁补偿系数向量；

8)使用总体最小二乘方法求解出磁补偿系数向量；

9)将求解出的磁补偿系数向量，用于飞机平台的实时补偿。

其中，在步骤6)中，建立单光系光泵磁力仪的方向误差，包括以下步骤：

对于单光系光泵磁力仪，主要用到的是铯光泵磁力仪，铯光泵磁力仪在转动时，其内部光轴与地磁场之间的夹角发生了改变，因此把方向误差写成关于飞机平台坐标系下地磁场的方向余弦u(t)的函数，经过泰勒函数展开得到的方向误差B_h(t)表达式为：

B_h(t)＝k₁*u₁(t)+k₂*u₂(t)+k₃*u₃(t)+k₄*u₁ ²(t)+k₅*u₂ ²(t)+k₆*u₃ ²(t)+k₇ *u₁(t)*u₂(t)+k₈*u₁(t)*u₃(t)+k₉*u₂(t)*u₃(t)

其中，k₁，k₂，k₃，k₄，k₅，k₆，k₇，k₈，k₉为9个待求解的方向误差系数。

在步骤6)中，建立多光系光泵磁力仪的方向误差，包括以下步骤：

多光系光泵磁力仪的方向误差B_h表达式为：

其中，z_i为方向误差系数，n为光泵原子气室个数，n＝3，分别沿着x，y，z轴方向；每个气室上的交变磁场线圈方向与光传播方向垂直，分别沿着y轴、x轴和x轴方向，h_i为与方向误差有关的方向余弦，表达式为：

其中，η_i为原子气室安装方向与地磁场的夹角，η_i＝X或Y或Z，ε_i为原子气室上加的交变磁场方向与地磁场的夹角，ε_i＝X或Y或Z。

在步骤8)中，使用总体最小二乘方法求解出磁补偿系数向量，包括以下步骤：

a)A和b分别为不可观测的无误差基函数数据矩阵和无误差飞机干扰磁场数据向量，实际观测的基函数数据矩阵A₀和飞机干扰磁场数据向量b₀分别为：

A₀＝A+ΔA，b₀＝b+Δb

其中，ΔA和Δb分别为误差数据矩阵和误差数据向量；根据式(7)得到飞机干扰数据向量和基函数数据矩阵的关系为：A*N＝b，N为待求的磁补偿系数向量，则若存在误差，该式的形式为：

即(A+ΔA)N＝b+Δb，等价为(B+D)q＝0；

式中，增广数掘矩阵B＝[A，b]和增广校正矩阵D＝[ΔA，Δb]均为m×(k+1)维矩

阵，q为总体最小二乘的解向量，而

为(k+1)×1向量，k为无误差基函数数掘矩阵A的列数，m为无误差基函数数掘矩阵A的行数；

b)计算增广数掘矩阵B的奇异值分解并存储右奇异矩阵V；

c)确定增广数掘矩阵B的有效秩p；

d)总体最小二乘的解向量

中，N为由p个线性独立的未知参数组成的列向量，因此总体最小二乘的求解就变成了k+1-p个总体最小二乘问题的求解： B_j ^(p)q＝0，j＝1，2...，k+1-p，B_j ^(p)为B的第j列到第p+j列组成的子矩阵；

e)根掘最小二乘原理，总体最小二乘解等价于求得解向量q使代价函数极小化，代价函数f(q)定义为：

其中

其中，u_j为增广数掘矩阵B左奇异矩阵第i行第j列上的元素，

为右奇异矩阵V第i行第j列上的元素；δ_j为增广数掘矩阵B的奇异值；

f)定义

则代价函数f(q)写成f(q)＝q^TS^(p)q的形式，经过整理得到：

g)代价函数f(q)的极小变化量是f(q)对q求导并令其导数等于0，求得S^(p)q＝λe₁，λ为常数表示误能量，e₁为单位向量，e₁＝[1，0，0...，0]^T；

h)求S^(p)的逆矩阵S^-(p)；

i)对上式S^(p)q＝λe₁求解，求得总体最小二乘的解向量q，因为

则磁补偿系数向量N＝[x_tls(1)，x_tls(2)，...，x_tls(p)]^T最优最小二乘近似解为：

其中，S^-(p)为S^(p)的逆矩阵，x_tls(i)为磁补偿系数向量N的第i个元素。

进一步，在飞机进行校准飞行前，在飞机上安装光泵磁力仪、三分量磁通门磁力仪、高度计和GPS；光泵磁力仪用来测量的总磁场，三分量磁通门磁力仪实时测量沿飞机平台坐标系x、y和z轴的地磁场测量值，高度计用来测量飞机飞行的高度，GPS用来测量飞机所处位置的经度和纬度。

本发明的优点：

本发明将原先的18项T-L模型扩展到50项甚至更多，解决三了轴磁通门磁力仪误差问题和飞机在非刚性情况下和极化电流产生的干扰磁场的数学模型问题；本发明大大地提高了航磁补偿的精度，减小了飞机平台本身产生的磁场对探测目标信号磁场的干扰。

附图说明

图1为本发明的考虑飞机非刚性以及极化电流因素的航空磁补偿方法的一个实施例的飞机飞行示意图；

图2为本发明的考虑飞机非刚性以及极化电流因素的航空磁补偿方法的流程图。

具体实施方式

下面结合附图，通过具体实施例，进一步阐述本发明。

在飞机进行校准飞行前，在飞机上安装光泵磁力仪、三分量磁通门磁力仪、高度计和GPS；光泵磁力仪用来测量的总磁场，三分量磁通门磁力仪实时测量沿飞机平台坐标系x、y和z轴的地磁场测量值，高度计用来测量飞机飞行的高度，GPS用来测量飞机所处位置的经度和纬度。

本实施例的考虑飞机非刚性以及极化电流因素的航空磁补偿方法，包括以下步骤：

1)建立经典的磁补偿T-L模型：

飞机在飞行过程中进行三种机动动作，俯仰、侧滚、偏航，每个机动动作四组，分为东西南北四个航向，如图1所示；

以飞机平台建立坐标系，建立经典的磁补偿T-L模型，得到飞机平台产生的干扰磁场B_TL(t)：

其中，a_i为飞机永久干扰磁场的补偿系数，b_i，j为飞机感应干扰磁场的补偿系数，c_i，j为飞机涡流干扰磁场的补偿系数，u_i(t)和u_j(t)为飞机平台坐标系与地磁场的方向余弦u(t)的分量，i＝1，2，3，j＝1，2，3，在飞机平台坐标系下分别为x、y和z轴，把它们展开写为u₁(t)＝cos X，u₂(t)＝cos Y，u₃(t)＝cos Z，其中X、Y和Z分别为飞机平台坐标系x、y和z轴与地磁场的夹角，cos X、cos Y和cos Z表达式分别如下：

其中，B₁(t)、B₂(t)和B₃(t)分别为在t时刻三分量磁通门磁力仪测量的沿飞机平台坐标系x、y和z轴的地磁场测量值；

2)建立飞机平台干扰磁场模型：

经典的T-L模型是在飞机平台为刚体这个基本假设条件建立的，但是飞机在转弯和做机动动作的时候，方向舵、升降舵和副翼这些活动部件发生了状态变化，机翼发生了弯曲，因此飞机在做机动动作的时候不能看成是一个理想刚体，当飞机在做侧滚、俯仰、偏航机动动作时，飞机的结构状况采用运动状态参数来描述，飞机平台坐标系下地磁场的方向余弦u(t)能够描述整体的飞机姿态，但是不能描述具体的各个活动部件的状态；

d)以方向余弦的微分

反映活动部件的运动状态，把式(1)中的补偿系数写成u(t)和

的泛函：

e)将补偿系数以近似的线性化表示为：

f)建立飞机平台干扰磁场模型，得到修改补偿系数之后飞机平台产生的干扰磁场B_coe(t)：

d_i，j，k为飞机主体感应磁场与活动部件永久磁场的耦合的补偿系数，g_i，j，k为感应磁场与活动部件涡流磁场的耦合的补偿系数；

3)建立三分量磁通门磁力仪误差模型的磁补偿模型：

三分量磁通门磁力仪存在非正交性和零漂问题，地磁场的真实值为：

其中，

为磁通门磁力仪测量的沿飞机平台坐标系x、y和z轴的地磁场真实值，

为修正矩阵，该矩阵是未知的，用来修正磁通门存在的非正交性问

题，

为误差向量，也是未知的向量，用来修正三分量磁通门磁力仪存在的零漂问题；

将式(4)写成简洁形式为：

因此在飞机平台坐标系下，修正后的地磁场的方向余弦

表示为：

其中，σ(t)为时间缓变误差函数，

将

代入式(1)，且

整理后得到包括修正磁通门数据建立的磁补偿模型，得到修正磁通门数据后飞机平台产生的干扰磁场B_cal(t)：

其中，a′_i为修正后的永久磁场的补偿系数，b′_i，j为修正后的感应磁场的补偿系数，c′_i，j为修正后的涡流磁场的补偿系数，h_i为由磁通门零漂和非正交性导致的磁干扰的补偿系数，β(t)是随时间变化的一次函数；

由于飞机的运动姿态也是随时间变化的，因此对式(5)中的补偿系数a′_i(t)、b′_i，j(t) 和c′_i，j(t)进行与式(2)形式一致的线性展开，再与式(5)合并之后得到通过校准磁通门数据和改变补偿系数形式而建立的飞机干扰磁场模型，得到通过校准磁通门数据和改变补偿系数形式之后飞机平台产生的干扰磁场B_interference(t)：

其中，h_i为由磁通门零漂和非正交性导致的磁干扰的补偿系数；

4)极化电流产生的干扰磁场：

在环境场中，不仅有磁场还有电场，经典的磁补偿T-L模型里并未考虑到环境电场产生的极化电流因素，因此忽略了极化电流产生的一部分干扰磁场；

海洋上空的环境电场E_env(r，t)为垂直指向海面的，则

E₀为环境电场的大小，e(t)为环境电场方向；

为由整体的飞机姿态组合成的矩阵，整体的飞机姿态以方向余弦u(t)表示，因此把

写成关于u(t)的函数，于是把e(t)用u(t)来表示：e(t)＝e(u(t))；

环境电场E_env(r，t)使飞机的金属机壳产生正比于

的极化电流，这个极化电流

产生的干扰磁场B_ele(t)表示为：

其中，

为环境电场方向e(t)的微分，

为环境电场方向的微分

的第j分量(j＝1，2，3)，

表达为：

其中，α_j为常数项系数，β_j，k为一次项系数，γ_j，k，l为二次项系数，M_i，j

为极化电流产生的干扰磁场系数，将该表达式与极化电流产生的干扰磁场结合，得到极化电流产生的干扰磁场B_ele(t)如下：

其中，σ_i，k为极化电流磁场干扰的补偿系数，τ_i，k，l为飞机感应磁场与极化电流磁场的耦合的干扰磁场的补偿系数；

5)对地磁场进行补偿：

采用国际地磁参考场(IGRF)，用IGRF模型算出地磁场值后，再用光泵测量值减去地磁场值，得到补偿后的地磁场值；

6)建立光泵磁力仪的方向误差：

在实际实验中，会用到光泵磁力仪测量的总磁场，光泵磁力仪测量的总磁场都存在方向误差磁场，导致求解系数会不准确，因此分别建立单光系或多光系光泵磁力仪的方向误差B_h(t)；

建立单光系光泵磁力仪的方向误差，包括以下步骤：

于单光系光泵磁力仪，主要用到的是铯光泵磁力仪，铯光泵磁力仪在转动时，其内部光轴与地磁场之间的夹角发生了改变，因此可把方向误差写成关于方向余弦u(t)的函数，经过泰勒函数展开得到的方向误差B_h(t)表达式为：

B_h(t)＝k₁*u₁(t)+k₂*u₂(t)+k₃*u₃(t)+k₄*u₁ ²(t)+k₅*u₂ ²(t)+k₆*u₃ ²(t) +k₇*u₁(t)*u₂(t)+k₈*u₁(t)*u₃(t)+k₉*u₂(t)*u₃(t)

其中，k₁，k₂，k₃，k₄，k₅，k₆，k₇，k₈，k₉为9个待求解的方向误差系数；

多光系光泵磁力仪的方向误差B_h表达式为：

其中，z_i为方向误差系数，n为光泵原子气室个数，n＝3，分别沿着x，y，z轴方向；每个气室上的交变磁场线圈方向与光传播方向垂直，分别沿着y轴、x轴和x轴方向， h_i为与方向误差有关的方向余弦，表达式为：

其中，η_i为原子气室安装方向与地磁场的夹角，η_i＝X或Y或Z，ε_i为原子气室上加的交变磁场方向与地磁场的夹角，ε_i＝X或Y或Z；

7)建立最终的飞机平台磁干扰模型：

将光泵磁力仪的方向误差加入到通过校准磁通门数据、改变补偿系数形式和增加极化电流磁场而建立的飞机干扰磁场模型中，建立最终的飞机平台磁干扰模型，得到最终的飞机干扰磁场B_f(t)：

B_f(t)＝B_interference(t)+B_ele(t)+B_h(t)；

上式简化为：

B_f(t)＝A*N (7)

其中，A为方向余弦组成的基函数矩阵，N为待求的磁补偿系数向量；

8)使用总体最小二乘方法求解出磁补偿系数向量：

a)A和b分别为不可观测的无误差基函数数据矩阵和无误差飞机干扰磁场数据向量，实际观测的基函数数据矩阵A₀和飞机干扰磁场数据向量b₀分别为：

A₀＝A+ΔA，b₀＝b+Δb

其中，ΔA和Δb分别为误差数据矩阵和误差数据向量；根据式(7)得到飞机干扰数据向量和基函数数据矩阵的关系为：A*N＝b，N为待求的磁补偿系数，则若存在误差，该式的形式为：

即(A+ΔA)N＝b+Δb，等价为(B+D)q＝0；

式中，增广数据矩阵B＝[A，b]和增广校正矩阵D＝[ΔA，Δb]均为m×(k+1)维矩阵，q为总体最小二乘的解向量

为(k+1)×1，k为无误差基函数数据矩阵A的列数，m为无误差基函数数据矩阵A的行数；

b)计算增广数据矩阵B的奇异值分解并存储右奇异矩阵V；

c)确定增广数据矩阵B的有效秩p；

d)总体最小二乘的解向量

中，N为由p个线性独立的未知参数组成的列向量，因此总体最小二乘的求解就变成了k+1-p个总体最小二乘问题的求解： B_j ^(p)q＝0，j＝1，2...，k+1-p，B_j ^(p)为B的第j列到第p+j列组成的子矩阵；

e)根据最小二乘原理，总体最小二乘解等价于使代价函数极小化，代价函数f(q)定义为：

其中

其中，u_j为增广数据矩阵B左奇异矩阵第i行第j列上的元素，

为右奇异矩阵V 第i行第j列上的元素；δ_j为增广数据矩阵B的奇异值；

f)定义

则代价函数f(q)写成f(q)＝q^TS^(p)q的形式，经过整理得到：

g)代价函数f(q)的极小变化量是f(q)对q求导并令其导数等于0，求得S^(p)q＝λe₁，λ为常数表示误能量，e₁为单位向量，e₁＝[1，0，0...，0]^T；

h)求S^(p)的逆矩阵S^-(p)；

i)对上式S^(p)q＝λe₁求解，求得总体最小二乘的解向量q，因为

则磁补偿系数向量N＝[x_tls(1)，x_tls(2)，...，x_tls(p)]^T最优最小二乘近似解为：

其中，S^-(p)为S^(p)的逆矩阵，x_tls(i)为磁补偿系数向量N的第i个元素。；

9)将求解出的补偿系数，用于飞机平台的实时补偿。

最后需要注意的是，公布实施例的目的在于帮助进一步理解本发明，但是本领域的技术人员可以理解：在不脱离本发明及所附的权利要求的精神和范围内，各种替换和修改都是可能的。因此，本发明不应局限于实施例所公开的内容，本发明要求保护的范围以权利要求书界定的范围为准。

Claims (5)

1.一种考虑飞机非刚性以及极化电流因素的航空磁补偿方法，其特征在于，所述航空磁补偿方法包括以下步骤：

1)建立经典的磁补偿T-L模型：

以飞机平台建立坐标系，建立经典的磁补偿T-L模型，得到飞机平台产生的干扰磁场B_TL(t)：

其中，a_i为飞机永久干扰磁场的补偿系数，b_i,j为飞机感应干扰磁场的补偿系数，c_i,j为飞机涡流干扰磁场的补偿系数，u_i(t)和u_j(t)为飞机平台坐标系与地磁场的方向余弦u(t)的分量，i＝1,2,3，j＝1,2,3，在飞机平台坐标系下分别为x、y和z轴，把它们展开写为u₁(t)＝cosX，u₂(t)＝cosY，u₃(t)＝cosZ，其中X、Y和Z分别为飞机平台坐标系x、y和z轴与地磁场的夹角，cosX、cos Y和cosZ表达式分别如下：

其中，B₁(t)、B₂(t)和B₃(t)分别为在t时刻三分量磁通门磁力仪测量的沿飞机平台坐标系x、y和z轴的地磁场测量值；

2)建立飞机平台干扰磁场模型：

经典的T-L模型是在飞机平台为刚体这个基本假设条件建立的，但是飞机在转弯和做机动动作的时候，方向舵、升降舵和副翼这些活动部件发生了状态变化，机翼发生了弯曲，因此飞机在做机动动作的时候不能看成是一个理想刚体，当飞机在做侧滚、俯仰和偏航机动动作时，飞机的结构状况采用运动状态参数来描述，飞机平台坐标系下地磁场的方向余弦u(t)能够描述整体的飞机姿态，但是不能描述具体的各个活动部件的状态；

a)以飞机平台坐标系下地磁场的方向余弦的微分

反映活动部件的运动状态，把式(1)中的补偿系数写成u(t)和

的泛函：

b)将补偿系数以近似的线性化表示为：

a_i ⁽⁰⁾，b_i,j ⁽⁰⁾和c_i,j ⁽⁰⁾分别表示补偿系数a_i、b_i,j和c_i,j的常数项分量，a_i,j ⁽¹⁾为补偿系数a_i关于u_j(t)的线性系数，b_i,j,k ⁽¹⁾为补偿系数b_i,j关于u_k(t)的线性系数，c_i,j,k ⁽¹⁾为补偿系数c_i,j关于u_k(t)的线性系数；a_i,j ⁽²⁾为补偿系数a_i关于

的线性系数，b_i,j,k ⁽²⁾为补偿系数b_i,j关于

的线性系数；

c)建立飞机平台干扰磁场模型，得到修改补偿系数之后飞机平台产生的干扰磁场B_coe(t)：

其中，d_i,j,k为飞机主体感应磁场与活动部件永久磁场的耦合的补偿系数，g_i,j,k为感应磁场与活动部件涡流磁场的耦合的补偿系数；

3)建立三分量磁通门磁力仪误差模型的磁补偿模型：

三分量磁通门磁力仪存在非正交性和零漂问题，地磁场的真实值为：

其中，

为磁通门磁力仪测量的沿飞机平台坐标系x、y和z轴的地磁场真实值，

为修正矩阵标记为w，该矩阵是未知的，用来修正磁通门存在的非正交性问题，

为误差向量，也是未知的向量，用来修正三分量磁通门磁力仪存在的零漂问题；

将式(4)写成简洁形式为：

因此在飞机平台坐标系下，修正后的飞机平台坐标系下地磁场的方向余弦

表示为：

其中，σ(t)为时间缓变误差函数，

W(t)表示为不同时刻t的修正矩阵；将

代入式(1)，且

是对W(t)的微分，整理后得到包括修正磁通门数据建立的磁补偿模型，得到修正磁通门数据后飞机平台产生的干扰磁场B_cal(t)：

其中，a′_i为修正后的永久磁场的补偿系数，b′_i,j为修正后的感应磁场的补偿系数，c′_i,j为修正后的涡流磁场的补偿系数，h_i为由磁通门零漂和非正交性导致的磁干扰的补偿系数，β(t)是随时间变化的一次函数；

由于飞机的运动姿态也是随时间变化的，因此对式(5)中的修正后的永久磁场的补偿系数a′_i、修正后的感应磁场的补偿系数b′_i,j和修正后的涡流磁场的补偿系数c′_i,j进行与式(2)形式一致的线性展开，再与式(5)合并之后得到通过校准磁通门数据和改变补偿系数形式而建立的飞机干扰磁场模型，得到通过校准磁通门数据和改变补偿系数形式之后飞机平台产生的干扰磁场B_interference(t)：

其中，h_i为由磁通门零漂和非正交性导致的磁干扰的补偿系数；

4)得到极化电流产生的干扰磁场：

在环境场中，不仅有磁场还有电场，经典的磁补偿T-L模型里并未考虑到环境电场产生的极化电流因素，因此忽略了极化电流产生的一部分干扰磁场；

海洋上空的环境电场E_env(r,t)为垂直指向海面的，则

E₀为环境电场的大小，e(t)为环境电场方向，

表示垂直指向海面的方向向量；

为由整体的飞机姿态组合成的矩阵，整体的飞机姿态以飞机平台坐标系下地磁场的方向余弦u(t)表示，因此把

写成关于u(t)的函数，于是把e(t)用u(t)来表示：e(t)＝e(u(t))；

环境电场E_env(r,t)使飞机的金属机壳产生正比于

的极化电流，这个极化电流产生的干扰磁场B_ele(t)表示为：

其中，

为环境电场方向e(t)的微分，

为环境电场方向的微分

的j分量，

为极化电流的系数，j＝1,2,3，

表达为：

其中，α_j为常数项系数，β_j,k为一次项系数，γ_j,k,l为二次项系数，M_i,j为极化电流产生的干扰磁场系数，将该表达式与极化电流产生的干扰磁场结合，得到极化电流产生的干扰磁场B_ele(t)：

其中，σ_i,k为极化电流磁场干扰的补偿系数，τ_i,k,l为飞机感应磁场与极化电流磁场的耦合的干扰磁场的补偿系数；

5)对地磁场进行补偿：

采用国际地磁参考场IGRF，用IGRF模型算出地磁场值后，再用光泵磁力仪测量的磁场值减去地磁场值，得到补偿掉地磁场后的剩余干扰磁场值；

6)建立光泵磁力仪的方向误差：

在实际实验中，会用到光泵磁力仪测量的总磁场，光泵磁力仪测量的总磁场都存在方向误差磁场，因此分别建立单光系或多光系光泵磁力仪的方向误差B_h(t)；

7)建立最终的飞机平台磁干扰模型：

将光泵磁力仪的方向误差加入到通过校准磁通门数据、改变补偿系数形式和增加极化电流磁场而建立的飞机干扰磁场模型中，建立最终的飞机平台磁干扰模型，得到最终的飞机干扰磁场B_f(t)：

B_f(t)＝B_interference(t)+B_ele(t)+B_h(t)

上式简化为：

B_f(t)＝A(t)*N (7)

其中，A(t)为飞机平台坐标系下地磁场的方向余弦组成的基函数矩阵，N为待求的磁补偿系数向量；

8)使用总体最小二乘方法求解出磁补偿系数向量；

9)将求解出的磁补偿系数向量，用于飞机平台的实时补偿。

2.如权利要求1所述的航空磁补偿方法，其特征在于，在步骤6)中，建立单光系光泵磁力仪的方向误差，包括以下步骤：

单光系光泵磁力仪采用铯光泵磁力仪，铯光泵磁力仪在转动时，内部光轴与地磁场之间的夹角发生了改变，因此把方向误差写成关于飞机平台坐标系下地磁场的方向余弦u(t)的函数，经过泰勒函数展开得到的方向误差B_h(t)表达式为：

B_h(t)＝k₁*u₁(t)+k₂*u₂(t)+k₃*u₃(t)+k₄*u₁ ²(t)+k₅*u₂ ²(t)+k₆*u₃ ²(t)+k₇*u₁(t)*u₂(t)+k₈*u₁(t)*u₃(t)+k₉*u₂(t)*u₃(t)

其中，k₁,k₂,k₃,k₄,k₅,k₆,k₇,k₈,k₉为9个待求解的方向误差系数。

3.如权利要求1所述的航空磁补偿方法，其特征在于，在步骤6)中，建立多光系光泵磁力仪的方向误差，包括以下步骤：

多光系光泵磁力仪的方向误差B_h表达式为：

其中，z_i为方向误差系数，n为光泵原子气室个数，n＝3，分别沿着x,y,z轴方向；每个气室上的交变磁场线圈方向与光传播方向垂直，分别沿着y轴、x轴和x轴方向，h_i为与方向误差有关的方向余弦，表达式为：

其中，η_i是原子气室安装方向与地磁场的夹角，η_i＝X或Y或Z，ε_i是原子气室上加的交变磁场方向与地磁场的夹角，ε_i＝X或Y或Z。

4.如权利要求1所述的航空磁补偿方法，其特征在于，在步骤8)中，使用总体最小二乘方法求解出磁补偿系数向量，包括以下步骤：

a)A和b分别为不可观测的无误差基函数数据矩阵和无误差飞机干扰磁场数据向量，实际观测的基函数数据矩阵A₀和飞机干扰磁场数据向量b₀分别为：

A₀＝A+ΔA，b₀＝b+Δb

其中，ΔA和Δb分别为误差数据矩阵和误差数据向量；根据式(7)得到飞机干扰数据向量和基函数数据矩阵的关系为：A*N＝b，N为待求的磁补偿系数向量，则若存在误差，该式的形式为：

即(A+ΔA)N＝b+Δb，等价为((B+D)q＝0；

式中，增广数据矩阵B＝[A,b]和增广校正矩阵D＝[ΔA,Δb]均为m×(k+1)维矩阵，q为总体最小二乘的解向量，而

为(k+1)×1向量，k为无误差基函数数据矩阵A的列数，m为无误差基函数数据矩阵A的行数；

b)计算增广数据矩阵B的奇异值分解并存储右奇异矩阵V；

c)确定增广数据矩阵B的有效秩p；

d)总体最小二乘的解向量

中，N为由p个线性独立的未知参数组成的列向量，因此总体最小二乘的求解就变成了k+1-p个总体最小二乘问题的求解：B_j ^(p)q＝0,j＝1,2...,k+1-p，B_j ^(p)为B的第j列到第p+j列组成的子矩阵；

e)根据最小二乘原理,总体最小二乘解等价于使代价函数极小化，代价函数f(q)定义为：

其中

其中，u_j为增广数据矩阵B左奇异矩阵第i行第j列上的元素，

为右奇异矩阵V第i行第j列上的元素；δ_j为增广数据矩阵B的奇异值；

f)定义

则代价函数f(q)写成f(q)＝q^TS^(p)q的形式，经过整理得到：

g)代价函数f(q)的极小变化量是f(q)对q求导并令导数等于0，求得S^(p)q＝λe₁，λ为常数表示误能量，e₁为单位向量，e₁＝[1,0,0...,0]^T；

h)求S^(p)的逆矩阵S^-(p)；

i)对上式S^(p)q＝λe₁求解，求得总体最小二乘的解向量q，因为

则磁补偿系数向量N＝[x_tls(1),x_tls(2)，...,x_tls(p)]^T最优最小二乘近似解为：

其中，S^-(p)为S^(p)的逆矩阵，x_tls(i)为磁补偿系数向量N的第i个元素。

5.如权利要求1所述的航空磁补偿方法，其特征在于，进一步，在飞机进行校准飞行前，在飞机上安装光泵磁力仪、三分量磁通门磁力仪、高度计和GPS；光泵磁力仪用来测量的总磁场，三分量磁通门磁力仪实时测量沿飞机平台坐标系x、y和z轴的地磁场测量值，高度计用来测量飞机飞行的高度，GPS用来测量飞机所处位置的经度和纬度。

Images