CN105549097A - 一种瞬变电磁信号工频及其谐波干扰消除方法及装置 - Google Patents

Abstract

本发明属于电磁信号处理，具体地来讲为一种瞬变电磁信号工频及其谐波干扰消除方法及装置，该方法包括，在第一时间段内接收环境噪声信号N(t)，经放大后，经低通滤波去除高频噪声，对其信号中的工频及其谐波噪声进行建模得到工频及其谐波干扰模型，从第二时间段开始，打开发射机，接收包括瞬变电磁感应电压信号和环境噪声信号，采用所述第一时间段内的工频及其谐波干扰模型计算所述第二时间段内及后续时间段的工频及其谐波估计值，将包括瞬变电磁感应电压信号和环境噪声信号经减法器去除工频及其谐波估计值得到去噪后的瞬变电磁信号。本方法可去除任意频率的工频及其谐波干扰。

Description

一种瞬变电磁信号工频及其谐波干扰消除方法及装置

技术领域

本发明属于电磁信号处理，具体地来讲为一种瞬变电磁信号工频及其谐波干扰消除方法及装置。

背景技术

瞬变电磁探测方法是一种基于电磁感应原理的地球物理探测方法，具有施工效率高、分辨能力强优势，广泛应用于矿产勘探、油气探测等领域。在瞬变电磁地球物理勘探中，经常引入各种干扰和噪声，而来自高压输送电网，供电变压器接地方式等原因产生的工频及其谐波干扰是十分普遍的，尤其在城市附近，工频及其谐波干扰十分强烈。

瞬变电磁探测系统接收线圈在接收瞬变电磁信号的同时，不可避免的接收到工频及其谐波噪声，引起瞬变电磁信号信噪比降低，导致探测精度下降。并且瞬变电磁信号随时间成指数衰减，晚期道信号十分微弱，工频及其谐波干扰严重影响晚期道数据质量，导致地球物理深部信息反演结果可信度降低，勘探深度减小。

CN101548885A公开了一种去除电生理信号中工频干扰的方法。该发明通过将采集到的生理信号中减去重构干扰信号的方法，得到纯净的电生理信号。由于电生理信号主要集中在30Hz以下，50Hz附近的能量非常小，该方法在重构过程中忽略了50Hz附近的信号，但是瞬变电磁为宽频信号，不能忽略50Hz附近的能量，因此该方法不适用于瞬变电磁信号去除工频干扰。

CN101666834A公开了一种抗工频干扰的信号采集方法及系统。该方法首先通过从输入的市电交流信号中提取工频信号，并对其进行同步跟踪，产生工频同步脉冲。其次根据工频同步信号锁定检测信号中工频干扰的过零点相位，在过零点相位对输入的检测信号进行同步采样。该方法虽然在实现消除工频干扰的过程中，未损失有用信号，但过零点相位检测时，过零点相位容易产生偏差，影响测量结果。

CN104702243A公开了一种基于模糊逻辑滤除工频干扰的自适应滤波器。通过数字信号的模糊逻辑分析，构建匹配的滤波器，对信号进行滤波，实现生物信号工频干扰的消除。该方法虽然使工频干扰得到衰减，但有用信号也遭到同样的衰减，导致测量保真度下降。

发明内容

本发明所要解决的技术问题在于提供一种瞬变电磁信号工频及其谐波干扰消除方法及装置，在不损失有用信号的前提下，去除工频及其谐波干扰对瞬变电磁数据的影响，以提高瞬变电磁信号的信噪比，尤其对于晚期道信号的提取，提高深部反演信息的准确性。

本发明是这样实现的，

一种瞬变电磁信号工频及其谐波干扰消除方法，该方法包括，

步骤1，在第一时间段内接收环境噪声信号N(t)，经放大后，经低通滤波去除高频噪声，对其信号中的工频及其谐波噪声进行建模得到工频及其谐波干扰模型，

步骤2，从第二时间段开始，打开发射机，接收包括瞬变电磁感应电压信号和环境噪声信号，采用所述第一时间段内的工频及其谐波干扰模型计算所述第二时间段内及后续时间段的工频及其谐波估计值，

