ARSTAN程序和R語言dplR擴(kuò)展包進(jìn)行樹輪年表分析的比較研究

趙守棟，江 源，焦 亮，王明昌，張凌楠，李文卿

北京師范大學(xué)資源學(xué)院， 北京 100875

ARSTAN程序和R語言dplR擴(kuò)展包進(jìn)行樹輪年表分析的比較研究

趙守棟，江 源*，焦 亮，王明昌，張凌楠，李文卿

北京師范大學(xué)資源學(xué)院， 北京 100875

在樹輪年代學(xué)領(lǐng)域，ARSTAN是去趨勢處理和建立年表方面應(yīng)用最為廣泛的程序，而新興的R語言dplR擴(kuò)展包實(shí)現(xiàn)了ARSTAN的主要功能，且具有源代碼公開、擴(kuò)展性強(qiáng)等優(yōu)點(diǎn)，是傳統(tǒng)程序的良好補(bǔ)充。使用賀蘭山青海云杉(PiceaCrassifolia)樹輪寬度數(shù)據(jù)，分析了ARSTAN和dplR進(jìn)行樹輪年代學(xué)分析所得結(jié)果的差異。結(jié)果顯示，兩種程序計(jì)算平均敏感度和一階自相關(guān)系數(shù)的平均誤差為0.005—0.008，但具有確定的轉(zhuǎn)換關(guān)系；兩種程序如果使用同種方法去趨勢，擬合曲線的參數(shù)相近，建立標(biāo)準(zhǔn)年表的平均誤差為0.002；擬合自回歸模型時(shí)差異較大，其中時(shí)域上表現(xiàn)為差值年表起始30a內(nèi)差異顯著，在頻域上表現(xiàn)為dplR的差值年表保留的低頻信息較少；年表統(tǒng)計(jì)量計(jì)算和公共區(qū)間分析中，不同程序計(jì)算樣本總體代表性和信噪比的差異較大。分析表明，兩程序在擬合生長趨勢和自回歸模型時(shí)存在算法上的較大差異，同時(shí)年表統(tǒng)計(jì)量和公共區(qū)間各指標(biāo)的算法也不盡一致，但存在較為確定的轉(zhuǎn)換關(guān)系。對開展不同來源數(shù)據(jù)的整合分析提出了建議，應(yīng)明確不同研究中樹輪數(shù)據(jù)的處理過程，在條件允許時(shí)使用同一程序或算法重新處理數(shù)據(jù)，確保結(jié)果的可比性。

R語言；dplR；去趨勢；公共區(qū)間分析；樹輪年代學(xué)

樹輪年代學(xué)研究樹木年輪的寬度、密度等性質(zhì)構(gòu)成的時(shí)間序列及其與環(huán)境因子的響應(yīng)關(guān)系[1]，在反演過去和預(yù)測未來氣候變化中扮演著重要作用[2- 5]。樹輪年代學(xué)中，去除樹木自身的生長趨勢并建立樹輪年表是開展研究的基礎(chǔ)步驟之一。在去趨勢處理和建立年表方面，Cook等編寫的ARSTAN程序是目前最權(quán)威的分析工具，在樹輪年代學(xué)的發(fā)展過程中發(fā)揮了巨大作用。然而ARSTAN基于Fortran語言編寫，同時(shí)缺乏系統(tǒng)而詳細(xì)的說明文檔，使用者不容易理清其內(nèi)部的算法，實(shí)現(xiàn)新的分析方法也比較困難，這在一定程度上限制了研究的深入和學(xué)科間的交流。

R語言是當(dāng)前發(fā)展較快的統(tǒng)計(jì)分析工具之一，在眾多研究領(lǐng)域得到了越來越廣泛的應(yīng)用。Bunn等使用R語言編寫了dplR擴(kuò)展包，實(shí)現(xiàn)了樹輪數(shù)據(jù)讀取、轉(zhuǎn)換、分析、繪圖等功能[6]。相比于經(jīng)典工具而言，R語言和dplR等相關(guān)擴(kuò)展包的源代碼完全公開，且容易理解，研究者可以根據(jù)自身需要修改或編寫程序；同時(shí)使用者眾多，網(wǎng)絡(luò)社區(qū)發(fā)達(dá)，有利于不同領(lǐng)域的相互交流。因此，R語言及dplR擴(kuò)展包是ARSTAN等傳統(tǒng)程序的良好補(bǔ)充，對于樹輪年代學(xué)的進(jìn)一步發(fā)展具有重要意義。

