Print
结合像元分解和STARFM模型的遥感数据融合
expand article info 谢登峰 , 张锦水 , 孙佩军 , 潘耀忠 , 云雅 , 袁周米琪
北京师范大学地表过程与资源生态国家重点实验室, 北京师范大学资源学院, 北京  100875

摘要

高空间、时间分辨率遥感数据在监测地表快速变化方面具有重要的作用。然而,对于特定传感器获取的遥感影像在空间分辨率和时间分辨率上存在不可调和的矛盾,遥感数据时空融合技术是解决这一矛盾的有效方法。本文利用像元分解降尺方法(Downscaling mixed pixel)和STARFM模型(Spatial and Temporal Adaptive Reflectance Fusion Model)相结合的CDSTARFM算法(Combination of Downscaling Mixed Pixel Algorithm and Spatial and Temporal Adaptive Reflectance Fusion Model)进行遥感数据融合。首先,利用像元分解降尺度方法对参与融合的MODIS数据进行分解降尺度处理;其次,利用分解降尺度的MODIS数据替代STARFM模型中直接重采样的MODIS数据进行数据融合;最后以Landsat 8和MODIS遥感影像数据对该方法进行了实验。结果表明:(1) CDSTARFM算法比STARFM和像元分解降尺度算法具有更高的融合精度;(2)CDSTARFM能够在较小的窗口下获得更高的融合精度,在相同的窗口下其融合精度也高于STARFM;(3)CDSTARFM融合的影像更接近真实影像,消除了像元分解降尺度影像中的"图斑"和STARFM模型融合影像中的"MODIS像元边界"。

关键词

像元分解;降尺度;STARFM;遥感;数据融合;CDSTARFM

1 引 言

高空间、时间分辨率的遥感数据在监测土地覆盖变化和作物生长以及物候参数反演等方面具有重要的作用(Zhang等,2013)。然而,目前获取到的卫星传感器数据存在空间分辨率、时间分辨率、光谱分辨率相互矛盾,即同时满足高空间、高时间和高光谱分辨率的卫星遥感数据是不现实的(Price,1994Emelyanova等,2013)。比如,L and sat卫星的多光谱影像空间分辨率是30 m,其适中的空间分辨率和易获取性在植被指数提取、监测土地覆盖动态变化和生态系统方面具有广泛的应用(Rees等,2003Masek等,2008Schroeder等,2011),但因其16天的重访周期和阴雨天气的影响,难以保证获得连续有效的影像和对农作物生长、物候变化和自然灾害进行有效的监测(Ju和Roy,2008González-Sanpedro等,2008);搭载在Terra/Aqua卫星上的中分辨率成像光谱仪MODIS能提供时间分辨率为1天、空间分辨率为250—1000 m的遥感影像,被广泛应用于大尺度区域的地表覆盖类型监测(Arvor等,2011Notarnicola等,2013Zhou等,2013),但MODIS 数据的低空间分辨率限制其在景观破碎和异质性较强区域的应用(Shabanov等,2003)。若能综合中高分辨率影像(如L and sat OLI)的高分辨率特性和MODIS时间分辨率的优势,将会大幅度提高植被遥感监测能力(Lunetta等,1998Arai等,2011)。

