基于能量運(yùn)算的磁共振信號(hào)尖峰噪聲抑制方法

萬(wàn)玲， 張揚(yáng)， 林君， 蔣川東， 林婷婷

吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院/地球信息探測(cè)儀器教育部重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)春　130026

?

基于能量運(yùn)算的磁共振信號(hào)尖峰噪聲抑制方法

萬(wàn)玲， 張揚(yáng)， 林君， 蔣川東， 林婷婷*

吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院/地球信息探測(cè)儀器教育部重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)春130026

摘要磁共振探測(cè)信號(hào)微弱，使用高靈敏度的核磁共振地下水探測(cè)儀，極易受到環(huán)境噪聲干擾.其中，工頻諧波噪聲和尖峰噪聲，是影響信號(hào)質(zhì)量最嚴(yán)重的兩類噪聲.國(guó)內(nèi)外研究表明，通過(guò)參考線圈的布設(shè)，依據(jù)探測(cè)線圈和參考線圈中噪聲相關(guān)性，利用自適應(yīng)參考對(duì)消算法，能夠有效抑制工頻諧波噪聲.然而，尖峰噪聲的存在嚴(yán)重影響了主通道與參考道的數(shù)據(jù)相關(guān)性，成為了參考對(duì)消算法應(yīng)用的難題與障礙.為解決這一問(wèn)題，本文提出磁共振信號(hào)中尖峰噪聲的抑制方法，推導(dǎo)了能量域磁共振信號(hào)表達(dá)式，通過(guò)計(jì)算信號(hào)能量，可有效檢測(cè)尖峰噪聲并突出不易識(shí)別的小幅度尖峰噪聲，采用基于中位數(shù)的絕對(duì)偏差法確定閾值，進(jìn)而剔除尖峰噪聲.為了驗(yàn)證去噪效果，與應(yīng)用較廣的統(tǒng)計(jì)疊加法進(jìn)行對(duì)比研究.仿真結(jié)果表明，對(duì)于干擾幅度較大、持續(xù)時(shí)間較長(zhǎng)的尖峰噪聲，能量運(yùn)算法和統(tǒng)計(jì)疊加法均能識(shí)別并剔除，且不損失有效的磁共振信號(hào)，標(biāo)準(zhǔn)差偏差控制在0.3%以內(nèi)，可以滿足實(shí)際應(yīng)用的要求；對(duì)于小于信號(hào)幅度1.5倍的尖峰噪聲，能量運(yùn)算法可有效識(shí)別并剔除，而統(tǒng)計(jì)疊加法無(wú)法實(shí)現(xiàn).針對(duì)多通道探測(cè)系統(tǒng)，使用能量運(yùn)算法剔除尖峰噪聲后，可明顯提高主通道和參考道的數(shù)據(jù)相關(guān)性，為后續(xù)自適應(yīng)參考對(duì)消算法的應(yīng)用奠定了基礎(chǔ).實(shí)測(cè)數(shù)據(jù)處理結(jié)果進(jìn)一步證明了本文方法的實(shí)用性.

關(guān)鍵詞磁共振信號(hào)； 參考對(duì)消算法； 能量運(yùn)算； 尖峰噪聲； 噪聲相關(guān)性

In the present paper, we provide a detailed insight into the technique of spike noise cancellation based on transferring the signal from time domain to energy domain and using the median absolute deviation method (MAD). First, we calculate the energy of MRS signal. After doing this, the signal is emphasized which is much more clear than in time domain. Then, a threshold principle based on MAD is provided to cancel all the noise above it.

In addition, this paper contains a comparison of the spike noise cancellation effects of Statistical Stacking and Energy cancellation. It is found that the two methods provide identical noise cancellation performance when the spike noise is strong and long enough. Statistical Stacking method is limited when the spike noise get weaker, especially weaker than 1.5 times of the signal. But Energy cancellation we proposed can still remove the spike noise and keep the MRS signal effectively. Moreover, from the noise properties analysis, correlation between the detection loop and the reference loop can be improved after spike noise cancellation. We first apply it on the synthetic data to see the correlation improvement and find that the correlation is improved from 0.1369 to 0.4941 after using the Energy cancellation.Then we obtain the true field data using the multi-channel MRS instrument in Changchun suburb, and apply the Energy cancellation method on the field data. The correlation is improved from 0.1575 to 0.2481. We also find that the harmonic noise is much more evident after spike noise cancellation, which is advantageous to use the adaptive noise cancellation later.

We anticipate that better noise cancelling results can be obtained with the Energy cancellation and adaptive noise cancellation methods used together and therefore a wider application of multi-channel MRS instrument can be made in future.

1引言

地下水中的氫原子被電磁場(chǎng)激發(fā)，會(huì)形成宏觀磁矩，在地面鋪設(shè)線圈可拾取宏觀磁矩進(jìn)動(dòng)產(chǎn)生的核磁共振信號(hào)(Legchenko and Valla, 2002; Legchenko et al.,2002；Lubczynski and Roy,2003；Roy and Lubczynski,2003；Schirov et al.,1991).信號(hào)的強(qiáng)弱與含水量大小有關(guān)，通過(guò)觀察拾取信號(hào)的強(qiáng)弱可判斷地下水含量.另外，不同強(qiáng)度的電磁場(chǎng)可激發(fā)地下不同深度的含水體.通過(guò)改變激發(fā)電磁場(chǎng)的強(qiáng)度，可獲得地下不同深度的含水量信息.與常規(guī)的物探找水方法相比，利用核磁共振原理探尋地下水，不打鉆就能確定地下含水層的深度和含水量大小等信息(林君等，2010；潘玉玲和張昌達(dá)，2000).因此，核磁共振測(cè)深(Magnetic Resonance Sounding, MRS)方法，作為一種直接探測(cè)地下水的地球物理方法(林君，2010；林君等，2012；陸其鵠等，2009；張榮等，2006)，已經(jīng)得到了世界上越來(lái)越多的關(guān)注.

