基于MODIS时序数据物候特征下的多源遥感玉米提取

2023-06-26刘惠楠王井利马运涛

江西农业学报 2023年4期

刘惠楠，王井利，周斌，马运涛

（1.沈阳建筑大学交通与测绘工程学院，辽宁沈阳 110168；2.辽宁省生态气象和卫星遥感中心，辽宁沈阳 110166）

玉米作为我国主要的粮食作物之一，准确掌握玉米的分布情况有助于推动我国精准农业的发展。目前，关于玉米种植结构的识别及提取的研究主要基于中高分辨率遥感影像的地块尺度方面［1-5］。中高分辨率的遥感影像虽然有较高的分类精度，但由于受到重访周期、图幅宽窄、云污染等影响，监测较大尺度农业种植面积的中高空间分辨率遥感数据难以满足精确研究农作物关键生育期识别的要求。

在实际应用中，玉米属于季节性作物，准确区分其物候期，可为实现玉米种植结构的准确提取提供重要依据［6］。其中，利用MODIS制作的时序数据最为广泛，它具有观测周期短、范围广等优点［7-8］。但在市(县)尺度的研究中，MODIS数据的空间分辨率较低，从而影响了对作物识别的精度，因此，引入较高空间分辨率的遥感影像联合应用可互相弥补不足之处，进而提高对作物识别的精度与效率［9-13］。

大量研究表明，多源数据联合应用的分类方法已成为遥感农作物识别的主要发展趋势。其中，MODIS数据的空间尺度较大，与其他高分辨率影像数据融合后极易造成混合像元，导致融合后影像数据的应用受限，对作物提取结果的精度并不理想。本文提出了不基于影像融合的多源数据联合应用提取作物种植结构的方法，为提高农作物遥感种植结构提取的准确度提供了一种新的思路。

1 资料与方法

1.1 研究区概况

沈阳市位于辽宁省中部，位于41°10′30″ N~43°0′42″N、122°23′22″E~123°46′38″E 之间，处于辽东山地与下辽河平原的过渡地带，属温带季风气候区湿润和半湿润大陆性气候，四季分明，降水集中，日照充足，雨热同季，年均降水量628 mm，适于农作物的生长发育。研究区地理位置和Landsat-8真彩色影像如图1所示。

图1 研究区位置示意图和Landsat-8真彩色影像

1.2 研究数据

1.2.1 多源遥感数据的获取与预处理本文应用Landsat-8、Sentinel-2、GF-1、MODIS等4类遥感影像，选取了沈阳区域2020年无云或云量较少的中高分辨率卫星遥感影像共28幅，2020年5月1日—10月1日低分辨率卫星遥感影像共154幅，多源遥感影像具体数据见表1。

表1 2020年5月1日—10月1日的多源遥感影像参数

