张智韬,于广多,吴天奎,张誉馨,白旭乾,杨 帅,周永财
(1. 西北农林科技大学水利与建筑工程学院,杨凌 712100;2. 西北农林科技大学旱区农业水土工程教育部重点实验室,杨凌 712100)
及时准确地获取作物水分胁迫状况对发展节水灌溉农业、提高水分利用效率有重要意义[1]。冠层温度由土壤—植物—大气连续体内热量的传输与转换决定[2],冠层温度及其相关指数可以指示作物水分胁迫状况[3-4],其中以作物水分胁迫指数(Crop Water Stress Index,CWSI)应用最为广泛。Idso 等[5]提出基于冠层温度的CWSI 经验公式,而后又有学者提出了CWSI 的理论公式与简化公式[6-8],通过蒸散模型计算CWSI 的理论公式考虑了诸多气象因素[9],理论基础较强,监测精确度较高。近年来,国内外低空无人机遥感技术快速发展,可借助热红外成像仪快速、大面积地采集冠层温度信息,为监测作物水分提供了一种可靠路径[10-12]。然而当作物未达到全覆盖时,无人机获取的试验区域图像中包含冠层、土壤等多种地物信息,并且不同地物的温度在一定范围内存在交集,这对提取冠层温度有较大的干扰,因此剔除无人机影像中的土壤背景可以提高作物水分胁迫的监测精度[13-14],快速准确地对热红外图像进行地物分类并提取其中的冠层温度是十分必要的。
按照是否借助可见光图像可将冠层温度提取方法分为两类:一种是依据冠层与土壤的温度信息差异,直接在热红外图像上分类,如经验公式法[15]、边缘检测法等[16-17],统称为直接法。经验公式法十分依赖所测定的大气温度的准确性,误差较大[18];边缘检测法要求图像具有较高的分辨率、不同地物温度差异较大且地物类别不宜过多[19],这些因素在一定程度上限制了此方法的应用。另一种是借助可见光图像的光谱信息实现热红外图像地物分离的间接法[20-22]。虽然高分辨率的可见光相机能够很好地区分出复杂地物信息,但在实际操作中,受相机自身参数、拍摄角度、气象因素以及人为因素等多个因素的干扰,两图像像元无法准确无误地一一对应[23],这直接影响冠层温度的提取精度。直接法与间接法各有优缺点,而现阶段对两类方法的分类效果进行比较,并将两类方法相结合来协同提取冠层温度的研究较少。
鉴于此,本研究以4 种不同水分处理的拔节期夏玉米为研究对象,利用无人机获取试验区域每日13:00 的热红外和可见光图像资料,分别采用 Otsu 算法、EXG-Kmeans 算法和Otsu-EXG-Kmeans 算法获取冠层区域图像,并从冠层轮廓图、冠层区域提取精度、温度直方图以及与实测冠层温度相关情况4 个方面对提取结果进行精度评价,以确定提取热红外图像冠层温度的最佳方法,而后求得相应CWSI,分析CWSI 同土壤含水率相关关系以及CWSI 日平均变化趋势,对玉米水分亏缺状况进行定性分析,以期为无人机遥感精准监测作物水分胁迫状况提供参考。
研究区域位于陕西省杨凌示范区的中国旱区节水农业研究院,34°17′N,108°04′E,海拔525 m。该区域属于温带大陆性季风气候,多年平均降水量约640 mm,主要集中在7-9 月,试验期间大气平均温度为33 ℃,总辐射较强。试验区域的土壤类型为中壤土,凋萎系数8.6%,0~1 m 平均田间持水量为23%~26%(均为质量含水率),平均土壤干容重1.44 g/cm3,地下水埋藏较深,其向上的补给量忽略不计。
试验设置严重水分胁迫(T1)、中度水分胁迫(T2)、轻度水分胁迫(T3)以及充分灌溉(T4)4 种水分处理梯度,土壤含水率分别维持在田间持水量的50%、65%、80%以及95%~100%[17],每个水分处理设有3 个重复,共计12 个小区(图1),小区为4 m×4 m,2 个小区之间布置保护行。玉米品种为“华农138”,播种时间为2020年6 月14 日,出苗时间为6 月20 日,播种深度约5 cm,种植行距60 cm,株距30 cm,行沿南北走向。试验区域设置小型气象站,并配有移动折叠雨棚以防止降雨干扰。
拔节期玉米水分胁迫程度直接影响玉米长势,且无人机相机在晴朗无风天气拍摄效果较好。考虑以上因素,本研究选择在2020 年7 月26 日、7 月27 日、7 月29 日、7 月31 日、8 月2 日的13:00 进行无人机图像以及地面数据的采集。
1)热红外与可见光图像采集
使用大疆公司研发的Matrice600Pro 六旋翼无人机,并搭载ZENMUSE XT 热红外相机与ZENMUSE X3 可见光相机进行图像采集。设置无人机飞行高度为20 m,镜头垂直地面进行拍摄,热红外相机工作波段为7.5~13.5μm,像素为640×512,水平视场角32°,垂直视场角26°,镜头焦距19 mm,温度分辨率0.05 ℃,20 m 图像分辨率为1.8 cm。可见光相机像素为4 000×3 000,视场角94°,镜头焦距20 mm,20 m 图像分辨率为1.0 cm。
2)地面数据采集
在各小区内选取可以代表小区平均长势且生长良好的3 株玉米植株进行标记[24],在图像采集完成后,及时利用测温枪获取标记植株承光面的冠层温度,同时利用土钻取土,取土深度为20 cm[25],并通过烘干法测定土壤含水率。
1.4.1 无人机图像预处理
首先采用PIX4D 软件对可见光与热红外图像进行拼接处理,而后以热红外图像为基准,利用ENVI Classic软件手动选取热红外与可见光图像中明显的地物特征点进行配准(特征点选取30 个以上,且保证每个小区包含至少2 个特征点,特征点均方根误差小于1),最后利用FLIR Tools 软件将热红外图像的灰度值转换为温度,并用实测水温校准温度。
1.4.2 热红外图像冠层温度提取
1)Otsu 算法
热红外相机获取的图像可分为两类:1)冠层像元,2)土壤像元[13]。Otsu 算法是灰度图像阈值分割的经典自适应阈值算法,该算法以像元类间方差为分割标准,将两种像元类间方差最大时的阈值作为冠层与土壤温度的分界线,其核心公式如下[26]:
式中σ2为类间方差;k为阈值;σ2(k)取最大值时的k值即为最佳阈值;w0为冠层像元出现的概率;w1为土壤像元出现的概率;μ为图像灰度均值;μ0为冠层灰度均值;μ1为土壤灰度均值。
2)EXG-Kmeans 算法
可见光图像的空间分辨率较高,且玉米冠层对绿波段反射强烈而土壤对绿波段反射较弱[27],故本研究采用EXG 指数[28]区分冠层与土壤信息,EXG 指数计算式如下:
式中R、G、B分别表示可见光图像红、绿、蓝波段灰度值。
同时为避免分类结果受主观因素的影响,采取非监督分类中的Kmeans 算法对EXG 指数图像进行分类,生成包含玉米冠层像素轮廓特征的矢量文件,而后对热红外图像进行掩膜处理获取冠层温度信息。
3)Otsu-EXG-Kmeans 算法
研究表明Otsu 算法可以有效剔除图像中阳光直射土壤,而对阴影覆盖土壤剔除效果较差[29];借助可见光图像可以准确确定冠层区域,对阴影土壤与阳光直射土壤的剔除效果较好,但由于存在配准偏差,导致少量土壤区域被误分到冠层区域中[23],影响EXG-Kmeans 算法提取冠层温度的精确度。因此本研究首先利用Otsu 算法剔除图像中阳光直射土壤,而后利用EXG-Kmeans算法剔除阴影土壤,称为Otsu-EXG-Kmeans 算法,冠层温度提取过程通过MATLAB 与ENVI 实现,计算流程如图2 所示。
采用Jones 等[7]提出并由张仁华[9]改良后的理论模型计算CWSI,公式如下:
式中Tc为图像冠层温度,℃;Ta为大气温度,℃;γ为湿度计常数,kPa/℃;Δ为饱和水汽压与温度曲线斜率,kPa/℃;d为空气饱和差,kPa;Rn为太阳净辐射,W/m3;G为土壤热通量,W/m3;ra为空气动力阻力,s/m;ρ为空气密度,kg/m3;Cp为空气比热,J/(kg·℃)。
本研究采用生产者精度(Producer's Accuracy,PA)和用户精度(User's Accuracy,UA)对热红外图像冠层区域提取结果进行量化评价,两者越接近100%,说明冠层区域提取精度越高,具体公式如下:
式中TP是预测为冠层的冠层样本个数;FP是预测为冠层的土壤样本个数;FN是预测为土壤的冠层样本个数。
本研究通过相关系数(r)评价两者的相关关系,r越接近1 两者相关性越高;通过均方根误差(Root Mean Square Error,RMSE)评价预测值的误差,RMSE 越接近0,预测值的误差越小。
2.1.1 不同算法冠层区域提取结果分析
以8 月2 日T1-1 小区中心区域为例,分别在RGB 影像与热红外影像中绘制不同算法确定的冠层区域图(图3)。对比RGB 影像可知,Otsu 算法确定的冠层区域包含极大部分玉米像元以及部分温度与冠层温度相近的土壤像元,而EXG-Kmeans 算法可在RGB 影像中准确剔除阴影覆盖土壤与阳光直射土壤。由于RGB 影像与热红外影像在配准时存在细微偏差,导致少量高于冠层温度的土壤像元被误划于冠层区域,对比热红外影像可知,EXG-Kmeans 算法无法避免此类现象发生,而Otsu-EXG-Kmeans 算法由于加入温度阈值条件,可更好地剔除高于冠层温度的土壤像元。
2.1.2 不同水分处理冠层区域提取结果分析
分别运用 Otsu 算法、 EXG-Kmeans 算法和Otsu-EXG-Kmeans 算法对热红外图像进行分类处理,每种方法的生产者精度和用户精度如表1 所示。
表1 不同提取方法分类精度Table 1 Classification accuracy of different extraction methods
依据图3 所示冠层区域并对比表1 中同一算法的生产者精度与用户精度可知:由于Otsu 算法确定的冠层区域包含大量土壤样本,而土壤区域只有极少数的冠层样本, 此算法的生产者精度远远大于用户精度(96.1%>66.6%);利用EXG-Kmeans 算法进行地物分类,误分到冠层区域的土壤样本与误分到土壤区域的冠层样本数量相近,生产者精度与用户精度相差较小;Otsu-EXG-Kmeans 算法可剔除误分到冠层区域的土壤样本,而无法将误分到冠层区域的土壤样本召回,故用户精度大于生产者精度(95.9%>85.6%)。
对比Otsu 算法对不同水分处理的分类效果,其中T1处理的用户精度高于T4 处理而生产者精度低于T4 处理(68.4%>64.0%,94.5%<97.2%),由此推断Otsu 算法在严重水分胁迫处理区域(T1)的适用性优于充分灌溉处理(T4),但是误分到土壤区域的冠层像元数量也要大于充分灌溉处理。EXG-Kmeans 算法中T4 处理的用户精度与生产者精度均高于T1 处理(89.4%>77.9%,88.7%>79.1%),故此算法对充分灌溉处理冠层区域提取效果要优于严重水分胁迫处理。而Otsu-EXG-Kmeans 算法的用户精度和生产者精度同水分胁迫情况呈负相关,且用户精度均值达95.9%,对冠层提取效果优于其他2 种算法。
水分胁迫处理越严重的小区冠层温度越高,温度极差越大[30],不同算法温度统计结果越明显。以8 月2 日T1-1 小区为例,绘制原始图像与不同分类方法获得的冠层温度频率直方图(图4),并统计冠层温度的最小值、最大值以及平均值(表2)。
表2 图像冠层温度统计Table 2 Image canopy temperature statistics ℃
地面实测数据表明,阳光直射土壤温度大于冠层温度,而阴影覆盖土壤温度与冠层温度接近,甚至低于冠层温度,故作物半覆盖区域的温度频率直方图大致为双峰形状。如图4a 所示,原始温度频率直方图左侧峰频率较大且温度较低,代表冠层集中区域,右侧峰频率较小且温度较高,代表土壤集中区域。原始图像与EXG-Kmeans 算法相比于Otsu 算法与Otsu-EXG-Kmeans算法,提取的冠层温度分布范围更广,进一步说明EXG-Kmeans 算法由于存在配准偏差的现象,少量高温的光照土壤像元无法有效剔除。由表 2 可知,Otsu-EXG-Kmeans 算法取得的冠层温度最小值与EXG-Kmeans 算法同为31.0℃,最大值与Otsu 算法同为46.4℃,故此方法有效地剔除了像元温度低于冠层温度的阴影覆盖土壤以及像元温度高于冠层温度的阳光直射土壤。对比冠层温度平均值,3 种方法均可有效降低土壤温度的干扰,且Otsu-EXG-Kmeans 算法效果最佳,Otsu 算法次之,EXG-Kmeans 算法最差(冠层温度:38.2℃<38.9℃<39.9℃<42.2℃)。
为进一步评估不同方法的提取效果,图像提取温度与地面实测温度进行相关性分析(n=60),如图5 所示,随着土壤像元的减少,图像冠层温度逐渐降低。对比3种处理方法,Otsu-EXG-Kmeans 算法提取的冠层温度与地面实测温度相关性最佳(r=0.788),误差最小(RMSE=1.901 ℃),同时Otsu-EXG-Kmeans 算法的拟合直线相比较其他方法,截距更小且斜率(回归系数)更接近1,整体趋势更接近1:1 线,故利用Otsu-EXG-Kmeans算法得到的冠层温度精确度更高。
2.4.1 冠层温度指标与土壤含水率相关性分析
上述研究表明,Otsu-EXG-Kmeans 算法对不同水分处理的适用性均较强,得到的冠层温度精确度也较高,故可借助Otsu-EXG-Kmeans 算法监测玉米水分胁迫状况。作物通过根系吸收土壤水分来进行蒸腾作用,从而维持冠层温度的平衡,因此土壤含水率与玉米缺水状况呈负相关[31],是监测玉米缺水状况的有效指标。Otsu-EXG-Kmeans 算法提取的冠层温度值代入式(3)计算出相应的作物水分胁迫指数,将冠层温度与作物水分胁迫指数分别同0~20 cm 土壤含水率进行相关性分析(n=60),以间接监测玉米缺水状况。结果表明(图6),冠层温度与作物水分胁迫指数均与土壤含水率负相关,即同玉米缺水状况呈正相关,相比于冠层温度,作物水分胁迫指数可以更为准确地指示作物的水分亏缺状况(r= -0.738)。
2.4.2 冠层温度与CWSI 日平均变化趋势对比
分别取各水分梯度3 个重复小区的平均值作为本梯度的冠层温度值与作物水分胁迫指数,由此绘制不同水分处理的冠层温度与作物水分胁迫指数日平均变化趋势线如图7 所示,不同水分处理的冠层温度与作物水分胁迫指数具有明显差异,并且T1>T2>T3>T4。根据试验记录显示,7 月28 日对T3、T4 处理进行了灌水处理,故2种处理在29 日出现拐点,而T1 与T2 梯度在试验期间一直未进行灌水处理,随着时间的推移,土壤水分不断被作物消耗,玉米受水分胁迫程度逐渐增加,相比于冠层温度,作物水分胁迫指数日平均变化趋势更好地描述了这一现象。
2.4.3 冠层温度与CWSI 分布图
为更直观地展示试验区域不同处理的作物水分胁迫程度,以8 月2 日为例,绘制Otsu-EXG-Kmeans 算法获取的冠层温度与作物水分胁迫指数分布图,如图8 所示,整体冠层温度在28~46.5 ℃之间分布,作物水分胁迫指数在0~1 之间分布。
本研究结果表明,Otsu 算法对水分亏缺严重的区域分类效果较好,对充分灌溉区域分类效果较差,然而,EXG-Kmeans 算法对不同水分胁迫处理的分类效果与Otsu 算法截然相反。通过探究可知,在供水充足的状态下,植物通过蒸腾作用散失大量辐射热,从而降低叶温和植物体温,避免高温伤害,同样地表也可通过蒸发作用延缓土壤温度上升,与此同时,作物体内的水分为叶片光合作用提供原料,叶绿体可正常工作,热红外图像中的冠层与土壤温度相差较小,而EXG-Kmeans 算法主要依据可见光图像中冠层与土壤像元色彩差异对图像进行分类,充分灌溉条件下,可见光图像中的冠层像元为绿色,与棕黄色土壤色差较大,配准法较为实用。作物水分胁迫区域,作物蒸腾作用与土壤蒸发作用减弱,温度均升高,而由于土壤的比热容相对水较小,升温比叶片速率快,同时缺水会导致叶片气孔导度下降,进入叶内的二氧化碳减少,光合产物在叶片中积累对光合作用产生反馈抑制作用,也会影响细胞伸长并抑制蛋白质的合成,叶片易出现暂时性萎蔫,使光合面积降低[32],严重缺水时甚至造成叶绿体类囊体结构被破坏,叶片对红光与蓝光吸收能力下降[33],热红外图像中的冠层与土壤温度相差较大,可见光图像中的冠层像元偏向黄色,与棕黄色土壤色差较小,而Otsu 算法主要依据热红外图像中冠层与土壤温度差异对图像进行分类,故在水分亏缺时Otsu 算法分类效果较好。上述情况是水分胁迫程度与冠层温度呈正相关的主要原因, 也是Otsu-EXG-Kmeans 算法在不同水分胁迫处理中适用性均较强的主要原因。
通过冠层温度提取值与实测值相关性分析发现,水分胁迫严重的小区,部分冠层温度提取值大于实测值,而充分灌溉的小区,部分冠层温度提取值小于实测值。这一现象出现的原因是实测值选取的是生长良好的承光面冠层区域,然而水分胁迫过重的区域玉米长势相对较差,且容易发生叶片卷曲现象(T1 处理中心区域),在提取过程中这部分温度较高的区域参与计算,导致冠层温度提取值高于实测值。反之充分灌溉区域玉米长势较好,冠层覆盖度较高,部分阴影投映在冠层区域,在提取过程中温度较低的阴影冠层参与计算,导致冠层温度提取值低于实测值。
本研究方法也存在不足之处:首先,在水分严重胁迫区域,部分玉米叶片变黄,其冠层温度与土壤温度接近,本研究提出的方法无法准确确定此部分冠层区域;其次,通过上述分析可知,由于剔除了可见光影像误分到冠层中的土壤样本,本研究方法确定的冠层区域小于仅借助可见光影像确定的冠层区域,即采用本研究方法计算的作物覆盖度略低;此外,不同作物的耐旱性不同导致温度分布不同,不同地区下垫面不同导致光谱信息不同,因此本研究提出的方法是否适用于其他地区的水分胁迫监测仍需进一步研究与验证。
本研究通过不同方法提取玉米冠层温度,计算了基于冠层温度的作物水分胁迫指数,探究了其在水分胁迫监测中的精确度, 得出以下结论: 1 ) 利用Otsu-EXG-Kmeans 算法获取的冠层区域图像最优,EXG-Kmeans 算法次之,Otsu 算法最差(用户精度:95.9%>83.2%>66.6%);Otsu-EXG-Kmeans 算法提取的冠层温度与实测温度相关关系接近1:1 线,相关性相比于其他方法更好(r=0.788),因此Otsu-EXG-Kmeans 算法是准确获取图像冠层温度的有效方法。2)相比于冠层温度,作物水分胁迫指数与土壤含水率的相关性更高,作物水分胁迫指数日平均变化趋势更符合实际情况,可更加精确地监测玉米缺水状况。