CN108345039A - 一种消除地面核磁共振数据中邻频谐波干扰的方法 - Google Patents

Abstract

本发明属于地面核磁共振数据处理领域，为一种消除地面核磁共振数据中邻频谐波的方法。该方法包括：录入地面核磁共振测量数据，测量实验地点地磁场，计算拉莫尔频率；搜索工频谐波基频，并对工频谐波进行建模，计算得到所有谐波系数；根据拉莫尔频率和每个谐波频率的频率差判断是否存在邻频谐波，若否，利用谐波基频和谐波系数计算得到谐波模型，在测量数据中减去谐波模型，若是，则去除领域谐波外的工频谐波后基于拟合方法去除领频谐波，本发明在不损失核磁共振信号的前提下，去除工频谐波中的邻频谐波对后续数据处理和反演的影响。尤其当实验地点的拉莫尔频率与工频谐波的高阶谐波相近时，可准确提取地面核磁共振信号特征，提高反演结果准确度。

Description

一种消除地面核磁共振数据中邻频谐波干扰的方法

技术领域

本发明属于地面核磁共振数据处理领域，具体地来讲为一种消除地面核磁共振数据中邻频谐波的方法。

背景技术

地面核磁共振是一种直接探测地下水赋存属性的地球物理方法，具有测量速度快，分辨率高和定量解释等优点，被广泛应用于水资源勘探、水文环境监测以及地下工程水害预警等领域。但是地面核磁共振方法获得的信号十分微弱，通常强度只有纳伏量级(nV，10^-9V)，易受到各种环境噪声的干扰，导致采集数据信噪比低，严重影响后续反演结果的准确性，制约了该方法在城镇、村庄以及地下工程等复杂环境中的应用。

地面核磁共振数据中主要包含三类噪声：尖峰噪声、工频及其谐波噪声和随机噪声，其中高压输送电网和电气设备等产生的工频谐波噪声最为强烈。对于工频谐波噪声，当存在其中一个谐波的频率与地面核磁共振信号的频率(称为拉莫尔频率)相近或相等时(小于等于3Hz)，我们称之为邻频谐波。

专利CN103913778A公开了一种多个近端参考线圈的核磁共振信号实时噪声抵消装置。该方法首先消除参考线圈内的地面核磁共振信号，并利用参考通道与主通道的频率相关性，实现工频谐波的消除。但是当存在邻频谐波时，相关系数无法正确计算，且铺设参考线圈增加了野外实验的复杂度。

专利CN104898172A公开了一种互相关的核磁共振全波信号噪声滤除方法。该方法通过构造一个与地面核磁共振信号频率相同的正弦信号作为参考信号，并基于随机噪声和尖峰噪声与参考信号相关性较弱，而地面核磁共振信号与参考线圈相关性较强的原则，实现地面核磁共振数据的噪声滤除。但是当存在邻频谐波时，该方法无法有效区分谐波噪声和地面核磁共振信号。

专利CN106772646A公开了一种地面核磁共振信号提取方法。对于工频谐波噪声，提出了利用谐波建模消噪方法，该方法无需额外的参考通道，仅通过采集数据建立工频谐波模型。存在邻频谐波噪声时，该方法无法同样无法区分谐波噪声和地面核磁共振信号，导致建模结果产生错误。

发明内容

本发明所要解决的技术问题在于提供一种消除地面核磁共振数据中工频及邻频谐波干扰的方法，在不损失核磁共振信号的前提下，去除工频谐波中的邻频谐波对后续数据处理和反演的影响。尤其当实验地点的拉莫尔频率与工频谐波的高阶谐波相近时，可准确提取地面核磁共振信号特征，提高反演结果准确度。

本发明是这样实现的，

一种消除地面核磁共振数据中邻频谐波干扰的方法，该方法包括：

步骤1，录入地面核磁共振测量数据V＝[V₁,V₂,L,V_K]，K为采样点数，测量实验地点地磁场B₀，根据式f_L＝0.04258×B₀计算拉莫尔频率f_L，

步骤2，搜索工频谐波基频f₀，并对工频谐波进行建模，计算得到所有谐波系数X；

步骤3，根据拉莫尔频率f_L和每个谐波频率f_m＝m·f₀的频率差df_m＝f_L-f_m判断是否存在邻频谐波，其中m＝1,2,L,M，M是谐波个数。若所有的频率差|df_m|≥3Hz，说明不存在邻频谐波，利用步骤2得到的谐波基频f_o和谐波系数X计算得到谐波模型V_harm，在测量数据V中减去谐波模型V_harm，实现工频及其谐波干扰的消除。

进一步地，步骤3中，若其中一个频率差|df_m|＜3Hz，则存在邻频谐波，进行步骤4；

步骤4，将谐波系数X中邻频谐波对应的系数置零，进行谐波建模，得到除邻频谐波外的谐波模型V′_harm，在测量数据V中减去谐波模型V′_harm，首先完成除邻频谐波外所有工频谐波的消除，得到数据V′，基于拟合的方法，得到邻频谐波V_m的参数，在数据V′中减去邻频谐波V_m，得到去除谐波噪声后的地面核磁共振信号。

