基于改進BP神經(jīng)網(wǎng)絡(luò)的頁巖地應(yīng)力預(yù)測模型

尚福華，王瑋卿，曹茂俊

(東北石油大學(xué) 計算機與信息技術(shù)學(xué)院，黑龍江 大慶 163318)

0 引 言

頁巖氣是中國能源勘探的重要領(lǐng)域和目標(biāo)。其中，地應(yīng)力分析技術(shù)是指導(dǎo)頁巖壓裂施工的核心技術(shù)[1]。目前，初始地應(yīng)力數(shù)據(jù)主要靠現(xiàn)場實測獲得，但是通常需要在現(xiàn)場利用水力壓裂法，應(yīng)力-應(yīng)變恢復(fù)法等方法進行測量，這些測量方法測試周期較長，成本昂貴，測試條件有限，所以目前地應(yīng)力場一般的演算方法為以實測地應(yīng)力數(shù)據(jù)為基礎(chǔ)，通過數(shù)值模擬或數(shù)值分析來求解[2]。然而，由于地層在漫長的時期中經(jīng)歷多次地質(zhì)作用及構(gòu)造演化，地質(zhì)構(gòu)造類型、作用程度和方向不同，導(dǎo)致構(gòu)造應(yīng)力場的大小和方向存在時間和空間的變異性，只能通過它的外部相關(guān)因素的狀態(tài)和方式來把握該事物的信息[3]。

神經(jīng)網(wǎng)絡(luò)具有自學(xué)習(xí)、自組織、可以高速尋找最優(yōu)解等優(yōu)點，可成功擬合多數(shù)非線性連續(xù)函數(shù)。因此，針對建立工程區(qū)域內(nèi)相關(guān)的參數(shù)與巖體初始地應(yīng)力非線性映射關(guān)系問題，采用神經(jīng)網(wǎng)絡(luò)是一種有效的解決手段[4]。目前，已有一些學(xué)者應(yīng)用BP神經(jīng)網(wǎng)絡(luò)分別對三維應(yīng)力場和巖石力學(xué)參數(shù)進行了反演分析。侯連浪等[5]建立并訓(xùn)練了BP神經(jīng)網(wǎng)絡(luò)對頁巖靜彈性模量進行預(yù)測，楊志強等[6]利用金川礦區(qū)實測地應(yīng)力樣本建立了人工神經(jīng)網(wǎng)絡(luò)模型進行訓(xùn)練，獲得了金川礦區(qū)地應(yīng)力分布規(guī)律。然而，由于BP神經(jīng)網(wǎng)絡(luò)的優(yōu)化算法為梯度下降法，目標(biāo)函數(shù)通常又極其復(fù)雜，因此，存在著網(wǎng)絡(luò)訓(xùn)練易過擬合、易陷入局部最優(yōu)、收斂速度較慢以及隱層節(jié)點數(shù)難于確定等缺點[7]。

目前，一些傳統(tǒng)的自適應(yīng)梯度下降算法可以通過自適應(yīng)地為各個參數(shù)分配不同學(xué)習(xí)率來解決該問題，但同樣存在隨著迭代次數(shù)增加，學(xué)習(xí)率越來越慢而導(dǎo)致提前收斂的問題[8]。因此，針對頁巖氣地應(yīng)力預(yù)測模型的構(gòu)建問題，該文以頁巖橫觀各向同性為基礎(chǔ)，針對傳統(tǒng)的自適應(yīng)梯度下降算法使用改進的梯度下降算法，并將改進后的算法應(yīng)用于神經(jīng)網(wǎng)絡(luò)中，以長寧工區(qū)若干口井的聲波測井資料與對應(yīng)深度的實測地應(yīng)力值數(shù)據(jù)為依據(jù)，反演出目標(biāo)層段的應(yīng)力邊界條件，最后代入到數(shù)值模擬模型中進行正演計算，進行該工區(qū)初始地應(yīng)力預(yù)測。

1 頁巖地應(yīng)力物理模型

1.1 橫觀各向同性物理模型

近年來，地應(yīng)力計算模型的建立大多數(shù)針對各向同性儲層，其中黃榮樽模型[9]的應(yīng)用比較廣泛：

