基于三維空間響應分布函數的脈沖中子孔隙度測井快速模擬

付新月 吳文圣 張恒榮 葛云龍 王 虎

1.中國石油大學（北京）油氣資源與探測國家重點實驗室 2.中海石油（中國）有限公司湛江分公司3.中國石油集團測井有限公司

0 引言

玻爾茲曼方程通常用來描述粒子在地層中的擴散運動。密度測井和中子孔隙度測井的物理過程可以用玻爾茲曼方程來描述，通過求解玻爾茲曼方程可以得到探測器的伽馬和中子通量[1-2]。然而直接用數學方法求解玻爾茲曼方程非常復雜，通常都是通過增加先驗條件來求解該方程[3]。與直接求解方法不同，蒙特卡洛方法可以通過跟蹤每個粒子的運動歷史來模擬探測器的粒子通量[4]。核測井響應的正演模擬大多用蒙特卡洛方法實現。這種方法精度高，可以設置三維地層和測井儀器模型并提供復雜的物質材料分布，能夠準確模擬得到探測器的中子通量并計算中子孔隙度值。對于中子孔隙度測井，該方法可以模擬中子測井響應，分析不同地層和環境對中子測井的影響[5-8]，同時也用于測井儀器的設計，如優化源距、屏蔽體、分析探測器靈敏度等[9-11]。然而這種方法存在計算速度較慢、需要大量計算機資源的缺陷。盡管計算機技術的最新進展大大減少了蒙特卡洛模擬所需的計算時間，但仍然無法滿足反演中所需要實時對比模擬測井曲線與實測曲線的要求[12]。

Watson[13]提出了用蒙特卡洛方法計算得到微分靈敏度方程來模擬巖性密度測井響應，成為快速模擬核測井響應的基礎。Couet等[14]用靈敏度函數方法分析了超熱中子測井儀器的空間響應，同時提出計算中子空間分布靈敏度函數的方法可以輔助測井儀器設計。Mendoza等[15-17]利用通量靈敏度函數和幾種近似方法來求解玻爾茲曼方程，實現了傳統化學源的核測井曲線的快速模擬。國外學者在核測井的快速模擬方法上做過一系列研究，國內對于快速模擬方法研究較少。葛云龍等[18-19]利用二維伽馬空間響應分布函數實現了補償密度測井曲線的快速模擬。目前國內在中子孔隙度測井曲線的快速模擬方面還沒有公開發表的文獻。同時在斜井的中子孔隙度測井中，不同方位地層不是各向同性，不同方位地層對探測器的貢獻不同[20]，利用二維空間響應分布函數來模擬中子孔隙度測井曲線會帶來一定誤差。

為此，筆者提出了用中子的三維空間響應分布函數快速模擬脈沖中子孔隙度測井曲線。首先，利用蒙特卡洛方法計算單位地層對探測器計數的貢獻即探測器的三維空間響應分布函數，考慮不同方位地層對探測器計數的影響，同時計算出不同地層條件下探測器的中子計數，建立不同地層的三維空間響應分布函數和探測器中子計數的數據庫。通過單一地層不同探測器的中子計數與探測范圍內地層對探測器計數的貢獻疊加得到不同組合地層下探測器的中子計數，根據不同探測器的中子計數比與地層孔隙度的關系求取地層的中子孔隙度值，實現中子孔隙度曲線的快速正演模擬。

1 脈沖中子孔隙度測井快速模擬原理

1.1 探測器計數計算原理

在一定體積內中子的通量隨時間變化率等于產生率減去泄漏率和吸收率，可以用玻爾茲曼方程來描述這一過程，如式（1）所示。等式左邊第一項是粒子的泄漏率，第二項是吸收率；等式右邊第一項是粒子碰撞轉換率，右邊第二項是粒子產生率[3]。

式中Ω、Ω分別表示中子的初始角度和碰撞后角度，（°）；φ表示能量為E、角度為Ω的中子通量，n/(cm2·s)；r表示空間中某一點與源的距離，cm；E、E分別表示中子的初始能量和碰撞后能量，MeV；∑t表示中子散射截面，1/cm；S表示源在單位時間內發射的粒子數，n/s；C表示中子發生散射的概率。

