基于无人机遥感图像纹理与植被指数的土壤含盐量反演

向友珍 李汪洋 台 翔 安嘉琪 王 辛 陈俊英

(1.西北农林科技大学旱区农业水土工程教育部重点实验室, 陕西杨凌 712100;2.西北农林科技大学水利与建筑工程学院, 陕西杨凌 712100)

0 引言

土壤盐渍化是指在自然和人为作用下,盐分在土壤表层不断积累的现象,是造成土壤退化、农业减产和生态环境恶化的重要因素[1]。河套灌区土壤盐渍化问题严重,对盐渍化土壤进行治理和开发利用是农业可持续发展与生态文明建设高质量发展的有力保障[2]。实时、高效地监测土壤盐渍化状况,是治理土壤盐渍化的重要前提。

人工定点采样法是土壤盐渍化监测的方式之一,该方法监测周期长、工作量大,且监测的精度与采样点的数量、布设方式有很大关系,很难反映盐渍化在空间上的实际分布。随着科学技术的不断发展,遥感被广泛应用于盐渍化研究和监测,为大面积时空分布动态化的盐渍化监测提供了新的途径。然而广泛可用的卫星图像无法提供高空间分辨率和数据采集时间的灵活性,在农业应用方面有很大的局限性[3]。

除卫星遥感外,利用无人机平台搭载的小型多光谱遥感设备对区域土壤盐渍化监测也是重要手段[4-5]。与卫星图像相比,基于无人机多光谱图像的盐渍化监测具有空间分辨率高、光谱分辨率强、波段连续性强等优点,可以获得多维、高精度的盐渍化检测信息,实现对土壤盐渍化的动态监测;同时,无人机传感器携带多光谱相机平台,具有实时、高分辨率、移动性和灵活性的优势,在高精度盐渍化监测方面具有优势[6]。刘旭辉等[7]利用无人机多光谱遥感影像建立的机器学习模型反演土壤含盐量,得出不同季节的土壤含盐量会有所不同。刘楠等[8]将局部区域的相关性分析和多维光谱变化特性采用二进制编码的形式来表征不同类型地物纹理特征,实现了区分和判别影像上不同纹理的目的,并证明该方法对提取地物纹理特性具有一定可行性。万亮等[9]利用无人机多光谱融合植被指数和纹理特征提高了水稻含水量预测结果;周聪等[10]通过实践证明,高分辨率遥感图像的纹理特征可以作为估算植被生长参数的更有效指标。陈鹏飞等[11]通过实践证实了基于无人机多光谱图像的剔除土壤背景和增加纹理特征处理可以提高棉花植株氮浓度的反演精度。但是基于无人机多光谱图像剔除土壤背景和纹理特征是否可以提高土壤含盐量精度方面研究较少。

本文以河套灌区沙壕渠灌域4块不同土壤盐渍化程度的典型盐渍化试验地为试验研究区,以灌区典型地物为研究对象,在灌区作物生长期内采集试验区盐渍土样本并获取土壤盐渍化信息。同时利用无人机多光谱采集遥感数据,对比研究基于盐渍化土壤和耐盐作物光谱的土壤含盐量反演模型的性能,建立基于无人机多光谱遥感平台的区域盐渍化监测体系,为区域盐渍化研究提供可靠依据。

1 材料与方法

1.1 研究区概况

本研究区域位于巴彦淖尔市沙壕渠灌域(40°52′～41°00′N,107°05′～107°10′E),隶属于河套灌区解放闸灌域,如图1所示。灌域因不合理的灌排方式,气候、土质和地貌等因素的综合影响,土壤盐渍化问题突出。区内种植作物以向日葵、玉米等耐盐粮油作物为主。沙壕渠灌域土壤类型为粉壤土、砂壤土和壤土。沙壕渠灌域面积约52.4 km2,南北跨度约15.1 km,东西横跨4.2 km,为典型的温带大陆性气候,多年平均气温3.7～7.6℃,多年平均降水量210～290 mm,多年平均年蒸发量2 100～3 080 mm,多年平均日照时长为3 000～3 200 h。

1.2 试验地布设

将试验地以土壤含盐量为基础从低到高依次分为1、2、3、4,每块试验地面积为15.5 hm2左右。4块试验地主要种植作物为向日葵和玉米。每块试验地均匀布设30个土壤采样点。采样点布设见文献[12]。

