CN104614778B - 基于ica的核磁共振地下水探测信号噪声消除方法 - Google Patents

Abstract

本发明公开了一种核磁共振地下水探测信号噪声消除方法，包括以下步骤：录入三组核磁共振响应数据，分别对这三组数据进行傅里叶变换，确定每组数据核磁共振中心频率附近所含工频谐波，然后构造与工频谐波同频率，与核磁共振相应数据同长度的正弦函数、余弦函数，并与核磁共振响应数据组成观测信号，采用独立分量分析算法对每组观测信号进行分离得到解混信号，进行数据重构以消除工频谐波的干扰，将三组去除工频谐波的核磁共振数据作为观测信号，再利用ICA算法处理，削弱剩余随机噪声干扰。本发明在消除工频谐波噪声干扰的同时，不会破坏信号的任何细节，不需要铺设参考线圈，操作简单，压制随机噪声不需要大量的数据，减少处理时间。

Description

基于ICA的核磁共振地下水探测信号噪声消除方法

技术领域：

本发明涉及一种核磁共振数据预处理方法，具体涉及一种基于独立分量分析(ICA)的核磁共振地下水探测信号噪声消除方法。

背景技术：

寻找和评价地下淡水资源的方法有很多种，常规的物探找水方法都是通过勘查含水构造的层位间接找水，不能解决何处有水、有多少水等一些与地下水紧密相关的基本问题，核磁共振地下水探测技术作为唯一的直接探测地下水的地球物理方法，具有信息量丰富、解唯一性等优点。但是核磁共振信号极其微弱，而且核磁共振探测仪器在接收核磁共振响应信号的同时不可避免地接收到复杂的噪声干扰，诸如雷电等自然干扰源引起的奇异噪声，由电力线、发电机和变电器等人为干扰源引起的工频谐波噪声，以及由其他噪声源引起的平稳噪声和时变噪声。复杂的噪声干扰使得核磁共振地下水探测信号基本被淹没，严重阻碍磁共振信号的提取，影响后期反演解释结果的准确性，因此核磁共振地下水探测信号消噪方法的研究意义重大。

目前已经应用到核磁共振地下水探测信号噪声消除的方法有很多种，其中利用非线性能量算子能够对尖峰噪声进行较好地消除，基于大量数据的统计叠加方法能够对随机噪声进行一定的压制，但是需要的数据量大，叠加次数多，工作时间长，效率较低，工频谐波噪声因为与核磁共振地下水探测信号的中心频率非常接近，难以消除。现有的工频谐波噪声消除的方法有陷波器、基于参考线圈和变步长自适应算法、多通道核磁响应信号的自适应抵消，但是陷波器在消除工频谐波噪声的同时也会破坏核磁共振信号中一些重要的细节信息，基于参考线圈和变步长自适应算法，以及多通道核磁响应信号的自适应抵消方法不仅需要铺设参考线圈，操作复杂，而且只能消除与参考通道相关的噪声，对参考通道的依赖性较大。

CN203759264U公开“一种基于工频整周期触发的核磁共振双极性叠加消噪装置”，是由双极性叠加消噪装置是由输出可调的大功率电源与大功率发射桥路连接，双极性主控制单元经发射桥路驱动和大功率发射桥路与谐振电容连接，双极性主控制单元分别与工频整周期同步触发单元、高压切换开关、信号调理电路和A/D采集单元连接，高压切换开关经信号调理电路和放大器电路与A/D采集单元连接构成。是通过分别采集两组方向相反的核磁共振信号和方向一致的噪声信号，通过相减的方法消除噪声，但是此方法不能保证两次采集的噪声完全一样，会有噪声残余，影响探测精度。

发明内容：

本发明的目的就在于针对上述现有技术的不足，结合ICA和核磁共振信号与噪声特性，提供一种基于ICA的核磁共振地下水探测信号噪声消除方法。

本发明的目的是通过以下技术方案实现的：

一种基于ICA的核磁共振地下水探测信号噪声消除方法，包括以下步骤：

A、录入三组核磁共振响应数据；

