馬 妮, 印興耀, 宗兆云, 孫成禹

(1.中國石油大學(華東) 地球科學與技術學院，青島 266580； 2.海洋國家實驗室 海洋礦產(chǎn)資源評價與探測技術功能實驗室，青島 266071)

0 引言

由于全球對能源的需求不斷增加且常規(guī)油氣資源的不斷耗竭，頁巖氣藏的勘探和開發(fā)逐漸成為專家學者關注的熱點，而開發(fā)頁巖氣藏的有效手段是對其進行水力壓裂。通過水力壓裂施工可以使儲層形成大量的裂縫，有利于油氣的儲存和運移。研究表明，裂縫密度與許多地層參數(shù)有關，如曲率、楊氏模量、地應力和各向異性參數(shù)等。通過對地層參數(shù)的研究，可以有效地預測儲層的裂縫密度，指導頁巖氣的勘探與開發(fā)。

曲線上某點的曲率定義為曲線方向的改變速度[1]。曲率屬性可以用于表示地層的彎曲程度，Murray等[2]得出曲率與裂縫密度存在一定關系的結論，曲率越大表明該區(qū)域的應變越大；Price等[3]推導了應力與曲率和彈性力學參數(shù)的定量關系，從而奠定了用曲率估算應力的理論基礎；Sheorey[4]提出用靜彈性熱應力模型估算地應力，該模型考慮了地層曲率和巖層彈性參數(shù)、密度與熱膨脹系數(shù)的變化；李志勇等[5]利用微分幾何學原理推導了主曲率的計算公式，并將趨勢面分析與差分法結合計算構造面主曲率預測裂縫密度；唐湘蓉[6]從古構造應力場出發(fā)，利用連續(xù)介質(zhì)有限元數(shù)值模擬法和巖石破裂準則預測古潛山儲層裂縫密度和裂縫方位；何英[7]推導了基于曲率分析的構造應力場計算公式，利用曲率估算地層的構造應力，克服了構建復雜構造應力場模型的難點；Hunt等[8]分析了與裂縫密度有關的地層參數(shù)，通過胡克定律推導應力與曲率和楊氏模量的關系，將曲率進行擴展得到裂縫預測參數(shù)預測裂縫密度。Starr[9]利用曲率和應變率的數(shù)學定義，得到了求取應變的簡化公式。

筆者根據(jù)Hunt[8]提出的預測裂縫密度的擴展曲率方法(即地應力近似等價于曲率與楊氏模量的乘積)進行了實際應用。首先利用基于tau-p變換的曲率屬性提取方法[10-11]提取地層的曲率屬性，再通過疊前地震反演得到地層的楊氏模量，將兩者進行結合得到應力-曲率屬性。曲率越大，該區(qū)域的應變越大，表明受到的應力越大。當給定曲率時，應力隨著楊氏模量的增加而增加，因此應力-曲率屬性高值表示該區(qū)域所受應力較大，增加了產(chǎn)生裂縫的可能性，有利于指導頁巖儲層的鉆井布置和水力壓裂開發(fā)等方面。

1 方法原理

1.1 基于局部tau-p變換的曲率屬性提取方法

曲率屬性是描述地層彎曲程度的一種地震幾何屬性，它能夠刻畫地層的撓曲、褶皺、斷裂以及河道等地質(zhì)構造特征，是間接預測裂縫密度的一種屬性[12-13]。王欣[10-11]提出了基于局部τ-p變換的曲率屬性提取方法，該方法是將(t,x)域的地震數(shù)據(jù)及其復地震數(shù)據(jù)根據(jù)一定的線性疊加關系變換至(τ,p)域，應用p和τ來描述波的軌跡。由于地震波場為點源波場，且遠源波場的局部小地震道集可以近似看作是局部平面波，因此對這種局部小道集進行τ-p變換稱為局部τ-p變換。利用局部τ-p變換可以將地震記錄上以時間延遲τ0為截距，以p0為斜率的同相軸(局部平面波)轉變成τ-p域中的一個點，其物理意義為遠源球面波場的局部平面波分解。對于一個三維數(shù)據(jù)體，其τ-p變換可表示為：

(1)

