基于多模型的黑河中游甘臨高地區土地利用情景模擬

蔣小芳 段翰晨 廖 杰 宋 翔 薛 嫻

(1.中國科學院西北生態環境資源研究院沙漠與沙漠化重點實驗室， 蘭州 730000； 2.中國科學院西北生態環境資源研究院干旱區鹽漬化研究站， 蘭州 730000)

0 引言

土地利用變化關系到生態系統的穩定和持續發展，是土地利用/覆被變化(Land use and cover change，LUCC)研究的熱點問題[1]。土地利用模型為這一研究提供了強有力的技術支撐。CA模型出現于20世紀40年代，TOBLER[2]于1979年首次將CA模型應用于地理建模。CA模型涉及人類系統與城市土地系統之間的相互作用，具有同質性、空間離散性和時空局部性等特點，空間模擬能力優越[3-5]。然而傳統CA模型中每個像元內的土地利用類型單一，而實際的土地利用類型分布情況復雜。MCCA模型利用的混合元胞包括了不同土地利用類型，符合現實中復雜的土地利用結構，在元胞空間、元胞狀態和轉化規則的元胞鄰域方面均與傳統的CA模型存在差異，實現了從定性、靜態模擬到定量、動態模擬的跨越[6]。

然而單一模型難以兼具數量模擬和空間格局預測的功能[7]。隨著研究的深入，綜合數量模型與空間模型的耦合模型出現[8-14]。目前已有的數量模擬模型主要有馬爾可夫鏈(Markov chain，MC)、灰色預測模型(Grey model，GM)和MOP模型等[5,15-20]。MOP模型可用于計算土地系統中多維目標函數，具有強大的計算功能，是一種自上而下的方法，但在空間優化方面存在缺陷[21]。CA和MCCA模型均為自下而上的空間格局模擬方法，能夠生成優化的土地利用結構[22]。因此有必要同時采用自頂向下和自底向上的耦合模型來優化土地利用數量-空間格局[23]。

部分學者采用包括小尺度土地利用變化及其空間效益模型(Conversion of land use and its effects at small regional extent，CLUE-S)在內的純凈元胞研究模型對黑河中游甘州區或張掖市等區域的土地利用結構進行模擬[24-25]，但是缺乏對混合元胞的探究。鑒于此，本文以黑河中游甘臨高地區為例，確定MCCA模型在研究區的適用性后，分別將MOP和普通線性回歸這兩個數量預測模型與MCCA空間模擬模型耦合，模擬得到2035年可持續發展(Sustainable development，SUD)情景和基本發展(Basic development，BAD)情景下土地利用空間結構，采用效益定量評價方法對比分析，優選最佳的土地配置方案。

1 研究區和研究方法

1.1 研究區域概況

黑河流域地處中國西北干旱氣候區。本文研究范圍為黑河流域中游的張掖市甘州區、臨澤縣和高臺縣(甘臨高地區，圖1)。該區域年均氣溫約為8.2℃，年均降水量在80～150 mm之間，生態環境較為脆弱。土壤類型以水稻土、風沙土、草甸土、石質土和棕鈣土為主。土地利用類型較多，主要包括耕地、草地、建設用地和荒漠等。黑河流域中游經濟發展速度較快，人口較為稠密，人地矛盾突出。研究該區土地利用變化對于西北地區的經濟發展和生態環境保護具有重要意義。

圖1 研究區簡圖Fig.1 Chart of study area

1.2 土地利用研究模型

1.2.1MCCA模型

(1)MCCA模型原理

MCCA模型主要包括3部分[6](圖2)。第1部分是混合元胞內多維土地利用變量連續變化的機理挖掘。MCCA模型基于起始年份和目標年份兩期土地利用數據得出不同地類的連續變化，采用隨機森林回歸(Random forest regression，RFR)算法挖掘不同地類與驅動力因素之間的映射關系。RFR有利于克服不同空間變量的多重共線性問題。RFR首先訓練各地類的變化量(即增加量)與驅動因素之間的關系，因變量為每個地類的變化量，自變量為各驅動因素；然后預測得到各地類中每個元胞的變化潛力。與以前的模型不同，RFR并未同時擬合多個地類，而是分別訓練每種地類的變化潛力；而某一地類的增加表示其他地類的減少，因此模型只關注增加的地類，減少表示不增加，避免重復計算。

