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

1 引言

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

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

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

2 基于雙重?cái)?shù)字濾波的TEM正演計(jì)算

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

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

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

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

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

根據(jù)階躍場和諧變場的關(guān)系，得到頻率域的階躍場的表達(dá)式為

對式（3）和式（4）獲得的頻率域電磁響應(yīng)，通過拉普拉斯逆變換可獲得時(shí)間域電磁場的表達(dá)式。顯然，這類方法計(jì)算量太大，采用 Gaver-Stehfest算法可以快速實(shí)現(xiàn)頻率域電磁場向時(shí)間域電磁場的轉(zhuǎn)換[12]：

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

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

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

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

3 模擬退火法反演步驟

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

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

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

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

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

4 仿真與分析

4.1 SA算法的改進(jìn)

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

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

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

此優(yōu)化策略具有以下特點(diǎn)：（1）ny根據(jù)溫度來調(diào)節(jié)每次迭代的步長，高溫時(shí)進(jìn)行大范圍的搜索，低溫時(shí)減小搜索范圍，僅在當(dāng)前反演模型附近進(jìn)行搜索，ny能根據(jù)nT 實(shí)現(xiàn)高效率的全局搜索；（2）P能根據(jù)步長被拒絕的次數(shù)而動(dòng)態(tài)變化，當(dāng)次數(shù)增大時(shí)，說明當(dāng)前模型遠(yuǎn)離真實(shí)值，要加大擾動(dòng)值P，進(jìn)而擴(kuò)大搜索范圍。每次內(nèi)部循環(huán)結(jié)束后，naN 被置零，SA算法重新進(jìn)行全局搜索。

4.2 3層模型的仿真結(jié)果與對比分析

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

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

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

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

表1 3層模型的反演結(jié)果對比

圖3 含有異常和沒有異常的階躍響應(yīng)對比圖

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

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

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

以下是對反演結(jié)果的分析：

（1）SA算法有效實(shí)現(xiàn)了對TEM方法中xE 仿真數(shù)據(jù)的反演，得到的模型參數(shù)與真實(shí)值基本一致。通過與Marquardt方法反演結(jié)果的對比，進(jìn)一步證實(shí)了SA算法搜索全局最優(yōu)解的能力。

（2）表1的反演結(jié)果表明，參數(shù)2ρ與其它參數(shù)相比，誤差率較大。這是因?yàn)樵趚E場值的正演計(jì)算中，參數(shù)2ρ對場值的貢獻(xiàn)較小，相關(guān)性較弱，所以xE的觀測數(shù)據(jù)對其不夠敏感，誤差相對較大。

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

4.3 5層模型的仿真結(jié)果與對比分析

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

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

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

以下是對反演結(jié)果的分析：

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

圖6 SA算法反演參數(shù)與真實(shí)參數(shù)的正演計(jì)算對比

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

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

表2 5層模型的反演結(jié)果對比

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

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

5 結(jié)論

本文提出了基于雙重?cái)?shù)字濾波的正演方法，實(shí)現(xiàn)了快速、精確的正演計(jì)算。將SA算法應(yīng)用于TEM方法的反演中，修改了算法的搜索策略，與傳統(tǒng)反演方法的對比實(shí)驗(yàn)驗(yàn)證了 SA算法的有效性。主要結(jié)論如下：（1）對于地下高阻性異常，xE 比其它場分量含有更豐富的地電信息，但其數(shù)據(jù)解釋難度較大。本文基于雙重?cái)?shù)字濾波方法的正演計(jì)算，簡化了計(jì)算過程，可快速獲得較高精度的電磁響應(yīng)，給xE場的數(shù)據(jù)反演提供了理論支撐。（2）SA 算法的搜索策略可以有效地實(shí)現(xiàn)層狀大地模型的反演，在搜索全局最優(yōu)解的過程中，目標(biāo)函數(shù)不是依次減小，而是波動(dòng)式減小，算法可以跳出局部極小值，實(shí)現(xiàn)全局優(yōu)化。（3）在正演計(jì)算過程中，參數(shù)對電磁場的貢獻(xiàn)越大，其反演的精度越高；加入隨機(jī)噪聲后，SA算法的目標(biāo)函數(shù)依然能收斂，表明 SA算法有一定的抗干擾能力。（4）隨著模型層數(shù)的增加，反演難度加大，但是，對于TEM方法常用的3層及其5層模型，SA算法都能獲得比較真實(shí)的地電模型參數(shù)。由于 SA算法進(jìn)行隨機(jī)搜索和步長的動(dòng)態(tài)調(diào)整，以避開局部極小值，SA算法需要的迭代次數(shù)較大，相應(yīng)的反演耗時(shí)較長，反演計(jì)算時(shí)間在數(shù)個(gè)小時(shí)之內(nèi)；但是，這仍然滿足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] 張建國, 武欣, 齊有政, 等. 時(shí)間域編碼電磁勘探方法研究[J].雷達(dá)學(xué)報(bào), 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].中國有色金屬學(xué)報(bào), 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] 湯井田, 任政勇, 化希瑞. 地球物理學(xué)中的電磁場正演與反演[J]. 地球物理學(xué)進(jìn)展, 2007, 22（4）: 1181-1194.

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

[8] 向陽, 于鵬, 陳曉, 等. 大地電磁反演中改進(jìn)的自適應(yīng)正則化因子選取[J]. 同濟(jì)大學(xué)學(xué)報(bào)（自然科學(xué)版）, 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] 孫娟娟, 黃宗潘, 楊大成. 基于模擬退火的移動(dòng)業(yè)務(wù)分布預(yù)測算法[J]. 電子與信息學(xué)報(bào), 2008, 30（9）: 2225-2228.

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

[11] 張宇, 張曉娟, 方廣有. 大尺度分層介質(zhì)電特性參數(shù)的反演方法研究[J]. 物理學(xué)報(bào), 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.