鄧馨卉，劉 財，郭智奇，劉喜武，劉宇巍

1.吉林大學地球探測科學與技術學院，長春 130026 2.頁巖油氣富集機理與有效開發(fā)國家重點實驗室，北京 100083 3.中國石化頁巖油氣勘探開發(fā)重點實驗室，北京 100083 4.中國石化石油勘探開發(fā)研究院，北京 100083

0 引言

近年來，隨著人們對油氣資源需求量的增加，以及常規(guī)油氣資源的不斷減少，頁巖油氣的研究與勘探逐漸成為國內油氣資源研究的熱門領域。位于渤海灣盆地的濟陽坳陷為富油氣坳陷，其中沾化凹陷羅家地區(qū)泥頁巖厚度大(1.0～1.5 m)，裂縫發(fā)育良好，沙河街組三段(簡稱沙三段)中下部的油泥頁巖為裂縫油氣藏主要賦存區(qū)域[1-10]。在地震分析的過程中，將井震數據結合對儲層特性進行研究十分重要。但由于測井資料與地震資料分辨尺度差異較大，將測井數據應用到地震響應計算中分析的方法并未得到充分的研究和應用。

為了將測井數據粗化到地震數據尺度，可采用Backus平均理論[11]對測井數據進行粗化，在長波長條件下將層狀介質等效為VTI(vertical transversely isotropy)介質。Schoenberg等[12]將該理論推廣到了一般的層狀各向異性介質中；Kumar[13]將Backus平均理論應用到具有低對稱各向異性的裂縫型頁巖儲層中。

地層中裂縫的發(fā)育引起地震波速度和反射振幅隨方位角的變化，因此可根據地震響應的方位各向異性特征識別和描述裂縫屬性。傳統(tǒng)的AVO(amplitude versus offset)研究以Zoeppritz方程[14]為基礎。由于Zoeppritz方程的復雜性，很多學者對其進行簡化并應用于實際研究當中。Banik[15]和Thomsen[16]提出弱各向異性參數可以反映TI(transeversely isotropy)介質中反射界面兩側時差速度和反射振幅的變化。Rüger[17]基于Banik和Thomsen的研究，給出了方位各向異性介質在任意方位角、大入射角情況下的縱波反射系數近似公式。AVAZ(amplitude versus azimuth)分析研究地震子波經過界面反射后振幅隨入射角和方位角的變化，通常是由一組垂直裂縫引起的各向異性。Sayers等[18]的研究表明，在地震勘探中，垂直裂縫的發(fā)育會引起更強更明顯的AVAZ異常響應。Schoenberg等[19]于1999年應用傳播矩陣算法[20]分析了裂縫儲層中地震反射波與方位角相關的調諧效應。

地震正演模擬可以有效地觀測地震波在各向異性介質中傳播的規(guī)律和特點，主要方法有反射率法、有限差分法和射線追蹤法，其中反射率法的精度高于其他兩種方法。Fuchs等[21]最初提出將反射率法應用于各向同性層狀介質地震響應的計算。Sheriff[22]對反射率法在各向同性介質中的應用進行了更詳盡的論述。Booth等[23]將反射率法推廣至各向異性層狀介質的計算。Fryer等[24]給出了存在彈性水平對稱面情況下的任意各向異性介質所對應的系統(tǒng)矩陣形式。Mallick等[25]通過反射率法計算得到方位各向異性層狀介質的合成地震記錄。

本文基于各向異性條件下的反射率法，結合Schoenberg傳播矩陣法中水平慢度的表達方式和VTI各向異性介質中的相速度公式，推導出了針對上覆層為各向異性介質的反射率法，并模擬得到全波場地震AVAZ響應。然后結合羅家地區(qū)裂縫型頁巖油儲層測井資料，應用各向異性Backus平均理論進行井震匹配，并通過Schoenberg線性滑動理論[26]引入垂直裂縫建立巖石物理模型，將儲層裂縫特征與地震響應較好地聯系起來；再應用各向異性反射率法進行各向異性地震AVAZ響應模擬并獲得方位振幅分布，由此來分析羅家地區(qū)頁巖油儲層的方位各向異性地震響應特征。

1 羅家地區(qū)各向異性頁巖油儲層測井數據分析

