沈振濤，趙 強，高 倩，武 瑋，劉玉玉

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

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

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

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

1 研究區(qū)概況

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

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

2.1 數(shù)據(jù)來源

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

2.2 研究方法

2.2.1 土地利用變化特征指標

凈轉移量是土地轉入與轉出面積之差，總轉移量是指各地類在某時間段的轉入和轉出面積之和，相關計算公式分別為

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

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

2.2.2 CA-Markov模型

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

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

3.1 土地利用結構變化特征

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

表1 大汶河流域不同時期土地利用類型面積對比

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

1995年，草地面積所占比例最大，達到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%，未利用地和水體面積變化相對較小。

2005年，各土地利用類型中所占比例最大的依然是草地，達到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%，水體面積變化相對較小。

2017年，耕地面積所占比例最大，達到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年來的土地利用變化，從不同土地利用類型在各年中所占比例(表1)可知：1985—2017年間，研究區(qū)建筑用地面積逐年增加，從667.78 km2擴大到1 398.63 km2；林地面積只在2005年有微弱增加，總體呈減少趨勢，從2 076.73 km2減少到1 022.22 km2；耕地面積先減少后增加，與1985年相比，2017年耕地面積增加了244.44 km2；草地面積變化與耕地相反，呈先增加后減少趨勢，與1985年相比，2017年草地面積仍增加了349.16 km2；未利用地面積逐漸減少，而水體面積無明顯變化。

3.2 土地利用動態(tài)轉移過程

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

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

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

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

3.3 土地利用驅動力分析

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

4 大汶河土地利用變化預測

4.1 驅動力因子選取

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

表4 元胞自動機-馬爾科夫(CA-Markov)模擬采用的驅動力因子

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

4.2 適宜性圖集

土地適宜性是土地在一定條件下對不同用途的適應程度[25]。在模擬土地利用轉化過程的適宜性時，選取了以上10個驅動力因子。利用邏輯回歸計算土地利用狀況跟驅動力因子之間的關系，然后依據(jù)回歸關系得到全部6種土地利用類型的適宜性圖集(見圖4)。該過程利用IDRISI Selva軟件進行模擬。

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

4.3 模型精度

Kappa系數(shù)可以用來檢驗模擬結果與實際數(shù)據(jù)的一致性，因此可用于2個圖件一致性評價和遙感解譯精度評價[27]。Kappa系數(shù)K的計算公式為

(9)

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

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

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

4.4 未來土地利用預測

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

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

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

圖6 1985—2030年大汶河流域土地利用類型變化

5 結論

本文中以大汶河流域為研究對象，通過分析1985—2017年土地利用變化情況、影響因子的驅動作用進行分析，并結合CA-Markov模型對2030年研究區(qū)土地利用變化趨勢進行了模擬預測，得到主要結論如下：

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

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

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