张永

（中国科学院地理科学与资源研究所，陆地表层格局与模拟重点实验室，北京 100101）

经验模态分解方法对祁连圆柏年轮宽度序列标准化初探

张永

（中国科学院地理科学与资源研究所，陆地表层格局与模拟重点实验室，北京 100101）

利用经验模态分解（EMD）方法对祁连圆柏树轮宽度序列进行标准化处理的初步探讨，分析该方法在圆柏宽度序列生长量订正中的应用潜力。利用EMD方法对祁连山地区的65个圆柏样芯宽度序列的逐级分解，较好地展示了各样芯由高频至低频的自然变化特征，其中最低频的趋势项虽包含了线性、负指数、类样条函数等多种形式，但末端基本趋于平稳，与树木的生物学遗传特性较为吻合，基本反映了与树木遗传相关的生物学趋势。与传统线性或负指数函数方法对祁连圆柏序列的标准化过程相比，传统处理中常常出现髓心端上翘或树皮端异常增大的问题，通过新的EMD方法的订正得到了改进，进一步分析了EMD方法在祁连圆柏轮宽序列生长量订正中的优势及不足。两种方法得到的祁连圆柏年表在整个时段上高频变化极为一致，仅在接近髓心端的100a内传统方法处理得到的指数值偏低，方差变化较大；而EMD处理得到的指数年表整个时段方差相对稳定。

树轮宽度序列；标准化方法；经验模态分解（EMD）

气候变化是影响世界环境、人类健康与福利和全球经济持续性发展的重要因素，了解气候变化的规律及提高对未来气候变化的预估水平成为现今地球科学最重要的研究课题之一。作为过去气候环境变化的天然信息库，树木年轮以其定年准确、分辨率高、连续性强、空间分布广、环境变化指示意义明确且可精确定量化等优势成为研究过去气候变化的一种重要代用指标[1-10]。

然而，树木年轮的形成与变异除了取决于气候和其它环境因子以外，还受树木本身的遗传因子控制。树轮的生长可以用线性聚类模式[11]表示如下：

Rt=At+Ct+△D1t+△D2t+Et（1）

其中Rt表示年轮宽度；At表示生长趋势；Ct表示气候信息；△D1t表示内部干扰；△D2t表示外部干扰；Et表示其它。

而我们主要研究环境因子中的气候因子对树木生长的影响，因此，对实际测量得到的年轮宽度值，在分析其能否反映出某些气候要素的变化之前，应设法消除树龄对年轮宽度的这种缓慢影响，即树木本身支配其生长的因素（遗传因子）应被剔除，从而得到仅由气候因子影响所决定的年轮序列变化量，该过程被称为生长量订正。然而，树龄的影响不像气候因子那样表现在逐年的轮宽相对增量上，而主要表现为缓慢的趋势性变化。一般来说，多数树木的年轮宽度均有下述变化趋势：幼龄期年轮宽度比较窄，后随树龄的进一步增加和树干直径增大的几何效应，到达树木生长最旺盛的时期时，平均的年轮宽度常为整个生长期内的极大值，此后，树木的年轮宽度随树龄的增加而趋于减小，减小的速度开始较快，后来越来越慢，到一定时期就大体趋于平稳，被称为树木生长趋势[12-13]。

如何去除树木生长趋势的影响，对树轮资料进行处理时存在的去趋势偏差、低频信号的提取和保留等也是树木年轮气候学研究中的一个重要部分[12]，这限制和影响着树轮资料在研究气候变化时的准确性。

1 传统的树轮生长量订正方法

精确地确定出树木的生长趋势并作出准确的生长量订正是一项十分困难的事情，这是由于随着树种和环境的不同，其树木的群体生长趋势也有较大差异。此外，生长量订正的另一目的是要把树轮宽度资料整理成可比较的形式。长期以来，为了达到上述两个目的，树木年轮气候学者提出了许多生长量订正的方法和手段[4]，设法消除这种影响，被称为年轮宽度序列标准化[13-14]。根据生长趋势曲线估算的思路和方法大体上可归为3类：树轮解析估计法、大样本估算法、多种函数拟合法[13-14]。

