周輝

（青島海洋大學海洋地球科學學院，青島　266003）

何樵登

（長春地質(zhì)學院地球物理系，長春　130026）

摘要　各向異性介質(zhì)參數(shù)反演通常為非線性優(yōu)化問題。非線性反演方法可以分為兩大類：隨機搜索方法，如Monte Carlo法、模擬退火和遺傳算法及基于非線性最小平方理論的非線性解析反演方法。遺傳算法能尋找到全局最優(yōu)解，但它為一種較費時的方法。非線性解析反演方法能給出一個與初始模型有關(guān)的局部最優(yōu)解。然而，這種方法具有較快的收斂速度。遺傳算法與非線性解析反演方法相結(jié)合的反演方法利用這兩種反演方法的優(yōu)點而克服其缺點。因此，結(jié)合反演方法既能快速收斂，又能尋找到全局最優(yōu)解。如何合理地將遺傳算法和非線性解析反演方法結(jié)合是十分重要的。本文提出一種結(jié)合方案，即在連續(xù)若干次遺傳算法迭代后作一次非線性解析反演。理論算例表明結(jié)合反演方法具有上述特點。

關(guān)鍵詞　遺傳算法　非線性解析反演　非線性結(jié)合反演　各向異性介質(zhì)

1　引言

遺傳算法為隨機搜索類方法之一，它以概率論為理論基礎(chǔ)，用于求解多極值復雜優(yōu)化問題[9]。遺傳算法不要求已知模型空間中后驗概率密度的形狀并能廣泛搜索模型空間。遺傳算法模擬自然選擇和遺傳規(guī)律，并遵循適者生存的原則。

遺傳算法由Holland在1975年提出[4]。Berg首先將遺傳算法應(yīng)用于地球物理優(yōu)化問題[1]。Stoffa等系統(tǒng)地研究了種群大小、交叉概率、選擇概率和變異概率對多參數(shù)優(yōu)化問題收斂性和收斂速度的影響[11]。Sen等討論了在選擇概率中引入溫度參數(shù)的作用并提出一些退火方案[10]。周輝等則研究了目標函數(shù)與收斂速度和解的精度的關(guān)系[16]。

基于最小平方優(yōu)化理論的非線性反演方法是兩大類反演方法之一。當給定的初始模型位于目標函數(shù)全局最優(yōu)解所在的峰谷附近時，這種下降類方法能給出正確解而與初始模型位置無關(guān)。下降類算法研究得較深入，應(yīng)用較廣。

Tarantola提出一種基于廣義最小二乘法的多維多偏移距聲波地震波形解釋的一般性非線性地震波形反演方法[12]。隨后，Tarantola將該理論推廣于各向同性介質(zhì)的彈性波反演[13]。Gauthier等用理論數(shù)據(jù)驗證了Tarantola提出的方法的正確性[2]。稍后，Tarantola研究非線性解析法反射波彈性反演的策略，指出以縱橫波的波阻抗和密度作為反演參數(shù)，才盡可能使反演參數(shù)之間相互獨立[14]。Pan用τ—P變換研究層狀聲學介質(zhì)中平面波地震記錄非線性解析反演的理論和可行性[6]。為了更多地利用地震數(shù)據(jù)中的信息，包括VSP資料中反射和轉(zhuǎn)換信息，Mora作了一些工作[5]。當僅用反射數(shù)據(jù)時反演主要解決引起反射的P波和S波的波阻抗突變。當利用轉(zhuǎn)換數(shù)據(jù)時，則能分辨大尺度的P波和S波速度變化。Sambridge等改進了修改模型的方法[8]。在子空間中，可同時得到P波、S波波阻抗和密度。周輝等將非線性梯度反演方法推廣于多維、多道、多分量任意彈性各向異性介質(zhì)參數(shù)的反演[17]。

非線性解析反演方法和遺傳算法結(jié)合的反演方法利用非線性解析反演和遺傳算法的優(yōu)點，克服它們的缺點。因此，結(jié)合反演方法不僅能搜索到全局最優(yōu)解，而且能較快地收斂。Porsani等在遺傳算法和廣義線性反演方法相結(jié)合方面作了一些研究[7]。

本文討論各向異性介質(zhì)的非線性解析反演方法和遺傳算法與非線性解析反演方法相結(jié)合的結(jié)合反演方法[17]。對于遺傳算法讀者可參考遺傳算法的相關(guān)文獻[3,9～11]。

2　各向異性介質(zhì)參數(shù)非線性解析反演方法