步骤3，将包括瞬变电磁感应电压信号和环境噪声信号经减法器去除工频及其谐波估计值得到去噪后的瞬变电磁信号。

方法进一步地，步骤1包括：

采用傅里叶级数拟合，对工频及其谐波噪声建立模型：其中k为采样时间，A_m为m次谐波分量对应的幅度，φ_m为m次谐波分量对应的相位，f₀为工频噪声频率，f_s为采样频率，M取小于等于200，N_p(k)表示工频及其谐波干扰模型函数，对式(1)进行正交分解，写成：

$N_p\left(k\right)=\underset{m=1}{\overset{M}{\Σ}}{\mathrm{\α}}_{m}cos\left(2\mathrm{\π}m\frac{f_0}{f_s}k\right)+\underset{m=1}{\overset{M}{\Σ}}{\mathrm{\β}}_{m}sin\left(2\mathrm{\π}m\frac{f_0}{f_s}k\right)---\left(2\right)$

其中

求出α_m、β_m及工频噪声频率f₀得到工频及其谐波干扰随时间变化的函数N_p(k)。

方法进一步地，根据工频噪声频率f₀的波动区间，采用二分法进行迭代计算，逐步减少搜索空间求取工频噪声频率f₀以及α_m和β_m。

方法进一步地，根据工频噪声频率f₀的波动区间，采用二分法进行迭代计算，逐步减少搜索空间求取工频噪声频率f₀以及α_m和β_m包括：

式(2)写成：

$\left[\begin{array}{ccccc}\cos (2{\mathrm{\πf}}_0/f_s)& \sin (2{\mathrm{\πf}}_0/f_s)& ...& \cos (2{\mathrm{\πmf}}_0/f_s)& \sin (2{\mathrm{\πmf}}_0/f_s)\\ \cos (4{\mathrm{\πf}}_0/f_s)& \sin (4{\mathrm{\πf}}_0/f_s)& ...& \cos (4{\mathrm{\πmf}}_0/f_s)& \sin (4{\mathrm{\πmf}}_0/f_s)\\ \·& \·& & \·& \·\\ \·& \·& ...& \·& \·\\ \·& \·& & \·& \·\\ \cos \left(2\mathrm{\π}\right(k-1)f_0/f_s)& \sin \left(2\mathrm{\π}\right(k-1)f_0/f_s)& ...& \cos \left(2\mathrm{\π}m\right(k-1)f_0/f_s)& \sin \left(2\mathrm{\π}m\right(k-1)f_0/f_s)\\ \cos (2{\mathrm{\πkf}}_0/f_s)& \sin (2{\mathrm{\πkf}}_0/f_s)& ...& \cos (2{\mathrm{\πmkf}}_0/f_s)& \sin (2{\mathrm{\πkf}}_0/f_s)\end{array}\right]\left[\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\β}}_1\\ \·\\ \·\\ \·\\ {\mathrm{\α}}_m\\ {\mathrm{\β}}_m\end{array}\right]=\left[\begin{array}{c}{Y}_1\\ Y_2\\ \·\\ \·\\ \·\\ Y_{k-1}\\ Y_k\end{array}\right]---\left(3\right)$

求解方程组(3)，得到系数X＝[α₁β₁…α_mβ_m]^T，

区间中点其中f_0a′表示波动区间的上端点、f_0b′表示波动区间的下端点，f₀′表示区间中点，分别计算波动区间的上端点f_0a′、波动区间的下端点f_0b′以及区间中点f₀′对应的工频及其谐波干扰模型，分别表示为：N′_pa(k)、N′_pb(k)、和N′_p0(k)；

将得到的f_0a′表示波动区间的上端点、f_0b′表示波动区间的下端点，f₀′表示区间中点的工频噪声分别与环境噪声N(t)相减得到剩余的随机噪声，并计算剩余随机噪声的均方差，分别表示为：σ_0a、σ_0b、和σ₀；

