999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

大興安嶺地區天然興安落葉松疏開木冠幅預測模型1

2017-05-17 11:40:49張樹森賈煒瑋王玉霞高慧淋
森林工程 2017年3期
關鍵詞:模型

張樹森,賈煒瑋,王玉霞,高慧淋

(1.黑龍江省尚志國有林場管理局,黑龍江 尚志 150600;2.東北林業大學 林學院 哈爾濱 150040)

大興安嶺地區天然興安落葉松疏開木冠幅預測模型1

張樹森1,賈煒瑋2*,王玉霞2,高慧淋2

(1.黑龍江省尚志國有林場管理局,黑龍江 尚志 150600;2.東北林業大學 林學院 哈爾濱 150040)

基于2016年大興安嶺地區十個林業局的天然興安落葉松林樣地調查數據,實測樣地內疏開木的胸徑與東、西、南、北四個方向的冠幅。分析疏開木胸徑與冠幅之間的關系,采用線性模型和非線性模型作為備選模型,分別構建了天然落葉松的冠幅預測模型,經過擬合優度指標與檢驗結果進行最優模型的選取。將所有樹木按照胸徑的變化范圍分為3個等級(即DBH≤20 cm、20 cm40 cm),基于選取的最優冠幅預測模型,采用啞變量的方法同時構建了3個等級天然興安落葉松疏開木冠幅的通用模型,并進一步分析3個等級冠幅模型之間的差別。結果表明:線性模型和非線性模型均能夠很好的擬合天然興安落葉松疏開木的冠幅模型,由擬合優度和檢驗指標最終選取非線性模型作為基礎模型。利用啞變量構建的3個等級疏開木冠幅的通用模型為CW=(0.353 3G1+0.225 5G2+1.537 6G3)DBH(0.587 5G1+0.738 5G2+0.190 4G3)。其中,DBH≤20 cm、20cm40 cm的兩個系數則沒有通過t檢驗。

天然興安落葉松林;冠幅預測模型;啞變量;大興安嶺地區

0 引言

在森林生態系統中,冠層是森林與外界環境進行相互作用最為直接的部分,同時也是林木進行營養循環的重要場所。樹冠能很好的反應樹木的長期競爭水平,林木通過樹冠進行光合作用和呼吸作用,產生并積累能量,供林木進行呼吸作用。樹冠通過林冠截留不僅能為林木的生長提供大量的水分,而且還能減緩雨水對地面的沖刷,防止水土流失。人們通常用材積三要素(形數、斷面積、樹高)計算林木材積,從而評價林木的經濟價值。斷面積的大小與林木樹冠的大小有著密不可分的關系。人們常用枝繁葉茂來評價一株林木的生長狀態。一般來說,同一立地條件下,同種林木、相同年齡,枝葉越茂密,冠幅值越大,林木生長就越好,獲得的單木經濟價值就越高。樹木在生長的過程中,如果遇到自然災害或者人為干擾,最先發現變化的是樹冠的變化,接著才會發現樹木的冠長、葉面積、樹冠表面積等也會發生變化。相比較而言,冠幅測量方法簡單、操作方便,而葉面積、樹冠表面積等需要大量的測量與計算,費時費力,比較麻煩。因此,人們通常用樹冠的大小來評價樹木是否生長健康,以及林木質量的優劣。胸徑和冠幅都是用來描述樹冠生長的重要指標,冠幅不但可以用來計算樹木的競爭指數,而且可以預測樹高或胸徑生長量等。

