侯志強(qiáng), 尹文筍, 李 鍵, 孫永壯, 劉 云

(中國(guó)海洋石油集團(tuán)有限公司上海分公司，上海 200335)

基于彈性理論的反射波地震勘探技術(shù)在能源、資源勘探開發(fā)和環(huán)境調(diào)查等領(lǐng)域發(fā)揮了重要作用，但目前隨著人們對(duì)深部能源勘探開發(fā)問題的日益重視和對(duì)地球深部地層結(jié)構(gòu)和巖性等問題的持續(xù)關(guān)注，業(yè)界都對(duì)深部地層的地震勘探精度提出了更高要求。由于地下巖石普遍具有粘彈性質(zhì)，常規(guī)基于彈性假設(shè)的地震勘探技術(shù)在解決深部地層的勘探問題時(shí)往往表現(xiàn)出不適應(yīng)性，表現(xiàn)在：(1)對(duì)于深部地層，由于地震反射波的傳播路徑很長(zhǎng)，介質(zhì)粘滯性對(duì)地震波傳播影響的累積效應(yīng)增大，地震波的實(shí)際傳播規(guī)律與彈性介質(zhì)假設(shè)情況嚴(yán)重不符，基于彈性介質(zhì)假設(shè)的地震資料處理與反演技術(shù)很難取得滿意處理效果；(2)巖石粘滯性造成深部地層反射波的高頻成分衰減嚴(yán)重，深部地層的地震反射波高頻成分缺失，頻帶變窄，倍頻程減小，分辨率降低；(3)巖石粘滯性造成深部地層的反射波能量減弱，信噪比降低，增加了成像難度；(4)地層粘滯性使得不同深度地層反射波的頻譜不一致，導(dǎo)致部分處理流程的參數(shù)選擇困難，同時(shí)將成像結(jié)果的垂向分辨率變成t0時(shí)的函數(shù)，導(dǎo)致地下淺、中、深層具有不同的分辨率，增加了解釋和反演難度。

研究并利用基于粘彈假設(shè)的地震勘探理論和方法可以克服上述缺陷，更好地解決深部地層的精確勘探問題。所謂粘彈介質(zhì)是指力學(xué)性質(zhì)介于完全彈性和完全粘性之間的介質(zhì)，這種介質(zhì)在外力作用下會(huì)同時(shí)表現(xiàn)出彈性性質(zhì)和粘性性質(zhì)，故地震波在這種介質(zhì)中傳播時(shí)具有不同于彈性介質(zhì)的傳播機(jī)理，這一獨(dú)特傳播機(jī)理是研發(fā)基于粘彈假設(shè)的地震勘探技術(shù)的理論基礎(chǔ)，因此研究粘彈介質(zhì)中的地震波傳播機(jī)理對(duì)于解決深部地層的精確勘探問題具有重要意義。

業(yè)界對(duì)粘彈理論的研究始于1945年的Stokes粘彈地震波動(dòng)方程[1]，該方程主要考慮了由質(zhì)點(diǎn)內(nèi)摩擦引起的地震波能量耗損，之后粘彈地震波傳播理論得到了快速發(fā)展[2-11]，目前已形成了Kelvin模型、Maxwell介質(zhì)模型、標(biāo)準(zhǔn)線性介質(zhì)和達(dá)朗貝爾模型等具有不同粘彈性質(zhì)的連續(xù)介質(zhì)力學(xué)模型，其中Kelvin模型是目前地震勘探領(lǐng)域應(yīng)用最多的粘彈性介質(zhì)模型。

上述模型對(duì)應(yīng)的地震波方程的解析解或數(shù)值解是研究不同類型粘彈介質(zhì)中地震波傳播規(guī)律的重要基礎(chǔ)，復(fù)雜模型條件下地震波方程解析解的求取極為困難，故業(yè)界往往采用波動(dòng)方程的數(shù)值解來研究波傳播機(jī)理或解決實(shí)際問題，粘彈介質(zhì)中地震波方程的數(shù)值求解是粘彈地震理論與方法的重要研究?jī)?nèi)容。

