土地利用变化对生态系统服务的影响机制极其复杂且具有较大时空异质性, 不同的研究方法或研究区域会使研究结果差别明显。因此, 需要对这些研究案例进行整合分析, 从而找到普遍结论并发现差异。在此情形下, 综合分析研究方法得到广泛使用, 尤其是Meta-analysis (整合分析)^[1]。

Meta-analysis又称元分析或荟萃分析, 相对于传统的文献综述, Meta-analysis扩大了样本例数, 提高统计学检验效能, 改善对效应量的估计, 因而具有全面、系统、客观和定量的特点, 在自然和社会科学各个学科得到广泛应用^[2]。诸多生态系统服务价值研究中, 最有影响力的是Daily^[3]编著的《Nature’s service: societal dependence on natural ecosystem》和Costanza等^[4]发表的“The value of the world’s ecosystem services and natural capital”。Bateman等^[5]从经济学视角研究土地利用变化引起的生态系统服务变动, 在土地利用规划方面提出相应的政策建议, 拓展了Meta-analysis在土地管理和政策制定方面的应用。Brander等^[6]应用meta-analysis研究湿地与其生态系统服务价值之间的关系, 使用价值转移研究范式, 建立生态系统服务价值回归方程, 并对其可行性进行探讨。谢高地等^[7]在Costanza等价值评价结果的基础上, 对国内200多位学者进行问卷调查, 获得中国陆地生态系统服务价值当量因子表。

中国生态系统服务价值评估工作起源于20世纪80年代起步的森林资源价值核算^[8-9]。欧阳志云等^[10-11]评述了生态系统服务及其生态经济价值评估方法, 并测算了中国陆地生态系统服务价值。陈仲新等^[12]按照Costanza等的单位价值量, 将中国生态系统服务的价值与世界生态系统服务总价值进行对比。宗跃光等^[13-14]根据不同土地利用类型, 计算区域生态系统服务价值, 并将Costanza等单纯的自然资本价值测算拓展到自然资本、经济资本和社会资本的综合测算。此外, 还有大量研究侧重于全国草地^[15]、森林^[16]、农田^[17]等不同类型生态系统服务价值的估算。

已有的生态系统服务价值研究大多只着眼于特定区域, 未考虑地理区位和社会经济等因素对生态系统服务价值的影响, 评估结果难以进行横向比较, 从而限制了其政策应用。针对这种情况, 本文选择各土地利用类型提供的生态系统服务价值作为因变量, 年代、地理分区、评估方法、人均GDP和人口密度等作为解释变量, 研究这些因素对各类土地利用生态系统服务价值的影响程度, 同时建立预测模型, 实现价值转移功能, 将其推广到不同地区以进行生态系统服务价值评估^[18]。

1 数据收集与样本遴选

为了使用Meta-analysis核算生态系统服务价值, 应尽可能收集多的研究案例。为此, 本研究共查阅国内327篇文献。其中, 许多案例直接引用了谢高地等^[7]的“中国生态系统单位面积生态服务价值当量因子”作为测算依据, 为了排除冗余, 本文将引用谢高地的文章作为一个方法案例。据此, 共筛选出38篇文献、197个数据作为案例数据, 样本的地理分布遍及中国23个省份(图 1)。

本研究将土地利用类型划分为6类, 即耕地、林地、草地、园地、水域和未利用地。其中, 耕地、林地、草地、园地分别与2007年国家标准《土地利用现状分类》中的地类相对应; 水域对应于国家分类中的水域及水利设施用地, 文献中没有细分的湿地类型也一并归入此类; 未利用地包括空闲地、沙地和裸地, 部分文献中的未利用地包括冰川及永久积雪, 导致其生态系统服务价值估算量较高。

