STARFM算法生成湿地类型TM反射率数据的应用评价

赵艳丽　李大成　贾琇明　崔鹏燕

(太原理工大学测绘科学与技术系　山西 太原 030024)



STARFM算法生成湿地类型TM反射率数据的应用评价

赵艳丽李大成贾琇明*崔鹏燕

(太原理工大学测绘科学与技术系山西 太原 030024)

当前数据获取条件下，很难直接获得兼具高时间与高空间分辨率的多光谱遥感数据,提出利用STARFM(SpatialandTemporalAdaptiveReflectionFusionModel)算法来合成高时间序列的高空间分辨率数据。该算法在我国地理区域的适用性与预测精度验证等工作尚未充分展开。为此，以内蒙古呼伦湖湿地自然保护区为研究样区，并借助于Landsat-5TM(ThematicMapper)与高时序MODIS反射率产品，利用STARFM算法生成具有高时序特征的TM数据，进而将其与真实TM数据进行对比验证分析。结果表明：STARFM算法能够在空间上保持一定预测精度的条件下，对湿地区域内不同地物类别随时相的变化特征具有较好的预测能力，尤其适用于对反射特征随时相变化较小的湿地区域进行时空拟合或数据预测研究。

呼伦湖湿地STARFM高时空分辨率Landsat-5TM

0　引　言

在当前遥感应用条件下，获取和利用高分辨率遥感数据具有重要的应用和发展前景，例如碳循环研究、土地覆盖和土地利用研究[1]，以及大区域制图[2]等均需要几十米级的地表覆盖数据。美国陆地资源卫星Landsat的TM/ETM+(EnhancedThematicMapperPlus)数据(30米)已广泛应用于土地覆盖、土地覆盖变化监测[3]以及生态系统动态变化监测中[4]。但较长的重访周期(16天)以及大气条件和云污染等因素的影响，严重制约了Landsat数据在快速地表变化监测方面的应用[5]。迄今，虽然很多具有高时间分辨率的遥感平台也已投入运行，但其在应用上的不足之处是空间分辨率很低，如MODIS(Moderate-resolutionImagingSpectroradiometer)数据的空间分辨率最高为250米。所以，如何借助这类数据在高时序特征上的优势来弥补LandsatTM/ETM+等中-高空间分辨率数据在捕捉时相特征能力上的不足就成为解决上述问题的有力途径。

时空融合的思想最早由Gao等在2006年提出，他们利用所发展的STARFM融合算法来生成高时间序列的地表反射率数据，并利用模拟数据与真实数据进行了初步验证[5]。该算法在局部尺度下的预测精度高，并得到了广泛的应用推广及验证[6-9]。例如，Schmidt等对利用STARFM算法生成长期密集时间序列数据集进行了分析，取得了很好的预测结果[10]；黄登成等则验证了STARFM算法在提高NPP估算精度中的应用潜力[11]；Hilker等[12]在进行深林干扰的制图研究时，基于STARFM算法提出一种STAARCH(SpatialTemporalAdaptiveAlgorithmformappingReflectanceChange)时空融合算法，使用该算法可以获取地表变化区域和发生变化的具体时相，通过选择最佳时相的数据进行融合来提高融合精度。但该方法在确定变化区域及其对应时相时需要借助多时序的高空间分辨率数据以及相应的低空间分辨率辅助数据，这极大限制了其应用范围。为提高STARFM算法对于地表反射率随时相变化较大的地物(如植被等)的预测精度，ZHU等[13]在STARFM算法的基础上提出了ESTARFM算法，当研究区有两幅以上清晰无云的高空间分辨率数据覆盖时，它具有比STARFM算法更为精确的融合能力[14]。然而，因该算法在推导过程中假设高空间分辨率图像内某一地物类别在两个特定时刻的反射率变化量是一个只与时相差异有关的增量，所以该算法更适用于那些起始融合时相与预测时相均分布在一个季相特征呈单向、线性变化的时段内的情形。

为明确处理二向性反射BRDF(Bidirectionalreflectancedistributionfunction)问题，Roy等[15]在2008年提出了一种基于半物理模型的时空融合算法，该算法在区域尺度具备一定的应用潜力；为解决该算法在计算复杂度、尺度差异等方面的限制，Li等[16]进一步提出了一种扩展的融合策略，该策略对处理大区域数据的时空融合问题具有效率高且光谱连续等优势，并通过引入像元的邻域光谱信息消除了斑块效应的影响[17]。但相比于STARFM算法，上述基于半物理模型的时空融合算法在预测精度以及模型的推广上都有待进一步的应用验证。

