陳前， 張鋒,2*， 梁啟軒， 刁書淵， 田立立，

范繼林1， 王才志4， 劉英明4， 遆永周5

0 引言

隨著油氣勘探開發(fā)理論、裝備的發(fā)展，油氣勘探開發(fā)不斷向陸上深層-超深層、海域深水和非常規(guī)三大領(lǐng)域挺進(jìn)(鄒才能等，2012).隨之而來的是地層巖性復(fù)雜、孔隙致密、構(gòu)造多變、測量環(huán)境惡劣等一系列問題，以隨鉆測井?dāng)?shù)據(jù)作為基礎(chǔ)的地質(zhì)導(dǎo)向及實(shí)時地層評價對于隨鉆水平井評價至關(guān)重要.當(dāng)前地質(zhì)導(dǎo)向過程中主要采用電阻率、伽馬、密度等測井技術(shù)，利用實(shí)時測量的數(shù)據(jù)與先導(dǎo)模型下的快速測井響應(yīng)進(jìn)行對比，以修正地層先導(dǎo)模型并確定地層參數(shù)，其中伽馬密度測井信息可以實(shí)現(xiàn)地層巖性、孔隙度、地層傾角等參數(shù)的計(jì)算.但是由于儀器與地層的空間位置及地層介質(zhì)等的變化，常用的蒙特卡羅(MCNP)數(shù)值模擬由于受自身耗時長、且建模受規(guī)則幾何的限制，不能滿足現(xiàn)場實(shí)時解釋評價及地質(zhì)模型參數(shù)快速準(zhǔn)確獲取的需求(Briesmeister, 2000)，因此需要開展密度快速計(jì)算研究.

國內(nèi)外學(xué)者針對測井方法的快速正演進(jìn)行了大量研究(Abubakar et al., 2008; Zhou, 2008; 邵才瑞等，2013; Qin et al., 2017; 胡旭飛等，2018; Wang and Fan, 2019).專門伽馬密度測井的快速正演，國外學(xué)者的研究較多，Both和Hendricks(1984)、Miller等(1990)分別對柵元重要性評估及柵元重要性自動生成技術(shù)進(jìn)行了研究；Liu和Gardner(1997)提出了可用于中子和伽馬射線的自適應(yīng)精細(xì)網(wǎng)格重要性圖模擬方法；Watson(1984)利用蒙特卡羅模擬了微分響應(yīng)函數(shù)；Mendoza等(2007,2010a,b)和Ijasan等(2012, 2013)引入通量靈敏度函數(shù)(Flux Sensitivity Function, FSF)，利用靈敏度函數(shù)結(jié)合格林函數(shù)實(shí)現(xiàn)了水平井和大斜度井條件下密度、中子孔隙度等的快速計(jì)算.相對而言，國內(nèi)學(xué)者關(guān)于伽馬密度快速計(jì)算的研究較少，袁超等(2018)引入空間響應(yīng)分布函數(shù)表征伽馬場分布及地層不同體積元對探測器計(jì)數(shù)貢獻(xiàn)，通過建立響應(yīng)函數(shù)數(shù)據(jù)庫進(jìn)而實(shí)現(xiàn)了補(bǔ)償密度的快速計(jì)算.其他的主要利用蒙特卡羅數(shù)值模擬獲取伽馬密度等放射性測井響應(yīng)(Liu et al., 2019; Zhang et al., 2019).

通常，蒙特卡羅數(shù)值方法因其模擬精度高與實(shí)際測量吻合程度好，被廣泛應(yīng)用于核測井儀器響應(yīng)獲取、圖版建立和核儀器結(jié)構(gòu)參數(shù)設(shè)計(jì).但因其計(jì)算速度問題，沒有辦法直接應(yīng)用于實(shí)際測井?dāng)?shù)據(jù)處理之中.

