基于Sentinel-2影像和机器学习算法的冬小麦秸秆覆盖度遥感估算

朱奇磊,梁 栋,徐新刚,安晓飞,陈立平,杨贵军,黄林生,许思喆

(1.安徽大学电子信息工程学院,安徽合肥 230601;2.北京市农林科学院信息技术研究中心,北京 100097;3.北京市农林科学院智能装备技术研究中心,北京 100097;4.安徽大学 农业生态大数据分析与应用技术国家地方联合工程研究中心,安徽合肥 230601)

作物秸秆是指收获后残留在农田地表的作物剩余物质。作物秸秆覆盖于农田土壤表面,具有减少土壤水分蒸发和减缓土壤风化的作用[1]。同时,秸秆还田也有助于提高农田土壤肥力,是保护性耕作的重要内容[2-3]。作物秸秆覆盖度(crop residue coverage,CRC)指单位面积田块中秸秆覆盖所占的比例,是表征作物秸秆数量与分布的重要参数。实现大面积、实时和准确的作物秸秆覆盖度监测,对评估农田生态环境和开展保护性耕作具有重要意义。

目前,作物秸秆覆盖度测取方法有拉绳法和照相法。拉绳法在实际操作过程耗时费力,工作效率低[4];照相法具有便携和易操作的优势[5-8],但对于不同区域农田秸秆覆盖度的监测,仍不能满足大范围、实时监测的需求。遥感技术具有大面积、快速和动态监测的特点,在作物秸秆覆盖度估算中扮演着越来越重要的角色。李志婷等利用地面高光谱数据,构建不同遥感影像的波段反射率,通过构建敏感光谱特征指数用于小麦秸秆覆盖度估测的分析比较,结果表明,Landsat-8 OLI数据在波段划分上更为精细,且利用其构建的光谱指数NDIOLI21估测精度优于Landsat-5 TM[9]。Cai等综合利用Sentinel-2光学与Sentinel-1微波数据,构建光谱与微波信息的复合指数对作物覆盖度的估测精度较线性模型得到了提升[10]。Daughtry等在复杂混合场景下获取的小麦、玉米、土壤等地物高光谱数据进行作CRC分析试验,结果发现,作物秸秆在2 100 nm处所特有的光谱吸收特征能够很好地将秸秆与土壤区分开来,进而构建纤维素吸收指数(cellulose absorption index,CAI)[11]。黄晋宇等运用哨兵Sentinel-2卫星遥感数据,结合典型的光谱指数与不同土壤类型因素实现玉米秸秆覆盖度的估测[12-13]。在众多CRC估测中运用光学特征居多。Jin等利用Landsat-8 OLI遥影像提取的植被指数和纹理特征估测玉米秸秆覆盖度,结果显示,影像纹理特征对玉米秸秆覆盖度进行估测具有可行性[14]。纹理特征在遥感数据中能提供丰富的空间信息。Li等运用机器学习算法结合GF-1的光谱特征与纹理特征与Sentinel-1微波特征实现森林生物量的估测,结果表明,利用GF影像提取的光谱与纹理特征并结合随机森林模型在森林生物量估测中精度最高[15]。

目前,在光谱植被特征信息的基础上加入纹理特征,在玉米秸秆覆盖度的遥感监测中是可行的,但小麦和玉米秸秆在农田地表结构和形态上具有较大的差异,在光谱植被特征基础上引入纹理信息是否能显著提高小麦秸秆覆盖的监测精度需要进一步分析探究。此外,随着遥感特征变量的增多,不同数据之间的特征融合,如何有效筛选对作物秸秆敏感的特征量,并尽可能提高估测精度是当前作物秸秆覆盖度遥感监测研究的重要内容。

本研究以当前广泛使用的Sentinel-2卫星影像为数据源,探讨利用典型光谱特征、植被指数和纹理特征及其组合在估测小麦秸秆覆盖度的可行性,通过灰色关联度与随机森林(GRA-RF)相结合的变量筛选方法,优选对小麦秸秆覆盖度敏感的特征变量,利用套索(LASSO)、岭回归(RR)、高斯过程回归(GPR)等机器学习算法,针对不同特征变量集,构建小麦秸秆覆盖度遥感估算模型并实现小麦CRC遥感估测,以期为大范围小麦秸秆覆盖度估测提供新的思路和方法。