1.3 光谱图像采集和处理

使用的无人机为深圳市大疆创新科技有限公司生产的M600型六旋翼无人机,其最大上升速度5 m/s,最大下降速度3 m/s,最大飞行速度18 m/s,飞行承载质量6 000 g,飞行高度2 500 m,单次飞行时间35～40 min。多光谱遥感相机采用美国Tetracam公司生产的6通道Micro-MCA多光谱相机,包括蓝光波段、绿光波段、红光波段、红边波段、近红外1波段、近红外2波段共6个遥感波段,波长分别为490、550、680、720、800、900 nm。试验时间为2022年7月16—20日。每次试验均在11:00—14:00进行,试验日晴朗无风,以确保充分的辐射强度,尽量减小植被阴影对光谱的影响。根据提前规划好的航线,设置无人机飞行高度120 m,对应多光谱相机分辨率为6.5 cm,相机拍摄速率为18～19幅/min,每次试验均设有白板进行图像标定。

1.4 光谱指数计算

光谱指数是综合考虑地物的各波段光谱特征,对不同波段反射率进行数学变换组合,以增强地物特定的信息[13]。分别选择8种植被指数和10种盐分指数,其计算公式如表1所示。

表1 光谱指数Tab.1 Spectral index

1.5 灰度共生矩阵法

灰度共生矩阵(GLCM)是一种基于统计数据的图像纹理特征提取方法,由HARALICK等[26]提出。

GLCM的元素是不同组合出现的频数,为了计算图像纹理的特征参数,对GLCM进行归一化。已知i、j分别为(x,y)、(x+Δx,y+Δy)的像素值,(i,j)的个数为V(i,j),将归一化GLCM记为P,则

(1)

(2)

式中Vi,j——矩阵元素N——矩阵元素个数

基于归一化灰度共生矩阵,计算图像纹理特征参数,本研究采用8个纹理特征参数,其计算公式和特性如表2所示。

表2 图像纹理特征参数Tab.2 Texture feature parameters of images

1.6 Otsu算法

Otsu算法是一种计算简单、自适应强且已得到最广泛使用的图像阈值自动选取方法[29]。设某一灰度级对应阈值为T,类间方差计算公式为

(3)

式中Wb——阈值T下背景占整幅图像的比重

Wf——阈值T下前景占整幅图像的比重

1.7 全子集筛选法

全子集筛选法是利用全子集回归分析,对自变量不同的组合,用最小二乘法进行建模分析,筛选最优的变量组合。选择最优模型的评价标准为:①似然函数最大化。②模型未知参数最小化[30]。本文通过R编程语言进行全子集筛选,利用决定系数R2和贝叶斯信息准则(Bayesian information criterion,BIC)来评价筛选结果,对比分析R2越大、BIC越小的筛选结果,为最优变量组合。

1.8 建模策略

建模策略分别为:未剔除土壤背景的光谱指数(策略1)、剔除土壤背景后的光谱指数 (策略2)、未剔除土壤背景的光谱指数+图像纹理特征(策略3)、 剔除土壤背景的光谱指数+图像纹理特征(策略4)。

1.9 模型方法和评价指标

极限学习机(Extreme learning machine, ELM)是一种基于最小二乘学习算法的隐层前馈网络[31]。ELM收敛速度比传统算法快,因为它无需迭代即可学习,同时,随机隐藏节点保证了全局逼近能力[32]。本文ELM模型采用R语言elmNNRcpp包构建。支持向量机(Support vector machines,SVM)是用于监督学习的强大计算工具[33]。本研究采用非线性支持向量机模型,利用R语言e1071包构建。通过决定系数R2、均方根误差(RMSE)、标准均方根误差(NRMSE)指标评价模型精度。R2越接近1,RMSE越接近0,说明模型效果越好。

本文利用R语言编程构建SVM和ELM土壤含盐量反演模型,其中实测土壤含盐量为因变量,光谱变量为自变量。将深度0～20 cm平均土壤含盐量按照比例2∶1随机划分建模集和验证集,通过调整参数,获得每个条件下最佳模型。建模流程如图2所示。

图2 建模流程图Fig.2 Model building flowchart

2 结果与分析

2.1 土壤含盐量

