劉朔, 楊建勇, 蔡凡隆

四川省林業(yè)和草原調(diào)查規(guī)劃院，四川 成都 610081

隨著遙感技術(shù)的發(fā)展，植被覆蓋研究成了分析生態(tài)環(huán)境、氣候演變和水文過程等眾多領(lǐng)域的基礎(chǔ)[1]。歸一化差值植被指數(shù)（NDVI）是反映植被覆蓋的一個重要參數(shù)，具有較強的植被檢測能力，正值表示有植被覆蓋，且數(shù)值隨覆蓋度增大而增大，NDVI不僅能更好地檢測植被的生長狀況，并且能消除部分輻射誤差[2-5]。像元二分模型是一種實用的遙感估算模型，是采用NDVI歸一化植被指數(shù)定量估算植被覆蓋度的模型，在植被覆蓋度研究和生態(tài)研究及評價等方面應(yīng)用廣泛，許多學(xué)者基于 NDVI 數(shù)據(jù)和像元二分模型做了多方面研究[6-10]，一般通過NDVI、植被覆蓋度與植被相關(guān)的指標(biāo)分析植被覆蓋變化和趨勢。國內(nèi)學(xué)者基于NDVI對毛烏素沙地、新疆、西藏高原等沙化廣布區(qū)進行了分析研究[11-14]，基礎(chǔ)數(shù)據(jù)主要為MODIS_NDVI 數(shù)據(jù)，該數(shù)據(jù)空間分辨率為250 m，適合大尺度、廣區(qū)域的研究區(qū)。基于NDVI對川西北沙區(qū)的研究鮮有報道，重要原因之一是川西北沙區(qū)沙化土地分布分散，主要呈小塊狀，點狀分布，MODIS_NDVI 數(shù)據(jù)精度難以滿足對該區(qū)域沙化土地的分析需求。

若爾蓋縣是川西北沙化最嚴峻的地區(qū)，掌握若爾蓋沙地動態(tài)變化趨勢有利于該區(qū)域生態(tài)環(huán)境保護與恢復(fù)，為該區(qū)域防沙治沙提供理論依據(jù)與數(shù)據(jù)支持。本文以若爾蓋縣第五次沙化監(jiān)測圖斑區(qū)為研究區(qū)，以2013年、2015年、2017年3個年度10月初的 Landsat 8影像作為年際對比基礎(chǔ)數(shù)據(jù)，主要原因是10月初植被已處于生長末期，更有利于提取低植被區(qū)沙化土地，減少將沙化土地植被過高判讀的可能，加之研究區(qū)以草地植被為主，植被單一，季節(jié)間生長規(guī)律同步性好。通過Landsat 8影像中的近紅外波段和紅光波段數(shù)據(jù)提取空間分辨率為30 m的NDVI，再利用像元二分模型反演植被覆蓋度，分析研究區(qū)植被變化特征，評估研究時段期間植被變化趨勢，為若爾蓋縣沙化治理和生態(tài)環(huán)境建設(shè)對策制定提供科學(xué)依據(jù)，應(yīng)用的技術(shù)方法可為川西北沙化監(jiān)測、沙區(qū)植被動態(tài)變化對比和評估提供科技支撐。

1 研究區(qū)概況

研究區(qū)為若爾蓋縣全國第五次沙化監(jiān)測的圖斑區(qū)，總面積79 851.78 hm2，主要分布在牧區(qū)，以草地沙化為主,平均海拔3 500 m，地質(zhì)母巖多為三疊系砂巖、板巖、灰?guī)r和第四紀的松散堆積物，其形成的土壤含沙量高。研究區(qū)屬大陸性季風(fēng)高原氣候區(qū)，氣候寒冷干燥，長冬無夏，日照強烈，溫差大，年均氣溫0.7℃ ，無絕對無霜期，年均日照2 389 h；年均降雨量656.8 mm，年均蒸發(fā)量1232 mm，年均濕度69%；平均風(fēng)速2.4 m·s-1，最大風(fēng)速40 m·s-1，年大風(fēng)日數(shù)50 d，年沙塵日數(shù)6 d。

2 研究方法

2.1 數(shù)據(jù)來源及處理

（1）研究采用Landsat 8遙感影像進行NDVI提取。Landsat 8遙感影像來自“地理空間數(shù)據(jù)云”網(wǎng)站，分別是2013年10月、2015年10月、2017年10月的Landsat 8遙感影像各1景（見表1），空間分辨率為30 m（僅全色波段為15 m）；（2）研究區(qū)范圍數(shù)據(jù)來源若爾蓋縣第五次全國沙化監(jiān)測數(shù)據(jù)。

