无人机热红外遥感反演玉米根域土壤含水率方法研究

2021-04-04陈俊英周永财崔文轩

节水灌溉 2021年3期

杨帅，陈俊英，周永财，崔文轩，杨宁

（1.西北农林科技大学水利与建筑工程学院，陕西杨凌712100；2.西北农林科技大学旱区农业水土工程教育部重点实验室，陕西杨凌712100）

0 引言

快速、准确地获取区域土壤含水率状况对于判断作物水分胁迫状况、指导农业生产具有重要意义[1]。传统的土壤含水率监测方法众多[2，3]，但都局限于点测范围，没有对区域土壤含水率状况做出整体反映。1963年，Tanner[4]首次发现冠层温度能够作为指示作物水分的有效指标并首次使用红外测温仪来量化冠层温度差异与作物水分胁迫之间的关系。随后，Idso[5，6]利用冠层温度与大气温度的差值建立了作物水分胁迫指数CWSI的经验公式。Jackson[7]和Jones[8]相继完善并补充了CWSI的定义公式。此后一大批国内外学者开始用CWSI反演土壤含水率[9-12]。近年来，随着无人机遥感的快速发展，利用无人机遥感平台搭载热红外成像仪快速获取区域面状温度进而判断区域水情和作物水分胁迫状况逐渐成为热点[13-15]，并为监测土壤含水率提供了新的思路。Agam 等[16]探究了CWSI对太阳辐射的敏感程度。孙圣等[17]利用无人机热红外技术测定了北方核桃园区的土壤含水率分布。张立元等[18]将无人机遥感与田间空气温湿度数据结合建立了监测大田玉米水分胁迫的经验模型。张智韬等[18，19]探究了基于无人机尺度构造的覆盖度和综合温度指数对不同深度土壤含水率的响应关系。边江等[20]以遥感热红外波段反演的不同农作物土壤体积含水率与灌区的土地分类影像为依据，快速推求灌区灌溉需水量。尚晓英等[21]将3 种CWSI与棉花根域土壤含水率建立模型，确定了最佳反演模型。

要实现从无人机热红外图像中准确提取作物冠层温度信息，剔除土壤背景干扰必不可少[18,23]。目前，从热红外图像中剔除土壤背景，提取特定地物信息的常见处理有两种，一种是直接从热红外图像中提取，主要思路为：①基于热红外图像灰度分割阈值来剔除干扰背景；②基于热红外图像特定地物形态边缘特征来提取[22-24]，这两种思路都对热红外图像的分辨率有较高要求[25]。另一种处理是结合可见光和热红外图像[26-28]，利用可见光图像中各波段对地物特征的不同反映来辅助提取热红外图像中的有效温度信息，这种方法考虑了地物的形态学特征，对低分辨率的图像效果较好[19]，但对图像的配准精度要求较高。两种处理都有其局限性。在实际运用中采用一种处理方法的较多，而同时使用两种处理方法并对提取效果进行精度评价的研究较少。

由于拔节期玉米枝叶细长且相互交错，在热红外图像上表现为线条特征繁多凌乱，且伴随有较大阴影覆盖，传统阈值分割方法没有考虑到叶片边缘特征以至于出现像素误分从而影响提取精度。本文以拔节期玉米为研究对象，采用可见光和热红外图像相结合的RGRI指数法，将可见光图像中玉米形状轮廓特征明显的优势和热红外图像的温度敏感优势结合，剔除热红外图像中的土壤背景，提取玉米冠层温度，同时设置Otsu阈值处理法和不剔除土壤背景处理作为对照，将基于3种方法提取的冠层温度与实测温度作对比，然后用基于这3种方法建立的CWSI指标反演不同深度的土壤含水率，分析并评价基于何种方法建立的CWSI指标反演效果最佳。

1 材料与方法

1.1 试验地概况

试验地位于陕西省杨凌示范区西北农林科技大学旱区节水农业研究院（108°4'E，34°42N'），地处关中平原，属于东亚暖温带半湿润半干旱气候区，具有夏热多雨、秋热凉爽多连阴雨等明显的大陆性季风气候特征。年均降水量635.1 mm，平均蒸发量为993.2 mm，土壤质地为中壤土。

1.2 试验设计

