李 邦，蔣川東,2,3，王 遠,2,3，田寶鳳,2，段清明,2,3，尚新磊,2,3

1.吉林大學儀器科學與電氣工程學院，長春 130026 2.地球信息探測儀器教育部重點實驗室(吉林大學)，長春 130026 3.國家地球物理探測儀器工程技術(shù)研究中心(吉林大學)，長春 130026

0 引言

地面磁共振探測(magnetic resonance sounding，MRS)是一種非侵入式直接探測地下水的技術(shù)。和其他地球物理探測方法相比，MRS技術(shù)具有直接高效、信息量豐富、定量準確及解釋唯一的優(yōu)點[1]。近年來，MRS技術(shù)還被應用到水污染監(jiān)測、水文環(huán)境評價[2]和地質(zhì)災害預警等領(lǐng)域[3]。

MRS信號微弱，僅為納伏級別，極易受到周圍環(huán)境中電磁干擾的影響，導致采集到的MRS數(shù)據(jù)質(zhì)量嚴重依賴于測量地點的電磁噪聲水平。尤其在城市和村鎮(zhèn)這些人類活動較多的區(qū)域進行磁共振探測時，MRS信號被噪聲完全淹沒，信噪比極低，這使得MRS數(shù)據(jù)反演解釋的準確度嚴重下降[4]。同時這也是制約目前MRS技術(shù)大規(guī)模應用的主要瓶頸。因此，研究磁共振信號的噪聲壓制具有非常重要的意義。

根據(jù)噪聲來源和信號特征的不同，地面磁共振數(shù)據(jù)中的噪聲可被分為隨機噪聲、工頻諧波噪聲和尖峰噪聲3種。去除工頻諧波噪聲的方法有自適應濾波[5]和工頻諧波建模[6]等，去除尖峰噪聲可以用壓縮小波變換方法[7]，通過上述方法可以較好地去除這兩種噪聲。隨機噪聲的特征為高斯白噪聲，一般通過疊加多次測量數(shù)據(jù)取平均值的方式減小[8]。2018年，王琦等[9]提出基于稀疏表示的隨機噪聲背景下多弛豫MRS信號的提取方法。林婷婷等[10]使用時頻峰值濾波(time-frequency peak filtering，TFPF)方法，對隨機噪聲起到了良好的壓制效果。但上述兩種方法處理過程較為繁瑣，計算耗費時間長，對于野外實驗取得的大量數(shù)據(jù)難以實現(xiàn)實時處理，且仍不能完全達到滿意的噪聲壓制效果。

針對MRS信號隨機噪聲壓制中存在的上述問題，本文提出使用卷積神經(jīng)網(wǎng)絡(convolutional neural networks，CNN)[11]對磁共振數(shù)據(jù)中的隨機噪聲進行壓制。本文首先介紹地面磁共振數(shù)據(jù)信號和噪聲的特征，以及基于CNN的MRS數(shù)據(jù)噪聲壓制的過程。然后，利用CNN對仿真MRS數(shù)據(jù)進行隨機噪聲壓制，并與TFPF方法處理的結(jié)果進行對比，證明CNN的優(yōu)勢。最后，通過野外實測數(shù)據(jù)的MRS噪聲壓制結(jié)果，驗證了基于CNN的MRS數(shù)據(jù)隨機噪聲壓制方法的有效性和實用性。

1 原理與方法

1.1 磁共振信號及其噪聲特征

在發(fā)射線圈中通入Lamor頻率電流，向地下發(fā)射電磁場，地下水中的氫質(zhì)子被電磁場激發(fā)，形成宏觀磁矩，這一宏觀磁矩在地磁場B0中以B0方向為軸進行旋進運動。自由感應衰減(free induction decay, FID)是核磁共振與磁共振成像中最簡單的信號形式。受激發(fā)的氫核對磁振頻譜儀或磁振造影掃描器的射頻線圈造成感應電流而產(chǎn)生信號，并因為發(fā)生弛豫而使信號強度逐漸衰減至0，這種逐漸衰減的信號即稱為FID。當激發(fā)場停止后，氫質(zhì)子旋進產(chǎn)生弛豫現(xiàn)象；地面接收線圈記錄到的宏觀磁矩進動產(chǎn)生的電磁信號就是FID信號VFID，可表示為

(1)

由于FID信號極為微弱，單位僅為納伏級別，要采集這樣的微弱信號就必須使用高靈敏度的接收器。因此，最終采集到的信號中包含了大量噪聲。

