基于程函方程快速行進(jìn)法的地震走時(shí)層析成像法

李天揚(yáng)，李桐林，孟相禹

吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

0 引言

地震走時(shí)層析成像法是一種利用不同地震臺(tái)站連續(xù)接收到的地震事件，利用地震走時(shí)斷層掃描反演地球內(nèi)部速度分布的方法。由1977年發(fā)展至今已有40多年的歷史，已成為如今對(duì)地下速度結(jié)構(gòu)探測(cè)的重要研究手段，在火山活動(dòng)、板塊運(yùn)動(dòng)學(xué)和地幔流等領(lǐng)域的研究中扮演著重要的角色[1]。

地震層析成像過(guò)程主要包括正演和反演兩個(gè)環(huán)節(jié)。早期的層析成像法正演采用的是傳統(tǒng)的射線追蹤方法，例如試射法、彎曲法，但由于其追蹤過(guò)程中存在多路徑問(wèn)題,致使反演不穩(wěn)定且多解，故針對(duì)該問(wèn)題相繼提出了基于程函方程數(shù)值解的地震波走時(shí)計(jì)算方法，主要包括有限差分法[2]、快速推進(jìn)法[3--4]和快速掃描法[5--7]。近幾十年來(lái)，最短路徑法[8--10]等方法發(fā)展起來(lái)。反演環(huán)節(jié)早期的地震層析成像法中應(yīng)用的阻尼最小二乘法、奇異值分解法，隨后發(fā)展的模擬退火法[11]、共軛梯度法、迭代重建法及子空間反演法等諸多方法都可以針對(duì)不同模型應(yīng)用于反演環(huán)節(jié)中。

地震層析成像中不同的波也具有不同的用處，其中最常見(jiàn)的有初至波、反射波。Rawlinson et al.[12]將FMM正演算法與Subspace反演算法相結(jié)合利用遠(yuǎn)震初至波資料研究了澳大利亞Murray盆地的上地幔速度結(jié)構(gòu)；2015年，付翠[13]正演利用旅行時(shí)線性插值射線追蹤法以及反演的最小二乘因式分解法將初至波與反射波聯(lián)合反演應(yīng)用到井間地震中；2017年,俞岱等[14]利用旅行時(shí)線性插值法以及初至波旅行時(shí)反演實(shí)現(xiàn)了初至波層析成像并行運(yùn)算；查樹(shù)貴[15]通過(guò)對(duì)反射波走時(shí)信息加約束條件提高成像的分辨率；Wang et al.[16]利用反射波振幅信息對(duì)模型速度的空間分布進(jìn)行反演計(jì)算；張?jiān)频萚17]在2019年基于多步快速行進(jìn)法和子空間反演法實(shí)現(xiàn)了對(duì)二維起伏層狀介質(zhì)中的P波和反射波的反演。國(guó)內(nèi)目前對(duì)于FMM算法以及利用其進(jìn)行射線追蹤不多見(jiàn)，孫章慶等[18--20]利用快速推進(jìn)法迎風(fēng)雙線性插值法實(shí)現(xiàn)了三維地震波走時(shí)計(jì)算以及在復(fù)雜山地下和復(fù)雜海底中地震波的射線追蹤和走時(shí)計(jì)算；王晗等[21]反演通過(guò)雙曲Radon變換實(shí)現(xiàn)了對(duì)地震數(shù)據(jù)的重建；肖漢等[22]利用快速匹配法計(jì)算得到了VTI介質(zhì)的走時(shí)。筆者在前人的工作基礎(chǔ)上，正演采用基于程函方程快速行進(jìn)法(FMM)進(jìn)行射線追蹤，實(shí)現(xiàn)了對(duì)入射波以及反射波的射線追蹤，反演采用了子空間反演法，并且設(shè)計(jì)了理論模型進(jìn)行反演，并對(duì)各反演結(jié)果進(jìn)行了分析。

1 正演理論

1.1 快速行進(jìn)法解程函方程

