本發(fā)明涉及地面核磁共振地下水探測信號噪聲濾除和參數(shù)提取技術領域，具體是利用基于工頻諧波建模和自相關的地面核磁共振信號提取方法。

背景技術：

地面核磁共振(Magnetic Resonance Sounding,MRS)作為一種地球物理新技術，近年來被廣泛應用到地下水探測中。與傳統(tǒng)的地球物理間接找水方法相比，地面核磁共振方法具有高分辨率，高效率，信息豐富和對水文參數(shù)具有唯一解釋的特點。地面核磁共振方法的基本原理是通過探測地下水中氫質子共振躍遷產生的MRS信號，提取其特征參數(shù)，再通過反演解釋來獲取地下水的水質信息。從一些國內外研究和應用中可以看出，MRS方法不僅在尋找地下水和地下水資源評價等方面取得了良好的應用效果，而且在獲取水文地質參數(shù)和定量化解釋方面也取得進展。

MRS方法雖然比其他傳統(tǒng)方法具有明顯優(yōu)勢，但也存在一些不足。實際探測到的MRS信號是非常微弱的納伏級信號，信號中噪聲干擾非常大，包括尖峰噪聲、隨機噪聲和工頻諧波噪聲等，其中由電力線干擾引起的工頻諧波噪聲對其影響最大。目前，針對濾除核磁共振信號中噪聲的理論和方法有很多。訾彥勇在著作《基于EMD的磁共振探測信號噪聲抑制方法》，提出一種新型的信號時頻處理方法，可以在沒有輸入信號任何先驗知識的情況下，自適應地將信號分解成若干個固有模態(tài)函數(shù)(Intrinsic Mode Functions,IMF)，實現(xiàn)信號趨勢的有效提取，但存在端點效應和模態(tài)混疊問題；田寶鳳等人在著作《基于參考線圈和變步長自適應的磁共振信號噪聲壓制方法》中直接利用觀測數(shù)據不斷遞歸更新處理參數(shù)，實現(xiàn)不同信噪比和信號強度下MRS信號噪聲的去除，但該方法需要遠端線圈同時采集噪聲數(shù)據做參考，對于單通道儀器無法實現(xiàn)，而且只能對相關噪聲處理，不能保證將多個工頻諧波噪聲全部消除。

CN104459809A公開了“基于獨立成分分析的全波核磁共振信號噪聲濾除方法”，該算法能夠有效地實現(xiàn)全波MRS信號的信噪分離，且數(shù)據擬合后初始振幅和弛豫時間的相對誤差小于±5.00％，相比于其他經典算法該算法更具有消噪性能的優(yōu)越性，但該方法只針對去除磁共振信號中的工頻諧波噪聲。CN104898172A公開了“一種基于互相關的核磁共振全波信號噪聲濾除方法”，該發(fā)明利用噪聲與拉莫爾頻率的正弦信號不相關，而MRS幅度衰減正弦信號與拉莫爾頻率的正弦信號具有相關性的特點，通過互相關運算濾除噪聲，然后擬合互相關波形的包絡，并重構不含噪聲的互相關波形，最后利用解卷積算法提取核磁共振全波數(shù)據中的MRS信號。該方法運算數(shù)據計算量小，可以同時壓制工頻及其諧波噪聲、隨機噪聲和尖峰噪聲，明顯提高核磁共振全波數(shù)據信噪比。但是在地磁場不穩(wěn)定時或信噪比較低時，拉莫爾頻率不能準確獲得，利用存在頻率偏差的信號進行互相關計算后，信號參數(shù)提取誤差大。CN104636609A中公開了“一種基于經驗模態(tài)分解與小波分析的信號聯(lián)合去噪方法”，根據信號的自相關性，對混有高斯白噪聲的信號進行EMD分解，由于EMD分解的性質，高斯白噪聲已不再是真正的白噪聲，但白噪聲的統(tǒng)計特性近似存在，即所述混有高斯白噪聲的信號的自相關函數(shù)在零點取得最大值，幅值隨著時間差的變化而變化，但其隨著時間的衰減很快。利用這種差異可以選取出噪聲起主導作用的IMF分量有效降低噪聲對信號的影響。但該方法主要是根據自相關特性對EMD分解后每個模態(tài)的屬性進行判斷，主要用于低信噪比下去除高斯白噪聲，所以其應用存在局限性。