判断若σ_0a-σ₀＞σ_0b-σ₀，则工频在[f₀′_,f_0b′]区间，令f_0a′＝f₀′后，继续迭代，直至找到剩余随机噪声的均方差最小值点，输出工频f₀′及对应系数矩阵X′₀；否则令f_0b′＝f₀′，继续迭代，直至找到随机噪声的均方差最小值点，此时f₀＝f₀′，X₀＝X′₀输出工频f₀及对应系数矩阵X₀；

确定α_m、β_m及f₀后，将其带入式(1)，得到工频及其谐波干扰模型。

一种瞬变电磁信号工频及其谐波干扰消除装置，该装置包括，接收线圈，接收第一时间的环境噪声信号N(t)以及第二时间的包括瞬变电磁感应电压信号和环境噪声；

前置放大器，接收的信号进行放大处理；

低通滤波器，将放大的信号去除高频噪声；

减法器，其第一输入接收去除高频噪声的信号；

处理器，存储通过第一时间计算的工频及其谐波干扰模型，根据第一时间计算的工频及其谐波噪声模型，计算出第二时间工频及其谐波噪声估计值，噪声经DA转换，转换为模拟量，至所述减法器的第二输入；

减法器将两路输入经减法计算后输出。

装置进一步地，所述处理器在所述第一时间内，采用傅里叶级数拟合，对工频及其谐波噪声进行建立包含工频噪声信号频率的模型：其中k为采样时间，A_m为m次谐波分量对应的幅度，φ_m为m次谐波分量对应的相位，f₀为工频噪声频率，f_s为采样频率，M取小于等于200，N_p(k)表示工频及其谐波干扰模型函数，对式(1)进行正交分解，写成：

$N_p\left(k\right)=\underset{m=1}{\overset{M}{\Σ}}{\mathrm{\α}}_{m}cos\left(2\mathrm{\π}m\frac{f_0}{f_s}k\right)+\underset{m=1}{\overset{M}{\Σ}}{\mathrm{\β}}_{m}sin\left(2\mathrm{\π}m\frac{f_0}{f_s}k\right)---\left(2\right)$

其中 $A_m=\sqrt{{{\mathrm{\α}}_m}^2+{{\mathrm{\β}}_m}^2},{\mathrm{tan\φ}}_m=-\frac{{\mathrm{\β}}_m}{{\mathrm{\α}}_m}.$

求出α_m、β_m及工频噪声频率f₀得到工频及其谐波干扰随时间变化的函数N_p(k)。

装置进一步地，所述处理器，根据工频噪声频率f₀的波动区间，采用二分法进行迭代计算，逐步减少搜索空间求取工频噪声频率f₀以及α_m和β_m。

装置进一步地，所述处理器，根据工频噪声频率f₀的波动区间，采用二分法进行迭代计算，逐步减少搜索空间求取工频噪声频率f₀以及α_m和β_m包括：

将式(2)写成：

$\left[\begin{array}{ccccc}\cos (2{\mathrm{\πf}}_0/f_s)& \sin (2{\mathrm{\πf}}_0/f_s)& ...& \cos (2{\mathrm{\πmf}}_0/f_s)& \sin (2{\mathrm{\πmf}}_0/f_s)\\ \cos (4{\mathrm{\πf}}_0/f_s)& \sin (4{\mathrm{\πf}}_0/f_s)& ...& \cos (4{\mathrm{\πmf}}_0/f_s)& \sin (4{\mathrm{\πmf}}_0/f_s)\\ \·& \·& & \·& \·\\ \·& \·& ...& \·& \·\\ \·& \·& & \·& \·\\ \cos \left(2\mathrm{\π}\right(k-1)f_0/f_s)& \sin \left(2\mathrm{\π}\right(k-1)f_0/f_s)& ...& \cos \left(2\mathrm{\π}m\right(k-1)f_0/f_s)& \sin \left(2\mathrm{\π}m\right(k-1)f_0/f_s)\\ \cos (2{\mathrm{\πkf}}_0/f_s)& \sin (2{\mathrm{\πkf}}_0/f_s)& ...& \cos (2{\mathrm{\πmkf}}_0/f_s)& \sin (2{\mathrm{\πkf}}_0/f_s)\end{array}\right]\left[\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\β}}_1\\ \·\\ \·\\ \·\\ {\mathrm{\α}}_m\\ {\mathrm{\β}}_m\end{array}\right]=\left[\begin{array}{c}{Y}_1\\ Y_2\\ \·\\ \·\\ \·\\ Y_{k-1}\\ Y_k\end{array}\right]---\left(3\right)$

