模擬退火方法在三維速度模型地震波走時反演中的應(yīng)用

李飛 張雪梅 陳宏峰 何少林

中國地震臺網(wǎng)中心，北京市西城區(qū)三里河南橫街5號 100045

0 引言

地震層析成像從醫(yī)學(xué)CT引入地球物理學(xué)領(lǐng)域用于地殼上地幔速度結(jié)構(gòu)研究（Aki et al，1976）以來，已發(fā)展成為研究地球內(nèi)部速度結(jié)構(gòu)（白志明等，2016）最重要的方法之一。利用地震層析成像得到地球內(nèi)部高分辨率的速度結(jié)構(gòu)模型，為探索全球構(gòu)造進(jìn)程，巖石圈的演化和構(gòu)造運動，全球板塊運動的規(guī)律和動力機(jī)制，以及地震、火山活動發(fā)生的深部構(gòu)造環(huán)境等提供了十分重要的科學(xué)依據(jù)（Zhao，2001；田有等，2009）。地震層析成像技術(shù)根據(jù)所使用的地震波、資料類型的不同，可分為體波層析成像（Aki et al，1977；Zhao，2004）和面波層析成像（Shapiro et al，2005）；根據(jù)研究區(qū)域尺度大小的不同，可分為全球?qū)游龀上窈蛥^(qū)域?qū)游龀上瘛Ｓ捎诿娌ū旧淼拈L波長性質(zhì)，面波層析成像主要用于全球范圍或大尺度區(qū)域的研究；而體波波長較短，具有較高的空間分辨率，因此，體波層析成像可用于局部區(qū)域的研究（Zhao，2001）。目前，應(yīng)用較多的依然是體波（P波）的走時層析成像，其所反演的參數(shù)也由單純的介質(zhì)速度結(jié)構(gòu)發(fā)展到介質(zhì)速度結(jié)構(gòu)和震源位置，或者同時反演介質(zhì)速度結(jié)構(gòu)和界面深度（Aki et al，1976、1977；Pavlis et al，1980；劉福田，1984；Bishop et al，1985；Zhao et al，1992；Zelt et al，1992；華標(biāo)龍等，1995；Rawlinson et al，2001；周龍泉等，2006）。Zhao（2001）、Rawlinson等（2010）分別對地震層析成像技術(shù)進(jìn)行了系統(tǒng)的總結(jié)。

地震走時層析成像是建立在一定的模型參數(shù)化和正演射線追蹤基礎(chǔ)之上的。在對地球介質(zhì)進(jìn)行模型參數(shù)化時，根據(jù)其空間分布特征可分為離散介質(zhì)和連續(xù)介質(zhì)2種類型。對于離散介質(zhì)類型一般將地球介質(zhì)進(jìn)行網(wǎng)格劃分，常見的為大小相等的矩形網(wǎng)格（三維為立方體），在網(wǎng)格內(nèi)部或節(jié)點上定義各自的地震波速度等屬性。網(wǎng)格結(jié)構(gòu)存在2個主要困難：一是難于描述復(fù)雜構(gòu)造的精細(xì)結(jié)構(gòu)；二是為了滿足描述最精細(xì)區(qū)域的要求，網(wǎng)格結(jié)構(gòu)通常需要密度較大的網(wǎng)格節(jié)點，而網(wǎng)格節(jié)點上的正演計算通常至少與網(wǎng)格節(jié)點數(shù)目成正比（Moser，1991），因此，受計算機(jī)內(nèi)存和運算速度的影響很大，對于三維空間尤其如此。而離散化模型適應(yīng)性強(qiáng)，在此基礎(chǔ)上的走時計算和射線追蹤方法健壯性好是其最大的優(yōu)點。

層狀結(jié)構(gòu)是常見的分段連續(xù)介質(zhì)描述方法，其以不同的分層界面來描述物性間斷面。層狀結(jié)構(gòu)描述簡單，模型基礎(chǔ)上的正演射線追蹤計算方便快捷（Zelt et al，1992；Zhang et al，2007、2013）。但層狀結(jié)構(gòu)在描述逆斷層等復(fù)雜三維地質(zhì)體時會變得十分困難，而塊狀結(jié)構(gòu)描述則容易的多。塊狀結(jié)構(gòu)模型中，三維地質(zhì)體被看做不同地質(zhì)塊組成的集合體。Gj?ystdal等（1985）首先提出了三維塊狀建模的概念，地質(zhì)塊采用集合運算的操作來構(gòu)造，顯示很不直觀。Pereyra（1996）發(fā)展了塊狀模型描述的技術(shù)，塊體更加直觀自然，但模型界面采用B樣條曲面描述，這使其受到限制，只能構(gòu)造尖滅和蘑菇云等典型地質(zhì)模型。我們在此基礎(chǔ)上采用三角形面片描述模型界面，在構(gòu)建復(fù)雜模型方面有了很大的提升，理論上可以構(gòu)建任意復(fù)雜構(gòu)造的地質(zhì)體（徐濤等，2004；Xu et al，2006、2008、2010；李飛等，2013）。

