三維似French模型高精度彈性波正演數值模擬及其波場特征

劉阿男,劉思彤,陳可洋,楊 微

(1.東北石油大學地球科學學院,黑龍江大慶 163318; 2.大慶油田有限責任公司勘探開發研究院,黑龍江大慶163712; 3.大慶油田有限責任公司第六采油廠,黑龍江大慶 163712)

三維似French模型高精度彈性波正演數值模擬及其波場特征

劉阿男1,劉思彤1,陳可洋2,楊 微3

(1.東北石油大學地球科學學院,黑龍江大慶 163318; 2.大慶油田有限責任公司勘探開發研究院,黑龍江大慶163712; 3.大慶油田有限責任公司第六采油廠,黑龍江大慶 163712)

基于三維一階雙曲型各向異性介質速度-應力彈性波動方程,給出各向異性介質向各向同性介質轉化的條件,采用高階交錯網格有限差分法和PML吸收邊界條件,實現高精度三維彈性波正演數值模擬.以均勻介質為例,研究彈性波在三維空間的波場特征;以三維似French模型為例,合成不同時刻的三分量波場快照和共炮點記錄,并對比分析三維各向同性介質和各向異性介質在該模型中的波場響應特征.數值結果驗證了方法的正確性和有效性,同時具有精度高和邊界吸收效果明顯的優點.

三維彈性介質;高精度正演數值模擬;PML吸收邊界條件;似French模型;高階交錯網格有限差分法

0 引言

彈性波正演數值模擬在多分量地震資料處理、觀測系統采集設計、處理參數試驗、解釋方案驗證等方面發揮越來越重要作用,是地震學和地震勘探學研究領域中的前沿課題之一,也是難度較大的課題之一[1].近些年,隨計算機計算能力的增強該技術得到飛躍發展.彈性波正演數值模擬技術大多是基于二維波場觀測研究地震勘探問題,原因在于二維彈性波場的計算速度快、算法易實現、剖面易顯示,而實際的地球介質是三維空間,二維波場觀測忽略三維構造體的側向效應等很多重要信息,無法準確、真實地刻畫三維空間的彈性波場響應特征.只有通過三維彈性波正演數值模擬,才能全面地認識彈性介質與彈性波場之間的耦合關系及彈性波場的空間變化特征,后者對于研究和認識地球介質的各向異性特征尤為重要,然而三維彈性波場正演問題研究困難,其原因在于:三維空間的計算變量較多;空間網格節點配置關系復雜;邊界吸收問題處理較難;剖面立體顯示困難;計算量和存儲量巨大,三維聲波正演數值模擬的算法要比彈性波正演數值模擬的算法更為簡單、易實現.

目前,對三維聲波正演問題已得到豐富的研究成果[2-7],如頻散問題[8-15](精度問題)、邊界吸收問題[16]、穩定性問題[17]等,而對三維彈性波正演問題的研究相對較少,且大多是基于均勻單層介質研究彈性波場的傳播特征(波場響應),如IgelH等[18]討論三維各向異性介質中彈性波數值模擬的有限差分實現方法;張文生等[19]對三維正交各向異性介質彈性波場進行有限差分模擬;裴正林[20-21]采用高階交錯網格有限差分法對三維各向異性介質中彈性波場進行高精度數值模擬;董良國等[17]提出三維任意各向異性介質彈性波動方程的穩定性條件;夏凡等[22-23]提出三維均勻介質彈性波動方程的吸收邊界條件及其高精度數值計算方法.

筆者借助Linux集群計算機的高速運算能力和大容量存儲能力,采用高階交錯網格有限差分法求解三維各向異性介質彈性波動方程,通過數值實例分析彈性波在三維空間經典似French速度模型中的三分量波場響應特征,給出三維各向異性介質向各向同性介質轉化的條件,同時使用任意方位波場切片技術(包括波場快照切片、模擬記錄切片等),實現三分量彈性波場全方位可視化顯示,研究三分量彈性波場的傳播規律,以此驗證方法的準確有效性.

1 基本理論

通常三維二階各向異性介質(VTI介質)彈性波動方程可表示為

