張?zhí)鞚桑?韓立國， 張 盼， 胡 勇， 毛 博

(吉林大學 地球探測科學與技術學院，長春 130012)

0 前言

全波形反演利用觀測地震記錄的波形信息來估計地下介質的物性參數(shù)(如密度，速度和各向異性參數(shù)等)，在油氣勘探，開發(fā)和地球動力學等領域有著很廣闊的應用前景[1]。隨著地下介質和大規(guī)?？茖W計算能力的提升，全波形反演的研究越來越受到重視，人們對反演所需要的優(yōu)化算法也有著越來越深入地研究。

20世紀80年代，TARANTOLA[2]提出了時間域全波形反演，并指出目標函數(shù)對參數(shù)模型的梯度可以通過計算殘差反傳波場與正傳波場的互相關得到,該方法被稱為伴隨狀態(tài)法，它有效地避免了Jacobi矩陣的計算。PRATT等[3-4]經(jīng)過推導得到FWI中的近似Hessian矩陣和精確Hessian矩陣，指出全波形反演可以利用基于Hessian矩陣的高斯-牛頓類優(yōu)化算法來進行反演。但當反演較大數(shù)據(jù)體時，高斯-牛頓類優(yōu)化算法的計算量很大，對存儲內存的要求很高，當Hessian矩陣高度病態(tài)或非正定時，反演存在不收斂的可能性。BROSSIER等[5]將l-BFGS算法應用到彈性波反演中。l-BFGS算法通過存儲固定數(shù)目的模型信息和梯度信息來計算Hessian矩陣，避免了高斯-牛頓類優(yōu)化算法需要求解Hessian矩陣逆的過程[6-8]，在各類優(yōu)化問題中得到了廣泛地應用。以上描述的優(yōu)化算法都要求在每次迭代中通過一定的步長搜索準則求解當前迭代所需要的步長，如Wolfe準則和Armijo準則等。在全波形反演的每次迭代中，這些準則要求程序不斷的試探步長[9]。當步長滿足特定的條件后才能進行模型的更新，而這些特定的條件通常要計算正演和求解梯度，所以依照準則的步長求解是全波形反演計算量大的重要原因之一。為了避免求解步長，王希云[10]提出了帶有固定步長的基于信賴域算法的步長計算公式。J.Sun,、X.D.Chen、Narushiama等[11-14]提出了梯度類優(yōu)化算法中的固定步長計算公式。馬巍和胡勇[15-16]提出了基于超記憶梯度優(yōu)化算法的固定步長計算公式。以上步長計算公式要求在反演的全過程固定步長進行模型迭代。在一般情況下，這些步長公式并不滿足全波形反演的要求。

為提高全波形反演的精度和計算效率，筆者將OBFGS(Online-BFGS)算法應用到時間域全波形反演。OBFGS算法不要求利用準則搜索步長，避免了求取步長時的梯度和目標函數(shù)計算。這里還提出了適合全波形反演的自適應步長計算公式，該計算公式避免了利用準則進行步長的搜索，簡化了步長求取的過程，能加快反演速率。筆者將這種自適應步長方法與OBFGS算法結合，應用到全波形反演中。通過數(shù)值反演，證明了基于自適應步長OBFGS優(yōu)化算法的全波形反演可以提高反演精度而且可以縮短反演所需要的時間，提高反演的計算效率。

1 方法以及原理

1.1 基于高斯-牛頓算法的全波形反演基本理論

時間域波動方程的離散形式為：

A(m)d=S

(1)

式中：d和S分別為波場和震源；m表示模型的速度參數(shù)；A(m)為正演矩陣。

二維時間域全波形反演的目標函數(shù)可以表示為：

(2)

式中：dcal為數(shù)值正演得到地震數(shù)據(jù)；dobs為實際觀測數(shù)據(jù)；‖g‖2表示L2范數(shù)。本文中的正演數(shù)據(jù)用時間二階空間八階有限差分算法計算得到[17-22]。

一般情況下全波形反演的模型更新公式定義為：

mt+1=mt+ηtpt

(3)

式中：t為迭代次數(shù)；pt和ηt分別為第次t迭代的更新方向和步長。對目標函數(shù)(式(2))進行二階泰勒展開:

(4)

