基于RFFLUS-InVEST-Geodetector耦合模型的平朔矿区生境质量时空演变及其影响因素
2021-11-24杨国婷郭东罡张霄羽
杨国婷，张 红，李 静，郭东罡，张霄羽
(1 山西大学 环境与资源学院，山西 太原 030000;2 中国科学院 地理科学与资源研究所，北京 100101)
生境质量是指生态系统保障生物适宜生存的能力，常被用来反映地域生物多样性状况[1]。人类对矿山进行露天开采，会导致生境的破损、退化与消失。土地是地表各类生态系统的载体[2]，土地利用变化是人类活动对生态环境作用的直接体现, 也是引起区域生境质量变化的重要因素[3]。因此，从土地利用的角度出发，评估并预测生境质量的时空演变规律，探究影响生境质量的自然和人为因素，可为区域生态保护提供参考。
目前，国内外生境质量评估方法主要有加权评价法[4]、层次分析法[5]、生态遥感指数法[6]等。随着研究的深入与技术的发展，生境质量的评估更趋向于定量化、空间化和精细化[7]。InVEST(integrated valuation of ecosystem services and trade-offs, InVEST)模型因具有数据量少、输入简单、评估结果精确、支持不同情景和不同尺度的评估及分析等优点，在国内外研究中被广泛应用[8]。众多学者利用InVEST模型对流域、自然保护区和城市等不同区域进行生境质量评估[9-11]，但对土地利用模拟预测、生境质量评估以及影响因素的综合分析鲜见报道，对矿区生境质量的研究也很少见。
土地利用变化是InVEST模型的重要参数，精确地对其进行模拟和预测，对于分析未来生境质量演变趋势十分重要。随机森林和元胞自动机相结合的土地利用模型RFFLUS(random forest future land use simulation，RFFLUS)，是在FLUS(future land use simulation)模型的基础上，用RF模型取代FLUS模型的ANN(artificial neural network，ANN)运行环节，可以更精准地模拟区域土地利用变化，具有数据可视化、精度高等优点。地理探测器(Geodetector)与地理加权回归、空间自相关指数等方法相比，具有交互探测的独特优势[12]。三个模型均适用于区域尺度的研究，RFFLUS模型可以精准预测土地利用变化并为生境质量的预测提供数据基础，InVEST模型能定量描述土地利用与生境质量的关系并实现生境质量的空间分布，Geodetector模型可以直接构建自然与社会驱动因子与生境质量之间的联系并明晰单因子和双因子对生境质量的作用程度。因此，将RFFLUS模型、InVEST模型和Geodetector模型相耦合，可以更好地研究区域生境质量的时空演变及归因分析。
矿区的土地利用变化剧烈，呈现阶段性变化特征，尤其是露天开采对区域生态系统造成严重破坏[13]。平朔矿区是我国目前规模最大、现代化程度最高的大型露采基地，历经30余年的开采，生境质量受到严重的损坏。本文基于长时间序列遥感影像，结合矿区环境特征、开采历史、社会经济指标，耦合RFFLUS模型、InVEST模型和Geodetector模型，对平朔矿区近30年土地利用变化和生境质量的时空演变及未来趋势进行分析、模拟和预测，探索生境质量时空演变规律，揭示影响生境质量的人为与自然驱动力因素，对矿区生态恢复实践和土地利用规划具有现实指导意义。
1 研究区概括与研究方法
1.1 研究区概况
平朔矿区是我国五大特大型露天煤矿之一，区内有安太堡、安家岭、东露天三大矿区，也是国家规划建设的晋北煤炭基地的重要组成部分。该矿区位于山西省北部、朔州市中西部，地理坐标为东经112°14′～112°32′、北纬39°22′～39°37′，总面积约442 km2(图1)；属温带半干旱大陆性季风气候区，年平均气温4.8～9.6℃，年平均降水量461.12～476.74 mm，为黄土丘陵地貌，地表水系比较丰富，以栗钙土为主，植被主要呈农业耕作景观，煤炭储量多[14]；主要含煤层为石炭系和二叠系，各可采煤层为气煤、肥煤，是良好的动力煤[15]。
图1 平朔矿区位置及高程图示意图Fig.1 Location and DEM of Pingshuo mining area
平朔矿区自1985年建矿以来，共经历了4个阶段：第1阶段为起步阶段(1985—1995年)，1987年安太堡露天矿建成投产，核定原煤生产能力为2 200万t/a，资源回收率达到95%以上[16]；第2阶段为发展阶段(1995—2005年)，2003年安家岭露天矿正式投产，建设用地进一步扩大；第3阶段为跨越阶段(2005—2015年)，2011年东露天煤矿正式投产，国家规划区的3座露天矿全部建成投产；第4阶段为转型阶段(2015年至今)，矿区进行土地复垦和生态重建工作，促进矿区生态可持续发展。
1.2 数据来源
本文使用的数据主要分为四大类(表1)：(1)InVEST模型。威胁源数据是描述威胁分布和强度的栅格数据，根据InVEST模型要求并参照以往研究经验，自然环境对外界威胁的影响更敏感，人工环境对威胁的影响相对更弱，本文将耕地、建设用地和其他用地设为威胁源。InVEST模型中每种地类的生境得分及其对威胁源的敏感度参数、威胁源权重和威胁源最大影响距离参数通过参考文献[17-20]并结合实际测定确定。(2)土地利用数据。土地利用数据为1985年、1995年、2005年、2015年的Landsat遥感影像数据，数据来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)，分辨率为30 m×30 m。通过人机交互的解译方式获取，并按照国家标准划分为六大类。(3)自然环境数据。包括高程、气温、降水数据，高程数据来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)，精度为30 m；NDVI由遥感影像计算获取；坡度、坡向数据从高程数据中提取；气温、降水数据来源于中国气象数据网(http://data.cma.cn)，精度为30 m。距离因子包括距公路距离、距铁路距离、距水系距离、距设施点距离、距居民点距离，距离因子数据由遥感影像通过欧氏距离处理得到。(4)社会经济数据。包括人口、GDP、煤炭生产能力数据，人口和GDP数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)，GDP数据和人口数据分辨率通过重采样提升至30 m，煤炭生产能力数据通过实地调查获取。
表1 数据来源与用途Tab.1 Sources and uses of data
1.3 研究方法
本文将RFFLUS模型、InVEST模型、Geodetector模型进行耦合(图2)，综合利用各模型的优势，预测生境质量时空演变规律并对生境质量的驱动力机制进行探测。在分析平朔矿区1985—2015年土地利用变化的基础上，利用RFFLUS模型精确地预测平朔矿区2025年土地利用变化，采用InVEST模型开展平朔矿区1985—2025年生境质量时空演变的探究，最后用Geodetector模型对平朔矿区生境质量影响因素进行归因分析。
图2 RFFLUS-InVEST-Geodetector耦合模型Fig.2 RFFLUS-InVEST-Geodetector coupling model
1.3.1 InVEST模型 InVEST模型的生境质量模块将输入的威胁源数据进行计算获取生境退化度，通过生境的适宜性数据和退化度得出生境质量[17]。生境退化程度计算公式为
(1)
式中：Dxj表示位于栅格x地类j的生境退化度；r为威胁因子；R为威胁因子的总数；y为威胁因子所在的栅格；Yr为威胁因子r的数量；Wr为威胁因子的权重；ry为位于栅格y威胁因子的值；βx为栅格x的可达性水平；Sjr为地类j受威胁因子r的敏感性；irxy为威胁因子ry对栅格x的影响。其中irxy有以下两种情况：
irxy=1-(dxy/drmax)(线性衰退),
(2)
irxy=exp[(-2.99/drmax)dxy](指数衰退)。
(3)
式中：dxy为x与y两个栅格的线性距离；drmax为威胁因子r的最大作用距离。
假定模型中，土地利用类型对威胁越敏感，其越容易退化。在得到生境退化程度后，利用半饱和函数将栅格退化分值解译成生境质量得分值，公式为
(4)
式中：Qxj指位于栅格x地类j的生境质量；Hj表示地类j的生境适宜程度；Dxj表示位于栅格x地类j的生境退化度；K为半饱和参数；z为尺度常数，一般取2.5。
该模块中需要输入的参数主要有矿区土地利用类型图、威胁源种类、土地利用类型对威胁因子的敏感度(表2)，威胁源权重和最大影响距离(表3)等。本研究在参考相关文献的基础上，结合研究区实际情况确定生境质量模块的输入参数。
表2 生境适宜性及各生境类型对威胁因子的相对敏感度Tab.2 Habitat suitability and relative sensitivity of habitat types to threat factors
表3 威胁源权重和最大影响距离Tab.3 Threats weight and maximum influencing distance
1.3.2 RFFLUS模型 土地利用变化的拟合与预测采用RFFLUS模型，该模型首先采用随机森林训练模块得到土地的初始适宜性概率；其次设计一种自适应惯性和竞争机制，解决不同土地利用类型之间的竞争和相互作用，得到土地利用类型的综合概率[21]；最后使用轮盘赌选择机制，确定何种土地利用类型占据网格单元，在达到土地利用需求下，得到模拟结果。
RF模型是Breiman在2001年提出的一种基于分类树的机器学习算法[22]，该模型可以将数据可视化并提高运行速度。本研究中的决策树个数为24，变量数设置为4，单个样本数据的2/3用于模型训练，其余1/3作为袋外数据(out of bag，OOB)用于验证每棵决策树的精度，OOB越小，精度越高[23]。其概率公式为
(5)
式中：P(p,t,k)为元胞p在t时刻转变为地类k的概率；θ为需要进行分类的数据集；hp(θ)表示元胞p的决策树运算结果；Yk为地类k的决策树分类结果；M为地类数目；I(*)代表分类结果指标函数[23]。
假定研究期间的土地利用规则不变，利用随机森林模型输入的空间因素计算各地类的初始适宜性概率。但是一个土地单元转变为特定的土地利用类型，不仅取决于其发生的概率还取决于预测期内不同发展状态的其他可变因素[21]。RFFLUS模型将初始适宜性概率与自适应惯性竞争机制、转移成本、领域效应相结合，计算综合适宜性概率。
自适应惯性竞争机制的核心思想是通过自适应惯性系数调整土地利用的发展趋势与宏观需求之间的矛盾[21]。自适应惯性系数定义为
(6)
领域效应定义为
(7)
综合适宜性概率为
(8)
1.3.3 地理探测器 地理探测器是利用空间方差分析来探测空间分异性，并揭示驱动因子的空间统计方法。该模型有4个模块，分别为因子探测模块、交互探测模块、风险探测模块与生态探测模块[12]，被广泛应用于生态、环境等领域。本文探测结果以2015年为例，以生境质量变化值为因变量，以自然和人为影响因素为自变量，通过Geodetector模型中的“因子探测器”和“交互作用探测器”两个模块，定量分析各影响因子解释力和影响因子之间的交互类型与大小。
(1)因子探测器模块通过比较子区域方差与总区域方差的关系，明确各个自变量对因变量的解释力[16]。q值是对自变量解释力的度量，q的值域为[0,1]，q值越大说明驱动因素对生境质量的变化解释力越强，反之越小。
(9)
(2)交互作用探测器模块通过比较两个影响因子对生境质量的q(X1)值和q(X2)值与两个影响因子交互q(X1∩X2)值之间的关系判断各自变量对生境质量的共同作用程度。若q(X1∩X2)
2 土地利用变化
2.1 土地利用演变分析
从平朔矿区的土地利用时空演变与各年份土地类型面积变化(图3)得知，1985—2015年平朔矿区土地利用类型以耕地、草地和建设用地为主，其占总面积的90%以上。从变化趋势分析，近30年来耕地和水域分别减少25.19%和0.74%；建设用地、草地、林地和其他用地分别增长16.23%、8.27%、0.56%和0.87%。
图3 1985—2015年平朔矿区土地利用时空演变示意图Fig.3 Spatial-temporal evolution of land use in Pingshuo mining area from 1985 to 2015
1985—1995年平朔矿区处于起步阶段，仅有安太堡矿区开采，开采强度较小。其中：耕地占矿区总面积的一半以上，增加527.94 hm2；草地广泛分布于整个研究区，减少751.23 hm2；林地主要分布在西北部，减少微弱，变化率为—0.12%；建设用地分布于中部，增幅较小，增加了351.63 hm2；水域面积极小，呈条状分布在西北部边缘，中部和南部零散分布，减少73.80hm2。
1995—2005年平朔矿区处于发展阶段，建设用地大幅度增加，占比由6.63%增至13.36%，这是由于安家岭煤矿的投产扩大了建设用地。林地面积增加1.22%、草地面积增加14.08%，新增其他用地191.25 hm2，约占区域面积0.43%；耕地大面积减少，由55.01%下降到33.14%，耕地复垦为林地、草地。
2005—2015年平朔矿区处于跨越阶段，建设用地不断扩大，其余地类均减小。建设用地增加8.70%，这是因为东露天煤矿的投产与扩建。林地减少0.54%，草地减少4.11%，耕地减少4.52%，矿区在林地、草地、耕地上进行了复垦。三十年来，露天开采使土地利用格局发生了巨大改变。
从平朔矿区土地利用转移矩阵(表4)可以看出，1985—2015年平朔矿区大量林地、草地、耕地转化为建设用地，转化面积分别为1 009.62、3 033.72、4 589.19 hm2，这表明近30年平朔矿区不断扩建，大量建设用地也向林地、草地和耕地转移，体现了平朔矿区“边开采边治理”的策略。林地、草地和耕地三者之间相互转化，其中耕地转为林地1 116.63 hm2，转为草地9 680.76 hm2，这说明矿区的复垦方向为林地、草地。
表4 1985—2015年平朔矿区土地利用转移矩阵Tab.4 Change of land use transfer in Pingshuo mining area from 1985 to 2015 单位：hm2
2.2 土地利用预测
本文采用RFFLUS模型预测平朔矿区2025年的土地利用状况，首先以2005年平朔矿区的土地利用数据为基础，运行一个10年的步长，模拟2015年平朔矿区的土地利用分布，并将模拟结果与2015年实际解译结果进行对比验证，kappa系数与总体精度均在0.9以上，证明利用RFFLUS模型可以精准预测未来土地利用分布，结果可靠。
图4表明，研究区2015—2025年土地利用结构发生显著变化，建设用地在矿区周边继续扩大，增幅明显，面积增加4 054.41 hm2，主要分布于研究区中部、西北部和东北部。草地、耕地显著减少，面积分别减少2 215.44、1 754.19 hm2，减少的原因主要是草地和耕地被附近的建设用地进一步挤压占用。林地和水域减少较小，变化幅度不大，在未来几年内保持稳定。
图4 平朔矿区2025年土地利用预测Fig.4 Land use forecast of Pingshuo mining area in 2025
3 平朔矿区生境质量时空演变分析
从空间格局上看(图5)，生境质量等级分布与土地利用分布具有一致性[24]。平朔矿区生境质量空间分布格局整体呈现由中心向外围逐渐增高的态势。低等级生境质量以建设用地为主，分布于研究区西北部、中部和东北部，采矿活动对生境破坏强烈，生境质量值低。较低等级生境质量分布最为广泛，在整个矿区夹杂分布，这些区域主要为耕地和部分草地，人类活动较多，人为干扰对生境质量造成破坏，生境质量处于偏低水平。中等生境质量分布较广，草地居多，人为干扰小，生态系统稳定，生境质量高。较高和高等生境质量以林地为主，零散分布于西北部，林地生境适宜性好，生境质量高。
图5 1985—2025年平朔矿区生境质量等级图Fig.5 Habitat quality grade of Pingshuo mining area from 1985 to 2025
从时间尺度上看(图6),1985年、1995年、2005年、2015年、2025年平均生境质量分别为0.48、0.47、0.53、0.48、0.42，生境质量呈现“降—增—降”趋势。平朔矿区开采初期生境质量降低；随着生态保护意识与复垦经验的提高，生境质量增加；三大露天煤矿区全部建成后给生境质量带来巨大压力，生境质量又出现降低。统计各等级生境的面积比，较低等级生境占比最高(近50%)，表明平朔矿区生境质量整体处于较低的水平。1985—2015年,各等级生境面积占比有所变化，低等级生境逐年升高，由5.83%提升至22.93%，低等级生境主要由较低等级生境转化而来；较低等级生境分布着大量草地和耕地，安太堡、安家岭、东露天三大矿区的建设侵占了耕地和草地，导致较低等级生境质量下降31.70%，较低等级生境向低等和中等级转换。2005—2015年，中等和较高等级明显增加，分别上升14.67%、1.35%;中等级生境主要为草地，较高等级生境为林地，表明生态恢复取得效果，生境质量有所改善。高等级生境质量占比不大，变化也不是很明显。分析2025年的生境变化可知，2025年平朔矿区平均生境质量为0.42，较2015年下降了0.06，这与2005—2015年生境质量的下降趋势一致。在未来几年内，低等级和较低等级生境将持续增长，分别增长9.47%、20.99%，建设用地有扩大的可能性，然而中等级、较高等级、高等级生境质量均有所下降，分别下降26.99%、2.72%和0.75%，其中中等级生境质量下降最为明显，耕地有很大可能转为建设用地，林地和草地基本保持不变。平朔矿区具有潜在的生境退化风险，仍需要坚持复垦工作，在发展经济的同时，更要注重生态保护。
图6 1985—2025年平朔矿区生境质量变化情况Fig.6 Changes of habitat quality in Pingshuo mining area from 1985 to 2025
将1985年、1995年、2005年、2015年与2025年生境质量分布图进行差值分析，以深入探究生境质量的空间分异特征(图7)。1985—1995年生境质量下降区域集中于中西部，整体生境无明显变化。1995—2005年，生境质量提升区域遍布整个矿区，这得益于矿区土地复垦工作的不断完善与推广；生境质量下降区域则位于中部，主要是因为建设用地的迅速扩张挤压、占用了土地，直接破坏生境空间，生境质量极速下降。2005—2015年，中部的安太堡矿区生境恢复，安太堡的西排土场、南排土场和内排土场在2005年生态复垦的基础上进行了强化，生境质量升高。2015—2025年，平朔矿区整体生境质量相对均衡，仍然存在潜在的下降风险，说明生态复垦需要持续坚持。
图7 1985—2025年平朔矿区生境质量空间变化Fig.7 Spatial change of habitat quality in Pingshuo mining area from 1985 to 2025
4 平朔矿区生境质量归因分析
采用地理探测器对平朔矿区的生境质量影响因素进行分析。从因子探测的结果可知，决定力(q值)从大到小依次排序为NDVI(0.157 5)、煤炭生产能力(0.069 9)、GDP(0.069 7)、降水(0.065 5)、高程(0.064 5)、气温(0.045 6)、坡度(0.0307)。其中NDVI对生境质量的影响最大，贡献率达0.157 5，明显高于其他因子的解释力，且NDVI与所有的因子均存在显著差异，NDVI是引起生境质量变化的重要驱动力。
从交互结果得知，因子之间的交互作用对生境质量的影响均高于单个因子的影响程度，交互类型均为增强。交互作用解释力最强的为NDVI与其他影响因子的协同作用，解释力均高于24%，其中NDVI与GDP的叠加作用最大，可解释28.66%的生境质量变化空间分布；其次为气温∩高程(0.208 4)、煤炭生产能力∩气温(0.201 7)、煤炭生产能力∩高程(0.201 3)、煤炭生产能力∩降水(0.198 5)、高程∩GDP(0.189 8)，坡度∩GDP(0.112 0)最低。NDVI作为反映植被健康态势的重要参数，与生境质量的高低有着密切的关系，因此NDVI与其他影响因子交互作用较大。此外，煤炭生产能力与自然因子的交互作用加大了对生境质量变化的解释力，说明人为因素与自然因素的共同作用对生境质量变化的影响力更强。
5 结论
RFFLUS-InVEST-Geodetector耦合模型能够有效地分析平朔矿区1985—2025年的土地利用变化，精准揭示生境质量的时空演变规律，从而能对生境质量的影响因素进行定量归因分析，得出：(1)1985—2015年耕地、水域大面积减少，分别减少了11 130.21、327.87 hm2，耕地主要转化为林地和草地。建设用地、草地分别增加了7 169.59、3 656.16 hm2。2015—2025 年土地利用结构发生显著变化，建设用地和其他用地增幅明显，草地、耕地显著减少，林地和水域减少较小。(2)生境质量空间分布格局整体呈现由中心向外围逐渐增高的态势。1985—2015年生境质量呈现“降—增—降”趋势，生境质量值总体下降0.06。2025年生境质量值持续下降了0.062 3，中等级、较高等级和高等级生境均有所下降，平朔矿区具有潜在的生境退化风险。(3)在地理探测器中，各影响因子对生境质量的解释力从大到小依次为：NDVI、煤炭生产能力、GDP、降水、高程、气温、坡度；各因子的交互表现类型为双因子增强或非线性增强，其中NDVI和GDP的交互作用最大。
本研究运用InVEST模型中的生境质量模块分析了平朔矿区生境质量变化，因该模型的参数设置存在主观性，生境质量评估受到一定影响，未来可根据矿区的实际情况进行调查并改进；此外，该模型输入的胁迫因子和土地利用类型种类较少，未来应利用决策树分类将土地利用类型细化，考虑矿区开采和矿区恢复下的土地利用和生境质量变化，从而帮助决策者更加清晰地了解矿区生境演变规律，促进矿区更加健康的发展。
