基于无人机遥感的青贮夏玉米水分亏缺指数反演研究

李星恕 程双飞 薛 志 熊秀芳 韩文霆,2 张立元

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

0 引言

我国是严重缺水国家，农业用水超过全国用水量的60%。玉米作为我国主要的粮食作物，其种植面积达到了0.2亿hm2[1]。在干旱半干旱地区，玉米种植主要依靠灌溉，而传统灌溉方式水资源利用率低、浪费严重。为了降低玉米生产过程中的灌溉用水，提高水资源利用效率，实时、准确监测大田玉米旱情非常必要。传统旱情监测方法一般通过地面传感器监测土壤含水率和植被生理参数，效率低、成本较高，监测较大面积的玉米旱情信息非常困难[2-4]。遥感技术提供了一种实时、无损的广域范围旱情监测方法，但是传统的卫星遥感时空精度较低，无法满足大田尺度的旱情信息监测[5]。而无人机遥感由于方便、快捷、精度高、成本低等优点，弥补了卫星遥感的不足，近年来在农业领域广泛应用[6-8]。

监测旱情的指标很多，其中冠层温度和以冠层温度为基础建立的多种干旱指数应用较多[9-12]，其中水分胁迫指数(Crop water stress index，CWSI)应用最为广泛。但CWSI干旱指数的计算需要精确获取作物的冠层温度[13-15]。遥感监测时，由于热红外相机分辨率较低，作物覆盖度较低的遥感影像中会存在大量的土壤-植被混合像元，从而影响作物冠层温度的提取[16]。虽然有学者采用基于RGB滤波等预处理方法消除混合像元对冠层温度提取的影响，但是这些方法需要主观选择阈值，存在耗时长、稳定性差等缺点，在实际应用中难以获得满意的水分胁迫监测效果[17]。因此，有学者在CWSI的基础上进一步提出了水分亏缺指数(Water deficit index，WDI)的概念[18]。WDI根据以植物的陆气温差(陆地表面混合温度与空气温度之差)和植被指数为横、纵坐标而形成的植被指数-温度梯形空间(Vegetation index-temperature trapezoid,VIT)的方法计算得到[19]；通过定义该梯形空间的4个边界点，在理论上将植物所有干旱情况包含在梯形空间内。与CWSI相比，WDI不需确定作物冠层温度，因此在植被覆盖度低的区域仍有较好的反演效果。近年来，有学者利用WDI监测矮化绿豆，发现WDI可监测产量和季节性作物水分变化状况[20]；利用卫星遥感反演WDI结合潜在蒸散量计算实际蒸散量，发现WDI可用于辅助计算实际蒸散量[21]；利用无人机遥感反演WDI研究水分胁迫对不同基因型苹果树的影响，结果表明，WDI作为干旱指标与果树茎水势密切相关，适用于不同覆盖度和冠层结构的苹果树旱情监测，可为苹果育种提供一种新的实时监测方法[22]；对不同生育期的大麦进行无人机遥感WDI反演，表明能够准确监测不同生育期的大麦旱情变化，且能用于及时发现大麦早熟[23]。近年来虽然利用无人机遥感、采用WDI指数对多种作物进行了干旱胁迫研究，但是利用无人机遥感技术实时监测拔节期青贮夏玉米旱情，研究不同时间尺度下不同水分胁迫对WDI干旱指数监测效果的影响，迄今为止研究较少。

为了研究水分胁迫对不同时间尺度下拔节期夏玉米WDI和陆气温差监测效果的影响，本文利用无人机遥感数据和地面采集数据，建立植被指数-温度梯形空间，计算WDI干旱指数并生成大田尺度的青贮夏玉米旱情分布图，对基于WDI与陆气温差的不同时间尺度下青贮夏玉米气孔导度和土壤含水率的相关性进行分析，并研究水分胁迫对其影响，为大田作物干旱信息的实时监测提供参考。

1 实验材料与方法

1.1 实验区概况

