CN111242022B - 基于低分辨率遥感产品降尺度的高分辨率fapar估算方法 - Google Patents

Abstract

基于低分辨率遥感产品降尺度的高分辨率FAPAR估算方法，先基于先验知识库，通过地表单元划分对MODIS FAPAR反演模型进行分层线性化，建立MODIS FAPAR与Landsat反射率间的尺度转换函数(即尺度自洽的FAPAR估算模型)；基于贝叶斯框架融合MODIS FAPAR和Landsat反射率新观测数据，更新MODIS FAPAR与Landsat反射率间的尺度转换函数；最后基于新观测Landsat数据实时估算FAPAR。解决了现有估算方法难以同时实现尺度自洽、空间连续、实时的高分辨率FAPAR估算的问题。

Description

基于低分辨率遥感产品降尺度的高分辨率FAPAR估算方法

技术领域

本发明属于遥感估算领域，涉及一种基于低分辨率遥感产品降尺度的高分辨率FAPAR估算方法。

背景技术

光合有效辐射吸收比率(Fraction of Absorbed Photosynthetically ActiveRadiation,FAPAR/FPAR)定义为植被冠层(绿色部分)吸收的光合有效辐射(Photosynthetically Active Radiation,PAR)占总入射PAR的比率，是表征植被光合作用的核心参数。一方面，FAPAR直接反映大气圈与生物圈间的物质交换过程，是建立碳循环模型、陆地生态系统模型的关键变量，也是植被生长状态与土地覆盖监测中的重要指标。另一方面，FAPAR反映植被冠层对光能的截获利用情况，是地气系统辐射平衡与能量交换的重要影响因素。因此，FAPAR的准确估算，是定量估算光合作用、提高全球生态系统碳估算精度的基础，在全球气候变化、碳排放、粮食安全等重大问题研究方面有着重要的应用价值。

尽管遥感估算已成为获取面状FAPAR信息的唯一有效手段，但目前全球及区域尺度上通用的FAPAR产品间精度差别较大，且都是在低分辨率尺度上生产的，丢失了地表大量的空间细节信息。公里级产品亚像元尺度上的地表异质性，不仅会影响FAPAR的反演精度，造成多尺度FAPAR产品间难以自洽，还直接影响生态系统生产力和植被生物量的估算。获取高空间分辨率FAPAR遥感产品，对于改进FAPAR遥感产品的精度和空间连续性，更好地满足生态学应用需求具有重要意义。

目前，获取高空间分辨率FAPAR有两种方式：一是采用基于植被指数的半经验方法或基于机理模型的物理反演方法，直接由高分辨率遥感数据估算FAPAR；二是以成熟的业务化FAPAR遥感产品为参考，通过降尺度转换得到高分辨率尺度上的FAPAR。然而，无论是采用基于植被指数的半经验方法或基于机理模型的物理反演方法，在较大范围上应用时都会面临算法普适性、参数可获取性、多尺度自洽性和真实性检验等多重挑战。相比直接反演方法，以MODIS(Moderate-resolution Imaging Spectroradiometer)为代表的业务化遥感产品建立在严格的物理模型基础上，在全球范围积累了大量的历史观测数据，并在不断产生新观测数据，算法及产品经过了数十年来的广泛检验，精度较为可靠。以成熟的业务化FAPAR遥感产品为参考进行降尺度转换，成为获取高分辨率FAPAR的一种可行方式。

T.Hwang,C.Song,P.V.Bolstad,and L.E.Band,"Downscaling real-timevegetation dynamics by fusing multi-temporal MODIS and Landsat NDVI intopographically complex terrain,"Remote Sensing of Environment,vol.115,pp.2499-2512,2011提出了一种MODIS FAPAR降尺度转换的方法，该方法假设MODIS FAPAR与Landsat NDVI(归一化植被指数，Normalized Difference Vegetation Index)间的关系具有尺不变性，逐像元计算MODIS FAPAR与Landsat NDVI聚合值的转换系数，再将其用于从Landsat NDVI中估算FAPAR(即基于NDVI的线性转换系数法)。由于转换系数是在MODIS像元尺度上逐日计算的，该方法考虑到了降尺度转换系数的季候差异，但却忽略了MODIS亚像元尺度上的空间异质性。另外，NDVI的尺度效应已被大量研究证实，FAPAR与NDVI间转换关系的尺不变假设显然会在异质性地表引入尺度误差，导致结果难以自洽。

为在降尺度转换过程中尽量避免尺度效应，F.Gao and W.P.Kustas,"Simplemethod for retrieving leaf area index from Landsat using MODIS leaf areaindex products as reference,"Journal of Applied Remote Sensing,vol.6,p.063554,2012提出了一种基于多尺度回归树模型的MODIS LAI(叶面积指数，Leaf AreaIndex)降尺度方法，也即一些现有专利提到的cubist方法。该方法利用多期Landsat反射率数据和MODIS LAI产品，选择纯像元样本在MODIS尺度上建立反射率与LAI的多尺度线性回归树模型，并将该模型用于Landsat尺度上LAI的估算。采用线性模型可以避免在降尺度转换中产生尺度效应。但该方法对全年所有日期数据采用相同的转换关系，忽略了LAI与反射率间的关系存在季候差异。这对于时间变化特征较强的植被参数(如FAPAR)，显然会引入额外误差。尤其是当需要对MODIS产品进行长时间序列上、实时的降尺度转换时，该方法难以对新观测数据进行融合更新，而是每次都需选择当年的多期数据来建立转换关系，难以达成实时的降尺度转换。

总结现有的两种典型算法，基于NDVI的线性转换系数法能够实现实时的降尺度转换，考虑到了植被参数的季候差异，但却是以牺牲MODIS亚像元尺度上的空间细节信息、忽略尺度效应为代价的。基于多尺度回归树模型的MODIS LAI降尺度方法能够保证结果的空间连续性和自洽性，但却忽略了植被参数的时变特征，难以达成实时的降尺度转换。从先验知识利用的角度分析，基于NDVI的线性转换系数法只利用了时间上的先验知识，而基于多尺度回归树模型的MODIS LAI降尺度方法仅利用空间上的先验知识。此外，基于NDVI的线性转换系数法是将新观测作为数据进行一次性的利用，基于多尺度回归树模型的MODIS LAI降尺度方法是建立在先验知识之上的，但却难以对先验知识进行更新。因此，需要在多尺度自洽的前提下，发展更为普适的、可综合利用各类先验知识并可有效更新的遥感数据降尺度转换方法。

发明内容

本发明实施例的目的在于提供一种基于低分辨率遥感产品降尺度的高分辨率FAPAR估算方法，以解决现有估算方法难以同时实现尺度自洽、空间连续、实时的高分辨率FAPAR估算的问题。

本发明所采用的技术方案是，基于低分辨率遥感产品降尺度的高分辨率FAPAR估算方法，按照以下步骤进行：

步骤S1，基于FAPAR遥感估算机理、遥感产品尺度转换理论、历史年限MODIS FAPAR产品与Landsat地表反射率数据、相关地理要素数据(土壤类型与地表覆盖类型)等先验知识，通过地表单元划分，对MODIS FAPAR反演模型进行分层线性化，建立MODIS FAPAR与Landsat地表反射率间的尺度转换函数(即尺度自洽的FAPAR先验估算模型)；

