基于遺傳退火進化算法的瞬變電磁非線性反演

2020-11-12 10:39:16司元雷鄭龍金

計算機應用與軟件 2020年11期

司元雷鄭龍金

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

0 引言

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

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

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

1 遺傳退火進化算法

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

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

(1)

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

2 具體實現

2.1 使用經典測試函數進行對比

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

-100≤xi≤100i=1,2

(2)

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

(3)

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

(a) Schaffer函數

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

2.1.1測試函數一

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

圖3 三種算法在Schaffer函數尋優過程的目標值圖

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

表1 三種算法的計算結果對比

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

2.1.2測試函數二

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

圖4 三種算法在Alpine函數尋優過程的目標值圖

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

表2 三種算法的計算結果對比

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

2.2 層狀地電模型反演

2.2.1反演目標函數的確定

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

(4)

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

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

2.2.3不同地電理論模型試算

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

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

表3 三層H型地電模型反演結果

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

圖6 反演結果柱狀圖對比

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