沈振濤，趙 強(qiáng)，高 倩，武 瑋，劉玉玉

(濟(jì)南大學(xué) 水利與環(huán)境學(xué)院，山東 濟(jì)南 250022)

土地利用/土地覆被變化(LUCC)作為社會(huì)、經(jīng)濟(jì)與自然環(huán)境相互作用的表現(xiàn)形式[1]，在全球環(huán)境變化、可持續(xù)發(fā)展研究中占有重要地位[2]。LUCC模型是探索未來(lái)土地利用演變過(guò)程的特征規(guī)律，揭示人類(lèi)活動(dòng)影響下的生態(tài)環(huán)境變化的重要方法[3]。

目前，用于土地利用變化模擬的方法有很多，常見(jiàn)的有利用變化及效應(yīng)(CLUE-S)模型[4]、馬爾可夫(Markov)預(yù)測(cè)模型[5]、元胞自動(dòng)機(jī)(CA)-Markov模型[6]、系統(tǒng)動(dòng)力學(xué)模型[7]、邏輯回歸(LR)模型[8]、人工神經(jīng)網(wǎng)絡(luò)(ANN)[9]等，其中的CA-Markov模型被公認(rèn)為對(duì)模擬和預(yù)測(cè)土地利用變化具有重要價(jià)值[10]。CA-Markov模型可以整合社會(huì)、經(jīng)濟(jì)、地形和道路數(shù)據(jù)并進(jìn)行模擬，從而反映社會(huì)經(jīng)濟(jì)發(fā)展的特征[11]，因此，許多學(xué)者使用該模型來(lái)預(yù)測(cè)和優(yōu)化未來(lái)的土地利用，例如：周浩等[12]基于1990、2013年撓力河流域的土地利用數(shù)據(jù)，對(duì)2025年土地利用格局進(jìn)行模擬預(yù)測(cè)；侯西勇等[13]基于1990、2000年河西走廊土地利用數(shù)據(jù)，模擬了20世紀(jì)90年代該地區(qū)土地利用變化；靳含等[14]運(yùn)用多時(shí)間跨度模擬四道河子鎮(zhèn)未來(lái)的土地利用變化。目前CA-Markov模型的應(yīng)用較為廣泛，且在土地利用演化模擬中發(fā)揮了重要作用。

本文中以大汶河流域?yàn)檠芯繉?duì)象，選取1985、1995、2005、2017年遙感影像，利用完整的遙感圖像處理平臺(tái)(ENVI)進(jìn)行解譯，結(jié)合土地利用變化的面積和土地利用動(dòng)態(tài)度指標(biāo)分析研究區(qū)內(nèi)1985—2017年土地利用的動(dòng)態(tài)變化特征，最后運(yùn)用CA-Markov模型對(duì)研究區(qū)未來(lái)的土地利用格局進(jìn)行模擬預(yù)測(cè)，為大汶河流域土地利用資源的合理運(yùn)用及生態(tài)環(huán)境保護(hù)提供參考。

1 研究區(qū)概況

大汶河又被稱(chēng)為古汶水，位于山東中部，是黃河下游水系中最大且最后的一條支流。大汶河發(fā)源地為山東旋崮山北麓，自東向西流經(jīng)山東省境內(nèi)7個(gè)縣(市)[15]。大汶河流域位于東經(jīng)116°11′15″—118°0′0″、北緯35°37′30″—36°32′30″之間，東部以魯山為界，西部以黃河流域?yàn)榻?，南部以蒙山及余脈和淮河流域?yàn)榻?，北部主要以泰山山脈和小清河為界，流域全長(zhǎng)209 km，面積8 634 km2。平均海拔為213 m，地勢(shì)北高南低，東高西低，南北兩面環(huán)山[16]。大汶河流域地理位置如圖1所示。

2 數(shù)據(jù)來(lái)源與研究方法

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

