基于梯形網格辛Runge-Kutta射線追蹤的角度域層析偏移速度分析

魏亦文，熊彬，唐國彬

(桂林理工大學a.地球科學學院;b.廣西礦冶與環境科學實驗中心，廣西桂林541004)

復雜構造區采集的低信噪比地震數據經過偏移后生成的共成像點道集(CIGs)比原始數據顯得簡單些，可將旅行時反演與偏移速度分析(MVA)組合為層析偏移速度分析方法(TMVA)，將波動方程層析反演和MVA組合成波動方程偏移速度分析法(WEMVA)。研究和實踐表明，在成像空間中，針對速度比不大的地下結構成像和偏移速度建模，TMVA比WEMVA的應用更普遍。在考慮構建TMVA的反演線性方程組問題時，剩余時差和偏導數矩陣的求取是最關鍵的問題。

疊前深度MVA是建立復雜介質速度模型的有效工具。Faye等[1]提出深度聚焦分析(DFA);Al-Yahya等提出剩余曲率分析(RCA)，利用剩余校正(RMO)來確定速度誤差，CIGs平直性來判斷速度模型的正確性[2];Berkhout等提出共聚焦速度分析(CFP)，然后對其進行推廣應用[3];Stork等基于RCA理論將反射旅行時反演應用到MVA中，提出TMVA[4]。

采用Kirchhoff偏移作為引擎的MVA，生成的CRP道集可能存在假象，利用DSRWDC或RTM提取的ADCIGs被證明是沒有假象的道集。Prucha等采用DSRWDC生成ADCIGs[5]，并將其應用于速度分析;Clapp等利用反射層位約束，基于ADCIGs實現了TMVA[6];Biondi等分析了ADCIGs的運動學特征，推導了求取RMO和旅行時殘差的公式[7];Sava等提出了時移成像條件，通過傾斜疊加得到時移角道集并用于偏移速度分析[8];劉守偉等基于時移ADCIGs，提出將DFA和RCA兩種方法進行組合求取偏移速度模型[9];張敏等對各種道集進行對比分析指出ADCIGs最適合MVA[10];秦寧等實現了自動拾取的成像空間走時反演速度模型[11]。研究表明，復雜構造中ADCIGs與偏移速度模型之間的關系比較復雜，確定兩者之間的關系是角度域層析反演的基礎。

本文實現了在角度域更新層速度模型的TMVA方法:①通過常規速度分析獲得矩形網格初始層速度模型，在此基礎上采用DSRWDC提取ADCIGs;②求取ADCIGs的RMO并計算旅行時殘差作為層析反演線性方程組的右端;③根據初始速度模型和偏移層位構建梯形網格速度模型，在此基礎上采用辛Runge-Kutta射線追蹤算法構建方程組左端的偏導數矩陣;④建立TMVA處理流程，對青海某工區實際地震資料進行處理，驗證該方法的可行性和有效性。

1 角度域層析偏移速度分析

1.1 角度域層析反演方程

角度域層析反演是在ADCIGs基礎上構建旅行時反演方程

其中:L是靈敏度矩陣;Δt是走時殘差;Δz是ADCIGs的剩余深度;Δs是慢度更新;α是地層傾角;γ是反射角。Δz可以通過ADCIGs剩余曲率選取控制點自動擬合求解，具體步驟見文獻[11];式(1)中Δt具體推導見文獻[7]。

為高效求解最小二乘法(LSQR)正規方程，數據量不大時直接采用GPU-LU或者遞歸LU[12－14]，數據量很大時采用矩陣壓縮存儲技術和并行算法[13]。

1.2 DSRWDC提取ADCIGs

單程波方程波場向下延拓包括炮域單平方根方程波場向下延拓和CMP域雙平方根方程波場向下延拓(DSRWDC)[8，15－16]，兩種方法提取ADCIGs效果相同。本文選取DSRWDC提取ADCIGs，將波場變換到頻率域沿深度方向向下延拓，對頻率求和相當于取零時間成像條件，得到局部半偏移距共成像點道集(ODCIGs)。速度模型正確時ODCIGs能量應該聚焦在零偏移距處，否則不聚焦。為利于RCA速度分析，利用文獻[16]中的方法將ODCIGs快速轉化為ADCIGs，算法的實現調用了美國Stanford大學勘探項目組(SEP)開發的SEPlib庫。

