馬友平

(湖北民族大學(xué) 林學(xué)園藝學(xué)院，湖北 恩施 445000)

石漠化是一種石質(zhì)荒漠化，有地球癌癥之稱[1]，是土壤受到侵蝕后的直接后果與巖溶地區(qū)生態(tài)系統(tǒng)逆向演替的頂級(jí)階段，一旦發(fā)生，恢復(fù)治理相當(dāng)困難，導(dǎo)致生態(tài)環(huán)境的崩潰，且恢復(fù)速率極慢，已成為我國三大生態(tài)自然災(zāi)害之一[2-4].在我國主要分布在以云貴高原為中心的西南八省、市、自治區(qū)的451個(gè)縣(市)中，這樣大范圍內(nèi)有關(guān)石漠化基礎(chǔ)數(shù)據(jù)的獲取一直是石漠化研究的核心問題之一，其中遙感影像已成為石漠化研究的主要數(shù)據(jù)源，遙感影像主要依據(jù)的是地物的波譜特性，易受成像時(shí)間、地點(diǎn)、地物類別等的影響較大，現(xiàn)還沒有一個(gè)大范圍的統(tǒng)一標(biāo)準(zhǔn)，為此作者擬在該方面進(jìn)行一些探討，豐富以遙感影像為主要數(shù)據(jù)源的石漠化研究.

圖1 恩施市位置圖Fig.1 Enshi city position map

在石漠化的研究中,多以特定的行政區(qū)劃范圍為單位研究石漠化的發(fā)生機(jī)理、評(píng)定方法和繪制石漠化的等級(jí)分布圖.如周忠發(fā)等[5]以貴州省清鎮(zhèn)市為例，應(yīng)用TM遙感影像制作了清鎮(zhèn)市石漠化等級(jí)分布圖;周濤[6]以貴州省關(guān)嶺縣為例，以ETM+、TM遙感影像為數(shù)據(jù)源，采用CA_MarkoV模型研究了2007年喀斯特石漠化分布情況，并對(duì)模型模擬情況進(jìn)行評(píng)價(jià);胡寶清等[7]以廣西都安瑤族自治縣為例，應(yīng)用TM影像，建立數(shù)據(jù)庫，分析影響因子，認(rèn)為石漠化產(chǎn)生的主導(dǎo)原因是人類不合理的經(jīng)營活動(dòng)；朱林富等[8]應(yīng)用TM影像和ArcGIS軟件平臺(tái)對(duì)重慶巫山縣喀斯特石漠化與土地覆被關(guān)系進(jìn)行了研究；這些都是以縣域?yàn)閱挝?然而還有以省域或更大區(qū)域?yàn)閷?duì)象的.如白曉永等[9]用三個(gè)時(shí)期TM影像對(duì)貴州省的石漠化演變過程進(jìn)行了分析和評(píng)價(jià)；陳起偉等[10]以TM、ASTER和ALOS 三期衛(wèi)星影像為遙感動(dòng)態(tài)監(jiān)測數(shù)據(jù)源對(duì)貴州省的石漠化進(jìn)行了遙感監(jiān)測；靖娟利等[11]利用MODIS遙感影像對(duì)廣西自治區(qū)的喀斯特石漠化演變特征進(jìn)行了研究；童立強(qiáng)等[12]應(yīng)用ETM、TM、CBERS-02B、ALOS和ASTER衛(wèi)星遙感數(shù)據(jù)對(duì)中國南方地區(qū)石漠化進(jìn)行了遙感調(diào)查與演變分析.上述研究其基礎(chǔ)數(shù)據(jù)的獲取主要來自遙感影像，以遙感影像的光譜特性為基礎(chǔ)，以不同波段之間的線性與非線性組合來提取能反映石漠化狀況的指標(biāo).為此本研究擬在借鑒上述成果時(shí)以湖北省恩施市為例，討論石漠化范圍的確定和等級(jí)的劃分標(biāo)準(zhǔn)；對(duì)于指標(biāo)提取時(shí)影像的“缺失”區(qū)通過插值法獲得，有研究證明遙感影像各點(diǎn)的灰度值(DN)，并不相互獨(dú)立，它們之間存在一定的相關(guān)性[13-14]，為“缺失”區(qū)石漠化信息的地統(tǒng)計(jì)學(xué)估計(jì)提供了可能.

1 研究地區(qū)概況及數(shù)據(jù)來源

1.1 研究區(qū)概況

