基于改進BP神經網絡的頁巖地應力預測模型

尚福華，王瑋卿，曹茂俊

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

0 引 言

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

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

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

1 頁巖地應力物理模型

1.1 橫觀各向同性物理模型

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

(1)

(2)

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

(3)

(4)

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

(5)

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

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

(6)

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

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

1.2 改進的頁巖地應力計算模型

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

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

(7)

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

(8)

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

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

2 BP神經網絡模型

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

2.1 BP神經網絡算法

BP神經網絡是一種按誤差逆傳播算法訓練的多層前饋網絡。學習過程中由信號的正向傳播與誤差的逆向傳播兩個過程組成。正向傳播時,各神經元從輸入層開始，接收前一級輸入，經激活函數計算后并輸入到下一級，直至輸出層。反向傳播時,將輸出層輸出值與期望值的誤差逐層返回,并不斷修改權值矩陣。權值不斷修改的過程,也就是網絡學習的過程，此過程一直進行到網絡輸出的誤差逐漸減少到可接受的程度或達到設定的學習次數為止。三層網絡結構如圖1所示，其中[X1,X2,…,Xn]為輸入樣本，[Y1,Y2,…,Ym]為輸出樣本。

圖1 三層BP神經網絡結構

令隱層權值矩陣為Wij，輸出層權值矩陣為Tjk，則：

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

Oj=f(∑WijXi-θj)

(9)

輸出層的輸出結果為：

Yj=f(∑TjkOj-θk)

(10)

式中，θ表示神經元閾值，f表示非線性函數。隨后按BP神經網絡的梯度下降準則，定義誤差函數ek后，進行權值更新，直到誤差收斂至期望時訓練結束：

(11)

Wjk=Wjk+ηHjek

(12)

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

2.2 改進的神經網絡算法

雖然BP神經網絡具有誤差反向傳遞的特征，但由于學習速率固定，仍存在以下顯著不足：

(1)若學習速率過快，有時會出現“鋸齒形現象”，導致網絡難以收斂；

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

鑒于常數型學習率的不足，目前為了解決此問題常采用自適應學習率算法。例如，Adagrad算法中網絡每代的學習率利用初始學習率除以歷代梯度平方和的均方根，根據歷史梯度的變化情況調整大小達到實現自適應效果[12]，即：

(13)

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

Gt=Gt-1+gt2

(14)

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

(15)

(16)

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

3 地應力預測模型構建

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

3.1 測井資料處理解釋

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

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

3.2 確定神經網絡訓練樣本

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

3.3 神經網絡特征選取與結構設計

(17)

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

(18)

式中，N表示樣本個數，y'表示神經網絡預測輸出樣本，y表示真實輸出。當誤差函數下降至期望時迭代結束。

3.4 神經網絡反演結果處理

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

實驗流程如圖2所示。

圖2 實驗流程

4 工程實例應用

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

4.1 樣本數據預處理

根據目標區域現場地質測試分析，可確定若干個測點的實測應力值與巖石力學參數，取5個均布的測點，根據巖石力學參數并反算出這些測點的模型中參數值(見表1)并組合，將巖石力學參數、模型中的參數按照正交設計表進行設計，共得到25組樣本(見表2)，并根據上述公式計算得到結果(見表3)。

表1 初始巖石力學參數

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

表3 部分輸入樣本

4.2 實驗及結果對比

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

圖3 誤差函數收斂曲線對比

由于200代之后誤差未發生顯著變化，則圖3只保留200代前的收斂情況。圖3顯示，普通的BP神經網絡的誤差函數由于學習率固定出現了提前收斂至10-2并不再下降的情況，而改進的自適應學習率算法的網絡收斂速度較快，均方誤差值穩定下降，當迭代過程進行到80代之后，誤差趨于穩定并達到期望值，隱含層最優中心權值已獲得。

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

表4 參數反演結果

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

表5 不同的BP神經網絡反演參數計算結果對比

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

5 結束語

利用BP神經網絡的學習能力與橫觀各向同性模型相結合，成功建立了頁巖地應力計算模型，并以長寧某工區的測井數據為基礎，進行了仿真實驗，得到以下結論：

(1)BP神經網絡可以建立巖石力學參數與應力值之間的非線性映射關系,輸出結果可用來計算初始地應力。

(2)自適應學習率的BP神經網絡具有更好的學習能力，收斂速度較快、結果誤差較小。

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