近年来,学者提出了一系列获取高时空分辨率数据的遥感数据融合模型。这些融合模型总体分为两大类:(1)像元分解降尺度算法(downscaling mixed pixel),该算法是基于线性光谱混合模型,通过从高分辨率数据中提取地物类别和丰度来分解低分辨率的像元获得类别的光谱值(Zhukov等,1999Keshava和Mustard,2002Zurita-Milla等,2008Bioucas-Dias等,2012)。该算法简单易行,不需要中、高分数据具有相同的波段,只需要高分辨率的地表覆盖类型数据就可以进行数据融合(Zurita-Milla等,2011),但是该算法获得的类别反射率是整个区域或局部区域的平均反射率,不能很好地体现地物光谱的空间差异性,会出现“图斑”现象(Busetto等,2008Wu等,2012邬明权等,2012);(2)自适应遥感图像时空融合方法STARFM(Spatial and Temporal Adaptive Reflectance Fusion Model)(Gao等,2006),该模型用于融合L and sat和MODIS数据构建每日的L and sat数据,在融合过程中不仅考虑了空间的差异性,而且也考虑了时间的差异性,是目前应用最广泛的时空融合模型之一(Hilker等,2009a2009bSingh,2011Emelyanova等,2013)。但STARFM模型存在两个问题:(1)若短暂或突变的地表变化信息没有被基期的L and sat影像记录下,那么融合的影像就不能捕捉到该地表变化信息(Gao等,2006Hilker等,2009a);(2)该模型构建的数据质量依赖于从MODIS数据中提取的纯净像元的时间变化信息,若在窗口内没有纯净的MODIS 像元,数据融合质量在一定程度上会降低(Gao等,2006Hilker等,2009b)。为解决问题(1)的不足,Hilker等人(2009a)基于STARFM模型提出了一种新的融合算法STAARCH(Spatial Temporal Adaptive Algorithm for Mapping Reflectance Change),分别从L and sat和MODIS数据中提取空间的变化和时间变化,通过选取最佳的基准期L and sat影像来提高融合算法的精度。对于问题(2),由于MODIS影像空间分辨率低,很难在地块破碎、异质性强的区域内找到纯净的MODIS像元,为了解决纯净MODIS 像元对STARFM算法精度的影响,Zhu等人(2010)提出了适合异质性区域的增强型STARFM(ESTARFM)模型,ESTARFM模型通过假设一段时间内地物反射率为线性变化以及混合像元值是由不同地物类别光谱值的线性组合,在模型中引入一个转换系数来提高破碎区域数据融合的精度,但ESTARFM模型至少需要两对同期的L and sat和MODIS基期影像,与STARFM相比增加了数据获取的难度(Zhu等,2010)。

针对STARFM模型因在破碎区域难搜寻到纯净像元而降低融合数据精度的问题,可以首先利用像元分解降尺度算法对低分辨率数据进行降尺度处理,然后将降尺度数据取代STARFM算法中重采样的低分辨率数据,最后再利用STARFM算法后续步骤进行数据融合。本文提出像元分解降尺度算法结合STARFM模型的CDSTARFM算法,利用L and sat 8和MODIS反射率数据构建具有L and sat 8 OLI空间分辨率(30 m)和MODIS时间分辨率(1 d)的高时空遥感数据,并验证该方法的适用性。

2 CDSTARFM介绍

CDSTARFM算法结合像元分解降尺度算法和STARFM各自的优势,以期解决像元分解降尺度融合数据出现的“图斑”和在一定程度上增加STARFM算法搜寻到纯净像元的概率,提高数据融合的精度。输入的数据包括一对基期(tk)的中、低分辨率遥感影像和预测期的低分辨率影像,输出的是预测时期的中分辨率影像。CDSTARFM算法流程见图 1,具体实现步骤如下。

图1 CDSTARFM算法流程图
Figure 1 Algorithm flow of CDSTARFM

2.1 类别定义及丰度提取

混合像元内的地物类别及其丰度是混合像元分解的主要参数。类别是通过传统方法聚类中分辨遥感影像获取,丰度是通过式(1)提取每个混合像元内各类别的丰度,并假设类别和丰度在一定时间内不随时间发生变化。利用传统非监督分类方法ISO-DATA对中分辨率遥感影像聚类获得分类图像,利用ArcGIS在分类图上创建大小为低分辨率像元尺度的网格提取格网内各类别的丰度(邬明权等,2012)。

$ {f_c}\left( {i,c} \right) = \frac{Q}{S} $ (1)

式中,fc(i,c)表示i位置低分辨率像元内c类地物的丰度,Q表示低分辨率像元内c类地物的像元数,S表示低分辨率像元内所有类别地物的像元数。

2.2 降尺度数据获取

混合像元分解的基础是线性光谱混合模型,在忽略不同传感器数据之间的地理配准误差和大气校正误差的条件下,混合像元的光谱值等于该像元内地物类别的光谱值与其丰度的线性组合(Settle和Drake,1993)(式(2))。然而,由于空间异质性的存在和最小二乘算法理论的需要,式(2)需要在一定大小的窗口内进行,但过大的窗口会削弱空间的异质性。同时,窗口的使用既可以消除中、高分辨率影像之间因配准误差带来的影响,又可以体现同类地物的空间差异性(Wu等,2012Zhang等,2013)。因此,选择合适的矩形窗口进行混合像元分解是必要的。为选择一个合适窗口,设置不同大小的窗口进行混合像元分解,对各窗口下降尺度的数据与真实数据进行比较,通过定量精度指标选择合适的窗口。在一定大小窗口下,利用式(1)求得窗口内每个混合像元内各类别丰度,组成到n2k列的丰度矩阵A(n2×k)(n2>k)(n表示矩形窗口大小,k表示窗口内类别数),利用最小二乘法(式(3))求解窗口内各类别的光谱值(Busetto等,2008Wu等,2012Zhang等,2013)。然后,把求得的值依照类别类型对应到窗口内相应像元的位置上,获得分解后的降尺度数据。

