2000-2017年中巴經濟走廊逐年荒漠化分布數據集

敏玉芳，馮克庭，康建芳，艾鳴浩

1. 中國科學院大學，北京 100049

2. 中國科學院西北生態(tài)環(huán)境資源研究院, 科學大數據中心，蘭州 730000

3. 國家特殊環(huán)境、特殊功能觀測研究臺站共享服務平臺, 平臺服務中心，蘭州730000

數據庫（集）基本信息簡介

引 言

荒漠化包括氣候變化和人類活動在內的種種因素造成的干旱、半干旱和干旱亞濕潤地區(qū)的土地退化，屬于破壞性的生物地理過程，通常會造成生物多樣性降低，土壤肥力下降乃至生態(tài)承載力喪失[1-2]。土地荒漠化是目前世界上最為嚴重的生態(tài)環(huán)境問題之一，據統(tǒng)計已有100多個國家正受到土地荒漠化的影響，荒漠化不僅破壞生態(tài)環(huán)境，也會削弱社會和經濟的發(fā)展。監(jiān)測土地荒漠化動態(tài)變化，掌握其變化規(guī)律對防治荒漠化有很重要的意義。

中巴經濟走廊起點在中國新疆喀什，終點在巴基斯坦瓜達爾港，北接“絲綢之路經濟帶”、南連“21世紀海上絲綢之路”，是貫通南北絲路的關鍵樞紐，是一條包括公路、鐵路、油氣和光纜通道在內的貿易走廊，也是“一帶一路”的重要組成部分。中巴經濟走廊全長約3000公里，其中1300余公里為嚴寒高原、干旱和荒漠區(qū)段。尤其是在南段，干旱和大面積荒漠是其主要的生態(tài)環(huán)境約束因素。

遙感技術為人們提供了一種全新的荒漠化監(jiān)測手段，它具有觀測范圍廣、信息量大、數據更新快和精度高的特點，通過遙感圖像解譯或者定量反演都可以準確及時地獲取荒漠化土地變化的信息?；哪谶b感影像上表現為裸地地表信息的增強和植被信息的減弱，可以采用地表反照率、地表溫度、地表濕度、植被指數、植被覆蓋度等指標因子表征。歸一化植被指數（NDVI）是反映地表植被狀態(tài)的重要生物物理參數，而地表反照率（Albedo）則是反映地表對太陽短波輻射反射特性的物理參量。隨著荒漠化程度的加重，地表植被遭受嚴重破壞，地表植被蓋度降低，生物量減少，地表粗糙度增加，在遙感圖像上表現為NDVI值相應減少，地表反照率得到相應的增加。Li[3]等經過定位觀測研究表明，當地表反照率達到一定數值時，會發(fā)生草地荒漠化，荒漠化發(fā)生的地表反照率閾值為0.3。Albedo 與 NDVI 這兩個指標因子分別通過表征裸露地表信息、植被覆蓋度信息，與荒漠化過程建立聯(lián)系，具有較好的相關性。曾永年等[4]提出不單獨依靠某一個參數，而是通過NDVI和Albedo的負相關關系構造“NDVI-Albedo”特征空間，并計算出荒漠化差值指數（DDI），利用DDI來進行荒漠化分級。該方法使用簡單，易于獲取指標，在荒漠化的分類及分級中，較僅使用遙感光譜信息進行分類的方法精度更高，所以在近年來不同地區(qū)的荒漠化定量評價中得到了一定的應用[4-7]。本文采用DDI評價中巴經濟走廊荒漠化程度，研究區(qū)如圖1所示。以NDVI與Albedo為監(jiān)測指標，通過構造Albedo-NDVI特征空間，并利用Albedo和NDVI之間負相關的關系，構建中巴經濟走廊DDI公式，并完成2000-2017年共18期中巴經濟走廊荒漠化分類專題數據集，為相關部門的治理和決策提供技術支持。

1 數據采集和處理方法

1.1 數據來源與預處理

中巴經濟走廊荒漠化分布數據的制備基于2種數據：MODIS植被指數產品MOD13A3數據和地表反照率產品MCD43A3數據，數據驗證采用Landsat 7數據，詳細列表如表1所示。

