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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

2 基本原理

2.1 Backus平均理論

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

(1)

式中：Cij和Cije分別為Backus平均前、后等效介質(zhì)的彈性剛度系數(shù)；〈·〉表示加權(quán)平均。

2.2 Schoenberg線性滑動理論

Schoenberg等[26]提出當?shù)貙又型瑫r存在水平層理和垂直裂縫時，在長波長條件下可等效為正交各向異性介質(zhì)。該正交各向異性模型由兩部分組成，第一部分為具有一個垂向?qū)ΨQ軸的橫向各向同性背景介質(zhì)，第二部分是一組平行垂直裂縫。正交各向異性介質(zhì)的剛度矩陣可以表示為

(2)

其中：

(3)

其中，無量綱量可定義為

(4)

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

(5)

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

2.3 各向異性反射率方法

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

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

G(px,py,ω)=

(6)

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

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

(7)

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

u=(ux,uy,uz)T，

(8)

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

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

(9)

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

在求解微分方程系統(tǒng)的過程中，可通過引入經(jī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處的響應(yīng)為b(z)=P(z,z1)b(z1)。通過系統(tǒng)矩陣A的特征值矩陣和特征向量矩陣，可求得當能量自上而下入射時，各向異性層狀介質(zhì)的底界面反射、透射系數(shù)矩陣R：

(12)

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

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

(13)

(14)

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

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

(15)

3 數(shù)據(jù)試算

3.1 理論模型

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

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

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

表1 理論模型參數(shù)

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

圖4 入射角為30°時PP波(a)、PSV波(b)、PSH波(c)地震AVAZ響應(yīng)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°

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

3.2 羅家地區(qū)L井測井數(shù)據(jù)應(yīng)用

3.2.1 儲層VTI各向異性

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

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

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

3.2.2 儲層正交各向異性

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

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

各向異性參數(shù)初始值：f. ε=0, γ=0, δ=0；g. ε=0.256, γ=0.481, δ=0.051。黑色曲線為原始測井曲線，灰色曲線為Backus平均后曲線。圖6 羅家地區(qū)L井及等效介質(zhì)測井曲線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井測井數(shù)據(jù)構(gòu)建的正交各向異性介質(zhì)模型中存在沿z方向的垂直裂縫，且VTI各向異性介質(zhì)中發(fā)育水平層理和沿y方向的水平裂縫，計算得到的縱橫波各方向速度大小符合正交各向異性特征。

3.2.3 全波場地震響應(yīng)模擬及分析

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

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

圖8為對L井進行的井震標定[34-36]示意圖。可以看出：應(yīng)用反射率法計算得到的AVAZ響應(yīng)顯示出的振幅變化與速度突變處一致；證明了對聲波數(shù)據(jù)建立的深-時轉(zhuǎn)換曲線的準確性，也證明了各向異性反射率法在研究區(qū)的適用性。因此，在測井曲線上表現(xiàn)不明顯的信息也可以通過地震資料來觀測。

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

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

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

入射角為30°。圖9 羅家地區(qū)L井反射PP波(a)、PSV波(b)、PSH波(c)地震AVAZ響應(yīng)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 結(jié)論

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

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

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

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