剔除土壤背景的棉花水分胁迫无人机热红外遥感诊断

张智韬 边 江 韩文霆 付秋萍 陈硕博 崔 婷

(1.西北农林科技大学旱区农业水土工程教育部重点实验室， 陕西杨凌 712100； 2.西北农林科技大学水利与建筑工程学院， 陕西杨凌 712100； 3.新疆农业大学水利与土木工程学院， 乌鲁木齐 830052)

0 引言

土壤水分是影响农作物生长的关键因素，也是指导灌溉的理论依据，利用热红外非接触性诊断作物水分状况成为当前热点[1-2]。1963年，TANNER[3]首先研究发现冠层温度可以指示作物水分胁迫，随后，IDSO等[4-6]利用空气温度Ta和冠层温度Tc之差与上下基线建立经验法水分胁迫指数CWSI。但经验法CWSI上下基线存在误差和不确定性，所以JACKSON等[7]根据冠层能量平衡推导出CWSI中上下基线的理论计算公式，用于小麦的水分诊断研究，取得较好的诊断精度，并一直作为水分诊断的可行方法[8-10]。由于理论公式所需参数较多，JONES等[11-12]将CWSI计算方法进行简化，利用干湿参考面确定上下温度极限值，使CWSI更加容易指导生产实践，而简化法CWSI评估棉花水分亏缺状况也取得了一定进展。MERON等[13]研究了冠层阳光直接照射叶片温度、冠层阴影叶片温度、土壤温度对热红外诊断棉花水分状况的影响。XU等[14]进一步将CWSI与蒸腾速率、净光合作用、气孔导度的关系作了深入的分析。近年来，作物水分胁迫指数CWSI在国内也逐渐成为研究的热点[15-18]，张立伟等[19]分析了冠层温度、冠气温度比、冠气温度差与玉米根系土壤含水率的相关关系，蔡甲冰等[20]将冠层温度作为自动灌溉决策系统的核心理论基础。

以上研究都是基于地面点测方式获取的冠层温度，然而将点状温度转为面状冠层温度具有一定误差和不确定性。随着低空无人机遥感平台的发展，快速获取高时间和空间分辨率的热红外图像成为可能,有效解决了面状冠层温度采集的问题。HOFFMANN等[21]利用无人机热红外技术采集冠层温度信息，并且结合多种植被指数共同诊断大麦的水分胁迫状况。BALUJA等[22]利用无人机热红外图像计算CWSI，用于评估葡萄园水分空间状况的差异性。RUD等[23]研究发现，无人机热红外图像计算的CWSI与马铃薯土壤含水率具有显著的相关关系，EGEA等[24]将CWSI应用在高密度的橄榄园及桃园的土壤水分监测中，并取得较好的研究成果。然而无人机热红外图像受到土壤背景的干扰，作物水分胁迫状况的监测精度普遍较低，HAN等[25]利用EM算法研究冠层温度的分布及土壤背景的干扰情况，但并没有将土壤背景进行有效剔除。Canny算子是一种图像边缘检测技术，可以提取不同物质类型的边界像元，PARK等[26]应用Canny边缘检测技术优化了水分胁迫指数CWSI干湿参考温度的计算，但利用Canny边缘检测算法剔除热红外图像土壤背景的研究较少。

本文以4种不同水分处理的花铃期棉花为研究对象，采用无人机高分辨率热红外成像仪连续监测，并分别采用二值化Ostu算法和Canny边缘检测算法对热红外图像进行掩膜处理，实现对土壤背景的剔除，然后分别计算两种热图像处理算法和包含土壤背景的作物水分胁迫指数CWSI，并研究3种方法计算的CWSI与棉花叶片气孔导度Gs的关系模型。

1 材料与方法

1.1 研究区概况及试验设计

1.1.1研究区概况

试验田位于陕西关中平原腹地，西北农林科技大学中国旱区节水农业研究院(108°4′20″E，34°17′42.17″N)，该地区平均年降雨量640 mm，平均蒸发量为993.2 mm，海拔525 m，地下水位埋藏较深，属于暖温带季风半湿润气候区，土壤质地为中壤，平均田间持水率23%～26%，凋萎含水率8.6%(均为质量含水率)，平均土壤干容重1.44 g/cm3。试验田开阔平坦配备遮雨架进行控水措施，并设有可移动式小型气象站，试验期间天气晴朗，总辐射较强,无阴雨天气干扰。

