2000–2020 年珠三角30 m 空間分辨率植被覆蓋度時空變化數據集

沈明潭，譚炳香*，侯瑞霞，黃逸飛

1.中國林業科學研究院資源信息研究所，北京 100091

引 言

珠三角在國家整體發展中具有重要的戰略地位，自改革開放以來經歷了快速的經濟增長、人口增長和城市化擴張，導致當地植被急劇退化和破壞，從而嚴重影響了生態系統的服務和功能，因此動態監測珠三角植被的時空變化和分布格局對當地的生態建設具有重要的參考意義。在已公開發布包含珠三角地區的植被覆蓋度數據集產品中，高吉喜等基于空間分辨率為250 m 的MODI3Q1 NDVI（Normalized differential vegetation index）數據，采用NDVI 最大合成法的方式進行中國區域的植被覆蓋反演[1]。與該數據集相比，本研究選用空間分辨率更高的30 m Landsat5 TM 和Landsat8 OLI 數據，利用中值合成法獲取研究區的影像數據，更能精確估測出珠三角地區的植被覆蓋度，更科學有效地反映出珠三角地區的植被覆蓋狀況。

植被覆蓋度（Fraction Vegetation Cover, FVC）是指地面植被（包括樹葉、莖和樹枝）的垂直投影面積占統計區域總面積的百分比，可以反映地表植被的生長狀況，是監測生態環境的重要基礎數據，在評估地表覆蓋變化、土地利用方式、景觀和生態系統功能等方面有重要指示作用[2]，被廣泛用于監測植被變化、區域環境質量評估、氣候預測、土壤侵蝕評估等。遙感技術因其獲取速度快、覆蓋范圍廣、能夠提取大量信息而廣泛應用于研究植被動態監測、生態環境評估和土地覆蓋分類等領域[3]，是監測區域和全球FVC 最高效的方法[4]，可以動態監測植被覆蓋變化及其影響因素。通過長時間序列和大規模的FVC 研究來評估特定區域的生態環境條件和質量，在揭示區域環境條件的演變過程和預測未來變化趨勢等方面具有重要意義。

1 數據來源和處理方法

本研究采用以下主要步驟對珠三角植被覆蓋度進行分析：（1）在GEE（Google Earth Engine）平臺調用Landsat 數據，基于NDVI 像元二分法估算珠三角2000、2005、2010、2015 和2020 年植被覆蓋度；（2）利用變異系數分析珠三角5 個時期植被覆蓋度變化的波動性；（3）利用Slope趨勢分析來探究植被覆蓋度變化的時空變化特征；（4）采用Hurst 指數預測研究區植被覆蓋度未來的發展趨勢。

1.1 數據來源

Landsat5 TM 和Landsat8 OLI 數據具有較高的空間分辨率（30 m），但由于研究區地處亞熱帶區域，受到長時間大面積高云量的影響，難以獲取研究區完整清晰的影像。因此本研究基于GEE 平臺，利用一年中云量低于20%的數據，按照成像時間進行排列，選取位于中間時段的像元，合成出質量較高的研究區影像（見表1）。調用“LANDSAT/LT05/C01/T1_SR”數據集，分別利用2000-01-01 至2000-12-31 期間的141 景Landsat5 TM 影像合成2000 年珠三角年度中值影像，2005-01-01 至2005-12-31 期間的119 景Landsat5 TM 影像合成2005 年珠三角年度中值影像，2010 年的數據量較少，需要擴大時間段才能得出質量較好的影像數據，因此利用2009-11-01 至2010-12-31 期間的101 景Landsat5 TM 影像合成2010 年珠三角年度中值影像。調用“LANDSAT/LC08/C01/T1_SR”數據集，分別利用2015-01-01 至2015-12-31 期間的162 景Landsat8 OLI 影像合成2015 年珠三角年度中值影像，2020-01-01 至2020-12-31 期間的152 景Landsat8 OLI 影像合成2020 年珠三角年度中值影像。由于按照時間段內云量低于20%的數據進行篩選，每個時間段內符合的景數不同，且能完整合成出研究區清晰無云的影像數據需要的景數也不同，因而每個期間所調用的影像景數會有差別。所用的Landsat 數據集都已經過輻射定標和大氣校正，分別裁剪出珠三角5 個時期的合成影像（見圖1）。

