CN105891755A - 航空器挂载式磁通门磁梯度张量仪的校正方法 - Google Patents

Abstract

本发明涉及一种航空器挂载式磁通门磁梯度张量仪的校正方法，一建立单个磁通门误差模型，二建立磁梯度张量分量误差校正模型，三在高空中磁场均匀的区域采集校正数据，四将上述张量分量的误差校正，五在地面磁场不均匀的区域采集校正数据，六求解张量分量的标度因子误差，七结合高空校正和地面校正所获得校正系数，八完整计算校正后的磁梯度张量。本发明的校正算法与现有的航空磁场梯度张量校正算法相比，不仅仅对单个张量分量进行校正，还包含了张量作为一个整体的校正，从而使得张量分量和整体的校正结果更加准确，同时其应用条件更符合实际的航空磁梯度张量探测区域的地质条件。

Description

航空器挂载式磁通门磁梯度张量仪的校正方法

技术领域：

本发明涉及一种航空物探磁测仪器校正，尤其是航空器挂载式磁通门磁梯度张量仪的校正方法。

背景技术：

磁梯度张量仪测量磁场三分量沿正交坐标系的空间变化率。磁梯度张量仪受地磁影响小，获取信息丰富，特别适合发现浅层潜伏矿，磁性运动目标等，是磁测仪器发展方向之一。在航空磁梯度张量测量中,为了避免张量仪受到航空器如直升机等的磁干扰，通常是把张量仪挂载到离航空器较远的地方，然后在飞行过程中测量目标的磁场梯度，本发明正是针对此类应用时张量仪的校正方法。

由于需要对磁场矢量敏感的传感器，故目前张量仪的实现主要有两种方式：一种是基于超导量子干涉器件的张量仪，如德国的航空超导磁梯度张量仪[R.Stolz et al,Magnetic full-tensor SQUID gradiometer system for geophysical applications,The Leading Edge,2006,25(2):178-180.]，澳大利亚的航空超导磁梯度张量仪[P.Schmidt et al,GETMAG-a SQUID magnetic tensor gradiometer for mineral andoil exploration,Exploration Geophysics,2004,35:297-305.]等。超导张量仪虽然具有较高的分辨率，但是价格昂贵，且不能在常温下使用，应用领域有较大的局限性。另一种是基于磁通门传感器构造的张量仪，如美国地质调查局的正四面体形磁通门张量仪[PJ.Brown et al,A case study of magnetic gradient tensor invariants appliedto the UXO problem,U.S.Geological Survey，2004：1-4.]。美国IBM公司研制的磁通门张量仪[RH.Koch et al,Room temperature three sensor magnetic field gradiometer,Review of Scientific Instruments,1996,64(1):230-235.]。

虽然磁通门构造的张量仪具有成本低，温度范围广，分辨率较高等优点，但是磁通门存在三个感应轴标度因子不完全一致造成的误差、三个感应轴不严格正交造成的误差、零点存在偏移、动态误差等问题。这些误差以及各个磁通门误差特性不一致造成磁通门张量仪误差大的问题，从而严重影响了航空磁梯度张量仪的应用效果。针对此问题，目前存在两种校正方法。第一种为实验室标定法，一般是在磁屏蔽环境中，利用亥姆赫兹线圈或麦克斯韦线圈，以及标准电流源等装置进行标定。这类方法除了标定环境难以建立以外，最大问题是标定所处的环境(如屏蔽环境中磁场很小，测量装置处于静止状态)和实际的航空磁梯度张量测量环境(如存在地球磁场，测量装置处于运动状态)相差甚远，直接应用实验室内获得的校正参数并不能保证航空测量应用中校正的正确性。第二种方法是在地面上找一块磁场均匀的区域，然后以多个磁通门测量的总场值或磁场分量值相等为条件求解线性或非线性方程以获得校正参数。这类方法有三个主要的问题：第一，在磁场均匀区域内，单个梯度分量为零是磁梯度张量仪校正中的必要条件，而非充要条件。既不能保证对非零梯度情况响应的正确性，也不能保证各个梯度分量形成张量后作为整体的正确性。第二，在地面难以复现飞行过程中磁通门的动态特性，故难以辨识出多个磁通门动态特性不一致性的问题。而航空磁梯度张量测量始终处于运动状态，动态特性的不一致同样也是形成运动噪声，淹没目标磁探测信号的原因之一。第三，实际航空磁场测量的目标区域中或附近，地面磁性异常往往较大，难以找到磁场均匀的区域进行校正。

航空挂载式磁通门磁梯度张量仪中磁通门的测量误差来源于以下几项：磁通门传感器三个感应轴标度因子不完全一致造成的误差、三个感应轴不严格正交造成的误差、零点存在偏移、由机械安装误差造成的非对准误差、动态特性不一致造成的误差。

发明内容：

本发明的目的是针对上述现有技术的不足，提供一种提高其测量准确性的航空器挂载式磁通门磁梯度张量仪的校正方法。

本发明的目的是通过以下技术方案实现的：

航空器挂载式磁通门磁梯度张量仪的校正方法，包括以下步骤：

步骤一、建立单个磁通门误差模型：

把惯性导航***和磁通门张量仪进行刚体连接，把各个磁通门传感器坐标系x,y,z轴分别向惯性导航坐标系的x,y,z轴方向对齐，则单个磁通门误差模型为：

$D_j\frac{d^{j}T}{{\mathrm{dt}}^j}+D_{j-1}\frac{d^{j-1}T}{{\mathrm{dt}}^{j-1}}+...+D_1\frac{dT}{dt}+D_{0}T={\mathrm{FEM}}^{-1}B+R+\mathrm{\Φ}---\left(1\right)$

式中：F为标度因子误差矩阵；E为非正交误差矩阵；R为零位误差矢量；M为磁通门坐标系与惯性导航坐标系非对准误差矩阵；D_j为磁通门测量值对时间的j阶导数的描述，即动态特性描述矩阵；T为磁通门在自身测量坐标系下的实测值；B为磁通门在惯导坐标系下的真实磁场值；Φ为磁通门的测量噪声；

