裂縫參數(shù)對縱波各向異性影響的數(shù)值模擬

段 茜 劉向君 梁利喜 熊 健

(①油氣藏地質(zhì)及開發(fā)工程國家重點實驗室(西南石油大學(xué))，四川成都 610500； ②西南石油大學(xué)理學(xué)院，四川成都 610500)

0 引言

裂縫描述是油藏開發(fā)的重要內(nèi)容，準(zhǔn)確預(yù)測裂縫發(fā)育帶的方向以及最大裂縫密度區(qū)域是勘探地球物理的一項重大技術(shù)難題。理論與實踐均已證實，當(dāng)?shù)叵陆橘|(zhì)中存在裂縫時，縱波響應(yīng)呈方位各向異性特征，即沿不同方位具有不同的動力學(xué)和運動學(xué)參數(shù)。因此，可根據(jù)裂縫型介質(zhì)中的縱波響應(yīng)方位各向異性預(yù)測裂縫參數(shù)。

自20世紀(jì)90年代以來，人們利用不同方位的速度、旅行時、振幅等多種屬性預(yù)測裂縫。Tsvanskin[1]研究了HTI介質(zhì)中縱波動校正速度隨方位的變化規(guī)律，提出采用縱波動校正速度估算橫波分裂系數(shù)，進而預(yù)測裂縫密度。Ruger[2]研究了HTI介質(zhì)中縱波反射系數(shù)隨炮檢距和方位角的變化，認(rèn)為振幅與AVO梯度隨方位角余弦函數(shù)的平方而變化，給出了反射振幅與各向異性系數(shù)之間的簡單關(guān)系。Perez等[3]采用多種方法預(yù)測委內(nèi)瑞拉西南部碳酸鹽巖儲層的裂縫走向。董淵等[4]設(shè)計了一個簡單的三層水平界面模型，采用交錯網(wǎng)格有限差分方法，在以45°夾角彼此相交的4條二維測線的縱波CDP道集上，通過拾取層間時差求取各向異性參數(shù)，并預(yù)測地下裂縫發(fā)育情況。Shen等[5]利用多信號分類技術(shù)，從二維縱波地震數(shù)據(jù)中提取振幅和頻率隨炮檢距的變化屬性檢測碳酸鹽巖儲層的裂縫方向。Li等[6]嘗試?yán)萌S縱波地震資料，在中國東部黃河三角洲地區(qū)檢測裂縫。張立勤等[7]將地震數(shù)據(jù)劃分為5個方位角數(shù)據(jù)，按不同的方位角范圍對地震道進行振幅補償、校正和疊前時間偏移，產(chǎn)生不同方位角的處理結(jié)果，從而預(yù)測儲層的各向異性特征。楊勤勇等[8]通過建立裂縫物理模型，得到了速度與振幅隨測線方位變化的關(guān)系曲線，結(jié)果表明，當(dāng)測線方位與裂縫方向平行時反射時間最小、振幅最大，隨著測線方位與裂縫走向間的夾角增大，反射時間逐漸增大、振幅逐漸減小，直至觀測方位與裂縫方向垂直時，反射時間最大、振幅最小。Zhu等[9]研究了橫觀各向同性樣品的P波衰減系數(shù)，認(rèn)為衰減各向異性遠遠超過速度各向異性，從而為裂縫探測和巖性識別提供了敏感屬性。齊宇等[10]運用物理模型模擬方法，在水中以50、100、150mm的炮檢距觀測了HTI介質(zhì)的地震響應(yīng)，模擬了裂縫介質(zhì)的反射波振幅、反射旅行時隨裂縫方位角變化的規(guī)律。王月英等[11]分析了兩組任意夾角的傾斜裂縫介質(zhì)模型的地震波場數(shù)值模擬結(jié)果，認(rèn)為在發(fā)育強度相同的情況下，傾斜裂縫較垂直裂縫的方位各向異性程度弱，裂縫發(fā)育趨于復(fù)雜時，各向異性特征趨于減弱。王洪求等[12]分析和對比了旅行時、旅行時差、振幅、AVO梯度的方位各向異性特征及其裂縫預(yù)測精度，在此基礎(chǔ)上優(yōu)選有利屬性進行融合，進一步提高了裂縫預(yù)測精度。王玲玲等[13]基于制作的多參數(shù)復(fù)雜裂縫儲層地震物理模型，利用相干、傾角、振幅、方差、曲率、衰減、頻譜分解等屬性識別裂縫及優(yōu)選敏感屬性。

