李陽何登科

1.中國礦業(yè)大學(xué)(北京)地球科學(xué)與測繪工程學(xué)院,北京 100083;2.煤炭資源與安全開采國家重點(diǎn)實(shí)驗(yàn)室,北京 100083

初至拾取技術(shù)是地震資料處理中的一項(xiàng)關(guān)鍵技術(shù)。 初至拾取的精度直接決定靜校正的精度,并對后續(xù)速度分析、偏移成像的質(zhì)量造成影響。 隨著地震勘探技術(shù)的不斷革新,地震勘探采集的數(shù)據(jù)體量也變得十分龐大。 使用傳統(tǒng)方法如長短時(shí)窗平均能量比法[1]、偏振分析法、時(shí)頻分析法等進(jìn)行初至拾取,需要耗費(fèi)大量時(shí)間進(jìn)行人工干預(yù),工作量巨大[2]。 因此,針對如何進(jìn)行高效的自動(dòng)拾取初至問題,諸多學(xué)者進(jìn)行了大量的討論研究。

近年來,得益于人工智能技術(shù)的迅猛發(fā)展,基于神經(jīng)網(wǎng)絡(luò)的深度學(xué)習(xí)算法被廣泛應(yīng)用到地球物理的各個(gè)領(lǐng)域[3-6],并取得了顯著成效。 在計(jì)算機(jī)視覺領(lǐng)域,初至拾取通常作為語義分割任務(wù)來處理。 全卷積網(wǎng)絡(luò)[7](Fully Convolutional Networks,FCN)已經(jīng)成為圖像語義分割的基本框架,分割效果較傳統(tǒng)神經(jīng)網(wǎng)絡(luò)在精度和速度上有明顯提高。后續(xù)出現(xiàn)的語義分割網(wǎng)絡(luò)如SegNet 網(wǎng)絡(luò)[8]和UNet 網(wǎng)絡(luò)[9],也都是基于FCN。

針對初至自動(dòng)拾取的準(zhǔn)確性、穩(wěn)定性和抗噪性等問題,研究者陸續(xù)提出了多種基于FCN 及其變體的新拾取方法。 這些方法主要分為兩類:一類是通過增加網(wǎng)絡(luò)深度或替換網(wǎng)絡(luò)模塊提高模型的精度,如劉佳楠等[10]構(gòu)建了3 種不同深度的FCN 網(wǎng)絡(luò)進(jìn)行初至拾取;丁建群等[11]提出通過引入空洞卷積改善拾取結(jié)果;楊雙瑜等[12]則通過增加Dense Block,利用深層網(wǎng)絡(luò)強(qiáng)化特征抽取,提高U-Net 在低信噪比環(huán)境下的穩(wěn)定性;陳德武等[13]通過引入殘差連接改進(jìn)U-Net 網(wǎng)絡(luò)的跳躍連接來提高U-Net的準(zhǔn)確率和效率。 另一類方法是通過增加后處理步驟改進(jìn)處理結(jié)果,如David Cova 等[14]提出通過視速度約束提高結(jié)果的連續(xù)性。 以上兩類方法雖然提高了拾取效果,但也增加了模型的復(fù)雜度和訓(xùn)練難度。

基于偏微分方程(Partial Different Equation,PDE)的方法也可用于圖像分割。 水平集方法是由Osher 和Sethian 于1988年提出的一種將二維閉合曲線演化問題轉(zhuǎn)化為高維曲面演化過程進(jìn)行處理的方法[15]。 傳統(tǒng)的水平集方法往往是利用求解能量方程最小化來分割圖像中的目標(biāo)區(qū)域[16],無法直接應(yīng)用到神經(jīng)網(wǎng)絡(luò)。 Kim 等[17]提出了一種水平集損失函數(shù),通過改造水平集方法,使其能應(yīng)用到卷積神經(jīng)網(wǎng)絡(luò)中。 這種基于水平集方法的水平集損失能更好地利用像素間的關(guān)系,而網(wǎng)絡(luò)中常用的損失函數(shù)如二元交叉熵?fù)p失函數(shù),都是獨(dú)立計(jì)算各像素的損失。 對于初至拾取任務(wù),這也是造成模型拾取結(jié)果連續(xù)性不好的一個(gè)原因。

