李勇德， 董良國， 劉玉柱

同濟大學海洋地質(zhì)國家重點實驗室， 上海　200092

?

基于改進的散射積分算法的初至波走時層析

李勇德， 董良國， 劉玉柱*

同濟大學海洋地質(zhì)國家重點實驗室， 上海200092

初至波走時層析是獲取近地表速度結構的一種常用方法.隨著采集技術的不斷發(fā)展，可使用的數(shù)據(jù)量迅速增多，傳統(tǒng)的基于射線追蹤和解方程組的地震走時層析成像方法面臨著內(nèi)存占用大、方程求解不穩(wěn)定等問題.為了解決這些問題，本文基于前人在波形反演研究中提出的一種改進的散射積分算法，提出了一種預條件最速下降法初至波走時層析.該方法無需存儲核函數(shù)矩陣與Hessian矩陣即可方便地實現(xiàn)目標函數(shù)梯度的計算與預條件，且該方法計算效率高、求解穩(wěn)定、易于并行.數(shù)值實驗結果表明，該方法可以獲得與傳統(tǒng)方法精度相當?shù)姆囱萁Y果，但所占用的內(nèi)存大幅減小.

地震層析成像； 散射積分法； 初至波走時； 近地表速度; 海森矩陣； 預條件； 最速下降法

1　引言

地震走時層析的目的是找到一個使理論合成走時與觀測走時之差在某種意義下達到極小的“最佳”模型(Woodward et al., 2008; Zhu et al., 1992).目前基于線性反演思想的方法主要有兩種實現(xiàn)方式.第一種是傳統(tǒng)的求解層析方程組的方法.這種方法求解的方程組系數(shù)矩陣規(guī)模龐大，且非常稀疏，病態(tài)嚴重，需要采用正則化保證反演的正確性與穩(wěn)定性(崔巖和王彥飛，2015)，因此反演精度嚴重依賴于方程組的求解方法.第二種是基于局部最優(yōu)化思想，通過對目標函數(shù)的局部線性化近似獲得速度的修改方向，然后通過迭代實現(xiàn)該非線性反演問題的線性化求解.但無論哪種實現(xiàn)方式，線性反演方法的主要缺點是對初始模型依賴性強，當初始模型選擇不合適時，很容易過早陷入局部極值.盡管如此，這些方法在實際應用中仍可以取得不錯的效果.

在第二種方法中，每一次迭代中速度修改方向的計算是關鍵.利用伴隨狀態(tài)法(Li et al., 2013; 謝春等, 2014)或散射積分法(Chen et al., 2007)可以獲得速度的修改方向.前者利用地震波的一次正傳和數(shù)據(jù)殘差的一次反傳即可獲得目標函數(shù)的梯度，無需存儲大型矩陣，但其預條件在射線走時層析成像中不易實現(xiàn).后者通過顯式地計算核函數(shù)，并與走時殘差向量相乘實現(xiàn)梯度的計算.但與傳統(tǒng)的層析成像方法一樣，散射積分法面臨內(nèi)存占用大的問題，特別是在利用Hessian矩陣時問題更加突出.

Hessian矩陣的利用在局部最優(yōu)化反演中對提高反演精度與多參數(shù)解耦有非常重要的作用.核函數(shù)表達的Hessian矩陣包含兩項，利用準確的Hessian矩陣構造的方向稱為牛頓方向(Pratt et al., 1998)，只保留Hessian第一項的矩陣稱為近似Hessian矩陣，相對應的方向為高斯-牛頓方向(Ma and Hale, 2012).實際應用中，無論是牛頓方向還是高斯-牛頓方向都涉及到Hessian矩陣求逆的操作，因此往往難以準確構造(Métivier et al., 2013).更通常的做法是利用近似Hessian矩陣的對角陣來構造Hessian矩陣的逆，所對應的方向是擬牛頓方向的一種(Shin et al., 2001).由于Hessian矩陣的逆是直接作用在梯度上面，因此這種做法又稱為梯度方向的預條件.正文中會展示，預條件的實施與否對層析有很大的影響.

