利用树轮资料和植被指数估测青海云杉林地上净初级生产力

范文武，凌飞龙，田 昕，汪小钦，闫 敏

(1. 福州大学空间数据挖掘与信息共享教育部重点实验室，福建 福州 350002； 2. 中国林业科学研究院资源信息研究所，北京 100091; 3. 重庆市地理信息中心，重庆 401120)



利用树轮资料和植被指数估测青海云杉林地上净初级生产力

范文武1, 2, 3，凌飞龙1，田 昕2，汪小钦1，闫 敏2

(1. 福州大学空间数据挖掘与信息共享教育部重点实验室，福建 福州 350002； 2. 中国林业科学研究院资源信息研究所，北京 100091; 3. 重庆市地理信息中心，重庆 401120)

利用MODIS产品MOD13A3、 样地调查数据和树轮资料等数据，构建回归模型以估测祁连山森林保护区的青海云杉林地上净初级生产力(ANPP). 结果表明： 经质量权重S-G滤波后，回归模型的估测精度有所提高； 2000-2013年祁连山森林保护区青海云杉林ANPP呈现波动上升的趋势，空间分布特征为东南多西北少，多年平均森林ANPP为51.34～274.24 g·m-2； 青海云杉林ANPP与降水和温度呈正相关，受温度变化影响较小，降水成为影响其变化的主要因子.

青海云杉； 地上净初级生产力； 增强植被指数； 树轮资料

0 引言

森林地上净初级生产力ANPP(aboveground net primary production)是地表碳循环的组成部分，准确估测森林ANPP的变化是全球碳循环研究的重要内容. 树木年轮数据具有分辨率高、 时间序列长等优势[1]，能获取长时间序列的森林ANPP，但在大尺度森林ANPP估测方面存在局限. MODIS数据具有宏观、 大面积监测等优点，其光谱信息与森林ANPP具有相关性. 其中，NDVI(normalized difference vegetation index)能表述森林碳通量的变化，但在浓密的常绿林中受到限制; 而EVI(enhanced vegetation index)避免了植被指数的饱和问题，减少了大气和土壤背景的影响.

以祁连山森林保护区为研究区，利用植被指数产品MOD13A3、 树轮资料等数据，构建iEVI(integrated EVI，集成增强植被指数)[2]和森林ANPP回归模型，并对模型可靠性进行验证，分析气候因子对青海云杉林空间分布的影响.

1 研究区概况与数据

1.1 研究区概况

黑河流域发源于祁连山北麓，位于河西走廊中部，流域范围介于96°42′～102°08′E，37°44′～42°42′N之间，流域面积约13×104km2. 本研究以黑河流域上游祁连山森林保护区为研究区，该地区林分是以青海云杉(Piceacrassifolia)为主的天然次生林，其中夹杂着少量的垂枝祁连圆柏(Sabinaprzewalsk). 青海云杉主要生长于祁连山北麓阴坡地带，海拔在2 500～3 300 m之间.

1.2 野外采样数据

2014年5月在祁连山森林保护区选取22块样地，大小为20 m×20 m，样地内胸径大于5 cm的青海云杉采用每木检尺的方式测量胸径、 树高. 样地内按6个径级(5 cm≤D<10 cm，10 cm≤D<15 cm，15 cm≤D<20 cm，20 cm≤D<25 cm，25 cm≤D<30 cm，30 cm≤D，D为胸径), 每径阶选取三株标准木，每株标准木在胸径处(1.3 m)沿东西、 南北方向钻取树芯，共钻取792颗树芯. 将树芯保存在塑料管中密封，在实验室用Lintab年轮分析仪器测量树芯的长度，树木年轮观测的精度为0.001 mm[3].

1.3 遥感数据及森林/非森林分类图

本研究采用的遥感数据为EOS/Terra卫星的植被指数产品MOD13A，该产品包含EVI和植被指数质量标记(QA)等数据，其空间分辨率为1 km，时间分辨率为30 d，每年12景(2000年11景)，时间序列为2000-2013年，共计167景影像. 本研究将2000-2013年每年12景EVI数据(2000年11景)分年累加求和得到iEVI(integrated EVI)数据.