技術實現(xiàn)要素：

本發(fā)明所要解決的技術問題在于提供一種地面核磁共振信號提取方法，解決了磁共振測深找水工作中由于強工頻諧波干擾和隨機噪聲造成的MRS信號有效提取的問題。

本發(fā)明是這樣實現(xiàn)的，一種地面核磁共振信號提取方法，該方法包括以下步驟：

步驟(1):利用地面核磁共振地下水探測儀器采集到一組觀測MRS含噪數(shù)據；

步驟(2):利用統(tǒng)計方法判斷是否存在尖峰噪聲，如果存在，去除尖峰噪聲并用插值結果代替，如果不存在，則測量數(shù)據保持不變；

步驟(3):將去除尖峰噪聲的數(shù)據，利用諧波建模的方法去除工頻諧波噪聲；

步驟(4):對步驟(3)的結果進行自相關和疊加處理減小隨機噪聲；

步驟(5):對步驟(4)后的結果進行MRS信號參數(shù)提取。

進一步地，步驟(3)中的建模方法的具體步驟為：

步驟3a：設定工頻諧波基頻的搜索范圍以及步長進行粗掃；

步驟3b：將工頻諧波模型轉化為Ax＝b的形式，解包含測量數(shù)據的線性方程組，得到工頻各諧波頻點的系數(shù)矩陣；

步驟3c：計算工頻諧波基頻搜索范圍內不同掃描值的模型估計值；

步驟3d：計算模型估計值與測量數(shù)據差值的2范數(shù)；

步驟3e：確定使誤差2范數(shù)取得最小值的頻點，從而獲得兩側相鄰的頻點；

步驟3f：將步驟3e中得到的兩側相鄰頻點作為下一次掃描范圍，選取適當步長再次掃描；

步驟3g：判斷是否掃描設定的次數(shù)，如果達到則根據誤差2范數(shù)最小的頻點進行工頻諧波建模；否則重復步驟3c～3g；

步驟3h:從步驟(2)的結果中減去工頻諧波模型，去除工頻諧波噪聲。

進一步地，步驟3f：步長為(f_m2-f_m1)/M，其中f_m1，f_m2分別為最小值頻點的兩側相鄰頻點，其中M＝3～7。

進一步地，步驟3g掃描設定的掃描次數(shù)為3-7次。

進一步地，所述的步驟(4)中對步驟(3)得到的結果進行自相關和疊加處理減小隨機噪聲的具體步驟為：

步驟4a：將MRS信號代入自相關函數(shù)公式，修改積分上下限，求得MRS信號的自相關表達式；

步驟4b：判斷采集時間是否大于1s，如果大于1s截取自相關結果的前半段用近似，其中e₀為初始振幅、為弛豫時間、f為拉莫爾頻率和τ為自相關時間間隔；如果小于1s自相0關結果用近似，其中t_max為時間最大值；

步驟4c：基于步驟4b的判斷后將多次測量數(shù)據的自相關結果進行疊加，進一步減小隨機噪聲。

進一步地，所述的步驟(5)中對步驟(4)得到的結果進行參數(shù)提取的具體步驟為：

步驟5a:對步驟(4)中疊加的數(shù)據進行希爾伯特變換，再通過低通濾波，轉化為兩個正交分量；

步驟5b:利用非線性擬合方法求解MRS信號的特征參數(shù)A₀和弛豫時間其中再計算初始振幅e₀結果，參數(shù)提取結果為初始振幅e₀和弛豫時間