1.1.2试验设计

试验设计4个水分梯度I1(50%田间持水量)、I2(65%田间持水量)、I3(80%田间持水量)、I4(对照组95%～100%田间持水量)，每个水分处理设置3个重复试验，共12个小区，小区面积4 m×5 m。小区配有水表，采用滴灌的方式进行灌溉，棉花每行放置一条滴灌带，滴头湿润半径20 cm，并且小区之间设有宽垄，以尽量避免各个小区土壤水分侧向入渗的影响(图1，2017年7月11日13:00采集)。图中小区编号I11代表1区的1个重复，其他编号类似。本研究选用西北农林科技大学农学院培育的新品种西农棉1008，于2017年4月19日人工覆膜播种，施用225 kg/hm2磷肥和375 kg/hm2尿素作为底肥，棉花苗期不做控水处理。试验时间为2017年7月11—15日，连续监测5 d，棉花处在花铃盛期，是棉花结铃生长发育的关键生育期。试验开始前棉花水分胁迫出现明显状态，参考设置合理的计划湿润层深度，严格按照4个水分处理进行灌溉，试验期间不进行任何灌水，每天于09:00、13:00、17:00进行无人机冠层温度采集。

图1 试验小区可见光影像Fig.1 Visible light image of experimental area

1.2 无人机热红外系统

1.2.1无人机系统

无人机系统采用大疆公司Matrice600搭载禅思XT热红外(TIR)相机(FLIR系统，Tau2机芯)获取棉花花铃期的冠层温度影像。Matrice600六旋翼无人机具有飞行姿态稳定且续航能力长等特点，最大起飞质量15.1 kg，最大可承受8 m/s风速(图2)。TIR相机工作的波段为7.5～13.5 μm，空间分辨率为640像素×512像素，温度灵敏度0.05℃，镜头焦距19 mm，视场角32°H×26°V(图3)。为了拍摄分辨率较高且质量较好的热红外影像，无人机悬停15 m高度采集棉花冠层的热红外影像，图像分辨率为0.011 m。

图2 大疆Matrice600无人机Fig.2 DJI M600 UAV

图3 大疆禅思XT热红外成像测温仪Fig.3 DJI Zenmuse XT thermal infrared imaging thermometer

1.2.2热红外图像的校准

热红外图像采用地面实测叶片温度与水温进行校准。无人机采集棉花冠层温度时，每个试验小区均匀选取阳光直接照射棉花叶片,地面同步利用手持热红外测温仪测量其叶片温度，最终取所测温度的平均值，并采集试验区放置的器皿中水的温度。将无人机采集的热红外图像输入FLIR Tools(禅思XT热红外影像处理软件)，设置辐射率为0.96，将叶片平均温度和水温作为参考温度，通过对比参考温度和热红外图像上相应像元温度的方法进行红外温度图像的校准，即可得到温度较为准确的无人机热红外温度图像。

1.3 生理指标数据采集

棉花叶片的气孔导度(Gs，mol/(m2·s))使用美国Licor公司LI-6400型便携式光合测定仪测定。气孔导度的测定与无人机采集热红外图像同步进行，每个小区选取3株棉花样本(倒三叶片)测量，每株重复3次。

1.4 气象数据采集

使用自行搭建的可移动式小型气象站自动收集气象数据，主要气象要素有：大气温度、相对湿度、风速、太阳净辐射、土壤热通量等，数据采集间隔为30 min，试验期间(2017年7月11—15日)的主要气象数据如表1所示。

表1 主要日平均气象因子Tab.1 Values of major daily average meteorological factors

1.5 热红外图像剔除土壤背景的方法

1.5.1二值化Otsu算法

无人机热红外成像测温仪获取的图像具有两类特征像元，一类为目标冠层像元，另一类为土壤背景像元。二值化Ostu算法(最大类间方差法)是灰度图像阈值分割的经典自适应阈值算法，在图像处理领域得到广泛的应用[27]，其核心理论公式为