傳統(tǒng)的運動學(xué)射線追蹤方法包括試射法（Langan et al，1985；Sun，1993；Sambridge et al，1995；徐濤等，2004）和彎曲法（Julian et al，1977；Thurber et al，1980；Um et al，1987；Xu et al，2006、2010、2014；李飛等，2013；Li et al，2014）。試射法在全局搜索接收器方面效果較好，而彎曲法的計算效率較高。Ceverny（2001）對射線追蹤方法進(jìn)行了很好的總結(jié)。在前期的工作中，基于塊狀建模方法我們已經(jīng)發(fā)展了適用于常速度（Xu et al，2006）、常梯度速度分布（Xu et al，2010）的逐段迭代射線追蹤算法，以及VTI介質(zhì)中的子三角形法射線追蹤算法（Xu et al，2008）。之后，我們將塊狀模型介質(zhì)擴(kuò)展到任意速度分布的非均勻介質(zhì)，發(fā)展了與之相適應(yīng)的逐段迭代射線追蹤算法，并結(jié)合偽彎曲法提出了射線追蹤的3點修正方案（李飛等，2013）。三維模型實驗顯示了該射線追蹤擾動修正方案的適用性，相對于試射法，其追蹤效率顯著提高。本文基于三維塊狀建模以及快速高效的逐段迭代射線追蹤，采用模擬退火非線性反演算法，開展了三維速度模型中的地震波走時反演研究。

1 模型參數(shù)化及正演射線追蹤

1.1 三維非均勻塊狀模型

塊狀結(jié)構(gòu)模型中，地質(zhì)體被看作由大小不等、形狀各異的地質(zhì)塊組成的集合體，每個地質(zhì)塊有各自的地質(zhì)屬性，如地震波速度、密度等物性結(jié)構(gòu)，并與其他的地質(zhì)塊存在相鄰的邊界（徐濤等，2004；Xu et al，2006）。我們對二維區(qū)域剖面采用“域-面—邊-點”的描述方式（徐果明等，2001；Xu et al，2010）。不同的地質(zhì)面元為封閉的，被邊界分隔開，邊界由離散點插值得到的3次樣條給定。三維地質(zhì)體的描述采用“體-塊—面-三角形-點”的層次結(jié)構(gòu)（徐濤等，2004；Xu et al，2006、2008、2010；李飛等，2013）。不同地質(zhì)塊之間的物性分界面由一系列的三角形面片拼接而成。

三維塊狀模型中界面的描述極為關(guān)鍵。與Coons、Bezier、B樣條等曲面相比，三角形面片有諸多優(yōu)勢（李飛等，2013）。在三維建模軟件GOCAD系統(tǒng)中，也是利用三角形來構(gòu)造模型界面的（Mallet，1989）。三角形內(nèi)部各點法向量是確定且唯一的，如果2個三角形面片不在1個平面上，則界面法向量在跨三角形的連接處會出現(xiàn)突變，因此，會造成反射、透射射線在跨連接處也會出現(xiàn)突變，這對于如擾動量修正射線軌跡的彎曲法射線追蹤存在極大的困難。針對該問題，我們對界面上各點的法向量進(jìn)行了光滑處理，重新定義界面上各點的法向量，使得法向量在整個界面內(nèi)連續(xù)變化（徐濤等，2004；Xu et al，2006、2008、2010；李飛等，2013）。

圖1為我們采用塊狀建模并結(jié)合三角形界面描述方式設(shè)計的一個典型三維模型，該模型比類伸展盆地中的塹壘構(gòu)造，該構(gòu)造廣泛存在于中國東部如渤海灣盆地、蘇北盆地以及松遼盆地等。模型由18個地質(zhì)塊構(gòu)成，6676個三角形描述模型界面，三角形則由2700個離散點控制。

圖1 復(fù)雜塹壘構(gòu)造地質(zhì)體

圖2 由矩形網(wǎng)格定義的非均勻速度場

