基于像元二分模型的植被覆蓋度反演-以北京市為例

1、基于像元二分模型的植被覆蓋度反演-以北京市為例 王玲（西北大學 城市與環(huán)境學院，陜西 西安 710127）摘要：采用遙感技術監(jiān)測植被覆蓋度具有重要意義。本文以北京市為例，基于2013年的Landsat8 OLI影像，選取NDVI值為參數，采用像元二分模型對植被覆蓋度進行反演，最終反演的結果與實際情況符合，說明采用該方法反演植被覆蓋度可行。關鍵詞：植被覆蓋度、像元二分模型、NDVI、植被指數引言植被覆蓋度（Vegetation fractional cover，簡稱fc）是指植被（包括葉、莖、枝）在地面的垂直投影面積占統計區(qū)總面積的百分比，即植土比。通常林冠稱郁閉度，灌草等植被稱覆蓋度1。它是衡

2、量地表植被覆蓋的一個最重要的指標，被覆蓋度及其變化是區(qū)域生態(tài)系統環(huán)境變化的重要指示，對水文、生態(tài)、全球變化等都具有重要意義2。根據監(jiān)測手段, 測量植被覆蓋度的方法可分為地面測量和遙感測量兩大類，測量常用于田間尺度，遙感估算常用于區(qū)域尺度。目前已經發(fā)展了很多利用遙感測量植被覆蓋度的方法，較為實用的地面測量單的方法就是目估法，缺點主要是主觀性太強?？陀^的測量方法有樣點法、樣方法、樣帶法等，借助于采樣儀器的測量方法，空間定量計、移動光量計等。這些方法雖然提高了測量精度，但野外操作不便，并且成本較高, 難以在大范圍內快速提取植被覆蓋度。而采用遙感技術為監(jiān)測大面積區(qū)域的植被覆蓋度，甚至全球的植被覆蓋度監(jiān)

3、測提供了可能3。目前已有許多利用遙感技術測量植被覆蓋度的方法，其中應用最廣泛的方法是利用植被指數近似估算植被覆蓋度，常用的植被指數為NDVI。一、數據源 本文選取兩景覆蓋北京市的Landsat8 OLI影像、土地覆蓋類型圖以及北京行政邊界矢量數據為數據源。其中，土地覆蓋類型圖是作為掩膜文件使用，其目的是為了便于植被覆蓋度的估算；北京行政邊界矢量數據是為了將兩景鑲嵌好的影像數據進行裁剪使用的，其目的是裁剪出北京市行政區(qū)內的范圍。另外，Landsat8 OLI影像是從地理空間數據云網站上下載得到的，其成像時間為2013年10月份。與Landsat7的ETM+成像儀相比，OLI成像儀獲取的遙感圖像輻

4、射分辨率達到12比特，圖像的幾何精度和數據的信噪比也更高。OLI成像儀包括9個短波譜段（波段1波段9），幅寬185km，其中全色波段地面分辨率為15m，其他譜段地面分辨率為30m4。表1 Landsat8 OLI陸地成像儀波段參數 二、研究方法本文反演植被覆蓋度所采用的是像元二分模型方法，像元二分模型是一種簡單實用的遙感估算模型，它假設一個像元的地表由有植被覆蓋部分地表(SV)與無植被覆蓋部分地表(SS)組成，而遙感傳感器觀測到的光譜信息(S)也由這2個組分因子線性加權合成，各因子的權重是各自的面積在像元中所占的比率，如其中植被覆蓋度可以看作是植被的權重。因此，像元二分模型的原理如下：遙感傳感

5、器觀測到的光譜信息(S)由有植被覆蓋部分地表(SV)與無植被覆蓋部分地表(SS)組成，可得出： S = SV + SS ·········公式1假設一個像元中有植被覆蓋的面積比例為fc , 即該像元的植被覆蓋度, 則裸土覆蓋的面積比例為1 -fc ，如果全由植被所覆蓋的純像元所得的遙感信息為Sveg , 則混合像元的植被部分所貢獻的信息Sv可以表示為Sveg與fc的乘積：Sv =fc·Sveg ········

6、83;公式2那么， Ss =(1 -fc )·Ssoil ·········公式3將公式2與公式3代入到公式1中，可得到：S =fc ·Sveg +(1 -fc)S soil ·········公式4對公式4進行變換, 可得以下計算植被覆蓋度的公式：fc =(S -Ssoil) (Sveg -Ssoil ) ·······

7、83;·公式5其中Ssoil 為純土壤像元的信息, Sveg 為純植被像元的信息, 因而可以根據公式5利用遙感信息來估算植被覆蓋度。將歸一化植被指數(NDVI)代入公式5可以被近似為:fc =(NDVI - NDVIsoil ) (NDVIveg -NDVIsoil) ·········公式6其中, NDVIsoil 為裸土或無植被覆蓋區(qū)域的NDVI值, 即無植被像元的NDVI 值;而NDVIveg 則代表完全被植被所覆蓋的像元的NDVI 值, 即純植被像元的NDVI 值2。當區(qū)域內可以近似取V

