基于半監督學習的井震聯合儲層橫向孔隙度預測方法

韓宏偉， 劉浩杰， 桑文鏡， 魏國華， 韓智穎， 袁三一*

1 中國石化集團公司勝利油田物探研究院， 東營 257000 2 中國石油大學(北京)油氣資源與探測國家重點實驗室， 北京 102249

0 引言

孔隙度是描述儲層巖性、物性和含油氣性特征的重要參數之一.準確預測孔隙度有利于開展甜點區域優選、油氣儲量估算、井網井位設計、油藏數值模擬等工作(Angeleri and Carpi，1982；Doyen，1988；Leite and Vidal，2011).根據不同的地球物理探測和觀測技術，孔隙度預測方法主要包括巖石物理實驗測量、測井解釋和井震聯合預測三類.巖石物理實驗測量常采用氦氣法和高壓壓汞法等測得巖心孔隙度(韓學輝等，2021).測井解釋主要采用經驗公式、三孔隙度模型、交會圖、巖石物理模型、神經網絡等方法構建孔隙度與密度、補償中子、聲波時差、自然伽馬、泥質含量、地層電阻率等敏感參數之間的計算模型(Archie，1942；Xu and White，1996；Khaksar and Griffiths，1998；李宏兵等，2019；安鵬等，2019).這兩類方法能實現精確測量巖心段或全井段孔隙度，但無法有效外推獲得井間孔隙度.

井震聯合儲層孔隙度預測的理論基礎是地震波在地層中的傳播速度與巖石孔隙度和泥質含量等有關(Wyllie et al.，1956；Klimentos and McCann，1990).因此，整合地震和測井數據在數據覆蓋范圍、縱橫向分辨率等方面的各自優勢，能較好地預測出儲層孔隙度的空間分布.井震聯合預測方法主要是基于不同的地震屬性直接或間接地預測孔隙度.根據利用的地震屬性和實現過程的差異，先后發展出疊后單屬性預測、疊后多屬性融合、疊前多屬性融合三類孔隙度預測方法.疊后單屬性預測方法一般是基于經驗公式和巖石模型等方式進一步線性或非線性地轉化疊后波阻抗反演結果(或速度及聲波時差)為儲層孔隙度(Angeleri and Carpi，1982；de Buyl et al.，1988).或考慮地震垂直走時和儲層孔隙度的橫向變化的空間相關性，采用協同克里金等地質統計學技術進行儲層建模并生成孔隙度模型(Doyen，1988；Moon et al.，2016).疊后多屬性融合預測方法采用神經網絡、隨機模擬等技術直接映射優選的地震屬性組合為儲層孔隙度.Hampson等(2001)通過多元逐步線性回歸分析優選與預測目標參數相關性最好的地震屬性組合，并結合兩個應用實例說明概率神經網絡相比于多元線性回歸能更好表達地震屬性與孔隙度之間的非線性關系，進而提高預測的可靠性.Leite和Vidal(2011)研究了一種基于神經網絡融合自然伽馬數據體和遞推反演得到的阻抗數據體的三維孔隙度預測方法.

相比于疊后地震數據，疊前地震數據的振幅隨偏移距變化(AVO)特征包含著更為豐富的儲層巖性、物性與流體信息.因此，理論上疊前數據預測儲層物性參數有利于提高預測準確性和可靠性，降低多解性.疊前多屬性融合預測方法依托疊前同步反演和多屬性融合兩大技術，通過疊前反演獲得與孔隙度有較好相關性、包含更多反映儲層物性特征的疊前彈性參數(如縱橫波速度、縱橫波速度比、縱橫波阻抗、密度、楊氏模量、泊松比等).隨后優選并組合疊前AVO屬性(如截距、梯度等)、疊后地震屬性(如振幅、頻率、波形三大類等)和疊前彈性參數中的一類或多類參數，利用神經網絡和隨機模擬等多屬性擬合技術建立基于優選屬性的孔隙度預測模型(AlMuhaidib et al.，2012).韓宏偉等(2021)系統地總結了地震儲層預測技術表征孔隙度等儲層參數的原理與應用案例.

