基于空間—波數(shù)域聯(lián)合深度學(xué)習(xí)的數(shù)值頻散壓制

張 巖 崔淋淇 宋利偉 董宏麗

(①東北石油大學(xué)計(jì)算機(jī)與信息技術(shù)學(xué)院，黑龍江大慶 163318； ②東北石油大學(xué)物理與電子工程學(xué)院，黑龍江大慶 163318； ③東北石油大學(xué)人工智能能源研究院，黑龍江大慶 163318； ④黑龍江省網(wǎng)絡(luò)與智能控制重點(diǎn)實(shí)驗(yàn)室，黑龍江 163318)

0 引言

地震波傳播的正演數(shù)值模擬是地震勘探開發(fā)中重要研究的課題。波動方程法能完整表征波在傳播過程中振幅、相位、頻率的變化，同時(shí)保持了地震波的運(yùn)動學(xué)與動力學(xué)特征，被廣泛應(yīng)用于地震正演模擬。有限差分法由于適用于大規(guī)模并行計(jì)算和易編程等優(yōu)勢，是求解波動方程的常用方法。然而用離散的差分算子代替連續(xù)的微分算子，產(chǎn)生的誤差會存在于整個(gè)數(shù)值模擬記錄中，導(dǎo)致數(shù)值頻散，影響正演模擬精度和逆時(shí)偏移成像效果。因此，頻散壓制是提高數(shù)值模擬精度的重要手段。目前數(shù)值頻散壓制方法從原理上大體分為以下三種。

(1)從數(shù)值頻散產(chǎn)生的機(jī)理出發(fā)，使用較長的差分算子或精細(xì)的差分網(wǎng)格進(jìn)行頻散壓制。Alford等[1]和Dablain[2]對聲波二階差分?jǐn)?shù)值頻散進(jìn)行分析，指出網(wǎng)格尺寸和地震波傳播方向是影響頻散的兩個(gè)因素； 董良國等[3-4]指出高階有限差分法可以較好解決數(shù)值頻散問題； 裴正林等[5]利用一階速度—應(yīng)力彈性波方程在計(jì)算域內(nèi)采用高階有限差分法高精度模擬了彈性波場； Boore[6]的研究表明，當(dāng)震源信號波長與網(wǎng)格間距的比值超過10時(shí)，頻散顯著減弱； 孫林潔等[7]根據(jù)頻散關(guān)系改變網(wǎng)格剖分方式，提高了有限差分法模擬的精度。此類方法會使計(jì)算復(fù)雜度和內(nèi)存需求顯著增加。

(2)通量校正傳輸(Flux-corrected Transport，F(xiàn)CT)法。Fei等[8]首次在彈性波的數(shù)值模擬中運(yùn)用FCT技術(shù)，發(fā)現(xiàn)該方法對數(shù)值頻散有較好壓制效果； Book等[9]在求解聲波方程時(shí)應(yīng)用FCT法有效壓制了數(shù)值頻散； 中國學(xué)者也對FCT法消除數(shù)值頻散現(xiàn)象方面進(jìn)行了研究、分析[10-13]。FCT技術(shù)壓制數(shù)值頻散的基本思想是先認(rèn)為全部的極值點(diǎn)是數(shù)值頻散形成的，然后通過對全部的網(wǎng)格節(jié)點(diǎn)都實(shí)行漫射校正，從而實(shí)現(xiàn)數(shù)值頻散的壓制。此類方法雖然減輕了計(jì)算量，但會降低有效信號的分辨率，造成振幅損失。

(3)優(yōu)化波動方程或差分系數(shù)。Etgen[14]使用高斯—牛頓迭代方法確定有限差分系數(shù)； 張志禹等[15]在常規(guī)差分方程中增加了頻散校正項(xiàng)，有效衰減高波數(shù)成分、抑制頻散； Liu[16]應(yīng)用最小二乘法獲取給定波數(shù)范圍內(nèi)的優(yōu)化有限差分系數(shù)。還有許多優(yōu)化差分系數(shù)求取方法[17-21]，在壓制數(shù)值頻散方面取得了較好的效果。此類方法是多參數(shù)優(yōu)化問題，需要憑借豐富的經(jīng)驗(yàn)選取所優(yōu)化的參數(shù)值，用較長計(jì)算時(shí)間確定系數(shù)且不能保證所得參數(shù)為全局最優(yōu)解。

上述頻散壓制方法在一定精度要求下，往往導(dǎo)致計(jì)算復(fù)雜度大幅度增加，耗時(shí)過長，因此探究新的數(shù)值頻散壓制方法具有重要意義。隨著硬件能力的提高和大數(shù)據(jù)技術(shù)的成熟，深度學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)廣泛應(yīng)用到各個(gè)領(lǐng)域。深度學(xué)習(xí)可以從大數(shù)據(jù)中自動學(xué)習(xí)樣本特征，從而表現(xiàn)出顯著優(yōu)勢，在地球物理領(lǐng)域也得到了廣泛的應(yīng)用[22-28]。Kaur等[29]首次提出利用cycleGAN網(wǎng)絡(luò)去除有限差分波外推中的頻散偽影。該方法通過深度學(xué)習(xí)創(chuàng)建一個(gè)傳播后置頻散校正濾波器，能夠以自適應(yīng)的方式消除頻散偽影，且與傳統(tǒng)方法相比效果明顯較好，但GAN網(wǎng)絡(luò)復(fù)雜度較高，訓(xùn)練難度大，同時(shí)僅考慮到了空間域特征提取，忽略了波場數(shù)據(jù)的稀疏特性，波場的主要特征提取不充分。

