高分辨率反射率時間序列重建算法研究

摘" 要：該文主要圍繞著高分辨率植被反射率時間序列重建過程中所遇到的數據噪聲多、數據缺失量大等問題，提出基于氣象數據和物候信息的高分辨率反射率時間序列重建算法（HTRA）。HTRA基于重建預測模型構建、時序數據去噪算法研究、重建策略建立為主要內容，以神經網絡和多源數據耦合為主要技術，利用真實值進行重建。以華北平原的一個區域為研究區，重建得到該區域農田的紅光反射率和近紅外反射率2016—2019年時空連續數據。通過對重建結果進行定性的分析驗證可以得到，HTRA算法具備很強的抗噪、抗缺失能力，重建出的反射率能夠很好的重建出植被本身的物候特征。

關鍵詞：反射率；時間序列重建；高分辨率；農田；時空連續

中圖分類號：P41" " " 文獻標志碼：A" " " " " 文章編號：2095-2945（2023）21-0001-08

Abstract： This paper focusing on the problems encountered in the reconstruction of high-resolution vegetation reflectance time series， such as excessive data noise and large amount of data missing， a high-resolution/reflectivity time-series reconstruction algorithm （HTRA） based on meteorological data and phenological information is proposed. HTRA is based on reconstruction prediction model construction， temporal data de-noising algorithm research， reconstruction strategy establishment as the main content， neural network and multi-source data coupling as the main technology， and uses the real value for reconstruction. Taking an area of the North China Plain as the study area， the spatial-temporal continuous data of the red light reflectance and near-infrared reflectance of the farmland in the area from 2016 to 2019 were reconstructed. Through qualitative analysis and verification of the reconstruction results， it can be concluded that the HTRA algorithm has strong anti-noise and anti-deletion ability， and the reconstructed reflectivity can well reconstruct the phenological characteristics of the vegetation itself.

Keywords： reflectivity; time series reconstruction; high resolution; farmland; spatio-temporal continuity

地表反射率是描述地球表面特性的重要參數，也是許多生物物理參數反演的關鍵參數，如葉面積指數（LAI）、植被光合有效輻射吸收系數（FAPAR）和地表反照率等[1]。根據中分辨率成像光譜儀（MODIS）、Landsat中高分辨率衛星等傳感器數據，已經生產許多陸地表面反射率產品[2-4]。地面反射率是大氣層頂反射率經大氣校正后計算而來的[5]。對于反射率的時間序列分析應用廣泛，可以用于變化檢測[6]、季節效應研究[7]、物候分析[8]等。但是由于云雨等因素影響導致地表反射率產品在時間和空間上不連續，為后續應用帶來挑戰[9]，因此重建地表反射率時間序列是當務之急。

中低分辨率衛星重訪周期短大多數在一天或者接近一天，具有較多的觀測值，可以利用插值方法生成規則采樣的時間序列。一些學者開展了低分辨率反射率數據時間序列重建的研究[10]。為了滿足精細化觀測、定量化描述的需求，高時空分辨率的產品生產技術迅速發展，高空間分辨率反射率數據的時間序列重建是其基礎。以Landsat系列為代表的中高分辨率衛星，重訪周期為16 d，有較長的重訪周期，稀疏的觀測數據為時間序列重建的插值方法帶來挑戰。Fisher等[11]和 Melaas 等[12]匯總多年的觀察資料，得到粗分辨率 MODIS 數據與 Landsat 數據融合的重建方法，并且已經顯示出較為良好的性能。但是在遇到在關鍵物候時期系統性缺失或低分辨率數據不存在的情況下該方法的應用受到限制。Zhu等[13]提出一種不依賴于其他傳感器的諧波建重建方法，通過 Landsat 進行了季節后云層篩選測試[14]，以及土地覆蓋監測測試[15]。但是獲得可靠結果的最低抽樣要求為12個樣本，這樣的要求在缺失較為嚴重的區域無法得到滿足。Pouliot等[16]通過分段諧波結合氣象站點和AVHRR數據進行建模，從而對Landsat的反射率和NDVI時間序列進行了重建。結果表明，在觀測稀疏的情況下，加入先驗知識可以改善模型預測。在±20 d內觀測值少于3次的情況下，比插補方法表現更好，平均反射率誤差減少0.001～2.5。

