江 山，石旭飛，郭常來，馮雨林，孫秀波，孫家全，倪 金

中國地質(zhì)調(diào)查局沈陽地質(zhì)調(diào)查中心，遼寧沈陽 110034

0 引言

土地利用和覆被變化（land use and cover change，LUCC）是全球環(huán)境變化的重要組成部分和驅(qū)動因素之一［1］，研究土地利用變化規(guī)律和發(fā)展趨勢，對于土地資源的合理開發(fā)與利用，促進區(qū)域協(xié)調(diào)可持續(xù)發(fā)展具有重要意義［2-3］.近年來，國內(nèi)外學者已對土地利用變化時空動態(tài)預測模型進行了深入研究，采用的預測模型主要有CLUE-S 模型、SLEUTH、多智能體系統(tǒng)模型（multi-agent system，MAS）、元胞自動機模型（cellular automata，CA）等［4-9］.其中，CA-Markov 模型耦合了元胞自動機和馬爾可夫鏈的各自優(yōu)勢，可以實現(xiàn)對長時間復雜土地利用類型時空演變過程的可視化模擬和預測，并具有較高的精度.國內(nèi)不少學者應用CA-Markov模型在新疆克州、昌化江流域、云南德宏州、三峽庫區(qū)、黃河流域、錫林河流域、喀什市、長沙市等地區(qū)開展土地利用變化分析與預測研究［10-16］.

大凌河是遼寧西部地區(qū)最大的河流，位于大凌河流域的錦州、阜新、朝陽三市是東接遼沈，西聯(lián)京津唐地區(qū)的重要城市.大凌河流域內(nèi)地形復雜，石漠化、水土流失等問題嚴重，其生態(tài)系統(tǒng)穩(wěn)定直接關(guān)系著遼西的經(jīng)濟發(fā)展［17-18］.一些學者在大凌河流域開展的土地利用變化研究有力揭示了流域生態(tài)環(huán)境變化特征［19-20］，只是研究數(shù)據(jù)較陳舊，不能反映流域新變化.鑒于此，本研究基于自然資源部新發(fā)布的土地覆蓋數(shù)據(jù)，從土地利用動態(tài)度、土地利用程度和土地利用轉(zhuǎn)移矩陣等不同角度，分析研究區(qū)過去20 年來來土地利用變化規(guī)律，并基于CA-Markov 模型對研究區(qū)2030 年土地利用格局進行預測，以期為大凌河流域生態(tài)環(huán)境保護與區(qū)域經(jīng)濟和社會發(fā)展提供科學依據(jù).

1 研究區(qū)簡介

研究區(qū)大凌河流域覆蓋的地理范圍為：118°53'—121°52'E、40°28'—42°38'N.河流全長435 km，上游分南、北兩支，于喀左縣大城子東南匯合后流經(jīng)朝陽、北票、凌海、義縣等地，最終匯入渤海（圖1）.流域總面積2.33×104km2，其中2.08×104km2位于遼寧省境內(nèi).地貌以山地丘陵為主，少量平原區(qū).氣候類型屬于中溫帶氣候，四季冷暖干濕分明，溫度變化較大，多年平均氣溫8.3 ℃，平均相對濕度53%，日照時數(shù)2 800 h，年均降雨量465 mm，年蒸發(fā)量1 974.4 mm，年均徑流量17.91×108m3［19］.

圖1 大凌河流域地理位置圖Fig.1 Geographical location map of Daling River Basin

2 數(shù)據(jù)來源及預處理

本研究主要基于自然資源部發(fā)布的30 m 全球地表覆蓋數(shù)據(jù)（GlobeLand30）.該數(shù)據(jù)從2014 年開始發(fā)布，目前已完成2020 年版更新.GlobeLand30 數(shù)據(jù)研制所使用的分類影像主要有美國陸地資源衛(wèi)星（Landsat）、中國環(huán)境減災衛(wèi)星（HJ-1）和高分一號（GF-1）等多光譜影像，優(yōu)選植被生長季少云影像［21］.GlobeLand30 數(shù)據(jù)共包括10 個一級類型，分別是：耕地、林地、草地、灌木地、濕地、水體、苔原、人造地表（建設(shè)用地）、裸地、冰川和永久積雪，研究區(qū)內(nèi)共有其中7類數(shù)據(jù)（如表1 所示）.GlobeLand30 數(shù)據(jù)精度評價由同濟大學和中國科學院空天信息創(chuàng)新研究院等牽頭完成，平均總體精度為84.61%，平均Kappa 系數(shù)為0.80.

