基于頻域?qū)ΨQ法的磁共振數(shù)據(jù)殘余噪聲消除方法

林婷婷，李玥，劉大震，萬玲

（1.吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院，吉林長春，130026；2.地球信息探測儀器教育部重點實驗室，吉林長春，130026）

磁共振測深法（magnetic resonance sounding，MRS）自20世紀80年代末期開始逐漸發(fā)展，成為最具競爭力的地球物理勘探方法之一，它具有非侵入式直接追蹤含水體分布信息的優(yōu)勢，被廣泛應(yīng)用于地下水資源探測和災(zāi)害水源勘查等領(lǐng)域[1-3]。然而，MRS 信號非常微弱（nV 級），探測過程易受環(huán)境噪聲干擾，導(dǎo)致反演獲取的地下水信息失真。目前的數(shù)據(jù)采集技術(shù)無法完全屏蔽環(huán)境噪聲，限制了MRS 方法在復(fù)雜環(huán)境中的探測潛力[4-6]。因此，如何最大限度地消除環(huán)境噪聲對MRS 信號的干擾是一個亟待解決的難題。

常見的MRS 環(huán)境噪聲主要分為3 類：尖峰噪聲、工頻諧波噪聲和隨機噪聲[7-9]。

1）尖峰噪聲通常由物體忽然放電產(chǎn)生，例如太陽磁暴或靜電放電。其特點是持續(xù)時間短，干擾幅度大。在數(shù)據(jù)處理過程中，大幅度尖峰噪聲所在的數(shù)據(jù)段會被丟棄；對于不易識別的小幅度尖峰噪聲，萬玲等[10]提出了基于能量運算的尖峰噪聲抑制方法進行消除。LARSEN[11]以2 個級聯(lián)二階濾波器的脈沖響應(yīng)建模，專用于消除來自電網(wǎng)的尖峰噪聲。

2）工頻諧波噪聲較恒定，其頻率固定在電力線基頻的整數(shù)倍處，通常來源于電力線、變壓器等工業(yè)及生活用電。WALSH[12]開發(fā)了多通道MRS地下水探測系統(tǒng)，基于通道間數(shù)據(jù)相關(guān)性提出了頻域自適應(yīng)消噪算法。這種方法可以自動調(diào)節(jié)濾波器系數(shù)跟蹤工頻諧波噪聲峰值位置，但其作用頻率范圍有限，且對參考源信號成分的單一性要求較高。之后，MüLLER-PETKE 等[13]研究發(fā)現(xiàn)，設(shè)置最佳濾波器寬度并延長數(shù)據(jù)記錄時間能夠顯著提高參考消噪方法的有效性。LARSEN等[14]提出了建模法與維納濾波法結(jié)合來消除工頻噪聲。HEIN等[15]提出頻域?qū)ΨQ法，實現(xiàn)了MRS信號中工頻噪聲的進一步消除。

3）隨機噪聲由外界環(huán)境和儀器自身的不可預(yù)知因素產(chǎn)生，在統(tǒng)計特征上具有隨機性。在實際應(yīng)用中常采用疊加法來抑制隨機噪聲，降低了探測效率[5,16]。為此，LIN等[17-19]提出了時頻峰值濾波法，可以有效抑制單次測量的隨機噪聲。BEHROOZMAND 等[1]總結(jié)了傳統(tǒng)MRS 消噪流程，其中包含4個步驟：1）去尖峰噪聲；2）消除工頻諧波噪聲；3）壓制隨機噪聲；4）包絡(luò)探測提取信號[20-22]。

采用傳統(tǒng)MRS 消噪流程處理數(shù)據(jù)可以在一定程度上削弱各類噪聲，提高信噪比。但受消噪方法原理限制，處理后的信號中仍然殘余小幅度工頻諧波噪聲和高斯隨機噪聲，影響信號特征參數(shù)提取準確度。針對這一問題，本文作者在頻域?qū)ΨQ法消除工頻噪聲的基礎(chǔ)上，分析該方法對殘余小幅度工頻諧波噪聲和高斯隨機噪聲的抑制效果，在數(shù)據(jù)處理過程中將該方法用于壓制隨機噪聲（第3 步），以MRS 包絡(luò)信號在頻域存在對稱性為理論基礎(chǔ)，對含噪MRS信號進行頻域解調(diào)與時頻變換，區(qū)分含噪信號的實虛分量；然后，利用工頻諧波噪聲的虛部分量估計實部分量，剔除后實現(xiàn)信噪分離。