B、利用ICA独立成分分析，依次去除每组数据核磁共振信号中心频率附近的工频谐波干扰，然后分别对每组处理后的数据进ICA逆变换后再进行数据重构，得到消除工频谐波干扰的数据；

C、将去除工频谐波噪声的三组数据作为观测信号，然后对其进行ICA处理，以消除剩余的随机噪声，并利用ICA逆变换进行数据重构，得到最终的消噪数据。

步骤A所述的三组，是指三个源信号。因为，ICA算法的使用前提是观测信号的个数大于等于源信号的个数，为了计算简单通常取相同个数，三组数据相对合理。核磁共振数据包含MRS信号和噪声，其中噪声又包括工频谐波干扰、随机噪声和尖峰干扰等。本发明仅针对工频谐波干扰和随机噪声，因此，源信号分为MRS信号、工频谐波干扰和随机噪声。

步骤B包括以下步骤：

a、对第一组核磁共振响应数据进行快速傅里叶变换，确定该核磁共振地下水探测信号中心频率附近的工频谐波干扰频率f₁,f₂,…,f_n；

b、结合所确定频率，分别构造与工频谐波同频率，与核磁共振响应数据同长度的正弦函数、余弦函数，例如工频谐波频率为f₁，则正弦函数为：sin(2πf₁k)，余弦函数为：cos(2πf₁k)，并与核磁共振响应数据构成观测信号；

c、对观测信号进行预处理；

d、利用快速独立分量分析算法对预处理后的观测信号进行分离，得到幅度有很大衰减的解混信号和解混矩阵w；

e、将步骤d分离出的工频谐波噪声屏蔽清零，利用ICA逆变换 进行数据重构，恢复核磁共振信号的幅度，其中独立成分y_j是分解出来的信号，只保留独立成分y_j，其他独立成分置零，是每个通道内所含的信号成分，由于核磁共振响应数据通道所含的核磁共振信号明显，正弦函数及余弦函数所在的通道基本不含核磁共振信号，因此选择核磁共振信号强的作为去工频谐波后的核磁共振信号，记为

f、针对第二组核磁共振响应数据，重复步骤a～e得到去除工频谐波的第二组核磁共振响应数据，记为

g、针对第三组核磁共振响应数据，重复步骤a～e得到去除工频谐波的第三组核磁共振响应数据，记为

步骤C包括以下步骤：

Ⅰ、将三组去除工频谐波噪声的数据组成观测信号，并对其进行预处理；

Ⅱ、利用快速独立分量分析算法对预处理后的观测信号进行分离，得到幅度有很大衰减的解混信号和解混矩阵A；

Ⅲ、将解混信号中噪声成分少、核磁共振信号明显的成分保留，其他成分屏蔽，利用ICA逆变换进行数据重构，恢复核磁共振信号幅度。

ICA算法包括以下步骤：

第一步、观测信号去均值，首先按照公式求解观测信号每一行的均值，然后利用公式对观测信号进行去均值处理，使数据中心化，满足零均值的假设；

第二步、对去均值后的数据进行白化处理，根据公式C_x＝E{xx^T}求观测信号的协方差矩阵，并求协方差矩阵的特征值d＝diag(d₀,d₁,…,d_2n)，特征向量e＝(e₀,e₁,…e_2n)，最后由公式x＝vx＝d^1/2e^Tx对零均值观测信号进行白化处理，使数据具有单位方差，以消除数据各分量之间的相关性；

第三步、求解混矩阵w，首先令j＝1，初始化解混向量w₁，利用根据负熵最大化独立判据和牛顿迭代优化算法推导出递推式：

求w₁，

利用公式对w₁进行正交化，在根据w_j＝w_j/||w_j||标准化w₁，当w₁收敛时，第一个独立成分对应的解混向量w₁求解完毕，判断w₁是否收敛，如果w₁不收敛重新利用公式求w₁，对w₁进行正交化和标准化，直到w₁收敛，第一个独立成分对应的解混向量w₁求解完毕，j＝j+1，判断j≤m是否成立，如果j≤m，按照上述的步骤求出w₂，直到j＞m时，独立成分对应的所有的解混向量w求解完毕，输出解混矩阵w；

第四步、利用公式y＝w^Tx求解独立成分y；

