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

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

0 引言

由于全球?qū)δ茉吹男枨蟛粩嘣黾忧页Ｒ?guī)油氣資源的不斷耗竭，頁巖氣藏的勘探和開發(fā)逐漸成為專家學(xué)者關(guān)注的熱點，而開發(fā)頁巖氣藏的有效手段是對其進(jìn)行水力壓裂。通過水力壓裂施工可以使儲層形成大量的裂縫，有利于油氣的儲存和運移。研究表明，裂縫密度與許多地層參數(shù)有關(guān)，如曲率、楊氏模量、地應(yīng)力和各向異性參數(shù)等。通過對地層參數(shù)的研究，可以有效地預(yù)測儲層的裂縫密度，指導(dǎo)頁巖氣的勘探與開發(fā)。

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

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

1 方法原理

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

曲率屬性是描述地層彎曲程度的一種地震幾何屬性，它能夠刻畫地層的撓曲、褶皺、斷裂以及河道等地質(zhì)構(gòu)造特征，是間接預(yù)測裂縫密度的一種屬性[12-13]。王欣[10-11]提出了基于局部τ-p變換的曲率屬性提取方法，該方法是將(t,x)域的地震數(shù)據(jù)及其復(fù)地震數(shù)據(jù)根據(jù)一定的線性疊加關(guān)系變換至(τ,p)域，應(yīng)用p和τ來描述波的軌跡。由于地震波場為點源波場，且遠(yuǎn)源波場的局部小地震道集可以近似看作是局部平面波，因此對這種局部小道集進(jìn)行τ-p變換稱為局部τ-p變換。利用局部τ-p變換可以將地震記錄上以時間延遲τ0為截距，以p0為斜率的同相軸(局部平面波)轉(zhuǎn)變成τ-p域中的一個點，其物理意義為遠(yuǎn)源球面波場的局部平面波分解。對于一個三維數(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ù)體通過傅氏變換從空間域轉(zhuǎn)換到波數(shù)域，將其乘以-ik因子相當(dāng)于其在空間域的一階導(dǎo)數(shù)的傅氏變換，并引入了分波數(shù)導(dǎo)數(shù)的概念[12]，可以得到多尺度曲率屬性的計算公式：

(3)

(4)

最大正曲率表示地層凸起構(gòu)造，因此選用最大正曲率作為地應(yīng)力預(yù)測的基礎(chǔ)。

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

宗兆云等[15-18]在Aki等[19]近似方程的基礎(chǔ)上，建立了縱波反射系數(shù)與楊氏模量反射系數(shù)、泊松比反射系數(shù)及密度反射系數(shù)的線性關(guān)系，稱之為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ù)地應(yīng)力預(yù)測奠定基礎(chǔ)。

1.3 基于曲率屬性的地應(yīng)力預(yù)測方法

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

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

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

(11)

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

(12)

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

(13)

將公式(13)簡化為：

σ=E·K

(14)

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

基于曲率屬性的地應(yīng)力預(yù)測方法分為3步：

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

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

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

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

2 應(yīng)用實例

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

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

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

圖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 應(yīng)力-曲率屬性層切片F(xiàn)ig.6 Horizontal slicer of stress-curvature

3 結(jié)論

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

參考文獻(xiàn)：

[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] 李志勇, 曾佐勛, 羅文強. 裂縫預(yù)測主曲率法的新探索[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] 唐湘蓉, 李晶. 構(gòu)造應(yīng)力場有限元數(shù)值模擬在裂縫預(yù)測中的應(yīng)用[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] 何英. 高精度曲率分析方法及其在構(gòu)造識別中的應(yīng)用[D]. 成都：成都理工大學(xué), 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]. 青島：中國石油大學(xué), 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]. 地球物理學(xué)報, 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].地球物理學(xué)報, 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]. 地球物理學(xué)報, 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] 印興耀，宗兆云，吳國忱. 巖石物理驅(qū)動下地震流體識別研究[J].中國科學(xué):地球科學(xué), 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ìn)展[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]. 青島：中國石油大學(xué), 2016.

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