基于雙向長短期記憶神經網絡的水平地應力預測方法

馬天壽，向國富，石榆帆，桂俊川，張東洋

1 西南石油大學“油氣藏地質及開發工程”國家重點實驗室, 成都 610500

2 西南石油大學工程學院, 南充 637001

3 中國石油西南油氣田公司頁巖氣研究院, 成都 610051

0 引言

準確預測地應力對于油氣勘探開發具有重要作用和意義[1]。地應力是井壁穩定分析的邊界條件，其大小決定了維護井壁穩定所需的當量泥漿密度[2-4]，同時地應力也是水力壓裂設計的關鍵參數，其大小直接約束了水力裂縫的起裂和擴展[5-6]。地應力由垂向地應力、最大和最小水平地應力構成，其中垂向地應力可直接由密度測井資料計算得到，而最大和最小水平地應力由于地質環境和技術水平的限制難以實現高效預測。目前，確定最大和最小水平地應力的方法包括直接測量法、測井解釋法、數值模擬法及地震預測法等[7]。其中，直接測量法和測井解釋法精度較高且技術相對成熟因此被廣泛采用。直接測量法包括聲發射測試、波速各向異性測試、水力壓裂測試、診斷裂縫注入測試等[8-10]，直接測量法精度最好，但是測試成本高、耗時長，而且只能獲得取芯深度點的地應力，所測得的數據十分有限；測井解釋法可以利用聲波測井、成像測井及鉆井誘導裂縫等資料確定地應力[11-12]，相較于直接測量法，測井解釋法精度偏低，但是具有縱向分辨率高、相對連續的特點，可以解釋得到相對連續的地應力剖面，因而被現場廣泛應用。然而，目前的測井解釋方法依賴于陣列聲波測井資料和實測地應力數據，嚴重限制了地應力的測井解釋。此外，由于地下情況復雜和非均質性較強，測井數據之間表現出極強的非線性關系[13]，使得傳統的解析模型和神經網絡模型難以表達測井數據與地應力之間的空間相關性。因此，亟需尋求一種簡單、高效的替代模型用于解決地應力預測問題。

近年來，隨著機器學習在科學和工程領域的廣泛應用，許多學者提出利用機器學習來解決地球科學問題[14-15]。張東曉等[15]提出了一種串級長短期記憶神經網絡，用于補全和生成人工測井曲線；Singh[16]提出了一種模糊邏輯模型進行地層識別；Mehrgini等[17]利用Elman神經網絡，通過訓練自然伽馬、電阻率、中子孔隙度、巖石體積密度及含水飽和度等參數來預測橫波速度；Korjani等[18]基于鄰井地質資料，構建了一種深度學習神經網絡預測巖石物理特征的油藏建模方法，用于構建缺少測井和巖心資料的虛擬測井數據；韓玉嬌等[19]采用AdaBoost機器學習算法，開展了大牛地氣田儲層流體智能識別研究。在地應力預測方面：馮鵬等[20]篩選了井徑、補償中子、自然伽馬、密度、深淺側向電阻率等測井參數作為訓練因子，利用支持向量回歸方法建立了煤層最小水平地應力預測模型；尚福華等[21]基于改進的BP神經網絡模型對巖石力學參數進行反演預測，進而進行地應力預測；HAN和YIN[22]利用人工神經網絡(ANN)，將井眼形變數據作為特征來預測地應力；IBRAHIM等[23]基于測井參數利用隨機森林、功能網絡和模糊推理系統等機器學習技術來預測地應力；GOWIDA等[24]利用ANN，通過訓練鉆壓、扭矩及機械鉆速等實鉆數據來預測地應力；TUNG[25]利用ANN，通過訓練垂深、孔隙壓力、垂向地應力等數據來預測最小水平地應力；PU等[26]通過決策樹，綜合考慮巖石彈性模量、泊松比及其他參數的影響來預測地應力；MA等[27]建立了一種混合神經網絡模型，用于單井地應力的預測。

國內外學者已經利用機器學方法對巖性識別、測井曲線重構、地應力預測等開展了大量研究，但地應力的機器學習方法大都采用傳統全連接神經網絡(FCNN)，這類方法只能構造點對點的映射，無法有效捕捉測井資料深度上的變異信息，而且，這類方法通常假設數據是低維和線性的，因而難以有效準確地預測地應力。此外，這類方法的泛化能力較差，很難用于新井地應力的預測。為此，采用一種先進的循環神經網絡，即雙向長短記憶神經網絡(BiLSTM)，利用常規測井參數來預測地應力。首先，選取工區兩口直井作為研究對象，分別對兩口井的測井資料進行地應力解釋，形成訓練和測試數據集；然后，研究數據預處理與樣本構造算法，分析測井參數與地應力相關性，建立基于BiLSTM的水平地應力預測模型，并闡述模型的核心算法和流程；最后，結合測井參數相關性與地質學含義，分析不同測井參數組合模式下水平地應力的預測效果。本文研究對于準確預測地應力、推廣機器學習在石油工程的應用具有重要的作用和意義。

