多源多時相生態(tài)脆弱性遙感定量評價方法研究

楊明芬吳旭

（1.西藏自治區(qū)科技信息研究所，西藏 拉薩 850000；2.成都理工大學(xué)，四川 成都 610059）

生態(tài)環(huán)境是人類生存和發(fā)展的重要物質(zhì)基礎(chǔ)之一，對生態(tài)環(huán)境脆弱性變化展開研究，是制定區(qū)域可持續(xù)發(fā)展戰(zhàn)略的重要前提。對生態(tài)環(huán)境進行檢測和評價分析，有助于及時了解環(huán)境質(zhì)量變化的規(guī)律并采取科學(xué)的應(yīng)對措施。生態(tài)環(huán)境脆弱性評價是區(qū)域生態(tài)綜合評價的重要內(nèi)容和生態(tài)規(guī)劃的重要依據(jù)，近年來在地理學(xué)、生態(tài)學(xué)和環(huán)境科學(xué)研究中廣泛應(yīng)用。

以RS和GIS技術(shù)支持的生態(tài)環(huán)境脆弱性研究被更多的學(xué)者采用。生態(tài)環(huán)境脆弱性評價可分為定性評價和定量評價，定量方法對環(huán)境的脆弱性、穩(wěn)定性、敏感性和外部環(huán)境脅迫對系統(tǒng)造成的影響進行定量描述。［1］利用生態(tài)數(shù)學(xué)語言提出了生態(tài)脆弱帶基礎(chǔ)判定的數(shù)學(xué)方法。［2］用統(tǒng)計模型對生態(tài)環(huán)境脆弱性密度分布進行了初步測算。［3］通過多源遙感數(shù)據(jù)下植被覆蓋變化研究對脆弱性評價進行了研究。由于生態(tài)環(huán)境脆弱性研究涉及地學(xué)、生態(tài)學(xué)等多個領(lǐng)域的研究，具有很強的綜合性，因此傳統(tǒng)評價方法使用的研究方法、評價指標(biāo)體系和分級標(biāo)準(zhǔn)呈現(xiàn)多樣化。各種研究方法的重點集中在確定生態(tài)脆弱性評價指標(biāo)、權(quán)重設(shè)定、指定各指標(biāo)穩(wěn)定度的標(biāo)準(zhǔn)方面。這些方法需要較長的時間來處理前期數(shù)據(jù)資料，研究周期較長，無法快速定期對生態(tài)環(huán)境脆弱性現(xiàn)狀和變化趨勢進行評價監(jiān)測，而且評價結(jié)果不夠全面。

文章在前人建立的生態(tài)環(huán)境脆弱性定量評價方程基礎(chǔ)上，研究一種新的區(qū)域態(tài)環(huán)境脆弱性定量評價模型?；诙嘣炊鄷r相的衛(wèi)星遙感影像數(shù)據(jù)提取植被指數(shù)、土地覆蓋分類、DEM坡度與坡向數(shù)據(jù)作為評價模型輸入的原始數(shù)據(jù)。采用多元非線性回歸方法提取生態(tài)環(huán)境脆弱性定量評價所需的模型參數(shù)。通過評價模型的計算實現(xiàn)研究區(qū)生態(tài)環(huán)境脆弱性的快速、定量評價，形成脆弱性評價分級圖，并對模型評價的精度和穩(wěn)定性進行了驗證。

1 研究方法

1.1 基本思路

文章采用生態(tài)環(huán)境脆弱性的因子相關(guān)分析方法，參考［4］劉燕華和王經(jīng)民的經(jīng)驗公式,構(gòu)建脆弱性定量評價模型：

式（1）中F為生態(tài)環(huán)境脆弱性值，U為土地覆蓋類型的變化梯度，f為變化區(qū)域的植被覆蓋度，其中：

式（2）中U為變化的土地類型面積；Si為該類型區(qū)域前一年土地總面積為土地覆蓋類型變化梯度；i為土地覆蓋類型，參考前人研究基礎(chǔ)并結(jié)合與植被有關(guān)的特性，劃分為林地、草地、濕地、水體、人工用地、沙地等六種類型；

