梁利喜, 高 陽, 劉向君, 戶海勝, 熊 健, 鐘自強

(1.西南石油大學(xué), 油氣藏地質(zhì)與開發(fā)工程國家重點實驗室, 成都 610500; 2.中國石油新疆油田分公司勘探開發(fā)研究院, 克拉瑪依 834000)

某地區(qū)百口泉組砂礫巖地層巖性油氣藏顯示出巨大的勘探潛力,成為西部某油田重要的增儲上產(chǎn)領(lǐng)域[1]。由于對該地區(qū)地應(yīng)力特征認識程度較低,在油氣開發(fā)過程中頻現(xiàn)儲層壓裂改造效果差異大、水力壓裂裂縫擴展延伸規(guī)律難以預(yù)測以及水平井井筒相繼出現(xiàn)砂橋而堵塞流體通道等問題,因此迫切需要進一步研究該地區(qū)的地應(yīng)力場狀態(tài),為完善井網(wǎng)部署、優(yōu)化壓裂設(shè)計、防治采油出砂等提供依據(jù)[2-3]。

地應(yīng)力是自然狀態(tài)下存在于巖體內(nèi)部的應(yīng)力,由巖體自重、地質(zhì)構(gòu)造作用及溫度壓力等共同作用下形成。隨著油氣勘探開發(fā)深度的加大,地質(zhì)環(huán)境越來越復(fù)雜,地應(yīng)力大小及方向?qū)τ蜌饩こ獭⒂蜌饪碧介_發(fā)等各個領(lǐng)域的重要影響越來越受到工業(yè)界的高度關(guān)注[4-5]。目前為止,中外學(xué)者研究地應(yīng)力的方法,主要有構(gòu)造形跡分析法、室內(nèi)巖心測試技術(shù)、礦場實測分析、測井曲線計算法、物理和數(shù)值模擬技術(shù)等,并朝著多種方法綜合應(yīng)用、相互印證方面發(fā)展完善[6]。隨著計算機科學(xué)的迅猛發(fā)展,數(shù)值模擬技術(shù)在地應(yīng)力場研究中的應(yīng)用愈加廣泛,其能夠較好地綜合考慮深部地層巖石物理力學(xué)性質(zhì)復(fù)雜多變、地層介質(zhì)各向異性、非均質(zhì)性、不連續(xù)性和復(fù)雜邊界條件等因素,該方法受到了學(xué)者們的青睞[7-10]。目前,大多數(shù)的地應(yīng)力場數(shù)值模擬反演文獻的技術(shù)路線都是以單點地應(yīng)力實測資料為約束與修正依據(jù),通過試算法或者人工智能優(yōu)化方法(BP神經(jīng)網(wǎng)絡(luò)、粒子群算法等)確定應(yīng)力或位移狀態(tài)邊界條件,進而獲得地應(yīng)力場。然而在非均質(zhì)性較強的砂礫巖地層中,地應(yīng)力受地層巖性影響,單點地應(yīng)力測試結(jié)果往往不能代表研究區(qū)域整體地應(yīng)力狀態(tài)。因此,以實測地應(yīng)力數(shù)據(jù)點為基礎(chǔ),建立合理的地應(yīng)力測井預(yù)測模型獲取地應(yīng)力剖面,從點到線對地應(yīng)力的縱向分布特征進行研究,減小巖性對地應(yīng)力結(jié)果的影響,增大了地應(yīng)力場反演約束數(shù)據(jù)量。在此基礎(chǔ)上,利用地震構(gòu)造模型構(gòu)建合理的三維空間數(shù)值模型,以單井地應(yīng)力剖面計算結(jié)果為基本約束,利用有限元法進行地應(yīng)力場的數(shù)值模擬反演,實現(xiàn)從線到體對工區(qū)地應(yīng)力場進行研究,從而獲取研究工區(qū)地應(yīng)力的橫向分布特征和縱向分布特征。

1 單井地應(yīng)力測試分析

原始地應(yīng)力測試就是確定存在于研究位置點及其周圍區(qū)域的未受擾動的三維應(yīng)力狀態(tài),包括地應(yīng)力大小以及地應(yīng)力方向兩個方面。工區(qū)目的儲層巖性為砂礫巖,物性條件差,一般都采用水力壓裂來改造儲層提高產(chǎn)能以達到工業(yè)要求。因此,采用資料豐富的壓裂改造資料及地球物理測井資料分析深部地應(yīng)力。