本發(fā)明與現(xiàn)有技術相比，有益效果在于：本發(fā)明提出了基于工頻諧波建模和自相關的地面核磁共振信號提取方法，針對單通道采集的全波MRS數(shù)據，不僅可以一次性去除所有諧波干擾，而且還結合了自相關方法壓制隨機噪聲，最后通過非線性擬合實現(xiàn)了有效MRS信號特征參數(shù)的提取。本發(fā)明方法解決了磁共振測深找水工作中由于強工頻諧波干擾和隨機噪聲造成的MRS信號有效提取的難題，獲取的MRS信號關鍵特征參數(shù)的擬合誤差較小，同時本發(fā)明突破了經典消噪方法需多通道探測等其他條件的限制，節(jié)省了大量的人力物力。

附圖說明

圖1為本發(fā)明地面核磁共振信號提取方法流程框圖；

圖2為工頻諧波建模消噪原理及效果圖，圖2(a)是其時域信號，MRS信號淹沒在噪聲之中；圖2(b)是其功率譜，在諧波的頻點上均存在不同幅度的諧波噪聲；圖2(c)諧波建模后，在測量數(shù)據b中減去V_harmonic完成消噪過程，2(d)諧波建模消噪后，MRS信號的輪廓呈現(xiàn)衰減趨勢，在功率譜中只存在MRS信號的頻點；

圖3為自相關消噪效果圖；圖3(a)中黑色線是理想的MRS信號，參數(shù)為e₀＝100nV，f＝2326Hz和灰色線是加入了100nV均方差的高斯白噪聲后的測量信號，此時數(shù)據的信噪比為0dB，MRS信號基本淹沒在噪聲中；圖3(b)是測量信號和MRS信號的頻譜，可見MRS信號的頻譜峰值是噪聲頻譜均值的10倍。3(c)經過自相關函數(shù)后，測量信號(灰色線)和MRS信號(黑色線)的結果；3(d)MRS信號自相關后的頻譜峰值與噪聲頻譜；

圖4為自相關后非線性擬合效果圖，圖4(a)是未經過自相關函數(shù)的曲線，圖4(b)是經過自相關函數(shù)的曲線；

圖5為仿真數(shù)據和數(shù)據處理結果圖，圖5(a)含噪MRS信號；圖5(b)從含噪MRS信號中減去工頻諧波模型，得到去工頻信號圖；5(c)自相關和疊加處理得到的消噪后信號圖，5(d)MRS信號的實部，5(e)MRS信號的虛部，5(f)MRS信號實部的疊加曲線和虛部見；5(g)MRS信號虛部的疊加曲線；

圖6為實測數(shù)據處理結果圖；圖6(a-1)為q＝1.4As的MRS信號；6(a-2)為q＝4.2As的MRS信號；6(a-3)為q＝8.6As的MRS信號；

圖6(b-1)為q＝1.4As的信噪比信號；6(b-2)為q＝4.2As的的信噪比信號；6(b-3)為q＝8.6As的的信噪比信號。

具體實施方式

為了使本發(fā)明的目的、技術方案及優(yōu)點更加清楚明白，以下結合實施例，對本發(fā)明進行進一步詳細說明。應當理解，此處所描述的具體實施例僅僅用以解釋本發(fā)明，并不用于限定本發(fā)明。

如圖1所示，一種基于工頻諧波建模和自相關的地面核磁共振信號提取方法，包括以下步驟：

步驟(1):利用地面核磁共振地下水探測儀器采集到一組觀測MRS含噪數(shù)據；

步驟(2):利用統(tǒng)計方法(如羅曼諾夫斯基準則和能量運算等算法)判斷是否存在尖峰噪聲，如果存在，去除尖峰噪聲并用插值結果代替，如果不存在，則測量數(shù)據保持不變；

步驟(3):利用諧波建模的方法去除工頻諧波噪聲；

步驟(4):對步驟(3)的結果進行自相關和疊加處理減小隨機噪聲；

步驟(5):對步驟(4)后的結果進行MRS信號參數(shù)提??；

如圖2所示，步驟(3)中的工頻建模方法的具體步驟：