(1)

(2)

由于頁巖物質(zhì)組成不均勻，節(jié)理發(fā)育不均勻，因此平行結(jié)構(gòu)面方向和垂直結(jié)構(gòu)面方向物理力學(xué)性質(zhì)不同，呈現(xiàn)出橫觀各向同性的特征。鑒于此特性，計算時仍采用上述各向同性模型顯然不夠準(zhǔn)確。目前，研究中一般采用Thiercelin模型[10]，具體如下：

(3)

(4)

式中，σv為垂直地應(yīng)力；σH為水平最大主應(yīng)力；σh為水平最小主應(yīng)力；Pp為地層孔隙壓力；α為Biot系數(shù)；ξ為滲流系數(shù)；Eh為層理面橫向的楊氏模量；EV為層理面法向的楊氏模量；vh為層理面橫向的泊松比；vv為層理面法向的泊松比；εH為水平最大構(gòu)造應(yīng)變；εh為水平最小構(gòu)造應(yīng)變。其中垂向主應(yīng)力的數(shù)值近似等于上覆巖層壓力σv，可以看作為地層密度和深度的函數(shù)，因此可用如下公式求出垂向地應(yīng)力：

(5)

式中，h為地層埋藏深度(單位m)；ρ(h)為地層密度隨地層深度變化的函數(shù)；g為重力加速度，9.8。

實際的地層密度ρ(h)隨深度的變化規(guī)律復(fù)雜，用連續(xù)函數(shù)不易表示，所以在實際計算中可以用分段加和的方法來計算：

(6)

式中，ρi為第i段地層平均體積密度(單位kg/m3)；ΔDi為第i段地層厚度(單位m)。

上覆巖層壓力計算完畢后，地層巖石的彈性參數(shù)可由鉆井取芯進行室內(nèi)試驗獲得，但由于鉆井取芯的不連續(xù)性及高成本，可以通過聲波測井資料、密度測井資料和自然伽馬測井資料來計算獲得。

1.2 改進的頁巖地應(yīng)力計算模型

由于大多數(shù)研究區(qū)構(gòu)造復(fù)雜，頁巖具有明顯的層理結(jié)構(gòu)，受層理分布的影響，地層自上而下傾角與傾向均有明顯變化，而地層傾角使得地層的上覆巖層壓力降低，增大了水平應(yīng)力，使水平地應(yīng)力的計算結(jié)果有一定誤差。并且，由重力、構(gòu)造應(yīng)力共同作用下的巖石應(yīng)變值在實際工程中不易獲得。為了提高地應(yīng)力計算精度，該文在Thiercelin模型的基礎(chǔ)上，考慮地層傾角因素，對上述地應(yīng)力計算模型進行改進：

(σv-αPp)sinφsin(β-β0)+αPp

(7)

(σv-αPp)sinφsin(β-β0)+αPp

(8)

式中，A為最大水平主應(yīng)力方向構(gòu)造運動(位移)系數(shù)；B為最小水平主應(yīng)力方向構(gòu)造運動(位移)系數(shù)；φ為地層傾角(單位°)；β為方位角(單位°)；β0為最大水平主應(yīng)力方向角(單位°)；K為側(cè)壓系數(shù)。

模型基于空間三角函數(shù)應(yīng)力分布關(guān)系進行計算，其中第一項首先通過構(gòu)造應(yīng)力系數(shù)反映水平構(gòu)造劇烈程度，其次反映由于傾角作用使得上覆巖層壓力對地層的壓實程度的減少情況；第二項考慮井斜方位角與構(gòu)造運動方向不同而產(chǎn)生的剪切應(yīng)力；最后，公式采用上覆巖層壓力產(chǎn)生的水平構(gòu)造應(yīng)力與原始水平構(gòu)造應(yīng)力加和的方式，代替利用地層的彈性應(yīng)變求彈性應(yīng)力的方式，避免使用因深度過大，取芯困難從而難以獲得的應(yīng)變數(shù)據(jù)。當(dāng)?shù)貙觾A角為0時，對地應(yīng)力的影響也為0，此時模型還原為黃氏模型。