現有的井震聯合儲層孔隙度預測方法要解決的核心問題是如何優選出與孔隙度相關的地震屬性以及提高彈性參數預測精度.人工智能技術可以間接解決上述兩個核心問題，為從地震數據中自動提取、優化和組合地震屬性來直接預測孔隙度提供了新的可行性途徑.近幾年來，深度學習在地球物理勘探領域的應用逐步從地震資料處理、解釋延伸到彈性參數反演、儲層參數預測等方面(Yuan et al.，2018，2022；Sang et al.，2021；Song et al.，2021；王迪等，2021).以卷積神經網絡、遞歸神經網絡和生成對抗網絡為主的不同類型深度神經網絡也開始逐漸應用于儲層孔隙度等物性參數預測任務.目前，這類智能建模方法在多測井曲線擬合井點位置的孔隙度方面研究較多，在井震聯合直接預測儲層橫向孔隙度方面研究相對較少.在預測孔隙度曲線方面，考慮到測井數據是反映不同地質時期地層沉積特征的局部相關的時間序列，眾多學者采用循環神經網絡建立多種敏感的測井參數轉換為孔隙度的非線性映射，以實現孔隙度建模與預測(安鵬等，2019；宋輝等，2019；王俊等，2020；Chen et al.，2020).相比于全連接神經網絡、多元線性回歸分析等方法只考慮來自同一深度點(或時間采樣點)的測井曲線之間的相互聯系，考慮測井曲線漸變沉積規律和前后關聯的循環神經網絡具有更高的預測精度和魯棒性.在井震聯合儲層橫向孔隙度直接預測方面，Das和Mukerji(2020)采用端到端卷積神經網絡進行有監督學習，從疊前地震數據預測出孔隙度.Feng等(2020)整合卷積神經網絡與褶積模型搭建物理約束的混合網絡，模擬地震反演和正演過程開展無監督孔隙度預測.該方法不需要大量的標簽數據，適合緩解地震測井數據不均衡引起的油氣小樣本問題.但是該方法依賴于準確的低頻信息和震源子波信息參與混合網絡的訓練.

為減少低頻模型對神經網絡預測孔隙度的影響，并緩解基于有監督學習范式下的智能化孔隙度預測結果橫向連續性差的問題，本文提出一種基于半監督學習的井震聯合儲層橫向孔隙度預測方法.該方法基于雙向門控遞歸單元(Bi-GRUs)構建由編碼網絡和解碼網絡組成的孔隙度預測混合網絡(Cho et al.，2014).其中，編碼網絡模擬地震反演及巖石物理建模過程，轉化輸入疊后地震數據為孔隙度；解碼網絡建立數據驅動的正演模型，實現孔隙度非線性映射為生成地震數據.此外，非井旁地震道在每次迭代過程中隨機引入混合網絡參與網絡訓練.通過整合井位置的地震數據和(偽)測井數據及非井位置的地震數據開展半監督學習，有望既保證預測孔隙度的精度，同時通過地震波形匹配間接保證孔隙度結果的橫向連續性，獲得高精度的儲層孔隙度模型.

1 方法原理

本小節首先回顧基于有監督學習構建孔隙度預測網絡的基本思想.之后介紹本文提出的半監督井震聯合孔隙度預測網絡的主要原理，包括網絡架構、目標函數和實現過程等.最后介紹合成數據集的生成過程，該數據集用于后續測試的網絡智能建模與結果定量評價.

1.1 有監督孔隙度預測網絡

隨著以深度學習為代表的人工智能技術逐漸用于求解高分辨率處理、波阻抗反演、層析反演、全波形反演、彈性參數反演等反問題，儲層物性參數智能反演也開始受到重視(Araya-Polo et al.，2018；Kim and Nakata，2018；Alfarraj and AlRegib，2018，2019；Biswas et al.，2019；Gu et al.，2019；Russell，2019；Sun and Alkhalifah，2020；Aleardi and Salusti，2021).目前，通過學習由地震數據和目標參數構成的訓練樣本對建立有監督儲層參數預測網絡是一種常用的儲層參數智能反演方法(Alfarraj and AlRegib，2019).在有監督學習框架下，本質上孔隙度等物性參數的智能反演與其他儲層參數反演類似.其實現的數學基礎是理論上神經網絡能直接建立地震數據到任意目標參數的非線性映射關系，通過梯度下降迭代優化目標函數而逼近目標解；實現的物理基礎是不同類型的地球物理數據從不同角度刻畫地下相同目標地質體，它們之間具有內在的物理關聯和相似的響應特征.

地層通常遵從漸變的沉積規律，局部地層性質是近似且有關聯的.因此，連續記錄的地震或測井數據在不同時刻或不同深度獲得的響應與其附近響應存在內在關聯和局部相似，綜合反映了周圍地層的性質.基于以上分析，地震記錄和測井曲線都可以視為內部關聯、局部依賴的序列數據.在語音識別和文本分類等自然語言處理場景中，循環神經網絡(RNN)已被廣泛應用于解決序列建模問題(Hochreiter and Schmidhuber，1997；Chorowski et al.，2014).在地震勘探領域的彈性參數反演和物性參數預測等方面，不少學者通過比較前饋神經網絡和循環神經網絡的預測性能發現，后者更適合準確預測儲層參數的變化趨勢和局部細節(安鵬等，2019；王俊等，2020).相比于前饋神經網絡，RNN能記憶上一時刻隱藏層學到的信息，并利用該歷史信息影響當前時刻隱藏層的輸出.但是，常規RNN由于不同時刻信息的長距離依賴容易出現梯度消失或梯度爆炸的問題, 導致其難以訓練獲得較精準的網絡模型.而雙向門控遞歸單元(Bi-GRUs)通過更新門和重置門靈活決定前一時刻隱層信息的記憶或遺忘，緩解常規RNN模型難以優化更新的問題.此外，其能避免常規RNN信息狀態的單向傳遞，采用兩個方向相反的單向門控遞歸單元(GRU)實現信息的正向和反向傳播，使得當前時刻的輸出由當前時刻之前和之后的隱層信息共同決定.這一工作機制使得其在預測某一時間點或深度點的孔隙度時能考慮局部相鄰的地震或測井模式，獲得更能反映地層變化規律的孔隙度結果.因此，本文采用Bi-GRUs作為構建孔隙度預測網絡的基本單元.以GRU為例，其內部細胞結構是由重置門、更新門和tanh激活函數三者構成的.重置門和更新門的輸出值都與當前時刻的輸入和上一時刻的隱層輸出有關，且更新門輸出值越小，歷史狀態信息對當前時刻GRU輸出結果的影響程度越小.由于Bi-GRUs的工作原理不是本文研究重點，更多的原理介紹參考Cho等(2014)、Alfarraj和AlRegib(2019)和Yuan等(2022).

