L-BFGS法時(shí)間域全波形反演中初始矩陣的選擇方法

王 義，董良國(guó)

(同濟(jì)大學(xué)海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092)

地震全波形反演(full waveform inversion，F(xiàn)WI)是近年來勘探地球物理領(lǐng)域的研究熱點(diǎn)之一。它利用疊前地震記錄的波形信息來重建高精度的地下介質(zhì)的物性參數(shù)(密度、速度和模量等)，可以為油氣勘探中的偏移成像提供更精確的背景速度[1]，可以應(yīng)用到儲(chǔ)層定量研究領(lǐng)域[2]，還可以應(yīng)用到全球局部和區(qū)域地殼和上地幔結(jié)構(gòu)等研究中[3]。因而，在油氣勘探開發(fā)和地球動(dòng)力學(xué)等領(lǐng)域有著廣闊的應(yīng)用前景。

20世紀(jì)80年代，Lailly[4]和Tarantola[5]利用最小二乘優(yōu)化的思想在時(shí)間域建立了FWI的核心框架。FWI希望找到最優(yōu)的地下介質(zhì)模型參數(shù)，使模擬地震記錄與觀測(cè)記錄達(dá)到最佳匹配。FWI涉及的數(shù)據(jù)量大、模型參數(shù)多，是典型的大規(guī)模優(yōu)化問題，通常采用局部迭代優(yōu)化方法求解[6]。目標(biāo)函數(shù)對(duì)模型參數(shù)的二階導(dǎo)數(shù)，即Hessian矩陣，對(duì)加速迭代方法的收斂起著重要作用[7]。但巨大的計(jì)算和存儲(chǔ)消耗使顯式獲取Hessian及其逆矩陣變得非常困難?？朔@個(gè)困難通常有三類方法:預(yù)條件共軛梯度法(preconditioned conjugate gradient method,PCG)、不精確Newton法和擬Newton法[6]。PCG法在FWI中獲得了大量應(yīng)用[8]，它通過使用預(yù)條件矩陣對(duì)梯度進(jìn)行校正，來加速迭代的收斂過程[9-13]。而預(yù)條件矩陣通常是先對(duì)Hessian矩陣做一些簡(jiǎn)單近似(如對(duì)角近似或塊對(duì)角近似)，然后求逆得到預(yù)條件矩陣。不精確Newton法在每步迭代中通過采用PCG法近似求解Newton方程來獲取本次迭代的更新方向。如果沒有好的預(yù)條件矩陣，每步迭代中的PCG 求解過程收斂很慢，導(dǎo)致不精確Newton法的整體計(jì)算效率偏低[6]。擬Newton法是求解無約束優(yōu)化問題的所有利用一階導(dǎo)數(shù)信息的方法中最有效的一類方法[14]，其中的L-BFGS法在FWI中得到了廣泛應(yīng)用[15-17]。

L-BFGS法不需要顯式形成Hessian及其逆矩陣，在非線性優(yōu)化方面表現(xiàn)出比CG法收斂快且穩(wěn)健的優(yōu)勢(shì)[18]。該方法在每步迭代中都需要提供Hessian逆矩陣的一個(gè)初始近似矩陣。這個(gè)矩陣可以隨迭代的進(jìn)行而更新，也可以固定不變。初始矩陣的選取和更新方式對(duì)L-BFGS法的收斂性能有重要影響。但初始矩陣的選取并沒有統(tǒng)一標(biāo)準(zhǔn)，需要根據(jù)具體優(yōu)化問題而定[6]。如前文所述，PCG法采用的預(yù)條件矩陣即是對(duì)Hessian逆矩陣的近似。因而，選擇初始矩陣時(shí)可以借鑒預(yù)條件矩陣的作法。此外，L-BFGS法需要存儲(chǔ)一定數(shù)目的梯度殘差和模型殘差向量對(duì)，存儲(chǔ)數(shù)目的多少也會(huì)對(duì)其收斂性能產(chǎn)生影響。

