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

基于小波能量系數和葉面積指數的冬小麥生物量估算

2021-02-14 01:56:02李長春李亞聰王藝琳馬春艷陳偉男
農業機械學報 2021年12期
關鍵詞:模型

李長春 李亞聰 王藝琳 馬春艷 陳偉男 丁 凡

(河南理工大學測繪與國土信息工程學院, 焦作 454000)

0 引言

生物量是指單位面積內作物在一定時間內累積的有機物[1]。生物量作為作物生長發育過程中的一個重要生理參數,其變化能夠直接表征作物進行凈光合作用積累有機物的能力,反映作物的生長情況[2]。因此,科學、快速、準確地獲取生物量信息,對監測冬小麥長勢以及產量估算具有重要的意義。

傳統生物量獲取方法主要是野外實測方法,該方法應用于小區域地區精度高,但對于大區域而言,工作過程繁瑣、費時費力,時效性差,且大面積的田間破壞性取樣會影響作物的生長[3]。遙感技術憑借其宏觀、動態和實時等特點在作物監測方面得到廣泛的關注,為估算作物生物量提供有效手段。高光譜遙感技術具有較高的光譜分辨率和較強的波段連續性,是對地觀測領域的重大成就,是當今遙感科學的前沿技術[4]。在利用遙感技術進行作物生物量估算方面,構建植被指數基于經驗統計的方法是目前生物量估算比較常用的方法,植被指數是由一個或多個光譜波段以一定的數學方式進行結合,能夠反映植被生長狀況。孫奇等[5]通過分析可見光-近紅外植被指數、短波紅外植被指數和8個融合植被指數與生物量的關系,結果顯示融合植被指數對冬小麥生物量的估算能力更高,其中MTVI2結合NDMI估算精度最高。李嵐濤等[6]利用波長400~950 nm處的冠層反射率,采用支持向量機和偏最小二乘回歸方法構建不同生育期的冬小麥生物量估算模型,結果顯示抽穗期效果最好,R2大于0.85。郭超凡等[7]利用單變量植被指數,基于多元線性和隨機森林方法對小麥地上部生物量進行估算,發現綠色葉綠素指數和歸一化差異水體指數與生物量具有最優擬合關系,隨機森林模型估算效果最優,R2達到0.74。陶惠林等[8]基于植被指數、紅邊參數結合多元線性回歸模型估算冬小麥生物量,結果表明植被指數耦合紅邊參數能夠提高模型估算精度。使用植被指數建立生物量估算模型方便簡單,但是光譜信息利用率較低,僅使用了其中的幾個波段,不能充分利用有效光譜信息,并且使用植被指數估算多生育期的生物量容易存在飽和現象。近年來,為了進一步提高生物量估算精度,一些學者通過對光譜進行變換,如小波變換[9-11]、光譜微分[12-14]、連續統去除法[15-16]等多種方法挖掘潛在的有效光譜信息,提取隱藏信息,增強光譜信息的敏感度,有效解決了模型的飽和問題。

然而,這些研究僅僅利用光譜及其變換信息估算作物生物量,未綜合考慮葉面積指數與生物量的關系。本文基于高光譜數據,通過連續小波變換得到的小波能量系數,耦合葉面積指數,分別利用支持向量回歸算法、隨機森林算法和高斯過程回歸算法,構建冬小麥生物量估算模型,分析不同模型生物量估算效果,尋求最佳估算方法,以期為基于遙感技術的作物高通量表型參數快速估算提供參考。

1 材料與方法

1.1 研究區概況

研究區位于北京市昌平區小湯山國家精準農業研究示范基地。基地位于北緯40°00′~40°21′,東經116°34′~117°00′,平均海拔為36 m,氣候類型是典型的暖溫帶半濕潤大陸季風氣候,春季干旱多風,夏季炎熱多雨,秋季涼爽,冬季寒冷干燥,四季分明。一般為每年9月末至10月初播種冬小麥,次年3月左右返青,經過拔節、孕穗、開花、灌漿4個生育期,6月中上旬收獲。研究區共設置48個小區,每個小區面積為48 m2(6 m×8 m),每16個小區為1個重復,共3個重復。設置4個氮肥水平:N1:0 kg/hm2,N2:15 kg/hm2,N3:30 kg/hm2,N4:45 kg/hm2。冬小麥管理模式參考一般的大田模式,具有優越的肥力、灌溉和排水條件。其他按照田間實際管理進行操作,研究區地理位置如圖1所示。

