地表附近宇宙射線(xiàn)中子測(cè)量土壤含水量刻度函數(shù)研究

李 會(huì)，侯英偉，李德源

(1.中國(guó)輻射防護(hù)研究院 核與輻射前沿技術(shù)研究中心，太原 030006；2.哈爾濱工程大學(xué) 核安全與仿真中心重點(diǎn)學(xué)科實(shí)驗(yàn)室，哈爾濱 150001)

0 引 言

宇宙射線(xiàn)中子法采用宇宙射線(xiàn)快中子在土壤中的慢化效果來(lái)反演土壤含水量。該方法測(cè)量區(qū)域平均土壤含水量，探測(cè)半徑約300 m，測(cè)量深度為12～75 cm[1-3]。該方法為被動(dòng)、非侵入式的中尺度土壤水分測(cè)量方法。可以實(shí)時(shí)監(jiān)測(cè)，無(wú)需大量取樣和實(shí)驗(yàn)室繁瑣的分析過(guò)程。實(shí)際應(yīng)用中的關(guān)鍵問(wèn)題是精確建立地表附近宇宙射線(xiàn)中子強(qiáng)度與土壤水分含量之間的刻度函數(shù)關(guān)系。通常采用蒙特卡洛模擬方法對(duì)刻度函數(shù)進(jìn)行研究，刻度函數(shù)的形式基本一致，而模型建立的正確與否決定了刻度函數(shù)適用性。

Zreda[4]等首先提出利用宇宙射線(xiàn)中子測(cè)量土壤水分的方法，采用MCNP對(duì)該測(cè)量方法的測(cè)量范圍和測(cè)量深度進(jìn)行了研究。宇宙射線(xiàn)中子從8 km[6]高度處入射，然后在大氣中輸運(yùn)。這種源設(shè)置方式存在以下問(wèn)題：①高能中子的作用截面不確定性高達(dá)50%，主要跟元素和反應(yīng)類(lèi)型有關(guān)。②大氣頂層專(zhuān)用于模擬的中子能譜無(wú)意的忽略了二次粒子，如質(zhì)子、π介子和μ介子產(chǎn)生的中子。③盡管H(氫)是中子慢化的主要成分，但大氣中水蒸汽在模擬中常常被忽略。

Desilets[5]采用MCNP建立模型用彌散分布在土壤中的人造中子源代替宇宙射線(xiàn)中子，計(jì)算得到有效測(cè)量半徑約300 m[5,6]。該方法存在的主要問(wèn)題：①高能中子在土壤中的指數(shù)衰減跟土壤類(lèi)型和地理位置有關(guān)。②土壤中的中子能譜有待驗(yàn)證。蒸發(fā)中子是主要成分，但不是能譜的全部成分。③從大氣入射的中子能譜包含低能成分和部分蒸發(fā)中子。

K?hli[7]采用URANOS蒙特卡洛程序利用Sato總結(jié)的宇宙射線(xiàn)中子能譜計(jì)算方法，將中子源設(shè)置在2～42 m高度處進(jìn)行模擬，計(jì)算的有效測(cè)量區(qū)域半徑遠(yuǎn)小于300 m[7,8]。將中子源設(shè)置在2～42 m的空間上，探測(cè)層設(shè)置在2 m處，雖然有利于與之前Zreda[4]等的模擬結(jié)果進(jìn)行對(duì)比，但實(shí)際上宇宙射線(xiàn)中子探測(cè)器可以安裝在飛行器上或其他高度下，而且通常情況下，在地形崎嶇復(fù)雜的地方，源設(shè)置在2～42 m大氣中不合理，忽略了土壤對(duì)地表附近宇宙射線(xiàn)中子的影響。

針對(duì)傳統(tǒng)模型射線(xiàn)源設(shè)置不合理，模擬結(jié)果存在差異問(wèn)題，本文利用Geant4軟件建立模型模擬宇宙射線(xiàn)粒子在地球磁場(chǎng)和大氣中的輸運(yùn)過(guò)程，研究地表附近宇宙射線(xiàn)中子對(duì)土壤容積含水量的響應(yīng)，并給出中子強(qiáng)度與土壤含水量之間的函數(shù)關(guān)系式。

1 地表附近宇宙射線(xiàn)中子測(cè)量土壤含水量原理

