宇宙射线中子法测量区域土壤水分准确性季节变化研究

武 强 黄 娜 罗孳孳 唐余学 朱玉涵 徐倩倩

(1.重庆市气象科学研究所, 重庆 401147; 2.重庆市江津现代农业气象试验站, 重庆 402260;3.中国农业大学资源与环境学院, 北京 100193; 4.合肥市气象局, 合肥 230041)

0 引言

土壤水分通过径流、蒸散发、植被水分供给等方式，在地表与大气的物质和能量交换中起着重要作用[1-4]。土壤水分也作用于植物光合作用以及土壤微生物群落功能，从而影响着生态系统初级生产力[5-7]。因此，准确获得土壤水分含量及其时空动态变化，对水土保持、水资源管理及利用等研究具有现实意义[8-11]。

宇宙射线中子法是一种介于原位点测量与遥感大尺度反演的中等尺度非入侵、被动获取式土壤水分测量方法[12]，可以反映测量区域土壤水分平均状况，且不受土壤结构异质性以及土壤化学性质影响，不破坏土壤结构、可连续观测[13-14]。其原理为宇宙射线初级粒子经过一系列级联反应以及与大气粒子碰撞后，由激发态至稳定态的过程中会释放出快中子，快中子在近地面层与氢原子碰撞发生慢化直至吸收[15]，而土壤水分是陆表的主要氢源[16]，宇宙射线中子法通过测量与氢原子碰撞慢化的快中子强度来反演计算测量源区内土壤水分含量。

宇宙射线中子法区域土壤水分监测系统由美国研发生产，并在其境内组建形成70多个站点的大型观测网络[17]，围绕宇宙射线中子法的反演方法、参数订正、结果修正以及不同观测方法的比较验证等开展了诸多研究[12,15-17]。2012年以来，我国学者开始将宇宙射线中子法应用于不同下垫面的土壤水分测量，并就其测量效果开展了评估，在农田、草原、山地等应用场景下，均取得了较好的测量效果[18-21]，但受标准样本采样时间间隔较长或研究区域冬季冻土等影响，完整数据集时间跨度多为非冻土时段的数月时长，一方面忽略了短期天气条件变化的影响，另一方面未能涵盖生态环境条件周期性变化下的全部情况。较长时间序列的研究多关注于下垫面生态系统特征，如地表生物量对测量结果准确性的影响以及土壤水分反演模型适用性等方面，宇宙射线中子法土壤水分测量结果准确性在典型气象条件的影响下是否存在差异以及普适性尚未得到证实[22-23]。土壤水分时空变异性受到气候[24-25]、地形[26]、植被[27-28]等诸多因素的影响，由于对不同控制因子和过程的强烈非线性反馈，从而表现出土壤水分的时空变异性[29-30]。宇宙射线中子法作为新型的土壤水分测量方法，应用其进行土壤水分的长期连续测定，需对其应用场景的下垫面因素、气象条件等的动态变化影响进行全面考量，从而对应用可行性有全面的掌握。本文研究区域的下垫面植被以灌木型茶树为主，伴随茶叶适龄期的采摘、打顶修整等农事活动以及冬季休眠停止生长等因素，植被生物量变化较小，对宇宙射线中子法的土壤水分测量结果影响基本恒定，而伴随季节变化的气象条件是影响土壤水分变化最具变异性的因素。本文拟研究不同季节气象条件特征以及典型天气条件下宇宙射线中子法测量结果准确性，以季节为时间划分尺度，并选取典型气象条件时段，对比分析宇宙射线中子法在不同气象条件下的土壤水分连续测量效果差异与影响因素，以期为宇宙射线中子法在同一下垫面不同气象条件的适用性以及典型天气条件下的土壤水分数据质量控制与结果订正提供科学依据。

1 研究区域与研究方法

1.1 研究区域

试验区域设置于重庆市巴南区白象山，为典型山地地形下垫面，地势自西向东倾斜，整体海拔落差约100 m，坡度约18°，土质以壤土和黏质土为主，局部区域为砂质土或岩性紫色土。下垫面植被主要为灌木型茶树，株高50～80 cm，间种高大乔木作为防护林与遮阴树，空地多为原生草本地被植物。参数率定于2018年6月22日—8月11日进行，对比研究时段为2019年3月—2020年2月，按照春季(3—5月)、夏季(6—8月)、秋季(9—11月)、冬季(12—2月)的标准进行时间划分。

