基于改進(jìn)短時(shí)傅里葉變換的磁共振隨機(jī)噪聲消減方法*

林婷婷 李玥 高興 萬玲?

1) (吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院, 長春 130026)

2) (吉林大學(xué), 地球信息探測儀器教育部重點(diǎn)實(shí)驗(yàn)室, 長春 130026)

磁共振測深法(magnetic resonance sounding, MRS)具有無需鉆探即可直接探測地下水含量的優(yōu)勢, 但是極低的信噪比(signal-to-noise ratio, SNR)限制了該方法的大范圍應(yīng)用, 目前的研究工作主要集中在消除MRS信號(hào)中的尖峰噪聲和工頻諧波噪聲上, 而隨機(jī)噪聲由于其無規(guī)律性導(dǎo)致難以消除, 但是它的影響不容忽視.目前MRS隨機(jī)噪聲的消減常采用疊加法, 需要重復(fù)采集數(shù)據(jù)疊加平均來達(dá)到消噪目的, 探測時(shí)間長且消噪效果有限.針對(duì)這一問題, 本文提出了一種改進(jìn)的短時(shí)傅里葉變換方法, 該方法通過處理單次采集的MRS包絡(luò)信號(hào)來降低數(shù)據(jù)量, 提高數(shù)據(jù)處理效率.改進(jìn)的短時(shí)傅里葉變換方法采用解析信號(hào)代替常規(guī)短時(shí)傅里葉變換中的實(shí)值信號(hào), 提高M(jìn)RS信號(hào)時(shí)頻域瞬時(shí)幅度的準(zhǔn)確度, 得到MRS信號(hào)的高精度時(shí)頻分布, 然后提取時(shí)頻域峰值幅度和峰值相位重構(gòu)信號(hào)來消除隨機(jī)噪聲.仿真實(shí)驗(yàn)和實(shí)測數(shù)據(jù)處理結(jié)果表明, 該方法能夠直接處理單次采集數(shù)據(jù), 在信噪比高于-17.21 dB的情況下可有效提取MRS信號(hào), 實(shí)現(xiàn)隨機(jī)噪聲的壓制, 且與傳統(tǒng)疊加法相比, 信噪比最多可提高27.88 dB, 均方根誤差最多縮小36.44倍, 參數(shù)估計(jì)值更加準(zhǔn)確.本文的研究結(jié)果為利用MRS獲取準(zhǔn)確的地下水分布情況奠定了良好的基礎(chǔ).

1 引 言

磁共振測深法(magnetic resonance sounding,MRS)是一種用于水文地質(zhì)調(diào)查的非侵入式地球物理方法[1-3], 由于其對(duì)質(zhì)子具有獨(dú)特敏感性,因此可以直接定量地估計(jì)地下水含量和孔隙結(jié)構(gòu)[4-6].近年來, 這種方法在測井[7]、實(shí)驗(yàn)室研究[8]和野外測量[9]中得到了快速發(fā)展.但是MRS在地磁場(0.05-0.06 mT)中進(jìn)行測量時(shí), 微弱的MRS信號(hào)易被環(huán)境噪聲淹沒, 導(dǎo)致信噪比(signal-tonoise ratio, SNR)降低, 水文參數(shù)估計(jì)錯(cuò)誤[10,11].