每个误差矩阵或矢量的具体表达式如下：

(1)标度因子误差矩阵：

$F=\left[\begin{array}{ccc}1+S_x& 0& 0\\ 0& 1+S_y& 0\\ 0& 0& 1+S_z\end{array}\right]---\left(2\right)$

其中S_x,S_y,S_z表示磁通门感应轴的标度因子误差，

(2)非正交误差矩阵：

假设实际磁通门坐标系x轴与理想磁通门坐标系x轴的方向一致，ρ为实际磁通门坐标系y轴与理想磁通门坐标系y轴的夹角；为实际磁通门坐标系z轴与理想磁通门坐标系z轴的夹角在磁通门xoz平面上的投影；λ为实际磁通门坐标系z轴与理想磁通门坐标系z轴的夹角在磁通门yoz平面上的投影，

(3)非对准误差矩阵：

$M=\left[\begin{array}{ccc}{c}_{\mathrm{\β}}{c}_{\mathrm{\γ}}& -c_{\mathrm{\α}}{s}_{\mathrm{\γ}}+s_{\mathrm{\α}}{s}_{\mathrm{\β}}{c}_{\mathrm{\γ}}& s_{\mathrm{\α}}{s}_{\mathrm{\γ}}+c_{\mathrm{\α}}{s}_{\mathrm{\β}}{c}_{\mathrm{\γ}}\\ c_{\mathrm{\β}}{s}_{\mathrm{\γ}}& c_{\mathrm{\α}}{c}_{\mathrm{\γ}}+s_{\mathrm{\α}}{s}_{\mathrm{\β}}{s}_{\mathrm{\γ}}& -s_{\mathrm{\α}}{c}_{\mathrm{\γ}}+c_{\mathrm{\α}}{s}_{\mathrm{\β}}{s}_{\mathrm{\γ}}\\ -s_{\mathrm{\β}}& s_{\mathrm{\α}}{c}_{\mathrm{\β}}& c_{\mathrm{\α}}{c}_{\mathrm{\β}}\end{array}\right]---\left(4\right)$

式中c_α＝cosα，c_β＝cosβ，c_γ＝cosγ，s_α＝sinα，s_β＝sinβ，s_γ＝sinγ,α、β和γ分别为实际磁通门坐标系相对惯导坐标系，绕实际磁通门坐标系x轴,y轴和z轴的非对准角。当非对准角非常小的时候可以将非对准误差矩阵近似表示为：

$M=\left[\begin{array}{ccc}1& -\mathrm{\γ}& \mathrm{\β}\\ \mathrm{\γ}& 1& -\mathrm{\α}\\ -\mathrm{\β}& \mathrm{\α}& 1\end{array}\right]---\left(5\right)$

(4)零位误差矢量：

$R=\left[\begin{array}{c}{\mathrm{\μ}}_x\\ {\mathrm{\μ}}_y\\ {\mathrm{\μ}}_z\end{array}\right]---\left(6\right)$

式中μ_x,μ_y,μ_z分别为磁通门三个轴分别存在的零位误差，

(5)动态特性描述矩阵：

$D_j=\left[\begin{array}{ccc}{d}_{11}& d_{12}& d_{13}\\ d_{21}& d_{22}& d_{23}\\ d_{31}& d_{32}& d_{33}\end{array}\right]---\left(7\right)$

当磁通门由一阶***描述时取j＝1，D₁为一阶导数前面的系数矩阵。

将所述三轴磁通门的误差模型改写为校正模型：

$B=\left[\begin{array}{c}{b}_x\\ b_y\\ b_z\end{array}\right]=C_{0}T+\underset{j=1}{\overset{h}{\Σ}}\frac{d^j{C}_{j}T}{{\mathrm{dt}}^j}+O---\left(8\right)$

式中：C₀为T第0阶导数的校正参数矩阵；C_j为Τ第j阶导数的校正参数矩阵；O为零位误差校正矢量；b_x,b_y,b_z分别表示磁通门在惯性导航坐标系x,y,z轴上的磁场值，则磁通门中单个磁场分量误差校正模型可表示为：

$b_u=P_{u0}T+\underset{j=1}{\overset{h}{\Σ}}\frac{d^j{Q}_{uj}T}{{\mathrm{dt}}^j}+O_u---\left(9\right)$

式中：u表示x,y,z中一个轴，P_u0为矩阵C₀中对应b_u的行向量校正参数；Q_uj为矩阵C_j中对应b_u的行向量校正参数；O_u为总偏差矢量O_n中对应b_u的校正参数；

步骤二、建立磁梯度张量分量误差校正模型：

磁场沿x,y,z三个方向的二阶导数可构成磁梯度张量G，尽管G有九个分量，由于在不包括场源的域内，磁场的散度和旋度为0，即对称，且迹为0，故只有五个分量是独立的，可表示为：

$G=\left(\begin{array}{ccc}\frac{\∂B_x}{\∂x}& \frac{\∂B_y}{\∂x}& \frac{\∂B_z}{\∂x}\\ \frac{\∂B_x}{\∂y}& \frac{\∂B_y}{\∂y}& \frac{\∂B_z}{\∂y}\\ \frac{\∂B_x}{\∂z}& \frac{\∂B_y}{\∂z}& \frac{\∂B_z}{\∂z}\end{array}\right)=\left(\begin{array}{ccc}{g}_{xx}& g_{yx}& g_{zx}\\ g_{yx}& -\left(g_{xx}+g_{zz}\right)& g_{yz}\\ g_{xz}& g_{yz}& g_{zz}\end{array}\right)---\left(10\right)$

磁通门张量仪由不同方式构建，主要有十字、四方体和三角形等结构，但都是利用差分近似微分的方式来表示磁梯度张量G各个分量。以十字型结构为例，在惯性导航坐标系下，张量分量表示方式为：

