[发明专利]结合矩阵补全和趋势滤波的植被指数时间序列重建方法

权利要求书

1.结合矩阵补全和趋势滤波的植被指数时间序列重建方法，其特征在于，包括：

S1：输入植被指数时间序列原始数据y以及与之对应的质量标记数据f，质量标记数据用以表征植被指数时间序列原始数据中每个像素的质量；

S2：针对每个像素的植被指数序列，通过矩阵变化从一维向量变成二维矩阵Y，其中，两个维度分别代表时间序列的年份覆盖长度和一年内的监测频次，用以表征植被的季节性变化特征及年际周期变化相似性，构建的二维矩阵为NDVI时间序列数据矩阵；然后基于质量标记数据，找出二维矩阵中的缺失数据，将对应位置记录为缺失；

S3：针对每个像素变换后的二维矩阵，建立低秩矩阵补全的最优化能量方程，通过非精确增广拉格朗日算法实现矩阵补全，得到初步不含数据缺失的时间序列补全矩阵X；

S4：将补全矩阵X向量化，得到不含缺失的NDVI时间序列数据x；

S5.基于不含缺失的NDVI时间序列数据和质量标记数据f，建立加权趋势滤波的优化方程，通过交替方向乘子法实现噪声滤除，得到平滑干净的高质量NDVI时间序列数据z，作为重建后的NDVI时间序列数据。

2.如权利要求1所述的植被指数时间序列重建方法，其特征在于，步骤S1中，输入的植被指数时间序列原始数据y有n年，每年具有k个时间序列点，y和f的总长度为kn。

3.如权利要求1所述的植被指数时间序列重建方法，其特征在于，当植被指数时间序列原始数据y为MODIS植被指数产品MOD13中的数据时，NDVI数据分为干净像元、可能干净的像元、云覆盖像元以及雪覆盖像元四类，质量标记分别为0,1,2,3。

4.如权利要求3所述的植被指数时间序列重建方法，其特征在于，步骤S2包括：

S2.1：将植被指数时间序列原始数据矩阵化，形成k行n列的矩阵其中每列为一年的NDVI时间序列数据，包括k个时间序列点数据，共有n列；

S2.2：根据质量标记数据f，找出其中被云覆盖以及雪覆盖的像元，并将找出的像元标记为缺失像元，矩阵Y中的对应元素记录为缺失元素，将其他像元标记为干净像元，矩阵Y中的对应元素记录为干净元素。

5.如权利要求1所述的植被指数时间序列重建方法，其特征在于，步骤S3包括：

S3.1：针对每个像素变换后的二维矩阵，建立低秩矩阵补全的最优化能量方程：

式中，为待求的补全矩阵，i表示矩阵的第i行，j表示第j列，用下标(i,j)表示元素的位置，将Y中干净元素的位置(i,j)的集合表示为Ω；

S3.2：采用矩阵的核范数最小化代替矩阵的秩最小化，将式(1)转化成凸优化问题，即：

式中，||X||_*为矩阵的核范数；

S3.3：利用非精确增广拉格朗日算法求解，首先通过引入辅助变量E，将式(2)转化成等价形式：

再构建式(3)的增广拉格朗日函数，即：

式中，M为拉格朗日乘子，μ为惩罚参数，式(4)通过下述交替迭代步骤获得最优解：

其中，上标k表示迭代次数，ρ表示一个大于1的常数，用于更新惩参数μ，P_Ω(·)为矩阵投影算子，为奇异值阈值算子，通过上述迭代求解获得的最优解X,为NDVI时间序列数据补全矩阵。

6.如权利要求1所述的植被指数时间序列重建方法，其特征在于，步骤S4包括：

将NDVI时间序列数据补全矩阵转变成时间序列向量NDVI时间序列数据补全矩阵中不含缺失数据。

7.如权利要求1所述的植被指数时间序列重建方法，其特征在于，步骤S5包括：

S5.1：根据质量标记数据f，将干净像元的权重设置为1，将缺失重建像元的权重设置为0.8，建立加权趋势滤波的优化方程为：

式中，为待求的去噪结果，λ为正则化参数，W为对角矩阵，对角元素为对应像元的权重，D为二阶差分矩阵；

S5.2：利用交替方向乘子法求解，首先将式(6)转化成带约束的模型：

再构建式(7)的增广拉格朗日函数，即：

式中，m为拉格朗日乘子，β为惩罚参数，式(6)通过下述交替迭代步骤获得最优解，

其中，θ表示一个大于1的常数，用于更新惩参数β，

为软阈值函数，上标k表示迭代次数通过上述迭代求解获得的最优解z，为最终的NDVI时间序列数据重建结果。