楊尚倍 白超英 ZHOU Bing

(①中國地質(zhì)調(diào)查局西安地質(zhì)調(diào)查中心，陜西西安 710054；②長安大學(xué)地質(zhì)工程與測繪學(xué)院地球物理系，陜西西安 710054；③哈利法科學(xué)與技術(shù)大學(xué)地球科學(xué)系，阿布扎比 2533)

0 引言

地震波場數(shù)值模擬不僅是了解地震波傳播規(guī)律和解釋地震數(shù)據(jù)的重要工具，也是地震成像的基礎(chǔ)。隨著計算機(jī)技術(shù)快速發(fā)展，各種地震波場數(shù)值模擬算法日趨完善，如有限差分法[1]、有限元法[2]、邊界元法[3]、偽譜法[4-5]、譜元法[6]、有限體積法[7]等。盡管3D波場數(shù)值解相對于2D更符合實際觀測記錄，然而巨大的內(nèi)存需求和運(yùn)算量使得實際3D大尺度波場數(shù)值模擬及波形反演面臨計算資源的巨大挑戰(zhàn)。因此，目前更多的實際應(yīng)用則在2D或小尺度的3D空間開展，如逆時偏移成像[8]和全波形反演[9]。

雖然2D數(shù)值模擬相較于3D高效，但引入了線源假設(shè)代替實際的點源激發(fā)。這種線源假設(shè)意味著震源激發(fā)大小沿走向不變，模擬的波場也是2D，其結(jié)果很難與實際觀測的3D地震記錄相匹配[10-11]。當(dāng)然可采用3D到2D地震數(shù)據(jù)轉(zhuǎn)換方法，將實際的點源數(shù)據(jù)(3D)轉(zhuǎn)換成線源數(shù)據(jù)(2D)，但目前這些轉(zhuǎn)換方法只適用于簡單的地質(zhì)模型，例如均勻介質(zhì)和水平層狀介質(zhì)，并且需要滿足遠(yuǎn)場近似和無波場重疊等苛刻條件[12]，不利于實際復(fù)雜介質(zhì)的成像。

鑒于此，為了克服2D波場模擬中線源假設(shè)帶來的失真和3D波場模擬中效率低下的不足，人們提出了2.5D波場模擬方法[13]。2.5D地震波場數(shù)值模擬方法假設(shè)地質(zhì)模型沿走向方向模型屬性不變，因此可以沿走向方向做傅里葉變換將3D問題簡化為重復(fù)的2D問題，從而大大提高了計算效率，且與3D解析解高度吻合[14]。如：Song等[13]在頻率域通過有限差分法重構(gòu)了3D地震波場；Zhou等[15]給出了精確的2.5D吸收邊界條件，并采用有限元法得到了頻率域2.5D波場數(shù)值解；Maokun等[16]則利用近似最佳正交有限差分算法進(jìn)行了2.5D彈性波數(shù)值模擬；Takenaka等[17]給出了時間域平面波入射下2.5D彈性波動力學(xué)方程的求解過程；Novais等[18]改進(jìn)了時間域2.5D地震波場有限差分?jǐn)?shù)值模擬方法的波數(shù)采樣方法。Xiong等[19]首先將空間域彈性波方程轉(zhuǎn)換為波數(shù)域彈性波方程，然后采用交錯網(wǎng)格有限差分求解波數(shù)域方程。

然而，上述時間域或頻率域2.5D地震波場數(shù)值模擬方法只針對單一的彈性、黏彈性或各向異性介質(zhì)，并沒有一種較為普適的算法，即可同時滿足不同介質(zhì)屬性(彈性、黏彈性、各向異性)和不同邊界條件(聲波自由地表，固體自由地表和固—液邊界)下的2.5D波場數(shù)值模擬算法。本文提出一種適合復(fù)雜多種介質(zhì)混合模型下時間域廣義2.5D一階波動方程及數(shù)值解法，適用于混合(聲波、彈性各向同性、彈性各向異性)介質(zhì)和各種邊界條件(聲波自由地表、固體自由地表和固—液邊界)的地震波場數(shù)值模擬。

另一方面，為了適應(yīng)復(fù)雜地表起伏模型中的數(shù)值模擬，本文采用曲網(wǎng)格進(jìn)行模型剖分，并采用低頻散、低耗散的同位網(wǎng)格MacCormack有限差分法格式[20]離散時間域廣義2.5D一階波動方程，從而得到高精度的數(shù)值解。最后，應(yīng)用數(shù)值模擬算例驗證了本文2.5D數(shù)值模擬方法的正確性和有效性。