1.1 地應(yīng)力大小

水力壓裂資料反演地應(yīng)力大小是深部地層應(yīng)力測量最有效的方法,也是國際巖石力學(xué)測試技術(shù)委員會推薦的巖體應(yīng)力測量的主要方法之一[11]?；趲r石為連續(xù)、均質(zhì)和各向同性的假設(shè),依據(jù)巖石力學(xué)分析和能量最低原則,水力壓裂縫的起裂發(fā)生在井壁切向應(yīng)力最小的部位,當(dāng)井軸與垂直主應(yīng)力方向一致時,裂縫的發(fā)育方位指示水平最大主應(yīng)力方位。壓裂過程中的裂縫閉合壓力是縫面相互接觸時的壓力,此時裂縫深度已達原巖應(yīng)力狀態(tài),故閉合壓力反映水平最小主應(yīng)力大小。從壓裂施工動態(tài)變化曲線中可以讀出破裂壓力和閉合壓力,根據(jù)式(1)可計算出水平最大主應(yīng)力,其中巖石抗張強度、biot系數(shù)及地層壓力均來自項目測試結(jié)果。

Pf=3σH2-σH1-αPP+St

(1)

式(1)中:Pf為破裂壓力,MPa；σH2為水平最大主應(yīng)力,MPa；σH1為水平最小主應(yīng)力,MPa；α為biot系數(shù)；PP為地層孔隙壓力,MPa；St為巖石抗張強度,MPa。

基于上述方法對研究區(qū)28個井段的壓裂資料進行分析,分別獲得壓裂深度點的破裂壓力和閉合壓力,計算出對應(yīng)的水平最大主應(yīng)力和水平最小主應(yīng)力。結(jié)果表明,水平最小主應(yīng)力分布在32.54～74.59 MPa,平均值為56.49 MPa；水平最大主應(yīng)力分布在37.48～99.91 MPa,平均值為70.17 MPa。

1.2 地應(yīng)力方位

在鉆井過程中,由于鉆具震動、原始地應(yīng)力釋放、鉆井液密度過大等原因會形成分布于井壁上的誘導(dǎo)裂縫,與地應(yīng)力有關(guān)的重泥漿壓裂縫和應(yīng)力釋放縫的方位都指示水平最大地應(yīng)力方向。根據(jù)彈性力學(xué)理論分析,高應(yīng)力作用下井壁附近在水平最小主應(yīng)力方向產(chǎn)生應(yīng)力集中,當(dāng)切向應(yīng)力超過井壁巖石抗剪強度時,井壁剪切破壞出現(xiàn)井眼應(yīng)力崩落現(xiàn)象,崩落方位指示水平最小主應(yīng)力方向。成像測井具有高密度采樣和高分辨率的特點,可以有效識別鉆井誘導(dǎo)縫和井眼崩落,獲取水平主應(yīng)力的方位[12-13]?；诘貙游㈦娮杪蕭呙璩上?formation micro-scanner image，F(xiàn)MI)測井資料對研究工區(qū)內(nèi)共6個井段的鉆井誘導(dǎo)縫和井眼崩落進行方位分析,某井部分井段分析結(jié)果如圖1所示。由圖1可見羽狀排列的應(yīng)力釋放縫組,指示該井水平最大主應(yīng)力方位為NWW向,優(yōu)勢方位約為105°。

圖1 誘導(dǎo)縫指示水平最大主應(yīng)力方位Fig.1 The direction of horizontal maximum principal stress indicated by conducted fracture

在各向異性地層中,橫波沿傳播方向?qū)⒎至殉蓛蓚€速度不同相互垂直的子波(快、慢橫波),這種現(xiàn)象叫橫波分裂。井周地層受不同方向的水平地應(yīng)力作用導(dǎo)致受壓程度不同,沿水平最大主應(yīng)力方向巖石受壓程度最高,快橫波的方位可以指示水平最大主應(yīng)力的方向[14]?；赟onic Scanner測井?dāng)?shù)據(jù),通過提取快、慢橫波方位信息,分析快、慢橫波時差大小及能量差等各向異性指標,對工區(qū)12個井段進行水平地應(yīng)力分析。圖2為某井部分井段分析結(jié)果,結(jié)果表明該井段快橫波指示的水平最大主應(yīng)力方位為E-W向,優(yōu)勢方位為90°。