1 基于頻域?qū)ΨQ法的MRS 探測信號消噪方法

1.1 磁共振探測地下水原理

地下水中氫質(zhì)子以地磁場強度B0決定的拉莫爾頻率fL繞背景場進動。磁共振方法在發(fā)射線圈中通以頻率為fL的交變電流，產(chǎn)生與地磁場方向垂直的交變磁場。當停止激發(fā)電流后，氫質(zhì)子的磁化強度不僅要圍繞激發(fā)電流產(chǎn)生的交變磁場進動，還要逐漸回落到地磁場B0的方向上，并產(chǎn)生MRS信號，通過地面鋪設(shè)的接收線圈獲取。MRS 信號經(jīng)過反演后，即可確定地下水體含水量、含水層結(jié)構(gòu)、含水層深度等信息[2,23]。MRS探測方法原理圖如圖1所示，感應(yīng)信號E的表達式為

圖1 MRS探測方法原理圖Fig.1 Schematic diagram of MRS detection method

其中：E0為初始振幅，E0反演后可得到含水量；T2*為氫原子的弛豫時間，表征含水層介質(zhì)孔隙度；φ0為初始相位，由地下電導(dǎo)率產(chǎn)生。噪聲會影響信號特征參數(shù)提取準確度，導(dǎo)致對地下水物理信息的錯誤判斷[24-26]。

1.2 MRS信號與工頻諧波噪聲頻域分析

在處理MRS 數(shù)據(jù)時，通?？珊雎猿跏枷辔坏挠绊懀处?=0）[1,26]，則磁共振信號E由式（2）和（3）所示的指數(shù)函數(shù)和余弦函數(shù)的乘積構(gòu)成。

根據(jù)傅里葉變換的性質(zhì)，式（1）中的傅里葉變換式等價于式（2）和（3）中的傅里葉變換形式（4）和（5）的卷積[27]。

式（4）所示為MRS包絡(luò)信號頻譜，在復(fù)平面上實部為偶函數(shù)，虛部為奇函數(shù)，經(jīng)傅里葉逆變換到時域后為純實數(shù)信號。式（5）中ω0=2πfL，是在f=±fL處高度為的沖激函數(shù)。式（4）和（5）卷積是對式（2）的頻譜調(diào)制，相當于將式（2）的頻譜進行尺度變換后由f=0搬移到f=±fL處。提取信號特征參數(shù)需要對式（1）進行頻譜解調(diào)（頻譜調(diào)制的逆過程）獲取MRS包絡(luò)信號。由于包絡(luò)信號頻譜隨|ω|增加不斷趨近于0，因此，可選擇有限帶寬內(nèi)的頻譜近似解調(diào)。

首先，選擇以拉莫爾頻率為中心的有限帶寬截取頻譜；然后，將所截頻譜沿頻率軸左移fL個單位；最后，將頻譜幅值放大2 倍，即可得到MRS包絡(luò)信號頻譜。設(shè)置MRS 信號的初始振幅E0=50 nV，弛豫時間T*2=100 ms，拉莫爾頻率fL=200 Hz。時域信號如圖2（a）中藍色曲線所示，解調(diào)前MRS信號頻譜峰值位于頻率f=±200 Hz處，如圖2（b）中虛線所示。解調(diào)后頻譜如圖2（b）中實線所示，可以看出該頻譜具有式（4）所述對稱性，經(jīng)傅里葉逆變換得到時域MRS 包絡(luò)信號為純實數(shù)，如圖2（a）中黑色實線所示，包含初始振幅E0和弛豫時間T2*的信息。

圖2 仿真MRS信號和解調(diào)后MRS信號的時域圖和頻域圖Fig.2 Time-domain and frequency-domain plots of a simulated MRS signal before and after demodulation