式（3）中f是估算的指標(biāo)覆蓋度，DNVI是歸一化植被指數(shù)，NDVImax和NDVImin分別是研究區(qū)同一年內(nèi)的最大值和最小值。

通過遙感監(jiān)督分類提取NDVI、土地覆蓋變化范圍，并通過編寫IDL程序計算土地變化的面積和總面積的最大值與最小值，求出變化梯度，根據(jù)NDVI求出植被覆蓋度，代入式（4）求出由于變化梯度導(dǎo)致的生態(tài)環(huán)境脆弱性值。在IDL環(huán)境中，對研究區(qū)的數(shù)字地形進行高程、坡度與坡位的綜合分析，選取與植被相關(guān)性最大的因素代入評價方程中，計算得出脆弱性評價結(jié)果。將數(shù)字地形模型與評價結(jié)果疊加，得出區(qū)域生態(tài)環(huán)境脆弱性綜合評價等級。實驗流程如圖1。

圖1 實驗流程圖

1.2 實驗數(shù)據(jù)

研究區(qū)位于四川省若爾蓋縣，四川、甘肅、青海三省交界處，地處青藏高原東部邊緣，海拔3400-3900m，東西與南北最大距離約150km，面積約10620km2，地理坐標(biāo)為 102°08′E 至 103°39′E，32°56′至34°19′N之間，是長江和黃河流域上游地區(qū)重要的水源涵養(yǎng)地，也是中國殘存面積最大的高原泥炭沼澤濕地，整個研究區(qū)形成一個相對完整的地理自然區(qū)，對于研究區(qū)域生態(tài)脆弱性變化極具代表性。

1.2.1 Landsat數(shù)據(jù)。Landsat1-8系列衛(wèi)星是由美國航空航天局 NASA（National Aeronautics and Space Ad?ministration）發(fā)射的對地觀測衛(wèi)星。文章采用Land?sat5 TM數(shù)據(jù)提取土地覆蓋分類信息，并計算土地類型變化率。TM的覆蓋范圍為之間的所有陸地區(qū)域，東西掃描范圍約183km，南北掃描范圍約為170km，數(shù)據(jù)更新周期為16天，空間分辨率30m。文章獲取2001-2010年研究區(qū)內(nèi)的Landsat5 TM遙感數(shù)據(jù)（每年4景，共40景數(shù)據(jù)）?？紤]到研究每一景影像可能存在多云覆蓋的問題，我們選擇的單景影像平均云量小于10%的影像數(shù)據(jù)，成像時間集中在8-9月。選擇出每一年內(nèi)成像質(zhì)量較好的影像數(shù)據(jù)，對各景數(shù)據(jù)進行鑲嵌，根據(jù)研究區(qū)矢量地圖裁剪得出研究區(qū)影像數(shù)據(jù)，對裁剪的研究區(qū)數(shù)據(jù)進行大氣校正、正射校正、去除薄云處理。其中投影為UTM、坐標(biāo)系為WGS-1984，如圖2。

圖2 研究區(qū)歷年影像

1.2.2 MODIS數(shù)據(jù)。文章采用MOD13A1 NDVI數(shù)據(jù)產(chǎn)品，共獲取2001、2003、2006、2008年4個年份，16天間隔的產(chǎn)品數(shù)據(jù)（每年23景，共計92景影像數(shù)據(jù)，空間分辨率500m）。在ENVI軟件環(huán)境中對每年數(shù)據(jù)及十年數(shù)據(jù)進行HANTS濾波處理，生成長時空序列數(shù)據(jù)集，采用波段運算工具，對NDVI的月均、年均值進行計算。對數(shù)據(jù)集進行Hants濾波處理，然后利用波段運算工具計算2001-2010年均NDVI的數(shù)據(jù)集，如圖3。

圖3 NDVI年均數(shù)據(jù)集