EVI=2.5×

式中：ρblue、ρred和ρNIR分别为MODIS影像蓝、 红、 近红外波段的反射率.

QA数据在像元尺度对植被指数产品MOD13A进行质量评价，结果见表1. 表1中第二字段“VI可用性”分16个等级，其中1代表植被指数数据质量最高，随着数值增大质量降低. 该字段16个等级是根据获取EVI时的气溶胶含量和云条件等因素评定的结果. 因此，本研究将该字段作为质量因子对EVI进行滤波.

表1 MODIS VI质量数据域

青海云杉林分布图采用郭云等[4]的研究结果(分类精度为90.39%)，并使用最近邻重采样方法将分类数据升尺度至1 km分辨率. 分类图中TM影像来源于美国地质勘测局，4景TM影像成像时间为2009年8月11日(1景)和2009年7月17日(3景).

1.4 气象数据

黑河流域2000-2013年大气驱动数据集来源于寒区旱区科学数据中心[5]，包含2000-2013年逐时0.05 ℃的2 m气温和累积降水等数据，采用WRF(weather research and forecasting)模式制备，经过22个站点观测资料进行逐日和逐时验证. 所使用气温和降水等数据为掩膜后的青海云杉林区域的数据.

2 研究方法

2.1 质量权重Savizky-Golay滤波

表2 QA可用性指数分级赋权

图1 研究区单一像元EVI时间序列及其对应的QA质量标记 Fig.1 A single pixel EVI time-series and the corresponding QA

Savizky-Golay(S-G)滤波是一种基于最小二乘的卷积运算以平滑时间序列为数据的滤波方法[6]. 质量权重S-G滤波方法是从MODIS QA数据域中提取出质量因子，将其转换为质量权重以替换自适应S-G滤波迭代过程中的距离权重系数. 权重设置见表2. 对MOD13A数据利用LDOPE质量分解软件得到QA数据，用MRT投影转换软件提取EVI数据和QA数据集中的16级EVI可用性指数，建立连续14 a的EVI和对应QA时间序列. 图1为单一像元EVI时序数据，横轴为EVI数据时间序列，每一期隔30 d，连续14 a的期数. 红圆点为该像元所对应的QA质量标记在滤波处理中的权重等级，圆半径越大，质量越高. 本研究质量权重S-G滤波方法所使用的系数为： 滤波半窗宽度m=6，迭代次数k=3，上包络线拟合强度为2，算法在timesat[7]环境下实现.

2.2 森林地上净初级生产力估测

表3 青海云杉各部分生物量与树高胸径的关系

森林ANPP是生物量的增长量(ΔW)、 凋落枯死量(L)和动物采食量(G)的总和[8]. 由于青海云杉位于祁连山国家自然保护区，地处高寒山区，受人为和动物干扰影响较小，而且森林冠层郁闭以后，凋落量基本恒定. 本研究忽略凋落枯死量和动物采食量，将地上生物量的年增长量近似表示森林ANPP，即ANPP≈ΔW，第N年青海云杉林的ANPPN≈WN-WN-1. 通过样地调查数据建立树高与胸径模型，借助树轮逐年宽度，得到逐年胸径和树高. 本研究结合文[9]对祁连山青海云杉各器官生物量生长方程W=a(D2H)b(表3)(D表示胸径，H表示树高，a和b为系数)，估算每径阶标准木的地上生物量，即青海云杉干、 枝、 叶三部分的生物量，然后分别乘以样地内各径阶青海云杉的株数，累加求和得到样地地上生物量，最终将青海云杉林地上生物量的增长量近似表示逐年的ANPP.

2.3 回归模型的建立与验证

