秦 越 ，劉 斌,2

(1.廣州海洋地質(zhì)調(diào)查局 自然資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 511458；2.南方海洋科學(xué)與工程實(shí)驗(yàn)室，廣東 廣州 511458)

1 引 言

海底地震儀(Ocean Bottom Seismometer, OBS)起初主要應(yīng)用于研究地球深部的構(gòu)造以及地球動(dòng)力學(xué)[1-4]。隨著采集技術(shù)的提高，OBS也越來越廣泛地應(yīng)用于油氣以及水合物的調(diào)查研究中[5-7]。OBS原始數(shù)據(jù)是通過在海底布設(shè)等間距或不等間距的海底地震儀，在海面激發(fā)氣槍震源或接收天然震源而獲得的。通過旅行時(shí)反演來獲得地下介質(zhì)的速度結(jié)構(gòu)是利用OBS數(shù)據(jù)的主要方式[8-12]。旅行時(shí)反演的第一步是拾取OBS數(shù)據(jù)上的初至旅行時(shí)或者反射波旅行時(shí)，并指定反射波對(duì)應(yīng)的反射界面，這一步也稱作震相識(shí)別[13,14]。在震相識(shí)別過程中，常常需要通過射線正演計(jì)算來輔助震相的確定。此外，由于OBS不僅記錄了縱波，還記錄了轉(zhuǎn)換橫波，波場(chǎng)比較復(fù)雜。為分析OBS數(shù)據(jù)，常常需要通過射線追蹤來確定波場(chǎng)的射線路徑。具體地，OBS的射線追蹤就是給定速度模型，計(jì)算指定炮檢點(diǎn)，指定界面和指定波型的射線路徑。

在給定的速度模型上計(jì)算指定的炮檢點(diǎn)所對(duì)應(yīng)的射線路徑屬于微分方程的兩點(diǎn)邊值問題，在數(shù)值計(jì)算領(lǐng)域是一個(gè)比較困難的問題。由于射線追蹤在地震旅行時(shí)反演以及基于射線的偏移方法中處于核心的地位，人們發(fā)展了非常多的方法，包括打靶法[15]，彎曲擾動(dòng)法[16-19]以及二步法[20]等。這些算法能夠適用于均勻和非均勻的速度模型，廣泛應(yīng)用于初至射線以及反射射線的計(jì)算[21]?；谶@些成熟的算法，已經(jīng)開發(fā)出了較多的開源軟件，比如rayinvr[22]，但這些算法的實(shí)現(xiàn)和應(yīng)用均比較復(fù)雜。另外，這些算法不能很好地支持任意界面以及任意模式轉(zhuǎn)換下的橫波射線追蹤。

針對(duì)層狀均勻介質(zhì)模型下OBS射線追蹤的問題，本文實(shí)現(xiàn)了一種簡(jiǎn)單快速，并且能夠適用于任意的縱波和橫波的算法。該算法基于費(fèi)馬原理和最優(yōu)化原理，首先建立旅行時(shí)的目標(biāo)函數(shù)，然后通過最小化目標(biāo)函數(shù)來獲得給定炮檢點(diǎn)的射線路徑。

2 本文射線追蹤方法

2.1 目標(biāo)函數(shù)的建立

對(duì)于給定的炮點(diǎn)和檢波點(diǎn)，以所有可能的射線路徑作為自變量，建立旅行時(shí)的目標(biāo)函數(shù)。對(duì)于一般的速度模型，射線路徑需要用一個(gè)空間的函數(shù)來表示，所建立的目標(biāo)函數(shù)實(shí)際是關(guān)于射線路徑函數(shù)的泛函。對(duì)于層狀均勻介質(zhì)模型，射線在一個(gè)地層內(nèi)是直線，所以射線路徑可以用射線與界面的交點(diǎn)的坐標(biāo)來表示，這樣，目標(biāo)函數(shù)就是關(guān)于一系列坐標(biāo)的多元函數(shù)。

假設(shè)VP，VS分別對(duì)應(yīng)地層的縱波速度和橫波速度，介質(zhì)模型如圖1(a)所示。此時(shí)可用一組光滑的函數(shù)來描述地層的界面：

zi=fi(x)

(1)

由于水平、傾斜層狀均勻介質(zhì)是彎曲層狀介質(zhì)的特例，因此式(1)中的函數(shù)fi(x)用一組特定的深度值描述，便可得到水平層狀均勻介質(zhì)模型(圖1b)：

zi=li(x)=zi0

(2)

其中，zi0為模型最左邊的深度。

相應(yīng)地，在式(1)中，函數(shù)fi(x)用界面最左邊的深度值以及斜率來描述，可以獲得傾斜層狀介質(zhì)模型(圖1c)：