实验田位于内蒙古自治区鄂尔多斯市达拉特旗昭君镇(40°26′0.29″ N，109°36′25.99″ E)，海拔1 010 m，面积1.13 hm2。实验田土壤为砂壤土，田间持水量16%。种植作物为青贮夏玉米(品种为钧凯918)；播种时间为2018年5月8日，收获时间为9月8日，全生育期共120 d，拔节期为6月26日—7月26日。播种深度5 cm，行距58 cm，株距25 cm，按照当地种植习惯进行田间管理。实验田采用自走中心支轴式喷灌机灌溉。

实验田划分出3个扇形灌溉区域。为了研究不同水分胁迫程度对青贮夏玉米旱情监测的影响，控制喷灌机灌水量使不同灌溉区域处于100%充分灌溉(田间持水量的95%，1区)、65%充分灌溉(2区)、40%充分灌溉(3区)，如图1所示。为了方便采集地面数据，在距离实验田中心点25 m的圆弧上均匀选择3个采样区、距离中心点47.5 m的圆弧上均匀选择6个采样区(图1中矩形)，每个采样区面积6 m×6 m；这样每个灌溉区内有3个均匀分布的采样区。为了提高遥感影像精度，采用 RTK x5plus测量系统(科利达，中国)在实验田内均匀布设5个像控点用于后期遥感影像拼接过程中的图像校正(图1中圆形)。在每个矩形采样区选择对角线中点和端点3个点作为玉米生理参数的测量点。

图1 实验地样区概况Fig.1 Experimental area

1.2 数据采集

数据采集主要包括地面数据和无人机遥感数据获取。地面数据采集用于建立WDI梯形空间边界。无人机获取遥感数据提取青贮夏玉米的NDVI和陆地表面混合温度用于计算WDI。

1.2.1无人机数据获取

采用自主研发的无人机多光谱影像采集系统和无人机热红外影像采集系统分别获取青贮夏玉米的多光谱和热红外影像。两套系统及相机具体参数见表1和表2。

表1 无人机多光谱影像采集系统主要参数Tab.1 Main parameters of UAV multispectral image acquisition system

表2 无人机热红外影像采集系统主要参数Tab.2 Main parameters of UAV thermal infrared image acquisition system

在青贮夏玉米拔节期间选择晴朗无风的天气进行无人机遥感监测。飞行时间为12：00—13:00，采用固定航线飞行，飞行高度70 m，航向和旁向重叠度为85%，多光谱影像分辨率为4.6 cm，热红外影像分辨率为7.8 cm。在后期遥感影像处理过程中，以遥感影像中的像控点为基准拼接多光谱和热红外影像，多光谱影像利用黑色和白色漫反射板进行反射率校正；热红外影像利用雷泰ST60+型手持测温仪(RAYTEK，美国)测量黑色漫反射板、白色漫反射板和水体温度对热红外影像进行温度校准；再利用重采样的方法将两类影像像素统一降至8 cm；然后从多光谱影像中提取归一化植被指数(Normalized difference vegetation index，NDVI)，从热红外影像中提取陆地表面混合温度。归一化植被指数计算式为

NDVI=(IR-R)/(IR+R)

式中IR——近红外波段R——红光波段

1.2.2地面数据获取

为了获取建立WDI模型和表征干旱的土壤含水率(Soil water content，SWC)、气孔导度(Stomatal conductance，Gs)等数据，需要采集以下地面数据：

(1)气象站数据：为了减小气象站点和实验田的气象数据差异，根据FAO-56建立标准的下垫面和气象站[24]，从而能够利用一个气象站监测更广范围内的气象数据。农业气象站(中国河北清易电子科技有限公司，中国)位于实验区南1 km处，下垫面为紫花苜蓿。可采集太阳净辐射、空气温度、相对湿度、风速数据，采集频率30 min，选择12:00—14:00的气象数据计算NDVI-陆气温差梯形空间边界。

(2)土壤含水率：采用干燥法测量实验田土壤含水率。在采样区中心点附近采集10 cm深度的土样，将干燥箱温度设置105℃干燥8 h至恒质量。

(3)气孔导度：在采样区的3个测量点附近选择有代表性的玉米植株，用AP4型气孔计(Delta-T公司，英国)测量被太阳直射的玉米第1片全展叶的叶片中间部分的气孔导度；为了确定凋萎状态下和充分灌溉下的玉米气孔导度，在玉米生育期内根据种植经验选择实验田接近凋萎状态和充分灌溉的玉米植株各10株，测量其气孔导度并取平均值。选择最小值作为凋萎状态的玉米气孔导度，选择最大值作为充分灌溉状态的玉米气孔导度。