有監督孔隙度預測網絡(圖1a，以下簡稱為有監督網絡)一般包括輸入層、隱藏層和輸出層三個部分.輸入層和輸出層分別是單道地震記錄和網絡預測的孔隙度曲線.為與后面的半監督孔隙度預測網絡形成對比，中間隱藏層由4個雙向門控遞歸單元和1個全連接層構成.有監督網絡主要通過成對的地震記錄和孔隙度曲線進行學習和訓練，該過程可以視為一種廣義的聯合疊后地震反演與巖石物理模型的建模流程，即模擬地震數據到彈性參數，彈性參數到物性參數的物理反演過程，實現地震數據直接轉化為孔隙度.有監督網絡(簡記為E)通過反向傳播算法不斷更新網絡參數，減小網絡預測結果與孔隙度標簽之間的差異，其差異使用如下目標函數LE計算測井監督損失：

(1)

其中，N為訓練樣本或標簽的數量，Si和Pori分別為第i個地震道曲線、第i個孔隙度曲線，θE為有監督網絡E的網絡參數，E(Si;θE)為網絡E預測的孔隙度結果.

通過迭代最小化目標函數，最終得到的有監督孔隙度預測模型能實現從單道地震記錄預測出最佳的單道孔隙度曲線.但是，由于該逐道反演過程沒有考慮儲層的橫向變化，特別是在預測的單道孔隙度結果出現較大偏差時，預測孔隙度剖面的橫向連續性差.

1.2 半監督孔隙度預測網絡

為緩解常見的有監督學習孔隙度預測方法存在著信息利用不足、預測結果連續性差等問題，本文提出一種基于半監督學習的孔隙度預測方法.該半監督孔隙度預測網絡(圖1b，以下簡稱半監督網絡)是由編碼網絡和解碼網絡組成的混合網絡.編碼網絡與上一小節有監督網絡的網絡結構保持一致(編碼網絡也簡記為E)，也是由4個雙向門控循環單元和1個全連接層構成；而解碼網絡D與編碼網絡E采用相同的網絡模塊.E和D的區別在于各自的輸入、輸出和網絡內部權重與偏置不同，E的輸入是單道地震記錄，輸出是預測的孔隙度曲線；D的輸入是E預測的孔隙度曲線，輸出是生成的地震記錄.以編碼網絡E為例，E中不同隱層的雙向門控遞歸單元通過更新候選隱層狀態和隱層狀態實現從輸入地震數據或上一隱層輸出進一步提取高維的多水平特征，不同水平特征最終通過線性回歸層加權輸出得到期望的孔隙度曲線.預測的孔隙度進一步進入解碼網絡D正演生成地震數據.

圖1 井震聯合儲層孔隙度預測網絡(a) 有監督孔隙度預測網絡； (b) 半監督孔隙度預測網絡.Fig.1 Reservoir porosity prediction network via seismic and well logs integration(a) Supervised porosity prediction network； (b) Semi-supervised porosity prediction network.

實際工區稀疏分布的測井數據和稠密采集的地震數據導致智能化儲層參數預測面臨著數據不均衡引起的油氣小樣本問題.即稀疏的井點位置可以生成地震測井標簽樣本對；而眾多的非井點位置只有無標簽的地震數據.為盡可能充分利用地震和測井信息保證井間孔隙度預測效果，本文提出的半監督網絡除使用有監督網絡中有限的成對井旁地震道和測井標簽外，還采用眾多的非井旁地震道共同構建孔隙度預測網絡的訓練集，即聯合少量的有標簽數據和眾多的無標簽數據開展半監督學習預測孔隙度.半監督網絡在每次迭代優化的過程中，首先利用少量的有標簽數據融合窄帶地震數據和寬帶測井數據這兩種不同尺度信息，同時建立地震數據到孔隙度的編碼網絡和孔隙度到生成地震數據的解碼網絡，并計算井旁地震道數據匹配損失和測井監督損失.接著進一步在當前迭代周期內隨機地從地震剖面上選擇選取N個非井位置的地震記錄作為無標簽數據依次通過編碼網絡和解碼網絡而參與網絡訓練與優化.非井旁地震道由于缺少真實測井標簽而無法采用測井監督損失評估其對應的孔隙度預測效果, 因此只采用非井旁地震數據匹配損失間接保證非井位置的孔隙度預測結果準確而且穩定，不易受到地震振幅特征變化而出現較大偏差.井旁地震道數據匹配損失、測井監督損失和非井旁地震數據匹配損失三項構成半監督混合網絡的總損失.每次迭代過程計算出總損失后再通過反向傳播實現梯度更新并進入到下一次迭代，當訓練集和驗證集損失都收斂到極小值后，完成半監督網絡的訓練.由于半監督網絡中編碼網絡與有監督網絡結構相同，其目標函數也都是公式(1)；解碼網絡的目標函數LD為：

