基于智能算法的油氣田地應力三維預測

袁 多，吳 超，盧運虎，張東清

1.頁巖油氣富集機理與有效開發(fā)國家重點實驗室，北京 昌平 102206

2.中國石化石油工程技術研究院，北京 昌平 102206

3.中國石油大學（北京）石油工程學院，北京 昌平 102249

引言

油氣田地應力三維預測是將已鉆井的地應力信息延拓外推至區(qū)域地層空間的運算過程，亦稱地應力空間建模，它有助于準確認識地應力在三維區(qū)域內的縱、橫向分布規(guī)律，從而在井壁失穩(wěn)控制、套管損傷和出砂預防、壓裂工藝優(yōu)化等方面發(fā)揮指導作用[1-2]。由于油氣儲層埋藏相對較深，鉆井與開發(fā)工程所涉及的地層空間范圍較廣，其多變的沉積與構造環(huán)境導致了地層具有復雜的巖石物理和力學特征，因此，針對油氣田地層空間進行地應力預測具有更大的難度。

在油氣田多種數(shù)據(jù)資料中，可用于較準確求取地應力的信息主要為實驗測試數(shù)據(jù)以及與應力存在直接關聯(lián)的聲波測井、密度測井等數(shù)據(jù)，這些信息都分布于工區(qū)內的完鉆井處。因此，油氣田巖層空間的地應力建模以完鉆井的測試和測井資料為基礎，在三維區(qū)域內的地質力學和地球物理信息約束下，運用數(shù)學手段將單井地應力數(shù)據(jù)推廣為三維空間的地應力數(shù)據(jù)體[3]。事實上，這種面向未鉆地層的地應力預測對于后續(xù)鉆井的指導意義更為重要。

通過巖芯和壓裂實驗可精準確定地應力大小[4-5]，但油氣田深層取芯和測試費用較高，且只能獲得井點個別深度的應力信息。因而利用測井資料進行分層地應力計算應用更為廣泛，其中基于地層各向同性假設的求取方法包括早期單軸應變法、黃氏方法、組合彈簧法、Holbrook 法和葛氏方法等[6]，近年各向異性模型也得到了完善[7-8]。由于資料來源制約，測井分析方法通常在鉆后進行，而且它只能反映井壁圍巖應力狀況，對于油氣工程尤其是鉆井的指導缺乏時效性和針對性。趙海峰等提出了利用實鉆資料反演地應力的方法，但這類井周應力實時監(jiān)測技術無法用于區(qū)域地應力預測[9]。對于油氣田地應力的三維預測問題，當前常規(guī)處理手段是借鑒礦業(yè)、水利等學科的地應力反演理論[10]，將測試、測井數(shù)據(jù)和有限元模擬等數(shù)值算法結合起來進行地應力預測[3，11-13]，但在建模計算中沒有使用三維地震信息進行約束控制，地應力預測精度和分辨率難以滿足精細勘探開發(fā)的需求。劉建偉等基于疊前彈性反演理論進行了頁巖油藏地應力的三維地震預測[14]，但資料處理和計算量的限制導致該技術在地層比例超過90%的非目的層難以推廣應用。近期有學者嘗試使用疊后地震數(shù)據(jù)預測三維地應力以提高運算時效[15-17]，但其操作步驟繁多，且用于應力解空間搜索的全隨機算法易陷入局部最優(yōu)，求解效率有待提高。斯倫貝謝公司的地震指導鉆井技術（SGD）中通過地震層析反演預測強度、應力、壓力等地層力學參數(shù)[18-20]，同樣由于運算量的制約，該方法僅針對鉆頭前方的小范圍空間。

本文通過考察巖石物理力學參數(shù)與地震信息之間的關系，針對不同的工況條件，運用BP 神經(jīng)網(wǎng)絡、模擬退火等智能算法建立了一套較全面的油氣田地應力場三維預測技術體系。該技術充分利用地震信息覆蓋面廣、信息量大的優(yōu)勢，以提高地應力預測精度和分辨率，同時通過優(yōu)化算法和處理流程來提高操作性能。