工頻諧波噪聲會破壞MRS 包絡(luò)信號的頻域?qū)ΨQ性，此時，傅里葉逆變換得到的時域數(shù)據(jù)為復(fù)數(shù)，實部由MRS 包絡(luò)信號和工頻諧波噪聲解調(diào)后的實時域分量共同構(gòu)成，而虛部完全由工頻諧波噪聲產(chǎn)生。在MRS 信號中加入1 個工頻噪聲成分（見圖3），幅度EN=10 nV，頻率fN=250 Hz。圖3（a）中藍色曲線為加入工頻諧波噪聲的MRS信號；黑色實線為解調(diào)后實部分量，具有衰減趨勢，是含噪信號的包絡(luò)曲線；紅色實線為工頻諧波噪聲解調(diào)后虛部分量。

對于頻譜中f=±|fL+x|處高度為λ+μi 的工頻諧波噪聲，解調(diào)之后噪聲峰值位于f=x處，f=-x處噪聲為0，如圖3（b）所示。由數(shù)學(xué)定理可知：任意函數(shù)都可以表示為1個奇函數(shù)和1個偶函數(shù)的和，且表達形式唯一[27]。則解調(diào)后的噪聲可以分解為Cre，Cro，Cie和Cio共4 部分，其表達 式分別為：

圖3 含1個諧波成分的MRS信號解調(diào)前后的時域圖和頻域圖Fig.3 Time-domain and frequency-domain plots of an MRS signal with a harmonic component before and after demodulation

Cre和Cro以及Cie和Cio之間的區(qū)別是在f=-x處的峰值符號相反，當上述4 個分量相加時，f=-x處的噪聲峰值會相互抵消，f=x處的峰值疊加產(chǎn)生工頻諧波噪聲λ+μi。Cre和Cio經(jīng)傅里葉逆變換到時域后是一個純實數(shù)序列，這就是與MRS 信號混疊在一起的工頻諧噪聲實時域分量，如圖4（a）中黑色曲線所示。Cro和Cie經(jīng)傅里葉逆變換到時域后得到純虛數(shù)序列，也就是工頻諧波噪聲解調(diào)后的虛時域分量，如圖4（a）中紅色曲線所示。

圖4 諧波噪聲解調(diào)后時域圖及其頻譜Fig.4 Time-domain and frequency-domain plots of a harmonic component after demodulation

1.3 頻域?qū)ΨQ法消噪原理

頻域?qū)ΨQ消噪理論是基于MRS 信號的頻域?qū)ΨQ性展開研究的。首先對含噪聲的MRS 數(shù)據(jù)進行傅里葉變換、頻域解調(diào)、傅里葉逆變換，得到時域復(fù)數(shù)數(shù)據(jù)。該復(fù)數(shù)虛部由工頻諧波噪聲產(chǎn)生，傅里葉變換后可得到已知分量Cro和Cie。復(fù)數(shù)實部中無法直接提取工頻諧波噪聲實部信息，因此，對噪聲虛部分量Cro和Cie在f=-x處的峰值乘-1，估計出未知的噪聲實部分量Cre和Cio。然后，在頻譜中f=±x處減去估計噪聲峰值并進行頻譜校正，再通過頻譜調(diào)制和傅里葉逆變換得到消噪后的MRS信號。

頻域?qū)ΨQ法處理MRS 數(shù)據(jù)時得到的虛時域分量中包含部分高斯隨機噪聲，這是由于頻域解調(diào)是在頻譜上有限連續(xù)帶寬內(nèi)進行，且高斯隨機噪聲的無規(guī)律性決定其在變換之后產(chǎn)生的時域數(shù)據(jù)也是復(fù)數(shù)。虛時域分量在估算噪聲實部分量之后被丟棄時，其中包含的高斯隨機噪聲也同時被消減。因此，頻域?qū)ΨQ法不僅可以消除工頻諧波噪聲的干擾，而且可以消除部分高斯噪聲。在實際應(yīng)用中，將頻域?qū)ΨQ法作為現(xiàn)有消噪流程的補充步驟，用在疊加法消除隨機噪聲之后，以期消除MRS 信號中仍殘余的工頻諧波噪聲和高斯隨機噪聲，進一步提高信噪比。頻域?qū)ΨQ法消噪流程圖如圖5所示。

