CN116401513A - 一种基于深度残差网络的磁共振工频谐波噪声抑制方法 - Google Patents

Abstract

本发明涉及磁共振信号处理领域，为一种基于深度残差网络的磁共振工频谐波噪声抑制方法。首先将混叠了工频谐波噪声的含噪信号与工频谐波噪声作为深度残差网络的输入和输出，通过前向传播和反向传播计算误差，并引入残差块建立恒等映射，重复训练使损失函数减少至期望值，即可确定去噪模型。通过模型提取噪声并与含噪信号作差，从而获得去除工频噪声的信号。与传统方法相比，本发明利用训练过的模型进行预测的时间非常短，并且不需要调整网络参数。与卷积网络相比，残差网络在高深度时解决了退化问题，精度更高，模型效果更好，更加节省算力。

Description

一种基于深度残差网络的磁共振工频谐波噪声抑制方法

技术领域

本发明属于核磁共振测深(Magnetic Resonance Sounding，MRS)信号噪声抑制方法领域，特别是基于深度残差网络(Deep Residual Network，DRN)的磁共振地下水探测工频谐波噪声抑制方法。

背景技术

核磁共振技术探测地下水是目前唯一的直接找水的地球物理方法。与传统的地球物理探测方法相比具有高分辨力、高效率、信息量丰富和解唯一性等优点。MRS信号的质量对于探测结果是至关重要的。信号参数估计的准确性依赖于测量地点的电磁噪声水平。由于通过接收线圈接收到的核磁共振信号幅度极小，只有nV级，环境噪声足以将其淹没，影响信号质量最严重的两类噪声是工频谐波噪声和尖峰噪声，其中工频及其谐波噪声来源于电力线、发电机和变电器等，频率固定为50Hz或60Hz的整数倍。工频谐波噪声对接收信号的影响最为严重，且拉莫尔频率附近的工频干扰尤为突出。而工频谐波的存在，会使接收装置中的信号调理模块提前达到饱和状态，影响了后续的数据反演效果，从而降低了核磁共振探测结果的准确性。

Legchenko,A等人在Journalof Applied Geophyics[2003，53(2-3)，103-120]上发表的论文“Removal of power-line harmonics from proton magnetic resonancemeasurements”中基于NUMIS进行了工频谐波消噪的研究，利用区块对消等数据处理方法进行工频谐波干扰的消除。Larsen,J等人在Geophysical Journal International[2014，196(2)，828-836]上发表的论文“Noise cancelling of MRS signals combining model-based removal of powerline harmonics and multichannel Wiener filtering”中提出基于工频谐波建模和多通道维纳滤波相结合的方法来进行工频谐波噪声抑制。李玥于2022年12月在论文“地面磁共振工频谐波噪声消减及全波信号提取方法研究”中提出一种工频谐波噪声处理的新方法，针对工频谐波噪声干扰严重且来源密集化的问题，提出了频域对称算法。利用工频谐波噪声经过傅里叶变换和解调等过程产生的虚部频谱，估计与MRS信号混叠在一起的实部频谱从而消除噪声。专利CN111860273A公开了一种“基于卷积神经网络的磁共振地下水探测噪声抑制方法”，其特点是利用卷积神经网络的学习性，使用模拟训练数据集来优化模型参数，待模型训练完成后即可用其进行实测数据的消噪处理。

上述发明的传统方法需要人工调整滤波参数，不利于非专业技术人员使用。基于卷积神经网络的消噪方法普适性强，但针对工频噪声的滤除效果有限，在强工频噪声环境下无法达到预期效果。

发明内容

本发明所要解决的技术问题在于提供一种基于深度残差网络的磁共振工频谐波噪声抑制方法，该方法利用深度残差网络的学习性训练出一套成型的消除工频谐波噪声的模型，向模型输入含工频谐波噪声的MRS信号，模型自动提取工频谐波噪声并将其消除。与传统方法相比，本发明利用训练过的模型进行预测的时间非常短，并且不需要调整网络参数。

本发明是这样实现的，

一种基于深度残差网络的磁共振工频谐波噪声抑制方法，该方法包括：