考慮到氣象因素對于植被生產的驅動機制，結合氣象數據輔助生成植被參數連續時間序列也有很多相關研究。Yu等[17]提出了一種基于氣候數據的間隙填充法（CGF），用于重建基于Landsat的NDVI時間序列。CGF模型考慮了NDVI時間序列和相鄰2年的氣候條件。利用太陽短波向下輻射、降水量和溫度等氣候變量來約束植被生長的年際變化。采用徑向基函數網絡（RBFN）描述NDVI年際序列之間的映射關系。研究結果表明，該方法能準確地重建和生成30 m連續8 d NDVI時間序列。當原始時間序列存在長時間缺失時，CGF方法優于傳統的時間序列重建方法（如Hants和SG濾波）。Zhu等[18]等提出了一種LAI時間序列重建算法MBPNN。通過在全連接神經網絡模型（FC-NN）中引入與增長相關的氣候信息和其他先驗知識來解決時間序列中存在連續缺口的問題。通過模擬不同的數據丟失場景，比較了MBPNN算法和HANTS算法的精度。結果表明，當時序間隙數從6（RMSE=0.184 1）上升到24（RMSE=0.807 3）時，HANTS算法的RMSE迅速增加，而MNNR算法的RMSE保持在0.2以下。李松澤[19]提出了一種全球LAI時間序列連續產品生產的MRA算法 ，通過考慮氣象數據的滯后效應以及對先驗知識數據進行去噪規范化處理的方法重建出全球植被的時空連續LAI。結果表明，該方法重建出的LAI能很好反映出植被生長物候規律。

已有研究表明，高分辨率反射率數據的重建算法主要基于濾波與融合。系統性缺失情況影響應用受限。在植被參數的重建中表明在模型中加入氣象數據，可以補充連續丟失的時間序列信息，輔助描述植被物候軌跡，有利于提高時間序列的重建精度、降低時間序列重建的不確定性，更加符合物候科學規律。本文基于此思路開展地表反射率的重建。考慮到紅光波段和近紅外波段對于植被參數的重要性，本文針對這2個波段的反射率重建開展研究分析。由于不同植被類型的物候特征差異大，且農田具備時空連續的特征對于作物生長研究具備重要意義，本文選擇以農田作為主要植被類型開展研究。

1" 研究區及數據

1.1" 研究區概況

本次研究的區域選在華北平原地區，位于北緯32°～40°，東經114°～121°，作為我國重要的糧棉油基地，華北平原在前一年年末至本年5、6月主要種植冬小麥，6—10月期間種植玉米，同時還種植了較大規模的棉花、甜菜等經濟作物。年平均含云量在50%～60%，這樣的研究區條件可以滿足我們對農作物反射率重建的研究。

1.2" 研究數據

考慮到Landsat數據具有質量標識文件，對于數據質量具有良好的控制效果。訓練數據集本研究采用USGS Landsat 8/7 Level2的紅和近紅外反射率波段。研究采用的土地覆蓋數據為GLC30土地覆蓋數據集，該數據集由中國國家地理信息中心（NGCC）基于像元-對象-知識（POK）的分類策略，將全球地表劃分為耕地、森林、草地、灌叢、濕地、開闊水域、人工地表、裸地和永久冰雪9類，總體分類準確率大于72.51%。

ERA5數據由歐洲中期天氣預報中心生成，由再分析網格化數據集組成，研究中使用了總降水量、地表太陽輻射下降量和地表溫度3種數據。將空間分辨率較低的ERA5數據投影并重采樣到30 m，以匹配幅Landsat ETM圖像。去除了無效和異常的觀測，并將每小時的數據聚合為8 d，以減少噪聲的影響。

2" 高分辨率反射率時間序列重建算法