針對(duì)上述問題，我們分析了幾種時(shí)間域可行的預(yù)條件矩陣的原理，對(duì)采用它們作為初始矩陣的L-BFGS法的反演性能進(jìn)行了對(duì)比分析，并與相應(yīng)的PCG法比較。此外，還分析了存儲(chǔ)向量對(duì)的數(shù)目對(duì)L-BFGS法的影響。最后通過時(shí)間域的Marmousi模型FWI反演試驗(yàn)，來說明不同的初始矩陣選擇以及更新方式對(duì)L-BFGS法收斂性能的影響。

1 時(shí)間域全波形反演理論簡(jiǎn)介

時(shí)間域波動(dòng)方程的離散形式為

(1)

式中：u,S分別為波場(chǎng)和震源;m為模型參數(shù)，維數(shù)M;A(m)為正演矩陣?？梢钥闯觯珹-1的每一列都對(duì)應(yīng)于特定震源位置和特定激發(fā)時(shí)刻的Green函數(shù)。

設(shè)最小二乘目標(biāo)函數(shù)為

(2)

式中：ns為炮數(shù)；ds是單炮觀測(cè)記錄；Rs為提取記錄的算子。Newton法的更新方向滿足

(3)

(4)

下面采用虛震源給出Jacobi矩陣的構(gòu)造原理。分別在(1)式兩邊對(duì)m求導(dǎo)可得

(5)

定義F為虛震源，它的每一列代表該位置處模型參數(shù)的虛震源。虛震源中除了包含波場(chǎng)的作用外，還涉及到傳播算子對(duì)模型參數(shù)的影響?A/?m，描述了數(shù)據(jù)對(duì)模型參數(shù)的隨方位角變化的敏感性[19]。以F的每列為震源分別做正演模擬可以得到偏導(dǎo)波場(chǎng)，進(jìn)而構(gòu)造出Jacobi矩陣。因而，Gauss-Newton的近似Hessian矩陣可以寫為

(6)

2 L-BFGS方法原理和初始矩陣的幾種選擇方法

L-BFGS法利用目標(biāo)函數(shù)的一階導(dǎo)數(shù)信息遞歸地計(jì)算近似Hessian逆矩陣與梯度的乘積，來得到迭代更新方向，克服了顯式求取Hessian矩陣及其逆的困難，在FWI中應(yīng)用廣泛。

2.1 L-BFGS優(yōu)化方法原理

每步迭代中近似Hessian逆矩陣與梯度向量乘積的具體計(jì)算公式為

(7)

(8)

這種初始矩陣在每次迭代中都可以更新[6]。以上2種方式并沒有具體結(jié)合FWI中Hessian矩陣的特點(diǎn)。下面我們討論幾種針對(duì)時(shí)間域FWI的初始矩陣的選擇方式。

2.2 L-BFGS法中初始矩陣的幾種選擇方式

預(yù)條件矩陣通常是由對(duì)Hessian矩陣做簡(jiǎn)單的近似，如對(duì)角近似或塊對(duì)角近似，再求逆后得到。因而，初始矩陣的選取可以借鑒預(yù)條件矩陣的做法。在頻率域也可以對(duì)Hessian矩陣做塊對(duì)角近似[21]，但在時(shí)間域計(jì)算塊對(duì)角Hessian矩陣的運(yùn)算量很大，所以一般并不可行。下面分析幾種在時(shí)間域可行的對(duì)角預(yù)條件矩陣方式。

2.2.1 Pseudo Hessian(PH方式)

Shin等人[9]定義FTF為Pseudo Hessian，并用對(duì)角矩陣diag(FTF)來近似Hessian矩陣。從虛震源F的定義可知，這種方式只考慮了震源波場(chǎng)的作用和傳播算子對(duì)模型參數(shù)的影響。它用單位矩陣替代(6)式的矩陣B，忽略了所有檢波點(diǎn)處的Green函數(shù)對(duì)目標(biāo)函數(shù)的影響。

2.2.2 New Pseudo Hessian(newPH方式)