樹冠研究指標一般分為兩種,一種是冠幅、葉面積指數等能夠直接測定的指標,另外一種是疏透度、圓滿度、樹冠率、樹冠競爭因子(CCF)等需要間接計算的指標。其中人們最常用的是冠幅指標,其次是樹冠競爭因子指標。在許多森林經營管理中,樹冠指標常常被用來作為制定營林決策的重要依據之一。大量研究表明,冠幅與胸徑有著明顯的相關關系。Curtis等是第一個使用線性回歸法,構建黃杉(Pseudotsugamenziesii(Mirb.)Franco)冠幅生長模型[1]。構建單木冠幅模型的方法相對比較多,但是線性回歸模型比較多[2-4]。國內外學者對單木冠幅模型進行了大量研究。例如,唐守正為了計算樹冠競爭因子(CCF),使用線性模型、冪函數和雙曲線函數模型來擬合杉木(Cunninghamialanceolata(Lamb.)Hook)和長白落葉松(LarixolgensisHerry)疏開木與胸徑關系模型,發現雙曲線函數模型效果略好一些[5]。符亞健等考慮單木冠長、每公頃株數和林木競爭壓力指數等因子,從9個常見的冠幅-直徑備選模型中選出一個擬合精度較高的三參數的邏輯斯蒂形式的冠幅-直徑模型作為構建華北落葉松(Larixprincipis-rupprechtiiMayr)單木冠幅模型的基礎模型,從而更好地模擬當地華北落葉松冠幅與胸徑的關系[6]。自21世紀初以來,對樹木冠幅的研究得到廣泛的重視,迄今為止,針對冠幅模型的構建已提出了多種建模方法,例如,Gill等利用加權最小二乘法擬合模型,構建加利福尼亞幾種針葉樹種樹冠半徑模型[7]。雷向東等利用SPSS軟件,采用逐步回歸分析法選取所需自變量,采用多元回歸模型建立了長白落葉松、冷杉(Abiesnephrolepis)、紅松(Pinuskoraiensis)、云杉(Piceakoraiensis)等9個樹種的冠幅預測模型[8]。符利勇等以湖南省黃豐橋國有林場調查的杉木為研究對象,把立地指數和嵌套在立地指數中的樣地作為隨機效應因子,利用2水平非線性混合效應模型構建混合效應模型單木冠幅預測模型[9]。

興安落葉松(Larixgmelinii(Rupr.)Kuzen),為松科、落葉松屬的落葉喬木,是中國東北林區的主要森林樹種,是大興安嶺優良的鄉土樹種,同時也是大興安嶺原始森林的主要組成樹種,主要分布于大、小興安嶺,蘇聯遠東地區也有分布。興安落葉松喜光性強,在各種不同環境均能生長,常組成大面積的單純林,興安落葉松木材蓄積豐富,也是東北地區進行荒山造林和更新的主要樹種。興安落葉松木材略重,硬度中等,易裂,邊材淡黃色,心材黃褐色至紅褐色,紋理直,結構細密,相對密度為0.32~0.52,有樹脂,耐久用等特點,可供房屋建筑、土木工程、電桿、舟車、細木加工及木纖維工業原料等用材,樹干可提取樹脂,樹皮可提取栲膠,是東北地區的主要速生豐產林木,也是我國重要的用材樹種。本次研究根據大興安嶺地區天然落葉松疏開木的胸徑與冠幅之間的關系,通過擬合與檢驗,構建以胸徑作為單一變量來預測冠幅的大小的模型。將樹木按照胸徑大小分為3個等級,即DBH≤20 cm,20 cm40 cm,采用啞變量的方法構建了天然落葉松疏開木3個等級冠幅的通用模型,并對3個等級的冠幅進行了比較,以期為今后大興安嶺地區森林可持續經營提供科學依據。

1 研究地區概況數據收集

1.1 研究地概況

大興安嶺地區地處我國最北部邊陲,東與小興安嶺毗鄰,西以大興安嶺山脈為界與內蒙古自治區接壤,南瀕廣闊的松嫩平原,北以黑龍江主航道中心線與俄羅斯為鄰。大興安嶺地區天然地處黑龍江省西北部,內蒙古自治區東北部,大興安嶺山脈的東北坡,處于根河與嫩江兩條地震大斷裂帶上,是黑龍江省下轄地區,屬黑龍江省十三個地級行政單位之一。大興安嶺地區位于東經121°12′到127°00′,東西橫跨6個經度,北緯50°10′到53°33′,南北縱越3個緯度,屬寒溫帶氣候區,年降水量偏低(大約為350~500 mm),年平均氣溫-2~-4 ℃,有的地方極端最低溫-52.3 ℃。大興安嶺地區是我國重點國有林區和天然林主要分布區之一,也是中國唯一的寒溫帶針葉林區和國內僅存的寒溫帶生物基因庫,森林覆蓋率79.83%。大興安嶺地區疆域廣闊,東西長410 km,南北寬386 km,行政區面積8.46萬km2,包括塔河、新林、西林、呼中、加格達奇、松嶺、韓家園、十八站、阿木爾和圖強十個林業局。

