基于MODIS影像喀斯特石漠化狀況研究
——以貴州省為例

吳林霖，官云蘭，李嘉偉，袁晨鑫，李 睿

(1.東華理工大學(xué)測繪工程學(xué)院，南昌 330013；2.流域生態(tài)與地理環(huán)境監(jiān)測國家測繪地理信息局重點實驗室，南昌 330013；3.廣東南方數(shù)碼科技股份有限公司，廣州 510665)

0 引言

喀斯特地貌在我國又稱為巖溶地貌，一般指碳酸鹽巖分布地區(qū)或存在流經(jīng)石灰?guī)r的地下水所形成的特有地貌現(xiàn)象[1]。貴州省是我國喀斯特地貌分布最廣的省份，比“喀斯特”名稱起源地的前南斯拉夫的喀斯特地貌分布面積占比還大，其石漠化問題也最為突出[2]。由袁道先院士首先提出了石漠化的概念，即含有植被、土地覆蓋的喀斯特地區(qū)經(jīng)各類因素使得巖石裸露的過程，并指出了石漠化是中國南方亞熱帶喀斯特地區(qū)嚴(yán)峻的生態(tài)問題[3]。我國對石漠化重視較早，取得了一定的治理進(jìn)展，研究方法也由最初的人工調(diào)查轉(zhuǎn)變?yōu)榛?S集成技術(shù)的石漠化時空演化與動態(tài)變化分析，以及石漠化信息提取和目標(biāo)識別。李水明等[4]基于K2結(jié)構(gòu)學(xué)習(xí)算法結(jié)合樣本分類信息獲取石漠化特征子集，并以此提取石漠化信息；藍(lán)安軍等[5]分析了植被覆蓋率與石漠化的負(fù)相關(guān)關(guān)系并引入石漠化動力指數(shù)來表征石漠化的分級和現(xiàn)狀；李松等[6]基于1992年、2001年2個時期的TM/ETM+影像，在進(jìn)行數(shù)據(jù)預(yù)處理后采用圖像差值法對石漠化進(jìn)行了變化檢測，并提取了變化信息。當(dāng)前對喀斯特石漠化的研究主要集中在石漠化分布狀況、動態(tài)監(jiān)測、預(yù)警分析、演化機(jī)制和生態(tài)治理等方面[7]，但大多是在小區(qū)域尺度范圍內(nèi)進(jìn)行的，僅依據(jù)遙感提取的石漠化信息做單因素的石漠化特征反演，獲得的成果有一定的局限性；并且由于選取的指標(biāo)各不相同，研究結(jié)果難以進(jìn)行相互間的分析、比較。

為此，本文利用MODIS影像大空間尺度的特點，綜合考慮人地因素對石漠化的影響，引入4類石漠化評價因子(植被覆蓋度(vegetation fractional cover，VFC)、基巖裸露度、坡度和人口密度(population density，PD))建立石漠化評價模型，研究貴州省近10 a間的石漠化分布狀況，提取石漠化信息并對石漠化進(jìn)行變化檢測，進(jìn)而分析石漠化變化特征與時空演化趨勢?？焖?、高效地了解大范圍的石漠化現(xiàn)狀，監(jiān)測石漠化變化情況是提供石漠化治理方向和改善生態(tài)環(huán)境的基礎(chǔ)，在石漠化研究中具有重大的意義。

1 研究區(qū)及數(shù)據(jù)源

1.1 研究區(qū)概況

貴州省地處長江和珠江兩大水系的分水嶺地區(qū)，云貴高原東部，氣候溫暖濕潤，屬亞熱帶濕潤季風(fēng)氣候[8]。全省總面積為176 167 km2，地勢西部最高，中部稍低，并分別向北部四川盆地、東部湖南低山丘陵和南部廣西盆地逐步過渡，東西呈三級階梯，南北有兩大斜坡，平均海拔為1 100 m，最大高差達(dá)2 763 m，形成了一個總體地勢較高，內(nèi)部分異較大，深受河流切割的亞熱帶喀斯特高原山區(qū)[9]。貴州省是一個碳酸鹽巖廣泛分布、喀斯特強(qiáng)烈發(fā)育的區(qū)域，碳酸鹽巖總厚度約17 000 m，占全省沉積巖總厚度的70%以上，其中純碳酸鹽巖厚度約為碳酸鹽巖總厚度的62%[10-11]。

