田亞軍, 高靜懷* , 王大興, 陳道雨

1 西安交通大學信息與通信工程學院， 西安 710049 2 海洋石油勘探國家工程實驗室， 西安 710049 3 中國石油長慶油田公司勘探開發(fā)研究院， 西安 710018 4 低滲透油氣田勘探開發(fā)國家工程實驗室， 西安 710018

0 引言

在儲層預測工作中常遇到強反射界面，這給儲層的精確描述工作帶來了極大困難.相鄰地層波阻抗差異大是導致強反射界面出現(xiàn)的直接原因.一般煤系地層、深水穩(wěn)定沉積形成的烴源巖層、風化層形成的不整合面等位置易形成強反射界面，而這類地層附近也容易產(chǎn)生油氣聚集.強反射的存在，會淹沒其附近儲層的弱反射信號，導致儲層預測工作無法正常進行.為此，開展強反射剝離研究對于儲層弱反射信號的識別、結(jié)構(gòu)的刻畫、含油性檢測并進一步輔助水平井設計提高鉆探成功率具有重要意義(Guo and Wang, 2004；李海山等，2014；秦雪霏和李巍，2014；杜玉山等，2016；蔡文杰等，2019).

針對強反射附近的儲層預測問題前人提出了諸多解決方案，主要包括三類思路：(1)去除或壓制強反射，其中基于匹配追蹤算法的研究最多，主要思想是根據(jù)強反射低頻強能量的特性匹配出強反射分量然后將其去除(Wang, 2010；李海山等，2014；許璐等，2019；何峰等，2019).多子波分解方法通過不同主頻的地震子波重構(gòu)原始信號，從重構(gòu)分量中去除反映強反射信號頻段的分量達到去除強反射的目的(趙爽等，2007；佘剛等，2013；陳人杰等，2014).類似的，Wang等(2016a)使用EMD分解方法分解出強反射所在的固有模態(tài)分量并將其分離達到壓制強反射的目的.(2)忽略強反射層的影響，通過調(diào)整時頻分析方法的參數(shù)突出弱反射信號的時頻特征(王寶江等，2012).(3)采用反Q濾波方法對強反射附近目標層反射波的相位和振幅進行補償校正來增強弱反射信號的能量(Guo and Wang, 2004).本研究屬于第(1)類.

在卷積模型的框架下可將強反射看成是強反射系數(shù)與地震子波卷積的結(jié)果(Robinson, 1985).構(gòu)成強反射的關(guān)鍵要素包括強反射系數(shù)的位置、大小以及地震子波，若能夠通過算法求解這三個要素，強反射就能被唯一的確定.然而，實際資料中的地震子波難以預測，加上強反射往往是由多個反射界面產(chǎn)生的反射波相互干涉疊加形成的復合波，導致原始地震記錄與強反射之間存在非常復雜的非線性映射關(guān)系，傳統(tǒng)方法難以求解這一映射關(guān)系.深度學習作為機器學習領(lǐng)域一個新的研究方向，不需要物理方程或理論知識就可以求解非常復雜的非線性映射問題(Kim and Nakata, 2018).更重要的是，當深度網(wǎng)絡訓練完成后可直接預測目標，不需要經(jīng)驗性調(diào)參.鑒于深度學習方法的強擬合性，大量的深度學習方法被應用于地震資料處理與解釋，例如波阻抗反演(Das et al., 2019; Wu et al., 2020)、鹽丘識別(Di and AlRegib, 2020; Shi et al., 2019)、初至拾取(Wu et al., 2019a；李薇薇等，2021)、層位解釋(Lou et al., 2020; Wu et al., 2019b)、斷層解釋(段艷廷等，2019；Wu et al., 2019c; Liu et al., 2020)、儲層預測(林年添等，2018；付超等，2018)、波形分類(趙明等，2019)、地震P波S波到時拾取(于子葉等，2018)以及地震成像(奚先和黃江清，2018，2020).

綜上所述，本研究通過分析強反射產(chǎn)生的原因，在卷積模型的框架下建立含有強反射地震記錄的生成模型；進一步將強反射預測問題分解成地震子波估計與強反射估計兩個子優(yōu)化問題，并引入交替迭代深度神經(jīng)網(wǎng)絡(AIDNN)和U-Net網(wǎng)絡分別求解；根據(jù)求解優(yōu)化模型的需求及訓練樣本的特點設計具體工作流程.最后通過模型數(shù)據(jù)及實際數(shù)據(jù)的試算，驗證方法的有效性.