第五步、数据重构，由于利用ICA算法求出的输出信号y具有幅度的不确定性，与源信号相比初始幅度明显减小，将第四步分离出的噪声屏蔽清零，利用ICA逆变换进行数据重构，恢复核磁共振信号的幅度，其中独立成分y_j是分解出来的信号，只保留独立成分y_j，其他独立成分置零，是每个通道内所含的信号成分。

有益效果：本发明对源信号和传输通道的先验知识没有要求，在试验过程中不需要铺设参考线圈，使工频噪声的消除效果不受参考线圈噪声采集情况的影响，在随机噪声消除的过程中不需要大量的数据进行叠加，操作简单、方便、效率高，由于采用了ICA算法，在分离噪声的过程中对信号的细节几乎没有破坏。

附图说明：

图1为基于ICA的核磁共振地下水探测信号工频谐波噪声消除方法流程图

图2为用ICA算法对去除工频噪声后随机噪声消除流程图

图3为ICA算法流程图

图4为对第一组数据利用ICA对工频谐波噪声进行消除时构造的观测信号；

图5为第一组数据ICA处理后的解混信号；

图6为第一组数据去除工频谐波后的核磁共振数据时域波形图；

图7为第一组数据去除工频谐波后频谱与原始数据频谱的对比图；

图8为第二组数据去除工频谐波后的核磁共振数据时域波形图；

图9为第二组数据去除工频谐波后频谱与原始数据频谱的对比图；

图10为第三组数据去除工频谐波后的核磁共振数据时域波形图；

图11为第三组数据去除工频谐波后频谱与原始数据频谱的对比图；

图12为三组去除工频谐波后的数据组成的观测信号时域波形图；

图13为第二次ICA处理后的解混信号时域波形图；

图14为第二次ICA处理后的重构信号时域波形图；

图15为未经处理的原始三组数据时域波形图；

图16为两次ICA消噪后核磁共振信号频谱与一个未经处理的原始信号频谱的对比图。

具体实施方式：

图1是基于ICA的核磁共振地下水探测信号工频谐波噪声消除方法的工作流程图。针对核磁共振地下水探测时受到的工频谐波的干扰，通过构造与工频谐波同频率的正弦函数、余弦函数，并和核磁共振响应信号组成观测信号，利用ICA算法对观测信号进行分离，并采用数据重构的方法解决ICA分离过程中产生的幅度衰减问题。

图2是利用ICA算法对去除工频噪声后的核磁共振地下水探测信号进行随机噪声消除的工作框图，将去工频后的三组核磁共振数据作为观测信号，用ICA算法对其进行处理，并采用数据重构的方法解决ICA分离过程中产生的幅度衰减问题，以压制核磁共振数据中的剩余噪声。

下面是基于ICA的核磁共振地下水探测信号噪声消除方法的具体步骤：

A、录入三组核磁共振响应数据，记为x₀(k)、x₁(k)、x₂(k)；

B、利用ICA算法依次去除x₀(k)、x₁(k)、x₂(k)中核磁共振信号中心频率附近的工频谐波干扰，具体步骤如下：

a、对数据x₀(k)进行快速傅里叶变换，确定该核磁共振地下水探测信号中心频率附近的工频谐波干扰频率f₁,f₂,…,f_n；

b、分别构造与工频谐波同频率，与核磁共振响应数据同长度的正弦函数(x_1sin＝sin(2πf₁k),x_2sin＝sin(2πf₂k),…,x_nsin＝sin(2πf_nk))、余弦函数(x_1cos＝cos(2πf₁k),x_2cos＝cos(2πf₂k),…,x_ncos＝cos(2πf_nk))，并与核磁共振响应数据构成观测信号x＝[x₀(k),x_1sin(k),x_1cos(k),x_2sin(k),x_2cos(k),…,x_nsin(k),x_ncos(k)]^T；

c、对观测信号进行预处理。首先按照公式依次求解观测信号每一行的均值，记为然后利用公式对观测信号进行去均值处理，在根据公式C_x＝E{xx^T}求观测信号的协方差矩阵，并求协方差矩阵的特征值d＝diag(d₀,d₁,…,d_2n)，特征向量e＝(e₀,e₁,…e_2n)，最后由公式x＝vx＝d^1/2e^Tx对零均值观测信号进行白化处理，得到预处理后的观测信号；