1 材料与方法

1.1 研究区概况

试验于2020年7月在北京市小汤山国家精准农业示范基地进行 (40°00′N～40°21′N,116°34′E～ 117°00′E,见图1)。研究区地势平坦,平均海拔高度约为36 m,四季分明,属于大陆性季风气候,春季干旱多风,夏季炎热多雨,年平均日照时长2 684 h,降水量500～600 mm。图1为研究区位置及Sentinel-2卫星影像采样点分布。

图1 小麦试验区域位置及样点分布

1.2 数据获取及预处理

1.2.1 田间实测小麦CRC获取

于2020年7月7日,采用手持多光谱相机(具有红、绿、蓝、近红、红边五个光谱波段)在固定高度对试验采样区域地面正上方垂直拍摄,共获取115个高分辨率秸秆覆盖数据及拍摄时间和经纬度信息。为了更接近秸秆还田真实场景,对试验区秸秆在自然环境中静置处理,以保证在下季作物种植期有更真实的秸秆覆盖数据。图2为原始多光谱数据真彩色合成图。

图2 田间地面实测图像

对田间试验样点采用人工目视解译的方法选取5个地物类别构建ROI感兴趣区域[土壤、秸秆、绿植、阴影(土壤阴影、非土壤阴影)]。目前常见的地物分类算法有支持向量机、随机森林、最小距离等,本研究通过试验分析选取分类精度较高的随机森林分类算法进行影像处理。地物分类后采用混淆矩阵与Kappa系数作为评价标准,其中RF的两种系数评估结果分别为88.33%、0.85。

1.2.2 遥感数据处理

借助Google Earth Engine(GEE)平台获取Sentinel-2在2020年7月6日本试验区域的遥感数据,并对遥感影像做相应预处理,包括辐射定标、大气校正、几何校正处理,遥感数据状态良好且无云[16-17]。

1.3 特征参数提取

1.3.1 遥感影像光谱特征提取

针对近地面高分辨率的多光谱数据,通过变量筛选的方法选取重要性程度高的光谱指数,对地面样点进行分类处理,以减少运算量和冗余信息,提高处理效率。以往的研究基于Sentinel-2数据利用短波区域合成的光谱指数对遥感影像处理,其中枯燃料指数(dead fuel index,DFI)、归一化差异耕作指数(normalized differential tillage index,NDTI)、简单耕作指数(simple tillage index)以及归一化差值指数NDI(x,y)等[18-19]对秸秆覆盖度估测效果较好。

DFI=100(1-SWIR2/SWIR1)×R/NIR

(1)

NDTI=(SWIR1-SWIR2)/(SWIR1+SWIR2)

(2)

STI=SWIR1/SWIR2

(3)

DNI(X,Y)=(Y-X)/(Y+X)

(4)

式中SWIR1、SWIR2分别为波段B11、B12,R、NIR分别为B4和B8,x,y为遥感影像的不同光谱波段。

1.3.2 遥感影像纹理特征提取

对Sentinel-2影像从不同波段分别获取灰度共生矩阵(gray-level co-occurrence matrix,GLCM)的8种纹理特征(表1),包括均值(mean)、方差(variance ,Var)、同质性(homogeneity,Hom)、对比度(contrast,Con)、差异性(dissimilarity,Dis)、熵(entropy,Ent)、二阶矩(second moment,Sec)、相关性(correlation,Cor)。随着这些纹理特征的加入,估测模型就可以反映空间特征。

表1 纹理特征表达式

1.4 数据分析方法

1.4.1 灰色关联-随机森变量筛选方法

