一种年度区域植被覆盖度计算方法

摘要

本发明公开了一种年度区域植被覆盖度计算方法，包括：对同年度区域所有遥感影像进行精确的辐射定标、大气校正和几何配准，得到遥感影像结果图集合；基于构建的云层指数、阴影指数、背景指数及植被指数，提取相应专题信息；对不同时相遥感影像分别进行去云算法和去阴影算法处理；对单时相植被指数信息进行空间合成；对年度内多时相信息进行融合，计算年度区域植被指数信息空间分布图；最后，计算年度区域植被覆盖度信息空间分布图。本发明消除了天气、地形、时相、植被类型、环境等因素对植被覆盖度的影响，解决了年度区域植被覆盖度计算的问题，提高了区域尺度植被覆盖度计算的准确度和可靠性，操作过程简单、灵活，易于在区域尺度中推广应用。

说明书

技术领域

本发明涉及一种年度区域植被覆盖度计算方法，属于生态建设、生态环境保护、灾害监测、生态功能评估等应用领域。

背景技术

植被覆盖度通常定义为植被(叶、茎、枝)在地面的垂直投影面积占统计区总面积的百分比，是衡量地表植被状况的主要指标之一，也是区域生态系统环境变化的重要指标，在大气圈、土壤圈、水圈和生物圈中占据着重要地位。植被覆盖度被认为是进行土地荒漠化、石漠化、盐渍化等土地退化评估中最有效的指标；它反映了生态系统对土壤侵蚀的影响，是控制土壤侵蚀的积极因素，是开展水土流失评价和综合治理的重要指标，也是水文生态模型中的一个重要变量，进而在区域生态功能评估和资源环境综合承载力评估中具有重要意义；同时，植被覆盖度与植被蒸散有着密切的联系，而植被蒸散是能量平衡与水分平衡的重要组成部分，是土壤-植被-大气系统水热通量传输中的一个过程，蒸散的监测需要根据不同植被覆盖度采用不同的方法，植被覆盖度是作为一个重要的控制因子而存在；另外，植被覆盖度也是全球变化研究中的一个敏感因子，是许多区域乃至全球气候模型、碳过程评估模型中必不可少一个变量，可为生态环境变化、区域植被生产力估算、土地覆盖变化监测提供有用的信息。因此，进行区域植被覆盖度的计算及其计算精度的提高是各个应用领域发展所需要的，对生态建设、环境保护、水文、农业、国土资源利用、灾害监测以及全球变化研究等都具有重要意义。

目前，植被覆盖度的计算方法大体可以分为两种：一种是传统的地面实测法，另一种是遥感监测法。

传统的地面实测法虽然测量精度较高，但该方法费时、费力，加上植被覆盖度本身具有显著的时空分异特征，该方法不适合区域尺度研究，仅适用于较小空间范围的研究应用中。

遥感技术的出现为监测区域植被覆盖度，乃至全球植被覆盖度变化提供了可能。现有遥感监测植被覆盖度的方法有很多，主要包括：经验模型法、混合像元分解法、物理模型法、光谱梯度法。经验模型法是基于实测数据，建立相应的统计模型或经验回归模型，进而计算植被覆盖度；该方法虽然简单、易于计算，但需要大量的实测数据以建立相关模型，同时，针对不同植被类型需要建立不同的经验模型，易受时间、地域条件、植被类型等制约，在大区域研究中具有较大的不确定性和局限性，不易推广应用。混合像元分解模型法也称为亚像元分解法，即在一定假设条件下，分解多个组分构成的像元中的遥感信息，建立像元分解模型，从而获得植被覆盖度。混合像元分解模型主要有线性模型、概率模型和模糊分析模型等，其中线性分解模型应用最为广泛。线性分解模型假设到达传感器的光子只与一个组分作用，不同组分相互独立，通过数值方法求解各组分在像元中的比例即是对应的植被覆盖度，其精度很大程度上取决于各端元的合理选取，因此它不适合低空间分辨率的遥感数据应用。常用的线性分解模型有像元二分法、人工神经网络法、线性光谱分析法、多端元光谱混合分析法等，此类方法要求空间分辨率较高的遥感数据。物理模型法是通过研究光与植被的相互作用以建立植被光谱信息与植被覆盖度物理关系的模型，理论上具有较广泛的应用，但此类模型需要大量的相关参数数据，而现有卫星遥感数据在保证模型精度的情况，很难满足参数要求，且相关模型计算也比较复杂，难易在区域尺度范围上推广。光谱梯度法是在分析植被和土壤反射基础上提出，并假设在有限波段范围内土壤反射率随波长线性变化，没有考虑植被土壤面积随波长的变化，模型误差较大。

