基于SD模型的晋北退耕还林工程与社会经济互馈关系研究

郑子潇,胡保安,韩海荣,康峰峰,程小琴,*

1 北京林业大学生态与自然保护学院,北京 100083

2 内蒙古七老图山森林生态系统国家定位观测研究站,赤峰 024400

退耕还林生态工程作为一项重要工程技术,具有明显的社会经济特征,具备保护生态环境作用的同时,也影响着社会经济的发展[1]。大规模的生态修复是复杂的系统工程,退耕还林工程的有效性和社会经济因素密切相关[2],社会经济又不断影响工程的实施质量[3]。退耕还林工程实施的最直接表现是植被的空间变化,研究表明,耕地面积、造林面积与生态工程关系密切,主要表示生态工程的“量”[4—6];生态工程带来生态环境的变化,“固碳”、“释氧”、“归一化植被指数(NDVI)”等成为生态工程“质”的代表[3,7];人口、地区经济、农村经济是影响许多发展中国家大规模植被恢复工程的主要社会经济因素[8—9]。多数研究表明退耕还林工程对社会经济产生积极影响,改善了生态环境,增加了地区生产总值(GDP)、林业产值等[4—6],而一些认为随着工程发展,退耕还林对粮食总产量、农村收入产生负面影响[4,10]。如何科学、全面地理解和把握退耕还林工程与社会经济间的关系是未来调整、改进退耕还林工程的关键。部分学者就不同退耕还林模式对社会经济的影响进行研究:退耕还林强度与粮食总产量成反比[11];退耕还林工程的持续推进对农村收入有显著正效应[12]。因此,迫切需求在不同退耕还林模式下探求退耕还林工程与社会经济的内在联系和互馈作用机制,以期为科学推进生态文明建设、促进区域生态-社会-经济协同发展提供理论支撑。目前,针对生态与社会经济关系的研究有以下几种方法:系统动力学(System Dynamics,SD)模型[13—14]、结构方程模型[6]、耦合协调度模型[15]等。学者们运用SD模型对水资源利用[16]、土地利用[17]、生态系统服务价值[18]与社会经济关系的研究较多:袁绪英等[16]针对水资源短缺、水污染过重、水生态平衡破坏的问题,构建滠水河流域SD模型,进行模拟预测,来寻求促使经济系统与环境系统协调发展的措施与途径,从而为滠水河流域可持续发展提供政策建议;易阿岚等[18]运用SD模型对上海市湿地生态系统服务价值的变化进行模拟与预测,结果可为上海市制定湿地保护或生态规划政策提供可行性分析。综上,SD模型在分析系统结构、模拟复杂系统研究中具备优势,可满足定量分析退耕还林工程与社会经济互馈关系和未来发展趋势的研究需要。生态稳态转换可理解为生态系统发生巨大、突然、持续改变时,往往对社会、经济系统产生相当大的影响[19],山西省晋北地区退耕还林工程实施的改变,会对当地GDP、第一产业产值、粮食总产量等社会经济变量产生一定影响,因此同样可视为一种稳态转换事件[20]。

因此,通过构建SD模型对山西省晋北地区退耕还林工程与社会经济系统进行仿真和预测,揭示其互馈作用机制,以造林与耕地之间的用地矛盾作为参数调整切入点,考察不同退耕还林模式对社会、经济、生态要素的影响,有助于为晋北地区退耕还林工程发展提供可操作性的参考建议,促进该地区退耕还林工程与社会经济和谐可持续发展。为此,本研究旨在(1)通过SD模型定量表征晋北地区2002—2035年退耕还林工程和社会经济互馈关系;(2)阐明不同退耕还林模式,对晋北地区社会、经济、生态要素的影响,提出退耕还林工程实施的最优情景。

1 晋北退耕还林工程与社会经济SD模型构建

