CN106997035A - 一种基于磁梯度不变量的磁梯度计校正方法 - Google Patents

Abstract

本发明涉及一种基于磁梯度不变量的磁梯度计校正方法，首先对构建张量仪时所存在的标度因子误差、非正交误差、非对准误差、零偏误差进行建模。然后借鉴磁矢量仪校正时的标量约束法，利用磁梯度张量两个独立的旋转不变量进行约束，获得最优的校正参数。直接从张量仪整体误差出发，对张量仪的误差进行建模，利用磁梯度张量两个独立的旋转不变量进行约束，获得最优的校正参数，更深入地对张量仪误差进行校正，使测量的数据更为可靠；相对于目前校正张量仪整体误差繁琐复杂的方法，直接对测量数据结果进行处理，不需要提供其它的数据，更为简便，此外本发明不受限于各种磁梯度测量原理既可在有地磁场的野外环境下使用，也可在磁屏蔽环境中使用。

Description

一种基于磁梯度不变量的磁梯度计校正方法

技术领域：

本发明涉及一种地球物理勘探中的磁测仪器的校正，尤其是基于磁梯度不变量的磁梯度计校正方法。

背景技术：

在地球磁场下进行磁场相关的科学研究和生产实践，已经从早期的磁场总场和磁场矢量测量发展到现在的磁梯度张量测量。磁梯度计是用来测量磁场空间梯度的装置[Pedersen L B,Rasmussen T M.The gradient tensor of potential field anomalies:Some implications on data collection and data processing of maps.Geophysics,1990,55(12):1558-1566.]。磁梯度张量仪具有以下优点：没有磁场测量中对环境磁场方向极度敏感的问题[Heath P,Heinson G,Greenhalgh S.Some comments on potentialfield tensor data.Exploration Geophysics,2003,34(2):57-62.]，具有更高的目标分辨能力和更加丰富的信息[Schmidt P W,Clark D A.The magnetic gradient tensor:Itsproperties and uses in source characterization.The Leading Edge,2006,25(1):75-78.]。因此被用于地面或水下未爆炸物检查[Pei Y H,Yeo H G,Kang X Y,etal.Magnetic gradiometer on an AUV for buried object detection,OCEANS2010.IEEE,2010:1-8.]、移动磁性目标探测[Nara T,Suzuki S,Ando S.A closed-formformula for magnetic dipole localization by measurement of its magnetic fieldand spatial gradients.IEEE Transactions on Magnetics,2006,42(10):3291-3293.]、航空或地面的矿产探查[Clark D A,Schmidt P W,Coward D A,et al.Remotedetermination of magnetic properties and improved drill targeting of magneticanomaly sources by Differential Vector Magnetometry(DVM).ExplorationGeophysics,1998,29(3/4):312-319.]，火山活动监测[Bracken R E,Grover T P,PuniwaiG S.Development and testing of a tensor magnetic gradiometer system withtrial monitoring near the Kilauea Volcano,Hawaii.US Dept.of the Interior,USGeological Survey,1998.]，地磁图的修正[Schiffler M,Queitsch M,Schneider M,etal.Calculation and Analysis of magnetic gradient tensor components of globalmagnetic models.AGU Fall Meeting Abstracts.2014,1:03.]，非屏蔽环境中磁纳米粒子的检查[Elrefai A L,Sasada I,Yoshida T.Fluxgate gradiometer for magneticnanoparticle magnetorelaxometry in unshielded environment,MagneticsConference(INTERMAG).2015IEEE.IEEE,2015:1-1.]，消除探测引力波的飞船的磁干扰[Griffin D K,Masseglia O,Hall M,et al.Design and calibration of a compactlow-noise magnetic gradiometer.Aerospace EMC,2012Proceedings ESA Workshopon.IEEE,2012:1-6.]，磁环境评估[Voigt J,Knappe-Grüneberg S,Gutkelch D,etal.Development of a vector-tensor system to measure the absolute magneticflux density and its gradient in magnetically shielded rooms.Review ofScientific Instruments,2015,86(5):055109.]等。

