王 輝 何兵壽* 邵祥奇

(①中國海洋大學(xué)海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東青島 266100；②青島海洋科學(xué)與技術(shù)國家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測技術(shù)功能實(shí)驗(yàn)室，山東青島 266100；③山東省煙臺(tái)市自然資源和規(guī)劃局，山東煙臺(tái) 264003)

0 引言

彈性波正演是研究縱橫波傳播規(guī)律的重要工具之一，在多波地震的資料采集、處理、解釋和反演等領(lǐng)域發(fā)揮著重要作用。彈性波方程及其模擬算法是多年來業(yè)界廣泛關(guān)注的熱點(diǎn)。

近年來，在彈性波方程的有限差分?jǐn)?shù)值模擬領(lǐng)域開展了大量的研究工作，取得了許多成果，包括彈性波偏微分方程組中空間偏導(dǎo)數(shù)的高階差分格式[1-5]、時(shí)間偏導(dǎo)數(shù)的高階差分格式[6-9]、各種差分格式的穩(wěn)定性條件[10-14]、數(shù)值頻散的壓制方法[15-19]、邊界條件[20-23]以及計(jì)算效率提升方法[24-25]等。這些成果對(duì)于推動(dòng)多波地震技術(shù)的發(fā)展具有重要意義。但總體而言，現(xiàn)有的彈性波正演技術(shù)仍存在諸多不足。①模擬精度與效率很難兼顧。小網(wǎng)格、高階差分格式可以提高模擬精度，但意味著計(jì)算量的指數(shù)級(jí)增加，導(dǎo)致模擬效率降低。雖然近似解析離散[26-27]等方法可以適應(yīng)大網(wǎng)格的差分計(jì)算，但其計(jì)算量仍遠(yuǎn)大于預(yù)期。②基于長方體網(wǎng)格的空間剖分方法無法準(zhǔn)確擬合復(fù)雜界面的形態(tài)，導(dǎo)致模擬結(jié)果中出現(xiàn)大量繞射波。近年出現(xiàn)的不規(guī)則網(wǎng)格技術(shù)[4，28-29]雖然緩解了這一問題，但距離根本解決還為時(shí)尚早；③現(xiàn)有算法能準(zhǔn)確模擬彈性波的傳播過程，模擬得到的三個(gè)分量均同時(shí)包含縱波和橫波[30]，必須借助額外的波場解耦算子[31-32]才能得到純縱波和純橫波的模擬結(jié)果，即模擬的縱、橫波記錄和快照同時(shí)受制于差分算法的精度和解耦算子的精度，增加了誤差來源[33-34]。

為解決上述第三個(gè)問題，本文通過對(duì)一階速度—脹縮—旋轉(zhuǎn)彈性波方程的差分離散實(shí)現(xiàn)彈性波的模擬。首先推導(dǎo)了一階速度—脹縮—旋轉(zhuǎn)彈性波方程在三維交錯(cuò)網(wǎng)格空間中的高階有限差分格式，給出了相應(yīng)的穩(wěn)定性條件；然后導(dǎo)出了適應(yīng)該方程的PML吸收邊界條件，在此基礎(chǔ)上實(shí)現(xiàn)了一階速度—脹縮—旋轉(zhuǎn)彈性波方程的正演模擬。本文的彈性波模擬方法的優(yōu)勢在于，不僅可以得到x、y、z三分量地震記錄，而且無需顯式地進(jìn)行縱、橫波解耦就可以得到橫波振動(dòng)速度矢量和縱波振動(dòng)速度矢量，實(shí)現(xiàn)了縱、橫波的保幅解耦。同時(shí)，數(shù)值求解該方程還可以直接得到各質(zhì)點(diǎn)的體應(yīng)變和旋轉(zhuǎn)矢量，為分析不同速度模型條件下縱、橫波的傳播規(guī)律提供更詳細(xì)的信息。

1 速度—脹縮—旋轉(zhuǎn)彈性波方程

各向同性介質(zhì)中的三維一階速度—脹縮—旋轉(zhuǎn)彈性波方程[35]為

(1)