盡管基于縱波屬性方位各向異性的裂縫預(yù)測方法較多，但在巖石物理實驗中，人工裂縫模型類型有限且裂縫多為均勻分布，研究結(jié)果具有一定的局限性。裂隙介質(zhì)等效理論在裂縫形狀、裂縫間的相互作用、流體作用等方面大多進行了近似或簡化處理，難以精細描繪復(fù)雜縫洞系統(tǒng)及其流體分布，計算結(jié)果通常與實驗結(jié)果存在較大差異。利用地震資料研究天然裂縫的縱波方位各向異性，目前還缺乏縱波屬性對裂縫的敏感性系統(tǒng)分析。為此，本文基于聲波波動理論，對隨機離散裂縫模型的聲波波場進行數(shù)值模擬，通過計算不同測線方位的聲波速度和衰減系數(shù)，分析聲學(xué)參數(shù)與裂縫分布、走向、密度及流體間的關(guān)系，以期為裂縫型儲層參數(shù)預(yù)測提供一定的理論基礎(chǔ)。

1 聲學(xué)參數(shù)的計算

1.1 計算聲波速度

本文對超聲波透射實驗進行數(shù)值模擬，以發(fā)射探頭的激發(fā)信號為震源，模擬其通過不同氣/水兩相裂縫型巖樣后接收探頭的波形圖，再拾取接收端的波形初至得到波的傳播時間，利用

(1)

即可計算聲波速度VP。式中：L為巖樣長度；TP為縱波傳播時間。

1.2 計算聲波衰減系數(shù)

采用標(biāo)準(zhǔn)樣品對比法，將鋁塊作為標(biāo)準(zhǔn)樣品，再取一塊與鋁塊長度相同的巖樣作為待測樣品，分別提取鋁塊和巖樣的聲波首波振幅值，即可利用

(2)

計算巖樣的聲波衰減系數(shù)α。式中：A0和A分別為鋁塊和巖樣的聲波首波振幅；α0為鋁塊的聲波衰減系數(shù)(實際計算時可近似為0)。

2 裂縫模型及觀測系統(tǒng)

2.1 單裂縫模型聲學(xué)參數(shù)試算

基于二階聲波波動方程，對氣/水兩相裂縫型巖樣的聲波波場進行有限差分?jǐn)?shù)值模擬。圖1、圖2分別為鋁塊、單裂縫巖樣聲波波場快照，圖3為鋁塊和巖樣接收端波形圖。由圖3可見：從接收端提取的巖樣首波初至?xí)r間為8.38 μs，由式(1)得到VP為5966.6m/s；從波形圖上提取鋁塊和巖樣的聲波首波振幅分別為1.0和0.1644，代入式(2)即可計算出該巖樣的α為36.11dB/m。

2.2 隨機離散裂縫模型的建立

前人已討論了各類連續(xù)隨機介質(zhì)模型的構(gòu)造方法[14-17]，但這些方法幾乎都使用一個均值為零的平穩(wěn)隨機過程表示介質(zhì)在空間小尺度的非均勻性，這并不適用于描述分布極不規(guī)則的裂縫型油氣藏。在地殼介質(zhì)中裂縫都呈離散且隨機分布，裂縫的存在使介質(zhì)的速度和密度發(fā)生不連續(xù)變化。因此，應(yīng)將裂縫型介質(zhì)看作裂縫和基質(zhì)的結(jié)合體，采用離散裂縫模型表征。離散裂縫模型可以表征裂縫性油藏在任意尺度的非均質(zhì)性，可清楚地顯示介質(zhì)中的每條裂縫，具有計算精度高、擬真性好的優(yōu)點?；跀?shù)字圖像處理技術(shù)，通過設(shè)置孔隙度φ、裂縫密度d、裂縫長寬比r及裂縫傾角θ(裂縫與水平線間的夾角)等裂縫參數(shù)，采用鄰點融合方法建立非均勻性的隨機裂縫介質(zhì)，構(gòu)建步驟如下。