1 時間域2.5D廣義一階速度—應(yīng)力 波動方程

首先建立廣義2.5D一階速度—應(yīng)力波動方程，并分析該方程的普適性問題。2.5D地震波場數(shù)值模擬是3D地震波場數(shù)值模擬的一種特殊情況，其假設(shè)地質(zhì)模型沿走向(y方向)具有不變的物理性質(zhì)，則2.5D地震波動方程可以通過對3D地震波動方程沿y方向做傅里葉變換得到[21]。

1.1 彈性各向異性介質(zhì)波動方程

三維彈性各向異性介質(zhì)中，地震波的控制方程[22]為

(1)

σ=cu

(2)

(3)

和

(4)

(5)

式中上標(biāo)“(S)”表示彈性各向同性或彈性各向異性的固體介質(zhì)。式(3)和式(4)可簡寫為矩陣形式

(6)

上式為任意彈性各向異性介質(zhì)中的2.5D地震波控制方程，系數(shù)矩陣A(S)、B(S)、C(S)的具體表達(dá)式見附錄A。

1.2 聲波介質(zhì)波動方程

在聲波介質(zhì)中，應(yīng)力滿足切向應(yīng)力為零和法向應(yīng)力的軸向分量等于液體壓力

σij=-δijP

(7)

P可表示為

P=-Kuj,j+s(t)δ(x-xs)

(8)

式中：K為體積模量；s(t)為xs處的震源函數(shù)。將式(7)代入式(1)，可得

(9)

對式(8)求時間導(dǎo)數(shù)，可得

(10)

根據(jù)式(9)，式(10)可以表示為

(11)

對式(9)和式(11)做關(guān)于y的傅里葉變換，可得

(12)

(13)

定義

(14)

可得式(12)和式(13)的矩陣形式

(15)

式中：上標(biāo)“(W)”表示聲波介質(zhì)；系數(shù)矩陣A(W)、B(W)、C(W)的具體表達(dá)式見附錄B。上式為聲波介質(zhì)中的2.5D地震波控制方程。

1.3 聲波自由地表

聲波自由地表邊界條件需要滿足地表處液體壓力P=0，代入式(12)和式(13)，可得

(16)

(17)

式中：上標(biāo)“(AW)”表示聲波自由地表；系數(shù)矩陣A(AW)、B(AW)、C(AW)分別為

(18)

(19)

C(AW)=0

(20)

對比式(18)～式(20)與附錄B，可以發(fā)現(xiàn)式(17)是式(15)的特例。

1.4 固體自由地表

固體自由地表邊界條件需要滿足地表處法向應(yīng)力為0，即

σ·n=0

(21)

展開為

(22)

式中n=(n1，n2，n3)表示界面的法向向量。將式(22)寫成矩陣形式

(23)

將上式寫為子矩陣格式

Γ2σ2=-Γ1σ1

(24)

式中

(25)

(26)

定義

(27)

式(23)可改寫為

σ2=Sσ1

(28)

由上式可知，應(yīng)力分量σ2可由σ1和表面法向分量n得到，只有應(yīng)力分量σ1未知。假設(shè)自由地表由z=z0(x)給出，其坡度為tanα=z′0(x)，法向分量為n=(-sinα,0,cosα)，則式(27)變?yōu)?/p>

(29)

根據(jù)式(28)和式(29)可得

(30)

將上式代入式(3)，可得

(31)

由式(4)可得

(32)

定義

(33)

則式(30)和式(31)可簡寫為

(34)

1.5 固—液邊界

在固—液邊界上，例如由z=z0(x)定義的二維海底界面，邊界條件需要滿足法向應(yīng)力和法向速度分量連續(xù)，并且切向應(yīng)力分量為0，即

(35)

(36)

(37)

(38)

將式(38)代入(37)，可得

(39)

(40)

(41)

為了表示簡單，假設(shè)固液邊界上沒有施加震源，將式(4)和式(13)代入(41)，可得

(42)

上式對于任意ky成立，可得

(43)

由式(35)可知

(44)

將式(42)代入式(4)，可得

(45)

將式(39)和式(40)代入式(3)，可得

(46)

式中ρS為固體密度。由式(12)可知

(47)

式中ρW為液體密度。定義

(48)

將式(45)～式(47)寫成矩陣形式

(49)

至此，獲得了2.5D廣義一階速度—應(yīng)力波動方程

(50)

式中的各個參量可根據(jù)不同介質(zhì)或不同邊界條件而有所不同。

2 曲線網(wǎng)格有限差分算法

2.1 曲線坐標(biāo)系與直角坐標(biāo)系的轉(zhuǎn)換關(guān)系

