張 冰，郭智奇，徐 聰，劉 財，劉喜武，劉宇巍

1.吉林大學地球探測科學與技術(shù)學院，長春 130026 2.中國電力工程顧問集團東北電力設(shè)計院有限公司，長春 130021 3.頁巖油氣富集機理與有效開發(fā)國家重點實驗室，北京 100083 4.中國石化頁巖油氣勘探開發(fā)重點實驗室，北京 100083 5.中國石化石油勘探開發(fā)研究院，北京 100083

0 引言

頁巖儲層中的非常規(guī)資源已經(jīng)成為油氣勘探的研究熱點。頁巖結(jié)構(gòu)復(fù)雜，內(nèi)部物質(zhì)的不均勻分布和定向排列使其通常表現(xiàn)出各向異性；可以通過研究頁巖儲層內(nèi)部結(jié)構(gòu)屬性和各向異性的聯(lián)系，為識別和勘探富含油氣區(qū)域提供依據(jù)。

巖石物理模型是建立頁巖微觀結(jié)構(gòu)與各向異性聯(lián)系的常用工具[1-2]。但頁巖各向異性受多種因素影響，包括黏土礦物的定向排列、非球形孔隙、微裂隙和裂縫屬性[3]。黏土礦物定向排列是由頁巖的壓實作用引起的，不同程度的壓實會引起黏土礦物不同的排列狀態(tài)，進而影響各向異性。Hornby等[4]和Sayers[5]使用方位分布函數(shù)(ODF)來模擬黏土礦物的方位分布狀態(tài)，然后使用各向異性微分等效介質(zhì)模型(DEM)來估計頁巖各向異性。Sayers[6]和Bandyopadhyay[3]指出非球形孔隙是頁巖各向異性的另一個來源。Kwon等[7]發(fā)現(xiàn)頁巖中裂隙的形狀、大小和定向排列影響各向異性。Chapman模型[8]使用了多尺度孔隙系統(tǒng)，包括球形孔隙、裂隙和裂縫，并假設(shè)各向異性只來源于水平排列的裂縫。Jiang等[9]建立了一個雙孔隙系統(tǒng)模型，即結(jié)合自相容等效介質(zhì)理論(SCM)和Chapman模型來克服Chapman模型中不包含非球形孔隙的缺點，并將該模型應(yīng)用于Haynesville頁巖。

對于裂縫型頁巖儲層，各向異性由裂縫屬性主導(dǎo)。本文在巖石物理建模過程中參照Jiang等[9]的工作，使用雙孔隙系統(tǒng)模型以準確描述頁巖的孔隙系統(tǒng)，并針對實際研究區(qū)域?qū)υ撃Ｐ瓦M行校正。由于模型的復(fù)雜性，反演需要確定非球形孔隙度、孔隙縱橫比和裂縫密度，因此反演過程使用模擬退火優(yōu)化粒子群算法(SA-PSO)[10]。作為一種智能尋優(yōu)算法，它能夠快速準確地尋找到各參數(shù)最優(yōu)解，適用于復(fù)雜的非線性反演問題。該算法被應(yīng)用到位于中國西南部的龍馬溪組頁巖，估計了該儲層物性參數(shù)，如非球形孔隙度、孔隙縱橫比、裂縫密度和各向異性參數(shù)。

1 巖石物理模型

橫向各向同性(TI)頁巖的各向異性屬性可以由Thomsen[11]各向異性參數(shù)表示：

(1)

式中：cij是剛度張量分量，由巖石內(nèi)部屬性決定；vP、vS和ρ分別是縱波速度、橫波速度和密度；ε、γ、δ是各向異性參數(shù)。巖石物理模型能夠描述巖石內(nèi)部屬性并得到等效彈性張量，然后使用實測縱、橫波速度對模擬結(jié)果進行校準。