本文結(jié)合波場數(shù)據(jù)的特點(diǎn)，參考卷積神經(jīng)網(wǎng)絡(luò)在地震資料處理方面的成功應(yīng)用，提出了一種基于卷積神經(jīng)網(wǎng)絡(luò)壓制頻散的網(wǎng)絡(luò)模型。該模型主要特征如下：

(1)對不同速度模型進(jìn)行正演模擬，分別使用低階和高階有限差分對波動方程求解，得到大量含頻散輸入數(shù)據(jù)集和干凈波場數(shù)據(jù)集作為網(wǎng)絡(luò)的訓(xùn)練輸入和標(biāo)簽；

(2)構(gòu)建聯(lián)合學(xué)習(xí)網(wǎng)絡(luò)模型，通過空間—波數(shù)域聯(lián)合損失誤差約束進(jìn)行多維度的特征提取和參數(shù)學(xué)習(xí)，增加描述波場特征的準(zhǔn)確性，提高網(wǎng)絡(luò)模型的波場特征學(xué)習(xí)能力；

(3)考慮到波場數(shù)據(jù)的稀疏特征，在網(wǎng)絡(luò)模型中引入稀疏正則化約束，將頻散壓制問題轉(zhuǎn)化為L1范數(shù)約束的最優(yōu)化問題求解，降低網(wǎng)絡(luò)模型復(fù)雜度，提高網(wǎng)絡(luò)模型的泛化能力；

(4)利用卷積神經(jīng)網(wǎng)絡(luò)壓制數(shù)值頻散，避免傳統(tǒng)方法的計(jì)算困難和復(fù)雜度，自適應(yīng)提取波場數(shù)據(jù)淺層的像素級特征和深層的語義級特征，與FCT方法相比，無需人為設(shè)定漫射系數(shù)等相關(guān)參數(shù)，技術(shù)人員無需具有豐富的地震波動力學(xué)理論基礎(chǔ)，可直接應(yīng)用預(yù)訓(xùn)練好的模型，實(shí)現(xiàn)“端到端”的處理模式。

本文首先通過均勻速度模型驗(yàn)證網(wǎng)絡(luò)內(nèi)部各結(jié)構(gòu)的合理性； 其次，通過Marmousi模型測試該算法在復(fù)雜模型中的頻散壓制能力； 最后，通過與遷移學(xué)習(xí)結(jié)合的方式驗(yàn)證該算法對于新模型正演數(shù)據(jù)的頻散壓制泛化性能。結(jié)果表明所提算法能夠保護(hù)有效信號，降低計(jì)算成本，有效地壓制有限差分正演模擬過程中產(chǎn)生的數(shù)值頻散。

1 數(shù)值頻散分析

1.1 數(shù)值頻散

在有限差分法求解波動方程的過程中，用離散的差分算子代替連續(xù)的微分算子，使波傳播的相速度減小，不同波數(shù)分量的波場以不同傳播速度進(jìn)行傳播，波前形狀發(fā)生變化并逐漸擴(kuò)散導(dǎo)致波形失真，造成頻散假象，影響數(shù)值模擬的精度。

以各向同性介質(zhì)中的二維波動方程

(1)

為例分析數(shù)值頻散產(chǎn)生的機(jī)制。式中：U為隨空間位置和時(shí)間變化的波場；v為速度；t為時(shí)間；x和z分別為空間上的水平、垂直方向坐標(biāo)。

式(1)的平面波解可表示為

p(x,z,t)=exp[i(ωt-kxcosφ-kzsinφ)]

(2)

式中：k為波數(shù)；φ為入射角；ω=kv0為角頻率，其中v0為數(shù)值相速度。通過推導(dǎo)可得v0與v的關(guān)系為

(3)

式中：M為有限差分階次的一半；cm為有限差分系數(shù)；h=Δx=Δz，為網(wǎng)格間距。

以均勻介質(zhì)中的單炮波場快照為例，監(jiān)測區(qū)域橫、縱方向均設(shè)置50個(gè)網(wǎng)格，地下介質(zhì)速度為2000m/s，在區(qū)域中心選取主頻為40Hz的Ricker子波作震源，網(wǎng)格尺寸為8m×8m，不同精度差分算子模擬的1000ms波場快照如圖1所示。由于差分階數(shù)較低，數(shù)值頻散造成的干擾波逐漸向有效波場內(nèi)部擴(kuò)散，使波場快照的分辨率和信噪比降低(圖1a)； 隨著差分算子精度的提高，數(shù)值頻散逐漸減弱(圖1b、圖1c)。地震有效波形所占檢測區(qū)域的比重較少，因此在空間域內(nèi)波場數(shù)據(jù)具有一定的稀疏性。

圖2為圖1在x=200m處抽取的單道波形數(shù)據(jù)。低階差分算子模擬的波形曲線中出現(xiàn)嚴(yán)重的振鈴效應(yīng)，干擾了有效波形信息； 隨著差分算子精度的提高，所得波形趨于光滑，數(shù)值頻散造成的偽振動被壓制，更接近于輸入的Ricker子波波形。可見通過提高差分算子精度的方式可有效壓制波場模擬中的數(shù)值頻散。

圖1 均勻模型不同精度差分算子模擬的波場快照(a)2階； (b)4階； (c)12階

1.2 波場數(shù)據(jù)的波數(shù)域特征

