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

段 茜 劉向君 梁利喜 熊 健

(①油氣藏地質(zhì)及開(kāi)發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室(西南石油大學(xué))，四川成都 610500； ②西南石油大學(xué)理學(xué)院，四川成都 610500)

0 引言

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

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

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

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

1.1 計(jì)算聲波速度

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

(1)

即可計(jì)算聲波速度VP。式中：L為巖樣長(zhǎng)度；TP為縱波傳播時(shí)間。

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

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

(2)

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

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

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

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

2.2 隨機(jī)離散裂縫模型的建立

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

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

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

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

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

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

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

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

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

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

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

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

2.3 實(shí)驗(yàn)觀測(cè)系統(tǒng)

圖4為縱波方位各向異性特征的實(shí)驗(yàn)觀測(cè)系統(tǒng)。

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

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

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

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

圖5 隨機(jī)離散裂縫模型巖心尺寸為50mm(長(zhǎng)度)×50mm(直徑)，時(shí)間采樣步長(zhǎng)為10ns，空間采樣步長(zhǎng)為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 測(cè)線方位與裂縫走向夾角不同時(shí)模型8的聲波位移波場(chǎng)快照(t=8 μs)(a)0°； (b)60°； (c)120°； (d)180°； (e)240°； (f)300°

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

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

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

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

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

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

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

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

(3)

(4)

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

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

圖12 ΔVP隨d的變化

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

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

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

圖14 裂縫模型

圖15 VP隨測(cè)線方位角的變化曲線

圖16 α隨測(cè)線方位角的變化曲線

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

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

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

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

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

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

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

4 結(jié)論

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