基于GF-1与Sentinel-2融合数据的地膜识别方法研究

罗 琪，刘晓龙,2*，史正涛，屈 冉，赵文智

(1.云南师范大学地理学部，云南 昆明 650500；2.内蒙古工业大学信息工程学院，内蒙古 呼和浩特 010051；3.生态环境部卫星环境应用中心，北京 100094；4.北京师范大学地理科学学部，北京 100875)

0 引言

地膜能够调节农田土壤光、热、水、气等条件，缓解作物遭受低温冻害、高温干旱、病虫害等的影响，进而促进作物生长发育、提高作物产量[1,2]，但地膜易破碎且难以回收和降解，大量残膜造成土壤肥力下降、作物减产及土壤污染等农业生产安全与生态环境问题[3-6]。例如，地膜残膜促使土壤中微塑性物质含量增加[7]，其吸收有机污染物/吸附病原体并通过空气传播被生物体吸收，进而进入食物网对人类健康构成威胁[8,9]。随着我国地膜使用面积的增加，准确获取地膜的空间分布信息并进行科学规划与管理对提高农业生产与土壤污染防治及保障人类安全具有重要意义[10]。

目前，国内外学者多采用遥感数据进行大范围的地膜识别研究。例如：Lu等[11]基于250 m低空间分辨率的MODIS-NDVI时间序列数据，利用光谱特征阈值法实现了对光谱特征比较稳定(田块尺度较大，作物均一性高)的新疆覆膜棉田的高精度识别，但目前我国大部分地区种植结构零碎，基于低分辨率数据的地膜识别方法存在一定的局限性。中分辨率数据则可以实现种植结构相对均一地区的地膜识别，如：Lanorte等[12]基于Landsat-8数据，利用支持向量机(SVM)分类法提取了意大利南部阿普利亚地区的地膜；沙先丽[13]基于Landsat-5 TM数据的光谱特征、指数特征和纹理特征，通过构建决策树对新疆覆膜棉区的识别精度达90%以上。但在种植结构破碎的局地地膜精细识别中，中分辨率遥感影像易受混合像元影响，难以准确捕捉地膜的细节信息，因此，多数学者利用高空间分辨率遥感影像提取纹理、形状等空间结构特征进行地膜识别研究。例如：高梦婕等[14]基于GF-2遥感影像探讨最佳分割尺度及最佳特征选取，并采用随机森林、CART决策树、支持向量机分类法提取塑料大棚，结果显示随机森林法的分类精度最高(89.65%)；赵璐等[15]基于GF-2遥感影像，利用面向对象最近邻方法提取农业大棚信息，精度达95.65%；吴锦玉等[16]基于GF-2遥感影像，通过对比分析不同纹理提取算法，发现利用影像纹理可显著提高大棚识别精度。

受限于目前的卫星成像技术，高分辨率遥感数据大都仅包含可见光波段和近红外波段，缺少对地物反射敏感的短波红外波段，而具有短波红外光谱信息的中分辨率数据存在分辨率不足。为此，一些学者尝试基于多源数据和多特征的组合进行地膜识别研究。如哈斯图亚[17]利用GF-1、Landsat-8和Radarsat-2数据的地膜识别精度显著高于单一GF-1数据的地膜识别精度，但Radarsat-2遥感影像重访周期长且成本高，限制了该方法的实用性[18]。考虑到GF-1系列卫星数增加及Sentinel-2数据的高重访(小于5 d)特性，本文尝试将两种数据融合，旨在充分利用GF-1数据的空间细节信息和Sentinel-2数据的光谱信息，探讨多源光学遥感数据融合在地膜识别中的可行性。

1 研究区与数据

1.1 研究区概况

本文研究区为河北省邯郸市邱县(图1)，区域高程为32.8～40.8 m，属北暖温带亚湿润大陆性季风气候，年均气温13 ℃，年均降雨量524.6 mm，气候温和，雨热同季，无霜期长，以棉花、蔬菜地膜种植为主，其中棉花种植面积达33万hm2[19]。

图1 研究区位置(GF-1影像RGB显示)Fig.1 Location of the study area (RGB display of GF-1 image)

1.2 数据及其处理

研究区地膜覆盖时间集中于4月中旬，揭(挑)地膜时间为6月中旬。本文选用的GF-1、Sentinel-2遥感影像获取时间均为2019年5月23日稳定覆膜期，是地膜光学遥感识别的理想时间窗口。