$g_{uv}=\frac{{\mathrm{\ΔB}}_u}{{\mathrm{\Δs}}_v}=\frac{b_{nu}-b_{mu}}{{\mathrm{\Δs}}_v}=\frac{\left(P_{nu0}{T}_n+\underset{j=1}{\overset{h}{\mathrm{\Σ}}}\frac{d^j{Q}_{nuj}{T}_n}{{\mathrm{dt}}^j}+O_{nu}\right)-\left(P_{mu0}{T}_m+\underset{j=1}{\overset{h}{\mathrm{\Σ}}}\frac{d^j{Q}_{muj}{T}_m}{{\mathrm{dt}}^j}+O_{mu}\right)}{{\mathrm{\Δs}}_v}---\left(11\right)$

式中：m,n表示磁通门标号，取1,2,3,4；u,v都代表坐标轴的方向，取x,y,z；g_uv为惯性导航坐标系下某个张量分量；ΔB_u为惯性导航坐标系下沿u轴方向的磁场变化量；Δs_v表示磁通门n与磁通门m之间的基线距离。

将所述张量分量误差模型改写为：

$g_{uv}=P_{nu0}\frac{T_n}{{\mathrm{\Δs}}_v}-P_{mu0}\frac{T_m}{{\mathrm{\Δs}}_v}+\underset{j=1}{\overset{h}{\mathrm{\Σ}}}\left(\frac{d^j{Q}_{nuj}\frac{T_n}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}-\frac{d^j{Q}_{muj}\frac{T_m}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}\right)+O_{uv}---\left(12\right)$

式中：O_uv为张量分量g_uv的总偏置误差；

步骤三、在高空中磁场均匀的区域采集校正数据：

将磁通门磁梯度张量仪挂载到航空器如直升机下方较远处。航空器飞行到高空，然后以方形，八字形等多种轨迹进行正常飞行，张量仪采集不同姿态下带误差的磁场测量数据；

步骤四、因为高空中地球磁场的梯度近似为0，设张量分量g_uv误差模型中磁通门n中沿u轴方向测量值前面的校正参数为p_uv，将上述张量分量的误差校正模型改写为：

$\frac{g_{uv}}{p_{uv}}=\frac{P_{nu0}}{p_{uv}}\frac{T_n}{{\mathrm{\Δs}}_v}-\frac{P_{mu0}}{p_{uv}}\frac{T_m}{{\mathrm{\Δs}}_v}+\underset{j=1}{\overset{h}{\mathrm{\Σ}}}\left(\frac{d^j\frac{Q_{nuj}}{p_{uv}}\frac{T_n}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}-\frac{d^j\frac{Q_{muj}}{p_{uv}}\frac{T_m}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}\right)+\frac{O_{uv}}{p_{uv}}=0---\left(13\right)$

故式(13)中有一个测量值的校正参数为1，这样以此项测量值为已知值，对式(13)列写线性方程组，通过最小二乘拟合方法解算线性方程组，可以获得式(13)中的各项参数。此时解算出各个磁梯度张量分量误差模型中参数的相对量， $\frac{P_{\mathrm{nu}0}}{p_{\mathrm{uv}}},\frac{P_{\mathrm{mu}0}}{p_{\mathrm{uv}}},\frac{Q_{\mathrm{nuj}}}{p_{\mathrm{uv}}},\frac{Q_{\mathrm{muj}}}{p_{\mathrm{uv}}},\frac{O_{\mathrm{uv}}}{p_{\mathrm{uv}}};$

步骤五、在地面磁场不均匀的区域采集校正数据：

在地面上有磁梯度的区域，把张量仪的传感器部分绕其中心点进行旋转，记录此过程的测量值；

步骤六、以不变量进行约束，求解张量分量的标度因子误差：

将高空校正之后所获得的各项校正系数重新带入到公式(13)中，可以获得各个张量分量的相对量值，表达为：

由于式(1)中误差矩阵的数值并不相同，因此每个分量所对应的p_uv也各不相同，故此时的张量分量可以认为存在张量分量标度因子误差。张量的不变量之一C_T，可以由即张量缩并运算获得，其特性是在坐标旋转时，一直保持值不变，当张量分量存在标度因子误差时，旋转时所获得的C_T会形成明显的波动，因此旋转过程中C_T可以用来校正张量分量的标度因子误差，

把张量分量的校正模型重新表达为：

则张量不变量校正模型表达为： $C_T^2=\underset{u,v=(x,y,z)}{\Σ}{g_{uv}}^2=\underset{u,v=(x,y,z)}{\Σ}{p_{uv}}^2{g_{{\mathrm{uv}}_{mid}}}^2---\left(16\right)$

将公式(16)中的p_uv都取为1，计算出此时的C_T值并进行平均处理。最后，C_T的平均值和作为公式(16)中已知量，采用广义最小二乘拟合可估计出p_uv；

步骤七、结合高空校正获得的各项系数，可以求得每个张量分量的最终校正参数：P_nu0,P_mu0,Q_nuj,Q_muj,O_uv；

步骤八、将三轴磁通门磁场测量值和校正参数绝对量输入到磁梯度张量分量校正模型中，计算校正后的磁梯度张量。

有益效果：1.与现有磁梯度张量仪实验室校正方法相比，本发明首先不用在野外构造磁屏蔽环境及应用昂贵的校正设备，其次本校正方法中数据采集的环境与实际航空磁测的环境完全一致，克服了实验室校正方法由于环境不一致造成的校正参数变化的问题。2.与现有磁梯度张量仪野外校正方法相比，本发明提出的校正方法多了在地面利用张量旋转不变量进行校正的过程。此过程是把张量各个分量形成一个整体后进行校正，不但使得单个分量更加准确，也保证了张量作为一个整体的正确性。3.本发明应用中不是在地面寻找磁场均匀区域，而是寻找磁场非均匀区域。由于航空磁测的主要应用是发现铁矿等磁性异常，因此在航空磁测的目标区域中或附近一般难以存在磁场均匀区域。故本发明的使用条件更符合实际的航空磁测探测区域的地质条件，即受环境限制的因素更少。未校正前，均方根误差为785.97，校正之后均方根误差为0.57，改善比为1380.58。其它张量分量的校正效果类似，从而验证了本校正方法的正确性和显著的效果。

