本發(fā)明涉及勘探地球物理領(lǐng)域，特別是涉及到一種正交介質(zhì)地震波場模擬頻散壓制方法。
背景技術(shù)：
：地震學(xué)研究的對象是地震波及其傳播的地球介質(zhì)，實(shí)際地球介質(zhì)是一種非均勻、非完全彈性、各向異性、多相態(tài)的介質(zhì)。地震學(xué)的發(fā)展歷程正是由簡單均勻、完全彈性、各向同性、單相態(tài)的波動理論向復(fù)雜的真實(shí)地球介質(zhì)的波動理論步步逼近的過程(梁鍇,2009)。開展地震各向異性研究對認(rèn)知地球介質(zhì)結(jié)構(gòu)、勘探開發(fā)復(fù)雜油氣藏和預(yù)報(bào)地質(zhì)災(zāi)害等方面均具有理論意義和實(shí)用價(jià)值，是認(rèn)知地球介質(zhì)歷史發(fā)展的必然(吳國忱，2006)。隨著油氣勘探開發(fā)的深入，勘探的目標(biāo)逐步由構(gòu)造油氣藏轉(zhuǎn)化為裂縫性油氣藏，裂隙是許多油氣藏中液體或氣體的重要通道，正確的識別裂隙發(fā)育位置對于油氣藏的勘探和開發(fā)均具有重要意義(秦海旭，2015)。裂隙中彈性波傳播表現(xiàn)為方位各向異性介質(zhì)。具有方位各向異性特征的介質(zhì)包括HTI介質(zhì)、正交各向異性介質(zhì)和單斜各向異性介質(zhì)等。Tsvankin在2001年指出可將發(fā)育在橫向各向同性介質(zhì)中的多組平行排列直立裂隙等效為正交各向異性(Tsvankin,2001)，圖1給出了正交各向異性介質(zhì)示意圖(Tsvankin,2001)。更多情況下前人研究的是由廣泛擴(kuò)容各向異性(EDA)(Crampin，1987)和旋回性薄互層(PTL)(Postma，1995)等效而成(何燕，2008)。本文研究的是兩組正交發(fā)育于各向同性介質(zhì)中的直立裂隙等效成的正交各向異性介質(zhì)(Bakulin，2000partII)。各向異性介質(zhì)三大裂隙理論主要包括Hudson裂縫等效理論(Hudson,1981)、Thomsen裂隙等效理論(Thomsen,1995)和線性滑動理論(Schoenberg,1995)。這里利用線性滑動理論將前文提到的裂隙模型等效為正交各向異性來進(jìn)行巖石物理建模，為地震波數(shù)值模擬提供模型。地震波正演模擬是模擬地震波在地球介質(zhì)中的傳播過程，并研究地震波的傳播特性與地球介質(zhì)參數(shù)的關(guān)系(孫成禹,2007)，開展地震波的正演模擬研究，對人們正確認(rèn)識地震波的傳播規(guī)律，驗(yàn)證所求地球模型的正確性，進(jìn)行實(shí)際地震資料的地質(zhì)解釋與儲層預(yù)測等，均具有重要的理論和實(shí)際意義(吳國忱，2006)。各向異性介質(zhì)中地震波正演模擬的數(shù)值方法主要有：射線追蹤法、有限差分法、反射率法、偽譜法和有限元法等，這些數(shù)值正演模擬方法各有優(yōu)缺點(diǎn)。有限差分法是一種應(yīng)用比較廣泛的正演模擬方法，能夠較精確地模擬任意非均勻介質(zhì)中的地震波場，并含有多次散射、轉(zhuǎn)換波與繞射波。有別于規(guī)則網(wǎng)格有限差分法，交錯(cuò)網(wǎng)格有限差分法采用一階速度-應(yīng)力彈性波方程，其無須對彈性常數(shù)進(jìn)行空間微分(Virieux,1984)，在相同儲存空間和計(jì)算量的情況下具有更高的模擬精度(董良國，2000)。傳統(tǒng)交錯(cuò)網(wǎng)格有限差分空間導(dǎo)數(shù)的高階差分系數(shù)一般是通過Taylor級數(shù)展開法求取的，而用該方法確定的差分系數(shù)來求解空間導(dǎo)數(shù)時(shí)，一般只是在一個(gè)較小的頻率范圍內(nèi)才能取得比較高的精度(楊蕾，2014)。劉洋在2013年采用最小二乘法求解二階空間導(dǎo)數(shù)差分系數(shù)(YangLiu,2013)，楊蕾在2014年使用最小二乘法優(yōu)化一階空間導(dǎo)數(shù)差分系數(shù)(楊蕾，2014)。技術(shù)實(shí)現(xiàn)要素：本發(fā)明的目的是提供一種可以壓制波場模擬中的數(shù)值頻散現(xiàn)象，明顯提高模擬精度的正交介質(zhì)地震波場模擬頻散壓制方法。本發(fā)明的目的可通過如下技術(shù)措施來實(shí)現(xiàn)：正交介質(zhì)地震波場模擬頻散壓制方法，該正交介質(zhì)地震波場模擬頻散壓制方法包括：步驟1，利用巖石物理理論和裂隙等效理論將裂隙物性參數(shù)轉(zhuǎn)化為正交介質(zhì)參數(shù)，為正交介質(zhì)正演模擬提供模型；步驟2，根據(jù)地震波動力學(xué)理論結(jié)合裂隙介質(zhì)剛度矩陣推導(dǎo)正交介質(zhì)彈性波一階速度-應(yīng)力方程；步驟3，利用有限差分法求解正交介質(zhì)彈性波方程模擬地震波在裂隙介質(zhì)中的傳播過程；步驟4，利用最小二乘優(yōu)化方法優(yōu)化正交介質(zhì)交錯(cuò)網(wǎng)格差分系數(shù)從而壓制頻散提高模擬精度；步驟5，進(jìn)行模擬測試。本發(fā)明的目的還可通過如下技術(shù)措施來實(shí)現(xiàn)：在步驟1中，借助線性滑動理論給出了各向同性介質(zhì)背景下的兩組正交直立裂隙的彈性系數(shù)矩陣：C=C11C12C13C12C22C23C13C23C33C44C55C66=C~100C~2---(1)]]>其中C為彈性系數(shù)矩陣，C11,C22,C33,C44,C55,C66,C12,C13,C23為彈性系數(shù)矩陣元素，為子矩陣，其又可表示為：C~1=1d(λ+2μ)l1m3λl1m1λl1m2λl1m1(λ+2μ)l3m1λl2m1λl1m2λl2m1(λ+2μ)(l1m3-l4)---(2)]]>C~2=μ(1-ΔT2)000μ(1-ΔT1)000μ(1-ΔT1)(1-ΔT2)(1-ΔT1ΔT2)---(3)]]>其中λ和μ分別為裂隙所在背景介質(zhì)的拉梅參數(shù)，ΔN和ΔT分別是法向柔度和切向柔度，它們與裂隙充填物有關(guān)，變量的下標(biāo)N1,N2,T1,T2分別表示第一組和第二組裂隙法向和切向，d＝1-r2ΔN1ΔN2,r＝1-2g,g＝μ/(λ+2μ)，l1,l2,l3,l4,m1,m2,m3,d為過渡參數(shù)，其可以表示為：l1=1-ΔN1,l2=1-r-ΔN1,l3=1-r2ΔN1,l4=4r2g2ΔN1ΔN2m1=1-ΔN2,m2=1-rΔN2,m3=1-r2ΔN2---(4)]]>切向柔量和法向柔量的表達(dá)式，其中K為體積模量：ΔNi=(λ+2μ)KNi1+(λ+2μ)KNi,ΔTi=μKTi1+μKTi---(5)]]>當(dāng)裂隙滿足Hudson理論的假設(shè)時(shí)，可以從Hudson理論中給出ΔT,ΔN的表達(dá)式：(1)當(dāng)裂縫中填充較小體積模量和剪切模量的固體時(shí)：ΔN=43g[1-g+(k′+4μ′3)/(πdμ)]e,ΔT=163[3-2g+4μ′πdμ]e---(6)]]>(2)當(dāng)裂縫為干裂縫時(shí)：ΔN=43g(1-g)e,ΔT=163(3-2g)e---(7)]]>(3)當(dāng)裂縫中填充無粘滯流體時(shí)：ΔN=0,ΔT=163(3-2g)e---(8)]]>其中，k'和μ'分別為填充介質(zhì)的體積模量和剪切模量，e和d分別為裂縫密度和裂縫的縱橫比。在步驟2中，根據(jù)應(yīng)力應(yīng)變關(guān)系：τxxτyyτzzτyzτxzτxy=c11c12c1300c16c12c22c2300c26c13c23c3300c36000c44c450000c45c550c16c26c3600c66exxeyyezzeyzexzexy---(9)]]>其中τxx,τyy,τzz,τxy,τxz,τyz為應(yīng)力張量，exx,eyy,ezz,exy,exz,eyz為應(yīng)變張量，C11,C22,C33,C44,C55,C66,C12,C13,C23為彈性系數(shù)矩陣元素，應(yīng)力位移關(guān)系：ρ∂2∂t2uxuyuz=∂∂x000∂∂z∂∂y0∂∂y0∂∂z0∂∂x00∂∂z∂∂y∂∂x0τxxτyyτzzτyzτxzτxy+ρFxFyFz---(10)]]>其中ρ為密度，ux,uy,uz為位移變量，F(xiàn)x,Fy,Fz為震源變量，位移應(yīng)變關(guān)系：exxeyyezzeyzexzexy=∂∂x000∂∂y000∂∂z0∂∂z∂∂y∂∂z0∂∂x∂∂y∂∂x0uxuyuz---(11)]]>即可推導(dǎo)三維正交介質(zhì)彈性波一階速度-應(yīng)力方程：ρ∂vx∂t=∂τxx∂x+∂τxz∂z+∂τxy∂y+ρFxρ∂vy∂t=∂τyy∂y+∂τyz∂z+∂τxy∂x+ρFyρ∂vz∂t=∂τzz∂z+∂τyz∂y+∂τxz∂x+ρFz]]>∂τxx∂t=c11∂vx∂x+c12∂vy∂y+c13∂vz∂z+c16(∂vx∂y+∂vy∂x)∂τyy∂t=c12∂vx∂x+c22∂vy∂y+c23∂vz∂z+c26(∂vx∂y+∂vy∂x)∂τzz∂t=c13∂vx∂x+c23∂vy∂y+c33∂vz∂z+c36(∂vx∂y+∂vy∂x)∂τyz∂t=c44(∂vy∂z+∂vz∂y)+c45(∂vx∂z+∂vz∂x)∂τxz∂t=c45(∂vy∂z+∂vz∂y)+c55(∂vx∂z+∂vz∂x)∂τxy∂t=c16∂vx∂x+c26∂vy∂y+c36∂vz∂z+c66(∂vx∂y+∂vy∂x)]]>其中，vx,vy,vz為速度變量，t表示時(shí)間變量。在步驟3中，利用泰勒展開得到時(shí)間2階和空間2N階差分形式：時(shí)間2階差分：∂u(t)∂t=u(t+Δt)-u(t-Δt)2Δt---(12)]]>空間2N階差分：∂u∂x=1ΔxΣn=1Ncn[u(x+(2n-1)Δx2)-u(x-(2n-1)Δx2)]+O(Δx2N)---(13)]]>其中O(Δx2N)為泰勒展開誤差項(xiàng)，N階差分系數(shù)cm可以通過以下方程求?。篴113151...(2N-1)1133353...(2N-1)3153555...(2N-1)5...............1(2N-1)3(2N-1)5(2N-1)...(2N-1)(2N-1)c1c2c3...cN=100...0---(14)]]>在步驟4中，將最小二乘優(yōu)化差分系數(shù)的方法應(yīng)用于三維正交介質(zhì)模擬中，一階空間導(dǎo)數(shù)Taylor展開：∂u∂x≈1hΣm=1Mcm[u(x+mh-0.5h)-u(x-mh+0.5h)]---(15)]]>其中h為空間采樣間隔，cm為空間差分系數(shù)，M為差分階數(shù)，將平面波解u(x)＝u0eikx帶入，其中k為波數(shù)，結(jié)合歐拉公式e±ix＝cosx±isinx，化簡可得：kh2≈Σm=1Mamsin[kh(2m-1)/2]---(16)]]>定義過渡函數(shù)β和為：可得：定義誤差函數(shù)：其中E表示誤差，aM為空間差分系數(shù)，根據(jù)最小二乘基本原理求解：∂E∂cm=0---(20)]]>寫成矩陣形式即：其中cm為空間差分系數(shù)，M為差分階數(shù)，算法和分別為：本發(fā)明中的正交介質(zhì)地震波場模擬頻散壓制方法，針對正交介質(zhì)波場模擬過程中數(shù)值頻散現(xiàn)象，和地震源子波主頻較高或模擬模型速度較低時(shí)頻散現(xiàn)象尤為明顯的情況，借助線性滑動理論，將裂隙密度、縱橫比等不同裂隙物性參數(shù)等效為正交介質(zhì)剛度矩陣，提出兩組正交直立裂隙介質(zhì)建模方法。利用交錯(cuò)網(wǎng)格高階有限差分法實(shí)現(xiàn)三維正交介質(zhì)彈性波波場模擬，并提出正交介質(zhì)中基于最小二乘法優(yōu)化的高精度交錯(cuò)網(wǎng)格高階有限差分法頻散壓制方法。模擬結(jié)果表明介質(zhì)各向異性強(qiáng)度與裂隙物性直接相關(guān)，相比普通交錯(cuò)網(wǎng)格，最小二乘優(yōu)化方法可以壓制波場模擬中的數(shù)值頻散現(xiàn)象，明顯提高模擬精度。附圖說明圖1為正交介質(zhì)模型的示意圖；圖2為本發(fā)明的正交介質(zhì)地震波場模擬頻散壓制方法的一具體實(shí)施例的流程圖；圖3為本發(fā)明的一具體實(shí)施例中正交介質(zhì)變量定義網(wǎng)格的示意圖；圖4為本發(fā)明的一具體實(shí)施例中優(yōu)化前后不同差分階數(shù)數(shù)值頻散隨頻率變化圖；圖5為本發(fā)明的一具體實(shí)施例中交錯(cuò)網(wǎng)格空間6階差分均勻介質(zhì)模型250MS正演波場優(yōu)化效果的示意圖；圖6為本發(fā)明的一具體實(shí)施例中交錯(cuò)網(wǎng)格空間6階差分均勻介質(zhì)模型500MS正演波場優(yōu)化效果的示意圖；圖7為本發(fā)明的一具體實(shí)施例中交錯(cuò)網(wǎng)格空間8階差分均勻介質(zhì)模型500MS正演波場優(yōu)化效果的示意圖。具體實(shí)施方式為使本發(fā)明的上述和其他目的、特征和優(yōu)點(diǎn)能更明顯易懂，下文特舉出較佳實(shí)施例，并配合附圖所示，作詳細(xì)說明如下。如圖2所示，圖2為本發(fā)明的正交介質(zhì)地震波場模擬頻散壓制方法的流程圖。步驟101：利用巖石物理理論和裂隙等效理論將裂隙物性參數(shù)轉(zhuǎn)化為正交介質(zhì)參數(shù)，為正交介質(zhì)正演模擬提供模型；該方法首先需要利用巖石物理理論和裂隙介質(zhì)等效理論建立正交介質(zhì)模型。Bakulin借助線性滑動理論給出了各向同性介質(zhì)背景下的兩組正交直立裂隙的彈性系數(shù)矩陣：C=C11C12C13C12C22C23C13C23C33C44C55C66=C~100C~2---(1)]]>其中C為彈性系數(shù)矩陣，C11,C22,C33,C44,C55,C66,C12,C13,C23為彈性系數(shù)矩陣元素，為子矩陣，其又可表示為：C~1=1d(λ+2μ)l1m3λl1m1λl1m2λl1m1(λ+2μ)l3m1λl2m1λl1m2λl2m1(λ+2μ)(l1m3-l4)---(2)]]>C~2=μ(1-ΔT2)000μ(1-ΔT1)000μ(1-ΔT1)(1-ΔT2)(1-ΔT1ΔT2)---(3)]]>其中λ和μ分別為裂隙所在背景介質(zhì)的拉梅參數(shù)，ΔN和ΔT分別是法向柔度和切向柔度，它們與裂隙充填物有關(guān)，變量的下標(biāo)N1,N2,T1,T2分別表示第一組和第二組裂隙法向和切向，d＝1-r2ΔN1ΔN2,r＝1-2g,g＝μ/(λ+2μ)，l1,l2,l3,l4,m1,m2,m3,d為過渡參數(shù)，其可以表示為：l1=1-ΔN1,l2=1-r-ΔN1,l3=1-r2ΔN1,l4=4r2g2ΔN1ΔN2m1=1-ΔN2,m2=1-rΔN2,m3=1-r2ΔN2---(4)]]>切向柔量和法向柔量的表達(dá)式，其中K為體積模量：：ΔNi=(λ+2μ)KNi1+(λ+2μ)KNi,ΔTi=μKTi1+μKTi---(5)]]>當(dāng)裂隙滿足Hudson理論的假設(shè)時(shí)，可以從Hudson理論中給出ΔT,ΔN的表達(dá)式：(1)當(dāng)裂縫中填充較小體積模量和剪切模量的固體時(shí)：ΔN=43g[1-g+(k′+4μ′3)/(πdμ)]e,ΔT=163[3-2g+4μ′πdμ]e---(6)]]>(2)當(dāng)裂縫為干裂縫時(shí)：ΔN=43g(1-g)e,ΔT=163(3-2g)e---(7)]]>(3)當(dāng)裂縫中填充無粘滯流體時(shí)：ΔN=0,ΔT=163(3-2g)e---(8)]]>其中，k'和μ'分別為填充介質(zhì)的體積模量和剪切模量，e和d分別為裂縫密度和裂縫的縱橫比。步驟102：根據(jù)地震波動力學(xué)理論結(jié)合裂隙介質(zhì)剛度矩陣推導(dǎo)正交介質(zhì)彈性波一階速度-應(yīng)力方程；利用地震波動力學(xué)理論結(jié)合裂隙各向異性理論推導(dǎo)正交介質(zhì)介質(zhì)彈性波一階速度-應(yīng)力方程。τxxτyyτzzτyzτxzτxy=c11c12c1300c16c12c22c2300c26c13c23c3300c36000c44c450000c45c550c16c26c3600c66exxeyyezzeyzexzexy---(9)]]>其中τxx,τyy,τzz,τxy,τxz,τyz為應(yīng)力張量，exx,eyy,ezz,exy,exz,eyz為應(yīng)變張量。應(yīng)力位移關(guān)系：ρ∂2∂t2uxuyuz=∂∂x000∂∂z∂∂y0∂∂y0∂∂z0∂∂x00∂∂z∂∂y∂∂x0τxxτyyτzzτyzτxzτxy+ρFxFyFz---(10)]]>其中ρ為密度，ux,uy,uz為位移變量，F(xiàn)x,Fy,Fz為震源變量。位移應(yīng)變關(guān)系：exxeyyezzeyzexzexy=∂∂x000∂∂y000∂∂z0∂∂z∂∂y∂∂z0∂∂x∂∂y∂∂x0uxuyuz---(11)]]>即可推導(dǎo)三維正交介質(zhì)彈性波一階速度-應(yīng)力方程：ρ∂vx∂t=∂τxx∂x+∂τxz∂z+∂τxy∂y+ρFxρ∂vy∂t=∂τyy∂y+∂τyz∂z+∂τxy∂x+ρFyρ∂vz∂t=∂τzz∂z+∂τyz∂y+∂τxz∂x+ρFz]]>∂τxx∂t=c11∂vx∂x+c12∂vy∂y+c13∂vz∂z+c16(∂vx∂y+∂vy∂x)∂τyy∂t=c12∂vx∂x+c22∂vy∂y+c23∂vz∂z+c26(∂vx∂y+∂vy∂x)∂τzz∂t=c13∂vx∂x+c23∂vy∂y+c33∂vz∂z+c36(∂vx∂y+∂vy∂x)∂τyz∂t=c44(∂vy∂z+∂vz∂y)+c45(∂vx∂z+∂vz∂x)∂τxz∂t=c45(∂vy∂z+∂vz∂y)+c55(∂vx∂z+∂vz∂x)∂τxy∂t=c16∂vx∂x+c26∂vy∂y+c36∂vz∂z+c66(∂vx∂y+∂vy∂x)]]>其中，vx,vy,vz為速度變量，t表示時(shí)間變量。步驟103：利用有限差分法求解正交介質(zhì)彈性波方程模擬地震波在裂隙介質(zhì)中的傳播過程；通常情況多利用泰勒展開可以得到時(shí)間2階和空間2N階差分形式：時(shí)間2階差分：∂u(t)∂t=u(t+Δt)-u(t-Δt)2Δt---(12)]]>空間2N階差分：∂u∂x=1ΔxΣn=1Ncn[u(x+(2n-1)Δx2)-u(x-(2n-1)Δx2)]+O(Δx2N)---(13)]]>其中O(Δx2N)為泰勒展開誤差項(xiàng)，方程不同變量在三維有限差分網(wǎng)格中的位置如圖3和表1所示。其中差分系數(shù)可以通過以下方程求取，如表2所示：113151...(2N-1)1133353...(2N-1)3153555...(2N-1)5...............1(2N-1)3(2N-1)5(2N-1)...(2N-1)(2N-1)c1c2c3...cN=100...0---(14)]]>表1正交介質(zhì)彈性參數(shù)變量的空間位置表表2普通交錯(cuò)網(wǎng)格差分系數(shù)步驟104：利用最小二乘優(yōu)化方法優(yōu)化正交介質(zhì)交錯(cuò)網(wǎng)格差分系數(shù)從而壓制頻散提高模擬精度；將最小二乘優(yōu)化差分系數(shù)的方法應(yīng)用于三維正交介質(zhì)模擬中。一階空間導(dǎo)數(shù)Taylor展開：∂u∂x≈1hΣm=1Mcm[u(x+mh-0.5h)-u(x-mh+0.5h)]---(15)]]>其中h為空間采樣間隔，cm為空間差分系數(shù)，M為差分階數(shù)，將平面波解u(x)＝u0eikx帶入，其中k為波數(shù)，結(jié)合歐拉公式e±ix＝cosx±isinx，化簡可得：kh2≈Σm=1Mamsin[kh(2m-1)/2]---(16)]]>定義過渡函數(shù)β和為：可得：定義誤差函數(shù)：其中E表示誤差，aM為空間差分系數(shù)，根據(jù)最小二乘基本原理求解：∂E∂cm=0---(20)]]>寫成矩陣形式即：其中cm為空間差分系數(shù)，M為差分階數(shù)，算法和分別為：若給定d＝0.75可求解新的差分系數(shù)如表3所示：表3優(yōu)化交錯(cuò)網(wǎng)格差分系數(shù)步驟105：進(jìn)行模擬測試，模型測試說明最小二乘數(shù)值頻散壓制技術(shù)對頻散的壓制效果和模擬精度的提高效果。圖4展示了不同子波主頻、不同差分階數(shù)情況下優(yōu)化差分系數(shù)前后理論模擬精度對比，由圖我們可以看出不同情況下模擬差分精度均有所提高。圖5、6、7展示了不同差分階數(shù)和不同旅行時(shí)差分系數(shù)優(yōu)化前后正交介質(zhì)正演波場情況，從圖中可以看出最小二乘優(yōu)化差分方法可以有效提高差分精度壓制數(shù)值頻散。本發(fā)明借助線性滑動理論，將裂隙密度、縱橫比等不同裂隙物性參數(shù)等效為正交介質(zhì)剛度矩陣，提出兩組正交直立裂隙介質(zhì)建模方法。利用交錯(cuò)網(wǎng)格高階有限差分法實(shí)現(xiàn)三維正交介質(zhì)彈性波波場模擬，并提出正交介質(zhì)中基于最小二乘法優(yōu)化的高精度交錯(cuò)網(wǎng)格高階有限差分法頻散壓制方法,使用最小二乘方法優(yōu)化正交介質(zhì)一階速度-應(yīng)力方程空間導(dǎo)數(shù)的差分系數(shù)，推導(dǎo)更大頻率范圍內(nèi)可以獲得高精度模擬結(jié)果的差分系數(shù)。模擬結(jié)果表明介質(zhì)各向異性強(qiáng)度與裂隙物性直接相關(guān)，相比普通交錯(cuò)網(wǎng)格，最小二乘優(yōu)化方法可以壓制波場模擬中的數(shù)值頻散現(xiàn)象，明顯提高模擬精度，子波主頻較高和模型波速較低時(shí)改善效果更為突出。當(dāng)前第1頁1 2 3