步驟3a:設定工頻諧波基頻f₀的搜索范圍49.9Hz～50.1Hz，f₀以步長為0.03Hz進行粗掃；

步驟3b:工頻諧波噪聲V_harmonic可表示為式(1)的形式

其中，A_n和分別是第n個諧波的幅度和相位，并有和t是時間。將表達式(1)整理成線性矩陣形式Ax＝b，

其中，b＝[V₁,V₂,…,V_P]^T中的V_p(p＝1,2,…P)是t_p時刻的接收數(shù)據；x＝[α₁,…,α_N,β₁,…,β_N]^T是諧波系數(shù)，N是諧波個數(shù)，根據接收信號的帶寬(1～3kHz)取N＝100。求解線性方程組(2)得到諧波系數(shù)矩陣x。

步驟3c:分別計算頻率搜索范圍內不同f₀的模型估計值V_harmonic(f₀)；

步驟3d:計算V_harmonic(f₀)與測量數(shù)據b的差值(誤差)的2范數(shù)：||b-V_harmonic(f₀)||₂；

步驟3e:確定使誤差2范數(shù)取得最小值時對應的頻點f_min，從而獲得兩側相鄰的頻點(f_m1,f_m2)；

步驟3f:確定下一次掃描范圍(f_m1,f_m2)，f₀選取步長為(f_m2-f_m1)/M，其中M＝3～7；

步驟g:判斷是否掃描5次(可以設置為3-7次，每次搜索的工頻基頻頻率相對變化量小于0.1％即可)，如果達到，則根據誤差2范數(shù)最小時對應的頻點f_min進行工頻諧波建模；否則重復步驟3c～3g；

步驟h:從步驟(2)的結果中減去工頻諧波模型V_harmonic(f_min)，從而去除工頻諧波噪聲；

如圖2所示的一組接收信號數(shù)據，圖2(a)是其時域信號，MRS信號淹沒在噪聲之中；圖2(b)是其功率譜，在諧波的頻點上均存在不同幅度的諧波噪聲。諧波建模后，在測量數(shù)據b中減去V_harmonic即完成了消噪過程，如圖2(c)和2(d)所示。諧波建模消噪后，MRS信號的輪廓呈現(xiàn)衰減趨勢，在功率譜中只存在MRS信號的頻點。

本發(fā)明提供的自相關和疊加處理的具體步驟為：

步驟4a：MRS信號用V_MRS表示：

其中，e₀表示初始振幅；表示平均弛豫時間，f表示拉莫爾頻率(與當?shù)氐卮艌鲇嘘P)，表示初始相位。

根據自相關原理，將V_MRS代入自相關函數(shù)的表達式中：

其中，τ為時間間隔。由于自相關函數(shù)是關于自變量τ的實偶函數(shù)，所以τ只取其正半軸。表達式(4)中的積分上下限分別取τ和t_max，t_max是測量時間長度，通過推導得到V_MRS的自相關表達式為：

R(τ)＝RC1(τ)+RC2(τ)+RC3(τ)+RC4(τ)

其中，

由表達式(5)～(7)可知，MRS信號自相關函數(shù)包含四個分量RC1～RC4，MRS信號自相關后的前半部分可以只用RC1代替；后半部分可以用RC1+RC2代替；RC3+RC4的結果很小，可以忽略。

步驟4b：判斷采集時間是否大于1s，如果大于1s截取自相關結果的前半段用RC1近似；如果小于1s自相關結果用RC1+RC2近似。

步驟4c：基于步驟4b的判斷后將多次測量數(shù)據b的自相關結果R進行疊加，即其中N是測量次數(shù)。N次疊加后信噪比提高倍，進一步減小隨機噪聲。