(2)

其中，S′i為第i個非井位置的地震記錄曲線，θD為解碼網絡D的網絡參數，D(E(Si;θE);θD)和D(E(S′i;θE);θD)分別為Si、S′i通過混合網絡后生成的地震數據.解碼網絡目標函數中的第一、二項分別用于計算井旁地震道和非井旁地震道的地震數據匹配損失.此外，不同迭代次數下的非井旁地震道是從原始地震數據中隨機選擇的，以保證網絡較好地適應不同波形特征預測出精確的孔隙度.混合網絡的總目標函數Ltotal為：

Ltotal=LE+λLD，

(3)

其中，λ為正則化參數，用于控制地震數據匹配程度對孔隙度預測效果的影響.當λ=0，公式(3)與公式(1)等價，即有監督學習范式是半監督學習范式的一種特例.本文提出的半監督孔隙度預測方法擯棄了彈性參數反演、屬性優選等復雜過程，編碼網絡和解碼網絡在地震數據匹配損失和孔隙度估計損失的聯合監督下同時迭代優化網絡參數，實現從地震數據直接定量預測地下儲層孔隙度參數的空間分布，而且保證預測孔隙度能通過解碼網絡“正演”生成逼近真實的地震數據.

1.3 合成數據集生成

為證明本文提出的方法預測孔隙度的有效性，選取常用的Marmousi模型(Martin et al.，2006)進行測試.圖2a為通過速度模型導出的孔隙度模型，其模擬地層的巖性以泥巖為主，還含有部分砂巖、泥灰巖和鹽巖.模型中部從淺到深依次經過三條大斷裂、一個泥灰巖背斜、一個不整合面、一個楔形真空鹽巖層和一個泥巖背斜地層.圖2b為阻抗導出的反射系數與震源子波正演生成的無噪合成地震數據.震源子波的主頻為35 Hz，采樣間隔為1 ms.圖2c為在圖2b基礎上添加隨機噪聲得到的含噪合成地震數據，噪信比為0.3(噪聲與地震信號的能量之比).

圖2 用于孔隙度預測的合成數據集(a) 孔隙度模型； (b) 無噪合成地震數據； (c) 含噪合成地震數據.Fig.2 The synthetic dataset for porosity prediction(a) Porosity model； (b) Noise-free synthetic seismic data； (c) Noisy synthetic seismic data.

選取孔隙度模型和無噪合成地震剖面中的部分數據構建訓練集和驗證集，無噪和含噪合成地震剖面作為測試數據，孔隙度模型作為參考標準評價有監督和半監督孔隙度預測方法的優劣.評價指標采用相對誤差(R)和相關系數(CC)兩種.相對誤差定義為：

(4)

其中，X和X′分別為真實孔隙度模型和預測孔隙度結果，m和n分別為孔隙度模型或預測孔隙度結果的時間采樣點數和道數.相關系數定義為：

(5)

其中，X為真實孔隙度模型(或觀測地震數據)，X′為預測孔隙度結果(或生成地震數據)，Cov(X,X′)、Var[X]和Var[X′]分別為X與X′之間的協方差、X的方差和X′的方差.

2 例子

2.1 數值模型測試

本小節采用Marmousi模型進行方法測試與評價.訓練有監督網絡時，訓練集由隨機抽取的8條地震記錄及其相應的孔隙度曲線(圖2b和圖2a中的黑線)構成，驗證集由等間隔抽取的6條地震記錄及其對應的孔隙度曲線(圖2b和圖2a中的紅線)組成，測試集由無噪地震數據(圖2b)、含噪地震數據(圖2c)和真實孔隙度模型(圖2a)形成.三個數據集中每條地震記錄或孔隙度曲線的大小都為750×1.訓練半監督網絡時，驗證集和測試集與有監督網絡使用的保持一致，訓練集為成對的地震記錄和孔隙度曲線(圖2b和圖2a中的黑線，即有監督網絡的訓練集)及非井旁地震道.在半監督網絡的不同迭代訓練周期內，其訓練集中的非井旁地震道是利用不同種子數隨機生成的，而有標簽數據對始終維持不變，即用于訓練半監督網絡的無標簽數據遠多于有標簽數據.此外，不同數據集中的地震數據和孔隙度都采用Z-score歸一化進行預處理：