為了提高測量結(jié)果的準(zhǔn)確性, 通常情況下會(huì)采用信號(hào)增強(qiáng)和噪聲抑制兩種途徑來提高信噪比.Davis等[12,13]采用善于檢測微弱信號(hào)的超導(dǎo)量子干涉器件(SQUID)代替?zhèn)鹘y(tǒng)接收線圈來提高磁共振儀器靈敏度; 2016年, Grunewald等[9]通過在滿足絕熱條件時(shí)增加扳倒角來提高地下水中氫質(zhì)子磁化強(qiáng)度的有效分量;2018年, Lin等[14]在磁共振測深法中使用直流預(yù)極化場增強(qiáng)氫質(zhì)子宏觀磁矩,測得了冰層下的流動(dòng)水信號(hào).這些方法均在一定程度上提高了磁共振信號(hào)強(qiáng)度, 但是環(huán)境噪聲的影響仍不容小覷.干擾MRS信號(hào)的環(huán)境噪聲主要包括尖峰噪聲、工頻諧波噪聲和隨機(jī)噪聲, 在消噪過程中需要對(duì)這些噪聲進(jìn)行分類逐步消除, 常規(guī)步驟為: 1)去尖峰; 2)消除工頻諧波噪聲; 3)抑制隨機(jī)噪聲[3,15].目前在MRS領(lǐng)域行之有效的降噪方法主要針對(duì)尖峰噪聲和工頻諧波噪聲, 因?yàn)樗鼈冇商囟ㄔ肼曉串a(chǎn)生, 具有一定的規(guī)律性, 可以總結(jié)出對(duì)應(yīng)的數(shù)學(xué)表達(dá)式從而進(jìn)行消除.Jiang等[16]于2011年提出了一種基于羅曼諾夫斯基準(zhǔn)則的統(tǒng)計(jì)疊加方法來丟棄尖峰噪聲; Dalgaard等[17]于2012年實(shí)現(xiàn)了基于非線性能量算子的尖峰噪聲檢測算法; Larsen[18]于2016年的一項(xiàng)研究表明, 由電圍網(wǎng)產(chǎn)生的尖峰噪聲可建模為可調(diào)諧四階帶通濾波器的脈沖激勵(lì)并進(jìn)行消減.陷波濾波法由Legchenko和Valla[19]于2003年提出, 是消除單通道信號(hào)中工頻諧波噪聲的最常見并有效的方法.之后, 基于Walsh[20]于2008年發(fā)明的多通道儀器, 遠(yuǎn)程參考技術(shù)被應(yīng)用于工頻諧波噪聲的消除[17,20,21].當(dāng)存在多個(gè)噪聲源時(shí), Larsen等[22]于2014年提出了一種基于建模的降噪策略, 與多通道維納濾波相結(jié)合來降低多源干擾.然而隨機(jī)噪聲是一種來源廣泛且不確定的前后獨(dú)立平衡隨機(jī)過程, 由于其無規(guī)律性、無法預(yù)測并且可能與有效信號(hào)發(fā)生頻譜重疊而無法得到有效抑制.目前最常見的隨機(jī)噪聲消除方法為疊加法[2,17,23], 通過對(duì)同一地點(diǎn)多次采集信號(hào)進(jìn)行疊加平均來提高信噪比, 但是該方法消耗時(shí)間長, 儀器探測效率低, 并且只能抑制部分隨機(jī)噪聲.林婷婷等[24]提出了一種分段時(shí)頻峰值濾波法來處理磁共振全波信號(hào)中的隨機(jī)噪聲.磁共振全波信號(hào)采樣率高, 能夠提供更加豐富的地下含水層信息,但在實(shí)際探測過程中, 為了實(shí)現(xiàn)實(shí)時(shí)觀測而對(duì)探測速度要求較高時(shí), 仍然需要采集包絡(luò)信號(hào)并對(duì)其中的隨機(jī)噪聲進(jìn)行抑制.

為了簡化MRS探測中的隨機(jī)噪聲消減方法,實(shí)現(xiàn)隨機(jī)噪聲的即時(shí)處理, 本文提出了一種改進(jìn)的短時(shí)傅里葉變換(modified short-time Fourier transform, MSTFT)方法, 應(yīng)用于包含地下水探測關(guān)鍵參數(shù)的MRS包絡(luò)信號(hào)中來消減隨機(jī)噪聲.采集包絡(luò)信號(hào)可以降低實(shí)測數(shù)據(jù)總量, 實(shí)現(xiàn)儀器的快速測量.該方法首先對(duì)MRS包絡(luò)信號(hào)作希爾伯特變換獲取解析信號(hào), 然后用解析信號(hào)代替MRS實(shí)值信號(hào)作對(duì)稱短時(shí)傅里葉變換, 得到高精度的MRS信號(hào)時(shí)頻分布, 最后提取峰值幅度和峰值相位重構(gòu)信號(hào), 消除隨機(jī)噪聲.本文結(jié)構(gòu)安排如下: 首先分析MRS信號(hào)中的隨機(jī)噪聲和傳統(tǒng)疊加法消除隨機(jī)噪聲的局限性; 然后介紹本文提出的改進(jìn)短時(shí)傅里葉變換算法的原理, 并闡述該算法提高M(jìn)RS信號(hào)時(shí)頻分布精度的原因; 最后通過仿真實(shí)驗(yàn)和處理實(shí)測數(shù)據(jù), 驗(yàn)證該算法能夠有效抑制隨機(jī)噪聲, 在不同信噪比下提取準(zhǔn)確的信號(hào)參數(shù).

2 磁共振信號(hào)中的隨機(jī)噪聲

原始MRS信號(hào)是由受交變磁場激發(fā)的氫質(zhì)子在拉莫爾頻率處震蕩產(chǎn)生的指數(shù)衰減信號(hào), 其計(jì)算公式為

其中t是采樣時(shí)間, E0是MRS信號(hào)初始振幅, 與含水量成正比,為橫向弛豫時(shí)間, 表征含水層介質(zhì)孔隙度, 拉莫爾頻率 fL由地磁場強(qiáng)度決定, 會(huì)隨地理位置的變換而稍有改變, φ0表示由導(dǎo)電介質(zhì)和儀器相移產(chǎn)生的初始相位信息, q為激發(fā)脈沖矩.

通常情況下, 測量(1)式中的MRS全波信號(hào)需要高達(dá)幾十kHz的采樣頻率, 記錄的數(shù)據(jù)量非常大.為了減少對(duì)數(shù)據(jù)采集的需求, 提高探測效率, 現(xiàn)有的儀器通常采用雙通道正交探測系統(tǒng)輔以低通濾波器來探測MRS包絡(luò)信號(hào)(JLMRS-I(Jiang et al.2011), NumisPlus和NumisLite(Bernard 2006))[16,25].

