基于无人机光谱遥感的田块尺度蒸散发空间分布估算

苏宝峰 王 琮 张茹飞 陈 山

(1.西北农林科技大学机械与电子工程学院，陕西杨凌 712100；2.农业农村部农业物联网重点实验室，陕西杨凌 712100)

0 引言

蒸散是农业水资源管理中灌溉管理和产量预报的重要参考指标[1-3]。准确的区域蒸散信息可以为区域水资源的管理和农田作物培养提供信息以及决策支持[4-6]。然而，目前缺少对农田蒸散量空间模式的有效量化方法，获取精准的农田蒸散量分布图对于开展精准农业研究、灌溉作业非常必要[4]。

目前，获取农田蒸散量的方法主要是涡动协方差法和卫星遥感估算法[7-8]。涡动协方差法结果较为准确，但数据不具备空间尺度差异性[8]；使用卫星数据，估算结果精度低、分辨率低与实时性差是制约卫星技术在农业中应用的重要因素[9-11]。无人机技术监测土壤和植被信息具有快速、宏观、方便等特点，已被广泛应用于农业领域[11-12]。因此，基于无人机技术可以建立适用于农田的蒸散模型，然而合适的模型和数据的获取是建立农田蒸散量分布图必须解决的问题。现存遥感蒸散模型大多基于卫星数据建立，根据无人机数据建立的蒸散模型较少，需要更深入的研究[13-15]。根据分析通量交换方式的不同，遥感蒸散量模型可分为两类[16-19]：一是将土壤和植被当做均匀单层进行研究的单层模型，二是考虑土壤和植被水热传输特性差异及相互作用的多层模型。近年来，适用于稀疏植被覆盖条件下土壤和植被能量交换的多层模型，尤其是双层模型被广泛运用[20]。SAMANI等[21]已经证明，RSEB模型可用于植被分布不均匀地区潜热通量或蒸散量的估算，这是一种适用于植被分布不均匀区域的蒸散模型。ORTEGA-FARAS等[14]使用RSEB模型估算了巴西橄榄园中的蒸散量，由于对多光谱数据利用不足，使用经验公式获取土壤热通量，导致适用性非常差，对于农田并不适用。

在数据获取方面，反演蒸散量需要温度数据和多光谱数据[22-24]，无人机上搭载小型热像仪和多光谱相机可以获取相关数据。在使用无人机的小型热像仪时，由于飞行中距离目标远等原因极易出现较大的测量误差，从而影响最终的估算结果，通过校正可以有效减小测量误差[25]。

本文在分析以上方法与模型优缺点的基础上，利用无人机搭载热像仪获取的地表温度(LST)以及多光谱相机获取的归一化植被指数(NDVI)进行蒸散量的估算，提出将RSEB模型的土壤热通量计算方式进行基于多光谱数据的改进，并对热像仪数据进行温度校正，建立农田蒸散量分布图，最后，利用涡度系统的通量观测数据验证结果的准确度，为开展大范围的农田水分监测和精准农业灌溉提供支持。

1 材料和方法

1.1 研究区概况

研究区位于中国旱区节水农业研究院的小麦种植试验场(34°17′58.47″N，108°4′2.93″E，海拔525 m)，年均蒸发量1 500 mm，土壤质地为中壤土，1 m土层田间持水率为23%～24%(质量含水率)，凋萎含水率为11%～12%(质量含水率)，平均干容重为1.4 g/cm3。试验期田间种植冬小麦，在2018年10月播种，试验区施肥和灌溉水平与该地区大田一致。试验区地形开阔，附近多为农田，气象监测站位于麦田中央。

1.2 硬件系统

基于蒸散模型估算蒸散量的数据需求进行硬件系统搭建，硬件系统平台为M100型大疆多旋翼无人机，无人机质量2 431 g(含电池)，最大起飞质量3 600 g，最大可承受风速10 m/s，最大航行速度22 m/s,无人机平台如图1所示。使用Pix4Dcapture进行无人机任务规划，飞行高度为15 m。在无人机底部通过云台搭载热像仪和多光谱相机，顶部为全球定位系统(GPS)模块和光照传感器。