圖2 技術流程圖Fig.2 Flow chart of this research

第2部分是耦合亞元胞尺度土地競爭與定量轉化機制的CA模型。傳統的CA模型基于純凈元胞進行模擬，每個元胞代表一個地類，這脫離了實際情況，現實中每個元胞內部含有多個地類；而MCCA模型基于混合元胞內部每種地類的占比進行模擬，即亞元胞。該模型通過驅動系數調節不同地類的變化方向，首先需要計算土地利用成分的總變化概率，計算公式為

(1)

式中t——迭代次數i——混合元胞

k——地類組分O——總變化率

P——變化潛力Ω——領域效應

D——需求反饋

第3部分是MCCA模型模擬結果評價體系，主要包括亞像元混淆矩陣(Sub-pixel confusion matrix，SCM)、平均相對熵(Relative entropy，RE)和混合元胞質量系數(Mixed-cell figure of merit，mcFoM)，分別表征模擬結果的整體精度、土地利用結構相似度和混合元胞的變化精度。SCM越低說明模擬結果的整體精度越高，RE越接近0說明模擬過程中信息丟失越少，mcFoM越接近1說明模擬結果越準確。除此之外還有總體分類精度(Overall accuracy，OA)、生產者精度(Producer’s accuracy，PA)和用戶精度(User’s accuracy，UA)，均是越接近1模擬精度越高。

MCCA模型采樣率設置為10%，回歸樹數目為100，領域為3；轉換矩陣中0表示兩種土地利用類型之間無法轉換，1表示可以轉換，本研究在模擬過程中未設置限制轉換區域；步長表示土地利用類型的轉化率，經多次測試發現步長均為1時各項精度預測指標值最大。

(2)MCCA模型數據源及預處理

2000年和2015年的土地利用數據來源于中國科學院地理科學與資源研究所的資源環境科學與數據中心(http:∥www.resdc.cn)，為純像元數據。該數據的每個像元表征一種地類，空間分辨率為30 m(表1)。通過ArcGIS 10.2軟件的分區統計制表工具將這一數據聚合成空間分辨率為250 m的混合像元數據，并計算得到每個像元內部耕地、林地、草地、水域、建設用地和未利用地的覆蓋比例數據。

表1 MCCA模型數據源信息Tab.1 Research data information of MCCA model

其他表征土地利用變化的驅動力數據主要有自然環境數據(高程、坡度、坡向、河流矢量數據、2015年的年均氣溫及年均降水量數據)和社會經濟驅動力數據(道路矢量數據、居民點矢量數據、2015年的人口和國內生產總值(Gross domestic product，GDP))，與土地利用數據的來源相同。ArcGIS 10.2軟件中的歐氏距離計算工具可處理得到每個像元與道路、河流、居民點的距離柵格數據，該數據的空間分辨率為30 m。MCCA模型提供類似于ArcGIS軟件重采樣的坐標對齊機制，不要求驅動力柵格數據具有統一的像元分辨率(見MCCA模型操作手冊，下載網址：https: ∥github.com/HPSCIL/Mixed_Cell_Cellullar_Automata)。考慮到大部分驅動力數據的空間分辨率為30 m，為提高MCCA模型運行效率，通過ArcGIS 10.2軟件中的重采樣工具將驅動力數據的空間分辨率統一為30 m。

1.2.2普通線性回歸模型

MCCA模型主要用于土地利用空間結構可視化模擬，該模型開發者研發的MCCA軟件同時具有普通線性回歸功能(圖2中簡稱為MCCA的線性預測模塊)，該功能基于一次函數進行預測[6]。除MCCA模型空間模擬需要的2000年和2015年的土地利用數據外，本文引入2005年和2010年各土地利用類型的面積數據對未來各地類面積進行線性回歸預測，其來源和處理方法與2000年和2015年土地利用數據一致。