DTEM_FAST為快橫波,μs/ft；DTSM_SLOW為慢橫波,μs/ft；FSH_P1N0為快慢橫波的方位,(°)圖2 快橫波指示最大主應(yīng)力方位Fig.2 The direction of horizontal maximum principal stress indicated by fast shear wave

分別采用FMI、Sonic Scanner兩種不同的成像測井資料,基于不同巖石物理力學(xué)理論分析得到不同井位的水平主應(yīng)力方位。結(jié)果表明,工區(qū)水平最大主應(yīng)力整體呈現(xiàn)E-W向。受局部構(gòu)造作用影響,不同井位處水平最大主應(yīng)力方位存在一定差異。

2 地應(yīng)力縱向剖面計算

以地應(yīng)力實測資料為基礎(chǔ),建立合理的地應(yīng)力計算模型,利用地球物理測井資料進行地應(yīng)力縱向剖面計算分析,是目前中外石油工業(yè)界研究地應(yīng)力的常用方法[15-16]。

2.1 垂向地應(yīng)力

地層垂向地應(yīng)力由上覆地層重力引起,數(shù)值大小取決于地層深度與巖石體積密度,采用密度測井的積分進行計算,公式為

(2)

式(2)中:σV為垂向地應(yīng)力,MPa；H0為測井起始點深度,m；ρ0(h)為未測井段深度為h點的密度,g/cm3；ρ(h)為深度為h點的測井密度,g/cm3；g為重力加速度,m/s2。

2.2 水平地應(yīng)力

利用測井資料計算水平地應(yīng)力依賴于合理的計算模型,目前常用的模型有Mohr-Columb模型、Matthews & Kelly模型、Anderson模型、Newberry模型、黃榮樽模型、組合彈簧模型等[17]。其中組合彈簧關(guān)系模型認為水平主應(yīng)力是上覆地層壓力和水平構(gòu)造應(yīng)力共同作用的結(jié)果,綜合考慮了地層巖石力學(xué)特性、孔隙壓力及構(gòu)造作用對地應(yīng)力的影響,在實際工程中應(yīng)用廣泛。模型計算公式為

(3)

式(3)中:μ為泊松比；E為巖石彈性模量,MPa；εH為沿水平最大主應(yīng)力方向構(gòu)造應(yīng)變系數(shù)；εh為沿水平最小主應(yīng)力方向構(gòu)造應(yīng)變系數(shù)。

根據(jù)地應(yīng)力計算模型確定地層巖石力學(xué)參數(shù)和地層孔隙壓力是獲取地應(yīng)力必需的基礎(chǔ)參數(shù)。巖石力學(xué)參數(shù)由巖石力學(xué)實驗數(shù)據(jù)建立模型,結(jié)合測井?dāng)?shù)據(jù)計算獲取；地層孔隙壓力獲取基于有效應(yīng)力理論,由實測資料及測井?dāng)?shù)據(jù)進行多元回歸分析建模獲得。構(gòu)造應(yīng)變系數(shù)εH、εh是開展地應(yīng)力剖面研究的關(guān)鍵性參數(shù),利用地應(yīng)力測試結(jié)果,根據(jù)式(3)反演即可得到。研究工區(qū)塊部分井構(gòu)造應(yīng)變系數(shù)計算結(jié)果如表1所示。

表1 部分井構(gòu)造應(yīng)變系數(shù)分析結(jié)果Table 1 Results of tectonic strain ratio in some wells

2.3 地應(yīng)力剖面

采用式(2)及式(3)的組合彈簧模型,結(jié)合分層構(gòu)造應(yīng)變系數(shù),利用全井巖石楊氏模量、泊松比、地層孔隙壓力等參數(shù),可計算得到全井的地應(yīng)力縱向剖面(圖3)。計算分析結(jié)果表明,研究區(qū)域垂向地應(yīng)力梯度大小為2.27～2.37 MPa/100 m,水平最大主應(yīng)力梯度大小為1.79～2.71 MPa/100 m,水平最小主應(yīng)力梯度大小為1.74～2.15 MPa/100 m,應(yīng)力狀態(tài)類型整體表現(xiàn)為潛在正斷型應(yīng)力場,即3個主應(yīng)力大小關(guān)系為σV>σH>σh。