羅家鼻狀構造帶北鄰渤南洼陷，南靠陳家莊凸起，研究井L位于濟陽坳陷沾化凹陷羅家鼻狀構造帶西翼。研究區(qū)域2 909.0～3 132.0 m為沙河街組三段中下部地層，發(fā)育厚層油泥巖、油頁巖，地層能量充足，鉆井過程中油氣顯示活躍。L井中礦物成分主要為方解石、黏土礦物和石英，同時含有少量白云石、黃鐵礦、斜長石和鉀長石等。圖1a為L井井區(qū)沙三段中下部地層礦物類型及體積分數示意圖。可以看出，隨著地層由老到新的變化，黏土礦物和石英的平均體積分數呈現逐漸增加的趨勢，方解石和白云石的平均體積分數則呈現逐漸降低的趨勢。

圖1b—f為L井頁巖儲層測井數據，在縱波速度、密度以及孔隙度曲線中存在明顯的突變處，即深度3 060.0 m處：此界面以上地層的縱波速度與密度明顯較低，孔隙度較高。這是由于富有機質頁巖儲層在常規(guī)測井曲線中具有一定的規(guī)律性，即高電阻率、高聲波時差、低密度等特征，與其他非泥頁巖層段之間具有較大區(qū)別，更易于進行識別。

圖2為L井沙河街組三段下亞段的全井眼地層微電阻掃描成像(FMI, formation micro-resistivity scanning imaging)[27]測井結果。通過放大的FMI圖像可以看出，該段地層內發(fā)育角度較低的層間縫和高角度的高導縫，且水平層理發(fā)育。地層的VTI各向異性主要由水平層理和層間縫引起；而高角度的高導縫為垂直裂縫，在VTI各向異性的背景下，垂直裂縫的存在使得頁巖儲層具有等效正交各向異性介質的特征，響應表現為方位各向異性特征。

圖1 羅家地區(qū)L井頁巖儲層測井數據Fig.1 Well log data from shale oil reservoir of Well L in Luojia area

圖2 羅家地區(qū)L井水平縫和垂直縫FMI成像結果Fig.2 FMI images of horizontal fractures and vertical fractures of Well L in Luojia area

2 基本原理

2.1 Backus平均理論

Backus平均理論常用于測井數據中估算長波長條件下薄層的彈性性質。在長波長的條件下，對VTI各向異性介質應用Backus平均理論，得到的等效介質可以由復數模量描述為

(1)

式中：Cij和Cije分別為Backus平均前、后等效介質的彈性剛度系數；〈·〉表示加權平均。

2.2 Schoenberg線性滑動理論

Schoenberg等[26]提出當地層中同時存在水平層理和垂直裂縫時，在長波長條件下可等效為正交各向異性介質。該正交各向異性模型由兩部分組成，第一部分為具有一個垂向對稱軸的橫向各向同性背景介質，第二部分是一組平行垂直裂縫。正交各向異性介質的剛度矩陣可以表示為

(2)

其中：

(3)

其中，無量綱量可定義為

(4)

式中：Zn為法向柔度；Zv為垂向正切柔度；Zh為水平正切柔度。假設δv=δh，并利用裂縫密度e來估計δv，δh，δn：

(5)

式中：g=(vS/vP)2，vS和vP分別為背景介質的S波和P波速度；孔隙縱橫比α=3φ/4πe，φ是裂縫的體積分數；μ是背景介質的剪切模量；k′和μ′分別為裂縫中流體成分的體積模量和剪切模量。

2.3 各向異性反射率方法

反射率方法可以模擬彈性或非彈性介質中地震波的反射透射問題，具有數值穩(wěn)定性高、精度高、計算成本低的特點。

首先假設介質僅在垂向變化，且不存在體力。為了對運動方程變量進行分離，將動量方程和本構方程進行三維傅里葉變換：

G(px,py,ω)=

(6)

式中：ω為角頻率；px、py為水平方向的慢度分量。G(x，y，t)表示笛卡爾坐標系下的動量方程和本構方程；t為時間。經過傅里葉變換后可得到微分方程系統(tǒng)[28]：

?zb(z)=iωA(z)b(z)。

(7)

式中：A為6×6的系統(tǒng)矩陣；b為應力位移矢量，表示連續(xù)穿過水平面的變量，

u=(ux,uy,uz)T，

(8)

