基于CA-Markov的鄖西縣城關(guān)鎮(zhèn)土地利用變化圖譜的模擬及預(yù)測(cè)

唐宇佳，王 浩，李 禹，金 鑫

(1.湖北大學(xué)，湖北 武漢 430062；2.湖北一二規(guī)劃設(shè)計(jì)有限公司，湖北 武漢 430062)

1 引言

目前，對(duì)LUCC課題的相關(guān)研究仍然很多[1～5]。在LUCC研究熱潮中，一部分人認(rèn)為土地變化不僅影響到生態(tài)、社會(huì)經(jīng)濟(jì)，而且也影響到了全球經(jīng)濟(jì)發(fā)展和生態(tài)保護(hù)[6，7]。眾多學(xué)者研究成果表明：土地利用的時(shí)空變化的研究已經(jīng)越來(lái)越重要，已然成為眾多學(xué)者的首選研究課題。

自該研究熱點(diǎn)研究以來(lái)，學(xué)者采用了動(dòng)力學(xué)模型[8～10]、CLUE-S模型[11，12]、Dyna-CLUES模型[13～15]、人工神經(jīng)網(wǎng)絡(luò)模型(Artificial Neural Network，ANN)[16，17]和元胞自動(dòng)機(jī)-馬爾科夫鏈(CA-Markov)模型[18～21]等來(lái)模擬并預(yù)測(cè)土地利用的時(shí)空變化情況。系統(tǒng)動(dòng)力學(xué)模型是通過(guò)系統(tǒng)論來(lái)研究土地變化的因素，但是其軟件和GIS的空間分析功能融合并不好，使其作用沒(méi)有完全發(fā)揮出來(lái)；CLUE-S模型可以對(duì)時(shí)間和空間進(jìn)行雙重模擬，但其系統(tǒng)中的參數(shù)主要來(lái)源于人主觀設(shè)置，缺乏客觀性和科學(xué)性；Dyna-CLUE模型是在CLUE、CLUE-S模型基礎(chǔ)上發(fā)展而來(lái)，該模型綜合了宏觀和微觀兩個(gè)方面的，對(duì)多尺度的土地利用變化有很好的適用性；CA-Markov模型通過(guò)發(fā)揮Markov精確的預(yù)測(cè)能力和CA強(qiáng)大的空間模擬能力，將其兩者結(jié)合來(lái)預(yù)測(cè)和模擬研究區(qū)土地利用變化趨勢(shì)，既可以模擬空間變化，又可以預(yù)測(cè)其準(zhǔn)確面積。相比其他模型，該模型克服了以往模型模擬空間能力的缺陷和預(yù)測(cè)精度問(wèn)題。目前，學(xué)者常用CA-Markov模型來(lái)模擬和預(yù)測(cè)土地利用變化。

利用圖譜變化來(lái)研究土地利用的時(shí)空變化逐漸受到廣大學(xué)者的關(guān)注[22～25]。本文對(duì)城關(guān)鎮(zhèn)9年土地變化情況進(jìn)行圖譜分析，然后通過(guò)CA-Markov模型來(lái)模擬及預(yù)測(cè)城關(guān)鎮(zhèn)土地的時(shí)空變化，預(yù)測(cè)結(jié)果可以給國(guó)土和規(guī)劃部門提供幫助。

鄖西縣城關(guān)鎮(zhèn)地形以山地為主，植被、水流資源豐富，是國(guó)家新階段扶貧開(kāi)發(fā)重點(diǎn)縣，南水北調(diào)中線工程水源保護(hù)區(qū)上游核心區(qū)。近幾年，城市范圍擴(kuò)大，占用大量耕地，城區(qū)的發(fā)展對(duì)河流、植被等生態(tài)造成了不同程度上的破壞。如何合理布局建設(shè)用地，減少耕地占有量是解決這一問(wèn)題的重中之重[26]。

