段晨陽(yáng)，王曉艷

（1.西安科技大學(xué) 測(cè)繪科學(xué)與技術(shù)學(xué)院，陜西 西安 710054）

植被作為生態(tài)系統(tǒng)的平衡調(diào)節(jié)器和生態(tài)變化的指示器，能夠?yàn)樯鷳B(tài)變化、環(huán)境演化研究提供科學(xué)依據(jù)。近年來(lái)基于時(shí)間序列的衛(wèi)星遙感植被覆蓋度變化分析逐漸增多，為植被監(jiān)測(cè)提供了數(shù)據(jù)支撐，歸一化植被指數(shù)（NDVI）對(duì)植物檢測(cè)靈敏度較高，能表示地表綠度和植被覆蓋特征，可廣泛應(yīng)用于大尺度地表植被活動(dòng)監(jiān)測(cè)與評(píng)估[1-3]。國(guó)內(nèi)外學(xué)者對(duì)植被覆蓋變化采用不同方法進(jìn)行研究，多是關(guān)于植被變化及產(chǎn)生因素、植被動(dòng)態(tài)變化監(jiān)測(cè)、物種分布等[4-13]?？臻g自相關(guān)是地理空間變量分布與相鄰位置間相關(guān)性的研究，通過(guò)監(jiān)測(cè)分析區(qū)域的變異是否依賴(lài)于臨近位置上的變異來(lái)判別是否存在空間相關(guān)性[8]，目前主要應(yīng)用于經(jīng)濟(jì)[9]、人口[10]、土地利用[11]等方面，也有一些學(xué)者利用空間自相關(guān)分析植被覆蓋在空間上的相關(guān)性[12-13]。MODIS NDVI在區(qū)域小比例尺范圍內(nèi)植被覆蓋動(dòng)態(tài)變化分析上具有較好的應(yīng)用。本文以山東省為研究區(qū)，進(jìn)行植被覆蓋的動(dòng)態(tài)變化以及空間自相關(guān)分析，首先對(duì)山東省植被覆蓋進(jìn)行等級(jí)劃分，將時(shí)間序列變化的MODIS影像與空間自相關(guān)分析相結(jié)合，利用一元線性回歸方法，擬合每一個(gè)像元的變化趨勢(shì)，并結(jié)合植被覆蓋重心遷移路徑分析，以揭示植被覆蓋在時(shí)間和空間上的變化規(guī)律，更好地表示植被覆蓋的動(dòng)態(tài)演變規(guī)律以及植被覆蓋的空間聚集程度，為山東省生態(tài)環(huán)境保護(hù)工作提供科學(xué)依據(jù)[14]。

1 研究區(qū)概況

山東省地處中國(guó)東部沿海地區(qū)，位于34°22′～38°23′N(xiāo)、114°19′～122°43′E之間，地域遼闊，受溫帶季風(fēng)氣候影響。圖1為山東省DEM圖。

圖1 研究區(qū)DEM圖

2 數(shù)據(jù)源與研究方法

2.1 數(shù)據(jù)來(lái)源及預(yù)處理

本研究采用美國(guó)地球資源觀測(cè)系統(tǒng)數(shù)據(jù)中心提供的植被指數(shù)產(chǎn)品（MOD13Q1），空間分辨率為250 m、時(shí)間分辨率為16 d。選取2008—2020年植被旺盛季節(jié)[15]即每年6—9月份影像，合計(jì)85景。利用MRT軟件對(duì)圖像進(jìn)行NDVI提取，投影到Albers投影坐標(biāo)系統(tǒng)下。為了去除云、大氣以及太陽(yáng)高度角等的影響，采用國(guó)際通用的最大值合成法（MVC）對(duì)數(shù)據(jù)進(jìn)行時(shí)間尺度擴(kuò)展計(jì)算[16]，利用山東省矢量邊界對(duì)圖像進(jìn)行裁剪。

利用ArcGIS地圖代數(shù)工具對(duì)影像進(jìn)行異常剔除以及歸一化處理[0，1]，利用焦點(diǎn)統(tǒng)計(jì)算法，選擇用于計(jì)算統(tǒng)計(jì)信息的每個(gè)單元格周?chē)鷧^(qū)域的形狀為矩形，統(tǒng)計(jì)類(lèi)型為平均值。隨機(jī)生成點(diǎn)文件，將焦點(diǎn)統(tǒng)計(jì)后的影像進(jìn)行重分類(lèi)，通過(guò)提取值到點(diǎn)工具，將柵格數(shù)據(jù)的NDVI值批量提取到對(duì)應(yīng)點(diǎn)的屬性表中，以便后期空間相關(guān)性分析。

采用Anselin等開(kāi)發(fā)的Geoda軟件進(jìn)行空間自相關(guān)分析，得到山東省2008—2018年NDVI全局Moran’s I指數(shù)，局部自相關(guān)分析得到空間自相關(guān)的LISA（LISA）聚集圖。

2.2 研究方法

采用MVC法進(jìn)行NDVI時(shí)間序列計(jì)算。

式中，NDVIi為第i年NDVI最大值；NDVIt為各像元植被16 d合成NDVI最大值；d為第i年16 d合成NDVI影像總數(shù)目。

采用一元線性回歸方法對(duì)每個(gè)像元NDVI值進(jìn)行時(shí)間序列變化分析，利用最小二乘法擬合像元不同年份的斜率來(lái)反映區(qū)域植被覆蓋變化的時(shí)空特征[15]。

式中，θslop為NDVI值變化斜率；n為時(shí)間序列值；i=1，2，3…，1為年號(hào)；Ni為第i年NDVI值。當(dāng)斜率為正值時(shí)，NDVI增加，表示植被覆蓋度增加，反之減少。

重心可以表示地理要素的時(shí)空分布特征，移動(dòng)方向指向高密度區(qū)域，移動(dòng)距離反映要素變化的空間差異程度。將研究區(qū)域劃分為若干子區(qū)域，每個(gè)子區(qū)域的重心坐標(biāo)為（Xi，Yi），Mi為該子區(qū)域某種意義下的重量，(Xˉ，Yˉ)為子區(qū)域重心坐標(biāo)在重量影響下的加權(quán)平均值。本文以植被NDVI作為研究區(qū)的重量，探索植被覆蓋在時(shí)空變化下的規(guī)律。

應(yīng)用空間自相關(guān)方法對(duì)山東省植被覆蓋空間分布進(jìn)行分析研究。通常采用全局和局部?jī)煞N指標(biāo)來(lái)衡量，全局指標(biāo)是用于探測(cè)某現(xiàn)象在整個(gè)研究區(qū)域的空間分布模式，分析其是否有聚集特性存在；局部指標(biāo)是反映整個(gè)大區(qū)域中，局部小區(qū)域單元上的屬性值與相鄰局部小區(qū)域單元上同一屬性值的相關(guān)程度[17]。本文運(yùn)用應(yīng)用較為廣泛的Moran’s I指數(shù)進(jìn)行研究。

首先確定空間權(quán)重矩陣，通常用一個(gè)二元對(duì)稱(chēng)空間權(quán)重矩陣W來(lái)表示n個(gè)區(qū)域位置的鄰近關(guān)系[18]，如公式（4）：

式中，n為空間單元個(gè)數(shù)；wij為區(qū)域i和的鄰接關(guān)系。全局Moran’s I指數(shù)為（5）：

式中，xi、xj為空間位置i和j的觀察值，全局Moran’s I的取值范圍為[-1，1]，樣本值越趨近于1空間差異性越小，反之空間差異性越大。正值表示樣本值趨于聚集，負(fù)值表示樣本值在空間上的聚集性逐漸減小或者消失，0表示樣本值在空間上相互獨(dú)立，沒(méi)有相關(guān)性。

對(duì)于全局Moran’s I指數(shù)，可以用標(biāo)準(zhǔn)化統(tǒng)計(jì)量Z來(lái)檢驗(yàn)，Moran’s I指數(shù)檢驗(yàn)的標(biāo)準(zhǔn)化統(tǒng)計(jì)量為：

式中，以正態(tài)分布95%置信區(qū)間雙側(cè)檢驗(yàn)閾值1.96為界限[6]，E(Ii)和Z(Ii)分別為理論期望和理論方差。當(dāng)Z(Ii)＞1.96時(shí)，空間關(guān)系為正相關(guān)；當(dāng)Z(Ii)＜-1.96時(shí)，空間關(guān)系為負(fù)相關(guān)；當(dāng)Z(Ii)介于-1.96～1.96之間時(shí)，觀測(cè)值呈隨機(jī)分布。

全局空間自相關(guān)假定空間同質(zhì)，但區(qū)域要素之間普遍存在異質(zhì)性，相同變量在不同空間尺度上的空間自相關(guān)程度還是有差別的，區(qū)域空間自相關(guān)可以更好地解釋植被空間異質(zhì)特性。

局部Moran’s I公式為：

3 結(jié)果與分析

3.1 NDVI的空間分布

NDVI是植物生長(zhǎng)狀態(tài)以及空間分布密度的最佳指示因子，與植被分布密度呈線性關(guān)系。將NDVI指數(shù)劃分為5個(gè)級(jí)別[2]：小于0.1的區(qū)域?yàn)榉侵脖粎^(qū)；0.1～0.3為植被覆蓋匱乏地區(qū)；0.3～0.6為植被覆蓋較匱乏地區(qū)；0.6～0.8為植被覆蓋良好地區(qū)；大于0.8為植被覆蓋優(yōu)秀地區(qū)。山東省植被類(lèi)型主要以農(nóng)業(yè)為主，主要分布在魯西北平原。落葉闊葉林、針葉林、灌叢、草地等主要分布在魯中南山地區(qū)陵、膠東丘陵地區(qū)。圖2為山東省2008—2020年的NDVI空間分布圖，近13 a來(lái)NDVI較高的地區(qū)在魯西南及魯西北等農(nóng)業(yè)區(qū)，NDVI均在0.6以上。中部山地地區(qū)以及半島丘陵地區(qū)NDVI指數(shù)處適中狀態(tài)。NDVI指數(shù)較低區(qū)域?yàn)轸斘髂系貐^(qū)的南四湖及東平湖、北部的黃河三角洲以及萊州灣地區(qū)以及城鎮(zhèn)建設(shè)地區(qū)。由此可見(jiàn)，植被覆蓋的空間分布受到土地利用類(lèi)型的影響。

圖2 NDVI空間分布

3.2 NDVI時(shí)間變化

時(shí)間序列的遙感影像圖能夠反映植被覆蓋的變化情況，采用一元線性回歸分析對(duì)每一個(gè)柵格NDVI的變化趨勢(shì)進(jìn)行模擬，斜率代表了該柵格的NDVI在2008—2018年間的變化趨勢(shì)。若斜率為正值，表示植被呈增長(zhǎng)趨勢(shì)；反之呈減少趨勢(shì)。

利用MATLAB處理計(jì)算，斜率最大值為0.102，最小值為-0.104。依據(jù)斜率的大小情況，劃分為5個(gè)等級(jí)[19]，-0.104～-0.02為顯著減少；-0.02～-0.005為輕微減少；-0.005～0.005為保持不變；0.005～0.02為輕微增加，0.02～0.102為顯著增加。由表1可知植被輕微增加的像元數(shù)最多占62.5%，相對(duì)來(lái)說(shuō)輕微減少的像元所占比重要最小，僅為2.1%，得出近年來(lái)山東省植被呈現(xiàn)出穩(wěn)定和增長(zhǎng)趨勢(shì)。由圖3可見(jiàn)，植被變化趨勢(shì)有明顯的空間分布差異，大部分地區(qū)呈穩(wěn)定增長(zhǎng)趨勢(shì)。其中植被顯著增長(zhǎng)地區(qū)主要集中在聊城、德州平原地區(qū)以及濰坊等農(nóng)業(yè)區(qū)，由于實(shí)施鄉(xiāng)村振興戰(zhàn)略，農(nóng)業(yè)得到很大發(fā)展，糧食產(chǎn)量創(chuàng)新高，造林綠化使得植被覆蓋增加。而北部的黃河三角洲及萊州灣、湖泊濕地及沿海部分地區(qū)有輕微導(dǎo)致植被覆蓋減少，而湖泊濕地周?chē)脖粯O不穩(wěn)定，受環(huán)境人為因素干擾大，出現(xiàn)減少趨勢(shì)。

表1 各變化等級(jí)統(tǒng)計(jì)

圖3 植被覆蓋變化趨勢(shì)

3.3 重心演化分析

為揭示山東省植被覆蓋在時(shí)間以及空間上的動(dòng)態(tài)變化，計(jì)算2008—2020年時(shí)間序列變化的植被空間重心坐標(biāo)。如圖4所示，取間隔2 a的NDVI重心坐標(biāo)生成散點(diǎn)圖，箭頭的指向即為重心遷移變化方向。由圖可知，2008—2020年間植被覆蓋的重心發(fā)展方向經(jīng)歷幾個(gè)過(guò)程，即東北方向、西南方向、西北方向，再到西南方向。結(jié)合圖3植被變化趨勢(shì)作對(duì)比分析，2008—2012年植被覆蓋重心坐標(biāo)往東北偏移，濰坊地區(qū)大力發(fā)展農(nóng)業(yè)，打造新型蔬菜大棚，著力于蔬菜出口貿(mào)易，使得植被覆蓋持續(xù)增加。2012—2014年植被覆蓋重心坐標(biāo)往西南偏移，菏澤、濟(jì)寧地區(qū)位于平原地區(qū)，農(nóng)業(yè)大力發(fā)展，使得植被覆蓋增加。2014—2018年植被覆蓋重心坐標(biāo)往西北偏移，聊城、德州地區(qū)農(nóng)業(yè)生產(chǎn)得到迅速發(fā)展，糧食產(chǎn)量大大提高，黃河三角洲地區(qū)實(shí)施生態(tài)保護(hù)政策，植被覆蓋度大大提高。2018—2020年往西南方向發(fā)展。植被覆蓋重心的遷移變化受到很多因素的影響，可能還受溫度降水等綜合影響。

圖4 NDVI重心遷移路徑

3.4 空間自相關(guān)分析

利用Geoda軟件進(jìn)行全局自相關(guān)分析，用Weight工具創(chuàng)建研究區(qū)空間權(quán)重矩陣文件，依照Threshold Distance規(guī)則創(chuàng)建空間權(quán)重矩陣，得到全局Moran’s I為0.898，通過(guò)了P＜0.05的顯著性檢驗(yàn)。Moran’s I系數(shù)反映了空間鄰近點(diǎn)NDVI的相似程度[20]，隨機(jī)點(diǎn)集中分布在第三象限，表明山東省表現(xiàn)為顯著的正空間自相關(guān)，植被覆蓋呈現(xiàn)出聚集狀態(tài)。山東省大部分地區(qū)植被類(lèi)型為1年一熟或者2年三熟的農(nóng)作物以及草地林地過(guò)渡區(qū)域，植被覆蓋季節(jié)性變動(dòng)態(tài)化小，植被覆蓋穩(wěn)定性強(qiáng)[21]，這些因素對(duì)于植被覆蓋在空間上的顯著相關(guān)性產(chǎn)生一定影響。

全局空間自相關(guān)分析是對(duì)整個(gè)研究區(qū)內(nèi)的地理空間要素整體綜合分析，來(lái)探求整體表現(xiàn)出來(lái)的規(guī)律，是一個(gè)整體的反映。但是對(duì)于內(nèi)部來(lái)說(shuō)，各個(gè)區(qū)域的局部空間自相關(guān)性并不是完全一樣的，經(jīng)常表現(xiàn)出不同程度的空間異質(zhì)性，因此對(duì)研究區(qū)進(jìn)行局部空間自相關(guān)分析，能夠發(fā)現(xiàn)全局分析沒(méi)有發(fā)現(xiàn)的規(guī)律，從而更好的揭示空間聚集的一般規(guī)律。由Moran’s I的數(shù)值看出NDVI在空間上表現(xiàn)出較強(qiáng)的正相關(guān)，但是沒(méi)有清晰地表明相關(guān)性在空間中的分布情況，而LISA圖能夠很好地彌補(bǔ)這一不足。

圖5為通過(guò)P＜0.05顯著性檢驗(yàn)的NDVI LISA聚集圖，除北部的黃河三角洲和萊州灣、南四湖、沿海以及零星地區(qū)NDVI表現(xiàn)為低-低自相關(guān)以外，山東省大部分地區(qū)呈高-高自相關(guān)趨勢(shì)。山東省為農(nóng)業(yè)大省，植被類(lèi)型大多以農(nóng)用地為主[22]，平原面積占全省面積的55.56%，主要分布在魯西北地區(qū)和魯西南局部地區(qū)，所以這些地區(qū)植被覆蓋總體較好，高-高聚集現(xiàn)象明顯；魯中南為山地區(qū)陵地區(qū)，主要為森林、灌叢和草地等，植被覆蓋總體良好，聚集現(xiàn)象不那么強(qiáng)；魯東沿海地區(qū)，旅游業(yè)發(fā)展迅速，開(kāi)發(fā)建設(shè)導(dǎo)致植被覆蓋整體較差，低-低聚集明顯；黃河三角洲以及萊州灣地區(qū)一直為NDVI低值聚集區(qū)，距離黃河河道越遠(yuǎn)，渤海海岸線越近，NDVI值越小。近海地區(qū)常年受到黃河水沙及人類(lèi)活動(dòng)的影響，土壤含鹽量較高，生態(tài)脆弱，植物種類(lèi)匱乏，覆蓋率低，低-低聚集現(xiàn)象明顯。

圖5 山東省NDVILISA圖

4 結(jié)語(yǔ)

本文以2008—2020年的MODIS影像為研究基礎(chǔ)，對(duì)山東省植被覆蓋的時(shí)間序列的變化情況，以及植被覆蓋在空間上的關(guān)聯(lián)情況做出研究，得到了以下結(jié)論：

1）近年來(lái)，山東植被覆蓋整體表現(xiàn)為穩(wěn)定上升的趨勢(shì)，75%的區(qū)域保持增長(zhǎng)趨勢(shì)，呈增加趨勢(shì)的區(qū)域遠(yuǎn)遠(yuǎn)大于呈減退趨勢(shì)的區(qū)域。

2）植被變化趨勢(shì)在空間上具有一定差異性，西部地區(qū)為增長(zhǎng)型，東部地區(qū)為穩(wěn)定性。植被覆蓋重心往魯西北、魯西南農(nóng)業(yè)區(qū)偏移。

3）全局Moran’s I均值為0.898，顯著性檢驗(yàn)P＜0.05，置信度較高，該結(jié)果表明山東省植被覆蓋有較強(qiáng)的空間聚集性。局部空間自相關(guān)分析來(lái)看，山東省大部分地區(qū)均表現(xiàn)為高-高自相關(guān)，植被覆蓋度較好；魯東沿海地區(qū)、黃河三角洲和萊州灣地區(qū)以及各地零星地區(qū)表現(xiàn)為低-低自相關(guān)，植被覆蓋度較差。

基于時(shí)間序列的遙感影像圖結(jié)合空間自相關(guān)方法，對(duì)山東省植被覆蓋進(jìn)行研究，能夠客觀地揭示植被覆蓋動(dòng)態(tài)變化規(guī)律以及空間聚集情況。統(tǒng)計(jì)學(xué)方法的數(shù)據(jù)指標(biāo)能夠反映一定的客觀規(guī)律，但是植被覆蓋的變化與諸多因素有關(guān)。在以后的研究工作中，應(yīng)綜合多種因素更深入的進(jìn)行探究。