為了使網(wǎng)格線與物理空間的邊界線重合，避免傳統(tǒng)矩形網(wǎng)格由于階梯狀近似產(chǎn)生的虛假散射問題，應(yīng)用曲線網(wǎng)格(貼體網(wǎng)格)對地質(zhì)模型進(jìn)行剖分，然后將實際的非規(guī)則物理空間變換到一個規(guī)則的計算空間[23]。假設(shè)物理空間的坐標(biāo)(x,z)和計算空間的坐標(biāo)(ξ,η)之間映射關(guān)系(圖1)為

圖1 曲線網(wǎng)格和計算網(wǎng)格的變換關(guān)系示意圖

(51)

由鏈?zhǔn)角髮?dǎo)法則可得

(52)

整理可得到坐標(biāo)轉(zhuǎn)換系數(shù)

(53)

式中J=x,ξz,η-x,ηz,ξ。為了避免網(wǎng)格非均勻引發(fā)的虛假震源項，x,ξ、x,η、z,ξ、z,η均采用數(shù)值方法計算，所采用的數(shù)值法與下文求解波動方程空間導(dǎo)數(shù)的數(shù)值算法相同。

2.2 曲線坐標(biāo)系中的廣義一階速度—應(yīng)力波動方程

根據(jù)求導(dǎo)的鏈?zhǔn)椒▌t，廣義一階速度—應(yīng)力彈性波方程的矩陣形式在曲線坐標(biāo)系中的表達(dá)式為

(54)

Hixon等[24]采用Tam等[25]提出的DRP (Dispersion Relation Preserving)方法對同位MacCormack格式的系數(shù)進(jìn)行了優(yōu)化，得到了低頻散、低耗散的DRP/opt MacCormack 格式，本文采用該格式有限差分算子離散方程。

同位網(wǎng)格MacCormack格式的兩個偏心差分格式——向前和向后差分格式分別為

(55)

(56)

差分系數(shù)為

(57)

采用DRP/opt MacCormack有限差分格式，式(54)空間導(dǎo)數(shù)可以離散為

(58)

在時間方向上采用四階Runge-Kutta法[26]，結(jié)合DRP/opt MacCormack有限差分算子，波場可以從時間步長第N步更新到第N+4步[20]，具體過程如下。

(59)

(60)

(61)

(62)

式中：α1=0.5；α2=0.5；α3=1.0；β1=1/6；β2=1/3；β3=1/3；β4=1/6。

3 震源、人工邊界和波數(shù)域采樣

3.1 震源和吸收邊界

本文在網(wǎng)格點上施加集中力源，并使用Ricker子波做為震源

(63)

式中：震源中心頻率fc=30Hz；延遲時間t0=30ms。為了消除人工邊界產(chǎn)生的邊界反射，本文采用新近提出的廣義剛度衰減法(The Generalized Stiffness Reduction Method，GSRM)[27]，吸收效果優(yōu)于最佳匹配層方法，適用于各種介質(zhì)，如聲波介質(zhì)、各向同性和各向異性介質(zhì)。

3.2 波數(shù)域采樣方法

(64)

(65)

式中vij表示施加i方向集中力源時得到的j方向速度分量。

(66)

4 2.5D與3D數(shù)值模擬的精度和效 率分析

表1 介質(zhì)彈性參數(shù)

圖2為三個地表水平均勻介質(zhì)模型的2.5D數(shù)值解的波場快照，從圖中可以清晰地看到波場的P、S波波前。為了驗證本文方法模擬得到的2.5D數(shù)值解的計算精度和正確性，將2.5D數(shù)值解與對應(yīng)的3D解析解及3D數(shù)值解進(jìn)行對比(圖3)，三者匹配良好，驗證了本文算法的計算精度和正確性。

圖2 地表水平均勻介質(zhì)模型2.5D模擬的波場快照(a)聲波介質(zhì)的P分量，0.15s；(b)、(c)分別為各向同性介質(zhì)的vx和vz分量，0.09s；(d)、(e)分別為VTI-1介質(zhì)的中vx和vz分量，0.08s

圖3 地表水平均勻介質(zhì)模型2.5D數(shù)值解與3D數(shù)值解、3D解析解的對比(a)聲波介質(zhì)中P分量；(b)、(c)分別為各向同性介質(zhì)的vx和vz分量；(d)、(e)分別為VTI-1介質(zhì)的vx和vz分量