(6)

在訓練數據、優化器、目標函數和反向傳播算法等的共同作用下，有監督網絡通過預先設定最大迭代次數完成初步學習與訓練，通過繪制損失曲線觀察訓練損失和驗證損失是否同時收斂到極小值，若沒有則調整最大迭代次數直至二者都滿足收斂條件.最終確定最大迭代次數設置為300次時可以獲得最優的有監督孔隙度預測模型，對應的損失曲線見圖3.之后推廣到歸一化后的無噪地震數據進行測試，網絡的輸出結果采用公式(6)中的均值和方差反歸一化處理化后得到最終預測的孔隙度(圖4a).孔隙度預測結果與真實孔隙度模型(圖2a)的相關系數和相對誤差分別為0.9616和2.59%，表明預測結果整體精度較高.網絡單道反演預測得到的孔隙度在淺部橫向連續性較好，這是因為孔隙度模型在淺部(小于0.25 s的部分)地層起伏較小，橫向變化較為緩慢，且淺部傾斜地層、褶皺和斷層內的孔隙度整體變化也較小.而孔隙度模型在大于0.25 s的中深部地層起伏較大，發育有陡傾角的斷層和背斜等，且孔隙度在中深度部的局部差異較大.導致有監督網絡在背斜等復雜構造位置推廣能力相對較差，因而預測結果的橫向連續性會變差，局部出現明顯的預測偏差.其原因可能是對應位置的低頻成分難以準確估計.可以考慮學習相對孔隙度、加入孔隙度初始模型參與網絡訓練、多道地震數據預測單道孔隙度等方式緩解該問題.此外，若采用不同隨機種子為有監督網絡設置不同的網絡參數初始化狀態，以此建立的多個有監督網絡預測出多個孔隙度模型，取不同結果的均值作為統計性孔隙度預測結果也可以緩解該問題.但是，該方式的計算成本較高.采用傳統模型驅動類方法可以獲得速度或阻抗參數，再根據孔隙度與速度或阻抗的巖石物理關系可以進一步獲得預測的孔隙度結果.參考She等(2018)的工作，圖4b為采用基于全變分約束阻抗反演方法導出的孔隙度結果.從定量預測角度來看，有監督方法預測出更符合真實模型的孔隙度結果，其橫向連續性較傳統方法預測效果要差一些.其原因是傳統方法使用了初始模型提供測井先驗信息.但是，該傳統方法由彈性參數到孔隙度的二次轉換及其方法本身對地下地質結構的塊狀假設導致預測的孔隙度整體出現較大偏差，表現為預測孔隙度值整體偏小.

圖3 有監督孔隙度預測網絡的損失曲線Fig.3 Loss curves of the supervised porosity prediction network

圖4 智能方法與傳統方法孔隙度預測結果對比(a) 基于有監督孔隙度預測網絡得到的孔隙度結果(迭代300次)； (b) 基于疊后波阻抗反演方法導出的孔隙度結果.Fig.4 Comparison of porosity prediction results between the intelligent method and the traditional method(a) Porosity prediction result via supervised porosity prediction network (300 iterations)； (b) Porosity prediction result derived from the post-stack acoustic impedance inversion method.