但是，由于MRS信號(hào)非常微弱，高靈敏度的接收裝置極易受到周圍環(huán)境噪聲的干擾(林君等，2010)，MRS信號(hào)常常被淹沒(méi)在噪聲中，主要包括工頻諧波噪聲、尖峰噪聲等.隨著國(guó)際多通道磁共振儀器的誕生，工頻諧波噪聲可以通過(guò)參考通道的布設(shè)，利用其與主通道的相關(guān)性，通過(guò)自適應(yīng)參考對(duì)消算法消除(Fabian,2010；Müller-Petke and Yaramanci,2010,2011a，2011b；Radic，2006；Walsh，2008；田寶鳳等，2012).然而，強(qiáng)尖峰噪聲干擾的存在，將影響自適應(yīng)參考對(duì)消算法對(duì)工頻噪聲的剔除效果，進(jìn)而阻礙MRS信號(hào)的有效提取，降低反演結(jié)果對(duì)水文地質(zhì)參數(shù)解釋的準(zhǔn)確度.

尖峰噪聲主要是由太陽(yáng)磁暴、雷暴或者任何物體突然放電等引起的.這種尖峰噪聲的主要特征是，時(shí)域上的分布和幅度具有偶然性和隨機(jī)性，并且幅度大于或遠(yuǎn)大于信號(hào)幅度；頻域上頻率不固定，且頻譜分布范圍廣，通常其覆蓋從幾赫茲到100 MHz以上，與MRS信號(hào)頻率混疊在一起(林君等，2010).針對(duì)尖峰噪聲對(duì)MRS信號(hào)干擾問(wèn)題，國(guó)內(nèi)外專家和學(xué)者進(jìn)行了相應(yīng)的研究.王中興等(2009)提出差閾值替代疊加方法削弱和抑制奇異噪聲，提高信號(hào)參數(shù)的提取精度；蔣川東等(2012)采用統(tǒng)計(jì)疊加方法去除尖峰噪聲干擾；Strehl(2006), Strehl等(2006)根據(jù)尖峰噪聲對(duì)MRS信號(hào)的干擾特點(diǎn)，提出了基于小波硬閾值變換方法去除MRS信號(hào)中的尖峰噪聲.差閾值和統(tǒng)計(jì)疊加的方法能夠判別出幅值較大的尖峰噪聲，然而當(dāng)尖峰噪聲幅度較小時(shí)，難以判別噪聲從而不能獲得好的消噪效果；采用小波變換的方法進(jìn)行噪聲去除，其去噪效果受到分解級(jí)數(shù)和閾值策略的限制，此外，當(dāng)尖峰噪聲持續(xù)時(shí)間較長(zhǎng)時(shí)，容易損失有效磁共振信號(hào).綜上，磁共振尖峰噪聲處理中需要一種對(duì)其進(jìn)行有效識(shí)別與去除的技術(shù).

本文借鑒了尖峰噪聲剔除方法在醫(yī)學(xué)領(lǐng)域(Mukhopadhyay and Ray,1998)的成功應(yīng)用，提出基于能量運(yùn)算的磁共振信號(hào)尖峰噪聲抑制方法.通過(guò)計(jì)算信號(hào)能量，突出不易識(shí)別的小幅度尖峰噪聲，選取適當(dāng)?shù)拈撝颠M(jìn)行尖峰識(shí)別并剔除，進(jìn)而有效提高主通道和參考道的數(shù)據(jù)相關(guān)性，為自適應(yīng)參考對(duì)消算法的廣泛應(yīng)用奠定基礎(chǔ).

2核磁共振探測(cè)噪聲分析及對(duì)策

2.1核磁共振地下水探測(cè)

將發(fā)射線圈鋪在地表，供入頻率為拉莫爾頻率的脈沖電流，形成激發(fā)磁場(chǎng)，當(dāng)脈沖終止后，利用接收線圈接收磁共振全波信號(hào)E(t)，表達(dá)式為：

(1)

圖1　發(fā)射的電流脈沖與接收的磁共振信號(hào)示意圖Fig.1　Excited pulse and received MRS signal

磁共振信號(hào)微弱，通常為nV級(jí).為提高M(jìn)RS信號(hào)有效獲取能力，可應(yīng)用參考對(duì)消算法的多通道磁共振探測(cè)系統(tǒng)被廣泛應(yīng)用.其實(shí)際工作時(shí)在距離探測(cè)線圈一定距離處，鋪設(shè)一個(gè)或多個(gè)參考線圈，滿足探測(cè)線圈接收含噪的磁共振信號(hào)，參考線圈只接收噪聲的要求.利用探測(cè)線圈和參考線圈中噪聲的相關(guān)性，采用自適應(yīng)參考對(duì)消算法可有效實(shí)現(xiàn)噪聲壓制(田寶鳳等，2012).圖2為帶一個(gè)參考線圈的多通道磁共振探測(cè)系統(tǒng)的布線示意圖.

圖2　帶一個(gè)參考線圈的多通道磁共振探測(cè)系統(tǒng)示意圖Fig.2　Multi-channel MRS instrumentation with a reference loop

2.2MRS信號(hào)中噪聲分析及消噪對(duì)策

分析發(fā)現(xiàn)，影響MRS信號(hào)質(zhì)量的噪聲主要包括工頻諧波噪聲、尖峰噪聲和自然噪聲(林君等，2010).其中，自然噪聲類似于均勻白噪聲或高斯白噪聲，通過(guò)多次測(cè)量，采用疊加算法即可消除，而工頻諧波噪聲和尖峰噪聲則需要選取一定的消噪對(duì)策進(jìn)行剔除.

多通道核磁共振地下水探測(cè)系統(tǒng)可以依據(jù)噪聲相關(guān)性，采用自適應(yīng)參考對(duì)消方法消噪，消噪效果的好壞主要受噪聲相關(guān)程度的影響(Müller-Petke and Yaramanci,2011b).研究發(fā)現(xiàn)，探測(cè)線圈和參考線圈中的工頻諧波噪聲相關(guān)性較好，而尖峰噪聲相關(guān)性很差，導(dǎo)致自適應(yīng)參考對(duì)消方法難以有效發(fā)揮作用.因此綜合考慮，MRS信號(hào)的消噪對(duì)策應(yīng)首先剔除尖峰噪聲，再通過(guò)自適應(yīng)參考對(duì)消方法壓制工頻諧波干擾，最后采用疊加算法剔除自然噪聲影響，圖3給出了MRS信號(hào)的處理流程.