采用1985、1995、2005年Landsat 4-5 TM數(shù)據(jù)和2017年Landsat 8 OLI_TIRS數(shù)據(jù)對(duì)研究區(qū)進(jìn)行土地分類(lèi)。為了便于進(jìn)行土地分類(lèi)，選取每年6月份且云量較少時(shí)的數(shù)據(jù)作為基礎(chǔ)數(shù)據(jù)，運(yùn)用ENVI軟件對(duì)遙感數(shù)據(jù)進(jìn)行輻射糾正、幾何校正等處理，參照GB/T 21010—2007《土地利用現(xiàn)狀分類(lèi)》，以及遙感影像分辨率、土地利用方式、土地利用實(shí)際情況和后期模型分析的需要，將研究區(qū)土地利用類(lèi)型分為6種，分別為建筑用地(URHD)、林地(FRST)、未利用地(BARR)、水體(WATR)、草地(PAST)、耕地(AGRL)。

2.2 研究方法

2.2.1 土地利用變化特征指標(biāo)

凈轉(zhuǎn)移量是土地轉(zhuǎn)入與轉(zhuǎn)出面積之差，總轉(zhuǎn)移量是指各地類(lèi)在某時(shí)間段的轉(zhuǎn)入和轉(zhuǎn)出面積之和，相關(guān)計(jì)算公式分別為

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

式中：Nc為轉(zhuǎn)移凈變化率，%；ΔUin為其他地類(lèi)轉(zhuǎn)化為該地類(lèi)的總面積，km2；ΔUout為該地類(lèi)轉(zhuǎn)化為其他地類(lèi)的面積，km2；Ua、Ub分別為該地類(lèi)在時(shí)間段內(nèi)初始、末期面積，km2；Tc為轉(zhuǎn)移總變化率，%；K為單一土地利用類(lèi)型的動(dòng)態(tài)變化率，%；T為時(shí)間，a；Ps為轉(zhuǎn)移狀態(tài)；Sa為總土地利用凈轉(zhuǎn)移變化率，%；ΔUin,i為其他土地利用類(lèi)型轉(zhuǎn)化為類(lèi)型i的總面積，km2；ΔUout,i為土地利用類(lèi)型i轉(zhuǎn)化為其他類(lèi)型的面積，km2；St為轉(zhuǎn)移總變化率，%；Lc為綜合土地利用動(dòng)態(tài)度；ΔUai,j為某時(shí)間段內(nèi)土地利用類(lèi)型i轉(zhuǎn)移為類(lèi)型j的面積，km2；ΔUai為土地利用類(lèi)型i的總面積，km2；n為土地利用類(lèi)型總數(shù)；Pt為某時(shí)間段內(nèi)的轉(zhuǎn)移狀態(tài)。

2.2.2 CA-Markov模型

Markov模型的原理是基于柵格的空間概率，根據(jù)現(xiàn)有狀況預(yù)測(cè)未來(lái)各個(gè)時(shí)刻的變化[17]，已被廣泛應(yīng)用于土地利用變化模擬[18]、景觀格局變化模擬[19]和森林和植被的可持續(xù)性模擬[20]。利用遙感與地理信息結(jié)合應(yīng)用系統(tǒng)(IDRISI Selva)中的CA-Markov模塊，以1995、2005年土地利用狀況為初始矩陣，得出1995—2005年間大汶河流域土地利用面積轉(zhuǎn)移矩陣并將其作為模型轉(zhuǎn)移規(guī)則，驗(yàn)證2017年土地利用狀況，然后以2005、2017年土地利用狀況為初始矩陣，得出2005—2017年間土地面積轉(zhuǎn)移矩陣，并預(yù)測(cè)研究區(qū)2030年土地利用狀況。

3 大汶河流域土地利用變化特征

3.1 土地利用結(jié)構(gòu)變化特征

根據(jù)大汶河流域1985、1995、2005、2017年共4期土地利用狀況，基于ArcGIS的空間分析功能，對(duì)研究區(qū)4期的土地利用類(lèi)型圖進(jìn)行處理，同時(shí)提取大汶河流域1985、1995、2005、2017年的4期土地利用數(shù)據(jù)，如圖2、表1所示。

表1 大汶河流域不同時(shí)期土地利用類(lèi)型面積對(duì)比