根据求解方程组(3)，得到系数X＝[α₁β₁…α_mβ_m]^T，

包括：计算区间中点其中f_0a′表示波动区间的上端点、f_0b′表示波动区间的下端点，f₀′表示区间中点，分别计算波动区间的上端点f_0a′、波动区间的下端点f_0b′以及区间中点f₀′对应的工频噪声，分别表示为：N′_pa(k)、N′_pb(k)、和N′_p0(k)；

将得到的f_0a′表示波动区间的上端点、f_0b′表示波动区间的下端点，f₀′表示区间中点的工频噪声分别与环境噪声N(t)相减得到剩余的随机噪声并计算剩余随机噪声的均方差，分别表示为：σ_0a、σ_0b、和σ₀；

判断若σ_0a-σ₀＞σ_0b-σ₀，则工频在[f₀′_,f_0b′]区间，令f_0a′＝f₀′后，继续迭代，直至找到剩余随机噪声的均方差最小值点，输出工频f₀′及对应系数矩阵X′₀；

否则令f_0b′＝f₀′，继续迭代，直至找到随机噪声的均方差最小值点，f₀＝f₀′，X₀＝X′₀输出工频f₀及对应系数矩阵X₀；

确定α_m、β_m及f₀后，将其带入式(1)，得到工频及其谐波干扰模型。

本发明与现有技术相比，有益效果在于：本发明公开的瞬变电磁信号工频及其谐波干扰实时消除方法及装置，通过直接对数据中的时域工频及其谐波噪声建立模型，可以有效的去除工频及其谐波干扰对瞬变电磁信号的影响，同时不损失有用信息。本方法可去除任意频率的工频及其谐波干扰。

经试验，本方法对瞬变电磁数据中工频及其谐波干扰的去除效率高，准确性好，提高了深部反演结果的可信度。

附图说明

图1为基于噪声模型实时去除工频及其谐波干扰原理框图；

图2为发射电流示意图；

图3为瞬变电磁中心回线工作示意图；

图4为工频及其谐波干扰噪声模型效果图，(a)为环境噪声信号波形图，(b)为工频及其谐波干扰模型波形图；

图5为利用噪声模型去除工频及其谐波干扰效果图；

图6为瞬变电磁数据实时去除工频及其谐波干扰的效果图。

具体实施方式

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

本发明实施例提供的基于噪声建模的瞬变电磁信号工频及其谐波干扰消除装置，原理框图如图1所示，本实施例采用基于中心回线的工作方式，如图3所示，发射线圈11铺设为方形回线，边长为100m，接收线圈1为边长为1m的方形回线。

接收线圈1接收环境噪声信号n(t)(图2中T0～T1时间段)，信号经前置放大器2放大后，通过RC低通滤波器3，去除高频噪声，输出环境噪声信号N(t)。信号N(t)经AD转换器4，将模拟量转化为数字量N(k)，输出信号如图4(a)所示。AD转换器4将数字量输出至处理器DSP5中，本实施例中，采用DSP，首先在DSP中建立环境噪声的工频及其谐波干扰模型，具体如下所述：

利用傅里叶级数拟合，对信号中的工频及其谐波噪声进行建模。

工频及其谐波干扰模型为：

$N_p\left(k\right)=\underset{m=1}{\overset{M}{\Σ}}{A}_{m}cos(2\mathrm{\π}m\frac{f_0}{f_s}k+{\mathrm{\φ}}_m)---\left(1\right)$

其中，k为采样时间，A_m为m次谐波分量对应的幅度，φ_m为m次谐波分量对应的相位，f₀为工频噪声频率，f_s为采样频率，M取200。

对式(1)进行正交分解，可以写成：

$N_p\left(k\right)=\underset{m=1}{\overset{M}{\Σ}}{\mathrm{\α}}_{m}cos\left(2\mathrm{\π}m\frac{f_0}{f_s}k\right)+\underset{m=1}{\overset{M}{\Σ}}{\mathrm{\β}}_{m}sin\left(2\mathrm{\π}m\frac{f_0}{f_s}k\right)---\left(2\right)$