磁梯度张量仪按测量原理主要可以分为三种：差分近似法、旋转调制法[TilbrookD L.Rotating magnetic tensor gradiometry and a superconductingimplementation.Superconductor Science and Technology,2009,22(7):075002.]、直线测量法[Sunderland A,Ju L,Blair D G,et al.Optimizing a direct string magneticgradiometer for geophysical exploration.Review of Scientific Instruments,2009,80(10):104705.]。其中差分近似法的测量原理最为简单，且可以测量全部的张量元素，因此成为磁梯度张量的主流测量方法。由于需要矢量敏感的传感器，因此主要有基于低温或高温超导量子干涉器件的张量仪[Stolz R,Zakosarenko V,Schulz M,etal.Magnetic full-tensor SQUID gradiometer system for geophysicalapplications.The Leading Edge,2006,25(2):178-180.]和以磁通门为核心张量仪[SuiY,Li G,Wang S,et al.Compact fluxgate magnetic full-tensor gradiometer withspherical feedback coil.Review of Scientific Instruments,2014,85(1):014701.]。前者在传感器级别已经形成梯度值，后者需要先测量磁场后再做差分。

由于地球磁场的存在，目前适合野外应用的磁梯度张量仪校正方法都是利用地球磁场来获取校正参数。其基本思想是如果磁梯度张量仪存在误差，那么在磁场均匀的区域，磁梯度张量仪的输出G_ik必然是磁场的函数[Schiffler M,Queitsch M,Stolz R,etal.Calibration of SQUID vector magnetometers in full tensor gradiometrysystems.Geophysical Journal International,2014,198(2):954-964.]。这样把磁梯度张量仪在均匀磁场中旋转，并且同时测量旋转时的磁场值，就可以利用最小二乘等方法获得校正系数。这类方法有两个主要的问题：

第一，在磁场均匀区域内，单个梯度分量为零是磁梯度张量仪校正方法中的必要条件，不能保证各个梯度分量形成张量后作为整体的正确性。特别是对那些能够直接测量梯度的张量仪，此方法仍停留在校正传感器梯度不平衡度的阶段上。

第二，作为参考的磁场矢量仪由于存在标度因子误差、非正交误差、安装误差和零偏等，同样需要严格地校正。现有的磁矢量校正往往又需要磁标量作为参考，这一方面造成校正过程复杂，参与设备众多，另一方面使得张量仪的性能受辅助测量仪器性能的制约。

目前总体上来说，磁梯度张量仪的校正可以划分为两个阶段：第一个阶段是校正梯度张量中的每个分量；第二个阶段是对梯度张量仪整体存在的误差进行校正。不同测量方法获得的磁梯度分量所存在的误差各不相同，难以统一表达。现有的技术大部分还停留在第一阶段，第二阶段的校正还没有什么比较简便的方法。

发明内容：

本发明的目的就在于针对现有技术的不足，提供一种直接从第二阶段的磁梯度张量分量出发获得最优的校正参数的一种基于磁梯度不变量的磁梯度计校正方法。

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

首先对构建张量仪时所存在的标度因子误差、非正交误差、非对准误差、零偏误差进行建模。然后借鉴磁矢量仪校正时的标量约束法，利用磁梯度张量两个独立的旋转不变量进行约束，获得最优的校正参数。

一种基于磁梯度不变量的磁梯度计校正方法，包括以下步骤：

步骤一、在二维空间对张量仪误差进行建模；

步骤二、由二维空间扩展到三维空间对张量仪误差进行建模；

步骤三、建立校正补偿模型；

步骤四、利用磁偶极子先进行旋转校正数据；

步骤五、利用磁梯度不变量特性通过最小二乘法获取校正系数；

步骤六、进行测线数据测量，将张量仪装配到飞行器上，让飞行器演设计好的测线进行飞行，张量仪采集测线上的磁场数据；

步骤七、进行校正，利用步骤五中的求得的校正系数对测线上测量的数据进行校正，得到校正值。

