青海大場金礦礦集區(qū)1988-2016年多時相植被覆蓋率估算及變化檢測研究

張 肅

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

1 土地植被覆蓋率計算的基本原理

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

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

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

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

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

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

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

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

1.1 分類前比較法

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

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

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

1.2 分類后比較法

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

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

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

2 研究區(qū)域概況

研究區(qū)域位于青海省大場金礦集礦區(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ù)

本次實驗主要是對植被覆蓋率估算，因此在實驗區(qū)域影像盡量選擇植被相對茂盛的季節(jié)。在數(shù)據(jù)源選取上，選擇易獲取的免費Landsat系列數(shù)據(jù)。下載的遙感影像是青海省大場金礦的四個時間點的數(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個時相的遙感影像略圖如圖2所示：

3.2 分類系統(tǒng)

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

4 研究方法

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

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

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

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

圖3 植被覆蓋率計算流程圖

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

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

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

4.2 NDVI計算

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

表2 Landsat8影像NDVI計算值

4.3 土壤線方程

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

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

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

4.4 土地植被覆蓋分類

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

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

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

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

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

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

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

4.5 植被的亞像元模型

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

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

表4 林地的NDVI植被覆蓋率

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

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

具體步驟如下：

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

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

通過以上公式，

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

當i=1時，DNs1是將r∞1、c1、DN∞1和DNi代入公式中，DNi是TM band3的DN值。

當i=2時，DNs2是將r∞2、c2、DN∞2和DNi代入公式中，DNi是TM band4的DN值。

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

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

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

5 研究結(jié)果

5.1 植被覆蓋率成圖

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

表7 植被覆蓋率分級

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

圖4 4個時相遙感影像的植被覆蓋率圖

圖5 4幅時相的遙感植被覆蓋率變化檢測圖

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

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

接上表

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

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

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

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

通過以上標準，便可以得到3幅變化檢測影像，分別為1996年～1988年，2004年～1996年，2016年～2004年。

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

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

6 結(jié)論

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

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