基于先驗約束的深度學習地震波阻抗反演方法

宋 磊 印興耀* 宗兆云 李炳凱 瞿曉陽 郗曉萍

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580； ②東方地球物理公司研究院地質研究中心，河北涿州 072751)

0 引言

地震波阻抗反演以地震資料為基礎，通過綜合利用地質和測井資料從有限頻寬的地震數據中恢復得到寬頻帶波阻抗數據，現已廣泛應用于油氣勘探階段的儲層定性、定量預測和油氣開發階段的井網部署、儲量計算、油藏動態監測等方面[1-2]。由于實際地球物理問題十分復雜，人們對其認知通常又十分模糊，這使得在反演波阻抗時構建的模型大多是近似的。而深度學習算法可從大量已有數據中學習到隱含知識，從而使所構建數學模型更有利于對那些知識背景不清晰、模型欠精確的問題進行求解，因此該算法適用于波阻抗反演。

深度學習源于神經網絡研究，此概念由Hinton等[3]率先提出。相比于提升方法、支持向量機、最大熵方法等通過人工抽取樣本特征且只能學習到沒有層次結構的單層特征的“淺層學習”方法而言，深度學習通過對輸入樣本逐層進行特征變換以學習得到數據的層次化特征表示[4-6]。深度學習模型的復雜度高，非線性層級數大，可實現對復雜函數的逼近，但通常復雜模型的訓練效率低、易陷入過擬合。而隨著云計算、大數據時代的到來，高性能計算集群有效地緩解了訓練的低效性，海量訓練數據降低了過擬合風險，促進了深度學習的研究與應用的逐步升溫[7]。目前深度學習已廣泛應用于計算機視覺[8]、圖像處理[9]、語音識別[10]和自然語言處理[11]等領域，并給這些領域帶來了巨大的變革。

神經網絡技術在油氣地球物理領域應用較早，印興耀等[12-14]在上世紀90年代初就開展了基于神經網絡的儲層預測技術研究，并取得良好應用效果。近年來隨著深度學習的興起，國內外眾多學者基于此開展了大量研究，如利用深度學習算法識別斷層和鹽體[15-17]，進行地震相分類[18]，提取彈性參數[19-20]，進行儲層預測等[21-22]； 在波阻抗反演方面，依據算法所采用的學習方式分為三類主要方法：基于監督學習的波阻抗反演[23]、基于無監督學習的波阻抗反演[24]和基于半監督學習的波阻抗反演[25-26]。

基于監督學習的波阻抗反演方法需耗費大量時間制作標簽數據(通常來源于測井波阻抗數據)集，且標簽數據的質量也會影響網絡的反演精度。無監督學習波阻抗反演方法將物理模型與神經網絡結合，在網絡學習過程中引入正演模型。這種無監督學習方式減少了網絡對標簽數據的依賴，但同時因缺少真實標簽數據的約束，故其反演精度通常低于監督學習地震反演方法。

半監督學習波阻抗反演方法結合監督學習反演方法與無監督學習反演方法的優勢，同時使用大量無標簽數據(指無對應測井數據的地震數據)和少量有標簽數據(有對應測井數據的地震數據)對網絡模型進行訓練。訓練時將地震數據輸入網絡模型，得到估計的波阻抗； 將該波阻抗再輸入到正演模型，得到合成記錄； 通過對比合成記錄與輸入記錄、標簽數據與網絡估計的波阻抗數據計算損失函數優化網絡模型。因此，半監督學習反演方法既減少了訓練時獲取大量波阻抗的難度，又將井中波阻抗作為標簽約束網絡反演，保證了反演精度。但此類方法仍存在一些問題： ①在構造復雜區，數據特征也十分復雜，在訓練網絡模型時，為了更精確地表征地震數據與波阻抗數據之間的映射關系，可能會丟失較多局部特征，使網絡反演結果喪失一些細節信息。盡管可通過增加網絡深度或延長訓練周期以便于網絡學習這些局部特征，但由于數據量限制，很容易陷入過擬合。②由于地震數據采集設備的限制及處理過程中的低頻損失，獲取的地震數據通常缺少低頻成分，由于深度學習的基礎是數據，若僅據帶限地震數據對網絡模型進行訓練，則網絡所預測的波阻抗也會缺少低頻成分。而低頻成分蘊含地質構造的基本信息，缺失該成分在一定程度上會降低基于波阻抗反演結果開展儲層預測的精度。③受現場環境及采集條件影響，所采集原始數據通常包含大量噪聲，盡管可在處理階段降噪，但這樣也會損害其中的有效信息，如何提高網絡模型對含噪地震數據的反演精度，也是亟待解決的問題之一。

