李 同,林 昊,田寶鳳,林 君
(吉林大學(xué)a.儀器科學(xué)與電氣工程學(xué)院;b.地球信息探測(cè)儀器教育部重點(diǎn)實(shí)驗(yàn)室,長(zhǎng)春130026)
地面核磁共振探測(cè)(SNMR:Surface Nuclear Magnetic Resonance)是20世紀(jì)80年代發(fā)展起來(lái)的一種非侵入性地下水資源探測(cè)方法。該方法通過(guò)產(chǎn)生激發(fā)磁場(chǎng)使水中氫質(zhì)子向高能級(jí)躍遷,激發(fā)磁場(chǎng)消失后氫質(zhì)子回歸平衡態(tài)的過(guò)程中將釋放能量,通過(guò)信號(hào)解析出地下環(huán)境中的一系列水文地質(zhì)參數(shù)[1]。但是,實(shí)際工程應(yīng)用中阻礙該方法得到大規(guī)模推廣的主要原因是核磁共振測(cè)得的信號(hào)很微弱,只有納伏級(jí),極易受環(huán)境噪聲中尖峰噪聲的干擾。尖峰噪聲的幅值往往為有效信號(hào)的數(shù)倍,混疊在一起會(huì)導(dǎo)致信號(hào)淹沒(méi)在其中無(wú)法辨識(shí)[2-4]。Dalgaard等[5]分析了核磁共振信號(hào)中各種噪聲的特征以及處理方案,指出數(shù)據(jù)處理過(guò)程中應(yīng)重視剔除尖峰干擾,否則在后續(xù)自適應(yīng)濾波過(guò)程中信號(hào)將產(chǎn)生異變,因此,在數(shù)據(jù)處理過(guò)程中應(yīng)首先對(duì)尖峰干擾進(jìn)行剔除。
針對(duì)尖峰噪聲的處理一般在探測(cè)時(shí)設(shè)置閾值,當(dāng)尖峰噪聲信號(hào)幅值超過(guò)閾值時(shí),則舍棄該組測(cè)量數(shù)據(jù),盡可能保留無(wú)干擾數(shù)據(jù)。然而,當(dāng)實(shí)驗(yàn)環(huán)境中有不間斷干擾源時(shí)需要大量重復(fù)實(shí)驗(yàn),因此,有必要進(jìn)行后期數(shù)據(jù)處理以解決該類問(wèn)題。針對(duì)此類問(wèn)題,Jiang等[6]提出了閾值疊加方法,對(duì)尖峰較少的數(shù)據(jù)效果較好,但其實(shí)現(xiàn)條件比較苛刻,尖峰越多計(jì)算復(fù)雜度越高,且只適用于數(shù)據(jù)處理的后期疊加過(guò)程。Strehl[7]提出了運(yùn)用小波變換的尖峰噪聲篩選方案,得到較好的效果,但算法的穩(wěn)定性較差。為此,筆者針對(duì)核磁共振探測(cè)應(yīng)用中的尖峰噪聲干擾問(wèn)題,提出數(shù)據(jù)處理前期的基于標(biāo)準(zhǔn)差中位數(shù)的尖峰篩選方案,經(jīng)較少次疊加采集后可有效剔除尖峰噪聲,以滿足核磁共振應(yīng)用的快速探測(cè)需要。
核磁共振地下水探測(cè)工作模式如圖1所示,地下水中的氫質(zhì)子經(jīng)過(guò)激發(fā)脈沖后在采集系統(tǒng)中形成自由感應(yīng)衰減信號(hào)(FID:Free Induction Decay)[8-10]通常表示為
其中E0為初始振幅;T*2為平均衰減時(shí)間常數(shù);ω0表示響應(yīng)信號(hào)的特征頻率;φ0為初始相位。
上述參數(shù)代表地下水的含水量,孔隙度等信息。其中E0一般為幾十到幾百納伏,而混入的噪聲干擾往往有幾千納伏,其中包括工頻干擾以及尖峰噪聲。工頻干擾主要來(lái)自電力線[11-13],尖峰噪聲是由雷電現(xiàn)象或工程機(jī)械運(yùn)轉(zhuǎn)時(shí)產(chǎn)生的瞬時(shí)電能,在儀器接收回路中產(chǎn)生的激發(fā)脈沖,主要特點(diǎn)是突發(fā)性,幅度高,隨機(jī)性強(qiáng)。圖2為一組采自農(nóng)安燒鍋鎮(zhèn)的混有尖峰噪聲的核磁共振信號(hào),可見(jiàn)其中混有兩處明顯的尖峰,相對(duì)微弱的信號(hào)被尖峰嚴(yán)重干擾,無(wú)法準(zhǔn)確顯示其衰減形態(tài),影響信號(hào)質(zhì)量。因此,筆者主要針對(duì)磁共振探測(cè)中的尖峰干擾去除問(wèn)題進(jìn)行研究。
圖1 核磁共振信號(hào)原理圖Fig.1 Surface nuclear magnetic resonance signal
圖2 實(shí)際采集混有尖峰干擾的信號(hào)Fig.2 Signal mixed with peak noise
以一組相同脈沖矩的N次疊加SNMR信號(hào)V(i,t)為例,其中i=1,2,3,…,N。首先通過(guò)Hilbert變換提取其包絡(luò)信號(hào)VENV(i,t)[14],再按疊加次數(shù)取其平均值得,進(jìn)行區(qū)域極大值判定,篩選信號(hào)中的尖峰噪聲。
當(dāng)t=T(n)時(shí),待處理的單次采集信號(hào)其中Tthreshold為閾值系數(shù),取值越小,則篩選條件越嚴(yán)苛;σ(t)為VENV(i,t)按采集次數(shù)求得的標(biāo)準(zhǔn)差;median{σ(t)}為取中值函數(shù)。
通過(guò)[V,C]=max{v(T)}找出尖峰噪聲的頂峰位置C,可得到尖峰的時(shí)域窗口 I,取值范圍為其中w為窗口寬度,用來(lái)估計(jì)尖峰噪聲的時(shí)域范圍。最后將鎖定的尖峰噪聲數(shù)據(jù)替換為疊加信號(hào)的均值
結(jié)合一般SNMR野外實(shí)驗(yàn)工作方法以及噪聲特點(diǎn)進(jìn)行仿真實(shí)驗(yàn)。依據(jù)式(1)生成一組32次采集的純凈SNMR信號(hào),其中單次采集信號(hào)如圖3所示。在此基礎(chǔ)上對(duì)每組數(shù)據(jù)同時(shí)混進(jìn)隨機(jī)噪聲及尖峰噪聲。圖4為32次疊加數(shù)據(jù)中的第7組數(shù)據(jù),該組信號(hào)共混有3處尖峰干擾,分別位于0.8 ms,110.3 ms和237.7 ms處,干擾幅值很大,位置無(wú)規(guī)律。
圖3 模擬的純凈核磁共振信號(hào)Fig.3 Synthetic signal
圖4 模擬的混有尖峰干擾的信號(hào)Fig.4 Synthetic signal mixed with peak noise
為消除尖峰噪聲的干擾,采用上述方法對(duì)數(shù)據(jù)集進(jìn)行處理,通過(guò)式(2)對(duì)多次采集的數(shù)據(jù)集中各時(shí)刻的比對(duì)篩選尖峰噪聲,其過(guò)程如圖5所示。圖5中黑色方塊分別表示圖4中3處尖峰時(shí)刻在各采集次數(shù)下的信號(hào)幅值,虛線的取值為式(2)不等式右側(cè)部分,代表篩選尖峰噪聲的下限參考值,所有大于該值的數(shù)據(jù)點(diǎn)被認(rèn)定為尖峰干擾,由此找出其時(shí)域坐標(biāo),進(jìn)而將其剔除,閾值Tthreshold的取值會(huì)影響識(shí)別尖峰干擾的準(zhǔn)確性。Tian等[15]討論了標(biāo)準(zhǔn)差系數(shù)與信噪比的適應(yīng)關(guān)系,認(rèn)為該系數(shù)設(shè)定取決于環(huán)境噪聲幅度。當(dāng)環(huán)境噪聲較弱時(shí),信噪比高,調(diào)整較低的閾值系數(shù)便可篩選尖峰干擾;當(dāng)環(huán)境噪聲較強(qiáng)時(shí),信噪比較低,此時(shí)閾值系數(shù)若不進(jìn)行調(diào)整,會(huì)將有效的核磁共振信號(hào)誤認(rèn)為尖峰干擾。因此,處理低信噪比數(shù)據(jù)時(shí)應(yīng)適當(dāng)增大閾值,根據(jù)噪聲水平,這里閾值系數(shù)取值5,可見(jiàn)在此閾值條件下能準(zhǔn)確地篩選各處尖峰干擾。
圖5 多次采集下閾值為5時(shí)的尖峰干擾篩選效果Fig.5 Recognizing peak noise when threshold is 5
該組數(shù)據(jù)經(jīng)過(guò)尖峰噪聲剔除的波形如圖6所示,圖7a,圖7b分別為該組數(shù)據(jù)處理前后的包絡(luò)信號(hào)。對(duì)比圖5可以看出,該組信號(hào)受到3組脈沖干擾明顯,幅值達(dá)到相應(yīng)時(shí)刻平均值的10倍。合理的閾值設(shè)定可以準(zhǔn)確地定位尖峰干擾,將其篩除,并很好地保留了其余信號(hào),可明顯看出尖峰剔除后信號(hào)的衰減形態(tài)。
圖6 尖峰剔除后的仿真信號(hào)Fig.6 Signal after deleting peak noise
上述仿真實(shí)驗(yàn)證實(shí)了尖峰剔除方法的有效性,可將該方法應(yīng)用于實(shí)際采集核磁共振信號(hào)的處理過(guò)程中。圖2是采集于農(nóng)安縣燒鍋鎮(zhèn)的核磁共振探測(cè)信號(hào),共進(jìn)行32次疊加采集,圖示為尖峰干擾較明顯的第6次疊加數(shù)據(jù),圖8、圖9給出了該組信號(hào)經(jīng)過(guò)尖峰剔除后的結(jié)果,原數(shù)據(jù)中存在兩處明顯尖峰,分別在19.92 ms以及75.92 ms處。由圖9可看出,經(jīng)過(guò)計(jì)算處理,兩處尖峰干擾均被剔除。
圖9 尖峰剔除前后的實(shí)測(cè)信號(hào)包絡(luò)曲線對(duì)比Fig.9 Envelope of field signal before and after filtering
通過(guò)包絡(luò)數(shù)據(jù)可以看出,尖峰干擾破壞了信號(hào)的連續(xù)衰減形態(tài),經(jīng)過(guò)計(jì)算處理后尖峰干擾可被有效濾除,證明了該算法可以解決實(shí)際采集信號(hào)中的尖峰干擾問(wèn)題,有較好的實(shí)用性。
由于尖峰干擾的存在,增大了后續(xù)數(shù)據(jù)處理的計(jì)算難度,而傳統(tǒng)的尖峰噪聲濾波方案在逐漸發(fā)展的核磁共振應(yīng)用環(huán)境中受到局限。筆者設(shè)計(jì)了基于重復(fù)采集數(shù)據(jù)的標(biāo)準(zhǔn)差中位數(shù)濾波方案,針對(duì)核磁共振信號(hào)中尖峰干擾進(jìn)行篩選剔除。該方案具有計(jì)算簡(jiǎn)單、準(zhǔn)確度高等優(yōu)勢(shì)。通過(guò)對(duì)仿真核磁共振信號(hào)的處理可以看出,該算法可準(zhǔn)確剔除混在信號(hào)中的3處尖峰干擾,同時(shí)保留真實(shí)信號(hào),尤其信號(hào)初段的尖峰幅值在真實(shí)信號(hào)的10倍以上,對(duì)其有效濾除后還原了真實(shí)信號(hào)的衰減形態(tài),對(duì)實(shí)際采集的核磁共振信號(hào)的處理也得到了同樣的效果,從而驗(yàn)證了該方法的有效性。
[1]林君.現(xiàn)代科學(xué)儀器及其發(fā)展趨勢(shì)[J].吉林大學(xué)學(xué)報(bào):信息科學(xué)版,2002,20(1):1-7.LIN Jun.Developing Trend of Modern Scientific Instrument[J].Journal of Jilin University:Information Science Edition,2002,20(1):1-7.
[2]YARAMANCI U.New Technologies in Groundwater Exploration Surface Nuclear Magnetic Resonance[J].Geologica Acta:An International Earth Science Journal,2004,2(2):109-120.
[3]JUAN PLATA,F(xiàn)ELIX RUBIO.MRS Experiments in a Noisy Area of a Detrital Aquifer in the South of Spain[J].Journal of Applied Geophysics,2002,50(1/2):83-94.
[4]UGUR YARAMANCI,GERHARD LANGE,MARIAN HERTRICH.Aquifer Characterisation Using Surface NMR Jointly with other Geophysical Techniques at the Nauen/Berlin Test Site[J].Journal of Applied Geophysics,2002,50(1/2):47-65.
[5]DALGAARD E,AUKEN E,LARSEN J J.Adaptive Noise Cancelling of Multichannel Magnetic Resonance Sounding Signals[J].Geophysical Journal International,2012,191(1):88-100.
[6]JIANG Chuandong,LIN Jun,DUAN Qingming,et al.Statistical Stacking and Adaptive Notch Filter to Remove High-Level Electromagnetic Noise from MRS Measurements[J].Near Surface Geophysics,2011,9(5):459-468.
[7]STREHL S.Development of Strategies for Improved Filtering and Fitting of Snmr Signals[D].Berlin:Institute of Applied Geosience,Berlin University of Technology,2006:57-64.
[8]林君,段清明,王應(yīng)吉,等.核磁共振找水儀原理與應(yīng)用[M].北京:科學(xué)出版社,2011:2-3.LIN Jun,DUAN Qingming,WANG Yingji,et al.Theory and Design of Magnetic Resonance Sounding Instrument for Groundwater Detection and Its Applications[M].Beijing:Science Press,2011:2-3.
[9]翁愛(ài)華,王雪秋,劉國(guó)興,等.導(dǎo)電性影響的地面核磁共振反演[J].地球物理學(xué)報(bào),2007,50(3):890-896.WENG Aihua,WANG Xueqiu,LIU Guoxing,et al.Nonlinear Inversion of Surface Nuclear Magnetic Resonance over Electrically Conductive Medium[J].Chinese Journal of Geophysics(in Chinese),2007,50(3):890-896.
[10]潘玉玲,張昌達(dá).地面核磁共振找水理論和方法[M].武漢:中國(guó)地質(zhì)大學(xué)出版社,2000:37.PAN Yuling,ZHANG Changda.The Theory and Method of SNMR [M].Wuhan:China University of Geosciences Press,2000:37.
[11]林婷婷,慧芳,蔣川東,等.分層多指數(shù)磁共振弛豫信號(hào)反演方法研究 [J].地球物理學(xué)報(bào),2013,56(8):2849-2861.LIN Tingting,HUI Fang,JIANG Chuandong,et al.Layered Multi-Exponential Inversion Method on Surface Magnetic Resonance Sounding Dataset[J].Chinese Journal Geophysics,2013,56(8):2849-2861.
[12]易曉峰,李鵬飛,林君,等.基于多匝環(huán)形線圈的核磁共振信號(hào)響應(yīng)計(jì)算與試驗(yàn)研究[J].地球物理學(xué)報(bào),2013,56(7):2484-2493.YI Xiaofeng,LI Pengfei,LIN Jun,et al.Simulation and Experimental Research of MRS Response Based on Multi-Turn Loop[J].Chinese Journal Geophysics,2013,56(7):2484-2493.
[13]田寶鳳,段清明.核磁共振信號(hào)工頻諧波的自適應(yīng)濾波方法 [J].吉林大學(xué)學(xué)報(bào):信息科學(xué)版,2009,27(3):223-228.TIAN Baofeng,DUAN Qingming.Removal Method of Industrial Frequency Harmonics in Nuclear Magnetic Resonance Signal Based on Adaptive Filter[J].Journal of Jilin University:Information Science Edition,2009,27(3):223-228.
[14]張旻,程家興,樊甫華,等.利用 Hilbert變換提取信號(hào)瞬時(shí)特征參數(shù)的問(wèn)題研究[J].電訊技術(shù),2003,43(4):44-48.ZHANG Min,CHENG Jiaxing,F(xiàn)AN Fuhua,et al.Study on the Problems in ExtractingInstantaneous Characters of Signals Based on Hilbert Transform [J].Telecommunication Engineering,2003,43(4):44-48.
[15]TIAN Baofeng,ZI Yanyong,ZHOU Yuanyuan,et al.Research on Pulse Noise Elimination Algorithm of Magnetic Resonance Sounding Full-Wave Signal [C]∥ 10th International Computer Conference on Wavelet Active Media Technology and Information Processing.Chengdu,China:IEEE,2013,6-10.