圖1 鋁塊聲波波場快照(t=8μs)設(shè)鋁塊尺寸為50mm(長度)×50mm(直徑)，時間采樣步長為10ns，將鋁塊的縱向剖面區(qū)域劃分成250×250的網(wǎng)格，因此空間采樣步長為0.2mm。震源為主頻100kHz的雷克子波，位于(25mm，0)處，接收探頭位于(25mm，50mm)處，采用吸收邊界條件[18]。鋁塊縱波速度為6350m/s?！锎碚鹪次恢?，▲代表接收探頭位置，下同

圖2 單裂縫巖樣聲波波場快照(t=8μs)巖樣的尺寸、采樣步長、縱向剖面區(qū)域劃分、震源、邊界條件等參數(shù)同圖1。裂縫位于模型中心，縫長為14mm，縫寬為3mm，裂縫與水平線之間的夾角為0°。巖樣骨架的縱波速度為6200m/s，裂縫孔隙流體為氣/水兩相，氣和水的密度分別為0.72、1000kg/m3，體積模量分別為0.1、2.1GPa。設(shè)孔隙的含水飽和度為50%，在流體替換過程中，按照Wood模型[19]計算出孔隙流體的密度為500.36kg/m3，體積模量為0.19GPa，因此裂縫孔隙流體的聲波速度為616.22m/s

圖3 鋁塊和巖樣接收端波形圖

(1)選定裂縫模型參數(shù)，包括φ、d、r和θ。

(2)根據(jù)設(shè)置的φ和d，計算給定模型范圍內(nèi)的裂縫的面積s和裂縫條數(shù)n。

(3)隨機生成n個像素點作為n條裂縫的中心Mn(xn,zn)，按如下步驟逐個融合其鄰域內(nèi)像素，直到模型區(qū)域內(nèi)的所有裂縫面積達到s為止。

1)裂縫的像素點構(gòu)成的圖像子集記為P，將P的鄰點構(gòu)成的圖像子集記為N′(P)，其中N′(P)與P無公共像素點。初始時，每條裂縫只包含一個像素點，即P={Mn}。

2)由r確定P中像素點的4個鄰域方向(上、下、左、右)的融合概率cj。

3)按照cj隨機將N′(P)的像素點M′加入P中P=P+{M′},M′∈N′(P)。

4)若P的面積達到s，則將圖像旋轉(zhuǎn)角度θ，否則轉(zhuǎn)步驟3)。

通過選擇4個裂縫參數(shù)，可按以上步驟構(gòu)造不同形式的隨機裂縫介質(zhì)模型。

2.3 實驗觀測系統(tǒng)

圖4為縱波方位各向異性特征的實驗觀測系統(tǒng)。

圖4 定炮檢距觀測系統(tǒng)示意圖發(fā)射探頭與接收探頭間的距離固定，并將兩探頭固定在巖心頂部某一水平面上同時繞中心點逆時針旋轉(zhuǎn)，觀測并記錄透過模型的超聲波速度與衰減特性。測線方位每旋轉(zhuǎn)15°采集1道數(shù)據(jù)，當(dāng)測線方位旋轉(zhuǎn)360°(共采集25道數(shù)據(jù))時采集1條測線。測線與裂縫走向平行時為0°或180°方位，測線與裂縫走向垂直時為90°或270°方位，以測線與裂縫走向平行時開始觀測

3 模擬結(jié)果與分析

3.1 裂縫傾角及密度對方位各向異性特征的影響

建立20個隨機離散裂縫模型(圖5)。由于采用360°觀測系統(tǒng)，在研究裂縫方位各向異性特征時將裂縫走向設(shè)置成某一固定方位。圖6為測線方位與裂縫走向夾角(下文簡稱測線方位角)不同時模型8的聲波位移波場快照。圖7～圖11分別為模型1～4、5～8、9～12、13～16、17～20的各向異性特征曲線。由模擬結(jié)果可見，VP和α隨測線方位呈各向異性特征，整個變化以180°為周期，隨著測線方位角增大，VP減小(圖7a、圖8a、圖9a、圖10a、圖11a)、α增大(圖7b、圖8b、圖9b、圖10b、圖11b)。如：當(dāng)測線方位角為0°或180°時，VP較大、α較??；當(dāng)測線方位角為90°或270°時，VP較小、α較大。