形成数据集：建立仿真磁共振地下水探测仪采集的磁共振信号数据集X(n)，加入多种频率的工频及其谐波噪声N(n)，得到含工频谐波噪声的磁共振信号Y(n)，形成训练数据集，其中每条信号的长度为N，每个数据集的信号个数为M；

构建网络：设计残差网络的结构，确定网络深度及残差块种类，将网络参数初始化；

训练模型：输入训练数据集，前向传播时计算误差，反向传播时利用梯度下降法更新网络参数，利用残差学习机制建立含噪信号与噪声之间的映射关系，使损失函数逐渐优化至最小；

测试模型：组建基于实际测量数据的测试数据集，利用训练好的模型进行降噪处理。

进一步地，残差网络包括：

卷积运算，对输入数据进行特征提取，第一层的卷积运算调整输入含噪信号的尺寸用于池化；

池化层，用于提取数据的特征信息，过滤冗余信息；

残差块，用于在深层网络中防止神经网络退化，使网络随着深度的增加获得更高的精度；

全连接层，设置在放在残差网络的末尾，使输出与输入尺寸相同。

进一步地，构建残差网络的具体步骤为：步骤a：选定残差块种类，确定残差网络结构的层数L、残差块个数n和每层隐藏层的类型及参数，定义卷积核的大小K、填充大小P以及步幅S；

步骤b:初始化网络参数θ＝{W,b}、学习速率β及停止迭代阈值ε，其中W、b分别表示各残差块内部的隐藏层与输出层的线性关系系数矩阵和偏倚向量。

进一步地，训练模型的具体步骤包括：

将含工频谐波噪声的磁共振信号Y(n)作为残差网络的输入，工频谐波噪声N(n)作为残差网络的输出，构造去工频谐波噪声的残差网络模型：

x＝y-RN(y；θ)

其中，RN表示搭建的残差网络，θ＝{W,b}为网络参数，网络参数会在反向传播中不断更新，将仿真的含工频谐波噪声的核磁信号Y(n)输入进残差网络，含噪的磁共振信号数据y经网络前向传播算法输出工频谐波噪声数据RN(y；θ)，将y减去噪声数据RN(y；θ)得到干净的信号x；

其中进行前向传播算法，将含噪数据作为残差网络的输入，引入残差块，所述残差块由多个卷积层加上一个恒等映射构成，设F为残差函数，h为恒等映射，xⁿ为第n个残差块的输入，yⁿ为第n个残差块的输出，Wⁿ为第n个残差块的内部关系矩阵，则有：

xⁿ⁺¹＝yⁿ＝h(xⁿ)+F(xⁿ,Wⁿ)；

设单个卷积层输入张量为a，对应内部系数矩阵为W_l，偏倚向量为b_l，则卷积层第L层的输入张量、内部系数矩阵、偏倚向量分别为a^L、W_l ^L和

经过卷积运算及激活操作得到第L+1层网络的输入张量：

其中，ReLU(x)＝max(0,x)表示激活函数，*表示卷积操作，在单个卷积层内，卷积运算和激活操作之间利用BN算法进行批量归一化，来达到加速网络训练的目的。

进一步地，所述残差块的层数为两层，若a^L为第n个两层残差块的输入张量，a^L＝xⁿ，根据残差块结构则有

F(xⁿ,Wⁿ)＝xⁿ⁺¹＝a^L+2

以残差块为单位进行前向传播：

建立损失函数

J(θ)越小，预测的噪声输出与实际噪声越接近；

进行反向传播算法：在残差网络中引入残差δ，利用损失函数J(θ)计算第n个残差块的残差：

则得到

其中bⁿ为第n个残差块对应的偏倚向量，利用梯度下降法更新网络参数：

(W^k)’←W^k-β·x^kδ^k+1

(b^k)’←b^k-β·δ^k+1，k＝1,2,...,n-1,θ＝{W,b}，β为学习速率；

当各隐藏层与输出层的线性关系系数矩阵W和偏倚向量b小于停止更新的阈值ε时，停止梯度下降法并将W和b输出，其中W、b分别表示各残差块内部的隐藏层与输出层的线性关系系数矩阵和偏倚向量。