1.2 數據獲取

從4月13日—5月27日,每隔13~15 d,在實驗區里選取樣本點進行小麥冠層高光譜數據、葉面積指數和生物量測量。

1.2.1冠層高光譜數據獲取

采用美國ASD(Analytica Spectra Deviecs,Inc)生產的FieldSpec便攜式地物光譜儀獲取冬小麥冠層高光譜數據。該光譜儀可測量的波長為350~2 500 nm,采樣間隔為1.4 nm(350~1 000 nm)和2 nm(1 000~2 500 nm),儀器內部重采樣后可達到1 nm,測量時間選擇在每天10:30—14:00。為減少土壤背景的影響,測量過程中,傳感器探頭垂直向下,距地表1.3 m處,視場角為25°。每個實驗小區均勻采集20條光譜數據,采集完成后利用與ASD光譜儀配套的光譜數據處理軟件ViewSpecPro導出無量綱的反射率,以每個實驗小區20條光譜數據的平均值作為各小區的冠層反射率。由于水分吸收帶嚴重,所以選取波段350~1 350 nm內的光譜數據。在測量過程中,當每個樣點測量完成后,應及時利用白板進行標準白板校準(標準白板反射率為1,獲得的目標光譜為相對反射率)。

1.2.2葉面積指數獲取

采集光譜數據的同時,在采樣點上獲取小麥植株,帶回室內,將取樣的冬小麥進行莖葉分離,利用CI-203型激光葉面積儀測定葉片葉面積,求其總和,進而計算出單莖葉面積,然后乘以單位土地面積單莖數得到冬小麥葉面積指數。4個生育期共獲取到192個冬小麥葉面積指數樣本數據。

1.2.3生物量獲取

利用收獲法獲取生物量數據。每個小區內,隨機選取10株長勢均勻冬小麥植株作為樣本,經莖葉分離、稱量等操作后編號放入紙袋中,將干燥溫度保持在105℃,運行30 min殺青,隨后將溫度設置為75℃干燥至質量恒定(48 h左右),測定各器官干質量,最后根據群體密度及樣本干質量,計算單位面積冬小麥地上生物量,最終獲取生物量。

1.3 數據處理與模型構建方法

1.3.1連續小波變換

小波變換是由傅里葉變換發展起來的一種數據處理方法,通過伸縮和平移等運算功能對函數或信號進行多尺度細化分析[17]。小波變換包括連續小波變換(Continuous wavelet transform,CWT)和離散小波變換(Discrete wavelet wransform,DWT),后者是前者的離散化。連續小波變換作為一種線性變換的方法,廣泛應用于高光譜數據處理[18-20]。本文采用CWT對ASD高光譜數據進行分解,其原理就是將冠層高光譜數據分解成不同尺度和不同波長的小波能量系數。將一維空間的高光譜反射率轉換為二維空間尺度和小波能量系數。利用CWT分別獲取4個生育期前10尺度下的小波能量系數,計算式為

(1)

(2)

式中f(λ)——冠層高光譜反射率

λ——波段350~1 350 nm范圍的光譜波段

Ψa,b——小波基函數

a——尺度因子b——平移因子

1.3.2模型構建方法

(1)高斯過程回歸

高斯過程回歸(Gaussian process regression,GPR)由RASMUSSEN和WILLIAMS在2005年提出[21],是將統計學習理論和貝葉斯理論相結合的一種機器學習方法,對于處理小樣本、高維度和非線性復雜關系的回歸和分類問題都有較好的效果,近年來成為作物參數估算最高效的機器學習算法之一[22-23]。

(2)支持向量回歸

支持向量回歸(Support vector regression,SVR)算法是由支持向量機算法演化而來,作為一種監督學習的機器學習算法,基本原理是通過建立一個最優決策超平面,使距樣本之間距離達到最小化,從而擬合樣本數據[24],與其他算法相比,模型準確性較高,對高維度、小樣本數據有較好的處理能力,具有良好的泛化能力和魯棒性[25]。

(3)隨機森林

