生态系统服务是指人类从生态系统中直接或间接获取的益处, 其对人类福祉的贡献是全球人均国内生产总值(GDP)的两倍以上^[1]。湿地是地球上初级生产力最高、最有价值的生态系统之一, 是重要的自然资源^[2]。湿地在生产食物、降解污染物、改善水质、减弱风暴潮、减少洪水、维持物种多样性、气候调节以及文化价值等方面发挥着重要作用^[3-4]。随着人们对湿地认识的不断加深, 湿地提供的生态系统服务越来越受到人们的重视^[5]。尽管中国政府出台了相关保护政策, 湿地退化问题仍然严峻, 因此需要将生态系统服务纳入到湿地保护政策中^[6]。对生态系统服务进行估值可以量化这些服务价值的变化对人类福祉所产生的收益或损失, 还可用于指导政策制定和土地管理, 以控制生态系统服务变化的人类驱动因素^[7]。此外, 研究湿地生态系统服务价值有助于探索资源和经济的可持续发展^[8]。

在城市化进程中, 湿地面积不断缩减, 湿地生态系统服务严重退化^[9]。湿地面积的变化和人口的变化会引起湿地生态系统服务的变化^[10]。湿地管理者需要对生态系统服务进行可靠、合理地评估, 即对生态系统服务进行全面的生物、物理以及经济评估才能辅助决策者做出科学决策^[11]。在过去20年中, 研究者们运用各种方法进行生态系统服务价值(ecosystem service value, ESV)评估, 例如, 实际市场法、替代市场法、虚拟市场法等^[12]。Costanza等提出的ESV估算方法以及谢高地提出的单位面积价值当量是目前ESV评估的基础并被广泛引用^{[1, 13]}, 但参数的时效性存在着不确定性。虽然现有的ESV评估方法各有特点, 但绝大多数研究都以千年生态系统评估的分类方法为依据将每种ESV进行单独评估, 很少能将生态系统服务与自然及人类影响因素结合起来进行综合评估。此外, 在评估ESV时应考虑单位面积价值当量的时效性。将生态系统服务纳入到政策议程中已成为一种趋势, 可以协调社会经济发展和环境保护之间的冲突^[14]。

上海市位于东海沿岸线和长江口沿岸, 该地区的沿海和河口湿地在动态变化中, 受到自然和人为的影响^[15]。2003—2013年, 上海市因经济发展和城镇建成区扩张损失了约5.05×10⁴hm²的湿地^[16]。上海市湿地面积的变化, 一方面是由于长江上游输送的沉积物导致河口湿地扩张; 另一方面是由于人口的增长、城市化进程的加快以及海平面上升所带来的海水侵蚀导致湿地面积不断缩减^[17]。湿地面积减少导致上海市ESV在1989—2010年间损失了约4.45亿美元^[7]。现有研究不能很好地解释湿地面积变化对生态系统服务的影响程度, 也未能量化不同影响因素对生态系统服务的影响程度, 有的甚至忽略了自然因素对生态系统服务的影响^[18-19]。河口湿地生态系统影响因素以及空间异质性大, 要对湿地生态系统进行评估, 需要综合湿地生态系统各个方面的评估指标^[20]。现阶段自然因素和人为因素对生态系统服务的影响研究大都处于定性水平或者仅研究两个生态系统服务之间的协同或权衡关系^{[10, 21-22]}, 某一特定时期的ESV价值估算众多, 对未来生态系统服务动态变化的预测很少^{[19, 23]}。综上, 虽然对ESV的评估研究众多, 但仍缺乏对城市内湿地生态系统的动态评估以及对未来生态系统服务变化的预测。本文构建的系统动力学(system dynamics, SD)模型, 借助Vensim PLE建模平台对上海市ESV进行仿真模拟及预测, 并考察各因素变化对生态系统服务的影响, 将有助于决策者在将来制定有效的政策。为此, 本文旨在：(1)通过构建SD模型对上海市1993—2025年湿地生态系统服务变化过程进行模拟与预测; (2)分析1993—2025年上海市湿地ESV构成的变化; (3)识别影响上海市ESV变化的主要驱动因素。

1 SD模型构建

