頁巖儲層地應力預測模型的建立和求解

鄧金根，陳崢嶸，，耿亞楠，劉書杰，朱海燕

(1.中國石油大學油氣資源與探測國家重點實驗室，北京102249;2.中海油研究總院，北京100027)

目前，常規能源的開發和利用已逐漸不能滿足現代工業的發展需要，勘探開發的重點投向了難以開采的頁巖氣。頁巖儲層是一種低孔隙度和特低滲透率的致密儲層，其生產完全依賴于壓裂效果。對頁巖進行壓裂施工必須對頁巖起裂的規律有足夠的認識，其中地應力模型是指導頁巖壓裂施工的極為重要的理論基礎。目前地應力預測模型大部分建立在地層均質各向同性的基礎上，而頁巖儲層具有較強的非均質性和各向異性。Thiercelin［1］討論了線彈性的橫觀各向同性地層地應力的計算模型，Higgins［2］在研究 Baxter頁巖地層的地應力時分別應用了各向同性和橫觀各向同性的地應力計算模型，并對二者進行了計算和比較。目前，多數學者的理論都只考慮上覆地層壓力和構造應力的影響，對地應力的描述并不全面，因此并不能完全描述橫觀各向同性地層地應力的特性。筆者在前人理論研究的基礎上，結合滲流和孔隙壓力的影響建立考慮頁巖儲層橫觀各向同性的地應力模型。

1 頁巖本構特性

頁巖由于其沉積歷史的特殊性［3］通常具有良好的分層結構，在構造和組成上呈連續變化，因此巖石的性質在平面和垂向上會有區別，頁巖巖石呈現出橫觀各向同性的特性［4-5］，其巖石力學性質也表現出橫觀各向同性［6］，即頁巖儲層水平方向是各向同性的而垂直方向則表現為各向異性(圖1)。因此，對頁巖儲層可通過橫觀各向同性模型進行研究。

圖1 頁巖橫觀各向同性單元體Fig.1 Element of shale representing transversely isotropy body

對于橫觀各向同性的頁巖，應力和應變之間的本構方程可通過修正后的胡克定律［7］表示為

式中，σij為二階應力張量的分量形式，MPa;εkl為二階應變張量的分量形式;Cijkl為四階彈性剛度張量的分量形式，MPa;i，j，k，l=1，2，3。

二階應力張量和二階應變張量實際分別都只有6個獨立的分量，因此四階彈性剛度張量可表示為6×6矩陣。根據地層的對稱性，剛度張量的獨立參數將減少，因此剛度張量C可表示為

對于橫觀各向同性地層，矩陣各參數有以下關系 C11=C22，C12=C21，C13=C31=C23=C32，C44=，其中 C66可以表示為因此可以用5個獨立的變量表示矩陣C［9-10］。要完全描述矩陣C需要求得以上5個獨立的參數，該5個參數［11］可以表示為

式中，ρ為體積密度，kg/m3;vPH和vPV分別為水平方向和垂直方向的P波波速，m/s;vSH和vSV分別為水平方向和垂直方向的S波波速，m/s;v45為與地層對稱軸成45°的P波(m=－1)或S波(m=1)波速，m/s。

經過聲波測井的測量，可求得橫觀各向同性的5個獨立參數。則垂直和水平方向上的彈性模量和泊松比［12-13］分別為

由彈性模量和泊松比表示的頁巖本構模型為

其中

式中，σ和τ分別為作用在單元體上的正應力和剪應力，MPa;ε和γ分別為單元體上的正應變和剪應變;Eh和Ev分別為水平方向和垂直方向上的彈性模量，MPa;νh和νv分別為施加水平和垂直應變時水平應變的泊松比;Gv為垂直平面的剪切模量，MPa。

地層力學參數的測量通常通過應力應變的關系或聲波速度測量得到，通過應力應變關系施加載荷測試，一般測定的是標準的力學彈性參數(如彈性模量、剪切模量和泊松比等)，聲波速度測試則是測量聲波的強度參數(如C11、C33和C12等)，而聲波實際上也是施加在巖石上的載荷，聲學性質與聲波傳播過程中巖石的應力應變有關。通過上述分析，這兩種力學參數的測量方法等價，通過以上兩種方法都可以解決巖石力學參數的測量問題，在應用過程中可根據實際情況擇優使用。

2 頁巖地應力預測模型的建立

地應力場中垂直地應力由上覆地層的重力產生，可通過地層的密度和厚度的積分計算得到。水平地應力的形成由地層構造運動產生的應變、不同地層之間地層性質的各向異性和非均質性及其他地質因素共同作用產生。為了預測水平地應力并構建理想的地應力模型，首先必須對所研究地層的物理力學參數的各向異性情況有足夠的了解。建立頁巖儲層地應力模型必須考慮頁巖的橫觀各向同性的特性。

