基于模擬退火法的磁共振測(cè)深多源諧波噪聲壓制方法

陳亮，付立恒，蔡凍，李凡，李振宇，魯愷

(1.中國(guó)電建集團(tuán) 江西省電力設(shè)計(jì)院有限公司，江西 南昌 330096； 2.中國(guó)地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院，湖北 武漢 430074)

0 引言

1 MRS方法的基本原理及多源諧波電磁噪聲特征

1.1 MRS方法的基本原理

MRS方法以NMR理論為基礎(chǔ)，即地下水中的氫核磁矩受到一定頻率的交變電流脈沖激發(fā)而發(fā)生偏轉(zhuǎn)，磁矩弛豫過(guò)程中產(chǎn)生核磁共振信號(hào)，位于地面的接收線圈接收到一個(gè)自由感應(yīng)衰減信號(hào)，即NMR信號(hào)，其過(guò)程如圖1所示。圖中：0～t1為儀器激發(fā)前的純?cè)肼暡杉瘯r(shí)間；t1～t2為儀器發(fā)射拉莫爾頻率的激發(fā)脈沖持續(xù)時(shí)間，脈沖包絡(luò)線為矩形；t2～t3為發(fā)射轉(zhuǎn)接收時(shí)間，稱為死區(qū)時(shí)間；t3～t4為接收NMR信號(hào)時(shí)間，NMR信號(hào)的包絡(luò)線按指數(shù)規(guī)律衰減。NMR信號(hào)ENMR為：

(1)

圖1 MRS方法激發(fā)/接收NMR信號(hào)時(shí)序圖Fig.1 Transmitting and receiving NMR signal in time domain of MRS

1.2 干擾噪聲的種類及其特征

1.2.1 噪聲種類

MRS方法接收到的NMR信號(hào)通常非常微弱，一般在10～1 000 nV之間，由于儀器接收靈敏度高易受到外部環(huán)境電磁噪聲的干擾，故接收到的信號(hào)信噪比通常較低，提取NMR信號(hào)的難度較大[17]。按照噪聲成分及來(lái)源，電磁噪聲分為尖峰噪聲、諧波噪聲和隨機(jī)噪聲等，可以將接收到的NMR信號(hào)和電磁噪聲表示為

Eobs=ENMR+Espike+Eharmonic+Erandom,

(2)

式中:Eobs為儀器接收到的信號(hào)數(shù)據(jù),ENMR為磁共振信號(hào),Espike為尖峰噪聲,Eharmonic為諧波噪聲,Erandom為隨機(jī)噪聲。

1.2.2 多源諧波電磁噪聲特征

在經(jīng)典諧波噪聲模型條件下，2013年Larsen等[14]提出假設(shè)諧波噪聲的基頻在信號(hào)采集持續(xù)期間內(nèi)是單一且穩(wěn)定的，由此寫出如下諧波模型式：

(3)

式中：Eharmonic(f0,n)表示采樣時(shí)間n處的諧波模型的振幅，Ak、φk分別為第k次諧波分量的振幅和相位，fs為采樣頻率，f0為諧波噪聲的基頻。

隨著MRS勘查工作逐步向城市靠近，面對(duì)的電磁噪聲情況也愈發(fā)復(fù)雜，在測(cè)量區(qū)域可能存在2個(gè)或2個(gè)以上的諧波干擾源。本文分析存在2個(gè)諧波源情形下的諧波噪聲特征，多源諧波與之相通[16]。

假設(shè)有2個(gè)基頻不同的諧波干擾源，且在測(cè)量時(shí)間內(nèi)不變，可以寫出2個(gè)諧波源的諧波模型式：

(4)

式中：f1和f2分別表示要搜索的2個(gè)不同的基頻，Ak和Bk分別表示每個(gè)基頻的第k次諧波的振幅，φk和φk分別表示每個(gè)基頻的第k次諧波的相位。