傅里葉變換是地震數(shù)據(jù)分析的主要數(shù)學(xué)工具。對有限差分法波場數(shù)值模擬結(jié)果，利用傅里葉變換將波場快照從空間域變換至波數(shù)域，降低了波場信號的相關(guān)性，能量集中在少數(shù)波數(shù)成分上，實(shí)現(xiàn)對信號的稀疏表示，使得在波數(shù)域中更易找出數(shù)值頻散的模式。二維傅里葉變換的形式為

(4)

圖3為不同階次差分算子模擬的波數(shù)域波場。由于波數(shù)較高的地震波相速度較慢，滯后于群速度形成“波拖尾”現(xiàn)象，因此當(dāng)不同差分精度算子所得波場快照變換至波數(shù)域時(shí)，含頻散波場數(shù)據(jù)傅里葉變換域內(nèi)的系數(shù)向高波數(shù)分量擴(kuò)散，而干凈波場數(shù)據(jù)波數(shù)域內(nèi)信息分布則更為集中。

圖2 不同階次差分算子模擬波場的單道對比

圖3 2(a)、4(b)、12階(c)差分算子模擬的波場快照波數(shù)域特征

2 數(shù)值頻散壓制網(wǎng)絡(luò)模型

2.1 網(wǎng)絡(luò)模型

本文應(yīng)用殘差學(xué)習(xí)的思想設(shè)計(jì)基于監(jiān)督學(xué)習(xí)的深度卷積網(wǎng)絡(luò)。訓(xùn)練數(shù)據(jù)以輸入—標(biāo)簽對的形式輸入網(wǎng)絡(luò)

{(x1,y1),(x2,y2),…,(xNd,yNd)}

(5)

式中：Nd為波場數(shù)據(jù)總量；x為含頻散波場數(shù)據(jù)；y為對應(yīng)的干凈波場數(shù)據(jù)，即標(biāo)簽數(shù)據(jù)。通過殘差學(xué)習(xí)策略提取數(shù)值頻散產(chǎn)生的偽波紋理及分散的特征。神經(jīng)網(wǎng)絡(luò)的輸出為

(6)

網(wǎng)絡(luò)模型結(jié)構(gòu)如圖4所示。本文網(wǎng)絡(luò)結(jié)構(gòu)由多層的卷積單元(Convolutional Layer，Conv)、批歸一化(Batch Normalization，BN)單元以及激活函數(shù)單元組成。選擇線性整流(Rectifier Linear Unit，ReLU)單元作為激活函數(shù)進(jìn)行非線性映射，可表示為

f(l)=max(0,l)

(7)

式中l(wèi)表示激活函數(shù)的輸入。ReLU函數(shù)具有單側(cè)飽和特性，在負(fù)半?yún)^(qū)的導(dǎo)數(shù)為零，減少激活的矩陣中非0元素個(gè)數(shù)，從而在網(wǎng)絡(luò)內(nèi)部引入稀疏性，避免無用信息的干擾。

第1層為輸入層，由Conv +ReLU單元組成，其中Conv用于提取含頻散數(shù)據(jù)的特征，ReLU用于執(zhí)行非線性映射，去除含頻散數(shù)據(jù)中的冗余特征； 第2～第19層為網(wǎng)絡(luò)的隱藏層，由交替的Conv +BN+ReLU單元組成，用于完成數(shù)據(jù)處理后特征提取的非線性映射； 前20層卷積核大小設(shè)為3×3，步長為1，每一層卷積后得到64張?zhí)卣鲌D； 最后一層為Conv單元，將上一層完成的非線性映射輸出作為輸入進(jìn)行特征提取，在前向傳遞過程中每進(jìn)行一次卷積后均進(jìn)行補(bǔ)零操作，保證訓(xùn)練過程中每個(gè)特征圖與輸入數(shù)據(jù)尺寸一致。

圖5為網(wǎng)絡(luò)壓制頻散的過程，在空間域訓(xùn)練含頻散波場數(shù)據(jù)后，對輸出的波場數(shù)據(jù)與標(biāo)簽分別進(jìn)行傅里葉變換，得到波場數(shù)據(jù)的波數(shù)域特征； 對空間域和波數(shù)域上的波場信息分別進(jìn)行誤差計(jì)算后得到聯(lián)合誤差函數(shù)，通過使聯(lián)合誤差函數(shù)最小化進(jìn)行網(wǎng)絡(luò)參數(shù)更新。

2.2 稀疏正則化聯(lián)合損失函數(shù)的建立

2.2.1 聯(lián)合損失函數(shù)

針對頻散特性，使用聯(lián)合損失函數(shù)從空間—波數(shù)聯(lián)合域上對預(yù)測波場與目標(biāo)波場的差異進(jìn)行多維度描述，聯(lián)合損失函數(shù)的形式為

(8)

(9)

式中sign(·)為符號函數(shù)。由于L1范數(shù)在0處不可導(dǎo)，設(shè)置sign(0)=0，代替不連續(xù)區(qū)域的梯度值。

圖4 網(wǎng)絡(luò)模型結(jié)構(gòu)

圖5 網(wǎng)絡(luò)頻散壓制流程

2.2.2 稀疏約束正則項(xiàng)

本文聯(lián)合損失函數(shù)僅可保證模型取得較小的訓(xùn)練誤差，網(wǎng)絡(luò)的復(fù)雜程度和泛化能力需要引入稀疏約束項(xiàng)。稀疏約束項(xiàng)可以理解為對網(wǎng)絡(luò)復(fù)雜度的約束，以減少網(wǎng)絡(luò)中的參數(shù)冗余，提高模型的泛化性能。由于波場數(shù)據(jù)在空間域和波數(shù)域均具有稀疏性，樣本中有效信號的比重較小，引入稀疏正則約束后則可以較好收縮模型、減小解空間，提高網(wǎng)絡(luò)預(yù)測的準(zhǔn)確性和魯棒性。因此本文通過在聯(lián)合損失函數(shù)中加入模型的稀疏約束，達(dá)到利用先驗(yàn)信息壓制波場頻散的目的。稀疏約束項(xiàng)為