d、利用快速独立分量分析算法对预处理后的观测信号进行分离，得到幅度有很大衰减的核磁共振信号和工频谐波噪声和解混矩阵w，具体步骤如下：

i.确定观测信号的维数m＝2n+1，即为独立分量的个数。令j＝1；

ii.初始化解混向量w_j；

iii.按照公式w_j＝w_j/||w_j||归一化解混向量；

iv.选择非线性函数

v.由公式求解解混向量；

vi.根据公式正交化解混向量w_j，并归一化w_j；

vii.判断，如果w_j收敛，j＝j+1，否则返回步骤v；

viii.判断，如果j<m，返回步骤ii，如果j＝m，输出解混矩阵w；

ix.由公式y＝w^Tx估计独立分量。

e、将步骤d分离出的工频谐波噪声屏蔽清零，利用ICA逆变换 进行数据重构，恢复核磁共振信号的幅度，其中独立成分y_j是分解出来的信号，只保留独立成分y_j，其他独立成分置零，是每个通道内所含的信号成分，由于核磁共振响应数据通道所含的核磁共振信号明显，正弦函数及余弦函数所在的通道基本不含核磁共振信号，因此选择核磁共振信号强的作为去工频谐波后的核磁共振信号，记为

f、针对第二组核磁共振响应数据，重复步骤a～e得到去除工频谐波的第二组核磁共振响应数据，记为

g、针对第三组核磁共振响应数据，重复步骤a～e得到去除工频谐波的第三组核磁共振响应数据，记为

C、将去除工频谐波噪声的三组数据作为观测信号，然后其进行ICA处理，以消除剩余的随机噪声，并利用ICA逆变换进行数据重构，得到最终的消噪数据，具体步骤如下：

a、将三组去除工频谐波噪声的数据组成观测信号并对其进行预处理，首先按照公式依次求解观测信号每一行的均值，然后利用公式对观测信号进行去均值处理，在根据公式求观测信号的协方差矩阵，并求协方差矩阵的特征值d＝diag(d₀,d₁,d₃)，特征向量e＝(e₀,e₁,e₃)，最后由公式对零均值观测信号进行白化处理，得到预处理后的观测信号；

b、利用快速独立分量分析算法对预处理后的观测信号进行分离，得到幅度有很大衰减的解混信号和解混矩阵A，具体步骤如下：

i.确定观测信号的维数m＝3，即为独立分量的个数。令i＝1；

ii.初始化解混向量A_i；

iii.按照公式A_i＝A_i/||A_i||归一化解混向量；

iv.选择非线性函数

v.由公式求解解混向量；

vi.根据公式正交化解混向量A_i，并归一化A_i；

vii.判断，如果A_i收敛，i＝i+1，否则返回步骤v；

viii.判断，如果i<m，返回步骤ii，如果j＝m，输出解混矩阵A；

ix.按照公式y＝A^Tx估计独立分量。

c、将解混信号中噪声成分少、核磁共振信号明显的成分保留，其他成分屏蔽，利用ICA逆变换进行数据重构，恢复核磁共振信号幅度。

应用示例：

以吉林省烧锅镇核磁共振地下水探测为例：根据当地的地磁场强度计算得到核磁共振的拉莫尔频率为2326Hz，也就是核磁共振响应数据的中心频率为2326Hz，数据长度为16000，采样频率为66666.7Hz，在matlab环境下，基于ICA算法对核磁共振地下水探测实测数据中的工频谐波噪声进行处理。

具体步骤如下：

A、录入三组核磁共振响应数据，记为x₀(k)、x₁(k)、x₂(k)；

B、利用ICA算法依次去除x₀(k)、x₁(k)、x₂(k)中核磁共振信号中心频率附近的工频谐波干扰，具体步骤如下：

a、对数据x₀(k)进行快速傅里叶变换，确定该核磁共振地下水探测信号中心频率附近的工频谐波干扰频率f1＝2300Hz，f2＝2350Hz，f3＝2450Hz；

