一种基于无人机高光谱影像的土壤墒情检测新方法

葛翔宇， 丁建丽*， 王敬哲， 孙慧兰， 朱志强

1. 新疆大学资源与环境科学学院， 新疆 乌鲁木齐 830046 2. 新疆大学绿洲生态教育部重点实验室， 新疆 乌鲁木齐 830046 3. 新疆大学智慧城市与环境建模自治区普通高校重点实验室， 新疆 乌鲁木齐 830046 4. 深圳大学海岸带地理环境监测自然资源部重点实验室， 广东 深圳 518060 5. 新疆师范大学地理科学与旅游学院， 新疆 乌鲁木齐 830054 6. 北京化工大学材料科学与工程学院， 北京 100029

引 言

土壤含水量(soil moisture content， SMC)影响土壤的理化过程， 参与全球生态、 环境、 水温和气候变化模式[1]。 同时， SMC是约束土壤养分状况的关键因素, 是影响精准农业发展的首要因素[2]。 绿洲农业承载着新疆维吾尔自治区的第一产业， 近年来的人类活动加剧干旱区内墒情失衡， 使绿洲生态环境暴露在干旱与盐渍化灾害中[3]。 在实施可持续土壤管理实践和精准农业时， 了解SMC的空间分布对于确定土壤墒情监测和土壤水盐运移的区域至关重要。 因此， 快速易行的获取SMC数据对农业监测、 产量估算以及合理灌溉具有必要的现实意义。

近十年来遥感技术的迅猛发展， 在地表观测的各领域得到了广泛应用[3-4]。 其中无人机(unmanned aerial vehicle， UAV)遥感技术的推广， 使得大尺度、 高效率地获取SMC信息成为可能[4]。 纵观相关研究的探索， 光学遥感能有效监测地物参量， 其中植被冠层可反映植被的生长状况和健康状况， 其光谱特征在不同土壤水分胁迫条件下表现出差异[5-6]。 因此， 使用UAV衍生的高光谱植被数据可作为准确评估SMC的替代方案。 然而高光谱会产生信息冗余和背景噪声等问题[7]。 光谱指数可以解决在高光谱定量估算研究中的类似问题， 易检测到敏感波长并且进一步增强目标的特定属性和光谱特性之间的相关性[8]。 通过冠层光谱指数模型， 讨论受到水胁迫的各类参数与SMC具有良好相关性[9]。 然而， 未经预处理的数据是几种复合信号与各种重叠数据的组合， 很难实现深度数据挖掘。 为弥补这一不足， 引入预处理来消除部分外部噪声， 增强非线性关系并提高地表参数估计模型的准确性。

1 实验部分

1.1 数据采集

采用DJI Matrice600 PRO®六旋翼UAV， 搭载Nano-Hyperspec®高光谱传感器。 高光谱传感器的波段范围为400～1 000 nm， 光谱分辨率为6 nm， 重采样间隔为2.2 nm。 采集时间为2018年4月17日15：00， 并对传感器进行暗电流校正及光谱定标。 后期利用Hyperspec Ⅲ及Headwall SpectralView软件处理校正。 布设70个50 cm×50 cm样方， 四点采样封装于铝盒， 并用GPS保存各点坐标， 后期在室内利用烘干法(105 ℃的恒温箱， 48 h)测定SMC。

图1 采样点分布、 高光谱传感器说明及UAV野外作业

1.2 数据预处理

在MATLAB R2016b环境中基于Savitzky-Golay(SG)滤波器平滑了高光谱图像。 以SG平滑的高光谱图像为基础影像(R)， 基于IDL + ENVI平台对原始数据进行一阶导数(FDR)， 二阶导数(SDR)， 连续体去除(CR)、 吸光度(A)、 吸光度一阶(FDA)和吸光度二阶(SDA)6种预处理方案。

1.3 适宜光谱指数构建

为了充分利用光谱数据， 根据以往的研究选择差值型指数(DI)， 比值型指数(RI)， 归一化型指数(NDI)和垂直型指数(PI)。 其中， PI的常数项是根据高光谱影像提取的土壤线的系数y=0.440 1x+0.330 8， 土壤线方程为。 4个光谱指数的数学表达式如式(1)—式(4)

DI(Ri，Rj)=Ri-Rj

(1)

RI(Ri，Rj)=Ri/Rj

(2)

NDI(Ri，Rj)=(Ri-Rj)/(Ri+Rj)

(3)

(4)

其中Ri和Rj是在高光谱传感器的工作范围内(400～1 000 nm)任意获得的i和j的光谱反射率。 使用MATLAB R2016b实现指数与SMC的相关性及指数的波段组合。