有益效果：本发明是直接从张量仪整体误差出发，对张量仪的误差进行建模，利用磁梯度张量两个独立的旋转不变量进行约束，获得最优的校正参数，更深入地对张量仪误差进行校正，使测量的数据更为可靠；相对于目前校正张量仪整体误差繁琐复杂的方法，本发明直接对测量数据结果进行处理，不需要提供其它的数据，更为简便，此外本发明不受限于各种磁梯度测量原理既可在有地磁场的野外环境下使用，也可在磁屏蔽环境中使用。

附图说明：

图1基于磁梯度不变量的磁梯度计校正方法流程图。

图2二维空间中测量张量分量G_xx和G_yy时梯度计的理想情况。

图3二维空间中测量张量分量G_xx和G_yy时张量仪存在非对准误差。

图4二维空间中测量张量分量G_xx时张量仪存在非正交误差的情况。

图5三维空间中非正交误差示意图。

图6测量值、校正值与理论值误差的RMSE对比。

图7测线CT测量值、校正值与理论值对比。

具体实施方式：

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

一种基于磁梯度不变量的磁梯度计校正方法，包括以下步骤：

步骤一、在二维空间对张量仪误差进行建模；

步骤二、由二维空间扩展到三维空间对张量仪误差进行建模；

步骤三、建立校正补偿模型；

步骤四、利用磁偶极子先进行旋转校正数据；

步骤五、利用磁梯度不变量特性通过最小二乘法获取校正系数；

步骤六、进行测线数据测量，将张量仪装配到飞行器上，让飞行器演设计好的测线进行飞行，张量仪采集测线上的磁场数据；

步骤七、进行校正，利用步骤五中的求得的校正系数对测线上测量的数据进行校正，得到校正值。

步骤一、在二维空间对张量仪误差进行建模

实际的仪器制造和安装过程中存在在非正交误差与非对准误差。在二维空间中对两个理论上测量磁场X分量在x轴梯度G_xx与磁场Y分量在y轴方向上梯度g_yy的梯度计进行分析，得到G_xx分量测量值R_xx与G_yy分量测量值R_yy的表示方式：先通过公式(1)和公式(2)分别计算矩阵Nx与Ny,然后取其中的N_x11为R_xx，N_y22为R_yy，其中，α为基线非对准角，γ为基线非正交角，β为梯度计所测磁场分量非对准角，θ为梯度计所测磁场分量非正交角。同理，二维空间中其他梯度的误差模型也可以通过类似方式表示。其中c*表示cos*，s*表示sin*。

步骤二、由二维空间扩展到三维空间对张量仪误差进行建模，并且建立补偿校正模型

假设R_xx和R_yy是G_xx和G_yy分量的测量值，其可通过公式(3)-(5)计算，分别取N_x和N_y中的N_x11和N_y22为R_xx and R_yy。

其中：

α_i,β_i,γ_i和δ_i表示误差角。

因此，每个测量的分量R_μυ(μ,υ＝x,y,和z)只含有非正交与非对准误差，可以表示为：

步骤三、可以建立校正补偿模型

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

而实际的测量值T_μυ还含有标度因子误差c,和零偏误差o,可以表示为：

T_μυ＝c_μυN_μυ+o_μυ. (7)

因此只含有五个独立分量的张量测量值可以表示为：

其中张量的校正模型可以写为:

步骤四、利用磁偶极子先进行旋转校正数据获取

在地面将仪器装配好，以磁偶极子作为磁源，让张量仪绕轴进行旋转，旋转要尽量覆盖到各个角度，张量仪采集不同姿态下带磁干扰的磁场测量数据。

步骤五、利用磁梯度不变量特性通过最小二乘法获取校正系数

不变量是指在张量进行坐标变换时其值可以通过一系列计算得到一个不变常数，张量具有两个非零独立不变量：

I₂＝|G|. (11)

第一步是用不变量I₁，进行校正：将公式(12)中的G_μυ用公式(9)中的T_μυ代入展开可以得到公式(12)

其中x₁-x₁₅是h_ij的多项展开式，x₁₆-x₂₀由h_ij和d_ij组成。