SD采取定性分析和定量分析相结合的方式来解决实际问题, 模型可以实现动态模拟^[24]。SD模型本质上是一个微分方程系统, 通过数值分析来模拟复杂系统的行为^[25], 且方程中的参数可根据实际情况调整, 从而使得模拟结果更加贴近实际情况并可进行多方案的比较分析^[24]。SD模型利用系统动力学的原理和方法对上海市快速城市化导致的湿地面积变化进行模拟分析, 并结合人口因素以及经济因素构建模型流程图(图 1)。

将表 1中的方程所涉及的参数输入Vensim PLE平台以搭建供给服务价值、调节服务价值以及文化服务价值计算子系统, 再根据式(1)建立ESV模型, 之后将湿地面积、湿地保护系数等参数纳入系统, 考虑到上海市湿地面积的变化并非线性, 湿地面积增长率采用表函数的形式输入模型。模型主要包含人口因素子系统、供给服务子系统、文化服务子系统、调节服务子系统中以及经济因素子系统。结合上海市人均GDP和上海市人均ESV可输出ESV/GDP的值, 湿地保护系数采用“IF THEN ELSE”函数, 即湿地保护系数的值由ESV/GDP值确定, 根据湿地保护系数就可以判断是否要采取适当的湿地保护政策来调节湿地面积增长率, 在本文的研究中, 湿地保护系数值取1, 即不采取任何保护政策的情况。通过ESV/GDP的值, 可以判断ESV在经济发展中的地位, 进而给环境规划与管理提供科学依据。此外, 上海市湿地生态系统的发展状况可根据ESV、人均ESV来判断。

2 材料与方法 2.1 研究区概况

上海市位于中国东部, 地处长江入海口(30°40′—31°53′N, 120°52′—122°12′E)。2018年末, 上海市行政区划面积为6340.5km², 除西南部有少数丘陵山脉外, 绝大多数土地为冲积平原, 平均海拔约2.19m, 全市常住人口2423.78万人。上海市被认为是中国现代经济的枢纽, 在全国经济建设和社会发展中具有十分重要的地位。该市拥有丰富的河口、河流、运河、溪流和湖泊等湿地资源, 包括两个国家级湿地站点, 四个自然保护区, 十三个市级湿地。上海市湿地主要分布于郊区, 近30年来, 湿地面积变化显著(图 2)。

2.2 数据来源

本研究中的土地利用/覆盖数据主要采用监督分类方法结合目视解译分别从Landsat TM/ETM遥感影像(1990年和2000年)、Landsat-TM遥感影像(1995年、2005年以及2010年)、Landsat 8遥感影像(2015年和2018年)获取。遥感影像无云, 空间分辨率为30m。在ENVI 5.3和ArcGIS 10.2软件中进行辐射定标和波段融合。根据《全国湿地资源调查与监测技术规程(试行)》将上海市土地利用/覆盖数据分为“湿地”和“其他用地”两种类型。“湿地”包含近海与海岸湿地、河流湿地、湖泊湿地、沼泽湿地以及人工湿地; “其他用地”则为“湿地”以外的土地(图 2)。上海市地处长江入海口, 长江流域上海段水下地形复杂, 故本文的研究区范围不涉及长江流域上海段的河口湿地。