在Vensim-PLE软件中设定2002年为初始年份,2035年为终止年份,以2002—2019年统计年鉴数据进行拟合,时间步长为1年。本研究基于对退耕还林工程、社会经济特征的理解并参考前人研究[2—10,14],主要关系如图1,其中退耕还林工程的直接表现有“质”、“量”两部分,植被固碳量、NDVI可作为“质”的代表,造林面积、耕地面积等可作为“量”的代表,在本研究中造林面积、耕地面积作为系统中的主要驱动器,为经济社会提供资源,经济具备资金投入的功能,社会消耗资源[21]。选取退耕还林工程变量包含造林面积、耕地面积、NDVI、植被固碳量;社会变量包含总人口、农村人口、农林牧渔人口、粮食总产量、粮食单产;经济变量包含GDP、人均GDP、第一产业产值、林业产值、固定资产投资、植树造林投资、农村总收入、农村人均纯收入(表1)。造林面积等变量不符合线性关系的运用Vensim-PLE中表函数进行输入,一般方程借助R 4.0.3软件进行线性回归,如粮食总产量、农村人口等,因篇幅有限,以2002—2019年粮食总产量与总人口、耕地面积的回归关系进行详细说明(表2)。系统动力学模型流图见图2。

表1 SD模型主要变量及方程Table 1 Main variables and equations of system dynamics model

表2 耕地面积、总人口与粮食总产量回归结果Table 2 Regression analysis results of cultivated land area,total population and grain yield

图1 山西晋北退耕还林工程、社会经济系统关系图 Fig.1 Relations between the Grain for Green Project,socio-economic system in northern Shanxi

图2 晋北退耕还林工程与社会经济系统动力学模型流图Fig.2 Dynamics model of Grain for Green Project and socio-economic system in northern Shanxi

2 数据与方法

2.1 研究区概况

本研究所选取的晋北地区位于山西省北部(38°42′—40°08′N,112°73′—113°30′E),共包含三个地级市,从北向南依次为大同市、朔州市、忻州市,总面积约为49836km2,地处恒山、太行山、吕梁山地,地势呈现东西两边高、中间低,土壤类型以棕壤和淋溶褐土为主,属大陆季风气候,平均海拔约1347m。该地区属温带落叶阔叶林带,植物种类丰富,主要乔木树种有华北落叶松(Larixprincipis-rupprechtii)、油松(Pinustabuliformis)等;灌木树种主要有忍冬(Lonicerajaponica)、绣线菊(Spiraeasalicifolia)等。晋北地区具备一定可研究性:2019年,晋北地区GDP达3359.56×108元,人均GDP达3.99×104元/人,农村总收入达633.52×108元,相比山西省其他地区经济发展稍显缓慢;该地区生态环境较脆弱,自2002年全面实施退耕还林工程,2000—2015年间晋北地区的耕地转林地、草地占比在整个山西省中占据较高的比例,其中广灵县(大同市)、右玉县(朔州市)、保德县(忻州市)为重点退耕县[22]。

2.2 数据来源及处理

