支持向量機(jī)在煤層地應(yīng)力預(yù)測(cè)中的應(yīng)用

馮 鵬，李 松，湯達(dá)禎，陳 博，鐘廣浩

(1. 中國地質(zhì)大學(xué)(北京) 能源學(xué)院，北京 100083；2. 煤層氣開發(fā)利用國家工程研究中心中國地質(zhì)大學(xué)(北京) 煤儲(chǔ)層物性實(shí)驗(yàn)室，北京 100083)

0 引 言

地應(yīng)力場(chǎng)的研究對(duì)煤層氣開發(fā)意義重大[1-2]。由于煤層本身的低孔低滲特征，決定了其只有經(jīng)過大規(guī)模壓裂改造才能獲得商業(yè)產(chǎn)能[3]，而地應(yīng)力是影響壓裂改造效果的關(guān)鍵因素，煤儲(chǔ)層中割理(裂隙)系統(tǒng)的開啟程度和增產(chǎn)改造壓裂縫的形態(tài)擴(kuò)展均受現(xiàn)今地應(yīng)力場(chǎng)的控制[4-7]，針對(duì)煤層開展地應(yīng)力研究是進(jìn)行煤儲(chǔ)層開采的必要環(huán)節(jié)。近年來，我國煤層氣勘探開發(fā)在地質(zhì)理論和開發(fā)技術(shù)方面已取得重大進(jìn)展，但在煤層地應(yīng)力方面的研究相對(duì)薄弱[8]。前人經(jīng)過大量研究，提出了數(shù)十種地應(yīng)力測(cè)量方法[9]。利用儀器測(cè)量得到的地應(yīng)力比較準(zhǔn)確，但受工程、技術(shù)及成本等因素的限制，特別是對(duì)深部地層而言，不可能進(jìn)行大規(guī)模的地應(yīng)力測(cè)量。為了降低成本，許多學(xué)者提出了地應(yīng)力計(jì)算模型[10-14]，如Mattews和Kelly模型、Anderson模型、Mohr模型、Coulomb模型、黃氏經(jīng)驗(yàn)?zāi)Ｐ汀⒏鹗辖?jīng)驗(yàn)?zāi)Ｐ秃徒M合彈簧模型等。還有學(xué)者結(jié)合技術(shù)手段對(duì)地應(yīng)力進(jìn)行預(yù)測(cè)：如黃繼新等建立了應(yīng)用成像測(cè)井資料進(jìn)行地應(yīng)力分析的方法[15]；劉景武、傅海成等基于多極子聲波測(cè)井探測(cè)巖層裂縫、判斷地應(yīng)力方位等進(jìn)而獲得巖石力學(xué)性質(zhì)參數(shù)[16-20]；樓一珊、林英松等首先利用聲波測(cè)井資料計(jì)算出巖體完整系數(shù)，進(jìn)一步基于靜態(tài)力學(xué)參數(shù)計(jì)算了地應(yīng)力的大小[21-22]；王曉鋒、陳書平等利用有限元方法建立了煤儲(chǔ)層的地質(zhì)模型以及模擬了各地質(zhì)時(shí)期的地應(yīng)力分布特征[23-24]；魏忠文等利用構(gòu)造應(yīng)力場(chǎng)數(shù)值模擬方法計(jì)算了構(gòu)造應(yīng)力場(chǎng)的分布特征[25]；趙冰等將有限元線法同巖土工程反分析理論相結(jié)合，探討了用巖土工程反演確定初始地應(yīng)力和圍巖參數(shù)方法的可行性[26]；張有天等根據(jù)少數(shù)測(cè)點(diǎn)數(shù)據(jù)及地表邊界條件對(duì)地應(yīng)力場(chǎng)進(jìn)行趨勢(shì)分析，提出了二維及三維地應(yīng)力場(chǎng)的計(jì)算方法[27]；劉志強(qiáng)基于三軸壓縮和凱撒實(shí)驗(yàn)，利用水力壓裂和測(cè)井?dāng)?shù)據(jù)，建立了一套煤層地應(yīng)力預(yù)測(cè)方法[28]。