表 1 Landsat 8影像主要信息表Tab. 1 Main information table of Landsat 8 image

2.2 數(shù)據(jù)處理及研究方法

數(shù)據(jù)總體處理流程為：（1）在ENVI5.3中主要完成NDVI和FC（植被覆蓋度）提取，主要步驟分別為：①輻射定標(biāo)；②大氣校正；③NDVI計算；④FC（采用像元二分模型反演植被覆蓋度）。（2）在Arcgis10.2中主要進行NDVI和FC的分級、植被變化的趨勢分析、轉(zhuǎn)移矩陣計算。主要步驟分別為：①通過無偏移天地圖衛(wèi)星影像將研究區(qū)矢量圖層、NDVI柵格、FC柵格校正為大地2000投影坐標(biāo)系；②按掩膜提取研究區(qū)的NDVI和FC柵格；③進行NDVI和FC的分級；④植被變化轉(zhuǎn)移矩陣計算；⑤在研究區(qū)創(chuàng)建生成1 000個隨機點矢量數(shù)據(jù)，用來提取各年度NDVI和FC的柵格數(shù)據(jù)值，并導(dǎo)出。（3）在Spss24中對導(dǎo)出的研究區(qū)1 000個矢量點導(dǎo)出數(shù)據(jù)進行數(shù)據(jù)類型分布探索的基礎(chǔ)上，進行各年度間NDVI和FC差異分析，以及與沙化程度的相關(guān)分析。經(jīng)Spss24數(shù)據(jù)探索，整體數(shù)據(jù)呈偏態(tài)分布，非正態(tài)分布（這是研究區(qū)域數(shù)據(jù)本身決定的，因為植被占絕大部分，而非植被的裸地等占比很小），因此選用Spss24 中“非參數(shù)檢驗”中“k個相關(guān)樣本”進行差異分析。

主要研究方法如下：

（1）NDVI提取

NDVI 是遙感監(jiān)測中反映植被在可見光、近紅外波段與土壤背景之間光學(xué)差異的指標(biāo)，是對地表植被活動的簡單、有效和經(jīng)驗的度量[15]。計算公式為：

NDVI=（Rnir - Rr ） /（Rnir + Rr）

式中：Rnir為近紅外波段反射率； Rr為紅波段反射率。

在ENVI5.3中必須先對影像進行輻射定標(biāo)和大氣校正，再通過 “band math” 計算工具中輸入公式或者使用軟件中專門的NDVI工具進行求算。在Arcgis10.2中根據(jù)NDVI數(shù)據(jù)取值及分布特征并結(jié)合經(jīng)驗成果，將其重分類為4個等級，分別為：低（-1～0.2）；較低（0.2～0.4）；中等（0.4～0.5）；較高（0.5～1）。

（2）反演植被蓋度

植被覆蓋度(vegetation fractional cover)簡稱FC，在眾多遙感測量植被覆蓋度的方法中，較為實用和普遍的方法是采用像元二分模型利用植被指數(shù)近似估算植被覆蓋度，通過一系列的模型運算[16]，可以得到：

式中：NDVIsoil為完全是裸土或無植被覆蓋區(qū)域的 NDVI 值，NDVIveg代表完全被植被所覆蓋像元的NDVI 值，即純植被像元的 NDVI 值。

受眾多因素影響， NDVIsoil并不是理論上的零值，通常在 -1～0.2 之間變化，NDVIveg也隨植被類型以及時間和空間而所差異，因此NDVIsoil和NDVIveg沒有固定值[17]。實際應(yīng)用中，在沒有實測數(shù)據(jù)的情況下，可根據(jù)經(jīng)驗取一定置信度范圍內(nèi)的NDVI max 和 NDVI min 作為和DVIveg和NDVIsoil則有：

FC在ENVI5.3中可通過 “band math” 計算工具中輸入公式進行求算，本研究通過多次測算和對比，將NDVI為-1～0之間的累計頻率作為NDVI min，將累計頻率為 90% 的 NDVI作為NDVI max。并將FC數(shù)據(jù)取值及分布特征并結(jié)合前人經(jīng)驗，將其重分類為4個等級：低（0～0.2）；較低（0.2～0.4）；中等（0.4～0.5）；較高（0.5～1）

（3）轉(zhuǎn)移矩陣分析

Arcgis10.2中將各期NDVI和FC柵格數(shù)據(jù)轉(zhuǎn)化為矢量數(shù)據(jù)，再分別進行疊加相交分析后導(dǎo)出數(shù)據(jù)，并在EXCEL中通過數(shù)據(jù)透視表整理提取面積轉(zhuǎn)移矩陣。

3 結(jié)果與分析

3.1 NDVI及FC年總體變化特征

