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

基于多源遙感數據的黃河三角洲人工刺槐林生物量估算

2023-12-29 00:00:00汪逸聰
湖北農業科學 2023年7期

摘要:利用哨兵影像、數字地形數據及森林實地樣方調查數據,分別構建K-近鄰(KNN)模型、隨機森林(RF)模型、極值梯度增強(XGBboost)模型、Stacking模型,實現對黃河三角洲人工刺槐(Robinia pseudoacacia)林生物量的估算。結果表明,相較于K-近鄰模型、隨機森林模型、極值梯度增強模型,集成學習Stacking模型明顯提高了生物量估測的精度(R2=0.61、RMSE=13.42 t/hm2)。

關鍵詞:哨兵;Stacking模型;刺槐(Robinia pseudoacacia)林;生物量;黃河三角洲

中圖分類號:S75" " " " "文獻標識碼:A

文章編號:0439-8114(2023)07-0143-06

DOI:10.14088/j.cnki.issn0439-8114.2023.07.025 開放科學(資源服務)標識碼(OSID):

Abstract: Using sentinel images, digital terrain data and forest field quadrat survey data, K-nearest neighbor (KNN) model, random forest (RF) model, extreme gradient enhancement (XGBboost) model and Stacking model were constructed respectively to estimate the biomass of artificial Robbin pseudoacacia forest in Yellow River Delta. The results showed that the integrated learning Stacking model significantly improved the accuracy of biomass estimation compared with K-nearest neighbor model, random forest model, and extreme gradient enhancement model (R2=0.61, RMSE=13.42 t/hm2).

Key words: sentinel; Stacking model; Robbin pseudoacacia forest; biomass; Yellow River Delta

為了應對氣候危機,全球178個締約方共同簽署了《巴黎協定》。該協定旨在將全球平均氣溫較前工業化時期上升幅度控制在2 ℃以內[1]。中國作為世界第一工業大國,碳排放量居世界首位,為積極應對氣候問題,彰顯負責任的大國形象與擔當,中國于2020年9月明確提出“雙碳”目標,即2030年實現“碳達峰”和2060年實現“碳中和”[2]?!疤歼_峰”是指二氧化碳排放量達峰值后逐步下降,“碳中和”是指采用減排和增匯措施抵消自身產生的二氧化碳排放,從而實現二氧化碳“零排放”的目標[3]。

森林生態系統作為陸地生態系統中最主要的組成部分,通過植被光合作用,吸收固定二氧化碳,是陸地生態系統中最大、最重要的貯碳庫,對于“雙碳”目標的實現起到至關重要的作用[3]。有學者研究表明,中國森林碳儲量增匯的主要貢獻來自人工林[4]。植樹造林被認為是最具生態效應的固碳減排方法,中國目前已廣泛開展了森林保護、植樹造林等生態工程,人工林面積已達7 954.28萬hm2,占世界人工林總面積的73%,居世界首位[5]。

黃河三角洲地區由于土地生產力低下、土壤鹽堿化等原因,該地區的原生喬木難以生存。刺槐(Robbin pseudoacacia)因具有耐旱和抗鹽堿的特性而被引入并廣泛種植,成為該地區最重要的生態保護林。然而,自20世紀90年代初,大面積刺槐出現枯梢、死亡[6],造成當地水土流失嚴重、土壤貧瘠以及植被破壞。為了研究黃河三角洲人工林碳循環和碳匯能力,有必要對刺槐生物量進行精確估算和動態變化分析。