針對地震初至拾取的主要問題,本文結(jié)合水平集方法對語義分割網(wǎng)絡(luò)U-Net 進(jìn)行了改進(jìn),即在原網(wǎng)絡(luò)中增加計(jì)算水平集損失。 網(wǎng)絡(luò)能夠利用地震道間初至的連續(xù)性特征識別地震數(shù)據(jù)體邊界,從而提高拾取結(jié)果的精度。 另外,為了克服訓(xùn)練數(shù)據(jù)不足問題,通過采用人工合成數(shù)據(jù)和實(shí)際數(shù)據(jù)交叉的方式制作訓(xùn)練數(shù)據(jù)集,既簡化了訓(xùn)練數(shù)據(jù)構(gòu)建標(biāo)簽難度,又保證了樣本質(zhì)量。 最后,為了驗(yàn)證該方法的有效性,對比測試了將該方法應(yīng)用到模擬數(shù)據(jù)和實(shí)際地震數(shù)據(jù)上的拾取效果,并對結(jié)果進(jìn)行了分析說明。

1 方法原理

1.1 U-Net 網(wǎng)絡(luò)

U-Net 網(wǎng)絡(luò)是一種經(jīng)典的基于FCN 的語義分割網(wǎng)絡(luò),其核心設(shè)計(jì)包括編碼器-解碼器結(jié)構(gòu)和跳躍連接[9]。 U-Net 網(wǎng)絡(luò)結(jié)構(gòu)呈U 型對稱,左側(cè)為編碼器,右側(cè)為解碼器。 編碼器實(shí)現(xiàn)圖像空間特征提取并對結(jié)果進(jìn)行下采樣,由4 個(gè)子模塊構(gòu)成,每個(gè)模塊包含2 個(gè)3×3 卷積操作和1 個(gè)2×2 的池化操作;解碼器對圖像特征進(jìn)行上采樣生成分割結(jié)果圖,與編碼器結(jié)構(gòu)相似,也由4 個(gè)子模塊構(gòu)成,每個(gè)模塊包含2 個(gè)3×3 卷積操作和1 個(gè)2×2 的反卷積操作。 跳躍連接則將解碼器上采樣結(jié)果與編碼器同層的具有相同分辨率的輸出進(jìn)行拼接,實(shí)現(xiàn)不同尺度圖像特征的融合[18],以提高圖像分割的準(zhǔn)確度。 網(wǎng)絡(luò)的最終輸出結(jié)果,通過1 個(gè)1×1 的卷積操作生成最終的分割結(jié)果圖。

U-Net 網(wǎng)絡(luò)的卷積塊使用修正線性單元(Rectified Linear Unit,ReLU)激活函數(shù),形成如下形式:

式(1)計(jì)算簡單,不存在梯度消失問題。 但ReLU 激活函數(shù)會(huì)造成部分神經(jīng)元失效,使網(wǎng)絡(luò)變得稀疏。 近年來出現(xiàn)一些ReLU 激活函數(shù)的改型,如隨機(jī)滲漏整流線性單元(Randomized Leaky ReLU,RReLU)、指數(shù)線性單元(Exponential Linear Unit,ELU)等,可克服ReLU 的不足。

1.2 水平集方法

使用水平集方法的卷積語義分割網(wǎng)絡(luò)[17],通過利用學(xué)習(xí)目標(biāo)的空間信息提升網(wǎng)絡(luò)學(xué)習(xí)性能,從而獲得了更加平滑、精確的分割效果。 而David Cova 等[14]設(shè)計(jì)的對比實(shí)驗(yàn)發(fā)現(xiàn),通過優(yōu)化超參數(shù)而不改變U-Net 網(wǎng)絡(luò)基本架構(gòu)也能夠獲得令人滿意的結(jié)果。 將水平集方法應(yīng)用到U-Net 網(wǎng)絡(luò),可使U-Net 學(xué)習(xí)到地震道之間初至的連續(xù)性特征。