1 地應力測井解釋

本文數據來自于四川盆地CL氣田的2口直井，工區位于川中隆起區的川西南低陡褶帶，該工區整體表現為由北西向南、東方向傾斜的寬緩單斜構造，且地層平緩，傾角小，斷裂不發育。為了獲得2口井的地應力測井解釋數據集，需要根據工區地層特征建立地質力學模型。

1.1 地質力學模型

地質力學模型是根據室內巖石力學實驗、礦場實測資料和測井解釋建立地層巖石力學參數、孔隙壓力、地應力一維剖面的常用方法[27-28]，其建模流程如圖1所示。首先，確定巖石力學參數，包括泊松比、彈性模量和Biot系數；其次，通過密度測井資料計算垂向地應力；然后，利用聲波測井資料確定孔隙壓力；進一步，根據室內測試或小型地破試驗結果，反演構造應力系數；最后，利用巖石力學參數、孔隙壓力、垂向地應力解釋結果和構造應力系數，確定最大和最小水平地應力。

圖1 地質力學模型建模流程圖[27]Fig. 1 Flowchart of geo-mechanical modeling[27]

大部分油氣測井作業中，由于成本限制，一般只進行縱波時差測井，缺少橫波時差測井資料[29]。因此，需要采用橫波重構來估算橫波時差，即利用工區已測井的縱橫波時差關系，建立橫波時差預測的經驗公式。對于目標工區，縱橫波時差之間的經驗關系為[27]：

式中：Δtp為縱波時差，μs/ft；Δts為橫波時差，μs/ft。

進一步，根據彈性介質中縱橫波傳播理論，建立動態彈性參數與縱波、橫波時差的關系[29]：

式中：Ed為動態彈性模量，GPa；μd為動態泊松比，無因次；ρb為巖石密度，g/cm3。

在實際工程應用中，需要采用靜態彈性模量和泊松比。因此，需要對彈性模量和泊松比進行動靜態轉換。對于目標工區，彈性模量和泊松比的動靜態轉換關系為[27]：

式中：Es為靜態彈性模量，GPa；μs為靜態泊松比，無因次。

而地層巖石Biot系數可表示為[29]：

式中：α為Biot系數，無因次；ρma為巖石骨架密度，g/cm3；Δtmap和Δtmas分別為巖石骨架的縱波時差和橫波時差，μs/ft。

垂向地應力也稱上覆巖層壓力，可由密度測井曲線直接積分求得：

式中：σv為垂向地應力，MPa；ρ為缺少密度測井資料井段的上覆巖層平均密度，g/cm3；TVD0為測井曲線起始點垂深，m；TVD為給定計算測點所對應的垂深，m。

通過對比孔隙壓力預測方法及預測效果，優選Eaton法計算目標工區地層孔隙壓力。Eaton法以泥巖正常壓實理論為基礎，因此，先要建立正常壓實趨勢線，然后才能計算地層孔隙壓力。對于目標工區，其正常壓實趨勢線可表示為[27]：

而Eaton法計算地層孔隙壓力的公式可表示為[30]：

式中：pp為地層孔隙壓力，MPa；pw為地層流體靜液柱壓力，MPa；Δtnp為正常趨勢線上對應的縱波時差，μs/ft；c為Eaton常數，無因次，目標工區Eaton常數取值為0.914。

最后，計算最大和最小水平地應力，水平地應力的計算模式主要分為兩種[31-32]：一種是單軸應變模式，認為水平方向無限大，地層在沉積過程中只發生垂向變形，水平方向的應力是由上覆巖層重量引起的，主要包括Terzaghi模型、Matthews-Kelly模型、Anderson模型和Newberry模型，這類模式沒有考慮構造應力，不符合實際地質條件；另一類是分層計算模式，主要包括組合彈簧模型、斯倫貝謝模型和黃氏模型，該模式考慮了構造應力對水平地應力的影響，因此被廣泛應用于油田地應力預測。由于目標工區整體呈低緩構造特征，因此，采用黃氏模型計算水平地應力，黃氏模型可表示為[32]：

式中：σH為最大水平地應力，MPa；σh為最小水平地應力，MPa；ω1和ω2分別為最大和最小水平地應力方向的構造應力系數，無因次。

在式(8)所示的黃氏模型中，構造應力系數的確定非常關鍵。對確定的工區，其構造應力系數可視為常數。利用室內地應力差應變測試結果，并結合孔隙壓力(pp)、靜態泊松比(μs)、Biot系數(α)測井解釋結果，可以利用式(8)反演出工區的構造應力系數，目標工區構造應力系數取值為ω1=0.458、ω2=0.364，進而再利用式(8)計算連續的水平地應力測井解釋剖面。