森林生物量是評價森林碳匯能力的關鍵指標,同時也反映森林的健康狀況和生產力水平。精確估算森林生物量是人工林科學管理和可持續經營的前提[7]。森林生物量的測定方法主要包括野外實地調查和基于遙感方法進行估算。對于小型林分,采用野外實地調查的結果更精確可靠。但這種方法過于耗時費力且對生態環境具有破壞性,難以推廣到大區域尺度森林生物量的測定[8]。而基于遙感方法可以實現大面積對地觀測,且獲取信息的效率高、時間周期短。遙感可以有效監測大區域尺度森林生物量的連續空間分布及動態變化[9]。雷達是一種主動遙感技術,具備穿透樹冠并直接與樹木生物量的主要部位(如樹干和樹枝)相互作用的能力[10]。雷達后向散射強度在一定范圍內隨著森林生物量的增加而增加,且波長越長對森林生物量的敏感度越高[11]。因此,長波段更適用于估算生物量[12]。然而,由于大部分長波雷達衛星是商業衛星,數據獲取和處理的費用相對較高,限制了長波雷達在生物量估算中的應用和推廣。歐洲航天局C波段哨兵1號任務是免費為全球用戶提供高空間分辨率的合成孔徑雷達數據。由于C波段的飽和點較低,所以僅靠雷達數據不足以估算森林生物量。除雷達數據外,光學遙感數據可以提取到大量與森林生物量密切相關的信息,如與植被光合作用直接相關的藍紫光以及紅光、植被指數等[13]。光學遙感目前已廣泛應用于區域尺度的森林生物量估算,但光學遙感穿透力弱,僅能與森林上部的冠層結構發生作用且信號易飽和,導致估測精度偏低[14]。各種遙感源在森林生物量估算中各具優勢,但也均存在局限性,因此,整合各數據源的優點并進行生物量估算是當今森林生物量估算的熱點。為了更好地對黃河三角洲人工刺槐林進行碳匯評估和科學可持續管理,本研究通過獲取野外實測樣方數據、哨兵1號C波段雷達數據、哨兵2號多光譜數據、數字地形數據和2017年課題組基于激光雷達反演得到孤島區域高精度生物量分布圖,利用集成學習方法構建森林生物量估算模型,從而獲得黃河三角洲人工刺槐林森林生物量空間分布情況并展開分析,為林業部門科學管理人工林提供依據。

1 研究區概況

黃河三角洲隸屬于山東省東營市,位于山東省東北部渤海南岸與萊州灣西岸,地處117°30′—119°20′E和36°55′—38°16′N,總面積約為5 400 km2。黃河三角洲海拔為0~15 m,地勢平緩,屬于暖溫帶大陸性季風氣候,夏季雨水充沛,年均降水量為530~630 mm,夏熱冬涼,年均氣溫為12 ℃。刺槐原產自北美洲,根系發達且適應力較強,在固沙保土方面表現突出,因此,在20世紀70年底刺槐被引入黃河三角洲區域并廣泛種植,形成了中國最大的人工刺槐林地(圖1)。

2 數據處理與方法

2.1 數據

1)野外實測樣方數據。2017年6月中旬在孤島進行野外實測調查,共設置72個10 m×10 m的樣方。樣方的生物量調查采用標準木法,即選擇1株樹作為標準木,使用激光測高儀測量其樹高,用皮尺測量其樹干直徑,根據刺槐林的異速生長方程計算標準木的生物量,并乘以樣方內的樹木總株數,計算每個樣方的生物量。

2)哨兵2號。哨兵2號數據來自哥白尼開放獲取中心(https://scihub.copernicus.eu/),分辨率為10 m×10 m,影像獲取時間為2017年6月21日,含云量低于2%。

3)哨兵1號。哨兵1號GRD數據來自哥白尼開放獲取中心(https://scihub.copernicus.eu/),獲取時間為2017年6月8日,分辨率為5 m×20 m,雙極化模式為VV、VH。

4)數字地形數據。數字高程模型是通過向山東省國土測繪院申請獲得,獲取時間為2013年,其空間分辨率為5 m,坐標系統為2000國家大地坐標系,為統一數據尺度將其重采樣至10 m分辨率。

5)高精度孤島生物量預測值。Lu等[15]基于無人機激光雷達和背包移動激光雷達,結合野外實測數據,對孤島不同健康等級的刺槐林進行生物量估測,結果表明,基于隨機森林方法3種不同健康等級的刺槐林估算精度分別為:健康林的決定系數(Coefficient of determination,R2)為0.92,均方根誤差(Root mean squared error,RMSE)為3.99 t/hm2,中度枯梢林的決定系數(R2)為0.95,均方根誤差(RMSE)為4.65 t/hm2,重度枯梢林的決定系數(R2)為0.91,均方根誤差(RMSE)為4.63 t/hm2。