為此，本文提出基于先驗約束的深度學習波阻抗反演方法，即根據待反演區地震相類型對其進行區域劃分，并為每個區域生成相應子網絡，使每個子網絡僅用于對應區域反演任務； 將初始模型用作一類特殊標簽以豐富反演結果的低頻信息、提高反演結果橫向連續性； 在模型中引入具有生物真實性的強抗噪神經元激活函數Rand Softplus[27]以提高其抗噪性能。最后將上述方法針對Marmousi 2模型進行測試，并將其成功地應用于M油田的實際地震數據的反演。

1 方法原理

1.1 反演流程

反演的具體流程主要包括五個方面。

(1)區域分割：將已取得的地震相分類結果作為最初的區域分割結果，為保證每個分割區域內具有較充足數據量并存在一定數量的井數據，合并相鄰的內部反射結構類似且面積較小的區域。

(2)數據預處理：首先對地震數據與測井數據進行標準化處理，再進行井震標定，構建時深關系，然后根據每個區域內的測井數據與地震數據為每個區域提取區域子波并構建區域子波集，最后根據測井數據與構造解釋結果構建初始模型。

(3)訓練反演模塊：反演模塊由一系列網絡模型構成，這些網絡模型是根據區域分割結果自動生成的，每個網絡模型唯一對應一個分割區域。網絡模型的輸入是地震數據，輸出是波阻抗數據，訓練時使用與網絡模型對應分割區域內的訓練數據，并采取半監督學習方式訓練這些網絡模型。

(4)波阻抗反演：將地震數據依據區域分割結果輸入到反演模塊的對應網絡模型中，再合并反演模塊中所有網絡模型預測的波阻抗以獲得最終反演結果，最后對反演結果進行反標準化處理。

(5)評價反演結果：利用均方誤差、相關系數和確定系數三項指標定量地評價反演結果。

1.2 網絡模型

循環神經網絡(RNN)具有一個隱藏狀態向量，可在序列數據之間傳遞，使其能捕獲序列數據中的長距離依賴關系。而卷積神經網絡(CNN)使用固定大小的卷積核并通過滑動窗口方式對輸入數據進行處理，故在局部特征提取方面有較大優勢。本文充分考慮了二者的優勢而構建了深度學習網絡模型框架(圖1)，包含全局特征提取層、局部特征提取層和回歸層三個主要部分。

圖1 深度學習網絡模型框架

1.2.1 全局特征提取層

全局特征提取層的功能是提取輸入地震數據中長時間依賴關系，即全局特征，由一系列GRU (Gate recurrent unit)[28]構成。GRU是RNN的一種，用于處理長期記憶與反向傳播時梯度消失的問題，它使用重置門和更新門控制信息流動。

在圖2所示結構中，前向傳播過程可表示為

圖2 GRU結構示意圖

(1)

式中：rt、ut分別為重置門和更新門的門控信號，rt門控信號用于控制是否忽略之前的隱藏狀態，而ut用于控制隱藏狀態變量的更新；xt、yt、ht、ht-1分別為當前時刻的輸入、輸出、狀態變量及前一時刻狀態變量；W和b分別是權值矩陣和偏置向量。

1.2.2 局部特征提取層