随着遥感技术的进步，遥感影像的空间分辨率已不再是精确获取植被覆盖度的限制因子，遥感的重访周期也大大缩短。与此同时，因天气、地形、传感器本身等因素，影像噪声(如云层、阴影、背景信息等)也随之增多，如何有效地降低噪声的影响，同时最大程度地利用遥感影像中的有用信息以提高植被覆盖度计算的准确性，已是当下急需解决的技术问题。另外，现有对植被覆盖度计算研究应用中，通常选择该地区植被生长比较好的一个时相的数据来代表该地区的植被覆盖度。而，现实中一个地区地表植被覆盖类型往往多种多样，每种类型的植被其生长季都有所不同；即使是同一类型的植被，因土壤条件、水分条件、种植时间等差异也会导致其生长状况的差异。因此，需要考虑区域范围内不同植被类型以及生长环境的影响，充分利用多个时相数据，以提高区域植被覆盖度计算的准确性、可靠性。现有应用不仅需要能真实代表该地区的植被覆盖度信息，尤其是随着研究范围区域的扩大，上述问题对植被覆盖度准确计算的影响会越来越显著。

发明内容

针对现有技术存在的不足，本发明提出了一种年度区域植被覆盖度计算方法，该方法充分考虑了云层、阴影、背景等噪声的干扰，能最大限度地挖掘和利用遥感影像信息，提高了现有植被覆盖度计算的准确度，操作过程简单、灵活，方便在较大的区域尺度上推广应用。

为解决上述问题，本发明所采取的技术方案是：

一种年度区域植被覆盖度计算方法，包括以下步骤：

S1，对同一年度植被生长季内区域所有遥感影像进行精确的辐射定标和大气辐射校正，得到不同时相代表地物真实反射率的遥感影像结果图；

S2，针对步骤S1所得遥感影像结果图，以任一时相遥感影像结果图为基准，对其他遥感影像结果图进行精确的几何配准，形成一个在地理空间坐标上相匹配的遥感影像结果图集合；

S3，针对步骤S2所得遥感影像结果图集合中每一时相遥感影像结果图分别进行以下处理：

S3-1，针对该时相遥感影像结果图分别计算云层指数、植被指数和阴影指数，得到相应指数信息空间分布图；

S3-2，针对该时相遥感影像结果图计算背景指数，并进行二值化处理，得到一个影像背景信息空间分布图；

S3-3，针对步骤S3-1所得植被指数信息空间分布图，进行去云算法处理，该算法包括以下子步骤：

S3-3-1，对该时相遥感影像结果图进行影像的增强算法处理，并在经过增强算法处理后的遥感影像结果图上随机采集云层区的像元样本为样本区；

S3-3-2，将该样本区与步骤S3-1所得云层指数信息空间分布图进行空间叠加，统计计算样本区对应云层指数信息数值从小到大的累计频率；

S3-3-3，根据该累计频率，设置一阈值范围，并制作两个与云层指数信息空间分布图相匹配的二值信息图；

S3-3-4，将制作的两个二值信息图与步骤S3-1所得植被指数信息空间分布图和步骤S3-2所得影像背景信息空间分布图进行空间运算，即可得到消除云层噪声后的植被指数信息结果图；

S3-4，针对步骤S3-1所得植被指数信息空间分布图，进行去阴影算法处理，该算法包括以下子步骤：

S3-4-1，对该时相遥感影像结果图进行影像的增强算法处理，并在经过增强算法处理后的遥感影像结果图上随机采集阴影区的像元样本为样本区；

S3-4-2，将该样本区与步骤S3-1所得阴影指数信息空间分布图进行空间叠加，统计计算样本区对应阴影指数信息数值从小到大的累计频率和相关统计变量；

S3-4-3，根据该累计频率及相关统计变量，设置一阈值范围，并制作两个与阴影指数信息空间分布图相匹配的二值信息图；