含有雙基頻的諧波噪聲頻譜相對(duì)于單一且穩(wěn)定基頻的諧波噪聲頻譜會(huì)出現(xiàn)2種特殊的峰值頻譜特征：“彎曲”和“雙峰”(圖2)，“雙峰”峰值頻譜特征出現(xiàn)的條件較“彎曲”峰值頻譜要更為嚴(yán)苛，要求頻率間隔大且幅值相近，因此“雙峰”峰值頻譜的出現(xiàn)概率要低于“彎曲”峰值頻譜。

2 同步刪除法和模擬退火法原理

2.1 同步刪除法原理

由雙基頻諧波噪聲的特征可知[17]，主要問(wèn)題是諧波的頻率彼此接近但不完全相同，使得其頻率內(nèi)容重疊但不完全一致。為解決雙基頻諧波噪聲的相互干擾問(wèn)題，可以拓展模型去噪的搜索方式，在一個(gè)二維頻率空間內(nèi)同步搜索2個(gè)基頻，求取觀測(cè)信號(hào)與雙基頻諧波模型的最小值：

‖Eobs-Eharmonic(f1,f2,n)‖→min 。

(5)

由于式(4)中的f1、f2、Ak、Bk、φk和φk的確定是一個(gè)非線性優(yōu)化問(wèn)題，可以通過(guò)標(biāo)準(zhǔn)化改寫式(4)中的余弦項(xiàng)：

圖2 雙基頻的諧波噪聲頻譜特征示意Fig.2 Spectrum characteristics of harmonic noise with double fundamental frequencies

(6)

式中變量αk、βk、γk和δk可通過(guò)以下關(guān)系與對(duì)應(yīng)諧波源的振幅和相位相關(guān)聯(lián)：

(7)

(8)

將式(6)中的余弦項(xiàng)展開(kāi)變換，假設(shè)f1和f2已知，觀測(cè)信號(hào)與諧波模型的差值最小值問(wèn)題就變成線性，可以用矩陣表示法寫成：

(9)

式中：[k0,k1,…,kN]表示構(gòu)建諧波模型的階數(shù)，[n0,n1,…,nf]為采樣時(shí)間樣本。同時(shí)，可將式(9)整理成矩陣形式：

Ax=b，

(10)

式中：x=[αk0,βk0,γk0,δk0,…,αkN,βkN,γkN,δkN]T,b=[Eobs(n0),Eobs(n1),…,Eobs(nf-1),Eobs(nf)]T。利用最小二乘法求解該標(biāo)準(zhǔn)線性方程可得：

x=(ATA)-1ATb。

(11)

每次代入1對(duì)假設(shè)基頻值f1和f2，可以得到1組諧波參數(shù)。利用式(5)尋找1組能使得該差值最小的假設(shè)的基頻值作為構(gòu)建諧波模型的基頻值，再利用測(cè)量數(shù)據(jù)Eobs減去構(gòu)建的Eharmonic，即可達(dá)到同步壓制2個(gè)干擾源諧波噪聲的效果。

2.2 模擬退火法原理

模擬退火法的發(fā)展受到金屬物理退火的啟發(fā)[18]：當(dāng)金屬?gòu)募t熱狀態(tài)緩慢冷卻時(shí)，金屬內(nèi)部形成有序的最小能量晶體結(jié)構(gòu)。最初，當(dāng)金屬處于紅熱狀態(tài)時(shí)，原子的運(yùn)動(dòng)完全由隨機(jī)的熱漲落控制，但是，隨著溫度的緩慢降低，原子間的作用力變得越來(lái)越重要，最后，原子變成了一個(gè)代表最小能量結(jié)構(gòu)的晶格[19]。在模擬退火算法中，參數(shù)T代表模擬溫度，誤差E代表模擬能量，T值很大時(shí)算法的行為類似于Monte Carlo搜索，T值很小時(shí)算法的行為更有方向性。在物理退火中，一開(kāi)始是一個(gè)大值的T，然后隨著越來(lái)越多的試驗(yàn)被檢驗(yàn)，T值慢慢降低；最初，大量的模型空間是隨機(jī)采樣的，但是隨著搜索的進(jìn)行，搜索變得越來(lái)越有方向性[20]。