式中：v=(vx，vy，vz)為質(zhì)點(diǎn)振動(dòng)速度矢量；vS=(vSx，vSy，vSz)為質(zhì)點(diǎn)的橫波振動(dòng)速度矢量；vP=(vPx，vPy，vPz)為質(zhì)點(diǎn)的縱波振動(dòng)速度矢量；cP、cS分別為縱、橫波傳播速度；θ為體應(yīng)變；ω=(ωx，ωy，ωz)為旋轉(zhuǎn)矢量。

式(1)不僅顯式地包含了常規(guī)彈性波方程中的質(zhì)點(diǎn)總振動(dòng)速度矢量，而且包含了由脹縮運(yùn)動(dòng)和剪切運(yùn)動(dòng)引起的質(zhì)點(diǎn)振動(dòng)速度矢量，同時(shí)還顯式包含了θ和ω，因此利用該方程進(jìn)行彈性波正演模擬不僅可以得到x、y、z三分量地震記錄，而且無需顯式地進(jìn)行縱、橫波解耦就可以得到縱波振動(dòng)速度矢量和橫波振動(dòng)速度矢量。同時(shí)，數(shù)值求解該方程還可以直接得到各質(zhì)點(diǎn)的體應(yīng)變和旋轉(zhuǎn)矢量，體應(yīng)變一般只與縱波有關(guān)，旋轉(zhuǎn)矢量一般只與橫波有關(guān)，這為分析不同速度模型條件下縱、橫波的傳播規(guī)律提供了更詳細(xì)的信息。

此外，由于目前的多波地震采集基本都采用縱波源激發(fā)，式(1)顯然十分便于解決彈性波正演中的縱波源設(shè)置問題，模擬時(shí)只需將震源函數(shù)加載到體應(yīng)變?chǔ)壬霞纯蓪?shí)現(xiàn)縱波源激發(fā)。而采用傳統(tǒng)的一階速度—應(yīng)力方程進(jìn)行彈性波模擬時(shí)，縱波源必須加載到三個(gè)正應(yīng)力分量上[9，13]，但縱波源的總能量按照什么樣的比例分配到三個(gè)正應(yīng)力分量上，目前并沒有一個(gè)公認(rèn)的結(jié)論。

2 交錯(cuò)網(wǎng)格數(shù)值模擬方法

2.1 交錯(cuò)網(wǎng)格高階差分格式

三維情況下，在圖1所示的交錯(cuò)網(wǎng)格空間[3-4]中求解式(1)，其中各波場分量在交錯(cuò)網(wǎng)格空間中的配置由表1給出。

表1 各波場分量在交錯(cuò)網(wǎng)格中的分布

圖1 三維交錯(cuò)網(wǎng)格示意圖

采用高階有限差分算法[6，8-9]對(duì)式(1)進(jìn)行差分離散，可得到三維一階速度—脹縮—旋轉(zhuǎn)彈性波方程正演模擬的時(shí)間2階、空間2N階有限差分格式。以θ分量為例，可表示為

--------------------

(2)

式中：Δt為時(shí)間延拓步長；Δx、Δy、Δz分別為直角坐標(biāo)系中x、y、z方向上的網(wǎng)格剖分間距；n為時(shí)間離散序號(hào)；i、j、k分別代表x、y、z方向上的離散網(wǎng)格點(diǎn)序號(hào)；Cm為2N階精度的交錯(cuò)網(wǎng)格差分系數(shù)，

其計(jì)算公式為[4，9，13]

(3)

表2為由式(3)得到的12階以內(nèi)各階差分的系數(shù)，一般來說，差分階數(shù)越高，模擬結(jié)果的精度越高。圖2為均勻介質(zhì)模型的空間2階和12階有限差分模擬的vx分量快照，其中縱、橫波速度分別為2800、1620m/s，空間三個(gè)方向的剖分間距均為5m，時(shí)間步長為0.5ms，縱波源是主頻為30Hz的Ricker子波，位于(500m，500m，500m)處。顯然，12階精度的差分格式能獲得更精確的模擬結(jié)果。雖然差分階數(shù)越高模擬精度越高[5]，但差分階數(shù)的提高也意味著計(jì)算量的迅速增大。數(shù)值實(shí)驗(yàn)表明，當(dāng)差分階數(shù)增高到一定程度后，繼續(xù)增高對(duì)模擬精度的提升十分有限。本文綜合考慮模擬精度和計(jì)算效率兩個(gè)因素后，采用12階精度的差分格式實(shí)現(xiàn)式(1)的數(shù)值模擬。