$ R\left( {i,t} \right) = \sum\limits_{c = 0}^k {{f_c}\left( {i,c} \right) \times \bar r\left( {c,t} \right) + \xi \left( {i,t} \right)} $ (2)

$ \left[ {\begin{array}{*{20}{c}} {R\left( {1,t} \right)}\\ \cdots \\ {R\left( {i,t} \right)}\\ \cdots \\ {R\left( {n,t} \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{f_c}\left( {1,1} \right)}& \cdots &{{f_c}\left( {1,c} \right)}& \cdots &{{f_c}\left( {1,k} \right)}\\ \cdots &{}& \cdots &{}& \cdots \\ {{f_c}\left( {i,1} \right)}& \cdots &{{f_c}\left( {i,c} \right)}& \cdots &{{f_c}\left( {i,k} \right)}\\ \cdots &{}& \cdots &{}& \cdots \\ {{f_c}\left( {n,1} \right)}& \cdots &{{f_c}\left( {n,c} \right)}& \cdots &{{f_c}\left( {n,k} \right)} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\bar r\left( {1,t} \right)}\\ \cdots \\ {\bar r\left( {c,t} \right)}\\ \cdots \\ {\bar r\left( {k,t} \right)} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\xi \left( {1,t} \right)}\\ \cdots \\ {\xi \left( {c,t} \right)}\\ \cdots \\ {\xi \left( {k,t} \right)} \end{array}} \right] $ (3)

约束条件: $ \sum\limits_{c = 0}^k {{f_c}\left({i,c} \right)= 1;0 < } \bar r\left({c,t} \right)< 1 $

式中,R(i,t)t时期i位置低分辨率像元的反射率,fc(i,c)i位置低分辨率像元内类别c占该像元的面积比,(c,t)t时期类别c的平均反射率,ξ(i,t)为残差,k为窗口内类别数。

2.3 融合影像生成

CDSTARFM利用t0tk两期低分辨率影像降尺度的数据和tk时期中分辨率遥感影像预测t0时期中分辨率遥感数据。STARFM算法实现的关键在于如何权重空间信息,且这个权重算法可以自适应于不同的区域(Gao等,2006)。Gao等人(2006)通过在一定窗口内筛选与中心像元光谱相似的相似性像元,利用相似性像元的光谱差异性,时间差异性和与中心像元的相对距离来加权计算预测期的中心像元值。由于CDSTARFM继承了STARFM算法,其实现过程也是通过加权窗口内相似性像元计算预测期的中心像元值。

窗口内中心像元的相似性像元由tk时期中分辨率影像各波段的阈值θb确定,窗口内与中心像元差值的绝对值小于各波段θb的像元作为相似性像元。

$ {\theta _b} = \sqrt {\frac{1}{N}{{\sum\limits_{i = 1}^N {\left( {{x_i} - \mu } \right)} }^2}} /m $ (4)

式中,θb表示窗口内b波段的阈值,N为窗口内中分辨率像元个数,xii位置的像元反射率,μ为窗口内像元反射率均值,m为类别数(Gao等,2006Zhu等,2010)。

相似像元的权重大小Wijk有3个指标来衡量:(1)光谱的差异性Sijk;(2)时间的差异性Tijk;(3)相似性像元与中心像元的相对距离Dijk

$ {S_{ijk}} = \left| {L\left( {{x_i},{y_i},{t_k}} \right) - M\left( {{x_i},{y_i},{t_k}} \right)} \right| $ (5)

$ {T_{ijk}} = \left| {M\left( {{x_i},{y_i},{t_k}} \right) - M\left( {{x_i},{y_i},{t_0}} \right)} \right| $ (6)

$ {D_{ijk}} = 1 + \sqrt {{{\left( {{x_{w/2}} - {x_i}} \right)}^2} + {{\left( {{y_{w/2}} - {y_j}} \right)}^2}} /A $ (7)

式中,A为表示空间距离相对于光谱差异性和时间差异性重要性的常量。A=750 m(Gao等,2006)

利用式(8)(9)加权窗口内相似像元的反射率计算预测期的中心像元反射率,生成最终的融合影像。

