趙海峰 ，楊昌增 ，馮 堃 ，甄懷賓

（1.中國石油大學(xué)（北京） 石油工程學(xué)院，北京 102249；2.中國石油塔里木油田分公司，新疆 庫爾勒 841901；3.中石油煤層氣有限責(zé)任公司，北京 100028；4.中聯(lián)煤層氣國家工程研究中心有限責(zé)任公司，北京 100095）

煤層氣是一種清潔高效的能源，加快對煤層氣的開發(fā)與利用能夠有效緩解國內(nèi)的能源供給壓力[1-2]?，F(xiàn)有煤層氣資源勘探調(diào)查結(jié)果顯示，煤層氣全球的可采儲量約為256.3×1012m3，我國的煤層氣資源含量十分豐富[3]。研究表明，韓城區(qū)塊埋深在1 500 m 以內(nèi)的煤層氣儲量達(dá)到2 007×108m3，埋 深1 500～2 000 m 之 間 的 煤 層 氣 儲 量 達(dá) 到 了934×108m3。這些數(shù)據(jù)表明，韓城區(qū)塊具有較大的煤層氣開采價值，因此，它是我國深層煤層氣開發(fā)的重點對象[4]。

研究顯示，現(xiàn)今地應(yīng)力的分布情況和煤層滲透率的大小決定著煤層氣勘探開發(fā)的成功與否[5-7]，對煤層氣的高效經(jīng)濟(jì)開采具有重要的指導(dǎo)意義。在實施煤層頂板壓裂的過程中，地應(yīng)力控制著壓裂裂縫從起裂到擴(kuò)展再到延伸的整個過程，從而決定了壓裂的成敗[8]。同時，地應(yīng)力對煤層的滲透率也有著較大的影響[9]。地應(yīng)力的確定是進(jìn)行煤層氣開采的基礎(chǔ)，目前通過鉆孔進(jìn)行地下應(yīng)力的測量是一種公認(rèn)的最為可靠的應(yīng)力測量方式，但由于受地質(zhì)環(huán)境和條件的限制，無法實現(xiàn)大面積的地下應(yīng)力測量，而在實際的工程應(yīng)用中需要的是整個區(qū)塊的地應(yīng)力分布[10-11]。

因此，根據(jù)研究區(qū)塊的地質(zhì)資料和實測數(shù)據(jù)，建立準(zhǔn)確地計算模型并進(jìn)行地應(yīng)力場反演，對后續(xù)的煤層氣開發(fā)具有重要的指導(dǎo)意義[12-13]。近年來，不同學(xué)者提出了多種用于地應(yīng)力反演分析的方法[14-19]，其中，多元線性回歸法原理簡單且解唯一，被眾多學(xué)者所沿用[20-23]。鑒于此，依據(jù)韓城北部某區(qū)域的地應(yīng)力實測數(shù)據(jù)，采用三維有限元回歸分析法對研究區(qū)進(jìn)行應(yīng)力場反演，進(jìn)一步分析，獲得了該研究區(qū)目前的地應(yīng)力場分布規(guī)律，為該區(qū)域的煤層氣的開采和增產(chǎn)提供了重要的數(shù)據(jù)基礎(chǔ)。

1 研究區(qū)域概況

韓城區(qū)塊地質(zhì)構(gòu)造為北東走向、北西傾向的單斜構(gòu)造。位于鄂爾多斯盆地渭北隆起東北部，渭北隆起自三疊紀(jì)后經(jīng)歷了印支運(yùn)動期NNW 向擠壓、燕山運(yùn)動中期NWW 向擠壓和喜馬拉雅運(yùn)動期NE 向擠壓與NW-SE 向伸展作用三期構(gòu)造作用。印支、燕山、喜馬拉雅三期構(gòu)造運(yùn)動疊加形成該區(qū)域構(gòu)造形態(tài)，其中燕山運(yùn)動奠定了本區(qū)的構(gòu)造基礎(chǔ)。