其中, Eenvelope(t) 表示MRS包絡(luò)信號(hào), Δ f 為拉莫爾頻率和激發(fā)電場之間的頻率偏移, 理想條件下Δf 的值趨近于零. σx(t) 和 σy(t) 分別為MRS信號(hào)的同相噪聲和異相噪聲.

MRS信號(hào)中的隨機(jī)噪聲來源于寬帶背景噪聲和儀器電子元件噪聲, 圖1所示為對(duì)一組野外實(shí)測記錄分別進(jìn)行1, 4, 8和16次疊加法處理的結(jié)果.常用的磁共振測深法工作方案為先測量空采噪聲,然后再發(fā)射電流采集信號(hào), 因此該組實(shí)測記錄包含噪聲數(shù)據(jù)和含噪信號(hào)兩部分, 分別如圖1中灰色曲線和黑色曲線所示.測量地點(diǎn)的拉莫爾頻率為2330 Hz, 尖峰噪聲和工頻諧波噪聲在使用疊加法之前已經(jīng)消除.如圖1(a)和圖1(b)所示, 1次疊加時(shí)信號(hào)的時(shí)間序列被噪聲破壞, 頻譜中幾乎看不到信號(hào)的頻率響應(yīng).從圖1(c)和圖1(d)可以看出,當(dāng)疊加次數(shù)為4時(shí), 隨機(jī)噪聲被抑制, 信號(hào)明顯增強(qiáng).隨著疊加次數(shù) N s 增大, 噪聲水平越來越小.8次疊加(圖1(e)和圖1(f))和16次疊加(圖1(g)和圖1(h))的結(jié)果對(duì)比表明, 雖然SNR隨著疊加次數(shù)的增大而增大了, 但是隨機(jī)噪聲仍然存在.

圖1 疊加法分別疊加1, 4, 8和16次時(shí)消除隨機(jī)噪聲效果時(shí)域和頻域圖Fig.1.Time-domain and frequency-domain diagrams after stacking 1, 4, 8 and 16 times to suppress random noise, respectively.

假設(shè)Si和σi分別為單次記錄的信號(hào)和噪聲,那么經(jīng)Ns次疊加后的信噪比為

3 基于MSTFT的隨機(jī)噪聲抑制方法

3.1 短時(shí)傅里葉變換

實(shí)值信號(hào) x (t) 的短時(shí)傅里葉變換形式為

F(τ,f) 為信號(hào)x (t)w(t-τ) 的頻譜, 包含信號(hào)x(t)w(t-τ) 的初始幅度和初始相位信息.窗函數(shù)w(t) 滑動(dòng)截取信號(hào), 分別在頻率域形成以f為中心的窄帶, 時(shí)間窗函數(shù)越短, 短時(shí)傅里葉變換對(duì)信號(hào)的時(shí)間分辨率越高, x (t)w(t-τ) 的初始幅度和初始相位可近似看作信號(hào) x (t) 在固定窗內(nèi)的瞬時(shí)幅度和瞬時(shí)相位, 則 F (τ,f) 可表示信號(hào) x (t) 隨時(shí)間變化的頻譜.但是相應(yīng)的頻域分辨率會(huì)降低, 因此由Heisenberg不確定準(zhǔn)則的限制, 窗函數(shù) w (t) 的面積要不小于2[26,27], 若定義

則 As(τ,f) 和 θs(τ,f) 可分別表示信號(hào)的時(shí)頻幅度譜和時(shí)頻相位譜.

3.2 改進(jìn)的短時(shí)傅里葉變換

為了準(zhǔn)確地描述非平穩(wěn)信號(hào), 獲得高精度的時(shí)頻振幅譜和時(shí)頻相位譜, 本文在STFT的基礎(chǔ)上,以解析信號(hào)代替實(shí)信號(hào), 提出了一種MSTFT方法來提取MRS信號(hào), 消除隨機(jī)噪聲, 其變換形式為

其中, h (t) 為 x (t) 的解析信號(hào), 由(8)式給出.

同樣地, G (τ,f) 表示信號(hào) h (t)w(t-τ) 的頻譜, 如果定義

則 A (τ,f) 和 θ (τ,f) 分別表示信號(hào)的高精度時(shí)頻幅度譜和高精度時(shí)頻相位譜, A (τ,fτ) 和 θ (τ,fτ) 為給定時(shí)間 τ 處信號(hào)頻譜的峰值幅度和峰值相位, fτ為取得峰值幅度和峰值相位處的頻率.由于峰值幅度和峰值相位受隨機(jī)噪聲影響小, 因此可以用來重構(gòu)信號(hào):

r (t) 為抑制隨機(jī)噪聲影響后的MRS信號(hào), 其中Re[] 表示選取復(fù)值信號(hào)的實(shí)數(shù)部分.

通常情況下傅里葉變換是在頻域內(nèi)計(jì)算的, 為了更好地理解MSTFT與常規(guī)短時(shí)傅里葉變換之間的區(qū)別, 將變換形式改寫為頻域形式:

其中, H (ω) , X (ω) 和 W (ω) 分別為 h (t) , x (t) 和w(t) 的傅里葉變換形式, 且窗函數(shù) W (ω) 滿足

