欒金凱, 劉登峰, 劉慧, 林木, 黃強(qiáng)

(1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710048; 2.中國(guó)科學(xué)院地理科學(xué)與資源研究所,陸地水循環(huán)及地表過程重點(diǎn)實(shí)驗(yàn)室,北京 100101; 3.中國(guó)科學(xué)院大學(xué),北京 100049; 4.中國(guó)水利水電科學(xué)研究院,北京 100038;5.中央財(cái)經(jīng)大學(xué) 統(tǒng)計(jì)與數(shù)學(xué)學(xué)院,北京 100081)

植被在地球生態(tài)系統(tǒng)平衡、氣候變化和水循環(huán)中起著協(xié)調(diào)的作用[1]。在流域的水文過程中,植被對(duì)蒸發(fā)、產(chǎn)流具有重要的影響。植被可在較短時(shí)間內(nèi)反映氣候變化和人類活動(dòng)的影響[2]。作為植被生長(zhǎng)及覆蓋的最佳指示因子,歸一化植被指數(shù)(Normalized Difference Vegetation Index,NDVI)是遙感數(shù)據(jù)的地表特征參量,反映植被分布特征及變化,是植被研究的基礎(chǔ)數(shù)據(jù)之一[3-6]。

在有關(guān)植被指數(shù)影響因素的研究中,研究較多的是氣溫和降水,它們與植被指數(shù)的變化存在著密切的關(guān)系。ZHAO Xia等[7]對(duì)新疆1982—2003年的NDVI分析發(fā)現(xiàn),新疆NDVI的增加與降水量和潛在蒸散發(fā)量的增加有關(guān)。FU Baihua等[8]分析了河岸植被NDVI與降水、氣溫、地表徑流、洪水及地下水位之間的關(guān)系,發(fā)現(xiàn)最高氣溫是影響NDVI值最主要的因素,其次是前28 d的降水量。趙舒怡等[9]基于像元尺度,應(yīng)用簡(jiǎn)單相關(guān)分析方法分析了2001—2013年植被指數(shù)與Palmer干旱指數(shù)(Palmer Drought Severity Index,PDSI)的相關(guān)關(guān)系。BIRTWISTLE Amy N等[10]研究發(fā)現(xiàn),可以利用植被指數(shù)數(shù)據(jù)為干旱和半干旱無(wú)雨量站觀測(cè)的偏遠(yuǎn)地區(qū)獲取降水量。GEORGANOS Stefanos等[11]運(yùn)用地理加權(quán)回歸方法,對(duì)NDVI和降水?dāng)?shù)據(jù)進(jìn)行了分析,發(fā)現(xiàn)可以使用降水?dāng)?shù)據(jù)對(duì)NDVI進(jìn)行較好的預(yù)測(cè)。WEN Zhaofei等[12]基于像元尺度,分析了三峽庫(kù)區(qū)1982—2011年的NDVI與降水、氣溫、太陽(yáng)輻射及人類活動(dòng)的相關(guān)關(guān)系。何云玲等[13]利用趨勢(shì)分析和線性相關(guān)分析方法對(duì)云南省植被的月變化趨勢(shì)、年際變化趨勢(shì)進(jìn)行了詳細(xì)分析,并利用相關(guān)分析方法分析了NDVI與主要?dú)夂蛩疅嵋蜃拥年P(guān)系。欒金凱等[14]應(yīng)用復(fù)直線回歸分析方法對(duì)榆林市近17年來(lái)的植被指數(shù)進(jìn)行分析發(fā)現(xiàn),在榆林市一半以上的區(qū)域內(nèi),人類活動(dòng)對(duì)植被生長(zhǎng)起到了促進(jìn)作用;同時(shí),還對(duì)漢江流域植被指數(shù)的影響因素進(jìn)行了研究[15-16]。

目前,對(duì)NDVI影響因素的研究主要包括3個(gè)方面:①對(duì)研究區(qū)域的植被指數(shù)取面平均值,并與降水、氣溫等影響要素進(jìn)行簡(jiǎn)單相關(guān)分析或回歸分析;②在像元尺度上進(jìn)行簡(jiǎn)單相關(guān)分析;③在研究區(qū)域內(nèi)選取氣象站點(diǎn)周圍的NDVI值或面平均值,與氣象因子和其他影響因子進(jìn)行簡(jiǎn)單相關(guān)分析和多元回歸分析。對(duì)NDVI進(jìn)行像元尺度上的多元(3個(gè)及以上影響因素)回歸分析的相關(guān)研究尚未見到。進(jìn)行NDVI像元尺度的分析,可以準(zhǔn)確地分析研究區(qū)域內(nèi)每個(gè)地點(diǎn)的影響因素對(duì)小范圍植被的影響,其空間連續(xù)性和異質(zhì)性能夠得到更好的體現(xiàn)和分析;利用多元回歸分析,可以同時(shí)考慮多個(gè)氣象因素對(duì)歸一化植被指數(shù)的影響,可以比較準(zhǔn)確地基于像元尺度預(yù)測(cè)未來(lái)植被的覆蓋狀況。因此,本研究以陜北地區(qū)的榆林市為研究區(qū)域,篩選出影響該區(qū)域歸一化植被指數(shù)的4個(gè)氣象因子,在像元尺度上對(duì)其歸一化植被指數(shù)進(jìn)行多元線性回歸分析,分析回歸系數(shù)的空間分布特征,評(píng)價(jià)回歸分析得到的歸一化植被指數(shù)值,為生態(tài)保護(hù)政策實(shí)施效果的監(jiān)測(cè)提供參考依據(jù)。

1 研究區(qū)域概況

榆林市位于陜西省最北部,東臨黃河與山西相望，西連寧夏、甘肅，北鄰內(nèi)蒙古，南接同省的延安市。經(jīng)、緯度范圍分別為107°28′E～111°15′E,36°57′N～39°35′N。地域東西長(zhǎng)385 km,南北寬263 km,總土地面積43 578 km2。地貌大體以長(zhǎng)城為界，北部為風(fēng)沙草灘區(qū)，占總面積的42%，南部為黃土丘陵溝壑區(qū)，占總面積的58%[17-18]。榆林市的位置和地形如圖1所示。

圖1 榆林市位置和地形圖

榆林市位于我國(guó)東部季風(fēng)氣候和西北干旱大陸性氣候的過渡地帶,也是毛烏素沙漠南緣與陜北黃土高原的過渡地帶,這也決定了榆林市生態(tài)環(huán)境的脆弱性。榆林市位于“三北防護(hù)林體系工程”區(qū)域內(nèi),當(dāng)前我國(guó)北方旱區(qū)森林覆蓋已經(jīng)接近氣候容量的極限,極端干旱事件正在導(dǎo)致我國(guó)北方旱區(qū)普遍出現(xiàn)森林死亡現(xiàn)象,所以考慮到未來(lái)極端氣候發(fā)生頻率增加,榆林市不宜進(jìn)一步增加森林覆蓋。在這樣的背景下,對(duì)榆林市植被的影響因素及未來(lái)發(fā)展趨勢(shì)進(jìn)行研究變得至關(guān)重要。

2 數(shù)據(jù)來(lái)源與處理

本文采用來(lái)自MODIS/Terra網(wǎng)站提供的NDVI遙感數(shù)據(jù),數(shù)據(jù)集全稱為MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid V005,簡(jiǎn)稱MOD13Q1。首先對(duì)下載的影像應(yīng)用MRT(MODIS Reprojection Tool)軟件進(jìn)行批量裁剪、投影及轉(zhuǎn)格式操作。本研究只對(duì)8月份的NDVI影像進(jìn)行分析,MOD13Q1數(shù)據(jù)8月份有2幅影像,將得到的影像進(jìn)行最大化合成(Maximum Value Composites,MVC)處理,最后每年得到一幅影像,所以2000—2017年共有18幅影像。

本文所用的8個(gè)氣象數(shù)據(jù)(日照時(shí)數(shù)、平均本站氣壓、平均最低氣溫、平均最高氣溫、平均氣溫、平均風(fēng)速、平均水汽壓、平均相對(duì)濕度)來(lái)自中國(guó)氣象數(shù)據(jù)共享服務(wù)網(wǎng)(http://data.cma.cn/),選取了榆林市內(nèi)及其周圍共17個(gè)氣象站點(diǎn)的2000—2017年共18 a的氣象數(shù)據(jù),17個(gè)站點(diǎn)(如圖1所示)分別為榆林、神木、定邊、靖邊、吳旗、橫山、綏德、延安、河曲、興縣、離石、鹽池、東勝、鄂托克旗、環(huán)縣、五寨、隰縣。利用反距離權(quán)重插值的方法來(lái)獲得每一個(gè)像元的氣象數(shù)據(jù),插值出的氣象數(shù)據(jù)的像元大小與NDVI像元大小相同且空間位置對(duì)應(yīng)。

3 研究方法

3.1 相關(guān)性分析法

分別采用Pearson相關(guān)分析法和學(xué)生t檢驗(yàn)法對(duì)NDVI與氣象數(shù)據(jù)進(jìn)行相關(guān)分析和顯著性檢驗(yàn),其計(jì)算公式如下[19]:

(1)

(2)

3.2 多元線性回歸分析法

本研究應(yīng)用多元線性回歸分析方法分析每個(gè)像元的NDVI與對(duì)應(yīng)像元的多個(gè)氣象因素之間的關(guān)系,多元線性回歸分析的表達(dá)式為[19]:

θ=A0+A1x1+A2x2+…+Anxn+ζ。

(3)

式中:θ為NDVI值;A0、A1、…、An為回歸系數(shù);x1、x2、…、xn為氣象因素,本文中n=4;ξ為殘差,即實(shí)際NDVI值與線性回歸分析結(jié)果的差值。

4 結(jié)果分析

4.1 氣象因子的選取

本文選取了日照時(shí)數(shù)、平均本站氣壓、平均最低氣溫、平均最高氣溫、平均氣溫、平均風(fēng)速、平均水汽壓和平均相對(duì)濕度共8個(gè)氣象因子。由于氣象因子存在空間差異性,選取了榆林市內(nèi)及其周圍共17個(gè)氣象站點(diǎn)資料,利用反距離權(quán)重插值的方法來(lái)獲得每一個(gè)像元的氣象數(shù)據(jù),插值出的氣象數(shù)據(jù)的像元大小與NDVI像元相同?？紤]到氣象因子對(duì)NDVI的影響存在滯后性,將榆林市2000—2017年8月份的NDVI分別與6月、7月、8月的各氣象因子及6—7月、7—8月、6—8月各氣象因子的平均值進(jìn)行相關(guān)分析,相關(guān)系數(shù)見表1。由表1可知,與榆林市8月份NDVI相關(guān)性最大的各氣象因子分別為7月日照時(shí)數(shù)、6—7月平均氣壓、7月平均最低氣溫、7月平均最高氣溫、7月平均氣溫、6月平均風(fēng)速、7月平均水汽壓和7月相對(duì)濕度。

表1 榆林市8月份NDVI與各時(shí)段氣象因子間的相關(guān)系數(shù)

對(duì)這8個(gè)與NDVI相關(guān)性最大的氣象因子進(jìn)行t分布檢驗(yàn),結(jié)果見表2。由表2發(fā)現(xiàn)，6—7月平均氣壓、7月平均最高氣溫、7月平均氣溫、7月平均相對(duì)濕度通過了95%的置信度檢驗(yàn)。因此,本文選擇6—7月平均氣壓、7月平均最高氣溫、7月平均氣溫和7月平均相對(duì)濕度作為NDVI的主要影響因子。

表2 部分氣象因子的t分布檢驗(yàn)值

4.2 多元線性回歸分析

選擇上述17個(gè)氣象站點(diǎn)2000—2017年6—7月平均氣壓、7月平均最高氣溫、7月平均氣溫、7月平均相對(duì)濕度,運(yùn)用反距離權(quán)重法進(jìn)行空間插值。插值后的柵格數(shù)據(jù)的像元大小和MODIS NDVI的像元大小相同,每幅影像為1 328×1 559的矩陣,去除周圍不在榆林市境內(nèi)的區(qū)域,共得到804 015個(gè)值,即榆林市共有804 015個(gè)像元。將2000—2014年共15 a的數(shù)據(jù)作為模擬數(shù)據(jù)用于建立模型,2015—2017年共3 a的數(shù)據(jù)用于驗(yàn)證模型。經(jīng)過回歸分析可得各像元的回歸方程，其表達(dá)式的一般形式為：

θ=A0+A1x1+A2x2+A3x3+A4x4。

(4)

式中:θ為NDVI的回歸分析值;x1為7月平均相對(duì)濕度;x2為6—7月平均氣壓;x3為7月平均氣溫;x4為7月平均最高氣溫;A0為回歸常數(shù);A1、A2、A3、A4分別代表7月平均相對(duì)濕度、6—7月平均氣壓、7月平均氣溫和7月平均最高氣溫的回歸系數(shù)。

各回歸參數(shù)(包括回歸常數(shù)和各回歸系數(shù))的空間分布及頻率分布分別如圖2和圖3所示。由圖2可知：回歸常數(shù)在大部分區(qū)域?yàn)檎?在榆陽(yáng)區(qū)和神木市部分區(qū)域?yàn)樨?fù)值,西部的定邊縣和東南部的部分區(qū)域?yàn)?;7月平均相對(duì)濕度的回歸系數(shù)在榆林大部分區(qū)域?yàn)檎?說(shuō)明在大部分區(qū)域,制約榆林植被生長(zhǎng)的主要因素是水。從空間分布來(lái)看,6—7月平均氣壓與7月平均相對(duì)濕度對(duì)植被的影響基本呈現(xiàn)相反的趨勢(shì),7月平均氣溫與7月平均相對(duì)濕度的回歸系數(shù)也呈現(xiàn)出相反的趨勢(shì),而7月平均最高氣溫的回歸系數(shù)與7月平均氣溫的回歸系數(shù)呈現(xiàn)出相反的趨勢(shì);不同的影響因子對(duì)同一地區(qū)的響應(yīng)不同,不同區(qū)域?qū)ν挥绊懸蜃拥捻憫?yīng)也不同。

圖2 各回歸參數(shù)的空間分布圖

圖3 各回歸參數(shù)的頻率分布圖

表3為各回歸參數(shù)(804 015個(gè)值)的基本統(tǒng)計(jì)特征值,包括平均值和標(biāo)準(zhǔn)差。

由表3可知,各影響因素的系數(shù)分布分散,說(shuō)明不同地點(diǎn)上影響因素的作用有較大差異。

表3 各回歸參數(shù)的基本特征統(tǒng)計(jì)值

4.3 NDVI模擬值的評(píng)價(jià)

榆林市2015—2017年各年NDVI模擬值與觀測(cè)值的空間分布如圖4所示。分析圖4發(fā)現(xiàn)：榆林市各年NDVI的模擬值與觀測(cè)值的空間分布特征基本一致,毛烏素沙漠邊緣NDVI值較低,榆林東部和東南部NDVI值較高;整體的模擬效果較好,其中2016年的模擬效果最好,2017年的模擬效果較差。

圖4 榆林市2015—2017年模擬NDVI與實(shí)際NDVI的空間分布圖

為了定量評(píng)價(jià)逐個(gè)像元的模擬效果,將實(shí)際NDVI值減去模擬NDVI值的結(jié)果記為殘差。殘差可視為是人類活動(dòng)對(duì)歸一化植被指數(shù)的影響。如果殘差為負(fù)值,說(shuō)明人類活動(dòng)對(duì)該區(qū)域的植被起到了抑制作用,可能出現(xiàn)了濫砍濫伐、植被破壞或者植被死亡等現(xiàn)象。如果殘差為正值,說(shuō)明人類活動(dòng)對(duì)植被的生長(zhǎng)起到了積極的促進(jìn)作用。殘差的絕對(duì)值與實(shí)際NDVI結(jié)果的比值記為相對(duì)誤差,是評(píng)價(jià)模擬效果的關(guān)鍵指標(biāo)。相對(duì)誤差越大,說(shuō)明人類活動(dòng)對(duì)植被的影響越大;反之越小。

榆林市2015—2017年各年NDVI模擬值的殘差和相對(duì)誤差的空間分布如圖5所示。從圖5可以看出：2015年殘差為負(fù)值的區(qū)域主要位于靖邊縣、橫山區(qū)、子洲縣、神木市的北部及府谷縣的南部,說(shuō)明在這些區(qū)域存在著植被破壞或死亡的現(xiàn)象;其余區(qū)域殘差值均為正值,說(shuō)明在這些區(qū)域內(nèi),人類活動(dòng)對(duì)植被生長(zhǎng)起到了積極的促進(jìn)作用,這些地區(qū)實(shí)施退耕還林、荒山荒地造林、封山育林等措施的效果較好。2016年殘差為負(fù)值的區(qū)域主要分布在神木市的北部和府谷縣的南部。2017年殘差為負(fù)值的區(qū)域主要位于靖邊縣的西部,佳縣、米脂縣和子洲縣的部分區(qū)域。

圖5 榆林市2015—2017年各年NDVI模擬值的殘差和相對(duì)誤差分布圖

榆林市2015—2017年各年NDVI模擬值的殘差及相對(duì)誤差的統(tǒng)計(jì)情況分別見表4和表5。由表4和表5可知：對(duì)于榆林市NDVI的模擬值,2015年的殘差在[-0.1,0.1)上的像元比例為66.50%,相對(duì)誤差在25%以內(nèi)的像元比例為68.68%;2016年的殘差在[-0.1,0.1)上的像元比例為77.62%,相對(duì)誤差在25%以內(nèi)的像元比例達(dá)到了88.93%;2017年的殘差在[-0.1,0.1)上的像元比例為72.82%,相對(duì)誤差在25%以內(nèi)的像元比例為89%。由此可以看出：2016年的模擬效果最好,該年人類活動(dòng)對(duì)植被的影響最小;2015年的模擬結(jié)果最差,該年人類活動(dòng)對(duì)植被的影響最為劇烈。

表4 榆林市2015—2017年NDVI模擬值的殘差統(tǒng)計(jì)

表5 榆林市2015—2017年NDVI模擬值的相對(duì)誤差統(tǒng)計(jì)

由表4和表5還可以得知,榆林市2015—2017年NDVI模擬值的殘差為正值的像元比例分別為60.14%、79.31%和85.73%,說(shuō)明最近幾年人類活動(dòng)對(duì)植被生長(zhǎng)和更新起到促進(jìn)作用的區(qū)域占到了很大部分。查閱榆林市林業(yè)局網(wǎng)站發(fā)現(xiàn)[20],榆林市2015年總造林面積45萬(wàn)畝,2016年和2017年分別造林73.3萬(wàn)畝和70萬(wàn)畝,這也導(dǎo)致2016和2017年殘差為正值的區(qū)域面積大于2015年的。另外,2016年橫山區(qū)造林面積15 300萬(wàn)畝,所以該區(qū)NDVI模擬值的殘差由2015年的負(fù)值變?yōu)榱?016年的正值。2016年春季,府谷縣存在嚴(yán)重的偷牧亂牧現(xiàn)象,所以2016年該區(qū)NDVI模擬值的殘差值小于0。2017年7月,綏德縣和子洲縣受到暴雨和洪災(zāi)的影響,當(dāng)?shù)亓謽I(yè)受損嚴(yán)重,導(dǎo)致2017年這兩個(gè)縣NDVI模擬值的殘差值呈現(xiàn)負(fù)值。由此可見,應(yīng)用本文所采用的方法,可以作為監(jiān)測(cè)封山育林、退耕還林、退牧還草等政策的落實(shí)和實(shí)施效果的參考;在現(xiàn)階段人類活動(dòng)變化幅度不大的情況下,還可以根據(jù)未來(lái)的氣象條件預(yù)測(cè)未來(lái)植被理論上的生長(zhǎng)狀況。

對(duì)榆林市2015—2017年NDVI的實(shí)測(cè)值與模擬值進(jìn)行逐像元相關(guān)系數(shù)分析,并作實(shí)際值與模擬值的散點(diǎn)圖,橫坐標(biāo)為榆林市804 015個(gè)像元的NDVI實(shí)測(cè)值,縱坐標(biāo)為其模擬值,如圖6所示。由圖6可知：2017年的模擬效果最好,相關(guān)系數(shù)達(dá)到了0.851 0;2015年的相關(guān)系數(shù)為0.616 6;2016年的相關(guān)系數(shù)為0.803 3。

圖6 榆林市2015—2017年NDVI實(shí)測(cè)值與模擬值的散點(diǎn)圖

5 結(jié)論

本文利用MODIS/Terra NDVI時(shí)間序列數(shù)據(jù),基于像元尺度對(duì)陜西省榆林市2000—2017年生長(zhǎng)季(8月份)的歸一化植被指數(shù)(NDVI)應(yīng)用多元線性回歸方法進(jìn)行影響因素分析研究,主要得到以下結(jié)論:

1)與榆林市8月份NDVI相關(guān)程度最大的因子為:6—7月平均氣壓、7月平均最高氣溫、7月平均氣溫和7月平均相對(duì)濕度,可視為榆林市8月份NDVI的主要影響因子。

2)不同的影響因子對(duì)同一地區(qū)的響應(yīng)不同,不同區(qū)域?qū)ν挥绊懸蜃拥捻憫?yīng)也不同。7月平均相對(duì)濕度與6—7月平均氣壓對(duì)植被的影響基本呈現(xiàn)相反的趨勢(shì)。7月平均氣溫與7月平均相對(duì)濕度的回歸系數(shù)呈現(xiàn)出相反的趨勢(shì),7月平均最高氣溫的回歸系數(shù)與7月平均氣溫的回歸系數(shù)也呈現(xiàn)出相反的趨勢(shì)。

3)對(duì)比NDVI的模擬值和實(shí)測(cè)值發(fā)現(xiàn):整體的模擬效果較好,2015年NDVI模擬值的殘差為負(fù)值的區(qū)域主要位于靖邊縣、橫山區(qū)、子洲縣、神木市的北部及府谷縣的南部,2016年NDVI模擬值的殘差為負(fù)值的區(qū)域主要分布在神木市的北部和府谷縣的南部,2017年NDVI模擬值的殘差值為負(fù)值的區(qū)域主要位于靖邊縣的西部,佳縣、米脂縣和子洲縣的部分區(qū)域。殘差為正值的區(qū)域,說(shuō)明人類活動(dòng)對(duì)該區(qū)域內(nèi)植被的生長(zhǎng)起到了促進(jìn)作用,這些地區(qū)的封山育林、退耕還林、退牧還草等措施的實(shí)施效果較好。殘差值為負(fù)值的區(qū)域,說(shuō)明人類活動(dòng)對(duì)該區(qū)域內(nèi)植被的生長(zhǎng)起到了抑制作用,可能存在亂砍亂伐、植被死亡等現(xiàn)象。應(yīng)用這種方法,在現(xiàn)階段人類活動(dòng)變化幅度不大的情況下,可以根據(jù)未來(lái)的氣象條件預(yù)測(cè)未來(lái)植被理論上的生長(zhǎng)狀況,可以為區(qū)域生態(tài)環(huán)境保護(hù)效果的評(píng)價(jià)提供參考。