基于氣象數據和物候信息的高分辨率反射率時間序列重建算法（以下稱為HRTA）是在低分辨率LAI重建算法（MRA）[19]的基礎上改進得來的。MRA算法主要利用氣象數據基于神經網絡模型對低分辨率LAI數據進行時間序列重建得到符合植被物候特征的時空連續LAI。與低分辨率的LAI重建不同的是：①高分辨率反射率數據由于重訪周期相對較長，缺失更多；②與LAI相比地表反射率對噪聲的響應更顯著；③地表反射率和LAI對氣象數據的響應不同。針對上述3點，研究受殘云影響顯著的反射率數據去噪方法，重新制定重建策略，構建了基于神經網絡和氣象數據的反射率時序列重建模型。具體重建流程如圖1所示。由圖1可知，①數據預處理：對樣本點的二次篩選與數據均衡化、待重建原始數據的有效值選取；②神經網絡預測模型構建：基于相關性分析的訓練數據集構建、神經網絡模型訓練；③利用重建模型和策略結合模型輸入得到重建結果。

2.1" 數據預處理

2.1.1" 樣本點選取與數據均衡化

神經網絡在時間序列重建中具有良好的應用[18-19]，因此在反射率重建中同樣也采用神經網絡的技術構建預測模型。訓練樣本的質量極大地影響著重建模型的應用效果，我們對樣本質量的控制主要體現在對樣本純度和樣本數據完整性的控制。樣本純度是關鍵因素之一，摻雜有其他樣本，尤其是生長模式不一致的訓練樣本會使得重建結果存在潛在偏差。樣本數據的完整性就是指樣本點在時間序列上數據的缺失情況，因為受到天氣影響有的數據往往會有系統性缺失，會導致原始樣本數據在不同物候時期上分布不均衡。

在訓練樣本選擇時，主要分為兩步，首先根據土地覆蓋產品GLC30進行隨機選點在中國境內選取了2013—2021年的4 086個農田覆蓋類型的樣本點時間序列數據。但是選點會由于土地類型突變、選擇數據缺失嚴重等各種原因導致所選點并不滿足要求；接下來進行第二步的篩選，考慮到農田NDVI的數值和時序圖的形狀具有規律性。于是逐點繪制農田NDVI時序圖，根據NDVI的數值大小和波動情況篩選訓練樣本。最后計算多年平均數據的缺失率，當缺失率超過10%則樣本點也不符合要求。

由于植被參數原始數據集存在系統性缺失，原始樣本數據在不同物候時期上不可避免會出現分布不均衡的情況。并且，大氣污染還會使原始反射率的時間序列數據中存在噪聲。MRA算法中運用HANTS濾波對植被參數進行了去噪均衡化處理得到了均衡的數據效果。HANTS的核心算法是最小二乘法和傅立葉變換，通過最小二乘法的迭代擬合去除時序反射率值中受云污染影響較大的點，借助于傅立葉在時間域和頻率域的正反變換實現曲線的分解和重構，從而達到時序數據去噪填補的目的。本研究也采用HANTS濾波對反射率數據進行均衡化處理。但是HANTS算法在缺失間隙較大的情況下容易算法失效。因此上文介紹的的樣本點篩選步驟就顯得更加重要。在缺失率較少情況下HANTS算法具備較好的去噪填補的效果，最終使本研究的原始反射率數據在不同物候期分布均衡，受到明顯噪聲影響更少。

2.1.2" 數據去噪

待重建的地表反射率數據，是反射率時序重建的原始數據，也是重建過程中模型參考的基準值。由于受到殘云及傳感器噪聲等影響，反射率數據會存在顯著的噪聲，算法的重建效果依賴于保留反射率值的觀測質量。因此，需要對原始反射率數據進行較嚴格的去噪。MRA算法利用包絡線方法進行篩選，但由于中高分辨率的反射率時間序列中缺失數據多、相鄰數據間隔大等原因，會使得包絡線所依賴的窗口間隔過大，出現錯誤選點、過度選點等問題，方法不再適用。本去噪算法以多年平均反射率數據為先驗知識，基于時間窗口判定、包絡線約束進行原始反射率數據有效值選取。

首先進行時間窗口大小判定，根據不同的缺失窗口大小，采用不同的有效值判定方法。其次為了時間序列上的連續分析，把無效值進行了有效化賦值參與后續篩選。

