基于GF-1与MODIS时空融合的南方丘陵区水稻提取研究*

王建勋，华 丽※，邓世超，王惠东， 2，陈家慧

(1.华中农业大学水土保持研究中心，湖北武汉 430070； 2.武汉大学测绘遥感信息工程国家重点实验室，湖北武汉 430079)

0 引言

水稻作为我国重要的粮食作物，其分布十分广泛。研究水稻空间种植分布可以为水稻精准估产与农业政策制定提供技术支持，对保障国家农业发展和粮食安全具有重要意义[1-3]。传统方法对水稻面积的研究多基于统计数据，该方法难以获得精确的空间分布信息。随着遥感技术的飞速发展，基于遥感技术的水稻提取研究取得了一定成果[4-8]。但这些研究区域多集中在地势平坦、地物类型单一的平原区，针对南方丘陵区地形复杂、地类斑块破碎的特点，如何提高水稻分布面积提取效率和精度成为当前研究的难点和热点。

当前水稻识别与分类主要包括两大类方法：基于单一高分辨率影像和基于多源影像数据的水稻分类提取。高分辨率影像主要包括IKONOS和SPOT等，其优势在于此类影像地物纹理、色调与形状等特征较为丰富，不同地类之间差异大且容易识别。国内外一些学者利用高分辨率数据对小区域的水稻种植区进行面积提取，达到了一定的提取精度[9， 10]。但是影像的时间分辨率低，无法充分利用作物识别所必需的物候信息。而利用多源数据进行水稻的提取，主要采用中高分辨率的Landsat影像与低分辨率的MODIS影像通过数据融合构建时空序列数据集进行水稻提取[11， 12]。该类方法优势在于影像的开源性且易于获取，但仅适合大尺度、地类单一区域的水稻提取，难以适用于地形复杂、地类斑块破碎的南方丘陵区。

自Gao F等[13]提出自适应遥感影像时空数据融合模型(STARFM)以来，许多学者利用该模型完成TM和MODIS数据融合，并进行地物分类与提取研究[14， 15]。但该模型存在无法解决混合像元的缺陷，Zhu等[16]针对这一缺陷提出了ESTARFM融合模型。该模型假定一段时间内地物反射率为线性变化，混合像元由不同地表覆盖线性组合而成，相较于STARFM方法利用单时相高空间分辨率影像与高时间、中低空间分辨率影像融合，ESTARFM利用两景高空间分辨率影像参与融合，可以减少融合过程中产生混合像元的问题，从而提高融合精度。该方法更能好地构建中高分辨率时空数据集并提供满足用户使用的多幅影像，从而为研究南方丘陵区水稻分布提供解决方法。

文章以湖南省衡阳县为研究区域，采用ESTARFM方法将MODIS分别与GF-1WFV和TM数据进行融合，构建高时空分辨率的NDVI时空数据集，再利用TIMESAT软件[17]对获取的NDVI数据集进行滤波处理并提取地物的物候参数，最后通过决策树分类方法提取水稻种植面积并对两类数据融合后的提取结果进行精度评价。

1 材料与方法

1.1 研究区概况

衡阳县地处湖南省衡阳市西北部，湘江中游，西起东经111°59′，东至东经112°46′，南起北纬26°51′，北至北纬27°26′，多为低山丘陵地貌，总面积为25.59万hm2(图1)。该县属亚热带季风气候，年均气温17.9℃左右，温暖湿润的气候成为水稻种植的天然优势。作为全国粮食生产先进县，衡阳县2016年水稻播种面积9.57万hm2，产量57万t，是湖南省主要水稻种植区。旱地作物与水稻交错分布现象普遍，水稻种植模式是单、双季稻皆有，以双季稻为主。单季稻每年5月开始播种育秧， 9月收割。双季稻中早稻4月播种， 7月中旬收割，接着播种晚稻， 10月中旬收割。

图1 研究区

1.2 数据及预处理

1.2.1 Landsat 8数据与GF-1WFV数据

