地震波動(dòng)方程方向行波波場(chǎng)分離正演數(shù)值模擬與逆時(shí)成像

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

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

0 引言

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

1 基本理論

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

式中：P為質(zhì)點(diǎn)的應(yīng)力波場(chǎ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為波場(chǎng)傳播時(shí)間，s；ρ為密度，kg/m3。

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

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

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

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

2.1均勻介質(zhì)模型

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

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

圖1 0.2 s時(shí)刻均勻介質(zhì)地震行波分離的波場(chǎng)快照及其波印廷矢量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，橫、縱向空間步長(zhǎng)均為5m。模型頂層①和底層⑦厚度分別為150m和250m，對(duì)應(yīng)的地震波速度分別為1 800m/s和3 500m/s。模型中間部分包含一組不同傾角的傾斜界面，而且傾斜界面自右向左的傾角依次為 30°，45°，60°和90°，模型中介質(zhì)②～⑥對(duì)應(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ā)，合成記錄長(zhǎng)度為1 s，時(shí)間步長(zhǎng)為0.5ms，滿足計(jì)算所需的穩(wěn)定性條件［7］。地面采集深度為100m，VSP采集位置距模型左側(cè)100m，波場(chǎng)快照記錄時(shí)間為0.2 s。

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

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

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

圖4 M armousi模型中地震方向行波分離的波場(chǎng)快照、模擬記錄及其波印廷矢量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è)檢波器。合成記錄長(zhǎng)度為2 s，時(shí)間步長(zhǎng)為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í)成像條件對(duì)采集到的所有炮記錄進(jìn)行逆時(shí)成像處理。

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

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

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

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

3 結(jié)論

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

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

（3）通過地震方向行波波場(chǎng)分離數(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維地震波場(chǎng)褶積微分算子法數(shù)值模擬［J］.地球物理學(xué)進(jìn)展，2008，23（4）：1099-1105.

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

［6］ 陳可洋.邊界吸收中鑲邊法的評(píng)價(jià)［J］.中國(guó)科學(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)格有限差分法縱橫波波場(chǎng)分離數(shù)值模擬［J］.地球物理學(xué)進(jìn)展，2010，25（3）：877-884.

［10］ 陳可洋.一階速度-應(yīng)力Biot雙相各向同性介質(zhì)彈性波波場(chǎng)分離數(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.