σ2(k)=ω0(μ-μ0)2+ω1(μ-μ1)2

(1)

式中ω0——C0出现的概率

ω1——C1出现的概率

μ0——C0的灰度均值

μ——图像总的灰度均值

μ1——C1的灰度均值

该算法以像元类间方差为分割标准，类间方差越大，说明两种物质的差别越大，因此可以将图像的背景与目标分割。设图像的任意灰度k，1≤k≤L(L为图像总灰度级)，将图像按灰度k分为C0和C1两类，σ2(k)取最大值的k即为最佳阈值。

1.5.2Canny边缘检测算法

无人机采集的高分辨率热红外影像冠层边缘具有明显的分界像元，因此可以对热红外图像进行边缘检测，得到棉花冠层边缘特征栅格图，将边缘特征栅格图进行重分类、线转面、裁剪等图像处理操作即可将冠层像元与土壤背景有效分离。边缘检测是图像预处理中的重要过程，其主要利用图像灰度阶跃变化的像素点集合，反映了图像DN值的突变性。传统的边缘检测算子有Sobel、Robert、Prewitt等[28]，但其都属于局部窗口的梯度算子，边缘检测效果并不理想[29]。Canny边缘检测算法具有边缘检验的3个判定准则[30]，因此对图像进行边缘提取具有较好的应用效果,并且Canny边缘检测算法在水平与垂直方向上的梯度变化检测具有较大优势[31]。Canny边缘检测算法借助高斯滤波器进行平滑处理，其高斯函数为

(2)

式中σ——高斯滤波器参数

应用一阶微分算子处理平滑后的影像I，并获得相应的梯度幅值G和梯度方向θ，点I(i，j)处x和y的偏导数分别为Gx(i，j)和Gy(i，j)[30]，即

(3)

(4)

其中

Gx(i,j)=(I(i,j+1)-I(i,j)+
I(i+1,j+1)-I(i+1,j))/2

(5)

Gy(i,j)=(I(i,j)-I(i+1,j)+I(i,j+1)-
I(i+1,j+1))/2

(6)

Canny边缘检测算法提取的边缘特征图像，首先将边缘特征图像在ArcGIS和ENVI中进行重分类和配准，再使用ArcToolbox工具中的栅格转线，线转面，拓扑检查及高级编辑功能将边缘线状图像进行闭合处理；最后将提取的棉花冠层面状矢量图导入ENVI软件，裁剪并使用建立掩膜工具和掩膜统计功能导出纯冠层像元的温度直方图和统计特征值，取本研究区具有代表性的局部，应用Canny边缘检测算法提取纯冠层像元的具体流程及操作见图4。

1.6 作物水分胁迫指数

冠层能量平衡公式[7]为

Rn=G+H+λE

(7)

其中

H=ρCp(Ta-Tc)/ra

(8)

(9)

式中Rn——太阳净辐射，W/m2

G——土壤热通量，W/m2

H——显热通量，W/m2

λE——潜热通量，W/m2

ρ——空气密度，kg/m3

Cp——空气热容量，J/(kg·K)

Ta——大气温度，℃

Tc——冠层温度，℃

ra——空气动力学阻力，s/m

ea——Ta时的水气压，Pa

rc——作物冠层阻力，s/m

γ——湿度计常数，Pa/℃

理论模式下的水分胁迫指数的计算公式为

(10)

(11)

式中rcp——潜在蒸散发下作物冠层阻力[32]，s/m

Δ——饱和水汽压与温度曲线的斜率，kPa/℃

2 结果与分析

2.1 剔除土壤背景图像结果分析

应用二值化Ostu算法和Canny边缘检测算法分别对研究区的热红外影像进行土壤背景剔除(图5)。