若以窗函數(shù) W1(ω) 代替 W (ω) , 則實(shí)值信號(hào)

表示信號(hào) x (t) 的零相帶通濾波響應(yīng), G (τ,f) 為R(τ,f) 的解析信號(hào), 因此采用MSTFT計(jì)算得到的 G (τ,f) 包含信號(hào) R (τ,f) 的瞬時(shí)幅度和瞬時(shí)相位信息.但對(duì)于常規(guī)的短時(shí)傅里葉變換 F (τ,f) , 由Hilbert變換和Fourier變換的性質(zhì)可知, 只有滿足

時(shí), 能得到

這時(shí)常規(guī)短時(shí)傅里葉變換能夠準(zhǔn)確獲取信號(hào) R(τ,f)的瞬時(shí)相位, 但是瞬時(shí)幅度僅為真實(shí)值的一半, 且在或時(shí), (4)式所描述的短時(shí)傅里葉變換獲取信號(hào) R (τ,f) 的瞬時(shí)幅度和瞬時(shí)相位均為近似值.因此, 在使用相同時(shí)間窗函數(shù)的情況下, MSTFT比常規(guī)短時(shí)傅里葉變換獲取的信號(hào)參數(shù)變化情況更為精確.

3.3 信號(hào)實(shí)例

為了說明MSTFT與常規(guī)傅里葉變換之間的區(qū)別, 以及MSTFT消除隨機(jī)噪聲提取MRS信號(hào)的過程, 仿真一組受隨機(jī)噪聲干擾的MRS信號(hào),信號(hào)參數(shù)E0= 200 nV,= 150 ms, fL= 2330 Hz,信號(hào)持續(xù)時(shí)間500 ms.加入的隨機(jī)噪聲均值為0,標(biāo)準(zhǔn)差為50 nV, 使信噪比為0.51 dB.以拉莫爾頻率對(duì)信號(hào)重采樣, 信號(hào)長度為1165.選取對(duì)稱窗函數(shù)為定義在 [ -L/2,L/2] 區(qū)間內(nèi)的海明窗, 其中L為窗長, 根據(jù)Heisenberg不確定準(zhǔn)則和MRS信號(hào)特征, 選取窗長 L =6 .對(duì)包絡(luò)信號(hào)和其解析信號(hào)分別進(jìn)行對(duì)稱窗短時(shí)傅里葉變換, 結(jié)果如圖2所示.圖2(a)和圖2(c)分別為時(shí)頻幅度譜 As(τ,f) 和高精度時(shí)頻幅度譜 A (τ,f) , 從圖中可知 As(τ,f) 和A(τ,f) 的變化趨勢基本一致, 但是 A (τ,f) 的變化范圍在0-300 nV之間, 與仿真信號(hào)的瞬時(shí)幅度一致, 而 As(τ,f) 在0-150 nV之間變化.圖2(e)描 述 了 As(τ,f) 和 A (τ,f) 之間 的 差異, 與 圖2(a)的變化趨勢相近, 說明 A (τ,f) 約為 As(τ,f) 的2倍.圖2(b)和圖2(d)為時(shí)頻相位譜 θs(τ,f) 和高精度時(shí)頻相位譜 θ (τ,f) , 二者之間的差異如圖2(f)所示, 可以看出 θs(τ,f) 和 θ (τ,f) 的差異幾乎為0, 與(16)式得出的結(jié)論相同.

圖2 MSTFT與常規(guī)短時(shí)傅里葉變換之間的區(qū)別 (a) 時(shí)頻幅度譜 A s(τ,f) ; (b) 時(shí)頻相位譜 θ s(τ,f) ; (c) 高精度時(shí)頻幅度譜A(τ,f) ; (d) 高精度時(shí)頻相位譜 θ (τ,f) ; (e) A s(τ,f) 和 A (τ,f) 的差異; (f) θ s(τ,f) 和 θ (τ,f) 的差異Fig.2.Difference between the MSTFT and the conventional short-time Fourier transform: (a) time-frequency amplitude spectrum As(τ,f) ; (b) time-frequency phase spectrum θ s(τ,f) ; (c) high-precision time-frequency amplitude spectrum A (τ,f) ; (d) high-precision time-frequency phase spectrum θ (τ,f) ; (e) difference between A s(τ,f) and A (τ,f) ; (f) difference between θ s(τ,f) and θ(τ,f) .

圖3 (a)為含噪信號(hào)的時(shí)域波形, 其瞬時(shí)幅度(圖3(b))和瞬時(shí)相位(圖3(c))受隨機(jī)噪聲影響嚴(yán)重, 無法準(zhǔn)確提取原信號(hào)信息.采用前文所述的MSTFT方法獲取高精度時(shí)頻幅度譜 A (τ,f) 和高精度時(shí)頻相位譜 θ (τ,f) 后, 提取頻譜峰值幅度 A(τ,fτ)和峰值相位 θ (τ,fτ) 分別如圖3(d)和圖3(e)所示,隨機(jī)噪聲的影響被消除, 信號(hào)的瞬時(shí)幅度展現(xiàn)出明顯的指數(shù)衰減趨勢, 瞬時(shí)相位維持在0°水平上.重構(gòu)MRS包絡(luò)信號(hào)由圖3(f)給出, 提取信號(hào)參數(shù)E0= 199.14 nV,= 149.2 ms, 信噪比提高為30.69 dB.