随着树轮学者对树轮去趋势理论认识的提高及经验的积累，多种函数逐步被尝试用来表达树木年轮随树龄变化的生长趋势曲线。其中主要包括线性函数、双曲线函数、指数函数、多项式函数及样条函数等。但每种函数也都有自己的局限性，其中目前以线性或负指数函数应用最为广泛，通常也被认为是比较适用于干旱、半干旱区树种生长趋势拟合的方法[15-16]。特殊情况会考虑样条函数，样条函数[15-16]一般用来拟合较湿润环境下树木生长的持续性以及由树木相互竞争产生的非同步扰动。应用中发现负指数和直线拟合在树轮曲线两端的拟合常出现较大偏差，样条函数一定程度上可以适应环境干扰信号的拟合，但当样条函数的步长选择不当时，虽然有可能看起来与原序列的拟合效果甚佳，但此时可能不仅消除了树木自身的生长趋势，也将受长期气候影响所造成的树木年轮宽度缓慢变化滤掉，造成“信号”或“噪声”重叠的低频部分的“过度拟合”，很难判定最终得出的气候信息是否完整可靠，所以在选用其进行生长量订正的过程中有一定争议。

近年来，研究中常有报道，即使采样点的树木生长对温度敏感，却依然很难重建出中世纪暖期（MWP）、小冰期（LIA）等近千年中典型气候事件，其主要原因就是去趋势的问题，如Mann等[8]利用树轮重建的北半球过去近千年来的温度距平序列，因中世纪暖期不明显而受到较大争议。就此问题，国外学者Briffa[18-19]，Esper[20]开展了关于树轮低频信息提取方面的研究，提出了区域曲线标准化（Regional Curve Standardization,RCS）方法，但局限性在于该方法必须保证对样本总体生长曲线（RC）有较好估计的情况下方可使用，否则会有较大偏差，另外对于单个样芯真实生理年龄的估计方法也存在较多争议。Melvinand Briffa[21]建立了“signal-free”方法用于解决十年至百年尺度的异常变化对传统的树木生长拟合曲线的影响，但也同样存在低频变化的缺失以及使用不同拟合方法导致的不确定性的问题；国内的树轮学者还提出了总体生长趋势曲线拟合法[22]，尝试进一步完善去趋势过程中低频信息的保留以及树木生长总体趋势的估计，但目前仍处于初步研究阶段，未推广应用。

2 经验模态分解（EMD）方法与圆柏生长量订正

2.1 EMD方法介绍

经验模态分解(Empirical mode decomposition, EMD)[14]，是20世纪末由Huang提出来的一个信号分析方法，EMD分解的主要思想是利用波动上、下包络的平均值去确定“瞬时平衡位置”，进而分解IMF分量，每次提取都将高频信号分离出，剩下低频信号，本质上是一个过滤器。按照它所确定的程序，通过逐步筛选，把一个混合信号序列（或其导数，视所需的分解精度而定）分解成有限个更为简单的反映某一个物理信号在窄波段上随时间变化的IMF分量。这些分量保持了原来信号的自然振荡特性，并且是相互独立的。当所有的IMF分量被提取完毕之后，剩余项被称为原始信号的趋势项。对于一个给定的时间序列x（t），可以通过下面的基本步骤完成它的分解：

（1）找出x（t）所有的极大值点并将其用三次样条函数拟合成原始数据的上包络线；类似的找出极小值点拟合成下包络线；并分别记为u（t）、v（t）和m（t），其中m（t）满足：

（2）计算x（t）和m（t）之差，并记为：h1（t）=x（t）-m（t）。

如果m（t）等于零，或者满足设定的精度要求，那么，h1（t）就给出了第一个固有模态函数，并记为IMF1；否则，重复步骤（1）和（2）。

（3）用残差r1（t）=x（t）-h1（t）代替原始信号x（t），并重复上面的筛选过程，可以得到第二个IMF2。

（4）重复步骤（3），依次得到ri（t）和hi+1（t）（i= 2，3，…，n）。