图1 无人机平台

试验所需的温度数据由热像仪获得，热像仪型号为FLIR VueProR，红外图像分辨率为640像素×512像素，测量范围-20～50℃，保存的数据为TIFF格式图像[26]。飞行前使用地面站设备通过蓝牙连接热像仪进行辐射校准。

试验所需的多光谱数据由多光谱传感器获得，多光谱传感器为美国Micasense RedEdge多光谱成像仪[27],镜头焦距5.5 mm，视场角47.2°。波段数5个,分别为蓝(中心波长475 nm)、绿(中心波长560 nm)、红(中心波长668 nm)、红边(中心波长717 nm)和近红外(中心波长840 nm)；单幅影像分辨率为1 280像素×960像素。另配有日光照度计和GPS模块,总质量仅170 g，适合小型无人机搭载，飞行前通过WiFi连接多光谱传感器，使用光谱标准面板的图像作为参考进行辐射校准。

1.3 数据获取和处理

1.3.1气象数据和冠层数据

在具有代表性的冬小麦夏玉米轮作旱田下垫面安装涡度相关系统(Open path eddy covariance，OPEC)自动观测仪器,如图2所示，该系统主要组成部分为：CAST3A型三维超声风速仪、LI75600A型开路CO2/H2O分析仪、CR1000型数据采集器、HMP-60型空气温湿度探头、LI200SZ型辐射量表和HFP01SC型热通量板等。

图2 涡度相关系统

SunScan冠层分析仪是一款简便测量和分析冠层中入射和透射光合有效辐射(PAR)的系统，提供关于影响田间作物生长限制因素的有价信息。研究期间使用SunScan在麦田中随机取4个1 m×1 m的方块，每个方块在不同位置测量9次取平均值，每7 d测量叶面积指数(LAI)。

1.3.2低空遥感数据获取

归一化差异植被指数(NDVI)的计算公式为[28]

(1)

式中NIR——近红外波段反射值

R——红光波段反射值

相机和无人机之间的时间需同步，以便将记录的GPS和旋转位置与每个图像链接，为后续图像拼接提供支持。使用Pix4Dcapture在移动设备上进行飞行任务设置，由于热红外相机不具备定位设备，热红外相机与多光谱相机设置相同的拍照频率和开启时间，使用同一个地理坐标数据，本研究中拍照间隔为1 s。

为了保证数据清晰并去除边缘影响，设置图像数据区域比研究区域边界多50 m，并在数据处理时将多余部分去除。对原始遥感数据坐标进行校正，通过GIS技术对数据进行配准处理[29]，选取同一时间拍摄的LST和NDVI图像，首先对两幅图像进行特征提取得到特征点；通过进行相似性度量找到匹配的特征点对；然后通过匹配的特征点对得到图像空间坐标变换参数；最后由坐标变换参数进行图像配准。

2019年夏季，每2 d取一组数据，其中在阴雨天气停止飞行任务。对重叠率进行对照试验，航向、行间重叠率在60%～95%时每隔5%取值，分析拼图后的准确度，取拼图效果最佳的值为试验标准。结果表明，在大田环境中，多光谱传感器在航向重叠率85%、行间重叠率75%时效果最优，热像仪航向重叠率85%、行间重叠率95%时效果最优。

每组数据包括不同重叠率的两次飞行，在温度变化比较明显的日期适当增加飞行任务。每次飞行都要采集温度与多光谱数据，在5月和6月共飞行80次。飞行时间如表1所示，表中时间采用中国标准时间。

表1 UAV飞行任务

1.3.3无人机数据校正

考虑到温度作为模型中重要的原始数据，必须保证其准确性以减小误差，参考METRIC模型中冷热点的校正方法，在有代表性的区域取点将实际值与热像仪返回的测量值进行对比，通过建立的模型对地表温度进行校正，减少由于设备和辐射等外界因素产生的系统误差。

