基于CA-Markov模型的哈爾濱市土地利用變化預測研究*

李志明，宋 戈.2※，魯 帥，王 蓓，徐四桂

(1.東北農業大學資源與環境學院，黑龍江哈爾濱 150030；2.東北大學土地管理研究所，遼寧沈陽 110819)

0 引言

土地利用變化是人類經濟社會發展對土地資源利用的直接反映，是全球環境變化和可持續發展研究的主要內容[1-2]。加速的土地利用變化產生更加復雜的土地利用格局，直接或間接引發生態環境變化[3]，如土壤質量[4-5]、氣候條件[6-7]、水文狀況[8]等。因此，科學合理地分析土地利用演變特征，探究未來土地利用變化趨勢，能夠為區域土地資源的合理開發利用、生態環境的改善和保護、經濟的可持續發展提供科學的參考依據[9-10]。目前，國內外土地利用變化預測模型主要包括側重時間維度分析的系統動力學模型、Markov模型、Logistic模型、神經網絡模型等和側重空間維度分析的CA模型、CLUE模型、CLUE-S模型等。CA-Markov模型由于結合了Markov模型的時間維度分析優勢和CA的空間維度分析能力[11]，得到了廣泛的應用和探討。研究尺度側重于中尺度的流域[11-13]和城市[14-15]，研究內容集中于模型指標選取對預測精度的影響[16]和未來土地利用趨勢對區域生態環境的影響[13]等方面，但多數研究中將Markov模型預測的條件概率圖像、較單一的影響因子或單生態方面的影響因子作為元胞自動機的轉換規則，轉換規則過于簡單，缺乏社會經濟因子對土地利用變化影響的考慮，使結果不能更好地貼合實際發展情況。因此，文章以哈爾濱市建成區為例，選取經濟和生態兩方面的GDP、坡度、水域、與道路距離、與水域距離、與建設用地距離共6種影響因子，采用MCE模型中的加權線性合并法單獨制作耕地和建設用地適宜性圖集參與到元胞自動機轉換規則中。應用GIS技術、Markov模型和土地利用重心遷移模型探究1992～2014年土地利用時空演變規律，最后基于CA-Markov模型對2025年研究區土地利用變化進行更加細致的預測。

1 研究區概況和數據來源

1.1 研究區概況

哈爾濱市位于黑龍江省南部，地跨125°42′E～130°10′E，44°04′N～46°40′N，地勢東南高，西北低，東南地區多山地丘陵，西北地區位于松嫩高平原上，平坦低洼，松花江貫穿其中。年平均降水570mm，四季分明，屬于溫帶大陸性季風氣候。全市總面積5.3萬km2， 2014年總人口達987萬，是重要的工業城市，城鎮化發展迅速。該文選取哈爾濱市的道外區、道里區、南崗區、阿城區、松北區、呼蘭區、香坊區和平房區8個主要市轄區為研究區。該區域位于哈爾濱市西北部，總面積約7 083km2，占哈爾濱市總面積的13.36%，土地利用類型豐富，以耕地和林地為主， 2014年土地利用情況與1992相比，耕地面積減少了10.91%，建設用地增加了74.37%，是研究土地利用變化的熱點區域。

1.2 數據來源

