CN109856689A - 一种超导航磁梯度张量数据抑噪处理方法和系统 - Google Patents

Abstract

本发明公开一种超导航磁梯度张量数据抑噪处理方法及系统。方法包括：获取在飞行测量平台坐标系下的大地背景参考磁场值、大地背景参考磁梯度场值和飞行测量线上的背景参考磁场值的时间变化率；根据获得的数据对三分量磁强计数据进行补偿校正处理，得到磁场校正补偿处理数据；根据磁场校正补偿处理数据、大地背景参考磁场值和飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型；根据磁梯度数据补偿模型，得到补偿后的磁梯度数据；对补偿后的磁梯度数据进行低通滤波处理，得到处理后的磁梯度数据；根据处理后的磁梯度数据对补偿效果进行评估。本发明能够提高实际测量数据的补偿有效性。

Description

一种超导航磁梯度张量数据抑噪处理方法和系统

技术领域

本发明涉及地球物理航空磁法勘探领域，特别是涉及一种超导航磁梯度张量数据抑噪处理方法和系统。

背景技术

航空磁法勘探是一种重要并广泛应用于矿产和环境勘测的地球物理勘探技术手段。近年来，得益于基于超导量子干涉仪(SQUID)的磁梯度测量系统的大力发展，已使得直接测量高精度的磁梯度数据成为可能，而超导全张量磁梯度仪因其高信息量、高磁场灵敏度、体积小等诸多优点，被认为是第三代航磁探测的发展方向。目前，我国已建设有基于SQUID的全张量磁梯度测量平台原型机，不同于传统的采用机载方式进行测量的航磁梯度测量飞行平台，该原型机将磁梯度测量仪及相关设备(GPS、惯导系统(INS)、读入读出电路等)装载于一个吊舱内，并通过直升机拖曳的方式进行飞行测量，该方式将大大降低直升机对磁测量干扰的影响，甚至在拖曳缆绳足够长的条件下直升机影响可以忽略不计，但吊舱内仍然存在大量会产生磁干扰的设备，因此有效地通过磁补偿处理提高飞行平台的测量精度，以实现高分辨率磁梯度测量具有重要的意义。

目前针对该问题的方法还十分有限，除了通过硬件设备压制干扰的方法外，已有的软补偿方法多数仅考虑传统干扰源对测量的影响，但在实际飞行测量过程中往往噪声系统更为复杂，仅通过一种模型及相关的单一流程很难将复杂的各类干扰充分滤除。

发明内容

本发明的目的是提供一种超导航磁梯度张量数据抑噪处理方法和系统，能够提高实际测量数据的补偿有效性。

为实现上述目的，本发明提供了如下方案：

一种超导航磁梯度张量数据抑噪处理方法，包括：

获取在飞行测量平台坐标系下的大地背景参考磁场值、大地背景参考磁梯度场值和飞行测量线上的背景参考磁场值的时间变化率；

根据所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对三分量磁强计数据进行补偿校正处理，得到磁场校正补偿处理数据；

根据所述磁场校正补偿处理数据、所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型；

根据所述磁梯度数据补偿模型，得到补偿后的磁梯度数据；

对所述补偿后的磁梯度数据进行低通滤波处理，得到处理后的磁梯度数据；

根据所述处理后的磁梯度数据对补偿效果进行评估。

可选的，所述获取在飞行测量平台坐标系下的大地背景参考磁场值、大地背景参考磁梯度场值和飞行测量线上的背景参考磁场值的时间变化率，具体包括：

获取在大地坐标系下的大地背景参考磁场值和大地背景参考磁梯度场值；

将所述大地坐标系下的大地背景参考磁场值和大地背景参考磁梯度场值转换到飞行测量平台坐标系中，得到在飞行测量平台坐标系下的大地背景参考磁场值和大地背景参考磁梯度场值；

根据所述飞行测量平台坐标系下的大地背景参考磁场值，通过数值差分或傅里叶变换的方式，得到飞行测量线上的背景参考磁场值的时间变化率。

可选的，所述根据所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对三分量磁强计数据进行补偿校正处理，得到磁场校正补偿处理数据，具体包括：

