張 肅

（四川省地質(zhì)調(diào)查院，四川 成都 610081）

1 土地植被覆蓋率計(jì)算的基本原理

土地及植被覆蓋是在全球環(huán)境變化，全球變暖研究中占據(jù)重要的地位。在分析人類活動(dòng)對(duì)于自然環(huán)境造成的影響以及自然環(huán)境的反映中，也具有重要意義。土地植被覆蓋分類為土地植被覆蓋變換監(jiān)測(cè)，氣候檢測(cè)，土壤理化性質(zhì)等提供強(qiáng)有力的數(shù)據(jù)支撐。

遙感能夠進(jìn)行空間連續(xù)覆蓋和長(zhǎng)時(shí)間持續(xù)觀測(cè)，因而成為土地植被覆蓋數(shù)據(jù)獲取的重要手段。中等空間分辨率的Landsat衛(wèi)星數(shù)據(jù)具有顯著表達(dá)地表覆蓋特征的能力，在土地植被利用分類與檢測(cè)領(lǐng)域發(fā)揮著重要作用[1]。在利用Landsat衛(wèi)星數(shù)據(jù)進(jìn)行土地植被制圖和變化監(jiān)測(cè)的研究中，Gong等人對(duì)8900多景TM和ETM+影像進(jìn)行分類，得到了全球土地覆蓋數(shù)據(jù)集[2]。王雪梅等人使用ETM+影像和副主數(shù)據(jù)對(duì)新疆于田綠洲進(jìn)行了土地植被覆蓋分類研究[3]。這些研究標(biāo)明，Landsat衛(wèi)星數(shù)據(jù)已經(jīng)成為土地植被資源調(diào)查的主要數(shù)據(jù)源之一，適合區(qū)域和全球等不同尺度的研究和應(yīng)用。

青海大場(chǎng)金礦是亞洲最大的金礦之一，對(duì)青海大場(chǎng)金礦地區(qū)的土地植被覆蓋研究具有重要意義，本研究時(shí)間跨度長(zhǎng)達(dá)30年，計(jì)劃選取4期不同時(shí)間的影像，影像間時(shí)間跨度約為30年，這樣的影像選擇可以全面的反映30年來研究區(qū)域內(nèi)土地植被覆蓋的變化規(guī)律。對(duì)于2013年之后的數(shù)據(jù)選擇，本研究項(xiàng)目選取Landsat8數(shù)據(jù)，對(duì)于2013年之前的數(shù)據(jù)選擇，本研究項(xiàng)目選取Landsat TM數(shù)據(jù)[4]。

目前，利用Landsat8衛(wèi)星影像進(jìn)行土地植被覆蓋分類的研究處于快速發(fā)展之中。Landsat8陸地成像儀（OLI）多光譜影像相對(duì)于TM/ETM+增加了新的特性。在土地植被覆蓋分類過程之中，衛(wèi)星影像預(yù)處理和分類方法會(huì)影響植被分類結(jié)果的經(jīng)度，因此對(duì)這些處理方法進(jìn)行研究和探討具有重要意義。本文從衛(wèi)星數(shù)據(jù)特性出發(fā)，結(jié)合處理中影響分類精度的因素，對(duì)青海大場(chǎng)金礦地區(qū)進(jìn)行土地植被覆蓋率分類方法研究，定量評(píng)價(jià)處理方法的實(shí)用性。

利用遙感資料估算植被覆蓋率的方法大致可歸納為兩種[5]：