如果rn-1（t）的极值点不超过2个，则整个分解过程结束。EMD分解最终把原序列分解成有限个（n个）IMF分量，不同层次得IMF分量都可能对应某一物理背景。因此，对于一个混沌序列，找出其主要特征量即对应主要物理背景的特征层次，对混沌序列的预测具有重要意义。而原始序列就可以由这些IMF分量及均值之和表示：

分解后所得到每个IMF分量具有如下特征：（1）从全局特性看，极值点数必须与过零点数一致或者至多相差一个；（2）在某一个局部点，极大值包络和极小值包络在该点的算术平均和为零。为了分辨分解结果中可能的虚假分量，可以通过计算IMF和细节分量与原序列的相关系数，得到与原序列显著相关的IMF和细节分量。

2.2 EMD对祁连圆柏轮宽序列的分解

由上文树轮宽度序列的聚类模式表达式可以看出，树轮宽度序列正是一个多信号混合的混沌序列，为了分析其中最主要的气候因子对轮宽变化的影响，借助于EMD分解的方法，将各信号逐级进行分离，去除与年龄相关的生物趋势及其他非气候噪音是一种可行的思路。利用EMD运算程序，对祁连山孔岗木祁连圆柏[24]（Kgm02）02a和09a样芯宽度测量序列进行逐级分解（图1和图2）。可以看出两个树芯的宽度序列均包含了7个本征函数分量，每个IMF分量对应于一个窄波段的信号，反映树轮序列记录的该尺度信号的自然波动过程，IMF7均为原始变量的最低频趋势分量，各分量之和可还原为原始变量序列（蓝线）。

图1所示Kgm0202a芯的7条IMF分量反映了该样芯由高频至低频不同准周期尺度上的瞬时波动特征，而本研究中主要关注EMD分解所得趋势项对树轮宽度序列生长趋势的拟合效果，趋势项IMF7分量1692—2000年间由高至低，逐渐趋于平稳，反映的生长趋势近似于负指数函数拟合曲线。

图2展示的Kgm0209a样芯的IMF7趋势分量变化过程与前者有较大差异，该序列的IMF7分量1651—1757年逐步加速至峰值，1757年之后开始逐步减慢，至19世纪末趋于平稳，与树木自幼龄期到达壮年的自然生长过程极为吻合。该曲线近似于样条函数拟合曲线，但不同之处在于样条函数的步长选取较为主观，研究中如果选取不当则会导致过度拟合的情形，而EMD分解则是遵循序列本身的自然震荡进行逐级分解，不存在人为假定和预先设定函数形式的问题。

Kgm0202a和Kgm0209a虽然均是按照同一原理进行分解，可因为序列本身包含的震荡过程不同，也就是说两条序列受起始年靠近髓心位置不同及局地小生境的不同，得到的序列生长趋势项也表现了较大差异。因此，靠一个模式或一个一成不变的公式是很难满足完全剔除掉所有样芯的生物学效应的。在实际应用中根据所采集样芯的生态环境及树种，参照多序列的综合分析，才是保证较合理剔除各样芯生长趋势的关键。

进一步对Kgm02样点的65个样芯分别进行了EMD分解，各序列的最低频分量趋势项见图3。由图3可以看出，EMD分解能较好地适应树轮序列趋势项多样不一的情况，针对各序列不同的生境及生长情况逐步分离各序列所包含的多尺度震荡分量，同时避免了选择样条函数拟合步长选取的主观性，有助于较准确的获取趋势项。

图1 Kgm0202a宽度序列EMD分解各分量

图2 Kgm0209a宽度序列EMD分解各分量

图3 Kgm02各样芯宽度序列EMD分解的趋势项分量

趋势项包含了传统使用的线性函数、负指数函数、类样条函数等各种形式，但总体而言，趋势项在末端均表现出趋于平稳的特点，该特点与树木生长晚期生理学效应减弱并趋于平稳的特征较为吻合，进一步证实了EMD在趋势项获取中的可行性。2.3EMD与传统方法对圆柏生长量订正的比较