表2 交錯(cuò)網(wǎng)格不同階次有限差分系數(shù)

圖2 均勻模型2階(a)和12階(b)有限差分模擬的波場快照

2.2 穩(wěn)定性條件

采用與文獻(xiàn)[13，36]相同的推導(dǎo)方法可得式(2)的穩(wěn)定性條件為

(4)

式中vmax為模型中的最大速度。

2.3 PML吸收邊界條件

本文引入地震波方程正演領(lǐng)域常用的PML邊界條件[21-23]吸收入射到模型邊界的外行波。采用圖3所示的鑲邊方案在模型周邊鑲嵌PML吸收層，圖3中區(qū)域Ⅰ為8個(gè)角點(diǎn)區(qū)，區(qū)域Ⅱ?yàn)?2條棱邊區(qū)，區(qū)域Ⅲ為剩余的6個(gè)面。

圖3 三維空間PML鑲邊方案示意圖

以θ分量為例，在鑲邊區(qū)域?qū)ⅵ确至糠纸鉃棣葂、θy和θz三個(gè)分量，分別計(jì)算各分量，再求和得到鑲邊區(qū)域的θ分量，即

(5)

式中dx、dy、dz分別為x、y、z三個(gè)方向上的衰減因子，其計(jì)算公式為[9，22]

(6)

式中：V為波速，一般取縱波速度cP；R為理論邊界反射系數(shù)，一般取經(jīng)驗(yàn)值0.0001[9]；δ為PML鑲邊層厚度；hx、hy、hz分別為鑲邊層中的各網(wǎng)格點(diǎn)與有效計(jì)算區(qū)域邊界的最短距離。

同樣在交錯(cuò)網(wǎng)格空間中對(duì)式(5)進(jìn)行差分離散，可得到式(5)的高階有限差分格式

(7)

同理可得出其他波場分量的PML吸收邊界及其高階差分格式。對(duì)比均勻模型加與不加PML吸收邊界的波場快照(圖4)可見，本文的PML邊界條件幾乎能完全吸收入射到截?cái)噙吔绲耐庑胁ā?/p>

圖4 不加(a)與加(b)PML吸收邊界的vx分量快照對(duì)比

3 模擬快照分析

對(duì)一個(gè)簡單的兩層水平層狀介質(zhì)模型分別采用一階速度—應(yīng)力彈性波方程[9]和一階速度—脹縮—旋轉(zhuǎn)彈性波方程[32]進(jìn)行數(shù)值模擬，利用模擬快照分析兩個(gè)方程的聯(lián)系和區(qū)別。模型尺寸為1000m×600m×1000m，上層縱、橫波速度分別為2250、1300m/s，下層縱、橫波速度分別為2900、1676m/s，界面埋深為350m，空間網(wǎng)格尺寸為5m×5m×5m，時(shí)間采樣間隔為0.5ms。兩個(gè)方程均采用空間12階精度的有限差分算法進(jìn)行模擬。經(jīng)驗(yàn)證，以上模擬參數(shù)同時(shí)滿足兩個(gè)方程的穩(wěn)定性條件。均采用主頻為30Hz的Ricker子波作為縱波激發(fā)源，炮點(diǎn)位于(500m，300m，0)處。對(duì)一階速度—應(yīng)力方程模擬時(shí)，縱波源加載到三個(gè)正應(yīng)力分量上，而對(duì)一階速度—脹縮—旋轉(zhuǎn)方程進(jìn)行模擬時(shí)，由于該方程沒有顯式包含應(yīng)力分量，故將縱波源加載到體應(yīng)變?chǔ)壬?。兩種算法模擬所需的計(jì)算量如表3所示(一階速度—應(yīng)力方程的模擬過程包含了基于散度和旋度算子的縱橫波分離運(yùn)算)，顯然，相同精度條件下，本文算法的運(yùn)算量低于一階速度—應(yīng)力彈性波方程數(shù)值模擬的運(yùn)算量。

表3 兩種方法計(jì)算量對(duì)比

