树轮年代学研究树木年轮的宽度、密度等性质构成的时间序列及其与环境因子的响应关系^[1]，在反演过去和预测未来气候变化中扮演着重要作用^{[2, 3, 4, 5]}。树轮年代学中，去除树木自身的生长趋势并建立树轮年表是开展研究的基础步骤之一。在去趋势处理和建立年表方面，Cook等编写的ARSTAN程序是目前最权威的分析工具，在树轮年代学的发展过程中发挥了巨大作用。然而ARSTAN基于Fortran语言编写，同时缺乏系统而详细的说明文档，使用者不容易理清其内部的算法，实现新的分析方法也比较困难，这在一定程度上限制了研究的深入和学科间的交流。

R语言是当前发展较快的统计分析工具之一，在众多研究领域得到了越来越广泛的应用。Bunn等使用R语言编写了dplR扩展包，实现了树轮数据读取、转换、分析、绘图等功能^[6]。相比于经典工具而言，R语言和dplR等相关扩展包的源代码完全公开，且容易理解，研究者可以根据自身需要修改或编写程序；同时使用者众多，网络社区发达，有利于不同领域的相互交流。因此，R语言及dplR扩展包是ARSTAN等传统程序的良好补充，对于树轮年代学的进一步发展具有重要意义。

目前在国外研究中，dplR被广泛应用于去趋势处理和建立年表^{[7, 8, 9]}，在交叉定年^[10]、小波分析^[3]等方面也发挥了作用；在其基础上又发展出交互式去趋势^[11]、响应函数分析^[12]、扰动事件重建^[13]等相关扩展包，逐渐形成了一套较为完整的树轮数据分析软件体系。而在国内树轮年代学领域，使用dplR开展的研究还较少，使用的功能也相对单一^[14]。作为一个新兴的树轮年表分析工具，dplR处理结果与传统程序一致与否，将直接影响到不同来源研究数据的可比性，因此有必要对ARSTAN和dplR在算法和结果上进行比较分析。

本文使用贺兰山青海云杉树轮宽度数据，分别借助ARSTAN和dplR进行原始序列统计量计算、去趋势处理、年表建立、公共区间分析等操作，通过相关分析、小波分析等手段比较了两者所得结果的差异，并从算法的角度探讨了差异的来源，为开展树轮年代学数据的整合分析提出了需要注意的问题及相应建议。

本文中使用的树轮资料为贺兰山青海云杉(Picea Crassifolia)树轮宽度数据^[15]，采样点位于贺兰山东坡林线附近，共包含来自25棵样树50根样芯的树轮宽度测量结果，数据精度为0.01 mm。使用COFECHA程序检验交叉定年质量，结果显示序列平均长度为83.9 a，平均敏感度为0.393，序列间相关系数为0.825。

ARSTAN程序目前分为Windows XP版和OS X版，早期的DOS版已不再更新。本文主要使用Windows XP版ARSTAN进行去趋势处理并建立年表，程序下载自美国哥伦比亚大学树木年轮实验室(www.ldeo.columbia.edu/tree-ring-laboratory)，版本号为44。选择的去趋势方法为修正的负指数函数，即：

点击浏览公式

R语言主程序和dplR扩展包均下载自综合R档案网络(cran.r-project.org)，其中R语言版本号为3.1.1 ，dplR版本号为1.6.1。使用“rwl.stats”函数计算原始树轮宽度序列和树轮宽度指数年表的基本统计量。使用“detrend.series”函数对原始树轮宽度序列逐一进行去趋势操作，确保各样芯使用的去趋势方法与ARSTAN中相同，使得结果具有可比性。使用“rwi.stats”函数对树轮宽度指数序列进行公共区间分析，时间段同样设定为1966—2009年，这也是程序默认的公共区间。使用“chron”函数采用双权重平均法建立年表，最终得到STD和RES年表。

