GEE环境下融合主被动遥感数据的冬小麦识别技术

2021-10-13马战林刘昌华薛华柱李静茹周俊利

农业机械学报 2021年9期

马战林刘昌华薛华柱李静茹房旭周俊利

(1.河南理工大学测绘与国土信息工程学院，焦作 454003； 2.河南省遥感测绘院，郑州 450003)

0 引言

随着我国粮食生产的市场化，冬小麦的种植面积呈现年际变化大的特点，及时获取冬小麦精确的种植面积信息及空间分布，对促进我国农业生产发展和粮食安全具有重要意义[1]。

遥感观测技术在大范围作物种植面积和结构的监测上已被广泛应用[2]。但由于技术条件限制和工作原理的不同，单一的传感器难以全面反映复杂的地表覆被特征。光学遥感数据受光谱分辨率、空间分辨率、光谱波段的因素及限制，导致同物异谱、异物同谱的现象存在，影响了地物识别的精度。已有学者通过不同方法提高地物分类精度[3-9]。上述方法需要足够的光学数据支持[10]。但在实际应用中，受云雨、光照等因素的影响，光学遥感数据源的质量和数量无法保障，一定程度上限制了地面作物信息的准确提取[11]。

合成孔径雷达(Synthetic aperture radar,SAR)不受云雨天气影响，具有全天时、全天候监测等优点。然而，SAR数据的信号易受相干斑噪声的干扰，影响对目标地物的提取，单一时相的SAR数据很难对地物进行精确提取。有研究表明，使用多时相或多极化SAR数据比单时相或单极化SAR数据能获得更好的分类结果[12]。若能将光学和微波遥感数据结合，利用两者数据的优点进行地物分类，则分类精度会有较大提升。文献[13-16]利用最大似然分类(Maximum likelihood classification,MLC)、人工神经网络(Artificial neural network, ANN)、支持向量机(Support vector machine,SVM)和监督分类等方法，对比光学遥感、微波遥感、光学遥感与微波遥感数据结合3种方式的作物分类精度，得到光学和微波遥感数据相结合的分类精度最高。

相比于单期光学数据和SAR数据的融合，时序数据进行融合应用时，作物分类精度进一步提高。文献[17]结合Landsat-8和Sentinel-1A数据构建时间序列遥感数据以提取作物物候特征，利用卷积神经网络算法进行不同作物的分类，并取得了较高的分类精度。文献[18]以多时相SAR(Sentinel-1A)和Landsat-8光学影像对南京市高淳区的冬小麦种植面积进行提取，结果表明添加纹理特征的时序SAR数据和光学数据结合时分类精度最高。

除数据源外，分类算法的选择也会直接影响分类结果[19]。近年来随着机器学习分类算法的发展，随机森林算法(RF)计算精度高，模型训练时间少，能够确定变量在模型中的相对重要性，同时对训练样本数量和质量的敏感度低，而被广泛应用于作物识别中[20-21]。

以上研究中均需将所需遥感数据下载至本地进行预处理，所需时间较长，而谷歌地球引擎(Google Earth Engine，GEE)是一个提供全球尺度地球观测数据储存和用户友好界面的数据访问平台，具有强大的数据处理能力，仅需内置的少数代码即可完成数据获取、预处理及算法的应用。因此，本文在GEE平台上，充分利用SAR极化数据所包含的地物结构信息和光学数据的光谱信息，结合RF算法，探究时间序列SAR数据、融合时间序列SAR和光学数据的不同特征值组合对冬小麦识别精度的影响。

1 研究区域与数据来源

1.1 研究区概况

河南省襄城县、鄢陵县和临颍县位于河南省中部地区(图1)，是中国北方典型的冬小麦种植区和河南省高标准冬小麦种植示范区。三县地跨113°2′～114°2′E，33°41′～34°15′N，总面积为2 612.60 km2，属北温带季风气候区，热量资源丰富，雨量丰沛，光照充足；年平均气温14.3～14.7℃，年平均降水量为671.1～736.0 mm，无霜期为217.5 d；地形地貌以平原为主，占75.81%，土壤肥沃，具有较高的自然肥力；农作物以冬小麦、夏玉米为主，其中，小麦播种面积高达夏粮总种植面积的90%以上。该区域的冬小麦10月中旬种植，到次年6月上旬收割。鄢陵县种植结构复杂，地块细碎，与襄城县和临颍县的连片种植冬小麦形成鲜明对比，在分类精度的评价上具有较好的区域对比性，可验证融合Sentinel主被动遥感数据在提取冬小麦方面的应用潜力。

