基于T-matrix的非線性參數(shù)估計方法

王本鋒, 吳如山, 陳小宏, 陸文凱

1 清華大學自動化系，智能技術與系統(tǒng)國家重點實驗室, 北京　100084　2 中國石油大學(北京)油氣資源與探測國家重點實驗室， 北京　102249　3 Modeling and Imaging Laboratory, University of California, Santa Cruz，95064, U.S.A.

?

基于T-matrix的非線性參數(shù)估計方法

王本鋒1,2,3, 吳如山3, 陳小宏2, 陸文凱1

1 清華大學自動化系，智能技術與系統(tǒng)國家重點實驗室, 北京1000842 中國石油大學(北京)油氣資源與探測國家重點實驗室， 北京1022493 Modeling and Imaging Laboratory, University of California, Santa Cruz，95064, U.S.A.

摘要全波形反演可提供高精度的地下介質(zhì)參數(shù)空間分布，但傳統(tǒng)的全波形反演方法建立在Born近似的基礎上，對初始模型具有一定的依賴性.為了擺脫Born近似的束縛，本文基于二維常密度聲波方程，在De Wolf近似的前提下，借助傳輸矩陣(T-matrix)方法，深入研究了逆薄板傳播算子(Inverse Thin-Slab Propagator, ITSP)，實現(xiàn)了速度擾動的非線性估計.ITSP方法避免了Born級數(shù)方法在擾動較強、擾動區(qū)域較大時的發(fā)散性問題，且只經(jīng)過一次掃描校正，計算效率較高.二維模擬數(shù)據(jù)分析驗證了本文方法的可行性以及有效性.

關鍵詞全波形反演； De Wolf 近似； 傳輸矩陣； 逆薄板傳播算子； Born級數(shù)

The Inverse Thin-Slab Propagator (ITSP) is studied in detail both theoretically and numerically based on Transmission Matrix (T-matrix) and De Wolf approximation using the 2D acoustic wave equation with constant density, which can overcome the limitations of Born approximation. The ITSP method considers all orders of scattering effects and has no initial model dependence, besides, it is a regularized method in each correction step and has no divergence effects. One correction sweep (half from upper part, half from lower part) is involved and is computationally efficient. Based on the ITSP method, the velocity perturbation can be achieved and the velocity distribution can be obtained finally with the background velocity.

We designed three numerical models including ellipse ball with positive velocity perturbations, Gaussian ball with negative velocity perturbations and ellipse ball with both positive and negative velocity perturbations. With the ITSP method, we achieve the velocity distribution based on the known T-matrix and the T-matrix estimation method will be left for future research. The reconstructed results are consistent with the true velocity distributions and the relative errors are minor enough for these three models, which demonstrate the validity of the proposed method. Using the ITSP method, only one correction sweep is involved to achieve the velocity distribution, and it is time efficient. Therefore, this nonlinear estimation method may have wide applications in the future.

The accurate velocity distribution is achieved based on the known T-matrix using the ITSP method. 2D numerical examples of the designed models including different perturbation types, demonstrate the validity of the proposed method. The ITSP method considers all scattering effects and overcomes the initial model dependence of Born approximation, besides, it is a regularized method in each correction step and overcomes the divergence effect of Born series for strong velocity perturbation or large perturbation volume. It is also time efficient because only one correction sweep is involved (half from upper part, half from lower part). This method is based on the known T-matrix, suitable for smooth media, then T-matrix estimation method and the improved ITSP method for complex media will be developed in future work.

1引言

