CN107957566B - 基于频率选择奇异谱分析的磁共振测深信号提取方法 - Google Patents

Abstract

本发明为一种基于频率选择的奇异谱分析的磁共振测深信号提取方法。首先利用核磁共振测深探水仪采集Larmor频率已知地区的MRS信号，通过宽频带带通滤波器进行部分噪声抑制后，基于功率谱分析在功率谱上找到Larmor频率对应的MRS信号的位置；然后，进行基于频率选择的奇异谱分析从而提取MRS信号。基于频率选择的奇异谱分析包括嵌入、RSVD分解、根据MRS信号幅值选择对应的奇异值进行矩阵重构和对角平均化四个步骤。本发明能够解决含噪MRS信号中随机噪声、尖峰噪声和工频谐波干扰的有效滤除，实现复杂强噪声干扰下MRS信号的有效提取，与传统MRS信号去噪方法相比，具有运算速度快、信噪比高、实用性强等优点。

Description

基于频率选择奇异谱分析的磁共振测深信号提取方法

技术领域

本发明涉及磁共振测深(又称地面核磁共振)地下水探测信号噪声滤除和参数提取技术领域，具体是利用基于频率选择奇异谱分析的磁共振测深信号提取方法。

背景技术

核磁共振测深(Magnetic Resonance Sounding,MRS)技术是唯一一种直接定量表征含水量和孔隙结构的方法，其基本原理是通过探测地下水中氢质子共振跃迁产生的MRS信号响应进行地下水探测。MRS信号的数量级是纳伏，十分微弱，极易被环境中的随机噪声、工频谐波和尖峰脉冲所干扰，导致探测到的MRS信号的质量受到影响，进而影响MRS信号特征参数的提取，降低了反演结果的准确性，造成了对区域水资源含量以及存储介质成分评定的不准确。

目前，针对磁共振测深信号噪声滤除和信号提取理论和方法研究方面，国内外专家学者开展了大量研究工作。Annette Hein等人在《Symmetry based frequency domainprocessing to remove harmonic noise from surface nuclear magnetic resonancemeasurements》(《Geophysical Journal International》2017年第208期：724-736页)提出了一种基于MRS信号在频域上对称的特点来消除工频噪声的方法，但是不能处理尖峰噪声和随机噪声；田宝凤等人在《Noise cancellation of a multi-reference full-wavemagnetic resonance sounding signal based on a modified sigmoid variable stepsize least mean square algorithm》(《Journal of Central South University ofTechnology》2017年第24期：900-911页)提出了一种基于改进s函数变步长LMS的多通道全波核磁共振信号参考对消方法，但是当主通道与参考通道相关性较差时，处理效果不好。张海如等人在《基于Bayes Bootstrap统计降噪方法的磁共振测深信号检测》(《中南大学学报(自然科学版)》2014年第45卷第9期：3144-3149页)中从各导电层MRS信号分量中提取2个最优估计点来重建MRS信号，但是处理后提取的信号参数幅值e₀和平均弛豫时间T₂ ^*的误差较大。

CN106772646A公开了“一种地面核磁共振信号提取方法”，该方法能自适应地找到基频为49.9Hz～50.1Hz的工频谐波，求得MRS信号的自相关表达式，并且快速有效地实现信号和噪声的分离，但该方法只针对处理核磁共振信号中的工频谐波噪声；CN104459809A公开了“一种基于独立成分分析的全波核磁共振信号噪声滤除方法”，采用独立成分分析算法对工频噪声进行消除，采用数字正交法构造虚拟输入通道信号解决欠定盲源分离问题，但该方法对强随机噪声和尖峰噪声干扰无能为力。CN106226407A公开了“一种基于奇异谱分析的超声回波信号在线预处理方法”，该方法将奇异谱分析用于超声在线检测技术中的回波信号预处理，可自动实现超声检测中回波信号的噪声去除和不同频段信号分量的分离提取；CN106404386A公开了“一种用于采集、提取及诊断齿轮箱早期故障特征信号的方法”，该方法将奇异谱分析用于故障诊断中。可见奇异谱分析已被成功应用到了信号处理的各个领域，但尚未见其应用于MRS信号的噪声滤除中。