模擬退火法從一個(gè)初始模型m(p)和其對(duì)應(yīng)的誤差E(m(p))開(kāi)始。然后生成一個(gè)在m(p)附近的測(cè)試模型m(*)，計(jì)算其對(duì)應(yīng)的誤差E(m(*))，m(*)可以通過(guò)向m(p)增加一個(gè)從高斯分布中提取的增量Δm。當(dāng)E(m(*))≤E(m(p))時(shí)，m(*)會(huì)替代m(p)成為新的比較值。但是,有時(shí)E(m(*))>E(m(p))，m(*)也可以接受。為了決定后面的情況，這個(gè)測(cè)試參數(shù)t可通過(guò)

(12)

計(jì)算。在區(qū)間[0,1]生成一個(gè)均勻分布的隨機(jī)數(shù)r，如果t>r，則接受m(*)。當(dāng)T值很大時(shí)，測(cè)試參數(shù)t趨近于1，故不管誤差值是多少，m(*)幾乎總是被接受的，這對(duì)應(yīng)著“熱運(yùn)動(dòng)”情況下，模型參數(shù)在空間上以無(wú)向的方式探索；當(dāng)T值很小時(shí)，測(cè)試參數(shù)t趨近于0，m(*)幾乎不被接受，這與定向搜索情況相對(duì)應(yīng)，只有使誤差E減小的測(cè)試模型才被接受。

3 仿真數(shù)據(jù)驗(yàn)證

3.1 網(wǎng)格搜索同步刪除法

前文論述過(guò)，若要式(5)成立需要多次求解線性問(wèn)題，同時(shí)優(yōu)化二維網(wǎng)格搜索，每個(gè)維度對(duì)應(yīng)不同的干擾源諧波基頻。為了加快搜索速度，可以將搜索分2步進(jìn)行：首先，覆蓋相對(duì)寬且粗糙的網(wǎng)格，每次計(jì)算模型減去后的剩余信號(hào)；然后，從在第一步中獲得的最小值開(kāi)始探索更細(xì)的網(wǎng)格，并且通過(guò)第二步獲得的最小剩余信號(hào)產(chǎn)生最佳頻率估計(jì)值。

圖3a是在一個(gè)相對(duì)粗糙的網(wǎng)格內(nèi)搜索2個(gè)基頻值，即圖中2個(gè)相互對(duì)稱的深藍(lán)色色塊;圖3b是選取圖3a中的一個(gè)色塊范圍進(jìn)一步精細(xì)化搜索2個(gè)基頻,可以明顯看出找到的基頻值與仿真設(shè)定的一致。搜索基頻頻率的精度要達(dá)到0.001 Hz需要進(jìn)行600多次的搜索，處理時(shí)間為60.36±1.27 s。圖4顯示，同步刪除法能夠壓制雙基頻諧波噪聲，經(jīng)過(guò)處理的信號(hào)與真實(shí)NMR信號(hào)相似程度較高，呈指數(shù)衰減趨勢(shì)。圖中縱坐標(biāo)變量A表示振幅。

圖3 頻率空間諧波基頻網(wǎng)格搜索Fig.3 Harmonic fundamental frequency grid search in frequency space

圖4 網(wǎng)格搜索同步刪除法處理雙基頻諧波時(shí)域Fig.4 Time domain diagram of double fundamental frequency harmonics processed by grid search simultaneous removal method

3.2 模擬退火同步刪除法

為了加快基頻的搜索和避免陷入局部極值點(diǎn)，利用模擬退火法來(lái)搜索干擾源諧波的基頻值。仿真條件同上，迭代次數(shù)200次，搜索路徑如圖5所示。圖5中黑色圓圈代表初始模型(50.000；50.000)，白色圓圈代表每次迭代被接受的模型，綠色圓圈代表迭代結(jié)束最終模型(50.031 7；49.984 6)，紅色線段代表基頻搜索路徑?？梢钥闯觯看伪唤邮艿哪Ｐ投荚谥鸩奖平鎸?shí)模型，雖然模擬退火法最終能找到最小值，但在有限的迭代次數(shù)時(shí)搜索路徑會(huì)存在一定的偏差。為了更好地評(píng)價(jià)去噪效果，特意引入均方根誤差(root mean square error，RMSE)這個(gè)指標(biāo)參數(shù)[21]：