本文即以内蒙古呼伦湖湿地自然保护区为例，以Landsat-5TM与MODIS反射率数据为基础研究数据，通过利用STARFM算法生成研究区高时序的TM反射率数据，并分别从空间与时间维度上分析STARFM算法的预测精度及其对于时相变化特征的捕捉能力。

1　研究区概况及数据预处理

1.1研究区概况

呼伦湖是我国第四大淡水湖，以其为核心的湿地自然保护区位于内蒙古呼伦贝尔市辖区内(图1中虚线框)，地理位置为东经116°50′10″～118°10′10″，北纬47°45′50″～49°20′20″，总面积达400 000公顷，平均海拔约545米。保护区内以沙地、盐碱地、湖泊、草地、湿地等分布为主，年平均降雨量约239毫米，主要集中于6至8月份。本文所选研究样区处于呼伦湖东侧的乌尔逊河上，其位置如图1中黑色方框所示。

图1　研究区地理位置

1.2数据预处理

本文所用实验数据为Landsat-5TM影像(空间分辨率30米，回访周期为16天)和MODIS每8天的500米地表反射率产品(MOD09A1)，并搜集了研究区夏季(5月至9月)的TM数据4幅和MODIS数据20幅 (如图2所示，粗线表示起始融合时相下的Landsat-5TM)。其中，TM数据的轨道号为124/26(列/行)，过境日期分别为2009年的6月25日、7月27日、8月12日、9月29日；MODIS数据的轨道号为25/4(列/行)，起终日期为2009年的5月1日和9月30日。在进行融合处理前，需对Landsat-5TM数据和MODIS反射率数据进行相应的数据预处理。

Landsat-5TM数据首先需进行辐射校正，从而将DN值转化为可用于大气校正处理的辐射亮度值。本文利用ENVI的FLAASH大气校正模块对经辐射校正处理后的TM数据进行大气校正，进而将辐射亮度值转换为地表反射率；为与TM数据保持一致，MODIS数据需利用相应的重投影工具(MODISReprojectionTool，MRT)将反射率数据(MOD09A1)的投影类型转换为UTM-WGS84投影坐标系，并采用最邻近法重采样到30米分辨率。经上述处理后的TM数据和MODIS数据还需进行精确的几何配准，从而使它们能够在空间上精确匹配。

图2　Landsat-5 TM和MODIS数据的获取日期

2　方　法

2.1STARFM算法理论基础

当影像间配准误差和大气校正误差影响较小时，某同质像元(xi,yj)在TM影像上t0时刻下的反射率L(xi,yj,t0)可以由其在MODIS影像上同时刻下的反射率M(xi,yj,t0)表示为[5]：

L(xi,yj,t0)=M(xi,yj,t0)+ε0

(1)

其中，εk为由传感器系统误差造成Landsat影像与MODIS影像的光谱差异值。

同理，该像元在时刻tk下的TM反射率数据与相应时刻的MODIS反射率数据间存在以下关系：

L(xi,yj,tk)=M(xi,yj,tk)+εk

(2)

其中，L(xi,yj,tk)和M(xi,yj,tk)分别为(xi,yj)位置上Landsat和MODIS的反射率值。

若Landsat影像与MODIS影像反射率差值不随时间发生变化，则有:

εk=ε0

(3)

因此，该像元在tk时刻下TM影像上的反射率值可表达为：

L(xi,yj,t0)=M(xi,yj,t0)+L(xi,yj,tk)-M(xi,yj,tk)

(4)

然而，由于实际中地表类型的复杂性，它们通常不是同质像元，因而无法满足上述条件下的推导过程。为了解决这一问题，STARFM算法通过借助目标像元周围的邻近像元的光谱信息来构建一个充分考虑了光谱、时间和空间等因素的权重函数Wijk来计算目标像元的反射率值。此权重函数决定了周围邻近像元对目标(中心)像元光谱信息的贡献大小，由光谱差异Sijk、时间差异Tijk以及空间差异Dijk三个因子组成：

(5)

(6)

(7)

其中：A是特定常数。上述三个参量的值越小，则该邻近像元对中心像元的反射率贡献越大，所被赋予的权重系数就越高。权重函数Wijk的计算公式为：

(8)

目标像元的反射率值的最终计算公式如下：

(M(xi,yj,t0)+L(xi,yj,tk)-M(xi,yj,tk))

(9)