雖然前人提出的地應(yīng)力預(yù)測(cè)模型具有顯著成效并得到了較好的應(yīng)用，但對(duì)于煤層地應(yīng)力的預(yù)測(cè)仍然存在一些不足。一方面，大多數(shù)地應(yīng)力模型都是經(jīng)驗(yàn)性的，存在一定的地域局限性；模型中還有大量參數(shù)需要利用其他方法或公式進(jìn)行計(jì)算，過程繁瑣，且參數(shù)的確定進(jìn)一步增大了誤差，實(shí)用性受到限制。另一方面，大多數(shù)的地應(yīng)力預(yù)測(cè)研究都僅針對(duì)常規(guī)砂巖儲(chǔ)層，而煤巖與砂巖存在顯著的性質(zhì)差異，導(dǎo)致上述方法不能直接應(yīng)用在煤層上[28]。因此，建立一個(gè)簡便且適用于煤層的地應(yīng)力預(yù)測(cè)模型十分必要。由于測(cè)井資料具有測(cè)量深度較深、信息量大、數(shù)據(jù)相對(duì)連續(xù)等特點(diǎn)，基于測(cè)井資料計(jì)算地應(yīng)力具有獨(dú)特的優(yōu)勢(shì)[9]。近年來在地球物理反演中應(yīng)用較廣泛的支持向量機(jī)算法是一種可用于對(duì)數(shù)據(jù)進(jìn)行高維特征分析的重要算法[29]。因此，本文基于測(cè)井資料，嘗試將支持向量機(jī)算法應(yīng)用于煤層地應(yīng)力的預(yù)測(cè)中，并進(jìn)一步借助三維地質(zhì)建模軟件完成地應(yīng)力分布模型的建立，為解決煤層地應(yīng)力的預(yù)測(cè)難題探索新的思路和方法。

1 實(shí)測(cè)地應(yīng)力計(jì)算方法

1.1 水平主應(yīng)力計(jì)算方法

根據(jù)地面垂直鉆孔水力壓裂測(cè)量地應(yīng)力方法[30]，閉合壓力Pc為最小水平主應(yīng)力σh，即

σh=Pc

(1)

最大水平主應(yīng)力σH的計(jì)算公式為：

σH=3Pc-Pf-Po+T

(2)

式中：T為煤或巖石的抗拉強(qiáng)度，MPa，本次計(jì)算中取定值0.8；Pf為煤層破裂壓力，MPa；Po為儲(chǔ)層壓力，MPa。

基于鄂爾多斯盆地典型煤層氣區(qū)塊注入/壓降試井?dāng)?shù)據(jù)，利用以上公式求出煤層的最大水平主應(yīng)力σH和最小水平主應(yīng)力σh，并對(duì)二者進(jìn)行了相關(guān)性分析(圖1)。最大與最小水平主應(yīng)力之間呈正相關(guān)關(guān)系，且相關(guān)性較強(qiáng)(R2=0.93)，故在已知最小水平主應(yīng)力σh的情況下，可利用如下模型計(jì)算最大水平主應(yīng)力σH：

σH=1.7394σh-3.4279

(3)

圖1 最大與最小水平主應(yīng)力擬合關(guān)系圖Fig.1 Fitting relation diagram of maximum and minimum horizontal principal stress

1.2 垂向主應(yīng)力計(jì)算方法

垂向主應(yīng)力σv是由上覆巖石的自重引起的，所以計(jì)算公式為：

(4)

其中：ρ為地層密度，kg/m3；g為重力加速度，取9.8 m/s2；h為深度，m。

實(shí)際計(jì)算時(shí)簡化積分，可以用如下公式計(jì)算[31]:

σv=0.027×h

(5)

2 最小水平主應(yīng)力預(yù)測(cè)方法

2.1 基于SVM回歸的煤層地應(yīng)力預(yù)測(cè)原理

支持向量機(jī)(SVM)是Cortes和Vapnik于1995年首先提出的[32]，它在解決小樣本、非線性及高維模式識(shí)別中有獨(dú)特優(yōu)勢(shì)[33]。研究所用的支持向量回歸模型，其基本思想包括升維和線性化。即當(dāng)數(shù)據(jù)集Xi在原始特征中不是線性可分時(shí)，通過非線性函數(shù)(Xi)把原始特征空間映射到一個(gè)高維特征空間，SVM算法將在這個(gè)高維特征空間中找到一個(gè)函數(shù)間隔最大的線性超平面來對(duì)數(shù)據(jù)進(jìn)行分類。最后把高維的分類面映射回低維，就可以在低維空間對(duì)數(shù)據(jù)分類(圖2)。

圖2 SVM結(jié)構(gòu)示意圖Fig.2 Structure schematic diagram of the Support Vector Machine (SVM)

在高維特征空間中構(gòu)造最優(yōu)超平面只需要計(jì)算支持向量與特征空間中向量的內(nèi)積，而SVM不能直接計(jì)算該內(nèi)積，因此需要用滿足Mercer定理的核函數(shù)K(x,y)來代替[34]。本文選取徑向基函數(shù)作為核函數(shù)：

(6)