图5a为原始无人机热红外影像，图中包括4种水分处理(I1、I2、I3、I4)，每种处理3个重复，总共12个小区，并且4种水分处理的小区形成不同覆盖度的棉花冠层。图5b为二值化Ostu算法剔除土壤后的图像，从图中可以看出，二值化Ostu算法对于低覆盖度的棉花冠层(I11、I12、I13)剔除土壤背景效果较为理想，其中I11小区的红外影像由于光线太弱的缘故，导致较多的冠层像元被剔除；然而，该算法不能较好地区分高覆盖度作物冠层(I41、I42、I43)的棉花与土壤像元，而对于中等覆盖度的局部冠层分割效果也存在混合像元干扰问题，因此二值化Ostu算法的分割效果总体不理想。图5c为利用Canny边缘检测算法进行热红外图像边缘检测并利用本文的图像处理操作后的分割效果，从图中可得，Canny边缘检测算法可以精准地识别棉花冠层与土壤背景的边界像元，并能较好地提取棉花纯冠层像元，其分割效果优于二值化Ostu算法。

图5 剔除热红外图像土壤背景的方法对比Fig.5 Comparison of different methods for removing soil pixels from thermal infrared images

2.2 剔除土壤背景温度直方图分析

图6分别为图4d研究区(2017年7月13日13:00采集)的原始热红外图像、二值化Otsu算法剔除土壤背景、Canny边缘检测算法剔除土壤背景后的温度直方图。

图6 不同处理方法的温度直方图Fig.6 Temperature histograms for different treatment methods

从图6a可知，包含土壤背景的原始热红外图像温度直方图为双峰形状[24, 28]，研究表明，双峰形状的直方图代表两种不同的物质类型，第1个峰占据大多数的温度像元，其代表的物质类型为棉花冠层，第2个峰占据较少的像元，并且其像元温度高于第1个峰的冠层像元温度，因此代表的是土壤背景；本文的二值化Ostu算法与Canny边缘检测算法最终目的是将第2个峰消除，通过热红外图像温度直方图的双峰或单峰形状可以判断土壤背景是否剔除。图6b为二值化Otsu算法处理后的冠层温度直方图，它的直方图形状与原始温度直方图的前半部分相似，该温度直方图粗糙，噪声多，不符合自然界同种物质正态分布的规律。图6c为Canny边缘检测算法剔除土壤背景后的冠层温度直方图，可以看出使用Canny边缘检测算法剔除土壤背景后的直方图呈现明显的单峰形状，由于原始图像温度直方图冠层像元部分具有凸起形状，所以应用Canny边缘检测算法剔除土壤背景后的冠层温度直方图下侧有小的凸起峰；该直方图大致服从偏态分布，轮廓光滑平整，噪声较少，因此本文的Canny边缘检测算法相比二值化Ostu算法具有较大优势。

通过以上分析可知，Canny边缘检测算法剔除热红外影像土壤背景干扰的效果较为理想，并且纯冠层温度直方图为单峰形状，为了进一步验证Canny边缘检测算法的背景剔除效果，本文绘制12个试验小区(2017年7月13日13:00采集)剔除土壤背景后的冠层温度直方图(图7)。

图7 不同试验处理小区的冠层温度直方图Fig.7 Canopy temperature histograms of different treatment plots

图8 不同图像处理CWSI日平均变化趋势Fig.8 Daily average trends of CWSI by using different image processing methods

从图7可以看出，利用Canny边缘检测算法剔除12个试验小区土壤背景后的冠层温度直方图都为单峰形状，轮廓光滑，效果较好。通过对比分析图7不同水分处理小区的冠层温度直方图可得，随着水分胁迫程度的递减，试验小区冠层温度逐渐降低，并且冠层温度的变化范围区间也随之减小，因此，Canny边缘检测算法可以较好地剔除热红外土壤背景干扰，是提取无人机热红外纯冠层温度像元的可行方法。

2.3 CWSI日平均变化趋势

利用式(10)分别计算Canny边缘检测算法与Ostu算法剔除土壤前后的水分胁迫指数CWSI，并绘制不同试验处理小区日平均CWSI变化过程(图8)。

图8a是包含土壤背景的热红外图像计算的水分胁迫指数CWSI，从图8a可以看出，不同处理小区CWSI由大到小依次为：I1、I2、I3、I4，4种水分处理具有明显的差异性，并且随着灌水日期的推移土壤水分不断被作物消耗，CWSI逐渐上升。图8a的水分胁迫指数CWSI明显高于图8b和图8c的CWSI，因此无人机高分辨率热红外图像的土壤温度像元对CWSI计算具有极大的影响，图中对照组I4小区的水分胁迫指数大约为0.27，I1小区的水分胁迫指数高达0.55。对比分析可知，土壤温度像元极大地高估了水分胁迫指数CWSI，严重影响无人机热红外诊断作物水分状况的精度。