宇宙射線(xiàn)中子測(cè)量土壤含水量實(shí)際應(yīng)用的關(guān)鍵問(wèn)題是如何建立土壤含水量θ與探測(cè)器的中子計(jì)數(shù)率N之間的關(guān)系。如圖1所示，要利用實(shí)時(shí)測(cè)量的宇宙射線(xiàn)中子注量率得到地表土壤含水量，需要3個(gè)過(guò)程：采集土壤樣品計(jì)算土壤平均含水量、采用蒙特卡洛方法計(jì)算得到刻度函數(shù)關(guān)系式、利用測(cè)量的中子注量率和土壤平均含水量代入刻度函數(shù)關(guān)系式中計(jì)算刻度函數(shù)中的N0，獲得完整的刻度函數(shù)(N0方法)。根據(jù)實(shí)時(shí)測(cè)量中子計(jì)數(shù)率N和刻度函數(shù)可計(jì)算實(shí)時(shí)土壤含水量θ[5]。

圖1 宇宙射線(xiàn)中子測(cè)量土壤含水量流程圖Fig. 1 Flowchart of the method for measuring the soil moisture content with cosmic ray neutrons

圖2 海拔600 m時(shí)中子注量率和能譜Fig.2 Neutron energy spectra at altitude 600 m

2 模型建立

2.1 入射宇宙射線(xiàn)中子能譜

本模型將宇宙射線(xiàn)中子源設(shè)置在海拔500～600 m高度的空間中，且中子源的角分布采用余弦分布。通過(guò)以下分析將宇宙射線(xiàn)源設(shè)置在500～600 m海拔更合理。利用基于Geant4的PLANETOCOSMICS軟件計(jì)算得到地表分別設(shè)為H2O(水)和SiO2(二氧化硅)時(shí)哈爾濱上空100～1 000 m海拔下(間隔100 m)宇宙射線(xiàn)中子能譜，統(tǒng)計(jì)不同海拔下中子注量率能譜差值的標(biāo)準(zhǔn)差。在海拔>600 m時(shí)，宇宙射線(xiàn)中子譜的差異不大；海拔600 m時(shí)，能譜差異如圖2。同時(shí)，根據(jù)Rosolem[9]等的計(jì)算結(jié)果，干空氣時(shí)能利用宇宙射線(xiàn)中子譜反映土壤含水量差異的最大高度為412 m，濕空氣時(shí)能利用宇宙射線(xiàn)中子譜反映土壤含水量差異的最大高度為265 m。故而本文采用500～600 m海拔處的宇宙射線(xiàn)中子譜作為入射中子譜合理。

本模型中根據(jù)PLANETOCOMISCS計(jì)算的600 m海拔向下的宇宙射線(xiàn)中子能譜和蒙特卡洛方法生成入射宇宙射線(xiàn)中子的能譜如圖3所示；中子入射方向的空間角分布按：F=cos(θ)概率函數(shù)生成，其中θ指入射方向與垂直地表方向的夾角，取值范圍為0°～90°。

圖3 本模型隨機(jī)生成的輸入譜Fig.3 The input spectrum randomly generated by this model

2.2 系統(tǒng)布局

土壤厚度設(shè)置為6 m，不考慮土壤垂直方向上的異質(zhì)；最大的幾何空間設(shè)置為2 km×2 km×1 km。本模型用于計(jì)算距離地表2 m處，土壤體積含水量分別為0%、1%、2%、3%、7%、10%、20%、30%、40%、50%和100%，空氣濕度為10 g/m3時(shí)，宇宙射線(xiàn)中子的輸運(yùn)過(guò)程。入射107個(gè)宇宙射線(xiàn)中子，中子源設(shè)置在500～600 m海拔的空間。

圖4 宇宙射線(xiàn)中子測(cè)量土壤水分含量模擬模型設(shè)置(單位：m)Fig.4 Simulation model setup for measuring soil moisture content by cosmic ray neutron

