基于CA-Markov的渭河流域NDVI時(shí)空變化模擬及預(yù)測

王麗霞, 張珈瑋, 孟妮娜, 隋立春,4, 張雙成, 劉 招

(1.長安大學(xué) 地質(zhì)工程與測繪學(xué)院, 西安 710054; 2.長安大學(xué) 地球科學(xué)與資源學(xué)院, 西安 710054;3.長安大學(xué) 環(huán)境科學(xué)與工程學(xué)院, 西安 710054; 4.地理國情監(jiān)測國家測繪地理信息局 工程技術(shù)研究中心, 西安 710054)

植被是全球生態(tài)系統(tǒng)的重要組成部分，是連接土壤、大氣和水分的中間環(huán)節(jié)[1-2]。植被覆蓋在很大程度上代表了整體的生態(tài)環(huán)境，地表植被覆蓋狀況如何，直接關(guān)系到區(qū)域生態(tài)環(huán)境的穩(wěn)定性和安全性[3]。遙感數(shù)據(jù)提供了關(guān)于植被類型/植被覆蓋動態(tài)變化過程的重要信息，目前一般通過遙感影像數(shù)據(jù)提取各類植被指數(shù)，在不同尺度下對植被覆蓋進(jìn)行長時(shí)間序列的動態(tài)監(jiān)測。其中NDVI是一種應(yīng)用最為廣泛的植被指數(shù)，在生態(tài)環(huán)境監(jiān)測、旱情監(jiān)測、植被生長力評估和土地利用評價(jià)等方面已有較多應(yīng)用。近年來，國內(nèi)外學(xué)者主要針對植被覆蓋現(xiàn)狀的時(shí)空分布進(jìn)行了較為深入的研究[4-9]。CA-Markov模型綜合了Markov模型能夠進(jìn)行長期預(yù)測的優(yōu)點(diǎn)與元胞自動機(jī)(CA)的模擬復(fù)雜系統(tǒng)空間變化的能力，具有較高的預(yù)測精度，在分析和模擬土地利用變化方面取得了良好的效果[10-15]。

渭河流域地處我國中心區(qū)域，地理位置介于104°00′—110°20′E，33°50′—37°18′N，橫跨陜西、甘肅、寧夏3省(自治區(qū))，面積13.5萬km2，是黃河最大的子流域。渭河流域地勢南北高，中間低，由西向東傾斜，地處半干旱、半濕潤地區(qū)的過渡地帶，屬溫帶大陸性季風(fēng)氣候，降水時(shí)空分布不均，南多北少，主要集中在夏季，年降水量為500～800 mm，年平均氣溫介于7.8～13.5℃[16]。流域受地理位置、地貌及氣候等因素的影響，植被類型多樣。關(guān)中平原、涇河、北洛河河谷地區(qū)以及甘肅隴東地區(qū)主要的植被類型為栽培植被，包括農(nóng)作物、耐寒經(jīng)濟(jì)作物以及落葉果樹等。秦嶺、六盤山、子午嶺地區(qū)主要分布著溫帶落葉闊葉林、溫帶針葉林、溫帶灌叢等自然植被，流域西北則以溫帶草原為主。該地區(qū)處于南北方的地形、氣候過渡區(qū)，是連通西北、西南的咽喉要地，水土流失嚴(yán)重，是我國生態(tài)環(huán)境脆弱的典型區(qū)域，綜合分析該流域的植被覆蓋變化狀況，對區(qū)域生態(tài)環(huán)境保護(hù)具有重要意義。

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

1.1 遙感數(shù)據(jù)來源及處理

本文選取的NDVI數(shù)據(jù)為MOD13Q1數(shù)據(jù)，數(shù)據(jù)來源于NASA網(wǎng)站(https:∥ladsweb.modaps.eosdis.nasa.gov/)。該數(shù)據(jù)的投影類型為正弦曲線投影，時(shí)間分辨率為16 d，空間分辨率為250 m，全年共23個(gè)時(shí)相。本文選用的數(shù)據(jù)的空間時(shí)間范圍為2000年2月—2017年12月。渭河流域在SIN投影系統(tǒng)中的空間編號為H26V05，H27V05，因此共有原始影像822景。

