竹林是森林生态系统的重要植被类型，我国竹林面积占森林面积的3.06%，竹林碳储量占森林碳储量的比例超过11%(陈先刚等, 2008)。针对竹林生态系统碳储量及碳汇功能，我国学者在全国尺度或典型区域上做了大量研究(王兵等, 2008; 黄从德等, 2007; 周国模等, 2004; 杨春花等, 2008)，但由于竹林研究基础数据的不完整(邓旺华等, 2008)和计算方法的差异，不同学者在全国尺度上估算的竹林生态系统碳储量差异较大(李正才等, 2003; 郭起荣等, 2005; 陈先刚等, 2008; 王兵等, 2008)，典型区域竹林碳储量也因竹种、人工抚育以及研究方法的不同差异显著(冯帅等, 2010)。因此开展中尺度竹林碳储量研究十分重要。

中国自20世纪90年代后期实施退耕还林工程，显著改变了森林面积、土地利用/覆盖类型(张志明等, 2009)、生态系统碳储量(Chen et al., 2009)和碳储量空间分配格局(胡亚林等, 2009)。研究发现，退耕竹林碳储量高于白桦(Betula platyphylla)林和落叶松(Larix olgensis)人工林等乔木林种(吴小山等, 2007; 李江等, 2006; 白雪爽等, 2008; 王春梅等, 2007)，其碳汇潜力巨大(蔡丽莎等, 2009; 陈先刚等, 2009)。随着遥感技术的发展及其所具有的局限性，以往对竹林动态变化的研究大多与土地利用类型变化有关(刘华等, 2008; 施拥军等, 2008), 退耕还林的实施对竹林碳储量影响的研究较少。本研究尝试利用3S技术，在极具退耕还林代表性的洪雅县研究中尺度退耕竹林碳储量分布格局。本研究选用1994—2007年3期Landsat卫星影像，利用遥感数据结合地面调查和竹林碳储量实测数据，构建竹林碳储量遥感模型估算各期竹林碳储量，在此基础上研究竹林面积、碳储量及碳密度的时空分布格局，并探讨退耕还林工程实施对竹林碳储量的影响，以期为正确评估退耕还林生态作用提供数据支撑，为退耕还林工程宏观决策提供依据。

1 研究区概况

研究区位于四川盆地西南边缘华西雨屏区核心区洪雅县(102°49′—103°32′ E，29°24′—30°00′ N)，总面积1 948.43 km²，耕地总面积19.52万hm²。海拔417~3 090 m，地形由西南向东北高低梯次变化形成高山、中山、深丘、浅丘、台地、河谷和平坝，地貌以山地丘陵为主。气候属亚热带湿润性山地气候，年均气温14~16 ℃，年降雨量1 300~2 000 mm，年均日照1 080 h，全年无霜期352天。全县地带性土壤为黄壤，兼有冲积土、紫色土，境内土壤侵蚀较为严重。地带性植被类型为亚热带常绿阔叶林。退耕竹林栽植于2000年，主要竹种为慈竹(Neosinocalamus affinis)、麻竹(Dendrocalamus latiflorus)、撑绿杂交竹(Bambusa pervariabilis×Dendrocalamopsis grandis)和苦竹(Pleioblastus amarus)(李贤伟等, 2009)。

2 研究方法

数据源于覆盖整个洪雅县区域的退耕还林以前(1994-06-26)、退耕还林后第4年(2004-04-18)和退耕还林后第7年(2007-09-18)的3期TM影像(30 m分辨率)、洪雅县地形图(比例尺1:50 000)、洪雅县林相图(比例尺1:10 000)、1999—2005年洪雅县各乡镇退耕还林工程验收图(比例尺1:10 000)、洪雅县各乡镇1999—2005年退耕还林年度总结报告及洪雅县2007年森林资源二类调查小班档案。

利用ENVI4.7对3期原始影像进行辐射定标、FLAASH大气校正及几何校正，结合实地调查和退耕还林小班验收图，对研究区3期TM遥感影像进行分类(总体分类精度分别为87.35%，86.56%和89.32%)并提取出洪雅县退耕竹林分布区域。在洪雅县范围内共布设474块30 m×30 m的样地，落入退耕竹林地的样地92块，其中可测样地87块(随机选取57块样地数据用于建模，剩余30块用于模型预测精度检验)。样地实测面积为30 m×30 m，与TM影像的像元大小一致，对样地每竹调查其胸径、地径、竹高、枝下高和秆龄，用GPS对样地中心点定位。