為了描述三維模型中非均勻的速度分布，我們對三維地質(zhì)體每個塊體內(nèi)部速度場單獨定義，分別建立1套獨立的矩形網(wǎng)格節(jié)點，并在節(jié)點上定義速度值（圖2）。如圖2所示，對于某塊體內(nèi)部的任意點P（x，y，z），首先找到該點在矩形網(wǎng)格中的位置，再由周圍8個節(jié)點位置處的速度三線性插值獲得P點的速度v（x，y，z）。

式中，v（i＋l，j＋m，k＋n）分別為8個節(jié)點的速度值；xi，yj，zk分別為矩形速度網(wǎng)格在x，y，z方向的位置坐標(biāo)。

綜上所述，我們在三維塊體內(nèi)部定義一系列立方體速度網(wǎng)格，網(wǎng)格內(nèi)部任意位置處速度由其所在速度網(wǎng)格頂點位置處的速度通過線性插值獲得，通過網(wǎng)格頂點速度值控制整個模型的速度分布。由此，三維塊體內(nèi)部可以定義為常速度、常梯度速度、任意非均勻速度等分布。相鄰塊體之間的分界面由一系列的三角形面片拼接而成，三角形頂點位置可以隨意改變，通過三角形頂點位置的抬升與下降來控制模型界面的起伏。

1.2 三維非均勻塊狀模型的射線追蹤

Um等（1987）提出了偽彎曲法射線追蹤用于解決連續(xù)介質(zhì)中的2點射線追蹤問題，但卻無法適應(yīng)存在強(qiáng)速度間斷面的情況。我們提出了逐段迭代射線追蹤算法，用于追蹤存在強(qiáng)速度間斷面模型，該算法屬于彎曲法的范疇，已經(jīng)在二維、三維常速度分布（Xu et al，2006）和常梯度速度分布（Xu et al，2010）中均取得了很好的射線追蹤效果。對于更普遍的非均勻介質(zhì)模型，我們發(fā)展了相適應(yīng)的逐段迭代射線追蹤算法（李飛等，2013）。下面簡述逐段迭代法修正界面路徑點的基本原理。

圖3為逐段迭代法修正中間點示意圖。如圖3所示，對于任意分布的非均勻介質(zhì)，假設(shè)Pk-1、Pk＋1為跨過界面的射線路徑的起始點和終止點；Pk為落在界面上的路徑點，滿足界面方程

圖3 逐段迭代法修正中間點示意圖

如果Pk－1、Pk＋1與Pk點間的距離很小，則射線軌跡可近似看作2個直線段Pk－1Pk和Pk Pk＋1。固定起始點Pk-1和終止點Pk＋1，則射線軌跡的走時T是中間點Pk坐標(biāo)xk（ξ，η）的函數(shù)

式（3）可近似寫為

其中，v0、v3為起始點Pk－1和終止點Pk＋1處的速度；v1、v2為界面路徑點Pk兩側(cè)的速度。

基于射線走時的費馬原理，假設(shè)修正后的路徑點位置為（ξ＋Δξ，η＋Δη），則有如下的偏導(dǎo)數(shù)方程成立

對式（5）進(jìn)行泰勒展開，并保留一階小量，最終可計算得到中間點的一階修正量（Δξ，Δη），詳細(xì)公式推導(dǎo)過程可參見文獻(xiàn)（李飛等，2013）。

針對存在速度間斷面的非均勻介質(zhì)，我們聯(lián)合偽彎曲法和逐段迭代法，提出了1種擾動修正方案：對于落在界面上，即速度間斷面上的路徑點，進(jìn)行逐段迭代法的擾動修正；而對于在模型連續(xù)介質(zhì)內(nèi)部的路徑點，用偽彎曲法來修正，這樣可適應(yīng)三維復(fù)雜非均勻介質(zhì)且存在速度間斷面的情況（李飛等，2013）。與Zhao等（1992）的射線追蹤過程相比，我們在逐段迭代過程中采用一階顯式修正的方式，而不是迭代的方法修正落在界面上的路徑點，因此，追蹤過程更高效省時。

圖4為我們設(shè)計的一個地層互切模型的射線追蹤結(jié)果，模型在x、y、z軸3個方向的尺度均為5km，定義z軸向上為正。地層互切模型（圖4（a））包含5個地質(zhì)塊、3635個三角形及1681個點。模型的x、y、z軸方向分別以紅、綠、藍(lán)3種顏色直線表示（圖4（a））。由圖4（a）可見，直線連接起始點S、終止點R以及與模型界面的交點P，作為初始射線路徑，圖中顯示了追蹤過程中的部分射線路徑變化及最終的射線追蹤結(jié)果。圖4（b）顯示了三維非均勻速度場在y＝2.5km處的速度切片。