1 強反射層剝離問題分析及模型構(gòu)建

地震記錄可以看成地震子波與反射系數(shù)卷積的結(jié)果(Robinson, 1985)：

s(t)=w(t)*r(t),

(1)

其中s(t)是合成地震記錄，*是卷積運算符，w(t)是地震子波，r(t)是反射系數(shù).同樣的，對于含有強反射的地震記錄可以表示為：

s(t)=w(t)*rw(t)+w(t)*rs(t),

(2)

其中rw(t)是弱反射系數(shù)，rs(t)是強反射系數(shù).圖1簡要闡明了帶有強反射的地震記錄的形成過程.顯然，若能預測出強反射(ss(t)=w(t)*rs(t))并將其從地震記錄s(t)中剝離，就能突出弱有效信號(sw(t)=w(t)*rw(t))能量.參照卷積模型框架下的反問題解法(Wang et al., 2016a)，本研究可以構(gòu)建如下優(yōu)化問題求解強反射：

(3)

式(3)中共有三個未知量：w(t)、rw(t)和rs(t)；這種含有三個未知量的反問題無法直接求解.

圖1 強反射剝離原理Fig.1 Principle of strong reflection removal

強反射可以看成是強反射系數(shù)與地震子波卷積的結(jié)果；實際工作中，在一定空間范圍內(nèi)地震子波可以看成是時不變的；本文首先假設地震子波在待處理區(qū)域時不變，這樣就能將式(3)中的地震子波分離出來單獨求解.參照Chen 等(2019)的工作將地震子波預測問題分解為地震子波與反射系數(shù)交替求解的兩個子問題，并采用AIDNN網(wǎng)絡訓練并求解.得到地震子波后，式(3)便剩下兩個未知量rw(t)和rs(t).建立如下優(yōu)化模型：

(4)

進一步假設在一定空間范圍內(nèi)反射系數(shù)序列同分布，選取同工區(qū)不包含強反射的地震數(shù)據(jù)作為弱反射地震記錄sw(t)來訓練式(4)中的優(yōu)化問題.式(4)可以表示為：

(5)

式(5)中僅包含一個未知量rs(t).

由于我們想得到的是強反射，為了方便訓練地震記錄與強反射的映射關(guān)系，將式(5)中的唯一未知量代換為ss(t)=w(t)*rs(t)；式(5)可表示為：

(6)

以相同工區(qū)相鄰位置的數(shù)據(jù)作為訓練樣本帶來的另一個好處是，樣本sw(t)本身的噪聲水平與分布規(guī)律可以看作與待處理地震記錄s(t)一致.因此，在訓練強反射ss(t)與地震記錄s(t)之間的映射關(guān)系時不必考慮噪聲.最后采用具有強大非線性映射能力的U-Net來求解這一問題.

2 基于AIDNN的地震子波預測

Adler和?ktem(2017)針對病態(tài)反問題提出了一種交替迭代深度神經(jīng)網(wǎng)絡(AIDNN)并將其應用于醫(yī)學高分辨成像.AIDNN以傳統(tǒng)最小二乘方法為先驗約束解空間，進一步利用深度網(wǎng)絡校正最小二乘解提高解的精度.Chen等(2019)將AIDNN應用于地震子波預測；本文借鑒Chen等人的工作建立地震子波預測方法.下面簡要介紹該思想：

根據(jù)卷積模型(Robinson, 1985)，地震記錄s(t)可以表示為：

s(t)=ω(t)*r(t)+n(t),

(7)

其中，*是卷積運算符，ω(t)是地震子波，r(t)是反射系數(shù)，n(t)是隨機噪聲.本研究的目標是求解w(t)，建立如下優(yōu)化模型：

+Ψ(w(t)),

(8)

其中Φ(r(t))是r(t)的凸正則化函數(shù)，Ψ(w(t))是w(t)的凸正則化函數(shù).受到Wang等(2016b)、Adler和?ktem(2017)的啟發(fā)，將問題(8)分解為兩個子問題，即：交替反演地震子波與反射系數(shù).具體可分為如下兩個步驟：

(9)

(10)

2.1 AIDNN網(wǎng)絡結(jié)構(gòu)

圖2 訓練地震數(shù)據(jù)與地震子波之間映射關(guān)系的AIDNN網(wǎng)絡結(jié)構(gòu)Fig.2 The structure of AIDNN network for training mapping relationship between seismic wavelets and seismic data