1.2 測井解釋結果

通過測井解釋得到了工區2口直井(X-1井和X-2井)的地應力剖面，包括垂向地應力、最大水平地應力和最小水平地應力，結果如圖2所示。由圖可知，每口井有10條測井曲線，包括自然伽馬(GR)、縱波時差(DTC)、井徑(CAL)、密度(DEN)、補償中子(CNL)、橫波時差(DTS)、測深(MD)、方位角(AZIM)、井斜角(DEVI)和垂深(TVD)。其中，X-1井測深為2809～3664.3 m，測段總長為855.5 m，；X-2井測深為3218.1～3693.0 m，測段總長為474.9 m。不難看出：該工區垂向應力＞最大水平地應力＞最小水平地應力，屬于潛在正斷層應力狀態，而且，地應力隨井深的增加呈增大趨勢。此外，為了驗證黃氏模型的準確性，將X-1井地應力解釋結果與3663.2 m、3657.3 m和3669.1 m處巖心差應變實驗結果進行對比發現，X-1井垂向地應力誤差為0.39%，最大水平地應力誤差為0.18%～0.64%，最小水平地應力誤差為0.29%。因此，測井解釋得到的地應力與實測地應力吻合較好，可以將測井解釋得到的地應力作為真實地應力。

圖2 地應力測井解釋成果圖Fig. 2 Logging interpretation results of in-situ stress

2 數據預處理與樣本構造

數據預處理和樣本構造是機器學習建模的重要環節，前者決定了機器學習模型預測和泛化能力的好壞，而后者確定了數據輸入模型的形式。本文分別將X-1井和X-2井作為訓練井和測試井，其中，訓練井用于開發機器學習模型，測試井用于驗證模型的預測和泛化能力。

2.1 數據處理

測井參數是對井下巖石物性特征的綜合響應，不同測井參數能夠反映地層巖性、泥質含量、含水飽和度、滲透率等特征[33]。而地應力取決于地質構造、沉積環境、地層巖性等因素，因此，測井參數與地應力之間存在著必然聯系，這也是利用測井數據預測地應力的理論依據。同時，測井參數與地應力之間關系復雜，不是所有測井參數都與地應力之間存在關聯關系。因此，需要對測井參數進行合理篩選，以消除冗余參數對地應力預測的影響。經典機器學習方法往往采用相關性分析來挖掘數據之間的關聯關系。

考慮到測井參數之間一般呈非正態分布，因此，采用Spearman秩相關系數來表征不同測井參數之間的相關性。Spearman秩相關系數是一種非參數的相關性度量準則，通過對變量X和Y求秩后再進行相關性分析，能夠滿足非正態分布數據的相關性計算，Spearman秩相關系數計算公式為[28]：

式中：RS為Spearman秩相關系數；n為樣本的數量；di為變量X和Y的秩次之差。

同時，根據式(9)計算得到的相關性系數，可以按照以下取值范圍判斷相關性強弱[34]：0.8≤|RS|≤1.0表 示 極 強 相 關；0.6≤|RS|≤0.79表 示 強 相 關；0.4≤|RS|≤0.59表示中等相關；0.2≤|RS|≤0.39表示弱相關；0.0≤|RS|≤0.19表示不相關。

圖3所示為X-1井測井參數相關系分析結果，不難看出：各測井參數與水平地應力之間具有極其相似的相關性系數，但不同測井參數與水平地應力之間的相關性存在較大差異：水平地應力與垂深、密度及自然伽馬等測井參數之間表現為正相關，與縱波時差、井徑、補償中子及橫波時差等測井參數之間表現為負相關；最大和最小水平地應力與垂深之間的相關性系數最大且均為0.98，其次是井徑、密度、補償中子和自然伽馬，分別為-0.38、0.37、-0.36和0.34，縱波時差和橫波時差相關性系數相同且最小為-0.27；根據相關性強度判別標準，水平地應力與垂深之間相關性極強，與井徑、密度、補償中子、自然伽馬、縱波時差、橫波時差之間相關性相對較弱。整體來看，X-1井各測井參數與水平地應力之間均存在一定程度的相關性。

圖3 X-1井測井參數相關性分析圖Fig. 3 Correlation analysis of logging parameters in X-1 well

原始測井數據受外部環境以及測井儀器自身性能的影響，使得原始測井數據具有較大噪聲[35]?？紤]到卡爾曼濾波具有無偏、穩定和概率意義上最優的特點[36]，采用卡爾曼濾波算法對原始測井數據進行去噪。常規卡爾曼濾波算法包含估計和修正兩個過程，估計過程主要依靠時間更新方程，即首先根據式(10)建立當前狀態的先驗估計值êk|k-1，然后根據式(11)推算當前狀態的誤差協方差估計值Pk|k-1：