目前，業(yè)界用于求取地震波方程數(shù)值解的算法主要包括反射率法[3-4]、有限元法[10]、虛譜法[6]、有限差分法[5,8-9,11]等，其中有限差分法由于具有計(jì)算速度快、精度高、易實(shí)現(xiàn)等優(yōu)點(diǎn)而得到了廣泛應(yīng)用。本文內(nèi)容屬粘彈性波方程的有限差分?jǐn)?shù)值求解范疇，首先基于交錯(cuò)網(wǎng)格技術(shù)[12-13]給出了基于Kelvin模型中的粘彈性波方程的高階有限差分格式、穩(wěn)定性條件和吸收邊界條件。其次針對(duì)粘彈介質(zhì)中縱橫波的解耦問題，本文借鑒完全彈性介質(zhì)中縱橫波的解耦技術(shù)[14]，從散度算子和旋度算子出發(fā)，通過理論分析給出了一種粘彈介質(zhì)中縱橫波的保幅分離算法，將粘彈介質(zhì)中的矢量彈性波場(chǎng)分解為矢量縱波場(chǎng)和矢量橫波場(chǎng)，再依據(jù)縱橫波傳播方向與偏振方向之間的關(guān)系，將三維矢量縱波與矢量橫波合成為標(biāo)量縱波與標(biāo)量橫波，以獲取具有實(shí)際意義的縱、橫波波場(chǎng)和單炮記錄，實(shí)現(xiàn)了基于縱橫波保幅分離的粘滯介質(zhì)彈性波正演模擬。

1 粘彈介質(zhì)中的彈性波方程及有限差分求解方法

1.1 三維粘彈介質(zhì)中的彈性波方程

三維Kelvin模型中的粘彈性波方程為：

(1)

其中：x,y,z為三個(gè)直角坐標(biāo)；t為時(shí)間；ρ為密度；vx、vy、vz分別為x,y,z三個(gè)方向的質(zhì)點(diǎn)震動(dòng)速度分量；σxx,σyyσzz,σxy,σxz,σyz為應(yīng)力分量；cp為縱波速度；cs為橫波速度；Qp為縱波品質(zhì)因子；Qs為橫波品質(zhì)因子；ω為圓頻率。

1.2 三維粘彈介質(zhì)中的彈性波方程的交錯(cuò)網(wǎng)格有限差分解法

交錯(cuò)網(wǎng)格法[12-13]是指在常規(guī)網(wǎng)格中引入半網(wǎng)格點(diǎn)，在半網(wǎng)格點(diǎn)上進(jìn)行空間導(dǎo)數(shù)的計(jì)算，把應(yīng)力分量和速度分量定義在兩套網(wǎng)格上。與常規(guī)網(wǎng)格相比，交錯(cuò)網(wǎng)格能夠有效解決一階彈性波方程速度分量和應(yīng)力分量的耦合關(guān)系，在不增加計(jì)算量的前提下提高計(jì)算精度和穩(wěn)定性。以式(1)中的σxx分量和vx分量為例，它在交錯(cuò)網(wǎng)格空間中的高階有限差分格式如式(2)、(3)所示：其他分量的差分格式可用類似方法導(dǎo)出。

(2)

(3)

(4)

差分計(jì)算方法為：

(5)

式(2)、(3)、(5)的穩(wěn)定性條件為：

(6)

其中cmax為模型中最大縱波速度。

1.3 吸收邊界條件

采用PML邊界條件[15-16]解決式(1)求解過程中的截?cái)噙吔鐔栴}。PML吸收邊界的基本思想是在計(jì)算區(qū)域增加吸收層，在吸收層內(nèi)設(shè)置衰減因子對(duì)波場(chǎng)進(jìn)行衰減。對(duì)計(jì)算區(qū)域鑲邊后的三維空間如圖1所示。以vx分量為例，依據(jù)PML的分裂思路[15-16]，可將其分解為x,y,z三個(gè)方向的分量vx_x，vx_y，vx_z，即：

vx=vx_x+vx_y+vx_z。

(7)

各分量的吸收方程如下：

(8)

其中d(x)、d(y)、d(z)分別為x、y、z三個(gè)方向上的衰減因子，其取值詳見文獻(xiàn)[15]。

在不同的邊界對(duì)不同的分量進(jìn)行吸收即可壓制截?cái)噙吔绲膫畏瓷洹Ｈ砸詖x分量為例，在圖1所示的三維吸收邊界示意圖中，各個(gè)邊界區(qū)域的衰減因子如下：

圖1 三維空間PML吸收邊界示意圖

2 縱橫波保幅解耦方法與矢量縱橫波的標(biāo)量合成

2.1 縱橫波保幅解耦方法

式(1)中vx、vy、vz的本質(zhì)是質(zhì)點(diǎn)的振動(dòng)速度矢量v在直角坐標(biāo)系三個(gè)坐標(biāo)軸上的投影，由于縱波與橫波均可引起這三個(gè)方向上的質(zhì)點(diǎn)振動(dòng)，因此vx、vy、vz分量都同時(shí)包含縱波與橫波，兩種波耦合在一起不便于分析縱橫波的傳播與衰減機(jī)理，因此有必要在粘彈介質(zhì)彈性波方程正演的過程中采用適當(dāng)方法對(duì)縱橫波進(jìn)行解耦，以得到合成縱波記錄和橫波記錄。