(10)

L0正則項(xiàng)從直觀上是對網(wǎng)絡(luò)非零參數(shù)的個(gè)數(shù)進(jìn)行約束，可實(shí)現(xiàn)網(wǎng)絡(luò)參數(shù)對模型的有效表達(dá)得到最稀疏的解，但L0范數(shù)非凸且不可微，很難對其進(jìn)行優(yōu)化求解。L1正則化是L0正則化的最優(yōu)凸近似，L1范數(shù)作為懲罰項(xiàng)在訓(xùn)練過程中對權(quán)重的更新進(jìn)行稀疏約束，使模型的大量參數(shù)近似為0，實(shí)現(xiàn)網(wǎng)絡(luò)模型參數(shù)的稀疏化和特征選擇優(yōu)化。對于卷積核的權(quán)重θ′={θ′1,θ′2,…，θ′J}，其中J為卷積核中參數(shù)的數(shù)量，稀疏約束表達(dá)式為

(11)

則添加正則約束的最終損失函數(shù)為

(12)

式中λ為稀疏約束項(xiàng)系數(shù)。當(dāng)λ→0時(shí)，表示不對網(wǎng)絡(luò)權(quán)值進(jìn)行約束，只通過樣本數(shù)據(jù)學(xué)習(xí)權(quán)值； 當(dāng)λ→∞時(shí)，表示只根據(jù)懲罰項(xiàng)確定網(wǎng)絡(luò)的權(quán)值。

為了平衡空間域與波數(shù)域損失的權(quán)重，設(shè)置μ∶v=1∶5，λ=5×10-5，采用自適應(yīng)矩估計(jì)(Adaptive Moment Estimation，Adam)算法對損失函數(shù)進(jìn)行最優(yōu)化學(xué)習(xí)。

3 模型實(shí)驗(yàn)及結(jié)果

分別選取均勻速度模型以及Marmousi速度模型對本文算法進(jìn)行訓(xùn)練和測試。對于每一個(gè)速度模型，設(shè)置不同的網(wǎng)格尺寸、主頻、采樣間隔等正演參數(shù)，選取Ricker子波作為震源，分別采用低階差分和12階差分進(jìn)行波場模擬，所得到的波場數(shù)據(jù)作為含頻散的訓(xùn)練數(shù)據(jù)集和標(biāo)簽數(shù)據(jù)集。

頻散壓制效果可采用信噪比(SNR)衡量，其數(shù)學(xué)定義為

(13)

式中Nx、Ny為波場切片的x和z方向上的樣點(diǎn)數(shù)。

實(shí)驗(yàn)選用GPU為NVIDA GTX-2080的平臺，操作系統(tǒng)為64位Ubuntu 18.04，軟件環(huán)境為Python 3.8，聯(lián)合深度學(xué)習(xí)模型使用Pytorch1.6框架搭建，CUDA版本為11.0。

3.1 均勻速度模型實(shí)驗(yàn)

為了驗(yàn)證模型的泛化性能，通過不同正演模擬參數(shù)數(shù)值模擬構(gòu)建數(shù)據(jù)集A：監(jiān)測區(qū)域橫向和縱向分別設(shè)置為200個(gè)網(wǎng)格，采用完全匹配層(Perfectly matched Layer，PML)吸收邊界條件在波場外圍布置100個(gè)網(wǎng)格的吸收衰減層避免邊界反射，時(shí)間步長為0.25ms。訓(xùn)練數(shù)據(jù)集A所用正演參數(shù)如表1所示，對每種正演參數(shù)數(shù)值模擬獲得含頻散數(shù)據(jù)和樣本標(biāo)簽數(shù)據(jù)，對每一個(gè)樣本、標(biāo)簽進(jìn)行傅里葉變換，得到波數(shù)域的樣本、標(biāo)簽。將數(shù)據(jù)集A分別按照8∶1∶1的比例劃分為訓(xùn)練集、驗(yàn)證集和測試集，對網(wǎng)絡(luò)模型各部分的有效性進(jìn)行驗(yàn)證。

表1 正演建模參數(shù)

3.1.1 學(xué)習(xí)率更新對網(wǎng)絡(luò)的影響

學(xué)習(xí)率是深度學(xué)習(xí)中重要的超參數(shù)之一，由于L1損失在網(wǎng)絡(luò)更新時(shí)的梯度始終相同，在接近最優(yōu)值處可能維持較大的梯度而錯(cuò)過最優(yōu)值，因此本文在訓(xùn)練迭代過程中自動調(diào)節(jié)學(xué)習(xí)率：當(dāng)驗(yàn)證集的損失函數(shù)在3次迭代后不再下降時(shí)，將學(xué)習(xí)率調(diào)整為當(dāng)前學(xué)習(xí)率的0.5倍。在訓(xùn)練過程中初始學(xué)習(xí)率設(shè)為0.003，批大小為8，迭代次數(shù)為400，學(xué)習(xí)率的更新與SNR的變化關(guān)系如圖6所示，分別在迭代64、92、124、156、200、254、285、316、348、380次時(shí)更新學(xué)習(xí)率，使目標(biāo)函數(shù)重新收斂到全局最小值，得到最優(yōu)解。

圖6 學(xué)習(xí)率更新與網(wǎng)絡(luò)性能的變化