(1)GF-1数据。时间分辨率为5 d，PMS传感器可获得空间分辨率为2 m的全色影像和8 m的多光谱(BlueGreenRedNIR波段)影像。首先对影像进行辐射定标，并利用FLAASH模块进行大气校正，然后利用GF-1卫星数据自带的RPC文件进行正射校正。通过Gram-Schmidt Spectral Sharpening算法对影像的多维特征进行正交化处理，可消除数据冗余且光谱保真度高，能明显提高影像的空间细节信息，融合效果好，尤其是在黑色背景区域较多的情况下[20-22]。为保证GF-1数据的光谱和空间分辨率，本文采用该算法融合GF-1多光谱和全色波段，最终获得空间分辨率为2 m的多光谱影像，用于对比分析本文方法的地膜识别精度。

(2)Sentinel-2数据。Sentinel-2当前在轨运行A、B两颗星，覆盖13个光谱波段，地面分辨率分别为10 m(BlueGreenRedNIR波段)、20 m(Vegetation Red EdgeSWIR波段)、60 m(Coastal aerosolWater vapourSWIR-Cirrus波段)，单颗卫星的重访周期为10 d，A、B双星互补后在赤道地区的重访周期达5 d。本文对下载的Sentinel-2-L1C产品用Sen2Cor插件进行大气校正，并根据研究需求，将数据的BlueGreenRedNIRSWIR波段合成并重采样为空间分辨率为10 m的多光谱影像。

2 研究方法

本文将GF-1全色影像与Sentinel-2多光谱影像进行NNDiffuse Pan Sharpening融合后，利用分层分类(HC)法进行地膜识别，并与基于GF-1遥感影像/多源融合影像的深度学习(DL)法、传统机器学习分类方法进行精度对比，验证方法的可行性。

2.1 融合算法

Gram-Schmidt Spectral Sharpening算法用于GF-1影像融合效果佳[23,24]，但该算法耗时长，且对不同数据源的影像配准精度要求高，更适用于同一平台传感器的影像融合[25],故本文仅用该算法对GF-1数据进行融合。NNDiffuse Pan Sharpening融合方法将单个像元光谱值作为融合运算的最小元素，通过混合模型产生分辨率增强的融合影像，该方法运算速度快，能较好地保留原影像的色彩和纹理，且融合后影像的空间特征和光谱特征不会发生较大变化[26]，但需注意的是，待融合遥感影像多光谱波段间的波长范围不能交叉覆盖，且全色波段能够基本覆盖多光谱所有波段的波长范围[27]，故本文利用该算法对GF-1、Sentinel-2数据进行融合。首先，以预处理后的GF-1全色影像为参考影像，Sentinel-2多光谱影像为待校正影像，采用6个地面控制点和一阶多项式变换模型进行几何精校正，减少非系统因素产生的误差，校正后影像的几何配准偏差小于0.5个像元；再以GF-1全色影像为高分辨率影像，Sentinel-2多光谱影像为低分辨率影像，进行NNDiffuse Pan Sharpening融合，得到空间分辨率为2 m的多光谱(BlueGreenRedNIRSWIR波段)影像。

2.2 地膜识别方法

