劉建華 張啟斌 YANG Di 岳德鵬 于 強 楊 斕

(1.北京林業(yè)大學精準林業(yè)北京市重點實驗室， 北京 100083； 2.弗羅里達大學地理系， 蓋恩斯維爾 FL 32611)

0 引言

《全國生態(tài)環(huán)境保護綱要(國發(fā)[2000]38號)》中提到了生態(tài)用地一詞，明確指出了生態(tài)用地的重要性[1]。眾多學者基于此不斷發(fā)展生態(tài)用地的概念和內(nèi)涵。一種觀點認為應該從土地類型上來劃分生態(tài)用地，即提供自然生態(tài)系統(tǒng)服務價值的土地都可以被視為生態(tài)用地。另一種觀點認為應該從土地的主體功能來區(qū)分，主張農(nóng)業(yè)用地不應作為生態(tài)用地[2-3]。根據(jù)前人理論成果，結合包頭市實地情況，本文將林地、草地、水體作為生態(tài)用地范圍。

生態(tài)用地作為自然生態(tài)系統(tǒng)服務的基本載體，是解決城市建設用地擴張與生態(tài)保護矛盾的綜合途徑[4]。在我國西北干旱地區(qū)，生態(tài)用地在防止土地沙漠化、水土保持等方面的生態(tài)功能更加突出，對于維持區(qū)域生態(tài)系統(tǒng)健康穩(wěn)定具有至關重要的作用[5-6]。對生態(tài)用地的空間變化進行模擬，可以為干旱地區(qū)生態(tài)用地科學規(guī)劃和生態(tài)安全格局構建提供重要參考價值。

生態(tài)用地的演變模擬是土地利用演變模擬的一部分，當前土地利用格局的演變模擬主要包括兩種類型，一種是宏觀回顧模型，如邏輯回歸(Logistic regression)模型[7]、馬爾可夫(Markov)模型[8]、系統(tǒng)動力學(System dynamic)模型等[9-10]，這些模型根據(jù)土地利用的時空數(shù)據(jù)，量化土地利用在宏觀時空尺度上的變化規(guī)律，對未來情景進行預測。另一種是微觀預測模型，包括多智能體模型(Multi-agent)[11]、元胞自動機(Cellular automata, CA)[12]等，這種模型從微觀尺度入手，通過研究土地利用與生態(tài)、人類活動等因子的相互作用，預測宏觀上土地利用的未來演變。將兩種模型結合，可得到二者的混合模型，既可以考慮全局宏觀效應，又可以考慮多個因子對土地利用的微觀作用機制，相比單一模型具有較大優(yōu)勢。將宏觀回顧模型與元胞自動機結合，是這種模型的一種常見形式，如CA-Markov模型[13]。在這種模型中，有兩個關鍵點對其模擬精度有較大影響，一是元胞鄰域內(nèi)不同土地利用類型間轉變規(guī)則，即鄰域規(guī)則；二是地形、生態(tài)環(huán)境等因子對土地利用演變影響，即適宜性規(guī)則[14]。當前傳統(tǒng)元胞自動機鄰域規(guī)則制定方法較多，如支持向量機、多目標決策、Markov模型等[15]，用戶往往難以選擇合適的規(guī)則定義方法，部分線性模型也難以根據(jù)鄰域特征精確模擬生態(tài)用地演化等非線性過程；基于邏輯回歸、多標準評價模型及層次分析法等制定的適宜性規(guī)則雖然考慮了多種因素對土地利用變化的影響，但多基于單個像元利用空間疊加手段獨立評價其適宜性，不能根據(jù)生態(tài)用地演化的過程量化研究區(qū)內(nèi)各景觀單元演化為生態(tài)用地的適宜性。

針對上述問題，選取西北干旱地區(qū)典型城市包頭市為研究區(qū)，基于源匯理論，利用最小累積耗費阻力(Minimum cumulative resistance, MCR)模型量化生態(tài)用地演化過程中，從“生態(tài)源地”到其他土地利用類型的適宜性，同時利用人工神經(jīng)網(wǎng)絡(Artificial neural network, ANN)提取CA鄰域轉換規(guī)則，構建MCR-ANN-CA模型對包頭市生態(tài)用地的演變過程進行模擬，以期為區(qū)域生態(tài)用地規(guī)劃及生態(tài)建設提供理論與方法支持。

1 材料與方法

1.1 研究區(qū)概況