本文基于時間獨(dú)立的玻耳茲曼輸運(yùn)方程，提出將微擾理論和蒙特卡羅模擬相結(jié)合的方法來實(shí)現(xiàn)伽馬-伽馬密度快速計(jì)算.在對三維儀器地層模型及密度響應(yīng)公式進(jìn)行合理簡化的基礎(chǔ)上，通過MCNP模擬建立不同地層密度條件下的擾動靈敏度函數(shù)數(shù)據(jù)庫，最終形成伽馬密度的快速計(jì)算方法，為后續(xù)中子、伽馬密度實(shí)時聯(lián)合反演提供基礎(chǔ)支撐.

1 伽馬密度快速計(jì)算方法

伽馬-伽馬密度測井是利用伽馬源照射地層介質(zhì)，伽馬射線與介質(zhì)發(fā)生康普頓散射來實(shí)現(xiàn)地層密度的測量.地層密度與探測器的能窗計(jì)數(shù)具有如下關(guān)系：

ρ=Aln(N)+B.

(1)

公式(1)中，ρ為地層密度，單位g·cm-3，A和B為常數(shù)系數(shù)，與密度儀器相關(guān)，N為遠(yuǎn)探測器的能窗計(jì)數(shù).

由于伽馬探測器類型不同，能量分辨率及探測效率存在差異，本研究以NaI伽馬探測器為例，經(jīng)地層介質(zhì)散射后的伽馬射線在探測器的響應(yīng)可以表示為:

(2)

公式(2)中，rR為源距，φ(rR,r,E,Ω)為不同地層位置r處到達(dá)探測器的伽馬通量，R(rR,r,E,Ω)為探測器響應(yīng)函數(shù)，與探測器類型和尺寸相關(guān)，與地層性質(zhì)無關(guān).到達(dá)探測器的伽馬射線通量φ(rR,r,E,Ω)取決于伽馬射線能量E、源距rR和地層衰減系數(shù)μ三個因子，對于給定的密度儀器，其源能量E和源距rR已知，故此時φ(rR,r,E,Ω)主要受地層伽馬衰減系數(shù)μ支配.

如圖1所示，由伽馬源發(fā)射的伽馬光子與地層介質(zhì)發(fā)生康普頓散射、光電效應(yīng)作用后探測器被探測到，其整個過程可以使用時間獨(dú)立的玻耳茲曼輸運(yùn)方程進(jìn)行描述.

圖1 密度測井測量示意圖

→EΩ)φ(r,E′,Ω′)+S(r,E,Ω),

(3)

公式(3)中，φ為能量角粒子通量，E為伽馬射線能量，Ω為粒子散射方向，r為位置向量，μ為總衰減系數(shù)，μc為康普頓衰減系數(shù)，S為放射性源強(qiáng)度.

當(dāng)介質(zhì)1變?yōu)榻橘|(zhì)2時，衰減系數(shù)由μ0變?yōu)棣?+Δμ，探測器的計(jì)數(shù)可以表示為

(5)

擾動靈敏度函數(shù)PSF主要反映基準(zhǔn)地層條件下地層和儀器信息，其含義表示地層介質(zhì)對伽馬探測器計(jì)數(shù)的貢獻(xiàn)權(quán)重.對于給定的密度儀器，PSF主要取決于基準(zhǔn)地層的衰減系數(shù)μ0，所以在快速計(jì)算之前需要模擬多種基準(zhǔn)地層條件下的擾動靈敏度函數(shù).

(7)

公式(7)中，C0，C1和C2為基準(zhǔn)地層相關(guān)的系數(shù)，可以通過最小二乘法擬合得到.由于伽馬-伽馬密度測井主要利用的是康普頓散射，故地層介質(zhì)衰減系數(shù)的變化可以近似為介質(zhì)康普頓衰減系數(shù)的變化，公式(7)可以改寫為：

(8)

聯(lián)立公式(1)和公式(8)，可以得到地層密度快速計(jì)算公式為：

(9)

