CN114814963B - 一种基于b样条抽道的时间域电磁数据甚低频噪声压制方法 - Google Patents

Abstract

本发明为一种基于B样条抽道的时间域电磁数据甚低频噪声压制方法。读取时间域电磁系统关发射采集的z分量噪声数据，频谱分析以获取甚低频噪声的中心频率；基于B样条基函数构建时间域电磁数据衰减曲线的抽道方案，根据衰减曲线特性设定模型阶数及节点向量；计算B样条抽道方案的幅频响应，优化B样条节点向量使得幅频响应的零点位于甚低频噪声的中心频率附近，通过降低抽道方案在甚低频噪声中心频率附近频率响应的幅度，有针对性地压制甚低频噪声。针对甚低频噪声对时间域电磁信号的干扰问题，设计频率响应可调控的抽道方案，既保证了时间域电磁信号的提取，又能提高对甚低频噪声的压制，同时避免了外部滤波器引入的电磁信号失真。

Description

一种基于B样条抽道的时间域电磁数据甚低频噪声压制方法

技术领域

本发明属于时间域电磁数据预处理领域，具体而言为一种基于B样条抽道的时间域电磁数据甚低频噪声压制方法。

背景技术

提升仪器采集效率、降低成本、提高分辨率是地球物理方法的重要发展趋势。特别是对于能够分析深部地下信息的时间域电磁系统而言，在大范围、高精度的探测要求下，探测结果往往会受到分辨率和噪声污染的制约。因此高效、可调控的噪声抑制方法对于高分辨率时间域电磁系统的发展至关重要。甚低频(Very low frequency,VLF)通信信号在时间域电磁系统中被视为噪声。对于时间域电磁测量中随时间和位置变化的噪声，甚低频噪声是15-25kHz频带内的主要噪声源。具体而言，甚低频噪声干扰甚低频波段时间域电磁信号的提取，相当于相干噪声，无法通过叠加消除。

通常，二阶巴特沃斯滤波器用于切断高频下的甚低频噪声。同时Allard(2007年)也介绍了陷波滤波器或平均处理有助于降低甚低频噪声。然而，陷波滤波或带宽限制等消除甚低频噪声的频域方法是不够的，会在早期衰减(例如在微秒或几十微秒范围内)数据中造成失真。这很大程度限制了额外的频域滤波器的应用。

Macnae(2015年)提出，通过对传输比特的数据分析和预测，可以从数据中减去甚低频通信信号。然而，通信发射机的比特率预测和甚低频通信信号的剪除仍将受到不同地区编码和调制方法的限制。

一般假设甚低频噪声以恒定频率传输，该频率的峰值提供了在数据预处理的频率响应中压制甚低频噪声的可能性。时间域电磁信号在早期以高振幅开始，迅速衰减至平坦，并覆盖较大的动态范围。通过对时间道中的信号进行积分，并对积分时间门的多次重复进行平均，同步检测可以显著减少时间域电磁信号在所需动态范围内的存储，增强信号并抑制非相干噪声。作为同步检测的一个步骤，抽道通常采用指数增加的窗宽来对应时间域电磁信号的衰减趋势，并且通常在这些栅极宽度的时间域电磁信号的平均中执行矩形。与矩形相比，平滑的抽道形状在抑制高频噪声方面具有优势。

Larsen等人(2021年)探索了一种平滑抽道策略，该策略提高了配备模拟boxcar积分器的时间域电磁系统的信噪比。然而，单一固定形状的平滑抽道对于不同频率的甚低频噪声缺乏灵活性。具体而言，该处理方法需要频率响应的可调零点，以满足某些特定高振幅甚低频无线电发射机的抑制。

发明内容

本发明所要解决的技术问题在于提供一种基于B样条抽道的时间域电磁数据甚低频噪声压制方法，利用B样条(B-spline)曲线的光滑性和局部支持特性，在其幅频响应(AFR)中减小带宽，增加旁瓣衰减，从而在不影响信号提取的情况下提高对平稳随机噪声的抑制。重要的是，优化了B样条的节点向量，以确保AFR的零点可以位于甚低频噪声的恒定频率附近。此外，B样条抽道增强了在抽道间隔内积分时间域电磁信号时的局部低通滤波效果，可以避免额外滤波器所引入的信号失真，提高同步检测的相干噪声去除效率，同时针对不同地区不同频率的甚低频噪声压制更具灵活性。

