CN109814163B - 一种基于扩展补偿模型的航磁张量数据抑噪方法及*** - Google Patents

Abstract

本发明公开一种基于扩展补偿模型的航磁张量数据抑噪方法及***。方法包括：建立磁梯度综合补偿模型；根据所述磁梯度综合补偿模型，得到扩展磁梯度综合补偿模型；对所述扩展磁梯度综合补偿模型进行求解，得到补偿参数；根据所述补偿参数，确定补偿后的磁梯度数据；根据所述补偿后的磁梯度数据对补偿效果进行评估。采用本发明能够有效提高补偿后磁梯度信号的信噪比，使得补偿算法能方便直接的适用于高灵敏度磁梯度传感器的数据补偿。

Description

一种基于扩展补偿模型的航磁张量数据抑噪方法及***

技术领域

本发明涉及地球物理航空磁法勘探领域，特别是涉及一种基于扩展补偿模型的航磁张量数据抑噪方法及***。

背景技术

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

由于磁梯度张量测量***采用SQUID传感器的高灵敏度特性，因此在低空勘察飞行中地磁场梯度对于传感器的影响不可以忽略，而利用传统的高空测试飞行获得的补偿系数则无法补偿低空飞行中磁梯度信号对于传感器的干扰影响，因此直接使用传统的补偿箱方法(即利用高空飞行数据获得补偿系数直接对低空勘察飞行数据进行补偿的方式)无法获得精度满足要求的磁梯度处理数据。

发明内容

本发明的目的是提供一种基于扩展补偿模型的航磁张量数据抑噪方法及***，可以有效提高最终补偿后磁梯度信号的信噪比，使得补偿算法能方便直接的适用于高灵敏度磁梯度传感器的数据补偿，并且补偿飞行操作可以在低空测量飞行高度直接进行。

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

一种基于扩展补偿模型的航磁张量数据抑噪方法，包括：

建立磁梯度综合补偿模型；

根据所述磁梯度综合补偿模型，得到扩展磁梯度综合补偿模型；

对所述扩展磁梯度综合补偿模型进行求解，得到补偿参数；

根据所述补偿参数，确定补偿后的磁梯度数据；

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

可选的，所述建立磁梯度综合补偿模型，具体包括：

获取相关参数数据，所述相关参数数据包括测量的磁梯度值、通道测量点处的真实磁梯度值、当地背景磁场、背景磁场时间变化率和三轴磁强计测量值；

根据所述相关参数数据建立磁梯度综合补偿模型Ax＝b；

其中，A为由相关参数数据

构成的N×M阶矩阵，

为测量的磁梯度值；