本文將城關(guān)鎮(zhèn)3年的土地現(xiàn)狀矢量數(shù)據(jù)根據(jù)國(guó)家國(guó)土資源部的分類標(biāo)準(zhǔn)分為耕地、林地、草地、水域、建設(shè)用地和未利用地6類，運(yùn)用ArcGIS10.2對(duì)幾種地類的空間布局和數(shù)量進(jìn)行分析統(tǒng)計(jì)；選擇CA-Markov模型模擬和預(yù)測(cè)鄖西縣城關(guān)鎮(zhèn)2017年土地利用變化情況，并將其和實(shí)際2017年土地利用情況進(jìn)行Kappa系數(shù)分析；結(jié)合土地轉(zhuǎn)移矩陣和適宜性圖集來(lái)模擬并預(yù)測(cè)城關(guān)鎮(zhèn)2020年的土地變化情況。預(yù)測(cè)結(jié)果可以給國(guó)土和規(guī)劃部門提供幫助。

2 研究區(qū)概況

鄖西縣城關(guān)鎮(zhèn)地處鄖西縣東部，北與土門鎮(zhèn)相鄰，南與觀音、河夾相鄰。秦嶺南坡漢水北岸，素有東西通“豫蜀”，南北連秦、楚之稱(圖1)。全鎮(zhèn)版圖面積127.83 km2，轄21個(gè)行政村、173個(gè)組、14個(gè)社區(qū)居委會(huì)、202個(gè)居民小組，總?cè)丝?48612人，是鄖西縣的政治、經(jīng)濟(jì)、文化、旅游及商貿(mào)中心。2017年末，土地總面積約為127.83 km2，土地分類情況如下：①耕地及林地比例較大，占土地總面積的82.76%，面積分別為17.63 km2，88.17 km2；②建設(shè)用地面積為12.29 km2，占土地總面積的9.61%；③水域和草地占比相當(dāng)，面積分別為2.12 km2和5.04 km2；④未利用地比例很小，只占0.09%，面積為0.11 km2。鄖西作為南水北調(diào)核心水源區(qū)，保證水源區(qū)周邊的生態(tài)環(huán)境至關(guān)重要，如何發(fā)揮“河長(zhǎng)”的角色，破解河道管理、水環(huán)境整治中長(zhǎng)期存在的困局，做到水系完整、河流暢通、水質(zhì)良好、生物多樣、岸線優(yōu)美。如何采取有效措施保證“一江清水送北京”。此外，國(guó)家最近實(shí)施的精準(zhǔn)扶貧政策工作已然拉起大幕，如何確保易地搬遷戶能夠種到好田，吃到無(wú)污染的蔬菜、糧食。如何合理安排布局安置房和保障房才能實(shí)現(xiàn)“真扶貧、扶真貧”的目的等一系列問(wèn)題值得思考。

3 數(shù)據(jù)源及研究方法

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

3.1.1 數(shù)據(jù)來(lái)源

本研究以城關(guān)鎮(zhèn)2009、2013和2017年土地利用現(xiàn)狀矢量數(shù)據(jù)和DEM數(shù)據(jù)作為主要數(shù)據(jù)源，其來(lái)源主要是國(guó)土規(guī)劃部門。與城關(guān)鎮(zhèn)相關(guān)部門有項(xiàng)目合作并簽訂了數(shù)據(jù)使用保密協(xié)議后才拿到數(shù)據(jù)。相比其他學(xué)者使用遙感影像解譯得到的土地利用現(xiàn)狀數(shù)據(jù)，直接使用土地利用現(xiàn)狀矢量數(shù)據(jù)可以減少監(jiān)督分類和人工目視解譯時(shí)的誤差，精度有所提高。

3.1.2 數(shù)據(jù)預(yù)處理

將3年的土地利用現(xiàn)狀數(shù)據(jù)根據(jù)國(guó)土部門發(fā)布的土地分類標(biāo)準(zhǔn)進(jìn)行劃分，并提取出需要的地類，對(duì)其進(jìn)行編號(hào)。在ArcGIS10.2中對(duì)各個(gè)地類面積進(jìn)行統(tǒng)計(jì)；將收集到的DEM數(shù)據(jù)處理后提取相關(guān)坡度值。