8、FCmax=100%，VFCmin=0%，VFC = (NDVI -NDVImin)/ ( NDVImax -NDVImin)，NDVImax 和NDVImin分別為區(qū)域內最大和最小的NDVI值。由于不可避免存在噪聲，NDVImax 和NDVImin一般取一定置信度范圍內的最大值與最小值，置信度的取值主要根據圖像實際情況來定；當區(qū)域內不能近似取VFCmax=100%，VFCmin=0%，當有實測數據的情況下，取實測數據中的植被覆蓋度的最大值和最小值作為VFCmax和VFCmin，這兩個實測數據對應圖像的NDVI作為NDVImax 和NDVImin。當沒有實測數據的情況下，取一定置信度范圍內的N

9、DVImax 和NDVImin。VFCmax和VFCmin根據經驗估算。三、 數據處理1、數據預處理本文使用的Landsat8 OLI為L1T級別數據，不需做幾何校正處理。而北京市需要兩景Landsat OLI數據覆蓋，因此首先要進行圖像鑲嵌和裁剪，然后進行大氣校正等預處理過程。（1） 輻射定標輻射定標是將傳感器記錄的電壓或數字值轉換成絕對輻射亮度的過程。目的是消除傳感器本身所產生的誤差,由于傳感器在不斷的運行中光學器件性能逐漸退化,因此定標的系數也隨之不同,這些定標系數也在不斷的更改,在用戶獲得數據的時候,這些定標系數也在影像的頭文件中同時提供給用戶。遙感數據輻射定標就是將傳感器得到的灰度值

10、轉換為星上的輻射亮度值或星上反射率,即表觀輻射度或表觀反射率。輻射定標主要校正由傳感器的靈敏度帶來的輻射誤差5。其目的是為FLAASH大氣校正準備數據：定標符合單位要求的輻射量數據、轉換數據儲存順序等。該處理過程在Envi5.2中實現，具體操作：在ENVIToolbox中，選擇Toolbox/Radiometric Correction/ Radiometric Calibration，選擇*_MultiSpectral多光譜組（7個波段），打開輻射定標工具，對兩景影像分別做輻射定標。（2）影像鑲嵌因本文所使用的影像數據源是兩景Landsat OLI影像，因此需進行影像鑲嵌，鑲嵌的目的是將不同

11、的影像文件無縫地拼接成一幅完整的包含研究區(qū)域的影像。該處理過程在Envi5.2中實現，具體操作：在Toolbox中，選擇/Mosaicking/Seamless Mosaic，打開無縫鑲嵌工具，然后進行相關參數設置，如下所示：（3）影像裁剪因本文所使用的影像數據包含了北京市行政區(qū)劃以外的部分地區(qū)，因此需進行影像裁剪，以將研究區(qū)裁剪出來，并且減小了數據量，加快了數據處理速度，本文使用北京行政邊界矢量裁剪圖像。過程在Envi5.2中的具體操作如下： 在Toolbox中，選擇/Regions of Interest/Subset Data from ROIs，打開裁剪工具：影像裁剪結果如下所示：（4

12、）Flaash大氣校正電磁波在大氣中的傳輸和遙感器觀測過程中受光照條件以及大氣作用等的影響,只有小部分(在0.85um波段80%,在0.45um波段50%)太陽輻射能反射到遙感器,導致遙感器的測量值與地物實際的光譜輻射率不一樣。輻射損失主要發(fā)生在大氣吸收和散射過程，因此地表參數的遙感定量反演研究中,必須糾正目標輻射的不確定性信息6。ENVI中的FLAASH模型是基于MODTRAN4+輻射傳輸模型，通過參數查找表來進行大氣校正的商業(yè)化軟件。FLAASH大氣校正模塊支持多種傳感器數據，其光譜處理范圍0.4m-2.5m，可以有效地去除水蒸氣/氣溶膠散射效應，同時該方法基于圖像像素級的校正，能夠解決目

13、標像元和鄰近象元的“鄰近效應”問題，校正結果精度高，簡單易行7。然后，對大氣校正前后同一地物的光譜曲線進行對比，這里以植被為例，光譜曲線如下圖所示：校正前的植被光譜曲線 校正后的植被光譜曲線2、植被覆蓋度估算（1）計算NDVI 本文選取NDVI值為參數，采用像元二分模型對植被覆蓋度進行反演，根據植被覆蓋度的計算公式可知，要求取植被覆蓋度，首先需要計算NDVI。在Envi5.2中的具體操作如下：在Toolbox中，選擇Spectral/Vegetation/NDVI，NDVI Calculation Input File面板中，選擇LC8_rad_beijing_ref.dat 圖像，求算NDV