1.2.3MOP模型

MOP模型包括目標函數和約束條件兩部分，其研究結果對于認識未來的土地利用結構和當前土地規劃具有參考價值[26]。研究區面積為1 073 078 hm2，本研究利用MOP算法，基于Lingo 12.0軟件，對2035年各地類面積進行了優化。

(1) 目標函數

土地利用結構的優化應該遵循動態性、統籌兼顧和生態優先原則。土地系統內部的各個要素處在動態變化過程中。土地利用結構優化時應基于系統論觀點，兼顧整體和局部、協調人地矛盾，保證人與自然和諧發展。

生態效益目標函數的確定能提高研究區的生態價值，增強生態服務功能，維持生態系統的穩定性。當前生態效益主要采用價值量估算，并用貨幣值量化。本研究采用謝高地等[27]的中國生態系統服務價值當量因子表來估算生態效益，采用黑河中游張掖地區的糧食單位面積產量與全國糧食單位面積產量的比值作為修訂系數。參考《中國統計年鑒》和《張掖年鑒》[28-29]，2015年全國糧食單位面積的產量為5 482.900 kg/hm2，全國糧食平均價格為2.370元/kg，張掖市同年的糧食單產為7 198.260 kg/hm2，計算得到修訂系數為1.313。單位面積耕地生態系統服務價值當量因子的經濟價值為本年度全國糧食單產價值的1/7，其計算結果為1 856元/hm2。校正后耕地、林地、草地、水域和未利用地的生態效益分別為16 840、53 250、17 650、112 040、1 020元/hm2，建設用地的生態效益較低，本研究按零值計算。

黑河中游甘臨高地區土地利用結構優化是為獲得最大經濟效益，經濟效益目標的實現有利于提升工農業發展水平。本文采用每種地類單位面積的國內生產總值作為經濟效益目標函數的系數。參考2015年相關統計年鑒的經濟指標數據[28-31]，耕地、林地、草地、水域和建設用地的經濟效益核算結果分別為30 400、400、1 700、200、379 000元/hm2。

社會效益主要包括糧食安全保障和農民生存保障，參考前人的研究，耕地、林地、草地、水域、建設用地和未利用地的社會價值分別設置為23 360、2 520、2 340、5 160、99 040、10元/hm2[32]。

(2) 約束條件

約束條件包括3部分：土地總量約束、社會經濟約束和生態環境約束。研究區各種土地利用類型的約束條件主要參考2015年實際土地類型面積比例、各土地面積變化趨勢以及普通線性回歸模型得到的2035年土地類型面積比例。社會經濟約束包括經濟總量約束、土地開發強度約束和土地利用率約束。生態環境約束包括生態承載力和生態綠當量兩方面。

科研數據管理過程具有明顯的周期性和階段性特征，基于科研活動不同階段的數據形態和數據處理活動，形成關于科研數據管理的相關生命周期理論。重點介紹了ICPSR社會科學數據存檔生命周期管理模型。

生態承載力是產量因子、均衡因子與各土地利用類型面積的乘積，用來衡量生態系統服務能力的上限[33-34]。參考相關研究，耕地、林地、草地、水域和建設用地的均衡因子分別為2.210、1.340、0.490、0.360和2.210[35]，產量因子分別為1.310、0.910、0.390、0.980和1.660[36-37]，研究區的生態承載力應高于2015年(表2)。

生態綠當量主要用于衡量某一區域生態環境的優劣。因地制宜確定本研究區的耕地、林地、草地均為擁有綠當量的土地。生態綠當量計算公式為

(2)

式中F——生態服務功能價值系數

i等于1為耕地，i等于2為林地，i等于3為草地。

2 結果與分析

2.1 MCCA模型適用性評價

2.1.1模型精度評價