圖1 城關(guān)鎮(zhèn)區(qū)域位置

3.2 研究方法

3.2.1 土地利用轉(zhuǎn)移矩陣

土地利用轉(zhuǎn)移矩陣是通過(guò)計(jì)算研究期初始和末的土地面積變化來(lái)得出的矩陣，通過(guò)對(duì)矩陣的分析，能夠研究各個(gè)地類在研究期內(nèi)的變化情況。其表達(dá)式為：

(1)

式(1)中A表示各種地類的面積，n表示當(dāng)前研究的地類，i表示研究初期地類，j表示研究期末地類。

3.2.2CA模型

研究某種地類在t+1時(shí)刻的地類情況，可通過(guò)研究其相鄰單元細(xì)胞和自身t時(shí)刻地類來(lái)計(jì)算得到。結(jié)合土地轉(zhuǎn)移矩陣可以計(jì)算出更復(fù)雜的土地變化情況。其表達(dá)式為：

S(t+1)=f(S(t),N)

(2)

式(2)中t和t+1為單元細(xì)胞所處的兩個(gè)時(shí)刻；S為單元細(xì)胞的狀態(tài)集合；f為單元細(xì)胞轉(zhuǎn)化函數(shù)；N為相鄰單元細(xì)胞。

3.2.3Markov模型

Markov土地利用變化模型是利用概率建立的一種隨機(jī)型的時(shí)序模型，使用該模型的步驟如下：①根據(jù)3年土地利用現(xiàn)狀數(shù)據(jù)推算出各種地類的轉(zhuǎn)出和轉(zhuǎn)入率，以此來(lái)構(gòu)建土地利用變化轉(zhuǎn)移矩陣；②統(tǒng)計(jì)預(yù)測(cè)初始年的土地利用情況；③通過(guò)IDRISI軟件中的馬爾科夫模型來(lái)預(yù)測(cè)2020年城關(guān)鎮(zhèn)土地利用情況。表達(dá)式為：

S(t+1)=Pij·S(t)

(3)

(4)

3.2.4 CA-Markov模型

本文將CA和Markov模型進(jìn)行集成，以2009年和2013年土地利用現(xiàn)狀數(shù)據(jù)為基礎(chǔ)，來(lái)模擬2017年土地利用情況，并計(jì)算其Kappa系數(shù)。在精度要求范圍內(nèi)，再以2017年為模擬初始年，將由驅(qū)動(dòng)因子制成的適宜性圖集和土地轉(zhuǎn)移概率矩陣輸入到模型中，來(lái)預(yù)測(cè)模擬2020年城關(guān)鎮(zhèn)土地利用情況。具體的步驟如下。

(1)轉(zhuǎn)移矩陣。將研究區(qū)3年的土地利用現(xiàn)狀矢量數(shù)據(jù)轉(zhuǎn)化為柵格形式，象元大小為100×100。時(shí)間周期設(shè)置為8a，比例系數(shù)設(shè)置為0.15，通過(guò)IDRISI軟件的Markov模塊來(lái)計(jì)算得到城關(guān)鎮(zhèn)這8年間的土地利用轉(zhuǎn)移概率矩陣。

(2)適宜圖集。結(jié)合研究區(qū)實(shí)際情況和主要研究地類，把水域和建設(shè)用地設(shè)置為主要限制因子，將坡度、距離等設(shè)置為限制條件。將土地利用柵格圖和限制條件、因子等參數(shù)輸入到COLLETION EDIT模塊中生成適宜性圖集。

(3)CA濾波器。參考其他研究學(xué)者的研究成果并結(jié)合研究區(qū)實(shí)際情況，該研究采取5×5的濾波器。

(4)確定CA循環(huán)次數(shù)。本文以 2009年土地利用格局為起始年，CA循環(huán)次數(shù)設(shè)為7。

3.2.5 Kappa系數(shù)預(yù)測(cè)精度檢驗(yàn)