对地膜识别的多特征信息进行可分离性分析有利于区分地膜和其他地物。本文通过目视解译获取190个样本，计算研究区不同土地类型(地膜、植被、水体、裸土和休耕地、建筑物)在融合数据RedGreenBlueNIRSWIR波段的反射率均值，并基于NDVI[28]、改进的归一化差异水体指数(MNDWI)[29]、归一化建筑指数(NDBI)[30]、建筑用地指数(IBI)[31]计算不同土地覆盖类型的反射率(图2)，进行地膜与其他地物的可分离性分析。由图2可知，地膜在Blue、Green、Red、NIR、SWIR波段的反射率均高于裸土和休耕地，且地膜在Blue、Green、Red、SWIR波段的反射率明显高于水体、植被，但其与建筑物的光谱反射率较为接近。由于利用归一化差异水体指数(NDWI)[32]提取水体时存在许多非水体信息，且建筑物等阴影在Green和NIR波段的波谱特征与水体相似[29]，故利用SWIR波段数据，采用MNDWI提取水体信息，再利用NDVI提取植被，实现地膜与水体、植被的分离；同时，利用NDBI和IBI提取建筑物，但这两种指数难以实现地膜与建筑物的分离。像元形状指数(PSI)[33]是利用中心像元及其相邻像元的关系，计算具有光谱相似性相邻像元组的维数，并拉大不同方向邻域像元间的特征值，以增强影像的空间细节信息。利用PSI能提取高分辨率影像的形状和结构特征，可弥补光谱特征的不足，在规则纹理人工地类提取中具有较强优势[16]，故本文利用PSI提取研究区建筑物信息。PSI的计算步骤为：1)定义围绕中心像元的方向线，不同方向线的长度通过相邻像元间的光谱同质性测度值(式(1))和阈值确定，每条方向线从中心像元向两边扩展，扩展条件包括光谱约束条件(PHi(x,y)

图2 不同土地覆盖类型的反射率与遥感指数值(均值)Fig.2 Reflectance and remote sensing index values of different land cover types (mean value)

(1)

(2)

(3)

分层分类(HC)法是根据各类地物的光谱特征，按照一定的准则对每一层级进行分解，以达到高精度分类或专门提取某一地类的目的。针对待分类地物的光谱特征选择不同的特征参数或分类方法，逐层提取地物信息并进行信息掩膜处理，从而消除对其他地物提取产生的影响，然后将逐层分类结果叠加得到最终分类结果[34]。首先，基于融合影像利用MNDWI、NDVI提取水体和植被，并对提取结果进行掩膜处理，然后利用PSI提取建筑物，并再次进行掩膜处理。本文是在覆膜稳定时期进行的地膜识别研究，适逢农作物生长发育期，裸土和休耕地分布较少。通过分析地膜与裸土和休耕地在Blue、Green、Red、NIR、SWIR波段的光谱反射率，可知地膜的反射率高于裸土和休耕地的反射率(图2)，利用反射率差异，基于最大类间方差法(OSTU)[35]实现地膜与裸土和休耕地的区分(图3)。OSTU计算简单，自适应性强，对类间方差为单峰的影像分割效果好[36-39]，故本文通过该方法确定水体、植被、建筑物及地膜的提取阈值。

图3 分层分类框架Fig.3 Framework of hierarchical classification

2.3 对比验证方法

为验证本文方法在地膜识别中的可行性，将其与基于深度学习、传统机器学习的分类方法进行比较。其中，深度学习法通过分层处理机制，构建含有多个隐藏层的学习模型，进行影像特征的自动提取和选择，具有很强的特征表达、容错能力和鲁棒性[40]，但其需要大量标记的训练样本，在小尺度遥感影像数据集情况下会降低分类性能，且需要对分类结果进行优化处理，以提高地物边界的细化程度，在遥感影像分类应用领域存在局限[41]；另外，该方法对计算机性能要求较高，需要使用GPU或TPU加快学习模型的训练速度，成本较高[42,43]。目前，基于深度学习的遥感影像分类以监督分类为主，标记样本数据(影像块、语义分割)的数量和质量也将影响模型应用的精度[44]。支持向量机能较好地解决小样本统计学习问题，对高维数据不敏感，具有较高适应性和泛化能力，且能避免过学习问题[45]。随机森林法的计算效率高，适合高维数据处理，对数据噪声和小样本数据不敏感，具有较好的鲁棒性和稳定性，并且在运算过程中人工干扰少，不易出现过拟合问题[46]。

3 结果与分析

3.1 融合结果及精度评价

本文采用定性和定量相结合的方法对GF-1与Sentinel-2的融合影像进行评价，以验证本文融合方法是否有益于地膜识别。1)为便于定性分析融合影像的颜色、清晰度和空间纹理特征等细节信息，在研究区裁剪出地膜分布集中且与其他地物边界明显的局部区域(图4)。由图4可知，相对于Sentinel-2多光谱影像，融合影像的空间分辨率明显提高，黑色圆圈中地物的纹理特征、边缘细节更清晰，但其相比GF-1影像清晰度有所降低。2)利用目前常用的平均梯度、信息熵和相关系数[47,48]指标对融合影像和GF-1、Sentinel-2影像的各波段数据进行评价(表1)。由表1可知，融合影像各波段的信息熵均接近7.0，保留了较多原始影像的信息量，且平均梯度与GF-1接近，但低于Sentinel-2。由于受多源数据获取方式以及卫星姿态、扰动等影响，多源遥感影像的几何精校正存在偏差，使融合方法存在不确定性[49]；另外GF-1数据与Sentinel-2数据的空间分辨率存在差异，且受地表复杂性等因素影响，难以获取理想参考影像，加之GF-1融合后的数据与Sentinel-2数据存在不可避免的偏差，导致利用GF-1、Sentinel-2数据和融合数据进行的指标评价存在较大的不确定性[50]。融合影像各波段相关系数高于0.8，且与Sentinel-2数据的相关系数达0.98以上，说明融合数据的保真度较高。