其中：W为局部移动窗口的大小，(xw/2,yw/2)为移动窗的中心像元位置，t0为预测时相，L(xw/2,yw/2,t0)为t0时刻中心像元的反射率值，为保证邻近像元光谱信息的有效性，只有在移动窗内且与中心像元光谱相似的清晰像元才参与计算过程。

该融合模型在局部尺度下的预测精度高，本文主要分析和探讨该模型在本研究区对不同随着时相变化的地物特征的预测能力。

2.2时空融合策略与评价方法

本文选取2009年8月12日的Landsat-5TM数据作为STARFM算法在起始融合时相下的数据，并选取与该时相最为接近的2009年8月13日的MODIS数据为起始融合时相下的MODIS数据，然后将2009年5月1日至9月30日间其他时刻的MODIS数据依次作为预测时相下的MODIS融合数据来生成相应时相下的TM反射率数据(共19幅)。

对于本实验所需生成的19个时相的预测结果中，6月26日、7月28日以及9月30日的TM反射率数据预测结果分别与6月25日、7月27日以及9月29日下真实的TM反射率数据进行融合精度的对比分析，由于这两组数据在时相上仅相差一天，因而可以认为两组数据间的反射率差异能够反映出预测数据与真实数据的反射率差异。本文首先利用基于反射率的散点分布图来分析上述两组数据分别在三个时相下反射率数据的相关性，然后进一步对两组数据中不同地物类别的反射率数据的相关性以及统计特征进行对比分析；最后分析短时期内(一个季相)不同地物类别在所预测的多时序TM数据中的NDVI变化特征与MODIS中NDVI数据的一致程度，从而验证不同地物随时相的变化特征。

考虑到研究区的分类精度会直接影响到本文后续的实验分析及评价过程，因而需选取分类精度高，并且处理性能稳定的监督分类方法来进行分类处理。在充分了解研究区地物分布特征的基础上，本文采用最大似然分类方法来对研究区在起始融合时相(2009年8月12日)下的Landsat-5TM数据进行分类处理。这里将研究区域按照草地、湿地、盐碱地、水体、沙地等共5种地物类别进行划分，其中每种类别至少选取10个以上的训练样本作为分类依据，样本可分离度均大于1.9。将分类结果中的小图斑剔除后，用分类结果叠加方法目视判断分类结果的准确性，将部分不准确的分类结果合并到相应类中。

对于本文所搜集的4个时相的TM数据， 2009/07/27、2009/08/12以及2009/09/29三个时相下的影像上有水体的存在，而2009年6月25日的TM影像上却不存在水体类别。因此，虽然上述4幅TM数据基本分布于一个季相范围内，但相对于其他三个时相，研究区在6月25日下的地物反射特征存在显著变化。

3　实验分析与讨论

本文拟分别从空间预测精度和时序数据拟合程度上分析STARFM算法的预测精度及其对不同地物类别随时相变化特征的捕捉能力。评价参数采用相关系数r2来分析实验结果的预测一致性。相关系数([0,1])表征着预测结果与真实值的接近程度，其值越大，相关程度越高，数据间的一致性就越好。

图3展示了三个预测时相下的融合结果与真实TM反射率数据在近红、红和绿三个波段上的散点分布图。其中，横坐标为真实的TM反射率数据，纵坐标为STARFM算法的预测结果；虚线为理想拟合线，实线为预测结果与真实数据的拟合线(其直线表达式已在图中标出)，预测结果与真实反射率数据间的相关系数r2标于图中左上方。由图可知，三个时相在近红外、红、绿这三个波段的相关系数均在0.7以上(r2>0.73)，尤其是在红波段和绿波段上的相关系数均在0.8以上，相关程度很高。若对比每个波段在不同时相下的相关性，可发现：9月30日、6月26日以及7月28日下的预测结果与相应真实数据间的相关程度依次增高。同时，对于上述任一时相而言，红波段和绿波段上的相关程度均比近红外波段上的相关程度高。由于7月28日的预测结果在时相上与起始融合时相(8月12日)相差较小，地表覆盖变化很小，因而其相关程度也较其他两个时相高。此外，研究区域在6月和7月间正处于植被生长旺盛期，而9月30日下的植被生长则趋于衰落阶段，因而6月26日的相关性较9月30日的相关性要好。

图3　真实Landsat-5 TM影像与其预测结果的相关性

表1展示了不同地物类别的真实反射率与其预测结果间的相关性统计信息。可以发现，在研究区的不同地物类别中，盐碱地、湿地以及沙地的真实反射率与其预测结果的相关程度较高，而每一种地类在三个目标时相(9月30日、6月26日以及7月28日)下的相关程度依次增高，且在红波段和绿波段上的相关系数整体上要高于近红外波段。

