各向異性彈性介質(zhì)方向行波波場分離正演數(shù)值模擬

陳可洋

（中國石油大慶油田有限責任公司勘探開發(fā)研究院，黑龍江大慶163712）

0 引言

多波多分量彈性波波動方程的地震波正演數(shù)值模擬技術(shù)一直是國內(nèi)外地球物理學界的熱點。隨著彈性波波動理論和計算機技術(shù)的不斷發(fā)展，自20世紀60年代以來這項技術(shù)便得到了飛速發(fā)展。彈性波正演數(shù)值模擬技術(shù)能夠基本保持彈性波的幾何學、運動學和動力學等特征，因此可以達到精確模擬彈性波波場傳播規(guī)律的目的［1］。目前，彈性波正演數(shù)值模擬技術(shù)的相關(guān)理論已基本成熟，并已形成了一系列經(jīng)典算法和研究成果。例如，在數(shù)值離散方面［2-5］，形成了有限差分法、有限元法和偽譜法等方法，數(shù)值頻散問題得到了有效壓制；在處理人為截斷邊界方面［6-8］，形成了旁軸近似吸收邊界、阻尼吸收邊界及最佳匹配層（PML）吸收邊界等各類吸收邊界條件，提高了計算結(jié)果的信噪比；在計算過程穩(wěn)定性方面［9-10］，形成了特征值分析法、矩陣分析法等平面波解的分析思路；在描述地球介質(zhì)方面［11-13］，形成了均勻各向同性介質(zhì)、各向異性介質(zhì)、黏滯介質(zhì)、雙相介質(zhì)、裂隙介質(zhì)及多相孔隙介質(zhì)等復雜的近似地球介質(zhì)的相關(guān)理論的假設。彈性波正演數(shù)值模擬技術(shù)的引入不僅大大降低了物理模擬分析和野外多分量地震勘探的成本，而且有效地指導了多分量地震資料的采集、處理和解釋，以及相關(guān)方法的驗證等，因此，在當前的勘探地震學領域發(fā)揮著至關(guān)重要的作用。

隨著地震波數(shù)值模擬技術(shù)的進步，如何在彈性波正演模擬過程中實現(xiàn)不同類型波場的自動分離逐漸得到廣泛重視。其中，最具代表性的為彈性波波場分離數(shù)值模擬方法，其次為基于散度算子計算純縱波和基于旋度算子計算純橫波的波場分離方法［14］，以及構(gòu)建等效彈性波方程的波場分離數(shù)值模擬方法［11，13］。 同時，研究的等效彈性介質(zhì)也已從目前的單相介質(zhì)延伸到了雙相孔隙介質(zhì)中，而方向行波分離的方法早已得到廣泛應用。例如，單程波疊前深度偏移方法需要對地震波波動方程進行算子分離處理，從而得到單程波偏移算子［15］；地震波數(shù)值模擬中的旁軸近似吸收邊界也是采用時域單程波算子解決某空間方向的邊界吸收問題［6］；還有上行波和下行波分離逆時成像條件可以對目的層有貢獻的波場進行選擇性成像［16］，從而有效壓制逆時偏移噪聲的形成。但是，目前尚未見到關(guān)于彈性波，特別是各向異性介質(zhì)中方向行波波場分離方面的相關(guān)報道。

在前人研究的基礎上，筆者以3個二維理論模型為例，在多波多分量各向異性介質(zhì)波場模擬中引入波印廷矢量，開展上行波、下行波、左行波和右行波4個方向行波波場分離的彈性波正演數(shù)值模擬實驗，以期提高對彈性波波場傳播規(guī)律的認識，并指導多分量彈性波地震資料的數(shù)字處理。

1 基本理論

根據(jù)彈性波波動力學的相關(guān)理論［1］，二維一階雙曲型各向異性介質(zhì)彈性波波動方程通常有如下的表達式［17］：

式中：vx和vz分別為質(zhì)點振動速度的水平分量和垂直分量（二維矢量 U={vx，vz}），m/s；τxx，τzz和 τxz分別為應力的3個分量（x和z空間方向的正應力和切應力），N；Cij為各向異性介質(zhì)彈性參數(shù)，109N/m3；ρ為介質(zhì)密度，kg/m3。式（1）是根據(jù)牛頓第二定律和胡克定律來構(gòu)建的［1］，其具體的數(shù)值離散方法、吸收邊界條件及其穩(wěn)定性條件可參見文獻［18-20］。

另外，Thomson 參數(shù)［21］符號定義如下：

式中：ε和δ均為Thomson參數(shù)，無量綱。

同時，泊松比σ公式定義如下：