從密度快速計(jì)算公式(9)中可以看出，最終計(jì)算的地層密度值主要取決于擾動靈敏度函數(shù)PSF和介質(zhì)康普頓衰減系數(shù)變化量Δμc，擾動靈敏度函數(shù)PSF模擬的準(zhǔn)確程度將直接影響伽馬-伽馬密度快速計(jì)算的精度.

2 擾動靈敏度函數(shù)

擾動靈敏度函數(shù)是儀器探測范圍內(nèi)，不同體積柵元對探測器的貢獻(xiàn)權(quán)重，在柱坐標(biāo)系下可以表示為：

(10)

以137Cs補(bǔ)償密度儀器結(jié)構(gòu)參數(shù)為基準(zhǔn)，建立如圖2所示的MCNP數(shù)值計(jì)算模型.其中，儀器直徑89 mm，近探測器尺寸為φ25.4×12.7 mm，遠(yuǎn)探測器尺寸為φ38.1×40 mm，且長短源距分別為18.5 cm和41 cm.地層為圓柱狀飽含水砂巖，直徑為80 cm，井眼中充填淡水，儀器貼井壁測量.為了保證足夠的靈敏度函數(shù)模擬精度，地層被劃分為厚度0.5 cm寬度0.25 cm的同心圓環(huán)狀柵元.模擬得到的近、遠(yuǎn)探測器擾動靈敏度函數(shù)如圖3所示.

圖2 密度儀器蒙特卡羅計(jì)算模型示意圖

從圖3可知，近、遠(yuǎn)探測器的擾動靈敏度函數(shù)高值集中于源和探測器附近區(qū)域的地層，且擾動靈敏度函數(shù)呈近似的對稱，同時源附近和探測器附近的擾動靈敏度函數(shù)與井軸呈一定角度，原因是密度測井儀器的源和探測器結(jié)構(gòu)設(shè)計(jì)上的開槽設(shè)置.

圖3 近、遠(yuǎn)探測器擾動靈敏度函數(shù)

圖4和圖5分別給出了不同基準(zhǔn)密度地層條件下的遠(yuǎn)探測器的軸向擾動靈敏度函數(shù)和徑向擾動靈敏度函數(shù).由圖4可知，在軸向上，擾動靈敏度函數(shù)在源附近和探測器附近出現(xiàn)極大值，且不同基準(zhǔn)密度地層的遠(yuǎn)探測軸向擾動靈敏度函數(shù)基本重合.從圖5中可以看出，徑向擾動靈敏度函數(shù)隨徑向距離的增大而先增加后減小，在距離井壁2.5 cm處達(dá)到最大值.除此之外，在距離井壁0～7.5 cm范圍內(nèi)，徑向擾動靈敏度函數(shù)隨基準(zhǔn)地層密度的增大而增大，但當(dāng)徑向距離大于7.5 cm時，徑向擾動靈敏度函數(shù)隨基準(zhǔn)地層密度的增大而減小.在實(shí)際快速計(jì)算過程中，靈敏度函數(shù)數(shù)據(jù)庫主要由孔隙度從0%到100%，間隔10%的砂巖構(gòu)成，在選取靈敏度函數(shù)進(jìn)行快速計(jì)算時遵循“就近”原則，選取距離地層密度較近的靈敏度函數(shù)來計(jì)算，理論上基準(zhǔn)地層的密度值與地層密度值的相對誤差不大于10%.

圖4 遠(yuǎn)探測器軸向擾動靈敏度函數(shù)

圖5 遠(yuǎn)探測器徑向擾動靈敏度函數(shù)

3 快速計(jì)算結(jié)果對比

3.1 不同巖性地層的快速計(jì)算結(jié)果對比