通过正则化最小二乘拟合或Huber范数拟合的方式求解方程得到第一校正补偿系数c₁；其中，A为直接获取的观测数据；

对测量数据B_m进行第一阶段校正补偿计算，得到初始处理数据

通过正则化最小二乘拟合或Huber范数拟合的方式求解方程得到第二校正补偿系数c₂；

对数据B1_m进行第二阶段补偿计算，通过公式得到磁场校正补偿处理数据。

可选的，所述根据所述磁场校正补偿处理数据、所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型，具体包括：

根据所述磁场校正补偿处理数据、所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型

其中，l为通道编号；为通道测量的磁梯度值；为通道l的补偿磁梯度值；H_e＝(H_ex,H_ey,H_ez)为参考背景磁场，dH_e＝(dH_ex,dH_ey,dH_ez)为参考背景磁场值H_e的时间变化率，B_r＝(B_rx,B_ry,B_rz)为测量磁场校正补偿处理后数据，t为惯导数据采样时间，fs为对应采样频率，p,q为干扰源信号滤波基长度，相应的，A为由各类观测数据及相关计算量构成的N×M阶矩阵，x＝(s1,o1,a_1,x,a_2,x,K,a_2p+1,z,b_1,x,K,b_2p+1,z,c_1,x,K,c_2q+1,z)^T为实际计算中使用的补偿参数，其中s1＝1/s^l，o1＝o^l/s^l，s^l为测量通道的尺度偏差；o^l为测量通道的偏差系数。

可选的，所述根据所述磁梯度数据补偿模型，得到补偿后的磁梯度数据，具体包括：

根据所述磁梯度数据补偿模型求解补偿系数；

根据所述补偿系数，确定补偿后的磁梯度数据。

可选的，所述对所述补偿后的磁梯度数据进行低通滤波处理，得到处理后的磁梯度数据，具体包括：

对所述补偿后的磁梯度数据采用巴特沃斯低通滤波方法进行低通滤波处理，得到处理后的磁梯度数据。

可选的，所述根据所述处理后的磁梯度数据对补偿效果进行评估，具体包括：

根据所述处理后的磁梯度数据，采用标准化参数对补偿效果进行评估，所述标准化参数包括均方根、标准差和改善比。

一种超导航磁梯度张量数据抑噪处理系统，包括：

获取模块，用于获取在飞行测量平台坐标系下的大地背景参考磁场值、大地背景参考磁梯度场值和飞行测量线上的背景参考磁场值的时间变化率；

磁场校正补偿处理数据确定模块，用于根据所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对三分量磁强计数据进行补偿校正处理，得到磁场校正补偿处理数据；

磁梯度数据补偿模型确定模块，用于根据所述磁场校正补偿处理数据、所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型；

补偿磁梯度数据确定模块，用于根据所述磁梯度数据补偿模型，得到补偿后的磁梯度数据；

磁梯度数据处理模块，用于对所述补偿后的磁梯度数据进行低通滤波处理，得到处理后的磁梯度数据；

补偿效果评估模块，用于根据所述处理后的磁梯度数据对补偿效果进行评估。

根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明提供一种超导航磁梯度张量数据抑噪处理方法，获取在飞行测量平台坐标系下的大地背景参考磁场值、大地背景参考磁梯度场值和飞行测量线上的背景参考磁场值的时间变化率；根据所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对三分量磁强计数据进行补偿校正处理，得到磁场校正补偿处理数据；根据所述磁场校正补偿处理数据、所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型；根据所述磁梯度数据补偿模型，得到补偿后的磁梯度数据；对所述补偿后的磁梯度数据进行低通滤波处理，得到处理后的磁梯度数据；根据所述处理后的磁梯度数据对补偿效果进行评估。本发明能够有效的解决多种测量平台系统噪声和环境噪声的干扰，提高了SQUID航磁梯度仪的磁梯度数据补偿的精度和有效性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例超导航磁梯度张量数据抑噪处理方法流程图；