模型如圖4所示，探測(cè)層布置在2 m高度處，探測(cè)層厚度為50 cm，探測(cè)層認(rèn)為是該高度下的透明大氣即理想探測(cè)器。選擇2 m高度目的是與Desilets[5]和K?hli[7]等的模擬結(jié)果進(jìn)行對(duì)比。探測(cè)層厚度50 cm是參考COSMOS宇宙射線(xiàn)中子探測(cè)系統(tǒng)實(shí)際中子探測(cè)器的長(zhǎng)度。采用探測(cè)層作為靈敏體積記錄輸運(yùn)的中子而不是采用探測(cè)器，這樣能夠在相同的計(jì)算量下獲得更多的統(tǒng)計(jì)數(shù)據(jù)。采用材料為透明的大氣是考慮到若設(shè)置成中子探測(cè)器的3He(氦三)或BF3(三氟化硼)材料，中子會(huì)與探測(cè)器材料反應(yīng)，這會(huì)導(dǎo)致從2 m高度以上本要入射到土壤的中子被吸收，計(jì)算結(jié)果不準(zhǔn)確。

2.3 大氣、土壤與水

2017年11月20日哈爾濱地區(qū)在0～1 km范圍內(nèi)的干空氣大氣氣壓及密度采用NRLMSISE00模型進(jìn)行計(jì)算得到，如表1

表1 2017年11月20日哈爾濱地表以上0～1 km干空氣密度Tab.1 Dry air density of 0～1 km above the surface of Harbin on November 20,2017

所示。地面附近飽和空氣濕度<30 g/m3，本次計(jì)算中空氣濕度選為10 g/m3[7]。

土壤由土和可調(diào)比例的水組成；其中土體積分?jǐn)?shù)為50%，孔隙水的體積分?jǐn)?shù)為0～50%，孔隙的體積分?jǐn)?shù)為50%。參考K?hli[7]對(duì)土壤建模方式，土壤成分由SiO2和Al2O3(三氧化二鋁)組成，體積比為75% SiO2和25% Al2O3，密度為2.86 g/cm3；土壤考慮到孔隙和孔隙水后，土壤容積含水量由0%增加到50%，土壤總密度從1.43 g/cm3可變化到1.93 g/cm3；本文的建模過(guò)程中不考慮土壤化學(xué)成分的差異，只考慮土壤體積含水量的差異。

3 模擬結(jié)果及分析

3.1 地表背散射中子能譜驗(yàn)證

地表背散射中子是指入射到土壤的宇宙射線(xiàn)中子與土壤相互作用后散射離開(kāi)地表回到大氣的宇宙射線(xiàn)中子；為與URANOS計(jì)算的地表背散射能譜歸一化后進(jìn)行比較，本文利用Geant4模擬土壤體積含水量為1%，空氣濕度為0.33 g/m3時(shí)中子的輸運(yùn)過(guò)程，統(tǒng)計(jì)所有離開(kāi)地表向上運(yùn)動(dòng)中子的動(dòng)能得到了圖5中紅色的中子能譜曲線(xiàn)；計(jì)算結(jié)果表明，Geant4計(jì)算的地表背散射宇宙射線(xiàn)中子能譜與URANOS計(jì)算結(jié)果共同點(diǎn)是級(jí)聯(lián)中子峰不明顯，蒸發(fā)中子峰明顯，不同點(diǎn)是Geant4計(jì)算的結(jié)果中熱中子峰明顯，慢化中子強(qiáng)度稍小。

圖5 Geant4計(jì)算得到的地表背散射中子譜與URANOS計(jì)算的背散射宇宙射線(xiàn)中子譜歸一化后的比較Fig.5 Comparison of normalized backscatter neutron spectra of cosmic rays calculated by Geant4 and URANOS

理論上，宇宙射線(xiàn)中子與土壤原子相互作用后，會(huì)有熱中子產(chǎn)生。由于H對(duì)中子慢化做主要貢獻(xiàn)，而Geant4計(jì)算的熱中子峰明顯，說(shuō)明Geant4計(jì)算的結(jié)果更準(zhǔn)確。Geant4計(jì)算中使用H的中子反應(yīng)截面與URANOS使用的不同；URANOS軟件中1H(氕)和2H(氘)的吸收截面和散射截面均采用的是2007年版的日本的評(píng)價(jià)數(shù)據(jù)庫(kù)JENDL-4.0，而Geant4在中子能量大于20 MeV時(shí)采用的是日本的評(píng)價(jià)數(shù)據(jù)庫(kù)JENDF，中子能量小于20 MeV時(shí)用的是美國(guó)的評(píng)價(jià)數(shù)據(jù)庫(kù)ENDF/B-VII。而JENDL中子截面庫(kù)在對(duì)高能區(qū)的中子輸運(yùn)計(jì)算比較準(zhǔn)確，在對(duì)低能區(qū)的中子輸運(yùn)計(jì)算并不如ENDF/B-VII計(jì)算準(zhǔn)確。而且在URANOS的計(jì)算中只考慮了中子的吸收截面和散射截面，在本模型計(jì)算中考慮了中子散射、吸收和裂變，使用了高精度HP物理模型。以上這些導(dǎo)致了基于Geant4的本模型計(jì)算的結(jié)果更準(zhǔn)確，在0.025eV左右的熱中子峰明顯存在。故而，本模型計(jì)算結(jié)果比URANOS計(jì)算結(jié)果更準(zhǔn)確。