METRIC模型中通过选择图像上的“极冷点”(植被覆盖度高,地表温度最低)和“极热点”(植被覆盖度低,地表温度最高)求解线性回归系数。通过假设地表温度与温差的线性关系进行参数标定。

本研究在温度分布线性关系的基础上进行温度校正。根据温度分布的线性关系，选取能代表温度分布情况的地点进行标定[30-32]，选取的点应当具有差异性、代表性，为温度校正提供支持。测量极干裸土地表的温度(裸土)、极干完全植被覆盖地表的温度(杂草)、土壤水分充足条件下的裸土(土壤)和完全植被覆盖地表的温度(冠层)这4种温度代表地区的温度值。使用UT302C型工业手持测温仪进行温度测量，并利用地表辐射温度和组分温度线性组合关系得出温度的变化规律。对实际温度与测量温度进行建模，进而校正通过热像仪获得的地表温度。

1.4 蒸散量模型

METRIC模型和RSEB模型估算地表瞬时蒸散量是通过计算地表能量平衡相关组分获得潜热通量LE，进而估算蒸散量，计算公式为[13,21]

(2)

式中Rn——净辐射通量，W/m2

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

H——显热通量，W/m2

LE——潜热通量，W/m2

ET——蒸散量，mm/h

ρw——水密度，kg/m3

δ——蒸发潜热，J/kg

1.4.1METRIC模型

METRIC模型优点在于利用地表辐射温度代替空气动力学温度,将显热通量H与地表辐射温度和空气温度之间的差值通过空气动力学阻抗建立联系，对温度进行内部调整，但是缺点在于只在植被覆盖均匀处才能有一个较好的结果，农田中植被和裸土交错分布对结果产生极大的影响。

METRIC模型中Rn计算公式为[14]

Rn=(1-α)Rs↓+RL↑-RL↓-(1-ε)RL↓

(3)

式中α——地表反射率

Rs↓——下行太阳短波辐射，W/m2

RL↑——上行太阳长波辐射，W/m2

RL↓——下行太阳长波辐射，W/m2

ε——地表比辐射率

土壤热通量在热量平衡中是一个相对较小的量，直接计算较为困难，METRIC模型估算土壤热通量和显热通量公式为[14]

(1-0.98NDVI4)Rn

(4)

(5)

(6)

(7)

式中Ts——地表温度，K

ρ——空气密度，kg/m3

CPair——空气定压比热容，J/(g·K)

z1——裸地地表粗糙度，m

z2——气象站的高度，m

ΔT——z1与z2处的空气温度差，K

rah——热量传输空气动力学阻抗，s/m

u*——像元的摩擦风速，m/s

k——冯卡门常数

u200——混合层风速，m/s

zom——像元地表粗糙度，m

1.4.2RSEB模型

RSEB模型[18]可用于估算地表为部分植被覆盖地区的潜热通量。RSEB模型为双层模型，和单层模型相比可以更准确描述植被分布不均匀地区的地表通量。RSEB模型的Rn计算公式为[14,22]

Rn=frRnc+(1-fr)Rns

(8)

其中

Rnc=(1-αc)Rse+Lin-Loutc-

(1-εc)Lin

(9)

Rns=(1-αs)Rse+Lin-Louts-

(1-εs)Lin

(10)

式中Rnc——冠层净辐射分量，W/m2

Rns——土壤净辐射分量，W/m2

fr——最低覆盖率

Rse——瞬时入射短波辐射，W/m2

αs——土壤反照率

αc——土壤和冠层反照率

Lin——瞬时传入长波辐射，W/m2

Loutc——冠层瞬时出射长波辐射，W/m2

Louts——土壤瞬时出射长波辐射，W/m2

εc——冠层表面热发射率

εs——土壤表面热发射率

G=0.323 6Rn-51.52

(11)

显热通量的计算公式为

H=frHc+(1-fr)Hs

(12)

(13)

(14)

式中Hc——冠层显热通量分量，W/m2

