基于自適應(yīng)純形模擬退火法一維大地電磁測深視電阻率和相位反演研究

孫歡樂，王世彪，郭榮文，周紹民，柳建新*

(1.中南大學(xué) a.地球科學(xué)與信息物理學(xué)院，b.有色資源與地質(zhì)災(zāi)害探查湖南省重點實驗室，長沙 410083；2.廣東有色工程勘察設(shè)計院，廣州 510080)

?

基于自適應(yīng)純形模擬退火法一維大地電磁測深視電阻率和相位反演研究

孫歡樂1a,1b，王世彪2，郭榮文1a,1b，周紹民1a,1b，柳建新1a,1b*

(1.中南大學(xué) a.地球科學(xué)與信息物理學(xué)院，b.有色資源與地質(zhì)災(zāi)害探查湖南省重點實驗室，長沙 410083；2.廣東有色工程勘察設(shè)計院，廣州 510080)

針對大地電磁測深反演中線性化方法容易陷入局部極值，而全局優(yōu)化方法收斂慢等問題，這里采用自適應(yīng)純形模擬退火綜合優(yōu)化方法進(jìn)行大地電磁測深數(shù)據(jù)反演。該優(yōu)化方法綜合了下降純形法和模擬退火法各自的優(yōu)點，已被證明具有全局搜所能力和收斂速度快的特點，并且實現(xiàn)了視電阻率和相位同時反演，減小多解性的同時提高了反演的分辨率。通過H、K和HKH型模型的數(shù)值計算，驗證了這種綜合優(yōu)化方法的搜索效率和全局收斂性；所有模型合成數(shù)據(jù)的反演結(jié)果都能較好地反映真實模型的結(jié)構(gòu)特征。對于中間層為相對高阻的模型，雖然相對其他層結(jié)構(gòu)恢復(fù)差些，但是反演后該層結(jié)構(gòu)基本得到恢復(fù)。

大地電磁； 視電阻率； 相位； 下降純形法； 模擬退火法； 自適應(yīng)純形模擬退火法

0 引言

大地電磁測深(MT)是發(fā)展較早的一種地球物理勘探方法，在探測地?zé)?、金屬礦勘探和地震預(yù)報等方面得到廣泛地應(yīng)用[1-2]。近些年來該方法發(fā)展迅速，許多學(xué)者一直致力于大地電磁測深反演方法的研究，對于一維MT，目前應(yīng)用最多的反演方法大都屬于線性或擬非線性方法(馬夸特方法和廣義逆矩陣法)[3-4]。但是大地電磁測深反演問題是非線性問題，而且具有非唯一性，線性或擬非線性反演方法依賴于初始模型的選擇，容易陷入局部極值，最終把局部最優(yōu)解當(dāng)作全局最優(yōu)解，導(dǎo)致后期資料的解釋不準(zhǔn)確[5]。由此，專家學(xué)者們提出了應(yīng)用非線性全局最優(yōu)方法求解地球物理反演問題[6]，其中，模擬退火法得到了許多地球物理學(xué)家的關(guān)注和重視[7]。模擬退火法受統(tǒng)計熱力學(xué)現(xiàn)象的啟發(fā)[8]，1953年Metropolis等人[9]將模擬退火算法的思想提出，但當(dāng)時并未得到廣泛地應(yīng)用，30年后，Kirkpatrick等人[10]將模擬退火法的思想成功地應(yīng)用于組合優(yōu)化問題，隨后模擬退火法得到了很大的發(fā)展。在地球物理領(lǐng)域，Rothman[11-12]在地震資料處理的剩余靜校正問題中首次應(yīng)用了模擬退火法的思想，其后Sen和Stoffa等人[13-14]將其用于一維地震波形反演、一維電測深反演。后來，模擬退火法在MT[15]、重力磁法勘探資料等反演問題中[16]得到了廣泛研究與應(yīng)用。

