地震數(shù)值模擬畢業(yè)設(shè)計(jì)

1、地震數(shù)值模擬研究及應(yīng)用摘要本論文是對(duì)數(shù)值模擬中波動(dòng)方程相移法及相移加插值方法的應(yīng)用。使用相移法對(duì)地下的地層模型進(jìn)行正演和偏移，但是相移法只針對(duì)縱向速度變化的模型。為了得到更好的地質(zhì)模型的效果，引入了相移加插值的方法對(duì)水平變化的復(fù)雜模型進(jìn)行延拓，最終跟相移法進(jìn)行比對(duì)，分析了兩種方法的優(yōu)勢(shì)與不足。在對(duì)相移法的實(shí)驗(yàn)過程中，采用了模型簡(jiǎn)化，以求取平均速度的方式采用深度延拓，然后對(duì)點(diǎn)脈沖進(jìn)行正演以及偏移，印證了實(shí)驗(yàn)的正確性后，又采用了這種方法對(duì)其他復(fù)雜模型進(jìn)行了實(shí)驗(yàn)，另外還在相移的基礎(chǔ)上使用了相移加插值的方法，對(duì)模型進(jìn)行了檢驗(yàn)和分析。通過一系列的試算后，發(fā)現(xiàn)相移加插值的方法比起相移法精度更高，并且在傾角

2、很大的速度模型上，相移加插值的方法更能夠識(shí)別地層的有效信息，并且在一些復(fù)雜的模型上也體現(xiàn)了出來，但是，相移加插值有計(jì)算量較大的缺點(diǎn)。關(guān)鍵詞：波場(chǎng) 偏移 脈沖 相移 插值The research and application of seismic numerical simulationAbstract: The research and application of the phase shift and interpolation of wave equation in numerical simulation is studied in this paper. The forward m

3、odeling and migration of the underground strata are carried out by using the phase shift method, but the phase shift method is only used to model the longitudinal velocity change. In order to get better geological models, the introduction of the reference velocity interpolation method on changes in

4、the level of complex model extension, eventually with phase shift method for comparison, analyses advantages and disadvantages of the two methods. In the phase shift method in the course of the experiment, the model is simplified, in order to obtain the mean velocity with depth extrapolation and to

5、point pulse forward modeling and migration, confirms the correctness of the experiment and uses this method to other complex model experiments were carried out, still another phase shift based on using the reference velocity of wave field interpolation, the model were test and analysis. Through a se

6、ries of trial and error analysis found higher than phase shift precision interpolation method of phase shift, and the inclination angle is very big velocity model and interpolation method can effective recognition of stratigraphic information, and in some complex models is also reflected in it. Howe

7、ver, phase shift interpolation calculation of large amount of defects.Keywords: wave field migration pulse phase shift interpolation II目錄摘要IABSTRACTII第一章：前言11.1 課題研究背景及意義11.2 研究現(xiàn)狀11.3主要研究?jī)?nèi)容3第二章 波動(dòng)方程偏移成像42.1偏移42.1.1 偏移的概念42.1.2 偏移脈沖響應(yīng)72.2 波動(dòng)方程偏移的成象原理72.2.1 爆炸反射界面成象原理82.3 相移波動(dòng)方程偏移92.3.1 基本概念及公式92.3.2

8、傅里葉變換的微分性質(zhì)92.3.3 單程聲波方程波場(chǎng)延拓公式102.3.4 相移加插值延拓13第三章 程序算法流程163.1 數(shù)字化163.1.1 頻率-波數(shù)域離散化163.1.2 邊界反射問題163.1.3 插值計(jì)算173.2 程序流程173.2.1 疊后相移正演流程173.2.2 疊后相移偏移流程183.2.3疊后相移插值正演流程183.2.4疊后相移插值偏移流程18第四章 疊后深度偏移算例19結(jié)論及建議25致謝26參考文獻(xiàn)27第一章：前言1.1 課題研究背景及意義地震數(shù)值模擬在當(dāng)前地球物理學(xué)勘探中有著廣泛的應(yīng)用和相當(dāng)?shù)姆萘浚?dāng)前的勘探技術(shù)正在往更加先進(jìn)，要求更加高的層面上發(fā)展，在這種情況下

9、就要求科研人員們對(duì)于一些特別的地質(zhì)模型，復(fù)雜的地質(zhì)要求下利用已經(jīng)學(xué)過的知識(shí)采用一些特別的成熟的勘探方法，還要有一些檢測(cè)工具來印證這些手段。用地震數(shù)值模擬是一個(gè)有效的方法，可以來解決現(xiàn)代開發(fā)工作、油氣勘探中所面臨的棘手問題，它也具有重要的意義。首先我們來回顧地震數(shù)值模擬的定義：即已知地下的參數(shù)及結(jié)構(gòu)，模擬波在地質(zhì)體中的傳播，并計(jì)算相對(duì)于地面對(duì)地下的觀測(cè)地震波記錄。地震數(shù)值模擬方法在地震勘探中得到了廣泛的應(yīng)用，地震數(shù)值模擬在地震勘探的各個(gè)階段都具有十分重要的作用。在勘探的采集過程之中，可以運(yùn)用這種方法對(duì)其觀測(cè)系統(tǒng)進(jìn)行系統(tǒng)的評(píng)估，可以依此對(duì)觀測(cè)系統(tǒng)進(jìn)行改進(jìn)。在物探數(shù)據(jù)的各個(gè)處理的過程中，可以對(duì)反演方

10、法的精度進(jìn)行有效的評(píng)估。在地震資料處理結(jié)果的解釋中，地震數(shù)值模擬可以檢驗(yàn)數(shù)據(jù)的結(jié)果。地震勘探的是為了研究地下地質(zhì)結(jié)構(gòu)而存在，并根據(jù)地震記錄在井或地面的不同地點(diǎn)所觀測(cè)到的數(shù)據(jù)，描述其狀態(tài)或巖性。這些都屬于反演的概念，但是，反演過程是建立在地震正演模型的基礎(chǔ)上。因此，地震數(shù)值模擬不僅可以進(jìn)行地震正演模擬，也是地震反演的基礎(chǔ)。1.2 研究現(xiàn)狀從過去到現(xiàn)在，地震波場(chǎng)數(shù)值模擬技術(shù)都是地震勘探領(lǐng)域的重要研究?jī)?nèi)容，波場(chǎng)數(shù)值模擬雖然沒有像使用物理模型那樣的真實(shí)，但它具有靈活的參數(shù)，方便的模型描述，較少的人力以及無需特定材料等特點(diǎn)而被廣泛地應(yīng)用，特別是計(jì)算機(jī)的迅速發(fā)展，各種軟件的交互化，建模技術(shù)的齊備，地震數(shù)值