图2为在飞行测量坐标系下对三分量磁强计测量值的校正补偿结果示意图，其中，(a)表示B_x分量校正前、校正后与相应参考场数据的对比结果；(b)表示B_y分量校正前、校正后与相应参考场数据的对比结果；(c)表示B_z分量校正前、校正后与相应参考场数据的对比结果；(d)表示磁场模量|B|校正前、校正后与相应参考场数据的对比结果；

图3为磁场三分量数据补偿校正前后，各磁场分量Fourier频率幅值谱分析的对比图，其中，(a)表示B_x分量校正前、校正后与相应参考场数据的对比结果；(b)表示B_y分量校正前、校正后与相应参考场数据的对比结果；(c)表示B_z分量校正前、校正后与相应参考场数据的对比结果；

图4为在飞行测量坐标系下对G1通道的磁梯度计测量值的补偿及滤波处理结果示意图，其中，(a)表示补偿前与补偿后数据的对比结果；(b)表示进行补偿处理后，滤波前与滤波后数据的对比结果；

图5为对G1通道的磁梯度计测量值进行补偿及滤波处理后，测量数据、补偿后数据、滤波后数据的Fourier频率幅值谱分析对比图；

图6为本发明实施例超导航磁梯度张量数据抑噪处理系统结构图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种超导航磁梯度张量数据抑噪处理方法和系统，能够提高实际测量数据的补偿有效性。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明实施例超导航磁梯度张量数据抑噪处理方法流程图。如图1所示，一种超导航磁梯度张量数据抑噪处理方法，包括：

获取在飞行测量平台坐标系下的大地背景参考磁场值、大地背景参考磁梯度场值和飞行测量线上的背景参考磁场值的时间变化率；

步骤101：根据所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对三分量磁强计数据进行补偿校正处理，得到磁场校正补偿处理数据；

具体包括：

获取在大地坐标系下的大地背景参考磁场值和大地背景参考磁梯度场值；

将所述大地坐标系下的大地背景参考磁场值和大地背景参考磁梯度场值转换到飞行测量平台坐标系中，得到在飞行测量平台坐标系下的大地背景参考磁场值和大地背景参考磁梯度场值；

根据所述飞行测量平台坐标系下的大地背景参考磁场值，通过数值差分或傅里叶变换的方式，得到飞行测量线上的背景参考磁场值的时间变化率。

首先通过飞行测量平台搭载的GPS系统获取各测量点处的经纬高数据(经度，纬度，高度)、通过全球地磁场参考模型IGRF-12以及区域地磁台网观测数据，可以计算出大地坐标系下(NEU),大地背景参考磁场值H_e0＝(H_e0x,H_e0y,H_e0z)^T、大地背景参考磁梯度场值H_ge0＝(H_0xx,H_0xy,H_0xz,H_0yy,H_0yz,H_0zz)。

然后，通过INS系统获取各测量点处的测量平台的飞行姿态数据(偏转角η、俯仰角λ、翻滚角ψ),可将大地坐标系下的H_e0及H_ge0转换到飞行测量平台坐标系下。其中，转换后的背景参考磁场值H_e及参考磁梯度场值H_ge可以分别表示为：

H_e＝(H_ex,H_ey,H_ez)^T＝T₂H_e (1)

其中，T₂为大地坐标系到飞行测量坐标系的坐标转换矩阵，即：

由于在实际中受飞行测量设备的限制，使得GPS、INS的数据采样率与磁强计B_m、磁梯度计G_i的采样率不一致(本实施例中，前者采样率要远低于后者)，因此造成相应计算得到的参考磁场数据H_e及H_ge的分辨率要低于磁传感器测量信号B_m及G_i的分辨率。对此，在插值精度允许的条件下，利用三次样条插值对H_e及H_ge进行数据加密，同时对B_m及G_i进行下采样处理，使得原来两个采样率不同的数据集经过处理后获得同样长度的数据样本，以便用于下面磁校正及补偿处理。为了简便表达式，经过插值或下采样处理的数据用原有同样的符号表示。

最后，将H_e看作时间序列H_e(t)，通过数值差分或FFT变换的方式，计算在飞行测量线上的背景磁场值H_e的时间变化率dH_e＝(dH_x,dH_y,dH_z)＝(H_e(t+1/fs)-H_e(t))·fs，其中fs为数据采样频率。