2 BP神經(jīng)網(wǎng)絡(luò)模型

在上述分析模型中，巖層巖石力學(xué)參數(shù)，包括各方向的楊氏模量、泊松比，井斜方位角等參數(shù)雖可用測井曲線資料計算獲得，但由于巖石在不同方向所受的構(gòu)造運動，溫度影響不同，使得模擬地層構(gòu)造環(huán)境過程復(fù)雜，不同水平方向的構(gòu)造應(yīng)力往往不同。因此模型中構(gòu)造應(yīng)力系數(shù)只能通過該地區(qū)的經(jīng)驗系數(shù)獲得，結(jié)果誤差較大。為了解決地層構(gòu)造應(yīng)力系數(shù)獲得較難的問題，該文利用BP神經(jīng)網(wǎng)絡(luò)的映射功能反演某區(qū)塊氣藏地層的水平構(gòu)造應(yīng)力系數(shù)，形成最終該地區(qū)的地應(yīng)力解釋模型。

2.1 BP神經(jīng)網(wǎng)絡(luò)算法

BP神經(jīng)網(wǎng)絡(luò)是一種按誤差逆?zhèn)鞑ニ惴ㄓ?xùn)練的多層前饋網(wǎng)絡(luò)。學(xué)習(xí)過程中由信號的正向傳播與誤差的逆向傳播兩個過程組成。正向傳播時,各神經(jīng)元從輸入層開始，接收前一級輸入，經(jīng)激活函數(shù)計算后并輸入到下一級，直至輸出層。反向傳播時,將輸出層輸出值與期望值的誤差逐層返回,并不斷修改權(quán)值矩陣。權(quán)值不斷修改的過程,也就是網(wǎng)絡(luò)學(xué)習(xí)的過程，此過程一直進行到網(wǎng)絡(luò)輸出的誤差逐漸減少到可接受的程度或達到設(shè)定的學(xué)習(xí)次數(shù)為止。三層網(wǎng)絡(luò)結(jié)構(gòu)如圖1所示，其中[X1,X2,…,Xn]為輸入樣本，[Y1,Y2,…,Ym]為輸出樣本。

圖1 三層BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

令隱層權(quán)值矩陣為Wij，輸出層權(quán)值矩陣為Tjk，則：

隱層各單元的輸出按下式計算：

Oj=f(∑WijXi-θj)

(9)

輸出層的輸出結(jié)果為：

Yj=f(∑TjkOj-θk)

(10)

式中，θ表示神經(jīng)元閾值，f表示非線性函數(shù)。隨后按BP神經(jīng)網(wǎng)絡(luò)的梯度下降準(zhǔn)則，定義誤差函數(shù)ek后，進行權(quán)值更新，直到誤差收斂至期望時訓(xùn)練結(jié)束：

(11)

Wjk=Wjk+ηHjek

(12)

式中，i=1,2,…,n；j=1,2,…,l；k=1,2,…,m。

2.2 改進的神經(jīng)網(wǎng)絡(luò)算法

雖然BP神經(jīng)網(wǎng)絡(luò)具有誤差反向傳遞的特征，但由于學(xué)習(xí)速率固定，仍存在以下顯著不足：

(1)若學(xué)習(xí)速率過快，有時會出現(xiàn)“鋸齒形現(xiàn)象”，導(dǎo)致網(wǎng)絡(luò)難以收斂；

(2)若學(xué)習(xí)速率過慢，即使網(wǎng)絡(luò)成功收斂，但首先收斂速度過慢，其次因為誤差平方和函數(shù)易陷入局部極小值，可能得到的只是局部最優(yōu)解[11]。

鑒于常數(shù)型學(xué)習(xí)率的不足，目前為了解決此問題常采用自適應(yīng)學(xué)習(xí)率算法。例如，Adagrad算法中網(wǎng)絡(luò)每代的學(xué)習(xí)率利用初始學(xué)習(xí)率除以歷代梯度平方和的均方根，根據(jù)歷史梯度的變化情況調(diào)整大小達到實現(xiàn)自適應(yīng)效果[12]，即：