灰色关联分析(GRA)是一种能够权衡特征变量与比较数列关联程度的无量纲统计分析方法,可通过计算变量之间的关联度,获得变量的排序情况,进而分析自变量与因变量之间的关系。关联度较大的参量表明其与比较数列的变化趋势更为相近[20-25]。随机森林(random forest,RF)算法属于一种集成算法,由许多决策树构成的分支集合而成,众多决策树的投票结果能够共同完成对数据属性的分类与预测任务,通过对每个特征在决策树上的平均贡献率进行评估,获得特征的重要性排序,同时所采用Bootstrap抽样法能有放回地随机生成多个决策树数据,而且随机性的引入使得模型不容易陷入过拟合[26-28]。灰色关联分析能对多目标灵活处理数据,但仅适用于线性关系的分析,随着数据信息量的增加,其分析过程也将变得复杂。而随机森林作为机器学习方法,既能实现多特征的非线性分析,也能应对大量的数据信息,在模型构建过程中对不同特征信息进行充分学习,并提供特征重要性排序。因此,本研究尝试对两者优势进行融合。首先运用两种分析方法各自运算获得的结果分别进行归一化处理,然后将特征重要性得分对应加和,获得特征信息排名。

(5)

(6)

式中,x0(j)为参考因子数列;xi(j)为比较因子数列;α为分辨系数,通常取0.5;i为比较因子的个数;j为数据的组数;用公式5计算关联系数后求取平均值,获得关联度ri。

1.4.2 机器学习算法

本研究应用偏最小二乘、岭回归、套索回归及高斯过程回归等代表性的4种机器算法,比较分析不同算法在冬小麦秸秆覆盖度遥感监测中的性能和作用。

偏最小二乘回归(Partial least squares regression,PLSR)属于数学统计的一种优化方法,能够用于构建多个因变量与多个自变量之间的模拟关系,结合试验参数进行主成分分析,突出对模型贡献较大的主成分信息,通过对多变量间的相关性检测,减少多重共线性下的冗余参量,最后进行多元线性回归分析[29-31]。

岭回归(Ridge regression,RR)是一种能够消除多重共线性所带来的回归误差的方法,是对最小二乘的一种改进算法,最先用来处理样本量少于输入特征参数的情况,通过对目标函数引入正则项的方式,使得无解的状况得到解决,减少了参数量,使得模型更稳定且更符合实际[24,32]。

套索回归(Least absolute shrinkage and selection operator,LASSO)由Robert Tibshirani在1996年首次提出。该算法集合了RR与最小二乘的优势,能在建模时自动将不重要回归系数自动剔除,降低了模型的复杂性即减少总体的回归系数,能够实现自动筛选变量的功能,在一定程度上降低共线性的影响[31,33-34]。

高斯过程回归(Gaussian process regression,GPR)是利用高斯过程对数据进行分析来寻找数据之间内在关联,在处理非线性问题上具有较好的表现,也适用于解决样本量小的回归问题,具有较好的泛化能力与适应性[35-38]。

1.5 精度评估

本研究在多种特征融合基础上构建机器学习回归模型,选取70%数据为训练集,30%数据为验证集,用10折交叉验证获取模型平均得分作为最终评估结果,采用决定系数(coefficient of determination,r2)和均方根误差(root mean square error,RMSE)来分析模型精度。其中,r2越接近于1,RMSE越接近于0,模型拟合效果越好。用方差膨胀因子(variance inflation factor,VIF)检测模型变量参数之间是否存在共线性问题,当VIF>10时,则认为变量之间存在共线性问题。

(7)

(8)

(9)

2 结果与分析

2.1 特征变量相关性与GRA-RF的敏感性分析

2.1.1 基于不同数据的特征相关性分析

对光谱波段、光谱指数、纹理特征与CRC进行相关性分析。从相关性由大到小的排序来看,排在前五的遥感光谱波段、光谱指数与纹理特征排序分别依次为B12>B4>B3>B8>B5、DFI>NDTI>STI>NDI(4,11)>NDI(3,11)和Mean_B12>Hom_B3>Sec_B3>Dis_B3>Ent_B3(图3)。由于不同特征变量之间具有一定的相关性,因而在多元线性回归模型中需要利用方差膨胀因子(VIF)检测来消除变量之间的多重共线性问题。基于这些特征变量与CRC线性回归建模,结果(表2)表明,光谱指数DFI、STI和NDTI的建模效果较好。

表2 CRC线性模型分析

图3 波段、指数、纹理特征相关系数分析

2.1.2 基于GRA-RF的敏感特征分析