步骤102：根据所述磁场校正补偿处理数据、所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型；

磁梯度张量数据测量仪是由6个不共面的一阶SQUID平面磁梯度计和3轴正交参考SQUID磁强计组成的9通道信号输出设备，由于超导芯片制备工艺的限制使得平面梯度计的两个线圈具有不一致性，因此每一个梯度计信号中还会包含寄生的三分量磁场信号。虽然多梯度计模块的几何构建可以降低自身固有的不平衡性，但该精度仍然无法满足在地磁背景场下条件的实际测量需求，因此需要在梯度计附近安装三轴磁强计，同时测量磁场三分量场，用于测量的磁梯度场数据的平衡度补偿处理。由于三轴磁强计与磁梯度计的测量环境相同，同样会受到多种干扰噪声的影响，因此在用于磁梯度计补偿前需要对其本身采集的磁场数据进行校正补偿处理。

首先建立基于飞行测量平台坐标系下的磁场三分量数据校正补偿综合模型。考虑磁测传感器的低频噪声干扰主要来源于两个部分：内部噪声(多个传感器敏感度差异性、传感器偏差、安装误差及软磁性干扰)及外部干扰(测量平台的剩磁、感磁、涡流磁场干扰)，结合三轴磁强计校正模型及Tolles-Lawson磁补偿模型，建立如下的三分量磁场数据校正补偿综合模型：

B_m＝SDB₀+O₁+O₂+M₁H_e+M₂dH_e+n (4)

其中，B₀＝(B_x,B_y,B_z)为在飞行测量坐标系下的真实磁场三分量值，在高空测量条件下应与H_e近似相等；

S＝diag(s_x,s_y,s_z)为尺度偏差因子；

为安装误差矩阵，偏差角a_xy为实际测量x’轴在直角坐标系xoy平面上的投影与x轴的小夹角，其他偏差角定义类同；

O₁＝diag(o_1x,o_1y,o_1z)为传感器偏差系数矩阵；

O₂＝diag(o_2x,o_2y,o_2z)为测量平台的剩磁干扰系数矩阵

M₁＝(m1_ij)_3×3为测量平台的感磁干扰系数矩阵；

M₂＝(m2_ij)_3×3为测量平台的涡流干扰系数矩阵；

n为传感器高斯白噪声。

因此，上述模型可以整理得到如下磁场校正补偿系数求解方程：

其中，为实际计算使用的校正补偿系数，其各个部分与原始具有物理含义的噪声干扰系数的关系为：F(B_m,1,dH_e)为由数据B_m,dH_e,1形成的方程系数矩阵，当采样数据量为N时，该矩阵维度为3N×21。

该校正补偿系数求解方程可以通过正则化最小二乘拟合或Huber范数拟合的方式进行求解，其求解精度由系数矩阵F的条件数决定。由理论模型测试分析可知，由B_m和dH_e构成的向量组成线性相关，因此在存在数据噪声的情况下直接求解该方程会造成严重的计算不稳定，因此实际计算中将校正补偿系数c分为两组及分别进行求解，并对磁场数据进行分阶段校正补偿处理，具体方法如下所述：

通过正则化最小二乘拟合或Huber范数拟合的方式求解方程得到第一校正补偿系数c₁；其中，A为直接获取的观测数据；

对测量数据B_m进行第一阶段校正补偿计算，得到初始处理数据

通过正则化最小二乘拟合或Huber范数拟合的方式求解方程得到第二校正补偿系数c₂；

对数据B1_m进行第二阶段补偿计算，通过公式得到磁场校正补偿处理数据。

步骤103：根据所述磁梯度数据补偿模型，得到补偿后的磁梯度数据；

首先，由于是对磁梯度计单一通道测量数据进行补偿，因此需要将位于飞行测量坐标系下的参考背景磁梯度场H_ge坐标旋转到与待补偿通道平行的方向，设坐标旋转矩阵为T₃，在实施例中，六个平面磁梯度计成六棱台分布，每一片在安装支架平面圆的投影成夹角为60度的均匀分布，具体的转化矩阵形式如下所示：

