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

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

0 前言

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

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

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

1 方法以及原理

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

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

A(m)d=S

(1)

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

二維時間域全波形反演的目標(biāo)函數(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迭代的更新方向和步長。對目標(biāo)函數(shù)(式(2))進(jìn)行二階泰勒展開:

(4)

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

H(m1)pt=-gt

(5)

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

(6)

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

(7)

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

pt=-H(mt)-1gt

(8)

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

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

(9)

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

1.2 標(biāo)準(zhǔn)BFGS算法

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

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

1.3 步長搜索準(zhǔn)則

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

以Wolfe準(zhǔn)則為例，Wolfe準(zhǔn)則要求步長滿足以下條件：

(10)

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

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算法同樣基于對真實海森矩陣進(jìn)行代替的思想，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算法的收斂性能緊密相關(guān)[6]。同BFGS算法一樣，基于Wolfe準(zhǔn)則的l-BFGS算法依舊在選擇步長時需要進(jìn)行梯度和目標(biāo)函數(shù)的計算。這種依賴準(zhǔn)則的步長求取方法會占用大量的內(nèi)存，是傳統(tǒng)全波形反演計算量大而且占用內(nèi)存多的重要原因之一。

1.5 OBFGS算法

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

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

1.6 自適應(yīng)步長設(shè)置策略

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

基于OBFGS算法的自適應(yīng)步長全波形反演的步長設(shè)置：

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

2)自適應(yīng)步長迭代公式：

(12)

κ∈[0.5,1],κ為步長縮進(jìn)因子。

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

1.7 對比總結(jié)

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

2 數(shù)值實驗

2.1 Marmousi模型反演

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

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

2.2 反演效率

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

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

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

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

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

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

優(yōu)化算法初始步長迭代次數(shù)計算時間/sCPU使用率/%最大內(nèi)存/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

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

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

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

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

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

3 結(jié)論與討論

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

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

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

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

[1] 劉國峰, 劉洪, 孟小紅,等. 頻率域波形反演中與頻率相關(guān)的影響因素分析[J].地球物理學(xué)報，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介質(zhì)聲波多參數(shù)全波形反演[J]. 地球物理學(xué)報, 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]王希云,仝建.帶有固定步長的非單調(diào)自適應(yīng)信賴域算法[J].應(yīng)用數(shù)學(xué),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(xué)，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.