實際中的地球物理反演是一個多參數(shù)、非線性優(yōu)化問題，所采用的目標(biāo)函數(shù)，即衡量由理論模型得出的正演響應(yīng)值與觀測值的吻合程度，常具有多個局部極值[17]，雖然模擬退火法隨機(jī)搜索，可以避免局部極小值的影響，最終得到全局最優(yōu)解，但是模擬退火法優(yōu)化方法是在可行解空間中隨機(jī)搜索，通過Metropolis法[18]來控制新解的取舍，所以其迭代速度相對較慢。此外，模擬退火法是否收斂于全局最優(yōu)解受多種因素的影響，比如初始溫度的選擇、溫度的下降模式、隨機(jī)擾動的方式等。溫度選擇越高，下降越慢，模型空間搜索的越徹底，可以避免陷入局部極小模型，但是耗費大量的計算時間，所以模擬退火法的全局收斂較慢。

針對上述問題，作者在研究地球物理各種優(yōu)化算法的基礎(chǔ)上，引進(jìn)了自適應(yīng)純形模擬退火綜合方法(Adaptive simplex simulated annealing，ASSA)，該方法是純形下降優(yōu)化方法和模擬退火法的有機(jī)結(jié)合，最初由Press等[19]提出解決NLP問題。純形下降法最早由 Spednley等人提出，Nelder 等[20]作了改進(jìn)，純形下降優(yōu)化方法可以實現(xiàn)直接快速地搜索到極小值，但是容易受到局部極小值的影響，所以容易陷入局部極值。因此，可以結(jié)合二者的優(yōu)點，協(xié)同快速地求解反演問題的全局最優(yōu)解，自適應(yīng)純形模擬退火優(yōu)化方法已在一維聲波和地聲學(xué)反演方面得到了成功地應(yīng)用[17,21]。一維MT反演通常只利用視電阻率或者阻抗幅值，往往忽略阻抗相位，其實阻抗相位包含著大量的地下電性信息。對于一維地電結(jié)構(gòu)來說，相位和視電阻率資料不是相互獨立的，其具有一定的相關(guān)性[2]。作者利用視電阻率數(shù)據(jù)和阻抗相位數(shù)據(jù)進(jìn)行同時反演，二者同時反演可以減小反演的多解性，提高分辨能力，還對壓制靜態(tài)位移具有積極的意義。這里通過對一維大地電磁數(shù)據(jù)測深反演，驗證了自適應(yīng)純形模擬退火綜合反演方法的全局搜索能力和效率，同時驗證了視電阻率數(shù)據(jù)和阻抗相位數(shù)據(jù)同時處理可以提高反演的精度。

1 一維大地電磁測深正演

1.1 正演公式

假設(shè)地下電性結(jié)構(gòu)是水平層狀分布的，如果地電結(jié)構(gòu)共有n層，則有2n-1個參數(shù)， hi(i=1,2,…,n-1)代表第i層的厚度；ρi(i=1,2,…,n)代表第i層的電阻率，最后一層為均勻半空間，則hn=∞。計算視電阻率ρs和阻抗相位φ的表達(dá)式如式(1)、式(2)所示[2-3]。

(1)

(2)

其中：ω表示角頻率；μ0表示導(dǎo)磁率；Z表示地面的波阻抗值。

公式(3)建立了視電阻率、阻抗相位和頻率及層參數(shù)間的關(guān)系[5]。

ρs=f(h1,ρ1,h2,ρ2,…,hn-1,ρn-1,ρn,f)

φ=g(h1,φ1,h2,φ2,…,hn-1,φn-1,φn,f)

(3)

1.2 目標(biāo)函數(shù)

目標(biāo)函數(shù)(能量函數(shù))[22]通過模型正演響應(yīng)值、觀測值及其數(shù)據(jù)噪聲的標(biāo)準(zhǔn)差來建立，具體表達(dá)式如式(4)所示。

(4)

其中：ρdi和φdi分別為頻點fi的觀測視電阻率值、觀測阻抗相位值；ρmi和φmi分別為頻點fi的模型視電阻率值、模型阻抗相位值；σρi和σφi分別表示相應(yīng)第i個視電阻率或相位數(shù)據(jù)噪聲的標(biāo)準(zhǔn)方差；n為頻點個數(shù)。

2 模擬退火算法