比较ARSTAN和dplR所得原始轮宽序列统计量的差异，使用Spearman相关系数考察误差与序列长度的关系。对比ARSTAN和dplR拟合生长趋势所得曲线的参数及拟合结果，使用Spearman相关系数考察两程序所得STD、RES年表的相关性，以ARSTAN所得结果为准，计算dplR所得年表的标准误；针对差异相对较大的RES年表，使用Morlet小波变换分析两程序所得RES年表在频域上的差异。最后比较不同程序在年表统计量和公共区间分析结果上的差异。相关分析中，显著性水平设置为0.05。

针对50个原始轮宽序列，ARSTAN和dplR计算所得均值和标准差完全相等；而dplR所得平均敏感度(mean sensitivity,MS)和一阶自相关系数(first-order autocorrelation,AC)均略小于ARSTAN所得结果，其中MS的平均误差为0.005，标准误为0.0002，AC的平均误差为0.008，标准误为0.0003，两者的误差与序列长度存在显著负相关关系(P < 0.05)。进一步计算发现，ARSTAN和dplR所得MS和AC满足：

点击浏览公式

ARSTAN和dplR拟合所得负指数函数的参数a、b完全相同，参数d的误差在0.0003以内，而线性函数参数完全相同。去趋势所得各样芯的两个标准化序列间存在差异的年份不超过10年，且误差均在0.003以内。在此基础上，两程序建立的STD年表相关系数为0.999(P < 0.05)，对应年份宽度指数误差在0.05以内，平均误差为0.005，标准误为0.0010，相位变化情况基本相同(图1)。

图 1 ARSTAN和dplR所得树轮宽度指数年表 Fig.1 Tree-ring chronologies produced by ARSTAN and dplR (a)为标准年表；(b)为差值年表；(c)为不同程序建立的年表之差

图选项

在建立差值年表时，ARSTAN对各标准化序列均拟合为一阶自回归模型，所得差值序列与原序列长度相等；而dplR则拟合成阶数各异的自回归模型，包括29个一阶模型和17个二阶或更高阶模型，序列初始损失的年份数与模型阶数相等，另有4个样芯未拟合成自回归模型。各样芯的两个差值序列间的平均误差为0.034，标准误为0.0011。两程序建立的RES年表间相关系数为0.987(P < 0.05)，在1914年前(即序列前30年)，两者平均误差为0.048，标准误为0.0096；而1914年后，两者平均误差为0.028，标准误为0.0021。两个RES年表仅在起始的10年内存在相位不一致的情况，例如1887—1889年(图1)。

对两个RES年表分别进行Morlet小波变换(图2)。结果显示，两个RES年表含有的高频信息基本相同，均表现出1940年前2—11a尺度的震荡能量较小，周期变化不明显，在1940年后2—11a尺度的能量增强，周期变化显著；而在低频信息方面，ARSTAN建立的RES年表在1940年后稳定存在着32a左右的周期，而dplR建立的RES年表32a左右的周期震荡能量始终较小。

图 2 差值年表小波功率谱 Fig.2 Wavelet power spectrum of the residual chronologies

图选项

年表统计量及公共区间分析结果见表1。在年表统计量方面，ARSTAN的STD、RES年表和dplR的STD年表的序列起止年相同，而dplR的RES年表与这3个年表相比少1年(1884年)。两个STD年表的均值、标准差均相差0.001；MS相差0.003，相对误差约为1%；AC相差0.005，相对误差约为1.5%。而两个RES年各项统计量的误差均大于STD年表，其中dplR的RES年表MS比ARSTAN低4%，而AC的绝对值低36%。

在公共区间分析方面，两个STD年表的所有样芯平均相关系数(all-series Rbar,Rtot)、树内平均相关系数(within-trees Rbar,Rwt)、树间平均相关系数(among-trees Rbar,Rbt)和年表有效信号(effective chronology signal,Reff)均完全相等，而ARSTAN的RES年表在这些指标上普遍略高于dplR的RES年表。dplR的STD和RES年表的样本总体代表性(expressed population signal,EPS)比ARSTAN计算结果低0.004，相对误差约为0.4%，而信噪比(signal-to-noise ratio,SNR)结果低30左右，相对误差高达30%—40%。