隨機森林(Random forest,RF)基本原理是利用多個決策樹算法對相同現象做重復的預測[26],作為一種組成式的有監督學習方法,是機器學習方法的分支之一,即集成學習方法,在處理大量輸入變量問題中優勢明顯,預測性能好、魯棒性強[27]。

本文以不同生育期小波能量系數、小波能量系數耦合LAI為自變量,生物量為因變量,在每個生育期,分別挑選2/3樣本數據(32個)作為建模集,1/3樣本數據(16個)作為驗證集,構建冬小麥生物量估算模型。

1.3.3相關系數

采用皮爾遜相關系數R來表示兩個隨機變量之間的相關程度,相關系數取值范圍為[-1,1],相關系數絕對值越大,相關性越強,相關系數絕對值越接近于0,相關度越弱[28]。相關系數計算式為

(3)

式中R(X,Y)——兩隨機變量相關系數

cov(X,Y)——兩隨機變量協方差

σ——標準差

1.3.4模型精度評價

采用決定系數(R2)、均方根誤差(RMSE)和標準均方根誤差(nRMSE)3個指標作為模型精度評價指標。R2越大,表明該模型的擬合效果越高,并且越穩定。nRMSE小于10%表示估算值和實測值一致性極好;nRMSE為10%~20%表示一致性較好;nRMSE為20%~30%表示一致性中等;nRMSE大于等于30%表示一致性較差。

2 結果與分析

2.1 冬小麥生物量與葉面積指數變化特征

為了研究冬小麥生物量與葉面積指數之間的關系,對4個生育期生物量和葉面積指數的結果進行統計,結果如表1所示。

表1 不同生育期實測生物量和葉面積指數統計分析Tab.1 Statistical analysis of measured biomass and leaf area index in different growth periods

由表1可知,從拔節期到灌漿期,冬小麥生物量總體呈先升高后下降的變化趨勢,拔節期生物量最小,平均值為2631.90 kg/hm2,開花期達到最大,平均值為5940.82 kg/hm2,開花期后生物量逐漸下降。變異系數在26.28%~31.87%之間,離散程度較穩定,表明在完整生育期內,冬小麥生物量值域分布隨著時間推移趨于集中。從拔節期到開花期,總生物量迅速上升,隨著葉片數量增加和葉面積擴大,光合速率提高,使得各構件生物量迅速積累;開花期到灌漿期,即成熟初期,由于葉片和莖稈逐漸老化變黃,光合作用能力下降,生成的光合產物低于自身消耗,乃至脫落,小麥生物量呈很小的下降趨勢,與侯學會等[29]的研究結果一致。

冬小麥葉面積指數在整個生育期內先增大后減小,呈拋物線型變化。在孕穗期之前,冬小麥營養生長和生殖生長較為活躍,葉片快速生長,葉面積指數迅速增加,且在拔節期與孕穗期之間達到峰值,孕穗期之后,隨著營養向生殖器官的快速轉移,葉片逐漸老化變黃、脫落,因而葉面積指數呈現下降趨勢。拔節期離散程度較小,變異系數為25.35%,而灌漿期葉面積指數離散程度較大,變異系數為55.58%。

2.2 相關性分析

2.2.1原始光譜與生物量相關性分析

對冬小麥拔節期、孕穗期、開花期和灌漿期去噪后的高光譜數據與相應的生物量數據進行相關性分析,計算相關系數R,繪制不同生育期相關系數曲線,結果如圖2所示。

由圖2可知,各個生育期的冬小麥冠層光譜反射率與生物量的相關系數變化趨勢表現出相似的特點。在波段500~700 nm附近,冬小麥冠層光譜反射率與生物量都有很高的負相關性,波段800~900 nm范圍內相關系數基本保持不變,冬小麥冠層光譜反射率與生物量都有很高的正相關性。在拔節期,波長691 nm處負相關程度最大,相關系數R為-0.71,波長823 nm處正相關程度最大,相關系數R為0.62;在孕穗期,波長695 nm處負相關程度最大,相關系數R為-0.74,波長868 nm處正相關性最大,相關系數R為0.64;在開花期,波長691 nm處負相關性最大,相關系數R為-0.73,波長819 nm處正相關性最大,相關系數R為0.84;在灌漿期,波長505 nm處負相關性最大,相關系數R為-0.80,波長877 nm處正相關性最大,相關系數R為0.85。