圖4 地層互切模型的射線追蹤結(jié)果

綜上所述，我們采用塊狀建模并結(jié)合三角形界面可以構(gòu)建任意復(fù)雜的三維地質(zhì)體，同時發(fā)展了適應(yīng)于非均勻速度分布且存在速度間斷面的逐段迭代法，并結(jié)合偽彎曲法，提出了復(fù)雜塊狀模型中射線追蹤的3點修正方案，該方案具有很強(qiáng)的適應(yīng)性，且相對于試射法追蹤效率顯著提高，這為地震體波走時反演研究打下了堅實的基礎(chǔ)。

2 模擬退火非線性反演算法

地球物理反演理論自誕生至今，已取得了長足的進(jìn)展，形成了一系列較為成熟的線性、非線性反演算法，如最小二乘（Miller，2006）、共軛梯度法（Fletcher et al，1964）、模擬退火（Bertsimas et al，1993）、遺傳算法等。本文采用模擬退火非線性反演算法，進(jìn)行了三維速度模型的地震波走時反演研究。

模擬退火法（Bertsimas et al，1993）是一種用于求解大規(guī)模優(yōu)化問題的隨機(jī)搜索算法。它以優(yōu)化問題求解過程與物理系統(tǒng)退火過程之間的相似性為基礎(chǔ)；優(yōu)化的目標(biāo)函數(shù)相當(dāng)于金屬的內(nèi)能；優(yōu)化問題的自變量組合相當(dāng)于金屬的內(nèi)能狀態(tài)空間；問題的求解過程就是尋找一個組合狀態(tài)，以使得目標(biāo)函數(shù)值最小。該方法利用Metropolis準(zhǔn)則并適當(dāng)?shù)乜刂茰囟鹊南陆颠^程實現(xiàn)模擬退火，從而達(dá)到在多項式時間內(nèi)求解全局優(yōu)化問題的目標(biāo)。

模擬退火算法來源于固體退火原理，將固體加溫至溫度充分高，再讓其徐徐冷卻，加溫時，固體內(nèi)部粒子隨溫度升高變?yōu)闊o序狀，內(nèi)能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達(dá)到平衡態(tài)，最后在常溫時達(dá)到基態(tài)，內(nèi)能減為最小。用固體退火模擬組合優(yōu)化問題，將內(nèi)能E模擬為目標(biāo)函數(shù)值f，溫度T演化成控制參數(shù)t，即得到解組合優(yōu)化問題的模擬退火算法。退火過程由冷卻進(jìn)度表控制，包括控制參數(shù)的初值T0及其衰減因子ΔT，每個T值時的迭代次數(shù)L和停止條件S。模擬退火算法的基本流程如下：

（1）隨機(jī)產(chǎn)生1個初始解x0，令xbest＝x0，并計算目標(biāo)函數(shù)值E（x0）。

（2）設(shè)置初始溫度T（0）＝T0，迭代次數(shù)i＝1。

（3）Do whileT（i）＞Tmin，其中，Tmin為設(shè)定的最低溫度

1）forj＝1～k；

2）對當(dāng)前最優(yōu)解xbest按照某一鄰域函數(shù)產(chǎn)生1個新的解xnew，并計算新的目標(biāo)函數(shù)值E（xnew），目標(biāo)函數(shù)值的增量ΔE＝E（xnew）－E（xbest）；

3）如果ΔE＜0，則xbest＝xnew；

4）如果ΔE＞0，則計算狀態(tài)轉(zhuǎn)移概率P＝exp（－ΔE／T（i）），隨機(jī)產(chǎn)生1個0～1之間的隨機(jī)數(shù)c＝random［0，1］，如果c＜P，則xbest＝xnew，否則xbest＝xbest；

5）End for。

（4）Ti＋1＝Ti－ΔT，i＝i＋1。

（5）End do。

（6）輸出當(dāng)前最優(yōu)點，計算過程結(jié)束。

狀態(tài)轉(zhuǎn)移概率是指從一個狀態(tài)xold（一個可行解）向另一個狀態(tài)xnew（另一個可行解）的轉(zhuǎn)移概率，可以通俗地理解為接受一個新解為當(dāng)前解的概率，它與當(dāng)前的穩(wěn)定參數(shù)T有關(guān)，并隨溫度下降而減小，一般采用Metropolis準(zhǔn)則