步骤S2，基于贝叶斯框架融合MODIS FAPAR和Landsat地表反射率新观测数据，更新MODIS FAPAR与Landsat地表反射率间的尺度转换函数，得到更新的FAPAR后验估算模型；

步骤S3，基于更新后的FAPAR后验估算模型，基于新观测Landsat数据实时估算FAPAR。

本发明实施例的有益效果是：通过地表覆盖类型和土壤类型对地表单元进行合理划分，针对每一地表单元建立FAPAR与多波段反射率间的线性估算模型，实现了对低分辨率FAPAR业务化遥感产品生产所用的复杂辐射传输模型的线性化，以避免尺度效应，实现了自洽的尺度转换；对每一地表单元，通过选取纯像元样本，建立MODIS FAPAR与Landsat地表反射率聚合值之间的尺度转换函数，采用统计分析方法求解模型参数，建立了尺度自洽的FAPAR估算模型，再将该模型用于高分辨率(即Landsat尺度)像元尺度，实现了空间连续的高分辨率FAPAR估算；同时，在遥感趋势面的框架下，利用多期历史数据构建FAPAR先验估算模型，再由贝叶斯定理融合新观测信息(即每一景实时获取的MODIS FAPAR和Landsat遥感影像)，对先验估算模型参数进行更新，得到融合了先验知识和新观测信息的FAPAR估算模型，其同时利用了时间、空间上的信息，是关于高分辨率FAPAR全部的知识，实现了实时的高分辨率FAPAR估算。通过上述方法同时实现了尺度自洽、空间连续、实时的高分辨率FAPAR估算，有效解决了现有估算方法难以同时实现尺度自洽、空间连续、实时的高分辨率FAPAR估算的问题。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例的MODIS FAPAR遥感数据的空间降尺度转换方法流程图。

图2(a)是本发明实施例采用不同回归模型进行回归分析R²指标拟合精度图。

图2(b)是本发明实施例采用不同回归模型进行回归分析RMSE指标拟合精度图。

图3是采用三种降尺度转换算法得到的2012年年积日为112的Landsat FAPAR结果与MODIS FAPAR对比图。

图4是采用三种降尺度转换算法得到的2012年年积日为176的Landsat FAPAR结果对比图。

图5是采用三种降尺度转换算法得到的2012年年积日为224的Landsat FAPAR结果对比图。

图6是采用三种降尺度转换算法得到的2012年年积日为272的Landsat FAPAR结果对比图。

图7是采用三种降尺度转换算法得到的2012年年积日为304的Landsat FAPAR结果对比图。

图3-7中，a为MODIS FAPAR，b为多尺度回归树算法结果，c为基于NDVI的线性转换系数法结果，d为本发明实施例的方法所得结果。

图8(a)是采用多尺度回归树方法所得结果与MODIS FAPAR的散点图。

图8(b)是采用基于NDVI的线性转换系数法所得结果与MODIS FAPAR的散点图。

图8(c)是采用本发明实施例的方法所得结果与MODIS FAPAR的散点图。

图9(a)是采用多尺度回归树方法所得结果与地面测量数据的散点图。

图9(b)是采用基于NDVI的线性转换系数法所得结果与地面测量数据的散点图。

图9(c)是采用本发明实施例的方法所得结果与地面测量数据的散点图。

图10(a)是2012年6月24日ASTER FAPAR真实参考数据的空间分布图。

图10(b)是2012年7月10日ASTER FAPAR真实参考数据的空间分布图。

图11(a)是采用多尺度回归树方法所得2012年6月24日FAPAR空间分布图。

图11(b)是采用多尺度回归树方法所得2012年7月10日FAPAR空间分布图。

图12(a)是采用基于NDVI的线性转换系数法所得2012年6月24日FAPAR空间分布图。

图12(b)是采用基于NDVI的线性转换系数法所得2012年7月10日FAPAR空间分布图。

图13(a)是采用本发明实施例的方法所得2012年6月24日FAPAR空间分布图。

图13(b)是采用本发明实施例的方法所得2012年7月10日FAPAR空间分布图。

图14(a)是采用多尺度回归树方法所得结果与ASTER FAPAR的散点图。

图14(b)是采用基于NDVI的线性转换系数法所得结果与ASTER FAPAR的散点图。

图14(c)是采用本发明实施例的方法所得结果与ASTER FAPAR的散点图。

图15(a)是农田典型像元时间序列曲线图。

图15(b)是林地典型像元时间序列曲线图。

图15(c)是稀疏植被典型像元时间序列曲线图。

图15(d)是裸地典型像元时间序列曲线图。

图16是各地面站点实测FAPAR与发明实施例的方法的反演FAPAR的时间趋势对比图。

图17(a)是研究区Landsat影像数据自动分类图。

图17(b)是研究区土地覆盖类型图。

图17(c)是研究区地表单元图，其中图例代表地表单元类型代码，1-90对应特定土壤类型与地表覆盖类型组合的地表单元。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例的基于现有业务化低分辨率FAPAR遥感产品降尺度的高空间分辨率FAPAR估算方法，可适用于各种高空间分辨率遥感数据的FAPAR估算，如对Landsat5/Thematic Mapper(TM),Landsat 7Enhanced Thematic Mapper Plus(ETM+),Landsat 8/Operational Land Imager(OLI),Advanced Spaceborne Thermal Emission andReflection Radiometer(ASTER),HJ-1/CDD,Sentinel-2MSI等空间分辨率不低于30m的遥感数据均具有普遍的适用性。

本实施例以MODIS FAPAR和Landsat遥感数据为例。但把Landsat换成以上这些高空间分辨率数据同样适用。

本发明实施例提出一种基于MODIS FAPAR降尺度的Landsat FAPAR估算方法，以对2012年的Landsat FAPAR进行估算为例进行说明，具体流程如图1所示，按照以下步骤进行：

步骤S1、基于先验知识库(FAPAR遥感估算机理、遥感数据尺度转换理论、2000～2011年的MODIS FAPAR历史数据和Landsat历史影像、地表覆盖类型数据、土壤类型数据)，通过地表单元划分对MODIS FAPAR反演模型进行分层线性化，建立MODIS FAPAR与Landsat地表反射率间的尺度转换函数(即尺度自洽的FAPAR先验估算模型)，具体步骤如下：

步骤S11、数据收集与预处理。

本发明实施例使用的数据主要为Landsat 5TM和Landsat 7ETM+地表反射率数据。Landsat TM/ETM+数据具有多个波段，影像空间分辨率30m。Landsat卫星重访周期为16天。受天气状况、质量因素等影响，实际可用的Landsat数据时间分辨率远低于16天。本发明实施例共筛选了2000年2月至2011年12月间云量小于10％的Landsat TM\ETM+地表反射率数据共103景，2012年高质量Landsat TM\ETM+地表反射率数据14景。其中，历史数据(2000-2011年)平均时间分辨率约为9景/年，用于从历史观测中归纳、总结提取先验知识；2012年的数据作为新观测信息。下载的Landsat TM\ETM地表反射率数据产品空间基准为WGS-84坐标系和UTM投影，对数据进行解压缩、裁剪和数值转换等预处理。

