地表附近宇宙射线中子测量土壤含水量刻度函数研究

李 会，侯英伟，李德源

(1.中国辐射防护研究院 核与辐射前沿技术研究中心，太原 030006；2.哈尔滨工程大学 核安全与仿真中心重点学科实验室，哈尔滨 150001)

0 引 言

宇宙射线中子法采用宇宙射线快中子在土壤中的慢化效果来反演土壤含水量。该方法测量区域平均土壤含水量，探测半径约300 m，测量深度为12～75 cm[1-3]。该方法为被动、非侵入式的中尺度土壤水分测量方法。可以实时监测，无需大量取样和实验室繁琐的分析过程。实际应用中的关键问题是精确建立地表附近宇宙射线中子强度与土壤水分含量之间的刻度函数关系。通常采用蒙特卡洛模拟方法对刻度函数进行研究，刻度函数的形式基本一致，而模型建立的正确与否决定了刻度函数适用性。

Zreda[4]等首先提出利用宇宙射线中子测量土壤水分的方法，采用MCNP对该测量方法的测量范围和测量深度进行了研究。宇宙射线中子从8 km[6]高度处入射，然后在大气中输运。这种源设置方式存在以下问题：①高能中子的作用截面不确定性高达50%，主要跟元素和反应类型有关。②大气顶层专用于模拟的中子能谱无意的忽略了二次粒子，如质子、π介子和μ介子产生的中子。③尽管H(氢)是中子慢化的主要成分，但大气中水蒸汽在模拟中常常被忽略。

Desilets[5]采用MCNP建立模型用弥散分布在土壤中的人造中子源代替宇宙射线中子，计算得到有效测量半径约300 m[5,6]。该方法存在的主要问题：①高能中子在土壤中的指数衰减跟土壤类型和地理位置有关。②土壤中的中子能谱有待验证。蒸发中子是主要成分，但不是能谱的全部成分。③从大气入射的中子能谱包含低能成分和部分蒸发中子。

Köhli[7]采用URANOS蒙特卡洛程序利用Sato总结的宇宙射线中子能谱计算方法，将中子源设置在2～42 m高度处进行模拟，计算的有效测量区域半径远小于300 m[7,8]。将中子源设置在2～42 m的空间上，探测层设置在2 m处，虽然有利于与之前Zreda[4]等的模拟结果进行对比，但实际上宇宙射线中子探测器可以安装在飞行器上或其他高度下，而且通常情况下，在地形崎岖复杂的地方，源设置在2～42 m大气中不合理，忽略了土壤对地表附近宇宙射线中子的影响。

针对传统模型射线源设置不合理，模拟结果存在差异问题，本文利用Geant4软件建立模型模拟宇宙射线粒子在地球磁场和大气中的输运过程，研究地表附近宇宙射线中子对土壤容积含水量的响应，并给出中子强度与土壤含水量之间的函数关系式。

1 地表附近宇宙射线中子测量土壤含水量原理

宇宙射线中子测量土壤含水量实际应用的关键问题是如何建立土壤含水量θ与探测器的中子计数率N之间的关系。如图1所示，要利用实时测量的宇宙射线中子注量率得到地表土壤含水量，需要3个过程：采集土壤样品计算土壤平均含水量、采用蒙特卡洛方法计算得到刻度函数关系式、利用测量的中子注量率和土壤平均含水量代入刻度函数关系式中计算刻度函数中的N0，获得完整的刻度函数(N0方法)。根据实时测量中子计数率N和刻度函数可计算实时土壤含水量θ[5]。

图1 宇宙射线中子测量土壤含水量流程图Fig. 1 Flowchart of the method for measuring the soil moisture content with cosmic ray neutrons

图2 海拔600 m时中子注量率和能谱Fig.2 Neutron energy spectra at altitude 600 m

2 模型建立

2.1 入射宇宙射线中子能谱

本模型将宇宙射线中子源设置在海拔500～600 m高度的空间中，且中子源的角分布采用余弦分布。通过以下分析将宇宙射线源设置在500～600 m海拔更合理。利用基于Geant4的PLANETOCOSMICS软件计算得到地表分别设为H2O(水)和SiO2(二氧化硅)时哈尔滨上空100～1 000 m海拔下(间隔100 m)宇宙射线中子能谱，统计不同海拔下中子注量率能谱差值的标准差。在海拔>600 m时，宇宙射线中子谱的差异不大；海拔600 m时，能谱差异如图2。同时，根据Rosolem[9]等的计算结果，干空气时能利用宇宙射线中子谱反映土壤含水量差异的最大高度为412 m，湿空气时能利用宇宙射线中子谱反映土壤含水量差异的最大高度为265 m。故而本文采用500～600 m海拔处的宇宙射线中子谱作为入射中子谱合理。