进一步地，构建训练数据集的具体步骤为：

步骤a：使用MATLAB构建仿真的磁共振信号，表达式为

其中t为时间，E₀是初始振幅，T₂ ^*是弛豫时间，f_L是拉莫尔频率，/>

是初始相位，得到X₁(n)，重复m次，得到信号个数为M，采集数据长度为N的核磁共振信号数据集X(n)＝[X₀(n),X₁(n),X₂(n),...,X_M-1(n)]；

步骤b：使用MATLAB构建仿真的工频谐波信号，表达式为

其中k为工频谐波的次数，f_N是基波频率，为50±0.1或60±0.1Hz，E_k和/>

是k次谐波的初始振幅和相位，改变谐波次数k以构建不同频率，不同振幅的工频谐波噪声N(n)＝[N₀(n),N₁(n),N₂(n),...,N_M-1(n)]，信号个数为M，采集数据长度为N；

步骤c：根据Y(n)＝X(n)+N(n)，通过改变初始振幅、初始相位、谐波次数来构造多组的核磁共振信号数据集X(n)、和工频谐波噪声信号N(n)，得到含噪的磁共振信号Y(n)。

本发明与现有技术相比，有益效果在于：本发明在模型训练完成后，无需调节滤波参数，可直接消除工频谐波噪声，且针对未训练的工频谐波噪声也有较好的去噪效果；与传统的去噪方法依赖人工经验、处理效率低相比，本发明消除工频谐波噪声的效率得到了很大的提高，可以很好的消除磁共振信号中的工频谐波噪声，达到理想的去噪效果；本发明能有效抑制工频谐波噪声造成的干扰，与传统的消噪方法相比具有更好的抗干扰能力，提高了可靠性，保障磁共振探测地下水的方法可以在有较强工频谐波噪声干扰的环境中使用。同时，与卷积神经网络相比，本发明引入了残差块，有效避免了神经网络的退化问题，针对工频噪声获得更高的精度。

附图说明

图1为基于深度残差网络的磁共振工频谐波噪声抑制方法流程图；

图2为深度残差网络结构图；

图3为深度残差网络中应用的双层残差块结构图。

具体实施方式

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

参见图1所示，一种基于深度残差网络的磁共振工频谐波噪声抑制方法，包括以下步骤：

a，形成数据集：建立仿真磁共振地下水探测仪采集的磁共振信号数据集X(n)，加入多种频率的工频及其谐波噪声N(n)，得到含工频谐波噪声的磁共振信号Y(n)，形成训练数据集,其中每条信号的长度为N，每个数据集的信号个数为M；

b，构建网络：如图2所示的深度残差网络结构，确定网络深度及残差块种类，将网络参数初始化；

c，训练模型：输入训练数据集，前向传播时计算误差，反向传播传播时利用梯度下降法更新网络参数，利用残差学习机制建立含噪信号与噪声之间的映射关系，使损失函数逐渐优化至最小；

d，测试模型：组建基于实际测量数据的测试数据集，利用训练好的模型进行降噪处理。

所述步骤a中形成数据集的基本步骤为：

步骤a1：构建仿真的磁共振信号

其中t＝(0:N-1)/fsfs为采样率，N＝2000，E₀∈(10,4000)nV，/>

f_L∈(1300,3700)Hz，

得到X₀(n)，重复1000次，得到信号个数为M＝1000，采集数据长度为N＝2000的核磁共振信号数据集X(n)＝[X₀(n),X₁(n),X₂(n),...,X₉₉₉(n)]；

步骤a2：构建仿真的工频谐波信号，其表达式为

其中t₀＝(0:N-1)/fs，fs为采样率，N＝2000，该信号为工频噪声的45次谐波与47次谐波的叠加，其频率恰好在拉莫尔频率的附近。45次谐波的初始振幅为50nV，初始相位为/>

47次谐波的初始振幅为40nV，初始相位为