根据87块样地调查资料，按竹高、胸径和秆龄每样地选取标准竹2~3株进行全竹取样，分别取竹叶、竹枝、竹秆、竹蔸、竹根和竹鞭，野外称其鲜质量，并各取1 kg鲜样室内测定其生物量和有机碳含量。有机碳含量的测定采用重铬酸钾-外加热法(鲁如坤，1999)。碳储量为生物量乘以碳转换系数即有机碳含量，表达式为C_i=C_c×W_i，C_i为碳储量(t), C_c为碳含量, W_i为生物量(t)。

碳储量遥感模型构建及优化：利用ArcGIS 9.3对地形图进行配准并矢量化，生成分辨率为30 m的数字高程模型(digital elevation model，简称DEM)。在ENVI4.7中将TM影像、DEM与87块竹林样地叠合提取影像各波段光谱反射率ρ因子、植被指数因子、主成分因子、穗帽变化因子以及地形因子作为遥感建模初选反演参数。对所选参数进行相关性和共线性分析(表 1)，选取相关性高、共线性小的因子构建碳储量遥感线性回归方程。综合考虑方程的相关系数、复相关系数、调整相关系数、显著水平、模型独立性检验指标值及方差膨胀因子值，初步遴选出模型C=73.042+0.011E－1.954ρ₂－0.633ρ₇+0.369P₃－15.826R+51.555M(表 2)，用地统计学软件对样点残差克里金插值得到最优模型C=82.869+0.011E+0.013A－1.602ρ₂－0.925ρ₇+0.837P₃+2.942M。C为碳储量(t), E为海拔(m), A为坡向, ρ₂为第2波段光谱反射率, ρ₇为第7波段光谱反射率, P₃为第3主成分, R为比值植被指数, M为土壤调整植被指数。

用未参加建模的30块样地实测碳储量与最优模型预测的该30块样地的碳储量进行拟合, 结果表明(图 1)，插值后实测碳储量与预测值的拟合度有一定程度的提高，预测精度为83.01%，较初选模型(81.56%)提高1.45%。本研究利用优化模型估算洪雅县1994，2004和2007年碳储量。

3 结果与分析 3.1 竹林面积动态

由图 2, 3和4可见，1994—2007年洪雅县竹林分布区域集中在中部和北部，1994年各竹林斑块面积较小，分布较零散，2004和2007年各竹林斑块呈现出片状分布。13年间竹林面积增加了0.77万hm²，其中1994—2004年增加了0.55万hm²，2004—2007年增加了0.21万hm²，年均增长率为4.39%。竹林面积的大幅增加主要与该时期内退耕还林工程的广泛开展有关，退耕还林工程实施后，全县竹林面积占县域面积的比重由1994年的5.19%提升至2007年的9.15%，约提升了4个百分点。

为获取竹林面积随各种地形因素的分异规律，在生成海拔、坡度、坡向专题数据基础上，为方便统计与分析，结合研究区各地形因子实际情况，采用ArcGIS 9.3对各因子进行了分级(李志林等, 2001)，分级标准见表 3。

将3期影像分类图与海拔、坡向、坡度分级专题图叠加分析发现，竹林空间分布向低海拔、平缓坡发展，受坡向影响不大。具体表现为：1994年竹林主要分布在海拔1 500 m以下，其中海拔1 000~1 500 m和500~1 000 m区域竹林面积分别占洪雅县竹林总面积的33.70%和64.12%。退耕还林工程实施后，竹林面积在＜1 000 m海拔区域增幅较大，2004年面积增加了0.60万hm²，为1994年的1.90倍，2007年面积增加0.23万hm²，为2004年的1.18倍, 在＞1 000 m海拔区域竹林面积逐年降低，2007年较1994年减少了17.24%，较2004年减少了4.96%。竹林随坡向的分布在各时期呈现相似规律，表现为阴坡＞半阴坡＞半阳坡＞阳坡＞无坡向。无坡向比重最少，约为4%，在其余4个坡向中的分布比重相差不大，约为阴坡29%、半阴坡25%、半阳坡22%、阳坡20%(表 4)。

