基于混合范数拟合的航磁梯度张量数据抑噪方法及系统
阅读：1051发布：2020-06-12
专利汇可以提供基于混合范数拟合的航磁梯度张量数据抑噪方法及系统专利检索，专利查询，专利分析的服务。并且本 发明 公开一种基于混合范数拟合的航磁梯度张量数据抑噪方法及系统。方法包括：建立磁梯度数据综合补偿模型；根据所述磁梯度综合补偿模型，得到观测数据的时间变化率 信号 ；根据所述时间变化率信号，构建基于信号变化率的补偿方程；根据所述基于信号变化率的补偿方程，建立数据拟合目标函数；对所述拟合目标函数进行求解，得到补偿参数；根据所述补偿参数，确定补偿后的磁梯度数据；根据所述补偿后的磁梯度数据对补偿效果进行评估。采用本发明的方法或装置能够提高最终补偿信号的 信噪比 。，下面是基于混合范数拟合的航磁梯度张量数据抑噪方法及系统专利的具体信息内容。
1.一种基于混合范数拟合的航磁梯度张量数据抑噪方法，其特征在于，包括：
建立磁梯度数据综合补偿模型；
根据所述磁梯度综合补偿模型，得到观测数据的时间变化率信号；
根据所述时间变化率信号，构建基于信号变化率的补偿方程；
根据所述基于信号变化率的补偿方程，建立数据拟合目标函数；
对所述拟合目标函数进行求解，得到补偿参数，具体包括：
基于Lp范数数据拟合方法，获取初始的补偿参数和初始补偿数据；
计算所述初始补偿数据的一阶时间变化率；
根据所述一阶时间变化率分析有效信号的区域特征，得到分析结果；
根据所述分析结果确定拟合Lp范数阶数p和q的取值；
根据所述p和q的取值，通过重加权Gauss-Newton迭代算法求解所述数据拟合目标函数，得到补偿参数；
根据所述补偿参数，确定补偿后的磁梯度数据；
根据所述补偿后的磁梯度数据对补偿效果进行评估，具体包括：
根据所述补偿后的磁梯度数据，采用标准化参数对补偿效果进行评估，所述标准化参数包括均方根、标准差和改善比。
2.根据权利要求1所述的基于混合范数拟合的航磁梯度张量数据抑噪方法，其特征在于，所述根据所述磁梯度综合补偿模型，得到观测数据的时间变化率信号，具体包括：
根据所述磁梯度综合补偿模型，构建补偿方程Ax＝b；
根据所述补偿方程，得到观测数据的时间变化率信号DA,Db；
其中，A为直接获取的观测数据，b为用于补偿梯度值的观测数据，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
3.根据权利要求1所述的基于混合范数拟合的航磁梯度张量数据抑噪方法，其特征在于，所述根据所述时间变化率信号，构建基于信号变化率的补偿方程，具体包括：
根据所述时间变化率信号，构建基于信号变化率的补偿方程DAx＝Db；
其中，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
4.一种基于混合范数拟合的航磁梯度张量数据抑噪系统，其特征在于，包括：
补偿模型建立模块，用于建立磁梯度数据综合补偿模型；
时间变化率信号获取模块，用于根据所述磁梯度综合补偿模型，得到观测数据的时间变化率信号；
补偿方程确定模块，用于根据所述时间变化率信号，构建基于信号变化率的补偿方程；
目标函数确定模块，用于根据所述基于信号变化率的补偿方程，建立数据拟合目标函数；
求解模块，用于对所述拟合目标函数进行求解，得到补偿参数；
磁梯度数据确定模块，用于根据所述补偿参数，确定补偿后的磁梯度数据；
评估模块，用于根据所述补偿后的磁梯度数据对补偿效果进行评估；
所述求解模块，具体包括：
基于Lp范数数据拟合方法，获取初始的补偿参数和初始补偿数据；
计算所述初始补偿数据的一阶时间变化率；
根据所述一阶时间变化率分析有效信号的区域特征，得到分析结果；
根据所述分析结果确定Lp范数阶数p和q的取值；
根据所述p和q的取值，通过重加权Gauss-Newton迭代算法求解自适应Lp范数数据拟合的目标函数；
所述评估模块，具体包括：
评估单元，用于根据所述补偿后的磁梯度数据，采用标准化参数对补偿效果进行评估，所述标准化参数包括均方根、标准差和改善比。
5.根据权利要求4所述的基于混合范数拟合的航磁梯度张量数据抑噪系统，其特征在于，所述时间变化率信号获取模块，具体包括：
第一补偿方程确定单元，用于根据所述磁梯度综合补偿模型，构建补偿方程Ax＝b；
时间变化率信号获取单元，用于根据所述补偿方程，得到观测数据的时间变化率信号DA,Db；
其中，A为直接获取的观测数据，b为用于补偿梯度值的观测数据，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
6.根据权利要求4所述的基于混合范数拟合的航磁梯度张量数据抑噪系统，其特征在于，所述补偿方程确定模块，具体包括：
第二补偿方程确定单元，用于根据所述时间变化率信号，构建基于信号变化率的补偿方程DAx＝Db；
其中，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
基于混合范数拟合的航磁梯度张量数据抑噪方法及系统
技术领域
[0001] 本发明涉及地球物理航空磁法勘探技术领域，特别是涉及一种基于混合范数拟合的航磁梯度张量数据抑噪方法及系统。
背景技术
[0002] 航空磁法勘探是一种重要并广泛应用于矿产和环境勘测的地球物理勘探技术手段。近年来，基于超导量子干涉仪(SQUID)的磁梯度测量系统大力发展，使得直接测量高精度的磁梯度数据成为可能，而超导全张量磁梯度仪因其高信息量、高磁场灵敏度、体积小等诸多优点，被认为是第三代航磁探测的发展方向。目前，我国已建设有基于SQUID的全张量磁梯度测量平台原型机，不同于传统的采用机载方式进行测量的航磁梯度测量飞行平台，该原型机将磁梯度测量仪及相关设备(GPS、惯导系统(INS)、读入读出电路等)装载于一个吊舱内，并通过直升机拖曳的方式进行飞行测量，该方式将大大降低直升机对磁测量干扰的影响，甚至在拖曳缆绳足够长的条件下直升机影响可以忽略不计，但吊舱内仍然存在大量会产生磁干扰的设备，因此有效地通过磁补偿处理提高飞行平台的测量精度，以实现高分辨率磁梯度测量具有重要的意义。
[0003] 由于磁梯度张量测量系统采用SQUID传感器的高灵敏度特性，因此在低空勘察飞行中地磁场梯度对于传感器的影响不可以忽略，而利用传统的高空测试飞行获得的补偿系数则无法补偿低空飞行中磁梯度信号对于传感器的干扰影响，因此直接使用传统的补偿方法(即利用高空飞行数据获得补偿系数直接对低空勘察飞行数据进行补偿的方式)无法获得精度满足要求的磁梯度处理数据。因此针对这一问题，需要改变原有的补偿操作方式，直接对低空飞行勘察数据进行补偿。但这一过程中需要面对的重要问题是，诱导产生干扰源的磁梯度信号同时也可能是针对勘察目的的有效信号，那么如何在压制干扰的前提下保留这些有效信号，提高最终重构信号的信噪比，成为磁梯度数据处理中的一个重要的问题。
发明内容
[0004] 本发明的目的是提供一种基于混合范数拟合的航磁梯度张量数据抑噪方法及系统，能够提高最终补偿信号的信噪比。
[0005] 为实现上述目的，本发明提供了如下方案：
[0006] 一种基于混合范数拟合的航磁梯度张量数据抑噪方法，包括：
[0007] 建立磁梯度数据综合补偿模型；
[0008] 根据所述磁梯度综合补偿模型，得到观测数据的时间变化率信号；
[0009] 根据所述时间变化率信号，构建基于信号变化率的补偿方程；
[0010] 根据所述基于信号变化率的补偿方程，建立数据拟合目标函数；
[0011] 对所述拟合目标函数进行求解，得到补偿参数；
[0012] 根据所述补偿参数，确定补偿后的磁梯度数据；
[0013] 根据所述补偿后的磁梯度数据对补偿效果进行评估。
[0014] 可选的，所述根据所述磁梯度综合补偿模型，得到观测数据的时间变化率信号，具体包括：
[0015] 根据所述磁梯度综合补偿模型，构建补偿方程Ax＝b；
[0016] 根据所述补偿方程，得到观测数据的时间变化率信号DA,Db；
[0017] 其中，A为直接获取的观测数据，b为用于补偿梯度值的观测数据，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
[0018] 可选的，所述根据所述时间变化率信号，构建基于信号变化率的补偿方程，具体包括：
[0019] 根据所述时间变化率信号，构建基于信号变化率的补偿方程DAx＝Db；
[0020] 其中，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
[0021] 可选的，所述对所述拟合目标函数进行求解，得到补偿参数，具体包括：
[0022] 基于Lp范数数据拟合方法，获取初始的补偿参数和初始补偿数据；
[0023] 计算所述初始补偿数据的一阶时间变化率；
[0024] 根据所述一阶时间变化率分析有效信号的区域特征，得到分析结果；
[0025] 根据所述分析结果确定Lp范数阶数p和q的取值；
[0027] 可选的，所述根据所述补偿后的磁梯度数据对补偿效果进行评估，具体包括：
[0028] 根据所述补偿后的磁梯度数据，采用标准化参数对补偿效果进行评估，所述标准化参数包括均方根、标准差和改善比。
[0029] 一种基于混合范数拟合的航磁梯度张量数据抑噪系统，包括：
[0030] 补偿模型建立模块，用于建立磁梯度数据综合补偿模型；
[0031] 时间变化率信号获取模块，用于根据所述磁梯度综合补偿模型，得到观测数据的时间变化率信号；
[0032] 补偿方程确定模块，用于根据所述时间变化率信号，构建基于信号变化率的补偿方程；
[0033] 目标函数确定模块，用于根据所述基于信号变化率的补偿方程，建立数据拟合目标函数；
[0034] 求解模块，用于对所述拟合目标函数进行求解，得到补偿参数；
[0035] 磁梯度数据确定模块，用于根据所述补偿参数，确定补偿后的磁梯度数据；
[0036] 评估模块，用于根据所述补偿后的磁梯度数据对补偿效果进行评估。
[0037] 可选的，所述时间变化率信号获取模块，具体包括：
[0038] 第一补偿方程确定单元，用于根据所述磁梯度综合补偿模型，构建补偿方程Ax＝b；
[0039] 时间变化率信号获取单元，用于根据所述补偿方程，得到观测数据的时间变化率信号DA,Db；
[0040] 其中，A为直接获取的观测数据，b为用于补偿梯度值的观测数据，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
[0041] 可选的，所述补偿方程确定模块，具体包括：
[0042] 第二补偿方程确定单元，用于根据所述时间变化率信号，构建基于信号变化率的补偿方程DAx＝Db；
[0043] 其中，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
[0044] 可选的，所述求解模块，具体包括：
[0045] 基于Lp范数数据拟合方法，获取初始的补偿参数和初始补偿数据；
[0046] 计算所述初始补偿数据的一阶时间变化率；
[0047] 根据所述一阶时间变化率分析有效信号的区域特征，得到分析结果；
[0048] 根据所述分析结果确定Lp范数阶数p和q的取值；
[0049] 根据所述p和q的取值，通过重加权Gauss-Newton迭代算法求解自适应Lp范数数据拟合的目标函数。
[0050] 可选的，所述评估模块，具体包括：
[0051] 评估单元，用于根据所述补偿后的磁梯度数据，采用标准化参数对补偿效果进行评估，所述标准化参数包括均方根、标准差和改善比。
[0052] 根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明提供一种基于混合范数拟合的航磁梯度张量数据抑噪方法，可以有效提高最终补偿后磁梯度信号的信噪比，使得补偿算法能方便直接的适用于高灵敏度磁梯度传感器的数据补偿。附图说明
[0053] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0054] 图1为本发明基于混合范数拟合的航磁梯度张量数据抑噪方法流程图；
[0055] 图2为本发明磁化率分布值图；
[0056] 图3为本发明磁梯度传感器G1通道测量信号及不同方法补偿数据值比较；
[0057] 图4为本发明磁梯度传感器G2通道测量信号及不同方法补偿数据值比较；
[0058] 图5为本发明磁梯度传感器G3通道测量信号及不同方法补偿数据值比较；
[0059] 图6为本发明磁梯度传感器G4通道测量信号及不同方法补偿数据值比较；
[0060] 图7为本发明磁梯度传感器G5通道测量信号及不同方法补偿数据值比较；
[0061] 图8为本发明磁梯度传感器G6通道测量信号及不同方法补偿数据值比较；
[0062] 图9为本发明基于混合范数拟合的航磁梯度张量数据抑噪系统结构图；
[0063] 图10通过本发明算法对实际G1通道的磁梯度计测量值的补偿及滤波处理结果示意图，其中,(a)表示补偿前与补偿后数据的对比结果；(b)表示进行补偿处理后，滤波前与滤波后数据的对比结果；
[0064] 图11对G1通道的磁梯度计测量值进行补偿及滤波处理后，测量数据、补偿后数据、滤波后数据的Fourier频率幅度谱分析对比图。
具体实施方式
[0065] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0066] 本发明的目的是提供一种基于混合范数拟合的航磁梯度张量数据抑噪方法及系统，能够提高最终补偿信号的信噪比。
[0067] 为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0068] 图1为本发明基于混合范数拟合的航磁梯度张量数据抑噪方法流程图。如图1所示，一种基于混合范数拟合的航磁梯度张量数据抑噪方法，包括：
[0069] 步骤101：建立磁梯度数据综合补偿模型；
[0070] 根据一般性假设，梯度计测量的信号由真实信号、低频干扰噪声与高频电磁噪声组成，其中低频干扰噪声主要来源于外部磁源干扰(测量平台的剩磁、感磁、涡流磁场)，以及未校准的剩余传感器系统偏差和寄生磁场等。因此，当考虑测量平台的干扰源主要来自吊舱内部，并且与磁传感器为刚性连接时，磁传感器的测量梯度值可以由如下模型进行表示：
[0071]
[0072] 其中，l为测量通道编号； 为测量的磁梯度值； 为通道测量点处的真实磁梯度值；He为当地背景磁场，可以由全球地磁场模型IGRF-12模型或局部地磁场模型近似得到；dHe＝(dHx,dHy,dHz)＝(He(t+1/fs)-He(t))·fs为背景磁场时间变化率信号，其中fs为数据采样频率；Br为补偿校正后的三轴磁强计测量值； 为未校准的剩余系统偏差， 为外部干扰源的剩磁干扰，两者由于在一定时间内对应单一梯度计测量值均为恒定不变的常值，可综合为一个偏差系数ol； 为外部干扰源感磁干扰的线性算子，由(2p+1)×3个系数组成； 为外部干扰源涡流干扰的线性算子，由(2p+1)×3个系数组成；h0为寄生磁场干扰的线性算子，由(2q+1)×3个系数组成；n为高频Gauss白噪声。通过对上述采集信号模型进行整理，可以得到用于磁梯度数据补偿的计算模型，如下所示：
[0073]
[0074] 其中，A为由各类观测数据及相关计算量 构成的N×M阶矩阵，N为梯度计采样点数； 为实
际计算中使用的补偿参数，共M＝12p+6q+10个。由于该模型考虑了各个传感器间的系统延迟特性(惯导系统、梯度计和磁强计)，因此对涉及不同传感器测量数据计算得到的外部干扰源的感磁、涡流干扰场，以及寄生磁场干扰进行线性滤波近似处理；当不考虑系统延迟特性时，可取p＝q＝0，则模型可以简化为具有M＝10个补偿系数x＝(ol,ax,ay,az,bx,by,bz,cx,T
cz,cz)的计算模型，可以表示为：
际计算中使用的补偿参数，共M＝12p+6q+10个。由于该模型考虑了各个传感器间的系统延迟特性(惯导系统、梯度计和磁强计)，因此对涉及不同传感器测量数据计算得到的外部干扰源的感磁、涡流干扰场，以及寄生磁场干扰进行线性滤波近似处理；当不考虑系统延迟特性时，可取p＝q＝0，则模型可以简化为具有M＝10个补偿系数x＝(ol,ax,ay,az,bx,by,bz,cx,T
cz,cz)的计算模型，可以表示为：
[0075]
[0076] 步骤102：根据所述磁梯度综合补偿模型，得到观测数据的时间变化率信号；
[0077] 根据所述磁梯度综合补偿模型，构建补偿方程Ax＝b；
[0078] 根据所述补偿方程，得到观测数据的时间变化率信号DA,Db；
[0079] 其中，A为直接获取的观测数据，b为用于补偿梯度值的观测数据，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
[0080] 步骤103：根据所述时间变化率信号，构建基于信号变化率的补偿方程；具体包括：
[0081] 根据所述时间变化率信号，构建基于信号变化率的补偿方程DAx＝Db；
[0082] 其中，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
[0083] 步骤104：根据所述基于信号变化率的补偿方程，建立数据拟合目标函数；
[0084] 图2为本发明磁化率分布值图。根据理论模型计算可知，如果目标信号源为多个不同形态特征的磁源组成，则通过一系列坐标转换后(大地坐标系-测量平台坐标系-单一梯度计传感器测量平行方向)，磁异常与磁梯度异常信号与其在大地坐标系下具有形态相似性，仅在磁异常幅值较大的中间部分有与飞行姿态变化相关的小幅值波动。因此磁异常信号可以基本划分为两种区域，第一类为非磁源区域或远磁源区域的相对光滑平坦区间，以及第二类为磁源区或磁源边界区域的磁场波动区。为了保持补偿后信号的具有较高的信噪比，需要在第一类区域内尽可能降低信号的无效波动，尤其是信号中的低频波动，因为高频波动可以通过低通滤波处理而不影响最后的补偿效果，而在第二类区域内则需要尽可能保持信号的波动幅值和形态，以及区域内与飞行姿态变化相关的小信号震荡。
[0085] 一般意义上，传统的L2范数数据拟合的方式，由于假设计算值与观测量之间的误差满足Gaussian分布，因此在有效信号幅值很低时可以很好的对干扰信号进行拟合，并且使得最终的补偿信号趋向于小幅值均匀震荡或常值；但是当计算值与观测量直接的残余量具有明显幅值变化时(不满足Gaussian分布)，由于该方法会对有效信号进行无差别拟合，最终使得补偿后数据在有效信号区域附近出现明显的无效信号震荡，因而降低了补偿信号的信噪比。而Lp范数(0≤p≤1)数据拟合不仅可以有效的保持信号的阶跃性，还能降低其中重建信号边界处的无效震荡。除此之外，信号的导数信息可以进一步显示信号的边界特性，因而根据上述分析，结合Lp范数和信号的时间导数信息，构建如下的目标函数：
[0086]
[0087] 其中，p,q均为Lp范数的阶数， 或是实际计算中会根据信号的特性来自适应选取各区域使用的p,q值；α为拟合平衡系数，用于平衡两部分拟合数据的权重，一般可以取
[0088] 由于在p,q＜2情况下目标函数不可微，因此为了便于计算通常将其上述目标函数离散形式近似为一种加权L2范数的目标函数形式，即为：
[0089]
[0090] 其中，Φ(x)为目标函数，W1(x),W2(x)为范数加权矩阵，其对角元为w1ii,w2ii，具体定义方式为如下公式所示：
[0091]
[0092]
[0093] 并且有εp,εq均为一取值足够小的常量，实际计算中可以取值为1e-8.用于计算时间变化来的梯度算子矩阵表示为如下形式：
[0094]
[0095] 上述问题可以通过重加权高斯-牛顿Gauss-Newton方法进行迭代求解。
[0096] 步骤105：对所述拟合目标函数进行求解，得到补偿参数；
[0097] 步骤105具体包括：
[0098] 基于Lp范数数据拟合方法，获取初始的补偿参数和初始补偿数据；
[0099] 计算所述初始补偿数据的一阶时间变化率；
[0100] 根据所述一阶时间变化率分析有效信号的区域特征，得到分析结果；
[0101] 根据所述分析结果确定Lp范数阶数p和q的取值；
[0102] 根据所述p和q的取值，通过重加权Gauss-Newton迭代算法求解所述数据拟合目标函数，得到补偿参数。
[0103] 步骤105的具体步骤为：
[0104] (一)步骤1:
[0105] 由于所述拟合目标函数对于测量数据和测量数据变化率在不同的数据区域内采用不同的Lp拟合范数阶数，以提高最终补偿信号的信噪比，因此需要一个初始的补偿数据作为算法输入，以判断磁异常源所在的初始可能性区域。本方法中使用一般L2范数拟合的结果给出初始补偿数据，求解方法可使用岭回归算法或基于主成分分析(PCA)的求解方法，具体计算方案如下所示：
[0106] (1)基于L2范数的数据拟合目标函数为：
[0107]
[0108] (2)通过岭回归的方法可以得到相应的补偿参数xβ，并计算得到补偿数据 基本计算方法为：
[0109]
[0110]
[0111] 其中，β为岭回归参数，I为M阶单位矩阵。
[0112] (3)通过基于PCA的求解方法，在系统矩阵线性相关性较为严重时计算更为稳定，具体的计算步骤为：
[0113] (i)将系数矩阵A逐列进行均值标准化处理为
[0114]
[0115] 其中，ai为矩阵A的列向量， 为向量ai的均值，N为采样点个数。
[0116] (ii)对 进行特征值分解C＝VΣVT，得到特征向量矩阵V＝(v1,v2,K,vM)，以及对应的特征值λ1＞λ2＞K＞λM,(Σ＝diag(λ1,λ2,K,λM))
[0117] (iii)选取特征值大于λmin的分量m(≤M)个，构成主成分变换补偿系数矩阵P＝AVm,Vm＝(v1,v2,K,vm)；其中，P为A的主成分矩阵，Vm＝(v1,v2,K,vm)为由前m个特征向量分量构成的子矩阵，vi为对应特征向量矩阵V中的第i个特征向量。
[0118] (iv)通过最小二乘法求解主成分变换补偿参数y及相应的补偿后数据Grl0：
[0119]
[0120]
[0121] (二)步骤2:
[0122] 根据初始L2范数数据拟合相关计算得到的补偿后数据 计算判断标定信号中的非信号区域、光滑信号区域和非光滑信号区域，具体计算方法如下：
[0123] (1)首先对补偿后信号 及其变化率值 进行归一化处理，得到归一化数据di,bi：
[0124]
[0125] (2)设定归一化信号d,b的阈值θ1,θ2，并且根据如下规则判断并取值数据及其变化率拟合范数的阶数p,q：
[0126]
[0127] 该参数将用于下面的数据拟合计算中。
[0128] (三)步骤3:
[0129] 根据上述方法确定的拟合范数的阶数p,q，通过重加权Gauss-Newton迭代求解基于Lp混合范数数据拟合的目标函数，具体导出过程如下所示。
[0130] 为了求得目标函数的极小值，其需要满足目标函数梯度值等于0，由于原目标函数中的加权项W1(x),W2(x)与补偿参数x相关，因此可以通过一阶Taylor展开获得原问题的近似目标函数，并且进一步可以得到相应的迭代计算公式：
[0131]
[0132] 根据上述目标函数极小点存在的条件：
[0133]
[0134] 可以进一步推导得到如下的迭代计算公式：
[0135]
[0136] 其中，该线性迭代公式中的系数矩阵和右端项可以具体表示为：
[0137]
[0138]
[0139] 公式中 为范数加权矩阵，其对角元为w1ii,w2ii，具体定义方式为如下公式所示：
[0140]
[0141]
[0142] 并且有εp,εq均为一取值足够小的常量，实际计算中可以取值为1e-8。
[0143] 因此根据上述的迭代计算公式，补偿参数x及相应的补偿后数据可以通过如下算法计算得到：
[0144] (i)给出初值x0，迭代最大步数maxstep，以及停滞误差error；
[0145] (ii)doi＝1,...,maxstep；
[0146] (ii-1)更新加权系数矩阵
[0147] (ii-2)构造Gauss-Hessian矩阵Hn，以及梯度向量gn：
[0148] (ii-3)通过迭代线性方程组求解算法(如PCG/GMRES/BiCG等)或直接方法(LU分解等)求解方程组HnΔx＝-gn，得到更新变量Δx，并更新计算值xn+1/2＝xn+Δx；
[0149] (ii-4)根据约束条件，对计算更新变量进行约束区间投影xn+1＝Π|Q(xn+1/2)；
[0150] (ii-5)检验迭代是否停滞，如果满足||xn+1-xn||/||xn||＜error则迭代停滞，否则转到(ii-1)继续进行计算。
[0151] (iii)通过补偿参数x计算得到补偿后数据
[0152] 计算中需要注意的是，用于构造磁补偿模型方程系数矩阵A的数据量以及数据类型不同(与干扰源模型的建模方式相关)，形成的系数矩阵A会具有严重的病态性(即条件数很大)，为了提高求解参数的稳定性和补偿参数的复用性，可以利用上述步骤1(3)中的PCA分析方法提取系数矩阵A的主成分矩阵P，并用同样的方式构建数据拟合目标函数并求解得到补偿参数y和补偿后数据 具体目标函数和补偿数据求解步如下所示：
[0153]
[0154] 步骤106：根据所述补偿参数，确定补偿后的磁梯度数据；
[0155] 步骤107：根据所述补偿后的磁梯度数据对补偿效果进行评估。
[0156] 一般对补偿后数据使用如下标准化参数对补偿结果进行评价：均方根RMS、标准差σc、改善比IR。参数值的具体定义如下所示：
[0157]
[0158]
[0159]
[0160] 其中，dc为补偿处理后的数据，d0为原始的采集数据，N为同一类型的数据个数。
[0161] 当进行理论模型测试的情况下，由于已知准确信号值dt，因此可以通过补偿后信号与准确信号的残量rc＝dt-dc来量化算法的有效性，具体定义如下三个评价参数：残量均方根RMS、残量标准差σc、信噪比SNR，以及峰值信噪比PSNR，参数值的具体定义如下所示：
[0162]
[0163]
[0164]
[0165]
[0166] 其中，残量均方根RMS、残量标准差σc用于度量参量误差，其值越小越好；信噪比SNR，以及峰值信噪比PSNR用于度量有效信号的重构精度，其值越大越好。
[0167] 图9为本发明基于混合范数拟合的航磁梯度张量数据抑噪系统结构图。如图9所示，一种基于混合范数拟合的航磁梯度张量数据抑噪系统，包括：
[0168] 补偿模型建立模块201，用于建立磁梯度数据综合补偿模型；
[0169] 时间变化率信号获取模块202，用于根据所述磁梯度综合补偿模型，得到观测数据的时间变化率信号；
[0170] 补偿方程确定模块203，用于根据所述时间变化率信号，构建基于信号变化率的补偿方程；
[0171] 目标函数确定模块204，用于根据所述基于信号变化率的补偿方程，建立数据拟合目标函数；
[0172] 求解模块205，用于对所述拟合目标函数进行求解，得到补偿参数；
[0173] 磁梯度数据确定模块206，用于根据所述补偿参数，确定补偿后的磁梯度数据；
[0174] 评估模块207，用于根据所述补偿后的磁梯度数据对补偿效果进行评估。
[0175] 所述时间变化率信号获取模块202，具体包括：
[0176] 第一补偿方程确定单元，用于根据所述磁梯度综合补偿模型，构建补偿方程Ax＝b；
[0177] 时间变化率信号获取单元，用于根据所述补偿方程，得到观测数据的时间变化率信号DA,Db；
[0178] 其中，A为直接获取的观测数据，b为用于补偿梯度值的观测数据，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
[0179] 所述补偿方程确定模块203，具体包括：
[0180] 第二补偿方程确定单元，用于根据所述时间变化率信号，构建基于信号变化率的补偿方程DAx＝Db；
[0181] 其中，x为实际计算中使用的补偿参数，DA为直接获取的观测数据的时间变化率信号，Db为用于补偿梯度值的观测数据的时间变化率信号。
[0182] 所述求解模块205，具体包括：
[0183] 基于Lp范数数据拟合方法，获取初始的补偿参数和初始补偿数据；
[0184] 计算所述初始补偿数据的一阶时间变化率；
[0185] 根据所述一阶时间变化率分析有效信号的区域特征，得到分析结果；
[0186] 根据所述分析结果确定Lp范数阶数p和q的取值。
[0187] 根据所述p和q的取值，通过重加权Gauss-Newton迭代算法求解自适应Lp范数数据拟合的目标函数。
[0188] 所述评估模块207，具体包括：
[0189] 评估单元，用于根据所述补偿后的磁梯度数据，采用标准化参数对补偿效果进行评估，所述标准化参数包括均方根、标准差和改善比。
[0190] 具体实施例1：
[0191] 下面通过理论模型数据和实测数据对本发明所磁梯度数据补偿方法的效果进行验证和说明。本发明中涉及的SQUID磁梯度测量仪共有6个磁梯度传感器输出通道：G1-G6。其中，在实测数据验证中，本实施例使用G1通道测量的数据对磁梯度数据补偿效果进行说明。
[0192] 首先通过仿真数据对本发明所提出的磁梯度数据补偿方法的效果进行验证和说明。假设基岩的磁化率值非常小，为1e-8(SI)，内嵌入两个磁化率值不同、形态不同的磁源体，左侧倾斜板状体为非均质源，其磁化率值从上到下依次为1e-3，2e-3，3e-3(SI)，右侧垂直管状体的磁化率值为5e-3(SI)，如图3所示。在一定观测高度下，通过惯导系统和GPS真实测量的姿态位置信息，将磁异常信号坐标转换到飞行测量坐标系下，获得近似大地背景磁场信号F0、磁异常场三分量Fv和磁梯度分量信号Gi，设定磁异常值为Fv，则相应的局部磁场数据可以表示为Fa＝F0+Fv。假设测量平台的干扰源等价于一单位偶极子，利用假设的剩磁、感磁及涡流参数以及不平衡度系数，则通过如下的模型构建可构建包含飞机机动噪声的实际测量梯度值：
[0193]
[0194] 其中，T为大地坐标系到传感器测量平面的坐标转换矩阵，其中使用的用于构建补偿模型系数矩阵的数据包括：F0，dF0。通过三种数据拟合目标函数求解方法：岭回归，PCA分析，本发明方法，获得补偿后的磁梯度数据如图5，6所示，其中图3为本发明磁梯度传感器G1通道测量信号及不同方法补偿数据值比较。图3中不同曲线分别表示不同方法计算得到的补偿数据，图4为本发明磁梯度传感器G2通道测量信号及不同方法补偿数据值比较；图5为本发明磁梯度传感器G3通道测量信号及不同方法补偿数据值比较；图6为本发明磁梯度传感器G4通道测量信号及不同方法补偿数据值比较；图7为本发明磁梯度传感器G5通道测量信号及不同方法补偿数据值比较；图8为本发明磁梯度传感器G6通道测量信号及不同方法补偿数据值比较；相应的补偿后数据的各类评价指标如表1-表3所示。通过一系列计算结果显示，新方法相较于传统方法可以有效保持有效信号的幅值和细节特征，并且评价参数各项值均要明显优于其他两种方法。
[0195] 表1基于岭回归数据拟合求解方案的磁梯度补偿方法的处理数据评价参数统计[0196]
[0197] 表2基于PCA分析数据拟合求解方案的磁梯度补偿方法的处理数据评价参数统计[0198]
[0199]
[0200] 表3基于新方法数据拟合求解方案的磁梯度补偿方法的处理数据评价参数统计[0201]
[0202] 然后，通过实测数据对本发明所磁梯度数据补偿方法的效果进行验证和说明。采样数据高度为1km，因此不包含有效磁梯度信号。图10通过本发明算法对实际G1通道的磁梯度计测量值的补偿及滤波处理结果示意图，其中,(a)表示补偿前与补偿后数据的对比结果；(b)表示进行补偿处理后，滤波前与滤波后数据的对比结果。图11对G1通道的磁梯度计测量值进行补偿及滤波处理后，测量数据、补偿后数据、滤波后数据的Fourier频率幅度谱分析对比图。通过图10和图11可知，该磁补偿方法可以有效的压制测量数据中的低频干扰，而补偿数据通过滤波处理后可以进一步消除其中的高频噪声，提高补偿精度。相应的，补偿和滤波处理结果对应的质量评价指标由表4给出，结合图9～11和表4可知，经过补偿和滤波处理后的磁梯度数据均方根小于20pT，改善比IR可达2.0e3，因此基于本发明提出的方法能够获得良好的补偿结果，可达到仪器分辨率的要求。
[0203] 表4基于综合磁梯度补偿方法的处理数据评价参数统计
[0204]
[0205] 至此，已经结合附图对本发明的基本方法和实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明的基于SQUID航磁梯度仪单一采集通道信号的磁梯度数据补偿方法有了清楚的认识。需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。
[0206] 本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0207] 本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。
| 标题 | 发布/更新时间 | 阅读量 |
|---|---|---|
| 一种对偶空间反投影射频层析成像方法、定位方法及装置 | 2020-10-02 | 0 |
| 一种基于多层次方法和赋权超图的层次式FPGA布局布线方法 | 2020-10-05 | 2 |
| 一种自动识别鼻咽癌原发肿瘤的方法及装置 | 2020-10-15 | 1 |
| 一种道路瓶颈路段识别及交通拥堵传播预警系统和方法 | 2020-12-17 | 1 |
| 一种基于支持向量机的无线网络吞吐量的评估方法及装置 | 2021-02-26 | 0 |
| 一种计及支路量测数据的短期负荷协调预测方法 | 2022-08-05 | 2 |
| 一种车辆出库管理方法、装置、计算机可读存储介质及设备 | 2021-12-13 | 1 |
| 一种混合液制备的控制方法、装置以及系统 | 2023-02-26 | 0 |
| 一种基于CFD数值模拟的流体设备阻力优化算法 | 2022-12-09 | 1 |
| 公交站点位置获取方法及装置 | 2020-12-31 | 2 |
高效检索全球专利专利汇是专利免费检索,专利查询,专利分析-国家发明专利查询检索分析平台，是提供专利分析，专利查询，专利检索等数据服务功能的知识产权数据服务商。
我们的产品包含105个国家的1.26亿组数据，免费查、免费专利分析。
分析报告
专利汇分析报告产品可以对行业情报数据进行梳理分析，涉及维度包括行业专利基本状况分析、地域分析、技术分析、发明人分析、申请人分析、专利权人分析、失效分析、核心专利分析、法律分析、研发重点分析、企业专利处境分析、技术处境分析、专利寿命分析、企业定位分析、引证分析等超过60个分析角度，系统通过AI智能系统对图表进行解读，只需1分钟，一键生成行业专利分析报告。
QQ群二维码
意见反馈
意见反馈