基于超宽带雷达和多光谱数据融合的土壤含水率检测

2021-10-13白清源郭文川

农业机械学报 2021年9期

郭交白清源郭文川,3

(1.西北农林科技大学机械与电子工程学院，陕西杨凌 712100；2.陕西省农业信息感知与智能服务重点实验室，陕西杨凌 712100；3.农业农村部农业物联网重点实验室，陕西杨凌 712100)

0 引言

土壤含水率是农业、生态、气象和水文等学科的重要参数之一[1-2]，对农作物生长具有重要作用，实时监测农作物的含水率情况，不仅对评价农作物健康水平有着重要意义，而且在节水灌溉、农作物光合作用和养分吸收以及产量预测等方面都具有重要作用[3-4]。

当前测定土壤含水率的方法，如烘干法[5]、中子法[6-7]、时域反射法(Time domain reflectometry, TDR)[8-9]、热脉冲法[10]等，能够较为准确地进行单个土样或者单点的小尺度测定。但数据时效性较差，很难保证测量数据的同步性，而且只能用一个或几个点的数据来表征大范围面的数据，若直接用于大面积测量，则耗费大量人力物力[11]。近年来随着星载遥感技术的迅速发展，基于卫星数据来反演土壤含水率已成为一种解决土壤含水率监测问题的新方法[12-16]，但受到当前卫星传感器的限制，在对局部农田进行反演时精度比较低，同时该方法受卫星重访周期和过境时天气影响比较大，很难有效指导精准农业生产活动[17-18]。低空扫描法是一种将雷达架设于地面一定高度处观测一定范围土壤含水率的方法[19-20]，这种方法操作简单，测量效率高并保持很高的精确性，比较适合于表层土壤含水率快速检测，但传统上与低空扫描法配合使用的雷达尺寸过大，测量难度和成本较高[21]。近些年来，由于超宽带(Ultra-wide band，UWB)雷达体积小、质量轻、穿透力强、抗干扰和低功耗的特性[22]，被越来越多地应用于土壤含水率相关的研究中，而且将超宽带雷达模块与低空扫描法相结合，可以简化低空扫描法中土壤信号采集的难度。

一些学者利用模糊逻辑系统、人工神经网络等手段对从超宽带雷达测量数据中得到土壤含水率相关信息进行了研究[23-27]，所选择的试验条件都是平坦裸土。但在实际的农田环境中，土壤上均会存在不同的植被覆盖且土壤中石子与秸秆残渣等干扰较多，这些都会影响土壤含水率监测的精度。而微波遥感不如光学遥感形象直观，易于理解。植被的雷达后向散射与系统参数(波长、入射角、极化方式)、植被参数以及地表参数均有关,电磁波与植被之间的相互作用是比较复杂的。为了对这个复杂的过程进行理解并提高地表参数的反演精度，研究人员通过对植被微波后向散射特性的研究，建立了各种类型的植被微波后向散射模型[28]。目前已有学者直接通过植被指数来消除植被覆盖产生的影响[29]。

基于以上分析，本文以土壤含水率为研究对象，在10块不同植被覆盖情况的试验田中进行不同土壤含水率的超宽带回波信号采集，试验中土壤含水率均采用体积含水率(Variable water content，VWC)计算。通过无人机多光谱遥感获得10块试验田的多光谱图像，进而计算其对应的归一化植被指数(NDVI)，以此来消除植被覆盖所产生的雷达后向散射影响。将无人机多光谱图像与超宽带雷达模块(PulsOn 440, P440)的低空扫描法相结合，融合多光谱与超宽带雷达数据，利用支持向量机(SVM)模型实现对不同程度植被覆盖土壤含水率进行分级预测。在此基础上，研究对比融合多源数据与只采用超宽带雷达回波数据的预测精度，并建立农田尺度上最优的土壤含水率预测模型。

1 研究区域与试验方法

1.1 研究区域

试验地点为陕西省咸阳市杨凌区某处试验田，土壤类型为陕西省关中地区的塿土，总面积约为100 m2，总体植被覆盖度约为42.10%，植被覆盖高度为0～30 cm。选择其中不同植被覆盖的10块试验地，植被覆盖类型为杂草，每块试验田的大小约为2.5 m×2.5 m。图1a为试验区域总体示意图，图1b、1c为具体地块试验过程中无人机定点拍摄的多光谱影像RGB波段融合后图像。

1.2 试验设备

(1)超宽带雷达模块

