李飛 張雪梅 陳宏峰 何少林

中國(guó)地震臺(tái)網(wǎng)中心，北京市西城區(qū)三里河南橫街5號(hào) 100045

0 引言

地震層析成像從醫(yī)學(xué)CT引入地球物理學(xué)領(lǐng)域用于地殼上地幔速度結(jié)構(gòu)研究（Aki et al，1976）以來(lái)，已發(fā)展成為研究地球內(nèi)部速度結(jié)構(gòu)（白志明等，2016）最重要的方法之一。利用地震層析成像得到地球內(nèi)部高分辨率的速度結(jié)構(gòu)模型，為探索全球構(gòu)造進(jìn)程，巖石圈的演化和構(gòu)造運(yùn)動(dòng)，全球板塊運(yùn)動(dòng)的規(guī)律和動(dòng)力機(jī)制，以及地震、火山活動(dòng)發(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ū)游龀上?。由于面波本身的長(zhǎng)波長(zhǎng)性質(zhì)，面波層析成像主要用于全球范圍或大尺度區(qū)域的研究；而體波波長(zhǎng)較短，具有較高的空間分辨率，因此，體波層析成像可用于局部區(qū)域的研究（Zhao，2001）。目前，應(yīng)用較多的依然是體波（P波）的走時(shí)層析成像，其所反演的參數(shù)也由單純的介質(zhì)速度結(jié)構(gòu)發(fā)展到介質(zhì)速度結(jié)構(gòu)和震源位置，或者同時(shí)反演介質(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）分別對(duì)地震層析成像技術(shù)進(jìn)行了系統(tǒng)的總結(jié)。

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

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

傳統(tǒng)的運(yùn)動(dò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）。試射法在全局搜索接收器方面效果較好，而彎曲法的計(jì)算效率較高。Ceverny（2001）對(duì)射線追蹤方法進(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點(diǎn)修正方案（李飛等，2013）。三維模型實(shí)驗(yàn)顯示了該射線追蹤擾動(dòng)修正方案的適用性，相對(duì)于試射法，其追蹤效率顯著提高。本文基于三維塊狀建模以及快速高效的逐段迭代射線追蹤，采用模擬退火非線性反演算法，開(kāi)展了三維速度模型中的地震波走時(shí)反演研究。

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

1.1 三維非均勻塊狀模型

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

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

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

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

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

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

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

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

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

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

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

圖3 逐段迭代法修正中間點(diǎn)示意圖

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

式（3）可近似寫(xiě)為

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

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

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

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

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

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

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

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

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

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

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

（1）隨機(jī)產(chǎn)生1個(gè)初始解x0，令xbest＝x0，并計(jì)算目標(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）對(duì)當(dāng)前最優(yōu)解xbest按照某一鄰域函數(shù)產(chǎn)生1個(gè)新的解xnew，并計(jì)算新的目標(biāo)函數(shù)值E（xnew），目標(biāo)函數(shù)值的增量ΔE＝E（xnew）－E（xbest）；

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

4）如果ΔE＞0，則計(jì)算狀態(tài)轉(zhuǎn)移概率P＝exp（－ΔE／T（i）），隨機(jī)產(chǎn)生1個(gè)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)點(diǎn)，計(jì)算過(guò)程結(jié)束。

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

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

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

3 模型實(shí)驗(yàn)

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

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

實(shí)際反演過(guò)程中采用的模型速度更新方式為

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

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

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

4 結(jié)論

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