圖3　MRS信號(hào)處理流程圖Fig.3　Flow chart of MRS signal processing

3基于能量運(yùn)算的去尖峰算法

3.1信號(hào)能量計(jì)算表達(dá)式的推導(dǎo)與定義

用xn代表采樣后離散的余弦信號(hào)，用式(2)表示：

(2)

式中，Ω=2πf/fs，fs是采樣頻率；A和φ分別是信號(hào)的振幅和初始相位.假定A、Ω和φ三個(gè)參數(shù)恒定，并選擇與xn相鄰的xn-1和xn+1組成方程組，得到表達(dá)式(3)：

(3)

根據(jù)三角函數(shù)的積化和差公式對(duì)表達(dá)式(3)進(jìn)行變換，可以得出：

(4)

再根據(jù)三角函數(shù)的二倍角公式對(duì)表達(dá)式(4)進(jìn)行變換，可以得出：

xn+1xn-1=A2cos2(Ω n+φ)-A2sin2Ω，

(5)

(6)

整理得：

(7)

當(dāng)Ω<π/4，即f/fs<1/8，在誤差允許的范圍內(nèi)，

(8)

由于余弦信號(hào)的能量正比于振幅的平方與頻率的平方的乘積(Miller,1937)，因此可用表達(dá)式(8)表示余弦信號(hào)的能量，得到：

(9)

式中，En即為余弦信號(hào)xn的能量.

MRS信號(hào)是呈指數(shù)規(guī)律衰減的余弦信號(hào)，因此MRS信號(hào)的振幅A是隨時(shí)間按指數(shù)規(guī)律變化的.假設(shè)A=A0·e-a n，Ω=Ω0，φ=φ0，可得到MRS信號(hào)為：

(10)

利用信號(hào)能量表達(dá)式(9)計(jì)算表達(dá)式(10)所代表的MRS信號(hào)的能量，可以得到

(11)

同樣，當(dāng)采樣頻率大于信號(hào)頻率8倍時(shí)，滿足

(12)

3.2尖峰噪聲的檢測(cè)原理

當(dāng)有尖峰噪聲干擾時(shí)，假設(shè)xn是接收系統(tǒng)采集的信號(hào)，x1n是不含尖峰噪聲的MRS信號(hào)，x2n是尖峰噪聲.可知，xn是x1n和x2n的線性組合，滿足如下關(guān)系：

(13)

假設(shè)使

(14)

成立，可以得出：

(15)如果x1n和x2n不相關(guān)，用數(shù)學(xué)期望值求信號(hào)的能量，滿足如下關(guān)系：

(16)

引入延遲時(shí)間為2τ的x(t)的自相關(guān)函數(shù)，表達(dá)式為

(17)

式中，Rx(t+τ,t-τ)的傅里葉變換為

(18)

式中，W(t,ω)就是該信號(hào)的譜估計(jì).由于W(t,ω)是角頻率w的偶函數(shù)

(19)

式中，B是通頻帶寬.因此，利用公式(15)至公式(19)對(duì)表達(dá)式(14)進(jìn)行變換，可以得出：

=Rx(n,n)-Rx(n+1,n-1)

(20)

式中，(1-cos2ωT)起高通濾波器的作用.假設(shè)

(21)

可知，W(n,ω)經(jīng)過(guò)高通濾波器后是W′(n,ω).由式(20)可以得到

(22)其中

(23)

(24)

經(jīng)過(guò)上述能量運(yùn)算后，尖峰噪聲被突出，此時(shí)需要選取適當(dāng)?shù)拈撝颠M(jìn)行尖峰識(shí)別并剔除.

圖4a為一組受尖峰噪聲干擾的MRS信號(hào)，可以看出MRS信號(hào)是按指數(shù)規(guī)律變化的衰減信號(hào)，而尖峰噪聲在數(shù)據(jù)列中的表現(xiàn)形式為異常值，其數(shù)值多大于信號(hào)值.經(jīng)過(guò)能量運(yùn)算后信號(hào)如圖4b所示，可見(jiàn)其依然按指數(shù)規(guī)律衰減.信號(hào)和尖峰噪聲的這一特點(diǎn)，使得選取的閾值既不能把信號(hào)最大值標(biāo)記成異常值，又必須識(shí)別不明顯的異常值，即幅度值小的尖峰噪聲.因此，常規(guī)采用基于數(shù)據(jù)列標(biāo)準(zhǔn)差的閾值確定方法(吳石林和張玘，2010)，當(dāng)數(shù)據(jù)列中存在數(shù)值較大的異常值時(shí)，閾值受其影響會(huì)相對(duì)較高，進(jìn)而降低了對(duì)不明顯異常的檢測(cè)能力.

基于中位數(shù)的絕對(duì)偏差法(Median Absolute Deviation, MAD)(Hoaglin et al.,2000)是一種非常有效的統(tǒng)計(jì)學(xué)方法.對(duì)于一系列測(cè)量值為x1,x2,…,xn的數(shù)據(jù)列，計(jì)算MAD，應(yīng)先求某測(cè)量值相對(duì)于其所在數(shù)據(jù)列中位數(shù)的絕對(duì)偏差，再對(duì)該絕對(duì)偏差求中位數(shù)，定義為：

(25)

圖4　峰值檢測(cè)實(shí)例(a) 受尖峰噪聲干擾的MRS信號(hào)； (b) 信號(hào)能量及閾值.Fig.4　Example of spike detection(a) MRS signal disturbed by spike noises； (b) Energy of signal and threshold.

將能量域中的信號(hào)與閾值進(jìn)行比較，超過(guò)閾值數(shù)據(jù)視作尖峰噪聲并剔除.由于實(shí)際采集的信號(hào)包含隨機(jī)噪聲，隨機(jī)噪聲中的高頻成分在能量域中容易被誤識(shí)別為“小尖峰”，導(dǎo)致增加無(wú)效的判斷次數(shù)及計(jì)算量.因此，可適當(dāng)?shù)赝ㄟ^(guò)低通濾波器，去除隨機(jī)噪聲中的高頻成分，以準(zhǔn)確識(shí)別并剔除尖峰噪聲.