1.2 數(shù)據(jù)源

本文使用的數(shù)據(jù)包括遙感影像、數(shù)字高程模型(digital elevation model，DEM)數(shù)據(jù)、人口數(shù)據(jù)及地理矢量數(shù)據(jù)等。遙感影像采用EOS系列衛(wèi)星上搭載的MODIS中分辨率成像光譜儀傳感器獲取的陸地標(biāo)準(zhǔn)產(chǎn)品MOD09A1 500M 地表反射率8 d合成產(chǎn)品，提供了1—7波段的500 m空間分辨率影像數(shù)據(jù)，參數(shù)如表1所示。選取2007年、2010年和2016年獲取時間相近的3期影像進(jìn)行后續(xù)處理。

表1 MODIS 500 m地表反射率8 d合成產(chǎn)品參數(shù)Tab.1 MODIS 500 m surface reflectance 8-days composition product

DEM數(shù)據(jù)為SRTM(shuttle Radar topography mission)90 m空間分辨率原始高程數(shù)據(jù)，由美國國家航空航天局和國防部國家圖像測繪局聯(lián)合測量。DEM數(shù)據(jù)與MODIS數(shù)據(jù)均來源于中國科學(xué)院計算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(http://www.gscloud.cn)。人口數(shù)據(jù)來源于貴州省統(tǒng)計局《貴州省統(tǒng)計年鑒》獲取的2007年、2010年和2016年各市(州)年末常住人口。地理矢量數(shù)據(jù)來源于國家基礎(chǔ)地理數(shù)據(jù)庫，包括省界和市界。

2 評價因子

目前并沒有確定的石漠化評價體系，而遙感石漠化監(jiān)測需要考慮遙感技術(shù)的特點和石漠化評價的應(yīng)用需求[12]?；鶐r裸露度和VFC是石漠化最直觀的地面表現(xiàn)特征，故為石漠化評價的關(guān)鍵指標(biāo)[13]。坡度會影響土壤侵蝕，是誘發(fā)石漠化的原因之一[14]。石漠化是自然因素和人為因素共同作用的結(jié)果，人類活動對自然環(huán)境有著不可磨滅的影響，人口劇增給土地帶來的壓力以及不合理的人類活動，會加速石漠化的進(jìn)程[15]。因此，本文選取VFC、基巖裸露度、坡度和PD作為評價因子，根據(jù)層次分析法(analytic hierarchy process,AHP)建立石漠化評價模型，提取石漠化信息。

2.1 VFC

VFC是指植被(包括葉、莖、枝)在地面的垂直投影面積占統(tǒng)計區(qū)總面積的百分比[16]。岳躍民等[17]已驗證植被指數(shù)對于提取石漠化信息的有效性，故本文選用歸一化差值植被指數(shù)(normalized difference vegetation index，NDVI)估算VFC，并以李苗苗[18]在像元二分法基礎(chǔ)上研究的VFC提取方法建立模型。由于巖溶地區(qū)土被分布淺薄且不連續(xù)，在提取VFC時可以暫不考慮土壤的影響[19]。本文假設(shè)研究區(qū)可近似取VFCmax=100%，VFCmin=0，則VFC的計算公式為

(1)

式中NDVImax和NDVImin分別表示研究區(qū)NDVI最大和最小值。由于噪聲的影響，NDVImax和NDVImin一般會取一定置信度范圍內(nèi)的最大和最小值。根據(jù)置信區(qū)間5%～95%，使用ENVI軟件查看不同時期NDVI直方圖，獲取不同時期NDVI最大和最小值如表2所示。

表2 NDVI最大和最小值Tab.2 Maximum and minimum of NDVI

將表中的NDVImax和NDVImin代入式(1)，可估算出研究區(qū)2007年、2010年和2016年的VFC。

2.2 基巖裸露度

基巖裸露度采用Rikimaru[20]提出的裸土指數(shù)(bare soil index，BSI)估算。由于裸土的光譜特征在中紅外波段反射率最高，BSI通過高反射率波段減去低反射率波段從而突出裸土信息，并且加入了藍(lán)光波段進(jìn)行歸一化，計算結(jié)果中裸土的亮度(DN值)最高，其次為建筑用地，易于區(qū)分出裸土、植被與水體。BSI公式為

(2)

式中r1，r2，r3和r6分別為MODIS數(shù)據(jù)紅光、近紅外、藍(lán)光和短波紅外波段的反射率值。

2.3 坡度

坡度是地面陡緩的程度，是衡量石漠化的重要因子[14]。采用地理空間數(shù)據(jù)云平臺DEM數(shù)據(jù)切割模型對DEM裁剪，再使用ENVI軟件對DEM拼接，最后生成坡度圖。坡度計算公式為

(3)

式中：Slope為坡度；fx和fy分別為水平方向和垂直方向高程變化率。

2.4 PD

PD是指單位面積土地上居住的人口數(shù)[15]。本研究的人口信息從貴州省統(tǒng)計年鑒獲得(2007年、2010年、2016年各市(州)年末常住人口)，通過ArcGIS軟件將信息添加入縣市級矢量數(shù)據(jù)屬性表中，再經(jīng)由核密度分析計算研究區(qū)各市級單位的人口分布情況，獲取研究區(qū)的PD。

2.5 評價因子歸一化

不同的評價因子數(shù)值大小不在同一量級，量綱不統(tǒng)一，直接對石漠化程度進(jìn)行評價較為困難。因此，需要對各評價因子進(jìn)行歸一化，使其具有相同的量綱進(jìn)而參與石漠化評價模型計算。將編碼值設(shè)為1—10，統(tǒng)一各評價因子的量值，以數(shù)值大小表示對石漠化的影響程度，影響越大，編碼值越大；影響越小，編碼值則越小。坡度、基巖裸露度和PD對石漠化起正向影響，數(shù)值從低到高分為1—10的等級。VFC值越高表示植被越茂密，對石漠化起反向影響；VFC值越高，石漠化程度越小，反之則越大。因此，將VFC值從高到低分為1—10的等級。

3 石漠化評價模型建立

石漠化評價模型是基于各類歸一化評價因子，通過比較評價因子重要性來綜合計算各評價因子權(quán)重而建立的，以Saaty[21]提出的AHP來確定評價因子權(quán)值。該方法是一種定性和定量相結(jié)合的目標(biāo)決策分析方法。首先，根據(jù)相對重要性的比例標(biāo)度評判各因子的重要性建立判斷矩陣，比例標(biāo)度如表3所示；其次，經(jīng)過一系列計算逐步剔除主觀元素，使計算的權(quán)值盡可能客觀公正；最后，所求出的權(quán)值還需要對其進(jìn)行一致性檢驗，計算出一致性檢驗系數(shù)CR，當(dāng)CR<0.1時，判斷矩陣通過一致性檢驗，否則不滿足一致性。

表3 比例標(biāo)度Tab.3 Proportion criteria

①8,6,4,2,1/2,1/4,1/6,1/8為以上重要性標(biāo)度的中間值。

具體計算步驟如下：

1)根據(jù)表3所示的比例標(biāo)度結(jié)合研究區(qū)實際情況將評價因子兩兩進(jìn)行比較確定相對重要性數(shù)值來建立判斷矩陣A，如表4所示。