上標T表示矩陣轉置，將式(4)進行最小化，可以得到高斯-牛頓法的更新方向所滿足的方程：

H(m1)pt=-gt

(5)

式中:gt表示目標函數(shù)對模型的一階導數(shù)，即梯度：

(6)

筆者利用正傳波場和殘差反傳波場的互相關求取梯度[2]。式(5)中H(mt)是目標函數(shù)對模型參數(shù)的二階導數(shù)，即Hessian矩陣：

(7)

所以更新方向可以表示為：

pt=-H(mt)-1gt

(8)

拓展到一般情況，模型的更新公式為：

mt+1=mt-αtH(mt)-1gt

(9)

式(9)為基于高斯-牛頓法的全波形反演模型迭代公式。

1.2 標準BFGS算法

由于對顯式Hessian矩陣及其逆矩陣的計算會消耗巨大內存。在反演時，為了避免直接計算Hessian矩陣及其逆矩陣，可以用一個近似矩陣來代替真實的Hessian矩陣。擬牛頓法就是應用了這種思想。在擬牛頓法中，通常用一個矩陣Bt來代替真實Hessian矩陣(t代表迭代次數(shù))。

標準BFGS算法是目前最流行，也是最有效的優(yōu)化算法之一。該算法由Broyden、Fletcher、Goldfarb和Shanno各自獨立提出[15]，隸屬于擬牛頓算法。標準BFGS算法的步驟如算法1(Algorithm1)所示。

1.3 步長搜索準則

在全波形反演的每次迭代中，求取正確的步長很關鍵。在傳統(tǒng)的全波形反演中通常依靠一定的步長搜索準則求取步長，這些準則通過逐步的試探步長求取最終步長。在得到一個步長之后，準則要求衡量目標函數(shù)和梯度以這個步長更新模型后，有沒有一個滿意的下降效果，若滿足這個下降效果，則用這個步長對模型進行更新，否則就把步長按照某種比率壓縮到滿足要求為止。

以Wolfe準則為例，Wolfe準則要求步長滿足以下條件：

(10)

式中:ρ∈(0,0.5);σ∈(ρ,1);ηt為當前步長。式(10)中第一個不等式判斷目標函數(shù)的下降情況，第二個不等式判斷梯度的下降情況。由式(10)可以看出，為確定合適的步長，每次算法都需要進行正演和梯度計算，這種依賴準則的步長搜索方法是全波形反演的計算量大的重要原因之一。

Algorithm1: Standard BFGS Method

Given:

·Current gradientgt)

·Line search line obeying Wolfe conditions

·Convergence tolerance ε>0

1.t:=0

2.B0=I

3. while ‖gt‖>ε

(a)pt=-Btf(m)

(b)ηt=line(f,m,pt)

(c)St=ηtpt

(d)mt+1=mt+St

(e)yt=f(mt+1)-f(mt+1)

(i)t:=t+1

4. return

--------------------------------------------Algorithm 2：Online BFGS Method

Given:

·Parameters 00

·Current gradientgt

·Sequence of step sizesηt>0

1.t:=0

2.Bo=?I

3. while not converged

(a)pt=-Btgt

(b)(no line search)

(c)St=(ηtpt)/c

(d)mt+1)=mt+St

(e)yt=gt+1-gt+λst

(i)t:=t+1

4. return

1.4 l-BFGS算法

l-BFGS算法同樣基于對真實海森矩陣進行代替的思想，l-BFGS算法利用前n次迭代過程中的梯度變化量和模型變化量去近似海森矩陣。在l-BFGS算法中 的更迭方程為：

(11)

l-BFGS算法可以在只存儲n+1個向量組的情況下，計算近似海森矩陣，可以加快反演的效率，但是l-BFGS算法犧牲了海森矩陣的精度。n的大小決定了在l-BFGS算法中近似海森矩陣B(n)的大小和精度[8]，因此n值的選取在一定程度上影響著l-BFGS 算法的最終反演精確度。此外，l-BFGS算法要求為海森矩陣提供一個近似的初始矩陣，初始矩陣的選擇同樣與l-BFGS算法的收斂性能緊密相關[6]。同BFGS算法一樣，基于Wolfe準則的l-BFGS算法依舊在選擇步長時需要進行梯度和目標函數(shù)的計算。這種依賴準則的步長求取方法會占用大量的內存，是傳統(tǒng)全波形反演計算量大而且占用內存多的重要原因之一。