基于激光雷達預測的高精度孤島生物量預測值可代替野外實測數據,本研究將生物量預測值作為訓練樣本,從而降低數據采集的成本和節約時間。此外,選用激光雷達提供的預測值作為訓練樣本,可以避免野外實測數據存在的偏差和誤差等問題,提高模型的精度和可靠性。

2.2 數據預處理

歐洲航天局發布的Sen2cor是針對哨兵2號L2C級數據預處理生成L2A級數據的處理器,它基于大氣輻射傳輸模型的方法對C級產品進行大氣、卷云以及地形校正,從而生成A級數據。由于哨兵2號波段分辨率分別為10、20、60 m,為提高生物量反演的精度,本研究將其他波段重采樣至10 m分辨率。

SNAP軟件是由歐洲航天局開發并用于哨兵數據的處理軟件,哨兵1號預處理的主要流程包括熱噪聲去除、軌道文件校正、輻射定標、濾波校正、多普勒地形校正,最終獲得圖像的后向散射強度信息,最后將圖像重采樣至10 m分辨率。

2.3 研究方法

2.3.1 集成學習 Stacking是一種集成學習策略,通過元模型將多個單一模型組合起來,這些單一模型被稱為基模型。Stacking集成學習可以有效提高模型的魯棒性和泛化能力,但基模型的表現會直接影響Stacking模型的最終效果。因此,在選擇基模型時,需要考慮學習器的充分性和多樣性,即基模型具有良好的學習能力,且各基模型之間相互獨立,以此實現模型間信息的有效互補[16]。

K-近鄰(K-nearest neighbors,KNN)模型、極值梯度增強(eXtreme gradient boosting,XGBboost)模型和隨機森林(Random forest,RF)模型在回歸問題上表現優異,具有較強的非線性擬合能力。KNN模型的理論基礎成熟且易于理解,XGBboost模型可根據弱學習器性能調整樣本分布,并賦予錯誤結果更大權重,通過不斷學習調整的樣本,對所有弱學習器進行加權組合以獲得最終結果[17]。RF模型是在Bagging模型基礎上進行改進,它通過自助聚集的方式對樣本進行隨機有放回抽樣,并在每個節點上選擇最佳特征進行劃分,從而實現高效訓練[18]。這3種回歸模型滿足基模型選擇的充分性和多樣性原則,因此被廣泛用于Stacking模型中。本研究將上述3種模型與線性回歸(Linear regression,LR)相結合,構建Stacking模型以估算黃河三角洲人工刺槐林的生物量。

2.3.2 特征參數提取 光學遙感技術可以通過不同波段的反射率來獲取地表特征的信息。由于Band1、Band9和Band10波段主要用于大氣校正和氣溶膠與水蒸氣的探測,因此,在本研究中將其剔除,選用其余10個波段的光譜信息進行生物量估測。植被指數是利用遙感技術對植被反射率和光譜特征進行計算,得到一個反映植被生長狀況的指標值。它是遙感技術在植被監測領域的一種應用。植被指數廣泛應用于農業、森林、水資源和環境等領域中的植被監測、植被變化檢測和生態環境評價等方面[19]。本研究對哨兵2號影像采用光譜和紋理2種指數進行分析。紋理指數是一種用于描述圖像中紋理信息的指數,通常應用于遙感影像分析。紋理可以被定義為圖像中局部灰度變化的表現,而紋理指數則可以用來描述不同尺度和方向的紋理信息,提取的特征如表1所示。