圖5 隨機離散裂縫模型巖心尺寸為50mm(長度)×50mm(直徑)，時間采樣步長為10ns，空間采樣步長為0.2mm。φ=4%，θ=65°，r均值為10，氣/水兩相裂縫孔隙流體的含水飽和度為50%。從左至右、從上至下依次為模型1～20，模型1～4、5～8、9～12、13～16、17～20的d分別為2000、4000、6000、8000、10000條/m2

圖6 測線方位與裂縫走向夾角不同時模型8的聲波位移波場快照(t=8 μs)(a)0°； (b)60°； (c)120°； (d)180°； (e)240°； (f)300°

圖7 模型1～4的各向異性特征曲線(a)VP隨測線方位角變化； (b)α隨測線方位角變化

圖8 模型5～8的各向異性特征曲線(a)VP隨測線方位角變化； (b)α隨測線方位角變化

圖9 模型9～12的各向異性特征曲線(a)VP隨測線方位角變化； (b)α隨測線方位角變化

圖10 模型13～16的各向異性特征曲線(a)VP隨測線方位角變化； (b)α隨測線方位角變化

圖11 模型17～20的各向異性特征曲線(a)VP隨測線方位角變化； (b)α隨測線方位角變化

以模型8為例。VP在測線方位角為0°或180°時較大，在測線方位角為60°和240°時出現(xiàn)極大值，在測線方位角為90°時并未達到最小值，而在測線方位角為105°和285°時值達到最小值(圖8a)。這是因為裂縫分布是隨機離散的，除了在測線方位角為0°或180°時，在其他測線方位角聲波繞射通過裂縫到達接收端的路程也可能較短，因此VP較高。VP變化曲線呈多峰值分布，且隨著裂縫密度d增大，多峰值分布現(xiàn)象越明顯。因此，VP最大值或最小值對應(yīng)的測線方位難以準(zhǔn)確地指示裂縫走向。

α在測線方位角為0°或180°時達到最小值，在測線方位角為60°和240°出現(xiàn)極小值，在測線方位角為120°和300°時達到最大值(圖8b)。因此，α最小值對應(yīng)的測線方位可準(zhǔn)確地指示裂縫走向，α最大值對應(yīng)的測線方位判定裂縫走向存在一定偏差。隨著d增加，α基本呈雙峰值分布，均在測線方位角為180°附近出現(xiàn)波谷。這是因為盡管在某些測線方位角聲波繞射通過裂縫到達接收端的路程較短，但裂縫嚴(yán)重阻擋了聲波能量傳播，即使VP較大，α也較大。本文對隨機離散裂縫模型縱波方位各向異性特征的數(shù)值模擬結(jié)果表明，VP不一定在測線方位角為0°或180°時出現(xiàn)最大值，也不一定在90°或270°時出現(xiàn)最小值。因此，只有α最小值對應(yīng)的測線方位才可較準(zhǔn)確地指示裂縫走向。本文進一步深化了楊勤勇等[8]有關(guān)裂縫介質(zhì)的縱波方位各向異性特征的認(rèn)識。

為了分析各向異性的差異，參考Thomsen[20]分析縱波各向異性程度的方法，定義VP相對變化量ΔVP和α相對變化量Δα表征各向異性程度

(3)

(4)

式中下標(biāo)max、min分別代表對應(yīng)變量的極大值、極小值。由式(3)、式(4)得到模型8的ΔVP和Δα分別為0.20和0.81，即Δα?ΔVP，因此衰減各向異性特征較速度各向異性特征更顯著。

圖12、圖13分別為ΔVP、Δα隨d的變化。由圖可見, ΔVP隨d的變化規(guī)律不明顯(圖12)， Δα隨d增加而減小(圖13)。造成上述現(xiàn)象的原因為：VP還受裂縫空間分布位置的影響，因此ΔVP與d不存在明顯的相關(guān)關(guān)系；當(dāng)孔隙度φ一定時，隨著d增加，由裂縫條數(shù)較少的大裂縫轉(zhuǎn)變?yōu)榱芽p條數(shù)眾多的細小裂縫，空間分布逐漸均勻，各向異性特征減小，因此Δα隨d增加而減小。楊勤勇[21]選用4種厚度的有機玻璃片制作不同裂縫密度的均勻裂縫巖心模型進行物理模型實驗，結(jié)果表明，隨著裂縫密度減小，測線方位角為0°和90°時的反射時差逐漸增大，即縱波方位各向異性隨測線方位角增大逐漸增大。本文對隨機離散裂縫模型縱波方位各向異性特征的數(shù)值模擬結(jié)果不支持文獻[21]的認(rèn)識。