当磁张量梯度仪绕其中心旋转，它可以测量n组带误差的张量值。根据不变的性质，x₁-x₂₀应该保证F_i(i＝1...n)不变.所以,我们设计公式(15)中的函数J利用最小二乘法计算出最佳的x₁-x₂₀。

将(9),(10),和(12)联立可以得到二次型：

(T-O)^TA(T-O)＝1. (14)

其中

所以零偏误差O可以通过x₁-x₂₀来计算。将零偏误差消除之后，(9)和(10)可以修改为：

同样的，y₁-y₁₅是h_ij组成的多项式可以通过最小二乘法进行计算。

以系数y₁为例进行展开有:

所以,15个系数可以得到15个方程，将这15个方程联立可以得到方程组E₁。由于方程个数小于未知量个数，因此还需要使用不变量I₂进行校正。

将(11)展开,我们得到了一个关于G_μυ的多项式：

将公式(18)中的用(15)中的T'代入可以得到公式(19)：

其中，z₁-z₃₅是h_ij组成的多项式可以通过最小二乘法进行计算。

因此35个等式可以组成方程组E₂，将E₁和E₂联立通过MATLAB的数值计算可以解出方程组的解得到校正系数，最后我们可以通过公式(9)来进行校正。

步骤六、进行测线数据测量

将张量仪装配到飞行器上，让飞行器演设计好的测线进行飞行，张量仪采集测线上的磁场数据。

步骤七、进行校正

利用步骤五中的求得的校正系数对测线上测量的数据进行校正，得到校正值。

本方法以单个磁梯度计为基础，由二维空间先对其构建张量仪时存在的非对准误差、非正交误差进行建模，再扩展到三维空间，在建模中加入标度因子误差、零偏误差建立校正补偿模型。现在地面进行旋转校正数据获取，利用磁梯度张量两个独立的旋转不变量进行约束，利用最小二乘法获得最优的校正参数，再将校正系数应用于测线数据校正，获取校正数据。以测量g_xx，g_xy，g_xz，g_yy和g_yz分量的张量仪进行分析，提供一种张量仪的校正补偿方法。

在具体应用实例中，可以通过仿真进行模拟。该方法通过以下步骤实现：

步骤一、在二维空间对张量仪误差进行建模

先在二维空间中分析非正交与非对准误差的情况，图2所示为测量磁场X分量在x轴方向与Y分量y轴方向上的梯度(G_xx,G_yy)时,两个梯度计的摆放示意图，其中虚线箭头表示梯度计的感应方向。

在制造和安装过程中存在非对准和非正交误差。图3示出了分别在基线和感应方向上具有非对准角度α和β的梯度计。原本应该测量G_yy分量的梯度计实际上测量的是y'轴方向上的磁场分量Y'的变化。Y'的感应方向等同于Y以角α和β为旋转角旋转两次的感测方向。

此外，x″轴和y'轴之间的非正交角γ以及感应方向X″和Y'之间的非正交角θ也在图4中给出，原本应该测量G_xx分量的梯度计实际测量值为磁场分量X″在x”轴方向上的变化。

二维空间中理论张量为：

因此如果实际测量值G_yy有图1中所示的误差，实际上Y'的感应方向等同于Y以角α和β为旋转角旋转两次，而y'轴相当于以角α旋转了一次。所以G_yy测量值可以通过公式(2)计算，用R_yy表示。R_yy值为N_y中的N_y22：

测量的分量G_xx的梯度计的基线存在非正交角γ。由于基线中的非正交误差，感测方向已经以角度γ进行旋转。此外，在X”和Y'的感测方向之间存在另一非正交角度θ。假设R_xx是G_xx的测量值，其可以经由(3)计算，R_xx值为N_x中的N_x11。

其中s_*和c_*表示正弦和余弦函数。二维空间中其他分量的表达方法与之类似。

步骤二、由二维空间扩展到三维空间对张量仪误差进行建模，并且建立补偿校正模型

我们以G_xx和G_yy分量为例推出三维空间对张量仪误差模型的表示。具体如下：