11、模擬技術(shù)得到了更加迅猛的發(fā)展。目前的地震數(shù)值模擬方法包括射線追蹤和波動(dòng)方程法。前者主要是將地震波理論簡(jiǎn)化為射線理論，并考慮地震波傳播的運(yùn)動(dòng)學(xué)特征?；诓▌?dòng)方程的常規(guī)成像方法主要有三種：傅立葉域法、基爾霍夫積分法和有限差分法。射線法具有概念清晰、操作方便、顯示直觀等優(yōu)點(diǎn)，而且其適應(yīng)性強(qiáng)，其缺點(diǎn)是還有一定的條件限制，計(jì)算還不是很精確，傳統(tǒng)的射線追蹤法有試射法和彎曲的方法，其局限性是處理高速度的變化很難，很難找到全局最小走時(shí)多值時(shí)間，計(jì)算效率低，還有陰影射線覆蓋密度問題。本論文主要研究的是波動(dòng)方程的傅氏域方法進(jìn)行正演和偏移，下面來介紹一下傅氏域法波場(chǎng)數(shù)值模擬的一些發(fā)展及應(yīng)用：根據(jù)前人的波動(dòng)方程的理論

12、，如果說地震在傳播之中均滿足于波動(dòng)方程的規(guī)律。如果說我們已經(jīng)知道了地面上的速度模型，我們可以通過延拓求解地震波場(chǎng)，再現(xiàn)波在地下介質(zhì)中傳播的過程。1972年JClearbout首先提出波動(dòng)方程有限差分偏移法，給我們留下了現(xiàn)在還在使用的15度有限差分的方法。即時(shí)這種方法對(duì)于一些構(gòu)造角度較大的地區(qū)并不能適用，但我們現(xiàn)在有很多的人仍然在采用他的這種偏移方法。1984年Jeno Gazdag使用了橫向均勻速度介質(zhì)的相移法外推的方法。中間結(jié)果是P波參考波場(chǎng)。他還使用了另一種方法，實(shí)際波場(chǎng)通過參考波場(chǎng)的計(jì)算再插值出來的，相移加插值的方法是很穩(wěn)定，方便的，適合于三維地震數(shù)據(jù)的偏移，這是使用這種方法研究的先例。

13、1990年Stoffa提出了分裂傅利葉的方法，但如果橫向速度過大，處理結(jié)果也會(huì)有明顯的偏差，1994年Ristow提出了傅利葉有限插分法，將兩種方法很聰明地結(jié)合起來，可以處理速度橫向變化較大的情況。下面介紹一些國內(nèi)的創(chuàng)新研究實(shí)例。1992年涂國田采用相移插值法中遞歸向下延拓方式，每一延拓層間的成像使用了給出的常速F-K再加上剩余速度F-K 偏移的方法，再與橫向速度變化有關(guān)插值來實(shí)現(xiàn)。1998年賀振華，熊高君等提出地震波傳播的定位原理，他們首次利用了聲波方程得到了地震勘探中人工合成共炮點(diǎn)記錄。2001年李錄明使用延拓和插值即實(shí)現(xiàn)了波場(chǎng)延拓模型的校正，實(shí)現(xiàn)了雙曲線校正時(shí)差。2003年李冰提出了一種

14、基于相移算子約束的離散光滑插值的構(gòu)造一維顯式短算子的優(yōu)化方法通過離散光滑插值法，在頻率一波數(shù)域中，以傳播區(qū)內(nèi)的相移算子作為約束，在傳播區(qū)外的算子的兩端處利用零點(diǎn)為約束，進(jìn)行離散的光滑插值，使得所得的算子能夠具有二階光滑可導(dǎo)性，那么它對(duì)應(yīng)的頻率一空間域中的算子就可以取很短。2005年王祥春在采用變換坐標(biāo)系中波動(dòng)方程波場(chǎng)延拓的方法來消除地表起伏的影響。2011年杜玉靜實(shí)現(xiàn)了時(shí)間 一空間域和頻率一波數(shù)域的地震波場(chǎng)外推。在時(shí)間一空間域?qū)崿F(xiàn)了角度較大的65度波動(dòng)方程有限差分的偏移，在頻率一波數(shù)域?qū)崿F(xiàn)了只適合地震波速度，僅隨深度變化的相移法偏移和適應(yīng)地震波速度橫向變化的分裂步有限差分偏移?？偟膩碚f，數(shù)值模

15、擬技術(shù)有很多自己的特點(diǎn)和優(yōu)勢(shì)，運(yùn)算速度和模擬精度是數(shù)值模擬技術(shù)發(fā)展的兩個(gè)主要的方向，簡(jiǎn)單快速的模擬方法可以快速地檢驗(yàn)各種處理解釋成果的正確性，應(yīng)用波場(chǎng)數(shù)值模擬技術(shù)不僅僅可以在勘探階段確定使用的處理解釋方法 并且可以利用其振幅對(duì)儲(chǔ)層的敏感性監(jiān)測(cè)油氣開采情況可以預(yù)測(cè)，隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展 波場(chǎng)方程的數(shù)值模擬技術(shù)的應(yīng)用范圍會(huì)越來越廣。1.3主要研究?jī)?nèi)容1.利用數(shù)值模擬的方法構(gòu)建一個(gè)具有地質(zhì)意義的地層模型，在模型剖面上做出一個(gè)脈沖，在實(shí)際工作中即為震源，脈沖處于中間位置，波速為4000m/s左右，道間距10米，采樣間隔4ms，延拓步長(zhǎng)5m。方便進(jìn)行波形延拓與插值計(jì)算，最開始使用水平地層模型，然后可

