一種基于互相關的核磁共振全波信號噪聲濾除方法
【技術領域】
[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] 信號經希爾伯特變換后，在頻域各頻率分量的幅度保持