式中:U、V、W分別為x、y、z方向質點的振動位移;C為三維彈性張量(矩陣組合形式C也可稱為物性矩陣,一般情況下沿主對角線方向是對稱正定的);ρ為彈性介質密度;t為時間.參數均為空間坐標的函數.

當各向異性介質對稱軸與觀測系統不一致時,需要進行觀測坐標旋轉或者物性矩陣旋轉,以達到合理觀測目的.此時需要用到Bond矩陣旋轉變換,對物性矩陣旋轉運算方程為C′=MCMT,其中M為正交正定矩陣,C′為變換后的物性矩陣.

根據Thomsen參數的定義,準縱波相速度Vp0=C33/ρ;準橫波相速度Vs0=C44/ρ;度量準縱波各向異性強度的參數ε=(C11-C33)/(2C33);度量準橫波各向異性強度的參數γ=(C66-C44)/(2C44);影響垂直于VTI介質對稱軸方向附加準縱波速度的參數δ=[(C13+C44)2-(C33-C44)2]/[2C33(C33-C44)].

為實現三維各向異性介質向各向同性介質轉化,借用各向同性介質中的泊松比ν概念,建立三維VTI介質中準縱波和準橫波相速度之間的關系式,即

同時可推得

為便于應用高階交錯網格有限差分法,將三維二階各向異性彈性波動方程(式(1))變換為一階雙曲型速度—應力彈性波動方程的微分形式,即

采用高階交錯網格有限差分法(時間2階、空間10階差分精度)對式(2)進行高精度差分離散(形成差分行式),對離散方程中的各變量進行交錯網格節點配置,可得到三維高精度彈性波正演遞推算子.在人工截斷邊界處,采用統一格式的最佳匹配層吸收邊界條件(吸收層厚度取為20個節點),消除或削弱截斷邊界處的反射波.要保證邊界計算過程穩定,還需滿足差分迭代所需的穩定性條件[17].

2 波場特征

2.1 均勻介質三維三分量

三維均勻各向同性介質模型尺寸為500m×500m×500m,3個方向空間網格長度為5m,采用純縱波震源(在3個方向正應力場中加入等量外力擾動)并置震源于模型中央,其最大頻率為40Hz(波形為雷克子波),介質的縱波速度為2000m/s;泊松比為0.25;密度為1g/cm3;時間步長取0.5ms.滿足穩定性條件,波場快照的記錄時間為0.125s.波場任意方向切片顯示方式的定義:(x=A,yoz)代表沿x方向深度A處且在yoz平面內的波場切片(其中y為圖片的橫坐標;z為圖片的縱坐標);(t=B,yoz)代表在時間B節點處且在yoz平面內的波場切片,其他情況依此類推.

三維任一分量彈性波場在3個相互垂直方向上的波場快照切片見圖1.由圖1可知:在純縱波震源激勵下,x方向質點速度場中縱波波形與yoz平面對稱;y方向質點速度場中縱波波形與xoz平面對稱;z方向質點速度場中縱波波形與yox平面對稱,且對稱面兩側極性相反.這驗證在均勻各向同性介質中,三維純縱波波場具有球對稱的特點.純橫波三分量彈性波場在任意方向的質點振動速度和與其相垂直的平面存在對稱關系,且極性相一致(見圖2,在x分量質點的振動速度處加入外力擾動).數值計算結果表明計算精度較高、頻散現象很弱.該三維三分量彈性波場特征分析方法同樣適合于表征三維各向異性介質準縱波和準橫波波場的特點,受各向異性參數的影響,準縱波和準橫波波形不再是球對稱關系,而是橢圓形甚至是具有不規則特征的對稱關系.

2.2 三維似French模型

圖1 0.125ms時刻縱波震源情況下三維各向同性介質彈性波波場快照

采用似French模型的速度模型,包含1個水平層界面、1個傾斜斷層、1個隆起構造、1個凹陷構造(見圖3).模型尺寸為500m×500m×500m,3個方向空間網格長度為5m,震源位于模型地表中央位置,其最大頻率為50Hz,每層密度均為1g/cm3,時間步長為0.5ms,滿足穩定性條件,接收時間為1s.對于三維各向同性介質模型,上層介質縱波速度為2000m/s,下層介質縱波速度為2500m/s,泊松比為0.25;對于三維各向異性介質模型,上層介質準縱波速度為2000m/s,下層介質準縱波速度為2500m/s,泊松比為0.25,ε為0.1,δ為0.2,γ為0.