式中：ux，uy，uz是位移的分量；τxz，τyz，τzz是應力張量的z方向分量。由此而得到6×6的系統(tǒng)矩陣A可由水平慢度px和py、密度ρ以及彈性參數cij(i,j=1, 2, 3, …, 6)來表示。如果彈性參數與角頻率無關，那么A與頻率無關。

當目標模型的上半空間為各向異性時，水平慢度px和py可表示為

(9)

式中：θ為入射角；Ω為方位角；vqP為上部介質的縱波相速度。

在求解微分方程系統(tǒng)的過程中，可通過引入經典的傳播算子法[28-29]來完成計算。傳播算子P(z,z1)為微分方程系統(tǒng)的唯一連續(xù)解：

?zP(z,z1)=iωA(z)P(z,z1)，

(10)

P(z,z1)=I。

(11)

式中：I是6×6的單位矩陣。當在深度z1處選擇b(z1)作為明確的邊界條件時，那么深度為z處的響應為b(z)=P(z,z1)b(z1)。通過系統(tǒng)矩陣A的特征值矩陣和特征向量矩陣，可求得當能量自上而下入射時，各向異性層狀介質的底界面反射、透射系數矩陣R：

(12)

式中：TU、TD分別為3×3上行波和下行波透射系數矩陣；RU、RD分別為3×3上行波和下行波反射系數矩陣。

通過Fryer等[30]引入的散射算子和散射矩陣計算并整理，可以得到層狀介質中第i+1層RDi+1、RUi+1、TDi+1、TUi+1和第i層RDi、RUi、TD、TUi之間的遞歸關系式[31]：

(13)

(14)

式中：rD、rU和tD、tU分別為界面深度z=zi+1處的反射矩陣和透射矩陣；h=zi+1-zi為深度差；λ為系統(tǒng)矩陣A的特征值。

計算反射響應是求取合成地震記錄的第一步，由計算得到RD中的頻率域反射響應rppD，與地震子波頻譜F(ω)相乘后進行逆傅里葉變換，可得到時間域合成地震記錄RppD：

(15)

3 數據試算

3.1 理論模型

根據Schoenberg等[19]的正交各向異性理論，設計一個三層理論模型進行分析，模型參數如表1所示，表中ε、γ、δ為各向異性參數。理論模型為三層正交各向異性介質，中間層復數剛度矩陣如圖3a所示。計算中用到的地震子波為40 Hz雷克子波。

理論模型中的正交各向異性介質由兩部分組成：VTI各向異性背景介質以及垂直裂縫系統(tǒng)，如圖3b所示。研究方位各向異性特征時，規(guī)定觀測面垂直于模型中的垂直裂縫發(fā)育面；當模型中僅具有VTI各向異性時，地震響應表現為方位各向同性，引入垂直裂縫后地震響應則具有方位各向異性。正交各向異性介質的地震響應不僅與背景介質的VTI各向異性參數有關，也與垂直裂縫密度有關，因此本文在計算中固定裂縫密度值為0.1。

依據各向異性反射率理論，大入射角處方位各向異性特征明顯，但臨界角處會發(fā)生全反射，因此選取入射角30°的地震響應進行分析。圖4為入射角為30°，反射PP波、PSV波、PSH波在方位角0°～360°時的地震AVAZ響應。3種反射波的AVAZ響應均表現出明顯的周期性和對稱性，且均可以觀察到明顯的地震方位各向異性特征，即介質中垂直裂縫的存在導致不同方位角處反射振幅產生明顯的變化。提取目標層振幅并繪制極坐標圖(圖5)，可直觀地觀察到反射振幅的特征；對三者的方位振幅分布擬合成橢圓，其中PP波和PSV波的長軸方向與裂縫方向一致，PSH波的長軸方向與裂縫方向垂直。

表1 理論模型參數

COrth=30.929.8612.620009.8643.3616.4700012.6216.4733.7900000012.100000009.3300000011.76aba.矩陣描述;b. 模型描述。圖3 正交各向異性介質示意圖Fig.3 Schematic of orthorhombic medium

圖4 入射角為30°時PP波(a)、PSV波(b)、PSH波(c)地震AVAZ響應Fig.4 Seismic AVAZ responses of PP wave(a), PSV wave(b) and PSH wave(c) at incident angle of 30°