式中：êk|k-1為先驗估計狀態向量；A為狀態轉移矩陣；Pk為后驗估計誤差協方差矩陣；Pk|k-1為先驗估計誤差協方差矩陣；Q為過程噪聲協方差矩陣。

然后，在獲得先驗估計的前提下，利用觀測值對當前狀態進行后驗估計：首先，根據式(12)計算當前狀態的卡爾曼增益矩陣Kk；然后，結合觀測向量zk，根據式(13)對先驗估計狀態向量êk|k-1進行修正，得到狀態估計向量êk，即為最優估計；最后，根據式(14)將誤差協方差Pk|k-1更新為Pk：

式中：êk為后驗狀態估計向量；zk為觀測向量；R為測量噪聲協方差矩陣；S為測量協方差矩陣；Kk為卡爾曼增益。

采用卡爾曼濾波對X-1井和X-2井的井徑、密度、自然伽馬、縱波時差、橫波時差以及最大和最小水平地應力等測井曲線進行了去噪處理，去噪前后的對比曲線如圖4所示。不難看出：卡爾曼濾波可以有效消除測井數據中儀器自身噪聲和隨機性噪聲信息，有助于避免后續訓練過程中可能會出現的噪聲擬合，從而避免噪聲引起的過擬合或欠擬合。因此，卡爾曼濾波去噪也是一種正則化手段。

圖4 測井參數去噪前后對比Fig. 4 Comparison of logging parameters before and after denoising

完成數據去噪后，需要對測井數據進行歸一化處理。歸一化處理的目的是將不同類型數據映射到一定范圍內，在不影響數據特性前提下，消除量綱對數據的影響，同時，加快機器學習模型的處理過程，從而提高訓練速度。本文采用最大最小歸一化，將測井數據歸一化到0和1之間，最大最小歸一化表示為[37]：

式中：x*為歸一化之后的數據；x為原始數據；xmax和xmin分別為數據的最大值和最小值。

所有測井參數在輸入模型前都要進行歸一化處理，值得注意的是最大和最小水平地應力作為輸出數據，即目標數據不需要進行歸一化，但為了使模型能夠快速收斂，將最大和最小水平地應力值縮小100倍。

2.2 樣本構造

為有效提取測井數據序列的深度特征信息，采用了一種滑動窗口的樣本構造方式[38]。如圖5所示，將歸一化的測井數據組成二維數據矩陣形式，其中，特征表示相關測井參數，標簽表示測井數據點所對應的地應力，n和l分別表示滑動窗口的步長和滑窗的大小，每一個樣本由輸入和輸出樣本組成。在樣本構造過程中，將左側彩色實線框中的數據作為輸入樣本，右側彩色實線框中的數據作為輸出樣本?；懊看窝刂疃确较蚧瑒右粋€步長，便可構造一個樣本，滑窗滑至最深處，即可完成所有樣本的構造?？紤]到本研究為一個深度序列預測任務，為了保證預測參數在深度上連續，取滑窗步長為1。

圖5 樣本構造方法Fig. 5 Sample construction method

3 雙向長短期記憶神經網絡模型

基于測井參數預測最大和最小水平地應力，可描述為已知某一井段測井參數X={x1, x2,…, xn}以及對應的地應力參數Y={y1, y2,…, yn}，其中，n表示樣本數量，通過構建模型，由訓練得到測井參數X與地應力參數Y之間的非線性映射關系。構建模型的目的在于，當給出新的且未進行地應力測井解釋的測井段時，模型能夠預測出最優地應力預測結果Y?={?1, ?2,…, ?n}。

3.1 模型理論

根據上述水平地應力預測問題的描述，本文使用一種特殊的循環神經網絡，即雙向長短期記憶神經網絡(BiLSTM)[39]，BiLSTM網絡結構如圖6所示。一個完整的BiLSTM網絡包括數據輸入層、前向LSTM、后向LSTM、激活函數層及數據輸出層。

圖6 BiLSTM網絡結構圖Fig. 6 Network structure diagram of BiLSTM

BiLSTM主要由前向和后向長短期記憶神經網絡(LSTM)組成[40]，通過兩個平行的隱藏單元層來捕捉前后兩個方向上的序列關系，進一步提高了特征提取的完整性和全局性，使得地應力預測更加準確。LSTM是BiLSTM的基本模塊，不同于普通的循環神經網絡(RNN)，LSTM由于特殊的結構設計使得其能夠有效避免時間序列模型在反向傳播過程中因為對序列產生長期依賴所導致的梯度消失和梯度爆炸問題[41]。相較于傳統的RNN，LSTM結構更為復雜，其包含了四個交互作用門層(遺忘門層、Tanh門層、輸入門層、輸出門層)[42-43]，從而形成了獨有的反饋循環結構，LSTM的網絡結構如圖7所示。

圖7 LSTM網絡結構圖Fig. 7 Network structure diagram of LSTM