表4 判斷矩陣及權(quán)值Tab.4 Judgment matrix and weight

2)計算判斷矩陣各行的乘積，即

(4)

式中：i和j分別為判斷矩陣的行和列號；n取值為4，因為選擇了4個評價因子建立模型。

3)計算Mi的n次方根，即

(5)

(6)

式中wi為第i行對應(yīng)的評價指標(biāo)權(quán)值。

5)將上述結(jié)果代入公式(7)求出最大特征根，即

(7)

式中(Aw)i為矩陣A第i行所有元素與向量w第i行對應(yīng)權(quán)值的乘積。

6)計算一致性系數(shù)CR，檢驗一致性。其表達(dá)式為

(8)

(9)

式中RI通過查表可得，當(dāng)n=4時，RI=0.89。當(dāng)CR<0.1時，通過一致性檢驗，否則修改判斷矩陣比例標(biāo)度，直到通過一致性檢驗。根據(jù)上述步驟計算出CR=0.060 726 704，通過一致性檢驗，說明判斷矩陣比例標(biāo)度合理，權(quán)值可行。根據(jù)計算的權(quán)值建立石漠化評價模型D，即

D=w1VFC+w2BSI+w3Slope+w4PD，

(10)

式中：w1=0.55，為VFC所對應(yīng)的權(quán)值；w2=0.28，為BSI所對應(yīng)的權(quán)值；w3=0.12，為坡度所對應(yīng)的權(quán)值；w4=0.05，為PD所對應(yīng)的權(quán)值。