1.5 OBFGS算法

OBFGS對標準BFGS算法中Bt的更新公式進行了修正，并在算法中引入了新的參數(shù)，達到了使算法不依賴準則搜索步長的效果。同時OBFGS算法引入了新的參數(shù)c、λ、?使得算法收斂，確保了Bt的正定性。OBFGS算法可以在合理選擇固定步長的前提下，使反演在這三個參數(shù)的控制下收斂。Nicol N. Schraudolph等[11]給出了這三個參數(shù)在OBFGS算法中的取值范圍，三個新參數(shù)的取值范圍已在算法2(Algorithm2)中表明，OBFGS算法和標準BFGS算法的區(qū)別已經(jīng)用橫線畫出。

在OBFGS算法中，如果?=1，λ=0，c=1，OBFGS算法中矩陣Bt的更新規(guī)則就會變成標準BFGS算法的更新法則?？梢钥闯鯫BFGS算法沒有依照準則進行步長的搜索，避免了在確定步長的過程中計算正演和梯度，因此可以提高反演效率。在本文的數(shù)值實驗中，新引入的三個參數(shù)的值為： ?=1，λ=1，c=1。

1.6 自適應步長設置策略

在本文中只有在目標函數(shù)不下降的情況下才會縮短步長，并以縮短后的步長直接進行模型更新進入下一次迭代，否則步長不變。

基于OBFGS算法的自適應步長全波形反演的步長設置：

1)設置步長的取值范圍：ηt∈ηmin,1)。

2)自適應步長迭代公式：

(12)

κ∈[0.5,1],κ為步長縮進因子。

限制步長小于“1”是為了防止步長過大使得算法不收斂。限定最小步長是因為在反演末期，如果采用過小的步長進行反演會使得反演進展得極為緩慢。將自適應步長公式帶入算法2的3(b)中就是基于自適應步長的OBFGS算法。

1.7 對比總結

通過以上對自適應步長OBFGS算法和BFGS算法的對比可以看出，自適應步長OBFGS算法不依賴步長收斂準則，可以避免在求取步長時計算梯度和目標函數(shù)，能降低反演的計算量，加快反演的計算效率。通過對自適應步長OBFGS算法和l-BFGS算法的對比可以看出，自適應步長OBFGS的海森矩陣精度更高，而且收斂效率不依賴初始海森矩陣，可以提高反演的精度。

2 數(shù)值實驗

2.1 Marmousi模型反演

筆者利用Marmousi模型進行反演測試，分別用OBFGS算法和l-BFGS算法進行對比實驗，突出本文算法的優(yōu)勢。真實模型如圖1(a)所示，初始模型由真實模型的線性平滑得到，如圖1(b)所示。模型的大小為50×200，網(wǎng)格的大小設置為40 m×40 m。震源函數(shù)為8 Hz主頻的雷克子波，設置20個炮點均勻的分布在模型的頂端，接收時間為1.9 s，采樣間隔為0.001 s，最大迭代次數(shù)為500次。

用基于l-BFGS算法的全波形反演進行對比，l-BFGS算法中的步長搜索方法為Wolfe準則(初始步長0.85)，OBFGS算法用自適應步長進行步長設置(初始步長設置為0.85，步長縮進因子設置為0.75)。兩種優(yōu)化算法反演的結果如圖2所示，由圖2可以看出，在不依賴準則搜索步長的情況下，在一定的誤差允許的范圍內，利用OBFGS算法得了更精確的反演結果。通過對比圖2中的黑色箭頭處的反演結果，可以看出OBFGS算法的反演結果更精確。圖3展示了反演結果第90道和100道處的速度反演曲線。從圖3中可以看出，基于OBFGS算法的全波形反演結果更接近真實速度。

2.2 反演效率

表1列出了可以評價算法效率的幾個參數(shù)[16]：(1)迭代次數(shù)；(2)進行全波形反演所需要的計算時間；(3)CPU平均使用率；(4)計算所需的最大內存；(5)最終的反演誤差。本文的反演在同一臺機器上

圖1 模擬實驗所用速度模型Fig.1 Model ased in numerical simulation test(a)真實模型；(b)初始模型