快速行進(jìn)法是射線追蹤方法的一種，主要用來(lái)求解非線性的程函方程來(lái)計(jì)算地震波走時(shí)場(chǎng)以及射線路徑。程函方程是指導(dǎo)波前傳播的物理現(xiàn)象表達(dá)式，其二維直角坐標(biāo)系下的程函方程為：

(1)

式中：t為地震波的旅行時(shí)間；s為地震波在介質(zhì)中傳播的慢度。程函方程的解即為對(duì)t的求解，采用迎風(fēng)差分方法對(duì)其離散并簡(jiǎn)化為：

(2)

(3)

Rawlinson and Sambridge將一、二、三階的算子進(jìn)行了精度上和計(jì)算效率上的比較，得到了一階算子精度達(dá)不到標(biāo)準(zhǔn)，二階、三階的精度差異很小，但是二階算子的計(jì)算效率遠(yuǎn)遠(yuǎn)高于三階算子的結(jié)論，故通常應(yīng)用時(shí)選擇二階快速行進(jìn)法。二階迎風(fēng)差分算子為：

(4)

走時(shí)t可通過(guò)對(duì)式(2)的離散程函方程解得，其流程是選取需要計(jì)算的節(jié)點(diǎn)以及該節(jié)點(diǎn)上下左右4個(gè)方向結(jié)點(diǎn)并向波前拓展的方向建立差分格式，解出該方程的黏性解即為該節(jié)點(diǎn)的走時(shí)。

1.2 快速行進(jìn)法解程函方程求取波前的方法

通過(guò)快速行進(jìn)法得到的程函方程的解即為波前。由于快速行進(jìn)法是一種基于網(wǎng)格的射線追蹤方法，因此要描述其在網(wǎng)格中如何將波前推進(jìn)。

如圖1a所示，網(wǎng)格中所有的點(diǎn)劃分為3個(gè)區(qū)域，分別為上風(fēng)區(qū)域、窄帶及下風(fēng)區(qū)域。其中上風(fēng)區(qū)域中的點(diǎn)是在波前的后面，稱為alive點(diǎn)，其走時(shí)是已知；窄帶中的點(diǎn)是在波前上，稱為close點(diǎn)，走時(shí)雖然已經(jīng)算出，但是還沒(méi)有進(jìn)行更新；下風(fēng)區(qū)域中的點(diǎn)在波前的前面，稱為far點(diǎn)，其走時(shí)還未進(jìn)行計(jì)算。

圖1b所示為窄帶原理在源點(diǎn)進(jìn)行波前演變的步驟。alive點(diǎn)的走時(shí)是已經(jīng)計(jì)算出來(lái)的，使用離散化公式計(jì)算出close點(diǎn)走時(shí)，然后選擇其中值最小的變成alive點(diǎn)，far點(diǎn)沒(méi)有計(jì)算值，窄帶是由計(jì)算出的close點(diǎn)不停地向前演變。在這個(gè)過(guò)程中，新的alive點(diǎn)不斷的被標(biāo)記，并持續(xù)更新鄰近的close點(diǎn)和far點(diǎn)，窄帶的傳播形狀近似為初至波前的形狀，快速行進(jìn)法使窄帶沿著網(wǎng)格傳播，計(jì)算結(jié)束的標(biāo)志是所有的點(diǎn)變成alive點(diǎn)。

1.3 射線求取以及反射走時(shí)雅克比矩陣的計(jì)算

從接收點(diǎn)開(kāi)始逐點(diǎn)尋找走時(shí)場(chǎng)梯度的最大下降方向，該方向點(diǎn)的連線即為接收點(diǎn)到源點(diǎn)之間的射線路徑。首先要找到接收點(diǎn)位置并記下，然后設(shè)置每一步射線增量路徑長(zhǎng)度進(jìn)行梯度的計(jì)算，計(jì)算完每一步之后得到的點(diǎn)記下并繼續(xù)上一步的計(jì)算直到設(shè)置的界面位置(此過(guò)程采用的是已經(jīng)計(jì)算過(guò)的反射波場(chǎng))，記下界面上的點(diǎn)并繼續(xù)以上操作直到源點(diǎn)為止(此過(guò)程采用的是已經(jīng)計(jì)算過(guò)的入射波場(chǎng))，連接此過(guò)程的點(diǎn)即為一條入射波和反射波的射線路徑。

