基于Landsat-7 ETM SLC-OFF影像的巖溶區(qū)石漠化研究

馬友平

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

石漠化是一種石質(zhì)荒漠化，有地球癌癥之稱[1]，是土壤受到侵蝕后的直接后果與巖溶地區(qū)生態(tài)系統(tǒng)逆向演替的頂級(jí)階段，一旦發(fā)生，恢復(fù)治理相當(dāng)困難，導(dǎo)致生態(tài)環(huán)境的崩潰，且恢復(fù)速率極慢，已成為我國(guó)三大生態(tài)自然災(zāi)害之一[2-4].在我國(guó)主要分布在以云貴高原為中心的西南八省、市、自治區(qū)的451個(gè)縣(市)中，這樣大范圍內(nèi)有關(guān)石漠化基礎(chǔ)數(shù)據(jù)的獲取一直是石漠化研究的核心問(wèn)題之一，其中遙感影像已成為石漠化研究的主要數(shù)據(jù)源，遙感影像主要依據(jù)的是地物的波譜特性，易受成像時(shí)間、地點(diǎn)、地物類別等的影響較大，現(xiàn)還沒(méi)有一個(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ù)庫(kù)，分析影響因子，認(rèn)為石漠化產(chǎn)生的主導(dǎo)原因是人類不合理的經(jīng)營(yí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ì)貴州省的石漠化演變過(guò)程進(jìn)行了分析和評(píng)價(jià)；陳起偉等[10]以TM、ASTER和ALOS 三期衛(wèi)星影像為遙感動(dòng)態(tài)監(jiān)測(cè)數(shù)據(jù)源對(duì)貴州省的石漠化進(jìn)行了遙感監(jiān)測(cè)；靖娟利等[11]利用MODIS遙感影像對(duì)廣西自治區(qū)的喀斯特石漠化演變特征進(jìn)行了研究；童立強(qiáng)等[12]應(yīng)用ETM、TM、CBERS-02B、ALOS和ASTER衛(wèi)星遙感數(shù)據(jù)對(duì)中國(guó)南方地區(qū)石漠化進(jìn)行了遙感調(diào)查與演變分析.上述研究其基礎(chǔ)數(shù)據(jù)的獲取主要來(lái)自遙感影像，以遙感影像的光譜特性為基礎(chǔ)，以不同波段之間的線性與非線性組合來(lái)提取能反映石漠化狀況的指標(biāo).為此本研究擬在借鑒上述成果時(shí)以湖北省恩施市為例，討論石漠化范圍的確定和等級(jí)的劃分標(biāo)準(zhǔn)；對(duì)于指標(biāo)提取時(shí)影像的“缺失”區(qū)通過(guò)插值法獲得，有研究證明遙感影像各點(diǎn)的灰度值(DN)，并不相互獨(dú)立，它們之間存在一定的相關(guān)性[13-14]，為“缺失”區(qū)石漠化信息的地統(tǒng)計(jì)學(xué)估計(jì)提供了可能.

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

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)過(guò)石漠化的專項(xiàng)治理，石漠化現(xiàn)狀有所改善，但問(wèn)題依然存在，治理任務(wù)依然十分艱巨.

1.2 數(shù)據(jù)來(lái)源

美國(guó)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ù)稱之為L(zhǎng)andsat-7 ETM SLC-OFF數(shù)據(jù)，該研究主要采用有條帶缺失的Landsat-7 ETM SLC-OFF影像，時(shí)間2015-10-22，軌道號(hào)/行號(hào)為126/39，空間分辨率30 m，云量0，而通過(guò)實(shí)際的圖片觀察恩施市上空有少量的云層；全市的DEM數(shù)據(jù)采用SRTM1(Shuttle Radar Topography Mission)，對(duì)應(yīng)的分辨率精度為30 m，來(lái)源地理空間數(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ù)問(wèn)題，以彌補(bǔ)在應(yīng)用過(guò)程中因云量、時(shí)段等導(dǎo)致的數(shù)據(jù)源缺乏的問(wèn)題，下面同時(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è)國(guó)土面積的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)過(guò)查閱文獻(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í)的劃分

石漠化主要分布在我國(guó)廣大的南方地區(qū)，面積大，地域差異明顯，由于研究方法和數(shù)據(jù)的差異，現(xiàn)在還沒(méi)有一個(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í)，分別為無(wú)石漠化(I)、輕度石漠化(II)、中度石漠化(III)和重度石漠化(IV)，詳見(jiàn)表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)行線性無(wú)偏最優(yōu)估計(jì)的一種方法，是地統(tǒng)計(jì)學(xué)的重要研究?jī)?nèi)容之一，以變異函數(shù)來(lái)研究空間的結(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)還沒(méi)有統(tǒng)一的評(píng)定方法和嚴(yán)格的定義，在大尺度研究時(shí)，多以遙感影像為數(shù)據(jù)源，通過(guò)光譜特征來(lái)計(jì)算，對(duì)于“缺失”區(qū)預(yù)測(cè)數(shù)據(jù)，無(wú)法同地形因子如高程等，與地面實(shí)測(cè)值相較，評(píng)價(jià)擬合結(jié)果，僅能將參與擬合的數(shù)據(jù)按隨機(jī)的方法抽選一部分作為檢驗(yàn)點(diǎn)，該研究中抽取5%的數(shù)據(jù)，而剩余95%的數(shù)據(jù)參與克里格法預(yù)測(cè)建模.建模精度擬采用三個(gè)指標(biāo)評(píng)價(jià)：①?zèng)Q定系數(shù)(r)，用以反映觀測(cè)值與估計(jì)值的關(guān)系緊密程度，值越大估測(cè)精度越高；②平均絕對(duì)誤差(MAE)，反映估測(cè)值與檢驗(yàn)值的偏離大小，其值越接近0效果越好；③均方根誤差(RMSE)，介于0與1之間，取值越接近0，估測(cè)越精確.MAE和RMSE的計(jì)算公式如下:

(3)

(4)

式中：O實(shí)測(cè)值，P預(yù)測(cè)值，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)分布，滿足克里格分析中簡(jiǎn)單克里格、普通克里格等要求數(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ù)測(cè).

從圖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，說(shuō)明像素之間的距離小于該值時(shí)，它們之間空間相關(guān)，反之則認(rèn)為無(wú)關(guān)，即相互獨(dú)立.

由于原始數(shù)據(jù)經(jīng)對(duì)數(shù)變換后服從正態(tài)分布，分別采用普通克里格法(Ordinary Kriging)、簡(jiǎn)單克立里格法(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%的沒(méi)有參與建模的數(shù)據(jù)，按照式(3)、(4)計(jì)算得圖4、圖5，可以發(fā)現(xiàn)無(wú)論決定系數(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ù)測(cè)石漠化指數(shù)精度圖 圖5 簡(jiǎn)單克里格法預(yù)測(cè)石漠化指數(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ù)值大，說(shuō)明石漠化指數(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,無(wú)重度,可以發(fā)現(xiàn)對(duì)于陰影區(qū)、云區(qū)二者結(jié)果差異較大.

4 結(jié)語(yǔ)

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

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

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