LSTM的核心思想是細胞態，其貫穿整個循環過程，在LSTM的處理過程中，通過精心設計的交互作用門層，可根據上一時刻的隱狀態和當前時刻的輸入，選擇性刪除上一時刻細胞態中的信息或者添加新信息到當前時刻的細胞態中。具體而言，第一個交互門層稱為遺忘門層，它決定了哪些信息該從細胞態中刪除，遺忘門層可表示為[39]：

式中：ft為遺忘門；Wf為遺忘門權重矩陣；ht-1為上一時刻的隱狀態；xt為當前時刻的輸入；ht為當前時刻的隱狀態；bf為遺忘門的偏置向量；σ為Sigmoid激活函數。

第二個交互作用門層稱為輸入門層，它決定了哪些信息將會被細胞態存貯，輸入門層可表示為式(17)。接下來便利用第三個交互門層，即Tanh門層生成一個新的候選值，用于描述當前時刻的細胞態，Tanh門層可表示為式(18)[40]：

式中：it為輸入門；t為新的候選值；Wi為輸入門權重矩陣；Wc為Tanh門權重矩陣；bi為輸入門偏置向量；bc為Tanh門偏置向量。

當經過上述步驟后，舊的細胞態Ct-1將被更新成為新的細胞態Ct。首先，將遺忘門的值ft與舊的細胞態Ct-1進行點乘以忘記部分信息，保留更為重要的信息；然后，將輸入門的值it與新的候選值t進行點乘以決定哪些輸入信息將會被保留；最后，將兩者的計算結果相加完成對細胞態的更新；該過程可表示為[41]：

式中：Ct為當前時刻的細胞態；Ct-1為上一時刻的細胞態。

第四個交互門層稱為輸出門層，決定細胞態的哪些信息將會被輸出。在這一階段，首先會將細胞態的值通過Tanh函數映射到-1和1之間，然后與輸出門層的計算結果進行點乘，輸出當前時刻的隱狀態，該過程可表示為：

式中：Ot為輸出門；Wo為輸出門權重矩陣；bo為輸出門偏置向量；ht為當前時刻隱狀態。

對于BiLSTM模型，其隱狀態包含兩部分，一部分來自于前向LSTM，一部分來自于后向LSTM，分別用和表示。因此，BiLSTM的隱狀態可以被定義為[41]：

式中：yt為BiLSTM的輸出；為當前時刻前向LSTM的隱狀態；為當前時刻后向LSTM的隱狀態；c1和c2為常數，且c1+c2=1。

3.2 模型構建

由BiLSTM構建的預測模型具有良好的適應性，在解決復雜地質問題方面表現出了卓越的預測能力[40]。本文所構建的地應力預測模型包括2個BiLSTM網絡層、1個丟棄層(Dropout)、1個拉直層(Flatten)和1個全連接層(Linear)，并分別在BiLSTM網絡層和Linear層后面加上雙曲正切(Tanh)激活函數層和Sigmoid激活函數層，模型結構和訓練流程如圖8所示。其中，BiLSTM網絡層是模型的核心層，用于捕捉測井參數的序列特征以及不同測井參數之間的內在聯系；Dropout層用于避免模型過擬合；Flatten層用于將高維數據特征變換成1維，實現BiLSTM網絡層到全連接層的過渡；Linear層則將所學習的特征映射到樣本標記空間形成預測值；Tanh和Sigmoid激活函數則完成數據的非線性變換。需要注意的是，在模型構建初期，暫時按照經驗將樣本構造中的窗口大小設置為30，數據批尺寸設置為16，學習率設置為0.01，迭代次數為50次；模型的訓練流程可分為以下幾個階段：

圖8 模型結構和訓練流程圖Fig. 8 Model structure and training flowchart

1) 將測井數據進行去噪和歸一化處理后按照2.2節所述樣本構造方式構造出特征和標簽樣本，并將樣本進行打亂，然后分批次喂入模型；

2) 樣本數據進入模型后首先經過BiLSTM網絡，隨著數據在網絡中的自由流動，預測地應力的關鍵特征被提取出來，在這一階段BiLSTM網絡的隱狀態設置為12維，并使用Tanh函數作為非線性激活函數；

3) 從BiLSTM網絡層提取的特征被送入Dropout層，Dropout層按照一定的概率對部分隱藏神經元進行丟棄操作，丟棄率設置為20%；

4) 經過Dropout層后的特征進入Flatten層，這時多維特征變被拉直成1維，使得特征能過渡到Linear層，最后特征經過Linear層，使用Sigmoid函數進行激活后映射成為預測值；

5) 計算預測值與標簽之間的誤差，并選擇自適應矩估計(Adam)算法來調整模型的權值參數，使得訓練過程中的誤差可以不斷減小直到達到迭代終止條件，達到迭代終止條件后保存模型直接用于后續測試集的預測。

4 地應力預測結果及討論

4.1 結果分析

