基于融合MODIS与Landsat数据的洞庭湖区水稻面积提取.docx
《基于融合MODIS与Landsat数据的洞庭湖区水稻面积提取.docx》由会员分享，可在线阅读，更多相关《基于融合MODIS与Landsat数据的洞庭湖区水稻面积提取.docx（26页珍藏版）》请在冰豆网上搜索。
基于融合MODIS与Landsat数据的洞庭湖区水稻面积提取
基于多时相Landsat数据融合的洞庭湖区
水稻面积提取
张猛1,2,曾永年1,2
（1.中南大学地球科学与信息物理学院，长沙410083；2.中南大学空间信息技术与可持续发展研究中心，长沙410083）
摘要：
洞庭湖区作为中国重要的商品粮基地，水稻种植面积的变化对国家粮食安全有重要的影响，准确获取水稻面积及其变化显的十分重要。
为解决数据缺失问题,该文利用STARFM（SpatialAndTemporalAdaptiveReflectanceFusionModel）模型融合高时间分辨率的MODIS数据与中等空间分辨率的Landsat数据,得到时序LandsatNDVI数据。
经S-G函数平滑处理，参考作物物候特征及可分离性分析（J-M距离）得到最佳时期的LandsatNDVI组合，结合Landsat8OLI影像对水稻种植面积进行提取。
结果显示，该方法能够有效的提取水稻种植面积，总体分类精度94.52%，Kappa系数为0.9128。
水稻分布几乎覆盖整个研究区，水稻种植总面积达7.88×105hm2。
双季稻种植面积为7.75×105hm2，主要集中于湖区北部及西北部，且分布较连续。
一季稻种植面积为1.3×104hm2，分布相对零散，有小范围集中于湖区中部及西北部。
关键词：
洞庭湖区；水稻；MODIS；Landsat；数据融合
中国分类号：
TP79；S127；F301.24文献标识码：
A
0引言
洞庭湖区作为中国的商品粮基地和湖南农业主产区，以种植粮食、棉花为主，分别占到全省的50.3%和89.3%[1-2]。
农作物面积的变化不仅关系着国家的粮食安全，也影响区域环境和气候变化，以及社会经济发展的决策。
随着洞庭湖区城市化速度的加快，耕地非农化现象突出，同时人口增长对粮食的需求量不断增大。
因此，洞庭湖区农作物种植面积，尤其是水稻种植面积及其变化信息尤为重要。
卫星遥感技术已被广泛应用于农作物分析[3-6]，MODIS数据由于其较高的时间分辨率，基于时间序列的MODIS数据的作物种植面积提取已开展了较多的研究。
Vintrou等利用时间序列的MODIS13Q1数据，采用景观分层的方法对非洲马里南部的农业用地进行分类，得到了与基础数据相类似的分类结果[7]。
Brown等利用时序MODIS植被指数对巴西MatoGrosso地区多年的农用地分类[8]。
然而，由于受空间分辨率、混合像元的影响，MODIS数据不能满足区域农作物分类以及精细提取的要求。
Landsat数据由于其较高的空间分辨率在区域土地利用/覆盖、农作物信息提取中得到应用。
Matejicek等利用Landsat数据对捷克西北部农用地的变化进行了研究[9]。
Vittek等利用LandsatMSS/TM数据对非洲西部土地利用变化进行了监测[10]。
但Landsat回访周期（16d）较长，加之阴雨天气的影响，难以获得时间序列的Landsat数据。
因此，基于物候特征Landsat数据的大面积农作物种植区信息提取受到数据的极大限制。
遥感技术及数据时空融合技术的发展，为获得时间序列Landsat数据提供有效的技术途径。
Gao等针对Landsat与MODIS数据的时空融合问题，提出了STARFM（SpatialAndTemporalAdaptiveReflectanceFusionModel）模型，试验验证了模型的有效性[11]。
Zhu等针对STARFM模型的不足，改进并提出ESTARFM（EnhancedSpatialAndTemporalAdaptiveReflectanceFusionModel）模型[12]。
邬明权等提出了基于混合像元分解的方法STDFM（SpatialAndTemporalDataFusionModel）用来融合MODIS与TM数据[13]。
Zhang等针对STDFM算法进行改进，提出了ESTDFM（EnhancedSpatialAndTemporalDataFusionModel）模型[14]。
对于上述几种模型在NDVI数据融合的效果方面，石月婵等以甘肃张掖市为例进行了对比分析，研究结果表明上述几种模型在NDVI数据融合效果方面基本相当[15]。
在数据数量要求方面STARFM模型具有一定的优势，除STARFM模型外，其他几种算法都需要在模型中输入2期Landsat数据。
STARFM模型提出后得到了较为广泛的应用，Hilker等利用STARFM模型融合得到时间序列的Landsat数据，融合Landsat数据的NDVI值也能够很好地反映研究区作物的物候特征[16]。
Walker等利用STARFM模型对Landsat数据和MODIS数据进行融合，以此分析了干旱地区森林的物候特征[17]。
Jia等利用STARFM模型融合的时序LandsatNDVI数据，对北京市土地覆盖进行分类制图研究[18]。
然而，基于融合MODIS与Landsat数据的水稻种植面积提取研究不多[19]，本文利用STARFM模型融合MODIS与Landsat数据对水稻种植面积进行提取具有一定意义。
洞庭湖区属于亚热带季风气候，云雨天气较多，难以获取水稻生长期的时间序列Landsat数据。
遥感数据时空融合技术为解决水稻生长期Landsat数据缺失问题提供了有效的技术途径。
为此，本文利用STARFM模型，融合不同类型数据，对研究区水稻种植面积进行提取。
1研究区及数据
1.1研究区概况
洞庭湖区位于长江中游荆江段南岸，跨湘、鄂两省，介于28°30′-29°31′N，111°40′-113°10′E之间。
湖区大部分海拔低于50m，属亚热带季风气候，年平均气温为15.8～17.4℃，年降水量在1200～1500mm。
本文选择的研究范围为环绕洞庭湖水域的丘陵和冲击平原地区，位于湖南省行政区内的洞庭湖区，地处湖南省东北部（如图1所示）。
研究范围涉及岳阳、常德、益阳3个市21个县（市、区），土地面积为4.56×104平方公里，占湖南省总面积的12.2%。
洞庭湖区是长江流域重要的商品粮基地，主要粮食作物为水稻（双季稻、一季稻），同期主要农作物为棉花。
图1研究区范围
Fig.1Studyarea
1.2数据与处理
本文所需遥感影像数据Landsat8OLI与MODIS13Q1均下载于USGS。
研究区横跨四幅Landsat8影像，编号分别为123/39、123/40、124/39、124/40；所需MODIS数据跨h27v05、h27v06两幅影像（如表1）。
表1遥感数据类型及获取日期
Table1Remotesensingdatatypesandacquisitiondata
数据类型
Datatypes
Landsat8OLI
MODIS13Q1
行列号
andnumberofRow/line
123/39
123/40
124/39
124/40
h27v05
h27v06
获取日期
Acquisitiondate
2013-05-12
2013-05-12
2013-08-7
2013-08-7
全年共23期
Atotalof23scenesofayear
全年共23期
Atotalof23scenesof
ayear
2013-06-13
2013-06-13
2013-10-10
2013-10-10
2013-07-31
2013-07-31
2013-09-17
2013-09-17
下载的Landsat8OLI数据分辨率为30m，质量较好，云覆盖均小于5%。
利用ENVI5.0图像处理软件的FLAASH模块进行了辐射定标与大气校正，以1：
50000的地形图为参考，采用二次多项式方法进行了图像几何校正。
将123/39与123/40、124/39与124/40分别进行影像镶嵌，得到123行4期、124行2期遥感影像，之后进行NDVI计算，并将NDVI值的取值范围设置为0～10000。
MODISNDVI产品数据MOD13Q1获取时间范围为2013-1-1～2013-12-19，共23期，分辨率为250m。
首先将MOD13Q1产品数据坐标系转换为高斯克吕格投影，与进行几何校正后的Landsa8OLI一致，并将Landsat8OLI与MODIS13Q1进行配准。
之后进行去无效值处理，并重采样使其空间分辨率与和Landsat8数据相一致（30m）。
以123行、124行Landsat8数据形成掩膜，分别对MODIS13Q1数据进行裁剪。
地面验证数据包括2013年研究区实地考察采样数据、2013年湖南省土地利用现状图及谷歌地球数据，验证样本像元数7597个，验证样本像元涵盖该研究所需的所有覆被类型，且较均匀分布于整个研究区。
2研究方法
2.1Landsat与MODIS数据时空融合
Landsat与MODIS数据的融合采用时空自适应反射率融合模型（STARFM），该模型基于t0时刻Landsat、MODIS数据，t1时刻的MODIS数据，结合不同的空间权重融合计算出t1时刻的Landsat数据[11,17]，STARFM模型的表达式如下：
（1）
式中：
L和M分别表示Landsat与MODIS像元反射率；W权重函数，权重函数决定了滑动窗口内各像元对预测值的贡献大小，利用光谱距离、时间距离与空间距离来确定函数权重；（xi,yj,t0）表示t0时刻位置为（xi,yj）处的像元；（xi,yj,t1）表示t1时刻位置为（xi,yj）处的像元；（xw/2,yw/2,t1）表示t1时刻移动窗口的中心像元。
本文利用时空自适应反射率融合模型（STARFM），由t0时刻的LandsatNDVI、MODISNDVI数据，以及16d间隔时间序列的MODISNDVI数据，通过式
（1）融合得出16d间隔时序的LandsatNDVI数据。
由于受天气及卫星回访周期的影响，只能得到123行4期、124行2期Landsat8数据。
利用STARFM模型进行影像融合时所需Landsat影像只需一期，且进行数据融合时尽量选择与目标日期相近的LandsatNDVI数据来预测目标日期的LandsatNDVI数据，以此来提高融合目标日期LandsatNDVI数据的效果。
将融合预测的123行与124行的LandsatNDVI数据进行镶嵌，并用研究区矢量范围图裁剪得到23期16d间隔时序的研究区LandsatNDVI影像数据。
2.2NDVI重构
为消除由云污染和大气变化引起的低值突变噪声对时序LandsatNDVI数据的影响，本文采用以Savitzky-Golay滤波法（S-G）、非对称高斯函数拟合法（AG,AsymmetricGaussianfunctionfittingmethod）及双Logistic曲线拟合法（D-L,DoubleLogisticfunctionfittingmethod）为内核[17-18]的Timesat软件包来对时间序列NDVI数据进行重构，并对3种曲线重构方法保持原始NDVI时间序列曲线整体特征及处理未受噪声污染点真实值的水平（保真性特征）进行比较，相关系数r表示2个样本组数间的相关程度，可以反映拟合后的NDVI时间序列曲线的保真性。
并且利用回归估计标准差来描述重建后的NDVI时间序列与原始值之间的平均差异程度，其值越小，拟合值的代表性越强。
Savitzky-Golay滤波方法首先根据云状态（像元可信度）对NDVI数列进行线性插值，然后利用S-G滤波器得到插值后曲线的模拟长趋势线，再根据上包络线得到新的NDVI曲线，将上述过程进行数次迭代并设置拟合影响系数作为迭代退出条件，最终得到较为平滑又能反映NDVI数值变化趋势的时间序列曲线。
其滤波过程可用公式
（2）表示：
（2）
式中：
Y是指NDVI原始值；
是NDVI拟合值；j是原始NDVI数组的系数；
是第i个NDVI值的卷积系数；N是滑动数组宽度（2m+1）；m是待拟合点左右两端各需的点数[20]。
非对称高斯函数拟合法是一种从局部拟合到整体拟合的方法，使用分段高斯函数来模拟植被生长过程，最终通过平滑连接各高斯拟合曲线时间序列重构[21]。
其中局部拟合公式为：
（3）
其中，
（4）
式中：
t表示t时刻的NDVI值；
和
为控制曲线的基准和幅度；
决定峰值和谷值的位置；
、
、
、
为控制曲线左、右半部分的宽度和陡峭度。
整体拟合函数为：
（5）
式中：
是NDVI的变化区间，
、
、
分别表示
区间内左边谷值、中间峰值及右边谷值所对应的局部函数；
、
为位于
之间的剪切系数。
双Logistic曲线拟合法是一种半局部拟合方法，其局部拟合方式与非对称高斯拟合方法类似（如式
（2））。
采用双Logistic函数（如式（6）），基于整体拟合函数（式（4））将各局部拟合函数的特征加以综合，重建新的NDVI时间序列曲线[22-24]。
（6）
2.3可分离性分析及最佳LandsatNDVI组合选择
J-M距离是基于特征计算不同类别样本间的距离，用来衡量类别间分离度的有效工具。
相对于欧式距离、巴氏距离等地表特征可分性判定方法，J-M距离更优[25]。
基于某一特征的2类样本的J-M距离计算公式如下：
（7）
式中：
B表示在某一特征维上的巴氏距离。
在样本对象满足正态分布的前提下，不同2类别间样本对象的巴氏距离（Bhattacharyyadistance,B）为：
（8）
式中，
表示某类特征的均值；
表示某类特征的方差，其中（k=1，2）。
表2选择计算J-M距离的LandsatNDVI数据
Table2LandatNDVIdataforJ-Mdistancecalculation
NDVI数据序号
NumberofNDVIdata
NDVI数据日期
DateofNDVIdata
一年中的天数
Dayofyear（DOY）
1
2013-04-23
113
2
2013-05-9
129
3
2013-05-25
145
4
2013-06-10
161
5
2013-06-26
177
6
2013-07-12
193
7
2013-07-28
209
8
2013-08-13
225
9
2013-08-29
241
10
2013-09-14
257
11
2013-09-30
273
12
2013-10-16
289
13
2013-11-1
305
J值在0～2之间，其大小代表样本间可分离程度。
当J=2时，表明2类在所选分类特征下完全分离；当J值较小时，表明分离性较差且会有较大数量的错分对象[26-28]。
本文结合Googleearth、湖南省土地利用现状图及部分实地调查数据，对林地、双季稻、一季稻、棉花等植被覆盖的J-M距离进行计算，选择最佳的LandsatNDVI日期组合。
结合研究区水稻生长期，确定用于J-M距离计算的LandsatNDVI数据及日期（表2）。
2.4基于时空特征的遥感分类与精度评价
支持向量机（SVM）可以自动寻找对分类有较大区分能力的支持向量，构造出类与类之间的间隔最大化的分类器，具有较高的分类准确率，因而在遥感分类中得到应用。
本文基于融合的时间序列LandsatNDVI影像，采用支持向量机（SVM）遥感影像分类方法，进行水稻种植区域信息的提取。
根据研究区具体情况，将研究区土地利用/覆盖划分为水体、双季稻、一季稻、棉花、林地和其他类共6类。
结合Googleearth、湖南省土地利用现状图及部分实地考察数据，在研究区内随机选取训练样本并利用支持向量机（SVM）进行分类（如表3）。
利用上文中地面验证数据生成验证感兴趣区并结合分类后数据生成混淆矩阵，采用总体精度（overallaccuracy）、生产者精度（produceraccuracy）、用户精度（useraccuracy）以及Kappa系数进行定量评价。
表3不同土地利用/覆盖类型训练样本和验证样本感兴趣区（ROI）及像元个数
Table3NumberofROIsandpixelsofdifferentlanduse/covertypesfortrainingandvalidating
ROI和像元个数
NumberofROIandpixels
水体
Water
双季稻
Doublecroppingrice
一季稻
Singleseasonrice
棉花
Cotton
林地
Forest
其他类
Others
训练样本
Trainingsamples
ROI个数
NumberofROIs
196
102
74
83
212
85
像元个数
Numberofpixels
3256
4585
2168
3749
4053
3821
验证样本Validationsamples
ROI个数
NumberofROIs
48
52
38
51
65
43
像元个数
Numberofpixels
1048
1204
850
1426
1987
967
3结果与分析
3.1遥感数据时空融合结果
利用公式
（1）融合预测出整个研究区23期LandsatNDVI数据，由于图像数量较多，本文只展示了研究区水稻生长关键期的融合LandsatNDVI数据，如图2，日期分别为2013-05-25（第145天）、2013-07-12（第193天）、2013-07-28（第209天）、2013-08-29（第241天）、2013-09-14（第257天）、2013-10-16（第289天）。
a．第145天，双季早稻抽穗b．第193天，双季早稻收割c．第209天，一季稻分蘖、抽穗
a．DOY145,Doublecroppingriceb．DOY193,Doublecroppingricec．DOY209,Singleseasonrice
（early）washeading（early）washarvestwastilleringandheading
d．第241天，双季晚稻抽穗e．第257天，一季稻收割f．第289天，双季晚稻收割
d．DOY241,Doublecroppingricee．DOY257,Singleseasonricef．DOY289,Doublecroppingrice
（late）washeadingwasharvest（late）washarvest
图2研究区水稻生长关键期的部分融合LandsatNDVI数据
Fig.2SomeoffusingLandsatNDVIdataduringcriticalstageofriceofstudyarea
3.2融合时序LandsatNDVI滤波结果
通过三种NDVI时间序列曲线重构方法及其保真性特征（相关系数）比较，得出S-G滤波法的整体保真性较好。
除双季稻外，其他3种植被类型时序NDVI经S-G滤波法重构后的保真性较其他两种算法要好（如图3a）。
以局部拟合为主的S-G滤波有较强的细节拟合能力，而AG与D-L拟合方法以曲线上包络线吻合为主要特征，在去噪的同时也会导致重建的NDVI偏离真实值。
与此同时，对各类型植被用不同方法重建的效果分析发现，除林地类型S-G滤波方法的重建前后差异程度（回归估计标准差）相对较小，即用该方法所得的拟合值整体代表性较好（3b）。
因此本文将选用经过S-G滤波后的时间序列NDVI数据进行后续研究。
a．三种算法下不同覆被类型NDVI重构后原始值与拟合值的相关系数
a．CorrelationcoefficientbetweenoriginalNDVIvalueandfittingvaluereconstructedbythreealgorithmsfordifferentlandcovertypes
b．三种算法下不同覆被类型NDVI重构后原始值与拟合值的回归估计标准差
b.RootmeansquareerrorbetweenoriginalNDVIvalueandfittingvaluereconstructedbythreealgorithmsfordifferentlandcovertypes
图3三种算法下不同覆被类型NDVI重构保真性比较
Fig.3ComparisonoftheabilityonkeepingthemaincharactersafterNDVIreconstructionbythreealgorithmsfordifferentlandcovertypes
拟合后的双季稻NDVI值呈现明显的“双峰”型，如图4。
在5月下旬（第145天）出现第1个峰值（0.69），此时早稻正处于抽穗阶段，故其NDVI值较高。
到7月上旬（第193天）出现谷值，因此时为早稻收割期，NDVI值会有所减小。
8月下旬（第241天）左右出现了第2个峰值（0.73），这段时期由于晚稻处于孕穗、抽穗状态，NDVI值相对较高。
10月上旬（第289天）晚稻相继收割，故NDVI出现第2个谷值。
后期NDVI值有上升趋势主要是由田间杂草或绿色肥料（满江红等），油菜的栽种等所导致的。
棉花的NDVI值呈“单峰”型。
6月下旬到7月上旬（第177-193天）棉花处于蕾期，NDVI值达到最大值（0.73）。
11月中旬之后，棉花相继落叶，故其NDVI值不断减小。
7月下旬至8月上旬（第193-209天）一季稻处于分蘖、孕穗期，NDVI值达到最大值0.72。
到8月下旬至9月上旬（第241-257天），一季稻收获，NDVI值较低。
之后杂草等会导致NDVI值小幅度升高。
由于研究区大都为常绿阔叶林，故林地的NDVI值表现最为平稳，除了2月至3月上旬由于气候影像，其NDVI值会稍有降低，林地NDVI值全年基本维持在0.8左右。
4种难以辨别覆被类型的NDVI值会随着时间发生不同的变化，这种NDVI的变化将有助于区分这四种覆被类型。
例如，5月下旬、7月上旬、9月上旬可以将早稻与一季稻进行区分，11月上旬可以区分棉花和水稻、棉花等。
图4LandsatNDVI拟合结果
Fig.4FittingresultofLandsatNDVIdata
3.3J-M距离及最佳LandsatNDVI组合
考虑到早稻移栽日期大概在四月中旬，以及晚稻收割大致在10月上旬，选择进行J-M距离计算的LandsatNDVI的日期为2013年4月23日～2013年11月1日（第113～305天），共13期。
根据公式（7）和（8），并结合滤波后的LandsatNDVI，按不同组合计算J-M距离。
进行J-M距离计算时，通过选取不同的地面验证数据进行对比分析，得到各植被类型间最佳的J-M距离（如表4）。
升级会员 