16、以構(gòu)造一個(gè)傾斜地層或者不規(guī)則地層。2.利用爆炸反射界面的原理，將一個(gè)點(diǎn)脈沖看做是由每個(gè)地層的脈沖波，而上一級(jí)是下面幾層的疊加，可以用付式變換在頻率波數(shù)域中用C語言，通過延拓公式和傅里葉變換，用來編寫點(diǎn)脈沖的波場(chǎng)延拓程序和偏移結(jié)果，得到正演與偏移后的脈沖響應(yīng)。3.由于需要研究水平速度變化的情況，采用幾個(gè)參考速度，用延拓公式和差值計(jì)算，將起始的參考速度向下延拓，經(jīng)過插值，計(jì)算，便可以計(jì)算研究每個(gè)地層的速度變化。最后每一個(gè)地層都有其對(duì)應(yīng)的水平速度變化，當(dāng)然，插值擬合的一些算法在地球物理學(xué)勘探中已經(jīng)被各種實(shí)踐證明中是可行的，在國內(nèi)外各種勘探中也已被廣泛應(yīng)用。第二章 波動(dòng)方程偏移成像2.1偏移2.1.1

17、 偏移的概念地震勘探中以激發(fā)地震波收到他們從地震的剖面信號(hào)作為基本方法，經(jīng)過地下地質(zhì)異常（彈性不均勻體或波阻抗界面）反射，然后利用儀器接收地震記錄，再推導(dǎo)成地下的地質(zhì)剖面，這樣做可以得到地下剖面或深度圖。這些地圖也可以用作發(fā)現(xiàn)地質(zhì)剖面研究上。地震記錄所定義的空間被稱為圖像空間或數(shù)據(jù)空間。地下地震地質(zhì)圖或地下埋深剖面稱為目標(biāo)空間或地質(zhì)空間。數(shù)據(jù)空間可以表示為：Px, y, t,其中x, y為空間坐標(biāo)，t為時(shí)間坐標(biāo)。目標(biāo)空間或地質(zhì)空間記為Px, y,z，其中P為表征地質(zhì)結(jié)構(gòu)的特性參數(shù)(如界面反射系數(shù)，彈性參數(shù)等)。一般來說，我們期望地震記錄可以直接作為地質(zhì)剖面，但因?yàn)榇嬖诘卣鹩涗涍x取方式的不同、觀

18、測(cè)系統(tǒng)的構(gòu)成形式有差別、地下界面的結(jié)構(gòu)和傳播的不同、影響因素、這些都存在一定的差異，因此用地震記錄來解釋地下的地質(zhì)剖面就成了一個(gè)有風(fēng)險(xiǎn)的事情。給出一個(gè)具體的例子來說明這種情況。xz/tAFOD 'DBB'M反射界面圖2-1 自激自首的情況下，反射同相軸B'D'相對(duì)它的反射界面BD沿下傾方向移動(dòng)了，長(zhǎng)度加大與地面夾角變小，并且存在關(guān)系式sin=tan，這些變化統(tǒng)稱為偏移如果說地下有一個(gè)反射界面BD(見圖2-1),我們?cè)O(shè)炮檢距為零， BD的反射同向軸為 B'D'，在前面我們已經(jīng)講到，BD、B'D'都在目標(biāo)域和象域之中。為了方

19、便起見，設(shè)波速在介質(zhì)中為一常數(shù)，通過h=Vt2，我們可以進(jìn)行時(shí)深轉(zhuǎn)換，則時(shí)間剖面(象域)和地質(zhì)界面(目標(biāo)域)就可以使用一張圖來顯示(圖2-1)，但是B'D'和BD在圖上看來并不是重合的。由B引BD的垂直線段交x于A點(diǎn)，則AB=h表示從A點(diǎn)激發(fā)信號(hào)再在A點(diǎn)觀測(cè)的來自B的反射路徑， B點(diǎn)反射回來的波應(yīng)于A點(diǎn)的下面的點(diǎn)B'，我們可以設(shè)波的單程反射時(shí)間是t0，那么就有公式AB'=t0V=AB。同理，從D引BD的垂直線段交水平面x于F點(diǎn)，得到了FD=FD'。當(dāng)然，我們說界面BD，它在象域中的B'D'(反射同相軸)和它自身的各個(gè)位置，狀態(tài)都是不一樣的

20、，簡(jiǎn)單的說，一般研究人員稱之為地震的偏移。在實(shí)際問題中，勘探的地質(zhì)地形一般都比較復(fù)雜，這就會(huì)使得野外施工所用的方法會(huì)存在許多誤差，讓我們對(duì)地下的構(gòu)成認(rèn)識(shí)十分困難，如果說我們要得到一個(gè)比較靠譜的結(jié)論的話，就要進(jìn)行它的矯正或者是偏移處理，也就是說將B'實(shí)現(xiàn)歸位到B處，D'歸位到D處。如果說要進(jìn)一步的研究，可以考慮一些偏移距或者角度的一些思想。偏移距是說一個(gè)反射的點(diǎn)(比如說B)和這個(gè)點(diǎn)反射的(B')之間的距離，BB'于水平的投影MB可以說是水平上的偏移距，那么在垂直方向上的投影MB'稱為垂直偏移距。若將BD和B'D'延長(zhǎng)，則比較于地面的O點(diǎn)，這

21、是因?yàn)镺點(diǎn)既是反射的點(diǎn)，它又是接收點(diǎn)(此時(shí)波的傳播時(shí)間為零)。設(shè)BDO與x軸的夾角為，B'D'O與x軸的夾角為，則按簡(jiǎn)單的幾何關(guān)系得sin=tan (2-1)此式描述了偏移處理前后同相軸視傾角和界面真傾角的關(guān)系。2. 波的傳播特性與偏移設(shè)均勻介質(zhì)中有一個(gè)自激自收的激發(fā)點(diǎn)(或散射點(diǎn))，波從繞射點(diǎn)出發(fā) (見圖2-2)。二維情況下波的傳播方程式為x-xd2+z-zd2=V2t2 (2-2)圖2-2 均勻介質(zhì)中繞射波及假想觀測(cè)面(zi=0，z1，z2，z3)式中xd，zd是繞射點(diǎn)坐標(biāo)，V，t是傳播的速度及旅行的時(shí)間。當(dāng)t固定時(shí)，其軌跡是一個(gè)圓，因?yàn)閂2t2=R2=常數(shù)時(shí)(2-2)變?yōu)閳A