仿真一組含有16000個(gè)等精度采樣點(diǎn)、采樣頻率為66667 Hz的四次疊加信號(hào)，單次采集的原始信號(hào)如圖5a所示，圖中的尖峰噪聲、工頻諧波噪聲和高斯白噪聲的干擾時(shí)刻以及干擾強(qiáng)度均是隨機(jī)的.使用能量運(yùn)算法剔除尖峰噪聲主要包括三個(gè)步驟：

(1) 計(jì)算單次采集信號(hào)的能量，將信號(hào)由時(shí)域轉(zhuǎn)化到能量域.再將計(jì)算后的信號(hào)的能量通過(guò)低通濾波器，如圖5b所示.

(2) 采用MAD方法確定閾值.

(3) 將低通濾波后信號(hào)的能量與閾值進(jìn)行比較，如圖5c所示，超出閾值的部分視作尖峰噪聲，用該時(shí)刻無(wú)尖峰噪聲的數(shù)據(jù)列的中位數(shù)進(jìn)行替代，以達(dá)到剔除尖峰噪聲的目的.圖5d給出了剔除尖峰噪聲后的信號(hào).

4消噪算法仿真分析

目前在磁共振信號(hào)處理中通常采用統(tǒng)計(jì)疊加方法進(jìn)行尖峰噪聲剔除(Jiang et al.,2011)，其基本思想是把尖峰噪聲看作統(tǒng)計(jì)學(xué)中的粗大誤差，首先提取一個(gè)可疑的測(cè)量值，然后按t分布檢驗(yàn)提取的測(cè)量值是否含有粗大誤差，再?zèng)Q定其是否需要剔除，最終剔除后的數(shù)據(jù)參與到疊加運(yùn)算中.

圖5　基于能量運(yùn)算的去尖峰算法消噪流程圖(a) 單次采集的原始信號(hào)； (b) 計(jì)算后的信號(hào)能量； (c) 中位數(shù)絕對(duì)偏差法確定的閾值； (d) 剔除尖峰噪聲后的信號(hào).Fig.5　The workflow of energy calculation(a) Original signals of each sampled； (b) Energy of signal calculated by energy calculation； (c) Threshold determined by median absolute deviation；(d) MRS signals of each sampled without spikes.

分別采用能量運(yùn)算法和統(tǒng)計(jì)疊加法，對(duì)含有不同干擾幅度以及干擾持續(xù)時(shí)間尖峰噪聲的信號(hào)進(jìn)行處理，驗(yàn)證能量運(yùn)算法的應(yīng)用效果，并觀察多通道系統(tǒng)采集中剔除尖峰噪聲后，主通道和參考道數(shù)據(jù)相關(guān)性的改善情況.

除了能從波形上觀察尖峰噪聲的剔除效果，本文引入信噪比和標(biāo)準(zhǔn)差進(jìn)行數(shù)值衡量.信噪比的定義為：

(26)

式中，信噪比SNR采用了dB的表示方式.

標(biāo)準(zhǔn)差的定義為：

(27)

4.1干擾幅度不同的尖峰噪聲去噪仿真

圖6　統(tǒng)計(jì)疊加法和能量運(yùn)算法對(duì)不同幅度尖峰噪聲的去噪結(jié)果(a) 仿真的含噪信號(hào)； (b) 應(yīng)用統(tǒng)計(jì)疊加法剔除尖峰噪聲； (c) 應(yīng)用能量運(yùn)算法剔除尖峰噪聲.虛線是區(qū)分不同幅度尖峰噪聲的標(biāo)識(shí)，幅度值超過(guò)該虛線的視作大幅度尖峰噪聲，反之，視作小幅度尖峰噪聲.Fig.6　Statistical stacking approach and energy calculation approach to remove spikes of different amplitude(a) Simulated signals with spikes； (b) Remove spikes by statistical stacking approach； (c) Remove spikes by energy calculation approach. The dash line is the sign of different amplitude，data above the dash line represent big spikes，conversely，represent small spikes.

圖中具有衰減趨勢(shì)的虛線，其幅度是信號(hào)幅度的1.5倍，將它作為區(qū)別不同幅度尖峰噪聲的標(biāo)識(shí).幅度值超過(guò)該虛線的視作大幅度的尖峰噪聲；反之，視作小幅度的尖峰噪聲.分別采用統(tǒng)計(jì)疊加法和能量運(yùn)算法對(duì)圖6a所示的仿真含噪信號(hào)進(jìn)行消噪處理，圖6b和6c為對(duì)應(yīng)的消噪結(jié)果.從消噪結(jié)果圖中可以明顯看出，應(yīng)用統(tǒng)計(jì)疊加法后，幅值較大的尖峰噪聲被剔除，但幅值較小的尖峰噪聲仍然存在；而經(jīng)過(guò)能量運(yùn)算法后，幅值較大的尖峰噪聲和幅值較小的尖峰噪聲均被剔除.應(yīng)用統(tǒng)計(jì)疊加法后，信噪比SNR=6.6617 dB，提高了5.7432 dB.另外，標(biāo)準(zhǔn)差STD=80.1300 nV，與不含尖峰噪聲的MRS信號(hào)相比，標(biāo)準(zhǔn)差偏差為3.4610%.而應(yīng)用能量運(yùn)算法后，信噪比SNR=12.6133 dB，提高了11.6948 dB.標(biāo)準(zhǔn)差STD=77.7195 nV，標(biāo)準(zhǔn)差偏差為0.3486%.

圖7　統(tǒng)計(jì)疊加法和能量運(yùn)算法對(duì)不同持續(xù)時(shí)間尖峰噪聲的去噪結(jié)果(a) 仿真的含噪信號(hào)； (b) 應(yīng)用統(tǒng)計(jì)疊加法剔除尖峰噪聲； (c) 應(yīng)用能量運(yùn)算法剔除尖峰噪聲.Fig.7　Statistical stacking approach and energy calculation approach to remove spikes of different duration time(a) Simulated signals with spikes； (b) Remove spikes by statistical stacking approach；(c) Remove spikes by energy calculation approach.