采用OA、Kappa系數、mcFoM、RE、PA和UA等參數評估MCCA模型的預測精度。OA、Kappa系數和mcFoM的計算結果分別為0.934、0.886和0.261。所有地類的PA和UA均高于0.700，且除林地和建設用地外其他土地利用類型的PA和UA均高于0.800。RE越低，模擬過程中損失的信息量越少，說明模型模擬結果與實際結果更相符。模擬結果和實際土地利用數據之間的RE均值為0.508，表明模擬過程中信息損失較少，二者之間的相似度高。各項評估指標值均較理想，這說明模型可用于預測黑河中游甘臨高地區未來的土地利用變化情況。

采用MCCA模型分別計算耕地、林地、草地、水域、建設用地和未利用地的2015年模擬結果和實際土地利用分布數據較2000年的變化情況，變化程度為-1～1，值為正表示土地利用類型的面積在擴張，值為負則表示減少。可以發現所有土地利用類型的變化區域基本吻合，變化地類的位置擬合效果較好，變化的斑塊數量差異不明顯。各土地利用成分的模擬結果與真實情況的變化規律基本一致，僅有少數像元在變化程度上存在差異(圖3，圖中左圖為實際的土地利用成分的變化程度，右圖為模擬的土地利用成分的變化程度)。

圖3 2000—2015年實際和模擬的土地利用成分變化程度Fig.3 Actual and simulated land use changes during 2000—2015

為更深入探究基于混合元胞的MCCA模型的優越性，使用PLUS模型對空間分辨率250 m的純像元土地利用結構數據進行模擬，模擬結果的FoM、Kappa系數、OA分別為0.117、0.853、0.916，而MCCA模型的OA、Kappa系數和mcFoM的計算結果分別為0.934、0.886和0.261。據此可以確認MCCA模型的模擬結果與地表實際狀況吻合度更高，能夠很好地模擬復雜土地利用結構。

2.1.3MCCA模型在研究區的最適模擬尺度

不同模型對于柵格尺度的要求存在差異。LIANG等[6]采用4種空間分辨率(250、500、750、1 250 m)的混合像元驗證MCCA模型的模擬精度，發現250 m的模擬精度最高。鑒于此本文將混合像元空間分辨率縮小至250～500 m之間進一步確定最佳的混合像元大小，結果表明在250、300、350、400、450、500 m的混合像元中，依然是250 m混合像元模擬精度最高。不同柵格尺度的模擬準確度評價指標值見圖4，6種柵格尺度模擬結果的OA均高于0.920，RE均值均低于0.550，Kappa系數均高于0.870，mcFoM均高于0.210。隨著混合像元的增大，RE均值波動上升，說明柵格尺度越大模擬過程中丟失的柵格信息越多；OA、Kappa系數和mcFoM則波動降低，說明MCCA模型的擬合準確度隨著柵格尺度的增大而趨于降低。

圖4 不同柵格尺度的模擬精度Fig.4 Simulation accuracy of different raster scales

2.2 2035年不同情景預測結果

預知未來的土地利用情況對維持生態系統的穩定性有重要意義，本文基于普通線性回歸模型的BAD情景中2015—2035年期間耕地、水域和建設用地逐年遞增(圖5)。說明隨著城市化進程的加快，二三產業的投資明顯增多，建設用地面積迅速擴張。尤其是人口老齡化導致二孩政策開放，人口數量的增加對糧食的需求量加大，耕地面積也在逐年上升。BAD情景中草地、林地和未利用地類型逐年遞減，表明生態環境在逐漸惡化。

圖5 研究區2000—2015年實際的和2015—2035年預測的土地利用類型覆蓋度Fig.5 Actual (2000—2015) and predicted (2015—2035) land use proportions in study area

BAD和SUD情景的耕地、水域和建設用地均較2015年有所增長，而草地則低于2015年。SUD情景中2035年的林地面積占比明顯高于BAD情景，草地略高于BAD情景，而耕地、建設用地和未利用地則略低于BAD情景。SUD情景占用的未利用土地較少，說明SUD情景的經濟發展狀況有利于保護脆弱的生態環境(表3)。

