唐 文 （中國石油大學(xué) （北京）CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京102249西安石油大學(xué)地球科學(xué)與工程學(xué)院，陜西 西安710065）

王尚旭 （中國石油大學(xué) （北京）CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京102249）

引起偏移成像精度不高的因素中，討論較多的是良好的成像算法、適當(dāng)?shù)钠瓶讖?、精確的速度模型、滿意的原始資料以及較好的成像條件[1～3]。對(duì)于子波延遲所造成的影響仍缺乏詳盡的討論。在介質(zhì)的速度參數(shù)為已知的條件下，確定反射圖像的任務(wù)就是求反射點(diǎn)的空間位置及其反射因數(shù)。由于無法求出確切的反射因數(shù)，成像的反射因數(shù)實(shí)際上就是用反射波與入射波振幅比來表示。偏移成像的目的就是把地面上觀測(cè)到的反射波和繞射波歸位到產(chǎn)生它的反射點(diǎn)上去，這需要進(jìn)行上行波場(chǎng)的反向外推。外推后求出的各點(diǎn)波場(chǎng)值，有的是來自本點(diǎn)的反射波，有的是該點(diǎn)下方許多點(diǎn)上的反射波疊加。因此，需要在外推波場(chǎng)中提取成像值。反射成像原則就是：反射面由反射點(diǎn)組成，在這些點(diǎn)上，其入射波的初至與反射波的產(chǎn)生時(shí)間相同，反射因數(shù)為反射波振幅u（x，z，td）（其中，x、z分別為水平坐標(biāo)、深度坐標(biāo)，m；td為反射時(shí)間，s）與入射波振幅d（x，z，td）的比值。因此，反射波成像可以用上下行波場(chǎng)進(jìn)行互相關(guān)提取成像值。下行波初始時(shí)間的確定決定了成像的位置和精度。通常假設(shè)下行波是最小相位子波或脈沖波時(shí)，其成像位置是精確的[4，5]；如果下行波是一個(gè)有短延續(xù)長度的子波時(shí)，它只是一個(gè)很好的近似成像公式，其成像位置有一定的誤差。為此，筆者以子波波峰出現(xiàn)時(shí)間作為初始時(shí)間，以一個(gè)繞射點(diǎn)模型為例，選用有限元法波動(dòng)方程疊前逆時(shí)偏移程序進(jìn)行模擬與偏移，利用互相關(guān)成像條件，就子波延遲時(shí)間的不同對(duì)成像位置精度產(chǎn)生的影響進(jìn)行了分析。

1 有限元逆時(shí)外推方程

時(shí)空域二維常密度聲波方程為：

式中，P ＝P（x，z，t）為聲壓，Pa；v＝v（x，z）為波速，m/s；x、z為空間變量，m；t為時(shí)間變量，s；q為震源項(xiàng)。采用部分離散，即只對(duì)空間域進(jìn)行離散，令是關(guān)于節(jié)點(diǎn)j的插值函數(shù)；αj（t）是節(jié)點(diǎn)j上的時(shí)間函數(shù)；n是節(jié)點(diǎn)數(shù)）。應(yīng)用加權(quán)余量法導(dǎo)出式（1）的等效積分弱形式的伽遼金有限元控制方程為：

式中，M為質(zhì)量矩陣；K為剛度矩陣；C為阻尼矩陣；t為場(chǎng)向量對(duì)時(shí)間t的二階導(dǎo)數(shù)；at為場(chǎng)向量對(duì)時(shí)間t的一階導(dǎo)數(shù)；at為場(chǎng)向量；Qt為t時(shí)刻節(jié)點(diǎn)上的場(chǎng)向量。應(yīng)用中心差分法得到顯式逆時(shí)外推公式[6，7]為：

式中，Δt為時(shí)間間隔，s。

2 互相關(guān)成像原理

成像條件的制定是地震偏移算法的關(guān)鍵之一，直接決定了成像效果和計(jì)算成本。選用目前普遍認(rèn)為成像效果較好的互相關(guān)成像條件公式[8～10]：