(13)

其中，Gt為累積歷史梯度的平方，即：

Gt=Gt-1+gt2

(14)

由于Adagrad算法的全局學(xué)習(xí)速率除以的是累加梯度的平方，到后面累加的比較大時，會導(dǎo)致每次的學(xué)習(xí)率逐漸減小，導(dǎo)致梯度更新緩慢，最終可能會導(dǎo)致訓(xùn)練中后期就停止迭代。為了克服其局限性，該文將Adagrad算法中的學(xué)習(xí)率通過對過去累積的梯度賦予不同的權(quán)重，對學(xué)習(xí)率改進的數(shù)學(xué)表達式如下：

(15)

(16)

該改進的好處在于，在充分考慮了為權(quán)重設(shè)計不同的學(xué)習(xí)率的基礎(chǔ)上，分子上利用冪指數(shù)函數(shù)保證學(xué)習(xí)率穩(wěn)定下降且計算量較小,分母上利用加權(quán)系數(shù)自適應(yīng)學(xué)習(xí)速率對最新的隨機梯度估計賦予了更多的權(quán)重，從而把握梯度變化的方向而加快收斂，避免震蕩，此外還可以避免Adagrad算法關(guān)于考慮歷史梯度的問題。

3 地應(yīng)力預(yù)測模型構(gòu)建

該文地應(yīng)力預(yù)測模型構(gòu)建包括頁巖地應(yīng)力計算模型的構(gòu)建和神經(jīng)網(wǎng)絡(luò)的反分析兩部分。神經(jīng)網(wǎng)絡(luò)的反分析即利用神經(jīng)網(wǎng)絡(luò)對初始地應(yīng)力進行反演得到邊界條件預(yù)測值，即模型中的參數(shù)。選擇目標(biāo)層段內(nèi)若干點的實測應(yīng)力值作為輸入，通過網(wǎng)絡(luò)反分析，得到區(qū)域的模型參數(shù)，然后建立上述地應(yīng)力預(yù)測模型，利用測井曲線獲得的聲波波速計算巖石力學(xué)參數(shù)，最后將巖石力學(xué)參數(shù)，模型參數(shù)代入模型內(nèi)進行計算，即可得到目標(biāo)區(qū)域內(nèi)任意一點的應(yīng)力值。具體步驟如下所示。

3.1 測井資料處理解釋

該文的初始數(shù)據(jù)來源為目標(biāo)工區(qū)裸眼井陣列聲波測井資料。首先對實際的波形數(shù)據(jù)進行了處理：根據(jù)儀器特性，選擇合適的時窗長度，使得在相關(guān)系數(shù)等值線圖上可有效區(qū)分各個模式波對應(yīng)的相關(guān)系數(shù)峰值；根據(jù)各模式波在全波波形上的持續(xù)時間，確定時窗尋峰的開始時間和結(jié)束時間以及時窗移動的步長；通過分析目的層段地質(zhì)、測井等資料，估算各模式波時差的大致范圍，從而確定尋峰的最小時差、最大時差以及時差變化的步長；根據(jù)上述確定的時窗長度、時窗尋峰的開始時間、時窗尋峰的結(jié)束時間等參數(shù)計算相應(yīng)的相關(guān)系數(shù)；由相關(guān)系數(shù)最大值對應(yīng)的時間和時差得到模式波的波至?xí)r間以及模式波的時差，最后得出聲波波速數(shù)據(jù)。

聲波波速確定后，接下來確定剛度矩陣中的C11,C13,C33,C44和C66五個參數(shù)。C33和C44可通過波速資料直接求取，由于缺少斯通利波數(shù)據(jù)，C66通過MANNIE3模型[13]的橫波和縱波各向異性參數(shù)之間的相關(guān)性進行計算，C11,C13用標(biāo)準(zhǔn)的ANNIE模型[14]進行計算，最終利用剛度矩陣元素根據(jù)相應(yīng)巖石力學(xué)參數(shù)經(jīng)驗公式[15-16]計算水平、豎直向的楊氏模量、泊松比等其他巖石力學(xué)參數(shù)。