圖5為0.35s時(shí)刻兩種方程模擬的vx、vy、vz分量快照，是野外采集中可以直接用速度檢波器接收并記錄的信息。由于波的傳播方向并非與坐標(biāo)軸完全平行，因此三個(gè)分量中均同時(shí)包含縱波與橫波。顯然，由兩個(gè)方程得到的質(zhì)點(diǎn)的振動(dòng)速度是完全一致的(因?yàn)槎叩奈锢硪饬x本就是相同的)。這說明，如果以得到三分量炮記錄為正演目標(biāo)，則式(1)可以獲得與一階速度—應(yīng)力方程完全相同的結(jié)果。

圖5 水平層狀模型兩種方程模擬的波前快照對(duì)比(一)

圖6a～圖6d為式(1)正演得到的ω的三分量快照和θ的快照。依據(jù)式(1)，θ和ω對(duì)時(shí)間的偏導(dǎo)數(shù)分別對(duì)應(yīng)v的散度和旋度，分別指示縱波和橫波。圖6e～圖6g為直接對(duì)圖5d～圖5f求旋度和散度的結(jié)果，它們幾乎與ω和θ的快照完全相同，表明θ和ω能代替一階速度—應(yīng)力彈性波方程中質(zhì)點(diǎn)振動(dòng)速度的散度和旋度。需要說明的是，雖然業(yè)界常通過直接求取v的散度和旋度分離縱、橫波，但這種方法得到的結(jié)果只在波形上與真實(shí)的縱、橫波相似，卻沒有明確的物理意義，因?yàn)橹挥挟?dāng)Helmholtz分解的對(duì)象是位移矢量時(shí)，分離結(jié)果才有明確的物理意義，因此圖6a～圖6d具有明確的物理意義，而圖6e～圖6g則不然。此外，對(duì)比圖6與圖5還可以看出，上述兩種方法得到的橫波的x分量和y分量與分離前混合波場中橫波的x分量和y分量存在90°的相位差，其極性反轉(zhuǎn)位置沒有與分離前混合波場中的各分量對(duì)應(yīng)。而分離后的橫波z分量為零，這顯然與實(shí)際情況不符。這說明基于散度和旋度算子的波場分離方法不能實(shí)現(xiàn)縱、橫波的保幅分離。圖6c和圖6g為零的原因是：在進(jìn)行旋度運(yùn)算時(shí)，ωz是vy在x方向的偏導(dǎo)數(shù)與vx在y方向的偏導(dǎo)數(shù)之差，由于本算例所用的模型在z方向上水平分層，vy在x方向的偏導(dǎo)與vx在y方向的偏導(dǎo)相等，故旋度的z分量為零。

圖6 水平層狀模型兩種方程模擬的波前快照對(duì)比(二)

圖7為在式(1)方程的正演直接得到的純vP和純vS的三分量快照，與圖5對(duì)比可見，圖7不僅實(shí)現(xiàn)了縱、橫波的完全分離，并且分離后各分量的縱、橫波的極性與混合波場中的極性完全對(duì)應(yīng)，相位也保持一致。同時(shí)，依據(jù)式(1)，由于混合波場v是兩個(gè)矢量波場vS和vP的簡單相加，因此利用式(1)得到的縱、橫波分離結(jié)果一定是保幅的。圖8為水平位置(250m，0)處vx、vPx、vSx分量垂向波形對(duì)比，證明了一階速度—脹縮—旋轉(zhuǎn)彈性波方程能夠?qū)崿F(xiàn)縱橫波的保幅分離。與現(xiàn)有的縱、橫波保幅分離算法相比，本文算法的優(yōu)勢在于：一階速度—脹縮—旋轉(zhuǎn)彈性波方程無需任何額外計(jì)算即可實(shí)現(xiàn)縱、橫波的完全分離，分離精度只與模擬算法的精度有關(guān)，而現(xiàn)有的縱、橫波保幅分離算法的精度同時(shí)受制于模擬算法的精度和解耦算子的精度。

圖7 水平層狀模型速度—脹縮—旋轉(zhuǎn)彈性波方程模擬的純vP、純vS各分量快照

圖8 本文方法模擬的vx、vPx、vSx分量垂向波形對(duì)比

4 數(shù)值算例