如圖2所示，求解問題(9)與問題(10)所使用的網(wǎng)絡結(jié)構(gòu)相同.除第一次迭代與最后一次迭代包含輸入輸出層外，每次迭代使用的網(wǎng)絡結(jié)構(gòu)相同.每次迭代過程中隱層數(shù)量相同，本文使用三個隱層且每個隱層使用線性整流函數(shù)(ReLU)作為激活函數(shù).卷積核大小為3×1.更多關(guān)于AIDNN網(wǎng)絡理論基礎可參見文獻Adler和?ktem(2017).

2.2 訓練樣本

對于地震子波的預測問題，本研究選用的訓練數(shù)據(jù)集包括地震子波、反射系數(shù)以及隨機選擇地震子波與反射系數(shù)并卷積生成的地震記錄.其中反射系數(shù)序列來自于圖3a所示的Marmousi模型，共包含737道數(shù)據(jù)，振幅范圍歸一化為[-1,1]，采樣間隔為2 ms、時間采樣點數(shù)為551.訓練使用的地震子波為3801個不同相位與主頻的雷克子波(相位范圍為-90°～90°，采樣間隔為1°；主頻范圍為20～40 Hz，

圖3 用于地震子波估計的訓練樣本(a) Marmousi模型反射系數(shù); (b) 雷克子波(主頻20 Hz，相位-90°～90°).Fig.3 Training sample for estimating seismic wavelets(a) Reflection coefficient of Marmousi model； (b) Ricker wavelets (phases are from -90° to 90°, and dominant frequency is 20 Hz).

采樣間隔為1 Hz；圖3b所示為部分主頻為20 Hz，相位介于-90°～90°之間的雷克子波樣本).訓練過程中隨機抽取地震子波與反射系數(shù)并卷積生成地震道.

2.3 訓練策略

3 基于U-Net的強反射層剝離

如式(6)所示，本研究將強反射預測問題構(gòu)建成一個端到端的優(yōu)化問題，并使用U-Net網(wǎng)絡求解.U-Net網(wǎng)絡首先由Ronneberger等(2015)提出并將其應用于醫(yī)學圖像分割.由于強大的非線性映射能力，U-Net被廣泛應用于許多領(lǐng)域.對于地震數(shù)據(jù)處理和解釋，其被應用到初至拾取(Hu et al., 2019)、斷層識別(Wu et al., 2019c; Liu et al., 2020)等問題上，并取得了良好效果.

3.1 U-Net的網(wǎng)絡結(jié)構(gòu)

圖4為U-Net的網(wǎng)絡結(jié)構(gòu)示意圖，該網(wǎng)絡的輸入是帶有強反射的地震記錄道，輸出為強反射.網(wǎng)絡除輸入與輸出層外是一個U形的對稱路徑(即下采樣路徑和上采樣路徑).每個下采樣步驟包含兩個重復的卷積核大小為3×1的卷積層，每個卷積層后面是一個線性整流激活函數(shù)(ReLU)和一個2×1的最大池化層，每次下采樣都將特征圖大小減半.每個上采樣步驟包括一個2×1的上采樣卷積層和兩個重復的卷積核大小為3×1卷積層，每個卷積層后面都有一個ReLU，每次上采樣都將特征圖大小翻一倍.更多關(guān)于U-Net的理論基礎可參見Ronneberger等(2015).

圖4 用于地震強反射預測的U-Net網(wǎng)絡架構(gòu)Fig.4 Framework of U-Net network used in strong seismic reflection prediction

3.2 訓練樣本

在地震子波在一定空間范圍內(nèi)是時不變的且反射系數(shù)序列在一定空間范圍內(nèi)是同分布的假設下，本文建立了式(6)所示的強反射預測優(yōu)化模型.預測強反射前，本文先確定強反射所在時窗范圍，并選擇強反射時窗附近的具有相同窗長的弱反射樣本.在強反射系數(shù)的大小和位置均是隨機的假設下，生成一系列振幅與位置均隨機變化的強反射系數(shù)并與前文基于AIDNN網(wǎng)絡求解的地震子波褶積獲得強反射樣本.隨機抽取弱反射樣本與強反射樣本并將二者相加便得到了含有強反射的地震記錄樣本.至此便獲得求解式(6)所需的所有訓練樣本.具體樣本制作細節(jié)與方法實施步驟見第4節(jié).

4 實施步驟