進行計算(32G內存，CORETMi3，4核處理器Ubuntu系統(tǒng))，程序用MATLAB進行編譯。

從對Marmousi模型進行反演的參數(shù)比較，可以看出當反演誤差為1.02%時，l-BFGS和OBFGS算法分別用時71 562 s和45 629 s(初始步長均為0.85)。OBFGS算法耗時更短，提高了反演的效率。

圖2 基于OBFGS和l-BFGS算法的全波形反演結果Fig.2 Full waveform inversion based on OBFGS methocl and L-BFGS method(a)基于OBFGS算法的全波形反演結果；(b)基于l-BFGS算法的全波形反演結果

圖3 基于OBFGS算法全波形反演和基于l-BFGS算法的全波形速度反演結果Fig.3 Velocity inversion result based on OBFGS and l-BFGS method(a)第90道；(b)第100道

優(yōu)化算法初始步長迭代次數(shù)計算時間/sCPU使用率/%最大內存/G誤差/%l-BFGS0.855007156251.56.21.02OBFGS0.854904562952.45.91.02OBFGS0.755414889851.75.51.02OBFGS0.55935513150.15.51.02OBFGS0.855004653452.15.90.99OBFGS0.755004696251.85.51.21

模型更新誤差計算公式為(|反演結果-真實模型|/真實模型)

從表1可以看出，在利用OBFGS算法進行反演時，如果固定反演的迭代次數(shù)，隨著初始步長的變小，最終的誤差會變大。利用不同初始步長進行基于OBFGS算法的全波形反演結果如圖4所示。由圖4可以看出，當初始步長設置較小時，最終得不到特別滿意的反演效果。這是因為在相同的迭代次數(shù)相同的前提下，初始步長設置太小模型將更新緩慢，從而使反演達不到理想的效果。

圖5展示了利用不同初始步長的基于OBFGS算法的全波形反演殘差曲線變化圖。由圖5可以看出，如果用不同的初始步長進行反演，殘差曲線在反演初期會呈現(xiàn)不同的下降趨勢。相對較大的步長殘差下降的更快，而小步長讓曲線下降的趨勢變慢。從圖5中可以看出，利用基于自適應步長OBFGS算法的全波形反演策略，可以使500次迭代反演全部收斂，是一種有效的反演算法。

圖4 不同初始步長反演結果Fig.4 Inversion result with different initial step(a)l-BFGS算法的反演結果；(b)初始步長設置為0.85的Marmousi模型反演結果;(c)初始步長設置為0.75的Marmousi模型反演結果；(d)初始步長設置為0. 5的Marmousi模型反演結果

圖5 基于不同初始步長的OBFGS算法的全波形反演殘差曲線Fig.5 The residual of FWI based on OBFGS with different initial step size

3 結論與討論

1)通過對OBFGS算法和l-BFGS與BFGS算法的理論對比可以看出，自適應步長OBFGS算法有著更高反演精度和計算效率。

2)通過Marmousi模型反演的結果分析表明，基于自適應步長的OBFGS算法的反演得到了更高精度的結果，提高了全波形反演的精度。

3)計算效率的對比結果表明，基于自適應步長的OBFGS算法計算效率更高。由于算法不依賴步長搜索準則進行步長求取，而是自適應的去計算步長，簡化了全波形反演中的步長求取過程，提高了反演的計算效率。

從表1可以看出，步長的設置過小會使反演進行的緩慢，所以初始步長會影響最終的反演結果。如何確定全波形反演優(yōu)化算法的眾多參數(shù)的取值范圍，從而使得反演的效率達到最高，依舊要進一步地研究。

[1] 劉國峰, 劉洪, 孟小紅,等. 頻率域波形反演中與頻率相關的影響因素分析[J].地球物理學報，2012, 55(4)：1345-1353. LIU G F, LIU H,MENG X H,et al. Frequency-related factors analysis in frequency domain waveform inversion [J]. Chinese Journal Geophysics, 2012, 55(4): 1345-1353.(In Chinese)

[2] TARANTOLA A. Inversion of seismic reflection data in the acoustic approximation [J]. Geophysics, 1984, 49(8):1259-1266.

[3] PRATT R G, SHIN C, HICKS. Gauss-newton and full newton methods in frequency-space seismic waveform inversion[J]. Geophysics. J. Int，1998(133): 341-362.