圖3 MSTFT方法消除隨機(jī)噪聲過程示例 (a) 仿真含噪數(shù)據(jù); (b)瞬時(shí)幅度; (c) 瞬時(shí)相位; (d) 峰值幅度; (e) 峰值相位; (f)重構(gòu)信號(hào)Fig.3.Example of the basic procedure for the random noise suppression using MSTFT: (a) Synthetic noisy signal; (b) instantaneous amplitude; (c) instantaneous phase; (d) peak amplitude; (e) peak phase; (f) reconstructed signal.

4 仿真分析

為了驗(yàn)證改進(jìn)的短時(shí)傅里葉方法消除隨機(jī)噪聲, 提取MRS信號(hào)的有效性, 在仿真隨機(jī)噪聲和實(shí)測隨機(jī)噪聲中分別加入模擬的單指數(shù)衰減信號(hào),除了在波形上觀察隨機(jī)噪聲的消噪效果外, 采用Müller-Petke等[15]于2016提出的MRSmatlab軟件中的同步檢測方法提取信號(hào)參數(shù), 并基于信噪比和均方根誤差(root mean square error, RMSE)來量化評(píng)估本文所提方法的消噪效果.均方根誤差定義為

其中K為信號(hào)長度, x (k) 為理想信號(hào), r (k) 為重構(gòu)信號(hào).

4.1 仿真噪聲與仿真信號(hào)的消噪示例

單指數(shù)信號(hào)的仿真參數(shù)E0= 200 nV,=150 ms, Δ f = 0 Hz, φ0= 60°, 信號(hào)持續(xù)時(shí)間為500 ms, 信號(hào)長度為1165.為了使仿真的隨機(jī)噪聲更接近野外實(shí)際環(huán)境中的隨機(jī)噪聲, 根據(jù)(18)式進(jìn)行仿真[28]:

其中, k為離散時(shí)間樣本, N(0,1)(k) 表示均值為0標(biāo)準(zhǔn)差為1 nV的高斯噪聲, Vbac(m) 表示背景噪聲, 考慮到結(jié)構(gòu)噪聲等非特定噪聲的存在, Vuni(m)為加入的均勻噪聲.

圖4給出了不同噪聲條件下MSTFT方法對(duì)隨機(jī)噪聲的消除效果.信號(hào)1中背景噪聲 Vbac(m)是均值為0標(biāo)準(zhǔn)差為50 nV的高斯噪聲, 均勻噪聲 Vuni(m) 是1%的理想信號(hào)值, 加入噪聲后信噪比為0.68 dB.然后應(yīng)用MSTFT算法處理噪聲,圖4(a)為含噪信號(hào)的高精度時(shí)頻幅度譜, 噪聲在時(shí)頻域內(nèi)隨機(jī)分布, 但是對(duì)信號(hào)峰值影響較弱, 峰值出現(xiàn)在0 Hz附近, 幅度在300 nV以內(nèi).沿頻率軸提取各個(gè)時(shí)刻幅度和相位的最大值后重構(gòu)信號(hào), 信號(hào)時(shí)域和頻域的處理結(jié)果分別如圖4(b)和圖4(c)所示, 其中灰色曲線表示含噪信號(hào), 黑色曲線為采用MSTFT方法處理數(shù)據(jù)后重構(gòu)的MRS信號(hào), 紅色曲線為仿真信號(hào), 可以看出, 隨機(jī)噪聲成分被消除, 信號(hào)得以保留.消噪之后提取的信號(hào)參數(shù)和SNR, RMSE如表1首行所列, 參數(shù)提取結(jié)果E0= (200.19 ± 5.01) nV,= (149.2 ± 2.8)ms,Δf = (-0.03 ± 0.01) Hz, 信號(hào)SNR提升為32.67 dB,RMSE為(1.03 ± 0.55) nV.為了測試該方法在較高噪聲水平下的有效性, Vbac(m) 的標(biāo)準(zhǔn)差增大為100 nV, Vuni(m) 增大為理想信號(hào)的3%, 使信號(hào)2的信噪比為-5.20 dB.表1的中間行為MSTFT方法對(duì)信號(hào)2的消噪結(jié)果, 提取參數(shù)和 Δf分別為(202.61 ± 7.90) nV, (152.0 ± 11.8) ms和(0.04 ± 0.03) Hz, SNR提高到24.01 dB, RMSE為(3.79 ± 1.89) nV.最后, 加入干擾程度嚴(yán)重的隨機(jī)噪聲, Vbac(m) 標(biāo)準(zhǔn)差為200 nV, Vuni(m) 為5%理想信號(hào), 信噪比降低至-11.22 dB.從圖4(h)的時(shí)間序列上看, 在隨機(jī)噪聲干擾較大的情況下,MSTFT方法仍然具有良好的消噪性能, 由圖4(g)可以看出, 雖然隨機(jī)噪聲水平增大了, 但是由于其隨機(jī)分布而不能在時(shí)頻域產(chǎn)生聚集性, 因此對(duì)信號(hào)峰值處幅度和相位影響較小, 可以直接提取而不損耗信號(hào)信息.消噪后MRS信號(hào)的參數(shù)提取結(jié)果E0為(204.12 ± 14.07) nV,為(154.4 ± 14.5) ms,Δf 為(0.04 ± 0.06) Hz, 信噪比SNR和均方根誤差RMSE分別為20.81 dB和(5.81 ± 2.42) nV.