恩施市位于湖北省的西南部,如圖1所示，是恩施土家族苗族自治州的首府，東經(jīng)109°4′48″～109°58′42″、北緯29°50′33″～30°39′30″，面積397 2 km2,有世界硒都之美稱，小氣候復(fù)雜多變，石灰?guī)r分布廣泛，是南方巖溶發(fā)育區(qū)的一部分；雖經(jīng)過石漠化的專項(xiàng)治理，石漠化現(xiàn)狀有所改善，但問題依然存在，治理任務(wù)依然十分艱巨.

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

美國Landsat系列數(shù)據(jù)現(xiàn)已免費(fèi)開放，其中1999-04-15發(fā)射了Landsat-7，有8個(gè)波段，在波段設(shè)置上與Landsat-5非常相似，不同的是波段數(shù)上增加了一個(gè)全色波段，空間分辨率為15 m，在第6波段熱紅外波段上空間分辨率從120 m提高到了60 m.但在2003-05-31發(fā)生了故障，導(dǎo)致了一景影像中約有25%的數(shù)據(jù)丟失，因此這個(gè)時(shí)段以后的Landsat-7遙感數(shù)據(jù)稱之為Landsat-7 ETM SLC-OFF數(shù)據(jù)，該研究主要采用有條帶缺失的Landsat-7 ETM SLC-OFF影像，時(shí)間2015-10-22，軌道號(hào)/行號(hào)為126/39，空間分辨率30 m，云量0，而通過實(shí)際的圖片觀察恩施市上空有少量的云層；全市的DEM數(shù)據(jù)采用SRTM1(Shuttle Radar Topography Mission)，對(duì)應(yīng)的分辨率精度為30 m，來源地理空間數(shù)據(jù)云(網(wǎng)址：http://www.gscloud.cn/).同時(shí)也收集一些有關(guān)恩施市的社會(huì)、經(jīng)濟(jì)、文化資料，主要為《恩施市統(tǒng)計(jì)年鑒》、《恩施州統(tǒng)計(jì)年鑒》.收集到恩施市的行政邊界和區(qū)劃圖后，利用GPS采集地面控制點(diǎn)坐標(biāo)進(jìn)行幾何校正，校正后坐標(biāo)系均為WGS_1984_UTM_Zone_49N，與遙感影像、DEM為同一坐標(biāo)系.

2 研究方法

對(duì)于大區(qū)域石漠化的研究，遙感影像作為基礎(chǔ)數(shù)據(jù)源是許多專家學(xué)者的選擇之一，該研究選擇了Landsat-7 ETM SLC-OFF數(shù)據(jù)源，探討在有云、陰影、條帶缺失區(qū)的情形下與石漠化研究相關(guān)的數(shù)據(jù)修復(fù)問題，以彌補(bǔ)在應(yīng)用過程中因云量、時(shí)段等導(dǎo)致的數(shù)據(jù)源缺乏的問題，下面同時(shí)結(jié)合其他數(shù)據(jù)源闡明遙感影像在石漠化應(yīng)用中的研究方法.

2.1 Landsat-7 ETM SLC-OFF數(shù)據(jù)處理

按行政邊界進(jìn)行裁剪，分別將第1波段B1(Bands 1)、第4波段B4(Bands 4)、第7波段B7(Bands 7)放入紅綠藍(lán)(RGB)通道進(jìn)行假彩色合成，提取影像數(shù)據(jù)的缺失區(qū)域，灰度級(jí)為1bit，正常區(qū)為1，缺失區(qū)數(shù)據(jù)為0，為機(jī)載掃描行校正器SLC(Scan Lines Corrector) 突然發(fā)生故障，導(dǎo)致獲取的數(shù)據(jù)丟失部分；云覆蓋的位置彩用梯度法[15-16]確定；陰影采用圖像RGB中的G通道，HIS變換后的I通道進(jìn)行閾值標(biāo)定后，進(jìn)行邏輯與運(yùn)算所得[17]；在獲取了圖像的掉線、陰影、云的位置后進(jìn)行彭脹一個(gè)象元即30 m的處理，如圖2所示.圖2(a)為第一、四、七波段合成的原圖，圖2(b)為去掉不正常區(qū)域后的一、四、七波段合成圖形，經(jīng)統(tǒng)計(jì) “缺失”區(qū)面積1 503.829 0 km2，占整個(gè)國土面積的37.860 8%.

(a) B147原圖 (b) B147“缺失”圖圖2 原圖與“缺失”圖Fig.2 The original image and the “missing” map

2.2 基于遙感影像的石漠化指數(shù)與基巖裸露率

石漠化信息的提取精度與方法，直接關(guān)系到研究區(qū)域石漠化的成果，經(jīng)過查閱文獻(xiàn)[18-20]、地面實(shí)際調(diào)查和比較，采用增強(qiáng)型植被指數(shù)法建立研究區(qū)域的石漠化指數(shù)模型，公式如下：

(1)

式(1)中：Di為第i像元的石漠化指數(shù)值，DNiB5、DNiB4分別為第i個(gè)像元在第5波段、第4波段中的灰度值，GB4、GB5分別為遙感影像第4、第5波段所有像元灰度值的平均值.

由式(1)提取了石漠化指數(shù)后，還要轉(zhuǎn)化為基巖裸露率，基巖裸露率是石漠化等級(jí)劃分的一個(gè)非常重的條件指標(biāo)，計(jì)算公式為[20-21]：

(2)

式(2)中:Dgi為研究區(qū)域中第i個(gè)像元的基巖裸露率，Dmax、Dmin分別研究區(qū)域中石漠化指數(shù)的最大值與最小值，該研究中不包含“缺失”區(qū)的數(shù)據(jù)值.

2.3 石漠化等級(jí)的劃分

石漠化主要分布在我國廣大的南方地區(qū)，面積大，地域差異明顯，由于研究方法和數(shù)據(jù)的差異，現(xiàn)在還沒有一個(gè)統(tǒng)一的劃分標(biāo)準(zhǔn)，作者在參考前輩的研究成果和在實(shí)際調(diào)查的基礎(chǔ)上[12-13，22]，石漠化等級(jí)的劃分標(biāo)準(zhǔn)為二項(xiàng)指標(biāo)，即基巖裸露率和坡度；等級(jí)劃分為四級(jí)，分別為無石漠化(I)、輕度石漠化(II)、中度石漠化(III)和重度石漠化(IV)，詳見表1.

表1 巖溶區(qū)石漠化分級(jí)標(biāo)準(zhǔn)Tab.1 Classification standard of karst rock desertification

圖3 圖2(b)石漠化指數(shù)分布圖Fig.3 Distribution map of karst rocky desertification index figure 2(b)

2.4 克里格法(Kriging)

克里格法(Kriging)又稱為空間局部插值法，是在有限區(qū)域內(nèi)對(duì)區(qū)域化變量的取值進(jìn)行線性無偏最優(yōu)估計(jì)的一種方法，是地統(tǒng)計(jì)學(xué)的重要研究內(nèi)容之一，以變異函數(shù)來研究空間的結(jié)構(gòu)性與隨機(jī)性[14，23-24]，不同于經(jīng)典統(tǒng)計(jì)學(xué)考慮了變量間的自相關(guān)性；而遙感影像為柵格數(shù)據(jù)結(jié)構(gòu)，像元在圖像上呈規(guī)則網(wǎng)狀排列，特征值(石漠化指數(shù))的提取只考慮了地物的灰度值，而灰度值存在一定的空間相關(guān)性，顯然應(yīng)用統(tǒng)計(jì)學(xué)中的克里格內(nèi)插法可對(duì)“缺失”區(qū)域的數(shù)據(jù)進(jìn)行預(yù)估，提高精度[25-26].

石漠化指數(shù)作為一項(xiàng)石漠化狀況評(píng)價(jià)的基礎(chǔ)指標(biāo)，現(xiàn)還沒有統(tǒng)一的評(píng)定方法和嚴(yán)格的定義，在大尺度研究時(shí)，多以遙感影像為數(shù)據(jù)源，通過光譜特征來計(jì)算，對(duì)于“缺失”區(qū)預(yù)測數(shù)據(jù)，無法同地形因子如高程等，與地面實(shí)測值相較，評(píng)價(jià)擬合結(jié)果，僅能將參與擬合的數(shù)據(jù)按隨機(jī)的方法抽選一部分作為檢驗(yàn)點(diǎn)，該研究中抽取5%的數(shù)據(jù)，而剩余95%的數(shù)據(jù)參與克里格法預(yù)測建模.建模精度擬采用三個(gè)指標(biāo)評(píng)價(jià)：①?zèng)Q定系數(shù)(r)，用以反映觀測值與估計(jì)值的關(guān)系緊密程度，值越大估測精度越高；②平均絕對(duì)誤差(MAE)，反映估測值與檢驗(yàn)值的偏離大小，其值越接近0效果越好；③均方根誤差(RMSE)，介于0與1之間，取值越接近0，估測越精確.MAE和RMSE的計(jì)算公式如下:

(3)

(4)

式中：O實(shí)測值，P預(yù)測值，n為檢驗(yàn)數(shù)據(jù)的個(gè)數(shù).

3 結(jié)果與分析

3.1 石漠化指數(shù)的統(tǒng)計(jì)特征

石漠化指數(shù)的提取以ERDAS軟件中的建模能功為工具，以圖2(b)為遙感數(shù)據(jù)源，由式(1)可得圖3.

圖3中白色顏表示數(shù)據(jù)“缺失”區(qū)，數(shù)據(jù)值為NoData，將其轉(zhuǎn)換為點(diǎn)狀的矢量數(shù)據(jù)結(jié)構(gòu)，屬性值為石漠化指數(shù)，對(duì)其進(jìn)行是否服從正態(tài)布檢驗(yàn).經(jīng)試驗(yàn)和檢驗(yàn)原始數(shù)據(jù)做對(duì)數(shù)變換后服從正態(tài)分布，滿足克里格分析中簡單克里格、普通克里格等要求數(shù)據(jù)服從正態(tài)分布的要求[27]，否則可能會(huì)使變異函數(shù)產(chǎn)生比例效應(yīng).

3.2 石漠化指數(shù)“缺失”的克里格估計(jì)

在使用遙感影像作為數(shù)據(jù)源時(shí)，由于地形起伏、成像時(shí)間、大氣狀況和成像傳感器自身等因素，導(dǎo)致所獲影像部分不是地物光譜特征的真實(shí)反映，有的甚至為數(shù)據(jù)缺失，如2003年后的Landsat-7影像，對(duì)這些區(qū)域的研究，統(tǒng)一稱為數(shù)據(jù)“缺失”區(qū)，由于灰度值的自相關(guān)性[13-14]，用灰度值推演出的石漠化指數(shù)也應(yīng)存在自相關(guān)性，因而對(duì)“缺失”區(qū)數(shù)據(jù)的估計(jì)采用地統(tǒng)計(jì)中的克里格法預(yù)測.