圖5 入射角為30°時PP波(a)、PSV波(b)、PSH波(c)方位振幅分布Fig.5 Azimuth versus amplitude in polar coordinate of PP wave(a), PSV wave(b) and PSH wave(c) at an incident angle of 30°

將各向異性反射率法應用于理論模型，得到的計算結果符合其構造和彈性特征，證明了該方法的有效性以及對于正交各向異性介質的適用性，同時也能夠較好地觀察到相應的裂縫特征。

3.2 羅家地區(qū)L井測井數據應用

3.2.1 儲層VTI各向異性

根據地層特性給定初始各向異性參數，應用各向異性Backus平均理論粗化測井數據，得到等效VTI各向異性介質，并求取尺度粗化后的各向異性參數。

圖6為L井測井曲線，原始測井曲線采樣間隔為0.1 m，應用各向異性Backus平均理論計算后曲線的計算長度為1 m[32]。巖石中存在兩種各向異性：與層狀結構有關的各向異性以及巖石的固有各向異性。當地層中僅存在與層狀結構有關的各向異性時，固有各向異性參數為ε=0、γ=0、δ=0。目前對于羅家地區(qū)的研究中沒有針對巖心的具體描述，因此當地層中同時具有兩種各向異性時，針對L井：自然伽馬數值大于55 API指示巖性為VTI各向異性的頁巖，可引入Sondergeld等[33]研究中給出的頁巖固有各向異性參數ε=0.256、γ=0.481、δ=0.051；自然伽馬數值小于等于55 API時，固有各向異性參數均為0。圖6f、g分別為兩種假設條件下計算得到的VTI等效介質的地震各向異性參數。

各向異性Backus平均理論在兩種假設條件下計算得到的等效介質彈性參數曲線沒有明顯變化，這與Backus提出的理論相符，即各向同性介質與VTI各向異性介質平均后得到的等效介質均為VTI各向異性，同時驗證了本文應用Backus平均理論計算的正確性。對比兩種假設的各向異性參數可以看出：ε和δ與縱波速度成相反的變化趨勢；而且地層中同時存在兩種各向異性時，各向異性參數數值變化更明顯，對各采樣點特征的反映更敏感，有助于精確描述地下結構特征。

3.2.2 儲層正交各向異性

利用Schoenberg線性滑動理論將垂直裂縫引入VTI各向異性等效介質，構建正交各向異性介質模型，計算中e=0.1。

圖7a正交各向異性介質模型中沿x、y、z方向的縱波速度分別為v11、v22和v33。v33

各向異性參數初始值：f. ε=0, γ=0, δ=0；g. ε=0.256, γ=0.481, δ=0.051。黑色曲線為原始測井曲線，灰色曲線為Backus平均后曲線。圖6 羅家地區(qū)L井及等效介質測井曲線Fig.6 Well log data from Well L in Luojia area and equivalent medium

圖7 羅家地區(qū)L井方位各向異性模型中各方向速度Fig.7 Velocities in different directions of the orthorhombic model of Well L in Luojia area

由以上分析可知，在羅家地區(qū)L井測井數據構建的正交各向異性介質模型中存在沿z方向的垂直裂縫，且VTI各向異性介質中發(fā)育水平層理和沿y方向的水平裂縫，計算得到的縱橫波各方向速度大小符合正交各向異性特征。

3.2.3 全波場地震響應模擬及分析

對理論模型的計算驗證了應用反射率法模擬正交各向異性介質地震響應的可行性。當入射角過大時會發(fā)生相位反轉；為了保證計算效果和分析的準確性，在接下來的地震AVAZ響應計算中選取入射角范圍為0°～40°。

前文對羅家地區(qū)L井建立巖石物理模型，并計算各方向上的縱、橫波速度，對儲層中的構造進行了初步分析?，F綜合測井、巖石物理和地震資料進行層位標定：首先對測井數據及各方向速度進行深-時轉換，然后應用各向異性反射率法模擬地震AVAZ響應，并將合成地震記錄上的標志層與測井曲線以及速度曲線進行標定。由于本文的聲波和密度測井數據不是全井段數據，模擬的合成地震記錄在地震剖面上存在浮動現象，因此需要通過比對測井數據和合成地震記錄中的標志性層位，來建立正確的深-時對應關系。