從相關(guān)數(shù)據(jù)可以看出，林地、耕地和草地是研究區(qū)主要的3種土地利用類(lèi)型，覆蓋了研究區(qū)的絕大部分地區(qū)。1985年，耕地面積所占比例最大，達(dá)到35.59%，面積為3 072.88 km2；之后是草地、林地，所占比例分別為28.12%、24.05%，面積分別為2 427.72、2 076.73 km2；建筑用地占比為7.73%，面積為667.78 km2；水體、未利用地所占比例相對(duì)較小，分別為1.21%、3.30%，面積分別為104.20、284.74 km2。在空間上，草地、耕地、林地交叉分布，林地遍布研究區(qū)各地，林地主要位于丘陵山區(qū)，建筑用地則主要分布于河流附近及平原地區(qū)，而耕地位于中間位置，主要依附水體且與建筑用地同存，未利用地主要零散分布于丘陵山區(qū)未開(kāi)墾處。

1995年，草地面積所占比例最大，達(dá)到37.01%，面積為3 195.31 km2；之后是耕地，所占比例為26.57%，面積為2 427.72 km2；林地、建筑用地占比分別為19.96%、10.69%，面積分別為1 723.63、923.33 km2；水體、未利用地所占比例分別為1.36%、4.40%，面積分別為117.81、379.71 km2。與1985年相比，草地面積增加了8.89%，耕地和林地面積分別減少了9.02%、4.09%，建筑用地面積增加了2.96%，未利用地和水體面積變化相對(duì)較小。

2005年，各土地利用類(lèi)型中所占比例最大的依然是草地，達(dá)到43.71%，面積為3 773.67 km2；林地和耕地次之，所占比例分別為20.95%、19.64%，面積分別為1 808.53、1 695.54 km2；建筑用地占比分別為13.91%，面積為1 201.01 km2；水體、未利用地所占比例分別為1.32%、0.48%，面積分別為114.27、41.04 km2。與1995年相比，草地面積持續(xù)增加，增加率為6.7%；耕地面積持續(xù)減少，減少率為6.93%；林地、建筑用地面積分別增加了0.99%、3.22%；未利用地面積減少了3.92%，水體面積變化相對(duì)較小。

2017年，耕地面積所占比例最大，達(dá)到38.42%，面積為3 317.32 km2；草地次之，所占比例為32.16%，面積為2 776.88 km2；之后是林地、建筑用地，所占比例分別為11.84%、16.20%，面積分別為1 022.22、1 398.63 km2；水體、未利用地所占比例較小，分別為1.34%、0.04%，面積分別為115.9、3.09 km2。與2005年相比，草地面積減少了11.55%，耕地面積增加了18.78%，林地面積減少了9.11%，建筑用地面積增加了2.29%，未利用地面積減少了0.44%，水體面積變化依然不大。

為了直觀分析大汶河流域近30年來(lái)的土地利用變化，從不同土地利用類(lèi)型在各年中所占比例(表1)可知：1985—2017年間，研究區(qū)建筑用地面積逐年增加，從667.78 km2擴(kuò)大到1 398.63 km2；林地面積只在2005年有微弱增加，總體呈減少趨勢(shì)，從2 076.73 km2減少到1 022.22 km2；耕地面積先減少后增加，與1985年相比，2017年耕地面積增加了244.44 km2；草地面積變化與耕地相反，呈先增加后減少趨勢(shì)，與1985年相比，2017年草地面積仍增加了349.16 km2；未利用地面積逐漸減少，而水體面積無(wú)明顯變化。

3.2 土地利用動(dòng)態(tài)轉(zhuǎn)移過(guò)程

