基于平均影响值算法和BP神经网络的根区土壤湿度估算

袁玲， 方秀琴,2*， 郭晓萌， 杨露露， 张晓祥， 任立良

(1.河海大学水文与水资源学院， 南京 210098； 2.河海大学海岸灾害及防护教育部重点实验室， 南京 210024)

土壤湿度(soil moisture, SM)是陆地表面系统的关键水文和生态变量，而根区土壤湿度(root zone soil moisture, RZSM)不仅影响降雨在入渗水和地表径流中的分配，是洪水灾害发生的重要孕灾环境，还影响着地表植被的水分吸收和蒸腾作用，进而影响植被的结构和组织[1]。此外，由于土壤含水量和土抗剪强度之间的内在关系，RZSM对土壤侵蚀的加速和产沙量的增加也有一定影响。因此，对RZSM的预测不仅可以为洪水风险分析过程提供依据，也是农业规划和干旱监测工作中的重要环节[2]。

在过去的几十年里，为了获取准确的土壤水分信息，发展了多种技术和仪器。现场土壤传感仪，如中子探头、电磁传感器和脉冲传感器被广泛应用于获取监测站点的土壤湿度数据[3]，但是这些监测方法需要耗费大量人力物力，而且局限于站点，得到的土壤湿度资料存在空间上密度低、时间上不连续，进而存在大量缺测值的问题[4]。此外，由于土壤质地、植被、地形等地表特征具有较高的空间异质性，因而利用点或场尺度的观测数据来反映土壤水分的空间分布具有相当大的局限性[5]。遥感技术因其能够获取大尺度连续区域的土壤水分资料而得到了广泛的关注，取得了重要进展，杨丹阳等[6]基于多时相高分一号影像建立反演模型对实验区土壤湿度进行了反演，达到了较高的精度，然而利用遥感技术获得土壤水分信息限制在地表5 cm以内，无法获取到根区的土壤湿度值。总体而言，现有的观测技术不足以满足RZSM获取的需求，在许多未测量区域和观测站密度较低的区域，RZSM观测在时间和空间上都存在大量的空白，而遥感技术由于测量深度有限，无法填补根区土壤湿度的空白。

近年来，研究人员提出了多种利用地表信息反演根区土壤湿度的技术。在众多方法中，将土壤-植被-大气传输模型和地表模型与数据同化技术相结合被证明是一种有效的RZSM估算方法，然而基于物理过程的模型的参数优化需要大量观测数据来确定，并且模型实现过程复杂，其中难免存在重要物理过程概化的缺陷。针对此问题，一些研究提出了基于数据驱动的RZSM估算技术。与传统水文模型不同，基于数据驱动的估算方法不需要通过明确的模型结构来表示所涉及的水文过程，仅需保证与目标值相关的训练数据的真实性和可靠性，刘立文等[7]研究发现土壤湿度通过各种水文过程与气象因素存在高度相关性，因此可以通过地表气象要素数据驱动来反演根区土壤湿度。

在各种类型的数据驱动方法中，人工神经网络被广泛应用于土壤湿度的估算。Lee等[8]利用研究区610 km2区域内的10个监测站点数据来训练人工神经网络，并对该区的一个独立站点的SM进行预测，获得了良好的估算结果；杨娜等[9]以甘肃省为例，对气象要素与土壤湿度之间的作用关系进行了研究，并建立BP神经网络完成了土壤湿度值的估算实验，获得了较好的预测结果；Meng等[10]提出了一种基于GF-3卫星数据和人工神经网络的农业地区土壤水分反演方法，使用从高级积分方程模型生成的训练样本数据集上对神经网络进行训练和测试，得到了较好的估算性能。虽然这些研究都显示了人工神经网络方法在土壤水分估计中的潜力，但都局限于小范围区域，鲜有关于模型在不同气候和土壤条件区域的适用性评价。在前人的研究基础上，现选择中国境内4个典型土壤湿度长期监测台站为研究对象，以地表信息为输入，根区土壤湿度值为输出构建BP神经网络模型，分别对不同站点的RZSM进行预测，并基于平均影响值(mean impact value，MIV)算法计算输入地表要素的重要性，探讨不同气候和土壤区域的根区土壤湿度的主要驱动因素和模型性能表现差异，对BP神经网络利用地表特征信息计算根区土壤湿度的性能和区域适用性进行研究。

1 实验区与数据

1.1 实验区概况