其中 $A_m=\sqrt{{{\mathrm{\α}}_m}^2+{{\mathrm{\β}}_m}^2},{\mathrm{tan\φ}}_m=-\frac{{\mathrm{\β}}_m}{{\mathrm{\α}}_m}.$

求出α_m、β_m及f₀即得到工频及其谐波干扰随时间变化的函数N_p(k)。

由于工频噪声频率f₀通常在49Hz～51Hz之间波动，因此首先需要确定工频噪声频率f₀，具体通过二分法进行迭代计算，逐步减小搜索空间，步骤如下：

a，计算区间中点分别计算工频区间端点f_0a′、f_0b′及区间中点f₀′对应的工频及其谐波干扰模型。

当工频确定时，式(2)可写成：

$\left[\begin{array}{ccccc}\cos (2{\mathrm{\πf}}_0/f_s)& \sin (2{\mathrm{\πf}}_0/f_s)& ...& \cos (2{\mathrm{\πmf}}_0/f_s)& \sin (2{\mathrm{\πmf}}_0/f_s)\\ \cos (4{\mathrm{\πf}}_0/f_s)& \sin (4{\mathrm{\πf}}_0/f_s)& ...& \cos (4{\mathrm{\πmf}}_0/f_s)& \sin (4{\mathrm{\πmf}}_0/f_s)\\ \·& \·& & \·& \·\\ \·& \·& ...& \·& \·\\ \·& \·& & \·& \·\\ \cos \left(2\mathrm{\π}\right(k-1)f_0/f_s)& \sin \left(2\mathrm{\π}\right(k-1)f_0/f_s)& ...& \cos \left(2\mathrm{\π}m\right(k-1)f_0/f_s)& \sin \left(2\mathrm{\π}m\right(k-1)f_0/f_s)\\ \cos (2{\mathrm{\πkf}}_0/f_s)& \sin (2{\mathrm{\πkf}}_0/f_s)& ...& \cos (2{\mathrm{\πmkf}}_0/f_s)& \sin (2{\mathrm{\πkf}}_0/f_s)\end{array}\right]\left[\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\β}}_1\\ \·\\ \·\\ \·\\ {\mathrm{\α}}_m\\ {\mathrm{\β}}_m\end{array}\right]=\left[\begin{array}{c}{Y}_1\\ Y_2\\ \·\\ \·\\ \·\\ Y_{k-1}\\ Y_k\end{array}\right]---\left(3\right)$

求解方程组(3)，得到系数X＝[α₁β₁…α_mβ_m]^T。将f_0a′＝49Hz、f_0b′＝51Hz及f₀′＝50Hz，带入方程(3)，求解方程组，得到对应的系数矩阵X′_0a、X′_0b、X′₀，分别将其带入式(1)即得到工频及其谐波干扰模型N′_pa(k)、N′_pb(k)、N′_p0(k)。

b，将得到的工频及其谐波干扰模型与环境噪声N(k)相减，得到剩余的随机噪声N_r(k)，计算随机噪声N_r(k)的均方差：

$\mathrm{\σ}=\sqrt{\frac{1}{K}\underset{k=1}{\overset{K}{\Σ}}\left(N_r\right(k)-{\mathrm{\μ}}_N)}---\left(4\right)$

其中均值K为采样点个数。分别计算步骤b中得到的N′_pa(k)、N′_pb(k)、N′_p0(k)的剩余随机噪声的均方差，得到σ_0a＝7.003，σ₀＝5.157，σ_0b＝6.172。

c，经判断σ_0a-σ₀＞σ_0b-σ₀，则工频在[50,51]Hz区间，令f_0a′＝50Hz，转至步骤a，继续迭代，反复进行上述步骤，直至f₀′＝50.406Hz时，σ₀＝1.057取得最小值，即输出工频f₀＝50.406及对应系数矩阵X′₀。确定α_m、β_m及f₀后，将其带入式(1)，得到工频及其谐波干扰模型如图4(b)所示，并保存至DSP中。利用噪声模型去除工频及其谐波干扰后效果如图5所示。