目前在國外研究中，dplR被廣泛應(yīng)用于去趨勢處理和建立年表[7- 9]，在交叉定年[10]、小波分析[3]等方面也發(fā)揮了作用；在其基礎(chǔ)上又發(fā)展出交互式去趨勢[11]、響應(yīng)函數(shù)分析[12]、擾動(dòng)事件重建[13]等相關(guān)擴(kuò)展包，逐漸形成了一套較為完整的樹輪數(shù)據(jù)分析軟件體系。而在國內(nèi)樹輪年代學(xué)領(lǐng)域，使用dplR開展的研究還較少，使用的功能也相對單一[14]。作為一個(gè)新興的樹輪年表分析工具，dplR處理結(jié)果與傳統(tǒng)程序一致與否，將直接影響到不同來源研究數(shù)據(jù)的可比性，因此有必要對ARSTAN和dplR在算法和結(jié)果上進(jìn)行比較分析。

本文使用賀蘭山青海云杉樹輪寬度數(shù)據(jù)，分別借助ARSTAN和dplR進(jìn)行原始序列統(tǒng)計(jì)量計(jì)算、去趨勢處理、年表建立、公共區(qū)間分析等操作，通過相關(guān)分析、小波分析等手段比較了兩者所得結(jié)果的差異，并從算法的角度探討了差異的來源，為開展樹輪年代學(xué)數(shù)據(jù)的整合分析提出了需要注意的問題及相應(yīng)建議。

1 數(shù)據(jù)與方法

1.1 樹輪寬度資料與氣候資料

本文中使用的樹輪資料為賀蘭山青海云杉(PiceaCrassifolia)樹輪寬度數(shù)據(jù)[15]，采樣點(diǎn)位于賀蘭山東坡林線附近，共包含來自25棵樣樹50根樣芯的樹輪寬度測量結(jié)果，數(shù)據(jù)精度為0.01 mm。使用COFECHA程序檢驗(yàn)交叉定年質(zhì)量，結(jié)果顯示序列平均長度為83.9 a，平均敏感度為0.393，序列間相關(guān)系數(shù)為0.825。

1.2 ARSTAN處理

ARSTAN程序目前分為Windows XP版和OS X版，早期的DOS版已不再更新。本文主要使用Windows XP版ARSTAN進(jìn)行去趨勢處理并建立年表，程序下載自美國哥倫比亞大學(xué)樹木年輪實(shí)驗(yàn)室(www.ldeo.columbia.edu/tree-ring-laboratory)，版本號為44。選擇的去趨勢方法為修正的負(fù)指數(shù)函數(shù)，即：

f(i)=a×exp(-b×i)+d

(1)

式中,f(i)為生長趨勢估計(jì)值，i為年份的序數(shù)。設(shè)定約束條件為d>0，擬合失敗則擬合斜率非正的線性函數(shù)。公共區(qū)間設(shè)定為所有樹芯均含有的共同時(shí)間段，即1966—2009年。ARSTAN最終建立標(biāo)準(zhǔn)年表(standard chronology, STD)、差值年表(residual chronology, RES)和自回歸年表(autoregressive standard chronology, ARS)，并給出年表的基本統(tǒng)計(jì)量和公共區(qū)間分析結(jié)果。

1.3 dplR處理

R語言主程序和dplR擴(kuò)展包均下載自綜合R檔案網(wǎng)絡(luò)(cran.r-project.org)，其中R語言版本號為3.1.1，dplR版本號為1.6.1。使用“rwl.stats”函數(shù)計(jì)算原始樹輪寬度序列和樹輪寬度指數(shù)年表的基本統(tǒng)計(jì)量。使用“detrend.series”函數(shù)對原始樹輪寬度序列逐一進(jìn)行去趨勢操作，確保各樣芯使用的去趨勢方法與ARSTAN中相同，使得結(jié)果具有可比性。使用“rwi.stats”函數(shù)對樹輪寬度指數(shù)序列進(jìn)行公共區(qū)間分析，時(shí)間段同樣設(shè)定為1966—2009年，這也是程序默認(rèn)的公共區(qū)間。使用“chron”函數(shù)采用雙權(quán)重平均法建立年表，最終得到STD和RES年表。

1.4 ARSTAN與dplR結(jié)果對比

比較ARSTAN和dplR所得原始輪寬序列統(tǒng)計(jì)量的差異，使用Spearman相關(guān)系數(shù)考察誤差與序列長度的關(guān)系。對比ARSTAN和dplR擬合生長趨勢所得曲線的參數(shù)及擬合結(jié)果，使用Spearman相關(guān)系數(shù)考察兩程序所得STD、RES年表的相關(guān)性，以ARSTAN所得結(jié)果為準(zhǔn)，計(jì)算dplR所得年表的標(biāo)準(zhǔn)誤；針對差異相對較大的RES年表，使用Morlet小波變換分析兩程序所得RES年表在頻域上的差異。最后比較不同程序在年表統(tǒng)計(jì)量和公共區(qū)間分析結(jié)果上的差異。相關(guān)分析中，顯著性水平設(shè)置為0.05。