发明内容

本发明所要解决的技术问题在于提供一种基于频率选择奇异谱分析的磁共振测深信号提取方法，解决由于环境中复杂尖峰噪声、工频谐波和随机噪声等影响带来的MRS信号提取问题。

本发明是这样实现的,基于频率选择奇异谱分析(Frequency chosen ofSingular-Spectrum-Analysis，F-SSA)的磁共振测深信号提取方法，该方法包括以下步骤：

步骤(1)：利用核磁共振测深(MRS)探水仪采集到一组已知Larmor频率的观测MRS信号X₁(t)；

步骤(2)：将采集的观测MRS信号X₁(t)通过宽频带带通滤波器得到X₂(t)；

步骤(3)：将通过宽频带带通滤波器得到的X₂(t)进行功率谱分析，将各频率的幅值进行降序排列，找到Larmor频率对应的信号幅值的排序位置；

步骤(4)：将通过宽频带带通滤波器得到的X₂(t)进行奇异谱分析提取MRS信号，奇异谱分析包括：嵌入，得到轨迹矩阵H；将轨迹矩阵H进行RSVD分解，得到按降序排列的奇异值和奇异向量U；根据步骤(3)得到的排序位置，选择与之对应的两个奇异值进行矩阵重构，得到矩阵C；将矩阵C进行对角平均化，得到提取的MRS信号Y(t)；

进一步地，步骤(4)所述的嵌入的具体步骤为：

信号X₂(t)为长度为N的一维实序列，X₂(t)＝(x₁,x₂,...,x_N)，正整数L为滑动窗口长度，1＜L＜N，L的取值由下式决定

通过嵌入操作原序列信号X₂(t)构成P个向量，每个向量可用h_i表示

h_i＝(x_i,x_i+1,…,x_i+L-1)^T

其中P＝N-L+1，i＝1,2,...,P，映射的结果形成轨迹矩阵H：

进一步地，步骤(4)所述的RSVD分解的具体步骤为：

1)、设置参数k和参数w，k为所取的用于近似重构矩阵的k个奇异值，w是用来保证重构矩阵条件成立的参数，其中w>k，w<L，w<P；

2)、构建0均值1方差的高斯随机矩阵G_L×w；

3)、计算轨迹矩阵H_L×P的采样矩阵M_P×w

4)、将采样矩阵M_P×w选取前k个奇异值进行SVD分解，得到其正交矩阵左奇异向量Q_P×k、对角矩阵Z_k×k，右奇异向量

5)、构建矩阵T_L×k：

T_L×k＝H_L×PQ_P×k

6)、对矩阵T_L×k进行SVD分解，得到其左奇异向量U_L×L、对角矩阵Σ_L×k和右奇异向量

7)、计算矩阵V_P×k：

V_P×k＝Q_P×kO_k×k；

8)、计算H_L×P的近似降秩矩阵

使用RSVD分解求得H的近似降秩矩阵，求取的近似降秩矩阵需满足以下条件：

λ_k+1为矩阵H的第k+1个奇异值；

对角阵Σ_L×k为L×k阶矩阵，其主对角元素为MRS信号近似降秩矩阵的k个奇异值；将k个奇异值进行降序排列(λ₁≥λ₂≥...≥λ_k)，得到与之对应的正交奇异向量集合U′＝(u₁,u₂,...u_k)，每个奇异值的贡献率为：

根据上式绘制出奇异谱图。

进一步地，步骤(4)所述的根据排序位置，选择与之对应的两个奇异值进行矩阵重构的具体步骤为：

根据步骤(3)找到MRS信号幅值大小的排序v,在奇异谱里面选择第2v-1和第2v个奇异值；

根据所选择的奇异值重构矩阵，重构步骤如下：

1)、计算重构信号矩阵的右奇异向量W:

其中，j＝2v-1,2v；

2)、重构信号矩阵C

进一步地，所述的步骤(4)中的对角平均化的具体步骤为：