本模型中根据PLANETOCOMISCS计算的600 m海拔向下的宇宙射线中子能谱和蒙特卡洛方法生成入射宇宙射线中子的能谱如图3所示；中子入射方向的空间角分布按：F=cos(θ)概率函数生成，其中θ指入射方向与垂直地表方向的夹角，取值范围为0°～90°。

图3 本模型随机生成的输入谱Fig.3 The input spectrum randomly generated by this model

2.2 系统布局

土壤厚度设置为6 m，不考虑土壤垂直方向上的异质；最大的几何空间设置为2 km×2 km×1 km。本模型用于计算距离地表2 m处，土壤体积含水量分别为0%、1%、2%、3%、7%、10%、20%、30%、40%、50%和100%，空气湿度为10 g/m3时，宇宙射线中子的输运过程。入射107个宇宙射线中子，中子源设置在500～600 m海拔的空间。

图4 宇宙射线中子测量土壤水分含量模拟模型设置(单位：m)Fig.4 Simulation model setup for measuring soil moisture content by cosmic ray neutron

模型如图4所示，探测层布置在2 m高度处，探测层厚度为50 cm，探测层认为是该高度下的透明大气即理想探测器。选择2 m高度目的是与Desilets[5]和Köhli[7]等的模拟结果进行对比。探测层厚度50 cm是参考COSMOS宇宙射线中子探测系统实际中子探测器的长度。采用探测层作为灵敏体积记录输运的中子而不是采用探测器，这样能够在相同的计算量下获得更多的统计数据。采用材料为透明的大气是考虑到若设置成中子探测器的3He(氦三)或BF3(三氟化硼)材料，中子会与探测器材料反应，这会导致从2 m高度以上本要入射到土壤的中子被吸收，计算结果不准确。

2.3 大气、土壤与水

2017年11月20日哈尔滨地区在0～1 km范围内的干空气大气气压及密度采用NRLMSISE00模型进行计算得到，如表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，本次计算中空气湿度选为10 g/m3[7]。

土壤由土和可调比例的水组成；其中土体积分数为50%，孔隙水的体积分数为0～50%，孔隙的体积分数为50%。参考Köhli[7]对土壤建模方式，土壤成分由SiO2和Al2O3(三氧化二铝)组成，体积比为75% SiO2和25% Al2O3，密度为2.86 g/cm3；土壤考虑到孔隙和孔隙水后，土壤容积含水量由0%增加到50%，土壤总密度从1.43 g/cm3可变化到1.93 g/cm3；本文的建模过程中不考虑土壤化学成分的差异，只考虑土壤体积含水量的差异。

3 模拟结果及分析

3.1 地表背散射中子能谱验证

地表背散射中子是指入射到土壤的宇宙射线中子与土壤相互作用后散射离开地表回到大气的宇宙射线中子；为与URANOS计算的地表背散射能谱归一化后进行比较，本文利用Geant4模拟土壤体积含水量为1%，空气湿度为0.33 g/m3时中子的输运过程，统计所有离开地表向上运动中子的动能得到了图5中红色的中子能谱曲线；计算结果表明，Geant4计算的地表背散射宇宙射线中子能谱与URANOS计算结果共同点是级联中子峰不明显，蒸发中子峰明显，不同点是Geant4计算的结果中热中子峰明显，慢化中子强度稍小。

图5 Geant4计算得到的地表背散射中子谱与URANOS计算的背散射宇宙射线中子谱归一化后的比较Fig.5 Comparison of normalized backscatter neutron spectra of cosmic rays calculated by Geant4 and URANOS

理论上，宇宙射线中子与土壤原子相互作用后，会有热中子产生。由于H对中子慢化做主要贡献，而Geant4计算的热中子峰明显，说明Geant4计算的结果更准确。Geant4计算中使用H的中子反应截面与URANOS使用的不同；URANOS软件中1H(氕)和2H(氘)的吸收截面和散射截面均采用的是2007年版的日本的评价数据库JENDL-4.0，而Geant4在中子能量大于20 MeV时采用的是日本的评价数据库JENDF，中子能量小于20 MeV时用的是美国的评价数据库ENDF/B-VII。而JENDL中子截面库在对高能区的中子输运计算比较准确，在对低能区的中子输运计算并不如ENDF/B-VII计算准确。而且在URANOS的计算中只考虑了中子的吸收截面和散射截面，在本模型计算中考虑了中子散射、吸收和裂变，使用了高精度HP物理模型。以上这些导致了基于Geant4的本模型计算的结果更准确，在0.025eV左右的热中子峰明显存在。故而，本模型计算结果比URANOS计算结果更准确。

3.2 地表附近宇宙射线中子能谱对土壤水分响应