1.4 模型的构建与评价

RF回归是较常见的机器学习算法之一， 具有理想的估计能力， 特别是对于高维数据集[9]。 它也是一种基于分类回归树(CART)的集成学习算法。 该算法关键参数有决策树的数量ntree(设置ntree=500)和节点数mtry(mtry=10， 即总变量数的1/3)。 GBRT是针对回归和分类问题的机器学习技术[12]。 每次迭代都适合对前一次遗留的残差的决定， 然后通过组合树来完成预测。 考虑到同RF和XGBoost进行比较， 相关参数与XGBoost保持一致。

XGBoost是一个可扩展且灵活梯度增强的优化实现， 有关XGBoost的详细说明， 请参阅文献[13]。 为确定最优迭代次数， 防止过拟合， 选用交叉验证的方法， 设定nfold=3， 测试结果如图2所示， 当迭代次数(nround)>300时验证集的精确度将不再提高， 为了便于同RF算法公平对比， 故本文选择nround=500作为模型参数。 (其他参数设置： 学习速率(eta)=0.01， 每棵树的最大深度(max_depth)=10)。

纳入标准：①年龄18岁及以上的近视和散光患者；②球镜度＞-6.00 D，柱镜度＜ 2.00 D；③屈光度数稳定2年以上，每年变化＜0.5 D；④角膜地形图检查正常；⑤软性角膜接触镜停戴1～2周，硬性角膜接触镜停戴3～4周，角膜塑形镜停戴3个月以上。

图2 XGBoost最优迭代次数选择

上述3个模型基于R 3.5.0平台实现。 为了量化3种集成学习估算模型的性能， 使用以下指标评估模型的效果[7]： 确定系数(R2)， 均方根误差(RMSE)和相对百分比偏差(RPD)。

2 结果与讨论

2.1 SMC及光谱特征

样本划分基于联合x-y距离 (simple set portioning based on jointx-ydistance， SPXY)[14]算法进行， 选取50个样点作为建模集， 20个样点作为验证集。 全集、 建模集和验证集具有以下描述性统计结果(图3)： 全集SMC的均值(Mean)为24.446%， 标准差(SD)为5.408%。 建模集(12.23%～37.63%)的Mean和验证集(14.95%～35.62%)的Mean分别为24.499%和24.374%。 SD和Mean相似表明建模集和验证集与全集SMC保持类似的统计分布， 在确保代表性样本的同时尽可能缩减了建模集和验证集中存在偏差的估计。

图3 SMC样本划分及其描述性统计

UAV高光谱影像进行了R和FDR， SDR， CR， A， FDA和SDA共6种预处理， 结果如图4， 红线代表平均光谱， 灰色区域代表平均光谱加减标准差。 随着导数的阶数增加， 处理后的光谱的强度也会降低， 而A和CR增强了某些区域的光谱， 尤其突显了红边信息。

2.2 适宜光谱指数构建

为了可视化表达SMC和光谱指数(DI， RI， NDI和PI)之间的相关性， 通过二维相关图来表达(图5)。 图5色柱表示SMC与光谱指数之间的相关系数的平方(r2)，x轴和y轴表示400～1 000 nm的波段。 暗红色描绘了SMC和光谱指数之间的高r2。 为了更好地描述， 下面中将r2转换为相关系数的绝对值(|r|)。

如表1所示， 使用SMC建立的28个光谱指数均通过0.01水平的显着性检验(阈值为±0.380)， 其中相关性最高的4个光谱指数DI， RI， NDI和PI的|r|分布范围为0.475～0.773。 在不同预处理方案下， 光谱指数也表现出明显差异。 其中， R_PI， A_PI， CR_NDI和CR_RI都在0.75以上， 这四个指数的效果比较显著。 基于FDR， SDR和SDA的3种预处理方案构建的光谱指数效果较差。 不同的预处理方案可以不同程度地改善光谱指数与SMC之间的相关性， 最优光谱指数为A_PI(|r|=0.773)。

2.3 构建模型与评价

图4 基于UAV高光谱影像的预处理结果

图5 SMC与最优光谱指数的r2图

Table 1 |r| between SMC and spectral indices based on different pretreatments