（1）經(jīng)驗(yàn)?zāi)Ｐ头?，通過建立實(shí)測(cè)植被覆蓋率數(shù)據(jù)與植被指數(shù)的經(jīng)驗(yàn)?zāi)Ｐ蛠砬笕〈竺娣e植被覆蓋率[6]。該方法需要到研究區(qū)域內(nèi)進(jìn)行野外調(diào)查，對(duì)研究區(qū)域內(nèi)進(jìn)行分散采樣，根據(jù)實(shí)地采集的植被覆蓋率數(shù)據(jù)，并結(jié)合遙感衛(wèi)星影像計(jì)算出的植被指數(shù)，建立統(tǒng)計(jì)的經(jīng)驗(yàn)?zāi)Ｐ?，然后將得出的?jīng)驗(yàn)?zāi)Ｐ蛻?yīng)用于整個(gè)研究區(qū)域，從而得到整個(gè)研究區(qū)域的土地植被覆蓋率[7]。該方法需要消耗大量的時(shí)間和資源，并且僅能對(duì)當(dāng)前時(shí)段的植被覆蓋率進(jìn)行野外調(diào)查，無法獲取以前的植被覆蓋率信息。故該方法多用于單時(shí)相的土地植被覆蓋率調(diào)查研究。

（2）植被指數(shù)轉(zhuǎn)換法：通過對(duì)各像元中植被類型及分布特征的分析，建立植被指數(shù)與植被覆蓋率的轉(zhuǎn)換關(guān)系來直接估算植被覆蓋率[8]。這種方法的優(yōu)點(diǎn)是可以根據(jù)遙感影像直接進(jìn)行土地植被覆蓋率的計(jì)算分析，而不需要進(jìn)行野外調(diào)查。因此可以對(duì)多時(shí)相的遙感影像進(jìn)行分析。本次研究采取的就是植被指數(shù)轉(zhuǎn)換法[9]。

植被覆蓋率計(jì)算得到結(jié)果后，就可以進(jìn)行變化檢測(cè)。多時(shí)相數(shù)據(jù)的變化檢測(cè)一般可以分為兩大類別：

1.1 分類前比較法

分類前比較法是先對(duì)不同時(shí)相的圖像進(jìn)行比較，比較的內(nèi)容可以是直接的光譜值，也可以是從光譜信息中間接提取的特征值，通過比較增強(qiáng)了變化信息和未變化信息的可分性，最后對(duì)比較的結(jié)果進(jìn)行分類得到變化檢測(cè)圖像。分類前比較方法又可以分為直接比較法和間接比較法。

直接比較法直接對(duì)光譜值進(jìn)行比較，簡(jiǎn)單直觀，但是由于多時(shí)相遙感數(shù)據(jù)往往存在同譜異物和同物異譜現(xiàn)象，檢測(cè)結(jié)果容易產(chǎn)生虛假變化。其次，直接比較法是基于像素的直接比較，對(duì)圖像預(yù)處理（如圖像配準(zhǔn)等）的精度要求很高。因此，直接比較法對(duì)圖像的適應(yīng)性較差[10]。

間接比較法是相對(duì)于直接比較法而言，它不是直接對(duì)光譜特征進(jìn)行比較，而是首先對(duì)光譜特征進(jìn)行必要的處理，提取有用信息，也是為了增加變化信息和未變化信息的可分性[11]。

1.2 分類后比較法

分類后比較法也是目前常用的變化檢測(cè)方法，它首先對(duì)配準(zhǔn)后的不同時(shí)相的圖像分別用相同的分類算法進(jìn)行分類，然后將分類結(jié)果逐像元比較以確定變化信息。其中對(duì)圖像中地物的準(zhǔn)確分類是這種方法的難點(diǎn)。因?yàn)榉诸惡蟊容^法中圖像分類是分別進(jìn)行的，所以它回避了分類前比較法所要求的不同時(shí)相的影像成像條件基本一致的條件及影像間輻射校正、匹配等問題，并且它可以同時(shí)獲得變化前后的類別信息和變化信息。但是它對(duì)分類方法的選擇以及分類器的精度要求較高。

本次研究采用的是分類后比較法，因?yàn)槊嫦驅(qū)ο蠓诸惙椒梢员ＷC土地植被分類的精度，產(chǎn)生的分類結(jié)果較為準(zhǔn)確可信，所以分類后比較法用于此研究較為合適。