1 間接預測法及其應用實例

根據(jù)地震勘探理論，疊后地震記錄數(shù)據(jù)可以表示為各套地層界面處的反射系數(shù)和地震子波的褶積，即

地震波反射系數(shù)通過反射界面兩側介質的波速和密度來計算

一般說來，地層縱波速度和密度之間存在良好的經(jīng)驗關系

由式（1）～式（3），可以將疊后地震記錄和層速度之間的關系表示為如下函數(shù)關系

眾所周知，巖石的彈性模量和泊松比數(shù)值由地層密度與聲波速度決定

大量巖石物理實驗結果表明，巖石縱波和橫波傳播速度之間存在一定的經(jīng)驗關系

目前國內外已提出多種油氣田地應力分層計算模式，各類模式中普遍涉及的影響因素為彈性模量、泊松比、孔隙壓力及構造應力系數(shù)，以得到廣泛使用的黃氏修正模式為例

利用密度測井數(shù)據(jù)積分可得到適用于目標工區(qū)的上覆壓力模式，構造應力系數(shù)通常由實驗室或現(xiàn)場測試結果擬合得到，孔隙壓力與聲速關系非常密切，常用的Eaton 計算模式

根據(jù)式（5）～式（8），層速度與地應力之間的函數(shù)關系可表示為

通過式（9），可以發(fā)現(xiàn)用于地應力建模的最關鍵、最基礎的數(shù)據(jù)是地層的縱波傳播速度。式（4）已表明層速度和疊后地震記錄之間存在定量關系。

因此，可以運用數(shù)學方法通過對地震數(shù)據(jù)反推求得速度信息，這在地球物理學中稱為地震速度反演。傳統(tǒng)的反演過程是基于Levenberg-Marguardt理論求取層速度的最小二乘解[21]，而式（4）在解析構建偏導數(shù)矩陣時巨大的工作量大大降低了運算時效，且初始解選取的嚴格性也難以保證穩(wěn)定的迭代過程。

考慮到BP 神經(jīng)網(wǎng)絡在模式識別、數(shù)值逼近、分類等智能學習方面具有很高的穩(wěn)定性與可靠性，因此筆者使用這種多層前饋網(wǎng)絡算法來識別層速度和地震記錄之間的函數(shù)映射關系[22]。

首先，收集工區(qū)完鉆井的聲波測井及井旁地震記錄，建立時深關系以將地震、測井數(shù)據(jù)統(tǒng)一標度，并從地震記錄中提取瞬時參數(shù)、頻譜參數(shù)、功率譜參數(shù)、自相關參數(shù)、自回歸參數(shù)、分形參數(shù)等特征參數(shù)等。

分層段運用聚類算法將地震特征參數(shù)和聲波速度建成學習樣本對，并將其分別導入BP 神經(jīng)網(wǎng)絡的輸入和輸出端。

基于激勵函數(shù)，利用式（10）依次計算網(wǎng)絡各節(jié)點（神經(jīng)元）的輸出

設置目標函數(shù)使網(wǎng)絡輸出聲波速度和學習樣本對中的實測速度值之間達最小誤差，選取目標函數(shù)E0，運用梯度下降法對網(wǎng)絡連接權值與閾值逐次迭代修改

重復以上運算直到某次計算的E0達到指定精度，可結束網(wǎng)絡學習過程。

基于上述分層學習成果，將工區(qū)內未鉆空間位置處的地震特征參數(shù)按照對應的地層輸入BP 神經(jīng)網(wǎng)絡，輸出得到相應位置的層速度預測值，即可建立層速度的三維空間數(shù)據(jù)體。

在預測三維層速度的基礎上根據(jù)工區(qū)實際地質情況選取力學模型進一步預測區(qū)域地應力[6-8]。上述思路是在層速度地震反演的基礎上進一步預測地應力，因此稱之為間接預測法。由于該方法所需的學習樣本和經(jīng)驗參數(shù)較多，其應用條件是工區(qū)具有豐富的實測資料用于參數(shù)統(tǒng)計回歸及BP 神經(jīng)網(wǎng)絡學習。為盡量避免BP 神經(jīng)網(wǎng)絡陷入局部極值，借鑒了Matlab 工具函數(shù)trainbfg 采用的基于擬牛頓法的訓練算法，實際收斂速度和精度均很理想。