$ \begin{array}{l}L\left({{x_{w/2}},{y_{w/2}},{t_0}} \right)= \\\sum\limits_{i = 1}^w {\sum\limits_{j = 1}^w {\sum\limits_{k = 1}^w {{W_{ijk}} \times \left({M\left({{x_i},{y_i},{t_0}} \right)+ L\left({{x_i},{y_i},{t_k}} \right)- M\left({{x_i},{y_i},{t_k}} \right)} \right)} } } \end{array} $ (8)

$ \begin{array}{l}{W_{ijk}} = \left({\frac{1}{{{S_{ijk}} \times {T_{ijk}} \times {D_{ijk}}}}} \right)/\\\sum\limits_{i = 1}^w {\sum\limits_{j = 1}^w {\sum\limits_{k = 1}^w {\left({\frac{1}{{{S_{ijk}} \times {T_{ijk}} \times {D_{ijk}}}}} \right)} } } \end{array} $ (9)

式(5)—(9)中,w为窗口大小,L(xi,yj,tk)M(xi,yj,tk)分别为tk时期给定位置(xi,yj)的中分辨率数据和降尺度数据的像元值,(xw/2,yw/2)为窗口的中心像元,光谱差异性Sijk值越小说明给定位置与邻近像元的光谱相似度越高,赋予较高的权重;时间差异性Tijk值越小说明该段时间内光谱变化越小,赋予较高的权重;相对距离Dijk值越小,赋予较高的权重(Gao等,2006)。

3 算法测试

3.1 测试数据与预处理

测试数据覆盖的区域(37.52°N— 37.73°N,115.44°E—115.71°E)位于河北省衡水市境内,面积约为25×25 km2(50×50个MODIS像元)。该区域以农业种植为主,作物类型以玉米为主,兼有棉花、大豆和花生等,影像右上区域为居民区,中间有一定面积的水体,地表类型复杂,在MODIS像元尺度下,纯净像元较少,满足本研究需要。

测试数据(表 1)为2景2014-08-19和2014-09-04 L and sat 8(OLI)数据,以及2景2014-08-19和2014-09-04 MODIS地表日反射率产品数据(MOD09GA),实验中选用绿、红和近红外波段进行算法测试(表 2)。测试数据来源于USGS网站。L and sat 8影像的坐标系为UTM-WGS84 Zone 50N,由于L and sat 8数据已进行过基于地形数据几何校正,可以直接使用。利用ENVI5.1提供的定标和FLAASH大气校正模块把DN值转换为地表真实反射率(Jia等,2014)。日反射率的MOD09GA标准产品数据采用的是HDF格式和正弦(Sinusoidal)坐标系统,由于该数据是经过大气校正和几何校正的2级标准产品数据并附带有详细的地理坐标信息,可以利用MODIS 重投影工具MRT(MODIS Re-projection Tool)转换为Geo-tif数据格式和L and sat数据的UTM-WGS84坐标系,并用最近邻域法重采样像元为480 m,便于后续混合像元分解(Wu等,2012)。

表 1 Landsat 8和MODIS数据主要特征
Table 1 Main information of Landsat 8 and MODIS data 下载CSV 
数据类型 获取时间 轨道号(Path/Row) 数据用途
Landsat 8(OLI) 2014-08-19 123/034 分类和筛选相似性像元
2014-09-04 精度评价
MOD09GA 2014-08-19 h27/v05 像元分解降尺度
2014-09-04 像元分解降尺度
表 2 Landsat 8和MOD09GA对应波段信息
Table 2 Band information of Landsat 8 and MOD09GA data 下载CSV 
数据波段波长/nm空间分辨率/m
Landsat 8 OLI Band 3530—59030
Landsat 8 OLI Band 4640—67030
Landsat 8 OLI Band 5850—88030
MOD09GA Band 4545—565500
MOD09GA Band 1620—670500
MOD09GA Band 2841—876500

3.2 实现流程

