一种全自动的航空电磁数据调平方法及装置

1.本发明属于航空电磁数据计算领域，特别涉及一种全自动的航空电磁数据调 平方法及装置。


背景技术：

2.航空电磁探测以飞机为载体，基于法拉第电磁感应定律的地球物理探测方 法，具有工作效率高、探测面积广以及探测深度深等优点，适合于地形地质复杂 地区的矿产探测、地质普查等。航空电磁探测技术为我国资源开发提供技术支撑， 已经被列为未来矿产勘查的重要手段之一，对国民经济的可持续发展具有战略意 义。航空电磁探测通常采用“s”型飞行方式，探测过程中测区各测线飞行状态 的不一致，将会导致测区存在沿测线方向的条带状调平误差，严重影响探测结果。 为去除测区条带状误差，提高探测数据质量，需要对航空电磁测量数据进行调平 处理。
3.中国专利cn106970426a公开了“一种基于测线差分与主成分分析的航空电 磁数据调平方法”。该方法首先选择基准测线，沿测线方向对测区数据做测线差 分；基于测区差分数据，沿切割线方向进行主成分分析，获取调平误差差分形式； 最后通过迭代求取各测线调平误差，得到调平结果。该方法通过选择基准测线作 为调平标准以完成调平，并未涉及航空电磁数据的全自动调平。
4.林振民(1980)、毛玉仙(1995)提出航磁数据的拟切割线调平，该方法需 要额外飞行切割线。李文杰(2007)提出基于一维、二维滤波处理的二维自动调 平算法，结合移动平均滤波，有效去除了频率域航空电磁数据调平误差。孙东明 等(2010)和陈雄等(2011)对拟切割线数据求导，判断其突变点以确定调平区 域，结合滤波器实现数据调平。郭涛等(2013)结合总区域背景场与局部区域背 景场，利用中值滤波实现区域控制调平。该类算法需要结合滤波器，设置滤波器 参数。至今仍未见航空电磁数据的全自动调平方法。


技术实现要素：