根據(jù)前文設計的強反射剝離優(yōu)化模型，結(jié)合AIDNN與U-Net兩個深度網(wǎng)絡，本研究制定了圖5所示的技術(shù)實施流程，具體內(nèi)容如下：

(1)強反射層層位Hs拾取.

(2)確定強反射層層位Hs所在的時窗范圍并將時窗范圍分別向上或向下延伸一個波長獲得待處理時窗范圍Wts-te，窗口中的數(shù)據(jù)記為si，其中i代表道號.

(3)地震子波預測訓練樣本制作，包括三個內(nèi)容：①選擇典型地質(zhì)模型的反射系數(shù)作為訓練樣本，本文參照Chen等(2019)的工作選擇Marmousi模型作為反射系數(shù)訓練樣本；②確定地震子波主頻大致范圍，并設置主頻與相位隨機變化的雷克子波為子波的訓練樣本；③隨機抽取反射系數(shù)與地震子波進行卷積獲得地震記錄樣本；本文使用的具體訓練樣本參見2.2節(jié).

(6)待處理數(shù)據(jù)si預處理：強反射與弱反射最大的差異在于能量強，為了更好地使U-Net網(wǎng)絡學習到這一知識，本研究對si的每一道分別做以強反射層Hs為中心的限幅處理使強反射層位置能量最強并對每道分別做歸一化得到sig：

圖5 強反射剝離工作流程圖Fig.5 Work flow of strong reflection removal

sig(t)=si(t)⊙Giw(t)./max(si(t)⊙Giw(t)),

(11)

其中⊙為哈達瑪積，代表對應元素相乘；max代表每道地震記錄的最大值；Giw(t)代表限幅處理使用的高斯窗函數(shù)：

Giw(t)=e-a(t-ti0)2,

(12)

其中a控制高斯窗寬，ti0是每道地震記錄強反射層所在時刻，i代表地震道道號.

(7)強反射訓練樣本制作：在卷積模型框架下，假設地震子波在一定空間范圍內(nèi)時不變，如果已知地震子波、強反射層反射系數(shù)的位置和大小，則強反射就唯一的確定下來了；據(jù)此，本文分三步構(gòu)建強反射樣本：①統(tǒng)計預處理后的數(shù)據(jù)sig(t)每道中強反射層峰值振幅與其次峰值振幅之比，確定所有道的比值范圍Rs-Re；②在時窗Wts-te范圍內(nèi)生成位置隨機且大小在Rs-Re范圍內(nèi)隨機變化的強反射系數(shù)；③與步驟(5)中地震子波w(t)卷積生成強反射樣本sjs，j代表第j次生成的強反射樣本.

(8)弱反射訓練樣本制作：對于實際數(shù)據(jù)，在一定地質(zhì)時期形成的地下地層之間存在必然聯(lián)系并且地下地層的巖性和物性常呈周期性出現(xiàn)；據(jù)此，為了使樣本訓練更具針對性，本研究假設在一定空間范圍內(nèi)的弱反射界面的反射系數(shù)概率同分布；從待處理數(shù)據(jù)所在的時窗Wts-te下部相鄰位置截取等寬的不含強反射的層段Wsts-te數(shù)據(jù)并做單道歸一化作為弱反射樣本siw.

(9)分別隨機抽取sjs和siw中的單道數(shù)據(jù)sks和skw并將二者相加獲得含強反射的地震記錄skr，k代表第k次抽取的訓練樣本.

(10)訓練樣本預處理：分別對sks、skw和skr使用與步驟(6)相同的高斯窗函數(shù)進行限幅處理并分別對三者做除以sks最大值的歸一化處理得到：

sksg(t)=sks(t)⊙Gkw(t)/max(sks(t)⊙Gkw(t)),

(13)

skwg(t)=skw(t)⊙Gkw(t)/max(sks(t)⊙Gkw(t)),

(14)

skrg(t)=skr(t)⊙Gkw(t)/max(sks(t)⊙Gkw(t)),

(15)

sksg、skwg和skrg分別是所需的強反射、弱反射和包含強反射的地震記錄訓練樣本.

(13)使用式(16)逐道恢復強反射層真振幅：

sispr(t)=

(16)

其中/號代表對應元素相除.

(14)逐道從原始記錄si(t)中減去強反射sispr(t)，獲得強反射剝離后地震數(shù)據(jù).

5 模型數(shù)據(jù)應用效果分析

為了驗證本研究提出的強反射剝離方案的可行性和有效性，首先將其應用于理論模型數(shù)據(jù)的驗算.