因此，韓城區(qū)塊擠壓構(gòu)造帶主要發(fā)育在東南邊淺部地層，雖然伸展構(gòu)造遍及全區(qū)，但主要集中在文家?guī)X隆起以南地區(qū)，而且發(fā)育規(guī)模及強(qiáng)度從南向北、從邊淺部向中深部慢慢變?nèi)?。南北區(qū)域構(gòu)造發(fā)育的差異性表現(xiàn)在南區(qū)以斷層為主，褶皺發(fā)育較少，北區(qū)則相反。在區(qū)域應(yīng)力擠壓的作用下，斷層較為發(fā)育，受3 條斷裂帶的影響，韓城區(qū)塊煤層氣的富集差異大。研究區(qū)塊構(gòu)造綱要如圖1。

圖1 研究區(qū)塊構(gòu)造綱要圖Fig.1 Outline diagram of research block structure

2 地應(yīng)力場實測數(shù)據(jù)分析

實測地應(yīng)力數(shù)據(jù)是現(xiàn)今地應(yīng)力場分析極其重要的數(shù)據(jù)基礎(chǔ)，在地應(yīng)力反演過程中實測數(shù)據(jù)亦是必不可少的重要參數(shù)。因此，在進(jìn)行地應(yīng)力反演前需進(jìn)行地應(yīng)力實測數(shù)據(jù)分析。通過鉆孔進(jìn)行地下應(yīng)力的測量是目前一種公認(rèn)的最為可靠的應(yīng)力測量方式，但是現(xiàn)場應(yīng)力實測需要耗費(fèi)巨額的資金，并且受地質(zhì)環(huán)境和技術(shù)條件的限制，無法實現(xiàn)大面積地下應(yīng)力測量。因此，利用水壓致裂法獲得的6 個測點的實測應(yīng)力數(shù)據(jù)對研究區(qū)塊的應(yīng)力場進(jìn)行反演分析。

2.1 水力壓裂法原理

向地層封隔段中注入高壓流體，記錄地層破裂的整個過程，通過壓裂曲線得出垂直于鉆孔平面的最小水平主應(yīng)力的大小，然后通過經(jīng)驗公式算出最大水平主應(yīng)力。彈性力學(xué)是水壓致裂法測量原理的基礎(chǔ)，須滿足以下3 個要求：①巖石是線性、均勻、各向同性的彈性體；②巖石是完整的且不具有滲透性，注入的流體符合達(dá)西定律；③鉆孔軸方向與垂向主應(yīng)力方向相同。

可以用3 個應(yīng)力分量來表示地殼中某一點的應(yīng)力狀態(tài)，而鉆孔中某處的應(yīng)力狀態(tài)可以通過壓裂曲線計算得出。在壓裂的過程中，封隔段內(nèi)的液柱壓力會不斷增大，壁面上的切向應(yīng)力將隨之變小，最后慢慢演變?yōu)閺垜?yīng)力。當(dāng)壓力超過巖石的抗拉強(qiáng)度 σt時，巖石會發(fā)生破裂。此時的液體壓力被稱之為破裂壓力，其表達(dá)式為[24]：

式中：pb為 破裂壓力；Sh為鉆孔破裂面上最小主應(yīng)力，MPa；SH為鉆孔破裂面上最大主應(yīng)力，MPa； σt為 巖石的抗拉強(qiáng)度；p0為孔隙壓力。

計算最大、最小主應(yīng)力的公式如下[24]：

式中：ps為 瞬時關(guān)閉壓力；pr為裂縫重壓力。

垂向應(yīng)力Sv可以通過地層的密度計算得出。

2.2 測點布置

研究區(qū)域的地應(yīng)力測量共選擇6 個勘探鉆孔作為應(yīng)力測試孔。所有測點均布置在煤層頂板砂巖層中。除6 號鉆孔外，其他鉆孔深度均未超過1 000 m，最大測點深度為1 003 m。在水壓致裂測量試驗結(jié)束之后進(jìn)行印模定向測試，所測得的破裂面的走向就是最大水平主應(yīng)力方向。測點布置如圖2，鉆孔水壓致裂法地應(yīng)力測量結(jié)果見表1。

表1 鉆孔水壓致裂法地應(yīng)力測量結(jié)果Table 1 Measurement results of in-situ stress by drilling hydraulic fracturing method

圖2 測點布置圖Fig.2 Measuring points layout