本发明是这样实现的，一种基于B样条抽道的时间域电磁数据甚低频噪声压制方法，该方法包括：

a、读取时间域电磁系统测量的时间域电磁z分量数据，截取off-time段衰减曲线，根据衰减曲线长度设定抽道窗宽及道数；

b、读取时间域电磁系统同一地点关闭发射机测量的z分量噪声数据，对其进行傅里叶变换获取甚低频噪声的中心频率；

c、利用B样条基函数构建单一道抽道初始模型，根据衰减曲线特性设定模型阶数；

d、定义初始模型的节点向量中首节点区间间隔和尾节点区间间隔与该道抽道窗宽的比值为自变量，设定比值的取值范围。初始模型中的每一个比值都对应一个B样条抽道方案；

e、计算步骤d中各比值对应B样条抽道方案的幅频响应及其零点；

f、根据步骤b得到的中心频率设定甚低频噪声的频率范围，统计幅频响应的零点落入甚低频噪声频率范围的比值个数；

g、对于比值个数不为零的抽道方案，选择在所有甚低频噪声中心频率下频率响应振幅总和减少最大的比值，该比值所对应抽道方案作为该道最终B样条抽道方案；

h、循环步骤c-g直至完成所有道数的B样条抽道方案的构建，对步骤a中截取的衰减曲线进行抽道处理。

进一步地，步骤a所述的截取off-time段衰减曲线，根据衰减曲线长度设定抽道窗宽及道数需要根据时间域电磁数据的发射电流的关断零点作为off-time段起始点开始截取衰减曲线，各抽道窗宽满足近似等对数间隔，相邻道包含重叠。

进一步地，步骤b所述的对z分量噪声数据进行傅里叶变换获取甚低频噪声的中心频率需要根据世界各国军用甚低频长波电台的频率信息来确定z分量噪声频谱中包含的能量较高的甚低频噪声中心频率。

进一步地，步骤c所述的利用B样条基函数构建单一道抽道初始模型计算过程如下：

设U是k+1个非递减数的集合，u₀≤u₁≤…≤u_i≤…≤u_k。u_i称为节点，集合U称为节点向量，半开区间[u_i,u_i+1)是第i个节点区间。本发明构建节点区间内u_i仅出现一次，属于简单节点。节点可认为是分隔点，将区间[u₀,u_k)细分为节点区间。所有B样条基函数被假设定义域在[u₀,u_k)上。根据Cox-de Boor递推公式，第i个p阶B样条基函数B_i,p(t)定义如下：

由上式可知，基函数B_i,p(t)将由p+1个0阶基函数构成，即B_i,p(t)是通过B_i,0(t),B_i+1,0(t),...,B_i+p,0(t)计算得到的，节点区间分别为[u_i,u_i+1),[u_i,u_i+1),…,[u_i+p,u_i+p+1)。由此可知，基函数B_i,p(t)在区间[u_i+p,u_i+p+1)上非零。这种局部支撑性确保了在平滑的抽道形状中的加权有效性。那么当i＝0且所有节点均为简单节点时，p阶B样条基函数B_i,0(t)定义于p+1个连续无重叠的节点区间[u₀,u₁),[u₁,u₂),…,[u_p,u_p+1)，此时k＝p+1。

时间域电磁探测中，信号的振幅变化率随衰减时间的增大而减小。为了在抽道过程中尽可能多的提取信号的有效信息，单个抽道时间间隔内采集的信号可实时分解为局部二阶多项式。对于某抽道区间内的N个采样点，其多项式系数计算如下：

式中，C_p为抽道区间内的p阶多项式系数，x[n]为抽道区间内的采样点。零阶多项式即为该道间隔内采集信号的平均值。同时平滑抽道形状可以增强以甚低频无线电传输为主的高频噪声分量的压制效果。基于上述时间域电磁信号的分解特点和平滑抽道要求，并结合B样条基函数特有的复杂多项式形式，设计采用2阶B样条基函数作为抽道窗函数，即p＝2，节点向量U由4个节点组成，U＝[u₀,u₁,u₂,u₃]^T。