表1 土地利用類型分類表Table 1 Classification of land use types

本次研究共收集了6 幅數(shù)據(jù)，數(shù)據(jù)列表為：N50_40_2000LC030、N51_40_2000LC030、N50_40_2010LC030、N51_40_2010LC030、N50_40_2020LC030、N51_40_2020LC030.數(shù)據(jù)處理與分析采用ArcGIS 和IDRISI 軟件完成.

3 研究方法

3.1 土地利用動態(tài)度

土地利用動態(tài)主要用于表征一段時間內(nèi)土地利用類型的變化情況，即土地利用的變化速度，分為單一土地利用動態(tài)度和綜合土地利用動態(tài)度［22-24］.計算公式如下：

式中，L 為研究時段內(nèi)單一土地利用類型的動態(tài)度；Ua和Ub分別表示該類土地利用類型研究始、末的面積；LC 為綜合土地利用動態(tài)度，LUi為第i 類土地類型在研究初期的面積；ΔLUi-j為研究時段內(nèi)第i 類土地利用類型轉(zhuǎn)變?yōu)榉莍 類（j 類）土地利用類型面積的絕對值；T 為研究時段，當該時段設(shè)定為年時，L 為該研究區(qū)某種土地利用類型年變化率，LC 為研究區(qū)土地利用年變化率［23］.

3.2 土地利用程度

土地利用程度是在某個時間內(nèi)由自然因素和人為因素雙重作用的結(jié)果，它反映了土地利用的廣度和深度.根據(jù)劉紀遠等［1］提出的土地利用程度的綜合分析方法，按照土地自然綜合體在社會因素影響下的自然平衡狀態(tài)將土地利用程度分為4 個級別，依次賦予分級指數(shù)（表2），在此基礎(chǔ)上建立土地利用程度綜合指數(shù)及變化模型［25］.

表2 土地利用類型及分級表Table 2 Land use types and grading indexes

1）土地利用程度綜合指數(shù)

土地利用程度綜合指數(shù)計算公式如下：

式中，Li表示土地利用程度綜合指數(shù)；Ai表示第i 級的土地利用程度分級指數(shù)；Ci表示第i 級的土地利用類型面積百分比；n 為土地利用程度分級數(shù).Li數(shù)值位于100～400 區(qū)間內(nèi)，數(shù)值大小反映土地利用程度的高低.

2）土地利用程度變化

土地利用程度的變化是指特定區(qū)域在一定的時段內(nèi)的綜合變化，通常由土地利用程度的變化量表征，計算公式如下：

式中，La和Lb分別為a、b 兩個時間的研究區(qū)土地利用程度綜合指數(shù)；Ai為第i 級的土地利用程度分級指數(shù)；Cia和Cib分別為研究區(qū)a 時間和b 時間第i 級土地類型利用程度面積百分比.如果ΔLb-a＞0，表示研究區(qū)土地利用處于發(fā)展時期；如果ΔLb-a＜0，表示該區(qū)域的土地利用處于衰退期.

3.3 土地利用轉(zhuǎn)移矩陣

土地利用變化不僅是面積數(shù)量上的增減，而且還存在類型之間的轉(zhuǎn)換關(guān)系.土地利用轉(zhuǎn)移矩陣能夠定量描述各種土地利用類型之間的變化，有助于了解期初各地類的流失去向和期末各地類的來源和構(gòu)成.通過對任意兩期土地利用數(shù)據(jù)進行空間疊加，獲取各研究時段土地利用類型轉(zhuǎn)移矩陣，其表達式為：

式中，S 為面積；n 為土地利用類型數(shù)量；i、j 分別為研究初期與末期的土地利用類型.

3.4 CA-Markov 模型

CA-Markov 模型結(jié)合了元胞自動機模型和馬爾可夫鏈的優(yōu)勢，常用于土地利用時空格局的模擬預測，具有高精度、可視化等優(yōu)點［26］.Markov 鏈描述的是土地利用從一個時刻（t）到另一個時刻（t＋1）的變化趨勢，并以此為規(guī)則預測將來的變化趨勢［26-27］.公式如下：

式中，St、St＋1分別為t、t＋1 時刻的系統(tǒng)狀態(tài)，Pij為狀態(tài)轉(zhuǎn)移概率.

