发明内容:
本发明的目的就在于针对上述现有技术的不足,结合ICA和核磁共振信号与噪声特性,提供一种基于ICA的核磁共振地下水探测信号噪声消除方法
本发明的目的是通过以下技术方案实现的:
一种基于ICA的核磁共振地下水探测信号噪声消除方法,包括以下步骤:
A、录入三组核磁共振响应数据;
B、利用ICA独立成分分析,依次去除每组数据核磁共振信号中心频率附近的工频谐波干扰,然后分别对每组处理后的数据进ICA逆变换后再进行数据重构,得到消除工频谐波干扰的数据;
C、将去除工频谐波噪声的三组数据作为观测信号,然后对其进行ICA处理,以消除剩余的随机噪声,并利用ICA逆变换进行数据重构,得到最终的消噪数据。
步骤A所述的三组,是指三个源信号。因为,ICA算法的使用前提是观测信号的个数大于等于源信号的个数,为了计算简单通常取相同个数,三组数据相对合理。核磁共振数据包含MRS信号和噪声,其中噪声又包括工频谐波干扰、随机噪声和尖峰干扰等。本发明仅针对工频谐波干扰和随机噪声,因此,源信号分为MRS信号、工频谐波干扰和随机噪声。
步骤B包括以下步骤:
a、对第一组核磁共振响应数据进行快速傅里叶变换,确定该核磁共振地下水探测信号中心频率附近的工频谐波干扰频率f1,f2,…,fn;
b、结合所确定频率,分别构造与工频谐波同频率,与核磁共振响应数据同长度的正弦函数、余弦函数,例如工频谐波频率为f1,则正弦函数为:sin(2πf1k),余弦函数为:cos(2πf1k),并与核磁共振响应数据构成观测信号;
c、对观测信号进行预处理;
d、利用快速独立分量分析算法对预处理后的观测信号进行分离,得到幅度有很大衰减的解混信号和解混矩阵w;
e、将步骤d分离出的工频谐波噪声屏蔽清零,利用ICA逆变换 进行数据重构,恢复核磁共振信号的幅度,其中独立成分yj是分解出来的信号,只保留独立成分yj,其他独立成分置零,是每个通道内所含的信号成分,由于核磁共振响应数据通道所含的核磁共振信号明显,正弦函数及余弦函数所在的通道基本不含核磁共振信号,因此选择核磁共振信号强的作为去工频谐波后的核磁共振信号,记为
f、针对第二组核磁共振响应数据,重复步骤a~e得到去除工频谐波的第二组核磁共振响应数据,记为
g、针对第三组核磁共振响应数据,重复步骤a~e得到去除工频谐波的第三组核磁共振响应数据,记为
步骤C包括以下步骤:
Ⅰ、将三组去除工频谐波噪声的数据组成观测信号,并对其进行预处理;
Ⅱ、利用快速独立分量分析算法对预处理后的观测信号进行分离,得到幅度有很大衰减的解混信号和解混矩阵A;
Ⅲ、将解混信号中噪声成分少、核磁共振信号明显的成分保留,其他成分屏蔽,利用ICA逆变换进行数据重构,恢复核磁共振信号幅度。
ICA算法包括以下步骤:
第一步、观测信号去均值,首先按照公式求解观测信号每一行的均值,然后利用公式对观测信号进行去均值处理,使数据中心化,满足零均值的假设;
第二步、对去均值后的数据进行白化处理,根据公式Cx=E{xxT}求观测信号的协方差矩阵,并求协方差矩阵的特征值d=diag(d0,d1,…,d2n),特征向量e=(e0,e1,…e2n),最后由公式x=vx=d1/2eTx对零均值观测信号进行白化处理,使数据具有单位方差,以消除数据各分量之间的相关性;
第三步、求解混矩阵w,首先令j=1,初始化解混向量w1,利用根据负熵最大化独立判据和牛顿迭代优化算法推导出递推式:
求w1,
利用公式对w1进行正交化,在根据wj=wj/||wj||标准化w1,当w1收敛时,第一个独立成分对应的解混向量w1求解完毕,判断w1是否收敛,如果w1不收敛重新利用公式求w1,对w1进行正交化和标准化,直到w1收敛,第一个独立成分对应的解混向量w1求解完毕,j=j+1,判断j≤m是否成立,如果j≤m,按照上述的步骤求出w2,直到j>m时,独立成分对应的所有的解混向量w求解完毕,输出解混矩阵w;
第四步、利用公式y=wTx求解独立成分y;
第五步、数据重构,由于利用ICA算法求出的输出信号y具有幅度的不确定性,与源信号相比初始幅度明显减小,将第四步分离出的噪声屏蔽清零,利用ICA逆变换 进行数据重构,恢复核磁共振信号的幅度,其中独立成分yj是分解出来的信号,只保留独立成分yj,其他独立成分置零,是每个通道内所含的信号成分。
有益效果:本发明对源信号和传输通道的先验知识没有要求,在试验过程中不需要铺设参考线圈,使工频噪声的消除效果不受参考线圈噪声采集情况的影响,在随机噪声消除的过程中不需要大量的数据进行叠加,操作简单、方便、效率高,由于采用了ICA算法,在分离噪声的过程中对信号的细节几乎没有破坏。
具体实施方式:
图1是基于ICA的核磁共振地下水探测信号工频谐波噪声消除方法的工作流程图。针对核磁共振地下水探测时受到的工频谐波的干扰,通过构造与工频谐波同频率的正弦函数、余弦函数,并和核磁共振响应信号组成观测信号,利用ICA算法对观测信号进行分离,并采用数据重构的方法解决ICA分离过程中产生的幅度衰减问题。
图2是利用ICA算法对去除工频噪声后的核磁共振地下水探测信号进行随机噪声消除的工作框图,将去工频后的三组核磁共振数据作为观测信号,用ICA算法对其进行处理,并采用数据重构的方法解决ICA分离过程中产生的幅度衰减问题,以压制核磁共振数据中的剩余噪声。
下面是基于ICA的核磁共振地下水探测信号噪声消除方法的具体步骤:
A、录入三组核磁共振响应数据,记为x0(k)、x1(k)、x2(k);
B、利用ICA算法依次去除x0(k)、x1(k)、x2(k)中核磁共振信号中心频率附近的工频谐波干扰,具体步骤如下:
a、对数据x0(k)进行快速傅里叶变换,确定该核磁共振地下水探测信号中心频率附近的工频谐波干扰频率f1,f2,…,fn;
b、分别构造与工频谐波同频率,与核磁共振响应数据同长度的正弦函数(x1sin=sin(2πf1k),x2sin=sin(2πf2k),…,xnsin=sin(2πfnk))、余弦函数(x1cos=cos(2πf1k),x2cos=cos(2πf2k),…,xncos=cos(2πfnk)),并与核磁共振响应数据构成观测信号x=[x0(k),x1sin(k),x1cos(k),x2sin(k),x2cos(k),…,xnsin(k),xncos(k)]T;
c、对观测信号进行预处理。首先按照公式依次求解观测信号每一行的均值,记为 然后利用公式对观测信号进行去均值处理,在根据公式Cx=E{xxT}求观测信号的协方差矩阵,并求协方差矩阵的特征值d=diag(d0,d1,…,d2n),特征向量e=(e0,e1,…e2n),最后由公式x=vx=d1/2eTx对零均值观测信号进行白化处理,得到预处理后的观测信号;
d、利用快速独立分量分析算法对预处理后的观测信号进行分离,得到幅度有很大衰减的核磁共振信号和工频谐波噪声和解混矩阵w,具体步骤如下:
i.确定观测信号的维数m=2n+1,即为独立分量的个数。令j=1;
ii.初始化解混向量wj;
iii.按照公式wj=wj/||wj||归一化解混向量;
iv.选择非线性函数
v.由公式 求解解混向量;
vi.根据公式正交化解混向量wj,并归一化wj;
vii.判断,如果wj收敛,j=j+1,否则返回步骤v;
viii.判断,如果j<m,返回步骤ii,如果j=m,输出解混矩阵w;
ix.由公式y=wTx估计独立分量。
e、将步骤d分离出的工频谐波噪声屏蔽清零,利用ICA逆变换 进行数据重构,恢复核磁共振信号的幅度,其中独立成分yj是分解出来的信号,只保留独立成分yj,其他独立成分置零,是每个通道内所含的信号成分,由于核磁共振响应数据通道所含的核磁共振信号明显,正弦函数及余弦函数所在的通道基本不含核磁共振信号,因此选择核磁共振信号强的作为去工频谐波后的核磁共振信号,记为
f、针对第二组核磁共振响应数据,重复步骤a~e得到去除工频谐波的第二组核磁共振响应数据,记为
g、针对第三组核磁共振响应数据,重复步骤a~e得到去除工频谐波的第三组核磁共振响应数据,记为
C、将去除工频谐波噪声的三组数据作为观测信号,然后其进行ICA处理,以消除剩余的随机噪声,并利用ICA逆变换进行数据重构,得到最终的消噪数据,具体步骤如下:
a、将三组去除工频谐波噪声的数据组成观测信号并对其进行预处理,首先按照公式依次求解观测信号每一行的均值,然后利用公式对观测信号进行去均值处理,在根据公式求观测信号的协方差矩阵,并求协方差矩阵的特征值d=diag(d0,d1,d3),特征向量e=(e0,e1,e3),最后由公式对零均值观测信号进行白化处理,得到预处理后的观测信号;
b、利用快速独立分量分析算法对预处理后的观测信号进行分离,得到幅度有很大衰减的解混信号和解混矩阵A,具体步骤如下:
i.确定观测信号的维数m=3,即为独立分量的个数。令i=1;
ii.初始化解混向量Ai;
iii.按照公式Ai=Ai/||Ai||归一化解混向量;
iv.选择非线性函数
v.由公式 求解解混向量;
vi.根据公式正交化解混向量Ai,并归一化Ai;
vii.判断,如果Ai收敛,i=i+1,否则返回步骤v;
viii.判断,如果i<m,返回步骤ii,如果j=m,输出解混矩阵A;
ix.按照公式y=ATx估计独立分量。
c、将解混信号中噪声成分少、核磁共振信号明显的成分保留,其他成分屏蔽,利用ICA逆变换 进行数据重构,恢复核磁共振信号幅度。
应用示例:
以吉林省烧锅镇核磁共振地下水探测为例:根据当地的地磁场强度计算得到核磁共振的拉莫尔频率为2326Hz,也就是核磁共振响应数据的中心频率为2326Hz,数据长度为16000,采样频率为66666.7Hz,在matlab环境下,基于ICA算法对核磁共振地下水探测实测数据中的工频谐波噪声进行处理。
具体步骤如下:
A、录入三组核磁共振响应数据,记为x0(k)、x1(k)、x2(k);
B、利用ICA算法依次去除x0(k)、x1(k)、x2(k)中核磁共振信号中心频率附近的工频谐波干扰,具体步骤如下:
a、对数据x0(k)进行快速傅里叶变换,确定该核磁共振地下水探测信号中心频率附近的工频谐波干扰频率f1=2300Hz,f2=2350Hz,f3=2450Hz;
b、分别构造与工频谐波同频率,与核磁共振响应数据同长度的正弦函数(x1sin=sin(2πf1k),x2sin=sin(2πf2k),x3sin=sin(2πf3k))、余弦函数(x1cos=cos(2πf1k),x2cos=cos(2πf2k),x3cos=cos(2πf3k)),并与核磁共振响应数据构成观测信号x=[x0(k),x1sin(k),x1cos(k),x2sin(k),x2cos(k),x3sin(k),x3cos(k)]T,观测信号如图4所示;
c、对观测信号进行预处理。首先按照公式依次求解观测信号每一行的均值,记为 然后利用公式对观测信号进行去均值处理,在根据公式Cx=E{xxT}求观测信号的协方差矩阵,并求协方差矩阵的特征值d=diag(d0,d1,d2,d3,d4,d5,d6),特征向量e=(e0,e1,e2,e3,e4,e5,e6),最后由公式x=vx=d1/2eTx对零均值观测信号进行白化处理,得到预处理后的观测信号;
d、利用快速独立分量分析算法对预处理后的观测信号进行分离,得到幅度有很大衰减的核磁共振信号和工频谐波噪声和解混矩阵w,分离后的独立分量估计如图5所示,具体步骤如下:
i.确定观测信号的维数m=7,即为独立分量的个数。令j=1;
ii.初始化解混向量wj;
iii.按照公式wj=wj/||wj||归一化解混向量;
iv.选择非线性函数
v.由公式 求解解混向量;
vi.根据公式正交化解混向量wj,并归一化wj;
vii.判断,如果wj收敛,j=j+1,否则返回步骤v;
viii.判断,如果j<7,返回步骤ii,如果j=7,输出解混矩阵w;
ix.按照公式y=wTx估计独立分量。
e、将步骤d分离出的工频谐波噪声屏蔽清零,利用ICA逆变换 进行数据重构,恢复核磁共振信号的幅度,其中独立成分y0是分解出来的信号,是每个通道内所含的信号成分,由于核磁共振响应数据通道所含的核磁共振信号明显,正弦函数及余弦函数所在的通道基本不含核磁共振信号,因此选择核磁共振信号强的作为去工频谐波后的核磁共振信号,记为重构后核磁共振信号如图6所示,消噪前后的频谱分析对比图如图7所示。
f、针对第二组核磁共振响应数据,重复步骤a~e得到去除工频谐波的第二组核磁共振响应数据,记为重构核磁共振信号如图8所示,消噪前后的频谱分析对比图如图9所示。
g、针对第三组核磁共振响应数据,重复步骤a~e得到去除工频谐波的第三组核磁共振响应数据,记为重构核磁共振信号如图10所示,消噪前后的频谱分析对比图如图11所示。
C、将去除工频谐波噪声的三组数据作为观测信号,并对其进行ICA处理,以消除剩余的随机噪声,得到最终的消噪数据,具体步骤如下:
a、将三组去除工频谐波噪声的数据组成观测信号观测信号如图12所示,并对其进行预处理,首先按照公式依次求解观测信号每一行的均值,然后利用公式对观测信号进行去均值处理,在根据公式求观测信号的协方差矩阵,并求协方差矩阵的特征值d=diag(d0,d1,d3),特征向量e=(e0,e1,e3),最后由公式对零均值观测信号进行白化处理,得到预处理后的观测信号;
b、利用快速独立分量分析算法对预处理后的观测信号进行分离,得到幅度有很大衰减的解混信号和解混矩阵A,分离后的独立分量估计如图13所示,具体步骤如下:
i.确定观测信号的维数m=3,即为独立分量的个数,令i=1;
ii.初始化解混向量Ai;
iii.按照公式Ai=Ai/||Ai||归一化解混向量;
iv.选择非线性函数
v.根据公式 求解解混向量;
vi.由公式正交化解混向量Ai,并归一化Ai;
vii.判断,如果Ai收敛,i=i+1,否则返回步骤v;
viii.判断,如果i<3,返回步骤ii,j=3,输出解混矩阵A;
ix.按照公式y=ATx估计独立分量。
c、将解混信号中噪声成分少、核磁共振信号明显的成分保留,其他成分屏蔽,利用ICA逆变换 进行数据重构,恢复核磁共振信号幅度,重构后的核磁共振信号如图14所示。
基于ICA算法的核磁共振地下水探测信号工频谐波噪声消除方法,不仅使核磁共振响应数据中的工频谐波被很好地消除,而且没有破坏核磁共振信号,对其他的噪声也有很好的保留,这样不会干扰后期对其他噪声的处理,基于ICA算法的核磁共振地下水探测信号随机噪声消除方法,很好地消除了随机噪声,消噪后的核磁共振信号的时域波形图很接近理想的呈e指数衰减的核磁共振信号,ICA算法对实测数据的结果分别如图14、图16。其中图14是两次消噪后核磁共振信号的时域波形图,而图15是消噪前核磁共振实测数据的时域波形图,对比图14、图15的结果,可以发现ICA消噪后的核磁共振信号e指数衰减的趋势更明显,图16是消噪后核磁共振信号与消噪前的其中一个实测数据频谱分析对比图,可以发现ICA消噪后核磁共振信号的频谱噪声成分更少。