NDVI時(shí)間變換一致性的貴州茂蘭植被覆蓋變化分析

陳炫熾，田鵬舉，陳 蓉，吳愈鋒，王躍躍

(1.貴州大學(xué)礦業(yè)學(xué)院，貴州 貴陽 550025；2.貴州省生態(tài)氣象和衛(wèi)星遙感中心，貴州 貴陽 550025)

植被是陸地生態(tài)系統(tǒng)的重要組成部分，是連接土壤、大氣和水分的紐帶，同時(shí)也是生態(tài)環(huán)境的一個(gè)重要影響因素。當(dāng)前，生態(tài)文明建設(shè)在全國(guó)各地廣泛開展，對(duì)貴州茂蘭植被覆蓋變化的研究，有助于掌握當(dāng)?shù)厣鷳B(tài)環(huán)境變化，加快推進(jìn)生態(tài)文明建設(shè)。

近年來，隨著科學(xué)技術(shù)的發(fā)展，遙感在植被覆蓋變化研究領(lǐng)域的應(yīng)用已不再受限于數(shù)據(jù)、資料來源匱乏、分析角度和處理手段單一的窘境：文獻(xiàn)[1]以Hyperion高光譜影像為數(shù)據(jù)源，利用混合像元分解估測(cè)了喀斯特地區(qū)植被覆蓋度；文獻(xiàn)[2]選取近30年中3時(shí)期Landsat TM/ETM/OLI數(shù)據(jù)，分析了梵凈山植被覆蓋時(shí)空變化特征和影響因素；文獻(xiàn)[3]利用MODIS NDVI數(shù)據(jù)，從時(shí)序演化、數(shù)量轉(zhuǎn)移和空間演化3個(gè)方面分析了淮南礦區(qū)植被覆蓋的時(shí)空演化特征。綜上，在GIMMS-NDVISPOT VGT/NDVIMODIS-NDVI等中低空間分辨率植被指數(shù)產(chǎn)品的支持下，大區(qū)域尺度的植被覆蓋變化研究已趨于成熟且成果斐然；但小區(qū)域尺度下長(zhǎng)時(shí)間序列的植被覆蓋變化研究因可用影像數(shù)據(jù)稀少、云覆蓋對(duì)影像質(zhì)量的影響等問題，難以構(gòu)建較為完整的植被覆蓋時(shí)間序列并進(jìn)行分析，故目前此方面的研究進(jìn)展大為受限。

因此，本文基于Landsat TM/ETM/OLI數(shù)據(jù)，以具有重要研究?jī)r(jià)值的貴州茂蘭國(guó)家級(jí)自然保護(hù)區(qū)為研究區(qū)，選擇已被廣泛用于表征地表植被覆蓋的歸一化植被指數(shù)(normalized difference vegetation index,NDVI)作為反映植被覆蓋變化的參數(shù)，使用基于NDVI時(shí)間變換一致性的方法消除云覆蓋的影響，建立較為完整的NDVI時(shí)間序列。在此基礎(chǔ)上通過一元線性回歸模型和相關(guān)分析法研究氣候?qū)ρ芯繀^(qū)植被覆蓋變化的影響，為增強(qiáng)生態(tài)文明建設(shè)對(duì)氣候變化的適應(yīng)能力提供科學(xué)依據(jù)，也為小區(qū)域尺度下的植被覆蓋研究提供技術(shù)支持。

1 研究區(qū)域與研究方法

1.1 研究區(qū)地理位置

貴州茂蘭國(guó)家級(jí)自然保護(hù)區(qū)位于貴州省黔南布依族苗族自治州荔波縣東南部，南與廣西壯族自治區(qū)接壤，毗鄰廣西木倫國(guó)家自然保護(hù)區(qū)，地理位置為東經(jīng)107°52′～108°05′，北緯25°09′～25°20′，東西寬為22.8 km，南北長(zhǎng)為21.8 km，狀似不規(guī)則的“凹”字形。保護(hù)區(qū)總面積為21 285 hm2，其中核心區(qū)為8305 hm2，緩沖區(qū)為8130 hm2，試驗(yàn)區(qū)為4850 hm2。

1.2 數(shù)據(jù)來源及預(yù)處理

