一种AVHRR地表反射率重建方法、系统与装置

摘要

本发明提供一种AVHRR地表反射率重建方法、系统与装置，所述方法包括：对原始地表反射率数据进行去除无效值处理；对处理获取的有效原始地表反射率数据进行指定时间分辨率的合成处理；基于合成处理获取的合成地表反射率数据，计算归一化植被指数NDVI，并基于所述NDVI，利用给定算法重建NDVI上包络线；基于所述NDVI和所述NDVI上包络线，对所述合成地表反射率数据进行云检测，去除受云影响的数据；基于去除受云影响数据所得的剩余数据，以及所述NDVI上包络线，通过函数拟合，获取指定时间内任意时刻的计算反射率数据，重建时间连续的地表反射率。本发明能够在不影响通用性的基础上有效去除地表反射率重建中的云干扰并进行缺失值填充。

说明书

技术领域

本发明涉及数据处理技术领域，更具体地，涉及一种AVHRR地表反射率重建方法、系统与装置。

背景技术

搭载在NOAA卫星上的AVHRR传感器每天获取全球影像，提供大量的监测大气、海洋、植被和陆表特性的光谱信息，形成一个每天分辨率的长时间全球影像序列。

基于AVHRR数据，利用不同数据处理方法，目前已生成了多个长时间序列数据集，这些数据集具有不同的时空分辨率和时间跨度。其中，最具代表性的一套长时间序列数据集是由LTDR项目生产的每天分辨率的全球NDVI和地表反射率。但是，LTDR AVHRR NDVI和地表反射率产品中包含大量的云干扰信息，而且一些时间的NDVI和地表反射率数据完全缺失。地表反射率和NDVI产品中云等的污染严重制约了这些产品在陆表监测中的应用，而且会导致遥感参数高级产品的时空不一致性。

针对上述问题，目前已有多种重建NDVI时间序列曲线的方法，总的来说这些方法可以分为时间域和频率域方法两大类。其中Tang等提出的一种TSCD算法，用于MODIS地表反射率数据的云检测和缺失值填充。TSCD算法在地表稳定或变化缓慢时，能够取得很好的效果，但该方法依赖于蓝光波段的反射率数据以及其他相关辅助信息。因此，TSCD算法很难直接用于AVHRR数据的云检测和数据重建。

发明内容

为了克服上述问题或者至少部分地解决上述问题，本发明提供一种AVHRR地表反射率重建方法、系统与装置，以达到有效去除AVHRR地表反射率重建中的云干扰并进行缺失值填充的目的。

第一方面，本发明提供一种AVHRR地表反射率重建方法，包括：S1，基于设定规则，对原始地表反射率数据进行去除无效值处理，获取有效原始地表反射率数据；S2，对所述有效原始地表反射率数据进行指定时间分辨率的合成处理，获取合成地表反射率数据；S3，基于所述合成地表反射率数据，计算归一化植被指数NDVI，并基于所述NDVI，利用给定算法重建NDVI上包络线；S4，基于所述NDVI和所述NDVI上包络线，利用给定判定条件，对所述合成地表反射率数据进行云检测，去除受云影响的数据，获取符合质量等级要求的地表反射率数据；S5，基于所述符合质量等级要求的地表反射率数据，以及所述NDVI上包络线，通过函数拟合，获取指定时间内任意时刻的计算反射率数据，重建时间连续的地表反射率。

其中，步骤S1中所述的设定规则包括：若任一像元点的红光波段反射率大于该像元点的近红外波段反射率，则判定该像元点的地表反射率为无效值；和/或，若任一像元点的红光波段与近红外波段增强植被指数大于NDVI，则判定该像元点的地表反射率为无效值。

其中，所述S2的步骤进一步包括：若在所述指定时间分辨率的合成窗口内，存在至少两个所述有效原始地表反射率数据，则采用观测角度约束的NDVI最大值合成法，合成所述有效原始地表反射率数据；或者，若在所述指定时间分辨率的合成窗口内，存在一个所述有效原始地表反射率数据，则将其作为所述合成地表反射率数据；或者，若在所述指定时间分辨率的合成窗口内，不存在所述有效原始地表反射率数据，则以指定年限内所述有效原始地表反射率数据的均值作为所述合成地表反射率数据。

其中，所述采用观测角度约束的NDVI最大值合成法，合成所述有效原始地表反射率数据的步骤进一步包括：在所述指定时间分辨率的合成窗口内，对所述有效原始地表反射率数据的观测天顶角按从小到大的顺序进行排列；基于所述顺序中最小的两个所述观测天顶角对应的所述有效原始地表反射率数据，分别计算对应的NDVI；选取二者中较大的所述对应的NDVI对应的所述有效原始地表反射率数据，作为所述合成地表反射率数据。