3.2 地表附近宇宙射線(xiàn)中子能譜對(duì)土壤水分響應(yīng)

統(tǒng)計(jì)進(jìn)入到探測(cè)層的宇宙射線(xiàn)中子動(dòng)能，得到如圖6所示的土壤體積含水量不同時(shí)，距地表2 m處的宇宙射線(xiàn)中子能譜。圖中縱坐標(biāo)為計(jì)數(shù)N，通過(guò)除以入射宇宙射線(xiàn)中子數(shù)107，即可得到相對(duì)中子強(qiáng)度。計(jì)算結(jié)果表明，H對(duì)中子慢化起主要貢獻(xiàn)；土壤水分體積含量增加時(shí)，慢化中子強(qiáng)度明顯減少，熱中子強(qiáng)度增加，蒸發(fā)中子峰降低，級(jí)聯(lián)中子峰強(qiáng)度基本不變。出現(xiàn)該現(xiàn)象是由于土壤含水量的增加，被土壤元素慢化和吸收中子數(shù)增加，導(dǎo)致大氣中蒸發(fā)中子峰、慢化中子峰和熱中子峰的總強(qiáng)度要減少。H對(duì)中子起主要慢化作用，當(dāng)H含量增加后，由于與熱中子的散射截面大，吸收截面小，熱中子強(qiáng)度相對(duì)慢化中子強(qiáng)度變大。

圖6 地表以上2 m高處土壤水分含量不同時(shí)，本模型計(jì)算的宇宙射線(xiàn)中子強(qiáng)度能譜Fig.6 The differences of cosmic ray neutron spectra at 2 m height from the earth surface for the different soil moisture contents

圖6中結(jié)果表明10-6～10-2MeV慢化中子對(duì)土壤水分含量有較好的分辨能力，慢化能區(qū)的中子強(qiáng)度隨土壤含水量的增加而減少。同時(shí)，在熱中子0.5×10-6MeV以下的中子對(duì)土壤水分在一定范圍內(nèi)變化時(shí)也具有分辨能力。目前的COSMOS系統(tǒng)[1,10]利用熱中子探測(cè)器與快中子探測(cè)器的計(jì)數(shù)率比值來(lái)確定地表積雪或植被水當(dāng)量，從而將測(cè)量的總水分含量減去地表積雪或植被水含量，更精確的得到土壤總含水量。

3.3 刻度函數(shù)的計(jì)算

為與Desilets[5]利用MCNP模擬計(jì)算的刻度函數(shù)關(guān)系式進(jìn)行對(duì)比，統(tǒng)計(jì)能量在10-8～10-4MeV范圍中子計(jì)數(shù)，該能量段包含了熱中子和部分快中子。當(dāng)土壤體積含水量為0%的干土壤時(shí)，探測(cè)層的計(jì)數(shù)為N0，其他土壤容積含水量下的，探測(cè)層的計(jì)數(shù)為N(假定探測(cè)層對(duì)所有能量的中子探測(cè)效率η=100%)。Geant4模擬計(jì)算得到的探測(cè)層計(jì)數(shù)如表2所示。

將容積含水量θvol和質(zhì)量含水量θw與N/N0一一對(duì)應(yīng)點(diǎn)繪制在圖7中，并與Desilets[5]利用MCNP軟件計(jì)算得到擬合形狀函數(shù)(圖中紅色曲線(xiàn))進(jìn)行對(duì)比。采用與Desilets[5]模擬的N0方法刻度函數(shù)類(lèi)似形式進(jìn)行擬合，本模型計(jì)算結(jié)果的擬合函數(shù)如式(1)所示。