(13)

式中:N為儀器采集到的信號(hào)長(zhǎng)度，Edenoised(i)是去噪后的第i個(gè)數(shù)據(jù)，ENMR(i)是真實(shí)NMR信號(hào)的第i個(gè)數(shù)據(jù)。

圖5 模擬退火同步刪除法的基頻搜索路徑Fig.5 Fundamental frequency search path of simulated annealing simulated annealing

復(fù)雜的電磁噪聲可能存在3個(gè)或更多基頻的多源諧波噪聲。為了驗(yàn)證模擬退火同步刪除法能夠適用于多源諧波電磁噪聲的情況，選取了3個(gè)干擾源進(jìn)行仿真，諧波基頻分別為49.950、50.010、50.050，迭代次數(shù)300次，初始基頻模型(50.000；50.000；50.000)，其他仿真條件同上，圖6為仿真結(jié)果。從圖6a可以看出模擬退火同步刪除法在3個(gè)基頻諧波的情況下能夠較好地壓制多源諧波噪聲，去噪后的藍(lán)色曲線與真實(shí)NMR信號(hào)的紅色曲線基本吻合。圖6b中黑色圓圈顯示了迭代初始模型的頻率空間位置，藍(lán)色圓圈是300次迭代結(jié)束最終模型的頻率空間位置(50.009 2；49.951 2；50.049 9)，紅色圓圈是3個(gè)真實(shí)諧波基頻不同排列順序的頻率空間位置，灰色曲線代表基頻搜索路徑；經(jīng)過(guò)一定迭代次數(shù)的搜索，最終結(jié)果會(huì)歸于6個(gè)紅色圓圈中任意一個(gè)的附近。

表1 不同迭代次數(shù)的模擬退火同步刪除法去噪效果

圖6 模擬退火同步刪除法壓制多源諧波噪聲Fig.6 Simulated annealing simultaneous removal method suppressing Multi- source harmonic noise

4 實(shí)測(cè)數(shù)據(jù)檢驗(yàn)

為了驗(yàn)證本文提出的噪聲壓制方法在實(shí)測(cè)數(shù)據(jù)中的應(yīng)用效果，選擇天門市多寶鎮(zhèn)漢江邊的一處場(chǎng)地進(jìn)行實(shí)驗(yàn)。該場(chǎng)地地下含水層已由場(chǎng)地內(nèi)的鉆孔所揭露，場(chǎng)地附近有多個(gè)人居村莊，符合多源諧波的噪聲特點(diǎn)，實(shí)測(cè)數(shù)據(jù)可以很好地應(yīng)用于去噪處理的實(shí)驗(yàn)工作中。數(shù)據(jù)采集使用法國(guó)IRIS公司的NUMISpoly儀器。本次測(cè)量采集技術(shù)參數(shù)：方形線圈邊長(zhǎng)100 m；拉莫爾頻率fL=2 135.2 Hz；激發(fā)脈沖矩范圍為0.12～10.866 A·s，分為16個(gè)；采樣率19 200 Hz；采樣時(shí)間1 s；死區(qū)時(shí)間0.04 s。

通過(guò)分析野外實(shí)際工作環(huán)境和噪聲，發(fā)現(xiàn)2個(gè)較強(qiáng)的諧波干擾源。根據(jù)數(shù)據(jù)處理流程對(duì)噪聲數(shù)據(jù)進(jìn)行了尖峰噪聲的處理后，采用100次迭代的模擬退火同步刪除法來(lái)壓制雙諧波噪聲(圖7)，再將去噪后的信號(hào)進(jìn)行疊加、希爾伯特變換和非線性擬合，結(jié)果見(jiàn)圖8。

圖7是利用模擬退火同步刪除法來(lái)壓制實(shí)測(cè)雙諧波噪聲。圖7a中經(jīng)過(guò)模擬退火同步刪除法壓制后的去噪曲線(紅色曲線)有明顯的指數(shù)衰減趨勢(shì)，僅剩余隨機(jī)噪聲的干擾；圖7b中相對(duì)于實(shí)測(cè)含噪信號(hào)(灰色曲線)，去噪曲線在50 Hz整數(shù)附近的峰值被較好壓制，拉莫爾頻率的信號(hào)保存較為完整。