GR為自然伽馬,API；CAL為井徑,in(1 in=25.4 mm)；AC為縱波時差,μs/ft；DEN為密度,g/cm3；CNL為中子,%；RT為深電阻率,Ω·m；RXO為淺電阻率,Ω·m；U為單軸抗壓強度,MPa；E為彈性模量,MPa；PP為地層孔隙壓力,MPa/100 m；SVC為垂向主應(yīng)力,MPa/100 m；SHMG為水平最大主應(yīng)力,MPa/100 m；SHNG為水平最小主應(yīng)力,MPa/100 m圖3 某井地應(yīng)力縱向剖面圖Fig.3 The profile of in-situ stress in a well

3 地應(yīng)力場數(shù)值模擬

3.1 數(shù)值模型

地應(yīng)力場與區(qū)域地質(zhì)構(gòu)造、地層巖石物理力學(xué)特性密切相關(guān),地質(zhì)模型的構(gòu)建和巖石力學(xué)參數(shù)的賦值是地應(yīng)力場反演的關(guān)鍵。由于地下深部構(gòu)造復(fù)雜多樣,考慮到計算工作量和計算精度,在深入認識研究區(qū)構(gòu)造格架的基礎(chǔ)上,對實際地質(zhì)模型進行合理的抽象簡化,基于已有的地球物理勘探資料地震解釋數(shù)據(jù)體,建立了圖4所示的三維數(shù)值模型。模型構(gòu)建過程中充分考慮目的層的展布特征和起伏特性,并對研究工區(qū)8條規(guī)模較大的斷裂破碎帶進行重點描述。模型范圍對應(yīng)的大地坐標為X(15 409 580～15 436 580)、Y(5 095 665～5 110 665)。模型面積為27×15=405 km2,采用八節(jié)點四面體單元進行網(wǎng)格離散劃分,經(jīng)離散后模型共有節(jié)點85 721個,單元696 886個。

圖4 研究工區(qū)目的層三維數(shù)值模型Fig.4 3D numerical model of the target formation

3.2 參數(shù)賦值

反演地應(yīng)力場采用的是線彈性本構(gòu)模型,需要對地質(zhì)體的巖石力學(xué)參數(shù)(楊氏模量、泊松比)和體積密度賦值。傳統(tǒng)數(shù)值模擬中力學(xué)參數(shù)的賦值往往將介質(zhì)假設(shè)為均值體,而這往往導(dǎo)致分析結(jié)果的誤差較大,故采用構(gòu)建測井預(yù)測模型獲取巖石力學(xué)參數(shù)的方法[18]。因此,通過室內(nèi)三軸壓縮試驗得到巖樣的靜態(tài)巖石力學(xué)參數(shù),通過縱、橫波測試實驗得到巖樣的聲波時差等物性參數(shù),再利用靜態(tài)力學(xué)參數(shù)與聲波時差等物性參數(shù)做數(shù)理統(tǒng)計法擬合獲得巖石力學(xué)參數(shù)測井預(yù)測模型,結(jié)合測井資料,獲得巖石力學(xué)參數(shù)的縱向剖面。在此基礎(chǔ)上,利用地震速度場三維數(shù)據(jù)體,使用井震聯(lián)合技術(shù),建立巖石力學(xué)參數(shù)的空間分布特性[19],為地應(yīng)力場反演提供巖石力學(xué)參數(shù)賦值基礎(chǔ)。由于研究工區(qū)內(nèi)斷裂帶巖石松散破碎,無法取樣測試,因此按照行業(yè)常用方法,將斷裂帶巖石力學(xué)參數(shù)按照一定比例降低[20],取用相應(yīng)空間位置點巖石的2/3。

3.3 數(shù)值模擬反演分析

有限元模型建立后,以實測地應(yīng)力方向為指導(dǎo),對模型東側(cè)和北側(cè)設(shè)置模型邊界荷載作用,并在模型西側(cè)、南側(cè)和底面設(shè)置法向位移約束以防止剛性漂移。合理確定施加邊界后,以工區(qū)地應(yīng)力測試及測井計算地應(yīng)力大小和方向作為約束和修正依據(jù),采用線性疊加原理反復(fù)修改邊界條件直到誤差最小,對地應(yīng)力場進行多次正演和反演試算,最終確定合理的應(yīng)力加載方式,并得出研究工區(qū)現(xiàn)今應(yīng)力場的真實情況。