表1 3組仿真隨機(jī)噪聲消噪后參數(shù)估計(jì)情況表Table 1.The parameter estimation after de-noising of 3 sets of simulated random noise.

圖4 3組仿真隨機(jī)噪聲消噪結(jié)果圖 (a)低噪聲強(qiáng)度下仿真數(shù)據(jù)高精度時(shí)頻域振幅; (b) 低噪聲強(qiáng)度下消噪結(jié)果時(shí)域圖; (c) 低噪聲強(qiáng)度下消噪結(jié)果頻域圖; (d) 中噪聲強(qiáng)度下仿真數(shù)據(jù)高精度時(shí)頻域振幅; (e) 中噪聲強(qiáng)度下消噪結(jié)果時(shí)域圖; (f) 中噪聲強(qiáng)度下消噪結(jié)果頻域圖; (g) 高噪聲強(qiáng)度下仿真數(shù)據(jù)高精度時(shí)頻域振幅; (h) 高噪聲強(qiáng)度下消噪結(jié)果時(shí)域圖; (i) 高噪聲強(qiáng)度下消噪結(jié)果頻域圖Fig.4.The de-noising results of 3 sets of random noise simulation: (a) High-precision time-frequency domain amplitude of simulated data under low noise intensity; (b) time domain results under low noise intensity; (c) frequency domain results under low noise intensity; (d) high-precision time-frequency domain amplitude of simulated data under moderate noise intensity; (e) time domain results under moderate noise intensity; (f) frequency domain results under moderate noise intensity; (g) high-precision timefrequency domain amplitude of simulated data under high noise intensity; (h) time domain results under high noise intensity; (i) frequency domain results under high noise intensity.

4.2 實(shí)測噪聲與仿真信號(hào)的消噪示例

為了進(jìn)一步驗(yàn)證MSTFT方法的有效性, 本節(jié)進(jìn)行了另一組仿真實(shí)驗(yàn), 由(2)式給出的仿真信號(hào)中E0= 200 nV,= 100 ms, Δ f = 0 Hz, φ0=57°, 信號(hào)持續(xù)時(shí)間為256 ms, 信號(hào)長度為596.加入在兩個(gè)不同地點(diǎn)進(jìn)行MRS實(shí)驗(yàn)采集的兩組實(shí)測噪聲數(shù)據(jù), 由于噪聲記錄中既包含尖峰噪聲和工頻諧波噪聲, 又包含隨機(jī)噪聲, 因此先采用統(tǒng)計(jì)疊加法和自適應(yīng)陷波法消除尖峰噪聲和工頻諧波噪聲[16], 只保留隨機(jī)噪聲, 并對(duì)隨機(jī)噪聲進(jìn)行重采樣,使其與仿真信號(hào)長度相等.然后采用MSTFT方法分別對(duì)兩個(gè)測點(diǎn)數(shù)據(jù)中的一次噪聲記錄進(jìn)行隨機(jī)噪聲壓制, 并將結(jié)果與疊加法處理同一脈沖矩?cái)?shù)據(jù)的結(jié)果進(jìn)行對(duì)比.由于測點(diǎn)2的噪聲水平略高于測點(diǎn)1, 因此疊加法處理隨機(jī)噪聲時(shí), 對(duì)測點(diǎn)1的每組脈沖矩?cái)?shù)據(jù)疊加16次, 測點(diǎn)2的每組脈沖矩?cái)?shù)據(jù)疊加32次.

圖5 給出了傳統(tǒng)疊加法和改進(jìn)短時(shí)傅里葉變換方法壓制隨機(jī)噪聲的對(duì)比結(jié)果, 從第一列圖中可以看出, 疊加法處理數(shù)據(jù)后的結(jié)果中仍殘留隨機(jī)噪聲成分, 而圖5中間列中, MSTFT方法處理單次噪聲記錄, 隨機(jī)噪聲抑制效果明顯, 信號(hào)衰減趨勢與仿真信號(hào)接近, 并且與傳統(tǒng)疊加法相比, 該方法在較低信噪比下仍能提取所需信號(hào).

圖5 疊加法和MSTFT方法消除隨機(jī)噪聲效果對(duì)比圖 (a) 疊加法消除隨機(jī)噪聲時(shí)域圖; (b) MSTFT方法消除隨機(jī)噪聲時(shí)域圖; (c) 疊加法和MSTFT方法消除隨機(jī)噪聲頻域?qū)Ρ葓DFig.5.Comparison of the de-noising results by using stacking and MSTFT methods: (a) Results of random noise elimination by stacking in time domain; (b) results of random noise elimination by MSTFT in time domain; (c) comparison of the de-noising results by using stacking and MSTFT methods in frequency domain.