以正東向為x軸正向，正北向為y軸的正向，豎直向上為z軸的正向，建立三維直角坐標(biāo)系。考慮到進(jìn)行地應(yīng)力場反演分析的過程中都是以應(yīng)力分量的形式來表征地應(yīng)力進(jìn)行回歸計算，而水壓致裂法測得的地應(yīng)力值是以主應(yīng)力的平面方位和傾角給出，所以為了便于計算，將根據(jù)方向余弦將實測地應(yīng)力值進(jìn)行坐標(biāo)變換。

首先，計算實測主應(yīng)力與建立的坐標(biāo)軸之間的方向余弦[25]經(jīng)進(jìn)一步推導(dǎo)得出：

式中：Li、Mi、Ni分別為 σi對x、y、z軸的方向余弦； σi為 實測主應(yīng)力； βi為 σi與水平面之間的夾角； αix為 σi與x軸正方向間的夾角。

根據(jù)所測得的主應(yīng)力的大小和方向余弦，再計算主應(yīng)力轉(zhuǎn)換坐標(biāo)系下的應(yīng)力分量[26]：

式中：σ1、σ2、σ3分 別為原始坐標(biāo)系下x、y、z方向的主應(yīng)力； σx、σy、σz分別為x、y、z方向的主應(yīng)力； τxy、τxz分別為與y向力相關(guān)的剪應(yīng)力和與z向力相關(guān)的剪應(yīng)力； τyz為σy與z向力相關(guān)的剪應(yīng)力。

坐標(biāo)系下實測地應(yīng)力分量值見表2。

表2 坐標(biāo)系下實測地應(yīng)力分量值Table 2 Measured geostress component values in the coordinate system MPa

綜合上述6 個測點的實測地應(yīng)力資料可以得出，地應(yīng)力狀態(tài)主要表現(xiàn)為Sv>SH>Sh，最小水平主應(yīng)力小于16 MPa，三向主應(yīng)力的這種分布關(guān)系表明現(xiàn)今地應(yīng)力處于拉伸應(yīng)力狀態(tài)，工程類比判斷韓城區(qū)塊的最大主應(yīng)力方向為N45°E，2 個鉆孔水壓致裂法實測均為N45°E 上下浮動3°～5°，最大主應(yīng)力方向基本一致。根據(jù)現(xiàn)場測井解釋結(jié)果顯示韓城區(qū)塊的平均地應(yīng)力梯度為2.56 MPa/hm；6 個測點埋深范圍為754～1 002 m 之間，測得的垂向應(yīng)力基本符合測井解釋的結(jié)果。

3 地應(yīng)力場三維有限元回歸分析

3.1 地應(yīng)力場回歸分析原理

為了全面掌握研究區(qū)域的現(xiàn)今地應(yīng)力場分布規(guī)律，引用有限元應(yīng)力回歸分析法進(jìn)行計算。該方法基于數(shù)值分析理論，通過計算機(jī)技術(shù)手段，把現(xiàn)有的地質(zhì)資料和現(xiàn)場實測地應(yīng)力資料及數(shù)理統(tǒng)計理論結(jié)合起來，充分考慮和分析影響地應(yīng)力的因素，保證反演出的應(yīng)力場最大限度地接近實際地層應(yīng)力。其基本原理是：

1）根據(jù)研究區(qū)域的地質(zhì)資料，建立起精確的有限元數(shù)值模型。

2）根據(jù)現(xiàn)今地應(yīng)力成因理論，選定垂直方向的自重因素和水平方向的2 個構(gòu)造因素作為待定因素，以單位載荷分別作用于先前建好的數(shù)值模型，得到各因素在單位載荷下的有限元計算值。

3）根據(jù)應(yīng)力線性疊加原理，將地應(yīng)力回歸值σ0n作為因變量，把有限元計算求得的單因素的地應(yīng)力計算值σin作為自變量，則回歸方程的形式為[27]：

最小二乘法殘差的平方和Q為[27]：

應(yīng)力分量m取1～6，分別對應(yīng)于6 個分量，根據(jù)最小二乘法的原理，求解使殘差的平方和取得最小值的方程式[27]，即可得到k個回歸系數(shù)D1、D2、···、Dk，將其代入式（13）即可求得模型區(qū)域任意位置的地應(yīng)力回歸計算值。

