時間域電磁勘探數據的模擬退火法反演研究

張建國 武 欣 趙海濤② 齊有政② 方廣有

1 引言

電磁勘探方法對高阻性異常體有較強的探測能力，已成為油氣資源勘探的一種有效方法[1,2]。利用接地長導線向地下發送電磁場（一次場），地下異常體中激勵的感應渦流將產生感應電磁場（二次場），利用接收機觀測感應電磁場，再通過數據解釋來推斷地下電性結構的空間分布，進而確定地下異常體的位置和規模。時間域電磁法（Time-domain ElectroMagnetic, TEM）在一次場的關斷間歇來觀測二次場，得到的觀測數據是沒有一次場“污染”的電磁響應，通過數據解釋可獲得更真實的地電模型參數[3,4]。

目前，電磁勘探中地電模型的精確反演是一大難題。雖然一些學者開展了電磁法的多維反演理論研究，但需要指出的是，這些方法應用于實踐還存在許多問題，能有效應用到工程中的仍然只是1維反演[5,6]。其中1維線性化反演方法仍然是電磁勘探中普遍采用的反演方法[7,8]。然而，由于其反演過程是采用梯度搜索方法對非線性模型進行線性化近似，這必然從理論上帶來了誤差；同時，由于地電模型的非惟一性，此類方法容易收斂于局部極小值，并高度依賴于初始模型的選擇，其數據解釋結果容易偏離真實的大地模型。

因此，本文開展 TEM 數據的非線性反演方法研究，將模擬退火（Simulated Annealing, SA）算法應用于TEM數據最優解的搜索[911]-。雖然SA算法已開始在地震勘探中應用，但是由于 TEM 地電模型的高度非線性和非惟一性，其在 TEM 數據反演中的應用鮮有文獻討論。本文首先采用雙重數字濾波方法提高TEM正演計算的速度，縮短SA算法反演的迭代時間。同時，改善 SA算法的搜索控制策略，使算法在全局最優解的搜索過程中能自適應地調節搜索步長，并逐步縮小搜索范圍，提高算法的搜索效率。與Marquardt反演方法（線性化反演方法）進行對比實驗，結果表明 SA算法能在不依賴于初始模型和梯度的條件下，有效地擺脫局部最優，搜索到全局最優解。

2 基于雙重數字濾波的TEM正演計算

如圖1所示，TEM方法通過接地電極AB向地下發射雙極性方波形式的激勵信號，雙極性方波的每個下降沿可視為地電模型的階躍激勵信號，同時利用多通道接收機采集具有地下異常信息的電磁響應數據。其中，ρn和hn分別為第n層介質的電阻率和厚度（n = 1 ,2,… , N ）, T和τ分別為雙極性方波的周期長度和占空比，角度Φ為發射點與接收點連線與x軸的夾角。

圖1 地電模型和TEM方法示意圖

地電模型的頻率域電磁場滿足的赫姆霍茲（Helmholtz）方程為

其中，k為波數。由于TEM方法是對時間域的階躍場進行觀測，因此正演的思路為：首先求解頻率域的諧變場，然后根據階躍場和諧變場的關系，得到頻率域的階躍場 G （ ω） = F （ ω）/（i ω） ，最后利用傅里葉逆變換或者拉普拉斯逆變換，獲得G（ω）對應的時間域的階躍場g（ t） 。

對于均勻各向同性的層狀大地模型，電偶極源激勵的諧變場（正演中只考慮水平電場 Ex和垂直磁場 Hz）的表達式為

根據階躍場和諧變場的關系，得到頻率域的階躍場的表達式為

對式（3）和式（4）獲得的頻率域電磁響應，通過拉普拉斯逆變換可獲得時間域電磁場的表達式。顯然，這類方法計算量太大，采用 Gaver-Stehfest算法可以快速實現頻率域電磁場向時間域電磁場的轉換[12]：

其中，Kn為預先確定的Gaver-Stehfest算法的數字濾波系數，Pn= n ? ln 2 /t 。ex（ t）,hz（ t） 為頻率域階躍場對應的時間域表達式，是TEM方法的正演結果。