包絡(luò)體理論是頁巖建模的一種常用方法。該理論把巖石中的組分當作不同形狀的包含物，形狀由縱橫比決定。SCM[12]是一種各向同性包絡(luò)體理論，能夠結(jié)合多種礦物和孔隙類型構(gòu)建等效介質(zhì)模型，等效體積模量Ksc和等效剪切模量μsc由下式給出：

(2)

式中：N是組分數(shù)目；fi是每種組分的體積分數(shù)；Ki和μi是每種組分的體積模量和剪切模量；Ui和Vi是幾何因子，由縱橫比和模量決定，具體公式由Berryman[13]給出。圖1顯示了SCM模型中不同孔隙縱橫比(α)和孔隙度(φ)下vP和vS的變化情況。所用基質(zhì)是黏土，顆?？v橫比為0.1，孔隙為水飽和。由圖1可見：隨著α的增加，vP和vS逐漸增大；而隨著φ的增加，vP和vS都變小。

Chapman多尺度孔隙模型[8]可用來描述流體噴射流動現(xiàn)象對巖石等效模量的影響。波長小于孔隙尺寸的波在巖石中傳播時會產(chǎn)生壓力梯度，噴射流動是指由此產(chǎn)生的流體相互作用[14]。該模型考慮了頻率和孔隙壓力的影響，適用于各向同性和各向異性介質(zhì)[9]。Chapman模型認為巖石包含球形孔隙、裂隙和裂縫。圖2是該模型的示意圖[15]，其中球形孔隙和裂隙隨機分布，裂縫水平排列，各向異性由裂縫引起。除了裂縫與裂縫之間，球形孔隙、裂隙和裂縫相互連通來保證流體的流動性[8]。Chapman模型中的孔隙度(φchap)可表示為球形孔隙度(φrp)、裂隙孔隙度(φc)和裂縫孔隙度(φf)三者之和：

φchap=φrp+φc+φf。

(3)

裂隙和裂縫被認為是理想橢球體，具有相同的縱橫比，因此裂隙密度(εc)和裂縫密度(εf)由下式得出：

(4)

Chapman模型的具體公式由Jiang等[9]總結(jié)給出。模型中參數(shù)都有相應(yīng)的物理意義，并且大部分可以直接測定或通過測定值估計得到。圖3a是孔隙縱橫比對各向異性參數(shù)的影響。由圖3a可見：vS、ε和δ隨著α的增大而增大，其中δ的變化最大；而vP和γ卻隨著α的增大而減小。圖3b為各向異性參數(shù)隨裂縫密度的變化。由圖3b可見：vP、vS和δ隨著εf的增大而減?。欢藕挺脜s隨著εf的增大而增大，其中γ變化最大。模型所用基質(zhì)是黏土，孔隙為水飽和。

Chapman模型中孔隙被認為是球形，但大部分頁巖包含非球形孔隙。Jiang等[9]將SCM和Chapman模型相結(jié)合構(gòu)建了雙孔隙系統(tǒng)模型來估計Haynesville頁巖的儲層屬性。在建模過程中，SCM模型用來描述礦物和非球形孔隙，非球形孔隙間相互不連通。Chapman模型包含球形孔隙、裂隙和裂縫，其中裂縫間相互不連通，裂隙和球形孔隙至少與一個裂縫相連通，用來描述流體噴射流動現(xiàn)象對巖石等效模量的影響。該雙孔隙系統(tǒng)模型將總孔隙度(φ)劃分為非球形孔隙(φnp)和Chapman模型孔隙系統(tǒng)，這兩個孔隙系統(tǒng)相互不連通，可表示為

a.φ=0.1；b. α=0.05。 圖1 SCM模型中縱、橫波速度隨孔隙縱橫比和孔隙度的變化情況Fig.1 Modeled P- and S- velocities as a function of pore aspect ratio and porosity by using SCM