为通道测量点处的真实磁梯度值；H_e为当地背景磁场，dH_e为背景磁场时间变化率，B_r为补偿校正后的三轴磁强计测量值，N为梯度计采样点数；x＝(o^l,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，x为实际计算中使用的补偿参数，共M＝12p+6q+10个，b为测量信号。

可选的，所述根据所述磁梯度综合补偿模型，得到扩展磁梯度综合补偿模型，具体包括：

根据所述测量信号b，得到测量信号的时间变化率信号Db；

根据所述N×M阶矩阵A，得到相应的系数变化率矩阵DA；

根据所述测量信号b、所述测量信号的时间变化率信号Db、所述N×M阶矩阵A和所述相应的系数变化率矩阵DA，构建扩展磁梯度综合补偿模型Kx＝d，其中，K＝[A；DA]为扩展综合补偿模型系数矩阵,d＝[b,Db]为扩展数据矢量。

可选的，所述对所述扩展磁梯度综合补偿模型进行求解，得到补偿参数，具体包括：

根据所述扩展磁梯度综合补偿模型，在时间域或变换域上建立拟合目标函数；

根据目标函数，得到补偿参数。

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

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

一种基于扩展补偿模型的航磁张量数据抑噪***，包括：

第一模型建立模块，用于建立磁梯度综合补偿模型；

第二模型建立模块，用于根据所述磁梯度综合补偿模型，得到扩展磁梯度综合补偿模型；

求解模块，用于对所述扩展磁梯度综合补偿模型进行求解，得到补偿参数；

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

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

可选的，所述第一模型建立模块，具体包括：

获取单元，用于获取相关参数数据，所述相关参数数据包括测量的磁梯度值、通道测量点处的真实磁梯度值、当地背景磁场、背景磁场时间变化率和三轴磁强计测量值；

第一模型建立单元，用于根据所述相关参数数据建立磁梯度综合补偿模型Ax＝b；

其中，A为由相关参数数据

构成的N×M阶矩阵，

为测量的磁梯度值；

为通道测量点处的真实磁梯度值；H_e为当地背景磁场，dH_e为背景磁场时间变化率，B_r为补偿校正后的三轴磁强计测量值，N为梯度计采样点数；x＝(o^l,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，x为实际计算中使用的补偿参数，共M＝12p+6q+10个，b为测量信号。

可选的，所述第二模型建立模块，具体包括：

时间变化率信号确定单元，用于根据所述测量信号b，得到测量信号的时间变化率信号Db；

系数矩阵确定单元，用于根据所述N×M阶矩阵A，得到相应的系数变化率矩阵DA；

第二模型建立单元，用于根据所述测量信号b、所述测量信号的时间变化率信号Db、所述N×M阶矩阵A和所述相应的系数变化率矩阵DA，构建扩展磁梯度综合补偿模型Kx＝d；

其中，K＝[A；DA]为扩展综合补偿模型系数矩阵,d＝[b,Db]为扩展数据矢量。

可选的，所述求解模块，具体包括：

目标函数拟合单元，用于根据所述扩展磁梯度综合补偿模型，在时间域或变换域上拟合目标函数；

补偿参数确定单元，用于根据目标函数，得到补偿参数。

可选的，所述评估模块，具体包括：

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

根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明提供一种基于扩展补偿模型的航磁张量数据抑噪方法及***，可以有效提高最终补偿后磁梯度信号的信噪比，使得补偿算法能方便直接的适用于高灵敏度磁梯度传感器的数据补偿，并且本发明可以在低空测量飞行高度直接进行。通过新的补偿参数求解算法，可以有效保持目标信号中存在多个较为尖锐的小磁梯度异常信息。

附图说明

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

图1为本发明基于扩展补偿模型的航磁张量数据抑噪方法流程图；

图2为用于数值验证的磁源理论模型对应的磁化率分布值图；

图3为在测量坐标系下，基于真实的测量的飞行姿态和GPS信息计算的背景磁场值(a)、磁异常值Bx,By,Bz(b)、磁异常梯度值G1-G6(c,d)；

图4为通过磁源理论模型和等价源方法计算得到G1通道方向的测量磁梯度数据和真实磁梯度数据以及相应一阶时间变化率：(a)测量磁梯度数据；(b)真实磁梯度数据；(c)测量磁梯度的一阶时间变化率数据；(d)真实磁梯度数据的一阶时间变化率数据；

图5为各类基本干扰源信号(F0，Fv，Fa)、模拟G1通道方向真实磁梯度数据(Gi)及其一阶时间变化率信号的频谱特征(dF0，dFv，dFa，dGi)，其中基本干扰源信号分为x,y,z三个分量；

图6为各类基本干扰源信号(A)、模拟G1通道方向真实磁梯度数据(Gi)、通过等价源方法模拟得到的G1通道测量磁梯度数据(Gm)、对应扩展数据(K,d)及其一阶时间变化率信号(DA,dGm,dGi)的小波变换数据对比；

图7为利用数值模拟得到的磁梯度传感器G1测量信号以及三种补偿算法计算得到的补偿信号；

图8为磁梯度传感器G2～G6通道测量信号通过三种补偿算法计算得到的补偿信号值比较：(a)G2；(b)G3；(c)G4；(d)G5；(e)G6；

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

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

图11为本发明基于扩展补偿模型的航磁张量数据抑噪***结构图。

具体实施方式

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

本发明的目的是提供一种基于扩展补偿模型的航磁张量数据抑噪方法及***，可以有效提高最终补偿后磁梯度信号的信噪比，使得补偿算法能方便直接的适用于高灵敏度磁梯度传感器的数据补偿，并且补偿飞行操作可以在低空测量飞行高度直接进行。

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

图1为本发明基于扩展补偿模型的航磁张量数据抑噪方法流程图。如图1所示，一种基于扩展补偿模型的航磁张量数据抑噪方法，包括：

步骤101：建立磁梯度综合补偿模型；

根据一般性假设，梯度计测量的信号由真实信号、低频干扰噪声与高频电磁噪声组成，其中低频干扰噪声主要来源于外部磁源干扰(测量平台的剩磁、感磁、涡流磁场)，以及未校准的剩余传感器***偏差和寄生磁场等。因此，当考虑测量平台的干扰源主要来自吊舱内部，并且与磁传感器为刚性链接时，则磁传感器的测量梯度值可以由如下模型进行表示：

其中，l为测量通道编号；

为测量的磁梯度值；

为通道测量点处的真实磁梯度值；H_e为当地背景磁场，可以由全球地磁场模型IGRF-12模型或局部地磁场模型近似得到；dH_e＝(dH_x,dH_y,dH_z)＝(H_e(t+1/fs)-H_e(t))·fs为背景磁场时间变化率，其中fs为数据采样频率；B_r为补偿校正后的三轴磁强计测量值；

为未校准的剩余***偏差，

为外部干扰源的剩磁干扰，两者由于在一定时间内对应单一梯度计测量值均为恒定不变的常值，可综合为一个偏差系数o^l；

为外部干扰源感磁干扰的线性算子，由(2p+1)×3个系数组成；

为外部干扰源涡流干扰的线性算子，由(2p+1)×3个系数组成；h₀为寄生磁场干扰的线性算子，由(2q+1)×3个系数组成；对应各个干扰源的基函数为基本干扰源信号；n为高频Gauss白噪声。通过对上述采集信号模型进行整理，可以得到用于磁梯度数据补偿的计算模型，如下所示：

其中，A为由各类观测数据及相关计算量

构成的N×M阶矩阵，N为梯度计采样点数；x＝(o^l,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，x为实际计算中使用的补偿参数，共M＝12p+6q+10个。由于该模型考虑了各个传感器间的***延迟特性(惯导***、梯度计和磁强计)，因此对涉及不同传感器测量数据计算得到的外部干扰源的感磁、涡流干扰场，以及寄生磁场干扰进行线性滤波近似处理；当不考虑***延迟特性时，取p＝q＝0，则模型可以简化为10个补偿系数x＝(o^l,a_x,a_y,a_z,b_x,b_y,b_z,c_x,c_z,c_z)^T的计算模型，可以表示为：

步骤102：根据所述磁梯度综合补偿模型，得到扩展磁梯度综合补偿模型；

图2为用于数值验证的磁源理论模型对应的磁化率分布值图。图3为在测量坐标系下，基于真实的测量的飞行姿态和GPS信息计算的背景磁场值(a)、磁异常值Bx,By,Bz(b)、磁异常梯度值G1-G6(c,d)。图4为通过磁源理论模型和等价源方法计算得到G1通道方向的测量磁梯度数据和真实磁梯度数据以及相应一阶时间变化率：(a)测量磁梯度数据；(b)真实磁梯度数据；(c)测量磁梯度的一阶时间变化率数据；(d)真实磁梯度数据的一阶时间变化率数据；

根据理论模型计算可知，如果目标信号源为一些列不同尺寸的块状磁源组成(如图2所示)，并且有明显的边界，则通过一系列坐标转换后(大地坐标系->测量平台坐标系->单一梯度计传感器测量平行方向)，磁异常与磁梯度异常信号与其在大地坐标系下具有形态相似性，仅在磁异常幅值较大的中间部分有与飞行姿态变化相关的小幅值波动，如图3(b～d)所示。具体原因分析如下：假设大地背景磁场平稳变化，因此在远离磁异常区域磁梯度值近似为0或常值，而该区域由于磁梯度幅值很小，所以通过一系列坐标转换后不会产生明显的与飞行姿态变化相关的幅值波动，因此磁异常场在经过坐标变化前后具有明显的局部区域特征相似性，受到飞行姿态影响小。而相对传感器的干扰场激发源主要来源于大地背景磁场，并且干扰信号在飞行坐标系下受到飞行姿态的明显影响，具有与姿态信号相似的频带特征，如图4(a,b)所示。但由于目标信号幅值(磁梯度值)与传感器测得的磁干扰幅值差值巨大，且在磁异常区域具有局部信号结构相似性，因此直接通过测量的原始信号构建补偿参数求解方程，一般的数据拟合方案很容易造成对有效信号的过度拟合，因此需要进一步加大干扰源信号和目标源信号的差异性，并且利用这一差异性进行有效信号建模，降低过度拟合的风险，才能有效提高最终补偿结果的信噪比。而通过比较测量信号和目标信号的时间变化率信号，如图4(c,d)所示，可见前者与姿态变换频率一致的信号中的高频信息特征得到很好的凸显，而目标信号中对应的高频信息明显较少，两种信号经过处理后差异性进一步加大。

因此，综合上述分析，可以建立如下基于时间变化率信息的扩展磁补偿模型，可以从基本模型上降低数据过度拟合的风险：

其中，

为扩展综合补偿模型算子，

为扩展数据矢量；α表示数据平衡系数，用于平衡两部分拟合数据的权重，一般可以取

表示计算时间变化率，其离散形式可以表示为如下矩阵：

其中，对于系数矩阵的计算等价于矩阵相乘的形式

通过对新方程进行拟合计算得到补偿参数x可以更好的保持重构信号中的有效信号。

步骤103：对所述扩展磁梯度综合补偿模型进行求解，得到补偿参数；具体包括：

根据所述扩展磁梯度综合补偿模型，在时间域或变换域上建立拟合目标函数；

根据目标函数，得到补偿参数。

在时间域上建立基于L2范数的数据拟合目标函数：

根据真实目标信号的原始数据和一阶时间变化率信号在时间域内均具有一定的光滑性，因此选取L2范数进行数据拟合，并可以通过岭回归或基于主成分分析(PCA)的方法进行补偿参数期间，具体目标函数以及求解算法如下所示：

(1)基于L2范数的数据拟合目标函数为：

(2)在补偿参数无约束条件下，通过一般岭回归方法进行计算；在有约束条件下采用带约束的线性最小二乘数据拟合的方法求解补偿系数x，求解算法可采用信赖域反射算法(Trust region reflective)，有效集算法(Active set)或内点法(Interior point)等。而相应补偿参数x求解的基本计算方法，以及利用计算得到的补偿参数重构的补偿数据

可以通过如下公式表示：

x_β＝(K^TK+βI)^-1K^Td,

其中，β为岭回归参数，I为M阶单位矩阵。

(3)通过基于PCA的求解方法，在***矩阵K线性相关性较为严重时计算更为稳定，具体的计算步骤为：

a.将系数矩阵A逐列进行均值标准化处理为

其中，a_i为矩阵A的列向量，

为向量a_i的均值，N为采样点个数。

b.对

进行特征值分解C＝VΣV^T，得到特征向量矩阵V＝(v₁,v₂,K,v_M)以及对应的特征值λ₁＞λ₂＞K＞λ_M,(Σ＝diag(λ₁,λ₂,K,λ_M))

c.选取特征值大于λ_min的分量m(≤M)个，构成主成分变换补偿系数矩阵P＝AV_m,V_m＝(v₁,v₂,K,v_m)；其中，P为A的主成分矩阵，V_m＝(v₁,v₂,K,v_m)为由前m个特征向量分量构成的子矩阵，v_i为对应特征向量矩阵V中的第i个特征向量。

d.通过最小二乘法求解主成分变换补偿参数y及相应的补偿后数据

在变换域上建立基于Huber范数的数据拟合目标函数：

图5为各类基本干扰源信号(F0，Fv，Fa)、模拟G1通道方向真实磁梯度数据(Gi)及其一阶时间变化率信号的频谱特征(dF0，dFv，dFa，dGi)，其中基本干扰源信号分为x,y,z三个分量。图6为各类基本干扰源信号(A)、模拟G1通道方向真实磁梯度数据(Gi)、通过等价源方法模拟得到的G1通道测量磁梯度数据(Gm)、对应扩展数据(K,d)及其一阶时间变化率信号(DA,dGm,dGi)的小波变换数据对比。如图5和图6所示，通过对一系列模拟基本干扰源信号和真实目标信号进行Fourier和小波变换后，可以明显发现，两者一阶时间变化率信号差异性更加凸显：在Fourier变换域内，目标信号的一阶时间变化率信号与干扰源信号的差异性非常明显，如图5所示；在小波变换域内，基本干扰源信号和测量信号的小尺度成分明显多于真实目标信号，尤其是对扩展干扰源信号的小波变换信号，在各个小波变换尺度内均有明显幅值，而目标信号仅在近似信号和大尺度区域内有明显幅值，如图6所示。因此，通过合适的数据变化可以有效增大信号间的差异性，进一步降低过度拟合的风险，因此利用变换域中有效信号的稀疏性，建立如下基于Huber范数拟合的目标函数，具体方法如下：

(1)对扩展系数矩阵K和数据d进行Fourier或多尺度小波变换，设变换矩阵为W，则有：

(2)建立基于Huber范数的新数据拟合目标函数：

s.t.x∈Q＝{x|x_min≤x≤x_max}

其中，

为Huber函数。

(3)利用交替方向乘子方法(ADMM)进行计算，求解补偿系数x，具体计算方法如下所述：

在该步骤中，利用ADMM算法进行求解时需要首先将原优化问题转化为ADMM形式：

其中，ρ₁,ρ₂为增广Language惩罚系数；f为区域Q的指示函数。具体的ADMM迭代方法可以表示成如下形式：

其中，

为软阈值算子，∏_Q为在区域Q上的投影函数。整个ADMM算法的迭代过程在迭代次数达到上限时，或是在如下定义的主变量和对偶变量残量值小于预设误差限时终止：

其中，绝对误差ε^abs和相对误差ε^rel的选择需要依赖应用问题自定义。

(4)通过如上计算方案求得补偿系数x之后，使用如下公式计算变换域中磁梯度信号补偿值，并将其做逆数据变换到时域内，通过截取前N个数据值得到最终的补偿信号

需要注意的是，用于构造磁补偿模型方程系数矩阵A的数据量以及数据类型不同(与干扰源模型的建模方式相关)，形成的系数矩阵A可能会具有严重的病态性，为了提高求解参数的稳定性和补偿参数的复用性，可以利用PCA分析提取系数矩阵A的主成分矩阵P，用同样的方式构建目标函数并求解得到补偿参数y和补偿后数据

具体目标函数和补偿求解步为：

步骤104：根据所述补偿参数，确定补偿后的磁梯度数据；

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

一般对补偿后数据使用如下标准化参数对补偿结果进行评价：均方根RMS、标准差σ_c、改善比IR。参数值的具体定义如下所示：

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

当进行理论模型测试的情况下，由于已知准确信号值d^t，因此可以通过补偿后信号与准确信号的残量r^c＝d^t-d^c来量化算法的有效性，具体定义如下几个评价参数：残量均方根RMS、残量标准差σ_c、信噪比SNR，以及峰值信噪比PSNR，参数值的具体定义如下所示：

其中，残量均方根RMS、残量标准差σ_c用于度量参量误差，其值越小越好；信噪比SNR，以及峰值信噪比PSNR用于度量有效信号的重构精度，其值越大越好。

步骤101，具体包括：

获取相关参数数据，所述相关参数数据包括测量的磁梯度值、通道测量点处的真实磁梯度值、当地背景磁场、背景磁场时间变化率和三轴磁强计测量值；

根据所述相关参数数据建立磁梯度综合补偿模型Ax＝b；

其中，A为由相关参数数据

构成的N×M阶矩阵，

为测量的磁梯度值；

为通道测量点处的真实磁梯度值；H_e为当地背景磁场，dH_e为背景磁场时间变化率，B_r为补偿校正后的三轴磁强计测量值，N为梯度计采样点数；x＝(o^l,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，x为实际计算中使用的补偿参数，共M＝12p+6q+10个，b为测量信号。

步骤102，具体包括：

根据所述测量信号b，得到测量信号的时间变化率信号Db；

根据所述N×M阶矩阵A，得到相应的系数变化率矩阵DA；

根据所述测量信号b、所述测量信号的时间变化率信号Db、所述N×M阶矩阵A和所述相应的系数变化率矩阵DA，构建扩展磁梯度综合补偿模型Kx＝d；其中，K＝[A；DA]为扩展综合补偿模型系数矩阵,d＝[b,Db]为扩展数据矢量。

本发明可以有效提高最终补偿后磁梯度信号的信噪比，使得补偿算法能方便直接的适用于高灵敏度磁梯度传感器的数据补偿，并且本发明可以在低空测量飞行高度直接进行。通过新的补偿参数求解算法，可以有效保持目标信号中存在多个较为尖锐的小磁梯度异常信息。

下面通过理论模型数据和实测数据对本发明所磁梯度数据补偿方法的效果进行验证和说明。本发明中涉及的SQUID磁梯度测量仪共有6个磁梯度传感器输出通道：G1-G6。其中，在实测数据验证中，本实施例使用G1通道测量的数据对磁梯度数据补偿效果进行说明。

首先通过仿真数据对本发明所提出的磁梯度数据补偿方法的效果进行验证和说明。假设基岩的磁化率值非常小，为1e-8(SI)，内嵌入三个磁化率值不同、形态不同的磁源体，从左到右磁化率值分为3e-4，4e-4，1e-3(SI),如图2所示。在一定观测高度下，通过惯导***和GPS真实测量的姿态位置信息，将磁异常信号坐标转换到飞行测量坐标系下，获得近似大地背景磁场信号F₀、磁异常场三分量F_v和磁梯度分量信号Gi同样如图4所示，设定磁异常值为F_v,则相应的局部磁场数据可以表示为F_a＝F₀+F_v。假设测量平台的干扰源等价于一单位偶极子，利用假设的剩磁、感磁及涡流参数以及不平衡度系数，则通过如下的模型构建可构建包含飞机机动噪声的实际测量梯度值：

其中，T为大地坐标系到传感器测量平面的坐标转换矩阵F_a＝(F_ax,F_ay,F_az)为局部磁场三分量数据，dF_a＝(dF_ax,dF_ay,dF_az)为局部磁场三分量数据的时间变化率数据。

图7为利用数值模拟得到的磁梯度传感器G1测量信号以及三种补偿算法计算得到的补偿信号。图8为磁梯度传感器G2～G6通道测量信号通过三种补偿算法计算得到的补偿信号值比较：(a)G2；(b)G3；(c)G4；(d)G5；(e)G6。模拟的G1通道测量数据如图7(a)实线所示。其中使用的用于构建补偿模型系数矩阵的数据包括：F₀,dF₀。通过三种数据拟合目标函数求解方法：岭回归，PCA分析，本发明方法，获得补偿后的磁梯度数据如图7,8所示，其中图7为G1通道的补偿结果比较，图8为磁梯度传感器G2～G6通道测量信号通过三种补偿算法计算得到的补偿信号值比较；相应的补偿后数据的各类评价指标如表1～表3所示。通过一系列计算结果显示，新方法相较于传统方法可以有效保持有效信号的幅值和细节特征，并且评价参数各项值均要明显优于其他两种方法,除了G5通道方向的PSNR值。但从补偿数据结果图(图8(d))的视觉观察可见，新方法的补偿数据结果明显与真实值更为接近，数据无效震荡更少。

表1基于岭回归数据拟合求解方案的磁梯度补偿方法的处理数据评价参数统计

表2基于PCA分析数据拟合求解方案的磁梯度补偿方法的处理数据评价参数统计

表3基于新方法数据拟合求解方案的磁梯度补偿方法的处理数据评价参数统计

然后，通过实测数据对本发明所磁梯度数据补偿方法的效果进行验证和说明。采样数据高度为1km，因此不包含有效磁梯度信号。图9为在飞行测量坐标系下通过本发明算法对G1通道的磁梯度计测量值的补偿及滤波处理结果示意图。其中,(a)表示补偿前与补偿后数据的对比结果；(b)表示进行补偿处理后，滤波前与滤波后数据的对比结果。图10为对G1通道的磁梯度计测量值进行补偿及滤波处理后，测量数据、补偿后数据、滤波后数据的Fourier频率幅值谱分析对比图。本发明的磁补偿方法可以有效的压制测量数据中的低频干扰，而补偿数据通过滤波处理后可以进一步消除其中的高频噪声，提高补偿精度。相应的，补偿和滤波处理结果对应的质量评价指标由表4给出，结合图9、图10和表4可知，经过补偿和滤波处理后的磁梯度数据均方根小于20pT,改善比IR可达2.0e3,因此基于本发明提出的方法能够获得良好的补偿结果，可达到仪器分辨率的要求。

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

图11为本发明基于扩展补偿模型的航磁张量数据抑噪***结构图。如图11所示，一种基于扩展补偿模型的航磁张量数据抑噪***，包括：

第一模型建立模块201，用于建立磁梯度综合补偿模型；

第二模型建立模块202，用于根据所述磁梯度综合补偿模型，得到扩展磁梯度综合补偿模型；

求解模块203，用于对所述扩展磁梯度综合补偿模型进行求解，得到补偿参数；

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

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

所述第一模型建立模块201，具体包括：

获取单元，用于获取相关参数数据，所述相关参数数据包括测量的磁梯度值、通道测量点处的真实磁梯度值、当地背景磁场、背景磁场时间变化率和三轴磁强计测量值；

第一模型建立单元，用于根据所述相关参数数据建立磁梯度综合补偿模型Ax＝b；

其中，A为由相关参数数据

构成的N×M阶矩阵，

为测量的磁梯度值；

为通道测量点处的真实磁梯度值；H_e为当地背景磁场，dH_e为背景磁场时间变化率，B_r为补偿校正后的三轴磁强计测量值，N为梯度计采样点数；x＝(o^l,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+10个，b为测量信号。

所述第二模型建立模块202，具体包括：

时间变化率信号确定单元，用于根据所述测量信号b，得到测量信号的时间变化率信号Db；

系数矩阵确定单元，用于根据所述N×M阶矩阵A，得到相应的系数变化率矩阵DA；

第二模型建立单元，用于根据所述测量信号b、所述测量信号的时间变化率信号Db、所述N×M阶矩阵A和所述相应的系数变化率矩阵DA，构建扩展磁梯度综合补偿模型Kx＝d；其中，K＝[A；DA]为扩展综合补偿模型系数矩阵,d＝[b,Db]为扩展数据矢量。

所述求解模块203，具体包括：

目标函数拟合单元，用于根据所述扩展磁梯度综合补偿模型，在时间域或变换域上建立拟合目标函数；

补偿参数确定单元，用于根据目标函数，得到补偿参数。

所述评估模块205，具体包括：

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

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

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

Claims (6)

1.一种基于扩展补偿模型的航磁张量数据抑噪方法，其特征在于，包括：

建立磁梯度综合补偿模型；

根据所述磁梯度综合补偿模型，得到扩展磁梯度综合补偿模型；

对所述扩展磁梯度综合补偿模型进行求解，得到补偿参数，具体包括：

根据所述扩展磁梯度综合补偿模型，在时间域或变换域上建立拟合目标函数；

根据目标函数，得到补偿参数；

根据所述补偿参数，确定补偿后的磁梯度数据；

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

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

根据所述扩展磁梯度综合补偿模型，在变换域上建立拟合目标函数，具体包括：

(1)对扩展系数矩阵K和数据d进行Fourier或多尺度小波变换，设变换矩阵为W，则有：

(2)建立基于Huber范数的新数据拟合目标函数：

s.t.x∈Q＝{x|x_min≤x≤x_max}

其中，

为Huber函数；

(3)利用交替方向乘子方法(ADMM)进行计算，求解补偿系数x，具体计算方法如下所述：

在该步骤中，利用ADMM算法进行求解时需要首先将原优化问题转化为ADMM形式：

其中，ρ₁,ρ₂为增广Language惩罚系数；f为区域Q的指示函数，具体的ADMM迭代方法表示成如下形式：

其中，

为软阈值算子，Π_Q为在区域Q上的投影函数，整个ADMM算法的迭代过程在迭代次数达到上限时，或是在如下定义的主变量和对偶变量残量值小于预设误差限时终止：

其中，绝对误差ε^abs和相对误差ε^rel的选择需要依赖应用问题自定义；

(4)通过如上计算方案求得补偿系数x之后，使用如下公式计算变换域中磁梯度信号补偿值，并将其做逆数据变换到时域内，通过截取前N个数据值得到最终的补偿信号

为了提高求解参数的稳定性和补偿参数的复用性，利用PCA分析提取系数矩阵A的主成分矩阵P，用同样的方式构建目标函数并求解得到补偿参数y和补偿后数据

具体目标函数和补偿求解步为：

2.根据权利要求1所述的基于扩展补偿模型的航磁张量数据抑噪方法，其特征在于，所述建立磁梯度综合补偿模型，具体包括：

获取相关参数数据，所述相关参数数据包括测量的磁梯度值、当地背景磁场、背景磁场时间变化率和三轴磁强计测量值；

根据所述相关参数数据建立磁梯度综合补偿模型Ax＝b；

其中，A为由相关参数数据

构成的N×M阶矩阵，

为测量的磁梯度值；H_e为当地背景磁场，dH_e为背景磁场时间变化率，B_r为补偿校正后的三轴磁强计测量值，N为梯度计采样点数；x为实际计算中使用的补偿参数，共M＝12p+6q+10个，b为测量信号。

3.根据权利要求2所述的基于扩展补偿模型的航磁张量数据抑噪方法，其特征在于，所述根据所述磁梯度综合补偿模型，得到扩展磁梯度综合补偿模型，具体包括：

根据所述测量信号b，得到测量信号的时间变化率信号Db；

根据所述N×M阶矩阵A，得到相应的系数变化率矩阵DA；

根据所述测量信号b、所述测量信号的时间变化率信号Db、所述N×M阶矩阵A和所述相应的系数变化率矩阵DA，构建扩展磁梯度综合补偿模型Kx＝d，其中，K＝[A,DA]为扩展综合补偿模型系数矩阵,d＝[b,Db]为扩展数据矢量。

4.一种基于扩展补偿模型的航磁张量数据抑噪***，其特征在于，包括：

第一模型建立模块，用于建立磁梯度综合补偿模型；

第二模型建立模块，用于根据所述磁梯度综合补偿模型，得到扩展磁梯度综合补偿模型；

求解模块，用于对所述扩展磁梯度综合补偿模型进行求解，得到补偿参数；

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

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

所述求解模块，具体包括：

目标函数拟合单元，用于根据所述扩展磁梯度综合补偿模型，在时间域或变换域上拟合目标函数；

补偿参数确定单元，用于根据目标函数，得到补偿参数；

根据所述扩展磁梯度综合补偿模型，在变换域上建立拟合目标函数，具体包括：

(1)对扩展系数矩阵K和数据d进行Fourier或多尺度小波变换，设变换矩阵为W，则有：

(2)建立基于Huber范数的新数据拟合目标函数：

s.t.x∈Q＝{x|x_min≤x≤x_max}

其中，

为Huber函数；

(3)利用交替方向乘子方法(ADMM)进行计算，求解补偿系数x，具体计算方法如下所述：

在该步骤中，利用ADMM算法进行求解时需要首先将原优化问题转化为ADMM形式：

其中，ρ₁,ρ₂为增广Language惩罚系数；f为区域Q的指示函数，具体的ADMM迭代方法表示成如下形式：

其中，

为软阈值算子，Π_Q为在区域Q上的投影函数，整个ADMM算法的迭代过程在迭代次数达到上限时，或是在如下定义的主变量和对偶变量残量值小于预设误差限时终止：

其中，绝对误差ε^abs和相对误差ε^rel的选择需要依赖应用问题自定义；

(4)通过如上计算方案求得补偿系数x之后，使用如下公式计算变换域中磁梯度信号补偿值，并将其做逆数据变换到时域内，通过截取前N个数据值得到最终的补偿信号

为了提高求解参数的稳定性和补偿参数的复用性，利用PCA分析提取系数矩阵A的主成分矩阵P，用同样的方式构建目标函数并求解得到补偿参数y和补偿后数据

具体目标函数和补偿求解步为：

所述评估模块，具体包括：

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

5.根据权利要求4所述的基于扩展补偿模型的航磁张量数据抑噪***，其特征在于，所述第一模型建立模块，具体包括：

获取单元，用于获取相关参数数据，所述相关参数数据包括测量的磁梯度值、当地背景磁场、背景磁场时间变化率和三轴磁强计测量值；

第一模型建立单元，用于根据所述相关参数数据建立磁梯度综合补偿模型Ax＝b；

其中，A为由相关参数数据

构成的N×M阶矩阵，

为测量的磁梯度值；H_e为当地背景磁场，dH_e为背景磁场时间变化率，B_r为补偿校正后的三轴磁强计测量值，N为梯度计采样点数；x为实际计算中使用的补偿参数，共M＝12p+6q+10个，b为测量信号。

6.根据权利要求5所述的基于扩展补偿模型的航磁张量数据抑噪***，其特征在于，所述第二模型建立模块，具体包括：

时间变化率信号确定单元，用于根据所述测量信号b，得到测量信号的时间变化率信号Db；

系数矩阵确定单元，用于根据所述N×M阶矩阵A，得到相应的系数变化率矩阵DA；

第二模型建立单元，用于根据所述测量信号b、所述测量信号的时间变化率信号Db、所述N×M阶矩阵A和所述相应的系数变化率矩阵DA，构建扩展磁梯度综合补偿模型Kx＝d；其中，K＝[A,DA]为扩展综合补偿模型系数矩阵,d＝[b,Db]为扩展数据矢量。

Images