用SEG/EAGE鹽丘模型檢驗(yàn)本文模擬算法的正確性，并進(jìn)一步分析一階速度—脹縮—旋轉(zhuǎn)方程在彈性波模擬與應(yīng)用方面的優(yōu)勢。模型尺寸為1000m×1000m×1000m，縱、橫波速度模型如圖9所示，空間網(wǎng)格尺寸為5m×5m×5m，時(shí)間步長為0.35ms，模擬參數(shù)滿足穩(wěn)定性條件?？v波震源是主頻為30Hz的Ricker子波，炮點(diǎn)坐標(biāo)為(500m，500m，0)。地震記錄長度為1.05s，8條多分量接收線位于地表，間距為125m，第1、第8條分別位于y=0、y=875m處，每條線200道，道間距為5m，單炮共有1600道。

圖9 三維鹽丘模型

圖10、圖11為正演得到的各分量單炮記錄。由于該模型包含了斷層、凹陷、隆起、傾斜界面和水平界面，使得彈性波的傳播過程極為復(fù)雜且混合波場中縱波與橫波耦合在一起，難以區(qū)分，不便于分析縱、橫波的傳播過程，也不便于評(píng)價(jià)觀測系統(tǒng)的合理性。而一階速度—脹縮—旋轉(zhuǎn)方程可以在波場延拓中直接獲得vS和vP的各個(gè)分量，使混合在一起的波場得到分離，為分別研究反射縱波或轉(zhuǎn)換橫波的傳播過程和機(jī)理提供了更準(zhǔn)確的信息。同時(shí)，上述結(jié)果也說明一階速度—脹縮—旋轉(zhuǎn)彈性波方程還有利于在延拓過程中直接求取純縱波或純橫波的某些屬性，如縱波的傳播方向、橫波的傳播方向及其極化方向等，而這些屬性信息又往往是偏移和反演等環(huán)節(jié)的關(guān)鍵信息。

圖10 部分三維鹽丘模型本文方法模擬的混合波場地震記錄

圖11 部分三維鹽丘模型本文方法模擬的純波地震記錄

5 結(jié)論與討論

本文給出了三維各向同性介質(zhì)中一階速度—脹縮—旋轉(zhuǎn)彈性波方程的高階有限差分格式及其穩(wěn)定性條件，以此為基礎(chǔ)實(shí)現(xiàn)了該方程的數(shù)值模擬，并利用模型模擬結(jié)果詳細(xì)分析了該方程的優(yōu)勢，得出如下結(jié)論。

(1)本文算法能實(shí)現(xiàn)一階速度—脹縮—旋轉(zhuǎn)彈性波方程的高精度數(shù)值模擬。

(2)除能獲得常規(guī)彈性波模擬中的三分量合成記錄和波場快照外，一階速度—脹縮—旋轉(zhuǎn)彈性波方程還能夠在不增加任何額外運(yùn)算的前提下直接實(shí)現(xiàn)彈性波的保幅解耦。解耦結(jié)果是一個(gè)縱波速度矢量和一個(gè)橫波速度矢量，除了便于分析縱、橫波在各個(gè)方向的振動(dòng)情況外，還特別有利于求取純縱波或純橫波的傳播方向和極化方向等屬性信息。同時(shí)，由于解耦結(jié)果是兩個(gè)矢量，避免了彈性波逆時(shí)偏移中進(jìn)行互相關(guān)成像時(shí)矢量和標(biāo)量無法相關(guān)的難題。

(3)對(duì)一階速度—脹縮—旋轉(zhuǎn)彈性波方程進(jìn)行數(shù)值模擬還可以直接得到一個(gè)體應(yīng)變?chǔ)群托D(zhuǎn)矢量ω，其物理意義與常規(guī)彈性波方程中位移矢量的散度和旋度算子的物理意義相同。θ與ω雖然與特定模型條件下的縱橫波傳播過程有關(guān)，但不能準(zhǔn)確描述縱、橫波傳播過程。

(4)本文模擬得到的vS和vP都是矢量，而事實(shí)上，各向同性介質(zhì)中的縱波是標(biāo)量，雖然vP可以準(zhǔn)確描述由脹縮運(yùn)動(dòng)導(dǎo)致的質(zhì)點(diǎn)振動(dòng)在各個(gè)方向的投影情況，但由于業(yè)界習(xí)慣于將縱波當(dāng)成一個(gè)獨(dú)立標(biāo)量處理和分析，因此將vP合成為一個(gè)標(biāo)量更易為業(yè)界接受，現(xiàn)有的矢量波場的標(biāo)量化處理技術(shù)完全可以解決這一問題。