第三步，第二个时间段内，打开瞬变电磁发射装置(图2中T1时刻)，接收线圈在off-time段(图2中T2～T3时间段)接收到信号包含感应电压信号v(t)＝s(t)+n(t)，如图6中黑色曲线所示，为双对数坐标系下显示；经前置放大器2放大后，通过RC低通滤波器3，去除高频噪声，输出v′(t)，至减法器7。

第四步，利用DSP5中存储的工频及其谐波干扰模型函数式(1)计算出T2～T3时间工频及其谐波干扰的估计值，噪声经DA转换器6，转换为模拟量N_p(t)，至减法器7。

第五步，信号v′(t)减去N_p(t)，得到去除工频及其谐波干扰的模拟瞬变电磁信号v_s(t)，经末端放大器8得到输出信号，双对数坐标系下如图6中灰色曲线所示，可以看出瞬变电磁数据经本发明的装置及方法去除工频及其谐波噪声后，信噪比明显提高，尤其对于反映地球物理深部信息的晚期数据。最后对去除工频干扰后的信号v_s(t)进行放大，得到满足数据采集系统要求的有效瞬变电磁信号。

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

Claims (8)

1.一种瞬变电磁信号工频及其谐波干扰消除方法，其特征在于，该方法包括，

步骤1，在第一时间段内接收环境噪声信号N(t)，经放大后，经低通滤波去除高频噪声，对其信号中的工频及其谐波噪声进行建模得到工频及其谐波干扰模型，

步骤2，从第二时间段开始，打开发射机，接收包括瞬变电磁感应电压信号和环境噪声信号，采用所述第一时间段内的工频及其谐波干扰模型计算所述第二时间段内及后续时间段的工频及其谐波估计值，

步骤3，将包括瞬变电磁感应电压信号和环境噪声信号经减法器去除工频及其谐波估计值得到去噪后的瞬变电磁信号。

2.按照权利要求1所述的瞬变电磁信号工频及其谐波干扰消除方法，其特征在于，步骤1包括：

采用傅里叶级数拟合，对工频及其谐波噪声建立模型：其中k为采样时间，A_m为m次谐波分量对应的幅度，φ_m为m次谐波分量对应的相位，f₀为工频噪声信号频率，f_s为采样频率，M取小于等于200，N_p(k)表示工频及其谐波干扰模型函数，对式(1)进行正交分解，写成：

$N_p\left(k\right)=\underset{m=1}{\overset{M}{\Σ}}{\mathrm{\α}}_{m}cos\left(2\mathrm{\π}m\frac{f_0}{f_s}k\right)+\underset{m=1}{\overset{M}{\Σ}}{\mathrm{\β}}_{m}sin\left(2\mathrm{\π}m\frac{f_0}{f_s}k\right)---\left(2\right)$

其中 $A_m=\sqrt{{{\mathrm{\α}}_m}^2+{{\mathrm{\β}}_m}^2},{\mathrm{tan\φ}}_m=-\frac{{\mathrm{\β}}_m}{{\mathrm{\α}}_m}.$

求出α_m、β_m及工频噪声频率f₀得到工频及其谐波干扰随时间变化的函数N_p(k)。

3.按照权利要求2所述的瞬变电磁信号工频及其谐波干扰消除方法，其特征在于，根据工频噪声频率f₀的波动区间，采用二分法进行迭代计算，逐步减少搜索空间求取工频噪声频率f₀以及α_m和β_m。

4.按照权利要求3所述的瞬变电磁信号工频及其谐波干扰消除方法，其特征在于，根据工频噪声频率f₀的波动区间，采用二分法进行迭代计算，逐步减少搜索空间求取工频噪声频率f₀以及α_m和β_m包括：

式(2)写成：