在二維直角坐標(biāo)系中，模型任一點(diǎn)P可表示為(x,z)，根據(jù)射線理論，由第i個(gè)源點(diǎn)經(jīng)反射點(diǎn)到第j個(gè)接收點(diǎn)的走時(shí)ti,j可表示為：

(5)

式中：ds為第i個(gè)震源到第j個(gè)檢波器的射線路徑Rij上的線積分元；v(x,z)為沿著該射線路徑上的速度函數(shù)，當(dāng)模型單元網(wǎng)格化后，沿著射線路徑Rij上的走時(shí)可以寫成求和的形式：

(6)

式中：N為射線穿過(guò)的單元總數(shù)；Rij,k、v(x,z)分別為射線穿過(guò)第k個(gè)單元的長(zhǎng)度和速度分布。在對(duì)速度進(jìn)行反演時(shí)，F(xiàn)rechet偏導(dǎo)數(shù)矩陣包含了走時(shí)關(guān)于速度變化的導(dǎo)數(shù)，其公式為：

(7)

當(dāng)射線穿過(guò)某一個(gè)網(wǎng)格單元時(shí)，走時(shí)對(duì)速度的一階偏導(dǎo)數(shù)可用射線兩端點(diǎn)導(dǎo)數(shù)的平均值代替：

(8)

式中的：

(9)

式中：c=ki,kj，ki,kj是穿過(guò)第k個(gè)單元射線段的兩個(gè)端點(diǎn)；BV(k)即為網(wǎng)格單元速度；fc(x,z)即為所求Frechet導(dǎo)數(shù)矩陣也就是雅克比矩陣。

2 反演理論以及算法流程

2.1 反演理論

子空間算法是Skilling et al.[23]在最大熵圖像重建中提出的，隨后被諸多學(xué)者將其概念推廣和改進(jìn)，并用其討論了地震多參數(shù)的同時(shí)反演問(wèn)題，主要用于對(duì)速度以及界面的反演中。

走時(shí)層析成像方法反演可歸結(jié)為如下目標(biāo)函數(shù)求極小值的問(wèn)題，即：

(10)

子空間反演法是在多個(gè)模型參數(shù)的子空間內(nèi)同時(shí)沿著多個(gè)方向去尋找目標(biāo)函數(shù)S(m)的最小值，算法的核心是將多個(gè)參數(shù)的模型擾動(dòng)量用p維子空間的一維基向量{aj},j=1,2,…，p來(lái)表示，即：

(11)

由公式(9)可推出：

(12)

(13)

(14)

②如果p≧2，則構(gòu)造基向量公式如下：

(15)

③采用奇異值分解法將①、②中構(gòu)建的基向量標(biāo)準(zhǔn)正交化。

相比于傳統(tǒng)的其他梯度類算法，其優(yōu)點(diǎn)是只需要對(duì)一個(gè)維數(shù)等于子空間維數(shù)的矩陣進(jìn)行反演，并正確選擇跨越模型的子空間的基向量，這個(gè)優(yōu)點(diǎn)也使其實(shí)用性高于其他算法。

2.2 算法流程

依據(jù)快速行進(jìn)法射線追蹤理論以及子空間反演法為基礎(chǔ)，在Windows環(huán)境下運(yùn)用Fortran語(yǔ)言編程實(shí)現(xiàn)。實(shí)現(xiàn)了對(duì)二維模型快速行進(jìn)法的射線追蹤以及對(duì)模型的反演。其中射線追蹤過(guò)程為：①給定初始模型并設(shè)置好網(wǎng)格范圍以及源點(diǎn)接收點(diǎn)；②在源點(diǎn)附近進(jìn)行源的細(xì)化計(jì)算，提高計(jì)算準(zhǔn)確度；③設(shè)置網(wǎng)格所有點(diǎn)為far點(diǎn)并開(kāi)始計(jì)算；④計(jì)算到網(wǎng)格底面使所有點(diǎn)計(jì)算成為alive點(diǎn)并備份計(jì)算所得網(wǎng)格時(shí)間場(chǎng)；⑤將計(jì)算結(jié)束的網(wǎng)格區(qū)域再次設(shè)置為far點(diǎn)，從保留了時(shí)間場(chǎng)的網(wǎng)格底面開(kāi)始向接收點(diǎn)進(jìn)行計(jì)算；⑥到接收點(diǎn)計(jì)算結(jié)束后完成對(duì)入射波以及反射波波場(chǎng)的計(jì)算；⑦如1.3所述完成對(duì)射線的追蹤。反演流程如圖2所示。正演反演結(jié)束后結(jié)果數(shù)據(jù)通過(guò)matlab和surfer繪制成圖。