3.2 三維有限元計算模型的建立

為了提高有限元力學(xué)模型的精度，充分利用現(xiàn)有測井資料，在Petrel 中建立地質(zhì)模型，并通過測井曲線將泊松比和彈性模量粗化到整個網(wǎng)格中，最終在Fracman 中完成地應(yīng)力反演過程。建模思路如圖3。

圖3 建模思路圖Fig.3 Modeling idea diagram

首先，從研究區(qū)域的井軌跡和分層數(shù)據(jù)中得出各個層段（頂板頂部、煤層頂部、煤層底部、底板底部）的各層面點數(shù)據(jù)；然后通過Petrel 的make surface 模塊進(jìn)行曲面擬合，形成三維空間的曲面，通過三維建模模塊生成三維網(wǎng)格模型；最后，將測井曲線中的泊松比和彈性模量粗化至每個網(wǎng)格中。將得到的模型導(dǎo)入Fracman 中，通過裂縫建模生成煤層的割理裂隙，并將割理裂隙通過Oda地應(yīng)力計算粗化至每個網(wǎng)格中以彈性模量和泊松比的變化體現(xiàn)出來[23]。該模型計算域為：東西向分別以19 433 040 m 和19 437 060 m 坐標(biāo)網(wǎng)格為邊界，南北向分別以3 933 605 m 和3 936 485 m坐標(biāo)網(wǎng)格為邊界，東西向?qū)? 020 m，南北向長2 880 m，頂板下表面深度取到海拔-245 m，上表面深度取-28 m，計算域內(nèi)共劃分178 640 個六面體單元，三維有限元計算模型如圖4。

圖4 三維有限元計算模型Fig.4 3D finite element calculation model

如果直接在模型加載應(yīng)力邊界條件，會產(chǎn)生邊界應(yīng)力集中，與實際應(yīng)力情況有很大差別，且煤層和頂?shù)装迳皫r層起伏較大，直接加載應(yīng)力也會產(chǎn)生較大誤差。

由于自重作用和構(gòu)造運(yùn)動是引起地應(yīng)力的主要因素，因此，模擬的邊界條件只考慮這2 個因素。構(gòu)造應(yīng)力場力學(xué)模型為：2 個側(cè)面分別施加x方向的I1的單位均布載荷和y方向的J1單位應(yīng)力載荷，來模擬水平方向構(gòu)造運(yùn)動作用；在垂直方向加載K1的單位應(yīng)力載荷來模擬上覆地層壓力；在相應(yīng)的未加載應(yīng)力的對側(cè)分別施加x、y、z方向的位移約束。

3.3 初始地應(yīng)力場反演回歸結(jié)果

基于6 個測點的測試數(shù)據(jù)和3 個地應(yīng)力場計算結(jié)果，最終得到對應(yīng)于自重應(yīng)力場 σz、水平方向構(gòu)造應(yīng)力場 σx和 σy的3 個回歸系數(shù)，分別為D1=20.2，D2=18.5，D3=12.5。因此，得到本次研究區(qū)的地應(yīng)力場σ0回歸方程為：

利用式（15）即可得到測點位置處的應(yīng)力回歸計算值，通過數(shù)學(xué)關(guān)系轉(zhuǎn)換得到對應(yīng)的最大和最小水平主應(yīng)力，并與應(yīng)力實測數(shù)據(jù)進(jìn)行對比分析。實測數(shù)據(jù)與回歸計算值對比見表3。

表3 實測數(shù)據(jù)與回歸計算值對比Table 3 Comparison of measured data and regression calculated values

