陳可洋，陳樹民，李來林，吳清嶺，范興才，劉振寬

（1.中國石油大慶油田有限責(zé)任公司 勘探開發(fā)研究院，黑龍江大慶163712；2.中國石油大慶油田有限責(zé)任公司勘探事業(yè)部，黑龍江 大慶163453）

0 引言

地震波正演數(shù)值模擬技術(shù)一直是國際地球物理學(xué)界的熱點(diǎn)內(nèi)容之一。隨著地震波動(dòng)理論和計(jì)算機(jī)技術(shù)的不斷發(fā)展，自20世紀(jì)60年代以來該項(xiàng)技術(shù)便得到了飛速發(fā)展，其中采用波動(dòng)方程的地震波數(shù)值模擬技術(shù)能更好地保持地震波的幾何學(xué)、運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)等特征，因此可達(dá)到精確模擬地震波傳播規(guī)律的目的［1］。目前，地震波正演數(shù)值模擬技術(shù)的相關(guān)理論已基本成熟，并已形成了一系列經(jīng)典算法和研究成果，如在數(shù)值離散方面［2-5］，以交錯(cuò)網(wǎng)格有限差分法最具代表性，有效壓制了數(shù)值頻散并提高了計(jì)算效率；在吸收邊界方面［6-8］，以最佳匹配層（PML）吸收邊界條件最具代表性，在減小計(jì)算量的同時(shí)減少了虛假波場的引入，有效提高了計(jì)算結(jié)果的信噪比。地震波正演數(shù)值模擬技術(shù)的引入不僅大大降低了物理模擬技術(shù)和實(shí)際野外地震勘探的成本，同時(shí)還可有效指導(dǎo)實(shí)際地震資料的采集、處理、解釋及方法驗(yàn)證等，因此，該技術(shù)在當(dāng)前勘探地震學(xué)領(lǐng)域中發(fā)揮著至關(guān)重要的作用。

隨著地震波數(shù)值模擬技術(shù)的進(jìn)步，如何在地震波正演模擬過程中實(shí)現(xiàn)不同類型波場的自動(dòng)識別與分離問題逐漸得到廣泛重視，有利于對復(fù)雜波場傳播規(guī)律的認(rèn)識和應(yīng)用。其中，最具代表性的是彈性波波場分離數(shù)值模擬方法，目前該方法包括基于散度算子計(jì)算純縱波和基于旋度算子計(jì)算純橫波的波場分離數(shù)值模擬方法［9］，也有構(gòu)建等效彈性波方程的縱橫波波場分離數(shù)值模擬方法［10］，它們均可實(shí)現(xiàn)純縱波和純橫波波場的自動(dòng)分離。同時(shí)，研究介質(zhì)也從單相介質(zhì)延伸到雙相孔隙介質(zhì)。截至目前，未見方向行波波場分離數(shù)值模擬方面的相關(guān)報(bào)道，但針對地震方向行波波場分離的相關(guān)方法早已得到應(yīng)用，如利用單程波偏移算子進(jìn)行地震成像［11］和利用單程波吸收邊界條件進(jìn)行地震波正演中的邊界吸收處理［12］等。在前人研究的基礎(chǔ)上，筆者采用3種模型實(shí)例，基于波印廷矢量，開展地震波場中上行波、下行波、左行波和右行波共4種類型方向行波的波場分離數(shù)值模擬實(shí)驗(yàn)，以期提高對地震波波場傳播規(guī)律的認(rèn)識，并指導(dǎo)實(shí)際地震資料數(shù)字處理。

1 基本理論

地震波波動(dòng)方程通?？杀硎鰹橐浑A雙曲型形式［6］

式中：P為質(zhì)點(diǎn)的應(yīng)力波場，無量綱；u和w分別為2個(gè)空間方向的質(zhì)點(diǎn)振動(dòng)速度，m/s；v為介質(zhì)速度，m/s；f為震源函數(shù)，無量綱；xs和zs分別為震源在模型x和z空間方向的坐標(biāo)，m；t為波場傳播時(shí)間，s；ρ為密度，kg/m3。

式（1）的具體數(shù)值離散方法、吸收邊界條件及其穩(wěn)定性條件可參見文獻(xiàn)［1］。

Yoon等［11］給出了用于地震波場的波印廷矢量計(jì)算公式，即

該波印廷矢量代表地震波在介質(zhì)中傳播時(shí)其波前面的法線傳播方向。筆者在此基礎(chǔ)上對波印廷矢量進(jìn)行了波場數(shù)值特征分析，并基于該矢量波場實(shí)現(xiàn)對原始波場的分離，最終實(shí)現(xiàn)了左行波、右行波、上行波和下行波的波場分離正演數(shù)值模擬，且計(jì)算過程簡單易行。