從圖1中可以發(fā)現(xiàn)恩施市的行政邊界極不規(guī)則，東西、南北、西北-東南、西南-東北方向的間距在60～100 km之間，石漠化指數(shù)的采樣間距除去“缺失”區(qū)外為30 m，可認(rèn)為是規(guī)則的采樣點(diǎn)，可將采樣間距或倍數(shù)設(shè)為滯后距離[28]，因此滯后距離設(shè)為3 km，滯后組數(shù)為12，變程為R=36 km，說明像素之間的距離小于該值時(shí)，它們之間空間相關(guān)，反之則認(rèn)為無關(guān)，即相互獨(dú)立.

由于原始數(shù)據(jù)經(jīng)對(duì)數(shù)變換后服從正態(tài)分布，分別采用普通克里格法(Ordinary Kriging)、簡單克立里格法(Simple Kriging)對(duì)“缺失”區(qū)進(jìn)行數(shù)據(jù)估算，變程均為36 km，滯后距離3 km，滯后組數(shù)12，數(shù)據(jù)按隨機(jī)原則從總體中抽取的訓(xùn)練數(shù)據(jù)，得半變異函數(shù)模型為指數(shù)模型.為了對(duì)其擬合精度檢驗(yàn)證，利用5%的沒有參與建模的數(shù)據(jù)，按照式(3)、(4)計(jì)算得圖4、圖5，可以發(fā)現(xiàn)無論決定系數(shù)(r)、平均絕對(duì)誤差(MAE)、均方根誤差(RMSE)均以普通克里格法最優(yōu)，塊金值C0=0.024 1，拱高C=0.008 1，變程R=36 km，半變異指數(shù)函數(shù)模型，如圖6所示，由此可得石漠化指數(shù)的Kriging插值，如圖7所示.

圖4 普通克里格法預(yù)測石漠化指數(shù)精度圖 圖5 簡單克里格法預(yù)測石漠化指數(shù)精度圖

圖6 石漠化指數(shù)半變異函數(shù)圖 圖7 石漠化指數(shù)的Kriging插值圖

圖8 石漠化等級(jí)分布圖Fig.8 Distribution map of karst rocky desertification(K.R.D)