根据原始轮宽序列统计量比较结果，ARSTAN与dplR在MS和AC的计算上存在着较为确定的换算关系，这很可能是由于两者使用了不同的计算公式。dplR计算MS的公式为^[16]：

点击浏览公式

dplR调用R中的“acf”函数计算AC，其计算公式可表达为^{[1, 17]}：

点击浏览公式

点击浏览公式

ARSTAN和dplR均具有多种去趋势方法，两者的对应关系如表2所示，其中在拟合负指数函数的程序实现方面，两者具有较大差异。在算法上，ARSTAN使用Fritts等专为负指数函数拟合问题设计的算法^[18]，该方法充分考虑到负指数函数的曲线特征，首先估计曲线的切线斜率b，通过调节b的值进行迭代。而dplR首先估计y轴截距a和渐近线d的初始值，b的初始值则根据经验直接设为0.01，再调用R内置的“nls”函数进行曲线参数估计。如果不设定约束条件，当树轮宽度序列的趋势较符合负指数函数曲线时，两种算法均能得到精确的拟合结果；而当树轮宽度序列的趋势近似于直线时，d为负值，这时dplR的初始值设定严重偏离最优解，导致迭代过程报错，转而使用线性函数进行拟合。如果设定约束条件d>0，ARSTAN首先在整个定义域范围内寻找各参数的最佳估计值，然后将参数带入约束条件中进行判断，如果不符合条件则拒绝拟合负指数函数，转而拟合线性函数；而dplR则是在“nls”函数中缩小参数的取值范围，最终拟合得到参数的局部最优解。当拟合方法出现分歧时，去趋势所得序列将出现较大差异，而本文在操作中参考ARSTAN的运行结果，使用dplR对每根样芯单独去趋势，确保两者使用的去趋势方法相同。

尽管拟合生长趋势的方法相同，ARSTAN和dplR建立的标准年表仍存在一定差异，这与缺失值的处理方式有关。为防止运算中可能出现的错误，dplR将序列中的缺失值(一般赋值为0)均替换成0.001，因此dplR所得标准化序列不存在树轮宽度指数为0的情况；而ARSTAN并未对缺失值进行特殊处理，因此标准化序列中仍存在0值。

建立差值年表，需要去除树轮宽度指数标准化序列中的自回归成分。ARSTAN首先用所有标准化序列拟合出一个p阶合并自回归模型，将拟合值作为树木群体所共有的持续性生长量。再对各标准化序列分别拟合限定为p阶的自回归模型，将拟合值作为整体生长模式下个体的持续性生长量，而将残差作为差值序列，并使用双权重平均法建立RES年表。如果将群体的持续性生长量加回到RES年表上，则得到ARS年表^[19]。dplR未考虑整体的生长模式，而是直接对各标准化序列分别拟合自回归模型，不同序列的模型可以具有不同的阶数，再将残差作为差值序列，建立RES年表。针对本文所用数据，dplR的计算过程对个体特有的持续性生长量剔除得更为彻底，低频信息保留得更少，这种差异在频域上更为明显(图2)。

拟合p阶自回归模型会造成序列前p年树轮宽度指数的损失。针对这一情况，ARSTAN使用序列均值将标准化序列向前延长p年^[19]，从而使自回归模型计算出的拟合值序列与原序列等长，避免了序列长度上的损失。而dplR未对初始年份的损失进行特殊处理，因此两个RES年表存在的差异相对STD年表更大，尤其是样本量较小的最初若干年。根据本文结果，起始值处理方式造成的差异会在约30a后消失，同时考虑到参与“树轮-气候”关系分析的年表区间通常远小于年表总长度，如果起始年份不参与后续的分析，则起始年份的不同处理方式对最终结果的影响较小^{[19, 20]}。

