陳可洋

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

0 引言

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

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

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

1 基本理論

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

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

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

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

同時(shí)，泊松比σ公式定義如下：

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

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

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

2 數(shù)值模擬算例

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

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

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

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

2.2 傾斜界面速度模型

為了驗(yàn)證文中方法在含不同傾角地質(zhì)情況下的準(zhǔn)確性和有效性，以含多個(gè)角度的傾斜界面速度模型（圖2）為例進(jìn)行分析。模型大小為2.5 km×1.0 km，縱橫向空間步長(zhǎng)均為5m；模型頂層Ⅰ左側(cè)部分和底層Ⅵ厚度均為250m，模型頂層Ⅰ右側(cè)部分最大厚度為750m。模型中間部分包含一組不同傾角的傾斜界面，傾斜界面自右向左的傾角依次為30°，45°，60°和 90°，對(duì)應(yīng)的彈性參數(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ā)，合成記錄的道長(zhǎng)為1 s，時(shí)間步長(zhǎng)為0.5ms，滿足計(jì)算所需的穩(wěn)定性條件。地面采集深度為100m，VSP采集位置距離模型左側(cè)100m處，波場(chǎng)快照記錄時(shí)間為0.2 s。

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

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

圖3 0.2 s時(shí)刻多分量各向異性介質(zhì)彈性波行波分離波場(chǎng)快照、模擬記錄及其波印廷矢量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模型

為了進(jìn)一步驗(yàn)證文中方法在更為復(fù)雜的彈性各向異性介質(zhì)模型中的適用性及其應(yīng)用效果，以Marmousi模型為例進(jìn)行了分析。模型總大小為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ā)，合成記錄的道長(zhǎng)為2s，時(shí)間步長(zhǎng)為0.2ms，滿足計(jì)算所需的穩(wěn)定性條件。地面采集深度為100m，波場(chǎng)快照記錄時(shí)間為0.6 s。

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

3 結(jié)論

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

（2）通過彈性波方向行波的波場(chǎng)分離數(shù)值模擬，進(jìn)一步提高了對(duì)彈性波波場(chǎng)傳播規(guī)律的認(rèn)識(shí)，同時(shí)還能夠?yàn)槎嗖ǘ喾至繌椥圆ㄙY料的方法驗(yàn)證提供指導(dǎo)和幫助。

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

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

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

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

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

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

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

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

［8］ 張會(huì)星，寧書年.彈性波動(dòng)方程疊前逆時(shí)偏移［J］.中國(guó)礦業(yè)大學(xué)學(xué)報(bào)，2002，31（5）：371-375.

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

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

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

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

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

［14］ 陳可洋.一階速度-應(yīng)力Biot雙相各向同性介質(zhì)彈性波波場(chǎng)分離數(shù)值模擬［J］.計(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í)偏移數(shù)值試驗(yàn)［J］.石油物探，2014，53（1）：8-16.

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

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

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

［20］ 陳可洋，吳清嶺，范興才，等.高階高密度三維多波多分量彈性波波場(chǎng)分離正演數(shù)值模擬［J］.油氣藏評(píng)價(jià)與開發(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］ 陳可洋.幾種地震觀測(cè)方式的逆時(shí)成像分析［J］.巖性油氣藏，2013，25（1）：95-101.

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

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