2 方向行波分離數(shù)值算例

2.1均勻介質(zhì)模型

模型大小為1 km×1 km，橫、縱向空間步長均為5m。采用最大頻率為60Hz的雷克子波作為震源，并置于模型中央位置處激發(fā)。地震波速度為2 000 m/s，介質(zhì)密度為 1.1×103k/m3，時(shí)間步長為0.5ms，滿足計(jì)算所需的穩(wěn)定性條件［7］，波場快照記錄時(shí)間為0.2 s。

分析圖1（a）可知，在各向同性介質(zhì)中，地震波波場快照的波前面為圓形。圖1（b）為水平方向波印廷矢量波場快照，藍(lán)色波場數(shù)值為正，紅色波場數(shù)值為負(fù)，且正好沿過震源的垂直線分為2瓣，即劃分為左行波波場和右行波波場；同理，圖1（c）為垂直方向波印廷矢量波場快照，正好沿過震源的水平線分為2瓣，即劃分為上行波波場和下行波波場。根據(jù)圖 1（b）和圖 1（c）對圖 1（a）完整的地震波波場快照進(jìn)行篩選，最終得到了完全分離出來的左行波波場［圖 1（d）］、右行波波場［圖 1（e）］、上行波波場［圖 1（f）］和下行波波場［圖 1（g）］。根據(jù)計(jì)算結(jié)果分析可知，在各向同性介質(zhì)中本文所提的方向行波波場分離數(shù)值模擬方法準(zhǔn)確可靠。

圖1 0.2 s時(shí)刻均勻介質(zhì)地震行波分離的波場快照及其波印廷矢量Fig.1 Seismic one-way wave separating snapshots in isotropic media and corresponding Poynting vector at 0.2 s

2.2 傾斜界面速度模型

圖2 復(fù)雜介質(zhì)模型Fig.2 Complex velocity model

圖2為含有多個(gè)角度的傾斜界面速度模型，模型大小為2.5 km×1.0 km，橫、縱向空間步長均為5m。模型頂層①和底層⑦厚度分別為150m和250m，對應(yīng)的地震波速度分別為1 800m/s和3 500m/s。模型中間部分包含一組不同傾角的傾斜界面，而且傾斜界面自右向左的傾角依次為 30°，45°，60°和90°，模型中介質(zhì)②～⑥對應(yīng)的地震波速度分別為2 000m/s，2 300m/s，2 600m/s，3 000m/s 和 2 700m/s，介質(zhì)密度均為1.1×103k/m3。采用最大頻率為60Hz的雷克子波作為震源，在模型（1.0 km，0.2 km）位置處激發(fā)，合成記錄長度為1 s，時(shí)間步長為0.5ms，滿足計(jì)算所需的穩(wěn)定性條件［7］。地面采集深度為100m，VSP采集位置距模型左側(cè)100m，波場快照記錄時(shí)間為0.2 s。

圖3（a）和圖3（b）分別為水平方向和垂直方向的波印廷矢量波場快照，其中藍(lán)色部分波場數(shù)值為正，紅色部分波場數(shù)值為負(fù)。分析圖3（a）和圖3（b）可知，在這種較為復(fù)雜的速度模型中，不同方向地震行波波場的傳播特征非常復(fù)雜，無法有效識別哪些波場屬左行波波場，哪些波場屬右行波波場，哪些波場屬上行波波場以及哪些波場屬下行波波場，而基于波印廷矢量可以對不同類型的地震方向行波波場實(shí)現(xiàn)自動(dòng)識別與分離。圖3（c）為0.2 s時(shí)刻的地震波全波場快照，圖 3（d）～（g）分別為對圖 3（c）的全波場快照進(jìn)行方向行波波場分離的結(jié)果，最終得到了完全分離的左行波、右行波、上行波和下行波波場快照。圖3（h）和圖3（j）分別為100m深度位置采集和VSP觀測方式采集的模擬記錄。從這2種模擬記錄中無法有效識別與分離不同類型的方向行波波場。 圖 3（i）和圖 3（k）分別為圖 3（h）和圖3（j）對應(yīng)的波印廷矢量數(shù)值模擬記錄，圖 3（l）～（o）分別為該方法分離出的左行波、右行波、上行波和下行波波場數(shù)值模擬記錄。綜合分析可知，該方法適合較為復(fù)雜的介質(zhì)模型，而且可以獲得準(zhǔn)確的左行波、右行波、上行波和下行波波場分離結(jié)果。

