基于植被指数的猕猴桃根域土壤水分反演影响因素研究

2022-02-08邓俊涛倪国威牛子杰潘时佳韩文霆

农业机械学报 2022年12期

张军邓俊涛倪国威牛子杰潘时佳韩文霆

(1.西北农林科技大学机械与电子工程学院，陕西杨凌 712100；2.西北农林科技大学中国旱区节水农业研究院，陕西杨凌 712100)

0 引言

水资源是农业生产中不可或缺的自然资源。植物水分胁迫状况反映了其水分需求状况，从而影响灌溉决策，但现有接触式测量方法难以大面积准确监测植株个体的水分胁迫状况。针对此问题，现阶段对于植物水分参数的部分研究已经确立了其与光谱、植被指数和温度的相关性，并建立了较为准确的反演模型。

植物表型监测的研究表明植物叶片形态具有判断植物胁迫状态的潜力，这种形态包括植物叶片多光谱反射率、叶片表面颜色、叶片热红外信息。基于植物叶片形态的无人机、卫星遥感作物监测模型也证明了其在遥感反演领域的应用潜力。在对监测模型的研究中，人工神经网络(Artificial neural network，ANN)算法的引入更进一步提升了遥感反演植物胁迫状态的能力[1-9]。

植物水分参数反演模型中植物的生理参数野外测量易受测量环境、测量方式影响，实验室测量方式耗费时间，但是对植物水分参数及其监测手段的研究为对土壤水分的反演提供了理论基础。研究证明在卫星遥感尺度上，土壤不同波长反射率也可以反映土壤干旱情况，这种反演方式也为低空无人机遥感反演土壤干旱情况提供了参考。相关研究已证明植物冠层反射率与土壤干旱程度具有一定关系，并建立了基于植被指数的土壤含水率反演模型，在此基础上有相关学者提出通过图像卷积进一步提取植物冠层反射率特征从而提升土壤含水率反演精度，但目前尚未对模型精度与影响因素进行研究[10-14]。

土壤含水率反演模型可以分析遥感图像上植被指数的空间异质性，对农田土壤水分状况进行评估，但是现有模型与研究依然面临以下问题：①果树相比于大田作物体积更大，根系与冠层分布更广。②冠层中存在不同年龄的叶片，这些叶片在反射率上有着较大的差异使得单个植株冠层植被指数具有空间异质性。③植株冠层之间会发生相互重叠，这种情况在非规范化果园中尤为明显，这使得在进行土壤含水率反演时不能像大田作物连续采集区域反射率，在评估植株根域土壤含水率(Root domain soil water content，RSWC)时需要讨论其数据采集的位置与大小。

针对以上问题，本文对猕猴桃冠层植被指数与RSWC进行相关性和显著性分析，并根据分析结果建立不同采样尺寸的植被指数-RSWC反演模型，最后通过绘制各个模型精度与采样尺寸的回归曲线研究采样尺寸对反演精度的影响并分析原因。

1 材料与方法

1.1 试验区域与试验设计

于2021年8月在陕西省宝鸡市眉县西北农林科技大学猕猴桃实验站(107°59′31″N，34°07′28″E，海拔643.22 m)进行试验，该地区属暖温带大陆性半湿润气候，年平均气温12.9℃，平均降水量609.5 mm，平均日照时数2 015.2 h，无霜期218 d。每年3—5月，回暖较快，秋季受冷空气影响，昼夜温差较为明显，是关中地区秋雨最多的区域之一；较为适合猕猴桃生长，但仍需人工对果园进行水分管理。该实验站种植有多种品种猕猴桃，本文试验植株品种为徐香，正处于果实膨大期(5—9月)，水分需求量大，树龄5～6年，株距3 m，行距4 m。试验地位置及试验区域标注如图1所示。

图1 试验地信息Fig.1 Experimental site information

试验设置了3组不同梯度(每个梯度包含4个植株)的试验组，统计不同处理下的土壤含水率空间变异性；并记录灌溉后5 d的数据，统计其时间变异性。3组试验组仅对土壤含水率进行控制，对应田间持水量的100%(充分灌溉)、90%(轻度亏缺)、80%(中度亏缺)，分别表示为T1、T2、T3，达到对应土壤含水率后使其土壤含水率自然下降，并在下降期间采集其土壤含水率与多光谱数据。