圖5 頻域?qū)ΨQ法消噪流程圖Fig.5 Flowchart of frequency-domain symmetry method

2 仿真結(jié)果及分析

為了驗證頻域?qū)ΨQ法消除殘余噪聲的效果，在仿真MRS 信號中加入不同幅度工頻諧波噪聲和高斯隨機噪聲，仿真信號參數(shù)E0=300 nV，T*2=150 ms，fL=2 330 Hz。采用Larsen 建模法和疊加法分別消除工頻諧波噪聲和隨機噪聲，然后用頻域?qū)ΨQ法消除殘余噪聲，通過對比頻域?qū)ΨQ法消噪前后信噪比、初始振幅提取誤差和弛豫時間提取誤差對算法進行量化評估。各項評價指標的定義如下。

信噪比σSNR定義為

式中：E和分別為仿真信號和信號估值。

初始振幅提取相對誤差σIAEE定義為

式中：為初始振幅的提取值。

弛豫時間提取相對誤差σRTEE定義為

式中：為弛豫時間的提取值。

2.1 不同幅度工頻諧波噪聲數(shù)據(jù)處理結(jié)果

為了驗證頻域?qū)ΨQ法對殘余工頻諧波噪聲的消除效果，在MRS 仿真信號中加入10 組30～300 nV 內(nèi)的工頻諧波噪聲AN，信噪比由23.38 dB下降到3.38 dB，每組噪聲對應(yīng)信噪比σSNR如表1所示。采用頻域?qū)ΨQ法處理各組數(shù)據(jù)后提取信號特征參數(shù)，所得結(jié)果均為=299.87 nV，=150.1ms，參數(shù)提取相對誤差σIAEE=0.04%，σRTEE=0.07%，信噪比提高到32.49 dB。實驗結(jié)果表明，即使工頻諧波噪聲殘余很小，頻域?qū)ΨQ法也能準確消除噪聲，并且在不同噪聲水平下信號提取結(jié)果相同，誤差很小，說明采用頻域?qū)ΨQ法消噪效果十分穩(wěn)定。

表1 不同幅度工頻諧波噪聲下的信噪比Table 1 SNR after adding powerline harmonics with different amplitudes

圖6所示為300 nV 工頻諧波噪聲消噪過程圖。在工頻諧波噪聲影響下，MRS 信號上重復(fù)出現(xiàn)高度約為300 nV 的峰值，信號衰減趨勢減弱。對含噪結(jié)果進行擬合，提取信號特征參數(shù)=247.38 nV，=259.0 ms，信噪比σSNR=3.38 dB，信號準確度嚴重受損。采用頻域?qū)ΨQ法處理數(shù)據(jù)。

首先，進行傅里葉變換、解調(diào)和逆傅里葉變換，由于MRS 包絡(luò)信號為純實數(shù)，因此，圖6（b）中實部曲線具有明顯的衰減趨勢，而圖6（b）中虛部曲線在水平直線附近上下震蕩。2條曲線中出現(xiàn)的重復(fù)性峰值說明時域數(shù)據(jù)的實部和虛部均含有工頻諧波噪聲。然后，對實時域序列和虛時域序列進行傅里葉變換，估算工頻諧波噪聲實時域分量，消減后獲取MRS 包絡(luò)信號頻譜如圖6（d）黑色曲線所示。將該頻譜逆傅里葉變換到時域后，包絡(luò)信號指數(shù)衰減特征明顯，原有重復(fù)性特征全部消除。最后，對消噪后的包絡(luò)信號頻譜調(diào)制，得到頻域特征如圖6（f）中黑色曲線所示，然后經(jīng)傅里葉逆變換到時域，圖6（e）中黑色曲線為消噪后MRS 全波信號。對處理后數(shù)據(jù)作非線性擬合，信號特征參數(shù)提取結(jié)果為：=299.87 nV，=150.1ms，準確度分別提高了438.5倍和1 038.1倍，信噪比提升至32.49 dB?？梢姡褂妙l域?qū)ΨQ法消除殘余工頻諧波噪聲可以提高提取信號準確度，其中對弛豫時間準確度的影響最顯著。

