基于空間響應(yīng)分布函數(shù)的水平井補(bǔ)償密度測井快速正演模擬

袁　超, 李潮流, 周燦燦, 葛云龍

(1.中國石油勘探開發(fā)研究院,北京 100083; 2.中國石油大學(xué)(北京)地球物理與信息工程學(xué)院,北京 102249)

隨著鉆井技術(shù)的發(fā)展,為了提高石油天然氣采收率,水平井的數(shù)量逐年增加[1-2]。井眼軌跡與地層關(guān)系是水平井測井解釋的關(guān)鍵和基礎(chǔ)[3],依據(jù)正演模擬曲線與實測曲線對比確定兩者關(guān)系時需要根據(jù)初始地層模型快速正演測井曲線。如果正演曲線與實測曲線吻合,則說明設(shè)定地層模型符合真實地層情況,否則需要重新調(diào)整地層模型[4]。目前,主要利用電阻率、自然伽馬和密度測井曲線確定井眼軌跡與地層關(guān)系,國內(nèi)外大量學(xué)者開展電阻率和自然伽馬測井快速正演模擬的研究[5-9]。國外學(xué)者對補(bǔ)償密度測井快速正演做了相關(guān)研究,Booth等[10]研究了正演模擬中的重要性估算;Watson[11]利用蒙特卡羅模擬探測器微分靈敏度函數(shù);Case等[12]研究了放射性測井建模中的靈敏度函數(shù);Liu等[13]基于精細(xì)網(wǎng)格研究了探測器響應(yīng)的重要性;Mendoza等[14]利用通量靈敏度函數(shù)和格林函數(shù)快速模擬核測井響應(yīng)。雖然國外學(xué)者對補(bǔ)償密度測井快速正演模擬做了一系列的研究,但是國內(nèi)學(xué)者主要采用蒙特卡羅研究測井響應(yīng)特性[15-16],未見關(guān)于補(bǔ)償密度測井快速正演模擬的公開發(fā)表文獻(xiàn)。蒙特卡羅模擬計算耗時長,不滿足井眼軌跡與地層關(guān)系確定中快速正演模擬的需求[17]。筆者提出一種補(bǔ)償密度測井快速正演模擬方法,利用蒙特卡羅程序建立快速正演模擬所需數(shù)據(jù)庫,計算不同地層區(qū)域?qū)﹂L、短源距探測器計數(shù)的貢獻(xiàn)近似值,進(jìn)而得到長、短源距探測器計數(shù),并將探測器計數(shù)轉(zhuǎn)換為地層密度,以快速獲取高精度的補(bǔ)償密度測井正演模擬結(jié)果。

1　補(bǔ)償密度測井快速正演模擬原理

伽馬源發(fā)射出的伽馬射線進(jìn)入地層后,在地層中的輸送過程可用玻爾茲曼方程來表示,伽馬射線在地層中輸送過程示意圖如圖1所示。利用玻爾茲曼方程可以得出探測器位置處伽馬射線計數(shù)為

(1)

式中,r為伽馬射線位置;E為伽馬射線能量;Ω為入射伽馬射線方向;ψ為伽馬射線角通量;rS為伽馬放射源位置;rD為探測器位置;R為探測器響應(yīng)函數(shù)。

從公式(1)可以看出,補(bǔ)償密度測井儀器在探測器位置處的伽馬射線計數(shù)與入射伽馬參數(shù)(包括入射伽馬能量和方向)、伽馬射線角通量、發(fā)生康普頓散射反應(yīng)的伽馬射線位置以及探測器響應(yīng)函數(shù)(與探測器位置、發(fā)生康普頓散射反應(yīng)的伽馬射線位置以及入射伽馬射線參數(shù)有關(guān))。在補(bǔ)償密度測井中,入射伽馬能量和方向等參數(shù)已知,因此在探測器位置處伽馬射線計數(shù)與伽馬射線角通量、發(fā)生康普頓散射的伽馬射線位置以及探測器響應(yīng)函數(shù)有關(guān)。因此探測器位置處伽馬射線計數(shù)受3個因素控制,即伽馬射線在地層中的空間通量分布、地層中不同位置處發(fā)生康普頓散射反應(yīng)的伽馬射線對探測器計數(shù)的貢獻(xiàn)以及探測器響應(yīng)函數(shù)有關(guān)[18]。