隨著勘探開發(fā)的不斷深入，地震勘探逐漸進入巖性油氣藏勘探階段，對地震資料處理以及反演提出了更高的要求.全波形反演可以精確描述地下介質(zhì)參數(shù)空間分布，但是其巨大的計算量以及對初始模型的依賴性，成為其廣泛應用的瓶頸(Virieux and Operto, 2009).為了克服其計算量大的缺陷，GPU加速技術、相位編碼技術以及震源編碼技術被相繼利用，達到提高計算效率的目的(Ben-Hadj-Ali et al., 2009; Larner et al., 1981; Luo et al., 2012; Moghaddam et al., 2013).由于計算量大的限制，大多數(shù)反演方法局限于二維情況，且屬于基于Fréchet導數(shù)的線性或擬線性方法，對初始模型具有一定的依賴性.為了減弱對初始模型的依賴性，許多學者做了大量的研究.Bunks等(1995)應用多尺度反演策略：大尺度反演的結果作為小尺度反演的初值，減弱了對初始模型的依賴性；Ren等(2014)應用第二代小波多尺度方法，實現(xiàn)了品質(zhì)因子以及速度的反演；Shin等(2008;2009)提出了Laplace域以及Laplace-Fourier域全波形反演，其可為波形反演提供初始速度模型；Wu等(2014b)提出了包絡反演方法，在低頻缺失的情況下，可為全波形反演提供可靠的初始速度模型；在包絡反演理論基礎上，Luo和Wu(2015)討論了包絡反演在降低目標函數(shù)局部極小值以及抗噪性方面的特性；另外，Chi等人也同時研究了包絡反演方法以及其為全波形反演方法提供初始模型的應用潛力(Chi et al., 2014).上述方法在減弱初始模型依賴性方面，均是設計對初始模型不敏感的目標函數(shù)，但是在更新迭代時，均基于Born近似，采用線性或擬線性方法，原因之一是其計算效率較高，因為基于蒙特卡洛的完全非線性方法，計算量大，目前無法在大規(guī)模計算中使用.

但是線性方法有其局限性，因此Wu和Zheng(2014)從理論上分析了Born近似(線性Fréchet導數(shù))的弊端，證明很多情況下高階Fréchet導數(shù)是不可以被忽略的，因此引入了非線性偏導數(shù)的概念.考慮到波傳播以及反演過程的非線性性質(zhì)，Weglein等(2003)研究了逆散射級數(shù)方法以及在地球物理勘探中的應用；但是Born級數(shù)以及逆Born級數(shù)的收斂性在很大程度上不能保證，因此為了改善收斂性，眾多學者做了大量研究.Jakobsen(2012)利用量子物理中重整化方法得到改善的Born級數(shù)，改善了原Born級數(shù)的收斂性；Kouri和Vijay(2003)以及Yao等(2014)利用重整化技術將Lippmann-Schwinger方程從Fredholm積分轉化為在1D介質(zhì)復平面具有絕對一致收斂性的Volterra級數(shù)；Kouri等(2004)基于Volterra逆散射級數(shù)，提出了若干策略以得到更好的聲波散射預測；Lesage等(2013)利用反射以及透射數(shù)據(jù)，用重整化的Volterra積分來描述逆聲散射級數(shù).該非線性模擬以及非線性反演，在重整化的框架下，改善了散射級數(shù)以及逆散射級數(shù)的收斂性.

本文基于Wu 和Zheng(2014)建議的非線性反演方法以及Wu等(2014a,2015)提出的對逆散射級數(shù)進行重整化的逆薄板算子，對逆?zhèn)鞑ニ阕永碚摷捌淇焖偎惴ㄗ鲞M一步研究.在二維常密度聲波正演方程中，首先引入傳輸矩陣(T-matrix)，得出擾動數(shù)據(jù)與T-matrix的線性關系；其次，根據(jù)解析關系，得到T-matrix與速度擾動的非線性關系，這種非線性關系可用Born級數(shù)表示；最后，在T-matrix已知的前提下，通過逆散射級數(shù)對速度擾動進行估計，為了保證收斂性，應用De Wolf近似，推導出了逆薄板傳播算子(Inverse Thin-Slab Propagator, ITSP)，以得到地下參數(shù)準確而高效的非線性估計.我們已對高斯球高速異常體模型進行了收斂性分析(Wu et al., 2014a, 2015)，本文對二維常密度聲波方程利用ISTP方法對橢球高速異常體、高斯球低速異常體以及既含高速異常體又含低速體異常體的雙擾動橢球模型進行測試，驗證本文方法的有效性.

