• <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é)報(bào) 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ò)展包實(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.

    猜你喜歡
    趨勢程序差異
    相似與差異
    音樂探索(2022年2期)2022-05-30 21:01:37
    趨勢
    找句子差異
    試論我國未決羈押程序的立法完善
    生物為什么會(huì)有差異?
    “程序猿”的生活什么樣
    初秋唇妝趨勢
    Coco薇(2017年9期)2017-09-07 21:23:49
    英國與歐盟正式啟動(dòng)“離婚”程序程序
    SPINEXPO?2017春夏流行趨勢
    創(chuàng)衛(wèi)暗訪程序有待改進(jìn)
    亚洲成人精品中文字幕电影 | 精品一区二区三卡| 精品乱码久久久久久99久播| 免费看十八禁软件| 久久久精品国产亚洲av高清涩受| 黄色a级毛片大全视频| 欧美日韩瑟瑟在线播放| 亚洲美女黄片视频| 午夜亚洲福利在线播放| 首页视频小说图片口味搜索| 亚洲精品av麻豆狂野| 九色亚洲精品在线播放| 日韩精品免费视频一区二区三区| 日韩人妻精品一区2区三区| 中出人妻视频一区二区| 曰老女人黄片| 免费av毛片视频| 琪琪午夜伦伦电影理论片6080| 国产精品免费一区二区三区在线| 无限看片的www在线观看| 色老头精品视频在线观看| 侵犯人妻中文字幕一二三四区| 女同久久另类99精品国产91| 每晚都被弄得嗷嗷叫到高潮| 欧美一区二区精品小视频在线| 亚洲成人精品中文字幕电影 | 亚洲av日韩精品久久久久久密| 日韩视频一区二区在线观看| 亚洲av熟女| 91国产中文字幕| 午夜老司机福利片| 成人三级做爰电影| xxx96com| 久久久久久久久中文| 日韩免费高清中文字幕av| 中文字幕色久视频| 精品福利永久在线观看| 丝袜人妻中文字幕| 99精品久久久久人妻精品| 欧洲精品卡2卡3卡4卡5卡区| xxx96com| 久久久久久久久中文| 欧美成人午夜精品| 婷婷精品国产亚洲av在线| av网站在线播放免费| netflix在线观看网站| 国产精品亚洲一级av第二区| 亚洲自拍偷在线| 亚洲成av片中文字幕在线观看| 一进一出好大好爽视频| 日本黄色日本黄色录像| av视频免费观看在线观看| 国产一区二区在线av高清观看| 国产成人精品在线电影| 狠狠狠狠99中文字幕| 99riav亚洲国产免费| 亚洲人成伊人成综合网2020| 中文字幕最新亚洲高清| 99国产精品99久久久久| 老熟妇乱子伦视频在线观看| 身体一侧抽搐| 亚洲国产精品sss在线观看 | 国产91精品成人一区二区三区| 视频在线观看一区二区三区| 人妻丰满熟妇av一区二区三区| 欧美成人午夜精品| 美女高潮到喷水免费观看| 操出白浆在线播放| 老司机深夜福利视频在线观看| 91成人精品电影| 国内毛片毛片毛片毛片毛片| 国产精品久久电影中文字幕| 国产精品98久久久久久宅男小说| 18禁国产床啪视频网站| 欧美人与性动交α欧美精品济南到| 久久午夜综合久久蜜桃| 欧美日韩黄片免| 宅男免费午夜| 男人操女人黄网站| 两个人看的免费小视频| 欧美不卡视频在线免费观看 | 欧美日韩视频精品一区| 国产精品一区二区在线不卡| 丰满饥渴人妻一区二区三| 免费av毛片视频| 久久久久国内视频| 亚洲欧美日韩高清在线视频| 一二三四在线观看免费中文在| 亚洲精品久久成人aⅴ小说| 满18在线观看网站| 中文字幕最新亚洲高清| 99久久久亚洲精品蜜臀av| 88av欧美| 很黄的视频免费| 看黄色毛片网站| 男人舔女人的私密视频| 国内毛片毛片毛片毛片毛片| 精品久久久久久成人av| 夜夜看夜夜爽夜夜摸 | 精品国产一区二区久久| 日日摸夜夜添夜夜添小说| 午夜视频精品福利| 黄色视频,在线免费观看| 夜夜夜夜夜久久久久| 亚洲国产欧美网| 日韩欧美国产一区二区入口| 性色av乱码一区二区三区2| 可以免费在线观看a视频的电影网站| 两性午夜刺激爽爽歪歪视频在线观看 | 久久久久久人人人人人| 国产免费现黄频在线看| 精品午夜福利视频在线观看一区| 曰老女人黄片| 91精品三级在线观看| 后天国语完整版免费观看| 久久人人爽av亚洲精品天堂| av天堂在线播放| 俄罗斯特黄特色一大片| 自拍欧美九色日韩亚洲蝌蚪91| 色婷婷久久久亚洲欧美| 男人舔女人下体高潮全视频| 香蕉国产在线看| 久久 成人 亚洲| 精品人妻1区二区| 精品免费久久久久久久清纯| 黄片播放在线免费| 国产野战对白在线观看| 大香蕉久久成人网| 搡老熟女国产l中国老女人| 日本欧美视频一区| 国产av又大| 真人做人爱边吃奶动态| 男人舔女人下体高潮全视频| 午夜精品久久久久久毛片777| 91在线观看av| 1024香蕉在线观看| 老司机午夜十八禁免费视频| 欧美不卡视频在线免费观看 | 国产熟女xx| avwww免费| 国产精品成人在线| 欧美性长视频在线观看| 亚洲久久久国产精品| 国产成人av教育| 一本综合久久免费| 美女高潮喷水抽搐中文字幕| 无人区码免费观看不卡| 一级a爱视频在线免费观看| 人人妻人人爽人人添夜夜欢视频| www国产在线视频色| 不卡一级毛片| 久久人妻熟女aⅴ| 日韩有码中文字幕| 免费在线观看日本一区| 精品一区二区三区视频在线观看免费 | 99国产极品粉嫩在线观看| 99精国产麻豆久久婷婷| 国产精品亚洲av一区麻豆| 两性午夜刺激爽爽歪歪视频在线观看 | 黄网站色视频无遮挡免费观看| 精品国产一区二区三区四区第35| 成年人免费黄色播放视频| 啪啪无遮挡十八禁网站| 在线av久久热| 亚洲情色 制服丝袜| 欧美 亚洲 国产 日韩一| 丝袜在线中文字幕| 免费av毛片视频| 日日干狠狠操夜夜爽| 老司机午夜十八禁免费视频| 国产精品电影一区二区三区| 国产精品免费视频内射| 99re在线观看精品视频| 国产三级黄色录像| tocl精华| xxx96com| 亚洲欧美激情在线| 国产又爽黄色视频| 亚洲精品粉嫩美女一区| 18禁裸乳无遮挡免费网站照片 | 成人免费观看视频高清| 在线观看午夜福利视频| 免费不卡黄色视频| 香蕉久久夜色| 午夜影院日韩av| 久久精品国产亚洲av高清一级| 97超级碰碰碰精品色视频在线观看| 热re99久久国产66热| 日本黄色视频三级网站网址| 男人舔女人下体高潮全视频| 欧美在线一区亚洲| 丝袜美腿诱惑在线| 丝袜美腿诱惑在线| 亚洲欧美激情综合另类| 好看av亚洲va欧美ⅴa在| 美女午夜性视频免费| 亚洲少妇的诱惑av| 欧美激情 高清一区二区三区| 亚洲第一青青草原| 高清黄色对白视频在线免费看| 黄色片一级片一级黄色片| 高清av免费在线| 久久久久久久久免费视频了| 亚洲成人免费电影在线观看| 欧美黄色片欧美黄色片| 成年人黄色毛片网站| 亚洲精品久久午夜乱码| 久久婷婷成人综合色麻豆| 国产欧美日韩精品亚洲av| 91国产中文字幕| 丰满人妻熟妇乱又伦精品不卡| 国产1区2区3区精品| 好男人电影高清在线观看| 国产欧美日韩一区二区三| 欧美激情高清一区二区三区| 后天国语完整版免费观看| 一级片'在线观看视频| 欧美 亚洲 国产 日韩一| 91成人精品电影| 亚洲欧美一区二区三区黑人| 国产精品美女特级片免费视频播放器 | 两个人免费观看高清视频| 高清毛片免费观看视频网站 | 激情视频va一区二区三区| 一区福利在线观看| 免费人成视频x8x8入口观看| 狂野欧美激情性xxxx| 9热在线视频观看99| 国产成人欧美| 亚洲国产看品久久| 国产成人精品在线电影| 国产欧美日韩一区二区精品| 啦啦啦在线免费观看视频4| 国产深夜福利视频在线观看| 国产午夜精品久久久久久| 午夜老司机福利片| 99国产精品一区二区三区| 亚洲第一av免费看| 中亚洲国语对白在线视频| 黑人猛操日本美女一级片| 欧美精品亚洲一区二区| 日韩 欧美 亚洲 中文字幕| 大香蕉久久成人网| 高清av免费在线| 亚洲中文字幕日韩| 国产蜜桃级精品一区二区三区| 久久久久九九精品影院| 脱女人内裤的视频| 99久久久亚洲精品蜜臀av| 露出奶头的视频| 可以免费在线观看a视频的电影网站| 欧美精品亚洲一区二区| 久久亚洲精品不卡| 制服诱惑二区| 9191精品国产免费久久| av天堂在线播放| 亚洲aⅴ乱码一区二区在线播放 | 午夜免费观看网址| 精品国产一区二区久久| 日本vs欧美在线观看视频| 国产精品99久久99久久久不卡| 91成人精品电影| 成人黄色视频免费在线看| www日本在线高清视频| 国产三级黄色录像| 日韩大尺度精品在线看网址 | 久久精品影院6| 久99久视频精品免费| 黄色成人免费大全| 日日干狠狠操夜夜爽| 亚洲中文日韩欧美视频| 亚洲av成人av| 成人18禁高潮啪啪吃奶动态图| 国产成人系列免费观看| 最新在线观看一区二区三区| 久久欧美精品欧美久久欧美| 日本一区二区免费在线视频| 久久婷婷成人综合色麻豆| 亚洲五月天丁香| 女人爽到高潮嗷嗷叫在线视频| www.999成人在线观看| 免费久久久久久久精品成人欧美视频| 天天影视国产精品| 最新在线观看一区二区三区| av中文乱码字幕在线| 麻豆国产av国片精品| 亚洲精品一卡2卡三卡4卡5卡| 国内毛片毛片毛片毛片毛片| 91九色精品人成在线观看| 亚洲熟妇中文字幕五十中出 | 国产精品日韩av在线免费观看 | 侵犯人妻中文字幕一二三四区| 欧美日韩乱码在线| 久久精品国产综合久久久| 国产欧美日韩精品亚洲av| 久久久久精品国产欧美久久久| 免费少妇av软件| 三级毛片av免费| 亚洲一区二区三区色噜噜 | 亚洲va日本ⅴa欧美va伊人久久| 女性生殖器流出的白浆| 免费一级毛片在线播放高清视频 | 欧美精品一区二区免费开放| 久久香蕉精品热| 日韩av在线大香蕉| 老司机深夜福利视频在线观看| 国产免费现黄频在线看| 国产99久久九九免费精品| 亚洲视频免费观看视频| √禁漫天堂资源中文www| 亚洲国产毛片av蜜桃av| 亚洲免费av在线视频| 大型黄色视频在线免费观看| 18禁黄网站禁片午夜丰满| 999久久久精品免费观看国产| 亚洲五月天丁香| 岛国在线观看网站| 狠狠狠狠99中文字幕| 女人爽到高潮嗷嗷叫在线视频| 美女高潮到喷水免费观看| 人人澡人人妻人| 1024香蕉在线观看| 中文欧美无线码| 国产精品99久久99久久久不卡| 亚洲av第一区精品v没综合| 午夜91福利影院| 国产亚洲欧美在线一区二区| 人妻丰满熟妇av一区二区三区| 少妇被粗大的猛进出69影院| 成人三级做爰电影| 国产成人系列免费观看| 久久久久久亚洲精品国产蜜桃av| 黑人巨大精品欧美一区二区蜜桃| 日本wwww免费看| 女人被躁到高潮嗷嗷叫费观| 国产精品美女特级片免费视频播放器 | 老熟妇乱子伦视频在线观看| 免费搜索国产男女视频| 在线十欧美十亚洲十日本专区| 黄色视频,在线免费观看| 国产精品野战在线观看 | 亚洲午夜理论影院| 亚洲人成77777在线视频| 看片在线看免费视频| 日韩三级视频一区二区三区| 国产精品一区二区免费欧美| www.999成人在线观看| 国产在线精品亚洲第一网站| 亚洲熟妇中文字幕五十中出 | 俄罗斯特黄特色一大片| 高清欧美精品videossex| 免费高清在线观看日韩| 亚洲精品av麻豆狂野| 国产欧美日韩精品亚洲av| 91老司机精品| 黄色视频不卡| 丝袜美足系列| 亚洲九九香蕉| 欧美人与性动交α欧美软件| 亚洲av熟女| 欧美日韩瑟瑟在线播放| 老鸭窝网址在线观看| 亚洲一区中文字幕在线| 夜夜爽天天搞| 亚洲欧美激情在线| 后天国语完整版免费观看| 欧美亚洲日本最大视频资源| 日韩有码中文字幕| 99国产极品粉嫩在线观看| 很黄的视频免费| 亚洲精品在线观看二区| 十分钟在线观看高清视频www| 亚洲少妇的诱惑av| 大码成人一级视频| 日韩欧美一区视频在线观看| 制服人妻中文乱码| 午夜福利,免费看| 一级片'在线观看视频| 亚洲三区欧美一区| 91在线观看av| 80岁老熟妇乱子伦牲交| 国产区一区二久久| 大型黄色视频在线免费观看| 国产黄色免费在线视频| 亚洲av片天天在线观看| 男女做爰动态图高潮gif福利片 | 久久精品影院6| 成熟少妇高潮喷水视频| 悠悠久久av| 人成视频在线观看免费观看| 99国产精品99久久久久| 在线永久观看黄色视频| 精品无人区乱码1区二区| 欧美成人午夜精品| 露出奶头的视频| 成人18禁在线播放| 亚洲国产欧美日韩在线播放| 亚洲黑人精品在线| 黄色毛片三级朝国网站| 一区福利在线观看| 人人妻人人澡人人看| 亚洲三区欧美一区| 国产成人精品无人区| 精品卡一卡二卡四卡免费| 精品乱码久久久久久99久播| 在线国产一区二区在线| 日韩欧美一区视频在线观看| 久久草成人影院| 午夜精品国产一区二区电影| 国产精华一区二区三区| 国产精品国产高清国产av| 五月开心婷婷网| 精品第一国产精品| 久久中文看片网| 丰满人妻熟妇乱又伦精品不卡| 国产不卡一卡二| 国产精品成人在线| 亚洲精品一卡2卡三卡4卡5卡| 欧美成人性av电影在线观看| 亚洲午夜理论影院| 在线天堂中文资源库| 免费av中文字幕在线| 一区福利在线观看| 中出人妻视频一区二区| 久久人人97超碰香蕉20202| 久久久久九九精品影院| 一级a爱视频在线免费观看| 日本a在线网址| 老司机深夜福利视频在线观看| 国产亚洲欧美在线一区二区| 国产成人av激情在线播放| 午夜福利欧美成人| 国产极品粉嫩免费观看在线| 亚洲人成网站在线播放欧美日韩| 久久人妻福利社区极品人妻图片| 一级毛片女人18水好多| 婷婷六月久久综合丁香| 性欧美人与动物交配| 老熟妇仑乱视频hdxx| 99riav亚洲国产免费| 热99国产精品久久久久久7| 国产精品日韩av在线免费观看 | 看免费av毛片| 国产精品久久久人人做人人爽| 国产91精品成人一区二区三区| 精品久久蜜臀av无| 日韩欧美一区视频在线观看| 亚洲欧美日韩高清在线视频| 极品教师在线免费播放| 亚洲精品一卡2卡三卡4卡5卡| 成年人黄色毛片网站| 欧美+亚洲+日韩+国产| aaaaa片日本免费| 欧美av亚洲av综合av国产av| 亚洲欧洲精品一区二区精品久久久| 日本精品一区二区三区蜜桃| 亚洲欧美一区二区三区久久| 免费搜索国产男女视频| 亚洲精品中文字幕在线视频| 国产97色在线日韩免费| 他把我摸到了高潮在线观看| 国产激情欧美一区二区| 午夜亚洲福利在线播放| 国产亚洲精品久久久久久毛片| 色精品久久人妻99蜜桃| 欧美日韩国产mv在线观看视频| 久久婷婷成人综合色麻豆| 亚洲av熟女| 国产成+人综合+亚洲专区| 美女高潮到喷水免费观看| 国产日韩一区二区三区精品不卡| 免费少妇av软件| 日韩一卡2卡3卡4卡2021年| 午夜视频精品福利| 午夜免费激情av| 日韩精品中文字幕看吧| 国产精品自产拍在线观看55亚洲| 黄色毛片三级朝国网站| 亚洲成人国产一区在线观看| www.精华液| 99re在线观看精品视频| 两人在一起打扑克的视频| 三上悠亚av全集在线观看| 啪啪无遮挡十八禁网站| 免费在线观看视频国产中文字幕亚洲| 久久亚洲真实| 每晚都被弄得嗷嗷叫到高潮| 精品欧美一区二区三区在线| 中文字幕人妻丝袜一区二区| 久久精品亚洲精品国产色婷小说| 高清毛片免费观看视频网站 | 99国产综合亚洲精品| 午夜激情av网站| 一区在线观看完整版| 精品人妻1区二区| 法律面前人人平等表现在哪些方面| 天堂动漫精品| 超碰成人久久| 欧美日韩视频精品一区| 美女大奶头视频| 亚洲情色 制服丝袜| 黄色视频,在线免费观看| 一级a爱视频在线免费观看| 久99久视频精品免费| 欧美在线黄色| xxxhd国产人妻xxx| 欧美日韩乱码在线| 在线观看www视频免费| 亚洲免费av在线视频| 亚洲激情在线av| 怎么达到女性高潮| 久久精品影院6| 欧美人与性动交α欧美精品济南到| 午夜成年电影在线免费观看| 美女大奶头视频| 精品欧美一区二区三区在线| 一级a爱片免费观看的视频| 人人妻人人澡人人看| 久久精品91蜜桃| 丝袜人妻中文字幕| 欧美日韩av久久| 亚洲精品中文字幕一二三四区| 久久精品国产综合久久久| 成人亚洲精品一区在线观看| 看片在线看免费视频| 欧美黑人欧美精品刺激| 黄片播放在线免费| 国产精品免费一区二区三区在线| 熟女少妇亚洲综合色aaa.| 亚洲第一av免费看| 桃色一区二区三区在线观看| ponron亚洲| 人人妻人人澡人人看| www.自偷自拍.com| 亚洲色图 男人天堂 中文字幕| 老汉色∧v一级毛片| 午夜精品在线福利| 丝袜在线中文字幕| 国产成人av教育| 亚洲精品中文字幕在线视频| 大码成人一级视频| 19禁男女啪啪无遮挡网站| 曰老女人黄片| 国产成人av教育| 另类亚洲欧美激情| 女人爽到高潮嗷嗷叫在线视频| 黄片大片在线免费观看| 亚洲国产精品sss在线观看 | 日韩欧美国产一区二区入口| 久久久久精品国产欧美久久久| 免费高清在线观看日韩| 亚洲中文字幕日韩| www国产在线视频色| 99国产精品免费福利视频| 色婷婷久久久亚洲欧美| 日韩有码中文字幕| 国产精品久久久人人做人人爽| 9191精品国产免费久久| a在线观看视频网站| 久久国产乱子伦精品免费另类| 精品一区二区三区av网在线观看| 欧美成人性av电影在线观看| 999久久久国产精品视频| 高清在线国产一区| 欧美日韩亚洲综合一区二区三区_| 亚洲精品在线美女| 成年女人毛片免费观看观看9| 91精品三级在线观看| 亚洲午夜精品一区,二区,三区| 黄色片一级片一级黄色片| 他把我摸到了高潮在线观看| 午夜福利影视在线免费观看| 欧美色视频一区免费| 在线观看免费视频网站a站| 亚洲av片天天在线观看| 国产精华一区二区三区| 精品一区二区三卡| 丁香六月欧美| 国产精品久久视频播放| av中文乱码字幕在线| 50天的宝宝边吃奶边哭怎么回事| 日日干狠狠操夜夜爽| 大型黄色视频在线免费观看| 日韩欧美在线二视频| 亚洲人成77777在线视频| 丰满饥渴人妻一区二区三| 可以免费在线观看a视频的电影网站| 在线天堂中文资源库| 久久国产精品影院| 国产区一区二久久| 久久久水蜜桃国产精品网| 在线观看66精品国产| 欧美中文日本在线观看视频| 日韩三级视频一区二区三区| 少妇 在线观看| 国产单亲对白刺激| 五月开心婷婷网| 黄网站色视频无遮挡免费观看| 精品日产1卡2卡| 欧美日韩中文字幕国产精品一区二区三区 | 亚洲午夜理论影院| 亚洲色图 男人天堂 中文字幕| 一区二区日韩欧美中文字幕| 午夜福利影视在线免费观看| 一级作爱视频免费观看| 亚洲视频免费观看视频| 色婷婷久久久亚洲欧美| 亚洲人成伊人成综合网2020| 国产精品一区二区三区四区久久 | 国产精品一区二区在线不卡| 日日干狠狠操夜夜爽| 久久欧美精品欧美久久欧美| 50天的宝宝边吃奶边哭怎么回事| 亚洲五月天丁香| 成人三级黄色视频|