2.5D與3D數(shù)值模擬計算機(jī)內(nèi)存占用和CPU運(yùn)行時間的對比如表2所示，可以看出，2.5D數(shù)值模擬所需CPU時間和計算機(jī)內(nèi)存都遠(yuǎn)小于3D，特別是計算機(jī)內(nèi)存。所以本文的2.5D波場數(shù)值模擬方法則具有較高的實際應(yīng)用價值，可以直接應(yīng)用于地震數(shù)據(jù)反演方法。

表2 2.5D與3D數(shù)值模擬用時和內(nèi)存占用對比

5 復(fù)雜模型的波場數(shù)值模擬

5.1 起伏地表均勻模型

應(yīng)用起伏地表聲波、彈性各向同性和彈性VTI-1均勻介質(zhì)模型(圖4)驗證本文施加自由地表邊界算法的正確性。網(wǎng)格剖分間距Δx=Δz=2m，震源位于(500m，400m)，檢波器位于(300m，500m)，模型參數(shù)見表1。

圖4 模型網(wǎng)格剖分示意圖紅星和藍(lán)三角分別為震源和檢波器

圖5、圖6、圖7分別為起伏地表聲波、彈性各向同性和彈性VTI-1均勻介質(zhì)模型的2.5D數(shù)值解的vx和vz分量波場快照，從中可以清楚地看到直達(dá)P波及自由地表產(chǎn)生的反射波。圖8、圖9、圖10為三個起伏地表模型2.5D與2D數(shù)值解的波形對比，其中2D數(shù)值解是令式(50)中ky=0后采用本文相同方法獲得。2.5D數(shù)值解和2D數(shù)值解之間不僅存在著明顯的振幅差異，還存在相移。

圖5 起伏地表聲波均勻介質(zhì)模型2.5D數(shù)值解vx(a)和vz(b)分量在0.4s的波場快照

圖6 起伏地表各向同性均勻彈性介質(zhì)模型2.5D數(shù)值解vx(a)和vz(b)分量在0.32s的波場快照

圖7 起伏地表均勻VTI-1介質(zhì)模型2.5D數(shù)值解vx(a)和vz(b)分量在0.32s的波場快照

圖8 起伏地表聲波均勻介質(zhì)模型2.5D與2D數(shù)值解對比(a)vx分量；(b)vz分量；(c)圖a的歸一化顯示；(d)圖b的歸一化顯示。下圖為上圖方框的放大顯示

圖9 起伏地表均勻各向同性彈性介質(zhì)模型2.5D與2D數(shù)值解的對比(a)vx分量；(b)vz分量；(c)圖a的歸一化顯示；(d)圖b的歸一化顯示。下圖為上圖方框的放大顯示

圖10 起伏地表均勻VTI-1介質(zhì)模型2.5D與2D數(shù)值解的對比(a)vx分量；(b)vz分量；(c)圖a的歸一化顯示；(d)圖b的歸一化顯示。下圖為上圖方框的放大顯示

5.2 固—液耦合模型

建立一個正弦起伏固—液耦合界面模型(圖11)驗證本文的2.5D數(shù)值模擬方法適用性。模型尺寸為1000m×1000m，網(wǎng)格間距為2m，第一層為聲波介質(zhì)，第二次為各向同性介質(zhì)(參數(shù)見表1)。震源位于(500m，300m)，5個檢波器等間距布置于(300m，0m)與(300m，400m)之間。圖11為2.5D數(shù)值解的0.28s時刻波場照，從圖中可以觀察到直達(dá)P波經(jīng)過固—液界面產(chǎn)生了清晰的透射P波和S波；圖12是2.5D與2D數(shù)值解的波形對比，可以看出2.5D與2D數(shù)值解之間同樣存在著明顯的振幅差異和相移現(xiàn)象。

圖11 固—液耦合模型2.5D數(shù)值解vx(a)和vz(b)分量在0.28s的波場快照圖中起伏曲線為正弦起伏界面

圖12 固—液耦合模型2.5D與2D數(shù)值解波形對比(a)vx分量；(b)vz分量

5.3 三層混合模型