模擬退火法(Simulated Annealing， SA)來源于熱力學(xué)的退火原理[23-25]：高溫物體具有一定的自由能，物體內(nèi)部晶體雜亂無序，處于十分不穩(wěn)定的狀態(tài)；不斷地降低物體的溫度，隨著溫度的降低物體的自由能也有所下降，物體內(nèi)部晶體從無序逐漸有規(guī)則的排列。當(dāng)溫度下降到最小值時，此時物體內(nèi)部的自由能達(dá)到最小，物體內(nèi)部的晶體有序排列，整個物體處于一種非常穩(wěn)定的狀態(tài)。溫度最低時，物體所具有的最小自由能就是優(yōu)化問題所求得最優(yōu)解。

模擬退火算法的實現(xiàn)步驟為：

1)設(shè)置全局可行解的搜索空間。在可行解的空間中初始化一個模型M0、設(shè)定起始溫度T0、最大迭代次數(shù)NT、擬合誤差閥值e、已迭代次數(shù)k等參數(shù)并令k=0。設(shè)定溫度下降方式，不同反演問題選擇不同的溫度下降方式(指數(shù)下降型和雙曲線下降型)。

2)對第k次迭代的模型參數(shù)值Mk，計算溫度Tk和目標(biāo)函數(shù)值E(Mk)，對模型Mk在可行解搜索空間里進(jìn)行隨機(jī)擾動形成新模型Mk+1，并計算目標(biāo)函數(shù)E(Mk+1)，E為目標(biāo)函數(shù)。

3)采用Metropolis抽樣準(zhǔn)則對模型選擇。令ΔE=E(Mk+1)-E(Mk)，若ΔE<0，新模型Mk+1被接受；否則，新模型Mk+1按照概率P(ΔE)=exp(-ΔE/T)被接受，產(chǎn)生隨機(jī)數(shù)η，η∈[0,1]，若η≤P，則新模型Mk+1被接受，否則，拒絕接受新模型Mk+1，轉(zhuǎn)步驟2)。

4)計算目標(biāo)函數(shù)E(Mk+1)，如果E(Mk+1)

3 下降純形法

下降純形法(Downhill Simplex， DS)是一種高效的近似梯度尋優(yōu)方法，該方法不需參數(shù)空間導(dǎo)數(shù)或靈敏度矩陣[26-27]參加運算，可以直接、快速地沿極小值的方向搜索，但容易陷入局部的極小值。其思想是在R維模型空間中，假定一個由R+1個點組成的初始純形空間(如圖1(a)R=3)，每個點都對應(yīng)一個函數(shù)值并按照其大小依次排列[21]，得到最高點(函數(shù)值最大)和最低點(函數(shù)值最小)，然后對純形進(jìn)行反射、收縮及延伸等步驟并計算對應(yīng)的函數(shù)值，舍棄惡化點，保留優(yōu)化點，最終構(gòu)造成新的純形空間，由此逐步逼近極小值點。

下降純形法步驟為[21,28]：

1)算法開始，具有最大函數(shù)值的點(高點)通過R維單純形的背向面移到另一側(cè)(圖1(b))以改善高點的函數(shù)值。

2)若反射后的函數(shù)值比最大值要小，那么這個點就被接受。若反射后的函數(shù)值在純形中最小，那么新的點將從高點的背向反射并延伸兩倍(圖1(c))。

圖1 下降純形法示意圖[19]Fig.1 Downhill simplex method diagram(a)初始純形； (b)反射；(c)反射和延伸；(d)收縮；(e)多向收縮

3)若新的點在純形中仍然具有最大的函數(shù)值，那么這個反射后的點將被丟棄，同時從高點沿反射面方向的收縮(圖1(d))。

4)若上述所有的步驟都無法使函數(shù)值降低到比最大值小，那么純形上R個點將沿低點方向進(jìn)行兩倍的收縮(圖1(e))。

重復(fù)以上四個步驟直到在純形中最高點函數(shù)值Ehigh和最低點函數(shù)值Elow之間滿足式(5)收斂標(biāo)準(zhǔn)[27]。

(5)

在反演中，這里采用ε=10-5。

4 自適應(yīng)純形模擬退火法