回归模型经常用于建立森林ANPP与各种实测树木生物物理参数之间的相关关系，回归方程是森林ANPP估算常用的统计方法. 研究中将样地观测数据按2/3和1/3比例随机选取，分别用于建立模型和精度验证. 在建立模型时，以样地的iEVI为自变量、 ANPP为因变量，建立线性或非线性模型，以模拟的ANPP与观测的ANPP拟合得到的R2(决定系数)和RMSE(均方根误差)为指标选择模型.

其中:xi表示第i个估测数据；yi表示第j个模型数据；R2表示表征依变数y的变异中有多少百分比可由控制的自变数x来解释； RMSE表示观测值与真实值的偏差.

3 结果与分析

3.1 滤波前后EVI数据对比

为了直观地对比质量权重S-G滤波前后对影像噪声的滤除程度，选取图1中第4个波谷对应的第36景(2003年1月1日-31日合成)噪声含量很高的影像进行对比评价(如图2所示). 由于大多噪声对应较低的EVI值，由图2(a)可知该时相影像含有大区域的噪声(黑色区域)，图2(b)清晰地显示出质量权重S-G滤波方法能够滤除绝大部分的噪声.

图2 高噪声EVI影像滤波前后对比Fig.2 An original strong noisy image and images filtered with quality-weighted S-G

为了定量分析滤波前后对不同质量因子EVI数据的影响，本研究选取2000年9月、 2004年1月、 2008年9月和2012年1月的影像进行统计. 图3显示了不同质量EVI数据在滤波前后的相关系数，质量因子为1的高质量EVI数据相关系数接近于0.9，说明对高质量EVI数据的规律性影响较小. 随着质量因子变大(噪声增强)，相关系数都有减小的趋势. 因此，质量权重S-G滤波方法能够有效地滤除噪声.

3.2 树高胸径模型

利用样地调查的树高和胸径数据进行曲线拟合，选取多个拟合模型进行分析评价，以最优决定系数R2、 最小均方根误差RMSE作为评价标准，确定胸径与树高最佳相关模型H=0.89D0.85，R2=0.88，RMSE=2.06 m(图4). 该拟合方程能较好地解释树高和胸径的关系，借助树轮宽度数据可反演样地标准木逐年树高.

图3 不同质量EVI数据滤波前后相关性

图4 树高和胸径的数据拟合曲线

3.3 iEVI与ANPP回归模型

利用树轮资料估算2000-2013年青海云杉22块样地ANPP，随机抽取14块样地数据用于构建ANPP与iEVI的回归模型，余下8块样地数据用来验证模型精度. 从建模样本来看，质量权重S-G滤波前后iEVI与ANPP建模精度R2分别为0.53(图5(a))和0.63(图5(b))，建模效果较好. 回归模型精度验证表明，滤波前回归模型的估测精度为R2=0.48，RMSE=36.39 g·m-2(图5(c))，滤波后回归模型的估测精度为R2=0.60，RMSE=34.82 g·m-2(图5(d))，质量权重S-G滤波后模型估测精度有所提高. 从点的分布来看(图5(d))，ANPP主要集中分布在140 g·m-2附近的低值区域，高值区域分布较少且零散，当实测值ANPP较低时，模型结果偏高，实测值ANPP较高时，模型结果偏低. 林分密度是影响森林地上生物量的重要因素，利用树轮数据反演过去林分地上生物量，尚无法估算林分密度在森林发育过程中的变化. 青海云杉林密度随林龄增加而下降，这与人工抚育和林分自疏有关，对于青海云杉林，以当前密度反演过去森林地上生物量可能会导致对林分恢复早期地上生物量的低估[2]，进而影响对森林ANPP的估算.

根据2000-2013年逐年的ANPP与质量权重S-G滤波后的iEVI数据建立对数回归模型，得到逐年的ANPP反演模型(图6). 各年份的iEVI和ANPP建模精度R2呈现差异，2008年iEVI和ANPP拟合效果较差(R2=0.49)，2005年iEVI和ANPP拟合效果最好(R2=0.68).

图6 2000-2013年iEVI与ANPP的相关性Fig.6 The correlation between iEVI and ANPP from 2000 to 2013

3.4 青海云杉林ANPP时空变化特征