其中，步骤S3中所述的基于所述NDVI，利用给定算法重建NDVI上包络线的步骤进一步包括：利用基于三维离散余弦变换的惩罚最小二乘回归法，通过迭代运算，最小化给定的包含时间序列NDVI值向量的第一代价函数，获取最优时间序列NDVI值向量估计；基于所述最优时间序列NDVI值向量估计，构建所述NDVI上包络线。

其中，所述S4的步骤进一步包括：基于所述合成地表反射率数据对应的NDVI，依次判断任一时刻NDVI是否满足如下给定判定条件：

|NDVI_i-NDVI_Env_i|＞α×NDVI_Env_i；

式中，NDVI_i表示第i时刻的NDVI值，NDVI_Env_i表示NDVI上包络线第i时刻的NDVI值，α表示设定阈值；

若判断获知第i时刻NDVI满足所述给定判断条件，则判定所述第i时刻NDVI对应的合成地表反射率数据为受云影响的数据；去除所述合成地表反射率数据中的所述受云影响的数据，获取所述符合质量等级要求的地表反射率数据。

其中，所述S5的步骤进一步包括：基于所述符合质量等级要求的地表反射率数据，获取时间序列红外及近红外波段反射率点集；基于所述红外及近红外波段反射率点集，以及所述NDVI上包络线，利用二次多项式函数，通过迭代运算最小化第二代价函数，拟合给定时间窗口内所述符合质量等级要求的地表反射率数据，获取指定时间内任意时刻的计算反射率数据；基于任意时刻的所述计算反射率数据，重建所述指定时间内时间连续的地表反射率。

其中，所述指定时间分辨率为8天分辨率。

第二方面，本发明提供一种AVHRR地表反射率重建系统，包括：数据预处理模块，用于基于设定规则，对原始地表反射率数据进行去除无效值处理，获取有效原始地表反射率数据；有效数据合成模块，用于对所述有效原始地表反射率数据进行指定时间分辨率的合成处理，获取合成地表反射率数据；NDVI上包络线重建模块，用于基于所述合成地表反射率数据，计算归一化植被指数NDVI，并基于所述NDVI，利用给定算法重建NDVI上包络线；云检测模块，用于基于所述NDVI和所述NDVI上包络线，利用给定判定条件，对所述合成地表反射率数据进行云检测，去除受云影响的数据，获取符合质量等级要求的地表反射率数据；地表反射率重建模块，用于基于所述符合质量等级要求的地表反射率数据，以及所述NDVI上包络线，通过函数拟合，获取指定时间内任意时刻的计算反射率数据，重建时间连续的地表反射率。

第三方面，本发明提供一种AVHRR地表反射率重建装置，包括：至少一个存储器、至少一个处理器、通信接口和总线；所述存储器、所述处理器和所述通信接口通过所述总线完成相互间的通信，所述通信接口用于所述AVHRR地表反射率重建装置与原始地表反射率数据存储设备之间的信息传输；所述存储器中存储有可在所述处理器上运行的计算机程序，所述处理器执行所述程序时实现如如上所述的地表反射率重建方法。

本发明提供的一种AVHRR地表反射率重建方法、系统与装置，在对每天时间分辨率的AVHRR地表反射率数据进行聚合的基础上，重建时间序列上连续平滑的NDVI上包络线，并基于时间序列的NDVI及其上包络线检测受云影响的反射率，最后以时间序列上连续平滑的NDVI上包络线为约束，进行基于高质量的地表反射率时间连续的地表反射率的重建，能够有效去除地表反射率重建中的云干扰并能有效进行缺失值填充。同时由于本发明不依赖于任何辅助信息，同样能够有效用于MODIS、VIIRS、FY、Landsat等不同传感数据的重建，通用性强，具有广阔的应用前景。

附图说明

图1为本发明实施例一种AVHRR地表反射率重建方法的流程图；

图2为本发明实施例一种采用观测角度约束的NDVI最大值合成法合成有效原始地表反射率数据的处理过程流程图；

图3为本发明实施例一种重建NDVI上包络线的处理过程流程图；

图4为本发明实施例一种云检测的处理过程流程图；

图5为本发明实施例一种基于符合质量等级要求的地表反射率数据和NDVI上包络线重建地表反射率的处理过程流程图；

图6为本发明实施例一种AVHRR地表反射率重建系统的结构示意图；

图7为本发明实施例一种AVHRR地表反射率重建装置的结构框图；

图8为根据本发明实施例重建AVHRR NDVI和地表反射率的仿真试验流程图；

图9(a)为根据本发明实施例重建的1992年至2005年Zhangbei地区植被类型站点LTDR AVHRR和GLASS AVHRR的NDVI及地表反射率时间序列仿真曲线以及SG AVHRR NDVI的时间序列仿真曲线；