圖6 300 nV工頻諧波噪聲消噪過程圖Fig.6 Process of eliminating powerline harmonic noise with an amplitude of 300 nV

2.2 不同幅度隨機噪聲數(shù)據(jù)處理結(jié)果

為了驗證頻域?qū)ΨQ法對殘余隨機噪聲的消除效果，在MRS 仿真信號中加入13 組0～60 nV 內(nèi)的高斯隨機噪聲，對比消噪前后的特征參數(shù)提取結(jié)果和相對誤差，判定頻域?qū)ΨQ法消除殘余高斯隨機噪聲的有效性，如圖7所示。從圖7可以看出：采用頻域?qū)ΨQ法消噪之后，特征參數(shù)提取結(jié)果更加接近仿真值，信號提取誤差變小。當殘余高斯隨機噪聲很小（5 nV）時，采用頻域?qū)ΨQ法消噪后，初始振幅E0由298.41 nV 提高到299.37 nV，初始振幅提取相對誤差由0.53%降低到0.21%，弛豫時間T*2由152.4 ms降低到151.1 ms，弛豫時間提取相對誤差由1.61%降低到0.71%，信噪比由21.34 dB提高到24.77 dB。說明在高斯隨機噪聲殘余較小的情況下，頻域?qū)ΨQ法也能將其消除，改善MRS 信號的提取質(zhì)量。

圖7 不同幅度高斯隨機噪聲水平下提取信號結(jié)果Fig.7 Signal fitting results after adding random noise with different amplitudes

2.3 仿真數(shù)據(jù)處理結(jié)果

為了驗證頻域?qū)ΨQ法作為現(xiàn)有消噪流程補充步驟的必要性，仿真9 組16 次采集的MRS 數(shù)據(jù)，假設(shè)數(shù)據(jù)中尖峰噪聲已經(jīng)被消除，工頻諧波噪聲和隨機噪聲分別在3 000 nV 和100 nV 以內(nèi)漸增式隨機產(chǎn)生。采用建模法、疊加法和頻域?qū)ΨQ法依次消除工頻諧波噪聲、隨機噪聲和殘余噪聲，其中1 組消噪結(jié)果如圖8所示。圖8中藍色曲線是采用建模法對單次數(shù)據(jù)直接消除工頻諧波噪聲后的處理結(jié)果，綠色曲線是建模法和疊加法共同消除噪聲之后的結(jié)果，紅色曲線是進一步采用頻域?qū)ΨQ法消除殘余噪聲后的結(jié)果。頻域?qū)ΨQ法消除殘余噪聲后，全波信號衰減趨勢更強，包絡(luò)信號擬合結(jié)果更接近仿真信號包絡(luò)曲線，信噪比由15.29 dB提升至19.33 dB，初始振幅提取相對誤差由2.03%減小為0.86%，弛豫時間提取相對誤差由5.53%減小為2.27%。因此，在傳統(tǒng)消噪流程之后增加頻域?qū)ΨQ法可以消除殘余噪聲，進一步提高信噪比和特征參數(shù)提取準確度，改善提取MRS信號的質(zhì)量。圖9所示為MRS 信號的三維時頻分布圖，經(jīng)建模法、疊加法和頻域?qū)ΨQ法處理后的數(shù)據(jù)，噪聲極大程度被消除，MRS信號成分有效保留。

圖8 仿真數(shù)據(jù)處理結(jié)果時域圖Fig.8 Time-domain plot of simulated processing results

圖9 仿真數(shù)據(jù)處理結(jié)果三維時頻圖Fig.9 Time-frequency domain of simulated processing results