對于式（3）和式（4）中 Bessel函數非線性積分的計算，使用快速Hankel變換方法可以獲得令人滿意的數值精度。將積分中Bessel函數的兩個變量（λ,r）通過指數代換，進而把非線性積分計算轉換成卷積形式。然后將卷積離散化，通過預先確定的數字濾波系數，獲得高精度的數值解[13,14]。

其中 H1,H2和 H3分別為預先確定的快速Hankel變換的數字濾波系數，Δ為濾波系數的坐標間隔。濾波系數越長，數值解的精度越高。

比較式（7）和式（8）可知，與 Hz相比， Ex包含的地下信息更為豐富。但是，由于反演解釋的復雜性不同，TEM 方法通常采用 Hz作為觀測數據。本文采用雙重數字濾波的正演方法，能有效地簡化正演計算過程，實現電磁場的快速、精確計算，本文反演計算針對 Ex場進行。在搭載 Intel Core i3-2120 CPU（主頻3.3 GHz），內存2 GB的PC機上，針對3層地電模型，20個時間點的正演計算耗時僅0.056 s。

3 模擬退火法反演步驟

層狀地電模型的參數包括每層的電阻率ρn和厚度 hn，把這些參數看作一個數組向量，則地電模型可記為x。由于地電模型是高度非線性的，所以通過電磁場的觀測值進行反演，推導地電模型的參數的過程屬于非線性最優化問題。將其轉化為一個迭代求解過程，且每次迭代過程中都通過一定的方法控制變量的變化，使目標函數S（x）隨著迭代的進行而逐步下降，從而達到逼近最優解的目的。地電模型的最優解 xopt應滿足： S （xopt）= min（S （x ） ） ，即

式中rx為真實模型，obv（）dx為觀測數據，（）d x表示根據真實模型進行正演計算得到的電磁場理論數據。

模擬退火（SA）算法應用于求解多元函數全局最小值，不依賴于初始模型的選擇，并能有效地擺脫局部最小值，同時不需要計算復雜的雅克比矩陣[15,16]。在地震勘探領域中，SA算法作為一種全局優化算法，得到越來越多的應用。但是，由于地電模型的高度非線性和非惟一性，采用 SA算法的TEM數據的反演，鮮有研究文章發表。本文基于雙重數字濾波方法的正演計算，通過改進的SA算法，實現了TEM方法的準確、快速反演。

當 SA算法的目標函數值較大時，溫度取較大值，那么模型迭代的步長可以取較大值；反之亦然。因此，反演不會陷入局部極小值，保證了反演過程的全局最優化。SA算法通過溫度來控制迭代過程，使用收斂閾值來判斷迭代結果是否被接受。SA算法的流程圖如圖2所示，其中：n為內部循環迭代次數，m為外部循環迭代次數，mT 為溫度，η為每次迭代中產生的一個隨機數。當目標函數（）S x小于收斂閾值或者m大于最大循環次數時，迭代過程結束。使（）S x小于收斂閾值的地電參數，即為SA算法的反演結果。

圖2 模擬退火算法的流程框圖

4 仿真與分析

4.1 SA算法的改進

傳統的 SA算法采用隨機擾動的搜索策略，效率較低，不能滿足 TEM 方法大數據量、高復雜度的非線性反演的需求。本文聯合柯西概率分布與“擾動因子”，來設定地電模型的步長：

其中，y= s ign（η - 0 .5）T [ （1 + 1 / T ）2ηn-1-1]為基

nnnn于柯西概率分布的修正參數，擾動因子為 P =exp（5.0 ?Nna/Ns），式中ηn為（0,1）之間的隨機數，Tn為溫度控制變量， Ns為內部循環的總次數， Nna為每個內部循環中步長未被接受的次數，系數 5.0是經驗參數。

此優化策略具有以下特點：（1）ny根據溫度來調節每次迭代的步長，高溫時進行大范圍的搜索，低溫時減小搜索范圍，僅在當前反演模型附近進行搜索，ny能根據nT 實現高效率的全局搜索；（2）P能根據步長被拒絕的次數而動態變化，當次數增大時，說明當前模型遠離真實值，要加大擾動值P，進而擴大搜索范圍。每次內部循環結束后，naN 被置零，SA算法重新進行全局搜索。

4.2 3層模型的仿真結果與對比分析