zi=li(x)=zi0+kix

(3)

其中，ki是界面的斜率。

下面以來自第三個(gè)反射界面的反射波為例(圖1)，推導(dǎo)上述層狀均勻介質(zhì)模型的射線路徑所對(duì)應(yīng)的旅行時(shí)目標(biāo)函數(shù)。對(duì)于其他情形的目標(biāo)函數(shù)可采用相同的方法推導(dǎo)。已知OBS的位置(xobs,zobs)和炮點(diǎn)的位置(xs,zs)，假設(shè)反射波的射線與海底的交點(diǎn)為(x31,z31)，與第二個(gè)反射界面的交點(diǎn)為(x21,z21)以及(x22,z22)，與第一個(gè)反射界面的坐標(biāo)為(x11,z11)?；谶@些坐標(biāo)，可以建立反射波旅行時(shí)的目標(biāo)函數(shù)。

對(duì)于簡(jiǎn)單的水平層狀均勻介質(zhì)模型，其對(duì)應(yīng)的旅行時(shí)目標(biāo)函數(shù)為：

(4)

對(duì)于傾斜層狀均勻介質(zhì)模型，射線路徑與界面的交點(diǎn)的縱坐標(biāo)z11，z21，z22，z31與橫坐標(biāo)x11，x21，x22，x31有關(guān)，其關(guān)系用如下公式表示：

z11=l1(x11)

(5-1)

z21=l2(x21)

(5-2)

z22=l2(x22)

(5-3)

z31=l3(x31)

(5-4)

用上述公式替換公式(1)中的z11，z21，z22，z31,則得到傾斜層狀均勻介質(zhì)的旅行時(shí)目標(biāo)函數(shù)T2：

(6)

對(duì)于更一般的層狀均勻介質(zhì)模型，射線路徑與界面交點(diǎn)的縱坐標(biāo)Z11，Z21，Z22，Z31與橫坐標(biāo)x11，x21，x22，x31之間的關(guān)系為：

z11=f1(x11)

(7-1)

z21=f2(x21)

(7-2)

z22=f2(x22)

(7-3)

z31=f3(x31)

(7-4)

把上述公式帶入公式(1)，得到一般的層狀均勻介質(zhì)模型射線追蹤的目標(biāo)函數(shù)T3：

(8)

2.2 最優(yōu)化求解射線路徑

根據(jù)費(fèi)馬原理，在所有可能的路徑中，使得旅行時(shí)最小的路徑就是真實(shí)的射線路徑。為獲得真實(shí)的路徑，需要求解目標(biāo)函數(shù)即式(4)，式(6)，式(8)的最小值。這樣射線追蹤的問題就轉(zhuǎn)化為數(shù)值最優(yōu)化的問題。對(duì)于水平層狀介質(zhì)的情形，目標(biāo)函數(shù)T1的最小值問題，可通過令函數(shù)的梯度等于零來求解。具體推導(dǎo)詳見附錄A部分。而對(duì)于一般的層狀介質(zhì)情形，也即T3， 在大部分情況下，f1(x)，f2(x)，f3(x)沒有解析表達(dá)式或者梯度的表達(dá)式過于復(fù)雜，此時(shí)目標(biāo)函數(shù)的最小值問題不能通過令梯度等于零的方式來求解，只能通過最優(yōu)化方法來求解。最優(yōu)化是一個(gè)重要的數(shù)學(xué)分支，它所研究的問題是討論在眾多的方案中什么樣的方案最優(yōu)以及怎樣找出最優(yōu)方案[23]。最優(yōu)化方法在交通運(yùn)輸、管理、電力、航天、通信等廣大工程學(xué)科中應(yīng)用非常廣泛。至今最優(yōu)化理論已發(fā)展出線性規(guī)劃、非線性規(guī)劃、幾何規(guī)劃、動(dòng)態(tài)規(guī)劃等多種方法。近年來隨著計(jì)算機(jī)技術(shù)的發(fā)展，遺傳算法[18,24]、粒子群優(yōu)化算法[25,26]、差分進(jìn)化算法[27]等進(jìn)化算法在最優(yōu)化領(lǐng)域也得到了廣泛的應(yīng)用。最近也涌現(xiàn)出了基于深度強(qiáng)化學(xué)習(xí)的最優(yōu)化新方法[28]。由于本文求解射線路徑最優(yōu)化問題的參數(shù)較少，因此采用基于梯度最優(yōu)化方法來求解目標(biāo)函數(shù)的極小值，其中梯度是關(guān)鍵。附錄B給出了水平層狀介質(zhì)情形下的Matlab代碼。

3 數(shù)值例子

