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

1.中國(guó)石油大學(xué)（北京）油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室 2.中海石油（中國(guó)）有限公司湛江分公司3.中國(guó)石油集團(tuán)測(cè)井有限公司

0 引言

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

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

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

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

1.1 探測(cè)器計(jì)數(shù)計(jì)算原理

在一定體積內(nèi)中子的通量隨時(shí)間變化率等于產(chǎn)生率減去泄漏率和吸收率，可以用玻爾茲曼方程來(lái)描述這一過(guò)程，如式（1）所示。等式左邊第一項(xiàng)是粒子的泄漏率，第二項(xiàng)是吸收率；等式右邊第一項(xiàng)是粒子碰撞轉(zhuǎn)換率，右邊第二項(xiàng)是粒子產(chǎn)生率[3]。

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

通過(guò)求解式（1），探測(cè)器的計(jì)數(shù)可以表示為每個(gè)從源發(fā)射的粒子對(duì)探測(cè)器貢獻(xiàn)的總和[21]，如式（2）：

式中N（rR）表示探測(cè)器的計(jì)數(shù)，n；rs表示空間中源的位置，無(wú)量綱；rR表示空間中探測(cè)器的位置，無(wú)量綱；ψ表示位于rs處的放射源發(fā)射的中子在r處的通量，n/(cm2·s)；D表示rR處的探測(cè)器響應(yīng)函數(shù)，與單位地層對(duì)探測(cè)器計(jì)數(shù)貢獻(xiàn)的重要性和中子散射截面有關(guān)。

對(duì)于脈沖中子孔隙度測(cè)井，脈沖中子源發(fā)射的快中子與原子核發(fā)生非彈性散射、彈性散射反應(yīng)，快中子能量降低變?yōu)闊嶂凶?，探測(cè)器主要探測(cè)熱中子。定義中子源發(fā)射的快中子經(jīng)過(guò)幾種作用后生成的熱中子對(duì)探測(cè)器計(jì)數(shù)的貢獻(xiàn)為中子的空間響應(yīng)分布函數(shù)FN，如式（3）：

式中FN表示中子的空間響應(yīng)分布函數(shù)，無(wú)量綱；w表示每個(gè)單位地層對(duì)探測(cè)器計(jì)數(shù)的重要性，即單位地層對(duì)探測(cè)器的貢獻(xiàn)率。

空間響應(yīng)分布函數(shù)定義為不同位置單位地層對(duì)探測(cè)器計(jì)數(shù)的貢獻(xiàn)。利用蒙特卡洛方法計(jì)算得到不同位置單位地層的熱中子通量和單位地層生成的熱中子對(duì)不同探測(cè)器計(jì)數(shù)的重要性值[22-23]。利用空間響應(yīng)分布函數(shù)結(jié)合參考地層探測(cè)器計(jì)數(shù)可以計(jì)算得到不同組合地層的探測(cè)器計(jì)數(shù)[18-19]。

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

圖1 脈沖中子孔隙度測(cè)井示意圖

式中n表示地層數(shù)，套；B1、B2…Bn表示探測(cè)范圍內(nèi)第1套到第n套地層；Nnear、Nfar分別表示探測(cè)范圍內(nèi)近、遠(yuǎn)探測(cè)器的中子計(jì)數(shù)，n；FNS、FNL分別表示近、遠(yuǎn)探測(cè)器的空間響應(yīng)分布函數(shù)，無(wú)量綱；分別表示模擬得到第n套單一地層的近、遠(yuǎn)探測(cè)器的中子計(jì)數(shù)，n。

1.2 中子孔隙度計(jì)算原理

脈沖中子孔隙度測(cè)井是通過(guò)建立近、遠(yuǎn)探測(cè)器熱中子計(jì)數(shù)比與地層孔隙度的關(guān)系來(lái)求取地層孔隙度。根據(jù)雙組擴(kuò)散理論[2]，距離源r處的熱中子計(jì)數(shù)如式（6）：

式中Dt表示熱中子擴(kuò)散系數(shù)，cm；Lf表示中子減速長(zhǎng)度，cm；Lt表示擴(kuò)散長(zhǎng)度，cm。

近、遠(yuǎn)探測(cè)器計(jì)數(shù)比如式（7）：

式中r1、r2分別表示近、遠(yuǎn)探測(cè)器與源的距離，cm。

對(duì)于熱中子來(lái)說(shuō)，擴(kuò)散長(zhǎng)度Lt很小，當(dāng)r較大時(shí)，Lt可以忽略不計(jì)，因此式（7）可以寫(xiě)為：

在儀器源距一定的情況下，中子計(jì)數(shù)比主要受減速長(zhǎng)度影響，地層的減速長(zhǎng)度主要與地層含氫指數(shù)有關(guān)?？紫吨械牧黧w如水和油含氫指數(shù)較高，因此地層孔隙度的變化對(duì)中子計(jì)數(shù)比影響很大。可以通過(guò)建立中子計(jì)數(shù)比與孔隙度的刻度關(guān)系來(lái)計(jì)算地層中子孔隙度。