對NDVI當相鄰有效值相差不超過5個單位（每個單位為8 d）的，進行基于滑動窗口的上包絡線求取，通過和滑動窗口內的最大值做比較判定值的有效性。當相鄰有效值相差超過5個單位則采用和對應時間點的多年平均數據的比較，當小于多年平均值一定范圍內判定無效。

對紅光反射率當相鄰有效值相差不超過5個單位（每個單位為8 d）的，進行基于滑動窗口的下包絡線求取，通過和滑動窗口內的最小值做比較判定值的有效性。當相鄰有效值相差超過5個單位則采用和對應時間點的多年平均數據的比較，當大于多年平均值一定范圍內判定無效。

2.2 時序預測神經網絡模型構建

為了得到一個良好的神經網絡模型構建具有強相關性的訓練數據集是相當重要的。根據Zhu等[18]、李松澤[19]等對LAI重建的研究，氣象數據是植被生長中的關鍵因素，并且多年歷史數據也在重建中發揮重要作用，通過差值預測差值可以更好地得到重建數據的物候特征。因此本研究在植被反射率重建中加入氣象因子、反射率多年平均的變化值、采用變化值作為模型的輸出。

氣候因子解釋了64%全球植被生長的變化，且考慮滯后效應會有更高的數據相關性。研究表明全球植被生長的首要驅動因素是氣溫，其次是降水和太陽輻射。因此得到這3種氣象數據對反射率數據的滯后情況，會有益于構建出一個和反射率變化量更加具有相關性的數據集。

通過訓練數據集中對氣象數據的相關性分析得到了具有強相關性的訓練數據集（表1、表2），然后需要進行神經網絡的搭建，全連接層神經網絡具有非常好的非線性映射能力、較高的并行性及全局優化的特點[20]，本研究選擇全連接層神經網絡作為模型訓練網絡。考慮到時序預測中往往具有雙向預測的特征，本文得到了以時間序列前向點為基準向后一個時相點進行預測的前向模型，和以時間序列后向點為基準向前一個時相點進行預測的后向模型。

2.3" 基于重建策略的時序缺失場景重建

基于前向模型和后向模型得到的重建結果需要針對不同重建情形進行具體分析，因此需要分情形制定重建策略進行應用。在應用中會主要遇到3種情況，第一種情況就是當缺失時相數據位于時間序列的最前面，沒有前向參考值；第二種情況就是當缺失時相位于時間序列中間，既有前向參考值也有后向參考值；第三種情況就是當缺失時相位于時間序列最后，沒有后向參考值。具體情況如圖2所示。圖中為模擬的一個反射率時間序列，并展示出了運用重建策略后的重建結果。

第一種情況下，當缺失時相數據位于時間序列的最前面，由于沒有前向參考值，利用后向模型進行重建。第二種情況下，當缺失時相位于時間序列中間，既有前向參考值也有后向參考值，為了利用好兩個模型在重建中的作用，本研究對前向模型和后向模型重建出的值基于前向點后向點到重建值的時間間隔，進行反距離加權得到時間影響的模型應用權重。得到前向模型和后向模型應用權重后，可以得到前向模型和后向模型結果加權后的重建反射率值。第三種情況下，當缺失時相數據位于時間序列的最末端，由于沒有后向參考值，利用后向模型進行重建。

3 結果與分析

3.1 重建精度評價

本研究尋找缺失較少的農田樣點，得到多年平均值時間序列（缺失率小于5%）并進行濾波去噪處理。得到的時序數據作為參考值進行評估，主要評估HRTA算法面對大量缺失場景以及大量噪聲場景時的重建能力。

在MuSyQ 高分 16 m空間分辨率 10 d合成的NDVI 植被指數產品中指出南方云雨天氣較多往往數據缺失量更大[21]。川黔地區年平均總云量均在80%以上，部分地區可達90%以上[22]，是中國區域云量最多的地區。因此本研究模擬了大量缺失的情況，即90%左右的缺失率，且包括關鍵物候期數據缺失。對于大氣帶來的噪聲影響，因為大氣對于紅光反射率具備正偏置效應，對于近紅外反射率負偏置效應更加突出。因此，在原有的大量缺失情景下加入本文所模擬的噪聲數據。把模擬好的2種情景數據作為本研究的待重建數據。