图4 原始影像与融合影像对比Fig.4 Comparison of original image and fusion image

表1 融合质量评价结果Table 1 Fusion quality evaluation results

3.2 地膜识别结果比较

本文设计不同分类方案(表2)，将分层分类法的识别结果与基于深度学习、支持向量机和随机森林方法的识别结果进行比较。为更好地体现不同分类方法的可比性，将研究区地类分为地膜和其他(建筑物、裸土和休耕地、植被、水体)两个类别。

表2 不同分类方案Table 2 Different classification schemes

本文利用混淆矩阵法，选用总体精度(OA)、Kappa系数、生产者精度(PA)、用户精度(UA)4个指标[51]对地膜识别结果进行定量分析(表3)，精度验证的参考数据为人工目视选取的173个样本。

表3 不同方案识别精度Table 3 Accuracy of different identification schemes

由表3可知，Fus+HC方案的分类精度最高，总体精度达94.87%，Kappa系数达0.89；其次是GF-1+DL和Fus+RF方案，总体精度分别为93.14%、92.49%，Kappa系数分别为0.85、0.83。与GF-1+DL方案相比，Fus+DL方案总体精度降低1.16%，Kappa系数降低0.03。多次实验表明，基于深度学习的识别方法易受影像分辨率的影响，用于混合像元较多的中、高分辨率影像时，其识别效果有限[41]。Fus+RF方案比GF-1+RF方案总体精度提高5.71%，Kappa系数提高0.11；Fus+SVM方案比GF-1+SVM方案总体精度提高1.59%，Kappa系数提高0.03。整体上，除基于深度学习的地膜识别方法外，基于多源融合数据的地膜识别精度均高于基于GF-1遥感数据的识别精度。

从识别结果的空间分布看(图5)：GF-1+SVM和GF-1+RF方案分类结果中地膜与建筑物存在较多混淆，尤其是在地物复杂区域；Fus+RF和Fus+SVM方案将建筑物区域内的植被误分为地膜，导致分类结果出现“椒盐”现象；GF-1+DL和Fus+DL方案对“椒盐”现象进行了改善，地膜识别精度均高于90%，但其识别结果中地物边界过于平滑，需引入影像后处理方法(如形态学滤波[52]、全连接CRFs[53]等)优化分类结果，使分类后的地物边界更接近真实地物边界。Fus+HC方案识别出的地物边界轮廓与真实地物边界基本一致，空间细节特征表现明显，分类效果更好，故本文利用Fus+HC方案提取研究区地膜的空间分布(图6)。

图5 不同方案分类结果Fig.5 Classification results of different schemes

图6 基于分层分类法的地膜识别结果Fig.6 Identification results of plastic mulch based on hierarchical classification

4 结论与讨论

针对单一遥感数据源在地膜识别中存在混合像元、光谱特征不足等问题，本文提出一种基于GF-1与Sentinel-2融合数据的地膜识别方法，并与机器学习、深度学习方法提取结果进行精度对比，结论如下：1)基于大样本空间的深度学习、传统机器学习分类方法虽然能够获得高地物识别精度，但其对于有限样本空间的地物识别不如阈值法，本文基于多源融合数据的分层分类法通过同时相GF-1数据和Sentinel-2数据融合，在地膜识别应用中获得了较为可靠的结果(总体精度为94.87%，Kappa系数为0.89)，明显高于Fus+DL、Fus+RF和Fus+SVM方案地膜识别的总体精度(分别为91.98%、92.49%和87.50%)和Kappa系数(分别为0.82、0.83、0.74)。2)地膜具有较强的方向反射特性，多角度成像是实现地膜识别的关键。本文未考虑不同卫星数据源多角度特征对识别精度的影响，后续将进一步分析识别精度与成像几何问题，并探讨不同空间分辨率(8 m、10 m、20 m)多光谱影像组合对识别精度的影响，以期获得更佳的识别效果。随着Sentinel-2数据和我国GF-1系列数据的发展，这两种数据的结合可实现多源遥感数据的优势互补，为解决地膜识别以及积雪和冰川识别等问题提供了新思路。