(4)株高：在采样区3个测量点附近随机选择3株玉米测量株高，取平均值。

(5)叶面积指数：在采样区内侧，选择距离采样区边界50 cm的玉米植株，用LAI-2200C型(LICOR，美国)冠层分析仪测量叶面积指数。

(6)温度：在采样区3个测量点附近，利用雷泰ST60+型手持测温仪在距离地表20 cm处垂直测量土壤表面温度，在玉米株高80%高度处沿水平方向测量冠层温度，取平均值。

1.3 基于无人机遥感的青贮夏玉米WDI反演方法

采用NDVI-陆气温差建立WDI梯形空间，WDI反演技术路线见图2。通过定义梯形空间的4个边界点：1点为水分充足的郁闭冠层、2点为凋萎状态下的郁闭冠层、3点为水分充足的裸土、4点为严重水分胁迫的裸土，再结合地面数据计算和确定青贮夏玉米的WDI梯形空间边界；线段13为WDI模型的湿边，线段24为WDI模型干边。通过多光谱影像中提取的NDVI和热红外影像中提取的陆地表面混合温度(与空气温度之差为陆气温差)，确定青贮夏玉米在梯形空间中的具体位置，见图3中B点。

图2 WDI反演技术路线图Fig.2 Technology road map for inversion WDI

图3 梯形空间示意图Fig.3 VIT trapezoid space

水分亏缺指数WDI计算公式为[15]

(1)

式中Ts——陆地表面混合温度，℃

Ta——空气温度，℃

WDI梯形空间的建立是以能量平衡为基础[19]。作物陆气温差的计算公式为[14]

Ts-Ta=

(2)

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

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

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

rc——冠层阻力，s/m

Cv——空气体积热容，J/(m3·℃)

Δ——饱和蒸气压与空气温度的斜率，kPa/℃

γ——计量常数，Pa/K

VPD——饱和蒸气压，kPa

其中Δ、γ、VPD参考FAO-56得到[24]。

利用式(2)计算WDI梯形空间4个顶点时，1、2点分别表示水分充足和凋萎状态下的郁闭冠层,因此冠层阻力rc分别用潜在蒸散量状态下的冠层阻力rcp和凋萎状态下的冠层阻力rex代替；3、4点分别为水分充足和水分胁迫时的裸土，因此可以定义3点和4点的冠层阻力rc分别为零和无穷大。裸土时(3、4点)，热量传递过程主要受到空气阻力ra和土壤阻力rs的影响，因此计算陆气温差时将ra替换为ra+rs。

确定梯形空间所需的参数计算公式如下：

冠层阻力rc[25]

(3)

式中Gs——叶面气孔导度，mol/(m2·s)

LAI——叶面积指数

空气动力学阻力ra[26]

(4)

其中

zo≈0.13hcd0=0.65hc

式中z——风速计高度，m

hc——植株最大高度，m

u——风速，m/s

zo——表面粗糙度，m

d0——零点高度，m

土壤阻力rs[27]

(5)

式中rs——土壤阻力，s/m

Tsoil——裸土表面温度，℃

Tc——冠层温度，℃b、c——常数

us——土壤表面上方0.2 m处的风速，m/s

其中c=0.002 5,b=0.012,us参考FAO-56利用2 m处风速计算得到[24]。

2 结果与分析

2.1 青贮夏玉米拔节期WDI梯形空间

青贮夏玉米拔节期建立WDI梯形空间边界所需数据随生长期的变化如表3所示。表中数据均利用地面测量的数据计算得到。

表3 青贮夏玉米拔节期梯形空间所需数据Tab.3 Silage summer maize jointing stage data of VIT trapezoids

根据表3数据计算生育期第49、57、65天(拔节期3 d)边界点，根据边界点的定义，1点和3点为郁闭冠层，NDVI值为1；2点和4点为裸土，NDVI值为0。边界点陆气温差如表4所示。

表4 梯形空间边界点陆气温差