2 結(jié)果

2.1 原始輪寬序列統(tǒng)計(jì)量

針對50個(gè)原始輪寬序列，ARSTAN和dplR計(jì)算所得均值和標(biāo)準(zhǔn)差完全相等；而dplR所得平均敏感度(mean sensitivity, MS)和一階自相關(guān)系數(shù)(first-order autocorrelation, AC)均略小于ARSTAN所得結(jié)果，其中MS的平均誤差為0.005，標(biāo)準(zhǔn)誤為0.0002，AC的平均誤差為0.008，標(biāo)準(zhǔn)誤為0.0003，兩者的誤差與序列長度存在顯著負(fù)相關(guān)關(guān)系(P<0.05)。進(jìn)一步計(jì)算發(fā)現(xiàn)，ARSTAN和dplR所得MS和AC滿足：

(2)

(2)

式中,MSr為dplR所得MS，MSa為ARSTAN所得MS，ACr為dplR所得AC，ACa為ARSTAN所得AC，n為序列長度。隨著序列長度的增加，ARSTAN與dplR所得結(jié)果的差異將逐漸減小，當(dāng)序列長度在100年以上時(shí)，相對誤差可以降至1%以下。

2.2 建立樹輪寬度指數(shù)年表

ARSTAN和dplR擬合所得負(fù)指數(shù)函數(shù)的參數(shù)a、b完全相同，參數(shù)d的誤差在0.0003以內(nèi)，而線性函數(shù)參數(shù)完全相同。去趨勢所得各樣芯的兩個(gè)標(biāo)準(zhǔn)化序列間存在差異的年份不超過10年，且誤差均在0.003以內(nèi)。在此基礎(chǔ)上，兩程序建立的STD年表相關(guān)系數(shù)為0.999(P<0.05)，對應(yīng)年份寬度指數(shù)誤差在0.05以內(nèi)，平均誤差為0.005，標(biāo)準(zhǔn)誤為0.0010，相位變化情況基本相同(圖1)。

在建立差值年表時(shí)，ARSTAN對各標(biāo)準(zhǔn)化序列均擬合為一階自回歸模型，所得差值序列與原序列長度相等；而dplR則擬合成階數(shù)各異的自回歸模型，包括29個(gè)一階模型和17個(gè)二階或更高階模型，序列初始損失的年份數(shù)與模型階數(shù)相等，另有4個(gè)樣芯未擬合成自回歸模型。各樣芯的兩個(gè)差值序列間的平均誤差為0.034，標(biāo)準(zhǔn)誤為0.0011。兩程序建立的RES年表間相關(guān)系數(shù)為0.987(P<0.05)，在1914年前(即序列前30年)，兩者平均誤差為0.048，標(biāo)準(zhǔn)誤為0.0096；而1914年后，兩者平均誤差為0.028，標(biāo)準(zhǔn)誤為0.0021。兩個(gè)RES年表僅在起始的10年內(nèi)存在相位不一致的情況，例如1887—1889年(圖1)。

對兩個(gè)RES年表分別進(jìn)行Morlet小波變換(圖2)。結(jié)果顯示，兩個(gè)RES年表含有的高頻信息基本相同，均表現(xiàn)出1940年前2—11a尺度的震蕩能量較小，周期變化不明顯，在1940年后2—11a尺度的能量增強(qiáng)，周期變化顯著；而在低頻信息方面，ARSTAN建立的RES年表在1940年后穩(wěn)定存在著32a左右的周期，而dplR建立的RES年表32a左右的周期震蕩能量始終較小。

圖1 ARSTAN和dplR所得樹輪寬度指數(shù)年表Fig.1 Tree-ring chronologies produced by ARSTAN and dplR(a)為標(biāo)準(zhǔn)年表；(b)為差值年表；(c)為不同程序建立的年表之差

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

2.3 年表統(tǒng)計(jì)量和公共區(qū)間分析結(jié)果

年表統(tǒng)計(jì)量及公共區(qū)間分析結(jié)果見表1。在年表統(tǒng)計(jì)量方面，ARSTAN的STD、RES年表和dplR的STD年表的序列起止年相同，而dplR的RES年表與這3個(gè)年表相比少1年(1884年)。兩個(gè)STD年表的均值、標(biāo)準(zhǔn)差均相差0.001；MS相差0.003，相對誤差約為1%；AC相差0.005，相對誤差約為1.5%。而兩個(gè)RES年各項(xiàng)統(tǒng)計(jì)量的誤差均大于STD年表，其中dplR的RES年表MS比ARSTAN低4%，而AC的絕對值低36%。