渤海灣某油田斷塊發(fā)育，地應力狀態(tài)復雜，井壁失穩(wěn)、出砂、套損等復雜故障頻繁發(fā)生。

由于該油田LZ 工區(qū)完鉆井實測信息豐富，所以在該工區(qū)應用了間接預測法進行地應力三維預測。通過實測數(shù)據(jù)的統(tǒng)計分析，確定了模型參數(shù)，其中工區(qū)分層構造應力系數(shù)見表1。

表1 LZ 工區(qū)分層構造應力系數(shù)Tab.1 The tectonic stress coefficients of different strata in LZ Area

收集到LZ 工區(qū)的三維疊后地震數(shù)據(jù)體，圖1顯示的剖面是其中的line1518，基于BP 神經(jīng)網(wǎng)絡算法預測三維層速度，圖2 顯示的速度剖面是其中的line1518，進一步預測三維地應力（圖3）。

圖1 LZ 工區(qū)的疊后地震數(shù)據(jù)體Fig.1 Post stack seismic data volume in LZ Area

圖2 基于BP 神經(jīng)網(wǎng)絡算法預測LZ 工區(qū)縱波速度數(shù)據(jù)體Fig.2 Predicted interval velocity data volume based on BP neural network algorithm in LZ Area

圖3 預測得到的LZ 工區(qū)最大水平地應力數(shù)據(jù)體Fig.3 Predicted maximum horizontal stress data volume in LZ Area

圖4 和圖5 為從地應力數(shù)據(jù)體中切片得到的分層段橫向展布等值線圖。

由等值線圖可發(fā)現(xiàn)LZ 工區(qū)地應力呈從北向南逐漸增大的區(qū)域特征，另外可通過觀測地應力等值線密集程度來大體確定工區(qū)斷層的走向[11]，可判斷L22、L16、L36、L32、L3-3 等井附近發(fā)育5 組斷塊（以上推斷與地質認識基本吻合）。以上分析成果可用于后續(xù)鉆井的井壁失穩(wěn)、出砂、套損預防及開發(fā)方案優(yōu)化設計。

2 直接預測法及其應用實例

在工區(qū)完鉆井數(shù)量與實測數(shù)據(jù)較少的情況下，運用間接預測法可能存在一些問題：（1）井點處可信的測試信息缺乏，式（7）中的構造應力系數(shù)難以合理確定。（2）用于建立學習樣本的測井資料較少，導致BP 神經(jīng)網(wǎng)絡模式識別及外推能力減弱。針對以上問題，可以考慮變換思路，由式（9）反推，建立以地應力及構造應力系數(shù)為自變量的函數(shù)關系式

將式（12）和式（4）結合，可以導出

式中：n每個地震道上地震記錄序號；

m每個地震道上待求解的地應力對應的反射界面序號；

Sn每個地震道上的各個地震記錄。

式（13）是一個地震道上的以地應力及構造應力系數(shù)為未知量的非線性方程組，可根據(jù)其通過地震數(shù)據(jù)直接反演求取三維區(qū)域內各地震道上的水平地應力解向量，也可同時解得構造應力系數(shù)。

理論上可以使用非線性優(yōu)化理論求解式（13），通過尋求合成與實際地震記錄的最優(yōu)匹配來直接求解地應力，避免經(jīng)驗系數(shù)不準確引起的誤差，因此稱之為直接預測法。

但在井點實測信息缺乏的情況下，若按照常規(guī)的非線性最小二乘法求解式（13），將會因為方程組規(guī)模大、非線性程度高、測井插值所得的初始解質量較差而導致求解過程穩(wěn)定性極差、運算時間過長，所以筆者考慮使用智能優(yōu)化理論中的模擬退火算法對其進行求解。