2方法原理

2.1Born級數(shù)以及Born近似

對于二維常密度聲波方程，如公式(1)所示：

(1)

其中p為壓力場，v為地下速度分布，s為震源項.方程的解可由Lippmann-Schwinger積分方程進行表征：p(x)=p0(x)+k2∫Vd3x′g0(x;x′)εv(x′)p(x′),

(2)

圖1　波傳播示意圖Fig.1　Cartoon of seismic wave propagation

公式(2)可以寫成矩陣、向量形式：

P=P0+G0VP,

(3)

其中P為總波場(從震源xs出發(fā)的向量)，P0為背景波場(向量)，G0為包含波數(shù)以及格林函數(shù)信息的矩陣，V為實值對角矩陣，包含速度擾動的信息.則利用循環(huán)迭代方式得總波場P的遞推公式：

(4)

將總波場P寫成級數(shù)的形式：

P=P0+G0VP0+G0VG0VP0+…,

(5)

公式(5)對應的級數(shù)即為Born級數(shù).右端第一項為入射場，第二項為一階散射，第三項為二階散射，后面的為高階散射.將所有高次項疊加便可以得到包含各階散射的合成地震記錄，但Born級數(shù)不能保證收斂性，在擾動強度較大或擾動區(qū)域較大時，得到的Born級數(shù)發(fā)散，不能得到期望的合成地震記錄.在擾動較弱時，通常對其進行線性近似，忽略高階項的影響，取其前兩項得Born近似，如公式(6)所示：

P=P0+G0VP0.

(6)

在擾動強度以及擾動范圍較小時，其具有較高的精度；擾動強度以及擾動范圍較大時，其精度較低，與有限差分得到的離散解相差較大.目前大多數(shù)反演方法都是基于Born近似理論，因此對初始模型具有一定的依賴性，只有在擾動相對較小時，才能得到合理的解估計.

2.2T-matrix方法理論

為了提高Born級數(shù)的收斂性，引入傳輸矩陣(T-matrix)的概念，使得：

TP0=VP,

(7)

其中T為復值傳輸矩陣，包含總波場P中各階散射的信息.由于傳輸矩陣T的引入，使得總波場P中的各階散射信息被轉移到傳輸矩陣T之中，基于傳輸矩陣T對速度擾動進行非線性估計.由于各階散射的影響，傳輸矩陣T的結構與速度擾動V不同，其非對角線元素表征了不同擾動點之間的散射信息.將公式(7)反代入公式(3)中得：

P=P0+G0VP=P0+G0TP0,

(8)

則總波場P與速度擾動V之間的非線性關系轉變?yōu)榕c傳輸矩陣T之間的線性關系，T的估計問題可以通過線性優(yōu)化方法進行求解；在正演模擬時，與Born級數(shù)方法相比，公式(8)具有較好的收斂性.傳輸矩陣T與速度擾動V之間的關系為

TP0=VP=V(P0+G0TP0)=VP0+VG0TP0

=(V+VG0T)P0,

(9)

由于P0是從震源xs出發(fā)到達散射點x′的向量，由于震源xs的任意性，使得向量P0是任意的，且滿足公式(9)，公式(9)兩端同時約去P0，得傳輸矩陣T與速度擾動V之間的非線性關系，隱式公式如(10)所示：

T=V+VG0T,

(10)

傳輸矩陣T與速度擾動V的相互關系可顯示表達為

(11)

其中H為傳播算子，將速度擾動V轉化為包含各階散射的傳輸矩陣T；H-為逆?zhèn)鞑ニ阕?，將包含各階散射的傳輸矩陣T恢復為速度擾動V.正演計算時，已知速度擾動V，利用公式(11)計算傳輸矩陣T，再基于公式(8)合成地震記錄；在反演計算時，可借助線性反演理論，基于公式(8)，估計復值傳輸矩陣T，再利用公式(11)進行速度擾動估計.但是當擾動區(qū)域較大時，利用矩陣求逆方法估計速度擾動計算量較大，與矩陣規(guī)模的三次方成正比，難以滿足大規(guī)模計算的要求(Wang et al., 2014).