后向散射系數(Backscatter coefficient)是微波雷達反演森林生物量的重要參數之一,是衡量微波雷達信號穿過植被后被反射回來的強度大小指標,通常表示為[σ0]。對于具有不同結構和密度的植被,其后向散射系數也不同。因此,通過對不同類型植被的后向散射系數進行測量和分析,可以建立植被生物量與后向散射系數之間的關系模型,從而實現對生物量的估算。本研究參考Sun等[20]的研究結果,選用哨兵1號數據不同極化方式下的后向散射系數,以及其線性組合和紋理特征,參與后續生物量建模,提取的特征如表2所示。

研究區內人工刺槐林的健康狀況與土壤含水量、土壤含鹽量具有高度相關性[21]。此外,微地形與地下水含鹽量密切相關。為了探究這些關系,本研究基于數字地形數據提取地形因子,包括高程、坡度、坡向。

2.3.3 特征篩選 特征篩選在機器學習過程中至關重要,其目的是從大量特征中選擇最有用的子集,并用于模型的訓練。模型的預測性能很大程度上取決于所選的特征。特征數量過多會導致維數災難,使模型過度擬合訓練數據;但缺少重要特征則可能導致模型欠擬合。通過特征篩選,可以選擇具有代表性和重要性的特征并消除噪聲和冗余,從而提高模型的預測性能。

隨機森林模型可以通過對每個特征的貢獻度進行評估,從而進行特征重要性排序,根據重要性大小實現特征篩選的目的。在隨機森林模型中可以使用以下2種方法計算特征重要性。平均不純度減少法(Mean decrease impurity):對于每個特征,計算在隨機森林模型中所有決策樹上分裂該特征所帶來的平均不純度減少程度,即該特征所帶來的平均信息增益。然后將所有特征的平均不純度減少進行排序,得出特征重要性排名。平均精確度減少法(Mean decrease accuracy):對于每個特征,在隨機森林模型中進行隨機置換,打亂該特征的順序,然后計算置換后模型在測試集上的準確度下降程度。特征的重要性評估可以通過計算所有特征的平均精度下降程度得出。然后將所有特征的平均精度下降程度進行排序,得出特征重要性排名。2種方法各有優缺點,特征數較少且特征之間較少交互作用,適用于平均不純度減少法[22]。特征數較多或特征之間存在較多交互作用,則適用于平均精確度減少法。

由于本研究選取參與生物量模型構建的特征較多,因此采用平均精確度減少法進行特征篩選,并從中選取特征重要性排名前20的特征,用于后續森林生物量建模。篩選后的特征如表3所示。

2.3.4 生物量模型的建立與評價 基于多源遙感影像的生物量估測大多是利用野外實測采樣得到的生物量與光學影像提取的指數建立統計模型,然后外推來估測生物量。針對研究區內不同健康狀況的刺槐進行生物量估測時,傳統方法需要大量的野外實測樣本,且需包含不同健康狀況的刺槐,樣本量的多少直接影響生物量模型的精度[23]。有學者基于無人機激光雷達和背包移動激光雷達,結合野外實測數據對孤島不同健康等級的刺槐林進行生物量估測(R2為0.91~0.95)[15]。

本研究采用激光雷達預測的高精度孤島刺槐林生物量數據作為連接野外實測數據和多源遙感數據的橋梁,對研究區刺槐人工林生物量進行建模估測。利用4種模型建立多源遙感數據生物量估測模型,訓練樣本采用36 315個激光雷達估測的孤島生物量作為真實值進行模型訓練,驗證樣本采用72個野外實測生物量進行模型驗證,模型精度評價通過R2和RMSE進行評估。

回歸模型的評價指標通常包括均方根誤差(RMSE)、平均絕對誤差(Mean absolute error,MAE)、決定系數(R2)等[24]。本研究采用決定系數(R2)、均方根誤差(RMSE)作為評價指標,以評估所構建的森林生物量模型的性能。均方根誤差(RMSE)是一種用于回歸模型評價的指標,用于度量模型預測值與真實值之間的偏差。決定系數(R2)用于衡量模型對數據的擬合程度,解釋了目標變量方差的比例,取值范圍在0~1,R2=1表示模型完美擬合數據,而R2=0則表示模型不能解釋目標變量的任何變化。