水平集方法通常用于分割前景和背景,多分類問題可以采用通道的方式轉(zhuǎn)換為多個(gè)二分類任務(wù)處理。 針對地震初至拾取任務(wù),為簡化網(wǎng)絡(luò),水平集方法處理中將其視作二分類任務(wù)。 本文采用的能量方程[17]形式如下:

式中,φ為符號距離函數(shù);G為待分割地震圖像的目標(biāo)區(qū)域;L為分割邊界集,即初至點(diǎn)集;Cl,1、Cl,2分別為分割邊界內(nèi)部和外部的灰度均值;Hε為階躍函數(shù)的近似表示;ε為常量,ε→0。

式(2)作為損失函數(shù),是可微的,能夠應(yīng)用于網(wǎng)絡(luò)進(jìn)行反向傳播更新參數(shù)。

1.3 改進(jìn)后的U-Net 網(wǎng)絡(luò)

應(yīng)用水平集方法后,改進(jìn)的U-Net 網(wǎng)絡(luò)結(jié)構(gòu)如圖1 所示。 首先將訓(xùn)練數(shù)據(jù)分批送入U(xiǎn)-Net 網(wǎng)絡(luò)模型中,經(jīng)過下采樣和上采樣步驟提取的融合圖像特征信息,經(jīng)由U-Net 網(wǎng)絡(luò)輸出一個(gè)像素分類的概率預(yù)測圖。 訓(xùn)練數(shù)據(jù)的標(biāo)簽使用3 種分類標(biāo)記,即背景、初至波和續(xù)至波。 由于采用3 種分類方式,水平集方法通常用于分割目標(biāo)和背景,因此需要合并計(jì)算初至和地震數(shù)據(jù),即水平集方法用于分割有效地震數(shù)據(jù)和噪聲數(shù)據(jù)。 將這些預(yù)測的結(jié)果代入式(2)中,就可以求出水平集能量泛函的項(xiàng)。 水平集損失通過式(3)獲得:

圖1 網(wǎng)絡(luò)結(jié)構(gòu)Fig.1 Network structure diagram

式中,yi為真實(shí)值;ei為水平集預(yù)測值。

由于訓(xùn)練樣本中初至點(diǎn)樣點(diǎn)比較少,90% 以上是標(biāo)簽為噪聲或地震數(shù)據(jù)的非初至樣點(diǎn),因此初至拾取存在類別不平衡、難易樣本不均衡問題。 針對該問題,損失選用Focal Loss 函數(shù)效果更好[20]。Focal Loss 形式如下:

式中,α,γ為權(quán)值;yi為真實(shí)值;pi為網(wǎng)絡(luò)預(yù)測值。

將U-Net 輸出的預(yù)測結(jié)果值與輸入數(shù)據(jù)的標(biāo)簽數(shù)據(jù)值代入到式(4)中,計(jì)算兩者的Focal Loss 損失。引入水平集方法之后,最終的損失函數(shù)定義為

式中,β為權(quán)值。

將兩種損失代入式(5)獲得總的損失,最后進(jìn)行誤差反向傳播,進(jìn)一步優(yōu)化校正網(wǎng)絡(luò)權(quán)重參數(shù)。

為了降低網(wǎng)絡(luò)模型參數(shù),抑制不相關(guān)背景區(qū)域的特征響應(yīng),引入注意力機(jī)制模塊[19],其結(jié)構(gòu)如圖2 所示。 首先將特征矩陣經(jīng)過1 個(gè)2×2 的卷積層、歸一化和激活函數(shù)Sigmoid,再與初始特征矩陣相乘,計(jì)算注意力向量提升特征學(xué)習(xí),可以在多尺度采樣后進(jìn)一步細(xì)化小像素點(diǎn)初至波的特征結(jié)構(gòu)[3]。

圖2 注意力機(jī)制模塊結(jié)構(gòu)[19]Fig.2 Attention model structure[19]

1.4 訓(xùn)練集制作方法