为了检验EMD方法在祁连圆柏生长量订正中的实际效果，遵循常规的树轮曲线标准化订正原理（原始曲线除以趋势项），分别按照传统拟合法和EMD方法对各序列进行趋势项分离，并对传统方法和EMD方法获得的65个轮宽序列去除生长趋势后得到的指数序列进行逐一对比。图4展示了四组具有代表性的序列细节对比图，上下两图分别为一组，共4组（a，b；c，d；e，f；g，h）。每组上图为轮宽曲线及传统线拟合（线性或负指数）及EMD分解的趋势项，下图为对应的去趋势后指数序列。

其中最为突出的差异之一表现为两种方法得到的指数序列在接近髓心端的差异。以祁连圆柏代表性样芯Kgm0211a（图4a和4b）为例可以看出，传统方法拟合的趋势项是斜率为负的直线（图4a红线），EMD分解得到的趋势项（图4a蓝线）自1600年开始先慢后快，至1850年左右趋于平稳。相对应两种方法得到的指数序列在髓心端（最后50a）差异较大，传统方法得到的序列值（图4b红线）高于EMD方法得到的指数值（图4b蓝线）大约0.2个标准差，导致髓心端异常上翘。

相对于髓心端的异常，两种方法对部分序列树皮端的订正也存在一定的差异。祁连圆柏代表性序列为Kgm0219b（图4c黑线），两种方法对原始宽度序列的拟合线（图4c红线和蓝线）较为相近，仅在末端直线拟合表现了持续性的下降，而EMD趋势项趋于水平值。相对应得到的指数序列前者在后5a出现了整个时段的异常极大值，而后者则较好地避免了因趋势线的误差导致的偏差。

结合以上两点不同，还有一类序列在两种方法的对比中，整个时段会表现为波动性偏差。以圆柏Kgm0222b（图4e黑线）为例可以看出，传统方法的拟合线仍为直线（图4e红线），而EMD得到的趋势项自序列始端（1787年）生长量的变化过程为由慢至快加速，1865年左右达到峰值，后逐渐减速至1990年趋于水平（图4e蓝线）。相对应的两条指数序列偏差可分为三段：1785—1830年（传统方法的指数序列偏低）、1830—1935年（传统方法的指数序列偏高）和1935—2000年（传统方法的指数序列偏低）。而这些偏差分析可知，正是由于传统拟合线对树木生物学效应低估的结果，EMD方法得到的指数序列则避免了这个偏差。

祁连圆柏Kgm0202a（图4g，4h）则代表了两种方法得到的较为一致的指数序列的例子。由图4g可以看出，传统方法拟合得到的趋势线与EMD方法分解得到的趋势项一致，相对应订正后的指数序列也几乎完全重合，高低频变化完全一致。对本文所选样本祁连圆柏的统计分析表明，当轮宽序列本身所包含生长趋势项较为单调，且生长过程中不存在异常小生境扰动时，两种方法得到的结果也会相对比较一致。

虽然整体上EMD对大部分样芯的订正都表现出了较大优势，但对祁连圆柏65个样芯的分析中也发现个别样芯由EMD分解得到的趋势项也会出现失效的情况。本文分析发现Kgm0215a和Kgm0218c两个样芯的EMD分解趋势项在生长量订正中需要特殊处理。为了更清楚地看出EMD方法与传统的负指数和线性函数拟合方法对多复本量树轮曲线订正的差异，本文还对比了两种方法最终建立的Kgm02研究点（65芯）的树轮年表，结果如图5。

由图5可以看出，两条序列整体所反映的高频信息极为一致，虽然由上面的分析可知少部分样芯在传统方法处理中会出现树皮端的变形，但所有序列平均之后树皮端变形被削弱，故两种方法得到的该样点年表在树皮端基本一致；然而，因为传统方法对大部分样芯接近髓心端的部分生物学效应剔除效果较差，因此两种方法得到的所有序列平均指数在髓心端的100a内出现了较大差异。EMD方法的年表在1430—1500年略高于传统方法得到的指数值，可能反映了该时期树木生长环境较为适宜的气候条件，然而，具体的气候指示意义需进一步分析。