2.1　共軛梯度法

反演的目的是利用地面或井中測得的位移場ui（xr,t）求取地下介質(zhì)密度分布ρ（x）和彈性參數(shù)分布Cijkl（x）。ρ（x）、Cijkl（x）稱為模型參數(shù)。x為研究介質(zhì)中或邊界上任一點，x=（x1,x2,x3）,xr為接收點。反演的目標是使目標函數(shù)

巖石圈構(gòu)造和深部作用

取極小值。其中Cd、Cm分別為數(shù)據(jù)（波場）和模型參數(shù)的協(xié)方差算子。m0為先驗?zāi)Ｐ蛥?shù)，m為反演過程中求得的模型參數(shù)。由于模型參數(shù)有多個，故用向量表示。ucal為給定m的波動方程正演記錄，uobs為觀測波場，上角標t表示轉(zhuǎn)置。地震記錄u和模型參數(shù)m之間的函數(shù)關(guān)系為

巖石圈構(gòu)造和深部作用

g為非線性算子，（2）式為波動方程的算子形式。記第n次迭代時的模型參數(shù)為mn，則有

巖石圈構(gòu)造和深部作用

及共軛梯度法的迭代公式[15]

巖石圈構(gòu)造和深部作用

其中Gn為g對mn的FrEChet導數(shù)，ηn為一常數(shù)，可由多種方法計算[5,8]。

梯度

為模型空間的對偶空間中的一個元素。模型空間和其對偶空間以模型參數(shù)的協(xié)方差算子Cm=Diag（Cp，Cc）由式（4d）相聯(lián)系。在后面將給出

和

的表達式。

式（4）為梯度反演方法的基本公式。當該公式中的每一量都已知時，迭代就可進行。在這些變量中，最關(guān)鍵的是梯度向量。

2.2　目標函數(shù)

在最小二乘理論中，權(quán)函數(shù)是協(xié)方差算子逆的核。假設(shè)數(shù)據(jù)集中的誤差是不相關(guān)的，它僅取決于時間或源和接收器的位置，那么有[14]

巖石圈構(gòu)造和深部作用

其中σ為數(shù)據(jù)的均方差。

2.3　各向異性介質(zhì)中的彈性波動方程

令fi（x,t;xs）是第s次激發(fā)的內(nèi)體力密度，Ti（x,t;xs）是地球表面S的應(yīng)力矢量分量，ni（x）是表面的單位法向分量。那么與第s次激發(fā)相應(yīng)的位移由以下微分方程組給出[15]

巖石圈構(gòu)造和深部作用

2.4　梯度向量

式（4）中梯度向量的分量為[17]

巖石圈構(gòu)造和深部作用

其中，T為地震記錄的長度，

為反向傳播場，滿足

巖石圈構(gòu)造和深部作用

其中，t∈[T,0]，

滿足終了時間條件。

3　結(jié)合反演方法

3.1　遺傳算法和非線性解析反演方法的優(yōu)缺點

遺傳算法是利用概率論來求解多極值復雜優(yōu)化問題的一種隨機搜索方法，由一組隨機選取的模型開始，不需要更多的先驗信息，廣泛而有效地對模型空間的最優(yōu)部分采樣。盡管遺傳算法是基于自然選擇、遺傳規(guī)律，搜索模型空間的最優(yōu)部分而求得最優(yōu)解，但它是一種計算量很大的方法。由于地震模型空間大，用全局最優(yōu)化方法估計各向異性介質(zhì)參數(shù)的地震波形反演十分費時。

目標函數(shù)的梯度信息是非線性解析反演方法修改模型參數(shù)的依據(jù)，它能給出一個接近初始模型的一個局部最優(yōu)解。如果初始模型選擇得合適，即當初始模型處在全局最優(yōu)解所在的目標函數(shù)低谷時，非線性解析反演方法能收斂于全局最優(yōu)解。然而，恰好給出一個接近全局最優(yōu)解的初始模型的概率是非常小的，尤其對沒有模型參數(shù)的任何先驗信息的情況。但應(yīng)強調(diào)的是，非線性解析反演方法具有較快的收斂速度。

發(fā)揮非線性解析反演方法快速收斂和遺傳算法能搜索到全局最優(yōu)解的優(yōu)點，而克服前者僅能尋找到局部最優(yōu)解和后者運算量大的缺點是很有意義的。非線性解析反演方法和遺傳算法相結(jié)合的反演方法可達到上述目的。在結(jié)合反演方法中，遺傳算法的作用是提供接近全局最優(yōu)解的模型，非線性解析反演的作用是盡快求出全局最優(yōu)解。因此，結(jié)合反演方法具有搜索到全局最優(yōu)解的能力和比遺傳算法收斂速度快的特點。