22、方程。當(dāng)觀測(cè)面固定(z=常數(shù))，(2-2)式可改寫為t=zi-zd2V2+x-xd2V2 (2-3)式中zi=0，z1，z2，zd，(2-3)是一雙曲線方程，對(duì)不同的zi，雙曲線形狀與位置不同，見圖2-3，這樣可以發(fā)現(xiàn)其中有一些共同的特點(diǎn)：在橫軸的坐標(biāo)點(diǎn)xd和激發(fā)點(diǎn)的水平坐標(biāo)xd其實(shí)都是一樣的的。另外，可以發(fā)現(xiàn)，當(dāng)z=zd時(shí)，xd的時(shí)間數(shù)據(jù)為0。這就可以說明：當(dāng)我們的觀測(cè)的位置是zd時(shí)，如果地震波的傳播時(shí)間是0的話就可以標(biāo)定xd, zd。事實(shí)上，在地下深度為zd處進(jìn)行觀測(cè)是不可能的事情，然而如果再經(jīng)過數(shù)學(xué)的運(yùn)算可以將在地面上的數(shù)據(jù)Px, z=0, t變到z=zd處，得到記錄Px, z=zd,

23、t，然后在令t=0，即得到Pxd, zd，這就是最終要求的激發(fā)點(diǎn)。按常理來說，由Px, z=0, t計(jì)算出Px, z, t這樣一個(gè)方式叫做延拓(外推)，由Px, z, t求Px, z叫做波場(chǎng)的成像，這兩個(gè)在地震偏移處理中是很重要的方法。上面所說的是單個(gè)點(diǎn)的繞射情況，其他的地震界面也可以用到這個(gè)方法，由于地震界面可以看做是點(diǎn)的合集。每一個(gè)繞射點(diǎn)都有其對(duì)應(yīng)的曲線，雙曲線簇漸近線是反射同相軸因?yàn)閠an=dzdx sin=dzdx=Vdtdx=VP (2-4)這個(gè)公式中P是射線的參數(shù)。因?yàn)閐tdx為我們雙曲線的斜率。這樣就將圖2-4和圖2-1聯(lián)系起來了。綜上所述，那么就可以對(duì)偏移的概念綜合起來：地震記

24、錄上的反射同相軸因?yàn)槭懿ǖ膫鞑ヌ匦缘挠绊懞陀涗浄绞降南拗仆c其相應(yīng)的反射地質(zhì)體在形態(tài)和位置上不一致，這種不一致性稱為偏移，設(shè)法消除偏移影響的方法叫偏移處理。xtoxdz=0xtoxdz=z1xtoxdz=z2xtoxdz=zd圖2-3 繞射點(diǎn)xd, zd處，觀測(cè)面在zi=0，z1，z2，zd時(shí)測(cè)得的時(shí)距圖，注意：z=zd時(shí)，時(shí)距曲線變?yōu)轫旤c(diǎn)在t=0的直線2.1.2 偏移脈沖響應(yīng)1. 一個(gè)脈沖在剖面的形成如果一個(gè)剖面(象域)中只有一個(gè)脈沖點(diǎn)，其他的都是0時(shí)通過偏移以后其對(duì)應(yīng)的地層中的空間(目標(biāo)空間)中的圖形被我們叫做輸入剖面偏移脈沖響應(yīng)。假定波速不變，剖面是自己激發(fā)自己接收的，那么它的偏移的圖

25、形為一個(gè)半圓形構(gòu)造(見圖2-4)。因?yàn)榈叵碌钠拭嫒羰且缘孛娴囊粋€(gè)點(diǎn)作為圓心的一個(gè)半圓，利用自激自收的方法來進(jìn)行測(cè)量，反射波就會(huì)形成一個(gè)脈沖波顯示在時(shí)間的剖面之上?？傮w來說，一個(gè)地震記錄從時(shí)間的剖面轉(zhuǎn)化為目標(biāo)空間我們就叫做偏移，從目標(biāo)空間變到時(shí)間剖面可以稱之為地質(zhì)模型或正演模擬。xtoxdxtoxdsRJ偏移處理V=常數(shù)正演模擬ba圖2-4a. 輸入的自激自收時(shí)間剖面(象空間)；b. V=常數(shù)時(shí)，偏移脈沖響應(yīng)為半圓，2. 一個(gè)脈沖響應(yīng)的方程如果說地下空間有一個(gè)點(diǎn)的脈沖響應(yīng)，地面上得到的觀測(cè)數(shù)據(jù)(象域，時(shí)間域)就叫做我們輸出的時(shí)間剖面的響應(yīng)。當(dāng)然，自激自收的話就是一個(gè)雙曲線方程(見圖2-2及圖2-

26、3)。其為t=4zd2V2+4x-xd2V2 (2-5)式中2zdV=t0是用實(shí)雙程旅行時(shí)。2.2 波動(dòng)方程偏移的成象原理如圖2-6所示，如果說(z=0)存在著激發(fā)與接收的兩個(gè)點(diǎn)，兩個(gè)的位置是一樣的，同時(shí)收到A點(diǎn)的反射回來的波，當(dāng)波的波速V在是一個(gè)定值時(shí)，射線為直線。那么反射的波就在S的下面B點(diǎn)處，時(shí)間為t=2SA/V，SB=Vt/2=SA，可以說B和反射的點(diǎn)之間的長(zhǎng)度為偏移距。假定我們可以把測(cè)線歸于z=z1處，要有A點(diǎn)的反射，那么兩個(gè)點(diǎn)的位置都要在S'處，這個(gè)時(shí)候波的旅行時(shí)間t'=2S'A/V,反射波可以放在B'處，且S'B'=AS'，

27、這個(gè)時(shí)候偏移距變?yōu)锳B'，同樣，要是觀測(cè)的位置在z=z2，那么AB''又是它的偏移距。由圖2-6可以發(fā)現(xiàn)，觀測(cè)的位置變化的過程中，有兩個(gè)特征：（1）偏移距越來越?。唬?）旅行時(shí)間越來越短。這也就是說波場(chǎng)向下延拓的時(shí)候，總可以讓偏移距變?yōu)?，實(shí)現(xiàn)研究的最終偏移的目標(biāo)。偏移距AB和波的旅行路徑AS，地層傾角的關(guān)系為 AB=2ASsin2 (2-6)當(dāng)觀測(cè)線向下延拓到A點(diǎn)時(shí)，路徑AS=0，偏移距AB=0。地層傾角為零時(shí)，偏移距AB恒為零，此時(shí)根本就不用進(jìn)行偏移處理。SS''S'ABB'B''VZ=0Z=Z1Z=Z2圖2-5 延拓