1.2 GEE及应用数据

GEE云平台存储了Sentinel数据、MODIS数据集、降水数据、The National Agriculture Imagery Program(NAIP)数据、海洋表面温度数据、Landsat数据、CHIRPS(Climate hazards group infraRed precipitation with station data)气候数据和海拔数据等海量数据，容量超过5PB，且每天增加约5 000幅影像，可以解决大面积土地覆盖制图方面最重要的数据存储下载问题[22]。用户可以使用基于网络的集成开发环境代码编辑器分析所有可用的遥感图像，而无需将这些数据下载到本地机器上。这样，用户可以轻松访问、选择和处理大研究区域的大量数据。除了快速处理之外，它提供了几个包含大量算法的软件包，可以简化专家和非专家对遥感工具的访问。GEE也允许用户上传自己的栅格和矢量数据(例如GeoTIFF或Shape文件)进行分析，完全控制访问[23]。

Sentinel-1主动微波遥感卫星由两颗极轨卫星A星和B星组成，搭载C波段的合成孔径雷达(SAR)传感器，重访周期小于10 d，本文采用分辨率为10 m、极化方式为VV和VH的后向散射系数数据。Sentinel-2由Sentinel-2A和Sentinel-2B两颗高分辨率卫星组成，单颗卫星的重访周期为10 d，两颗互补，重访周期为5 d。其搭载的MSI传感器有可见光、近红外和短波红外波段，空间分辨率在不同的波段有10、20、60 m。本研究主要选取该卫星分辨率10 m的蓝光、绿光、红光和近红外波段。选用的Sentinel-1微波遥感数据和Sentinel-2 MSI多光谱遥感数据均在GEE云平台上进行调用、处理。在GEE云平台上使用函数直接获取不同极化方式的后向散射系数，减轻处理SAR数据的难度。根据图像的大气校正状态，Sentinel-2数据可分为两个不同的级别:1C级和2A级。1C级为大气顶部图像，需要大气校正。2A级数据为大气底部反射数据，已进行大气校正。选用2A级数据。同时考虑云雾对数据质量的影响，在GEE中设置云量低于10%。遥感影像参数如表1所示。

表1 遥感影像参数Tab.1 Remote sensing image parameters

本文使用高分六号(GF-6)数据，该卫星配备2 m全色/8 m型多光谱高分辨率相机以及16 m型多光谱中分辨率宽幅相机。应用2020年4月6日和2020年5月17日两景GF-6融合后的2 m多光谱数据，分析评价分辨率10 m Sentinel主被动数据在碎部地区的分类效果。

1.3 特征变量与数据集

结合研究区植被的生态环境和物候特点，选取Sentinel-2影像的光谱反射率及相关植被指数和熵值纹理(Entropy,ENT)共11个特征变量。Sentinel-1 SAR数据的特征变量包含极化特征变量和纹理特征变量，为避免多纹理特征带来的特征信息冗余[24]，在保留最大信息量、便于计算和较高分类精度的前提下，选取灰度共生矩阵(Gray-level co-occurrence matrix,GLCM)生成的角二矩阵(Angular second moment,ASM)、对比度(Contrast,CONTRAST)、相关性(Correlation,CORR)和熵值4个纹理特征变量，共6个特征变量[25]。研究中共17个特征变量，如表2所示。

表2 特征变量Tab.2 Characteristic variable