1.2 數據來源

本次研究的數據是來自2016年大興安嶺地區的十個林業局,研究區域內有天然分布的興安落葉松林。在每一個林業局內選擇有地域代表性,體現不同立地條件差異的不同徑階的興安落葉松疏開木為調查對象,利用胸徑尺測量距離地面1.3處的直徑,稱為胸徑(DBH),胸徑范圍為3.6~50.2 cm,精確到0.1 cm。利用鋼尺以樣木為中心,分別對其正東、正西、正南、正北4個方向的最大樹冠半徑進行測量,精確到0.1 m。統計不同徑階樹木的基本信息(表1)。統計東、西、南、北 4 個方向的冠幅半徑,求出單木冠幅的平均值,其計算公式為

式中:CW東、CW西、CW南、CW北分別代表以樣木為中心所測量的東、南、西、北4 個方向的最大樹冠半徑,CW是指單木冠幅,m。

將所有數據按照3∶1的比例隨機的分為擬合樣本和檢驗樣本,擬合樣本用于模型參數的估計,檢驗樣本用于對模型的預測能力進行評價。

表1 不同徑階樹木變量基本信息統計Tab.1 Statistical information for the individual tree variables from different DBH class

2 研究方法

2.1 備選模型的選取

Krajicek等認為冠幅(CW)與胸高直徑(DBH)呈線性相關關系[10]。后來研究者發現胸徑生長與冠幅大小的相關規律不受立地條件與林齡差異的影響[11]。國內外對冠幅和胸徑關系進行了大量的研究,其中有9個常用的冠幅-胸徑模型[12-15],本文選取了幾個常用的模型進行了研究和比較。

基于以往研究成果,通過對天然落葉松疏開木胸徑(DBH)與冠幅(CW)的關系進行分析,以樣木胸徑(DBH)為自變量,實測冠幅(CW)為因變量,做DBH-CW的散點圖。根據散點圖的分布趨勢對散點圖進行分析,可看出隨著DBH的越大,CW也呈現增大的趨勢。本文主要構建以下3個備選模型,然后用SAS9.3軟件PROC Model模塊進行模型的擬合。備選模型的形式如下所示:

CW=a1+a2DBH。

(1)

CW=a1+a2DBH+a3DBH2。

(2)

CW=a2DBHa2。

(3)

式中:CW為單木冠幅,DBH為樣木胸徑,a1、a2、a3為模型參數。

2.2 基于啞變量的不同等級樹木通用冠幅模型的建立

近年來,啞變量模型的研究越來越受到人們的重視。在使用最優模型的基礎上,在不增加自變量的前提下,可以使用啞變量方法構建模型,從而提高模型的預測精度,這在以前就有所研究。例如,符利勇等對東北地區興安落葉松和長白落葉松兩個主要樹種構建生物量方程時,考慮林分不同起源能提高方程預測,因此把林分起源差異當作固定效應,通過啞變量方法構建不同樹種生物量模型方法,增加自變量能提高方程預測效果[16]。

常用啞變量模型方法進行各種回歸分析和建模實踐[17-19]。啞變量為虛擬變量,啞變量是對定性因子或分類變量進行處理的一種常用方法,一般取值為0 或1。啞變量的定義為:對于等級性(定性)數據x,用變量δ(x,i)表示,當x取值為第i等級時,δ(x,i)=1,否則為0[19]。

2.3 模型擬合及檢驗指標

利用均方根誤差(RMSE)和調整后的決定系數(R2)兩個指標對模型的擬合優度進行評價。當R2的值越大,而RMSE的值越小時,說明模型的擬合效果越高,模型愈適用。再分別對模型進行獨立性檢驗,計算指標平均誤差(ME)、平均絕對誤差(MAE)、平均誤差絕對值(MA%E)的值。ME、MAE、MA%E的值越小,說明模型的預測效果越好。各指標計算公式如下:

(4)

(5)

(6)

(7)

(8)

3 結果

3.1 最優模型的選取

表2 各個備選模型的擬合優度及檢驗結果比較Tab.2 Goodness-of-fit and validation results for the basic models

3.2 不同徑階的冠幅通用模型的建立

3.2.1 構建啞變量模型