為了驗(yàn)證CA-Markov模型預(yù)測(cè)精度，在分類的精度評(píng)價(jià)中，不同的精度評(píng)價(jià)方法有不同的劃分標(biāo)準(zhǔn)。筆者主要采用IDRISI軟件中CROSSTAB工具計(jì)算的overall Kappa系數(shù)來(lái)評(píng)價(jià)預(yù)測(cè)圖和原始圖之間的一致性。因此直接借用Cohen提出的Kappa系數(shù)分類評(píng)價(jià)標(biāo)準(zhǔn)。如表1所示。

4 結(jié)果分析

4.1 城關(guān)鎮(zhèn)土地利用分析

4.1.1 2009～2017年城關(guān)鎮(zhèn)土地利用變化分析

利用ArcGIS10.2統(tǒng)計(jì)得到城關(guān)鎮(zhèn)三期土地利用面積如表2所示。其土地利用現(xiàn)狀如圖2所示。

表1 Kappa系數(shù)分類標(biāo)準(zhǔn)

表2 城關(guān)鎮(zhèn)2009、2013和2017年土地利用變化

圖2 城關(guān)鎮(zhèn)2009、2013及2017年土地利用現(xiàn)狀

由表2可知，2009～2017年城關(guān)鎮(zhèn)耕地和林地減少較多，水域和建設(shè)用地增加較多。究其原因，南水北調(diào)工程水源區(qū)工程的完工，導(dǎo)致水域面積增加，易地扶貧搬遷導(dǎo)致建設(shè)用地增加。

從數(shù)字上分析，3年林地面積占比分別為69.03%、69.01%和68.97%，其面積呈現(xiàn)出減少趨勢(shì)，原因是作為南水北調(diào)核心水源區(qū)，城關(guān)鎮(zhèn)擔(dān)任著保護(hù)好水源水質(zhì)的重?fù)?dān)，修建堤壩、防風(fēng)、加固壩等工程均會(huì)占用林地。耕地面積由2009年的1787.70 hm2減少到2017年的1762.58 hm2，占比從13.99%減少到13.79%，面積減少25.12 hm2，占比減少0.2%。8a期間草地面積減少8.52 hm2，占比減少0.07%。由于蓄水量逐年增加，水域面積由2009年的212.75 hm2增加到2017年的214.91 hm2，占比增加0.02%。建設(shè)用地面積由2009年的1196.81 hm2增加到2017年的1228.87 hm2，增加了32.06 hm2。

4.1.2 2009～2017年城關(guān)鎮(zhèn)土地利用變化圖譜分析

根據(jù)式(2)，借助于ArcGIS10.2求出城關(guān)鎮(zhèn)2009～2017年土地利用變化圖譜如表3所示。

由表3可知，城關(guān)鎮(zhèn)2009～2017年期間土地利用變化類型中，穩(wěn)定變化型面積最大，其變化圖斑為25756個(gè)，變化面積為18882.15 hm2，后期變化型次之，變化圖斑為1307個(gè)，變化面積為160.3 hm2，持續(xù)性變化型面積最小，其變化圖斑為17個(gè)，變化面積為0.19 hm2。2009～2017年城關(guān)鎮(zhèn)土地利用變化面積由大到小排序?yàn)榉€(wěn)定型、后期型、前期型、反復(fù)型、持續(xù)型。

表3 城關(guān)鎮(zhèn)2009～2017年土地利用變化圖譜分析

4.1.3 2009～2017年城關(guān)鎮(zhèn)土地利用流向分析

本節(jié)利用ArcGIS10.2空間分析工具計(jì)算土地利用轉(zhuǎn)移矩陣，城關(guān)鎮(zhèn)2009～2013年及2013～2017年的土地利用變化轉(zhuǎn)移面積矩陣如表4、表5所示。

由表3、4可知，在2009～2013年城關(guān)鎮(zhèn)土地利用轉(zhuǎn)移面積矩陣中，轉(zhuǎn)移面積最大的是水域向建設(shè)用地的轉(zhuǎn)移，為84.03 hm2，其次是耕地向建設(shè)用地的轉(zhuǎn)移，為15.27 hm2。究其原因是南水北調(diào)初期，堤壩等加固設(shè)施沒(méi)有完善，再加上城市擴(kuò)展所致。