图9(b)为根据本发明实施例重建的1992年至2005年Yucheng地区植被类型站点LTDR AVHRR和GLASS AVHRR的NDVI及地表反射率时间序列仿真曲线以及SG AVHRR NDVI的时间序列仿真曲线；

图9(c)为根据本发明实施例重建的1992年至2005年Puechabon地区植被类型站点LTDR AVHRR和GLASS AVHRR的NDVI及地表反射率时间序列仿真曲线以及SG AVHRR NDVI的时间序列仿真曲线；

图9(d)为根据本发明实施例重建的1992年至2005年Larose地区植被类型站点LTDR AVHRR和GLASS AVHRR的NDVI及地表反射率时间序列仿真曲线以及SG AVHRR NDVI的时间序列仿真曲线；

图9(e)为根据本发明实施例重建的1992年至2005年Wankama地区植被类型站点LTDR AVHRR和GLASS AVHRR的NDVI及地表反射率时间序列仿真曲线以及SG AVHRR NDVI的时间序列仿真曲线；

图9(f)为根据本发明实施例重建的1992年至2005年Turco地区植被类型站点LTDRAVHRR和GLASS AVHRR的NDVI及地表反射率时间序列仿真曲线以及SG AVHRR NDVI的时间序列仿真曲线；

图10为根据本发明实施例构建的2010年1月9日和7月12日地表反射率的全球仿真分布图。

具体实施方式

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

作为本发明实施例的一个方面，本实施例提供一种AVHRR地表反射率重建方法，参考图1，为本发明实施例一种AVHRR地表反射率重建方法的流程图，包括：

S1，基于设定规则，对原始地表反射率数据进行去除无效值处理，获取有效原始地表反射率数据。

可以理解为，运行中的AVHRR传感器对每天的地表影像数据进行采集和存储，基于这些影像数据，利用一定的处理方法生成的长时间序列地表反射率数据集即为原始地表反射率数据。由于存在干扰等各种因素，这些原始地表反射率数据中会存在一些无效值。考虑到各种干扰因素，设定对应的无效值判断规则，根据该设定规则筛选出原始地表反射率数据中的无效值，其余的原始地表反射率数据中即为有效原始地表反射率数据。

其中可选的，步骤S1中所述的设定规则包括：若任一像元点的红光波段反射率大于该像元点的近红外波段反射率，则判定该像元点的地表反射率为无效值；和/或，若任一像元点的红光波段与近红外波段增强植被指数大于NDVI，则判定该像元点的地表反射率为无效值。

可以理解为，对每天时间分辨率的原始地表反射率数据进行筛选，确定有效原始地表反射率数据，筛选规则包括：对于每一个像元点，若其红光波段反射率大于近红外波段反射率，则判定该像元点的地表反射率为无效值；另外，若一个像元点的红光波段与近红外波段增强植被指数均大于NDVI，则判定该像元点的地表反射率也为无效值。

S2，对所述有效原始地表反射率数据进行指定时间分辨率的合成处理，获取合成地表反射率数据。

可以理解为，在上述步骤去除原始地表反射率数据中的无效值时，并不能消除有效原始地表反射率数据中的云干扰。为了尽可能减小云和大气对AVHRR地表反射率的影响，将每天时间分辨率的AVHRR地表反射率进行合成。为了达到相应的分辨率需求，设定数据合成的时间分辨率，并按照该时间分辨率对有效原始地表反射率数据进行聚合，合成后的数据即为合成地表反射率数据。在一个实施例中，所述指定时间分辨率为8天分辨率。

S3，基于所述合成地表反射率数据，计算归一化植被指数NDVI，并基于所述NDVI，利用给定算法重建NDVI上包络线。

可以理解为，合成的地表反射率仍然包含大量云的信息。因此，根据合成地表反射率数据计算得到的NDVI时间序列上不连续。考虑到受云影响的地表反射率计算得到的NDVI值偏低，采用相应方法筛选受云影响的地表反射率。具体根据上述步骤获取的合成的地表反射率数据，计算数据对应的归一化植被指数NDVI。为了方便筛选出其中受云影响较小或未受云影响的数据，即NDVI值较大的数据，根据计算获取的NDVI，利用一定的数据重建方法重建NDVI上包络线，如基于三维离散余弦变换的惩罚最小二乘回归法。

S4，基于所述NDVI和所述NDVI上包络线，利用给定判定条件，对所述合成地表反射率数据进行云检测，去除受云影响的数据，获取符合质量等级要求的地表反射率数据。

可以理解为，云等对地表反射率的影响一般会导致NDVI值的负偏置噪声，本步骤中采用时间序列的NDVI及其上包络线对受残余云影响的反射率进行检测。通过一一检测合成地表反射率数据中任意时刻对应的NDVI和NDVI上包络线是否满足给定判断条件，判断该时刻的合成地表反射率数据是否受到云的影响。然后去除合成地表反射率数据中通过该方法筛选出的受到云影响的数据，获取到的剩余数据即为符合质量等级要求的地表反射率数据。