Hs——土壤显热通量分量，W/m2

rac——冠层上方空气动力学阻抗，s/m

raa——冠层与参考水平之间的空气动力学阻抗，s/m

ras——土壤与冠层之间的空气动力学阻抗，s/m

Tc——冠层温度，K

Ta——空气温度，K

1.4.3土壤热通量的RSEB模型改进

无人机多光谱遥感获得的数据量大而且光谱分辨率高，可以有效表现出地物特性。ORTEGA-FARAS等[14]使用RSEB模型估算蒸散量时没有引入多光谱数据，G通过经验公式得到，由于使用环境只种植单一的植被，所以不具有普遍适用性，只适用于特定环境。

为了使估算方法更好地适用于大田试验环境，分别考虑METRIC和RSEB模型的优点，因为METRIC模型对土壤热通量的计算是基于多光谱数据进行的，可以反映不同环境下的下垫面差别，所以选择式(4)进行土壤热通量的计算[24]，使土壤热通量的估算能够适用于不同的环境。

引入NDVI数据计算土壤热通量可以更好地体现空间的差异性，更好地适用于农田尺度。基于多光谱数据改进的RSEB模型比原始模型适用于更复杂的环境。

1.5 精度验证

通过平均绝对误差(MAE)、均方根误差(RMSE)和决定系数(R2)3种指标验证数据的可靠性，分别分析不同蒸散量模型和温度校正的结果。MAE和RMSE越小，R2越接近1，说明结果越准确。

2 结果与分析

剔除飞行任务中由于环境和设备影响而产生的异常值，选取其中26组飞行测量数据进行对比试验，选取数据的时间序列是在试验期间均匀分布的。将得到的遥感数据进行拼接裁剪等预处理，之后导入GIS中进行计算，得到如下数据：NDVI分辨率约为0.01 m，LST分辨率约为0.1 m，最终净辐射、土壤热通量、显热通量和潜热通量的空间分布图以0.1 m的分辨率获得。

2.1 无人机估算与OPEC测量值对比

2.1.1METRIC模型与RSEB模型比较

试验区为大量冬小麦覆盖的规则田地，在最初的METRIC模型中对地表温度的计算是采用热红外辐亮度通过对比辐射率的校正所得[14]，本研究中使用热图像直接得到LST信息。

图3 G的测量值与估算值比较结果

由于G在两种模型中采用的计算方式相同[24]，通过与OPEC系统测量值进行对比，验证结合多光谱数据计算农田中G的准确度，对比结果如图3所示。对数据分析后得出，RMSE为5.192 W/m2，MAE为4.511 W/m2，R2为0.94，G的估算结果接近测量值，这种对于G的估算方式适用于无人机数据和大田环境。

然后对比H和LE，LE可以直观地反映蒸散量的大小[13]，将单层模型的计算结果与双层模型计算结果分别与涡度系统测量数据对比，结果如图4所示。

图4 LE和H的测量值与估算值比较结果

分别对METRIC模型和RSEB模型进行分析，RMSE和MAE如表2所示，可以发现，METRIC模型相对于RSEB模型在估算LE时估算值平均超过实际值20 W/m2，这是因为单层模型在估算时没有区分土壤和植被的差异。

表2 METRIC模型和RSEB模型计算结果误差

在估算LE时，METRIC模型的R2为0.779，RSEB模型的R2为0.788，说明模型估算结果与真实值一致且RSEB模型的结果更为接近，误差分析如表2所示。

双层模型的结果一定程度上优于单层模型，这是因为双层模型考虑了冠层与土壤的差别，但总体结果还是有一定的误差，除了环境因素产生的误差外主要误差来自热像仪的测温误差，飞行的高度和速度都会对其造成影响[26]。通过校正对原始数据进行优化，以减小误差，使估算结果能够更好地接近测量结果。

2.1.2温度校正对结果的影响