圖8中灰色曲線是去噪后的信號(hào)經(jīng)過(guò)疊加處理和希爾伯特變換得到的包絡(luò)曲線，紅色曲線是灰色曲線經(jīng)過(guò)非線性擬合得到的擬合曲線。包絡(luò)曲線和擬合曲線重合性較好，表明模擬退火同步刪除法能夠較好地壓制多諧波噪聲，也能看出殘余諧波噪聲對(duì)包絡(luò)曲線的影響在可接受的范圍，且對(duì)擬合曲線的參數(shù)提取影響較小。

圖7 模擬退火同步刪除法在實(shí)測(cè)數(shù)據(jù)處理中的效果Fig.7 Effect of simulated annealing simultaneous removal method in measured data processing

圖8 去噪后各脈沖矩的包絡(luò)線及其擬合曲線Fig.8 Envelope and fitting curve of each pulse moment after denoising

圖9 天門河灘MRS方法某測(cè)深點(diǎn)實(shí)測(cè)數(shù)據(jù)及其反演結(jié)果與鉆探結(jié)果對(duì)比Fig.9 Measured data, inversion results and drilling results of a sounding point by MRS method in River Beach, Tianmen

5 結(jié)論

研究與實(shí)踐表明，對(duì)于單一且規(guī)律性的諧波電磁噪聲，利用模型去噪的方法能夠有效地壓制，而變化、多源的諧波噪聲的壓制是磁共振測(cè)深方法信號(hào)處理中的難點(diǎn)和研究熱點(diǎn)，壓制效果的好壞直接影響后續(xù)反演解釋工作的準(zhǔn)確性。根據(jù)多源諧波電磁噪聲的特點(diǎn)出發(fā)，驗(yàn)證了網(wǎng)格搜索同步刪除法的可行性，為了節(jié)約計(jì)算時(shí)間成本，提出了模擬退火同步刪除法這一新的磁共振測(cè)深方法信號(hào)處理手段，通過(guò)仿真數(shù)據(jù)驗(yàn)證和實(shí)測(cè)數(shù)據(jù)檢驗(yàn)，本文論述的信號(hào)處理手段取得了較好的噪聲壓制效果。

1)針對(duì)磁共振測(cè)深信號(hào)中模型去噪方法無(wú)法有效壓制多源諧波噪聲的問(wèn)題，利用網(wǎng)格搜索同步刪除法和模擬退火同步刪除法均能有效找到較為準(zhǔn)確的基頻值，達(dá)到壓制多源諧波的目的。

2)通過(guò)模擬仿真實(shí)驗(yàn)，證明網(wǎng)格搜索同步刪除法和模擬退火同步刪除法均能壓制多諧波噪聲。在雙基頻諧波情況下比較不同迭代次數(shù)下的壓制效果，模擬退火同步刪除法較網(wǎng)格搜索同步刪除法的計(jì)算效率提高了2.35倍，大大降低了時(shí)間成本，且在多源諧波噪聲壓制方面也能有較好的效果，為開(kāi)展更加復(fù)雜的諧波情況的壓制提供了可能。

3)通過(guò)野外MRS方法的實(shí)測(cè)資料證明了模擬退火同步刪除法壓制多源諧波噪聲的有效性，反演得到的含水量與鉆探結(jié)果基本相符，極大地拓展了磁共振測(cè)深方法在人居區(qū)等電磁情況復(fù)雜地區(qū)作業(yè)的有效性，也為多通道磁共振多維探測(cè)提供新的數(shù)據(jù)處理思路。

仿真和實(shí)測(cè)結(jié)果表明，模擬退火同步刪除法能夠有效壓制磁共振測(cè)深方法中的多源諧波噪聲，為后期反演資料的解釋提供了技術(shù)保證和選擇方法。本文主要討論單通道情況下多源諧波噪聲的壓制，未來(lái)將繼續(xù)研究二維三維多通道變化諧波噪聲的壓制方法。