3.1.2 稀疏約束正則項(xiàng)效果分析

將不含稀疏約束項(xiàng)的網(wǎng)絡(luò)模型記為算法Ⅰ，本文算法和算法Ⅰ對驗(yàn)證集的SNR變化曲線如圖7所示，可以看出兩個(gè)網(wǎng)絡(luò)模型在驗(yàn)證集性能上存在較大差別，原因在于算法Ⅰ在訓(xùn)練過程中網(wǎng)絡(luò)過分?jǐn)M合訓(xùn)練數(shù)據(jù)的特征，導(dǎo)致在未知樣本集上表現(xiàn)不佳； 本文算法在驗(yàn)證集上的SNR接近于訓(xùn)練集，提升了網(wǎng)絡(luò)的泛化能力。

圖7 引入稀疏約束前、后的網(wǎng)絡(luò)性能隨迭代次數(shù)的變化

在神經(jīng)網(wǎng)絡(luò)提取特征的學(xué)習(xí)過程中，淺層卷積核提取紋理等細(xì)節(jié)特征，較深層次卷積核提取波場數(shù)據(jù)的輪廓、形狀等抽象語義特征。為了進(jìn)一步說明稀疏約束項(xiàng)的效果，對比兩種算法隱藏層中較深層次卷積核權(quán)重分布如圖8所示。由圖8a可以看出，結(jié)合稀疏約束項(xiàng)的卷積核所學(xué)習(xí)的參數(shù)特征分布基本一致，表明本文算法與算法Ⅰ相比，在特征學(xué)習(xí)時(shí)可較為穩(wěn)定地提取波場數(shù)據(jù)的深層語義信息； 算法Ⅰ中卷積核學(xué)習(xí)的權(quán)重較為分散(圖8b)，由于對卷積核施加L1正則約束，本文算法卷積核權(quán)重分布服從拉普拉斯分布，更多較小權(quán)重被懲罰近似為0，實(shí)現(xiàn)特征優(yōu)化選擇。綜上分析可知加入稀疏約束的網(wǎng)絡(luò)模型具有一定的合理性。

3.1.3 空間—波數(shù)域聯(lián)合誤差效果分析

將只含空間域誤差與稀疏約束項(xiàng)的模型記為算法Ⅱ，只含波數(shù)域誤差與稀疏約束項(xiàng)的模型記為算法Ⅲ。分別使用算法Ⅱ、算法Ⅲ和本文算法對波場數(shù)據(jù)進(jìn)行頻散壓制，均進(jìn)行400次迭代，不同算法網(wǎng)絡(luò)性能隨迭代次數(shù)的變化如圖9所示。從SNR隨迭代次數(shù)變化曲線可知，算法Ⅱ由于未結(jié)合波數(shù)域誤差，在網(wǎng)絡(luò)學(xué)習(xí)過程中未能充分提取特征，結(jié)果地SNR低于本文算法。算法Ⅱ、算法Ⅲ和本文算法進(jìn)行400次迭代后SNR的值分別為22.87dB、6.72dB與34.82dB，本文算法的SNR高于使用單一空間域或波數(shù)域誤差作為損失函數(shù)的算法(圖9)。從算法Ⅲ中第400次迭代中任取一個(gè)樣本展示其壓制頻散效果，由于僅含波數(shù)域計(jì)算目標(biāo)函數(shù)，能量沒有被恢復(fù)，所提取特征有限導(dǎo)致無法根據(jù)現(xiàn)有特征重構(gòu)出完整波形，頻散紋理仍存在于有效波形附近且波場快照分辨率較低(圖10)。

圖8 本文算法(a)與算法Ⅰ(b)部分卷積核權(quán)值統(tǒng)計(jì)左、中、右分別為第15、16、17層

圖9 是否使用聯(lián)合誤差的網(wǎng)絡(luò)性能隨迭代次數(shù)的變化

3.1.4 殘差學(xué)習(xí)效果分析

本文所提出模型在訓(xùn)練過程中對頻散形成的偽波紋理及偽振動所造成的能量損失進(jìn)行學(xué)習(xí)，之后由輸入數(shù)據(jù)與殘差數(shù)據(jù)相減得到壓制頻散后的干凈波場。

將本文方法與不含殘差學(xué)習(xí)策略的模型(記為算法Ⅳ)進(jìn)行對比。由二者SNR隨迭代次數(shù)變化的曲線(圖11a)可見，本文算法比算法Ⅳ更容易優(yōu)化，且網(wǎng)絡(luò)收斂更快(圖11b)。從總體壓制頻散效果來看，引入殘差學(xué)習(xí)策略的模型優(yōu)于直接映射訓(xùn)練的網(wǎng)絡(luò)模型。

圖10 算法Ⅲ頻散壓制結(jié)果

圖11 是否引入殘差學(xué)習(xí)的網(wǎng)絡(luò)性能隨迭代次數(shù)的變化(a)SNR曲線對比； (b)損失函數(shù)的對比

3.1.5 不同損失函數(shù)效果分析

損失函數(shù)的選擇對網(wǎng)絡(luò)性能有很大影響。將使用L2范數(shù)為損失函數(shù)的模型記為算法Ⅴ，將該算法與本文算法進(jìn)行對比測試。訓(xùn)練過程中其余參數(shù)保持不變，兩種算法均迭代400次。由圖12可以看出，算法Ⅴ訓(xùn)練時(shí)的SNR曲線較為平滑，收斂時(shí)的SNR遠(yuǎn)低于本文模型，這是由于算法Ⅴ使用L2損失度量預(yù)測值與真實(shí)值差異時(shí)將誤差平方化，在進(jìn)行頻散壓制的任務(wù)中會賦予偽波紋理較大的權(quán)重，以犧牲有效信號的精度為代價(jià)，朝著減少偽波紋理的方向更新網(wǎng)絡(luò)，降低了模型的整體性能。