通過對模型進行擬合與檢驗,得出模型(3)為最優基礎模型。本次研究收集的胸徑數據,考慮胸徑值之間的差異性,將胸徑分為0~20 cm、20~40 cm、>40 cm3個徑級,將這3個不同的徑級帶入基礎模型中,引入啞變量,將不同徑級范圍的樣木用定性代碼表示,從而整合成一個模型來構建。當G1=1、G2=0、G3=0時,代表徑階為0~20c m的冠幅預測模型;當G1=0、G2=1、G3=0時,代表徑階為20~40 cm的冠幅預測模型;當G1=0、G2=0、G3=1時,代表>40 cm徑階??紤]到落葉松的胸徑大小不同,預測的冠幅數值可能會存在一定的差異。在所構建的冠幅預測模型的基礎上,通過引入啞變量的方法,構建了不同徑階的大興安嶺地區天然落葉松疏開木冠幅的通用模型,模型的具體形式如下所示。

CW=(a1G1+a2G2+a3G4)DBH(b1G1+b2G2+b3G3)。

(9)

表3 基于啞變量人工落葉松不同徑階大小林木冠幅模型的參數估計結果Tab.3 Estimates for the parameters of the crown width model for the different classes of natural Larix gmelini plantation based on dummy variable

3.2.2 消除異方差,繪制殘差圖

做各徑階冠幅預測模型的殘差圖,通過觀察模型殘差分布圖,發現殘差值的分布均勻,殘差值分布的離散程度有擴大的趨勢,表明預測模型存在一定程度的異方差性,因此需要消除異方差,使殘差值均勻分布。消除異方差,常用方法有對數回歸法和加權回歸法[20-23]。本文通過加權的方法消除異方差,得到分布相對均勻的殘差分布圖(圖1)。

(a)

(b)

(c)圖1 各徑階冠幅預測模型消除異方差后的殘差圖 (a、b、c分別代表≤20 cm、20 cm<40 cm、>40 cm 3個徑階的殘差圖)Fig.1 Residual plots for the predicted crown width model of different classes after the eliminating for the heteroscedasticity(a,b,c represents the residual plots for ≤20 cm、20 cm<40 cm、>40 cm,respectively)

4 結論與討論

本研究是收集了2016年大興安嶺地區十個林業局分布的天然興安落葉松的胸徑與冠幅數據。通過分析胸徑與冠幅之間的關系,選取了3個備選模型,經過擬合與檢驗,得出冪函數為最優基礎模型。然后將最優基礎模型與啞變量結合起來,將胸徑分為0~20 cm、20~40 cm、>40 cm3個徑級,利用啞變量構建了大興安嶺地區天然落葉松冠幅的通用預測模型。采用SAS9.3軟件PROC Model模塊進行參數的求解,得出啞變量最優模型:CW=(0.353 3G1+0.225 5G2+1.537 6G3)DBH0.587 5G1+0.738 5G2+0.190 4G3。利用加權法消除冠幅預測值的殘差,得到比較均勻的殘差圖。由啞變量參數的估計結果可知,不同徑階的冠幅模型存在顯著的不同,因此建議在今后的研究中將林分按照徑階進行分組建模,可以在一定程度上提高模型的擬合效果和預測精度。

雖然胸徑生長與冠幅大小呈正相關關系的規律不受立地條件、林分密度、林齡等的影響,但是具體的模型及其參數還是會有或多或少的不同。大興安嶺地區地域遼闊,溫度、降雨量等差別比較大,而且各個林場的立地條件也大不相同,每個林場、每個林班的密度也不盡相同,甚至每棵樹的林齡以及所在的坡度、坡位、坡向也有一定的不同,而且樹高、枝下高、冠長比等因子與冠幅也存在顯著相關性。若要考慮以上相關因素,就需要收集與之相關的數據,然而本次研究只有林場及冠幅與胸徑的數據,因此如何綜合考慮冠幅生長相關的各種因子,構造更準確的冠幅模型,仍需要進一步的研究。

[1]Curtis R O,Reukema D L.Crown Development and Site Estimates in a Douglas-Fir Plantation Spacing Test[J].Forest Science,1970,16(3):287-301.

[2]Grote R.Estimation of Crown Radii and Crown Projection Area from Stem Size and Tree Position[J].Annals of Forest Science,2003,60(5):393-402.

[3]田曉筠.林木直徑與冠幅的灰色模型及應用[J].科技資訊,2008(25):86.

[4]韋雪花,王佳,馮仲科.北京市13個常見樹種胸徑估測研究[J].北京林業大學學報,2013,35(5):56-63.

[5]唐守正,杜紀山.利用樹冠競爭因子確定同齡間伐林分的斷面積生長過程[J].林業科學,1999,35(6):35-41.

[6]符亞健,呂飛舟,朱光玉,等.華北落葉松天然次生林單木冠幅模型構建[J].林業資源管理.2016(5):65-70.

[7]Gill S J,Biging G S,Murphy A C.Modeling Conifer Tree Crown Radius and Estimating Canopy Cover[J].Forest Ecology and Management,2000,126(3):405-416.

[8]雷相東,張則路,陳曉光.長白落葉松等幾個樹種冠幅預測模型的研究[J].北京林業大學學報,2006,28(6):75-79.

[9]符利勇,孫華.基于混合效應模型的杉木單木冠幅預測模型[J].林業科學,2013,49(8):65-74.

[10]Krajicek J E,Brinkman K,Gingrich S F.Crown Competition Factor-a Measure of Density[J].Forest Sci,1961,7(1):35-42.

[11]張金友,姜毛孩,劉德祥.小興安嶺人工樟子松林立木胸徑與冠幅直徑的相關性[J].林業勘察設計,1997(4):51.

[12]羅玲,廖超英.榆林沙區樟子松冠幅與胸徑的相關關系分析[J].安徽農學通報,2007,13(24):92-97.

[13]盧昌泰,李吉躍,康強,等.馬尾松胸徑與根徑和冠徑的關系研究[J].北京林業大學學報,2008,30(1):58-63.

[15] S?nmez T.Diameter at Breast Height-crown Diameter Prediction Models for Picea orientalis[J].African Journal of Agricultural Research,2009,4(3):214-219.

[16]符利勇,唐守正,張會儒,等.東北地區兩個主要樹種地上生物量通用方程構建[J].生態學報,2015,35(1):150-157.

[17]李希菲,洪玲霞.用啞變量法求算立地指數曲線族的研究[J].林業科學研究,1997,10(2):215-219.

[18]李河,麥勁壯,肖敏,等.啞變量在 Logistic 回歸模型中的應用[J].循證醫學,2008,8(1):42-45.

[19]高東啟,鄧華鋒,王海賓,等.基于啞變量的蒙古櫟林分生長模型[J].東北林業大學學報,2014,42(1):61-64

[20]曾偉生,唐守正,夏忠勝,等.利用線性混合模型和啞變量模型方法建立貴州省通用性生物量方程[J].林業科學研究,2011,24(3):285-291.

[21]曾偉生,唐守正.東北落葉松和南方馬尾松地下生物量模型研建[J].北京林業大學學報,2011,33(2):1-6.

[22]楊英,冉啟香,陳新云,等.啞變量在云杉地上生物量模型中的應用研究[J].林業資源管理,2015(6):71-76.

[23]張甜,朱玉杰,董希斌.撫育間伐對大興安嶺天然用材林冠層結構及光環境特征的影響[J].東北林業大學學報,2016,44(10):1-7.

Modeling Crown Width for the Opening-grown Trees of NaturalLarixGmeliniPlantation in Daxing’an Mountains1

Zhang Shusen1,Jia Weiwei2*,Wang Yuxia2,Gao Huilin2

(1.Heilongjiang Shangzhi National Forest Farm Management Bureau,Shangzhi 150600;2.School of Forestry,Northeast Forestry University,Harbin 150040)