由圖1和圖2和表2可知，研究區(qū)NDVI和FC從2013年、2015年、2017年均呈顯著增加（差異均達到極顯著），NDVI中位數(shù)值從2013年的0.477 5，增加到2015年0.537 4，至2017年達0.592 1；FC中位數(shù)值從2013年的0.530 5，增加到2015年0.597 1，至2017年達0.657 8，說明研究區(qū)植被覆蓋度明顯增加。主要原因是2007年起若爾蓋縣逐步開始了持續(xù)沙化治理工程，包括：省級沙化治理試點工程以及后續(xù)的成果鞏固工程、 2013年開始實施的國家川西藏區(qū)沙化治理工程、國家重點生態(tài)功能區(qū)轉(zhuǎn)移支付資金項目等，工程中相關(guān)圍欄封育、植灌種草、人工種草、生物沙障等措施成效顯著[18-20]。同時，近年來若爾蓋縣持續(xù)開展畜牧業(yè)提檔升級，推行畜種改良，以草定畜等，草地超載情況得到明顯改善，2012年末全縣牲畜存欄114.1萬混合頭[21]，到2017年降低為95.96萬混合頭[22]。

圖 1 NDVI各年度中位數(shù)變化圖Fig. 1 Annual median change chart of NDVI

圖 2 FC各年度中位數(shù)變化圖Fig. 2 Annual median change chart of FC

圖 3 NDVI 變化趨勢圖Fig. 3 NDVI change trend

圖 4 FC變化趨勢圖Fig. 4 FC change trend

3.2 NDVI及FC動態(tài)變化轉(zhuǎn)移矩陣

由表3、表4、表5可知，研究區(qū)分級后的NDVI從2013年、2015年、2017年均呈顯著增加。其中：2013年至2015年，較高等級面積占比從7.26%提高到30.69%；中等級面積占比從75.97%降低到59.31%；較低等級面積占比從13.56%降低到7.32%；低等級面積占比從3.2%降低到2.67%。其具體轉(zhuǎn)化方向主要是：中等級向較高等級轉(zhuǎn)化，轉(zhuǎn)化面積占自身面積31.9%；中等級向較高等級轉(zhuǎn)化，轉(zhuǎn)化面積占自身面積56.9%；低等級向較低等級轉(zhuǎn)化，轉(zhuǎn)化面積占自身面積30.2%。

表 2 研究區(qū)各年度NDVI和FC差異顯著性系數(shù)（Bonferroni調(diào)整）Tab. 2 Significant coefficient of NDVI and FC in each year in the study area (Bonferroni adjustment)

表 3 研究區(qū)NDVI和FC分級差異顯著性系數(shù)（Bonferroni調(diào)整）Tab. 3 Significant coefficient of NDVI and FC classification grades in the study area (Bonferroni adjustment)

表 4 研究區(qū)2013—2015年NDVI動態(tài)變化矩陣Tab. 4 Dynamic change matrix of NDVI in the study area from 2013 to 2015

表 5 研究區(qū)2015—2017年NDVI動態(tài)變化矩陣Tab. 5 Dynamic change matrix of NDVI in the study area from 2015 to 2017

2015年至2017年，較高等級面積占比從30.69%提高到46.22%；中等級面積占比從59.31%降低到45.87%；較低等級面積占比從7.32%降低到4.94%；低等級面積占比從2.67%增加到2.96%。其具體轉(zhuǎn)化方向主要是：中等級向較高等級轉(zhuǎn)化，轉(zhuǎn)化面積占自身面積34.05%；中等級向較高等級轉(zhuǎn)化，轉(zhuǎn)化面積占自身面積49.39%；低等級向較低等級轉(zhuǎn)化，轉(zhuǎn)化面積占自身面積27.88%。

由表3、表6、表7可知，研究區(qū)分級后的FC從2013年、2015年、2017年均呈顯著增加。其中：2013年至2015年，較高等級面積占比從64.92%提高到82.14%；中等級面積占比從25.45%降低到10.71%；較低等級面積占比從6.88%降低到4.85%；低等級面積占比從2.74%降低到2.29%。其具體轉(zhuǎn)化方向主要是：中等級向較高等級轉(zhuǎn)化，轉(zhuǎn)化面積占自身面積71.41%；中等級向較高等級轉(zhuǎn)化，轉(zhuǎn)化面積占自身面積30.08%；低等級向較低等級轉(zhuǎn)化，轉(zhuǎn)化面積占自身面積29.42%。

表 6 研究區(qū)2013—2015年FC動態(tài)變化矩陣Tab. 6 Dynamic change matrix of FC in the study area from 2013 to 2015