表1 研究采用數據列表

Albedo是遙感反演中的重要參數，指單位時間、單位面積內地表全部反射輻射通量與入射太陽總輻射通量之比。隨著荒漠化程度增加，地表形態(tài)發(fā)生顯著變化，集中表現在地表植被覆蓋率下降，土壤有機質含量降低，土壤水分減少，地表粗糙度增加，反照率上升。MCD43A3是綜合了Terra和Aqua衛(wèi)星的數據產品，為天尺度L3級產品，時空分辨率是day/500 m。MCD43A3產品數據包含白空反照率和局地正午太陽角的黑空反照率，且白空反照率和黑空反照率都包含1-7波段共7個窄波段，以及3個寬波段可見光（0.3～0.7 μm）、近紅外（0.7～5.0 μm）、短波（0.3～5.0 μm）。數據來源于美國國家航空航天局（NASA）陸地過程分布式數據檔案中心（LPDAAC）。數據以HDF-EOS格式存儲，采用正弦曲線投影。覆蓋整個中巴走廊需要6景6個波段的數據。本研究獲取2000-2017年18年MCD43A3數據共39 090景。首先使用MRT工具對數據進行鑲嵌和重采樣，將500 m分辨率的數據重采樣為 1 km，并將所有數據正弦曲線投影轉為 WGS84。然后利用中巴走廊矢量邊界和Gdal的warp工具進行數據裁剪。

圖1 研究區(qū)范圍地理位置示意圖

NDVI是反應植物量和植被檢測的重要指標，主要應用于監(jiān)測植被生長狀態(tài)、植被覆蓋度和消除部分輻射誤差等。本文中使用1km的MODIS月合成植被指數產品MOD13A3。MOD13A3數據包括1 km月合成的NDVI、增強植被指數EVI、NDVI質量文件、EVI質量文件、可見光紅光波段反射率、近紅外波段反射率、可見光藍光波段反射率、中紅外波段反射率、平均觀測天頂角、平均太陽天頂角、平均方位角。數據來源于LPDAAC。數據以HDF-EOS格式存儲，采用Sinusoidal投影。本研究獲取2000-2017年18年MOD13A3L3數據共1296景。首先使用MRT工具對數據進行鑲嵌和重采樣，將所有數據Sinusoidal坐標轉為WGS84。然后利用中巴走廊矢量邊界和Gdal的warp工具進行數據裁剪，最后得到216景研究區(qū)NDVI數據。

本文中還使用Landsat7 ETM+數據用于數據質量評估。Landsat 7 ETM+影像包括8個波段，band 1-band 5和band 7的空間分辨率為30 m，band 6的空間分辨率為60 m，band 8的空間分辨率為15 m，時間分辨率為16天。本研究中獲取2010年7-8月，云量小于5%的Landsat 7影像8景。利用ENVI軟件對數據進行幾何校正和大氣校正，NDVI計算與統(tǒng)計，最后利用像元二分模型進行植被覆蓋度的遙感估算，得到研究區(qū)的植被覆蓋度數據。

NDVI和Albedo原始數據有異常值和像元缺失的情況，質量不是非?？煽?。需要進行無效值剔除和空間插補。NDVI數據的有效范圍為-2000～10000，無效值用-3000填充。MCD43A3產品中所使用的1-6白空反照率波段的有效值范圍為0～1000，無效值用32766填充。本文采用反距離加權（IDW）方法實現數據的空間差值，IDW方法參數設置如下：

（1）距離指數：通常取值范圍在0.5～3之間，本文選擇值為2；

（2）搜索半徑：定義對缺失像元值進行插值的輸入點，選擇可變搜索半徑方式，插值的最鄰近輸入采樣點數量為12。

1.2 數據處理步驟

1.2.1 NDVI和Albedo數據年尺度重建

植被覆蓋度的變化是荒漠化最直觀的表現，需要在一年中植被最茂盛的情況下評價這一年的荒漠化程度。所以，基于年度NDVI和年度Albedo構建的荒漠化差值指數，需要以NDVI的年度最大值，Albedo的年度最小值作為基礎數據。計算過程包括以下幾個方面：

（1）Albedo數據計算