S5，基于所述符合质量等级要求的地表反射率数据，以及所述NDVI上包络线，通过函数拟合，获取指定时间内任意时刻的计算反射率数据，重建时间连续的地表反射率。

可以理解为，在根据上述步骤去除了无效数据和受云影响的数据之后，根据获取的符合质量等级要求的地表反射率数据和NDVI上包络线数据，利用一定的数据拟合方法，拟合指定时间段内对应的符合质量等级要求的地表反射率数据。然后通过迭代运算最小化第二代价函数，计算获取该指定时间段内任意时刻的反射率数据，即计算反射率数据。最后根据该指定时间段内所有时刻的计算反射率数据重建时间连续的地表反射率。

本发明实施例提供的一种AVHRR地表反射率重建方法，根据合成地表反射率数据的NDVI重建NDVI上包络线，并基于NDVI和NDVI上包络线对合成地表反射率数据进行云检测，有效剔除受云影响的数据，最终基于剩余未受云影响的数据和NDVI上包络线，实现时间连续的地表反射率的有效重建，能够有效去除地表反射率重建中的云干扰并能有效进行缺失值填充；同时由于本发明不依赖于任何辅助信息，同样能够有效用于MODIS、VIIRS、FY、Landsat等不同传感数据的重建，通用性强，具有广阔的应用前景。

在另一个实施例中，所述S2的步骤进一步包括：若在所述指定时间分辨率的合成窗口内，存在至少两个所述有效原始地表反射率数据，则采用观测角度约束的NDVI最大值合成法，合成所述有效原始地表反射率数据；或者，若在所述指定时间分辨率的合成窗口内，存在一个所述有效原始地表反射率数据，则将其作为所述合成地表反射率数据；或者，若在所述指定时间分辨率的合成窗口内，不存在所述有效原始地表反射率数据，则以指定年限内所述有效原始地表反射率数据的均值作为所述合成地表反射率数据。

可以理解为，如上述实施例所述，为了达到相应的分辨率需求，设定数据合成的时间分辨率，数据合成在该时间分辨率的合成窗口内进行。为方便理解，以8天时间分辨率为指定时间分辨率举例进行说明。由于在前序步骤中对原始地表反射率数据进行了去除无效值处理，因此在8天的时间窗口内，剩余的有效原始地表反射率数据个数存在三种情况，即0个、1个和至少2个。

对于8天时间窗口内剩余的有效原始地表反射率数据为至少2个的情况，采用观测角度约束的NDVI最大值合成法，对这至少2个有效原始地表反射率数据进行合成，合成结果作为本合成窗口内的合成地表反射率数据；对于8天时间窗口内剩余的有效原始地表反射率数据为1个的情况，直接将该有效原始地表反射率数据作为本合成窗口内的合成地表反射率数据；对于8天时间窗口内剩余的有效原始地表反射率数据为0个的情况，计算指定年限内有效原始地表反射率数据的平均值，并将该平均值作为本合成窗口内的合成地表反射率数据。

其中可选的，所述采用观测角度约束的NDVI最大值合成法，合成所述有效原始地表反射率数据的进一步处理步骤参考图2，为本发明实施例一种采用观测角度约束的NDVI最大值合成法合成有效原始地表反射率数据的处理过程流程图，包括：

S21，在所述指定时间分辨率的合成窗口内，对所述有效原始地表反射率数据的观测天顶角按从小到大的顺序进行排列。

仍以8天时间分辨率为指定时间分辨率举例进行说明，可以理解为，对于8天时间窗口内的至少两个有效原始地表反射率数据，获取每个有效数据的观测天顶角，并根据观测天顶角从小到大的顺序，对对应的有效原始地表反射率数据进行排序。

S22，基于所述顺序中最小的两个所述观测天顶角对应的所述有效原始地表反射率数据，分别计算对应的NDVI。

可以理解为，根据上述步骤的排序，选取其中最小的两个观测天顶角对应的有效原始地表反射率数据，分别计算对应的NDVI。

S23，选取二者中较大的所述对应的NDVI对应的所述有效原始地表反射率数据，作为所述合成地表反射率数据。

可以理解为，根据上述步骤计算获取的两个最小观测天顶角分别对应的两个NDVI值，选取NDVI值较大者对应的有效原始地表反射率数据，将该数据作为本合成窗口内的合成地表反射率数据。

本发明实施例提供的一种AVHRR地表反射率重建方法，能够按需进行指定时间分辨率的地表反射率数据合成，从而能够有效减小云和大气对AVHRR地表反射率的影响。

在上述实施例的基础上，其中步骤S3中所述的基于所述NDVI，利用给定算法重建NDVI上包络线的进一步处理步骤参考图3，为本发明实施例一种重建NDVI上包络线的处理过程流程图，包括：

