齊紅宇 傅紅筍 楊 露

(大連海事大學理學院， 遼寧大連 116026)

0 引言

全波形反演(Full Waveform Inversion，F(xiàn)WI)充分利用地震記錄中包含的豐富信息，可獲得更精細的深度域速度模型，成為油氣勘探領域的研究熱點[1]。20世紀80年代，Lailly[2]和Tarantola[3]率先在時域引入基于廣義最小二乘優(yōu)化原理的FWI理論與方法，其開創(chuàng)性的工作對反演理論的發(fā)展產(chǎn)生了重要而深遠的影響。隨后，Pratt[4]將FWI推廣到頻率域，并提出由低頻到高頻的反演格式，只需要幾個離散頻率就能獲得較高精度的反演結果。然而，由于觀測數(shù)據(jù)中低頻信息的缺失和實際觀測中大量噪聲的存在，F(xiàn)WI通常是嚴重不適定問題，主要表現(xiàn)為多解性和對噪聲的敏感性[5-7]。為此，必須采用有效的正則化技術求取其穩(wěn)定的近似解。

經(jīng)典的Tikhonov正則化[8]選取L2范數(shù)作為罰項，其作用在于通過壓制大的分量而產(chǎn)生光滑的近似解。全變差(Total Variation，TV)正則化[9]能在抑制噪聲的同時保持尖銳邊緣。近十幾年來，該方法已廣泛應用于求解地球物理反問題，實現(xiàn)了對層間尖銳界面的有效識別，以及高對比度鹽體的有效成像[10-11]。但TV正則化在反演過程中會導致圖像結構細節(jié)信息的丟失，并產(chǎn)生“階梯偽影”[12-13]。

由于在復雜地質(zhì)背景下，速度模型中可能同時包含平滑結構和銳利邊界，單一的正則化技術很難提供較好的反演結果。為此，Lin等[14]建立了基于修正TV正則化的FWI模型，將FWI分解為兩個交錯的子問題：一是具有Tikhonov正則化項的傳統(tǒng)FWI問題，二是基于一階TV正則化去噪問題，得到較理想的反演結果，但階梯狀偽影依然存在。Peng等[15]將TV正則化與箱約束相結合，提升模型的物理可解釋性，并提出改善反演精度的自適應原始對偶梯度法。Aghamiry等[16]將TV正則化項與Tikhonov正則化項加權耦合，并利用交替方向乘子法求解相應的目標泛函，數(shù)值試算結果表明所提耦合方法優(yōu)于單一正則化方法。

本文將速度模型擾動量的L2范數(shù)與TV正則化相融合，構建能克服FWI的不適定性的優(yōu)化目標泛函，針對目標泛函的不可微性，提出修正正交有限內(nèi)存擬牛頓算法(modified Orthant-Wise Limited- Memory Quasi-Newton，mOWL-QN)進行迭代計算；該算法基于有限內(nèi)存擬牛頓方法[17]，在下降的過程中引入正交函數(shù)尋找下降方向以避免不可微項的求導。將mOWL-QN算法引入FWI，不僅提高了FWI的計算效率，而且可引導反演結果更貼近真實速度?；诰哂袕碗s構造的修正Marmousi模型和BG Compass模型的數(shù)值試驗，且與無正則化FWI及鄰近有限內(nèi)存擬牛頓(proximal Limited-memory Broyden-Fletcher-Goldfarb-Shanno，proxL-BFGS)方法[18]進行比較，驗證了mOWL-QN算法求解混合正則化FWI的有效性。

1 原理方法

1.1 頻域FWI的混合正則化模型

從觀測數(shù)據(jù)與計算數(shù)據(jù)相吻合的角度出發(fā)，F(xiàn)WI可歸結為如下的非線性二次優(yōu)化問題

(1)

為了克服FWI的不適定性，獲得一個既包含光滑內(nèi)部結構，又具有銳利邊緣的速度模型，將Tikho-nov與TV正則化項線性加權后添入式(1)，得到

(2)