通过iEVI与ANPP逐年的回归模型，估算祁连山森林保护区2000-2013年的青海云杉林ANPP平均值，得到近14 a来森林ANPP的变化情况(图7). 从图7中可以看出，2000-2013年平均ANPP在106.51～161.96 g·m-2之间; 2000-2005年间除ANPP在2001年到达最小106.51 g·m-2，其它年份均值在125.76～134.35 g·m-2之间小幅摆动； 2006-2010年ANPP则呈现出明显的上升趋势，2009年ANPP达到最大值161.96 g·m-2； 2011-2013年ANPP则呈现下降的趋势. 总体上，青海云杉林ANPP呈缓慢波动增加趋势.

如图8所示青海云杉林主要分布在祁连山区北麓阴坡地带，总面积约为777 km2，2000-2013年，年平均森林ANPP为51.34～274.24 g·m-2. 森林ANPP值在180～274.24 g·m-2之间的面积为129 km2，约占总面积的16.60%，分布在肃南县和祁连县东南部、 民乐县和山丹县地区； 森林ANPP值在138～180 g·m-2之间的面积为212 km2，约占总面积的27.28%，分布在肃南县中部和东南部、 祁连县东南部以及民乐县地区； 森林ANPP值在102～138 g·m-2之间的面积为257 km2，约占总面积的33.08%，分布较为零散； 森林ANPP值在51.34～102 g·m-2之间的面积为179 km2，约占总面积的23.04%，分布在肃南县中部和西北部地区； 森林ANPP总体空间分布特征大致为东南多西北少.

图7 祁连山森林保护区青海云杉逐年平均ANPP

图8 祁连山森林保护区青海云杉年平均ANPP空间分布

3.5 青海云杉林ANPP与气候因子的相关性

2000-2013年青海云杉林区域年降水总量在278.82～671.01 mm之间(图9(a))，降水的年际变化与森林ANPP的年际变化总体趋势较为相似，2000-2010年的增减具有同步性，2011-2013年呈现相反的趋势，但年降水总量整体呈现上升的趋势. 2000-2013年平均温度在-4.53～4.78 ℃之间(图9(b))，气温的年际变化与森林ANPP的年际变化在2000-2004年表现出相反的趋势，2005-2013年表现为一致的趋势，年平均温度整体呈现下降的趋势，整体相关性较小.

图9 青海云杉林ANPP与降水、 气温的年际关系Fig.9 Correlations between annual forest ANPP of P.crassifolia and temperature and precipitation

对青海云杉林区域2000-2013年年降水总量、 年平均气温以及平均森林ANPP进行相关性分析，森林ANPP分别与年降水总量(R=0.49，P<0.05)和年平均气温(R=0.22，P>0.05)呈正相关. 青海云杉林地处高寒地区，温度变化对森林ANPP影响较小，降水成为影响森林ANPP变化的主要因子. 2000-2004年气温较高，热量水平的提高加大了区域水分的蒸发，且降水量较少，在一定程度上限制了森林ANPP的增长. 森林ANPP整体上随降水量的增加而波动增长，表明降水有利于加强森林的固碳能力.

4 结论

本研究在考虑MODIS数据质量的前提下，对EVI数据进行质量权重S-G滤波，使得高质量EVI数据保真性较好，黑色噪声区域得到有效滤除. S-G滤波前后模型反演精度有所提高，但从模型验证精度来看，ANPP主要集中在140 g·m-2附近的低值区域，高值区域分布较少且零散，模型预测结果偏小.

通过长时间序列的树轮资料反演森林ANPP，并建立植被指数iEVI和森林ANPP的回归模型，实现样地尺度到区域尺度的转化. 从回归模型验证精度(R2=0.60，RMSE=34.82 g·m-2)来看，该方法是可行的. 但受到空间异质性等影响，后期需对MODIS EVI进行混合像元分解或使用30 m分辨率的Landsat TM数据计算EVI，与样地数据更好的空间匹配，提高模型反演精度.