將X-1井常規測井數據和測井解釋得到的水平地應力數據作為訓練集，在訓練過程中隨機選取訓練集30%的數據作為驗證集，將X-2井常規測井參數和測井解釋得到的水平地應力數據作為測試集，以驗證該方法的泛化能力。其中，訓練集和驗證集用于建立和評估模型，測試集用于測試BiLSTM模型預測精度。為了分析不同測井參數組合對水平地應力預測的影響，根據相關性分析結果，按照相關性系數的大小順序，設計了7種不同的測井參數組合方式：組合1僅使用1種測井參數：TVD；組合2使用2種測井參數：TVD和CAL；組合3使用3種測井參數：TVD、CAL和DEN；組合4使用4種測井參數：TVD、CAL、DEN和CNL；組合5使用5種測井參數：TVD、CAL、DEN、CNL和GR；組合6使用6種測井參數：TVD、CAL、DEN、CNL、GR和DTC；組合7使用7種測井參數：TVD、CAL、DEN、CNL、GR、DTC和DTS。

為了評價模型預測效果，采用均方根誤差(RMSE)、平均絕對誤差(MAE)、平均絕對百分比誤差(MAPE)作為評價指標，其計算公式如下[33,38]：

式中：N為樣本數量；yi為真實地應力值，MPa；?i為預測地應力值，MPa；RMSE為均方根誤差，MPa；MAE為平均絕對誤差，MPa；MAPE為平均絕對百分比誤差，%。

利用X-1井訓練機器學習模型，然后按照7種測井參數組合方式，分別對X-2井3218.1～3693.0 m井段地層的最大和最小水平地應力進行預測，同時進行了超參數優化，優化前后預測結果如圖9所示。圖9a為實際和預測最大水平地應力對比，圖9b為實際和預測最小水平地應力對比，其中，紅色實線代表實際水平地應力，綠色實線代表預測得到的水平地應力，藍色實線代表超參數優化后得到的水平地應力。

不同測井參數組合下預測結果的評價指標如表1所示，結合圖9分析不難發現：1) X-2井預測井段實際最大水平地應力介于78.5～98.3 MPa，平均值為84.4 MPa；最小水平地應力介于74.1～93.7 MPa，平均值為79.7 MPa。2) 7種不同測井參數組合模式下，組合6的預測結果與實際水平地應力最為接近，最大和最小水平地應力的平均絕對百分比誤差(MAPE)分別為1.39‰和1.76‰；組合2的預測效果最差，最大和最小水平地應力的平均絕對百分比誤差(MAPE)分別為4.26%和4.33%。3) 組合4～7的預測效果明顯優于組合1～3；從組合3到組合4，在增加了補償中子作為輸入參數之后，預測誤差大幅下降，最大和最小水平地應力的平均絕對百分比誤差(MAPE)分別由38.48‰和40.93‰下降至10.37‰和11.35‰，說明在增加了補償中子后，模型學習到了預測地應力的關鍵信息，分配給了補償中子更大的權重；另一方面，補償中子用于指示巖石的孔隙度，對地層孔隙壓力極為敏感，說明水平地應力與孔隙壓力關系較為密切。

圖9 不同參數組合下超參數優化前后X-2井最大和最小水平地應力預測結果對比Fig. 9 Comparison of prediction results of maximum and minimum horizontal in-situ stress in X-2 well before and after hyperparameter optimization under different parameter combinations

表1 不同參數組合下X-2井水平地應力預測效果評價指標Table 1 Evaluation indices of horizontal in-situ stresses prediction for X-2 well under different parameter combinations

由于機器學習模型超參數的選擇與估計對模型的預測與泛化性能具有重要影響[45]，很多研究學者在模型超參數優化過程中將模型學習率、批尺寸、迭代次數等超參數作為主要影響因素，利用網格搜索、隨機搜索和貝葉斯優化等方法對模型進行優化[46]。為此，將模型學習率、批尺寸、迭代次數三種超參數作為主要影響因素來研究其對模型性能的影響，此外，考慮到樣本構造方式對的特殊性，將滑窗的窗口大小也作為主要影響因素進行研究。針對多因素的交互影響，普通實驗設計方案需要大量重復實驗且耗時長，為此，采用正交設計實驗方法，該方法適合于多因素、多水平的研究，能夠從眾多實驗數據中選取適當具有代表性的數據進行測試，這些實驗數據具有分布均勻、齊整可比的特點[44]。此處以學習率(A)、批尺寸(B)、迭代次數(C)以及窗口大小(D)等四個超參數為因子，選取各因子的典型值作為水平值，每個因子選取四個水平值，其中，學習率分為0.005、0.01、0.015、0.02四個水平值，批尺寸分為16、32、64、128四個水平值，迭代次數分為50、100、150、200四個水平值，窗口大小分為10、20、30、40四個水平值，從而形成了四因子四水平的16組正交實驗方案，如表2所示。此外，正交實驗采用的輸入特征參數組合為組合6，計算了不同實驗方案預測得到的最大和最小水平地應力的評價指標RMSE、MAE和MAPE，結果如表2和圖10所示。