在近35 a間，大汶河流域存在幾十種土地利用轉(zhuǎn)移關(guān)系，其中絕大部分來(lái)自建筑用地、林地、耕地以及草地的轉(zhuǎn)移(見(jiàn)表2)。由表可見(jiàn)：1985—2017年間，建筑用地主要轉(zhuǎn)化為耕地，轉(zhuǎn)化面積為186.43 km2，轉(zhuǎn)化比例為28%；林地主要轉(zhuǎn)化為耕地，轉(zhuǎn)化面積為670.14 km2，轉(zhuǎn)化比例達(dá)到32%；林地向草地轉(zhuǎn)化的面積為354.04 km2，轉(zhuǎn)化比例為17%；未利用地主要轉(zhuǎn)化為草地、耕地，轉(zhuǎn)化面積分別為161.41、87.34 km2，轉(zhuǎn)化比例分別為58%、31%；耕地主要轉(zhuǎn)化為草地，其次轉(zhuǎn)化為建筑用地,轉(zhuǎn)化面積分別為1 102.33、536.76 km2，轉(zhuǎn)化比例分別為36%、17%；草地則主要轉(zhuǎn)化為耕地，轉(zhuǎn)化面積為1 038.96 km2，轉(zhuǎn)化比例為43%。水體面積的轉(zhuǎn)化量相對(duì)較小。

表2 1985—2017年大汶河流域土地利用變化轉(zhuǎn)移矩陣 km2

根據(jù)式(1)—(8)計(jì)算可得到大汶河流域1985—2017年土地利用動(dòng)態(tài)度和土地利用轉(zhuǎn)移態(tài)勢(shì)，如表3所示。由表中數(shù)據(jù)可知，大汶河流域在1985—2017年間，土地利用動(dòng)態(tài)度Lc為0.91，土地利用轉(zhuǎn)移狀態(tài)Pt為0.27，流域整體土地利用屬于弱平衡型轉(zhuǎn)換；水體、耕地以及草地Ps分別為0.13、0.07和0.11，說(shuō)明水體、耕地和草地呈緩慢擴(kuò)張趨勢(shì)；建筑用地Ps為0.52，說(shuō)明建筑用地呈擴(kuò)張趨勢(shì)，且變化趨勢(shì)較為明顯。林地Ps為-0.73，說(shuō)明大汶河流域林地呈顯著減少趨勢(shì)；未利用地Ps為-0.98，說(shuō)明未利用地顯著減少，呈逆向轉(zhuǎn)換。

表3 1985—2017年大汶河流域土地利用動(dòng)態(tài)度和轉(zhuǎn)移狀態(tài)

3.3 土地利用驅(qū)動(dòng)力分析

人口與經(jīng)濟(jì)都是城市化發(fā)展進(jìn)程的主要影響因素。在人類(lèi)發(fā)展過(guò)程中，為了滿(mǎn)足生產(chǎn)、生活的需要，人們對(duì)土地的開(kāi)發(fā)和利用改變了當(dāng)?shù)赝恋乩妙?lèi)型和分布。大汶河自東向西流經(jīng)萊蕪、泰安市，最后注入東平湖。2005年泰安市總?cè)丝跒?49.95萬(wàn)，到2017年全市總?cè)丝谶_(dá)到570.70萬(wàn)，13年間人口增加20.75萬(wàn)；2005年萊蕪市總?cè)丝跒?24.29萬(wàn)，2017年全市總?cè)丝谶_(dá)到129.68萬(wàn)，13年間人口增加5.39萬(wàn)。流域內(nèi)人口數(shù)量的增加使大汶河流域土地利用格局發(fā)生變化，建筑用地面積不斷擴(kuò)大。經(jīng)濟(jì)發(fā)展過(guò)程實(shí)際上是土地利用結(jié)構(gòu)的優(yōu)化調(diào)整與土地利用集約度提高的過(guò)程。2017年，泰安、萊蕪市固定資產(chǎn)分別達(dá)到2 992.5、568.14億元，經(jīng)濟(jì)的快速發(fā)展，使得市區(qū)建筑面積不斷擴(kuò)大，計(jì)算結(jié)果顯示，2005—2017年研究區(qū)建筑用地增加16.45%，草地和未利用地轉(zhuǎn)化為建筑用地的比例分別為11.98%和7.58%。政策因素是短期內(nèi)推動(dòng)土地利用類(lèi)型變化的關(guān)鍵因素之一。1987年的《泰安市林業(yè)區(qū)劃》中規(guī)定，山區(qū)建設(shè)要以水源涵養(yǎng)林為主，使林地的種植面積得到了增加。2006年的《萊蕪市土地利用總體規(guī)劃》中規(guī)定，要優(yōu)化土地利用空間格局，加強(qiáng)中心城區(qū)土地利用控制，增強(qiáng)農(nóng)用地保護(hù)和基本農(nóng)田建設(shè)，推進(jìn)建筑用地節(jié)約和集約利用等，對(duì)萊蕪市土地利用變化起到了控制作用。