1993—2017年的总人口、农林牧渔业从业人口、其他人口(除了湿地农从业人口之外的人口)、商品零售价指数、GDP、人均GDP、农林牧渔业总产值数据来源于中国知网上海市经济社会发展统计数据库(http://data.cnki.net/area/Yearbook/Single/N2019040068?z=D09)。滨海湿地面积变化的主要影响因素是农林牧渔业和人口^[31]。因缺失湿地农业产值数据, 而湿地农业是指在天然湿地基础上改造成以稻田、苇塘、鱼塘、小型水库为主体的农、林、牧、副、渔综合发展的人工农业复合生态系统^[32], 因此采用上海市农林牧渔业从业人口年平均产值、农林牧渔业从业人口、农林牧渔业总产值分别近似代替湿地农业从业人口年平均产值、湿地农业从业人口以及湿地农业总产值。

上海市1995—2017年的市辖区旅游总收入数据来源于上海市人民政府年鉴(http://www.shanghai.gov.cn/nw2/nw2314/nw24651/nw45010/index.html)。2000—2017年河流年输沙量数据来源于中国人民共和国水利部《中国河流泥沙公报》(http://www.mwr.gov.cn/sj/tjgb/zghlnsgb/), 1993—1999年大通站的河流年输沙量因无具体统计数据, 故以公报中大通站1953—1999年平均年输沙量4.36亿t代替1993—1999年的河流年输沙量, 依据水库泥沙容重1.5t/m³换算单位^[28]。

2.3 研究方法

上海市湿地ESV的计算方法表达如下：

(1)

式中, T_ESV为湿地ESV, V_p、V_r、V_c分别为湿地生态系统提供的供给服务价值、调节服务价值和文化服务价值。∑V_pi、∑V_ri、∑V_ci则分别表示各服务指标的总价值。

上海市湿地ESV核算模型的部分公式及输入参数参考商慧敏等^[10]、张绪良等^[28]、钟水映等^[24]的核算方法(表 1)。其中为简化供给服务价值的计算, 本研究参考崔丽娟等^[27]的核算方法将湿地的供给服务价值等价为湿地农业总产值和水资源价值。设定湿地农业总产值等于湿地农业从业人口乘以湿地农业从业人口年平均产值。湿地年平均净初级生产力取自然湿地年平均净初级生产力和人工湿地年平均净初级生产力的算术平均值25.58t/hm^2[28]。

2010年的碳税率为702.95元/t^[33]。2002年的单位围垦工程单价为8.00元/t^[28]。2002年单位面积生物多样性维持服务价值为2786.90元hm^-2 a^-1, 1998年单位自然面积湿地氮吸收能力为0.30t/kg, 1998年建设人工湿地吸收单位质量氮成本5750元/t^[28]。2002年单位水资源价值为3.52元/t^[30]。1997年的单位面积科研教育价值为5880元/hm^2[34]。2010年单位释氧价值为352.93元/t^[33]。由于1993—2017年的25年间物价变化较大, 故1993—2017年期间的碳税率、单位围垦工程价格、单位面积湿地生物多样性维持服务价值、单位自然面积湿地氮吸收能力、建设人工湿地吸收单位质量氮成本、单位水资源价值(取单位地表水价值和单位地下水价值的算术平均值)、单位面积科研教育价值以及单位释氧价值根据文献提供的数据, 按照张绪良等^[28]的方法采用全国零售物价总指数进行换算, 保证计算结果的可比性。

湿地光合作用每形成1g有机质需要吸收1.63g二氧化碳量并释放1.19g氧气量^[28]。年单位湿地水面蒸发量为2177.78m³/hm², 单位体积水量转化为蒸汽耗电量约为125kWh^[27]。电价参考国家电网提供的上海市销售电价表, 取0.61元/kWh。单位面积水源涵养量为744.48m³/hm², ^[35]。

湿地生态系统模型方程中既有线性方程, 也存在一些非线性关系。Vensim PLE软件针对非线性关系可以采取表函数的形式输入, 即可以用图表的方式来反映两个变量之间的非线性关系, 当部分自变量缺失时, 系统可以用线性插值法取其近似值, 从而弥补因部分数据缺失而无法估算的不足。表函数的输入数据均来源于已知的统计年鉴数据或经过全国零售物价总指数换算后的数据。市辖区旅游总收入的统计数据缺失1993—1994年的数据, 遥感数据难以有连续年份的数据, 所以在该模型输入数据时市辖区旅游总收入和湿地面积都采取表函数的形式输入。此外, 为了更好地反映自变量和因变量之间的特殊非线性关系, 在该模型中, 其他人口增长率、湿地农业从业人口增长率、湿地农业从业人口年平均产值、GDP、湿地面积增长率、单位面积湿地生物多样性维持服务价值、单位水资源价值、碳税率、单位释氧价值、年单位自然面积湿地氮吸收能力、建设人工湿地吸收单位质量氮成本、单位围垦工程价格以及河流年输沙量共13个变量也采用表函数的方式建模。表函数的运用可以使模型更加接近实际情况, 减小实际值与模拟值之间的误差。现阶段生态系统支持服务核算方法尚未成熟, 且支持服务通常又是其他三类服务的形成条件, 所以学界已达成共识, 不再核算支持服务的价值量, 以避免重复核算, 因而在本文中不涉及生态系统支持服务的计算。根据模型的主要方程(表 1), 在Vensim PLE平台绘制成系统流图(图 1)。输入相应的参数和经验数据, 设定模型的初始年为1993年, 终止年为2017年, 步长为1年。模拟结果使用Origin 2019软件处理。

3 模型检验和情景模拟 3.1 SD模型验证

模型验证的作用一方面通过模型检验评估其是否能够准确反映现实状态, 另一方面通过灵敏度分析评估模型系统在运行中是否具有稳定的置信度。

3.1.1 模型检验

本研究从人口因素子系统、供给服务子系统、文化服务子系统、调节服务子系统、以及经济因素子系统中共选取10个变量, 分别为总人口、河流年输沙量、湿地农业从业人口、其他人口、湿地农业从业人口年平均产值、人均GDP、GDP、湿地农业总产值、市辖区旅游总收入以及湿地面积。以1993—2017年间上海市统计年鉴、中国河流泥沙公报以及中国统计年鉴等数据为基础, 对10个变量的模型模拟值与历史实际值进行相对误差检验, 因湿地遥感数据没有连续年份的数据, 故湿地面积只采用1995年、2000年、2005年、2010年以及2015年的历史实际值与相应年份的模拟值进行相对误差检验。结果表明(图 3)：1993—2017年间所有检验指标的平均相对误差均小于10%, 模型的仿真值与实际值拟合较好, 因此可以认为SD模型具有可靠性, 可以进行实际仿真操作。

3.1.2 灵敏度分析

灵敏度分析可以评估关键参数的变化对模型输出结果的影响程度。根据顾朝林等的方法, 通过增加或减少10%关键参数值来确定其影响程度^[29]。若模型对大部分参数的变化表现不灵敏, 则说明模型具有很好的稳健性。

SD模型的灵敏度分析分别从人口因素子系统、供给服务子系统、文化服务子系统、调节服务子系统、以及经济因素子系统中选取15个变量, 考察这15个变量对上海市湿地生态系统服务的影响程度即变化幅度均值。这15个变量分别为湿地农业从业人口增长率、其他人口增长率、单位围垦工程价格、湿地面积增长率、单位释氧价值、碳税率、单位面积科研教育价值、湿地年平均净初级生产力、单位面积湿地蒸发量、市辖区旅游总收入、年单位面积湿地氮吸收能力、河流年输沙量、单位水资源价值、湿地农业从业人口年平均产值以及单位面积湿地生物多样性维持服务价值。灵敏度检验结果表明(表 2)：除了3个变量(单位面积湿地蒸发量、市辖区旅游总收入以及湿地农业从业人口年平均产值)的灵敏度较高外, 其他12个变量的灵敏度值均未超过10%, 这说明SD模型具有稳健性, 能对实际情况进行可靠的预测。

灵敏度的大小代表着各变量对湿地ESV的影响程度。灵敏度较高的变量既是对系统影响最大的因素, 也是日后影响湿地生态系统服务变化的主要动力。灵敏度最高的4个变量分别为单位面积湿地蒸发量、市辖区旅游总收入、湿地农业从业人口年平均产值以及湿地面积增长率(表 2), 因此可以确定这4个因素是影响上海市湿地生态系统服务的主要因素, 对整个模型的影响最大。其他人口增长率和单位围垦工程价格这两个变量的灵敏度最小, 即这两个变量的波动对上海市湿地生态系统服务的影响甚微。

3.2 情景模拟 3.2.1 情景模拟参数设定

2017年12月上海市出台《上海市湿地保护修复制度实施方案》(简称《上海实施方案》), 该方案明确, 到2020年, 上海湿地面积不少于46.4×10⁴hm²。利用Vensim PLE软件进行仿真模拟, 设置单位时限为1年, 以1993年为起始年, 终止年为2025年, 进行为期33年的模拟, 初始参数值为1993年所对应的实际值。模型仿真的时长涵盖“十四五”规划期(2021—2025年), 可以为“十四五”规划期的湿地保护决策提供参考。

本文研究湿地的变化对上海市湿地生态系统服务的影响, 上海市2017年湿地面积年增长率为-1.30%, 因此通过设置不同的湿地面积增长率来分析不同情景下上海市湿地生态系统服务变化趋势。在4种情景中, 2018—2025年的湿地面积增长率分别设置为：将维持现状背景下的发展模式定义为情景1, 即湿地面积增长率为-1.30%;将优先发展经济, 湿地总面积维持不变的模式设定为情景2, 即湿地面积增长率为0;假设从2017年后开始逐渐保护湿地生态环境, 划定湿地面积生态红线, 将这种湿地面积逐渐扩大的发展模式设定为情景3, 设定湿地面积增长率为2017年湿地面积增长率的绝对值, 即湿地面积增长率为1.30%;若在情景3中的模式下, 每年湿地增长的面积不足以弥补因自然因素损失的湿地, 需要继续加大湿地恢复力度的模式设定为情景4, 设定湿地面积增长率为2017年湿地面积增长率的绝对值的2.5倍, 即湿地面积增长率为3.50%。

将设定的参数输入模型, 根据表 3中的结果可以看出, 到2020年底, 按照情景1、情景2、情景3、情景4预测, 上海市湿地面积将分别为40.44×10⁴、42.04×10⁴、43.70×10⁴、46.61×10⁴hm²。4种情景中, 只有按照情景4的发展模式, 上海市湿地面积在2020年才能达到《上海实施方案》中的要求, 即2018—2020年年湿地面积增长率设定在3.50%。在情景1和情景2中, ESV和人均ESV均呈下降趋势, 即维持湿地面积不变或减少湿地面积都会导致ESV和人均ESV的下降, 但湿地面积不管是逐年增长、逐年减少还是年际间保持不变, 供给服务价值并未增长, 反而逐年下降。1993—2025年期间, 在生态系统服务构成中, 随着上海市湿地面积逐年减少(图 2), 调节服务价值在ESV中的比例总体上呈下降趋势, 供给服务价值在ESV中的比例总体上呈下降趋势, 而文化服务价值在ESV中的比例总体上呈上升趋势(图 4)。

3.2.2 情景模拟结果及分析

综合4种情景下的模拟结果, 可以发现, 在2017—2025年上海市湿地生态系统服务构成中, 除本文未涉及的生态系统支持服务外, 文化服务价值最高(占总价值的64%以上), 研究时段内呈增长趋势; 调节服务价值次之(占ESV的30%以上), 呈增长趋势; 供给服务价值最低(占ESV的6%以下), 呈下降趋势。依据中国知网上海市经济社会发展统计数据, 在1993—2017上海市三大产业产值中, 第三产业产值增长最快, 第一产业最慢, 到2017年, 第三产业产值占比最大(>69%)。这是上海市湿地生态系统文化服务价值占比最高且不断增长、供给服务价值占比最低且不断下降的原因之一。

1993—2025年, ESV/GDP的值总体上呈下降趋势(图 4)。尽管生态系统服务的估值与经济活动或GDP之间没有必要的关系, 但ESV/GDP的值对于理解生态系统对人类的福祉非常有用, 2010年全国人均ESV为28569元(以2010年美元平均汇率6.77计算), ESV/GDP为0.87^[26]。2017年—2025年, ESV/GDP的值除了在情景1和情景2中逐年下降外(表 3), 其他两种情景中ESV/GDP的值均逐年增大, 即ESV占上海市GDP的份额逐年增大。

到2025年, 情景1、情景2、情景3、情景4中上海市湿地ESV/GDP的值分别为0.20、0.22、0.24以及0.25, 相比2010年全国ESV/GDP的值, 上海市ESV/GDP值远低于全国平均水平。在上海市今后的发展中, 应该逐步增加湿地面积, 恢复湿地的生态系统服务功能, 朝着情景2, 3的模式发展。此外, 模型也可根据ESV/GDP的值来判断是否应设置适当的湿地保护系数来促进湿地面积的增长。例如可以在Vensim PLE软件中将“湿地保护系数”设置为“IF THEN ELSE((“生态系统服务价值/国内生产总值”≥1), (1), (1.2))”, 这表示如果ESV/GDP < 1时, 系统会将“湿地保护系数”设定为1.2, 即需要加大湿地保护力度; 若ESV/GDP>1时, 系统则将“湿地保护系数”设定为1, 即维持湿地面积现状不变。

4 讨论

近30年来, 上海市经济发展突飞猛进, 城市化不断推进。在经济发展的过程中上海市湿地面积不断被开垦为其他的用地, 以满足经济发展的需求。湿地的ESV日益受到人们的重视, ESV受自然因素和人为因素的影响, 各种影响因素动态变化, 这使得在实际估算模拟当中, 模拟方程往往呈非线性。而vensim PLE中的表函数的出现, 使得非线性方程的使用成为现实, 它允许用户通过列举自变量值和因变量值的方式来构建模型。且当某个自变量为非已知点或数据缺失时, 可以用线性插值的方法取得其近似值, 这样就能让模型更加接近实际情况, 同时, 表函数的使用可以让复杂的生态系统服务核算方程变得更为简便。关于ESV的核算, 现阶段的研究主要采用价值当量因子方法计算, 但这种方法存在依赖主观判断的风险, SD模型在ESV的核算方面, 多采用较为精细的计算方式, 在很大程度上减少了依赖主观判断的风险。

目前, ESV的评估众多, 但往往都是基于某一时期的估算, 对于ESV的变化难以持续模拟。这使得决策者难以通过某一时期或几个时期的ESV来判断该地区的生态系统服务变化趋势。SD模型可以同时实现多变量多年份连续动态模拟, 这使得ESV的变化趋势变得更为直观, 也可量化ESV构成之间的变化, 具有实际参考价值。本文通过设定不同的湿地面积增长率分别模拟上海市在采取不同的政策下湿地ESV未来发展情况, 这种模拟可以为划定生态红线以及其他政策的制定提供超前评估, 预防政府决策失误。

“湿地保护系数”采用“IF THEN ELSE”函数, 可以理解为模型能够依据期望设定ESV/GDP的阈值, 再把实际的ESV/GDP值与期望阈值作比较, 根据比较结果再采取不同的湿地保护系数, 这样就使得这个函数实际上具有了两种含义。它能提前发现未来ESV/GDP的变化, 从而为构建生态安全格局提供参考依据。

本文所构建的SD模型虽然部分参数采用了上海市尺度之外的参数, 但理论上各要素之间的相互关系可在本文所确定的概念框架之上通过结合研究区的特点而灵活调整参数, 从而使得模型更加贴近实际情况。此外, 本文采用全国零售物价总指数对不同年份的单位价值量进行换算, 从而减小通货膨胀等因素所带来的估算偏差, 使得模拟结果更具有时效性。但本文所构建的SD模型仍有待将土地因素、气候因素、社会经济情况等因素嵌入模型之中, 并尝试将其量化, 以便在其他尺度的区域进行推广。

5 结论

本文通过构建SD模型对上海市湿地ESV动态变化过程进行仿真模拟, 可以得出如下结论：

(1) 基于1993—2017年的经验数据, 通过对该模型进行系统存流量检验和系统灵敏度分析, 显示该模型有效且具有实际仿真的可操作性, 可根据政策要求设定参数, 对上海市湿地ESV变化进行模拟与预测。

(2) 4种基于不同的“湿地面积增长率”情景模拟显示了未来经济发展在不同增速条件下的湿地面积及ESV变化。按此4种情景预测, 到2020年底上海市湿地面积将分别为40.44×10⁴hm²、42.04×10⁴hm²、43.70×10⁴hm²和46.61×10⁴hm²。无论哪种情景, 在湿地生态系统服务构成中, 文化服务价值最高, 约占总价值的64%以上, 研究时段内呈增长趋势; 调节服务价值次之, 呈增长趋势; 供给服务价值最低, 呈下降趋势。在影响上海市湿地ESV变化的因素中, 最重要的4个影响因素分别为单位面积湿地蒸发量、市辖区旅游总收入、湿地农业从业人口年平均产值、湿地面积增长率。对比4种情景, 只有情景4满足当下对于湿地面积的要求, 适合当前的发展方向, 即年湿地面积增长率设定在3.50%。

(3) 到2025年, 即“十四五”规划末期, 按照现行的政策(情景1), 上海市湿地面积将减少至37.92×10⁴hm², 湿地ESV将减少至4711.56×10⁸元, 其中文化服务价值为3085.67×10⁸元, 调节服务价值为1441.27×10⁸元, 供给服务价值为184.62×10⁸元。人均ESV为1.95×10⁴元, 湿地ESV占上海市GDP的20%。