求解式(2)的主要計算困難在于L1范數(shù)的不可微性。為了解決這一問題，Gong等[19]提出mOWL-QN算法。mOWL-QN方法先利用TV正則項內(nèi)變量的每個坐標都不變號的特點計算正則項的導數(shù)；再使用當前正交對象所求得的導數(shù)構造一個搜索方向，并將其限制在該正交對象上[20]。mOWL-QN算法在高效計算不可微優(yōu)化問題的同時繼承了LBFGS方法求解大規(guī)模問題的優(yōu)點，所以本文利用mOWL- QN算法求解式(2)，以達到權衡計算效率與反演精度的目的。

1.2 mOWL-QN算法

首先，給出mOWL-QN算法中必需的正交函數(shù)π(x，y)的定義[21]

函數(shù)π(x，y)主要目的是限制x與y的象限相同； sign(·)為符號函數(shù)。

為了簡化符號，令

將式(2)簡寫成

(3)

式中：g(m)為f(m)的梯度，即gi(m)=?f/?mi，mi表示參數(shù)m的第i個元素；wi表示δm向量化后的第i個元素。

對于大規(guī)模非線性優(yōu)化問題，LBFGS方法常被用來計算近似的Hessian陣的逆。其主要思想是通過存儲向量的方式避免存儲大量的矩陣，提升計算效率。為此，式(3)的下降方向為

dk=-Hk◇g(mk)

(4)

式中：Hk是由LBFGS方法計算的目標泛函Φ(m)的Hessian陣的逆； ◇g(mk)是目標函數(shù)的偽梯度；k為迭代次數(shù)。同時，將正交函數(shù)作用在下降方向dk上，得到新的下降方向

pk=π[dk，-◇g(mk)]

并利用正交函數(shù)π(x，y)對m進行更新，得到

mk+1=π(mk+α1pk，ξk)

(5)

但由于式(3)中的不可微項的存在，導致證明式(5)的收斂性存在一定困難。因此，Gong等[19]提出mOWL-QN算法，主要的修正體現(xiàn)在原算法中加入最速下降步。最速下降方向是將梯度的負方向作為模型下一步的迭代方向，并且利用線搜索方法尋找合適的步長。在最速下降法中，對mk+1進行迭代更新的表達式為

mk+1=mk-α2◇g(mk)

(6)

在迭代中，mk+1的更新策略是由集合Ik決定的。集合Ik定義如下

式中ε1是足夠小的正數(shù)。如果集合Ik是空集，那么mk+1=π(mk+αpk，ξk)； 否則，利用最速下降法式(6)進行更新。

表1 mOWL-QN算法的偽代碼

1.3 鄰近有限內(nèi)存擬牛頓算法

使用傳統(tǒng)LBFGS算法求解式(3)具有挑戰(zhàn)性，其主要原因在于TV正則化的不可微性。因此，本文選擇proxL-BFGS作為對比算法求解式(3)。該算法將傳統(tǒng)LBFGS與鄰近算子結合，克服目標泛函的不可微性，得到反演模型參數(shù)的近似解。

(7)

根據(jù)凸優(yōu)化的一階條件可得式(7)的最優(yōu)解

proxG(z)=[I+?G(m)]-1(z)

(8)

式中： ?G(m)為凸泛函G(m)的次梯度；I為單位矩陣。

在LBFGS算法中，迭代更新公式為

(9)

下面給出利用鄰近有限內(nèi)存擬牛頓算法(表2)的偽代碼。

表2 proxL-BFGS算法的偽代碼

為了保證對比實驗的公平性，上述算法步驟中的參數(shù)與本文所提算法mOWL-QN的參數(shù)設置保持一致。

2 數(shù)值算例

數(shù)值試算基于具有復雜構造的修正Marmousi模型及BG Compass模型進行。考慮到算法的實用性，實驗數(shù)據(jù)中添加了1%的高斯噪聲。本文選擇結構相似性指數(shù)(SSIM)、峰值信噪比(PSNR)及均方根誤差(RMSE)三項指標對算法進行評估。其中：SSIM和PSNR越高表明反演結果越接近真實值； RMSE越低表明反演結果與真實值之間的誤差越小。為了保證實驗的公平性，所有實驗均使用Matlab 2018a進行反演試算，并且在內(nèi)存為64GB、八核處理器的計算機環(huán)境中運行。