在與有監督網絡的最大迭代次數保持一致的情況下，使用提前準備的數據集訓練半監督網絡，其使用(偽)測井數據和(偽)井位置及非(偽)井位置的地震數據匹配共同監督和質控網絡的訓練過程.目標函數(式(3))中地震數據匹配項的正則化參數λ會影響半監督網絡的迭代優化方向和孔隙度預測效果.下面對比分析不同正則化參數情況下訓練得到的半監督網絡推廣到無噪地震數據(圖2b)預測的孔隙度結果與有監督網絡預測孔隙度結果之間的差異，并討論λ大小對孔隙度預測及地震數據匹配效果的影響.當λ設置為0時，半監督網絡預測的孔隙度結果(圖5a)與有監督網絡預測的孔隙度結果(圖4a)完全一致，表明此時兩種網絡建立的孔隙度預測模型是基本相同的.由于此時的地震數據匹配項完全沒有起到任何作用，所以預測孔隙度結果(圖5a)的精度及橫向連續性并沒有得到改善，并且其進一步通過半監督網絡中的解碼部分生成的地震數據(圖5f)與觀測地震數據(圖2b)無法匹配，二者的殘差(圖5k)極大.當λ設置較小時(λ=0.1)，半監督網絡預測得到的孔隙度結果(圖5b)相較于圖3的橫向連續性得到一定改善，但是地震數據匹配損失占總損失的比例較小，導致半監督網絡生成的地震數據(圖5g)與觀測地震數據(圖2b)之間的殘差仍然較大(圖5l).當λ設置較大時(λ=1.5)，地震數據匹配過多地影響孔隙度預測效果，使得孔隙度結果(圖5c)中出現了與觀測地震數據(圖2b)形態相似的預測假象，并且在淺部(0～0.1 s)表現尤為明顯.此時，生成地震數據(圖5h)相比于圖5g進一步接近觀測地震數據.當λ設置合理時(λ=0.37)，在地震波形控制和測井監督的共同作用下，不僅可以獲得比有監督方法預測精度更高且橫向連續性更好的孔隙度預測結果(圖5d)，而且預測孔隙度通過數據驅動的“正演模型”(即解碼網絡)可以進一步實現準確的地震數據匹配(圖5i).以上4組對照測試與有監督網絡設置的迭代次數一致，都為300次.在確定出合適的正則化參數后，繼續增加迭代次數(如增大到1000次)可以進一步提高孔隙度預測結果(圖5e)的質量，同時進一步減小生成地震數據(圖5j)與觀測地震數據之間的差異(圖5o).對比圖5a、b及圖4a可以看到，引入地震數據匹配約束會改善孔隙度預測結果的橫向連續性，但需要設置合理的正則化參數才能保證預測準確性和連續性的同時提升.對比圖5f—j及圖5k—o不難發現，地震數據匹配的好壞可以在一定程度上反映預測孔隙度剖面的精度，并且通過生成地震數據也可以間接判斷預測孔隙度剖面的橫向連續性.這是因為當孔隙度預測結果出現較大偏差時也會導致生成的地震記錄不準確，即預測孔隙度剖面與生成地震剖面的橫向連續性具有較強的一致性(如圖5a、f，圖5b、g等).

圖5 基于半監督孔隙度預測網絡測試得到的孔隙度結果(a)—(e) 孔隙度預測結果； (f)—(j) 孔隙度預測結果通過解碼網絡D生成的地震數據； (k)—(o) 生成地震數據(f)—(j)與觀測地震數據(圖2b)之間的殘差. (a)—(e)對應的正則化參數分別為0、0.1、1.5、0.37、0.37.此外，(a)—(d)對應的迭代次數都設置為300次，(e)對應的迭代次數設置為1000次.Fig.5 Porosity results via the semi-supervised porosity prediction network(a)—(e) Porosity prediction results； (f)—(j) Generated seismic data by inputting the predicted porosity of (a)—(e) into the decoder D，respectively； (k)—(o) Data residual between the generated seismic data of (f)—(j) and the observed seismic data of Fig.2b，respectively. The regularization parameter corresponding to (a)—(e) is set to 0，0.1，1.5，0.37 and 0.37， respectively. In addition, the iterations for (a)—(d) are set to 300，and for (e) are set to 1000.

在直觀對比圖4和圖5中的預測孔隙度剖面差異基礎上，進一步對比有監督與半監督兩種網絡在盲井位置預測的單道孔隙度曲線的差異.圖6每個子圖中的黑線、藍線和紅線分別表示從真實孔隙度模型(圖2a)、有監督網絡預測的孔隙度剖面(圖4a)和半監督網絡預測的孔隙度剖面(圖5d)中提取出的來自同一CDP位置的孔隙度曲線.圖6a—c對應的CDP位置分別是16、499和688.在經過楔形體及下伏地層(圖6a的橢圓內)時，孔隙度曲線出現了急劇增大到固定值后又急劇減小的現象.兩種網絡都沒有很好地適應這種縱向跨度比較大的孔隙度突變，但是半監督網絡相比于有監督網絡預測效果略好.在經過斷層及背斜等(圖6b、c橢圓內)區域時，半監督網絡比有監督網絡能更好地適應局部孔隙度的復雜變化，說明前者更適用于實際非均質性強的復雜儲層孔隙度預測.從圖6可以看出，半監督網絡相比于有監督網絡能提高不同位置的單道孔隙度預測精度，最終在孔隙度預測剖面上表現為減小數據驅動的逐道反演引起的橫向不連續性.

圖6 有監督與半監督網絡盲井測試結果對比基于第16道(a)、第499道(b)、第668道(c)地震記錄預測的孔隙度曲線.Fig.6 Comparisons of blind well test via supervised network and semi-supervised networkPredicted porosity curves based on the 16th (a)，499th (b)， and 668th (c) seismic record. Black，blue and red lines represent true supervised network inverted，and semi-supervised network inverted porosity curves，respectively.

實際地震數據在野外采集時往往由于人為或環境因素的干擾而受到噪聲污染，為進一步說明半監督網絡相比于有監督網絡具有更好的抗噪性.采用與圖4和圖5d分別對應的無噪數據訓練得到的有監督網絡和半監督網絡直接推廣到含噪合成地震數據(圖2c，即模擬的實際數據).此時，半監督網絡預測孔隙度剖面(圖7b)比有監督網絡預測孔隙度剖面(圖7a)更為干凈和連續，說明前者比后者更適用于低信噪比地震資料.表1統計了圖4—圖7中所有測試預測孔隙度剖面與真實孔隙度剖面之間的相對誤差和相關系數，及半監督網絡生成地震數據與對應的輸入地震數據之間的相關系數.