1.2.3 地形數(shù)據(jù)。本文選取的DEM數(shù)據(jù)為SRTMDE?MUTM（90m分辨率），提取高程數(shù)據(jù)；坡度數(shù)據(jù)為SRTMSLOPE（90m分辨率）；坡向數(shù)據(jù)為SRTMAS?PECT（90m分辨率），如圖4所示。這三種數(shù)據(jù)用于分析其與NDVI的相關(guān)性，估算未知區(qū)域NDVI的取值，為定量評價模型提供所需參數(shù)。

圖4 地形數(shù)據(jù)

1.3 實驗步驟

1.3.1 計算土地覆蓋類型的變化梯度。研究區(qū)典型地物類型包括：河流、湖泊、沼澤、沼澤化草甸、濕草甸、灌叢、灌木林、森林、裸地、裸巖、沙地和建筑用地等［5］。其中河流、湖泊、泥炭沼澤、沼澤化草甸和濕草甸是研究區(qū)的主要濕地類型［7］。文章參考《全國遙感監(jiān)測土地利用/覆蓋分類體系》、前人研究成果，并考慮模型計算的簡便性，把土地覆蓋類型分為：草地、林地、裸地（裸土、沙地、裸巖、積雪）、水體（河流、湖泊）、沼澤、人工用地（居民地、耕地、道路）六類［6］。為保證分類精度和信息提取的快捷，文章采用ENVI軟件下的支持向量機分類（Support Vector Machine，SVM）方法對遙感圖像進行監(jiān)督分類，再通過人工操作修改錯誤分類。分類過程中結(jié)合野外調(diào)查的先驗知識選擇訓(xùn)練樣本。影像中難免存在不可去除的云層干擾等因素，由此產(chǎn)生錯誤分類，文章通過與影像觀測時間最接近的數(shù)據(jù)進行比較人工修改分類。把相同年份的遙感預(yù)處理數(shù)據(jù)與分類結(jié)果數(shù)據(jù)疊加，在遙感預(yù)處理影像中隨機選取各類型純凈像元，采用ENVI的分類后處理工具中的誤差混合矩陣進行分類精度評價，分類結(jié)果精度如表1所示（以2001年為例）。

表1 土地覆蓋類型分類精度表

從各時相分類精度可以看出，總體精度在85%以上，可以滿足應(yīng)用需求。其中沼澤與水體、沼澤與草地、裸地與草地誤分的現(xiàn)象明顯，這主要是由于高植被覆蓋沼澤與草地的光譜信息比較接近，地覆蓋草地與裸露地表、沙地等裸地在圖像上比較相似，從而導(dǎo)致誤分；同時少量區(qū)域由于厚云根本上無法提出，對林地、草地也存在誤分。

文章在獲取2001、2003、2006、2008年土地覆蓋類型的分類數(shù)據(jù)基礎(chǔ)上，采用的分類統(tǒng)計工具計算每種覆蓋類型相鄰年份的變化值，計算出土地覆蓋類型的變化梯度，對土地覆蓋類型的動態(tài)變化進行分析。通過計算結(jié)果獲得變化的類型、面積、像元數(shù)量，并統(tǒng)計出土地覆蓋類型變化的百分比。

1.3.2 植被覆蓋度估算。然后基于多源遙感數(shù)據(jù)的生態(tài)環(huán)境脆弱性評價，評價模型參數(shù)組要依賴于植被覆蓋度和土地覆蓋類型的變化梯度。但一個地區(qū)的生態(tài)環(huán)境脆弱性往往是多種因素共同的結(jié)果。研究植被指數(shù)、地形高度、坡度以及坡向變化的規(guī)律有利于區(qū)域生態(tài)環(huán)境脆弱性定量評價的準(zhǔn)確性。文章對獲取的2001-2010年NDVI數(shù)據(jù)集與研究區(qū)地形數(shù)據(jù)的相關(guān)性分析，步驟如下：

將NDVI數(shù)據(jù)集與研究區(qū)高程、坡度、坡向數(shù)據(jù)疊加，并配準(zhǔn)；

編寫程序，在研究區(qū)內(nèi)指定的矩形窗口范圍內(nèi)生成隨機采樣點，獲取每一年份同一經(jīng)緯度的NDVI值、海拔高程、坡度值、坡向值；

