秦 越 ，劉 斌,2

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

1 引 言

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

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

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

2 本文射線追蹤方法

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

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

假設(shè)VP，VS分別對應(yīng)地層的縱波速度和橫波速度，介質(zhì)模型如圖1(a)所示。此時可用一組光滑的函數(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是界面的斜率。

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

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

(4)

對于傾斜層狀均勻介質(zhì)模型，射線路徑與界面的交點的縱坐標z11，z21，z22，z31與橫坐標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ù)T2：

(6)

對于更一般的層狀均勻介質(zhì)模型，射線路徑與界面交點的縱坐標Z11，Z21，Z22，Z31與橫坐標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ì)模型射線追蹤的目標函數(shù)T3：

(8)

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

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

3 數(shù)值例子

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

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

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

4.2 算法實現(xiàn)

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

(9)

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

5 結(jié) 論

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

附錄A

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

(A1)

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

(A2)

求解上述線性方程，可得到使得T1取最小值時的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;

%對炮點循環(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