局部特征提取層的功能是提取輸入地震數據的局部特征。由于卷積層中的膨脹系數(dilation)可控制卷積核的感受野大小，所以首先利用一組膨脹系數不同的并行卷積塊提取地震數據中不同尺度的局部特征，再使用一個全連接層和一個卷積塊組合局部特征。這里的卷積塊由一維卷積層(Conv1d)、批歸一化層(Batch normalization)[29]和Rand Softplus激活函數組成。

卷積層由包含多個神經元的特征面組成，每個神經元通過卷積核與前一層特征面局部區域相連以提取輸入特征[30]。卷積層的前向傳播過程表示為

yl=yl-1*Wl+bl

(2)

式中：yl-1、yl分別為第l-1、第l層輸出；Wl、bl是第l層卷積核和偏置向量； “*”是卷積運算符。

批歸一化層用于處理中間層數據分布變化問題，加速網絡的收斂速度。其定義為

(3)

式中：μ、σ分別為輸入數據的均值和標準差；ε是常量，用于保證分式中的分母不為0；λ、ψ均為可學習參數。

激活函數是網絡模型的重要組成部分，是網絡非線性的來源，使網絡能逼近高度非線性函數。常用的激活函數有Tanh、Relu、Softplus等。由于這些激活函數只是對生物神經元輸出相應特性的高度簡化與模擬，使用這些激活函數的網絡模型在抗噪性、不確定性信息處理及功耗等方面與人腦仍存在巨大差距。Rand Softplus激活函數通過引入反映生物神經元隨機性的參數，將Relu與Softplus函數有機結合，使其具有強抗噪能力，定義為

f(x)=(1-η)max(0,x)+ηln(1+ex)

(4)

式中：x為神經元輸入值；η為生物神經元隨機性參數，其值與輸入數據的噪聲水平有關，詳細的計算步驟[27]如下。

(1)計算神經網絡輸入數據所含噪聲標準差

(5)

(2)調整噪聲標準差：由于同一層內部分標準差遠大于該層均值，為避免歸一化操作導致其分布不均勻，將大于2倍均值的標準差設定為2倍均值。

(3)噪聲標準差歸一化處理：隨著網絡深度增加噪聲標準差呈逐漸減小趨勢，為此對噪聲標準差做歸一化處理，從而得到參數η。

(4)η置零操作：為控制網絡的稀疏性，設置比例系數閾值ηth，將η中小于ηth的部分置為0。

1.2.3 回歸層

回歸層的功能是將數據從特征域映射到目標域，即首先使用一個全連接層來組合全局特征提取層和局部特征提取層所提取的特征。由于在特征提取時改變了輸入數據尺寸，所以再使用一組串行反卷積塊對合并后特征進行上采樣，使其具有與標簽數據相同的采樣率； 再使用一個GRU和一個全連接層將上采樣后數據從特征域映射到目標域。其中反卷積塊由反卷積層、批歸一化層和Rand Softplus激活函數組成。反卷積層與卷積層思路類似，僅是運算不同，可將其理解為卷積層的逆過程，其前向傳播過程與卷積層的反向傳播過程相同，反向傳播過程與卷積層的前向傳播過程相同[31]。

1.3 反演框架

本文采用半監督學習方式，并構建圖3所示的基于先驗信息約束的半監督學習反演框架。

圖3 本文方法反演框架

該反演框架由數據分類、反演、正演和數據合并等四個主要模塊組成。數據分類模塊的功能是將輸入地震數據按區域分割結果進行分類。為消除分塊訓練后的拼接痕跡，在對區域邊界處數據進行分類時，對數據向上、向下延拓半個子波長度； 當可延拓長度小于半個子波長度時，延拓至整個二維剖面邊界處； 最后將分類后數據依次輸入到反演模塊的相應網絡模型中。反演模塊由n個網絡模型組成，每個網絡模型(圖1)對應一個分割區域，網絡模型的輸入是地震數據，輸出是波阻抗，訓練時使用相同分割區域內的訓練數據訓練相應網絡模型。數據合并模塊的功能是根據區域分割結果將各分割區域內的波阻抗數據或地震數據進行合并，從而得到與輸入地震數據對應的波阻抗數據或合成地震數據。正演模塊的功能是將波阻抗進行正演，從而得到合成地震記錄。即首先根據