圖2 0.125ms時刻混合源情況下三維各向同性介質彈性波波場快照(x分量)

圖3 三維似French模型(單位:m)

在0.20,0.25s時,三維各向同性介質和各向異性介質波場快照中的x、y、z方向的中間位置,且沿平行于yoz、xoz、xoy平面的波場快照切片(只給出三維垂直分量的波場快照)見圖4和圖5.考慮三維似French速度模型與三維各向同性介質彈性波場快照,經過對比容易識別直達P波、傾斜界面的反射PP波和反射PS波、水平層界面引起的反射PP波和反射PS波、隆起構造的反射PP波和反射PS波、凹陷構造引起的回轉反射PP波和反射PS波;彈性波傳播至下層介質時形成的透射PP波和透射PS波,以及彈性波在隆起構造和凹陷構造之間相互作用形成的復雜波(包括轉換波能量).與三維各向異性介質彈性波場快照對比,彈性波場更為復雜,波場快照中包含由各向異性特性引起的直達S波,模型中構造形成的反射SP波、反射SS波、透射SP波、透射SS波,以及這些波在構造之間相互作用形成的復雜波形,直達S波的波場能量比其他反射波能量要強.由圖4可知,在不同的切片中,凹陷構造和隆起構造的側向傳播效應較為明顯.三維各向同性介質和各向異性介質數值模擬記錄的在不同位置處的波場切片(見圖5),兩者最明顯的特征在于各向異性介質中存在直達S波及其轉換波,其波場情況明顯比各向同性介質更為復雜.由此可知,三維彈性波在三維空間中進行傳播,僅考慮二維空間的彈性波傳播不夠準確和全面.另外,地震波到達模型邊界時無邊界反射波形成,三維彈性波場的傳播過程中無明顯的頻散現象,說明數值計算方法具有較高精度.

圖4 三維彈性波場快照切片顯示

圖5 三維彈性波模擬記錄切片顯示

3 結論

(1)采用高階交錯網格有限差分法求解三維一階雙曲型各向異性介質速度—應力彈性波動方程及統一格式的PML吸收邊界條件方程構建三維高精度彈性波正演遞推算子,給出三維VTI介質與三維各向同性介質之間的轉換關系,實現了三維均勻介質和似French模型的彈性波高精度正演數值模擬.

(2)使用任意方位彈性波場切片技術實現三維三分量彈性波場全方位可視化,研究三分量彈性波場在三維空間中的波場特征及其傳播規律,縱波或準縱波具有平面對稱、極性相反的特點;橫波或準橫波具有平面對稱、極性相同的特點.

(3)三維彈性波場正演數值模擬技術可以更加全面地表征彈性波在全三維空間中的傳播過程,計算方法準確有效.

[1] 裴正林,牟永光.地震波傳播的數值模擬[J].地球物理學進展,2004,19(4):933-941.

[2] 劉金平,張向君,董建平.波動方程顯式短算子三維正演模擬[J].吉林大學學報:地球科學版,2004,34(1):138-141.

[3] 王潤秋,張明,費建博,等.精細積分法三維地震波正演模擬[J].勘探地球物理進展,2006,29(6):394-397.

[4] 何兵壽,魏修成,劉洋.三維波動方程的數值頻散關系及其疊前和疊后數值模擬[J].石油大學學報:自然科學版,2001,25(1):67-71.

[5] 熊曉軍,賀振華,黃德濟.三維波動方程正演及模型應用研究[J].石油物探,2005,44(6):554-556.

[6] 熊高君,張學工,賀振華,等.三維定位原理與三維反射波場模擬[J].礦物巖石,2002,22(3):93-97.

[7] 熊高君,賀振華,張毅祥.三維混合延拓一步法波動方程正演模擬[J].成都理工學院學報,1999,26(1):48-51.

[8] Alford R M,Kelly K R,Boore D M.Accuracy of finite-difference modeling of the acoustic wave equation[J].Geophysics,1974,39 (6):834-842.