圖3 復(fù)雜介質(zhì)模型中地震方向行波波場分離的波場快照、模擬記錄及其波印廷矢量Fig.3 Seismic one-way wave separating snapshots, numerical records and corresponding Poyntingvector in complex velocity model

2.3 M armousi模型

為了進(jìn)一步驗(yàn)證方向行波波場分離數(shù)值模擬方法在更為復(fù)雜的介質(zhì)模型中的應(yīng)用效果，以Marmousi模型進(jìn)行模擬。模型大小為3.4 km×1.4 km，橫、縱向空間步長均為5m。地震波最小速度為1 028m/s，最大速度為4 670m/s，采用最大頻率為60Hz的雷克子波作為震源，在（1.70 km，0.15 km）處激發(fā)，合成記錄長度為2 s，時(shí)間步長為0.2ms，滿足計(jì)算所需的穩(wěn)定性條件。地面采集深度為100m，VSP采集位置距模型左側(cè)1 km，波場快照記錄時(shí)間為0.6 s。

圖 4（a）和圖 4（b）分別為 Marmousi模型水平方向和垂直方向的波印廷矢量波場快照。分析圖4（a）和圖4（b）可知，由于介質(zhì)的復(fù)雜性，在該介質(zhì)的波場快照中左行波、右行波、上行波和下行波波場的傳播特征變得更為復(fù)雜，而基于波印廷矢量的地震波波場仍可有效地分離出不同類型的方向行波波場。圖 4（c）為 0.6 s時(shí)刻的地震波全波場快照，圖 4（d）～（g）分別為對圖4（c）的完整地震波波場快照進(jìn)行方向行波波場的分離處理，最終得到了完全分離的左行波、右行波、上行波和下行波波場快照。圖4（h）和圖4（j）分別為Marmousi模型100m深度位置采集和VSP觀測方式采集的地震波正演模擬記錄。從這2種模擬記錄中也無法有效識別和分離不同類型的方向行波波場。圖4（i）和圖4（k）分別為圖4（h）和圖 4（j）對應(yīng)的波印廷矢量模擬記錄，圖 4（l）～（o）分別為該方法分離出的左行波、右行波、上行波和下行波的波場模擬記錄。Marmousi模型數(shù)值分析表明，該方法的應(yīng)用效果也準(zhǔn)確有效。

圖4 M armousi模型中地震方向行波分離的波場快照、模擬記錄及其波印廷矢量Fig.4 Seismic one-way wave separating snapshots, numerical records and the corresponding Poynting vector in Marmousi model

2.4 行波分離方法在逆時(shí)成像中的應(yīng)用

以圖2速度模型為例，震源采用最大頻率為60 Hz的雷克子波，從模型最左側(cè)5m處開始向右側(cè)激發(fā)，震源深度均為10m，炮間距為50m，共激發(fā)50炮，檢波器埋置深度為10m，道間距為5m，共501個(gè)檢波器。合成記錄長度為2 s，時(shí)間步長為0.5ms，滿足計(jì)算所需的穩(wěn)定性條件［7］。采用交錯(cuò)網(wǎng)格有限差分法，其差分精度分別為時(shí)間2階和空間10階，在邊界處采用20個(gè)網(wǎng)格點(diǎn)的PML吸收邊界條件來壓制邊界反射波，可保證有效模擬區(qū)域內(nèi)具有較高的信噪比和人工截?cái)嗵幘哂休^好的邊界吸收效果，然后采用相關(guān)法逆時(shí)成像條件對采集到的所有炮記錄進(jìn)行逆時(shí)成像處理。

圖5（a）為炮點(diǎn)和檢波點(diǎn)處上行波和下行波傳播方向示意圖，S代表炮點(diǎn)波場，R代表檢波點(diǎn)波場，原因在于100m深度位置檢波器僅能采集這2種方向行波波場。分析波場傳播機(jī)理可知，檢波點(diǎn)的下行波波場（逆向傳播）其實(shí)與炮點(diǎn)上行波波場的傳播路徑一致，檢波點(diǎn)的上行波波場與炮點(diǎn)下行波波場的傳播路徑一致，將這些波波場兩兩作互相關(guān)法逆時(shí)成像得到的是逆時(shí)偏移噪聲，而檢波點(diǎn)的下行波波場與炮點(diǎn)下行波波場的傳播路徑以及檢波點(diǎn)的上行波波場與炮點(diǎn)上行波波場的傳播路徑不同，因此將這些波波場兩兩作互相關(guān)法逆時(shí)成像得到的是真實(shí)的地層有效信號［16］。圖 5（b）為常規(guī)逆時(shí)偏移成像剖面，其成像結(jié)果為

式中：*代表相關(guān)成像。

