基于縱橫波保幅分離的粘彈介質(zhì)彈性波正演模擬*

侯志強, 尹文筍, 李 鍵, 孫永壯, 劉 云

(中國海洋石油集團有限公司上海分公司，上海 200335)

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

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

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

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

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

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

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

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

(1)

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

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

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

(2)

(3)

(4)

差分計算方法為：

(5)

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

(6)

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

1.3 吸收邊界條件

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

vx=vx_x+vx_y+vx_z。

(7)

各分量的吸收方程如下：

(8)

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

在不同的邊界對不同的分量進行吸收即可壓制截斷邊界的偽反射。仍以vx分量為例，在圖1所示的三維吸收邊界示意圖中，各個邊界區(qū)域的衰減因子如下：

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

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

2.1 縱橫波保幅解耦方法

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

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

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

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

v=vp+vs。

(9)

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

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

(10)

由于縱波為無旋場，橫波為無散場，有：

(11)

(7)式可寫為：

(12)

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

(13)

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

(14)

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

(15)

對上式作傅里葉反變換，可得

(16)

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

(17)

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

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

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

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

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

圖3 散度和旋度算子的波場解耦快照

圖4 李振春等算法的波場解耦快照

圖5 本文算法的波場解耦快照

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

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

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

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

圖6 不同分離方法得到的波場快照比較

圖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，記錄長度為1.05 s。采用Ricker子波作為震源，主頻為35 Hz，震源位于(760 m，250 m,0 m)處，三線接收，線距50 m，每條測線300道接收。圖10為基于本文算法的合成縱波記錄和轉(zhuǎn)換橫波單炮記錄，圖11為該模型完全彈性情況下的合成反射縱波記錄和轉(zhuǎn)換橫波記錄，對比圖10，11可以看出粘彈介質(zhì)情況下，由于受到地層的粘滯吸收作用，反射縱波和轉(zhuǎn)換橫波的能量均弱于彈性情況，圖12為第1條線100道縱、橫波分量中的375～900 ms時間段波場對比圖，圖中無論是縱波還是轉(zhuǎn)換橫波，介質(zhì)完全彈性情況下的振幅明顯高于粘彈情況，且隨著時間的增大，這種差別也逐漸增大，其原因為：傳播時間的增大往往意味著傳播距離的增加，即傳播的波長數(shù)增大，介質(zhì)粘滯性的累積效果增加。

圖9 縱橫波速度模型

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

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

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

4 結(jié)論

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

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

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