在评估方法上, 引用Costanza估算值的称为Costanza系数法, 引用谢高地估算值的称为谢高地系数法。这些研究利用了价值转移, 而不是初级评估手段; 利用当地粮食价格和谢高地当量因子计算当地生态系统服务价值量的称为公式法; 市场价值法、机会成本法、影子工程法、替代花费法等归为一类, 称为其他方法。

一般认为, 单位生态系统服务价值随土地利用类型、评估方法、年代等变化较大^[19](图 2)。Bran-der^[6]在对世界范围的湿地生态系统服务价值研究中发现, 年代、地域和方法的选择对生态系统服务价值影响很大。事实上, 生态系统服务价值量也受其所在地区的社会经济特性制约(图 2), 但这方面的研究略有欠缺。本研究中的影响因素包括这些社会经济因子, 然而一些研究案例中并不包含诸如人口密度和人均GDP等社会经济信息。为此, 本研究从各地统计年鉴中补录了上述数据。

从图 2可以看出, 生态系统服务价值由高到低的顺序分别是水域、林地、园地、耕地、草地和未利用地。

为了减少整合分析中的不确定性, 对所有样本中的各个土地利用类型的价值进行奇异值检查(图 3)。从图 3看出, 草地、水域和未利用地有1~2个奇异值。对这些奇异值溯源可以发现, 草地中第10个样本, 由造林成本法核算的固定CO₂价值较高, 由机会成本法核算的土壤保持价值也较高, 这是方法选择引起的生态系统服务价值偏高。水域中, 第13个样本评估得出的长白山天池的水源涵养、调节服务价值很高。未利用地中, 第12和第14个样本将冰川积雪包含到未利用地中, 使难利用地与水体的面积比例约为3:1, 导致其水源涵养和保持、废物处理价值很高, 取它们的加权平均值作为未利用地的生态系统服务价值后, 整体上仍偏高。

2 Meta回归分析 2.1 Meta回归模型

为了获取相关因素之间的定量关系, 本文使用Meta回归模型来同时评估所有潜在相关因素的相对重要性^[20]。根据土地利用类型划分方案, 共建立12个回归方程式, 其中, 解释模型和预测模型各6个。两类模型的回归参数可能相差甚远且符号不一致, 因为预测模型仅用来预测, 不具备解释能力; 同理, 解释模型仅用来解释, 不具备预测能力^[21]。回归方程中, Y_n(n取1, 2, …, 6)分别表示耕地、林地、草地、园地、水域和未利用地的单位面积生态系统服务价值(元/(hm² ·年))。解释变量共有5个, 其中, X_s表示年代, X_p表示人均GDP, X_e表示人口密度, X_f表示方法(X_f1: Costanza系数法; X_f2:谢高地系数法; X_f3:公式法; X_f4:其他方法), X_i表示区位(X_i1:东部; X_i2:中部; X_i3:西部)。回归模型如下:

${Y_n} = {b_s}{X_s} + {b_p}{X_p} + {b_e}{X_e} + {b_f}{X_f} + {b_i}{X_i},$

(1)

式(1)中, 年代、人均GDP和人口密度都是数值型变量, 方法和分区则是类别型变量。在解释模型中, 为了使回归系数具有可比性, 对变量都用标准计分(standard score)方法^[21]做了标准化处理。由于引入全部分类变量, 默认常数为0, b是相应变量的回归系数, 检验排除了多重共线性, DW值通过检验。

2.2 土地利用类型回归结果

2.2.1 耕地

运行土地利用Meta回归模型, 得到耕地回归方程如下。

1)解释模型:

Y₁=0.142X_s+0.063X_p+0.019X_e-1.759X_f1-0.611X_f2-0.669X_f3+0.052X_i1+0.535X_i2-0.586X_i3。

2)预测模型:

Y₁′=591.98X_s+0.01X_p+0.27X_e+4796.15X_f2+4555.5X_f3+7349.73X_f4-499.58X_i1+1516.51X_i2-3165.47X_i3。