图8b为二值化Ostu算法剔除土壤背景后的CWSI日平均变化趋势。从CWSI的整体趋势来看，二值化Ostu算法的CWSI相对包含土壤的CWSI整体偏小，其最大值为0.42，最小值为0.22。图8c是利用Canny边缘检测算法剔除土壤背景后计算的水分胁迫指数CWSI，图中对照组I4小区CWSI的曲线从0.2增大至0.25，其CWSI变化范围较小，说明作物充分供水时水分胁迫指数相对偏小，I1为缺水处理小区，其CWSI最大值接近0.4，并小于包含土壤背景的水分胁迫指数CWSI。

Canny边缘检测算法与二值化Ostu算法剔除土壤背景后，相对于包含土壤背景的水分胁迫指数CWSI都相应的减小，说明高分辨率无人机热红外图像的土壤背景对水分胁迫指数具有较大影响。然而，Canny边缘检测算法小于二值化Ostu算法的CWSI，对比两种算法可知，二值化Ostu算法剔除土壤背景局部效果较差，而Canny边缘检测算法可以完全剔除土壤背景，对于无人机热红外诊断棉花水分状况具有重要的意义。

2.4 CWSI与棉花叶片气孔导度的关系

棉花叶片的气孔导度(Gs，mol/(m2·s))反映了棉花水分的真实情况，可以作为棉花是否缺水的判定标准。为了验证剔除土壤背景后水分胁迫指数CWSI的准确性，本文利用棉花叶片的气孔导度Gs与CWSI建立关系模型(表2)。

表2 CWSI与气孔导度的相关性Tab.2 Correlation between CWSI and stomatal conductance

注：表中Y为气孔导度Gs，X为水分胁迫指数CWSI。下同。

从表2可得，水分胁迫指数CWSI与叶片气孔导度Gs具有极高的线性关系，在13:00时，包含土壤背景水分胁迫指数CWSI与气孔导度Gs的模型决定系数R2为0.697，达到显著相关，但09:00与17:00的模型相关性较低。二值化Ostu算法剔除土壤背景后的CWSI与气孔导度Gs相关关系并没有提高，可能由于Ostu算法对于低覆盖度的棉花冠层剔除土壤效果理想，而高覆盖度的棉花热红外冠层剔除土壤效果较差，因此最终计算的水分胁迫指数偏差更大，导致CWSI与气孔导度Gs的关系较差。Canny边缘检测算法剔除土壤后的CWSI与气孔导度Gs模型关系在13:00时的决定系数高达0.84，并且09:00与17:00的模型相关性也高于包含土壤背景处理。对比可知，土壤背景温度对水分胁迫的计算具有较大的干扰，应用Canny边缘检测算法剔除无人机高分辨率热红外影像的土壤背景可以有效提高水分胁迫指数CWSI的计算精度。

2.5 CWSI与棉花叶片气孔导度的模型验证

利用棉花叶片实测气孔导度Gs数据，验证表2中的关系模型，均方根误差RMSE和决定系数R2共同作为模型检验精度的评判标准(表3)。

从表3可知，3种热红外图像处理方法的模型验证精度均在13:00最高，并且包含土壤背景、二值化Ostu算法、Canny边缘检测算法的水分胁迫指数CWSI与棉花叶片气孔导度Gs验证精度在13:00的决定系数R2分别为0.734、0.775、0.893，均方根误差RMSE分别为0.086、0.065、0.041，可得Canny边缘检测算法的模型验证精度最高，二值化Ostu算法次之，包含土壤背景的验证精度最差；在09:00与17:00的模型验证精度，Canny边缘检测算法同样高于包含土壤背景和二值化Ostu算法的模型验证精度；从模型验证的决定系数R2与均方根误差RMSE来看，Canny边缘检测算法提高了CWSI计算精度，可以更准确地诊断棉花水分状况。

表3 CWSI与气孔导度关系的验证Tab.3 Verification of relationship between CWSI and stomatal conductance