S31，利用基于三维离散余弦变换的惩罚最小二乘回归法，通过迭代运算，最小化给定的包含时间序列NDVI值向量的第一代价函数，获取最优时间序列NDVI值向量估计。

可以理解为，令向量X中包含时间序列的NDVI值x_i,i＝1,2,…,n，W为对角矩阵，包含对应NDVI值x_i的权重w_i。基于三维离散余弦变换的惩罚最小二乘回归法，通过最小化如下第一代价函数，获取X的最优平滑估计

式中，X表示时间序列NDVI值向量，W表示时间序列NDVI值的权重矩阵，表示X的最有估计，D表示拉普拉斯算子，s表示标量常量，能够决定的平滑程度。

通过迭代运算，最小化上述第一代价函数，获取第一代价函数取最小值时的估计值利用Ⅱ类离散余弦变换，第k+1次迭代得到的可表示为如下形式：

式中，和分别表示第k次和第k+1次迭代运算获取的X的估计值，X表示时间序列NDVI值向量，W表示时间序列NDVI值的权重矩阵，Γ为对角矩阵。

上式中Γ的各分量可通过如下公式计算：

Γ_i,i＝(1+s(2-2cos((i-1)π/n))²)^-1；

式中，Γ_i，i表示对角矩阵Γ的第i行第i列元素，s表示标量常量，能够决定的平滑程度，i表示对角矩阵Γ的第i行或第i列，n表示对角矩阵Γ的行数或列数。

由于云等因素的影响，导致NDVI值偏低。因此，在迭代过程中这些偏低的NDVI值被赋予了较小的权重，而高质量的NDVI值则被赋予了较大的权重。权重矩阵W的各元素计算公式如下：

式中，w_i表示第i个NDVI值x_i的对应权重，u_i表示学生化残差。

上式中学生化残差u_i根据如下公式计算：

式中，u_i表示学生化残差，表示第i个观测残差，s表示标量常量，能够决定的平滑程度，MAD表示中位数绝对偏差。

S32，基于所述最优时间序列NDVI值向量估计，构建所述NDVI上包络线。

可以理解为，在上述步骤获取了时间序列NDVI值向量X的最优估计之后，根据该最优估计向量的元素，构建NDVI上包络线。

本发明实施例提供的一种AVHRR地表反射率重建方法，利用离散余弦变换法，能够实现NDVI上包络线的快速重建，而且整个重建过程全自动，尤其适合长时间序列的全球数据的处理。

在上述实施例的基础上，其中所述S4的进一步处理步骤参考图4，为本发明实施例一种云检测的处理过程流程图，包括：

S41，基于所述合成地表反射率数据对应的NDVI，依次判断任一时刻NDVI是否满足如下给定判定条件：

|NDVI_i-NDVI_Env_i|＞α×NDVI_Env_i；

式中，NDVI_i表示第i时刻的NDVI值，NDVI_Env_i表示NDVI上包络线第i时刻的NDVI值，α表示设定阈值。

可以理解为，如上述实施例所述，云干扰会影响表反射率数据的NDV的取值，在根据上述实施例获取合成地表反射率数据的NDVI和NDVI上包络线之后，根据NDVI和NDVI上包络线筛选受云影响的合成地表反射率数据。具体为，首先一一判断给定时间段内任意时刻NDVI是否满足如下给定判定条件：

|NDVI_i-NDVI_Env_i|＞α×NDVI_Env_i；

式中，NDVI_i表示第i时刻的NDVI值，NDVI_Env_i表示NDVI上包络线第i时刻的NDVI值，α表示设定阈值。

S42，若判断获知第i时刻NDVI满足所述给定判断条件，则判定所述第i时刻NDVI对应的合成地表反射率数据为受云影响的数据。

可以理解为，对于第i时刻的NDVI值NDVI_i，如果满足上述步骤中给定判定条件，则认为该时刻的地表反射率受云等的影响。否则，认为该时刻的地表反射率的质量好，不受云等的影响。

S43，去除所述合成地表反射率数据中的所述受云影响的数据，获取所述符合质量等级要求的地表反射率数据。

可以理解为，在根据上述步骤筛选出受云影响的数据之后，将这些数据从合成地表反射率数据中剔除，其余数据即为符合质量等级要求的地表反射率数据。

本发明实施例提供的一种AVHRR地表反射率重建方法，根据NDVI和NDVI上包络线，对地表反射率数据中的残余云进行检测，剔除所有受云影响的反射率，只保留质量好的地表反射率数据，以便融合NDVI的上包络线，重建时间序列上连续的地表反射率。

在上述实施例的基础上，其中所述S5的进一步处理步骤参考图5，为本发明实施例一种基于符合质量等级要求的地表反射率数据和NDVI上包络线重建地表反射率的处理过程流程图，包括：