退耕还林以前竹林广泛分布于平缓坡(5°~25°)区域，其中10°~25°分布最多，达到0.57万hm²，占总面积的56.28%。退耕还林工程使得坡度＞25°的农耕地大量退耕，1994—2004年，竹林面积增加了0.55万hm²，其中坡度＞25°区域竹林面积在2004年增加了0.43万hm²，为退耕还林前的2.87倍，其占总面积的比例由1994年的22.74%提升到42.25%，为退耕前的1.86倍。可见退耕还林的实施使陡坡竹林面积大幅增加。2007年竹林面积在＞35°坡度段略有减少，其他坡度段均呈现增加趋势。总体表现为＜25°区域增加量大于＞25°区域，竹林分布呈现出逐渐向平缓坡发展趋势(图 5)。

3.2 竹林碳储量动态

13年间竹林碳储量呈现持续增长的特征，共增加了25万t，年均增长4.33%，前期增速略大于后期。1994年洪雅县竹林碳储量为34万t，2004年竹林碳储量较1994年增加了18万t。这主要得益于退耕还林工程的实施，退耕还林开展后洪雅县大量发展竹林，竹林面积的大幅增加促进了竹林碳储量的增加。2007年竹林碳储量为59万t，较2004年增加了7万t。这些数据说明13年间洪雅县竹林均表现为碳汇，碳汇量为2万t·a^-1。竹林碳密度在1994—2004年和2004—2007年这2个时段呈现先降后升的趋势。1994年竹林碳密度为33.76 t C·hm^-2。退耕还林初期，由于大量的退耕还竹，未成熟竹林比重增加，2004年竹林碳密度有所下降，为33.25 t C·hm^-2。到2007年，竹林碳密度略微增加，达到33.67 t C·hm^-2，尚未恢复到1994年的水平，低于全国森林平均水平41.10 t C·hm^-2(Fang et al., 2007)，说明洪雅县竹林碳吸存能力还比较低。究其原因，研究期间退耕还林新栽竹林多处于中、幼龄阶段，人工抚育措施和管理水平不高，从而导致碳密度水平较低。随着竹林的生长、成熟和经营管理水平的提高，洪雅县竹林碳密度、碳储量尚有增加潜力。

对比分析竹林碳储量随海拔的分布发现(图 6)，1994—2007年竹林碳储量主要分布于海拔2 000 m以下，其中高海拔(＞1 000 m)区域竹林碳储量逐年减少，中低海拔(＜1 000 m)区域竹林碳储量逐渐增加，且增加趋势明显。2004年高海拔区域竹林碳储量较1994年(11.65万t)减少1.66万t，2007年较2004年减少2.97万t；在500~1 000 m海拔区域，2004年竹林碳储量较1994年(21.90万t)净增加16.70万t，增加了1.76倍，2007年较2004年增长量有所下降，增加了6.09万t；在＜500 m海拔区域竹林碳储量持续增加，且增速加快，2004年较1994年(0.61万t)增加2.84万t，年均增长19.03%，2007年较2004年增加4.38万t，年均增长31.44%。

对比分析各时期竹林碳储量随坡度的分布发现(图 7)，退耕还林以前洪雅县竹林碳储量主要分布于＜25°坡度区域，为26.39万t，占总碳储量的77.26%。其中以15°~25°坡度区域分布最多，为12.48万t，占总碳储量的36.54%，＞35°坡度区域分布最少，仅占6.54%。2004年竹林碳储量在各个坡度段均有一定增加(10°~15°坡度区域除外)，其中＞25°坡度区域的增加最为明显，增加了14.22万t，为退耕还林前的2.83倍，这主要得益于退耕还林工程的实施。而在＜25°坡度区域竹林碳储量也有一定增长，增加了3.66万t，为1994年的1.14倍，这部分竹林碳储量的增加主要是瓦屋山旅游沿途竹林景观带的建设，使得平缓坡竹林面积增加所致。随竹林年龄结构变化和竹林分布改变，竹林碳储量在2007年表现出在＜15°坡度区域持续增加，在25°~35°坡度区域增速减慢，在15°~25°和＞35°坡度区域略有减少。说明随着时间推移，竹林碳储量分布区逐渐向平缓坡发展，在陡坡区域竹林碳储量呈逐渐减少趋势。