表 7 研究區(qū)2015—2017年FC動態(tài)變化矩陣Tab. 7 Dynamic change matrix of FC in the study area from 2015 to 2017

2015年至2017年，較高等級面積占比從82.14%提高到89.41%；中等級面積占比從10.71%降低到4.07%；較低等級面積占比從4.85%降低到3.92%；低等級面積占比從2.29%增加到2. 60%。其具體轉(zhuǎn)化方向主要是：中等級向較高等級轉(zhuǎn)化，轉(zhuǎn)化面積占自身面積77.41%；中等級向較高等級和高等級轉(zhuǎn)化，轉(zhuǎn)化面積分別占自身面積23.24%和25.03%；低等級向較低等級轉(zhuǎn)化，轉(zhuǎn)化面積占自身面積28.31%。

3.3 NDVI及FC與沙化程度相關(guān)性分析

將研究區(qū)NDVI 及FC等級與研究區(qū)2015年監(jiān)測的沙化程度（輕度、中度、重度、極重度）進行相關(guān)分析，發(fā)現(xiàn)它們均呈顯著正相關(guān)（見表8）。主要原因是沙化程度的劃分主要依據(jù)之一就是植被覆蓋度，而NDVI和FC均是反映植被生長情況的指標(biāo)。另一方面也說明，本次NDVI 及FC重分類等級的劃分基本符合研究區(qū)沙化植被覆蓋度現(xiàn)狀。由于全國沙化監(jiān)測工作是每5年開展一次，時間跨度長，利用NDVI及其反演的植被覆蓋度的方法理論上可在不同時段（如年際、月際），對研究區(qū)沙化植被動態(tài)變化進行分析評估，是一種非常靈活高效的方法。

表 8 研究區(qū)NDVI 及FC等級與沙化程度相關(guān)性Tab. 8 Correlation between NDVI, FC grades and the desertification degree in the study area

4 討論

（1）Landsat 8作為免費的影像數(shù)據(jù)，自2013年來已被廣泛應(yīng)用于自然資源調(diào)查、生態(tài)環(huán)境監(jiān)測、地質(zhì)災(zāi)害響應(yīng)、教育科研和政府管理等領(lǐng)域，其空間分辨率為30 m，較MODIS影像（空間分辨率為250 m）高很多，可進行更加細致和精確的研究。但是，其拍攝周期間隔16 d，相對較長，因此影像可選擇的數(shù)量有限，加之?dāng)?shù)據(jù)若受到氣候、天氣等因素的影響，云量較大的現(xiàn)象也較普遍。筆者搜索了研究區(qū)2013—2020年4—10月的Landsat 8影像，再排除受云量影響較大影像后，結(jié)合研究內(nèi)容特點，最終篩選確定了2013年、2015年、2017年共3個年度10月的影像數(shù)據(jù)進行研究。不同于MODIS可直接下載合成好的NDVI數(shù)據(jù)，Landsat 8需要根據(jù)紅光波段和近紅外波段求算NDVI值，在此過程中，輻射定標(biāo)和大氣校正是必須進行預(yù)處理工作，這是基于Landsat 8準確求算NDVI的前提，未進行大氣校正測算的NDVI值要普遍小于大氣校正后測算的NDVI值。

（2）研究區(qū)2013年、2015年、2017年的NDVI和FC均呈顯著增加，由于NDVI和FC指標(biāo)變化趨勢一致，都可反映植被蓋度變化趨勢，實際生產(chǎn)過程中可根據(jù)研究需要靈活選擇指標(biāo)進行研究分析。

（3）在實際生產(chǎn)工作中往往需要對植被蓋度進行分級，確定各等級間的區(qū)劃閾值最為關(guān)鍵，但區(qū)劃等級間的閾值沒有統(tǒng)一標(biāo)準，要綜合影像、研究目的、研究方法和研究區(qū)本底情況或結(jié)合已有研究而確定。本研究中，研究區(qū)NDVI 及FC等級與研究區(qū)2015年監(jiān)測的沙化程度（輕度、中度、重度、極重度）均呈顯著正相關(guān)，說明本次NDVI 及FC重分類等級的劃分基本符合研究區(qū)沙化植被覆蓋度現(xiàn)狀。在沙化監(jiān)測及成效對比評估工作中，充分利用遙感數(shù)據(jù)的同步性、時效性、周期性和成本優(yōu)勢，深入研究NDVI與工作區(qū)植被蓋度關(guān)聯(lián)特征，將有效提升監(jiān)測、評估工作的時效，降低因調(diào)查時間差異或調(diào)查者個體差異產(chǎn)生的判讀誤差，為若爾蓋縣沙化治理和生態(tài)環(huán)境建設(shè)對策制定提供科學(xué)依據(jù)。