4 石漠化狀況分析與評價

4.1 石漠化分級標(biāo)準(zhǔn)

對石漠化程度分級，基于不同的影像數(shù)據(jù)和評價因子有不同的分級標(biāo)準(zhǔn)。劉芳等[22]在Landsat8熱紅外影像基礎(chǔ)上，確定不同石漠化程度的亮溫差，并依據(jù)亮溫差建立了評價標(biāo)準(zhǔn)；李麗等[23]基于NDVI估算VFC，并依據(jù)研究區(qū)的VFC建立了石漠化評價標(biāo)準(zhǔn)。本文沿用現(xiàn)行的分級方法[24]，將石漠化程度劃分為4個等級，即無石漠化、輕度石漠化、中度石漠化和重度石漠化。根據(jù)歸一化編碼值(1—10)分為4類，如表5所示。

表5 石漠化分級標(biāo)準(zhǔn)Tab.5 Hierarchy standard of rocky desertification extent

4.2 3個時期石漠化狀況分析

根據(jù)上述方法及分級標(biāo)準(zhǔn)，得到2007年、2010年及2016年3期貴州省不同等級的石漠化空間分布狀況，如圖1所示。

(a)2007年 (b)2010年 (c)2016年

圖1 石漠化分級結(jié)果

Fig.1Classificationofrockydesertificationextent

1)2007年石漠化分布。2007年全省大部分區(qū)域為輕度石漠化，其次為中度石漠化。輕度石漠化占全省總面積的50.80%；中度石漠化與重度石漠化占全省總面積的比例分別為42.16%和6.99%，其中凱里市和遵義市南部存在較少的石漠化現(xiàn)象，中度石漠化主要分布在遵義市東北部、畢節(jié)市和六盤水市，興義市、貴陽市也有少量分布；重度石漠化主要分布在畢節(jié)市，少數(shù)分布于安順市與六盤水市的交界區(qū)、遵義市和貴陽市。

2)2010年石漠化分布。相較于2007年而言，全省范圍的重度石漠化區(qū)域有所減少，約占全省總面積的3.69%；中度石漠化約占全省面積的40.53%；而輕度石漠化則占全省面積55.67%。其中畢節(jié)市依舊表現(xiàn)為大面積的中度石漠化，重度石漠化相較2007年有部分減少；六盤水市和興義市西部中度、重度石漠化分布范圍增加；都勻市與遵義市主要為中度石漠化和極小片的重度石漠化；安順市與凱里市基本為輕度石漠化和無石漠化。