Tab.4 Boundary pointTs-Taof VIT trapezoid℃

生育期1点2点3点4点第49天-7.31.93.068.3第57天-4.11.43.851.9第65天-3.51.43.747.9

根据梯形空间边界点建立WDI梯形空间边界，从多光谱遥感影像中提取NDVI，从热红外遥感影像中提取地表温度。青贮夏玉米旱情的WDI梯形空间如图4所示。

图4 青贮夏玉米拔节期梯形空间Fig.4 Silage summer maize jointing stage VIT trapezoid

由图4可以看出，梯形空间的右侧干边可以很好地包含右侧所有的像素点，但是左侧湿边只有生育期第49天的梯形空间很好地包含了左侧所有像素点，原因是由于生育期第49天空气湿度明显低于另外两天。在3 d的梯形空间中，生育期第57天的梯形空间的值整体偏左，这是由于生育期第55天发生了降雨事件，表明梯形空间对短期降雨反应敏感,可以反映农田整体的干旱情况。

2.2 陆气温差和WDI分布图

陆气温差为陆地表面混合温度与空气温度之差，也反映作物的干旱程度。图5、6为拔节期青贮夏玉米的陆气温差和WDI分布图。由图5、6可知，生育期第49、57、65天的陆气温差和WDI分布图中均有5条明显的分界线，是为了方便田间管理在实验田种植过程中预留的裸地。在分布图中可以明显地区分裸地和玉米，且随着青贮夏玉米的生长，裸地与玉米的差异更加明显。由于生育期第55天下雨的原因，可以明显看出生育期第57天陆气温差和WDI分布图整体偏暗，证明陆气温差和WDI分布图对短期降水事件反应敏感。生育期第65天的干旱梯度最明显，可以看出1区为充分灌溉，WDI值最小。由图5、6可知，1区和3区有明显偏暗的区域，这是由于灌溉不均匀导致局部土壤含水率过高，在陆气温差和WDI分布图中形成了明显偏暗的区域。

图5 青贮夏玉米拔节期陆气温差Ts-Ta分布图Fig.5 Silage summer maize jointing stage Ts-Ta maps

图6 青贮夏玉米拔节期 WDI分布图Fig.6 Silage summer maize jointing stage WDI maps

2.3 日间尺度下WDI与陆气温差对土壤含水率的相关性分析

图7为拔节期青贮夏玉米WDI和陆气温差与土壤含水率的线性拟合，关系式如表5所示。由图7和表5可知，随着玉米生长，WDI和陆气温差与土壤含水率拟合关系式的斜率整体呈减小趋势；日间尺度下，WDI和陆气温差与土壤含水率的相关性无明显差异。生育期第65天，WDI和陆气温差与土壤含水率的相关性明显高于其他两天。原因是生育期第65天干旱严重，干旱梯度明显。当干旱越严重干旱梯度越明显，则WDI和陆气温差与土壤含水率的相关性越显著。

图7 WDI和Ts-Ta与土壤含水率(SWC)的关系Fig.7 Relationship between WDI or Ts-Ta and soil water content

表5 WDI和Ts-Ta与土壤含水率(SWC)的关系Tab.5 Relationship between WDI or Ts-Ta and SWC

2.4 日间尺度下WDI与陆气温差对气孔导度的相关性分析

图8为拔节期间青贮夏玉米的WDI和陆气温差与气孔导度的线性拟合，关系式如表6所示。由图8和表6可知，随着玉米生长，WDI和陆气温差与气孔导度拟合关系式的斜率呈减小趋势。与陆气温差相比，WDI与气孔导度拟合关系式的斜率变化更加明显。由表6可知，日间尺度下，WDI和陆气温差与气孔导度的相关性没有明显差异；对比拔节期3d的数据可知，干旱越严重，干旱梯度越明显，则WDI和陆气温差与气孔导度相关性越显著。

2.5 旬间尺度下WDI与陆气温差对土壤含水率、气孔导度的相关性分析