b、分别构造与工频谐波同频率，与核磁共振响应数据同长度的正弦函数(x_1sin＝sin(2πf₁k),x_2sin＝sin(2πf₂k),x_3sin＝sin(2πf₃k))、余弦函数(x_1cos＝cos(2πf₁k),x_2cos＝cos(2πf₂k),x_3cos＝cos(2πf₃k))，并与核磁共振响应数据构成观测信号x＝[x₀(k),x_1sin(k),x_1cos(k),x_2sin(k),x_2cos(k),x_3sin(k),x_3cos(k)]^T，观测信号如图4所示；

c、对观测信号进行预处理。首先按照公式依次求解观测信号每一行的均值，记为然后利用公式对观测信号进行去均值处理，在根据公式C_x＝E{xx^T}求观测信号的协方差矩阵，并求协方差矩阵的特征值d＝diag(d₀,d₁,d₂,d₃,d₄,d₅,d₆)，特征向量e＝(e₀,e₁,e₂,e₃,e₄,e₅,e₆)，最后由公式x＝vx＝d^1/2e^Tx对零均值观测信号进行白化处理，得到预处理后的观测信号；

d、利用快速独立分量分析算法对预处理后的观测信号进行分离，得到幅度有很大衰减的核磁共振信号和工频谐波噪声和解混矩阵w，分离后的独立分量估计如图5所示，具体步骤如下：

i.确定观测信号的维数m＝7，即为独立分量的个数。令j＝1；

ii.初始化解混向量w_j；

iii.按照公式w_j＝w_j/||w_j||归一化解混向量；

iv.选择非线性函数

v.由公式求解解混向量；

vi.根据公式正交化解混向量w_j，并归一化w_j；

vii.判断，如果w_j收敛，j＝j+1，否则返回步骤v；

viii.判断，如果j<7，返回步骤ii，如果j＝7，输出解混矩阵w；

ix.按照公式y＝w^Tx估计独立分量。

e、将步骤d分离出的工频谐波噪声屏蔽清零，利用ICA逆变换 进行数据重构，恢复核磁共振信号的幅度，其中独立成分y₀是分解出来的信号，是每个通道内所含的信号成分，由于核磁共振响应数据通道所含的核磁共振信号明显，正弦函数及余弦函数所在的通道基本不含核磁共振信号，因此选择核磁共振信号强的作为去工频谐波后的核磁共振信号，记为重构后核磁共振信号如图6所示，消噪前后的频谱分析对比图如图7所示。

f、针对第二组核磁共振响应数据，重复步骤a～e得到去除工频谐波的第二组核磁共振响应数据，记为重构核磁共振信号如图8所示，消噪前后的频谱分析对比图如图9所示。

g、针对第三组核磁共振响应数据，重复步骤a～e得到去除工频谐波的第三组核磁共振响应数据，记为重构核磁共振信号如图10所示，消噪前后的频谱分析对比图如图11所示。

C、将去除工频谐波噪声的三组数据作为观测信号，并对其进行ICA处理，以消除剩余的随机噪声，得到最终的消噪数据，具体步骤如下：

a、将三组去除工频谐波噪声的数据组成观测信号观测信号如图12所示，并对其进行预处理，首先按照公式依次求解观测信号每一行的均值，然后利用公式对观测信号进行去均值处理，在根据公式求观测信号的协方差矩阵，并求协方差矩阵的特征值d＝diag(d₀,d₁,d₃)，特征向量e＝(e₀,e₁,e₃)，最后由公式对零均值观测信号进行白化处理，得到预处理后的观测信号；

b、利用快速独立分量分析算法对预处理后的观测信号进行分离，得到幅度有很大衰减的解混信号和解混矩阵A，分离后的独立分量估计如图13所示，具体步骤如下：

i.确定观测信号的维数m＝3，即为独立分量的个数，令i＝1；

ii.初始化解混向量A_i；

iii.按照公式A_i＝A_i/||A_i||归一化解混向量；

iv.选择非线性函数

v.根据公式求解解混向量；

vi.由公式正交化解混向量A_i，并归一化A_i；