CA 模型同時具備時間和空間預測功能，模型變量包括狀態(tài)、領(lǐng)域空間和時間等離散變量，公式如下：

式中，S 為元胞有限、離散的狀態(tài)集合；t 和t＋1 表示不同時刻；N 為元胞的鄰域；f 為局部空間元胞轉(zhuǎn)換規(guī)則［28］.

4 結(jié)果分析

4.1 土地利用空間格局

研究區(qū)土地利用以耕地和草地為主.其中耕地占總面積的50%左右，主要集中分布在大凌河中、下游，少量位于上游丘陵地區(qū)谷地；草地占比超過30%，主要分布在上游丘陵區(qū)；林地占比12%左右；建設(shè)用地約占比5%；濕地和水體面積較小，約占總面積的1%；裸地面積可忽略不計（圖2，表3），表明研究區(qū)土地開發(fā)程度處于較高水平.

表3 大凌河流域土地利用面積統(tǒng)計表Table 3 Statistics of land use area in Daling River Basin

圖2 大凌河流域2000、2010、2020 年土地利用空間格局圖Fig.2 Land use spatial patterns of Daling River Basin in 2000，2010 and 2020

4.2 土地利用變化分析

1）土地利用結(jié)構(gòu)變化

2000、2010、2020 年研究區(qū)土地利用結(jié)構(gòu)有較大變化.2000—2020 年期間，研究區(qū)內(nèi)耕地持續(xù)減少，其中前10 年減少了132.85 km2，后10 年減少了1 360.13 km2，減少速度在加快.近20 年間，耕地共減少了1 492.97 km2，面積占比由2000 年的53.23%減少為2020 年的46.02%.與此同時，林地、草地和建設(shè)用地面積不斷增加，其中，草地增加了616.17 km2，建設(shè)用地增加492.23 km2，林地共增加346.97 km2（表3）.

2）土地利用動態(tài)度

從土地利用動態(tài)度計算結(jié)果（表4）來看，2000—2010 年間，研究區(qū)綜合土地利用動態(tài)度為43.53%，單地類中動態(tài)度最大的是濕地，平均每年增加86.32%，動態(tài)度最小的是耕地，平均每年減少0.12%.2010—2020 年間，研究區(qū)綜合土地利用動態(tài)度為10.36%，單地類中動態(tài)度最大的是水體，平均每年增加5.00%，動態(tài)度最小的是濕地，平均每年減少7.22%.總體來看，近20 年間，研究區(qū)綜合土地利用動態(tài)度為15.60%，處于不斷減少的狀態(tài).單地類中動態(tài)度最大的是濕地，平均每年增加8.38%，動態(tài)度最小的是耕地，平均每年減少0.68%.

表4 大凌河流域土地利用動態(tài)度統(tǒng)計表Table 4 Dynamic degree statistics of land use in Daling River Basin

3）土地利用程度

根據(jù)式（3）（4）計算了研究區(qū)2000、2010 和2020年土地利用綜合指數(shù)，從數(shù)據(jù)（表5）來看，研究區(qū)平均土地利用程度綜合指數(shù)為261.5，土地利用程度處于中等偏高水平.由2000 年的262.53 變?yōu)?020 年259.91，土地利用變化量在持續(xù)降低，表明大凌河流域近20 年土地處于衰退期，且后半期（2010—2020 年）衰退速度大于前半期（2000—2010）.

表5 大凌河流域土地利用綜合指數(shù)統(tǒng)計表Table 5 Comprehensive land use index of Daling River Basin

4）土地利用變化矩陣

通過建立土地利用轉(zhuǎn)移矩陣（圖3，表6），可知2000—2010 年間，耕地減少了132.85 km2，主要轉(zhuǎn)為草地、林地和建設(shè)用地.其中轉(zhuǎn)為草地和林地的耕地面積分別為130.74 和54.96 km2，主要分布于大凌河下游義縣-錦州一帶；共有33.34 km2耕地轉(zhuǎn)為建設(shè)用地，主要分布于朝陽市和錦州市等地.2010—2020 年間，耕地減少面積達1 360.13 km2，其中有723.68 km2耕地轉(zhuǎn)為草地，主要發(fā)生在大凌河上游和下游；共有403.15 km2耕地轉(zhuǎn)為林地，主要位于大凌河上游地區(qū)；另有473.15 km2耕地轉(zhuǎn)為建設(shè)用地，集中分布在阜新市一帶.與此同時，有少量的草地和林地轉(zhuǎn)化為耕地，面積分別為143.18 和77.74 km2（如圖3，表7）.