图9为拔节期青贮夏玉米在旬间尺度下WDI和陆气温差与土壤含水率和气孔导度的线性拟合。由图9可知，旬间尺度下监测拔节期青贮夏玉米旱情的连续变化时， WDI和土壤含水率、气孔导度的相关性优于陆气温差。原因是玉米在拔节期覆盖度和生理状态发生明显改变，由表3可知，3 d的气象因素也有明显的差异，WDI在计算过程中，考虑到了植物生理参数和气象因素的变化。证明在旬间尺度下，对拔节期青贮夏玉米的旱情进行连续监测时，WDI比陆气温差更具优势。

图8 WDI和Ts-Ta与气孔导度的关系Fig.8 Relationship between WDI or Ts-Ta and stomatal conductance

表6 WDI和Ts-Ta与叶片气孔导度(Gs)的关系Tab.6 Relationship between WDI or Ts-Ta and Gs

图9 WDI和陆气温差与土壤含水率(SWC)和气孔导度(Gs)的关系Fig.9 Relationships between WDI or Ts-Ta and soil water content or Gs

2.6 旬间尺度下不同水分胁迫的WDI和陆气温差与土壤含水率和气孔导度的相关分析

图10为拔节期青贮夏玉米在旬间尺度不同水分胁迫下数据的拟合，关系式如表7所示。充分灌溉下(1区100%充分灌溉)，WDI和陆气温差与气孔导度和土壤含水率的相关性都不高。而在不同水分胁迫下(2区65%充分灌溉，3区40%充分灌溉)WDI和陆气温差与气孔导度和土壤含水率都有较强的相关性。对比2区和3区，发现在水分胁迫不同的情况下，WDI与土壤含水率和气孔导度的相关性无明显差异，而陆气温差与土壤含水率和气孔导度的相关性却出现了较大差异。在不同的水分胁迫下，WDI能较稳定地监测青贮夏玉米的旱情，而陆气温差对青贮夏玉米的旱情监测效果却出现了较大的波动。

3 结论

(1)植被指数-温度梯形空间、WDI分布图和陆气温差分布图对短期降雨反应敏感，WDI和陆气温差分布图可明显区分青贮夏玉米拔节期的裸地和玉米，以及灌溉不均匀导致的含水率异常点。

(2)日间尺度下WDI和陆气温差与土壤含水率和气孔导度的相关性无明显差异。日间尺度下监测青贮夏玉米拔节期旱情在空间上的分布时，易于获取的陆气温差比WDI更具优势。

图10 不同水分胁迫下WDI和陆气温差(Ts-Ta)与土壤含水率(SWC)和气孔导度(Gs)的关系Fig.10 Relationship of WDI or Ts-Ta and soil water content or stomatal conductance under different water stresses

水分胁迫自变量x因变量y拟合关系式决定系数R2100%充分灌溉WDITs-TaSWCy=-0.0300x+0.07780.0889Gsy=-0.4889x+0.65590.1139SWCy=-0.0007x+0.07790.0978Gsy=-0.0112x+0.65520.108965%充分灌溉WDITs-TaSWCy=-0.0593x+0.07150.7593Gsy=-1.0105x+0.63460.8200SWCy=-0.0023x+0.07770.6939Gsy=-0.0411x+0.75450.807440%充分灌溉WDITs-TaSWCy=-0.0976x+0.07910.7283Gsy=-1.4127x+0.72200.7725SWCy=-0.0260x+0.07650.3566Gsy=-0.0412x+0.71520.4662

(3)旬间尺度下WDI与土壤含水率和气孔导度的相关性明显优于陆气温差。在监测旬间尺度下青贮夏玉米拔节期旱情的连续变化时，WDI的监测效果更好；旬间尺度下对不同水分胁迫梯度的数据进行相关性分析发现，充分灌溉下WDI和陆气温差都不能很好地监测青贮夏玉米拔节期旱情。WDI在不同的水分胁迫下可以较为稳定地监测青贮夏玉米拔节期旱情，而陆气温差的监测效果却出现了较大的波动，这表明与陆气温差相比，WDI监测不同水分胁迫的青贮夏玉米拔节期旱情有着更好的稳定性。

(4)卫星遥感数据、无人机遥感数据、地面监测数据3种数据融合进行监测，可提高监测范围和监测精度。采用无人机遥感反演WDI监测青贮夏玉米旱情，可为后期“星-机-地”数据融合的新型遥感监测方法提供了研究基础。