這種方式是對(duì)上一種方式的改進(jìn)，由Choi等人在頻率域提出[10]。代替矩陣B的仍是對(duì)角矩陣，對(duì)角元素的值則是所有震源位置的Green函數(shù)在該空間點(diǎn)處的波場(chǎng)的振幅之和。如果觀測(cè)系統(tǒng)中炮點(diǎn)和檢波點(diǎn)有重合的位置，則這種方式能考慮重合部分的檢波點(diǎn)處的Green函數(shù)的作用。在時(shí)間域可以采用類似的做法，用所有震源波場(chǎng)在地下同一點(diǎn)的振幅之和，作為代替矩陣B的對(duì)角矩陣在該位置的對(duì)角元素。

2.2.3 采用震源正傳波場(chǎng)和記錄殘差反傳波場(chǎng)的能量矩陣(RcdIW方式)

該方式采用震源正傳波場(chǎng)與以模擬和觀測(cè)記錄的殘差作為震源的反傳波場(chǎng)的能量分布來構(gòu)造梯度的預(yù)條件矩陣，目的是消除正、反傳播過程中的幾何擴(kuò)散效應(yīng)[11]。記u(x;xs)為xs處的震源的正傳波場(chǎng)，λ(x;xr)為記錄殘差反傳的波場(chǎng)，那么，正傳波場(chǎng)的能量矩陣為

(9)

反傳波場(chǎng)的能量矩陣為

(10)

對(duì)梯度進(jìn)行預(yù)條件:

(11)

嚴(yán)格地講，這種方式并不是對(duì)Hessian逆矩陣的近似，而是從波場(chǎng)能量分布的角度來對(duì)梯度進(jìn)行校正。它考慮了記錄殘差反傳波場(chǎng)的作用，但也攜帶了散射點(diǎn)的信息。采用它來預(yù)條件梯度會(huì)破壞模型擾動(dòng)的聚焦。另外，該方式?jīng)]有考慮傳播算子對(duì)模型參數(shù)的影響。

2.2.4 采用震源正傳波場(chǎng)和模擬記錄反傳波場(chǎng)的能量矩陣(ResIW方式)

該方法是對(duì)上一種方式的改進(jìn)[12]，不同之處在于采用的是模擬記錄反傳波場(chǎng)的能量分布，不會(huì)對(duì)模型擾動(dòng)的聚焦造成破壞。但該方式同樣沒有考慮傳播算子對(duì)模型參數(shù)的影響。它需要增加一次波場(chǎng)模擬來計(jì)算模擬記錄的反傳波場(chǎng)v(x;xr)。正傳波場(chǎng)的能量矩陣仍采用(9)式。反傳波場(chǎng)的能量矩陣以及預(yù)條件后的梯度分別為

2.2.5 用Hessian-vector 乘積構(gòu)造對(duì)角矩陣

Korta等[13]在雙參數(shù)FWI中采用這種方式構(gòu)造塊狀對(duì)角的Hessian近似。在單參數(shù)FWI中，可以使用有限差分法[6]、二階伴隨狀態(tài)法[22]或基于Born核函數(shù)的矩陣分解方法計(jì)算出Hessian-vector乘積Hv，其中v=(1,1,…,1,1)T。以該乘積為對(duì)角元素構(gòu)造一個(gè)對(duì)角矩陣。此時(shí)，每個(gè)對(duì)角元素均為Hessian矩陣該行所有元素的疊加。如

果Hessian矩陣是嚴(yán)格對(duì)角占優(yōu)的話，即對(duì)角線元素的絕對(duì)值大于同一行所有其它元素的絕對(duì)值之和，經(jīng)過按行疊加得到的元素會(huì)偏離Hessian矩陣的主對(duì)角元素的精確值。如果不是嚴(yán)格對(duì)角占優(yōu)，計(jì)算出的乘積很有可能是零或負(fù)值。即使取絕對(duì)值的話，數(shù)值試驗(yàn)表明其反演過程容易出現(xiàn)不穩(wěn)定的現(xiàn)象。所以，在單參數(shù)FWI中，這種方式并不是合適的預(yù)條件矩陣。

3 數(shù)值反演試驗(yàn)

