各向異性地層深度域速度建模

李漢青,李向陽,侯 嵩

(1.中國石油大學(xué)(北京) CNPC物探重點實驗室，北京 102249；2.CGG公司英國技術(shù)中心,英國 雷德希爾 RH11PP)

偏移成像是常規(guī)地震數(shù)據(jù)處理的最后一步，也是十分關(guān)鍵的步驟。但在實際處理中地下構(gòu)造復(fù)雜且目標儲層深度較大，疊后時間偏移成像方法效果往往不佳。疊前深度偏移(Pre-stack Depth Migration，PSDM) 是解決速度橫向變化劇烈的復(fù)雜地區(qū)的地震資料成像的關(guān)鍵技術(shù)。它可以使散射波收斂，并將地震反射波移動到其在深度域中的正確位置。它可以提供比其他成像方法更好的復(fù)雜地下結(jié)構(gòu)圖像[1-2]。

當前對于縱波和轉(zhuǎn)換波的多分量偏移處理，是通過波場分離將數(shù)據(jù)分離為縱波和橫波場，然后基于聲學(xué)假設(shè)對分離后的波場進行縱波和橫波標量偏移成像[3-4]。這種處理方法沒有考慮彈性矢量波場的特性，在波場分離的過程中，不同波形的能量不能完全分離。剩余的非純波能會在偏移結(jié)果中產(chǎn)生大量的串擾噪聲，成像效果嚴重降低[5]。

基于矢量波場延拓的彈性波偏移成像技術(shù)為解決此類問題提供了有效方法。Pao和Varatharajulu[6]推導(dǎo)出了彈性波Kirchhoff-Helmholtz積分形式。秦福浩等[7]在此基礎(chǔ)上推出了彈性波Kirchhoff積分偏移公式。榮駿召等[8]從多分量地震波矢量場的特征出發(fā)，基于均勻各向同性彈性波波動方程，利用彈性波波動方程的矢量特性，在偏移過程中將泄漏的地震波信息還原到PP波與PS波分量中，實現(xiàn)了偏移過程中的波場分離。

與其他偏移方法相比，PSDM需要更精確的速度模型。傳統(tǒng)的速度模型更新方法可能不夠精確。疊前深度偏移本身對速度模型的敏感性高，是建立速度模型的有效工具。偏移速度分析有兩個主要原則：深度聚焦分析(Depth-Focusing Analysis，DFA)和剩余曲率分析(Residual-Curvature Analysis，RCA)。深度聚焦分析是在利用疊加功率測量速度誤差的基礎(chǔ)上進行的，剩余曲率分析是基于剩余時差來測量速度誤差的。在DFA中，如果偏移深度和聚焦深度的差為零，偏移速度是可以接受的。在RCA中，如果不同偏移量的成像深度差為零，偏移速度是可以接受的[9-10]。

再者，由于實際地下介質(zhì)具有彈性各向異性，因此在深度偏移過程中需要對介質(zhì)中的各向異性進行補償。Crampin等[11-14]的研究發(fā)現(xiàn)，地下介質(zhì)中存在廣泛的各向異性。巖石的層狀結(jié)構(gòu)、定向排列的裂隙都會引起各向異性。橫波對于各向異性的敏感，使其在多波地震資料的處理和解釋中有了重要的價值。轉(zhuǎn)換波速度建模比縱波速度建模難度更大。因為，對于P-S轉(zhuǎn)換波，下行波為P波，上行波為S波，這使得速度分析可靠性降低。在各向異性介質(zhì)中，波的傳播速度不僅是空間坐標的函數(shù)，而且與入線方向有關(guān)。如果不考慮各向異性的走時，共成像點(Common Image Point，CIP)道集在大偏移量下將不會被拉平，從而影響成像質(zhì)量。為此提出了一種估計地下各向異性信息并將其應(yīng)用于疊前深度偏移的方法。

