• <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ì)四種桉樹生長的影響
    一进一出抽搐动态| 国产91精品成人一区二区三区| 一级毛片精品| 亚洲国产色片| 制服丝袜大香蕉在线| 美女大奶头视频| 九九在线视频观看精品| 人妻丰满熟妇av一区二区三区| 久9热在线精品视频| 欧美在线一区亚洲| 长腿黑丝高跟| 一级a爱片免费观看的视频| 国内精品美女久久久久久| 韩国av一区二区三区四区| 亚洲一区高清亚洲精品| 精品欧美国产一区二区三| 国内毛片毛片毛片毛片毛片| www.999成人在线观看| 搡老岳熟女国产| 欧美一区二区精品小视频在线| 国产成人av教育| 亚洲精品美女久久av网站| 老汉色av国产亚洲站长工具| 男女午夜视频在线观看| 精品一区二区三区视频在线观看免费| 色老头精品视频在线观看| 日韩大尺度精品在线看网址| 最新美女视频免费是黄的| 后天国语完整版免费观看| a级毛片a级免费在线| 国产一区二区在线av高清观看| 日本精品一区二区三区蜜桃| 好男人电影高清在线观看| 国产精品综合久久久久久久免费| 国产极品精品免费视频能看的| 国产精品98久久久久久宅男小说| 亚洲专区国产一区二区| 老司机深夜福利视频在线观看| 黑人欧美特级aaaaaa片| 国产综合懂色| 国产成人一区二区三区免费视频网站| 午夜精品在线福利| 两性夫妻黄色片| 少妇裸体淫交视频免费看高清| 亚洲av成人av| 欧美黄色片欧美黄色片| av天堂在线播放| 动漫黄色视频在线观看| 国产熟女xx| av片东京热男人的天堂| 欧美色视频一区免费| 精品99又大又爽又粗少妇毛片 | 国产在线精品亚洲第一网站| 一卡2卡三卡四卡精品乱码亚洲| 在线免费观看不下载黄p国产 | 国产私拍福利视频在线观看| av天堂在线播放| 嫩草影视91久久| 精品国产亚洲在线| 白带黄色成豆腐渣| 亚洲av五月六月丁香网| 国产黄片美女视频| 国产精品亚洲美女久久久| 精品久久蜜臀av无| av片东京热男人的天堂| 中出人妻视频一区二区| 国产69精品久久久久777片 | 99热这里只有是精品50| 中文字幕久久专区| 亚洲 欧美一区二区三区| 熟女少妇亚洲综合色aaa.| 老司机午夜十八禁免费视频| 免费在线观看成人毛片| 99热6这里只有精品| 中文字幕精品亚洲无线码一区| 1024手机看黄色片| 禁无遮挡网站| 久久欧美精品欧美久久欧美| 日日摸夜夜添夜夜添小说| 无遮挡黄片免费观看| 97碰自拍视频| 日日夜夜操网爽| 高潮久久久久久久久久久不卡| 99热这里只有精品一区 | 小蜜桃在线观看免费完整版高清| 看片在线看免费视频| 天堂av国产一区二区熟女人妻| 国产 一区 欧美 日韩| 级片在线观看| 欧美日韩亚洲国产一区二区在线观看| 丰满人妻熟妇乱又伦精品不卡| 99久国产av精品| 99久久精品国产亚洲精品| 别揉我奶头~嗯~啊~动态视频| 欧美日韩一级在线毛片| 母亲3免费完整高清在线观看| 天天躁日日操中文字幕| 国产精品久久视频播放| 国产成人福利小说| 色综合亚洲欧美另类图片| 丰满人妻一区二区三区视频av | 麻豆一二三区av精品| 色综合婷婷激情| 久久久久久久久中文| 亚洲中文日韩欧美视频| 一区二区三区国产精品乱码| 亚洲av片天天在线观看| 91在线观看av| 黑人巨大精品欧美一区二区mp4| 亚洲专区国产一区二区| 男女之事视频高清在线观看| 琪琪午夜伦伦电影理论片6080| 亚洲av日韩精品久久久久久密| 日韩欧美三级三区| 国产91精品成人一区二区三区| 一区福利在线观看| 国产精品久久久久久精品电影| 99re在线观看精品视频| 日本与韩国留学比较| 国产不卡一卡二| 91麻豆精品激情在线观看国产| 丰满的人妻完整版| 美女午夜性视频免费| 午夜免费观看网址| 久久久久久久久久黄片| 亚洲国产欧美人成| 黄色视频,在线免费观看| 夜夜爽天天搞| 亚洲精品久久国产高清桃花| 久久久国产精品麻豆| 久久精品国产清高在天天线| 久久亚洲精品不卡| 亚洲国产看品久久| 一级a爱片免费观看的视频| 18禁黄网站禁片免费观看直播| 草草在线视频免费看| 18禁国产床啪视频网站| 国产av不卡久久| 嫩草影院精品99| 男女之事视频高清在线观看| 国产精品久久视频播放| 免费观看的影片在线观看| 亚洲欧美日韩高清专用| bbb黄色大片| 国产成人一区二区三区免费视频网站| 成人亚洲精品av一区二区| 精品一区二区三区av网在线观看| 一区二区三区国产精品乱码| 久久精品国产综合久久久| 一个人看视频在线观看www免费 | 两个人看的免费小视频| 丁香六月欧美| 国产精品一区二区精品视频观看| 亚洲专区国产一区二区| 两个人的视频大全免费| 一进一出抽搐gif免费好疼| 午夜福利在线在线| 亚洲国产看品久久| 床上黄色一级片| 国内毛片毛片毛片毛片毛片| 极品教师在线免费播放| 久久精品91无色码中文字幕| 国产成人av教育| 亚洲国产中文字幕在线视频| 亚洲欧美日韩高清在线视频| 搡老熟女国产l中国老女人| 国内精品一区二区在线观看| 国产极品精品免费视频能看的| 日本黄大片高清| 无遮挡黄片免费观看| 最新中文字幕久久久久 | 亚洲av电影在线进入| а√天堂www在线а√下载| 毛片女人毛片| 欧美丝袜亚洲另类 | 精品一区二区三区四区五区乱码| 国产精品亚洲av一区麻豆| 色av中文字幕| 99久久精品一区二区三区| 老司机深夜福利视频在线观看| 麻豆成人av在线观看| 亚洲专区中文字幕在线| 精品久久久久久,| 国产欧美日韩一区二区三| 宅男免费午夜| 嫩草影院精品99| 18美女黄网站色大片免费观看| 国产淫片久久久久久久久 | 一级毛片女人18水好多| 99久久99久久久精品蜜桃| 岛国在线观看网站| 亚洲国产高清在线一区二区三| 国产三级在线视频| 搡老岳熟女国产| 美女高潮喷水抽搐中文字幕| 琪琪午夜伦伦电影理论片6080| 99国产精品一区二区蜜桃av| 最近在线观看免费完整版| 国产单亲对白刺激| 国产综合懂色| 中文字幕av在线有码专区| 亚洲午夜精品一区,二区,三区| 全区人妻精品视频| 波多野结衣高清无吗| 国产真人三级小视频在线观看| 国产私拍福利视频在线观看| 日韩av在线大香蕉| 99在线人妻在线中文字幕| 露出奶头的视频| 欧美黑人巨大hd| 中文字幕精品亚洲无线码一区| 成人三级做爰电影| 精品国内亚洲2022精品成人| 在线国产一区二区在线| 国产精品国产高清国产av| 亚洲精品在线观看二区| 精品乱码久久久久久99久播| 69av精品久久久久久| av女优亚洲男人天堂 | 欧美又色又爽又黄视频| 国产99白浆流出| 人人妻人人澡欧美一区二区| 国产91精品成人一区二区三区| 欧美一区二区国产精品久久精品| 日韩高清综合在线| e午夜精品久久久久久久| 国产精品永久免费网站| 成人精品一区二区免费| 日本成人三级电影网站| 亚洲精品在线观看二区| 91麻豆精品激情在线观看国产| 制服丝袜大香蕉在线| 老司机午夜十八禁免费视频| 久久久久久久久免费视频了| 757午夜福利合集在线观看| 欧美激情在线99| 美女cb高潮喷水在线观看 | 黑人巨大精品欧美一区二区mp4| 波多野结衣高清无吗| 亚洲av成人av| 18美女黄网站色大片免费观看| 美女高潮的动态| 脱女人内裤的视频| 757午夜福利合集在线观看| 97碰自拍视频| 国产精品一区二区三区四区久久| 亚洲五月天丁香| 久久久色成人| 最近最新中文字幕大全电影3| 国内精品一区二区在线观看| 精品人妻1区二区| 亚洲av免费在线观看| 午夜亚洲福利在线播放| 国产精品美女特级片免费视频播放器 | 午夜精品在线福利| 国产久久久一区二区三区| 国产免费男女视频| 大型黄色视频在线免费观看| 日本精品一区二区三区蜜桃| 午夜亚洲福利在线播放| 精品熟女少妇八av免费久了| 国产又黄又爽又无遮挡在线| 不卡av一区二区三区| 国产精品综合久久久久久久免费| 岛国视频午夜一区免费看| 天天躁日日操中文字幕| 亚洲av电影不卡..在线观看| 熟妇人妻久久中文字幕3abv| 午夜免费成人在线视频| 麻豆国产97在线/欧美| 男人的好看免费观看在线视频| 亚洲国产精品久久男人天堂| 成人高潮视频无遮挡免费网站| 伊人久久大香线蕉亚洲五| 男女床上黄色一级片免费看| 国产精品久久久久久人妻精品电影| 性欧美人与动物交配| 9191精品国产免费久久| ponron亚洲| 国产精品av视频在线免费观看| 亚洲熟妇熟女久久| 精品欧美国产一区二区三| 精品久久久久久久久久久久久| 少妇熟女aⅴ在线视频| 国产男靠女视频免费网站| 88av欧美| 日本 欧美在线| 亚洲av电影不卡..在线观看| 日韩欧美在线乱码| 大型黄色视频在线免费观看| 午夜福利免费观看在线| 欧美乱妇无乱码| 久久国产精品人妻蜜桃| 国内精品一区二区在线观看| 97人妻精品一区二区三区麻豆| av在线天堂中文字幕| 三级毛片av免费| 免费一级毛片在线播放高清视频| 麻豆久久精品国产亚洲av| 久久婷婷人人爽人人干人人爱| 麻豆国产97在线/欧美| 欧美色视频一区免费| 精品国产乱子伦一区二区三区| 成人国产综合亚洲| 日韩 欧美 亚洲 中文字幕| 51午夜福利影视在线观看| 欧美成人免费av一区二区三区| www日本在线高清视频| 欧美av亚洲av综合av国产av| 国产精品日韩av在线免费观看| 国产真实乱freesex| 成人一区二区视频在线观看| 我的老师免费观看完整版| 国产高清视频在线观看网站| 高清毛片免费观看视频网站| 国产69精品久久久久777片 | 亚洲熟妇中文字幕五十中出| 一本久久中文字幕| 男女那种视频在线观看| 一级黄色大片毛片| 99久久精品热视频| 99riav亚洲国产免费| 真实男女啪啪啪动态图| 夜夜看夜夜爽夜夜摸| 亚洲精品国产精品久久久不卡| 最近最新中文字幕大全免费视频| 成人鲁丝片一二三区免费| 国产精品野战在线观看| 久久久久久九九精品二区国产| 国产一区二区激情短视频| 麻豆av在线久日| 伦理电影免费视频| 国产精品亚洲一级av第二区| 香蕉av资源在线| 丝袜人妻中文字幕| 国产成人av激情在线播放| 国产一级毛片七仙女欲春2| 亚洲精品色激情综合| 真人一进一出gif抽搐免费| 可以在线观看毛片的网站| 999久久久精品免费观看国产| 亚洲av免费在线观看| 国产精华一区二区三区| 久久久久久久久免费视频了| 日韩精品青青久久久久久| 天天添夜夜摸| 又爽又黄无遮挡网站| 少妇熟女aⅴ在线视频| 国内揄拍国产精品人妻在线| 香蕉国产在线看| 午夜日韩欧美国产| 国产成人系列免费观看| 国模一区二区三区四区视频 | 午夜影院日韩av| 一夜夜www| 精品久久久久久成人av| 999久久久精品免费观看国产| 国产黄a三级三级三级人| 高清在线国产一区| 亚洲av中文字字幕乱码综合| 男人舔奶头视频| 成人特级av手机在线观看| 美女 人体艺术 gogo| 特级一级黄色大片| 一级a爱片免费观看的视频| www.www免费av| 日韩大尺度精品在线看网址| 身体一侧抽搐| 琪琪午夜伦伦电影理论片6080| 日本黄色片子视频| 婷婷亚洲欧美| 日本 av在线| 给我免费播放毛片高清在线观看| 国产1区2区3区精品| 男女午夜视频在线观看| 欧美黑人巨大hd| 91老司机精品| 国产三级中文精品| 国产精品 国内视频| 制服人妻中文乱码| 欧美三级亚洲精品| 久久久水蜜桃国产精品网| 中文字幕精品亚洲无线码一区| 国产成+人综合+亚洲专区| 久久精品aⅴ一区二区三区四区| 免费观看的影片在线观看| 亚洲午夜理论影院| 久久这里只有精品中国| 欧美成狂野欧美在线观看| 搡老熟女国产l中国老女人| 国产精品av视频在线免费观看| 18禁裸乳无遮挡免费网站照片| 啦啦啦观看免费观看视频高清| 日韩免费av在线播放| 亚洲美女视频黄频| 一本一本综合久久| 一区二区三区国产精品乱码| 五月伊人婷婷丁香| 亚洲aⅴ乱码一区二区在线播放| 真实男女啪啪啪动态图| 婷婷精品国产亚洲av在线| 在线观看美女被高潮喷水网站 | 国产99白浆流出| 91av网站免费观看| 日韩欧美国产在线观看| 国产午夜福利久久久久久| 国产欧美日韩精品亚洲av| 婷婷丁香在线五月| 动漫黄色视频在线观看| 欧美性猛交黑人性爽| a在线观看视频网站| 久久久精品大字幕| 午夜福利18| 午夜福利视频1000在线观看| 久久久久免费精品人妻一区二区| 国产欧美日韩一区二区三| 久久久久性生活片| 麻豆成人av在线观看| 国产精品爽爽va在线观看网站| 99久久综合精品五月天人人| 成年女人看的毛片在线观看| 亚洲国产欧美人成| 在线播放国产精品三级| 99久久久亚洲精品蜜臀av| 男女之事视频高清在线观看| 国产 一区 欧美 日韩| 久久久久久大精品| 国产淫片久久久久久久久 | 国产主播在线观看一区二区| 真人一进一出gif抽搐免费| 一本综合久久免费| 在线免费观看的www视频| 国产精品亚洲av一区麻豆| 噜噜噜噜噜久久久久久91| x7x7x7水蜜桃| 色老头精品视频在线观看| 无遮挡黄片免费观看| 亚洲av日韩精品久久久久久密| 精品免费久久久久久久清纯| 国产成人系列免费观看| 国产伦一二天堂av在线观看| 亚洲精华国产精华精| 舔av片在线| 午夜免费激情av| 99国产精品99久久久久| 欧美在线黄色| 好男人在线观看高清免费视频| 国语自产精品视频在线第100页| 国产av不卡久久| 午夜成年电影在线免费观看| 99久久久亚洲精品蜜臀av| 亚洲五月天丁香| 欧美日本视频| 免费高清视频大片| 欧美3d第一页| 夜夜看夜夜爽夜夜摸| 国产黄色小视频在线观看| 99精品在免费线老司机午夜| 精品国产美女av久久久久小说| 成年免费大片在线观看| www.999成人在线观看| 午夜精品久久久久久毛片777| 欧美日韩瑟瑟在线播放| 久久久水蜜桃国产精品网| 成人18禁在线播放| 日韩欧美国产一区二区入口| 中国美女看黄片| 深夜精品福利| 国产成人欧美在线观看| 熟女电影av网| 中文字幕熟女人妻在线| 国产探花在线观看一区二区| avwww免费| 两人在一起打扑克的视频| 黄色 视频免费看| a级毛片在线看网站| 亚洲一区二区三区不卡视频| 国产成年人精品一区二区| 人妻久久中文字幕网| 亚洲成人久久性| 一本久久中文字幕| 亚洲国产精品合色在线| 成人特级av手机在线观看| 九色国产91popny在线| 草草在线视频免费看| 夜夜夜夜夜久久久久| 日本免费一区二区三区高清不卡| 午夜福利在线观看吧| 亚洲最大成人中文| 90打野战视频偷拍视频| 亚洲在线自拍视频| 美女午夜性视频免费| 午夜激情福利司机影院| 99久久久亚洲精品蜜臀av| 免费观看的影片在线观看| 麻豆成人av在线观看| 91在线精品国自产拍蜜月 | 国产高清videossex| 亚洲成人久久性| 十八禁人妻一区二区| 激情在线观看视频在线高清| 国产伦人伦偷精品视频| 99久久精品一区二区三区| 欧美黄色片欧美黄色片| www国产在线视频色| 日韩欧美国产一区二区入口| 真人一进一出gif抽搐免费| 久久国产精品影院| 嫩草影院入口| 少妇熟女aⅴ在线视频| 久久中文字幕人妻熟女| 日本一本二区三区精品| 欧美午夜高清在线| 法律面前人人平等表现在哪些方面| 亚洲精品在线观看二区| 18禁黄网站禁片午夜丰满| 成年女人看的毛片在线观看| 天天躁日日操中文字幕| 超碰成人久久| 久久中文字幕人妻熟女| 一级毛片高清免费大全| 国产精品一区二区三区四区久久| 国产又黄又爽又无遮挡在线| 色噜噜av男人的天堂激情| 日本撒尿小便嘘嘘汇集6| 每晚都被弄得嗷嗷叫到高潮| 国产亚洲av嫩草精品影院| 国产极品精品免费视频能看的| 国产精品 国内视频| 精品福利观看| 最近最新免费中文字幕在线| 国产精品一区二区三区四区免费观看 | 香蕉国产在线看| av福利片在线观看| 中亚洲国语对白在线视频| 国产黄片美女视频| 亚洲黑人精品在线| 两性午夜刺激爽爽歪歪视频在线观看| 国产精品爽爽va在线观看网站| 丁香六月欧美| 长腿黑丝高跟| 精品一区二区三区av网在线观看| 美女黄网站色视频| 精品久久久久久成人av| 巨乳人妻的诱惑在线观看| 成年免费大片在线观看| 村上凉子中文字幕在线| 国产69精品久久久久777片 | 成人性生交大片免费视频hd| 亚洲av五月六月丁香网| 亚洲av第一区精品v没综合| 真人做人爱边吃奶动态| 亚洲av美国av| 九色国产91popny在线| 午夜免费观看网址| 久久午夜综合久久蜜桃| 在线永久观看黄色视频| 国产av一区在线观看免费| 黄色女人牲交| 级片在线观看| 国产精品,欧美在线| av黄色大香蕉| 亚洲自拍偷在线| av片东京热男人的天堂| 欧美+亚洲+日韩+国产| 三级毛片av免费| 国产三级中文精品| 国产亚洲欧美98| 亚洲国产精品999在线| 99久久综合精品五月天人人| 啦啦啦免费观看视频1| 熟女少妇亚洲综合色aaa.| 最新中文字幕久久久久 | 亚洲最大成人中文| 黄频高清免费视频| 国产精品99久久久久久久久| 99国产综合亚洲精品| 韩国av一区二区三区四区| 免费看美女性在线毛片视频| 91九色精品人成在线观看| 露出奶头的视频| 亚洲国产欧美网| 国产伦人伦偷精品视频| 又黄又粗又硬又大视频| 69av精品久久久久久| av福利片在线观看| 国产精华一区二区三区| 我的老师免费观看完整版| 狂野欧美白嫩少妇大欣赏| 亚洲av成人不卡在线观看播放网| 非洲黑人性xxxx精品又粗又长| 午夜免费成人在线视频| 欧美中文日本在线观看视频| 搡老岳熟女国产| 久久伊人香网站| 男女床上黄色一级片免费看| 日韩大尺度精品在线看网址| 亚洲av成人av| 免费在线观看日本一区| a级毛片a级免费在线| 特级一级黄色大片| 91久久精品国产一区二区成人 | 亚洲av成人精品一区久久| 日韩欧美三级三区| 一级毛片精品| 久久中文字幕人妻熟女| 成人高潮视频无遮挡免费网站| 免费人成视频x8x8入口观看| 一二三四在线观看免费中文在| 给我免费播放毛片高清在线观看| 国产爱豆传媒在线观看| 97超级碰碰碰精品色视频在线观看| 99re在线观看精品视频| 亚洲天堂国产精品一区在线| 久久精品影院6|