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

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

0 引言

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

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

1 基本理論

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

式中：P為質點的應力波場，無量綱；u和w分別為2個空間方向的質點振動速度，m/s；v為介質速度，m/s；f為震源函數，無量綱；xs和zs分別為震源在模型x和z空間方向的坐標，m；t為波場傳播時間，s；ρ為密度，kg/m3。

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

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

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

2 方向行波分離數值算例

2.1均勻介質模型

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

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

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

2.2 傾斜界面速度模型

圖2 復雜介質模型Fig.2 Complex velocity model

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

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

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

2.3 M armousi模型

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

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

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

2.4 行波分離方法在逆時成像中的應用

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

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

式中：*代表相關成像。

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

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

3 結論

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

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

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

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

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

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

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

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

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

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

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

［9］ 劉洪林，陳可洋，楊微，等.高階交錯網格有限差分法縱橫波波場分離數值模擬［J］.地球物理學進展，2010，25（3）：877-884.

［10］ 陳可洋.一階速度-應力Biot雙相各向同性介質彈性波波場分離數值模擬［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，Formel 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，Formel 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］ 陳可洋，吳清嶺，范興才，等.地震波疊前逆時偏移脈沖響應研究與應用［J］.石油物探，2013，52（2）：163-170.