首先,利用2014-08-19的L and sat 8多光谱数据进行ISODATA聚类,子类合并后类别为10类,结合聚类数据和最小二乘算法分别在窗口大小w(MODIS像元)为5×5、7×7、11×11、15×15、21×21、31×31和41×41对2014-08-19的MODIS数据进行像元分解。利用分解后的降尺度数据与真实L and sat 8数据之间的相关系数r、偏差Bias、均方根误差RMSE和ERGAS(the Erreur Relative Globale Adimensionalle de Synthèse)定量指标对降尺度数据进行精度评价,选出最适的窗口大小。其次,在该最佳窗口下对2014-09-04的MODIS数据进行像元分解降尺度。最后,利用CDSTARFM算法和2014-08-19和2014-09-04两期最佳窗口下的降尺度数据以及基期2014-08-19 L and sat 8数据进行数据融合。CDSTARFM 算法的搜寻窗口w(L and sat OLI像元)设为7×7、11×11、31×31、61×61、101×101和151×151。为了便于该算法与STARFM算法比较,STARFM算法利用与CDSTARFM算法相同的参数进行数据融合。通过融合数据与真实数据之间的r、Bias、RMSE和ERGAS等指标对融合数据进行精度评价(Amorós-López等,2013Gevaert和García-Haro,2015)。

$r = \frac{{\sum\limits_{j = 1}^N {\left( {{x_j} - \bar x} \right)\left( {{y_j} - \bar y} \right)} }}{{\sqrt {\sum\limits_{j = 1}^N {{{\left( {{x_j} - \bar x} \right)}^2}} \sum\limits_{j = 1}^N {{{\left( {{y_j} - \bar x} \right)}^1}} } }}$ (10)

$ Bais = \frac{1}{N}\sum\limits_{j = 1}^N {\left({{x_j} - {y_i}} \right)} $ (11)

$ RMSE = \sqrt {\frac{{{{\sum\limits_{j = 1}^N {\left({{x_j} - {y_i}} \right)} }^2}}}{N}} $ (12)

式中,N表示像元个数,j表示融合数据和真实数据的像元位置,xjyj分别表示j位置融合数据和真实数据的像元值,$ {\bar x} $和$ {\bar y} $分别表示融合数据和真实数据的平均值。

$ ERGAS = 100\frac{H}{L}\sqrt {\frac{1}{{{N_{ban}}}}\sum\limits_{i = 1}^{{N_{ban}}} {{{\left({RMS{E_i}/{M_i}} \right)}^2}} } $ (13)

式中,H表示高分辨率影像的像元大小(30 m),L表示低分辨影像的像元大小(480 m),Nban表示波段的数量,RMSEi表示i波段融合数据与真实数据对应波段的均方根误差,Mi表示真实L and sat 8数据i波段的平均值。ERGAS值越小说明融合数据与真实数据在空间上的差异就越小,融合数据的质量就越高(Amorós-López等,2013Gevaert和García-Haro,2015)。

4 结果分析

4.1 定量分析

表 3表 4表 5分别表示像元分解降尺度算法、STARFM模型和CDSTARFM算法在不同窗口下融合的Green、Red和NIR波段数据与同期真实L and sat 8数据对应波段的相关系数r、偏差Bias、均方根误差RMSE和ERGAS指标。