自適應(yīng)純形模擬退火法(ASSA)是一種綜合優(yōu)化算法，綜合了下降純形法(DS)和模擬退火法(SA)[29]。模擬退火法是一種利用Metropolis準(zhǔn)則在可行解空間隨機(jī)搜索的優(yōu)化方法，它可以避免局部極小值的干擾，最終可以搜索到全局最優(yōu)解，但其耗費大量時間、收斂速度慢[28]。下降純形法能直接、快速地搜索到極小值，但當(dāng)目標(biāo)函數(shù)有不止一個極小值時，不同的初值選擇會導(dǎo)致不同的反演結(jié)果，容易導(dǎo)致局部最優(yōu)解代替全局最優(yōu)解[30]。因此，綜合二者的優(yōu)點，既可以提高搜索效率，又避免了陷入局部極小值。

自適應(yīng)純形模擬退火法的基本思想為對于初始模型M0，首先采用純形下降法沿著梯度方向搜索到一個極小值點，然后利用模擬退火法在可行解空間內(nèi)隨機(jī)搜索，使其跳出該局部極小值點，當(dāng)搜索到一個比該極小值點更小的點時，立刻以新的極小值點為初始點再次采用純形下降法搜索新的極小值點附近的另外一個極小值點，如此交叉進(jìn)行純形下降法和模擬退火法，直到滿足精度終止條件，算法結(jié)束。

自適應(yīng)純形模擬退火優(yōu)化算法的實施步驟[28]：

1)設(shè)定初始值。設(shè)置控制參數(shù)T(溫度)的初值T0、初始步長V0、控制精度τ、最大迭代次數(shù)NT、已迭代次數(shù)k等參數(shù)并令k=0。隨機(jī)產(chǎn)生R+1個點的初始模型，其中能量值最小模型為M0，令EOPT=E(M0)，MOPT=M0。

2)下降純形法搜索。以M0為初始點，采用下降純形算法搜索該點周圍的一個極小值Mmin，并令Mk=Mmin，E(Mk)=E(Mmin)，MOPT=Mmin，EOPT=E(Mmin)。

3)模擬退火法。以Mk為初始點采用模擬退火算法，計算溫度Tk，在可行解搜索空間隨機(jī)產(chǎn)生新點Mk+1=Mk+rV0，其中r為隨機(jī)數(shù)，r∈(-1,1)，計算E(Mk+1)。

4)Metropolis接受準(zhǔn)則。若E(Mk+1)≤E(Mk)，新解Mk+1被接受；否則，新解Mk+1按照概率P(ΔE)=exp(-Δ/T)被接受：產(chǎn)生隨機(jī)數(shù)η，η∈[-1,1]，若η≤P，則新解Mk+1被接受。

5)重復(fù)步驟3)和步驟4)，直到E(Mk+1)

6)改變溫度Tk+1=update(Tk)，并令k=k+1。

自適應(yīng)純形模擬退火優(yōu)化當(dāng)中，下降純形法反射、延伸和收縮等每個步驟的函數(shù)值計算之前，對該模型的各參數(shù)設(shè)置Cauchy分布的擾動[21,26]。

(6)

c=λ1tan(λ2π/2)

(7)

式中:λ1和λ2是隨機(jī)分布在[-1,1]上的變量。

溫度的下降方式影響到自適應(yīng)純形模擬退火法的收斂性，ASSA中溫度利用指數(shù)下降模式，比較符合退火的實質(zhì)[21]，即

Tk=αkT0

(8)

式中:k為迭代次數(shù)(降溫次數(shù))；α是小于1的常數(shù)；T0為初始溫度，初始溫度應(yīng)該使80%以上的隨機(jī)擾動都被接受；Tk為迭代k次時的溫度；這里α取值0.998。

5 數(shù)值模擬計算實例

為了驗證自適應(yīng)純形模擬退火法優(yōu)化方法的搜索能力和全局收斂性，作者設(shè)計H型、K型和HKH共3個地電模型，模型參數(shù)設(shè)置如表1所示。頻率點數(shù)為40個，本文中的電阻率的單位為Ω·m，層厚度單位m。

表1 模型參數(shù)

5.1 計算實例對比

我們選擇一個H模型，參數(shù)設(shè)置如表1所示，分別利用下降純形法、模擬退火法和自適應(yīng)純形模擬退火法進(jìn)行反演，試比較三者的反演速度和擬合精度。實際采集的野外資料都含有一定的噪聲，這里采用理論合成模型加入適當(dāng)?shù)脑肼暷M野外實測數(shù)據(jù)，對該模型的正演響應(yīng)(視電阻率和相位)加入2%的高斯噪聲后再反演，三種反演方法結(jié)果如表2所示。

