全波形反演正則化方法對比

李昕潔 王維紅 郭雪豹 張庭俊

(東北石油大學地球科學學院，黑龍江大慶 163318)

0 引言

精確的速度模型是復雜介質疊前深度偏移準確成像的前提。傳統的旅行時層析和偏移速度分析等方法已無法滿足當前高精度成像的需求。相比傳統方法，全波形反演能夠充分利用地震數據中的信息，理論上可以實現更高精度的速度建模[1-4]。早在二十世紀八十年代，Tarantola[5-6]就給出了時間域全波形反演的理論框架，但受制于當時計算機水平，該理論未得到充分發展。隨后，Pratt[7]實現了頻率域全波形反演，并提出多尺度反演的概念。二維情況下，頻率域方法計算效率更高，單一頻率波場模擬僅需一次矩陣分解，即可通過替換震源項實現多炮并行。同時，在多尺度反演框架下，只需由低到高幾個頻率成分即可實現反演過程，不僅降低了計算成本，而且降低了反演的非線性。由此，全波形反演開始得以迅速發展，并逐步應用于實際數據。

全波形反演是高度非線性和不適定的反問題。當初始模型遠離真實模型時，會存在周波跳躍現象，極易陷入局部極值。雖然全局優化方法在理論上可以得到全局最優解，但高昂的計算成本使其難以應用于實際。在目標函數中添加正則化項約束反演過程，是緩解全波形反演非線性和不適定性的一種處理方式，如在模型域中通過將先驗信息應用于模型空間等[8-10]。典型的Tikhonov正則化方法是基于L2范數[11]，能緩解反問題的不適定性，但通常會使反演結果變得光滑，不利于刻畫地層邊界。相較于L2范數正則化，L1范數能使反演結果具有稀疏表征的效果，更符合實際地質情況。

全變分(Total variation，TV)正則化方法[12]通過求解全變分約束的非線性最小化問題實現噪聲壓制，同時保留圖像原本的結構特征。Askan[13]將TV正則化應用于各向異性介質的全波形反演，以同時反演速度和黏彈性參數，在壓制反演結果高振蕩現象的同時，保留了結果的不連續變化。Burs-tedde等[14]以一維全波形反演為例，證明了TV正則化的反演結果比Tikhonov正則化方法更符合實際層狀介質。Anagaw[15]將TV正則化應用于地震成像，詳細給出了TV正則化的構造公式。隨后，TV正則化得到了進一步發展，出現了TV正則化的初對偶實現方法、L2范數+加權L2范數正則化的兩步法[16]等。Anagaw等[17]給出基于TV正則化的全波形反演方法，提高了反演結果分辨地層邊界的能力。

近年來，人們開始將Tikhonov正則化和TV正則化相結合，以期同時具備兩者的優勢。 Lin等[18]提出了改進的TV正則化方法，將其解耦為基于L2范數的Tikhonov正則化和標準L2范數問題；Sun等[19]融合Tikhonov正則化和TV正則化，提出了雙參數整形正則化全波形反演方案，提高了邊界精度和收斂速度[20]；Aghamiry等[21]構建兩類基于Tikhonov和TV正則化的混合正則項：第一類為二者的凸組合(CC)，第二類則是二者的卷積下確界(IC)，認為基于IC的復合正則化反演的誤差最小。

Hu等[22]基于L2范數和加權L2范數正則化方法，提出了一種結合高斯—牛頓方法和乘法正則化的同時多頻反演算法。Guitton等[23]提出了一種L2范數預處理全波形反演，利用傾角估計，通過定向Laplacian濾波器對梯度施加具有地質意義的約束。隨后，Guitton[24]又提出了一種塊狀正則化方法，并比較了L1范數和Cauchy函數正則化項的差異。Ma等[25]研究了圖像引導梯度，在稀疏模型空間中計算梯度，用一種改進的圖像引導共軛梯度方法求解稀疏全波形反演，通過約束梯度與地下傾角補充低頻成分。Asnaashari等[26]提出了一種全波形反演方案，其中有兩個模型懲罰項，Tikhonov項和從聲波測井、鉆井數據或野外地質信息中導出的先驗模型項。Van Leeuwen等[27]提出通過擴展搜索空間減少全波形反演中的局部極小值。Yu等[28]建立了基于L2范數和非光滑L1范數組合的混合雙參數正則化，該方法可以快速收斂到真實模型，反演精度較高。Yan等[29]提出了一種基于Lp-Lq范數的最小化模型，突出了模型結構的邊緣信息。Xue等[30]通過使用Seislet正則化實現全波形反演，Seislet正則化可以在多震源同時激發導致的強串擾偽像和數據存在強隨機噪聲的情況下，恢復高精度速度模型，極大地提高全波形反演的穩健性。利用條件Lipschitz穩定性估計，Shi等[31]設計了迭代正則化的全波形反演，通過將縮放后梯度投影到與參數化相關的子空間實現，保證了Lipschitz穩定性。