表2列出了實(shí)測數(shù)據(jù)中選定脈沖矩的單指數(shù)信號(hào)擬合參數(shù)結(jié)果, 并且計(jì)算了每種情況下信號(hào)的SNR和RMSE.可以看出, 疊加法處理隨機(jī)噪聲后, 提取參數(shù)仍具有較大的擬合誤差, 信噪比水平較低, 且每個(gè)信號(hào)T2*參數(shù)的估計(jì)誤差較大.采用MSTFT方法處理單次記錄數(shù)據(jù)提取的信號(hào)參數(shù)與仿真信號(hào)參數(shù)更為接近.實(shí)驗(yàn)結(jié)果表明在信噪比高于-17.21 dB時(shí)均可有效提取信號(hào), 與疊加法消噪結(jié)果相比, 信號(hào)的信噪比SNR最多可提高27.88 dB, 均方根誤差RMSE最多縮小36.44倍,信號(hào)擬合結(jié)果更精確.

表2 疊加法和MSTFT消除隨機(jī)噪聲提取參數(shù)結(jié)果對(duì)比Table 2.Comparison of the parameter estimation after random noise elimination by stacking and MSTFT methods.

5 實(shí)測結(jié)果

本節(jié)分別采用MSTFT方法和疊加法對(duì)同一組實(shí)測MRS數(shù)據(jù)進(jìn)行處理, 對(duì)比兩種方法下該組實(shí)測數(shù)據(jù)的消噪結(jié)果、參數(shù)估計(jì)和MRS反演情況.所有數(shù)據(jù)均使用吉林大學(xué)自主開發(fā)的JLMRS-I系統(tǒng)在吉林省長春市郊區(qū)進(jìn)行采集, 該地的地磁場強(qiáng)度為54720 nT, 拉莫爾頻率為2330 Hz.接收線圈為邊長100 m的正方形, 16個(gè)發(fā)射脈沖矩分布在240-8350 A·ms范圍內(nèi), 每組脈沖矩?cái)?shù)據(jù)采集16次.

圖6所示為不同脈沖矩下原始數(shù)據(jù)及經(jīng)疊加法和MSTFT方法處理后的數(shù)據(jù)隨時(shí)間的變化情況.從圖6(a)中可以看出, 原始數(shù)據(jù)受噪聲干擾嚴(yán)重, 無法辨別信號(hào)的存在.因此采用預(yù)處理方法預(yù)先消除尖峰噪聲和工頻諧波噪聲, 只保留隨機(jī)噪聲, 并對(duì)含噪數(shù)據(jù)進(jìn)行重采樣.然后采用傳統(tǒng)疊加法消除隨機(jī)噪聲, 結(jié)果如圖6(b)所示, 可以看出不同脈沖矩下信號(hào)幅度的變化情況和衰減趨勢, 但未被完全消除的隨機(jī)噪聲仍會(huì)對(duì)信號(hào)產(chǎn)生影響.采用MSTFT方法消除單次記錄的預(yù)處理數(shù)據(jù)中的隨機(jī)噪聲, 結(jié)果如圖6(c)所示, 隨機(jī)噪聲抑制效果更加明顯, 信號(hào)更為突出.

圖6 實(shí)測數(shù)據(jù)處理結(jié)果圖 (a) 野外探測原始數(shù)據(jù); (b) 疊加法消除隨機(jī)噪聲后數(shù)據(jù); (c) MSTFT方法消除隨機(jī)噪聲后數(shù)據(jù)Fig.6.De-noising results of field data: (a) Original field data; (b) random noise elimination by stacking; (c) random noise elimination by MSTFT.

圖7為16個(gè)發(fā)射脈沖矩中任意4個(gè)脈沖矩信號(hào)的消噪結(jié)果.其中黑色曲線為疊加法消噪后的MRS包絡(luò)信號(hào), 紅色曲線為MSTFT方法處理單次數(shù)據(jù)記錄后得到的MRS包絡(luò)曲線.圖7(a)和圖7(c)中可以觀察到指數(shù)衰減趨勢, 但是圖7(e)和圖7(g)中SNR仍然很低.由時(shí)域和頻域消噪結(jié)果可以看出, 與疊加法相比, 經(jīng)MSTFT方法壓制后信號(hào)中隨機(jī)噪聲水平更低.并且如圖7(f)和圖7(h)的內(nèi)置圖所示, 當(dāng)信噪比越低時(shí), MSTFT方法的優(yōu)勢越明顯.另外, 雖然實(shí)際測量環(huán)境中的隨機(jī)噪聲不是純白噪聲數(shù)據(jù), 但MSTFT方法仍然可行.

圖7 疊加法和MSTFT方法處理4個(gè)脈沖矩信號(hào)結(jié)果對(duì)比圖 (a) 時(shí)域?qū)Ρ冉Y(jié)果; (b) 頻域?qū)Ρ冉Y(jié)果Fig.7.Comparison of the de-noising results of 4 pulse moments by using stacking and MSTFT methods: (a) Time domain comparison; (b) frequency domain comparison.