從表2可以看出，對于H型曲線三種反演方法都可以找到一個合適的結(jié)果，但在個別參數(shù)上存在一些誤差。DS反演中最佳模型參數(shù)和真值相差最大不到10.09% ，擬合差為40.57 ，反演時間相對較短；SA反演中最佳模型參數(shù)和真值相差最大不到9.89% ，擬合差為39.36 ，反演時間相對較長；ASSA反演中最佳模型參數(shù)和真值相差最大不到6.95% ，擬合差為30.17，反演時間相對適中?？梢钥闯?，綜合考慮反演速度和擬合精度，自適應(yīng)純形模擬退火法相對于下降純形法和模擬退火法具有一定的優(yōu)越性。作者著重討論自適應(yīng)純形模擬退火反演算法。

表2 對H型理論曲線反演結(jié)果

5.2 ASSA反演(加2%高斯噪聲)

分析研究自適應(yīng)純形模擬退火法，我們對表1的2個三層模型和1個五層模型的正演響應(yīng)(視電阻率和相位)分別加入2%的高斯噪聲后再反演。H、K、HKH模型分別經(jīng)過11 189、10 164、10 506次迭代，最終的擬合差依次為：30.17、36.19、32.74。反演結(jié)果如表3所示，反演模型基本反映了真實模型結(jié)構(gòu)，最佳模型參數(shù)和真值相差最大不到13%。從表3可以看出，對于三層模型，反演后結(jié)果顯示，K型的中間層分辨率不是很高，H型模型反演結(jié)果與真實模型基本吻合。五層模型的反演結(jié)果顯示第三層的深度和電阻率反演效果相對略差，但基本反映出該高阻異常體，反演結(jié)果顯示出地電結(jié)構(gòu)的變化。圖2(a) ～圖2(c)分別對應(yīng)H、K、HKH型地層反演后模型響應(yīng)和合成數(shù)據(jù)。從圖2可以看出，反演后的響應(yīng)曲線與合成數(shù)據(jù)擬合較好，阻抗相位的擬合度比視電阻率的擬合度略差。在反演中，為了檢驗該方法的全局性，我們采用了不同的隨機(jī)初始模型，進(jìn)行多次的反演計算，所有反演計算結(jié)果都收斂于同一模型的最優(yōu)解，整個算法在計算機(jī)上的運行時間為60 s左右，時間很短，反演效率比較高。

表3 加入2%高斯噪聲反演結(jié)果

5.3 ASSA反演(加5%高斯噪聲)

對加入2%高斯噪聲電磁數(shù)據(jù)的反演效果良好，我們將噪聲干擾增加到5%的高斯誤差。H、K、HKH模型分別經(jīng)過9 514、8 290、8 559次迭代，最終的擬合差依次為：37.52、43.81、40.49。反演模型參數(shù)如表4所示，反演模型參數(shù)值和真值相差不大，最佳模型參數(shù)和真值相差最大不到23%。圖3(a)～圖3(c)為合成數(shù)據(jù)曲線與反演的響應(yīng)曲線的擬合對比圖，擬合效果相對于加2%高斯噪聲的情況略差，但整體的數(shù)據(jù)變化趨勢基本一致。對于K型

圖2 合成曲線(加2% 高斯噪聲)與反演的響應(yīng)曲線對比圖Fig.2 Synthetic(2% gaussian noise) and inversion responses(a)H型；(b)K型；(c)HKH型

參數(shù)ρ1ρ2ρ3ρ4ρ5h1h2h3h4搜索范圍[11000][11000][11000][11000][11000][110000][110000][110000][110000]H型107.0710.499.01//1014.883897.84//K型21.05103.5220.41//1045.654573.40//HKH型253.1722.1093.4611.2826.67614.911087.214033.855875.55

模型的中間頻率段擬合相對于低頻和高頻段比較差，不過電性結(jié)構(gòu)還是很清楚的體現(xiàn)出來了；對于HKH型模型，低頻部分(深部)的擬合情況相對比較差，其原因為從地下深部返回的地電信息相對較弱，當(dāng)觀測資料含有與高頻相同水平的噪聲時，使得這些信息的信噪比大大降低，造成深部的反演結(jié)果變差。但是整體的擬合情況還是很好的，說明了該反演方法對于噪聲有一定的抑制作用，所有的反演都基本搜索到了模型的全局最優(yōu)解。此外，隨著噪聲的增加，收斂時擬合值增大，但是不影響該反演方法的全局收斂性。