2 空間響應(yīng)分布函數(shù)計(jì)算

地層的空間響應(yīng)分布函數(shù)為單位地層的熱中子通量與單位地層對(duì)探測(cè)器計(jì)數(shù)的重要性的乘積。利用蒙特卡洛方法分別計(jì)算單位地層的熱中子通量與地層的重要性，進(jìn)而得到不同地層的空間響應(yīng)分布函數(shù)數(shù)據(jù)庫(kù)。

利用蒙特卡洛方法建立三維地層模型，設(shè)置地層為圓柱體，高度為110 cm、半徑為60 cm，井眼半徑設(shè)置為10 cm，井眼流體為淡水、密度為1.00 g/cm3，地層中烴類(lèi)氣和油的密度分別為0.10 g/cm3和0.80 g/cm3。放射源采用脈沖中子源，發(fā)射能量為14 MeV，儀器包含兩個(gè)He-3探測(cè)器，探測(cè)器長(zhǎng)、短源距分別為65 cm和38 cm，源與探測(cè)器之間采用理想屏蔽體。地層分為9個(gè)扇區(qū)，靠近儀器的地層每30°設(shè)置一個(gè)扇區(qū)，儀器后部的地層每60°設(shè)置一個(gè)扇區(qū)。地層?xùn)旁趶较蚝涂v向上的間隔為1 cm。三維模型如圖2所示，圖2-a為圓柱體地層的剖面、圖2-b為圓柱體地層的橫截面。

圖2 蒙特卡洛地層模型圖

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

圖3 三維空間響應(yīng)分布函數(shù)圖

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

圖4 不同地層條件下的空間響應(yīng)分布函數(shù)圖

3 中子孔隙度刻度公式

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

圖5 淡水石灰?guī)r中子孔隙度測(cè)井刻度關(guān)系圖

對(duì)于淡水石灰?guī)r地層，擬合得到地層孔隙度與中子計(jì)數(shù)比的關(guān)系，如式（9）：

式中φ表示地層孔隙度；R表示中子計(jì)數(shù)比，無(wú)量綱。

計(jì)算中子孔隙度時(shí)，首先根據(jù)空間響應(yīng)分布函數(shù)，利用式（4）、（5）計(jì)算得到近、遠(yuǎn)探測(cè)器計(jì)數(shù)比，然后根據(jù)式（9）計(jì)算出脈沖中子孔隙度值。

4 脈沖中子孔隙度測(cè)井正演模擬結(jié)果

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

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

圖6 直井快速模擬結(jié)果圖

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

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

圖7 斜井模型及其三維空間響應(yīng)分布函數(shù)圖

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

圖8 斜井快速模擬結(jié)果圖

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

4.3 多組地層模擬

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

表1 多組地層模型中物質(zhì)組分和性質(zhì)參數(shù)表

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

5 結(jié)論

1）利用中子的三維空間響應(yīng)分布函數(shù)可以實(shí)現(xiàn)脈沖中子孔隙度測(cè)井曲線的快速模擬。通過(guò)蒙特卡洛方法模擬得到不同地層條件下的中子三維空間響應(yīng)分布函數(shù)以及不同地層近、遠(yuǎn)探測(cè)器的熱中子計(jì)數(shù)，建立一個(gè)包含三維空間響應(yīng)分布函數(shù)和不同條件地層探測(cè)器計(jì)數(shù)的數(shù)據(jù)庫(kù)。

2）對(duì)空間響應(yīng)分布函數(shù)在不同地層區(qū)域進(jìn)行積分得到不同地層區(qū)域?qū)μ綔y(cè)器計(jì)數(shù)的貢獻(xiàn)。利用不同地層區(qū)域?qū)μ綔y(cè)器計(jì)數(shù)的貢獻(xiàn)，結(jié)合預(yù)先模擬得到的單一地層近、遠(yuǎn)探測(cè)器熱中子計(jì)數(shù)，計(jì)算得到不同組合地層探測(cè)器的熱中子計(jì)數(shù)。然后，利用中子

計(jì)數(shù)比與孔隙度的關(guān)系得到地層的中子孔隙度曲線。3）與蒙特卡洛方法模擬結(jié)果相比，快速模擬方法結(jié)果在直井和斜井中最大誤差分別為2.5 PU和3.0 PU。對(duì)于多組地層，快速模擬方法結(jié)果與蒙特卡洛方法模擬結(jié)果基本一致，可以滿(mǎn)足中子孔隙度曲線正演需求。與蒙特卡洛方法花費(fèi)時(shí)間幾十小時(shí)相比，快速模擬方法僅用十幾秒，大大提高了計(jì)算速度，為后續(xù)反演等應(yīng)用打下了堅(jiān)實(shí)的計(jì)算基礎(chǔ)。