本次研究使用植被指數(shù)轉(zhuǎn)換法計(jì)算植被覆蓋率，數(shù)據(jù)預(yù)處理和NDVI計(jì)算使用ENVI及ERDAS軟件實(shí)現(xiàn)，土地植被分類采用ecognition軟件實(shí)現(xiàn)，植被覆蓋率計(jì)算采用MATLAB軟件編程實(shí)現(xiàn)。

2 研究區(qū)域概況

研究區(qū)域位于青海省大場(chǎng)金礦集礦區(qū)，位于96°00'E～96°25'E,35°10'N～35°30'之間的區(qū)域，平均海拔約3600米左右。主要的植被類型為草地、裸地、雪地、灌叢等四種?？偯娣e約為1034平方公里。研究區(qū)域范圍如圖1所示：

圖1 研究區(qū)域范圍

3 數(shù)據(jù)源及植被分類系統(tǒng)

3.1 遙感數(shù)據(jù)

本次實(shí)驗(yàn)主要是對(duì)植被覆蓋率估算，因此在實(shí)驗(yàn)區(qū)域影像盡量選擇植被相對(duì)茂盛的季節(jié)。在數(shù)據(jù)源選取上，選擇易獲取的免費(fèi)Landsat系列數(shù)據(jù)。下載的遙感影像是青海省大場(chǎng)金礦的四個(gè)時(shí)間點(diǎn)的數(shù)據(jù)，分別來自于Landsat TM數(shù)據(jù)以及Landsat 8影像數(shù)據(jù)。數(shù)據(jù)下載來自于地理空間數(shù)據(jù)云網(wǎng)站(http://www.gscloud.cn/)。

表1 遙感影像數(shù)據(jù)信息

研究區(qū)域的4個(gè)時(shí)相的遙感影像略圖如圖2所示：

3.2 分類系統(tǒng)

參考全國(guó)遙感監(jiān)測(cè)土地利用/土地覆蓋分類體系及遙感數(shù)據(jù)土地利用/土地覆蓋分類體系，該體系對(duì)于不同空間尺度和對(duì)應(yīng)遙感數(shù)據(jù)源，都具有相對(duì)應(yīng)的分類，且類型逐漸細(xì)化。根據(jù)該分類體系，結(jié)合Landsat8數(shù)據(jù)特性及研究對(duì)象地物分布的特點(diǎn)，定義五大類土地利用覆蓋類型，各類型如下：①水域；②灌從；③草地；④冰雪；⑤林地；⑥裸地。

4 研究方法

研究方法如圖3所示，首先對(duì)覆蓋區(qū)域內(nèi)的Landsat8影像進(jìn)行預(yù)處理，包括輻射校正,大氣校正，正射校正及幾何校正。并按照研究區(qū)域的范圍進(jìn)行影像裁剪。其次，將預(yù)處理之后的影像進(jìn)行分類處理，并評(píng)價(jià)分類精度，得到土地植被覆蓋分類圖。同時(shí)，計(jì)算覆蓋區(qū)域的NDVI值，根據(jù)分類情況和NDVI值，通過植被指數(shù)轉(zhuǎn)換法得出研究區(qū)域內(nèi)的植被覆蓋率。

圖2 研究區(qū)域4個(gè)時(shí)相遙感影像

4.1 數(shù)據(jù)預(yù)處理

（1）輻射校正。下載的Landsat8/TM L1T級(jí)數(shù)據(jù)，元數(shù)據(jù)信文件（124028_MTL）有詳細(xì)影像參數(shù)，由于研究使用數(shù)據(jù)為L(zhǎng)1T級(jí)數(shù)據(jù)，經(jīng)過系統(tǒng)輻射校正的數(shù)據(jù)。由于1級(jí)產(chǎn)品的DN值是由輻射亮度線性變換得到的，因此從1級(jí)產(chǎn)品計(jì)算輻射亮度只需利用相關(guān)參數(shù)（Gain和Bias）進(jìn)行線性反變換即可，計(jì)算過程比較簡(jiǎn)單。采用ENVI4.7軟件進(jìn)行輻射校正操作。

圖3 植被覆蓋率計(jì)算流程圖

（2）大氣校正。6S模型是常用的大氣校正方法，利用6S模型將原始多光譜影像大氣表觀反射率（TOA）校正為地表反射率，6S校正使用Vermote等人編寫的6S軟件完成。

（3）正射校正。使用基于DEM的正射模型對(duì)6S校正后的影像進(jìn)行地形校正，去除地形效應(yīng)對(duì)地物反射率的影響。DEM的數(shù)據(jù)為30m分辨率的ASTER全球DEM，可以在于地理空間數(shù)據(jù)云網(wǎng)站(http://www.gscloud.cn/)下載。正射校正較好消除了6S校正后影像地形的3D效果，陰坡和陽坡的亮度差異減少，地形效應(yīng)減弱。

（4）幾何校正。本次研究使用經(jīng)過正射校正的影像以及地形圖進(jìn)行幾何校正，地形圖為1:50000比例尺，采用ENVI軟件進(jìn)行幾何校正，校正結(jié)果的RMSE為0.5個(gè)像素，對(duì)于該研究的尺度，該精度已經(jīng)可以滿足要求。

4.2 NDVI計(jì)算

使用大氣校正后得到的反射率數(shù)據(jù)，直接用ENVI軟件計(jì)算NDVI，具體計(jì)算公式為NDVI=(NIR-R)/(NIR+R)。NDVI取值范圍為[-1,1]，負(fù)值表示地面覆蓋為云、水、雪等，對(duì)可見光高反射；0表示有巖石或裸土等，NIR（近紅外波段）和R（紅色波段）近似相等；正值，表示有植被覆蓋，且隨覆蓋度大而增大。

表2 Landsat8影像NDVI計(jì)算值

4.3 土壤線方程

土壤線指土壤的光譜值在近紅外波段和紅波段的反射率或亮度值所構(gòu)成的二維平面上的線性關(guān)系，是對(duì)大量土壤光譜信息的綜合描述。

土壤線方程為：NIR=a R+b其中：

NIR—近紅外波段的反射率或DN值；R—紅色波段的反射率或DN值；a和b—分別為土壤線的斜率和截距。TM影像的band 3是紅色波段，band 4是近紅外波段。用R軟件做土壤線方程的線性回歸計(jì)算。上表通過回歸求得a=0.666，b=-15.871

4.4 土地植被覆蓋分類

對(duì)于研究區(qū)域內(nèi)土地植被覆蓋分類，基于像素的方法很難精確地提取各類地物的邊界，并且容易出現(xiàn)錯(cuò)分和漏分的情況。近年來，面向?qū)ο蟮倪b感圖像分類方法逐漸開始發(fā)展起來，其中以econgnition軟件最為成熟。面向?qū)ο蟮姆诸惙椒ㄏ鄬?duì)于傳統(tǒng)的基于像素的分類方法的優(yōu)點(diǎn)是：模擬人腦解譯方式，充分利用高分辨率遙感影像豐富的紋理和結(jié)構(gòu)多種特征，提高信息分類精度；更適合處理空間尺度、空間分析等問題；面向?qū)ο蠓椒梢栽诓煌叨瓤臻g提取特定主題的信息，提高分類精度。