5.1 二維速度模型的建立與正演模擬

Chen等(2019)在論文中詳細的對比了AIDNN網(wǎng)絡應用于卷積模型制作的合成地震記錄上預測地震子波的效果.為了使模擬結(jié)果更加接近真實情況，本研究使用能夠更好的保證地震波傳播過程中幾何學與動力學特性的波動方程進行正演模擬.圖6a為本研究使用的二維速度模型，紅色箭頭指示將產(chǎn)生強反射的反射界面，藍色箭頭指示需要識別的弱反射界面.子波采用圖7中藍線所示的主頻30 Hz，相位-45°的雷克子波，采樣間隔為2 ms；使用波動方程模擬炮集記錄(Zhu and Harris, 2014; Zhang and Shi, 2019)并抽取零偏移距地震道作為待處理數(shù)據(jù)(圖6b).不難看出，由于強反射層反射信號的存在，圖6b所示的原始地震剖面下部目標層的反射信號淹沒其中，無法識別(即圖6b中無法識別圖6a中藍色箭頭指示的反射界面反射的信號).

圖6 二維地質(zhì)模型及其模擬地震剖面(a) 二維地質(zhì)模型，紅色與藍色箭頭分別指示強反射界面和弱反射界面； (b) 模擬地震剖面，紅框內(nèi)為待處理部分，籃框內(nèi)為弱反射訓練樣本，綠色曲線代表強反射層T0.Fig.6 Two-dimensional geological model and its simulated seismic profile(a) Two-dimensional geological model. The red and blue arrows indicate strong and weak reflection interface, respectively； (b) Simulated seismic profile. The red rectangle denotes the data to be processed. The blue rectangle denotes the training samples of weak reflection. The green curve indicates strong reflector T0.

5.2 應用效果分析

為了剝離強反射凸顯目的層弱反射，本研究按照第4節(jié)中的實施步驟進行逐步處理.需要注意的是，圖6b中紅框內(nèi)部為待處理數(shù)據(jù)，藍框中為弱反射信號樣本所在時窗.圖7a所示為使用AIDNN網(wǎng)絡預測的地震子波與實際子波對比，可見二者波形匹配較好；由圖7b真實子波與預測子波之差可見，誤差最大的位置是子波第二個旁瓣的位置為0.3，進一步說明了地震子波的預測效果；圖8為模型數(shù)據(jù)強反射剝離前后結(jié)果對比分析，可以看到，相比于圖8a所示的原始地震剖面，圖8b所示的強反射剝離后的剖面中藍色箭頭指示的目標層弱反射信號凸顯出來.圖8c為預測出的強反射，其橫向能量變化與波形變化與原始剖面(圖8a)中的強反射一致，這說明了本方法良好的保真保幅性.

圖7 實際地震子波與估計地震子波對比(a) 實際地震子波(藍色曲線)與估計地震子波(紅色曲線)； (b) 預測結(jié)果與真實結(jié)果之差.Fig.7 Comparison of real and estimated seismic wavelets(a) Real (blue curve) and estimated (red curve) seismic wavelets； (b) Difference between estimated and real seismic wavelets.

圖8 模型數(shù)據(jù)的強反射剝離(a) 原始地震剖面； (b) 強反射層剝離后的地震剖面，藍色箭頭指示弱反射信號； (c) 強反射.Fig.8 Strong reflection removal of synthetic data(a) Original seismic profile； (b) Seismic profile after strong reflections removed. Blue arrow indicates weak reflection signal； (c) Strong reflection.

6 實際數(shù)據(jù)應用效果分析

本節(jié)選擇儲層在強反射下方的典型疊后地震剖面作為研究對象，通過本研究提出的強反射剝離方法輔助儲層砂體結(jié)構(gòu)表征.

6.1 實際數(shù)據(jù)及其地質(zhì)背景

待處理數(shù)據(jù)選自中國西北鄂爾多斯盆地古峰莊地區(qū)，目的層是河流沉積體系演化為大型三角洲沉積體系過程中形成的優(yōu)質(zhì)儲層；其上部層段為厚達40 m的低速烴源巖，形成于深水湖盆時期；這為下部儲層提供了良好的生油與封蓋條件，二者形成了最佳的成藏組合關(guān)系.但由于烴源巖與圍巖波阻抗差大導致了強反射同向軸T0的形成，加之烴源巖與儲層相鄰導致下部儲層反射信號被淹沒在強反射之中，難以有效識別，這嚴重影響了儲層預測工作(屈雪峰等，2017).因此，本節(jié)針對該工區(qū)實際疊后地震資料展開強反射剝離工作，以輔助儲層砂體展布形態(tài)的刻畫.