28、與偏移關(guān)系圖，當(dāng)測(cè)線向下延拓時(shí)，A點(diǎn)反射波的旅行路徑AS逐步減，偏移距AB也逐步減小2.2.1 爆炸反射界面成象原理這個(gè)原理是數(shù)值模擬中很常用的一種方法。即將地下的反射界面的所有狀態(tài)都跟地面上的狀態(tài)一致，其產(chǎn)生的地震波也都是一致的。并且假定在t=0時(shí)刻，所有的界面都激發(fā)，到達(dá)地面上，其上行波被地面所接受到，波的傳播速度為Ve。如果我們利用波動(dòng)方程將波場(chǎng)(地震記錄)向下延拓到t=0時(shí)的波場(chǎng)值就實(shí)現(xiàn)了其自動(dòng)歸位，自動(dòng)實(shí)現(xiàn)偏移成像。oxzoxzVVe=V2反射界面爆炸反射界面ab圖2-6a. 自激自收反射界面模型，上行波和下行波路徑相同，波的傳播速度為V；b. 爆炸反射界面模型，t=0時(shí)，波自界面發(fā)

29、出(只有上行波)，波的傳播速度Ve=V/2，由于Ve比V小一倍，因此模型a和b中波到達(dá)地面的時(shí)間相同此成象原理適用于水平疊加后的地震資料的偏移處理，因?yàn)榀B加后資料相當(dāng)于零炮檢距地震剖面，自炮點(diǎn)發(fā)出的下行波到達(dá)反射點(diǎn)的路徑與自該點(diǎn)反射返回地面的上行波的路徑完全一致，這樣，可以考慮上行波而不必考慮下行波。但是實(shí)際記錄上的反射波到達(dá)時(shí)間是雙程時(shí)間，若僅考慮上行波，波的到達(dá)時(shí)間將減少一倍，即Ve=V/2。爆炸反射界面成象原理概況在圖2-6中。2.3 相移波動(dòng)方程偏移相移波動(dòng)方程偏移是由J.Gazdag(1978)提出的，對(duì)于只有縱向速度變化的情況，它的結(jié)果也是準(zhǔn)確的。2.3.1 基本概念及公式1. 波

30、場(chǎng)：pressurePz, x, t設(shè)x, z為空間坐標(biāo)，t為時(shí)間，地震波傳播速度為Vx, z，則二維介質(zhì)中任意位置、任意時(shí)刻的地震波場(chǎng)為Pz, x, t。2. 頻率f：?jiǎn)挝粫r(shí)間內(nèi)包含的地震波的周期T數(shù) f=1T , =2f=2/T (2-7)3. 波數(shù)(k)：?jiǎn)挝豢臻g(距離)內(nèi)包含的波長(zhǎng)()數(shù) k=2/ (2-8)4. 頻率-空間域波場(chǎng)：Pz, x, ，Pz, x, t只對(duì)時(shí)間t做傅里葉變換到頻率域 Pz, x, =-Pz, x, te-itdt (2-9)5. 頻率-波數(shù)域波場(chǎng)：頻率域波場(chǎng)Pz, x, 對(duì)x做傅里葉變換到波數(shù)k域 Pk, x, =-Pz, x, e-ikxdx (2-10)

31、2.3.2 傅里葉變換的微分性質(zhì)1. 傅里葉變換的微分性質(zhì)若P=-Pte-itdt -傅里葉變換 (2-11)P=12-Peitdt-反傅里葉變換 (2-12)則 iP=-dPtdte-itdt (2-13)i2P=-d2Ptdt2e-itdt (2-14)2. 利用傅里葉變換微分特性求波場(chǎng)的各偏導(dǎo)數(shù)的傅里葉變換波場(chǎng)空間偏導(dǎo)數(shù)對(duì)時(shí)間的傅里葉變換：-2Pz, x, tx2e-itdt=2Pz, x, x2 (2-15)-2Pz, x, tz2e-itdt=2Pz, x, z2 (2-16)波場(chǎng)時(shí)間偏導(dǎo)數(shù)對(duì)時(shí)間傅里葉變換：由式(3-8)得i2Px,z,=-d2Px,z,tdt2e-itdt (2-

32、17)波場(chǎng)空間偏導(dǎo)數(shù)對(duì)空間的傅里葉變換：由式(3-8)同理得ik2Pk=-d2Pxdx2e-ikxdx (2-18)則ik2Pk,z,t=-d2Px,z,tdx2e-ikxdx (2-19)速度隨深度變換介質(zhì)的波場(chǎng)的時(shí)間傅里葉變換：-d21Vz2Px,z,tdt2e-itdt=ik2Vz2Px,z, (2-20)2.3.3 單程聲波方程波場(chǎng)延拓公式二維各向同性均勻介質(zhì)中的地震波傳播方程為：2Px,z,tx2+2Px,z,tz2=1V2x,z2Px,z,tt2 (2-21)若速度只隨深度變化，則可化為：2Px,z,tx2+2Px,z,tz2=1V2z2Px,z,tt2 (2-22)對(duì)t和x做傅里

33、葉變換： 先對(duì)t做傅里葉變換的頻率域-2Px,z,tx2+2Px,z,tz2e-itdt=-1V2z2Px,z,tt2e-itdt (2-23)得到頻率-空間域波動(dòng)方程：2Px, z, x2+2Px, z, z2=i2V2zPx,z, (2-24) 再對(duì)x做傅里葉變換-2Px, z, x2+2Px, z, z2e-ikxdx=-i2V2zPx,z,e-ikxdx (2-25)得到頻率-波數(shù)域波動(dòng)方程：ik2Pk,z,+2Pk, z, z2=i2V2zPk,z, (2-26)由頻率-波數(shù)域方程得：2Pk, z, z2=-2V2z-k2Pk,z, (2-27)令kz=2V2z-k2 (2-28)得