由于飞行任务耗时较长，考虑到温度随时间变化，所以飞行前后对不同位置选取点分别取3次数据的平均值作为实测值，对热像仪与测温仪的值进行比较，平均误差为1.73℃，部分温度数据示例如表3所示。

表3 实际温度与热像仪测量值比较结果

在试验区域内对经过校正后的LST随机取点与实测值进行比较，平均误差为0.65℃，说明校正方法对于热像仪数据是可靠的。温度分布图如图5所示，中心部分设备为涡度系统。

图5 热像仪测量的地表温度和校正后地表温度

根据温度校正改进模型,先估算显热通量H，进而推算潜热通量LE，将经过校正的数据与未经过校正的数据结果分别与涡度系统测量数据对比，结果如图6所示。

图6 瞬时LE和H的测量值与估计值比较结果

经过温度校正后，估算LE时METRIC模型的R2为0.848，RSEB模型的R2为0.940，说明经过温度校正后RSEB模型的结果与真实值相关性是最好的，误差分析如表4所示。

对比数据明显看出，校正温度数据后的H与LE的RMSE比校正温度数据前平均降低了9 W/m2，MAE平均降低了6 W/m2，R2平均提升了0.12，校正后的估算结果更为合理准确，说明校正方法对于无人机的热红外数据非常有效，经过校正后的温度数据更贴合实际温度。本研究校正方法可以减小无人机热像仪数据获取过程中的误差，使结果更为准确。研究结果表明，温度校正后的RSEB模型估算得到的通量值和监测得到的通量值有很高的一致性，R2为0.940，RMSE为40.202 W/m2，MAE为26.017 W/m2, 可以通过本研究方法获取农田中的LE。

表4 校正后的METRIC和RSEB计算结果误差

2.2 蒸散量空间分布图

将获取的数据代入模型对遥感数据进行计算，NDVI如图7a所示，得到净辐射通量Rn、显热通量H、土壤热通量G数据分别如图7b～7d所示。

图7 NDVI和通量的空间分布图

根据式(1)和求得的通量数据，进而求得潜热通量LE如图8所示，通过估算得到的潜热通量与测量的潜热通量非常一致，生成热通量的空间分布图进而得到蒸散量的空间分布图，如图8b所示。

图8 潜热通量和蒸散量的空间分布图

通过现有其他类型技术很难得到相同结果，涡动协方差系统测量不能得到空间尺度差异性的数据[8]，使用卫星数据的遥感估算模型无法得到空间分辨率较高数据，并且受到时间条件约束[9]。通过无人机数据最终得到ET如图8b所示，图中可以清楚地看到不同位置蒸散量的差异，本研究结果在空间尺度上表示的农田蒸散量可以达到分米级。

3 结论

(1)搭建了获取蒸散模型基础数据的无人机平台，并将无人机数据代入蒸散模型进行研究，比较METRIC模型与RSEB模型在农田环境下的结果，对RSEB模型的土壤热通量计算方式进行基于多光谱数据的改进，校正无人机温度数据，最终得到农田蒸散量分布图。

(2)对比典型单层模型和典型双层模型在农田中的适用性发现，双层模型的结果更优。在通量计算方面，对RSEB模型的土壤热通量计算方式进行基于多光谱数据的改进，使模型的适用性更强。在LST获取方面，对无人机测温方式进行优化，基于FLIR VueProR热像仪多次试验，取最佳飞行重叠率与高度，并进行温度校正，即选取有代表性的温度点，分别位于裸地、冠层、杂草、土壤，随机采集5个点的真实温度，与无人机全局测温图进行建模，获得校正后温度。温度校正有效减小了测量中产生的误差，使热像仪提供的数据更为准确。

(3)将无人机遥感、涡度相关技术、地理信息技术相结合，通过无人机获取温度和多光谱数据，GIS软件进行数据处理，建立有效的蒸散模型，获取蒸散量结果。将结果与涡度系统测量真实数据进行比较，结果表明，基于无人机与气象数据建立的蒸散模型适用于大田作物，空间分辨率较高(分米级)，能够实现蒸散量的反演。