樊文杰, 戴曉愛*, 謝一茹, 高怡凡

(1.成都理工大學(xué)地球科學(xué)學(xué)院， 成都 610059； 2.北京師范大學(xué)環(huán)境演變與自然災(zāi)害教育部重點實驗室， 北京 100875)

土地是人類賴以生存和發(fā)展的物質(zhì)基礎(chǔ)，是人類存在、生活和生產(chǎn)的源泉[1-2]。土地利用是人類根據(jù)其自身發(fā)展需求，對土地加以改造開發(fā)利用的過程[3]。隨著人口數(shù)量的不斷增加以及社會的進(jìn)步，人地矛盾日益尖銳，因此對土地利用進(jìn)行科學(xué)的模擬和規(guī)劃至關(guān)重要。土地利用變化模型通過分析以往的土地利用數(shù)據(jù)來理解土地利用系統(tǒng)的機(jī)理[4]，并評估土地政策對自然生態(tài)和經(jīng)濟(jì)發(fā)展的影響，為土地利用提供決策支持，以達(dá)到人類社會的可持續(xù)發(fā)展。

土地利用變化模型層出不窮，其中使用較為廣泛的模型包括：元胞自動機(jī)模型、馬爾科夫鏈模型、灰色系統(tǒng)預(yù)測模型、CLUE-S(conversion of land use and its effects at small regional extent)模型等。中國學(xué)者已經(jīng)應(yīng)用不同的土地利用變化模型對不同尺度下的研究區(qū)土地利用做了豐富的研究[5-12]。其中， CLUE-S模型使用廣泛，得到了非常好的驗證。此外，與其他模型相比，CLUE-S模型能夠?qū)Σ煌恋乩梅绞竭M(jìn)行同步模擬。

素有天府之國美譽的四川省，其GDP總量在中國位列前茅，生態(tài)系統(tǒng)服務(wù)多樣[13]。然而四川省作為人口大省，伴隨著近年來人口的持續(xù)增加，該地區(qū)的生態(tài)保護(hù)壓力、耕地安全壓力和經(jīng)濟(jì)發(fā)展壓力日益增加，經(jīng)濟(jì)發(fā)展與糧食安全、生態(tài)保護(hù)方面的協(xié)調(diào)發(fā)展至關(guān)重要。土地利用變化模擬能夠在政策目標(biāo)約束下實現(xiàn)對四川省未來土地利用變化空間格局的預(yù)測，預(yù)測結(jié)果對當(dāng)前四川省面臨的問題將具有重要的參考意義。當(dāng)前關(guān)于四川省的土地利用變化模擬研究較為豐富，但大部分研究區(qū)域僅針對市、縣一級[14-16]，未針對四川省全域進(jìn)行探索，且模擬情景單一，多為對現(xiàn)狀的分析[17-18]。

因此，現(xiàn)使用土地利用變化模型CLUE-S，對四川省全域范圍進(jìn)行土地利用時空格局動態(tài)演變模擬，探究未來數(shù)年的土地利用變化格局，探尋緩和人地矛盾、調(diào)控區(qū)域生態(tài)安全的方法，以期為四川省實現(xiàn)高質(zhì)量、可持續(xù)發(fā)展提供科學(xué)參考。

1 研究區(qū)域

1.1 研究區(qū)域概況

四川省(26°03′～34°19′N、97°21′～108°12′E)位于中國西南部，轄21個市(州)，總?cè)丝诩s8 375萬人，總面積48.67萬km2；位于青藏高原地區(qū)向長江中下流平原的過渡地帶，整體地勢呈西高東低。川西地區(qū)為高原山地氣候，以高原和山地為主，森林草地分布范圍廣，是主要的生態(tài)系統(tǒng)服務(wù)供給區(qū)域；川東地區(qū)為亞熱帶季風(fēng)氣候，以丘陵和盆地為主，城市化水平高、經(jīng)濟(jì)發(fā)達(dá)，是經(jīng)濟(jì)開發(fā)重點區(qū)域。四川省2015年各土地利用類型面積從大到小依次為：草地(17.53萬km2)、林地(16.56萬km2)、耕地(11.93萬km2)、未利用地(1.89萬km2)、建筑用地(0.47萬km2)、水域(0.29萬km2)，空間分布如圖1所示。