统计进入到探测层的宇宙射线中子动能，得到如图6所示的土壤体积含水量不同时，距地表2 m处的宇宙射线中子能谱。图中纵坐标为计数N，通过除以入射宇宙射线中子数107，即可得到相对中子强度。计算结果表明，H对中子慢化起主要贡献；土壤水分体积含量增加时，慢化中子强度明显减少，热中子强度增加，蒸发中子峰降低，级联中子峰强度基本不变。出现该现象是由于土壤含水量的增加，被土壤元素慢化和吸收中子数增加，导致大气中蒸发中子峰、慢化中子峰和热中子峰的总强度要减少。H对中子起主要慢化作用，当H含量增加后，由于与热中子的散射截面大，吸收截面小，热中子强度相对慢化中子强度变大。

图6 地表以上2 m高处土壤水分含量不同时，本模型计算的宇宙射线中子强度能谱Fig.6 The differences of cosmic ray neutron spectra at 2 m height from the earth surface for the different soil moisture contents

图6中结果表明10-6～10-2MeV慢化中子对土壤水分含量有较好的分辨能力，慢化能区的中子强度随土壤含水量的增加而减少。同时，在热中子0.5×10-6MeV以下的中子对土壤水分在一定范围内变化时也具有分辨能力。目前的COSMOS系统[1,10]利用热中子探测器与快中子探测器的计数率比值来确定地表积雪或植被水当量，从而将测量的总水分含量减去地表积雪或植被水含量，更精确的得到土壤总含水量。

3.3 刻度函数的计算

为与Desilets[5]利用MCNP模拟计算的刻度函数关系式进行对比，统计能量在10-8～10-4MeV范围中子计数，该能量段包含了热中子和部分快中子。当土壤体积含水量为0%的干土壤时，探测层的计数为N0，其他土壤容积含水量下的，探测层的计数为N(假定探测层对所有能量的中子探测效率η=100%)。Geant4模拟计算得到的探测层计数如表2所示。

将容积含水量θvol和质量含水量θw与N/N0一一对应点绘制在图7中，并与Desilets[5]利用MCNP软件计算得到拟合形状函数(图中红色曲线)进行对比。采用与Desilets[5]模拟的N0方法刻度函数类似形式进行拟合，本模型计算结果的拟合函数如式(1)所示。

表2 土壤含水量不同时，探测层记录的10-8～10-4 MeV能区的中子计数Tab.2 Neutrons in the 10-8～10-4 MeV energy region recorded by the detection layer for the different soil moisture contents

(1)

当θ为θvol时，拟合参数a1=0.028 632 5，a2=0.567 39,a3=0.399 761；

当θ为θw时a1=0.021 725，a2=0.041 394,a3=0.378 991。

图7 土壤含水量与N/N0的关系Fig.7 N/N0 as a function of soil moisture

显然，采用Geant4计算得到的拟合函数与Desilets[5]采用MCNP计算得到的拟合函数比较接近，且对干燥土壤含水量变化更为灵敏。采用MCNP计算给出的是形状函数，同时适用于土壤容积含水量和质量含水量。Desilets在建模过程中用土壤中的弥散人造中子源近似实际入射到土壤的宇宙射线中子，这导致计算结果不准确，前言中已有说明，在此不再赘述。在实际应用过程中，通常会对Desiltes[5]给出的刻度函数进行修正，要么采取对N进行修正，要么对θ进行修正。

前文为了与Desilets计算结果对比，故而中子计数能量区间选择了10-8～10-4MeV。根据图6所示，显然能量区间10-6～10-2MeV中子计数与土壤含水量的响应关系更直观。不同能量区间的中子计数与土壤含水量的刻度函数拟合曲线如图8所示，两能量区间刻度曲线随着土壤含水量的增大差异逐渐增大。导致差异的主要原因是10-8～10-4MeV包含了热中子。然而热中子能量区间对土壤含水量变换并非线性响应，故能量区间10-6～10-2MeV的刻度函数结果具有更好的灵敏度。其拟合参数值为：

当θ为θvol时，a1=0.021 597;a2=0.172 635;a3=0.028 315；

当θ为θw时，a1=0.028 511;a2=0.399 549;a3=0.056 440。

图8 不同中子能区土壤含水量与N/N0的关系Fig.8 N/N0 as a function of the soil moisture for the neutrons within different energy ranges

4 结 论

本文针对常用的基于MCNP和URANOS的土壤含水量刻度函数计算模型中中子源设置不合理的问题，采用Geant4建立了新的刻度函数计算模型，中子源的入射高度和能谱是本模型不同于传统模型之处。在不同地表成分和地理位置以及海拔情形下，采用宇宙射线粒子输运程序PLANETOCOMICS分别计算了宇宙射线中子能谱并分析差异，发现地表对海拔600 m及以上的宇宙射线中子强度无影响，中子源设置在600 m海拔处，计算结果优于传统MCNP和URANOS模型的计算结果，且能量区间10-6～10-2MeV的中子对土壤水分的响应比10-8～10-4MeV能量区间更灵敏。基于以上结论给出了合适能量区间内的刻度函数拟合形式和拟合参数，结合图1流程计算N0可用于宇宙射线中子实时测量土壤含水量。