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

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

?

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

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

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

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

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

0 引言

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

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

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

1 一維大地電磁測(cè)深正演

1.1 正演公式

假設(shè)地下電性結(jié)構(gòu)是水平層狀分布的，如果地電結(jié)構(gòu)共有n層，則有2n-1個(gè)參數(shù)， hi(i=1,2,…,n-1)代表第i層的厚度；ρi(i=1,2,…,n)代表第i層的電阻率，最后一層為均勻半空間，則hn=∞。計(jì)算視電阻率ρ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]通過(guò)模型正演響應(yīng)值、觀測(cè)值及其數(shù)據(jù)噪聲的標(biāo)準(zhǔn)差來(lái)建立，具體表達(dá)式如式(4)所示。

(4)

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

2 模擬退火算法

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

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

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

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

3)采用Metropolis抽樣準(zhǔn)則對(duì)模型選擇。令Δ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)計(jì)算目標(biāo)函數(shù)E(Mk+1)，如果E(Mk+1)

3 下降純形法

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

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

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

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

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

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

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

重復(fù)以上四個(gè)步驟直到在純形中最高點(diǎn)函數(shù)值Ehigh和最低點(diǎn)函數(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)解，但其耗費(fèi)大量時(shí)間、收斂速度慢[28]。下降純形法能直接、快速地搜索到極小值，但當(dāng)目標(biāo)函數(shù)有不止一個(gè)極小值時(shí)，不同的初值選擇會(huì)導(dǎo)致不同的反演結(jié)果，容易導(dǎo)致局部最優(yōu)解代替全局最優(yōu)解[30]。因此，綜合二者的優(yōu)點(diǎn)，既可以提高搜索效率，又避免了陷入局部極小值。

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

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

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

2)下降純形法搜索。以M0為初始點(diǎn)，采用下降純形算法搜索該點(diǎn)周?chē)囊粋€(gè)極小值Mmin，并令Mk=Mmin，E(Mk)=E(Mmin)，MOPT=Mmin，EOPT=E(Mmin)。

3)模擬退火法。以Mk為初始點(diǎn)采用模擬退火算法，計(jì)算溫度Tk，在可行解搜索空間隨機(jī)產(chǎn)生新點(diǎn)Mk+1=Mk+rV0，其中r為隨機(jī)數(shù)，r∈(-1,1)，計(jì)算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)中，下降純形法反射、延伸和收縮等每個(gè)步驟的函數(shù)值計(jì)算之前，對(duì)該模型的各參數(shù)設(shè)置Cauchy分布的擾動(dòng)[21,26]。

(6)

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

(7)

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

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

Tk=αkT0

(8)

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

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

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

表1 模型參數(shù)

5.1 計(jì)算實(shí)例對(duì)比

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

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

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

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

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

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

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

對(duì)加入2%高斯噪聲電磁數(shù)據(jù)的反演效果良好，我們將噪聲干擾增加到5%的高斯誤差。H、K、HKH模型分別經(jīng)過(guò)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)曲線的擬合對(duì)比圖，擬合效果相對(duì)于加2%高斯噪聲的情況略差，但整體的數(shù)據(jù)變化趨勢(shì)基本一致。對(duì)于K型

圖2 合成曲線(加2% 高斯噪聲)與反演的響應(yīng)曲線對(duì)比圖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

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

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

6 結(jié)論

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

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

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

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

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

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

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

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

[3] 陳樂(lè)壽.大地電磁測(cè)深資料處理與解釋[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] 劉云峰,曹春蕾.一維大地電磁測(cè)深的遺傳算法反演[J].浙江大學(xué)學(xué)報(bào)，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é)-中國(guó)地質(zhì)大學(xué)學(xué)報(bào),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é)報(bào), 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)性觀測(cè)數(shù)據(jù)條件下的大地電磁貝葉斯反演[D].長(zhǎng)沙:中南大學(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ù)值并行算法(第一冊(cè))模擬 退火算法[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é)-中國(guó)地質(zhì)大學(xué)學(xué)報(bào),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é)報(bào),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].長(zhǎng)沙:中南大學(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é)報(bào),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ì)算機(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]曹治國(guó),汪勇.基于模擬退火-單純形法的目標(biāo)函數(shù)的優(yōu)化[J].華中科技大學(xué)學(xué)報(bào): 自然科學(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

國(guó)家科技基礎(chǔ)專(zhuān)項(xiàng)(2013FY110800)；國(guó)家自然科學(xué)基金(41174103)；青年自科基金(41204081)

孫歡樂(lè)(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