表1 年表統(tǒng)計(jì)量及公共區(qū)間分析結(jié)果Table 1 Chronologies statistics and common interval analysis

3 討論

3.1 ARSTAN與dplR計(jì)算序列統(tǒng)計(jì)量的算法差異

根據(jù)原始輪寬序列統(tǒng)計(jì)量比較結(jié)果，ARSTAN與dplR在MS和AC的計(jì)算上存在著較為確定的換算關(guān)系，這很可能是由于兩者使用了不同的計(jì)算公式。dplR計(jì)算MS的公式為[16]：

(4)

式中,n為序列長度，wi為第i年的樹輪寬度。根據(jù)結(jié)果推測，ARSTAN在計(jì)算MS時(shí)，將公式中的n-1換成了n-2。

dplR調(diào)用R中的“acf”函數(shù)計(jì)算AC，其計(jì)算公式可表達(dá)為[1,17]：

(5)

式中,s為序列的標(biāo)準(zhǔn)差。ARSTAN在計(jì)算AC時(shí)，可能將s處理成了總體標(biāo)準(zhǔn)差：

(6)

盡管計(jì)算方法不同，由于兩指標(biāo)相對n的取值具有單調(diào)性，因此同一個(gè)程序的計(jì)算結(jié)果具有可比性，而不同程序的計(jì)算結(jié)果需要經(jīng)過換算后才可相互比較。

由于不帶蓄電池備電，所以在沒有市電的情況下，這種供電方式不能為信號源提供電能。綜上所述，直流輸入不帶備電的方案比較適合于信號源為直流輸入，對服務(wù)要求不高的場景。

3.2 ARSTAN與dplR建立年表的算法差異

ARSTAN和dplR均具有多種去趨勢方法，兩者的對應(yīng)關(guān)系如表2所示，其中在擬合負(fù)指數(shù)函數(shù)的程序?qū)崿F(xiàn)方面，兩者具有較大差異。在算法上，ARSTAN使用Fritts等專為負(fù)指數(shù)函數(shù)擬合問題設(shè)計(jì)的算法[18]，該方法充分考慮到負(fù)指數(shù)函數(shù)的曲線特征，首先估計(jì)曲線的切線斜率b，通過調(diào)節(jié)b的值進(jìn)行迭代。而dplR首先估計(jì)y軸截距a和漸近線d的初始值，b的初始值則根據(jù)經(jīng)驗(yàn)直接設(shè)為0.01，再調(diào)用R內(nèi)置的“nls”函數(shù)進(jìn)行曲線參數(shù)估計(jì)。如果不設(shè)定約束條件，當(dāng)樹輪寬度序列的趨勢較符合負(fù)指數(shù)函數(shù)曲線時(shí)，兩種算法均能得到精確的擬合結(jié)果；而當(dāng)樹輪寬度序列的趨勢近似于直線時(shí)，d為負(fù)值，這時(shí)dplR的初始值設(shè)定嚴(yán)重偏離最優(yōu)解，導(dǎo)致迭代過程報(bào)錯(cuò)，轉(zhuǎn)而使用線性函數(shù)進(jìn)行擬合。如果設(shè)定約束條件d>0，ARSTAN首先在整個(gè)定義域范圍內(nèi)尋找各參數(shù)的最佳估計(jì)值，然后將參數(shù)帶入約束條件中進(jìn)行判斷，如果不符合條件則拒絕擬合負(fù)指數(shù)函數(shù)，轉(zhuǎn)而擬合線性函數(shù)；而dplR則是在“nls”函數(shù)中縮小參數(shù)的取值范圍，最終擬合得到參數(shù)的局部最優(yōu)解。當(dāng)擬合方法出現(xiàn)分歧時(shí)，去趨勢所得序列將出現(xiàn)較大差異，而本文在操作中參考ARSTAN的運(yùn)行結(jié)果，使用dplR對每根樣芯單獨(dú)去趨勢，確保兩者使用的去趨勢方法相同。

盡管擬合生長趨勢的方法相同，ARSTAN和dplR建立的標(biāo)準(zhǔn)年表仍存在一定差異，這與缺失值的處理方式有關(guān)。為防止運(yùn)算中可能出現(xiàn)的錯(cuò)誤，dplR將序列中的缺失值(一般賦值為0)均替換成0.001，因此dplR所得標(biāo)準(zhǔn)化序列不存在樹輪寬度指數(shù)為0的情況；而ARSTAN并未對缺失值進(jìn)行特殊處理，因此標(biāo)準(zhǔn)化序列中仍存在0值。

