本发明属于遥感估算领域，涉及一种基于低分辨率遥感产品降尺度的高分辨率fapar估算方法。
背景技术：
：光合有效辐射吸收比率(fractionofabsorbedphotosyntheticallyactiveradiation,fapar/fpar)定义为植被冠层(绿色部分)吸收的光合有效辐射(photosyntheticallyactiveradiation,par)占总入射par的比率，是表征植被光合作用的核心参数。一方面，fapar直接反映大气圈与生物圈间的物质交换过程，是建立碳循环模型、陆地生态系统模型的关键变量，也是植被生长状态与土地覆盖监测中的重要指标。另一方面，fapar反映植被冠层对光能的截获利用情况，是地气系统辐射平衡与能量交换的重要影响因素。因此，fapar的准确估算，是定量估算光合作用、提高全球生态系统碳估算精度的基础，在全球气候变化、碳排放、粮食安全等重大问题研究方面有着重要的应用价值。尽管遥感估算已成为获取面状fapar信息的唯一有效手段，但目前全球及区域尺度上通用的fapar产品间精度差别较大，且都是在低分辨率尺度上生产的，丢失了地表大量的空间细节信息。公里级产品亚像元尺度上的地表异质性，不仅会影响fapar的反演精度，造成多尺度fapar产品间难以自洽，还直接影响生态系统生产力和植被生物量的估算。获取高空间分辨率fapar遥感产品，对于改进fapar遥感产品的精度和空间连续性，更好地满足生态学应用需求具有重要意义。目前，获取高空间分辨率fapar有两种方式：一是采用基于植被指数的半经验方法或基于机理模型的物理反演方法，直接由高分辨率遥感数据估算fapar；二是以成熟的业务化fapar遥感产品为参考，通过降尺度转换得到高分辨率尺度上的fapar。然而，无论是采用基于植被指数的半经验方法或基于机理模型的物理反演方法，在较大范围上应用时都会面临算法普适性、参数可获取性、多尺度自洽性和真实性检验等多重挑战。相比直接反演方法，以modis(moderate-resolutionimagingspectroradiometer)为代表的业务化遥感产品建立在严格的物理模型基础上，在全球范围积累了大量的历史观测数据，并在不断产生新观测数据，算法及产品经过了数十年来的广泛检验，精度较为可靠。以成熟的业务化fapar遥感产品为参考进行降尺度转换，成为获取高分辨率fapar的一种可行方式。t.hwang,c.song,p.v.bolstad,andl.e.band,"downscalingreal-timevegetationdynamicsbyfusingmulti-temporalmodisandlandsatndviintopographicallycomplexterrain,"remotesensingofenvironment,vol.115,pp.2499-2512,2011提出了一种modisfapar降尺度转换的方法，该方法假设modisfapar与landsatndvi(归一化植被指数，normalizeddifferencevegetationindex)间的关系具有尺不变性，逐像元计算modisfapar与landsatndvi聚合值的转换系数，再将其用于从landsatndvi中估算fapar(即基于ndvi的线性转换系数法)。由于转换系数是在modis像元尺度上逐日计算的，该方法考虑到了降尺度转换系数的季候差异，但却忽略了modis亚像元尺度上的空间异质性。另外，ndvi的尺度效应已被大量研究证实，fapar与ndvi间转换关系的尺不变假设显然会在异质性地表引入尺度误差，导致结果难以自洽。为在降尺度转换过程中尽量避免尺度效应，f.gaoandw.p.kustas,"simplemethodforretrievingleafareaindexfromlandsatusingmodisleafareaindexproductsasreference,"journalofappliedremotesensing,vol.6,p.063554,2012提出了一种基于多尺度回归树模型的modislai(叶面积指数，leafareaindex)降尺度方法，也即一些现有专利提到的cubist方法。该方法利用多期landsat反射率数据和modislai产品，选择纯像元样本在modis尺度上建立反射率与lai的多尺度线性回归树模型，并将该模型用于landsat尺度上lai的估算。采用线性模型可以避免在降尺度转换中产生尺度效应。但该方法对全年所有日期数据采用相同的转换关系，忽略了lai与反射率间的关系存在季候差异。这对于时间变化特征较强的植被参数(如fapar)，显然会引入额外误差。尤其是当需要对modis产品进行长时间序列上、实时的降尺度转换时，该方法难以对新观测数据进行融合更新，而是每次都需选择当年的多期数据来建立转换关系，难以达成实时的降尺度转换。总结现有的两种典型算法，基于ndvi的线性转换系数法能够实现实时的降尺度转换，考虑到了植被参数的季候差异，但却是以牺牲modis亚像元尺度上的空间细节信息、忽略尺度效应为代价的。基于多尺度回归树模型的modislai降尺度方法能够保证结果的空间连续性和自洽性，但却忽略了植被参数的时变特征，难以达成实时的降尺度转换。从先验知识利用的角度分析，基于ndvi的线性转换系数法只利用了时间上的先验知识，而基于多尺度回归树模型的modislai降尺度方法仅利用空间上的先验知识。此外，基于ndvi的线性转换系数法是将新观测作为数据进行一次性的利用，基于多尺度回归树模型的modislai降尺度方法是建立在先验知识之上的，但却难以对先验知识进行更新。因此，需要在多尺度自洽的前提下，发展更为普适的、可综合利用各类先验知识并可有效更新的遥感数据降尺度转换方法。技术实现要素：本发明实施例的目的在于提供一种基于低分辨率遥感产品降尺度的高分辨率fapar估算方法，以解决现有估算方法难以同时实现尺度自洽、空间连续、实时的高分辨率fapar估算的问题。本发明所采用的技术方案是，基于低分辨率遥感产品降尺度的高分辨率fapar估算方法，按照以下步骤进行：步骤s1，基于fapar遥感估算机理、遥感产品尺度转换理论、历史年限modisfapar产品与landsat地表反射率数据、相关地理要素数据(土壤类型与地表覆盖类型)等先验知识，通过地表单元划分，对modisfapar反演模型进行分层线性化，建立modisfapar与landsat地表反射率间的尺度转换函数(即尺度自洽的fapar先验估算模型)；步骤s2，基于贝叶斯框架融合modisfapar和landsat地表反射率新观测数据，更新modisfapar与landsat地表反射率间的尺度转换函数，得到更新的fapar后验估算模型；步骤s3，基于更新后的fapar后验估算模型，基于新观测landsat数据实时估算fapar。本发明实施例的有益效果是：通过地表覆盖类型和土壤类型对地表单元进行合理划分，针对每一地表单元建立fapar与多波段反射率间的线性估算模型，实现了对低分辨率fapar业务化遥感产品生产所用的复杂辐射传输模型的线性化，以避免尺度效应，实现了自洽的尺度转换；对每一地表单元，通过选取纯像元样本，建立modisfapar与landsat地表反射率聚合值之间的尺度转换函数，采用统计分析方法求解模型参数，建立了尺度自洽的fapar估算模型，再将该模型用于高分辨率(即landsat尺度)像元尺度，实现了空间连续的高分辨率fapar估算；同时，在遥感趋势面的框架下，利用多期历史数据构建fapar先验估算模型，再由贝叶斯定理融合新观测信息(即每一景实时获取的modisfapar和landsat遥感影像)，对先验估算模型参数进行更新，得到融合了先验知识和新观测信息的fapar估算模型，其同时利用了时间、空间上的信息，是关于高分辨率fapar全部的知识，实现了实时的高分辨率fapar估算。通过上述方法同时实现了尺度自洽、空间连续、实时的高分辨率fapar估算，有效解决了现有估算方法难以同时实现尺度自洽、空间连续、实时的高分辨率fapar估算的问题。附图说明为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。图1是本发明实施例的modisfapar遥感数据的空间降尺度转换方法流程图。图2(a)是本发明实施例采用不同回归模型进行回归分析r2指标拟合精度图。图2(b)是本发明实施例采用不同回归模型进行回归分析rmse指标拟合精度图。图3是采用三种降尺度转换算法得到的2012年年积日为112的landsatfapar结果与modisfapar对比图。图4是采用三种降尺度转换算法得到的2012年年积日为176的landsatfapar结果对比图。图5是采用三种降尺度转换算法得到的2012年年积日为224的landsatfapar结果对比图。图6是采用三种降尺度转换算法得到的2012年年积日为272的landsatfapar结果对比图。图7是采用三种降尺度转换算法得到的2012年年积日为304的landsatfapar结果对比图。图3-7中，a为modisfapar，b为多尺度回归树算法结果，c为基于ndvi的线性转换系数法结果，d为本发明实施例的方法所得结果。图8(a)是采用多尺度回归树方法所得结果与modisfapar的散点图。图8(b)是采用基于ndvi的线性转换系数法所得结果与modisfapar的散点图。图8(c)是采用本发明实施例的方法所得结果与modisfapar的散点图。图9(a)是采用多尺度回归树方法所得结果与地面测量数据的散点图。图9(b)是采用基于ndvi的线性转换系数法所得结果与地面测量数据的散点图。图9(c)是采用本发明实施例的方法所得结果与地面测量数据的散点图。图10(a)是2012年6月24日asterfapar真实参考数据的空间分布图。图10(b)是2012年7月10日asterfapar真实参考数据的空间分布图。图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)是采用多尺度回归树方法所得结果与asterfapar的散点图。图14(b)是采用基于ndvi的线性转换系数法所得结果与asterfapar的散点图。图14(c)是采用本发明实施例的方法所得结果与asterfapar的散点图。图15(a)是农田典型像元时间序列曲线图。图15(b)是林地典型像元时间序列曲线图。图15(c)是稀疏植被典型像元时间序列曲线图。图15(d)是裸地典型像元时间序列曲线图。图16是各地面站点实测fapar与发明实施例的方法的反演fapar的时间趋势对比图。图17(a)是研究区landsat影像数据自动分类图。图17(b)是研究区土地覆盖类型图。图17(c)是研究区地表单元图，其中图例代表地表单元类型代码，1-90对应特定土壤类型与地表覆盖类型组合的地表单元。具体实施方式下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。本发明实施例的基于现有业务化低分辨率fapar遥感产品降尺度的高空间分辨率fapar估算方法，可适用于各种高空间分辨率遥感数据的fapar估算，如对landsat5/thematicmapper(tm),landsat7enhancedthematicmapperplus(etm ),landsat8/operationallandimager(oli),advancedspacebornethermalemissionandreflectionradiometer(aster),hj-1/cdd,sentinel-2msi等空间分辨率不低于30m的遥感数据均具有普遍的适用性。本实施例以modisfapar和landsat遥感数据为例。但把landsat换成以上这些高空间分辨率数据同样适用。本发明实施例提出一种基于modisfapar降尺度的landsatfapar估算方法，以对2012年的landsatfapar进行估算为例进行说明，具体流程如图1所示，按照以下步骤进行：步骤s1、基于先验知识库(fapar遥感估算机理、遥感数据尺度转换理论、2000～2011年的modisfapar历史数据和landsat历史影像、地表覆盖类型数据、土壤类型数据)，通过地表单元划分对modisfapar反演模型进行分层线性化，建立modisfapar与landsat地表反射率间的尺度转换函数(即尺度自洽的fapar先验估算模型)，具体步骤如下：步骤s11、数据收集与预处理。本发明实施例使用的数据主要为landsat5tm和landsat7etm 地表反射率数据。landsattm/etm 数据具有多个波段，影像空间分辨率30m。landsat卫星重访周期为16天。受天气状况、质量因素等影响，实际可用的landsat数据时间分辨率远低于16天。本发明实施例共筛选了2000年2月至2011年12月间云量小于10％的landsattm\etm 地表反射率数据共103景，2012年高质量landsattm\etm 地表反射率数据14景。其中，历史数据(2000-2011年)平均时间分辨率约为9景/年，用于从历史观测中归纳、总结提取先验知识；2012年的数据作为新观测信息。下载的landsattm\etm地表反射率数据产品空间基准为wgs-84坐标系和utm投影，对数据进行解压缩、裁剪和数值转换等预处理。本发明实施例使用的modisfapar产品为最新发布的c6(collection6)版本。相较于前一版本，c6版本的主要改进在于空间分辨率从1km提高到500m。本发明实施例收集了2000年2月至2012年12月的583景mod15a2hc6产品数据，空间分辨率约500m，所使用的参数主要包括fapar及相关质量标识(lai\faparqc,faparstd)等。参照modislai\fapar产品手册，对mod15a2h进行数据层提取、裁剪，投影转换(至wgs-84坐标系utm投影)，重采样(至480m)、数值转换等处理。由于landsat影像和modisfapar产品具有不同的时间分辨率，首先筛选出与landsat影像观测日期对应的modisfapar产品。步骤s12、按照植被生长旺盛期、无云或少云、成像质量较好等标准，选择研究区具有代表性的一景影像(本实施例选用2012年7月10日的landsat影像)，采用k均值聚类算法将研究区的landsat影像自动划分为林地、农田、稀疏植被、裸地和其他(如云)等五种地类类别，将landsat分类结果按照众数规则升尺度至1km尺度，即得到了分辨率1km的地表覆盖类型数据。基于世界土壤数据库(harmonizedworldsoildatabase,hwsd)的中国土壤数据集(v1.1)(空间分辨率1km)，使用土壤类型栅格数据，对土壤类型数据和地表覆盖类型数据进行空间上的交集操作，将研究区划分为相应的地表单元，每一地表单元具有特定的土壤类型和地表覆盖类型。即对土壤类型数据(共18个类别)和地表覆盖类型(共5种)进行空间上的交集操作，理论上可将研究区划分为18×5种地表单元。由于某些土壤类型单元并不覆盖所有地表覆盖类型，研究区实际上共划分为73个地表单元，如图17(a)-(c)所示。土地覆盖类型代表地表覆盖的一般状态。在地表覆盖类型未发生剧烈变化(如人为干扰、自然灾害等)的情况下，只需要对一景代表性的landsat影像进行地表覆盖分类即可。如果地表覆盖类型发生了剧烈变化，可以有针对性地、在发生变化的日期后选择一景landsat影像分析、重新划分地表单元；在这种情况下，按照发生变化前和发生变化后两个时期，可有两个不同的地表单元数据，在数据筛选与分析中，按照成像时间(即是否已经发生了变化)将landsat数据和modis数据对应至相应的地表单元数据。步骤s13、筛选modisfapar高质量纯像元样本。(1)，筛选modisfapar高质量反演值：modisfapar产品生产基于主算法和备用算法两种算法。其中，主算法基于复杂三维辐射传输模型的查找表算法，缺少显性表达式，但具有严格的物理基础，精度较高。本发明实施例的目的也即：将具有高精度的modisfapar主算法进行分层线性化，得到其显性表达式。因此，只选用modis主算法反演的fapar高质量反演数据建立样本。筛选标准为modisfapar质量控制标识qc＝0。(2)，modis纯像元筛选：modis尺度上的纯像元是建立modisfapar与landsat地表反射率间尺度转换函数的关键。本发明实施例采用升尺度地类频率和变异系数cv(coefficientofvariation)两种指标来度量modis亚像元尺度上地表异质性，进行纯像元筛选。变异系数cv定义为标准差与均值的比值，cv值表示地表反射率的均匀程度，值越小，像元越均质。升尺度地类频率代表优势地类在亚像元尺度上所占的面积比例，值越大，像元越均质。为保证筛选出的像元尽可能均质，且能够包含足够多的样本，本实施例按照cv≤0.2且升尺度地类频率≥0.95的标准筛选纯像元。同时，为避免在非植被生长季选入过多的裸地像元作为样本，影响样本的代表性，在纯像元筛选中仅使用植被生长期(4至10月)的数据。具体地，利用landsat数据计算纯像元筛选的指标。对每景landsat地表反射率数据进行如下处理：a)，将landsat地表反射率升尺度至modis尺度(480m分辨率)，并计算变异系数cv(coefficientofvariation)。变异系数定义为标准差与均值的比值。具体计算过程为：首先，对landsat可见光-近红外波段的每一个波段，计算每一modis像元的波段cv值，即该modis像元空间范围包含的16*16个landsat像元波段反射率的标准差与均值的比值；其次，计算每一个modis像元各个波段(即landsat可见光-近红外的6个波段)cv值的均值，作为该像元的cv值。b)，将landsat地表分类数据升尺度至480m，按照众数规则得出modis像元尺度上的优势地类，即升尺度像元对应的地类。同时，计算该优势地类在modis亚像元尺度上所占的面积比例，即升尺度地类频率。(3)，对2000-2011年间的每一景landsat地表反射率数据及modisfapar数据，按照步骤(1)-(2)中的过程，筛选出具有modisfapar高质量反演值的纯像元。然后，根据s12中划分的地表单元，筛选出每一类地表单元对应的modisfapar高质量反演值的纯像元；提取所有筛选像元对应的modisfapar值和聚合至480m分辨率的landsat可见光和近红外波段的地表反射率值，即构成每一地表单元的筛选样本。步骤s14、建立modisfapar与landsat地表反射率间自洽的尺度转换函数，对每一地表单元，通过对筛选样本进行回归分析得出fapar先验估算模型的参数。具体过程如下：modisfapar产品是基于冠层三维辐射传输模型，采用查找表反演策略生产的。modisfapar值与地表反射率间的关系可表示为：其中，下标m代表modis尺度；faparm代表modisfapar值；sm代表modis尺度的地表反射率；fm表示modis尺度的fapar估算模型，既可以代表简单的经验线性模型，也可以代表复杂的物理模型，对于modis产品，fm代表一个复杂的物理模型；下标l代表landsat尺度，i表示给定modis像元内对应landsat像元的索引，sl(i)表示modis像元内第i个landsat像元的地表反射率；n代表modis像元对应landsat像元的总数。本发明实施例使用的modisfapar数据空间分辨率480m，一个modis像元内部对应16*16个landsat像元。modisfapar空间降尺度的问题可表述为：如何在已知faparm和sl(i)的情况下，求出fm，且fm能够适用于landsat尺度。即fm＝fl，fl为landsat尺度fapar估算模型。由尺度效应产生的机理分析[1,2]，若fm为非线性模型，则将fm不加修正地用于landsat尺度时会产生尺度效应，则fm≠fl，无法实现自洽的尺度转换。解决这个问题的困难还在于，对modisfapar产品，fm代表一个复杂的物理模型，不具备显性的数学表达式。因此，需要由fapar的遥感反演机理入手，将modisfapar产品所代表的物理反演模型进行线性化，表示为地表反射率的线性组合。文献[1]h.wuandz.li,"scaleissuesinremotesensing:areviewonanalysis,processingandmodeling,"sensors,vol.9,pp.1768-1793,2009.文献[2]z.huands.islam,"aframeworkforanalyzinganddesigningscaleinvariantremotesensingalgorithms,"ieeetransactionongeoscienceandremotesensing,vol.13,pp.747-755,1997.modisfapar产品算法考虑了植被结构、地表类型和土壤背景等多种因素间的相互作用，且现有研究表明在建立线性回归模型时考虑地表类型可以减小地表参数估算中的不确定性。因此，fapar与地表反射率间的关系与地表类型、土壤背景相关。本实施例结合土壤类型和地表覆盖类型进行地表单元划分，对每一地表单元建立特定的fapar估算模型，形式为：其中，下标soil代表土壤类型，下标lc代表地表覆盖类型。表示特定土壤类型和地表覆盖类型条件下modis尺度fapar估算模型。由上述分析可知，当为线性模型时，有表示特定土壤类型和地表覆盖类型条件下landsat尺度fapar估算模型。建立在modis尺度上的fapar估算模型可适用于landsat尺度fapar的估算，即建立了尺度自洽的fapar估算模型，则landsat尺度fapar的估算faparl可由下式计算：其中，sj表示第j波段的地表反射率，表示线性函数第j波段的系数，k为波段总数，为线性函数的常数项。landsat尺度上，对同类地表单元，fapar估算采用相同的线性模型则modis纯像元上的fapar可由下式计算：式中，si,j为modis像元内对应第i个landsat像元第j波段的地表反射率；为该modis像元内所有landsat像元第j波段地表反射率的均值，也就是landsat像元第j波段地表反射率聚合至modis480m分辨率的均值。因此，对每一地表单元，通过筛选具备modisfapar高质量反演值的modis纯像元，即可在已知modisfapar和landsat多波段的地表反射率聚合值(即聚合至modis480m分辨率的地表反射率均值)的情况下，通过回归分析，求解该地表单元对应的fapar线性估算模型的系数和得出该地表单元对应的fapar先验估算模型。根据采用的波段组合不同，降尺度转换函数可有多种形式，本实施例选用三种模型形式进行对比分析，依次为：基于landsat2-7波段反射率的多元回归模型、基于landsat3-4波段反射率的二元回归模型，以及基于ndvi的一元回归模型。分别采用三种模型对提取的modisfapar纯像元高质量反演样本进行回归分析，以决定系数r2和均方根误差rmse(rootmeansquarederror)为指标对比模型拟合精度，如图2(a)～(b)所示，可知，基于landsat2-7波段反射率的多元回归模型拟合r2最高，rmse最低，精度最高。基于landsat3、4波段反射率的二元回归模型拟合精度次之，且明显优于基于landsatndvi的一元线性回归模型。相较而言，基于landsat3、4波段反射率的二元回归模型在rmse上与基于landsat2-7波段反射率的多元回归模型差异很小，二者间最大差异为0.0095。landsat3-4波段对应红光、近红外波段，与植被特征关系最为密切，且基于landsat3、4波段反射率的二元回归模型模拟精度高于基于landsatndvi的一元线性回归模型，模型复杂度低于基于landsat2-7波段反射率的多元回归模型，因此综合考虑模型模拟精度、模型复杂度和模型尺度效应等因素，本实施例采用基于landsat3、4波段反射率的二元回归模型作为fapar估算模型的基本形式，表示为：式中，s代表地表反射率；s3和s4分别为landsat第3(红光)、4(近红外)波段的地表反射率；和分别表示线性函数fsoil,lc的常数项、变量s3和s4的系数，即模型系数。由于该模型适用于modis和landsat两个尺度，这里不再用下标m或l区分，代表fapar估算的通用模型。fsoil,lc表示尺度自洽的fapar估算模型，即同时适用于modis和landsat两个尺度的fpar估算模型。对应地，fapar、s、s3和s4分别代表与fapar估算模型运行尺度相对应的fapar值、地表反射率和landsat第3、4波段的地表反射率。如fapar估算模型运行在modis尺度上，fapar代表modisfapar值；s、s3和s4分别代表landsat地表反射率聚合至modis尺度(480m)的值；如fapar估算模型运行在landsat尺度，则通过landsat地表反射率s3和s4可计算出landsat尺度fapar。对每个地表单元分别建立modisfapar与landsat地表反射率聚合值的线性模型，即降尺度转换函数，以此来实现对modisfapar反演模型分层线性化；通过对每一地表单元筛选样本进行回归分析，求解上式中的模型系数，即可建立尺度自洽的(即适用于modis和landsat两个尺度)fapar先验估算模型。基于步骤s1中由历史数据建立的筛选样本，采用回归分析方法，可以求出每一类地表单元的模型系数，得出对应地表单元的fapar先验估算模型。步骤s2、基于贝叶斯框架融合modisfapar和landsat反射率新观测数据，更新modisfapar与landsat反射率间的尺度转换函数，得到更新的fapar后验估算模型。由于fapar先验估算模型建立在历史数据基础上，大量数据的平均效应会导致该模型在估算单个日期fapar时产生偏差。此外，fapar具有高度的时间动态性。已有的研究表明，fapar与地表反射率间的关系会随太阳天顶角、大气状况、土壤背景发生变化，即fapar与地表反射率间的关系会随时间变化。已建立的fapar估算模型在用于不同日期fapar估算时会产生偏差。历史数据为我们提供了先验知识，仍需要融合新观测样本数据，才能实现fapar估算模型的实时更新，实现实时的降尺度转换。因此，本发明实施例提出一种基于贝叶斯融合新观测的样本数据，对fapar先验估算模型进行实时更新的方法。把fapar估算模型参数x(即降尺度转换函数参数)看成是具有某种先验分布的随机变量，可通过贝叶斯理论融合新观测的样本数据。具体步骤是：步骤s21、对每一景新观测的landsat地表反射率数据，采用k均值聚类算法，将地表自动划分为林地、农田、稀疏植被、裸地和其他等五类地类。选取其对应日期的modisfapar数据。基于每一对新观测的landsat地表反射率数据和modisfapar数据，对每一个地表单元，按照步骤s13中的方法筛选具有fapar高质量反演值的modis纯像元，提取其对应的modisfapar值与landsat第3、4波段地表反射率聚合值，构成新观测数据样本。同时提取新观测modisfapar数据中提供的fapar反演标准差(即faparstd字段)。步骤s22、将步骤s14中建立的fapar先验估算模型作为先验知识，按照贝叶斯公式融合新观测数据，更新fapar先验估算模型系数，得到fapar后验估算模型。具体过程如下：将fapar估算模型参数x看成是具有某种先验分布的随机变量，则modis尺度上fapar的估计值fm可由modis尺度反射率sm和fapar估算模型参数x计算得到，有f＝smx；fapar估算模型参数x分布可由概率分布函数描述，也即先验分布p(x)。融合modis尺度上新观测信息fm，由先验分布得到模型参数后验概率的关键是贝叶斯公式，表示为：其中，p(x)为fapar估算模型参数x的先验概率密度函数，p(fm|x)为给定fapar估算模型参数x时得到modisfapar新观测数据fm的条件概率；p(fm)＝∫m′p(fm|x)p(x)dx，m′为fapar估算模型参数x的取值空间，p(fm)表示得到新观测值fm的全部概率；p(x|fm)为fapar估算模型参数x的后验概率，即给定modisfapar新观测数据fm时得到模型参数x的条件概率，融合了先验知识与新观测信息，是关于x的全部知识。假设模型参数误差与观测数据误差均服从高斯分布，根据贝叶斯定理，可通过最小化代价函数得到fapar估算模型参数x的最优估计，即fapar估算模型参数x后验分布的均值μx和协方差矩阵cx，表示为：其中，cd为观测误差协方差矩阵，cp是fapar估算模型参数x先验分布的协方差矩阵(即p(x)的协方差矩阵)，xp为先验知识对fapar估算模型参数x的最优估计(即fapar估算模型参数先验估计值)。μx即为融合新观测信息得到后验估算模型参数。cx也即μx的误差协方差矩阵。为简单起见，使用球形协方差矩阵定义cp和cd，表示为：其中，代表fapar估算模型参数先验估计值xp的方差，代表modisfapar新观测数据的方差，i(r)为r×r的单位矩阵，i(m)为m×m的单位矩阵，其中，r值为fapar估算模型参数个数，m为新观测值数量，即对某一地表单元，基于新观测数据筛选出具有modisfapar高质量反演值的modis纯像元数量(即筛选样本所含的像元数)。由modisfapar数据中给出的fapar反演标准差(即faparstd字段，代表σd)计算得出。本实施例建立的fapar估算模型为landsat数据红光、近红外波段反射率构成的线性函数，共有3个参数(含常数项)，即r＝3。是由筛选样本以fapar估算模型形式进行线性回归分析中，各模型参数估计的标准差(即σp)计算得出。回归分析中模型参数的标准差，决定了模型参数的95％置信区间，表示模型参数估计的不确定性。由于fapar估算模型是在由地表覆盖类型和土壤类型划分的地表单元(即像元集合)上构建的，fapar估算模型的参数更新也同样基于地表单元。对每一地表单元，按照上式，通过融合新观测数据对fpar先验估算模型参数进行更新。步骤s3、基于更新的fapar后验估算模型，基于新观测landsat数据实时估算fapar。根据降尺度转换原理，更新后的fapar后验估算模型可直接应用于新观测landsat地表反射率数据上相应的地表单元，进而估算得到landsat尺度fapar数据。转换结果分析：1)与modisfapar数据对比为验证本发明实施例提出的降尺度转换方法的精度，选择两种典型降尺度算法进行对比验证，包括多尺度回归树方法和基于ndvi的线性转换系数法。多尺度回归树方法可直接移植至modisfapar降尺度转换，具体实现过程为：筛选modis尺度纯像元上的高质量反演值作为样本，建立modisfapar与landsat聚合反射率间的线性回归树模型，进而将该模型用于估算landsatfapar值。在样本构建方面，该方法与本发明实施例的降尺度转换方法存在一定差异。该算法不考虑地类和土壤类型，直接对选取的样本，按照反射率进行树划分并对叶结点建立fapar线性回归模型。该算法中的纯像元筛选标准仅考虑cv。为保证算法对比的基准一致，具体实现时，其和本发明实施例的方法都采用了相同的样本，即同时考虑cv和升尺度地类频率提取出的样本。基于ndvi的线性转换系数法的具体实现过程为：计算landsatndvi，并将其聚合至modis尺度；对每个观测日期，逐像元计算聚合ndvi到modisfapar的转换系数；将该转换系数用于将landsatndvi转换至landsatfapar。将三种降尺度转换算法得到的landsatfapar与modisfapar进行对比，分析二者的一致性，作为衡量降尺度算法精度和多尺度自洽性的指标之一。首先，选择植被生长过程中的代表性日期(年积日分别为：112,176,224,272,304)，通过目视判读对比降尺度转换结果与modisfapar数据，如图3～7所示。图3～7中a均为modisfapar，b均为采用多尺度回归树方法的降尺度转换结果，c均为基于ndvi的线性转换系数法的降尺度转换结果，d均为本实施例方法的降尺度结果。由图可知，基于ndvi的线性转换系数法的结果与modisfapar空间分布的整体趋势最为一致，但空间连续性较差，表现出明显的斑块状特征，且存在较多的数据空洞。尤其是在植被生长初期和植被衰退期(如年积日为112时)，研究区西南地形起伏区域的植被分布在空间上呈现出不连续、不规则的特征。这表明，基于ndvi的线性转换系数法尽管能在理论上与modisfapar数据高度吻合，但却丢失了大量的空间细节信息，降尺度转换结果空间连续性较差，且精度直接依赖modisfapar的反演精度。本实施例方法与多尺度回归树方法得出的结果，均表现出较好的空间连续性。但多尺度回归树方法结果在植被生长期(如年积日为176和224时)表现出对fapar的明显低估，而本实施例算法得到的fapar在分布趋势和值域范围上都与modisfapar更为接近。这是由于多尺度回归树方法建立在历史数据的回归分析上，大量数据的平均效应虽稳定、可靠，但如不考虑fapar的年际变化和季候差异，会引入额外误差。进一步对三种算法结果与modisfapar间进行散点图分析。在产品的直接对比与验证中，通常认为取3×3像元窗口可有效减小配准误差的影响。因此，将modisfapar和landsatfapar同时升尺度至1.5km，选取modis高质量反演数据对比二者量上的差异，结果如图8(a)～(c)所示。多尺度回归树方法的结果在fapar高值部分存在明显低估。基于ndvi的线性转换系数法理论上可以得到与modisfapar完全吻合的结果，但散点图显示仍有不少像元存在明显偏离。这是由于其忽略了modis亚像元尺度上的异质性。尤其是当modis亚像元尺度上存在非植被地类时，landsat尺度上ndvi<0，而聚合至modis尺度ndvi接近0，会产生较大的转换系数，导致fapar异常值的产生，这是一种典型的尺度效应。本发明实施例的降尺度转换方法结果与modisfapar整体吻合较好，在fapar高值部分略有低估。以平均绝对误差mae(meanabsoluteerror)和均方根误差rmse(rootmeansquarederror)为指标，定量分析landsatfapar与modisfapar间的差异，如表1所示。表1landsatfapar与modisfapar差值统计qc＝0像元多尺度回归树方法基于ndvi的线性转换系数法本发明实施例的方法mae0.03570.00430.0275rmse0.05860.02410.0454由表1可看出，基于ndvi的线性转换系数法结果与modisfapar相比差异最小，这与其算法原理相关。其次是本发明实施例的降尺度转换方法得出的结果，而多尺度回归树方法结果与modisfapar相比差异最大。综合看来，本发明实施例的降尺度转换方法得出的结果既与modisfapar在空间趋势和量值分布上都具有较好的一致性，又表现出较好的时空连续性，具有明显优势。2)与地面验证数据对比地面测量数据基于独立观测手段，是验证降尺度算法精度最为有效而可靠的指标。本实施例采用的地面验证数据是经过预处理及时间归一化校正的fapar地面测量数据，即校正至卫星过境时刻的fapar测量数据。地面测量数据的可用日期分别为2012年6月8日、6月24日和7月10日。与地面实测数据直接对比验证：首先对三种算法结果与地面测量数据进行散点图分析，如图9(a)～(c)所示。由于modisfapar在农田反演精度较高，以modisfapar为参考的三种算法结果与地面实测数据的吻合都较好。但多尺度回归树方法在fapar高值部分(7月10日)表现出明显低估，基于ndvi的线性转换系数法在fapar高值部分(7月10日)有高估的倾向。本发明实施例的方法结果与地面实测fapar吻合较好，仅在fapar低值部分(6月8日)略有高估。进一步分析三种算法结果与地面实测数据量值上的差异，如表2所示。以mae和rmse指标定量评价。表2landsatfapar与地面实测fapar差值统计评价指标多尺度回归树方法基于ndvi的线性转换系数法本发明实施例的方法mae0.06930.06690.0546rmse0.08940.07980.0710由表2可看出，本发明实施例的方法精度最高，明显优于多尺度回归树方法和基于ndvi的线性转换系数法。与asterfapar真实参考数据对比：首先对比三种算法结果与asterfapar真实参考数据的空间分布，如图10(a)(b)-13(a)(b)所示，多尺度回归树方法结果在6月24日时与asterfapar空间分布吻合较好，但在7月10日植被生长旺盛时表现出明显低估；基于ndvi的线性转换系数法表现出明显的斑块状特征，空间分布连续性较差，尤其是在荒漠与农田交界处空间不连续特征尤为明显；本发明实施例的方法在两个日期都与asterfapar空间分布趋势非常一致，具有很好的空间连续性。进一步由散点图对比三种算法结果与asterfapar真实参考数据量值上的差异。为避免配准误差的影响，将landsatfapar和asterfapar同时升尺度至90m。如图14(a)～(c)所示，多尺度回归树方法在fapar低值部分有所高估，在fapar高值部分有明显低估。基于ndvi的线性转换系数法整体上对fapar存在明显高估。受尺度效应影响，大量像元在尺度转换中产生了异常的转换系数导致fapar值远大于1，表现出fapar值截断至1的特征。本发明实施例的方法在fapar低值部分存在一定程度上的高估，整体上与asterfapar具有高度的一致性。分别将landsatfapar和asterfapar升尺度至90m和480m，分析三种算法结果与地面实测数据量值上的差异，以mae和rmse指标定量评价，结果如表3所示。表3landsatfapar与asterfapar差值统计由表3可知，本发明实施例的方法精度最高，多尺度回归树方法次之，基于ndvi的线性转换系数法精度最低。结合表2，尽管基于ndvi的线性转换系数法与modisfapar高度吻合，但由于未考虑modis亚像元尺度上的异质性，导致其与地面验证数据相比精度并不理想。本发明实施例的方法既与modisfapar吻合较好，与地面验证数据相比精度也较高，且空间分布连续性较好，具有明显优势。3)时间趋势从时间趋势上，通过与modisfapar和地面站点实测数据对比，检验本发明实施例的方法的结果的精度。首先检验降尺度fapar与modisfapar时间趋势上的一致性。选择农田、林地、稀疏植被、裸地等地类的典型像元绘制时间序列曲线，如图15(a)～(d)所示。典型像元是以变异系数cv和升尺度地类频率为指标，筛选出modis各地类具有代表性的纯像元。在纯像元基础上，landsatfapar与modisfapar吻合越好，说明降尺度转换结果与modisfapar的自洽性越好。由图15(a)～(d)可知，在modis尺度上，各地类像元内的异质性随时间推移呈增大趋势。但无论像元内部异质性高或者低，landsatfapar升尺度均值能够与modisfapar在时间趋势上很好地吻合，说明了本发明实施例的方法精度可靠。其次，比较13个地面站点fapar实测数据与对应landsat像元上降尺度fapar的时间趋势。如图16所示，在所有站点上，本发明实施例的方法降尺度转换得到的fapar均能很好地展现出作物的生长变化规律，曲线平滑、连续，且都能与地面站点实测fapar的时间趋势很好地吻合。上述结果均表明：多尺度回归树算法能够保证降尺度转换结果的空间连续性和自洽性，但却未融合植被参数的时变特征，估算得到的landsatfapar结果与modisfapar吻合程度最低，尤其是对植被生长旺盛时期的fapar存在明显低估；基于ndvi的线性转换系数法能够实现实时的降尺度转换，考虑到了植被参数的季候差异，但却是以牺牲modis亚像元尺度上的空间细节信息、忽略尺度效应为代价的。该方法降尺度转换结果的空间连续性较差，存在明显的斑块状特征，与地面真实参考数据相比精度也较低。本发明实施例提出的方法，既对尺度转换过程中的尺度效应进行了充分考虑，又通过融合新观测数据对fapar的时变特征进行了充分考虑，在此基础上实现了实时的降尺度转换。本发明实施例的降尺度转换方法得到的结果精度最高、时空连续性最好，与modisfapar和地面验证数据相比均表现出很好的一致性，且本发明实施例方法降尺度转换得到的fapar，在升尺度至modis尺度时精度(表3：rmse＝0.0264)更优于modisfapar产品([3]：rmse＝0.054)。既验证了发明实施例方法降尺度转换方法的优越性，也说明了构建遥感趋势面的重要意义：提高遥感产品的时空分辨率，有望提高遥感产品的精度。[3]y.wang,d.xie,s.liu,r.hu,y.li,andg.yan,"scalingoffaparfromthefieldtothesatellite,"remotesensing,vol.8,p.310,2016.modisfapar产品生产用到的算法没有显性的表达式，是由复杂的三维辐射传输模型反演而来的，但是它用途广泛、精度较高。因此，本实施例方法重点是建立modisfapar与landsat反射率间的关系，得到landsat尺度上fapar的显性表达式。为了尺度自洽，通过地表单元划分实现了对二者关系的分层线性化。方法与多尺度回归树方法有所类似，即都是建立了多元线性回归模型。但建立方法不同：第一：多尺度回归树方法等是基于样本自然划分，缺少机理，且样本难以代表所有像元，而本实施例的方法是基于图像划分地表单元，建立在fapar估算机理和尺度转换理论基础上的。第二：多尺度回归树方法是基于多期数据建立一个模型后，对任意一天的影像都采用同样的模型。这种方法没有考虑到二者关系的复杂性，即随时间会发生变化，所以在对某一景影像fapar进行估算时会产生误差，并且难以实现fapar的实时估算。本实施例的方法采用了贝叶斯实时融合新的数据，考虑到fapar估算模型随时间的变化。基于ndvi的线性转换系数法可以实时进行尺度转换，但是由于ndvi存在尺度效应，它无法实现自洽的尺度转换。并且，基于ndvi的线性转换系数法忽略了modis亚像元尺度上的异质性，对modis像元内部所有landsat像元采用相同的转换系数，导致了结果误差较大、空间连续性较差。与本实施例的方法区别也较大。本发明实施例基于遥感趋势面的主要思想，在尺度自洽的前提下，发展了一种基于先验知识的、实时融合新观测数据的modisfapar降尺度转换方法，通过对研究区modisfapar降尺度转换的对比验证表明，本发明实施例提出的方法在算法精度、多尺度自洽性和时空连续性等方面，均优于多尺度回归树法和基于ndvi的线性转换系数法等两种典型算法。在先验知识库支持下，利用土壤类型和地表覆盖类型数据对地表单元进行合理划分，将modisfapar与landsat地表反射率间的非线性关系在特定的地表单元转换为线性函数，避免了尺度转换过程中的尺度效应，保证了降尺度转换结果的多尺度自洽性。且在大量历史观测数据的基础上，建立了modisfapar降尺度转换函数作为先验知识，在贝叶斯方法框架下通过融合新观测数据，实现了对landsatfapar的实时估算。以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。当前第1页1 2 3
技术特征：1.基于低分辨率遥感产品降尺度的高分辨率fapar估算方法，其特征在于，按照以下步骤进行：
步骤s1，基于先验知识库，通过地表单元划分对modisfapar反演模型进行分层线性化，建立modisfapar与landsat地表反射率间的尺度转换函数，得到fapar先验估算模型；
步骤s2，基于贝叶斯框架融合modisfapar和landsat地表反射率新观测数据，更新modisfapar与landsat地表反射率间的尺度转换函数，得到更新的fapar后验估算模型；
步骤s3，基于更新后的fapar后验估算模型，基于新观测landsat数据实时估算fapar。
2.根据权利要求1所述的基于低分辨率遥感产品降尺度的高分辨率fapar估算方法，其特征在于，所述步骤s1按照以下步骤进行：
步骤s11，数据收集与预处理；
使用landsat地表反射率数据，影像空间分辨率30m，采用可见光、近红外波段，筛选历史年限和当年云量小于10％的高质量地表反射率数据；历史年限数据用于先验知识构建，当年的数据作为新观测信息；地表反射率数据产品空间基准为wgs-84坐标系和utm投影，对数据进行解压缩、裁剪和数值转换；
使用modisfapar产品，在modis分辨率尺度上，筛选出与landsat影像观测日期对应的modisfapar产品，对产品数据进行数据层提取、裁剪、投影转换、重采样、数值转换；
步骤s12，选择研究区具有代表性的一景影像，基于landsat影像分析，采用k均值聚类算法将研究区的landsat影像自动划分为林地、农田、稀疏植被、裸地和其他五种地类类别，将landsat分类结果按照众数规则升尺度至1km尺度，得到分辨率1km的地表覆盖类型数据；使用土壤类型栅格数据，对土壤类型数据和地表覆盖类型进行空间上的交集操作，将研究区划分为相应的地表单元，每一地表单元具有特定的土壤类型和地表覆盖类型；
步骤s13，筛选modisfapar高质量纯像元样本；
包括(1)，筛选modisfapar高质量反演值：选用modis主算法反演的fapar高质量反演数据建立样本，筛选标准为modisfapar质量控制标识qc＝0；
(2)，modis纯像元筛选：采用升尺度地类频率和变异系数cv来度量modis亚像元尺度上地表异质性，进行纯像元筛选；变异系数定义为标准差与均值的比值；按照cv≤0.2且升尺度地类频率≥0.95的标准筛选纯像元；在纯像元筛选中仅使用植被生长期的数据；
对每景landsat地表反射率数据进行如下处理：
a)，将landsat地表反射率升尺度至modis480m分辨率尺度，并计算变异系数cv：首先，对landsat可见光-近红外波段的每一个波段，计算每一modis像元的波段cv值，即该modis像元空间范围包含的16×16个landsat像元波段反射率的标准差与均值的比值；其次，计算每一个modis像元各个波段cv值的均值，作为该像元的cv值；
b)，将landsat地表分类数据升尺度至480m，按照众数规则得出modis像元尺度上的优势地类，即升尺度像元对应的地类；同时，计算该优势地类在modis亚像元尺度上所占的面积比例，即升尺度地类频率；
(3)，对历史年限的每一景landsat地表反射率数据及modisfapar数据，按照步骤(1)-(2)中的过程，筛选出具有modisfapar高质量反演值的纯像元；然后，根据划分的地表单元，筛选出每一类地表单元对应的modisfapar高质量反演值的纯像元；提取所有筛选像元对应的modisfapar值和聚合至480m分辨率的landsat可见光和近红外波段的地表反射率值，即构成每一地表单元的筛选样本；
步骤s14，对每一地表单元，基于筛选样本，建立modisfapar与landsat地表反射率间的尺度转换函数，通过回归分析建立fapar先验估算模型。
3.根据权利要求2所述的基于低分辨率遥感产品降尺度的高分辨率fapar估算方法，其特征在于，所述步骤s14按照以下步骤进行：
modisfapar值与地表反射率间的关系可表示为：
其中，下标m代表modis尺度；faparm代表modisfapar值；sm代表modis尺度的地表反射率；fm表示modis尺度的fapar估算模型，对于modis产品，fm代表一个复杂的物理模型；下标l代表landsat尺度，i表示给定modis像元内对应landsat像元的索引，sl(i)表示modis像元内第i个landsat像元的地表反射率；n代表modis像元对应landsat像元的总数；
结合土壤类型和地表覆盖类型进行地表单元划分，对每一地表单元建立特定的fapar估算模型，形式为：
其中，下标soil代表土壤类型，下标lc代表地表覆盖类型，表示特定土壤类型和地表覆盖类型条件下modis尺度fapar估算模型；
当为线性模型时，有表示特定土壤类型和地表覆盖类型条件下landsat尺度fapar估算模型；建立在modis尺度上的fapar线性估算模型可适用于landsat尺度fapar的估算，即建立了尺度自洽的fapar估算模型，则landsat尺度fapar的估算faparl可由下式计算：
其中，sj表示第j波段的地表反射率，表示线性函数第j波段的系数，k为波段总数，为线性函数的常数项；
landsat尺度上，对同类地表单元，fapar估算采用相同的线性模型则modis纯像元上的fapar可由下式计算：
式中，si,j为modis像元内对应第i个landsat像元第j波段的地表反射率；为该modis像元内所有landsat像元第j波段地表反射率的均值，也就是landsat像元第j波段地表反射率聚合至modis480m分辨率的均值；由上式可建立modisfapar与landsat地表反射率之间的联系，即在已知faparm和si,j的条件下，求出尺度自洽的fapar估算模型
对每一地表单元，筛选具备modisfapar高质量反演值的modis纯像元，即可在已知modisfapar和landsat多波段的地表反射率聚合值的情况下，通过回归分析，求解该地表单元对应的fapar线性估算模型的系数和得出该地表单元对应的fapar先验估算模型。
4.根据权利要求3所述的基于低分辨率遥感产品降尺度的高分辨率fapar估算方法，其特征在于，所述步骤s14中，可采用基于landsat红光、近红外波段反射率的二元回归模型建立尺度转换函数，也即尺度自洽的fapar估算模型；模型形式表示为：
式中，s代表多波段地表反射率；s3和s4分别为landsat红光、近红外波段的地表反射率；和分别表示线性函数fsoil,lc的常数项、变量s3和s4的系数，即模型系数；fsoil,lc表示尺度自洽的fapar估算模型，即同时适用于modis和landsat两个尺度的fpar估算模型；fapar表示与估算模型运行尺度对应尺度即modis或landsat尺度上的fapar值；
对每个地表单元分别建立modisfapar与landsat地表反射率聚合值的线性模型，即降尺度转换函数，以此来实现对modisfapar反演模型分层线性化；通过对每一地表单元筛选样本进行回归分析，求解上式中的模型系数，即可建立尺度自洽的fapar先验估算模型。
5.根据权利要求1-4任一项所述的基于低分辨率遥感产品降尺度的高分辨率fapar估算方法，其特征在于，所述步骤s2的具体步骤是：
步骤s21，对每一景新观测的landsat地表反射率数据，采用k均值聚类算法，将地表自动划分为林地、农田、稀疏植被、裸地和其他；
选取其对应日期的modisfapar数据；
基于每一对新观测的landsat地表反射率数据和modisfapar数据，对每一个地表单元，按照步骤s13，筛选具有fapar高质量反演值的modis纯像元，提取其对应的modisfapar值与landsat红光、近红外波段地表反射率聚合值，构成新观测数据样本；同时提取新观测modisfapar数据中提供的fapar反演标准差；
步骤s22，将步骤s14中建立的fapar先验估算模型作为先验知识，按照贝叶斯公式融合新观测数据，更新fapar先验估算模型系数，得到fapar后验估算模型。
6.根据权利要求5所述的基于低分辨率遥感产品降尺度的高分辨率fapar估算方法，其特征在于，所述步骤s22的具体步骤是：
将fapar估算模型参数x看成是具有某种先验分布的随机变量，则modis尺度上fapar的估计值fm可由modis尺度反射率sm和fapar估算模型参数x计算得到，有
fm＝smx；
fapar估算模型参数x分布可由概率分布函数描述，也即先验分布p(x)；融合modis尺度上新观测信息fm，由先验分布得到模型参数后验概率的关键是贝叶斯公式，表示为：
其中，p(x)为fapar估算模型参数x的先验概率密度函数，p(fm|x)为给定fapar估算模型参数x时得到modisfapar新观测数据fm的条件概率；p(fm)＝∫m′p(fm|x)p(x)dx，m′为fapar估算模型参数x的取值空间，p(fm)表示得到新观测值fm的全部概率；p(x|fm)为fapar估算模型参数x的后验概率，即给定modisfapar新观测数据fm时得到模型参数x的条件概率，融合了先验知识与新观测信息，是关于x的全部知识；
假设模型参数误差与观测数据误差均服从高斯分布，根据贝叶斯定理，可通过最小化代价函数得到fapar估算模型参数x的最优估计，即fapar估算模型参数x后验分布的均值μx和协方差矩阵cx为：
其中，cd为观测误差协方差矩阵，cp是fapar估算模型参数x先验分布的协方差矩阵，xp为先验知识对fapar估算模型参数x的最优估计，即fapar估算模型参数先验估计值；
使用球形协方差矩阵定义cp和cd，表示为：
其中，代表fapar估算模型参数先验估计值xp的方差，代表modisfapar新观测数据的方差，i(r)为r×r的单位矩阵，i(m)为m×m的单位矩阵，其中，r值为fapar估算模型参数个数，m为新观测值数量，即对每一地表单元，基于新观测数据筛选出具备高质量modisfapar反演值的modis纯像元数量；由modisfapar数据中给出的fapar反演标准差σd计算得出；建立的fapar估算模型为landsat数据红光、近红外波段反射率构成的线性函数，共有3个参数，即r＝3；是由筛选样本以fapar估算模型形式进行线性回归分析中，各模型参数估计的标准差σp计算得出；基于土壤类型和新观测地表覆盖类型划分的地表单元，按照上式，通过融合新观测数据对fpar先验估算模型参数进行更新。
技术总结基于低分辨率遥感产品降尺度的高分辨率FAPAR估算方法，先基于先验知识库，通过地表单元划分对MODIS FAPAR反演模型进行分层线性化，建立MODIS FAPAR与Landsat反射率间的尺度转换函数(即尺度自洽的FAPAR估算模型)；基于贝叶斯框架融合MODIS FAPAR和Landsat反射率新观测数据，更新MODIS FAPAR与Landsat反射率间的尺度转换函数；最后基于新观测Landsat数据实时估算FAPAR。解决了现有估算方法难以同时实现尺度自洽、空间连续、实时的高分辨率FAPAR估算的问题。
技术研发人员：王祎婷;阎广建;谢东辉
受保护的技术使用者：西安科技大学
技术研发日：2020.01.10
技术公布日：2020.06.05