选择分别位于中国境内东北黑土区、中南红壤区、西南紫色土区和西北半干旱土区的四大整合农田土壤湿度监测场地为实验区域，分别为辽宁沈阳(123.37°E, 41.52°N)、湖南桃源(111.47°E, 28.92°N)、四川盐亭(105.46°E, 31.27°N)和新疆阿拉尔(80.85°E, 40.62°N)，站点土地利用类型和气象条件情况如表1所示。

1.2 数据资料

1.2.1 土壤湿度数据

土壤湿度数据来自中国土壤数据库(http://vdb3.soil.csdb.cn/)提供的站点土壤水分动态数据，该数据中土壤湿度用体积含水率来表示，即土壤中水分占有的体积和土壤总体积的比值。土壤湿度获取方法为中子管法，数据时间步长为5 d，从每月1日开始，5 d一采样，对每个监测站点在不同采样地得到的土壤湿度取平均值作为该监测站点当日的土壤湿度值。5个站点提供10～150 cm每隔10 cm深度的土壤湿度监测值，综合考虑数据的完整性，最终共选择7层有效深度的土壤湿度数据为研究对象，对应深度和时间序列如下：盐亭站第1层至第7层垂直深度分别为10、20、30、40、50、60、100 cm，时间范围为2007年1月5日—2010年12月31日；沈阳站为10、20、30、40、50、70、90 cm，时间范围为2004年4月30日—2011年10月4日；桃源站为10、20、30、40、50、70、90 cm，时间范围为2004年1月11日—2011年12月20日；阿拉尔站为10、20、30、40、50、70、90 cm，时间范围为2008年1月5日—2011年11月30日。其中，除沈阳站每年11月—次年3月存在数据缺测外，其余站点均为连续完整数据。将10 cm土壤湿度视为表层土壤湿度数据，其余层的土壤湿度视为对应不同深度的RZSM。

表1 实验站点概况Table 1 Overview of experimental sites

对于SM时间序列中的少量缺失数据依据K邻近插值算法(K-nearest neighbors，KNN)的基本思想进行插补，即在K个最近邻对象中通过设定函数值来填充缺失值。研究设定K为10，使用10个最邻近对象的加权平均值来填补缺失数据，得到4个野外台站共1 310组完整时间序列的土壤湿度观测值，其中盐亭站287组，沈阳站306组，桃源站414组，阿拉尔站303组, 将每个站点的数据按照6∶2∶2的比例分为训练集、验证集和测试集，其中训练集用于模型构建，验证集用于调整模型的超参数，测试集用来评估模型最终的泛化能力。

1.2.2 气象数据

气象要素是土壤湿度计算中重要的地面特征变量，能在一定程度上反映土壤水分的变化情况，从中国气象数据网(http://data.cma.cn)提供的观测站点中提取与4个土壤湿度观测站地理位置一致的站点，获取对应站点的地面气候资料日值数据集，并将其整理为模型所需的地表气象输入数据，包括日平均气温(daily mean temperature，TEM)、日累计降水量(daily cumulative precipitation，PRE)、日照时长(sunshine duration，SSD)、蒸发量(evaporation，EVP)、相对湿度(relative humidity，RHU)、风速(wind speed，WIN)6种气象数据，并提取与土壤湿度数据的时间序列对应的各气象因子的日值。

1.2.3 植被指数

近年来，归一化植被指数(normalized difference vegetation index, NDVI)和归一化差异红外指数(normalized difference infrared index, NDII)被逐渐用于对土壤湿度的研究。Schnur等[11]利用美国一些土壤湿度站点实测数据和MODIS(moderate-resolution imaging spectroradiometer)的NDVI和EVI产品，研究发现植被指数与每个深度的土壤湿度都高度相关，当植被指数滞后5～10 d时，相关性达到最高，且将回归模型应用到相似的远处站点时，估算值能够解释观测值的50%～88%；Sriwongsitanon等[12]发现在流域尺度上NDII与根区土壤湿度显著相关，两者的R2在旱季可以达到0.87；Joiner等[13]发现NDVI和NDII在几天到几周的时间内与根区土壤湿度有显著的相关性；Castelli等[14]用Landsat提取的NDII序列与突尼斯干旱区的多个土壤湿度观测站点数据进行分析，发现所有站点上两者都有显著相关(考虑整年数据的R2在0.62～0.67)，Castelli等[15]也将NDII成功用于非洲其他地区的根区土壤湿度估算研究；Pradhan[16]在美国西北部山区的研究证明NDVI能够反映干旱/半干旱区的根区土壤水分的动态变化。Liu等[17]在中国华北小麦地的实验结果表明太阳诱导叶绿素荧光(solar-induced chlorophyll fluorescence，SIF)和修正的NDVI都与根区土壤水分有显著正相关。这些研究都证实植被的光谱指数(NDVI、NDII、SIF等)可能成为大尺度根区土壤湿度监测的有效途径。

基于此，利用MODIS地表反射率产品MOD09GA，根据红波段(620～670 nm)、近红外波段(841～876 nm)和中红外波段(1 628～1 652 nm)的平均反射率ρ，按式(1)和式(2)计算4个研究站点对应时间序列的NDVI和NDII值，即

NDVI=(ρ2-ρ1)/(ρ2+ρ1)

(1)

NDII=(ρ2-ρ6)/(ρ2+ρ6)

(2)

式中：ρ1、ρ2、ρ6分别为红波段、近红外波段和中红外波段的平均反射率。

2 研究方法

2.1 输入变量选择

与根区土壤湿度相关的表层土壤湿度、气象驱动数据和植被指数共9个模型作为输入变量，为了使模型更加符合真实情况，考虑变量存在的时间滞后问题，将变量分为三组，其中，表层土壤湿度(SSM)、蒸发量(EVP)、相对湿度(RHU)、风速(WIN)、日照时常(SSD)为当日记录值，反映日增量；日累计降水量(PRE)和日平均气温(TEM)为前a天的累计值，为模型提供对过去值的记忆；归一化植被指数(NDVI)和植被水分指数(NDII)为滞后b天数据，这是由于植被生长状况对土壤含水量的响应需要一定时间，且Pradhan[16]研究表明，NDVI和NDII与土壤水分在5～10 d时间滞后时表现出最为显著的相关性。a和b的具体值通过比较每个站点变量与土壤湿度之间皮尔森相关系数值的大小及显著性来确定，最终确定4个站点的输入变量为：表层土壤湿度(SSM)、蒸发量(EVP)、相对湿度(RHU)、风速(WIN)、日照时常(SSD)、前45 d降水量累计值(PRE)、前45 d气温累计值(TEM)、滞后7 d NDVI值(NDVI)、滞后7 d NDII值(NDII)。

2.2 BP神经网络

近20年来，人工神经网络在水文学研究领域得到了普及与应用[18-19]，在多类神经网络中，BP神经网络因其较强的非线性映射能力而被广泛使用。

2.2.1 模型结构和参数

采用三层BP神经网络，包含输入层、隐含层和输出层。输入层神经元个数为9，对应模型的9个输入变量；输出层神经元个数为1，对应站点当前垂直深度的根区土壤湿度值；隐含层神经元个数依据经验公式[式(3)]并结合试错法确定，公式为

(3)

式(3)中：M、P、N分别为输入层、隐含层和输出层神经元的数目；L∈[1，10]。

隐含层神经元个数的具体确定方法为：以每个站点20 cm深度土壤湿度值的估算为例，将隐含层神经元个数从4增加到14，计算测试集中预测值与目标值之间的均方根误差(root mean squared error，RMSE)，结果如图1所示。从图1可以看出，RMSE并没有大幅度上升或下降的趋势，但总体来说，当隐含层神经元个数为6时，4组模型的综合性能达到最优，故最终确定隐含层神经元个数为6。

图1 均方根误差RMSE与隐含层神经元个数的关系Fig.1 The relationship between the RMSE and the number of neurons in the hidden layer

以测试集预测精度为依据，设置网络其他参数如下：学习率依据站点不同略有差异，盐亭站为0.1、沈阳站为0.25、桃源站为0.45、阿拉尔站为0.15；最大迭代次数均为1 000次；激活函数为sigmoid函数；优化算法为梯度下降法。

2.2.2 模型性能评价指标

使用以下3个指标来反映估计值与目标值之间的拟合优度。

(1)决定系数R2。R2为相关系数的平方，反映输入变量对结果的解释程度，其值越大表示模型拟合优度越好。

(4)

(2)均方根误差RMSE。RMSE为预测值与真值偏差的平方和与观测次数比值的平方根，反映预测的精密度和模型的稳定性，其值越小表示模型仿真性能越好。

(5)

(3)平均绝对误差MAE。MAE为所有预测值与真实值之间绝对误差的平均值，反映预测值实际误差的大小，其值越小表示模拟结果越接近真实值。

(6)

2.3 输入变量重要性计算

为了计算9个输入变量对4个站点不同深度土壤湿度预测结果的影响，更加直观地反映地表各特征在根区土壤湿度预测中的重要程度，使用平均影响值(mean impact value，MIV)算法分别计算4个站点各输入变量的影响值来反映其重要性。

平均影响值算法是Dombi等[20]提出的一种反映神经网络中权重矩阵变化情况的方法，被广泛应用于非线性算法中的变量排序，徐磊等[21]将BP神经网络和MIV算法相结合，确定了公交场站污染物暴露水平的主要影响因素。

假设神经网络共有l个输入变量， MIV算法步骤如下。

步骤1依次将其中一个自变量Xi(i=1,2,…,l)的数值增加10%，其他l-1个变量不变，得到新的样本后将其输入训练好的网络进行仿真得到l组模拟结果Pi(i=1，2，…,l)。

步骤2依次将其中一个自变量Xi(i=1，2，…,l)的数值减小10%，其他l-1个变量不变，得到新的样本后将其输入训练好的网络进行仿真得到l组模拟结果Qi(i=1，2，…,l)。

步骤3增加自变量得到的模拟结果减去减小自变量的模拟结果求平均，得到变量Xi的MIV值。

3 结果与分析

3.1 模型性能

使用验证集来评估模型性能，为了避免偶然误差，提升模型的可信程度，对4个站点每个深度的模型各训练500次，取500次R2、RMSE、MAE的平均值来衡量模型性能，结果如表2～表5所示。从站点来看，4个站点中桃源站的整体模型性能最好，而从土层深度来看，浅层土壤湿度的模拟精度高于要深层土壤。对4个站点的评估结果求平均(图2)，结果如表6所示，结合表6和图2可以看出，随着土层深度的增加，模型的决定系数R2不断减小，而均方根误差(RMSE)和平均绝对误差(MAE)则呈现明显的增加趋势。这与Pan等[5]的研究结果一致，即仅通过地表资料足以反演根区土壤水分状况，但反演精度会随着土壤深度的增加而减弱。

表2 盐亭站不同土层深度的模型性能Table 2 Model performance in different depths at Yanting

表3 沈阳站不同土层深度的模型性能Table 3 Model performance in different depths at Shenyang

表4 桃源站不同土层深度的模型性能Table 4 Model performance in different depths at Taoyuan

表5 阿拉尔站不同土层深度的模型性能Table 5 Model performance in different depths at Alaer

表6 4个站点不同土层深度下模型评价指标均值Table 6 Averages of model evaluation indexes under different depths at four stations

图2 4个站点不同土层深度下模型评价指标均值Fig.2 Averages of model evaluation indexes under different depths at four station

3.2 地表特征信息的重要性

使用平均影响值(MIV)算法分别计算4个站点各地表特征变量对RZSM计算模型的贡献，结果如表7～表10所示，其中红色背景代表变量对RZSM为正影响，且颜色越深影响值越大；绿色背景代表变量对RZSM为负影响，且颜色越深影响越大。

(1)表层土壤湿度(SSM)。在众多地表特征变量中，无论对于哪个站点，SSM对于根区土壤水分计算模型来说都是最重要的输入变量，这是因为，表层土壤水通过下渗进入深层土壤，是根区土壤水分的直接来源，也是RZSM计算的重要依据。

(2)降水(PRE)、归一化植被指数(NDVI)和归一化差异红外指数(NDII)。这3个特征变量的影响值在4个站点多为正值，表明这3个地表特征在不同气候和土壤区对不同深度的RZSM具有正面影响。这是因为，在自然状态下，降水是土壤水分最重要的来源之一，降水量的多少直接影响土壤中水分含量，而地表植被对土壤湿度既存在正面影响也存在负面影响，当植被根系层对土壤水的保蓄作用大于冠层截留时，NDVI和NDII与土壤含水量之间表现为正相关。但从桃源和阿拉尔结果显示，对于土地利用形式为水田和水浇地的两个站点来说，PRE、NDVI和NDII指标的影响值都很小，这可能是由于外界因素的干扰，改变了土壤水分的主要来源途径，使得降雨和植被的变化不能很好地反映土壤含水量的变化。

(3)蒸发量(EVP)。在沈阳和阿拉尔两个位于温带地区的站点，EVP对RZSM表现为负影响，而在位于亚热带的两个站点，EVP与浅层(20 cm)土壤水分存在较为明显的负相关关系，而在深层(30 cm以下)对RZSM的影响不明显。这可能是因为亚热带气候雨热同期的特点，当夏季地面蒸发量大的时候雨水也充沛，导致EVP这一地面特征并不能很好地反映深层根区土壤湿度的变化情况。

(4)日照(SSD)。在盐亭、沈阳和桃源3个位于湿润地区的站点，SSD对RZSM没有表现出明显的正影响或负影响，而在处于干旱区的阿拉尔站，SSD在浅层(50 cm及以上)对RZSM表现为明显的负影响，这是因为在湿润地区，土壤水分不仅具有多重来源，而且水分含量也受蒸发量、植被冠层截留等多重因素的影响，不确定性大，而在位于干旱区的阿拉尔站，土壤水分的主要来源不是自然降水，而是人工灌溉，日照时长会影响土壤水分的蒸发速度，对干旱区土壤湿度存在明显的负影响。

(5)相对湿度(RHU)、风速(WIN)和气温(TEMP)。这3个特征变量对RZSM的影响在本研究中的4个实验站点没有表现出明显的规律。

总的来说，根区土壤湿度值计算的主要依据是地表土壤湿度，多个地面特征之间的相互作用为辅助条件，且不同特征在不同气候区和土壤类型区对RZSM的影响情况也不一致。

表7 盐亭站不同土层深度变量的MIV值Table 7 MIV values of variables in different depths at Yanting

表8 沈阳站不同土层深度变量的MIV值Table 8 MIV values of variables in different depths at Shenyang

表9 桃源站不同土层深度变量的MIV值Table 9 MIV values of variables in different depths at Taoyuan

表10 阿拉尔站不同土层深度变量的MIV值Table 10 MIV values of variables in different depths at Alaer

3.3 模型应用

为了验证模型的有效性，将训练好的神经网络应用到4个站点，可用时间序列的全部数据来预测不同深度的根区土壤湿度，并将其与观测值进行比较，结果如图3所示。

从整体上看，几乎所有预测结果都能较好地反映根区土壤湿度值的变化趋势，具体到预测精度上来说，桃源和阿拉尔两个站点的预测结果与真实值最为接近，沈阳站模拟结果较好，而盐亭站预测值与真实值之间的偏差较大。

同一站点不同深度的土壤湿度值略有差异但变化趋势基本一致，说明以表层土壤水分数据为依据，结合气象要素预测深层土壤湿度值这一方法是可行的。

4个站点的预测结果普遍存在高值偏高和低值偏低的问题，即模型对极值的预测还不够准确。

模型对于浅层土壤湿度的估算性能普遍比深层土壤湿度要好，50 cm以上土层的拟合精度高于50 cm以下土层。

总的来说，构建BP神经网络模型，依据表层土壤湿度和气象要素来估算根区土壤湿度，在中国不同气候区的4个站点都取得了较为理想的精度，证实了这一方法的可行性。

由于存在数据缺测，沈阳站缺失每年11月—次年3月数据值图3 4个站点不同土层深度的模型预测值与真实值比较Fig.3 Comparison of predicted value with true value in different depths at the four sites

4 结论

选择中国境内不同气候区的4个土壤湿度监测站点为研究对象，构建了以地表信息为输入，根区土壤湿度值为输出的BP神经网络模型，对不同站点的RZSM进行了预测，并基于MIV算法计算了输入变量的重要性，对BP神经网络的性能和区域适用性进行了探讨。结论如下。

(1)在选择合理参数的前提下，利用BP神经网络模型来预测根区土壤湿度能到达较好的精度，在20～90 cm和20～100 cm深度，模型在验证集上的平均R2值在4个站点分别为0.79、0.69、0.66、0.56、0.51和0.47；RMSE为1.91%、2.17%、2.51%、2.71%、2.82%和3.08%；MAE为1.44%、1.61%、1.75%、1.89%、2.04%和2.35%。说明构建的BP神经网络模型能应用于中国不同气候和土壤区的四大整合农田土壤湿度监测场根区土壤湿度的估算，且越接近表层，土壤水分估算模型的性能越好。

(2)从地表特征的重要性计算结果来看，SSM是模型最重要的输入变量，此外，对位于不同气候区的站点来说，气象输入要素对根区土壤湿度的影响都不相同。其中，PRE、NDVI和NDII 3个特征的影响值在4个站点多为正值；EVP在两个温带地区的站点对RZSM表现为负影响，在位于亚热带的两个站点，对浅层(20 cm)SM表现为负影响，而对深层(30 cm以下)SM的影响不明显；SSD在3个位于湿润地区的站点，没有对RZSM表现出明显的正影响或负影响，而在处于干旱区的站点，SSD对浅层(50 cm及以上)SM表现为明显的负影响；RHU、WIN和TEMP对4个实验站点RZSM的影响没有表现出明显的规律。

(3)从BP神经网络的预测结果来看，该方法在所有站点几乎都能较好地反映出根区土壤湿度值的变化趋势，但在4个站点都存在极值预测不够准确的问题，此外，模型对于浅层土壤湿度的估算性能普遍比深层土壤湿度要好，50 cm以上土层的拟合精度高于50 cm以下土层。