通过对角平均化将奇异谱分析重构信号矩阵{c₁,c₂,…,c^q,…,c^P}转化为对应的重构序列{g¹,g²,…,g^q,…,g^P}，其中序列g^q表示第q个奇异谱分析重构序列，过程如下所示：

其中，c^q表示第q个奇异谱分析重构信号矩阵， 为矩阵c^q中第m行第n列元素，表示重构序列g^q的第d个元素，q＝1,2,...,P，d＝1,2,...,N；将重构序列{g₁,g₂,…,g^P}累加求和，得到去噪后的MRS信号序列其中Y(t)＝{y(t₁),y(t₂),…,y(t_N)}＝{y₁,y₂,…,y_N}。

本发明与现有技术相比，有益效果在于本发明提出了基于频率选择的奇异谱分析的磁共振测深信号的提取方法，针对单通道采集的探测数据，通过寻找MRS信号幅值在功率谱中的位置排序，采用所对应的奇异值进行矩阵重构和信号恢复，可以一次性去除尖峰噪声、工频谐波和随机噪声的干扰，实现了MRS信号的有效提取。本发明方法解决了磁共振测深找水工作中由于尖峰噪声、工频谐波和随机噪声造成的信号难以有效提取的难题，同时本发明突破了经典消噪方法需多通道探测等其他条件的限制，节省了大量的财力物力，开辟了奇异谱分析在核磁共振信号消噪领域的新天地。

附图说明

图1为本发明基于频率选择的奇异谱分析的磁共振测深信号的提取方法的流程框图；

图2为本发明RSVD分解算法流程框图；

图3为本发明含噪MRS信号与纯净MRS信号时域及其功率谱，其中(a)为时域谱，(b)为功率谱；

图4为本发明含噪MRS信号带通滤波器处理前后时域及其功率谱，其中(a)为时域谱，(b)为功率谱；

图5为本发明实例1选用带通滤波器曲线；

图6为本发明实例1奇异谱图；

图7为本发明仿真MRS信号F-SSA处理前后时域及其功率谱，其中(a)为时域谱，(b)为功率谱；

图8为本发明实测MRS信号带通滤波器处理前后时域及其功率谱，其中(a)为时域谱，(b)为功率谱；

图9为本发明实例2选用带通滤波器曲线；

图10为本发明实例2奇异谱图；

图11为本发明实测MRS信号F-SSA处理前后时域及其功率谱，其中(a)为时域谱，(b)为功率谱。

具体实施方式

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

如图1所示，基于频率选择奇异谱分析的磁共振测深信号提取方法，包括以下步骤：

步骤(1)：利用核磁共振测深(MRS)探水仪采集到一组已知Larmor频率的观测MRS信号X₁(t)；

步骤(2)：将采集的观测MRS信号X₁(t)通过宽频带带通滤波器得到X₂(t)；

步骤(3)：将通过宽频带带通滤波器得到的X₂(t)进行功率谱分析，将各频率的幅值进行降序排列，找到Larmor频率对应的信号幅值的排序位置；

步骤(4)：将通过宽频带带通滤波器得到的X₂(t)进行奇异谱分析提取MRS信号。奇异谱分析包括：嵌入，得到轨迹矩阵H；将轨迹矩阵H进行RSVD分解，得到按降序排列的奇异值和奇异向量U；根据步骤(3)得到的排序位置，选择与之对应的两个奇异值进行矩阵重构，得到矩阵C；将矩阵C进行对角平均化，得到提取的MRS信号Y(t)；

一种基于频率选择奇异谱分析的磁共振测深信号提取方法所述的嵌入的具体步骤为：

信号X₂(t)为长度为N的一维实序列，X₂(t)＝(x₁,x₂,...,x_N)，正整数L为滑动窗口长度，1＜L＜N，L的取值由下式决定

通过嵌入操作原序列信号X₂(t)构成P个向量，每个向量可用h_i表示

h_i＝(x_i,x_i+1,…,x_i+L-1)^T

其中P＝N-L+1，i＝1,2,...,P，映射的结果形成轨迹矩阵H：

如图2所示，本发明提供的RSVD分解的具体步骤为：