(6)

將輸入的波阻抗轉換為反射系數； 然后從區域子波集中選取與輸入波阻抗數據屬同一分割區域的區域子波； 最后將反射系數與區域子波進行褶積，得到合成地震記錄。式中IP為縱波阻抗序列。

訓練時，反演模塊中各個網絡模型內部可學習參數的優化受以下三個過程的綜合影響：

(1)無標簽地震數據輸入到數據分類模塊后，地震數據按所屬區域類別輸入到反演模塊的相應網絡模型中，得到多組預測波阻抗數據，這些波阻抗數據經過正演模塊后又得到多組合成地震數據，通過計算合成地震數據與分類后輸入地震數據之間的均方誤差(lseismic)，更新反演模塊中相應網絡模型內部參數。

(2)將井旁道地震數據輸入數據分類模塊進行分類，同樣將分類后地震數據按類別輸入到反演模塊的相應網絡模型中，得到多組預測波阻抗數據，通過計算這些預測波阻抗數據與相應井段之間的均方誤差(lwell)，更新反演模塊中相應網絡模型內部參數。

(3)將所有地震數據輸入到數據分類模塊中進行分類，再按類別將地震數據輸入到反演模塊的相應網絡模型中，得到多組預測波阻抗數據，計算這些預測波阻抗與對應初始模型數據之間的均方誤差(lmodel)，從而更新反演模塊中相應網絡模型內部參數。

在訓練初期，初始模型可使網絡快速收斂到某一局部低值，但由于初始模型與地下真實波阻抗存在一定差距，若在訓練后期仍將其作為標簽對網絡進行訓練，則會使反演模塊的預測結果過度模型化，因此應隨著訓練的進行對初始模型進行更新。本次采取每隔一定訓練周期將反演模塊預測的多組波阻抗輸入到數據合并模塊中進行合并，再用合并后波阻抗數據更新初始模型。

綜上，構建如下損失函數

(7)

(8)

反演模塊中網絡模型優化完成后，便可利用上述反演框架開展波阻抗反演，其具體過程是：將待反演地震數據輸入到數據分類模塊進行分類，即將地震數據按類別輸入到反演模塊的相應網絡模型中； 通過反演模塊得到一組預測波阻抗數據； 再將這組預測波阻抗輸入到數據合并模塊，得到反演結果； 最后再對該反演結果進行定量評價。

1.4 反演結果的評價

反演模塊經訓練后，地震數據通過數據分類、反演、合并模塊可得到預測波阻抗數據，地震數據通過數據分類、反演、正演、合并模塊可得到合成地震數據。通過計算真實波阻抗與預測波阻抗、輸入地震數據與合成地震數據的均方誤差(MSE)、相關系數(PCC)和確定系數(R2)，定量評價反演結果。

(1)MSE是預測數據與真實數據對應點誤差的平方和，其值越小說明預測數據與原始數據擬合越好，其定義如式(8)所示。

(9)

(3)R2用于評價變量之間的擬合優度，值域為[0，1]，其值越大，表明變量之間擬合優度越大，自變量對因變量的解釋程度越高，自變量引起的變動占總變動的百分比越高。其定義為

(10)

2 模型數據測試

2.1 無噪聲條件下的波阻抗反演

選取Marmousi 2模型(圖4a)，生成波阻抗數據(圖4b)，再根據下式對波阻抗數據做標準化處理

(11)

式中x是任意輸入數據。將標準化后的波阻抗數據與主頻為30Hz的Ricker子波輸入到正演模塊，得到相應的正演地震數據； 再對該地震數據進行濾波，使其缺失5Hz以下頻帶(圖4a)。將地震相劃分結果作為區域分割結果的初值，對相鄰的內部反射結構類似且面積較小的區域進行合并，如將剖面左側的透鏡體與周圍區域進行合并(圖4c區域Ⅰ)、將剖面中部的兩個斷層結構進行合并(圖4c區域Ⅱ)，得到最終的區域分割結果(圖4c)。均勻地抽取20道的波阻抗作為偽井數據。然后進行層位解釋及斷層解釋，據此將二維剖面分割為多個區塊。針對具有測井數據的區塊，將測井數據在該區塊內進行內插外推； 而對于無測井數據的區塊，則根據相鄰區塊內的數據對該區塊做插值或延伸，生成圖4d所示的初始模型。