圖1 珠三角Landsat8 標準假彩色影像Figure 1 Standard false colour images of Landsat8 in the Pearl River Delta

表1 影像合成信息Table 1 Image synthesis information

1.2 數據處理方法

1.2.1 植被覆蓋度估測

本研究基于GEE 平臺生成珠三角年度中值植被覆蓋度影像[5-6]。主要的計算過程為：（1）輸入珠三角范圍；（2）構建Landsat 遙感影像的時間序列數據集。分別調用LANDSAT/LT05/C01/T1_SR、LANDSAT/LC08/C01/T1_SR 數據集，根據時間段進行篩選數據，然后進行去云處理；（3）生成珠三角年度中值遙感影像數據集。對構建后的影像數據集按不同的年份進行分組，逐像元計算生成每個時期的年度中值影像；（4）計算NDVI。調用normalized Difference 函數計算NDVI（公式1）；（5）去除水體。為了避免水體對植被覆蓋度估測的影響，調用mask 函數將水體進行掩膜處理；（6）計算珠三角FVC。選用基于NDVI 的像元二分法（公式2），分別計算每個時期的植被覆蓋度；（7）導出計算結果。

1.2.2 歸一化差值植被指數

Kriegler 等人[7]提出了一種簡單的波段變換：用近紅外光反射率減去紅光反射率，再除以近紅外光反射率加紅光反射率，稱為歸一化差異植被指數（NDVI）。NDVI 對植被的生長狀態和空間分布特征非常敏感[8]，可以在很大程度上消除地形和大氣的干擾，并與植被覆蓋度有顯著的線性相關[9]，在大規模植被動態監測等方面具有明顯優勢，NDVI 值越高，植被覆蓋度越高，是目前應用最廣泛的植被指數[10]。計算公式為：

式中，NDVI是歸一化差異植被指數，NIR和R分別為傳感器在近紅外光和紅光區域記錄的反射率。

1.2.3 像元二分模型

植被覆蓋度與NDVI 之間具有顯著的線性相關[11]，利用植被指數與像元分解模型相結合的方法，具有一定的物理意義、反演算法簡單[12]，而且在很大程度上弱化了大氣、土壤背景等因素的影響，被廣泛應用于區域和全球范圍內植被覆蓋度提取[13-14]。因此本研究利用基于NDVI 的像元二分模型[15]對研究區植被覆蓋度進行估算：

式中FVC為植被覆蓋度，本研究選取置信區間為[5%, 95%]的NDVI 值作為NDVIsoil和NDVIveg[16]。NDVIsoil表示無植被覆蓋像元NDVI 值，NDVIveg表示純植被覆蓋像元NDVI 值。

1.2.4 變異系數

對研究區5 個時期的植被覆蓋度數據逐像元計算變異系數，探究研究區植被覆蓋度的穩定性[17]。CV 值越大，植被變化的波動性越大，將波動性劃分為4 種類型，見表2。計算公式如下所示：

表2 植被覆蓋度波動變化類型劃分表Table 2 Classification of FVC fluctuation types

式中CV植被覆蓋度的變異系數，n 是研究期數，FVCi第i 年的植被覆蓋度值，F???VC???是研究期內植被覆蓋度均值。

1.2.5 Slope 趨勢分析

采用一元線性回歸分析法中的斜率公式計算要素（y）的變化趨勢（Slope）并進行顯著性檢驗，根據Slope值和P 值對植被覆蓋度的變化趨勢進行劃分，將變化趨勢分為7 種類型[18]，見表3。其中Slope的計算公式如下：

表3 植被覆蓋度變化趨勢類型劃分表Table 3 Classification of the trend types of FVC changes

式中，Slope表示要素的年際變化趨勢；i為時期變量；n是研究期數，本研究為5 個時期；yi為第i年的要素。Slope>0 和Slope<0 分別表示要素序列隨著時間增加和減少，Slope的絕對值越大表明要素變化越大。

采用F 檢驗對變化趨勢進行顯著性檢驗，當檢驗結果P<0.05 為顯著變化，0.050.1 為不顯著變化。統計量計算公式如下：

式中，U為回歸平方和，Q為殘差平方和，yi為第i年的要素，y?i為要素的回歸值，y?為要素的平均值，n表示期數。