針對這些不足，提出了基于矢量波波場延拓的彈性波疊前深度偏移速度分析和各向異性參數(shù)的估計方法，更新了速度結(jié)構(gòu)，建立了縱波和轉(zhuǎn)換波速度模型。該方法基于RCA的最基本原理，避免了數(shù)學(xué)逼近，可用于復(fù)雜速度結(jié)構(gòu)的建立。

1 彈性波Kirchhoff積分法

鑒于傳統(tǒng)多分量處理方法沒有考慮彈性矢量波場的特征，剩余能量會造成偏移結(jié)果中的噪聲干擾。本文采用矢量波波場延拓的彈性波偏移成像技術(shù)。彈性波Kirchhoff-Helmholtz積分方程，可以同時向下延拓P波和PS波的波場。Pao等[6]推導(dǎo)出了彈性波Kirchhoff-Helmholtz積分形式[7]。Kuo等[15]將此方法應(yīng)用于單界面模型數(shù)據(jù)的彈性波Kirchhoff疊前深度偏移，同時輸入X、Y、Z三分量的單炮數(shù)據(jù)，分別得到X、Y、Z三個方向的偏移剖面。在自由表面邊界條件下，地表應(yīng)力為零，可以得到簡化的偏移方程：

(1)

下標標注方式滿足愛因斯坦求和約定。

根據(jù)射線出射角將矢量位移分別投影到PP分量和PS分量：

(2)

其中，θ為位移水平分量與x軸夾角；φ為位移與軸夾角?？紤]二維多分量勘探，將uy略去，得到二維彈性波Kirchhoff積分法矢量偏移公式：

(3)

(4)

實際介質(zhì)是各向異性的，地下繞射點到地表接收點的PP波與PS波的射線路徑不再是直線，兩者的射線路徑、走時和出射角均不相同。因此需要對均勻各向同性介質(zhì)偏移公式進行改進。走時t，tP，tS，路徑出射角θP、θS和等效縱橫波速度VP、VS通過各向異性介質(zhì)的射線追蹤和速度分析得到。求得方程(3)、(4)寫為：

(5)

(6)

式中，tp為上行P波走時，ts為上行SV波走時；SP、SS分別為上行P、SV波的射線路徑長度；θP、θS分別為P、SV波射線出射角；VP、VS分別為P、SV波速度。

與基于聲波方程的Kirchhoff積分偏移算法比較，彈性波Kirchhoff積分矢量偏移算法具有很多優(yōu)點：具有波場分離的特性，偏移時同時分離縱波和轉(zhuǎn)換波；偏移算法具有保幅性；各向異性射線追蹤和速度分析使算法適用于各向異性介質(zhì)。

2 深度偏移速度更新方法

利用疊前深度偏移對速度模型的敏感性，將它作為一種有力的速度分析工具。在進行速度建模時，可以通過疊前時間偏移結(jié)果和Dix公式得到一個初始的速度模型。在這個模型基礎(chǔ)上給定一個初始速度范圍和步長，通過疊前深度偏移得到不同速度下的CIP道集。

剩余曲率分析是根據(jù)剩余旅行時差測量速度誤差，將剩余曲率分析原理應(yīng)用于共成像點(Common Imaging Point，CIP)道集，得到共成像點道集拉平準則，即當偏移速度準確時共成像點道集上不同炮檢距的成像深度相同；當偏移速度偏低時，同相軸上翹；反之，偏移速度偏高時，同相軸下垂。

對于各向異性地層，Sarkar和Tsvankin[16-18]給出各向異性近似下縱波成像深度與炮檢距關(guān)系

(7)

同理，上述將方程拓展到SV波，方程形式與縱波方程相同[2]，可以表示為

(8)

首先按照介質(zhì)為各向同性時進行速度分析，此時r2為0。本文采用的速度更新方法就是在此基礎(chǔ)上提出的。速度更新步驟如下：

(1)在一定的速度范圍內(nèi)，通過疊前深度偏移構(gòu)建CIP道集。