1.2 数据采集

试验数据采集于2021年8月3、4、6、9、10日的12：00—14：00，该时段猕猴桃果树较易出现水分胁迫现象且无人机获取的正射影像阴影较少。

遥感数据与地面数据采集同时进行，使用大疆Phantom 4 V2.0型无人机搭载Survy3型多光谱相机(MAPIR，美国)在高度15 m采集试验区域不同波长反射率的像元亮度(Digital number，DN)值。获取的图像使用MAPIR多光谱校准软件(MAPIR camera controler，MCC)与多光谱标定板校正；图像地面分辨率为0.007 m，根据果园的株行距试验区域按照600像素×400像素(约为4 m×3 m)进行分割，并作为初始采样尺寸对区域内的DN值进行计算，如图2所示。

图2 采集设备与采集点Fig.2 Acquisition equipment and acquisition points

为了获得正射影像，无人机在采集试验区图像时会向外拓展一部分距离，保证Pix4D Mapper软件在进行正射拼接时获得足够的重叠率。

土壤含水率采集深度为40 cm，处于水分活跃层(20～80 cm)，受到大气影响较小，但依然活跃，是猕猴桃根系的主要生长深度[15]。灌溉后在猕猴桃主干附近采集40 cm处的土壤剖面，并使用PR-3001-TRREC-N01型土壤水分检测仪(仁智测控技术公司)重复测量3次取平均值作为被测植株RSWC，并对异常值进行筛选去除。

1.3 数据处理

图像处理首先使用MCC对图像进行校正；其次使用Pix4D Mapper软件拼接试验区域正射影像；使用Python的PIL库对图像进行不同尺寸的中心剪裁，最后提取不同区域内的DN值并计算植被指数(Vegetation index，VI)，评估的植被指数包括绿色指数(GI)、修改型土壤调整植被指数(MSAVI)[16]、绿度归一化植被指数(gNDVI)[17]、归一化植被指数(NDVI)[18]、优化型土壤调节植被指数(OSAVI)[19]、重归一化植被指数(RDVI)[20]、土壤调整植被指数(SAVI)[21]、简单比值指数(SR)[22]。

1.4 统计分析

使用SPSS软件在95%置信水平(P<0.05)下使用单向单因素方差分析。为了比较不同处理的植株之间的植被指数，将5 d数据分为60组，RSWC为25.4%～30.9%，并使用Pearson系数对植被指数与RSWC共计9个变量进行相关性分析。

使用不同采样尺寸下的植被指数建模，评估采样尺寸对建模精度的影响。使用多层感知机(Multi-layer perceptron，MLP)算法进行建模，MLP的隐节点采用输入向量与权向量的内积作为激活函数的自变量，激活函数采用ReLU函数。各参数对网络的输出具有同等地位的影响，因此MLP是对非线性映射的全局逼近，该模型以冠层植被指数作为输入，使用MLP模型预测RSWC。

最后，采集试验植株新生与成熟叶片叶绿素SPAD值，统计分析同一植株的叶绿素含量方差，讨论各个参数对反演精度的影响原因。

2 结果与分析

2.1 植被指数与根域含水率变化

于2021年8月3日进行灌溉，灌溉后让植株自然生长。根据周边气象站数据记录，试验期内日平均空气相对湿度(Mean relative humidity，MRH)与测量时空气相对湿度(Relative humidity，RH)总体呈上升趋势；测量时空气相对湿度最低值出现在8月6日，为46%，试验期间空气相对湿度变化如图3所示。

图3 2021年8月日平均空气相对湿度与采集时间点空气相对湿度Fig.3 Daily mean relative humidity in August, 2021 and relative humidity at collection time point