在本研究中，面向?qū)ο蟮耐恋刂脖桓采w分類提取的基本步驟為。

（1）根據(jù)圖像的光譜特征和幾何特征，將整幅影像進(jìn)行分割，將圖像整體分割為大小不一的小塊，這些小塊具有相似的光譜特征和幾何特征。

（2）分割出小的分塊之后，根據(jù)一系列的分類規(guī)則對(duì)這些小的分塊進(jìn)行分類。分類的參數(shù)選擇上，可以選取每一個(gè)波段的均值和方差，以及形態(tài)因子作為分類的依據(jù)。在設(shè)定好總的地物類別之后（包括裸地，灌從，河流，冰雪等類別），對(duì)每一個(gè)小分塊賦予一個(gè)類別。

（3）最后，對(duì)于剩下的錯(cuò)誤類別，可以采用手動(dòng)分類進(jìn)行修正。給每一個(gè)分錯(cuò)的小塊直接指定正確的類別；這對(duì)于最后修正某些類別的錯(cuò)誤邊界有很好的效果。至此，Landsat8影像分類得到了完整的結(jié)果，包括灌從，草地、林地、裸地、冰雪、河流等地物類別。

（4）分類結(jié)果評(píng)價(jià)，以2016年的Landsat8影像為例，通過人工判讀方式在研究區(qū)域隨機(jī)選擇均勻分布的樣本點(diǎn)3456個(gè)，樣本點(diǎn)的選取參考谷歌地球的worldview2影像。得到的精度統(tǒng)計(jì)表格如下，可以看出，采用面向?qū)ο蠓诸惸艿玫捷^好的分類結(jié)果。

