江 山,石旭飞,郭常来,冯雨林,孙秀波,孙家全,倪 金
中国地质调查局沈阳地质调查中心,辽宁沈阳 110034
土地利用和覆被变化(land use and cover change,LUCC)是全球环境变化的重要组成部分和驱动因素之一[1],研究土地利用变化规律和发展趋势,对于土地资源的合理开发与利用,促进区域协调可持续发展具有重要意义[2-3].近年来,国内外学者已对土地利用变化时空动态预测模型进行了深入研究,采用的预测模型主要有CLUE-S 模型、SLEUTH、多智能体系统模型(multi-agent system,MAS)、元胞自动机模型(cellular automata,CA)等[4-9].其中,CA-Markov 模型耦合了元胞自动机和马尔可夫链的各自优势,可以实现对长时间复杂土地利用类型时空演变过程的可视化模拟和预测,并具有较高的精度.国内不少学者应用CA-Markov模型在新疆克州、昌化江流域、云南德宏州、三峡库区、黄河流域、锡林河流域、喀什市、长沙市等地区开展土地利用变化分析与预测研究[10-16].
大凌河是辽宁西部地区最大的河流,位于大凌河流域的锦州、阜新、朝阳三市是东接辽沈,西联京津唐地区的重要城市.大凌河流域内地形复杂,石漠化、水土流失等问题严重,其生态系统稳定直接关系着辽西的经济发展[17-18].一些学者在大凌河流域开展的土地利用变化研究有力揭示了流域生态环境变化特征[19-20],只是研究数据较陈旧,不能反映流域新变化.鉴于此,本研究基于自然资源部新发布的土地覆盖数据,从土地利用动态度、土地利用程度和土地利用转移矩阵等不同角度,分析研究区过去20 年来来土地利用变化规律,并基于CA-Markov 模型对研究区2030 年土地利用格局进行预测,以期为大凌河流域生态环境保护与区域经济和社会发展提供科学依据.
研究区大凌河流域覆盖的地理范围为:118°53'—121°52'E、40°28'—42°38'N.河流全长435 km,上游分南、北两支,于喀左县大城子东南汇合后流经朝阳、北票、凌海、义县等地,最终汇入渤海(图1).流域总面积2.33×104km2,其中2.08×104km2位于辽宁省境内.地貌以山地丘陵为主,少量平原区.气候类型属于中温带气候,四季冷暖干湿分明,温度变化较大,多年平均气温8.3 ℃,平均相对湿度53%,日照时数2 800 h,年均降雨量465 mm,年蒸发量1 974.4 mm,年均径流量17.91×108m3[19].
图1 大凌河流域地理位置图Fig.1 Geographical location map of Daling River Basin
本研究主要基于自然资源部发布的30 m 全球地表覆盖数据(GlobeLand30).该数据从2014 年开始发布,目前已完成2020 年版更新.GlobeLand30 数据研制所使用的分类影像主要有美国陆地资源卫星(Landsat)、中国环境减灾卫星(HJ-1)和高分一号(GF-1)等多光谱影像,优选植被生长季少云影像[21].GlobeLand30 数据共包括10 个一级类型,分别是:耕地、林地、草地、灌木地、湿地、水体、苔原、人造地表(建设用地)、裸地、冰川和永久积雪,研究区内共有其中7类数据(如表1 所示).GlobeLand30 数据精度评价由同济大学和中国科学院空天信息创新研究院等牵头完成,平均总体精度为84.61%,平均Kappa 系数为0.80.
表1 土地利用类型分类表Table 1 Classification of land use types
本次研究共收集了6 幅数据,数据列表为:N50_40_2000LC030、N51_40_2000LC030、N50_40_2010LC030、N51_40_2010LC030、N50_40_2020LC030、N51_40_2020LC030.数据处理与分析采用ArcGIS 和IDRISI 软件完成.
土地利用动态主要用于表征一段时间内土地利用类型的变化情况,即土地利用的变化速度,分为单一土地利用动态度和综合土地利用动态度[22-24].计算公式如下:
式中,L 为研究时段内单一土地利用类型的动态度;Ua和Ub分别表示该类土地利用类型研究始、末的面积;LC 为综合土地利用动态度,LUi为第i 类土地类型在研究初期的面积;ΔLUi-j为研究时段内第i 类土地利用类型转变为非i 类(j 类)土地利用类型面积的绝对值;T 为研究时段,当该时段设定为年时,L 为该研究区某种土地利用类型年变化率,LC 为研究区土地利用年变化率[23].
土地利用程度是在某个时间内由自然因素和人为因素双重作用的结果,它反映了土地利用的广度和深度.根据刘纪远等[1]提出的土地利用程度的综合分析方法,按照土地自然综合体在社会因素影响下的自然平衡状态将土地利用程度分为4 个级别,依次赋予分级指数(表2),在此基础上建立土地利用程度综合指数及变化模型[25].
表2 土地利用类型及分级表Table 2 Land use types and grading indexes
1)土地利用程度综合指数
土地利用程度综合指数计算公式如下:
式中,Li表示土地利用程度综合指数;Ai表示第i 级的土地利用程度分级指数;Ci表示第i 级的土地利用类型面积百分比;n 为土地利用程度分级数.Li数值位于100~400 区间内,数值大小反映土地利用程度的高低.
2)土地利用程度变化
土地利用程度的变化是指特定区域在一定的时段内的综合变化,通常由土地利用程度的变化量表征,计算公式如下:
式中,La和Lb分别为a、b 两个时间的研究区土地利用程度综合指数;Ai为第i 级的土地利用程度分级指数;Cia和Cib分别为研究区a 时间和b 时间第i 级土地类型利用程度面积百分比.如果ΔLb-a>0,表示研究区土地利用处于发展时期;如果ΔLb-a<0,表示该区域的土地利用处于衰退期.
土地利用变化不仅是面积数量上的增减,而且还存在类型之间的转换关系.土地利用转移矩阵能够定量描述各种土地利用类型之间的变化,有助于了解期初各地类的流失去向和期末各地类的来源和构成.通过对任意两期土地利用数据进行空间叠加,获取各研究时段土地利用类型转移矩阵,其表达式为:
式中,S 为面积;n 为土地利用类型数量;i、j 分别为研究初期与末期的土地利用类型.
CA-Markov 模型结合了元胞自动机模型和马尔可夫链的优势,常用于土地利用时空格局的模拟预测,具有高精度、可视化等优点[26].Markov 链描述的是土地利用从一个时刻(t)到另一个时刻(t+1)的变化趋势,并以此为规则预测将来的变化趋势[26-27].公式如下:
式中,St、St+1分别为t、t+1 时刻的系统状态,Pij为状态转移概率.
CA 模型同时具备时间和空间预测功能,模型变量包括状态、领域空间和时间等离散变量,公式如下:
式中,S 为元胞有限、离散的状态集合;t 和t+1 表示不同时刻;N 为元胞的邻域;f 为局部空间元胞转换规则[28].
研究区土地利用以耕地和草地为主.其中耕地占总面积的50%左右,主要集中分布在大凌河中、下游,少量位于上游丘陵地区谷地;草地占比超过30%,主要分布在上游丘陵区;林地占比12%左右;建设用地约占比5%;湿地和水体面积较小,约占总面积的1%;裸地面积可忽略不计(图2,表3),表明研究区土地开发程度处于较高水平.
表3 大凌河流域土地利用面积统计表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
1)土地利用结构变化
2000、2010、2020 年研究区土地利用结构有较大变化.2000—2020 年期间,研究区内耕地持续减少,其中前10 年减少了132.85 km2,后10 年减少了1 360.13 km2,减少速度在加快.近20 年间,耕地共减少了1 492.97 km2,面积占比由2000 年的53.23%减少为2020 年的46.02%.与此同时,林地、草地和建设用地面积不断增加,其中,草地增加了616.17 km2,建设用地增加492.23 km2,林地共增加346.97 km2(表3).
2)土地利用动态度
从土地利用动态度计算结果(表4)来看,2000—2010 年间,研究区综合土地利用动态度为43.53%,单地类中动态度最大的是湿地,平均每年增加86.32%,动态度最小的是耕地,平均每年减少0.12%.2010—2020 年间,研究区综合土地利用动态度为10.36%,单地类中动态度最大的是水体,平均每年增加5.00%,动态度最小的是湿地,平均每年减少7.22%.总体来看,近20 年间,研究区综合土地利用动态度为15.60%,处于不断减少的状态.单地类中动态度最大的是湿地,平均每年增加8.38%,动态度最小的是耕地,平均每年减少0.68%.
表4 大凌河流域土地利用动态度统计表Table 4 Dynamic degree statistics of land use in Daling River Basin
3)土地利用程度
根据式(3)(4)计算了研究区2000、2010 和2020年土地利用综合指数,从数据(表5)来看,研究区平均土地利用程度综合指数为261.5,土地利用程度处于中等偏高水平.由2000 年的262.53 变为2020 年259.91,土地利用变化量在持续降低,表明大凌河流域近20 年土地处于衰退期,且后半期(2010—2020 年)衰退速度大于前半期(2000—2010).
表5 大凌河流域土地利用综合指数统计表Table 5 Comprehensive land use index of Daling River Basin
4)土地利用变化矩阵
通过建立土地利用转移矩阵(图3,表6),可知2000—2010 年间,耕地减少了132.85 km2,主要转为草地、林地和建设用地.其中转为草地和林地的耕地面积分别为130.74 和54.96 km2,主要分布于大凌河下游义县-锦州一带;共有33.34 km2耕地转为建设用地,主要分布于朝阳市和锦州市等地.2010—2020 年间,耕地减少面积达1 360.13 km2,其中有723.68 km2耕地转为草地,主要发生在大凌河上游和下游;共有403.15 km2耕地转为林地,主要位于大凌河上游地区;另有473.15 km2耕地转为建设用地,集中分布在阜新市一带.与此同时,有少量的草地和林地转化为耕地,面积分别为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
本研究利用IDRISI 软件进行土地利用的模拟预测分析.第一步在Markov 模块设置相对误差为0.15,分别计算研究区2000—2010 年、2010—2020 年的土地利用转移概率矩阵;第二步在CA-Markov 模块以2000—2010 年的土地利用转移概率矩阵为条件,设置30 m×30 m 元胞、5×5 邻域结构、元胞自动机循环次数为10 等参数,模拟预测2020 年研究区土地利用分布图;第三步在CROSSTAB 模块对2020 年预测值和实际值进行交叉验证,得到Kappa 系数为0.8835,满足精度要求.然后重复第二步,以2020 年实际土地利用数据为基础,模拟预测研究区2030 年土地利用格局[29](如图4a 所示).
根据预测结果,建立2020—2030 年土地利用转移矩阵(图4b,表8).结果表明2030 年研究区耕地将进一步减少1 699.92 km2,其中有630.42 km2耕地转为草地,主要位于大凌河上游和下游;共有571.65 km2耕地转为建设用地,主要发生在朝阳市和义县-锦州一带;另有454.19 km2转为林地,主要发生在大凌河中上游地区.
表8 大凌河流域2020—2030 年土地利用变化矩阵Table 8 Land use change matrix of Daling River Basin during 2020-2030
通过分析和预测大凌河流域土地利用格局,得出如下结论:
1)2000—2020 年大凌河流域土地利用以耕地和草地为主,二者占比超过80%.近20 年间,研究区内林地、草地和建设用地面积不断增加,耕地面积持续减少,且减少的速度在加快,前后两个10 年分别减少了132.85 km2、1 360.13 km2,耕地面积占比由2000 年的53.23%,减少为2020 年的46.02%.林地和草地增加主要位于大凌河上游和下游,建设用地增加主要位于朝阳市和锦州市周边.这些变化表明区域内“退耕还林还草”工程实施效果明显,建设用地的扩张规模基本可控.
2)从土地利用动态度分析来看,近20 年来,大凌河流域综合土地利用动态度从2000—2010 年期间的43.53%,下降到2010—2020 年期间为10.36%,处于不断减少的状态.单地类来看,湿地动态度最大,平均每年增加8.38%.表明区域土地利用方式趋于稳定,湿地等生态系统正在恢复.
3)2000—2020 年间,大凌河流域土地利用程度综合指数由2000 年的262.53 变为2020 年259.91,表明区域土地处于衰退期,土地开发利用程度在逐步降低,自然生态系统正在恢复之中.
4)基于2010—2020 年土地利用变化概率矩阵,利用CA-Markov 模型预测了大凌河流域在自然条件下2030 年土地利用格局,模拟精度Kappa 系数为0.8835,说明结果可靠.统计表明,2030 年,大凌河流域耕地面积将进一步减少1 699.92 km2,主要变为草地、建设用地和林地.当然,土地利用变化预测是一个系统工程,除了历史土地变化规律外,还受到人口、经济、环境和土地调控政策等多因素的影响,今后研究中需综合考虑其他因素,进一步优化预测模型.