本文假設傳輸矩陣T已知，估計速度擾動V.由于矩陣求逆的計算量較大，當‖G0T‖<1時，公式(11)中的矩陣求逆可通過級數(shù)的形式進行計算：

(12)

公式(12)實際上為逆Born級數(shù)，可以通過迭代的方式進行計算，遞推公式為

V0=T,

(13)

其中Vn是速度擾動第n次迭代的估計值.在收斂的前提下，隨著迭代的進行，速度擾動的估計值Vn逐漸退化為實值對角矩陣，因此，我們主要關注估計值的對角元素，即

(14)

當速度擾動強度較大或擾動范圍較大時，逆Born級數(shù)(公式(12)—(14))并不能保證收斂性，因此引入DeWolf近似方法(DeWolf, 1971;WuandHuang, 1995;WuandZheng, 2014;Wuetal., 2007)，進行收斂性的改善.

2.3De Wolf 近似

公式(11)在‖G0T‖<1的前提下，可以展開為逆Born級數(shù)；但是在較強擾動或擾動范圍較大時，‖G0T‖<1假設不成立，導致逆Born級數(shù)式(12)以及迭代公式(13)的收斂性不能得到保證，這是逆Born級數(shù)的主要弊端.Wu和Zheng(2012;2014)指出可以將散射算子分解為前向散射以及后向散射算子，再對前向散射算子進行正則化處理，可以將關于傳輸矩陣T的逆Born級數(shù)轉變?yōu)榫哂惺諗啃缘腄e Wolf級數(shù).因此，將傳輸矩陣T分解為前向散射(透射波)以及所有的后向散射(包含反射波在內(nèi)的所有多次反射波)：

(15)

其中Tf是對應于前向散射的傳輸矩陣，Tb是對應于所有后向散射的傳輸矩陣.為簡化問題，僅僅考慮只產(chǎn)生前向散射的光滑介質(zhì)，應用于強反差復雜介質(zhì)的T-matrix方法將是下一步的研究目標.則

(16)

前向散射傳輸矩陣Tf按擾動點(圖2中紅色十字星)所在層的位置可以分成三部分：上半空間速度擾動產(chǎn)生的干涉、下半空間速度擾動產(chǎn)生的干涉以及當前層速度擾動產(chǎn)生的干涉，如圖2所示，則

Tf(x,x′)=Tf-(x,x′)+Tf+(x,x′)+Tf0(x,x′)

=Tu(x,x′)+Td(x,x′)+Tz(x,x′),

(17)

然后根據(jù)Tu, Tz和Td對逆?zhèn)鞑ニ阕親-進行分解得

(18)

其中I為單位算子；Cu,Cd,Cz分別為對應于Tu,Td和Tz的前向散射校正算子，且有

(19)

由于觀測孔徑的影響，Tz一般比較小，甚至不可能得到，因此在數(shù)值試驗時忽略當前層的影響.Cu的第n次校正項為

=(-1)n(GdTu)(xn,xn-1)(GdTu)(xn-1,xn-2)…

(GdTu)(x1,x0),z=zn>zn-1>…>z1>z0

(20)

則公式(18)—(19)可以總結為

(21)

公式(18)轉化為

(22)

2.4逆薄板傳播算子(Inverse Thin-Slab Propagator，ITSP)

對公式(18)進行直接計算，計算量較大.基于重整化的薄板法(Thin-Slab Propagator, TSP)思想(Wu and Wu, 2006)，對于逆散射級數(shù)式(18)，我們推導出了相應的薄板法公式，稱之為逆薄板傳播算子(Inverse Thin-Slab Propagator, ITSP).由公式(22)知

=

.

