防风固沙型重点生态功能区防风固沙服务的评估与受益区识别

徐 洁, 肖 玉, 谢高地, 王洋洋, 江 源, 陈文辉

1 北京林业大学自然保护区学院, 北京 100083 2 中国科学院地理科学与资源研究所, 北京 100101 3 中国科学院大学资源与环境学院, 北京 100049 4 北京师范大学地理科学学部, 北京 100875 5 浙江农林大学信息工程学院, 临安 311300 6 浙江省林业智能监测与信息技术研究重点实验室, 临安 311300

2010年12月国务院正式发布关于印发全国主体功能区规划的通知,按照生态脆弱性和生态重要性两个指标,划定了涵养水源、保持水土、防风固沙、生物多样性维护4种类型25个国家重点生态功能区。其中,防风固沙型重点生态功能区主要指沙漠化敏感性高、土地沙化严重、沙尘暴频发并影响较大范围的区域,其沙尘暴影响范围涉及华北、东北及周边的国家和地区[1-2]。因此,防风固沙服务的跨区流动效应显著,防风固沙型重点生态功能区的防风固沙服务功能提升对中国和东北亚地区的人类福利和社会经济发展具有重要的保障作用[3]。

防风固沙服务是生态系统植被对风沙的抑制和固定作用[4],是风蚀地区自然生态系统提供的一项重要防护型服务,为区域生产生活的可持续发展创造条件。退耕(牧)还林(草)、控制放牧强度、禁止开矿等生态工程建设都是主要的防风固沙措施。这些措施的实行限制了当地社会经济的发展,降低了当地居民的经济收入,增加了区域的发展机会成本,实际上这部分损失应该由防风固沙服务的受益区来补偿。但是已有的补偿措施主要集中在中央政府转移支付等纵向的生态补偿,如京津风沙源治理工程、三北防护林建设等,缺少区域间的生态补偿政策。区域间生态补偿的关键在于建立补偿方(生态系统服务受益区)与受偿方(生态系统服务供给区)的空间对应关系。但是已有的防风固沙型重点生态功能区研究主要集中在防风固沙服务功能的评估[5- 8]、生态系统变化状况分析[9]、土地利用格局演变特征的分析[10]、生态安全评价[11]等方面。很少涉及防风固沙服务的跨区流动分析,因而不能界定防风固沙服务的受益区范围,难以建立防风固沙服务功能与受益区之间的时空对应关系。

防风固沙服务功能的评估是防风固沙服务流动分析的基础,可以通过风力侵蚀模型计算植被引起的风蚀减少量获得。目前修正风蚀方程(Revised Wind Erosion Equation,RWEQ)得到广泛应用,其综合考虑了气候条件、植被覆盖状况、土壤可蚀性、土壤结皮、地表粗糙度、社会驱动力等因素,具有简单的模拟过程和较少的输入参数,可用于识别风蚀区域、评估风蚀风险[12]。国内学者通过参数修正和公式调整将该模型应用到我国的风蚀评估工作中[13- 17]。防风固沙服务的受益区识别可以通过模拟没有防风固沙服务条件下(裸露地表下)沙尘的扩散影响范围获得,防风固沙服务的存在避免了该部分沙尘由沙尘源区(防风固沙服务供给区)流向沙尘沉降区(防风固沙服务受益区),受益区沙尘沉降量的减少是防风固沙服务效应的最直观体现。已有的沙尘流动路径模拟模型多为空气质量模型,主要有天气预报模式(The Weather Research and Forecasting Model,WRF)[18]、全球环境多尺度模型(Global Environmental Multiscale,GEM)[19]、通用多尺度空气质量模型(Community Multiscale Air Quality,CMAQ)[20]、混合单粒子拉格朗日积分轨迹模型(Hybrid Single-Particle Lagrangian Integrated Trajectory,HYSPLIT)[21]等。其中WRF、GEM、CMAQ模型模拟精度较高,但是需要气象学、流体力学等专业知识背景,且对气象、沙尘源等数据精度的要求过于苛刻,可操作性低,不适用于生态系统服务扩散的跨学科研究领域。HYSPLIT模型可操作性较强,是由美国国家海洋和大气管理局(the National Oceanic and Atmospheric Administration,NOAA)的空气研究实验室(Air Research Laboratory,ARL)研发的一种用于计算和分析大气污染物输送、扩散轨迹的专业模型,已经被广泛地应用于多种污染物的扩散研究[22- 24],包括沙尘[25- 27]、苯[28-29]、火山喷发[30]、PM2.5[31-32]、霾污染[33]等的扩散模拟。