MODIS窄波段反照率向寬波段反照率的轉化使用的是Liang[8-9]研究的算法，選用短波反照率計算方法，公式如下：

其中αshort是短波反照率，α1-α6分別代表MCD43A3中1-6波段，本研究中選用了白空反照率來計算Albedo，因為白空反照率是各個入射角的積分，更接近一般意義上的地表反照率。

（2）NDVI年度數據計算

采用最大值合成法（Maximum Value Compositing，MVC）。該方法在遙感數據處理方面，主要用于在一定時間段內像元的數據分析重構。計算公式：

其中，i為像元名，j為[1,n]時間區(qū)間的時間點，NDVIij指第i個像元在第j時間的NDVI值。

（3）Albedo年度數據計算

Albedo數據的年度重建采用最小值合成法。該方法是在某一指定的時間段內，選取像元的最小值作為新的像元值。年度Albedo的計算公式為：

其中，i為像元名，j為[1,n]時間區(qū)間的時間點，Albedoij指第i個像元在第j時間的Albedo值。

1.2.2 基于Albedo與NDVI構建荒漠化指數

NDVI值與植被覆蓋度成正比，而Albedo與植被覆蓋度成反比。根據NDVI與植被覆蓋度的正相關關系，以及NDVI與Albedo的負相關關系可以得到Albedo-NDVI二維空間特征圖。本文以2010年數據為例，先對NDVI和Albedo數據進行歸一化處理，歸一化處理公式如下：

然后在研究區(qū)域均勻選擇1500個樣本點，以NDVI作為X軸，Albedo作為Y軸，構建NDVIAlbedo特征空間。特征空間散點圖如圖2所示，根據線性統(tǒng)計回歸分析結果，Albedo和NDVI間的負相關關系可以用以下線性關系式表示：

圖2 中巴經濟走廊2010年Albedo-NDVI二維特征空間圖

如圖2所示，由Albedo和NDVI構造的特征空間中，Albedo與NDVI的負相關關系式的不同位置，代表荒漠化不同階段的狀態(tài)和程度，荒漠化程度隨著NDVI值的減少而增加，隨著Albedo值的增加而增加，即Albedo與NDVI的負相關線性表達式可以反映荒漠化的變化趨勢。根據Verstraete和 Pinty[10]的研究結論，如果在代表荒漠化變化趨勢的垂直方向上劃分 Albedo-NDVI特征空間，可以將不同的荒漠化土地有效區(qū)分開來，如圖3所示，用荒漠化差值指數模型DDI來表示：

圖3 NDVI-Albedo空間中荒漠化差值指數（DDI）圖形表達

在具體應用中，常數α可根據Albedo-NDVI特征空間中斜率來確定，k為特征方程的斜率。

根據上面的計算結果，Albedo-NDVI的負相關關系斜率k=-0.2303，則α=4.3422，DDI的表達式如下：

1.2.3 荒漠化等級劃分

為荒漠化評價與制圖的需要，1984年聯(lián)合國糧農組織（FAO）和聯(lián)合國黃金規(guī)劃署（UNEP）在《荒漠化評價與制圖方案》中，從植被退化、風蝕、水蝕、鹽堿化等4個方面，提出了荒漠化現狀、發(fā)展速率、內在危險性評價的具體定量植被，并把荒漠化按其發(fā)展程度的不同，分為輕度、中度、重度、極重度。

通過上一節(jié)構建的 DDI模型，得到 2000-2017年的中巴經濟走廊荒漠化指數專題數據?；贒DI進行荒漠化等級劃分時，大部分研究人員都采用 ArcGIS重分類中的自然斷裂法[5,7,11]。自然斷裂法是基于統(tǒng)計學的Jank最優(yōu)法得出的分級點[12]，這種分類法使各類別中差異最小，類之間差異最大。本文也采用自然斷裂法和荒漠化程度四分法，將DDI值劃分為5個區(qū)間，確定5級分級指標（表2），用戶也可以結合實地調查數據，進一步精細化調整DDI分級表。

表2 DDI分級表

1.2.4 荒漠化指數專題數據模擬流程