Thiercelin討論了線彈性各向同性和橫觀各向同性的地層地應力計算模型，對于各向同性地層，在平面應變的假設下，考慮上覆地層壓力和構造應力的影響，水平地應力模型［1］可表示為

式中，σv為垂直地應力，MPa;σH和σh分別為水平最大和最小地應力，MPa;pp為地層孔隙壓力，MPa;α為有效應力系數，通常取1;εH和εh分別為水平最大和最小構造應變。

對于橫觀各向同性的地層，Thiercelin在各向同性地應力模型的基礎上考慮了地層垂直和水平方向上力學性質的差別，建立了橫觀各向同性地層地應力模型，即

式中，ξ為滲流系數，通常取0。

其他地質條件也可能產生水平地應力，如孔隙壓力的變化。結合上述因素，基于線彈性假設建立考慮地層橫觀各向同性、孔隙壓力變化、構造載荷和滲流的有效水平地應力模型，即

式中，σh，g，σh，p和 σh，t分別為上覆地層壓力、孔隙壓力變化和構造應力產生的水平地應力，MPa。

圖2所示為頁巖儲層中影響水平地應力的各種因素。

圖2 頁巖儲層中影響水平地應力的各種因素Fig.2 Factors of influencing horizontal in-situ stress in shale formation

如果忽略構造應力的影響，所產生的水平地應力將相同(σH=σh)，如果某一方向上存在構造作用產生的變形，則將導致水平地應力在不同方向上的應力差(σH＞σh)。相應的地應力為

其中

孔隙壓力的變化通過系數(1+K0)改變水平地應力。地質上的構造應變則通過壓縮和拉伸使地應力相應的增加和減少，二者之間的比例關系可用地層的各向異性參數K1和K2表示。對于橫觀各向同性地層，如果不存在構造應力產生的變形，則兩個水平地應力相等，而一旦考慮了構造應力，兩個水平地應力將產生差別。除以上因素外，還有一些影響地應力預測的地質力學因素。該模型相對簡單且所需的數據可通過測井解釋方法獲取，便于實際應用。

3 算 例

圖3 Eh/Ev對最大和最小水平地應力的影響Fig.3 Effect of Eh/Evon the maximum and minimum horizontal in-situ stress

根據上述理論分析，編寫計算機程序模擬Eh/Ev、νv/νh、孔隙壓力變化對地應力的影響。模擬參數:地層平均密度為2.61 g/cm3，有效應力系數為1，滲流系數為0，εH和 εh分別為0.3和0.1。

3.1 Eh/Ev對頁巖儲層地應力的影響

模擬參數:垂直和水平方向的泊松比均為0.25，孔隙壓力變化為0.5 MPa。圖3為Eh/Ev對最大和最小水平地應力的影響。隨著Eh/Ev的增大，地層的各向異性更強，水平地應力將整體變大。當Eh/Ev=4時，最大和最小水平地應力分別是各向同性(Eh/Ev=1)情況下的2.30和2.57倍，平均增加144%，因此忽略橫觀各向同性對地層地應力的影響，預測結果會產生較大誤差，各向同性的預測結果明顯偏低。

3.2 νv/νh對頁巖儲層地應力的影響

模擬參數:垂直和水平方向的彈性模量均為16 GPa，孔隙壓力變化為0.5 MPa。圖4 為 νv/νh對最大和最小水平地應力的影響。

圖4 vv/vh對最大和最小水平地應力的影響Fig.4 Effect of vv/vhon the maximum and minimum horizontal in-situ stress

從圖4中可以看出，最大和最小水平地應力隨νv/νh值增大而增大。其中，νv/νh=2.2 的最大和最小水平地應力分別是νv/νh=1時的1.15和1.18倍，平均增加16.5%。因此垂直和水平方向上的泊松比對地應力的預測亦會產生較大的影響。

3.3 孔隙壓力變化對頁巖儲層地應力的影響

圖5 孔隙壓力變化對最大和最小水平地應力的影響Fig.5 Effect of pore pressure on the maximum and minimum horizontal in-situ stress

模擬參數:垂直和水平方向的彈性模量均為16 GPa，垂直和水平方向的泊松比均為0.2。圖5為孔隙壓力變化對最大和最小水平地應力的影響。預測地應力考慮孔隙壓力的變化時，最大和最小水平地應力隨著孔隙應力變化的增大而增大。其中，孔隙壓力變化值為2.4 MPa時的最大和最小水平地應力分別是孔隙壓力不變化時的1.13和1.15倍，平均增加13.8%。因此，孔隙壓力的變化對地應力的預測也是一個不可忽視的因素。

3.4 地層升沉對頁巖儲層地應力的影響

一些學者在建立地應力預測模型時還考慮了地層的升沉作用(地層上升或者沉降)對地應力預測的影響［14-15］。該模型表達式為

對應的最大和最小水平地應力分別為