S3-4-4，将制作的两个二值信息图与步骤S3-1所得植被指数信息数据空间分布图和步骤S3-2所得影像背景信息空间分布图进行空间运算，即可得到消除阴影噪声后的植被指数信息结果图；

S3-5，基于步骤S3-2所制作的影像背景信息空间分布图以及步骤S3-3和步骤S3-4所制作的二值信息图，制作两个新的二值信息图，进而利用新的二值信息图，对步骤S3-3和步骤S3-4所得植被指数信息结果图进行单时相信息的空间合成，生成该时相的最终植被指数信息结果图；

S4，基于步骤S3所得各时相的最终植被指数信息结果图，将区域内年度所有时相的最终植被指数信息结果图进行融合，形成一个包含多个波段的影像数据；

S5，在步骤S4基础之上，对覆盖相同区域的最终植被指数信息结果图，采用最大值合成法，进行最大值合成处理，得到该年度区域植被指数信息空间分布图；

S6，对步骤S5所得区域植被指数信息空间分布图进行数值运算处理，消除云层、阴影和背景信息对统计分析结果的影响；为下一步进行统计分析提供准确可靠的数据源；

S7，基于步骤S6所得结果，计算植被覆盖度，即得该年度该区域准确的植被覆盖度空间分布图。

进一步的，所述步骤S3-1中构建的云层指数CCI，其计算公式为：CCI＝(ρ_R+ρ_G+ρ_B)*ρ_NIR，ρ_R为红波段的光谱反射率值，ρ_G为绿波段的光谱反射率值，ρ_B为蓝波段的光谱反射率值，ρ_NIR为遥感影像近红外波段光谱反射率值。

进一步的，所述步骤S3-1中植被指数为归一化植被指数NDVI，该指数计算公式为：ρ_NIR为遥感影像近红外波段光谱反射率值，ρ_R为红波段的光谱反射率值。

进一步的，所述步骤S3-1中所述阴影指数YYI计算公式为：其中ρ_NIR为遥感影像近红外波段光谱反射率值，ρ_SWIR为短波红外波段的光谱反射率值,NDVI为归一化植被指数。

进一步的，所述步骤S3-2中构建的背景指数BGI，其计算公式为：其中ρ_NIR为遥感影像近红外波段光谱反射率值,ρ_R为红波段的光谱反射率值，ρ_G为绿波段的光谱反射率值。

进一步的，所述步骤S3-2中二值化处理是指将背景指数BGI中背景信息像元赋值为-99，其他像元赋值为0，得到一个影像背景信息空间分布图BGII。

进一步的，所述步骤S3-3-1中影像的增强算法处理是指2％的线性拉伸增强；随机采集样本像元是指采用2×2的采样窗口随机采集云层信息像元，为避免错选上非云层区像元而造成误差，要求样本像元数量大于200。

进一步的，所述步骤S3-3-3中阈值范围的设置是指将累计频率达到2％时的云层指数信息数值作为阈值的下限，该阈值范围上限阈值为无穷大。

进一步的，所述步骤S3-3-3中两个二值信息图为CCIB1和CCIB2，CCIB1是针对云层指数CCI，将阈值范围内的数值赋值为0，其他值赋值为1，计算而得；CCIB2是针对云层指数CCI，将阈值范围内的数值赋值为-99，其他值赋值为0，计算而得。

进一步的，所述步骤S3-3-4中空间运算的计算公式为：NDVI_CCI＝NDVI*CCIB1+CCIB2+BGII，其中NDVI_CCI为消除云层噪声后的植被指数信息结果图，BGII是影像背景信息空间分布图，CCIB1和CCIB2是步骤S3-3-3计算所得二值信息图。

进一步的，所述步骤S3-4-1中所述影像的增强算法处理是指2％的线性拉伸增强；随机采集样本像元是指采用2×2的采样窗口随机采集影像阴影信息像元，为避免错选上非云层区像元而造成误差，要求样本像元数量大于200。

进一步的，所述步骤S3-4-3中所述阈值范围的设置是指将累计频率达到95％时的阴影指数信息数值作为阈值的上限，最小值统计变量(样本区对应阴影指数数值)作为阈值范围的下限。

进一步的，所述步骤S3-4-3中所述两个二值信息图为YYIB1和YYIB2，YYIB1是针对阴影指数YYI，将阈值范围内的数值赋值为0，其他值赋值为1，计算而得；YYIB2是针对阴影指数YYI，将阈值范围内的数值赋值为-99，其他值赋值为0，计算而得。