圖1 四川省2015年土地利用分布Fig.1 Land use distribution in Sichuan Province in 2015

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

本文中使用了中國科學(xué)院地理科學(xué)與資源研究所(http://www.resdc.cn)的2015年中國省級行政邊界數(shù)據(jù)、2015年中國GDP空間分布公里網(wǎng)格數(shù)據(jù)集、全國1 km的數(shù)字高程模型和4期(2000年、2005年、2010年、2015年)土地利用數(shù)據(jù)；國家地球系統(tǒng)科學(xué)數(shù)據(jù)中心(http://www.geodata.cn)的2015年中國1 km分辨率平均氣溫數(shù)據(jù)集和平均降水量數(shù)據(jù)集；其他空間數(shù)據(jù)主要包括：Global Human Settlement(https://ghsl.jrc.ec.europa.eu/ghs_pop2019.php)的全球人口空間分布數(shù)據(jù)；Weiss等[19]提供的城市可達(dá)性數(shù)據(jù)；OpenStreetMap(http://download.geofabrik.de)的河流、湖泊和公路網(wǎng)數(shù)據(jù)；國際土壤參考資料中心(https://data.isric.org)的土壤有機(jī)碳密度和土壤質(zhì)地數(shù)據(jù)。社會經(jīng)濟(jì)數(shù)據(jù)來自中國國家統(tǒng)計局和四川省十四五規(guī)劃相關(guān)文件。研究中使用到的坡度數(shù)據(jù)是使用ArcGIS中Slope工具將高程數(shù)據(jù)作為輸入生成的；可達(dá)性柵格數(shù)據(jù)是使用ArcGIS中Euclidean Distance將對應(yīng)的矢量文件作為輸入得到的。所有使用到的柵格數(shù)據(jù)經(jīng)重采樣處理后統(tǒng)一成1 km×1 km的分辨率。

2 研究方法

2.1 模型

研究中采用由荷蘭瓦赫寧根大學(xué)的Verburg等[20]開發(fā)的CLUE-S模型。該模型的前身是CLUE模型，相比于CLUE模型而言更適合省級區(qū)域范圍的研究。CLUE-S模型是以空間位置適宜性分析為基礎(chǔ)，結(jié)合土地系統(tǒng)時空動態(tài)競爭和相互作用，來模擬土地利用變化。

CLUE-S模型可以分為需求量模塊和空間分配模塊[20-22]。需求量模塊根據(jù)需求情景在整體層面上計算不同土地利用類型的總需求量，該模塊通常是通線性規(guī)劃或者馬爾可夫過程等方法處理后作為CLUE-S模型的輸入?？臻g分配模塊將上述土地需求量反映為空間上某一位置土地利用類型的變化。CLUE-S模型需要土地政策及限制、土地利用類型轉(zhuǎn)換、土地利用類型需求量和土地位置屬性四類信息。土地政策及限制表明根據(jù)相關(guān)法律政策哪些區(qū)域的土地利用類型不允許轉(zhuǎn)換，比如永久基本農(nóng)田、國家級自然保護(hù)區(qū)、國家森林公園等。土地利用類型轉(zhuǎn)換有轉(zhuǎn)換規(guī)則和轉(zhuǎn)換彈性兩項，前者表明不同的各種土地利用類型能否相互轉(zhuǎn)換，后者反映土地利用類型維持自身狀態(tài)的能力。土地利用需求量是第一個模塊得到的結(jié)果。土地位置屬性信息用于將一系列相關(guān)因子作為輸入使用邏輯回歸計算出某個位置最佳的土地利用類型。

2.2 發(fā)展情景設(shè)計

基于研究區(qū)域自身特點和其發(fā)展戰(zhàn)略，綜合考慮生態(tài)效益、經(jīng)濟(jì)效益，設(shè)計了三種情景。第1種是生態(tài)效益型情景，該情景在滿足政府預(yù)期經(jīng)濟(jì)增長速度的前提下，能夠達(dá)到生態(tài)效益最大化。第2種是經(jīng)濟(jì)效益型情景，該情景以不減少生態(tài)效益為基礎(chǔ)，最大化經(jīng)濟(jì)效益。第3種是綜合效益性情景，該情景為情景1和情景2的折中實現(xiàn)，該情景首先要求生態(tài)效益增量為情景1下生態(tài)效益增加量的50%，然后最大化經(jīng)濟(jì)效益。

選取線性規(guī)劃來計算土地利用需求量。根據(jù)模擬情景，可以確定目標(biāo)函數(shù)[式(1)]、面積約束條件[式(2)～式(4)]、價值量約束條件[式(5)]為

X=argmaxX(KX)

(1)

s.t.Xi≥Xmin

(2)

Xj≤Xmax

(3)

(4)

(5)

式中：X為通過線性規(guī)劃求解得到的各土地利用類型面積；Xi表示第i種土地利用類型的所占面積,km2；Xmin是第i種土地利用類型面積的最小值,km2；Xmax是第i種土地利用類型面積的最大值,km2；Varea是四川省的總面積,km2。

在情景1中，式(1)中K表示各土地利用類型的生態(tài)效益系數(shù)，式(5)中r表示到模擬終止年份GDP的增長率，Ci表示各土地利用類型的經(jīng)濟(jì)效益系數(shù)。在情景2和情景3中，式(1)中K表示經(jīng)濟(jì)效益系數(shù)，Ci表示生態(tài)效益系數(shù)。在情景2中，式(5)中r=0；而在情景3中，式(5)中r的值等于情景1中生態(tài)效益增長率的50%。

生態(tài)效益系數(shù)法、經(jīng)濟(jì)效益系數(shù)法[11,23]用來度量土地的生態(tài)效益和經(jīng)濟(jì)效益。經(jīng)濟(jì)效益系數(shù)是根據(jù)國家統(tǒng)計局公布的研究區(qū)域農(nóng)業(yè)總產(chǎn)值、林業(yè)總產(chǎn)值、牧業(yè)總產(chǎn)值和漁業(yè)總產(chǎn)值以及研究區(qū)域全年GDP除以每種土地利用類型相應(yīng)的面積計算得到的。此外，四川省內(nèi)的未利用地主要是由裸巖石質(zhì)地的土地組成，所以將未利用地的經(jīng)濟(jì)效益系數(shù)設(shè)置為0。生態(tài)效益系數(shù)根據(jù)四川省生態(tài)系統(tǒng)服務(wù)價值量[24]得到。計算得到的系數(shù)見表1。

表1 土地利用類型的生態(tài)效益系數(shù)和經(jīng)濟(jì)效益系數(shù)Table 1 Weight of ecosystem services value and economic benefit coefficients by land use types

根據(jù)四川省“十四五”規(guī)劃文件，本文中設(shè)定在所有的模擬情景中耕地面積不再減少，草地的面積不少于當(dāng)前草地面積的95%，未利用地的面積不少于當(dāng)前未利用地面積的90%。水域的面積較短時期內(nèi)往往難以發(fā)生驟變，因此采用GM(1,1)模型預(yù)測其面積。根據(jù)上述方程使用MATLAB軟件中的線性規(guī)劃相關(guān)函數(shù)求解得到2030年各土地利用類型的面積預(yù)測值。由于CLUE-S模型需要輸入逐年的面積預(yù)測值，故依據(jù)2015年和2030年的數(shù)據(jù)使用線性插值法得到中間各年面積的預(yù)測值。

2.3 模型參數(shù)設(shè)定

2.3.1 Logistic回歸分析

在CLUE-S模型中土地利用轉(zhuǎn)換往往會發(fā)生在“適宜性”最佳的位置，這種適宜性依賴于研究區(qū)域的自然地理環(huán)境和社會經(jīng)濟(jì)發(fā)展?fàn)顩r。為了量化這種適宜性，根據(jù)研究區(qū)特性，選取了人均GDP、人口密度、氣溫、降水、高程、坡度等13種驅(qū)動因子，使用Logistic回歸分析得到某一柵格單元轉(zhuǎn)變?yōu)槟撤N土地利用類型適宜性值，即

