司元雷 鄭龍金

1(中國礦業(yè)大學資源與地球科學學院 江蘇 徐州 221116) 2(江蘇建筑職業(yè)技術(shù)學院信電工程學院 江蘇 徐州 221116) 3(江西有色地質(zhì)勘查五隊 江西 九江 332000)

0 引 言

瞬變電磁法(簡稱TEM)屬于時間域電磁法，它采用不接地的回線或接地線源向地層發(fā)射一次脈沖磁場[1](也稱一次場)，在發(fā)射波形關(guān)斷的期間，通過線圈或者感應探頭來接收地下目標體被激發(fā)而產(chǎn)生的二次磁場，以此來解決與之相關(guān)的地球物理探測問題[2]。TEM自20世紀被提出后，經(jīng)過幾十年的發(fā)展，目前已在多個勘測領(lǐng)域得到了廣泛的應用，例如能源開發(fā)、礦產(chǎn)勘探、水文地質(zhì)、工程物探、考古調(diào)查、軍事和環(huán)境保護等領(lǐng)域均有所涉及[3]。

目前，瞬變電磁數(shù)據(jù)的高維精準反演是擺在廣大研究人員面前的一道大難題。盡管TEM的高維反演理論研究工作早已展開，但由于高維正反演的數(shù)值計算量非常之巨大，因此在這一方面一直未曾有較大的突破。正演是反演的基礎(chǔ)，瞬變電磁高維反演還是更多地停留在理論研究階段，大規(guī)模應用于工業(yè)生產(chǎn)還有一段很長的路要走，瞬變電磁法在目前階段的反演解釋主要還是使用較為實用的一維反演方法[4]。實測資料的反演一直是TEM的核心問題之一。為此，許多地球物理學家也提出了非常多的反演方法。例如煙圈反演法、浮動薄板解釋法、人機對話自動反演法、成像類反演算法、馬奎特法和廣義逆反演法等，而且這些方法在國內(nèi)外也都得到了不同程度的應用。然而，傳統(tǒng)的反演方法大多基于最小方差原理[5]，其反演結(jié)果依賴于初始模型的選擇，容易陷于局部極值，并且巨大的矩陣運算量在反演過程中需要耗費很多時間，但是反演結(jié)果卻仍然不穩(wěn)定。從這些反演方法的實質(zhì)上來講，它們都是將非線性的問題進行線性化來處理[5]，并非真正意義上的非線性方法。因此，采用完全非線性方法反演瞬變電磁數(shù)據(jù)勢在必行，對于提升反演結(jié)果的穩(wěn)定性和精度也同樣有著非常重要的實際意義[6]。

鑒于大多數(shù)地球物理反演問題具有高度的非線性，因此，自非線性反演提出以來，便引起了相關(guān)研究學者的高度重視。它們通過各種方式非間接地求解非線性問題，將數(shù)據(jù)空間映射到模型空間，進而反演求解問題被歸化為求一個泛函極值問題的解，并采用各種優(yōu)化方法來解決對應的問題，從而在求解非線性地球物理反演問題的過程中能夠快速地得到與之對應的模型參數(shù)[1]，最終更加精準快速地解譯野外實測物探資料。凡事都具有兩面性，非線性反演方法也有不足之處。它雖然可以在全局中自動尋找最優(yōu)化的解，但是對算法參數(shù)的要求卻較為嚴格，一旦相關(guān)參數(shù)選取不合適或?qū)偈狗囱莩霈F(xiàn)早熟，或出現(xiàn)無法收斂的情況[1]。目前，現(xiàn)有的非線性反演算法和線性反演算法都依賴于初始模型的設(shè)置，若設(shè)置不得當，很有可能得到的并不是最優(yōu)解。除此之外，也存在著容易陷入局部極值、收斂速度緩慢等一系列問題[7]。例如：遺傳算法(Genetic Algorithms，GA)可以從概率的角度隨機找到最優(yōu)解，且把握全局搜索過程的能力很強[8]，但在實際運用中也存在一些問題，例如早熟現(xiàn)象和較差的局部優(yōu)化性能；模擬退火(Simulated Annealing，SA)算法在局部搜索方面具有優(yōu)勢，可以避免在搜索過程當中落入局部最優(yōu)解，但它在把握全局搜索空間方面力不從心，搜索過程不能進入最有希望的尋優(yōu)區(qū)域，運行效率低下[6]。本文首先通過測試函數(shù)對遺傳退火進化算法(Genetic Annealing Evolutionary Algorithm，GAEA)進行數(shù)值實驗，然后將其應用于瞬變電磁地電模型理論數(shù)據(jù)的反演，從而驗證反演算法的有效性，以期提升瞬變電磁反演解釋工作的準確性。