表1　湿地区域不同地物类别预测结果的相关性统计

此外，在2009年6月26日的预测结果中，水体的预测精度最差(相关性最低)。参照6月25日的TM真实数据可知，该时相下研究区内的水体已接近干涸，地物类别发生了明显转变，而其它两个时相(9月30日和7月28日)下的水体分布特征则与起始融合时相一致。不同地物类别在不同预测时相下的融合结果与相应的TM真实图像的差异数据的统计信息(均值与标准差)如图4所示。横坐标表示不同地物类别，纵坐标为预测结果影像与真实反射率数据的差异大小。十字符、米字符、三角符分别代表近红外、红、绿三个波段，实线和虚线分别表示均值和标准差的大小。可见，每种地物类别在所有时相和波段上的差异大小均在零值附近波动，其中沙地的统计特征最接近零值，而水体的拟合结果最差。总体上看，所有类别的差异程度从9月30日、6月26日至7月28日基本上依次降低，即融合精度依次增大，即与图5的结论一致。

图5为不同地物类别(盐碱地-saline，湿地-wetland，沙地-sandyland，草地-grass，水体-water)在预测的时序数据上及其在真实MODIS数据上的NDVI曲线。易发现，所预测的研究区各地类的NDVI变化曲线与MODIS数据的NDVI参考曲线基本一致。具体而言，5月到6月间，植被生长相对缓慢，两者的NDVI曲线较为接近；7月至8月间，植被处于快速生长时期，因而两者的NDVI数据的差异逐渐增大，且对于草地的而言，这种差异最为突出；而在九月底时，植被生长季结束，两者在NDVI数据上的差异程度也相应降低。

总体而言，STARFM算法能够在保持空间预测精度的基础上，对不同地物类别随时相的变化特征有较好的拟合能力；但对于那些反射特征随时相变化较大的地物，其预测精度稍差。

图4　STARFM对湿地区域不同地物在不同时相下的反射率预测精度

图5　STARFM对湿地区域不同地物在一个季相下的NDVI预测曲线

4　结　语

本文采用STARFM算法来预测湿地研究区下多时序的TM反射率数据，并从单时相数据的空间预测精度与多时相数据的时序变化特征上对预测结果进行充分的应用分析与评价。依据本文的实验及分析结果，可初步得到以下结论：

(1) 在一个季相分布内，起始融合时相与预测时相的差异越小，STARFM算法对于TM反射率数据的预测精度就越高(如实验中2009年7月28日的预测结果)。而研究区在两个时相间的地物类别发生改变会显著影响该算法的融合效果，如本实验中因季相特征变化引起的部分湖泊干涸等。

(2)STARFM算法对于湿地区域不同地物类别的预测精度也不相同，对于反射特征随时相变化较小的地物，其预测精度较高；反之，其预测效果稍差。因此，该算法适用于地物反射特征随时相变化较小的地理区域进行时空融合研究。

(3)STARFM算法对一个季相条件下湿地区域不同地物类别随时相的变化或分布特征具有很好的拟合能力，其本质上是能够很好地保持融合前后各波段反射率数据间的相对关系，因而该算法在模拟或预测地物反射的变化特征或规律上具有一定的应用潜力。

[1]IGBP,HDP.Land-useandland-coverchange：science/researchplan[R].IGBPReportNo.35andHDPReportNo.7.Stockholm, 1995.

[2]CihlarJ.Landcovermappingoflargeareasfromsatellites:statusandresearchpriorities[J].InternationalJournalofRemoteSensing, 2000,21(6-7):1093-1114.

[3]WulderMA,WhiteJC,GowardSN,etal.Landsatcontinuity:Issuesandopportunitiesforlandcovermonitoring[J].RemoteSensingofEnvironment, 2008,112(3):955-969.

[4]HealeySP,CohenWB,YangZQ,etal.ComparisonofTasseledCap-basedLandsatdatastructuresforuseinforestdisturbancedetection[J].RemoteSensingofEnvironment, 2005,97(3):301-310.

[5]GaoF,MasekJ,SchwallerM,etal.OntheblendingoftheLandsatandMODISsurfacereflectance:PredictingdailyLandsatsurfacereflectance[J].GeoscienceandRemoteSensing,IEEETransactionson, 2006,44(8):2207-2218.

[6]HilkerT,WulderMA,CoopsNC,etal.GenerationofdensetimeseriessyntheticLandsatdatathroughdatablendingwithMODISusingaspatialandtemporaladaptivereflectancefusionmodel[J].RemoteSensingofEnvironment, 2009,113(9):1988-1999.