各向同性介質(zhì)中縱波是無旋場(chǎng)，橫波是無散場(chǎng)，因此可通過求取彈性波場(chǎng)的散度與旋度得到縱波場(chǎng)與橫波場(chǎng)，Aki和Richards以此為基礎(chǔ)，利用v的散度和旋度算子實(shí)現(xiàn)了彈性介質(zhì)中的縱橫波分離[17]，這種方法實(shí)現(xiàn)簡(jiǎn)單，計(jì)算量小，但會(huì)使波場(chǎng)的相位和振幅信息產(chǎn)生畸變[18-19]，且分離后的波場(chǎng)在極性反轉(zhuǎn)位置上無法與分離前混合波場(chǎng)各分量對(duì)應(yīng)，因此基于散度和旋度算子的波場(chǎng)解耦方法是不保幅的。

由于本文研究的Kelvin模型仍屬于各向同性介質(zhì)范疇，因此彈性各向同性介質(zhì)中基于散度與旋度算子的波場(chǎng)解耦方法同樣無法解決Kelvin模型的縱橫波保幅分離問題。本文的主要目標(biāo)就是研究新的方法實(shí)現(xiàn)Kelvin粘彈模型的縱橫波保幅分離。

假設(shè)Kelvin粘彈模型中的矢量波場(chǎng)v由vp和vs兩個(gè)矢量場(chǎng)組成：

v=vp+vs。

(9)

其中vp為由縱波引起的質(zhì)點(diǎn)振動(dòng)速度矢量；vs為由橫波引起的振動(dòng)速度矢量，這兩個(gè)矢量在笛卡爾坐標(biāo)系中的表達(dá)形式為vp=(vp_x,vp_y,vp_z)，vs=(vs_x,vs_y,vs_z)。

對(duì)式(9)分別求散度和旋度可得：

(10)

由于縱波為無旋場(chǎng)，橫波為無散場(chǎng)，有：

(11)

(7)式可寫為：

(12)

將上式變換到波數(shù)域，有：

(13)

縱橫波的波數(shù)域單位傳播矢量Kp、Ks與縱橫波速度及圓頻率ω之間滿足以下關(guān)系：

(14)

聯(lián)立式(10)和式(11)可得：

(15)

對(duì)上式作傅里葉反變換，可得

(16)

上式即為各向同性粘彈介質(zhì)中的縱橫波保幅分離算子，它可以在時(shí)間空間域利用有限差分來實(shí)現(xiàn)，具體計(jì)算公式為：

(17)

和Aki等的方法[17]相比，本文方法的優(yōu)勢(shì)在于：不會(huì)引起縱橫波振幅和相位畸變，且將縱波當(dāng)作矢量進(jìn)行處理，在物理意義和波場(chǎng)的極性反轉(zhuǎn)位置上能與分離前混合波場(chǎng)各分量對(duì)應(yīng)；和李振春等的方法[20]相比，本文方法的優(yōu)勢(shì)在于：首先利用地層中縱橫波的傳播速度對(duì)分離后的矢量波場(chǎng)進(jìn)行振幅補(bǔ)償，對(duì)補(bǔ)償結(jié)果再沿時(shí)間方向進(jìn)行積分使得分離結(jié)果更具保真性。

利用圖2a所示的水平層狀模型驗(yàn)證本文算法的保幅性，兩層介質(zhì)的縱波速度分別為2 500、3 000 m/s，橫波速度分別為1 443、1 764 m/s，密度分別為2 000、2 500 kg/m3，縱、橫波品質(zhì)因子為常數(shù)80，界面埋深400 m。波場(chǎng)模擬所用的參數(shù)為：震源為主頻f0=35 Hz的雷克子波，震源置于地表，其水平位置為(500 m，250 m)，空間網(wǎng)格大小5 m×5 m，時(shí)間步長(zhǎng)0.5 ms。由于本文是在時(shí)間域?qū)κ?1)進(jìn)行求解，故假定圓頻率ω為常數(shù)，其值為2πf0。圖2b～2d為正演過程中記錄t=350 ms時(shí)的三分量波場(chǎng)快照，由圖可見，每個(gè)分量快照中都同時(shí)包含縱波與橫波，縱橫波耦合在一起，互為串?dāng)_，必須將之分解為相對(duì)獨(dú)立的縱波分量和橫波分量才便于分析粘彈介質(zhì)中的縱橫波傳播規(guī)律。