按照DL/T 573-95《電力變壓器檢修導則》規定，變壓器一般在投入運行后的5年內和以后每間隔10年大修一次。我廠#6B高廠變為N.M.G變壓器廠生產的、接線組別為△/YO/YO-1-1的變壓器 （額定電壓:22/6.3V）。1996年投運，當時對變壓器內部進行了全面檢查，未發現問題，修后運行狀況良好。根據我廠檢修計劃的安排，決定在2015年#6機組大修中對該變壓器進行全面預試。

圖4 Marmousi2模型數據

使用地震數據、偽井、初始模型和子波訓練反演模塊中的網絡模型，訓練時將迭代周期設置為1000，每個迭代周期的批訓練量為50，損失函數中的α、β、γ分別設置為1.0、1.0、0.1，訓練時每隔50個迭代周期更新一次初始模型。隨后用訓練好的反演模塊進行波阻抗反演(圖5a，設定為Result1)。

基于上述反演框架，首先測試無初始模型約束條件下反演效果，即將γ設置為0，得到了圖5c所示反演結果(設定為Result2)； 然后測試無初始模型和區域分割結果約束條件下的反演效果，即將γ設置為0且不進行區域分割(認為整個地震剖面屬于同一區域)，得到圖5e所示反演結果(設定為Result3)。表1顯示不同先驗約束條件下網絡在訓練階段與預測階段的計算時長(Intel i7 9700、Nvidia GeForce GTX 1050Ti)，發現引入各種先驗信息約束后，網絡的計算時間仍保持在可接受范圍內。

表1 不同先驗約束下網絡的訓練和預測階段運行時長

為了更直觀地對比不同先驗約束條件下的反演效果，根據

(12)

分別計算Result1、Result2與Result3與真實波阻抗的相對差值，結果如圖5b、圖5d和圖5f所示。

由于事先已知模型的波阻抗，可通過計算Result1、Result2和Result3與真實波阻抗的MSE、PCC和R2平均值定量評價不同約束條件下網絡模型的反演效果。表2、表3和表4分別為Result1、Result2和Result3與真實波阻抗在整個剖面上以及在圖5中紅色矩形1、2、3標示區域的MSE、PCC和R2平均值。

分析、對比圖5中不同先驗約束條件下的反演結果及表2～表4中定量評價結果，可直觀地發現，充分利用初始模型和區域分割結果等先驗信息能較大程度上消除反演結果中縱向抖動現象，提高反演結果的橫向連續性； 且在一些構造復雜部位，這種效果改善更明顯。如圖5a、圖5c和圖5e中紅色矩形1、2、3標示區域，這主要是由于該區域地下構造較復雜，存在構造突變和斷層，此時初始模型和區域分割結果向網絡提供基本構造信息，并對網絡預測結果做出明確的空間約束。而且反演模塊中每個網絡模型只負責構造形態相似區域的反演，使這些網絡模型能從訓練數據中捕捉更多局部特征，在使用訓練好的網絡模型進行反演時，其反演結果的細節信息也會更豐富，即反演精度更高。

表2 Result1的定量評價

表3 Result2的定量評價

表4 Result3的定量評價

但在圖5a中紅色箭頭所指區域的反演效果卻出現了退化，這主要是由于在進行區域分割時該處(圖4cⅠ區域)縱向長度較短，由于地震數據的自相關性，當將該處地震數據輸入到網絡模型時，網絡易將該處地震數據與其他某些道集數據相混淆，所以呈現反演精度降低的現象。

圖5 不同先驗信息約束條件下的反演結果(左)及其與真實波阻抗的相對差值(右)

2.2 含噪聲條件下的波阻抗反演

