• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

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

    2015-02-05 06:03:50趙守棟王明昌張凌楠李文卿
    生態(tài)學(xué)報 2015年22期
    關(guān)鍵詞:趨勢程序差異

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

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

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

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

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

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

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

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

    R語言是當(dāng)前發(fā)展較快的統(tǒng)計分析工具之一,在眾多研究領(lǐng)域得到了越來越廣泛的應(yīng)用。Bunn等使用R語言編寫了dplR擴(kuò)展包,實現(xiàn)了樹輪數(shù)據(jù)讀取、轉(zhuǎn)換、分析、繪圖等功能[6]。相比于經(jīng)典工具而言,R語言和dplR等相關(guān)擴(kuò)展包的源代碼完全公開,且容易理解,研究者可以根據(jù)自身需要修改或編寫程序;同時使用者眾多,網(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]、擾動事件重建[13]等相關(guān)擴(kuò)展包,逐漸形成了一套較為完整的樹輪數(shù)據(jù)分析軟件體系。而在國內(nèi)樹輪年代學(xué)領(lǐng)域,使用dplR開展的研究還較少,使用的功能也相對單一[14]。作為一個新興的樹輪年表分析工具,dplR處理結(jié)果與傳統(tǒng)程序一致與否,將直接影響到不同來源研究數(shù)據(jù)的可比性,因此有必要對ARSTAN和dplR在算法和結(jié)果上進(jìn)行比較分析。

    本文使用賀蘭山青海云杉樹輪寬度數(shù)據(jù),分別借助ARSTAN和dplR進(jìn)行原始序列統(tǒng)計量計算、去趨勢處理、年表建立、公共區(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程序檢驗交叉定年質(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é)樹木年輪實驗室(www.ldeo.columbia.edu/tree-ring-laboratory),版本號為44。選擇的去趨勢方法為修正的負(fù)指數(shù)函數(shù),即:

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

    (1)

    式中,f(i)為生長趨勢估計值,i為年份的序數(shù)。設(shè)定約束條件為d>0,擬合失敗則擬合斜率非正的線性函數(shù)。公共區(qū)間設(shè)定為所有樹芯均含有的共同時間段,即1966—2009年。ARSTAN最終建立標(biāo)準(zhǔn)年表(standard chronology, STD)、差值年表(residual chronology, RES)和自回歸年表(autoregressive standard chronology, ARS),并給出年表的基本統(tǒng)計量和公共區(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ù)計算原始樹輪寬度序列和樹輪寬度指數(shù)年表的基本統(tǒng)計量。使用“detrend.series”函數(shù)對原始樹輪寬度序列逐一進(jìn)行去趨勢操作,確保各樣芯使用的去趨勢方法與ARSTAN中相同,使得結(jié)果具有可比性。使用“rwi.stats”函數(shù)對樹輪寬度指數(shù)序列進(jìn)行公共區(qū)間分析,時間段同樣設(shè)定為1966—2009年,這也是程序默認(rèn)的公共區(qū)間。使用“chron”函數(shù)采用雙權(quán)重平均法建立年表,最終得到STD和RES年表。

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

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

    2 結(jié)果

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

    針對50個原始輪寬序列,ARSTAN和dplR計算所得均值和標(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)一步計算發(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年以上時,相對誤差可以降至1%以下。

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

    ARSTAN和dplR擬合所得負(fù)指數(shù)函數(shù)的參數(shù)a、b完全相同,參數(shù)d的誤差在0.0003以內(nèi),而線性函數(shù)參數(shù)完全相同。去趨勢所得各樣芯的兩個標(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)。

    在建立差值年表時,ARSTAN對各標(biāo)準(zhǔn)化序列均擬合為一階自回歸模型,所得差值序列與原序列長度相等;而dplR則擬合成階數(shù)各異的自回歸模型,包括29個一階模型和17個二階或更高階模型,序列初始損失的年份數(shù)與模型階數(shù)相等,另有4個樣芯未擬合成自回歸模型。各樣芯的兩個差值序列間的平均誤差為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。兩個RES年表僅在起始的10年內(nèi)存在相位不一致的情況,例如1887—1889年(圖1)。

    對兩個RES年表分別進(jìn)行Morlet小波變換(圖2)。結(jié)果顯示,兩個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)計量和公共區(qū)間分析結(jié)果

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

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

    3 討論

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

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

    (4)

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

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

    (5)

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

    (6)

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

    盡管擬合生長趨勢的方法相同,ARSTAN和dplR建立的標(biāo)準(zhǔn)年表仍存在一定差異,這與缺失值的處理方式有關(guān)。為防止運(yùn)算中可能出現(xiàn)的錯誤,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)化序列擬合出一個p階合并自回歸模型,將擬合值作為樹木群體所共有的持續(xù)性生長量。再對各標(biāo)準(zhǔn)化序列分別擬合限定為p階的自回歸模型,將擬合值作為整體生長模式下個體的持續(xù)性生長量,而將殘差作為差值序列,并使用雙權(quán)重平均法建立RES年表。如果將群體的持續(xù)性生長量加回到RES年表上,則得到ARS年表[19]。dplR未考慮整體的生長模式,而是直接對各標(biāo)準(zhǔn)化序列分別擬合自回歸模型,不同序列的模型可以具有不同的階數(shù),再將殘差作為差值序列,建立RES年表。針對本文所用數(shù)據(jù),dplR的計算過程對個體特有的持續(xù)性生長量剔除得更為徹底,低頻信息保留得更少,這種差異在頻域上更為明顯(圖2)。

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

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

    公共區(qū)間分析是指選擇包含特定時間范圍的樹輪寬度序列,通過計算年表信號,評價最終所得年表對森林總體生長的代表性,其中年表信號是衡量序列中包含的共同變化的統(tǒng)計量[21]。DOS版ARSTAN在導(dǎo)入測量數(shù)據(jù)時會提供一個推薦的公共區(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)的時間范圍。

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

    (7)

    (8)

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

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

    目前的樹輪年代學(xué)研究中,一棵樣樹往往取兩個甚至多個樣芯,因此dplR采用Reff計算EPS和SNR,計算結(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é)果時應(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具有更多的低頻信息;在時域上,起始年份處理方式的不同,造成年表在最初的20—30年內(nèi)存在較明顯的差異。不同程序在MS、AC等指標(biāo)上可能存在計算方法差異;ARSTAN使用Rtot計算EPS和SNR等指標(biāo),相對于文獻(xiàn)中基于Reff的算法而言,可能高估了年表中包含的共同信息[25- 26]。

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

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

    (2)如果僅能獲得文獻(xiàn)資料,建議記錄其處理過程的詳細(xì)信息,包括軟件的版本信息、去趨勢的具體方法等;當(dāng)信息足夠時,建議根據(jù)同一標(biāo)準(zhǔn)設(shè)定公共區(qū)間,并重新計算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)計量估計,完成了去趨勢處理并建立了標(biāo)準(zhǔn)年表和差值年表,通過公共區(qū)間分析評價了年表的質(zhì)量。分析表明,ARSTAN與dplR在擬合生長趨勢和自回歸模型時存在算法上的較大差異,同時平均敏感度、一階自相關(guān)系數(shù)、樣本總體代表性和信噪比等指標(biāo)的算法也不一致,但存在較為確定的轉(zhuǎn)換關(guān)系。整合分析不同來源的樹輪資料時,應(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é)報, 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é)報, 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é)報: 自然科學(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é)基金資助項目(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é)報,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.

    猜你喜歡
    趨勢程序差異
    相似與差異
    音樂探索(2022年2期)2022-05-30 21:01:37
    趨勢
    找句子差異
    試論我國未決羈押程序的立法完善
    生物為什么會有差異?
    “程序猿”的生活什么樣
    初秋唇妝趨勢
    Coco薇(2017年9期)2017-09-07 21:23:49
    英國與歐盟正式啟動“離婚”程序程序
    SPINEXPO?2017春夏流行趨勢
    創(chuàng)衛(wèi)暗訪程序有待改進(jìn)
    另类精品久久| 国产亚洲av片在线观看秒播厂| 80岁老熟妇乱子伦牲交| 你懂的网址亚洲精品在线观看| 国产色视频综合| 国产精品欧美亚洲77777| 好男人视频免费观看在线| 激情视频va一区二区三区| 人妻 亚洲 视频| 久久99一区二区三区| 黄色片一级片一级黄色片| 欧美日韩视频高清一区二区三区二| 一二三四在线观看免费中文在| 亚洲欧美色中文字幕在线| 免费高清在线观看日韩| 成年人黄色毛片网站| 国产人伦9x9x在线观看| 国产精品一二三区在线看| 午夜免费观看性视频| 啦啦啦 在线观看视频| 久久亚洲国产成人精品v| 麻豆av在线久日| 日韩人妻精品一区2区三区| 成人黄色视频免费在线看| 十分钟在线观看高清视频www| 在线观看免费日韩欧美大片| 国产熟女午夜一区二区三区| 狂野欧美激情性xxxx| 又大又黄又爽视频免费| 别揉我奶头~嗯~啊~动态视频 | 日本wwww免费看| 国产成人av教育| 午夜激情久久久久久久| 日韩av免费高清视频| 在线观看国产h片| 久久久久久久久免费视频了| 亚洲精品成人av观看孕妇| 亚洲三区欧美一区| 新久久久久国产一级毛片| 国产欧美日韩精品亚洲av| 久久久亚洲精品成人影院| bbb黄色大片| 亚洲国产中文字幕在线视频| 七月丁香在线播放| 91麻豆av在线| 9191精品国产免费久久| 国产一级毛片在线| 男人舔女人的私密视频| 亚洲视频免费观看视频| 免费观看av网站的网址| 国产成人精品无人区| 亚洲国产av新网站| 精品第一国产精品| 91精品伊人久久大香线蕉| 国产日韩欧美视频二区| 国产又爽黄色视频| 美女中出高潮动态图| 99久久人妻综合| 亚洲精品久久成人aⅴ小说| 亚洲成人国产一区在线观看 | 老司机在亚洲福利影院| 一级a爱视频在线免费观看| 欧美亚洲日本最大视频资源| 91字幕亚洲| 日韩电影二区| 国产精品一区二区在线观看99| 亚洲美女黄色视频免费看| 亚洲欧洲精品一区二区精品久久久| 欧美中文综合在线视频| 一本久久精品| 亚洲三区欧美一区| 国产精品国产av在线观看| 欧美日韩亚洲综合一区二区三区_| 欧美日韩亚洲综合一区二区三区_| 满18在线观看网站| 成年人黄色毛片网站| 久久久精品94久久精品| 在线 av 中文字幕| 国产精品麻豆人妻色哟哟久久| 精品久久久久久久毛片微露脸 | 18禁国产床啪视频网站| 亚洲欧美精品自产自拍| 只有这里有精品99| 一本—道久久a久久精品蜜桃钙片| 亚洲精品久久午夜乱码| 一边摸一边做爽爽视频免费| 新久久久久国产一级毛片| 国产欧美日韩一区二区三 | 国产成人一区二区在线| 精品久久久久久久毛片微露脸 | 国产有黄有色有爽视频| 国产男女超爽视频在线观看| 黄色片一级片一级黄色片| 成人18禁高潮啪啪吃奶动态图| 国产精品免费视频内射| 精品国产一区二区三区四区第35| 中文字幕人妻丝袜制服| 只有这里有精品99| 一区二区av电影网| 午夜福利视频精品| 飞空精品影院首页| 尾随美女入室| 天天躁夜夜躁狠狠久久av| 欧美日韩视频精品一区| 日本a在线网址| 少妇粗大呻吟视频| 免费黄频网站在线观看国产| 欧美黑人精品巨大| 久久女婷五月综合色啪小说| 老鸭窝网址在线观看| 少妇人妻 视频| 伦理电影免费视频| 91九色精品人成在线观看| 亚洲成人免费电影在线观看 | 丝袜人妻中文字幕| 欧美黄色淫秽网站| 亚洲国产精品一区三区| 色播在线永久视频| 女性被躁到高潮视频| 亚洲精品在线美女| 国产xxxxx性猛交| 十八禁网站网址无遮挡| 国产精品 国内视频| 亚洲成国产人片在线观看| av国产精品久久久久影院| 五月天丁香电影| 午夜免费男女啪啪视频观看| 精品国产一区二区久久| 国产精品99久久99久久久不卡| 蜜桃在线观看..| 麻豆国产av国片精品| 丝袜美足系列| 18禁国产床啪视频网站| 七月丁香在线播放| 久久久久久免费高清国产稀缺| 黄色视频不卡| 熟女少妇亚洲综合色aaa.| 欧美成人午夜精品| 成人亚洲欧美一区二区av| 国产男人的电影天堂91| 国产免费福利视频在线观看| 一级毛片黄色毛片免费观看视频| 丝袜人妻中文字幕| 女警被强在线播放| 精品一区在线观看国产| 午夜免费鲁丝| 中文字幕精品免费在线观看视频| 黄色视频在线播放观看不卡| 欧美 日韩 精品 国产| 自线自在国产av| 在线观看免费视频网站a站| 国产在线免费精品| 日韩av不卡免费在线播放| 操出白浆在线播放| 欧美久久黑人一区二区| 热99国产精品久久久久久7| 亚洲精品国产av蜜桃| 亚洲,欧美精品.| 麻豆av在线久日| 两人在一起打扑克的视频| 9色porny在线观看| 亚洲午夜精品一区,二区,三区| 亚洲av在线观看美女高潮| 国产一区有黄有色的免费视频| 成年人免费黄色播放视频| 婷婷色麻豆天堂久久| 成人国语在线视频| 黄色毛片三级朝国网站| 日韩一区二区三区影片| 精品人妻一区二区三区麻豆| 人妻一区二区av| 99国产精品一区二区蜜桃av | 国产亚洲精品久久久久5区| 国产成人av教育| 色视频在线一区二区三区| 狂野欧美激情性bbbbbb| 亚洲国产精品国产精品| 国产精品久久久人人做人人爽| 国产高清不卡午夜福利| 嫁个100分男人电影在线观看 | 免费黄频网站在线观看国产| 国产免费现黄频在线看| 成人黄色视频免费在线看| 18禁裸乳无遮挡动漫免费视频| 大片免费播放器 马上看| 亚洲一区中文字幕在线| 青春草视频在线免费观看| 成在线人永久免费视频| 99久久精品国产亚洲精品| 亚洲情色 制服丝袜| 欧美精品啪啪一区二区三区 | 日韩一区二区三区影片| 无遮挡黄片免费观看| 波野结衣二区三区在线| 菩萨蛮人人尽说江南好唐韦庄| 日日夜夜操网爽| 亚洲欧美日韩另类电影网站| 久久久精品免费免费高清| 最新的欧美精品一区二区| 9色porny在线观看| 操美女的视频在线观看| 香蕉国产在线看| 国产成人av教育| 国产精品一区二区免费欧美 | 亚洲精品国产av成人精品| 99久久人妻综合| 老汉色av国产亚洲站长工具| 日韩av在线免费看完整版不卡| 777米奇影视久久| 一级毛片我不卡| 男人舔女人的私密视频| 久久精品人人爽人人爽视色| 汤姆久久久久久久影院中文字幕| 一本色道久久久久久精品综合| netflix在线观看网站| 亚洲精品第二区| 80岁老熟妇乱子伦牲交| 日韩人妻精品一区2区三区| 欧美国产精品一级二级三级| 成年人午夜在线观看视频| 日韩 欧美 亚洲 中文字幕| 一级毛片黄色毛片免费观看视频| av福利片在线| 久热爱精品视频在线9| 性高湖久久久久久久久免费观看| 久久久精品94久久精品| 国产xxxxx性猛交| 亚洲人成电影免费在线| 激情视频va一区二区三区| 中文字幕人妻丝袜制服| 国产免费视频播放在线视频| 国产成人一区二区在线| 欧美国产精品一级二级三级| 午夜免费成人在线视频| 久久性视频一级片| 欧美日韩精品网址| 视频区欧美日本亚洲| 欧美国产精品一级二级三级| 午夜福利视频在线观看免费| 亚洲美女黄色视频免费看| 一本大道久久a久久精品| 亚洲精品一二三| 性少妇av在线| 91字幕亚洲| 18在线观看网站| 老司机深夜福利视频在线观看 | 香蕉国产在线看| 国产精品熟女久久久久浪| 亚洲国产日韩一区二区| 国产人伦9x9x在线观看| 国产成人系列免费观看| 亚洲国产av影院在线观看| 麻豆乱淫一区二区| 伊人亚洲综合成人网| 亚洲国产欧美网| 国产91精品成人一区二区三区 | 免费不卡黄色视频| 纯流量卡能插随身wifi吗| 伊人亚洲综合成人网| 咕卡用的链子| 波野结衣二区三区在线| 亚洲,欧美,日韩| 亚洲成国产人片在线观看| 久久久久久久大尺度免费视频| 欧美在线一区亚洲| 9191精品国产免费久久| 五月开心婷婷网| 精品第一国产精品| 欧美97在线视频| 亚洲欧洲国产日韩| 免费在线观看影片大全网站 | 国产精品久久久久成人av| 久久鲁丝午夜福利片| 久久久久精品国产欧美久久久 | 黄片小视频在线播放| 国产精品人妻久久久影院| 国产高清视频在线播放一区 | 国产一区二区 视频在线| 国产在线观看jvid| 在线av久久热| 久久久久久免费高清国产稀缺| 中文字幕人妻丝袜一区二区| 在线观看免费午夜福利视频| 男男h啪啪无遮挡| 午夜福利视频在线观看免费| 伊人亚洲综合成人网| 亚洲av电影在线观看一区二区三区| 亚洲精品国产区一区二| 久久精品人人爽人人爽视色| 中文字幕色久视频| 水蜜桃什么品种好| 久久av网站| 国产国语露脸激情在线看| 精品久久久精品久久久| 亚洲精品久久成人aⅴ小说| 成年人黄色毛片网站| 91九色精品人成在线观看| 两个人免费观看高清视频| 久久九九热精品免费| 在线观看免费高清a一片| videos熟女内射| 香蕉丝袜av| 亚洲欧洲日产国产| 晚上一个人看的免费电影| 肉色欧美久久久久久久蜜桃| 夫妻午夜视频| 亚洲国产欧美日韩在线播放| 色精品久久人妻99蜜桃| 亚洲专区国产一区二区| 中国美女看黄片| www日本在线高清视频| 美女脱内裤让男人舔精品视频| av在线播放精品| 亚洲精品一二三| 国产精品 欧美亚洲| 国产一区有黄有色的免费视频| 菩萨蛮人人尽说江南好唐韦庄| 女人被躁到高潮嗷嗷叫费观| 69精品国产乱码久久久| 国产一级毛片在线| 69精品国产乱码久久久| 久久久久久久大尺度免费视频| 成人手机av| 一边摸一边做爽爽视频免费| 免费人妻精品一区二区三区视频| 国产日韩欧美视频二区| 超色免费av| 国产精品一区二区精品视频观看| 欧美精品一区二区大全| 热99国产精品久久久久久7| 久久中文字幕一级| 亚洲国产欧美网| 五月开心婷婷网| 亚洲色图综合在线观看| 久久鲁丝午夜福利片| 久久久久久久久免费视频了| www.999成人在线观看| 久久热在线av| 午夜福利,免费看| 多毛熟女@视频| 日韩视频在线欧美| 成人三级做爰电影| 久久狼人影院| 午夜久久久在线观看| 日本欧美视频一区| 美女扒开内裤让男人捅视频| 精品一区二区三卡| 夫妻午夜视频| 国产91精品成人一区二区三区 | 777米奇影视久久| 欧美av亚洲av综合av国产av| 亚洲免费av在线视频| 亚洲少妇的诱惑av| 欧美黄色淫秽网站| 国产成人系列免费观看| 熟女少妇亚洲综合色aaa.| 精品国产一区二区三区四区第35| 中文字幕色久视频| 国产精品香港三级国产av潘金莲 | 亚洲国产精品999| 久久影院123| 1024香蕉在线观看| 9色porny在线观看| 亚洲国产日韩一区二区| 亚洲成人手机| 免费观看人在逋| 一本综合久久免费| 高清黄色对白视频在线免费看| 又黄又粗又硬又大视频| 中文字幕人妻丝袜制服| 一边摸一边抽搐一进一出视频| 国产成人精品久久久久久| 热99国产精品久久久久久7| 亚洲欧美精品自产自拍| 水蜜桃什么品种好| 亚洲五月婷婷丁香| 蜜桃国产av成人99| av福利片在线| 精品亚洲乱码少妇综合久久| 日韩大片免费观看网站| 国产淫语在线视频| 久久久精品94久久精品| 亚洲av美国av| 久久青草综合色| 黄色一级大片看看| 热re99久久国产66热| 国产精品一区二区在线观看99| 99久久99久久久精品蜜桃| 90打野战视频偷拍视频| 好男人视频免费观看在线| 成在线人永久免费视频| 老司机亚洲免费影院| 欧美在线一区亚洲| 午夜两性在线视频| 欧美日韩一级在线毛片| 乱人伦中国视频| 女人精品久久久久毛片| 国产亚洲精品第一综合不卡| 国产亚洲欧美精品永久| 少妇被粗大的猛进出69影院| 50天的宝宝边吃奶边哭怎么回事| 看十八女毛片水多多多| 欧美人与性动交α欧美软件| 欧美在线一区亚洲| 亚洲熟女毛片儿| 欧美日韩视频精品一区| 久久精品国产a三级三级三级| 亚洲国产av新网站| 免费少妇av软件| 大片电影免费在线观看免费| kizo精华| 久久青草综合色| 久久久国产精品麻豆| 午夜日韩欧美国产| 男女免费视频国产| 国产精品免费视频内射| 免费人妻精品一区二区三区视频| 亚洲精品国产av蜜桃| 日本vs欧美在线观看视频| 丝袜喷水一区| 自拍欧美九色日韩亚洲蝌蚪91| 在线观看免费视频网站a站| 欧美精品av麻豆av| 亚洲久久久国产精品| 久久青草综合色| 人人妻人人澡人人爽人人夜夜| av国产久精品久网站免费入址| 成年美女黄网站色视频大全免费| 色播在线永久视频| 美女高潮到喷水免费观看| xxxhd国产人妻xxx| 波野结衣二区三区在线| 成年人午夜在线观看视频| 亚洲精品自拍成人| 蜜桃国产av成人99| 大香蕉久久成人网| 亚洲五月色婷婷综合| 久久99精品国语久久久| 热re99久久精品国产66热6| 母亲3免费完整高清在线观看| 无限看片的www在线观看| 欧美在线一区亚洲| 免费看av在线观看网站| 日韩一卡2卡3卡4卡2021年| 在线观看一区二区三区激情| www.精华液| 国产一卡二卡三卡精品| 少妇猛男粗大的猛烈进出视频| 大型av网站在线播放| 欧美精品啪啪一区二区三区 | 女警被强在线播放| 国产精品久久久av美女十八| 男女高潮啪啪啪动态图| 亚洲精品美女久久久久99蜜臀 | 日韩一本色道免费dvd| 精品国产超薄肉色丝袜足j| 美女主播在线视频| av在线老鸭窝| 免费av中文字幕在线| 宅男免费午夜| 久久 成人 亚洲| 国产av精品麻豆| 国产老妇伦熟女老妇高清| 1024视频免费在线观看| 悠悠久久av| 久久久久久亚洲精品国产蜜桃av| 欧美变态另类bdsm刘玥| 亚洲精品久久午夜乱码| 欧美激情高清一区二区三区| 亚洲欧美精品综合一区二区三区| 如日韩欧美国产精品一区二区三区| 每晚都被弄得嗷嗷叫到高潮| 桃花免费在线播放| 另类亚洲欧美激情| 少妇猛男粗大的猛烈进出视频| 国产伦人伦偷精品视频| 又紧又爽又黄一区二区| 操出白浆在线播放| 日韩视频在线欧美| 黄片播放在线免费| 亚洲九九香蕉| 欧美成人午夜精品| 国产一区二区在线观看av| 青春草视频在线免费观看| 黑人欧美特级aaaaaa片| 最近中文字幕2019免费版| 亚洲午夜精品一区,二区,三区| 两性夫妻黄色片| 韩国精品一区二区三区| 国产成人啪精品午夜网站| 国产欧美日韩精品亚洲av| 一本久久精品| 中国美女看黄片| 18禁观看日本| 2018国产大陆天天弄谢| 亚洲精品美女久久久久99蜜臀 | 18在线观看网站| 伊人久久大香线蕉亚洲五| 啦啦啦在线观看免费高清www| 亚洲男人天堂网一区| 欧美日韩国产mv在线观看视频| cao死你这个sao货| 每晚都被弄得嗷嗷叫到高潮| xxx大片免费视频| 一本一本久久a久久精品综合妖精| 9191精品国产免费久久| 亚洲一区二区三区欧美精品| 18禁观看日本| 午夜久久久在线观看| 国产一区二区三区综合在线观看| 国产精品熟女久久久久浪| 青草久久国产| 精品亚洲成国产av| av天堂久久9| 免费高清在线观看日韩| 又紧又爽又黄一区二区| 亚洲伊人久久精品综合| 国产又色又爽无遮挡免| 精品国产国语对白av| 麻豆av在线久日| 国产深夜福利视频在线观看| www.精华液| 美女高潮到喷水免费观看| 美女国产高潮福利片在线看| 日本色播在线视频| 欧美av亚洲av综合av国产av| 黄色a级毛片大全视频| 国产亚洲精品久久久久5区| 一区二区三区激情视频| 亚洲av男天堂| 国产伦理片在线播放av一区| 亚洲欧美清纯卡通| 日本av免费视频播放| 欧美精品高潮呻吟av久久| 青青草视频在线视频观看| 免费观看av网站的网址| 亚洲精品中文字幕在线视频| 色综合欧美亚洲国产小说| 日本黄色日本黄色录像| 老司机影院成人| 亚洲一卡2卡3卡4卡5卡精品中文| 在线观看一区二区三区激情| 免费在线观看日本一区| 在线看a的网站| 少妇粗大呻吟视频| 国产精品免费大片| 久久久国产精品麻豆| 一级a爱视频在线免费观看| 久久这里只有精品19| 国产又色又爽无遮挡免| 国语对白做爰xxxⅹ性视频网站| 久久精品成人免费网站| 精品亚洲乱码少妇综合久久| 亚洲av电影在线观看一区二区三区| 久久中文字幕一级| 婷婷色av中文字幕| 一级片'在线观看视频| 国产在线免费精品| 日韩 亚洲 欧美在线| 老司机午夜十八禁免费视频| 尾随美女入室| 欧美国产精品va在线观看不卡| 亚洲天堂av无毛| 久久精品熟女亚洲av麻豆精品| av在线播放精品| 日韩人妻精品一区2区三区| 久久精品久久久久久久性| 欧美日韩视频精品一区| 色精品久久人妻99蜜桃| 9色porny在线观看| av有码第一页| 欧美乱码精品一区二区三区| 亚洲av欧美aⅴ国产| 久久精品久久久久久久性| 亚洲国产av新网站| 久久国产精品男人的天堂亚洲| 久久午夜综合久久蜜桃| 视频区欧美日本亚洲| 热99国产精品久久久久久7| 日本a在线网址| 美女国产高潮福利片在线看| 国产精品成人在线| 日本五十路高清| 国产精品三级大全| 91精品三级在线观看| 男男h啪啪无遮挡| 永久免费av网站大全| 激情视频va一区二区三区| 日本wwww免费看| 国产主播在线观看一区二区 | av在线app专区| 国产精品三级大全| 99re6热这里在线精品视频| 亚洲国产最新在线播放| cao死你这个sao货| 久热爱精品视频在线9| 观看av在线不卡| 晚上一个人看的免费电影| 日韩一卡2卡3卡4卡2021年| 日韩电影二区| xxxhd国产人妻xxx| 丝袜美腿诱惑在线| 欧美精品亚洲一区二区| 国产一区二区三区av在线| 91字幕亚洲| 国产成人欧美在线观看 | 久久亚洲国产成人精品v| 国产1区2区3区精品| 免费观看av网站的网址| 亚洲精品成人av观看孕妇| 99热国产这里只有精品6| 国产成人系列免费观看| 中文字幕最新亚洲高清| 九草在线视频观看| 久久精品亚洲熟妇少妇任你|