張繁昌 桑凱恒 路亞威
(中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院, 山東青島 266580)
天然裂縫型油氣藏近年來在勘探和生產(chǎn)中引起了較大的關(guān)注。天然裂縫控制儲層的滲透率,決定油氣的分布以及采收率,因此描述波在裂縫介質(zhì)中的傳播特征對于識別裂縫以及提高采收率具有重要意義。
疊前方位反演是目前最常用的識別裂縫方向、估計裂縫密度以及預(yù)測流體分布的方法[1-3],其中的反射系數(shù)方程由準(zhǔn)靜態(tài)等效裂縫模型所得,以HTI介質(zhì)反射系數(shù)公式[4]、傅里葉級數(shù)展開[5]為代表?,F(xiàn)階段最常用、最直觀的裂縫介質(zhì)等效理論模型分別為線性滑移模型[6]、Hudson[7]剛度擾動硬幣型裂縫模型和Thomson[8]各向異性參數(shù)裂縫模型,通過忽略裂縫微觀結(jié)構(gòu)以及孔縫間的流體交換模擬巖石彈性模量。上述方法在甚高頻率下模擬的飽和巖石彈性模量適用于超聲頻帶。對于低頻,一般先求得干燥裂縫介質(zhì)的等效模量,然后利用各向異性Gassmann方程進行流體替換[9]。雖然人們探討了單組、兩組及多組裂縫等效模型及裂縫間的相互影響[10-13],但仍不能準(zhǔn)確地描述波在裂縫介質(zhì)中的傳播特征。
波在介質(zhì)中傳播時發(fā)生的頻散和衰減是表征裂縫介質(zhì)彈性特征的重要屬性。除介質(zhì)本身固有的衰減機制外,波在介質(zhì)中傳播時引起的孔隙間、縫隙間以及孔縫間的流體流動也是一種重要的衰減模式。常規(guī)準(zhǔn)靜態(tài)裂縫模型只適用于弛豫和非弛豫狀態(tài),因而裂縫介質(zhì)的多尺度、全頻帶衰減特征研究成為近年的研究熱點。Chapman等[14-15]假定每個硬幣型裂縫都與固定數(shù)量的球形孔隙相連,提出微觀到介觀尺度頻變裂縫裂隙模型?;赪hite周期層狀等效模型[16],Brajanovski等[17-18]將裂縫等效為無限長平面薄層,借助等效傳播矩陣的求解方法,得到了裂縫的法向頻變彈性模量,進一步結(jié)合層狀介質(zhì)的Backus平均理論,分別得到高、低頻極限裂縫介質(zhì)的頻變彈性模量。Galvin等[19-21]將裂縫的影響視為各向同性背景介質(zhì)中的擾動,將裂縫視為薄硬幣形狀,利用多重散射定理,給出了隨機裂縫介質(zhì)的法向頻變彈性模量。Brajanovski模型和Galvin模型的衰減在低頻范圍與頻率成正比,在高頻范圍與頻率的平方根成正比,與Johnson[22]的認識一致。Gurevich等[23]將Brajanovski模型和Galvin模型統(tǒng)一,給出了全頻帶裂縫介質(zhì)彈性模量的統(tǒng)一表達式;Guo等[24]進一步將其拓展到有限厚度平面裂縫模型。另外,借助Krzikalla 等[25]提出的頻變張量近似式,Galvin等[26]進一步分析了硬幣型裂縫VTI介質(zhì)的相速度頻散和衰減,Yang等[27]和Kong等[28]研究了平面裂縫介質(zhì)模型的相速度頻散、衰減以及黏彈性AVO響應(yīng)。
本文利用離散積分方法給出了硬幣型裂縫法向彈性模量的數(shù)值模擬過程,并結(jié)合線性滑移理論,分別得到了高、低頻彈性模量,進一步分析了油、氣、水飽和情況下的頻散、衰減以及黏彈性反射響應(yīng)。
對于硬幣型裂縫介質(zhì)的頻變性質(zhì),Galvin等[21]提出了數(shù)值模擬的Fredholm方程
(1)
式中:x、y為積分變量;S(x)=(2/π)[sin(ax)-axcos(ax)]/x2;R(x,y)=sin[a(x-y)]/(x-y)-sin[a(x+y)]/(x+y),a為硬幣型裂縫長度;B為待求函數(shù);p0為流體壓力;T(y)為中間函數(shù)
(2)
硬幣型裂縫介質(zhì)在垂直方向的頻變彈性模量p33與B的關(guān)系為[21]
(3)
式中:ρ為介質(zhì)密度;ε為裂縫密度;k1為快縱波的波數(shù)。
式(1)、式(3)即為Galvin等[21]的硬幣型裂縫數(shù)值模擬表達式,利用常規(guī)積分方法求解難度較大。為此,下面深入分析數(shù)值求解方法。
由文獻[20]可知,除了直接由式(1)求解B外,還可以使用方程
(4)
求取。式中:θ(ζ)為待求函數(shù);Γ(ζ,ξ)為核函數(shù),滿足
(5)
式中J為Bessel函數(shù)。θ(ζ)與B(x)的關(guān)系為
(6)
利用式(4)~式(6)也可以計算B,但相比式(1),式(4)~式(6)在形式上更復(fù)雜。在實際求解中,發(fā)現(xiàn)式(6)可以避免奇異點的影響,而且求取極限的步驟可以合并,進而基于Gauss-Lobatto離散積分,通過數(shù)值求解得到頻變彈性模量。
首先,取Bessel函數(shù)的高階項
對核函數(shù)Г化簡
(7)
對于區(qū)間[0,a]的積分函數(shù)f(x)可以離散為
(8)
式中:xj∈[0,a]為求積節(jié)點;Aj為求積系數(shù)。式(1)可以離散為
(9)
令x=x1,…,xn, 并記θi=θ(xi),Γij=Γ(xi,xj),Fi=F(xi), 則式(9)可寫為
(10)
將上式表示為矩陣形式
(11)
可以看出,求解式(11)可以得到每一個節(jié)點的θi值。需要注意,Г函數(shù)是T(y)sin(ζy)sin(ξy)在區(qū)間(0,∞)積分的震蕩函數(shù),當(dāng)積分變量y的值較大后積分函數(shù)快速衰減。圖1為當(dāng)頻率為10、100、10000Hz時積分函數(shù)隨lgy的變化曲線。由圖可見,當(dāng)積分變量y>104后,積分函數(shù)趨于0,因此為加快計算速度,積分區(qū)間一般取(0,104)。
進一步
(12)
上式將求極限問題轉(zhuǎn)變?yōu)榍蠓e分。利用Gauss-Lobatto方法對上式離散化
(13)
只要將式(11)得到的每一個節(jié)點處的θ(xj)值乘以對應(yīng)的節(jié)點值和積分系數(shù),最后求和就可以得到極限值。
圖1 當(dāng)頻率為10(a)、100(b)、10000Hz(c)時積分函數(shù)隨lgy的變化曲線
Galvin等[21]只給出裂縫法向的彈性模量,其他方向的彈性模量可以通過Krzikalla 等[25]推導(dǎo)的頻變彈性模量關(guān)系得到
(14)
式中Cij為剛度系數(shù),下標(biāo)i,j=1,…,6,上標(biāo)“high”代表高頻,“l(fā)ow”代表低頻??梢?,除了已知法向彈性模量p33,還需知道高、低頻極限條件的其他各個彈性參數(shù)。
研究低頻條件的流體充填孔隙介質(zhì)的等效彈性模量時,一般利用Gassmann方程[29]對干巖石骨架進行流體置換分析。
對于干巖石的描述,Galvin和Hudson裂縫模型給出的結(jié)果相同。干硬幣型裂縫的法向弱度ΔN和切向弱度ΔT分別為[30]
(15)
(16)
結(jié)合線性滑移模型,將干巖石和裂縫參數(shù)代入Gassmann方程,可以得到其他彈性參數(shù)
(17)
式中λ、μ為干巖石的拉梅參數(shù)。
當(dāng)波在介質(zhì)中傳播時,會引起應(yīng)力變化。在高頻條件下慣性占據(jù)主導(dǎo),因此流體沒有足夠的時間在孔隙和縫隙間流動使介質(zhì)中的應(yīng)力達到平衡狀態(tài),此時流體流動產(chǎn)生的影響可以忽略,可將裂縫和孔隙視為獨立存在、互不影響。巖石的高頻彈性模量可模擬為飽和背景巖石與飽和裂縫彈性模量的疊加,流體飽和裂縫法向弱度ΔN由Hudson模型模擬[30]
(18)
式中:μf為流體的剪切模量(近似為零);c為裂縫厚度。
飽和裂縫的切向弱度等于干巖石裂縫的切向弱度ΔT(式16)。
將飽和背景巖石與飽和裂縫的彈性模量疊加,得到裂縫介質(zhì)的高頻模量
(19)
VTI介質(zhì)中復(fù)相速度VP的頻散表達式為[31]
(20)
其中
D={[(p11-p44)sin2θ-(p33-p44)cos2θ]2+
黏彈性VTI介質(zhì)縱、橫波反射、透射系數(shù)滿足下式[32]
(21)
式中:R為反射系數(shù);T為透射系數(shù);下標(biāo)P代表縱波,S代表橫波,1代表界面上方,2代表界面下方。其中
(22)
(23)
W=p55(γsx+βsz)
(24)
Z=βp13sx+γp33sz
(25)
(26)
式中sz為垂直方向的慢度,和水平方向慢度sx及剛度的關(guān)系為
(27)
(28)
(29)
首先利用離散積分法模擬硬幣型裂縫介質(zhì)飽含水的法向頻變彈性模量,VTI裂縫介質(zhì)參數(shù)如表1所示,流體參數(shù)如表2所示,其頻散與衰減如圖2所示。
表1 VTI裂縫介質(zhì)參數(shù)
由圖可見:①在弛豫狀態(tài),由于波長遠大于裂縫尺度,由波動引起的應(yīng)力變化在很短時間達到平衡,速度幾乎不變;在非弛豫狀態(tài),應(yīng)力變化加快,流體的慣性逐漸成為主導(dǎo),速度趨于穩(wěn)定;在過渡狀態(tài),因為巖石中的流體流動減弱,整體柔度減小,對應(yīng)剛度增加,因而隨頻率增加,速度不斷變大(圖2左)。②當(dāng)裂縫密度增加時,巖石柔度增加,速度變小,對應(yīng)的衰減整體變大,但衰減峰值頻率基本不變(圖2右)。③如果裂縫密度不變,裂縫長度由小變大,速度與衰減曲線的趨勢基本不變,只是整體左移。
圖2 不同裂縫密度(上)、裂縫長度(下)的硬幣型裂縫的頻散(左)與衰減(右)
表2 流體參數(shù)
圖3為不同頻率的含水介質(zhì)速度和衰減隨入射角的變化,圖4為不同入射角的含水介質(zhì)頻散和衰減。由圖可見:在垂直于裂縫方向(θ=0°)速度和衰減變化最大,在平行裂縫方向(θ=90°)趨于定值;隨入射角變大,各向異性程度減小,速度的頻散趨勢減緩,衰減峰值也逐漸減小。 圖5為含氣介質(zhì)的速度與能量隨入射角和頻率的變化。由圖可見:由于含氣介質(zhì)的密度較含水介質(zhì)小,速度整體變小(圖5左);由于氣體黏度較低,當(dāng)頻率較高時,在短時間內(nèi)通過流體交換達到平衡,因而弛豫和非弛豫之間的過渡區(qū)向高頻移動,振幅峰值也向高頻移動(圖5右)。
圖6為含油介質(zhì)的速度與能量隨入射角和頻率的變化。由圖可見:速度變化與飽水介質(zhì)相似(圖6左);由于油的黏度較高,在較低頻帶就達到非弛豫狀態(tài),因而峰值頻率向低頻移動(圖6右)。
圖7為含水、含油、含氣介質(zhì)頻變反射系數(shù)與流度。由圖可見:①由于上層波阻抗大于裂縫介質(zhì),垂直反射系數(shù)為負,因此在弛豫與非弛豫過渡區(qū),隨頻率增加,裂縫剛度增加、巖石的波阻抗增加,導(dǎo)致反射系數(shù)幅值減小(圖7左)。②當(dāng)裂縫介質(zhì)飽含氣時,在入射角較小范圍,反射系數(shù)幅值明顯大于含水、含油(圖7下左)。③流度的峰值頻率與能量的峰值頻率相對應(yīng),即含油時峰值流度最大(圖7中右),含氣時峰值流度最小(圖7下右);隨入射角增大,各向異性程度減弱,整體幅值減小。
圖3 不同頻率的含水介質(zhì)速度(左)和衰減(右)隨入射角的變化
圖4 不同入射角的含水介質(zhì)頻散(左)和衰減(右)
圖5 含氣介質(zhì)的速度(左)與能量(右)隨入射角和頻率的變化
圖6 含油介質(zhì)的速度(左)與能量(右)隨入射角和頻率的變化
本文通過巖石物理分析,研究了VTI硬幣型裂縫介質(zhì)頻散、衰減規(guī)律的數(shù)值模擬方法,并分析了不同流體飽和介質(zhì)的頻散、衰減現(xiàn)象,得到以下認識:
(1)通過離散積分和高階近似方法,將第二類Fredholm積分方程與多重散射方程的遠場極限值求解相統(tǒng)一,為數(shù)值模擬硬幣型裂縫模型的法向頻變彈性模量提供快速、準(zhǔn)確的方法;結(jié)合線性滑移理論和Hudson裂縫模型,得到了頻變彈性模量的解析式。
(2)分析VTI介質(zhì)相速度與黏彈性反射、透射系數(shù)方程可知,裂縫長度決定弛豫向非弛豫過渡區(qū)的頻帶,裂縫密度影響衰減峰值;當(dāng)裂縫介質(zhì)中飽含不同流體時,流體的黏度成為決定弛豫向非弛豫過渡的又一因素,飽氣時反射系數(shù)幅值大于飽油、飽水時,流度的變化趨勢與衰減的變化趨勢相似。