3.2 確定神經(jīng)網(wǎng)絡(luò)訓(xùn)練樣本

按上述步驟計算后，模型中還缺少最大、最小水平主應(yīng)力方向的構(gòu)造運動系數(shù)，此時需建立實測應(yīng)力值與模型參數(shù)的映射，從而獲得模型參數(shù)，即地應(yīng)力的反演。首先選擇5個地應(yīng)力實測點，根據(jù)提取的聲波資料值計算相應(yīng)的巖石力學(xué)參數(shù)，其次根據(jù)目標(biāo)區(qū)塊的地質(zhì)分析，確定兩個模型中待反演的參數(shù)A,B的上下限，并取5個均布的水平。最后將巖石力學(xué)參數(shù)條件按5×5正交與均勻設(shè)計原理進行組合，通過正演計算計算出相應(yīng)條件下的地應(yīng)力值，從而建立樣本。

3.3 神經(jīng)網(wǎng)絡(luò)特征選取與結(jié)構(gòu)設(shè)計

(17)

式中，h為隱層節(jié)點個數(shù)，m為輸入層節(jié)點個數(shù)，n為輸出層節(jié)點個數(shù)，α為調(diào)節(jié)常數(shù)。網(wǎng)絡(luò)訓(xùn)練采用sigmoid函數(shù)；由于輸出層的神經(jīng)元個數(shù)與實際輸出的個數(shù)相對應(yīng)，取純線性傳遞函數(shù)purelin根據(jù)隨機抽取的方法將樣本數(shù)據(jù)分為訓(xùn)練集和測試集。訓(xùn)練集用于應(yīng)力值與巖石力學(xué)參數(shù)的擬合，測試集用于對網(wǎng)絡(luò)模型計算均方誤差和篩選活動神經(jīng)元。應(yīng)用以上樣本，初始化權(quán)值、閾值矩陣、最大迭代次數(shù)，學(xué)習(xí)率，對網(wǎng)絡(luò)進行訓(xùn)練，采用均方誤差函數(shù)作為回歸損失函數(shù)，即：

(18)

式中，N表示樣本個數(shù)，y'表示神經(jīng)網(wǎng)絡(luò)預(yù)測輸出樣本，y表示真實輸出。當(dāng)誤差函數(shù)下降至期望時迭代結(jié)束。

3.4 神經(jīng)網(wǎng)絡(luò)反演結(jié)果處理

將最優(yōu)化構(gòu)造運動參數(shù)施加在已建立的橫觀各向同性模型建立正分析模型上進行正分析計算，最后與所選測點的實測應(yīng)力值進行對比。

實驗流程如圖2所示。

圖2 實驗流程

4 工程實例應(yīng)用

針對西南油氣田某區(qū)塊內(nèi)的龍馬溪組頁巖進行地應(yīng)力分布預(yù)測。龍馬溪組頁巖為脆性層理頁巖，應(yīng)力敏感性系數(shù)差異較大,各向異性明顯，巖體較為致密，巖石力學(xué)強度較高，頁巖氣含量為5.09～7.48 m3/t，并且隨深度增大而增大。該地區(qū)歷經(jīng)的多期次構(gòu)造活動對區(qū)域頁巖氣的保存產(chǎn)生了重大影響，是頁巖氣的有利保存區(qū)，因此具有一定的研究價值。

4.1 樣本數(shù)據(jù)預(yù)處理

根據(jù)目標(biāo)區(qū)域現(xiàn)場地質(zhì)測試分析，可確定若干個測點的實測應(yīng)力值與巖石力學(xué)參數(shù)，取5個均布的測點，根據(jù)巖石力學(xué)參數(shù)并反算出這些測點的模型中參數(shù)值(見表1)并組合，將巖石力學(xué)參數(shù)、模型中的參數(shù)按照正交設(shè)計表進行設(shè)計，共得到25組樣本(見表2)，并根據(jù)上述公式計算得到結(jié)果(見表3)。

表1 初始巖石力學(xué)參數(shù)

表2 正交設(shè)計計劃表(部分輸出樣本)

表3 部分輸入樣本

4.2 實驗及結(jié)果對比