本发明实施例使用的MODIS FAPAR产品为最新发布的C6(Collection 6)版本。相较于前一版本，C6版本的主要改进在于空间分辨率从1km提高到500m。本发明实施例收集了2000年2月至2012年12月的583景MOD15A2H C6产品数据，空间分辨率约500m，所使用的参数主要包括FAPAR及相关质量标识(LAI\FAPAR QC,FAPAR STD)等。参照MODIS LAI\FAPAR产品手册，对MOD15A2H进行数据层提取、裁剪，投影转换(至WGS-84坐标系UTM投影)，重采样(至480m)、数值转换等处理。

由于Landsat影像和MODIS FAPAR产品具有不同的时间分辨率，首先筛选出与Landsat影像观测日期对应的MODIS FAPAR产品。

步骤S12、按照植被生长旺盛期、无云或少云、成像质量较好等标准，选择研究区具有代表性的一景影像(本实施例选用2012年7月10日的Landsat影像)，采用K均值聚类算法将研究区的Landsat影像自动划分为林地、农田、稀疏植被、裸地和其他(如云)等五种地类类别，将Landsat分类结果按照众数规则升尺度至1km尺度，即得到了分辨率1km的地表覆盖类型数据。基于世界土壤数据库(Harmonized World Soil Database,HWSD)的中国土壤数据集(v1.1)(空间分辨率1km)，使用土壤类型栅格数据，对土壤类型数据和地表覆盖类型数据进行空间上的交集操作，将研究区划分为相应的地表单元，每一地表单元具有特定的土壤类型和地表覆盖类型。即对土壤类型数据(共18个类别)和地表覆盖类型(共5种)进行空间上的交集操作，理论上可将研究区划分为18×5种地表单元。由于某些土壤类型单元并不覆盖所有地表覆盖类型，研究区实际上共划分为73个地表单元，如图17(a)-(c)所示。

土地覆盖类型代表地表覆盖的一般状态。在地表覆盖类型未发生剧烈变化(如人为干扰、自然灾害等)的情况下，只需要对一景代表性的Landsat影像进行地表覆盖分类即可。如果地表覆盖类型发生了剧烈变化，可以有针对性地、在发生变化的日期后选择一景Landsat影像分析、重新划分地表单元；在这种情况下，按照发生变化前和发生变化后两个时期，可有两个不同的地表单元数据，在数据筛选与分析中，按照成像时间(即是否已经发生了变化)将Landsat数据和MODIS数据对应至相应的地表单元数据。

步骤S13、筛选MODIS FAPAR高质量纯像元样本。

(1)，筛选MODIS FAPAR高质量反演值：MODIS FAPAR产品生产基于主算法和备用算法两种算法。其中，主算法基于复杂三维辐射传输模型的查找表算法，缺少显性表达式，但具有严格的物理基础，精度较高。本发明实施例的目的也即：将具有高精度的MODIS FAPAR主算法进行分层线性化，得到其显性表达式。因此，只选用MODIS主算法反演的FAPAR高质量反演数据建立样本。筛选标准为MODIS FAPAR质量控制标识QC＝0。

(2)，MODIS纯像元筛选：MODIS尺度上的纯像元是建立MODIS FAPAR与Landsat地表反射率间尺度转换函数的关键。本发明实施例采用升尺度地类频率和变异系数CV(Coefficient of Variation)两种指标来度量MODIS亚像元尺度上地表异质性，进行纯像元筛选。变异系数CV定义为标准差与均值的比值，CV值表示地表反射率的均匀程度，值越小，像元越均质。升尺度地类频率代表优势地类在亚像元尺度上所占的面积比例，值越大，像元越均质。

为保证筛选出的像元尽可能均质，且能够包含足够多的样本，本实施例按照CV≤0.2且升尺度地类频率≥0.95的标准筛选纯像元。同时，为避免在非植被生长季选入过多的裸地像元作为样本，影响样本的代表性，在纯像元筛选中仅使用植被生长期(4至10月)的数据。

具体地，利用Landsat数据计算纯像元筛选的指标。对每景Landsat地表反射率数据进行如下处理：

a)，将Landsat地表反射率升尺度至MODIS尺度(480m分辨率)，并计算变异系数CV(Coefficient of Variation)。变异系数定义为标准差与均值的比值。具体计算过程为：首先，对Landsat可见光-近红外波段的每一个波段，计算每一MODIS像元的波段CV值，即该MODIS像元空间范围包含的16*16个Landsat像元波段反射率的标准差与均值的比值；其次，计算每一个MODIS像元各个波段(即Landsat可见光-近红外的6个波段)CV值的均值，作为该像元的CV值。

b)，将Landsat地表分类数据升尺度至480m，按照众数规则得出MODIS像元尺度上的优势地类，即升尺度像元对应的地类。同时，计算该优势地类在MODIS亚像元尺度上所占的面积比例，即升尺度地类频率。

(3)，对2000-2011年间的每一景Landsat地表反射率数据及MODIS FAPAR数据，按照步骤(1)-(2)中的过程，筛选出具有MODIS FAPAR高质量反演值的纯像元。然后，根据S12中划分的地表单元，筛选出每一类地表单元对应的MODIS FAPAR高质量反演值的纯像元；提取所有筛选像元对应的MODIS FAPAR值和聚合至480m分辨率的Landsat可见光和近红外波段的地表反射率值，即构成每一地表单元的筛选样本。

步骤S14、建立MODIS FAPAR与Landsat地表反射率间自洽的尺度转换函数，对每一地表单元，通过对筛选样本进行回归分析得出FAPAR先验估算模型的参数。

具体过程如下：

MODIS FAPAR产品是基于冠层三维辐射传输模型，采用查找表反演策略生产的。MODIS FAPAR值与地表反射率间的关系可表示为：

其中，下标M代表MODIS尺度；FAPAR_M代表MODIS FAPAR值；S_M代表MODIS尺度的地表反射率；f_M表示MODIS尺度的FAPAR估算模型，既可以代表简单的经验线性模型，也可以代表复杂的物理模型，对于MODIS产品，f_M代表一个复杂的物理模型；下标L代表Landsat尺度，i表示给定MODIS像元内对应Landsat像元的索引，S_L(i)表示MODIS像元内第i个Landsat像元的地表反射率；n代表MODIS像元对应Landsat像元的总数。本发明实施例使用的MODISFAPAR数据空间分辨率480m，一个MODIS像元内部对应16*16个Landsat像元。

MODIS FAPAR空间降尺度的问题可表述为：如何在已知FAPAR_M和S_L(i)的情况下，求出f_M，且f_M能够适用于Landsat尺度。即f_M＝f_L，f_L为Landsat尺度FAPAR估算模型。由尺度效应产生的机理分析[1,2]，若f_M为非线性模型，则将f_M不加修正地用于Landsat尺度时会产生尺度效应，则f_M≠f_L，无法实现自洽的尺度转换。解决这个问题的困难还在于，对MODIS FAPAR产品，f_M代表一个复杂的物理模型，不具备显性的数学表达式。因此，需要由FAPAR的遥感反演机理入手，将MODIS FAPAR产品所代表的物理反演模型进行线性化，表示为地表反射率的线性组合。