(2)通過參數(shù)r1拉平CIP道集，并拾取每個速度值所對應(yīng)的r1的值(圖1)。

圖1 不同偏移速度得到的CIP道集的r1參數(shù)譜Fig.1 r1parametric spectrum of CIP gathers obtained at different migration velocities

(3)通過線性回歸的方法計算r1=0時的速度值，用這個速度作為估計速度更新速度模型。

3 深度偏移各向異性參數(shù)反演

在各向同性理論下，地震波速度不隨地下介質(zhì)空間坐標變化；而在各向異性介質(zhì)中，地震波速度隨著射線角度的變化而改變。

Thomsen[13]用5個各向異性參數(shù)VP0、VS0、ε、δ、γ描述了速度的弱各向異性，VP0、VS0是垂直縱波和橫波速度，橫向各向同性介質(zhì)(Vertical Transverse Isotropy，VTI)中地震波速度與這5個參數(shù)的關(guān)系式：

VP(θ)=VP0(1+δsin2θcos2θ+εsin4θ)

(9)

(10)

VSH(θ)=VS0(1+γsin2θ)

(11)

上一節(jié)方程(7)給出了各向異性地層偏移速度、偏移深度和真實速度的關(guān)系。由方程可以看出CIP道集上的剩余時差是非雙曲函數(shù)，其大小由參數(shù)Vnm0和η決定。其中Vnm0決定近偏移距時差(二次)，而η對時差的影響在遠偏移距變得更大(4次)。如果偏移值和真實值一致，則偏移深度z(h)與半跑間距h無關(guān)，即同相軸變平。雖然方程(7)是假設(shè)反射層水平時推導(dǎo)而來，但Sarkar證明即使在地層傾角較大時也成立[18]。

為了定量分析方程(7)中各向異性參數(shù)ε、δ對偏移深度的影響，建立含有2個反射層的各向異性介質(zhì)速度模型，第1層是傾斜層，第2層為水平層。第1層和第2層速度以及各向異性參數(shù)相同。VP0=2 000 m/s，ε=0.3，δ=0.1，則η=ε-δ/1+2δ=0.167。運用CWPSU程序susynlvfti正演合成數(shù)據(jù)。再運用第一節(jié)所述基于矢量波波場延拓的彈性波偏移技術(shù)進行偏移得到CIP道集。由于各向異性參數(shù)的存在，偏移的關(guān)鍵步驟是各向異性射線追蹤。這里運用SU程序ray2dan進行各向異性射線追蹤計算走時和出射角，該程序依據(jù)旁軸近似將每條射線計算的旅行時外推到相鄰的網(wǎng)格點上。

圖2給出了在VP0=2 000 m/s正確時，各向異性參數(shù)ε、δ的錯誤取值和正確取值時的CIP道集。當選用正確的參數(shù)ε=0.3,δ=0.1進行偏移時，CIP道集同相軸是平直的(圖2(c))。圖2(a)是假設(shè)地層為各向同性時(即ε=δ=η=0)的偏移結(jié)果，由于忽略了各向異性參數(shù)的影響，同相軸上翹。圖2(b)是各向異性參數(shù)欠估計時的偏移結(jié)果(ε=0.3,δ=0.05，η=0.136)，與之對應(yīng)圖2(d)是各向異性參數(shù)過估計時的偏移結(jié)果(ε=0.4,δ=0.15，η=0.192)，同相軸下傾。

圖2 各項異性參數(shù)對CIP道集同相軸的影響Fig.2 Influence of anisotropy parameters on CIP gathers

這種速度各向異性在CIP道集上表現(xiàn)為遠偏移距的同相軸無法拉平。由方程(9)、(10)可知，P波和SV波在傳播過程中是耦合在一起的，同時由Thompson參數(shù)ε、δ決定。因為P波剖面有較強的能量和較高的信噪比，所以本文通過縱波反演Thompson參數(shù)ε、δ。