其中，分别为与6个平面磁梯度计通道方向平行的参考背景磁梯度场值，α为平面梯度计倾斜面与安装支架平面的夹角(该采样设备倾角值α＝63°)。

然后，建立针对单一磁梯度计通道测量数据的综合补偿模型：

根据平面磁梯度计的工艺说明及建立三分量磁场校正补偿综合模型的原理可知，磁梯度计中单一通道的低频噪声干扰主要来源于三个部分：内部偏差(多个传感器敏感度差异性、传感器偏差、安装误差及软磁性干扰)、外部干扰(测量平台的剩磁、感磁、涡流磁场)以及寄生磁场。

综合考虑，外部干扰源等价的磁矩矢量M与数据H_e,dH_e的关系可以表示为：

其中，M_r,M_ind,M_v分别表示外部干扰源等价的剩磁、感磁、涡流磁矩矢量；k＝(k_x,k_y,k_z)为感磁系数，e＝(e_x,e_y,e_z)为涡流系数，其分别与干扰源的磁性结构及几何形态分布有关。

而磁矩矢量M与磁梯度张量分量B_ij的关系可以表示为如下形式，也同样为一种线性关系：

其中，G_ij＝(G_ijx,G_ijy,G_ijz)对应于磁梯度分量的核函数向量，其具体公式为：

(1)单位球形源：

(2)正六面体源：

r_p(x_p,y_p,z_p)为测量点位置，r_Q(x_Q,y_Q,z_Q)为干扰源位置；Ω_Q为干扰源的几何形状，如将该区域干扰源离散近似为一系列等价球形偶极子源，则外部干扰源产生的磁梯度场可以近似表示为如下形式：

其中，N表示离散的磁球个数，V表示离散的磁球体积。

结合上述的飞行测量坐标系到测量通道平面的转换公式T₃(公式(6))，以及外部干扰源的磁矩模型(公式(7))，可以将外部干扰源沿单一通道方向的磁梯度测量干扰场表示为如下形式：

其中，以及g₀,g₁是由坐标转换变换T₃ ^l、磁梯度核函数G_ij、离散近似磁球体积V、以及感磁系数k，涡流系数e等常值参数构成的线性算子。经过一系列等价转换后，初始的复杂干扰场方程被转化为一个关于向量1,H_e,dH_e的简单线性组合形式，这种表示可以理解为外部干扰场空间的局部可以表示为如下一个线性空间：

V_out＝{G_out|G_out∈span{1,H_e,dH_e}}。 (13)

因此，基于类似信号线性空间理论，类比建立了磁梯度测量的寄生磁场空间：

V_balance＝{G_balance|G_balance∈span{B_r}}。 (14)

进一步，考虑数据采集中的系统误差及数据预处理过程造成的数值误差，比如不同采集设备、同一采集设备不同采集通道之间的信号延迟，以及步骤101中对于采样频率不一致的数据采用插值或下采样的处理方式，均会造成不同采样数据间无法做到完全同步。因此，如若干扰场的基信号仅选取当前时间点的一个值进行线性估计则很可能造成不可忽略的偏差量。为了更好的处理这一问题，本发明借鉴了线性滤波的思想对基本模型进行了优化：选取基信号会首先利用临近点信息对当前点的基信号进行线性估计，然后再利用估计量作为基向量对干扰场进行线性重构，此时需要将采集信号看作时间序列。据此，修正后的外部干扰场空间和寄生磁场空间可以表示为：

其中，2p+1,2q+1为两组基向量的估计域长度。综上所述，结合采集通道信号的内部偏差校正量，具体实施中采用如下模型对单一通道采集信号进行表示：

其中，l为通道编号；为通道测量的磁梯度值；G^l为通道测量点处的真实磁梯度值；s^l为测量通道的尺度偏差；为测量通道的系统偏差，为外部干扰源的剩磁干扰，两者由于在一定时间内对应单一通道测量值均为恒定不变的常值，可综合为一个偏差系数o^l；为外部干扰源感磁干扰的线性算子，由(2p+1)×3个系数组成；为外部干扰源涡流干扰的线性算子，由(2p+1)×3个系数组成；h₀为寄生磁场干扰的线性算子，由(2q+1)×3个系数组成；n为高频Gauss白噪声。通过对上述采集信号模型进行整理，可以得到用于磁梯度数据补偿的计算模型，如下所示：