1.2 数据来源

研究区域中心(29°28′30″N，106°44′7″E，海拔558.60 m)为4 m×4 m观测场，观测场内安装有自动气象站，设置空气温度、空气相对湿度、气压、降水量、风向、风速、土壤温度等观测要素。观测数据用于研究时段内气象要素分析以及宇宙射线中子法参数率定。

土壤水分测定设置宇宙射线中子法(Cosmic-ray neutron sensing，CRNS)、频域反射法(Frequency domain reflectometry，FDR)、烘干称重法(Oven-drying method，ODM)3种不同时间、空间尺度的测量方法。参数率定试验采用烘干称重法多点采样获取区域土壤水分标准值，率定得到CRNS区域土壤水分反演模型。准确性季节变化特征研究由FDR得到多点连续观测数据计算结果作为区域土壤水分连续观测标准值。ODM与FDR均结合CRNS足迹，采取同样的垂直足迹线性加权、水平足迹等权重加权处理，得到对照数据集，与CRNS进行对比分析。

宇宙射线中子法区域土壤水分观测系统(Cosmic-ray soil moisture observing system，COSMOS)设置于研究区域中心观测场内，数据采样间隔设定为1 h，原始数据中子数经参数修正后，采用12 h滑动平均滤除随机误差，反演计算获得区域土壤含水率。

频域反射型土壤水分自动站两套，分别位于COSMOS所在的观测场内及COSMOS西南方向约100 m处(29°26′35″N，106°45′2″E，海拔580.70 m)，测量深度0～10 cm、10～20 cm、20～30 cm、30～40 cm、40～50 cm共5个梯度。数据采样间隔1 h，与COSMOS同步。数据集时间为2019年3月—2020年2月，期间每月1日、16日采用便携式土壤水分测定仪(FDR)按照图1烘干称重法取样位置进行补充测定，校准FDR区域土壤水分连续观测标准值。

图1 宇宙射线中子法参数率定样点分布示意图Fig.1 Distribution of CRNS parameter calibration sample points

烘干称重法为国际上沿用的土壤含水率测定标准方法，以多点、多土层土壤含水率计算得到区域土壤含水率作为标准值代表。土样于105℃下干燥至质量恒定，土壤失水质量与原土样质量的比值即为质量含水率，结合土壤容重换算得到体积含水率。

1.2.1测量足迹

CRNS水平测量足迹为以COSMOS为中心的面圆区域[16]，与气压呈反比关系，其关系式为

(1)

式中Ps,o——参考气压,取标准大气压101.325 kPa

Ps——实际气压，kPa

Rs,o——水平足迹参考半径，取300 m

Rs——实际气压Ps下的水平足迹半径，m

垂直测量足迹即探测深度，主要与测量区域土壤含水率有关，探测有效深度可表示为

(2)

式中Z——垂直足迹有效测量深度，cm

ρbd——测量源区土壤容重，g/cm3

ρw——液态水密度，g/cm3

τ——晶格水占矿质颗粒和束缚水质量总和的比例，%

θ——测量源区土壤体积含水率，cm3/cm3

CRNS水平足迹采用等权重加权方法进行处理，垂直足迹采用线性加权方法进行处理，从而获得测量足迹内的区域土壤含水率标准值。等权重加权方法即测量足迹内水平距离多点均匀采样，计算平均值。垂直足迹线性加权方法[31-32]即认为不同土层深度的土壤含水率对区域土壤含水率反演呈线性权重分配，线性加权函数为

(3)

(4)

可以得到

(5)

式中wt(h)——垂直测量足迹内的深度权重系数

h——土层深度，cm

1.2.2数据处理

CRNS是通过测量与氢原子碰撞而慢化了的快中子的强度来反演计算测量源区内土壤水分含量[15]。土壤水分以外的含氢物质如空气中的水汽也会起到慢化快中子的作用[33]；气压的变化改变高能粒子与空气粒子碰撞而丢失能量发生级联反应的几率[16]；太阳长周期变化如太阳黑子周期、昼夜变化等也会直接影响高能粒子密度[15]，从而影响土壤水分测量结果准确性。基于上述因素，需对CRNS快中子数原始数据进行相应的校正。校正函数为