表3 不同情景下的2035年土地利用面積占比預測結果Tab.3 Land use percentage under different scenarios in 2035

MCCA模型對2035年地類結構進行空間可視化的結果顯示(圖6，水域變化程度較小，故未列出)，BAD情景中研究區的耕地不斷向外擴張，侵占林地和草地；SUD情景的耕地和草地變化趨勢與BAD情景一致，但是占用林地的現象較少。IPCC報告指出2016—2035年全球氣溫和降水會逐年增加[38]，研究區的水域范圍的增加驗證了這一結論。隨著城市化進程的加快，兩種情景中城鎮用地和農村用地均在原有基礎向外擴張，甘州區的市區范圍擴展尤為明顯；高臺縣和臨澤縣的城市發展依托于交通網絡，未來在城市周邊的交通集散地將形成較大的衛星城。兩種情景中的建設用地和耕地在未來20年將向未利用地延伸，現有的沙漠和耕地交界處在未來有極大可能轉變為農耕區，這可能是經濟技術發展以及政策引導的結果。但是并不是所有的建設用地均會向外擴展，MCCA模型的模擬結果可能與現實情況存在少數不一致之處，比如臨澤北部的工礦用地未來可能會萎縮，而模擬的結果卻顯示該區域將會擴張，這可能是因為MCCA模型模擬過程中受到驅動因素時間滯后性的影響。

圖6 2015年實際土地利用結構、2035年模擬土地結構和2015—2035年土地利用成分變化程度Fig.6 Actual land use structure in 2015, simulated land use structure in 2035, and change of land use percentage during 2015—2035

2.3 2015年與2035年的土地資源價值核算

通過量化2035年BAD情景下的生態效益、經濟效益、社會效益(表4)可以發現，2035年綜合效益達3 504 943萬元，而2015年僅為2 930 312萬元。BAD情景下耕地、水域和建設用地的3種效益均增加，且耕地及建設用地的效益增速較水域快。西部大開發戰略的實施帶動經濟發展，企業投資增加，建設用地隨之擴張。國家人口政策的實施帶動了未來總人數的上漲，糧食需求量大增，農用地面積增加。BAD情景雖然經濟效益快速增加，但是生態效益增速較慢，社會效益增速居中。BAD情景下林地生態效益、經濟效益和社會效益均明顯降低。

表4 研究區不同情景效益對比Tab.4 Benefit comparison of different scenario in study area 萬元

2015年研究區建設用地的經濟效益、社會效益分別為854 144、223 213萬元，而2035年SUD情景下則分別增至1 077 746、281 647萬元。黑河中游甘臨高地區的經濟總量持續增長，建設用地的范圍大幅擴張。2035年耕地面積得到了優化，耕地的各項效益高于2015年。但是耕地并未盲目擴張，占用林地的現象較少，增加的耕地主要是未利用地和草地演變而來，這有利于保障資源的可持續利用。2035年SUD情景下草地和未利用地的效益略低于2015年，而林地與2015年持平。這說明政府在發展經濟的同時注重保護環境。

2035年SUD情景和BAD情景中林地的生態效益分別為36 057、22 752萬元，SUD情景的生態效益較BAD情景增速更快。SUD情景較BAD情景更好地協調了生態、經濟和社會三者之間的關系。SUD情景建設用地和耕地各效益的增速低于BAD情景，而林地的三大效益高于BAD情景，工農業用地的適度增加有助于維持社會穩定、保障糧食安全和保護生態環境。尤其是西北干旱區生態脆弱，環境問題突出，保護環境勢在必行。

3 討論

3.1 基于混合元胞的MCCA模型模擬精度和時間跨度分析