归一化差异植被指数(Normalized difference vegetation index,NDVI)研究较高吸收和叶绿素反射的波段来识别光合作用植被。归一化水体指数(Normalized difference water index,NDWI)用来识别水体，NDVI和NDWI的结合可用于提取一些沿河畔、坑塘生长的植物的信息[26]。增强型植被指数(Enhanced vegetation index,EVI)包括近红外、红色和蓝色波段[27]，可在高生物量区域获得更好的灵敏度，并通过解耦冠层背景对不同物种的冠层结构变化做出更好的响应[28]。绿色归一化差异植被指数(Green normalized difference vegetation index,GNDVI)已被开发用于估计叶片叶绿素含量[29]。绿叶指数(Green leaf index,GLI)为颜色植被指数，反映植被的生长至衰老颜色变化信息。优化土壤调节植被指数(Optimization soil-adjust vegetation index,OSAVI)能够反映冬小麦的生长状况信息[30]。

受制于研究区作物种植结构和田块零碎的影响，研究中未使用Sentinel-2数据中分辨率为20 m的红边波段。纹理特征的提取能够较好地描述纹理的细节和随机性，适于地物分布复杂的SAR影像[31]。

为获取襄城县、临颍县、鄢陵县主要地物类型分布，结合Google Earth Engine，在2020年3月15日到2020年5月12日冬小麦生长期间进行实地采样。依据研究区的农田、道路、房屋、水体占土地覆被总面积90%左右，将覆被类型分为建筑用地、冬小麦、水体、林地、其他植被、其他用地6大类别。建筑用地包括房屋道路等建设用地，其他植被包含蔬菜、花卉、瓜果等，其他用地包括裸土、裸岩以及采石场等。应用中海达V30 型GPS仪器进行实地采样，共获得1 051个样本点(表3)，其中田地采样点均在大面积(边长100 m以上)田块中心采集。验证数据集应用GPS在研究区内随机均匀地对不同类别的数据进行区域采集，生成矢量数据，应用矢量范围数据对遥感图像进行裁剪得到验证数据集像元数量。

表3 样本数据集数量Tab.3 Sample data sets

2 研究方法

2.1 随机森林算法

随机森林算法的预测结果是在对组成森林的多棵决策树的决策结果进行众数求解得出的。每次采用Bootstrp采样技术从原始的训练数据集中有放回地抽取大约2/3的样本用于对当前决策树模型进行训练，故原始训练集中的部分样本可能被同时用来对决策树进行训练，也有可能从未被任意一棵决策树使用。抽取剩余的约1/3的样本构成对随机森林模型进行性能评估的袋外数据，计算该模型的预测错误率，该错误率被称为袋外误差。随机森林分类的基本思想：使用Bootstrap抽样从原始训练集中提取k个样本，每个样本的容量与原始训练集相同；为k个样本建立k个决策树模型，以获得k个分类结果；对每个记录进行表决以确定其最终分类。虽然随机森林分类器计算速度比单个决策树慢，需要自身算法推断出超出范围的独立变量或非独立变量，但其能处理高维特征，不易产生过拟合，模型训练速度比较快，可以产生高准确度的分类器，分类效果比较好。但随机森林决策树的数量对随机森林算法的效率有重要影响。当决策树数量偏少时，分类精度降低；决策树偏多时，分类精度会趋于稳定，但运算速度缓慢。有研究表明，在1～100范围内，当决策树大于50时，再继续增大随机森林决策树数量已无意义，因此本文选取决策树数量为50[32]。

2.2 精度评价

为有效评估不同数据对不同覆被类型的提取精度，本文从分类总体精度、混淆矩阵、Kappa系数、制图精度和用户精度5方面进行比较分析[33]。

3 实验结果与分析

3.1 融合月平均时序Sentinel-1 SAR数据分类

在气候湿润、多雾多云的地方，光学数据难以获取。Sentinel-1卫星搭载的C波段波长远高于云粒子的典型尺寸，且自身发出能量脉冲，故可穿透大部分云雾并不受光照条件的影响，被广泛应用于农作物提取研究。单利用某一时相的SAR极化特征数据有极大的局限性。有研究表明，将月尺度上观测数据进行均值合成后，能够降低降水对分类精度的影响，较大程度提高禾谷类作物的分类精度[12]。因此本研究选取冬小麦生长期内2019年11月到2020年5月Sentinel-1 SAR数据，在月尺度上进行均值合成，用于冬小麦识别研究。