我們通過Marmousi模型合成數(shù)據(jù)的FWI反演數(shù)值試驗(yàn)，對(duì)比采用上述前4種方式作為初始矩陣的L-BFGS法的反演性能，同時(shí)與對(duì)應(yīng)的預(yù)條件共軛梯度法進(jìn)行比較。

從原始Marmousi模型中抽取了5940m×2980m的模型，實(shí)際模型(局部)和初始模型如圖1 所示。網(wǎng)格剖分為298×150，網(wǎng)格間距20m。共使用56炮，間距為100m,檢波器298個(gè)，間距為20m。時(shí)間采樣間隔2ms，記錄長(zhǎng)度3s。震源采用主頻7Hz的Ricker子波。采用空間4階、時(shí)間2階有限差分法求解常密度聲波方程，四周均為PML吸收邊界[23]。

圖1 實(shí)際Marmousi模型(a)和初始模型(b)

3.1 殘差向量對(duì)的數(shù)目對(duì)L-BFGS性能的影響

L-BFGS法構(gòu)造近似Hessian逆矩陣時(shí)，需要存儲(chǔ)r對(duì)梯度殘差和模型殘差。r越大，構(gòu)造近似Hessian矩陣用到的信息越豐富，但所需內(nèi)存也隨之增大。假設(shè)模型向量的長(zhǎng)度為M，則存儲(chǔ)這些殘差向量需要2M×r個(gè)浮點(diǎn)類型(float類型)的內(nèi)存。首先，我們測(cè)試了L-BFGS法中需要存儲(chǔ)的向量對(duì)的數(shù)目r對(duì)其性能的影響。

分別取r為 4，8，12，16和20，初始矩陣采用公式(8)的方式，對(duì)應(yīng)的目標(biāo)函數(shù)收斂曲線如圖2。全波形反演的主要計(jì)算消耗為正演模擬，所以用目標(biāo)函數(shù)隨正演次數(shù)的變化來描述優(yōu)化方法的收斂速度。可以看出，隨著r的增大，L-BFGS法的收斂速度也加快。r=4時(shí)收斂最慢，r為8～20時(shí)收斂速度逐漸加快，但提升并不明顯。當(dāng)r=20時(shí)收斂最快，但存儲(chǔ)殘差向量的內(nèi)存是r=8時(shí)的2.5倍。圖3分別給出了r=8時(shí)的第324代(972次正演)和r=20時(shí)的第291代(873次正演)反演結(jié)果。此時(shí)，二者的目標(biāo)函數(shù)均下降到1.01×10-3附近。圖4分別對(duì)比了圖3中反演結(jié)果在地表3km和5km處的縱向速度曲線。從圖3 和圖4可以看出，2個(gè)反演結(jié)果非常接近。因而，綜合考慮收斂速度和內(nèi)存消耗，參數(shù)r的取值無需過大。在此選取r=8。

圖2 存儲(chǔ)殘差向量對(duì)的數(shù)目r對(duì)L-BFGS法目標(biāo)函數(shù)收斂的影響

圖3 存儲(chǔ)不同向量對(duì)數(shù)目r的L-BFGS法的反演結(jié)果a r=8,第324代; b r=20,第291代

圖4 對(duì)應(yīng)圖3中反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對(duì)比

3.2 不同初始矩陣情況下的L-BFGS法的反演性能對(duì)比

我們分別選擇前4種方式作為L(zhǎng)-BFGS法的初始矩陣,即依次為PH,newPH,RcdIW和ResIW。初始矩陣在迭代過程中可以固定不變，也可以每次迭代都更新(Up1)。另外，一并測(cè)試單位矩陣(Identity)和公式(8)對(duì)應(yīng)的矩陣(Nocedal)的情況。圖5給出了不同初始矩陣、不同更新方式情況下的L-BFGS法的目標(biāo)函數(shù)收斂曲線。

從圖5可以看出，不同初始矩陣、不同更新方式的L-BFGS法的反演性能存在較大差異。按收斂速度快慢大致可以分為3類。