結(jié)合統(tǒng)計(jì)疊加法和能量運(yùn)算法的消噪原理對(duì)其消噪差異進(jìn)行分析，可知，統(tǒng)計(jì)疊加法先剔除一個(gè)可疑的測(cè)量值，然后檢驗(yàn)剔除的測(cè)量值是否含有粗大誤差.在某些情況下，幅度值較小的干擾被認(rèn)為不含有粗大誤差，不予剔除.因此，統(tǒng)計(jì)疊加法能剔除幅度值較大的尖峰噪聲，卻常常忽略幅度值較小的尖峰噪聲.然而，經(jīng)過(guò)能量運(yùn)算后，原來(lái)時(shí)域中幅度值較小的尖峰噪聲在能量域中其幅度值會(huì)被增加，明顯地超過(guò)信號(hào)的幅度，很容易被識(shí)別并去除.

4.2持續(xù)時(shí)間不同的尖峰噪聲去噪仿真

分別采用統(tǒng)計(jì)疊加法和能量運(yùn)算法對(duì)圖7a所示的仿真含噪信號(hào)進(jìn)行消噪處理，圖7b和7c為對(duì)應(yīng)的消噪結(jié)果.應(yīng)用統(tǒng)計(jì)疊加法后，信噪比SNR=18.3899 dB，提高了18.1638 dB.另外，標(biāo)準(zhǔn)差STD=77.6652 nV，與不含尖峰噪聲的MRS信號(hào)相比，標(biāo)準(zhǔn)差偏差為0.2785%.而應(yīng)用能量運(yùn)算法后，信噪比SNR=18.3921 dB，提高了18.1660 dB.標(biāo)準(zhǔn)差STD=77.6744 nV，標(biāo)準(zhǔn)差偏差為0.2904%.從消噪結(jié)果圖以及消噪后的信噪比和標(biāo)準(zhǔn)差中可以明顯看出，統(tǒng)計(jì)疊加法和能量運(yùn)算法均可剔除干擾幅度較大，持續(xù)時(shí)間較長(zhǎng)的尖峰噪聲，且不損失有效的磁共振信號(hào)，標(biāo)準(zhǔn)差偏差控制在0.3%以內(nèi)，滿足實(shí)際應(yīng)用的需求.

4.3多通道系統(tǒng)數(shù)據(jù)去噪后提高信號(hào)相關(guān)性驗(yàn)證

模擬探測(cè)線圈采集含噪聲的磁共振信號(hào)，參考線圈只采集周圍環(huán)境噪聲，二者共同構(gòu)成多通道采集系統(tǒng).參考對(duì)消算法與通道間數(shù)據(jù)相關(guān)性有關(guān)，因此需引入相關(guān)性的概念.設(shè)Pxy(f)是信號(hào)x和信號(hào)y的互相關(guān)功率譜密度，Pxx(f)是信號(hào)x的自相關(guān)功率譜密度，Pyy(f)是信號(hào)y的自相關(guān)功率譜密度.相關(guān)性的定義為：

圖8　使用統(tǒng)計(jì)疊加方法和使用基于能量運(yùn)算的尖峰噪聲抑制方法剔除主通道中的尖峰噪聲(a) 統(tǒng)計(jì)疊加方法； (b) 能量運(yùn)算法； (c) 頻譜結(jié)果.Fig.8　Remove spikes by statistical stacking approach and energy calculation approach in detection loop(a) Statistical stacking approach； (b) Energy calculation approach； (c) Results in frequency domain.

(28)

式中，r介于0和1之間，當(dāng)r=1時(shí)，信號(hào)x和信號(hào)y完全相關(guān)；當(dāng)r=0時(shí)，信號(hào)x和信號(hào)y不相關(guān).

仿真模型參數(shù)如下：MRS信號(hào)的中心頻率fL=2326 Hz，工頻諧波干擾頻率f1=2250 Hz，f2=2350 Hz，參考道中工頻干擾的幅度是主通道中工頻干擾幅度的3倍，相位相同.主通道中隨機(jī)加入幅度為150 nV到550 nV、持續(xù)時(shí)間為4 ms到15 ms的五個(gè)尖峰噪聲干擾；同時(shí)，參考道中隨機(jī)加入幅度為300 nV到1000 nV、持續(xù)時(shí)間為5 ms到20 ms的八個(gè)尖峰噪聲干擾.使得主通道和參考道的數(shù)據(jù)相關(guān)性r=0.1369.

圖9　使用統(tǒng)計(jì)疊加方法和使用基于能量運(yùn)算的尖峰噪聲抑制方法剔除參考道中的尖峰噪聲(a) 統(tǒng)計(jì)疊加方法； (b) 能量運(yùn)算法； (c) 頻譜結(jié)果.Fig.9　Remove spikes by statistical stacking approach and energy calculation approach in reference loop(a) Statistical stacking approach； (b) Energy calculation approach； (c) Results in frequency domain.

分別使用統(tǒng)計(jì)疊加法和能量運(yùn)算法剔除主通道以及參考道中尖峰噪聲.經(jīng)過(guò)計(jì)算，應(yīng)用統(tǒng)計(jì)疊加法剔除尖峰噪聲后，主通道與參考道信號(hào)相關(guān)性r=0.4502，提高了0.3133.應(yīng)用能量運(yùn)算法剔除尖峰噪聲后，主通道與參考道信號(hào)相關(guān)性r=0.4941，提高了0.3572.圖8和圖9分別給出了統(tǒng)計(jì)疊加法和能量運(yùn)算法剔除尖峰噪聲前后主通道以及參考道信號(hào)的時(shí)域波形及頻譜.消噪前后信號(hào)的時(shí)域波形進(jìn)一步驗(yàn)證了4.1和4.2節(jié)得出的結(jié)論：統(tǒng)計(jì)疊加法和能量運(yùn)算法均可剔除干擾幅度較大，持續(xù)時(shí)間較長(zhǎng)的尖峰噪聲；另外，能量運(yùn)算法可以識(shí)別并剔除干擾幅度較小的尖峰噪聲.從消噪前后信號(hào)頻譜圖以及消噪后的相關(guān)性中可以明顯看出，采用能量運(yùn)算法剔除尖峰噪聲后，能夠顯著提高主通道與參考道信號(hào)相關(guān)性，為后續(xù)應(yīng)用自適應(yīng)噪聲對(duì)消算法更好地剔除工頻諧波噪聲奠定了基礎(chǔ).