其中，A为由各类观测数据及相关计算量构成的N×M阶矩阵，x＝(s1,o1,a_1,x,a_2,x,K,a_2p+1,z,b_1,x,K,b_2p+1,z,c_1,x,K,c_2q+1,z)^T为实际计算中使用的补偿参数，共M＝12p+6q+11个；需要注意的是，在实际计算中，由于补偿系数具有一定的物理含义，如s^l，因此需要根据实际传感器性能定义求解值域约束条件，如实施例中设约束区间s^l∈[0.5,1]。

根据上述磁梯度数据综合补偿模型，建立如下的约束优化问题求解补偿系数：

(1)对于测量区域内梯度值很小的情况下，采用带约束的线性最小二乘数据拟合的方法求解补偿系数x，求解算法可采用信赖域反射算法(Trust region reflective)，有效集算法(Active set)或内点法(Interiorpoint)等。实施例中采用的是常用的有效集算法求解如下问题：

(2)对于测量区域内存在局部较大梯度值的情况下，为了保证求解的稳定性，采用Huber数据拟合的方式求解补偿系数x，求解算法采用交替方向乘子方法(ADMM)：

通过如上计算方案求得补偿系数x之后，使用如下公式计算信号通道l采集的磁梯度补偿值：

其中，l为通道编号；为通道测量的磁梯度值；为通道l的补偿磁梯度值；H_e＝(H_ex,H_ey,H_ez)为参考背景磁场，dH_e＝(dH_ex,dH_ey,dH_ez)为参考背景磁场值H_e的时间变化率，B_r＝(B_rx,B_ry,B_rz)为测量磁场校正补偿处理后数据，t为惯导数据采样时间，fs为对应采样频率，p,q为干扰源信号滤波基长度，相应的，A为由各类观测数据及相关计算量构成的N×M阶矩阵，x＝(s1,o1,a_1,x,a_2,x,K,a_2p+1,z,b_1,x,K,b_2p+1,z,c_1,x,K,c_2q+1,z)^T为实际计算中使用的补偿参数，其中s1＝1/s^l，o1＝o^l/s^l，s^l为测量通道的尺度偏差；o^l为测量通道的偏差系数。

综上所述，新建立的磁梯度数据补偿综合模型不仅可以有效的滤除飞行测量平台在实际航磁测量时可能存在的多种类型的干扰噪声，而且基于信号线性空间的思想，该模型也有利于进一步的扩展和修正，只要能分析得到与干扰信号频带相同的特征信号即可用于扩展磁测干扰场空间域，通过求解相关干扰场的滤波补偿系数，消除干扰噪声而对测量数据进行有效补偿。

步骤104：对所述补偿后的磁梯度数据进行低通滤波处理，得到处理后的磁梯度数据；

通过上述一系列的数据补偿处理后，补偿后的数据还存在一部分在仪器输出条件下无法滤除的高频噪声，因此需要根据实际应用问题的信号频带特征，以及已知的传感器芯片的本底噪声频带特征，对补偿后的数据进一步进行低通滤波处理。在本发明的实施例中，已知航空测量的地磁场的有效信号频率应小于10Hz，因此最后通过巴特沃斯低通滤波方法对补偿数据进行滤波处理，获得最终的磁梯度数据处理结果。

步骤105：根据所述处理后的磁梯度数据对补偿效果进行评估。

完成上述磁补偿及数据滤波处理后，还需要对补偿质量进行评价。实施例中对补偿后数据使用如下标准化参数对补偿结果进行评价：均方根RMS、标准差σ_c、改善比IR。参数值的具体定义如下所示：

其中，d^c为补偿及滤波处理后的数据，d⁰为原始的采集数据，N为同一类型的数据个数。

图6为本发明实施例超导航磁梯度张量数据抑噪处理系统结构图。如图6所示，一种超导航磁梯度张量数据抑噪处理系统，包括：