S51，基于所述符合质量等级要求的地表反射率数据，获取时间序列红外及近红外波段反射率点集。

可以理解为，令(t_i,I_i),i＝1,2,…,m为根据上实施例获取的符合质量等级要求的地表反射率数据的时间序列点集，其中t_i为时间，I_i为不同波段的反射率。本实例以红光和近红外波段反射率为例进行处理，因此其中和分别表示t_i时刻红光和近红外波段反射率。

S52，基于所述红外及近红外波段反射率点集，以及所述NDVI上包络线，利用二次多项式函数，通过迭代运算最小化第二代价函数，拟合给定时间窗口内所述符合质量等级要求的地表反射率数据，获取指定时间内任意时刻的计算反射率数据。

可以理解为，在根据上述步骤获取红光和近红外波段反射率之后，对每一数据点，利用二次多项式函数f(t)＝at²+bt+c拟合合成窗口内2n+1个地表反射率数据，即利用迭代算法最小化如下第二代价函数：

式中，J(X_i)表示代价函数，X_i＝(a_red,b_red,c_red,a_NIR,b_NIR,c_NIR)^T表示二次多项式函数的系数，f_red(t_i)和f_NIR(t_i)分别表示红光和近红外波段的二次多项式函数，和分别表示t_i时刻红光和近红外波段反射率，NDVI_Sim_i表示利用二次多项式函数拟合的反射率数据计算获取的NDVI，NDVI_Env_i表示NDVI上包络线第i时刻的NDVI值，N表示合成窗口的大小。

在通过上述迭代运算获取上述拟合二次多项式函数的系数之后，利用二次多项式函数f_red(t_i)和f_NIR(t_i)分别计算t_i时刻的红光和近红外波段反射率，即获取计算反射率数据。

S53，基于任意时刻的所述计算反射率数据，重建所述指定时间内时间连续的地表反射率。

可以理解为，在基于上述步骤获取任意时刻的红光和近红外波段反射率之后，利用这两个波段反射率分别进行对应波段的时间连续的地表反射率的重建。

本发明实施例提供的一种AVHRR地表反射率重建方法，利用符合质量等级要求的地表反射率数据和NDVI上包络线，通过二次多项式函数拟合合成窗口内地表反射率数据，实现红光和近红外波段时间连续的地表反射率的重建，无需设置任何参数，重建地表反射率的过程稳定、快速且自动化程度高。

作为本发明实施例的另一个方面，本实施例提供一种AVHRR地表反射率重建系统，参考图6，为本发明实施例一种AVHRR地表反射率重建系统的结构示意图，包括：数据预处理模块1、有效数据合成模块2、NDVI上包络线重建模块3、云检测模块4和地表反射率重建模块5。

其中，数据预处理模块1用于基于设定规则，对原始地表反射率数据进行去除无效值处理，获取有效原始地表反射率数据；有效数据合成模块2用于对所述有效原始地表反射率数据进行指定时间分辨率的合成处理，获取合成地表反射率数据；NDVI上包络线重建模块3用于基于所述合成地表反射率数据，计算归一化植被指数NDVI，并基于所述NDVI，利用给定算法重建NDVI上包络线；云检测模块4用于基于所述NDVI和所述NDVI上包络线，利用给定判定条件，对所述合成地表反射率数据进行云检测，去除受云影响的数据，获取符合质量等级要求的地表反射率数据；地表反射率重建模块5用于基于所述符合质量等级要求的地表反射率数据，以及所述NDVI上包络线，通过函数拟合，获取指定时间内任意时刻计算反射率数据，重建时间连续的地表反射率。

可以理解为，组成本实施例地表反射率重建系统的处理程序模块至少包括数据预处理模块1、有效数据合成模块2、NDVI上包络线重建模块3、云检测模块4和地表反射率重建模块5，各程序模块通过相互间的数据传递和程序模块的计算处理，实现地表反射率的重建。具体体现为：

由于存在干扰等各种因素，原始地表反射率数据中会存在一些无效值。考虑到各种干扰因素，由数据预处理模块1设定对应的无效值判断规则，并根据该设定规则筛选出原始地表反射率数据中的无效值，获取其余的原始地表反射率数据中即为有效原始地表反射率数据。

由于数据预处理模块1并不能消除有效原始地表反射率数据中的云干扰。为了尽可能减小云和大气对AVHRR地表反射率的影响，有效数据合成模块2将每天时间分辨率的AVHRR地表反射率进行合成。为了达到相应的分辨率需求，有效数据合成模块2设定数据合成的时间分辨率，并按照该时间分辨率对有效原始地表反射率数据进行聚合加成，合成后的数据即为合成地表反射率数据。