針對油氣資源勘探，設置一個含有高阻油氣層的1維層狀大地模型，以驗證本文反演方法的有效性。地電模型共含有3層水平各向同性介質，表層介質的電阻率 ρ1= 5 0.0 Ω ?m ，厚度 h1= 2 00.0 m；中間層介質表征油氣層的電學特性，設置其電阻率為 ρ2= 1 000.0 Ω ?m ，厚度 h2= 5 0.0 m；底層為均勻半空間，介質的電阻率 ρ3= 1 00.0 Ω ?m ，這是TEM方法典型的地電模型。采用 Ex作為觀測數據，Ex比其它場分量含有更多的地下異常信息。發射電流I= 100A，雙極性方波的頻率 f = 1 Hz（占空比τ為 0.25），發射線長度AB = 2 00m，接收線長度MN = 1 00 m，偏移距 r = 2 km。為了檢驗TEM方法對高阻異常體的探測能力，針對以上模型，采用基于雙重數字濾波的正演方法，在有、無異常條件下分別進行計算，得到的觀測數據如圖3所示。由圖3可知，兩條曲線在0.1～10 ms時間段內的差異非常明顯，證明了 TEM 方法對高阻異常體有較強的探測能力。

將上述含有高阻異常體的正演仿真結果（0.1～10 ms時間段）作為反演對象，分別使用改進的 SA方法和傳統線性化反演方法的代表Marquardt方法，對大地模型的參數（ρ1, ρ2, ρ3,h1和 h2）進行反演，確定地電模型的參數值，并對比二者結果。基于上述模型參數和正演仿真數據，先使用Marquardt方法進行1維反演。設定地電模型的初始值為 x0=（80, 300, 800, 100, 150），進行10次反演計算，取10次結果的平均值作為地電模型的參數估計值。目標函數隨著迭代次數的增加而逐漸減小，收斂曲線如圖4所示，目標函數最終收斂于 2 × 1 0-4，得到的反演結果平均值如表1（電阻率單位為Ω?m，厚度單位為m）所示。由表1結果可知，Marquardt方法的反演結果最終陷于局部最小值，沒有獲得全局最優解，偏離真實模型參數。

采用相同的模型參數和正演仿真數據，再采用SA算法進行反演計算。仍然設定地電模型的初始值為 x0=（80, 300, 800, 100, 150），模型空間為真實值的±40%，內部循環次數N= 2 00，外部循環最大次數 M = 5 00。根據反演的經驗，外部循環的溫度控制變量的初始值設定為 T0= 1 0，并按照 Tm=T0?0. 9 7m的指數形式變化；溫度控制變量的設置要多次嘗試，以適應模型的擾動和目標函數的遞減趨勢。仍然進行10次反演計算，取結果的平均值。

目標函數S（x）的收斂閥值設為S= 1 0-6，即當0S（x）＜S0時，SA算法停止搜索，認為當前的模型是滿足反演收斂條件的真實模型。圖5為目標函數收斂曲線，相應的反演結果平均值如表1所示。需要指出的是，基于梯度的反演方法的目標函數隨著迭代次數的增加而逐漸減小（圖4），SA算法的目標函數是波動式的逐步減小（圖5），正是由于波動式的收斂曲線使得SA算法能跳出局部極小值。將SA算法的反演結果代入正演計算中，并與真實模型參數的正演計算結果對比（如圖 6所示），二者曲線基本一致，證實了反演結果的有效性。

表1 3層模型的反演結果對比

圖3 含有異常和沒有異常的階躍響應對比圖

圖4 3層模型Marquardt方法的收斂曲線

圖5 3層模型SA方法的收斂曲線

在搜索過程中，SA算法通過改進的搜索策略能成功跳出局部極小值。在跳出局部極小值的過程中，目標函數是增大的；但隨著迭代的進行，搜索會逐步向全局極小值靠近，目標函數也隨之逐步減小，最終在全局極小值處收斂。為了驗證算法的有效性，還對反演數據加入了5%的隨機噪聲，進行10次反演計算，同樣取平均值得到反演結果如表1所示。

以下是對反演結果的分析：

（1）SA算法有效實現了對TEM方法中xE 仿真數據的反演，得到的模型參數與真實值基本一致。通過與Marquardt方法反演結果的對比，進一步證實了SA算法搜索全局最優解的能力。

