具体实施方式

下面结合附图与实施例对本发明作进一步说明。

应该指出，以下详细说明都是例示性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

实施例一

本实施例提供了一种非全时半航空瞬变电磁数据运动噪声去除方法，如图2所示，包括以下步骤：

步骤1：确定时域二次场样本数据的采样频率、样本点个数、样本数据周期和每个周期样本点个数。

步骤2：确定所有样本点对应的真实时刻，将观测到的二次场样本数据延拓至全时长。首先将响应剖面各样本点根据其所属的半周期进行分割，假设共a个样本点，b个半周期，则每个半周期有a/b个样本点。各个样本点之间时间间隔一致，每两个相邻半周期之间存在一个半周期的时间间隔。假设一个半周期时间间隔为单位“1”，将实测第1、2、3……b个半周期样本数据分别移动到第1、3、5……2b-1半周期位置处，第2、4、6……2b位置不做任何处理。

步骤3：根据所有半周期数据叠加曲线确定晚期数据范围，通过勒让德多项式对延拓至全时长的二次场晚期数据进行初步拟合，对拟合得到的曲线进行频谱分析，将得到的频谱进行能量累加，累加范围0Hz～1kHz，以80％累加能量和为界限确定期望频率(含运动噪声的频率)最大值。

步骤4：根据上一步确定的期望频率，基于傅里叶级数构造半航空瞬变电磁晚期数据运动噪声超定线性方程组，通过最小二乘反演求解全时长的运动噪声，并将所获得运动噪声从原数据中剔除。

所述步骤1中，以广西达墨隧道半航空瞬变电磁勘探一点的观测数据为例。采样频率为30000Hz，每个半周期300个采样点，共32个半周期，9600个采样点。

所述步骤2中，如图3中(a)图所示，将二次场数据延拓至全时长，将采集到的数据分段(每半个周期)对应到采集的真实时刻，上一个半周期未点与下一个半周期起点之间不做任何插值，再进行后续去噪处理。更具体地说，将第1～300、第301～600……第9301～9600样本点分别移动到第1～300、第601～900……第18901～19200样本点位置。

如图4中(a)图所示，观测到的瞬变电磁响应信号呈指数衰减，可以根据其特性将其衰减曲线划分为三个阶段：早期、中期、晚期，各阶段信号的特点如下：早期，信号强、能量大，衰减很快，有效信号和噪声重叠；中期，有效电磁响应信号逐渐减小，能量变弱，噪声影响逐渐变大；晚期，有效电磁响应信号很小，能量微弱，几乎被噪声淹没。在实测过程中，当一次脉冲场突然关断后，二次场的产生需要一定的过程。因此，实测每个半周期数据二次场尖峰前存在能量增强段，最前端低能量平滑区为噪声部分，与衰减曲线晚期信号特点相似，参见图4中(b)图。因此，采集到的二次场数据时间上不连续，每两个相邻半周期终点和起点之间存在明显的阶跃，导致运动噪声以及有效电磁数据的连续性消失。

所述步骤3中，根据二次场信号衰减曲线分三个阶段的特性，认为衰减晚期有效信号微弱，几乎被噪声淹没。因此，拟合运动噪声的基本点仅为每个半周期二次场衰减晚期数据点，早期、中期的数据不用于拟合运动噪声，保证早、中期的衰减特性不对拟合运动噪声的发展趋势产生影响。

首先圈定晚期数据的范围：将所有半周期数据叠加，观察各样本点位置的能量大小，将低能量平滑区定为晚期数据。然后，读取每个半周期衰减晚期数据以及脉冲前端部分数据并存储在特定矩阵中，利用这些数据点作为已知噪声点拟合全时长的运动噪声，减小早、中期数据衰减趋势对拟合运动噪声产生的影响。图4中(b)图给出所有半周期叠加后曲线，认为第11～50样本点为衰减早、中期数据。将各个半周期对应位置数据剔除，得到晚期数据，参见图3中(b)图。

初步拟合运动噪声须确定勒让德多项式的最高阶数。勒让德多项式最高阶数的选择，以衰减晚期数据去除多项式拟合曲线后剩余噪声平均能量小于去噪前相对于均值的平均能量的10％为准。

初步拟合基线后，对拟合基线进行频谱分析。运动噪声频率主要集中在0～1000Hz，将此范围内的频率能量值从小到大依次进行累加，当累加值超过0～1000Hz频率能量值总和的80％时，将当前频率记为期望频率的最大值。

如图5所示，采用10阶勒让德多项式初步拟合运动噪声。对得到的基线进行频谱分析及能量累加，极限频谱分析结果和0～1000Hz能量累加结果分别如图6中(a)和(b)所示，确定期望频率的最大值为132Hz，频率间隔为二分之采样频率除以样本点总个数，即(30000/2)/9600Hz。将0-132Hz的期望频率作为后续超定线性方程组的参数。

值得注意的是，此方法同样适用于全时数据。

