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

1.吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院，長(zhǎng)春 130026 2.地球信息探測(cè)儀器教育部重點(diǎn)實(shí)驗(yàn)室(吉林大學(xué))，長(zhǎng)春 130026 3.國(guó)家地球物理探測(cè)儀器工程技術(shù)研究中心(吉林大學(xué))，長(zhǎng)春 130026

0 引言

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

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

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

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

1 原理與方法

1.1 磁共振信號(hào)及其噪聲特征

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

(1)

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

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

VMRS=VFID+Vrandom+Vspike+Vharmonic。

(2)

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

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

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

(3)

(4)

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

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

(5)

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

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

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

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

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

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

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

(6)

a. 純凈MRS信號(hào)的實(shí)部；b. 純凈MRS信號(hào)的虛部；c. 僅含隨機(jī)噪聲MRS信號(hào)的實(shí)部；d. 僅含隨機(jī)噪聲MRS信號(hào)的虛部。

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

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

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

CNN大多采用池化層對(duì)上一層的輸出特征進(jìn)行下采樣以及對(duì)網(wǎng)絡(luò)引入平移不變性。本文研究?jī)?nèi)容為由輸入含噪時(shí)頻圖像得出純凈MRS信號(hào)的時(shí)頻圖像，由于訓(xùn)練集中信號(hào)的頻率與測(cè)試MRS信號(hào)的頻率較為接近，訓(xùn)練集和驗(yàn)證集的時(shí)頻圖像出現(xiàn)位置相差不大，因此不需要用到池化層的平移不變性?；谝陨蟽牲c(diǎn)原因，本文構(gòu)建了不含有池化層的CNN。

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

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

(7)

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

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

3 結(jié)果

3.1 實(shí)驗(yàn)過(guò)程

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

設(shè)置每次訓(xùn)練的批次大小為512，初始學(xué)習(xí)率設(shè)置為10-3，每迭代4輪之后，學(xué)習(xí)率衰減為之前的1/2，一共迭代24輪。

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

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

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

3.3 噪聲影響

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

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

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

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

3.4 CNN與TFPF對(duì)比

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

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

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

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

(8)

式中，k為時(shí)頻譜矩陣中元素?cái)?shù)目。

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

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

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

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

3.5 野外實(shí)測(cè)數(shù)據(jù)驗(yàn)證

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

4 結(jié)論

1)為了克服隨機(jī)噪聲給磁共振反演解釋帶來(lái)的問(wèn)題，本文設(shè)計(jì)了用于壓制MRS信號(hào)中隨機(jī)噪聲的CNN，并使用仿真數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)做了驗(yàn)證。

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

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