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

1. 中國(guó)科學(xué)院大學(xué)，北京 100049

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

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

數(shù)據(jù)庫(kù)（集）基本信息簡(jiǎn)介

引 言

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

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

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

1 數(shù)據(jù)采集和處理方法

1.1 數(shù)據(jù)來(lái)源與預(yù)處理

中巴經(jīng)濟(jì)走廊荒漠化分布數(shù)據(jù)的制備基于2種數(shù)據(jù)：MODIS植被指數(shù)產(chǎn)品MOD13A3數(shù)據(jù)和地表反照率產(chǎn)品MCD43A3數(shù)據(jù)，數(shù)據(jù)驗(yàn)證采用Landsat 7數(shù)據(jù)，詳細(xì)列表如表1所示。

表1 研究采用數(shù)據(jù)列表

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

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

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

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

NDVI和Albedo原始數(shù)據(jù)有異常值和像元缺失的情況，質(zhì)量不是非常可靠。需要進(jìn)行無(wú)效值剔除和空間插補(bǔ)。NDVI數(shù)據(jù)的有效范圍為-2000～10000，無(wú)效值用-3000填充。MCD43A3產(chǎn)品中所使用的1-6白空反照率波段的有效值范圍為0～1000，無(wú)效值用32766填充。本文采用反距離加權(quán)（IDW）方法實(shí)現(xiàn)數(shù)據(jù)的空間差值，IDW方法參數(shù)設(shè)置如下：

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

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

1.2 數(shù)據(jù)處理步驟

1.2.1 NDVI和Albedo數(shù)據(jù)年尺度重建

植被覆蓋度的變化是荒漠化最直觀的表現(xiàn)，需要在一年中植被最茂盛的情況下評(píng)價(jià)這一年的荒漠化程度。所以，基于年度NDVI和年度Albedo構(gòu)建的荒漠化差值指數(shù)，需要以NDVI的年度最大值，Albedo的年度最小值作為基礎(chǔ)數(shù)據(jù)。計(jì)算過(guò)程包括以下幾個(gè)方面：

（1）Albedo數(shù)據(jù)計(jì)算

MODIS窄波段反照率向?qū)挷ǘ畏凑章实霓D(zhuǎn)化使用的是Liang[8-9]研究的算法，選用短波反照率計(jì)算方法，公式如下：

其中αshort是短波反照率，α1-α6分別代表MCD43A3中1-6波段，本研究中選用了白空反照率來(lái)計(jì)算Albedo，因?yàn)榘卓辗凑章适歉鱾€(gè)入射角的積分，更接近一般意義上的地表反照率。

（2）NDVI年度數(shù)據(jù)計(jì)算

采用最大值合成法（Maximum Value Compositing，MVC）。該方法在遙感數(shù)據(jù)處理方面，主要用于在一定時(shí)間段內(nèi)像元的數(shù)據(jù)分析重構(gòu)。計(jì)算公式：

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

（3）Albedo年度數(shù)據(jù)計(jì)算

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

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

1.2.2 基于Albedo與NDVI構(gòu)建荒漠化指數(shù)

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

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

圖2 中巴經(jīng)濟(jì)走廊2010年Albedo-NDVI二維特征空間圖

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

圖3 NDVI-Albedo空間中荒漠化差值指數(shù)（DDI）圖形表達(dá)

在具體應(yīng)用中，常數(shù)α可根據(jù)Albedo-NDVI特征空間中斜率來(lái)確定，k為特征方程的斜率。

根據(jù)上面的計(jì)算結(jié)果，Albedo-NDVI的負(fù)相關(guān)關(guān)系斜率k=-0.2303，則α=4.3422，DDI的表達(dá)式如下：

1.2.3 荒漠化等級(jí)劃分

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

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

表2 DDI分級(jí)表

1.2.4 荒漠化指數(shù)專題數(shù)據(jù)模擬流程

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

圖4 荒漠化指數(shù)專題數(shù)據(jù)模擬流程

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

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

圖5 2010年中巴經(jīng)濟(jì)走廊荒漠化分級(jí)圖

3 數(shù)據(jù)質(zhì)量控制和評(píng)估

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

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

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

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

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

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

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

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

致 謝

感謝LPDAAC提供的MODIS數(shù)據(jù)，感謝USGS提供的Landsat 7數(shù)據(jù)。感謝國(guó)家特殊環(huán)境、特殊功能觀測(cè)研究臺(tái)站共享服務(wù)平臺(tái)給予的觀測(cè)和數(shù)據(jù)方面指導(dǎo)與支持。