进一步的，所述步骤S3-4-4中空间运算的计算公式为：NDVI_YYI＝NDVI*YYIB1+YYIB2+BGII，其中NDVI_YYI为消除阴影信息后的植被指数信息结果图，BGII是步骤S3-2计算所得影像背景信息空间分布图，YYIB1和YYIB2是步骤S3-4-3计算所得二值信息图。

进一步的，所述步骤S3-5中两个新的二值信息图为NDVI_B1和NDVI_B2，其计算公式分别如下：

NDVI_B2＝(NDVI_B1+99)*(NDVI_B1-1)

其中BGII是步骤S3-2计算所得背景信息空间分布图，CCIB1是步骤S3-3计算所得二值信息图，YYIB1是步骤S3-4计算所得二值信息图。

进一步的，所述步骤S3-5中所述空间合成方法，其计算公式如下：

其中NDVI_CCI为步骤S3-3所得消除云层噪声后的植被指数信息结果图，NDVI_YYI为步骤S3-4所得消除阴影信息后的植被指数信息结果图，NDVI_B1和NDVI_B2为步骤S3-5中两个新的二值信息图。

进一步的，步骤S4中所述融合方法是指将多个单时相植被指数信息通过波段合成得到一个在空间上相匹配的影像。

进一步的，所述步骤S5中最大值合成法是指针对同一个像元，该像元在不同时相会有不同的数值，取这些值中的最大值，即为最大值合成。

进一步的，步骤S6中所述数值运算处理是指将步骤S5所得区域植被指数信息空间分布图中数值为-99的区域赋值为NAN，为下一步进行统计分析提供准确可靠的数据源。

进一步的，所述步骤S7中计算植被覆盖度可以利用任何一种基于植被指数信息来估算植被覆盖度的方法或模型。

进一步的，所述任何一种估算植被覆盖度的方法或模型包括混合像元分解法中常用的像元二分法、神经网络模型。

进一步的，所述像元二分法计算公式如下：

其中VFC为像元内的植被覆盖度，NDVI为该像元上的归一化植被指数，NDVI_soil为无植被覆盖像元的NDVI值，NDVI_veg为完全被植被覆盖的像元NDVI值；在本发明中，以0.5％置信度截取区域NDVI的上下阈值分别代表NDVI_veg和NDVI_soil，进而计算区域植被覆盖度VFC。

本发明的有益效果在于：

1)本发明的优点是从分析天气、地形、时相特征等条件对区域植被覆盖度计算的影响入手，构建了云层指数、阴影指数、背景指数等模型算法，最大程度地消除了环境背景噪声对植被覆盖度计算结果的影响，解决了年度区域植被覆盖度计算的问题，提高了区域尺度植被覆盖度计算的准确度和可靠性，操作过程简单、灵活，易于在区域尺度应用中推广。

2)相比现有技术方法，本发明充分利用了现今遥感数据的高时间分辨率、多时相特征，综合考虑区域内植被类型、生长环境、生长季的差异，提出了一种年度植被覆盖度计算方法，提高了植被覆盖度作为衡量地表植被状况、区域生态系统环境变化主要指标的准确度、可靠性和可信度，以及作为生态、水文、气候模型必要输入参数的准确度，其具有较强的普适性和适应性，易推广应用。

附图说明

图1为本发明的流程示意图；

图2为四个不同时相的遥感影像图；

图3为去云和去阴影的空间运算原理示意图；

图4为四个时相经过去云和去阴影处理后的植被指数信息结果图；

图5为年度区域植被指数信息空间分布图；

图6为年度区域准确的植被覆盖度空间分布图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本实施例以某区域同年度、四个不同时相的遥感影像为例进行说明。

参照图1，图1为本发明一种年度区域植被覆盖度计算方法的流程示意图，包括以下步骤：

S1，收集本区域同一年度生长季内所有遥感影像，分别对其进行精确的辐射定标和大气辐射校正，得到能代表地物真实反射率的遥感影像结果图。如图2为本区域四个时相的遥感影像图，图2中第一幅图(第一行左图)为6月份的遥感影像图，图2中第二幅图(第一行右图)为7月份的遥感影像图，图2中第三幅图(第二行左图)为8月份的遥感影像图，图2中第四幅图(第二行右图)为9月份的遥感影像图，除了季节性的差异外，同时都有不同程度的云层、云层阴影、山体阴影、背景噪声等信息。要获取该年度本区域的植被覆盖度，任一时相遥感数据均不能完全代表该区域的植被状况。