在大量缺失情景下，對其進行重建最終得到重建結果。如圖3（a）和圖3（c）依次為某一樣點的紅光反射率、近紅外反射率的重建結果。由圖中可以得到在年中的關鍵物候期數據缺失情況下依然可以重建出符合其生長的物候特征。在存在大量噪聲情景下，對其進行重建最終得到重建結果。如圖3（b）和圖3（d）依次為該樣點的紅光反射率、近紅外反射率的重建結果。由圖中可以得到在存在大量噪聲的情景下利用HRTA算法重建，重建結果依然可以和參考反射率具備較好的一致性。

3.2" 時間序列重建結果

本研究基于HRTA算法對華北平原的一個地區的農田地類像進行紅光反射率和近紅外反射率時間序列重建，得到圖4的空間分布重建結果。其中圖4（a）和圖4（c）分別為反射率時間序列重建前紅光反射率、近紅外反射率的原始空間分布，圖4（b）和圖4（d）分別為經過時間序列重建后的兩者的空間分布。圖中的黑色代表缺失或者區域內的非農田地類。圖中001—361是在該地區2018年的第1天到第361天以8 d為時間間隔的連續時間分布。

通過圖4中重建前后紅光反射率、近紅外反射率空間分布的對比，可以發現重建后的數據不僅保持了和原有數據高度一致的空間分布，并且在原有數據具有明顯噪音、大量甚至全部缺失情況下仍然可以重建得到符合農田空間分布規律及物候生長規律的數據。

為了直觀體現出重建過程中HRTA算法對單點的重建效果，本研究進行了點的時間序列重建結果展示。圖5（a）和圖5（b）分別代表對紅光反射率、近紅外反射率的重建結果。圖中包括重建的原始參考點、重建選取的有效值點以及重建前后的時間序列曲線。

由圖5可以觀察到重建后的紅光反射率和近紅外反射率具備良好的抵抗噪聲的能力。也可以明顯觀察到本文得到的紅光反射率具備良好且穩定的物候特征，并在儒略日2016185—2017001及2019185—2019361之間在關鍵物候點缺失的情況下重建出的結果和未缺失的以往年份得到的結果仍然保持趨勢相同大小一致，也證明本研究的重建結果具備良好的物候恢復能力以及抗缺失能力。

4" 結束語

本文對農田地類的30 m分辨率紅光反射率和近紅外反射率數據進行了時間序列重建。利用了氣象數據和歷史數據的物候信息結合相關性分析得到訓練數據集并結合全連接層神經網絡訓練得到神經網絡預測模型。并利用歷史數據為先驗知識，結合時序滑動窗口判定以及包絡線篩選進行時間序列數據去噪幫助本文得到重建中的有效點。利用前后向模型結合時間間隔隱含的時間信息得到了重建策略。基于神經網絡預測模型、去噪選點算法和重建策略構建了反射率時間序列重建模型，結合待重建區域的原始數據最終得到了建反射率。

重建的方法需要在更多地類進行測試和完善，進而建立全地類高分辨率反射率時間序列重建算法體系。對于算法精度的評價還需要針對抗噪能力進行測試，以及對在不同缺失梯度下算法的重建精度變化情況進行分析。

參考文獻：

[1] LIANG S，ZHANG X，XIAO Z， et al. Global LAnd Surface Satellite （GLASS） Products[M]. Springer International Publishing， 2014.

[2] VERMOTE E F ， SALEOUS N Z E ， JUSTICE C O .Atmospheric correction of MODIS data in the visible to middle infrared： first results[J]. Remote Sensing of Environment， 2002，83（1-2）：97-111.

[3] CLAVERIE M，VERMOTE E F，FRANCH B，et al. Evaluation ofthe Landsat - 5 TM and Landsat - 7 ETM + surface reflectance products[J]. Remote Sensing of Environment，2015（169）：390-403.

[4] YAN P ，HE G ，ZHANG Z ， et al. Study on atmospheric correction approach of Landsat-8 imageries based on 6S model and look-up table[J]. Journal of Applied Remote Sensing， 2016， 10（4）：045006.