fCWV=1+0.005 4(ρ-ρ0)

(6)

(7)

(8)

式中fCWV——水汽修正系数

fp——气压修正系数

fi——太阳活动修正系数

ρ——地表绝对水汽密度，g/cm3

ρ0——参考条件下地表绝对水汽密度，一般取0

P0——气压参考值，kPa

P——试验期内平均气压，kPa

L——高能中子质量衰减长度，按其随纬度变化规律，试验地点取137 g/cm3

Io——中子基准强度，使用海拔3 560 m瑞士少女峰站宇宙线数据作为中子基准强度

Im——探测器测得的中子强度

得到CRNS中子数校正函数为

NCORR=NRAWfCWVfpfi

(9)

式中NCORR——CRNS校正后中子数

NRAW——CRNS原始中子数

CRNS的反演函数[34]，即中子量与土壤含水率之间的函数关系为

(10)

式中θ(N)——测量源区平均土壤体积含水率

N0——测量源区土壤不含水条件下的中子数，其取值受地理纬度、海拔高度、土壤质地影响，对于同一测量源区，N0不变[35-37]

其中修正参数a0=0.080 8，a1=0.372，a2=0.115。

1.3 评价指标

CRNS区域土壤水分反演结果的评价采用1∶1线、偏差均值(Mean value of difference，d)、相关系数(Pearson correlation coefficient，R)、决定系数(Coefficient of determination,R2)、均方根误差(Root mean square error，RMSE)、纳什效率系数(Nash-Sutcliffe efficiency coefficient，NSE)、Kling-Gupta系数(Kling-Gupta efficiency coefficient，KGE)。其中d越小，一致性越高；RMSE表示测量值和标准值的偏离程度，RMSE越趋近于0，测量值与标准值越接近，误差越小；NSE越趋近于1，反演模型精度越高，反演结果与标准值一致性越好。KGE为综合考虑了模拟值与实际观测值的相关系数、变异性以及偏差量的一致性评价指标，KGE越趋近于1，表示模拟值与实际观测值综合一致性越好。

2 结果与分析

2.1 参数率定

2018年6月22日—8月11日开展参数率定试验，以ODM区域土壤含水率标准值作为基准进行CRNS参数率定。样点分布如图1所示，以COSMOS为圆心，垂直方向4个方位，50、100、200、250 m处共计16个取样点，取样深度0～10 cm、10～20 cm、20～30 cm和30～40 cm，结合试验期间CRNS垂直足迹实际情况，确定以0～20 cm土层测量结果，按照垂直足迹线性加权、水平足迹等权重加权方法进行处理，获得ODM区域土壤含水率标准值。试验期间每5 d取样1次，累计11次。计算CRNS、ODM获得的区域土壤含水率，结果比较如图2所示，计算得到RMSE为0.014 m3/m3， NSE为0.995，KGE为0.745，回归直线斜率与1∶1线偏差仅0.211，可见经参数率定后，CRNS对区域土壤水分的反演结果与ODM得到的标准值一致性良好，CRNS区域土壤水分计算模型可信度高。

图2 宇宙射线中子法参数率定结果Fig.2 Comparison of CRNS parameter calibration results

2.2 测量足迹季节变化

2019年3月—2020年2月CRNS测量足迹如图3所示，试验期内水平足迹较为稳定，年平均值为320.5 m，最小值314.5 m，最大值326.0 m，变化率为±2%；其中春季平均值321.2 m，夏季平均值323.5 m，秋季平均值319.5 m，冬季平均值318.4 m。垂直足迹变化较大，年平均值16.9 cm，最小值12.3 cm，最大值25.3 cm，变化率为-37%～50%，尤其夏季、秋季受降水及持续高温天气等影响，土壤持续失墒或降水补充等导致垂直测量足迹出现相应变化，具体表现为与土壤含水率呈负相关性，即土壤含水率越低，垂直足迹越大，土壤含水率越高，垂直足迹越小；其中春季平均值19.0 cm，夏季平均值17.1 cm，秋季平均值15.7 cm，冬季平均值16.0 cm。可见在本试验区域，CRNS水平足迹基本稳定，对时间、气象条件以及土壤含水率等的变化响应不明显；垂直足迹受到土壤含水率变化的影响而波动，但是基本维持在15.7～19.0 cm之间，受降水影响出现过极小值12.3 cm，但在雨后随着地表水的径流与渗透，可恢复到平均水平。耕作层有效深度大约20 cm，是土壤水分和养分含量发生波动的主要土层。因此，在山地地形下垫面，CRNS水平足迹基本稳定，随季节变化不明显，垂直足迹在不同季节随降水差异而有一定波动，但均集中在耕作层范围内变化，即CRNS能够较为稳定地反映测量源区内耕作层土壤水分的变化情况。