式中，σh，u為升沉作用影響產生的水平地應力，MPa;h為當前地層深度，m，Δh為當前地層深度與上升或者沉降之前地層深度的差值，即升沉高度，m。

模擬參數:垂直和水平方向的彈性模量均為16 GPa，垂直和水平方向的泊松比均為0.2，孔隙壓力變化為0.5 MPa，其他參數同上。圖6為升沉高度對最大和最小水平地應力的影響。最大和最小水平地應力隨著升沉高度的增大而增大，但升沉作用對地應力的預測影響不大。

圖6 升沉高度對最大和最小水平地應力的影響Fig.6 Effect of subsidence or uplift on the maximum and minimum horizontal in-situ stress

4 模型驗證

選取Baxter頁巖的地層參數進行模擬計算說明式(2)的準確性，將計算結果與地層瞬間閉合壓力(ISIP)進行比較(圖7)，可以看出本文模型和Thiercelin模型的結果都與ISIP基本一致，本文模型結果誤差最大僅為1.2%，所得地應力符合Baxter頁巖地應力的計算情況［2］，驗證了所建立的地應力模型的準確性。相同條件下，各向同性模型計算的地應力結果與實際誤差較大，并不準確，因此采用橫觀各向同性地應力模型更符合現場實際情況。

圖7 各向同性模型、Thiercelin模型和本文模型計算的最小水平地應力與ISIP數據的比較Fig.7 Comparison of ISIP and the minimum horizontal in-situ stress computed by isotropy，Thiercelin and transversely isotropy models

5 結論

(1)水平與垂直彈性模量之比、垂直與水平泊松比的比值和孔隙壓力的變化對地應力的預測影響較大，且以上3種因素增大，地應力的預測結果也將增大。升沉作用對地應力的影響較小。

(2)對于頁巖儲層，橫觀各向同性的地應力模型要比各向同性的地應力模型更能表示地層地應力的真實情況。準確計算地應力的分布對頁巖儲層現場的完井壓裂增產設計有重要的指導作用。

［1］THIERCELIN M J，PLUMB R A.A core-based prediction of lithologic stress contrasts in east Texas formations［J］.SPE Formation Evaluation，1994，9(4):251-258.

［2］HIGGINS S，GOODWIN S，DONALD A，et al.Anisotropic stress models improve completion design in the Baxter shale［R］.SPE 115736，2008.

［3］RIVERA S R，HANDWERGER D，KIESCHNICK J，et al.Accounting for heterogeneity provides a new perspective for completions in tight gas shales［R］.American Rock Mechanics Association 05-758，2005.

［4］LINGS M L，PENNIGNTON D S，DASH D F T.Anisotropic stiffness parameters and their measurement in a stiff natural clay［J］.Geotechnique，2000，50(2):109-125.

［5］GRAHAM J，CROOKS J H A，BELL A L.Time effects on the stress-strain behavior of natural soft clays ［J］.Geotechnique，1983，33(3):327-340.

［6］SAYERS C M.Seismic anisotropy of shales［J］.Geophysical Prospecting，2005，53(5):667-676.

［7］ GAUTAM R.Anisotropy in deformations and hydraulic properties of Colorado shale［D］.Calgary:Civil Engineering，University of Calgary，2004.

［8］NYE J F.Physical properties of crystals［M］.Oxford:Oxford University Press，1985.

［9］NORRIS A N，SINHA B K.Weak elastic anisotropy and the tube wave［J］.Geophysics，1993，58(8):1091-1098.

［10］WALSH J，SINHA B，DONALD A.Formation anisotropy parameters using borehole sonic data［R］.Society of Petrophysicists and Well-Log Analysis 2006-TT，2006.

［11］HORNBY B E，HOWIE J M，INCE D W.Anisotropy correction for deviated-well sonic logs:application to seismic well tie［J］.Geophysics，2003，68(2):464-471.

［12］SAROUT J，MOLEZ L，GUEGUEN Y，et al.Shale dynamic properties and anisotropy under triaxial loading:experimental and theoretical investigations［J］.Physical and Chemistry of the Earth，2007，32(8/14):896-906.

［13］RIVERA S R，DEENDAYALU C，YANG Y K.Unlocking the unconventional oil and gas reservoirs:the effect of laminated heterogeneity in wellbore stability and completion of light gas shale reservois［R］.Offshore Technology Conference 20269，2009.

［14］SAFDAR K，SAJJAD A，HAN H X，et al.Importance of shale anisotropy in estimating in-situ stresses and wellbore stability analysis in Horn river basin［R］.SPE 149433，2011.

［15］PRIOUL R，KARPFINGER F，DEENADAYALU C，et al.Improving fracture initiation predictions on arbitrarily oriented wells in anisotropic shales［R］.SPE 147462，2011.

［16］RIVERA S R，DEENADAYALU C，CHERTOV M，et al.Improving horizontal completions on heterogeneous tight shales［R］.SPE 146998，2011.