1.2.6 Hurst 指數

基于R/S 分析法的Hurst 指數[19]，能夠用于定量描述時間序列信息的長期相關性，目前已廣泛運用于分析區域植被覆蓋未來變化趨勢，Hurst 指數取值范圍存在以下3 種情況：（1）若0.5

時間序列{FVC(t)},t=1, 2, 3,…,n，對于任意一個正整數τ=1，定義均值系列：

累計離差：

極差：

標準差：

考慮比值R(τ)/S(τ)?R/S，若存在R/S∝τH，則說明時間序列{φ(t)}, t=1, 2, 3,…,n存在Hurst 現象，其中H 值稱為Hurst 指數。參考劉明霞等[21]將植被覆蓋度未來變化趨勢分為7 類，見表4。

表4 植被覆蓋度未來發展趨勢類型劃分表Table 4 Classification of future development trends of FVC

2 數據樣本描述

2.1 數據組成

本數據集包括：（ 1 ） PRD_FVC_2000.tif 、 PRD_FVC_2005.tif 、 PRD_FVC_2010.tif 、PRD_FVC_2015.tif、PRD_FVC_2020.tif，分別為2000 年、2005 年、2010 年、2015 年和2020 年5 個時期的年度植被覆蓋度數據；（2）PRD_FVC_CV.tif 為植被覆蓋度變異系數數據；（3）PRD_FVC_Slope_F.tif 為植被覆蓋度變化趨勢數據；（4）PRD_FVC_Hurst.tif 為植被覆蓋度未來變化趨勢數據。上述數據空間分辨率為 30 m，地理坐標系為 WGS1984，投影坐標系為WGS_1984_UTM_Zone_49N，數據保存為GeoTIFF 格式，整個數據集解壓后大小約為2.75 GB。

2.2 樣本描述

基于GEE 平臺，利用公式1 和公式2 進行珠三角植被覆蓋度反演，得到2000 年、2005 年、2010年、2015 年和2020 年5 個時期的年度植被覆蓋度影像（見圖2）?；谥槿? 個時期的植被覆蓋度反演結果，利用公式3 對珠三角植被覆蓋度變化的波動性進行分析，根據表2 將變異系數劃分為4 種等級，結果見圖2(a)；利用公式4 計算出珠三角5 個時期的植被覆蓋度的時空變化趨勢，采用F檢驗（公式5、6 和7）對變化趨勢進行顯著性檢驗，根據表3 的將變化趨勢劃分為7 個等級，結果見圖2(b)；利用公式4 計算出珠三角植被覆蓋度的變化趨勢，再利用Hurst 指數（公式8、9、10 和11）來反映植被覆蓋度變化趨勢的持續性，從而揭示珠三角植被覆蓋度未來的發展趨勢，根據表4 將珠三角發展趨勢劃分為7 種類型，結果見圖2(c)。