選取研究區1992年、2003年和2014年3期遙感影像(1992年和2003年為Landsat 5 TM影像， 2014年為Landsat 8 OLI影像，分辨率均為30m)為基礎數據(數據下載于http：//www.gscloud.cn/)，利用ENVI 5.1對遙感影像進行波段融合、圖像校正等預處理后，應用ArcGIS10.2對數據進行解譯，解譯過程借助地形圖等相關材料進行校準，參考《第二次全國土地調查土地分類》標準和研究區實際情況，將土地利用類型分為耕地、林地、草地、水域、建設用地和未利用地6個一級分類(圖1)。其他數據包括研究區30m分辨率的DEM數據，哈爾濱市社會經濟統計數據等。

2 研究方法

基于研究區3期土地利用解譯數據，運用GIS空間分析技術，利用Markov模型和土地利用重心遷移模型分析1992～2014年土地利用動態變化特征，選取一定約束條件和適宜性因子，結合MCE模型，建立適宜性圖集，采用CA-Markov模型對土地利用變化進行預測研究，并應用土地利用動態度對預測結果進行評價。

2.1 土地利用重心遷移模型

土地利用重心是在人口分布重心理論的背景下提出，用于表征土地開發利用過程的重心空間位置以及變化情況。計算方法如下：

(1)

式(1)中，Xk、Yk分別表示k種土地利用類型重心X、Y坐標值；Aki表示k種土地利用類型i圖斑的面積；Xki、Yki分別表示k種土地利用類型i圖斑的幾何中心X、Y坐標值；Ak表示k種土地利用類型總面積。

2.2 MCE模型

MCE是為實現一個具體目標，而設定一系列的評價準則的過程[17]，其最基本的問題是如何整合不同的決策準則信息，即如何整合約束條件和適宜性因子兩類決策準則信息。該文采用MCE中的加權線性合并法，能夠較好的表達連續性的適宜性因子，其運算過程如下：

S=∑wixi∏cj

(2)

式(2)中，S為適宜性指數；wi為第i個適宜性因子的權重指數；xi是第i個適宜性因子的得分。該文將適宜性因子統一模糊處理至0～255范圍內，數值越大，適宜度越高；cj是約束條件j的數值，其值為0或1，表示不流轉或可流轉。

結合研究區經濟、生態條件，選取GDP、坡度、水域、與道路距離、與水域距離、與建設用地距離6種影響因子，針對建設用地和耕地，單獨設置約束條件和適宜性因子用于制作適宜性圖集，其余土地利用類型由于限制因素較少，故選擇Markov模型生成的條件概率矩陣作為適宜性圖集。具體如下：

(1)建設用地。設定約束條件為： ①水域，考慮到水資源短缺、水域污染等生態環境問題，出于對研究區的環境保護，將水域設定為不發生用地類型轉換，賦值為0，其余賦值為1； ②坡度，根據《水土保持工作條例》和研究區自然條件，設定坡度大于25°為禁止開發區域，賦值為0，其余賦值為1。設定適宜性因子為： ①與建設用地距離，建設用地對周邊土地的轉入有較高的吸引，且隨著距離的增大，吸引力不斷弱化，直至形成無差別影響。故設定500m內適宜流轉，適宜度隨距離的增加單調遞[12]； ②與水域距離，出于對環境的保護，不考慮距水域100m范圍內土地的流轉， 100m外隨距離增加，轉變為建設用地適宜度增高[13]； ③與道路距離，道路對經濟發展的影響顯著，設定研究區主要道路100m內，適宜度不斷減弱， 100m外適宜度無差別[10]； ④GDP，由于建設用地對GDP貢獻最大，故設定GDP越高，建設用地轉入適宜度越高。

(2)耕地。設定的約束條件與建設用地相同。設定適宜性因子為： ①與建設用地距離， 500m內適宜性最弱， 500m外適宜度不斷增加； ②與水域距離，考慮到水源對耕地的重要性，設置100m內，適宜度最高， 100m外適宜度逐漸減弱； ③與道路距離，與主要道路100m內適宜度最弱， 100m外適宜度逐漸增強； ④GDP，耕地對GDP的貢獻相對較小，故設定適宜度隨GDP的降低而降低。

2.3 CA-Markov模型及預測方法

Markov模型是用不同時間的系統狀態給出一個從一種土地利用類型到另一種類型的轉換概率矩陣，預測未來某一時刻的系統狀態，是一種空間概率模型，強調的是從時間維度上分析土地利用的變化趨勢。CA模型是基于元胞的前一個狀態及周邊鄰域狀態按一定轉換規則預測未來元胞狀態的模型，是局部網絡動力學模型，其主要強調在空間維度上進行預測。CA-Markov模型結合Markov模型時間維度分析的優勢和CA模型在空間維度分析的優勢，降低了制定轉換規則的難度，減少了人為因素的干擾。具體預測過程如下：

(1)計算土地轉移矩陣。將1992年和2003年土地利用解譯數據進行疊置分析，利用Markov模型計算土地轉移概率矩陣和面積矩陣，并作為轉換規則參與后面預測。

(2)適宜性圖集的建立。應用MCE整合不同準則標準，得到不同土地類型的轉移適宜性圖，合并成適宜性圖集，參與土地利用變化的預測。

(3)參數設定。該文以2003年為預測時間起點，設定5×5的元胞濾波器，對2014年土地利用空間分布進行預測。在預測精度較好的情況下，以2014年為起始年，進行2025年的土地利用空間分布預測。

2.4 土地利用動態度

土地利用動態度是反映土地利用變化程度的重要指標，其動態模型可分為單一和綜合土地利用動態度，該文選擇單一土地利用動態度指標進行預測結果的評價，可以反映某一土地利用類型變化的劇烈程度。

(3)

式(3)中，i表示土地利用類型；Ki表示一定時間段內i土地利用類型的動態度；At2和At1分別表示變化目標年和變化基準年土地利用面積；t2-t1表示研究期時間間隔。

3 結果分析

3.1 土地利用結構變化分析

基于ArcGIS10.2分析了1992、2003和2014年3期解譯的土地利用數據，得到土地利用情況(表1)。

由表1可見，研究區內耕地占主體地位，所占比例分別為72.35%、70.45%和64.45%，在此期間，耕地以每年0.35個百分點的速率不斷減少，共減少5.592 7萬hm2，耕地壓力不斷增加。草地和未利用地也呈現減少態勢，減少速率較為緩和。林地和水域先減少后增加， 2014年相較于1992年所占比例有所提升，體現出近年來哈爾濱對生態環境保護工作的重視及取得的成果。建設用地不斷增加， 22年期間增加5.64個百分點，增加面積為3.996 8萬hm2，反映出在城鎮化水平不斷提高、人口壓力增大、經濟不斷發展的背景下，對建設用地擴展的需求。

表1 1992年、2003年和2014年研究區各類土地利用面積及所占比例

土地利用類型199220032014面積(萬hm2)百分比(%)面積(萬hm2)百分比(%)面積(萬hm2)百分比(%)耕地51.248072.3549.900670.4545.655364.45林地11.377616.0611.178915.7811.415016.12草地0.48340.680.44290.630.33650.48水域2.14213.022.01392.844.01495.67建設用地5.36977.587.179310.149.363513.22未利用地0.21280.300.11790.170.04850.07

3.2 土地利用變化的空間轉換分析

單純的土地利用類型面積的增減不能較好地反映各類型間的轉換情況，該文通過IDRISI中的空間分析功能，得到1992～2003年和2003～2014年土地利用面積轉移矩陣(表2)，從而充分了解土地利用空間變化過程。

表2 1992～2014年研究區土地利用面積轉移矩陣 hm 2

1992～2003年間，耕地凈轉出1.350 4萬hm2，主要轉出為建設用地2.354 7萬hm2和林地1.294 0萬hm2； 林地凈轉出1 436hm2，轉入轉出基本平衡； 建設用地凈轉入1.810 7萬hm2，耕地轉入的2.354 7萬hm2是建設用地增加的主要原因； 水域少量減少，凈轉出1 811hm2，水域和耕地之間的相互轉換明顯； 草地和未利用地由于面積較小，總體轉入轉出不明顯。但空間位置變化明顯。2003～2014年間，耕地持續減少，凈轉出4.195 2萬hm2，主要轉出為建設用地2.932 2萬hm2，水域2.142 7萬hm2和林地1.298 3萬hm2，主要轉入未建設用地8 797hm2和林地7 399hm2； 林地由凈轉出變為凈轉入4 297hm2，主要由于轉出為耕地的面積減少和耕地的轉入增多； 建設用地不斷增加，增加速度高于上時段，凈增加2.184 6萬hm2，耕地仍然是建設用地重要來源，轉入2.932 2萬hm2，轉出主要為耕地8 797hm2； 水域較上時段明顯不同，不降反增加1.756 5萬hm2，轉入主要為耕地2.142 7萬hm2，總轉出由上時段的1.197 8萬hm2減少到4 942hm2，可見環境保護效果明顯。

圖1 1992、2003和2014年研究區土地利用情況

圖2 土地利用重心點空間位置

3.3 土地利用重心變化

從上述分析可知， 1992～2014年各土地利用類型中，耕地和建設用地變化最為顯著，同時也是區域經濟發展和糧食安全的重要影響因子，故該文以耕地和建設用地兩種地類為代表，分析其重心變化特征，借助GIS空間分析技術和式(1)得到圖2，并將2025年預測結果展布到圖2，用于后期預測分析。1992～2003年，耕地重心向東南方向移動356m，移動速度較慢； 建設用地向西移動2 307m，表明建設用地向四周擴張的同時，更偏向于哈爾濱建成區西南方向發展。2003～2014年，耕地重心向北移動886m，較上時間段速度明顯加快，未來幾年有加快移動的趨勢； 建設用地向南移動1 791m，移動速度減緩。整體來看， 1992～2014年間，受研究區北部糧食主產區的影響，耕地重心整體向北方向移動，向呼蘭區、巴彥縣方向靠攏，建設用地由于松花江流域的阻隔和氣候的限制，整體向有利于經濟發展的西南方向移動。

3.4 CA-Markov模型預測精度檢驗

模型的一個重要性步驟是檢驗。一般情況下，利用某一土地覆蓋狀況已知時段的預測，估算對一個過程的認識和模仿的功效，這就是用于驗證的測試。Kappa統計是較為常用的精度檢驗模型(式3)，其判別標準為Kappa值從0～1表示兩幅圖像從完全不吻合到完全吻合的轉變，其中Kappa≥0.75表明模型預測精度較高，Kappa≤0.4預測精度較差。該文應用CA-Markov模型預測的2014年土地利用類型預測結果(圖3b)與2014年已知土地利用數據(圖3a)進行Kappa驗證，結果為Kappa=0.878 0，檢驗精度較高，能夠進行后續預測研究，公式表達如下[18]。

(3)

式(3)中，p0表示土地利用分布預測正確比例；pc表示隨機情形下預測正確比例；pp表示理想狀態下預測完全正確的比例。

圖3 2014年研究區土地利用現狀圖(a)和預測圖(b)

圖4 2014年研究區土地利用現狀圖和2025年研究區土地利用預測

表3 2014～2025年研究區土地利用轉移概率矩陣

土地類型耕地林地草地水域建設用地未利用地耕地0.7410.0520.0020.0870.1180.000林地0.1710.7760.0000.0200.0330.000草地0.3390.0240.5360.0070.0950.000水域0.3280.0150.0000.6410.0150.000建設用地0.2410.0170.0020.0030.7360.001未利用地0.3820.0810.0030.0760.2100.249 注:橫向代表變化基準年2014年,縱向代表變化目標年2025年

表4 2003～2025年研究區土地利用變化情況

土地利用類型20252003～20142014～2025預測面積(萬hm2)面積變化(萬hm2)動態度(%)面積變化(萬hm2)動態度(%)耕地39.4527-4.2453-8.51-6.2026-13.59林地11.63790.23602.110.22301.95草地0.2895-0.1065-24.03-0.0470-13.96水域6.75432.001099.362.739568.23建設用地12.77642.184230.423.412936.45未利用地0.0398-0.0694-58.87-0.0087-17.96

3.5 CA-Markov模型預測結果分析

在模型預測精度檢驗檢驗合格的基礎上，以2014年土地利用數據為基礎，制作用于2025年預測的適宜性圖集，并與Markov模型預測的2014～2025年土地面積轉移矩陣(表3)同時作為轉換規則進行2025年研究區土地利用變化的預測，結果見圖3、圖4和表4。

預測結果表明， 2025年研究區耕地面積39.452 7萬hm2，延續1992年以來減少的趨勢，較2014年共減少6.202 6萬hm2，且減少速率與2003～2014年時段相比不降反增，耕地保護壓力加大，耕地重心加速向北移動，向糧食主產區靠近， 2014～2025年間預計移動1 780m； 林地面積11.637 9萬hm2，較2014年穩中有增，變化不顯著，空間分布多集中于東南部； 草地面積減少到2 895hm2，與2014年相比減少41hm2，相比于上一時間段減少速度降低，但動態度反映草地變化仍然活躍； 水域面積6.254 3萬hm2，增加2.739 5萬hm2，空間分布變化體現為在2014年原有水域的基礎上擴寬、延長，水域的大量增加將使生態環境得以改善； 建設用地面積12.776 4萬hm2，自1992年以來不斷增加， 2014～2025年間增加3.412 9萬hm2，增速高于2003～2014時段。土地利用變化活躍，空間分布結果表明，建設用地增加過程集聚效應顯著，增加的面積多分布于原建設用地周邊，且占用耕地現象明顯，建設用地重心繼續向西南方向移動， 2014～2025年預計移動1 114m，移動速度減緩； 未利用地面積398hm2，比2014年減少87hm2。

4 結論與討論

(1)該文基于Markov模型和GIS、RS技術，分析哈爾濱市建成區1992～2014年間土地利用類型演變特征。選取經濟生態方面的6個影響因子，水域和坡度兩個約束條件和GDP、與水域距離、與建設用地距離、與道路距離4個適宜性因子參與到轉換規則的制定中，假設預測時間段內未發生突發性事件，應用CA-Markov模型預測了2025年土地利用情況。試驗得到的2014年土地利用結果和實際解譯數據Kappa系數達0.878 0，表明轉換規則的制定可行，模型可信度較高。轉換規則的制定提高了預測結果的準確合理性，使預測結果能夠為土地利用規劃的修編和生態環境的保護提供決策支持。

(2)土地利用變化演變特征主要集中于： 1992～2014年間，耕地面積持續減少，耕地重心向北加速移動，向糧食主產區呼蘭、巴彥縣方向靠近； 建設用地由于耕地的大量轉入使其在原有空間位置上不斷擴張，建設用地重心向西南方向移動； 林地在數量和空間分布上相對穩定； 草地和未利用地面積變化微弱，但空間分布變化明顯； 水域面積先減后增，反映出環境保護政策中水域保護效果明顯，耕地和水域之間大量轉換。

(3)2025年預測表明，耕地繼續減少，且減少速率不見緩和，耕地保護壓力較大，耕地重心持續向北移動，且速度加快； 林地延續前期變化趨勢，穩中有升； 水域面積明顯增加，土地利用變化活躍，對未來環境改善有積極促進作用； 建設用地不斷增加，占用大量耕地，擴展明顯，建設用地重心持續向西南移動，移動速度減慢； 草地和未利用地總面積較小，呈零星分布。

研究中，對模型參數選擇采用統一的30m×30m元胞大小， 5×5濾波器，未考慮元胞大小以及濾波器對預測結果的影響； 土地利用變化受到社會、經濟、生態等多方面影響，該文由于數據的不足和部分因子定量化困難的原因，只考慮了經濟生態方面的6個影響因子，對影響土地利用的社會因子如永久基本農田、人口密度等分析不足。在未來的研究中，加大對上述問題的探索力度，實現更好的空間預測效果。

[1] 后立勝，蔡運龍.土地利用/覆被變化研究的實質分析與進展評述.地理科學進展， 2004， 23(6)： 96～104

[2] 王兵， 臧玲.我國土地利用/土地覆被變化研究近期進展.地域研究與開發， 2006, 25(2): 86～91

[3] Weber C，Puissant A.Urbanization pressure and modeling of urban growth：Example of the tunis metropolitan area.Remote Sensing of Environment， 2003， 86(1/2)： 341～352

[4] 李志， 劉文兆，鄭粉莉.基于CA-Markov模型的黃土塬區黑河流域土地利用變化.農業工程學報， 2010， 26(1)： 346～352， 391

[5] Bormann H，Breuer L，Graff T，et al.Analysing the effects of soil properties changes associated with land use changes on the simulated water balance：a comparison of three hydrological catchment models foe scenario analysis.Ecol Model， 2007， 209(1)； 29～40

[6] 劉紀遠， 邵全琴，延曉冬，等.土地利用變化對全球氣候影響的研究進展與方法初探.地球科學進展， 2011, 26(10): 1015～1022

[7] 趙麗， 張蓬濤，周智.環首都貧困區土地利用變化對生態系統的影響及驅動因素分析——以河北淶水縣為例.中國農業資源與區劃， 2013, 34(5): 74～81

[8] Brath A，Montanan A，Moretti G.Assessing the effect on flood frequency of land ues change via hydrological simulation(with uncertainty).Journal of Hydrology， 2006， 324(1/2/3/4)： 141～153

[9] 汪佳莉， 吳國平，范慶亞，等.基于CA-Markov模型的山東省臨沂市土地利用格局變化研究及預測.水土保持研究， 2015, 22(1): 212～216

[10]趙冬玲， 杜萌，楊建宇，等.基于CA-Markov模型的土地利用演化模擬預測研究.農業機械學報， 2016, 47(3): 278～285

[11]王友生， 余新曉，賀康寧，等.基于CA-Markov模型的藉河流域土地利用變化動態模擬.農業工程報， 2011， 27(12)： 330～336， 442

[12]肖明， 吳季秋，陳秋波，等.基于CA-Markov模型的昌化江流域土地利用動態變化.農業工程學報， 2012, 28(10): 231～238

[13]何丹， 周璟，高偉，等.基于CA-Markov模型的滇池流域土地利用變化動態模擬研究.北京大學學報(自然科學版)， 2014, 50(6): 1095～1105

[14]龔文峰， 袁力，范文義.基于CA-Markov的哈爾濱市土地利用變化及預測.農業工程學報， 2012, 28(14): 216～222

[15]許月卿， 崔麗，孟繁盈.大城市邊緣區土地利用變化與社會經濟發展關系分析——以北京市平谷區為例.中國農業資源與區劃， 2008, 29(4): 16～21

[16]杜萌， 趙冬玲，楊建宇，等.基于元胞自動機復合模型的土地利用演化模擬——以北京市海淀區為例.測繪學報， 2015， 44(S1)： 68～74

[17]Carver，S.J.Integrating multi-criteria evaluation with geographical information systems，International Journal of Geographical Information Systems， 1991, 5(3): 321～339

[18]劉潔， 李宏，馬勇剛.基于CA-Markov模型的中亞典型城市土地利用變化預測分析.水土保持研究， 2014, 21(3): 51～56