根据记录，灌溉后的RSWC自然下降，第3天下降至25.4%，为了防止水分亏缺对猕猴桃植株造成不可逆损伤，在8月6日后对猕猴桃果园进行少量连续灌溉使其逐渐恢复；之后3组试验组RSWC持续回升。T3(中度亏缺)组RSWC在整个试验期间一直较低，最高为27.4%(8月10日)，同日T1、T2组大致为30.0%。T3组RSWC最低为25.4%，同日T1、T2组RSWC也达到试验期间最低，大致为28.0%。8月6日和10日，T1和T2组RSWC十分接近，这是由于RSWC受到作物自身生理作用、冠层结构以及气候因素多重影响，这使得灌溉后自然下降的试验组间可能出现相近的测量数据，这也是果园内RSWC时间与空间异质性的一种体现。同时，也正是复杂因素导致的水分异质性，才有必要进行果园内水分精细管理。试验期间不同组RSWC变化如图4所示。

图4 2021年8月不同处理试验区域平均RSWC日变化Fig.4 Daily variations of average RSWC in different test areas in August, 2021

根据不同采样尺寸下获取的冠层DN值，计算不同的植被指数。统计5 d内3种不同处理的12个植株，共获取60组RSWC与植被指数样本，每组样本包含8个植被指数；60组样本的植被指数与RSWC变化如图5所示。从图5可以看出，除RDVI、SR、GI外，其他的植被指数均处于0～0.4之间；在RDVI、SR、GI中，SR与GI的变化幅度较小，分别处于1.0～1.5与0.5～1.0区间；而RDVI则处于1.0～2.1之间，变化跨度大，且不同处理与时间的样本之间变化幅度较大，这种现象表明RDVI对土壤含水率的变化较为敏感，将有利于反演模型的建立；植被指数GI与其他植被指数变化趋势相反。

图5 所有样本RSWC与8个植被指数变化曲线Fig.5 Changes of RSWC and eight vegetation indices for all samples

2.2 不同植被指数与RSWC关系

为了研究不同植被指数对RSWC反演的影响，使用SPSS软件对正则化后的不同植被指数与RSWC进行Pearson相关性分析。

不同植被指数与RSWC的Pearson相关系数如表1所示，其中只有RDVI的Pearson相关系数大于0.6，表示其与RSWC呈强相关性；GI、SR与RSWC的相关系数绝对值小于0.5，其中SR在所有植被指数中与RSWC相关系数最小，呈弱相关性；其余植被指数的相关系数相互之间差异不明显。

表1 不同植被指数与RSWC的Pearson相关系数Tab.1 Pearson correlation coefficient between different vegetation indices and RSWC

为了排除植被指数之间的互相关，通过进行两两参数之间的Pearson相关性分析，获得Pearson相关性矩阵，如图6所示。横纵轴交叉处即为Pearson相关系数。GI与RSWC及其他植被指数的Pearson相关系数小于0，呈负相关；其中GI、SR、gNDVI与RDVI的互相关性相较于其他植被指数较弱(Pearson相关系数低于0.85)。

图6 Pearson相关系数矩阵Fig.6 Pearson correlation coefficient matrix

通过SPSS的单因素方差分析计算植被指数与RSWC之间的显著性，结果如表2所示；除SR外其他植被指数均与RSWC呈现出较强的显著性(P<0.05)，其中RDVI与RSWC的显著性最强，其次是gNDVI、MSAVI、GI，这些植被指数的P值均小于0.04。

表2 植被指数与RSWC显著性Tab.2 Significance result of vegetation index and RSWC

根据相关性与方差分析结果，考虑到植被指数与RSWC的相关性与显著性以及植被指数之间的互相关性，本文使用RDVI作为反演模型的主要植被指数,使用RDVI与RSWC建立基于多层感知机反演模型。

2.3 不同采样尺寸与RSWC相关性

使用RDVI作为输入(两列输入向量均为RDVI)建立多层感知机反演模型时，发现不同冠层采样尺寸获得的RDVI数据组对同一参数的人工神经网络模型反演精度会造成影响。在每个植株冠层初始采样尺寸(600像素×400像素)的基础上对采样范围进行变换，如图7所示，图中b为采样区域宽度，h为采样区域高度，L为采样区域对角线长度，黑框为初始采样尺寸。采样区域宽度计算公式为