由表3 可知，最大、最小水平主應(yīng)力和垂向應(yīng)力都表現(xiàn)出隨著埋深的增加而逐漸增大的趨勢，且三者的大小關(guān)系基本保持為Sv>SH>Sh。現(xiàn)今地應(yīng)力處于拉伸狀態(tài)，本次反演結(jié)果與現(xiàn)場實測結(jié)果的應(yīng)力特征基本保持一致。實測值與回歸計算值之間的誤差范圍在0～5 MPa 之間。其中，垂向應(yīng)力的誤差較小基本在0～1 MPa 之間浮動。這是因為垂向應(yīng)力與地層的深度基本呈線性變化，擬合結(jié)果較為準(zhǔn)確；水平應(yīng)力的誤差偏大，在2～5 MPa 之間波動。由于水平方向的應(yīng)力成因復(fù)雜，單純用一種線性載荷來模擬是不夠嚴(yán)謹(jǐn)?shù)?，因此誤差偏大，但也在工程應(yīng)用可接受最大不超過25%的范圍之內(nèi)。因此，本次反演結(jié)果能夠滿足工程精度的要求，反演得到的地應(yīng)力場可表達(dá)研究區(qū)地應(yīng)力場的主要特征。

煤層頂面水平最大主應(yīng)力大小云圖如圖5，煤層頂面水平應(yīng)力差云圖如圖6，煤層頂面水平最大主應(yīng)力方向如圖7。

圖5 煤層頂面水平最大主應(yīng)力大小云圖Fig.5 Cloud diagram of horizontal maximum principal stress on top surface of coal seam

圖6 煤層頂面水平應(yīng)力差云圖Fig.6 Cloud diagram of horizontal stress difference on top surface of coal seam

圖7 煤層頂面水平最大主應(yīng)力方向Fig.7 Direction of maximum horizontal principal stress on top surface of coal seam

從圖5 可以看出，煤層的水平最大主應(yīng)力整體趨勢還是隨深度的增加而增大，變化范圍在10.7～25 MPa 之間，基本中部、北西向高，東部低的趨勢；在研究區(qū)中部有高應(yīng)力區(qū)的存在，這是可能是由于此處處于地層的褶皺區(qū)，在頂部受到了較大的拉伸應(yīng)力造成的。

從圖6 可以看出，煤層的整體應(yīng)力差處于1～5 MPa 之間，在北西、西南和中南部存在較大的應(yīng)力差，結(jié)合前面的應(yīng)力分析，其原因可能是此處的地層起伏較大，構(gòu)造應(yīng)力較為復(fù)雜，導(dǎo)致水平最大主應(yīng)力的大小超過了垂向應(yīng)力，出現(xiàn)了高應(yīng)力區(qū)。所以，建議在前期布井時應(yīng)該盡量避免這3 個高應(yīng)力差的地區(qū)，盡量選擇在東北部布井。

從圖7 可以看出，除了中部的高應(yīng)力區(qū)外，整體的應(yīng)力方向趨勢呈N45°E。在這2 個高應(yīng)力區(qū)的應(yīng)力方向由水平向突變?yōu)樨Q直向，初步分析認(rèn)為是此處的最大水平應(yīng)力超過了垂向應(yīng)力，而繪制這個應(yīng)力方向矢量圖用的是水平方向和傾角，此時的中間應(yīng)力由原本的最大水平應(yīng)力方向變?yōu)榱舜瓜驊?yīng)力，所以會發(fā)生應(yīng)力突變。

4 結(jié) 語

1）通過分析研究區(qū)域的構(gòu)造特點、實測數(shù)據(jù)和已有的研究成果，經(jīng)多次的反演計算之后得出，在垂向施加20 MPa、東西向施加18.5 MPa、南北向施加12.5 MPa 應(yīng)力，所取得的反演效果最好。

2）測點位置處的三向應(yīng)力大小均表現(xiàn)出隨著埋深的增加而逐漸增大的趨勢，而現(xiàn)場實測應(yīng)力由部分與回歸計算值有較大差異。這由于是測點的位置在煤層與頂板的交界處，該處的地質(zhì)參數(shù)不穩(wěn)定從而影響地應(yīng)力大小產(chǎn)生不規(guī)律的變化。3）煤層的水平最大主應(yīng)力處于10～25 MPa 之間，呈現(xiàn)中部、北西向高，東部低的趨勢；煤層的整體應(yīng)力差處于1～5 MPa 之間，在西、西南和中南部有較大應(yīng)力差存在；究區(qū)整體應(yīng)力方向趨勢呈N45°E。建議在前期布井階段應(yīng)盡量避免這3 個高應(yīng)力差的地區(qū)，盡量選擇在東北部布井。