Liu等(2015)在波形反演研究中提出了一種改進的散射積分算法.該方法將大規(guī)模的核函數(shù)-向量乘運算表示為具有明確物理含義的向量-標量乘的累加運算，因此避免了大規(guī)模矩陣的存儲，使得傳統(tǒng)散射積分法波形反演更加實用化.同時，該方法可以方便地實現(xiàn)Hessian矩陣主對角線元素的計算，進而可以方便實現(xiàn)梯度的預條件.本文將這種在波形反演中提出的改進的散射積分算法應用于初至波走時層析成像中，避免了傳統(tǒng)層析成像中的大規(guī)模矩陣的存儲和預條件實現(xiàn)不方便的問題，進而實現(xiàn)了帶預條件的最速下降法地震走時層析成像.

2　方法原理

2.1散射積分法初至波走時層析方法原理

建立如下最小二乘目標函數(shù)：

(1)

其中，t表示在當前模型下計算出的初至波理論走時，t0表示實際觀測到的走時，s表示介質(zhì)的慢度.在高頻射線假設下，地震波走時t可以表示為

(2)

其中，積分上下限表示沿著射線路徑的曲線積分，s(r)表示空間位置r處的慢度，dl表示射線路徑上的微元.假設共有n個初至波走時數(shù)據(jù)，m個空間未知數(shù)，那么公式(2)離散化后得到方程組為

(3)

其中，sj表示模型空間中某一點的慢度，ti表示某一炮檢對所對應的初至波到達時.kij表示第i條射線在第j個模型網(wǎng)格內(nèi)的長度，即系數(shù)矩陣中的每一行表示一條射線路徑.公式(3)用矩陣符號表示為

(4)

這里，K即為走時層析中的敏感核函數(shù)(Fréchet導數(shù)矩陣)，其規(guī)模為n×m.模型參數(shù)的修改量可以表示為修改方向p與步長α的乘積.通過對目標函數(shù)(1)式進行二階泰勒級數(shù)展開，可以求得慢度修改方向p：

(5)

這里，Δt表示實際觀測走時與當前模型下理論合成走時之差，H表示準確的Hessian矩陣，是目標函數(shù)對慢度參數(shù)的二階導數(shù)，對于它的不同近似，可以得到不同的修改方向.步長α可以用下式求得

(6)

根據(jù)先驗信息給定初始模型s0，對于第l次迭代，可利用公式(5)求得模型參數(shù)修改方向，用公式(6)計算出步長，然后利用下式對模型進行更新：

(7)

2.2基于改進的散射積分算法的初至波走時層析

高斯-牛頓方向(用符號Ha表示，H≈Ha=KTK)可以用矩陣表示為

(8)

(8)式中矩陣KT規(guī)模龐大，不易存儲，利用Liu等(2015)的改進算法，(8)式右端的大型核函數(shù)-向量乘可以表示為如下多個向量-標量乘的累加運算：

(9)

根據(jù)前面的分析可知，公式(9)中等式右端每一個列向量(ki1…kij…kim)T即為某一炮檢對所對應的射線路徑.因此，不必再存儲完整的矩陣KT，僅需要存儲單一炮檢對所對應的射線路徑即可實現(xiàn)梯度的計算，大大減少了內(nèi)存的占用.

(10)

這樣就避免了存儲大型矩陣KT和K，僅需將單炮檢對核函數(shù)中的每一個元素平方，再對所有的炮檢對累加求和即可.為了避免由于H0對角元素中存在很小的數(shù)而產(chǎn)生計算的不穩(wěn)定，需要引入正則化項，所以預條件后的最速下降方向為

(11)

其中，I表示單位矩陣，λ表示一個較小的數(shù).將式(11)代入式(7)，即可實現(xiàn)預條件最速下降法初至波走時層析.

3　數(shù)值實驗

為了驗證本文方法的有效性，選取了一個二維簡單模型和一個二維復雜起伏地表模型進行理論模型測試.觀測方式分別為全方位觀測和地表觀測.基于相同的走時正演方法(劉洪等，1995)，同一組“觀測數(shù)據(jù)”，利用改進的散射積分法與傳統(tǒng)解方程組的方法分別進行初至波走時層析(為方便，后文分別將二者稱為散射積分法層析與傳統(tǒng)層析).需要注意的是，為了提高計算效率與減少內(nèi)存占用量，本文傳統(tǒng)層析中大型病態(tài)方程組的求解采用的是SIRT法，且沒有進行任何正則化處理(劉玉柱和董良國,2007；劉玉柱等，2007).

3.1全方位觀測情況下的簡單模型實驗