综上所述，EMD方法与传统方法对祁连圆柏生长量订正的改进主要体现在序列的前后两端，EMD方法分解得到的趋势项与树木生长自幼龄期开始慢—快—慢的生长特征较为吻合，因此对于序列两端生长量的订正极为适应，效果远高于线性或负指数函数拟合。在实际应用中可重点分析去趋势序列前后两端方差的异常变化，进而评估树轮生长量的订正效果。然而，对于一种新方法的应用，仍然需要更大量样本实验的对比验证，尤其要根据不同的树种及不同的采样环境进行深入分析，确定最适宜的生长量订正方法和手段。

3 结论

图4 Kgm02各样芯宽度序列传统方法拟合（红线）和EMD分解趋势分量拟合（蓝线）对比

图5 Kgm02样点传统年表与EMD年表对比

EMD是对一个时间序列多尺度变化特征全面分析的新方法，可以将时间序列包含的不同尺度的自然振荡和趋势逐级分离，有利于更好地认识目标序列的多尺度变化特征和总体趋势。本文针对祁连圆柏宽度序列生长量的订正进行了传统方法和EMD分析手段的全面对比分析，初步探讨EMD分解趋势项在气候变化研究中对生长量订正的改进效果及可应用性。分析初步表明，与传统方法相比，EMD方法在生长量订正的应用在序列髓心端和树皮端表现出最为明显的提高。传统线性或负指数的拟合在树皮端会引起少数样芯序列5～10a的异常极大值；而在髓心端对大部分样芯先慢后快再慢至平稳的生长趋势的反映较差，因此在髓心端大部分样芯生长量订正失效，出现接近髓心端100a内异常上翘或降低的现象，因此在研究中树轮年表的后50～100a常常出现较大的方差而被遗弃，EMD方法在一定程度上缓解了这一问题。

总体而言，EMD方法在祁连圆柏宽度序列生长量订正中的适用性得到了初步证实，为生长量的订正提供了一种新的可选的方法。与传统方法相结合，对不同样芯灵活选取最合适的拟合方法将会得到更加准确的气候代用序列。该方法在改进树轮宽度生长量订正及标准化的应用中有较大优势和潜力。然而，针对不同的树种、不同的研究环境，对该方法的应用仍需谨慎，需要更多的实验进行验证和分析，尤其是对每一个样芯的情况均需要仔细研究，综合分析，才可能得到真正能反映气候信息的可靠序列，得到对气候变化的合理解释和分析结果。

[1]ICSU.The Initial core projects[M].IGBP“Global Change”Report No.12,Berne,Switzerland.1990.330.

[2]张永香.巴丹吉林南缘近两百年来年降水重建及初步分析[J].沙漠与绿洲气象，2015，9（1）：12-16.

[3]尚华明，魏文寿，袁玉江，等.帕米尔东北部昆仑圆柏850年树轮宽度年表的建立及其气候意义[J].沙漠与绿洲气象，2015，9（1）：6-11.

[4]Zhu h F,Zheng Y h,Shao X m,et al.Millennial

temperature reconstruction based on tree-ring widths of Qilian juniper from Wulan,Qinghai Province,China[J]. Chinese Science Bulletin，2008，53：3914-3920.

[5]袁玉江，喻树龙，穆桂金，等.天山北坡玛纳斯河355a来年径流量的重建与分析[J].冰川冻土，2005，27：411-417.

[6]Liang E,Liu X,Yuan Y,et al.The 1920s drought recorded by tree rings and historical documents in the semi-arid areas of Northern China[J].Climatic Change，2006，79：403-432.

[7]Yang B,Qin C,Wang J,Hem,et al.A 3，500-year treering record of annual precipitation on the northeastern Tibetan Plateau[J].PNAS，2014，111：2903-2908.

[8]Liu Y,Sun J,Song H,et al.Tree-ring hydrologic reconstructions for the Heihe River watershed,western China since AD 1430[J].Water Research，2010，44：2781-2792.

[9]Shao X,Huang L,Liu H,et al.Reconstruction of precipitation variation from tree rings in recent 1000 years in Delingha,Qinghai[J].Science in China Seriesd：Earth Sciences，2005，48：939-949.

[10]Zhang,Q B,Michael N E,Lyu L X.Moisture dipole over the Tibetan Plateau during the past five and a half centuries[J].Nature Communication.2015，6：8062，doi：10.1038/ncomms9062.