[9] Marfurt KJ.Accuracy of finite difference and finite element modeling of the scalar and elastic equations[J].Geophysics,1984,49(5): 533-549.

[10] Dablain M A.The application of high-differencing to the scalar wave equation[J].Geophysics,1986,51(1):54-66.

[11] Fei T,Larner K.Elimination of numerical dispersion in finite-difference modeling and migration by flux-corrected transport[J].Geophysics,1995,60(6):1830-1842.

[12] 董良國,李培明.地震波傳播數值模擬中的頻散問題[J].天然氣工業,2004,24(6):53-56.

[13] 吳國忱,王華忠.波場模擬中的數值頻散分析與校正策略[J].地球物理學進展,2005,20(1):58-65.

[14] 寧剛,熊章強,陳持.波動方程有限差分正演模擬誤差來源分析[J].物探與化探,2008,32(2):203-206.

[15] 董良國,馬在田,曹景忠,等.一階彈性波方程交錯網格高階差分解法[J].地球物理學報,2000,43(3):411-419.

[16] 陳可洋,楊微.優化的三維地震波旁軸近似吸收邊界條件[J].勘探地球物理進展,2009,32(3):179-181,206.

[17] 董良國,馬在田,曹景忠.一階彈性波方程交錯網格高階差分解法穩定性研究[J].地球物理學報,2000,43(6):856-864.

[18] Igel H,Morap,Riolletb.Anisotropic wave propagation through finite-difference grids[J].Geophysics,1995,60(4):1203-1216.

[19] 張文生,宋海斌.三維正交各向異性介質三分量高精度有限差分正演模擬[J].石油地球物理勘探,2001,36(4):422-432.

[20] 裴正林.三維各向異性介質中彈性波方程交錯網格高階有限差分法數值模擬[J].石油大學學報:自然科學版,2004,28(5):23-29.

[21] 裴正林.三維各向同性介質彈性波方程交錯網格高階有限差分法模擬[J].石油物探,2005,44(4):308-315.

[22] Xia F,Dong L G,Ma Z T.Absorbing boundary conditions in 3D elastic-wave numerical modeling[J].Chinese J.Geophys,2004,47 (1):132-136.

[23] Xia F,Dong L G,Ma Z T.The numerical modeling of 3-D elastic wave equation using a high-order,staggered-grid,finite difference scheme[J].Applied Geophysics,2004,1(1):38-41.

High-precision elastic w ave forw ard numerical simulation and its w ave field characteristic analysis on 3D like French model/2010,34(4):45-50

LIU A-nan1,LIU Si-tong1,CHEN Ke-yang2,YANG Wei3
(1.Geoscience College,Northeast Petroleum University,Daqing,Heilongjiang163318,China;2. Ex ploration and Development Institute ofDaqing Oilf ield Company L td.,Daqing,Heilongjiang 163712,China;3.Oil Recovery Plant N o.6,Daqing Oilf ield Corp.L td.,Daqing,Heilongjiang 163712,China)

3D elastic wave forward numerical simulation is the key means for studying seismic wave’s wave field characteristic and its propagating rules in complex 3D spatial elastic wave medium.Based on 3D one-order hyperbolic anisotropic medium velocity-stress elastic wave equation,the paper gives out the transformation condition from anisotropic medium to isotropic medium,and uses high-order staggered-grid finite-difference method and PML absorbing boundary condition to achieve high-precision 3D elastic wave forward numerical simulation.Take an isotropic medium as an example to study the wave field characteristic in 3D space,and take 3D like French model as another model to obtain three-component snapshots with different time and common-shot record,compare and analyze the wave field response in this model of 3D isotropic medium and anisotropic medium.The numerical result verifies the method’s correctness and effectiveness,and shows the advantages of high-precision and obvious perfect boundary absorption effect.

3D elastic medium ; high-precision forward numerical simulation ; PML absorbing boundarycondition ; French like model ; high-order staggered-grid finite-difference met hod

book=4,ebook=364

TE19

A

1000-1891(2010)04-0045-06

2010-06-22;審稿人:王云專;編輯:任志平,張兆虹

國家“973”重大基礎研究項目(2001CB209104)

劉阿男(1973-),女,碩士,工程師,主要從事油氣成藏方面的研究.