采集到的MRS數(shù)據(jù)可表示為

VMRS=VFID+Vrandom+Vspike+Vharmonic。

(2)

式中：Vrandom為隨機噪聲；Vspike為尖峰噪聲；Vharmonic為工頻諧波噪聲。由于尖峰噪聲和工頻諧波噪聲已經(jīng)有了多種成熟有效的去除方式，因此本文僅對MRS信號中隨機噪聲的壓制展開研究。

1.2 MRS數(shù)據(jù)預處理

對采集到的數(shù)據(jù)進行尖峰壓制和工頻諧波建模消噪后，得到僅包含MRS信號和隨機噪聲的數(shù)據(jù)，對其進行希爾伯特變換，得到包絡信號(以單指數(shù)信號為例)的實部Vr和虛部Vi分別為：

(3)

(4)

式中：df為頻率偏量；εr與εi分別為噪聲的實部和虛部分量。

CNN善于提取二維信息中的特征[11-12]。為了使用MRS數(shù)據(jù)訓練CNN，首先利用式(5)對MRS信號E=Er+iEi(Er和Ei分別為MRS信號的實部和虛部)進行短時傅里葉變換 (short-time Fourier transform，STFT)：

(5)

式中：τ和f分別為時頻譜的時間點和頻率點；ω(t-τ)為窗函數(shù)；X為神經(jīng)網(wǎng)絡的輸入時頻譜。

圖1展示了純凈MRS信號和含噪MRS包絡數(shù)據(jù)經(jīng)過STFT后的時頻譜[12-13]?？梢钥闯?，STFT的結(jié)果可以同時反映信號在時域和頻域上的分布。含有更多信息的二維時頻信號可以使神經(jīng)網(wǎng)絡學習到更多MRS信號特征，提高噪聲壓制結(jié)果的準確度。

2 基于卷積神經(jīng)網(wǎng)絡的噪聲壓制方法

CNN是一種多層的監(jiān)督學習神經(jīng)網(wǎng)絡，運行的機制是輸入層的數(shù)據(jù)在多個隱藏層中依次進行特征提取，最后在輸出層輸出，若輸出與理想結(jié)果偏差超過可接受的范圍，則調(diào)整隱藏層神經(jīng)元的權(quán)重和偏置[13]。這樣根據(jù)輸入和輸出結(jié)果之間的關(guān)系，訓練得到一個最優(yōu)的模型的方法，稱為監(jiān)督學習。本文利用如圖2所示CNN結(jié)構(gòu)。

網(wǎng)絡主要包含輸入層、卷積層、批歸一化(batch normalization，BN)層、激活函數(shù)層和輸出層。

卷積層通過將卷積核與輸入數(shù)據(jù)作卷積運算，提取輸入數(shù)據(jù)不同位置的局部特征。而通過設置平移步長，可以令卷積核分別提取數(shù)據(jù)不同位置的特征，最終得到輸入數(shù)據(jù)的特征圖(feature map)。卷積操作的公式為

C(i,j)=(X*W)(i,j)=

(6)

a. 純凈MRS信號的實部；b. 純凈MRS信號的虛部；c. 僅含隨機噪聲MRS信號的實部；d. 僅含隨機噪聲MRS信號的虛部。

含噪時頻譜和結(jié)果時頻譜圖片前面的是實部，后面的是虛部(示意圖)。

式中：i、j、m、n分別為位置坐標；W為卷積核；*為卷積運算；C為卷積結(jié)果。

為了避免神經(jīng)網(wǎng)絡的內(nèi)協(xié)變量偏移，每層使用批歸一化方法[14]。批歸一化層在神經(jīng)網(wǎng)絡層的中間進行預處理，即上一層的輸入經(jīng)歸一化處理后再進入網(wǎng)絡的下一層，這樣可有效地防止訓練時發(fā)生梯度爆炸，并加速網(wǎng)絡訓練。每層卷積層均使用線性修正單元(rectified linear unit，ReLU)函數(shù)[15-16]作為激活函數(shù)。

CNN大多采用池化層對上一層的輸出特征進行下采樣以及對網(wǎng)絡引入平移不變性。本文研究內(nèi)容為由輸入含噪時頻圖像得出純凈MRS信號的時頻圖像，由于訓練集中信號的頻率與測試MRS信號的頻率較為接近，訓練集和驗證集的時頻圖像出現(xiàn)位置相差不大，因此不需要用到池化層的平移不變性。基于以上兩點原因，本文構(gòu)建了不含有池化層的CNN。