該文基于Matlab軟件對三層BP神經(jīng)網(wǎng)絡(luò)進行建模。將網(wǎng)絡(luò)的最大迭代次數(shù)設(shè)置為1 000，誤差期望值設(shè)置為10-4。隱藏節(jié)點個數(shù)取輸入單元數(shù)的(2n+1)倍，該文根據(jù)折中的思想，將隱含節(jié)點個數(shù)分別設(shè)置為13、15、17、19、21分別進行仿真實驗，經(jīng)比較后隱層選擇17個節(jié)點進行網(wǎng)絡(luò)訓(xùn)練，網(wǎng)絡(luò)的誤差函數(shù)收斂情況如圖3所示。

圖3 誤差函數(shù)收斂曲線對比

由于200代之后誤差未發(fā)生顯著變化，則圖3只保留200代前的收斂情況。圖3顯示，普通的BP神經(jīng)網(wǎng)絡(luò)的誤差函數(shù)由于學(xué)習(xí)率固定出現(xiàn)了提前收斂至10-2并不再下降的情況，而改進的自適應(yīng)學(xué)習(xí)率算法的網(wǎng)絡(luò)收斂速度較快，均方誤差值穩(wěn)定下降，當(dāng)?shù)^程進行到80代之后，誤差趨于穩(wěn)定并達到期望值，隱含層最優(yōu)中心權(quán)值已獲得。

網(wǎng)絡(luò)訓(xùn)練結(jié)束后，將第21～25組樣本作為測試樣本代入數(shù)值模擬模型進行正演計算，將計算結(jié)果作為模擬實測值，再代入網(wǎng)絡(luò)進行反分析，即可模擬出實測值下的巖石力學(xué)參數(shù)與模型參數(shù)。反演結(jié)果如表4所示。

表4 參數(shù)反演結(jié)果

將上述得到的反演結(jié)果代入到模型中，重新進行正演計算，并進行5個測點的預(yù)測值與現(xiàn)場實測數(shù)據(jù)的結(jié)果對比與誤差計算，結(jié)果如表5所示。

表5 不同的BP神經(jīng)網(wǎng)絡(luò)反演參數(shù)計算結(jié)果對比

表5顯示，兩種訓(xùn)練方法得到的結(jié)果是接近的，利用BP神經(jīng)網(wǎng)絡(luò)反演結(jié)果的最大相對誤差為16.58，最小相對誤差為0.01；利用自適應(yīng)學(xué)習(xí)率網(wǎng)絡(luò)進行反演的最大相對誤差為11.4，最小相對誤差為0.1。根據(jù)國內(nèi)外資料顯示，初始地應(yīng)力值測量的允許誤差達到25%～30%[17]。雖然2種方法的計算精度都能夠較好地滿足工程要求，但改進BP神經(jīng)網(wǎng)絡(luò)的相對誤差較小(除個別測點外)，精度較高，計算收斂代數(shù)少，網(wǎng)絡(luò)性能得到很大程度的提高。所以該改進有一定的實際意義。

5 結(jié)束語

利用BP神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)能力與橫觀各向同性模型相結(jié)合，成功建立了頁巖地應(yīng)力計算模型，并以長寧某工區(qū)的測井?dāng)?shù)據(jù)為基礎(chǔ)，進行了仿真實驗，得到以下結(jié)論：

(1)BP神經(jīng)網(wǎng)絡(luò)可以建立巖石力學(xué)參數(shù)與應(yīng)力值之間的非線性映射關(guān)系,輸出結(jié)果可用來計算初始地應(yīng)力。

(2)自適應(yīng)學(xué)習(xí)率的BP神經(jīng)網(wǎng)絡(luò)具有更好的學(xué)習(xí)能力，收斂速度較快、結(jié)果誤差較小。

(3)所建立的網(wǎng)絡(luò)模型對地應(yīng)力邊界條件的預(yù)測誤差較低，將該模型用于初始地應(yīng)力場預(yù)測，在誤差允許的范圍內(nèi)可以代替地應(yīng)力實測。因此，該模型在初始地應(yīng)力場預(yù)測方面具有良好的適用性與可行性，具有一定的使用價值。