文献[1]H.Wu and Z.Li,"Scale Issues in Remote Sensing:A Review onAnalysis,Processing and Modeling,"Sensors,vol.9,pp.1768-1793,2009.

文献[2]Z.Hu and S.Islam,"A framework for analyzing and designingscale invariant remote sensing algorithms,"IEEE Transaction on Geoscience andRemote Sensing,vol.13,pp.747-755,1997.

MODIS FAPAR产品算法考虑了植被结构、地表类型和土壤背景等多种因素间的相互作用，且现有研究表明在建立线性回归模型时考虑地表类型可以减小地表参数估算中的不确定性。因此，FAPAR与地表反射率间的关系与地表类型、土壤背景相关。本实施例结合土壤类型和地表覆盖类型进行地表单元划分，对每一地表单元建立特定的FAPAR估算模型，形式为：

其中，下标soil代表土壤类型，下标lc代表地表覆盖类型。

表示特定土壤类型和地表覆盖类型条件下MODIS尺度FAPAR估算模型。

由上述分析可知，当

为线性模型时，有

表示特定土壤类型和地表覆盖类型条件下Landsat尺度FAPAR估算模型。建立在MODIS尺度上的FAPAR估算模型可适用于Landsat尺度FAPAR的估算，即建立了尺度自洽的FAPAR估算模型，则Landsat尺度FAPAR的估算FAPAR_L可由下式计算：

其中，s_j表示第j波段的地表反射率，

表示线性函数

第j波段的系数，k为波段总数，

为线性函数

的常数项。

Landsat尺度上，对同类地表单元，FAPAR估算采用相同的线性模型

则MODIS纯像元上的FAPAR可由下式计算：

式中，s_i,j为MODIS像元内对应第i个Landsat像元第j波段的地表反射率；

为该MODIS像元内所有Landsat像元第j波段地表反射率的均值，也就是Landsat像元第j波段地表反射率聚合至MODIS 480m分辨率的均值。

因此，对每一地表单元，通过筛选具备MODIS FAPAR高质量反演值的MODIS纯像元，即可在已知MODIS FAPAR和Landsat多波段的地表反射率聚合值(即聚合至MODIS 480m分辨率的地表反射率均值)的情况下，通过回归分析，求解该地表单元对应的FAPAR线性估算模型的系数

和

得出该地表单元对应的FAPAR先验估算模型。

根据采用的波段组合不同，降尺度转换函数可有多种形式，本实施例选用三种模型形式进行对比分析，依次为：基于Landsat 2-7波段反射率的多元回归模型、基于Landsat3-4波段反射率的二元回归模型，以及基于NDVI的一元回归模型。分别采用三种模型对提取的MODIS FAPAR纯像元高质量反演样本进行回归分析，以决定系数R²和均方根误差RMSE(Root Mean Squared Error)为指标对比模型拟合精度，如图2(a)～(b)所示，可知，基于Landsat 2-7波段反射率的多元回归模型拟合R²最高，RMSE最低，精度最高。基于Landsat3、4波段反射率的二元回归模型拟合精度次之，且明显优于基于Landsat NDVI的一元线性回归模型。相较而言，基于Landsat 3、4波段反射率的二元回归模型在RMSE上与基于Landsat 2-7波段反射率的多元回归模型差异很小，二者间最大差异为0.0095。Landsat 3-4波段对应红光、近红外波段，与植被特征关系最为密切，且基于Landsat 3、4波段反射率的二元回归模型模拟精度高于基于Landsat NDVI的一元线性回归模型，模型复杂度低于基于Landsat 2-7波段反射率的多元回归模型，因此综合考虑模型模拟精度、模型复杂度和模型尺度效应等因素，本实施例采用基于Landsat 3、4波段反射率的二元回归模型作为FAPAR估算模型的基本形式，表示为：

式中，S代表地表反射率；s₃和s₄分别为Landsat第3(红光)、4(近红外)波段的地表反射率；

和

分别表示线性函数f_soil,lc的常数项、变量s₃和s₄的系数，即模型系数。由于该模型适用于MODIS和Landsat两个尺度，这里不再用下标M或L区分，代表FAPAR估算的通用模型。f_soil,lc表示尺度自洽的FAPAR估算模型，即同时适用于MODIS和Landsat两个尺度的FPAR估算模型。对应地，FAPAR、S、s₃和s₄分别代表与FAPAR估算模型运行尺度相对应的FAPAR值、地表反射率和Landsat第3、4波段的地表反射率。如FAPAR估算模型运行在MODIS尺度上，FAPAR代表MODIS FAPAR值；S、s₃和s₄分别代表Landsat地表反射率聚合至MODIS尺度(480m)的值；如FAPAR估算模型运行在Landsat尺度，则通过Landsat地表反射率s₃和s₄可计算出Landsat尺度FAPAR。

对每个地表单元分别建立MODIS FAPAR与Landsat地表反射率聚合值的线性模型，即降尺度转换函数，以此来实现对MODIS FAPAR反演模型分层线性化；通过对每一地表单元筛选样本进行回归分析，求解上式中的模型系数，即可建立尺度自洽的(即适用于MODIS和Landsat两个尺度)FAPAR先验估算模型。

基于步骤S1中由历史数据建立的筛选样本，采用回归分析方法，可以求出每一类地表单元的模型系数，得出对应地表单元的FAPAR先验估算模型。

步骤S2、基于贝叶斯框架融合MODIS FAPAR和Landsat反射率新观测数据，更新MODIS FAPAR与Landsat反射率间的尺度转换函数，得到更新的FAPAR后验估算模型。

由于FAPAR先验估算模型建立在历史数据基础上，大量数据的平均效应会导致该模型在估算单个日期FAPAR时产生偏差。此外，FAPAR具有高度的时间动态性。已有的研究表明，FAPAR与地表反射率间的关系会随太阳天顶角、大气状况、土壤背景发生变化，即FAPAR与地表反射率间的关系会随时间变化。已建立的FAPAR估算模型在用于不同日期FAPAR估算时会产生偏差。历史数据为我们提供了先验知识，仍需要融合新观测样本数据，才能实现FAPAR估算模型的实时更新，实现实时的降尺度转换。

因此，本发明实施例提出一种基于贝叶斯融合新观测的样本数据，对FAPAR先验估算模型进行实时更新的方法。

把FAPAR估算模型参数X(即降尺度转换函数参数)看成是具有某种先验分布的随机变量，可通过贝叶斯理论融合新观测的样本数据。

具体步骤是：

步骤S21、对每一景新观测的Landsat地表反射率数据，采用K均值聚类算法，将地表自动划分为林地、农田、稀疏植被、裸地和其他等五类地类。

选取其对应日期的MODIS FAPAR数据。

基于每一对新观测的Landsat地表反射率数据和MODIS FAPAR数据，对每一个地表单元，按照步骤S13中的方法筛选具有FAPAR高质量反演值的MODIS纯像元，提取其对应的MODIS FAPAR值与Landsat第3、4波段地表反射率聚合值，构成新观测数据样本。同时提取新观测MODIS FAPAR数据中提供的FAPAR反演标准差(即FAPARStd字段)。

