林婷婷 李玥 高興 萬玲?

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

2) (吉林大學, 地球信息探測儀器教育部重點實驗室, 長春 130026)

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

1 引 言

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

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

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

2 磁共振信號中的隨機噪聲

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

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

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

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

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

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

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

3 基于MSTFT的隨機噪聲抑制方法

3.1 短時傅里葉變換

實值信號 x (t) 的短時傅里葉變換形式為

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

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

3.2 改進的短時傅里葉變換

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

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

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

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

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

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

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

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

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

時, 能得到

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

3.3 信號實例

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

圖2 MSTFT與常規(guī)短時傅里葉變換之間的區(qū)別 (a) 時頻幅度譜 A s(τ,f) ; (b) 時頻相位譜 θ s(τ,f) ; (c) 高精度時頻幅度譜A(τ,f) ; (d) 高精度時頻相位譜 θ (τ,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)為含噪信號的時域波形, 其瞬時幅度(圖3(b))和瞬時相位(圖3(c))受隨機噪聲影響嚴重, 無法準確提取原信號信息.采用前文所述的MSTFT方法獲取高精度時頻幅度譜 A (τ,f) 和高精度時頻相位譜 θ (τ,f) 后, 提取頻譜峰值幅度 A(τ,fτ)和峰值相位 θ (τ,fτ) 分別如圖3(d)和圖3(e)所示,隨機噪聲的影響被消除, 信號的瞬時幅度展現(xiàn)出明顯的指數(shù)衰減趨勢, 瞬時相位維持在0°水平上.重構MRS包絡信號由圖3(f)給出, 提取信號參數(shù)E0= 199.14 nV,= 149.2 ms, 信噪比提高為30.69 dB.

圖3 MSTFT方法消除隨機噪聲過程示例 (a) 仿真含噪數(shù)據(jù); (b)瞬時幅度; (c) 瞬時相位; (d) 峰值幅度; (e) 峰值相位; (f)重構信號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 仿真分析

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

其中K為信號長度, x (k) 為理想信號, r (k) 為重構信號.

4.1 仿真噪聲與仿真信號的消噪示例

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

其中, k為離散時間樣本, N(0,1)(k) 表示均值為0標準差為1 nV的高斯噪聲, Vbac(m) 表示背景噪聲, 考慮到結構噪聲等非特定噪聲的存在, Vuni(m)為加入的均勻噪聲.

圖4給出了不同噪聲條件下MSTFT方法對隨機噪聲的消除效果.信號1中背景噪聲 Vbac(m)是均值為0標準差為50 nV的高斯噪聲, 均勻噪聲 Vuni(m) 是1%的理想信號值, 加入噪聲后信噪比為0.68 dB.然后應用MSTFT算法處理噪聲,圖4(a)為含噪信號的高精度時頻幅度譜, 噪聲在時頻域內(nèi)隨機分布, 但是對信號峰值影響較弱, 峰值出現(xiàn)在0 Hz附近, 幅度在300 nV以內(nèi).沿頻率軸提取各個時刻幅度和相位的最大值后重構信號, 信號時域和頻域的處理結果分別如圖4(b)和圖4(c)所示, 其中灰色曲線表示含噪信號, 黑色曲線為采用MSTFT方法處理數(shù)據(jù)后重構的MRS信號, 紅色曲線為仿真信號, 可以看出, 隨機噪聲成分被消除, 信號得以保留.消噪之后提取的信號參數(shù)和SNR, RMSE如表1首行所列, 參數(shù)提取結果E0= (200.19 ± 5.01) nV,= (149.2 ± 2.8)ms,Δf = (-0.03 ± 0.01) Hz, 信號SNR提升為32.67 dB,RMSE為(1.03 ± 0.55) nV.為了測試該方法在較高噪聲水平下的有效性, Vbac(m) 的標準差增大為100 nV, Vuni(m) 增大為理想信號的3%, 使信號2的信噪比為-5.20 dB.表1的中間行為MSTFT方法對信號2的消噪結果, 提取參數(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.最后, 加入干擾程度嚴重的隨機噪聲, Vbac(m) 標準差為200 nV, Vuni(m) 為5%理想信號, 信噪比降低至-11.22 dB.從圖4(h)的時間序列上看, 在隨機噪聲干擾較大的情況下,MSTFT方法仍然具有良好的消噪性能, 由圖4(g)可以看出, 雖然隨機噪聲水平增大了, 但是由于其隨機分布而不能在時頻域產(chǎn)生聚集性, 因此對信號峰值處幅度和相位影響較小, 可以直接提取而不損耗信號信息.消噪后MRS信號的參數(shù)提取結果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組仿真隨機噪聲消噪后參數(shù)估計情況表Table 1.The parameter estimation after de-noising of 3 sets of simulated random noise.

圖4 3組仿真隨機噪聲消噪結果圖 (a)低噪聲強度下仿真數(shù)據(jù)高精度時頻域振幅; (b) 低噪聲強度下消噪結果時域圖; (c) 低噪聲強度下消噪結果頻域圖; (d) 中噪聲強度下仿真數(shù)據(jù)高精度時頻域振幅; (e) 中噪聲強度下消噪結果時域圖; (f) 中噪聲強度下消噪結果頻域圖; (g) 高噪聲強度下仿真數(shù)據(jù)高精度時頻域振幅; (h) 高噪聲強度下消噪結果時域圖; (i) 高噪聲強度下消噪結果頻域圖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 實測噪聲與仿真信號的消噪示例

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

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

圖5 疊加法和MSTFT方法消除隨機噪聲效果對比圖 (a) 疊加法消除隨機噪聲時域圖; (b) MSTFT方法消除隨機噪聲時域圖; (c) 疊加法和MSTFT方法消除隨機噪聲頻域對比圖Fig.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ù)據(jù)中選定脈沖矩的單指數(shù)信號擬合參數(shù)結果, 并且計算了每種情況下信號的SNR和RMSE.可以看出, 疊加法處理隨機噪聲后, 提取參數(shù)仍具有較大的擬合誤差, 信噪比水平較低, 且每個信號T2*參數(shù)的估計誤差較大.采用MSTFT方法處理單次記錄數(shù)據(jù)提取的信號參數(shù)與仿真信號參數(shù)更為接近.實驗結果表明在信噪比高于-17.21 dB時均可有效提取信號, 與疊加法消噪結果相比, 信號的信噪比SNR最多可提高27.88 dB, 均方根誤差RMSE最多縮小36.44倍,信號擬合結果更精確.

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

5 實測結果

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

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

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

圖7 疊加法和MSTFT方法處理4個脈沖矩信號結果對比圖 (a) 時域對比結果; (b) 頻域對比結果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個脈沖矩信號的參數(shù)估計情況如圖8所示.可以看出, 采用傳統(tǒng)疊加法和MSTFT方法處理數(shù)據(jù)后提取的信號初始振幅差距很小, 但是MSTFT方法消噪后提取的信號弛豫時間比疊加法消噪后提取的弛豫時間更為準確, 兩種方法消噪后信號的頻率偏移在零附近波動并且變化較為平滑, 因此可以認為采用MSTFT方法抑制隨機噪聲后MRS信號的提取結果更為準確.

圖8 疊加法和MSTFT處理實測數(shù)據(jù)提取參數(shù)結果對比 (a) 初始振幅; (b) 弛豫時間; (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.

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

圖9 疊加法和MSTFT處理實測數(shù)據(jù)提取MRS信號反演結果與鉆孔結果對比 (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 結 論

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

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

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