图3 宇宙射线中子法测量足迹变化Fig.3 Diagram of CRNS footprint

2.3 测量结果准确性季节变化

统计研究区域气象数据得到，春季平均气温17.0℃，极端最高气温34.3℃，极端最低气温5.7℃。雨日数53 d，降水总量431.5 mm，雨量分布情况为3月20%，4月30%，5月50%。春季CRNS与FDR的测量结果对比如图4a所示，总体变化趋势基本一致，两种方法测量偏差均值为0.019 m3/m3。对两种方法得到的区域土壤含水率日值进行拟合(图4a)，得到线性方程y=1.638x-0.177，R2=0.442。具体地，3月气温较为稳定，降水较少，偏差均值为0.004 m3/m3；4月气温回升，土壤水分蒸腾量增加，且在4月中下旬有两次较强降水过程，土壤水分含量变化较大，偏差均值为0.017 m3/m3；5月雨日数较多，雨量较大，土壤增墒明显，偏差均值为0.041 m3/m3。

图4 宇宙射线中子法区域土壤含水率季节变化Fig.4 Seasonal variations of regional soil moisture by CRNS

夏季平均气温25.6℃，极端最高气温37.7℃，极端最低气温17.6℃。雨日数43 d，降水总量431.3 mm，雨量分布情况为6月42%，7月53%，8月5%。夏季CRNS与FDR的测量结果对比如图4b所示，两种方法的变化趋势一致性较好，两种方法测量偏差均值为0.020 m3/m3。对两种方法得到的区域土壤含水率日值进行拟合(图4b)，得到线性方程y=1.299x-0.083，R2=0.777。在降水集中的6、7月，受降水影响气温波动大，土壤短期失墒与水分补充变化较大，偏差均值分别为0.039、0.023 m3/m3，雨量较少的8月，气温较高，平均值26.3℃，土壤持续失墒，8月中下旬，日最高气温均在35℃以上，8月下旬，土壤含水率状况逐渐由轻旱发展为中度干旱，整个8月土壤含水率测量偏差均值仅0.001 m3/m3。

秋季平均气温17.3℃，极端最高气温35.4℃，极端最低气温5.9℃。雨日数54 d，降水总量367.9mm，雨量分布情况为9月25%，10月57%，11月18%。秋季CRNS与FDR的测量结果对比如图4c所示，两种方法测量偏差均值为0.034 m3/m3，但是在几次主要降水过程时，表现出较为明显的测量差异。对两种方法得到的区域土壤含水率日值进行拟合(图4c)，得到线性方程y=1.407x-0.103，R2=0.648。9月雨日数较少，降雨过程较为集中，偏差均值为0.024 m3/m3；10月雨日数多，但多为较小量级的降水，偏差均值为0.062 m3/m3；11月雨日数较少，雨量较小，偏差均值为0.034 m3/m3。

冬季平均气温8.2℃，极端最高气温16.8℃，极端最低气温0.3℃。雨日数43 d，降水总量86.4 mm，雨量分布情况为12月30%，1月34%，2月36%。冬季CRNS与FDR的测量结果对比如图4d所示，可见两种方法的变化趋势一致性较好，两种方法测量偏差均值为0.001 m3/m3。对两种方法得到的区域土壤含水率日值进行拟合(图4d)，得到线性方程y=1.929x-0.349，R2=0.182。冬季降水过程分布较为平均，且均为小雨量级降水，12、1、2月偏差均值分别为0.018、0.011、0.021 m3/m3。