GRA-RF法分别对光谱、指数、纹理特征进行敏感特征筛选,每种特征关联度较大的前五个特征排序结果与相关性分析结果具有一定的相似性(图4),即关联度最高的特征相同。

图4 不同特征的GRA-RF特征重要性排序

2.1.3 基于机器学习的CRC遥感估测

为探究不同数据特征及其相互组合对CRC模型精度的影响,根据重要性排序结果,结合多种机器学习算法来构建CRC估测模型。结果发现,多个波段特征组合的模型精度较线性模型明显提升,其中较稳定的模型有LASSO、RR、PLSR模型(表3)。在用光谱指数构建的回归模型中RR模型较为稳定且与DFI模型精度相近,光谱指数特征及其组合特征在所选机器学习模型中r2均稳定在0.5以上。纹理特征作为GPR模型特征参数时获得最稳定且精度最高的CRC估测结果。多种特征的融合在不同模型中各有突出的表现。与线性回归模型相比,运用机器学习算法与特征组合方式使得CRC模型总体估算精度得到了 改善。

表3 多种特征组合回归模型精度评估

2.2 基于不同光谱特征变量数的CRC遥感估算趋势分析

为进一步探究特征数量对模型精度影响,基于以上分析结果,通过逐步增加特征数量来观察不同特征及其组合在机器学习建模过程中的精度变化趋势。其中,对光谱波段、光谱指数以及纹理特征通过关联分析选取前10个排序结果,不同特征之间组合各选取前5个作为模型训练参数,其中模型参数最多的为波段、纹理、光谱三者组合,训练时最多含有15个变量且初始变量随机。

通过建模分析绘制CRC模型趋势图(图5),对其进行分析发现,含有光谱指数的所有模型的初始精度较高且整体较为稳定,对CRC具有较好的估算效果。而含有波段特征的精度曲线中,回归模型随着特征数的增加,其精度首先迅速提升后趋于平稳,其中在模型趋于稳定前达到最高的两个特征波段为B12与B11,说明机器学习算法结合敏感变量筛选构建的多元模型可实现小麦秸秆覆盖度较好估测。纹理特征在前5个特征中估测精度相对不高,但在GPR模型中随着特征数增加,模型最终趋于稳定且精度优于其他模型。LASSO与RR模型总体趋势较为相似,仅在少数变量变化中出现较小的上下浮动,而LASSO是基于RR回归的一种改进算法,随着变量数的增加,能自动剔除对模型精度贡献不大的特征。综合分析后,从两种算法中可选取精度较好且相对稳定的LASSO构建最优模型。表4为部分特征与LASSO算法建模参数及精度,所有特征变量均通过方差膨胀因子检验,模型各参数的VIF均小于10。图6对应表4中CRC模型在训练集与测试集的精度评估结果。依据三种模型分析得到研究区CRC的分布状况(图7),与实际观察较为相符。

图5 CRC模型精度趋势图(r2为三角,RMSE为圆点)

a:波段特征;b:波段与纹理;c:光谱指数、波段与纹理特征。

图7 不同特征组合下LASSO模型估测CRC分布

表4 不同特征与LASSO模型精度及公式

3 讨 论

本研究利用不同光谱特征及其组合构建机器学习模型,实现小麦秸秆覆盖度的估测,依据训练集与验证集的决定系数差异最小的原则,筛选稳定性最优的CRC估测模型。从相关性分析结果看,在遥感影像波段特征中B12与小麦CRC敏感性,高于其他光谱波段,主要是因为作物秸秆在2 100 nm处光谱反射率与土壤有较大的区别[11];光谱指数与CRC的相关性优于波段特征和纹理特征,其中利用DFI、STI、NDI(x,y)构建的CRC估算模型具有较好的精度,与已有研究结果相近[18-19]。当前,通过纹理特征实现作物秸秆估测研究较少,与玉米秸秆的纹理特征相比,小麦秸秆更为细长密集且反射面积小,空间分辨率越高,纹理特征可能会越明显。在纹理特征中Mean_B12与CRC的相关性表现较好,但构建多特征模型时,其他纹理特征对CRC精度的提升优于Mean_B12,这可能是B12的光谱特征或指数特征与其本身的纹理信息存在一定的共线性。在以往的多元线性回归模型中,随着变量特征数的增加,模型精度会得到一定的提升,但不同变量特征之间可能存在一定的共线性问题。因此,本研究利用敏感特征筛选和模型趋势分析,探究机器学习在CRC