34、2Pk, z, z2=-kz2Pk,z,= ikz2Pk,z, (2-29)回顧常系數(shù)二階微分方程的解：d2fxdx2=-a2fx 或 d2fxdx2=ia2fx (2-30)其解為：fx=c1e-iax+c2e+iax (2-31)比較常系數(shù)二階微分方程(2-26)與頻率-波數(shù)域方程(2-27)，可知頻率-波數(shù)域波動(dòng)方程為變系數(shù)二階偏微分方程，其中變化的系數(shù)為kz=2V2z-k2。所以要想按照常系數(shù)二階微分方程來解頻率-波數(shù)域方程，需要將頻率-波數(shù)域方程轉(zhuǎn)換為常系數(shù)微分方程。下面使用常系數(shù)問題的解決方案來求解頻率-波數(shù)域波動(dòng)方程的解：速度Vz只隨深度變化，但假定在一定深度間隔內(nèi)是常數(shù)，即在z

35、i和zi+1間隔內(nèi)傳播的速度Vz不隨z變，是z的常數(shù)，即如圖2-7的層狀介質(zhì)。v1v2vnzizi+1圖2-7 層狀介質(zhì)模型那么，式(3-23)對(duì)z的微分中，令：Pk,z,= fz，得d2fzdz2=ikz2fz (2-32)由于在zi和zi+1間隔為zi內(nèi)傳播的速度Vz不隨z變，是z的常數(shù)，則式(2-32)的解為：fx=c1e-ikzzi+c2e+ikzzi (2-33)由于Pk,z,= fz，所以Pk,zi,=c1e-ikzzi+c2e+ikzzi (2-34)在相位移波動(dòng)方程偏移中，使用下行波方程Pk,zi,=c1e-ikzzi (2-35)或Pk,zi,=c1e-ikzzi+1-zi

36、(2-36)式(2-35)或(2-36)就是頻率-波數(shù)域波動(dòng)方程的單程波場(chǎng)延拓公式。比較常系數(shù)二階微分方程(2-27)與頻率-波數(shù)域方程(2-28)，可知頻率-波數(shù)域波動(dòng)方程為變系數(shù)二階偏微分方程，本文討論的是零炮檢距(疊加后)二維地震剖面的偏移問題。如果說是水平層狀的模型，且假定在zizzi+zi的深度間隔內(nèi)波的傳播速度Vi不變，在zi+1zzi+1+zi+1的深度間隔內(nèi)，波速度Vi+1也是常數(shù)，但是不要求速度Vi和Vi+1必須相等，這樣，在各個(gè)深度間隔內(nèi)，頻率-波數(shù)域中標(biāo)量波動(dòng)方程上行波的解如式(2-36),可寫為PkX,zi+zi,=PkX,zi,e-ikzizi (2-37)式中zi為

37、波場(chǎng)從深度z=zi的水平面上延拓到z=zi+1水平面的距離，zi=zi+1-zi。而PkX,zi,等于PkX,zi,=-Px,zi,te-it-kXxdxdt (2-38)kzi=Vzi1-VzikX212 (2-39)當(dāng)zi=0時(shí)，Px,0,t和PkX,0,，即為地面上的疊加記錄和疊加記錄的二維付氏變換。將式(2-38)對(duì)kX做一維付氏反變換得Px,zi+zi,=12-PkX,zi,e-ikzizi+kXxdkX (2-40)上式對(duì)再做一次一維付氏反變換Px,zi+zi,t=12-Px,zi+zi,eitd (2-41)按照爆炸反射界面成象原理，令t=0，同時(shí)在式(3-33)中將速度Vzi減

38、半，得到偏移結(jié)果Px,zi+zi=12-Px,zi+zi,d (2-42)實(shí)際計(jì)算的時(shí)候，因?yàn)榭偸怯邢薜?，且采用離散值，上式變?yōu)镻x,zi+zi=12n=-N2N2Px,zi+zi,n (2-43)波動(dòng)方程相位移偏移處理使用單頻率延拓，這就需要在每個(gè)延拓的步長(zhǎng)之內(nèi)，對(duì)波的各個(gè)頻率成分重復(fù)進(jìn)行。頻率分量的多少是與計(jì)算的速度成正比。因此，我們可以利用實(shí)信號(hào)的傅里葉正負(fù)頻率的共軛特性將(2-43)式進(jìn)行化簡(jiǎn)，則可以節(jié)約一倍的運(yùn)算工作量。具體如下。設(shè)Pt的頻譜為P，如果Pt為實(shí)函數(shù)，則P=P- (2-44)式中的橫線代表復(fù)共軛。又設(shè)P=A+iB (2-45)則P-=A-iB (2-46)式中A和B均為

39、實(shí)數(shù)。于是(2-43)式化為Px,zi+zi=n=-N2N2Px,zi+zi,n=2n=0N2Ax,zi+zi,n=2n=0N2Re<Px,zi+zi,n> (2-47)最后的表達(dá)式中，設(shè)=1。如果將n的變化范圍限制在比0<n<n2的一個(gè)更窄一些的范圍內(nèi)，不但會(huì)使計(jì)算時(shí)間進(jìn)一步減少，還對(duì)記錄做了帶通濾波。2.3.4 相移加插值延拓設(shè)地震波的傳播速度在深度間隔z內(nèi)是橫向可變、縱向不變的，即V=V(x)，在這種情況下將z=zi深度處的波場(chǎng)Px,zi,t延拓到z=zi+z=zi+1深度，由于V(x)的變化，無法在波數(shù)域中進(jìn)行延拓。因此，設(shè)想用幾個(gè)不同的常速V1,V2,V3對(duì)波

40、場(chǎng)Px,zi,t分別延拓，如果在z間隔內(nèi)，波的傳播速度不變，任設(shè)Vx=V1，顯然，利用上節(jié)所述的方法，可求得zI+1處頻率、波數(shù)域延拓后的波場(chǎng)P1kx,zi+1,=Pkx,zi,e-ikzizi (2-48)同樣，如果分別假設(shè)z間隔內(nèi)地震波的傳播速度為Vx=V1, Vx=V2, Vx=V3，可求得P2kx,zi+1,，P3kx,zi+1,。對(duì)P1kx,zi+1,|l=1,2,3關(guān)于kx，做二維傅氏反變換可得P1kx,zi+1,|l=1,2,3。如果三個(gè)參考速度選擇得當(dāng)，根據(jù)與它們對(duì)應(yīng)的三個(gè)波場(chǎng)值P1,P2,P3，按插值方法可求出所有的以不同速度傳播的延拓波場(chǎng)值Pkx,zi+1,。從而近似地解決