3)2016年石漠化分布。2016年的石漠化情況明顯改善了許多，輕度石漠化占全省面積的79.71%，中度石漠化和重度石漠化分別減少到了17.35%和2.87%。其中畢節(jié)市由重度石漠化大范圍分布轉(zhuǎn)變?yōu)檩p度石漠化，使輕度石漠化所占比重明顯增加，重度石漠化仍主要分布在西部地區(qū)；遵義市、銅仁市和都勻市的石漠化情況較為相似，有少許中度石漠化，而大范圍分布輕度石漠化；貴陽市重度石漠化分布沒有明顯變化，說明貴陽市重度石漠化程度較穩(wěn)定；興義市西北部中度石漠化分布較多；凱里市主要為輕度石漠化，中度石漠化有零散分布。

4.3 石漠化總體變化

通過圖像差值法，將對應(yīng)像元的灰度值相減，獲取不同時相的差異以達(dá)到對2景不同時期影像進(jìn)行石漠化變化檢測的目的。對2007年、2010年和2016年3個時期的石漠化評價結(jié)果兩兩進(jìn)行差值，將差值圖像閾值設(shè)為：-1為石漠化減少、1為石漠化增加，結(jié)果如圖2所示。

(a)2007—2010年 (b)2010—2016年

圖2 石漠化變化檢測

Fig.2ChangedetectionofKarstrockydesertification

由圖2(a)可以看出，2007—2010年間，石漠化減少區(qū)域明顯大于石漠化增加區(qū)域。在貴州省西部(六盤水市、興義市)，中部(貴陽市、都勻市)和東北部(遵義市、銅仁市)等有少量石漠化增加現(xiàn)象，而貴州西北部(畢節(jié)市)石漠化明顯減少。如圖2(b)所示，2010—2016年間，石漠化趨勢主要表現(xiàn)為減少，大部分減少區(qū)域位于貴州省西北部(畢節(jié)市、六盤水市)，東北部(遵義市、銅仁市)，中部、西南部有少許減少趨勢；增加區(qū)域主要在貴州省西南部(興義市、安順市和六盤水市)，與2007—2010年相比，呈現(xiàn)持續(xù)增加趨勢。

4.4 石漠化分級變化

通過比較分級區(qū)域的變化，可以檢測出輕度石漠化、中度石漠化和重度石漠化相互間的轉(zhuǎn)移變化情況，2007—2010年和2010—2016年期間石漠化轉(zhuǎn)移情況如圖3所示。

(a)2007—2010年 (b)2010—2016年

圖3 石漠化轉(zhuǎn)移情況

Fig.3Conversiondistributionofrockydesertification

由圖3可以看出，在2007—2010年間，輕度石漠化總體增加趨勢最明顯，由中度和重度石漠化轉(zhuǎn)移為輕度石漠化的比例分別為8.15%和0.25%，但仍有6.14%的輕度石漠化轉(zhuǎn)化為中度石漠化；重度石漠化轉(zhuǎn)為中度石漠化的比例僅有2.29%，表明這期間對石漠化的治理有一定成效，但還是存在石漠化程度加重的情況。在2010—2016年期間，石漠化現(xiàn)象得到明顯改善，其中中度石漠化向輕度石漠化的轉(zhuǎn)移比例為14.61%，轉(zhuǎn)換為重度石漠化的比例僅為0.58%；輕度石漠化轉(zhuǎn)化為中度石漠化的比例也只有2.59%。

對比2個階段的石漠化轉(zhuǎn)移情況，發(fā)現(xiàn)輕度轉(zhuǎn)向中度石漠化的比例遠(yuǎn)遠(yuǎn)大于輕度轉(zhuǎn)為重度石漠化的比例，重度石漠化轉(zhuǎn)為中度石漠化的比例大于轉(zhuǎn)為輕度石漠化的比例。該現(xiàn)象說明中度石漠化為重度石漠化與輕度石漠化相互轉(zhuǎn)化的中間過程，從石漠化治理角度可以先由重度石漠化改善為中度石漠化，再由中度石漠化改善為輕度石漠化；中度石漠化與輕度石漠化相互之間的轉(zhuǎn)化占比較大，說明這2種石漠化類型具有跳躍性；重度石漠化轉(zhuǎn)移面積比較少，說明重度石漠化具有穩(wěn)定性，難以改善。這與白曉永等[9]、趙麗蘋[25]發(fā)現(xiàn)石漠化轉(zhuǎn)移及演替的特征基本一致。