1.3 初始速度模型梯形網格參數化

常規MVA方法可以獲得旅行時反演需要的初始速度模型，主要包括疊加速度分析、Dix層速度反演、Kirchhoff MVA[17－18]。矩形網格速度模型通過層位控制插值后可以轉換成梯形網格[12]，圖1是Zelt C等提出速度模型的梯形網格參數化方式，每個梯形單元的上下兩條邊具有不同的斜率k1、k2和截距b1、b2，但左右兩條邊必須是垂直的，x1、x2分別表示梯形網格左側和右側邊界橫坐標。可以利用梯形單元的4個節點速度或慢度(速度的倒數)插值得到梯形單元內任何位置的速度或慢度值。本文計算中梯形單元4個節點是慢度s1、s2、s3、s4，梯形單元內的任意點慢度值沿著梯形的4個邊都是線性變化的，因此在梯形單元中存在水平慢度梯度或者垂直慢度梯度。地層界面可以用梯形網格上下界面連接組成，為計算界面的斜率需要對界面節點進行三次樣條插值平滑處理[14]，以滿足射線追蹤需要。梯形網格速度建模利用了構造信息，因此它更符合地質規律。

圖1 梯形網格與矩形網格慢度插值示意圖Fig.1 Slowness interpolation of trapezoid and rectangular grid

圖1中梯形網格任意位置(x，z)處的慢度s(x，z)計算公式為

其中，系數c1、c2、c3、c4、c5、c6、c7表達式為

當s1、s2、s3、s4的4個頂點坐標分別為(0，0)，(1，0)，(0，1)，(1，1)時，利用式(3)求取系數c1～c7并代入式(2)得

式(4)就是矩形網格雙線性插值的基本公式，顯然它是梯形網格雙線性插值的特例。把慢度插值函數看成是4個頂點慢度si的函數，從式(3)可以看出，系數c1～c7中除了c6、c7與s1、s2、s3、s4無關外，其余均有關，由此可以得到梯形網格慢度場對4個頂點慢度的偏導數為

其中系數ej1、ej2、ej3、ej4、ej5的取值為

對第j個頂點來說，sign、k、xb、b是常數，見表1。

表1 偏導數中的參數值Table 1 Parameter values of partial derivative

當s1、s2、s3、s4的4個頂點坐標分別為(0，0)，(1，0)，(0，1)，(1，1)時，根據式(2)和式(3)及表1容易證明梯形網格慢度函數s(x，y)對各節點慢度sj的導數與式(4)慢度函數對各節點慢度導數相同。

1.4 角度域辛Runge-Kutta射線追蹤

Zelt C開發了RAYINVR軟件，針對寬角反射射線追蹤，在梯形網格中使用Runge-Kutta算法求解簡化的射線追蹤方程[12]。本文并不對射線方程進行簡化，根據射線追蹤方程一般形式，引入弧長τ得到射線追蹤方程組[19－20]

其中:x=(x1，x2)為坐標;s(x1，x2)為介質慢度分布;T為射線走時;是慢度向量，τ(T)=τ(T0)+∫Tc2dT。

式(7)為Hamilton系統，Hamilton函數為

考慮Hamilton系統

可以證明式(7)為Hamilton量H(x，p)≡0的正則方程，聯合式(8)有

其中，h為步長;j為迭代次數，j=1，2，3，4。給定當前點(k(0)，x(0))利用式(11)經過4次迭代計算得到下一個點射線參數和位置(k(4)，x(4))。aj，bj是最優化系數參見文獻[22]。