1 遺傳退火進化算法

GA是基于自然遺傳學進化思想的啟發(fā)式搜索算法，它基于對隨機搜索的智能化開發(fā)，將其轉(zhuǎn)移到獲得問題最優(yōu)解的方向。遺傳算法的概念是約翰·霍蘭德(John Holland)提出的，它是一種解決優(yōu)化問題并克服傳統(tǒng)優(yōu)化方法所帶來的局限性的方法。搜索過程從一組稱為種群的染色體開始，每個染色體表示為問題解的編碼。每一代種群通過使用遺傳算子的重組過程隨著時間的推移而進化，并朝著更好的個體進化以產(chǎn)生新種群。為了量化個體的適應度，使用了適應度函數(shù)，基于適者生存原理的自然選擇和遺傳操作生成新一代個體。GA是一個迭代過程，需要按周期進行操作，直到滿足收斂或停止標準為止[7]。SA算法被認為是一種模擬金屬退火過程的迭代隨機搜索算法。SA算法的基本思想是搜索當前個體的鄰居以找到新的個體，它具有避免陷入局部最優(yōu)狀態(tài)的能力，這是因為SA算法在每次迭代中都會將新個體與當前個體進行比較。如果新的個體更好(基于其適應性值判斷)，那么它將被選作下一次迭代的基礎(chǔ)；否則，根據(jù)收受概率，SA算法可能會采用此個體[9-10]。

GAEA是一種混合遺傳模擬退火算法[11]，它通過組合GA和SA算法的優(yōu)勢來求解最優(yōu)化問題[12-22]。GA是此混合算法的主要框架，而SA算法被用作局部搜索策略，以幫助GA跳出局部最優(yōu)值。本文中的遺傳退火進化算法主要參考邢文訓等[10]及何則干等[16]給出的方法，這里僅給出該算法的流程圖，如圖1所示，其余部分將不再贅述。對于瞬變電磁測深的反演建模，首先選擇一組(或總體)參數(shù)模型，并定義了模型的電阻率和層厚度等參數(shù)，然后設(shè)置每個模型參數(shù)的搜索范圍[23]：

(1)

式中：mk,j是參數(shù)向量；n是參數(shù)個數(shù)。在本文使用的遺傳算法中，模型的所有參數(shù)都以二進制形式進行編碼。

2 具體實現(xiàn)

2.1 使用經(jīng)典測試函數(shù)進行對比

在進行瞬變電磁非線性反演之前，為了檢驗這種GAEA的有效性，這里首先選用兩個經(jīng)典的測試函數(shù)：Schaffer函數(shù)及Alpine函數(shù)。其中Schaffer函數(shù)是一種含有無數(shù)個極小值點的多維函數(shù)，由于其分布呈現(xiàn)強烈震蕩形態(tài)，因而尋找全局最優(yōu)值的難度非常大；而Alpine函數(shù)是一種多峰最小化測試函數(shù)。當函數(shù)趨于無窮大時，該函數(shù)將沿著自變量的方向生成大量不同的局部極值，這給優(yōu)化帶來了很大的困難。因此，可將GA、SA和GAEA三種算法進行比較。兩個測試函數(shù)分別為：

-100≤xi≤100i=1,2

(2)

-10≤xi≤10i=1,2,…,n

(3)

式中：Schaffer函數(shù)在其定義域內(nèi)只有一個全局最小點f(0,0)=0；x*為當Alpine函數(shù)取得全局最優(yōu)值時xi的值，當x*=0時，Alpine函數(shù)取得全局最優(yōu)值f(x*)=0。函數(shù)圖形如圖2所示。

(a) Schaffer函數(shù)

(b) Alpine函數(shù)圖2 兩個測試函數(shù)的立體圖

2.1.1測試函數(shù)一

將種群設(shè)置為60，精度設(shè)置為1e-20，使用GA、SA和GAEA三種算法分別求解Schaffer函數(shù)的最小值點f(0,0)=0，圖3是通過三種算法計算所得到的每一代種群中最小目標值的尋優(yōu)過程。可以看出GAEA精度最高，SA算法精度最低。另外，相對于GA，GAEA不僅尋優(yōu)能力更強且收斂速度更快。

