大气中二氧化碳含量增加是国际社会广泛关注的环境问题, 如何减少碳排放、提升碳储存能力成为了研究热点和难点。陆地生态系统通过释放和吸收CO₂、N₂O等温室气体调节区域气候, 是全球碳循环过程的重要环节^[1]。陆地生态系统碳储量提高能够有效减少大气中CO₂的含量, 是缓解温室效应、减缓全球气候变化的重要方式之一。人类活动改变土地利用格局, 而土地利用变化进一步影响陆地生态系统结构、功能, 改变陆地生态系统固碳能力, 是陆地生态系统碳储量变化的主导因素^[2]。因此, 预测不同情景下土地利用变化, 定量评估土地利用变化所致的陆地生态系统碳储量, 为未来土地利用结构优化方向和陆地生态系统碳储量提升途径提供了指导作用, 为区域可持续发展和气候变化减缓做出贡献^[3]。

近年来, 众多学者基于土地利用格局变化对单一或综合性陆地生态系统碳储量进行了研究。研究结果表明, 全球范围内热带地区森林植被减少、森林生态系统功能退化造成了全球陆地生态系统碳储量的下降^[4-5]；陈耀亮等运用簿记(Bookkeeping)模型评估中亚地区森林生态系统碳储量, 指出森林砍伐造成了碳储量下降且土地利用/覆盖变化对碳储量有显著影响^[6]；刘亚男等研究了中国主要湿地碳储量并指出土地利用变化造成了湿地生态系统功能的退化, 严重影响湿地的“碳汇”功能^[7]；柯新利等人基于InVEST模型研究湖北省陆地生态系统碳储量, 认为其下降的主要原因是城市建设用地扩张和耕地保护^[8]。这些研究大多通过InVEST模型、Bookkeeping模型评估碳储量, 表明了土地利用变化对陆地生态系统碳储量有显著影响, 为土地利用结构优化提供了新的科学依据。一些学者运用模型模拟了未来土地利用格局与碳储量变化, 却没有考虑不同土地利用情景下的差异^[9-11]。流域作为一个复杂系统, 包括各种自然资源、生态系统等要素, 具有明显的层次结构和整体特征, 从流域尺度分析并解决生态环境问题是切实可行的有效途径^[12]。然而, 针对流域尺度, 结合时间和空间尺度预测不同土地利用情景与生态系统碳储量的研究并不多见, 预测未来土地利用格局, 模拟不同土地利用情景, 评估生态系统碳储量, 有助于区域生态环境发展, 促进土地利用优化配置。

本文选择汾河上游作为研究区, 该区域位于山西省中北部, 汾河发源地, 对山西省整体生态环境有重要影响, 近年来人类活动对土地利用格局产生巨大影响, 生态系统碳储存功能下降。因此, 本文以2007年、2011年、2017年汾河上游土地利用数据为基础构建了SD_CLUE-S复合模型, 预测2030年生态保护情景和自然增长情景下各土地利用类型面积, 利用InVEST模型对2007—2017年及不同情景下2030年生态系统碳储量进行估算并分析碳储量空间分布情况, 为汾河上游土地利用规划提供科学依据, 促进区域生态系统碳储存功能提升, 对山西省整体生态环境可持续发展具有重要意义。

1 研究区概况

汾河流域(38°7′57″—39°8′41″N, 111°42′43″—112°36′55″E)位于山西省忻州市中部, 东依原平市、忻府区, 西隔管涔山与五寨、岢岚县为邻、南连太原市的娄烦县、古交市, 北靠神池县、朔州市的朔城区(图 1), 总面积达到了399737.47hm²。空间上总体呈现“两山一川”大格局, 以芦芽山、云中山为两翼, 以汾河川为中心从南到北贯穿整个区域, 自然生态要素齐备, 生态功能完整。该地区位于北温带大陆性季风气候, 多年平均气温6.30℃, 最高和最低均出现在河川阶地区；降水多集中在7—9月份, 历年平均降水量420—704mm, 区域内平均年蒸发量为1812.60mm；境内山地多, 整体海拔较高, 地势呈中部低、东西两边高, 境内草地、林地、耕地所占面积较大, 建设用地面积逐渐扩大, 水体面积多年来持续下降。研究区是山西省重要的水源涵养区、华北平原重要生态屏障, 全省汾河水源地保护与治理重点区域, 生态功能系统碳储存功能脆弱是制约经济发展和社会进步的突出因素。

2 材料与方法 2.1 数据来源及处理

研究区2007、2011、2017年遥感影像来自地理数据空间云(http://www.gscloud.cn), 运用ENVI软件进行遥感解译且经过野外验证及Google Earth Pro验证精度达到了89.20%、91.45%, 90.87%, 根据《土地利用现状分类(GB/T21010—2007)》将研究区3期土地利用数据分为草地、耕地、建设用地、林地、未利用地、水体6种类型。

研究区年降雨量及年均气温来自《忻州市统计年鉴》；DEM数据来源于ASTER GDEM, 驱动因子中海拔因子和坡度因子数据来自DEM数据的ArcGIS处理；土壤质地因子来源于联合国粮农组织(FAO)和维也纳国际应用系统研究所(IIASA)所构建的世界土壤数据库以及第二次全国土地调查南京土壤所所提供的1:100万土壤数据；距离因子中省道、县道、乡道、河流、铁路、农村居民点以及县域中心的数据均来自全国地理信息资源目录服务系统, 运用欧式距离法得到研究区每个栅格到某空间驱动因子的距离；人口因子数据来源于地理国情监测云平台中全国1km×1km的人口栅格数据。

2.2 研究方法 2.2.1 系统动力学

系统动力学基于信息反馈系统认识并解决系统问题, 认为系统的行为模式与特性主要取决于系统内部的动态结构与反馈机制^[13], 一般用系统结构框图、因果关系图以及流图来构建特定模型进行长期预测, 以满足实际应用需要, 在研究数据不够充分的情况下尤其适用。

本文基于2007—2017年研究区域土地利用转移矩阵信息构建了土地利用变化的系统动力学模型(图 2), 各土地利用类型设置为水平变量, 土地转移面积设为流速变量^[14], 转移概率设为影响转换面积的常量, 并根据实际情况对转移概率进行修正和调整, 以便更精确地预测2030年土地利用格局。

2.2.2 CLUE-S模型

CLUE-S模型分为非空间需求和空间分配模块, 其中, 非空间需求模块确定土地利用变化驱动因子及其与土地利用变化分布概率之间的关系^[15], 本文通过二元Logistic回归模型及ROC曲线验证了驱动因子选取的精度, 同时根据土地利用的SD模型预测2030年研究区各土地利用类型面积。空间分配模块分为输入和输出两个阶段, 输入阶段包括土地利用转换规则, 限制性区域、空间特征以及非空间土地需求模块求得的土地需求结果^[16], 通过迭代将2030年土地利用类型的空间分配结果输出。

根据《忻州市土地利用总体规划(2006—2020年)调整方案》将各生态修复区域及基本农田区域选定为限制区域(图 3), 设定生态保护情景下草地、耕地、建设用地、林地、未利用地、水体的转移弹性值为0.5、0.6、0.7、0.3、0.7、0.6, 并设置土地利用变化转移矩阵(表 1)；自然增长情景下各土地利用类型ELSA值设为0.5、0.5、0.7、0.3、0.9、0.9, 同时将土地利用变化转移矩阵中水体转移为未利用地设为0, 其余皆为1, 并且不设置限制区域。

2.2.3 InVEST模型

InVEST模型根据不同土地利用类型对应的地上生物碳库、地下生物碳库、土壤碳库和凋落物碳库估算碳储量^[17], 并生成平均碳密度空间分布图。根据研究区具体实际情况以及数据获取难度, 本文选取对碳储量评估影响较大的地上生物碳库、地下生物碳库以及土壤碳库进行估算, 计算方法如下：

(1)

(2)

式中, C_kt表示土地利用类型k的碳密度总量(t/hm²)；C_kabove表示土地利用类型k的地上生物碳密度(t/hm²)；C_kbelow表示土地利用类型k的地下生物碳密度(t/hm²)；C_ksoil表示土地利用类型k的土壤碳密度(t/hm²)；C表示研究区碳储量(t); S_k表示土地利用类型k的面积(hm²)。

根据李克让、解宪丽等人的研究得到全国六种地类的碳密度数据^[18-19], 基于S. A. Alam、陈光水等人提出的公式^[20-21]对碳密度进行修正(表 2), 由于死亡有机物碳密度难以获得且所占比例不大, 因此将其填为0。

(3)

(4)

(5)

式中, C_SP表示年降雨量修正下得到的土壤碳密度(kg/m²)；MAP表示年均降雨量(mm)；C_BP和C_BT分别表示年降雨量和年均气温修正下得到的生物量碳密度(kg/m²)；MAT指的是年均气温(℃)。

(6)

(7)

(8)

(9)

式中, K_BP表示生物量碳密度年降雨量因子修正系数；K_BT表示生物量碳密度年均气温因子修正系数；K_B表示生物量碳密度修正系数；K_S表示土壤碳密度修正系数；C′和C″分别表示研究区及全国相应数据。

3 研究结果 3.1 2007—2017年汾河上游地区土地利用变化分析 3.1.1 土地利用变化幅度分析

2007年与2011年草地是研究区主要土地利用类型, 2011年与2007年相比草地面积下降且降幅为3%, 耕地、建设用地、林地以及未利用地面积增加, 其中林地增加幅度最大达到了2.29%, 其次是耕地增加幅度2.03%, 最小的是未利用地。此外, 2011年水体面积减少了7209.89 hm²(如图 4、表 3所示)。

2011—2017年草地仍然是土地利用的主要类型, 但是面积仍旧在下降且降幅达到了2.49%, 此外未利用地与水体面积存在不同程度下降, 其中未利用地面积减少了4080.85 hm², 主要原因是人们加大了对未利用地的开发程度；水体面积下降了0.24%, 近年来水体面积的持续下降对研究区域的生态环境产生了恶劣影响。耕地、建设用地、林地继续保持上升趋势, 其中耕地面积增加了9468.66 hm², 建设用地增加幅度达到了0.49%, 林地面积增加了3547.33 hm²(如图 4、表 3所示)。

2007—2017年, 耕地、建设用地、林地面积持续上涨, 增加幅度分别为4.40%、0.74%、2.18%, 草地、未利用地、水体面积下降, 其中草地面积降幅最大为5.49%, 其次是水体为2.04%, 减少面积达到了8144.77 hm², 未利用地面积则下降了3137.40 hm²(如图 5所示)。

3.1.2 土地利用类型转移矩阵分析

草地转出主要流向耕地, 其次是林地, 建设用地、水体与未利用地所占比重较小；耕地主要转移成了草地, 林地、建设用地、水体、未利用地次之；建设用地主要流向耕地, 转出比例达到了41.11%, 其次是草地；林地转化为其他土地利用类型的程度最小, 其中转为草地的比重最大；未利用地转移成林地比例最大, 其次是草地、耕地与建设用地；对水体的占用主要是草地与耕地, 两者的转移概率总和将近80%, 其次是建设用地、林地, 未利用地没有占用水体(如表 4、表 5所示)。

3.2 土地利用模拟预测结果及分析 3.2.1 基于SD模型的土地利用仿真及预测

根据系统动力学模型得到2007—2017年每一年各个土地利用类型面积的仿真结果(表 6), 将其与2017年解译结果进行对比, 确定SD模型构建的可行性和准确性。

根据仿真值与解译值对比可以看出, 所有土地利用类型仿真结果的相对误差都在1%以下(表 7), 说明SD模型仿真结果与实际情况基本一致, 可以进行下一步的预测, 即输入2017年各土地利用类型面积作为初始值, 设置步长为13年, 预测2030年研究区土地利用情况。

根据系统动力学模型预测结果可以看到(表 8), 研究区域在未来13年间各土地利用类型变化面积较大, 其中耕地变化面积最大, 减少了1117.00 hm²；林地增加量次之, 达到了1045.46 hm²；建设用地和未利用地面积分别减少了181.10 hm²和101.69 hm²；草地和水体分别增加了234.00 hm²和118.16 hm²。

3.2.2 土地利用变化驱动因子分析

本文选取了11个驱动因子进行分析, 其中距离因子包括距省道的距离、距县道的距离、距乡道的距离、距铁路的距离、距农村居民点的距离、距河流水系的距离和距县中心的距离；地形因子包括高程和坡度；土壤因子与社会经济因子分别为土壤类型和人口密度因子。

本文通过空间数据分析、二元Logistic回归分析和ROC曲线测算, 得到驱动因子的B值和Exp(B)(如表 9所示), 其中B值是各驱动因子的回归系数, Exp(B)表示发生比率, 即驱动因子每变化一个单位, 各土地利用类型发生比的变化；ROC曲线能够有效验证土地利用模型拟合度, 一般来看ROC值大于0.5即可表明驱动因子对土地利用类型空间分布具有较好的解释能力, 值越大表明拟合度越高, 根据计算得到ROC值均大于0.65, 拟合度较高, 符合要求。

研究结果表明, 研究区草地分布与县道、乡道、农村居民点、河流水系距离以及坡度呈正相关, 与海拔因子呈负相关, 其中坡度对其影响最大, 坡度每增加1°, 草地分布概率就增加2.9%；耕地与距农村居民点的距离、与县中心的距离、高程、坡度呈负相关, 与距河流水系的距离及人口密度因子无关, 与其他因子为正相关, 其中县域中心距离对耕地分布影响最大, 距离每增加1m, 耕地分布概率就减少62.5%, ；建设用地分布与县域中心距离相关性最强；土壤类型对林地分布影响最大, 将土壤类型设置为壤土-0, 砂土-1, 壤砂-2, 砂壤-3, 粉粘-4, 因此土壤类型对应编码每增加一个单位, 林地分布概率就减少14.96%；未利用地分布仅受到海拔与坡度的正相关影响, 其中海拔对其分布概率影响最大, 海拔每增加1m, 未利用地分布概率就增加36.89%；水体分布概率与县域中心距离相关性最强。

3.2.3 基于生态保护情景的土地利用变化模拟

将2007年数据作为起始年份数据, 仿真模拟2017年土地利用数据验证Kappa精度, 保证2030年土地利用数据预测准确性, 如图 6所示, 模拟正确的栅格有59973个, 总栅格数有63704个, 则P0为0.94, PC为0.30, PP为1, 得到Kappa系数为91.58%, 即2017年研究区域模拟精度较高, 完全满足进一步的研究需要, 输入自然情景以及生态保护下预测2030年土地利用数据的相关基础数据, 运行模型并输出结果(图 7)。

自然增长情景下, 林地面积大幅下降, 减少了近52908.26hm², 相较于生态保护情景减少了近54925.56hm²；自然增长情景下2030年草地面积达到了198100.00hm², 而生态保护情景下仅有145325.00hm², 主要原因是ELSA值、转移矩阵及空间分布概率影响；生态保护情景下退耕还林、还草等生态保护措施实施耕地面积下降, 与之相比自然增长情景下耕地面积增加了467.73hm²；建设用地在生态保护情景下发生增长而在自然增长情景下减少, 两种情景下面积相差较小；生态保护情景下大部分未利用地转化为林地、草地及水体, 且利用面积相较于自然增长情景多了1230.00hm²；研究区作为汾河的重要源头, 水体占比具有重要意义, 2007—2017年水体面积不断下降, 水资源环境持续恶化, 在自然增长情景下水体仍旧呈现下降趋势, 而生态保护情景下实现了142.16hm²的增长。

3.3 汾河上游地区碳储量评估

生态保护情景下, 2007—2030年研究区碳储量呈稳步上升趋势(图 8), 碳密度持续增长, 其中2007—2011年碳储量增幅最大, 达到了868423.72t, 2011—2017年碳储量增量下降为368719.31t, 2017—2030年研究区碳储量增加了341217.04t；2007年碳密度为144.45t/hm², 2030年碳密度达到了147.95t/hm², 且2007—2011年的碳密度变化最大；从土地利用类型来看, 草地碳储量在2007—2017年间持续下降, 2030年与2017年相比发生了51714.05t的增长, 2007—2017年耕地碳储量持续增长, 2017—2030年下降且降幅达到了160758.96t；建设用地碳储量在2011—2017年间增加量最大, 其中增幅最小为8550.58t；未利用地在2007—2011年间有小幅增长, 在2011—2030年持续下降, 且2011—2017年下降幅度较大；水体在2007—2017年间碳储量持续下降, 在2017—2030年发生了增长, 其中2007年下降幅度最大, 达到了11455.27t(图 9)。

自然增长情景下, 2017—2030年碳储量下降176917.86t, 碳密度变化0.41t/hm², 相较于生态保护情景碳储量减少341217.04t, 碳密度下降0.85t/hm², 从土地利用类型分析, 草地碳储量发生极大增长, 比生态保护情景多8176915.15t, 同时耕地与未利用地分别增长228058.52t和47338.88t；反之, 林地大幅度下降且减少了8771373.06 t, 建设用地下降21579.67t, 水体碳储量变化较小, 有576.87t下降量。

本文将碳储量变化量在≤-150t/hm²范围设定为下降, 处于-150t/hm²到150t/hm²设定为基本不变, ≥150t/hm²表示增加。由此可知, 2007—2011年碳储量下降区域主要以带状分布在静乐县汾河流域周围以及宁武县阳方口工矿镇区域, 其余较为零散地分布在研究区林地范围内, 增长主要发生在静乐县县城、宁武县芦芽山自然保护区、县内汾河及汾河支流附近。2011—2017年静乐县碳储量增长区域集聚分布在以汾河流域为中心的范围内, 而下降则发生在汾河流域内部；宁武县碳储量沿芦芽山自然保护区边界、恢河流域及县城周围村庄下降, 碳储量增长主要在县城中心范围、东马坊乡及荷叶坪区域(图 10)。

生态保护情景下2030年碳储量下降区域集中在研究区汾河流域, 增长区域零散分布在宁武县芦芽山自然保护区、恢河流域以及东马坊乡；自然增长情景下碳储量基本没有变化, 仅在宁武县恢河流域附近发生增长(图 11)。

4 结论和讨论 4.1 讨论

(1) 总体来看, 2007—2017年汾河上游社会经济发展, 生产活动频繁, 人们对耕地、建设用地的需求增加, 导致耕地、建设用地面积持续增长；宁武县工矿企业增加造成了水体污染, 同时人们生产生活中不合理利用水资源导致了汾河水量、水质不断下降, 水体面积逐年萎缩；汾河上游拥有大面积的原始森林, 随着政府实行林地管控及退耕还林政策, 林地面积发生增长。

(2) 系统动力学(SD)能够反馈系统内部要素相互作用, 借助模型对复杂问题实行运算, 能够有效预测不同情况下土地利用变化、土地利用需求等, 系统仿真提高了SD模型预测精度, 但是SD无法获得空间信息和空间分布预测结果^[22]。CLUE-S模型的空间模块能够模拟土地利用空间分布, 但是无法获得土地利用需求量, 在应用中人们往往会结合其他模型进行多领域的研究, 高志强等和将其与Dinamica EGO模型结合, 分析了中国未来土地利用格局和土地利用变化驱动力因子^[23]；周锐等结合Markov模型预测了辛庄镇在3种不同情景下未来土地利用格局的变化情况^[24]；刘静怡等结合灰色线性规划方法预测了嘉兴北部地区不同情景下的土地利用格局, 并进行土地利用结构的优化^[25]。然而这些方法主要运用数学方法进行土地利用需求预测, 忽略了人为决策因素, 或者需要大量数据支撑, 这在实际研究中可能难以获得, SD模型可以很好地解决这些问题。因此本文运用SD_CLUE-S复合模型消除了单一模型的不足, 避免了其他模型的劣势, 两者相辅相成, 提高了土地利用格局预测精度, 为进一步研究提供技术支撑, 但是不足之处在于SD模型构建中缺少了社会经济因素、政策因素考量, 这些影响因素的量化和空间表达是未来研究重点；此外空间尺度对土地利用格局影响较大, 本文分析一种空间尺度下研究区土地利用变化, 今后应当设置不同空间尺度对比分析。

(3) 目前陆地生态系统碳储量评估大多通过采样或植被、土壤清查数据计算获得, 具体方法可分为地表测量、模型模拟和遥感估算三类^[26], 其中地表测量大多是通过采样、实地勘测等, 结果精确性高, 但实施难度较大；遥感估算是当前研究的新方向和新思路, 通过遥感影像人们不仅可以获得地表信息, 还能对检测地下生物碳储量, 但遥感影像精度会影响研究结果；模型模拟方法中人们普遍采用InVEST模型, 数据获取方便, 结果准确性较高, 但是模型忽略了植被种类、植被生长情况等对碳密度的影响, 因此在研究成果基础上对研究区各土地利用类型碳密度进行修正, 这比以往直接引用全国碳密度值准确性高, 不足之处是缺乏实地测量数据对修正结果进行验证, 只能通过相近地区研究结果进行对比, 在今后研究中可以采用实验和预测模型或实地调研和预测模型相结合的方法, 提高碳储量估算精度。

4.2 结论

通过SD_CLUE-S模型模拟汾河上游生态保护情景和自然增长情景下未来土地利用格局, 并运用InVEST模型估算生态系统碳储量。结果表明, 2007年汾河上游碳储量和碳密度分别为57740767.95t和144.45t/hm², 2007—2017年汾河上游碳储量持续增长, 但增长幅度明显减缓。2030年汾河上游自然增长情景下碳储量和碳密度出现较为明显的下降, 与2017年相比分别减少176917.86t和0.44t/hm²。生态保护情景改变了这种趋势, 碳储量和碳密度显著增加, 增幅达到了2164299.18t和0.41t/hm², 增强了汾河上游固碳能力。未来对汾河上游土地利用结构优化, 应基于生态保护情景, 对汾河进行河道治理, 提高河流岸坡生态环境, 保证汾河水质、水量, 积极保护林地, 合理规划芦芽山保护区的旅游资源, 限制工矿用地过度扩张, 确保基本农田数量, 这样在提高碳储量的同时也可以保证汾河上游地区社会经济可持续发展。