對于獲取的NDVI數(shù)據(jù)，為便于ENVI，ArcGIS軟件讀取處理，以及與現(xiàn)有地圖數(shù)據(jù)實(shí)現(xiàn)空間上的匹配，需對其進(jìn)行拼接、投影轉(zhuǎn)換、格式轉(zhuǎn)換等預(yù)處理。本文利用MRT(Modis Reprojection Tools)工具進(jìn)行預(yù)處理。同時(shí)，由于影像數(shù)據(jù)會受到大氣、土壤以及高度角等因素的影響，導(dǎo)致NDVI值出現(xiàn)噪音。因此選擇最大值合成法(MVC,maximum Value Coposites)，對每月的兩期NDVI數(shù)據(jù)進(jìn)行合成，以獲取渭河流域2000—2017年月時(shí)序NDVI數(shù)據(jù)集[17]。

1.2 其他數(shù)據(jù)來源及處理

研究所用的各省省界區(qū)劃、模型計(jì)算公式來自其他地理信息網(wǎng)站。氣象數(shù)據(jù)來源于國家氣象信息中心(http:∥data.cma.cn/)，包括渭河流域及其周邊的27個(gè)氣象站點(diǎn)。選取的氣象因子為月平均氣溫和月降水量。通過對各種插值方法的實(shí)證比對研究，發(fā)現(xiàn)反距離權(quán)重插值法(IDW,Inverse Distance Weighing)精度較高，插值結(jié)果相對較好。因此選用反距離權(quán)重插值法將上述氣象要素插值為空間分辨率與NDVI一致的柵格數(shù)據(jù)。

2 研究方法

2.1 CA-Markov模型

元胞自動機(jī)(CA,Cellular Automata)是一種時(shí)間、空間、狀態(tài)均離散的局部網(wǎng)格動力學(xué)模型[18]。該方法具有強(qiáng)大的模擬復(fù)雜系統(tǒng)時(shí)空演化過程的能力，因而被廣泛應(yīng)用于土地變化、人口遷移、城市化等領(lǐng)域的研究[19-20]。元胞自動機(jī)模型可以表示為：

(1)

式中：Sij為元胞ij的狀態(tài)；f為轉(zhuǎn)換函數(shù)，定義了元胞從時(shí)刻t到t+1的轉(zhuǎn)換規(guī)則；n為元胞的鄰域，屬于f的一個(gè)輸入變量。應(yīng)用定義好的轉(zhuǎn)換規(guī)則就可以實(shí)現(xiàn)對局部空間范圍內(nèi)的元胞演化。

馬爾可夫模型(Markov Model)是具有無后效性的一種特殊隨機(jī)運(yùn)動過程，在許多地理現(xiàn)象的模擬與預(yù)測中得到應(yīng)用[21-23]。其基本原理是利用觀察系統(tǒng)中離散狀態(tài)之間的經(jīng)驗(yàn)傳遞概率以確定系統(tǒng)各狀態(tài)的變化趨勢，從而預(yù)測未來的狀態(tài)。在本文中，單像元尺度的NDVI狀態(tài)被看作是一個(gè)隨機(jī)過程，像元的NDVI值的等級被看作是Markov鏈的狀態(tài)。

CA模型和Markov模型均為時(shí)間離散、狀態(tài)離散的動力學(xué)模型，但Markov模型無法預(yù)測植被變化的空間變量，CA模型具有較強(qiáng)的模擬復(fù)雜空間系統(tǒng)時(shí)空演變的能力[24]。IDRISI 17.0軟件中的CA-Markov模塊將二者結(jié)合到一起，可以對植被動態(tài)變化的時(shí)空格局進(jìn)行模擬及預(yù)測。

2.2 Kappa系數(shù)

Kappa系數(shù)常用于評價(jià)遙感數(shù)據(jù)分類的精度、分析兩個(gè)圖件的一致性[25]。Kappa系數(shù)從空間位置和數(shù)量角度定量地闡明景觀變化過程中數(shù)量、位置和綜合信息的變化，其計(jì)算公式如下[26-27]：

(2)

式中：P0為正確模擬的比例；Pc為隨機(jī)情況下正確模擬比例的期望值。Kappa系數(shù)通常介于0～1。當(dāng)Kappa系數(shù)大于0.75時(shí)，模擬結(jié)果與真實(shí)分布之間的一致性較高，當(dāng)Kappa系數(shù)介于0.4～0.75時(shí)，一致性一般，當(dāng)Kappa系數(shù)小于0.4時(shí)，一致性較差。本文利用IDRISI軟件中提供的交叉表(Crosstable)模塊計(jì)算NDVI模擬結(jié)果和實(shí)際結(jié)果的Kappa系數(shù)，以評價(jià)模擬結(jié)果的精度。

3 結(jié)果與分析

3.1 NDVI空間分布格局模擬