第一類收斂最快，包含以下5種方式。初始矩陣選取PH和newPH且固定不變時(shí)(藍(lán)色和粉色虛線)收斂最快，二者非常接近。選RcdIW且固定不變時(shí)(綠色虛線)和選PH每次迭代都更新時(shí)(藍(lán)色點(diǎn)劃線)也比較接近，相比于前兩種方式則稍微慢一些。選newPH且每代都更新時(shí)(粉色點(diǎn)劃線)在開始階段較慢，之后收斂加快。

選RcdIW且每代都更新時(shí)(綠色點(diǎn)劃線)與選ResIW固定不變時(shí)(黑色虛線)收斂速度比較接近，屬于第二類收斂較快的方式。

選ResIW且每代都更新(黑色點(diǎn)劃線)，Nocedal方式(紅色實(shí)線)和Identity(青色實(shí)線)屬于收斂比較慢的第三類。其中，單位矩陣并未提供有關(guān)Hessian的信息，收斂最慢。

綜合以上對(duì)比結(jié)果可知，PH和newPH作初始矩陣且固定不變時(shí)，都可以快速收斂。RcdIW方式稍慢于前兩者，可能與其未考慮正演算子對(duì)模型參數(shù)的影響有關(guān)。因存在不利于模型擾動(dòng)聚焦的因素，ResIW在4種方式中收斂最慢。另外一個(gè)結(jié)論是，每種方式固定不變時(shí)(虛線)比每次迭代都更新(點(diǎn)劃線)的收斂要快。這說明每次迭代更新初始矩陣所增加的計(jì)算消耗并未能帶來目標(biāo)函數(shù)的快速下降。

下面對(duì)比反演所得的速度模型。表1列出了

目標(biāo)函數(shù)(歸一化后)下降到1.01×10-3左右時(shí)，圖5中各種方式的L-BFGS法所需要的迭代次數(shù)和正演次數(shù)的對(duì)比，其中單位矩陣(Identity)方式收斂最慢，只下降到1.202×10-3。圖6是PH方式第88代和RcdIW方式第97代的反演結(jié)果，對(duì)應(yīng)的目標(biāo)函數(shù)值分別為1.016×10-3和1.011×10-3,正演次數(shù)則為264次和292次。圖7給出了這兩個(gè)模型在地表3km和5km處的縱向速度曲線的對(duì)比?？梢钥闯?，兩種方式的反演結(jié)果非常接近，與真實(shí)模型也比較吻合。右下方包含高速體的區(qū)域則都存在偏差，對(duì)真實(shí)Hessian矩陣逼近的程度不夠可能是原因之一。圖7中還給出了PH方式第34代(102次正演)反演的速度曲線(粉色虛線)，對(duì)應(yīng)的目標(biāo)函數(shù)值為1.03×10-3。此時(shí)的反演模型與真實(shí)模型有一定偏差，說明仍需迭代更新。從這些對(duì)比可以看出，當(dāng)目標(biāo)函數(shù)下降到基本相同的水平時(shí)，PH和RcdIW方式的反演結(jié)果也非常接近。而且，PH方式充當(dāng)初始矩陣時(shí)在計(jì)算效率上占有少量?jī)?yōu)勢(shì)。

圖5 不同初始矩陣和不同更新方式的L-BFGS法的目標(biāo)函數(shù)收斂曲線

NocedalPHnewPHRcdIWResIW目標(biāo)函數(shù)值1.017×10-31.016×10-31.039×10-31.011×10-31.018×10-3正演次數(shù)966264250292507迭代次數(shù)322888397169IdentityPH-Up1newPH-Up1RcdIW-Up1ResIW-Up1目標(biāo)函數(shù)值1.202×10-31.015×10-31.007×10-31.019×10-31.011×10-3正演次數(shù)999336336464849迭代次數(shù)33311284116283

圖6 PH方式第88代(a)和RcdIW方式第97代(b)反演的速度模型