在此基礎上［22-25］對波印廷矢量進行波場數(shù)值特征分析，并采用基于該矢量波場對原始多分量彈性波進行波場分離計算，最終實現(xiàn)了多分量各向異性介質(zhì)彈性波左行波、右行波、上行波和下行波的波場分離數(shù)值模擬，并成功地應用于多個數(shù)值模擬算例中。

式中：σ為泊松比，無量綱。

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

2 數(shù)值模擬算例

2.1 均勻彈性各向異性介質(zhì)模型

為了驗證文中方法的準確性，設計了均勻彈性各向異性介質(zhì)模型。模型總大小為1 km×1 km，縱橫向空間步長均為5m。采用雷克子波波形的漲縮震源，其最大頻率為40Hz，并置于模型中央位置處激發(fā)。采用的彈性參數(shù)C33為4×109N/m3，ε，δ和σ均為無量綱，其值分別為 0.2，0.5和 0.25，密度為103kg/m3，時間步長為0.5ms。滿足計算所需的穩(wěn)定性條件，波場快照記錄時間為0.2 s。

圖 1（a1）和圖 1（a2）分別為采用公式（1）得到的水平分量和垂直分量各向異性介質(zhì)彈性波波場快照。經(jīng)分析可知，波場快照中存在2種波場，傳播較快的是縱波，其波前面為橢圓形，傳播較慢的是由各向異性引起的橫波，其波前面形狀不規(guī)則。在水平分量波場快照中，其左右兩側(cè)的波場極性相反；在垂直分量波場快照中，上下兩側(cè)的波場極性也相反。圖1（b1）為水平方向的波印廷矢量波場快照，藍色部分其數(shù)值為正，紅色部分其數(shù)值為負，且正好沿著過震源的垂直線分為兩瓣，即劃分出左行波波場和右行波波場；圖1（b2）為垂直方向的波印廷矢量波場快照，正好沿著過震源的水平線分為兩瓣，即劃分出上行波波場和下行波波場。于是根據(jù)圖1（b1）和圖 1（b2）對圖 1（a1）和圖 1（a2）的多分量彈性波波場快照進行分離，最終得到了分離后的左行波波場［圖 1（c1）和圖 1（c2）］、右行波波場［圖 1（d1）和圖 1（d2）］、上行波波場［圖 1（e1）和圖 1（e2）］和下行波波場［圖 1（f1）和圖 1（f2）］，且分離結(jié)果準確可靠。圖1（g1）和圖1（g2）分別為正應力和剪應力的波場快照。在以震源為原點的4個象限中，正應力波場不存在極性反轉(zhuǎn)現(xiàn)象，而剪應力波場存在極性反轉(zhuǎn)現(xiàn)象，對角線方向極性相一致。綜合分析認為，文中方法準確實現(xiàn)了均勻、多波多分量彈性各向異性介質(zhì)中方向行波的波場分離處理。

圖1 0.2 s時刻均勻介質(zhì)彈性波行波波場分離的波場快照及其波印廷矢量Fig.1 Elastic one-waywave field separating snapshots in isotropicmedia and the corresponding Poynting vector at0.2 s

2.2 傾斜界面速度模型

為了驗證文中方法在含不同傾角地質(zhì)情況下的準確性和有效性，以含多個角度的傾斜界面速度模型（圖2）為例進行分析。模型大小為2.5 km×1.0 km，縱橫向空間步長均為5m；模型頂層Ⅰ左側(cè)部分和底層Ⅵ厚度均為250m，模型頂層Ⅰ右側(cè)部分最大厚度為750m。模型中間部分包含一組不同傾角的傾斜界面，傾斜界面自右向左的傾角依次為30°，45°，60°和 90°，對應的彈性參數(shù) C33依次是：Ⅰ層為 3.24×109N/m3，Ⅱ?qū)訛?5.29×109N/m3，Ⅲ層為6.76×109N/m3，Ⅳ層為 9×109N/m3，Ⅴ層為 7.29×109N/m3，Ⅵ層為 12.25×109N/m3，密度均為 103kg/m3，ε，δ和σ均為無量綱，分別為0.2，0.5和0.25。采用雷克子波波形的漲縮震源，其最大頻率為40Hz，在模型橫向1 km和縱向0.2 km處激發(fā)，合成記錄的道長為1 s，時間步長為0.5ms，滿足計算所需的穩(wěn)定性條件。地面采集深度為100m，VSP采集位置距離模型左側(cè)100m處，波場快照記錄時間為0.2 s。

圖2 各向異性彈性介質(zhì)速度模型Fig.2 Anisotropic elasticmedium velocitymodel