本文使用的研究数据主要包括山西省晋北地区2002—2019年退耕还林工程、社会、经济数据。其中社会经济数据来源于各年的《山西统计年鉴》,造林面积、耕地面积数据来源于《山西统计年鉴》、《中国林业统计年鉴》、山西省统计局(http://tjj.shanxi.gov.cn),植被固碳量数据来源于中国碳核算数据库(https://www.ceads.net.cn),NDVI基础数据来源于美国国家航空航天局(NASA)MOD13Q1全球地理信息数据产品并运用ArcGIS 10.2对NDVI进行分区统计处理。其中晋北地区人均GDP、农村人均纯收入等数据采用加权平均方法计算;对于个别缺失的统计年鉴,根据相邻年份的数据,利用插值法或灰色预测方法[23]补全。运用Microsoft Excel 2019进行数据的整理。

2.3 研究方法

系统动力学(System Dynamics,SD)于1956年由Forrester教授提出,是发展最早、最常用的动态模拟技术[24]。SD模型可以充分反映系统之间的相关关系和相互作用,达到定性与定量结合的效果,用于处理和研究发生巨大、突然和持续变化的系统,模拟复杂的系统过程,是识别稳态转换的模型方法[20],其重要特点是关系可视化和模型仿真[25]。本研究借助Vensim-PLE 8.2.1软件构建山西省晋北地区退耕还林工程的变化与发展的SD模型,并结合社会经济因素构建模型流图。运用Origin 2022进行作图。

3 SD模型检验

SD模型的检验实质上是一个证伪的过程,一般情况有以下几种常见检验方法:运行检验、历史检验、灵敏度检验、极端情况测试、积分误差测试、心智模型测试等[26]。本研究选取前三种方式进行模型检验。

对模型进行运行检验,在软件中显示“Model is OK.”,说明该模型通过了运行检验。此外,本研究从社会、经济、退耕还林工程三方面中选取GDP、第一产业产值、林业产值、固定资产投资、农村总收入、总人口、粮食总产量、耕地面积、造林面积九个变量进行模型历史检验,以2002—2019年统计年鉴数据为基础,对以上九个变量的模拟值和实际历史值进行误差计算。结果表明(图3),2002—2019年间关键变量误差均小于10%,符合模拟的一致性检验要求。对模型进行灵敏度检验,灵敏度检验是通过增加或减少某一变量参数的10%,探究模型运行后其他变量的变化情况[27]。本文选取十个主要变量进行灵敏度检验,除粮食总产量对总人口变化灵敏度约为10.44%左右外,其他变量灵敏度均小于10%(表3)。

表3 2002—2019年关键变量系统灵敏度分析结果Table 3 System sensitivity analysis of key variables in the system dynamics model from 2002 to 2019

图3 2002—2019年关键变量检验结果Fig.3 The system test of key variables in the system dynamics model from 2002 to 2019

综合以上对模型进行运行检验、历史检验和灵敏度检验的结果,认为该模型具备一定可信度,可进一步进行仿真操作。

4 情景模拟与结果分析

4.1 情景模拟参数设定

为了解不同退耕还林模式下晋北社会经济的情况,根据控制变量原则,结合相关文献[28—30]对退耕还林参数进行情景设定。参考2002—2019年以来的耕地面积变化率、造林面积的变化趋势及数据进行后续的计算[28—29]。山西省人民政府表示晋北地区在保护耕地的同时,也要持续进行生态造林的实施,晋北地区2020年耕地面积应不低于9942.53km2,预计在2035年山西省人均GDP至少达到2万美元;Pearson相关分析结果显示,造林面积与耕地面积变化量呈现显著负相关关系(相关系数=-0.825,P<0.01)。针对造林与耕地之间的矛盾,本研究设计四种方案(表4),模型仿真时长涵盖“十四五”和“2035远景目标”规划期。

表4 山西晋北地区退耕还林工程-社会经济模型仿真情景设计Table 4 Scenarios design of Grain for Green Project and socio-economic system in northern Shanxi

4.2 情景模拟结果与分析

经过耕地面积变化率的调整,持续退耕还林型(情景1)耕地面积持续降低,低速退耕还林型(情景2)、保护耕地型(情景3)、间断退耕还林型(情景4)耕地面积均呈现上升的状态,耕地面积上升速度情景3>情景4>情景2>情景1。造林面积在情景1、2、3分别保持在2002—2019年造林面积平均值、二分之一平均值、四分之一平均值左右,情景4 2025、2030、2035年呈现较高状态,其他年份则保持最低造林面积。情景模拟结果显示植被固碳量与NDVI值变化趋势类似,情景2、3、4植被固碳量与NDVI值均呈现不同速率的上升态势,只有情景1植被固碳量与NDVI值呈现缓慢下降趋势。在2024年前,植被固碳量大小表现为情景1>情景2>情景3,情景4受到造林面积参数变化的影响在51.62—57.87万t间波动上升,整体来看,情景1以-0.89%的速率缓慢下降,情景2未来增长趋势较平缓,情景3以1.18%的年增长速率快速上升。在2020年前,NDVI值呈现波动上升态势,2028年前情景1NDVI值较高,情景3NDVI值最低,情景4NDVI值围绕0.610左右波动,结果表明情景3NDVI值年均增长率最高,为0.85%,呈现快速上升态势,并且在2032年之后赶超情景1、2。

情景模拟结果显示,经济发展均呈现不同程度的上升态势(图4)。GDP 、第一产业产值、人均GDP、农村总收入在2020—2035年间呈现逐渐上升的趋势,均在2019年后呈现情景2>情景3>情景4>情景1,2035年GDP分别达18412.90×108元、17367.80×108元、16742.50×108元、15480.50×108元(表5),年均增长率分别为9.34%、8.93%、8.67%、8.10%;林业产值呈现情景1>情景2>情景4>情景3,其中情景4受到造林面积参数变化的影响,呈现波动上升,各情景从大到小年均增长率分别为7.57%、6.91%、6.81%、6.52%。

图4 2002—2035年晋北退耕还林工程、经济社会的相关预测Fig.4 Predictions of Grain for Green Project and socio-economic system in northern Shanxi from 2002—2035

时间Time 情景SceneGDMJ/km2ZLMJ/km2NDVIZBGT/104tGDP/108元LYCZ/108元DYCYCZ/108元ZRK/104人LSZC/104 t2002年10787.502003.400.59149.26502.193.9518.26750.79183.102020年110547.701022.300.61457.384812.9923.47265.21860.78467.90210547.70515.900.60253.324812.9923.47265.21860.78467.90310547.70262.700.59651.304812.9923.47265.21860.78467.90410547.701022.300.61457.384812.9923.47265.21860.78467.902025年110358.101072.470.61357.107126.3428.68358.99895.31558.09210554.20564.870.60353.827557.8426.15394.45895.66612.95311498.10310.940.59550.977416.6924.88393.08895.56617.45410738.101105.080.61958.887320.7525.68371.09895.38577.052030年110194.401131.540.61356.9210519.8034.16474.63932.01653.82210580.70622.640.60554.3911821.0029.09530.62932.78717.55312541.70368.410.60855.1811376.4026.55515.44931.96707.80410950.001201.100.62360.4911092.4028.25492.84931.98672.802035年110036.501205.090.61356.8715480.5039.96619.86971.06755.51210607.50694.870.60755.0718412.9032.34710.43972.20828.62313665.10440.880.62159.8117367.8028.53674.92969.87801.79411163.601317.820.62962.2816742.5031.27649.73970.62774.00

社会各指标同样呈现上升态势。总人口差距不明显,到2035年情景2可达到最高972.20万人。在2027年前,情景3粮食总产量最高,其次为情景2、情景4、情景1,2027年之后情景2粮食总产量超过情景3,且2020—2035年间情景2粮食总产量平均增长率始终高于情景3;在2023—2035年间情景2粮食单产均保持最高,上升了近54.18%,情景3粮食单产上升速率小于其他情景。

从图4中可看出,由于造林面积、耕地面积、植树造林投资的变化,情景2粮食总产量、粮食单产社会变量的增速较高,说明情景2有利于粮食总产量的提升,而情景1实施反而会限制粮食总产量,情景3由于保护耕地使得耕地面积快速增加,粮食总产量、粮食单产在模拟前期往往呈现较高水平,但是整体粮食总产量、粮食单产年增率却小于情景2。情景1经济变量呈现最低,情景2更利于经济的发展,更高的经济水平会带来更高的植树造林投资,从而形成良性循环。

5 讨论与结论

5.1 讨论

森林和耕地的土地利用竞争已经进行了数千年[31],近年来多项政策中均提到未来发展既要满足对耕地的保护,又要不断提升生态环境的保护,以上就涉及到耕地面积与造林面积之间的矛盾,这是本研究选择变化参数的依据,抓住耕地与造林之间的矛盾,通过改变退耕还林模式探求对社会经济的影响。本研究方程借助多元线性回归方程探究退耕还林工程与社会经济之间的关系,绝大多数方程关系拟合较好,而呈现非线性关系的运用表函数进行构建,使拟合结果更接近实际情况[16,18]。

与相关学者研究类似,晋北地区退耕还林工程直接对经济产生一定积极影响[15],在整体上对经济结构优化起到提升作用[2,8—9],然而并非越大规模的退耕还林带来越好经济发展,持续快速的退耕还林可能会给经济发展带来制约作用[32],情景1下各经济指标均低于其他退耕还林模式,情景2各经济指标最高。退耕还林工程直接增加了林业产值,造林面积是退耕还林工程“量”的重要指标,林业产值主要受到造林面积的影响,情景1林业产值最高,此时造林面积的提升大大增加了林业产值,这种正相关关系会随着工程实施年数的增加而越显著[33]。耕地作为土地利用方式之一,其面积是评价退耕还林工程“量”的指标,在粮食安全中发挥着巨大的作用,从而粮食总产量将是判断社会发展的重要指标[31],粮食单产作为粮食总产量的基本要素,是提升粮食总产量的重要方式[4]。本研究中情景1由于其耕地的快速下降粮食总产量一直呈现最低,2027年前情景3最高,而在之后情景2最高,情景1的实施在初期不会造成显著的粮食安全问题,但是长期可能造成负面影响[34—35],2027年之后情景2粮食单产和粮食总产量均保持最高,表明低速退耕还林的实施在长期内未对粮食安全产生危害[36],随着时间增长农业的进步可以保证在退耕还林过程中的粮食安全[4],在保证粮食安全的情况下,应及时调整退耕还林模式。

2024年前,情景1植被固碳量和NDVI值均高于其他情景,是因为植被固碳量和NDVI值受到造林面积影响多[4,37],但是情景1植被固碳量、NDVI呈现缓慢减少的趋势,因为相比较植被自然恢复,盲目大规模增加植树造林可能引起更严重的沙漠化进程[38],另外,植树坑会加重土壤侵蚀,进而导致环境恶化[39]。申强等[40]对比高退耕态势的吴起县和低退耕态势的米脂县,发现较低退耕态势的米脂县产生更多间接农地资源生态服务价值,本研究显示到2032年情景3植被固碳量和NDVI值均超过情景1,印证了低退耕还林态势的情景3更利于生态长期积极发展[40]。综上,在社会、经济、生态方面按照原计划持续退耕还林往往弊大于利[41],应该采取低速退耕还林型。

退耕还林对社会经济的影响有时空差异。实施退耕还林工程应分区域讨论,充分考虑地域差异,刘璨等[42]发现退耕还林造林面积对我国东部和中部地区的粮食产量产生的负面影响要高于西部地区;王子婷等[43]选取陇中近郊、远郊村落探讨退耕还林工程对各自农村经济的影响研究,发现退耕还林对远郊地区的社会经济影响较小;刘晓琳等[44]探究了黄土高原退耕还林对粮食可持续性影响的空间演变,发现包含了晋北的土石地区粮食可持续性指数在黄土高原中最高,这印证了本研究晋北地区退耕还林利于粮食总产量的发展。退耕还林可以在短时间内快速影响社会经济的发展,而对于植被固碳量等自然变量的影响会有滞后效应[45]:在模拟时间中各社会经济变量随着时间的增长呈现较快速的增长,而植被固碳量在相同时间段内的变化在47×104t—63×104t之间,其变化幅度远小于社会经济变量变化幅度。同时,社会经济同样影响了退耕还林工程,本研究通过植树造林投资的增长,短期内在政策上对造林、耕地的变化产生直接影响,增加了植被固碳能力[46]。

2035年,四种情景下GDP年均增长率分别达8.10%、9.34%、8.93%、8.67%,且四种情景均满足2035年人均GDP达到2万美元/人的目标(按现行汇率计算),退耕还林工程实施的最终目的是在保证良好生态环境的条件下丰富和优化人民生活,因此本研究选取2035年社会经济指标最高的情景2作为最优发展情景。生态工程与社会经济间的互馈关系依旧处于探索阶段,退耕还林工程影响了GDP、粮食总产量等社会经济变量,但是还存在时空差异,探讨气候、生物、土壤等对生态工程的影响和适应策略将是进一步研究的课题[47];此外对于2035年的定量规划数据尚不明确,此后可依据充足政策数据进一步探讨山西省不同区域退耕还林工程的最优情景。

5.2 结论

本文通过构建山西省晋北地区退耕还林工程和社会经济之间的系统动力学模型,以该地耕地和造林的矛盾作为参数调整的切入点,通过调整退耕还林模式,设计四种发展情景,阐明了不同退耕还林情景对社会经济的影响。基于2002—2019年历史数据,构建的退耕还林工程和社会经济之间的SD模型通过了运行检验、历史检验、灵敏度检验,实现了较好的模型性能,具备一定可操作性,可通过调整参数进行模拟预测。对比四种情景,按原计划持续实施退耕还林会制约未来社会经济发展,情景2利于退耕还林工程与社会经济的和谐发展,情景3在后期促进NDVI、植被固碳量的增长:2020—2035年间GDP、第一产业产值、农村总收入、人均GDP均为情景2>情景3>情景4>情景1。林业产值表现为情景1>情景2>情景4>情景3。研究初期情景3粮食总产量较高,而后期情景2最高,年均增长率达4.66%。情景3 NDVI、植被固碳量的增长率最高。到2035年,山西省晋北地区以情景2作为最优发展情景,使得耕地面积年变化率保持在0.35%状态,造林面积保持在694.87 km2左右,相比其他情景此时各社会经济变量较高,能较好地统筹耕地与造林用地的矛盾,在保护生态建设的同时为社会经济可持续发展提供保证。