表 3 不同窗口2014-08-19和2014-09-04 MODIS降尺度数据的精度指标
Table 3 Accuracy indicators of downscaled data on August 19 and September 4, 2014 at different window size 下载CSV 
时间 窗口 w (n×n) r/% RMSE ERGAS Bias
Green Red NIR Green Red NIR Green Red NIR Green Red NIR
注:粗体数值表示实验中的最佳值。
2014-08-19 5 51.32 58.71 79.73 0.0261 0.0310 0.0623 2.6171 3.5138 1.2301 0.0106 0.0055 0.0074
7 64.26 69.94 85.77 0.0218 0.0253 0.0511 2.1808 2.8769 1.0087 0.0107 0.0056 0.0072
11 78.42 81.98 89.76 0.0186 0.0205 0.0436 1.8602 2.3223 0.8598 0.0106 0.0055 0.0070
15 82.12 86.02 91.56 0.0179 0.0189 0.0399 1.7921 2.1475 0.7883 0.0106 0.0056 0.0067
21 84.73 88.68 92.14 0.0174 0.0179 0.0387 1.7420 2.0291 0.7635 0.0106 0.0054 0.0067
31 85.63 89.59 91.96 0.0171 0.0173 0.0391 1.7106 1.9628 0.7713 0.0104 0.0052 0.0068
41 84.86 89.47 92.11 0.0173 0.0174 0.0387 1.7308 1.9765 0.7643 0.0105 0.0054 0.0070
2014-09-04 31 82.18 85.79 89.95 0.0182 0.0222 0.0401 1.8052 2.3020 0.8737 0.0067 0.0033 0.0095
表 4 STARFM算法在不同窗口融合的2014-09-04数据的精度指标
Table 4 Accuracy indicators of prediction data by STARFM on September 4, 2014 at different window size 下载CSV 
窗口 w (n×n) r/% RMSE ERGAS Bias
Green Red NIR Green Red NIR Green Red NIR Green Red NIR
注:粗体数值表示实验中的最佳值。
7 88.22 89.26 94.16 0.0130 0.0172 0.0373 1.2823 1.8061 0.8088 -0.0020 -0.0033 0.0182
11 88.80 89.87 93.94 0.0130 0.0172 0.0351 1.2865 1.8129 0.7611 -0.0025 -0.0027 0.0134
31 88.95 90.00 94.89 0.0127 0.0170 0.0300 1.2559 1.7916 0.6507 -0.0011 -0.0016 0.0090
61 88.44 89.48 94.90 0.0131 0.0176 0.0302 1.2968 1.8498 0.6553 -0.0019 -0.0019 0.0103
101 88.04 89.31 94.74 0.0133 0.0177 0.0309 1.3117 1.8673 0.6702 -0.0014 -0.0017 0.0114
151 87.92 89.21 94.75 0.0133 0.0179 0.0310 1.3147 1.8782 0.6735 -0.0025 -0.0021 0.0119
表 5 CDSTARFM算法不同窗口融合的2014-09-04数据的精度指标
Table 5 Accuracy indicators of prediction data by CDSTARFM on September 4, 2014 at different window size 下载CSV 
窗口 w (n×n) r/% RMSE ERGAS Bias
Green Red NIR Green Red NIR Green Red NIR Green Red NIR
注:粗体数值表示实验中的最佳值。
7 91.16 92.26 96.00 0.0118 0.0151 0.0260 1.1678 1.5876 0.5654 -0.0030 -0.0013 0.0068
11 91.29 92.29 96.31 0.0116 0.0151 0.0249 1.1502 1.5850 0.5416 -0.0018 -0.0002 0.0067
31 91.21 91.92 96.50 0.0117 0.0154 0.0245 1.1550 1.6224 0.5317 -0.0022 -0.0008 0.0073
61 91.06 91.71 96.50 0.0117 0.0156 0.0245 1.1564 1.6437 0.5326 -0.0020 -0.0003 0.0076
101 90.94 91.58 96.50 0.0117 0.0158 0.0246 1.1615 1.6572 0.5334 -0.0019 -0.0002 0.0077
151 90.83 91.45 96.50 0.0118 0.0159 0.0246 1.1671 1.6700 0.5341 -0.0019 -0.0003 0.0078

表 3中2014-08-19 MODIS降尺度数据的各项评价指标可以看出,当滑动窗口w为31×31 个MODIS像元时,评价指标几乎均达到了最佳值(表中的下划线粗体数值)。因此,选择该窗口作为分解降尺度的最佳窗口,并在该窗口下对2014-09-04的MODIS数据进行降尺度处理,其分解降尺度数据的精度指标见表 3

表 4表 5中可以看出,当窗口大小为31×31个L and sat 8 OLI像元时,STARFM算法融合的各波段数据的r、Bias、RMSE及ERGAS指标几乎同时达到了最佳。STARFM算法的精度依赖于在搜寻窗口内能否找到纯净的MODIS像元,MODIS像元的纯度是通过同期同位置的L and sat像元与MODIS像元之间的光谱值差异来衡量,差异越小说明MODIS像元越纯净(Gao等,2006)。CDSTARFM算法融合数据的精度(表 5)在窗口为11×11 L and sat OLI像元时达到最佳,其最佳窗口大小明显减小,并且也提高了融合数据的精度。CDSTARFM算法融合的green、red和NIR各波段数据与STARFM算法和降尺度算法相比,其相关系数r分别平均提高了0.02和0.06; Bias分别平均提高了0.003和0.002; RMSE分别平均提高了0.003和0.008;ERGAS分别平均提高了0.2和0.6。这是因为CDSTARFM算法利用的MODIS降尺度数据比STARFM中重采样的MODIS数据更准确地反映地表信息,能够在较小的窗口内搜寻到纯净像元的概率增加,从而提高融合数据的精度。同时,相同窗口大小下,CDSTARFM算法的融合数据各波段的r、RMSE和ERGAS均高于STARFM算法。由此,CDSTARFM算法可以获得精度更高的融合数据。