2.2.2小波能量系數與生物量相關性分析

對于冬小麥4個生育期,利用連續小波變換對冬小麥冠層高光譜數據進行1~10尺度分解,不同波段數與不同分解尺度對應的相關系數,表征小波能量系數與冬小麥生物量之間的相關性,通過連續小波變換間接對目標光譜信息進行優化。在每個分解尺度下,各波段生成特定的小波能量系數,然后分別與生物量進行相關性分析,生成相關性矩陣,結果如圖3所示。

由圖3看出,4個生育期小波能量系數與生物量的相關性隨著分解尺度的增加,呈現先升高后降低的趨勢;通過0.01極顯著水平檢驗的光譜波段數量不斷增加,當分解尺度為10,4個生育期均達到1 001個波段。由圖3a可知,在拔節期當分解尺度為1~8時,小波能量系數與生物量的最大相關系數絕對值均在0.71以上,分解尺度為2時,在波長948 nm處|R|最大可達0.77;當分解尺度為9和10時,|R|最大分別為0.64、0.52,小波能量系數和生物量相關性較小,分解得到的小波能量系數已經無法表征冬小麥的生物量信息。由圖3b可知,在孕穗期當分解尺度為1~8時,小波能量系數與生物量的最大相關系數絕對值均在0.80以上,二者呈高度相關,分解尺度為2時,在波長768 nm處|R|最大可達0.84;當分解尺度為9和10時,|R|最大分別為0.64、0.44,小波能量系數和生物量中低度相關,無法表征冬小麥的生物量信息。由圖3c可知,在開花期當分解尺度為1~8時,小波能量系數與生物量的最大相關系數絕對值均在0.89以上,二者呈高度相關,當分解尺度為3時,在912 nm處|R|最大可達0.91;當分解尺度為9和10時,|R|最大分別為0.74、0.64,小波能量系數和生物量中度相關,不能充分表征冬小麥的生物量信息。由圖3d可知,在灌漿期當分解尺度為1~9時,小波能量系數與生物量的最大相關系數絕對值均在0.86以上,二者呈高度相關,分解尺度為7時,在波長1 001 nm處|R|最大可達0.88;當分解尺度為10時,|R|最大為0.67,小波能量系數和生物量中度相關,不能充分表征冬小麥的生物量信息。

2.2.3葉面積指數與生物量相關性分析

將實測得的葉面積指數與生物量數據利用式(3)計算相關系數。4個生育期葉面積指數與生物量相關系數有所差異,4個生育期|R|分別為0.91、0.82、0.93和0.87,均大于0.8,且達到0.01顯著性水平,說明葉面積指數與生物量具有很高的相關性。因此,在各生育期生物量預測模型中,可將葉面積指數與小波能量系數耦合共同作為模型輸入變量以提高冬小麥生物量估算的精度。

2.3 敏感波段篩選

利用圖3對每個生育期生物量與小波能量系數進行相關性分析,篩選出相關系數絕對值較高的前10個波段作為輸入變量,估算冬小麥生物量,篩選結果如表2所示。

由表2可知,不同生育期不同尺度下最優波段有所不同,在拔節期,最優波段尺度在1、2、3、4、6、7;在孕穗期,最優波段尺度在1、2、3、4、5、6;在開花期,最優波段尺度在3、4、5、6;在灌漿期,最優波段尺度在4、5、6、7、8。綜合分析來看,4個生育期中,最優波段尺度主要集中在3、4、5、6,并且每種尺度下篩選的敏感波段與生物量均達到0.01顯著水平。

2.4 模型構建與精度評價

基于以上相關性分析結果,每個生育期以冬小麥生物量作為因變量,以小波能量系數、小波能量系能數耦合LAI為自變量,分別采用支持向量回歸算法(SVR)、隨機森林算法(RF)、高斯過程回歸(GPR)算法構建冬小麥4個生育期的生物量高光譜估算模型,并且進行驗證,結果見表3、4。

由表3、4可知:

(1)利用相關性較好的小波能量系數為輸入變量,分別利用支持向量回歸算法、隨機森林回歸算法、高斯過程回歸算法進行冬小麥生物量估算。拔節期和灌漿期利用隨機森林算法進行生物量估算效果最好,拔節期模型驗證R2、RMSE、nRMSE分別為0.40、505.83 kg/hm2、19.74%,灌漿期模型驗證R2、RMSE、nRMSE分別為0.89、954.50 kg/hm2、19.67%。孕穗期和開花期利用高斯過程回歸算法對生物量估算效果最佳,孕穗期模型驗證R2、RMSE、nRMSE分別為0.83、974.78 kg/hm2、21.69%,開花期模型驗證R2、RMSE、nRMSE分別為0.93、537.82 kg/hm2、9.27%。綜合考慮模型建立和模型驗證的評價指標結果以及模型的估算能力,開花期利用高斯過程回歸算法構建模型具有最優估算精度。

表2 前10敏感波段及對應尺度Tab.2 Top ten sensitive bands and corresponding scales

表3 不同生育期不同方法冬小麥生物量的建模精度Tab.3 Modeling results of winter wheat biomass in different growth periods and different methods

表4 不同生育期不同方法估算冬小麥生物量的驗證精度Tab.4 Validation results of different growth periods and different methods to estimate winter wheat biomass

(2)將葉面積指數信息融合相關性較好的小波能量系數,分別利用支持向量回歸算法、隨機森林回歸算法、高斯過程回歸算法進行冬小麥生物量估算,模型建立和模型驗證散點圖如圖4所示。4個生育期均利用高斯過程回歸算法進行生物量估算效果最好,拔節期模型驗證R2、RMSE、nRMSE為0.77、363.73 kg/hm2、13.87%,孕穗期模型驗證R2、RMSE、nRMSE為0.84、669.61 kg/hm2、13.86%,開花期模型驗證R2、RMSE、nRMSE為0.94、382.32 kg/hm2、6.34%,灌漿期模型驗證R2、RMSE、nRMSE分別為0.91、613.97 kg/hm2、11.61%。綜合考慮模型建立和模型驗證的評價指標結果以及模型的估算能力,表明在開花期利用高斯過程回歸構建的模型具有最優估算精度。通過對比分析,葉面積指數耦合小波能量系數估算生物量,在一定程度上能夠克服冠層光譜飽和現象,提高模型估算精度。

(3)與僅僅基于小波能量系數相比,耦合葉面積指數數據,利用支持向量回歸、隨機森林回歸、高斯過程回歸方法分析進行生物量估算時,拔節期驗證的R2分別提高0.21、0.33、0.38,孕穗期驗證的R2分別提高0.01、0.02、0.01,開花期驗證的R2分別提高0.03、0.02、0.01,拔節期驗證的R2分別提高0.01、0.01、0.06,表明小波能量系數耦合葉面積指數數據,可以提高生物量估算精度。

3 討論

將4個生育期原始光譜、小波能量系數分別與生物量進行相關性分析,經對比可知,原始光譜和小波能量系數與生物量的相關性絕對值均呈先升高后降低的變化趨勢,開花期的相關系數絕對值最大。可能的原因是從拔節期到開花期,葉片數量增加和葉面積擴大,光合速率大大提高,使得各構件生物量迅速積累和增加,冬小麥生物量在開花期達到最大,使得原始光譜和小波能量系數與生物量的關系由弱變強。開花期到灌漿期,即成熟初期,葉片和莖稈逐漸老化變黃,植被覆蓋度明顯降低,光合作用能力下降,生成的光合產物低于自身消耗,乃至脫落,小麥生物量下降,使得上述2個變量與生物量的相關性由強變弱。

每個生育期提取的小波能量系數與生物量的相關性都高于原始高光譜與生物量的相關性,小波能量系數與生物量的相關系數絕對值最大值較原始光譜與生物量的相關系數絕對值最大值分別提高了0.06、0.1、0.07和0.03,這是因為高光譜數據采集過程中,由于受土壤背景、環境因素等影響,導致原始光譜存在噪聲,影響敏感光譜信息提取,利用小波變換技術可以細化光譜信息,深度挖掘潛在的敏感光譜信息,凸顯出更多與生物量相關的有效信息,使得篩選出的小波能量系數與生物量的相關性更高。