(23)

首先分析δVu的估計問題.由公式(23)知

=-T GdTu+(-1)2T GdTuGdTu+…,

(24)

以積分或加和的形式將上式寫成核算子的形式，上部散射對第m個薄板的校正值為

(25)

核函數(shù)表達式為

(26)

利用ITSP方法計算第m層校正值為

(27)

(28)

可以看出，在逆薄板法(ITSP)每一步中，均對之前薄板(第m個薄板之前)的多次逆散射進行了部分疊加，該重整化策略消除了逆散射級數(shù)的發(fā)散性問題，圖2顯示了用于高斯球高速異常體模型測試(Wuetal., 2014a;Wuetal., 2015)的ITSP的示意圖：第n層的紅色十字表示要恢復的速度擾動點，該點對應的T-matrix(主頻20Hz)的實部為該示意圖的背景圖(紅色表示正值，藍色表示負值).為了恢復給定點的速度擾動，需要已知其所對應的T-matrix，其分布于整個速度擾動區(qū)域.

圖2　逆薄板法示意圖Fig.2　Schematic cartoon for ITSP (inverse thin-slab propagator)

與計算上部散射校正值一樣，基于ITSP方法得到下部散射校正值：

(29)

另外，當前層的校正值為

(30)

逆薄板傳播算子(ITSP)是一個消除散射的算子，它逐漸將T-matrix中包含的散射信息消除，逐漸退化為速度擾動對應的對角矩陣V.因此，為了恢復某一層的速度擾動，需要知道T-matrix的各個元素Tu(i,m),Td(i,m)和Tz(m,m)，其估計問題將是下一步的研究內(nèi)容.Tz一般比較小，甚至不可能得到，因此數(shù)值試驗時忽略當前層的影響.

3數(shù)值算例

我們已對高斯球高速異常體模型進行了收斂性測試，該試驗指出了逆Born級數(shù)的缺陷，且驗證了ITSP方法的可行性.為了進一步驗證ITSP方法的有效性，本文設計了橢球高速異常體、高斯球低速異常體以及既包含高速異常體又包含低速異常體的雙擾動橢球模型.一系列數(shù)值試驗分析驗證了ITSP方法的可行性以及有效性.

3.1橢球高速異常體模型

橢球高速異常體模型如圖3所示：橢球的長軸為500 m，短軸為300 m，背景速度為2000 m·s-1，高速異常體位于介質(zhì)中間，模型總規(guī)模201×201，異常體的規(guī)模31×51，空間分辨率為10 m.T-matrix是一個頻率依賴的復值矩陣，其規(guī)模為N×N；由于異常體主要分布于模型中間，僅考慮異常體區(qū)域即可，N=1581.反演中利用的T-matrix為頻率20 Hz的精確T-matrix信息，其對應全孔徑測量時的情形，包含全孔徑測量數(shù)據(jù)的全部信息.一般情況下，T-matrix可從有限孔徑的觀測數(shù)據(jù)中，利用線性反演方法得到，至于采集孔徑對T-matrix估計精度的影響將在后續(xù)的研究中進行討論.

圖3　光滑橢球擾動模型(50%速度擾動)Fig.3　The model of smooth ellipse ball with 50% perturbation

圖4　不同擾動強度橢球高速異常體速度估計結果以及相對誤差(a) 15%； (b) 20%； (c) 50%.Fig.4　The estimated velocity distributions and their relative errors for the high velocity anomalous ellipse ball model with different perturbation strengths

利用ITSP方法對速度模型的整個橢球擾動區(qū)域進行處理，估計得到速度分布以及它們的相對誤差，如圖4(a—c)所示，分別為15%、20%和50%速度擾動時估計得到的速度分布以及相對誤差，其中相對誤差定義為絕對誤差除以最大擾動量.可以看出，ITSP方法克服了逆Born級數(shù)的發(fā)散性問題，且精度較高，只經(jīng)過一次掃描校正(一半來自上部，一半來自下部，示意圖見圖2)，效率高.與高斯球高速異常體模型試驗結果(Wu et al., 2014a; Wu et al., 2015)類似，在精確T-matrix已知的前提下，進行速度擾動估計，克服了逆Born級數(shù)的發(fā)散性問題，在全波形反演中具有較好的應用前景.