圖3給出了含噪MRS信號經自相關和疊加處理的結果。圖3(a)中黑色線是理想的MRS信號，參數(shù)為e₀＝100nV，f＝2326Hz和灰色線是加入了100nV均方差的高斯白噪聲后的測量信號，此時數(shù)據的信噪比為0dB，MRS信號基本淹沒在噪聲中。圖3(b)是測量信號和MRS信號的頻譜，可見MRS信號的頻譜峰值是噪聲頻譜均值的10倍。經過自相關函數(shù)后，測量信號(灰色線)和MRS信號(黑色線)的結果見圖3(c)所示。此時，灰色線與黑色線類似，也呈現(xiàn)明顯的衰減趨勢，只出現(xiàn)了少量的隨機噪聲，信噪比為23.75dB，與自相關之前相比大幅提高。從圖3(d)中的頻譜分析也可以得出，MRS信號自相關后的頻譜峰值約為噪聲頻譜均值的100倍，與自相關之前提高了一個數(shù)量級。因此，自相關函數(shù)能夠有效壓制隨機噪聲，提高信噪比。

如圖4所示，本發(fā)明提供的參數(shù)提取方法的具體步驟為：

步驟5a:對疊加的數(shù)據進行希爾伯特變換，希爾伯特變換表達式為：可將實信號R(t)轉化為復信號其包絡信號為再通過低通濾波，得到兩個正交分量V_RC＝V_real+i·V_imag；

步驟5b:利用非線性擬合方法對步驟a得到的兩個正交分量進行參數(shù)提取如式(8)，

獲得MRS信號的特征參數(shù)A₀和其中再計算e₀結果，參數(shù)提取結果為初始振幅e₀和弛豫時間

圖3中的數(shù)據經過Hilbert變換后的兩個分量見圖4所示，其中圖4(a)是未經過自相關函數(shù)的曲線，圖4(b)是經過自相關函數(shù)的曲線。分別利用非線性擬合方法求解得到的兩個分量的擬合曲線見圖4中黑色實線和黑色虛線所示。未經過自相關函數(shù)的數(shù)據信噪比低，參數(shù)提取結果分別為e₀＝108.3nV，f＝2326.1Hz；而經過自相關函數(shù)的數(shù)據信噪比提高明顯，參數(shù)提取結果為e₀＝100.5nV，f＝2326.0Hz，誤差分別降低了7.8％、10.0％和0.01％。因此，自相關函數(shù)能夠提高數(shù)據的信噪比，從而提高了參數(shù)提取的準確度。

實施例1

基于工頻諧波建模和自相關的地面核磁共振信號提取方法，本發(fā)明模擬了一組包含16次獨立采集的測量數(shù)據集，參照圖1，包括以下步驟：

步驟(1)：利用式(3)構造含噪MRS信號如圖5(a)所示，其參數(shù)為e₀＝100nV，f＝2326Hz和每次采集的諧波噪聲基頻在49.9～50.1Hz隨機產生，諧波個數(shù)為100，諧波幅度在200nV內隨機分布；每次采集的隨機噪聲為200nV的高斯白噪聲；采集數(shù)據中未加入尖峰噪聲，假設所有的尖峰噪聲已通過已有的方法消除。

步驟(2)：利用諧波建模的方法得到工頻諧波模型其中f₀＝50.018125Hz，從含噪MRS信號中減去工頻諧波模型，得到去工頻信號，如圖5(b)所示，可見工頻諧波消噪后，數(shù)據的信噪比大幅提升，平均提高了17.03dB，數(shù)據幅度減小到±500nV范圍內；

步驟(3):對步驟(2)的結果進行自相關和疊加處理，減小隨機噪聲，得到消噪后信號，如圖5(c)所示，可見自相關處理后，數(shù)據呈現(xiàn)明顯的衰減趨勢，除了信號初期存在脈沖(隨機噪聲的自相關特性)，其他時刻的數(shù)據幅度比自相關之前的數(shù)據明顯提高而隨機噪聲的影響較小，因此信噪比有較大程度的提升，為15.64～16.43dB，信噪比平均提高了16.10dB；