由于有效数据合成模块2合成的地表反射率仍然包含大量云的信息，根据合成地表反射率数据计算得到的NDVI时间序列上不连续。考虑到受云影响的地表反射率计算得到的NDVI值偏低，NDVI上包络线重建模块3采用相应方法筛选受云影响的地表反射率。具体根据有效数据合成模块2获取的合成地表反射率数据，NDVI上包络线重建模块3计算数据对应的归一化植被指数NDVI。为了方便筛选出其中受云影响较小或未受云影响的数据，根据计算获取的NDVI，利用一定的数据重建方法重建NDVI上包络线，如基于三维离散余弦变换的惩罚最小二乘回归法。

云检测模块4采用时间序列的NDVI及其上包络线对受残余云影响的反射率进行检测。通过一一检测合成地表反射率数据中任意时刻对应的NDVI和NDVI上包络线是否满足给定判断条件，云检测模块4判断该时刻的合成地表反射率数据是否受到云的影响。然后云检测模块4去除合成地表反射率数据中通过该方法筛选出的受到云影响的数据，获取到的剩余数据即为符合质量等级要求的地表反射率数据。

在上述程序模块去除了无效数据和受云影响的数据之后，地表反射率重建模块5根据获取的符合质量等级要求的地表反射率数据和NDVI上包络线数据，利用一定的数据拟合方法，拟合指定时间段内符合质量等级要求的地表反射率数据。即通过迭代运算最小化第二代价函数，计算获取该指定时间段内任意时刻的反射率数据，即计算反射率数据。最后地表反射率重建模块5根据该指定时间段内所有时刻的计算反射率数据重建时间连续的地表反射率。

本发明实施例提供的一种AVHRR地表反射率重建系统，能够有效去除地表反射率重建中的云干扰并能有效进行缺失值填充，同时由于本发明不依赖于任何辅助信息，同样能够有效用于MODIS、VIIRS、FY、Landsat等不同传感数据的重建，通用性强，具有广阔的应用前景。

作为本发明实施例的又一个方面，本实施例提供一种AVHRR地表反射率重建装置，参考图7，为本发明实施例一种AVHRR地表反射率重建装置的结构框图，包括：至少一个存储器701、至少一个处理器702、通信接口703和总线704。

其中，存储器701、处理器702和通信接口703通过总线704完成相互间的通信，通信接口703用于所述地表反射率重建装置与原始地表反射率数据存储设备之间的信息传输；存储器701中存储有可在处理器702上运行的计算机程序，处理器702执行所述程序时实现如上述实施例所述的地表反射率重建方法。

可以理解为，所述的地表反射率重建装置中至少包含存储器701、处理器702、通信接口703和总线704，且存储器701、处理器702和通信接口703通过总线704形成相互之间的通信连接，并可完成相互间的通信。

通信接口703实现地表反射率重建装置与原始地表反射率数据存储设备之间的通信连接，并可完成相互间信息传输，如通过通信接口703实现原始地表反射率数据的读取等。

装置运行时，处理器702调用存储器701中的程序指令，以执行上述各方法实施例所提供的方法，例如包括：基于所述NDVI和所述NDVI上包络线，利用给定判定条件，对所述合成地表反射率数据进行云检测，去除受云影响的数据，获取符合质量等级要求的地表反射率数据等。

本发明另一个实施例中，提供一种非暂态计算机可读存储介质，所述非暂态计算机可读存储介质存储计算机指令，所述计算机指令使所述计算机执行如上述实施例所述的地表反射率重建方法。

可以理解为，实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。

以上所描述的地表反射率重建方法的实施例仅仅是示意性的，其中作为分离部件说明的单元可以是或者也可以不是物理上分开的，既可以位于一个地方，或者也可以分布到不同网络单元上。可以根据实际需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性的劳动的情况下，即可以理解并实施。

通过以上实施方式的描述，本领域的技术人员可以清楚地了解，各实施方式可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件。基于这样的理解，上述技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在计算机可读存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令，用以使得一台计算机设备(如个人计算机，服务器，或者网络设备等)执行上述各方法实施例或者方法实施例的某些部分所述的方法。

本发明实施例提供的一种AVHRR地表反射率重建装置和一种非暂态计算机可读存储介质，能够有效去除地表反射率重建中的云干扰并能有效进行缺失值填充，同时由于本发明不依赖于任何辅助信息，同样能够有效用于MODIS、VIIRS、FY、Landsat等不同传感数据的重建，通用性强，具有广阔的应用前景。

为了更清楚的说明本发明的技术方案和有益效果，本发明实施例进行了仿真试验，仿真过程和仿真结果如下：

参考图8，为根据本发明实施例重建AVHRR NDVI和地表反射率的仿真试验流程图，图中在对每天时间分辨率的AVHRR地表反射率数据进行8天时间分辨率聚合的基础上，重建时间序列上连续平滑的NDVI上包络线，并基于该NDVI上包络线及时间序列的NDVI检测受云影响的反射率数据，最后以时间序列上连续平滑的NDVI上包络线为约束，进行基于高质量的地表反射率时间连续的地表反射率的重建。