圖3 最小水平主應(yīng)力預(yù)測(cè)流程圖Fig.3 Prediction flow chart of minimum horizontal principal stress

2.2 煤層地應(yīng)力主要關(guān)聯(lián)因子的篩選

測(cè)井?dāng)?shù)據(jù)類型多樣，為了避免引入某些噪聲信息，在建立預(yù)測(cè)模型之前，需要對(duì)輸入的測(cè)井?dāng)?shù)據(jù)類型進(jìn)行篩選?；疑P(guān)聯(lián)分析方法是一種多因素統(tǒng)計(jì)分析方法，以各因素的樣本數(shù)據(jù)為依據(jù)，用灰色關(guān)聯(lián)度描述因素間關(guān)系的強(qiáng)弱、大小和次序[35]。前人提出了多種關(guān)聯(lián)度模型，本文采用鄧氏關(guān)聯(lián)度模型進(jìn)行分析，Xi與X0的關(guān)聯(lián)度可表示為：

(7)

其中：

(8)

通過對(duì)鄂爾多斯盆地東緣典型煤層氣區(qū)塊的測(cè)井和注入/壓降試井?dāng)?shù)據(jù)的整理，將包含9種不同類型測(cè)井?dāng)?shù)據(jù)(密度、自然伽馬、中子、自然電位、聲波時(shí)差、兩個(gè)方向井徑、深側(cè)向和淺側(cè)向電阻率)及對(duì)應(yīng)的最小水平主應(yīng)力的34個(gè)樣本數(shù)據(jù)(表1)導(dǎo)入Matlab軟件中，作為34×10的矩陣進(jìn)行分析。由于各組數(shù)據(jù)的量綱不同，計(jì)算關(guān)聯(lián)系數(shù)之前，先將數(shù)列進(jìn)行初值化處理。通過公式(7)計(jì)算得到如下結(jié)果(圖4)：聲波時(shí)差(DT)和自然電位(SP)測(cè)井值與最小水平主應(yīng)力關(guān)聯(lián)度相對(duì)較低，中子測(cè)井(CNL)與最小水平主應(yīng)力關(guān)聯(lián)度最高。以關(guān)聯(lián)度大于0.7為界，可得出影響最小水平主應(yīng)力的主要因子為兩個(gè)方向的井徑測(cè)井(CAL1和CAL2)、深淺側(cè)向電阻率測(cè)井(LLD和LLS)、補(bǔ)償中子測(cè)井(CNL)、自然伽馬測(cè)井(GR)和密度測(cè)井(DEN)共7個(gè)因素。由于CAL1和CAL2以及LLD和LLS的值基本一致，分別與最小水平主應(yīng)力的關(guān)聯(lián)度也基本相等，為了進(jìn)一步減少影響因子的數(shù)量，以CAL1和CAL2的平均值(CAL)作為一個(gè)影響因子，以LLD和LLS的平均值(R)作為一個(gè)影響因子。由此得到最小水平主應(yīng)力的5個(gè)主要影響因子：CAL、R、CNL、GR和DEN。

圖4 各類測(cè)井?dāng)?shù)據(jù)的關(guān)聯(lián)度柱狀圖Fig.4 Correlation histogram of various types of logging data

2.3 基于SVM的回歸模型的建立

按照預(yù)測(cè)樣本和訓(xùn)練樣本1∶3的比例從34個(gè)樣本中隨機(jī)選取了9個(gè)預(yù)測(cè)樣本，其余25個(gè)樣本用于訓(xùn)練(表1)。其中井徑測(cè)井(CAL)、補(bǔ)償中子測(cè)井(CNL)、自然伽馬測(cè)井(GR)、密度測(cè)井(DEN)和深淺側(cè)向電阻率測(cè)井值的平均值(R)作為輸入數(shù)據(jù)，最小水平主應(yīng)力作為輸出數(shù)據(jù)。

在對(duì)所有數(shù)據(jù)利用離差標(biāo)準(zhǔn)化方法進(jìn)行歸一化處理之后，將輸入數(shù)據(jù)和學(xué)習(xí)目標(biāo)導(dǎo)入編寫好的程序中，對(duì)樣本進(jìn)行多次迭代訓(xùn)練。通過不斷調(diào)整，最終得出參數(shù)c和μ的值分別為20.06和1.94時(shí)模型最優(yōu)。對(duì)模型進(jìn)行檢驗(yàn)，實(shí)測(cè)值(由式(1)計(jì)算所得)與預(yù)測(cè)值平方相關(guān)系數(shù)可達(dá)到0.94，平均相對(duì)誤差為0.07，具有較好的預(yù)測(cè)能力(圖5)。