為了驗(yàn)證CA-Markov模型在預(yù)測渭河流域植被動態(tài)變化中的可行性，本文以2005年和2010年的NDVI空間分布圖作為基礎(chǔ)數(shù)據(jù)，模擬2015年的NDVI分布，并與真實(shí)值作Kappa精度分析，以檢驗(yàn)?zāi)Ｐ驮谥脖蛔兓瘎討B(tài)預(yù)測中的可靠性。

利用CA-Markov模型進(jìn)行NDVI空間分布模擬的過程如下：

(1) NDVI數(shù)據(jù)分級。由于CA-Markov模型中的輸入數(shù)據(jù)需為空間、狀態(tài)均離散的柵格數(shù)據(jù)，因此，根據(jù)渭河流域NDVI的統(tǒng)計(jì)特征以及《土壤侵蝕分類分級標(biāo)準(zhǔn)》(SL190—2007)中土壤侵蝕強(qiáng)度分級指標(biāo)表，將其分為5個(gè)級別：高度植被覆蓋區(qū)(NDVI>0.6)、較高植被覆蓋區(qū)(0.5

圖1 渭河流域2005年、2010年NDVI分級

(2) 創(chuàng)建Markov狀態(tài)轉(zhuǎn)移概率矩陣。以2005年和2010年的平均NDVI空間分布圖為基礎(chǔ)，利用Markov模塊生成2005—2010年NDVI概率轉(zhuǎn)移矩陣(表1)和狀態(tài)轉(zhuǎn)變適宜性圖集(圖2)。轉(zhuǎn)移概率矩陣記錄了2005—2010年NDVI狀態(tài)發(fā)生轉(zhuǎn)移的概率(包括向自身轉(zhuǎn)移)，狀態(tài)轉(zhuǎn)變適宜性圖集顯示了在每一像元上找到不同的NDVI狀態(tài)的可能性。

表1定量分析了各等級NDVI之間的轉(zhuǎn)移方向以及不同等級NDVI之間的轉(zhuǎn)換概率。將植被覆蓋度由低到高的演變過程認(rèn)定為進(jìn)化演變，反之則認(rèn)定為退化演變[28]。由表1可以看出，2005—2010年，除高度植被覆蓋區(qū)之外，其他等級NDVI均屬于進(jìn)化演變。低度植被覆蓋區(qū)、較低植被覆蓋區(qū)、中度植被覆蓋區(qū)、較高植被覆蓋區(qū)進(jìn)化演變的概率分別為0.385 7，0.332 3，0.232 3，0.143 4。高度植被覆蓋區(qū)主要向較高植被覆蓋區(qū)退化演變，概率為0.178 5。

表1 渭河流域2005-2010年NDVI狀態(tài)轉(zhuǎn)移概率矩陣

(3) 應(yīng)用CA-Markov模型。在CA-Markov模塊中，以2010年的NDVI空間分布圖為基期數(shù)據(jù)，以2005—2010年NDVI狀態(tài)轉(zhuǎn)移概率矩陣和概率空間分布圖集作為轉(zhuǎn)換規(guī)則，并指定模型運(yùn)行周期數(shù)為5(即模擬的時(shí)間間隔)，運(yùn)行CA-Markov模型可以得到2015年NDVI空間分布模擬圖(圖3)。根據(jù)模擬結(jié)果(表2)，低度植被覆蓋區(qū)、較低植被覆蓋區(qū)面積有所減少，占比分別為13.03%，26.59%，空間上主要分布于流域西部和北部；中度植被覆蓋區(qū)、較高植被覆蓋區(qū)和高度植被覆蓋區(qū)面積有所增加，其中，中度植被覆蓋區(qū)主要分布在涇河和北洛河河谷地區(qū)，較高植被覆蓋區(qū)主要分布在關(guān)中平原、子午嶺等地，高度植被覆蓋區(qū)主要分布在秦嶺、六盤山。

表2 渭河流域2015年NDVI模擬結(jié)果

3.2 模擬結(jié)果精度驗(yàn)證

渭河流域2015年NDVI真實(shí)分布狀況如圖3B所示?？梢?，模擬結(jié)果和真實(shí)分布狀況具有較高的一致性。本文利用Kappa系數(shù)和ArcGIS疊置分析對模擬結(jié)果的精度進(jìn)行驗(yàn)證。

(1) Kappa系數(shù)。利用IDRISI軟件生成渭河流域2015年NDVI模擬結(jié)果和真實(shí)分布之間的混淆矩陣(單位為像元個(gè)數(shù))(表3)，同時(shí)計(jì)算二者之間的Kappa系數(shù)。結(jié)果表明，各等級NDVI模擬精度由大到小依次為：高度植被覆蓋區(qū)>低度植被覆蓋區(qū)>較低植被覆蓋區(qū)>中度植被覆蓋區(qū)>較高植被覆蓋區(qū)，其正確模擬的像元所占百分比分別為95.17%，73.92%，73.71%，67.08%，66.34%。計(jì)算得到渭河流域2015年NDVI模擬結(jié)果與真實(shí)分布之間的Kappa系數(shù)為0.785 0，表明模擬結(jié)果達(dá)到了較高的精度。

圖2 渭河流域2005-2010年NDVI狀態(tài)轉(zhuǎn)移適宜性圖集

圖3 渭河流域2015年NDVI模擬、實(shí)測結(jié)果

表3 渭河流域2015年NDVI模擬結(jié)果與真實(shí)分布之間的混淆矩陣 個(gè)

(2) 疊置分析。利用ArcGIS柵格計(jì)算器功能，將渭河流域2015年NDVI模擬結(jié)果與真實(shí)分布做柵格相減運(yùn)算。結(jié)果為0即為模擬準(zhǔn)確。經(jīng)統(tǒng)計(jì)，計(jì)算結(jié)果為0的像元共有1 555 990個(gè)，占像元總數(shù)的72.0%。利用隨機(jī)樣點(diǎn)生成工具(Create random points)在研究區(qū)內(nèi)生成1 000個(gè)空間位置隨機(jī)分布的樣點(diǎn)，并提取樣點(diǎn)所在位置的NDVI模擬的等級與真實(shí)值，通過作差運(yùn)算，得出模擬結(jié)果的準(zhǔn)確度為73.4%。表明模擬方法及過程可靠度較高。

3.3 NDVI空間分布格局預(yù)測

精度驗(yàn)證結(jié)果證明了CA-Markov模型在植被動態(tài)變化預(yù)測中具有較高的可信度。本文繼續(xù)以2010年、2015年渭河流域NDVI空間分布圖為基礎(chǔ)，建立CA-Markov模型，預(yù)測渭河流域2020年NDVI空間分布格局，進(jìn)而以2015年NDVI空間分布圖和2020年NDVI預(yù)測結(jié)果為基礎(chǔ)，預(yù)測渭河流域2025年NDVI空間分布格局。

渭河流域2020年和2025年NDVI空間分布格局預(yù)測結(jié)果見圖4和表4。根據(jù)預(yù)測結(jié)果，到2020年，高度植被覆蓋區(qū)和較高植被覆蓋區(qū)面積將繼續(xù)增加，分別比2015年增加6.52%，2.08%。空間上，高度植被覆蓋區(qū)擴(kuò)張的區(qū)域主要在子午嶺、六盤山等區(qū)域；較高植被覆蓋區(qū)擴(kuò)張的區(qū)域主要在關(guān)中平原。低植被覆蓋區(qū)、較低植被覆蓋區(qū)和中度植被覆蓋區(qū)面積有所減少，分別減少了4.45%，3.92%，0.23%；低度植被覆蓋區(qū)發(fā)生轉(zhuǎn)移的區(qū)域主要在隴中高原，主要流向是較低植被覆蓋區(qū)；較低植被覆蓋區(qū)發(fā)生轉(zhuǎn)移的區(qū)域主要在涇河流域和北洛河流域，主要流向是中度植被覆蓋區(qū)；中度植被覆蓋區(qū)的面積變化較小，其南部被高度植被覆蓋區(qū)和較高植被覆蓋區(qū)侵占，而北部又侵占了較低植被覆蓋區(qū)的面積。

圖4 渭河流域2020年、2025年NDVI預(yù)測結(jié)果

渭河流域2025年NDVI空間分布格局預(yù)測結(jié)果顯示，各等級NDVI的演變趨勢基本不變。到2025年，高度植被覆蓋區(qū)和較高植被覆蓋區(qū)的面積將分別增加4.14%，3.25%?？臻g上，高度植被覆蓋區(qū)擴(kuò)張的主要區(qū)域是關(guān)中平原，較高植被覆蓋區(qū)主要擴(kuò)張的區(qū)域是涇河和北洛河河谷地區(qū)。低度植被覆蓋區(qū)、較低植被覆蓋區(qū)和中度植被覆蓋區(qū)面積分別減少了1.84%，3.86%，1.69%?？臻g上，低度植被覆蓋區(qū)主要分布在流域北部和渭河沿岸的城鎮(zhèn)地區(qū)；較低植被覆蓋區(qū)一方面向北侵占低植被覆蓋區(qū)的面積，另一方面其南部被中度植被覆蓋區(qū)侵占；中度植被覆蓋區(qū)所在的區(qū)域整體上向北遷移。

表4 渭河流域2020年和2025年NDVI預(yù)測結(jié)果

4 討論與結(jié)論

(1) 渭河流域NDVI空間分布存在一定差異。低度植被覆蓋區(qū)、較低植被覆蓋區(qū)主要分布于隴中和陜北黃土高原地區(qū)，中度植被覆蓋區(qū)主要分布在涇河和北洛河河谷地區(qū)，較高植被覆蓋區(qū)主要分布在關(guān)中平原、子午嶺等地，高度植被覆蓋區(qū)主要分布在秦嶺、六盤山。

(2) CA-Markov模型對于模擬植被動態(tài)變化具有較好的適用性。利用CA-Markov模型對渭河流域植被的空間分布進(jìn)行模擬，結(jié)果表明模擬結(jié)果與NDVI真實(shí)分布之間Kappa系數(shù)達(dá)到0.785 0，疊置分析和隨機(jī)取樣檢驗(yàn)的精度分別為72.0%，73.4%，與原影像分級圖之間具有高度的一致性，說明在政策環(huán)境和氣候變化相對穩(wěn)定的條件下，利用CA-Markov模型定量地預(yù)測某一區(qū)域的植被覆蓋分布狀況的變化是可行的。

(3) 渭河流域NDVI時(shí)間變化規(guī)律明顯。整體上，2015—2025年，流域植被覆蓋狀況將進(jìn)一步改善。低度植被覆蓋區(qū)、較低植被覆蓋區(qū)、中度植被覆蓋區(qū)面積逐漸減少，向更高一級植被覆蓋區(qū)轉(zhuǎn)移，低度植被覆蓋區(qū)發(fā)生轉(zhuǎn)移的區(qū)域主要在流域西北部和渭河沿岸的城鎮(zhèn)地區(qū)；較低植被覆蓋區(qū)發(fā)生轉(zhuǎn)移的區(qū)域主要在涇河流域和北洛河流域，其一方面向北侵占低植被覆蓋區(qū)的面積，另一方面南部被中度植被覆蓋區(qū)侵占；中度植被覆蓋區(qū)整體向北遷移；較高植被覆蓋區(qū)和高度植被覆蓋區(qū)面積進(jìn)一步增加，空間上向北擴(kuò)張。

傳統(tǒng)的空間模擬及預(yù)測依據(jù)大量數(shù)據(jù)作回歸分析，以回歸方程估算屬性特征，其分析過程及模擬結(jié)果存在較多不確定性。本文利用CA-Markov模型，其工作原理是以預(yù)測基期的植被覆蓋空間分布為初始狀態(tài)，以基期和之后各級植被覆蓋區(qū)轉(zhuǎn)移面積及適宜性圖集表述的像元適宜NDVI為依據(jù)，對植被覆蓋空間分布進(jìn)行重新匹配，從而得到預(yù)測的各級植被覆蓋區(qū)面積，有效地模擬和預(yù)測了植被的時(shí)空分布特征。

需要指出的是，目前CA-Markov模型主要用于土地利用格局的模擬和預(yù)測，對于植被覆蓋的模擬和預(yù)測研究甚少。植被作為地表覆蓋的一種用地類型，其動態(tài)變化實(shí)質(zhì)上也是土地覆蓋的動態(tài)變化，因此研究嘗試?yán)肅A-Markov模型進(jìn)行植被覆蓋分布的模擬預(yù)測。CA-Markov模型的參數(shù)涉及時(shí)間間隔、數(shù)據(jù)尺度、濾波器等，對模型精度均有影響。本文選用最大值合成法對每月的兩期NDVI數(shù)據(jù)進(jìn)行合成，減少NDVI數(shù)據(jù)的噪音；另外針對不同濾波器參數(shù)下的模擬結(jié)果做了精度比較，發(fā)現(xiàn)5×5型濾波器模擬結(jié)果的準(zhǔn)確性最高，因此選用5×5型濾波器，提高模型模擬精度。

同時(shí)，植被覆蓋變化是一個(gè)復(fù)雜的過程，不僅受到氣候變化、自然災(zāi)害等諸多自然因素的影響，而且受到土地利用政策、社會經(jīng)濟(jì)發(fā)展及其他人類活動等不確定因素的影響，因此，如何綜合考慮多種因素設(shè)置CA-Markov模型的各參數(shù)是今后需要深入探究的科學(xué)問題。