基于梯形網(wǎng)格辛Runge-Kutta射線追蹤的角度域?qū)游銎扑俣确治?/h1>

2015-04-08 05:10:58魏亦文熊彬唐國彬

桂林理工大學(xué)學(xué)報訂閱 2015年3期 收藏

關(guān)鍵詞：分析模型

魏亦文，熊彬，唐國彬

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

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

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

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

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

1 角度域?qū)游銎扑俣确治?/h2>

1.1 角度域?qū)游龇囱莘匠?/h3>

角度域?qū)游龇囱菔窃贏DCIGs基礎(chǔ)上構(gòu)建旅行時反演方程

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

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

1.2 DSRWDC提取ADCIGs

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

1.3 初始速度模型梯形網(wǎng)格參數(shù)化

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

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

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

其中，系數(shù)c1、c2、c3、c4、c5、c6、c7表達(dá)式為

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

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

其中系數(shù)ej1、ej2、ej3、ej4、ej5的取值為

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

表1 偏導(dǎo)數(shù)中的參數(shù)值Table 1 Parameter values of partial derivative

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

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

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

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

式(7)為Hamilton系統(tǒng)，Hamilton函數(shù)為

考慮Hamilton系統(tǒng)

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

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

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

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

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

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

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

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

圖4 TMVA流程圖Fig.4 TMVA flowchart

2.1 初始速度模型

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

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

圖5 矩形網(wǎng)格初始速度模型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初始偏移剖面拾取層位，并根據(jù)圖5c速度模型建立梯形網(wǎng)格速度模型(圖9a)，其均勻采樣輸出速度模型見圖9b，表明速度模型符合地質(zhì)規(guī)律，然后選擇成像點進(jìn)行射線追蹤(圖10)。

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

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

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

圖11 4次迭代結(jié)果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 結(jié)論

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

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

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

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

[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，F(xiàn)omel S.Angle-Gathers by Fourier Transform[R].Stanford Exploration Project，Report 103，2000:123－133.

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

[10]張敏，李振春，張凱，等.疊前偏移共成像點道集的對比分析[J].地球物理學(xué)進(jìn)展，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].地球物理學(xué)報，2006，49(2):540－545.

[14]魏亦文，王彥春.混合插值法重構(gòu)近地表模型[J].計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報，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]劉守偉，程玖兵，王華忠，等.偏移距域/角度域共成像點道集與偏移速度的關(guān)系[J].地球科學(xué)——中國地質(zhì)大學(xué)學(xué)報，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]高亮，李幼銘，陳旭榮，等.地震射線幸?guī)缀嗡惴ǔ跆絒J].地球物理學(xué)報，2000，43(3):402－410.

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

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

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

猜你喜歡

分析模型

一半模型

童話王國·奇妙邏輯推理(2024年5期)2024-06-19 16:03:38

隱蔽失效適航要求符合性驗證分析

民用飛機(jī)設(shè)計與研究(2020年4期)2021-01-21 09:15:02

重要模型『一線三等角』

中學(xué)生數(shù)理化·七年級數(shù)學(xué)人教版(2020年10期)2020-11-26 08:24:50

重尾非線性自回歸模型自加權(quán)M-估計的漸近分布

數(shù)學(xué)物理學(xué)報(2020年2期)2020-06-02 11:29:24

電力系統(tǒng)不平衡分析

電子制作(2018年18期)2018-11-14 01:48:24

電力系統(tǒng)及其自動化發(fā)展趨勢分析

山東工業(yè)技術(shù)(2016年15期)2016-12-01 05:31:22

3D打印中的模型分割與打包

光學(xué)精密工程(2016年6期)2016-11-07 09:07:19

FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究

核科學(xué)與工程(2015年4期)2015-09-26 11:59:03

中西醫(yī)結(jié)合治療抑郁癥100例分析

中國中醫(yī)藥現(xiàn)代遠(yuǎn)程教育(2014年11期)2014-08-08 13:23:44

在線教育與MOOC的比較分析

終身教育研究(2014年5期)2014-02-28 01:23:06

桂林理工大學(xué)學(xué)報2015年3期

桂林理工大學(xué)學(xué)報的其它文章

基于改進(jìn)熵值賦權(quán)法和TOPSIS模型的綜合評價應(yīng)用

基于振動試驗要求的汽車組合開關(guān)夾具設(shè)計與模態(tài)分析

響應(yīng)面法優(yōu)化柚皮基活性炭對CR(Ⅵ)的吸附條件

基于遺傳算法的楊赤中插值權(quán)系數(shù)解算

基于IPA方法的民族村寨旅游開發(fā)社區(qū)居民滿意度分析——以三江程陽橋景區(qū)為例

麥芽糖－三聚氰胺－甲醛共縮聚樹脂的熱性能分析