表4 城關(guān)鎮(zhèn)2009～2013年土地利用變化轉(zhuǎn)移面積矩陣

表5 城關(guān)鎮(zhèn)2013～2017年土地利用變化轉(zhuǎn)移面積矩陣

4.2 城關(guān)鎮(zhèn)土地利用演變模擬

4.2.1 城關(guān)鎮(zhèn)土地利用轉(zhuǎn)移概率矩陣

將三年的土地利用現(xiàn)狀矢量數(shù)據(jù)加載到ArcGIS10.2軟件中，計(jì)算出8年各個(gè)研究地類的轉(zhuǎn)入率和轉(zhuǎn)出率，然后根據(jù)式(3)(4)，得到2009～2013和2013～2017年土地利用轉(zhuǎn)移概率矩陣，如表6，表7所示。

4.2.2 2017年城關(guān)鎮(zhèn)土地利用適宜性圖集準(zhǔn)備

鄖西縣地處鄂西北，地形以山地為主，導(dǎo)致城關(guān)鎮(zhèn)的建設(shè)只能以組團(tuán)方式發(fā)展；自南水北調(diào)動(dòng)工以來(lái)，核心水源區(qū)加固工程將建成并開(kāi)始蓄水，該地區(qū)的耕地和林地急劇減少，水域及建設(shè)用地急劇增加。結(jié)合研究區(qū)實(shí)際情況和主要研究地類，把水域和建設(shè)用地設(shè)置為主要限制因子，將坡度、距離等設(shè)置為限制條件。選擇基于布爾運(yùn)算的MCE制作各地類適宜性圖集。

4.2.3 2017年城關(guān)鎮(zhèn)土地利用演變模擬

將2013年土地利用情況作為研究初年，將土地利用轉(zhuǎn)移概率矩陣、各個(gè)參數(shù)和適宜性圖集輸入到CA-Markov模塊中，循環(huán)次數(shù)設(shè)置為4，然后運(yùn)行該模塊得到2017年土地利用的空間布局。如圖3。將模擬的城關(guān)鎮(zhèn)2017年土地利用數(shù)據(jù)和實(shí)際2017年土地利用數(shù)據(jù)進(jìn)行Kappa系數(shù)分析。分析統(tǒng)計(jì)表如表8所示。

由于城關(guān)鎮(zhèn)地處山區(qū)，地形比較復(fù)雜，海拔較高，另外由于南水北調(diào)工程的影響和易地扶貧搬遷政策等因素決定了未來(lái)多年內(nèi)該區(qū)域土地利用的供給情況和空

表6 2009～2013年土地利用變化轉(zhuǎn)移概率矩陣

表7 2013～2017年土地利用變化轉(zhuǎn)移概率矩陣

間分布情況。從表8可知，模擬的各個(gè)地類中林地、耕地和建設(shè)用地面積較大，這和實(shí)際的2017年城關(guān)鎮(zhèn)土地分布情況是一致的。為了檢驗(yàn)?zāi)M精度，文章利用IDRISI軟件中Kappa精度來(lái)檢測(cè)模擬圖和實(shí)際圖之間的模擬精度，經(jīng)統(tǒng)計(jì)，實(shí)際圖和模擬圖的Kappa精度為93%，證明了用該模型模擬城關(guān)鎮(zhèn)土地利用變化是合適的。

圖3 城關(guān)鎮(zhèn)2017年土地利用現(xiàn)狀模擬圖及實(shí)際圖

土地利用類型2017年實(shí)際土地利用結(jié)構(gòu)2017年模擬土地利用結(jié)構(gòu)面積對(duì)比面積/hm2占比/%面積/hm2占比/%面積差/hm2精度/%耕地1762.5813.791460.5811.6530294.57林地8816.6368.977947.6363.3986991.24草地504.233.94326.2326.0617893.56水域214.911.68237.9155.33-2394.78建設(shè)用地1228.879.613151.8725.14-192395.88未利用地11.020.0915.023.94-491.67

4.3 城關(guān)鎮(zhèn)土地利用演變預(yù)測(cè)