其中：m(τ,p,q)為三維τ-p域記錄；d(τ-pxj-qyj,xj,yj)為三維地震信號；J為分析窗口內(nèi)的地震道數(shù)；xj和yj分別是第j道相對分析窗口中心道之間的距離；p和q為不同方向的射線參數(shù)，表示地層的視傾角屬性。

將得到的三維τ-p域數(shù)據(jù)體帶入相似性掃描算法中，計算公式為式(2)。

(2)

通過相似性掃描可以得到視傾角數(shù)據(jù)體p和q以及地層的相干屬性。相對于常規(guī)中心窗算法而言，基于局部tau-p變換中心窗算法具有較高的穩(wěn)定性和識別精度，但其計算的相干屬性空間分辨率較低，且得到的傾角屬性體包含一定的干擾噪音。而采用基于局部tau-p變換中心窗-離心窗算法則解決了上述存在的問題，不僅使計算的相干屬性具有較高的空間分辨率，還降低了傾角屬性存在的干擾噪音[10]。最后，將得到的傾角數(shù)據(jù)體用于體曲率屬性計算中。Al-Dossary等[14]提出了分波數(shù)曲率提取法，該方法可以將傾角數(shù)據(jù)體通過傅氏變換從空間域轉換到波數(shù)域，將其乘以-ik因子相當于其在空間域的一階導數(shù)的傅氏變換，并引入了分波數(shù)導數(shù)的概念[12]，可以得到多尺度曲率屬性的計算公式：

(3)

(4)

最大正曲率表示地層凸起構造，因此選用最大正曲率作為地應力預測的基礎。

1.2 基于疊前地震反演的楊氏模量提取方法

宗兆云等[15-18]在Aki等[19]近似方程的基礎上，建立了縱波反射系數(shù)與楊氏模量反射系數(shù)、泊松比反射系數(shù)及密度反射系數(shù)的線性關系，稱之為YPD近似方程。

(6)

將公式(6)代入公式(5)可得式(7)。

(7)

(8)

其中:E0、σ0和ρ0為E、σ和ρ的平均值，為了求解方便，將其線性化，對公式(8)兩邊取對數(shù)，可得式(9)。

(9)

從公式(9)可以看出，提取楊氏模量、泊松比以及密度需要三個不同入射角的彈性阻抗體，將三個不同入射角度值代入公式(9)可得到方程組(10)。

(10)

利用井旁道彈性阻抗數(shù)據(jù)體和井曲線，計算出常系數(shù)a(θ1)、a(θ2)、a(θ3)、b(θ1)、b(θ2)、b(θ3)、c(θ1)、c(θ2)和c(θ3)，將它們代入公式(10)即可得到任意采樣點的楊氏模量[15-18]，為后續(xù)地應力預測奠定基礎。

1.3 基于曲率屬性的地應力預測方法

曲率是表示地層彎曲程度的量化參數(shù)，當構造曲面的曲率增加時，說明地層彎曲程度變大，該區(qū)域所受的應力增加，會增加產(chǎn)生裂縫的可能性。楊氏模量是表征巖石脆性的巖石力學參數(shù)，楊氏模量越高，巖石的脆性越強，越容易產(chǎn)生裂縫。Hunt等[8]提出將曲率與楊氏模量的乘積近似表示為地層應力，這是一種預測裂縫密度的擴展曲率方法。Roberts[1]描述了一個簡單褶皺層的應變，如圖1所示。

圖1 簡單巖層彎曲示意圖Fig.1 A simple layer bending diagram

假設層的原始長度為L0，層的頂部長度為L1，層的底部長度為L2，層的厚度為h，彎曲巖層的半徑為R，則該彎曲巖層的頂部應變e可表示為式(11)。

(11)

由曲率的定義可知，曲率是半徑R的倒數(shù)，可得應變與曲率的關系為式(12)。

(12)

其中K代表最大曲率[2]。由于最大正曲率是描述裂縫、撓曲、褶皺和斷層最有效的曲率屬性，直接反映地層隆起結構特征，因此筆者選用最大正曲率代替最大曲率用于計算地層應力。根據(jù)胡克定律可以得到曲率與應力的關系為式(14)。

(13)

將公式(13)簡化為：

σ=E·K

(14)