由上面的分析和圖6可以發(fā)現(xiàn)基臺(tái)值C0+C=0.032 2，C0/(C0+C)=74.844 7%，表示隨機(jī)因素所引起的空間異質(zhì)性占總變異的比例，該數(shù)值大，說明石漠化指數(shù)的空間相關(guān)性弱，反之則相關(guān)性強(qiáng)，由結(jié)構(gòu)性引起的空間變異度較大.經(jīng)研究表明該數(shù)值小于25%，空間相關(guān)性強(qiáng)，25%～75%空間相關(guān)性中等，大于75%中間相關(guān)性弱[13]，顯然74.844 7%小于75%，屬于中等空間相關(guān).結(jié)構(gòu)性因素巖性、氣候和隨機(jī)性因素土地耕作、林木采伐等對(duì)石漠化指數(shù)變異的影響差異不大.

3.3 石漠化等級(jí)的空間分布

圖7為石漠化指數(shù)分布圖，依據(jù)式(2)將其轉(zhuǎn)變基巖裸露率分布圖；坡度的獲取從SRTM1數(shù)據(jù)中提取，然后據(jù)石漠化等級(jí)劃分標(biāo)準(zhǔn)(表1)，可得恩施市2015年的石漠化等級(jí)分布(圖8)，統(tǒng)計(jì)得輕度石漠化面積850.531 5 km2,中度石漠化面積27.267 3 km2，重度石漠化面積0.223 2 km2，石漠化的總面積為878.022 0 km2.

2015-10-22的Landsat-7遙感影像由于傳感器自身的原因存在條帶缺失，如果僅對(duì)條帶缺失區(qū)域的數(shù)據(jù)采用克里格插值法進(jìn)行完善，不考慮云、影陰的影響，可分為：情形①，不對(duì)影陰區(qū)進(jìn)行克里格插值，用原影像的光譜信息分析則輕度石漠化為108.895 5 km2，中度石漠化為0.082 8 km2，而處理后的結(jié)果為輕度石漠化為130.309 2 km2，中度為6.585 3 km2；情形②，對(duì)于云區(qū)，不進(jìn)行克里格處理，結(jié)果為輕度石漠化為4.471 2 km2,中度石漠化為0.234 0 km2，處理后輕度石漠化為2.686 5 km2，中度石漠化為0.042 3 km2,無重度,可以發(fā)現(xiàn)對(duì)于陰影區(qū)、云區(qū)二者結(jié)果差異較大.

4 結(jié)語

該研究選取了傳感器存在缺陷的Landsat-7數(shù)據(jù)源，條帶缺失明顯，存在大面積的數(shù)據(jù)缺失，云層受大氣狀況的影響，所選影像全景云量很低近似為零；陰影由于地形起伏、太陽高度角、成像時(shí)間等的影響，作者在研究時(shí)對(duì)云、陰影的位置和范圍進(jìn)行了確定，作為數(shù)據(jù)“缺失”區(qū)處理.

地統(tǒng)計(jì)學(xué)分析時(shí)，得半變異函數(shù)為指數(shù)模型，變程為36 km，說明在36 km的范圍內(nèi)石漠化指數(shù)是存在空間相關(guān)性的，可以認(rèn)為所有“缺失”區(qū)的圖斑均在該范圍內(nèi)，應(yīng)用克里格空間插值是存在意義的；對(duì)其預(yù)測精度的評(píng)估，現(xiàn)多采用在建模時(shí)預(yù)留一部分檢證數(shù)據(jù)和參與建模的數(shù)據(jù)與預(yù)測值的誤差大小來評(píng)價(jià)，而實(shí)際調(diào)查只能得到石漠化的狀況和面積大小，而無法獲得石漠化指數(shù)值.

對(duì)于條帶缺失區(qū)的數(shù)據(jù)修復(fù)是必做的，云區(qū)、陰影區(qū)作者做了一個(gè)對(duì)比試驗(yàn)，即使用原有的光譜輻射值和克里格空間插值兩種方法對(duì)比，結(jié)果波動(dòng)性較大.對(duì)石漠化指數(shù)的空間結(jié)構(gòu)分析表明，由于受氣候、巖性、坡度和人們經(jīng)營措施等的影響，研究區(qū)域內(nèi)以石漠化指數(shù)為指標(biāo)的空間結(jié)構(gòu)呈中等空間自相關(guān)，恩施市的石漠化以輕度石漠化為主，面積850.531 5 km2，占國土總面積的21.413 2%，主要分布在研究區(qū)域的西北和東南地區(qū)，若采取有效的措施石漠化將得到有效的控制，轉(zhuǎn)為無石漠化地區(qū)，反之則相反.