表2 ARSTAN和dplR去趨勢方法的對應(yīng)關(guān)系Table 2 Corresponding relationships between detrending methods of ARSTAN and dplR

建立差值年表，需要去除樹輪寬度指數(shù)標(biāo)準(zhǔn)化序列中的自回歸成分。ARSTAN首先用所有標(biāo)準(zhǔn)化序列擬合出一個(gè)p階合并自回歸模型，將擬合值作為樹木群體所共有的持續(xù)性生長量。再對各標(biāo)準(zhǔn)化序列分別擬合限定為p階的自回歸模型，將擬合值作為整體生長模式下個(gè)體的持續(xù)性生長量，而將殘差作為差值序列，并使用雙權(quán)重平均法建立RES年表。如果將群體的持續(xù)性生長量加回到RES年表上，則得到ARS年表[19]。dplR未考慮整體的生長模式，而是直接對各標(biāo)準(zhǔn)化序列分別擬合自回歸模型，不同序列的模型可以具有不同的階數(shù)，再將殘差作為差值序列，建立RES年表。針對本文所用數(shù)據(jù)，dplR的計(jì)算過程對個(gè)體特有的持續(xù)性生長量剔除得更為徹底，低頻信息保留得更少，這種差異在頻域上更為明顯(圖2)。

擬合p階自回歸模型會(huì)造成序列前p年樹輪寬度指數(shù)的損失。針對這一情況，ARSTAN使用序列均值將標(biāo)準(zhǔn)化序列向前延長p年[19]，從而使自回歸模型計(jì)算出的擬合值序列與原序列等長，避免了序列長度上的損失。而dplR未對初始年份的損失進(jìn)行特殊處理，因此兩個(gè)RES年表存在的差異相對STD年表更大，尤其是樣本量較小的最初若干年。根據(jù)本文結(jié)果，起始值處理方式造成的差異會(huì)在約30a后消失，同時(shí)考慮到參與“樹輪-氣候”關(guān)系分析的年表區(qū)間通常遠(yuǎn)小于年表總長度，如果起始年份不參與后續(xù)的分析，則起始年份的不同處理方式對最終結(jié)果的影響較小[19- 20]。

3.3 ARSTAN與dplR進(jìn)行公共區(qū)間分析的算法差異

公共區(qū)間分析是指選擇包含特定時(shí)間范圍的樹輪寬度序列，通過計(jì)算年表信號，評價(jià)最終所得年表對森林總體生長的代表性，其中年表信號是衡量序列中包含的共同變化的統(tǒng)計(jì)量[21]。DOS版ARSTAN在導(dǎo)入測量數(shù)據(jù)時(shí)會(huì)提供一個(gè)推薦的公共區(qū)間范圍，其選擇標(biāo)準(zhǔn)是公共區(qū)間內(nèi)包含盡可能多的測量值，即滿足區(qū)間內(nèi)完整樣芯數(shù)與區(qū)間長度乘積取得最大值；而Windows XP版ARSTAN不再提供這一結(jié)果。dplR中對數(shù)據(jù)集進(jìn)行去趨勢使用“detrend”函數(shù)，其默認(rèn)公共區(qū)間是所有樣芯共有的公共區(qū)間，該區(qū)間長度一般小于DOS版ARSTAN提供的公共區(qū)間長度。dplR中的“common.interval”函數(shù)可以給出各種標(biāo)準(zhǔn)對應(yīng)的時(shí)間范圍。

在進(jìn)行公共區(qū)間分析時(shí)，最初使用方差分析對年表信號進(jìn)行估計(jì)[1]，后來的研究發(fā)現(xiàn)使用相關(guān)矩陣也可以估計(jì)年表信號的大小[22]。以序列間的平均相關(guān)系數(shù)(mean series intercorrelation, Rbar)作為年表信號的衡量指標(biāo)，可以通過計(jì)算EPS和SNR來評價(jià)年表的質(zhì)量[17]：

(7)

(8)

式中,t為樣樹數(shù)量。對于某樹輪年表，如果每棵樣樹均取一根樣芯，則Rbar的取值為Rbt；如果每棵樣樹取多根樣芯，則Rbar的取值應(yīng)當(dāng)為Rbt與Rwt加權(quán)平均結(jié)果，即Reff，如果仍以Rbt作為年表信號則會(huì)低估樣本中的公共信息[21]。

根據(jù)運(yùn)算結(jié)果推測，ARSTAN可能使用樣芯數(shù)量m和Rtot計(jì)算EPS和SNR。這種算法將所有樣芯等同看待，沒有剔除同一樣樹的不同樣芯包含著的重復(fù)信息，因此可能高估了年表質(zhì)量。EPS除了用于評價(jià)年表整體質(zhì)量外，還用于計(jì)算子樣本信號強(qiáng)度(subsample signal strength, SSS)[21]，而EPS和SSS均是衡量年表可靠長度的重要指標(biāo)。由于ARSTAN計(jì)算得出的EPS偏高，可能會(huì)令EPS或SSS高于閾值(一般設(shè)為0.8或0.85)的最小樣本量偏低，造成對年表可靠長度的高估。