基于上述模擬得到的不同密度地層的康普頓靈敏度函數(shù)數(shù)據(jù)庫，以砂巖為基準(zhǔn)地層，針對砂巖、碳酸巖、白云巖三種巖性和純水，密度變化范圍為1.0～2.87 g·cm-3的地層進(jìn)行快速計(jì)算，結(jié)果如圖6a所示.Δρ為快速計(jì)算密度值與模型密度值的差值絕對值，其結(jié)果如圖6b所示.由圖6a可知，快速計(jì)算的密度值與模型密度值吻合較好，均分布于45°線左右.由表1可知，砂巖密度快速計(jì)算平均誤差Δρ為0.0032 g·cm-3以內(nèi)，白云巖和碳酸巖密度快速計(jì)算平均誤差Δρ分別為0.0091 g·cm-3和0.0077 g·cm-3，說明密度快速計(jì)算仍然受光電效應(yīng)的影響，同時快速計(jì)算的速度是MCNP模擬速度的106倍以上.

圖6 快速計(jì)算結(jié)果與MCNP模擬結(jié)果對比

表1 單個數(shù)據(jù)點(diǎn)快速計(jì)算與MCNP模擬效率對比

3.2 不同地層條件下的密度快速計(jì)算與MCNP模擬結(jié)果對比

3.2.1 水平層狀地層

建立如圖7所示的不同密度及厚度的層狀水平地層模型，模型為半徑30 cm、長7.1 m的圓柱體，地層由多種密度的含水飽和砂巖構(gòu)成，密度變化范圍為1.99 g·cm-3到2.5675 g·cm-3.含水砂巖夾層厚度分別為30 cm、30 cm、30 cm、30 cm、30 cm、100 cm、40 cm、40 cm和40 cm，且夾層傾角為0°.井眼中充填淡水，伽馬密度儀器貼井壁測量，不考慮井眼尺寸和泥漿侵入影響，模擬儀器上提測量的整個過程.圖8給出了快速正演響應(yīng)(FF)、MCNP模擬響應(yīng)結(jié)果及地層模型界面參數(shù).其中密度快速正演響應(yīng)(FF)是選用2.155 g·cm-3、2.32 g·cm-3和2.5675 g·cm-3等多種砂巖地層作為基準(zhǔn)地層，同時選擇對應(yīng)基準(zhǔn)地層的擾動靈敏度函數(shù)，通過公式(9)計(jì)算得到.

圖7 水平層狀地層模型示意圖

圖8 水平層狀地層的伽馬密度正演響應(yīng)結(jié)果對比

由圖8可知，MCNP模擬密度響應(yīng)結(jié)果與通過快速計(jì)算得到的密度響應(yīng)結(jié)果一致，快速計(jì)算與MCNP模擬的誤差在0.003 g·cm-3以內(nèi)，由于該儀器的遠(yuǎn)探測器縱向分辨率為40 cm，故在30 cm的夾層計(jì)算的密度值與模型參數(shù)存在差異.與圖6結(jié)果相比，多種密度和厚度的層狀水平地層密度響應(yīng)更加復(fù)雜，曲線變化更大，但利用快速正演密度響應(yīng)曲線的半幅點(diǎn)確定的地層界面位置與模型界面位置吻合.當(dāng)?shù)貙雍穸刃∮趦x器縱向分辨率時，測量的密度值受圍巖和夾層的密度影響會偏大或偏小，但利用曲線確定的界面位置不受影響.當(dāng)圍巖與夾層密度差值較小時，在界面位置密度曲線變化較緩；當(dāng)圍巖與夾層密度差值較大時，在界面位置密度曲線急劇變化.

3.2.2 傾斜地層

建立如圖9所示的不同密度、傾角和厚度的地層模型，模型為半徑30 cm、長7.1 m的圓柱體，地層由多種密度的含水飽和砂巖構(gòu)成，密度變化范圍為1.99 g·cm-3到2.65 g·cm-3.含水砂巖夾層厚度分別為30 cm、40 cm、50 cm、60 cm、30 cm和30 cm，對應(yīng)的夾層傾角分別為135°、135°、120°、120°、0°、45°和45°(向右為正方向).圖10給出了快速正演響應(yīng)(FF)、MCNP模擬響應(yīng)結(jié)果及地層模型界面參數(shù).其中密度快速正演響應(yīng)(FF)是選用2.155 g·cm-3、2.32 g·cm-3和2.5675 g·cm-3等多種砂巖地層作為基準(zhǔn)地層，同時選擇對應(yīng)基準(zhǔn)地層的擾動靈敏度函數(shù)，通過公式(9)計(jì)算得到.