為測試網絡模型對噪聲數據的適應能力，向原地震數據中分別加入4%、8%、12%的高斯隨機噪聲，然后在網絡模型中分別使用Rand Softplus和Tanh激活函數對含噪地震數據進行反演。圖6顯示0、4%、8%和12%含噪條件下使用Rand Softplus和Tanh兩種激活函數的反演結果。圖7為上述反演結果與真實波阻抗的相對差值。表5和表6歸納在上述噪聲條件下所得反演結果的定量評價。

表5 基于Rand Softplus激活函數的反演結果定量評價

表6 基于Tanh激活函數的反演結果定量評價

圖6 不同噪聲條件下使用Rand Softplus(左)和Tanh(右)激活函數的網絡模型的反演結果

圖7 圖6所得反演結果與真實波阻抗的相對差值

3 實際數據應用

將本文反演方法應用于M油田實際三維地震數據。所選測線(圖8)穿過3口井(W1、W2、W3)，包含427個CDP，縱向時間長度為378ms、采樣間隔為2ms。區內Hor3界面(紅線)的上部與下部地震數據的反射結構有較大差異，下部振幅強度大、連續性好、內部反射結構有序； 而上部振幅強度小、連續性差、內部反射結構相對雜亂，因此以Hor3為分界線將該剖面分為上、下兩個區域。本次將W1和W3設定為訓練井，W2作為測試井，訓練井用作網絡標簽，而測試井不參與子波提取、初始模型構建與反演模塊的訓練。

圖8 地震數據、測井數據與地震層位解釋結果

首先對數據進行預處理。

(1)據式(11)對井數據、地震數據做標準化處理。

(2)井震標定，構建時深關系； 再從井旁地震道分別提取Hor3層上、下相鄰區域的子波(圖9)。

圖9 Hor3層上、下相鄰區域的子波

(3)當標簽數據與輸入地震數據的頻譜范圍差異較大時，不利于網絡快速收斂，還會降低網絡預測精度[33]，因此需濾去測井數據中的高頻成分。

(4)根據地震層位解釋結果，并綜合利用測井數據，構建初始模型(圖10)。

圖10 本文方法構建的初始模型

然后用區域子波、井數據、地震數據與初始模型訓練反演模塊中的網絡模型，訓練周期為500，每個周期的批訓練量為20，損失函數中的α、β、γ分別設定為0.2、1.0、0.1，訓練時每隔30個周期更新一次初始模型。

從上述反演結果可知，本文反演方法不僅在訓練數據集中有標記數據上具有較好反演效果，且在無標記數據上仍能取得較高反演精度。合成地震記錄與輸入地震記錄的對比結果(圖11d)充分證明了本文方法反演結果的合理性與準確性； 與稀疏脈沖反演法結果(圖11b)相比，本文方法反演結果中橫向連續性更好，且層與層之間的粘連效應較弱。

圖11 本文方法在實際數據中的應用

4 結論

在深度學習地震反演算法中，本文嘗試引入不同的先驗約束條件，由此提出一種基于先驗約束的深度學習地震反演方法。為增強算法的穩定性與抗噪性，并減少算法對標簽數據的依賴性，在網絡模型中選用了強抗噪性的Rand Softplus激活函數，且采取半監督學習方式對網絡模型進行訓練。經過理論分析與實驗測試取得以下認識和結論：

(1)采用區域分割策略后，每個網絡模型只需負責數據特征相似區域的反演，這使網絡模型能從訓練數據中提取更多有效信息，從而豐富了反演結果的細節信息，提高了反演精度。但由于地震數據的自相關性，在一些縱向長度較短的地震道，反演精度會有所退化。

(2)將初始模型作為一種特殊標簽來約束網絡的反演過程，可使網絡充分利用初始模型中含有的低頻分量和構造細節信息，從而豐富了反演結果的低頻信息，提高了反演結果的橫向連續性。

(3)Rand Softplus激活函數通過引入反映數據隨機性的參數，提高了網絡對噪聲數據的適應能力。尤其是在強噪聲條件下，該激活函數對反演精度的提升更為明顯。