$\left[\begin{array}{ccccc}\cos \left(2{\mathrm{\πf}}_0/f_s\right)& \sin \left(2{\mathrm{\πf}}_0/f_s\right)& ...& \cos \left(2\mathrm{\π}m{f}_0/f_s\right)& \sin \left(2\mathrm{\π}m{f}_0/f_s\right)\\ \cos \left(4{\mathrm{\πf}}_0/f_s\right)& \sin \left(4{\mathrm{\πf}}_0/f_s\right)& ...& \cos \left(4\mathrm{\π}m{f}_0/f_s\right)& \sin \left(4\mathrm{\π}m{f}_0/f_s\right)\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ ...\\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ \cos \left(2\mathrm{\π}\left(k-1\right)f_0/f_s\right)& \sin \left(2\mathrm{\π}\left(k-1\right)f_0/f_s\right)& ...& \cos \left(2\mathrm{\π}m\left(k-1\right)f_0/f_s\right)& \sin \left(2\mathrm{\π}m\left(k-1\right)f_0/f_s\right)\\ \cos \left(2\mathrm{\π}k{f}_0/f_s\right)& \sin \left(2\mathrm{\π}k{f}_0/f_s\right)& ...& \cos \left(2\mathrm{\π}mk{f}_0/f_s\right)& \sin \left(2{\mathrm{\πkf}}_0/f_s\right)\end{array}\right]\left[\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\β}}_1\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{\α}}_m\\ {\mathrm{\β}}_m\end{array}\right]=\left[\begin{array}{c}{Y}_1\\ Y_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ Y_{k-1}\\ Y_k\end{array}\right]---\left(3\right)$

求解方程组(3)，得到系数X＝[α₁β₁…α_mβ_m]^T，

区间中点其中f_0a′表示波动区间的上端点、f_0b′表示波动区间的下端点，f′₀表示区间中点，分别计算波动区间的上端点f_0a′、波动区间的下端点f_0b′以及区间中点f′₀对应的工频及其谐波干扰模型，分别表示为：N′_pa(k)、N′_pb(k)、和N′_p0(k)；

将得到的f_0a′表示波动区间的上端点、f_0b′表示波动区间的下端点，f′₀表示区间中点的工频噪声分别与环境噪声N(t)相减得到剩余的随机噪声，并计算剩余随机噪声的均方差，分别表示为：σ_0a、σ_0b、和σ₀；

判断若σ_0a-σ₀＞σ_0b-σ₀，则工频在[f′₀,f_0b′]区间，令f_0a′＝f′₀后，继续迭代，直至找到剩余随机噪声的均方差最小值点，输出工频f′₀及对应系数矩阵x′₀；否则令f_0b′＝f′₀，继续迭代，直至找到随机噪声的均方差最小值点，此时f₀＝f′₀，x₀＝x′₀输出工频f₀及对应系数矩阵x₀；

确定α_m、β_m及f₀后，将其带入式(1)，得到工频及其谐波干扰模型。

5.一种瞬变电磁信号工频及其谐波干扰消除装置，其特征在于，该装置包括，接收线圈，接收第一时间的环境噪声信号N(t)以及第二时间的包括瞬变电磁感应电压信号和环境噪声；

前置放大器，接收的信号进行放大处理；

低通滤波器，将放大的信号去除高频噪声；

减法器，其第一输入接收去除高频噪声的信号；

处理器，存储通过第一时间计算的工频及其谐波干扰模型，根据第一时间计算的工频及其谐波干扰模型，计算出第二时间工频及其谐波估计值，噪声经DA转换，转换为模拟量，至所述减法器的第二输入；

减法器将两路输入经减法计算后输出。

6.按照权利要求5所述的瞬变电磁信号工频及其谐波干扰消除装置，其特征在于，所述处理器在所述第一时间内，采用傅里叶级数拟合，对工频及其谐波噪声进行建立包含工频噪声频率的模型：其中k为采样时间，A_m为m次谐波分量对应的幅度，φ_m为m次谐波分量对应的相位，f₀为工频噪声频率，f_s为采样频率，M取小于等于200，N_p(k)表示工频及其谐波干扰模型函数，对式(1)进行正交分解，写成：

$N_p\left(k\right)=\underset{m=1}{\overset{M}{\Σ}}{\mathrm{\α}}_m\cos \left(2\mathrm{\π}m\frac{f_0}{f_s}k\right)+\underset{m=1}{\overset{M}{\Σ}}{\mathrm{\β}}_m\sin \left(2\mathrm{\π}m\frac{f_0}{f_s}k\right)---\left(2\right)$