选用的超宽带雷达模块是美国Time Domain公司的PulsOn 440单基站雷达模块，实物图如图2所示。它是一种工作频率在3.1～4.8 GHz之间的超宽带雷达，能够采用双向飞行时间方式在2个以上的模块之间进行精确测距与通信，也可以作为单基地雷达、双基地雷达或者多基地雷达工作。P440单基站雷达模块功耗低，供电电压为直流4.5～48 V，小型移动电源就可以满足其长时间稳定的供电需求，同时其尺寸小、质量轻的特点使其便携性很强，同时具有很强的穿透力。这些特性使其可以很好地适应农田间复杂的实际农情环境，方便田间试验时数据的采集，完成长时间高效率测量任务。

(2)雷达天线

雷达天线选用P440标配的BroadSpec天线，实物图如图2所示。该天线是一种椭圆形的平面偶极子天线，通过标准SMA公接口连接到P440模块的两个天线端口，能够提供全方向发射和接收模式。其时域特性和频域特性十分稳定，适用于P440模块的信号传输要求，有效射频带宽范围为3.1～5.3 GHz，天线效率约为90%。

(3)土壤含水率测定仪

使用顺科达TR-6D型土壤温湿度计测量土壤的体积含水率。TR-6D是一种便携式土壤温湿度测量仪，具有体积小、质量轻、能耗低的优点，可以很好地适应农田间的复杂实际情况。含水率测量范围为0～100%，精度为±2%。TR-6D采用接触式的测量法，所用探针长度为70 mm，垂直安装时传感器可测量近60 mm长感应区间的土壤体积含水率。

(4)无人机与多光谱相机

采用搭载多光谱相机的无人机获取不同植被覆盖情况试验地块的多光谱图像。所用无人机型号为四旋翼的大疆M100，搭载的多光谱相机为MicaSense的Red Edge五通道多光谱相机，其光谱参数如表1所示。该相机带有光强传感器和灰板，可以在无人机航拍飞行过程中矫正外界光线变化对光谱造成的影响，并对航拍影像进行反射率矫正[30]。

表1 Red Edge的光谱参数Tab.1 Spectral parameters of Red Edge nm

1.3 试验方法

本文采用低空扫描法进行土壤回波信号的采集。低空扫描法是利用分析接收的微波信号与土壤样品含水率之间的关系来反演土壤含水率。基于低空扫描法的数据测量系统由P440超宽带雷达、发射天线、土壤样品、接收天线、计算机和其余辅助装置构成。待测土壤区域处于发射天线和接收天线中间正下方，通过P440超宽带雷达模块完成回波信号的发射与接收，计算机用于回波信号的储存和数据分析，如图3所示。

1.4 试验步骤

(1)考虑到雷达天线耦合噪声的影响，在采集测量过程中 P440雷达应高于地面0.8 m以上[31]。选取周围比较空旷的试验田进行试验，搭建好雷达固定支架后，将P440雷达放置在距离地面1.0 m的高度位置进行回波信号的采集，通过USB接口与计算机连接，利用MRM软件控制回波信号数据采集的开始和停止，并将接收信号数据直接储存在计算机中以便进行后期数据处理。

(2)使用P440超宽带雷达在10块不同植被覆盖情况的试验田中进行不同含水率的土壤雷达回波数据采集。共采集了40个不同含水率的土壤雷达回波信号，土壤含水率区间为14.10%～28.59%，符合实际农田土壤含水率的情况。将土壤含水率从14.01%开始按照2%的差值分为8级，其中第1级为14.01%～16.00%，第2级为16.01%～18.00%，以此类推第8级为28.01%～30.00%。每个含水率含有1 000个回波信号样本，按照3∶1的比例将750个样本作为训练集，250个样本作为测试集。

(3)使用TR-6D型土壤温湿度计测量土壤的体积含水率，测量区域为P440雷达下方的1.5 m×1.5 m矩形区域。测量时选择区域内不同位置的12个点进行测量并记录数据，去除最低值和最高值后对10个点的体积含水率取平均值作为此次试验土壤的实际体积含水率。

(4)使用搭载多光谱相机的无人机对试验田进行拍照，获得试验区域的无人机多光谱遥感图像，用于计算不同试验地块对应的植被指数。

2 数据处理

2.1 雷达回波信号提取

P440超宽带雷达采集的回波信号中，每个回波的离散时间长度为480 ns，首先需要从回波信号中提取出含有土壤含水率信息的有用信号。根据雷达波时延公式，可以得到最先接收到含有土壤表层信息的回波信号的时间，从而在所有信号中提取出含有土壤信息的回波信号的起始点，公式为

(1)

(2)

(3)

式中V——雷达回波信号传播速度

R——雷达波传播距离，取1.0 m

