一种基于互相关的核磁共振全波信号噪声滤除方法
一种基于互相关的核磁共振全波信号噪声滤除方法
【技术领域】
[0001] 本发明涉及一种核磁共振测深（Magnetic Resonance Sounding, MRS)信号噪声滤 除方法，具体是基于互相关噪声研制原理的核磁共振全波信号噪声的滤除方法。
【背景技术】
[0002] 核磁共振地下水探测利用人工激发的电磁场使地下水中的氢核形成宏观磁矩，这 一宏观磁矩在地磁场中产生旋进运动，使用线圈接收宏观磁矩进动产生的电磁信号，该电 磁信号的幅度随时间按指数形式衰减，成为自由感应衰减（Free Induction Decay，FID)信 号，通过对FID进行分析可以探测地下是否存在地下水。核磁共振找水是一种快速、直接的 地下水探测方法。FID信号的表达式为：
[0004] 上式中，表示拉莫尔角频率，由当地的地磁场强度决定，一般为1000-3000HZ ; FID信号的三个特征参数E。（初始振幅）、T/(平均弛豫时间）和9。（初始相位）分别与含 水层的分布、厚度、平均含水率、渗透率和导电性等信息有着密切的关系，准确测量FID信 号是进行核磁共振地下水探测的关键，测量FID信号时需要在野外勘探过程中将线圈接收 信号的时域波形记录下来，该时域波形成为核磁共振全波信号。进行地下水探测时，核磁共 振信号非常微弱，一般为nV级，而周围电磁场噪声却非常强烈，如工频及其谐波噪声、随机 噪声和尖峰噪声，即使在信号检测电路中使用了带通滤波器，核磁共振全波信号的信噪比 依然很低，一般低于OdB，这严重限制了核磁共振地下水探测的应用。
[0005] 目前，核磁共振全波信号一般采用数据叠加方法，该方法能抑制大部分的随机噪 声，但对工频噪声和尖峰噪声抑制效果差，且需要多次采集核磁共振信号，测量效率低。 CN103823244A公开了一种磁共振三分量消噪装置，该发明基于参考通道噪声与主通道含 噪信号中噪声的相关性来实现信噪分离，但由于某些噪声的无规律性和不稳定性以及混合 机制的不确定性，限制了算法的应用，而且增加参考线圈和参考通道使仪器系统庞大复杂。 CN104459809A公开了一种基于独立成分分析的全波核磁共振信号噪声滤除方法，首先利用 核磁共振测深探水仪采集FID信号，通过频谱分析获得采集信号中含有的工频谐波干扰频 率，采用数字正交法构造输入通道信号解决欠定盲源分离问题；然后，将构造的输入通道信 号与采集的全波信号一并作为输入信号进行独立成分分析，得到分离的FID信号；最后采 用频谱校正法解决分离后MRS信号幅值不定问题，进而提取去噪后的FID信号。该发明主要 针对核磁共振全波信号中工频谐波干扰或某一单频干扰，不能滤除随机噪声和尖峰噪声， 该方法算法过程复杂，计算量大，而且输出结果不稳定。
【发明内容】
[0006] 本发明的目的就在于针对上述现有技术的不足，提供一种基于互相关的核磁共振 全波信号噪声滤除方法。
[0007] 本发明的目的是通过以下方式实现的：
[0008] -种基于互相关的核磁共振全波信号噪声滤除方法，包括以下步骤：
[0009] A、读取核磁共振接收机野外采集的全波数据s (k)及其采样率fs;
[0010] B、使用快速傅里叶变换求取全波数据的频谱S(?);
[0011] C、读取频谱S (?)中FID信号的频率f0所对应的相位0 0 ;
[0012] D、读取频谱S (?)中，拉莫尔频率f0左右两侧的两个工频谐波频率n和f2,及其 所对应的幅度al，a2,相位0 1和0 2 ;
[0013] E、重构两个工频谐波噪声的波形nl (k)和n2 (k);
[0014] F、在全波数据s(k)中减去nl (k)和n2(k)，得到x(k);
[0015] G、将x(k)与参考信号cos(2itkfO/fs)进行互相关运算，得到R(m);
[0016] H、利用希尔伯特变换获取R(m)的包络信号，En(m);
[0017]I、将包络信号En(m)分为长度相等的两段，并分别使用双指数函数拟合包络曲 线；
[0018] J、使用双指数函数的拟合结果重构包络信号，Env (m) = a ? exp (bmfO/ fs)+c ? exp(dmfO/fs)，a、b、c、d 为拟合参数；
[0019] K、构造不含噪声的互相关波形 rf (m) = Env (m) ? cos (2itmf0/fs+0 0);
[0020] 使用参考信号C〇S(2 3ikf0/fS)对rf(m)进行解卷积运算，得到的解卷积结果 ff(k)即为滤除噪声的核磁共振全波信号。
[0021] 有益效果：本发明提供的基于互相关核磁共振全波信号噪声滤除方法，是利用噪 声与拉莫尔频率的正弦信号不相关，而fid幅度衰减正弦信号与拉莫尔频率的正弦信号具 有相关性的特点，通过互相关运算滤除噪声，然后拟合互相关波形的包络，并重构不含噪声 的互相关波形，最后利用解卷积算法提取核磁共振全波数据中的FID信号。该方法运算数 据计算量小，可以同时压制工频及其谐波噪声、随机噪声和尖峰噪声，明显提高核磁共振全 波数据信噪比，有利于扩展核磁共振找水仪的应用范围和探测深度，而且该发明不需要额 外的参考消噪装置，节约了成本，一次采集的全波数据通过互相关即可达到理想的滤波效 果，节省测量时间。
【附图说明】
[0022] 附图1是基于互相关的核磁共振全波信号噪声滤除方法的流程图。
[0023] 附图2是野外采集的不含噪声的核磁共振全波数据的波形图。
[0024] 附图3是野外采集的含有噪声的核磁共振全波数据的波形图。
[0025] 附图4是含有噪声的核磁共振全波数据的频谱图。
[0026] 附图5是核磁共振全波数据与参考信号进行互相关后的波形图。
[0027] 附图6是互相关波形的包络拟合曲线。
[0028] 附图7是重构后得到的不含噪声的互相关波形图。
[0029]附图8是解卷积得到的滤除噪声的核磁共振全波数据波形图。
【具体实施方式】
[0030] 下面结合附图和实施例作进一步详细说明：
[0031] 如图1所示，一种基于互相关的核磁共振全波信号噪声滤除方法，包括以下步骤：
[0032] A、读取核磁共振接收机野外采集的全波数据，核磁共振找水仪接收机野外工作时 以恒定的采样率f s采集数据，读取接收机记录的全波数据s (k)，k = 1，2, 3…N ;
[0033] B、使用快速傅里叶变换求取全波数据的频谱S(?)，在Matlab中使用fft(s，N)命 令可得到全波的频谱S(?)，S(?)包含N个频率所对应的复数，频率间隔为fs/N;
[0034]C、读取频谱S (?)中FID信号的频率fQ所对应的相位0。，S (?)包含N个频率所 对应的复数代表了该频率的幅度和相位信息，即
[0036] 上式中，abs表示取模运算符，arctan为反正切运算符，im为取虚部运算符，re为 取实部运算符；
[0037] D、读取频谱S(?)中，拉莫尔频率&左右两侧的两个工频谐波频率f :和f 2，全波 数据与参考信号作互相关运算时，如果在频谱有上接近拉莫尔频率的周期信号，该周期信 号会影响互相关运算的结果，这是由于当周期噪声的频率接近信号频率时，周期噪声与FID 信号的相关性增强，互相关噪声抑制效果变差。因此，需要将频率最接近拉莫尔频率的两个 工频谐波滤除掉。读取这两个工频谐波噪声的幅度apa 2，相位^和9 2;
[0038] E、重构两个工频谐波噪声的波形ni (k)和n2 (k)，由于这两个工频谐波为频率已知 的周期信号，当其幅度和相位已知时，其波形就可唯一确定，可使用如下公式重构这两个工 频谐波噪声的波形：
[0040] F、在全波数据s (k)中减去(k)和n2 (k)，得到x (k)
[0041] x (k) = s (k) -rii (k) _n2 (k)
[0042] 得到的x(k)不含有频率最接近拉莫尔频率的两个工频谐波噪声，但依然含有其 他频率的工频噪声、随机噪声和尖峰噪声；
[0043] G、将x(k)与参考信号cos(2 Jr 进行互相关运算，得到R(k)，计算公式为：
[0045] 相关检测是一种噪声抑制的有效方法，它基于信号与噪声的统计特性进行检测， 互相关函数是两个时域信号相似性的一种度量，被检测信号与参考信号具有相关性，而与 工频及其谐波噪声、随机噪声和尖峰噪声的相关性较差或没有相关性，利用这一差异可 把被检测信号与噪声区分开了。虽然互相关检测一般应用于周期的正弦信号，但是核磁 共振的FID信号是幅度按指数规律衰减正弦信号，与频率为拉莫尔频率的参考正弦信号 cos (2 31 具有较强相关性，而与随机噪声和尖峰噪声没有相关性，与频率不等于拉莫 尔频率的周期噪声相关性较差，因此，互相关检测适用于核磁共振全波数据噪声滤除。两个 长度均为N的数据进行互相关运算后得到的互相关数组长度为2N-1 ;
[0046] H、利用希尔伯特变换获取R(m)的包络信号En(m)，互相关数据R(m)的希尔伯特变 换为：
[0048] 信号经希尔伯特变换后，在频域各频率分量的幅度保持 不变，但相位将出现90° 相移。得到解解析信号为：z(m) =R(m)+jR' (m)
[0049] 解析信号z (k)的瞬时振幅
'即为互相关R(m)波形的包
[0050]络 En(m);
[0051]I、将 En(m)分为两段，每段长度均为 N，即 Enjk) = En(l:N)，En2(k)= En (N: 2N-1)，分别使用双指数函数拟合Eni (k)和En2 (k)，采用方法为最小二乘拟合，即
[0054] ai、VVcpdi为拟合得到的参数；
[0055] J、使用双指数函数的拟合结果重构包络信号，Env(1:N)= ?exp (hkf。/ fJ+Ci ? exp^kf。/;^)，Env(N:2N_l) = a2 ? expO^kfo/f^+q ? 6叉口((121^。/;^)，拟合后的包 络信号将不含有噪声；
[0056]K、构造不含噪声的互相关波形 rf (m) = Env(m) ? cos(2 3imfQ/fs+0 Q)，Env(m)为 拟合得到的包络，cos (2 Jr 0 J为FID信号相位相同的正弦信号，其频率为拉莫尔频 率；
[0057] 使用参考信号cosOnkf^/f；)对1^〇〇进行解卷积运算，互相关函数可表示为卷 积运算：
[0058] R12= fj(t) *f2 (t),
[0059] 表示卷积运算，可见将f2 (t)反褶（变量取符号）与⑴做卷积即得到(t) 与4的互相关函数，同时，由于参考信号的反褶为其本身，R(m)可以看作 x(k)与参考信号cosQjikfVf；)的卷积运算，重构的互相关波形rf(m)是不含噪声的核 磁全波数据与参考信号的卷积。因此，通过反卷积可以得到不含噪声的核磁全波数据。在 Matlab是使用解卷积命令deconv即可得到ff (k) :ff = deconv (rf, cos (2 Jr fQt))，得到的 解卷积结果ff(k)即为滤除噪声的核磁共振全波信号。
[0060] 野外采集的核磁共振全波数据含有FID的波形为：
[0062] 上式中，FID信号的初始幅度为100nV，平均弛豫时间T/为300ms，拉莫尔频率f0 为2326Hz，初始相位为30°，接收机采样率为100kHz，信号采集时间为500ms，采样点数为 50000。不含噪声的FID信号波形见图2所示。
[0063] 采集到的核磁共振全波数据中除了微弱的FID信号外，还包含工频及其谐波噪 声、随机噪声和尖峰噪声，噪声的表达式为：
[0064] n (k) = lCT6square (t，50)+wag (1，N，-140)+lCT6pulstran (t，0? 05: 0? 05: 0? 5, ' gau spuls'）
[0065] 上式中，Matlab中square用于产生50Hz方波，来模拟工频工频及其谐波噪声，方 波的幅度为1 U V ;wgn(l，N，-140)用于产生标准差为100nV的零均值随机噪声；pulstran 函数用于产生尖峰噪声，尖峰噪声的幅度为1 UV ;这些噪声为加性噪声叠加在FID信号上， 核磁共振全波数据s(k)的表达式为：
[0066] s (k) = f (k) +n (k)
[0067] 该核磁共振全波数据的信噪比为_30dB。包含噪声的核磁共振全波数据波形如图 3所示。
[0068] 在Matlab中使用函数fft函数对全波数据x(k)进行快速傅里叶变换，得到的频 谱图如图4所示，频谱数据S(?)包含50000个频率所对应的复数，频率间隔为2Hz。
[0069] 从频谱数组S(?)中找到第1164个元素为，对此复数取相角得到0^*31°。
[0070] 根据拉莫尔频率2326Hz，得到最接近拉莫尔频率的两个工频谐波频率分别为 2250Hz和2350Hz，从频谱数组S(?)中找到第1126个元素为1. 3216e-09-l. 3472e-08i，对 此复数取模得到叫为27.731^，；从频谱数组S(?)中找到第1176个元素为 6.5767e-10-1.3849e-08i，对此复数取模得到a 2S 27.07nV，0 ^为-1.4°。重构得到的这 两个工频谐波波形为：
[0072] 从核磁全波数据s (k)中减去这两个工频谐波噪声得到x (k)，计算式如下：
[0073] x (k) = s (k) -rii (k) -n2 (k)
[0074] 由于减去的这两个谐波噪声相对总体噪声而言很小，核磁共振全波数据的波形几 乎不会发生明显变化。将x(k)与参考信号cos (2 31 按照如下公式进行互相关运算， 在Matlab中使用xcorr命令即可得到两个数组的互相关函数，得到的R(m)波形如图5所 示，可见互相关运算后的波形与原始的FID信号存在差异，波形可分为两段：前一段呈现幅 度按指数规律增加的正弦信号，后一段呈现幅度按指数规律衰减的正弦信号。
[0075] 利用希尔伯特变换获取R(m)的包络信号En(m)，将En(m)分为两段，每段长度均 为 50000,即 E& (k) = En (1:50000)，En2 (k) = En (50000:99999)，分别使用双指数函数拟合 Eni (k)和En2 (k)，采用方法为最小二乘拟合，得到的拟合结果为：
[0076] a:= 0. 001364, b := 0. 1022, 〇!=-〇. 001367, d := -3. 625 ；a 2= -〇. 0001759, b 2 =0? 4319, c2= 0? 001402, d2= -3. 628。
[0077] 使用双指数函数的拟合结果重构包络信号，Env (1 :N) = ? exp fJ+Ci ? exp^kf。/;^)，Env(N:2N_l) = a2 ? expO^kfo/f^+q ? 6叉口((121^。/;^)，拟合后的包 络信号将不含有噪声，图6是互相关波形的包络拟合曲线，包络拟合曲线分为两段。
[0078] 构造不含噪声的互相关波形 rf (m) = Env(m) ? cos (2 Jr mfQ/fs+9 Q)，Env(m)为拟 合得到的包络，图7是重构后得到的不含噪声的互相关波形图，重构后的互相关波形比原 互相关波形光滑。
[0079] 使用参考彳目号cos (2 Jr kfd/fj对rf (k)进彳丁解卷积运算，在Matlab是使用解卷积 命令deconv即可得到ff (k) :ff = deconv(rf, cos(2 Jr fQt))，得到的解卷积结果ff (k)即 为滤除噪声的核磁共振全波信号，图8是解卷积得到的滤除噪声的核磁共振全波数据波形 图，对比图2和图8发现二者波形一致，说明基于互相关的核磁共振全波信号噪声滤除方法 能够恢复淹没在噪声中的FID信号。
【主权项】
1. 一种基于互相关的核磁共振全波信号噪声滤除方法，其特征在于，包括以下步骤： A、 读取核磁共振接收机野外采集的全波数据s(k)及其采样率fs; B、 使用快速傅里叶变换求取全波数据的频谱S(U); C、 读取频谱S(?)中FID信号的频率fO所对应的相位0O; D、 读取频谱S(?)中，拉莫尔频率fO左右两侧的两个工频谐波频率fl和f2，及其所对 应的幅度al，a2,相位0 1和0 2 ; E、 重构两个工频谐波噪声的波形nl(k)和n2 (k); F、 在全波数据s(k)中减去nl(k)和n2 (k)，得到X(k); G、 将x(k)与参考信号C〇S(2 3ikfO/fS)进行互相关运算，得到R(m); H、 利用希尔伯特变换获取R(m)的包络信号，En(m); I、 将包络信号En(m)分为长度相等的两段，并分别使用双指数函数拟合包络曲线； J、 使用双指数函数的拟合结果重构包络信号，Env(m) =a?exp(bmfO/ fs)+c?exp(dmfO/fs)，a、b、c、d为拟合参数； K、 构造不含噪声的互相关波形rf(m) =Env(m) ?cosQjimfO/fs+QO);使用参考信号 cos(2 3ikf0/fs)对rf(m)进行解卷积运算，得到的解卷积结果ff(k)即为滤除噪声的核磁 共振全波信号。
【专利摘要】本发明涉及一种基于互相关的核磁共振全波信号噪声滤除方法，是利用噪声与拉莫尔频率的正弦信号不相关，而FID幅度衰减正弦信号与拉莫尔频率的正弦信号具有相关性的特点，通过互相关运算滤除噪声，然后拟合互相关波形的包络，并重构不含噪声的互相关波形，最后利用解卷积算法提取核磁共振全波数据中的FID信号。该方法运算数据计算量小，可以同时压制工频及其谐波噪声、随机噪声和尖峰噪声，明显提高核磁共振全波数据信噪比，有利于扩展核磁共振找水仪的应用范围和探测深度，而且该发明不需要额外的参考消噪装置，节约了成本，一次采集的全波数据通过互相关即可达到理想的滤波效果，节省测量时间。
【IPC分类】G01V3/14
【公开号】CN104898172
【申请号】CN201510256496
【发明人】林君, 刘立超
【申请人】吉林大学
【公开日】2015年9月9日
【申请日】2015年5月19日
【技术领域】
[0001] 本发明涉及一种核磁共振测深（Magnetic Resonance Sounding, MRS)信号噪声滤 除方法，具体是基于互相关噪声研制原理的核磁共振全波信号噪声的滤除方法。
【背景技术】
[0002] 核磁共振地下水探测利用人工激发的电磁场使地下水中的氢核形成宏观磁矩，这 一宏观磁矩在地磁场中产生旋进运动，使用线圈接收宏观磁矩进动产生的电磁信号，该电 磁信号的幅度随时间按指数形式衰减，成为自由感应衰减（Free Induction Decay，FID)信 号，通过对FID进行分析可以探测地下是否存在地下水。核磁共振找水是一种快速、直接的 地下水探测方法。FID信号的表达式为：
[0004] 上式中，表示拉莫尔角频率，由当地的地磁场强度决定，一般为1000-3000HZ ; FID信号的三个特征参数E。（初始振幅）、T/(平均弛豫时间）和9。（初始相位）分别与含 水层的分布、厚度、平均含水率、渗透率和导电性等信息有着密切的关系，准确测量FID信 号是进行核磁共振地下水探测的关键，测量FID信号时需要在野外勘探过程中将线圈接收 信号的时域波形记录下来，该时域波形成为核磁共振全波信号。进行地下水探测时，核磁共 振信号非常微弱，一般为nV级，而周围电磁场噪声却非常强烈，如工频及其谐波噪声、随机 噪声和尖峰噪声，即使在信号检测电路中使用了带通滤波器，核磁共振全波信号的信噪比 依然很低，一般低于OdB，这严重限制了核磁共振地下水探测的应用。
[0005] 目前，核磁共振全波信号一般采用数据叠加方法，该方法能抑制大部分的随机噪 声，但对工频噪声和尖峰噪声抑制效果差，且需要多次采集核磁共振信号，测量效率低。 CN103823244A公开了一种磁共振三分量消噪装置，该发明基于参考通道噪声与主通道含 噪信号中噪声的相关性来实现信噪分离，但由于某些噪声的无规律性和不稳定性以及混合 机制的不确定性，限制了算法的应用，而且增加参考线圈和参考通道使仪器系统庞大复杂。 CN104459809A公开了一种基于独立成分分析的全波核磁共振信号噪声滤除方法，首先利用 核磁共振测深探水仪采集FID信号，通过频谱分析获得采集信号中含有的工频谐波干扰频 率，采用数字正交法构造输入通道信号解决欠定盲源分离问题；然后，将构造的输入通道信 号与采集的全波信号一并作为输入信号进行独立成分分析，得到分离的FID信号；最后采 用频谱校正法解决分离后MRS信号幅值不定问题，进而提取去噪后的FID信号。该发明主要 针对核磁共振全波信号中工频谐波干扰或某一单频干扰，不能滤除随机噪声和尖峰噪声， 该方法算法过程复杂，计算量大，而且输出结果不稳定。
【发明内容】
[0006] 本发明的目的就在于针对上述现有技术的不足，提供一种基于互相关的核磁共振 全波信号噪声滤除方法。
[0007] 本发明的目的是通过以下方式实现的：
[0008] -种基于互相关的核磁共振全波信号噪声滤除方法，包括以下步骤：
[0009] A、读取核磁共振接收机野外采集的全波数据s (k)及其采样率fs;
[0010] B、使用快速傅里叶变换求取全波数据的频谱S(?);
[0011] C、读取频谱S (?)中FID信号的频率f0所对应的相位0 0 ;
[0012] D、读取频谱S (?)中，拉莫尔频率f0左右两侧的两个工频谐波频率n和f2,及其 所对应的幅度al，a2,相位0 1和0 2 ;
[0013] E、重构两个工频谐波噪声的波形nl (k)和n2 (k);
[0014] F、在全波数据s(k)中减去nl (k)和n2(k)，得到x(k);
[0015] G、将x(k)与参考信号cos(2itkfO/fs)进行互相关运算，得到R(m);
[0016] H、利用希尔伯特变换获取R(m)的包络信号，En(m);
[0017]I、将包络信号En(m)分为长度相等的两段，并分别使用双指数函数拟合包络曲 线；
[0018] J、使用双指数函数的拟合结果重构包络信号，Env (m) = a ? exp (bmfO/ fs)+c ? exp(dmfO/fs)，a、b、c、d 为拟合参数；
[0019] K、构造不含噪声的互相关波形 rf (m) = Env (m) ? cos (2itmf0/fs+0 0);
[0020] 使用参考信号C〇S(2 3ikf0/fS)对rf(m)进行解卷积运算，得到的解卷积结果 ff(k)即为滤除噪声的核磁共振全波信号。
[0021] 有益效果：本发明提供的基于互相关核磁共振全波信号噪声滤除方法，是利用噪 声与拉莫尔频率的正弦信号不相关，而fid幅度衰减正弦信号与拉莫尔频率的正弦信号具 有相关性的特点，通过互相关运算滤除噪声，然后拟合互相关波形的包络，并重构不含噪声 的互相关波形，最后利用解卷积算法提取核磁共振全波数据中的FID信号。该方法运算数 据计算量小，可以同时压制工频及其谐波噪声、随机噪声和尖峰噪声，明显提高核磁共振全 波数据信噪比，有利于扩展核磁共振找水仪的应用范围和探测深度，而且该发明不需要额 外的参考消噪装置，节约了成本，一次采集的全波数据通过互相关即可达到理想的滤波效 果，节省测量时间。
【附图说明】
[0022] 附图1是基于互相关的核磁共振全波信号噪声滤除方法的流程图。
[0023] 附图2是野外采集的不含噪声的核磁共振全波数据的波形图。
[0024] 附图3是野外采集的含有噪声的核磁共振全波数据的波形图。
[0025] 附图4是含有噪声的核磁共振全波数据的频谱图。
[0026] 附图5是核磁共振全波数据与参考信号进行互相关后的波形图。
[0027] 附图6是互相关波形的包络拟合曲线。
[0028] 附图7是重构后得到的不含噪声的互相关波形图。
[0029]附图8是解卷积得到的滤除噪声的核磁共振全波数据波形图。
【具体实施方式】
[0030] 下面结合附图和实施例作进一步详细说明：
[0031] 如图1所示，一种基于互相关的核磁共振全波信号噪声滤除方法，包括以下步骤：
[0032] A、读取核磁共振接收机野外采集的全波数据，核磁共振找水仪接收机野外工作时 以恒定的采样率f s采集数据，读取接收机记录的全波数据s (k)，k = 1，2, 3…N ;
[0033] B、使用快速傅里叶变换求取全波数据的频谱S(?)，在Matlab中使用fft(s，N)命 令可得到全波的频谱S(?)，S(?)包含N个频率所对应的复数，频率间隔为fs/N;
[0034]C、读取频谱S (?)中FID信号的频率fQ所对应的相位0。，S (?)包含N个频率所 对应的复数代表了该频率的幅度和相位信息，即
[0036] 上式中，abs表示取模运算符，arctan为反正切运算符，im为取虚部运算符，re为 取实部运算符；
[0037] D、读取频谱S(?)中，拉莫尔频率&左右两侧的两个工频谐波频率f :和f 2，全波 数据与参考信号作互相关运算时，如果在频谱有上接近拉莫尔频率的周期信号，该周期信 号会影响互相关运算的结果，这是由于当周期噪声的频率接近信号频率时，周期噪声与FID 信号的相关性增强，互相关噪声抑制效果变差。因此，需要将频率最接近拉莫尔频率的两个 工频谐波滤除掉。读取这两个工频谐波噪声的幅度apa 2，相位^和9 2;
[0038] E、重构两个工频谐波噪声的波形ni (k)和n2 (k)，由于这两个工频谐波为频率已知 的周期信号，当其幅度和相位已知时，其波形就可唯一确定，可使用如下公式重构这两个工 频谐波噪声的波形：
[0040] F、在全波数据s (k)中减去(k)和n2 (k)，得到x (k)
[0041] x (k) = s (k) -rii (k) _n2 (k)
[0042] 得到的x(k)不含有频率最接近拉莫尔频率的两个工频谐波噪声，但依然含有其 他频率的工频噪声、随机噪声和尖峰噪声；
[0043] G、将x(k)与参考信号cos(2 Jr 进行互相关运算，得到R(k)，计算公式为：
[0045] 相关检测是一种噪声抑制的有效方法，它基于信号与噪声的统计特性进行检测， 互相关函数是两个时域信号相似性的一种度量，被检测信号与参考信号具有相关性，而与 工频及其谐波噪声、随机噪声和尖峰噪声的相关性较差或没有相关性，利用这一差异可 把被检测信号与噪声区分开了。虽然互相关检测一般应用于周期的正弦信号，但是核磁 共振的FID信号是幅度按指数规律衰减正弦信号，与频率为拉莫尔频率的参考正弦信号 cos (2 31 具有较强相关性，而与随机噪声和尖峰噪声没有相关性，与频率不等于拉莫 尔频率的周期噪声相关性较差，因此，互相关检测适用于核磁共振全波数据噪声滤除。两个 长度均为N的数据进行互相关运算后得到的互相关数组长度为2N-1 ;
[0046] H、利用希尔伯特变换获取R(m)的包络信号En(m)，互相关数据R(m)的希尔伯特变 换为：
[0048] 信号经希尔伯特变换后，在频域各频率分量的幅度保持 不变，但相位将出现90° 相移。得到解解析信号为：z(m) =R(m)+jR' (m)
[0049] 解析信号z (k)的瞬时振幅
'即为互相关R(m)波形的包
[0050]络 En(m);
[0051]I、将 En(m)分为两段，每段长度均为 N，即 Enjk) = En(l:N)，En2(k)= En (N: 2N-1)，分别使用双指数函数拟合Eni (k)和En2 (k)，采用方法为最小二乘拟合，即
[0054] ai、VVcpdi为拟合得到的参数；
[0055] J、使用双指数函数的拟合结果重构包络信号，Env(1:N)= ?exp (hkf。/ fJ+Ci ? exp^kf。/;^)，Env(N:2N_l) = a2 ? expO^kfo/f^+q ? 6叉口((121^。/;^)，拟合后的包 络信号将不含有噪声；
[0056]K、构造不含噪声的互相关波形 rf (m) = Env(m) ? cos(2 3imfQ/fs+0 Q)，Env(m)为 拟合得到的包络，cos (2 Jr 0 J为FID信号相位相同的正弦信号，其频率为拉莫尔频 率；
[0057] 使用参考信号cosOnkf^/f；)对1^〇〇进行解卷积运算，互相关函数可表示为卷 积运算：
[0058] R12= fj(t) *f2 (t),
[0059] 表示卷积运算，可见将f2 (t)反褶（变量取符号）与⑴做卷积即得到(t) 与4的互相关函数，同时，由于参考信号的反褶为其本身，R(m)可以看作 x(k)与参考信号cosQjikfVf；)的卷积运算，重构的互相关波形rf(m)是不含噪声的核 磁全波数据与参考信号的卷积。因此，通过反卷积可以得到不含噪声的核磁全波数据。在 Matlab是使用解卷积命令deconv即可得到ff (k) :ff = deconv (rf, cos (2 Jr fQt))，得到的 解卷积结果ff(k)即为滤除噪声的核磁共振全波信号。
[0060] 野外采集的核磁共振全波数据含有FID的波形为：
[0062] 上式中，FID信号的初始幅度为100nV，平均弛豫时间T/为300ms，拉莫尔频率f0 为2326Hz，初始相位为30°，接收机采样率为100kHz，信号采集时间为500ms，采样点数为 50000。不含噪声的FID信号波形见图2所示。
[0063] 采集到的核磁共振全波数据中除了微弱的FID信号外，还包含工频及其谐波噪 声、随机噪声和尖峰噪声，噪声的表达式为：
[0064] n (k) = lCT6square (t，50)+wag (1，N，-140)+lCT6pulstran (t，0? 05: 0? 05: 0? 5, ' gau spuls'）
[0065] 上式中，Matlab中square用于产生50Hz方波，来模拟工频工频及其谐波噪声，方 波的幅度为1 U V ;wgn(l，N，-140)用于产生标准差为100nV的零均值随机噪声；pulstran 函数用于产生尖峰噪声，尖峰噪声的幅度为1 UV ;这些噪声为加性噪声叠加在FID信号上， 核磁共振全波数据s(k)的表达式为：
[0066] s (k) = f (k) +n (k)
[0067] 该核磁共振全波数据的信噪比为_30dB。包含噪声的核磁共振全波数据波形如图 3所示。
[0068] 在Matlab中使用函数fft函数对全波数据x(k)进行快速傅里叶变换，得到的频 谱图如图4所示，频谱数据S(?)包含50000个频率所对应的复数，频率间隔为2Hz。
[0069] 从频谱数组S(?)中找到第1164个元素为，对此复数取相角得到0^*31°。
[0070] 根据拉莫尔频率2326Hz，得到最接近拉莫尔频率的两个工频谐波频率分别为 2250Hz和2350Hz，从频谱数组S(?)中找到第1126个元素为1. 3216e-09-l. 3472e-08i，对 此复数取模得到叫为27.731^，；从频谱数组S(?)中找到第1176个元素为 6.5767e-10-1.3849e-08i，对此复数取模得到a 2S 27.07nV，0 ^为-1.4°。重构得到的这 两个工频谐波波形为：
[0072] 从核磁全波数据s (k)中减去这两个工频谐波噪声得到x (k)，计算式如下：
[0073] x (k) = s (k) -rii (k) -n2 (k)
[0074] 由于减去的这两个谐波噪声相对总体噪声而言很小，核磁共振全波数据的波形几 乎不会发生明显变化。将x(k)与参考信号cos (2 31 按照如下公式进行互相关运算， 在Matlab中使用xcorr命令即可得到两个数组的互相关函数，得到的R(m)波形如图5所 示，可见互相关运算后的波形与原始的FID信号存在差异，波形可分为两段：前一段呈现幅 度按指数规律增加的正弦信号，后一段呈现幅度按指数规律衰减的正弦信号。
[0075] 利用希尔伯特变换获取R(m)的包络信号En(m)，将En(m)分为两段，每段长度均 为 50000,即 E& (k) = En (1:50000)，En2 (k) = En (50000:99999)，分别使用双指数函数拟合 Eni (k)和En2 (k)，采用方法为最小二乘拟合，得到的拟合结果为：
[0076] a:= 0. 001364, b := 0. 1022, 〇!=-〇. 001367, d := -3. 625 ；a 2= -〇. 0001759, b 2 =0? 4319, c2= 0? 001402, d2= -3. 628。
[0077] 使用双指数函数的拟合结果重构包络信号，Env (1 :N) = ? exp fJ+Ci ? exp^kf。/;^)，Env(N:2N_l) = a2 ? expO^kfo/f^+q ? 6叉口((121^。/;^)，拟合后的包 络信号将不含有噪声，图6是互相关波形的包络拟合曲线，包络拟合曲线分为两段。
[0078] 构造不含噪声的互相关波形 rf (m) = Env(m) ? cos (2 Jr mfQ/fs+9 Q)，Env(m)为拟 合得到的包络，图7是重构后得到的不含噪声的互相关波形图，重构后的互相关波形比原 互相关波形光滑。
[0079] 使用参考彳目号cos (2 Jr kfd/fj对rf (k)进彳丁解卷积运算，在Matlab是使用解卷积 命令deconv即可得到ff (k) :ff = deconv(rf, cos(2 Jr fQt))，得到的解卷积结果ff (k)即 为滤除噪声的核磁共振全波信号，图8是解卷积得到的滤除噪声的核磁共振全波数据波形 图，对比图2和图8发现二者波形一致，说明基于互相关的核磁共振全波信号噪声滤除方法 能够恢复淹没在噪声中的FID信号。
【主权项】
1. 一种基于互相关的核磁共振全波信号噪声滤除方法，其特征在于，包括以下步骤： A、 读取核磁共振接收机野外采集的全波数据s(k)及其采样率fs; B、 使用快速傅里叶变换求取全波数据的频谱S(U); C、 读取频谱S(?)中FID信号的频率fO所对应的相位0O; D、 读取频谱S(?)中，拉莫尔频率fO左右两侧的两个工频谐波频率fl和f2，及其所对 应的幅度al，a2,相位0 1和0 2 ; E、 重构两个工频谐波噪声的波形nl(k)和n2 (k); F、 在全波数据s(k)中减去nl(k)和n2 (k)，得到X(k); G、 将x(k)与参考信号C〇S(2 3ikfO/fS)进行互相关运算，得到R(m); H、 利用希尔伯特变换获取R(m)的包络信号，En(m); I、 将包络信号En(m)分为长度相等的两段，并分别使用双指数函数拟合包络曲线； J、 使用双指数函数的拟合结果重构包络信号，Env(m) =a?exp(bmfO/ fs)+c?exp(dmfO/fs)，a、b、c、d为拟合参数； K、 构造不含噪声的互相关波形rf(m) =Env(m) ?cosQjimfO/fs+QO);使用参考信号 cos(2 3ikf0/fs)对rf(m)进行解卷积运算，得到的解卷积结果ff(k)即为滤除噪声的核磁 共振全波信号。
【专利摘要】本发明涉及一种基于互相关的核磁共振全波信号噪声滤除方法，是利用噪声与拉莫尔频率的正弦信号不相关，而FID幅度衰减正弦信号与拉莫尔频率的正弦信号具有相关性的特点，通过互相关运算滤除噪声，然后拟合互相关波形的包络，并重构不含噪声的互相关波形，最后利用解卷积算法提取核磁共振全波数据中的FID信号。该方法运算数据计算量小，可以同时压制工频及其谐波噪声、随机噪声和尖峰噪声，明显提高核磁共振全波数据信噪比，有利于扩展核磁共振找水仪的应用范围和探测深度，而且该发明不需要额外的参考消噪装置，节约了成本，一次采集的全波数据通过互相关即可达到理想的滤波效果，节省测量时间。
【IPC分类】G01V3/14
【公开号】CN104898172
【申请号】CN201510256496
【发明人】林君, 刘立超
【申请人】吉林大学
【公开日】2015年9月9日
【申请日】2015年5月19日
最新回复(0)