6.2 應用效果分析

圖9為待處理的二維過井疊后地震剖面，共901道，道間距20 m，采樣間隔2 ms.其中紅框內(nèi)為待處理數(shù)據(jù)，綠線代表待剝離強反射層T0，T0下部為河道砂體儲層，探測這部分信號是該工區(qū)重要任務.圖9中藍框內(nèi)的數(shù)據(jù)為弱反射層訓練樣本所在時間窗；另外，為了增加樣本數(shù)量，我們截取臨近30條CMP線上的相同時窗內(nèi)的共27030個地震道作為弱反射樣本.按照第4節(jié)中實施步驟進行逐步處理.圖10所示為地震子波預測結(jié)果.

圖9 待處理的原始地震剖面紅框內(nèi)為待處理數(shù)據(jù)；藍框內(nèi)為弱反射訓練樣本；綠色曲線指示強反射層T0.Fig.9 Original seismic profile to be processedThe red rectangle denotes the data to be processed. The blue rectangle denotes the training samples of weak reflection. The green curve indicates strong reflector T0.

圖10 基于AIDNN估計的實際地震子波Fig.10 Seismic wavelet of real seismic data estimated by AIDNN

圖11為實際資料強反射剝離效果對比分析.從圖11a所示的原始剖面中的自然伽馬曲線(紅色曲線)與聲波速度曲線(藍色曲線)不難看出，強反射下還有一個層位(藍色箭頭指示)；但在地震剖面上該層位反射波淹沒在其上部強反射之中無法識別.圖11b是使用本文方法剝離強反射后的地震剖面，可見強反射下部藍色箭頭指示的弱反射軸顯現(xiàn)出來且其位置與井曲線上指示的位置一致.圖11c為剝離出來的強反射，與圖11a中強反射層對比可見，二者橫向上的能量與波形的變化一致且圖11c中強反射之外無任何能量，這進一步說明了本文方法的有效性.

圖11 實際地震數(shù)據(jù)強反射剝離(a) 原始地震剖面，藍色箭頭指示井曲線上的弱反射層； (b) 強反射層剝離后的地震剖面，藍色箭頭指示弱反射信號； (c) 強反射；另外，(a)—(c)中綠色地震道為合成記錄，紅色曲線和藍色曲線分別是自然伽馬曲線與聲波速度曲線；水平方向的黑色曲線和紅色曲線之間是待預測儲層砂體所在區(qū)間.Fig.11 Strong reflection removal of real data(a) Original seismic profile. The blue arrow indicates the weak reflector on the well curve； (b) Seismic profile with strong reflections removed. Blue arrow indicates weak reflection signal； (c) Strong reflection. In (a) — (c), green seismic traces are synthetic seismic records, red and blue curves are gamma ray curves and acoustic velocity curves, respectively. Between black and red curves in horizontal direction is the sand body to be predicted.

為了進一步說明效果，將實際剖面中井2與井4的測井曲線與井旁道抽取出來進行對比分析(如圖12所示).圖12中黃線位置是強反射層位置，從測井曲線上可以看到其下部藍色箭頭指示的目標層位，強反射軸剝離之前在井旁道上無法顯示這一層位產(chǎn)生的反射信號，將強反射剝離后該層位弱反射能量凸顯出來且位置與井點上顯示的位置一致(藍色箭頭所指).

圖12 在井旁道分析強反射剝離結(jié)果(a) 井2； (b) 井4.Fig.12 Analysis of strong reflection removal results in borehole-side trace(a) Well 2； (b) Well 4.