由于式(11)表示的辛Runge-Kutta算法直接在τ域計算射線路徑，不需數值積分求旅行時，因此能高效求取反演線性方程組左側的偏導數矩陣(靈敏度矩陣)。設梯形網格節點個數為N，射線有M條，式(1)中偏導數矩陣L的元素為第i(i=1，2，…，M)條射線旅行時對第j(j=1，2，…，N)個節點慢度sj的偏導數。首先，找到以sj為頂點的所有梯形單元，單元個數為K(1≤K≤4)，求取每個梯形單元內射線路徑，然后利用式(5)求取每個單元慢度函數對頂點sj的梯度作為權系數，最后對K個射線路徑加權平均得到頂點sj對應的偏導數矩陣元素。

角度域射線追蹤是從成像點向地表進行的，對于圖2所示的理論速度模型，采用層位控制梯形網格參數化方式，應用式(11)進行射線追蹤(圖3)。

圖2 理論層速度模型Fig.2 Theoretic interval velocity model

圖3 辛Runge-Kutta射線追蹤Fig.3 Symplectic Runge-Kutta seismic ray tracing

2 青海2-D疊前資料處理

針對青海的某2-D工區低信噪比疊前地震數據，通過前期預處理，再利用疊加速度分析、Dix反演和Kirchhoff偏移速度分析建立初始速度模型，然后應用TMVA求偏移速度模型，本文采用的速度分析處理流程(圖4)。

圖4 TMVA流程圖Fig.4 TMVA flowchart

2.1 初始速度模型

首先采用疊加速度分析求取均方根速度(RMS)，然后通過最小二乘Dix反演層速度模型，再應用常規MVA求偏移速度，最后利用文獻[13]中混合插值法進行圓滑得到初始速度模型(圖5)。

采用Dix反演層速度模型(圖5b)和Kirchhoff疊前深度偏移提取的地面偏移距共成像點道集(CIGs)，CMP范圍215～221的CIGs，其同相軸向下彎曲(圖6a)，表明偏移速度過大[9－11];圖6b為偏移距范圍100～1 000 m的疊加剖面。

圖5 矩形網格初始速度模型Fig.5 Initial velocity model of rectanguler grids

圖6 Kirchhoff疊前深度偏移Fig.6 Kirchhoff pre-stack depth migration

2.2 提取ADCIGs

利用初始速度模型和DSRWDC偏移提取的ADCIGs(－30°～30°)見圖7，CMP位置從722開始，間隔75，圖7同相軸向上彎曲較大，表明ADCIGs中存在明顯的剩余深度。對ADCIGs疊加，見圖8。

圖7 初始速度模型DSRWDC提取的ADCIGsFig.7 ADCIGs extracted by DSRWDC with initial velocity model

圖8 初始ADCIGs疊加剖面Fig.8 Initial stack of ADCIGs

2.3 層析偏移速度分析

利用圖6b初始偏移剖面拾取層位，并根據圖5c速度模型建立梯形網格速度模型(圖9a)，其均勻采樣輸出速度模型見圖9b，表明速度模型符合地質規律，然后選擇成像點進行射線追蹤(圖10)。

利用圖9速度模型和TMVA迭代4次后的速度模型見圖11。用DSRWDC偏移后，與圖7相同CMP位置處的ADCIGs顯示結果見圖12，其ADCIGs小角度(－10°～10°)同相軸比圖7更平直，10°以上同相軸有彎曲，為大角度數據缺失后偏移所致。對ADCIGs進行疊加(圖13)，與圖8角度域初始疊加剖面相比，成像精度明顯得到了提高，說明TMVA迭代更新后速度模型精度也得到了提高。

圖9 初始速度模型Fig.9 Initial velocity model

圖10 角度域辛Runge-Kutta射線追蹤Fig.10 Angle domain symplectic Runge-Kutta seismic ray tracing

圖11 4次迭代結果Fig.11 Result after 4 times iteration

圖12 最終速度模型偏移ADCIGsFig.12 ADCIGs migrated by final velocity model

圖13 TMVA后角度域疊加剖面Fig.13 Angle domain stack after TMVA

3 結論

(1)對速度模型采用了梯形網格參數化方式，與矩形網格相比，可靈活選擇層位控制點及速度控制點位置，且更符合地質規律。

(2)根據ADCIGs視反射角范圍從成像點向地表進行射線追蹤，辛Runge-Kutta算法求解τ域射線追蹤方程，不需要數值積分求射線走時，只需要計算射線路徑可快速求取層析反演線性方程組左側的偏導數矩陣。