τ——P440的信号采样频率，取61 p/s

Δt——传播时间m——时间指数

C——光速ξr——介电常数

解得m约为109，表明每个回波的离散时间长度480 ns中，从109 ns起才包含有与土壤含水率有关信息，因此舍弃前108 ns，将样本的离散时间长度缩小为372 ns。雷达波在土壤中的穿透深度约为0.5 m[32]，其在土壤中的传播速度如式(3)所示,其中土壤的介电常数介于4～40之间，由经验公式可得含水率30%土壤的介电常数约为20[33]，可以通过式(3)计算出包含土壤信息的有效回波信号采样离散时间长度约为244 ns，本试验中取250 ns，即一个回波信号样本中109～359 ns的离散时间长度上的信息是包含土壤信息的有效回波信号。以土壤体积含水率为22.28%为例，如图4所示。

本研究中植被覆盖高度约为0.3 m，即雷达与植被距离约为0.7 m。根据雷达波延时公式可以计算出每个回波的离散时间长度480 ns中，约从76 ns起包含有与植被有关信息。将8号试验田24.31%与9号试验田24.40%的雷达回波样本进行比较分析，结果如图5所示。从图中可以看出，约70 ns开始两种雷达回波产生较为明显的变化，这是因为采集的雷达回波中开始包含植被信息，与计算结果一致，后续采集到的回波中均包含植被对回波散射产生的噪声信息。

2.2 雷达回波时域特征提取

由于农田中环境比较复杂，超宽带回波信号采集过程中噪声比较严重且无法很好地直接剔除，而且UWB信号是非平稳信号，幅值较大，信息冗余较多，若直接采用提取出来的含有土壤信息的回波信号进行建模，消耗时间过长且无法避免噪声带来的影响，因此本文对每一个采集的回波信号进行时域特征提取，选用平均值、方差、最大幅值、最小幅值、峰-峰值、均方根、偏斜度、峰值因子和峭度9个时域特征。用提取出来的9个回波时域特征参数来表征原始超宽带回波信号。

2.3 植被指数计算

根据植被的光谱特性，将多光谱相机拍摄的可见光和近红外波段进行组合，计算出植被指数NDVI。NDVI可以消除大部分与仪器定标、太阳角、地形、云阴影和大气条件有关辐照度的变化，是目前已有的40多种植被指数中应用最广的一种[34-35]。

计算得出10块试验田的NDVI分别为0.002 3、0.000 3、0.037 8、0.268 7、0.349 6、0.364 1、0.006 8、0.123 4、0.572 8、0.109 0，本研究中1、2、3、7号试验田植被覆盖十分稀疏，4、5、6、8、9、10号试验田有不同程度明显植被覆盖。

2.4 支持向量机模型

支持向量机(Support vector machines，SVM)可用于模式分类和非线性回归，主要思想是建立一个分类超平面作为决策曲面，使得正例和反例之间的隔离边缘被最大化。其特点是使用核函数、稀疏解、边界追踪控制边界和支持向量数目。SVM的计算复杂度不依赖于输入空间的维数，同时具有良好的泛化能力，分类精度比较高。在SVM中，输入样本x首先通过非线性映射Φ(x)映射到一个高维的特征空间，然后在这个特征空间中建立一个线性模型来估计回归函数，公式为

(4)

式中ω——超平面的法向量

b——超平面的截距

通过引入拉格朗日函数将约束优化问题转换成对偶问题，通过结对偶函数得到式(4)的解为

(5)

nsv——支持向量的个数

K(xi,x)——核函数

2.5 模型评价指标

本文采用混淆矩阵对不同特征输入情况下的预测结果进行精度分析比较。基于混淆矩阵计算总体精度(Overall accuracy，OA)和Kappa系数。总体精度和Kappa系数都可用于对比预测结果的精度，Kappa系数可以分为5组来表示不同级别的一致性：0～0.20，极低一致性；0.21～0.40，一般一致性；0.41～0.60，中等一致性；0.61～0.80，高度一致性；0.81～1.00，几乎完全一致。

3 结果与分析

3.1 回波信号时域特征相关性分析

将超宽带雷达回波提取到的9个时域特征分别求取平均值后按照含水率从低到高进行绘图，结果如图6所示。第1组土壤回波信号实际含水率最低，第40组土壤回波信号实际含水率最高。