包頭市位于內(nèi)蒙古自治區(qū)西部(東經(jīng)109°13′～111°26′，北緯40°13′～42°44′)，面積27 768 km2。包頭市深處內(nèi)陸，氣候為典型的溫帶干旱、半干旱大陸性氣候，冬季寒冷干燥、夏季炎熱多雨[16]，年平均氣溫2.0～7.7℃，年均降水量175～400 mm，年均蒸發(fā)量為2 100～2 700 mm。包頭市可利用地表水總量為9×108m3，地下水補給量為8.6×109m3。黃河流經(jīng)包頭境內(nèi)214 km，水面寬130～458 m，最大流量6 400 m3/s，年平均徑流量為2.6×1011m3，是包頭市主要用水來源[17]。包頭市海拔976～2 317 m，地勢中間高南北低，北部丘陵、中部山地、南部平原分別占總面積的14.49%、75.51%和10%[18]。干旱的氣候條件與起伏的地貌特征使得包頭市生態(tài)環(huán)境較為脆弱，面臨較高的水土流失與土地沙漠化、荒漠化風險。近年來，包頭市城市建設規(guī)模擴張迅速，房地產(chǎn)開發(fā)項目、工業(yè)園區(qū)等建設項目不斷涌現(xiàn)，導致草地、林地、濕地等生態(tài)用地遭到占用與破壞，生態(tài)風險有所升高[19]。

1.2 數(shù)據(jù)來源與處理

本研究主要數(shù)據(jù)源包括包頭市土地利用數(shù)據(jù)(2006、2011、2016年)、歸一化植被指數(shù)(Normalized difference vegetation index, NDVI)、數(shù)字地面高程模型(Digital elevation model, DEM)、各區(qū)縣旗人口密度及工業(yè)園區(qū)分布數(shù)據(jù)。其中土地利用數(shù)據(jù)由Landsat-8遙感數(shù)據(jù)解譯得到，工業(yè)園區(qū)分布數(shù)據(jù)來自《包頭市土地利用總體規(guī)劃》和《包頭市城市總體規(guī)劃》；歸一化植被指數(shù)(NDVI)、數(shù)字地面高程模型(DEM)來自地理空間數(shù)據(jù)云(http:∥www.gscloud.cn)；各區(qū)人口密度數(shù)據(jù)來自包頭市統(tǒng)計年鑒。

根據(jù)研究需要，從遙感影像中提取耕地、林地、草地、建設用地、水體、未利用地并通過實地驗證確保其精度，其中林地、草地、水體為本文中生態(tài)用地范圍；利用DEM數(shù)據(jù)提取高程、坡度數(shù)據(jù)。

1.3 最小累積耗費阻力模型

MCR模型由KANPPEN提出，最早應用于物種遷徙過程研究，之后在物種保護、景觀格局分析等方面取得了廣泛應用[20]。該模型主要考慮3個因素，即“源”、阻力和累積代價，通過對3個因素的分析，對“源”克服阻力向外傳播所耗費的代價或者所做的功進行描述[21]。MCR模型的一般形式為

R=f∑DijRi

(1)

式中R——最小累積阻力

f——未知負函數(shù)，表示最小累積阻力與生態(tài)適宜性的負相關關系

Dij——從源j到景觀單元i的空間距離

Ri——景觀單元i處的阻力

本文將生態(tài)用地的演化過程看作生態(tài)用地對其他景觀的競爭性控制過程，且這種演化必須通過克服阻力實現(xiàn)，這樣生態(tài)用地的演化過程就可以抽象為從源(現(xiàn)有生態(tài)用地斑塊)到匯(其他景觀單元)克服阻力做功的水平過程[22]。由于區(qū)域下墊面差異，不同空間位置的土地演化為生態(tài)用地的阻力是不同的，通過由“源”到當前像元的累積阻力可量化當前像元演化為生態(tài)用地的概率，即阻力越大，該像元演化為生態(tài)用地的概率越小，演化為其他用地類型的概率越大。本文利用MCR模型，綜合考慮土地利用、NDVI、坡度、政府規(guī)劃工業(yè)園區(qū)、人口密度、水體距離、高程因子構建累積耗費阻力面，利用該累積耗費阻力面構建CA模型的適宜性規(guī)則，對上述過程進行模擬，提高模型預測精度。

1.4 元胞自動機模型

CA模型是一種時間、空間、狀態(tài)都離散的動力學模型，具有明顯的時間與空間特征，CA模型的一般形式為[23]

Si,t+1=f(Si,t,SN,t)

(2)

式中Si,t+1——元胞i在t+1時刻的狀態(tài)

Si,t——元胞i在t時刻的狀態(tài)

SN,t——元胞i的鄰域集合在t時刻的狀態(tài)

f——轉換規(guī)則