3 結果與分析

3.1 森林生物量估算結果

由表4可知,綜合利用光學和微波雷達數據與地形數據進行生物量建模時,KNN模型、RF模型、XGBoost模型的R2分別為0.51、0.55、0.52,RMSE分別為14.88、13.85、14.72 t/hm2。采用Stacking模型進行集成學習可以進一步提升預測的準確度,R2=0.61、RMSE=13.42 t/hm2。為了直觀地展示各模型的預測表現,結合樣地生物量預測值與實際值繪制散點圖,由圖2可知,4種模型均存在生物量實際值高值被低估、低值被高估的現象,但這種現象在Stacking模型中有所改善。

3.2 森林生物量分布

本研究以前人基于激光雷達估測的36 315個高精度孤島刺槐生物量為訓練樣本,以2017年野外實測的72個數據為驗證數據,以多源遙感數據和數字地形數據提取的地形特征為模型自變量,建立生物量估測模型,估測研究區人工刺槐林生物量。比較KNN模型、RF模型、XGBoost模型和Stacking模型在生物量估算中的精度,結果表明,Stacking模型表現最優,因此選用該模型進行黃河三角洲人工刺槐林生物量估算??紤]到非參數模型具有一定隨機性,本研究采用5次預測結果的平均值進行生物量制圖,研究結果如圖3所示。

4個林場的生物量空間分布存在差異,其中孤島和大汶流自然保護區的生物量高于黃河故道和軍馬場2處刺槐林的生物量。孤島東北部生物量相對較低,主要是因為該區域的刺槐枯死情況比較嚴重,許多刺槐在生長未完全時就已經死亡,且刺槐林密度較低。黃河故道生物量普遍偏低的主要原因是由于20世紀90年代海水入侵和極端洪澇頻繁發生,導致土壤鹽漬化,嚴重影響了林地的健康生長[25]。此外,黃河故道上有農場進行水稻種植,水稻生長需要適合的土壤條件,為了調節土壤酸堿值,當地農民挖掘了灌溉渠,從附近河流引水泡田;該行為導致整個地區包括刺槐林遭受水淹,刺槐樹根淹沒在水中,缺氧容易導致根部腐爛,進而加速刺槐的死亡。

軍馬場與黃河故道的刺槐有相似的特征,主要受到人為活動的影響,也導致該地區刺槐的死亡率較高。前人的實地調查結果表明,當地居民對刺槐的砍伐和在刺槐林附近改種果樹、玉米、水稻等作物是導致刺槐死亡的重要原因之一[26]。大汶流自然保護區的整體生物量相較于黃河故道和軍馬場更高,這主要是由于該地區位于自然保護區,刺槐的管理和經營得到較好的保障,人為干擾較少,刺槐的主要影響因素是自然因素;此外,大汶流自然保護區刺槐的種植時間晚于其他3個林場[27]。

4 小結與討論

人工林是森林生態系統的重要組成部分,實現對人工林生物量連續空間分布及動態變化的快速、實時、精確估算是實現人工林可持續經營的前提,也是研究人工林碳循環的基礎。本研究基于多源遙感數據與數字地形數據,分別構建K-近鄰(KNN)模型、隨機森林(RF)模型、極值梯度增強(XGBboost)模型、Stacking模型,實現對黃河三角洲人工刺槐林生物量的估算,并得到對應的森林生物量空間分布。相較于K-近鄰(KNN)模型、隨機森林(RF)模型、極值梯度增強(XGBboost)模型,集成學習Stacking模型明顯提高了生物量估測的精度(R2=0.61、RMSE=13.42 t/hm2)。

孤島和大汶流自然保護區的生態環境相對穩定,森林生物量較高,軍馬場與黃河故道的森林生物量較低,這是由于軍馬場和黃河故道受到了人類農業活動的影響,例如灌溉和泡田壓鹽等。因此,這些區域的刺槐林生態環境相對脆弱。

參考文獻:

[1] 馬 欣.“雙碳”背景下零碳園區建設研究[J].合作經濟與科技,2023(8):4-7.

[2] 姜克雋,賀晨旻,莊 幸,等.我國能源活動CO2排放在2020—2022年之間達到峰值情景和可行性研究[J].氣候變化研究進展,2016,12(3):167-171.

[3] 劉菁婕.“雙碳”目標下公眾生態意識培養研究[J].湖北經濟學院學報(人文社會科學版),2023,20(3):28-32.

[4] BROWN S. Measuring carbon in forests: Current status and future challenges[J]. Environmental pollution, 2002, 116(3): 363-372.

[5] 方精云,陳安平.中國森林植被碳庫的動態變化及其意義[J].植物學報,2001(9):967-973.

[6] WANG H, PU R, ZHU Q, et al. Mapping health levels of Robinia pseudoacacia forests in the Yellow River Delta, China, using IKONOS and landsat 8 OLI imagery[J]. International journal of remote sensing, 2015, 36(4): 1114-1135.

[7] 陳幸良,巨 茜,林昆侖.中國人工林發展現狀、問題與對策[J].世界林業研究,2014,27(6):54-59.

[8] 李 蘭,陳爾學,李增元,等.合成孔徑雷達森林樹高和地上生物量估測研究進展[J].遙感技術與應用,2016,31(4):625-633.

[9] 張加龍,胥 輝.基于遙感的森林生物量估測樣地調查方法的研究動態[J].西南林業大學學報(自然科學),2019,39(4):166-173.

[10] LE TOAN T, BEAUDOIN A, RIOM J, et al. Relating forest biomass to SAR data[J]. IEEE transactions on geoscience and remote sensing, 1992, 30(2): 403-411.

[11] BEAUDOIN A, LE TOAN T, GOZE S, et al. Retrieval of forest biomass from SAR data[J]. International journal of remote sensing, 1994, 15(14): 2777-2796.

[12] ANTI E, PALOSCIA S, PETTINATO S, et al. The potential of multifrequency SAR images for estimating forest biomass in Mediterranean areas[J]. Remote sensing of environment, 2017, 200: 63-73.

[13]王長委,胡月明,沈德才,等.多源光學遙感數據估算桉樹森林生物量[J].測繪通報,2014(12):20-23.

[14] ZHENG D,RADEMACHER J, CHEN J, et al. Estimating aboveground biomass using Landsat 7 ETM+ data across a managed landscape in northern Wisconsin, USA[J]. Remote sensing of environment, 2004, 93(3): 402-411.

[15] LU J, WANG H, QIN S, et al. Estimation of aboveground biomass of Robinia pseudoacacia forest in the Yellow River Delta based on UAV and Backpack LiDAR point clouds[J]. International journal of applied earth observation and geoinformation,2020,86.DOI:10.1016/j.jag.2019.102014.

[16] 吳 彤,李 勇,葛 瑩,等.利用Stacking集成學習估算柑橘葉片氮含量[J].農業工程學報,2021,37(13):163-171.

[17] OGUNLEYE A, WANG Q G. XGBoost model for chronic kidney disease diagnosis[J]. IEEE/ACM transactions on computational biology and bioinformatics, 2019, 17(6): 2131-2140.

[18] BELGIU M, DRǎGU? L. Random forest in remote sensing: A review of applications and future directions[J]. ISPRS journal of photogrammetry and remote sensing, 2016, 114: 24-31.

[19] FOODY G M, BOYD D S, CUTLER M E J. Predictive relations of tropical forest biomass from landsat TM data and their transferability between regions[J]. Remote sensing of environment, 2003," " 85(4): 463-474.

[20] SUN G, RANSON K J, KHARUK V I. Radiometric slope correction for forest biomass estimation from SAR data in the Western Sayani Mountains, Siberia[J]. Remote sensing of environment, 2002, 79(2-3): 279-287.

[21] 宋 音,王 紅,路開宇,等.基于CCA方法的黃河三角洲不同健康刺槐林的土壤屬性研究[J].江西農業學報,2017,29(10):48-53.