為檢驗(yàn)上述算法的有效性，建立水平層狀，傾斜層狀以及一般層狀三個(gè)層狀均勻介質(zhì)模型，并用波動(dòng)方程計(jì)算相應(yīng)的地震記錄(圖2～圖7)。其中，三個(gè)模型的大小一樣，橫向大小為5 000 m，縱向大小為2 000 m(圖2、圖4和圖6)。以第三個(gè)反射界面上發(fā)射波的射線路徑為例，對(duì)每一個(gè)模型采用本文提出的算法進(jìn)行OBS射線追蹤。在顯示射線路徑時(shí)，每隔500 m的炮點(diǎn)顯示一條射線。計(jì)算相應(yīng)的地震記錄，通過在地震記錄上疊合顯示射線追蹤計(jì)算得到的旅行時(shí)來檢驗(yàn)算法的有效性(圖3，圖5和圖7)。

圖2 水平層狀均勻速度模型Fig.2 Horizontally layered homogeneous velocity model

圖3 水平層狀均勻速度模型的射線路徑及旅行時(shí)(黑色)在地震數(shù)據(jù)上的疊合顯示 Fig.3 Ray-paths for the horizontally layered homogeneous and the travel time stacked on the seismic data

圖4 傾斜層狀均勻速度模型Fig.4 Homogeneous velocity model with tilted layers

圖5 傾斜層狀均勻速度模型的射線路徑及旅行時(shí)(黑色)在地震數(shù)據(jù)上的疊合顯示Fig.5 Ray-paths for the tilted layered homogeneous and the travel time stacked on the seismic data

圖6 一般層狀均勻速度模型Fig.6 Homogeneous velocity model with curved layers

圖7 一般層狀均勻速度模型的射線路徑及旅行時(shí)(黑色)在地震數(shù)據(jù)上的疊合顯示 Fig.7 Ray-paths for the layered homogeneous with curved layers and the travel time stacked on the seismic data

4 討 論

4.1 適用性

射線追蹤是地震學(xué)領(lǐng)域中非常重要的一個(gè)方面，在射線類偏移和反演算法中起著非常重要的作用，發(fā)展了很多成熟的算法。本文針對(duì)層狀均勻模型下OBS射線追蹤這一問題，實(shí)現(xiàn)了一種簡(jiǎn)單快速、適用于任意的縱波和橫波的算法。當(dāng)射線追蹤用于輔助OBS震相識(shí)別或OBS數(shù)據(jù)波場(chǎng)分析時(shí)，所使用的速度模型一般是層狀均勻的速度模型，這也是本文算法的基礎(chǔ)。該算法基于費(fèi)馬原理，通過求解旅行時(shí)函數(shù)的最小值問題來獲得指定炮點(diǎn)和檢波點(diǎn)，指定界面以及指定波型(縱波或橫波)對(duì)應(yīng)的射線路徑。它能夠穩(wěn)定、快速地計(jì)算出射線路徑。數(shù)值例子表明了該算法的有效性。該算法不僅適用于縱波射線追蹤，還適用于橫波。在計(jì)算轉(zhuǎn)換波的射線路徑時(shí)，只需要在旅行時(shí)公式中把相應(yīng)的速度替換成橫波速度即可。該算法也適用于其他觀測(cè)系統(tǒng)，比如海面拖纜地震。此外，該算法的應(yīng)用并不局限于震相識(shí)別和波場(chǎng)分析，還可以用于特定界面射線的正演研究，一個(gè)典型的應(yīng)用就是計(jì)算與莫霍面相關(guān)的射線。

4.2 算法實(shí)現(xiàn)

通過2.1節(jié)方法部分的公式推導(dǎo)可知，對(duì)于某一類特定的射線，首先需要定義一個(gè)目標(biāo)函數(shù)，然后通過極小化該目標(biāo)函數(shù)來獲得射線路徑。特別地，對(duì)于最簡(jiǎn)單的水平層狀均勻介質(zhì)模型，由于其梯度較為簡(jiǎn)單，可通過令梯度等于零得到一個(gè)線性方程組，通過求解該方程來獲得射線與界面的交點(diǎn)。而對(duì)于其他的情形，則需要通過非線性最優(yōu)化方法來實(shí)現(xiàn)目標(biāo)函數(shù)的極小化。其中，基于梯度下降的優(yōu)化方法是最常用的。對(duì)于該類優(yōu)化方法，梯度計(jì)算是最為關(guān)鍵的部分。當(dāng)射線路徑涉及的地層數(shù)目較少時(shí)，計(jì)算梯度的公式比較簡(jiǎn)單。當(dāng)?shù)貙虞^多，或者射線與界面的交點(diǎn)比較多時(shí)，梯度計(jì)算的公式非常繁瑣。為避免梯度公式推導(dǎo)的繁瑣，可以采用數(shù)值方法來計(jì)算梯度，然后用非線性最優(yōu)化方法來求解極小值。即使對(duì)于水平層狀均勻介質(zhì)模型，也采用這種方式來處理。用這種方式處理時(shí)，對(duì)于不同的射線類型，只需要定義相應(yīng)的目標(biāo)函數(shù)即可，而不再需要去推導(dǎo)其梯度。數(shù)值方式計(jì)算梯度的公式如下：