步骤S22、将步骤S14中建立的FAPAR先验估算模型作为先验知识，按照贝叶斯公式融合新观测数据，更新FAPAR先验估算模型系数，得到FAPAR后验估算模型。具体过程如下：

将FAPAR估算模型参数X看成是具有某种先验分布的随机变量，则MODIS尺度上FAPAR的估计值F_M可由MODIS尺度反射率S_M和FAPAR估算模型参数X计算得到，有

F＝S_MX；

FAPAR估算模型参数X分布可由概率分布函数描述，也即先验分布P(X)。融合MODIS尺度上新观测信息F_M，由先验分布得到模型参数后验概率的关键是贝叶斯公式，表示为：

其中，P(X)为FAPAR估算模型参数X的先验概率密度函数，P(F_M|X)为给定FAPAR估算模型参数X时得到MODIS FAPAR新观测数据F_M的条件概率；P(F_M)＝∫_M′P(F_M|X)P(X)dX，M′为FAPAR估算模型参数X的取值空间，P(F_M)表示得到新观测值F_M的全部概率；P(X|F_M)为FAPAR估算模型参数X的后验概率，即给定MODIS FAPAR新观测数据F_M时得到模型参数X的条件概率，融合了先验知识与新观测信息，是关于X的全部知识。

假设模型参数误差与观测数据误差均服从高斯分布，根据贝叶斯定理，可通过最小化代价函数得到FAPAR估算模型参数X的最优估计，即FAPAR估算模型参数X后验分布的均值μ_X和协方差矩阵C_x，表示为：

其中，C_d为观测误差协方差矩阵，C_p是FAPAR估算模型参数X先验分布的协方差矩阵(即P(X)的协方差矩阵)，X_p为先验知识对FAPAR估算模型参数X的最优估计(即FAPAR估算模型参数先验估计值)。μ_x即为融合新观测信息得到后验估算模型参数。C_x也即μ_x的误差协方差矩阵。

为简单起见，使用球形协方差矩阵定义C_p和C_d，表示为：

其中，

代表FAPAR估算模型参数先验估计值X_p的方差，

代表MODIS FAPAR新观测数据的方差，I(r)为r×r的单位矩阵，I(m)为m×m的单位矩阵，其中，r值为FAPAR估算模型参数个数，m为新观测值数量，即对某一地表单元，基于新观测数据筛选出具有MODISFAPAR高质量反演值的MODIS纯像元数量(即筛选样本所含的像元数)。

由MODIS FAPAR数据中给出的FAPAR反演标准差(即FAPARStd字段，代表σ_d)计算得出。本实施例建立的FAPAR估算模型为Landsat数据红光、近红外波段反射率构成的线性函数，共有3个参数(含常数项)，即r＝3。

是由筛选样本以FAPAR估算模型形式进行线性回归分析中，各模型参数估计的标准差(即σ_p)计算得出。回归分析中模型参数的标准差，决定了模型参数的95％置信区间，表示模型参数估计的不确定性。由于FAPAR估算模型是在由地表覆盖类型和土壤类型划分的地表单元(即像元集合)上构建的，FAPAR估算模型的参数更新也同样基于地表单元。对每一地表单元，按照上式，通过融合新观测数据对FPAR先验估算模型参数进行更新。

步骤S3、基于更新的FAPAR后验估算模型，基于新观测Landsat数据实时估算FAPAR。

根据降尺度转换原理，更新后的FAPAR后验估算模型可直接应用于新观测Landsat地表反射率数据上相应的地表单元，进而估算得到Landsat尺度FAPAR数据。

转换结果分析：

1)与MODIS FAPAR数据对比

为验证本发明实施例提出的降尺度转换方法的精度，选择两种典型降尺度算法进行对比验证，包括多尺度回归树方法和基于NDVI的线性转换系数法。

多尺度回归树方法可直接移植至MODIS FAPAR降尺度转换，具体实现过程为：筛选MODIS尺度纯像元上的高质量反演值作为样本，建立MODIS FAPAR与Landsat聚合反射率间的线性回归树模型，进而将该模型用于估算Landsat FAPAR值。在样本构建方面，该方法与本发明实施例的降尺度转换方法存在一定差异。该算法不考虑地类和土壤类型，直接对选取的样本，按照反射率进行树划分并对叶结点建立FAPAR线性回归模型。该算法中的纯像元筛选标准仅考虑CV。为保证算法对比的基准一致，具体实现时，其和本发明实施例的方法都采用了相同的样本，即同时考虑CV和升尺度地类频率提取出的样本。

基于NDVI的线性转换系数法的具体实现过程为：计算Landsat NDVI，并将其聚合至MODIS尺度；对每个观测日期，逐像元计算聚合NDVI到MODIS FAPAR的转换系数；将该转换系数用于将Landsat NDVI转换至Landsat FAPAR。将三种降尺度转换算法得到的LandsatFAPAR与MODIS FAPAR进行对比，分析二者的一致性，作为衡量降尺度算法精度和多尺度自洽性的指标之一。首先，选择植被生长过程中的代表性日期(年积日分别为：112,176,224,272,304)，通过目视判读对比降尺度转换结果与MODIS FAPAR数据，如图3～7所示。图3～7中a均为MODIS FAPAR，b均为采用多尺度回归树方法的降尺度转换结果，c均为基于NDVI的线性转换系数法的降尺度转换结果，d均为本实施例方法的降尺度结果。由图可知，基于NDVI的线性转换系数法的结果与MODIS FAPAR空间分布的整体趋势最为一致，但空间连续性较差，表现出明显的斑块状特征，且存在较多的数据空洞。尤其是在植被生长初期和植被衰退期(如年积日为112时)，研究区西南地形起伏区域的植被分布在空间上呈现出不连续、不规则的特征。这表明，基于NDVI的线性转换系数法尽管能在理论上与MODIS FAPAR数据高度吻合，但却丢失了大量的空间细节信息，降尺度转换结果空间连续性较差，且精度直接依赖MODIS FAPAR的反演精度。本实施例方法与多尺度回归树方法得出的结果，均表现出较好的空间连续性。但多尺度回归树方法结果在植被生长期(如年积日为176和224时)表现出对FAPAR的明显低估，而本实施例算法得到的FAPAR在分布趋势和值域范围上都与MODISFAPAR更为接近。这是由于多尺度回归树方法建立在历史数据的回归分析上，大量数据的平均效应虽稳定、可靠，但如不考虑FAPAR的年际变化和季候差异，会引入额外误差。

进一步对三种算法结果与MODIS FAPAR间进行散点图分析。在产品的直接对比与验证中，通常认为取3×3像元窗口可有效减小配准误差的影响。因此，将MODIS FAPAR和Landsat FAPAR同时升尺度至1.5km，选取MODIS高质量反演数据对比二者量上的差异，结果如图8(a)～(c)所示。多尺度回归树方法的结果在FAPAR高值部分存在明显低估。基于NDVI的线性转换系数法理论上可以得到与MODIS FAPAR完全吻合的结果，但散点图显示仍有不少像元存在明显偏离。这是由于其忽略了MODIS亚像元尺度上的异质性。尤其是当MODIS亚像元尺度上存在非植被地类时，Landsat尺度上NDVI<0，而聚合至MODIS尺度NDVI接近0，会产生较大的转换系数，导致FAPAR异常值的产生，这是一种典型的尺度效应。本发明实施例的降尺度转换方法结果与MODIS FAPAR整体吻合较好，在FAPAR高值部分略有低估。