圖3 三種算法在Schaffer函數(shù)尋優(yōu)過程的目標值圖

鑒于上述三種算法的搜索方式具有隨機性，為了更客觀地評價算法的性能，這里將算法的容許誤差設(shè)置為1e-20，三種算法獨立地運行60次，得到的尋優(yōu)結(jié)果如表1所示。

表1 三種算法的計算結(jié)果對比

通過對比表1中三種算法的計算結(jié)果，GAEA相對GA而言精度更高，GAEA得到的最優(yōu)解與Schaffer函數(shù)的全局最小值點最相近，遠比SA算法優(yōu)越。

2.1.2測試函數(shù)二

將種群設(shè)置為60，精度設(shè)置為1e-20，使用GA、SA和GAEA三種算法分別求解Alpine函數(shù)的最小值點f(x*)=0，圖4是通過三種算法計算所得到的每一代種群中最小目標值的尋優(yōu)過程?？梢钥闯鯣AEA精度最高，SA算法精度最低。另外，相對于GA，GAEA不僅尋優(yōu)能力更強且收斂速度更快。

圖4 三種算法在Alpine函數(shù)尋優(yōu)過程的目標值圖

同樣地，這里將算法的容許誤差設(shè)置為1e-20，三種算法獨立地運行60次，得到的尋優(yōu)結(jié)果如表2所示。

表2 三種算法的計算結(jié)果對比

從運算結(jié)果中可知，GAEA有12次得到全局最優(yōu)解0；遠遠好于SA算法和GA。綜合兩個函數(shù)的測試結(jié)果可知，GA和SA算法均以一定的概率逼近全局最優(yōu)解。從比較結(jié)果不難發(fā)現(xiàn)，GAEA逼近全局最優(yōu)解的概率較大，落入局部極值的概率較小。

2.2 層狀地電模型反演

2.2.1反演目標函數(shù)的確定

通過反演算法對瞬變電磁數(shù)據(jù)進行的分析包括使用感應電動勢或視電阻率的值來估計被分析介質(zhì)的真實電阻率、厚度或深度值。對于一維情況，水平分層介質(zhì)被廣泛使用。該模型的自由參數(shù)是每層的電阻率ρn和厚度hn，由向量m表示，同時地電模型參數(shù)可被記為x。在這種情況下，反演問題即可被轉(zhuǎn)化為一個求解多目標非線性優(yōu)化問題。因此，構(gòu)建目標函數(shù)S(x)，層狀地電模型的最優(yōu)解xopt需滿足S(xopt)=min(S(x))，即：

(4)

2.2.2遺傳退火進化算法反演參數(shù)的選擇

基于上述的目標函數(shù)，針對算法參數(shù)進行了詳細的試算與篩選，初始群體由隨機概率產(chǎn)生，種群的大小設(shè)置為60，最大迭代次數(shù)設(shè)置為30 000，種群中個體的交叉概率設(shè)置為0.7，每一位個體的變異概率設(shè)置為0.1，個體的選擇方法使用輪盤賭法，采用浮點法進行編碼，允許誤差設(shè)置為1e-04。模擬退火初始溫度T0為1 500℃，降溫方式T=T0×0.98k-1，模擬退火的終止步數(shù)k設(shè)置為20。另外，在模擬退火算法部分，隨機生成小的擾動。

2.2.3不同地電理論模型試算