訓(xùn)練數(shù)據(jù)集通常采用人工標(biāo)注制作標(biāo)簽,工作量大,效率不高。 為了解決這個(gè)問題,本文在訓(xùn)練數(shù)據(jù)集的標(biāo)簽構(gòu)建過程中,先采用人工合成地震記錄制作訓(xùn)練數(shù)據(jù)集,訓(xùn)練模型;再使用訓(xùn)練后的模型處理實(shí)際地震數(shù)據(jù)集,獲得初步的初至記錄;最后,實(shí)際地震數(shù)據(jù)樣本經(jīng)人工校正后,和合成地震記錄混合制作訓(xùn)練數(shù)據(jù)集。

為了便于獲得地震記錄的初至,采用褶積模型生成模擬地震數(shù)據(jù),同時(shí)記錄初至到時(shí)數(shù)據(jù)。 地震記錄褶積模型的時(shí)間域表達(dá)式如下:

式中,s(t)為有效波;ω(t)為地震子波;r(t)為反射系數(shù);n(t)為隨機(jī)環(huán)境噪聲。

在進(jìn)行合成前,先根據(jù)實(shí)際工區(qū)地層情況,設(shè)置地震波速度和地層數(shù)量隨機(jī)范圍,隨機(jī)生成地層的地震波速度模型。 然后,根據(jù)反射系數(shù)和波阻抗的對應(yīng)關(guān)系計(jì)算得到反射系數(shù),再分別與不同頻率的Ricker 子波褶積得到相應(yīng)的人工合成地震數(shù)據(jù)。

訓(xùn)練數(shù)據(jù)集構(gòu)建標(biāo)簽使用三分類方式,初至上方像素標(biāo)記為0,將初至上下4 像素標(biāo)記為初至波值為1,初至下方地震數(shù)據(jù)體標(biāo)記為2。 如果采用二值分類,目標(biāo)區(qū)域和背景區(qū)域像素?cái)?shù)量相差不大,初至位置的細(xì)微改變對損失函數(shù)的影響不明顯。 而U-Net 輸出的結(jié)果代表著每個(gè)像素劃分為目標(biāo)區(qū)域或者背景區(qū)域的概率,該結(jié)果本身存在著一定的誤差,最終造成網(wǎng)絡(luò)靈敏度不夠,降低分割算法的準(zhǔn)確率。

2 數(shù)據(jù)準(zhǔn)備與模型訓(xùn)練

2.1 制作人工合成地震記錄

首先使用褶積模型制作人工合成地震記錄。為了使模擬數(shù)據(jù)更接近實(shí)際數(shù)據(jù)情況,對其添加不同級別的高斯白噪聲,生成不同信噪比的模擬地震數(shù)據(jù)。 所使用信噪比公式可表示為

式中,s(t)為原始信號;n(t)為高斯白噪聲。

按照褶積模型合成的地震記錄如圖3 所示。 圖3(a)為不含噪聲的地震記錄,圖3(b)為添加隨機(jī)噪聲的合成地震記錄,藍(lán)色標(biāo)記為初至數(shù)據(jù)記錄。

圖3 合成地震數(shù)據(jù)Fig.3 Synthetic data

2.2 實(shí)際地震數(shù)據(jù)及預(yù)處理

實(shí)際地震數(shù)據(jù)為某礦區(qū)三維地震資料,包含243 個(gè)炮集。 道間距為20 m,最大炮檢距為3 200 m,每個(gè)炮集包含475 道,采樣間隔為2 ms。

為保證網(wǎng)絡(luò)學(xué)習(xí)效果,在學(xué)習(xí)前對實(shí)際地震數(shù)據(jù)進(jìn)行必要的數(shù)據(jù)預(yù)處理。

(1) 去噪處理。 為了保留圖像邊緣細(xì)節(jié)信息,采用中值濾波算法處理地震數(shù)據(jù)。 中值濾波作為一種非線性濾波技術(shù),在一定條件下可以克服線性濾波帶來的信號邊界模糊等問題[22]。

(2) 振幅補(bǔ)償處理。 使用球面擴(kuò)散補(bǔ)償和地層吸收能量補(bǔ)償,對地震振幅進(jìn)行補(bǔ)償和校正。