(3)在青海地震資料處理中，使用的角度域多次迭代TMVA方法建立的圓滑層速度模型比常規疊前MVA方法建立的速度模型精度更高，表明該方法是可行、有效的。

中國地質大學(北京)王彥春教授提供了地震數據，并在處理方法上給予了很好的建議，在此表示衷心的感謝。

[1]Faye J P，Jeannot J P.Prestack migration velocities from focusing depth analysis[C]//SEG 56th Annual Meeting，Expanded Abstracts，Houston:Society of Exploration Geophysicist，1986:438－440.

[2]Al-Yahya K.Velocity analysis by iterative profile migration[J].Geophysics，1989，54(6):718－729.

[3]Berkhout A J，Rietveld W E.Determination of macro models for prestack migration:Part 1 Estimation of macro velocities[C]//SEG 64th Annual Meeting，Expanded Abstract，Los Angeles:Society of Exploration Geophysicist，1994:1330－1333.

[4]Stork C，Clayton R W.Linear aspects of tomographic velocity analysis[J].Geophysics，1991，56(4):483－495.

[5]Prucha M L，Biondi B L，Symes W W.Angle-domain common-image gathers by wave-equation migration[C]//SEG 69th Annual Meeting，Expanded Abstracts，Houston:Society of Exploration Geophysicist，1999:824－827.

[6]Clapp R G，Biondi B L，Claerbout J F.Incorporating geologic information into reflection tomography[J].Geophysics，2004，69(2):533－546.

[7]Biondi B，Symes W W.Angle-domain common-image gathers for migration velocity analysis by wavefield-continuation imaging[J].Geophysics，2004，69(5):1283－1298.

[8]Sava P，Fomel S.Angle-Gathers by Fourier Transform[R].Stanford Exploration Project，Report 103，2000:123－133.

[9]劉守偉，王華忠，程玖兵，等.時空移動成像條件及偏移速度分析[J].地球物理學報，2008，51(6):1883－1891.

[10]張敏，李振春，張凱，等.疊前偏移共成像點道集的對比分析[J].地球物理學進展，2011，26(1):220－228.

[11]秦寧，李振春，楊曉東，等.自動拾取的成像空間域走時層析速度反演[J].石油地球物理勘探，2012，47(3):392－398.

[12]Zelt C A，Smith R B.Seismic traveltime inversion for 2-D crustal velocity structure[J].Geophysical Journal International，1992，108(1):16－34.

[13]劉勁松，劉福田，劉俊，等.地震層析成像LSQR算法的并行化[J].地球物理學報，2006，49(2):540－545.

[14]魏亦文，王彥春.混合插值法重構近地表模型[J].計算機輔助設計與圖形學學報，2012，24(4):466－470.

[15]Rickett J，Sava P.Offset and angle domain common image gathers for shot profile migration[R].Stanford Exploration Project，Report 108，2001:1－8.

[16]劉守偉，程玖兵，王華忠，等.偏移距域/角度域共成像點道集與偏移速度的關系[J].地球科學——中國地質大學學報，2007，32(4):575－582.

[17]Liu Z Y.An analytical approach to migration velocity analysis[J].Geophysics，1997，62(4):1238－1249.

[18]Bloor R，Gonzalez A，Alertin U，et al.Tomographic velocity model updating for prestack depth migration[C]//SEG 69th Annual Meeting，Expanded Abstracts，Houston:Society of Exploration Geophysicist，1999:1255－1258.

[19]高亮，李幼銘，陳旭榮，等.地震射線幸幾何算法初探[J].地球物理學報，2000，43(3):402－410.

[20]李川，王有學，何曉玲，等.基于二維三次卷積插值算法的辛幾何射線追蹤[J].地球物理學報，2014，57(4):1235－1240.

[21]韓果花，李川，黃明曦，等.二維復雜速度結構中的射線追蹤[J].桂林理工大學學報，2013，33(3):425－429.

[22]McLachlan R I，Atela P.The accuracy of symplectic integrators[J].Nonlinearity，1992，5(2):541－562.