式中, 解释模型R²=0.523, 预测模型R²=0.877, 调整的R²=0.774。与数据预分析得出的结论一致, 表明人均GDP和年代对耕地的生态系统服务价值有正向作用。人口密度与其关系不显著, 评估方法对耕地生态系统服务价值有负向作用。

2.2.2 林地

运行土地利用Meta回归模型, 得到林地回归方程如下。

1)解释模型:

Y₂=0.191X_s+0.224X_p+0.239X_e-1.175X_f1-0.878X_f2-0.366X_f3+0.100X_i1+0.149X_i2-0.816X_i3。

2)预测模型:

Y₂′=1688.13X_s+0.08X_p+7.32X_e+2639.49X_f2+7166.53X_f3+10412.73X_f4-4374.81X_i1+3935.40X_i2-1974.43X_i3。

式中, 解释模型R²=0.540, 预测模型R²=0.927, 调整的R²=0.838。人均GDP、人口密度、年代都对林地的生态系统服务价值有正向作用, 地域和评估方法与林地生态系统服务价值呈负相关。

2.2.3 草地

运行土地利用Meta回归模型, 得到草地回归方程如下。

1)解释模型:

Y₃=0.316X_s+0.327X_p+0.1X_e-0.932X_f1-2.240X_f2-2.086X_f3+0.804X_i1+0.222X_i2-0.169X_i3。

2)预测模型:

Y₃′=1620.31X_s+0.06X_p+1.9X_e-4782.51X_f1-11498.58X_f2-10706.65X_f3+9011.88X_i1+6022.24X_i2+5751.55X_i3。

式中, 解释模型R²=0.653, 预测模型R²=0.879, 调整的R²=0.722。人均GDP、人口密度、年代和地域都对草地的生态系统服务价值有正向作用, 而方法皆为负影响, 说明文中列举的方法计算草地生态系统服务价值有偏低的可能。

2.2.4 园地

运行土地利用Meta回归模型, 得到园地回归方程如下。

1)解释模型:

Y₄=0.239X_s+0.435X_p+0.371X_e-1.223X_f2-0.449X_f3+0.281X_f4-0.463X_i1-0.877X_i3,

2)预测模型:

Y₄′=1556.52X_s+0.09X_p+7.82X_e-120.95X_f1+4911.71X_f2+9663.1X_f3+0.09X_f4-3017.13X_i1-5704.67X_i3,

式中, 解释模型R²=0.648, 预测模型R²=0.938, 调整的R²=0.694。人均GDP、人口密度和年代对园地的生态系统服务价值有正向作用。地域上, 东部和西部与园地生态系统服务价值呈负相关。方法上, 除“其他方法”对园地生态系统服务价值是正影响外, 其余两种方法皆为负影响(园地中没有用Costanza系数评估的样本案例), 说明用谢高地系数、公式法计算园地生态系统服务价值有偏低的可能性。

2.2.5 水域

运行土地利用Meta回归模型, 得到水域回归方程如下。

1)解释模型:

Y₅=-0.111X_s-0.16X_p+0.409X_e-2.669X_f2-2.136X_f3+2.775X_i1+2.719X_i2+2.584X_i3。

2)预测模型:

Y₅′=-2398.31X_s-0.14X_p+34.03X_e-57696.29X_f2-46167.17X_f3+120986.1X_i1+119763.87X_i2+116841.87X_i3。

式中, 解释模型R²=0.951, 预测模型R²=0.813, 调整的R²=0.673。人均GDP、年代和方法与水域生态系统服务价值呈负相关。随着时间的推进和经济的发展, 人们在源源不断地获取水源涵养与洪水控制的生态系统服务的同时, 欠缺对水域合理的保护; 文中列举的评估方法可能低估了水域的生态系统服务价值。人口密度和地域则与其有正向作用。

2.2.6 未利用地

运行土地利用Meta回归模型, 得到未利用地回归方程如下。

1)解释模型:

Y₆=0.075X_s-0.236X_p-0.253X_e-0.348X_f1-0.594X_f2-0.064X_f4+0.372X_i1-0.512X_i2-0.265X_i3,

2)预测模型:

Y₆′=188.12X_s-0.03X_p-2.36X_e-713.33X_f1-1331.65X_f2-+161.1X_f3+1.491X_f4+3318.13X_i1+1096.55X_i2+1716.84X_i3。

式中, 解释模型R²=0.218, 预测模型R²=0.349。人均GDP、人口密度和方法与未利用地生态系统服务价值呈负相关, 经济发达或者人口密集的地方, 未利用地提供的生态系统服务种类减少, 质量降低。年代则对其有正向作用。地域上, 东部与未利用地生态系统服务价值呈正相关, 西部和中部与其呈负相关。

为了更清晰地反映各土地利用类型提供的生态系统服务价值与其影响因素之间的关系, 表 1将各地类的Meta回归结果进行汇总。

3 价值转移 3.1 可行性分析

将Meta回归分析结果用于价值转移, 即利用一个地区的自然和社会经济属性数据预测和估计当地土地利用类型生态系统服务价值^[19]。价值转移有直接价值转移和价值转移方程^[4]两种方法。前者简单地将某几个样本区的研究结果应用到研究区, 理想情况下, 研究区的特征参数应当与样本区相似或可调整到相似^[22]。后者基于已知的属性数据, 利用价值回归方程核算研究区的生态系统服务价值, 一般采用meta回归方法^[23]。应用价值转移方程需满足以下条件^[24]: 1)有一个评估方程式, 可将生态系统服务的价值量与研究区的时间、空间等属性联系起来; 2)研究区与样本区之间的差异可以通过价值体现出来; 3)价值在时间变化下是稳定的, 或者考虑了时间的影响; 4)所取用样本的核算价值是准确的。

一般情况下, 价值转移方程的效果比直接价值转移稍有优势。原因如下:第一, 价值转移方程中包含了更多的信息量; 第二, 价值转移方程控制了样本估算方法之间的差别; 第三, 类别变量也可量化, 更精确地核算价值^[19]。

在价值转移过程中, 可能会出现一些“转移误差”。首先, 获取评估方程式中需要的研究区自然社会属性数据有一定的不确定性; 其次, 有些地方土地利用类型和生态系统服务类型划分界限不明确, 因此类别变量对社会经济因素的评估不能准确捕捉真实的变化量^[24]。本研究也不可避免地面临这样的问题, 应对措施是采用的社会经济因素尽量集中在数值变量上, 例如人均GDP和人口密度等。这些因素与生态系统服务价值越相关, 就越能增加价值转移评估的正确性; 最后, 回归所采用的案例样本的核算误差会导致价值转移的误差。

在实现价值转移之前, 先分析模型预测的效能。这里引入一个检验指标——平均完全百分误差(MAPE)^[4]:

$MAPE = |({Y_{obs}} - {Y_{pred}})/{Y_{obs}}|{\rm{ }},$

(2)

其中, Y_obs表示案例中直接给出的生态系统服务价值, Y_pred表示使用回归方程核算的生态系统服务价值。下面以耕地和林地为例给出误差示意图(图 4)。经过计算, 在38个案例中, 平均转移误差为5%, 68%以上的样本数据转移误差低于10%, MAPE在可接受范围内^[18], 因此可以进行价值转移。

3.2 基于价值转移的京津冀地区土地利用类型生态系统服务价值量

京津冀地区包括北京市、天津市以及河北省的唐山、石家庄、秦皇岛、廊坊、保定、沧州、衡水、邢台、邯郸、张家口和承德等市(图 5)。使用2012和2013年京津冀13个城市的人均GDP和人口密度, 代入回归方程。因不同回归方程避免多重共线性所舍弃的方法变量不同^[25], 且此处是同等条件下分析比较, 所以一致舍弃方法选择类变量。价值转移结果如表 2和3所示。将价值转移后的京津冀地区土地利用类型生态系统服务价值进行空间化表达, 并分成4个等级(图 6)。