MCCA模型能夠自下而上生成優化的土地利用空間格局，改進了傳統的純像元CA模型，混合元胞的轉化規則是定量的，其鄰域和元胞空間的狀態具有混合、連續和多維的特點，與傳統元胞土地利用成分離散、一維和定性相反，是從定性的靜態模擬到定量的動態模擬的跨越[6]。因此，MCCA模型較基于純凈元胞的PLUS模型的模擬準確度更高，本文的研究結果證實了這一點。而LIANG等[6]將MCCA模型與基于純凈像元的未來土地利用(Future land use simulation model，FLUS)模型進行對比分析同樣得出基于混合元胞的MCCA模型模擬準確度更高。

受到現實情況的制約，MCCA模型未統一兩期現有數據的時間跨度與現有數據-預測年份數據的時間跨度，PLUS、FLUS和CLUE-S等空間模擬模型也存在類似問題。這類模型通常是基于土地利用柵格數據以及驅動力柵格數據采用回歸算法(如隨機森林回歸、人工神經網絡和邏輯斯蒂回歸等算法)挖掘不同地類與驅動力因素之間的關系[6,39-40]；然后基于現有的土地利用數據中各地類的面積，采用數量預測模型(如馬爾科夫鏈、灰色預測模型、多目標線性優化模型和系統動力學模型等)預測未來指定年份的地類面積，因此時間跨度的不統一在地類面積數值預測方面存在合理性[41]；最后基于基期土地利用柵格數據、地類變化可能性和未來指定年份的地類面積進行預測，獲得未來的土地利用柵格數據。當研究中可獲得柵格數據之間時間跨度過短時，考慮到社會發展的慣性及預測年份數據的參考作用，本研究未統一時間跨度。

3.2 不同情境效益分析和MOP模型優勢

多目標空間優化問題較為復雜，需要同時兼顧多方面因素[42]。MOP模型采用自上而下的方式解決復雜土地系統問題，具備強大的數據處理能力，能夠在綜合考慮生態效益、經濟效益和社會效益的基礎上，從多個方面優化土地利用結構[43]。而普通線性回歸模型僅從各地類面積的數量變化方面進行分析，難以實現土地優化配置的目的。將基于普通線性回歸的BAD情景與基于MOP算法的SUD情景進行對比，計算不同情景下各土地利用類型的生態效益、經濟效益和社會效益，發現2035年土地優化配置的SUD情景綜合效益遠高于2015年，雖然其經濟效益和社會效益略低于BAD情景，但生態效益明顯高于2015年，這有利于緩解當地脆弱的生態環境。而將MOP和MCCA進行耦合模擬2035年研究區的土地利用空間結構，能夠充分利用數量模擬模型和空間模擬模型的優勢，實現未來土地利用結構的可視化，為相關部門進行土地規劃和生態修復提供參考。

4 結論

(1)采用Kappa系數、OA、PA、UA、RE和mcFoM等指標進行精度驗證，各項指標的計算結果均說明MCCA模型的模擬準確度高，且MCCA模型的模擬精度明顯優于基于純凈元胞的PLUS模型，這驗證了MCCA模型在黑河中游甘臨高地區具有高度適用性。為進一步確定最佳像元尺度，選擇分辨率250、300、350、400、450、500 m的混合元胞數據進行模擬，結果表明分辨率為250 m時模型模擬精度最高，Kappa系數、mcFoM和平均RE分別為0.886、0.261和0.508，因此采用250 m的混合像元數據預測2035年的土地利用結構。

(2)在基于普通線性回歸模型的BAD情景中，2000—2035年期間的耕地、水域和建設用地范圍逐年增加，林地、草地和未利用地則相反，存在工農業用地占用林草地的現象。SUD情景的建設用地、草地和未利用地的變化趨勢與BAD情景基本一致。但SUD情景耕地的增加速率較BAD情景偏低，其2035年林地的范圍大于BAD情景，該情景中工農業用地占用林草地的情況較少。SUD情景的經濟效益和社會效益略低于BAD情景，而生態效益明顯高于BAD情景。這說明基于多目標線性規劃(MOP)模型與MCCA模型的SUD情景能夠實現土地優化配置，在發展工農業的同時保護生態環境。