3.2高斯球低速異常體模型

高斯球低速異常體模型如圖5所示：高斯球的半徑為500 m，背景速度為2000 m·s-1，低速異常體位于介質(zhì)中間，模型總規(guī)模201×201，異常體的規(guī)模51×51，空間分辨率為10 m.T-matrix是一個頻率依賴的復值矩陣，其規(guī)模為N×N；由于異常體主要分布于模型中間，僅考慮異常體區(qū)域即可，N=2601.與橢球高速異常體類似，反演中利用的T-matrix為頻率20 Hz的精確T-matrix信息，其誤差對后續(xù)速度擾動估計的影響將在后續(xù)的研究中進行討論.

圖5　高斯球低速異常體擾動模型(50%速度擾動)Fig.5　The model of the low velocity anomalous Gaussian ball with 50% perturbation

與橢球高速異常體模型類似，利用ITSP方法對速度模型的整個擾動區(qū)域進行處理，得到速度分布以及它們的相對誤差.如圖6(a—c)所示，分別為15%、20%和50%速度擾動時估計得到的速度分布以及相對誤差，其中相對誤差依然定義為絕對誤差除以最大擾動量.反演結果表明ITSP方法精度較高，且只經(jīng)過一次掃描校正，效率高.與高速異常體試驗結果相比，反演的誤差相對較大，特別是擾動較大時(圖6c)，但相對誤差依然在允許范圍之內(nèi).

3.3雙擾動橢球模型

以上模型試驗均為單一類型速度異常體模型(高速異常體或低速異常體)，本小節(jié)設計了既含有高速異常體，又含有低速異常體的雙擾動橢球模型，如圖7所示，每個橢球的長軸為500 m，短軸為300 m，背景速度為2000 m·s-1，雙擾動異常體位于介質(zhì)中間，模型總規(guī)模201×201，異常體規(guī)模為31×101，空間分辨率為10 m.T-matrix依然是一個頻率依賴的復值矩陣，大小為N×N；由于異常體主要分布于模型中間，僅考慮異常體區(qū)域N=3131.類似于以上模型測試，反演中利用精確T-matrix，T-matrix的估計誤差對速度估計的影響，將作為后續(xù)的研究課題.

利用ITSP方法估計得到的速度分布及其對應的相對誤差如圖8(a—c)所示，分別為速度擾動15%、20%以及50%對應的結果.與之前測試結果吻合，ISTP方法克服了逆Born級數(shù)的發(fā)散性問題，估計得到的解精度較高，且具有較高的計算效率.該測試基于精確T-matrix信息，ITSP方法的精度以及計算效率高，應用前景較好.

圖6　不同擾動強度高斯球低速異常體速度估計結果及相對誤差(a) 15%； (b) 20%； (c) 50%Fig.6　The estimated velocity distributions and their relative errors for the low velocity anomalous Gaussian ball model with different perturbation strengths

圖7　雙擾動橢球模型(15%速度擾動)Fig.7　The model of the high and low velocity anomalous ellipse ball with 15% perturbation

圖8　不同擾動強度雙擾動橢球異常體速度估計結果以及相對誤差(a) 15%； (b) 20%； (c) 50%.Fig.8　The estimated velocity distributions and their relative errors for the high and low velocity anomalous ellipse ball model with different perturbation strengths

4結論