考虑到抽道过程对抽道窗宽内的信号进行积分所起到的低通滤波效果，抽道过程可以用离散时间加权平均滤波器的形式表示如下：

其中，滤波器的输出y[n]表示抽道结果，输入x[n]表示抽道间隔内时间域电磁数据的采样点，w_k表示用于平均的加权系数，N_w是归一化系数，定义为

因此对于B样条抽道，选用二阶B样条基函数如下式：

为满足时间域电磁探测off-time段衰减曲线的e指数衰减规律，节点向量U＝[u₀,u₁,u₂,u₃]^T中，四个节点分割形成三个节点区间，其中第一个节点区间[u₀,u₁)和最后一个节点区间[u₂,u₃)对数等间隔。节点向量或节点序列为非均匀序列。以单一道为例，对数等间隔可以表示为：

Δu＝lnu₁-lnu₀＝lnu₃-lnu₂

其中，Δu为首节点区间间隔和尾节点区间间隔。

进一步地，步骤d所述的首节点区间间隔和尾节点区间间隔与该道抽道窗宽的比值计算如下：

比值r的选取范围为[1％,50％)。

进一步地，步骤e所述的B样条抽道方案的幅频响应表达式如下：

其中，角频率ω表示为ω＝2πf/f_s，f为频率，f为采样率。则幅频响应的零点表示为：

z＝findpeaks(-|H(e^jω)|)

其中z表示比值r对应幅频响应的零点，findpeaks表示MATLAB函数，该函数会返回一个带有输入信号局部最大值(峰值)的向量。通过对幅频响应的幅度绝对值取负数并返回局部峰值来获得零点。

进一步地，步骤f所述的根据步骤b得到的甚低频噪声中心频率，设[f_i-500,f_i+500]为第i个甚低频噪声频率范围，当零点在该频率范围内时，将其记录为1，并计算其总和：

其中f_i表示第i个甚低频噪声的中心频率，i＝1,2,…,N_v；N_v是探测得知的中心频率的数量。

进一步地，步骤g所述的对于比值个数不为零的抽道方案，计算其在甚低频中心频率处的频率响应幅值的累加和，选取累加和绝对值最大比值，表达式如下：

当不同的比值具有相同的幅度缩减时，选择最小比值作为B样条抽道的最终参数，以保持B样条在对数尺度上的对称性。这有助于确定抽道中心时间，获得更稳定的抽道值。

进一步地，步骤h所述的循环道数由步骤a中衰减曲线长度和抽道窗宽决定，各道独立完成。

本发明与现有技术相比，有益效果在于：

本发明利用B样条基函数构建时间域电磁数据的抽道方案，通过优化B样条基函数的节点向量来调控抽道方案的幅频响应零点使其位于甚低频噪声的频率范围内，能够有针对性的降低甚低频噪声经过抽道系统后的幅值。本发明充分利用预处理过程中的抽道技术，结合时间域电磁数据的衰减规律和B样条基函数的多项式特性，基于B样条抽道方案来压制甚低频噪声，既保证了抽道处理对时间域电磁信号的提取，同时相比于传统的矩形抽道增强了甚低频噪声的压制能力。因本发明可依据当地实测的甚低频噪声的中心频率来优化B样条抽道方案，对比单一的平滑抽道方案在压制甚低频噪声方面更具有灵活性和针对性。本发明在数据预处理过程中压制噪声，避免了由于采用常用的外部频域滤波器而引入的信号失真问题。

附图说明

图1为本发明实施例提供的基于B样条抽道的时间域电磁数据甚低频噪声压制方法的流程图；

图2为本发明实施例提供的HTEM系统关闭发射机采集的z分量噪声数据的频谱；

图3为本发明实施例提供的第17道B样条抽道方案的节点区间比值；

图4(a)为第17道矩形、平滑(Larsen等，2021年)和B样条三种抽道形状；

图4(b)为第17道矩形、平滑和B样条三种抽道方案对应的幅频响应曲线。浅灰线是甚低频噪声的频谱。灰色三角形为B样条抽道位于甚低频噪声频率范围的幅频响应零点；