式中，Map（x，z）為像空間；Pr（x，z，t）為接收波場(chǎng)聲壓；Ps（x，z，t）為震源波場(chǎng)聲壓。被積函數(shù) Pr（x，z，t）、Ps（x，z，t）分別在t時(shí)刻對(duì)整個(gè)波場(chǎng)作一次成像運(yùn)算，積分則說明像空間Map（x，z）中的像是各時(shí)間步互相關(guān)成像值的疊加。

逆時(shí)偏移互相關(guān)成像條件的成像步驟是：

1）首先要完成一次波動(dòng)方程正演計(jì)算，并保存每一時(shí)間步的波場(chǎng)信息。

2）逆時(shí)外推觀測(cè)波場(chǎng)。

3）在每一時(shí)間步，讀取剛才保存的相同時(shí)刻的震源波場(chǎng)，再與當(dāng)前的外推觀測(cè)波場(chǎng)作互相關(guān)（即Pr（x，z，t）＊Ps（x，z，t）），再將得到的像加入像空間中。

4）重復(fù)步驟2）、3），直到t＝0。

在已進(jìn)入上市輔導(dǎo)且IPO申請(qǐng)未被受理的395家新三板企業(yè)中，輔導(dǎo)期在半年以內(nèi)的企業(yè)共有54家，占比13.67%。也就是說，除已申報(bào)和已進(jìn)入輔導(dǎo)期的新三板企業(yè)外，新三板還有數(shù)量龐大的IPO后備軍。這些企業(yè)一部分因?yàn)闃I(yè)績暫時(shí)無法IPO，一部分出于新三板政策紅利的考慮，則處于觀望搖擺階段。

應(yīng)用式（1）、（2）要涉及到選擇上、下行波的初始時(shí)間，筆者選擇雷克子波，以其波峰出現(xiàn)時(shí)間作為初始時(shí)間。

3 子波延遲引起成像誤差分析

設(shè)計(jì)一個(gè)點(diǎn)繞射模型來分析子波延遲對(duì)成像精度的影響，如圖1所示，A是地下一繞射點(diǎn)，繞射點(diǎn)周圍速度是v。以自激自收為例，S是炮點(diǎn)，R是接收點(diǎn)，炮點(diǎn)S到反射點(diǎn)A的傳播時(shí)間為tSA，反射點(diǎn)A到接收點(diǎn)R傳播時(shí)間為tRA。雙程波疊前逆時(shí)偏移成像值是震源波場(chǎng)和接收波場(chǎng)每次外推到同一時(shí)刻就要進(jìn)行互相關(guān)并提取成像值，如果震源波場(chǎng)子波和接收波場(chǎng)子波為脈沖波或最小相位子波時(shí)，則波場(chǎng)外推到A點(diǎn)能準(zhǔn)確成像。如果震源子波波峰相對(duì)于接收波場(chǎng)子波有一個(gè)延遲時(shí)間tBA，則震源波場(chǎng)和接收波場(chǎng)延拓到A點(diǎn)時(shí)，震源波場(chǎng)Ps子波波峰在B點(diǎn)，而接收波場(chǎng)Pr子波波峰到達(dá)了A點(diǎn)，在tA時(shí)刻相關(guān)不能成像，只有當(dāng)震源波場(chǎng)延拓經(jīng)過A點(diǎn)到達(dá)C點(diǎn)即外推時(shí)間為時(shí)，震源波場(chǎng)波峰到達(dá)D點(diǎn)深度（D是BA中點(diǎn)），即接收波場(chǎng)外推到時(shí)刻到達(dá)C點(diǎn)時(shí)，其波峰也到達(dá)C點(diǎn)深度，此刻D、C點(diǎn)在同一深度，這時(shí)互相關(guān)提取的成像值深度在DC水平線深度上，相對(duì)于正確成像位置A點(diǎn)，深度減小了為繞射點(diǎn)上方介質(zhì)速度，m/s）。

圖1 震源子波延遲tBA的互相關(guān)成像示意圖