利用原始高光譜數據估算生物量的困難在于提取原始高光譜數據中對生物量較敏感的光譜波段,以往基于高光譜數據的生物量估算模型大多基于植被指數。本研究采用SVR、RF和GPR 3種算法結合小波能量系數、小波能量系數耦合葉面積指數構建生物量估算模型,實現對冬小麥生物量的估算。研究結果顯示,4個生育期利用小波能量系數耦合葉面積指數基于GPR算法構建的模型對生物量估算效果最好,建模和驗證精度均較高,4個生育期建模R2分別為0.86、0.89、0.93、0.91,驗證R2分別為0.77、0.84、0.94、0.91,這與增加葉面積指數作為輸入變量有關,葉面積指數包含作物冠層結構中非常重要的信息,與生物量高度相關,是決定作物光合作用強弱的重要變量之一,研究結果表明耦合葉面積指數可以提高生物量估算精度,與文獻[30-32]研究結果相一致。另外,利用GPR算法估算冬小麥生物量精度較高,原因可能是GPR作為一種概率方法,訓練模型通過擬合均值和協方差函數找到所有訓練數據的單獨函數,與其他算法學習方法相比,更適合訓練小樣本模型,且受生育期的影響較其他兩種算法小。

4 結束語

基于小波能量系數,同時耦合葉面積指數數據,利用支持向量回歸、隨機森林回歸、高斯過程回歸方法,構建冬小麥生物量估算模型,并對模型精度進行驗證。結果表明:4個生育期,開花期生物量估算精度最高;利用支持向量回歸、隨機森林回歸、高斯過程回歸方法估算生物量,高斯過程回歸算法估算生物量效果最好;小波能量系數耦合葉面積指數與僅使用葉面積指數為模型因子相比,生物量估算精度得到提高,構建的驗證模型擬合性和穩定性都較好。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲swag精品自拍一区| 久草网视频在线| 国产精品流白浆在线观看| 88av在线播放| 精品欧美日韩国产日漫一区不卡| 精品视频一区二区三区在线播| 国产精品免费电影| 亚洲综合香蕉| 中文字幕日韩视频欧美一区| 高清久久精品亚洲日韩Av| 尤物国产在线| 欧美伦理一区| 极品私人尤物在线精品首页| 91色在线观看| 国产精品自拍露脸视频| 高潮毛片无遮挡高清视频播放| 日本欧美一二三区色视频| 奇米精品一区二区三区在线观看| 五月婷婷综合网| 亚洲精品国产乱码不卡| 欧美无专区| 中文字幕日韩丝袜一区| 亚洲精品在线影院| 天天躁日日躁狠狠躁中文字幕| 2048国产精品原创综合在线| 亚洲中文在线视频| 亚洲无码精品在线播放| 亚洲aⅴ天堂| 丁香六月综合网| 无码国产伊人| 四虎亚洲精品| 538国产视频| 凹凸国产分类在线观看| 热久久综合这里只有精品电影| 免费无码AV片在线观看中文| 国产迷奸在线看| 波多野结衣一区二区三区四区视频 | 国产福利小视频在线播放观看| 国产xx在线观看| 国内嫩模私拍精品视频| 超碰91免费人妻| 国产小视频免费| 欧美亚洲日韩中文| 亚洲男人在线| 青青国产成人免费精品视频| 久久久精品久久久久三级| 免费啪啪网址| 99热最新在线| 亚洲色欲色欲www网| 亚洲成人在线免费| 久久久波多野结衣av一区二区| 国产网站免费观看| 午夜老司机永久免费看片 | 99re这里只有国产中文精品国产精品 | 久久人午夜亚洲精品无码区| 97超级碰碰碰碰精品| 中文字幕久久波多野结衣| 亚洲精品波多野结衣| 精品国产成人a在线观看| 欧美精品一区在线看| 亚洲三级电影在线播放| 日韩精品久久久久久久电影蜜臀| 97se亚洲综合在线| www.狠狠| 欧美第一页在线| 91国语视频| 91免费在线看| 国产精品55夜色66夜色| 亚洲国产成人久久精品软件| 亚洲经典在线中文字幕| 免费又黄又爽又猛大片午夜| 免费A级毛片无码免费视频| 天天色天天综合网| P尤物久久99国产综合精品| 欧美影院久久| 中文字幕啪啪| 久久精品只有这里有| 99re热精品视频国产免费| 久久不卡国产精品无码| 国产视频 第一页| 成人日韩欧美| 999精品视频在线|