圖12 使用不同損失函數(shù)的網(wǎng)絡(luò)性能隨迭代次數(shù)的變化

3.1.6 不同網(wǎng)絡(luò)結(jié)構(gòu)模型測試

將上述訓(xùn)練好的網(wǎng)絡(luò)模型保存，使用與訓(xùn)練數(shù)據(jù)不同震源位置的波場切片對上述網(wǎng)絡(luò)模型進(jìn)行測試。圖13為不同算法測試后的波場快照SNR統(tǒng)計(jì)，可見本文所提出的算法在測試數(shù)據(jù)上具有更高的信噪比，使用結(jié)合稀疏約束的聯(lián)合損失函數(shù)有效降低了傳統(tǒng)損失函數(shù)對波場重構(gòu)帶來的負(fù)面影響，在頻散壓制效果上有明顯優(yōu)勢。

Kaur等[29]使用生成對抗網(wǎng)絡(luò)(GANS)中的cycleGAN網(wǎng)絡(luò)壓制有限差分法中的數(shù)值頻散。將含頻散的波場數(shù)據(jù)記為源域，無頻散偽影的波場數(shù)據(jù)記為目標(biāo)域。該網(wǎng)絡(luò)的目的是將數(shù)據(jù)由源域的特征轉(zhuǎn)換為目標(biāo)域特征，采用了循環(huán)一致性損失等多種損失函數(shù)來增強(qiáng)由源域至目標(biāo)域的映射約束。為了比較兩種網(wǎng)絡(luò)對數(shù)值頻散的壓制效果，選用與本文均勻速度模型相同的數(shù)據(jù)集進(jìn)行對比實(shí)驗(yàn)。cy-cleGAN網(wǎng)絡(luò)參數(shù)為：將cycleGAN網(wǎng)絡(luò)生成器部分調(diào)整為9個(gè)殘差塊，網(wǎng)絡(luò)初始學(xué)習(xí)率設(shè)為0.0002，每迭代50次調(diào)整學(xué)習(xí)率，批大小數(shù)量設(shè)為8。兩個(gè)網(wǎng)絡(luò)模型均迭代400次。為了進(jìn)行兩個(gè)模型泛化性能的對比，將訓(xùn)練好的網(wǎng)絡(luò)模型保存，使用與訓(xùn)練數(shù)據(jù)震源位置不同的波場切片對兩種網(wǎng)絡(luò)模型進(jìn)行測試，結(jié)果如圖14所示。

圖13 不同算法壓制頻散后的SNR對比

由兩種網(wǎng)絡(luò)模型測試結(jié)果(圖14)可以看出，頻散現(xiàn)象均被壓制，但cycleGAN網(wǎng)絡(luò)壓制結(jié)果的波形邊緣不清晰，影響了波場快照質(zhì)量。由于cycleGAN網(wǎng)絡(luò)在壓制過程中僅僅關(guān)注于波場數(shù)據(jù)的空間域特征，無法獲取更多維度特征去加強(qiáng)對目標(biāo)域的映射，轉(zhuǎn)化后的波場快照質(zhì)量與樣本標(biāo)簽存在一定差距； 從運(yùn)行時(shí)間上看，由于GAN網(wǎng)絡(luò)復(fù)雜度較高，因此訓(xùn)練難度較大，迭代400次所需時(shí)間為137h42min，而本文算法訓(xùn)練耗時(shí)為6h3min，本文網(wǎng)絡(luò)模型的運(yùn)行時(shí)間遠(yuǎn)小于cycleGAN。

3.1.7 模型的適用性分析

為了驗(yàn)證本文模型的適用性能，選用不同正演參數(shù)的波場數(shù)據(jù)進(jìn)行測試(表2)。SNR統(tǒng)計(jì)結(jié)果表明本文網(wǎng)絡(luò)對未經(jīng)訓(xùn)練過的數(shù)據(jù)頻散壓制效果明顯。測試數(shù)據(jù)集d和測試數(shù)據(jù)集e的處理結(jié)果如圖15所示，可見，數(shù)值頻散形成的偽波紋理被充分壓制，得到較為清晰的波形邊緣，說明本文網(wǎng)絡(luò)具有一定的泛化能力和適用性，不僅僅在訓(xùn)練數(shù)據(jù)上有著很好的擬合效果，當(dāng)模型不變，正演數(shù)值模擬參數(shù)在一定范圍內(nèi)發(fā)生變化時(shí)，無需重新訓(xùn)練即可獲取高質(zhì)量的波場數(shù)據(jù)。

圖14 cycleGAN網(wǎng)絡(luò)與本文網(wǎng)絡(luò)壓制頻散效果對比(a)與訓(xùn)練數(shù)據(jù)震源位置不同的測試數(shù)據(jù)； (b)cycleGAN網(wǎng)絡(luò)的輸出結(jié)果； (c)本文模型的輸出結(jié)果

表2 測試數(shù)據(jù)集正演參數(shù)及頻散壓制前、后信噪比統(tǒng)計(jì)

3.2 復(fù)雜Marmousi模型