各季节CRNS与FDR统计指标如表1所示，二者平均值μ差异较小，但是CRNS标准差σ明显偏大，主要与二者测量原理差异有关，CRNS通过测定与氢原子碰撞慢化的中子数反演计算得到土壤含水率，快中子计数在一定范围内变化，较为离散；而FDR传感器设置于土层中，通过电磁波在不同土壤含水率即不同介电常数介质中的变化来反演土壤含水率，除蒸腾与降水渗透以外，不受其他地上因素干扰，而土层中的水分变化相对较为稳定连续，因而测得土壤含水率也较为稳定。相关系数R夏季、秋季均达到0.8以上，冬季最低，为0.427。RMSE春季、夏季均为0.035 m3/m3，秋季0.055 m3/m3，冬季0.031 m3/m3。NSE与KGE均在夏季最高，分别达到0.604、0.508，其余季节均较低。主要与各季节的降水、土壤水分蒸散以及下垫面特征有关。结合研究区域下垫面蒸散量计算得到春季有效降水量为40.8 mm，夏季200.0 mm，秋季136.4 mm，冬季13.7 mm。由于研究区域下垫面主要为茶树，种植较为密集，树冠对于降水截留作用明显，较小量级的降水受到茶树冠层的截留削减，到达土层的入渗水量过小不足以引起频域反射传感器的明显响应，而CRNS反演结果涵盖了测量区域内所有含氢物质的总和，包括地表植被，尤其是降水期间冠层的截留水量。茶树下垫面主要为密植灌木，生长较为缓慢，加之采茶以及打顶等农管措施对幼嫩芽叶的削减，地表生物量总体较为稳定，因而主要是降水引起的冠层截留水量会导致测量结果的偏差。尤其在冬季，天气条件以阴天为主，总日照时数90.3 h，日照时数在1 h以上日数仅22 d，期间有效降水量极小，即水分蒸腾与补充基本接近平衡，埋入式的FDR传感器得到的土壤含水率在整个冬季也十分稳定，维持在0.36～0.39 m3/m3，而历次降水过程日降水量均不超过10 mm，主要为冠层截留水与土表吸收，该部分氢源的增加导致CRNS明显的响应，尤其2月有几次较为集中的降水过程，在FDR土壤含水率基本不变的情况下，CRNS表现出较为明显的结果差异。

表1 宇宙射线中子法与频域反射法区域土壤水分测量结果统计指标Tab.1 Statistical index values of regional soil moisture measurement results by CRNS and FDR

2.4 典型土壤水分条件下宇宙射线中子法测量准确性变化

进一步分析典型土壤水分条件下CRNS测量准确性变化，选取2019年3月—2020年2月自然条件下多雨偏湿以及高温干旱时段作比较。7月作为多雨偏湿的代表时段，平均气温25.7℃，极端最高气温37.5℃，极端最低气温19.1℃。雨日数18 d，降水总量228.4 mm，其中7月8日、7月19日、7月29日达到大雨量级，7月31日达大暴雨量级，其余降水时段为小雨量级，在整个多雨偏湿期，土壤表现为逐渐增墒效应，土壤含水率约增加0.06 m3/m3。对比多雨偏湿期CRNS与FDR的测量结果，两种方法测量偏差均值为0.002 m3/m3。对几次强降水过程，CRNS与FDR区域土壤含水率测量值均表现出明显的瞬时响应。对两种方法得到的区域土壤含水率日数据进行拟合(图5a)，得到线性方程y=1.366x-0.110，R2=0.272。8月作为高温干旱的代表时段，平均气温28.2℃，极端最高气温37.7℃，极端最低气温18.1℃。雨日数9 d，降水总量20.1 mm，其中除8月31日达到中雨量级，日降水量为12.2 mm，其余降水时段以零星小雨为主，日降水量不超过2 mm，在整个干旱期，土壤持续失墒，土壤含水率降低约0.13 m3/m3。对比高温干旱期CRNS与FDR的测量结果，两种方法测量偏差均值为0.001 m3/m3。对两种方法得到的区域土壤含水率日数据进行拟合(图5b)，得到线性方程y=1.031x-0.011，R2=0.926。

图5 典型土壤水分条件下宇宙射线中子法测量结果Fig.5 Result by CRNS in typical soil moisture conditions