如圖2所示，A是地下一繞射點(diǎn)，D為淺層成像點(diǎn)。當(dāng)接收波場(chǎng)Pr子波波峰相對(duì)于震源波場(chǎng)Ps子波波峰延遲時(shí)間為tBA，tRD＝tBA。震源波場(chǎng)開始傳播時(shí)刻，接收波場(chǎng)子波波峰已經(jīng)在D點(diǎn)，震源波場(chǎng)向下延拓tSE達(dá)到E點(diǎn)，接收波場(chǎng)向上延拓tDE也到達(dá)E點(diǎn)，2種波場(chǎng)在此刻相關(guān)能得到一個(gè)成像值，其深度為。當(dāng)波場(chǎng)延拓到A點(diǎn)的時(shí)刻即tSA，震源波場(chǎng)子波波峰到達(dá)A點(diǎn)，接收波場(chǎng)子波波峰因?yàn)檠舆t則波前已到達(dá)B點(diǎn)，此刻波場(chǎng)互相關(guān)，則在A點(diǎn)不能成像，只有當(dāng)接收波場(chǎng)從B向上傳播時(shí)間為到達(dá)C點(diǎn)，震源波場(chǎng)經(jīng)過A點(diǎn)傳時(shí)間到達(dá)C點(diǎn)（C是AB中點(diǎn)），即時(shí)刻2種波場(chǎng)互相關(guān)在C點(diǎn)成像，C點(diǎn)深度相對(duì)于真實(shí)位置A點(diǎn)深度增加了為繞射點(diǎn)下方介質(zhì)速度，m/s）。所以接收波場(chǎng)子波相對(duì)于震源波場(chǎng)子波有一個(gè)正向延遲，則會(huì)多出現(xiàn)一個(gè)假成像點(diǎn)。

通過類似分析，當(dāng)震源波場(chǎng)和接收波場(chǎng)的子波延遲時(shí)間相同時(shí)，對(duì)成像精度沒有影響即成像點(diǎn)在準(zhǔn)確位置。

圖2 接收記錄子波延遲tBA的互相關(guān)成像示意圖

4 模型試算

為了驗(yàn)證上述理論分析的正確性，選取零相位雷克子波波峰出現(xiàn)時(shí)間為初始時(shí)間，設(shè)計(jì)點(diǎn)繞射模型圖進(jìn)行測(cè)試。圖3是速度模型，在深度500m處有一個(gè)繞射點(diǎn)，其周圍速度是4000m/s。在繞射點(diǎn)正上方布置震源，兩邊布置檢波器，道間距10m，共101個(gè)檢波器接收。

圖4是震源子波沒有延遲時(shí)模擬的地震記錄。采用互相關(guān)成像條件成像，使用有限元法波動(dòng)方程疊前逆時(shí)偏移程序進(jìn)行偏移。當(dāng)選擇的震源波場(chǎng)子波延遲時(shí)間為0.05s時(shí)，成像位置在深度為400m的地方，比實(shí)際深度少100m（如圖5所示）。當(dāng)選擇的震源波場(chǎng)子波波峰延遲時(shí)間為0.1s時(shí)，成像位置在深度為300m的地方，比實(shí)際深度少200m（如圖6所示）。

圖7是震源波場(chǎng)子波波峰延遲0.05s時(shí)生成的地震記錄。選用震源波場(chǎng)子波波峰延遲0.05s的正向延拓波場(chǎng)與該地震記錄逆時(shí)延拓波場(chǎng)進(jìn)行互相關(guān)成像，其成像位置仍在深度為500m的地方 （如圖8所示）。當(dāng)選用震源波場(chǎng)子波波峰無延遲時(shí)的正向傳播波場(chǎng)與圖7地震記錄反向傳播的波場(chǎng)進(jìn)行互相關(guān)成像時(shí)，則成像位置在深度600m的地方，比實(shí)際位置深100m （如圖9所示）。根據(jù)前面的理論分析，成像剖面在深度為100m的位置有一個(gè)假成像點(diǎn)，由于是單炮逆時(shí)偏移，淺部成像點(diǎn)由于能量弱未能顯示出來。