圖8對(duì)比了PH方式第88代和newPH方式第83代反演結(jié)果的速度曲線。此時(shí)它們分別運(yùn)行了264次和250次正演，效率相當(dāng)，反演結(jié)果也非常相似。圖9是PH方式第88代和ResIW方式第169代反演結(jié)果的速度曲線對(duì)比，可以看出反演結(jié)果比較接近。但所需正演次數(shù)分別為264次和507次。所以ResIW方式的計(jì)算效率偏低。

下面分析初始矩陣每次迭代都更新時(shí)的反演結(jié)果。圖10是PH-Up1第112代(336，表示正演次數(shù)，下同)和PH第88代(264)反演結(jié)果的速度曲線對(duì)比。圖11是RcdIW-Up1第116代(464)和RcdIW 第97代(292)的對(duì)比。可以看出，對(duì)于每種方式來說，當(dāng)目標(biāo)函數(shù)下降到基本相同的水平時(shí)，初始矩陣固定和每次迭代都更新這兩種做法的反演結(jié)果很接近。但初始矩陣固定時(shí)計(jì)算效率更高。對(duì)于newPH和ResIW兩種方式，結(jié)論相同。

綜上可知，初始矩陣不變時(shí)，PH和newPH方式收斂最快，計(jì)算效率和反演結(jié)果均相當(dāng);RcdIW方式效率稍慢一些，ResIW方式則收斂較慢。以上4種方式作初始矩陣且每次迭代都更新時(shí)，L-BFGS方法的收斂速度要慢于初始矩陣固定不變的情況。因而，選好初始矩陣后，固定不變即可。初始矩陣選Nocedal方式和單位矩陣時(shí)收斂很慢。

圖7 對(duì)應(yīng)圖6中反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對(duì)比

圖8 PH方式和newPH方式反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對(duì)比

圖9 PH方式和ResIW方式反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對(duì)比

圖10 PH-Up1方式和PH方式反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對(duì)比

圖11 RcdIW-Up1方式和RcdIW方式反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對(duì)比

此外，L-BFGS法構(gòu)造近似Hessian逆矩陣時(shí)使用的是梯度殘差向量。分別嘗試采用PH，newPH，RcdIW和ResIW 4種方式預(yù)條件后的梯度殘差向量來構(gòu)造近似Hessian逆矩陣。如此進(jìn)行的L-BFGS法FWI試驗(yàn)的收斂速度慢于上述作為初始矩陣的情況，而且容易導(dǎo)致反演不穩(wěn)定。所以，不推薦對(duì)L-BFGS法中的梯度進(jìn)行額外的預(yù)條件。

3.3 不同預(yù)條件CG法的反演性能對(duì)比

分別采用上述4種方式的預(yù)條件共軛梯度(CG)法進(jìn)行同樣的全波形反演試驗(yàn)。圖12給出了4種預(yù)條件CG法與之前不同初始矩陣(固定)L-BFGS法的目標(biāo)函數(shù)收斂曲線的對(duì)比?？梢钥闯?，在優(yōu)化過程前期4種預(yù)條件CG法能達(dá)到和L-BFGS法比較接近的收斂速度，但都逐漸變慢。其中，ResIW預(yù)條件CG法收斂最慢。其它3種預(yù)條件CG法也都慢于初始矩陣固定的PH，newPH和RcdIW的L-BFGS法，但優(yōu)于初始矩陣為單位矩陣和Nocedal的方式。所以，在全波形反演中選取適當(dāng)?shù)某跏季仃嚕梢允筁-BFGS法的性能優(yōu)于預(yù)條件CG法。這與Brossier等在頻率域使用PH方式進(jìn)行FWI的結(jié)論一致[16]。

表2列出了不做預(yù)條件的CG法和上述4種預(yù)條件CG法的反演過程中的一組數(shù)據(jù)對(duì)比，主要包括目標(biāo)函數(shù)值、所需迭代次數(shù)以及正演次數(shù)。圖13a 給出了PH方式的預(yù)條件CG法第258代反演的速度模型，此時(shí)正演次數(shù)為774次，目標(biāo)函數(shù)值約為1.012×10-3。圖13b則給出了RcdIW方式預(yù)條件CG法第203次迭代的結(jié)果，目標(biāo)函數(shù)下降為1.014×10-3，正演次數(shù)為812次。