通過求解式（1），探測器的計數可以表示為每個從源發射的粒子對探測器貢獻的總和[21]，如式（2）：

式中N（rR）表示探測器的計數，n；rs表示空間中源的位置，無量綱；rR表示空間中探測器的位置，無量綱；ψ表示位于rs處的放射源發射的中子在r處的通量，n/(cm2·s)；D表示rR處的探測器響應函數，與單位地層對探測器計數貢獻的重要性和中子散射截面有關。

對于脈沖中子孔隙度測井，脈沖中子源發射的快中子與原子核發生非彈性散射、彈性散射反應，快中子能量降低變為熱中子，探測器主要探測熱中子。定義中子源發射的快中子經過幾種作用后生成的熱中子對探測器計數的貢獻為中子的空間響應分布函數FN，如式（3）：

式中FN表示中子的空間響應分布函數，無量綱；w表示每個單位地層對探測器計數的重要性，即單位地層對探測器的貢獻率。

空間響應分布函數定義為不同位置單位地層對探測器計數的貢獻。利用蒙特卡洛方法計算得到不同位置單位地層的熱中子通量和單位地層生成的熱中子對不同探測器計數的重要性值[22-23]。利用空間響應分布函數結合參考地層探測器計數可以計算得到不同組合地層的探測器計數[18-19]。

如圖1所示，可以將探測器計數視為探測區域內多套地層（地層B1、B2…Bn）對于探測器計數貢獻的疊加。模擬時首先利用蒙特卡洛方法計算探測區域內各單一地層的近、遠探測器空間響應分布函數作為背景空間響應分布函數，并分別計算得到各單一地層的近、遠探測器熱中子計數；然后在每個地層內對背景空間響應分布函數進行積分，得到每個地層對近、遠探測器計數的貢獻量。分別用每個地層的貢獻量除以探測區域內所有地層貢獻量的和得到不同地層貢獻率，隨后用相應地層的貢獻率與模擬得到的單一地層近、遠中子計數相乘后疊加得到綜合地層的近、遠探測器的計數，如式（4）、（5）所示。

圖1 脈沖中子孔隙度測井示意圖

式中n表示地層數，套；B1、B2…Bn表示探測范圍內第1套到第n套地層；Nnear、Nfar分別表示探測范圍內近、遠探測器的中子計數，n；FNS、FNL分別表示近、遠探測器的空間響應分布函數，無量綱；分別表示模擬得到第n套單一地層的近、遠探測器的中子計數，n。

1.2 中子孔隙度計算原理

脈沖中子孔隙度測井是通過建立近、遠探測器熱中子計數比與地層孔隙度的關系來求取地層孔隙度。根據雙組擴散理論[2]，距離源r處的熱中子計數如式（6）：

式中Dt表示熱中子擴散系數，cm；Lf表示中子減速長度，cm；Lt表示擴散長度，cm。

近、遠探測器計數比如式（7）：

式中r1、r2分別表示近、遠探測器與源的距離，cm。

對于熱中子來說，擴散長度Lt很小，當r較大時，Lt可以忽略不計，因此式（7）可以寫為：

在儀器源距一定的情況下，中子計數比主要受減速長度影響，地層的減速長度主要與地層含氫指數有關?？紫吨械牧黧w如水和油含氫指數較高，因此地層孔隙度的變化對中子計數比影響很大?？梢酝ㄟ^建立中子計數比與孔隙度的刻度關系來計算地層中子孔隙度。

2 空間響應分布函數計算

地層的空間響應分布函數為單位地層的熱中子通量與單位地層對探測器計數的重要性的乘積。利用蒙特卡洛方法分別計算單位地層的熱中子通量與地層的重要性，進而得到不同地層的空間響應分布函數數據庫。

利用蒙特卡洛方法建立三維地層模型，設置地層為圓柱體，高度為110 cm、半徑為60 cm，井眼半徑設置為10 cm，井眼流體為淡水、密度為1.00 g/cm3，地層中烴類氣和油的密度分別為0.10 g/cm3和0.80 g/cm3。放射源采用脈沖中子源，發射能量為14 MeV，儀器包含兩個He-3探測器，探測器長、短源距分別為65 cm和38 cm，源與探測器之間采用理想屏蔽體。地層分為9個扇區，靠近儀器的地層每30°設置一個扇區，儀器后部的地層每60°設置一個扇區。地層柵元在徑向和縱向上的間隔為1 cm。三維模型如圖2所示，圖2-a為圓柱體地層的剖面、圖2-b為圓柱體地層的橫截面。