[7]LiuF,WangZ.SyntheticLandsatdatathroughdataassimilationforwinterwheatyieldestimation[C]//Geoinformatics, 2010 18thInternationalConferenceon.IEEE, 2010:1-6.

[8]WattsJD,PowellSL,LawrenceRL,etal.Improvedclassificationofconservationtillageadoptionusinghightemporalandsyntheticsatelliteimagery[J].RemoteSensingofEnvironment,2011,115(1):66-75.

[9]SinghD.Evaluationoflong-termNDVItimeseriesderivedfromLandsatdatathroughblendingwithMODISdata[J].Atmósfera,2012,25(1):43-63.

[10]SchmidtM,UdelhovenT,GillT,etal.LongtermdatafusionforadensetimeseriesanalysiswithMODISandLandsatimageryinanAustralianSavanna[J].JournalofAppliedRemoteSensing,2012,6(1):063512-1-063512-18.

[11] 黄登成,张丽,尹晓利,等.数据融合技术在提高NPP估算精度中的应用[J/OL].计算机工程与应用,2013.http://www.cnki.net/kcms/detail/11.2127.TP.20130418.1618.016.html.

[12]HilkerT,WulderMA,CoopsNC,etal.Anewdatafusionmodelforhighspatial-andtemporal-resolutionmappingofforestdisturbancebasedonLandsatandMODIS[J].RemoteSensingofEnvironment, 2009,113(8):1613-1627.

[13]ZhuX,ChenJ,GaoF,etal.Anenhancedspatialandtemporaladaptivereflectancefusionmodelforcomplexheterogeneousregions[J].RemoteSensingofEnvironment, 2010,114(11):2610-2623.

[14] 朱孝林.高时空分辨率遥感数据生成技术研究[D].北京师范大学，2010.

[15]RoyDP,JuJ,LewisP,etal.Multi-temporalMODIS-Landsatdatafusionforrelativeradiometricnormalization,gapfilling,andpredictionofLandsatdata[J].RemoteSensingofEnvironment, 2008,112(6):3112-3130.

[16]LiDC,TangP,HuCM,etal.Spatialtemporalfusionalgorithmbasedonanextendedsemi-physicalmodelanditspreliminaryapplication[J].JournalofRemoteSensing, 2014,18(2):307-319.

[17] 李大成.时空融合方法生成高时空分辨率区域地表反射率数据集的技术研究[D].中国科学院大学，2013.

APPLICATIONEVALUATIONOFSTARFMALGORITHMINGENERATINGWETLAND-TYPETMREFLECTANCEDATA

ZhaoYanliLiDachengJiaXiuming*CuiPengyan

(Department of Mapping Science and Technology, Taiyuan University of Technology, Taiyuan 030024, Shanxi, China)

Multispectralremotesensingdatawithbothhightemporalandhighspatialresolutionscanhardlybeobtaineddirectlyundercurrentdataacquisitionconditions.WeproposedtousSTARFM(spatialandtemporaladaptivereflectancefusionmodel)tosynthesisethehightimeseriesdatawithhighspatialresolution.Nevertheless,theapplicabilityofthismodelforthegeographicareasinChinaanditspredictionprecisionvalidationhavenotyetbeensubstantiallydone.Accordingly,theHulunWetlandNatureReserveinInnerMongoliaisselectedasthestudyarea,andwiththehelpofLandsat-5TMdataandhightimeseriesMODISreflectanceproducts,weusedSTARFMalgorithmtogeneratetheTMdatawithhightimeseriesfeature,andthenmadethecomparativevalidationanalysisonitwiththeactualTMreflectancedata.ResultshowedthattheSTARFMalgorithmhasahigherpredictioncapabilityonthephasefeaturesofdifferentterraincategoriesinwetlandareavaryingalongwithtimephasesundertheconditionofpreservingcertainpredictionaccuracyinspace,itisparticularlysuitableforthespatiotemporalfittingordatapredictionresearchinthosewetlandareaswherethereflectioncharacteristicschangelittlealongwiththetimephases.

HulunwetlandSTARFMHighspatiotemporalresolutionLandsat-5TM

2014-07-24。国家高技术研究发展计划项目(2009AA122002)。赵艳丽，硕士，主研领域：遥感信息处理以及遥感与地理信息系统的一体化集成。李大成，讲师。贾琇明，教授。崔鹏燕，硕士。

TP751

ADOI:10.3969/j.issn.1000-386x.2016.03.063