所述步骤4中，如图7所示，通过离散傅立叶逆变换构建针对离散信号的傅立叶级数，进而构造超定线性方程组，求解期望频率频率域系数值。然后，将求解得到的期望频率系数向量x中各个元素放到在频率域对应的位置，利用离散傅里叶逆变换反演得到全时长的时间域运动噪声，并将其从对应位置的非全时电磁数据中去掉以达到去除运动噪声的目的，去除运动噪声后的结果如图7中(b)图所示。去噪数据与原始数据叠加曲线对比结果如图8所示。

基于傅里叶级数进行基线拟合，如图7中(a)图所示，可以用其表示任意信号，对表示周期性信号或者近似周期信号具有天然的优势。

基于傅里叶级数的超定线性方程组的构建：

IDFT基本公式：

式中，f(m)是在第m位置的时间域幅值，F(k)是在第k位置的频率域系数值，N是样本点的总个数，w＝e^2πi/N，

将IDFT基本公式展开为矩阵形式：

假设等式中f(x)为运动噪声的连续采样，则F(y)为运动噪声对应的离散频率域系数，但事实上只有部分运动噪声是准确知晓的，即瞬变电磁晚期数据，而在早期，有效信号与运动噪声重叠是无法直接分离的，运动噪声一般以中低频为主，能量亦主要集中于中低频段，因此刻画运动噪声并不需要全部频率域系数，只需要部分频率域系数即可，尤其是低频段对应频率域系数，如果噪声在时间域已知点的数目大于频率域未知系数的数目，则可以通过构建超定方程组求解噪声的频率域系数。

如式(3)所示，不妨设运动噪声在频率域中有两个期望频率(实际可能有更多的频率)，即四个频率域系数，则可以通过选择时间域中任意四个值准确已知的点来构造一个适定的线性方程组求得唯一解。“期望频率”是指那些包含线圈运动噪声的频率。

式中，右边变量F(1)、F(2)、F(N-2)和F(N-1)为期望频率，即运动噪声对应频率。左边变量f(a)、f(b)、f(m)和f(n)为在时间域中所选中数据点。

由于实际数据中除运动噪声外仍存在高斯噪声等其他噪声，所以需要在时间域中选择更多的点构造一个超定线性方程组来求解运动噪声。假设选择了六个时间域数据点，有四个待求的频率系数，超定线性方程组如式(4)所示：

式中，F(k₁)是的复共轭，m₁……m₆代表样本点在时间域的不同位置，k₁、k₂、代表期望频率在频率域的不同位置。

将式(4)改写为一个矩阵向量公式，如式(5)所示：

Ax＝b (5)

其中

在式(5)两边左乘A的转置矩阵A^T，得到

A^TAx＝A^Tb (7)

在式(7)两边左乘(A^TA)^-1，利用最小二乘求解期望频率系数向量x，得到

x＝(A^TA)^-1A^Tb (8)

然后，将求解得到的期望频率系数向量x中各个元素放到在频率域对应的位置，利用离散傅里叶逆变换反演得到全时长的时间域运动噪声基线，并将其从对应位置的非全时电磁数据中去掉以达到去除运动噪声的目的，运动噪声构造方法如式(9)所示，即不同4个频率系数构造整个时间序列噪声：

更一般的，当实测瞬变电磁数据均值明显偏移二次场能量0位置，可以在上述式(6)矩阵A最右侧加入一列一阶勒让德多项式，修正曲线的偏移，以得到更准确的结果。

一阶勒让德多项式公式如下：

p₁(x)＝x (10)

式中，x∈[-1，1]，x个数等于N，各点间隔2/(N-1)。求解矩阵如下式(11)

式中p₁(m₁)……p₁(m₆)代表样本点在时间域对应位置的一阶勒让德多项式，P₁(k)代表一阶勒让德多项式系数。

实施例二

本实施例的目的是在实施例一所述方法的基础上，提供一种非全时半航空瞬变电磁数据运动噪声去除系统，所述系统包括：

采样数据获取模块，用于获取时域二次场样本数据，及其采样频率和样本数据周期；

采样数据处理模块，用于将所述二次场样本数据包含的多个半周期样本数据延拓至全时长；

期望频率获取模块，用于对全时长的运动噪声基线进行初步拟合，并通过频谱分析得到期望频率；

运动噪声拟合模块，用于根据所述期望频率，利用离散傅里叶逆变换反演得到全时长的时域运动噪声；

运动噪声剔除模块，用于从二次场样本数据中将所获得运动噪声剔除。

实施例三

本实施例的目的是提供一种电子设备。

一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现如实施例一所述的方法。

实施例四

本实施例的目的是提供一种计算机可读存储介质。

一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现如实施例一所述的方法。

以上实施例二至四中涉及的各步骤与方法实施例一相对应，具体实施方式可参见实施例一的相关说明部分。术语“计算机可读存储介质”应该理解为包括一个或多个指令集的单个介质或多个介质；还应当被理解为包括任何介质，所述任何介质能够存储、编码或承载用于由处理器执行的指令集并使处理器执行本发明中的任一方法。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。