圖中包括隨機分布的球形孔隙、隨機分布并且隨機方位的裂隙和水平定向排列的裂縫。裂隙和裂縫的縱橫比相同，孔隙、裂隙和礦物顆粒的尺寸相同。據(jù)文獻[15]。圖2 Chapman模型示意圖Fig.2 Schematic of Chapman’s model

φ=φnp+φchap。

(5)

使用SCM模型時需要考慮各包絡(luò)體的縱橫比，本文中黏土礦物縱橫比是0.10，脆性礦物，例如石英、方解石、黃鐵礦等是1.00，干酪根是0.01。SCM模型中影響速度的主要因素是孔隙形狀，在建模流程中，SCM模型中的非球形孔隙與Chapman模型孔隙系統(tǒng)中的裂隙和裂縫有相同的縱橫比，且非球形孔隙尺寸與礦物顆粒以及球形孔隙、裂隙尺寸一致[9]。裂縫尺寸對噴射流動現(xiàn)象影響很大，通?？赏ㄟ^巖心數(shù)據(jù)確定。

雙孔隙系統(tǒng)頁巖巖石物理建模流程由圖4給出。建模過程中，首先使用SCM模型構(gòu)建包含礦物和非球形孔隙的各向同性基質(zhì)，礦物和孔隙在該基質(zhì)中均勻并且隨機方位分布且孔隙互不連通。然后使用Chapman模型向該基質(zhì)中加入球形孔隙、裂隙和裂縫的影響，其中裂隙的方位隨機，裂縫呈水平排列，最終得到各向異性頁巖模型。在該巖石物理模型中：各向異性來源于水平裂縫，主要受孔隙縱橫比和裂縫密度影響；而各向同性基質(zhì)的縱、橫波速度受非球形孔隙度和孔隙縱橫比影響；Chapman模型中其他參數(shù)可根據(jù)實際數(shù)據(jù)給定。圖5是巖石物理模型中非球形孔隙度、縱橫比、裂縫密度這3個重要儲層參數(shù)對各向異性參數(shù)的影響，所用基質(zhì)是黏土，孔隙是水飽和。其中：φnp對vS和γ的影響較大；α較小時對各向異性參數(shù)的影響大，但隨著α增加，影響逐漸減?。沪舊的取值對各向異性參數(shù)的影響始終較大。

a. εf=0.1，φ=0.1，εc=0.025；b. α=0.05，φ=0.1，εc =0.025。圖3 Chapman模型得到的各向異性參數(shù)隨孔隙縱橫比和裂縫密度的變化情況Fig.3 Modeled anisotropy parameters as a function of pore aspect ratio and fracture density by using Chapman’s model

圖4 頁巖巖石物理模型Fig.4 Rock physics model for shale

a. α =0.05，φ=0.1，εf =0.1，εc=0.025；b. φnp=0.02，φ=0.1，εf=0.1，εc=0.025；c. φnp=0.02，α =0.05，φ=0.1，εc=0.025。圖5 巖石物理模型中各向異性參數(shù)隨非球形孔隙度、孔隙縱橫比和裂縫密度的變化情況Fig.5 Modeled anisotropy parameters as a function of non-spherical porosity, pore aspect ratio and fracture density by using the rock physics model

2 實際數(shù)據(jù)分析

研究對象是位于中國西南四川盆地的龍馬溪組頁巖。測井數(shù)據(jù)在圖6中顯示，包括伽馬射線、縱橫波速度、密度、孔隙度和組分信息。目標層段頁巖致密，孔隙度較低，干酪根和黏土含量高。巖心資料顯示該層段水平裂縫發(fā)育，裂縫尺寸平均值是100 μm。

目標層段是包含水平裂縫的橫向各向同性(TI)頁巖。鄧繼新等[16]對龍馬溪組頁巖微觀結(jié)構(gòu)的研究發(fā)現(xiàn)：該頁巖的孔隙類型包含原生孔隙、次生孔隙和裂隙，各種孔隙尺度不同，但大部分小于2 μm；孔隙大部分是非球形孔隙，且呈孤立分布，連通性差。

