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

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

1.中國林業(yè)科學(xué)研究院資源信息研究所，北京 100091

引 言

珠三角在國家整體發(fā)展中具有重要的戰(zhàn)略地位，自改革開放以來經(jīng)歷了快速的經(jīng)濟(jì)增長、人口增長和城市化擴(kuò)張，導(dǎo)致當(dāng)?shù)刂脖患眲⊥嘶推茐?，從而?yán)重影響了生態(tài)系統(tǒng)的服務(wù)和功能，因此動(dòng)態(tài)監(jiān)測珠三角植被的時(shí)空變化和分布格局對當(dāng)?shù)氐纳鷳B(tài)建設(shè)具有重要的參考意義。在已公開發(fā)布包含珠三角地區(qū)的植被覆蓋度數(shù)據(jù)集產(chǎn)品中，高吉喜等基于空間分辨率為250 m 的MODI3Q1 NDVI（Normalized differential vegetation index）數(shù)據(jù)，采用NDVI 最大合成法的方式進(jìn)行中國區(qū)域的植被覆蓋反演[1]。與該數(shù)據(jù)集相比，本研究選用空間分辨率更高的30 m Landsat5 TM 和Landsat8 OLI 數(shù)據(jù)，利用中值合成法獲取研究區(qū)的影像數(shù)據(jù)，更能精確估測出珠三角地區(qū)的植被覆蓋度，更科學(xué)有效地反映出珠三角地區(qū)的植被覆蓋狀況。

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

1 數(shù)據(jù)來源和處理方法

本研究采用以下主要步驟對珠三角植被覆蓋度進(jìn)行分析：（1）在GEE（Google Earth Engine）平臺(tái)調(diào)用Landsat 數(shù)據(jù)，基于NDVI 像元二分法估算珠三角2000、2005、2010、2015 和2020 年植被覆蓋度；（2）利用變異系數(shù)分析珠三角5 個(gè)時(shí)期植被覆蓋度變化的波動(dòng)性；（3）利用Slope趨勢分析來探究植被覆蓋度變化的時(shí)空變化特征；（4）采用Hurst 指數(shù)預(yù)測研究區(qū)植被覆蓋度未來的發(fā)展趨勢。

1.1 數(shù)據(jù)來源

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

圖1 珠三角Landsat8 標(biāo)準(zhǔn)假彩色影像Figure 1 Standard false colour images of Landsat8 in the Pearl River Delta

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

1.2 數(shù)據(jù)處理方法

1.2.1 植被覆蓋度估測

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

1.2.2 歸一化差值植被指數(shù)

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

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

1.2.3 像元二分模型

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

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

1.2.4 變異系數(shù)

對研究區(qū)5 個(gè)時(shí)期的植被覆蓋度數(shù)據(jù)逐像元計(jì)算變異系數(shù)，探究研究區(qū)植被覆蓋度的穩(wěn)定性[17]。CV 值越大，植被變化的波動(dòng)性越大，將波動(dòng)性劃分為4 種類型，見表2。計(jì)算公式如下所示：

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

式中CV植被覆蓋度的變異系數(shù)，n 是研究期數(shù)，F(xiàn)VCi第i 年的植被覆蓋度值，F(xiàn)???VC???是研究期內(nèi)植被覆蓋度均值。

1.2.5 Slope 趨勢分析

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

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

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

采用F 檢驗(yàn)對變化趨勢進(jìn)行顯著性檢驗(yàn)，當(dāng)檢驗(yàn)結(jié)果P<0.05 為顯著變化，0.050.1 為不顯著變化。統(tǒng)計(jì)量計(jì)算公式如下：

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

1.2.6 Hurst 指數(shù)

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

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

累計(jì)離差：

極差：

標(biāo)準(zhǔn)差：

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

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

2 數(shù)據(jù)樣本描述

2.1 數(shù)據(jù)組成

本數(shù)據(jù)集包括：（ 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 個(gè)時(shí)期的年度植被覆蓋度數(shù)據(jù)；（2）PRD_FVC_CV.tif 為植被覆蓋度變異系數(shù)數(shù)據(jù)；（3）PRD_FVC_Slope_F.tif 為植被覆蓋度變化趨勢數(shù)據(jù)；（4）PRD_FVC_Hurst.tif 為植被覆蓋度未來變化趨勢數(shù)據(jù)。上述數(shù)據(jù)空間分辨率為 30 m，地理坐標(biāo)系為 WGS1984，投影坐標(biāo)系為WGS_1984_UTM_Zone_49N，數(shù)據(jù)保存為GeoTIFF 格式，整個(gè)數(shù)據(jù)集解壓后大小約為2.75 GB。

2.2 樣本描述

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

圖2 珠三角5 期植被覆蓋度（2000、2005、2010、2015、2020）、變異系數(shù)(a)、變化趨勢(b)、Hurst 指數(shù)(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 數(shù)據(jù)質(zhì)量控制和評估

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

3.1 精度驗(yàn)證

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

3.2 一致性檢驗(yàn)

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

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

4 數(shù)據(jù)價(jià)值

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

5 數(shù)據(jù)使用方法和建議

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

數(shù)據(jù)作者分工職責(zé)

沈明潭，男，云南文山人，碩士研究生，研究方向?yàn)榱謽I(yè)遙感技術(shù)與應(yīng)用。主要承擔(dān)工作：數(shù)據(jù)處理，數(shù)據(jù)論文撰寫。

譚炳香，女，北京市人，博士，研究員，研究方向?yàn)樯仲Y源遙感監(jiān)測與森林信息遙感提取等。主要承擔(dān)工作：數(shù)據(jù)論文修改。

侯瑞霞，女，山西翼城人，博士，副研究員，研究方向?yàn)镚IS 與數(shù)據(jù)分析。主要承擔(dān)工作：數(shù)據(jù)論文撰寫建議。

黃逸飛，女，湖南邵陽人，碩士研究生，研究方向?yàn)榱謽I(yè)遙感技術(shù)與應(yīng)用。主要承擔(dān)工作：基礎(chǔ)數(shù)據(jù)收集。