41、了橫向變速時(shí)波場(chǎng)的延拓問題。為了近一步理解波場(chǎng)插值的物理意義，下面用一個(gè)簡(jiǎn)單的例子說明，設(shè)Px,zi,t=Pzi,t為一正弦波，暫時(shí)不考慮x方向的變化。，zi+1 -zi=z=10m參考速度V1=400m/s, V2=500m/s, V3=700m/s求波速V=600m/s時(shí)延拓波場(chǎng)Pzi,t。所謂向下延拓，實(shí)際上等于波的反向傳播，因此波場(chǎng)P1zi+1,t等于將Pzi,t反向傳播t1=zv1=25ms時(shí)所取得的值，同理，可求得t2=zv2=20ms，P2和t3=zv3=14ms， P3，如圖2-8所示。圖2-8波場(chǎng)插值原理圖一般情況下，采用拉格郎日插值公式即能滿足精度要求。如果選用三個(gè)參考速度

42、，應(yīng)取三點(diǎn)插值公式,得到Pzi+1,t為Pzi+1,t= t-t2t-t3t1-t2t1-t3P1zi+1,t+t-t1t-t3t2-t1t2-t3P2zi+1,t+ t-t1t-t2t3-t1t3-t2P3zi+1,t (2-49)由于與速度成反比，式(2-2-7)還可表示為Pzi+1,t= 1V-1V21V-1V31V1-1V21V1-1V3P1zi+1,t+(1V-1V1)(1V-1V3)(1V2-1V1)(1V2-1V3)P2zi+1,t+ (1V-1V1)(1V-1V2)(1V3-1V1)(1V3-1V2)P3zi+1,t (2-50)如果只選用了兩個(gè)參考速度，應(yīng)使用兩點(diǎn)線性插值公式

43、Pzi+1,t=1-qP1zi+1,t+qP2zi+1,t (2-51)式中，q=1V-1V11V2-1V1 (2-52) (2-46)和(2-47)式對(duì)t做傅氏變換，且用慢度S 代替速度的倒數(shù)，得到相應(yīng)的頻率域插值公式：Pzi+1,t= S-S2S-S3S1-S2S1-S3P1zi+1,t+(S-S1)(S-S3)(S2-S1)(S2-S3)P2zi+1,t+ (S-S1)(S-S2)(S3-S1)(S3-S2)P3zi+1,t (2-53)其中kz=v11-v2kx2=v11-k12,k10.95式中S, S1, S2, S3分別為1V,1V1,1V2,1V3。 第三章 程序算法流程3.1

44、 數(shù)字化根據(jù)前面的疊后偏移成像的概念以及其意義，要實(shí)現(xiàn)其算法流程還是需要對(duì)其進(jìn)行一個(gè)數(shù)字化的處理，將它的速度及反射系數(shù)進(jìn)行離散化和初始化。首先將初始速度與波場(chǎng)定義成公式（3-1）所示。Px,t=PIx×x,It×t (3-1-a)Rx,z=RIx×x,Iz×z (3-1-b)Vx,z=VIx×x,Iz×z (3-1-c)在程序中，可以使用數(shù)組的方式來實(shí)現(xiàn)每一層每一道的算法。3.1.1 頻率-波數(shù)域離散化按照時(shí)間和頻率的抽樣定理，一個(gè)受限的頻率信號(hào)在不大于12tm頻率的間隔進(jìn)行采樣，那么抽樣得到的離散信號(hào)可以用來代表原信號(hào)，那么我們可以

45、使頻率和波數(shù)的采樣間隔分別為：=mNt，kx=kxmNx，由于傅里葉變換的共軛對(duì)稱性，所以時(shí)間序列與頻率域可以先取至一半進(jìn)行計(jì)算。3.1.2 邊界反射問題由上述步驟，可以做任意復(fù)雜構(gòu)造的正演模擬,但有缺陷。最重要的缺陷是邊界反射問題，即將實(shí)際無窮空間區(qū)域中求解波場(chǎng)的問題化為有窮區(qū)域求解時(shí)，左右兩邊使用零邊界條件。物理上假設(shè)探區(qū)距與兩個(gè)端點(diǎn)很遠(yuǎn)，在兩個(gè)端點(diǎn)上收到的反射波很弱。但是，上述條件在實(shí)際中不能成立，造成零邊界條件反而成為絕對(duì)阻止波通過的強(qiáng)反射面。在正演模擬的剖面上出現(xiàn)了邊界假反射干涉正常界面的反射，反射現(xiàn)象如圖3-1所示。圖3-1 疊后偏移正演（未鑲邊）如果說在剖面的兩端加上若干道記錄，

46、雖然無疑增加了計(jì)算難度。在這個(gè)繞射點(diǎn)的實(shí)例中采用削波的方法進(jìn)行處理，即波的兩端乘以一個(gè)越來越小的相移因子進(jìn)行振幅減弱處理。3.1.3 插值計(jì)算在這里，采用（2-50）的計(jì)算公式來計(jì)算水平速度變化的模型，在中間使用它的最大速度與最小速度，分別使用最大速度與最小速度計(jì)算相移因子并延拓，再使用拉格朗日公式進(jìn)行插值，使用的延拓因子的公式仍然是正演模擬中的公式。3.2 程序流程根據(jù)以上算法與流程我們可以對(duì)算法進(jìn)行一個(gè)算法梳理，以便算法的執(zhí)行。3.2.1 疊后相移正演流程1. 參數(shù)初始化，對(duì)每層速度求取平均值。2. 形成邊界削波數(shù)據(jù)。3. 波場(chǎng)初始化。4. Zmax層波場(chǎng)延拓到深度Zmax-1 。5. Z

47、i+1層波場(chǎng)延拓到深度Zi 。6. 重復(fù)5，從Iz=Nz-1開始，直到Iz=1，得測(cè)線上的頻率空間域波場(chǎng)。7. 頻率-空間域波場(chǎng)對(duì)頻率做反傅里葉變換，得時(shí)間-空間波場(chǎng)。8. 換取一個(gè)頻率，重復(fù)4到7的流程。3.2.2 疊后相移偏移流程1. 對(duì)輸入的時(shí)間剖面波場(chǎng)對(duì)時(shí)間進(jìn)行一維傅氏變換得到空間頻率域波場(chǎng)。2. 對(duì)波場(chǎng)以為參數(shù)在x方向進(jìn)行吸收衰減。3. 對(duì)波場(chǎng)進(jìn)行傅氏變換，得到波數(shù)頻率域的波場(chǎng)。4. 計(jì)算相移因子，注意速度減半。5. 使用相移因子乘以當(dāng)前波場(chǎng)。6. 對(duì)當(dāng)前波場(chǎng)的進(jìn)行一維傅里葉反變換，向深度域進(jìn)行延拓。7. 返回步驟2，繼續(xù)深度循環(huán)，直到最大深度時(shí)為止。8. 再取一個(gè)頻率值，重復(fù)步驟2

