[發明專利]結合矩陣補全和趨勢濾波的植被指數時間序列重建方法

在審

• 申請號: 202210037286.5
• 申請日: 2022-01-13
• 公開（公告）號: CN114463202A
• 公開（公告）日: 2022-05-10
• 發明（設計）人: 管小彬;儲棟;李星華;沈煥鋒
• 申請（專利權）人: 武漢大學
• 主分類號: G06T5/00
• 分類號: G06T5/00;G06T5/50
• 代理公司: 武漢科皓知識產權代理事務所(特殊普通合伙) 42222
• 代理人: 羅飛

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時間序列數據重建結果。