目前的樹輪年代學(xué)研究中，一棵樣樹往往取兩個(gè)甚至多個(gè)樣芯，因此dplR采用Reff計(jì)算EPS和SNR，計(jì)算結(jié)果更符合現(xiàn)有理論。使用基于Reff的公式重新評估本文中ARSTAN建立的年表，則STD年表EPS為0.986，SNR為72.276，而RES年表EPS為0.989，SNR為90.741，這些結(jié)果均與dplR所得結(jié)果接近。

3.4 對比不同來源結(jié)果時(shí)應(yīng)當(dāng)注意的問題

隨著樹木年輪基礎(chǔ)數(shù)據(jù)日漸豐富，對前期樣點(diǎn)尺度的研究成果進(jìn)行整合，開展區(qū)域、大陸甚至全球尺度上森林對于全球氣候變化的響應(yīng)分析已成為研究熱點(diǎn)之一[23- 24]。根據(jù)本文所得結(jié)果，在樹木年輪研究中如果使用不同程序處理原始數(shù)據(jù)，所得年表及其質(zhì)量評估需要考慮算法差異可能造成的影響。在標(biāo)準(zhǔn)年表建立方面，即使ARSTAN與dplR擬合生長趨勢的參數(shù)設(shè)定相同，擬合結(jié)果可能差異很大。而對于差值年表，在頻域上，ARSTAN建立的RES年表比dplR具有更多的低頻信息；在時(shí)域上，起始年份處理方式的不同，造成年表在最初的20—30年內(nèi)存在較明顯的差異。不同程序在MS、AC等指標(biāo)上可能存在計(jì)算方法差異；ARSTAN使用Rtot計(jì)算EPS和SNR等指標(biāo)，相對于文獻(xiàn)中基于Reff的算法而言，可能高估了年表中包含的共同信息[25- 26]。

綜合以上討論，對開展不同來源樹輪數(shù)據(jù)的整合分析提出如下建議：

(1)如果能夠在國際樹輪數(shù)據(jù)庫等共享平臺上獲取原始數(shù)據(jù)，建議使用同一程序完成年表建立和公共區(qū)間分析等處理步驟，并注意檢查各樣芯去趨勢方法的具體參數(shù)。

(2)如果僅能獲得文獻(xiàn)資料，建議記錄其處理過程的詳細(xì)信息，包括軟件的版本信息、去趨勢的具體方法等；當(dāng)信息足夠時(shí)，建議根據(jù)同一標(biāo)準(zhǔn)設(shè)定公共區(qū)間，并重新計(jì)算MS、EPS等指標(biāo)，確保結(jié)果的可比性。

4 結(jié)論

在樹輪年代學(xué)研究中，新興的R語言dplR擴(kuò)展包是對ARSTAN等傳統(tǒng)程序的良好補(bǔ)充，其平臺開放性可以促進(jìn)研究者對程序算法的認(rèn)識和跨學(xué)科的相互交流。本文對比了ARSTAN程序和dplR擴(kuò)展包，分別使用兩種程序?qū)R蘭山青海云杉樹輪寬度數(shù)據(jù)進(jìn)行原始序列統(tǒng)計(jì)量估計(jì)，完成了去趨勢處理并建立了標(biāo)準(zhǔn)年表和差值年表，通過公共區(qū)間分析評價(jià)了年表的質(zhì)量。分析表明，ARSTAN與dplR在擬合生長趨勢和自回歸模型時(shí)存在算法上的較大差異，同時(shí)平均敏感度、一階自相關(guān)系數(shù)、樣本總體代表性和信噪比等指標(biāo)的算法也不一致，但存在較為確定的轉(zhuǎn)換關(guān)系。整合分析不同來源的樹輪資料時(shí)，應(yīng)明確其數(shù)據(jù)處理過程，在條件允許的情況下使用同一程序或算法進(jìn)行去趨勢處理和公共區(qū)間分析。

[1] Fritts H C. Tree Rings and Climate. London: Academic Press, 1976:10-11,254-260.

[2] 彭劍峰, 勾曉華, 陳發(fā)虎, 劉普幸, 張永, 方克艷. 阿尼瑪卿山地不同海拔青海云杉(Piceacrassifolia)樹輪生長特性及其對氣候的響應(yīng). 生態(tài)學(xué)報(bào), 2007, 27(8): 3268- 3276.