预处理光谱指数|r|预处理光谱指数|r|预处理光谱指数|r|预处理光谱指数|r|RDI(R479, R619)0.724NDI(R431, R446)0.748RI(R431, R446)0.747PI(R446, R471)0.772FDRDI(R435, R746)0.662NDI(R702, R726)0.674RI(R702, R724)0.668PI(R435, R744)0.693SDRDI(R710, R753)0.551NDI(R417, R753)0.487RI(R444, R895)0.475PI(R653, R753)0.554CRDI(R400, R446)0.737NDI(R431, R446)0.755RI(R431, R446)0.755PI(R446, R466)0.746ADI(R431, R446)0.748NDI(R431, R619)0.725RI(R461, R619)0.720PI(R446, R471)0.773FDADI(R435, R744)0.742NDI(R513, R726)0.616RI(R420, R726)0.624PI(R435, R713)0.738SDADI(R579, R753)0.577NDI(R677, R753)0.561RI(R440, R446)0.424PI(R753, R946)0.569

表2 不同SMC预测模型建模效果

Table 2 Calibration and validation results for SMC estimation based on different modeling strategies

模型建模集验证集R2calRMSECR2valRMSEPRPDGBRT0.8972.4190.8872.1032.002RF0.8921.9930.9052.0582.109XGBoost0.9252.1980.9261.9432.556

图6 不同模型的验证结果

此外， 集成学习能通过重要性来演绎各膺选变量对模型结果的影响[12]。 如图7所示， 选取了前12个变量经行重要性表达， A_PI指数在3种模型中均位列第一， 可见A_PI指数的在与SMC非线性关系中的重要性。 在4类适宜光谱指数(DI， RI， NDI和PI)中， 由于PI在构建指数时考虑到土壤背景对植被光谱的干扰， 故而指数的效果在众多指数中名列前茅。 DI， NDI和RI指数效果与前人研究具有相似性[5]， 相较而言， PI对UAV高光谱影像消除背景噪声方面略胜一筹。

简单光谱指数仅考虑光谱与物体之间的相互作用， 而不考虑反射光谱之间的相互作用[9]。 由于受到高光谱传感器所限， 主要光谱集中在400～1 000 nm。 农田植被覆盖度高， 难以直接通过土壤光谱信息进行估算。 由于干旱区农业植物在不同程度的水分胁迫下， 作物冠层叶绿素随干旱程度波动， 故而土壤水分含量与叶绿素之间存在强烈的正相关关系。 基于实际效用考虑， 充分利用植被冠层光谱信息作为媒介， 间接监测土壤含水量。 根据适宜光谱指数所得到的波段组合(表1)较集中的处于420， 440， 460和700 nm附近。 420～460 nm附近是叶绿素、 类胡萝卜素和水分的强吸收带[5]， 由此佐证了指数构建的合理性。 本研究中确定的信息光谱带和区域可为传感器的光谱通道设计提供参考， 以监测干旱区的农业土壤墒情。

利用6种预处理方式(FDR， SDR， CR， A， FDA和SDA)对UAV高光谱影像的处理达到了一定效果。 由于一维光谱表达的局限性， 引入由Noda[16]开发二维同步相关光谱， 以便检测到更多难以在一维光谱中发现的光谱信息。 显然在图8中， 这些二维同步光谱图的对角线上存在一些自相关峰。 他们是证明两者间对外部扰动有协同响应。 对比不同预处理下的自相关峰， 一阶导数和二阶导数虽能剔除大量无关信息， 并确定更小范围的自相关峰， 但同时也损失了更多协同响应的光谱信息。 在R， CR和A中比较， R对光谱信息表达能力不如CR； A方案的光谱信息表达最佳。 在A中， 出现四个自相关峰， 分别位于450 nm附近， 670 nm附近， 740 nm附近以及980 nm附近， 说明在四个区域植被光谱发生显著变化。 这与前文光谱指数合理性的结果相似， 在证明预处理效果的同时， 佐证了本文光谱指数的响应机制。

图7 不同集成学习的变量重要性

图8 不同预处理方案的二维同步相关图

图9 SMC空间数字制图

机器学习算法被广泛地应用在地表属性的估算中， XGBoost估算SMC模型也在空谱结合的空间数字制图中得到较好的结果(图9)， 图9(a)与(b)表现出相似的空间分布特征， 残差最大为4.347%。 这表明该模型可以很好的对SMC可视化表达。

利用UAV高光谱影像数据， 在不同SMC估算模型下， 收获理想的估算结果， 在“星-空-地”遥感系统应用中尝试一种可替代方案， 并为低空遥感监测土壤墒情进行了有益探索。 因条件所限， 未能根据植被物候及土壤养分进行多时相影像数据采集， 所建立的SMC的机器学习估算模型迁移能力需进一步完善。 因此， 进一步的科研工作将着手于时间和区域尺度的SMC与植被光谱的研究， 为应对日益恶化干旱区生态环境提供精准农业的技术指导， 为相关部门提供作物长势监测、 估算产量、 合理灌溉及旱情监测的科学方案。

4 结 论