圖 3（a1）和圖 3（a2）分別為 0.2 s時刻水平分量和垂直分量各向異性介質(zhì)彈性波波場快照。圖3（b1）和圖3（b2）分別為水平方向和垂直方向的波印廷矢量波場快照，藍色部分波場其數(shù)值為正，紅色部分波場其數(shù)值為負。分析圖3可知，在含傾斜界面的復雜速度模型中，左行波、右行波、上行波和下行波波場的傳播特征較為復雜，但基于波印廷矢量可以對不同類型的行波波場加以自動識別和區(qū)分。圖 3（c）～（f）分別為根據(jù)圖 3（b1）和圖 3（b2）對圖 3（a1）和圖3（a2）的多分量彈性波波場快照進行的分離處理結(jié)果，最終得到了完全分離的水平分量和垂直分量的左行波［圖 3（c1）和圖 3（c2）］、右行波［圖3（d1）和圖 3（d2）］、上行波［圖 3（e1）和圖 3（e2）］和下行波［圖 3（f1）和圖 3（f2）］波場快照。 分析分離后的多分量彈性波波場快照可知，文中方法準確實現(xiàn)了較為復雜介質(zhì)情況下的上行波、下行波、左行波和右行波的波場分離處理。圖3（g1）和圖3（g2）分別為正應力和剪應力的波場快照。圖3（h1）和圖3（h2）分別為地面采集的水平分量和垂直分量彈性波模擬記錄。圖3（i）～（l）分別為對應的水平分量和垂直分量左行波［圖 3（i1）和圖 3（i2）］與右行波［圖 3（j1）和圖 3（j2）］。 圖 3（k1）和圖 3（k2）分別為正應力和剪應力的彈性波模擬記錄。圖3（l1）和圖3（l2）分別為VSP采集的水平分量和垂直分量彈性波數(shù)值模擬記錄。圖3（m）～（n）分別為對應的水平分量和垂直分量上行波［圖 3（m1）和圖 3（m2）］與下行波［圖 3（n1）和圖 3（n2）］。 圖 3（o1）和圖 3（o2）分別為VSP采集的正應力和剪應力的數(shù)值模擬記錄。綜上所述，文中方法在較為復雜的多分量各向異性介質(zhì)彈性波數(shù)值模擬過程中準確實現(xiàn)了上行波、下行波、左行波和右行波的波場分離處理。

圖3 0.2 s時刻多分量各向異性介質(zhì)彈性波行波分離波場快照、模擬記錄及其波印廷矢量Fig.3 Elastic one-way wave separating snapshots and numerical records in multi-components anisotropic mediaand the corresponding Poynting vector at 0.2 s

2.3 M armousi模型

為了進一步驗證文中方法在更為復雜的彈性各向異性介質(zhì)模型中的適用性及其應用效果，以Marmousi模型為例進行了分析。模型總大小為3.4 km×1.4 km，縱橫向空間網(wǎng)格均為5m；最小彈性參數(shù) C33為 1.06×109N/m3，最大彈性參數(shù) C33為21.81×109N/m3，密度均為 103kg/m3，ε，δ和 σ 均為無量綱，分別為0.2，0.5和0.25。采用雷克子波波形的漲縮震源，其最大頻率為40Hz，在橫向1.7 km和縱向0.15 km處激發(fā)，合成記錄的道長為2s，時間步長為0.2ms，滿足計算所需的穩(wěn)定性條件。地面采集深度為100m，波場快照記錄時間為0.6 s。

圖4 0.6 s時刻M armousi模型多分量各向異性介質(zhì)彈性波行波分離波場快照、模擬記錄及其波印廷矢量Fig.4 Elastic one-way wave separating snapshots and numerical records in multi-components anisotropic mediaMarmousi model and the corresponding Poynting vector at 0.6 s