[3] Pederson G T, Gray S T, Woodhouse C A, Betancourt J L, Fagre D B, Littell J S, Watson E, Luckman B H, Graumlich L J. The unusual nature of recent snowpack declines in the North American Cordillera. Science, 2011, 333(6040): 332- 335.

[4] Liu H Y, Park Williams A, Allen C D, Guo D L, Wu X C, Anenkhonov O A, Liang E Y, Sandanov D V, Yin Y, Qi Z H, Badmaeva N K. Rapid warming accelerates tree growth decline in semi-arid forests of Inner Asia. Global Change Biology, 2013, 19(8): 2500- 2510.

[5] 蘆曉明, 梁爾源. 灌木年輪學(xué)研究進(jìn)展. 生態(tài)學(xué)報(bào), 2013, 33(5): 1367- 1374.

[6] Bunn A G. A dendrochronology program library in R (dplR). Dendrochronologia, 2008, 26(2): 115- 124.

[7] Bader M K F, Leuzinger S, Keel S G, Siegwolf R T W, Hagedorn F, Schleppi P, K?rner C. Central European hardwood trees in a high-CO2future: synthesis of an 8-year forest canopy CO2enrichment project. Journal of Ecology, 2013, 101(6): 1509- 1519.

[8] Boden S, Kahle H P, Wilpert K V, Spiecker H. Resilience of Norway spruce (Piceaabies(L.) Karst) growth to changing climatic conditions in Southwest Germany. Forest Ecology and Management, 2014, 315: 12- 21.

[9] Rodríguez-González P M, Campelo F, Albuquerque A, Rivaes R, Ferreira T, Pereira J S. Sensitivity of black alder (Alnusglutinosa[L.]Gaertn.) growth to hydrological changes in wetland forests at the rear edge of the species distribution. Plant Ecology, 2014, 215(2): 233- 245.

[10] Martin A R, Caspersen J P, Fuller M M, Jones T A, Thomas S C. Temporal dynamics and causes of postharvest mortality in a selection-managed tolerant hardwood forest. Forest Ecology and Management, 2014, 314: 183- 192.

[11] Campelo F, García-González I, Nabais C. detrendeR -A Graphical User Interface to process and visualize tree-ring data using R. Dendrochronologia, 2012, 30(1): 57- 60.

[12] Zang C, Biondi F. Dendroclimatic calibration in R: The bootRes package for response and correlation function analysis. Dendrochronologia, 2013, 31(1): 68- 74.

[13] Altman J, Fibich P, Dolezal J, Aakala T. TRADER: A package for tree ring analysis of disturbance events in R. Dendrochronologia, 2014, 32(2): 107- 112.

[14] 王蔚蔚, 張軍輝, 戴冠華, 王秀秀, 韓士杰, 張寒松, 王云. 利用樹木年輪寬度資料重建長白山地區(qū)過去240年秋季氣溫的變化. 生態(tài)學(xué)雜志, 2012, 31(4): 787- 793.

[15] 趙守棟, 何新, 張冰琦, 劉琦, 王輝, 江源. 賀蘭山東坡高山林線青海云杉(Piceacrassifolia)徑向生長對氣候因子的響應(yīng). 北京師范大學(xué)學(xué)報(bào): 自然科學(xué)版, 2013, 49(5): 501- 505.

[16] Biondi F, Qeadan F. Inequality in paleorecords. Ecology, 2008, 89(4): 1056- 1067.

[17] Cook E R, Pederson N. Uncertainty, emergence, and statistics in dendrochronology//Hughes M K, Swetnam T W, Diaz H F, eds.Dendroclimatology: Progress and Prospects. Netherlands: Springer, 2011: 77- 112.

[18] Fritts H C, Mosimann J E, Bottorff C P. A revised computer program for standardizing tree-ring series. Tree-Ring Bulletin, 1969, 29(1- 2): 15- 20.

[19] Cook E R. A Time Series Analysis Approach to Tree Ring Standardization[D]. Tucson, Arizona, USA: University of Arizona, 1985:63-73.

[20] Cook E R, Shiyatov S, Mazepa V. Estimation of the mean chronology // Cook E R, Kairiukstis L A. Methods of Dendrochronology: Applications in the Environmental Sciences. Dordrecht: Kluwer Academic Publishers, 1990: 123- 132.

[21] Briffa K R, Jones P D. Basic chronology statistics and assessment // Cook E R, Kairiukstis L A. Methods of Dendrochronology: Applications in the Environmental Sciences. Dordrecht: Kluwer Academic Publishers, 1990: 137- 152.