为了消除转移支付政策对防风固沙型重点生态功能区的影响,更为客观地评价防风固沙型重点生态功能区防风固沙服务的跨区流动效应,本研究选取主体功能区规划正式实施之前的2010年为节点,基于RWEQ模型评估2010年防风固沙型重点生态功能区的防风固沙量,基于HYSPLIT模型模拟各个防风固沙型重点生态功能区的防风固沙服务流动路径,识别出各个功能区防风固沙服务的受益区范围,在此基础上得到整个防风固沙型重点生态功能区的防风固沙服务流动路径和受益区范围,建立防风固沙型重点生态功能区的防风固沙服务功能量及其受益区的时空联系,为防风固沙型重点生态功能区的区域间生态补偿政策制定提供科学的参考依据。

1 研究区与研究方法

1.1 研究区

防风固沙型重点生态功能区包括塔里木河荒漠化防治生态功能区(38.14×104km2,以下简称塔里木)、阿尔金草原荒漠化防治生态功能区(33.70×104km2,以下简称阿尔金)、阴山北麓草原生态功能区(9.68×104km2,以下简称阴山)、浑善达克沙漠化防治生态功能区(16.74×104km2,以下简称浑善达克)、科尔沁草原生态功能区(11.12×104km2，以下简称科尔沁)、呼伦贝尔草原草甸生态功能区(4.49×104km2,以下简称呼伦贝尔),均位于我国北方的干旱、半干旱地区(图1),风沙活动频繁,属于沙尘活动的沙源区,易对区域生态环境产生较大影响。

图1 防风固沙型重点生态功能区分布及其2010年土地覆被类型Fig.1 Distribution of the national key ecological function zone of windbreak and sand fixation type in China and its land cover types in 2010防风固沙型重点生态功能区代码说明：TLM：塔里木河荒漠化防治生态功能区;AEJ：阿尔金草原荒漠化防治生态功能区;YS：阴山北麓草原生态功能区;HSDK：浑善达克沙漠化防治生态功能区;KEQ：科尔沁草原生态功能区;HLBR：呼伦贝尔草原草甸生态功能区

1.2 研究方法1.2.1 防风固沙服务功能量的评估

当风经过地表时,会受到来自植被的阻挡,使得风力减弱,风蚀量降低,由植被作用引起的风蚀减小量为防风固沙服务功能量,用潜在风蚀量与实际风蚀量的差值表示。潜在风蚀量指没有植被的裸土条件下的土壤风蚀量,实际风蚀量指现实中地表植被覆盖条件下的土壤风蚀量。

(1)

Qrmax=109.8·(WF·EF·SCF·K′)

(2)

Sr=150.71·(WF·EF·SCF·K′)-0.3711

(3)

(4)

Qmax=109.8·(WF·EF·SCF·K′·C)

(5)

S=150.71·(WF·EF·SCF·K′·C)-0.3711

(6)

SR=SLr-SL

(7)

式中,SLr表示潜在风蚀量(kg/m2)；Qrmax表示潜在风力的最大输沙能力(kg/m)；Sr表示潜在关键地块长度(m)；SL表示实际风蚀量(kg/m2)；Qmax表示风力的最大输沙能力(kg/m)；S表示关键地块长度(m)；SR表示防风固沙量(kg/m2)；z表示所计算的下风向距离(m),本次计算取50 m；WF为气候因子(kg/m)；EF为土壤可蚀性百分比(%)；SCF为土壤结皮因子(无量纲)；K′表示土壤糙度因子(无量纲)；C表示植被因子(无量纲)。