48、7，直到所有頻率循環(huán)完畢為止。3.2.3疊后相移加插值正演流程1. 在插值的方法中，我們的流程與上訴的流程不同，需要先將波場(chǎng)轉(zhuǎn)化至?xí)r間波數(shù)域場(chǎng)。2. 對(duì)波場(chǎng)進(jìn)行深度的延拓。3. 對(duì)波場(chǎng)以為參數(shù)在x方向進(jìn)行吸收衰減。4. 求取每層最大速度與最小速度，計(jì)算波場(chǎng)相移因子并延拓。5. 利用拉格朗日插值進(jìn)行波場(chǎng)的插值6. 對(duì)當(dāng)前波場(chǎng)的進(jìn)行一維傅里葉反變換，向深度域進(jìn)行延拓。7. 返回步驟2，繼續(xù)深度循環(huán)，直到最大深度時(shí)為止。8. 再取一個(gè)頻率值，重復(fù)步驟27，直到所有頻率循環(huán)完畢為止。3.2.4疊后相移加插值偏移流程1. 讀取初始速度文件和波場(chǎng)文件。2. 將正演波場(chǎng)轉(zhuǎn)化為空間頻率域，進(jìn)行角頻率的循環(huán)。3

49、. 深度循環(huán)，利用成像公式進(jìn)行延拓。4. 求取最大速度與最小速度，分別計(jì)算相移因子，計(jì)算波場(chǎng)與相移因子的乘積。5. 轉(zhuǎn)為頻率波數(shù)域，用最大與最小速度計(jì)算延拓因子并延拓。6. 對(duì)當(dāng)前波場(chǎng)的進(jìn)行一維傅里葉反變換，向深度域進(jìn)行延拓。7. 返回步驟2，繼續(xù)深度循環(huán)，直到最大深度時(shí)為止。8. 轉(zhuǎn)為空間域進(jìn)行拉格朗日的插值計(jì)算，得到結(jié)果波場(chǎng)。第四章 疊后深度偏移算例如圖4-1所示，模型剖面上的一個(gè)脈沖圖，二維數(shù)據(jù)體中點(diǎn)脈沖的位置為（128，128），波速為3000m/s，道間距為10m，采樣間隔為4ms，延拓步長(zhǎng)5m。圖4-2為常速介質(zhì)中存在單個(gè)繞射點(diǎn)時(shí)的相移法疊后正演記錄雙曲線，從偏移的角度來說，這張圖

50、表示了繞射點(diǎn)到地面?zhèn)鞑r(shí)間的不同，此雙曲線稱為正演脈沖響應(yīng)，圖4-3為其相應(yīng)的偏移剖面。從這幾個(gè)結(jié)果可以更加清楚地知道地震記錄和正演的概念性上的不同，圖4-2上顯示了此繞射點(diǎn)的時(shí)距曲線，即其到達(dá)地面各觀測(cè)點(diǎn)的時(shí)間，而經(jīng)過偏移之后得到了其真實(shí)的地質(zhì)剖面：一個(gè)繞射點(diǎn)，相反，如果對(duì)輸入的一個(gè)點(diǎn)記錄作為偏移記錄，那么得到的是一個(gè)雙曲線的剖面，在這里，也對(duì)其正演記錄做出了褶積，使其記錄及剖面變得更加明顯，同樣的，在對(duì)復(fù)雜地質(zhì)體使用F-K偏移時(shí)，其每一個(gè)點(diǎn)均可以看作一個(gè)繞射點(diǎn)，由此也可以判斷出程序的正確性與實(shí)用性。xz圖4-1 初始繞射點(diǎn)模型xt圖4-2 繞射點(diǎn)正演響應(yīng)zx圖4-3 一個(gè)繞射點(diǎn)的偏移剖面對(duì)

51、于復(fù)雜的地質(zhì)模型來說如下圖生物礁模型所示，相移正演的方法是使用求取每層平均速度，對(duì)其進(jìn)行深度的延拓，最終得到圖4-5的正演圖形，而圖4-7得到的是偏移之后的圖形，從這里可以看出記錄所得到的時(shí)距曲線以及偏移之后的剖面顯示的十分清晰。而圖4-6及圖4-8所顯示的是相移加插值之后的正演和偏移的結(jié)果，對(duì)比未經(jīng)插值采用平均速度方法的模型，其顯示的更加清晰與完整，從正演及偏移后的剖面的細(xì)節(jié)比起相移法更加完好，但如果說選用未采用平均速度使用相移法，剖面會(huì)變得及其不準(zhǔn)確，選用平均速度來簡(jiǎn)化整個(gè)界面是一個(gè)十分有效的方法。在地質(zhì)模型的正演及偏移之中，地質(zhì)模型斷層外的能量團(tuán)顯示的不是很清晰，分析其中存在的原因可能是

52、斷層上下的速度差的差別本身不是很大的緣故，也可能是相移的方法對(duì)于水平速度的劇烈變化不是非常敏感。zx圖4-4生物礁斷層模型txtx 圖4-5生物礁的相移正演模型 圖4-6同一地質(zhì)模型相移加插值正演結(jié)果zzxx 4-7 相移偏移剖面圖 圖4-8 生物礁模型相移加插值偏移結(jié)果為了更好地比對(duì)兩種方法的不同，假定一個(gè)水平速度變化劇烈的情況，即采用一個(gè)橫向雙層速度模型的相同的點(diǎn)脈沖響應(yīng)，顯然，使用經(jīng)典的相移算法可以使用平均速度的方法得到圖4-1的繞射點(diǎn)記錄，可以看出，相移加插值對(duì)于剖面的響應(yīng)比起相移法更加敏感，準(zhǔn)確，在界面之上發(fā)生了偏折，如圖4-10所示，所以相移加插值的優(yōu)勢(shì)更加地突顯出來，使用平均速度簡(jiǎn)化地質(zhì)模型的方法也出現(xiàn)了缺陷，而相移加插