CN114463202A - 结合矩阵补全和趋势滤波的植被指数时间序列重建方法 - Google Patents
结合矩阵补全和趋势滤波的植被指数时间序列重建方法 Download PDFInfo
- Publication number
- CN114463202A CN114463202A CN202210037286.5A CN202210037286A CN114463202A CN 114463202 A CN114463202 A CN 114463202A CN 202210037286 A CN202210037286 A CN 202210037286A CN 114463202 A CN114463202 A CN 114463202A
- Authority
- CN
- China
- Prior art keywords
- matrix
- data
- vegetation index
- time sequence
- completion
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 239000011159 matrix material Substances 0.000 title claims abstract description 125
- 238000000034 method Methods 0.000 title claims abstract description 50
- 238000001914 filtration Methods 0.000 title claims abstract description 36
- 230000008859 change Effects 0.000 claims abstract description 22
- 230000003190 augmentative effect Effects 0.000 claims abstract description 16
- 238000005457 optimization Methods 0.000 claims abstract description 16
- 239000013598 vector Substances 0.000 claims abstract description 12
- 238000012217 deletion Methods 0.000 claims description 7
- 230000037430 deletion Effects 0.000 claims description 7
- 230000001932 seasonal effect Effects 0.000 claims description 7
- 230000009466 transformation Effects 0.000 claims description 7
- 150000001875 compounds Chemical class 0.000 claims description 6
- 238000012544 monitoring process Methods 0.000 claims description 6
- 230000000737 periodic effect Effects 0.000 claims description 4
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000006870 function Effects 0.000 description 11
- 230000008569 process Effects 0.000 description 4
- 238000012545 processing Methods 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 239000000443 aerosol Substances 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 238000004177 carbon cycle Methods 0.000 description 1
- 230000019771 cognition Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000010485 coping Effects 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000001308 synthesis method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
- G06T2207/30188—Vegetation; Agriculture
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
本发明提供了一种结合矩阵补全和趋势滤波的植被指数时间序列重建方法，包括如下步骤：首先通过植被指数产品中的质量标记数据确定时间序列数据中的缺失位置，然后针对单个像素的植被指数序列，通过矩阵变化从一维向量变成二维矩阵，接着针对每个像素变换后的矩阵，建立低秩矩阵补全的最优化能量方程，通过非精确增广拉格朗日算法实现矩阵补全，得到初步不含数据缺失的时间序列补全矩阵。最后再将该补全矩阵进行向量化，在一维向量的基础上建立加权趋势滤波的能量优化方程，通过交替方向乘子法实现模型的求解，从而进一步滤除残留的噪声，得到平滑干净的高质量植被指数时间序列数据，实现长时间遥感植被指数序列的高精度重建。
Description
技术领域
本发明涉及遥感技术和植被生态技术领域，尤其涉及一种结合矩阵补全和趋势滤波的植被指数时间序列重建方法。
背景技术
遥感植被指数是研究陆地生态***植被状态和变化的重要指标，能够反应大尺度区域植被结构和功能的长时间动态变化。在过去的几十年里，随着大量的光学遥感卫星和传感器成像仪的升空，如Landsat、AVHRR、MODIS、SPOT等，形成了系列具备长达十余年甚至几十年的遥感植被指数时间序列，包括归一化植被指数NDVI、增强植被指数EVI等。目前这些标准的NDVI时间序列产品，在植被动态监测、农业估产、土地覆盖变化、陆地碳循环领域得到了广泛应用，对于生态环境建设、全球变化应对与认知起到了重要支撑。然而，光学传感器成像过程中，不可避免的会受到传感器故障、云、气溶胶、及大气条件等的影响，从而在植被指数时间序列中存在较多的噪声和信息缺失。虽然目前MODIS、AVHRR、SPOT等这些标准植被指数产品都会使用最大值合成方法来减少噪声和缺失的影响，但仍会存在着较多的偏差，极大的影响了其后续应用。
为此，国内外学者发展了系列时域滤波方法来去除遥感植被指数时间序列中的噪声和缺失问题，主要可以分为以下几类：首先是基于函数拟合方法，通过预先设定的函数形式，对整个植被指数时间序列进行拟合，从而去除噪声的影响，典型的方法如非对称高斯函数(AG)、双逻辑曲线(DL)等方法。这类方法在表征植被生长的季节性规律变化模式方面具有独特的优势，从而被广泛应用于物候信息的提取，但在季节性变化不规律以及复杂的植被区域应用能力有限。第二个类别是基于局部滑动窗口的滤波方法，通过局部开窗口对目标时间序列进行平滑滤波，如IDR迭代插值、Savitzky-Golay(SG滤波)等方法，这两者有更高的计算效率，但其不同参数的选择对结果影响较大。第三类是频率域的方法，包括谐波分析(HANTS)、小波变换等方法，他们的参数确定也是一个较大的问题。此外，还有其他不能归到以上三类方法，如基于数据同化的方法法、Whittaker滤波等。
这些时域滤波方法因模型简单、效率较高，在植被相关研究中得到了广泛应用，但它们通常存在难以处理时间上连续缺失的问题，从而无法得到高质量NDVI时间序列数据。
发明内容
本发明提出了一种结合矩阵补全和趋势滤波的植被指数时间序列重建方法，用以解决现有技术中无法得到高质量NDVI时间序列数据的技术问题。由于本发明提出的方法考虑了NDVI时间序列数据物候周期性和平滑性特点，因而可以得到高质量NDVI时间序列数据。
本发明所采用的技术方案是：结合矩阵补全和趋势滤波的植被指数时间序列重建方法，该方法包括以下步骤：
S1：输入植被指数时间序列原始数据y以及与之对应的质量标记数据f，质量标记数据用以表征植被指数时间序列原始数据中每个像素的质量；
S2：针对每个像素的植被指数序列，通过矩阵变化从一维向量变成二维矩阵Y，其中，两个维度分别代表时间序列的年份覆盖长度和一年内的监测频次，用以表征植被的季节性变化特征及年际周期变化相似性，构建的二维矩阵为NDVI时间序列数据矩阵；然后基于质量标记数据，找出二维矩阵中的缺失数据，将对应位置记录为缺失；
S3：针对每个像素变换后的二维矩阵，建立低秩矩阵补全的最优化能量方程，通过非精确增广拉格朗日算法实现矩阵补全，得到初步不含数据缺失的时间序列补全矩阵X；
S4：将补全矩阵X向量化，得到不含缺失的NDVI时间序列数据x；
S5.基于不含缺失的NDVI时间序列数据和质量标记数据f，建立加权趋势滤波的优化方程，通过交替方向乘子法实现噪声滤除，得到平滑干净的高质量NDVI时间序列数据z，作为重建后的NDVI时间序列数据。
在一种实施方式中，步骤S1中，输入的植被指数时间序列原始数据y有n年，每年具有k个时间序列点，y和f的总长度为kn。
在一种实施方式中，当植被指数时间序列原始数据y为MODIS植被指数产品MOD13中的数据时，NDVI数据分为干净像元、可能干净的像元、云覆盖像元以及雪覆盖像元四类，质量标记分别为0,1,2,3。
在一种实施方式中，步骤S2包括：
S2.1：将植被指数时间序列原始数据
矩阵化，形成k行n列的矩阵
其中每列为一年的NDVI时间序列数据，包括k个时间序列点数据，共有n列；
S2.2：根据质量标记数据f，找出其中被云覆盖以及雪覆盖的像元，并将找出的像元标记为缺失像元，矩阵Y中的对应元素记录为缺失元素，将其他像元标记为干净像元，矩阵Y中的对应元素记录为干净元素。
在一种实施方式中，步骤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Ω(·)为矩阵投影算子，S1/μk(·)为奇异值阈值算子，通过上述迭代求解获得的最优解X,为NDVI时间序列数据补全矩阵。
在一种实施方式中，步骤S4包括：
将NDVI时间序列数据补全矩阵
转变成时间序列向量
NDVI时间序列数据补全矩阵中不含缺失数据。
在一种实施方式中，步骤S5包括：
S5.1：根据质量标记数据f，将干净像元的权重设置为1，将缺失重建像元的权重设置为0.8，建立加权趋势滤波的优化方程为：
式中，
为待求的去噪结果，λ为正则化参数，W为对角矩阵，对角元素为对应像元的权重，D为二阶差分矩阵；
S5.2：利用交替方向乘子法求解，首先将式(6)转化成带约束的模型：
再构建式(7)的增广拉格朗日函数，即：
式中，m为拉格朗日乘子，β为惩罚参数，式(6)通过下述交替迭代步骤获得最优解，
其中，θ表示一个大于1的常数，用于更新惩参数β，
本申请实施例中的上述一个或多个技术方案，至少具有如下一种或多种技术效果：
(1)NDVI时间序列数据具有明显的物候周期性特点，因此NDVI时间序列数据构建成的矩阵具有显著的低秩性特点，从而通过低秩矩阵补全可以实现对缺失值的填，本发明能有效地填补NDVI时间序列中的缺失值，特别是对于时序连续缺失值有很好的补全效果。
(2)此外，考虑NDVI时间序列数据具有平滑性特点，因此通过加权趋势滤波进一步地滤除噪声，本发明能有效地去除NDVI时间序列中的噪声，获得平滑的NDVI时间序列结果，同时保持最大值点、最小值点等关键点。
总之，本发明提出的结合低秩矩阵补全和加权趋势滤波的NDVI时间序列数据重建方法能有效提高NDVI时间序列数据的质量，有很强的稳定性和通用性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例中基于低秩矩阵补全和加权趋势滤波的时序NDVI重建方法流程图；
图2为本发明实施例中时间序列NDVI矩阵的奇异值变化图；
图3为本发明实施例中时间序列NDVI数据原始数据和重建结果对比图。
具体实施方式
本申请发明人通过大量的研究与实践发现：现有技术中的时域滤波方法因模型简单、效率较高，在植被相关研究中得到了广泛应用，但它们通常存在难以处理时间上连续缺失的问题。这是因为这些方法是对每个像素的植被指数序列进行单独处理，因此当时间序列中存在时间连续缺失时候，目标像素很难有充分的参考进行有效的重建，从而难以获取理想的结果。鉴于此，有必要发展新的时域滤波方法，在保证处理效率的同时，能够有效的保证不同缺失情况的处理精度。通过结合低秩矩阵补全和加权趋势滤波，对植被指数时间序列进行整体建模，同时通过矩阵变化有效利用时间序列中其他年份的参考信息，能够实现植被指数时间序列的高质量、高效处理。
本发明提供了一种结合低秩矩阵补全和加权趋势滤波的时间序列NDVI数据重建方法，考虑了NDVI时间序列数据物候周期性和平滑性特点，能有效填补NDVI时间序列数据中的缺失数据并滤除噪声，获得高质量的NDVI时间序列数据。
为了达到上述目的，本发明的主要构思如下：
针对长时序遥感植被指数序列，首先通过植被指数产品中的质量标记数据确定时间序列数据中的缺失位置，然后针对单个像素的植被指数序列，通过矩阵变化从一维向量变成二维矩阵，两个维度分别代表时间序列的年份覆盖长度以及一年内的监测频次，即表征了植被的季节性变化特征及年际周期变化相似性，从而能够在除了时间邻域外给缺失像元提供更多的参考信息，实现不同缺失情景下的高精度重建。针对每个像素变换后的矩阵，建立低秩矩阵补全的最优化能量方程(见公式(1))，通过非精确增广拉格朗日算法实现矩阵补全，得到初步不含数据缺失的时间序列补全矩阵。最后再将该补全矩阵进行向量化，在一维向量的基础上建立加权趋势滤波的能量优化方程(见公式(6))，通过交替方向乘子法实现模型的求解，从而进一步滤除残留的噪声，得到平滑干净的高质量植被指数时间序列数据，实现长时间遥感植被指数序列的高精度重建。
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明实施例提供了结合矩阵补全和趋势滤波的植被指数时间序列重建方法，包括：
S1：输入植被指数时间序列原始数据y以及与之对应的质量标记数据f，质量标记数据用以表征植被指数时间序列原始数据中每个像素的质量；
S2：针对每个像素的植被指数序列，通过矩阵变化从一维向量变成二维矩阵Y，其中，两个维度分别代表时间序列的年份覆盖长度和一年内的监测频次，用以表征植被的季节性变化特征及年际周期变化相似性，构建的二维矩阵为NDVI时间序列数据矩阵；然后基于质量标记数据，找出二维矩阵中的缺失数据，将对应位置记录为缺失；
S3：针对每个像素变换后的二维矩阵，建立低秩矩阵补全的最优化能量方程，通过非精确增广拉格朗日算法实现矩阵补全，得到初步不含数据缺失的时间序列补全矩阵X；
S4：将补全矩阵X向量化，得到不含缺失的NDVI时间序列数据x；
S5.基于不含缺失的NDVI时间序列数据和质量标记数据f，建立加权趋势滤波的优化方程，通过交替方向乘子法实现噪声滤除，得到平滑干净的高质量NDVI时间序列数据z，作为重建后的NDVI时间序列数据。
如图1所示，为本发明提供的一种结合低秩矩阵补全和加权趋势滤波的NDVI时间序列数据重建方法的流程图。
具体实施过程中，植被指数时间序列原始数据y为NDVI时间序列数据，质量标记数据f标记了y中每个像素的质量情况。
S2：通过植被指数产品中的质量标记数据确定时间序列数据中的缺失位置，然后针对单个像素的植被指数序列，通过矩阵变化从一维向量变成二维矩阵(NDVI时间序列数据矩阵)，其中的两个维度分别代表时间序列的年份覆盖长度以及一年内的监测频次，即表征了植被的季节性变化特征及年际周期变化相似性，从而能够在除了时间邻域外给缺失像元提供更多的参考信息，实现不同缺失情景下的高精度重建。
步骤S3针对每个像素变换后的矩阵，建立低秩矩阵补全的最优化能量方程(见公式(1))，通过非精确增广拉格朗日算法实现矩阵补全，得到初步不含数据缺失的时间序列补全矩阵。再通过步骤S4将该补全矩阵进行向量化，最后通过步骤S5在一维向量的基础上建立加权趋势滤波的能量优化方程(见公式(6))，通过交替方向乘子法实现模型的求解，从而进一步滤除残留的噪声，得到平滑干净的高质量植被指数时间序列数据，实现长时间遥感植被指数序列的高精度重建。
请参见图2，为本发明实施例中时间序列NDVI矩阵的奇异值变化图。
从图2可以看出，矩阵Y的奇异值曲线下降非常快，因此矩阵Y具有显著的低秩特性，从而可以通过低秩矩阵补全重建缺失数据。建立低秩矩阵补全的优化方程，通过非精确增广拉格朗日算法实现矩阵补全，得到NDVI时间序列数据补全矩阵X。
在一种实施方式中，步骤S1中，输入的植被指数时间序列原始数据y有n年，每年具有k个时间序列点，y和f的总长度为kn。
NDVI时间序列数据以年为周期呈现出周期相似性特点。
在一种实施方式中，当植被指数时间序列原始数据y为MODIS植被指数产品MOD13中的数据时，NDVI数据分为干净像元、可能干净的像元、云覆盖像元以及雪覆盖像元四类，质量标记分别为0,1,2,3。
在一种实施方式中，步骤S2包括：
S2.1：将植被指数时间序列原始数据
矩阵化，形成k行n列的矩阵
其中每列为一年的NDVI时间序列数据，包括k个时间序列点数据，共有n列；
S2.2：根据质量标记数据f，找出其中被云覆盖以及雪覆盖的像元，并将找出的像元标记为缺失像元，矩阵Y中的对应元素记录为缺失元素，将其他像元标记为干净像元，矩阵Y中的对应元素记录为干净元素。
在一种实施方式中，步骤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Ω(·)为矩阵投影算子，S1/μk(·)为奇异值阈值算子，通过上述迭代求解获得的最优解X,为NDVI时间序列数据补全矩阵。
具体来说，由于步骤3.1建立的优化方程是非凸的，因此通过S3.2用矩阵的核范数最小化代替矩阵的秩最小化，将式(1)转化成凸优化问题。
S3.3中的步骤1～6为求解X的算法。
在一种实施方式中，步骤S4包括：
将NDVI时间序列数据补全矩阵
转变成时间序列向量
NDVI时间序列数据补全矩阵中不含缺失数据。
在一种实施方式中，步骤S5包括：
S5.1：根据质量标记数据f，将干净像元的权重设置为1，将缺失重建像元的权重设置为0.8，建立加权趋势滤波的优化方程为：
式中，
为待求的去噪结果，λ为正则化参数，W为对角矩阵，对角元素为对应像元的权重，D为二阶差分矩阵；
S5.2：利用交替方向乘子法求解，首先将式(6)转化成带约束的模型：
再构建式(7)的增广拉格朗日函数，即：
式中，m为拉格朗日乘子，β为惩罚参数，式(6)通过下述交替迭代步骤获得最优解，
其中，θ表示一个大于1的常数，用于更新惩参数β，
由于优化方程(6)不能直接获得解析解，因此利用交替方向乘子法求解，首先将式(6)转化成带约束的模型：
请参见图3，为本发明实施例中时间序列NDVI数据原始数据和重建结果对比图。
应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。
Claims (7)
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的常数，用于更新惩参数β，
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210037286.5A CN114463202A (zh) | 2022-01-13 | 2022-01-13 | 结合矩阵补全和趋势滤波的植被指数时间序列重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210037286.5A CN114463202A (zh) | 2022-01-13 | 2022-01-13 | 结合矩阵补全和趋势滤波的植被指数时间序列重建方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114463202A true CN114463202A (zh) | 2022-05-10 |
Family
ID=81410122
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210037286.5A Pending CN114463202A (zh) | 2022-01-13 | 2022-01-13 | 结合矩阵补全和趋势滤波的植被指数时间序列重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114463202A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114936202A (zh) * | 2022-05-16 | 2022-08-23 | 中国地质大学（武汉） | 一种极地反照率遥感数据的重构方法、装置及计算机设备 |
| CN117388872A (zh) * | 2023-09-05 | 2024-01-12 | 武汉大学 | 一种北斗地基增强***参考站坐标框架维持方法和*** |
-
2022
- 2022-01-13 CN CN202210037286.5A patent/CN114463202A/zh active Pending
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114936202A (zh) * | 2022-05-16 | 2022-08-23 | 中国地质大学（武汉） | 一种极地反照率遥感数据的重构方法、装置及计算机设备 |
| CN114936202B (zh) * | 2022-05-16 | 2023-07-11 | 中国地质大学（武汉） | 一种极地反照率遥感数据的重构方法、装置及计算机设备 |
| CN117388872A (zh) * | 2023-09-05 | 2024-01-12 | 武汉大学 | 一种北斗地基增强***参考站坐标框架维持方法和*** |
| CN117388872B (zh) * | 2023-09-05 | 2024-03-19 | 武汉大学 | 一种北斗地基增强***参考站坐标框架维持方法和*** |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109102477B (zh) | 一种基于非凸低秩稀疏约束的高光谱遥感图像恢复方法 | |
| CN110378858B (zh) | 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法 | |
| CN109102469B (zh) | 一种基于卷积神经网络的遥感图像全色锐化方法 | |
| CN114463202A (zh) | 结合矩阵补全和趋势滤波的植被指数时间序列重建方法 | |
| CN108133465B (zh) | 基于空谱加权tv的非凸低秩松弛的高光谱图像恢复方法 | |
| CN106897707B (zh) | 基于多源中分的特征影像时间序列合成方法及装置 | |
| CN107341837B (zh) | 基于影像金字塔的栅格-矢量数据转换及连续尺度表达方法 | |
| CN109859110A (zh) | 基于光谱维控制卷积神经网络的高光谱图像全色锐化方法 | |
| CN114119444A (zh) | 一种基于深度神经网络的多源遥感图像融合方法 | |
| CN104504740A (zh) | 压缩感知框架下的图像融合方法 | |
| CN109657081A (zh) | 高光谱卫星遥感数据的分布式处理方法、***及介质 | |
| CN102542296A (zh) | 一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法 | |
| CN114820741B (zh) | 一种高光谱影像全波段超分重建方法 | |
| CN113706418B (zh) | 一种基于光谱分离的长波红外遥感图像恢复方法 | |
| CN111598822A (zh) | 一种基于gfrw与iscm的图像融合方法 | |
| CN113421198B (zh) | 一种基于子空间的非局部低秩张量分解的高光谱图像去噪方法 | |
| CN111399041A (zh) | 一种小波紧框架自适应稀疏三维地震数据重构方法 | |
| Guo et al. | A flexible object-level processing strategy to enhance the weight function-based spatiotemporal fusion method | |
| CN111062888B (zh) | 一种基于多目标低秩稀疏及空谱全变分的高光谱影像去噪方法 | |
| CN113569772A (zh) | 遥感影像耕地实例掩膜提取方法及***、设备、存储介质 | |
| CN116720041A (zh) | 一种基于地表覆盖类型的精细化ndvi数据重建方法 | |
| Banothu et al. | High-order total bounded variation approach for gaussian noise and blur removal | |
| CN115034978A (zh) | 基于子空间表示结合图嵌入约束的高光谱图像去噪方法 | |
| CN113191958B (zh) | 一种基于鲁棒张量低秩表示的图像去噪方法 | |
| CN115439345A (zh) | 跨轨宽幅扫描相机的条带噪声抑制方法及*** |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |