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

IBIS模型驗證與東北東部森林NPP季節變化模擬研究

2010-06-21 06:12:18國慶喜劉經偉
森林工程 2010年4期
關鍵詞:模型

劉 曦,國慶喜,劉經偉

(東北林業大學,哈爾濱 150040)

生態系統凈第一性生產力NPP(Net Primary Production)指綠色植物在單位時間和單位面積上所積累的有機干物質總量,是生物地球化學循環的關鍵環節,不僅是碳循環的原動力,而且是判定碳源、碳匯以及調節生態過程的主要因子[1]。其大小反映了植物固定和轉化光合產物的效率[2],同時也反映了環境因子對森林生長的綜合影響[3]。目前對大尺度的NPP動態變化的估計仍有很大的不確定性,因此更準確的估算NPP、理解它與氣象要素之間的相互關系,對于了解氣候變化對生態系統功能和結構的影響及判定碳匯具有重要作用,同時也是揭示生態系統對氣候變化響應的有效工具。

目前估算植被凈第一性生產力的模型很多,Ruimy[4]將模型概括為3類,即:統計模型 (statistical model)、參數模型 (parameter model)和過程模型 (process model)。IBIS模型屬于生態系統過程模型,該模型通過一種獨立、自然連續的模型框架結構將陸地表面生物物理、陸地碳通量和全球植被動態表達出來。模型集成了大范圍的生物物理、生理以及生態過程,并且這種模式框架能夠直接與大氣環流模式 (AGCMs)進行耦合,能夠模擬變化復雜的、時間跨度從秒到數百年的生物圈過程,屬于新一代全球生物圈模型DGVMs(dynamic global vegetation models),代表著全球碳循環模擬的研究方向。

本文首先基于實際觀測數據對IBIS模型有效性進行全面驗證,而后估算2004年和2005年NPP季節變化情況,旨在為氣候變化影響下大尺度NPP估算研究提供可靠的方法和數據支持。

1 研究方法

1.1 研究區域概況與樣點選取

研究地點位于黑龍江和吉林兩省境內,地跨40.9°~47.09°N 和125.0°~133.9°E,氣候為溫帶大陸性氣候,冬季漫長而寒冷,夏季短而多雨。年平均氣溫-4~6℃,1月和7月平均氣溫變化范圍分別是-22~-12℃,9~23℃。年平均降水變化范圍460~1200 mm,降水多集中于6~9月,夏季約占全年的60%以上;秋季次之,占年降水量的10%~25%;春季占10%~20%;冬季不足5%。

本研究黑龍江省境內樣點的選取,是基于2002年黑龍江省森林資源Ⅰ類清查數據,每隔8 km布點,共計542個樣點。吉林境內339樣點則是基于中國數字化植被圖集每隔15 km布點,各樣點所對應植被類型一目了然。這些樣地的位置如圖1所示。

圖1 研究用881樣點空間分布圖Fig.1 Study areas and spatial distribution of 881sample points

1.2 IBIS模型介紹

IBIS模型是美國威斯康星大學Foley教授等(1996年)研制的面向生物圈和區域尺度的景觀過程模型,發展至今已有最新版本IBIS 2.6[5]。IBIS模型本身包括地表、生物地球化學、植被物候和植被動態4個模塊,以氣候變量作為模型的驅動數據,不僅能對植被的長期動態進行模擬,還可以用于生態系統的碳、氮、水和能量平衡的模擬。在本模型中根據溫帶森林呼吸消耗碳的特點,NPP按照冠層凈光合與維持呼吸Rm(maintenance respiration)差值的70%計算。其中,光合計算采用Farquhar等,1980模型結構,受控于氣象條件 (溫度和降水)和植被生態生理學參數 (LAImax和Rubisco酶光合限制速率)。Rm即葉呼吸 (Rleaf)、根呼吸(Rroot)和干呼吸 (Rstem)的總稱,分別為參考溫度15℃時的Rubisco酶最大羧化速率 (Vm)、Jackson(1996,1997)根的擴散方程和LAImax的函數。具體計算方程如下:

式中:Ag為冠層凈光合速率,mol CO2/(m2·s1),η為生長呼吸消耗碳的比例,取值為0.3。

式中:Je表示為光限速率,與葉片吸收的光合有效輻射 (PAR)(W/m2)、植被CO2吸收的內在量子效率、葉片細胞間隙CO2濃度和光合作用的光補償點有關;Jc代表Rubisco酶光合限制速率,Js是當在具有高胞間二氧化碳濃度和強輻射條件下,才考慮的影響因子,Jc和Js都是關于Vm(Rubisco酶最大羧化速率)的方程。

式中:γ為葉片呼吸系數,Vm為葉片在15℃的Rubisco酶最大羧化速率,mol CO2/(m2·s1)。γ和Vm均因植被類型不同,取值不同。

式中:SLA代表單位重量的葉面積 (m2/kg),LAImax為每種植被類型的最大葉面積指數 (m2/m2),Troot為根部溫度 (K),根據Jackson(1996,1997)根的擴散方程計算。

式中:Tstem為上冠層干表面溫度 (K),Ps(sapwood)代表邊材所占比例,針對所有植被類型初始值均為0.05。

有關模型詳細描述見http://www.sage,wisc.edu/pages.datamodels.html或文獻[5-7]。

1.3 模型改進

本研究結合我國東北地區情況,將IBIS進行了相應改造。改造內容主要包括:

(1)將原有全球生物圈尺度的12種植物功能型劃分為適合我國東北地區的8種,即:溫帶常綠針葉、溫帶落葉闊葉、溫帶落葉針葉、寒溫帶常綠針葉、寒溫帶落葉針葉、寒溫帶落葉闊葉、落葉灌木和草本。

(2)原有的15種陸地植物群落類型修改成7種植被類型,即:溫帶常綠針葉林 (Temperate evergreen coniferous forest)、溫帶針闊葉混交林(Temperate coniferous and broadleaved mixed forest)、溫帶硬闊葉林 (Temperate hardwood forest)、溫帶軟闊葉林 (Temperate soft-leaved forest)、蒙古櫟林 (Mongolia quercus forest)、溫帶雜木林 (Temperate mixed deciduous forest)、寒溫帶興安落葉松林(Cold temperate dahurian larch forest)。

(3)土壤質地劃分同樣沿用美國土壤質地三角形,但土壤最大深度設置為1 m,因根據研究地實測數據及文章記載1 m以下土壤微生物活性幾乎為0。

(4)植被功能型參數,如Vmax,比葉面積,葉、根、干轉換時間常數,葉碳分配、根碳分配;植被生理學參數,如葉呼吸系數,氣孔導度m、b系數,絕對最小氣孔導度等均做了修正。

2 數據來源與數據分析

2.1 數據來源

在本研究中,IBIS模型的空間分辨率為0.5°*0.5°,模型運行需要輸入的數據為氣象數據、土壤特征數據和植被特征3方面的數據。

2 .1 .1 氣象數據

本研究所采用的氣象基礎數據來源于國家氣象中心氣象資料室發布的194氣象臺站信息,時間序列為1971年1月~2000年12月氣象數據,數據內容包括累年各月平均 (本站氣壓、氣溫、極端)、最高 (低)氣溫、相對濕度、總云量、風速、地溫、降水量、日照時數和蒸發量)以及氣象站點的經緯度和海拔高度。根據研究地樣點坐標,依據就近原則,選取對應氣象站信息。

2 .1 .2 土壤特征數據

包括土壤類型的顆粒組成數據,即沙粒含量、粘土含量、粉粒含量、土壤田間持水量和萎蔫含水量和容重等,來源于文獻[7]。

2 .1 .3 植被特征數據