融合7个月均值Sentinel-1 SAR数据的VV、VH极化雷达后向散射系数和其通过GLCM计算的纹理特征，应用RF算法进行分类。在GLCM计算纹理特征共生矩阵大小的选择上，选取4、8、16邻域数值进行计算，分类精度最高为4邻域。此时，融合多时相Sentinel-1极化特征和纹理特征的分类总体精度为91.00%，Kappa系数为0.84。冬小麦的制图精度为95.66%，用户精度为96.58%，结果如表4和图2所示。从碎部地区分类结果与GF-6的融合多光谱分辨率2 m数据(图3)对比发现，加入纹理特征导致覆被结构复杂的区域呈现斑块现象，碎部覆被细节体现不明显，这是共生矩阵计算纹理特征时应用周围像素参与计算和破碎地块像素变化不规律导致的结果。

表4 多时相Sentinel-1 SAR数据的极化特征和纹理特征分类结果Tab.4 Classification results of multi-temporal Sentinel-1 SAR data combining texture features and polarization features

在遥感影像分类中，特征变量权重能够增强对特征变量之间的理解。在RF算法中，基于袋外数据，可以对输入数据的权重进行验证，获得每个特征的权重值，用平均精度减少值来表示，其原理是将某特征参数数值变为随机数，计算其对模型准确率的影响，根据得到的精度减少值来计算此参数的权重，该值越大表明该变量的权重越大[24]。GEE中通过explain()函数可输出各特征的权重，特征变量权重的总和与变量的数目相关，并非固定值，其权重在自身变量组中具有相对意义。图4表示在分类过程中各变量的特征值权重，各特征变量在不同时期的重要性差异较大。对比VH、VV与纹理特征的特征值，在每个月平均值中，除去2019年11月和2020年2月，极化数据特征对分类的贡献率相对纹理特征较大。有研究表明，分类精度和特征重要性与特征值之间的差异性呈正比例关系[13]。从图5中发现，2019年11月和2020年2月重要性降低的原因是冬小麦与其他用地、其他植被、林地的极化值差异较小。

去除纹理特征，仅使用月平均的VH和VV极化数据进行冬小麦提取，结果如表5和图6所示。此时，融合月平均VH和VV极化数据的分类总体精度达到85.93%,Kappa系数为0.75，相对于添加纹理特征，总体精度降低5.07个百分点，Kappa系数降低10.71%。由此表明，综合极化数据和纹理特征使分类精度较高，说明多源特征量的综合有利于总体分类结果的改善。受冬小麦地块连片种植的影响，冬小麦的制图精度和用户精度变化仅为0.01、0.4个百分点，识别精度受纹理特征影响较小。对于其他覆被类型，虽单使用融合极化月平均数据识别精度较低，但破碎地块细节体现明显，其原因有待进一步研究。

表5 多时相Sentinel-1 SAR的极化特征分类结果Tab.5 Classification results of multi-temporal Sentinel-1 SAR data polarization features

3.2 融合Sentinel-1 SAR和Sentinel-2光学数据分类

SAR数据的地物后向散射特性异于光学遥感影像。光学数据反映目标体的光谱特性，SAR数据的穿透性不仅能够获取植被的表面信息，对植被的叶、茎、枝干等信息也有一定的反映，获取的是不同于光学数据的地物信息[34]。图7为不同覆被在冬小麦生长期间的Savizky-Golay滤波NDVI变化图。越冬期过后，2—3月温度上升，冬小麦进入返青期和拔节期，叶绿素含量逐渐升高，而林草与其他植被在此时生长较为缓慢，导致NDVI特征值差异较大。3—4月为冬小麦拔节期/孕穗期，整体植被指数较高，NDVI差异值减小。在5月中旬以后，冬小麦在灌浆期、收获期叶片含水率、叶绿素含量降低，叶黄素含量升高，植被指数呈陡然下降趋势，但其他植被的植被指数开始升高，此时冬小麦NDVI会与其他覆被的NDVI存在交叉。基于特征值差异值越大，RF算法分类精度越高的原则，选取返青期作为光学数据的输入，数据日期为2020年2月22日。在GEE中，通过get(‘CLOUDY_PIXEL_PERCENTAGE’)函数，可获取该地区2020年2月22日Sentinel-2 光学影像的含云量，输出结果为0.45%。