表2 超參數正交設計實驗表及其對應的誤差結果Table 2 Orthogonal design experiment table of hyperparameters and the corresponding errors

圖10 不同實驗方案預測結果的誤差對比Fig. 10 Comparison of prediction errors in different experimental schemes

結合表2和圖10不難看出：1) 不同超參數條件預測的最大和最小水平地應力誤差總體相似，且最小水平地應力的預測誤差略大于最大水平地應力；2) 除了第4組實驗外，其余實驗方案均能夠將水平地應力預測誤差控制在較小的范圍內，其中，最大水平地應力預測結果的RMSE介于0.06～0.16 MPa、MAE介于0.04～0.13 MPa、MAPE介于0.48～1.49‰，最小水平地應力預測結果的RMSE介于0.06～0.20 MPa、MAE介于0.04～0.16 MPa、MAPE介于0.50～1.86‰，說明超參數對水平地應力的影響相對較??；3) 第13組實驗預測得到的水平地應力誤差最小，其中，最大水平地應力預測結果的RMSE、MAE、MAPE分別為0.06 MPa、0.04 MPa和0.48‰，最小水平地應力預測結果的RMSE、MAE、MAPE分 別 為0.06 MPa、0.04 MPa和0.50‰，因此，推薦采用學習率、批尺寸、窗口大小和迭代次數分別為0.02、16、20和150。

根據正交設計實驗計算結果，分析了不同超參數對水平地應力預測誤差的影響規律，結果如圖10子圖所示，不難看出：1) 隨著學習率的增加，MAPE逐漸降低，當學習率為0.02時候，MAPE小幅增加；2) 隨著批尺寸的增加，MAPE逐漸增加，與之相反，隨著迭代次數的增加，MAPE逐漸下降；3) 窗口大小對MAPE的影響無明確規律，在模擬參數條件下窗口大小10或30的MAPE更低；4) 不同超參數影響下最大和最小水平地應力預測結果的MAPE變化趨勢基本一致，但最小水平地應力預測結果的MAPE略高于最大水平地應力。

此外，按照方案13的超參數組合，對表1所述的7種輸入特征參數組合的水平地應力進行了預測，預測結果如圖9所示。圖9a為參數優化前后最大水平地應力預測結果與測井解釋結果的對比，圖9b為參數優化前后最小水平地應力預測結果與測井解釋結果的對比，其中，紅色實線代表實際水平地應力，綠色實線代表預測得到的水平地應力，藍色實線代表超參數優化后得到的水平地應力，通過分析不難發現：1) 按照第13組超參數對模型進行優化后，不同輸入特征組合方式預測得到的最大和最小水平地應力的精度均得到了不同程度的提升，尤其是輸入特征參數組合1～4提升最為明顯，說明第13組超參數優化效果好且非常穩定；2) 超參數優化后，輸入參數組合6預測的水平地應力誤差最小，即輸入參數組合6的預測結果仍然是最優的，而且，優化前后最大水平地應力預測結果的RMSE、MAE、MAPE分別降低了0.13 MPa、0.08 MPa和0.91‰，最小水平地應力預測結果的RMSE、MAE、MAPE分別降低了0.23 MPa、0.11 MPa和1.26‰，說明最小水平地應力的預測精度提升更加明顯；3) 不同輸入參數組合方式下，在增加或減少某一測井參數后的預測效果與超參數優化前的結論基本一致，也進一步說明輸入參數組合6是最優的。

4.2 討論

首先，分析了7種不同測井參數組合對預測水平地應力的影響。對于組合1，只加入了垂深作為輸入參數，其預測結果表現為1條斜線，這是可以預見的，因為垂深本身呈線性，所以模型也只學到了這種線性關系；對于組合2，在增加了井徑作為輸入參數之后預測結果反而變差，常規測井解釋中井徑用于指示地層巖性和判別地應力方向，但對于最大和最小水平地應力的大小而言并無直接關系，因此，井徑的增加給予了模型錯誤的權重配比，反而使預測效果變差；對于組合3，在增加了密度測井曲線作為輸入參數之后，模型預測效果好轉，因為密度測井曲線與最大和最小水平地應力的計算有直接關系，所以增加密度測井曲線有利于模型進行預測；對于組合5，在增加了自然伽馬作為輸出參數之后，預測誤差進一步減小，說明自然伽馬對水平地應力的預測也起到了一定作用，通常自然伽馬用于識別地層巖性和確定泥質含量，顯然也構成了計算地應力的積極要素；對于組合6，在增加了縱波時差之后，預測效果有了進一步提升，事實上縱波時差與地層構造應力之間有直接關系，通常當地層受到額外的構造應力時，縱波時差會偏離正常趨勢線，因此，有學者利用這種關系建立起了預測水平地應力的經驗關系式[47]；對于組合7，當增加橫波時差之后誤差反而增加，可能是因為橫波時差主要反映地層的巖性、孔隙度特征，與補償中子、自然伽馬、縱波時差等測井參數互為同質參數，說明過多的同質參數會使數據維度增加，一定程度上降低模型的學習效果，從而增加預測誤差。綜上所述，結合測井參數實際含義，當增加或減少某些測井參數后，由于不同參數對地應力的指示作用不同，使模型在學習過程中發生權重偏移，進而對不同參數組合模式下的預測結果表現出差異性。因此，綜合考慮測井參數的實際含義和參數獲取的難易程度，推薦使用垂深、井徑、密度、補償中子、自然伽馬和縱波時差6種測井參數作為水平地應力的輸入參數組合。