5實(shí)測(cè)數(shù)據(jù)處理

為了驗(yàn)證算法的實(shí)用性，在長(zhǎng)春市郊區(qū)燒鍋鎮(zhèn)進(jìn)行了野外試驗(yàn).使用核磁共振地下水探測(cè)儀(林君等，2010)進(jìn)行信號(hào)采集，當(dāng)?shù)氐睦獱栴l率fL=2326 Hz，圖10a中藍(lán)色線為儀器采集的原始信號(hào)時(shí)域波形圖.從圖中可以看出：該信號(hào)含有較多的尖峰噪聲，并且從信號(hào)的時(shí)域波形上不能明顯看出MRS信號(hào)指數(shù)型的衰減趨勢(shì).經(jīng)計(jì)算，原始信號(hào)的信噪比SNR=-37.7090 dB.

依次操作3.3節(jié)中使用能量運(yùn)算法的三個(gè)步驟，對(duì)圖10a所示實(shí)測(cè)數(shù)據(jù)進(jìn)行尖峰噪聲剔除，得到消噪后的時(shí)域信號(hào)，如圖10a中紅色線所示.可見(jiàn)，原始信號(hào)中的尖峰噪聲被明顯剔除，消噪后的信號(hào)呈現(xiàn)指數(shù)形式衰減.再?gòu)念l域上進(jìn)行比較，如圖10b所示，消噪后尖峰噪聲的頻率分量被壓制，突出了信號(hào)的頻率.經(jīng)計(jì)算，消噪后信號(hào)的信噪比SNR1=-29.1380 dB，剔除尖峰噪聲后信噪比提高了ΔSNR1=8.5710 dB.

圖10　基于能量運(yùn)算的尖峰噪聲抑制方法剔除實(shí)測(cè)數(shù)據(jù)中的尖峰噪聲(a) 剔除尖峰噪聲前后的時(shí)域圖； (b) 剔除尖峰噪聲前后的頻域圖.Fig.10　Spikes in MRS experimental data removed by energy calculation approach(a) Remove spikes in time domain；(b) Remove spikes in frequency domain.

為了進(jìn)一步說(shuō)明方法的有效性，分別對(duì)不同信噪比下的另外兩組實(shí)測(cè)數(shù)據(jù)進(jìn)行處理，處理結(jié)果如圖11所示.由圖11消噪前后的時(shí)域信號(hào)及頻譜可以看出，對(duì)于不同信噪比實(shí)測(cè)數(shù)據(jù)，基于能量運(yùn)算的尖峰噪聲抑制方法可以獲得較好的消噪效果.

為了驗(yàn)證算法在多通道系統(tǒng)數(shù)據(jù)處理中的實(shí)用性，在長(zhǎng)春市凈月潭進(jìn)行了野外試驗(yàn).使用帶參考線圈的核磁共振地下水探測(cè)儀進(jìn)行信號(hào)采集，當(dāng)?shù)氐睦獱栴l率fL=2326 Hz，受到2250 Hz和2350 Hz兩個(gè)明顯的工頻噪聲干擾.使用能量運(yùn)算法剔除主通道和參考道中的尖峰噪聲，結(jié)果如圖12和圖13所示.從圖中可以看出：信號(hào)中的尖峰噪聲被剔除，并且去除尖峰噪聲后，工頻噪聲的頻率成分更加明顯.經(jīng)計(jì)算，剔除尖峰噪聲前，主通道和參考道的信號(hào)相關(guān)性r=0.1575；剔除尖峰噪聲后，主通道和參考道的信號(hào)相關(guān)性r=0.2481.可見(jiàn)，對(duì)多通道磁共振探測(cè)系統(tǒng)，先剔除主通道和參考道中的尖峰噪聲，提高信號(hào)相關(guān)性后，再應(yīng)用自適應(yīng)噪聲對(duì)消算法去除工頻諧波噪聲將會(huì)取得更好的去噪效果.

6結(jié)論

在強(qiáng)尖峰噪聲干擾下如何更好地實(shí)現(xiàn)微弱MRS信號(hào)可靠提取是核磁共振地下水探測(cè)抗干擾技術(shù)中的關(guān)鍵問(wèn)題，它將直接影響后續(xù)應(yīng)用自適應(yīng)參考對(duì)消算法對(duì)工頻諧波噪聲的有效濾除.能量運(yùn)算法可有效檢測(cè)尖峰噪聲并突出不易識(shí)別的小幅度尖峰噪聲.本文根據(jù)能量運(yùn)算法的這一特點(diǎn)，針對(duì)強(qiáng)尖峰噪聲干擾問(wèn)題，推導(dǎo)了能量域信號(hào)表達(dá)式，通過(guò)計(jì)算信號(hào)能量，有效識(shí)別尖峰噪聲，然后采用基于中位數(shù)的絕對(duì)偏差法確定閾值，實(shí)現(xiàn)對(duì)尖峰噪聲的有效剔除.與應(yīng)用效果較好的統(tǒng)計(jì)疊加法對(duì)比研究，通過(guò)數(shù)值仿真和實(shí)測(cè)數(shù)據(jù)處理，得出如下結(jié)論：

圖11　不同信噪比下實(shí)測(cè)數(shù)據(jù)處理結(jié)果Fig.11　Processing results of MRS signals on different SNRs

圖12　基于能量運(yùn)算的尖峰噪聲抑制方法剔除主通道中的尖峰噪聲(a) 剔除尖峰噪聲前后的時(shí)域圖； (b) 剔除尖峰噪聲前后的頻域圖.Fig.12　Spikes in detection loop removed by energy calculation approach(a) Remove spikes in time domain； (b) Remove spikes in frequency domain.

圖13　基于能量運(yùn)算的尖峰噪聲抑制方法剔除參考道中的尖峰噪聲(a) 剔除尖峰噪聲前后的時(shí)域圖； (b) 剔除尖峰噪聲前后的頻域圖.Fig.13　Spikes in reference loop removed by energy calculation approach(a) Remove spikes in time domain； (b) Remove spikes in frequency domain.