Based on the investigation data for the naturalLarixgmeliniplantationfrom ten forest bureaus in Daxing’an Mountains,the diameter at the breast height(DBH)and the crown width from four directions(i.e.,east,west,south and north)for all the opening-grown trees in the sample plots were measured.The relationship between theDBHand crown width was analyzed.The linear equation and nonlinear equation were used as the basic models to develop the predicted crown width model for the naturalLarixgmeliniplantation and the optimal model was selected based on the goodness-of-fit and validation results.All the sample trees were classified into three categories(DBH≤20cm,20cm40cm)based on the diameter class.Using the optimal model selected by the goodness-of-fit and validation statistics,the dummy variable approach was employed to develop the general common crown width model for the three classes.The differences of crown width between the three categories were further analyzed.The results showed that both of linear and nonlinear equations were suitable to model the crown width of naturalLarixgmeliniplantation in Daxing’an Mountains.Considering the goodness-of-fit and validation statistics,the nonlinear equation was finally selected as the best model.The general common crown width model for the three categories of the naturalLarixgmeliniplantation wasCW=(0.353 3G1+0.225 5G2+1.537 6G3)DBH(0.587 5G1+0.738 5G2+0.190 4G3).TheparametersofDBH≤20cm,20cm40cmwerenotpassedthettest.

NaturalLarix gmeliniplantation;crownwidthpredictedmodel;dummyvariable;Daxing’anmountains

2017-1-20

橫向課題(黑龍江省市縣林區森林植被空間分布及多樣性研究);黑龍江省林業廳項目(黑龍江省森林可持續經營試驗示范區建設);大興安嶺地區科技計劃項目(大興安嶺林區林分密度調控技術研究及示范)

張樹森,學士,高級工程師。研究方向:森林經營。E-mali:zssymp@126.com

張樹森,賈煒瑋,王玉霞,等.大興安嶺地區天然興安落葉松疏開木冠幅預測模型[J].森林工程,2017,33(3):33-38.

U270.1

A

1001-005X(2017)03-0033-06

*通信作者:賈煒瑋,博士,副教授。研究方向:林分生長模型。E-mali:Jiaww2002@163.com

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 欧美一级特黄aaaaaa在线看片| 中文字幕资源站| 久草视频中文| 久久国产黑丝袜视频| 日韩精品久久无码中文字幕色欲| 少妇高潮惨叫久久久久久| 中文字幕无码电影| 亚洲免费毛片| 99国产在线视频| 亚洲三级片在线看| 日本手机在线视频| 四虎国产精品永久一区| 无码精品国产dvd在线观看9久 | 亚洲国产日韩在线观看| 在线观看网站国产| 青青青视频91在线 | 91九色国产porny| 亚洲欧美日韩中文字幕在线| 国产成人免费高清AⅤ| 亚洲一区二区精品无码久久久| 日韩毛片免费观看| 成人福利在线视频| 久久国产拍爱| 日韩精品无码免费一区二区三区| 精品伊人久久大香线蕉网站| 欧美曰批视频免费播放免费| 成人亚洲天堂| 国内嫩模私拍精品视频| 无码AV动漫| 91在线播放国产| 免费国产无遮挡又黄又爽| 性欧美精品xxxx| 国模视频一区二区| 免费在线成人网| 一区二区三区在线不卡免费| 国产网站免费观看| 国产成人综合在线视频| 亚洲欧洲日产国码无码av喷潮| 国产成人a毛片在线| 91成人免费观看| 成人综合网址| 日韩欧美网址| 国产自产视频一区二区三区| 亚洲色图另类| 亚洲人成网18禁| 欧美日韩一区二区在线播放| 国产日产欧美精品| 狠狠久久综合伊人不卡| 爆操波多野结衣| 久久精品这里只有国产中文精品 | 在线欧美国产| 欧美成a人片在线观看| 偷拍久久网| 欧美精品一区在线看| 国产永久在线视频| a毛片基地免费大全| 国产精品免费福利久久播放| 国产精品专区第1页| 免费一级无码在线网站| 91无码国产视频| 婷婷综合缴情亚洲五月伊| 亚洲一级无毛片无码在线免费视频 | 亚洲成人播放| 国产精品页| 亚洲中文精品久久久久久不卡| 国内嫩模私拍精品视频| 亚洲香蕉在线| 狠狠亚洲婷婷综合色香| 99免费在线观看视频| 香蕉蕉亚亚洲aav综合| 国产在线观看高清不卡| 亚洲欧美自拍一区| 99尹人香蕉国产免费天天拍| 日韩成人在线网站| 国产精品午夜电影| 国产一二三区在线| 日本精品一在线观看视频| 91热爆在线| 99热最新网址| 欧美一级99在线观看国产| 国产精品成人一区二区不卡 | 在线另类稀缺国产呦|