如上分析，當?shù)貙哟嬖诟飨虍愋詴r，即使VP0準確，在偏移獲得的CIP道集中，成像深度仍然隨炮檢距變化。由于各向異性的存在，不同偏移距應(yīng)當采用不同的偏移速度。因為成像點Zh和Z0的走時相同，有

(12)

這里Vθ即公式(9)中的VP(θ)，將方程(9)與(11)聯(lián)立，得

(13)

方程(13)可改寫為

Aδ+Bε=C

(14)

(15)

(16)

(17)

A、B、C的值與成像深度和炮檢距有關(guān)，成像深度和相對應(yīng)的炮檢距可以從CIP道集上拾取，如圖3、圖4所示。

圖3 成像深度為炮檢距的函數(shù)Fig.3 Relation between imaging depths and offsets

圖4 拾取成像深度與炮檢距Fig.4 Picking imaging depths and offsets

然后，根據(jù)公式(17)計算得到列向量A、B、C。最后通過求解廣義逆矩陣的方法，通過公式(16)求解Thompson參數(shù)ε和δ。

同理，對于Thompson參數(shù)γ，可以在純橫波CIP道集上進行各向異性參數(shù)分析。

(18)

γ=A-1·B

(19)

計算得到列向量A、B(圖4)。再通過求解廣義逆矩陣的方法，通過式(19)求解Thompson參數(shù)。

4 實際工區(qū)模型測試

以勝利油田羅家工區(qū)Ken71模型進行試驗。Ken71模型為復(fù)雜的地質(zhì)構(gòu)造模型。該地區(qū)具有斷裂系統(tǒng)復(fù)雜、微幅構(gòu)造難以識別、薄互層發(fā)育等地質(zhì)特征，是典型的河流相與三角洲相沉積類型油藏[19]。針對薄互層地質(zhì)特點，建立了典型的Ken71地質(zhì)模型并進行了正演模擬。各向異性速度建模和各向異性參數(shù)分析流程如圖5所示，P波和S波速度場模型如圖6所示。

圖5 各向異性速度建模和各向異性參數(shù)分析流程Fig.5 Anisotropic velocity modeling and anisotropic parameter analysis process

圖6 Ken71速度場模型Fig.6 Velocity model of Ken71

圖中紅色括號標示出目標薄互層位置，1 480～1 630 m，其各向異性參數(shù)為ε=0.20,δ=0.05，γ=0.10。網(wǎng)格大小為2 766×901，橫向的采樣間隔為10 m，而縱向的采樣間隔為2 m。正演彈性波記錄共有75炮，每炮202道接收，道間距為20 m，此模型時間采樣點共有3 000個，采樣點間時間的間隔為1 ms[20],從中可以看出地質(zhì)構(gòu)造的復(fù)雜性。

正演彈性波場的X分量和Z分量單炮記錄如圖7所示，可以看出由于地下地質(zhì)構(gòu)造的復(fù)雜性導(dǎo)致地震記錄信息也非常復(fù)雜。

圖7 Ken71模型正演記錄Fig.7 Synthetic shot gathers of Ken71

運用圖5提出的速度建模和各向異性參數(shù)分析流程，得到縱波和轉(zhuǎn)換波速度模型，如圖8所示。在速度更新時為簡化層剝過程，只拾取了7個反射界面，逐層計算了縱波和轉(zhuǎn)換波的速度。

圖8 速度更新得到的速度模型Fig.8 Velocity model after velocity analysis