进一步地，步骤4中基于拟合的方法包括对数据V′以邻频谐波的频率f_m进行希尔伯特变换，并经过低通滤波器后得到数据V′包络的实部V_Re和虚部V_Im；

进一步地，对数据V′包络的实部V_Re和虚部V_Im进行非线性拟合获得6个未知参数x₁～x₆，其中四个参数对应于地面核磁共振信号参数，两个参数为邻频谐波的参数，获得邻频谐波的幅度和相位，计算得到邻频谐波模型V_m，在数据V′中减去邻频谐波模型V_m完成邻频谐波的消除。

经过低通滤波器后得到V′包络的实部V_Re和虚部V_Im分别为：

其中e₀、和θ分别为地面核磁共振信号的初始幅度，弛豫时间和相位，df_m＝f_L-f_m是拉莫尔频率f_L与邻波谐波频率f_m的频率差。A_m和分别为邻频谐波的幅度和相位。ε_r和ε_i分别是随机噪声的实部和虚部分量。

进一步地，根据表达式(1)和式(2)对数据V′包络的实部V_Re和虚部V_Im进行非线性拟合获得6个未知参数x₁～x₆，

其中，分别为经拟合得到的实部和虚部，x₁～x₄对应于地面核磁共振信号参数，x₅和x₆为邻频谐波的参数。

进一步地，步骤2包括：利用搜索得到的工频谐波基频f₀，根据式(3)对工频及其谐波噪声建立模型V_harm：

其中，A_m和分别是第m个谐波的幅度和相位，M是谐波个数，t为采样时间，t＝[t₁,t₂,L,t_K]，对式(3)中的第m项进行线性化为：

其中，α_m和β_m为第m阶谐波对应系数，和求解线性方程组，

得到谐波系数X＝[α₁,L,α_M,β₁L,β_M]^T，其中T是转置符号。

本发明与现有技术相比，有益效果在于：本发明公开的去除地面核磁共振测量数据中邻频谐波的方法，不需要参考线圈额外获取噪声数据。根据地面核磁共振信号和邻频谐波的幅度随时间变化的特征不同，利用非线性拟合方法同时估计两者的特征参数，进而实现对邻频谐波的建模和去除，且不损失有用信息。使得地面核磁共振方法可应用在拉莫尔频率和工频谐波的高阶谐波相近甚至重合的地点，拓宽了地面核磁共振方法的应用范围。

经试验，本发明方法对地面核磁共振数据中邻频谐波干扰的去除效率高，准确性好，提高了拉莫尔频率和工频谐波的高阶谐波相近的地点的探测结果的可信度。

附图说明

图1为地面核磁共振数据邻频谐波干扰去除方法的流程图；

图2地面核磁共振数据波形图和消除除邻频谐波外的其他谐波噪声的结果，(a)时域波形(b)频域波形；

图3希尔伯特变换和低通滤波后，地面核磁共振数据包络及拟合结果(a)实部(b)虚部；

图4地面核磁共振数据邻频谐波处理效果图(a)时域波形(b)频域波形。

具体实施方式

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

以一组初始振幅为50nV，弛豫时间为0.3s，拉莫尔频率为2300Hz，初始相位为π/3rad的地面核磁共振信号；工频谐波基频为50.010Hz，100个谐波的幅度和相位呈随机分布；随机噪声为10nV的数据为例进行说明。

参见图1，步骤1，录入地面核磁共振数据V＝[V₁,V₂,L,V_K]，K为采样点数，其时域和频域图见图2(a)和图2(b)黑色曲线所示。测量实验地点地磁场B₀为54016nT，计算拉莫尔频率f_L＝0.04258×B₀＝2300Hz。

步骤2，利用搜索得到的工频谐波基频f₀＝50.010Hz，根据式(3)对工频及其谐波噪声建立模型：

其中，A_m和分别是第m个谐波的幅度和相位，M是谐波个数，t为采样时间，t＝[t₁,t₂,L,t_K]。对式(3)中的第m项进行线性化为：

其中，和将f₀＝50.010Hz带入式(4)，求解线性方程组

得到谐波系数X＝[α₁,L,α_M,β₁L,β_M]^T，其中T是转置符号。

步骤3，根据拉莫尔频率f_L和每个谐波频率f_m＝m·f₀的频率差df_m＝f_L-f_m判断是否存在邻频谐波，其中第46阶谐波的频率为f₄₆＝n₄₆×f₀＝2300.46Hz，|df₄₆|＝0.46Hz，小于3Hz，即存在邻频谐波。

步骤4，将谐波系数X中邻频谐波对应的系数置零，即α₄₆＝0，β₄₆＝0。带入式(3)得到V′_harm，在测量数据V中减去V′_harm得到数据V′，即实现除邻频外的工频及其谐波噪声消除，数据V′的时域和频域结果见图2中圆圈曲线。