1)、设置参数k和参数w，k为所取的用于近似重构矩阵的k个奇异值，w是用来保证重构矩阵条件成立的参数，其中w>k，w<L，w<P；

2)、构建0均值1方差的高斯随机矩阵G_L×w；

3)、计算轨迹矩阵H_L×P的采样矩阵M_P×w

4)、将采样矩阵M_P×w选取前k个奇异值进行SVD分解，得到其正交矩阵左奇异向量Q_P×k、对角矩阵Z_k×k，右奇异向量

5)、构建矩阵T_L×k：

T_L×k＝H_L×PQ_P×k

6)、对矩阵T_L×k进行SVD分解，得到其左奇异向量U_L×L、对角矩阵Σ_L×k和右奇异向量

7)、计算矩阵V_P×k：

V_P×k＝Q_P×kO_k×k；

8)、计算H_L×P的近似降秩矩阵H～_L×P：

使用RSVD分解求得H的近似降秩矩阵，求取的近似降秩矩阵需满足以下条件：

λ_k+1为矩阵H的第(k+1)个奇异值；

对角阵Σ_L×k为L×k阶矩阵，其主对角元素为MRS信号近似降秩矩阵的k个奇异值；将k个奇异值进行降序排列(λ₁≥λ₂≥...≥λ_k)，得到与之对应的正交奇异向量集合U′＝(u₁,u₂,...u_k)，每个奇异值的贡献率为：

根据上式绘制出奇异谱图。

一种基于频率选择奇异谱分析的全波核磁共振信号提取方法所述的的矩阵重构的具体步骤为：

根据步骤(3)找到MRS信号幅值大小的排序v,在奇异谱里面选择第2v-1和第2v个奇异值；

根据所选择的奇异值重构矩阵，重构步骤如下：

1)、计算重构信号矩阵的右奇异向量W

其中，j＝2v-1,2v；

2)、重构信号矩阵C

一种基于频率选择奇异谱分析的全波核磁共振信号提取方法所述的的对角平均化的具体步骤为：

通过对角平均化将奇异谱分析重构信号矩阵{c₁,c₂,…,c^q,…,c^P}转化为对应的重构序列{g¹,g²,…,g^q,…,g^P}，其中序列g^q表示第q个奇异谱分析重构序列，过程如下所示：

其中，c^q表示第q个奇异谱分析重构信号矩阵， 为矩阵c^q中第m行第n列元素，表示重构序列g^q的第d个元素，q＝1,2,...,P，d＝1,2,...,N。将重构序列{g₁,g₂,…,g^P}累加求和，得到去噪后的MRS信号序列其中Y(t)＝{y(t₁),y(t₂),…,y(t_N)}＝{y₁,y₂,…,y_N}。

实施例1

本实施例是在MATLAB 2015a编程环境下开展的本发明方法的仿真实验。基于频率选择奇异谱分析的磁共振测深信号提取方法的仿真算法，参照图1，包括以下步骤：

步骤(1)：构造拉莫尔频率为2114Hz，幅值e₀为180nV，弛豫时间T₂ ^*为0.1s，相位为1.03的纯净MRS信号，信号采样率为10kHz，信号长度为500ms，数据点数为5000，如图3(a)和图3(b)所示。在该信号的基础上，仿真16组含噪信号，每组信号在0Hz～5000Hz之间添加100个幅值为100nV，相位随机，频率为50Hz整数倍的工频干扰；幅值为100nV，持续时间为10ms的尖峰噪声和幅值为60nV的随机噪声。将16组噪声进行统计叠加形成信噪比为-13.3930dB的观测MRS信号X₁(t)(为行向量)，时域图如图3(a)所示，功率谱图如图3(b)所示；

步骤(2)：将观测MRS信号X₁(t)经过如图5所示的带通滤波器，所采用的带通滤波器为切比雪夫滤波器，通带左边界为1900Hz，通带右边界为2300Hz，阻带左边界为1700Hz，阻带右边界为2500Hz，通带边带区衰减0.1dB，阻带截止区衰减30dB，得到信号X₂(t)，如图4(a)所示。