对试验地总计87个采样点的土壤含盐量进行统计,将计算表层和深度10～20 cm实测土壤含盐量的平均值作为深度10～20 cm的平均土壤含盐量。将深度0～20 cm平均土壤含盐量按照比例2∶1随机划分建模集和验证集,结果如表3所示。将各总集、建模集、验证集样本点土壤含盐量划分为4个等级:非盐土(D1,SSC(土壤含盐量)小于等于0.2%)、轻度盐渍化(D2,SSC为(0.2%、0.5%])、重度盐渍化(D3,SSC为(0.5%、1%])和盐土(D4,SSC大于1.0%)[34]。非盐土、轻度盐渍化、重度盐渍化和盐渍土占比分别为34.1%、55.7%、9.1%和1.1%。含盐量变异系数均处于中等差异(变异系数CV反映样点值的离散程度,CV<0.1为弱变异性;0.11为强变异性)[34]。

表3 土壤含盐量特征统计分析Tab.3 Statistical analysis of soil salt content characteristics

对各盐渍化等级下建模集、验证集和总集的含盐量进行统计分析,如图3(图中1.5IOR表示1.5倍的四分位距)所示。从图3可以看出,D1、D2、D3等级的建模集、验证集和总集的含盐量分布、值域和均值相近,确保建模集和验证集数据的代表性。

图3 土壤含盐量特征统计Fig.3 Statistical map of soil salinity characteristics

2.2 不同植被指数分类结果

基于近红外和可见光波段构建4个植被指数,利用对土壤背景敏感的4个植被指数(NDVI、MSAVI、DVI、CRSI)进行图像分类,验证不同指数的图像分类精度。为了对比4个植被指数图像分类的精度,从87个采样点中随机选取40个样本,进行如图4所示,4个植被指数的精度存在显著差异,其中NDVI的图像分类精度最高,40个样本的总体精度均大于93.4%,均值为96.9%;Kappa系数在0.84～0.99之间,均值为0.92。MSAVI的分类精度略低于NDVI,总体精度和Kappa系数均值分别为95.8%和0.91;DVI和CRSI 2个植被指数的分类效果相对较差,总体精度均值小于91%,Kappa系数均值均小于0.88。基于以上分析,NDVI图像分类结果表现最优。

图4 基于不同植被指数的图像分类结果Fig.4 Evaluation based on image classification results of different vegetation indices

分类结果评价。分别对40个样本的多光谱图像进行2.1节的操作,获得分类结果,并在ENVI 5.3中对分类结果进行评价,得到4个植被指数的总体精度和Kappa系数。

2.3 图像纹理计算

灰度共生矩阵具有丰富的特征参数,能从不同的角度对纹理进行细致刻画[35]。由盐渍化土壤覆盖植被冠层的光谱特征分析可知,植被冠层光谱中红光波段与土壤含盐量的相关性最高,对盐渍化的变化有最显著的响应,故选择波段Band3计算多光谱图像纹理。为了更好地体现植株间冠层图像的纹理差异,纹理特征参数计算的窗口选择尺度为11像素×11像素,尺寸为71.5 cm;每个特征参数都有4个不同方向的值,取其平均值作为方向无关的特征值。将采样点的多光谱图像输入ENVI 5.3软件,利用二阶统计滤波工具计算2个波段灰度共生矩阵的8个特征参数。如图5、6所示,分别为剔除土壤背景的植被冠层图像(E1)和原始多光谱图像(E2)处理的可见光示意图,均值(MEA)、方差(VAR)、均匀性(HOM)、对比度(CON)、差异(DIS)、熵(ENT)、二阶矩(SEC)和相关性(COR)特征参数灰度图。从图中可看出,灰度共生矩阵的特征参数表现出图像丰富的纹理特征,剔除土壤背景前后的多光谱图像纹理特征存在明显差异,统计灰度图的纹理特征参数数值,对纹理特征进行进一步分析。

图5 剔除土壤背景多光谱图像灰度共生矩阵特征参数灰度图Fig.5 Gray scale of characteristic parameters of co-occurrence matrix of multispectral image with soil background eliminated

图6 原始多光谱图像灰度共生矩阵特征参数灰度图Fig.6 Gray level co-occurrence matrix characteristic parameters gray level

2.4 图像纹理特征