選取Marmousi速度模型(圖16)驗(yàn)證本文方法在復(fù)雜波場情況下的頻散壓制能力。模型網(wǎng)格數(shù)為298×300，PML吸收邊界厚度設(shè)為100個(gè)網(wǎng)格點(diǎn)。網(wǎng)格尺寸為8m×8m，時(shí)間步長為0.25ms。由于Marmousi模型中包含豐富的反射界面與速度信息，為使訓(xùn)練數(shù)據(jù)包含波場信息更豐富，分別在(800m，800m)和(1200m，1200m)處設(shè)置主頻為40Hz的Ricker子波作為震源，兩處震源的波場快照按1∶1的比例制作成訓(xùn)練數(shù)據(jù)集B，并按照8∶1∶1的比例劃分為訓(xùn)練集、驗(yàn)證集與測試集。

使用本文算法對訓(xùn)練數(shù)據(jù)集B進(jìn)行壓制頻散處理，400次迭代后一個(gè)樣本的壓制頻散結(jié)果如圖17所示，可以看出：訓(xùn)練輸入樣本出現(xiàn)很多同相軸，頻散形成的偽波與反射波混疊，波形邊緣模糊且分辨率很低，數(shù)值頻散現(xiàn)象嚴(yán)重； 經(jīng)過頻散壓制后波場快照的分辨率有所提升，去除了大量的偽波紋理的同時(shí)保留了有效信號，具有較高的保真性。

圖15 本文網(wǎng)絡(luò)模型對不同參數(shù)正演波場的頻散壓制結(jié)果(a)測試數(shù)據(jù)集d； (b)測試數(shù)據(jù)集e。左：原始數(shù)據(jù)； 右：網(wǎng)絡(luò)輸出結(jié)果

圖16 Marmousi速度模型

抽取圖17在x=720m處波形曲線(圖18a)，并計(jì)算波數(shù)譜(圖18b)。由圖18a可以看出，頻散波波形經(jīng)過網(wǎng)絡(luò)頻散壓制后，輸出的波形與樣本標(biāo)簽接近一致，有效去除了頻散干擾。由圖18b可以看出，波場中低波數(shù)部分的能量對應(yīng)著真正的波場值，數(shù)值頻散形成的能量集中在高波數(shù)區(qū)域，輸入數(shù)據(jù)頻散現(xiàn)象較為嚴(yán)重，在沒有振動的高波數(shù)成分上也產(chǎn)生振動，形成了虛假繞射波； 頻散壓制后數(shù)據(jù)的高波數(shù)分量上的幅值低于輸入數(shù)據(jù)，意味著頻散能量被有效壓制。

圖17 Marmousi模型訓(xùn)練樣本頻散壓制結(jié)果(a)訓(xùn)練樣本； (b)樣本標(biāo)簽； (c)第400次迭代后的網(wǎng)絡(luò)輸出，SNR=24.19dB

圖18 Marmousi模型模擬數(shù)據(jù)頻散壓制前、后單道數(shù)據(jù)(a)及其波數(shù)譜(b)對比

將訓(xùn)練好的網(wǎng)絡(luò)模型保存，保持其他參數(shù)不變，使用與訓(xùn)練數(shù)據(jù)不同時(shí)刻的波場切片對模型進(jìn)行測試，為對比測試前后的SNR，同樣利用正演記錄獲得測試數(shù)據(jù)的樣本標(biāo)簽，頻散壓制前、后的空間、波數(shù)域如圖19所示。

由圖19a可以看出，頻散現(xiàn)象較為嚴(yán)重，信噪比較低； 經(jīng)過本文模型頻散壓制后(圖19b)，波場切片變得清晰，對于差分精度低造成的能量損失也得到了較精確的補(bǔ)償，消除了頻散帶來的干擾，在空間域內(nèi)頻散所形成的偽波被去除、在波數(shù)域內(nèi)高波數(shù)分量的系數(shù)被壓制，能量更為集中。

為進(jìn)一步對比本文所提方法的壓制頻散效果，將網(wǎng)絡(luò)模型與本文相同、使用MSE為損失函數(shù)且不含波數(shù)域誤差和稀疏約束項(xiàng)的模型記為算法Ⅵ。圖20a為算法Ⅵ的壓制頻散結(jié)果，圖20b為本文算法的輸出結(jié)果，圖20c和20d分別為兩種算法輸出結(jié)果與標(biāo)簽數(shù)據(jù)的殘差。由于算法Ⅵ僅提取了波場的空間域特征進(jìn)行學(xué)習(xí)，且未引入稀疏約束項(xiàng)，因此頻散偽波紋理學(xué)習(xí)不夠充分，少量頻散偽波仍存在于壓制后的數(shù)據(jù)中，波形整體邊緣模糊，頻散壓制效果較差； 從殘差剖面上看，本文壓制頻散結(jié)果與標(biāo)簽數(shù)據(jù)的殘差更小，對頻散造成的波形和能量損失學(xué)習(xí)得較為充分。抽取兩種方法頻散壓制結(jié)果在z=600m處進(jìn)行單道分析(圖21)。由圖21a可以看出，測試樣本的振動曲線波動劇烈，算法Ⅵ壓制結(jié)果仍存在鋸齒狀的波動，本文算法頻散壓制后可得到完整的Ricker子波波形。由圖21b可以看出，本文算法頻散壓制后，與樣本標(biāo)簽相比振幅誤差約為0.04，遠(yuǎn)小于算法Ⅵ的預(yù)測誤差，說明本文算法對振幅補(bǔ)償效果較好。由圖21c可以看出，本文算法壓制后的波數(shù)譜曲線高波數(shù)能量低于算法Ⅵ壓制結(jié)果，且充分保留了低波數(shù)分量能量，說明本文方法有效壓制了數(shù)值頻散且沒有破壞有效信號，驗(yàn)證了本文算法處理復(fù)雜波場數(shù)據(jù)的有效性。