步骤(3)：求观测MRS信号X₂(t)的功率谱，如图4(b)所示。确定含噪信号中的Larmor信号幅值排在第一位；

步骤(4)：将通过宽频带带通滤波器得到的X₂(t)进行嵌入，选择的窗口长度为得到轨迹矩阵H；

将轨迹矩阵H进行RSVD分解，令k＝20，w＝200， 满足近似条件，得到按降序排列的奇异值和奇异向量U，绘制出奇异谱如图6；

根据步骤(3)得到的排序位置，Larmor信号幅值排在第一位，因此选择第一个和第二个两个奇异值进行矩阵重构，得到矩阵X；

将矩阵X进行对角平均化，得到提取的MRS信号Y(t)，时域图和功率谱如图7(a)和图7(b)所示；

为了验证本发明方法的实用性，将去噪后MRS信号Y(t)进行了信噪比(SNR)估计。经计算，其SNR＝9.2147dB，较分离前的SNR提高了22.6076dB；接着对Y(t)进行了包络提取和数据拟合，以获得分离信号的关键参数初始振幅e₀和弛豫时间T₂ ^*，计算可得，相对误差分别为1.6024％、3.0142％，均控制在±5％以内，满足应用要求。

实施例2

本实施例以长春市烧锅镇(该地拉莫尔频率约为2332.5Hz)采集的MRS信号作为本发明方法的处理对象。如图1所示，基于频率选择奇异谱分析的磁共振测深信号提取方法，包括以下步骤：

步骤(1)：利用核磁共振测深(MRS)探水仪采集到一组观测MRS信号X₁(t)(为行向量)，如图8(a)和图8(b)所示，采样率25kHz，数据点数为12476，计算其信噪比为SNR＝-7.1453dB；

步骤(2)：将观测MRS信号X₁(t)经过如图9所示的带通滤波器，所采用的带通滤波器为切比雪夫滤波器，通带左边界为2100Hz，通带右边界为2500Hz，阻带左边界为1800Hz，阻带右边界为2800Hz，通带边带区衰减0.1dB，阻带截止区衰减30dB，得到信号X₂(t)，如图8(a)所示。

步骤(3)：求观测MRS信号X₂(t)的功率谱，如图8(b)蓝色所示。确定含噪信号中的Larmor信号幅值排在第二位；

步骤(4)：将通过宽频带带通滤波器得到的X₂(t)进行嵌入，选择的窗口长度为得到轨迹矩阵H；

将轨迹矩阵H进行RSVD分解，令k＝20，w＝200， 满足近似条件，得到按降序排列的奇异值和奇异向量U，绘制出奇异谱如图10；

根据步骤(3)得到的排序位置，Larmor信号幅值排在第二位，因此选择第三个和第四个两个奇异值进行矩阵重构，得到矩阵X；

将矩阵X进行对角平均化，得到提取的MRS信号Y(t)，时域图和功率谱如图11(a)和图11(b)所示；

为了验证本发明方法的实用性，将去噪后MRS信号Y(t)进行了信噪比(SNR)估计。经计算，其SNR＝11.9542dB，较分离前的SNR提高了19.0995dB；接着对Y(t)进行了包络提取和数据拟合，以获得分离信号的关键参数初始振幅e₀和弛豫时间T₂ ^*，计算可得，e₀＝63.7791nV，T₂ ^*＝0.1890s，与实际水文地质钻孔资料结果相符合。

Claims (5)

1.一种基于频率选择奇异谱分析的磁共振测深信号提取方法，其特征在于，包括以下步骤：

步骤(1)：利用核磁共振测深探水仪采集到一组已知Larmor频率的观测MRS信号X₁(t)；

步骤(2)：将采集的观测MRS信号X₁(t)通过宽频带带通滤波器得到信号X₂(t)；

步骤(3)：将通过宽频带带通滤波器得到的信号X₂(t)进行功率谱分析，将各频率的幅值进行降序排列，找到Larmor频率对应的信号幅值的排序位置；

步骤(4)：将通过宽频带带通滤波器得到的信号X₂(t)进行奇异谱分析提取MRS信号，其中：奇异谱分析包括：