二維簡單模型(圖1a)是一個包含101×101個離散網(wǎng)格的正方形模型，空間離散間隔為10 m×10 m.背景速度為1000 m·s-1，中間有一個半徑為100 m的圓形高速異常體，速度為1160 m·s-1.在模型每一邊均勻放置50個炮點，共計200個炮點.檢波器均勻放置在模型四周，每一邊50個，共200個.基于真實模型的高精度射線追蹤得到的走時作為觀測走時，初始模型為均勻背景模型，經(jīng)過20次迭代后，反演結果如圖1所示.

可以看出，散射積分法層析第一次迭代的梯度(圖1b)主要聚焦在速度異常附近，在速度無異常處梯度盡管不為0(這是照明角度不均勻造成的)，但是相對較弱.為了更好地對比兩種反演方法的精度，抽取地下500 m處的速度曲線(圖2).可以看出，散射積分法層析與傳統(tǒng)層析的反演精度基本相當.表1為兩種方法的計算時間與內(nèi)存占用量對比，可以看出，散射積分法層析內(nèi)存占用僅為傳統(tǒng)層析的4.6%，計算效率也比傳統(tǒng)層析高.

表1　二維簡單模型不同方法內(nèi)存占用量與 一次迭代所需時間對比Table 1　Comparisons of the memory requirements and computation time between different methods

圖2　地表下500 m處速度曲線對比初始模型(粉色)、真實模型(黑色)、散射積分法層析反演結果(紅色)、傳統(tǒng)層析反演結果(藍色).Fig.2　Profile comparisons of the results using SI method and traditional method at depth 500 mStarting model (pink); True model (black); Result using SI method (red) and the result using traditional method (blue).

圖1　二維簡單模型反演結果對比圖(a) 真實模型； (b) 散射積分法層析第一次迭代時的梯度； (c) 散射積分法層析20次迭代后的反演結果； (d) 傳統(tǒng)層析20次迭代后的反演結果.Fig.1　Comparisons of the inversion results of a 2D simple model(a) True model; (b) Gradient of the first round computed by SI method; (c) The inversion result using SI method after 20 iterations; (d) The inversion result using the traditional method after 20 iterations.

3.2地表觀測情況下的復雜模型實驗

為了進一步驗證該方法的反演能力，本文針對二維復雜起伏地表模型進行實驗(圖3).該模型地表起伏落差達450 m，近地表速度橫向變化劇烈.模型被離散為4001×151個網(wǎng)格，空間離散間隔為10 m×10 m.利用彈性波方程模擬地表地震記錄，并拾取Z分量記錄初至波走時作為觀測數(shù)據(jù).共模擬了619炮數(shù)據(jù)，第一炮位于5 km處，炮點間隔40 m.檢波器均勻、對稱地分布在炮點兩端，水平間隔為20 m.每炮201個檢波器，最小炮檢距為0 m，最大炮檢距為2.0 km.水平位置9 km和26 km處的兩個Z分量炮記錄如圖4所示.由于初至波走時層析僅能反演得到近地表的速度結構，因此為了更好地對比反演結果，下文僅顯示0～0.7 km深度范圍內(nèi)的速度結構.該范圍的真實模型如圖5a所示，初始模型如圖5b所示，散射積分法層析與傳統(tǒng)層析經(jīng)過40次迭代后的反演結果分別如圖5c、5d所示.為了定量地對比不同方法的反演結果，抽取了地下10 m、100 m、200 m處的速度曲線(圖6).由于模型左右兩端并沒有炮點覆蓋，無法得到正確結果，因此僅顯示了水平方向10～30 km的速度模型.對比不同深度的速度曲線可知，散射積分法層析可以獲得與傳統(tǒng)層析精度相當?shù)姆囱萁Y果.表2為兩種方法的內(nèi)存占用量與計算效率對比.可以看出，散射積分法層析的內(nèi)存占用量僅為傳統(tǒng)層析的8.2%，計算效率也略高.

圖3　真實模型(深度范圍0～1.5 km)Fig.3　True model (0～1.5 km in depth)

圖4　彈性波模擬地表Z分量觀測記錄((a)9 km處；(b)26 km處)Fig.4　Z component records of the elastic wave simulation at (a) 9 km and (b) 26 km

圖5　0.7 km深度范圍內(nèi)的(a)真實模型、(b)初始模型、(c)散射積分法層析40次迭代反演結果與(d)傳統(tǒng)層析40次迭代反演結果Fig.5　Zoomed views of the (a) true model, (b) starting model and inversion results using (c) SI method and (d) traditional method after 40th iterations, from 0 to 0.7 km in depth