获取模块201，用于获取在飞行测量平台坐标系下的大地背景参考磁场值、大地背景参考磁梯度场值和飞行测量线上的背景参考磁场值的时间变化率；

磁场校正补偿处理数据确定模块202，用于根据所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对三分量磁强计数据进行补偿校正处理，得到磁场校正补偿处理数据；

磁梯度数据补偿模型确定模块203，用于根据所述磁场校正补偿处理数据、所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型；

补偿磁梯度数据确定模块204，用于根据所述磁梯度数据补偿模型，得到补偿后的磁梯度数据；

磁梯度数据处理模块205，用于对所述补偿后的磁梯度数据进行低通滤波处理，得到处理后的磁梯度数据；

补偿效果评估模块206，用于根据所述处理后的磁梯度数据对补偿效果进行评估。

具体实施例1：

下面通过实测数据对本发明所提出的基于SQUID磁梯度计单通道测量数据补偿方法的效果进行验证和说明。本发明中涉及的SQUID磁梯度测量仪共有9个磁传感器输出通道，ch0-ch2为三分量磁强计输出通道，ch3-ch8为磁梯度计输出通道，本实施例使用ch0-ch3通道测量的数据对磁梯度数据(ch3)补偿效果进行说明。

图2为在飞行测量坐标系下对三分量磁强计测量值的校正补偿结果示意图，其中，(a)表示B_x分量校正前、校正后与相应参考场数据的对比结果；(b)表示B_y分量校正前、校正后与相应参考场数据的对比结果；(c)表示B_z分量校正前、校正后与相应参考场数据的对比结果；(d)表示磁场模量|B|校正前、校正后与相应参考场数据的对比结果；图3为磁场三分量数据补偿校正前后，各磁场分量Fourier频率幅值谱分析的对比图，其中，(a)表示B_x分量校正前、校正后与相应参考场数据的对比结果；(b)表示B_y分量校正前、校正后与相应参考场数据的对比结果；(c)表示B_z分量校正前、校正后与相应参考场数据的对比结果；图4为在飞行测量坐标系下对G1通道的磁梯度计测量值的补偿及滤波处理结果示意图，其中，(a)表示补偿前与补偿后数据的对比结果；(b)表示进行补偿处理后，滤波前与滤波后数据的对比结果；图5为对G1通道的磁梯度计测量值进行补偿及滤波处理后，测量数据、补偿后数据、滤波后数据的Fourier频率幅值谱分析对比图。

相应的，补偿和滤波处理结果对应的质量评价指标由表1给出，结合图4、图5和表1可知，经过补偿和滤波处理后的磁梯度数据均方根小于20pT,改善比IR可达2.36e3,因此基于本发明提出的方法能够获得良好的补偿结果，可达到仪器分辨率的要求。

表1基于综合磁梯度补偿方法的处理数据评价参数统计

至此，已经结合附图对本发明的基本方法和实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明的基于SQUID航磁梯度仪单一采集通道信号的磁梯度数据补偿方法有了清楚的认识。需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (8)

1.一种超导航磁梯度张量数据抑噪处理方法，其特征在于，包括：

获取在飞行测量平台坐标系下的大地背景参考磁场值、大地背景参考磁梯度场值和飞行测量线上的背景参考磁场值的时间变化率；

根据所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对三分量磁强计数据进行补偿校正处理，得到磁场校正补偿处理数据；

根据所述磁场校正补偿处理数据、所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型；

根据所述磁梯度数据补偿模型，得到补偿后的磁梯度数据；

对所述补偿后的磁梯度数据进行低通滤波处理，得到处理后的磁梯度数据；

根据所述处理后的磁梯度数据对补偿效果进行评估。

2.根据权利要求1所述的超导航磁梯度张量数据抑噪处理方法，其特征在于，所述获取在飞行测量平台坐标系下的大地背景参考磁场值、大地背景参考磁梯度场值和飞行测量线上的背景参考磁场值的时间变化率，具体包括：

获取在大地坐标系下的大地背景参考磁场值和大地背景参考磁梯度场值；

将所述大地坐标系下的大地背景参考磁场值和大地背景参考磁梯度场值转换到飞行测量平台坐标系中，得到在飞行测量平台坐标系下的大地背景参考磁场值和大地背景参考磁梯度场值；