圖7 有監督與半監督網絡含噪測試結果對比(迭代300次)(a) 有監督網絡基于含噪地震數據(圖2c)得到的孔隙度預測結果； (b) 半監督網絡基于含噪地震數據(圖2c)得到的孔隙度預測結果.Fig.7 Comparisons of noise-tolerance test results via supervised network and semi-supervised network (300 iterations)(a) Predicted porosity result by inputting noisy seismic data of Fig.2c into the supervised network； (b) Predicted porosity result by inputting noisy seismic data of Fig.2c into the semi-supervised network.

表1 不同方法預測孔隙度結果對比Table 1 Comparisons of porosity prediction results using different methods

Marmousi模型數據測試結果表明：(1)有監督或半監督孔隙度預測方法相比于傳統的聯合地震反演與巖石物理模型的反演方法能避免彈性參數到物性參數的二次轉換引起誤差累積和放大，實現直接從地震數據預測出孔隙度，且預測結果更為準確.(2)相比于有監督網絡，半監督網絡通過非井位置的地震數據匹配間接保證井間孔隙度的預測精度.井間地震數據匹配程度一定程度上可以間接反映孔隙度預測質量.當非井位置的孔隙度預測結果也較好而不出現較大偏差時，其通過解碼網絡才能合成與觀測地震數據接近的生成地震數據.半監督網絡通過學習反映局部地層變化的大量地震波形信息而增強了刻畫薄層孔隙度細節及薄厚層孔隙度突變的能力，提高了孔隙度預測的精度和橫向連續性.(3)在使用的有標簽數據一致的情況下，半監督網絡與有監督網絡的預測效果存在一定關聯.即當公式(3)中的正則化參數越小時，二者預測效果越接近.半監督網絡需要引入合適的地震數據匹配約束數據驅動反演過程，以提高孔隙度預測精度.

2.2 實際數據應用

最后，采用來自中國東部某地區的疊后地震數據對本文提出的半監督孔隙度預測方法進行檢驗.由于工區直井缺少實際解釋的孔隙度曲線，本文首先采用工區內已解釋的斜井孔隙度曲線與不同的孔隙度敏感參數進行交會分析.如圖8所示，在統計工區111口斜井在儲層段(1300～1400 m)的孔隙度與速度關系發現，使用線性回歸即可較好地從速度曲線擬合出孔隙度曲線，擬合結果與已知孔隙度值的相關系數較高，為0.8385.在明確孔隙度曲線生成方法的基礎上，下面針對某連井剖面(圖9a)開展儲層孔隙度預測.該連井剖面有141道，每道地震記錄有111個時間采樣點，道間距為12.5 m，時間采樣間隔為2 ms.連井剖面的時間范圍為1.1～1.32 s，所在儲層埋深較淺，自上而下依次為曲流河沉積和辮狀河沉積，主要的沉積微相包括心灘壩、辮狀河道、河道邊緣和河間灘等.儲層物性較好，平均孔隙度為30%左右，各小層的平均滲透率在(600～1100)×10-3μm2左右.由于河流頻繁改道及多條大斷層的發育使得儲層非均質性強，巖性縱橫向變化大，呈現出多個砂泥巖薄互層縱向交替疊置.在地震剖面上表現為整體較破碎，地震同相軸不連續且變化快.

圖8 儲層段孔隙度與速度的交會圖Fig.8 Cross-plot of porosity versus velocity in the reservoir intervals

如圖9a所示，該連井地震剖面從左到右依次經過6口井(命名為A1—A6).6口井經過井震標定合成的地震數據與對應的井旁地震道匹配程度高，平均相關系數為0.79.利用圖8的巖石物理關系擬合出的6條孔隙度曲線，選取前5條孔隙度曲線及其對應的地震數據訓練有監督網絡和半監督網絡.訓練的迭代次數都設置為3000次，正則化參數設置為1，其他設置與合成數據例子相同.訓練好的有監督和半監督孔隙度預測模型推廣到整個地震剖面(圖9a)，預測孔隙度剖面分別為圖9b、c.自然電位(SP)曲線在該工區能很好地刻畫河道砂體，通過插入井位置的SP曲線(圖9b、c中的虛線)可以看出，整體上有監督網絡和半監督網絡預測的孔隙度具有一定的相似性，但是后者在井旁預測的孔隙度與SP曲線變化趨勢更為吻合，更好地表征出不同河道砂體孔隙度的差異.有監督或半監督網絡預測的孔隙度模型在剖面上部(1.1～1.25 s)以相對低孔地層包裹高孔地層為主，符合曲流河出現的“泥包砂”現象；在剖面下部(1.25～1.32 s)相對高孔地層包裹低孔地層為主，符合辮狀河呈現的“砂包泥”現象.但是，半監督方法比有監督方法預測出的孔隙度結果縱橫向分辨率更高，連續性更好.此外，半監督孔隙度預測結果(圖9c)對應的生成地震數據(圖9d)與實際地震數據(圖9a)接近，二者殘差小(圖9e)，表明預測的孔隙度剖面是滿足地震數據匹配的.