4 大汶河土地利用變化預(yù)測(cè)

4.1 驅(qū)動(dòng)力因子選取

驅(qū)動(dòng)力因子是影響土地利用模擬的重要因素，因此選取適宜的驅(qū)動(dòng)力因子對(duì)模型的模擬非常重要[21]。結(jié)合研究區(qū)經(jīng)濟(jì)、生態(tài)條件，參考相關(guān)研究結(jié)果[22-24]，本文中選取了自然及社會(huì)因素共10種驅(qū)動(dòng)力因子(見(jiàn)表4)。A1、A2分別為高程和坡度，代表自然因素；A3—A10代表社會(huì)因素，其中A4—A9

表4 元胞自動(dòng)機(jī)-馬爾科夫(CA-Markov)模擬采用的驅(qū)動(dòng)力因子

表示道路到達(dá)的可能性，A10代表與河流的距離，表示植被獲取水資源的可能性。利用ArcGIS軟件將各類(lèi)影響因子與土地利用進(jìn)行統(tǒng)計(jì)分析，對(duì)土地利用的影響進(jìn)行歸一化處理，結(jié)果見(jiàn)圖3。

4.2 適宜性圖集

土地適宜性是土地在一定條件下對(duì)不同用途的適應(yīng)程度[25]。在模擬土地利用轉(zhuǎn)化過(guò)程的適宜性時(shí)，選取了以上10個(gè)驅(qū)動(dòng)力因子。利用邏輯回歸計(jì)算土地利用狀況跟驅(qū)動(dòng)力因子之間的關(guān)系，然后依據(jù)回歸關(guān)系得到全部6種土地利用類(lèi)型的適宜性圖集(見(jiàn)圖4)。該過(guò)程利用IDRISI Selva軟件進(jìn)行模擬。

適宜性圖集中的因變量是展開(kāi)的像元，以二進(jìn)制柵格表示，其中1表示特定像元的變化，說(shuō)明可轉(zhuǎn)化為該類(lèi)地；0表示這一段時(shí)間內(nèi)沒(méi)有變化，說(shuō)明可轉(zhuǎn)化為其他類(lèi)地[26]。從圖4中可以看出，耕地和草地轉(zhuǎn)化的適宜性變化最大，而建筑用地、林地及未利用地轉(zhuǎn)化的適宜性變化較小，水體以保持原樣為主。

4.3 模型精度

Kappa系數(shù)可以用來(lái)檢驗(yàn)?zāi)M結(jié)果與實(shí)際數(shù)據(jù)的一致性，因此可用于2個(gè)圖件一致性評(píng)價(jià)和遙感解譯精度評(píng)價(jià)[27]。Kappa系數(shù)K的計(jì)算公式為

(9)

式中：Po為模擬正確的柵格比例；Pc為期望模擬正確的柵格比例；Pp理想狀態(tài)下模擬正確的柵格比例；N為柵格總數(shù)；N1為模擬正確的柵格數(shù)。

基于1995、2005年土地利用數(shù)據(jù)，利用CA-Markov模型預(yù)測(cè)大汶河流域2017年土地利用狀況，并將2017年土地利用類(lèi)型圖與模型的模擬圖對(duì)比，進(jìn)行CA-Markov模型精度檢驗(yàn)。大汶河流域2017年土地利用類(lèi)型圖共有柵格個(gè)數(shù)為9 586 156。

利用IDRISI軟件中的CROSSTAB模塊計(jì)算模型精度，得到Kappa系數(shù)為0.705，表明模型可以較好地模擬大汶河流域2017年土地利用狀況，且精度較高，可以用來(lái)模擬研究區(qū)未來(lái)土地利用狀況。