根据所述飞行测量平台坐标系下的大地背景参考磁场值，通过数值差分或傅里叶变换的方式，得到飞行测量线上的背景参考磁场值的时间变化率。

3.根据权利要求1所述的超导航磁梯度张量数据抑噪处理方法，其特征在于，所述根据所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对三分量磁强计数据进行补偿校正处理，得到磁场校正补偿处理数据，具体包括：

通过正则化最小二乘拟合或Huber范数拟合的方式求解方程得到第一校正补偿系数c₁；其中，A为直接获取的观测数据；

对测量数据B_m进行第一阶段校正补偿计算，得到初始处理数据

通过正则化最小二乘拟合或Huber范数拟合的方式求解方程得到第二校正补偿系数c₂；

对数据B1_m进行第二阶段补偿计算，通过公式得到磁场校正补偿处理数据。

4.根据权利要求1所述的超导航磁梯度张量数据抑噪处理方法，其特征在于，所述根据所述磁场校正补偿处理数据、所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型，具体包括：

根据所述磁场校正补偿处理数据、所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型

其中，l为通道编号；为通道测量的磁梯度值；为通道l的补偿磁梯度值；H_e＝(H_ex,H_ey,H_ez)为参考背景磁场，dH_e＝(dH_ex,dH_ey,dH_ez)为参考背景磁场值H_e的时间变化率，B_r＝(B_rx,B_ry,B_rz)为测量磁场校正补偿处理后数据，t为惯导数据采样时间，fs为对应采样频率，p,q为干扰源信号滤波基长度，相应的，A为由各类观测数据及相关计算量构成的N×M阶矩阵，x＝(s1,o1,a_1,x,a_2,x,K,a_2p+1,z,b_1,x,K,b_2p+1,z,c_1,x,K,c_2q+1,z)^T为实际计算中使用的补偿参数，其中s1＝1/s^l，o1＝ol/s^l，s^l为测量通道的尺度偏差；o^l为测量通道的偏差系数。

5.根据权利要求1所述的超导航磁梯度张量数据抑噪处理方法，其特征在于，所述根据所述磁梯度数据补偿模型，得到补偿后的磁梯度数据，具体包括：

根据所述磁梯度数据补偿模型求解补偿系数；

根据所述补偿系数，确定补偿后的磁梯度数据。

6.根据权利要求1所述的超导航磁梯度张量数据抑噪处理方法，其特征在于，所述对所述补偿后的磁梯度数据进行低通滤波处理，得到处理后的磁梯度数据，具体包括：

对所述补偿后的磁梯度数据采用巴特沃斯低通滤波方法进行低通滤波处理，得到处理后的磁梯度数据。

7.根据权利要求1所述的超导航磁梯度张量数据抑噪处理方法，其特征在于，所述根据所述处理后的磁梯度数据对补偿效果进行评估，具体包括：

根据所述处理后的磁梯度数据，采用标准化参数对补偿效果进行评估，所述标准化参数包括均方根、标准差和改善比。

8.一种超导航磁梯度张量数据抑噪处理系统，其特征在于，包括：

获取模块，用于获取在飞行测量平台坐标系下的大地背景参考磁场值、大地背景参考磁梯度场值和飞行测量线上的背景参考磁场值的时间变化率；

磁场校正补偿处理数据确定模块，用于根据所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对三分量磁强计数据进行补偿校正处理，得到磁场校正补偿处理数据；

磁梯度数据补偿模型确定模块，用于根据所述磁场校正补偿处理数据、所述大地背景参考磁场值和所述飞行测量线上的背景参考磁场值的时间变化率对磁梯度计单一通道测量数据进行补偿，得到磁梯度数据补偿模型；

补偿磁梯度数据确定模块，用于根据所述磁梯度数据补偿模型，得到补偿后的磁梯度数据；

磁梯度数据处理模块，用于对所述补偿后的磁梯度数据进行低通滤波处理，得到处理后的磁梯度数据；

补偿效果评估模块，用于根据所述处理后的磁梯度数据对补偿效果进行评估。