嵌入，得到轨迹矩阵H；

将轨迹矩阵H进行RSVD分解，得到按降序排列的奇异值和奇异向量U；

根据步骤(3)得到的排序位置，选择与之对应的两个奇异值进行矩阵重构，得到矩阵C；

将矩阵C进行对角平均化，得到提取的MRS信号Y(t)。

2.按照权利要求1所述的一种基于频率选择奇异谱分析的磁共振测深信号提取方法，其特征在于，所述的步骤(4)中的嵌入的具体步骤为：

信号X₂(t)为长度为N的一维实序列：X₂(t)＝(x₁,x₂,...,x_N)，正整数L为滑动窗口长度，1＜L＜N，L的取值由下式决定：

通过嵌入操作原序列信号X₂(t)构成P个向量，每个向量用h_i表示

h_i＝(x_i,x_i+1,…,x_i+L-1)^T

其中P＝N-L+1，i＝1,2,...,P，映射的结果形成轨迹矩阵H：

3.按照权利要求1所述的一种基于频率选择奇异谱分析的磁共振测深信号提取方法，其特征在于，所述的步骤(4)中的RSVD分解的具体步骤为：

1)、设置参数k和参数w，k为所取的用于近似重构矩阵的k个奇异值，w是用来保证重构矩阵条件成立的参数，其中w>k，w<L，w<P；

2)、构建0均值1方差的高斯随机矩阵G_L×w；

3)、计算轨迹矩阵H_L×P的采样矩阵M_P×w：

4)、将采样矩阵M_P×w选取前k个奇异值进行SVD分解，得到其正交矩阵左奇异向量Q_P×k、对角矩阵Z_k×k，右奇异向量表示为：

5)、构建矩阵T_L×k：

T_L×k＝H_L×PQ_P×k；

6)、对矩阵T_L×k进行SVD分解，得到其左奇异向量U_L×L、对角矩阵Σ_L×k和右奇异向量表示为：

7)、计算矩阵V_P×k：

V_P×k＝Q_P×kO_k×k；

8)、计算H_L×P的近似降秩矩阵

使用RSVD分解求得H的近似降秩矩阵，求取的近似降秩矩阵需满足以下条件：

λ_k+1为轨迹矩阵H的第k+1个奇异值；

对角矩阵Σ_L×k为L×k阶矩阵，其主对角元素为MRS信号近似降秩矩阵的k个奇异值；将k个奇异值进行降序排列(λ₁≥λ₂≥...≥λ_k)，得到与之对应的正交奇异向量集合U′＝(u₁,u₂,...u_k)，每个奇异值的贡献率为：

根据上式绘制出奇异谱图。

4.按照权利要求3所述的一种基于频率选择奇异谱分析的磁共振测深信号提取方法，其特征在于，所述的步骤(4)中的根据排序位置，选择与之对应的两个奇异值进行矩阵重构包括：

根据步骤(3)找到MRS信号幅值大小的排序v，在奇异谱里面选择第2v-1和第2v个奇异值；

根据所选择的奇异值重构矩阵，重构步骤如下：

1)、计算重构信号矩阵的右奇异向量W

其中，j＝2v-1,2v；

2)、重构信号矩阵C

5.按照权利要求1所述的一种基于频率选择奇异谱分析的磁共振测深信号提取方法，其特征在于，所述的步骤(4)中的对角平均化的具体步骤为：

通过对角平均化将奇异谱分析重构信号矩阵{c₁,c₂,…,c^q,…,c^P}转化为对应的重构序列{g¹,g²,…,g^q,…,g^P}，其中序列g^q表示第q个奇异谱分析重构序列，过程如下所示：

其中，c^q表示第q个奇异谱分析重构信号矩阵， 为矩阵c^q中第m行第n列元素，表示重构序列g^q的第d个元素，q＝1,2,...,P，d＝1,2,...,N；将重构序列{g₁,g₂,…,g^P}累加求和，得到去噪后的MRS信号序列其中Y(t)＝{y(t₁),y(t₂),…,y(t_N)}＝{y₁,y₂,…,y_N}。