根據(jù)Kappa系數(shù)可知模擬和實(shí)際的模擬精度值為93%，說(shuō)明該模型在模擬和預(yù)測(cè)該區(qū)域的土地利用情況是合適的。將土地利用轉(zhuǎn)移概率矩陣、各個(gè)參數(shù)和適宜性圖集輸入到CA-Markov模塊中，循環(huán)次數(shù)設(shè)置為4，然后運(yùn)行該模塊得到2017年土地利用的空間布局。如圖4和表9所示。

由表9可知，到2020年，城關(guān)鎮(zhèn)土地利用繼續(xù)呈現(xiàn)出耕地與林地的大量減少、水域與建設(shè)用地的迅速增加、以及草地與未利用地不明顯變化等特點(diǎn)。耕地面積由2017年的1762.58 hm2減少到2020年的1678.58 hm2，所占百分比由13.79%下降到13.39%，面積減少84 hm2，占比減少0.4%。為了防止耕地?zé)o限制減少，政府應(yīng)當(dāng)嚴(yán)格執(zhí)行耕地保護(hù)紅線，牢牢保護(hù)好劃定的永久基本農(nóng)田，控制和阻止盲目開(kāi)墾耕地的行為，以此促進(jìn)該區(qū)域的生態(tài)環(huán)境的保護(hù)和改善。林地面積由2017年的8816.63 hm2減少到2020年的6781.63 hm2，所占百分比由68.97%下降到54.09%，面積減少2035 hm2，占比減少14.88%。林地的變化強(qiáng)度大于2009～2017年，林地是城關(guān)鎮(zhèn)主要的用地類型，為了保護(hù)該區(qū)域的生態(tài)環(huán)境及社會(huì)經(jīng)濟(jì)的可持續(xù)發(fā)展，政府應(yīng)嚴(yán)格控制林地的減少量。水域面積由2014年的214.91 hm2增加到2020年的456.91 hm2，所占百分比由1.68%上升到3.64%，面積增加242 hm2，占比增加1.96%。南水北調(diào)竣工以后，水源區(qū)不斷蓄水使得該區(qū)域社會(huì)效益與經(jīng)濟(jì)效益都有一定程度的增加，但是由于水深增加，水流速度減慢，水域由一開(kāi)始的自凈化能力較強(qiáng)的天然河流變?yōu)樯鷳B(tài)環(huán)境功能及其脆弱的人工湖泊，干支流納污能力減弱，水體富營(yíng)養(yǎng)化嚴(yán)重。該區(qū)域應(yīng)關(guān)注水體質(zhì)量，避免水生態(tài)環(huán)境遭受破壞，污染水源。建設(shè)用地面積由2017年的1228.87 hm2增加到2020年的2346.87 hm2，所占百分比由9.61%上升到18.72%，面積增加1118 hm2，占比增加9.11%。該區(qū)域在制定土地利用總體規(guī)劃時(shí)，應(yīng)當(dāng)嚴(yán)格控制建設(shè)用地的增加量，嚴(yán)格禁止一切違法建設(shè)行為。

圖4 城關(guān)鎮(zhèn)2020年土地利用現(xiàn)狀預(yù)測(cè)

土地利用類型2017年解譯數(shù)據(jù)2017年預(yù)測(cè)數(shù)據(jù)面積/hm2占比/%面積/hm2占比/%耕地1762.5813.791678.5813.39林地8816.6368.976781.6354.09草地504.233.94506.234.04水域214.911.68456.913.64建設(shè)用地1228.879.612346.8718.72未利用地11.020.0931.020.25

本文選擇的研究區(qū)域具有一定的代表性，能夠?yàn)樵搮^(qū)域的國(guó)土規(guī)劃部門提供一定的理論支撐。但由于選擇的驅(qū)動(dòng)因子不完整，選擇的因子并不能夠真實(shí)影響到土地利用變化情況。在以后的研究中應(yīng)當(dāng)根據(jù)實(shí)際情況選擇具有說(shuō)服力、具有真實(shí)性的因子。