韓城區塊西北部三維地應力場反演分析

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

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

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

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

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

1 研究區域概況

韓城區塊地質構造為北東走向、北西傾向的單斜構造。位于鄂爾多斯盆地渭北隆起東北部，渭北隆起自三疊紀后經歷了印支運動期NNW 向擠壓、燕山運動中期NWW 向擠壓和喜馬拉雅運動期NE 向擠壓與NW-SE 向伸展作用三期構造作用。印支、燕山、喜馬拉雅三期構造運動疊加形成該區域構造形態，其中燕山運動奠定了本區的構造基礎。

因此，韓城區塊擠壓構造帶主要發育在東南邊淺部地層，雖然伸展構造遍及全區，但主要集中在文家嶺隆起以南地區，而且發育規模及強度從南向北、從邊淺部向中深部慢慢變弱。南北區域構造發育的差異性表現在南區以斷層為主，褶皺發育較少，北區則相反。在區域應力擠壓的作用下，斷層較為發育，受3 條斷裂帶的影響，韓城區塊煤層氣的富集差異大。研究區塊構造綱要如圖1。

圖1 研究區塊構造綱要圖Fig.1 Outline diagram of research block structure

2 地應力場實測數據分析

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

2.1 水力壓裂法原理

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

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

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

計算最大、最小主應力的公式如下[24]：

式中：ps為 瞬時關閉壓力；pr為裂縫重壓力。

垂向應力Sv可以通過地層的密度計算得出。

2.2 測點布置

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

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

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

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

首先，計算實測主應力與建立的坐標軸之間的方向余弦[25]經進一步推導得出：

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

根據所測得的主應力的大小和方向余弦，再計算主應力轉換坐標系下的應力分量[26]：

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

坐標系下實測地應力分量值見表2。

表2 坐標系下實測地應力分量值Table 2 Measured geostress component values in the coordinate system MPa

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

3 地應力場三維有限元回歸分析

3.1 地應力場回歸分析原理

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

1）根據研究區域的地質資料，建立起精確的有限元數值模型。

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

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

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

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

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

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

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

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

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

如果直接在模型加載應力邊界條件，會產生邊界應力集中，與實際應力情況有很大差別，且煤層和頂底板砂巖層起伏較大，直接加載應力也會產生較大誤差。

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

3.3 初始地應力場反演回歸結果

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

利用式（15）即可得到測點位置處的應力回歸計算值，通過數學關系轉換得到對應的最大和最小水平主應力，并與應力實測數據進行對比分析。實測數據與回歸計算值對比見表3。

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

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

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

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

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

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

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

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

從圖7 可以看出，除了中部的高應力區外，整體的應力方向趨勢呈N45°E。在這2 個高應力區的應力方向由水平向突變為豎直向，初步分析認為是此處的最大水平應力超過了垂向應力，而繪制這個應力方向矢量圖用的是水平方向和傾角，此時的中間應力由原本的最大水平應力方向變為了垂向應力，所以會發生應力突變。

4 結 語

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

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