以平均绝对误差MAE(Mean Absolute Error)和均方根误差RMSE(Root MeanSquared Error)为指标，定量分析Landsat FAPAR与MODIS FAPAR间的差异，如表1所示。

表1 Landsat FAPAR与MODIS FAPAR差值统计

QC＝0像元	多尺度回归树方法	基于NDVI的线性转换系数法	本发明实施例的方法
				MAE	0.0357	0.0043	0.0275
RMSE	0.0586	0.0241	0.0454

由表1可看出，基于NDVI的线性转换系数法结果与MODIS FAPAR相比差异最小，这与其算法原理相关。其次是本发明实施例的降尺度转换方法得出的结果，而多尺度回归树方法结果与MODIS FAPAR相比差异最大。综合看来，本发明实施例的降尺度转换方法得出的结果既与MODIS FAPAR在空间趋势和量值分布上都具有较好的一致性，又表现出较好的时空连续性，具有明显优势。

2)与地面验证数据对比

地面测量数据基于独立观测手段，是验证降尺度算法精度最为有效而可靠的指标。本实施例采用的地面验证数据是经过预处理及时间归一化校正的FAPAR地面测量数据，即校正至卫星过境时刻的FAPAR测量数据。地面测量数据的可用日期分别为2012年6月8日、6月24日和7月10日。

与地面实测数据直接对比验证：

首先对三种算法结果与地面测量数据进行散点图分析，如图9(a)～(c)所示。由于MODIS FAPAR在农田反演精度较高，以MODIS FAPAR为参考的三种算法结果与地面实测数据的吻合都较好。但多尺度回归树方法在FAPAR高值部分(7月10日)表现出明显低估，基于NDVI的线性转换系数法在FAPAR高值部分(7月10日)有高估的倾向。本发明实施例的方法结果与地面实测FAPAR吻合较好，仅在FAPAR低值部分(6月8日)略有高估。

进一步分析三种算法结果与地面实测数据量值上的差异，如表2所示。以MAE和RMSE指标定量评价。

表2 Landsat FAPAR与地面实测FAPAR差值统计

评价指标	多尺度回归树方法	基于NDVI的线性转换系数法	本发明实施例的方法
				MAE	0.0693	0.0669	0.0546
RMSE	0.0894	0.0798	0.0710

由表2可看出，本发明实施例的方法精度最高，明显优于多尺度回归树方法和基于NDVI的线性转换系数法。

与ASTER FAPAR真实参考数据对比：

首先对比三种算法结果与ASTER FAPAR真实参考数据的空间分布，如图10(a)(b)-13(a)(b)所示，多尺度回归树方法结果在6月24日时与ASTER FAPAR空间分布吻合较好，但在7月10日植被生长旺盛时表现出明显低估；基于NDVI的线性转换系数法表现出明显的斑块状特征，空间分布连续性较差，尤其是在荒漠与农田交界处空间不连续特征尤为明显；本发明实施例的方法在两个日期都与ASTER FAPAR空间分布趋势非常一致，具有很好的空间连续性。

进一步由散点图对比三种算法结果与ASTER FAPAR真实参考数据量值上的差异。为避免配准误差的影响，将Landsat FAPAR和ASTER FAPAR同时升尺度至90m。如图14(a)～(c)所示，多尺度回归树方法在FAPAR低值部分有所高估，在FAPAR高值部分有明显低估。基于NDVI的线性转换系数法整体上对FAPAR存在明显高估。受尺度效应影响，大量像元在尺度转换中产生了异常的转换系数导致FAPAR值远大于1，表现出FAPAR值截断至1的特征。本发明实施例的方法在FAPAR低值部分存在一定程度上的高估，整体上与ASTER FAPAR具有高度的一致性。

分别将Landsat FAPAR和ASTER FAPAR升尺度至90m和480m，分析三种算法结果与地面实测数据量值上的差异，以MAE和RMSE指标定量评价，结果如表3所示。

表3 Landsat FAPAR与ASTER FAPAR差值统计

由表3可知，本发明实施例的方法精度最高，多尺度回归树方法次之，基于NDVI的线性转换系数法精度最低。结合表2，尽管基于NDVI的线性转换系数法与MODIS FAPAR高度吻合，但由于未考虑MODIS亚像元尺度上的异质性，导致其与地面验证数据相比精度并不理想。本发明实施例的方法既与MODIS FAPAR吻合较好，与地面验证数据相比精度也较高，且空间分布连续性较好，具有明显优势。

3)时间趋势

从时间趋势上，通过与MODIS FAPAR和地面站点实测数据对比，检验本发明实施例的方法的结果的精度。首先检验降尺度FAPAR与MODIS FAPAR时间趋势上的一致性。选择农田、林地、稀疏植被、裸地等地类的典型像元绘制时间序列曲线，如图15(a)～(d)所示。典型像元是以变异系数CV和升尺度地类频率为指标，筛选出MODIS各地类具有代表性的纯像元。在纯像元基础上，Landsat FAPAR与MODIS FAPAR吻合越好，说明降尺度转换结果与MODISFAPAR的自洽性越好。由图15(a)～(d)可知，在MODIS尺度上，各地类像元内的异质性随时间推移呈增大趋势。但无论像元内部异质性高或者低，Landsat FAPAR升尺度均值能够与MODIS FAPAR在时间趋势上很好地吻合，说明了本发明实施例的方法精度可靠。

其次，比较13个地面站点FAPAR实测数据与对应Landsat像元上降尺度FAPAR的时间趋势。如图16所示，在所有站点上，本发明实施例的方法降尺度转换得到的FAPAR均能很好地展现出作物的生长变化规律，曲线平滑、连续，且都能与地面站点实测FAPAR的时间趋势很好地吻合。

上述结果均表明：多尺度回归树算法能够保证降尺度转换结果的空间连续性和自洽性，但却未融合植被参数的时变特征，估算得到的Landsat FAPAR结果与MODIS FAPAR吻合程度最低，尤其是对植被生长旺盛时期的FAPAR存在明显低估；基于NDVI的线性转换系数法能够实现实时的降尺度转换，考虑到了植被参数的季候差异，但却是以牺牲MODIS亚像元尺度上的空间细节信息、忽略尺度效应为代价的。该方法降尺度转换结果的空间连续性较差，存在明显的斑块状特征，与地面真实参考数据相比精度也较低。本发明实施例提出的方法，既对尺度转换过程中的尺度效应进行了充分考虑，又通过融合新观测数据对FAPAR的时变特征进行了充分考虑，在此基础上实现了实时的降尺度转换。本发明实施例的降尺度转换方法得到的结果精度最高、时空连续性最好，与MODIS FAPAR和地面验证数据相比均表现出很好的一致性，且本发明实施例方法降尺度转换得到的FAPAR，在升尺度至MODIS尺度时精度(表3：RMSE＝0.0264)更优于MODIS FAPAR产品([3]：RMSE＝0.054)。既验证了发明实施例方法降尺度转换方法的优越性，也说明了构建遥感趋势面的重要意义：提高遥感产品的时空分辨率，有望提高遥感产品的精度。