附图说明：

图1为航空器挂载式磁通门磁梯度张量仪的校正方法流程图

图2十字全张量仪安装示意图

图3张量分量之一g_zz在校正前后与真实值的对比图

图4测线上校正前和校正后g_zx和g_xz的对比图。

建立“十字”型磁通门张量仪，选取第一个磁通门1坐标系为o₀x₀y₀z₀6，第二个磁通门2坐标系为o₁x₁y₁z₁7，第三个磁通门3坐标系为o₂x₂y₂z₂8，第四个磁通门4坐标系为o₃x₃y₃z₃9，惯性导航***(INS)5的坐标系为oxyz 10，磁通门1和磁通门3放置在坐标系10的x轴上关于原点o对称，间距为Δx，采用无磁性材料连接。磁通门2和磁通门4放置在坐标系10的z轴上关于原点o对称，间距为Δz，采用无磁性材料连接。惯性导航***5放置在y轴上，距离原点o间距为Δy，采用无磁性材料与其它部件连接。安装时尽量保证磁通门1坐标系o₀x₀y₀z₀6、磁通门2坐标系o₁x₁y₁z₁7、磁通门3坐标系o₂x₂y₂z₂8和磁通门4坐标系o₃x₃y₃z₃9与惯性导航***5坐标系oxyz 10相同。上述连接构成“十字”型磁通门张量仪。

具体实施方式：

下面结合附图和实施例对本发明做进一步的详细说明。

本发明利用高空中地球磁场的梯度张量近似为零的特点进行校正，然后根据地面磁梯度非零区域中张量仪旋转时不变量的可以指示各个梯度分量标度因子不一致的特点，计算出高空校正后未确定的各个梯度分量的标度因子，最终达到校正磁通门张量仪的目的。

本发明以单个磁通门误差模型为理论基础，推导出磁通门张量仪的误差模型，然后利用高空中地球磁场的梯度张量近似为零的特点，以及张量形成的不变量可以指示各个梯度分量标度因子不一致的特点，针对四个磁通门按“十”字结构排列所构成的磁梯度张量仪，提供一种由航空器挂载后进行测量应用时的校正方法。

航空挂载式磁通门磁梯度张量仪中磁通门的测量误差来源于以下几项：磁通门传感器三个感应轴标度因子不完全一致造成的误差、三个感应轴不严格正交造成的误差、零点存在偏移、由机械安装误差造成的非对准误差、动态特性不一致造成的误差。

在本实施例中以“十字”型磁通门梯度张量仪为例，图2为“十字”型磁通门梯度仪安装结构图，其中包括安装方式(四个磁通门，惯性导航***，刚性连接)和坐标系。

在具体应用实例中，可以通过仿真进行模拟。仿真中四个磁通门按十字形排列，基线距离为0.1米。地球磁场强度为50000nT，磁偏角为60°，磁倾角为-9°。

航空器挂载式磁通门磁梯度张量仪的校正方法，包括以下步骤：

步骤一、建立单个磁通门误差模型：

把惯性导航***和磁通门张量仪进行刚体连接，把各个磁通门传感器坐标系x,y,z轴分别向惯性导航坐标系的x,y,z轴方向对齐，则单个磁通门误差模型为：

$D_j\frac{d^{j}T}{{\mathrm{dt}}^j}+D_{j-1}\frac{d^{j-1}T}{{\mathrm{dt}}^{j-1}}+...+D_1\frac{dT}{dt}+D_{0}T={\mathrm{FEM}}^{-1}B+R+\mathrm{\Φ}---\left(1\right)$

式中：F为标度因子误差矩阵；E为非正交误差矩阵；R为零位误差矢量；M为磁通门坐标系与惯性导航坐标系非对准误差矩阵；D_j为磁通门测量值对时间的j阶导数的描述，即动态特性描述矩阵；T为磁通门在自身测量坐标系下的实测值；B为磁通门在惯导坐标系下的真实磁场值；Φ为磁通门的测量噪声；

每个误差矩阵或矢量的具体表达式如下：

(1)标度因子误差矩阵：

$F=\left[\begin{array}{ccc}1+S_x& 0& 0\\ 0& 1+S_y& 0\\ 0& 0& 1+S_z\end{array}\right]---\left(2\right)$

其中S_x,S_y,S_z表示磁通门感应轴的标度因子误差，

(2)非正交误差矩阵：

假设实际磁通门坐标系x轴与理想磁通门坐标系x轴的方向一致，ρ为实际磁通门坐标系y轴与理想磁通门坐标系y轴的夹角；为实际磁通门坐标系z轴与理想磁通门坐标系z轴的夹角在磁通门xoz平面上的投影；λ为实际磁通门坐标系z轴与理想磁通门坐标系z轴的夹角在磁通门yoz平面上的投影，

(3)非对准误差矩阵：

$M=\left[\begin{array}{ccc}{c}_{\mathrm{\β}}{c}_{\mathrm{\γ}}& -c_{\mathrm{\α}}{s}_{\mathrm{\γ}}+s_{\mathrm{\α}}{s}_{\mathrm{\β}}{c}_{\mathrm{\γ}}& s_{\mathrm{\α}}{s}_{\mathrm{\γ}}+c_{\mathrm{\α}}{s}_{\mathrm{\β}}{c}_{\mathrm{\γ}}\\ c_{\mathrm{\β}}{s}_{\mathrm{\γ}}& c_{\mathrm{\α}}{c}_{\mathrm{\γ}}+s_{\mathrm{\α}}{s}_{\mathrm{\β}}{s}_{\mathrm{\γ}}& -s_{\mathrm{\α}}{c}_{\mathrm{\γ}}+c_{\mathrm{\α}}{s}_{\mathrm{\β}}{s}_{\mathrm{\γ}}\\ -s_{\mathrm{\β}}& s_{\mathrm{\α}}{c}_{\mathrm{\β}}& c_{\mathrm{\α}}{c}_{\mathrm{\β}}\end{array}\right]---\left(4\right)$