估算中最优的模型方法及最佳特征变量。

本研究通过敏感特征筛选,利用排序后的特征进行建模,并分析特征变量的变化对CRC估算精度的影响,结果表明,随着光谱变量数的增加,机器学习算法模型如PLSR模型对CRC的估测精度没有实质性的提升,这与前人研究结果相似[15]。LASSO等算法依靠自身具有的正则化约束作用,在一定程度上消除了共线性问题[39]。从特征信息的增加对CRC模型精度的影响看,随着变量数的增加,利用机器学习算法构建的估算模型的精度有上升或下降,表明不同特征的加入对模型估算精度的影响不同。通过趋势线分析,经GPA-RF筛选的5个特征在单特征建模中,模型精度不随变量的增加而增加,且部分特征组合在最高精度具有相一定的相似性。相比仅利用波段特征信息,波段与纹理的组合模型中,随着纹理信息的加入,GPR、PLSR模型得到稳定,而LASSO与RR模型精度在稳定后得到提升,在光谱指数模型中加入纹理也有相似的表现。在多种特征组合情况,RR与LASSO模型精度趋于稳定甚至整体精度得到提升,而GPR与PLSR模型在个别特征组合中最高精度并无明显差异,说明建模信息过度冗余,后续需要通过对曲线精度最高点的特征进行提取,在原有基础上构建新的CRC估测模型以降低模型的复杂性。

本研究通过GRA-RF方法筛选特征变量,再结合多种机器学习算法构建CRC估测模型,相比传统光谱指数模型,混合光谱特征模型能显著提高整体精度,为实现大面积小麦CRC估测提供新的研究思路。但本研究仅基于秸秆还田试验基地,且部分田块采用轮作的种植方式,故从采集的图片数据中可以看出,在静置处理过程中出现新一季作物幼苗,这在秸秆还田试验中需要考虑,在后续的研究中可以结合作物物候期与多时相遥感数据对研究区秸秆还田时期进行判定。此外,本研究区局部地形的变化对模型影响较小,在大区域尺度进行CRC研究时需考虑不同地区在土壤含水量上的时空分布差异。同时,不同遥感影像在空间分辨率上有所区别,而较高分辨率影像包含更多的纹理信息。相比遥感光谱信息,纹理特征能够反映影像数据的表观性质,在CRC估算中能反映秸秆与土壤纹理上的空间差异,弥补了仅利用单一光谱特征对CRC估算的局限性。随着国产遥感卫星技术的不断发展,利用较高的空间分辨率与光谱分辨率优势实现作物秸秆覆盖度估算具有重要研究意义。综上所述,在作物秸秆实地研究中,田间场景的复杂程度、地表物质的多样性以及土壤的含水量与有机质含量等需要进行采样分析。在后续研究中,将考虑区域范围更广、秸秆类型更多、分布更为均匀的野外覆盖度估测分析,并为不同区域场景构建相应的地物特征库,进而构建具有普适意义的秸秆覆盖度模型。

4 结 论

基于Sentinel-2遥感数据,通过GRA-RF方法对波段、光谱指数、纹理特征进行敏感性分析,将不同特征信息的组合并结合多种机器学习算法构建冬小麦秸秆覆盖度估测模型。结果表明,基于单一特征的模型中,利用光谱指数的估测精度较高,其中DFI表现最好(r2=0.54,RMSE= 10.26%);单一纹理特征对冬小麦CRC的敏感性较差,但利用多个纹理信息构建机器学习估算模型时精度得到改善(r2=0.43～0.46)。通过构建特征变量与精度趋势变化曲线,对变量变化拐点处的特征进行提取,筛选出小麦CRC最优模型特征,其中3种特征(B11,B8,B4_B11,STI,Ent_B3,Hom_B3,Mean_B2)混合的模型精度最优,r2与RMSE分别为0.65%和9.25%。因此,与单特征变量构建的模型相比,多特征信息的组合显著提高了小麦CRC的估算精度。