圖6　地表下(a)10 m、(b)100 m、(c)200 m深度處的速度曲線對比初始模型(粉色)、真實模型(黑色)、散射積分法層析反演結果(紅色)、傳統(tǒng)層析反演結果(藍色).Fig.6　Comparisons of the horizontal profiles of the inversion results using SI method and traditional method at depths of (a) 10 m, (b) 100 m and (c) 200 mStarting model (pink); True model (black); Results using SI method (red) and traditional method (blue).

圖7　預條件對梯度的影響第一次迭代時的(a)無預條件最速下降方向與(b)預條件最速下降方向.Fig.7　Influence of preconditioning on gradientThe gradients of the first iteration (a) without preconditioning and (b) with preconditioning.

圖8　0.7 km深度范圍內(nèi)的最速下降法74次迭代后的結果Fig.8　Zoomed views of the inversion results using the SI method without preconditioning after 74th iterations, above 0.7 km

內(nèi)存占用量/Mb一次迭代反演所需時間/s傳統(tǒng)層析23339.92散射積分法層析1937.79

圖9　不同方法目標函數(shù)(自然對數(shù)值)變化曲線不實施預條件(藍色)、實施預條件(紅色)情況下的目標函數(shù).Fig.9　Comparisons of the object functionWithout preconditioning (blue) and with preconditioning (red).

4　討論

公式(11)采用預條件是該方法可以取得與傳統(tǒng)方法相當?shù)姆囱菥?圖6)的關鍵.在走時層析中，預條件的實現(xiàn)去除了射線密度的影響，起到了“照明補償”的效果.因此，與無預條件的梯度(圖7a)相比，預條件(圖7b)使得梯度修改范圍更深，速度修改方向更加合理，反演精度更高.如圖8所示，在不實施預條件的情況下，即使迭代次數(shù)達到74次依舊無法得到預條件情況下40次的反演精度.目標函數(shù)(自然對數(shù)值)的收斂情況(圖9)也有力地說明了預條件對于初至波走時層析的重要性.與無預條件的最速下降法相比，實施預條件可以使得目標函數(shù)收斂更快，并且避免過早陷入局部極值.

本文介紹的散射積分法層析對H0的實現(xiàn)有天然的優(yōu)勢(公式10).在計算梯度的同時便可以獲得H0，幾乎不需要增加額外的計算量.該思路可以直接推廣到三維層析中，為大規(guī)模三維近地表速度建模提供了一種簡單易行的方法.

5　結論

本文將改進的散射積分算法應用于地震走時層析成像中，將龐大的核函數(shù)-向量乘運算表示為具有明確物理含義的向量-標量乘的累加運算.該方法同時方便地實現(xiàn)了Hessian對角元素的累加運算，進而實現(xiàn)了預條件的最速下降法初至波走時層析.二維簡單模型和二維復雜起伏地表理論模型測試表明，該方法可以獲得與傳統(tǒng)解方程組的地震層析成像方法精度相當?shù)姆囱萁Y果，但內(nèi)存占用量大大降低.同時，該方法還具有計算效率高、物理意義明確、反演穩(wěn)定、易于并行等優(yōu)點.這為以后大規(guī)模野外觀測數(shù)據(jù)，尤其是三維地震層析成像處理提供了一條新的思路.

Chen P, Jordan T H, Zhao L. 2007. Full three-dimensional tomography: a comparison between the scattering-integral and adjoint-wavefield methods.GeophysicalJournalInternational, 170(1): 175-181. Cui Y, Wang Y F. 2015. Tikhonov regularization and gradient descent algorithms for tomography using first-arrival seismic traveltimes.ChineseJ.Geophys. (in Chinese), 58(4): 1367-1377, doi: 10.6038/cjg20150423.

Li S W, Vladimirsky A, Fomel S. 2013. First-break traveltime tomography with the double-square-root eikonal equation.Geophysics, 78(6): U89-U101.

Liu H, Meng F L, Li Y M. 1995. The interface grid method for seeking global minimum travel-time and the correspondent raypath.ChineseJ.Geophys. (in Chinese), 38(6): 823-832.