為了與采用對(duì)應(yīng)方式作初始矩陣的L-BFGS法比較(圖6)，我們抽取了相應(yīng)的速度曲線。圖14 顯示了PH方式的預(yù)條件CG法和作為L(zhǎng)-BFGS法初始矩陣時(shí)反演結(jié)果的速度對(duì)比曲線。圖15則是RcdIW方式的對(duì)比?？梢钥闯觯瑘D6與圖13對(duì)應(yīng)的反演速度模型是非常接近的。從計(jì)算消耗上來說，初始矩陣為PH和RcdIW的L-BFGS法分別進(jìn)行了264次和292次正演?？梢娖溆?jì)算效率高于相應(yīng)的預(yù)條件CG法。對(duì)于newPH方式也可以得到同樣的結(jié)論，這里不再給出具體的圖像和曲線。不做預(yù)條件的CG法和ResIW方式的預(yù)條件CG法的目標(biāo)函數(shù)分別下降到1.604×10-3和1.948×10-3時(shí)，已經(jīng)消耗掉約1000次正演，可見其效率偏低。

綜合以上對(duì)比分析可以看出，盡管使用預(yù)條件矩陣可以提高CG法的效率，但仍不及采用同樣方式作初始矩陣的L-BFGS法。

圖12 不同預(yù)條件CG法與不同初始矩陣(固定)L-BFGS法的目標(biāo)函數(shù)收斂曲線對(duì)比

CGPreCG-PHPreCG-newPHPreCG-RcdIWPreCG-ResIW目標(biāo)函數(shù)值1.604×10-31.012×10-31.038×10-31.014×10-31.948×10-3正演次數(shù)999774698812999迭代次數(shù)333258174203333

圖13 預(yù)條件CG法的反演結(jié)果a PH方式第258代； b RcdIW方式第203代

圖14 PH方式預(yù)條件CG法和充當(dāng)H0且固定不變時(shí)的L-BFGS法的反演結(jié)果在地表3km(a)和5km(b)處的vP曲線

圖15 RcdIW方式預(yù)條件CG法和充當(dāng)H0且固定不變的L-BFGS法的反演結(jié)果在地表3km(a)和5km(b)處的vP曲線

4 結(jié)束語

我們討論了全波形反演中L-BFGS法的有效使用問題，主要分析了初始矩陣的選取方式和更新方式對(duì)L-BFGS法性能的影響。通過時(shí)間域Marmousi模型反演試驗(yàn)的比較，可以得出幾點(diǎn)結(jié)論。

1) 初始矩陣固定不變時(shí)，只采用虛震源的PH方式和增加了震源波場(chǎng)振幅的newPH方式的L-BFGS方法收斂速度相當(dāng)且最快。采用震源正傳波場(chǎng)和模擬記錄反傳波場(chǎng)能量的RcdIW方式收斂稍慢。采用震源正傳波場(chǎng)和記錄殘差反傳波場(chǎng)能量的ResIW方式收斂明顯慢于前3種方式，而使用Nocedal方式和單位矩陣的L-BFGS法收斂最慢。另外，在單參數(shù)FWI中，采用Hessian向量乘積構(gòu)造對(duì)角矩陣的方式并不適合作為初始矩陣。

2) PH，newPH，RcdIW和ResIW4種方式作為初始矩陣且固定不變的L-BFGS法的收斂速度快于同種方式但初始矩陣每次迭代都更新的情形。

3) 與預(yù)條件CG法相比，采用同種方式作初始矩陣的L-BFGS方法收斂更快。

4) L-BFSG法的收斂速度會(huì)隨存儲(chǔ)的殘差向量對(duì)的數(shù)目的增大而有所提升，但所需內(nèi)存也會(huì)明顯增多。因而綜合計(jì)算效率和內(nèi)存消耗的考慮，存儲(chǔ)數(shù)目無需過大(如8左右)。

參 考 文 獻(xiàn)