(3) 歸一化處理。 為了避免振幅數(shù)據(jù)值對訓(xùn)練結(jié)果的影響,在訓(xùn)練網(wǎng)絡(luò)之前需要對數(shù)據(jù)進(jìn)行歸一化。 針對一道地震數(shù)據(jù){x1,x2,…,xm},歸一化計(jì)算公式如下:

式中,xi為該道當(dāng)前樣點(diǎn)的值;xmin為該道所有樣點(diǎn)振幅的最小值;xmax為該道所有樣點(diǎn)振幅的最大值。

2.3 制作訓(xùn)練數(shù)據(jù)集

使用人工合成地震記錄訓(xùn)練U-Net 網(wǎng)絡(luò),生成初始模型。 應(yīng)用初始模型對預(yù)處理后的實(shí)際地震數(shù)據(jù)進(jìn)行拾取,拾取結(jié)果如圖4 所示。

圖4 初始模型拾取結(jié)果Fig.4 The results picked by the initial model

使用地震資料處理軟件校正初至記錄,并更新到樣本對應(yīng)的道頭數(shù)據(jù)庫。 校正完成后,再將樣本及其初至數(shù)據(jù)保存至文件,完成實(shí)際地震數(shù)據(jù)集的初至拾取。 數(shù)據(jù)及對應(yīng)標(biāo)簽如圖5 所示。

圖5 地震數(shù)據(jù)及對應(yīng)標(biāo)簽數(shù)據(jù)Fig.5 Seismic data and labeled data

根據(jù)實(shí)際地震資料情況,每39 道合成1 個(gè)訓(xùn)練數(shù)據(jù)。 隨機(jī)選取人工合成地震數(shù)據(jù)和實(shí)際地震數(shù)據(jù)。 訓(xùn)練數(shù)據(jù)集共1 000 個(gè)訓(xùn)練數(shù)據(jù),其中60% 用作訓(xùn)練集,20% 用作測試集,20% 用作驗(yàn)證集。

2.4 訓(xùn)練環(huán)境與模型訓(xùn)練

網(wǎng)絡(luò)模型使用Python 語言基于PyTorch 深度學(xué)習(xí)框架實(shí)現(xiàn),并在CPU 為Intel Xeon Silver 4114、主頻2.2 GHz、內(nèi)存為64 G 的計(jì)算機(jī)上訓(xùn)練。 在訓(xùn)練過程中,采用小批量梯度下降方法,即把整個(gè)訓(xùn)練數(shù)據(jù)集隨機(jī)分為若干個(gè)批次,每次迭代僅使用其中的1 個(gè)批次的樣本參與參數(shù)更新。 為了獲得該網(wǎng)絡(luò)最優(yōu)的結(jié)果,分別對激活函數(shù)(Activation Function)、優(yōu)化器(Optimizer)、學(xué)習(xí)率(Learning Rate)進(jìn)行超參數(shù)優(yōu)化,所使用方法和值見表1。

表1 超參數(shù)優(yōu)化參數(shù)Table 1 Hyperparametric values

采用網(wǎng)絡(luò)訓(xùn)練的權(quán)重(Weight) 系數(shù)設(shè)為(0.005,0.015,0.98)。 訓(xùn)練的初始學(xué)習(xí)率設(shè)為0.01,epoch 設(shè)為500,小批量梯度下降的批次大小(Batch Size)設(shè)置為20。 每輪迭代50 步,每輪用驗(yàn)證集驗(yàn)證1 次。 訓(xùn)練過程耗時(shí)6 h 29 min,第140次迭代時(shí),網(wǎng)絡(luò)逐漸收斂。

3 模型應(yīng)用

3.1 模擬數(shù)據(jù)拾取

為了檢驗(yàn)?zāi)Ｐ驮诓煌肼晽l件下的表現(xiàn),使用Devito 框架[23]通過二階交錯(cuò)網(wǎng)格有限差分波動(dòng)方程正演獲得地震模擬數(shù)據(jù)進(jìn)行測試。 地層模型和生成的地震數(shù)據(jù)如圖6 所示,震源使用Ricker 子波,采用2 ms 采樣,地震記錄長度3 s,道間距25 m,共201 道。