隱藏層中包含27個卷積層，每個卷積層后都設置一個批量歸一化層與激活函數(shù)層，共81層。不同卷積層所對應的神經(jīng)元節(jié)點數(shù)在300～1 000不等。第一個卷積層使用了大小為9×8的卷積核，后續(xù)卷積層使用的卷積核大小為9×1與5×1。就如前文所述，卷積核大小為9×1的卷積層代替了池化層的作用。在隱藏層的最后，使用全連接(full connection，F(xiàn)C)層獲得神經(jīng)網(wǎng)絡的輸出。

本文將純凈FID數(shù)據(jù)和去噪后數(shù)據(jù)的時頻譜之間的均方誤差(mean squared error，EMS)作為損失函數(shù)來訓練網(wǎng)絡參數(shù)：

(7)

最后，我們采用Adam算子[17]最小化損失函數(shù)對網(wǎng)絡進行優(yōu)化。

第二天一早，一看表，七點多啦！完蛋了，看來又得挨批了。我趕緊起床洗漱，一抬頭，哇！鏡子里那只“熊貓”是誰呀！都是蚊子惹的禍！唉，蚊子呀蚊子，你犯下的滔天大罪我會永遠記住的！

3 結(jié)果

3.1 實驗過程

對訓練集數(shù)據(jù)的實部和虛部分別進行STFT，得到其時頻譜。實驗中STFT的窗寬設置為256，重疊率為0.75。為了降低神經(jīng)網(wǎng)絡所需的運算資源，在進行訓練時，只取-700～700 Hz的時頻信息作為神經(jīng)網(wǎng)絡的輸入和輸出。每一次訓練分別是以含噪信號的時頻矩陣作為網(wǎng)絡輸入，對應的純凈MRS信號的時頻譜作為輸出。

設置每次訓練的批次大小為512，初始學習率設置為10-3，每迭代4輪之后，學習率衰減為之前的1/2，一共迭代24輪。

本文訓練CNN所用的計算機配置為：CPU i5 8500，RAM 32 Gb。該計算機使用給定的訓練集對CNN進行訓練，總用時538 min。神經(jīng)網(wǎng)絡模型一次訓練完成后，后續(xù)每次對MRS信號做噪聲壓制時只需把信號做STFT后，得到的時頻矩陣與神經(jīng)網(wǎng)絡的參數(shù)矩陣做線性運算，這個過程僅需0.32 s，而TFPF方法耗時為6.89 s，相較之下極大節(jié)省了數(shù)據(jù)處理所需的時間。

3.2 基于CNN的MRS隨機噪聲壓制結(jié)果

使用訓練完成的CNN模型對隨機挑選的一組驗證MRS信號的噪聲壓制后，數(shù)據(jù)的信噪比為5 dB。從圖3a、b可以看出，經(jīng)噪聲壓制之后的MRS信號與純凈MRS信號包絡能夠較好地吻合；由圖3c、d,3e、f與圖3g、h可以看出，經(jīng)CNN噪聲壓制后的時頻譜與純凈信號的時頻譜基本相同，證明了CNN方法可以對MRS信號中的隨機噪聲起到良好的壓制作用且不存在MRS信號的能量損失。

3.3 噪聲影響

a. MRS數(shù)據(jù)實部的時域圖；b. MRS數(shù)據(jù)虛部的時域圖；c. 含噪數(shù)據(jù)實部的時頻譜；d. 含噪數(shù)據(jù)虛部的時頻譜；e. 純凈信號實部的時頻譜；f. 純凈信號虛部的時頻譜；g. 噪聲壓制后實部的時頻譜；h. 噪聲壓制后虛部的時頻譜。

a. SNR=5 dB,實部；b. SNR=5 dB，虛部；c. SNR=0 dB，實部 ；d. SNR=0 dB，虛部；e. SNR=-5 dB，實部；f. SNR=-5 dB，虛部；g. SNR=-10 dB，實部；h. SNR=-10 dB，虛部。

由圖4可以看出，當信噪比大于等于-5 dB時，經(jīng)CNN噪聲壓制后，實部和虛部的包絡曲線與純凈FID信號包絡(圖中紅色曲線)基本重合。當信噪比為-10 dB時，消噪后包絡曲線略有扭曲。