(6)

式(6)中：Pi是當(dāng)柵格單元是第i種土地利用類型時的適宜性，X是驅(qū)動因子在該柵格單元的值。逐步回歸可以使得與第i種土地利用類型相關(guān)性較小的驅(qū)動因子從回歸方程種剔除出去，也就是將驅(qū)動因子對應(yīng)的β系數(shù)設(shè)置為0。由于ROC(receiver operation characteristic)方法具有不受樣本不均衡影響等一系列優(yōu)點，故使用ROC方法對回歸結(jié)果的好壞做出評價。ROC曲線下面積稱為AUC，AUC的值介于0.5～1。AUC的值越大，說明擬合結(jié)果越好，驅(qū)動因子對柵格單元土地利用類型的解釋性就越強(qiáng)。

由于要對每種土地利用類型進(jìn)行回歸分析，故把土地利用類型根據(jù)中國科學(xué)院土地利用分類系統(tǒng)將土地利用類型重分類成六大類[25]，再次使用重分類工具制作對應(yīng)的二值圖。CLUE-S模型沒有內(nèi)置Logistic回歸分析功能，所以先使用ArcGIS中Raster to ASCII工具把13種驅(qū)動因子對應(yīng)的柵格數(shù)據(jù)以及土地利用二值圖像轉(zhuǎn)換成ASCII格式，再使用模型自帶的FileConvert工具把ASCII格式的數(shù)據(jù)轉(zhuǎn)換成統(tǒng)計分析軟件SPSS接受的格式，最后運行SPSS中的二元Logistic回歸分析功能，得到的結(jié)果如表2所示。表2表明，所選驅(qū)動因子對耕地、草地和未利用地解釋性非常好，其AUC值分別為0.900 512、0.872 227和0.845 733；對建設(shè)用地和林地的解釋性較好，其AUC值分別為0.789 206和0.705 876。