重复1000次，得到信号个数为M＝1000，采集数据长度为N＝2000的工频谐波噪声信号数据集N(n)＝[N₀(n),N₁(n),N₂(n),...,N₉₉₉(n)]。

步骤a3：根据Y(n)＝X(n)+N(n)，通过改变初始振幅、初始相位、谐波次数等参数来构造多组的核磁共振信号数据集X(n)、和工频谐波噪声信号N(n)，得到含噪的磁共振信号Y(n)。

所述步骤b中构建残差网络的具体步骤为：

步骤b1：选用图3所示的双层残差块，确定残差网络结构的层数L和每层隐藏层的类型及参数，定义卷积核的大小K，填充大小P，步幅S；

步骤b2:初始化网络参数θ＝{W,b}、学习速率β＝0.05及停止迭代阈值ε＝50nV。

所述步骤c中训练模型的具体步骤为：

步骤c1：含工频谐波噪声的磁共振信号Y(n)作为残差网络的输入，工频谐波噪声N(n)作为残差网络的输出，构造去工频谐波噪声的残差网络模型

x＝y-RN(y；θ)

其中，RN表示搭建的残差网络，θ＝{W,b}为网络参数。

进行前向传播算法，将含噪数据作为残差网络的输入。设F为残差函数，h为恒等映射，xⁿ为第n个残差块的输入，yⁿ为第n个残差块的输出，Wⁿ为第n个残差块的内部关系矩阵，则有

xⁿ⁺¹＝yⁿ＝h(xⁿ)+F(xⁿ,Wⁿ)；

以残差块为单位进行前向传播：

步骤c2：建立损失函数

J(θ)越小，预测的噪声输出与实际噪声越接近；

步骤c3：进行反向传播算法：利用梯度下降法更新网络参数：

(W^k)’←W^k-β·x^kδ^k+1

(b^k)’←b^k-β·δ^k+1，k＝1,2,...,n-1,θ＝{W,b},β为学习速率；

当各隐藏层与输出层的线性关系系数矩阵W和偏倚向量b小于停止更新的阈值ε时，停止梯度下降法并将W和b输出。

所述步骤d具体实施过程为：

在某MRS探测点进行实际数据采集获得数据集Y_test(n)，此地点的拉莫尔频率应与网络训练的拉莫尔频率相近。将Y_test(n)输入进训练好的残差网络RN(y；θ)中，得到预测的工频谐波噪声输出为N_test(n)，利用Y_test(n)-N_test(n)＝X_test(n)，得到消去工频谐波噪声的信号X_test(n)。

Claims (6)

1.一种基于深度残差网络的磁共振工频谐波噪声抑制方法，其特征在于，该方法包括：

形成数据集：建立仿真磁共振地下水探测仪采集的磁共振信号数据集X(n)，加入多种频率的工频及其谐波噪声N(n)，得到含工频谐波噪声的磁共振信号Y(n)，形成训练数据集，其中每条信号的长度为N，每个数据集的信号个数为M；

构建网络：设计残差网络的结构，确定网络深度及残差块种类，将网络参数初始化；

训练模型：输入训练数据集，前向传播时计算误差，反向传播时利用梯度下降法更新网络参数，利用残差学习机制建立含噪信号与噪声之间的映射关系，使损失函数逐渐优化至最小；

测试模型：组建基于实际测量数据的测试数据集，利用训练好的模型进行降噪处理。

2.按照权利要求1所述的方法，其特征在于，残差网络包括：

卷积运算，对输入数据进行特征提取，第一层的卷积运算调整输入含噪信号的尺寸用于池化；

池化层，用于提取数据的特征信息，过滤冗余信息；

残差块，用于在深层网络中防止神经网络退化，使网络随着深度的增加获得更高的精度；

全连接层，设置在放在残差网络的末尾，使输出与输入尺寸相同。

3.按照权利要求2所述的方法，其特征在于，构建残差网络的具体步骤为：步骤a：选定残差块种类，确定残差网络结构的层数L、残差块个数n和每层隐藏层的类型及参数，定义卷积核的大小K、填充大小P以及步幅S；