巖石物理模型中參數(shù)眾多，部分參數(shù)可由前人研究給定。Guo等[17]直接使用了Chapman[8]給出的頻率和弛豫時間參數(shù)值。Jiang等[9]使用雙孔隙系統(tǒng)模型對Haynesville頁巖儲層參數(shù)進行了估計。龍馬溪組頁巖類比于Haynesville頁巖，兩者主要礦物組分都是黏土、石英和方解石，孔隙度和滲透性都很低，并且都具有橫向各向同性，因此認為兩者具有相似性。據(jù)此，我們?yōu)閹r石物理建模過程中部分參數(shù)賦予了相同值。其中：Chapman模型中頻率(f)取10 kHz，接近測井頻率；非球形孔隙、球形孔隙、裂隙的尺寸(ξ)相同，取1 μm；裂隙密度由裂縫密度得到，

(6)

水飽和裂隙的弛豫時間(τm)由Chapman[8]給出；而裂縫的弛豫時間(τf)由它們的尺寸(af)決定，

(7)

礦物和流體的模量和密度由表1給出。最終，巖石物理模型中的未知參數(shù)可確定為非球形孔隙度、孔隙縱橫比和裂縫密度。通過使用測量得到的縱、橫波速度可反演得到這3個儲層參數(shù)，進而估計頁巖各向異性。

GR. 自然伽馬。圖6 頁巖測井數(shù)據(jù)Fig.6 Shale information from the well data

體積模量/GPa剪切模量/GPa密度/(g/cm3)黏土干酪根石英方解石黃鐵礦水26.02.938.077.0147.02.816.02.744.032.0132.002.551.302.652.714.931.09

3 反演方法和結(jié)果

反演目標函數(shù)通?？蓪懽?/p>

F(φnp,α,εf)=min(‖vP-real-vP-model‖2+

‖vS-real-vS-model‖2)。

(8)

式中：vP-real和vS-real是縱、橫波速度的測量值；vP-model和vS-model是由巖石物理模型計算的縱、橫波速度。

傳統(tǒng)反演方法基于遍歷搜索，即將各參數(shù)取值空間劃分成n個區(qū)間，遍歷各參數(shù)的各取值區(qū)間來尋找滿足目標函數(shù)的值[18]。對于本文中的三參數(shù)反演問題，傳統(tǒng)方法需要迭代n3次。巖石物理模型復(fù)雜，若區(qū)間數(shù)過多則計算量大，過少則無法找到最優(yōu)解。而且，多參數(shù)反演存在多解性，需要加入更多反演約束條件來減少多解性。地震反演方法通常在目標函數(shù)中加入約束項來減少多解性，李坤等[19]使用平滑模型作為約束項來減少地震反演中低頻趨勢的多解性。本文考慮到實際儲層屬性存在垂向連續(xù)性、各深度點儲層參數(shù)并不完全獨立，在目標函數(shù)中加入平滑約束項后為

Fs(φnp,α,εf)=min(‖vP-real-vP-model‖2+

(9)

基于巖石物理模型和模擬退火優(yōu)化粒子群算法的頁巖儲層參數(shù)反演流程(圖7)如下。

1)初始化參數(shù)。從參數(shù)取值區(qū)間抽取M個隨機點，被稱作粒子：

s=[0, 0.5(Su-Sd)]。

(10)

圖7 基于巖石物理模型和模擬退火優(yōu)化粒子群算法的頁巖儲層參數(shù)反演流程Fig.7 Workflow for rock properties inversion by using shale rock physics model and SA-PSO

本問題中3個儲層參數(shù)的取值范圍分別是：

Sφnp=[10-4,0.8φ]，Sα=[10-6,0.2]，

Sεf=[10-4,0.2]。