[3]Y.Wang,D.Xie,S.Liu,R.Hu,Y.Li,and G.Yan,"Scaling of FAPAR from theField to the Satellite,"Remote Sensing,vol.8,p.310,2016.

MODIS FAPAR产品生产用到的算法没有显性的表达式，是由复杂的三维辐射传输模型反演而来的，但是它用途广泛、精度较高。因此，本实施例方法重点是建立MODIS FAPAR与Landsat反射率间的关系，得到Landsat尺度上FAPAR的显性表达式。为了尺度自洽，通过地表单元划分实现了对二者关系的分层线性化。方法与多尺度回归树方法有所类似，即都是建立了多元线性回归模型。但建立方法不同：第一：多尺度回归树方法等是基于样本自然划分，缺少机理，且样本难以代表所有像元，而本实施例的方法是基于图像划分地表单元，建立在FAPAR估算机理和尺度转换理论基础上的。第二：多尺度回归树方法是基于多期数据建立一个模型后，对任意一天的影像都采用同样的模型。这种方法没有考虑到二者关系的复杂性，即随时间会发生变化，所以在对某一景影像FAPAR进行估算时会产生误差，并且难以实现FAPAR的实时估算。本实施例的方法采用了贝叶斯实时融合新的数据，考虑到FAPAR估算模型随时间的变化。基于NDVI的线性转换系数法可以实时进行尺度转换，但是由于NDVI存在尺度效应，它无法实现自洽的尺度转换。并且，基于NDVI的线性转换系数法忽略了MODIS亚像元尺度上的异质性，对MODIS像元内部所有Landsat像元采用相同的转换系数，导致了结果误差较大、空间连续性较差。与本实施例的方法区别也较大。

本发明实施例基于遥感趋势面的主要思想，在尺度自洽的前提下，发展了一种基于先验知识的、实时融合新观测数据的MODIS FAPAR降尺度转换方法，通过对研究区MODISFAPAR降尺度转换的对比验证表明，本发明实施例提出的方法在算法精度、多尺度自洽性和时空连续性等方面，均优于多尺度回归树法和基于NDVI的线性转换系数法等两种典型算法。在先验知识库支持下，利用土壤类型和地表覆盖类型数据对地表单元进行合理划分，将MODIS FAPAR与Landsat地表反射率间的非线性关系在特定的地表单元转换为线性函数，避免了尺度转换过程中的尺度效应，保证了降尺度转换结果的多尺度自洽性。且在大量历史观测数据的基础上，建立了MODIS FAPAR降尺度转换函数作为先验知识，在贝叶斯方法框架下通过融合新观测数据，实现了对Landsat FAPAR的实时估算。

以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。

Claims (5)

1.基于低分辨率遥感产品降尺度的高分辨率FAPAR估算方法，其特征在于，按照以下步骤进行：

步骤S1，基于先验知识库，通过地表单元划分对MODIS FAPAR反演模型进行分层线性化，建立MODIS FAPAR与Landsat地表反射率间的尺度转换函数，得到FAPAR先验估算模型；

步骤S2，基于贝叶斯框架融合MODIS FAPAR和Landsat地表反射率新观测数据，更新MODIS FAPAR与Landsat地表反射率间的尺度转换函数，得到更新的FAPAR后验估算模型；

步骤S3，基于更新后的FAPAR后验估算模型，基于新观测Landsat数据实时估算FAPAR；

所述步骤S2的具体步骤是：

步骤S21、对每一景新观测的Landsat地表反射率数据，采用K均值聚类算法，将地表自动划分为林地、农田、稀疏植被、裸地和其他；

选取其对应日期的MODIS FAPAR数据；

基于每一对新观测的Landsat地表反射率数据和MODIS FAPAR数据，对每一个地表单元，筛选具有FAPAR高质量反演值的MODIS纯像元，提取其对应的MODIS FAPAR值与Landsat红光、近红外波段地表反射率聚合值，构成新观测数据样本；同时提取新观测MODIS FAPAR数据中提供的FAPAR反演标准差；

步骤S22、将FAPAR先验估算模型作为先验知识，按照贝叶斯公式融合新观测数据，更新FAPAR先验估算模型系数，得到FAPAR后验估算模型，具体是：

将FAPAR估算模型参数X看成是具有某种先验分布的随机变量，则MODIS尺度上FAPAR的估计值F_M由MODIS尺度反射率S_M和FAPAR估算模型参数X计算得到，有

F_M＝S_MX；

FAPAR估算模型参数X分布由概率分布函数描述，也即先验分布P(X)；融合MODIS尺度上新观测信息F_M，由先验分布得到模型参数后验概率的关键是贝叶斯公式，表示为：

其中，P(X)为FAPAR估算模型参数X的先验概率密度函数，P(F_M|X)为给定FAPAR估算模型参数X时得到MODIS FAPAR新观测数据F_M的条件概率；P(F_M)＝∫_M′P(F_M|X)P(X)dX，M′为FAPAR估算模型参数X的取值空间，P(F_M)表示得到新观测值F_M的全部概率；P(X|F_M)为FAPAR估算模型参数X的后验概率，即给定MODIS FAPAR新观测数据F_M时得到模型参数X的条件概率，融合了先验知识与新观测信息，是关于X的全部知识；

假设模型参数误差与观测数据误差均服从高斯分布，根据贝叶斯定理，通过最小化代价函数得到FAPAR估算模型参数X的最优估计，即FAPAR估算模型参数X后验分布的均值μ_X和协方差矩阵C_x为：

其中，C_d为观测误差协方差矩阵，C_p是FAPAR估算模型参数X先验分布的协方差矩阵，X_p为先验知识对FAPAR估算模型参数X的最优估计，即FAPAR估算模型参数先验估计值；

使用球形协方差矩阵定义C_p和C_d，表示为：

其中，

代表FAPAR估算模型参数先验估计值X_p的方差，

代表MODIS FAPAR新观测数据的方差，I(r)为r×r的单位矩阵，I(m)为m×m的单位矩阵，其中，r值为FAPAR估算模型参数个数，m为新观测值数量，即对每一地表单元，基于新观测数据筛选出具备高质量MODISFAPAR反演值的MODIS纯像元数量；

由MODIS FAPAR数据中给出的FAPAR反演标准差σ_d计算得出；建立的FAPAR估算模型为Landsat数据红光、近红外波段反射率构成的线性函数，共有3个参数，即r＝3；

是由筛选样本以FAPAR估算模型形式进行线性回归分析中，各模型参数估计的标准差σ_p计算得出；基于土壤类型和新观测地表覆盖类型划分的地表单元，按照上式，通过融合新观测数据对FPAR先验估算模型参数进行更新。

2.根据权利要求1所述的基于低分辨率遥感产品降尺度的高分辨率FAPAR估算方法，其特征在于，所述步骤S1按照以下步骤进行：