3.2　結(jié)合方案

遺傳算法在優(yōu)化過程中連續(xù)不斷地搜索整個模型空間。在每次迭代結(jié)束后，得到一個本代的最優(yōu)模型。根據(jù)遺傳算法的數(shù)學原理[3]，最優(yōu)模型的數(shù)量在下一代中得以增加，同時經(jīng)交叉和變異作用又有新的模型產(chǎn)生。在下一代種群中，最優(yōu)模型可能與前一代的相同，也有可能劣于前一代的最優(yōu)模型。所有這些最優(yōu)模型可能在目標函數(shù)的同一低谷處，也有可能在其它低谷處。遺傳算法尋找最優(yōu)模型要經(jīng)過多次迭代才能確定一個極值。遺傳算法的隨機性導致遺傳算法是一種費時的方法。然而正是遺傳算法的這種隨機性保證了它能搜索到全局最優(yōu)解。

如果將每次遺傳算法迭代的最優(yōu)解作為非線性解析反演的初始模型，非線性解析反演可以找出與初始模型毗鄰的局部最優(yōu)解。由于非線性解析反演是一種確定性的方法，它按目標函數(shù)的梯度方向修改模型，所以非線性解析反演方法只需幾次迭代即可收斂。非線性解析反演求得的解是否為全局最優(yōu)解，非線性解析反演方法本身是無法得以保證的。只有當遺傳算法提供接近全局最優(yōu)解的初始模型時，非線性解析方法反演才能收斂到全局最優(yōu)解。

結(jié)合反演方法中遺傳算法和非線性解析反演方法的匹配方式是十分重要的。非線性解析反演方法得到接近遺傳算法提供的初始模型的局部最優(yōu)解后，在以后若干代中因遺傳算法的隨機性而使其最優(yōu)解與該局部最優(yōu)解相同。如果每次遺傳算法迭代后作非線性解析反演，那么結(jié)合反演的結(jié)果在幾代內(nèi)都是相同的。顯然其中的一些非線性解析反演是沒有必要的。因此，結(jié)合方式應(yīng)為在連續(xù)多次遺傳算法迭代后作一次非線性解析反演，然后將非線性解析反演的結(jié)果作為下一代種群中的一個母本模型。圖1為結(jié)合反演的框圖。

圖1　結(jié)合反演框圖

4　算例

為了驗證結(jié)合反演方法的優(yōu)越性，對一維多層橫向各向同性介質(zhì)參數(shù)的反演理論實例作了分析。

圖2是目標函數(shù)值與迭代次數(shù)的關(guān)系圖。在該結(jié)合反演算例中每次遺傳算法迭代后就作一次非線性解析反演迭代。結(jié)合反演的誤差在開始幾次迭代中下降很快，尤其在前3次。結(jié)合反演方法在第10次迭代達到的較小誤差，遺傳算法在第42次迭代才達到。結(jié)合反演的誤差比遺傳算法的跳躍得嚴重。這是因為非線性解析反演得到的模型在遺傳算法中作為母代參加繁衍。這個模型因遺傳算法的隨機性常常被新的模型替代。這兩個模型可能位于目標函數(shù)兩個不同的低谷中，因此非線性解析反演的結(jié)果不同。

盡管結(jié)合反演的目標函數(shù)有些振蕩，但也存在連續(xù)幾次迭代目標函數(shù)幾乎不變的現(xiàn)象。這意味著這幾次迭代的最優(yōu)模型是很接近的。在這種情況下非線性解析反演不能提供較大的改進。所以，此時的非線性解析反演是沒有必要的，否則只能增加計算量。

圖2　結(jié)合反演（實線）和遺傳算法（虛線）的誤差與迭代次數(shù)的關(guān)系

結(jié)合反演中每次遺傳算法迭代后作一次非線性解析反演迭代

圖3是另一個例子。在該結(jié)合反演例子中，每五次遺傳算法迭代作一次非線性解析反演。在這里遺傳算法占主要地位。此時結(jié)合反演的誤差函數(shù)明顯比遺傳算法的小。結(jié)合反演的誤差在第5次迭代末突然下降，并在第10次迭代時的小誤差，遺傳算法在42代才達到。遺傳算法始終沒有到達結(jié)合反演的最小誤差。結(jié)合反演的誤差在后期迭代過程中平建筑聲學測量穩(wěn)下降，這是遺傳算法占主導地位的原因。