试验地划分为12 个小区，设置4 种水分梯度处理（T1、T2、T3、T4），分别以95%～100%、80%、65%、50%田间持水量为灌水上限，以80%、65%、50%、40%田间持水量为灌水下限，每个处理设置3 个重复试验。每个小区面积为4 m×4 m，小区间距为1.5 m，每个小区采用行播的方式播种7行玉米种子，并以滴灌方式进行灌溉。试验开始前2 d 内对小区进行控水处理，以确保水分梯度。

1.3 数据采集

在玉米拔节期阶段，试验分别于2019年7月27日（拔节前期）、7月31日（拔节中期）、8月2日（拔节后期）三日13：00进行，通过大疆M600六旋翼无人机平台搭载可见光相机和大疆禅思XT 热红外测温成像仪进行图像信息采集，飞行高度为20 m，期间地面同步进行手持测温枪采集玉米冠层、水桶和辐射定标板温度信息。飞行任务结束后，通过打土钻的方式获取各小区内10、20、30、40、60 cm 深度的土样，放入铝盒内采用烘干法来测定各深度土样的质量含水率。

1.4 无人机图像处理

1.4.1 热红外图像校准和灰度值转换

获取的原始热红外图像，首先要导入Flir tools软件并根据地面手持测温枪实测的水温和辐射定标板温度进行校准，设置辐射率为0.96。将校准后的温度信息导至ENVI 软件中与原图像进行波段叠加。再通过ENVI 拓展工具对叠加图像中灰度波段和温度波段进行统计，得到灰度值与温度值的线性转换关系式。

1.4.2 基于RGRI指数法的图像处理

可见光图像中的绿波段对于植物绿反射较为敏感，红波段穿透性强，能较好减轻阴影带来的干扰，因此通过二者比值构造RGRI指数能较好区分玉米前景和土壤、阴影等背景，对于RGRI指数的计算见下式。

式中：bandred为可见光图像中的红波段；bandgreen为可见光图像中的绿波段。

对于获取的可见光图像，其视场角较广，包含多个小区，首先选取出具有各试验小区正射视角的可见光图像作为各小区的代表图像进行裁剪，裁剪尺寸为512×512像素。将裁剪后的可见光图像和校准后的热红外图像导入Arcgis软件中定义相同地理坐标系，接着在ENVI 软件中根据地物特征进行配准，并将配准后的可见光图像进行波段运算，得到RGRI指数图。通过RGRI统计直方图（图1）可以看到玉米像素和土壤像素二者有较为明显的RGRI分界值，由此确定玉米和土壤的分类阈值。

通过Arcgis 软件根据确定的分类阈值对RGRI指数图进行二值化操作，把玉米像素定义为1，把土壤像素定义为0，并对玉米像素进行提取，接着通过栅格转面操作得到具有玉米像素边界特征的矢量图层（图2）。最后将该矢量图层与配准后的热红外图像导入到ENVI软件中进行掩膜处理（图3），统计得到小区内玉米像素灰度值表，根据前面算出的灰度温度转换关系式即可得到相应的玉米冠层温度值表。

本研究为了减少配准和RGRI分类时误将土壤像素划归到玉米像素带来的误差，将得到的玉米冠层温度直方图剔除前后1%的温度值（图4），然后统计得到玉米冠层平均温度。

1.4.3 基于Otsu阈值法的图像处理

Otsu 阈值法又称最大类间方差法，它可以对灰度图像进行阈值分割，使得前景与背景图像的类间方差最大，从而达到较为理想的分割效果。

本研究对经过校准的热红外图像按小区进行裁剪，裁剪尺寸为480×480 像素。将裁剪后的各小区热红外图像导入MATLAB 软件中，通过编程实现玉米像素与土壤像素的分割，得到玉米像素温度矩阵。为减少图像阈值分割过程中误将土壤像素识别为玉米像素带来的干扰，对于获得的玉米温度采取剔除最大1%温度的处理，由此得到的平均温度作为小区玉米冠层平均温度。

1.4.4 不剔除土壤背景的图像处理

不剔除土壤背景将会导致图像提取温度偏高，本研究将其作为对照组。将校准后的热红外图像导入至Flir tools 软件中，通过软件自带的框选工具，框选整个小区，利用软件的均值统计工具得到框选区内的平均温度作为小区玉米冠层的平均温度。

1.4.5 水分胁迫指数CWSI

本研究作物水分胁迫指数（CWSI）的计算参考JONES等[8]的简化公式并借鉴了张智韬等[19]的研究。本研究中选取12 个小区中冠层温度均值最大值Tmax加上5 ℃作为干参考面，以小区中冠层温度均值最小值Tmin减去2℃作为湿参考面。具体计算见下式。