式中c_α＝cosα，c_β＝cosβ，c_γ＝cosγ，s_α＝sinα，s_β＝sinβ，s_γ＝sinγ,α、β和γ分别为实际磁通门坐标系相对惯导坐标系，绕实际磁通门坐标系x轴,y轴和z轴的非对准角。当非对准角非常小的时候可以将非对准误差矩阵近似表示为：

$M=\left[\begin{array}{ccc}1& -\mathrm{\γ}& \mathrm{\β}\\ \mathrm{\γ}& 1& -\mathrm{\α}\\ -\mathrm{\β}& \mathrm{\α}& 1\end{array}\right]---\left(5\right)$

(4)零位误差矢量：

$R=\left[\begin{array}{c}{\mathrm{\μ}}_x\\ {\mathrm{\μ}}_y\\ {\mathrm{\μ}}_z\end{array}\right]---\left(6\right)$

式中μ_x,μ_y,μ_z分别为磁通门三个轴分别存在的零位误差，

(5)动态特性描述矩阵：

$D_j=\left[\begin{array}{ccc}{d}_{11}& d_{12}& d_{13}\\ d_{21}& d_{22}& d_{23}\\ d_{31}& d_{32}& d_{33}\end{array}\right]---\left(7\right)$

当磁通门由一阶***描述时取j＝1，D₁为一阶导数前面的系数矩阵。

将所述三轴磁通门的误差模型改写为校正模型：

$B=\left[\begin{array}{c}{b}_x\\ b_y\\ b_z\end{array}\right]=C_{0}T+\underset{j=1}{\overset{h}{\Σ}}\frac{d^j{C}_{j}T}{{\mathrm{dt}}^j}+O---\left(8\right)$

式中：C₀为T第0阶导数的校正参数矩阵；C_j为Τ第j阶导数的校正参数矩阵；O为零位误差校正矢量；b_x,b_y,b_z分别表示磁通门在惯性导航坐标系x,y,z轴上的磁场值，则磁通门中单个磁场分量误差校正模型可表示为：

$b_u=P_{u0}T+\underset{j=1}{\overset{h}{\Σ}}\frac{d^j{Q}_{uj}T}{{\mathrm{dt}}^j}+O_u---\left(9\right)$

式中：u表示x,y,z中一个轴，P_u0为矩阵C₀中对应b_u的行向量校正参数；Q_uj为矩阵C_j中对应b_u的行向量校正参数；O_u为总偏差矢量O_n中对应b_u的校正参数；

步骤二、建立磁梯度张量分量误差校正模型：

磁场沿x,y,z三个方向的二阶导数可构成磁梯度张量G，其中只有五个量是独立，利用差分近似微分的方式表示磁梯度张量各个分量。对于“十字”型的磁通门磁梯度张量仪在惯性导航坐标系下磁梯度张量G可表示为：(公式也换了)

$G=\left(\begin{array}{ccc}{g}_{xx}& g_{yx}& g_{zx}\\ g_{xy}& g_{yy}& g_{zy}\\ g_{xz}& g_{yz}& g_{zz}\end{array}\right)=\left(\begin{array}{ccc}{g}_{xx}& g_{yx}& g_{zx}\\ g_{yx}& -\left(g_{xx}+g_{zz}\right)& g_{yz}\\ g_{xz}& g_{yz}& g_{zz}\end{array}\right)=\left(\begin{array}{ccc}\frac{B_{0x}-B_{2x}}{\mathrm{\Δ}x}& \frac{B_{0y}-B_{2y}}{\mathrm{\Δ}x}& \frac{B_{0z}-B_{2z}}{\mathrm{\Δ}x}\\ \frac{B_{0y}-B_{2y}}{\mathrm{\Δ}x}& -\left(\frac{B_{0x}-B_{2x}}{\mathrm{\Δ}x}+\frac{B_{1z}-B_{3z}}{\mathrm{\Δ}z}\right)& \frac{B_{1y}-B_{3y}}{\mathrm{\Δ}z}\\ \frac{B_{1x}-B_{3x}}{\mathrm{\Δ}z}& \frac{B_{1y}-B_{3y}}{\mathrm{\Δ}z}& \frac{B_{1z}-B_{3z}}{\mathrm{\Δ}z}\end{array}\right)---\left(10\right)$

张量分量表示方式为：

$g_{uv}=\frac{{\mathrm{\ΔB}}_u}{{\mathrm{\Δs}}_v}=\frac{b_{nu}-b_{mu}}{{\mathrm{\Δs}}_v}=\frac{\left(P_{nu0}{T}_n+\underset{j=1}{\overset{h}{\mathrm{\Σ}}}\frac{d^j{Q}_{nuj}{T}_n}{{\mathrm{dt}}^j}+O_{nu}\right)-\left(P_{mu0}{T}_m+\underset{j=1}{\overset{h}{\mathrm{\Σ}}}\frac{d^j{Q}_{muj}{T}_m}{{\mathrm{dt}}^j}+O_{mu}\right)}{{\mathrm{\Δs}}_v}---\left(11\right)$

式中：m,n表示磁通门标号，取1,2,3,4；u,v都代表坐标轴的方向，取x,y,z；g_uv为惯性导航坐标系下某个张量分量；ΔB_u为惯性导航坐标系下沿u轴方向的磁场变化量；Δs_v表示磁通门n与磁通门m之间的基线距离。

将所述张量分量误差模型改写为：

