工业园区是以工业为主的复杂经济系统，在相对较小的空间内聚集了大量企业，资源消耗量大，污染物排放集中[1]，园区的生态化发展对实现产业绿色转型[2]、节能减排[3]、减污增效[4]具有重要意义。
目前，针对工业园区绿色发展的研究主要集中在绩效评价方面：国家层面先后推出了国家生态工业示范园区创建、园区循环化改造、国家低碳工业园区试点、绿色园区建设等模式探索[5]，并出台了相应的指标体系对园区的发展水平进行评价；评价模型方面，数据包络（DEA）[6-9]、层次分析[10-12]、综合评价[13-14]等方法也被广泛应用于园区发展水平的定量评价中。绩效评价方法能很好地反映园区的发展现状[15]，但不能预测中长期发展战略下园区的发展趋势，很难从系统角度提出园区发展目标与优化路径。近些年，多目标优化法建立产业结构调整模型[16]、多目标优化与遗传算法结合建立产业模型[17]、动态投入产出法建立生态化决策模型[18]、系统动力学建立园区发展模型[19-21]等方法被越来越多地应用在园区的趋势预测和路径选择研究中。不同于多目标优化和投入产出模型依据数学逻辑推演得出最佳解，系统动力学方法依据系统实际观测信息，从整体出发得出多变量输出解，以其能容纳较多变量且准确便捷等优点越来越多地被应用于工业园区决策中。例如，Zhao等[19]采用系统动力学和灰色聚类方法对长春经济技术开发区的生态化发展进行系统评价和情景优化；瞿庆玲等[20]建立了人口–产业–环境系统动力学模型，对徐州经济技术开发区未来产业进行预测；郑斯瑞等[21]建立了人口–经济–环境系统动力学模型，提出某化工园区的发展对策。但是目前大部分研究缺乏对发展趋势的绩效评价，鲜少考虑以工业为主导的园区尺度下经济、能源与环境的相互制约关系，较难适应新形势下减污降碳协同的环保要求，且在指导工业园区生态化发展方面还需进一步优化和完善。
笔者综合运用园区趋势预测和绩效评价方法，从剖析园区经济–能源–环境（3E）系统结构出发，构建了工业园区3E系统动力学（SD）模型。以长三角某国家级生态工业示范园区为案例，设置不同情景进行趋势预测，预测结果对比HJ 274—2015《国家生态工业示范园区标准》[22]（以下简称《标准》）进行分析，提出在满足《标准》前提下的经济增速、污染减排、能耗强度、能源结构调整比例。
1. 园区3E系统动力学模型构建
1.1 系统动力学方法
系统动力学作为一种综合的仿真模型，不同于其他黑箱模型，它是通过揭示系统行为的内在结构，构造系统的基本结构和信息反馈机制，进而模拟与分析系统的动态行为[23]，是处理高阶、非线性、多重反馈的复杂时变系统的一种有效工具，目前已普遍运用于国家[24-25]、区域[26-27]以及园区[28]、行业[29]等不同尺度的发展模式研究中。
工业园区系统动力学研究主要是针对园区内部子系统之间的复杂关系而建立系统的结构模型，用因果回路描述系统内部结构框架，然后绘出表示系统要素之间逻辑关系的存量–流量图，进而用方程建立要素之间的数学关系，最终用计算机仿真软件（VENSIM）进行模拟分析。
1.2 工业园区3E系统因果关系分析
工业园区发展表现为经济的高速增长，资源的集约利用以及生态环境的美化，经济、资源和环境子系统之间存在着多重非线性和复杂的相互作用[30]，如何摆脱资源与环境对经济发展的限制作用，实现经济–资源–环境（ERE）系统的可持续协调发展[31]，是园区发展需要解决的关键问题。本研究主要考虑系统中因工业活动带来的经济增长、资源消耗和环境污染，暂不考虑由居民生活带来的资源消耗和环境污染。在资源系统中，暂不考虑工业水耗与废水中污染物排放关系，仅考虑由能源消耗带来的二氧化碳排放及大气污染物排放过程，故将资源系统简化为能源系统，主要针对工业园区3E系统开展研究。
在工业园区3E系统中，能源系统是经济发展的物质基础，环境是能量的载体，也是工业活动的空间支撑，经济发展是整个系统的核心要素[32]。经济、能源、环境3个子系统相互联系，相互制约，共同构成了工业园区3E系统，保障园区绿色可持续发展（图1）。在园区3E系统关键参数选取过程中，根据工业园区现状，参考《标准》中的部分指标。
经济子系统属于整个工业园区的动力子系统，是整个工业园区的核心部分。经济的增长一方面带来能源消耗和污染物排放，对环境子系统形成正向压力；另一方面增加环保投资和提高技术水平，对环境子系统形成反馈，进而降低能耗水平和污染水平[33]。经济子系统选取工业增加值作为园区经济发展的关键衡量参数。
能源子系统为工业活动提供能量基础，也造成了环境系统污染物的排放[34]。能源子系统选取能耗强度与能源结构作为关键参数，通过单位工业增加值能耗降低与能源结构调整减少对环境系统的压力。
环境子系统在整个工业园区中，属于约束系统。环境子系统对经济子系统起到支撑和制约的作用，一方面给经济子系统提供资源保障，另一方面也对经济的发展有限制作用。工业园区环境子系统主要考虑水、气和固体废物的排放，园区固体废物综合利用率一般较高，故暂不考虑固体废物指标。废水排放方面，选取废水COD及氨氮排放，通过单位工业增加值污染物排放强度与经济子系统建立联系。废气方面，二氧化硫的排放与能源系统中煤炭消耗有密切关系，选取二氧化硫排放建立能源与环境系统的联系；氮氧化物排放与能源消耗的关系暂不研究。减污降碳协同方面，考虑能源消耗引起的二氧化碳排放指标，通过能耗强度降低和能源结构的调整，达到减少大气污染物与温室气体排放的目的。在经济增长对污染物排放的负反馈方面，考虑将污染物减排因子以及固碳技术减排比例作为衡量指标。
1.3 工业园区3E系统动力学模型构建
通过对园区经济、能源、环境子系统之间关系的分析，结合3个子系统的相互联系以及各变量间的反馈关系，利用VENSIM软件建立园区3E系统动力学模型（图2）。该模型实际上是对3个子系统的耦合，其中用单位工业增加值综合能耗、单位工业增加值污染物排放、煤炭消耗二氧化硫排放系数、不同能源二氧化碳排放系数、减排因子等指标将3个子系统有机结合在一起。
工业园区经济–能源–环境SD模型包含3个存量、5个流量、8个内生变量、12个外生变量（表1），通过对外生变量的调控可实现对模型中所有参数的仿真模拟。
2. 模型模拟与校验
2.1 模拟对象
以长三角区域某园区为研究对象，该园区于2015年被评为国家级生态工业示范园区，园区内主导行业为电子电器、机械制造、轻工纺织等，园区发展属于同类园区的中上游水平。该园区内没有燃煤电厂，已实现集中供热，建有污水处理厂，基础设施完备。经过生态工业园区的创建，该园区2015—2019年实现工业增加值年均11%的增速，单位工业增加值综合能耗维持在0.3 t/万元（以标准煤当量计），单位工业增加值新鲜水耗维持在6 m3/万元，单位产值污染物产生量持续降低，符合《标准》指标要求。
模型的空间边界为该工业园区边界，模型基准年为2019年，主要历史数据为2015—2019年，模拟时间为2020—2030年，时间步长设置为1年，采用欧拉（Euler）算法。调控参数为不同能源种类占比、单位工业增加值年均增长率、单位工业增加值综合能耗、减排因子、固碳技术减排比例，模拟参数为工业增加值、综合能耗、废水COD排放、废水氨氮排放、二氧化硫排放、二氧化碳排放。
2.2 模型主要参数与特征方程建立
系统的特征方程主要用来描述变量之间的数量关系，可以借助趋势外推法得到方程式；对于不存在线性相关关系的变量，可以通过表函数的形式参与模型构建，模型中的关键参数来源主要为资料调研和文献参考。主要有以下几种：1）工业增加值、废水COD最终排放量、废水氨氮最终排放量等存量为流量的积分，主要与存量的初始值有关。2）工业增加值年均增长率、单位工业增加值COD排放量、单位工业增加值氨氮排放量等参数用均匀随机函数表示，根据园区该参数历史数值趋势外推得到。3）单位工业增加值综合能耗、不同能源种类占比等预测数据与预测年份不存在线性关系的，用表函数来表示各预测年份的变量取值。单位工业增加值综合能耗满足累计下降要求，不同能源占比根据园区目前的能源结构，结合碳达峰碳中和战略、园区的用能情况以及各种能源的发展潜力，并参考相关文献[35]，调整园区煤、电等能源占比。4）固碳技术发展与减排因子设置为减排量的百分比[36]。固碳技术结合国家碳达峰碳中和战略要求，拟从2025年开始投产使用，用阶梯函数表示。5）其他关键变量由各变量之间的相互关系确定。
2.3 模型不确定性分析与校验
该模型在构建时为了确定变量之间的数学函数关系，做了一些简化处理，给模型带来一定的不确定性。模型暂未考虑废气氮氧化物排放与能源消耗的关系，未考虑水资源消耗与废水污染物排放的复杂关系；在园区二氧化碳排放中，暂未考虑因森林生态系统变化引起的碳储量变化，未考虑电力二氧化碳排放系数的年际变化。污染物排放方面，将污染物排放因子简化为常数，通过减排因子来调控其减量。
模型校验时运用Vensim PLE软件，以2015年为基准年，对2015—2019年数据进行历史校验。仿真计算园区2015—2019年工业增加值、二氧化碳排放量、各污染物排放量，并与同期实际数据对比进行一致性检验，根据模拟结果，历年模拟结果与历史数据拟合误差均小于10%，在允许误差范围内[37]。因此，该仿真模型符合系统动力学模型要求，能够对园区3E系统进行模拟，且具有较高的可信度。
3. 结果与讨论
3.1 发展情景设定与模拟仿真
按照国家、长三角区域以及园区的发展规划及相关文件设定工业园区发展的不同情景（表2）。工业园区层面，以单位工业增加值综合能耗下降、单位工业增加值二氧化碳下降代表单位GDP综合能耗下降、单位GDP二氧化碳下降指标。
情景一：基准情景。经济子系统工业增加值年均增长率按照2015—2019年的平均增速10%增长，不对园区的能源结构做出调整，不进行减排和固碳技术的研发应用。
情景二：节能减碳情景。工业增加值年均增长率维持10%，同时单位工业增加值综合能耗5年累计下降率参考《中华人民共和国国民经济和社会发展第十四个五年规划和2035年远景目标纲要》中单位GDP能源消耗降低5年累计13.5%的指标要求，对接碳达峰碳中和战略，从2025年开始固碳技术的碳减排量为能源碳排放量的2%。
情景三：减污情景。在情景二的基础上，增加污染物的减量因子，由于技术进步或者政策约束，每年污染物的减排量约为总排放量的5%。
情景四：能源结构调整情景。在情景三的基础上，调整各能源种类占能源消费总量的比例。该工业园区内没有燃煤电厂，已全部淘汰燃煤锅炉，目前煤炭用量主要为工业原料，占比不大，未来由于园区产业结构由重转轻的调整，到2030年煤炭使用比例可逐步降到0；预计石油占比10年累计下降20%，天然气占比10年累计提升15%；产生二氧化碳的电力占比（不包含绿电）维持不变，可再生能源（计入绿电）占比10年累计提升10%。
情景五：经济适度增长情景。适度放缓经济增长水平，根据园区所在省市《国民经济和社会发展第十四个五年规划和2035年远景目标纲要》设置园区的工业增加值年均增长率为7%，适度放缓经济增长速度，节能减排、能源结构调整的力度与情景四相同。
情景六：加大节能减碳力度情景。根据国家碳达峰碳中和战略，在情景五的基础上，加大节能力度，使单位工业增加值综合能耗5年累计下降18%，并提高固碳技术的研发，从2025年开始，通过固碳技术减少二氧化碳量约为排放量的4%。
模型模拟结果主要关注该园区2019—2030年工业增加值、综合能耗、废水COD、废水氨氮、二氧化硫和二氧化碳的排放结果数据（图3）。工业增加值由经济增速决定，综合能耗总量受经济增速和能耗强度影响，废水COD及氨氮排放总量与经济增速和减排因子有关，二氧化硫排放量与经济增速、能源结构和能耗强度有关，二氧化碳排放量与经济增速、能耗强度、能源结构及固碳技术发展有关。由于经济的发展，带来了能源的增加以及污染物排放量的增长；随着能源结构调整，二氧化硫排放量逐步减少；随着能源结构调整和能耗强度下降，二氧化碳排放量缓慢增长，但依然没有达峰。
3.2 不同情景下园区发展水平评价
针对园区不同情景下的经济、能源、环境子系统的模拟预测结果，对比《标准》中的综合能耗弹性系数、主要污染物排放弹性系数、单位工业增加值二氧化碳排放年均削减率3个必选指标，进行园区的生态化发展水平评价。
按照《标准》的计算方法，以2019年为基准，不同情景下综合能耗弹性系数如表3所示。园区保证工业增加值10%年均增长率，达到单位工业增加值综合能耗5年累计下降14%（情景二、三、四），但不能达到指标要求（建设期年均增加率≥0，综合能耗弹性系数≤0.6）。情景五适度降低工业增加值年均增长率，在园区工业增加值年均增长率维持在7%，单位工业增加值综合能耗5年累计下降14%条件下，恰好可以达到综合能耗弹性系数≤0.6的要求。
不同情景下单位工业增加值二氧化碳年均削减率如表4所示，除基准情景外，其他情景通过能耗强度的调整，均能满足单位工业增加值二氧化碳年均削减率≥3%的指标要求。不同情景下单位工业增加值二氧化碳5年累计下降率结果如表5所示，情景二、三、四、五在2025年固碳技术减排2%的前提下依然不能满足《中华人民共和国国民经济和社会发展第十四个五年规划和2035年远景目标纲要》中单位GDP二氧化碳排放降低5年累计18%的约束性指标要求。通过加大能源强度下降比例，情景六可以初步满足国家对二氧化碳强度的减排要求。在目前的情景中，园区二氧化碳排放均未达峰，后期还需要加强固碳技术研发以及深度能源结构调整。
污染物排放方面，模型中未考虑废气氮氧化物的排放量，能源结构调整情景下，废气二氧化硫的排放量大幅下降，弹性系数为负值，对污染物的弹性系数影响较大，暂不考虑二氧化硫排放系数。以2019年为基准，不同情景下废水COD和氨氮弹性系数平均值如表6所示。在情景二、三、四中，工业增加值保持年均10%的增长率，废水污染物减量为排放量的5%时，并不能满足《标准》要求（建设期年均增加率≥0，主要污染物弹性系数≤0.3）。可以通过增加减排因子，或者降低经济增长速度来实现主要污染物弹性系数的指标值。废水COD和氨氮的去除技术已经相对比较成熟，且园区执行比较严格的水污染物控制政策，目前暂不增加因技术和政策变化引起的减排比例，情景五和情景六通过适度降低经济增长速度来达到指标要求。但在2026年之后，与2019年基准年相比，还需要加大减排技术应用和废水污染物减量政策调控才能满足《标准》要求。
减污降碳协同方面，主要考虑因为经济发展、能源结构调整、能耗强度降低带来的二氧化碳与二氧化硫协同减排力度。综合对比单位工业增加值二氧化碳（表4）与二氧化硫（表7）的年均削减率，情景二和三情景下，通过能耗强度的调控，单位工业增加值二氧化碳与单位工业增加值二氧化硫均有削减，且年均削减率相同。情景四、五、六，在能源结构调整的情景下，二氧化硫排放强度与二氧化碳排放强度均下降，二氧化硫的减排力度远大于二氧化碳的减排力度。这表明，能源结构的调整会带来温室气体与大气污染物的协同减排，在能源结构调整情景下，由能源使用带来的二氧化硫减排力度更大。
综上所述，模拟结果表明，园区在能源结构初步调整前提下（2019—2030年，煤炭占比逐步降到0，石油占比10年累计下降20%，天然气占比10年累计提升15%，电力占比维持不变，可再生能源占比10年累计提升10%），工业增加值每年保持7%的增长速度，单位工业增加值综合能耗5年累计下降率达到14%，废水污染物减排量为产生量的5%，可以满足《标准》中必选指标的要求，但不能满足《中华人民共和国国民经济和社会发展第十四个五年规划和2035年远景目标纲要》中对于单位GDP二氧化碳5年累计下降率达18%的约束性指标要求。适度增加节能力度，保证单位工业增加值综合能耗5年累计下降18%，可以满足国家规划要求。这为园区未来发展提供了一条满足《标准》的生态化路径选择，但该选择并不是唯一的，可以继续调控外生变量，求得多种解。
4. 结论
以某国家级生态工业园区为案例进行动态仿真，通过未来10年的仿真结果可知，在能源结构初步调整前提下，园区保持7%的增速，能耗强度5年累计降低18%，污染物减排5%，可以持续推进国家生态工业园区建设。但是在这种发展情景下，园区的二氧化碳排放量还未达峰。通过园区发展趋势预测，可为该园区绿色发展提供有效路径：1）园区目前经济的高速增长（10%）给能源和环境系统带来一定压力，园区可通过适度推广循环经济，延伸现有产业链，升级产业结构，将经济升级到高效率、低能耗、低污染的结构状态。2）能源结构调整与能耗强度降低直接影响大气污染物与二氧化碳的排放。园区层面可通过减少化石能源使用，发展屋顶光伏，提升能源密集产业的能源利用效率，推广节能减排技术等方式进一步调整能源结构，提高能源效率。3）科学技术是减污降碳协同发展的动力。未来园区要加大污染物减量技术研发、提高工业生产效率、加大固碳技术的研发与应用，尽早实现二氧化碳的排放达峰。
该模型在构建和模拟过程中，由于部分指标间关系难以构建，部分指标数据缺失，在参数选择时受到一定限制。继续研究能源与大气污染物排放关系，构建能源系统与大气氮氧化物排放方程，考虑园区森林生态系统的固碳作用，增加二氧化碳排放达峰情景判定，建立园区经济系统产业结构变化与能源和环境系统的关系，都是有待深入研究的方向。