[5] LIANG S， FANG H， CHEN M. Atmospheric correction of Landsat ETM+ land surface imagery. I. Methods[J]. IEEE Transactions on Geoscience amp; Remote Sensing， 2001， 39（11）：2490-2498.

[6] POULIOT， LATIFOVIC， FERNANDES， et al. Evaluation of annual forest disturbance monitoring using a static decision tree approach and 250 m MODIS data[J]. REMOTE SENS ENVIRON，2009，113（8）：1749-1759.

[7] ZHAO H， FERNANDES R. Daily snow cover estimation from Advanced Very High Resolution Radiometer Polar Pathfinder data over Northern Hemisphere land surfaces during 1982-2004[J].Journal of Geophysical Research Atmospheres，2009，114（D5）：113.

[8] POULIOT D， LATIFOVIC R， FERNANDES R， et al. Evaluation of compositing period and AVHRR and MERIS combination for improvement of spring phenology detection in deciduous forests[J]. Remote Sensing of Environment， 2011，115（1）：158-166.

[9] ZHU X， HELMER E H. An automatic method for screening clouds and cloud shadows in optical satellite image time series in cloudy regions[J]. Remote Sensing of Environment，2018（214）：135-153.

[10] ZHIQIANG X， SHUNLIN L， TONGTONG W， et al. Reconstruction of Satellite-Retrieved Land-Surface Reflectance Based on Temporally-Continuous Vegetation Indices[J]. Remote Sensing，2015，7（8）：9844-9864.

[11] FISHER J I， MUSTARD J F， VADEBONCOEUR M A. Green leaf phenology at Landsat resolution： Scaling from the field to the satellite[J]. Remote Sensing of Environment， 2006， 100（2）：265-279.

[12] MELAAS E K， FRIEDL M A， ZHU Z. Detecting interannual variation in deciduous broadleaf forest phenology using Landsat TM/ETM+ data[J].Remote Sensing of Environment： An Interdisciplinary Journal，2013：132.

[13] ZHU Z， WOODCOCK C E， HOLDEN C， et al. Generating synthetic Landsat images based on all available Landsat data： Predicting Landsat surface reflectance at any given time[J]. Remote Sensing of Environment，2015（162）：67-83.

[14] ZHU Z， WOODCOCK C E. Automated cloud， cloud shadow， and snow detection in multitemporal Landsat data：An algorithm designed specifically for monitoring land cover change[J].Remote Sens. Environ，2014（152）：217-234.

[15] ZHE Z， WOODCOCK C E. Continuous Change Detection and Classification of Land Cover Using All Available Landsat Data[J]. Remote Sensing of Environment，2013， 144（1）：152-171.

[16] POULIOT D， LATIFOVIC R. Reconstruction of Landsat time series in the presence of irregular and sparse observations： Development and assessment in north-eastern Alberta， Canada[J].Remote Sensing of Environment，2018（204）：979-996.

[17] YU W， LI J， LIU Q， et al. Gap Filling for Historical Landsat NDVI Time Series by Integrating Climate Data[J]. Remote Sensing， 2021， 13（3）：484.

[18] ZHU X，" LI J，" LIU Q， et al. Use of a BP Neural Network and Meteorological Data for Generating Spatiotemporally Continuous LAI Time Series[J]. IEEE Transactions on Geoscience and Remote Sensing， 2021（99）：1-14.

[19] 李松澤.結合氣象數據的全球植被葉面積指數遙感產品時間序列重建算法研究[D].北京：中國科學院空天信息創新研究院，2022.

[20] 張馳，郭媛，黎明.人工神經網絡模型發展及應用綜述[J].計算機工程與應用，2021，57（11）：57-69.

[21] 李松澤，李靜，于文濤，等.MuSyQ高分16米空間分辨率10天合成的NDVI植被指數產品（2018—2020年中國01版）［J］.中國科學數據（中英文網絡版），2022，7（1）：237-246.

[22] 林之光.貴州和四川盆地云量的氣候研究[J].地理學報，1986（4）：289-30.