步驟(4):對步驟(3)的結果進行MRS信號參數(shù)提取，利用希爾伯特變換(參考頻率2325Hz)和低通濾波(截至頻率200Hz)，得到MRS信號的實部和虛部見圖5(d)和5(e)所示。建模消噪處理后，由于還存在大量的隨機噪聲，數(shù)據的實部和虛部曲線區(qū)分不明顯，只能觀察到曲線的整體趨勢稍有起伏(1Hz的震蕩)。相比自相關處理后，數(shù)據的實部和虛部能夠清晰地區(qū)分開，兩者都呈現(xiàn)出1Hz的震蕩衰減過程，再對16次采集和處理的數(shù)據進行疊加，得到MRS信號實部和虛部的疊加曲線見圖5(f)和(g)所示。疊加處理能夠進一步降低隨機噪聲，16次疊加信噪比提高4倍(12.04dB)，利用非線性擬合方法對兩組數(shù)據進行MRS信號參數(shù)提取，通過建模消噪和自相關處理后的參數(shù)提取結果為：e₀＝100.45nV，f＝2326.0Hz，相對誤差分別為：0.45％，-2.14％和0，均控制在±3％以內，滿足應用要求。

實施例2

本實施例以長春市農安縣燒鍋鎮(zhèn)的太平池水庫旁(該地地磁場為54720nT，對應的Larmor頻率在2330Hz，變化范圍小于2Hz)利用自主研制的地面核磁共振地下水探測儀器進行了野外數(shù)據采集與數(shù)據處理實驗。為了從淺到深地對含水層進行探測，在0.2～8.5As的范圍內由小到大按對數(shù)分布設置了20組發(fā)射脈沖矩，每組脈沖矩重復疊加16次，采樣率為50kHz，采集時間為1s。將第10、16和20個脈沖矩采集的MRS信號作為本發(fā)明方法的處理對象。如圖1所示，基于工頻諧波建模和自相關的地面核磁共振信號提取方法，包括以下步驟:

步驟(1):利用地面核磁共振地下水探測儀器采集到一組觀測MRS含噪數(shù)據；

步驟(2):利用統(tǒng)計方法判斷是否存在尖峰噪聲，如果存在，去除尖峰噪聲并用插值結果代替，如果不存在，則測量數(shù)據保持不變；

步驟(3)：利用諧波建模的方法得到模型其中f₀＝50.018125Hz，從步驟(2)的結果中減去工頻諧波模型，得到去工頻信號，如圖6(a-1)為q＝1.4As的MRS信號；6(a-2)為q＝4.2As的MRS信號；6(a-3)為q＝8.6As的MRS信號所示，可見實測數(shù)據中存在大量的諧波噪聲，噪聲水平在2500～3000nV范圍內，而MRS信號幅度在80～375nV，因此信噪比較低。經過工頻諧波消噪處理后，大多數(shù)諧波噪聲顯著降低，諧波基頻在49.975～50.041Hz之間，MRS信號(2330Hz)的峰值比較明顯，其中q＝1.4As的信號較大，而q＝4.2As和8.6As的信號較小，信噪比提高了近20dB；

步驟(4):對步驟(3)的結果進行自相關和疊加處理，減小隨機噪聲，得到消噪后信號；

步驟(5):對步驟(4)后的結果進行MRS信號參數(shù)提取，如圖6(b-1)，6(b-2)，6(b-3))所示，3個脈沖矩下信噪比增量分別為12.05dB、2.98dB和7.71dB，此時經過非線性擬合后，參數(shù)提取結果為：

q＝4.2As時，e₀擬合結果為80.7nV，不準確度為15.9nV，而擬合結果為235.9ms，不準確度為38.2ms；q＝1.4As時，e₀擬合結果為396.8nV，不準確度為17.4nV，而擬合結果為248.6ms，不準確度為10.8ms；q＝8.6As時，e₀擬合結果為175.6nV，不準確度為13.4nV，而擬合結果為263.1ms，不準確度為33.1ms，說明SNR的增加，能夠提高非線性擬合結果的準確度。

以上所述僅為本發(fā)明的較佳實施例而已，并不用以限制本發(fā)明，凡在本發(fā)明的精神和原則之內所作的任何修改、等同替換和改進等，均應包含在本發(fā)明的保護范圍之內。