圖12 ΔVP隨d的變化

3.2 裂縫流體對方位各向異性特征的影響

圖14為裂縫模型，圖15、圖16分別為VP、α隨測線方位角的變化曲線。由模擬結(jié)果可見： 當(dāng)測線方位角為0°或180°時，裂縫及充填流體對聲波首φ=4%，θ=30°，d=4800條/m2，r均值為20。裂縫中流體為氣/水兩相，裂縫內(nèi)孔隙流體的含水飽和度分別為15%、30%、45%、60%、75%、90%、95%和100%波影響很小，因此不同飽和度間的VP(圖15)和α(圖16)的變化差異很?。划?dāng)測線方位角較大時，裂縫嚴(yán)重阻擋聲波首波傳播，因此不同飽和度間的VP(圖15)和α(圖16)的變化差異明顯。

圖13 Δα隨d的變化R為相關(guān)系數(shù)

圖14 裂縫模型

圖15 VP隨測線方位角的變化曲線

圖16 α隨測線方位角的變化曲線

圖17、圖18分別為測線方位角為90°或270°時VP、α隨含水飽和度的變化曲線。由圖可見：①當(dāng)含水飽和度較小時，隨著含水飽和度增加，VP先減小并出現(xiàn)最低值(圖17)，α逐漸增大并出現(xiàn)最大值(圖18)，這與氣/水兩相孔隙流體的微觀分布模式有關(guān)[22]；隨著含水飽和度進一步增加，VP逐漸增大(圖17)，α逐漸減小(圖18)。②當(dāng)含水飽和度高于90%時，VP迅速增大(圖17)，α迅速減小(圖18)。這是因為隨著含水飽和度增加，聲波透射能量增強，聲波透射通過裂縫的時間越短，透射波對首波初至的影響也越明顯，導(dǎo)致VP增加、α減小。

圖19、圖20分別為ΔVP、Δα隨含水飽和度的變化曲線。由圖可見，隨著含水飽和度增加，ΔVP、Δα均先增大后減小。由于當(dāng)測線方位角為0°或180°時VP、α的變化很小，因此ΔVP和Δα主要由測線方位角較大時VP和α的變化決定。由圖15可見，含水飽和度分別為15%、60%、100%的VP最大值幾乎相同，不同含水飽和度的VP最小值先減小后增加，因此ΔVP先增大后減小。同理，由圖16可見，含水飽和度分別為15%、60%、100%的α最小值幾乎相同，不同含水飽和度的α最大值先增加后減小，因此Δα先增大后減小。

圖17 測線方位角為90°或270°時VP隨含水飽和度的變化曲線

圖18 測線方位角為90°或270°時α隨含水飽和度的變化曲線

圖19 ΔVP隨含水飽和度的變化曲線

圖20 Δα隨含水飽和度的變化曲線

丁拼搏等[23]在不同飽和流體條件下利用超聲波透射法測試含裂縫分布的巖樣，結(jié)果表明，巖樣在飽氣條件下的縱波各向異性系數(shù)遠高于巖樣飽水。受基質(zhì)滲透性和裂縫分布狀態(tài)的影響，實驗中往往難以控制裂縫流體的部分飽和分布狀態(tài)，本文的數(shù)值模擬結(jié)果支持丁拼搏等的實驗結(jié)果，還討論了氣/水兩相部分飽和狀態(tài)下隨機離散裂縫模型的縱波各向異性特征，并進一步總結(jié)了縱波速度和衰減系數(shù)相對變化量隨含水飽和度的變化規(guī)律。

4 結(jié)論

本文對隨機離散裂縫模型的聲波波場的數(shù)值模擬結(jié)果表明：裂縫參數(shù)變化對聲波衰減系數(shù)的影響遠遠大于對聲波速度的影響，采用聲波衰減系數(shù)最小值對應(yīng)的測線方位可較準(zhǔn)確地判定裂縫走向；隨著裂縫密度增加，衰減系數(shù)相對變化量減小；隨著含水飽和度增加，速度和衰減系數(shù)的相對變化量均先增大后減小。上述認(rèn)識為利用聲學(xué)參數(shù)檢測裂縫提供了理論基礎(chǔ)。