在完鉆井資料缺乏的工況下，利用該算法求解式（13），具有兩個突出優(yōu)點：（1）對初始解要求不嚴格，概率重要性采樣算法可以保證穩(wěn)定的尋優(yōu)運算過程。（2）在地震道上進行全局尋優(yōu)，減少得到局部最優(yōu)解的可能性。模擬退火屬于智能算法，它模仿固體的退火過程，在緩慢降溫時每個溫度下內部粒子達平衡態(tài)，最終降至低溫時內能最小[23]。根據(jù)式（13）的求解要求，設定目標函數(shù)。

同時設置模擬退火過程的初始溫度T0、最低溫度Tf，溫度降低策略采用常規(guī)的Tj+1=λTj（0 ＜λ ＜1，j為降溫次數(shù)）。具體運算過程如下：

計算并分析工區(qū)完鉆井的水平地應力數(shù)據(jù)，統(tǒng)計其所服從分布類型（一般為正態(tài)分布），確定各參數(shù)候選解搜索鄰域空間并建立各自的分布密度函數(shù)，以隨機方式產(chǎn)生某地震道上的初始解σ0，并在此后迭代過程中對解空間的隨機采樣，逐次獲得候選解。

在初始溫度下T0，計算初始解目標函數(shù)Φ0，然后在鄰域內隨機產(chǎn)生一個候選解σ1，計算其目標函數(shù)Φ1。

如果Φ1≤Φ0，則接受σ1為當前最優(yōu)解；否則需計算（K--Boltzmann 常數(shù)），并產(chǎn)生一個隨機數(shù)R∈[0，1]，并運用Metropolis 準則：若P＞R，則也接受σ1為當前最優(yōu)解，否則不接受σ1。

在初始溫度下反復產(chǎn)生候選解，并按上述方式不斷更新當前最優(yōu)解。

當連續(xù)若干步目標函數(shù)值無明顯變化或候選解搜索次數(shù)達到設定馬爾可夫鏈長度L，該溫度下的迭代運算結束。

按照降溫策略逐次降溫，每個溫度Tj下均執(zhí)行上一步的操作（需注意計算P時需更新Tj）。

直至溫度降低至最低溫度Tf或連續(xù)兩次搜尋到的最優(yōu)解之差小于指定容差ε，全部迭代尋優(yōu)運算結束，得到本地震道上的水平地應力最終優(yōu)化解，進入下一地震道重復以上求解過程，直至解得三維巖層空間上的地應力值。

必須說明的是，式（13）的求解涉及到地層密度、上覆壓力、孔隙壓力計算模式中的經(jīng)驗系數(shù)，求取這些系數(shù)相對容易，且通用性好，因此不再作為未知量參與求解。

但由于在直接預測法中沒有事先預測層速度，即無法通過層速度進一步預測孔隙壓力，而孔隙壓力又是式（13）求解的必需參數(shù)，所以可參照上述地應力直接求解的思路，利用地震數(shù)據(jù)先期求取孔隙壓力。

與式（13）的推導過程類似，直接通過地震信息求取孔隙壓力

式（14）的求解過程和式（13）基本一致，這里不再贅述。在求得三維孔隙壓力數(shù)據(jù)后，將其作為求解方程組（13）的已知信息。

式（13）和式（14）的求解需要分地層進行，對于同套地層，構造應力系數(shù)作為固定值參與求解。

油田的NP-2 工區(qū)屬新探區(qū)，巖芯與現(xiàn)場測試及測井數(shù)據(jù)較少，前期使用有限元建模得到的地應力平均計算誤差約22%，因此，在該工區(qū)引入了地震信息進行區(qū)域地應力計算分析。

圖6 為NP-2 工區(qū)的三維疊后地震數(shù)據(jù)體（以其中201 地震測線為例對其二維剖面進行顯示）。

圖6 NP-2 工區(qū)的疊后地震數(shù)據(jù)體Fig.6 Post stack seismic data volume in NP-2 Area

在運算設計時，充分考慮了快速模擬退火策略[24-25]，以期提高其現(xiàn)場應用時的運算速度。