表2所示為使用不同方法處理仿真數(shù)據(jù)后的特征參數(shù)提取情況和信噪比。從表2可以看出：頻域?qū)ΨQ法參數(shù)提取結(jié)果優(yōu)于建模法參數(shù)提取結(jié)果，信噪比平均提高4.23 dB，但在建模法中加入疊加法消除隨機噪聲之后，提取信號準確度高于單獨使用頻域?qū)ΨQ法消噪時的信號準確度，信噪比平均提高7.81 dB，若在疊加法之后加入頻域?qū)ΨQ法，信號質(zhì)量會得到進一步改善，信噪比平均提高3.78 dB。因此，本文提出的頻域?qū)ΨQ法不宜單獨使用，而應(yīng)當作為傳統(tǒng)消噪流程的補充步驟，消除MRS 信號中殘余的工頻諧波噪聲和隨機噪聲，進一步提高信噪比。

表2 使用不同方法處理仿真數(shù)據(jù)后的特征參數(shù)提取情況和信噪比Table 2 Parameter fitting and signal-to-noise ratio after processing simulated data using different methods

3 實測數(shù)據(jù)處理

為了驗證本文提出算法的實用性，采用自主研制的MRS 地下水探測儀[1]在長春市燒鍋鎮(zhèn)采集數(shù)據(jù)。當?shù)乩獱栴l率fL=2 337 Hz。3 組發(fā)射脈沖矩Q在0.99～1.84 A·s 之間均勻分布。采樣頻率fs=25kHz，采集時間256 ms，每組脈沖矩疊加16次。

對實測數(shù)據(jù)進行尖峰噪聲處理后，得到3個脈沖矩信號的處理結(jié)果如圖10中灰色曲線所示，信號完全淹沒在工頻諧波噪聲和隨機噪聲中。采用建模法和疊加法分別消除工頻噪聲和隨機噪聲，數(shù)據(jù)處理結(jié)果如圖10中紅色曲線所示，可以看出實測數(shù)據(jù)中大部分噪聲被壓制，信號初見衰減趨勢，但數(shù)據(jù)中仍有部分噪聲殘留，Q為0.99，1.46 和1.84 A·s 時脈沖矩信號的信噪比σSNR分別為-3.95，-4.14和-4.31 dB。信號特征參數(shù)E0分別為382.13，383.30 和397.62 nV，T*2分別為259.9，380.0 和348.4 ms。進一步采用頻域?qū)ΨQ法消除殘余噪聲，結(jié)果如圖10中黑色曲線所示。時域信號衰減趨勢更明顯，頻域中紅色殘余噪聲被消除，信號峰譜更突出，計算Q為0.99，1.46 和1.84 A·s時脈沖矩信號的信噪比σSNR分別為-1.36，-0.61和-0.83 dB。信號特征參數(shù)E0分別為252.63，330.12和286.49 nV，T*2分別為169.9，168.5 和174.9 ms。相比于頻域?qū)ΨQ法消噪之前，3個脈沖矩信號的信噪比增量分別為2.59，3.53 和3.48 dB，初始振幅和弛豫時間提取準確度均有所提高。

圖10 實測數(shù)據(jù)處理結(jié)果Fig.10 Processing results of field data

4 結(jié)論

1）頻域?qū)ΨQ法不宜單獨使用，而應(yīng)當作為傳統(tǒng)消噪流程的補充步驟，進一步削弱環(huán)境噪聲對探測信號的影響，提高特征參數(shù)提取準確度。仿真結(jié)果表明，加入頻域?qū)ΨQ法后信號信噪比平均提高了3.78 dB。

2）即使殘余工頻諧波噪聲和高斯隨機噪聲幅度很小，頻域?qū)ΨQ法也能準確消除噪聲，提高獲取信號的準確性和穩(wěn)定性，為核磁共振測深法的廣泛應(yīng)用提供技術(shù)支持。實測數(shù)據(jù)處理結(jié)果表明，頻域?qū)ΨQ法消除殘余噪聲后，3個脈沖矩信號的信噪比增量分別為2.59，3.53和3.48 dB。

3）由于隨機噪聲成分復(fù)雜，對MRS 信號的影響無法完全消除，下一步工作將繼續(xù)研究隨機噪聲干擾嚴重情況下的MRS消噪方法。