[22] Wigley T M, Briffa K R, Jones P D. On the average value of correlated time series, with applications in dendroclimatology and hydrometeorology. Journal of Climate and Applied Meteorology, 1984, 23(2): 201- 213.

[23] 鄭景云, 邵雪梅, 郝志新, 葛全勝. 過去2000年中國氣候變化研究. 地理研究, 2010, 29(9): 1561- 1570.

[24] PAGES 2k Consortium. Continental-scale temperature variability during the past two millennia: Supplementary Information. Nature Geoscience, 2013, 6(5): 339- 346.

[25] Yang B, He M H, Melvin T M, Zhao Y, Briffa K R. Climate control on tree growth at the upper and lower Treelines: A case study in the Qilian Mountains, Tibetan Plateau. PLoS ONE, 2013, 8(7): e69065.

[26] Ahmed M, Palmer J, Khan N, Wahab M, Fenwick P, Esper J, Cook E. The dendroclimatic potential of conifers from northern Pakistan. Dendrochronologia, 2011, 29(2): 77- 88.

A comparative analysis of ARSTAN and the dplr package of R language in analyses of tree-ring chronologies

ZHAO Shoudong, JIANG Yuan*, JIAO Liang, WANG Mingchang, ZHANG Lingnan, LI Wenqing

CollegeofResourcesScienceandTechnology,BeijingNormalUniversity,Beijing100875,China

Dendrochronology plays an important role in estimating past climatic conditions and predicting future climate change. Detrending and chronology development are the fundamental steps of the study of dendrochronology. ARSTAN is the most popular program used to accomplish this step, and it has played an important role in the development of dendrochronology. However, ARSTAN uses the Fortran programming language, so users find it difficult to understand and revise the algorithm of the source program to meet their needs. An emerging package of the R language named dplR provides similar functions to ARSTAN. R and dplR’s source code is fully open to the public; thus, it has numerous users. When scholars from different domains communicate and share the methods and results of dendrochronology, it can help them improve those chronologies. In addition, R and dplR have become a good supplement to traditional analysis software. This paper compares the different dendrochronological analysis algorithms and results provided by ARSTAN and dplR with tree ring width data fromPiceacrassifoliaon Helan Mountain, Ningxia Hui Autonomous Region, China. The results show that the two programs calculated exactly the same means and standard deviations. The mean error of the mean sensitivities (MS) and first-order autocorrelations (AC) were 0.005 and 0.008, respectively, but they had a clear conversion relationship. When using the same method for detrending with both types of software, the parameters of fitting curves were generally equal, and the corresponding standard chronologies developed by the two programs had a mean error of only 0.002. However, the residual chronologies were very different. In the time domain, a significant difference was observed in the residual chronologies in the first 20—30 years. In the frequency domain, the residual chronologies created using ARSTAN showed more low frequency information than that created using dplR. For example, the former showed periods of 32 years with higher power than those of dplR. In the common interval analysis, ARSTAN gave a higher expressed population signal (EPS) and signal-to-noise ratio (SNR) of chronologies than dplR. EPS error was 0.4% and SNR error was 30%—40%. By comparing the algorithms of the two programs, we found that ARSTAN and dplR have different initial value setting rules and nonlinear fitting methods to choose the best fitting model during detrending. When fitting an autoregression model, ARSTAN used a pooled algorithm to find the integral growing pattern and used the same fitting order for different sequences. However, dplR directly used different optimal fitting models for different sequences. In addition, the two programs used different, but similar, formulas for calculating MS, AC, EPS, and SNR. Although the absolute value of the results was different, calculation results of the same program using different data were comparable. In conclusion, this paper offers two suggestions for the meta-analysis of tree ring data from different sources. First, if the source data are available, researchers should choose a single program for statistical calculation, detrending, and common interval analysis based on their needs. Second, if the source data are not available, information related to the chronologies is sufficient; researchers should use only a single program to calculate EPS and SNR chronological statistics to ensure that the results will be comparable.

R language; dplR; detrending; common interval analysis; dendrochronology

國家自然科學(xué)基金資助項(xiàng)目(41171067)

2014- 03- 30;

日期:2015- 04- 20

10.5846/stxb201403300597

*通訊作者Corresponding author.E-mail: jiangy@bnu.edu.cn

趙守棟，江源，焦亮，王明昌，張凌楠，李文卿.ARSTAN程序和R語言dplR擴(kuò)展包進(jìn)行樹輪年表分析的比較研究.生態(tài)學(xué)報(bào),2015,35(22):7494- 7502.

Zhao S D, Jiang Y, Jiao L, Wang M C, Zhang L N, Li W Q.A Comparative analysis of ARSTAN and the dplr package of R language in analyses of tree-ring chronologies.Acta Ecologica Sinica,2015,35(22):7494- 7502.