本文利用二維常密度聲波方程，討論了基于T-matrix的非線性速度擾動估計方法，在De Wolf近似的前提下，深入研究了逆薄板法(Inverse Thin-Slab Propagator, ITSP)，該算子擺脫了Born近似的束縛；由于考慮了多次散射項，克服了對初始模型的依賴性.在T-matrix已知的前提下，ITSP方法克服了Born級數(shù)的發(fā)散性問題，可對強速度擾動進行高精度的估計；由于其只經(jīng)過一次掃描校正(一半來自上部散射的貢獻、一半來自下部散射的貢獻)，計算效率較高.不同擾動類型以及擾動強度的模擬數(shù)據(jù)分析驗證了本文方法精度高、效率高的優(yōu)點，但其依賴于T-matrix的估計精度，具體量化分析將在后續(xù)研究中進行討論.本文方法基于二維常密度聲波方程，但其可以推廣至三維常密度聲波方程中，與二維方法類似，三維方法也是逐層校正；不同的是，三維中的層面為平面，二維中的層面為線條.

References

Ben-Hadj-Ali H, Operto S, Virieux J. 2009. Three-dimensional frequency-domain full waveform inversion with phase encoding.∥ 2009 SEG Annual Meeting. Houston, Texas: SEG, 2288-2291.

Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion.Geophysics, 60(5): 1457-1473.

Chi B X, Dong L G, Liu Y Z. 2014. Full waveform inversion method using envelope objective function without low frequency data.JournalofAppliedGeophysics, 109: 36-46.

De Wolf D. 1971. Electromagnetic reflection from an extended turbulent medium: cumulative forward-scatter single-backscatter approximation.IEEETransactionsonAntennasandPropagation, 19(2): 254-262. Jakobsen M. 2012. T-matrix approach to seismic forward modelling in the acoustic approximation.StudiaGeophysicaetGeodaetica, 56(1): 1-20. Kouri D J, Vijay A. 2003. Inverse scattering theory: Renormalization of the Lippmann-Schwinger equation for acoustic scattering in one dimension.PhysicalReviewE, 67(4): 046614. Kouri D J, Vijay A, Hoffman D K. 2004. Inverse scattering theory: Strategies based on the volterra inverse series for acoustic scattering.TheJournalofPhysicalChemistryB, 108(29): 10522-10528.

Larner K, Gibson B, Rothman D. 1981. Trace interpolation and the design of seismic surveys.Geophysics, 46(4): 407-415.

Lesage A C, Yao J, Eftekhar R, et al. 2013. Inverse acoustic scattering series using the volterra renormalization of the lippmann-schwinger equation.∥ 2013 SEG Annual Meeting. Houston, Texas: SEG, 4645-4649.

Luo J R, Gao J H, Wang B L. 2012. Multi-GPU based two-level acceleration of full waveform inversion.JournalofSeismicExploration, 21(4): 377-394.

Luo J R, Wu R S. 2015. Seismic envelope inversion: reduction of local minima and noise resistance.GeophysicalProspecting, 63(3): 597-614.

Moghaddam P P, Keers H, Herrmann F J, et al. 2013. A new optimization approach for source-encoding full-waveform inversion.Geophysics, 78(3): R125-R132. Ren Z M, Liu Y, Zhang Q S. 2014. Multiscale viscoacoustic waveform inversion with the second generation wavelet transform and adaptive time-space domain finite-difference method.GeophysicalJournalInternational, 197(2): 948-974.

Shin C, Cha Y H. 2008. Waveform inversion in the Laplace domain.GeophysicalJournalInternational, 173(3): 922-931.

Shin C, Cha Y H. 2009. Waveform inversion in the Laplace-Fourier domain.GeophysicalJournalInternational, 177(3): 1067-1079.

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

Wang B F, Chen X H, Li J Y, et al. 2014. A stable and efficient

attenuation compensation method based on inversion.ChineseJournalofGeophysics(in Chinese), 57(4): 1265-1274, doi: 10.6038/cjg20140423.

Weglein A B, Araújo F V, Carvalho P M, et al. 2003. Inverse scattering series and seismic exploration.InverseProblems, 19(6): R27. Wu R S, Huang L J. 1995. Reflected wave modeling in heterogeneous acoustic media using the De Wolf approximation.∥ Proc. SPIE 2571, Mathematical Methods in Geophysical Imaging III. SPIE, 176-186.Wu R S, Xie X B, Wu X Y. 2007. One-way and one-return approximations (de Wolf approximation) for fast elastic wave modeling in complex media.AdvancesinGeophysics, 48: 265-322.