圖2 珠三角5 期植被覆蓋度（2000、2005、2010、2015、2020）、變異系數(a)、變化趨勢(b)、Hurst 指數(c）空間分布Figure 2 Spatial distribution of vegetation coverage in five phases (2000, 2005, 2010, 2015, 2020), coefficient of variation (a), trend (b) and Hurst index (c) in the Pearl River Delta

3 數據質量控制和評估

本數據集參考植被覆蓋度的相關研究成果，在GEE 平臺利用中值合成法得到研究區清晰無云的影像，提高反演數據源的質量和精度。植被覆蓋度的實際采樣點觀測值與像元二分法的估測值有較好的一致性，利用像元二分法對植被覆蓋度進行反演是科學可靠的方法[22-23]。因此選擇基于NDVI的像元二分法進行珠三角地區植被覆蓋度反演，根據研究區的氣候條件和植被覆蓋狀況，在5%的置信水平下，采用累計頻率在5%～95%的NDVI 值作為上下閾值，在一定程度上減少遙感影像的噪點影響[24]?；贜DVI 的二分像元模型對圖像輻射校正的影響不敏感，具有較高的精度和易操作性，不需要實際的植被覆蓋度數據建模，僅需簡單計算，易于解釋[18]，目前廣泛應用于區域、省、流域、國家或全球范圍等不同尺度的植被覆蓋變化研究。參考Song 等學者的研究[25]，利用相關系數R2和均方根誤差RMSE 進行產品的精度驗證和一致性檢驗。

3.1 精度驗證

本研究在珠三角范圍內隨機生成100 個采樣點，通過world imagery 2020-07-22 高分辨率影像進行精度驗證，將相同坐標位置的2020 年植被覆蓋度遙感反演值和高分辨率影像目測值進行回歸分析，兩者相關系數R2=0.90，均方根誤差RMSE=0.03，表明本數據集的精度基本滿足區域性植被覆蓋度研究的需要。

3.2 一致性檢驗

關于產品的一致性檢驗，目前公開發布的包含珠三角地區的植被覆蓋度產品都為粗分辨率的，比如GEOV2 的2000 m 植被覆蓋度產品，GLASS 的500 m 植被覆蓋度產品，高吉喜團隊的中國區域250 m 植被覆蓋度數據集（2000–2022）。本數據集以高吉喜團隊2020 年植被覆蓋度數據為對比產品，在珠三角范圍內隨機生成1000 個點，提取本研究2020 年植被覆蓋度數據和對比產品相同位置的值進行一致性檢驗，相關系數R2=0.69，均方根誤差RMSE=0.09，表明兩個產品整體上有較好的一致性，但根據world imagery 2020-07-22 影像可以看出，本研究的FVC 數據相較于對比產品表現出了更多的細節（見圖3）。一致性不高的原因有：（1）雖然都采用了像元二分法，但NDVI 的取值方式不一樣，對比數據集采用的是NDVI 最大值合成法，而本研究利用中值法合成影像，然后再計算NDVI；（2）兩個數據集的空間分辨率不一樣；（3）數據源不一樣。

圖3 world image 影像（a）、本研究2020 年FVC（b）、對比產品2020 年7 月FVC（c）Figure 3 world image data (a), 2020 FVC (b) of this study, July 2020 FVC (c) of comparative products

4 數據價值

本數據集提供了珠三角2000、2005、2010、2015 和2020 年5 個時期植被覆蓋度的反演結果、變異系數、變化趨勢和Hurst 指數，可用于定量研究珠三角植被覆蓋度的時空變化特征和預測未來變化趨勢，探究珠三角植被覆蓋度變化的機理，更好地了解珠三角植被變化的動態過程，對植被資源的管理及生態環境的保護具有重要的指導意義。在珠三角快速的城市化進程中，本數據集具有重要的實際應用價值。在城市規劃方面，可以評估城市化對植被生態系統的影響，有助于指導城市的規劃和發展。在生態保護建設方面，植被覆蓋度研究有助于科學地認識生態系統及生態過程，可以指導生態環境的保護和恢復工作。在氣候變化研究方面，監測植被變化可以深入了解氣候變化對生態系統的影響，為可持續發展提供關鍵信息。

5 數據使用方法和建議

本數據集以GeoTIFF 格式進行存儲，在常用的GIS 軟件（ArcGIS, QGIS 等）及其他數據分析軟件（R, MATLAB 等）都可進行查看和再分析，可以實現2000、2005、2010、2015 和2020 年珠三角植被覆蓋度的定量評價和時空變化特征分析，可以應用于珠三角植被覆蓋度動態變化與影響因子之間的相關性探究，以及定量分析自然因子和人為因子對珠三角植被覆蓋度變化的相對影響程度等方面的研究。植被覆蓋度的取值范圍為[0,1]，為了監測地表植被覆蓋細微的時空變化，同時參考相關的植被覆蓋度數據集[1]，本研究FVC 值保留6 位有效數字，用戶在使用過程中可根據實際需要自行決定保留幾位有效數字。

數據作者分工職責

沈明潭，男，云南文山人，碩士研究生，研究方向為林業遙感技術與應用。主要承擔工作：數據處理，數據論文撰寫。

譚炳香，女，北京市人，博士，研究員，研究方向為森林資源遙感監測與森林信息遙感提取等。主要承擔工作：數據論文修改。

侯瑞霞，女，山西翼城人，博士，副研究員，研究方向為GIS 與數據分析。主要承擔工作：數據論文撰寫建議。

黃逸飛，女，湖南邵陽人，碩士研究生，研究方向為林業遙感技術與應用。主要承擔工作：基礎數據收集。