与已有的研究对比, 本文经价值转移得到的北京市生态系统服务价值比赖瑾瑾等^[26]在顺义区的出的结果略高, 与两个研究的时空尺度不一致有关。与马程等^[27]得出的京津冀地区不同土地利用类型生态系统服务价值相比, 耕地、水域和未利用地的生态系统服务价值相近, 而林地和草地则较之略高, 与时间尺度与方法选择不同有关。实际上, 北京市作为一国之首都, 人口密度大, 其担负的提供生态系统服务责任重大, 而北京市及周围地区的生态系统服务价值偏高也符合实际情况。其未利用地价值偏低与北京及周边的未利用地得到大力开发有关。

将土地利用类型按照生态系统服务价值高低排序, 依次为:水域、林地、草地、园地、耕地和未利用地。研究发现, 京津冀地区不同市区的不同土地利用类型提供的生态系统服务价值高低不同, 按照生态系统服务价值由高到低的顺序汇总到表 4。

耕地和林地的生态系统服务价值空间格局较相似, 前3个城市分别是北京、天津和唐山, 后3个城市分别是衡水、承德和张家口。天津的草地生态系统服务价值最高, 剩下的城市排名与耕地一致。城市按园地生态系统服务价值量高低排名与林地完全一致。北京、天津和邯郸的水域生态系统服务价值最高, 秦皇岛、张家口和承德的最低。按照未利用地生态系统服务价值高低与按照其他地类排序得到的结果完全相反, 最高的是张家口、承德和衡水, 唐山、天津和北京最低。总生态系统服务价值最高的是北京, 其次是天津、邯郸, 张家口和承德最低(图 7)。

4 结论与讨论

本文收集了近年中国土地利用类型生态系统服务价值评估的相关案例, 建立Meta回归模型, 通过价值转移计算京津冀地区2012和2013年土地利用类型的生态系统服务价值。研究结果如下。

1)从前研究常被忽视的社会经济因素, 如人均GDP和人口密度等, 在解释生态系统服务价值变化中起到了重要作用。除未利用地外, 人口密度与其他地类的生态系统服务价值都呈正相关, 人均GDP与水域和未利用地的生态系统服务价值都呈负相关; 除水域外, 年代与各地类生态系统服务价值呈正相关, 表明随着时间的推移, 整体上各地类生态系统服务价值量有所提高, 水域价值量则略有下降。

2)在评估方法上, 应用Costanza系数会高估耕地的生态系统服务价值, 而低估草地和林地的生态系统服务价值。除耕地外, 谢高地提出的系数及公式与生态系统服务价值呈负相关, “其他方法”则高估了园地的生态系统服务价值。

3)利用平均完全百分误差(MAPE)评估了Meta回归模型价值转移的可行性。价值转移在京津冀地区各城市的结果表明:耕地、林地、草地和园地的生态系统服务价值空间格局较相似, 前三名分别是北京、天津和唐山, 后三名分别是衡水、承德和张家口; 水域生态系统服务价值最高的是北京、天津和邯郸, 最低的是秦皇岛、承德和张家口; 未利用地生态系统服务价值的空间格局和其他地类完全相反, 张家口、承德和衡水最高, 唐山、天津和北京最低。

基于Meta回归模型的生态系统服务价值方程纳入了社会经济变量, 给出更全面的核算结果。价值转移过程节省了资金和时间的投入, 较为高效、准确地核算生态系统服务价值量, 进而为政策制定提供支撑。在应用Meta回归模型的价值转移方程时, 应当格外注意所采用的社会经济数据的准确性^[1]。未来的研究工作中, 需要范围更广、更高质量的初始评估研究进行Meta回归, 进一步提高结果精确性。