表2 Logistic回歸的分析結(jié)果與檢驗結(jié)果Table 2 Results of logistic regression and ROC test

2.3.2 轉(zhuǎn)換規(guī)則和轉(zhuǎn)換彈性

通過使用2000年、2005年、2010年、2015年四期四川省土地利用柵格數(shù)據(jù)進(jìn)行統(tǒng)計分析，可以得到近些年來土地利用類型轉(zhuǎn)換情況，然后使用上述數(shù)據(jù)作為土地利用轉(zhuǎn)換矩陣，如表3所示。土地利用轉(zhuǎn)換矩陣是一個二維數(shù)字矩陣，矩陣中的取值有0和1兩種。如果第i行第j列的值為0，則表示的是不允許第i行所代表的土地利用類型向第j列代表的土地利用類型發(fā)生轉(zhuǎn)換，否則可以發(fā)生轉(zhuǎn)換。

表3 土地利用轉(zhuǎn)換矩陣Table 3 Matrix of land use conversion

同樣使用上述四期數(shù)據(jù)進(jìn)行統(tǒng)計分析，可以得到近些年來各土地利用類型向其他土地利用類型轉(zhuǎn)出的比率。比率越小，說明該土地利用類型維持類型不變的能力越強(qiáng)，彈性也就越大。結(jié)合轉(zhuǎn)出比率數(shù)據(jù)以及專家經(jīng)驗，得到轉(zhuǎn)出彈性系數(shù)設(shè)置結(jié)果，如表4所示。

表4 土地利用類型轉(zhuǎn)換彈性Table 4 Conversion elasticity of land use types

2.4 模型驗證

為了驗證CLUE-S模型運行結(jié)果的準(zhǔn)確性，在正式模擬之前，首先以四川省2010年土地利用圖為起點，模擬四川省2015年土地利用情況。將模擬結(jié)果與真實結(jié)果作對比。如果相同或高度相似，則說明模擬效果好，反之亦然。在評價相似性時，采用Kappa系數(shù)法。具體而言，Kappa系數(shù)法需要輸入兩幅影像，綜合考慮像元數(shù)量相似性以及位置相似性，來判斷這一對影像是否具有空間分布一致性。本文中共采用以下三種Kappa系數(shù)：標(biāo)準(zhǔn)Kappa、反映數(shù)量相似性的KHis、反映位置相似性的KLoc，這些系數(shù)的計算公式為