Wu R S, Zheng Y C. 2012. Nonlinear Fréchet derivative and its De Wolf approximation.∥ 2012 SEG Annual Meeting. Las Vegas, Nevada: SEG, 1-6.

Wu R S, Zheng Y C. 2014. Non-linear partial derivative and its De Wolf approximation for non-linear seismic inversion.GeophysicalJournalInternational, 196(3): 1827-1843.

Wu R S, Hu C H, Wang B F. 2014a. Nonlinear sensitivity operator and inverse thin-slab propagator for tomographic waveform inversion.∥ 2014 SEG Annual Meeting. Denver, Colorado, USA: SEG, 928-933.

Wu R S, Luo J R, Wu B Y. 2014b. Seismic envelope inversion and modulation signal model.Geophysics, 79(3): WA13-WA24.

Wu R S, Wang B F, Hu C H. 2015. Renormalized nonlinear sensitivity kernel and inverse thin-slab propagator in T-matrix formalism for wave-equation tomography.InverseProblems, 31(11): 115004.

Wu X Y, Wu R S. 2006. Wavefield and AVO modeling using elastic thin-slab method.Geophysics, 71(5): C57-C67.

Yao J, Lesage A C, Bodmann B G, et al. 2014. One dimensional acoustic direct nonlinear inversion using the Volterra inverse scattering series.InverseProblems, 30(7): 075006.

附中文參考文獻

王本鋒, 陳小宏, 李景葉等. 2014. 基于反演的穩(wěn)定高效衰減補償方法. 地球物理學報, 57(4): 1265-1274, doi: 10.6038/cjg20140423.

(本文編輯胡素芳)

基金項目國家自然科學基金項目(U1262207)，中國博士后科學基金(2016M590102)和國家科技重大專項課題(2016ZX05024001-005)聯(lián)合資助.

作者簡介王本鋒，男，1986年生，博士后，主要從事地震資料處理以及反演方面的研究工作. E-mail:wbf1232007@126.com

doi:10.6038/cjg20160628 中圖分類號P631

收稿日期2015-07-18，2016-04-19收修定稿

Non-linear parameter estimation method based on T-matrix

WANG Ben-Feng1,2,3, WU Ru-Shan3, CHEN Xiao-Hong2, LU Wen-Kai1

1StateKeyLaboratoryofIntelligentTechnologyandSystems,DepartmentofAutomation,TsinghuaUniversity,Beijing100084,China2StateKeyLaboratoryofPetroleumResources&Prospecting,ChinaUniversityofPetroleum,Beijing102249,China3ModelingandImagingLaboratory,UniversityofCalifornia,SantaCruz, 95064,U.S.A.

AbstractWith the development of seismic exploration and exploitation, it is necessary to study the accurate parameter distribution finely describing the subsurface. Full Waveform Inversion (FWI) can help us achieve this goal, however, traditional FWI theory is based on the Born approximation and the performance is dependent on the initial model. This is one of the main factors which prevent FWI being used widely in practice. In order to study the accurate parameter distribution, we studied nonlinear estimation method based on T-matrix which overcomes the initial model dependence.

KeywordsFull waveform inversion; De Wolf approximation; Transmission matrix (T-matrix); Inverse Thin-Slab Propagator (ITSP); Born series

王本鋒, 吳如山, 陳小宏等. 2016. 基于T-matrix的非線性參數(shù)估計方法.地球物理學報，59(6):2257-2265,doi:10.6038/cjg20160628.

Wang B F, Wu R S, Chen X H, et al. 2016. Non-linear parameter estimation method based on T-matrix.ChineseJ.Geophys. (in Chinese),59(6):2257-2265,doi:10.6038/cjg20160628.