S2，基于步骤一所得结果，以能覆盖本区域的任一时相遥感影像结果图为基准，对其他三个时相的遥感影像进行精确的几何配准，使这四个时相的遥感影像数据在地理空间坐标上相匹配，得到一个标准的遥感影像结果图集合。

S3，针对步骤S2所得遥感影像结果图集合中每一时相遥感影像结果图分别进行以下处理：

S3-1，针对该时相遥感影像结果图分别计算云层指数、植被指数和阴影指数，得到相应指数信息空间分布图；

由于卫星遥感传感器受大气密度、云层变化等影响，同一区域很多影像都存在云层遮挡问题。云层会使获得的遥感地物信息衰减，甚至损失。在本实施例中云层遮挡会降低或掩盖植被指数信息，致使最终所计算植被覆盖度的准确度、可信度大大下降。本发明通过深入云层信息在卫星遥感传感器不同光谱波段的响应特征，以及与植被、土壤、水体、裸地、城镇等在相同波段位置的响应异同，选择卫星遥感传感器常用的、且共有的光谱波段，构建了一个云层指数CCI。基于该云层指数可以有效地提取云层信息，进而消除云层信息对植被覆盖度计算的影响。该云层指数CCI计算公式为：CCI＝(ρ_R+ρ_G+ρ_B)*ρ_NIR，ρ_R为红波段的光谱反射率值，ρ_G为绿波段的光谱反射率值，ρ_B为蓝波段的光谱反射率值，ρ_NIR为遥感影像近红外波段光谱反射率值。

目前，归一化植被指数NDVI被广泛应用于遥感监测植被生长状况中，同时也是本实施例计算植被覆盖度的基础。另外，在对不同地物光谱响应特性分析时，归一化植被指数NDVI的引入能有效地把云层阴影和山体阴影与遥感影像中水体信息等低亮度值的地物区分开来。因此，需首先对归一化植被指数NDVI进行计算，计算方法主要是利用非线性归一化的方法对较为简单的比值植被指数进行处理，其计算公式为：

ρ_NIR为遥感影像近红外波段光谱反射率值，ρ_R为红波段的光谱反射率值。

云层除了对地物信息有遮挡和削弱作用外，同时，云层本身阴影对阴影区的地物植被指数信息也有一定程度的降低作用。如果能将阴影区的地物信息被另一时相无阴影的遥感影像信息所替代，就能较好地提高最终计算植被覆盖度的准确度和可信度。除了云层会有阴影存在，因太阳高度角、方位角以及地形因素，山体阴影也是不容忽略的一部分。在分析不同地物光谱特性时，云层阴影和山体阴影的反射率值很低，如果将近红外波段光谱反射率值(ρ_NIR)与短波红外波段的光谱反射率值(ρ_SWIR)进行加和，能很好地拉大阴影信息与其他地物信息间的距离，但水体信息与阴影信息间距离较近，容易造成错判，两者的和不能将阴影与水体较好的分开；通过引入归一化植被指数(NDVI)则能进一步显著地拉大水体信息与阴影信息间的距离，进而构建了一个阴影指数(YYI)，其计算公式为：

其中，YYI为阴影指数，ρ_NIR为遥感影像近红外波段光谱反射率值，ρ_SWIR为短波红外波段的光谱反射率值,NDVI为归一化植被指数。

S3-2，针对该时相遥感影像结果图计算背景指数，并进行二值化处理，得到一个影像背景信息空间分布图；

在实际工作中，有很多遥感影像并不能完全覆盖整个地区，未覆盖的区域会有背景信息存在(见图2所示)，如何处理这些信息将直接影响整个技术方法的实施效率；另外，背景信息在经过大气校正后被默认赋值为0，这会造成跟其他地物的植被指数信息数值一样，处理不当将会影响最终植被覆盖度计算的准确度。本发明通过分析不同地物在各个波段的光谱响应特征及数量级，构建了背景指数BGI，其计算公式为：