3 讨论

本文选取二值化Ostu算法和Canny边缘检测算法分别对无人机热红外影像的土壤背景进行剔除。二值化Ostu作为一种优化的算法，其可以直接对图像进行阈值分割，而Canny边缘检测算法可以对图像的边缘特征栅格图像处理后间接实现图像分割处理，并且Canny边缘检测算法对无人机热红外图像的阈值分割鲜有报道。研究表明，Canny边缘检测算法剔除无人机热红外图像的土壤背景效果较为理想，而二值化Ostu算法只适用于低覆盖度的棉花冠层。目前剔除土壤背景提取纯冠层温度像元的主要方法有：RGB图像阈值法[33]、热红外图像直方图阈值法(Tair-10℃

从无人机监测作物水分胁迫状况的最佳时间分析，本文研究了无人机监测棉花水分胁迫状况的3个时刻。通过对比分析可知，包含土壤背景、二值化Ostu算法和Canny边缘检测算法3种处理方法13:00的模型精度最高。并且，在13:00利用Canny边缘检测算法处理得到的水分胁迫指数CWSI与棉花叶片气孔导度Gs的决定系数R2高达0.84，无人机热红外监测棉花水分胁迫精度最高。棉花叶片气孔导度Gs是反映作物叶片蒸腾强度和水分消耗规律的指标，随着日出后太阳净辐射的增强，大气温度回升，冠层温度趋近于大气温度[36]，作物水分消耗速率和气孔导度Gs逐渐增大，土壤可提供作物根系一定量的水分，以保证作物的生理活动。到13:00，如果土壤不能供给作物根系充足的水分，作物则表现出水分胁迫状态，气孔导度Gs变小，反之亦然。而此时，作物水分胁迫现象最为明显，不同水分处理小区水分胁迫差异性最大，是作物水分胁迫诊断的最佳时刻，该时刻大气温度较高，作物耗水多，比较符合作物的生理耗水规律，这为无人机监测作物水分状况提供了实践应用基础。因此，13:00的作物水分胁迫指数CWSI与叶片气孔导度Gs的相关关系高于其他时刻。

然而，本研究也存在不足，棉花冠层理论上可以分为阳光直接照射叶片和阴影叶片，并且两种冠层叶片的温度差别较大，将两种冠层叶片分别进行分析，可以更深入地了解棉花不同冠层温度对水分响应的变化规律，并能进一步提高无人机热红外诊断作物水分的精度。但是，从热红外图像处理角度分析，区分两种冠层叶片成分具有较大的难度，需要研究更优化可行的图像分割算法；从无人机热红外技术角度分析，热红外图像空间分辨率较低仍是瓶颈问题。因此，无人机热红外诊断作物水分状况所面临的难点需要更深入的研究和探讨。

4 结论

(1)Canny边缘检测算法剔除高分辨率热红外影像的土壤背景优于二值化Ostu算法，对于不同覆盖度的作物冠层，Canny边缘检测算法具有较大的适用性和理想的分割效果；原始热红外影像的温度直方图为双峰形状，其主峰为作物冠层，占据大部分像元，偏锋为土壤背景像元，剔除土壤背景后的热红外图像冠层温度直方图呈单峰形状的偏态分布。

(2)Canny边缘检测算法剔除土壤背景后的水分胁迫指数CWSI最小，二值化Ostu算法次之，包含土壤背景的最大，土壤背景温度像元对水分胁迫指数CWSI的影响较大，因此应用无人机热红外技术诊断作物水分状况应该考虑土壤背景的影响。

(3)Canny边缘检测算法剔除土壤背景的水分胁迫指数CWSI与棉花叶片气孔导度Gs的模型相关关系最高，优于二值化Ostu算法和包含土壤背景的模型精度；通过对无人机热红外3个时刻监测的水分胁迫指数CWSI分析可知，Canny边缘检测算法、二值化Ostu算法、包含土壤背景3种处理均在13:00的模型精度最高，09:00与17:00的模型精度相近，因此本研究的最佳无人机诊断作物水分状况的时间为13:00，并且Canny边缘检测算法明显提高棉花水分胁迫诊断的精度。