三北地区是我国最重要的自然资源储备区和生态屏障区，也是我国荒漠化最严重的地区^[1]，荒漠化土地面积占全国荒漠化土地面积的55%以上^[2]。地表径流和地下水资源不足，地面植被覆盖度不高、生产量偏低，生态环境极其脆弱^[3]，风沙危害和水土流失十分严重^{[4, 5]}。为了从根本上改善三北地区的生态环境状况，我国于1978年启动了三北防护林体系建设工程，至2010年完成了三北防护林四期工程。2000年以来，我国在三北地区还实施了退耕还林(还草)、退牧还草、天然林保护等重大生态保护与建设工程。

许多学者开展了三北地区生态环境状况相关的研究，集中在三北地区防护林体系区域分异与开发潜力^[3]、区域气候变化及防护林的气候生态适应性^{[5, 6, 7]}，区域蒸散发变化^[8]，乔、灌木树种最适分布范围^[9]，防护林空间分布信息提取及遥感监测^{[10, 11]}，多年NDVI与植被覆盖度时空变化^{[12, 13, 14]}，土地退化及趋势^[15]，水土流失和土壤侵蚀估算^{[16, 17]}等。此外，探究三北地区生态问题亦成为众多学者的一个主要研究方向，比如违背植被地带性分布规律，重乔轻灌草、土壤水分亏缺、土壤肥力不足和不平衡等^{[18, 19, 20]}。然而，上述研究主要针对三北地区生态环境单因素开展分析，缺少对该地区生态系统整体状况的全面、科学、准确地把握。

人们迫切需要了解，三北地区过去10年实施了大量重大生态建设工程后，区域生态系统如何变化？发生变化的主要驱动因素是什么？哪些因素产生正面作用，哪些产生负面作用？因此，本文试图通过分析2000—2010年生态系统类型、质量以及防风固沙、土壤保持生态系统服务的时空格局变化，了解气候变化与生态工程背景下三北地区近10年生态系统变化状况，并分析气候因子与人类扰动对三北地区生态系统变化的驱动作用，为科学地评估生态工程的生态效应以及后续生态工程的滚动实施提供理论依据。

1 研究方法与数据 1.1 研究区概况

中国三北地区系指昆仑-秦岭-大巴山以北的部分地区，横跨我国西北、华北和东北三大地域^[2]，辖13省(市)，面积约占国土面积的42.4%。该区气候差异显著，大部分属于干旱、半干旱气候带，降水量由东向西、从南到北递减，年均降水量在20—450mm^[4]。受降水量的影响，天然植被类型从东到西依次为森林、森林草原、草原和荒漠。

1.2 生态系统类型及其时空变化分析

三北地区生态系统类型主要划分为林地、草地、湿地、耕地、人工表面、其它等6个类型。林地即木本为主的植物群落，郁闭度不低于0.2，高度在0.3m以上，包括自然、半自然植被及集约化经营和管理的人工木本植被。草地是一年或多年生草本植被为主的植物群落，覆盖度大于20%，高度在3m以下。湿地包括沼泽、水域、永久性冰川、滩地等。耕地是以收获为目的、有耕犁活动的人工种植草本植物覆盖表面。人工表面即人工建造用于城乡居民点、工矿、交通等的陆地表面。其它包括年内最大植被覆盖度小于20%的地表及冰雪。基于2000、2005、2010年TM/ETM+、环境小卫星等30m空间分辨率遥感影像数据，经图像精纠正和拉伸处理后，通过遥感解译判读，生成几期生态系统类型空间分布数据。通过1°×1°经纬度交叉点的野外调查信息进行结果的校正与验证，类型精度达到85%以上。进而对生态系统类型空间数据进行统计分析，综合评价6大生态系统类型的时空变化态势。

1.3 生态系统质量及其变化分析

收集2000年至2010年MODIS的16天250m分辨率NDVI数据(MOD13Q1)，数据进行格式转换、重投影、图像的空间拼接、重采样和滤波处理，采用最大合成法(MVC)得到年NDVI数据。利用NDVI数据计算植被覆盖度，根据像元二分模型理论，可以认为一个像元的NDVI值是由绿色植被部分贡献的信息与无植被覆盖部分贡献的信息组合而成，最大植被覆盖度根据以下公式计算：