式中i——比例系数

bi——经过比例系数i裁剪后的采样区域宽度不同尺寸下的植被指数将被用来训练同一个MLP，以比较不同植被指数下对反演模型精度的影响，MLP的参数分别为：层结构为(2，2，1)，优化器使用L-BFGS算法，最大训练世代为2000代，激活函数为ReLU函数，学习率使用动态更新策略。

图7 采样尺寸变换过程示意图Fig.7 Sampling size transformation process

模型训练集包括所有不同处理与观测时间数据的60%，剩余40%作为测试集。其中训练集每一组数据由两个植被指数与对应土壤含水率组成，用于神经网络模型训练，测试集由具有未知标签的植被指数组成，用于评估模型预测效果。

对使用不同采样尺寸下的植被指数建立的反演模型决定系数R2进行统计，获得不同采样尺寸与模型精度，如表3所示。其中，采样面积A=hb。

表3 采样尺寸与模型精度Tab.3 Sampling size and model accuracy

从R2在L、A上的分布可以看出，模型结果的决定系数随着采样尺寸的减小，先增后减，但是不同采样尺寸建模下的模型均方根误差(RMSE)没有显著差别，这说明会存在最佳采样尺寸使得模型精度达到最大。分别建立采样区域对角线长度与模型决定系数二次、三次拟合曲线和采样区域面积与模型决定系数二次、三次拟合曲线L-R2和A-R2，如图8所示，可以看到L-R2与A-R2散点图的不同阶次多项式拟合曲线获得的最高决定系数分别为0.993与0.991，说明采样尺寸对模型结果的决定系数有很大影响。

图8 采样区域参数与模型决定系数二次与三次拟合曲线Fig.8 Quadratic and cubic fitting curves of sampling area parameters and model determination coefficient

在L-R2散点图中，模型最大决定系数出现在L为360像素左右，此时b大致为300像素；在A-R2散点图的二次拟合曲线中，模型最大决定系数出现在A为50 000像素左右，三次拟合曲线则A为64 000像素左右，此时b为273～309像素。根据以上分析，b在273～300像素之间时建模精度达到最大，由于地面分辨率为0.007 m，所以这个范围的宽度在实际中应该在1.953～2.170 m之间，面积在2.540～3.038 m2之间。

3 讨论

猕猴桃作为喜水植物对水分的需求量比其他果树多[23]，最适宜猕猴桃种植地区的降雨量要达到800～1 200 mm，而陕西关中地区的降雨量不能满足这一要求；在自然条件下，猕猴桃植株会逐渐出现水分胁迫，这个现象在果实膨大期尤为明显。这种由于土壤含水率下降导致的胁迫可以通过叶水势来进行评估，其中正午叶水势反映植株水分最大亏缺程度。众多关于植被指数与水分参数相关性的研究[4,6,24-26]以及叶水势与土壤含水率相关性的研究表明，植物在水分胁迫时的植被指数、叶水势以及土壤含水率三者之间存在较强的相关性[27]，这也使得通过植株冠层对土壤含水率进行反演成为了可能；本文将8种植被指数与土壤含水率进行了相关性与显著性分析，确定了RDVI为表征RSWC的最佳指标(相关系数为0.744，P<0.05)，结果与张智韬等[13]的剔除土壤背景的RSWC研究结论类似，通过冠层植被指数可以较好地反映出一定深度土壤的含水率；SEO等[14]通过卷积神经网络建立冠层图像与土壤含水率的研究也说明了这一点，但是针对于不同根系深度的植株还有待研究。

猕猴桃作为藤本果树与大田作物有着较大的区别，这主要体现在：①单株经济价值更高。②植株体型庞大，根系更深覆盖面积更广。③叶片宽大相互交叠，新老叶片的水分状态与颜色差异较大。在对猕猴桃冠层的DN值进行采集时，发现采集区域不同将对结果造成较大的影响，从而影响植被指数与建模精度。对猕猴桃冠层的NDVI图像反归一化后进行了假彩色处理，如图9所示，图中黑色圆圈代表植株的主干位置。