（2）表1的反演結果表明，參數2ρ與其它參數相比，誤差率較大。這是因為在xE場值的正演計算中，參數2ρ對場值的貢獻較小，相關性較弱，所以xE的觀測數據對其不夠敏感，誤差相對較大。

（3）加入噪聲后，反演誤差增大，但是反演誤差小于5%，仍然在可接受的范圍內，通過多次反演取平均值的方法，提高了計算結果的精度。

4.3 5層模型的仿真結果與對比分析

對于 TEM 數據的反演方法驗證，一般只采用典型的3層地電模型。為了進一步測試本文改進的SA算法對TEM數據反演的有效性，將3層模型擴展至5層模型，在高阻油氣層的上方和下方各加 1層電阻率為300 Ω?m的介質，地電模型的初始值為x0=（70, 150, 200, 600, 700, 70, 200, 700, 150），仍然采用3層模型中的收發參數。

采用Marquardt方法對上述模型的正演數據進行反演操作，同樣取10次反演結果的平均值作為地電模型參數的估計值。目標函數的收斂曲線如圖 7所示，Marquardt方法在迭代21次后，目標函數最終收斂于 1 .6 × 1 0-4，得到的反演結果平均值如表2（電阻率單位為Ω?m，厚度單位為m）所示。與3層模型的反演結果相比，5層模型的反演結果誤差更大，Marquardt方法的搜索過程陷于局部最小值，反演結果偏離全局最優解。

采用SA算法對相同的正演數據進行反演操作，地電模型的初始值也與Marquardt方法相同，模型空間為真實值的±40%，內部循環次數N= 2 00，外部循環最大次數 M = 1 000。溫度變量的初始值設置為 T0= 1 0，并按照 Tm= T0?0. 9 7m的指數形式變化。仍然進行10次反演計算，取結果的平均值。目標函數S（x）的收斂閾值仍然為S= 1 0-6，圖8為目標函0數收斂曲線，相應的反演結果平均值如表2所示。SA算法的目標函數仍然波動式減小，在迭代過程中避開多個局部極小值，進而逼近全局最優，獲得比較真實的地電模型估計值。SA算法在外部循環迭代474次時，獲得低于收斂閥值的目標函數值，此時的反演結果是TEM數據對應的全局最優解。

以下是對反演結果的分析：

（1）SA算法基本實現了5層模型TEM數據的反演，獲得的模型參數估計值能反映真實模型的特性，模型中介質電阻率由上而下為“低-高-低”的變化趨勢，在反演結果中得到了反映；與Marquardt方法相比，SA方法對TEM數據反演更準確、更有效。

圖6 SA算法反演參數與真實參數的正演計算對比

圖7 5層模型Marquardt方法的收斂曲線

圖8 5層模型SA算法的收斂曲線

表2 5層模型的反演結果對比

（2）由于5層模型的參數較多，相應的反演難度較大，因此，反演結果的誤差增大；但是，SA方法仍然能有效獲得模型的關鍵信息，如第3層（高阻油氣層）和第1層（表層）的電阻率及其厚度信息。

（3）隨著模型參數的增多，反演難度加大，在加入5%噪聲條件下，反演結果與真實值有一定偏差，SA方法對噪聲的抑制能力降低。

5 結論

本文提出了基于雙重數字濾波的正演方法，實現了快速、精確的正演計算。將SA算法應用于TEM方法的反演中，修改了算法的搜索策略，與傳統反演方法的對比實驗驗證了 SA算法的有效性。主要結論如下：（1）對于地下高阻性異常，xE 比其它場分量含有更豐富的地電信息，但其數據解釋難度較大。本文基于雙重數字濾波方法的正演計算，簡化了計算過程，可快速獲得較高精度的電磁響應，給xE場的數據反演提供了理論支撐。（2）SA 算法的搜索策略可以有效地實現層狀大地模型的反演，在搜索全局最優解的過程中，目標函數不是依次減小，而是波動式減小，算法可以跳出局部極小值，實現全局優化。（3）在正演計算過程中，參數對電磁場的貢獻越大，其反演的精度越高；加入隨機噪聲后，SA算法的目標函數依然能收斂，表明 SA算法有一定的抗干擾能力。（4）隨著模型層數的增加，反演難度加大，但是，對于TEM方法常用的3層及其5層模型，SA算法都能獲得比較真實的地電模型參數。由于 SA算法進行隨機搜索和步長的動態調整，以避開局部極小值，SA算法需要的迭代次數較大，相應的反演耗時較長，反演計算時間在數個小時之內；但是，這仍然滿足TEM方法的要求。