從公式(14)可以看出，應力可以近似等價為楊氏模量和曲率的乘積，這表明大曲率和高楊氏模量的區(qū)域具有較高的應力，且該區(qū)域容易產(chǎn)生裂縫，易于形成油氣聚集的空間和運移的通道。

基于曲率屬性的地應力預測方法分為3步：

1)根據(jù)基于局部tau-p變換的曲率屬性提取方法，得到儲層的曲率屬性體：先將地震數(shù)據(jù)進行希爾伯特變換得到它的復地震數(shù)據(jù)，然后利用得到的地震數(shù)據(jù)和復地震數(shù)據(jù)通過基于局部τ-p變換的傾角屬性提取方法得到三維數(shù)據(jù)體的傾角數(shù)據(jù)體p和q，最后將傾角數(shù)據(jù)體用于提取曲率屬性。由于公式(14)是基于凸起的褶皺層推導得到，且最大正曲率屬性是描述地層隆起構造最有效的曲率屬性，因此估算地應力的曲率屬性選擇最大正曲率屬性。

2)根據(jù)疊前地震反演方法提取儲層的楊氏模量數(shù)據(jù)體：通過對三個不同入射角的疊前地震數(shù)據(jù)分別進行反演得到對應的不同角度的彈性阻抗數(shù)據(jù)體，結合測井數(shù)據(jù)和AVO分析結果，利用公式(10)以及反演得到的各個角度的彈性阻抗數(shù)據(jù)體得到地層的楊氏模量。

3)根據(jù)楊氏模量、曲率與地應力的計算公式，將前面兩步得到的曲率數(shù)據(jù)體和楊氏模量數(shù)據(jù)體相乘得到地應力數(shù)據(jù)體。具體流程如圖2所示。

圖2 基于曲率屬性的地應力預測方法流程圖Fig.2 The flow chart of in-situ stress prediction method based on curvature attribute

2 應用實例

實際資料選自中國東部某泥頁巖裂縫型儲層工區(qū)，該工區(qū)儲層段中發(fā)育著以灰質(zhì)泥巖、油泥巖、油頁巖等巖性的泥質(zhì)烴源巖，其中泥頁巖富含灰質(zhì)成分和脆性成分，容易受構造活動的影響形成大量的裂縫，是油氣賦存最為有利的區(qū)域。工區(qū)內(nèi)部有三組斷層呈近東西方向發(fā)育，由于近北東向、北西向斷層的存在，它們被切割成縱多斷塊[20]。鑒于該工區(qū)復雜的斷裂構造，為了能夠尋找利于產(chǎn)生裂縫的區(qū)域，采用能夠刻畫斷層系統(tǒng)的曲率屬性和表征地層脆性特征的楊氏模量結合的方法得到地應力，由該方法得到的地應力屬性既能展現(xiàn)地層的構造特征又能指示地層中容易壓裂破碎的區(qū)域，能夠較好地評價儲層中容易產(chǎn)生裂縫的區(qū)域。

利用該工區(qū)的地震數(shù)據(jù)做實際應用，首先利用基于局部tau-p變換的傾角屬性提取方法，得到該工區(qū)的視傾角數(shù)據(jù)體p和q(圖3)。其次，利用得到的視傾角數(shù)據(jù)體p和q提取地層的曲率屬性(圖4)。然后通過基于彈性阻抗的疊前地震反演方法提取楊氏模量(圖5)。最后將楊氏模量與曲率相乘(應力-曲率屬性)近似表征地層應力，用于評價地層中易于產(chǎn)生裂縫的區(qū)域，預測地層的裂縫密度，指導頁巖儲層的水力壓裂開發(fā)(圖6)。