圖8為對L井進行的井震標定[34-36]示意圖?？梢钥闯觯簯梅瓷渎史ㄓ嬎愕玫降腁VAZ響應顯示出的振幅變化與速度突變處一致；證明了對聲波數據建立的深-時轉換曲線的準確性，也證明了各向異性反射率法在研究區(qū)的適用性。因此，在測井曲線上表現不明顯的信息也可以通過地震資料來觀測。

圖9為入射角為30°的反射PP波、PSV波和PSH波的全波場地震AVAZ響應，顯示方位角范圍為0°～360°?？梢钥闯?，固定入射角度時，方位各向異性特征隨著方位角的變化有著明顯的規(guī)律性變化：在PP波和PSV波的AVAZ響應中這種特征存在周期性，180°為一個周期，在一個周期內振幅隨方位角的變化以90°為軸兩側具有對稱性；PSH波的AVAZ響應則是以180°為軸兩側對稱，而且可以觀察到振幅隨方位角的變化發(fā)生了極性反轉。另外，在AVAZ響應中可以明確觀察到：PP波和PSV波觀測方向與垂直裂縫發(fā)育方向垂直時(方位角為0°，180°，360°)，地震響應中的振幅值應較小，當觀測方向與垂直裂縫發(fā)育方向平行時(方位角為90°，270°)振幅達到峰值；PSH波與之相反。同時，深層位置的方位各向異性特征比淺層更明顯。

將速度突變界面處的3種反射波振幅提取出來，觀測其隨方位角的變化規(guī)律，結果如圖10所示。PP波的反射振幅隨著方位角變化的分布可擬合成橢圓形，長軸方向與模型中的垂直裂縫發(fā)育方向一致；PSV波的反射振幅隨方位角變化成數字“8”分布，垂向對稱軸方向與垂直裂縫方向一致；PSH波的反射振幅成“∞”形狀分布，水平向對稱軸與垂直裂縫方向垂直。在方位振幅分布上，可以看出3種反射波的圖像在垂直和水平兩個方向上均存在對稱性，且反射PSV波和PSH波均存在零值點；但是結合全波場地震AVAZ響應可以發(fā)現，反射PSH波的反射振幅產生了極性反轉。因此，結合地震AVAZ響應和反射波方位振幅分布，可以清晰地觀察到方位各向異性特征，并以此判斷頁巖儲層中垂直裂縫的走向。

a.縱波速度；b. 反射PP波全波場AVAZ響應。圖8 羅家地區(qū)L井的井震標定結果Fig.8 Well-to-seismic calibration of Well L in Luojia area

入射角為30°。圖9 羅家地區(qū)L井反射PP波(a)、PSV波(b)、PSH波(c)地震AVAZ響應Fig.9 Seismic AVAZ responses of PP wav(a)e, PSV wave(b) and PSH wave(c) of Well L in Luojia area

圖10 羅家地區(qū)L井反射PP波(a)、PSV波(b)、PSH波(c)方位振幅分布Fig.10 Azimuth versus amplitude in polar coordinate of PP wave(a), PSV wave(b), PSH wave(c) of Well L in Luojia area

4 結論

針對羅家地區(qū)頁巖儲層構造特征，本文提出并實現了基于測井數據建立地震尺度非均勻的正交各向異性速度模型，并利用各向異性反射率法對儲層反射PP波、PSV波和PSH波的地震AVAZ響應進行模擬，主要結論如下：

1)通過在地震尺度非均勻各向異性條件下構建正交各向異性等效介質模型，模擬了水平層理發(fā)育垂直裂縫的儲層地質特征和彈性特性。

2)針對羅家地區(qū)L井建立正交各向異性模型，并應用各向異性反射率法，計算其反射PP波、PSV波和PSH波地震AVAZ響應，清晰展示了由于裂縫的發(fā)育，地震響應的方位各向異性特征及其變化規(guī)律。相比于淺層位置，深層的方位各向異性特征更明顯。

3)提取標志層反射PP波、PSV波和PSH波反射振幅繪制全波場方位振幅分布圖，同時結合地震AVAZ響應，可以觀察到裂縫相對于觀測方向的發(fā)育角度，PP波和PSV波的分布圖中，擬合成橢圓的長軸方向與裂縫方向一致，PSH波的分布圖擬合成橢圓的長軸方向則與裂縫方向相垂直。因此，全波場地震響應可作為判斷裂縫方向的指示工具，為頁巖儲層裂縫識別和預測提供依據。