[1] Zhang Jian-guo, Wu Xin, Qi You-zheng, et al.. Research on 3D marine electromagnetic interferometry with synthetic sources for suppressing the airwave interference[J]. Applied Geophysics, 2013, 10（4）: 373-383.

[2] 何繼善. 海洋電磁法原理[M]. 北京: 高等教育出版社, 2012:3-4 .

[3] 張建國, 武欣, 齊有政, 等. 時間域編碼電磁勘探方法研究[J].雷達學報, 2014, 3（2）: 158-165.Zhang Jian-guo, Wu Xin, Qi You-zheng, et al.. Investigation on time domain coded electromagnetic exploration method[J].Journal of Radars, 2014, 3（2）: 158-165.

[4] Ziolkowski A, Hobbs B A, and Wright D. Multitransient electromagnetic demonstration survey in France[J].Geophysics, 2007, 72（4）: F197-F209.

[5] 劉穎, 李予國, 柳建新, 等. 海洋可控源電磁場的一維反演[J].中國有色金屬學報, 2013, 23（9）: 2551-2556.Liu Ying, Li Yu-guo, Liu Jian-xin, et al.. One-dimentional inversion of marine controlled-source electromagnetic fields[J].The Chinese Journal of Nonferrous Metals, 2013, 23（9）:2551-2556.

[6] 湯井田, 任政勇, 化希瑞. 地球物理學中的電磁場正演與反演[J]. 地球物理學進展, 2007, 22（4）: 1181-1194.

[7] 毛立峰, 王緒本, 陳斌. 直升機航空瞬變電磁自適應正則化一維反演方法研究[J]. 地球物理學進展, 2011, 26（1）: 300-305.

[8] 向陽, 于鵬, 陳曉, 等. 大地電磁反演中改進的自適應正則化因子選取[J]. 同濟大學學報（自然科學版）, 2013, 41（9）:1429-1434.Xiang Yang, Yu Peng, Chen Xiao, et al.. An improved adaptive regularized parameter selection in magnetotelluric inversion[J]. Journal of Tongji University（Natural Science）,2013, 41（9）: 1429-1434.

[9] 孫娟娟, 黃宗潘, 楊大成. 基于模擬退火的移動業務分布預測算法[J]. 電子與信息學報, 2008, 30（9）: 2225-2228.

[10] 葛建軍, 張春城. 基于模擬退火算法的機載脈沖多普勒雷達中重復頻率選擇研究[J]. 電子與信息學報, 2008, 30（3）:573-575.

[11] 張宇, 張曉娟, 方廣有. 大尺度分層介質電特性參數的反演方法研究[J]. 物理學報, 2013, 62（4）: 044204.Zhang Yu, Zhang Xiao-juan, and Fang Guang-you. A data inversion method for electromagnetic scattering from large-scale layered medium[J]. Acta Physica Sinica, 2013,62（4）: 044204.

[12] Antonijevic S and Poljak D. A novel time-domain reflection coefficient function: TM case[J]. IEEE Transactions on Electromagnetic Compatibility, 2013, 55（6）: 1147-1153.

[13] Ding Ping-ping, Qiu Cheng-wei, Zouhdi S, et al.. Rigorous derivation and fast solution of spatial-domain green’s functions for uniaxial anisotropic multilayers using modified fast hankel transform method[J]. IEEE Transactions on Microwave Theory and Techniques, 2012, 60（2）: 205-217.

[14] Kong Fan-nian. Hankel transform filters for dipole antenna radiation in a conductive medium[J]. Geophysical Prospecting,2007, 55（1）: 83-89.

[15] Freethy S J, Shevchenko V F, and Vann R G L. Optimization of wide field interferometric arrays via simulated annealing of a beam efficiency function[J]. IEEE Transactions on Antennas and Propagation, 2012, 60（11）: 5442-5446.

[16] Wang Ruo, Yin Chang-chun, Wang Miao-yue, et al..Simulated annealing for controlled-source audio-frequency magnetotelluric data inversion[J]. Geophysics, 2013, 77（2）:E127-E133.