图9 NDVI反归一化后的假彩色图像Fig.9 False color images after inverse normalization of NDVI

从图9可以发现，猕猴桃冠层的NDVI分布并不均匀，如图9冠层顶部的植被指数普遍大于冠层边缘的植被指数，且呈向中心聚集的趋势；部分冠层并没有完全覆盖土壤，且土壤的NDVI与冠层的植被指数有较大差异；如图9b果园内部猕猴桃植株之间会发生藤蔓缠绕与冠层重叠，重叠部分NDVI较高，但是无法判断属于哪一个植株。

同时在8月10日，使用SPAD502(Konica Minolta Japan)随机对每个植株的3片叶子的上中下部分进行了3次测量，获得每个植株3片叶子的平均SPAD值，该测量值表示叶片叶绿素浓度，如图10 所示。记录叶片包括冠层顶部成熟叶片与边缘新生叶片，统计数据表明，不同处理植株之间，叶片SPAD值差距最大为21.67；即使是同样处理的植株，叶片SPAD值也有很大差异，最大差值也达到21.17，标准差为6.403；同一个植株不同区域叶片SPAD最大差值为12.26。通过观察，成熟叶片与新生叶片在颜色，质地上也存在较大差异。研究表明植物叶片通过叶绿素和叶黄素等植物色素吸收可见光波段的大部分辐亮度[24]，植物受外界胁迫时会导致叶绿素色素的损害，光合作用效率下降，这也使得叶片吸光度和反射率的变化[28]，改变了自身冠层的反射率格局，使可见光波段反射率增加，近红外波段反射率降低，从而改变植被指数。

图10 所有植株叶片SPAD平均值Fig.10 Mean values of SPAD in leaves of all plants

根据对假彩色图像中NDVI的分布与叶片SPAD值的统计，推测以下原因使得采样尺寸对模型精度造成了影响：①新生与成熟叶片的反射率差别及两者在冠层上的不均匀分布导致了冠层植被指数分布不均匀，新生叶片植被指数较低，而成熟叶片的相对较高，这也使得在变换采样尺寸时获得的植被指数平均值不同；越是向冠层中心集中，植被指数均值也会更加趋近于成熟叶片均值。②土壤像元植被指数与冠层像元植被指数有着较大差异，部分猕猴桃冠层没有完全遮盖土地，与同组遮盖完全的植株相比，提取到的植被指数均值会更小；越是接近冠层中心，裸露的土地越少，从而也减弱了土壤像元对冠层植被指数的影响。③猕猴桃植株属藤本类果树，果园内棚架相连使得猕猴桃藤蔓相互缠绕重叠，在大采样尺寸下，小冠层植株的植被指数均值可能包括周围其他处理植株冠层的植被指数；采样范围接近冠层中心时，将避免其他冠层对目标植株的影响。④当采样尺寸过小时实际上是对冠层顶部中心的几片叶子进行植被指数采集，叶子自身状态(例如病虫害等)将会对植被指数均值造成极大的影响，这也是在缩小采样尺寸到一定程度后模型精度下降的原因。

经过试验，果实膨大期徐香猕猴桃的植被指数采集面积在2.540～3.038 m2之间较为合适，本研究中并没有考虑采样形状与植株品种的影响。

ROMERO等[6]结合机载高光谱获得的高分辨率图像与人工神经网络，建立了葡萄园的叶水势反演模型，该方法取代了传统的压力室法测量叶水势，研究证明了通过无人机高光谱监测农田含水率的可行性与ANN在遥感数据挖掘中的巨大潜力。本文通过MLP神经网络建立了植被指数-RSWC反演模型，使用2.540～3.038 m2之间采集的RDVI数据训练的模型表明，基于MLP网络的土壤含水率反演模型可以在所有数据拟合上取得较高的精度(R2为0.638，RMSE为0.016)，表明了模型的精确性与泛用性。人工神经网络等算法仍处于发展阶段，随着算法更替将会使得模型的精度与泛化能力进一步提升，这也将是未来的一个研究方向。