CA模型中，中心元胞在下一時刻的狀態(tài)是其在當前時刻狀態(tài)及其鄰域集合狀態(tài)的函數(shù)，準確定義該函數(shù)對于CA模型的模擬精度具有關鍵作用[24-27]。由于生態(tài)用地的演化過程是一種非線性的復雜動力學過程，當前元胞及其鄰域狀態(tài)影響中心元胞過程很難用簡單的規(guī)則定義，因此本文采用人工神經(jīng)網(wǎng)絡模型提取CA模型轉換規(guī)則，對土地利用模擬過程中的現(xiàn)有規(guī)則進行改進[28-29]，以提高CA模型的模擬精度。

本研究中，元胞形狀為正方形的柵格像元，尺寸為30 m×30 m，采用3×3經(jīng)典摩爾鄰域定義鄰域空間，模型初始柵格數(shù)據(jù)中的元胞狀態(tài)定義為研究區(qū)景觀格局類型，依次為耕地、林地、草地、建設用地、水體、未利用地，根據(jù)中心元胞及鄰域元胞狀態(tài)，利用ANN模型提取元胞轉換規(guī)則，利用該規(guī)則計算元胞演化為生態(tài)用地的概率。

1.5 人工神經(jīng)網(wǎng)絡模型

ANN模型是20世紀80年代以來人工智能領域興起的研究熱點[30]。ANN模型在復雜的非線性系統(tǒng)的模擬方面具有明顯優(yōu)勢，其自組織、自學習、聯(lián)想以及記憶的優(yōu)勢能夠有效簡化CA模型，從原始訓練數(shù)據(jù)中提取CA轉換規(guī)則，避免主觀因素影響，提高模擬精度[31]。

由于ANN模型在解決此類非線性問題中的明顯優(yōu)勢，國內(nèi)外學者很早就開始了相關研究，嘗試利用ANN模型提取CA模型的轉換規(guī)則[32]。此類研究往往把鄰域規(guī)則和適宜性規(guī)則統(tǒng)一放入ANN模型進行規(guī)則提取，然而在模擬生態(tài)用地演化的CA模型中，相比適宜性規(guī)則，鄰域規(guī)則部分更難定義也更難被人類理解，因此本文利用ANN模型提取CA模型鄰域部分轉換規(guī)則，適宜性規(guī)則采用MCR模型構建。

采用經(jīng)典BP神經(jīng)網(wǎng)絡算法，網(wǎng)絡結構共3層，分別是輸入層、隱含層、輸出層。模型的輸入層共7個節(jié)點，隱含層共5個節(jié)點，采用tansing激勵函數(shù)，輸出層共1個節(jié)點，采用sigmoid激勵函數(shù)，其輸出為中心元胞轉換為生態(tài)用地的概率，值域為0～1。輸入變量及其取值范圍如表1所示。

表1 ANN模型輸入變量及取值范圍Tab.1 Input variables and their ranges of ANN model

1.6 MCR-ANN-CA模型

以CA模型為基礎框架，耦合MCR模型與ANN模型，構建MCR-ANN-CA模型用于包頭市生態(tài)用地的演化模擬，技術路線如圖1所示。模型模擬的基本步驟為：

圖1 技術路線Fig.1 Technical roadmap

(1)利用MCR模型，綜合考慮多種土地適宜性因子，構建研究區(qū)范圍內(nèi)各像元演化為生態(tài)用地的累積阻力面。

(2)利用歷史數(shù)據(jù)的轉換情況及鄰域特征，構建訓練數(shù)據(jù)，訓練ANN模型，提取CA模型鄰域規(guī)則。

(3)以當前土地利用為CA模型的輸入數(shù)據(jù)，針對任一元胞，利用訓練好的ANN模型，根據(jù)其鄰域狀態(tài)判斷其演化為生態(tài)用地的概率P。

圖2 阻力因子的阻力分布Fig.2 Spatial distributions of resistance values of resistance factors

(4)將P與累積阻力面中像元最大值相乘，其結果為S(取值范圍為0到累積阻力面最大值)，若S大于當前元胞所對應的累積阻力，則該元胞演化為生態(tài)用地，否則，該元胞演化為非生態(tài)用地。

(5)將模擬結果轉換為柵格圖像并輸出結果。

2 結果與分析

2.1 基于MCR模型的適宜性規(guī)則構建

對土地利用、NDVI、坡度、政府規(guī)劃工業(yè)園區(qū)、人口密度、水體距離、高程等多種因子進行空間化與標準化處理，結果如圖2所示。其中將NDVI、坡度、人口密度、高程因子數(shù)據(jù)原始值進行歸一化處理作為其阻力。通過對不同用地類型阻力賦值，并對賦值結果歸一化處理得到土地利用數(shù)據(jù)。政府規(guī)劃工業(yè)園區(qū)阻力賦值為1，其他區(qū)域賦值為0，將數(shù)據(jù)柵格化后作為一個阻力因子。