5.本发明的目的就在于针对上述现有技术的不足，提供一种航空电磁数据的全 自动调平方法及装置，解决了传统调平方法需要切割线数据辅助、需要设置滤波 器参数或者需要选择调平基准的技术问题。
6.为实现上述目的，一种全自动的航空电磁数据调平方法包括以下步骤：
7.a、录入航空电磁测区数据，为保持数据本身的特征，无需对测区数据进行 网格化；
8.b、对各测线数据进行镜像周期延拓，去除测线数据连接处的调平误差效应；
9.c、利用格拉布斯准则对异常区域数据进行剔除处理，避免异常区域数据干 扰调平过程；
10.d、对预处理后数据进行变分模态分解，根据调平误差在测区与测线的时频 特性，引入自适应信号处理的变分模态分解，通过迭代计算搜寻本征模态模型的 最优解，确定每
个本征模态分量的中心频率及带宽，将信号的频率实现自适应剖 分和各分量的有效分离，从而精确提取低频变化的调平误差分量；
11.e、取第一个变分模态分量作为调平误差分量，把去均值后的第一个变分模 态分量作为调平误差；
12.f、去除镜像周期延拓部分的调平误差，并从原始数据中去除对应的调平误 差作为调平结果；
13.g、输出数据调平结果并成图。
14.步骤a所述为录入航空电磁测区数据，为保持数据本身的特征，无需对测区 数据进行网格化。对于待调平的测区航空电磁数据，设测区有l条测线，各测线 包含测点为n＝[n1,n2,
…
,nl]，测区数据表示为无需对测区数据进行网格化。
[0015]
步骤b所述为对各测线数据进行镜像周期延拓，去除测线数据连接处的调平 误差效应。由于测区存在调平误差，各测线数据水平不一致，若把测区各测线数 据首尾相连，将会在两条测线连接处出现跳变，为避免数据分解时的端点效应， 在测线两端分别进行镜像周期延拓。设定各测线两端的延拓测点数为此时各测线包含测点为：
[0016][0017]
各测线数据经镜像周期延拓后，测区数据表示为：
[0018][0019]
步骤c所述为利用格拉布斯准则对异常区域数据进行剔除处理，避免异常区 域数据干扰调平过程。航空电磁的异常区域数据会影响调平误差的提取，为实现 全自动的调平过程，需要对测区数据进行预处理，剔除测区的异常区域数据。以 第一条测线为例，延拓后的第一条测线数据为当测点数据与该测线数据均值 的残差满足格拉布斯准则时，认为该测点数据会影响调平误差的提取，该测点数 据由该测线的中值数据代替：
[0020][0021]
其中，为该测线数据的均值，g(num,a)为临界值系数，取决于测量次数 num和显著性水平a，在本发明中取g(num,a)为1.148，为测线数据的标准 差。测区第一条测线剔除异常区域数据后，可表示为，
[0022][0023]
按照步骤c，依次对测区各测线数据进行异常区域数据的剔除处理，并将测 区各测线数据进行首尾相连，此时，预处理后的测区数据可记为测区数据可记为
[0024]
步骤d所述为对预处理后数据进行变分模态分解。变分模态分解通过迭代计 算搜寻变分模态模型的最优解，确定每个本征模态分量的中心频率及带宽，从而 将信号的频率实现自适应剖分和各分量的有效分离。
[0025]
为了提取测区沿测线低频变化的调平误差，本发明中设定测区数据d
ext_gru
经变分模态分解后，将分解得到2个模态分量u1(t)和u2(t)。为保证变分模态分 解的良好收敛过程，引入二次惩罚项和拉格朗日乘子，具体的步骤如下：
[0026]
(1)变分模态分解参数设置：令延拓后的测区数据d
ext_gru
为待分解数据，二 次惩罚因子α取为待分解数据长度的10倍，收敛标准τ取为10-6
。
[0027]
初始化变分模态分量u(t)、变分模态分量的中心频率ω和拉格朗日乘子λ，即 令：
[0028][0029]
(2)根据设置的变分模态分解参数，循环迭代更新变分模态分量u(t)、变分 模态分量的中心频率ω和拉格朗日乘子λ，在第m次迭代时，变分模态分量计算公 式如下：
[0030][0031]
其中，δ(t)为狄拉克函数，*表示卷积。根据帕塞瓦尔定理，函数平方的和 等于其傅里叶转换式平方之和，将式(4)变换到谱域，得到下式：
[0032][0033]
其中，sgn为符号函数。利用实际信号的埃尔米特矩阵对称性，可将式(5) 中的两项写为非负频率上的半空间积分，同时令ω＝ω ω1：
[0034][0035]
根据式(6)可得到第m次迭代的变分模态分量：
[0036][0037]
变分模态分量的中心频率计算公式如下：
[0038][0039]
同样，将式(8)变换到谱域，从而得到第m次迭代的变分模态分量的中心 频率：
[0040][0041]
拉格朗日乘子更新迭代公式如下：
[0042][0043]
(3)设置变分模态分解的迭代终止条件，当达到式(11)中的一项时，迭代 终止。
[0044][0045]
(4)对变分模态分量和进行傅里叶反变换，得到时域的分解结 果u1(t)和u2(t)，至此完成测区数据d
ext_gru
的变分模态分解过程。
[0046]
步骤e所述为取去均值的第一个变分模态分量作为调平误差。在该方法中， 第一个变分模态分量取得低频的调平误差分量，对该分量进行去均值处理，得到 预处理后的测区调平误差：
[0047][0048]
步骤f所述为去除延拓部分的调平误差，并从原始数据中去除调平误差作为 调平结果。由于步骤e求得的调平误差d
ext_error
为原始数据经过镜像周期延拓后 数据d
ext_gru
所对应的调平误差，为求取原始测区数据d所对应的调平误差，从 d
ext_error
中去除延拓部分：
[0049][0050]
本发明的调平结果为从测区原始数据中去除该部分对应的调平误差(式 (13))：d
result
＝d-d
error
。
[0051]
步骤g为输出数据调平结果并成图。输出测区数据调平结果d
result
，并根据 各测点
所在坐标进行成图。
[0052]
与现有技术相比，本发明的有益效果为：本发明公开的一种全自动的航空电 磁数据调平方法可以快速、有效的去除测区数据的条带状调平误差，如图5所示； 本发明对于航空磁数据和电磁数据均适用；本发明不需要测区切割线的辅助校 正，节省了大笔费用；本发明中所使用的参数均为固定的，不需要再进行调整， 是一种完全的自动调平方法。另外，从图2和图4中可以看出，由于本发明结合 了调平误差在时域和频域的特征，调平结果中保留了数据本身的特征，调平误差 去除较为完整。
附图说明
[0053]
图1为全自动的航空电磁数据调平方法流程图；
[0054]
图2为航空磁数据测区平面图；
[0055]
图3为航空磁数据测区数据曲线图；
[0056]
图4为航空磁数据测区预处理数据曲线图；
[0057]
图5为航空磁数据测区调平结果平面图；
[0058]
图6为航空磁数据测区调平误差平面图。
具体实施方式
[0059]
为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例 中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述 的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的 实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实 施例，都应当属于本发明保护的范围。
[0060]
需要说明的是，本发明的说明书和权利要求书及上述附图中的术语“第一终 端”、“第二终端”等是用于区别类似的对象，第一、第二为本发明所详细公开的 技术方案并非顺次执行。应该理解这样使用的数据在适当情况下可以互换，以便 这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实 施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他 的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必 限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过 程、方法、产品或设备固有的其它步骤或单元。
[0061]
下面结合附图和实施例作进一步的详细说明：
[0062]
本发明一种全自动的航空电磁数据调平方法是以变分模态分解为基础的调 平方法。以一组geotech公司探测的一组航空磁数据为例，测区包含117条测线 (记为l10160-l11320)，测线间距为200m，测区测线最小测点数为569，最大 测点数为733。
[0063]
a、录入航空电磁测区数据。对于待调平的测区航空电磁数据，测区有117 条测线，各测线包含测点为n＝[687,733,
…
,614,606]，测区数据表示为 如图2所示，无需对测区数据进行网格化。
[0064]
b、对各测线数据进行镜像周期延拓。由于测区存在调平误差，各测线数据 水平不一致，若把测区各测线数据首尾相连，将会在两条测线连接处出现跳变， 为避免数据分解时的端点效应，在测线两端分别进行镜像周期延拓。设定各测线 两端的延拓测点数为n
add
＝[343,366,
…
,307,303]，此时各测线包含测点为 n
ext
＝[1373,1465,
…
,1228,1212]。各测线数据经镜像周期延拓后，测区数据表 示为：
[0065][0066]
c、利用格拉布斯准则进行异常区域数据的剔除处理。航空电磁的异常区域 数据会影响调平误差的提取，为实现全自动的调平过程，需要对测区数据进行预 处理，剔除测区的异常区域数据。以第一条测线为例，延拓后的第一条测线数据 为当测点数据与该测线数据均值的 残差满足格拉布斯准则时，认为该测点数据会影响调平误差的提取，该测点数据 由该测线的中值数据代替：
[0067][0068]
其中，为该测线数据的均值，为测线数据的标准差。测区第一条 测线剔除异常区域数据后，可表示为，
[0069][0070]
按照步骤c，依次对测区各测线数据进行异常区域数据的剔除处理，并将测 区各测线数据进行首尾相连，此时，预处理后的测区数据可记为测区数据可记为如图4所示。
[0071]
d、对预处理后数据进行变分模态分解。变分模态分解通过迭代计算搜寻变 分模态模型的最优解，确定每个本征模态分量的中心频率及带宽，从而将信号的 频率实现自适应剖分和各分量的有效分离。
[0072]
为了提取测区沿测线低频变化的调平误差，本发明中设定测区数据d
ext_gru
经变分模态分解后，将分解得到2个模态分量u1(t)和u2(t)。为保证变分模态分 解的良好收敛过程，引入二次惩罚项和拉格朗日乘子，具体的步骤如下：
[0073]
(1)变分模态分解参数设置：令延拓后的测区数据d
ext_gru
为待分解数据，二 次惩罚因子α取为待分解数据长度的10倍，收敛标准τ取为10-6
。
[0074]
初始化变分模态分量u(t)、变分模态分量的中心频率ω和拉格朗日乘子λ，即 令：
[0075][0076]
(2)根据设置的变分模态分解参数，循环迭代更新变分模态分量u(t)、变分 模态分量的中心频率ω和拉格朗日乘子λ，在第m次迭代时，变分模态分量计算公 式如下：
(t)和u2(t)，至此完成测区数据d
ext_gru
的变分模态分解过程。
[0093]
e、取去均值的第一个变分模态分量作为调平误差。在该方法中，第一个变 分模态分量取得低频的调平误差分量，对该分量进行去均值处理，得到预处理后 的测区调平误差：
[0094][0095]
f、去除延拓部分的调平误差，并从原始数据中去除调平误差作为调平结果。 由于步骤e求得的调平误差d
ext_error
为原始数据经过镜像周期延拓后数据 d
ext_gru
所对应的调平误差，为求取原始测区数据d所对应的调平误差，从 d
ext_error
中去除延拓部分：
[0096][0097]
本发明的调平结果为从测区原始数据中去除该部分对应的调平误差(式 (13))：d
result
＝d-d
error
。
[0098]
g、输出数据调平结果并成图。输出测区数据调平结果d
result
，并根据各测 点所在坐标进行成图，如图5所示，对应的调平误差d
error
如图6所示。
[0099]
实施例二
[0100]
提供了一种全自动的航空电磁数据调平装置
[0101]
处理器，适于实现如上所述的一种全自动的航空电磁数据调平方法；
[0102]
处理器采用，集成式pc端专用处理器或移动设备所专用的手机端或者平板 电脑端处理器。
[0103]
以及
[0104]
存储单元，适于存储实现如上所述的任一项所述的全自动的航空电磁数据调 平方法，所述实现如上所述的任一项全自动的航空电磁数据调平方法由处理器加 载并执行。
[0105]
本方法可以继承在便携式智能设备中，如笔记本电脑、手机和平板电脑等， 也可以集成在某些专用设备和载具中，比如专用测绘仪、飞机和航拍设备等。
[0106]
在本技术所提供的几个实施例中，应该理解到，所揭露的技术内容，可通过 其它的方式实现。其中，以上所描述的装置实施例仅仅是示意性的，例如所述单 元的划分，可以为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如 多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或 不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是 通过一些接口，单元或模块的间接耦合或通信连接，可以是电性或其它的形式。
[0107]
所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单 元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可 以分布到多个单元上。可以根据实际的需要选择其中的部分或者全部单元来实现 本实施例方案的目的。
[0108]
另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也 可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。 上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实 现。
[0109]
所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或 使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的 技术方案本质
上或者说对现有技术做出贡献的部分或者该技术方案的全部或部 分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中， 包括若干指令用以使得一台计算机设备(可为个人计算机、服务器或者网络设备 等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括： u盘、只读存储器(rom，read-only memory)、随机存取存储器(ram，randomaccess memory)、移动硬盘、磁碟或者光盘等各种可以存储程序代码的介质。
[0110]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技 术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些 改进和润饰也应视为本发明的保护范围。