3 實(shí)例應(yīng)用

韓城區(qū)塊位于鄂爾多斯盆地東緣南部，渭北煤田東北緣，整體表現(xiàn)為北東走向、北北西傾向的單斜構(gòu)造。位于韓城區(qū)塊的H3井組煤層埋深跨度較大，受構(gòu)造活動(dòng)影響較大，導(dǎo)致煤層的地應(yīng)力場(chǎng)在平面和垂向上均存在較大差異，而地應(yīng)力是煤層氣開采的一個(gè)重要的影響因素[36]。因此，開展該區(qū)H3井組11號(hào)煤層地應(yīng)力的預(yù)測(cè)，同時(shí)借助Petrel 2015軟件完成該區(qū)地應(yīng)力場(chǎng)的三維地質(zhì)建模，厘定研究區(qū)煤層地應(yīng)力在平面和垂向上的展布特征，為該區(qū)煤層氣勘探開發(fā)提供理論依據(jù)。

表1 訓(xùn)練樣本和測(cè)試樣本原始數(shù)據(jù)

圖5 最小水平主應(yīng)力實(shí)測(cè)值與預(yù)測(cè)值對(duì)比圖Fig.5 Comparison between actual and predicted values of the minimum horizontal principal stress

3.1 研究區(qū)地應(yīng)力分布特征

三維地質(zhì)建模是用模型來表征儲(chǔ)層和物性參數(shù)的三維空間分布以及變化規(guī)律的一種手段[37-39]。本次研究基于SVM回歸模型計(jì)算韓城區(qū)塊H3井組11號(hào)煤層最小水平主應(yīng)力，進(jìn)而利用最小水平主應(yīng)力和煤層埋深分別計(jì)算得到最大水平主應(yīng)力和垂向應(yīng)力。以計(jì)算得到的地應(yīng)力數(shù)據(jù)作為基礎(chǔ)數(shù)據(jù)，運(yùn)用序貫高斯模擬法構(gòu)建煤層地應(yīng)力場(chǎng)的地質(zhì)模型，厘定煤層地應(yīng)力場(chǎng)的三維展布特征。計(jì)算表明，研究區(qū)11號(hào)煤層最小水平主應(yīng)力為4.86～27.12 MPa，最大水平主應(yīng)力為5.02～43.74 MPa，垂向主應(yīng)力15.93～24.96 MPa。陳世達(dá)等通過232組實(shí)測(cè)地應(yīng)力數(shù)據(jù)統(tǒng)計(jì)分析得出[2]，我國煤儲(chǔ)層最小水平主應(yīng)力為2.14～30.26 MPa，最大水平主應(yīng)力為2.70～45.20 MPa，垂向主應(yīng)力3.67～37.99 MPa。由此進(jìn)一步說明，本文通過支持向量機(jī)回歸建立的最小水平主應(yīng)力預(yù)測(cè)模型，反演得到的各應(yīng)力大小均在合理范圍。

地應(yīng)力場(chǎng)三維地質(zhì)建模結(jié)果顯示，從平面展布來看，研究區(qū)水平最小主應(yīng)力的高應(yīng)力范圍較小，且主要集中在研究區(qū)東部，西部次之，研究區(qū)中部和北部為水平主應(yīng)力的低應(yīng)力區(qū)；從垂向分布來看，該煤層頂?shù)装甯浇淖钚∷街鲬?yīng)力整體偏大，而煤層中部最小水平主應(yīng)力較小(圖6(a))。研究區(qū)最大水平主應(yīng)力平面展布和垂向分布規(guī)律與最小水平主應(yīng)力基本一致(圖6(b))。垂向應(yīng)力在平面展布上分布較為均勻，高應(yīng)力區(qū)以西南部為主，低應(yīng)力區(qū)以東北部為主，這是由煤層?xùn)|北高、西南低的地勢(shì)以及區(qū)域地質(zhì)構(gòu)造所造成；垂向上由于煤層較薄，沒有明顯變化規(guī)律(圖6(c))。

圖6 韓城礦區(qū)H3井組11號(hào)煤層地應(yīng)力分布柵狀圖(地應(yīng)力單位：MPa)Fig.6 Grid diagram of in-situ stress distribution in No.11 coal seam of well group H3, Hancheng mining area

3.2 不同深度的煤層地應(yīng)力分析

為了分析研究區(qū)不同深度煤層的應(yīng)力場(chǎng)分布狀態(tài)，基于前文建立的地應(yīng)力預(yù)測(cè)模型，反演獲得研究區(qū)295個(gè)地應(yīng)力數(shù)據(jù)點(diǎn)進(jìn)行應(yīng)力與埋深相關(guān)性分析(圖7)。