图5(a)为HTEM实测衰减曲线和三种抽道方案的抽道结果对比图；

图5(b)为三种抽道结果对比的晚期道局部放大图；

图6为以矩形抽道结果的标准误为归一化参数，三种抽道结果的标准误比。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明实施例提供的基于B样条抽道的时间域电磁数据甚低频噪声压制方法，首先读取HTEM系统的z分量数据并依据发射电流的关断点进行off-time段衰减曲线截取，根据衰减曲线长度设定抽道窗宽及道数；读取同一探测地点关发射采集的z分量噪声数据并进行频谱分析，查阅各国军用平台甚低频长波电台的频率信息以获取甚低频噪声的中心频率；以单一道为例，利用B样条基函数构建抽道初始模型，根据衰减曲线特性设定模型阶数；定义初始模型的节点向量中首节点区间间隔和尾节点区间间隔与该道抽道窗宽的比值为自变量，设定比值的取值范围；计算各比值对应B样条抽道方案的幅频响应及其零点；根据甚低频噪声的中心频率设定其频率范围，统计幅频响应零点落入甚低频噪声频率范围的比值个数；对于比值个数不为零的抽道方案，选择在所有甚低频噪声中心频率下频率响应振幅总和减少最大的比值，该比值所对应抽道方案作为该道B样条抽道方案；根据设定道数依次循环对应B样条抽道方案的构建过程，并对衰减曲线进行抽道处理。

a、读取时间域电磁系统测量的时间域电磁z分量数据，截取off-time段衰减曲线，根据衰减曲线长度设定抽道窗宽及道数；

b、读取时间域电磁系统同一地点关闭发射机测量的z分量噪声数据，对其进行傅里叶变换获取甚低频噪声的中心频率；

c、利用B样条基函数构建单一道抽道初始模型，根据衰减曲线特性设定模型阶数；

d、定义初始模型的节点向量中首节点区间间隔和尾节点区间间隔与该道抽道窗宽的比值为自变量，设定比值的取值范围。初始模型中的每一个比值都对应一个B样条抽道方案；

e、计算步骤d中各比值对应B样条抽道方案的幅频响应及其零点；

f、根据步骤b得到的中心频率设定甚低频噪声的频率范围，统计幅频响应的零点落入甚低频噪声频率范围的比值个数；

g、对于比值个数不为零的抽道方案，选择在所有甚低频噪声中心频率下频率响应振幅总和减少最大的比值，该比值所对应抽道方案作为该道最终B样条抽道方案；

h、循环步骤c-g直至完成所有道数的B样条抽道方案的构建，对步骤a中截取的衰减曲线进行抽道处理。

进一步地，步骤a所述的截取off-time段衰减曲线，根据衰减曲线长度设定抽道窗宽及道数需要根据时间域电磁数据的发射电流的关断零点作为off-time段起始点开始截取衰减曲线，各抽道窗宽满足近似等对数间隔，相邻道包含重叠。

进一步地，步骤b所述的对z分量噪声数据进行傅里叶变换获取甚低频噪声的中心频率需要根据世界各国军用甚低频长波电台的频率信息来确定z分量噪声频谱中包含的能量较高的甚低频噪声中心频率。

进一步地，步骤c所述的利用B样条基函数构建单一道抽道初始模型计算过程如下：

设U是k+1个非递减数的集合，u₀≤u₁≤…≤u_i≤…≤u_k。u_i称为节点，集合U称为节点向量，半开区间[u_i,u_i+1)是第i个节点区间。本发明构建节点区间内u_i仅出现一次，属于简单节点。节点可认为是分隔点，将区间[u₀,u_k)细分为节点区间。所有B样条基函数被假设定义域在[u₀,u_k)上。根据Cox-de Boor递推公式，第i个p阶B样条基函数B_i,p(t)定义如下：