式中，F_C是植被覆盖度，NDVI_veg是纯植被像元的NDVI值；NDVI_soil是完全无植被覆盖像元的NDVI值，纯植被像元和完全无植被像元根据生态系统类型数据进行识别。

采用CASA模型^[21]计算植被净初级生产力，输入参数包括气象数据插值的1km平均温度、蒸散量、日照时数，MODIS1km分辨率反照率与植被指数产品，数字化的1 ∶ 100万植被类型分布等，计算过程如下：

从资料文档、气象数据中得到太阳总辐射量及日照时数等信息，然后结合像元经纬度计算得到PAR。利用NDVI计算比值指数SR，通过FPAR与SR之间存在关系得到FPAR。

式中，FPAR_min和FPAR_max的取值与植被类型无关，分别取值为0.001和0.95，SR_min和SR_max与植被类型有关，为对应植被类型NDVI的5%和95%的下侧百分位数。NIR和RED分别表示近红外波段和红波段的反射率。

植被将吸收的光合有效辐射(APAR)通过光合作用转化为有机碳的效率计算如下：

式中，ε^*是最大光利用率(g/MJ),T₁和T₂表示环境温度对光利用的抑制影响，W为水分影响胁迫系数，均为无量纲参数。

叶面积指数(LAI)是利用遥感反射率数据，结合生态系统分类数据得到调查区域内不同植被类型的像元的NDVI，然后采用查找表的方法，反演区域尺度上植被冠层LAI。

1.4 生态系统服务及其变化分析 1.4.1 防风固沙服务量变化分析

在充分考虑气候条件、植被状况、地表土壤的粗糙度、土壤可蚀性、土壤结皮的情况下，利用修正的土壤风蚀方程(RWEQ)定量评估土壤风蚀量。