5 討論

在熊康寧等[26]的研究中先將研究區(qū)劃分為喀斯特區(qū)與非喀斯特區(qū)。本文考慮到貴州省碳酸鹽巖幾乎分布全境，石漠化區(qū)域也比較大，所以提取石漠化信息時沒有先剔除非喀斯特區(qū)，而是將整個研究區(qū)作為喀斯特區(qū)域。由于選擇的研究區(qū)范圍較大，選擇MODIS中分辨率成像光譜儀產(chǎn)品較為適合，并且選用MODIS地表反射率8 d合成產(chǎn)品，避免復(fù)雜的數(shù)據(jù)預(yù)處理過程。但由于MODIS影像的空間分辨率較低，不可避免的存在混合像元，加之研究區(qū)域范圍較大，且沒有實地考察來驗證本文研究的提取精度。因此，模型計算結(jié)果只能宏觀表示研究區(qū)石漠化概況與變化趨勢。評價因子計算中，李霞等[27]曾基于SPOT影像，采用歸一化土壤指數(shù)(normalized difference soil index，NDSI)和不透水面指數(shù)(normalized difference impervious surface index，NDISI)雙指數(shù)法確定一定經(jīng)驗閾值區(qū)分裸土與建筑用地。但是，MODIS影像尚不具備區(qū)分建筑用地與裸土的空間分辨率，故本文只采用BSI單一指數(shù)估算了基巖裸露度。在今后的研究中可以嘗試結(jié)合其他方法提高估算精度。

6 結(jié)論

基于MODIS地表反射率8 d合成產(chǎn)品，選取植被覆蓋度、基巖裸露度、坡度和人口密度作為評價因子，采用層次分析法計算各因子權(quán)值，通過一致性檢驗后建立石漠化評價模型，確定分級標(biāo)準(zhǔn)劃分石漠化程度，得到了不同時期的石漠化評價結(jié)果；對2007年、2010年和2016年3個時期的石漠化評價結(jié)果兩兩進(jìn)行變化檢測，獲得石漠化的總體變化以及石漠化分級變化狀況；完成了對貴州省石漠化狀況的評價，并分析了石漠化變化趨勢及時空演化特征。取得的主要結(jié)論如下：

1)研究區(qū)石漠化在近10 a間已得到了較大的改善，2010—2016年期間石漠化治理效果尤為突出。由2010年輕度石漠化與中度石漠化分別占全省總面積的55.67%和40.53%，到2016年輕度石漠化與中度石漠化分別占全省總面積的79.71%和17.35%，中度石漠化同比減少了23.18%，輕度石漠化增加了24.04%。重度石漠化有明顯的改善，從2007年占全省總面積的6.99%降到了2016年的2.87%。

2)在石漠化等級轉(zhuǎn)化過程中，中度石漠化與輕度石漠化相互之間的轉(zhuǎn)化率較高，說明這2種石漠化類型具有跳躍性；中度石漠化為重度石漠化與輕度石漠化相互轉(zhuǎn)化的中間過程，從石漠化治理角度可以先由重度石漠化改善為中度石漠化，再由中度石漠化改善為輕度石漠化。

3)對比2007—2010年和2010—2016年期間的石漠化轉(zhuǎn)移矩陣，重度石漠化轉(zhuǎn)移面積最小，表明重度石漠化較為穩(wěn)定，改善難度較大。2007—2010年間由重度石漠化轉(zhuǎn)化為中度石漠化和輕度石漠化所占比例分別為0.25%和2.29%，2010—2016年間由重度石漠化轉(zhuǎn)化為中度石漠化和輕度石漠化所占比例分別為0.59%和0.45%。