(7)

(8)

(9)

式中：po為觀測圖和模擬圖中像元值相同的像元個數(shù)與總像元數(shù)的比值，即總體精度；pe為服從觀測圖數(shù)據(jù)分布下的一致性期望概率；pmax為服從觀測圖數(shù)據(jù)分布下的一致性最大概率。標(biāo)準(zhǔn)Kappa的值等于KHis與KLoc的乘積。通常情況下Kappa值介于0～1，0表示觀測圖與模擬圖之間具有極低的一致性，1表示觀測圖與模擬圖之間具有高度一致性。

本文中使用了根據(jù)每一類別像元總數(shù)計算得到的KHis、根據(jù)每一類別像元空間分布計算得到的KLoc和標(biāo)準(zhǔn)Kappa。將上述兩幅影像作為 Map Comparison Kit 軟件的輸入，可以得到KHis系數(shù)的值為0.999 77、KLoc系數(shù)的值為0.896 76、Kappa系數(shù)的值為0.896 55，表明CLUE-S模型參數(shù)的設(shè)定能較好地反映和預(yù)測研究區(qū)域土地利用情況。

3 結(jié)果與分析

3.1 不同情景下土地需求結(jié)果分析

通過線性規(guī)劃計算出3種模擬情景下的土地利用需求面積，如表5所示。在3種情景下耕地的需求面積與2015年耕地面積基本相同。林地面積在生態(tài)效益情景下最大，經(jīng)濟(jì)效益情景下最小，綜合效益型居于兩者之間，這也體現(xiàn)出林地具有較高的生態(tài)效益，并且三種情景中林地面積與現(xiàn)狀值相比都略有增長。從表5種可以看出草地、水域和未利用地的面積在3種情景中分別相同，這也符合本文的需求設(shè)定。建設(shè)用地的面積在經(jīng)濟(jì)效益型情景中最大，在生態(tài)效益型情景中最小，綜合效益型居于兩種情景之間，這也體現(xiàn)出建設(shè)用地具有較高的經(jīng)濟(jì)效益。

3.2 不同情景下土地模擬結(jié)果分析

運行CLUE-S模型可以得到不同模擬情景下的四川省2030年不同土地利用類型面積的大小及其空間分布，如表6和圖2所示。通過比較表5和表6，發(fā)現(xiàn)兩者對應(yīng)數(shù)值的大小幾乎保持一致，說明模型結(jié)果能夠較為理想地反映設(shè)計的情景需求。

表5 不同模擬情景預(yù)測的土地利用需求面積Table 5 Area of land use demand predicted by different simulation scenarios

表6 不同模擬情景模擬的土地利用需求面積Table 6 Area of land use demand simulated by different simulation scenarios

圖2 不同情景下土地利用分布Fig.2 Land use distribution under different scenarios

根據(jù)圖2，從整體上看，耕地主要分布在成都平原、四川省的東北部地區(qū)、四川省的南部地區(qū)、涼山州和攀枝花市。林地分布在除了成都平原及其周邊地區(qū)外的地區(qū)。草地和未利用地集中分布在甘孜州和阿壩州。水域分布在四川省的東北部以及涼山州和攀枝花市。建設(shè)用地分布在成都周邊的城市、四川省的東北部以及攀西地區(qū)。這與當(dāng)前的土地利用分布現(xiàn)狀基本相同。