進一步分別對圖11中強反射剝離前后的地震剖面用優(yōu)化三參數(shù)小波變換(Tian et al., 2020; 高靜懷等，2006；Liu et al., 2017; 高靜懷等，2018)分解為圖13與圖14所示的20 Hz、30 Hz、40 Hz和50 Hz的譜分量來分析儲層砂體展布情況(Lu and Li, 2013; Liu et al., 2019).圖13為強反射層剝離前譜分解結(jié)果，不難看出，低頻部分(圖13a、b)目的層無法和強反射分離，高頻部分(圖13c、d)受強反射影響目的層內(nèi)部能量弱，無法準確分析砂體分布情況.圖14為強反射剝離后的譜分解結(jié)果，可見在目的層范圍內(nèi)不同頻率成分譜分解結(jié)果均不受強反射的影響.從圖14d高頻50 Hz譜分解結(jié)果來看，由淺至深存在三個砂體結(jié)構(gòu)且譜分解結(jié)果能夠清晰的描述三個砂體展布情況(藍箭頭、紅箭頭和黑箭頭分別指示第一至第三個砂體分布情況).此外，觀察五口井的自然伽馬曲線，不難看出，自然伽馬曲線低值位置箭頭指示的砂體位置與50 Hz譜分解結(jié)果一致；相對而言圖13d所示的含強反射的原始剖面的譜分解結(jié)果無法清晰刻畫著三個砂體，尤其是與強反射較近的第一個砂體的展布情況；這進一步說明本研究構(gòu)建的強反射剝離方法的效果及其對這類儲層預測的重要性.

圖13 強反射剝離前譜分解結(jié)果(a)—(d)分別是由三參數(shù)小波變換計算的20 Hz、30 Hz、40 Hz、50 Hz頻譜分量.井點位置黑色曲線代表自然伽馬曲線.Fig.13 Spectral decomposition results before strong reflection removal(a)—(d)are 20 Hz, 30 Hz, 40 Hz and 50 Hz spectral components calculated by the three-parameter wavelet transform, respectively. The black curves on the well points are the gamma measurements.

圖14 強反射剝離后剖面譜分解結(jié)果(a)—(d)分別是由三參數(shù)小波變換計算的20 Hz、30 Hz、40 Hz、50 Hz頻譜分量.井點位置黑色曲線代表自然伽馬曲線.Fig.14 Spectral decomposition results after strong reflection removal(a)—(d)are 20 Hz, 30 Hz, 40 Hz, and 50 Hz spectral components calculated by three-parameter wavelet transform, respectively. The black curve on the well point are the gamma measurements.

進一步對上述工區(qū)55 km2的三維疊后地震數(shù)據(jù)體目的層強反射進行強反射剝離處理驗證本文方法對三維實際數(shù)據(jù)的適用性和實用性.待處理數(shù)據(jù)Inline范圍為150～400間隔為1，Crossline范圍為60～600間隔為1，共135791道地震數(shù)據(jù)，每道351個采樣點，采樣間隔2 ms，數(shù)據(jù)大小212 M.本研究在TensorFlow框架下編寫Python語言程序，在一塊i7-8700 CPU上進行計算，訓練地震記錄與地震子波之間映射關(guān)系需要3 h 31 min，訓練地震記錄與強反射之間的映射關(guān)系需要44 min；對目標區(qū)域的預測及剝離過程使用了30 min.這里地震記錄與地震子波之間映射關(guān)系的訓練結(jié)果可以推廣應用至其他工區(qū)；地震記錄與強反射之間的映射關(guān)系訓練結(jié)果僅能應用于目標工區(qū)，這與選擇的弱反射訓練樣本有關(guān).圖15分別為強反射剝離前后對20 Hz與50 Hz譜分解結(jié)果做的沿層切片，其中20 Hz單頻切片表征厚砂體平面展布范圍，50 Hz單頻切片表征薄砂體平面展布范圍.不難看出，強反射剝離前20 Hz(圖15a)與50 Hz(圖15b)的單頻切片結(jié)果差異小，刻畫的信息大部分為強反射的頻率成分，無法精細刻畫目的層不同厚度砂體平面展布范圍.利用本文提出的強反射剝離方法處理后地震數(shù)據(jù)計算的20 Hz(圖15c)與50 Hz(圖15d)單頻切片高低頻特征差異明顯，能夠精細刻厚砂體與薄砂體平面展布特征.

圖15 強反射剝離前后砂體平面展布預測效果對比(a)和(b)分別為強反射剝離前20 Hz、50 Hz沿層切片； (c)和(d)分別為強反射剝離后20 Hz、50 Hz沿層切片.Fig.15 Comparison of prediction results of sand-body distribution before and after strong reflection removal(a) and (b) are 20 Hz and 50 Hz horizon slices before strong reflection removing, respectively；(c) and (d) are 20 Hz and 50 Hz horizon slices after strong reflection removing, respectively.

7 結(jié)論與討論