在假設(shè)速度模型準確的基礎(chǔ)上對第六層深度為1 400～1 680 m處當作各向異性地層(VTI)處理，反演得到各向異性參數(shù)平均結(jié)果約為ε=0.22,δ=0.05，γ=0.12。深度偏移得到X=3000 m處縱波CIP道集如圖9所示，其中圖9(a)為忽略各向異性(即各向同性ε=δ=γ=0，偏移得到CIP道集，圖9(b)為利用各向異性參數(shù)ε=0.22,δ=0.05、偏移得到CIP道集。對比二圖看到各向異性建模偏移的CIP道集拉平效果較好。

圖9 深度偏移得到的縱波CIP道集Fig.9 P-wave CIP gathers from depth migration

類似圖9，利用各向同性和各向異性速度建模最終偏移的縱波成像結(jié)果和轉(zhuǎn)換波成像結(jié)果如圖10和圖11所示。

圖10 深度偏移得到的疊加剖面Fig.10 Stacked profile from depth migration

圖11 轉(zhuǎn)換波深度偏移得到的疊加剖面Fig.11 Stacked profiles obtained by depth migration of converted waves

深度偏移得到的轉(zhuǎn)換波CIP道集如圖12所示，其中圖12(a)為各向同性偏移得到的CIP道集，圖12(b)為利用各向異性參數(shù)γ=0.12偏移得到的CIP道集，在VTI層(上部紅色箭頭處)，二者拉平效果接近，但對深層區(qū)別較大(下部紅色箭頭處)?？v橫波成像深度基本一致。各向異性參數(shù)對實際資料拉平效果沒有前述模擬資料明顯，主要原因可能是實際資料各向異性參數(shù)較小，另外實際資料的信噪比較低和Ken71地層復(fù)雜的結(jié)構(gòu)有關(guān)。但總體拉平效果有所改善，估計的參數(shù)值可以接受。

圖12 深度偏移得到的轉(zhuǎn)換波CIP道集Fig.12 Converted wave CIP gathers obtained by depth migration

圖12中圓圈處顯示偏移結(jié)果得到了提高，對于縱波，同相軸連續(xù)性變得更好，對于轉(zhuǎn)換波由于CIP道集拉平，疊加結(jié)果使疊加剖面同相軸更清晰。圖13是各向異性速度建模最終偏移的縱波成像結(jié)果和轉(zhuǎn)換波成像結(jié)果對比。

圖13 Ken71地震數(shù)據(jù)縱波和轉(zhuǎn)換波深度偏移剖面對比Fig.13 Comparison of P-wave and converted-wave depth migration profiles of Ken71 seismic data

從圖13中可以看出，兩者對地下復(fù)雜的地質(zhì)構(gòu)造例如起伏地表、斷層等都具有不錯的成像效果，結(jié)合地震解釋資料，深度剖面上對含油透鏡體、油水薄互層、油頂油水同層以及油氣水同層等主要層位都進行了準確的刻畫，轉(zhuǎn)換波的分辨率高于縱波?？梢钥吹皆谏疃? 400～1 680 m處正好是油水薄互層，薄互層一般等效為VTI地層，與前文設(shè)定的各向異性層對應(yīng)，層位偏移歸位到真實位置。

5 結(jié)論

(1)本文利用Kirchhoff積分疊前深度偏移公式，考慮各向異性介質(zhì)的射線追蹤，實現(xiàn)了各向異性地層矢量偏移，輸入多分量地震數(shù)據(jù)，同時分離得到縱波和轉(zhuǎn)換波的成像剖面。

(2)之后分析了縱波和轉(zhuǎn)換波深度域地層速度模型建模方法。本方法基于RCA原理，在CIP道集上進行速度更新，并將偏移速度與真實速度的差異參數(shù)化，通過拾取能量團的方式估計真實速度。然后對目標層進行疊前深度偏移，拾取速度結(jié)構(gòu)，實現(xiàn)速度建模。在速度模型準確的基礎(chǔ)上，本文提出了新的各向異性參數(shù)分析方法，推導(dǎo)出了各向異性參數(shù)、偏移距與成像深度之間的關(guān)系式，通過求解廣義逆矩陣的方法反演各向異性參數(shù)，實際數(shù)據(jù)測試取得了較好的效果，驗證了方法的可行性。