由上式可知，基函数B_i,p(t)将由p+1个0阶基函数构成，即B_i,p(t)是通过B_i,0(t),B_i+1,0(t),…,B_i+p,0(t)计算得到的，节点区间分别为[u_i,u_i+1),[u_i,u_i+1),…,[u_i+p,u_i+p+1)。由此可知，基函数B_i,p(t)在区间[u_i+p,u_i+p+1)上非零。这种局部支撑性确保了在平滑的抽道形状中的加权有效性。那么当i＝0且所有节点均为简单节点时，p阶B样条基函数B_i,0(t)定义于p+1个连续无重叠的节点区间[u₀,u₁),[u₁,u₂),…,[u_p,u_p+1)，此时k＝p+1。

时间域电磁探测中，信号的振幅变化率随衰减时间的增大而减小。为了在抽道过程中尽可能多的提取信号的有效信息，单个抽道时间间隔内采集的信号可实时分解为局部二阶多项式。对于某抽道区间内的N个采样点，其多项式系数计算如下：

式中，C_p为抽道区间内的p阶多项式系数，x[n]为抽道区间内的采样点。零阶多项式即为该道间隔内采集信号的平均值。同时平滑抽道形状可以增强以甚低频无线电传输为主的高频噪声分量的压制效果。基于上述时间域电磁信号的分解特点和平滑抽道要求，并结合B样条基函数特有的复杂多项式形式，设计采用2阶B样条基函数作为抽道窗函数，即p＝2，节点向量U由4个节点组成，U＝[u₀,u₁,u₂,u₃]^T。

考虑到抽道过程对抽道窗宽内的信号进行积分所起到的低通滤波效果，抽道过程可以用离散时间加权平均滤波器的形式表示如下：

其中，滤波器的输出y[n]表示抽道结果，输入x[n]表示抽道间隔内时间域电磁数据的采样点，w_k表示用于平均的加权系数，N_w是归一化系数，定义为

因此对于B样条抽道，选用二阶B样条基函数如下式：

为满足时间域电磁探测off-time段衰减曲线的e指数衰减规律，节点向量U＝[u₀,u₁,u₂,u₃]^T中，四个节点分割形成三个节点区间，其中第一个节点区间[u₀,u₁)和最后一个节点区间[u₂,u₃)对数等间隔。节点向量或节点序列为非均匀序列。以单一道为例，对数等间隔可以表示为：

Δu＝lnu₁-lnu₀＝lnu₃-lnu₂

其中，Δu为首节点区间间隔和尾节点区间间隔。

进一步地，步骤d所述的首节点区间间隔和尾节点区间间隔与该道抽道窗宽的比值计算如下：

比值r的选取范围为[1％,50％)。

进一步地，步骤e所述的B样条抽道方案的幅频响应表达式如下：

其中，角频率ω表示为ω＝2πf/f_s，f为频率，f为采样率。则幅频响应的零点表示为：

z＝findpeaks(-|H(e^jω)|)

其中z表示比值r对应幅频响应的零点，findpeaks表示MATLAB函数，该函数会返回一个带有输入信号局部最大值(峰值)的向量。通过对幅频响应的幅度绝对值取负数并返回局部峰值来获得零点。

进一步地，步骤f所述的根据步骤b得到的甚低频噪声中心频率，设[f_i-500,f_i+500]为第i个甚低频噪声频率范围，当零点在该频率范围内时，将其记录为1，并计算其总和：

其中f_i表示第i个甚低频噪声的中心频率，i＝1,2,…,N_v；N_v是探测得知的中心频率的数量。

进一步地，步骤g所述的对于比值个数不为零的抽道方案，计算其在甚低频中心频率处的频率响应幅值的累加和，选取累加和绝对值最大比值，表达式如下：

当不同的比值具有相同的幅度缩减时，选择最小比值作为B样条抽道的最终参数，以保持B样条在对数尺度上的对称性。这有助于确定抽道中心时间，获得更稳定的抽道值。

进一步地，步骤h所述的循环道数由步骤a中衰减曲线长度和抽道窗宽决定，各道独立完成。

本发明以一组HTEM系统在中国内蒙古自治区阿拉善盟的时间域航空电磁数据为例。对上述提到的方法进行说明。