圖2 模型示意圖及其三分量快照

在正演模擬過程中分別利用Aki等的方法[17]、李振春等的方法[20]和本文算法進(jìn)行波場(chǎng)解耦。圖3、4和5分別為上述三種方法的縱橫波解耦結(jié)果快照，由圖可見，這三種方法都能在波場(chǎng)模擬過程中實(shí)現(xiàn)縱橫波的解耦，其中Aki的方法[17]將彈性波場(chǎng)分解為標(biāo)量縱波和矢量橫波，后兩種方法則將彈性波場(chǎng)分解為矢量縱波和矢量橫波，由于標(biāo)量可以看作一種特殊的矢量，因此解耦結(jié)果中的縱波無論是標(biāo)量還是矢量，在理論上都是正確的。但圖3中矢量橫波各個(gè)分量中的極性反轉(zhuǎn)位置與原波場(chǎng)不一致，同時(shí)分離結(jié)果中橫波的z分量(見圖3d)在inline方向的波場(chǎng)值為零，這與vz分量中的橫波存在明顯差異，這些現(xiàn)象說明散度算子和旋度算子對(duì)橫波場(chǎng)具有改造作用，它無法得到地下真實(shí)的橫波場(chǎng)，只能得到改造后的橫波場(chǎng)。

圖4、圖5為采用后兩種方法得到的縱橫波分離快照，由圖可見，這兩種方法都能實(shí)現(xiàn)縱橫波的完全解耦，且解耦前后縱、橫波各分量的極性反轉(zhuǎn)位置能夠準(zhǔn)確對(duì)應(yīng)，這表明后兩種方法的解耦精度高于散度和旋度算子。但對(duì)于縱波的三個(gè)分量，李振春等的方法[20]得到的結(jié)果與分離前的數(shù)據(jù)存在90°的相位差，而本文方法與原始數(shù)據(jù)之間不存在相位變化，說明本文算法的解耦精度高于第二種方法。

圖3 散度和旋度算子的波場(chǎng)解耦快照

圖4 李振春等算法的波場(chǎng)解耦快照

圖5 本文算法的波場(chǎng)解耦快照

為證明本文算法的保幅優(yōu)勢(shì)，從不同方法分離前后的快照結(jié)果中選取一道數(shù)據(jù)進(jìn)行比較，該道在地表的投影位置為：(650，50)，圖6a為不同方法分離前后的縱波z分量顯示，其中第一道數(shù)據(jù)為波場(chǎng)分離前的z分量混合波場(chǎng)，第二道為利用散度算子[17]得到的標(biāo)量縱波，第三、第四道為分別利用李振春等[20]和本文算法得到的矢量縱波z分量，由圖可見，本文算法得到的縱波z分量結(jié)果與原波場(chǎng)中的縱波的振幅和相位完全一致，而另外兩種方法得到的結(jié)果的振幅比原波場(chǎng)小1～2個(gè)數(shù)量級(jí)，且存在90°相位差。圖6b為不同方法分離前后的橫波x分量顯示，其中第一道數(shù)據(jù)為場(chǎng)分離前的x分量混合波場(chǎng)，其余三道分別為利用旋度算子、李振春等的方法和本文算法得到的矢量橫波的x分量，三種方法對(duì)橫波的分量結(jié)果都不存在相位畸變，但前兩種分離算法分離結(jié)果的振幅比原始數(shù)據(jù)小一個(gè)數(shù)量級(jí)，本文方法的分離結(jié)果則與分離前的橫波完全一致。圖7為該位置處用不同方法得到的正演單道記錄的比較，分析該圖可以得到與圖6相同的結(jié)論，這表明本文給出的粘彈介質(zhì)縱橫波分離方法是保幅的。

2.2 矢量縱、橫波的標(biāo)量合成方法

本文的縱橫波保幅解耦方法能為研究粘彈介質(zhì)縱、橫波的傳播規(guī)律提供保真的數(shù)據(jù)，還能為基于點(diǎn)積互相關(guān)的彈性波逆時(shí)偏移成像[21]提供保真的矢量縱波和矢量橫波數(shù)據(jù)。但在常規(guī)逆時(shí)偏移技術(shù)[22-23]中往往需要用標(biāo)量的縱波與橫波進(jìn)行互相關(guān)成像，同時(shí)，工業(yè)界也傾向于利用更具明確地球物理意義的標(biāo)量縱波與標(biāo)量橫波來解決地質(zhì)問題，在這種情況下就需要對(duì)分離后的矢量縱波與矢量橫波進(jìn)行標(biāo)量合成。