圖2 反演流程圖Fig.2 Inversion flow chart

3 模型試算

采用理論模型進(jìn)行試算，網(wǎng)格剖分情況為26 km×30 km。異常體情況如圖3、圖4所示，低速異常體和高速異常體大小形狀位置均相同，其中低速異常體背景速度為5 km/s，低速模型為速度3 km/s；高速異常體背景速度為5 km/s，高速模型速度為8 km/s。

圖3 低速異常體理論模型Fig.3 Theoretical model of low velocity abnormal body

圖4 高速異常體理論模型Fig.4 Theoretical model of high velocity abnormal body

圖5、6、7為分別為低速異常體和入射波場(chǎng)、反射波場(chǎng)及射線路徑圖，由正演計(jì)算所得，在圖中可以看出對(duì)應(yīng)異常體的位置；圖8、9、10分別為高速異常體的入射波場(chǎng)、反射波場(chǎng)及射線路徑圖，圖中可以看出對(duì)應(yīng)位置異常體的位置。

圖5 低速異常體入射波場(chǎng)示意圖Fig.5 Schematic diagram of incident wave field of low velocity abnormal body

圖6 低速異常體反射波場(chǎng)示意圖Fig.6 Schematic diagram of reflection wave field of low velocity abnormal body

圖7 低速異常體射線路徑示意圖Fig.7 Schematic diagram of low velocity abnormal body ray path

圖8 高速異常體入射波場(chǎng)示意圖Fig.8 Schematic diagram of incident wave field of high velocity abnormal body

圖9 高速異常體反射波場(chǎng)示意圖Fig.9 Schematic diagram of reflection wave field of high velocity abnormal body

圖11和圖12為低速異常體和高速異常體的反演結(jié)果，由反演結(jié)果可以看出無(wú)論是低速或者是高速異常體的空間位置與給定模型基本吻合，但是異常體的形狀和大小有些許差異，這是由于射線的密度和射線角度的覆蓋率所導(dǎo)致。此外異常體的速度與背景速度的差異也會(huì)對(duì)反演結(jié)果產(chǎn)生影響。

圖10 高速異常體射線路徑示意圖Fig.10 Schematic diagram of high velocity abnormal body ray path

圖11 低速異常體反演結(jié)果Fig.11 Inversion results of low velocity abnormal body

圖12 高速異常體反演結(jié)果Fig.12 Inversion results of high velocity abnormal body

4 結(jié)論

(1)快速行進(jìn)法(FMM)不僅能追蹤初至波，而且能追蹤反射波。

(2)通過(guò)對(duì)初始模型以及異常體模型的反射波射線追蹤，反射波能夠反映速度結(jié)構(gòu)的變化。

(3)子空間算法反演得到的反射波成像結(jié)果，能夠準(zhǔn)確的恢復(fù)出異常體的空間位置，但是由于射線密度和射線角度的原因?qū)е庐惓ｓw的形狀與大小出現(xiàn)了些許差異。

(4)本文所設(shè)計(jì)固定界面進(jìn)行射線追蹤以及異常體反演，在實(shí)際地質(zhì)條件下界面往往更為復(fù)雜，后續(xù)還需對(duì)復(fù)雜界面下的速度進(jìn)行射線追蹤和反演以及對(duì)界面進(jìn)行反演。