(9)

其中，Δx是一個(gè)比較小的值?？紤]到T(x11,x21,x22,x31)是簡(jiǎn)單的函數(shù)，采用數(shù)值方法計(jì)算梯度的效率非常高。此外，在最優(yōu)化階段，也可采用共軛梯度、擬牛頓等方法來提高收斂的速度。

5 結(jié) 論

射線追蹤廣泛用于輔助OBS數(shù)據(jù)的震相識(shí)別以及波場(chǎng)分析。基于層狀均勻介質(zhì)模型，本文實(shí)現(xiàn)了一種快速簡(jiǎn)單，并且能夠適用于任意縱波和橫波的算法。該算法基于費(fèi)馬原理和最優(yōu)化方法。對(duì)于指定的炮點(diǎn)和檢波點(diǎn)，指定的界面和指定的波型，首先定義一個(gè)目標(biāo)函數(shù)，然后通過求取該目標(biāo)函數(shù)的最小值來獲得射線路徑。與其他射線追蹤實(shí)現(xiàn)方法相比，該方法具有兩個(gè)優(yōu)勢(shì)：①能夠很好地實(shí)現(xiàn)任意縱波和橫波的射線路徑追蹤；②實(shí)現(xiàn)簡(jiǎn)單，效率高。不同復(fù)雜程度均勻介質(zhì)模型的數(shù)值例子表明了該方法的有效性。

附錄A

對(duì)于水平層均勻介質(zhì)模型，為獲得真實(shí)的射線路徑，求解公式(1)的最小值問題?？紤]到對(duì)于正數(shù)，根號(hào)函數(shù)與根號(hào)的平方單調(diào)性是一致的。公式(1)的最小值問題可以轉(zhuǎn)化為如下公式的最小值問題：

(A1)

令T1的梯度等于0可獲得T1的極小值，得到線性方程組：

(A2)

求解上述線性方程，可得到使得T1取最小值時(shí)的x11,x21,x22,x31，也即得到射線路徑。

附錄 B Matlab 代碼

代碼1：T12332_flat.m

function [f,g] = T12332_flat(X,xs,zs,xobs,zobs,z1,z2,z3,v1,v2,v3,v4,v5)

% the reflection 12332

% OBS position, xobs,zobs,

% source position, xs,zs,

% v1: the velocity for the first segment

% v2: the velocity for the second segment

% ....

x11=X(1);

x21=X(2);

x31=X(3);

x22=X(4);

seg1=sqrt((xs-x11)^2+(zs-z1)^2);

seg2=sqrt((x11-x21)^2+(z1-z2)^2);

seg3=sqrt((x21-x31)^2+(z2-z3)^2);

seg4=sqrt((x31-x22)^2+(z2-z3)^2);

seg5=sqrt((x22-xobs)^2+(z2-zobs)^2);

f=seg1/v1+seg2/v2+seg3/v3+seg4/v4+seg5/v5;

g1=(x11-xs)/seg1/v1+(x11-x21)/seg2/v2;

g2=(x21-x11)/seg2/v2+(x21-x31)/seg3/v3;

g3=(x31-x21)/seg3/v3+(x31-x22)/seg4/v4;

g4=(x22-x31)/seg4/v4+(x22-xobs)/seg5/v5;

g=[g1,g2,g3,g4];

End

代碼2：example_code.m

%#########################################################################

%水平層狀介質(zhì)最優(yōu)化射線追蹤代碼

zs=0;

xobs=2 500; zobs=500;

z1=500;

z2=1 000;

z3=1 500;

v1=1 500;

v2=2 000;

v3=2 500;

v4=2 500;

v5=2 000;

%對(duì)炮點(diǎn)循環(huán)

for xs=500:25:4 500

fh=@(X)T12332_flat(X,xs,zs,xobs,zobs,z1,z2,z3,v1,v2,v3,v4,v5);

X0=[1 025,1 500,1 725,2 000];

XF = fminsearch(fh,X0);

x11=XF(1);

x21=XF(2);

x31=XF(3);

x22=XF(4);

if mod(xs,500)==0

line_X=[xs,x11,x21,x31,x22,xobs];

line_Z=[zs,z1,z2,z3,z2,zobs];

line(line_X,line_Z,'Color',[.8.8.8]);

end

flipy;

end