圖2 蒙特卡洛地層模型圖

設置地層為砂巖，地層孔隙度為10%，孔隙內為水，利用式（3）計算得到不同地層條件下的三維空間響應分布函數。計算得到的三維空間響應分布函數如圖3所示。圖3所示為空間響應分布函數沿x軸和井眼的切面，不同顏色代表單位地層對探測器計數貢獻的大小。圖3-a為遠探測器的空間響應分布函數，圖3-b為近探測器的空間響應分布函數。從圖3可以看出，中子探測器計數主要由探測器附近地層所貢獻，對近探測器有貢獻的地層分布范圍稍小，說明近探測器探測范圍較小、探測深度較淺，對遠探測器有貢獻的地層分布范圍廣、探測深度大。井周向不同方位的地層對探測器計數貢獻不同，探測器前部地層的貢獻最大，探測器后部地層對探測器計數貢獻很小。

圖3 三維空間響應分布函數圖

通過圖2的模型計算得到不同地層條件下單位地層的熱中子通量與重要性。設置地層為砂巖，地層孔隙度為0、5%、10%、20%、30%、40%，孔隙流體為淡水和甲烷氣體，模擬得到不同地層條件下的三維空間響應分布函數。不考慮不同方位角度地層的影響，將圓柱體橫切面地層視為一個整體，對不同地層條件下的空間響應分布函數在縱向和徑向上進行積分，結果如圖4所示。從圖中可以看出，當地層性質如孔隙度或孔隙內流體發生改變時，空間響應分布函數同時也發生變化?？紫督橘|為水的地層空間響應分布函數值較大，孔隙度變化，含水地層的空間響應分布函數變化量較大；孔隙介質為氣的地層空間響應分布函數值與純砂巖地層的空間響應分布函數值相接近，孔隙度變化，含氣地層空間響應分布函數變化量較小。不同地層條件下，空間響應分布函數的歸一化幾何因子也不同。利用不同的空間響應分布函數計算的近、遠探測器計數有區別。快速模擬計算前，需要先建立不同地層條件下的空間響應分布函數的數據庫以便于后續計算。

圖4 不同地層條件下的空間響應分布函數圖

3 中子孔隙度刻度公式

利用圖2中的地層模型，模擬地層設置為單一地層，改變地層的孔隙度，分別模擬得到不同孔隙度地層的熱中子計數比，建立中子計數比與孔隙度的刻度關系。設置地層為石灰巖，孔隙流體為淡水。地層孔隙度設定為0、5%、10%、20%、30%、40%。模擬中子計數比與孔隙度的關系如圖5所示。

圖5 淡水石灰巖中子孔隙度測井刻度關系圖

對于淡水石灰巖地層，擬合得到地層孔隙度與中子計數比的關系，如式（9）：

式中φ表示地層孔隙度；R表示中子計數比，無量綱。

計算中子孔隙度時，首先根據空間響應分布函數，利用式（4）、（5）計算得到近、遠探測器計數比，然后根據式（9）計算出脈沖中子孔隙度值。

4 脈沖中子孔隙度測井正演模擬結果

4.1 直井中單一巖性地層的正演模擬

設置儀器探測范圍內地層只有B1和B2地層，地層巖性均為砂巖且地層孔隙中為水，B2地層孔隙度為30%，B1地層孔隙度為5%，儀器自上而下運行，地層傾角為0°，以砂巖地層且孔隙度為10%時模擬得到的空間響應分布函數為參考，利用式（4）、（5）進行積分得到近、遠探測器的計數，在通過式（9）計算出脈沖中子孔隙度值，得到的模擬結果與蒙特卡洛方法模擬的結果進行對比如圖6所示，最大誤差為 2.5 PU。

圖6 直井快速模擬結果圖

4.2 斜井中單一巖性地層的正演模擬

