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

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

(1.中國電建集團 江西省電力設計院有限公司，江西 南昌 330096； 2.中國地質大學(武漢) 地球物理與空間信息學院，湖北 武漢 430074)

0 引言

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

1.1 MRS方法的基本原理

MRS方法以NMR理論為基礎，即地下水中的氫核磁矩受到一定頻率的交變電流脈沖激發而發生偏轉，磁矩弛豫過程中產生核磁共振信號，位于地面的接收線圈接收到一個自由感應衰減信號，即NMR信號，其過程如圖1所示。圖中：0～t1為儀器激發前的純噪聲采集時間；t1～t2為儀器發射拉莫爾頻率的激發脈沖持續時間，脈沖包絡線為矩形；t2～t3為發射轉接收時間，稱為死區時間；t3～t4為接收NMR信號時間，NMR信號的包絡線按指數規律衰減。NMR信號ENMR為：

(1)

圖1 MRS方法激發/接收NMR信號時序圖Fig.1 Transmitting and receiving NMR signal in time domain of MRS

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

1.2.1 噪聲種類

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

Eobs=ENMR+Espike+Eharmonic+Erandom,

(2)

式中:Eobs為儀器接收到的信號數據,ENMR為磁共振信號,Espike為尖峰噪聲,Eharmonic為諧波噪聲,Erandom為隨機噪聲。

1.2.2 多源諧波電磁噪聲特征

在經典諧波噪聲模型條件下，2013年Larsen等[14]提出假設諧波噪聲的基頻在信號采集持續期間內是單一且穩定的，由此寫出如下諧波模型式：

(3)

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

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

假設有2個基頻不同的諧波干擾源，且在測量時間內不變，可以寫出2個諧波源的諧波模型式：

(4)

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

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

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

2.1 同步刪除法原理

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

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

(5)

由于式(4)中的f1、f2、Ak、Bk、φk和φk的確定是一個非線性優化問題，可以通過標準化改寫式(4)中的余弦項：

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

(6)

式中變量αk、βk、γk和δk可通過以下關系與對應諧波源的振幅和相位相關聯：

(7)

(8)

將式(6)中的余弦項展開變換，假設f1和f2已知，觀測信號與諧波模型的差值最小值問題就變成線性，可以用矩陣表示法寫成：

(9)

式中：[k0,k1,…,kN]表示構建諧波模型的階數，[n0,n1,…,nf]為采樣時間樣本。同時，可將式(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。利用最小二乘法求解該標準線性方程可得：

x=(ATA)-1ATb。

(11)

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

2.2 模擬退火法原理

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

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

(12)

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

3 仿真數據驗證

3.1 網格搜索同步刪除法

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

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

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

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

3.2 模擬退火同步刪除法

為了加快基頻的搜索和避免陷入局部極值點，利用模擬退火法來搜索干擾源諧波的基頻值。仿真條件同上，迭代次數200次，搜索路徑如圖5所示。圖5中黑色圓圈代表初始模型(50.000；50.000)，白色圓圈代表每次迭代被接受的模型，綠色圓圈代表迭代結束最終模型(50.031 7；49.984 6)，紅色線段代表基頻搜索路徑?？梢钥闯?，每次被接受的模型都在逐步逼近真實模型，雖然模擬退火法最終能找到最小值，但在有限的迭代次數時搜索路徑會存在一定的偏差。為了更好地評價去噪效果，特意引入均方根誤差(root mean square error，RMSE)這個指標參數[21]：

(13)

式中:N為儀器采集到的信號長度，Edenoised(i)是去噪后的第i個數據，ENMR(i)是真實NMR信號的第i個數據。

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

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

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

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

4 實測數據檢驗

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

通過分析野外實際工作環境和噪聲，發現2個較強的諧波干擾源。根據數據處理流程對噪聲數據進行了尖峰噪聲的處理后，采用100次迭代的模擬退火同步刪除法來壓制雙諧波噪聲(圖7)，再將去噪后的信號進行疊加、希爾伯特變換和非線性擬合，結果見圖8。

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

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

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

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

圖9 天門河灘MRS方法某測深點實測數據及其反演結果與鉆探結果對比Fig.9 Measured data, inversion results and drilling results of a sounding point by MRS method in River Beach, Tianmen

5 結論

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

1)針對磁共振測深信號中模型去噪方法無法有效壓制多源諧波噪聲的問題，利用網格搜索同步刪除法和模擬退火同步刪除法均能有效找到較為準確的基頻值，達到壓制多源諧波的目的。

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

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

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