[1] Boonyasiriwat C,Schuster G T,Paul V,et al.Applications of multiscale waveform inversion to marine data using a flooding technique and dynamic early-arrival windows[J].Geophysics,2010,75(6):R129-R136

[2] 單蕊，卞愛飛，於文輝，等.利用疊前全波形反演進(jìn)行儲(chǔ)層預(yù)測(cè)[J].石油地球物理勘探,2011,46(4):629-633

Shan R,Bian A F,Yu W H,et al.Pre-stack full waveform inversion for reservoir prediction[J].Oil Geophysical Prospecting,2011,46(4):629-633

[3] Fichtner A,Kennett B L N,Igel H,et al.Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods[J].Geophysical Journal International,2009,179(3):1703-1725

[4] Lailly P.The seismic inverse problem as a sequence of before stack migrations[C]∥Conference on inverse scattering:theory and application.NewYork:SIAM,1983:206-220

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

[6] Nocedal J,Wright S J.Numerical optimization[M].2nd ed.New York:Springer,2006:101-184

[7] Pratt R G,Shin C,Hicks G J.Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J].Geophysical Journal International,1998,133(2):341-362

[8] Virieux J,Operto S.An overview of full-waveform inversion in exploration geophysics[J].Geophysics,2009,74(6):WCC1-WCC26

[9] Shin C,Jang S,Min D J.Improved amplitude preservation for prestack depth migration by inverse scattering theory[J].Geophysical Prospecting,2001,49(5):592-606

[10] Choi Y,Min D J,Shin C.Frequency-domain elastic full waveform inversion using the new pseudo-Hessian matrix:experience of elastic Marmousi-2 synthetic data[J].Bulletin of the Seismological Society of America,2008,98(5):2402-2415

[11] Zhang Z,Lin Y,Huang L.Full-waveform inversion in the time domain with an energy weighted gradient[J].Expanded Abstracts of 81stAnnual Internet SEG Mtg,2011,2772-2776

[12] Zhang Z,Lin Y,Huang L.A wave-energy-based precondition approach to full-waveform inversion in the time domain[J].Expanded Abstracts of 82ndAnnual Internet SEG Mtg,2012,1-5

[13] Korta N,Fichtner A,Sallaers V.Block-diagonal approximate hessian for Preconditioning in full waveform inversion[J].Expanded Abstracts of 75thEAGE Conference,2013,Tu-P04-16

[14] 袁亞湘.非線性優(yōu)化計(jì)算方法[M].北京:科學(xué)出版社，2008:82-88

Yuan Y X.Computation Methods for Nonlinear Optimization[M].Beijing:Science Press,2008:82-88

[15] Guitton A,Symes W W.Robust inversion of seismic data using the Huber norm[J].Geophysics,2003,68(4):1310-1319

[16] Brossier R,Operto S,Virieux J.Seismic imaging of complex structures by 2D elastic frequency domain full-waveform inversion[J].Geophysics,2009,74(6):WCC105-WCC118

[17] Prieux V,Brossier R,Operto S,et al.Multi-parameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field.Part 1:imaging compressional wavespeed,density and attenuation[J].Geophysical Journal International,2013,194(3):1640-1664

[18] Kelley C T.Iterative methods for optimization[M].New York:SIAM,1999:71-83

[19] Gholami Y,Ribodetti A,Brossier R,et al.Sensitivity analysis of full waveform inversion in VTI media[J].Expanded Abstracts of 72ndEAGE Conference,2010,P371

[20] Plessix R E,Mulder W A.Frequency-domain finite-difference amplitude preserving migration[J].Geophysical Journal International,2004,157(3):975-987

[21] Hu W,Abubabak A,Habashy T M,et al.Preconditioned non-linear conjugate gradient method for frequency domain full waveform seismic inversion[J].Geophysical prospecting,2011,59(3):477-491

[22] Fichtner A,Trampert J.Hessian kernels of seismic data functionals based upon adjoint techniques[J].Geophysical Journal International,2009,185(2):775-798

[23] Berenger J P.A perfectly matched layer for the absorption of electromagnetic waves[J].Journal of Computational Physics，1994,114(2):185-200