圖3 合成曲線(加5% 高斯噪聲)與反演的響應(yīng)曲線對比圖Fig.3 Synthetic(5% gaussian noise) and inversion responses(a)H型；(b)K型；(c)HKH型

6 結(jié)論

作者采用了模擬退火法和下降純形法的非線性綜合優(yōu)化方法，綜合二者的優(yōu)點，實現(xiàn)快速有效地搜索全局最優(yōu)解。為了減小大地電磁測深反演的多解性和提高分辨率，進(jìn)行了視電阻率和阻抗相位的同時反演。

首先對比了下降純形法、模擬退火法和自適應(yīng)純形模擬退火法，突出了自適應(yīng)純形模擬退火法的反演速度和擬合精度的優(yōu)越性。然后采用自適應(yīng)純形模擬退火反演方法對三層H、K和五層HKH等模型進(jìn)行了數(shù)值試驗，檢驗了不同隨機(jī)初始模型對反演結(jié)果的影響，我們對每個模型的合成數(shù)據(jù)進(jìn)行多次的反演計算，所有的反演計算最終都收斂至同一最優(yōu)解，不受初始模型選擇的影響；加入2%高斯噪聲的情況下，H型頻段中間頻段和HKH型的中間頻段反演效果略差，其他各層電性結(jié)構(gòu)得到較好的恢復(fù)；在加入5%高斯噪聲的情況下，同樣H型頻段中間頻段和HKH型的中間頻段的反演結(jié)果相對較差，說明數(shù)據(jù)對這些結(jié)構(gòu)反映不敏感，其它各層結(jié)構(gòu)都得到較好的恢復(fù)。正演曲線和反演的響應(yīng)曲線基本吻合，反演結(jié)果基本反映模型的地電結(jié)構(gòu)變化。

數(shù)值試驗過程中，反演后模型的分辨率隨著噪聲強(qiáng)度的增大都有所降低，其中H和HKH型的分辨率降低地比較明顯，H型中頻段和HKH型中頻段的分辨率隨著噪聲的增強(qiáng)逐步降低，這是視電阻率和相位對中頻段高阻層變化的不靈敏，這也是物理規(guī)律本身決定的。

這里所利用的反演方法是綜合利用局部和全局最優(yōu)的搜索方法，采用阻抗相位和視電阻率的同時反演，大大降低了解的非唯一性，提高了反演的精度和可靠性，是一種值得推廣的混合優(yōu)化算法，同時也為后續(xù)貝葉斯反演理論的研究奠定了良好的基礎(chǔ)。

[1] 孫正江,王華俊.地電概論[M].北京:地震出版社,1984.

SUN Z J,WANG H J. Overview of geoel- ectric[M].Beijing: Earthquake Press,1984.(In Chinese)

[2] 柳建新.大地電磁測深法勘探-資料處理、反演與解釋[M].北京:科學(xué)出版社，2011.

LIU J X.Exploration of magnetotelluric sounding method[M].Beijing:Science Press,2011.(In Chinese)

[3] 陳樂壽.大地電磁測深資料處理與解釋[M].北京:石油工業(yè)出版社,1989.

CHEN L S. Magnetotelluric sounding data processing and interpretation[M].Beijing: Oil industry Press,1989. (In Chinese)

[4] 王家映.地球物理反演理論[M].北京:高教出版社,2007.

WANG J Y.Geophysical inversion theory [M]. Beijing:HighereducationPress,2007.(In Chinese)

[5] 劉云峰,曹春蕾.一維大地電磁測深的遺傳算法反演[J].浙江大學(xué)學(xué)報，1997,31(3):300-305.

LIU Y F,CAO C L. Genetic algorithm inversion of one-dimensional magnetote- lluric sounding[J].Journal of zhejiang university.1997,31(3):300-305. (In Chinese)

[6] 楊文釆.非線性地球物理反演方法回顧 與展望[J].地球物理學(xué)進(jìn)展,2001，17(2):255-261.