通過程序自動提取對應(yīng)坐標(biāo)點的植被指數(shù)、值高程值、坡度和坡向值，生成NDVI-高程、NDVI-坡度、NDVI-坡向的數(shù)據(jù)記錄文件，以便進一步分析它們之間的關(guān)系。

1.3.3 NDVI與地形數(shù)據(jù)的相關(guān)性分析。將采集到的10個時相NDVI數(shù)據(jù)與地形高程進行研究：利用IDL編程繪制出對應(yīng)坐標(biāo)點的NDVI隨地形高度變化的散點圖，對每一時相2000多個采樣點進行統(tǒng)計分析。從繪制出的散點圖可以觀察出：NDVI值隨著高度的增加，是先增加后減少，其變化規(guī)律和植被分布的垂直地帶性是相符合。坡度為0的時候，ND?VI的值在坡度0.4-0.5區(qū)間分布均勻，這表明在草地這樣的平坦地形，NDVI的值是沒有反映坡度的變化規(guī)律。因此NDVI隨坡度的變化表現(xiàn)平直，沒有明顯的升降變化。同時，NDVI與坡向的變化分布均勻，沒有明顯的規(guī)律性變化，數(shù)據(jù)與各種擬合函數(shù)相關(guān)性極低。NDVI與坡度和坡向的相關(guān)系數(shù)大都小于0.2，因此跟研究中只針對NDVI與高度值進行回歸分析和函數(shù)擬合。根據(jù)散點圖顯示，NDVI與高度之間有明顯的規(guī)律性變化，但無明顯的線性或非線性函數(shù)關(guān)系，考慮對NDVI值進行某種轉(zhuǎn)換，研究分析轉(zhuǎn)換后的值與高程之間的相關(guān)關(guān)系，如表2所示。

表2 NDVI與地形數(shù)據(jù)的相關(guān)性分析

文章嘗試對NDVI與NDVI均值差的絕對值取自然對數(shù)，然后分析轉(zhuǎn)換后的值與高程的關(guān)系、轉(zhuǎn)換值，通過編寫程序生成與高程值的散點圖，如圖5所示。根據(jù)散點圖的分布情況進一步分析擬合函數(shù)的形態(tài)。

圖5 NDVI與高程的非線性分析

1.3.4 基于多元非線性回歸的NDVI估值。對多年份數(shù)據(jù)的轉(zhuǎn)換值Y與高度值的散點分布進行分析，發(fā)現(xiàn)散點的分布特征十分接近e的冪函數(shù)類型，因此設(shè)計一種多元非線性回歸方程對數(shù)據(jù)點進行擬合，并通過IDL編程計算回歸參數(shù)：Y=C*eax+b，其中x為高度值，a、b、c為待定參數(shù)。通過以上計算阿對未知點的NDVI進行估值，從而估算出研究區(qū)的植被覆蓋度。圖6描述了計算的待定參數(shù)形成的擬合函數(shù)與Y的關(guān)系。

擬合函數(shù)對NDVI估計值的精度檢驗如表2所示。

表2 NDVI估值精度檢驗

圖6 NDVI非線性函數(shù)擬合

1.4 定量評價模型

非線性回歸模型表示為：

其中f為期望函數(shù)，xn為第n個響應(yīng)的回歸向量或自變量向量。模型的期望響應(yīng)是關(guān)于參數(shù)的非線性函數(shù)。θ表示未知參數(shù)，當(dāng)分析一個特定的數(shù)據(jù)集時，假定向量xn，，n=1，…N，N是固定的，并幾種考慮期望響應(yīng)關(guān)于θ的相依關(guān)系。構(gòu)造N維向量η（θ），其第n個向量為：ηn（θ）=f（xn，θ），n=1，2…，N；記非線性回歸模型為：

對于非線性的情形，利用迭代的方法求解最小二乘估計量。利用期望函數(shù)的線性近似，通過迭代改善的初始值θ0，直到迭代無變化為止。

將期望函數(shù)f（xn，θ）在θ0處進行一階泰勒展開：