生態(tài)效益型情景模擬結(jié)果如圖2(a)所示。該情景下的生態(tài)效益為108 405.5億元，經(jīng)濟(jì)效益為71 913.6億元，與其他兩種情景相比其生態(tài)效益達(dá)到最大值。從地域分布來看：攀西地區(qū)的耕地面積減少，川東北地區(qū)北部耕地面積也減少，川西平原西南部地區(qū)的耕地面積增加；川西平原西南部地區(qū)的林地減少，攀西地區(qū)的林地會有所增加；四川盆地北部邊沿和西部邊沿的草地會減少，攀西地區(qū)西部的草地也會減少，川西北地區(qū)草地有所增加；四川盆地內(nèi)水域面積會減少，攀西地區(qū)和川東北地區(qū)水域面積會增加；新的建設(shè)用地主要集中于攀西地區(qū)和川東北地區(qū)；未利用地主要在川西北地區(qū)內(nèi)發(fā)生變遷。經(jīng)濟(jì)效益型情景模擬結(jié)果如圖2(b)所示。該情景下的生態(tài)效益為107 857.9億元，經(jīng)濟(jì)效益為78 346.3億元，與其他兩種情景相比其經(jīng)濟(jì)效益達(dá)到最大值。綜合效益性情景模擬結(jié)果如圖2-c所示。該情景下的生態(tài)效益為108 068.2億元，經(jīng)濟(jì)效益為75 480.7億元，其生態(tài)效益和經(jīng)濟(jì)效益都處于中間水平。在情景2和情景3中，各土地利用類型的動態(tài)變化與情景1基本保持一致。

三種情景中土地利用分布在成都平原經(jīng)濟(jì)區(qū)、川東北經(jīng)濟(jì)區(qū)和攀西經(jīng)濟(jì)區(qū)差異最為顯著。與情景2相比，情景1在成都平原經(jīng)濟(jì)區(qū)有更多的林地(+484 km2)，在川東北經(jīng)濟(jì)區(qū)和攀西經(jīng)濟(jì)區(qū)建設(shè)用地面積會減少(分別為-371 km2和-612 km2),體現(xiàn)出了該情景對于生態(tài)效益增加的需求；與情景3相比，情景2在成都平原經(jīng)濟(jì)區(qū)林地面積會減少(-220 km2)，在川東北經(jīng)濟(jì)區(qū)和攀西經(jīng)濟(jì)區(qū)有更多的建設(shè)用地(分別為+151 km2和+276 km2),體現(xiàn)出了該情景對于經(jīng)濟(jì)效益增加的需求；與情景1相比，情景3在成都平原經(jīng)濟(jì)區(qū)的林地也會減少(-264 km2)，在川東北經(jīng)濟(jì)區(qū)和攀西經(jīng)濟(jì)區(qū)建設(shè)用地面積會增加(分別為+220 km2和+336 km2)。

4 結(jié)論

根據(jù)研究區(qū)的發(fā)展設(shè)定了3種發(fā)展情景，使用CLUE-S模型對四川省2030年土地利用變化進(jìn)行了模擬研究。

(1)研究結(jié)果表明生態(tài)效益型情景能夠最大化生態(tài)系統(tǒng)服務(wù)價值，生態(tài)系統(tǒng)服務(wù)價值總量為108 405.5億元；經(jīng)濟(jì)效益型情景能夠維持生態(tài)系統(tǒng)服務(wù)水平現(xiàn)狀，還可以最大程度地滿足經(jīng)濟(jì)發(fā)展對土地的需求，經(jīng)濟(jì)效益為78 346.3億元；綜合效益型情景綜合考慮生態(tài)和經(jīng)濟(jì)發(fā)展，實現(xiàn)了兩者的平衡，生態(tài)系統(tǒng)服務(wù)價值總量和經(jīng)濟(jì)效益分別為108 068.2億元和75 480.7億元。

(2)綜合考慮3種情景，從整體的發(fā)展趨勢上來看，與現(xiàn)狀相比到2030年時，林地增加部分主要分布在攀西地區(qū)；水域增加部分主要分布在攀西地區(qū)；建筑用地增加部分主要分布在攀西地區(qū)和四川盆地。

(3)研究結(jié)果可以為四川省土地利用規(guī)劃和管理提供參考，本文的方法也可以應(yīng)用于我國的其他省份土地利用研究。在未來的研究中，可以通過收集更加豐富的土地利用數(shù)據(jù)，完善模型中的轉(zhuǎn)換規(guī)則；使用不同的方法計算預(yù)測年份的各土地利用類型面積大?。粚⒛Ｐ椭惺褂玫倪壿嫽貧w替換為隨機(jī)森林來完成“適宜性”的計算。