圖5(a)～圖5(c)分別為數(shù)值模擬計算得到的研究工區(qū)砂礫巖目的層的水平最大主應(yīng)力、水平最小主應(yīng)力以及垂向主應(yīng)力的三維空間分布圖。由圖5可見,研究工區(qū)整體處于壓應(yīng)力作用下,垂向主應(yīng)力分布在35～90 MPa范圍內(nèi)，水平最大主應(yīng)力分布在40～80 MPa范圍內(nèi)，水平最小主應(yīng)力分布在40～70 MPa范圍內(nèi)。受地層埋深和空間起伏影響,地應(yīng)力整體呈現(xiàn)由西北向東南逐漸增大的趨勢,地應(yīng)力狀態(tài)從局部走滑型(σH>σV>σh)逐漸轉(zhuǎn)變過渡為整體正斷型(σV>σH>σh)。工區(qū)斷裂破碎帶3個主應(yīng)力較低,而斷裂帶附近垂向主應(yīng)力、水平最小主應(yīng)力變化較大且存在局部高應(yīng)力帶,此特征在C5、C6、C1井旁的北東向斷裂尤為顯著,這與斷裂帶空間位置和形態(tài)有關(guān)。圖6所示為數(shù)值模擬計算得到的工區(qū)地層水平最大主應(yīng)力方位,可見整體水平最大主應(yīng)力方位呈近E-W,范圍為75°～95°。在斷裂附近,地應(yīng)力向斷裂延伸方向偏轉(zhuǎn)的趨勢明顯,表現(xiàn)為NE、NWW走向斷裂附近地應(yīng)力方位呈NE向，NNW走向斷裂附近地應(yīng)力方位趨于垂直斷裂延伸方向。

圖5 三維空間分布Fig.5 3D space distribution

圖6 水平最大主應(yīng)力方位分布Fig.6 The azimuth of horizontal maximum principal stress

4 地應(yīng)力場合理性驗證

為驗證地應(yīng)力場模擬結(jié)果的合理性,使用研究區(qū)未參與約束和修正的5口井單點地應(yīng)力測試資料與數(shù)值模擬計算結(jié)果進行對比,結(jié)果如表2所示。由表2可見二者的相對誤差較小,說明數(shù)值模擬結(jié)果可以較好地反映本區(qū)塊現(xiàn)今地應(yīng)力場的三維展布特征。

表2 單井地應(yīng)力測試結(jié)果與數(shù)值模擬結(jié)果對比Table 2 Comparison between experimental measurement results and numerical simulation results

研究區(qū)塊構(gòu)造格局形成于白堊紀早期。從區(qū)域構(gòu)造演化歷史來看,海西構(gòu)造運動和印支構(gòu)造運動對于研究區(qū)構(gòu)造格局的形成起了至關(guān)重要的作用[21]。海西運動、印支運動造成西北緣地殼局部隆起,燕山期后沖斷帶格局基本形成,喜馬拉雅運動后成為碰撞前陸型沉積坳陷。故三疊紀的印支運動產(chǎn)生S-N向的擠壓應(yīng)力,是研究區(qū)目的層現(xiàn)今水平主應(yīng)力E-W方向形成的主要原因。工區(qū)主要發(fā)育兩條NE向斷裂,是次級斷裂,與主斷裂呈平行狀態(tài),其主要形成于海西構(gòu)造運動時期,后在印支運動形成現(xiàn)今應(yīng)力場的作用下,導(dǎo)致斷裂附近地應(yīng)力場方向沿著斷裂偏轉(zhuǎn)。因此,研究獲取的工區(qū)地應(yīng)力方向與區(qū)域構(gòu)造背景狀態(tài)相符。

5 結(jié)論

(1)基于豐富的壓裂改造和地球物理測井資料,分析了研究工區(qū)已鉆井目的層的實測地應(yīng)力大小與方向,建立了適用于砂礫巖儲層的地應(yīng)力評價模型并計算了地應(yīng)力單井縱向剖面。

(2)在地應(yīng)力縱向剖面研究基礎(chǔ)上,利用地震構(gòu)造模型構(gòu)建合理的三維空間數(shù)值模型,以單井地應(yīng)力為基本約束,利用有限元法對工區(qū)現(xiàn)今三維地應(yīng)力場進行反演。研究表明,工區(qū)地應(yīng)力整體處于壓應(yīng)力作用下,表現(xiàn)為西北向東南方向逐漸增大的趨勢。受斷裂發(fā)育及次級構(gòu)造影響,斷裂內(nèi)3個主應(yīng)力顯著降低,斷裂邊緣和端部存在應(yīng)力激增帶。地應(yīng)力方位整體近E-W,受斷裂和局部構(gòu)造發(fā)育的影響,局部地應(yīng)力方向發(fā)生偏轉(zhuǎn)。