進一步局部放大圖9a—c，比較有監督方法與半監督方法描述辮狀河儲層段孔隙度的差異.圖10a—c分別為儲層段的過井地震剖面，有監督和半監督方法預測的儲層段孔隙度.其中，圖10a—c中2條紅線表示儲層段上下界面位置的兩個目標層位，圖10b—c中的黑框部分為井位置的真實孔隙度，井間為預測孔隙度.相比于有監督方法，半監督方法預測的孔隙度橫向變化在盲井位置(A6)附近更為合理.從圖10b—c中的黑框可以看到，A5井位置的孔隙度值比其余5口井的要大.而有監督方法建立的地震到孔隙度的非線性映射關系很大程度依賴于測井標簽，導致A5井參與有監督網絡訓練后，A4井與A5井之間孔隙度預測結果(圖10b的白框)出現成片的高孔特征，這與辮狀河橫向砂體變化快的地質認識不符.半監督方法比有監督方法提高了地震數據在孔隙度預測中的參與程度，緩解了特殊測井標簽(A5)對網絡建模的影響.通過地震波形匹配約束降低了數據驅動逆過程的多解性，間接保證不同空間位置孔隙度的差異性.因此，在A4井與A5井之間的半監督孔隙度結果(圖10c的白框)更符合地質認識，表明半監督方法在儲層段刻畫不同尺度河道砂體的孔隙度的能力更強.最后，對比一口訓練井(A1)和一口測試盲井(A6)的單道孔隙度預測結果可以看到，半監督孔隙度預測方法(圖11中的紅線)和有監督孔隙度預測方法(圖11中的藍線)都能較好地實現訓練井位置的孔隙度曲線擬合，但是前者的盲井擬合效果更符合實際孔隙度曲線(圖11b中的黑線)的整體趨勢，刻畫出更為準確的小層孔隙度變化細節.在盲井位置，半監督方法和有監督方法預測的孔隙度結果與真實孔隙度的相關系數分別為為0.8151和0.6949.

圖11 基于有監督與半監督網絡的訓練井與測試井預測孔隙度結果對比(a) 訓練井(A1)位置的孔隙度預測結果； (b) 測試井(A6)位置的孔隙度預測結果. (a)、(b)中的黑線、藍線和紅線分別代表真實、有監督網絡預測、半監督網絡預測的孔隙度曲線.Fig.11 Comparisons of predicted porosity curves between supervised network and semi-supervised network at the training well and test well locations(a) Predicted porosity curves at the training well (A1)； (b) Predicted porosity curves at the test well (A6). Black，blue and red lines in (a) and (b) represent true，supervised network inverted，and semi-supervised network inverted porosity curves，respectively.

實際數據測試結果表明：半監督方法比有監督方法更適用于表征實際非均質性儲層的孔隙度.它通過構建地震測井雙約束目標函數降低了測井標簽對孔隙度建模的決定性作用，地震數據波形匹配控制緩解了特殊測井數據對建立的孔隙度預測網絡的負面影響，并間接保證橫向孔隙度預測精度，實現更充分地利用地震測井多尺度信息估計更為準確的孔隙度模型.

3 結論

本文通過模擬傳統的地震正反演與巖石物理建模過程，提出了一種基于雙向門控遞歸單元神經網絡的半監督井震聯合儲層橫向孔隙度預測方法.通過充分利用少量的(偽)測井數據及對應的井旁地震道和大量的非井旁地震道建立地震數據到孔隙度、孔隙度到地震數據的閉環映射，實現地震數據直接預測孔隙度.合成數據和實際數據預測結果表明，相比于僅使用成對的地震數據和測井標簽建立的有監督孔隙度預測方法，半監督孔隙度預測方法具有更好的泛化性和魯棒性.其在保證井位置預測精度的同時通過非井旁地震數據匹配進一步提高井間孔隙度預測精度，進而改善預測孔隙度剖面的橫向連續性和準確性.通過學習不同位置不同信噪比的實際地震數據，半監督方法更好地適應地震波形變化，并降低對測井標簽的過度依賴，緩解特殊測井數據引起的孔隙度建模誤差.因此，半監督孔隙度預測方法更適用于實際非均質性儲層，實現合理地刻畫不同辮狀河河道砂體形態及其內部的孔隙度差異.

地震數據匹配項的正則化參數直接影響半監督孔隙度預測方法的應用效果，如何自動優選出合理的正則化參數有待進一步研究.此外，半監督孔隙度預測方法本身受到地震測井數據質量的固有限制，提高地震數據的信噪比和分辨率，制作數量充足、質量可靠、更貼近實際的(相控)孔隙度曲線也是下一步需要研究的工作.