在目標函數中添加正則化項是約束全波形反演結果的一種重要途徑，而對比分析不同正則化項對反演結果的影響有助于理解不同正則化的特點。因此，本文首先分別簡要介紹了Tikhonov、TV、雙參數整形、混合雙參數、稀疏結構約束正則化的基本原理，而后通過對背斜—超覆模型、Marmousi模型測試，分析五種正則化方法的特點，為正則化方法的應用提供參考。

1 不同正則化方法的基本原理

基于L2范數的全波形反演目標函數[29]為

(1)

式中：m是速度模型；dcal為模擬數據；dobs為觀測數據；δd為模擬數據與觀測數據的殘差；r(m)為正則化項。

本文在頻域中實現全波形反演，采用伴隨狀態法求取梯度，應用L-BFGS優化算法(Limited-memory Broyden-Fletcher-Goldfarb-Shanno Algorithm)更新模型。

1.1 Tikhonov正則化

Tikhonov正則化使用模型的二次函數作為懲罰項，以獲得具有平滑效果的穩定近似解，其形式為[11,32-33]

(2)

式中：m0是初始速度模型；α是正則化項權重。正則化項權重α保證正則化項占目標函數的5%～10%。α增大會使反演結果更平滑，并降低反演精度。

r(m)的梯度可寫為

?r(m)=2α(m-m0)

(3)

1.2 TV正則化

TV正則化是稀疏約束全波形反演方法的一種，該方法能夠突出并保留反演結果的邊界信息，其形式[34-35]為

(4)

式中：ε1通常被選擇為一個極小的正數；Ω為模型空間范圍。正則化項權重α可以通過實驗獲得，以保證正則化項占目標函數的5%～10%。α越大，邊界越清晰，但會降低反演精度。

r(m)的梯度為

(5)

其離散形式為

(6)

式中：i=1，2，…，Nx；j=1，2，…，Nz。其中Nx、Nz分別為模型在水平和垂直方向的坐標點總個數。

1.3 雙參數整形正則化

基于Tikhonov和TV正則化的雙參數整形正則化項為[19]

(7)

r(m)的梯度為