该研究选取云量少于10%以下的Landsat 8和GF-1WFV数据，Landsat 8影像获取时间是2016-03-01，2016-11-28，2016-12-30，该数据由美国地质勘探局提供(https：//earthexplorer.usgs.gov/)。GF-1WFV影像获取时间分别是2016-02-29， 2016-03-28和2016-05-12，数据来源自遥感集市(http：//www.rscloudmart.com/)。Landsat 8 陆地成像仪(OLI)数据包含蓝、绿、红、近红外和2个短波红外波段，GF-1WFV数据包含蓝、绿、红、近红外4个波段，波段信息如表1所示。对GF-1WFV与Landsat影像的红波段和近红外波段进行重投影、大气校正、几何校正等预处理操作。

1.2.2 MODIS数据

该研究选取8d合成的250m分辨率的地表反射率产品数据(MOD09Q1)，该数据是由NASA提供(https：//earthexplorer.usgs.gov/)的标准2级产品，包括红波段(Red)和近红外波段(NIR)的地表反射率，波段信息如表1所示。

数据预处理过程包括：首先使用MODIS Reprojection Tool工具将MODIS数据重投影为UTM-WGS84坐标系； 其次，将数据重采样至30m和16m分辨率，分别与Landsat 8和GF-1WFV数据进行融合处理。用于实验处理的数据特征如表2所示。

表1 Landsat8、GF-1WFV和MODIS影像波段信息

Landsat 8(OLI)GF1-WFVMODIS(Terra)波段号波长宽度(nm)波段号波长宽度(nm)波段号波长宽度(nm)4(Red)630～6803(Red)630～6903(Red)620～6705(NIR)845～8504(NIR)770～8904(NIR)841～876

表2 Landsat8、GF-1WFV和MODIS影像波段信息

Landsat 8(OLI)GF1-WFVMODIS(Terra)获取时间数据用途获取时间数据用途获取时间数据用途2016-03-01融合验证2016-02-29融合输入2016-01-012016-11-28融合输入2016-03-28融合验证…融合输入2016-12-30融合输入2016-05-12融合输入2016-12-27

1.3 研究方法

1.3.1 NDVI计算

NDVI指数是由Rouse[18]等提出，其本质是利用植物叶片在红色与近红外波段的光谱反差组合来设计的，是目前应用最广泛的植被指数。该指数已被证明可以真实地反映地表植被的生长状况。NDVI的计算公式为：

NDVI=(Nir-Red)/(Nir-Red)

(1)

式(1)中，Nir和Red各代表影像的红色与近红外波段。

1.3.2 ESTARFM模型

在忽略大气校正和误差配准影响的前提下，同一分辨率下的MODIS和Landsat数据地表反射率之间存在线性关系(式2)。

L(x，y，tk，B)=a×M(x，y，tk，B)+b

(2)

式(2)中，L、M分别代表Landsat和MODIS数据； (x，y)表示像元位置，tk为获取影像的时间，B表示影像的波段，a、b为线性方程的系数且会随着像元位置的不同发生改变。

当存在2个时期的2对影像，且地物类型在2个时期都没有明显变化，则根据式(2)可得到式(3)[19]。

L(x，y，tp，B)=L(x，y，t0，B)+a×(M(x，y，tp，B)-M(x，y，t0，B))

(3)

使用ESTARFM方法能有效获取30m的融合结果，但该分辨率应用在南方丘陵山区的水稻分布提取显然是不够的，因此该研究尝试使用GF-1数据的16m分辨率替换Landsat数据进行融合并构建时空数据集。由于GF-1数据与MODIS数据在波谱响应和Landsat数据存在差异，因此试验中还加入Landsat融合结果作为对比分析。

ESTARFM模型是一种基于移动窗口的数据融合模型，在进行融合时，考虑从邻近的光谱相似的像元获取辅助信息来提高精度。在一个给定的移动的窗口内，中心像元的预测反射率可以表述为：

(4)

式(4)中，N为包括预测像元在内的相似像元的数目；ω为滑动窗口大小； (xi，yi)是第i个相似像元的位置；Wi是第i个相似像元的权重大小，可描述为对预测中心像元的贡献程度，是依据空间距离、时间距离、光谱距离3项来确定的；Vi是转换系数。

为提高结果的准确率，可使用两期数据计算(xω/2，yω/2)点的反射率，公式如(5)(6)。

L(xω/2，yω/2，tp，B)=μm×L(xω/2，yω/2，tm，B)+μn×L(xω/2，yω/2，tn，B)

(5)

(6)

式(5)、(6)中，L(xω/2，yω/2，tm，B)和L(xω/2，yω/2，tn，B)分别为tm和tn时刻的观测数据融合得到的tp时刻的反射率和二者的权重。

1.4 物候参数提取

植物物候是指植物受生物因子与非生物因子如气候、水文、土壤等影响而出现的以年为周期的自然现象，它包括各种植物的发芽、展叶、开花、叶变色、落叶等现象[20]。利用时间序列NDVI可以提取出与植被作物长势有关的物候参数。该研究使用TIMESAT软件提取6个物候参数以获取季节性植被生长规律信息：作物生长期开始(SOS)、生长期结束(EOS)、生长期长度(LOS)、生长期NDVI最大值(MON)、NDVI基准值(BON)和NDVI幅值(AON)。其中，SOS、EOS及LOS以天为单位计算，涉及NDVI的MON、BON及AON等参数指标的范围为-1～1。

1.5 决策树分类

决策树分类是一种常见的分类方法，利用二叉树的逻辑性质对研究对象进行逐步剔除。遥感应用中经常利用该方法对地物进行提取分类[21-24]。该方法最大的特点和优势是能够综合运用遥感特征量和非遥感特征量等多源数据进行研究物的提取分类[25]。主要思想是：按照一定规则把遥感数据集逐层次由上向下细分从而得到具有不同属性的各个子集。决策树由3部分组成：根节点、内部节点(简称为分支)以及最终节点(简称为叶)，在每一个内部节点处，根据定义的规则，可以将该节点处的数据集进行划分得到满足条件的两个子集，如此针对所有需要分类的对象进行相同处理，直至所有数据被分为预期设定的不同类别，即最终节点。这些规则一般基于专家知识和经验积累进行设置，也可以定义一些算法自动提取[26]。

2 结果与分析

2.1 融合结果

该文采用ESTARFM模型分别对Landsat 8和GF-1WFV数据进行融合构建时空数据集。分别将Landsat 8的2016-11-28和2016-12-30以及GF-1WFV的2016-02-29和2016-05-12四景影像作为输入数据，与对应分辨率的MODIS数据融合得到1年期的多景融合结果。选取Landsat 8的2016-03-01和GF-1WFV的2016-03-28两景影像的子区域与对应时刻融合结果的子区域进行对比分析(图2、图3)。

图2 融合影像与真实Landsat8影像比较及对应波段的相关性系数

图3 融合影像与真实GF-1WFV影像比较及对应波段的相关性系数

对比两种影像融合结果，使用ESTARFM方法融合得到的Landsat和GF-1WFV数据与真实的红光、近红外波段之间仍存有少许差异，且Landsat数据融合结果与真实值之间的相关性系数略高于GF-1数据结果，主要原因是MODIS传感器与Landsat传感器的波谱信息较为相似。但通过对两种影像融合后的NDVI进行相关性分析，发现两类NDVI之间的差距较小(图4)。另外针对南方丘陵山区地形，使用16m的GF-1数据在细节表现上优于30m的Landsat数据，这对后续水稻种植区域提取有很大优势。因此选用质量好的GF-1数据替代Landsat数据进行融合构建时序数据集是可行的。

图4 融合影像NDVI与真实影像NDVI比较及相关性系数

图5 不同地物的标准NDVI时序曲线图

2.2 NDVI拟合及分类结果

分别将预测的44景GF-1和Landsat 8 NDVI数据与原始的NDVI数据构成1年期的时间序列数据导入TIMESAT软件中，选取S-G滤波算法进行NDVI时序平滑。以GF-1WFV数据为例，该研究将研究区地类分为单季稻、双季稻、林地、居民地和水体5类，针对不同地类各随机选取50个样点，以各类样点的平均时序曲线作为各地类标准时序曲线，统计了林地、水体、居民地和单季稻的物候信息，以提取的物候信息作为分类参数，并通过决策树分类法区分5种地类。各地物的NDVI拟合曲线如图5所示。

根据不同地类的物候参数，利用决策树方法建立分类规则并进行水稻识别与提取。以GF-1WFV数据为例，选取质量较好的第17景融合后的影像作为输入影像，从NDVI统计图可以看出，此时的水体生长期NDVI最大值(MON)处在0值附近，且大多都为负值，故用NDVI<0.15和NIR<1100建立第一个规则提取水体； 再利用统计图及物候参数MON、AON建立第二个规则NDVI<0.25和MON<0.25和AON<0.2提取居民地； 利用统计图以及物候参数LOS、MON建立第三个规则： 0.4

表3 基于 GF-1WFV数据的不同地物物候信息参数

SOS(d)EOS(d)LOS(d)BONMONAON林地1153162010.4 0.8 0.3 水体134290155-0.2 0.2 0.3 居民地103268 165 0.1 0.2 0.1 单季稻148 2691210.3 0.6 0.3

图6 分类结果

表4 两种数据分类结果面积统计

类别GF-1WFVLandsat 8象元个数面积(万km2)象元个数面积(万km2)双季稻331.734 4849.2493.212 2838.91水体15.636 740.034.542 240.88居民地19.718 850.485.465 649.19单季稻38.707 099.099.615 686.54林地593.812 51 520.16171.497 81 543.48

对两种分类结果进行像元个数统计得到各地类的面积，结果如表4所示。对比两种数据源融合提取结果得出：居民地和水体面积统计结果比较接近，而在水稻和林地的统计结果上存在一些差异，GF-1WFV数据和MODIS融合数据提取的水稻面积大于Landsat 8融合结果，而在林地面积上，GF-1 WFV则较小。

2.3 分类精度评价

2.3.1 混淆矩阵精度验证

结合实地采样点与Google Earth高分辨率影像样本点，在ENVI软件中结合野外调查建立的解译标志，随机挑选出感兴趣区域ROI作为标准。通过混淆矩阵对感兴趣区与分类结果进行精度验证，由于GF-1WFV数据与Landsat 8数据分辨率不同，结合两者分辨率比例系数关系，该研究共挑选GF-1WFV验证图斑22.418 5万个，Landsat 8验证图斑共5.303 8万个，比例大致为4： 1关系，与两者分辨率比例保持一致。其中GF-1WFV数据总体分类精度为86.37%，kappa系数为0.80，其中双季稻的分类精度为86.70%，单季稻分类精度为62.26%，林地分类精度为93.29%，水体分类精度为82.98%，居民地精度为70.85%，混淆矩阵见表5。Landsat 8数据总体分类精度为80.96%，Kappa系数为0.72，各地类除林地较好外，其他分类结果与GF-1WFV分类均存在差距。根据衡阳市统计局衡阳县2016年国民经济和社会发展统计公报的数据[27]： 2016年衡阳县稻谷播种面积9.566万hm2，该结果与GF-1WFV数据统计结果十分接近，因此可以表明基于GF-1WFV数据提取的水稻分布面积为9.483万hm2，较Landsat 8的9.254万hm2更为准确。

表5 GF-1WFV数据混淆矩阵

类别水体居民地双季稻单季稻林地总像元数用户精度(%)水体22 232475987922 89397.11居民地1 03120 2881 828203623 20387.44双季稻2 5247 90846 1122 4915 60464 63971.34单季稻2011622 8285 6381 50110 33054.58林地8052301 81890099 367103 12096.36总像元数26 79328 63553 1849 056106 517224 185制图精度(%)82.9870.8586.7062.2693.29 总体精度: 86.37%;Kappa系数: 0.80

表6 Landsat 8数据混淆矩阵

类别水体居民地双季稻单季稻林地总像元数用户精度(%)水体5 219421222945 87588.83居民地2223 986616104 28593.02双季稻7931 5949 4794512 01114 32866.16单季稻15445881118561 82161.39林地53883 04838223 05826 62986.59总像元数6 3026 13313 3981 96625 13953 038制图精度(%)81.5264.9970.7456.8791.72 总体精度: 80.96%;Kappa系数: 0.72

衡阳县三面环山，山地区域大面积林地可以被有效提取，研究区是双季稻的主要种植区，耕地基本以种植双季稻为主，因此双季稻提取精度较高。而单季稻提取精度较低则因为单季稻种植面积较少且零星分布于双季稻周围，从而与双季稻混杂在影像上构成混合像元，易被误分为双季稻。研究区湖泊与河流较多，水体识别比较容易，但受空间分辨率影响，少许面积较小的坑塘难以进行有效提取。居民地分为主城区和农村居民地，其中主城区能进行有效区分，而农村居民地会因为混合象元效应造成分类效果不理想。

2.3.2 地形复杂度对水稻提取影响

研究区地形复杂多样，地块破碎，为水稻提取增加了难度。坡度是评价地形复杂度的一个有效而且简单的指标。该研究以5°为区分地形复杂度高低的阈值，将研究区分为两类：一类为地形复杂度较低的区域，该区域地势较为平坦，地形起伏小，田块大都呈现团块聚集且较为规整的特点； 另一类为地形复杂度较高的区域，该类区域主要位于海拔较高或地形起伏较大的低山丘陵区，水田呈现零星分布特点。针对此类特点，对两种不同地形复杂度区域中的水稻提取精度进行了对比分析。在两种区域内分别选取300个单季稻和600个双季稻的验证像元点进行了精度验证，统计结果见表7。结果显示：GF-1 WFV融合后的单季稻在地形复杂度低与地形复杂度高的区域提取精度分别为57.6%与56.3%，略低于Landsat 8融合后的数据的提取精度(63.3%与56.3%)，主要原因是单季稻在衡阳县分布较少，样本数量对精度计算有一定的影响； GF-1 WFV融合后的双季稻在地形复杂度低与地形复杂度高的区域提取精度分别为73.7%与68.7%，明显优于Landsat 8融合后的提取精度(68.5%与63.2%)，不同地形下单双季稻的提取精度与混淆矩阵中的提取精度具有一致性。地形复杂度较低的区域水田主要围绕河流与居民点进行分布，单双季稻混种比较常见，容易造成分类混淆； 而地形复杂度较高区域的内水田主要位于河流或交通道路一侧，周围遍布灌木、果园和小片林地等其他植被，受混合像元的影响大，因此单双季稻的提取精度低于平缓区。根据对比研究结果，说明GF-1 WFV数据应用于地形较为复杂的南方丘陵区的作物提取具有一定的优势。

表7 平缓区与丘陵区单双季稻提取精度对比统计

GF-1 WFVLandsat 8验证像元正确像元错误像元正确率(%)验证像元正确像元错误像元正确率(%)地形复杂度低单季稻30017312757.630019011063.3双季稻60044215873.760041118968.5地形复杂度高单季稻30016913156.330016913156.3双季稻60041218868.760037922163.2

3 结论

该文基于ESTARFM模型用GF-1WFV、Landsat 8两种数据分别与MODIS数据进行时间序列融合重构，计算2016年衡阳县NDVI时序数据，并获取物候参数，通过决策树分类方法进行水稻种植面积提取。得到以下结论。

(1)在地块破碎、地形复杂多样的南方丘陵区，GF-1 WFV数据的空间分辨率高于Landsat 8，故对水稻提取有一定的优势。采用ESTARFM模型分别做了Landsat 8 与MODIS数据、GF-1WFV与MODIS数据的两组融合处理，虽然Landsat 8数据融合结果与原始影像的红、近红外波段真实值的相关性优于GF-1WFV结果，但经过NDVI计算后的结果相关性差异较小。因此可以使用GF-1WFV与MODIS数据构建NDVI时序数据。

(2)利用ESTARFM模型分别构建16m和30m分辨率NDVI时空数据集，对生成的NDVI数据集进行拟合，消除噪声干扰的同时提取各地物物候参数。利用决策树对物候信息进行分类得到两种融合数据的分类精度。其中16m数据集的总体分类精度为86.37%，单季稻和双季稻的提取精度分别为62.26%和86.70%，Kappa系数为0.80。而30m数据集的总体分类精度为80.96%，单季稻和双季稻的提取精度分别为56.87%和70.74%，Kappa系数为0.72。对比两者分类精度，可以看出16m数据集提取的总体精度和水稻精度均优于30m数据集； 此外，对不同地形复杂度下的水稻提取精度进行对比，发现GF-1 WFV数据在地形复杂的区域提取精度优于Landsat 8数据。因此，针对南方丘陵盆地区域复杂地貌和破碎斑块特征下的水稻提取，使用GF-1WFV数据进行融合处理提取的结果是可以信赖的，该方法的适用性可以推广应用于同类典型区域进行作物提取。