經(jīng)過迭代運算后通過地震信息直接得到了最大和最小水平地應力三維數(shù)據(jù)體（圖7、圖8）。

圖7 NP-2 工區(qū)預測最大水平地應力數(shù)據(jù)體剖面Fig.7 Predicted maximum horizontal stress data volume in NP-2 Are a

圖8 預測得到的NP-2 工區(qū)最小水平地應力數(shù)據(jù)體Fig.8 Predicted minimum horizontal stress data volume in NP-2 Area

圖9 和圖10 分別為從三維地應力數(shù)據(jù)體中切片得到的分層段橫向展布等值線圖。

圖9 NP-2 區(qū)塊Nd1 段最大水平地應力切片等值線圖Fig.9 Slice contour map of maximum horizontal stress of Nd1 Section in NP2 Area

圖10 NP-2 區(qū)塊Ng2 段最小水平地應力切片等值線圖Fig.10 Slice contour map of minimum horizontal stress of Ng2 Section in NP-2 Area

從圖上可發(fā)現(xiàn)，NP-2 工區(qū)地應力具有北部大、南部小的整體特征，同時從等值線疏密程度判斷出NP2-13、NP2-2、NP2-10 井附近發(fā)育3 組次級斷層，且大體呈NE45°走向，實際地質勘探成果驗證了以上推斷。上述間接和直接預測方法在該油田5 個工區(qū)根據(jù)實際工況進行了選擇應用，在LZ、MYL、GSP 等老工區(qū)應用了間接預測法；而NP-1、NP-2 是新探構造，適用直接預測法。為了全面驗證預測結果，從各工區(qū)的實際測試結果中取出部分樣本，它們不參與經(jīng)驗參數(shù)統(tǒng)計擬合，僅作為驗證數(shù)據(jù)，將其與相同空間位置處的預測值進行對比。從表2 可以看到，地應力整體預測精度基本可滿足工程需要，其中直接預測法的計算精度總體上低于間接預測法，因直接預測所使用的約束信息相對較少，但其平均預測精度為90.73%，仍明顯高于單純力學建模方法（78.00%），同時剖面和切片圖也可達到較高分辨率。

表2 多個試驗工區(qū)預測與實測地應力數(shù)據(jù)對比Tab.2 Comparison between predicted and measured in-situ stress data of several test areas

另外，通過與前期試用過的疊前地震預測算法對比，本技術合理利用了疊后地震信息，在現(xiàn)場應用中的資料處理和計算量大大減少，優(yōu)化求解運算過程穩(wěn)定，操作流程更為簡便，不僅可用于鉆前地質力學特征綜合分析，還可用于隨鉆井壁穩(wěn)定預測。以上情況表明，該預測技術在油氣田現(xiàn)場應用具備較高的適用性和可靠性。

3 結論

（1）地層聲波傳播速度是進行油氣田地應力定量分析的基礎數(shù)據(jù)，而疊后地震記錄與層速度之間存在著非線性函數(shù)關系，因此，可運用智能算法解析這種關系，從而利用區(qū)域地震資料預測三維地應力。

（2）對于完鉆井數(shù)量較多、實測信息較豐富的工區(qū)，其區(qū)域地應力計算適用間接預測法，該方法運用BP 神經(jīng)網(wǎng)絡算法識別完鉆處的聲波測井和地震信息之間的關系，并根據(jù)其外推預測三維層速度；同時利用測試數(shù)據(jù)統(tǒng)計相關巖石物理和力學經(jīng)驗參數(shù)，基于以上工作進行地應力三維預測。

（3）在工區(qū)測試和測井資料較缺乏的情況下，應使用直接預測法進行地應力三維建模，該方法以地應力和地震記錄之間的映射關系為基礎，運用模擬退火算法直接搜尋合成與實際地震記錄達最優(yōu)匹配下的地應力解向量，尋優(yōu)運算過程穩(wěn)定，易于得到全局最優(yōu)解。對多數(shù)探區(qū)來說，本方法更具意義。