圖9 傾斜地層模型示意圖

由圖10可知，不同傾角、密度和厚度的MCNP模擬密度響應(yīng)結(jié)果與通過快速計(jì)算得到的密度響應(yīng)結(jié)果一致，且二者計(jì)算的密度絕對誤差小于0.0035 g·cm-3.

圖10 傾斜地層的伽馬密度正演響應(yīng)結(jié)果對比

與層狀水平地層相比，傾斜地層的密度曲線在薄層位置呈現(xiàn)明顯的軸非對稱性，并且利用密度曲線半幅點(diǎn)確定的層界面位置與模型設(shè)置參數(shù)存在差別.當(dāng)?shù)貙觾A斜角度不同時，地層界面與水平方向的夾角越大，密度曲線變化越緩慢，且利用密度曲線半幅點(diǎn)確定的地層厚度明顯偏大.對于厚度小于儀器縱向分辨率40 cm的傾斜薄層，其測量的密度值與模型密度值差距最大會達(dá)到0.1 g·cm-3.同時對于不同傾角的夾層，其密度曲線達(dá)到最大值的位置不同，主要是由于響應(yīng)函數(shù)分布和地層幾何形態(tài)的綜合影響，傾斜夾層對探測器的貢獻(xiàn)達(dá)到最大值的時間要延后，傾角越大，在密度曲線上表現(xiàn)為推遲的距離越長.傾斜地層和水平地層的密度曲線差異特征可以為后續(xù)利用伽馬密度進(jìn)行地層參數(shù)反演提供基礎(chǔ).

4 結(jié)論

密度快速正演作為中子和密度聯(lián)合反演模型的輸入，直接決定了實(shí)時反演確定地層幾何結(jié)構(gòu)及屬性參數(shù)的準(zhǔn)確程度.本文從時間獨(dú)立的玻耳茲曼輸運(yùn)方程入手，通過微擾理論結(jié)合蒙特卡羅模擬計(jì)算，形成伽馬-伽馬密度快速正演方法.同時快速正演了不同地層密度、厚度及傾角條件下的連續(xù)地層密度響應(yīng)并和MCNP模擬結(jié)果進(jìn)行了對比.具體結(jié)論如下：

(1)砂巖、白云巖和碳酸巖的密度快速正演結(jié)果與MCNP模擬結(jié)果絕對誤差小于0.01 g·cm-3.以砂巖作為基準(zhǔn)地層時，白云巖和碳酸巖的密度快速正演誤差明顯增大，表明快速計(jì)算的密度仍然受光電效應(yīng)的影響.

(2)水平層狀地層條件下，快速計(jì)算的密度響應(yīng)與MCNP模擬結(jié)果吻合，快速計(jì)算密度的絕對誤差在0.003 g·cm-3以內(nèi).當(dāng)?shù)貙雍穸刃∮趦x器的縱向分辨率時，計(jì)算的地層密度與模型參數(shù)有差異，且受圍巖和夾層的雙重影響，但不影響利用密度曲線半幅點(diǎn)確定的地層界面位置的準(zhǔn)確性.

(3)傾斜地層條件下，快速計(jì)算的密度響應(yīng)與MCNP模擬結(jié)果吻合且絕對誤差小于0.0035 g·cm-3.與層狀水平地層相比，傾斜地層的密度曲線在薄層位置呈現(xiàn)明顯的軸非對稱性，并且利用密度曲線半幅點(diǎn)確定的層界面位置與模型設(shè)置參數(shù)存在差異.同時對于不同傾角的夾層，其密度曲線達(dá)到最大值的位置不同，傾斜地層和水平地層的密度曲線差異特征可以為后續(xù)利用伽馬密度進(jìn)行地層參數(shù)反演提供基礎(chǔ).