公共区间分析是指选择包含特定时间范围的树轮宽度序列，通过计算年表信号，评价最终所得年表对森林总体生长的代表性，其中年表信号是衡量序列中包含的共同变化的统计量^[21]。DOS版ARSTAN在导入测量数据时会提供一个推荐的公共区间范围，其选择标准是公共区间内包含尽可能多的测量值，即满足区间内完整样芯数与区间长度乘积取得最大值；而Windows XP版ARSTAN不再提供这一结果。dplR中对数据集进行去趋势使用“detrend”函数，其默认公共区间是所有样芯共有的公共区间，该区间长度一般小于DOS版ARSTAN提供的公共区间长度。dplR中的“common.interval”函数可以给出各种标准对应的时间范围。

在进行公共区间分析时，最初使用方差分析对年表信号进行估计^[1]，后来的研究发现使用相关矩阵也可以估计年表信号的大小^[22]。以序列间的平均相关系数(mean series intercorrelation,Rbar)作为年表信号的衡量指标，可以通过计算EPS和SNR来评价年表的质量^[17]：

点击浏览公式

根据运算结果推测，ARSTAN可能使用样芯数量m和Rtot计算EPS和SNR。这种算法将所有样芯等同看待，没有剔除同一样树的不同样芯包含着的重复信息，因此可能高估了年表质量。EPS除了用于评价年表整体质量外，还用于计算子样本信号强度(subsample signal strength,SSS)^[21]，而EPS和SSS均是衡量年表可靠长度的重要指标。由于ARSTAN计算得出的EPS偏高，可能会令EPS或SSS高于阈值(一般设为0.8或0.85)的最小样本量偏低，造成对年表可靠长度的高估。

目前的树轮年代学研究中，一棵样树往往取两个甚至多个样芯，因此dplR采用Reff计算EPS和SNR，计算结果更符合现有理论。使用基于Reff的公式重新评估本文中ARSTAN建立的年表，则STD年表EPS为0.986，SNR为72.276，而RES年表EPS为0.989，SNR为90.741，这些结果均与dplR所得结果接近。

随着树木年轮基础数据日渐丰富，对前期样点尺度的研究成果进行整合，开展区域、大陆甚至全球尺度上森林对于全球气候变化的响应分析已成为研究热点之一^{[23, 24]}。根据本文所得结果，在树木年轮研究中如果使用不同程序处理原始数据，所得年表及其质量评估需要考虑算法差异可能造成的影响。在标准年表建立方面，即使ARSTAN与dplR拟合生长趋势的参数设定相同，拟合结果可能差异很大。而对于差值年表，在频域上，ARSTAN建立的RES年表比dplR具有更多的低频信息；在时域上，起始年份处理方式的不同，造成年表在最初的20—30年内存在较明显的差异。不同程序在MS、AC等指标上可能存在计算方法差异；ARSTAN使用Rtot计算EPS和SNR等指标，相对于文献中基于Reff的算法而言，可能高估了年表中包含的共同信息^{[25, 26]}。

综合以上讨论，对开展不同来源树轮数据的整合分析提出如下建议：

(1)如果能够在国际树轮数据库等共享平台上获取原始数据，建议使用同一程序完成年表建立和公共区间分析等处理步骤，并注意检查各样芯去趋势方法的具体参数。

(2)如果仅能获得文献资料，建议记录其处理过程的详细信息，包括软件的版本信息、去趋势的具体方法等；当信息足够时，建议根据同一标准设定公共区间，并重新计算MS、EPS等指标，确保结果的可比性。

在树轮年代学研究中，新兴的R语言dplR扩展包是对ARSTAN等传统程序的良好补充，其平台开放性可以促进研究者对程序算法的认识和跨学科的相互交流。本文对比了ARSTAN程序和dplR扩展包，分别使用两种程序对贺兰山青海云杉树轮宽度数据进行原始序列统计量估计，完成了去趋势处理并建立了标准年表和差值年表，通过公共区间分析评价了年表的质量。分析表明，ARSTAN与dplR在拟合生长趋势和自回归模型时存在算法上的较大差异，同时平均敏感度、一阶自相关系数、样本总体代表性和信噪比等指标的算法也不一致，但存在较为确定的转换关系。整合分析不同来源的树轮资料时，应明确其数据处理过程，在条件允许的情况下使用同一程序或算法进行去趋势处理和公共区间分析。