模擬退火過程中的冷卻進(jìn)度表是指從某一個高溫狀態(tài)T0向低溫狀態(tài)冷卻時的降溫管理表，時刻t的溫度用T（t）表示。經(jīng)典模擬退火算法的降溫方式為，而快速模擬退火算法的降溫方式為。2種降溫方式都能使得模擬退火算法收斂于全局最小點。其中，初始溫度T0越大，獲得高質(zhì)量解的幾率越大，但耗費的計算時間也會隨之增大。因此，初溫的確定應(yīng)折中考慮優(yōu)化質(zhì)量和優(yōu)化效率。

可以通過增加某些環(huán)節(jié)以實現(xiàn)對模擬退火算法的改進(jìn)，主要的改進(jìn)方式包括：在算法進(jìn)程中的適當(dāng)時機(jī)，增加升溫和重升溫過程，將溫度適當(dāng)提高，可激活各狀態(tài)的接受概率，進(jìn)而調(diào)整搜索進(jìn)程中的當(dāng)前狀態(tài)，避免算法在局部極小解處停滯不前；為避免搜索過程中由于執(zhí)行概率接受環(huán)節(jié)而遺失當(dāng)前遇到的最優(yōu)解，可通過增加存儲環(huán)節(jié)，將當(dāng)前遇到的最優(yōu)解記憶下來；在退火過程結(jié)束后，以搜索到的最優(yōu)解為初始解，再次執(zhí)行模擬退火過程或者局部性搜索，即增加補(bǔ)充搜索過程；還可以結(jié)合其他搜索機(jī)制的算法，如遺傳算法等。

3 模型實驗

運用模擬退火算法，我們進(jìn)行了層狀介質(zhì)中均勻速度反演的模型實驗，理論模型采用圖5所示的簡單層狀模型。理論模型中第1、2兩層介質(zhì)速度分別為4.8、5.4km／s，2層界面深度分別為－1.35、－2.4km，第2層介質(zhì)底界面為地震波反射界面。建立如圖5所示的單炮記錄觀測系統(tǒng)，接收器均勻排列于地表。初始模型中固定界面深度值為理論值，且保持不變，介質(zhì)速度初始值在真實值附近某一鄰域范圍內(nèi)隨機(jī)產(chǎn)生。

圖5 水平層狀均勻介質(zhì)理論模型及射線追蹤結(jié)果示意圖

模擬退火算法的目標(biāo)函數(shù)設(shè)為

優(yōu)化的目標(biāo)即為尋找1組速度值，以使得目標(biāo)函數(shù)值最小。設(shè)定退火初始溫度為T0＝10000，初始溫度越大，則獲得高質(zhì)量解的概率越大。

降溫方式采用指數(shù)衰減形式

其中，β為給定的溫度衰減常數(shù)，如β＝0.2；k為外層循環(huán)次數(shù)。

降溫過程中，在某一確定的溫度狀態(tài)T（k）下，對當(dāng)前最優(yōu)速度解按照某一鄰域函數(shù)產(chǎn)生1組新的解進(jìn)行，并重新計算目標(biāo)函數(shù)。在迭代進(jìn)程中根據(jù)需要可適當(dāng)縮小鄰域的步長。

實際反演過程中采用的模型速度更新方式為

其中，ξ為［－1，1］之間的隨機(jī)數(shù)；k的定義同式（8）。

在迭代進(jìn)程中，如果搜索到的最優(yōu)值連續(xù)若干步保持不變，則認(rèn)為迭代終止，此時獲得的速度參數(shù)即為通過模擬退火非線性反演算法獲得的最優(yōu)解。采用模擬退火非線性反演算法進(jìn)行的水平層狀均勻模型介質(zhì)速度反演結(jié)果如圖6所示，其中，2層介質(zhì)速度真實值分別為v1＝4.8km／s，v2＝5.4km／s。

圖6 模擬退火速度反演結(jié)果隨迭代次數(shù)的變化

4 結(jié)論

塊狀結(jié)構(gòu)結(jié)合三角形界面可描述復(fù)雜的三維模型，在很大程度上克服了層狀結(jié)構(gòu)描述復(fù)雜模型的困難。正演過程中采用的逐段迭代射線追蹤方法可適應(yīng)三維復(fù)雜非均勻介質(zhì)的塊狀模型。模擬退火非線性反演算法實現(xiàn)了三維速度模型中的地震波走時反演。模型測試結(jié)果表明了反演工作的高效性，下一步將進(jìn)一步考慮三維復(fù)雜模型中速度與界面的聯(lián)合反演。