Liu Y Z, Dong L G. 2007. Analysis of influence factor of first-breaks tomography.OilGeophysicalProspecting(in Chinese), 42(5): 544-553. Liu Y Z, Dong L G, Xia J J. 2007. Normalized approach in tomographic imaging of first breaks travel time.OilGeophysicalProspecting(in Chinese), 42(6): 682-685, 698.Liu Y Z, Yang J Z, Chi B X, et al. 2015. An improved scattering-integral approach for frequency-domain full waveform inversion.GeophysicalJournalInternational, 202(3): 1827-1842.Ma Y, Hale D. 2012. Quasi-Newton full-waveform inversion with a projected Hessian matrix.Geophysics, 77(5): R207-R216.

Métivier L, Brossier R, Virieux J, et al. 2013. Full waveform inversion and the truncated newton method.SIAMJournalonScientificComputing, 35(2): B401-B437.

Pratt G, Shin C, Hicks. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.GeophysicalJournalInternational, 133(2): 341-362.Shin C, Jang S, Min D J. 2001. Improved amplitude preservation for prestack depth migration by inverse scattering theory.GeophysicalProspecting, 49(5): 592-606.

Woodward M J, Nichols D, Zdraveva O, et al. 2008. A decade of tomography.Geophysics, 73(5): VE5-VE11.

Xie C, Liu Y Z, Dong L G, et al. 2014. First arrival traveltime tomography based on the adjoint state method.OilGeophysicalProspecting(in Chinese), 49(5): 877-883. Zhu X H, Sixta D P, Angstman B G. 1992. Tomostatics: turning-ray tomography+static corrections.TheLeadingEdge, 11(12): 15-23.

附中文參考文獻

崔巖, 王彥飛. 2015. 基于初至波走時層析成像的Tikhonov正則化與梯度優(yōu)化算法. 地球物理學報, 58(4): 1367-1377, doi: 10.6038/cjg20150423.

劉洪, 孟凡林, 李幼銘. 1995. 計算最小走時和射線路徑的界面網(wǎng)全局方法. 地球物理學報, 38(6): 823-832.

劉玉柱, 董良國. 2007. 初至波層析影響因素分析. 石油地球物理勘探, 42(5): 544-553.

劉玉柱, 董良國, 夏建軍. 2007. 初至波走時層析成像中的正則化方法. 石油地球物理勘探, 42(6): 682-685, 698.

謝春, 劉玉柱, 董良國等. 2014. 伴隨狀態(tài)法初至波走時層析. 石油地球物理勘探, 49(5): 877-883.

(本文編輯汪海英)

First arrival traveltime tomography based on an improved scattering-integral algorithm

LI Yong-De, DONG Liang-Guo, LIU Yu-Zhu*

StateKeyLaboratoryofMarineGeology,TongjiUniversity,Shanghai200092,China

First arrival traveltime tomography is a commonly used method at present to obtain the near-surface velocity structure. With the development of acquisition technology, the data quantity is growing rapidly. Thus, traditional ray-based equations-solving method faces the challenges of memory consuming and instability. To solve these problems, an improved scattering-integral (SI) algorithm proposed in a former study on full-waveform inversion (FWI) is employed here to develop a new preconditioned steepest-descent ray-based first arrival traveltime tomography. This method can compute the gradient of the objective function and implement the precondition conveniently without storing the sensitivity kernel matrix or Hessian matrix in advance. In addition, the inversion process is efficient, stable and easy to be paralleled. Numerical experiments show that this method can achieve accurate inversion results comparable to traditional methods, while the memory requirement is significantly reduced.

Seismic tomography；Scattering-integral algorithm；First-arrival traveltime；Near surface velocity; Hessian matrix; Precondition; Steepest-descent method

10.6038/cjg20161026.

國家自然科學基金(41274116,41474034,41674115,41630964),中央高校基本科研基金(20123197,1350219162)共同資助.

李勇德，男，1992年生，碩士在讀，主要從事初至波走時層析方法與應用研究. E-mail: tju_lyd@126.com

劉玉柱，男，1979年生，博士，副教授、博士生導師，主要從事地震波正反演理論、方法與應用研究.E-mail: liuyuzhu@#edu.cn

10.6038/cjg20161026

P631

2015-12-05，2016-08-26收修定稿

李勇德，董良國，劉玉柱. 2016. 基于改進的散射積分算法的初至波走時層析. 地球物理學報，59(10):3820-3828,

Li Y D, Dong L G, Liu Y Z. 2016. First arrival traveltime tomography based on an improved scattering-integral algorithm.ChineseJ.Geophys. (in Chinese),59(10):3820-3828,doi:10.6038/cjg20161026.