包括植物功能型、每種植被類型、各樣點對應植被類型及枯落物儲量數據。植物功能型根據植被外貌 (喬木、灌木或草本)、葉行為 (常綠或落葉)、葉型 (針葉或闊葉)和光合途徑等4個指標將東北東部地區植被標劃分為8種植物功能型,每種植被類型葉面積指數 (LAI)及枯落物儲量數據來源于羅天祥博士論文 (http://www.geodata.cn或http://www.cern.ac.cn)。

2.2 數據分析

本文首先對IBIS進行了全面的驗證,這也是在區域模型研究方面的難點之一。本文采取和他人實測資料模擬結果相比較的方法,在一定程度上說明模型的有效性。所用模型驗證數據來源于2004年和2005年帽兒山生態站開展的NPP測定,共包括7種林型:楊樺林、硬闊葉林、蒙古櫟林、雜木林、紅松人工林和落葉松人工林。

模型估測結果準確性采用平均相對誤差(MRE)計算方法,計算公式如下:式中:n代表觀測點數量,Simulated(i)代表第i個觀測點模擬值,Observed(i)代表第i個觀測點的實測值。

3 結果與分析

3.1 模型驗證

3 .1 .1 氣溫與降水

本研究將2004年和2005年月平均溫度和降水的模擬值與氣象站點的平均觀測數據進行比較(如圖2和圖3所示),從圖中可以看出模型能夠很好的模擬溫度月變化趨勢,呈對稱鐘型,溫度1~7月逐漸上升,并在7~8月間達到最大值而后逐漸降低。模擬值與實測值之間R為0.997(2004),0.996(2005)(P<0.01)。模型模擬2004年和2005年1~4月氣溫均低于實測數據,此時段平均MRE分別是20.9%和47%,其中1~2月誤差較大。模型對生長季節 (5~10月)溫度模擬較好,MRE僅有1%和1.2%。其中,溫度最大值分別出現在8月 (2004年)和7月 (2005年),最大值溫度為21.4℃和21.5℃,模擬最大值與實測值之間的MRE分別為5%和5.6%。最低值均出現在1月份,溫度分別為-19.1℃和-16.2℃,與實測值之間的MRE分別為8.3%和11.3%。2004年平均溫度為3.7℃ ,高于2005年6.7%,其中生長季節平均溫度為16.4℃,高于2005年9.9%。

圖2 IBIS模擬2004年氣溫月動態與實測值比較Fig.2 Comparison between monthly dynamic simulation using IBIS and observed temperature in 2004

圖3 IBIS模擬2005年氣溫月動態與實測值比較Fig.3 Comparison between monthly dynamic simulation using IBIS and observed temperature in 2005

如圖4和圖5所示為降水模擬值與實測值對比圖,從圖中可以看出模型能夠較好的模擬降水量月變化趨勢,R為 0.979(2004),0.960(2005)(P<0.01),曲線變化呈草帽型。降水量1~7月逐漸增加,并在7月達到最大值,而后逐漸降低。1~4月降水模擬值均高于實測值,此時段MRE也較大,均超過60%,其中以3~4月降雨量估算誤差最大。生長季節降水模擬較好,MRE僅有23%(2004)和3%(2005)。其中,降水量最大值均出現在7月,最大降雨量分別為149.1 mm和149.9 mm,模擬最大值與實測值之間的MRE分別為 -11.7%和 -35%。最低值出現在 2月份(2004)和1月 (2005),分別為5 mm和6.3 mm,與實測值之間的MRE分別為285%和42%。2004年總降雨量為616.3 mm,低于2005年2%,其中生長季節降雨量為545 mm,低于2005年1%。模型模擬2004年和2005年兩年間的全年降水均低于實測值11.7%和35%。

IBIS對溫度和降雨的模擬結果顯示,模型能夠較好的模擬生長季溫度和降水,但對非生長季尤其是初春和冬季溫度和降雨 (雪)模擬不足,應是今后模型需要改進和完善的地方。

圖4 IBIS模擬2004年降水量月動態與實測值比較Fig.4 Comparison between monthly dynamic simulation using IBIS and observed precipitation in 2004

圖5 IBIS模擬2005年降水量月動態與實測值比較Fig.5 Comparison between monthly dynamic simulation using IBIS and observed precipitation in 2005

3 .1 .2年NPP模擬驗證

見表1為IBIS模擬NPP與觀測值比較,研究發現IBIS能夠很好的模擬不同林型NPP,其模擬值與觀測值MRE介于3.8%~17.96%之間,R值高達0.955。不同林型NPP誤差分別為紅松3.8%、針闊葉混交林5.87%、硬闊葉林6.55%、落葉松林7.73%、蒙古櫟林9.35%、楊華林15.56%和雜木林17.96%。根據NPP模擬結果,發現NPP值主要與林分類型[8]和所處的生長階段有關,針闊葉混交林作為溫帶森林生態系統中主要的林分類型,具有較高的NPP,其值介于0.769~0.799 kg C/(m2·a),落葉松林次之,NPP平均值為 0.737 kg C/m2,這與文獻[8,9]研究結果相似。硬闊葉林具有較低水平的NPP,年均值僅有0.311 kg C/(m2·a),這可能與所處的地理環境干燥和林齡較低有關。

森林NPP的形成是一個十分復雜的過程,影響因子也十分復雜,既有森林類型本身的生理生態學特性,也包括大量的環境因子,區域尺度的NPP估算模型無法考慮所有影響因子和生態過程,都進行了一定的簡化,因此模擬結果與實測值之間存在誤差在所難免[10],本研究有關NPP的模擬誤差均在文獻[11-14]范圍之內,因此IBIS模型的模擬結果總體上反映了東北東部森林生態系統的NPP實際情況,表明了該模型在本區域內的適用性。

表1 IBIS模擬東北東部2004年和2005年7種植被類型NPP與觀測值比較Tab.1 Comparison between the simulated using IBIS and observed of seven forest types NPP in the eastern part of northeast in 2004 and 2005

3.2 不同森林植被類型NPP的季節變化

如圖6和圖7所示為東北東部主要森林生態系統2004年和2005年NPP月變化曲線。從圖中可以看出NPP月變化呈對稱的鐘形,NPP由4月開始逐漸增加,至6~7月間達到最大值,這主要是由于夏季,太陽輻射豐富,水熱條件優良,才出現了峰值。從圖形來看,針葉樹種6月達到最大,闊葉樹種7月達到最大,闊葉樹種的NPP最大值均高于針葉樹種。不同植被類型在2004年NPP最大值分別是常綠針葉林0.15 kg C/(m2·month),針闊葉混交林0.18 kg C/(m2·month),硬闊葉林0.16 kg C/(m2·month),軟闊葉林 0.17 kg C/(m2·month),蒙古櫟林 0.18 kg C/(m2·month),雜木林0.13 kg C/(m2·month),興安落葉松林0.13 kg C/(m2·month)。2005年NPP最大值分別是常綠針葉林0.15 kg C/(m2·month),針闊葉混交林0.19 kg C/(m2·month),硬闊葉林0.17 kg C/(m2·month),軟闊葉林 0.17 kg C/(m2·month),蒙古櫟林0.18 kg C/(m2·month),雜木林0.16 kg C/(m2·month),興安落葉松林0.16 kg C/(m2·month)。2005年NPP最大值比2004年略有升高,升高幅度為5.6%~23% ,以雜木林和興安落葉松林NPP升高最明顯。2004和2005年最大值分別占全年凈初級生產力的20.2%~26.3% ,19.7%~23.4%,均是針葉林所占比例最小,蒙古櫟林所占比例最大,夏季NPP總量是全年NPP最大的季節,占全年NPP的53%~70%。

圖6 IBIS模擬2004年不同植被類型NPP月動態變化Fig.6 Monthly dynamic variation Simulation of NPP for different vegetation types using IBIS in 2004

圖7 IBIS模擬2005年不同植被類型NPP月動態變化Fig.7 Monthly dynamic variation simulation of NPP for different vegetation types using IBIS in 2005

7月末開始NPP逐漸減低,NPP在11月份接近于0,最低值出現在冬季12月份,這是因為冬季是我國氣溫最低、太陽輻射量最少的季節,NPP也是最小的一個月,大部分地方植被停止生長,NPP都幾乎為0,12月至次年3月NPP變化趨勢平穩。2004年NPP最低值分別是常綠針葉林0.000 58 kg C/(m2·month),針闊葉混交林 0.013 kg C/(m2·month),硬闊葉林0.009 8 kg C/(m2·month),軟闊葉林0.015 kg C/(m2·month),蒙古櫟林0.006 6 kg C/(m2·month),雜 木 林 0.033 kg C/(m2·month),興安落葉松林0.005 8 kg C/(m2·month)。2005年NPP最低值分別是常綠針葉林0.001 4 kg C/(m2·month),針闊葉混交林 0.019 kg C/(m2·month),硬闊葉林0.003 kg C/(m2·month),軟闊葉林0.007 kg C/(m2·month),蒙古櫟林0.001 3 kg C/(m2·month),雜 木 林 0.001 4 kg C/(m2·month),興安落葉松林0.001 7 kg C/(m2·month)。2005年不同植被類型NPP最低值與2004年相比,除硬闊葉林、雜木林和蒙古櫟林有較大幅度降低(71%~79%),以雜木林降低最明顯。其余樹種均有升高,升高幅度為6% ~137%。以常綠針葉林升高最顯著。

不同的季節內,不同植被類型NPP隨氣候變化的反應程度是各異的,這與Braswell等,1997的結論是一致的[15]。見表2和表3說明的是不同植被類型NPP隨季節變化的幅度差別,2004年和2005年NPP均以春季變化最劇烈,這是因為春季升溫促進植被萌發,各項生理活動增強。其中尤以針闊葉混交林的NPP變化幅度最顯著,這與文獻孫睿等[16,17]的研究截然相反,接下來依次為常綠針葉林和落葉松林,闊葉林的變化幅度一致。

秋季NPP是繼春季之后又一變化明顯的季節,仍然是針闊葉混交林NPP變化最為明顯。比較兩年間的季節變化 (見表4),2005年春季和秋季NPP較2004年變化劇烈,NPP在春季和秋季的平均值也高于2004年,這主要是因為2004年冬季低溫次年溫度升溫明顯,2005年夏季低溫多雨緩解區域干旱情況以及秋季降水比2004平均高7%,這些因素疊加造成了2005年春秋季NPP較高。

表3 2005年東北東部典型森林生態系統季節NPP統計特征值Tab.3 Season statistical eigenvalue of NPP for typical forest ecosystems in the eastern part of northeast in 2005(kg C·m-2·month-1)

表4 2004-2005年東北東部典型森林生態系統季節溫度 (℃)與降水 (mm)統計特征值Tab.4 Season statistical eigenvalue of temperature(℃)and precipitation(mm)for typical forest ecosystems in the eastern part of northeast in 2004-2005

4 結論與討論

4.1 模型驗證

IBIS模型能夠較好的模擬生長季溫度和降水,但對非生長季尤其是初春和冬季溫度和降雨 (雪)模擬不足,應是今后模型需要改進和完善的地方。對不同植被類型NPP的模擬誤差均小于眾文獻,因此IBIS模型的模擬結果總體上反映了東北東部森林生態系統的NPP實際情況,表明了該模型在本區域內的適用性。

4.2 NPP的季節變化

NPP的季節變化中,夏季NPP是全年最高的季節,占全年NPP的53%~70%。這主要是由于夏季,太陽輻射豐富,水熱條件優良所致。NPP最大值均出現在夏季6月和7月間,闊葉樹種NPP最大值大于針葉樹種,2005年NPP最大值比2004年略有升高,升高幅度為5.6% ~23%,這主要是由2005年夏季低溫多雨所致。

NPP最低值出現在冬季12月份,2005年不同植被類型NPP最低值與2004年相比,除硬闊葉林、雜木林和蒙古櫟林有較大幅度降低 (71%~79%),其余樹種均有升高,升高幅度為6% ~137%。以常綠針葉林升高最顯著。

2004年和2005年NPP均以春季變化最劇烈,針葉林NPP變化幅度大于闊葉樹種,這說明了針葉樹種對春季積雪融化,溫度升高反應很敏感,而闊葉樹種此時還沒有展葉,因此光合作用無法進行,NPP變化并不顯著。秋季NPP是繼春季之后又一變化明顯的季節,仍然是針闊葉混交林NPP變化最為明顯。

4.3 研究展望

本文僅運用森林生態系統兩年的NPP數據對IBIS模型進行了檢驗,系統分析了生態系統NPP的季節變化,但要深入分析生態系統NPP的季節變化規律及其環境控機制還需要更長時間的定位觀測與模擬研究,在長期實測資料的支持下,分析NPP年際變化及控制因子,評價森林生態系統在區域和全球碳循環中的作用。

從本研究結果看,NPP模擬與實測值較為一致,但數值上還存在一定的偏差。氣孔導度隨季節的變化仍需要進一步深入研究,模型應基于輔助試驗,注重對典型森林生態系統的長期觀測,加強光合和凈初級生產力等生物要素的同步觀測,增進生物因子對各呼吸組分作用的影響過程與機制的理解,完善生理模型,增強模型在不同尺度上運行的精度。

[1]Field C B,Behrenfeld M J,Randerson J T,et al.Primary Production of the Biosphere:Integrating Terrestrial and Oceanic Components[J].Science 1998,281:237-240.

[2]蘇宏新.全球氣候變化條件下新疆天山云杉林生長的分析與模擬 [D].北京:中國科學院,2005.

[3]張冬有.黑龍江省森林植被凈初級生產力遙感估算研究[D].北京:北京林業大學,2009.

[4]Ruimy S,Saugier B,Dedieu G.Methodology for Estimation of Terrestrial Net Primary Production from Remotely Sensed Data[J].Journal of Geophysics Research,1999,99:5263-5383.

[5]Foley J A,Prentice C,Ramankutty N,et al.An Integrated Biosphere Model of Land Surface Process,Terrestrial Carbon Balance and Vegetation Dynamics [J].Global Biogeochemical Cycles,1996,10:603-628.

[6]Foley J A,Kucharik C J,Polzin D.Integrated Biosphere Simulator Model(Ibis),Version 2.5.Model Product[EB,OL].A-vailable via http://daac.ornl.gov.2008,May 15.

[7]Kucharik C J,Foley J A,Delire C,et al.Testing the Performance of a Dynamic Global Ecosystem Model:Water Balance,Carbon Balance,and Vegetation Structure [J].Global Change Biology,2000,14:796 -825.

[8]Jiang H,Apps M J,Zhang Y,et al.Modelling the Spatial Pattern of Net Primary Productivity in Chinese Forests[J].Ecological Modelling,1999,122:275-288.

[9]Matsushita B,Tamura M.Integrating Remotely Sensed Data with an Ecosystem Model to Estimate Net Primary Productivity in East A-sia[J].Remote Sensing of Environment,2002,81:58-66.

[10]李 慧.福建省森林生態系統NPP和NEP時空模擬研究[D].福州:福建師范大學,2008.

[11]Esprey L J,Sands P J,Smith C W.Understanding 3-Pg Using a Sensitivity Analysis [J].Forest Ecology and Management,2004,193:235-250.

[12]Laurent J M,Francois L,Hen A B.European Bioclimatic Affinity Groups:Data-Model Comparisons[J].Global and Planetary Change,2008,61:28-40.

[13]Saltelli A,Tarantola S,Chan K.A Quantitative Model Independent Method for Global Sensitivity Analysis of Model Output[J].Technometrics,1999,41:39 -56.

[14]Thornton P E,Law B E,Gholz H L.Modeling and Measuring the Effects of Disturbance History and Climate on Carbon and Water Budgets in Evergreen Needle-Leaf Forests[J].Agricultural and Forest Meteorological,2002,113:185 -222.

[15]Mohamed M A A,Babiker I S,Chen Z M,et al.The Role of Climate Variability in the Interannual Variation of Terrestrial Net Primary Production(NPP)[J].Science of the Total Environment,2004,332:123 -137.

[16]孫 睿,朱啟疆.中國陸地植被凈第一性生產力及季節變化研究 [J].地理學報,2000,55(1):36-45.

[17]沈 微,王立海,孟 春.小興安嶺天然針闊混交林擇伐后土壤呼吸動態變化[J].森林工程,2009,25(3):1-4.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲福利视频一区二区| 欧美日韩精品在线播放| 国产又色又刺激高潮免费看| 免费观看无遮挡www的小视频| 久久大香伊蕉在人线观看热2| 亚洲精品国产综合99| 日本91视频| 成人噜噜噜视频在线观看| 日韩小视频在线播放| 91视频区| yy6080理论大片一级久久| 亚洲91精品视频| 亚洲中文无码av永久伊人| 亚洲永久色| 91网红精品在线观看| 中日韩一区二区三区中文免费视频 | 中国黄色一级视频| 香蕉99国内自产自拍视频| 亚洲天堂区| 久久国产热| 久久久四虎成人永久免费网站| 国产精品.com| 亚洲国产欧美国产综合久久| 性欧美久久| 国产本道久久一区二区三区| 自拍偷拍一区| 亚洲精品欧美日本中文字幕| 玖玖精品视频在线观看| 国产亚洲欧美在线专区| 青青网在线国产| 好紧太爽了视频免费无码| 欧美人人干| 亚洲国产高清精品线久久| 久久国产精品电影| 精品一区二区久久久久网站| 国内精自视频品线一二区| 伊人久久综在合线亚洲91| 国产精品xxx| 亚洲中文精品久久久久久不卡| 5388国产亚洲欧美在线观看| 亚洲自偷自拍另类小说| 婷婷激情亚洲| 91精品人妻互换| aⅴ免费在线观看| 亚洲欧美不卡视频| 久久青草热| 久久公开视频| 亚洲91精品视频| 99re这里只有国产中文精品国产精品 | a毛片在线| 波多野结衣无码中文字幕在线观看一区二区 | 亚洲欧美综合在线观看| 国产成人亚洲精品蜜芽影院 | 国产黑人在线| 精品欧美一区二区三区在线| 毛片网站观看| 日韩精品一区二区三区大桥未久| 成年看免费观看视频拍拍| 国产麻豆另类AV| 久久亚洲黄色视频| 伊人欧美在线| 国产欧美高清| 精品国产成人av免费| 久久精品国产亚洲麻豆| 91美女视频在线| 精品国产香蕉伊思人在线| 亚洲综合九九| 香蕉国产精品视频| 真人免费一级毛片一区二区| 中文字幕资源站| 国产美女久久久久不卡| 人妻出轨无码中文一区二区| 午夜丁香婷婷| 91亚洲视频下载| 欧美成人A视频| 国产精品污视频| 2021国产精品自产拍在线| 国产一线在线| 国产成人综合亚洲欧美在| 欧美性色综合网| 97se亚洲综合在线| 国产精品一区不卡|