張 巖 崔淋淇 宋利偉 董宏麗

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

1 數(shù)值頻散分析

1.1 數(shù)值頻散

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

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

(1)

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

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

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

(2)

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

(3)

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

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

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

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

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

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

(4)

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

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

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

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

2.1 網(wǎng)絡模型

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

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

(5)

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

(6)

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

f(l)=max(0,l)

(7)

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

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

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

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

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

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

(8)

(9)

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

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

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

2.2.2 稀疏約束正則項

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

(10)

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

(11)

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

(12)

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

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

3 模型實驗及結(jié)果

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

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

(13)

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

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

3.1 均勻速度模型實驗

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

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

3.1.1 學習率更新對網(wǎng)絡的影響

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

圖6 學習率更新與網(wǎng)絡性能的變化

3.1.2 稀疏約束正則項效果分析

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

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

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

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

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

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

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

3.1.4 殘差學習效果分析

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

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

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

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

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

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

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

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

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

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

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

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

3.1.7 模型的適用性分析

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

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

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

3.2 復雜Marmousi模型

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

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

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

圖16 Marmousi速度模型

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

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

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

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

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

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

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

圖20 本文算法與算法Ⅵ壓制頻散效果分析(a)算法Ⅵ壓制頻散結(jié)果； (b)本文算法壓制頻散結(jié)果； (c)圖a與標簽數(shù)據(jù)的差； (d)圖b與標簽數(shù)據(jù)的差圖a、圖b中大方框為小方框內(nèi)的放大顯示

3.3 斷陷模型

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

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

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

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

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

圖22 斷陷模型

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

4 結(jié)論

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

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

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

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

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

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

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