vii.判断，如果A_i收敛，i＝i+1，否则返回步骤v；

viii.判断，如果i<3，返回步骤ii，j＝3，输出解混矩阵A；

ix.按照公式y＝A^Tx估计独立分量。

c、将解混信号中噪声成分少、核磁共振信号明显的成分保留，其他成分屏蔽，利用ICA逆变换进行数据重构，恢复核磁共振信号幅度，重构后的核磁共振信号如图14所示。

基于ICA算法的核磁共振地下水探测信号工频谐波噪声消除方法，不仅使核磁共振响应数据中的工频谐波被很好地消除，而且没有破坏核磁共振信号，对其他的噪声也有很好的保留，这样不会干扰后期对其他噪声的处理，基于ICA算法的核磁共振地下水探测信号随机噪声消除方法，很好地消除了随机噪声，消噪后的核磁共振信号的时域波形图很接近理想的呈e指数衰减的核磁共振信号，ICA算法对实测数据的结果分别如图14、图16。其中图14是两次消噪后核磁共振信号的时域波形图，而图15是消噪前核磁共振实测数据的时域波形图，对比图14、图15的结果，可以发现ICA消噪后的核磁共振信号e指数衰减的趋势更明显，图16是消噪后核磁共振信号与消噪前的其中一个实测数据频谱分析对比图，可以发现ICA消噪后核磁共振信号的频谱噪声成分更少。

Claims (3)

1.一种基于ICA的核磁共振地下水探测信号噪声消除方法，其特征在于：包括以下步骤：

A、录入三组核磁共振响应数据；

B、利用ICA独立成分分析，依次去除每组数据核磁共振信号中心频率附近的工频谐波干扰，然后分别对每组处理后的数据进ICA逆变换后再进行数据重构，得到消除工频谐波干扰的数据；

C、将去除工频谐波噪声的三组数据作为观测信号，然后对其进行ICA处理，以消除剩余的随机噪声，并利用ICA逆变换进行数据重构，得到最终的消噪数据。

2.按照权利要求1所述的基于ICA的核磁共振地下水探测信号噪声消除方法，其特征在于，步骤B包括以下步骤：

a、对第一组核磁共振响应数据进行快速傅里叶变换，确定该核磁共振地下水探测信号中心频率附近的工频谐波干扰频率f₁，f₂，…，f_n；

b、结合所确定频率，分别构造与工频谐波同频率，与核磁共振响应数据同长度的正弦函数、余弦函数，工频谐波频率为f_n，n＝1，2，…，n，则正弦函数为：sin(2πf_nk)，余弦函数为：cos(2πf_nk)，并与核磁共振响应数据构成观测信号；

c、对观测信号进行预处理；

d、利用快速独立分量分析算法对预处理后的观测信号进行分离，得到幅度有很大衰减的解混信号和解混矩阵w；

e、将步骤d分离出的工频谐波噪声屏蔽清零，利用ICA逆变换

进行数据重构，恢复核磁共振信号的幅度，其中独立成分y_j是分解出来的信号，只保留独立成分y_j，其他独立成分置零， 是每个通道内所含的信号成分，由于核磁共振响应数据通道所含的核磁共振信号明显，正弦函数及余弦函数所在的通道基本不含核磁共振信号，因此选择核磁共振信号强的作为去工频谐波后的核磁共振信号，记为

f、针对第二组核磁共振响应数据，重复步骤a～e得到去除工频谐波的第二组核磁共振响应数据，记为

g、针对第三组核磁共振响应数据，重复步骤a～e得到去除工频谐波的第三组核磁共振响应数据，记为

3.按照权利要求1所述的基于ICA的核磁共振地下水探测信号噪声消除方法，其特征在于，步骤C包括以下步骤：

Ⅰ、将三组去除工频谐波噪声的数据组成观测信号，并对其进行预处理；

Ⅱ、利用快速独立分量分析算法对预处理后的观测信号进行分离，得到幅度有很大衰减的解混信号和解混矩阵A；

Ⅲ、将解混信号中噪声成分少、核磁共振信号明显的成分保留，其他成分屏蔽，利用ICA逆变换进行数据重构，恢复核磁共振信号幅度。