步骤S11，数据收集与预处理；

使用Landsat地表反射率数据，影像空间分辨率30m，采用可见光、近红外波段，筛选历史年限和当年云量小于10％的高质量地表反射率数据；历史年限数据用于先验知识构建，当年的数据作为新观测信息；地表反射率数据产品空间基准为WGS-84坐标系和UTM投影，对数据进行解压缩、裁剪和数值转换；

使用MODIS FAPAR产品，在MODIS分辨率尺度上，筛选出与Landsat影像观测日期对应的MODIS FAPAR产品，对产品数据进行数据层提取、裁剪、投影转换、重采样、数值转换；

步骤S12，选择研究区具有代表性的一景影像，基于Landsat影像分析，采用K均值聚类算法将研究区的Landsat影像自动划分为林地、农田、稀疏植被、裸地和其他五种地类类别，将Landsat分类结果按照众数规则升尺度至1km尺度，得到分辨率1km的地表覆盖类型数据；使用土壤类型栅格数据，对土壤类型数据和地表覆盖类型进行空间上的交集操作，将研究区划分为相应的地表单元，每一地表单元具有特定的土壤类型和地表覆盖类型；

步骤S13，筛选MODIS FAPAR高质量纯像元样本；

包括(1)，筛选MODIS FAPAR高质量反演值：选用MODIS主算法反演的FAPAR高质量反演数据建立样本，筛选标准为MODIS FAPAR质量控制标识QC＝0；

(2)，MODIS纯像元筛选：采用升尺度地类频率和变异系数CV来度量MODIS亚像元尺度上地表异质性，进行纯像元筛选；变异系数定义为标准差与均值的比值；按照CV≤0.2且升尺度地类频率≥0.95的标准筛选纯像元；在纯像元筛选中仅使用植被生长期的数据；

对每景Landsat地表反射率数据进行如下处理：

a)，将Landsat地表反射率升尺度至MODIS 480m分辨率尺度，并计算变异系数CV：首先，对Landsat可见光-近红外波段的每一个波段，计算每一MODIS像元的波段CV值，即该MODIS像元空间范围包含的16×16个Landsat像元波段反射率的标准差与均值的比值；其次，计算每一个MODIS像元各个波段CV值的均值，作为该像元的CV值；

b)，将Landsat地表分类数据升尺度至480m，按照众数规则得出MODIS像元尺度上的优势地类，即升尺度像元对应的地类；同时，计算该优势地类在MODIS亚像元尺度上所占的面积比例，即升尺度地类频率；

(3)，对历史年限的每一景Landsat地表反射率数据及MODIS FAPAR数据，按照步骤(1)-(2)中的过程，筛选出具有MODIS FAPAR高质量反演值的纯像元；然后，根据划分的地表单元，筛选出每一类地表单元对应的MODIS FAPAR高质量反演值的纯像元；提取所有筛选像元对应的MODIS FAPAR值和聚合至480m分辨率的Landsat可见光和近红外波段的地表反射率值，即构成每一地表单元的筛选样本；

步骤S14，对每一地表单元，基于筛选样本，建立MODIS FAPAR与Landsat地表反射率间的尺度转换函数，通过回归分析建立FAPAR先验估算模型。

3.根据权利要求2所述的基于低分辨率遥感产品降尺度的高分辨率FAPAR估算方法，其特征在于，所述步骤S14按照以下步骤进行：

MODIS FAPAR值与地表反射率间的关系表示为：

其中，下标M代表MODIS尺度；FAPAR_M代表MODIS FAPAR值；S_M代表MODIS尺度的地表反射率；f_M表示MODIS尺度的FAPAR估算模型，对于MODIS产品，f_M代表一个复杂的物理模型；下标L代表Landsat尺度，i表示给定MODIS像元内对应Landsat像元的索引，S_L(i)表示MODIS像元内第i个Landsat像元的地表反射率；n代表MODIS像元对应Landsat像元的总数；

结合土壤类型和地表覆盖类型进行地表单元划分，对每一地表单元建立特定的FAPAR估算模型，形式为：

其中，下标soil代表土壤类型，下标lc代表地表覆盖类型，

表示特定土壤类型和地表覆盖类型条件下MODIS尺度FAPAR估算模型；

当

为线性模型时，有

表示特定土壤类型和地表覆盖类型条件下Landsat尺度FAPAR估算模型；建立在MODIS尺度上的FAPAR线性估算模型

适用于Landsat尺度FAPAR的估算，即建立了尺度自洽的FAPAR估算模型，则Landsat尺度FAPAR的估算FAPAR_L由下式计算：

其中，s_j表示第j波段的地表反射率，a_jsoil,lc表示线性函数

第j波段的系数，k为波段总数，

为线性函数

的常数项；

Landsat尺度上，对同类地表单元，FAPAR估算采用相同的线性模型

则MODIS纯像元上的FAPAR由下式计算：

式中，s_i,j为MODIS像元内对应第i个Landsat像元第j波段的地表反射率；

为该MODIS像元内所有Landsat像元第j波段地表反射率的均值，也就是Landsat像元第j波段地表反射率聚合至MODIS 480m分辨率的均值；由上式建立MODIS FAPAR与Landsat地表反射率之间的联系，即在已知FAPAR_M和s_i,j的条件下，求出尺度自洽的FAPAR估算模型

对每一地表单元，筛选具备MODIS FAPAR高质量反演值的MODIS纯像元，即可在已知MODIS FAPAR和Landsat多波段的地表反射率聚合值的情况下，通过回归分析，求解该地表单元对应的FAPAR线性估算模型的系数

和

得出该地表单元对应的FAPAR先验估算模型。

4.根据权利要求3所述的基于低分辨率遥感产品降尺度的高分辨率FAPAR估算方法，其特征在于，所述步骤S14中，采用基于Landsat红光、近红外波段反射率的二元回归模型建立尺度转换函数，也即尺度自洽的FAPAR估算模型；模型形式表示为：

式中，S代表多波段地表反射率；s₃和s₄分别为Landsat红光、近红外波段的地表反射率；

和

分别表示线性函数f_soil,lc的常数项、变量s₃和s₄的系数，即模型系数；f_soil,lc表示尺度自洽的FAPAR估算模型，即同时适用于MODIS和Landsat两个尺度的FPAR估算模型；FAPAR表示与估算模型运行尺度对应尺度即MODIS或Landsat尺度上的FAPAR值；

对每个地表单元分别建立MODIS FAPAR与Landsat地表反射率聚合值的线性模型，即降尺度转换函数，以此来实现对MODIS FAPAR反演模型分层线性化；通过对每一地表单元筛选样本进行回归分析，求解上式中的模型系数，即可建立尺度自洽的FAPAR先验估算模型。

5.根据权利要求2-4任一项所述的基于低分辨率遥感产品降尺度的高分辨率FAPAR估算方法，其特征在于，步骤S21中，按照步骤S13，筛选具有FAPAR高质量反演值的MODIS纯像元，提取其对应的MODIS FAPAR值与Landsat红光、近红外波段地表反射率聚合值，构成新观测数据样本。

Images