$g_{uv}=P_{nu0}\frac{T_n}{{\mathrm{\Δs}}_v}-P_{mu0}\frac{T_m}{{\mathrm{\Δs}}_v}+\underset{j=1}{\overset{h}{\mathrm{\Σ}}}\left(\frac{d^j{Q}_{nuj}\frac{T_n}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}-\frac{d^j{Q}_{muj}\frac{T_m}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}\right)+O_{uv}---\left(12\right)$

式中：O_uv为张量分量g_uv的总偏置误差；

步骤三、在高空中磁场均匀的区域采集校正数据：

首先，根据地磁场强度、磁偏角和磁倾角获得地磁场三个分量，把地磁场作为四个磁通门在大地坐标系下的真实值。其次，将装载张量仪的吊舱的姿态进行随机变化来模仿高空飞行过程中张量仪的姿态变化，计算四个磁通门相应的空间位置及姿态，结合地磁场值计算四个磁通门在磁通门坐标系下的真实值。然后，按照公式(1)及四个磁通门的误差系数获得每个磁通门测量值。误差系数如表1所示，其中动态特性由一阶导数的系数控制，这是由于实际磁通门可以等效为一个一阶***。最后，用差分近似微分的方法，可以获得高空中带误差的张量数据，并以此作为高空中采集的所需校正数据。

表1仿真的误差系数

步骤四、因为高空中地球磁场的梯度近似为0，设张量分量g_uv误差模型中磁通门n中沿u轴方向测量值前面的校正参数为p_uv，将上述张量分量的误差校正模型改写为：

$\frac{g_{uv}}{p_{uv}}=\frac{P_{nu0}}{p_{uv}}\frac{T_n}{{\mathrm{\Δs}}_v}-\frac{P_{mu0}}{p_{uv}}\frac{T_m}{{\mathrm{\Δs}}_v}+\underset{j=1}{\overset{h}{\mathrm{\Σ}}}\left(\frac{d^j\frac{Q_{nuj}}{p_{uv}}\frac{T_n}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}-\frac{d^j\frac{Q_{muj}}{p_{uv}}\frac{T_m}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}\right)+\frac{O_{uv}}{p_{uv}}=0---\left(13\right)$

故式(13)中有一个测量值的校正参数为1，这样以此项测量值为已知值，对式(13)列写线性方程组，通过最小二乘拟合方法解算线性方程组，可以获得式(13)中的各项参数。此时解算出各个磁梯度张量分量误差模型中参数的相对量， $\frac{P_{\mathrm{nu}0}}{p_{\mathrm{uv}}},\frac{P_{\mathrm{mu}0}}{p_{\mathrm{uv}}},\frac{Q_{\mathrm{nuj}}}{p_{\mathrm{uv}}},\frac{Q_{\mathrm{muj}}}{p_{\mathrm{uv}}},\frac{O_{\mathrm{uv}}}{p_{\mathrm{uv}}};$

步骤五、在地面磁场不均匀的区域采集校正数据：

在地面上有磁梯度的区域，把张量仪的传感器部分绕其中心点进行旋转，记录此过程的测量值。地面张量数据模拟方法和高空类似，只有两点区别，一个是把磁偶极子产生磁场(公式17所示)叠加到地磁场上，模仿地面具有磁场梯度的情况。另一个是张量旋转不变量是通过随机改变旋转角后，再利用张量缩并公式计算而得。(此处为了和前面的公式编号保持一致，故编号未做更改)

$B(r,M)=\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}\frac{3(M\·r)r-{\mathrm{Mr}}^2}{r^5}---\left(17\right)$

式中，r为从磁源到某个位置的方位矢量，M为磁矩矢量。

步骤六、以不变量进行约束，求解张量分量的标度因子误差：

将高空校正之后所获得的各项校正系数重新带入到公式(13)中，可以获得各个张量分量的相对量值，表达为：

把张量分量的校正模型重新表达为：

则张量不变量校正模型表达为： $C_T^2=\underset{u,v=(x,y,z)}{\Σ}{g_{uv}}^2=\underset{u,v=(x,y,z)}{\Σ}{p_{uv}}^2{g_{{\mathrm{uv}}_{mid}}}^2---\left(16\right)$

将公式(16)中的p_uv都取为1，计算出此时的C_T值并进行平均处理。最后，C_T的平均值和作为公式(16)中已知量，采用广义最小二乘拟合可估计出p_uv。如表2所示，通过实际和求出的张量分量标度因子误差系数对比，两者的相对误差从0.137％到0.528％，验证不变量校正方法的正确性。

表2实际和求出的张量分量标度因子误差系数对比

	pxx	pyx	pzx	pxz	pyz	pzz
							理论值	0.9336	1.0201	0.9248	1.0522	0.9888	1.0976
解算值	0.9291	1.0187	0.9215	1.0551	0.9934	1.0918
							相对误差	0.482％	0.137％	0.357％	0.276％	0.465％	0.528％

步骤七、结合高空校正获得的各项系数，可以求得每个张量分量的最终校正参数：P_nu0,P_mu0,Q_nuj,Q_muj,O_uv；

步骤八、将三轴磁通门磁场测量值和校正参数输入到磁梯度张量分量校正模型中，计算校正后的磁梯度张量。为模拟校正前后的数据，首先按航空磁测的方式以地下10米的一个磁偶极子作为探测目标，在目标正上方40米高度，形成一条400米的测线。然后利用公式(18)计算出磁偶极子在测线上形成的理论张量。测量的张量数据由类似高空数据模拟方法生成，只是叠加上公式(17)计算的磁偶极子的磁场。

$G_{ij}=-3\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}\frac{\left(M\·r\right)\left(5r_i{r}_j-r^2{\mathrm{\δ}}_{ij}\right)-r^2\left(r_i{M}_j+r_j{M}_i\right)}{r^7}---\left(18\right)$

式中，δ_ij为Kronecker’s delta，i，j＝1,2,3表示在笛卡尔坐标系中x,y,z。

本发明方法的具体实施步骤如图1所示，实施结果如图3所示为测线上的理想、未校正和校正后张量分量之一g_zz的对比。

Claims (1)