(11)

2)計算粒子適應(yīng)度。

fit(xi)=-Fs(xi)。

(12)

式中：fit是粒子適應(yīng)度函數(shù)，由目標函數(shù)值決定；xi為粒子位置。

3)移動粒子。

(13)

(14)

式中：Tk是當前溫度；r3是屬于[0,1]的隨機數(shù)。

4)降溫操作。

Tk+1=cTk。

(15)

式中，c是降溫系數(shù)，取0.93。

5)重復(fù)步驟2)至4)直到滿足收斂條件。

對目標層段進行反演，使用100個粒子，迭代60次。圖8展示了一個深度點處的反演迭代收斂過程。由圖8可以看出，粒子的適應(yīng)度由低變高，并且逐漸聚集到最優(yōu)解處。最優(yōu)解滿足目標函數(shù)公式(9)，此時粒子在三維參數(shù)空間的位置也就是各參數(shù)的取值。如果傳統(tǒng)方法中每個參數(shù)區(qū)間劃分50份，則需要迭代125 000次，而模擬退火優(yōu)化粒子群算法共需迭代6 000次，計算量得以大大降低。

圖9是目標層段使用模型計算出的縱、橫波速度值和實測值，兩者擬合很好，表明了本文方法的準確性。圖10是反演得到的儲層參數(shù)曲線，可清楚地看出儲層參數(shù)隨深度的變化趨勢，并首次提供了該研究區(qū)域的非球形孔隙度。Qian等[20]對龍馬溪組頁巖的研究表明，該層段孔隙縱橫比較小，平均低于0.05，與本文反演值(圖10)一致。Jiang等[9]對Haynesville頁巖的研究給出它的裂縫密度是0.04，由于本文目標層裂縫發(fā)育，反演值較高(圖10)。圖11是估計的各向異性參數(shù)，其中ε和γ具有強相關(guān)性，而γ和δ在部分區(qū)域出現(xiàn)負相關(guān)；該現(xiàn)象與Bandyopadhyay[3]及Sayers[6]的研究一致。目標層段各向異性參數(shù)γ總體大于ε，并且δ為負值；Qian等[20]的研究也給出了相似的結(jié)果，但目標層各向異性強度由于裂縫發(fā)育而更大。

a. 初始化后粒子狀態(tài)；b. 迭代20次后；c. 迭代40次后；d. 迭代60次后。圖8 模擬退火優(yōu)化粒子群算法收斂過程Fig.8 Optimization procedure of SA-PSO

圖9 縱、橫波速度計算值(實線)和實測值(虛線)Fig.9 Modeled (solid) and measured (dashed) P- and S-wave velocities

圖10 目標層段儲層參數(shù)反演結(jié)果Fig.10 Inversion result of the formation properties

圖11 目標層各向異性參數(shù)反演結(jié)果Fig.11 Inversion result of the anisotropy parameters

4 結(jié)論與展望

本文使用改進的巖石物理模型和模擬退火優(yōu)化粒子群算法對頁巖裂縫屬性和各向異性進行了反演。巖石物理建模為使用等效介質(zhì)模型和Chapman模型相結(jié)合建立的雙孔隙系統(tǒng)模型，以準確地描述頁巖微觀孔隙結(jié)構(gòu)，反演過程使用智能算法準確求解復(fù)雜巖石物理模型下三參數(shù)反演問題。將方法應(yīng)用到龍馬溪組頁巖氣儲層，反演結(jié)果與實測數(shù)據(jù)和已有研究結(jié)果相匹配，說明了本文方法的有效性。

本文的研究為頁巖儲層的描述提供了更多的有效信息。頁巖結(jié)構(gòu)復(fù)雜，影響各向異性的因素很多，進一步的研究需要使用更加準確復(fù)雜的巖石物理模型進行模擬。而求解復(fù)雜巖石物理模型，反演方法同樣需要進一步提高。