為了檢驗本文2.5D波場數(shù)值模擬方法對于復(fù)雜介質(zhì)的適用性，選擇一個三層模型(圖13)進(jìn)行實驗。該模型第一層為聲波介質(zhì)(水)，第二層為彈性各向同性介質(zhì)，第三層為VTI彈性介質(zhì)(VTI-2)，各層彈性參數(shù)見表1。震源位于(1000m，450m)；三個檢波器分別位于(500m，200m)、(500m，500m)、(500m，800m)，每層各有一個。圖14為三層模型2.5D數(shù)值解vx分量在不同時刻的波場快照，可以看出清晰的反射波、轉(zhuǎn)換波、透射波和多次波。由圖14a、圖14b可以看出，當(dāng)P波和S波到達(dá)固—液邊界時，透射和轉(zhuǎn)換后的P波出現(xiàn)在第一層(聲波介質(zhì))中，而反射和轉(zhuǎn)換后的P波和SV波在第二層(各向同性介質(zhì))中向下傳播；在第三層(VTI-2介質(zhì))中出現(xiàn)P波和S波穿過固—固界面產(chǎn)生的透射和轉(zhuǎn)換波。由圖14c～圖14f可以發(fā)現(xiàn)，當(dāng)波傳播到自由表面產(chǎn)生了反射的P波，在第二層介質(zhì)中出現(xiàn)了清晰的多次波。圖15、圖16、圖17分別為三個檢波器的2.5D與2D模擬記錄對比，可以看出，無論檢波器位于何種介質(zhì)中，二者存在明顯的振幅差異和相移。

圖13 三層模型示意圖星和三角分別表示震源和接收器點位置

圖14 三層模型2.5D數(shù)值解vx分量在不同時刻的波場快照(a)0.12s；(b)0.20s；(c)0.28s；(d)0.36s；(e)0.44s；(f)0.52s

圖15 三層模型第一個檢波器的2.5D與2D數(shù)值模擬記錄對比(a)vx分量；(b)vz分量；(c)圖a的歸一化顯示；(d)圖b的歸一化顯示。下圖為上圖方框的放大顯示

圖16 三層模型第二個檢波器的2.5D與2D數(shù)值模擬記錄對比(a)vx分量；(b)vz分量；(c)圖a的歸一化顯示；(d)圖b的歸一化顯示。下圖為上圖方框的放大顯示

圖17 三層模型第三個檢波器的2.5D與2D數(shù)值模擬記錄對比(a)vx分量；(b)vz分量；(c)圖a的歸一化顯示；(d)圖b的歸一化顯示。下圖為上圖方框的放大顯示

三層混合模型模擬結(jié)果表明，本文方法適用于復(fù)雜介質(zhì)2.5D地震波數(shù)值模擬。

6 結(jié)論

本文通過推導(dǎo)、整理得到了一個矩陣形式的2.5D廣義一階速度—應(yīng)力波動方程，該方程不僅適用于聲波介質(zhì)和一般各向異性介質(zhì)，而且適用于不同的邊界條件(聲波自由地表、固體自由地表和固—液邊界條件)。通過定義不同的波場矢量，并根據(jù)模型離散點介質(zhì)參數(shù)的變化給出了三個系數(shù)矩陣，在壓力或應(yīng)力矢量上施加點源后采用數(shù)值方法求解，可以用一個計算機(jī)程序?qū)崿F(xiàn)復(fù)雜介質(zhì)中的2D或2.5D地震波場數(shù)值模擬。為了更好地貼合起伏地形，本文選擇曲線網(wǎng)格有限差分法求解波動方程，數(shù)值模擬實驗表明，在不同的介質(zhì)中本文方法模擬得到的2.5D數(shù)值解與3D解析解非常匹配，并且適用于不同的邊界條件。除此之外，本文驗證了2.5D數(shù)值方法的計算效率遠(yuǎn)優(yōu)于3D數(shù)值方法。這種2.5D數(shù)值模擬方法可以直接用于實際點源地震數(shù)據(jù)的逆時偏移成像。

附錄A 彈性各向異性介質(zhì)中波動方 程的系數(shù)矩陣

將式(3)、式(4)寫成矩陣形式，則式(6)中的系數(shù)矩陣可表示為

(A-1)

(A-2)

(A-3)

式中

(A-4)

(A-5)

(A-6)

(A-7)

(A-8)

(A-9)

(A-10)

(A-11)

(A-12)

(A-13)

(A-14)

附錄B 聲波介質(zhì)中波動方程的系數(shù) 矩陣

將式(12)和式(13)寫成矩陣形式，則式(15)的系數(shù)矩陣為

(B-1)

(B-2)

(B-3)

附錄C 固體自由地表中波動方程的 系數(shù)矩陣

將式(31)和式(32)寫成矩陣形式，則式(34)系數(shù)矩陣可表示為

(C-1)

(C-2)

(C-3)

式中：

(C-4)

(C-5)

(C-6)

(C-7)

(C-8)

(C-9)

附錄D 固—液邊界中波動方程的系 數(shù)矩陣

將式(45)～式(47)聯(lián)合寫成一個矩陣形式，可導(dǎo)出式(49)的系數(shù)矩陣，為

(D-1)

(D-2)

(D-3)

式中：

(D-4)

(D-5)

(D-6)

(D-7)

(D-8)

(D-9)