其次，采用正交設計實驗方法設計了16組實驗，研究了不同超參數對BiLSTM模型預測性能的影響。結果發現在輸入參數組合6條件下16組超參數組合的精度都比較高，一方面是因為輸入參數組合6為最優輸入組合，可以確保模型訓練效果和預測精度；另一方面，超參數取值范圍通過經驗設置在常用范圍內，避免了因為梯度爆炸或者梯度消失而引起的精度降低。將輸入參數組合6條件下優化獲得得超參數用于其它6種輸入參數組合，結果發現不同輸入特征組合方式預測得到的最大和最小水平地應力的精度均得到了不同程度的提升，而且輸入參數組合6的預測結果仍然是最優的。從預測結果來看，適當減小批尺寸和增加迭代次數可以增加模型預測精度，但減小批尺寸將導致模型學習效率的下降，增加迭代次數不僅增加訓練成本而且增大了過擬合的可能性，因此，需要針對不同數據集，按照實際情況進行合理調參。對于本研究推薦學習率、批尺寸、窗口大小和迭代次數等超參數分別取值為0.02、16、20和150。

訓練數據集數量限制，數據集樣本獲取途徑單一，都是造成預測精度低的重要原因。考慮到工程中實際條件的限制，通常無法獲取足夠的測井數據集進行訓練，因此，采用滑窗采樣的樣本構造方式和特征組合篩選實驗來優化輸入特征，并引入BiLSTM這一特殊循環神經網絡來優化模型結構，另外設計了正交實驗對超參數取值進行優化，以提高模型的預測精度和泛化能力。但是，本研究尚存在一定的局限性：首先，本研究僅對CL氣田兩口直井的測井數據進行測井解釋、訓練和驗證，所建立的模型數據來源較為單一，因此有必要加入地震、巖心實驗、礦場實測等多源、多態、多尺度的數據擴充訓練樣本集，從而構建出普適性和穩健性更好的地應力預測模型。其次，神經網絡模型是基于經驗風險最小化原則，在小樣本預測問題中存在過度擬合的風險，因此在后續建模過程中可以嘗試以結構風險為原則或以結構和經驗風險相結合為原則的機器學習模型，強化模型在理論方面的引導，避免因純數據驅動造成的模型可解釋性和泛化能力降低的問題。

5 結論

1) X-1和X-2井地應力測井解釋結果表明，地應力排序為垂向地應力＞最大水平地應力＞最小水平地應力，屬于潛在正斷層應力狀態；對比測井解釋和巖心差應變測試結果發現，垂向地應力解釋誤差為0.39%，最大水平地應力解釋誤差為0.18%～0.64%，最小水平地應力解釋誤差為0.29%，說明測井解釋得到的地應力與實際地應力吻合較好。

2) 水平地應力與垂深、密度及自然伽馬等測井參數之間表現為正相關，與縱波時差、井徑、補償中子及橫波時差等測井參數之間表現為負相關；水平地應力與垂深之間相關性極強，與井徑、密度、補償中子、自然伽馬、縱波時差、橫波時差之間相關性相對較弱。

3) 不同測井參數蘊含了不同的地質信息，不同測井參數對地應力預測起到了不同的指示作用。通過對比7種不同測井參數組合模式發現，組合4～7的預測效果明顯優于組合1～3，最優輸入參數組合為組合6，即包括垂深、井徑、密度、補償中子、自然伽馬、縱波時差；進一步，通過正交設計實驗對超參數取值進行調優，最佳的超參數取值為學習率0.02、批尺寸16、窗口大小20，迭代次數150，預測得到的最大和最小水平地應力平均絕對百分比誤差分別為0.48‰和0.50‰。

4) BiLSTM模型是一種可以實現水平地應力精準預測的機器學習方法，得益于BiLSTM對時間序列出色的處理能力，使得其能夠有效捕捉測井參數隨深度的變化趨勢以及測井參數的前后關聯信息。因此，推薦使用BiLSTM模型進行水平地應力的預測。