[4] PRATT R G. Inverse theory applied to multi-source cross-hole tomography I: acoustic wave-equation method [J]. Geophysics Prospecting, 1990, 38(3):287-310.

[5] BROSSIER R, OPERTOS, VIREUX J. Seismic imagine of complex onshore structures by 2D elastic frequency-domain full-waveform inversion [J]. Geophysics, 2009, 74(6):WCC105-WCC118.

[6] 王義,董良國. L-BFGS法時間域全波形反演中初始矩陣的選擇方法[J]. 石油物探, 2014, 53(5): 545-555. WANG Y, DONG L G. Selection strategy of the initial matrix for L-BFGS method in time domain full waveform inversion[J]. Geophysical Prospecting for Petroleum, 2014, 53(5): 545-555.(In Chinese)

[7] 王義, 董良國.基于截斷牛頓法的VTI介質聲波多參數(shù)全波形反演[J]. 地球物理學報, 2015,58(8): 2873-2885. WANG YI, DONG Liang-Guo.Multi-Parameter full waveform inversion foe acoustic VTI media using the truncated Newton method[J]. Chinese Journal Geophysics,2015,58(8):2873-2885.(In Chinese)

[8] 苗永康. 基于L-BFGS算法的時間域全波形反演[J]. 石油地球物理勘探, 2015, 50(3):469-474. MIAO Y K.Full waveform inversion in time domain based o limited-memory BFGS algorithm [J]. OGP,2015, 50(3):469-474.(In Chinese)

[9] KELLEY C T. Iterative methods for optimization [M]. New York: SIAM, 1999.

[10]王希云,仝建.帶有固定步長的非單調自適應信賴域算法[J].應用數(shù)學,2009,22(03):496-500. WANG X Y, TONG J.A nonmonotone adaptive trust region algorithm with fixed stepsize for unconstrained optimization problems[J]. Mathematica Applicata. 2009,22(03):496-500. (In Chinese)

[11]SCHRAUDOLPH N N, YU J, GüNTER S. A stochastic Quasi-newton method for online convex optimization [J]. Journal of Machine Learning Research, 2007(2):436-443.

[12]CHEN X, SUN J. Global convergence of a two-parameter family of conjugate gradient methods without line search[J]. Journal of Computational & Applied Mathematics, 2002, 146(1):37-45.

[13]YU Z. Global convergence of a memory gradient method without line search[J]. Journal of Applied Mathematics & Computing, 2008, 26(1-2):545-553.

[14]NARUSHIMA Y. A memory gradient method without line search for unconstrained optimization[J]. Sut J Math, 2006(42):191-206.

[15]馬巍. 無約束優(yōu)化問題的超記憶梯度法的若干研究 [D]. ?？冢汉Ｄ洗髮W，2013. MA W. Unconstrained optimization problem of super memory gradient method of research[D]. Haikou: Hainan University, 2013. (In Chinese)

[16]胡勇,韓立國,張盼,等. 混合超記憶梯度法多尺度全波形反演[J]. 石油物探, 2016, 55(4): 559-567. HU Y, HAN L G, ZHANG P, et al. Multi-scale full waveform inversion with hybrid super memory gradient method[J]. Geophysical Prospecting for Petroleum, 2016, 55(4): 559-567. (In Chinese)

[17]MOCZO P, KRISTEK J, GLIS M. The finite-difference modelling of earthquake motions. Waves and ruptures [J]. Cambridge Univ Pr, 2014:1-365.

[18]GAZDAG J. Modeling of the acoustic wave equation with transform methods[J]. Geophysics, 1981, 46(6): 854-859.

[19]KOSLOFF D, BAYSAL E. Forward modeling by a Fourier method [J]. Geophysics, 1982, 47(10): 1402-1412.

[20]RESHEF M, KOSLOFF D, EDWARDS M, et al. Three-dimensional elastic modeling by the Fourier method [J]. Geophysics, 1988, 53(9): 1184-1193.

[21]FORNBERG B. The pseudo spectral method; comparisons with finite differences for the elastic wave equation [J]. Geophysics, 1987, 52(4): 483-501.

[22]FORNBERG B. The pseudo spectral method; accurate representation of interfaces in elastic wave calculations [J]. Geophysics, 1988, 53(5): 625-637.