植被一般在9—10月時(shí)生長(zhǎng)速度逐漸降低直至停止生長(zhǎng)，這段時(shí)間植被的整體變化較小，因此根據(jù)影像質(zhì)量情況，收集研究區(qū)近18年Landsat TM/ETM+/OLI影像共18幅，見表1，其中包括成像時(shí)間在9—10月的質(zhì)量良好無需處理的影像6幅，需要進(jìn)行NDVI時(shí)間變換一致性處理的目標(biāo)影像6幅及對(duì)應(yīng)年份的5—12月基準(zhǔn)影像6幅。數(shù)據(jù)來源于美國(guó)地質(zhì)勘探局(United States Geological Survey，USGS)的數(shù)據(jù)共享網(wǎng)站(http:∥glovis.usgs.gov/)。氣象數(shù)據(jù)來源于中國(guó)氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)。

由于傳感器、成像日期等的不同，需用ENVI 5.3對(duì)影像進(jìn)行裁剪、大氣校正、輻射定標(biāo)等預(yù)處理。

表1 不同時(shí)期的影像數(shù)據(jù)

1.3 植被覆蓋表示方法

文獻(xiàn)[4—6]研究表明，NDVI與植被覆蓋度(fractional vegetation cover,FVC)有很好的相關(guān)性，故本文取NDVI平均值作為反映研究區(qū)植被覆蓋變化的指標(biāo)。

NDVI計(jì)算公式為

NDVI=(NIR-Red)/(NIR+Red)

(1)

式中，NIR和Red分別代表近紅外和紅光波段的反射率。

1.4 NDVI時(shí)間變換一致性方法

針對(duì)云覆蓋對(duì)影像質(zhì)量的影響，文獻(xiàn)[7]采用填補(bǔ)法，即用成像時(shí)間相近的另一張影像的相同區(qū)域?qū)υ聘采w區(qū)進(jìn)行鑲嵌填補(bǔ)，但其精度尚未得到驗(yàn)證。為解決這一問題，本文提出了一種基于相同地物NDVI時(shí)間變化一致性的處理方法。

NDVI時(shí)間變化一致性的基本思路如圖1所示，核心思想是對(duì)于相同地物而言，NDVI的變化趨勢(shì)在短期內(nèi)是一致的。

具體實(shí)現(xiàn)方法為：首先以成像時(shí)間在9—10月附近，受云覆蓋影響無法直接用于本文研究的一幅遙感影像為目標(biāo)，尋找與其成像時(shí)間接近且影像質(zhì)量較好的另一影像為基準(zhǔn)。經(jīng)預(yù)處理后，在同時(shí)保證處理效率和處理結(jié)果質(zhì)量的前提下，經(jīng)過多次試驗(yàn)，選擇基準(zhǔn)影像中計(jì)算NDVI所需的紅光和近紅外波段進(jìn)行ISODATA非監(jiān)督分類，將研究區(qū)地物分為若干類別。再在目標(biāo)影像云覆蓋區(qū)域S附近根據(jù)不同地物分布情況設(shè)置各類別采樣點(diǎn)數(shù)量為50～200，采樣點(diǎn)中取部分用于精度驗(yàn)證，剩余的用于采集不同地物在基準(zhǔn)影像和目標(biāo)影像上的NDVI值并進(jìn)行線性擬合，得到變換規(guī)則。利用變換規(guī)則以基準(zhǔn)影像中范圍S的NDVI為基礎(chǔ)，變換得到目標(biāo)影像中范圍S的NDVI值。處理結(jié)果的精度用均方根誤差(RMSE)表示，均方根誤差s越小，表示處理結(jié)果精度越高。

1.5 一元線性回歸模型

使用一元線性回歸模型[8]進(jìn)行研究區(qū)植被覆蓋變化趨勢(shì)分析和植被覆蓋變化對(duì)氣候變化響應(yīng)的分析。該模型的原理如下

對(duì)于一元線性回歸方程y=a+bx，式中，x、y分別為自變量和因變量，a為直線的縱截距，b為直線斜率。

設(shè)有n對(duì)x與y的數(shù)值，令

W(a,b)=∑(y-a-bx)2

(2)

當(dāng)W(a,b)值最小時(shí),根據(jù)最小二乘原理計(jì)算可得a、b為

(3)

1.6 相關(guān)分析法

使用相關(guān)系數(shù)r計(jì)算自變量x與因變量y之間的相關(guān)性[9]，設(shè)有n對(duì)x與y的數(shù)值，則有

(4)

當(dāng)rxy>0時(shí)，x與y呈正相關(guān)；當(dāng)rxy<0時(shí)，x與y呈負(fù)相關(guān)。

|rxy|≤1，值越大說明x與y的相關(guān)性越高，計(jì)算得到的回歸直線越符合數(shù)據(jù)變化的趨勢(shì)：|rxy|<0.3，為弱相關(guān)；0.3≤|rxy|<0.5，為低相關(guān)；0.5≤|rxy|<0.8，為顯著相關(guān)；0.8≤|rxy|≤1，為極顯著相關(guān)。