式中，SL表示土壤风蚀模数，x表示地块长度，Q_x表示地块长度x处的沙通量(kg/m)；Q_max表示风力的最大输沙能力(kg/m)；s表示关键地块长度(m)；WF表示气象因子；EF表示土壤可蚀性成分；SCF表示土壤结皮因子；K′表示土壤糙度因子；COG表示植被因子，包括平铺、直立作物残留物和植被冠层。气候因子中的风和土壤湿度根据中国气象科学数据共享服务网(http://cdc.cma.gov.cn)下载的日均风速、降水、温度、日照时数等计算；雪盖因子利用从中国西部环境与生态科学数据中心(http://westdc.westgis.ac.cn)下载的中国雪深长时间序列数据集来计算。土壤可蚀性因子根据Fryear等^[22]方程计算。土壤数据来源于西部环境与生态科学数据中心1 ∶ 100万土壤类型图所附的土壤属性表和空间数据。根据文献收集地面测定的我国不同地区的风蚀模数结果^{[23, 24, 25, 26]}对本估算结果进行验证。

通过对生态系统防风固沙服务量进行定量分析来衡量生态系统防风固沙的能力。防风固沙服务量为裸土条件和地表覆盖植被条件下的土壤风蚀量的差值：

式中,SL_SV表示防风固沙服务量，SL_S表示裸土条件下的潜在土壤风蚀量，SL_V表示植被覆盖条件下的现实土壤风蚀量。

1.4.2 土壤保持服务量变化分析

基于修正的通用水土流失方程(RUSLE)计算单位面积土壤流失量，即土壤侵蚀模数。RUSLE包含降雨侵蚀力因子(R)、土壤可蚀性因子(K)、坡长因子(L)、坡度因子(S)、覆盖和管理因子(C)等:

基于国家气象台站的日降雨量数据以及日平均气温等气象资料，采用日降雨量拟合模型^[27]来估算降雨侵蚀力。利用1 ∶ 100万土壤数据库，根据Nomo图法估算土壤可蚀性因子值。基于SRTM3 V4.1数据加工制作得到的90m分辨率DEM数据，根据McCool等^[28]和Liu等^[29]的核心算法计算坡度坡长因子。盖度和管理因子采用蔡崇法^[30]提出的利用植被覆盖度计算C值的方法。通过查找文献搜集到多个地区的土壤侵蚀模数资料，将估算得到的土壤侵蚀模数与其进行对比验证。

通过对生态系统土壤保持量进行定量分析来衡量生态系统保育土壤的能力。土壤保持量为生态系统在极度退化状况下的土壤流失量与现实状况下土壤流失量的差值。

式中，A_C为土壤保持量，A_D和A_R分别表示生态系统在极度退化状况下的土壤流失量与现实状况下土壤流失量，可以发现A_D和A_R的差异仅仅体现在C因子。

1.5 驱动力因子分析

采用最小二乘法分析气温、降水等气象要素的年际变化趋势：

式中，S为变化斜率，X_i为气温、降水量观测值，i＝1,2,3,…,n,m_i为年份序数，m₁=1,m₂=2,m₃=3,…,m_n=n。

此外，根据赵国松等^[31]，通过人类扰动指数评价人类活动对自然生态系统的影响程度，扰动指数值越高，人类扰动程度越高，当同一区域人类扰动程度下降时，认为人类对自然生态系统的威胁程度减小。只考虑有自然植被分布的生态系统类型的人类扰动，而对稀少植被或无植被分布生态系统类型的人类扰动不作考虑。将生态系统分级赋值，得到4种扰动程度的分级指数，如表 1所示。

根据不同类型所占比例进行加权求和，形成一个0—3之间的值，标准化后扰动指数值域范围0—1，综合反映了某一地区的扰动程度，计算方法如下：

式中，A_i表示第i级生态系统扰动程度分级指数，P_i表示第i级扰动程度分级面积所占百分比，D为人类扰动指数。为便于多期综合人类扰动指数的对比，直接扣除无植被或稀疏植被分布生态系统类型占比超过95%的区域。

2 结果与分析 2.1 三北地区生态系统宏观结构及其变化

结果表明(图 1，表 2)，2010年，三北地区以荒漠等其它类型与草地生态系统为主，其它类型面积约占全区总面积的44%，草地占全区总面积的22.9%。其次是耕地和林地，分别占全区的15.8%和13%。此外，湿地面积占2.6%，人工表面占1.6%。2000—2010年，三北地区林地面积净增加了0.69万km²(表 3)，转入林地面积11312.4 km²，以草地和耕地转入林地为主，林地转出面积4438.6 km²，以转出到耕地为主。草地面积总体减少了0.13万km²，其中，转入草地面积6848.8 km²，以耕地和其它转入草地为主，草地转出面积8132.5 km²，主要转出为耕地。湿地面积净减少0.21万km²，其中转入湿地面积2360km²，湿地转出面积4415.7 km²，主要转出为耕地和其它。耕地面积减少了0.87万km²，其中转入耕地面积7209.2km²，以林地、草地转入耕地为主，耕地转出面积15906.1km²，主要转出为人工表面。人工表面增加了0.57万km²，其它类型面积持续减少了496.5km²。

2.2 三北地区生态系统质量状况及其变化

2000—2010年，三北地区NDVI平均值约0.39，总体呈现轻微弱增加趋势(图 2，表 4)，年增加速率0.008，特别是黄土高原丘陵沟壑区，年增速为0.077，仅内蒙中东部呈现显著减少趋势。三北地区多年平均最大植被覆盖度为37.66%，总体呈现微弱减少趋势，年减少速率约为0.07%，其中东北华北区平均最大植被覆盖度89.94%，近10年减少速率最大，约为0.196%/a；黄土高原丘陵沟壑区平均最大植被覆盖度61.02%，近10年表现为大面积覆盖度增加趋势，增加速率0.611%/a。多年平均植被LAI约为1.13，近10年总体呈现微弱增加趋势，增加速率0.009/a，特别是黄土高原丘陵沟壑区和东北华北平原农区，年增加速率分别为0.041和0.024。近10年，三北地区多年平均地上植被生产力为236.65 gC m^-2 a^-1，总体呈现增加趋势，年增速为1.345gC/m²，其中，黄土高原区年增速高达8.828gC/m²。

2.3 三北地区生态系统关键服务及其变化 2.3.1 生态系统防风固沙服务及其变化

2000—2010年，三北地区风蚀区多年平均单位面积土壤风蚀模数为20.69 t/hm²(表 5，图 3)，其中西北荒漠区最高，约30.71 t/hm²，黄土高原丘陵沟壑区和东北华北平原农区最低，分别约为1.07和2.01 t/hm²。近10年，全区单位面积土壤风蚀模数呈现减少趋势，减少速率为每年1.13 t/hm²，年平均风蚀模数下降了25.25%，其中减速最高的是西北荒漠区。甘肃西部与蒙新衔接区域、柴达木盆地等则呈现明显增加趋势。三北地区多年平均生态系统防风固沙服务量为32.17 t/hm²，其中风沙区最高，约53.31t/hm²，黄土高原丘陵沟壑区最低，约为6.06 t/hm²。然而，近10年全区年平均防风固沙服务量则呈现减弱趋势，这与风蚀力逐年减弱有关，风蚀驱动力的减少使得裸土和植被覆盖条件下的土壤风蚀量明显减少，服务量也相对减少。但这并不能表明生态系统的防风固沙服务能力在降低，而仅仅表示在风蚀力减小的情况下，土壤防治的风蚀量也相应减少。

2.3.2 生态系统土壤保持服务及其变化

2 000—2010年，三北地区水蚀区多年平均土壤水蚀模数为10.37 t/hm²(表 5，图 4)，其中黄土高原丘陵沟壑区最高，约31.23 t/hm²，东北华北平原农区和风沙区最低，分别约为2.73和5.79 t/hm²。近10年，全区土壤水蚀模数总体呈现微弱增加趋势，年增加速率为0.08 t/hm²。而黄土高原区和风沙区的土壤水蚀模数则呈现降低趋势，特别是黄土高原区，年减少1.93 t/hm²。三北地区水蚀区多年平均土壤保持服务量为12.31 t/hm²，其中黄土高原丘陵沟壑区最高，约54.76 t/hm²，东北华北平原农区和西北荒漠区最低，约为8.7和9.3 t/hm²。近10年，单位面积生态系统土壤保持服务量呈现增加趋势，年增速0.28 t/hm²，其中黄土高原区年增速达到1.23 t/hm²。

2.4 三北地区生态系统变化的驱动力分析

2000—2010年，三北地区气温以增温为主，年平均气温变化斜率为0.02℃/a，特别是西北荒漠区。年降水量以增加趋势为主，变化斜率为3.00mm/a，特别是东北华北平原农区(图 5，表 6)。适度增温、降水量增加有助于增加植被高度和生产力，加快植物物候进程，延长生长季，促进植物生长发育。近10年，三北地区人类扰动指数增加了0.003(图 6，表 6)，特别是风沙区和西北荒漠区，分别增加了0.309和0.007，而东北华北平原农区和黄土高原丘陵沟壑区则呈现相反的减少趋势，分别减少了0.151和0.005，说明人类活动对这两个区域自然生态系统的扰动有所遏制，有利于植被恢复。

3 讨论与结论

本文分析了2000—2010年我国三北地区生态系统宏观结构、质量及关键服务变化状况及其驱动因素，得到几点结论：

首先，三北地区林地面积持续增加，草地面积持续减少。根据林业统计年鉴，2000—2010年，三北防护林体系工程完成人工造林面积达4.61万km²。而本文结果显示10年间，林地面积净增加了0.69万km²，由其他类型转变为林地的面积为1.13万km²，仅占统计面积的24.51%，主要原因是三北防护林中的农田防护林、公路大绿化和村旁植树难以遥感识别。

1 0年间，三北地区生态系统质量和服务能力有所提升，NDVI、LAI、植被生产力总体呈现轻微增加趋势，特别是黄土高原丘陵沟壑区，然而植被覆盖度总体呈现微弱减少趋势，部分区域的植被覆盖度呈现增长趋势，特别是陕甘宁地区，而新疆北部、内蒙古东部等植被覆盖度呈下降趋势。近10年，三北地区由于风蚀力逐年减弱导致土壤风蚀模数降低，同时生态系统防风固沙服务量也呈现减弱趋势。水蚀区土壤水蚀模数总体呈现微弱增加趋势，而生态系统土壤保持服务量呈现增加趋势。

从过去10年三北地区生态系统时空变化可以看出，近40年的生态工程在三北地区生态建设中发挥了极为重要的作用，但也存在较多问题，本区多为干旱半干旱区，很多地方抽取地下水喷灌、滴灌造林，甚至移土造林、在山地草甸上造林，不仅成本高、成活率也低，破坏了草地原生植被，而且不可持续。因此，未来三北地区生态建设应该基于“宜乔则乔、宜灌则灌、宜草则草、宜荒则荒”原则，优先保护优良、原生生态系统，以自然恢复为主。此外，如何科学正确地评估三北地区生态工程的生态效应，是调整后续工程措施的主要依据，也是未来的研究重点。