首先，对返青期的光学数据进行分类研究，结果如表6和图8所示。返青期Sentinel-2光学数据的分类总体精度为94.01%,Kappa系数为0.89，冬小麦制图精度和用户精度分别为96.34%、97.67%。相对于加纹理特征的融合Sentinel-1 SAR数据分类精度，返青期Sentinel-2的总体精度提高3.01个百分点，Kappa系数提高5.95%。但冬小麦的制图精度和用户精度差距为1个百分点左右。单期光学影像相对于时序SAR数据，在识别冬小麦过程中优势明显。

表6 冬小麦返青期Sentinel-2光学数据分类结果Tab.6 Sentinel-2 optical data classification results of winter wheat at greenback period

图9为冬小麦返青期Sentinel-2光学数据特征值的权重。光学数据熵值纹理特征此时为0，植被指数特征值权重相对均衡，主要是在选取各类指数过程中，不同指数反映了作物生长阶段的一些生长特征。但遥感信息反映的只是地表或作物群体瞬间的物理状况，而作物生长是一个随时间变化的连续过程。在进一步研究中，应考虑光学遥感信息与物候信息的结合，突出作物生长-衰老交替引起的季节性植被变化，降低同物异谱和异物同谱现象，特别在混合像元及破碎地块对冬小麦的提取上具有很强的优势。但本文主要研究时序SAR数据和光学数据的互补性，开发一种应用少量光学数据或无光学数据的分类方法。

融合时序Sentinel-1 SAR、SAR纹理特征和返青期Sentinel-2光学数据，进行冬小麦提取研究。表7和图10结果表明，当时序Sentinel-1 SAR加入纹理特征时，与返青期Sentinel-2光学数据融合的分类总体精度为95.31%，Kappa系数为0.92，冬小麦的制图精度及用户精度为97.71%和97.23%。融合Sentinel主被动遥感数据的总体精度和Kappa系数，相对于加入纹理特征的融合月平均Sentinel-1 SAR数据提高4.31个百分点和9.52%，冬小麦制图精度和用户精度提高2个百分点左右，相对于返青期的Sentinel-2光学数据提高1.3个百分点和3.37%。融合Sentinel主被动遥感数据能够提高分类总体精度和冬小麦的识别精度。

表7 多时相Sentinel-1 SAR数据的极化特征、纹理特征融合单期光学数据分类结果Tab.7 Classification results of multi-temporal Sentinel-1 SAR data polarization features and texture features integrating a single optical data

在融合Sentinel主被动遥感数据特征值权重方面，如图11所示，光学数据的权重高于SAR数据的权重，这也是应用RF算法时，单基于光学数据高于融合月平均时序SAR数据的原因。纹理特征值的权重相对较弱。时序Sentinel-1 SAR不加入纹理特征时，与返青期Sentinel-2光学数据融合的分类总体精度为95.78%，Kappa系数为0.92，相对于未添加纹理特征的融合月平均Sentinel-1 SAR数据提高9.85个百分点和22.67%，冬小麦的制图精度及用户精度为98.56%和97.92%，如表8和图12所示。相对于添加纹理特征的联合数据总体精度增加0.47个百分点，两者精度差别不明显。但在冬小麦及其他地物提取方面，添加纹理特征的融合数据冬小麦识别精度降低0.8个百分点左右，特别是其他用地的制图精度降低18.86个百分点。综合对比表4和表5，表7和表8，纹理特征对冬小麦的识别精度影响程度较小。可见，融合Sentinel-1 SAR和Sentinel-2光学数据分类与融合Sentinel-1 SAR数据同样表现为，添加纹理特征的融合数据在覆被类型复杂地块的分类精度上，并不随特征值的增加而占优势。

表8 多时相Sentinel-1 SAR极化特征融合单期光学数据分类结果Tab.8 Classification results of multi-temporal Sentinel-1 SAR data polarization features integrating a single optical data