矢量波場(chǎng)的標(biāo)量合成問題一般由振幅計(jì)算和極性確定兩部分組成。對(duì)于標(biāo)量波的振幅計(jì)算問題，不管是縱波還是橫波，都可以通過求取矢量波場(chǎng)的模來完成，問題的難點(diǎn)在于如何確定標(biāo)量化后的波場(chǎng)的極性。本文采用質(zhì)點(diǎn)振動(dòng)法求取標(biāo)量橫波的極性[24]，對(duì)于標(biāo)量縱波的極性問題，本文規(guī)定質(zhì)點(diǎn)振動(dòng)方向與z軸夾角小于90°時(shí)為負(fù)，反之為正，由此可以將縱波的極性求取問題轉(zhuǎn)換為質(zhì)點(diǎn)振動(dòng)方向的求取問題，由于縱波的傳播方向與質(zhì)點(diǎn)的振動(dòng)方向相同，因此可利用縱波的傳播方向確定質(zhì)點(diǎn)振動(dòng)方向，進(jìn)而確定標(biāo)量縱波的極性?？v波的傳播方向信息可利用坡印廷矢量得到，彈性波坡印廷矢量的求取方法已有多人做過研究[25-26]，本文不贅述。圖8為利用上述原理對(duì)圖5所示的波場(chǎng)快照進(jìn)行標(biāo)量合成的結(jié)果，標(biāo)量合成后的縱、橫波場(chǎng)具有更為明確的物理意義且更便于實(shí)際應(yīng)用。

圖6 不同分離方法得到的波場(chǎng)快照比較

圖7 不同分離方法得到的合成記錄單道比較

圖8 矢量縱、橫波標(biāo)量化后的快照

3 模型算例

模型的縱橫波速度如圖9所示，縱、橫波品質(zhì)因子均取常數(shù)80，正演所用的參數(shù)如下：模型大小為1 500 m×250 m×1 050 m，空間網(wǎng)格大小為5 m×5 m×5 m，采用間隔為0.35 ms，記錄長(zhǎng)度為1.05 s。采用Ricker子波作為震源，主頻為35 Hz，震源位于(760 m，250 m,0 m)處，三線接收，線距50 m，每條測(cè)線300道接收。圖10為基于本文算法的合成縱波記錄和轉(zhuǎn)換橫波單炮記錄，圖11為該模型完全彈性情況下的合成反射縱波記錄和轉(zhuǎn)換橫波記錄，對(duì)比圖10，11可以看出粘彈介質(zhì)情況下，由于受到地層的粘滯吸收作用，反射縱波和轉(zhuǎn)換橫波的能量均弱于彈性情況，圖12為第1條線100道縱、橫波分量中的375～900 ms時(shí)間段波場(chǎng)對(duì)比圖，圖中無論是縱波還是轉(zhuǎn)換橫波，介質(zhì)完全彈性情況下的振幅明顯高于粘彈情況，且隨著時(shí)間的增大，這種差別也逐漸增大，其原因?yàn)椋簜鞑r(shí)間的增大往往意味著傳播距離的增加，即傳播的波長(zhǎng)數(shù)增大，介質(zhì)粘滯性的累積效果增加。

圖9 縱橫波速度模型

圖10 粘彈條件下的合成縱橫波記錄

圖11 完全彈性條件下的合成縱橫波記錄

圖12 不同條件下第1條線100道地震記錄對(duì)比

4 結(jié)論

(1) 本文將標(biāo)量縱波看作一種特殊的矢量，推導(dǎo)了三維粘滯彈性波的縱橫波保幅解耦公式，給出了差分求解方法。本文算法的解耦結(jié)果能夠?qū)崿F(xiàn)與原三分量波場(chǎng)的對(duì)應(yīng)，解耦結(jié)果可方便的用于波場(chǎng)分析，且具有很高的保幅性。

(2) 本文算法解耦后的縱、橫波三分量數(shù)據(jù)可直接用于基于點(diǎn)積互相關(guān)的逆時(shí)偏移成像；縱、橫波的標(biāo)量合成結(jié)果可直接用于常規(guī)的彈性波逆時(shí)偏移成像。

(3) 本文基于縱橫波保幅分離的粘彈介質(zhì)彈性波正演模擬算法既可以獲得常規(guī)三分量合成地震記錄，也可以獲得波場(chǎng)解耦后的縱波合成記錄與轉(zhuǎn)換橫波合成記錄，還可以獲得三分量矢量縱波記錄和矢量橫波記錄。