?r(m)=α[γ1(m-m0)+

(8)

式中：γ1為0到1之間的常數。當γ1趨近于0時，反演結果趨向于TV正則化；當γ1等于1時，接近Tikhonov正則化結果，本文取γ1=0.4。α增大會降低反演精度。

1.4 混合雙參數正則化

混合雙參數正則化反演方法結合了L2范數和L1范數的優勢，使用L2范數擬合實際數據，并通過L1范數約束減少了非唯一性和離散值，形式為[28]

(9)

式中：r1(m)是L2范數正則化項；r2(m)是L1范數正則化項；參數γ2調控反演精度。

r1(m)的梯度與式(3)相同。r2(m)中的L1范數用Huber范數fHuber近似，即

(10)

式中hε2(·)定義為非線性分段函數

(11)

其中ε2通常選擇為一個極小的正數。

r2(m)的梯度可用Huber函數的導數近似

(12)

h′ε2(·)可以寫為

(13)

1.5 稀疏結構約束正則化

地下模型結構具有隱含的稀疏特征，因此可以用L1范數來描述。通過差分算子提取結構信息，利用不同的p和q參數進行Lp-Lq范數約束。

定義由三部分組成的結構化正則化項[29]為

=r1(m)+r2(m)

(14)

式中：λ1、λ2和γ3是正則化項權重；r1(m)為差分算子正則化項；r2(m)為L1范數正則化項。

對于p=q=2，r1(m)的梯度可寫為

(15)

當p=q=1時，L1范數用Huber范數近似，L1范數正則化項的梯度可用Huber函數的導數近似。

在不同的反演階段，權重的選擇不相同。r1和r2的最大值不應超過E(m)的百分之十。本文不同方向的正則化項權重取相同值，即λ1=λ2，也可根據實際需求或地下構造選取不同的參數組合。λ1和λ2在增強反演結果的連續性、去噪方面效果明顯。

2 不同正則化反演結果分析

本文采用背斜—超覆模型和Marmousi模型對比不同正則化方法對反演結果的影響。

2.1 背斜—超覆模型

背斜—超覆模型(圖1a)橫向有386個網格點，縱向有200個網格點，空間采樣間隔為10m。共設置42炮，均勻分布于地表水平方向0.4～3.6km范圍內。每炮126道接收，檢波點均勻分布于地表水平方向0.4～3.6km范圍內，檢波器固定。四周均采用完全匹配層吸收邊界條件。采用真實模型平滑后的結果(圖1b)作為初始模型。

圖1 背斜—超覆速度(v)模型(a)真實；(b)初始

圖2a為背斜—超覆模型未加正則化的反演結果，速度分界面不清晰，底層高速體(最底層紅色的高速體)存在層狀假象，深層(深度1～2km范圍)均衡性差。圖2b為背斜—超覆模型未加正則化的反演結果與真實模型的殘差(Δv)，可見深層反演精度較低，速度分界面殘差較大。

圖2 背斜—超覆速度模型未加正則化的反演結果(a)及其殘差(b)

Tikhonov正則化反演結果(圖3a)有一定的平滑效果，底層高速體層狀假象相對于未正則化反演結果有所減弱，箭頭所指的速度分界面在箭頭右側清晰度明顯降低。TV正則化反演結果(圖3b)底層高速體速度反演精度大幅提高，層狀假象相對于其他反演結果最弱，深層反演均衡性最好，邊界刻畫最清晰，但不能清晰反演出箭頭位置右側的速度分界面。

圖3 背斜—超覆模型不同正則化方法反演結果對比(a)Tikhonov；(b)TV；(c)雙參數整形；(d)混合雙參數；(e)稀疏結構約束(p=q=2)；(f)稀疏結構約束(p=q=1)

圖3c為雙參數整形正則化反演結果，邊界刻畫清晰，層位清晰，底層高速體層狀假象比TV正則化反演結果略強，深層反演均衡性好，可較清晰分辨速度分界面至箭頭位置。相比未正則化反演，混合雙參數正則化反演(圖3d)結果底層高速體層狀假象變弱，深層反演均衡性有所提高。

當p=q=2時，稀疏結構約束正則化反演結果(圖3e)速度連續性最強，層位清晰，底層高速體層狀假象比未正則反演結果弱，深層反演均衡性好。當p=q=1時，稀疏結構約束正則化反演結果(圖3f)底層高速體層狀假象比未正則反演結果弱。

圖4為背斜—超覆模型不同正則化方法反演結果與真實模型殘差。Tikhonov正則化反演結果殘差(圖4a)在矩形框內殘差值小于未正則化，箭頭所指的速度分界面處殘差也小于未正則化反演結果。TV正則化(圖4b)底層高速體與真實模型殘差最?。患^所指的速度分界面殘差值最小(與其他反演結果相比)。雙參數整形正則化(圖4c)在矩形框內與真實模型的殘差小于TV正則化?；旌想p參數正則化(圖4d)在矩形框內和底層高速體殘差比未正則化殘差小。當p=q=2時，稀疏結構約束正則化(圖4e)在淺層矩形框內殘差小于未正則化，整體上相比其他反演結果更接近真實模型速度。當p=q=1時，稀疏結構約束正則化(圖4f)，底層高速體與真實模型殘差小于未正則化。

圖4 背斜—超覆模型不同正則化方法反演結果與真實模型的殘差(a)Tikhonov；(b)TV；(c)雙參數整形；(d)混合雙參數；(e)稀疏結構約束(p=q=2)；(f)稀疏結構約束(p=q=1)

圖5為背斜—超覆模型不同正則化方法反演目標目標函數的下降曲線，可以發現：Tikhonov正則化、TV正則化、雙參數整形正則化目標函數收斂更快；混合雙參數正則化、稀疏結構約束正則化反演目標函數前30次下降速度快，后30次迭代目標函數下降趨勢變緩。

圖5 背斜—超覆模型不同正則化反演目標函數下降曲線

圖6為背斜—超覆模型不同正則化方法反演結果的模擬數據與觀測數據的相對殘差曲線，可以看出，Tikhonov正則化、TV正則化殘差值會略大于未正則化，因為兩側的反演精度降低；雙參數整形、混合雙參數、稀疏結構約束正則化(p=q=2、p=q=1)的反演結果，殘差值均小于未加正則化反演結果的殘差值，說明這三種方法整體上的反演精度高于未加正則化，稀疏結構約束(p=q=2)正則化的反演精度最高。

圖6 背斜—超覆模型不同正則化反演的數據相對殘差橫坐標1～7分別代表：未加正則化、Tikhonov正則化、TV正則化、雙參數整形正則化、混合雙參數正則化、稀疏結構約束正則化(p=q=2)、稀疏結構約束正則化(p=q=1)

圖7為背斜—超覆模型不同正則化方法反演結果在x=1.78km處的單道對比，可見：Tikhonov正則化和混合雙參數正則化速度反演結果略大于未正則反演結果；TV正則化和雙參數整形正則化能更準確地反映分界面的速度突變；稀疏結構約束正則化(p=q=2)反演速度整體上更接近真值。

圖7 背斜—超覆模型反演結果的單道對比

2.2 Marmousi模型

Marmousi模型(圖8a)橫向有630個網格，縱向有200個網格，空間采樣間隔為10m。共設置60炮，均勻分布于地表水平方向0.4～5.8km范圍內。每炮180道接收，均勻分布于地表水平方向0.4～5.8km范圍內，四周均采用完全匹配層吸收邊界條件。采用真實模型的平滑結果作為反演的初始模型(圖8b)。

圖8 Marmousi模型(a)真實速度；(b)初始速度；圖a中三個方框為重點分析區

2.2.1 整體反演結果分析

由于本文采用較為理想的初始模型及多尺度反演方法，因此，在未加正則化的情況下即可較真實反演出與真實速度相匹配的速度場(圖9)，但整體反演結果都有噪聲分布，底部高速體和速度分界面反演精度較低，速度分界面細節刻畫不清晰。Tikhonov正則化反演結果(圖10a)較未正則的反演結果更平滑，速度分界面精度有所提高，箭頭所示的深部高速體形態反映較為真實。TV正則化反演結果(圖10b)層位清晰，界面刻畫、去噪效果很好，在箭頭所示處，速度分界面刻畫更精確，幾乎完全消除了噪聲，但深層(縱向1.6～2km)反演精度有所下降。

圖9 Marmousi模型未加正則化反演結果

圖10c為雙參數整形正則化反演結果，融合了Tikhonov正則化和TV正則化的優勢，不但精確刻畫了界面、壓制了噪聲，而且提高了深層的反演精度。與TV正則化相比，雙參數整形正則化反演的深層高速體形態更真實(箭頭所示)?；旌想p參數正則化反演結果(圖10d)壓制了部分噪聲，速度分界面細節刻畫相對清晰，淺層(縱向0～1km)低速體保真度較高，深層反演精度略有下降(與未加正則反演相比)。

圖10 Marmousi模型不同正則化方法反演結果對比(a)Tikhonov；(b)TV；(c)雙參數整形；(d)混合雙參數；(e)稀疏結構約束(p=q=2)；(f)稀疏結構約束(p=q=1)

稀疏結構約束正則化(p=q=2)反演結果(圖10e)能有效地增強地層的連續性，去噪效果明顯，提高了反演精度；與本文的其他正則化方法相比，圖中箭頭所指的深層高速體，速度連續性最強，速度分界面細節刻畫最清晰。稀疏結構約束正則化(p=q=1)的反演結果(圖10f)淺層低速體保真度較高，與其他正則化方法相比，深層反演精度最高(箭頭所示)。

圖11為Marmousi模型不同正則化方法反演結果的模擬數據與觀測數據的相對殘差曲線，可以發現：Tikhonov、TV、雙參數整形、混合的雙參數、稀疏結構約束(p=q=2)正則化的反演殘差均小于未加正則化，說明這五種方法的反演精度要高于未加正則化；僅稀疏結構約束正則化(p=q=1)反演殘差大于未加正則化；稀疏結構約束正則化(p=q=2)反演殘差最小，反演精度最高。

圖11 Marmousi模型不同正則化反演的數據相對殘差橫坐標1～7分別代表：未加正則化、Tikhonov正則化、TV正則化、雙參數整形正則化、混合雙參數正則化、稀疏結構約束正則化(p=q=2)、稀疏結構約束正則化(p=q=1)

圖12為Marmousi模型不同正則化反演目標函數下降曲線。TV、混合雙參數正則化目標函數趨于穩定的速度更快。所有正則化反演方法，在前20次時目標函數值下降快，后20次迭代目標函數趨于平穩。

圖12 Marmousi模型不同正則化反演目標函數下降曲線

2.2.2 局部反演結果分析

將圖8方框1處Marmousi模型不同正則化反演結果進行局部放大，如圖13和圖14所示。未正則化的反演結果(圖13)的深層高速體地層連續性較差，速度分界面刻畫不夠清晰；Tikhonov正則化反演結果(圖14a)中深層高速體地層連續性有所增強但精度略有降低；TV正則化反演結果(圖14b)的速度分界面刻畫清晰，去噪效果明顯，深層反演精度明顯降低；雙參數整形正則化反演結果(圖14c)在TV正則化反演凸顯速度分界面和去噪的基礎上，提高了深層反演精度；混合雙參數正則化反演結果(圖14d)增強了深層高速體地層的連續性，深層反演精度降低較小。

圖13 Marmousi模型未加正則化反演結果的局部放大(圖8a方框1處)

稀疏結構約束正則化(p=q=2)反演結果(圖14e)，相比于前幾種方法，深層高速體地層連續性最好，保真度最高，速度分界面細節刻畫最清晰，去噪效果明顯，幾乎沒有降低深層反演精度。稀疏結構約束正則化(p=q=1)反演結果(圖14f)提高了深層的反演精度，與未加正則化的反演結果相比速度分界面更清晰。

圖14 Marmousi模型不同正則化反演結果的局部放大(圖8a方框1處)(a)Tikhonov；(b)TV；(c)雙參數整形；(d)混合雙參數；(e)稀疏結構約束(p=q=2)；(f)稀疏結構約束(p=q=1)

將圖8方框2處Marmousi模型不同正則化反演結果的殘差進行局部放大，如圖15和圖16所示，可以直觀地顯示反演結果在不同位置處的失配程度。Tikhonov正則化反演結果(圖16a)比未正則化(圖15)在橫向3.5～5.0km處殘差變小，在橫向5.0～6.0km處殘差略有增大。與其他方法相比在速度分界面處TV正則化反演結果的殘差(圖16b)最小。與TV正則化相比，雙參數整形正則化反演結果的殘差(圖16c)比TV正則化在橫向5.0～5.5km處明顯減小。混合雙參數正則化反演結果的殘差(圖16d)小于未正則化。稀疏結構約束正則化(p=q=2)反演結果的殘差(圖16e)最小，增強了速度的連續性，使反演結果保真度有所提高。稀疏結構約束正則化(p=q=1)反演結果的殘差(圖16f)小于未正則化。

圖15 Marmousi模型未加正則化反演結果殘差的局部放大(圖8a方框2處)

圖16 Marmousi模型不同正則化反演結果殘差的局部放大(圖8a方框2處)(a)Tikhonov；(b)TV；(c)雙參數整形；(d)混合雙參數；(e)稀疏結構約束(p=q=2)；(f)稀疏結構約束(p=q=1)

圖17a是Marmousi模型不同正則化反演結果在x=3km處的單道(圖8a方框3處)對比，可明顯看出，在第一個峰值兩側的速度界面，TV和雙參數整形正則化反演結果更接近真值，其次是稀疏結構約束正則化(p=q=2)的反演結果較接近真實值。圖17b是Marmousi模型不同正則化反演結果殘差在x=3km處的單道(圖8a方框3處)對比，TV正則化和雙參數整形正則化的殘差更小。

圖17 Marmousi模型不同正則化反演結果(a)及其殘差(b)在x=3km處的單道對比

3 結論

本文詳細對比分析了Tikhonov正則化等五種正則化方法在全波形反演中的作用。背斜—超覆模型和Marmousi模型測試表明：

(1)Tikhonov正則化增強反演結果的連續性，使反演結果更平滑；TV正則化在刻畫邊界的同時，去噪效果明顯，但這兩種典型的正則化方法的深層反演精度都有所降低；

(2)雙參數整形正則化方法融合了Tikhonov正則化和TV正則化的優點，并提高了深層的反演精度；

(3)混合雙參數正則化增強了反演結果的地層連續性，對淺層低速體的反演精度有一定程度的提高；

(4)稀疏結構約束正則化(p=q=2)有效增強了反演結果的地層連續性，速度分界面細節刻畫更清晰，去噪效果良好，且深層反演精度降低較??；

(5)稀疏結構約束正則化(p=q=1)，對深層反演精度有所提高。

總之，稀疏結構約束正則化(p=q=2)的反演結果兼顧了各方法在去噪、增強地層連續性以及保留邊緣結構等方面的優勢。