1.一种航空器挂载式磁通门磁梯度张量仪的校正方法，其特征在于，包括以下步骤：

步骤一、建立单个磁通门误差模型：

把惯性导航***和磁通门张量仪进行刚体连接，并把各个磁通门传感器坐标系x,y,z轴分别向惯性导航坐标系的x,y,z轴方向对齐，则单个磁通门误差模型为：

$D_j\frac{d^{j}T}{{\mathrm{dt}}^j}+D_{j-1}\frac{d^{j-1}T}{{\mathrm{dt}}^{j-1}}+...+D_1\frac{dT}{dt}+D_{0}T={\mathrm{FEM}}^{-1}B+R+\mathrm{\Φ}---\left(1\right)$

式中：F为标度因子误差矩阵；E为非正交误差矩阵；R为零位误差矢量；M为磁通门坐标系与惯性导航坐标系非对准误差矩阵；D_j为磁通门测量值对时间的j阶导数的描述，即动态特性描述矩阵；T为磁通门在自身测量坐标系下的实测值；B为磁通门在惯导坐标系下的真实磁场值；Φ为磁通门的测量噪声；

每个误差矩阵或矢量的具体表达式如下：

(1)标度因子误差矩阵：

$F=\left[\begin{array}{ccc}1+S_x& 0& 0\\ 0& 1+S_y& 0\\ 0& 0& 1+S_z\end{array}\right]---\left(2\right)$

其中S_x,S_y,S_z表示磁通门感应轴的标度因子误差，

(2)非正交误差矩阵：

假设实际磁通门坐标系x轴与理想磁通门坐标系x轴的方向一致，ρ为实际磁通门坐标系y轴与理想磁通门坐标系y轴的夹角；为实际磁通门坐标系z轴与理想磁通门坐标系z轴的夹角在磁通门xoz平面上的投影；λ为实际磁通门坐标系z轴与理想磁通门坐标系z轴的夹角在磁通门yoz平面上的投影，

(3)非对准误差矩阵：

$M=\left[\begin{array}{ccc}{c}_{\mathrm{\β}}{c}_{\mathrm{\γ}}& -c_{\mathrm{\α}}{s}_{\mathrm{\γ}}+s_{\mathrm{\α}}{s}_{\mathrm{\β}}{c}_{\mathrm{\γ}}& s_{\mathrm{\α}}{s}_{\mathrm{\γ}}+c_{\mathrm{\α}}{s}_{\mathrm{\β}}{c}_{\mathrm{\γ}}\\ c_{\mathrm{\β}}{s}_{\mathrm{\γ}}& c_{\mathrm{\α}}{c}_{\mathrm{\γ}}+s_{\mathrm{\α}}{s}_{\mathrm{\β}}{s}_{\mathrm{\γ}}& -s_{\mathrm{\α}}{c}_{\mathrm{\γ}}+c_{\mathrm{\α}}{s}_{\mathrm{\β}}{s}_{\mathrm{\γ}}\\ -s_{\mathrm{\β}}& s_{\mathrm{\α}}{c}_{\mathrm{\β}}& c_{\mathrm{\α}}{c}_{\mathrm{\β}}\end{array}\right]---\left(4\right)$

式中c_α＝cosα，c_β＝cosβ，c_γ＝cosγ，s_α＝sinα，s_β＝sinβ，s_γ＝sinγ,α、β和γ分别为实际磁通门坐标系相对惯导坐标系，绕实际磁通门坐标系x轴,y轴和z轴的非对准角。当非对准角非常小的时候可以将非对准误差矩阵近似表示为：

$M=\left[\begin{array}{ccc}1& -\mathrm{\γ}& \mathrm{\β}\\ \mathrm{\γ}& 1& -\mathrm{\α}\\ -\mathrm{\β}& \mathrm{\α}& 1\end{array}\right]---\left(5\right)$

(4)零位误差矢量：

$R=\left[\begin{array}{c}{\mathrm{\μ}}_x\\ {\mathrm{\μ}}_y\\ {\mathrm{\μ}}_z\end{array}\right]---\left(6\right)$

式中μ_x,μ_y,μ_z分别为磁通门三个轴分别存在的零位误差，

(5)动态特性描述矩阵：

$D_j=\left[\begin{array}{ccc}{d}_{11}& d_{12}& d_{13}\\ d_{21}& d_{22}& d_{23}\\ d_{31}& d_{32}& d_{33}\end{array}\right]---\left(7\right)$

当磁通门由一阶***描述时取j＝1，D₁为一阶导数前面的系数矩阵。

将上述三轴磁通门的误差模型改写为校正模型：

$B=\left[\begin{array}{c}{b}_x\\ b_y\\ b_z\end{array}\right]=C_{0}T+\underset{j=1}{\overset{h}{\Σ}}\frac{d^j{C}_{j}T}{{\mathrm{dt}}^j}+O---\left(8\right)$

式中：C₀为T第0阶导数的校正参数矩阵；C_j为Τ第j阶导数的校正参数矩阵；O为零位误差校正矢量；b_x,b_y,b_z分别表示磁通门在惯性导航坐标系x,y,z轴上的磁场值，则磁通门中单个磁场分量误差校正模型可表示为：

$b_u=P_{u0}T+\underset{j=1}{\overset{h}{\Σ}}\frac{d^j{Q}_{uj}T}{{\mathrm{dt}}^j}+O_u---\left(9\right)$

式中：u表示x,y,z中一个轴，P_u0为矩阵C₀中对应b_u的行向量校正参数；Q_uj为矩阵C_j中对应b_u的行向量校正参数；O_u为总偏差矢量O_n中对应b_u的校正参数；

步骤二、建立磁梯度张量分量误差校正模型：

磁场沿x,y,z三个方向的二阶导数可构成磁梯度张量G，尽管G有九个分量，由于在不包括场源的域内，磁场的散度和旋度为0，即对称，且迹为0，故只有五个分量是独立的，表示为：