(1) 對(duì)于干擾幅度較大，持續(xù)時(shí)間較長(zhǎng)的尖峰噪聲，能量運(yùn)算法和統(tǒng)計(jì)疊加法均能識(shí)別并剔除，且不損失有效的磁共振信號(hào)，仿真數(shù)據(jù)去噪后信噪比提高了18 dB，標(biāo)準(zhǔn)差偏差控制在0.3%以內(nèi)，滿足實(shí)際應(yīng)用的需求.

(2) 對(duì)于小于信號(hào)幅度1.5倍的小幅度尖峰噪聲，能量運(yùn)算法可有效識(shí)別并剔除，而統(tǒng)計(jì)疊加法無(wú)法實(shí)現(xiàn).同一組仿真數(shù)據(jù)應(yīng)用統(tǒng)計(jì)疊加法剔除尖峰噪聲后，信噪比提高了5.7432 dB；應(yīng)用能量運(yùn)算法剔除尖峰噪聲后，信噪比提高了11.6948 dB.

(3) 針對(duì)多通道系統(tǒng)，仿真數(shù)據(jù)經(jīng)過(guò)能量運(yùn)算法剔除尖峰噪聲后，相關(guān)性提高了0.3572，可見(jiàn)剔除尖峰噪聲后可明顯提高主通道和參考道的數(shù)據(jù)相關(guān)性.

(4) 本文所述能量運(yùn)算法不受尖峰噪聲干擾幅度以及干擾持續(xù)時(shí)間限制.隨機(jī)抽取三組不同信噪比的實(shí)測(cè)數(shù)據(jù)，采用能量運(yùn)算法剔除尖峰噪聲后，信噪比均得到大幅度提升，獲得了很好的實(shí)際應(yīng)用效果.

(5) 針對(duì)受工頻噪聲干擾的多通道磁共振探測(cè)系統(tǒng)的實(shí)測(cè)數(shù)據(jù)，使用本文所述能量運(yùn)算法剔除主通道和參考道的尖峰噪聲后，可提高數(shù)據(jù)相關(guān)性.

在核磁共振地下水實(shí)際探測(cè)應(yīng)用中，除尖峰噪聲之外，MRS數(shù)據(jù)還會(huì)受到工頻諧波噪聲和其他未知噪聲等多種影響，尤其在城市周邊和村莊附近探測(cè)時(shí)，高強(qiáng)度的工頻噪聲干擾會(huì)將MRS信號(hào)淹沒(méi).鑒于此，帶有參考線圈的多通道探測(cè)系統(tǒng)被推出，自適應(yīng)參考對(duì)消算法成為MRS消噪算法的主流方向.在運(yùn)用本文能量運(yùn)算法對(duì)MRS信號(hào)進(jìn)行尖峰噪聲剔除后，主通道和參考道噪聲相關(guān)性的提高為后續(xù)自適應(yīng)參考對(duì)消算法的有效應(yīng)用奠定了基礎(chǔ).下一步將研究基于自適應(yīng)參考對(duì)消算法對(duì)工頻諧波噪聲進(jìn)行濾除的問(wèn)題.

References

Fabian N. 2010. Processing of full time series, multichannel surface NMR signals [Master′s thesis]. Switzerland: ETH Zürich.

Hoaglin D C, Mosteller F, Tukey J W. 2000. Understanding Robust and Exploratory Data Analysis. New York: Wiley.

Jiang C D, Lin J, Duan Q M, et al. 2011. Statistical stacking and adaptive notch filter to remove high-level electromagnetic noise from MRS measurements.NearSurf.Geophys., 9(5): 459-468.

Legchenko A, Valla P. 2002. A review of the basic principles for proton magnetic resonance sounding measurements.J.Appl.Geophys., 50(1-2): 3-19.

Legchenko A, Baltassat J M, Beauce A, et al. 2002. Nuclear magnetic resonance as a geophysical tool for hydrogeologists.J.Appl.Geophys., 50(1-2): 21-46.

Lin J, Duan Q M, Wang Y J, et al. 2010. Theory and Design of Magnetic Resonance Sounding Instrument for Groundwater Detection and Its Applications (in Chinese). Beijing: Science Press.

Lin J. 2010. Situation and progress of nuclear magnetic resonance technique for groundwater investigations.ProgressinGeophysics(in Chinese), 25(2): 681-691, doi: 10.3969/j.issn.1004-2903.2010.02.043. Lin J, Jiang C D, Duan Q M, et al. 2012. The situation and progress of magnetic resonance sounding for groundwater investigations and underground applications.JournalofJilinUniversity(EarthScienceEdition) (in Chinese), 42(5): 1560-1570.

Lu Q H, Wu T B, Lin J. 2009. A research report on development of instrument science for geophysics.ProgressinGeophysics(in Chinese), 24(2): 750-758.

Lubczynski M, Roy J. 2003. Hydrogeological interpretation and potential of the new magnetic resonance sounding (MRS) method.JournalofHydrology, 283(1-4): 19-40.

Mukhopadhyay S, Ray G C. 1998. A new interpretation of nonlinear energy operator and its efficacy in spike detection.IEEET.Biomed.Eng., 45(2): 180-187.

Müller-Petke M, Yaramanci U. 2010. Improving the signal-to-noise ratio of surface-NMR measurements by reference channel based noise cancellation.∥ 16thEAGE European Meeting of Environmental and Engineering Geophysics. EAGE.Müller-Petke M, Yaramanci U. 2011a. Noise Cancellation for surface NMR: Derivation of time and frequency domain approaches.∥ 17thEAGE European Meeting of Environmental and Engineering Geophysics. EAGE.

Müller-Petke M, Yaramanci U. 2011b. Noise cancellation for surface NMR-application of time and frequency domain approaches.∥ 17thEAGE European Meeting of Environmental and Engineering Geophysics. EAGE.

Miller D C. 1937. Sound Waves: Their Shape and Speed. New York: Macmillan Company Press.

Pan Y L, Zhang C D. 2000. Theories and Methods of Surface Nuclear Magnetic Resonance (in Chinese). Beijing: China University of Geosciences Press.Radic T. 2006. Improving the signal-to-noise ratio of surface NMR data due to the remote reference technique.∥ 12thEAGE European Meeting of Environmental and Engineering Geophysics. EAGE.Roy J, Lubczynski M. 2003. The magnetic resonance sounding technique and its use for groundwater investigations.Hydrogeol.J., 11(4): 455-465. Schirov M, Legchenko A, Creer G. 1991. A new direct non-invasive groundwater detection technology for Australia.ExplorationGeophysics, 22(2): 333-338.