图 2表示3种方法的最佳融合数据与真实L and sat 8数据的散点图。由于像元分解降尺度算法削弱了地物的空间差异性,得到的是窗口内地物类别的平均反射率,使得散点图上会出现“条带”。而STARFM和CDSTARFM算法是利用中心像元周围相似性像元的信息加权进行数据融合很好地体现了地物的空间差异性,不会出现 “条带”现象,对应波段的散点较好地分布在1∶1对角线两侧,r也明显高于像元分解降尺度算法。但与STARFM相比,CDSTARFM融合数据的r提高了约2%,其散点图也更集中在对角线两旁。由此可见,CDSTARFM算法融合的数据与真实数据之间具有更高的相关性。

图2 3种算法在最佳窗口融合的2014-09-04数据与真实数据的散点图
Figure 2 Scatterplots of reference images and the fusion images on September 4, 2014, predicted by the three methods

4.2 融合图像分析

像元降尺度算法、STARFM算法和CDSTARFM算法的最佳融合影像分别见图 3(c)(e)。3种方法融合的影像与2014-08-19和2014-09-04真实L and sat 8影像进行目视比较,图 3(c)(e)图 3(b)相比较可以看出,3种方法对面积较大的地物(水体、居民区和较宽的道路)均可清晰地识别。但从放大的局部细节图(图 3(h)(j))来看,在较破碎的区域,像元降尺度算法的融合影像(图 3(h))有“图斑”,不能识别较小的地物及地物内部的纹理信息,而图 3(i)(j)能够反映出细小地物和纹理特征;由于STARFM算法使用的是重采样的MODIS数据,使得同一MODIS像元内重采样的像元之间具有很强的均质性,相邻的MODIS像元之间在光谱上具有一定的差异,尤其在破碎区域其差异性更大,这是导致融合影像出现“MODIS像元边界”(图 3(i)中黑色框内所示)的主要原因,而CDSTARFM算法的融合影像可以有效消除MODIS边界。图 3(g)中黄色椭圆框是图 3(f)发生变化的地物,但3种融合方法均没有捕捉到该变化信息,这是因为用于融合的基期L and sat 8影像(图 3(f))没有体现该变化信息(Gao等,2006),由于CDSTARFM算法大部分继承了STARFM算法,所以也没有解决该问题;但Hilker等提出的STAARCH方法通过选择合适的基期L and sat数据解决了这个问题(Hilker等,2009a)。总体来说,CDSTARFM融合影像的目视效果更接近真实影像。

图3 3种方法在最佳窗口融合的2014-09-04影像与真实影像的比较(NIR-Red-Green组合)
Figure 3 Comparison of reference image and fusion images predicted by the downscaling method, STARFM, and CDSTARFM on September 4, 2014 at the best window size, respectively (NIR-Red-Green combination)

5 结论

时空融合技术是一种有效解决遥感数据的时间分辨率和空间分辨率相互矛盾的方法。本文利用像元分解降尺度算法和STARFM模型相结合的CDSTARFM算法融合中、低分辨率遥感数据,以L and sat 8和MODIS数据对CDSTARFM算法进行了实验,得出如下结论:

(1)CDSTARFM算法利用的MODIS降尺度数据比STARFM模型中的MODIS重采样数据,更真实地反应地表情况,能够增加“纯净”像元的数量,提高融合数据的精度。CDSTARFM算法融合的green、red和NIR各波段数据与STARFM算法和降尺度算法相比,其相关系数r分别平均提高了0.02和0.06; Bias分别平均提高了0.003和0.002;RMSE分别平均提高了0.003和0.008;ERGAS分别平均提高了0.2和0.6,其各波段的散点图也更为集中,与真实影像具有更高的相关性。

(2)CDSTARFM算法精度最高的融合窗口大小w(11×11个OLI像元)小于STARFM算法的最佳窗口w(31×31个OLI像元)。同时,在相同的窗口大小下,CDSTARFM算法融合数据的精度也均好于STARFM算法融合数据的精度。

(3)在破碎区域,CDSTARFM算法可以消除像元分解降尺度融合影像中出现的“图斑”和STARFM算法融合影像中出现的“MODIS像元边界”,使得CDSTARFM算法的融合影像与真实影像的目视效果最接近。

(4)CDSTARFM算法不仅适合L and ssat 8和MODIS数据的融合,还适合其他中/高分辨率数据(如高分一号16 m数据、SPOT-5数据)与MODIS数据融合构建高时空遥感数据。