典型土壤含水率条件下CRNS与FDR统计指标值如表1所示，二者平均值μ在高温干旱期基本一致，多雨偏湿期较大。在多雨偏湿期、高温干旱期，R分别为0.522、0.962；RMSE分别为0.034、0.014 m3/m3。NSE与KGE在多雨湿润期均为负值，在高温干旱期则高达0.925、0.919。可见在降水干扰较小，即下垫面氢源较为稳定的时段，CRNS与FDR的测量结果可表现出良好的一致性，测量结果偏差极小，CRNS区域土壤含水率反演模型准确度高，反演值与标准值匹配度良好。在降水较多的土壤湿润期，即便受降水过程干扰，CRNS反演区域土壤含水率结果准确性受到一定影响，但是在土壤含水率的变化趋势以及量级较大降水过程时土壤含水率的响应等基本一致，测量结果偏差均值仅为0.023 m3/m3，测量结果准确性较好。

3 讨论

不同土层往往表现不同的水文过程和生态功能[38]，浅层土壤受降水入渗和蒸散发影响大，同时也是植被生长常用水分来源，深层土壤水分相对较为稳定，形成“土壤水库”。山地地形下垫面条件下，CRNS的水平测量足迹基本稳定，不随时间、气象条件以及土壤含水率等的变化而发生明显变化，因而在代表性样点的选取时，应注意如岩石、建筑、水源等水平足迹范围内干扰项的规避；在测量结果的订正，如水平足迹范围内植被含水率的排除计算等方面，则可根据水平足迹较为稳定的特点而给出理想试验方案。垂直足迹则主要受到土壤含水率变化的影响而波动，但测量足迹基本在耕作层范围内变化，因此，CRNS可应用于耕作层土壤水分反演，或地表生物量、植物水分截留等相关研究。

比较CRNS在不同季节的土壤水分测量结果，总体上CRNS与FDR的测量结果一致性较好，在土壤含水率较为稳定(冬季)表现为与FDR极小的偏差，在少雨持续失墒(夏季)以及多雨水分得到补充(秋季)的土壤含水率剧烈变化的时段，CRNS与FDR表现出一致的变化趋势。但是由于研究区域下垫面主要为密植的灌木型茶树，有明显的降水截留效果，被截留降水无法参与地面径流或地下水分运动过程[39]，但是截留水作为氢源可以被CRNS捕捉；通常雨量较大时，植被截留率较低，但是雨量较小时，植被截留则显得尤为重要，研究显示，在湿润森林区，植被截留损失占年降水量10%～30%[40-42]，而在降水偏少的地区，植被截留损失可以达到40%～50%[43]。因而在试验过程中，存在FDR几乎无变化，但CRNS已出现明显波动的情况。对于较小量级的短期降水(春季)与量级较大的持续降水(夏季)，CRNS较FDR能够更加灵敏地响应，就土壤水分测定而言，结果包含了冠层截留水等干扰，但对于水资源评估、作物水分平衡等研究具有重要的指示意义。

进一步分析典型土壤水分条件下的CRNS测量结果准确性变化，在高温干旱期，CRNS表现出最佳测量结果，与FDR表现出极高的一致性。在多雨偏湿的增墒时段，则受降水过程影响，相较于高温干旱时段表现为较大偏差。因此在应用CRNS对测量区域的土壤水分整体评价中，合理排除降水时段的土壤水分测量偏差尤为重要。

4 结束语

CRNS水平测量足迹的季节变化较为稳定，垂直足迹在耕作层范围内有一定波动。不同季节降水过程的量级差异是引起CRNS测量结果偏差的主要因素，CRNS对不同量级降水均可作出明显响应。夏秋两季，在高温少雨导致的持续失墒以及较大量级降水过程等土壤水分变化较为剧烈的时段，CRNS与FDR测量结果一致性较好。春季多较小量级降水，作物冠层截留作用导致小量级降水不足以引起探入式FDR的响应，但CRNS则可以灵敏捕捉到水分的变化，二者一致性略差。典型土壤水分条件下的CRNS测量准确性变化，进一步说明降水过程对CRNS准确性的干扰，而在持续高温干旱时段则一致性良好，RMSE、NSE与KGE分别为0.014 m3/m3、0.925、0.919。