最大正曲率主要描述儲層凸起特征，圖4中高值的區(qū)域代表地層隆起的區(qū)域，該區(qū)域應變較大表明所受的應力較大，會增加裂縫生成的可能性，但缺少對巖石性質(zhì)變化的考慮。圖5是反演得到的楊氏模量層切片，主要展示地層巖石性質(zhì)的變化。從公式(14)可知，當曲率一定時，楊氏模量與應力成正比的關系，楊氏模量越高所受的應力越大，其高值表示該區(qū)域的巖石脆性較強，容易發(fā)生破裂，形成裂縫。當巖石的脆性不變時，曲率越大地層越彎曲，受到的應力越大，越容易產(chǎn)生裂縫，裂縫發(fā)育程度與彎曲程度成正比。用楊氏模量擴展曲率得到的應力-曲率屬性層切片，可以近似表征地層應力，既結合了地層的連續(xù)性又體現(xiàn)了地層的巖性變化，其高值區(qū)域表示該區(qū)域容易產(chǎn)生裂縫，低值區(qū)域表示該區(qū)域不易發(fā)生破裂，可作為預測儲層裂縫密度的可靠依據(jù)，指導頁巖儲層的水力壓裂開發(fā)。通過地應力的分布可以找到泥巖裂縫的有利分布區(qū)域，圖6中北部中間部分區(qū)域應力-曲率屬性值較大且較為集中，該區(qū)域所受的應力較大，地層容易受到擠壓作用產(chǎn)生裂縫，在此區(qū)域進行水力壓裂能夠獲得較多的裂縫，容易形成裂縫密集區(qū)域，而西部應力-曲率屬性較小，巖石不容易破裂，應力較小，不易產(chǎn)生裂縫，應避開類似區(qū)域進行水力壓裂。

圖3 視傾角屬性p和q層切片F(xiàn)ig.3 Horizontal slicer of apparent dip p and q (a)p層切片；(b)q層切片

圖4 最大正曲率屬性層切片F(xiàn)ig.4 Horizontal slicer of the most positive curvature

圖5 反演得到的楊氏模量層切片F(xiàn)ig.5 Horizontal slicer of Young’s modulus inverted result

圖6 應力-曲率屬性層切片F(xiàn)ig.6 Horizontal slicer of stress-curvature

3 結論

筆者根據(jù)預測裂縫密度的擴展曲率方法對某工區(qū)實際資料做了實際應用。利用基于tau-p變換的曲率屬性提取方法得到的最大正曲率與疊前地震反演方法得到的楊氏模量相乘得到應力-曲率屬性，將其用于評價因構造彎曲作用而產(chǎn)生的裂縫發(fā)育情況，指示儲層中容易產(chǎn)生裂縫的區(qū)域，實現(xiàn)曲率屬性在地應力地震預測中的應用。從實際資料可以看出,應力-曲率屬性不僅考慮了巖石性質(zhì)變化，且能夠較好地反映地層的構造形態(tài)，為裂縫預測提供了借鑒，地層的北部中間區(qū)域高值較多且分布較為集中，該區(qū)域的應力較大，容易產(chǎn)生較多的裂縫，易于形成裂縫密集區(qū)域，而地層的西部則所受應力較小，該區(qū)域的巖石不容易發(fā)生破裂，其在布置井位置和進行水力壓裂時可作為參考，但在判定裂縫的密集區(qū)域時仍需要其他資料進行驗證，綜合應用多種地質(zhì)資料判斷裂縫的發(fā)育區(qū)域，進行油氣的勘探與開發(fā)。

參考文獻：

[1] ROBERTS A. Curvature attributes and their application to 3D interpreted horizons[J]. First break, 2001, 19(2):85-100.

[2] MURRAY G. H..Quantitative fracture study—Spanish Pool, McKenzie County, North Dakota[J]. AAPG Bulletin, 1968, 52(1):57-65.

[3] PRICE N.J.，COSGROVE J.W. Analysis of geological structures[M]. Cambridge: Cambridge University Press, 1990.

[4] SHEOREY P R. A theory for in situ, stresses in isotropic and transverseley isotropic rock[J]. International Journal of Rock Mechanics & Mining Science & Geomechanics Abstracts, 1994, 31(1):23-34.

[5] 李志勇, 曾佐勛, 羅文強. 裂縫預測主曲率法的新探索[J]. 石油勘探與開發(fā), 2003, 30(6):83-85.

LI Z Y, ZENG Z X, LUO W Q. A new approach for predicting fractures using principal curvature[J]. Petroleum Exploration and Development, 2003, 30(6):83-85. (In Chinese)

[6] 唐湘蓉, 李晶. 構造應力場有限元數(shù)值模擬在裂縫預測中的應用[J]. 特種油氣藏, 2005, 12(2):25-28.

TANG X R, LI J. Application of finite element numerical simulation of tectonic stress field in fracture prediction[J]. Special Oil and Gas Reservoirs, 2005, 12(2):25-28. (In Chinese)