其中ρ_NIR为遥感影像近红外波段光谱反射率值,ρ_R为红波段的光谱反射率值，ρ_G为绿波段的光谱反射率值。该指数只能有效地把背景信息与影像中地物信息区精确地区分开来，还需要做进一步的处理。本技术方案，对其做进一步的二值化处理，将背景指数BGI中背景信息像元赋值为-99，其他像元赋值为0，得到一个影像背景信息空间分布图BGII，以便作为一个信息因子参与下一步计算。

S3-3，针对步骤S3-1所得植被指数信息空间分布图，进行去云算法处理，该算法包括以下子步骤：

S3-3-1，在处理之前，首先对该时相遥感影像结果图进行影像的增强算法处理，并在经过增强算法处理后的遥感影像结果图上随机采集云层区的像元样本为样本区，所述影像的增强算法处理是指2％的线性增强算法处理，以增加不同地物类型色彩对比度，方便进行下一步像元样本的选择，同时，选择2％线性拉伸算法也能降低环境背景噪声的影响，保真度较高。为避免错选非云层覆盖区像元，在经过增强算法处理后的影像图上随机采集云层区的像元采用2×2的采样窗口随机采集云层信息像元样本为样本区，为保证数据精度，样本像元数量大于200。

S3-3-2，将该样本区与步骤S3-1所得云层指数信息空间分布图进行空间叠加，统计计算样本区对应云层指数信息数值从小到大的累计频率；

S3-3-3，根据该累计频率，设置一阈值范围，并制作两个与云层指数信息空间分布图相匹配的二值信息图；所述阈值范围的下限为累计频率达到2％时的云层指数信息数值，所述阈值范围的上限为无穷大。所述两个与云层指数信息空间分布图相匹配的二值信息图分别命名为CCIB1和CCIB2，CCIB1是针对步骤三中云层指数CCI，将阈值范围内的数值赋值为0，其他值赋值为1，计算而得，见图3中②；而CCIB2是针对云层指数CCI，将阈值范围内的数值赋值为-99，其他值赋值为0，计算而得，见图3中③。

S3-3-4，将制作的两个二值信息图与步骤S3-1所得植被指数信息空间分布图和步骤S3-2所得影像背景信息空间分布图进行空间运算，即可得到消除云层噪声后的植被指数信息结果图；该空间运算的计算公式为：

NDVI_CCI＝NDVI*CCIB1+CCIB2+BGII，其中NDVI_CCI为消除云层噪声后的植被指数信息结果图，BGII是影像背景信息空间分布图，CCIB1和CCIB2是步骤S3-3-3计算所得二值信息图。为更清楚地解释该空间运算的原理，有助于对本发明计算方案的理解，如图3为该空间运算过程原理的示意图。其中①为NDVI，②为CCIB1，③为CCIB2，④为BGII，⑤为NDVI_CCI。

S3-4，针对步骤S3-1所得植被指数信息空间分布图，进行去阴影算法处理，该算法包括以下子步骤：

S3-4-1，对该时相遥感影像结果图进行影像的增强算法处理，并在经过增强算法处理后的遥感影像结果图上随机采集阴影区的像元样本为样本区；

在处理之前，首先对该时相遥感影像结果图进行影像的增强算法处理，并在经过增强算法处理后的遥感影像结果图上随机采集阴影区的像元样本为样本区，所述影像的增强算法处理是指2％的线性增强算法处理，以增加不同地物类型色彩对比度，方便进行下一步像元样本的选择，同时，选择2％线性拉伸算法也能降低环境背景噪声的影响，保真度较高。为避免错选非阴影区像元，在经过增强算法处理后的影像图上随机采集阴影区的像元采用2×2的采样窗口随机采集阴影信息像元样本为样本区，为保证数据精度，样本像元数量大于200。

S3-4-2，将该样本区与步骤S3-1所得阴影指数信息空间分布图进行空间叠加，统计计算样本区对应阴影指数信息数值从小到大的累计频率和相关统计变量(样本区对应阴影指数数值)；

S3-4-3，根据该累计频率及相关统计变量，设置一阈值范围，将累计频率达到95％时的阴影指数信息数值作为阈值的上限，最小值统计变量作为阈值范围的下限，并制作两个与阴影指数信息空间分布图相匹配的二值信息图，分别命名为YYIB1和YYIB2，YYIB1是针对阴影指数YYI，将阈值范围内的数值赋值为0，其他值赋值为1，计算而得；YYIB2是针对阴影指数YYI，将阈值范围内的数值赋值为-99，其他值赋值为0，计算而得。