其中,气候因子WF表征了在考虑降雨、温度、日照及雪盖等因素下风力对土壤颗粒的搬运能力,其计算公式如下：

WF=Wf·(ρ/g)·SW·SD

(8)

Wf=u2·(u2-u1)2·Nd

(9)

式中,WF为气象因子(kg/m)；Wf为风力因子(m/s3)；g为重力加速度(m/s2)；ρ为空气密度(kg/m3)；SW为土壤湿度因子(无量纲)；SD为雪盖因子(无积雪覆盖天数/研究总天数,定义雪盖深度<25.4 mm为无积雪覆盖)；u1为2 m处起沙风速,本次计算取5 m/s[13,14,34]；u2为气象站2 m处月均风速值(m/s)；Nd为各月2 m处日均风速大于5 m/s的天数[13]。

土壤湿度因子SW的表达式如下：

(10)

式中,ETo为潜在蒸散发量(mm),利用彭曼公式计算得到,P为降水量(mm),I为灌溉量(mm),本文计算取0[35],Rd为降雨次数,N为试验天数。

土壤可蚀性因子EF是一定土壤理化条件下土壤受风蚀影响的大小,其表达式如下：

(11)

式中sa为土壤砂粒含量(%)；si为土壤粉砂含量(%)；cl为土壤粘粒含量(%)；OM为土壤有机质含量(%)；CaCO3为碳酸钙含量(%),本次计算取0[13]。采用对数正态分布模型进行土壤粒径转化[34]。对于土壤数据库中部分土壤类型的属性值缺失情况,采用相近土壤类型的属性值进行替代(附表)。

土壤表层的坚硬结皮能有效防止风蚀的发生。土壤结皮因子SCF即在一定土壤理化条件下土壤结皮抵抗风蚀能力的大小,其表达式如下：

(12)

式中cl为土壤粘粒含量(%)；OM为土壤有机质含量(%)。

植被覆盖因子C表示一定植被条件对风蚀的抑制程度,其表达式如下：

C=e-0.0483(SC)

(13)

SC=(NDVI-NDVImin)/(NDVImax-NDVImin)

(14)

式中,SC为植被覆盖度(%),NDVI、NDVImax、NDVImin分别代表NDVI实际值、最大值及最小值。

地表糙度K′是地形所引起的地表粗糙程度对风蚀影响的反映,其表达式如下[35]：

K′=cosα

(15)

式中,α为坡度,由1 km的DEM数据经过ArcGIS的slope模块计算得到。

1.2.2防风固沙服务空间流动路径模拟

防风固沙服务是一项防护型服务,假设防风固沙型重点生态功能区内没有植被覆盖或其他固沙措施,全部为沙化土地或沙漠,一旦当地风速高于起沙风速,则沙尘将向下风向扩散,对下风向区域造成损害。因此,为了确定防风固沙型重点生态功能区防风固沙服务的受益区范围,首先需要明确防风固沙型重点生态功能区在风速高于起沙风速时的气团运行轨迹[36]。起沙风速受土壤粒径、土壤水分、地表粗糙度、土壤含盐量等因素影响[37-39]。不同地区土壤沙粒的起沙风速存在较大差异,根据已有文献梳理出各个防风固沙型重点生态功能区代表站点的起沙风速(表1),分别筛选出各站点风速高于起沙风速的时间点,从而模拟相应时间点的气团扩散路径。

表1 防风固沙型重点生态功能区代表站点的平均起沙风速

Table 1 The average threshold wind velocity of the selected stations in the national key ecological function zone of windbreak and sand fixation type