时间特征上，2000-2013年青海云杉林ANPP在2001年为最低值106.51 g·m-2，2009年达到最高值161.96 g·m-2，在14 a的时间中增加趋势较为缓慢，这与其它学者对祁连山ANPP变化研究得出的结论基本一致[10]. 空间特征上，青海云杉主要分布在祁连山区北麓阴坡亚高寒山区，总体空间分布特征为东南多西北少. 青海云杉林ANPP受温度变化影响较小，而降水成为影响森林ANPP变化的主要因子. 降水量的增加，为植被的生长提供了必需的水分，使得青海云杉林固碳能力有所加强.

[1] 何吉成. 用树轮宽度资料重建长白山森林净初级生产力[D]. 北京： 中国科学研究院地理科学与资源研究所，2005.

[2] CAMPOS G E P, MORAN M, HUETE A,etal. Ecosystem resilience despite large-scale altered hydroclimatic conditions[J]. Nature, 2013, 494(7 437): 349-353.

[3] HOLMES R L. Computer-assisted quality control in tree-ring dating and measurement[J]. Tree-Ring Bulletin,1983, 43(3): 69-75.

[4] 郭云，李增元，陈尔学，等. 甘肃黑河流域上游森林地上生物量的多光谱遥感估测[J]. 林业科学，2015，51(1)： 40-149.

[5] 赵俊芳，延晓冬，朱玉洁. 陆地植被净初级生产力研究进展[J]. 中国沙漠，2007，27(5)： 780-786.

[6] 周增光，唐娉. 基于质量权重的Savitzky-Golay时间序列滤波方法[J]. 遥感技术与应用，2013，28(2)： 232-239.

[7] JÖNSSON P, EKLUNDH L. TIMESAT: a program for analyzing time-series of satellite sensor data[J]. Computer Geoscience,2004, 30(8): 833-845.

[8] PASTOR J, POST W M. Response of northern forests to CO2-induced climate change[J]. Nature,1988,334(6 177): 55-58.

[9] 王金叶，车克钧，傅辉恩，等. 祁连山水源涵养林生物量的研究[J]. 福建林学院学报，1998，18(4)： 319-323.

[10] 张福平，冯起，李旭谱，等. 黑河流域NPP遥感估算及其时空变化特征[J]. 中国沙漠，2014，34(6)： 1 654-1 664.

(责任编辑： 林晓)

Estimation ofPiceacrassifoliaforest aboveground net primary productivity from tree ring data and vegetation index

FAN Wenwu1, 2, 3, LING Feilong1, TIAN Xin2, WANG Xiaoqin1, YAN Min2

(1. Key Laboratory of Spatial Data Mining & Information Sharing of Ministry Education, Fuzhou University, Fuzhou, Fujian 350002,China； 2. Institute of Forest Resources Information Techniques, Chinese Academy of Forestry, Beijing 100091, China; 3. Chongqing Geomatics Center, Chongqing 401120, China)

The MOD13A3 EVIproducts, forest inventory data and tree ring data were used to estimateANPP ofPiceacrassifoliain nature reserve of Qilian Mountainby the regression model.The performance of iEVI in estimating ANPP was improved after the quality weighted S-G method. The annual forest ANPP showed a fluctuant increasing trend from 2000 to 2013 in the nature reserve of Qilian Mountain. Furthermore, the ANPP increased from the northwest to the southeastern and the annual average ANPP was 51.34～274.24 g·m-2. The ANPP ofPiceacrassifoliawaspositively correlated with precipitationand temperature. Compared with the temperature, the precipitation is main factor affecting the change of ANPP.

Piceacrassifolia; aboveground net primary productivity; enhanced vegetation index; tree ring data

10.7631/issn.1000-2243.2017.03.0348

1000-2243(2017)03-0348-07

2015-10-19

田昕(1979 -)，副研究员，主要从事遥感技术与应用方面研究，tianxin@ifrit.ac.cn

中央级公益性科研院所基本科研业务费专项资金资助项目(IFRIT201302)； 国家重点基础研究发展计划资助项目(2013CB733404)

S757

A