首先重建NDVI和Albedo年數據，用最大值合成法計算年度NDVI值，用最小值合成法計算年度Albedo值；然后通過擬合Albedo和NDVI數據，構建Albedo-NDVI特征空間，得到擬合斜率k；最后構建荒漠化指數，進行荒漠化等級劃分，完成基于荒漠化指數的荒漠化分級圖。數據處理流程如圖4所示。

圖4 荒漠化指數專題數據模擬流程

2 數據樣本描述

中巴經濟走廊逐年荒漠化分布數據集有2000-2017年18個文件，每個文件夾包括年度NDVI、Albedo、DDI柵格數據和2010年荒漠化分級圖示例。全部數據的空間分辨率均是1 km，年度NDVI、年度Albedo、DDI柵格數據保存格式為tif，荒漠化分級圖采用jpg文件格式保存，用戶可以根據DDI數據制作荒漠化分級數據。所有數據地理坐標系為WGS1984。數據結果展示如圖5。

圖5 2010年中巴經濟走廊荒漠化分級圖

3 數據質量控制和評估

本文采用高分辨率數據來評價荒漠化分級數據的質量。以2010年數據為例，選取8景Landsat 7 影像數據在小范圍上進行驗證。8景影像空間分辨率為30 m，平均云量小于5%，數據質量良好。DDI是通過NDVI年度最大值和Albedo年度最小值計算而來，所以本文用處于植被生長旺季的8月份Landsat數據作為驗證數據，反演年度植被覆蓋度最大值。運用ENVI軟件對數據進行預處理，主要包括配準、幾何校正、圖像增強和FLAASH大氣校正。首先計算8景影像的NDVI值，然后進行植被蓋度的估算。植被覆蓋度計算公式如下：

式中：vfc 為研究區(qū)的植被覆蓋度，NDVIsoil為本研究區(qū)域純沙漠的NDVI值，NDVIveg為完全植被覆蓋區(qū)域的NDVI值。一般以5%置信度截取NDVI的上下閾值分別近似代表NDVIveg和NDVIsoil進行計算。

本文參考其他研究工作，將中巴經濟走廊地區(qū)的荒漠化程度按植被覆蓋度大概劃分為：非荒漠化（高植被覆蓋，vfc>60%）、輕度荒漠化（中等植被蓋度，30%～60%）、中度荒漠化（中低植被覆蓋，15%～30%）、重度荒漠化（低植被覆蓋，<15%）。驗證時，在每個影像上選取各級別荒漠化驗證點30個，由于Landsat跟Modis數據分辨率不一致，所以在選取驗證點時，在Landsat上選取小范圍且植被覆蓋度差別不大的區(qū)域作為一個驗證點，與DDI分級結果進行比較（表3），并計算Kappa系數。表4是基于此方法驗證的數據評價精度表，總體評價精度達到81.67%，Kappa系數為75.42%，兩種分辨率的數據評價結果基本一致，反映出本文的數據在區(qū)域尺度上評價中巴經濟走廊荒漠化具有可行性。

表4 單一類別精度及總體精度分析

4 數據價值

雖然荒漠化的影響因素復雜，使用目的多樣化，造成指標因子繁雜，但所有的指標最終都可以歸結到最基本的影響因子——植被。植被既是土地覆蓋質量的綜合體現，也是影響生態(tài)環(huán)境中水土保持等方面的主體因素。本文基于易獲取的遙感影像數據，使用反映植被信息的NDVI與反映土壤、水分等信息的Albedo，構建中巴經濟走廊的DDI模型，直觀反映中巴經濟走廊荒漠化程度。為定量評價荒漠化嚴重程度提供參考，為中巴經濟走廊沿線國家進行荒漠化防治工作以及宏觀決策的制定提供基礎資料。

5 數據使用方法和建議

2000-2017年中巴經濟走廊荒漠化分布相關數據保存為柵格tif格式。ArcGIS、QGIS、ENVI、ERDAS等常用的GIS與遙感軟件均支持該數據的讀取和操作。

致 謝

感謝LPDAAC提供的MODIS數據，感謝USGS提供的Landsat 7數據。感謝國家特殊環(huán)境、特殊功能觀測研究臺站共享服務平臺給予的觀測和數據方面指導與支持。