圖1　補(bǔ)償密度測井中伽馬射線輸送過程示意圖Fig.1　Schematic diagram of gamma ray transportation in compensated density logging

定義伽馬放射源發(fā)射出的伽馬射線對探測器計數(shù)貢獻(xiàn)為空間響應(yīng)分布函數(shù)FS:

FS=Φη.

(2)

式中,FS為空間響應(yīng)分布函數(shù),無量綱;Φ為伽馬射線的空間通量分布;η為伽馬射線與空間不同位置體積元發(fā)生康普頓散射后對探測器計數(shù)的貢獻(xiàn)。

空間響應(yīng)分布函數(shù)FS的物理意義為地層中不同位置體積元對補(bǔ)償密度測井伽馬探測器計數(shù)的貢獻(xiàn)。利用蒙特卡羅方法可以模擬單一密度地層中的長、短源距探測器處的空間響應(yīng)分布函數(shù)FS以及長、短源距探測器的伽馬計數(shù)。

對于如圖2所示的兩種密度組合地層模型中,可認(rèn)為探測器計數(shù)是區(qū)域Z1和Z2對探測器計數(shù)貢獻(xiàn)的線性疊加?？焖僬菽M中,選取密度接近組合地層密度平均值的地層作為參考地層,在不同密度地層區(qū)域Z1和Z2,對參考地層的空間響應(yīng)分布函數(shù)積分得到每個區(qū)域?qū)﹂L、短源距探測器計數(shù)的貢獻(xiàn)近似值,將不同地層密度區(qū)域的貢獻(xiàn)近似值歸一化后與相應(yīng)地層單一密度條件下記錄的長、短源距探測器計數(shù)相乘并累加,計算得到長、短源距探測器計數(shù),即

(3)

(4)

式中,NS和NL分別為快速正演模擬的短源距和長源距探測器計數(shù);FSRS和FSRL分別為短源距和長源距探測器位置處利用蒙特卡羅方法模擬計算的參考地層的空間響應(yīng)分布函數(shù);NρS1和NρS2分別為區(qū)域Z1和Z2的單一密度地層中短源距探測器計數(shù);NρL1和NρL2分別為區(qū)域Z1和Z2的單一密度地層中長源距探測器計數(shù)。

圖2　兩種密度組合地層模型橫截面示意圖Fig.2　Cross section of formation model with two density zones

(5)

(6)

則公式(3)和(4)變?yōu)?/p>

NS?GrSNρS1+(1-GrS)NρS2,

(7)

NL?GrLNρL1+(1-GrL)NρL2.

(8)

模擬計算得到探測器計數(shù)后,通過短源距和長源距探測器計數(shù)與地層密度轉(zhuǎn)換關(guān)系將探測器計數(shù)轉(zhuǎn)換為地層密度。

2　快速正演模擬所需數(shù)據(jù)庫建立

利用蒙特卡羅方法在單一密度地層中建立補(bǔ)償密度測井快速正演模擬所需要的數(shù)據(jù)庫,包括單一密度地層與伽馬射線發(fā)生康普頓散射反應(yīng)對探測器計數(shù)貢獻(xiàn)的空間響應(yīng)分布函數(shù)、探測器計數(shù)以及探測器計數(shù)與地層密度的轉(zhuǎn)換關(guān)系。

2.1　蒙特卡羅計算模型

利用蒙特卡羅方法[19-20]建立儀器-井眼-地層的三維柱狀模型,計算模型設(shè)定參數(shù)如下:

(1)地層模型高度和半徑分別為60 和35 cm,對模擬地層進(jìn)行網(wǎng)格化以提高計算精度,每個網(wǎng)格尺寸為0.5 cm×0.5 cm。

(2)井眼半徑為10 cm,井眼內(nèi)充滿淡水。

(3)放射源采用Cs137,發(fā)射伽馬射線能量為662 keV。

(4)短源距探測器和長源距探測器的長度分別為5 和10 cm,距離放射源的距離分別為18 和42 cm。

(5)放射源和短源距探測器之間放置鎢屏蔽體,屏蔽體長度為3 cm。

2.2　空間響應(yīng)分布函數(shù)