根據(jù)包頭市2016年土地利用數(shù)據(jù)，提取包頭市生態(tài)用地范圍，如圖3a所示，由提取結果可知，包頭市生態(tài)源地主要分布在北部草原帶、大青山一帶及黃河沿岸，其他區(qū)域也有零星生態(tài)用地分布，大致形成了三屏多點的格局。

利用GIS疊加分析，將圖2中各因子進行空間疊加，得到包頭市各個空間位置演化為生態(tài)用地的阻力面，結果如圖3b所示。基于該阻力面，采用MCR模型，對源地進行累積耗費阻力計算，量化由近及遠的空間范圍內(nèi)各像元演化為生態(tài)用地面臨的阻力，得到如圖3c所示的累積生態(tài)阻力面。由圖3可知，包頭市南北部生態(tài)阻力較低，中部累積生態(tài)阻力較高，尤其圖3c中A區(qū)域及B區(qū)域，累積生態(tài)阻力達到包頭市最高水平。

圖3 包頭市生態(tài)用地演化累積生態(tài)阻力面構建Fig.3 Construction of ecological land transition minimum accumulative resistance surface of Baotou City

2.2 基于ANN模型的CA鄰域規(guī)則提取

基于研究區(qū)2006、2011年景觀格局分布數(shù)據(jù)，通過Matlab隨機抽取3 000個景觀像元，統(tǒng)計2006年景觀格局分布數(shù)據(jù)中，以3 000個像元為中心的3×3鄰域內(nèi)的各輸入變量大小，形成ANN模型訓練數(shù)據(jù)庫。選其中70%進行模型訓練，15%作為驗證數(shù)據(jù)集，15%作為測試數(shù)據(jù)集。

采用量化共軛梯度法對ANN進行訓練，模型經(jīng)57次迭代后均方根誤差達到最小，為0.12，如圖4所示，圖中圓圈處為最優(yōu)擬合點，此時驗證集的決定系數(shù)R2為0.90，模型擬合精度較好地滿足了本研究要求。

圖5 MCR-ANN-CA與CA-Markov模擬結果對比Fig.5 Comparison of simulation results of MCR-ANN-CA model and CA-Markov model

圖4 ANN模型均方根誤差Fig.4 Root mean square error of ANN model

2.3 MCR-ANN-CA生態(tài)用地模擬

將2011年景觀格局數(shù)據(jù)輸入訓練后的ANN模型，判斷各像元在鄰域規(guī)則影響下是否會演變?yōu)樯鷳B(tài)用地，進而根據(jù)MCR模型生成的累積耗費阻力面，判斷其最終演變方向,將模型模擬結果輸出為柵格圖像。包頭市2016年實際生態(tài)用地8 409.32 km2，MCR-ANN-CA模型的模擬結果中包頭市生態(tài)用地面積為8 670.01 km2，如圖5a、5b所示。

2.4 模型模擬精度對比分析

為驗證模型模擬精度，利用CA-Markov模型對包頭市2016年生態(tài)用地進行模擬，結果如圖5c所示。以2016年生態(tài)用地實際分布為參照，利用Idrisi Selva軟件中的Cross Tab模塊對2個模型的模擬結果進行分析，定量分析模型的模擬精度，結果如表2所示。表中KIA指數(shù)為卡帕一致性指數(shù)。

表2 模型模擬精度評價Tab.2 Accuracy assessment of simulation result

由表2可知，MCR-ANN-CA模型與CA-Markov模型模擬結果中，生態(tài)用地的面積均大于實際值，但總體上與2016年生態(tài)用地的實際分布保持了較高的一致性，其中MCR-ANN-CA模型的模擬精度略高于CA-Markov模型，二者的相對誤差分別為3.10%與5.31%，KIA指數(shù)分別為0.89和0.87。

相比CA-Markov模型，MCR-ANN-CA模型通過ANN模型提取了元胞自動機鄰域內(nèi)的轉換規(guī)則，同時利用MCR模型構建累積耗費阻力面，對不同空間位置元胞演化為生態(tài)用地的阻力進行了量化，因此模擬精度得到了進一步提高。

3 結論

(1)利用ANN模型提取了元胞自動機的鄰域規(guī)則，同時利用MCR模型構建累積耗費阻力面，基于MCR-ANN-CA模型對包頭市生態(tài)用地演化情況進行模擬，模擬精度較高。

(2)將MCR-ANN-CA模型與CA-Markov模型模擬結果進行對比，KIA指數(shù)分別為0.89和0.87，相對誤差分別為3.10%和5.31%，MCR-ANN-CA模型對包頭市生態(tài)用地的演化過程具有更高的模擬精度。