表3 分類精度比較結(jié)果

4.5 植被的亞像元模型

（1）等密度模型(林地)。假設(shè)像元中植被類型教委單一并且植被垂直密度足夠高，即LAI→∞，對(duì)應(yīng)的，因此，在等密度模型假設(shè)之下，植被覆蓋率：

根據(jù)公式（1），首先求得NDVI的最大值和最小值，在ENVI中，利用分類后的影像和NDVI影像做波段運(yùn)算得到林地類別的NDVI值。其次，根據(jù)公式（1）求解林地植被覆蓋率。在ENVI中，利用分類后影像和NDVI影像做波段運(yùn)算，得到林地的NDVI值如表4所示：

表4 林地的NDVI植被覆蓋率

（2）非等密度模型（草地、灌從、裸地）。非密度模型與等密度模型類似，該模型也假設(shè)像元中植被類型較為單一，但此時(shí)植被垂直密度較小，即LAI＜＜∞。因此植被覆蓋部分的NDVIg值，即，需要由Bear'定律確定，k=1,如下面公式所示：

其中LAI可以通過兩向近似法進(jìn)行的低垂直密度植被的LAI估算。

具體步驟如下：

①求草地、灌從的NDVI最大值和最小值，方法與林地中第一步相同。

②求解LAI（葉面積指數(shù)）。使用如下公式求解LAI。

通過以上公式，

參數(shù)取值如表5所示：

當(dāng)i=1時(shí)，DNs1是將r∞1、c1、DN∞1和DNi代入公式中，DNi是TM band3的DN值。

當(dāng)i=2時(shí)，DNs2是將r∞2、c2、DN∞2和DNi代入公式中，DNi是TM band4的DN值。

將上述的DNs1和DNs2代入土壤線方程，DNs2=a＊DNs1+b。這樣就得到一個(gè)只有LAI一個(gè)未知數(shù)的方程。但是這個(gè)方程中的LAI不能直接求解，因?yàn)閕=1時(shí)，e2c1LAI=e1.2LAI；i=2時(shí)，e2c2LAI=e0.42LAI。

采用牛頓迭代法求解LAI：使用Matlab編程用逼近的方式求解LAI。

求得LAI之后，可以計(jì)算出植被覆蓋率fg，使用MATLAB編程實(shí)現(xiàn)。具體數(shù)據(jù)如表6所示：

5 研究結(jié)果

5.1 植被覆蓋率成圖

將水體、裸地以及云的植被覆蓋率賦值為0，即可得到植被覆蓋率圖像。為了方便解譯，將植被覆蓋率分為5個(gè)等級(jí)。

表7 植被覆蓋率分級(jí)