YANG W C. Retrospect and prospect non-linear geophysical inversion methods[J]. Earth physics reviewed,2001，17(2): 255-261. (In Chinese)

[7] 姚姚.蒙特卡洛非線性反演方法及應(yīng)用[M].北京:冶金工業(yè)出版社,1997.

YAO Y. The monte carlo nonlinear invers-ion method and its application[M].Beijing：Metallurgical industry Press,1997.(In Chinese)

[8] KIRKPATRICK S,GELATT JR CD,VECCHI M P.Optimization by simulated annealing[J]. Science,1983,220:671-680.

[9] METROPOLIS N,ROSENBLUTH A,R OSENBLUTH M,et.al. Equation of state calculations by fast computing machines[J].J .Chem Phys，1953,21:1087-1092 .

[10]KIR KPATRICK S,CELATTCD,VECCHI M P.Optimization by simulated annealing[J].Sciences,1983,220:671-680.

[11]ROTHMAN D H .Nonlinear inversion,stati-stical mechanics,and residual statics esti-mation[J]. Geophysics,1985,50(12):2784 -2796.

[12]ROTHMAN DH.Automatic estimation of large residual statics correction[J]. Geop-hysics ,1986,51(2):337-346.

[13]SEN M K , STOFFA P L .Nonlinear one-dimensional seismic wave form inversion using simulated annealing[J].Geophysics,1991,56(10):1624-1638.

[14]SEND M K ,BHATTACHARYA B B,STOFFA P L. Nonlinear inversion of resistivity sounding data[J]. Geophysics, 1993,58(4): 496-507 .

[15]師學(xué)明,王家映.一維層狀介質(zhì)大地電磁模擬退火反演法[J].地球科學(xué)-中國地質(zhì)大學(xué)學(xué)報,1998,23(5):543-545.

SHI X M,WANG J Y. one-dimensional layered media magnetotelluric simulated annealing inversion method[J]. Journal of China university of geosciences,1998,23(5): 543-545. (In Chinese)

[16]于鵬,王家林,吳健生,等.重力與地震資料的模擬退火約束聯(lián)合反演[J].地球物理學(xué)報, 2007, 50(2):529-538.

YU P,WANG J L,WU J S,et,al. Gravity and simulated annealing constrained joint inversion of seismic data[J]. J. geophys，2007, 50(2):529-538. (In Chinese)

[17]劉鵬程,紀(jì)晨.改進(jìn)的模擬退火法-單純形綜合反演方法[J].地球物理學(xué),1995, 38(2): 199-205.

LIU P C.JI C. Improved simulated annea-ling method and simplex inversion meth-od[J].Geophys，1995,38(2):199-205. (InChinese)

[18]HOLLAND,J.H.Adaptation in natural and artificial systems[M].University of Michigan Press,Ann Horbor, Michigan,1975.

[19]PRESS, W. H., TEUKOLSKY, S.A..Simulated annealing optimization over continuous spaces[J]. Comput Phys, 1991(5): 426-450.

[20]NELDER, J. A., MEAD, R. A simplex method for function minimization[J]. computer Journal,1965(7):308-313.

[21]DOSSO SE,WILMUT M J,LAPINSKI A L S.An adaptive-hybrid algorithm for geoacoustic inversion[J]. IEEE Journal of Oceanic Engineering,2001, 26(3): 324-336.

[22]劉文劼.頻率和空間相關(guān)性觀測數(shù)據(jù)條件下的大地電磁貝葉斯反演[D].長沙:中南大學(xué).2012.

LIU W J. Frequency and spatial correla-tion under the condition of the earth obs-ervation data and electromagnetic bayesi-an inversion[D]. Changsha:Central south university,2012. (In Chinese)

[23]康立山.非數(shù)值并行算法(第一冊)模擬 退火算法[M].北京:科學(xué)出版社,1997.

KANG L S. Numerical parallel algorithm of simulated annealing algorithm[M]. Beijing: Science Press,1997. (In Chinese)

[24]師學(xué)明,王家映.一維層狀介質(zhì)大地電磁模擬退火反演法[J].地球科學(xué)-中國地質(zhì)大學(xué)學(xué)報,1998,23(5):543-545.