圖3　蒙特卡羅計算模型橫截面示意圖Fig.3　Cross section of Monte Carlo simulation model

利用如圖3所示的蒙特卡羅計算模型,改變地層孔隙度為0、5%、10%、15%、20%、25%、30%、35%和40%,設(shè)定不同的地層密度為2.65、2.567 5、2.485、2.402 5、2.32、2.237 5、2.155、2.072 5和1.99 g/cm3,模擬伽馬射線的空間分布以及伽馬射線與空間不同位置體積元發(fā)生康普頓散射反應(yīng)對伽馬探測器計數(shù)的貢獻(xiàn),兩者乘積為空間響應(yīng)分布函數(shù)。對于未模擬地層密度的空間響應(yīng)分布函數(shù)需要利用插值的方法獲取。地層密度設(shè)定為2.485 g/cm3時,短源距和長源距探測器的空間響應(yīng)分布函數(shù)如圖4所示。

圖4　空間響應(yīng)分布函數(shù)Fig.4　Spatial response distribution function

從圖4可以看出,伽馬探測器計數(shù)的貢獻(xiàn)主要來自于放射源以及探測器附近地層,其他地層區(qū)域?qū)ゑR探測器計數(shù)影響較小。另外,當(dāng)計算模型發(fā)生變化時,空間響應(yīng)分布函數(shù)也會發(fā)生相應(yīng)變化,因此需要根據(jù)測井儀器類型以及不同井眼和地層條件,建立補(bǔ)償密度測井快速正演模擬所需要的空間響應(yīng)分布函數(shù)的數(shù)據(jù)庫。

2.3　地層密度刻度關(guān)系

在模擬空間響應(yīng)分布函數(shù)的同時記錄短源距和長源距探測器計數(shù),建立不同地層密度條件下短源距和長源距探測器計數(shù)的數(shù)據(jù)庫,并建立短源距和長源距探測器計數(shù)與地層密度的轉(zhuǎn)換關(guān)系數(shù)據(jù)庫。由于實際儀器測量結(jié)果為補(bǔ)償密度,尤其是當(dāng)測量結(jié)果受圍巖影響時,長、短源距密度與實際儀器響應(yīng)有較大誤差,通過脊肋圖的方法獲得補(bǔ)償密度,通過模擬地層模型與實際儀器響應(yīng)對比分析,獲得地層真實密度值。在設(shè)定的計算模型條件下,短源距和長源距探測器計數(shù)與地層密度的轉(zhuǎn)換關(guān)系如圖5所示。

如圖5所示的地層密度轉(zhuǎn)換關(guān)系的脊肋圖,在設(shè)定的計算模型條件下,脊線和肋線方程分別為

lnNL=3.261 5lnNS-2.453,

(9)

lnNL=0.996 2lnNS+b.

(10)

把長、短源距探測器計數(shù)NL和NS代入公式(9)肋線方程,得到截距b,與脊線方程(10)聯(lián)立,得到交點的長、短源距探測器計數(shù)NLC和NSC,密度校正量為

Δρ=k(lnNLC-lnNL).

式中,k為長源距密度斜率。則補(bǔ)償密度為

ρ=ρL+Δρ.

圖5　地層密度刻度關(guān)系Fig.5　Formation density calibration relationship

3　補(bǔ)償密度測井快速正演模擬及誤差分析

利用建立的數(shù)據(jù)庫,開展補(bǔ)償密度測井快速正演模擬,并對比快速正演模擬結(jié)果與蒙特卡羅模擬結(jié)果,分析快速正演模擬結(jié)果的誤差。

3.1　直井模型

在如圖2所示的地層模型中,區(qū)域Z1和Z2的地層密度分別為2.567和2.155 g/cm3,井軸與Z軸的夾角(相對地層傾角)為0°,測井儀器不斷上提模擬在直井中的測井過程。選取區(qū)域Z2為參考地層,根據(jù)在密度為2.155 g/cm3的單一密度地層中模擬的空間響應(yīng)分布函數(shù),在每個深度點近似計算區(qū)域Z1和Z2對探測器計數(shù)的貢獻(xiàn):

(11)

(12)

(13)

(14)