根據(jù)不同的級(jí)別，對(duì)不同的覆蓋率級(jí)別不同的顏色，4個(gè)時(shí)相的植被覆蓋率分級(jí)圖如圖4所示。

圖4 4個(gè)時(shí)相遙感影像的植被覆蓋率圖

圖5 4幅時(shí)相的遙感植被覆蓋率變化檢測(cè)圖

4個(gè)時(shí)相的遙感影像的植被覆蓋率統(tǒng)計(jì)表如表8所示：

表8 1988-2016年植被覆蓋率統(tǒng)計(jì)表1988年

接上表

從表7可以看出，類別4和類別5百分比之和，近30年來，開始呈現(xiàn)減小趨勢(shì)，而到了2016年，又轉(zhuǎn)而增加。說明了該研究區(qū)域的植被顯示出“減少-增加”的反轉(zhuǎn)趨勢(shì)，這與30年來的氣候變化和人類工程活動(dòng)，環(huán)境治理等多方面因素有關(guān)。而研究區(qū)域內(nèi)的絕大部分土地植被覆蓋率小于40%，表明區(qū)域內(nèi)植被覆蓋不理想。

5.2 植被覆蓋率分類后變化檢測(cè)結(jié)果

植被覆蓋率被分為5類，根據(jù)4幅不同時(shí)相的植被覆蓋率分類圖，可以對(duì)這4幅圖進(jìn)行變化檢測(cè)。本研究采用的方法是分類后比較法，將植被覆蓋率分類結(jié)果逐像元比較以確定變化信息，該程序在MATLAB中實(shí)現(xiàn)。變化檢測(cè)的定量分析標(biāo)準(zhǔn)如下：

①前后兩幅影像同一像元級(jí)別從低變化至高，將其劃分為增加。②前后兩幅影像同一像元級(jí)別從沒有變化，將其劃分為不變。③前后兩幅影像同一像元級(jí)別從高變化至低，將其劃分為減少。

通過以上標(biāo)準(zhǔn)，便可以得到3幅變化檢測(cè)影像，分別為1996年～1988年，2004年～1996年，2016年～2004年。

從圖5可以看出，1996你那～1988年時(shí)間段，中部地區(qū)植被覆蓋率較少，而下方部分地區(qū)植被覆蓋率增加。2004-1996年，左邊部分地區(qū)植被覆蓋率減小，而中上部部分地區(qū)植被覆蓋率增加。而2016年～2004年之間，周圍大部分植被覆蓋率增加。

由于變換檢測(cè)方法采用的分類之后比較的方法，由于面向?qū)ο蟮姆椒ūＷC了分類的精度，植被覆蓋率的計(jì)算也相應(yīng)的保證了精度。所以變化檢測(cè)的結(jié)果是真實(shí)可信的。同時(shí)，將分析結(jié)果與4個(gè)時(shí)間點(diǎn)的遙感影像相對(duì)比，可以發(fā)現(xiàn)，植被覆蓋率的變化趨勢(shì)和遙感影像的變化趨勢(shì)是一致的。

6 結(jié)論

本研究通過對(duì)青海大場(chǎng)金礦礦集區(qū)4個(gè)時(shí)相的遙感影像數(shù)據(jù)進(jìn)行分析處理，得到了4個(gè)時(shí)相的土地植被覆蓋率分類影像，其分類的準(zhǔn)確性達(dá)到了較好的水平。

同時(shí)，根據(jù)4個(gè)時(shí)相的植被覆蓋率的統(tǒng)計(jì)結(jié)果，對(duì)4個(gè)時(shí)相的植被覆蓋率進(jìn)行了變化檢測(cè)，分析了近30年來植被覆蓋率的變化趨勢(shì)，并與遙感影像進(jìn)行對(duì)比，證明了植被覆蓋率變化趨勢(shì)和遙感影像的變化趨勢(shì)相同，從而對(duì)結(jié)果的可靠性進(jìn)行了驗(yàn)證。