SHI X M,WANG J Y. One-dimensional layered media magnetotelluric simulated annealing inversion method[J].Journal of China university of geosciences,1998, 23(5):543-545. (In Chinese)

[25]師學(xué)明,王家映.地球物理資料非線性反演方法講座(三)模擬退火法[J].工程地球物理學(xué)報,2007,4(3):165-174.

SHI X M,WANG J Y. Geophysical data nonlinear inversion method lecture - the simulated annealing method[J].Journal of engineering geophysics,2007, 4(3): 165-174. (In Chinese)

[26]郭榮文.貝葉斯MT反演的非線性和不確定度分析[D].長沙:中南大學(xué), 2011.

GUO R W. Bayesian MT inversion of the nonlinear and uncertainty analysis[D].Chang sha:Central south university,2011. (In Chinese)

[27]范鳴玉,張瑩. 最優(yōu)化技術(shù)基礎(chǔ)[M].北京:清華大學(xué)出版社,1982.

FAN M Y,ZHANG Y. Optimization technology base[M].Beijing: Tsinghua university Press,1982.

[28]席自強(qiáng).單純形-模擬退火法[J].湖北工學(xué)院學(xué)報,2000,15(1):27-29.

XI Z Q. Simplex simulated annealing method[J].Journal of hubei institute of technology,2000,15(1):27-29.(In Chinese)

[29]劉麗,蘇敏. 基于改進(jìn)模擬退火算法的醫(yī)學(xué)圖像配準(zhǔn)[J].計算機(jī)應(yīng)用及軟件, 2008,25(5):223-259.

LIU L,SU M. Based on improved simu-lated annealing algorithm of medical imageregistration[J].Computer applicat-ions and software,2008,25(5):223-259.(In Chinese)

[30]曹治國,汪勇.基于模擬退火-單純形法的目標(biāo)函數(shù)的優(yōu)化[J].華中科技大學(xué)學(xué)報: 自然科學(xué)版, 2005,33(6):67-68.

CAO Z G,WANG Y. Based on simulated annealing, the objective function of opti-mization of simplex method [J] .Journal of huazhong universityof science and technology,2005,33(6):67-68. (In Chinese)

One dimensional magnetotelluric sounding apparent resistivity and phase inversion of adaptive simplex simulated annealing study

SUN Huan-le1a,b, WANG Shi-biao2, GUO Rong-wen1a,b，ZHOU Shao-min1a,b, LIU Jian-xin1a,b*

(1.Central South University a.School of Geosciences and Info-Physics, b.Key Laboratory of Non-ferrous Resources and Geological Detection,Ministry of Hunan province,Changsha 410083,China; 2.Guang Dong Nonferrous Metals Engineering Investigation Design Institute,Guangzhou 510080,China)

In order to avoid the linear inversion of MT falling into local extremum,and slow convergence rate of global optimization method, the authors adopt a comprehensive optimization method of adaptive simplex simulated annealing to the inversion of MT.On the one hand,this method could combine the advantages of both downhill simplex and simulate anneal arithmetic,which is proved to be a good way with global search capability and high-speed convergence.On the other hand, it achieved the simultaneously inversion of apparent resistivity and phase, which leads to the result of reducing the multi solutions and increasing resolution. Through numerical calculation of H,K and HKH models,this method shows high search efficiency and global convergence of comprehensive optimization method.The structure characteristics of real model could be reflected by inversion result of all the model synthetic data.As for some intermediate layers are relative high resistance models,the structures still could recover after inversion even if they are relative bad to other layers.

magnetotelluric; apparent resistivity; phase; downhill simplex; simulated annealing; adaptive simplex simulated annealing

2015-08-11 改回日期：2015-09-24

國家科技基礎(chǔ)專項(2013FY110800)；國家自然科學(xué)基金(41174103)；青年自科基金(41204081)

孫歡樂(1990-)，男，碩士，從事電磁法正反演研究與應(yīng)用， E-mail：miaomis@163.com。

*通信作者：柳建新(1962-),男，教授，博士生導(dǎo)師，從事電磁法理論與應(yīng)用研究，E-mail：ljx6666@126.com。

1001-1749(2016)05-0584-09

P 631.4

A

10.3969/j.issn.1001-1749.2016.05.02