Strehl S. 2006. Development of strategies for improved filtering and fitting of SNMR-signals. Germany: Technical University of Berlin, Institute of Applied Geosciences, Department of Applied Geophysics.

Strehl S, Rommel I, Hertrich M, et al. 2006. New strategies for filtering and fitting of MRS signals.∥ Proceedings 3rdInternational MRS Workshop. Madrid, Spain, 65-68. Tian B F, Lin J, Duan Q M, et al. 2012. Variable step adaptive noise cancellation algorithm for magnetic resonance sounding signal with a reference coil.ChineseJ.Geophys. (in Chinese), 55(7): 2462-2472, doi: 10.6038/j.issn.0001-5733.2012.07.030.

Walsh D O. 2008. Multi-channel surface NMR instrumentation and software for 1D/2D groundwater investigations.J.Appl.Geophys., 66(3-4): 140-150. Wang Z X, Rong L L, Lin J. 2009. Spike noise elimination from surface nuclear magnetic resonance signal for underground water.JournalofJilinUniversity(EngineeringandTechnologyEdition) (in Chinese), 39(5): 1282-1287.

Wu S L, Zhang Q. 2010. Error Analysis and Data Processing (in Chinese). Beijing: Tsinghua University Press.

Zhang R, Hu X Y, Yang D K, et al. 2006. Review of development of surface nuclear magnetic resonance.ProgressinGeophysics(in Chinese), 21(1): 284-289.

附中文參考文獻(xiàn)

林君, 段清明, 王應(yīng)吉等. 2010. 核磁共振找水儀原理與應(yīng)用. 北京: 科學(xué)出版社.

林君. 2010. 核磁共振找水技術(shù)的研究現(xiàn)狀與發(fā)展趨勢(shì). 地球物理學(xué)進(jìn)展, 25(2): 681-691, doi: 10.3969/j.issn.1004-2903.2010.02.043. 林君, 蔣川東, 段清明等. 2012. 復(fù)雜條件下地下水磁共振探測(cè)與災(zāi)害水源探查研究進(jìn)展. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 42(5): 1560-1570.

陸其鵠, 吳天彪, 林君. 2009. 地球物理儀器學(xué)科發(fā)展研究報(bào)告. 地球物理學(xué)進(jìn)展, 24(2): 750-758.

潘玉玲, 張昌達(dá). 2000. 地面核磁共振找水理論和方法. 北京: 中國(guó)地質(zhì)大學(xué)出版社.

田寶鳳, 林君, 段清明等. 2012. 基于參考線圈和變步長(zhǎng)自適應(yīng)的磁共振信號(hào)噪聲壓制方法. 地球物理學(xué)報(bào), 55(7): 2462-2472, doi: 10.6038/j.issn.0001-5733.2012.07.030.

王中興, 榮亮亮, 林君. 2009. 地面核磁共振找水信號(hào)中的奇異干擾抑制. 吉林大學(xué)學(xué)報(bào)(工學(xué)版), 39(5): 1282-1287. 吳石林, 張玘. 2010. 誤差分析與數(shù)據(jù)處理. 北京: 清華大學(xué)出版社. 張榮, 胡祥云, 楊迪琨等. 2006. 地面核磁共振技術(shù)發(fā)展述評(píng). 地球物理學(xué)進(jìn)展, 21(1): 284-289.

(本文編輯何燕)

基金項(xiàng)目國(guó)家重大科學(xué)儀器設(shè)備開(kāi)發(fā)和應(yīng)用專項(xiàng)(2011YQ030113)，國(guó)家自然科學(xué)基金面上項(xiàng)目(41374075)，吉林省科技重點(diǎn)攻關(guān)項(xiàng)目(20150519008JH，20140204022GX),中國(guó)博士后科學(xué)基金面上項(xiàng)目(2014M561296)聯(lián)合資助.

作者簡(jiǎn)介萬(wàn)玲，女，1986年生，2013年畢業(yè)于吉林大學(xué)，講師，主要從事磁共振信號(hào)處理方法研究.E-mail：wanling@jlu.edu.cn *通訊作者林婷婷，女，1983年生，2010年畢業(yè)于吉林大學(xué)，博士生導(dǎo)師，主要從事磁共振地下水探測(cè)方法研究.E-mail：ttlin@jlu.edu.cn

doi:10.6038/cjg20160631 中圖分類號(hào)P631

收稿日期2015-02-02，2016-05-03收修定稿

Spikes removal of magnetic resonance sounding data based on energy calculation

WAN Ling， ZHANG Yang， LIN Jun， JIANG Chuan-Dong， LIN Ting-Ting*

CollegeofInstrumentationandElectricalEngineering/KeyLaboratoryofGeo-ExplorationandInstrumentation,MinistryofEducation,JilinUniversity,Changchun130026,China

AbstractMagnetic Resonance Sounding (MRS) signal is extremely easy corrupted by the noise, especially by the harmonic noise and spike noise. Harmonic noise cancellation is often based on remote references and the adaptive noise cancellation(ANC)algorithm to increase the signal-to-noise ratio (SNR). However, ANC algorithm cannot play an effective role when spike noise exists. Because the spike noise often lead to a bad correlation between the detection loop and the remote reference loop, it is necessary to remove the spike noise before using the adaptive noise cancellation algorithm to cancel the harmonic noise.

KeywordsMagnetic Resonance Sounding signal； Noise cancellation algorithm； Energy calculation； Spike noise； Noise correlation

萬(wàn)玲， 張揚(yáng)， 林君等. 2016. 基于能量運(yùn)算的磁共振信號(hào)尖峰噪聲抑制方法.地球物理學(xué)報(bào)，59(6):2290-2301,doi:10.6038/cjg20160631.

Wan L, Zhang Y, Lin J, et al. 2016. Spikes removal of magnetic resonance sounding data based on energy calculation.ChineseJ.Geophys. (in Chinese),59(6):2290-2301,doi:10.6038/cjg20160631.