式中：TC为作物冠层平均温度；Tmax为小区作物冠层平均温度最大值；Tmin为小区作物冠层平均温度最小值。

1.5 评价指标

本研究选取决定系数R2和均方根误差RMSE作为回归评价指标，R2越接近1，RMSE越接近0，模型反演效果越显著。采用F检验来统计样本显著性，若F0.05≤F≤F0.01，即0.01＜P≤0.05，则各处理间差异显著；若F≥F0.01，即P≤0.01，表示各处理间差异极显著。

2 结果与分析

2.1 基于热红外图像提取温度与实测温度相关性分析

将基于图像提取的冠层温度与实测温度进行线性拟合，如图5所示。从相关性来分析，在整个拔节期3 个阶段，基于RGRI指数法提取的温度与实测温度拟合效果最好，R2分别为0.909、0.891、0.828，其次是基于Otsu 阈值法，R2分别为0.862、0.865、0.875，而不剔除土壤背景提取温度表现最差，R2为0.787、0.645、0.676。对比均方根误差RMSE，按照RGRI指数法、Otsu 阈值法、不剔除土壤背景的顺序，在拔节前期，三者RMSE分别为0.65、1.10、2.46 ℃，在拔节中期，三者RMSE分别为0.58、0.57、1.00 ℃，在拔节后期三者RMSE分别为0.54、0.89、1.31 ℃。整体来看，RGRI指数法表现最好，误差最小，其次是Otsu 阈值法，而不剔除土壤背景误差最大。对比温度拟合线趋势，RGRI指数法和Otsu 阈值法二者趋势线较低且互有交叉，提取温度与实测温度较为接近，而不剔除土壤背景趋势线较高，由于土壤背景的干扰使得提取温度与实测温度相差较远，温度偏高。

2.2 不同水分梯度下CWSI指标变化趋势

根据式（2）计算出基于3 种方法的CWSI值，为了对比3种方法提取CWSI指标的差异性，本研究取每种水分梯度下3个重复小区的CWSI均值作为该水分梯度下的CWSI值。

如图6所示，在玉米拔节期3 个阶段，基于3 种方法得到的CWSI指标均在整体上表现出随着土壤含水率的减少而呈现增大的趋势，这表明CWSI指标与土壤含水率整体呈现出负相关的响应趋势。在拔节前期，T1 和T2 水分梯度下3 种方法得到的CWSI值较低这可能是由于试验前部分梯度小区要进行补水处理，而作物吸收水分时出现滞后效应导致的。随着土壤水分的消耗，到了拔节中期和拔节后期，T1 和T2 水分梯度下3 种方法得到的CWSI值均较拔节前期有了增加。整体来看，CWSI在T3 和T4 水分梯度下增长较快，而在T1 和T2 水分梯度下增长较缓，这表明相比于土壤水分充足状态，CWSI对于水分亏缺条件下作物的水分胁迫状况较为敏感，对于该状态下土壤含水率的变化应有较好的指示作用。从CWSI指数在不同水分梯度下的变化幅度来看，与拔节前期相比，拔节中期和拔节后期CWSI指数的下界增大，上界减小。这可能是因为随着玉米生长，玉米对于水分的依赖性逐渐减弱。

2.3 CWSI与不同深度土壤含水率之间的相关关系

本研究将基于3 种方法得到的各小区CWSI值与玉米拔节前、中、后期不同深度土壤含水率数据在SPSS 软件中采用线性回归方法构建模型，判定哪种方法得到的CWSI反演土壤含水率效果最好，同时探究玉米不同时期下CWSI对土壤含水率的最佳响应深度，对应相关关系见表1～表3（表中Smc代表土壤含水率）。

表1 拔节前期CWSI与土壤含水率的回归模型

表2 拔节中期CWSI与土壤含水率的回归模型