S3-4-4，将制作的两个二值信息图与步骤S3-1所得植被指数信息数据空间分布图和步骤S3-2所得影像背景信息空间分布图进行空间运算，即可得到消除阴影噪声后的植被指数信息结果图。该空间运算的计算公式为：

NDVI_YYI＝NDVI*YYIB1+YYIB2+BGII，其中NDVI_YYI为消除阴影信息后的植被指数信息结果图，BGII是步骤S3-2计算所得影像背景信息空间分布图，YYIB1和YYIB2是步骤S3-4-3计算所得二值信息图。该空间运算的原理与去云处理算法中的一样，可参见图3所示。

S3-5，基于步骤S3-2所制作的影像背景信息空间分布图以及步骤S3-3和步骤S3-4所制作的二值信息图，制作两个新的二值信息图，分别命名为NDVI_B1和NDVI_B2，其计算公式分别如下：

NDVI_B2＝(NDVI_B1+99)*(NDVI_B1-1)

其中BGII是步骤S3-2计算所得背景信息空间分布图，CCIB1是步骤S3-3计算所得二值信息图，YYIB1是步骤S3-4计算所得二值信息图。

进而利用新的二值信息图，对步骤S3-3所得消除云层噪声后的植被指数信息结果图和步骤S3-4所得消除阴影噪声后的植被指数信息结果图进行单时相信息的空间合成，生成该时相的最终植被指数信息结果图。该空间合成方法计算公式如下：

其中NDVI_CCI为步骤S3-3所得消除云层噪声后的植被指数信息结果图，NDVI_YYI为步骤S3-4所得消除阴影信息后的植被指数信息结果图，NDVI_B1和NDVI_B2为步骤S3-5中两个新的二值信息图。

如图4为四个时相经过去云和去阴影处理后的最终植被指数信息结果图，其中图4中第一幅图(第一行左图)为6月份的最终植被指数信息结果图，图4中第二幅图(第一行右图)为7月份的最终植被指数信息结果图，图4中第三幅图(第二行左图)为8月份的最终植被指数信息结果图，图4中第四幅图(第二行右图)为9月份的最终植被指数信息结果图。

S4，基于步骤S3所得4个时相的最终植被指数信息结果图，将区域内年度4个时相的最终植被指数信息结果图进行融合，形成一个包含多个波段的影像数据。在本实施例中是利用ENVI中的LayerStacking将步骤S3所得四个时相的最终植被指数信息结果图合成在一个空间上相匹配的影像。

S5，在步骤S4基础之上，对覆盖相同区域的四个时相的最终植被指数信息结果图，采用最大值合成法，进行最大值合成处理，得到该年度区域植被指数信息空间分布图，如图5为最大值合成结果图。

S6，对步骤S5所得区域植被指数信息空间分布图进行数值运算处理，消除云层、阴影和背景信息对统计分析结果的影响，为下一步进行统计分析提供准确可靠的数据源。具体的为在上述运算过程中云层信息、阴影信息和背景信息都被赋值为-99，为避免这些像元信息参与统计计算，避免统计误差，在计算植被覆盖度之前需将数值为-99的区域赋值为NAN。

S7，基于步骤S6所得结果，利用基于植被指数的像元二分法计算本区域植被覆盖度，即得该年度该区域准确的植被覆盖度空间分布图，如图6所示。所述像元二分法计算公式如下：

其中VFC为像元内的植被覆盖度，NDVI为该像元上的归一化植被指数，NDVI_soil为无植被覆盖像元的NDVI值，NDVI_veg为完全被植被覆盖的像元NDVI值；在本发明中，以0.5％置信度截取区域NDVI的上下阈值分别代表NDVI_veg和NDVI_soil，进而计算区域植被覆盖度VFC。

尽管这里参照本发明的多个解释性实施例对本发明进行了描述，但是，应该理解，本领域技术人员可以设计出很多其他的修改和实施方式，这些修改和实施方式将落在本申请公开的原则范围和精神之内。更具体地说，在本申请公开、附图和权利要求的范围内，可以对主题组合布局的组成部件和/或布局进行多种变型和改进。除了对组成部件和/或布局进行的变形和改进外，对于本领域技术人员来说，其他的用途也将是明显的。