功能区代码Code站点名称Station地理坐标Geographic coordinates平均起沙风速Average threshold wind velocity/(m/s)观测高度Observation height/m位置和文献来源Location and sourceHLBR新巴尔虎右旗(48.67°N, 116.82°E)6.0010—12(气象站风速观测高度)中国北方干旱、半干旱地区[40]KEQ扎鲁特(44.57°N, 120.90°E)HSDK苏尼特左旗(43.85°N, 113.63°E)5.472浑善达克[41]YS乌拉特中旗(41.57°N, 108.52°E)5.650.6内蒙古自治区武川县[42]AEJ若羌(39.03°N, 88.17°E)52塔里木盆地[43]TLM塔什库尔干(37.77°N, 75.23°E)4.62塔中[44]

在风速筛选的基础上,利用HYSPLIT模型模拟沙尘扩散的后向轨迹[21]。HYSPLIT模型使用拉格朗日法计算气团轨迹,假设空气中的粒子随风飘动,那么它的移动轨迹就是其在时间和空间上位置矢量的积分,最终的位置由初始位置和第一次推测位置的平均速率得到。本研究将防风固沙型重点生态功能区对应的6个代表气象站点(表1)分别作为HYSPLIT模型模拟的起始点,模拟各站点风速在高于起沙风速时的气团5日后向扩散路径,之后将各站点路径合并得到防风固沙型重点生态功能区防风固沙服务的空间流动路径。

1.2.3防风固沙服务受益区识别

沙尘天气的减少是防风固沙服务最直接的效益。沙尘天气的危害主要包括[45,46]：(1)风沙流吹蚀或磨蚀耕地土壤,造成根系外露、农作物受损；(2)风沙流搬运沙粒沉积掩埋灌渠、居民区、工矿、道路等基础设施,加剧土地沙化；(3)与沙尘暴相伴的强风破坏公共设施、工农业设施、树木和花果,伤害人畜；(4)沙尘暴导致大气颗粒物浓度增加,对人体健康造成损害；(5)风沙流向海洋输送和沉降,对海洋中营养盐的供给和海洋生态系统造成影响。因此,防风固沙的效益包括：(1)防止人体健康受损；(2)防止交通、通讯、电力、工矿、建筑、灌渠等设施受损；(3)防止耕地受损；(4)防止草场和牲畜受损；(5)防止果园及其他森林受损；(6)防止海洋生态系统破坏；(7)缓解沙漠化的出现与恶化。对应的受益区包括：(1)人类活动区域,如建设用地、耕地、草地、林地、水域与湿地等土地覆被类型；(2)基础设施分布区域,如建设用地、内陆水体、湿地等；(3)耕地；(4)草地；(5)林地；(6)海域；(7)沙地和裸地。

基于HYSPLIT模型模拟生成的多条轨迹,利用HYSPLIT的频率分布模块插值生成1°×1°栅格的轨迹分布频率图,栅格轨迹分布频率计算方法为：

(16)

式中,pi为第i个栅格风沙轨迹分布频率；Li为第i个栅格风沙轨迹通过数；L为起始点的风沙总轨迹数。每个栅格内的数值为通过该栅格的轨迹分布频率。插值范围即为防风固沙型重点功能区防风固沙服务的受益区范围,受益区内的林地、草地、耕地、建设用地、内陆水体、湿地、沙地和裸地等即为防风固沙服务的受益土地覆被类型。

1.3 数据来源

为评估防风固沙服务及其受益区范围,本研究需要以下数据：