本文提出了一種基于深度學習網(wǎng)絡的地震強反射剝離方法，并用于輔助被強反射壓制的儲層弱反射信號的檢測工作.利用深度學習網(wǎng)絡我們可以全自動地實現(xiàn)端到端的復雜非線性映射，在提升處理效果的基礎上省去專業(yè)處理和解釋人員大量的經(jīng)驗性調(diào)參過程，計算速度快，適用于大規(guī)模三維地震數(shù)據(jù)處理，可以大大提高該類儲層的解釋精度與工作效率.得到如下結(jié)論：

在卷積模型的框架下構(gòu)造的包括地震子波、強反射系數(shù)和弱反射系數(shù)三個未知量的強反射預測反問題的多解性非常強，即便是能夠求解復雜非線性映射問題的深度學習方法也難以直接求解.本文假設地震子波在一定空間范圍內(nèi)時不變，設計先求解地震子波后估計強反射的工作流程，將強反射預測問題分解為地震子波預測與強反射預測兩個端到端的子優(yōu)化問題，并分別采用AIDNN與U-Net兩個深度網(wǎng)絡進行求解.模型資料試算結(jié)果表明，本文提出的方法在不傷害強反射周圍有效反射信號的前提下能夠準確剝離強反射，且保真保幅性好.實際資料的強反射剝離結(jié)果及在該結(jié)果基礎上進行的儲層砂體展布預測結(jié)果在井點處與測井解釋結(jié)果一致.以上效果得益于：(1)假設地震子波在一定空間范圍內(nèi)時不變，將難以求解的強反射預測問題分解為子波預測與強反射預測兩個端到端的優(yōu)化子問題；(2)假設強反射附近的不含強反射的資料與強反射預測時窗內(nèi)的弱反射近似同分布，選擇強反射附近的資料作為弱反射訓練樣本，保證了訓練樣本的代表性且避免了訓練過程中噪聲的影響；(3)通過對訓練樣本與待處理數(shù)據(jù)的預處理，放大強反射特征，使得深度網(wǎng)絡更容易學到強反射的特征.綜上所述，本文提出的方法能夠有效的應用于強反射剝離，這為該類儲層的預測工作提供了重要技術(shù)支撐.

本研究雖然取得了重要進展，但仍然存在以下兩點問題需要進一步探討：(1)訓練樣本：訓練樣本設置的合理性對于深度網(wǎng)絡的訓練效果影響很大，理想數(shù)據(jù)集的每一類別都需要充分的樣本數(shù)量，并且能夠覆蓋地震波形所具有的各種復雜特征；本文在模型數(shù)據(jù)應用部分預測的地震子波并不理想，導致這一結(jié)果的原因有以下幾點：①在卷積模型框架下制作的合成地震記錄訓練樣本不能很好的模擬地震波傳播過程中幾何學與動力學特性；但本文為了增加難度，驗證算法所使用的模型數(shù)據(jù)是采用波動方程正演得到的.②地震子波訓練樣本只選用了雷克子波，但實際資料中子波的基本形態(tài)，它的頻寬，僅用瑞克子波并不能很好地描述.③訓練地震數(shù)據(jù)與地震子波的映射關(guān)系時沒有加入待研究工區(qū)的實際地下地質(zhì)規(guī)律，僅選用了Marmousi模型的反射系數(shù).針對訓練樣本的三點問題，擬提出以下三點探索計劃：①考慮地震波傳播過程中幾何學與動力學特性，在波動方程框架下模擬得到地震記錄訓練樣本.②進一步擴大地震子波樣本庫，在雷克子波的基礎上進一步加入俞氏子波、帶通子波、統(tǒng)計子波等作為地震子波訓練樣本.③結(jié)合測井數(shù)據(jù)設計更加匹配待處理工區(qū)的地質(zhì)模型，并獲得反射系數(shù)訓練樣本.(2)泛化性：本研究假設在一定空間范圍內(nèi)地震子波時不變，且反射系數(shù)同分布；超出一定范圍便不滿足這些假設；因此，處理不同工區(qū)資料時都需要花費時間重新訓練實際地震道與強反射之間的映射關(guān)系；為了使訓練結(jié)果泛化性能更好，應當進一步總結(jié)不同地震資料反射層分布特征，設計更具代表性的樣本，以增加訓練樣本的復雜性及普適性，從而使得所提出算法的訓練結(jié)果普適性更好.

致謝感謝審稿專家對本文提出的寶貴意見和建議，感謝中國石油長慶油田公司勘探開發(fā)研究院提供的實際地震數(shù)據(jù).