為驗證上述反演算法在瞬變電磁資料解釋中的可行性，使用不同瞬變電磁理論模型進行了反演試算。其中瞬變電磁正演計算部分采用文獻[24-27]中的方法。反演過程中先建立n層介質(zhì)模型，其參數(shù)為λ=(ρ1,ρ2,…,ρn,h1,h2,…,hn-1)，ρ和h分別為各層介質(zhì)的電阻率和厚度；使用式(4)構(gòu)建反演算法的目標函數(shù)，設(shè)定初始模型后進行迭代計算直至達到擬合終止條件。計算中，發(fā)射線圈和接收線圈的邊長分別設(shè)置為100 m和1 m，匝數(shù)均為1，供電電流為10 A。為了與傳統(tǒng)線性方法進行比較，除采用GA、SA算法、GAEA方法外，還采用傳統(tǒng)線性化反演方法的代表馬奎特方法(Marquardt)進行了反演計算[28]，每種算法獨立運行20次，然后取20次反演結(jié)果的均值作為模型參數(shù)的估計值。其中：GA的種群的大小設(shè)置為60，最大迭代次數(shù)為30 000，個體的交叉概率為0.7，變異概率為0.1，選擇方法使用輪盤賭法，采用浮點法編碼，允許誤差設(shè)置為1e-04。在交叉與變異運算中，進行交叉運算的個體和交叉的位置及變異位，均通過產(chǎn)生的隨機數(shù)與交叉概率與變異概率的大小比較來判定。SA算法的模型群體數(shù)量為60，最大迭代次數(shù)為30 000，初始溫度T0為1 500 ℃，降溫方式T=T0×0.98k-1，模擬退火的終止步數(shù)k設(shè)置為20，隨機產(chǎn)生小的擾動。下面將分別針對選取的三層低阻模型及五層地電模型進行反演試算，并分析相應的反演結(jié)果。

(1) 三層H型地電模型。表3為三層H型地電模型的參數(shù)及多種算法的反演結(jié)果，同時，GAEA的反演結(jié)果經(jīng)過正演計算得到感應電位之后，與真實模型的正演響應進行了對比(如圖5所示)，二者的感應電動勢曲線基本重合，證明反演結(jié)果是有效的。另外，將表3中的反演結(jié)果繪制成柱狀圖可以更加直觀地看出每種算法反演的效果，如圖6所示，可以看出GAEA反演效果最好，其各參數(shù)反演結(jié)果的誤差均控制在2%以內(nèi)，可以滿足實際需要。

表3 三層H型地電模型反演結(jié)果

圖5 原始模型響應和GAEA反演結(jié)果響應的對比

圖6 反演結(jié)果柱狀圖對比

(2) 五層HKH型地電模型。表4為五層HKH型地電模型的參數(shù)及多種算法的反演結(jié)果，GAEA的反演結(jié)果經(jīng)過正演計算得到感應電位之后，與真實模型的正演響應進行了對比，如圖7所示，二者的感應電動勢曲線基本重合，證明反演結(jié)果是有效的。從圖8中的直觀結(jié)果易知，與Marquardt法、GA、SA算法相對比，GAEA在反演TEM數(shù)據(jù)方面更為準確和有效。因為五層模型參數(shù)相較于三層更多，所以對應的反演難度加大，從而導致了結(jié)果的誤差也相應增大，但是該模型的關(guān)鍵信息仍可以由GAEA獲取，其大部分參數(shù)的反演結(jié)果的誤差控制在6%以內(nèi)，僅個別參數(shù)誤差超過10%(第三層電阻率參數(shù)ρ3的誤差為12.9%)，這在瞬變電磁法勘探的實際生產(chǎn)中，仍能滿足實際需要。

表4 五層HKH型地電模型反演結(jié)果

圖7 原始模型響應和GAEA反演結(jié)果響應的對比

圖8 反演結(jié)果柱狀圖對比

3 結(jié) 語

本文將一種新的GAEA應用于瞬變電磁法的反演中，與傳統(tǒng)線性反演方法的對比實驗驗證了GAEA的有效性和可行性。主要結(jié)論如下：

(1) 經(jīng)融合后得到的遺傳退火進化算法，結(jié)合了GA的優(yōu)秀全局搜索性能和SA算法的局部搜索能力強的特點。與此同時，種群在模擬退火操作中還能充分利用融合算法所帶來的全局信息。

(2) GAEA方法在1D合成數(shù)據(jù)中的實現(xiàn)和應用對于更好地理解這些方法的優(yōu)缺點具有重要意義。合成數(shù)據(jù)的反演結(jié)果表明，基于遺傳退火進化算法的瞬變電磁非線性反演比線性反演更有效，為2D反演算法的發(fā)展奠定了基礎(chǔ)。

(3) 將GAEA應用于兩個測試函數(shù)，即Schaffer函數(shù)及Alpine函數(shù)，該技術(shù)在尋找測試函數(shù)的最優(yōu)極值方面具有非常好的效果；將GAEA應用于瞬變電磁地電模型的合成數(shù)據(jù)，模型的關(guān)鍵信息仍可以由GAEA獲取，但是需要大量的處理時間。這驗證了該算法的適用性，其缺點可以在將來的實施中加以改進。