由表1～表3可知，3 种方法得到的CWSI值与土壤含水率的拟合系数均为负值，呈现负相关关系。从3 种CWSI与各深度土壤含水率相关性的角度来分析。在拔节前期，CWSIRGRI与0～10、0～20、0～30、0～40 cm 土壤含水率相关关系较好，R2分别为0.513、0.632、0.761、0.675，CWSIOtsu效果次之，R2分别为0.452、0.496、0.592、0.640，CWSIsoil表现效果最差，R2为0.329、0.405、0.536、0.629，而在0～60 cm 深度内，CWSIsoil的R2达到0.648，高于CWSIOtsu的0.623 和CWSIRGRI的0.591。对于CWSIRGRI，其在拔节中期和拔节后期五个深度内的表现优于CWSIOtsu和CWSIsoil，R2均大于0.5，效果显著(P＜0.01)。对于CWSIOtsu和CWSIsoil，在拔节中期0～20、0～30 cm深度内，CWSIOtsu表现较好（R2为0.683、0.732），在0～10、0～40、0～60 cm则是CWSIsoil拟合效果更优（R2为0.687、0.595、0.530），分析原因为在玉米拔节中期，部分地块热红外图像中的玉米像素与土壤像素灰度值较为接近，此时依靠Otsu 方法并不能很好的确定二者分类阈值，导致提取的玉米像素中掺杂了较多土壤像素，因此拟合效果较差。在拔节后期5个深度内，CWSIOtsu的R2整体高于CWSIsoil。因此整体来看，CWSIRGRI对各深度土壤含水率反演效果最好，CWSIOtsu次之，CWSIsoil效果最差。

从CWSI对不同深度土壤含水率的响应程度来看，选取CWSIRGRI作为最优CWSI指标，其在玉米拔节前、中、后期五个深度内的R2均表现出先增大后减小的趋势，且减小幅度大于上升幅度，并在0～30 cm 处达到最大值，这可能是因为此时玉米根系主要集中在0～30 cm 深度内，并且此时玉米根系毛管大部分应集中在土壤20～30 cm 处，因此该深度对土壤含水率的响应最为敏感，而分布于土壤30 cm 深度下的玉米根系由于较少且分布远不及0～30 cm 深度土壤内的根系分布，故而CWSIRGRI在0～40、0～60 cm 深度下的R2出现了明显的下降趋势。总体来说，CWSI对0～30 cm 深度土壤含水率反演效果最好，对0～60 cm土壤含水率反演效果最差。

2.4 模型验证评价

本研究对3 种处理下的CWSI与不同深度的土壤含水率模型进行验证，以拔节前期和拔节中期数据建立函数关系，以拔节后期数据进行验证，通过对比决定系数R2以及均方根误差RMSE对建立的3种模型进行评价（表4）。

从模型验证表可以看出，整体R2均在0.5 以上。对于同一深度，RGRI指数法的决定系数最大，Otsu 阈值法次之，不剔除土壤背景效果最差。3 种方法在5 个土壤深度内的R2呈现单峰分布且都在0～30 cm处达到最大值（R2分别为0.801、0.658、0.552）。基于RGRI指数法的均方根误差RMSE最小，均在0.1附近。3 种处理方法在不同土壤深度的RMSE表现呈现单峰波谷分布且都在0～30 cm 处达到最小值，这可能间接表明在玉米拔节期构建的CWSI指标能较好地反演0～30 cm 土壤深度内玉米根系活动层的土壤含水率。

3 讨论

本研究采用的RGRI指数法很大程度上依赖可见光图像与热红外图像的配准精度。在试验过程中，无人机平台无法同时搭载可见光相机与热红外相机，这就造成可见光图像与热红外图像的获取过程会出现短暂的时差，在这时段内光照条件的突然变化以及风力扰动带来玉米叶片的变化都会给配准精度带来很大干扰。为了减少干扰，下一步研究将把研究区选定为大尺度区域以期通过增大尺度的方式来均摊干扰，提高整体精度。

CWSI的建立有很多种方法[5-8]，张智韬等[20]通过对棉花叶片涂凡士林和洒水的方式来营造干湿参考面，进而确定CWSI模型的上下基线。本研究建立的CWSI模型是基于小区平均冠层温度来确定的干、湿参考面[27]，究竟哪种更适合还有待进一步探究。本次研究选取的数据样本只选取了13：00一个时刻，且有一定的时间间隔，这也可能会对试验结论的精确性产生一定的影响，下一步研究将考虑对连续时间内不同时刻的数据样本进行分析。

随着玉米生长阶段的变化，玉米根系分布也会不同，导致CWSIRGRI的最佳土壤含水率反演深度也会发生变化，具体深度是多少还有待进一步探究。针对不同水分梯度对反演精度的影响，本研究也缺乏考虑，在今后的研究中还需要加以探究。