仿真试验中根据本发明实施例提供的方法重建1982年至2015年的全球NDVI和地表反射率。为方便起见，将重建的NDVI和地表反射率记为GLASS AVHRR，将AVHRR原始反射率数据记为LTDR AVHRR。在不同植被类型站点上，将GLASS AVHRR NDVI和地表反射率与LTDRAVHRR NDVI和地表反射率进行比较，同时考虑到使用最为广泛的SG滤波方法，将GLASSAVHRR NDVI与利用SG滤波重建的NDVI(记为SG AVHRR)进行比较。另外，在空间上将GLASSAVHRR地表反射率与LTDR AVHRR地表反射率进行比较。仿真试验中，对SG滤波设置大小为5的平滑窗口，图9和图10为本发明实施例的仿真结果示意图。

参考图9，为根据本发明实施例重建的1992年至2005年不同植被类型站点LTDRAVHRR和GLASS AVHRR的NDVI及地表反射率时间序列仿真曲线以及SG AVHRR NDVI的时间序列仿真曲线。图9(a)至图9(f)分别为6个不同类型站点1992年至2005年的GLASS AVHRRNDVI和地表反射率、LTDR AVHRR NDVI和地表反射率以及SG AVHRR NDVI的时间序列仿真曲线，其中曲线1表示GLASS AVHRR，曲线2表示LTDR AVHRR，曲线3表示SG AVHRR。

由图9可见，在2000年和1994年下半年，LTDR AVHRR地表反射率数据缺失，而根据本发明实施例提供的重建方法，在这些站点LTDR AVHRR地表反射率缺失的时间段，都提供了合理的NDVI和地表反射率数据。在这些站点，大多数LTDR AVHRR红光波段反射率值都在0.2以上，对应的NDVI值都接近0，这主要是由于云等污染导致。根据本发明实施例提供的重建方法成功的剔除了这些受云影响的地表反射率，并重建了时间序列上连续的NDVI和地表反射率数据。

同时，重建的GLASS AVHRR NDVI与LTDR AVHRR NDVI的上包络线具有非常好的一致性，而重建的GLASS AVHRR地表反射率与LTDR AVHRR地表反射率的下包络线具有很好的一致性；除Yucheng站点外，GLASS AVHRRNDVI和SG AVHRR NDVI的时间序列曲线具有非常好的一致性。SG滤波在Zhangbei,Puechabon,Larose,Wankama和Turco站点能够取得很好的NDVI重建效果。但是由于Yucheng站点的植被类型为具有两个生长季节的作物，很显然SG滤波的窗口设置过大，导致在该站点SG AVHRR NDVI不能很好的重建两个生长季节之间较小的NDVI值。

不同于SG滤波法需要确定合适的平滑窗口，本发明提供的方法、系统和设备，无需设置任何参数，因此重建地表反射率的过程稳定、快速且全自动，非常适合长时间序列的全球数据的处理。同时该方法不依赖于任何辅助信息，可以用于MODIS、VIIRS、FY、Landsat等不同传感数据的重建，具有广阔的应用前景。另外，本发明提供的方法、系统和设备，可以同时重建不同波段的地表反射率数据，能够有效避免不同波段数据分别重建时导致的物理意义上的不一致性。

参考图10，为根据本发明实施例构建的2010年1月9日和7月12日地表反射率的全球仿真分布图。其中，图10(a)和图10(b)分别为2010年1月9日的GLASS AVHRR和LDTR AVHRR地表反射率RGB影像，图10(c)和图10(d)分别为2010年7月12日的GLASS AVHRR和LDTRAVHRR地表反射率RGB影像。

由图10(b)可见，由于LDTR AVHRR地表反射率被云污染，图中大部分区域都被云覆盖。10(a)中应用本发明实施例的重建方法，有效剔除了这些残余云，得到空间上完整的GLASS AVHRR地表反射率。同样的，由于LDTR AVHRR地表反射率被云污染，10(d)中存在大片的云，10(c)中应用本发明实施例的重建方法，有效剔除了所有受云影响的反射率值，并重建了合理的地表反射率值。

综上，本发明提供的一种AVHRR地表反射率重建方法、系统与装置，在对每天时间分辨率的AVHRR地表反射率数据进行聚合的基础上，重建时间序列上连续平滑的NDVI上包络线，并基于时间序列的NDVI及其上包络线检测受云影响的反射率，最后以时间序列上连续平滑的NDVI上包络线为约束，进行基于高质量的地表反射率时间连续的地表反射率的重建，能够有效去除地表反射率重建中的云干扰并能有效进行缺失值填充。同时由于本发明不依赖于任何辅助信息，同样能够有效用于MODIS、VIIRS、FY、Landsat等不同传感数据的重建，通用性强，具有广阔的应用前景。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。