[11]Cook E R,Kairiukstis La.Methods of Dendrochronology[M].Netherlands:Kluwer Academic，1990

[12]Fritts h C.Tree Ringsand Climate[M].Academic Press, London，1976.

[13]吴祥定.树木年轮与气候变化[M].气象出版社，1990.

[14]Phipps R L.Annual growth of suppressed chestnut oak and red maple,a basis for hydrologic inference[M],U.S.geological Survey Professional Paper 485-c，Washington D.C.1967.

[15]Cook E R,Briffa K R,Shiyatov S，et al.Tree-ring standardization and growth-trend estimation[M].In:Cook ER,Kairiukstis L A（eds.）.Methods of Dendrochronology.Kluwer Academic Publishers，1990，104-123.

[16]Cook E R.A time series analysis approach to tree-ring standardization[M].Disscrtation for the Doctoral Degree. Tucson:University of Arizona，1985，1-171.

[17]Mann M E,Bradley R S,Hughes M K.North hemisphere temperature during the last millennium:inferences, uncertainties and limitations[J].Geophysical Research Letters，1999，26：759-762.

[18]Briffa K R,Jones Pd,Schweingruber F H.Tree-ring reconstructions of summer temperature patterns across western North America since 1600[J].Journal of Climate，1992，5：735-754.

[19]Briffa K R,Jones Pd,Schweingruber F H，et al.Tree-ring variables as proxy-climate indicators:Problems with low frequency signals[J].Springer-Verlag，Berlin.、1996.

[20]Esper J,Cook E R,Schweingruber F H.Low-frequency signals in long tree-ring chronologies for reconstructing past temperature variability[J].Science，2002，295：2250-2253.

[21]Melvin，T M and Briffa,K R.A"Signal-Free"approach to Dendroclimatic Standardisation[J].Dendrochronologia 26，2008，71-86.doi:10.1016/j.dendro.2007.12.001.

[22]徐岩，邵雪梅.柴达木盆地东缘祁连圆柏轮宽序列标准化的方法研究[J].地理学报，2006，9：919-928.

[23]Huang N E,Shen Z,Long R S，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary series analysis[J].Proc.R. Soc.Lond.A.1998，454：899-995.

[24]田沁花，周秀骥，勾晓华，等.祁连山中部近500年来降水重建序列分析[J].中国科学：地球科学，2012，42（4）：536-544.

Study of Standardization of the Qilian Juniper Ring-width Series Based on Empirical Mode Decomposition Method

ZHANG Yong
（Key Laboratory of Land Surface Pattern and Simulation,Institute of geographic Sciences and Natural Resources Research,Chinese academy of Sciences，Beijing 100101，China）

Standardization of tree-ring measurements is one of the most important procedures in dendroclimatology.In this paper,applying Empirical Mode decomposition（EMD）method to the Qilian Juniper（Sabina przewalskii kom）ring-width samples,the potential and limitations of the EMD method are investigated.Tree-ring series are firstly decomposed by EMD to produce a number of Intrinsic mode Functions（IMFs）in different frequency band with different physical meanings. The lowest-frequency signals mainly reflect the age-dependent growth trend,which are employed to correct the raw ring-width series to get the standardization chronology.Both the chronologies from the traditional and EMD method preserve the same high-frequency,but the different low-frequency variations in the 50-100yrs near to the pith.In this period,the variations of tree-ring index based on traditional method are large in general,while those based on EMD method are more stable.EMD method is a suitable tool for preserving more reliable climate information than traditional method, especially for the pith periods and young periods of tree-ring series,which can be applied to fit the growth trend of Qilian Juniper.

tree ring index;standardized method;Empirical Mode Decomposition（EMD）

P468

A

1002-0799（2015）05-0009-07

张永.经验模态分解方法对祁连圆柏年轮宽度序列标准化初探[J].沙漠与绿洲气象，2015，9（5）：9-15.

10.3969/j.issn.1002-0799.2015.05.002

2015-04-02；

2015-06-30

国家自然科学面上基金（41471087和41175066）资助。

张永（1979-），男，助理研究员，主要从事古气候重建与气候变化研究。E-mail：zhangyong@igsnrr.ac.cn