从图6中可以看出，随着土壤含水率的增加，9个回波时域特征均有不同程度起伏变化。灰色关联分析是对一个系统发展变化态势的定量描述和比较的方法，它反映了曲线间的关联程度，通常用于分析各个因素对于结果的影响程度。本文使用灰色关联分析对从雷达回波信号中提取出来的9个时域信号特征分别与土壤含水率进行相关性分析，相关性从高到底排序依次为峰值因子、峭度、均方根、峰-峰值、最大幅值、偏斜度、方差、平均值、最小幅值。峰值因子与土壤含水率关联度最高，为0.822 9，而最小幅值与含水率的关联度最低，为0.756 2，图6中最小幅值随着含水率升高的变换曲线相对较平缓，波动不明显也可以验证灰色关联分析结果的正确性。

将9个时域特征按照相关度由高到低排列，前3个特征组合为a1，即峰值因子、峭度和均方根的组合，中间3个特征组合为a2，即峰-峰值、最大幅值和方差的组合，后3个特征组合为a3，即偏斜度、平均值和最小幅值的组合。

3.2 不同特征组合的模型精度分析

根据灰色关联分析，按照相关度从高到低，本文分别将最高的3个时域特征、6个时域特征和全部9个时域特征组合在一起作为SVM分类的输入特征，同时加入植被指数NDVI。输入特征组合分别为a1、a1+NDVI、a1+a2、a1+a2+NDVI、a1+a2+a3和a1+a2+a3 +NDVI，不同参数组合下SVM预测结果如表2所示，其中特征组合为a1、a1+NDVI、a1+a2+a3和a1+a2+a3 +NDVI的预测结果如图7所示。

表2 不同特征组合下SVM模型预测结果Tab.2 Classification results of SVM model under different feature combinations

3种时域特征组合中，a1特征组合的预测效果最差，总体精度为64.88%，Kappa系数为0.594 8，最优组合为a1+a2+a3特征组合，总体精度为95.59%，Kappa系数为0.949 2。a1+a2+a3特征组合作为特征输入比a1特征组合的预测精度提高了30.71个百分点，Kappa系数提高了0.354 4，比a1+a2特征组合的预测精度提高了21.56个百分点，Kappa系数提高了0.247 6。引入植被指数NDVI后，a1+NDVI组合作为特征输入总体精度为77.24%，Kappa系数为0.737 2。a1+a2+a3+ NDVI组合作为特征输入总体精度为98.09%，Kappa系数为0.978 0，在6种特征输入组合中总体精度和Kappa系数均为最高，预测效果最佳。加入NDVI后的3种不同时域特征组合分类模型比未加入NDVI的总体精度分别提高了12.36、13.57、2.50个百分点，Kappa系数分别提高了0.142 4、0.155 4、0.028 8，效果与未加入NDVI相比提升显著。不同的时域特征组合结合NDVI共同作为SVM特征输入后，总体精度和Kappa系数均有明显提高。

3.3 不同植被覆盖试验田的模型精度分析

本文中将稀疏植被覆盖的1、2、3、7号试验田分为一组，将有不同程度植被覆盖的4、5、6、8、9、10号试验田分为一组，对比将9个雷达回波信号时域特征作为SVM模型输入特征时，引入植被指数NDVI对于模型精度的影响，结果如表3所示。由表3可得：

表3 不同植被覆盖数据的模型精度Tab.3 Model accuracy of different vegetation coverage data

(1)植被覆盖十分稀疏的试验田随着太阳自然暴晒、土壤含水率降低的同时，不同试验田地表存在各不相同、不同程度的龟裂情况且石块裸露于地表。由于数据采集周期的原因，植被覆盖稀疏的试验田中地表龟裂程度变化对于雷达回波影响较大。对1、2、3、7号试验田植被指数归一化之后方差为0.218 6，而4、5、6、8、9、10号试验田植被指数归一化之后方差为0.139 1，极其稀疏植被覆盖的4块地之间差别大于有明显植被覆盖的试验田，因此a1+a2+a3特征组合作为模型输入时模型精度较差，略低于有明显植被覆盖的试验田。

(2)对于植被覆盖十分稀疏的试验田，通过多光谱遥感数据获得NDVI指数引入SVM模型中可以提高模型精度，本研究中4块近乎稀疏植被覆盖的试验田引入NDVI后土壤含水率预测精度提高了1.62个百分点。对于有明显植被覆盖的试验田，引入植被指数NDVI后土壤含水率预测精度提高了3.05个百分点。

(3)通过NDVI指数的引入，有明显植被覆盖试验田采集的数据进行分类建模后精度提升大于稀疏植被覆盖的试验田采集的数据，说明NDVI这一特征对于土壤含水率分类中植被的影响较为敏感，能够一定程度上降低由植被覆盖等因素所带来的噪声，提高模型精度。