式中,CS1和CS2分別為地層區(qū)域Z1和Z2對短源距探測器計數(shù)貢獻(xiàn)近似值;CL1和CL2分別為地層區(qū)域Z1和Z2對長源距探測器計數(shù)貢獻(xiàn)近似值;FSRS和FSRL分別為數(shù)據(jù)庫中密度為2.155 g/cm3的參考地層與伽馬射線發(fā)生康普頓散射反應(yīng)對短源距和長源距探測器計數(shù)貢獻(xiàn)的空間響應(yīng)分布函數(shù)。

根據(jù)上述計算的地層區(qū)域Z1和Z2對短源距和長源距探測器計數(shù)貢獻(xiàn)的近似值,以及數(shù)據(jù)庫中在地層區(qū)域Z1和Z2的單一密度地層中記錄的探測器計數(shù),可得到在該地層模型中短源距和長源距探測器計數(shù):

(15)

(16)

式中,NS和NL分別為快速正演模擬計算的短源距和長源距探測器計數(shù);NS1和NS2分別為數(shù)據(jù)庫中區(qū)域Z1和Z2的單一密度地層中記錄的短源距探測器計數(shù);NL1和NL2分別為數(shù)據(jù)庫中在區(qū)域Z1和Z2的單一密度地層中記錄的長源距探測器計數(shù)。

在相同地層模型條件下,選取10個位置點,利用蒙特卡羅方法模擬短源距和長源距探測器計數(shù)。根據(jù)如圖5建立的短源距和長源距探測器計數(shù)與地層密度的轉(zhuǎn)換關(guān)系,將快速正演模擬和蒙特卡羅模擬的探測器計數(shù)轉(zhuǎn)換為地層密度,結(jié)果如圖6所示。

3.2　斜井模型

在如圖2所示的地層模型中,區(qū)域Z1和Z2的地層密度分別為2.567和2.155 g/cm3,井軸與Z軸的夾角(相對地層傾角)為70°,測井儀器不斷上提模擬在斜井中的測井過程。利用與在直井中快速正演模擬的方法,模擬計算每個深度點的短源距和長源距探測器計數(shù)。在相同地層模型條件下,選取10個位置點,利用蒙特卡羅方法模擬短源距和長源距探測器計數(shù)。根據(jù)如圖5建立的短源距和長源距探測器計數(shù)與地層密度的轉(zhuǎn)換關(guān)系,將快速正演模擬和蒙特卡羅模擬的探測器計數(shù)轉(zhuǎn)換為地層密度,結(jié)果如圖7所示。

圖6　直井中快速正演模擬與蒙特卡羅模擬結(jié)果對比Fig.6　Comparison between fast forward simulation and Monte Carlo simulation results in vertical wells

圖7　斜井中快速正演模擬與蒙特卡羅模擬結(jié)果對比Fig.7　Comparison between fast forward simulation and Monte Carlo simulation results in deviated wells

密度成像測井是在測井過程中記錄不同方位的地層密度值,實現(xiàn)井周地層密度的方位測量,并將地層密度值按方位更加直觀的顯示[21]。本文中還實現(xiàn)了密度成像測井的快速正演模擬,如圖8(a)所示,傾斜地層界面與井軸的交點為b,但在不同測量方位上地層界面與井軸的交點不同,利用柱坐標(biāo)下的地層界面與井筒方程可確定不同測量方位下地層界面與井軸的交點。按如圖8(b)所示的記錄扇區(qū)示意圖,在每個采樣點模擬計算16個扇區(qū)方位的地層密度測井值,得到的密度成像如圖9所示。

圖8　斜井地層界面與井筒關(guān)系及不同記錄扇區(qū)示意圖Fig.8　Relationship between formation interface and wellbore in deviated wells, and schematic diagram of different recording sectors

圖9為測井儀器從密度為2.567 g/cm3的地層上提到密度為2.155 g/cm3的地層中模擬測量的方位密度成像圖,模型的地層界面相對傾角α為70°。利用模擬的方位密度成像圖計算得出的地層界面相對傾角為68.5°,與地層真實相對傾角的相對誤差為2.1%。

圖9　密度成像測井模擬結(jié)果Fig.9　Results of density imaging logging simulation

3.3　模擬結(jié)果誤差分析