從該例可知，若遺傳算法與非線性解析反演方法比較合理地結(jié)合，結(jié)合反演方法比遺傳算法具有快得多的收斂速度。

5　結(jié)論

非線性結(jié)合反演方法揚遺傳算法和非線性解析反演方法之長，抑其之短，它是一種具有較快收斂速度的全局反演方法。

在結(jié)合反演中遺傳算法和非線性解析反演方法的結(jié)合方式是重要的。從算例可得出，五次遺傳算法迭代后作一次非線性解析反演的結(jié)合反演的效果明顯優(yōu)于每次遺傳算法迭代后都作非線性解析反演的結(jié)合反演的效果。但是在結(jié)合反演中連續(xù)作多少次遺傳算法迭代及連續(xù)迭代次數(shù)在整個迭代過程中的可變性還有待于進一步研究。

圖3　結(jié)合反演（實線）和遺傳算法（虛線）的誤差與迭代次數(shù)的關(guān)系

結(jié)合反演中每五次遺傳算法迭代后作一次非線性解析反演迭代

在結(jié)合反演中遺傳算法的作用是提供接近全局最優(yōu)解的初始模型。結(jié)合反演的運算速度主要取決于遺傳算法的運算速度。均勻設(shè)計理論可以應(yīng)用于遺傳算法以加快隨機搜索的速度。

與遺傳算法相同，其它隨機搜索方法也可用來與非線性解析反演方法形成結(jié)合反演方法。

參考文獻

[1]E.Berg.Simple convergent genetic algorithm for inversion of multiparameter data.SEG60 Expanded Abstracts,1990，Ⅱ，1126～1128.

[2]O.Gauthier,J.Virieux and A.Tarantola.Two-dimensional nonlinear inversion of seismic waveforms:Numerical results.Geophysics,1986,51,1387～1403.

[3]D.E.Goldberg.Genetic Algorithms in Search,Optimiztion,and Machine Learning.Addison-Wesley,Reading,MA,1989.

[4]J.H.Holland.Adaptation in Natural and Artifical Systems.The University of Michigan Press,Ann Arbor,1975.

[5]P.Mora.2D elastic inversion of multi-offset seismic data.Geophysics,1988,52,2031～2050.

[6]G.S.Pan,R.A.Phinney,and R.I.ODOm.Full-waveform inversion of plane-wave seismograms in stratified acoustic media:Theory and feasibility.Geophysics,1988,53,21～31.

[7]M.J.Porsani,P.L.Stoffa,M.K.Sen,et al..A combined Genetic and linear inversion algorithm for seismic wave-form inversion.SEG63 Expanded Abstracts,1993,692～695.

[8]M.S.Sambridge,A.Tatantola and Kennet.An alternative strategy for nonlinear inversion of seismic waveforms.Geophysical Prospecting,1991,39,723～736.

[9]M.Sambridge,and G.Drijkoningen.Genetic algorithms in seismic waveform inversion.Geophys.J.Int.,1992,109,323～342.

[10]M.K.Sen,P.L.Stoffa.Rapid sampling of model space using genetic algorithms:examples from seismic waveform inversion.Geophys.J.Int.,1992,109,323～342.

[11]P.L.Stoffa,M.K.Sen.Nonlinear multiparametre optimization using genetic algorithms:Inversion of plane-wave seismograms.Geophysics,1991,56,1794～1810.

[12]A.Tarantola.Inversion of seismic reflection data in the acoustic approximation.Geophysics,1984（a）,49,1259～1266.

[13]A.Tarantola.The seismic reflection inverse problem.In:F.Santosa,Y.-H.Pao,W.W.System,and C.Holland Eds.Inverse problems of acoustic and elastic waves.Soc.Industr.Appl.Math.,1984（b）,104～181.

[14]A.Tarantola.A strategy for nonlinear elastic inversion of seismic reflection data.Geophysics,1986,51,1893～1903.

[15]A.Tarantola.Inverse problem theory:Methods for data fitting and model parameter estimation.Elsevier Science Publ.Co.Inc.,1987.

[16]周輝，何樵登.遺傳算法在各向異性介質(zhì)參數(shù)反演中的應(yīng)用.長春地質(zhì)學院學報，1995,25，增刊1,62～67.

[17]周輝.各向異性介質(zhì)波動方程正演及其非線性反演方法研究.長春地質(zhì)學院博士論文，1995.