圖19 Marmousi模型不同時(shí)刻測試樣本(a)及其頻散壓制結(jié)果(b)左：空間域； 右：波數(shù)域

圖20 本文算法與算法Ⅵ壓制頻散效果分析(a)算法Ⅵ壓制頻散結(jié)果； (b)本文算法壓制頻散結(jié)果； (c)圖a與標(biāo)簽數(shù)據(jù)的差； (d)圖b與標(biāo)簽數(shù)據(jù)的差圖a、圖b中大方框?yàn)樾》娇騼?nèi)的放大顯示

3.3 斷陷模型

為了測試本文方法對不同地下介質(zhì)速度模型壓制頻散的能力，設(shè)計(jì)一個(gè)斷陷速度模型(圖22)。該速度模型橫向網(wǎng)格數(shù)為400，縱向網(wǎng)格數(shù)為300，網(wǎng)格尺寸為10m×10m。在(2000m，1500m)處激發(fā)40Hz的Ricker子波作為震源，采用時(shí)間2階、空間2階和12階的有限差分?jǐn)?shù)值算法完成該斷層速度模型地震波場正演數(shù)值模擬，得到新的數(shù)據(jù)集C。

結(jié)合遷移學(xué)習(xí)的思想應(yīng)用該斷陷模型數(shù)據(jù)集對網(wǎng)絡(luò)進(jìn)行再訓(xùn)練，具體實(shí)現(xiàn)方式如下：首先，利用多種地下介質(zhì)速度模型所得正演波場數(shù)據(jù)作為源域進(jìn)行訓(xùn)練，得到預(yù)訓(xùn)練網(wǎng)絡(luò)模型； 選取數(shù)據(jù)集C中的20%作為目標(biāo)域的訓(xùn)練集，將預(yù)訓(xùn)練模型隱層參數(shù)賦值到新數(shù)據(jù)集模型相應(yīng)層中進(jìn)行權(quán)重初始化，訓(xùn)練迭代300次； 使用目標(biāo)域的80%剩余數(shù)據(jù)進(jìn)行測試，驗(yàn)證經(jīng)過遷移學(xué)習(xí)訓(xùn)練的網(wǎng)絡(luò)模型效果。

圖23為測試數(shù)據(jù)中一個(gè)樣本及經(jīng)網(wǎng)絡(luò)頻散壓制后的結(jié)果，測試樣本的SNR為0.9357dB，經(jīng)頻散壓制后SNR達(dá)到26.5088dB，可見經(jīng)頻散壓制后，有效波形突顯出來，可清晰觀察到反射波波形。

將在其他正演數(shù)據(jù)集上預(yù)訓(xùn)練過的網(wǎng)絡(luò)模型作為初始網(wǎng)絡(luò)，從而避免了模型隨機(jī)初始化參數(shù)時(shí)的逐步調(diào)優(yōu)過程，使用少量的新的模型數(shù)據(jù)對網(wǎng)絡(luò)進(jìn)行訓(xùn)練，就可得到高質(zhì)量波場數(shù)據(jù)，可見通過遷移學(xué)習(xí)的方式節(jié)約了訓(xùn)練成本、提升了頻散壓制的效率。

圖21 本文算法與算法Ⅵ的散壓制效果單道分析(a)波形對比； (b)殘差對比； (c)波數(shù)域?qū)Ρ?/p>

圖22 斷陷模型

圖23 斷陷模型測試樣本(a)及頻散壓制結(jié)果(b)

4 結(jié)論

本文提出了一種基于聯(lián)合學(xué)習(xí)與稀疏約束的波場數(shù)值頻散壓制方法，該方法通過對波場數(shù)據(jù)的空間域和波數(shù)域特征分析定義了聯(lián)合損失函數(shù)，并在網(wǎng)絡(luò)中增加稀疏約束，達(dá)到頻散壓制目的。通過均勻模型、復(fù)雜的Marmousi模型以及斷陷模型驗(yàn)證了本文所提方法的有效性、適用性與魯棒性。獲得的結(jié)論如下：

(1)與傳統(tǒng)數(shù)值頻散壓制方法相比，可以較低的計(jì)算代價(jià)完成頻散壓制，避免了FCT方法所帶來的振幅損失或計(jì)算量過大的問題；

(2)聯(lián)合空間域和波數(shù)域的信息定義損失函數(shù)，能夠更充分抑制頻散形成的高波數(shù)分量，使頻散特征提取更充分、對頻散壓制更徹底；

(3)選取合理的稀疏正則化約束參數(shù)對卷積核權(quán)重進(jìn)行約束，使網(wǎng)絡(luò)頻散壓制效果提升且對于未經(jīng)訓(xùn)練的未知樣本頻散壓制效果較好，提升網(wǎng)絡(luò)模型的泛化性能；

(4)由于邊緣保持能力及魯棒性，L1范數(shù)損失函數(shù)更適用于頻散壓制，結(jié)果更清晰，邊緣更明顯；

(5)本文算法與遷移學(xué)習(xí)結(jié)合，對新模型的正演波場進(jìn)行頻散壓制可取得較好的實(shí)驗(yàn)結(jié)果，在節(jié)約訓(xùn)練成本的同時(shí)使頻散壓制效率有所提升。

三個(gè)不同速度模型的實(shí)驗(yàn)結(jié)果表明，本文方法在有效壓制頻散的同時(shí)對有效波形損傷較小，有效保護(hù)了波形的完整性和邊緣的光滑性。