16個(gè)脈沖矩信號(hào)的參數(shù)估計(jì)情況如圖8所示.可以看出, 采用傳統(tǒng)疊加法和MSTFT方法處理數(shù)據(jù)后提取的信號(hào)初始振幅差距很小, 但是MSTFT方法消噪后提取的信號(hào)弛豫時(shí)間比疊加法消噪后提取的弛豫時(shí)間更為準(zhǔn)確, 兩種方法消噪后信號(hào)的頻率偏移在零附近波動(dòng)并且變化較為平滑, 因此可以認(rèn)為采用MSTFT方法抑制隨機(jī)噪聲后MRS信號(hào)的提取結(jié)果更為準(zhǔn)確.

圖8 疊加法和MSTFT處理實(shí)測數(shù)據(jù)提取參數(shù)結(jié)果對(duì)比 (a) 初始振幅; (b) 弛豫時(shí)間; (c) 頻率偏移Fig.8.Parameter estimation after random noise elimination of field data by using stacking and MSTFT methods: (a) Initial amplitude; (b) relaxation time; (c) frequency offset.

對(duì)兩種方法消噪后提取的MRS信號(hào)分別進(jìn)行初始振幅反演, 反演結(jié)果和實(shí)驗(yàn)地點(diǎn)附近的鉆孔結(jié)果對(duì)比如圖9所示.可以看出, 兩組數(shù)據(jù)的反演結(jié)果得到的含水層信息均與鉆孔結(jié)果顯示的地層結(jié)構(gòu)一致.該地點(diǎn)存在兩個(gè)含水層, 其中主含水層位于地下12.5-24 m的細(xì)-中-粗砂孔隙介質(zhì)中, 最大含水量在11%-12%之間, 次含水層位于地下38-40 m的中砂孔隙介質(zhì)中, 平均含水量約為4%.但是相比于疊加法, 由MSTFT方法消噪后的數(shù)據(jù)進(jìn)行反演, 反演結(jié)果得到的含水層信息與實(shí)際更為接近, 且實(shí)驗(yàn)數(shù)據(jù)采集中只需進(jìn)行一次探測, 提高了野外工作效率.

圖9 疊加法和MSTFT處理實(shí)測數(shù)據(jù)提取MRS信號(hào)反演結(jié)果與鉆孔結(jié)果對(duì)比 (a) 疊加法; (b) MSTFT方法Fig.9.Comparison of the inversion of the field data processd by stacking and MSTFT methods respectively with the borehole log:(a) Stacking result; (b) MSTFT result.

6 結(jié) 論

隨機(jī)噪聲由于其不確定性、無法預(yù)測而成為MRS信號(hào)中難以抑制的噪聲種類之一.在強(qiáng)噪聲干擾下有效抑制隨機(jī)噪聲, 準(zhǔn)確提取MRS信號(hào)是磁共振測深法能夠廣泛應(yīng)用的基礎(chǔ).針對(duì)現(xiàn)有隨機(jī)噪聲消除方法存在的探測效率低且只能抑制部分隨機(jī)噪聲的問題, 本文研究了一種改進(jìn)的短時(shí)傅里葉變換方法, 通過對(duì)單次采集的MRS包絡(luò)信號(hào)進(jìn)行Hilbert變換求取解析信號(hào), 然后代替實(shí)值信號(hào)進(jìn)行短時(shí)傅里葉變換后提取時(shí)頻峰值幅度和峰值相位重構(gòu)信號(hào), 實(shí)現(xiàn)MRS信號(hào)中隨機(jī)噪聲的有效抑制.通過理論推導(dǎo)和分析數(shù)值仿真和實(shí)測數(shù)據(jù)處理結(jié)果, 可以得出以下結(jié)論:

1)在使用相同時(shí)間窗函數(shù)的情況下, 與常規(guī)短時(shí)傅里葉變換方法相比, MSTFT方法獲取信號(hào)的瞬時(shí)相位不變, 但是瞬時(shí)幅度更為精確, 且由于時(shí)頻域中峰值幅度和峰值相位受隨機(jī)噪聲影響小,可以用來重構(gòu)信號(hào)抑制隨機(jī)噪聲;

2)與傳統(tǒng)疊加法相比, MSTFT方法可以在較低信噪比下提取所需信號(hào), 且信噪比越低時(shí), MSTFT方法的優(yōu)勢越明顯, 仿真實(shí)驗(yàn)表明, MSTFT方法在信噪比高于-17.21 dB時(shí)可有效提取信號(hào),與疊加法相比, 消噪后信噪比最多可提高27.88 dB,均方根誤差RMSE最多縮小到原來的1/36.44.另外, MSTFT方法處理單次記錄數(shù)據(jù)提取的信號(hào)參數(shù)更為準(zhǔn)確, 尤其能提高參數(shù)估計(jì)的穩(wěn)定性.采用MSTFT方法消除隨機(jī)噪聲可以優(yōu)化消噪效果,提高探測效率, 為MRS方法快速獲取地下水信息提供可行性.