表6 大凌河流域2000—2010 年土地利用變化矩陣Table 6 Land use change matrix of Daling River Basin during 2000-2010

表7 大凌河流域2010—2020 年土地利用變化矩陣Table 7 Land use change matrix of Daling River Basin during 2010-2020

圖3 大凌河流域土地利用分類變化圖Fig.3 Land use variation map of Daling River Basin

4.3 土地利用預測分析

本研究利用IDRISI 軟件進行土地利用的模擬預測分析.第一步在Markov 模塊設(shè)置相對誤差為0.15，分別計算研究區(qū)2000—2010 年、2010—2020 年的土地利用轉(zhuǎn)移概率矩陣；第二步在CA-Markov 模塊以2000—2010 年的土地利用轉(zhuǎn)移概率矩陣為條件，設(shè)置30 m×30 m 元胞、5×5 鄰域結(jié)構(gòu)、元胞自動機循環(huán)次數(shù)為10 等參數(shù)，模擬預測2020 年研究區(qū)土地利用分布圖；第三步在CROSSTAB 模塊對2020 年預測值和實際值進行交叉驗證，得到Kappa 系數(shù)為0.8835，滿足精度要求.然后重復第二步，以2020 年實際土地利用數(shù)據(jù)為基礎(chǔ)，模擬預測研究區(qū)2030 年土地利用格局［29］（如圖4a 所示）.

根據(jù)預測結(jié)果，建立2020—2030 年土地利用轉(zhuǎn)移矩陣（圖4b，表8）.結(jié)果表明2030 年研究區(qū)耕地將進一步減少1 699.92 km2，其中有630.42 km2耕地轉(zhuǎn)為草地，主要位于大凌河上游和下游；共有571.65 km2耕地轉(zhuǎn)為建設(shè)用地，主要發(fā)生在朝陽市和義縣-錦州一帶；另有454.19 km2轉(zhuǎn)為林地，主要發(fā)生在大凌河中上游地區(qū).

表8 大凌河流域2020—2030 年土地利用變化矩陣Table 8 Land use change matrix of Daling River Basin during 2020-2030

5 結(jié)論

通過分析和預測大凌河流域土地利用格局，得出如下結(jié)論：

1）2000—2020 年大凌河流域土地利用以耕地和草地為主，二者占比超過80%.近20 年間，研究區(qū)內(nèi)林地、草地和建設(shè)用地面積不斷增加，耕地面積持續(xù)減少，且減少的速度在加快，前后兩個10 年分別減少了132.85 km2、1 360.13 km2，耕地面積占比由2000 年的53.23%，減少為2020 年的46.02%.林地和草地增加主要位于大凌河上游和下游，建設(shè)用地增加主要位于朝陽市和錦州市周邊.這些變化表明區(qū)域內(nèi)“退耕還林還草”工程實施效果明顯，建設(shè)用地的擴張規(guī)?；究煽?

2）從土地利用動態(tài)度分析來看，近20 年來，大凌河流域綜合土地利用動態(tài)度從2000—2010 年期間的43.53%，下降到2010—2020 年期間為10.36%，處于不斷減少的狀態(tài).單地類來看，濕地動態(tài)度最大，平均每年增加8.38%.表明區(qū)域土地利用方式趨于穩(wěn)定，濕地等生態(tài)系統(tǒng)正在恢復.

3）2000—2020 年間，大凌河流域土地利用程度綜合指數(shù)由2000 年的262.53 變?yōu)?020 年259.91，表明區(qū)域土地處于衰退期，土地開發(fā)利用程度在逐步降低，自然生態(tài)系統(tǒng)正在恢復之中.

4）基于2010—2020 年土地利用變化概率矩陣，利用CA-Markov 模型預測了大凌河流域在自然條件下2030 年土地利用格局，模擬精度Kappa 系數(shù)為0.8835，說明結(jié)果可靠.統(tǒng)計表明，2030 年，大凌河流域耕地面積將進一步減少1 699.92 km2，主要變?yōu)椴莸?、建設(shè)用地和林地.當然，土地利用變化預測是一個系統(tǒng)工程，除了歷史土地變化規(guī)律外，還受到人口、經(jīng)濟、環(huán)境和土地調(diào)控政策等多因素的影響，今后研究中需綜合考慮其他因素，進一步優(yōu)化預測模型.