圖 4（a1）和圖 4（a2）分別為 0.6 s時刻水平分量和垂直分量各向異性介質(zhì)彈性波波場快照。圖4（b1）和圖4（b2）分別為水平方向和垂直方向的波印廷矢量波場快照，藍色部分波場其數(shù)值為正，紅色部分波場其數(shù)值為負。分析圖4可知，在基于Marmousi速度模型的波場快照中，左行波、右行波、上行波和下行波的波波場傳播特征更為復雜，存在各種類型的彈性波場轉(zhuǎn)換關(guān)系，但基于波印廷矢量可對不同類型的行波波場加以識別和分離。圖4（c）～（f）分別為根據(jù)圖 4（b1）和圖 4（b2）對圖 4（a1）和圖4（a2）的多分量彈性波波場快照進行的分離處理結(jié)果，最終得到了完全分離的水平分量和垂直分量左行波［圖 4（c1）和圖 4（c2）］、右行波［圖 4（d1）和圖 4（d2）］、上行波［圖 4（e1）和圖 4（e2）］與下行波［圖 4（f1）和圖 4（f2）］波場快照。 分析分離后的多分量各向異性介質(zhì)彈性波波場快照可知，文中方法也準確實現(xiàn)了上行波、下行波、左行波和右行波的波場分離處理。圖4（g1）和圖4（g2）分別為正應力和剪應力的波場快照。圖4（h1）和圖4（h2）分別為地面采集的水平分量和垂直分量彈性波模擬記錄。圖4（i）～（j）分別為對應的水平分量和垂直分量左行波［圖 4（i1）和圖 4（i2）］與右行波［圖 4（j1）和圖 4（j2）］波場的數(shù)值模擬記錄。 ［圖 4（k1）和圖 4（k2）］分別為地面采集的正應力和剪應力的數(shù)值模擬記錄。綜合分析可知，本文方法準確實現(xiàn)了復雜彈性各向異性介質(zhì)上行波、下行波、左行波和右行波的波場分離處理。

3 結(jié)論

（1）基于多分量彈性波波印廷矢量，實現(xiàn)了多分量各向異性介質(zhì)彈性波上行波、下行波、左行波和右行波的方向行波波場分離數(shù)值模擬。該方法計算量小，計算過程簡單，易實現(xiàn)。通過3個不同的理論模型與數(shù)值模擬，驗證了該方法的準確性和有效性。

（2）通過彈性波方向行波的波場分離數(shù)值模擬，進一步提高了對彈性波波場傳播規(guī)律的認識，同時還能夠為多波多分量彈性波資料的方法驗證提供指導和幫助。

（3）由于多波多分量彈性波逆時偏移理論仍處于理論研究和探索階段，該方法目前尚未在逆時成像方面進行應用，但這種方向行波分離的思路構(gòu)建的優(yōu)化逆時成像條件在聲波波動方程中已取得了較好的應用效果，因此，該方法在多波多分量地震資料處理中具有重要的應用價值。

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

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

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

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

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

［6］ 張秉銘.各向異性介質(zhì)中彈性波數(shù)值模擬與偏移研究［D］.北京：中國科學院地球物理研究所，1997.

［7］ 楊微，陳可洋.加權(quán)吸收邊界條件的優(yōu)化設計［J］.石油物探，2009，48（3）：244-246.

［8］ 張會星，寧書年.彈性波動方程疊前逆時偏移［J］.中國礦業(yè)大學學報，2002，31（5）：371-375.

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

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

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

［12］ 寧剛，熊章強，陳持遜.波動方程有限差分正演模擬誤差來源分析［J］.物探與化探，2008，32（2）：203-206.

［13］ 陳可洋.高精度地震純波震源數(shù)值模擬［J］.巖性油氣藏，2012，24（1）：84-91.

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

［15］ DenliH，Huang Lianjie.Elastic-wave reverse-timemigrationwith a wavefield-separation imaging condition［C］.SEG Las Vegas2008 AnnualMeeting，2008：2346-2350.

［16］ 陳可洋，陳樹民，李來林，等.彈性波聯(lián)合疊前逆時偏移數(shù)值試驗［J］.石油物探，2014，53（1）：8-16.

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

［18］ 裴正林，牟永光.地震波傳播的數(shù)值模擬［J］.地球物理學進展，2004，19（4）：933-941.

［19］ 陳可洋.地震波數(shù)值模擬中差分近似的各向異性分析［J］.石油物探，2010，49（1）：19-22.

［20］ 陳可洋，吳清嶺，范興才，等.高階高密度三維多波多分量彈性波波場分離正演數(shù)值模擬［J］.油氣藏評價與開發(fā)，2013，3（2）：6-14.

［21］ Thomsen L.Weak elastic anisotropy［J］. Geophysics，1986，51（10）：1954-1966.

［22］ Yoon K，Marfurt Kurt J，Starr W. Challenges in reverse-time migration［G］. SEG 74th Annual International Meeting Expanded Abstracts，2003：1057-1060.

［23］ 陳可洋.幾種地震觀測方式的逆時成像分析［J］.巖性油氣藏，2013，25（1）：95-101.

［24］ 陳可洋，吳沛熹，楊微.擴散濾波方法在地震資料處理中的應用研究［J］.巖性油氣藏，2014，26（1）：117-122.

［25］ 陳可洋，楊微，吳清嶺，等.地震反射波與散射波波場分離方法初探［J］.巖性油氣藏，2013，25（2）：76-81.