圖5 方向行波示意圖及其波場分離方法在逆時(shí)成像中的應(yīng)用Fig.5 Directional one-way wave maps and its application in reverse-time migration

分析圖 5（b）和圖 5（c）可知，應(yīng)用方向行波波場分離后的逆時(shí)成像剖面，其偏移噪聲得到了有效壓制，信噪比更高。由此可見，基于方向行波的地震波場分離方法在實(shí)際地震資料復(fù)雜構(gòu)造成像處理中具有較好的應(yīng)用效果。

3 結(jié)論

（1）根據(jù)波印廷矢量在地震波場中的特點(diǎn)，將其應(yīng)用于地震方向行波的波場分離正演數(shù)值模擬中，實(shí)現(xiàn)對常規(guī)地震波場中的左行波、右行波、上行波和下行波波場進(jìn)行自動(dòng)識別與分離，而且該方法計(jì)算量小，計(jì)算過程簡單易行。

（2）通過3種不同模型的數(shù)值模擬，驗(yàn)證了該方法的準(zhǔn)確有效性，并且得到了完全分離的左行波、右行波、上行波和下行波4種方向行波的波場快照和數(shù)值模擬記錄。通過地震方向行波的波場分離處理，可進(jìn)一步提高對地震波場的傳播過程及規(guī)律的認(rèn)識，同時(shí)還能為實(shí)際地震資料和VSP資料的方法驗(yàn)證提供指導(dǎo)與幫助。

（3）通過地震方向行波波場分離數(shù)值模擬在逆時(shí)成像中的應(yīng)用，其成像效果得到顯著改善，偏移噪聲得到有效壓制，因此，該方法在實(shí)際地震資料成像等處理中具有重要的研究和應(yīng)用價(jià)值。

［1］ 陳可洋.高階彈性波波動(dòng)方程正演模擬及逆時(shí)偏移成像研究［D］.大慶：大慶石油學(xué)院，2009.

［2］ Dablain M A.The application of high-order differencing to the scalar waveequation［J］.Geophysics，1986，51（1）：54-66.

［3］ 張中杰，騰吉文，楊頂輝.聲波與彈性波數(shù)值模擬中的褶積微分算子法［J］.地震學(xué)報(bào)，1996，18（1）：63-69.

［4］ 程冰潔，李小凡.2.5維地震波場褶積微分算子法數(shù)值模擬［J］.地球物理學(xué)進(jìn)展，2008，23（4）：1099-1105.

［5］ 朱生旺，魏修成.波動(dòng)方程非規(guī)則網(wǎng)格任意階精度差分法正演［J］.石油地球物理勘探，2005，40（2）：149-153.

［6］ 陳可洋.邊界吸收中鑲邊法的評價(jià)［J］.中國科學(xué)院研究生院學(xué)報(bào)，2010，27（2）：170-175.

［7］ 廖振鵬，周正華，張艷紅.波動(dòng)數(shù)值模擬中透射邊界的穩(wěn)定實(shí)現(xiàn)［J］.地球物理學(xué)報(bào)，2002，45（4）：533-545.

［8］ 馬在田.地震成像技術(shù)有限差分法偏移［M］.北京：石油工業(yè)出版社，1989.

［9］ 劉洪林，陳可洋，楊微，等.高階交錯(cuò)網(wǎng)格有限差分法縱橫波波場分離數(shù)值模擬［J］.地球物理學(xué)進(jìn)展，2010，25（3）：877-884.

［10］ 陳可洋.一階速度-應(yīng)力Biot雙相各向同性介質(zhì)彈性波波場分離數(shù)值模擬［J］.計(jì)算物理，2011，28（3）：404-412.

［11］ YYoon K，Marfurt K J，Starr W. Challenges in reverse-time migration［C］. 74th SEG Annual International Meeting，Denver，Colorado，2004：1057-1060.

［12］ Sava P C，F(xiàn)ormel S. Angle-domain common-image gathers by wavefieldcontinuation methods［J］. Geophysics，2003，68（14）：1065- 1074.

［13］ Fomel S. Theory of 3-D angle gathers in wave-equation seismic imaging［J］. Journal of Petroleum Exploration and Production Technology， 2011，1：11-16.

［14］ Sava P，F(xiàn)ormel S. Time-shift imaging condition in seismic migration ［J］. Geophysics，2006，71（6）：209-217.

［15］ Xu Sheng，Zhang Yu，Tang Bing. 3D angle gathers from Reversetime migration［J］. Geophysics，2011，76（2）：77-92.

［16］ 陳可洋，吳清嶺，范興才，等.地震波疊前逆時(shí)偏移脈沖響應(yīng)研究與應(yīng)用［J］.石油物探，2013，52（2）：163-170.