與η（θ0）相比，現(xiàn)在點y。因此取更好的參數(shù)值θ1=θ0+δ0，開始另一次迭代，繼續(xù)計算新的殘差z1=y-η（θ），新的導(dǎo)數(shù)矩陣V1和新的增量。重復(fù)上述過程，直至達到收斂為止。

2 實驗結(jié)果與分析

在生態(tài)環(huán)境評價中，由于評價單元不同，評價結(jié)果會存在差異,因此評價單元的選擇與研究區(qū)范圍密切相關(guān)［7］。文章以遙感影像中的柵格數(shù)據(jù)為載體，以每一格柵格單元為評價的最小單位，在模型計算的基礎(chǔ)上，按照指數(shù)的數(shù)值大小進行等級劃分，確定其優(yōu)劣程度。把脆弱性指標(biāo)值劃分為5個等級：無脆弱性、輕微脆弱性、一般脆弱性、中度脆弱性和嚴(yán)重脆弱性。文章采用此方法對2003、2006、2008三年的評價結(jié)果進行分級，以此分級閾值作為評價，最終得到研究區(qū)的生態(tài)環(huán)境脆弱性定量評價結(jié)果，如圖7所示。

圖7 區(qū)域生態(tài)環(huán)境脆弱性定量評價計算結(jié)果

由計算結(jié)果可知研究區(qū)評價值主要集中在0.2-0.8之間，生態(tài)環(huán)境脆弱性嚴(yán)重的區(qū)域主要集中在北部和中部部分地方。這里主要是裸露地表、沼澤邊緣和建筑用地、耕地等受人類活動影響較為頻繁的區(qū)域。從圖7可以看出，各評價等級相間分布，各等級間沒有明顯的界線。

3 結(jié)論

文章從少量因子出發(fā)，先抽取和植被生態(tài)方面有關(guān)的因子（土地覆蓋類型變化率、植被覆蓋度），研究生態(tài)環(huán)境脆弱性和它們之間的關(guān)系。再研究植被指數(shù)和地形變化的規(guī)律，逐個因子添加不僅可以使研究的復(fù)雜問題變的簡單而且有助于揭示它們之間的聯(lián)系。主要研究成果包括：實現(xiàn)了生態(tài)環(huán)境脆弱性的遙感定量研究。解決了如何從遙感影像中提取出反映生態(tài)環(huán)境脆弱性的因子，并把這種參數(shù)帶入到生態(tài)環(huán)境脆弱性的方程中，從而實現(xiàn)生態(tài)環(huán)境脆弱性的遙感定量研究是問題的關(guān)鍵。從植被的角度出發(fā)把造成生態(tài)環(huán)境脆弱性的因子劃分為：植被覆蓋度、土地覆蓋類型變化梯度、地形高度、坡度以及坡向。從遙感影像中提取出土地覆蓋類型變化的梯度，利用NDVI和植被覆蓋度的關(guān)系，定量表達植被覆蓋度。通過模型求出生態(tài)環(huán)境脆弱度，實現(xiàn)生態(tài)環(huán)境脆弱性的遙感定量研究，研究具有創(chuàng)新意義。揭示了研究區(qū)植被指數(shù)隨地形高度、坡度、坡向變化的規(guī)律。研究中對每一時相采集到2000多個具有典型意義的點，數(shù)據(jù)豐富，資料可靠。通過研究發(fā)現(xiàn)，NDVI和高程存在著明顯的相關(guān)關(guān)系，通過IDL編程求出了具體的定量表達關(guān)系式。實驗表明植被指數(shù)隨地形坡度和坡向的變化規(guī)律不明顯，只簡單地分布在均值周圍?；诙嘣蔷€性回歸方法，建立了生態(tài)環(huán)境脆弱性植被遙感定量研究的模型，通過改造傳統(tǒng)意義上的脆弱度方程，實現(xiàn)了參數(shù)的替代。該模型中脆弱度被定義為土地覆蓋類型變化的梯度和植被覆蓋度的商。通過反映生態(tài)環(huán)境脆弱性本底值的方法來評價劃分生態(tài)環(huán)境脆弱性等級。