用偽譜法合成零偏移距地震記錄

1、第28卷第4期長春科技大學學報V o l.28N o.41998年10月JOU RNAL O F CHAN GCHUN UN I V ER S IT Y O F SC IEN CE AND T ECHNOLO GYO ct.1998用偽譜法合成零偏移距地震記錄張文生1何樵登2李樹勇1(1中國科學院計算數(shù)學與科學工程計算研究所,北京100080;2長春科技大學,長春130026摘要零偏移距地震記錄是疊后偏移的基礎,常用射線追蹤法來制作,但當介質復雜時,射線會出現(xiàn)交叉現(xiàn)象,對三維情況計算量更大。從二維標量波動方程出發(fā),用偽譜法研究了零偏移距地震記錄的合成,給出了不同情況下的數(shù)值計算公式,并對邊界吸

2、收問題作了有效處理。數(shù)值計算結果表明了方法的正確性和有效性。該方法可用于偏移中正演模擬的快速計算,且易推廣到三維情形。關鍵詞偽譜法正演零偏移距記錄反周期擴展法中圖分類號P631.40引言在波動方程偏移中常需制作零偏移距地震記錄,零偏移距記錄通常用射線追蹤法和波動方程法來制作。射線追蹤法依據(jù)的是幾何地震學,不易反映波場的動力學特點,另外對復雜介質會出現(xiàn)射線交叉現(xiàn)象,對三維情況,因需作走時計算,計算量更大。波動方程法基于波在介質中傳播的波動性,通過數(shù)值求解波動方程來實現(xiàn),能較全面地反映波場特征。波動方程數(shù)值求解通常有有限差分法、有限元法和偽譜法三種。有限差分法是最常用的一種有效方法,但精度低,網(wǎng)格

3、步長必須明顯小于最短波長的一半,頻散現(xiàn)象也較嚴重,如采用高階差分法,能提高精度,但勢必會增加計算量,目前該方法向高階交錯網(wǎng)格法、不規(guī)則網(wǎng)格法及大網(wǎng)格快速差分法等方向發(fā)展。有限元法也是正演模擬的有效手段,能有效地求解不規(guī)則區(qū)域,但計算量相對較大,目前該方法向多重網(wǎng)格的區(qū)域分解法和混合元方向發(fā)展。偽譜法是一種較新和發(fā)展較快的方法,不同于傳統(tǒng)的有限差分法,其空間差分用快速傅氏變換來實現(xiàn),精度可認為是無限階的1,且理論上已證明,偽譜法在一個空間波長上最少僅需兩個網(wǎng)格點,因此可采用較大的空間步長,減少計算量。有限差分法,有限元法和偽譜法在地震波動方程數(shù)值模擬中都得到了較廣泛應用,但通常都用來合成共炮點地

4、震記錄及作波場傳播瞬時切片圖,如要得到零偏移距即自激自收記錄,固然可再作抽道集等處理得到,但顯然是費時和不必要的。本文用偽譜法研究了零偏移距地震記錄的合成,因其計算速度快,故可應用于C laerbou t偏移特別是三維偏移處理中。中國博士后科學基金資助第一作者簡介張文生男31歲博士后地球物理專業(yè)已發(fā)表“最小二乘法與反問題誤差最優(yōu)模的選擇”等論文收稿日期199711141方法原理1.1計算公式零偏移距地震記錄的制作要用到爆炸反射界面原理2。該原理把地下反射界面假設成具有爆炸性的物質,計算時將震源函數(shù)放于界面上,并且假定在t =0時所有的爆炸反射界面同時起爆。當介質的速度用半速度代替后,由界面上各

5、散射源發(fā)出的初至波被地面檢波器接收,其結果就是近似的CD P 疊加剖面即零偏移距記錄。假定地震波滿足如下二維標量波動方程52u 5t 2=v 2(52u 5x 2+52u 5z 2(1其中u (x ,z ,t 是位移,t 是時間,v (x ,z 是介質速度,x ,z 分別是水平和垂直距離。將(1式左端用二階時間差商代替,右端用偽譜法計算,可得u (x ,z ,t +t 2u (x ,z ,t -u (x ,z ,t -t -t 242(k 2x +k 2z v 2(x ,z u (k x ,k z ,t e i (k x x +k z z d k x d k z (2其中t 為時間采樣,k x

6、 、k z 分別為x 、z 方向的波數(shù),u (k x ,k z ,t 為u (k x ,k z ,t =u (x ,z ,t e -i (k x x +k z z d x d z (3特別地,如速度v (x ,z 為常數(shù)c ,則由(1式可導出u (x ,z ,t 的解析解3。先對(1式作關于x 、t 的二維傅氏變換,整理得52u (k x ,z ,z 2=-(2c 2-k 2x u (k x ,z ,=-k 2z u (k x ,z ,(4其中u (k x ,z ,=u (x ,z ,t e -i (k x x +t d x d t (5因速度為常數(shù),故在z 方向上的波數(shù)也為常數(shù),由(4式并利

7、用爆炸反射界面條件,經(jīng)推導可得u (k x ,z ,=u (k x ,k z ,t =0e ik z z c 1-(ck x 2(6再對(6式作關于k x 、的二維反傅氏變換,得u (x ,z ,t =142u (k x ,k z ,01c 1-(ck x 2e i (k x x +k z z +t d k x d (7由(7式即得零偏移距剖面u (x ,0,t =142u (k x ,k z ,0c 1-(ck x2e i (k x x +t d k x d (8其中k 2z =2c2-k 2x 。824長春科技大學學報1998年上面(2式是變速時的零偏移距剖面的數(shù)值計算公式,而(8式是常速

8、時的零偏移距剖面的解析表達式,這兩式均使用了全波波動方程。對某些模型或復雜模型,記錄中會出現(xiàn)多次波,為避免多次波的干擾,可采用無層間反射的單程波波動方程4。5u (k x ,k z ,t 5t=iv (x ,z k z 1+k 2x k 2z u (k x ,k z ,t (9該式也易從(1式導出(從略,其中u (k x ,k z ,t 如(3式所示。同推導(2式相類似的方法,采取時間中心差商格式,最后可得u (x ,z ,t +t u (x ,z ,t -t +ti 22u (k x ,k z ,t v (x ,z k z 1+k 2x k 2z e i (k x x +k z z d k

9、x d k z (10由(10式可算出波場u (x ,z ,t ,令z =0即得零偏移距剖面u (x ,z =0,t 。1.2震源、穩(wěn)定性條件及吸收邊界條件用方程(2、(8、(10作數(shù)值計算時,要用到初始條件,這里采用時間為雷克子波的源函數(shù)形式s (x ,z ,t ,即s (x ,z ,t =A f (x ,z co s2f 0(t -t 0e -(t -t 02f 20(11其中A 為振幅,為衰減系數(shù),f 0為主頻,f (x ,z 為界面函數(shù)。偽譜法計算的穩(wěn)定性條件可從矩陣穩(wěn)定性分析得到5,不難推得,(2式、(8式和(10式計算的穩(wěn)定性條件為t 2d v m ax (12其中空間網(wǎng)格取d =

10、x =y ,v m ax 取成最大速度值。數(shù)值模擬中另一必須解決而又難以解決的問題是吸收邊界條件。通常有限差分法采用粘滯型(或稱海綿型吸收條件,它在擴充區(qū)域內逐漸將波場衰減吸收掉,從而不產(chǎn)生邊界反射;該方法適用性廣,但難以選擇合適的吸收系數(shù)(一般用試驗方法得到,而不合適的吸收系數(shù)仍會導致很強的反射。另一種要推導邊界吸收方程,有限差分時將其結合在差分方程中計算,但它僅對一定角度的入射波有理想吸收效果。這些方法都難以應用于偽譜法中。在偽譜法中,根據(jù)傅氏變換的周期性條件,人們提出了一種反周期擴展的方法6,該方法已在各向異性介質共炮點地震記錄正演中得到了成功應用5,本文采用了這種方法。它主要通過擴展計

11、算區(qū)域增加計算量來消除記錄中的邊界反射,效果良好。2數(shù)值計算前面導出了不同情況的用于零偏移距地震記錄的公式(2式、(8式和(9式,數(shù)值計算時可根據(jù)情況利用。根據(jù)這些公式,給出如下三個模型的計算結果,加以說明和驗證。模型1模型由一個點繞射體和一個水平界面組成,其空間采樣大小為N x =64,N z =64。繞射點位于(32,7處,水平界面位于N z =20處,介質速度為3000m s 。圖1是所得的零偏移距記錄剖面,它由一個繞射雙曲線和一個水平界面構成,圖中的采樣x =z =20m ,t =2m s ,時間樣點數(shù)為250,所示單位均以空間或時間方向的采樣點數(shù)來表示(下同,其長度可由給定的采樣步長

12、換算得到。根據(jù)射線幾何路徑分析,不難驗證,圖1波場924第4期張文生,何樵登,李樹勇:用偽譜法合成零偏移距地震記錄 圖1模型1的零偏移距記錄F ig .1Zero -off set syn thetic record of m odel 1是正確的。模型2在模型1的基礎上討論一個稍復雜的模型(圖2。它模擬了一系列向斜構造(45°,速度仍假定為常數(shù)3000m s 。圖3是用(2式算得的未加邊界吸收的記錄,圖4是對應的加邊界吸收后的結果。比較可知,吸收效果明顯。圖5是用射線追蹤法所得的記錄,比較圖4與圖5可知,從另一側面說明了本方法及算法的可行性和正確性,圖中的采樣x =z =20m ,

13、t =2m s ,時間樣點數(shù)為400,空間樣點數(shù)為64。模型3如圖6,該模型由三個界面四層構成,四層的層速度分別為2000m s 、2500m s 、3000m s 、2500m s 。采用無層間反射的波動方程計算,所得的零偏移距地震記錄如圖7所示。由于受第二界面中嵌型構造的影響。第三界面中央記錄中出現(xiàn)彎曲;第二個界面所導致的繞射波形清晰可見。圖中x =z =25m ,t =2m s ,時間樣點為600,空間樣點為64,所用CPU 時間約為0.42s ,所用計算時間約為2.86m in 。圖2類似向斜構造的常速模型F ig .2A con stan t veloc ity m odel ana

14、 logued a se -r ies of syncli nes 圖3模型2(圖2的未加邊界吸收的零偏移距記錄F ig .3Zero -off set syn thetic record of m odel 2(f ig .2with no boundary absorbi ng con -d ition s034長春科技大學學報1998年 圖4模型2(圖2加邊界吸收的零偏移距記錄F ig .4Zero -off set syn thetic record of m odel 2(f ig .2af ter boundary absorbi ng 圖5模型2(圖2用射線追蹤法合成的零偏移距記

15、錄F ig .5Zero -off set syn thetic record of m odel 2(f ig .2usi ng the ray trac i ng m ethod圖6層變速模型F ig .6A m odel with layer var ied veloc ity 圖7模型3(圖6的零偏移距記錄F ig .7Zero -off set syn thetic record of m odel 3(f ig .63結論快速有效地合成零偏移距記錄是進行疊后偏移的基礎。本文給出了用全波波動方程和無層間反射的波動方程制作零偏移距地震記錄的方法,對邊界吸收用反周期擴展法作134第4期張

16、文生,何樵登,李樹勇:用偽譜法合成零偏移距地震記錄432 長春科技大學學報 1998 年 了有效處理 ( 當然不排除使用衰減型吸收邊界條件 。 模型計算表明了所用方法及算法的 正確性和有效性, 效果良好。 由于傅氏積分采用了快速離散傅氏變換來實現(xiàn), 提高了計算 速度, 盡管在邊界吸收中采用了反周期擴展法擴展了計算域增加了計算量, 但計算速度還 是快的, 這些為推廣到三維情形奠定了良好基礎。 本文得到了張關泉老師的大力幫助和支持, 同韓立國副教授在傅氏變換方面也曾作 過有益討論, 在此均表示衷心感謝。 參 考 文 獻 1 Fo rnberg B. T he p seudo sp ectra l

17、m ethod: Com p a rison s w ith fin ite d ifferences fo r the ela stic w ave equa tion. Geop hysics, 1987, 52 ( 4 : 483 501 2 L oew en tha l D , L u L , Roberson R , et a l T he w ave equa tion app lied to m ig ra tion. Geop hysica l . P ro sp ecting, 1976, 24 ( 2 : 380 390 3賀振華主編 1 反射地震資料偏移處理與反演方法 1

18、 重慶: 重慶大學出版社, 1989 ( 6 : 854 859 1997 1 2 Zhang W en sheng , H e Q iaodeng , L i Shuyong 4 Gazdag J. M odeling of the acou stic w ave equa tion w ith tran sfo rm m ethod s Geop hysics, 1981, 46 . 5張文生 1 各向異性介質彈性參數(shù)反演理論方法與數(shù)值方法研究: 學位論文 1 長春: 長春科技大學, w ith the p seudo sp ect ra l m ethod. N o t on ly w

19、e g ive the num erica l com p u ta t ion fo rm u lae . ca l ca lcu la t ion s show the co rrectness and effect iveness of m ethod s So w e supp ly a rap id m odelling m ethod fo r m ig ra t ion, w e a lso can ex tend it to th ree d i en sion s ca se ea sily. m p eriod ic ex ten sion m ethod ca ted m

20、 ed ia, and com p u ta t ion is t i e 2con sum ing fo r a 32 ca se. In the p ap er, ba sing on m D © 1995-2005 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. fo r d ifferen t ca ses, bu t a lso so lve the bounda ry ab so rb ing p rob lem effect ively. N um eri2 6 Fu rum u ra T , T akenaka H. A nw rap a round eli ina tion techn ique fo r the p seudo sp ectra l w ave syn 2 m ly, it is aecom p lished by the ray t racing m ethod, how ever, ray s m ay cro ss fo r a com p li2 the 22 sca la r