⑴、假设在基线中存在非对准角α₁,β₁,和γ₁，则梯度计可以被认为以这三个角为欧拉角绕轴旋转，变换矩阵可以由M定义。同样的在感应方向中若存在非对准角α₂,β₂,和γ₂，则其变换矩阵可以由L定义。

⑵、假设在基线中存在三个非正交角α₃,β₃,和γ₃，新产生的坐标轴为x',y',和z。其中理论上的z轴与实际z轴是重合的；y'oz平面与yoz也是重合的，δ₃是y'轴与y轴之间的角度；p是x'轴在xoy平面上的投影，其中α₃是x'轴和p轴之间的夹角，β₃是x轴和p轴之间的夹角。所以非正交变换矩阵可以由W₃定义。同时感应方向受到了基线非正交的影响产生变化：X'相当于X方向分别绕z轴和y轴以β₃角和α₃角依次旋转两次得到，变换矩阵可以由W₂定义；Y'相当于Y方向绕x轴以δ₃角旋转得到，变换矩阵可以由W₅定义。

⑶、假设在感应方向中存在三个非正交角α₄,β₄,和γ₄，其变换可以由相类似的方法推导，非正交矩阵由W₁和W₄定义。

假设R_xx和R_yy是G_xx和G_yy分量的测量值，其可通过公式(4)-(5)计算，分别取N_x和N_y中的N_x11和N_y22为R_xx and R_yy。

其中：

α_i,β_i,γ_i和δ_i表示误差角。

因此，每个测量的分量R_μυ(μ,υ＝x,y,和z)只含有非正交与非对准误差，可以表示为：

步骤三、可以建立校正补偿模型

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

而实际的测量值T_μυ还含有标度因子误差c,和零偏误差o,可以表示为：

T_μυ＝c_μυN_μυ+o_μυ. (8)

因此只含有五个独立分量的张量测量值可以表示为：

其中张量的校正模型可以写为:

其中包含了30个待求系数,H中有25个,写作h_ij,D中有5个写作d_ij.

步骤四、利用磁偶极子先进行旋转校正数据获取

在地面线将仪器装配好，以磁偶极子作为磁源，让张量仪绕轴进行旋转，旋转要尽量覆盖到各个角度，张量仪采集不同姿态下带磁干扰的磁场测量数据。

在仿真中模拟有5个安装好的梯度计分别测量磁梯度张量中的五个独立的分量G_xx,G_xy,G_xz,G_yy,和G_yz。在地面先进行旋转校正，再在低空的一条测线上用旋转校正得到的校正系数观测校正效果，具体的仿真有如下设定：

①、仿真中的磁偶极子在某一位置形成的磁梯度张量可以由磁场强度、磁偏角、磁倾角以及其具体位置通过公式(11)计算。同样，在测线上计算每个采样点的梯度值便可以得到整条测线上的梯度曲线。

其中，μ₀是真空中磁导率，r是方向向量，M是磁偶极子磁矩，δ_ij是克罗内克δ，i,j＝1,2,3表示笛卡尔坐标系下的x,y,z。

②、给张量仪加入随机的旋转角，计算得到不同姿态下张量仪理论输出值，用于得到校正系数。

③、张量仪的输出可以由各种误差由上述模型导出，具体的误差设置如表1所示。

④、将磁矩矢量为(250,000,250,000,-353,553.39)A·m²的磁偶极子作为测量目标埋入地下5米，采集磁偶极子正上方距离地面10米高处的测线数据，测线长100米，测线上均匀设置了500个采样点，得到测线上的测量数据。

表1仿真中张量仪误差设定

标度因子误差	在-0.01—0.01之间均匀分布
		轴间非正交误差	在-1°—1°之间均匀分布
轴间非对准误差	在-3°—3°之间均匀分布
		感应非对准误差	在-1°—1°之间均匀分布
零偏	在-15nT/m—15nT/m之间均匀分布
		噪声	振幅为0.1nT/m的均匀分布噪声。

步骤五、利用磁梯度不变量特性通过最小二乘法获取校正系数