4.4 未來(lái)土地利用預(yù)測(cè)

以2017年土地利用狀況為基準(zhǔn)，利用2005—2017年土地利用轉(zhuǎn)移矩陣以及完成的適宜性圖集，利用CA-Markov模型，循環(huán)13次，預(yù)測(cè)大汶河流域2030年土地利用狀況，結(jié)果如圖5所示。

將1985—2030年間共5期土地利用類(lèi)型圖中各類(lèi)地所占比例進(jìn)行對(duì)比，結(jié)果如圖6所示。與2017年相比，大汶河流域草地面積有所增加，耕地、未利用地以及水體面積基本不變，林地面積大幅度減少，而建筑用地有所增加。未來(lái)研究區(qū)草地面積呈先增加、后減少的趨勢(shì)，到2030年又有所增加；耕地面積則正好相反，先減少、后增加，但2017—2030年間建筑用地和耕地面積逐漸增加；未利用地面積在1995—2005年間大幅度減少，直到2030年變化都不大；未來(lái)大汶河流域水體面積基本保持不變。預(yù)測(cè)結(jié)果顯示，未來(lái)一段時(shí)間內(nèi)研究區(qū)建筑用地增加有所放緩，說(shuō)明通過(guò)政府的政策影響對(duì)萊蕪市土地利用變化起到了控制作用，因此要嚴(yán)格貫徹執(zhí)行《萊蕪市土地利用總體規(guī)劃》，積極響應(yīng)國(guó)家號(hào)召，加強(qiáng)中心城區(qū)土地利用控制，加大農(nóng)用地保護(hù)和基本農(nóng)田建設(shè)、保護(hù)生態(tài)環(huán)境可持續(xù)發(fā)展的力度，推進(jìn)建筑用地節(jié)約和集約利用。

大汶河流域地圖從國(guó)家測(cè)繪地理信息局標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載，地圖審批編號(hào)為GS(2019)1697(http://bzdt.ch.mnr.gov.cn/browse.html?picId=%224o28b0625501ad13015501ad2bfc0258%22)，經(jīng)過(guò)Arcgis 10.7軟件數(shù)字化處理后得到。圖5 大汶河流域2030年土地利用狀況預(yù)測(cè)

圖6 1985—2030年大汶河流域土地利用類(lèi)型變化

5 結(jié)論

本文中以大汶河流域?yàn)檠芯繉?duì)象，通過(guò)分析1985—2017年土地利用變化情況、影響因子的驅(qū)動(dòng)作用進(jìn)行分析，并結(jié)合CA-Markov模型對(duì)2030年研究區(qū)土地利用變化趨勢(shì)進(jìn)行了模擬預(yù)測(cè)，得到主要結(jié)論如下：

1)通過(guò)對(duì)1985年以來(lái)大汶河流域土地利用變化進(jìn)行分析，水體、耕地和草地呈緩慢擴(kuò)張趨勢(shì)，建筑用地呈明顯擴(kuò)張趨勢(shì)，林地呈顯著減少趨勢(shì)，未利用地呈逆向轉(zhuǎn)換。究其原因，人口數(shù)量增加和經(jīng)濟(jì)發(fā)展使得大汶河流域土地利用格局發(fā)生變化，城市化進(jìn)程處于不斷擴(kuò)張階段。

2)對(duì)大汶河流域2017年土地利用進(jìn)行驗(yàn)證，Kappa系數(shù)為0.705，說(shuō)明利用CA-Markov模型可以較好地模擬大汶河流域2017年土地利用狀況，且精度較高，可以用來(lái)模擬未來(lái)土地利用狀況。

3)運(yùn)用CA-Markov模型對(duì)大汶河流域2030年土地利用變化進(jìn)行模擬，結(jié)果顯示，草地面積逐漸增加，耕地、未利用地以及水體面積基本沒(méi)有變化，林地面積大幅度減少，而建筑用地有所增加，表明政府的政策調(diào)控措施優(yōu)化了中心城區(qū)土地利用格局，對(duì)城市土地利用變化起到了一定的控制作用。