14、I，如下：NDVI求算結果如下：由于大氣校正后的結果有部分像元為負值，主要集中在陰影地區(qū)，這部分區(qū)域計算得到的NDVI在-1,1之外，為了便于后面的分析，我們這里統一將這部分像元進行處理，即NDVI值大于1的變?yōu)?，小于-1的變成-1。在Bandmath中的表達式為：-1>b1<1，其中b1代表的是NDVI，得到去除異常值文件：NDVI_去除異常值.dat（2） 掩膜文件制作該過程主要是為了計算NDVI的最大值、最小值所服務的，根據土地利用分類圖（共5類，林地、農業(yè)用地、城市用地、水體與其他）制作各種土地利用類型的掩膜文件，在Envi5.2中的具體操作如下：在Toolbox中選擇/

15、Raster Management/Masking/Apply Mask，打開制作掩膜工具：采用該方法制作林地、農業(yè)用地、城市用地、水體與其他的掩膜文件，其中林地與耕地的掩膜文件制作結果如下：（2）獲取閾值這一步就是求解NDVImax和NDVImin，使用上一步獲取的掩膜文件分別對NDVI圖像文件進行統計，在一定置信度范圍內獲取每個掩膜文件（也就是土地覆蓋類型）對應的最大和最小NDVI值。在Toolbox中，選擇/Statistics/Compute Statistics，進行統計，然后在統計結果中，取一定的置信度獲取最大和最小的NDVI值。如這里的林地覆蓋區(qū)域的統計結果（如下圖），這個過程帶

16、有很大的主觀性，我們需要根據統計學原理自己制定一套規(guī)則（比如5%的置信度），這里我就以NDVI值對應像元數量增加到5位數字為置信區(qū)間，選擇NDVImin=0.3804，NDVImax=0.8667。同樣的方法得到其他地物覆蓋類型的NDVI閾值，其中 ，水體沒有植被（水藻不屬于植被），認為這部分區(qū)域的植被覆蓋度為0，如下表：（3）生成參數文件植被覆蓋度的計算公式：fc =(NDVI - NDVIsoil ) (NDVIveg -NDVIsoil)，該過程是根據上面得到的NDVI閾值分別生成NDVIsoil和NDVIveg參數文件，也即NDVImin與NDVImax。該過程主要使用Envi5.2的

17、bandmath工具，并且：NDVIsoil：b1*0.3804+b2*0.1604+b3*0.0404+b4*0+b5*0.0946其中，b1：林地掩膜文件b2：農業(yè)用地掩膜文件b3：城市用地掩膜文件b4：水體掩膜文件b5：其他用地掩膜文件NDVIveg：b1*0.8667+b2*0.7794+b3*0.4789+b4*0+b5*0.5347其中，b1:林地掩膜文件b2：農業(yè)用地掩膜文件b3：城市用地掩膜文件b4：水體掩膜文件b5：其他用地掩膜文件最終，生成的參數文件如下所示：（4）植被覆蓋度估算利用上一步得到的NDVIsoil和NDVIveg參數文件帶入公式：fc =(NDVI - NDV

18、Isoil ) (NDVIveg -NDVIsoil)，該過程也是利用Envi5.2中的Bandmath工具來實現的，其表達式為：(b1-b2)/(b3-b2)，其中，b1為NDVI(對應的文件名為“NDVI_去除異常值.dat”)、B2為 NDVIsoil參數文件、B3:為NDVIveg參數文件，植被覆蓋度估算結果如下：我們分析下結果，會發(fā)現有一些異常值，即值在0,1之外，這些異常值是在NDVI置信度之外的那部分像元產生的（也包括NDVI異常像元）。這些像元數量不多，大約占3.7%左右。還有背景和水體區(qū)域的植被覆蓋度的值為-NaN，即無效值，因為分母為0造成的。第一種異常值可以將小于 0的值

19、變成0，大于1的值變成1，用 bandmath工具即可， Bandmath 表達式為： 0.0>b1<1.0，其中b1為植被覆蓋度 ；-NaN 可以用掩膜進行處理，即在Build Mask中用 -NaN生成掩膜。去掉異常值之后，并對其進行分類顯示，最終得到的植被覆蓋度圖如下：（5）結果驗證目前業(yè)內植被覆蓋度的驗證方法主要是通過野外調查結果進行驗證，如下為一種方法：以與地面垂直的角度用數碼相機拍攝采樣點的地面照片，使用GPS定位獲得采樣點的經緯度坐標。為獲得準確的植被覆蓋度實測數據減少像片邊緣變形誤差，將數碼相機得到的數字圖像截取長、寬各三分之二的中心地帶，采用非監(jiān)督分類法為10類，并將分類結果分為植被、非植被兩類，以此來計算出植被覆蓋度。為了保證驗證時能夠正確定位，野外一般選取3*3個像元大小的樣方，即9