[22] BIAU G, SCORNET E. A random forest guided tour[J]. Test, 2016, 25: 197-227.

[23] HYDE P, NELSON R, KIMES D, et al. Exploring LiDAR-RaDAR synergy—predicting aboveground biomass in a southwestern ponderosa pine forest using LiDAR, SAR and InSAR[J]. Remote sensing of environment, 2007, 106(1): 28-38.

[24] CHAI T, DRAXLER R R. Root mean square error (RMSE) or mean absolute error (MAE)?-Arguments against avoiding RMSE in the literature[J]. Geoscientific model development,2014,7(3):1247-1250.

[25] 趙 玉,王 紅,張珍珍.基于遙感光譜和空間變量隨機森林的黃河三角洲刺槐林健康等級分類[J].遙感技術與應用,2016,31(2):359-367.

[26] 張珍珍,王 紅.基于衛星IKONOS影像的黃河三角洲人工刺槐林健康狀況分類[J].科學技術與工程,2014,14(33):73-79.

[27] WANG H, ZHONG Y, PU R, et al. Dynamic analysis of Robinia pseudoacacia forest health levels from 1995 to 2013 in the Yellow River Delta, China using multitemporal landsat imagery[J]. International journal of remote sensing, 2018, 39(12): 4232-4253.

主站蜘蛛池模板: 久久久久久久97| 久久精品午夜视频| 综合网久久| 91网在线| 亚洲成人精品久久| a级毛片免费播放| 91久久夜色精品国产网站| 美女被躁出白浆视频播放| 色欲国产一区二区日韩欧美| 九九热免费在线视频| 2021国产乱人伦在线播放| 国产在线高清一级毛片| 国产91精选在线观看| 国产精品亚洲综合久久小说| 久久九九热视频| 亚洲欧美成人综合| 天堂成人在线| a毛片免费在线观看| 国产女人在线观看| 五月天久久婷婷| 国产精品开放后亚洲| 老司机精品一区在线视频| 中文字幕人成乱码熟女免费| jizz国产在线| 日本欧美一二三区色视频| 久久精品人妻中文系列| 久久综合色天堂av| 久久综合九九亚洲一区| 91麻豆国产在线| 国产一级裸网站| 91网在线| 波多野结衣中文字幕一区| 欧美日韩免费在线视频| 亚洲天堂高清| 国国产a国产片免费麻豆| 久久99国产乱子伦精品免| 暴力调教一区二区三区| 日韩欧美中文字幕在线韩免费| 她的性爱视频| 美女啪啪无遮挡| 国产成人亚洲日韩欧美电影| 久久精品无码一区二区日韩免费| 制服丝袜一区| 高潮毛片免费观看| 国产美女主播一级成人毛片| 在线看片中文字幕| 亚洲手机在线| 尤物精品视频一区二区三区| 国产国语一级毛片在线视频| 久久成人18免费| 毛片久久网站小视频| 久久久波多野结衣av一区二区| 亚洲性日韩精品一区二区| 欧美日本中文| 成人精品免费视频| 伊人五月丁香综合AⅤ| 人妻熟妇日韩AV在线播放| 99在线国产| 国产爽妇精品| 美女免费黄网站| 国产一级毛片在线| 中国国产A一级毛片| 人人妻人人澡人人爽欧美一区| 中文字幕乱妇无码AV在线| 国产亚洲精品自在久久不卡| 亚洲成a人片77777在线播放| 色九九视频| 欧美中文字幕一区| 亚洲人妖在线| 欧美中文字幕第一页线路一| 国产精品xxx| a级毛片毛片免费观看久潮| 国产精品观看视频免费完整版| 亚洲成A人V欧美综合天堂| 一区二区三区国产精品视频| 国产H片无码不卡在线视频| 免费不卡视频| 久久久久久久蜜桃| 亚洲AV成人一区二区三区AV| 自拍偷拍欧美日韩| 日韩在线视频网站| 国产激爽爽爽大片在线观看|