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

    利用經(jīng)驗(yàn)-過程混合建模方法優(yōu)化華山松過程模型的參數(shù)*

    2021-12-01 01:50:42薛海連田相林曹田健
    林業(yè)科學(xué) 2021年9期
    關(guān)鍵詞:樹高胸徑樣地

    薛海連 田相林 曹田健

    (1.西北農(nóng)林科技大學(xué)理學(xué)院 楊凌 712100; 2.赫爾辛基大學(xué)林學(xué)系 赫爾辛基 FI-00014;3.西北農(nóng)林科技大學(xué)生態(tài)仿真優(yōu)化實(shí)驗(yàn)室 楊凌 712100)

    森林生長模型是森林生態(tài)系統(tǒng)動(dòng)態(tài)模擬與森林經(jīng)營規(guī)劃的重要工具,依據(jù)模型構(gòu)建方式,可分為經(jīng)驗(yàn)?zāi)P?、過程模型和混合模型(Kimmins, 2004)。傳統(tǒng)的經(jīng)驗(yàn)?zāi)P徒R蕾嚧罅繕拥卣{(diào)查數(shù)據(jù),采用統(tǒng)計(jì)回歸技術(shù)預(yù)估模型參數(shù),準(zhǔn)確可靠地預(yù)測林分生長和收獲,但外推至建模數(shù)據(jù)之外的環(huán)境條件時(shí)應(yīng)用受到限制。相比之下,過程模型的參數(shù)通過生理學(xué)觀測獲得,模型形式基于理論性植物與環(huán)境的關(guān)系,可外推計(jì)算環(huán)境變化對(duì)樹木生長的影響(Weiskitteletal., 2011),但其結(jié)構(gòu)復(fù)雜、參數(shù)繁多且參數(shù)值難以實(shí)測確定?;旌夏P筒捎冒虢?jīng)驗(yàn)半機(jī)理的建模方式,一方面在機(jī)理建模的基礎(chǔ)上結(jié)合部分經(jīng)驗(yàn)?zāi)K,另一方面對(duì)部分機(jī)理模塊進(jìn)行簡化,如基于生理學(xué)過程的3-PG(physiological principles in predicting growth)模型(Landsbergetal., 1997)和CABALA(CArbon BALAnce)模型(Battagliaetal., 2004)均是通過集成經(jīng)驗(yàn)的異速生長模塊來簡化樹木生長的機(jī)制模擬。雖然混合模型可簡化模型結(jié)構(gòu)且能夠降低建模數(shù)據(jù)需求,但仍包含動(dòng)輒數(shù)十或上百個(gè)參數(shù)。無論是過程模型還是混合模型,高維變量條件下的參數(shù)預(yù)估方法設(shè)計(jì)與選取一直是模型建模研究的焦點(diǎn)問題。

    過程模型參數(shù)預(yù)估通常涉及樹木生理、結(jié)構(gòu)、功能性狀等信息,這些信息往往具有一定的不確定性,有些參數(shù)甚至無法測量,只能間接推斷。目前,過程模型參數(shù)預(yù)估方法主要包括同時(shí)擬合法(Siev?nenetal., 1993; Vanclayetal., 1997)、Monte Carlo法(Hornbergeretal., 1985)、貝葉斯法(van Oijenetal., 2005; Siviaetal., 2006)和試錯(cuò)法(M?kel?, 1997; 張廷龍等, 2011)等,其中,同時(shí)擬合法、Monte Carlo法和貝葉斯法均需先給定參數(shù)的參考值或可能的變化區(qū)間等先驗(yàn)信息,而后搜索計(jì)算似然值(M?kel?etal., 2000; M?kel?, 1988; van Oijenetal., 2005; Ewen, 2013)。如Minunno等(2019)研究CROBAS(a growth model based on crown and base dynamics)模型參數(shù)的不確定性時(shí),基于芬蘭國家森林資源連續(xù)清查數(shù)據(jù),采用貝葉斯法校正歐洲赤松(Pinussylvestris)、挪威云杉(Piceaabies)和垂枝樺(Betulapendula)的CROBAS模型; 然而,為收斂至參數(shù)的完整聯(lián)合后驗(yàn)分布,貝葉斯法在處理不可解析的高維模型時(shí)消耗大量時(shí)間,并需進(jìn)行收斂性的反復(fù)測試(van Oijenetal., 2005; Hartigetal., 2012)。試錯(cuò)法主要基于實(shí)測數(shù)據(jù)、擬合數(shù)據(jù)、文獻(xiàn)數(shù)據(jù)和專家經(jīng)驗(yàn)等信息系統(tǒng)測試模型中的參數(shù),是多數(shù)過程模型進(jìn)行參數(shù)本地化的主要方法,如3-PG(Gonzalez-Beneckeetal., 2014; Xieetal., 2017)、BIOME-BGC(何麗鴻等, 2016; 曾慧卿等, 2008)和Tree-BGC(Koroletal., 1995)模型等。M?kel?(1997)基于碳平衡原理開發(fā)的歐洲赤松CROBAS原模型,部分參數(shù)也采用試錯(cuò)法確定,并與林分生長經(jīng)驗(yàn)?zāi)P皖A(yù)測結(jié)果進(jìn)行對(duì)照。試錯(cuò)法相對(duì)簡單易行,但變量過多時(shí)難以同時(shí)降低多個(gè)變量預(yù)測值與實(shí)測值的離差,且預(yù)測值累積誤差會(huì)隨著預(yù)測周期增加而逐漸增大。

    模型預(yù)測準(zhǔn)確性決定模型的應(yīng)用價(jià)值,3-PG和CROBAS等過程模型之所以應(yīng)用較廣,也正是因其預(yù)測準(zhǔn)確性較高(Weiskitteletal., 2011)。森林生物統(tǒng)計(jì)模型的可靠性還常通過不同模型的結(jié)果對(duì)比來評(píng)價(jià)。相對(duì)于Monte Carlo法和貝葉斯法,參數(shù)優(yōu)化方法的優(yōu)點(diǎn)是在處理復(fù)雜模型參數(shù)化時(shí)的收斂速度更快。不同于經(jīng)典的非線性規(guī)劃,種群算法(population-based algorithms)無需設(shè)定決策變量的經(jīng)驗(yàn)初始值,對(duì)目標(biāo)函數(shù)無可導(dǎo)可微的要求,特別是種群算法中的差分演化(differential evolution, DE)算法(Stornetal., 1997),結(jié)構(gòu)簡單、控制參數(shù)少、具有較強(qiáng)的全局收斂性和魯棒性,更適合于解決復(fù)雜環(huán)境中的優(yōu)化問題(盛萬興等, 2017)。目前,差分演化算法已被應(yīng)用于包括森林經(jīng)營規(guī)劃(Pukkala, 2009; Arias-Rodiletal., 2015)在內(nèi)的諸多領(lǐng)域(Lpez Cruzetal., 2003; Onwuboluetal., 2006),但其是否適用于過程模型建模仍有待證明。

    本研究針對(duì)過程模型普遍存在的結(jié)構(gòu)變量復(fù)雜、生理學(xué)數(shù)據(jù)獲取困難的問題,以華山松過程模型CROBAS-PA(PA=Pinusarmandii)的參數(shù)優(yōu)化為例,提出一種結(jié)合經(jīng)驗(yàn)?zāi)P团c過程模型的混合建模方法,探索在建模數(shù)據(jù)不足情況下預(yù)估復(fù)雜過程模型參數(shù)的有效途徑。

    1 研究區(qū)概況

    華山松林分生長數(shù)據(jù)源于秦嶺南坡中段的火地塘林區(qū)。該林區(qū)地處陜西省寧陜縣境內(nèi)(108°21′—108°39′E,33°18′—33°28′N), 海拔1 420~2 474 m,年均氣溫6~13 ℃,年降水量700~1 200 mm,屬北亞熱帶和暖溫帶過渡地帶。土壤主要為火成巖、變質(zhì)巖、石灰?guī)r母質(zhì)上發(fā)育的弱酸性山地棕壤,厚度50~70 cm?;鸬靥亮謪^(qū)華山松林面積占林地面積的22%(李悅黎等, 1993),主要分布在海拔1 600~2 200 m之間。

    2 研究方法

    2.1 數(shù)據(jù)來源

    經(jīng)驗(yàn)?zāi)P徒?shù)據(jù)來源于457塊角規(guī)繞測臨時(shí)樣地(角規(guī)系數(shù)為1)(表1),包括火地塘教學(xué)試驗(yàn)林場1990(244塊)、2005(143塊)和2013(70塊)森林二類調(diào)查樣地。2014年采集該林區(qū)31塊20 m × 20 m固定樣地的每木檢尺數(shù)據(jù)(表1),按地位級(jí)(Ⅰ級(jí)13塊,Ⅱ級(jí)14塊,Ⅲ級(jí)4塊)(吳恒等, 2015)分類,用于參數(shù)優(yōu)化后的過程模型檢驗(yàn)。

    表1 華山松樣地概況Tab.1 Stand attributes of P. armandii plots

    2.2 經(jīng)驗(yàn)?zāi)P图捌鋮?shù)預(yù)估

    本研究建立的華山松經(jīng)驗(yàn)?zāi)P蚎UASSI 1.0系統(tǒng)由式(1)~(6)構(gòu)成。因缺乏優(yōu)勢木樹高數(shù)據(jù),故選取基準(zhǔn)年齡時(shí)的林分平均高,即地位級(jí)指數(shù)SCI(site class index)作為立地質(zhì)量評(píng)價(jià)指標(biāo)。華山松林木材積與樹高連年生長量在林齡30年前已明顯下降(陳存根等, 1994),本研究以30年為基準(zhǔn)年齡。式(1)為立地級(jí)導(dǎo)向曲線,采用Hossfeld方程,已知林分初始年齡和樹高后,可由該式推算SCI并折算不同地位級(jí)指數(shù)下的樹高生長過程(McDilletal., 1992);式(2)為林分密度指數(shù)(Stand density index, SDI)模型(Reineke, 1933),基于林分中立木度最高的部分樣地通過參數(shù)擬合得到;式(3)和(4)是林分胸高斷面積(G, m2·hm-2)模型和平均胸徑(D,cm)模型,屬于Schumacher模型形式(Schumacher, 1939);林分蓄積量(V, m3·hm-2)模型(式6)為華山松材積方程(陳存根等, 1996):

    (1)

    (2)

    (3)

    (4)

    (5)

    (6)

    式中:H為樹高(m);T為林齡(a);a1、a2、a3、b1、b2、b3、b4、c1、c2、c3、c4和β均為系數(shù);N為林分密度(trees·hm-2)。

    本模型構(gòu)建數(shù)據(jù)來自表1描述的樣地,參數(shù)估計(jì)值與檢驗(yàn)見表2。

    為保證檢驗(yàn)樣本獨(dú)立性,按調(diào)查年份將數(shù)據(jù)分為3組(表1),采用交叉驗(yàn)證方法測試模型表現(xiàn),最終以合并數(shù)據(jù)所獲得的模型參數(shù)(表2)模擬林分生長。本研究篩選樣地均為人工起源的華山松相對(duì)純林,林分內(nèi)部個(gè)體年齡差異基本在一個(gè)齡級(jí)范圍內(nèi),按同齡林進(jìn)行建模。式(1)為樹高導(dǎo)向曲線,未考慮除年齡以外引起樹高數(shù)據(jù)變異的因素,導(dǎo)致較低的決定系數(shù)。而胸徑預(yù)測模型(式4)決定系數(shù)較低,是因?yàn)閮H考慮了年齡、立地和林分密度3個(gè)因素,多數(shù)樣地?zé)o法獲取撫育、間伐、擇伐的歷史信息,導(dǎo)致數(shù)據(jù)中較高比例變異未被解釋。

    表2 華山松經(jīng)驗(yàn)?zāi)P蚎UASSI 1.0參數(shù)估計(jì)值Tab.2 Parameter estimates of the empirical model QUASSI 1.0 for P. armandii

    2.3 過程模型

    華山松過程模型CROBAS-PA的原模型為M?kel?(1997)開發(fā)的CROBAS,該模型在模擬單木各器官碳分配過程時(shí)提出3個(gè)假設(shè)(式7—9): 1) 樹冠表面積(Ac, m2)與葉生物量(Wf, kg)異速相關(guān); 2) 葉和細(xì)根生物量(Wr, kg)之比為常數(shù); 3) 葉生物量和冠基處邊材面積之比為常數(shù)(αq, m2·kg-1)(Shinozakietal., 1964)。

    (7)

    Wr=αrWf;

    (8)

    Aq=αqWf,(q=s,b,t)。

    (9)

    式中:Aq為各部分的邊材面積;s、b、t分別代表樹干、樹枝、運(yùn)輸根;ξ′和z為參數(shù);αr為細(xì)根與葉生物量比值。

    模型鏈接了光合作用、呼吸作用、自然整枝、枯損等模塊。華山松過程模型CROBAS-PA以初始葉生物量Wf、初始第一活枝高度Hs、林分初始密度N0為輸入變量,計(jì)算模型的狀態(tài)變量(亦稱模型輸出變量),包括林木各器官在20年間的邊材面積和生物量、冠長、第一活枝高度等。

    本研究通過改變細(xì)根與葉生物量比值αr和細(xì)根衰老速率sr(表3,Sr=srWf,Sr為單樹細(xì)根每年的衰老量,單位為kg·a-1)反映立地條件變化,從而影響林木各器官的碳分配(M?kel?etal., 2016; M?kel?, 1997; Persson, 1980),即以第Ⅱ地位級(jí)參數(shù)為基準(zhǔn)(αr=0.2,sr=1),計(jì)算CROBAS-PA第Ⅰ地位級(jí)(αr=0.1,sr=0.5)和第Ⅲ地位級(jí)(αr=0.25,sr=1.5)的參數(shù)。該方法的依據(jù)是細(xì)根與葉生物比值αr隨立地質(zhì)量提高而變小(Vanninenetal., 1999),即立地較差的林分會(huì)有更高比例的碳分配至細(xì)根,以提高養(yǎng)分水分吸收量。

    2.4 參數(shù)優(yōu)化

    參數(shù)優(yōu)化模型(式10—11)的目標(biāo)函數(shù)設(shè)為過程模型CROBAS-PA與經(jīng)驗(yàn)?zāi)P蚎UASSI 1.0在樹高和生物量預(yù)測上的離差。

    續(xù)表3 Continued

    min{x1,x2,…,x10|Z0}f=∑dRRMSEd;

    (10)

    s.t.X=(x1,x2,…,x10),xj∈[Lj,Uj],

    j=1,2,…,10 。

    (11)

    式中:Z0指林分初始狀態(tài);X=(x1,x2,…,x10)為決策變量x1、x2、…、x10所組成的向量,向量X中的10個(gè)元素是表3中的最后10個(gè)參數(shù);Lj與Uj分別指xj的下界與上界;d為樹高或樹木各部分生物量;RRMSEd指各變量的相對(duì)均方根誤差(relative root mean square error, RRMSE)。

    分別利用經(jīng)驗(yàn)?zāi)P蚎UASSI 1.0與過程模型CROBAS-PA模擬特定初始年齡、立地、密度條件下90塊華山松樣地生長過程,通過優(yōu)化過程模型參數(shù),最小化二者間的差異。樣地具體信息為: 初始年齡20、30、40年,每個(gè)齡組中分別令地位級(jí)指數(shù)SCI=8、9、10、11、12、13 m,林分密度指數(shù)SDI=500、600、700、800、900。將初始年齡、SCI和SDI輸入全林模型后,計(jì)算90塊樣地的林分初始狀態(tài)及后續(xù)20年間的生長過程。20年間的葉、枝和根生物量變化可由CROBAS-PA計(jì)算,并與經(jīng)驗(yàn)生物量公式(陳存根等, 1996)對(duì)比。

    參數(shù)優(yōu)化模型的決策變量選取過程模型CROBAS中基于生理學(xué)和物候?qū)W設(shè)定的10個(gè)參數(shù),即樹葉衰老率、樹冠樹葉分形維數(shù)、比葉面積、最大光合速率、消光系數(shù)、葉表面積密度、自然整枝參數(shù)、樹枝邊材率、樹干邊材衰老率和樹枝邊材衰老率。約束條件為這10個(gè)決策變量的可行域(表3)。

    采用Storn等(1997)提出的差分演化算法(DE),該算法是基于“貪婪競爭”機(jī)制進(jìn)行全局尋優(yōu)的隨機(jī)算法,通過父代變異、雜交的方式產(chǎn)生子代,而后由競爭機(jī)制選擇最優(yōu)解。DE的變異操作充分利用其他個(gè)體信息,增加了種群多樣性,同時(shí)減少了單純個(gè)體變異帶來的隨機(jī)性和盲目性。利用最優(yōu)個(gè)體加隨機(jī)個(gè)體差分法的變異方法,既提高了搜索速度,也降低了陷入局部最優(yōu)的危險(xiǎn);但是若種群分布密度過高,也可能使其陷入局部最優(yōu)。差分演化算法的初始群體從10個(gè)決策變量的可行域中隨機(jī)抽取,10個(gè)決策變量的值和表3中的前26個(gè)參數(shù)一起輸入CROBAS-PA模型進(jìn)行優(yōu)化計(jì)算,從中選取使目標(biāo)函數(shù)最小的一組參數(shù)值。通過以下6個(gè)步驟執(zhí)行優(yōu)化搜索。

    1) 初始化。產(chǎn)生初始種群,種群個(gè)體Xi=(xi1,xi2,…,xijD),其中xij∈[Lj,Uj],i=1,2,…,Np,j=1,2,…,D,D為自變量的維數(shù)。種群大小Np=45,縮放因子F=0.6,雜交概率CR=0.8。迭代次數(shù)Gm=50。

    2) 評(píng)價(jià)初始群體。計(jì)算初始群體中各個(gè)體的適應(yīng)值f(Xi),并記錄最優(yōu)值與最優(yōu)個(gè)體。

    3) 差分變異操作。在種群中隨機(jī)選取不同的個(gè)體Xa、Xc、Xd、Xe、Xf,進(jìn)行變異操作取得新的種群個(gè)體Vi:

    rand/1Vi=Xi+F(Xa-Xc);

    rand/2Vi=Xa+F1(Xc-Xd)+F2(Xe-Xf);

    best/1Vi=Xb+F(Xa-Xc);

    best/2Vi=Xb+F1(Xa-Xc)+F2(Xd-Xe);

    current to best/1Vi=Xi+α(Xc-Xi)+F(Xa-Xc)。

    式中:α為D維向量,其元素是[0,1]之間的隨機(jī)數(shù);F1、F2均為縮放因子(Liuetal., 2010)。

    4) 雜交操作(二項(xiàng)式雜交)。若αj[0,1]

    5) 選擇操作。若f(Ui)

    6) 檢驗(yàn)是否符合條件。如果當(dāng)前的迭代次數(shù)達(dá)到預(yù)定設(shè)置的最大迭代次數(shù)Gm,則停止迭代,輸出最優(yōu)解,否則轉(zhuǎn)到步驟3。

    本研究在優(yōu)化算法選用時(shí),比較了差分演化算法(DE)的6種變異方式(圖1)。其中rand/2變異方式收斂速度最慢,best/1收斂速度最快,但變異方式rand/2與best/1之間的最優(yōu)值竟相差5.8%。通過設(shè)置不同的初始種群,并測試其收斂結(jié)果,最后根據(jù)算法的穩(wěn)定性和收斂速度選用current to best/1變異方式。

    圖1 算法DE的迭代過程Fig.1 Process of iteration of algorithms DE

    2.5 模型檢驗(yàn)及敏感性分析

    為檢驗(yàn)CROBAS-PA模型的準(zhǔn)確性,以31塊固定樣地(表1)作為一組獨(dú)立數(shù)據(jù)對(duì)參數(shù)優(yōu)化后的過程模型進(jìn)行驗(yàn)證。此外,設(shè)計(jì)不同立地(18種)、不同初始密度(5種)的90塊樣地(n=18×5=90),分別采用CROBAS-PA和QUASSI 1.0預(yù)測樹高和胸徑生長過程,而后對(duì)二者在20年間的樹高和胸徑預(yù)測差異進(jìn)行統(tǒng)計(jì)學(xué)驗(yàn)證。使用經(jīng)驗(yàn)?zāi)P蚎UASSI 1.0生成的90塊虛擬樣地可以避免過程模型參數(shù)優(yōu)化因抽樣非均衡帶來的誤差。具體虛擬樣地包括5個(gè)林分密度指數(shù)(SDI=500、600、700、800、900)及3個(gè)地位級(jí)共18個(gè)地位級(jí)指數(shù)(第Ⅰ地位級(jí)SCI=12.3、12.7、13.1、13.5、13.9、14.3 m; 第Ⅱ地位級(jí)SCI=9.9、10.3、10.7、11.1、11.5、11.9 m; 第Ⅲ地位級(jí)SCI=7.7、8.1、8.5、8.9、9.3、9.7 m)。所用的統(tǒng)計(jì)量有平均誤差(mean error, ME)、平均絕對(duì)誤差(mean absolute error, MAE)和平均相對(duì)誤差(mean relative error, MRE),3個(gè)統(tǒng)計(jì)量的絕對(duì)值越小,說明CROBAS-PA預(yù)測的樹高和胸徑越接近QUASSI 1.0 預(yù)測的樹高和胸徑,預(yù)測精度也越高。

    本研究基于生理學(xué)意義選取12個(gè)參數(shù)(樹冠樹葉分形維數(shù)、消光系數(shù)、比葉面積、最大光合速率、樹葉衰老率、葉表面積密度、自然整枝參數(shù)、樹枝邊材率、樹干邊材衰老率、樹枝邊材衰老率、細(xì)根與葉生物量比值和細(xì)根衰老速率)進(jìn)行全局敏感度分析,該方法考慮了單個(gè)參數(shù)及多個(gè)參數(shù)間的相互作用對(duì)模型輸出變量的影響。

    應(yīng)用Sobol法(Cukieretal., 1975)將過程模型全局靈敏度的目標(biāo)函數(shù)設(shè)為CROBAS-PA與QUASSI 1.0在預(yù)測樹高或胸徑變化上的Nash-Sutcliffe效率(NSE)(Nashetal., 1970)。NSE常用于復(fù)雜模型靈敏度分析的目標(biāo)函數(shù)(Yang, 2011; Nossentetal., 2011)。NSE的理論范圍為-∞到1,通常0

    本研究用于靈敏度分析計(jì)算的樣地共36塊(n=3×4×3=36,由QUASSI 1.0模擬),初始年齡分別設(shè)為20、30、40年,地位級(jí)指數(shù)SCI=6.5、8.5、10.5和12.5 m, 林分密度指數(shù)SDI=500、800、1 100。NSE具體計(jì)算公式如下(式12):

    (12)

    3 結(jié)果與分析

    3.1 過程模型參數(shù)本地化

    3.2 過程模型與經(jīng)驗(yàn)?zāi)P捅容^

    過程模型CROBAS-PA與經(jīng)驗(yàn)?zāi)P蚎UASSI 1.0的20~40年樹高生長預(yù)測對(duì)比(圖2)顯示,CROBAS-PA預(yù)測的樹高相較QUASSI 1.0低估0.1~1.0 m,且這種差異在地位級(jí)高的林分中更明顯。在初始密度為2 000和1 500株·hm-2的樣地中,CROBAS-PA預(yù)測的胸徑比QUASSI 1.0低估0~1.1 cm,而在初始密度1 000株·hm-2中,CROBAS-PA預(yù)測的胸徑比QUASSI 1.0高估0.6~1.2 cm(圖3)。

    圖3 基于兩種模型的不同立地與林分初始密度條件下林分胸徑預(yù)測值Fig.3 Comparisons between two models for DBH developments under different site and initial density conditions

    圖4 經(jīng)驗(yàn)生物量公式(陳存根等, 1996)與過程模型CROBAS-PA預(yù)測的葉生物量、枝生物量和根生物量Fig.4 Comparison between empirical biomass equations(Chen et al.,1996) and CROBAS-PA for foliage biomass, branch biomass, and root biomass predictions

    圖5 第Ⅱ地位級(jí)不同初始密度下的林分密度、第一活枝高度和冠長率的比較Fig.5 Comparisons of stand density, crown base height and crown ratio between stands of different initial densities in site classⅡ

    本研究以中等立地與密度條件的林分為例(圖4,第Ⅱ地位級(jí)、初始密度1 500株·hm-2),比較CROBAS-PA和經(jīng)驗(yàn)生物量公式(陳存根等, 1996)預(yù)測的各器官生物量,結(jié)果表明,林齡在20~40年間,CROBAS-PA模型的葉生物量預(yù)測值偏低,枝生物量和根生物量的預(yù)測結(jié)果分別在林齡35年之前和34年之前偏低,但之后偏高。

    經(jīng)驗(yàn)?zāi)P屯ㄟ^林分密度控制來描述競爭,過程模型可通過第一活枝高和冠長率的變化反映林木競爭(圖5)。林分枯損模塊是基于林分密度指數(shù)進(jìn)行構(gòu)建的(Reineke, 1933)即密度高于給定SDI值的林分會(huì)自然枯損。CROBAS-PA基于林冠結(jié)構(gòu)和光截留計(jì)算林冠內(nèi)的蔭蔽環(huán)境和自然整枝過程。密度較大的林分對(duì)空間和光照的競爭更激烈,第一活枝高度較高且冠長率較小(圖5)。

    3.3 預(yù)測誤差與參數(shù)敏感性

    固定樣地實(shí)測值與過程模型CROBAS-PA預(yù)測值的相關(guān)分析表明,CROBAS-PA模型的樹高和胸徑預(yù)測值與實(shí)測值基本呈線性相關(guān)(圖6)。即便31個(gè)用于對(duì)比的固定樣地林分狀況差別較大,如林齡從20到52年,林分密度從289株· hm-2到4 616株· hm-2,但樹高和胸徑的實(shí)測值和模型預(yù)測值的R2仍分別達(dá)到了0.714 8和0.643 3。

    虛擬樣地20年間的生長統(tǒng)計(jì)檢驗(yàn)結(jié)果顯示,CROBAS-PA模型在樹高和胸徑預(yù)測上存在低估,但誤差較小(表4)。CROBAS-PA模型樹高預(yù)測值的平均誤差與平均絕對(duì)誤差均低于1.03 m,平均相對(duì)誤差小于5.59%,其中第Ⅰ地位級(jí)的預(yù)測誤差略高,第Ⅲ地位級(jí)的預(yù)測誤差則較低(表4)。CROBAS-PA胸徑預(yù)測值的平均誤差、平均絕對(duì)誤差和平均相對(duì)誤差分別小于0.36 cm、1.19 cm和2.59%,其中平均相對(duì)誤差最高(2.59%),出現(xiàn)在第Ⅰ地位級(jí); 最高平均誤差(0.36 cm)發(fā)生于第Ⅱ地位級(jí); 而平均絕對(duì)誤差最高(1.19 cm)的為第Ⅲ地位級(jí)(表4)。

    全局靈敏度Sobol分析表明(表5),最大光合速率P、樹冠樹葉分形維數(shù)z、消光系數(shù)k和比葉面積An對(duì)樹高生長影響的總靈敏度系數(shù)排序最高,分別為0.444 3、0.361 0、0.312 9和0.301 0。胸徑的總靈敏度系數(shù)排序最前的4個(gè)參數(shù)與樹高一致。樹高和胸徑對(duì)各參數(shù)的一階靈敏度系數(shù)之和分別為0.646 6、0.588 5,均小于1,說明該模型是不可加的,即非線性的。

    圖6 樹高、胸徑的實(shí)測值和預(yù)測值比較Fig.6 Comparison between observed and simulated values for tree height and DBH

    表4 CROBAS-PA樹高、胸徑預(yù)測值的檢驗(yàn)結(jié)果Tab.4 Model validation of CROBAS-PA on tree height and DBH predictions

    表5 樹高和胸徑的一階靈敏度系數(shù)和總靈敏度系數(shù)Tab.5 Sobol’s first-order indices and total-effect indices for height and DBH

    4 討論

    參數(shù)優(yōu)化后的過程模型CROBAS-PA可以實(shí)現(xiàn)華山松林木各器官生物量的預(yù)估,并預(yù)測不同立地與密度條件下林分胸徑和樹高的生長。靈敏度分析顯示胸徑和樹高生長對(duì)比葉面積An、最大光合速率P、消光系數(shù)k和樹冠樹葉分形維數(shù)z較為敏感。

    CROBAS-PA模型的樹高、胸徑預(yù)測誤差與模型結(jié)構(gòu)簡化直接相關(guān)。Shcherbinina(2012)基于CROBAS預(yù)測短葉松(Pinusbanksiana)生長時(shí),亦出現(xiàn)類似本研究中樹高被低估的結(jié)果,其原因是細(xì)根與葉生物量比值αr的變化對(duì)樹高生長的影響為間接過程(M?kel?, 1997; Shcherbinina, 2012)。另外,本研究顯示CROBAS-PA模型預(yù)測的樹高生長受林分初始狀態(tài)影響較小,但林分初始密度對(duì)CROBAS-PA模型預(yù)測的胸徑精度有一定影響,林分密度越大時(shí)胸徑預(yù)測偏差越大。由于本研究采用異速生長方程(陳存根等, 1996)約束華山松過程模型CROBAS-PA的生物量分配,簡化了林分密度對(duì)葉生物量的影響,而根據(jù)葉生物量計(jì)算胸徑生長,誤差由此而生(圖3)。

    本研究的靈敏度分析證明比葉面積An、最大光合速率P、消光系數(shù)k對(duì)樹高和胸徑生長影響較高(表5),與Ewen(2013)的研究結(jié)果一致。Weiskittel 等(2011)指出過程模型對(duì)光合作用、光截留方面的參數(shù)是非常敏感的。增加An、P、k意味著光合作用增強(qiáng),進(jìn)而促進(jìn)了胸徑和樹高生長。本研究對(duì)葉表面積密度x的靈敏度分析結(jié)果表明,葉表面積密度x對(duì)胸徑和樹高生長影響較低,樹冠樹葉分形維數(shù)z的影響則高于葉表面積密度x(表5),與歐洲赤松的CROBAS模型靈敏度分析結(jié)果類似(M?kel?, 1997)。因?yàn)闃涔跇淙~分形維數(shù)z是異速生長方程中的指數(shù),其不僅影響葉生物量曲線,還影響該曲線的形狀。

    本研究的CROBAS模型是與年齡無關(guān)(age-independent)的林分生長模型,因而可用于模擬復(fù)層或混交林生長動(dòng)態(tài)(Minunnoetal., 2019)。然而,這種與年齡無關(guān)的模型,通常無法描述樹木衰老引起的生理組織代謝速率和碳利用效率的改變。Vanninen等(1996)對(duì)歐洲赤松林木(18~212年生)生物量的研究顯示,細(xì)根與葉生物量在樹木整體生物量的相對(duì)占比隨年齡增加而減小,而與此同時(shí),細(xì)根與葉生物量比值會(huì)隨年齡增大而增大,本研究中參數(shù)細(xì)根與葉生物量比值和細(xì)根衰老速率對(duì)胸徑與樹高預(yù)測的全局靈敏度值相對(duì)較低(表5),因而未設(shè)定其與年齡的關(guān)系,但這種與年齡無關(guān)(age-independent)的模型簡化會(huì)導(dǎo)致對(duì)不同齡級(jí)林木細(xì)根生物量估計(jì)的誤差。因此,在構(gòu)建與年齡相關(guān)(age-dependent)的生長模型時(shí),應(yīng)補(bǔ)充考慮這種變化。

    5 結(jié)論

    本研究通過參數(shù)優(yōu)化模型調(diào)試CROBAS過程模型參數(shù),基于經(jīng)驗(yàn)?zāi)P秃瓦^程模型的混合建模方法構(gòu)建了華山松過程模型CROBAS-PA。研究表明,經(jīng)過生物量異速生長方程和全林經(jīng)驗(yàn)?zāi)P图s束的華山松過程模型CROBAS-PA,可準(zhǔn)確預(yù)測林分樹高和胸徑生長及林木各器官中的生物量分配。全局靈敏度分析顯示,在CROBAS-PA模型參數(shù)中,比葉面積An、最大光合速率P、消光系數(shù)k和樹冠樹葉分形維數(shù)z的敏感度較高。研究證明,在有限實(shí)測數(shù)據(jù)和文獻(xiàn)信息條件下,可在林分水平上通過經(jīng)驗(yàn)?zāi)P图s束保證過程模型預(yù)測的準(zhǔn)確性。這意味著對(duì)于數(shù)據(jù)缺失情境下的復(fù)雜過程模型的參數(shù)預(yù)估,這種經(jīng)驗(yàn)-過程混合建模方法有著特定的優(yōu)勢與應(yīng)用潛力。

    猜你喜歡
    樹高胸徑樣地
    不同造林撫育方式對(duì)木荷林生長的影響
    白城山新1號(hào)楊育苗密度研究
    武漢5種常見園林綠化樹種胸徑與樹高的相關(guān)性研究
    額爾古納市興安落葉松中齡林植被碳儲(chǔ)量研究
    昆明市主要綠化樹種閾值測定與分析
    基于角尺度模型的林業(yè)樣地空間結(jié)構(gòu)分析
    福建省森林資源監(jiān)測體系抽樣調(diào)查中胸徑測量精度范圍的精準(zhǔn)確定
    人工福建柏胸徑與樹高關(guān)系的研究
    不同種源馬尾松樹高與胸徑生長相關(guān)模型研建
    綠色科技(2017年1期)2017-03-01 10:17:01
    不同坡位和造林密度對(duì)四種桉樹生長的影響
    麻豆久久精品国产亚洲av | 黑丝袜美女国产一区| 欧美在线一区亚洲| 国产97色在线日韩免费| 热re99久久国产66热| 亚洲五月色婷婷综合| 久久久国产欧美日韩av| 欧美激情高清一区二区三区| 精品无人区乱码1区二区| 级片在线观看| 成人精品一区二区免费| www.自偷自拍.com| 极品教师在线免费播放| 少妇的丰满在线观看| 国产精品电影一区二区三区| 纯流量卡能插随身wifi吗| 午夜福利,免费看| 一级,二级,三级黄色视频| 成年女人毛片免费观看观看9| 黑人操中国人逼视频| 18禁观看日本| 亚洲av日韩精品久久久久久密| 十八禁人妻一区二区| 国产精品久久电影中文字幕| 又紧又爽又黄一区二区| 宅男免费午夜| 丝袜人妻中文字幕| 黄色a级毛片大全视频| 国产亚洲精品久久久久久毛片| 国产精品香港三级国产av潘金莲| 亚洲国产欧美一区二区综合| www.熟女人妻精品国产| 在线观看66精品国产| 在线观看免费视频网站a站| 日本a在线网址| 久久久国产成人精品二区 | 免费搜索国产男女视频| 久久久久精品国产欧美久久久| 亚洲国产精品合色在线| 色在线成人网| 99国产极品粉嫩在线观看| 纯流量卡能插随身wifi吗| 亚洲午夜精品一区,二区,三区| 精品国产国语对白av| 免费在线观看视频国产中文字幕亚洲| 色老头精品视频在线观看| 身体一侧抽搐| 色在线成人网| 午夜福利,免费看| 亚洲狠狠婷婷综合久久图片| 亚洲黑人精品在线| 麻豆国产av国片精品| 性欧美人与动物交配| 老司机午夜十八禁免费视频| 一级片免费观看大全| 亚洲性夜色夜夜综合| 少妇裸体淫交视频免费看高清 | 国产成人精品无人区| 免费搜索国产男女视频| 久久午夜亚洲精品久久| 麻豆国产av国片精品| 视频区图区小说| 首页视频小说图片口味搜索| www.精华液| 黑人欧美特级aaaaaa片| 精品国产国语对白av| 热99国产精品久久久久久7| 欧美精品亚洲一区二区| 午夜精品国产一区二区电影| 久久人人97超碰香蕉20202| 最好的美女福利视频网| 国产一卡二卡三卡精品| 久久人妻福利社区极品人妻图片| cao死你这个sao货| 国产极品粉嫩免费观看在线| 免费久久久久久久精品成人欧美视频| 黄色毛片三级朝国网站| e午夜精品久久久久久久| 欧美日韩福利视频一区二区| 男女下面插进去视频免费观看| 国产欧美日韩精品亚洲av| 久久久精品欧美日韩精品| 亚洲九九香蕉| 啪啪无遮挡十八禁网站| 久热爱精品视频在线9| 成人三级黄色视频| 亚洲精品国产精品久久久不卡| 大型av网站在线播放| 无人区码免费观看不卡| 亚洲精品国产色婷婷电影| 日日爽夜夜爽网站| 9热在线视频观看99| 男人舔女人的私密视频| 日韩欧美一区视频在线观看| 色综合站精品国产| 色精品久久人妻99蜜桃| aaaaa片日本免费| 亚洲情色 制服丝袜| 在线观看舔阴道视频| 久久精品91蜜桃| 一区福利在线观看| 国产av一区二区精品久久| 久久香蕉国产精品| cao死你这个sao货| 人人妻人人澡人人看| 男女高潮啪啪啪动态图| 免费一级毛片在线播放高清视频 | 99国产精品一区二区三区| 身体一侧抽搐| 一区二区日韩欧美中文字幕| 90打野战视频偷拍视频| 亚洲人成电影观看| 欧美成狂野欧美在线观看| 久久久久国产一级毛片高清牌| 日韩大尺度精品在线看网址 | 亚洲 欧美 日韩 在线 免费| 亚洲国产精品999在线| 亚洲第一青青草原| 久久久久久久午夜电影 | 欧美乱色亚洲激情| 国产无遮挡羞羞视频在线观看| 丁香欧美五月| 色综合站精品国产| 男人的好看免费观看在线视频 | 搡老熟女国产l中国老女人| 超碰97精品在线观看| 在线观看一区二区三区激情| 长腿黑丝高跟| 中亚洲国语对白在线视频| 精品日产1卡2卡| 久久精品国产综合久久久| 亚洲精品在线观看二区| 19禁男女啪啪无遮挡网站| 亚洲精品中文字幕一二三四区| 亚洲精品国产一区二区精华液| 国产精品1区2区在线观看.| 国产欧美日韩一区二区三| 久久九九热精品免费| 啦啦啦在线免费观看视频4| 老熟妇仑乱视频hdxx| 啦啦啦免费观看视频1| 亚洲第一av免费看| 国产又爽黄色视频| 亚洲成人精品中文字幕电影 | 亚洲一卡2卡3卡4卡5卡精品中文| 最新美女视频免费是黄的| 一二三四在线观看免费中文在| 精品一区二区三卡| 亚洲熟妇中文字幕五十中出 | 国产高清视频在线播放一区| 午夜视频精品福利| 午夜精品国产一区二区电影| 色综合婷婷激情| 色婷婷久久久亚洲欧美| 免费在线观看影片大全网站| 久久久国产欧美日韩av| 欧美日韩中文字幕国产精品一区二区三区 | 99热国产这里只有精品6| 性欧美人与动物交配| 亚洲精品国产色婷婷电影| 色老头精品视频在线观看| 搡老岳熟女国产| 极品教师在线免费播放| 国产国语露脸激情在线看| 久久人妻av系列| 国产精品1区2区在线观看.| 美女午夜性视频免费| 精品国产乱码久久久久久男人| 久久久精品国产亚洲av高清涩受| 丁香六月欧美| 久久婷婷成人综合色麻豆| 日韩欧美三级三区| www.www免费av| 男女高潮啪啪啪动态图| 国产精品亚洲一级av第二区| 国产一卡二卡三卡精品| 最近最新免费中文字幕在线| 国产亚洲欧美98| 精品国内亚洲2022精品成人| 老熟妇乱子伦视频在线观看| 国产av在哪里看| 久久久久久大精品| 成人免费观看视频高清| 波多野结衣高清无吗| 黄色女人牲交| 高清黄色对白视频在线免费看| 两人在一起打扑克的视频| 两个人看的免费小视频| 黄色成人免费大全| 国产又色又爽无遮挡免费看| 中文字幕人妻丝袜制服| 色婷婷av一区二区三区视频| 国产成人精品久久二区二区免费| 波多野结衣高清无吗| 亚洲人成77777在线视频| 精品国产一区二区久久| 中文字幕高清在线视频| 动漫黄色视频在线观看| 亚洲成人国产一区在线观看| 嫁个100分男人电影在线观看| 国产人伦9x9x在线观看| 国产亚洲av高清不卡| 宅男免费午夜| 欧美成狂野欧美在线观看| 久久久久久久久中文| 国产不卡一卡二| 午夜精品久久久久久毛片777| 757午夜福利合集在线观看| 男女床上黄色一级片免费看| 欧洲精品卡2卡3卡4卡5卡区| 韩国精品一区二区三区| 如日韩欧美国产精品一区二区三区| 国产精品爽爽va在线观看网站 | 日韩欧美一区二区三区在线观看| 国产亚洲欧美98| 日本黄色视频三级网站网址| 大型av网站在线播放| 一级毛片精品| 一进一出好大好爽视频| 日本撒尿小便嘘嘘汇集6| 国产亚洲av高清不卡| 俄罗斯特黄特色一大片| 黄色毛片三级朝国网站| 另类亚洲欧美激情| 久久草成人影院| 色哟哟哟哟哟哟| 亚洲国产精品一区二区三区在线| 国产精品野战在线观看 | 亚洲色图综合在线观看| 少妇被粗大的猛进出69影院| 亚洲欧美精品综合久久99| 丝袜在线中文字幕| 亚洲黑人精品在线| 老司机深夜福利视频在线观看| 国产欧美日韩精品亚洲av| 校园春色视频在线观看| 男女午夜视频在线观看| 中文字幕人妻熟女乱码| 丰满饥渴人妻一区二区三| 黄频高清免费视频| 黑人操中国人逼视频| 日韩欧美三级三区| 男女高潮啪啪啪动态图| 久久国产乱子伦精品免费另类| 欧美乱妇无乱码| www日本在线高清视频| 欧美激情极品国产一区二区三区| 两个人看的免费小视频| 黄色女人牲交| 亚洲avbb在线观看| 在线观看一区二区三区激情| 桃色一区二区三区在线观看| 国产深夜福利视频在线观看| 一级作爱视频免费观看| 热re99久久精品国产66热6| 大型av网站在线播放| 久久久久久久久中文| 免费在线观看完整版高清| 精品欧美一区二区三区在线| 久久人人精品亚洲av| 欧美人与性动交α欧美精品济南到| 美女国产高潮福利片在线看| 久久中文字幕一级| 精品高清国产在线一区| 免费高清视频大片| √禁漫天堂资源中文www| 日本五十路高清| 真人一进一出gif抽搐免费| av免费在线观看网站| 大陆偷拍与自拍| 脱女人内裤的视频| 男女做爰动态图高潮gif福利片 | 国产激情欧美一区二区| 国产精品一区二区三区四区久久 | 国产黄色免费在线视频| 九色亚洲精品在线播放| 大型av网站在线播放| 国产亚洲精品综合一区在线观看 | 老司机深夜福利视频在线观看| 美女高潮到喷水免费观看| 最近最新中文字幕大全电影3 | 曰老女人黄片| 亚洲精品中文字幕一二三四区| 亚洲中文日韩欧美视频| 国产黄a三级三级三级人| 在线播放国产精品三级| 搡老岳熟女国产| 很黄的视频免费| 色哟哟哟哟哟哟| 视频在线观看一区二区三区| 国产精品99久久99久久久不卡| 中文字幕另类日韩欧美亚洲嫩草| 好看av亚洲va欧美ⅴa在| 亚洲熟妇熟女久久| 我的亚洲天堂| 男女床上黄色一级片免费看| 精品第一国产精品| 国产高清激情床上av| 十分钟在线观看高清视频www| 国产精品综合久久久久久久免费 | 女人精品久久久久毛片| 最好的美女福利视频网| 操美女的视频在线观看| 嫩草影视91久久| 一区在线观看完整版| 亚洲欧美日韩另类电影网站| 视频在线观看一区二区三区| 麻豆成人av在线观看| 国产亚洲精品久久久久5区| 精品高清国产在线一区| 91成年电影在线观看| 亚洲色图综合在线观看| 久久香蕉国产精品| 真人做人爱边吃奶动态| 91成年电影在线观看| 国产激情欧美一区二区| 国产一区二区在线av高清观看| 精品久久久久久久久久免费视频 | 欧美丝袜亚洲另类 | 亚洲av电影在线进入| 中文字幕av电影在线播放| 天堂动漫精品| 亚洲在线自拍视频| 免费高清视频大片| 国产在线精品亚洲第一网站| 涩涩av久久男人的天堂| 一边摸一边抽搐一进一小说| 亚洲精品国产一区二区精华液| www国产在线视频色| 99国产精品99久久久久| 久久国产精品男人的天堂亚洲| 麻豆成人av在线观看| 男女做爰动态图高潮gif福利片 | 80岁老熟妇乱子伦牲交| 一进一出抽搐动态| 亚洲av五月六月丁香网| 亚洲精品粉嫩美女一区| 成熟少妇高潮喷水视频| 麻豆久久精品国产亚洲av | 精品卡一卡二卡四卡免费| av视频免费观看在线观看| 精品福利观看| 波多野结衣高清无吗| 色精品久久人妻99蜜桃| 欧美激情高清一区二区三区| 亚洲自偷自拍图片 自拍| 久久人妻福利社区极品人妻图片| 中文字幕精品免费在线观看视频| 久久人妻熟女aⅴ| 国产亚洲精品一区二区www| 亚洲av电影在线进入| 精品高清国产在线一区| 黑丝袜美女国产一区| 亚洲国产欧美网| 久久久久久大精品| 精品高清国产在线一区| 丝袜美腿诱惑在线| 欧美乱色亚洲激情| 99精品欧美一区二区三区四区| 丰满人妻熟妇乱又伦精品不卡| 男女下面进入的视频免费午夜 | 交换朋友夫妻互换小说| 亚洲一区中文字幕在线| 午夜福利在线观看吧| 午夜激情av网站| 99国产精品一区二区蜜桃av| 一级a爱视频在线免费观看| 色精品久久人妻99蜜桃| 黄网站色视频无遮挡免费观看| 日韩高清综合在线| 亚洲精品在线美女| 精品熟女少妇八av免费久了| 老汉色∧v一级毛片| 久久午夜亚洲精品久久| 亚洲精品在线美女| 99热国产这里只有精品6| 18禁观看日本| 亚洲免费av在线视频| 久久人人精品亚洲av| 首页视频小说图片口味搜索| 国产精品香港三级国产av潘金莲| 高清在线国产一区| 男女床上黄色一级片免费看| 国产1区2区3区精品| 在线天堂中文资源库| 久久 成人 亚洲| 久久久国产精品麻豆| 黑人操中国人逼视频| 国内毛片毛片毛片毛片毛片| 亚洲一区中文字幕在线| 99riav亚洲国产免费| 丰满饥渴人妻一区二区三| tocl精华| 999久久久国产精品视频| 日本wwww免费看| 国产亚洲精品第一综合不卡| 婷婷丁香在线五月| 一进一出好大好爽视频| 免费观看精品视频网站| 91成年电影在线观看| 国产亚洲精品一区二区www| 天堂中文最新版在线下载| 岛国在线观看网站| 日本a在线网址| 亚洲精品久久午夜乱码| 久久久国产欧美日韩av| 美女大奶头视频| 久久精品亚洲av国产电影网| 亚洲精品成人av观看孕妇| 亚洲aⅴ乱码一区二区在线播放 | 18禁观看日本| 国产av又大| 午夜两性在线视频| 这个男人来自地球电影免费观看| 欧美不卡视频在线免费观看 | 久久久久九九精品影院| 老司机深夜福利视频在线观看| av有码第一页| 久久热在线av| 男女下面进入的视频免费午夜 | 一二三四在线观看免费中文在| 变态另类成人亚洲欧美熟女 | 成人三级黄色视频| 亚洲第一av免费看| av在线天堂中文字幕 | 97人妻天天添夜夜摸| 精品少妇一区二区三区视频日本电影| 在线观看免费视频网站a站| 极品人妻少妇av视频| 高清欧美精品videossex| 午夜影院日韩av| 这个男人来自地球电影免费观看| 夜夜躁狠狠躁天天躁| 国产真人三级小视频在线观看| 青草久久国产| av免费在线观看网站| 国产成人精品无人区| 欧美久久黑人一区二区| 亚洲 欧美 日韩 在线 免费| 久久精品人人爽人人爽视色| 人人妻人人添人人爽欧美一区卜| 精品一区二区三卡| 精品无人区乱码1区二区| 男人的好看免费观看在线视频 | 国产av在哪里看| 精品熟女少妇八av免费久了| 天堂俺去俺来也www色官网| 亚洲国产欧美一区二区综合| 99久久久亚洲精品蜜臀av| 美女高潮喷水抽搐中文字幕| 国产av一区在线观看免费| 麻豆久久精品国产亚洲av | 这个男人来自地球电影免费观看| 看黄色毛片网站| 99国产精品免费福利视频| 日韩大码丰满熟妇| 亚洲欧美激情综合另类| 欧美日韩黄片免| 亚洲av第一区精品v没综合| 午夜91福利影院| 19禁男女啪啪无遮挡网站| 久久精品影院6| 99国产精品一区二区蜜桃av| 岛国在线观看网站| 香蕉丝袜av| 不卡av一区二区三区| 国产成年人精品一区二区 | 在线视频色国产色| x7x7x7水蜜桃| 三上悠亚av全集在线观看| 久久久久国内视频| 中文字幕精品免费在线观看视频| 少妇裸体淫交视频免费看高清 | 午夜福利影视在线免费观看| 黑丝袜美女国产一区| netflix在线观看网站| 久久久精品欧美日韩精品| 国产伦一二天堂av在线观看| 亚洲九九香蕉| 亚洲专区字幕在线| 日本a在线网址| 岛国在线观看网站| 亚洲国产欧美一区二区综合| 欧美人与性动交α欧美软件| 久久人妻熟女aⅴ| 女性被躁到高潮视频| av超薄肉色丝袜交足视频| 国产精品98久久久久久宅男小说| av超薄肉色丝袜交足视频| 久久香蕉精品热| 国产成人免费无遮挡视频| 一级a爱片免费观看的视频| 男女做爰动态图高潮gif福利片 | 欧美日韩中文字幕国产精品一区二区三区 | 亚洲一卡2卡3卡4卡5卡精品中文| 丝袜在线中文字幕| tocl精华| 亚洲精品中文字幕在线视频| 91国产中文字幕| 亚洲av五月六月丁香网| 51午夜福利影视在线观看| 看黄色毛片网站| 在线永久观看黄色视频| 欧美成人免费av一区二区三区| 国产91精品成人一区二区三区| 久久久久久久久中文| 免费女性裸体啪啪无遮挡网站| 一本综合久久免费| 日韩欧美免费精品| 极品人妻少妇av视频| 黄色怎么调成土黄色| 午夜福利,免费看| 国产精品 欧美亚洲| 日韩 欧美 亚洲 中文字幕| 亚洲成a人片在线一区二区| 亚洲av成人不卡在线观看播放网| 午夜老司机福利片| 99久久人妻综合| 女人精品久久久久毛片| 日韩成人在线观看一区二区三区| 十分钟在线观看高清视频www| 精品久久久久久,| 精品欧美一区二区三区在线| 精品国产一区二区三区四区第35| 国产区一区二久久| 在线视频色国产色| 久久久久国产一级毛片高清牌| 巨乳人妻的诱惑在线观看| 国产亚洲欧美精品永久| av欧美777| 亚洲欧洲精品一区二区精品久久久| 91麻豆精品激情在线观看国产 | 91成人精品电影| 性欧美人与动物交配| 一a级毛片在线观看| 国产精品永久免费网站| a级片在线免费高清观看视频| 最近最新中文字幕大全电影3 | 免费看a级黄色片| 成人影院久久| 欧美激情高清一区二区三区| 国产无遮挡羞羞视频在线观看| 国产单亲对白刺激| 天堂√8在线中文| 女性被躁到高潮视频| 精品熟女少妇八av免费久了| 国产高清激情床上av| 丝袜美足系列| 日韩高清综合在线| 亚洲成人精品中文字幕电影 | 丰满的人妻完整版| 女人爽到高潮嗷嗷叫在线视频| a在线观看视频网站| 国产三级在线视频| 丰满迷人的少妇在线观看| 少妇裸体淫交视频免费看高清 | 黄片大片在线免费观看| 久久人人精品亚洲av| 黄色a级毛片大全视频| 高清黄色对白视频在线免费看| 午夜精品国产一区二区电影| 亚洲精品国产一区二区精华液| 欧美日韩中文字幕国产精品一区二区三区 | 在线国产一区二区在线| 一边摸一边抽搐一进一小说| 国产av在哪里看| 久久久久久久久免费视频了| 亚洲av熟女| 三上悠亚av全集在线观看| 国产精品 欧美亚洲| 99精国产麻豆久久婷婷| 精品国产超薄肉色丝袜足j| 自拍欧美九色日韩亚洲蝌蚪91| 国产欧美日韩精品亚洲av| 亚洲av片天天在线观看| 性色av乱码一区二区三区2| 国产精品1区2区在线观看.| 国产精品成人在线| 丰满人妻熟妇乱又伦精品不卡| 亚洲伊人色综图| 亚洲第一欧美日韩一区二区三区| 一级a爱片免费观看的视频| 国产成人啪精品午夜网站| 搡老岳熟女国产| 久久中文字幕人妻熟女| 午夜福利影视在线免费观看| 日韩av在线大香蕉| 久久香蕉精品热| 黄色片一级片一级黄色片| 欧美日本亚洲视频在线播放| 女人精品久久久久毛片| 男女午夜视频在线观看| 国产精品电影一区二区三区| 在线天堂中文资源库| 欧美乱码精品一区二区三区| 国产成人精品无人区| 侵犯人妻中文字幕一二三四区| 两性夫妻黄色片| 亚洲少妇的诱惑av| 一级片'在线观看视频| 欧美激情 高清一区二区三区| 亚洲午夜理论影院| 国产精品久久久久久人妻精品电影| 亚洲欧美激情在线| 99热国产这里只有精品6| 日本欧美视频一区| 两个人看的免费小视频| 日韩视频一区二区在线观看| 女性被躁到高潮视频| 久久九九热精品免费| 日本五十路高清| 高清在线国产一区| 视频区欧美日本亚洲| 曰老女人黄片| 亚洲一卡2卡3卡4卡5卡精品中文|