其中 $A_m=\sqrt{{{\mathrm{\α}}_m}^2+{{\mathrm{\β}}_m}^2},{\mathrm{tan\φ}}_m=-\frac{{\mathrm{\β}}_m}{{\mathrm{\α}}_m}.$

求出α_m、β_m及工频噪声频率f₀得到工频及其谐波干扰随时间变化的函数N_p(k)。

7.按照权利要求6所述的瞬变电磁信号工频及其谐波干扰消除装置，其特征在于，所述处理器，根据工频噪声频率f₀的波动区间，采用二分法进行迭代计算，逐步减少搜索空间求取工频噪声频率f₀以及α_m和β_m。

8.按照权利要求6所述的瞬变电磁信号工频及其谐波干扰消除装置，其特征在于，

所述处理器，根据工频噪声频率f₀的波动区间，采用二分法进行迭代计算，逐步减少搜索空间求取工频噪声频率f₀以及α_m和β_m包括：

将式(2)写成：

$\left[\begin{array}{ccccc}\cos \left(2{\mathrm{\πf}}_0/f_s\right)& \sin \left(2{\mathrm{\πf}}_0/f_s\right)& ...& \cos \left(2\mathrm{\π}m{f}_0/f_s\right)& \sin \left(2\mathrm{\π}m{f}_0/f_s\right)\\ \cos \left(4{\mathrm{\πf}}_0/f_s\right)& \sin \left(4{\mathrm{\πf}}_0/f_s\right)& ...& \cos \left(4\mathrm{\π}m{f}_0/f_s\right)& \sin \left(4\mathrm{\π}m{f}_0/f_s\right)\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ ...\\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ \cos \left(2\mathrm{\π}\left(k-1\right)f_0/f_s\right)& \sin \left(2\mathrm{\π}\left(k-1\right)f_0/f_s\right)& ...& \cos \left(2\mathrm{\π}m\left(k-1\right)f_0/f_s\right)& \sin \left(2\mathrm{\π}m\left(k-1\right)f_0/f_s\right)\\ \cos \left(2\mathrm{\π}k{f}_0/f_s\right)& \sin \left(2\mathrm{\π}k{f}_0/f_s\right)& ...& \cos \left(2\mathrm{\π}mk{f}_0/f_s\right)& \sin \left(2{\mathrm{\πkf}}_0/f_s\right)\end{array}\right]\left[\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\β}}_1\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{\α}}_m\\ {\mathrm{\β}}_m\end{array}\right]=\left[\begin{array}{c}{Y}_1\\ Y_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ Y_{k-1}\\ Y_k\end{array}\right]---\left(3\right)$

根据求解方程组(3)，得到系数X＝[α₁β₁…α_mβ_m]^T，

包括：计算区间中点其中f_0a′表示波动区间的上端点、f_0b′表示波动区间的下端点，f′₀表示区间中点，分别计算波动区间的上端点f_0a′、波动区间的下端点f_0b′以及区间中点f′₀对应的工频噪声，分别表示为：N′_pa(k)、N′_pb(k)、和N′_p0(k)；

将得到的f_0a′表示波动区间的上端点、f_0b′表示波动区间的下端点，f′₀表示区间中点的工频噪声分别与环境噪声N(t)相减得到剩余的随机噪声并计算剩余随机噪声的均方差，分别表示为：σ_0a、σ_0b、和σ₀；

判断若σ_0a-σ₀＞σ_0b-σ₀，则工频在[f′₀,f_0b′]区间，令f_0a′＝f′₀后，继续迭代，直至找到剩余随机噪声的均方差最小值点，输出工频f′₀及对应系数矩阵x′₀；

否则令f_0b′＝f′₀，继续迭代，直至找到随机噪声的均方差最小值点，f₀＝f′₀，x₀＝x′₀输出工频f₀及对应系数矩阵x₀；

确定α_m、β_m及f₀后，将其带入式(1)，得到工频及其谐波干扰模型。