另外,CDSTARFM算法精度在一定程度上依赖MODIS降尺度数据的精度,高精度的降尺度数据必定会提高CDSTARFM算法的精度。而本文实验仅使用最简单的混合像元分解模型,这对CDSTARFM算法的精度会产生一定影响。除此之外,CDSTARFM算法仅验证了窗口大小对融合数据的影响,没有考虑地表聚类数量对融合数据精度的影响以及 “极度”破碎区域(如中国南方地区)的地物类型、地块异质性等地表景观特征对该算法的影响;仅仅使用相关指标对融合数据进行评价,没有验证该算法的融合数据在反演地表参数(如NDVI)等方面的准确性和适用性,这些需要下一步的实验来验证。

参考文献

Remote sensing data fusion by combining STARFM and downscaling mixed pixel algorithm
expand article info XIE Dengfeng , ZHANG Jinshui , SUN Peijun , PAN Yaozhong , YUN Ya , YUAN Zhoumiqi
State Key Laboratory of Earth Surface Processes and Resource Ecology, College of Resources Science and Technology, Beijing Normal University, Beijing 100875, China

Abstract

High spatial and temporal resolution remote sensing data play an important role in monitoring the rapidly changing information regarding the earth's surface. However, the spatial and temporal resolutions of a remote sensing image for a specific sensor contain irreconcilable conflict. Fusing the spatial and temporal features of different remote sensing images to generate high-spatial-temporal remote sensing data is an effective method to solve the contradiction. This paper aims to improve the fusion data performance of STARFM in heterogeneous areas by downscaling rather than directly resampling a coarse pixel.The combination of the downscaling mixed pixel method and the spatial and temporal adaptive reflectance fusion model (STARFM) (CDSTARFM) was proposed. In this approach, the downscaling mixed pixel method first reduces the MODIS data that participated in the fusion. Then, the downscaled MODIS data replaces the direct resample MODIS data that appeared in the original STARFM. Finally, the subsequent steps of the STARFM are completed to predict Landsat-like images. The proposed algorithm produced high-resolution temporal synthetic Landsat 8 data based on Landsat 8 and MODIS remote sensing images.Results show that the CDSTARFM was more accurate than STARFM and downscaling methods. The downscaled data used in the CDSTARFM can more fairly reflect surface information than the resampled data of MODIS, which increase the probability of "pure pixels" and allow the CDSTARFM method to more easily determine the "pure" similar pixels in the search window. Therefore, the three indicators[i.e., correlation coefficient (r), root mean square error (RMSE), and the Erreur Relative Globale Adimensionalle de Synthèse (ERGAS)] as well as the scatterplots and visual effect of synthetic images for the CDSTARFM are better than those obtained by STARFM and downscaling methods. Moreover, the optimal window size (11×11 OLI pixels) of the CDSTARFM is smaller than that of the STARFM (31×31 OLI pixels). The accuracy of the CDSATRFM at the same window size is higher than that of the STARFM. The predicted NIR band is evaluated as an example in this study. The correlation coefficients (r) for the CDSTARFM, STARFM, and downscaling methods were 0.96, 0.95, 0.90, respectively; RMSEs were 0.0245, 0.0300, 0.0401, respectively; and ERGAS were 0.5416, 0.6507, 0.8737, respectively. Moreover, synthetic images effectively eliminated the "homogeneous spot" that appeared in the fusion images predicted by the downscaling algorithm and the "boundary of MODIS pixel" that appeared in the synthetic images produced by the STARFM.The temporal-spatial fusion of remote sensing data is an effective approach to solve the conflict of temporal and spatial resolutions of a sensor. This paper used the CDSTARFM algorithm, which combines downscaling method and STARFM to fuse remote sensing data. The CDSTARFM algorithm was verified by the experimental data of Landsat 8 and MODIS images, and it was compared with the STARFM and downscaling methods. The conclusions are listed bellow:(1) CDSTARFM using the downscaled data presented improved results to replace the directly resampled data used in STARFM. The three indicators (i.e., r, RMSE, and ERGAS) and the scatterplots for the CDSTARFM are the best compared with those for STARFM and downscaling methods. (2) The window size at which the best synthetic images were predicted by CDSTARFM is smaller than that by STARFM. In addition, CDSTARFM exhibits better accuracy than STARFM and downscaling methods at the same window size. (3) The synthetic images produced by CDSTARFM are more visually similar to the reference images, especially in the fragmented regions. CDSTARFM can eliminate the homogeneous "plot" in the synthetic images that were generated by downscaling methods and the "MODIS pixel boundary" in the prediction images generated by STARFM in fragmented landscapes.

Key words

decomposition of mixed pixel;downscaling;STARFM;remote sensing;data fusion;CDSTARFM