斜井的地層模型如圖7-a所示。設置地層為含水砂巖地層，地層傾角為75°，地層孔隙度為10%，利用蒙特卡洛方法計算出傾斜角為75°的斜井的空間響應分布函數如圖7-b、c所示。從圖7中可以看出，斜井的空間響應分布函數與直井有較大差異，儀器后部地層對儀器的測量響應也有一定的影響。對于斜井的快速模擬需根據斜井的空間響應分布函數進行計算。

圖7 斜井模型及其三維空間響應分布函數圖

根據圖7-a所示地層模型，開展斜井模型的脈沖中子孔隙度測井快速模擬。如圖7-a所示，地層界面將儀器的探測區域分割成兩個部分，設置上部地層孔隙度為5%、下部地層孔隙度為30%，地層為砂巖且孔隙中均為淡水。模擬儀器自上到下運行，參考的空間響應分布函數為75°斜井中10%孔隙度的砂巖地層的空間響應分布函數，利用式（4）、（5）計算不同探測器的中子計數，再通過式（9）將中子計數比轉化為中子孔隙度值。在相同的地層模型下利用蒙特卡洛方法模擬得到長、短源距探測器的計數并利用式（9）轉化成中子孔隙度值。對比蒙特卡洛方法模擬結果與快速模擬結果，得到結果如圖8所示。

圖8 斜井快速模擬結果圖

從圖6、8中可以看出，不管是直井還是斜井，快速模擬結果與蒙特卡洛方法計算結果基本一致，其中直井最大誤差為2.5 PU，斜井最大誤差為3.0 PU。在20核的服務器中蒙特卡洛方法模擬需要12 h，而快速模擬方法只需要十幾秒，大大縮短了模擬時間。

4.3 多組地層模擬

設置每組地層厚度相等，地層巖性分別為砂巖、石灰巖、白云巖、泥巖混層，不同組地層模型的孔隙度不同，孔隙內流體為水、氣和油，具體地層模型如表1所示，其中地層序號1～11、12～22和23～33分別代表砂巖、石灰巖和白云巖層。利用快速模擬方法和蒙特卡洛方法分別模擬不同地層條件下直井和75°斜井的脈沖中子孔隙度測井曲線，并對比快速模擬方法和蒙特卡洛方法的計算結果。模擬結果如表2所示。從表2中可以看出，快速模擬結果與蒙特卡洛模擬結果基本一致。含水石灰巖地層中的中子孔隙度計算結果與真孔隙度相等，含水砂巖地層的中子孔隙度小于真孔隙度，含水白云巖地層的中子孔隙度大于真孔隙度。當孔隙中含氣時，由于存在挖掘效應，含氣地層的中子孔隙度接近于0 PU；當孔隙中含油時，由于油的含氫指數低于水，油層的中子孔隙度略小于真實孔隙度；由于泥巖含氫指數較高，泥巖層的中子孔隙度大于真實孔隙度。

表1 多組地層模型中物質組分和性質參數表

表2 不同巖性在直井、斜井中的地層快速模擬結果表

5 結論

1）利用中子的三維空間響應分布函數可以實現脈沖中子孔隙度測井曲線的快速模擬。通過蒙特卡洛方法模擬得到不同地層條件下的中子三維空間響應分布函數以及不同地層近、遠探測器的熱中子計數，建立一個包含三維空間響應分布函數和不同條件地層探測器計數的數據庫。

2）對空間響應分布函數在不同地層區域進行積分得到不同地層區域對探測器計數的貢獻。利用不同地層區域對探測器計數的貢獻，結合預先模擬得到的單一地層近、遠探測器熱中子計數，計算得到不同組合地層探測器的熱中子計數。然后，利用中子

計數比與孔隙度的關系得到地層的中子孔隙度曲線。3）與蒙特卡洛方法模擬結果相比，快速模擬方法結果在直井和斜井中最大誤差分別為2.5 PU和3.0 PU。對于多組地層，快速模擬方法結果與蒙特卡洛方法模擬結果基本一致，可以滿足中子孔隙度曲線正演需求。與蒙特卡洛方法花費時間幾十小時相比，快速模擬方法僅用十幾秒，大大提高了計算速度，為后續反演等應用打下了堅實的計算基礎。