[7] 何英. 高精度曲率分析方法及其在構造識別中的應用[D]. 成都：成都理工大學, 2011.

HE Y. High precision curvature analysis and its application of structural identification[D]. Chengdu: Chengdu University of Technology, 2011. (In Chinese)

[8] HUNT L, REYNOLDS S, HADLEY S, et al. Causal fracture prediction: Curvature, stress, and geomechanics[J]. The Leading Edge, 2011, 30(11):1274-1286.

[9] STARR J.Modified curvature analysis to quantify strain within the Marcellus shale[C].SEG Technical Program Expanded Abstracts,2014:2373-2376.

[10] 王欣. 基于局部tau-p變換的地震幾何屬性提取及融合方法研究[D]. 青島：中國石油大學, 2016.

WANG X. The study of seismic geometry attribute and image fusion methods based on local tau-p transformation[D]. Qingdao: China University of Petroleum, 2016. (In Chinese)

[11] 王欣，印興耀，楊繼東，等. 基于局部tau-p變換的地震曲率屬性提取方法[C].中國地球物理年會第三十一屆年會，2015:1449-1451.

WANG X, YIN X Y, YANG J D,et al. The study of seismic attribute extraction methods based on local tau-p transformation[C]. CGU, 2015:1449-1451. (In Chinese)

[12] 印興耀，高京華，宗兆云. 基于離心窗傾角掃描的曲率屬性提取[J]. 地球物理學報, 2014,57(10):3411-3421.

YIN X Y, GAO J H, ZONG Z Y. Curvature attribute based on dip scan with eccentric window[J]. Chinese J. Geophys, 2014, 57(10):3411-3421. (In Chinese)

[13] 印興耀，周靜毅. 地震屬性優(yōu)化方法綜述[J]. 石油地球物理勘探, 2005,40(4):482-489.

YIN X Y, ZHOU J Y. Summary of optimum methods of seismic attributes[J]. Oil Geophysical prospectiong, 2005,40(4):482-489. (In Chinese)

[14] Al-DOSSARY S, MARFURT K J. 3D volumetric multispectral estimates of reflector curvature and rotation[J]. Geophysics, 2006, 71(5):41-51.

[15] 宗兆云，印興耀，張峰，等. 楊氏模量和泊松比反射系數(shù)近似公式及疊前地震反演[J].地球物理學報, 2012, 55(11):3786-3794.

ZONG Z Y, YIN X Y, ZHANG F, et al. Reflection coefficient equation and pre-stack seismic inversion with Young’s modulus and Poisson ratio[J]. Chinese J. Geophys, 2012, 55(11):3786-3794. (In Chinese)

[16] 宗兆云，印興耀，吳國忱. 基于疊前地震縱橫波模量直接反演的流體檢測方法[J]. 地球物理學報, 2012, 55(1):284-292.

ZONG Z Y, YIN X Y, WU G C. Fluid identification method based on compressional and shear modulus direct inversion[J]. Chinese J. Geophys, 2012, 55(1):284-292. (In Chinese)

[17] 印興耀，宗兆云，吳國忱. 巖石物理驅動下地震流體識別研究[J].中國科學:地球科學, 2015, 45(1):8-21.

ZONG Z Y, YIN X Y, WU G C. Research on seismic fluid identification driven by rock physics[J]. Science China: Earth Sciences, 2015, 45(1):8-21. (In Chinese)

[18] 印興耀，曹丹平，王保麗，等. 基于疊前地震反演的流體識別方法研究進展[J]. 石油地球物理勘探, 2014, 49(1):22-34,46.

YIN X Y, CAO D P, WANG B L, et al. Research progress of fluid discrimination with pre-stack seismic inversion[J]. Oil Geophysical prospectiong, 2014, 49(1):22-34,46. (In Chinese)

[19] AKI K, RICHARDS P G. Quantitative Seismology(2 ed)[C]. San Francisco: W. H. Freeman, 1980:158-164.

[20] 羅輯. 方位各向異性疊前地震反演方法研究[D]. 青島：中國石油大學, 2016.

LUO J. Study on methodology of pre-stack seismic inversion for azimuthal anisotropic medium[D]. Qingdao: China University of Petroleum, 2016. (In Chinese)