a、读取HTEM系统的z分量数据，系统主要参数如下：峰值发射电流190A，基频25Hz，采样率为192kHz。该采样率下的时间域电磁数据在15-25kHz频带内将受到甚低频噪声的严重干扰。根据发射电流截取衰减曲线长度为700个采样点共3650μs，设定抽道道数为18道，抽道窗宽满足近似等对数间隔如表1所示。所有衰减曲线均经过叠加和运动噪声去除处理。

表1第10-18道HTEM实测数据三种抽道形状的抽道窗宽。其中，矩形抽道窗宽采用传统的1/2重叠，平滑抽道和B样条抽道窗宽采用2/3重叠。

b、读取时间域电磁系统同一探测地点关闭发射机测量的z分量噪声数据，对其进行傅里叶变换如图2所示，通过查阅世界各国军用甚低频长波电台可得本组时间域航空电磁数据中甚低频噪声的中心频率为24.1kHz；

c、以第17道为例，利用B样条基函数构建抽道初始模型如下：

设U是k+1个非递减数的集合，u₀≤u₁≤…≤u_i≤…≤u_k。u_i称为节点，集合U称为节点向量，半开区间[u_i,u_i+1)是第i个节点区间。本发明构建节点区间内u_i仅出现一次，属于简单节点。节点可认为是分隔点，将区间[u₀,u_k)细分为节点区间。所有B样条基函数被假设定义域在[u₀,u_k)上。根据Cox-de Boor递推公式，第i个p阶B样条基函数B_i,p(t)定义如下：

由上式可知，基函数B_i,p(t)将由p+1个0阶基函数构成，即B_i,p(t)是通过B_i,0(t),B_i+1,0(t),…,B_i+p,0(t)计算得到的，节点区间分别为[u_i,u_i+1),[u_i,u_i+1),…,[u_i+p,u_i+p+1)。由此可知，基函数B_i,p(t)在区间[u_i+p,u_i+p+1)上非零。这种局部支撑性确保了在平滑的抽道形状中的加权有效性。那么当i＝0且所有节点均为简单节点时，p阶B样条基函数B_i,0(t)定义于p+1个连续无重叠的节点区间[u₀,u₁),[u₁,u₂),…,[u_p,u_p+1)，此时k＝p+1。

时间域电磁探测中，信号的振幅变化率随衰减时间的增大而减小。为了在抽道过程中尽可能多的提取信号的有效信息，单个抽道时间间隔内采集的信号可实时分解为局部二阶多项式。对于某抽道区间内的N个采样点，其多项式系数计算如下：

式中，C_p为抽道区间内的p阶多项式系数，x[n]为抽道区间内的采样点。零阶多项式即为该道间隔内采集信号的平均值。同时平滑抽道形状可以增强以甚低频无线电传输为主的高频噪声分量的压制效果。基于上述时间域电磁信号的分解特点和平滑抽道要求，并结合B样条基函数特有的复杂多项式形式，设计采用2阶B样条基函数作为抽道窗函数，即p＝2，节点向量U由4个节点组成，U＝[u₀,u₁,u₂,u₃]^T。

考虑到抽道过程对抽道窗宽内的信号进行积分所起到的低通滤波效果，抽道过程可以用离散时间加权平均滤波器的形式表示如下：

其中，滤波器的输出y[n]表示抽道结果，输入x[n]表示抽道间隔内时间域电磁数据的采样点，w_k表示用于平均的加权系数，N_w是归一化系数，定义为

因此对于B样条抽道，选用二阶B样条基函数如下式：

为满足时间域电磁探测off-time段衰减曲线的e指数衰减规律，节点向量U＝[u₀,u₁,u₂,u₃]^T中，四个节点分割形成三个节点区间，其中第一个节点区间[u₀,u₁)和最后一个节点区间[u₂,u₃)对数等间隔。节点向量或节点序列为非均匀序列。以单一道为例，对数等间隔可以表示为：

Δu＝lnu₁-lnu₀＝lnu₃-lnu₂

其中，Δu为首节点区间间隔和尾节点区间间隔。

d、首节点区间间隔和尾节点区间间隔与该道抽道窗宽的比值计算如下：

比值r的选取范围为[1％,50％)，步长为1％。

e、B样条抽道方案的幅频响应表达式如下：

其中，角频率ω表示为ω＝2πf/f_s，f为频率，f为采样率。则幅频响应的零点表示为：