2 結(jié)果與分析

2.1 NDVI時(shí)間變換一致性處理結(jié)果與精度評(píng)價(jià)

選擇成像時(shí)間為2003年11月14日的影像作為基準(zhǔn)影像，成像時(shí)間為2003年9月27日的影像作為目標(biāo)影像，經(jīng)預(yù)處理后計(jì)算得到兩張NDVI影像，如圖2所示。

建立一元線性、多元線性、冪、指數(shù)、對(duì)數(shù)等方程進(jìn)行回歸分析，經(jīng)過多次試驗(yàn)，發(fā)現(xiàn)各類地物使用一元線性方程擬合效果最好。對(duì)基準(zhǔn)時(shí)刻影像的紅光、近紅外波段進(jìn)行ISODATA非監(jiān)督分類，分為5種地物，經(jīng)采樣計(jì)算后得到各類地物的變換規(guī)則見表2。

利用表1中的變換規(guī)則，以基準(zhǔn)影像中范圍S的NDVI為基礎(chǔ)，變換得到目標(biāo)影像中范圍S的NDVI值，并用處理結(jié)果代替目標(biāo)影像中范圍S處的NDVI進(jìn)行鑲嵌，得到圖3(a)中NDVI時(shí)間變換一致性處理結(jié)果；用基準(zhǔn)影像中范圍S處的NDVI直接代替目標(biāo)影像中范圍S處的NDVI進(jìn)行鑲嵌，得到圖3(b)中填補(bǔ)法處理結(jié)果。

表2 NDVI時(shí)間變換一致性各類別規(guī)則與處理精度

注：x為采樣點(diǎn)在基準(zhǔn)影像上的NDVI值，y為以x為基準(zhǔn)進(jìn)行變換后的NDVI值；s1為NDVI時(shí)間變換一致性處理結(jié)果的均方根誤差；s2為填補(bǔ)法的均方根誤差。

從表2可以看出，NDVI時(shí)間變化一致性處理得到各類別變換規(guī)則均通過了P<0.001的顯著性檢驗(yàn)，說明研究區(qū)相同地物短期內(nèi)的確存在一致的變化趨勢(shì)，各個(gè)規(guī)則模型的R2均大于0.5，說明模型擬合效果較好；NDVI時(shí)間變化一致性處理后數(shù)據(jù)的均方根誤差均小于0.05，滿足使用要求，且除類別C3外，其余地物類型的均方根誤差與填補(bǔ)法相比有明顯改善。

從圖3中處理結(jié)果來看，成像日期雖相差近50 d，不過經(jīng)填補(bǔ)法處理后進(jìn)行鑲嵌的視覺效果尚可，而經(jīng)NDVI時(shí)間變換一致性處理后的NDVI與研究區(qū)目標(biāo)時(shí)刻下的NDVI鑲嵌效果更好。

以上結(jié)果表明，NDVI時(shí)間變換一致性處理后得到的研究區(qū)NDVI，不僅大幅度削弱、消除了云層覆蓋的影響，達(dá)到了本文研究的預(yù)期目的，而且其精度和效果明顯優(yōu)于填補(bǔ)法，滿足大部分研究的使用要求。

2.2 植被覆蓋變化情況

經(jīng)NDVI時(shí)間變換一致性處理后得到的研究區(qū)近18年植被覆蓋年際變化情況如圖4所示。

由圖4可知，研究區(qū)目前植被覆蓋較好，近18年中植被覆蓋呈顯著的緩慢上升趨勢(shì)，且每年的9—10月前后的平均NDVI在0.69～0.81之間上下波動(dòng)，增速為0.401(%)/a，相關(guān)系數(shù)r=0.667。此外，從圖4中還可以看出，2005、2010、2011和2015年研究區(qū)平均NDVI處于“波谷”，相應(yīng)的，研究區(qū)在這幾年都有不同程度的旱災(zāi)發(fā)生[10-11]，這從側(cè)面表明了氣候變化可能會(huì)對(duì)研究區(qū)植被生長(zhǎng)產(chǎn)生一定的影響，有必要對(duì)植被覆蓋變化與氣候的響應(yīng)關(guān)系作進(jìn)一步分析。

2.3 植被覆蓋變化與氣候因子相關(guān)分析

目前在分析植被覆蓋對(duì)氣候的響應(yīng)關(guān)系時(shí)主要趨向于使用溫度和降水兩個(gè)因子，貴州省多為陰冷濕雨天氣，因此本文除著重研究月平均溫度和降水量之外，還加入了平均相對(duì)濕度進(jìn)行分析。考慮氣候?qū)χ脖桓采w的影響可能存在滯后性[12-13]，因此，在分析時(shí)除當(dāng)月(2007年和2016年為9月，其余為10月)外還需考慮更早月份的氣候。