圖6 地層模型與模擬數(shù)據(jù)Fig.6 Stratigraphic model and simulated data

應(yīng)用U-Net 網(wǎng)絡(luò)模型分別對無噪聲和添加了不同信噪比的隨機(jī)噪聲模擬地震數(shù)據(jù)進(jìn)行初至拾取測試,拾取結(jié)果如圖7 所示。 可以看出,改進(jìn)后的網(wǎng)絡(luò)在較低信噪比的模擬數(shù)據(jù)上拾取結(jié)果與無噪聲數(shù)據(jù)基本重合,在信噪比較低的情況下,拾取結(jié)果出現(xiàn)較大偏差,但在近中炮檢距的范圍內(nèi)拾取結(jié)果依然準(zhǔn)確。

圖7 不同信噪比下模擬數(shù)據(jù)拾取結(jié)果Fig.7 First break picking results of simulated data with different SNR noise

3.2 實(shí)際地震數(shù)據(jù)拾取

為了驗(yàn)證改進(jìn)的網(wǎng)絡(luò)應(yīng)用于實(shí)際地震數(shù)據(jù)的適應(yīng)性,使用未參與訓(xùn)練的一炮地震數(shù)據(jù)進(jìn)行初至拾取。 圖8 所示為實(shí)際地震數(shù)據(jù)的拾取結(jié)果。 為了清楚展示驗(yàn)證效果,單獨(dú)提取出中間600 ～790道數(shù)據(jù)的拾取結(jié)果。 完整炮集拾取結(jié)果見表2?？梢钥闯?本方法應(yīng)用于實(shí)際礦區(qū)地震勘探數(shù)據(jù)時(shí),結(jié)果連續(xù)性相對較好,并可以獲得精度較高的拾取結(jié)果。

圖8 實(shí)際數(shù)據(jù)拾取結(jié)果Fig.8 The picking result using actual data

表2 實(shí)際數(shù)據(jù)拾取結(jié)果統(tǒng)計(jì)Table 2 The statistical result picked using actual data

為了驗(yàn)證模型在復(fù)雜地形條件下的性能,選取地表起伏較大的某礦區(qū)實(shí)際地震數(shù)據(jù)進(jìn)行測試。實(shí)際拾取結(jié)果如圖9 所示,圖中藍(lán)色標(biāo)記為初至拾取結(jié)果,中間子圖為檢波器高程曲線。 檢波器最大高差20 m,整體與實(shí)際初至基本吻合。

圖9 拾取結(jié)果與相對高程Fig.9 The picking result and relative elevation

4 結(jié) 論

為解決使用U-Net 網(wǎng)絡(luò)拾取地震初至存在的拾取結(jié)果不連續(xù)等問題,在不增加網(wǎng)絡(luò)復(fù)雜度的情況下,本文結(jié)合U-Net 和水平集方法提出了一種網(wǎng)絡(luò)模型,利用地震道間初至的連續(xù)性特征識別地震數(shù)據(jù)體邊界,并進(jìn)行了模擬數(shù)據(jù)與實(shí)際地震數(shù)據(jù)測試。 通過驗(yàn)證有如下結(jié)論:

(1) 通過采用水平集方法可以較好地解決UNet 網(wǎng)絡(luò)拾取結(jié)果不連續(xù)問題。

(2) 通過采用Focal Loss 損失函數(shù)可以改善訓(xùn)練集樣本類別不均衡問題。

(3) 當(dāng)算法不依賴地層信息時(shí),可以通過采用預(yù)拾取和混合訓(xùn)練數(shù)據(jù)的方法,緩解實(shí)際數(shù)據(jù)標(biāo)簽集不易制作和訓(xùn)練集數(shù)據(jù)不足的問題。

(4) 對于較低信噪比的地震數(shù)據(jù),應(yīng)用結(jié)合U-Net 和水平集方法的改進(jìn)模型,相對于原始UNet 模型能較準(zhǔn)確地拾取地震波初至。