图像纹理是用于识别图像中感兴趣的对象或区域的重要特征之一[36],为了对比剔除土壤背景以及不同盐渍化等级下多光谱图像纹理特征的差异,对E1和E2处理的多光谱图像Band3的纹理特征参数进行统计,计算2个处理下3个盐渍化等级样本点多光谱图像的纹理特征参数的平均值,如表4所示。

表4 多光谱图像纹理特征参数Tab.4 Multispectral image texture feature parameters

分别对比2个处理下3个盐渍化等级的纹理特征参数平均值。整体上,相比剔除土壤背景前的E2,剔除土壤背景后的E1处理的纹理特征值中,SEC、CON、MEA特征值减小,DIS、VAR、ENT、COR、HOM特征值增大。其中CON、ENT、COR 3个参数的变化最显著,CON减小表示图像纹理沟纹变浅,图像清晰度下降;ENT变大,表示图像纹理增多;COR增大,表示图像纹理变粗。

2.5 基于全子集回归筛选光谱自变量

通过R语言全子集回归模型算法,分别对E1和E2处理的 a (植被指数)和b (植被指数和图像纹理特征参数) 数据集进行全子集变量筛选,筛选评价指标选择R2和 BIC 两个指标。综合衡量R2最大和 BIC最小的变量组合,得到每个数据集中最优的变量组,结果如表5所示。

表5 全子集筛选最优变量组合统计Tab.5 Optimal variable combination statistics screened by full subset method

对比E1处理下2个变量种类的筛选结果,变量组a筛选得到的敏感变量组含有NDVI、SRVI、NDSI 3个植被指数;变量组b的敏感变量组中包含NDVI、

SRVI、ARVI、TVI 4个植被指数,筛选掉了纹理特征参数。对比E2处理下2个变量种类的筛选结果,变量组a筛选得到的敏感变量组包含DVI、EVI、SRVI 3个植被指数;变量组b敏感变量组中包含RVI、NDSI两个植被指数和ENT、COR、SEC 3个图像纹理特征参数。

2.6 基于图像纹理特征和植被指数的SSC反演模型

通过不同植被指数分类结果分析可以得出NDVI可以达到良好的图像分类效果。因此,利用NDVI对多光谱图像进行分类,并对图像分类得到的土壤像元进行掩膜处理,剔除土壤像元,获得纯净度更高的植被冠层图像。分别将E1和E2输入ENVI 5.3软件提取6个波段的反射率。将多光谱可见光和近红外1波段的反射率代入植被计算公式(表1),得到相应的植被指数,涉及的植被指数有NDVI、DVI、RVI、MSAVI、ARVI、EVI、CRSI。并基于全子集回归变量模型筛选结果,构建基于敏感变量组的SSC反演模型。

基于E1和E2处理下a、b两种变量组全子集法筛选的变量结果,构建基于敏感变量的SVM和ELM土壤含盐量反演模型,各模型建模集和验证集的R2和NRMSE如表6所示。

表6 基于不同变量组的SVM、ELM反演模型精度Tab.6 SVM and ELM inversion model based on different variable groups

基于全子集筛选的2种敏感变量组,利用SVM和ELM两种机器学习方法,构建E1和E2处理共8个土壤含盐量反演模型,验证集实测值和模拟值如图7所示。

图7 土壤含盐量实测值和模拟值对比Fig.7 Comparison of measured and simulated soil salt contents

E1、E2处理下2个变量组的SVM和ELM机器学习土壤含盐量反演模型中,ELM算法的E1-a、E1-b、E2-a、E2-b模型的R2分别为0.602、0.657、0.584、0.681,均大于对应的SVM模型。ELM模型的标准均方根误差分别为0.181、0.156、0.201、0.191,均小于SVM模型验证集的标准均方根误差。综上,ELM模型的表现更好。

3 讨论

表7 剔除土壤背景前后纹理特征参数平均值Tab.7 Average value of image texture characteristic parameters before and after removing soil background

本研究仅基于无人机遥感多光谱反射率、图像纹理特征与植被指数的土壤含盐量反演,对于作物的耕种方式、作物不同生育期受土壤盐渍化的影响,灌水量对土壤含盐量的影响尚未考虑。本研究建立了特定情形下的土壤含盐量反演模型,模型结果较好,但模型是否具有广泛的适用性,不仅需要更加深入的研究,还需要综合考虑更多因素,进一步提高反演模型精度。

4 结论