如表3所示，當(dāng)月平均溫度與植被覆蓋變化為顯著正相關(guān)，影響最強(qiáng)，到上1—2月時(shí)相關(guān)性降低，但仍有一定影響，上3月時(shí)只存在弱相關(guān)；植被覆蓋變化與當(dāng)月和上3月的降水量?jī)H為弱相關(guān)，與上1月降水量的相關(guān)性為顯著正相關(guān)，上前2月變?yōu)榈拖嚓P(guān)；植被覆蓋變化與當(dāng)月平均相對(duì)濕度響應(yīng)最弱，與上1月平均相對(duì)濕度顯著正相關(guān)，隨后逐漸降低，在上3月變?yōu)榈拖嚓P(guān)。

表3 植被覆蓋與各氣象因子擬合結(jié)果

注：x為各氣象因子，y為植被覆蓋。

總的來說，各氣候因子對(duì)研究區(qū)植被覆蓋變化均有顯著影響，平均溫度的影響在當(dāng)月最強(qiáng)，而降水量和平均相對(duì)濕度的影響則有一定的滯后性。

3 討論與結(jié)論

本文選用1999—2016年的Landsat TM/ETM+/OLI影像，建立了一種基于NDVI時(shí)間變換一致性的處理方法，在解決了數(shù)據(jù)來源和影像質(zhì)量的問題后，構(gòu)建出較為完整的植被指數(shù)時(shí)間序列，并基于此分析了茂蘭自然保護(hù)區(qū)的植被覆蓋變化趨勢(shì)，實(shí)現(xiàn)了對(duì)研究區(qū)植被覆蓋變化與氣象因子之間響應(yīng)關(guān)系的研究，也為小區(qū)域尺度下的植被覆蓋研究提供了技術(shù)支持。

(1)基于NDVI時(shí)間變化一致性的處理方法原理簡(jiǎn)單明了，便于操作，且很好地解決了研究區(qū)遙感影像中云覆蓋的影響，經(jīng)驗(yàn)證處理結(jié)果中各類地物的均方根誤差均小于0.05，處理精度和效果較好，滿足大部分研究的使用要求。

(2)研究發(fā)現(xiàn)，近18 a貴州茂蘭自然保護(hù)區(qū)的植被覆蓋較高且呈顯著的緩慢上升趨勢(shì)，r=0.667，每年9—10月的平均NDVI在0.69～0.81之間上下波動(dòng)，增速為0.401(%)/a，這與貴州省和保護(hù)區(qū)管理局對(duì)國(guó)家森林資源管理和生態(tài)保護(hù)政策的堅(jiān)決執(zhí)行是分不開的。

(3)經(jīng)分析，平均溫度、降水量和平均相對(duì)濕度這三種氣候因子對(duì)研究區(qū)植被覆蓋變化均存在顯著正相關(guān)，其中平均溫度的影響在當(dāng)月最強(qiáng)，相關(guān)系數(shù)r=0.60，而降水量和平均相對(duì)濕度與研究區(qū)平均NDVI的相關(guān)系數(shù)r在當(dāng)月最低，分別為-0.14和0.2，僅為弱相關(guān)；在上1月時(shí)最強(qiáng)，分別為0.67和0.70，隨后逐漸減弱，說明降水量和平均相對(duì)濕度對(duì)研究區(qū)植被覆蓋變化的影響具有滯后性。

根據(jù)本文研究，研究區(qū)植被覆蓋變化會(huì)受到氣候因子如溫度、降水和相對(duì)濕度的影響，目前在貴州省旱情頻繁、部分地區(qū)水資源短缺等種種壓力下，植被覆蓋和降水與平均相對(duì)濕度的正相關(guān)關(guān)系可能會(huì)影響森林植被的正常生長(zhǎng)；雖然目前全球變暖也會(huì)作用于平均溫度對(duì)植被覆蓋的正相關(guān)影響，但若溫度一直增長(zhǎng)直至超出某一臨界值，其造成的干旱等連鎖反應(yīng)仍會(huì)對(duì)植被生長(zhǎng)造成惡劣影響[14-16]。總而言之，如何在保持經(jīng)濟(jì)快速增長(zhǎng)的同時(shí)進(jìn)行生態(tài)環(huán)境建設(shè)、解決環(huán)境危機(jī)，這一問題亟待更多的研究和突破。