步骤b:初始化网络参数θ＝{W,b}、学习速率β及停止迭代阈值ε，其中W、b分别表示各残差块内部的隐藏层与输出层的线性关系系数矩阵和偏倚向量。

4.按照权利要求1所述的方法，其特征在于，训练模型的具体步骤包括：

将含工频谐波噪声的磁共振信号Y(n)作为残差网络的输入，工频谐波噪声N(n)作为残差网络的输出，构造去工频谐波噪声的残差网络模型：

x＝y-RN(y；θ)

其中，RN表示搭建的残差网络，θ＝{W,b}为网络参数，网络参数会在反向传播中不断更新，将仿真的含工频谐波噪声的核磁信号Y(n)输入进残差网络，含噪的磁共振信号数据y经网络前向传播算法输出工频谐波噪声数据RN(y；θ)，将y减去噪声数据RN(y；θ)得到干净的信号x；

其中进行前向传播算法，将含噪数据作为残差网络的输入，引入残差块，所述残差块由多个卷积层加上一个恒等映射构成，设F为残差函数，h为恒等映射，xⁿ为第n个残差块的输入，yⁿ为第n个残差块的输出，Wⁿ为第n个残差块的内部关系矩阵，则有：

xⁿ⁺¹＝yⁿ＝h(xⁿ)+F(xⁿ,Wⁿ)；

设单个卷积层输入张量为a，对应内部系数矩阵为W_l，偏倚向量为b_l，则卷积层第L层的输入张量、内部系数矩阵、偏倚向量分别为a^L、W_l ^L和

经过卷积运算及激活操作得到第L+1层网络的输入张量：

其中，ReLU(x)＝max(0,x)表示激活函数，*表示卷积操作，在单个卷积层内，卷积运算和激活操作之间利用BN算法进行批量归一化，来达到加速网络训练的目的。

5.按照权利要求4所述的方法，其特征在于，所述残差块的层数为两层，若a^L为第n个两层残差块的输入张量，a^L＝xⁿ，根据残差块结构则有

F(xⁿ,Wⁿ)＝xⁿ⁺¹＝a^L+2

以残差块为单位进行前向传播：

建立损失函数

J(θ)越小，预测的噪声输出与实际噪声越接近；

进行反向传播算法：在残差网络中引入残差δ，利用损失函数J(θ)计算第n个残差块的残差：

则得到

其中bⁿ为第n个残差块对应的偏倚向量，利用梯度下降法更新网络参数：

(W^k)’←W^k-β·x^kδ^k+1

(b^k)’←b^k-β·δ^k+1，k＝1,2,...,n-1,θ＝{W,b}，β为学习速率；

当各隐藏层与输出层的线性关系系数矩阵W和偏倚向量b小于停止更新的阈值ε时，停止梯度下降法并将W和b输出，其中W、b分别表示各残差块内部的隐藏层与输出层的线性关系系数矩阵和偏倚向量。

6.按照权利要求1所述的方法，其特征在于，构建训练数据集的具体步骤为：

步骤a：使用MATLAB构建仿真的磁共振信号，表达式为

其中t为时间，E₀是初始振幅，T₂ ^*是弛豫时间，f_L是拉莫尔频率，/>

是初始相位，得到X₁(n)，重复m次，得到信号个数为M，采集数据长度为N的核磁共振信号数据集X(n)＝[X₀(n),X₁(n),X₂(n),...,X_M-1(n)]；

步骤b：使用MATLAB构建仿真的工频谐波信号，表达式为

其中k为工频谐波的次数，f_N是基波频率，为50±0.1或60±0.1Hz，E_k和/>

是k次谐波的初始振幅和相位，改变谐波次数k以构建不同频率，不同振幅的工频谐波噪声N(n)＝[N₀(n),N₁(n),N₂(n),...,N_M-1(n)]，信号个数为M，采集数据长度为N；

步骤c：根据Y(n)＝X(n)+N(n)，通过改变初始振幅、初始相位、谐波次数来构造多组的核磁共振信号数据集X(n)、和工频谐波噪声信号N(n)，得到含噪的磁共振信号Y(n)。

Images