表1 CNN與TFPF單指數(shù)信號噪聲壓制結(jié)果參數(shù)擬合效果對比

3.4 CNN與TFPF對比

圖5為這1組多指數(shù)MRS信號加入不同水平噪聲后分別使用CNN和TFPF兩種方法進行噪聲壓制的結(jié)果。

由圖5可以得出，對于多指數(shù)MRS信號，當信噪比大于0 dB時，CNN和TFPF兩種方法均可以對噪聲實現(xiàn)有效的抑制。當信噪比小于-5 dB時，TFPF方法得到的包絡信號產(chǎn)生嚴重的失真，此時CNN方法的效果顯著優(yōu)于TFPF方法。

針對不同信噪比，分別使用CNN和TFPF對不同信噪比情況下的100組仿真MRS數(shù)據(jù)進行噪聲壓制。單指數(shù)信號處理結(jié)果的信噪比和均方根誤差(root mean squard error，ERMS)見表2，多指數(shù)信號處理結(jié)果的信噪比和均方根誤差見表3。其中均方根誤差采用式(8)計算：

a. SNR=5 dB，實部；b. SNR=5 dB，虛部；c. SNR=0 dB，實部；d. SNR=0 dB，虛部；e. SNR=-5 dB，實部；f. SNR=-5 dB，虛部；g. SNR=-10 dB，實部；h. SNR=-10 dB，虛部。

(8)

式中，k為時頻譜矩陣中元素數(shù)目。

表2 CNN與TFPF單指數(shù)信號噪聲壓制結(jié)果信噪比和均方根誤差對比

表2和表3列出的數(shù)據(jù)表明：經(jīng)CNN方法進行噪聲壓制后的信號擬合誤差相比TFPF方法的結(jié)果擬合誤差更低；在進行實驗對比的信噪比下，CNN的擬合誤差明顯低于TFPF。

另外可以看出，基于CNN和TFPF的兩種MRS隨機噪聲壓制方法，在信噪比高于0 dB時，均可有效地壓制噪聲，其ERMS均小于20 nV。在信噪比低于0 dB時，CNN方法的ERMS仍小于20 nV，而TFPF方法的ERMS已經(jīng)高于20 nV。整體來看，CNN方法噪聲壓制結(jié)果的信噪比提升量在各個場景下都比TFPF方法高出8～11 dB，而CNN處理結(jié)果的ERMS值比TFPF方法降低了59%～82%。由此可以得出：CNN進行噪聲壓制后的MRS仿真信號信噪比獲得了提升，且信噪比提升效果明顯優(yōu)于TFPF方法處理的效果。以ERMS作為評估指標時，CNN處理結(jié)果擁有更低的ERMS值，效果同樣優(yōu)于TFPF方法。

表3 CNN與TFPF多指數(shù)信號噪聲壓制結(jié)果信噪比和均方根誤差對比

3.5 野外實測數(shù)據(jù)驗證

為了驗證本文設計的CNN對實際采集MRS信號進行隨機噪聲壓制的有效性，在吉林省長春市燒鍋鎮(zhèn)進行了MRS探測實驗，探測點附近均為農(nóng)田和樹林，地勢平坦，且臨近水庫，具有豐富的地下水資源。野外實驗數(shù)據(jù)使用JLMRS-IV型磁共振地下水探測儀獲得，測量地點的地磁場強度為54 908 nT，Larmor頻率為2 338 Hz，地磁傾角為60°，發(fā)射/接收一體線圈的尺寸為100 m×100 m，使用不同發(fā)射脈沖矩測得16組MRS信號。以實測磁共振信號作為研究對象，并與TFPF方法作比較。圖6展示其中一組測量信號的噪聲壓制結(jié)果。

4 結(jié)論

1)為了克服隨機噪聲給磁共振反演解釋帶來的問題，本文設計了用于壓制MRS信號中隨機噪聲的CNN，并使用仿真數(shù)據(jù)和實測數(shù)據(jù)做了驗證。

2)使用CNN對單指數(shù)/多指數(shù)MRS信號進行噪聲壓制時均能取得良好的效果。與TFPF方法進行對比，CNN的噪聲壓制效果更優(yōu)秀。

3)通過野外實測數(shù)據(jù)驗證，得到了和仿真實驗相同的結(jié)論，證明了此方法在實際工程應用中的價值。