不变量是指在张量进行坐标变换时其值可以通过一系列计算得到一个不变常数，张量具有两个非零独立不变量：

I₂＝|G|. (13)

第一步是用不变量I₁，进行校正：将公式(12)中的G_μυ用公式(10)中的T_μυ代入展开可以得到公式(14)

其中x₁-x₁₅是h_ij的多项展开式，x₁₆-x₂₀由h_ij和d_ij组成。

当磁张量梯度仪绕其中心旋转，它可以测量n组带误差的张量值。根据不变的性质，x₁-x₂₀应该保证F_i(i＝1...n)不变.所以,我们设计公式(15)中的函数J利用最小二乘法计算出最佳的x₁-x₂₀。

将(10),(12),和(14)联立可以得到二次型：

(T-O)^TA(T-O)＝1. (16)

其中

所以零偏误差O可以通过x₁-x₂₀来计算。将零偏误差消除之后，公式(10)和公式(12)修改为：

同样的，y₁-y₁₅是h_ij组成的多项式可以通过最小二乘法进行计算。

以系数y₁为例进行展开有:

所以,15个系数可以得到15个方程，将这15个方程联立可以得到方程组E₁。由于方程个数小于未知量个数，因此还需要使用不变量I₂进行校正。

将(13)展开,我们得到了一个关于G_μυ的多项式：

将公式(18)中的用(17)中的T'代入可以得到公式(21)：

其中，z₁-z₃₅是h_ij组成的多项式可以通过最小二乘法进行计算。

因此35个等式可以组成方程组E₂，将E₁和E₂联立通过MATLAB的数值计算可以解出方程组的解得到校正系数，最后我们可以通过公式(10)来进行校正。

步骤六、进行测线数据测量

将张量仪装配到飞行器上，按表1中的设定仿真飞行器沿设计好的测线进行飞行，张量仪采集测线上的磁场数据。

步骤七、进行校正

利用步骤五中的求得的校正系数对测线上测量的数据进行校正，得到校正值。具体的校正效果如下：所有仿真系数根据表1生成，仿真重复1000次，均方根误差用来评估理论和未校正G_xy和校正G_xy之间的误差,1000次仿真中每一对均方根误差仿真结果如图6所示。其他分量校正结果结果如表2所示，校正的效果为38.58倍到60.83倍.

表2五个分量1000次仿真中平均RMSEs

	Gxx	Gxy	Gxz	Gyy	Gyz
						校正前(nT/m)	7.30	7.50	7.50	7.33	7.54
校正后(nT/m)	0.12	0.13	0.18	0.19	0.15
						改善幅度	60.83	57.70	41.67	38.58	50.27

图7展示了在某一次仿真中，在测区中间两条测线上理论CT与校正前后CT的对比结果，RMSE在校正之后从17.78nT/m减小到0.09nT/m。五个分量校正结果类似于CT，具体如表3所示:五个分量校正结果类似于CT如表3所示.校正改善效果为30.15倍到99.89倍.

表3一次仿真中五个分量的RMSE

	Gxx	Gxy	Gxz	Gyy	Gyz
						校正前(nT/m)	5.00	3.92	9.77	6.98	8.99
校正后(nT/m)	0.08	0.13	0.25	0.14	0.09
						改善幅度	62.50	30.15	39.08	49.86	99.89

Claims (1)

1.一种基于磁梯度不变量的磁梯度计校正方法，其特征在于，包括以下步骤：

步骤一、在二维空间对张量仪误差进行建模；

步骤二、由二维空间扩展到三维空间对张量仪误差进行建模；

步骤三、建立校正补偿模型；

步骤四、利用磁偶极子先进行旋转校正数据；

步骤五、利用磁梯度不变量特性通过最小二乘法获取校正系数；

步骤六、进行测线数据测量，将张量仪装配到飞行器上，让飞行器演设计好的测线进行飞行，张量仪采集测线上的磁场数据；

步骤七、进行校正，利用步骤五中的求得的校正系数对测线上测量的数据进行校正，得到校正值。