圖3 速度模型 （繞射點(diǎn)在深度500m處，繞射點(diǎn)周圍速度v＝4000m/s）

圖4 震源波場(chǎng)子波沒有延遲時(shí)的單炮地震記錄

圖5 震源波場(chǎng)子波延遲0.05s時(shí)的圖4逆時(shí)偏移結(jié)果 （成像深度400m）

圖6 震源波場(chǎng)子波延遲0.1s時(shí)的圖4偏移結(jié)果 （成像深度300m）

圖7 震源波場(chǎng)子波延遲0.05s的單炮地震記錄

5 結(jié) 論

圖8 震源波場(chǎng)子波延遲0.05s時(shí)的圖7逆時(shí)偏移結(jié)果 （成像深度500m）

圖9 震源波場(chǎng)子波波峰無延遲時(shí)的圖7逆時(shí)偏移結(jié)果 （成像點(diǎn)深度600m）

1）如果震源外推波場(chǎng)子波與觀測(cè)波場(chǎng)子波波峰延遲時(shí)間相同，則互相關(guān)成像位置與實(shí)際界面位置符合，成像準(zhǔn)確。

2）震源外推波場(chǎng)所用子波波峰相對(duì)于觀測(cè)波場(chǎng)子波波峰有延遲時(shí)間t，則互相關(guān)成像深度比實(shí)際深度小h，用公式表示為

3）觀測(cè)波場(chǎng)子波波峰相對(duì)于震源外推波場(chǎng)子波波峰延遲時(shí)間t時(shí)，則會(huì)產(chǎn)生2個(gè)成像點(diǎn)，一個(gè)成像點(diǎn)深度在真實(shí)深度下方處；另一個(gè)在地表下方處（其中，v表示淺部成像點(diǎn)上方速度）。c

使用互相關(guān)成像條件的疊前逆時(shí)偏移方法處理實(shí)際資料時(shí)，如果以子波波峰出現(xiàn)時(shí)間作為初始時(shí)間，成像誤差與震源波場(chǎng)子波波峰和地震記錄子波波峰相對(duì)延遲時(shí)間以及反射點(diǎn)周圍速度有關(guān)，子波波峰相對(duì)延遲時(shí)間越大，反射點(diǎn)周圍速度越大，成像點(diǎn)位置誤差就越大。因此進(jìn)行偏移時(shí)，選取正向傳播震源波場(chǎng)子波波峰延遲時(shí)間與地震資料子波波峰延遲時(shí)間相同的子波，有利于確定成像點(diǎn)的準(zhǔn)確位置。

[1]Parkes G E，Hatton L.Towards a systematic understanding of the effects of velocity model errors on depth and time migration of seismic data [A].1984SEG Annual Meeting [C].Atlanta，1984－12－2～6.

[2]Zhu Jinming，Lines L，Gray S.Smiles and frowns in migration velocity analysis [J].Geophysics，1998，63 （4）：1200～1209.

[3]Liu Zhenyue，Bleistein N.Migration velocity analysis：Theory and an iterative algorithm [J].Geophysics，1995，60 （1）：142.

[4]Claerbout J F.Toward a unified theory of reflection mapping [J].Geophysics，1971，36 （3）：467～481.

[5]Claerbout J F.Imaging the Earth＇s Interior [M].Oxford：Blackwell Scientific Publications，1985.

[6]薛東川，王尚旭 .起伏地表復(fù)雜介質(zhì)波動(dòng)方程有限元數(shù)值模擬方法 [J].地球物理學(xué)進(jìn)展，2007，22（2）：522～529.

[7]薛東川，王尚旭 .波動(dòng)方程有限元疊前逆時(shí)偏移 [J].石油地球物理勘探，2008，43（1）：17～21.

[8]Jrg S，JesséC C，Amélia N.A comparison of imaging conditions for wave－equation shot－profile migration [J].Geophysics，2008，73 （6）：219～227.

[9]Chattopadhyay S，McMechan G A.Imaging conditions for prestack reverse－time migration[J].Geophysics，2008，73 （3）：81～89.