表2 土壤含水量不同時(shí)，探測(cè)層記錄的10-8～10-4 MeV能區(qū)的中子計(jì)數(shù)Tab.2 Neutrons in the 10-8～10-4 MeV energy region recorded by the detection layer for the different soil moisture contents

(1)

當(dāng)θ為θvol時(shí)，擬合參數(shù)a1=0.028 632 5，a2=0.567 39,a3=0.399 761；

當(dāng)θ為θw時(shí)a1=0.021 725，a2=0.041 394,a3=0.378 991。

圖7 土壤含水量與N/N0的關(guān)系Fig.7 N/N0 as a function of soil moisture

顯然，采用Geant4計(jì)算得到的擬合函數(shù)與Desilets[5]采用MCNP計(jì)算得到的擬合函數(shù)比較接近，且對(duì)干燥土壤含水量變化更為靈敏。采用MCNP計(jì)算給出的是形狀函數(shù)，同時(shí)適用于土壤容積含水量和質(zhì)量含水量。Desilets在建模過(guò)程中用土壤中的彌散人造中子源近似實(shí)際入射到土壤的宇宙射線(xiàn)中子，這導(dǎo)致計(jì)算結(jié)果不準(zhǔn)確，前言中已有說(shuō)明，在此不再贅述。在實(shí)際應(yīng)用過(guò)程中，通常會(huì)對(duì)Desiltes[5]給出的刻度函數(shù)進(jìn)行修正，要么采取對(duì)N進(jìn)行修正，要么對(duì)θ進(jìn)行修正。

前文為了與Desilets計(jì)算結(jié)果對(duì)比，故而中子計(jì)數(shù)能量區(qū)間選擇了10-8～10-4MeV。根據(jù)圖6所示，顯然能量區(qū)間10-6～10-2MeV中子計(jì)數(shù)與土壤含水量的響應(yīng)關(guān)系更直觀。不同能量區(qū)間的中子計(jì)數(shù)與土壤含水量的刻度函數(shù)擬合曲線(xiàn)如圖8所示，兩能量區(qū)間刻度曲線(xiàn)隨著土壤含水量的增大差異逐漸增大。導(dǎo)致差異的主要原因是10-8～10-4MeV包含了熱中子。然而熱中子能量區(qū)間對(duì)土壤含水量變換并非線(xiàn)性響應(yīng)，故能量區(qū)間10-6～10-2MeV的刻度函數(shù)結(jié)果具有更好的靈敏度。其擬合參數(shù)值為：

當(dāng)θ為θvol時(shí)，a1=0.021 597;a2=0.172 635;a3=0.028 315；

當(dāng)θ為θw時(shí)，a1=0.028 511;a2=0.399 549;a3=0.056 440。

圖8 不同中子能區(qū)土壤含水量與N/N0的關(guān)系Fig.8 N/N0 as a function of the soil moisture for the neutrons within different energy ranges

4 結(jié) 論

本文針對(duì)常用的基于MCNP和URANOS的土壤含水量刻度函數(shù)計(jì)算模型中中子源設(shè)置不合理的問(wèn)題，采用Geant4建立了新的刻度函數(shù)計(jì)算模型，中子源的入射高度和能譜是本模型不同于傳統(tǒng)模型之處。在不同地表成分和地理位置以及海拔情形下，采用宇宙射線(xiàn)粒子輸運(yùn)程序PLANETOCOMICS分別計(jì)算了宇宙射線(xiàn)中子能譜并分析差異，發(fā)現(xiàn)地表對(duì)海拔600 m及以上的宇宙射線(xiàn)中子強(qiáng)度無(wú)影響，中子源設(shè)置在600 m海拔處，計(jì)算結(jié)果優(yōu)于傳統(tǒng)MCNP和URANOS模型的計(jì)算結(jié)果，且能量區(qū)間10-6～10-2MeV的中子對(duì)土壤水分的響應(yīng)比10-8～10-4MeV能量區(qū)間更靈敏?；谝陨辖Y(jié)論給出了合適能量區(qū)間內(nèi)的刻度函數(shù)擬合形式和擬合參數(shù)，結(jié)合圖1流程計(jì)算N0可用于宇宙射線(xiàn)中子實(shí)時(shí)測(cè)量土壤含水量。