本文所用Landsat-8 OLI多光谱影像数据来源于美国地质勘查局网站(https://glovis.usgs.gov/)，Sentinel-2A/B MSL多光谱影像数据来源于欧空局哥白尼数据中心(https://scihub.copernicus.eu/)，GF-1 WFV数据来源于辽宁省高分数据中心，MODIS数据来源于NASA网站(http://modis.gsfc.nasa.gov/)提供的MODIS1B陆地带MOD09(MODIS表面反射率)的MOD09GA数据产品。

遥感数据均经过辐射定标、大气校正、投影转换、波段合成、正射校正、图像配准、镶嵌、裁剪及光谱指数计算而来。各光谱指数计算公式为:

增强型植被指数：

归一化水体指数：

归一化建筑指数：

1.2.2 样本的建立与数据的验证基于现有的东北地区10 m作物类型分布产品数据集，制作玉米(157处)、水稻(100处)、大豆(37处)的样本区域数据。计算Landsat-8、Sentinel-2、GF-1、MODIS遥感影像中玉米157处样本区域的EVI平均值，并以此作为转化模型的样本数据。结合实地调查(图1)与研究区Google Earth影像制作玉米(86处)、水稻(50处)、大豆(60处)、林地(50处)、水体(37处)、其他植被(40处)、非植被与水体区域(100处)为结果验证数据。

1.3 研究方法

1.3.1 Savitzky-Golay滤波 S-G滤波是一种在域内基于多项式，通过移动窗口利用最小二乘法最佳拟合的方法［14］。具体计算公式为：

式(4)中，Yj+i为原始EVI的时序数据，′为EVI的拟合数据，Ci为第i个EVI值的滤波系数，N为卷积数目，等于组数的宽度(2m+1)，j为原始EVI数据的系数，m是平滑窗口大小的一半。

该研究利用2020年5月1日—10月1日共计154 d的每日MOD09GA数据中玉米、水稻、大豆各处样本区域合计的日平均EVI值，制作作物完整物候期的时间序列。由于EVI数据受云的影响明显，导致部分日期数据存在明显偏差。利用S-G滤波，滑动窗口设置为5，使得保留细节特征的同时更接近EVI时序的上包络线，对EVI数据的噪声进行有效剔除。

1.3.2 各传感器间EVI指数转化模型的构建利用中高分辨率EVI数据来弥补MODIS-EVI时间序列分辨率较低的问题。选取与中高分辨率影像对应同一天的MODIS数据，可以较大程度地避免卫星数据获取过程中由地物变化导致影像辐射值变化的问题，影像预处理可消除或尽量减小数据获取时太阳高度角、气溶胶等因素的影响。在尽可能地减小由数据来源差异引起的误差后，分别计算MODIS与Landsat-8、Sentinel-2、GF-1各个影像157处样本区域的EVI平均值，使之一一对应构建数据对，利用线性回归方法，以中高分辨率EVI数据作为自变量、对应的MODISEVI作为因变量构建EVI指数转化模型，并利用决定系数R2评价模型的拟合效果。

1.3.3 分类方法 CART决策树、随机森林、最大似然法为传统监督分类方法，利用已知分类样本像元对未知像元进行判别。该研究利用ENVI软件，输入待分类数与样本数据后计算机自动输出分类结果。

决策树分类技术在遥感分类方面具有很大的优势，可逐步建立决策树，手动设置阈值，将目标地物逐步从原始数据中分层掩膜出来，能有效地避免在提取过程中相似地物的相互干扰，通过共同运用大量不同类别的输入数据，从而达到多种类数据复合应用分类的目的，且分类结构清晰并容易理解［15-16］。

1.3.4 混合像元分解模型玉米常与其他作物间作，在遥感解译中极易产生混合像元，导致地物分类精度下降。该研究选择Heinz等［17-18］提出的全约束最小二乘法混合像元分解模型(Fully Constrained Least Squares,FCLS)，在决策树分类结果的基础上进行混合像元分解，进而提高对数据的解译精度。

1.3.5 精度评价方法该研究利用R2评价模型的拟合优度，以及标准差、均方根误差评价经模型转换前、后各个样本点的离散程度。利用验证数据的混淆矩阵，计算各地类的分类精度，使用Kappa系数和总体精度评价整体分类结果的精度，利用各地类的制图精度与用户的精度评价各地类的漏分与错分情况。

2 结果与分析

2.1 构建时间序列与分析物候特征

S-G滤波重构后的各类农作物时序曲线如图2所示，图2中2020年5月1日对应时间为0 d，2020年10月1日对应时间为153 d。

图2 作物全生长周期EVI时序曲线

玉米EVI在5月中下旬开始增加，此时玉米开始出苗生长；7月上旬的EVI增长最快，此时玉米处于拔节期，其生长态势最旺盛；7月14日(第75天)达到峰值，为0.660，此时由于玉米顶端雄穗与玉米植株颜色的反差导致各通道反射率发生变化，使得EVI开始下降；8月8日(第100天)到达波谷，为0.512，此时玉米处于吐丝期；随后植株继续生长，EVI经过小幅度增长后，玉米进入乳熟期；9月EVI迅速下降，玉米则逐渐成熟。

水稻于5月初出苗；EVI在5月中下旬缓慢减小，稻苗叶片收缩，此时水稻处于返青期；6月中下旬EVI迅速增长，此时水稻处于分蘖期；7月17日(第78天)EVI达到峰值，为0.683，此时水稻处于拔节期；随后EVI开始逐渐下降，水稻进入抽穗期，虽然此时稻穗颜色发生了变化，并由此导致各通道的反射率也随之发生了变化，但稻穗颜色的变化比玉米雄穗与植株的少，因此，此时的时序曲线未产生波谷；9月EVI迅速下降，水稻进入逐渐成熟。

大豆EVI在5月中旬开始增加，此时大豆开始出苗生长；6月中下旬EVI增长速率最快，生长态势最旺盛，此时大豆处于分枝期；7月20日(第81天)EVI为0.571时开始缓慢增长，大豆进入开花期；8月20日(第112天)达到峰值0.583，此时大豆处于结荚期；9月EVI迅速下降，大豆进入鼓粒期后逐渐成熟。由上述分析获取的作物物候期如表2所示。

表2 沈阳市主要作物物候期

2.2 各传感器间EVI指数转化

各传感器间EVI指数的转化模型及相关系数(R2)如图3所示。各转化模型的R2＞0.85，各转化模型的拟合优度较高。

图3 EVI指数转换模型

模型转换前，MODISEVI的数据相较于其他类型的标准差更小，离散程度较小，数据较稳定。经模型转换后，GF-1EVI、Landsat-8EVI与Sentinel-2EVI的标准差与均方根误差均有一定程度的减小，EVI离散范围明显缩小，具体标准差与均方根误差如表3所示。

表3 模型转换前、后各类EVI样本的标准差与均方根误差比较

转换后的数据能够与同日期的MODISEVI数据相匹配，具体如图4所示。将研究区域的Landsat-8、Sentinel-2、GF-1的EVI影像数据进行对应模型转换后，得到了各个日期的MODISLandsat-8EVI、MODISSentinel-2EVI、MODISGF-1EVI。

图4 模型转换前、后数据的对比图

2.3 作物种植结构的提取与精度评价

2.3.1 传统分类方法的玉米种植结构提取采用CART决策树、随机森林、最大似然法分别对中高分辨率影像进行地物提取，各类方法的总体精度、Kappa系数如表4所示。总体分类精度均小于86.5%，Kappa 系数均小于0.8，由此可知，这3种分类方法结果的精度并不理想。

表4 传统遥感分类方法的结果精度比较

2.3.2 基于决策树的多源遥感数据种植结构提取经反复实验发现，EVI在水体和建筑用地的分类中会将部分水体与建筑用地混淆，故利用9月30日的Sentinel-2 NDWI数据，将NDWI值＞0.5的区域判定为水体。利用7月22日的Landsat-8 NDBI数据，将NDBI值＞-0.1236且NDWI值＜0.5的区域判定为非植被与水体区域。利用7月22日的Landsat-8 EVI数据，将EVI值＜0.5755、NDWI值＜0.5、NDBI值＜-0.1236的区域判定为其他植被。

利用7月15日的GF-1影像、8月7日的Landsat-8 EVI数据、9月20日的Sentinel-2 EVI数据，根据作物的物候期将植被细分为林地、玉米、水稻、大豆。图2是由MODIS-EVI制作，故其只能作为各地类间植被指数大小比较的参考，其具体数值不能作为未经模型转换的数据分类阈值使用。基于此，经反复实验确定，由于8月7日的水稻处于抽穗期，此时的EVI值大于玉米和大豆的，故将8月7日的Landsat-8 EVI值在1.6213~1.8525范围内的区域判定为水稻；9月20日的大豆处于鼓粒期，此时的EVI值大于玉米的，故将不满足上述条件且9月20日的Sentinel-2 EVI值在-0.1021~0.0662范围内的区域判定为大豆。7月15日的玉米与水稻的EVI值远大于大豆的，玉米处于吐丝期，故将水稻与大豆提取后剩余区域内7月15日的EVI值＞0.5112的区域判定为玉米，余下区域被划分为林地。分辨率为30 m的分类结果如图5所示。

图5 基于决策树的多源遥感数据种植结构提取图

基于多源遥感数据的分类总体精度为87.68%，Kappa系数为0.8162，相较于传统监督分类方法，其总体精度与Kappa系数均有明显提高，各地类具体的用户精度与制图精度如表5所示。

表5 基于决策树的多源遥感数据种植结构分类结果的混淆矩阵

2.3.3 基于决策树的多源转化遥感数据种植结构提取经模型转换后的各类传感器EVI数据可根据MODIS-EVI时序曲线特征阈值，因此对植被细化分类的阈值进行更改。经反复实验确定将8月7日MODISLandsat-8EVI值＞0.5788的区域判定为水稻；将不满足上述条件且9月20日MODISSentinel-2EVI值在0.3006~0.3435范围内的区域判定为大豆；将水稻与大豆提取后，剩余区域内7月15日MODISGF-1EVI值＞0.634的区域判定为玉米，剩余区域被划分为林地。各时相EVI样本平均值变化见表6。

表6 模型转换后各时相EVI样本平均值

对比可知，利用时序曲线提取物候特征对决策树准确阈值能够有效区分玉米、水稻和大豆，并纠正研究区南部水稻与玉米混淆错分的问题，分类结果分辨率为30 m，具体如图6所示。

图6 基于决策树的多源转化遥感数据种植结构提取图

基于多源转化遥感数据的分类总体精度为92.27%，Kappa系数为0.8825，相较于数据未经模型转换的决策树分类方法，玉米制图精度、用户精度、总体精度分别提高了6.25、0.80、4.59个百分点，Kappa系数提高了0.0663，各地类具体的用户精度与制图精度如表7所示。

表7 基于决策树的多源转化遥感数据种植结构分类结果的混淆矩阵

2.3.4 基于决策树的混合像元分解模型将MODISGF-1EVI、MODISLandsat-8EVI、MODISSentinel-2EVI与2.3.2节的分类结果叠加，在玉米分类结果范围内选取玉米的端元，利用FCLS进行线性混合像元分解以实现玉米反演。决策树与混合像元分解结合提取玉米种植结构的总体精度为95.98%，玉米地类提取结果的边界清晰，能够有效地解决决策树分类结果出现的地块边缘错分、细碎图斑的问题(图7)。