对数据V′以邻频谐波频率f₅₆＝2300.46Hz进行希尔伯特变换，并经过低通滤波后，得到数据V′包络的实部和虚部(见图3(a)和图3(b)灰色曲线)分别为：

其中，e₀、和θ分别为地面核磁共振信号的初始幅度，弛豫时间和相位，df₄₆＝f_L-f₄₆＝-0.46Hz是拉莫尔频率f_L与邻波谐波频率f₅₆的频率差。A₅₆和分别为邻频谐波的幅度和相位。ε_r和ε_i分别是随机噪声的实部和虚部分量。

根据表达式(1)和(2对测量数据的实部V_Re和虚部V_Im进行非线性拟合(拟合结果见图3中黑色圆圈曲线)，

得到：x₁＝50.15(nV)，x₂＝0.305(s)，x₃＝-0.462(Hz)，x₄＝1.048(rad)，对应的地面核磁共振信号见图4(a)和图4(b)中黑色方框曲线，x₅＝43.35(nV)和x₆＝24.82(nV)，对应的邻频谐波见图3中黑色虚线。得到邻频谐波的幅度和相位分别为：A₅₆＝49.95nV和带入式(3)得到邻频谐波模型V₅₆，在数据V′中减去V₅₆即完成了邻频谐波的消除，时域和频域结果见图4(a)和图4(b)中黑色曲线。

结果显示本发明方法对地面核磁共振数据中邻频谐波干扰的去除效率高，准确性好，提高了拉莫尔频率和工频谐波的高阶谐波相近的地点的探测结果的可信度。

Claims (6)

1.一种消除地面核磁共振数据中邻频谐波干扰的方法，其特征在于，该方法包括：

步骤1，录入地面核磁共振测量数据V＝[V₁,V₂,L,V_K]，K为采样点数，测量实验地点地磁场B₀，根据式f_L＝0.04258×B₀计算拉莫尔频率f_L；

步骤2，搜索工频谐波基频f₀，并对工频谐波进行建模，计算得到所有谐波系数X；

步骤3，根据拉莫尔频率f_L和每个谐波频率f_m＝m·f₀(m＝1,2,L,M，M是谐波个数)的频率差df_m＝f_L-f_m判断是否存在邻频谐波，若所有的频率差|df_m|≥3Hz，说明不存在邻频谐波，利用步骤2得到的谐波基频f_o和谐波系数X计算得到谐波模型V_harm，在测量数据V中减去谐波模型V_harm，实现工频及其谐波干扰的消除。

2.按照权利要求1所述的方法，其特征在于，步骤3中，若其中一个频率差|df_m|＜3Hz，则存在邻频谐波，进行步骤4；

步骤4，将谐波系数X中邻频谐波对应的系数置零，进行谐波建模，得到除邻频谐波外的谐波模型V_h′_arm，在测量数据V中减去谐波模型V_h′_arm，首先完成除邻频谐波外所有工频谐波的消除，得到数据V′，基于拟合的方法，得到邻频谐波V_m的参数，在数据V′中减去邻频谐波V_m，得到去除谐波噪声后的地面核磁共振信号。

3.按照权利要求2所述的方法，其特征在于，步骤4中基于拟合的方法包括对数据V′以邻频谐波的频率f_m进行希尔伯特变换，并经过低通滤波器后得到数据V′包络的实部V_Re和虚部V_Im；

对数据V′包络的实部V_Re和虚部V_Im进行非线性拟合获得6个未知参数x₁～x₆，其中四个参数对应于地面核磁共振信号参数，两个参数为邻频谐波的参数，获得邻频谐波的幅度和相位，计算得到邻频谐波模型V_m，在数据V′中减去邻频谐波模型V_m完成邻频谐波的消除。

4.按照权利要求3所述的方法，其特征在于，经过低通滤波器后得到V′包络的实部V_Re和虚部V_Im分别为：

其中e₀、和θ分别为地面核磁共振信号的初始幅度，弛豫时间和相位，

df_m＝f_L-f_m是拉莫尔频率f_L与邻波谐波频率f_m的频率差。A_m和分别为邻频谐波的幅度和相位。ε_r和ε_i分别是随机噪声的实部和虚部分量。

5.按照权利要求4所述的方法，其特征在于，根据表达式(1)和式(2)对数据V′包络的实部V_Re和虚部V_Im进行非线性拟合获得6个未知参数x₁～x₆，

其中，分别为经拟合得到的实部和虚部，x₁～x₄对应于地面核磁共振信号参数，x₅和x₆为邻频谐波的参数。

6.按照权利要求2所述的方法，其特征在于，步骤2包括：利用搜索得到的工频谐波基频f₀，根据式(3)对工频及其谐波噪声建立模型V_harm：

其中，A_m和分别是第m个谐波的幅度和相位，M是谐波个数，t为采样时间，t＝[t₁,t₂,L,t_K]，对式(3)中的第m项进行线性化为：

其中，α_m和β_m为第m阶谐波对应系数，和求解线性方程组，

得到谐波系数X＝[α₁,L,α_M,β₁L,β_M]^T，其中T是转置符号。