$G=\left(\begin{array}{ccc}\frac{\∂B_x}{\∂x}& \frac{\∂B_y}{\∂x}& \frac{\∂B_z}{\∂x}\\ \frac{\∂B_x}{\∂y}& \frac{\∂B_y}{\∂y}& \frac{\∂B_z}{\∂y}\\ \frac{\∂B_x}{\∂z}& \frac{\∂B_y}{\∂z}& \frac{\∂B_z}{\∂z}\end{array}\right)=\left(\begin{array}{ccc}{g}_{xx}& g_{yx}& g_{zx}\\ g_{yx}& -\left(g_{xx}+g_{zz}\right)& g_{yz}\\ g_{xz}& g_{yz}& g_{zz}\end{array}\right)---\left(10\right)$

磁通门张量仪由不同方式构建，主要有十字、四方体和三角形等结构，但都是利用差分近似微分的方式来表示磁梯度张量G的各个分量，以十字型结构为例，在惯性导航坐标系下，张量分量表示方式为：

$g_{uv}=\frac{{\mathrm{\ΔB}}_u}{{\mathrm{\Δs}}_v}=\frac{b_{nu}-b_{mu}}{{\mathrm{\Δs}}_v}=\frac{(P_{nu0}{T}_n+\underset{j=1}{\overset{h}{\Σ}}\frac{d^j{Q}_{nuj}{T}_n}{{\mathrm{dt}}^j}+O_{nu})-(P_{mu0}{T}_m+\underset{j=1}{\overset{h}{\Σ}}\frac{d^j{Q}_{muj}{T}_m}{{\mathrm{dt}}^j}+O_{mu})}{{\mathrm{\Δs}}_v}---\left(11\right)$

式中：m,n表示磁通门标号，取1,2,3,4；u,v都代表坐标轴的方向，取x,y,z；g_uv为惯性导航坐标系下某个张量分量；ΔB_u为惯性导航坐标系下沿u轴方向的磁场变化量；Δs_v表示磁通门n与磁通门m之间的基线距离。

将所述张量分量误差模型改写为：

$g_{uv}=P_{nu0}\frac{T_n}{{\mathrm{\Δs}}_v}-P_{mu0}\frac{T_m}{{\mathrm{\Δs}}_v}+\underset{j=1}{\overset{h}{\Σ}}\left(\frac{d^j{Q}_{nuj}\frac{T_n}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}-\frac{d^j{Q}_{muj}\frac{T_m}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}\right)+O_{uv}---\left(12\right)$

式中：O_uv为张量分量g_uv的总偏置误差；

步骤三、在高空中磁场均匀的区域采集校正数据：

将磁通门磁梯度张量仪挂载到航空器如直升机下方较远处。航空器飞行到高空，然后以方形，八字形等多种轨迹进行正常飞行，张量仪采集不同姿态下带误差的磁场测量数据；

步骤四、因为高空中地球磁场的梯度近似为0，设张量分量g_uv误差模型中磁通门n中沿u轴方向测量值前面的校正参数为p_uv，将上述张量分量的误差校正模型改写为：

$\frac{g_{uv}}{p_{uv}}=\frac{P_{nu0}}{p_{uv}}\frac{T_n}{{\mathrm{\Δs}}_v}-\frac{P_{mu0}}{p_{uv}}\frac{T_m}{{\mathrm{\Δs}}_v}+\underset{j=1}{\overset{h}{\Σ}}\left(\frac{d^j\frac{Q_{nuj}}{p_{uv}}\frac{T_n}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}-\frac{d^j\frac{Q_{muj}}{p_{uv}}\frac{T_m}{{\mathrm{\Δs}}_v}}{{\mathrm{dt}}^j}\right)+\frac{O_{uv}}{p_{uv}}=0---\left(13\right)$

故式(13)中有一个测量值的校正参数为1，这样以此项测量值为已知值，对式(13)列写线性方程组，通过最小二乘拟合方法解算线性方程组，可以获得式(13)中的各项参数。此时解算出各个磁梯度张量分量误差模型中参数的相对量， $\frac{P_{nu0}}{p_{uv}},\frac{P_{mu0}}{p_{uv}},\frac{Q_{nuj}}{p_{uv}},\frac{Q_{muj}}{p_{uv}},\frac{O_{uv}}{p_{uv}};$

步骤五、在地面磁场不均匀的区域采集校正数据：

在地面上有磁梯度的区域，把张量仪的传感器部分绕其中心点进行旋转，记录此过程的测量值；

步骤六、以不变量进行约束，求解张量分量的标度因子误差：

将高空校正之后所获得的各项校正系数重新带入到公式(13)中，可以获得各个张量分量的相对量值，表达为：

由于式(1)中误差矩阵的数值并不相同，因此每个分量所对应的p_uv也各不相同，故此时的张量分量可以认为存在张量分量标度因子误差。张量的不变量之一C_T，可以由即张量缩并运算获得，其特性是在坐标旋转时，一直保持值不变，当张量分量存在标度因子误差时，旋转时所获得的C_T会形成明显的波动，因此旋转过程中C_T可以用来校正张量分量的标度因子误差，

把张量分量的校正模型重新表达为：

则张量不变量校正模型表达为： $C_T^2=\underset{u,v=(x,y,z)}{\Σ}{g_{uv}}^2=\underset{u,v=(x,y,z)}{\Σ}{p_{uv}}^2{g_{{\mathrm{uv}}_{\mathrm{mi}d}}}^2---\left(16\right)$

将公式(16)中的p_uv都取为1，计算出此时的C_T值并进行平均处理。最后，C_T的平均值和作为公式(16)中已知量，采用广义最小二乘拟合可估计出p_uv；

步骤七、结合高空校正获得的各项系数，可以求得每个张量分量的最终校正参数：P_nu0,P_mu0,Q_nuj,Q_muj,O_uv；

步骤八、将三轴磁通门磁场测量值和校正参数输入到磁梯度张量分量校正模型中，计算校正后的磁梯度张量。