z＝findpeaks(-|H(e^jω)|)

其中z表示比值r对应幅频响应的零点，findpeaks表示MATLAB函数，该函数会返回一个带有输入信号局部最大值(峰值)的向量。通过对幅频响应的幅度绝对值取负数并返回局部峰值来获得零点。

f、根据步骤b得到的甚低频噪声中心频率，设[f_i-500,f_i+500]为第i个甚低频噪声频率范围，当零点在该频率范围内时，将其记录为1，并计算其总和：

其中f_i表示第i个甚低频噪声的中心频率，i＝1,2,…,N_v；N_v是探测得知的中心频率的数量。

g、对于比值个数不为零的抽道方案，计算其在甚低频中心频率处的频率响应幅值的累加和，选取累加和绝对值最大比值，表达式如下：

如图3所示，以第17道为例，灰色点表示甚低频噪声频率范围内零点个数为零的比值及其对应频率响应的幅度，蓝色点表示甚低频噪声频率范围内零点个数不为零的比值及其对应频率响应的幅值，红色星号则在表示甚低频噪声频率范围内零点个数不为零，且频率响应幅值最低时对应的比值，即r＝r_B。可以得出，当首、尾节点间隔与第17道窗宽的比值为35％时，其对应B样条抽道方案在甚低频噪声中心频率24.1kHz处幅频响应幅值最小，为-114.3365dB。

为对比B样条与其他抽道形状的幅频特性，以第17道为例，绘制矩形、平滑形状和B样条的波形(图4(a))及其对应幅频响应曲线(图4(b))。。由图4(b)可以看出，矩形抽道对甚低频噪声几乎不存在压制效果，而B样条抽道相比于平滑抽道在甚低频噪声中心频率处降低约20dB。

h、重复步骤c-g直至完成所有18道B样条抽道方案的构建，并对时间域电磁衰减曲线进行抽道。B样条抽道与矩形抽道、平滑抽道的抽道对比结果如图5(a)所示。由于抽道的去噪效果受实测数据中较大的随机噪声和采样率的限制，因此在1-1.5ms内，三种抽道方案的errorbar之间的差异不明显；然而，如图5(b)所示的对比结果的局部放大图中可以看出，第17和第18道抽道结果仍然可以证明B样条抽道的标准误差相对较小，这表明在随机噪声干扰下，它在抑制甚低频噪声方面仍然优于矩形抽道和平滑抽道。图6中给出三种抽道结果的标准误对比结果。以矩形抽道结果的标准误为归一化参数。该结果充分说明了在第13-18道中B样条抽道的标准误降低更大，可以抑制甚低频噪声，同时确保有用信号的提取，并遵循时间域电磁信号的正常衰减规律。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于B样条抽道的时间域电磁数据甚低频噪声压制方法，其特征在于，该方法包括：

a、读取时间域电磁系统测量的时间域电磁z分量数据，截取off-time段衰减曲线，根据衰减曲线长度设定抽道窗宽及道数；

b、读取时间域电磁系统同一地点关闭发射机测量的z分量噪声数据，对z分量噪声数据进行傅里叶变换获取甚低频噪声的中心频率；

c、利用B样条基函数构建单一道抽道初始模型，根据衰减曲线特性设定模型阶数，计算过程如下：

设U是k+1个非递减数的集合，u₀≤u₁≤…≤u_i≤…≤u_k，u_i称为节点，集合U称为节点向量，半开区间[u_i，u_i+1)是第i个节点区间；所有B样条基函数被假设定义域在[u₀，u_k)上，根据Cox-de Boor递推公式，第i个p阶B样条基函数B_i，p(t)定义如下：

基函数B_i，p(t)将由p+1个0阶基函数构成，即B_i，p(t)是通过B_i，0(t)，B_i+1，0(t)，...，B_i+p，0(t)计算得到的，节点区间分别为[u_i，u_i+1)，[u_i，u_i+1)，...，[u_i+p，u_i+p+1)，基函数B_i，p(t)在区间[u_i+p，u_i+p+1)上非零，当i＝0且所有节点均为简单节点时，p阶B样条基函数B_i，0(t)定义于p+1个连续无重叠的节点区间[u₀，u₁)，[u₁，u₂)，...，[u_p，u_p+1)，此时k＝p+1；