(1)气象数据：用于计算土壤湿度和气象因子的气象数据来源于中国气象科学数据共享服务网(http://cdc.cma.gov.cn)的国家台站数据(图2),包括日均风速、最大风速、降水、平均温度、最高温度、最低温度、平均相对湿度、日照时数等数据。用于模拟防风固沙服务前向运行轨迹的数据包括风速数据、NCEP再分析资料。风速数据为2010年1月1日0：00到12月31日18：00的1小时内10分钟最大风速,时间间隔为6小时,数据来源为华云信息技术工程有限公司(http://www.huaxin-hitec.com/)。NCEP再分析资料是从NOAA的ARL获取的2010年1月至2011年1月每日6小时的GBL资料,空间水平分辨率为2.5°×2.5°。(2)土壤数据：来源于第二次全国土地调查中国科学院南京土壤研究所提供的1∶100万土壤数据,用于计算可蚀性因子与土壤结皮因子。(3)雪盖数据来源于旱区寒区科学数据中心(http://westdc.westgis.ac.cn)下载的“中国雪深长时间序列数据集(1979—2016年)”,空间分辨率为0.25°。(4)土地覆被数据：来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)的2010年土地覆被数据,空间分辨率为1 km。(5)数字高程模型(DEM)：来源于中国科学院资源环境科学数据中心(http://www.resdc.cn),空间分辨率为30 m。(6)NDVI数据：来源于中国科学院计算机网络信息中心国际科学数据镜像网站(http://www.gscloud.cn)2010年的月尺度数据,空间分辨率为250 m。(7)防风固沙型重点生态功能区范围：利用1∶100万全国县域数据和2010年《全国主体功能区规划》中国家重点生态功能区名录生成。

图2 防风固沙型重点生态功能区防风固沙服务评估相关国家气象站点的空间分布Fig.2 Spatial distribution of the selected national weather stations used to evaluate the wind erosion prevention service by the ecosystems in the national key ecological function zone of windbreak and sand fixation type in China

2 结果与分析

2.1 防风固沙型重点生态功能区的防风固沙服务功能评估

图3 2010年防风固沙型重点生态功能区防风固沙量空间分布格局Fig.3 The spatial pattern of wind erosion prevented amount of the national key ecological function zone of windbreak and sand fixation type in China in 2010

从各个防风固沙型重点生态功能区的防风固沙总量来看,2010年浑善达克的防风固沙总量最高,可达1.91×1012kg,其次为科尔沁(1.38×1012kg)、塔里木(0.79×1012kg)、阴山(0.73×1012kg)、阿尔金(0.54×1012kg)、呼伦贝尔的防风固沙总量最低,为0.20×1012kg。从单位面积防风固沙量来看,2010年科尔沁的单位面积防风固沙量最高,可达12.51 kg/m2,其次为浑善达克(11.49 kg/m2)、阴山(7.61 kg/m2)、呼伦贝尔(4.86 kg/m2)、塔里木(2.17 kg/m2),阿尔金的单位面积防风固沙量最低,为1.72 kg/m2。2010年防风固沙型重点生态功能区的防风固沙总量为5.55×1012kg,平均单位面积防风固沙量为6.72 kg/m2。从空间分布来看(图3),科尔沁南部、浑善达克中部和西北部、阴山中部的单位面积防风固沙量相对较高,塔里木西南部、阿尔金南部、呼伦贝尔的单位面积防风固沙量相对较低。刘璐璐等[6]基于RWEQ模型计算得到2000—2010、2010—2015年防风固沙型重点生态功能区年均防风固沙总量分别为12.52×1012kg、9.19×1012kg,年均单位面积防风固沙量为4.52 kg/m2、3.32 kg/m2。黄麟等[7]基于RWEQ模型计算得到2010年防风固沙型重点生态功能区平均单位面积防风固沙量为7.35 kg/m2,阴山、浑善达克、科尔沁的单位面积防风固沙量相对较高,分别为12.49 kg/m2、12.23 kg/m2、11.69 kg/m2,呼伦贝尔、阿尔金、塔里木的单位面积防风固沙量相对较低,分别为5.17 kg/m2、2.10 kg/m2、0.38 kg/m2。虽然防风固沙型重点生态功能区的划分范围因数据源和年限的不同有所差异,但是本研究的防风固沙量评估结果介于以上两者之间,具有可比性。

2.2 防风固沙型重点生态功能区的防风固沙服务空间流动路径分析

2010年各防风固沙型重点生态功能区的代表站点均有1460条每六小时的风速观测记录,浑善达克、呼伦贝尔、塔里木、阴山、科尔沁、阿尔金分别有194、183、141、109、73、55条风速记录超过当地起沙风速,分别对应194、183、141、109、73、55条防风固沙服务流动路径(图4)。从防风固沙服务流动路径的空间分布来看,塔里木和阿尔金的流动路径属于西北路径和偏西路径,西北路径部分随西北风向东南方向移动,后随气旋收缩北上转向东北方向移动,偏西路径则主要向偏东方向移动[47],两者主要流经中国境内地区,包括中国西北、华北、东北的小部分地区,在境外主要涉及俄罗斯远东地区、朝鲜半岛、日本等地。阴山、浑善达克、科尔沁、呼伦贝尔的防风固沙服务流动路径方向相似,均属于偏北路径,一般起源于蒙古国乌兰巴托以南的广大地区,主要影响西北地区东部、华北和东北的大部分区域[47],其跨境效应明显,境外主要流经朝鲜半岛、日本、俄罗斯远东地区,甚至会影响到美国阿拉斯加和菲律宾等东南亚地区。浑善达克、呼伦贝尔的路径分布密度相对较高。基于各个功能区防风固服务流动路径的空间叠加得到防风固沙型重点生态功能区的防风固沙服务流动路径(图4),共有755条。整体来看,路径集中分布在中国的西北、华北、东北的广大区域,境外主要涉及朝鲜半岛、日本、俄罗斯远东地区和北太平洋的广大海域,路径密度随着下风向传输距离的增加而降低。

2.3 防风固沙型重点生态功能区的防风固沙服务受益区识别

根据防风固沙服务流动路径的空间插值得到2010年各重点生态功能区的受益区范围及路径分布频率(图5)。从各个功能区来看,2010年呼伦贝尔的受益区面积最大,高达15.71×106km2,其次为浑善达克(15.66×106km2)、科尔沁(13.25×106km2)、阴山(11.50×106km2)、塔里木(8.64×106km2),阿尔金的受益区面积最小,为6.18×106km2。其中中国境内受益区占比最大的为阿尔金,占阿尔金受益区总面积的76.14%,其次为塔里木(57.05%)、阴山(33.10%)、浑善达克(20.74%)、科尔沁(15.68%)、呼伦贝尔(15.41%)。境内受益区总面积最大的为塔里木,占中国总面积的51.37%,其次为阿尔金(49.00%)、阴山(39.66%)、浑善达克(33.82%)、呼伦贝尔(25.21%)、科尔沁(21.65%)。从防风固沙服务流动路径的分布频率来看,各个功能区的防风固沙服务流动路径集中分布在路径分布频率低于2%的受益区内(表2)。防风固沙服务流动路径越密集,受益的防风固沙服务流动效应就越高。从防风固沙服务受益高频区域的空间分布来看(图5),塔里木的主要位于新疆西南部地区,阿尔金的主要位于新疆东南部、青海西北部地区,阴山的主要位于阴山东北部地区,包括内蒙古中部、京津冀地区西南部、山西北部,浑善达克的主要位于浑善达克东北部地区、辽宁与吉林西部地区,科尔沁的主要位于科尔沁沙地东南部与吉林、辽宁的边界地区,呼伦贝尔的集中在内蒙古东北部地区、吉林西部、黑龙江西南部。综合来看,各功能区防风固沙服务受益区范围内的受益频率均呈现明显的圈层状递减特征。将各个功能区的受益区范围进行叠加,得到防风固沙型重点生态功能区防风固沙服务的受益区范围(图5),总面积为32.16×106km2,86.26%的受益区范围内的受益频率低于2%。受益高频地区主要位于新疆西南部、内蒙古中部地区、京津冀、东北地区西部、山西北部地区,并以此为中心呈圈层状递减。

表2 2010年防风固沙型重点生态功能区防风固沙服务受益区路径分布频率统计/(106km2)

Table 2 The areas of the service beneficiary areas (SBAs) of the wind erosion prevention service (WEP) by the ecosystems in the national key ecological function zone of windbreak and sand fixation type (NKEFZ-WSF) in China in 2010 analyzed by frequency of trajectories

防风固沙服务流动路径分布频率Frequency of flow paths各防风固沙型重点生态功能区防风固沙服务受益区面积Areas of the SBAs of the WEP of each NKEFZ-WSF in ChinaTLMAEJYSHSDKKEQHLBR防风固沙型重点生态功能区受益区面积Areas of the SBAs of the WEP of the entire NKEFZ-WSF in China<2%5.80 3.16 7.71 10.41 7.23 10.97 27.74 2%—5%1.92 1.23 2.20 3.49 3.94 3.13 3.29 5%—10%0.50 1.09 0.95 1.26 1.75 1.07 0.99 10%—20%0.31 0.48 0.44 0.42 0.28 0.42 0.12 20%—30%0.06 0.11 0.13 0.05 0.03 0.07 0.03 ≥30%0.06 0.11 0.07 0.03 0.04 0.05 0.00 总面积Total area8.64 6.18 11.50 15.66 13.25 15.71 32.16

防风固沙型重点生态功能区代码说明：TLM：塔里木河荒漠化防治生态功能区;AEJ：阿尔金草原荒漠化防治生态功能区;YS：阴山北麓草原生态功能区;HSDK：浑善达克沙漠化防治生态功能区;KEQ：科尔沁草原生态功能区;HLBR：呼伦贝尔草原草甸生态功能区

图5 2010年防风固沙型重点生态功能区防风固沙服务受益区分布Fig.5 Service beneficiary areas of the wind erosion prevention service by the ecosystems in the national key ecological function zone of windbreak and sand fixation type in China in 2010

综合受益区范围内的防风固沙服务受益频率,海域及中国境外受益区的受益频率较低,防风固沙服务的流动效应不明显,因此本研究仅对中国境内陆地部分受益区的防风固沙服务流动效应进行分析。2010年防风固沙型重点生态功能区的境内受益面积为7.76×106km2,受益草地的占比最大,占境内受益区总面积的31.12%,占中国草地总面积的80.30%,主要位于防风固沙服务流动路径分布频率低于0.5%的受益区范围内(表3),基本涵盖了全国主要的草地分布区,包括内蒙古呼伦贝尔草原和锡林郭勒草原、新疆伊犁草原、西藏那曲高寒草原(图6)。其次为未利用地,占境内受益区总面积的23.24%,占中国未利用地总面积的93.96%,主要位于西北荒漠区(图6)。大部分受益荒漠的路径经过频率介于0.5%至1%(表3),防风固沙服务的存在能够降低沙尘流动路径经过荒漠区的频率,从而减少荒漠区沙尘的扩散。受益耕地和林地的面积相当,分别占境内受益区总面积的19.73%、19.59%,分别占中国耕地和林地总面积的85.40%、67.48%。受益耕地主要位于路径分布频率在2%—5%的受益区范围内,包括华北平原、东北平原、湖北的江汉平原、湖南的洞庭湖平原(图6),涉及全国重要的粮食主产区。受益林地主要位于路径分布频率低于0.5%的受益区范围内(表3),涉及我国东北地区、东南地区的主要林区(图6)。受益水体和湿地的面积占境内受益区总面积的3.97%,占全国水体和湿地总面积的85.96%,主要位于路径分布频率低于0.5%的受益区范围内(表3),涉及黄河、长江中下游水系、多数内陆湖泊等(图6)。受益的建设用地虽然仅占境内受益区总面积的2.36%,但是占中国建设用地总面积的92.21%,主要位于路径分布频率介于2%和5%的受益区内(表3),基本涵盖了全国主要的城市群(图6)。由于建设用地承载着人类主要的生活、生产活动,该部分受益区对保障人类福祉具有重要意义,虽然在受益区中占比较小,但是受益的效益最为显著。

图6 2010年防风固沙型重点生态功能区防风固沙服务受益土地覆被类型分布Fig.6 Spatial distribution pattern of the benefiting land cover of the wind erosion prevention service by the ecosystems in the national key ecological function zone of windbreak and sand fixation type in China in 2010

3 结论

本研究分析了转移支付前(2010年)防风固沙型重点生态功能区防风固沙量的空间格局,从沙尘传输路径的角度识别出各个防风固沙型重点生态功能区的防风固沙服务流动路径与受益区范围,进而得到整个防风固沙型重点生态功能区的防风固沙服务流动路径、受益区范围、受益的土地覆被类型,防风固沙服务功能效应通过流动路径在受益区内得以实现。主要结论如下：(1)2010年防风固沙型重点生态功能区的防风固沙总量为5.55×1012kg,其中浑善达克的防风固沙总量最高,占34.36%。平均单位面积防风固沙量为6.72 kg/m2,其中科尔沁的单位面积防风固沙量最高,为12.51 kg/m2。(2)2010年防风固沙型重点生态功能区共有755条防风固沙服务流动路径,主要流经中国的西北、华北、东北的广大区域,朝鲜半岛,日本,俄罗斯远东地区和北太平洋的广大海域,路径密度随着下风向传输距离的增加而降低。(3)2010年防风固沙型重点生态功能区防风固沙服务的受益区面积为32.16×106km2,大部分受益区的受益频率均低于2%,受益高频地区主要位于新疆西南部、内蒙古中部地区、京津冀、东北地区西部、山西北部地区,并以此为中心呈圈层状递减。(4)2010年防风固沙型重点生态功能区防风固沙服务的境内受益区面积为7.76×106km2,其中草地占比最大,其次为未利用地、耕地、林地,水体和湿地、建设用地的占比较小,但是受益建设用地占中国建设用地总面积的比例最高,高达93.96%,基本涵盖了中国主要的城市群区域,防风固沙服务流动效益明显。本研究能够为防风固沙型重点生态功能区的区域间生态补偿标准核算提供科学依据,对健全防风固沙型重点生态功能区生态补偿机制、弥补转移支付政策的不足具有重要的现实意义。

4 讨论

本研究明确了防风固沙型重点生态功能区防风固沙服务的受益区范围,主要位于东部经济发展水平较高的地区,防风固沙服务流动为保障这些地区的经济发展具有贡献作用。防风固沙服务的维持与提升涉及防风固沙型重点生态功能区的多项生态保护政策,但是主要集中在国家政府的纵向生态补偿,区域间生态补偿机制缺乏。因此本研究中受益区范围的确定有助于明确补偿主体,制定区域间生态补偿政策,弥补功能区中农牧民为保护植被、增强防风固沙服务而牺牲的发展机会成本,完善我国的生态补偿制度。但是,本研究仍然存在一些不确定性和限制因素。(1)防风固沙服务量的计算中：没有考虑不同土壤类型临界起沙风速的差异；由于数据限制,风力因子计算使用的风速数据为月均尺度,“平滑化效应”导致防风固沙量的估算结果偏小[48]；地表糙度因子计算过程中忽略了地势起伏变化的尺度效应。(2)起沙机制的简化：假定只要在风速大于等于临界起沙风速的条件下,有充足的沙源可供流动。然而沙尘的传输需要同时满足以下三个条件：充足的沙源、强风和大气的不稳定。沙尘的扩散机制复杂,不同粒径的沙粒在传输过程中沉降的先后顺序与多种因素有关,且沙尘的扩散与季节关系密切,不同季节的植被覆盖度、土壤湿度等影响起沙的条件差异显著,本研究均未作考虑。(3)气象数据的空间异质性：由于数据的限制,采用距离各个防风固沙型重点生态功能区中心点最近的气象站风速数据模拟各防风固沙型重点生态功能区的气流运动,没有考虑气象数据的空间异质性,从而影响模拟结果的准确度。今后的研究中应结合实地观测数据,提高防风固沙量的时空模拟精度,将沙尘流动过程的沉降、传输量计算与防风固沙量相结合,明确防风固沙服务流动的物质流过程,进一步建立防风固沙服务供给区与受益区之间服务流动的时空定量关系。同时,结合防风固沙服务的价值核算,确定防风固沙服务流动的价值流过程,建立防风固沙服务价值流与区域间生态补偿之间的时空定量关系,为防风固沙服务的区域间生态补偿提供更为直接的科学依据。