蒙特卡羅方法是放射性測井中常用的數(shù)值模擬方法,可以精確設(shè)定計算模型,得到高精度的數(shù)值模擬結(jié)果,所以將蒙特卡羅方法模擬結(jié)果作為理論精確值。在直井和斜井條件下的快速正演模擬和蒙特卡羅模擬結(jié)果的交會圖如圖10所示,圖中藍(lán)色線為對角線,計算結(jié)果越靠近藍(lán)色線表明精度越高。

圖10　快速正演模擬與蒙特卡羅模擬結(jié)果交會圖Fig.10　Cross plot of fast forward simulation and Monte Carlo simulation results

為驗證快速正演模擬結(jié)果的精確性,利用如下公式計算快速正演模擬與蒙特卡羅模擬結(jié)果的絕對誤差平均值EAA、相對誤差平均值EAR及皮爾遜相關(guān)系數(shù)Pr:

(17)

(18)

(19)

式中,n為采樣數(shù)據(jù)點的個數(shù);SFi為在第i采樣點的快速正演模擬結(jié)果;SMi為在第i采樣點的蒙特卡羅方法模擬的結(jié)果。

快速正演模擬和蒙特卡羅模擬結(jié)果的絕對誤差平均值、相對誤差平均值及相關(guān)系數(shù)計算值以及計算速度如表1所示。

表1　快速正演模擬和蒙特卡羅模擬結(jié)果誤差分析和計算速度Table 1　Error analysis and calculation speed of fast simulation and Monte Carlo simulation results

根據(jù)圖10所示的交會圖以及表1所示誤差分析和計算速度結(jié)果,可以看出快速正演模擬與蒙特卡羅模擬結(jié)果基本一致,快速正演模擬可提供與蒙特卡羅模擬精度一致的補(bǔ)償密度測井正演模擬結(jié)果。但是,本文快速正演模擬方法的計算速度非?？?利用快速正演模擬方法計算本文中直井或斜井實例的密度測井曲線可在2 s內(nèi)完成,而利用蒙特卡羅方法則需要大約60 h,計算速度提高了十萬倍以上。本文方法可在極短時間內(nèi)得到與蒙特卡羅方法精度基本一致的補(bǔ)償密度測井正演模擬結(jié)果,可以滿足水平井井眼軌跡與地層關(guān)系確定補(bǔ)償密度測井實時正演模擬計算的需求。

4　快速正演模擬應(yīng)用實例

圖11為在四川頁巖氣水平井X井中快速正演模擬實例,第1道為自然伽馬曲線,第2道為深度道,第3道為相位電阻率曲線,第4道為密度測井曲線及提取的地層模型,第5道為實際測量的隨鉆方位密度成像圖像,第6道為根據(jù)密度測井曲線提取的地層模型快速正演模擬計算的方位密度成像圖像。該模擬實例中30 m的井段正演模擬計算時間約需要20 s,快速正演模擬計算的結(jié)果與實際測井資料基本一致。

圖11　頁巖氣X井中快速正演模擬應(yīng)用Fig.11　Application of fast forward simulation in shale gas X wells

5　結(jié)　論

(1)提出一種補(bǔ)償密度測井快速正演模擬方法,引入空間響應(yīng)分布函數(shù)表征伽馬場分布及地層不同位置體積元對探測器計數(shù)貢獻(xiàn)。采用蒙特卡羅方法預(yù)先建立快速正演模擬所需數(shù)據(jù)庫,利用數(shù)據(jù)庫中的空間響應(yīng)分布函數(shù)計算不同地層區(qū)域?qū)μ綔y器計數(shù)貢獻(xiàn),并結(jié)合數(shù)據(jù)庫中單一密度地層中模擬記錄的探測器計數(shù),得到長、短源距探測器計數(shù),進(jìn)而根據(jù)數(shù)據(jù)庫中的轉(zhuǎn)換關(guān)系將探測器計數(shù)轉(zhuǎn)換為地層密度,可獲取高精度的補(bǔ)償密度測井快速正演模擬結(jié)果。

(2)本文快速正演模擬方法的計算精度與蒙特卡羅方法基本一致,但計算速度提高十萬倍以上,能夠滿足水平井井眼軌跡與地層關(guān)系確定中快速正演計算的需求。

(3)補(bǔ)償密度測井快速正演模擬結(jié)果的計算精度和計算速度與實際測井資料基本一致,驗證了新方法的實用性。