基于微擾理論的伽馬-伽馬密度測井快速正演方法

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

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

0 引言

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

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

通常，蒙特卡羅數值方法因其模擬精度高與實際測量吻合程度好，被廣泛應用于核測井儀器響應獲取、圖版建立和核儀器結構參數設計.但因其計算速度問題，沒有辦法直接應用于實際測井數據處理之中.

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

1 伽馬密度快速計算方法

伽馬-伽馬密度測井是利用伽馬源照射地層介質，伽馬射線與介質發生康普頓散射來實現地層密度的測量.地層密度與探測器的能窗計數具有如下關系：

ρ=Aln(N)+B.

(1)

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

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

(2)

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

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

圖1 密度測井測量示意圖

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

(3)

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

當介質1變為介質2時，衰減系數由μ0變為μ0+Δμ，探測器的計數可以表示為

(5)

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

(7)

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

(8)

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

(9)

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

2 擾動靈敏度函數

擾動靈敏度函數是儀器探測范圍內，不同體積柵元對探測器的貢獻權重，在柱坐標系下可以表示為：

(10)

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

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

從圖3可知，近、遠探測器的擾動靈敏度函數高值集中于源和探測器附近區域的地層，且擾動靈敏度函數呈近似的對稱，同時源附近和探測器附近的擾動靈敏度函數與井軸呈一定角度，原因是密度測井儀器的源和探測器結構設計上的開槽設置.

圖3 近、遠探測器擾動靈敏度函數

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

圖4 遠探測器軸向擾動靈敏度函數

圖5 遠探測器徑向擾動靈敏度函數

3 快速計算結果對比

3.1 不同巖性地層的快速計算結果對比

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

圖6 快速計算結果與MCNP模擬結果對比

表1 單個數據點快速計算與MCNP模擬效率對比

3.2 不同地層條件下的密度快速計算與MCNP模擬結果對比

3.2.1 水平層狀地層

建立如圖7所示的不同密度及厚度的層狀水平地層模型，模型為半徑30 cm、長7.1 m的圓柱體，地層由多種密度的含水飽和砂巖構成，密度變化范圍為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給出了快速正演響應(FF)、MCNP模擬響應結果及地層模型界面參數.其中密度快速正演響應(FF)是選用2.155 g·cm-3、2.32 g·cm-3和2.5675 g·cm-3等多種砂巖地層作為基準地層，同時選擇對應基準地層的擾動靈敏度函數，通過公式(9)計算得到.

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

圖8 水平層狀地層的伽馬密度正演響應結果對比

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

3.2.2 傾斜地層

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

圖9 傾斜地層模型示意圖

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

圖10 傾斜地層的伽馬密度正演響應結果對比

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

4 結論

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

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

(2)水平層狀地層條件下，快速計算的密度響應與MCNP模擬結果吻合，快速計算密度的絕對誤差在0.003 g·cm-3以內.當地層厚度小于儀器的縱向分辨率時，計算的地層密度與模型參數有差異，且受圍巖和夾層的雙重影響，但不影響利用密度曲線半幅點確定的地層界面位置的準確性.

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