单个抽道时间间隔内采集的信号实时分解为局部二阶多项式，对于某抽道区间内的N个采样点，其多项式系数计算如下：

式中，C_p为抽道区间内的p阶多项式系数，x[n]为抽道区间内的采样点，零阶多项式为该道间隔内采集信号的平均值，采用2阶B样条基函数作为抽道窗函数，即p＝2，节点向量U由4个节点组成，U＝[u₀，u₁，u₂，u₃]^T；

抽道过程用离散时间加权平均滤波器的形式表示如下：

其中，滤波器的输出y[n]表示抽道结果，输入x[n]表示抽道间隔内时间域电磁数据的采样点，w_k表示用于平均的加权系数，N_w是归一化系数，定义为

二阶B样条基函数如下式：

d、定义单一道抽道初始模型的节点向量中首节点区间间隔和尾节点区间间隔与该道抽道窗宽的比值为自变量，设定比值的取值范围，单一道抽道初始模型中的每一个比值都对应一个B样条抽道方案，其中抽道宽度为步骤1中的设定的抽道宽度；

e、计算步骤d中各比值对应、B样条抽道方案的幅频响应及其零点；

f、根据步骤b得到的中心频率设定甚低频噪声的频率范围，统计幅频响应的零点落入甚低频噪声频率范围的比值个数；

g、对于比值个数不为零的抽道方案，选择在所有甚低频噪声中心频率下频率响应振幅总和减少最大的比值，该比值所对应抽道方案作为该道最终B样条抽道方案；

h、循环步骤c-g直至完成所有步骤1中设定的道数的B样条抽道方案的构建，对步骤a中截取的衰减曲线进行抽道处理。

2.按照权利要求1所述的方法，其特征在于，步骤a所述的截取off-time段衰减曲线，根据衰减曲线长度设定抽道窗宽及道数需要根据时间域电磁数据的发射电流的关断零点作为off-time段起始点开始截取衰减曲线，各抽道窗宽满足近似等对数间隔，相邻道包含重叠。

3.按照权利要求1所述的方法，其特征在于，步骤b所述的对z分量噪声数据进行傅里叶变换获取甚低频噪声的中心频率具体为：根据世界各国军用甚低频长波电台的频率信息来确定z分量噪声频谱中包含的能量较高的甚低频噪声中心频率。

4.按照权利要求1所述的方法，其特征在于，步骤d所述的首节点区间间隔和尾节点区间间隔与该道抽道窗宽的比值计算如下：

比值r的选取范围为[1％，50％)。

5.按照权利要求4所述的方法，其特征在于，步骤e所述的B样条抽道方案的幅频响应表达式如下：

其中，角频率ω表示为ω＝2πf/f_s，f为频率，f_s为采样率，则幅频响应的零点表示为：

z＝findpeaks(-|H(e^jω)|)

其中z表示比值r对应幅频响应的零点，findpeaks表示MATLAB函数，该函数会返回一个带有输入信号局部最大值的向量，通过对幅频响应的幅度绝对值取负数并返回局部峰值来获得零点。

6.按照权利要求1所述的方法，其特征在于，步骤f所述的根据步骤b得到的甚低频噪声中心频率，设[f_m-500，f_m+500]为第m个甚低频噪声频率范围，当零点在该频率范围内时，将其记录为1，并计算其总和：

其中f_m表示第m个甚低频噪声的中心频率，m＝1，2，...，N_v；N_v是探测得知的中心频率的数量。

7.按照权利要求6所述的方法，其特征在于，步骤g所述的对于比值个数不为零的抽道方案，计算其在甚低频中心频率处的频率响应幅值的累加和，选取累加和绝对值最大比值，表达式如下：

当不同的比值具有相同的幅度缩减时，选择最小比值作为B样条抽道的最终参数，以保持B样条在对数尺度上的对称性。

8.按照权利要求1所述的方法，其特征在于，步骤h所述的循环道数由步骤a中衰减曲线长度和抽道窗宽决定，各道独立完成。

Images