2.1 修正Marmousi模型

首先，對修正Marmousi模型進行試算，其真實模型如圖1b所示，尺寸為Nz×Nx=260×891。本文選擇平滑處理后的初始模型(圖1a)，并將網(wǎng)格間距設置為3.8m。在模型的表層放置77個點源及231個接收器(接收點源發(fā)射的完整數(shù)據(jù))。為了避免周期跳躍性，在[2Hz，26Hz]頻率范圍選取12個重復頻段進行反演。實驗中，最大迭代次數(shù)為60，正則化參數(shù)分別為λ=10-3、β=10-6，設置誤差限ε1=ε2=10-7，選擇α0=1作為初始步長。

圖1 修正Marmousi模型

三種算法得到的反演結果如圖2所示。仔細觀察對比這三種反演結果，無正則化FWI的反演結果(圖2a)存在較多偽影且無法恢復模型的精細結構；本文算法(圖2c)更好地保留了地質(zhì)模型的細節(jié)；與proxL-BFGS算法(圖2b)相比，本文算法還可抑制噪聲的影響，且只存在較少偽影。圖3展示的是三種算法得到的反演結果在x=1880m處的速度曲線。雖然proxL-BFGS算法與mOWL-QN算法的反演速度模型都很貼合真實速度，但mOWL-QN算法得到的反演速度(紅色實線)與模型真實速度(藍色虛線)擬合效果更好。算法的評估指標結果如表3所示，本文mOWL-QN算法以少量的計算時間增加(比無正則化FWI多耗時約14%)在三種定量評估數(shù)值上表現(xiàn)最佳，且全方面完勝proxL-BFGS算法。

圖2 修正Marmousi模型的三種方法反演結果

圖3 修正Marmousi模型的三種方法反演的縱向速度曲線對比

表3 修正Marmousi模型三種方法反演的指標及耗時統(tǒng)計

2.2 BG Compass模型

再對BG Compass模型進行試算。設定網(wǎng)格間隔為10m，真實模型網(wǎng)格點數(shù)為205×701(圖4b)； 無任何橫向信息的初始模型如圖4a。在靠近地表的水平面上放置64個點源和192個接收器。同樣選取12個重復頻段進行反演，頻率范圍是2.6～28Hz。在實驗中，設置正則化參數(shù)λ=10-2、β=10-5，其他參數(shù)同前述實驗。

圖4 BG Compass 模型

最終反演結果如圖5所示，三種算法都可得到較滿意的反演效果。mOWL-QN算法(圖5c)能更多地保留地質(zhì)結構細節(jié)，具有更好的視覺效果。圖6為三種算法所得反演結果在x=1500m處的速度曲線?？梢钥闯?，相對于無正則化FWI及proxL- BFGS算法，本文算法的擬合程度更高，尤其在模型的下半部分，mOWL-QN算法得到的反演速度(紅色實線)更接近模型的真實速度(藍色虛線)。表4展示的評估指標統(tǒng)計結果表明本文算法在定量分析上具有優(yōu)越性。該實驗結論與前述算例保持一致，再次驗證了mOWL-QN算法的有效性。

圖5 BG Compass模型的三種方法反演結果

圖6 BG Compass模型三種方法反演的縱向速度曲線對比

表4 BG Compass模型三種方法反演的評價指標及耗時統(tǒng)計

3 結束語

FWI是一種非線性的、不適定的優(yōu)化問題，通常使用正則化技術緩解其不適定性。應用單一的Tikhonov正則化會產(chǎn)生過度光滑的解，將Tikhonov正則化與TV正則化相融合，能有效改進反演成像的效果。為了克服混合正則化目標泛函的不可微性，本文提出mOWL-QN反演優(yōu)化迭代算法，該算法繼承了擬牛頓方法優(yōu)點，且在一定程度上提升了大規(guī)模反演問題的計算效率，得到高分辨率成像結果。基于修正Marmousi模型和BG Compass模型進行的數(shù)值仿真和對比試驗，驗證了所提方法的有效性。