圖7 韓城礦區(qū)H3井組煤層地應(yīng)力與埋深擬合關(guān)系Fig.7 Fitting relations between in-situ coal seam stress and burial depth of H3 well group in Hancheng mining area

圖8 韓城礦區(qū)H3井組煤層側(cè)壓系數(shù)與埋深擬合關(guān)系Fig.8 Fitting relations between coal seam lateral pressure coefficient and burial depth of H3 well group in Hancheng mining area

由地應(yīng)力隨埋深的變化可以看出，該礦區(qū)H3井組煤層埋深為500～950 m，煤層所受三個(gè)方向的主應(yīng)力都隨著埋深的增加而增大，而垂向應(yīng)力隨埋深的增加而增加的幅度明顯大于水平主應(yīng)力，三個(gè)方向的應(yīng)力之間的相互關(guān)系隨著埋深不同也在發(fā)生變化，即在垂向上地應(yīng)力場(chǎng)的類型也存在分異(圖8)：

(1)500～600 m，以σH>σv>σh型應(yīng)力場(chǎng)為主，最小水平主應(yīng)力介于9.32～19.28 MPa之間，平均為13.10 MPa，最大水平主應(yīng)力占主導(dǎo)作用，地應(yīng)力表現(xiàn)為大地動(dòng)力型。水平方向與垂直方向的應(yīng)力大小關(guān)系可進(jìn)一步用側(cè)壓系數(shù)來表示，側(cè)壓系數(shù)即最大與最小水平主應(yīng)力的平均值與垂直主應(yīng)力的比值[4]。該深度范圍內(nèi)側(cè)壓系數(shù)較大，煤儲(chǔ)層主要處于擠壓帶，部分處于過渡帶。

(2)600～750 m，以σH≈σv>σh型應(yīng)力場(chǎng)為主，為地應(yīng)力場(chǎng)轉(zhuǎn)化的過渡階段，各地應(yīng)力間大小關(guān)系變化較為頻繁。在此深度范圍內(nèi)，大地動(dòng)力場(chǎng)型和大地靜力場(chǎng)型兩種應(yīng)力場(chǎng)共存。該深度范圍煤儲(chǔ)層主要處于過渡帶，部分處于伸張帶。

(3)750～950 m，以σv>σH>σh型應(yīng)力場(chǎng)為主，最小水平主應(yīng)力介于7.05～22.76 MPa之間，平均為13.85 MPa，垂向主應(yīng)力占主導(dǎo)作用，地應(yīng)力表現(xiàn)為大地靜力型。該深度范圍煤儲(chǔ)層主要處于伸張帶，部分處于過渡帶。

4 結(jié) 論

(1)基于灰色關(guān)聯(lián)法模型分析，井徑測(cè)井(CAL)、補(bǔ)償中子測(cè)井(CNL)、自然伽馬測(cè)井(GR)、密度測(cè)井(DEN)和深淺側(cè)向電阻率測(cè)井值的平均值(R)與煤層最小水平主應(yīng)力關(guān)聯(lián)度最好。以這5個(gè)測(cè)井參數(shù)作為訓(xùn)練因子，利用基于小樣本理論的支持向量機(jī)回歸方法建立煤層最小水平主應(yīng)力預(yù)測(cè)模型，實(shí)測(cè)值與預(yù)測(cè)值平方相關(guān)系數(shù)可達(dá)到0.94，因此該模型可以較為高效、精準(zhǔn)地對(duì)最小水平主應(yīng)力進(jìn)行預(yù)測(cè)。

(2)以鄂爾多斯盆地東緣韓城區(qū)塊H3井組煤層為研究對(duì)象，通過建立煤層地應(yīng)力場(chǎng)三維地質(zhì)模型以及地應(yīng)力與埋深的擬合關(guān)系綜合分析，研究區(qū)煤層埋深范圍在500～950 m，地應(yīng)力隨埋深的增加呈現(xiàn)遞增的趨勢(shì)。煤層埋深處于500～600 m時(shí)，以σH>σv>σh的大地動(dòng)力型應(yīng)力場(chǎng)為主；埋深處于600～750 m時(shí)，以σH≈σv>σh型應(yīng)力場(chǎng)為主，為地應(yīng)力場(chǎng)轉(zhuǎn)化的過渡階段；埋深處于750～950 m時(shí)，以σv>σH>σh的大地靜力型應(yīng)力場(chǎng)為主?？傮w而言，隨著埋深的增加，煤儲(chǔ)層所處的應(yīng)力環(huán)境由擠壓帶過渡為伸張帶。