对比分析竹林碳储量随坡向的分布发现(图 8)，各时相竹林碳储量随坡向的分布规律与竹林面积随坡向的分布规律一致，表现为阴坡(比重为29.10%)＞半阴坡(比重为25.17%)＞半阳坡(比重为21.93%)＞阳坡(比重为20.27%)＞无坡向(比重为3.56%)。随着退耕还林工程的实施和推进，阴坡、半阴坡、半阳坡、无坡向竹林碳储量都呈现持续增加的趋势，阳坡碳储量表现为先增加后降低。2004年较1994年竹林碳储量大幅增加，以半阴坡、半阳坡和阳坡为主，分别增加了4.81万，5.55万和4.40万t。无坡向竹林碳储量增加最少，仅为0.64万t。2007年较2004年半阴坡增加量最多，为2.26万t，其次为无坡向，增加量为2.17万t，阳坡减少0.16万t。对比各年份不同坡向竹林碳储量所占比重发现，1994—2007年期间半阴坡、无坡向竹林碳储量呈现增加趋势，半阳坡、阳坡呈现先增后降趋势，阴坡则逐年降低，各坡向比重差异逐渐缩小，表明竹林碳储量呈现各坡向均匀发展趋势。

4 结论与讨论

13年间洪雅县退耕还林后竹林面积和碳储量均呈现逐年上升趋势。具体表现为：1994—2007年竹林碳储量呈现逐年上升的趋势；退耕还林后高海拔区域竹林碳储量逐年减少，中低海拔区域竹林碳储量逐渐增加，且增加趋势明显；退耕还林的实施使得＞25°坡度区域竹林碳储量大幅增加，增加了14.22万t，随着退耕还林的深入，竹林碳储量呈现逐渐向平缓坡发展趋势；随坡向的分布差异不大，总体表现为阴坡＞半阴坡＞半阳坡＞阳坡＞无坡向；退耕还林工程实施后，洪雅县竹林斑块连通性增大，面积增幅较大，1994—2007年竹林面积增加了0.76万hm²；竹林分布呈现向低海拔区域发展的趋势，在500~1 000 m海拔区域分布最多，分别占各年洪雅县竹林总面积的64.12%，74.17%和75.82%；随着退耕还林工程的实施与推进，竹林面积在＞25°坡度区域增加速度减慢或面积减小；竹林面积随坡向的分布表现为：无坡向分布最少，其余4个坡向分布比重相差不大，半阴坡、无坡向呈现增加趋势，半阳坡和阳坡呈现先增后降趋势，阴坡则逐年降低。

洪雅县竹林碳密度表现出1994—2007年先降后升，总体略微下降的特点，以1994年最高为33.76 t C·hm^-2，2004年下降到33.25 t C·hm^-2，2007年又恢复到33.67 t C·hm^-2，整体上略高于毛竹(Phyllostachys edulis)林碳密度(29.46~32.99 t C·hm^-2)(周国模等, 2004; 王兵等, 2009)，低于四川省(38.04 t C·hm^-2)和全国(41.10 t C·hm^-2)森林平均碳密度(黄从德等，2009; Fang et al., 2007)。究其原因主要是植被类型不同和人类干扰所致，退耕还林工程大量造林更新使得洪雅县竹林年龄结构发生了改变，幼龄林、中龄林和近成熟林比重加大，成熟林、过熟林比重下降，这种林龄结构导致竹林碳密度出现了一定程度的降低。随着时间推移，竹林整体成熟度不断增加，同时竹林经营管理水平提高，竹林碳密度将不断增大，并且这种回升势头将在一定时间内得到保持，然后稳定在一个较高的水平，碳储量尚有较大增长空间。

本研究基于模型实用性考虑，利用地统计学建立碳储量遥感回归模型，对洪雅县退耕竹林进行了碳储量估算，模型变量少，简单实用，但模型的估测精度不高仅为83.01%，低于非线性模型(Gilabert et al., 1996; 张连义等, 2008)以及非参数模型(王淑君等, 2007; 张元元, 2009; 刘恩斌等, 2010)；对于所建立的竹林碳储量遥感模型，在进行后续研究时应考虑更多的非生物因子(王敏等, 2010; 徐小军等, 2008)，如区域气候、竹龄和立地指数等。同时，由于现有影像数据的限制，不能有效区分竹种，因此本研究没有分竹种进行建模，这对模型的精度也会造成一定的影响，如何有效提取竹种信息，以修正和优化竹林碳储量遥感模型有待进一步深入研究。