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

無人機高光譜遙感估算冬小麥葉面積指數

2021-01-19 04:59:48陳曉凱李粉玲王玉娜史博太侯玉昊常慶瑞
農業工程學報 2020年22期
關鍵詞:模型

陳曉凱,李粉玲,2※,王玉娜,史博太,侯玉昊,常慶瑞,2

(1.西北農林科技大學資源環境學院,楊凌 712100;2.農業部西北植物營養與農業環境重點實驗室,楊凌 712100)

0 引 言

葉面積指數(Leaf Area Index,LAI)被定義為作物在單位土地面積上的葉片的面積之和[1]。快速、無損、精準地監測冬小麥關鍵生育期的葉面積指數是準確掌握作物冠層結構、長勢信息、地上生物量、產量以及病蟲害監測等的重要途徑,對田間生產管理具有重要意義[2]。遙感技術能以較低的成本獲得 LAI的時空變化信息,克服了傳統田間地面測量只能獲取單點數據的不足。相對于衛星遙感技術,無人機遙感技術體積小、機動靈活、成本低、時空分辨率高,逐漸成為小范圍內農業遙感作物參數反演的重要手段[3-4],在水稻、小麥、大豆、棉花等作物的葉綠素、花青素、葉面積指數等理化參數的遙感估算中取得了較好的效果[5-8]。

利用無人機高光譜進行作物理化參數遙感估算主要是通過構建光譜指數,包括典型植被指數和任意兩波段組合的窄波段光譜指數,利用回歸分析方法進行估算模型構建。對于典型植被指數,高林等[9]以比值植被指數、歸一化植被指數、土壤調整植被指數、差值植被指數、三角植被指數等 5種傳統植被指數,采用經驗模型法對大豆葉面積指數進行反演,結果表明歸一化植被指數線性回歸模型精度最高;田明璐等[10]基于比值植被指數、差值植被指數、歸一化植被指數和綠波段歸一化植被指數構建極值植被指數,研究顯示以多個極值植被指數參與構建的偏最小二乘回歸模型對棉花的葉面積指數預測精度最高。而大多研究表明通過任意兩波段組合構建窄波段光譜指數更能挖掘出作物理化參數的敏感位置,顯著提高模型的預測精度[11-12]。王偉東等[13]通過構建任意兩波段組合的比值光譜指數、差值光譜指數、歸一化光譜指數分析了冬小麥花青素的空間分布,結果表明以窄波段光譜指數為自變量的花青素估算模型精度顯著高于典型植被指數模型;秦占飛等[14]通過任意兩波段組合的比值光譜指數、二次修正土壤調節光譜指數、差值光譜指數、歸一化光譜指數探尋反演水稻LAI的最優波段組合,結果表明以比值光譜指數為自變量建立的水稻LAI估算模型精度最高,但這些窄光譜指數均基于原始冠層光譜提取。姚霞等[15]對地面冠層高光譜進行數學變換并探索最佳波段組合對冬小麥氮素含量的響應能力,結果表明一階導數變換下的任意兩波段組合比值光譜指數、歸一化光譜指數、土壤調節光譜指數所構建的模型較原始光譜模型的決定系數明顯提高;李粉玲等[16-17]研究表明對冠層原始光譜進行一階導數變換、吸光度變換和連續統去除變換均能不同程度的提高氮素的遙感估算效果。而基于變換光譜提取窄波段光譜指數能否提高LAI的預測精度值得探討。

綜上,本研究以關中平原冬小麥為研究對象,以無人機高光譜成像儀影像數據為基礎,系統構建基于原始光譜、一階導數光譜和連續統去除光譜的任意兩波段間的差值光譜指數(Difference Spectral Index,DSI)、比值光譜指數(Ratio Spectral Index,RSI)和歸一化光譜指數(Normalized Spectral Index,NDSI),探討基于變換光譜窄波段光譜指數的LAI估算能力,以期為冬小麥的田間管理提供科學的技術指導。

1 材料與方法

1.1 試驗設計

冬小麥田間試驗布置在陜西省乾縣梁山鄉齊南村(108°07′E,34°38′N),地處陜北黃土高原南緣與關中平原的交界地帶,屬溫帶大陸性季風氣候,以種植冬小麥和夏玉米為主。本試驗共設置36個小區試驗和 4個大田試驗,供試品種均為當地常見的小偃22號。36個小區面積均為9 m×10 m,小區試驗進行氮(N)、磷(P)、鉀(K)3種肥料的脅迫試驗,分別設置 6個水平的施肥處理,分別為N(0、30、60、90、120、150 kg/hm2)、P (0、22.5、45、67.5、90、112.5 kg/hm2)和 K(0、22.5、45、67.5、90、112.5 kg/hm2),每個處理設置重復2次。4個大田面積均為480 m2,設置4個水平的氮(N)施肥處理0、60、120、180 kg/hm2,小區和大田的管理方式均與當地常規冬小麥相同。本研究在拔節期(2018年3月29日)進行冬小麥冠層LAI和無人機高光譜影像獲取。

1.2 數據采集與預處理

1.2.1 冬小麥葉面積指數(LAI)采集

冠層分析儀SunScan(Delta Co.,英國)能同時獲得直射光合有效輻射(Photosynthetically Active Radiation,PAR)、散射PAR、透過PAR和LAI等多種指標。其中LAI是通過傳感器測量冠層的輻射值和橢球體葉傾角分布參數,利用冠層投射光學模型的方法獲取,具有無損性和可重復性觀測的特點[18]。研究表明SunScan冠層分析儀測定的葉面積指數與傳統方法測定的葉面積指數有很好的相關性,在多種作物的LAI研究中得到廣泛應用[19-20]。本研究基于 SunScan冠層分析儀,在各小區和大田地塊分別選取1~2個樣點進行冠層LAI測定。測量時,每個樣點分別在東、西、南、北 4個方位上測定 4次,取其平均值作為該樣點上的LAI值。測量LAI的同時用差分全球定位系統(Global Positioning System,GPS)記錄該樣點的坐標信息。共獲取LAI有效樣本數據66個,將LAI值按從小到大的順序進行排序,按照2∶1的比例抽樣獲取建模集(44個)和驗證集(22個),樣本LAI值統計特征如表1所示。

圖1 研究區位置及采樣點分布Fig.1 Location of study area and distribution of sampling sites

表1 拔節期冬小麥葉面積指數特征統計Table 1 Feature statistics of winter wheat leaf area index at jointing stage

1.2.2 無人機高光譜影像獲取

試驗采用八旋翼無人機搭載 Cubert UHD185(簡稱UHD185)高光譜成像儀在地面LAI測定的同時獲取冬小麥近地冠層高光譜影像。UHD185成像儀的光譜覆蓋范圍為450~950 nm,光譜分辨率和采樣間隔均為4 nm,共125個探測通道。數據獲取之前,首先使用參考板對成像光譜儀進行輻射校正,并規劃飛行航線,設定1 ms采樣間隔進行數據采集,飛行高度為 100 m,飛行速度為6 m/s,光譜儀鏡頭垂直向下,視場角 30°,并保證 80%的航向重疊和 60%的旁向重疊,獲取的高光譜影像地面分辨率為0.3 m,每幅影像幅寬約16 m。在無人機影像專業處理軟件Cube-Pilot1.4.4下對高光譜影像進行拼接,依據Google Earth 影像和地面控制點在地理信息系統系列軟件ArcGIS 10.6下對研究區拼接影像進行地理配準;在遙感圖像處理軟件 ENVI(Environment for Visualizing Images)中采用 3×3均值濾波對配準后的圖像進行卷積運算,提取3×3卷積核范圍內的平均值作為中心像元值。依據地面上66個LAI觀測樣點位置,提取樣點所對應的高光譜反射率。

1.2.3 冠層高光譜變換

在獲取高光譜數據時,由于光子效應等原因,難免會受到噪聲的干擾,在遙感圖像處理軟件ENVI中對無人機高光譜數據進行平滑濾波處理,剔除依附在冠層高光譜影像上的噪聲信息,得到 66個樣本點的原始光譜(Original Spectrum,OS)。對去噪光譜進行一階導數變換和連續統去除變換,變換后的冠層光譜分別標記為一階導數光譜(First Derivative Spectrum,FDS)、連續統去除光譜(Continuum Removal Spectrum,CRS)。

1.2.4冠層高光譜窄波段光譜指數

高光譜提供了豐富的光譜信息,針對高光譜任意波段進行兩兩組合時,便可構建高光譜窄波段光譜指數。為了便于研究作物冠層高光譜窄波段光譜指數與LAI的相關關系,本研究在450~950 nm的波段范圍內,選取了計算最為簡單且最常用的歸一化植被指數、差值植被指數和比值植被指數,系統構造了原始光譜、一階導數光譜和連續統去除光譜的任意兩波段間的差值光譜指數(Difference Spectral Index,DSI)、比值光譜指數(Ratio Spectral Index,RSI)和歸一化光譜指數(Normalized Spectral Index,NDSI),其計算如式(1)~式(3)所示:

式中i、j為高光譜波長,nm;Ri、Rj分別表示i、j波長所對應的高光譜反射率(本研究分別取原始光譜、一階導數光譜和連續統去除光譜)。

1.3 模型構建與精度檢驗

1.3.1 模型構建

本研究基于450~950 nm波段范圍內的無人機冠層高光譜原始光譜反射率、一階導數光譜和連續統去除光譜分別提取3種光譜指數(差值光譜指數、比值光譜指數和歸一化光譜指數),分析任意兩波段交叉提取的光譜指數與LAI的相關性,篩選最優的光譜指數進行拔節期冬小麥LAI的估算模型研究。采用的建模方法包括一元回歸(Simple Regression,SR)、多元線性回歸(Multiple Linear Regression,MLR)、BP神經網絡(Back Propagation Neural Network,BP)和隨機森林回歸(Random Forest Regression,RFR)。其中 BP神經網絡和隨機森林回歸作為機器學習算法,在農作物生理生化參數的定量估算中表現突出[4,11,21-24]。

BP神經網絡模型由一個輸入層、至少一個隱含層和一個輸出層構成,是一種基于梯度下降思維、采用反向傳播算法不斷調節隱含層神經元個數對數據進行反復訓練獲取最優模型參數的機器學習算法[25]。本研究首先對建模集和驗證集數據進行歸一化處理,在MATLAB 2017a中調用BP神經網絡工具箱訓練數據,經過反復調試參數對數據進行訓練,最終確定隱含層節點數為4,輸出層節點數為1時獲取到最佳模型。分別將建模集和驗證集輸入模型得出各自預測值并輸出,再將輸出結果反歸一化完成BP神經網絡模型構建。

Breiman[26]在 2001年提出隨機森林回歸,這是一種基于分類樹的回歸算法。其實質是將樣本不斷放回并進行多次取樣形成訓練集,通過決策樹的組合對預測結果求平均來進行預測,這種算法主要考慮決策樹數量以及分割結點的數量這2個參數。本研究在統計軟件R中調用RandomForsest程序包進行隨機森林回歸分析。通常情況下,隨著決策樹數目的增加,預測模型的誤差會逐漸下降并趨于平穩。本研究通過繪制決策樹數目與模型預測誤差之間的曲線圖來確定最優決策樹數目。經過不斷調參,最終確定決策樹數量為500,分割節點數量為1[27-28]。

1.3.2 模型精度的評價指標

本研究選取決定系數(coefficient of determination,R2)、均方根誤差(Root Mean Square Error,RMSE)和相對預測偏差(Residual Prediction Deviation,RPD)進行模型精度檢驗。R2表示模擬值與實測值的擬合程度,R2越高越好,越接近于1,表明模型的估測效果越好;RMSE反映了模擬值與實測值的偏離度,RMSE越低,表明在該模型下的預測值與實測值越接近,模型效果越好,Rc2、RMSEc表示建模集決定系數和均方根誤差;Rp2、RMSEp表示驗證集決定系數和均方根誤差;RPD是通過衡量預測值和實測值的相對偏差程度來評價一個模型預測能力的重要指標。通常認為RPD<1.5 模型不具備預測能力;1.5<RPD<2.0模型可以粗略的估測樣本,預測能力尚可;RPD>2.0 表明模型具有極好的預測能力。RPD計算如式(4)所示[29]:

式中SD為驗證集標準偏差,RMSEp為驗證集均方根誤差。

2 結果與分析

2.1 窄波段光譜指數與葉面積指數(LAI)的相關性分析

本研究基于無人機冠層高光譜(450~950 nm),根據數學表達式(1)~式(3)在MATLAB中構建3種變換光譜基于任意兩波段的窄波段光譜指數與 LAI的相關性等勢圖(圖2)。結果表明窄波段光譜指數與LAI的相關性顯著提高,各變換光譜與 LAI相關性較高的波段組合范圍比較集中,主要分布在 680~750 nm 與 818~922 nm之間,9種窄波段光譜指數與LAI的相關系數最大值均在 0.65~0.85之間。根據相關系數最大的原則,篩選出估算冬小麥葉面積指數的 9種光譜指數的最佳波段組合如表2所示。通過查閱相關系數顯著性表,所選9種窄波段光譜指數均通過了0.01水平的顯著性檢驗。

2.2 冠層光譜與葉面積指數(LAI)的相關性分析

冬小麥冠層光譜以及變換光譜與LAI值的相關性關系如圖3所示。原始冠層光譜與LAI值在450~734 nm波段范圍內呈負相關關系,即LAI值隨著冠層光譜反射率的增加而減少;在734~950 nm波段范圍內,冬小麥LAI與反射率存在正相關關系,LAI值隨著冠層光譜反射率的增加而增加,其中在750~950 nm波段范圍內達到顯著正相關(P<0.01),并在918 nm處,原始光譜反射率與LAI值的相關系數達到最大值0.66。冠層一階導數光譜與 LAI的相關性較為復雜,波動性較大,在 730~770 nm波段范圍與LAI達到0.01水平上的顯著正相關,且相關性高于原始冠層光譜,對應最敏感位置在754 nm處,相關系數達到0.70。連續統去除光譜在450~778 nm范圍內與LAI值呈負相關關系,在778~950 nm波段范圍內呈正相關關系,其中482~762 nm波段范圍內,兩者呈現顯著性負相關關系,并在738 nm處相關系數到達最大值為-0.81。整體上來看,LAI與各變換光譜在可見光部分呈現負相關關系,在近紅外部分呈現正相關關系,且變換光譜與LAI的相關性略高于原始冠層光譜。

圖2 窄波段光譜指數與冬小麥葉面積指數的相關系數等勢圖Fig.2 Equipotential diagram of correlation coefficient between narrow band spectral index and winter wheat Leaf Area Index (LAI)

表2 窄波段光譜指數的最佳波段組合Table 2 Best band combination of narrow band spectral index

圖3 冬小麥葉面積指數與不同變換光譜反射率的相關性Fig.3 Correlation between winter wheat Leaf Area Index (LAI)and reflectance of different transformed spectra

2.3 基于窄波段光譜指數的葉面積指數(LAI)估算

以冬小麥葉面積指數為因變量,分別以原始光譜、一階導數光譜、連續統去除光譜下的差值光譜指數、比值光譜指數、歸一化光譜指數的最佳波段組合作為自變量,建立冬小麥葉面積指數的線性或非線性回歸模型,篩選出每種光譜變換下的最優單變量預測模型如表3所示。結果表明最優窄波段光譜指數與LAI的關系均表現為非線性,更適合用一元二次函數進行擬合。在3種光譜變換下,原始光譜和連續統去除光譜均以歸一化光譜指數(NDSI)為自變量構建的LAI預測模型最優,而導數變換光譜以比值光譜指數(RSI)為自變量構建的預測模型最優,3種最優單變量模型均到達了 0.01的顯著性檢驗,其R2均超過了0.55,其中基于原始光譜所構建的預測模型在幾種單變量預測模型中表現相對最佳,R2為0.75,RMSE為0.29,與基于連續統去除光譜所建模型精度較為接近。

表3 基于最優光譜指數的冬小麥葉面積指數含量預測模型分析(建模集)Table 3 Analysis results of winter wheat Leaf Area Index(LAI)content prediction model based on the best spectral index(Modeling set)

2.4 基于多光譜指數的葉面積指數(LAI)估算

眾多農作物生理生化參數的遙感估算研究表明,通過多光譜指數構建的農學參數估算模型精度顯著高于單光譜指數構建的模型,能有效提高提高農學參數的預測精度[30-31]。本研究嘗試基于窄波段多光譜指數構建以下模型,1)對各變換光譜分別進行基于最優差值光譜指數、比值光譜指數和歸一化光譜指數的多元線性回歸模型構建;2)以各變換光譜下的最優差值光譜指數、比值光譜指數和歸一化光譜指數合計9個光譜指數為自變量,構建基于BP神經網絡的LAI估算模型;3)構建基于以上9個光譜指數的隨機森林回歸估算模型。各模型的預測方程和精度如表4所示。

在 3種光譜變換下,采用多元線性回歸建模對冬小麥LAI進行預測,其中基于原始光譜建模所得決定系數最高(Rc2= 0.75),均方根誤差最低(RMSEc= 0.29);相對于多元線性回歸建模,2種機器學習算法模型精度均顯著提高。其中,對9個窄波段光譜指數進行LAI的隨機森林回歸建模優于BP神經網絡模型,其決定系數最高(Rc2= 0.92),均方根誤差最低(RMSEc= 0.17),且預測過程沒有出現過擬合和欠擬合現象。

2.5 模型精度對比

由表3和表4所構建的LAI預測模型可知,多元線性回歸模型建模精度相對傳統的一元回歸模型并無顯著提高,而基于多光譜指數的機器學習算法建模效果明顯優于傳統回歸模型。基于以上構建的冬小麥LAI估算模型,計算得到驗證集LAI的預測值,并將LAI預測值與實測值進行線性擬合。依據模型精度的評價指標對不同模型的預測精度進行比較分析,結果如表5所示。預測精度較建模精度整體上有所下降,基于連續統去除光譜的最佳窄波段光譜指數LAI估算模型(CRS-NDSI)精度略高于OS-NDSI和FDS-RSI,多元線性回歸模型驗證精度依然是連續統去除光譜(CRS-LAI-MLR)的預測精度略高于OS-LAI-MLR和FDS-LAI-MLR,但多元線性回歸模型驗證精度相對一元回歸模型并無顯著提高。基于FDS-RSI所構建的最優光譜指數驗證模型不具備預測樣本的能力(RPD < 1.5),其他7種驗證模型中除RFR模型外,均具有粗略的估測樣本能力(1.5 < RPD < 2),且RPD 由小到大依次為 FDS-LAI-MLR、OS-NDSI、OS-LAI-MLR、CRS-LAI-MLR、CRS-NDSI和LAI-BP,隨機森林回歸模型(LAI-RFR)RPD>2,驗證模型的決定系數為 0.77,均方根誤差為 0.27,在所有驗證模型中綜合表現最好。相對預測偏差RPD>1.5的所有模型的LAI實測值與預測值的空間分布如圖4所示,可以看出基于隨機森林構建的冬小麥LAI驗證模型,實測值與預測值的擬合分布更接近于 1∶1,說明估算模型具有極好的預測能力。

表4 基于多光譜指數的冬小麥葉面積指數含量預測模型(建模集)Table 4 Prediction model of winter wheat Leaf Area Index(LAI)content based on multispectral index (Modeling set)

表5 基于光譜指數的冬小麥葉面積指數含量驗證模型分析結果(驗證集)Table 5 Analysis results of winter wheat Leaf Area Index(LAI)content verification model based on spectral indices (Validation set)

圖4 冬小麥葉面積指數實測值和預測值分布Fig.4 Distribution of measured and predicted values of winter wheat Leaf Area Index(LAI)

3 討 論

葉面積指數是一種描述作物冠層結構的重要參數,快速準確地獲取LAI能精確地掌握農作物生長狀況[32-33]。冬小麥拔節期氮素積累量幾乎占總的氮積累量的一半,拔節期的生長狀況是后期產量和品質形成的關鍵,對拔節期冬小麥的葉面積指數進行精確監測可以為后期冬小麥的田間管理提供科學的理論依據[34-35]。基于高光譜相機(UHD 185)的作物理化參數反演得到了廣泛的驗證,田明璐等[10]對無人機平臺搭載 UHD高光譜相機所測的地物反射光譜曲線與地面對應目標的反射光譜進行比對,認為UHD高光譜影像具有較好的穩定性和光譜一致性。本研究基于無人機平臺和UHD高光譜影像進行拔節期LAI的估算研究,結果表明冬小麥UHD原始冠層高光譜、一階導數光譜和連續統去除光譜與LAI分別在750~950、730~770和 482~762 nm范圍內達到極顯著相關(P<0.01),通過分析原始冠層光譜、一階導數光譜和連續統去除光譜所構建的任意兩波段差值光譜指數、比值光譜指數和歸一化光譜指數與 LAI的相關關系,發現與 LAI相關性較高的波段組合中均包含了 680~750 nm的波段(表2),這與前人研究結果基本一致[14,36-39],本研究中基于連續統去除光譜的最佳波段組合為NDSI(738,822)。劉偉東等[40]研究發現水稻的葉面積光譜敏感波段集中在740~760 nm;蘇偉等[41]發現基于哨兵2號遙感影像構建有紅邊參與的 NDSI(705,783)光譜指數能夠較為準確地估算玉米LAI值,這都反映了紅邊位置對于LAI估算的重要性。因此,可以預見當無人機多光譜或者衛星多光譜在攜帶紅邊波段的情況下,將會提升LAI的估算能力。

本研究中基于一階導數變換光譜和連續統去除變換光譜與LAI的相關性較原始光譜顯著提高(圖2),與王偉東等[13]、姚霞等[15]、李粉玲等[16]、武旭梅等[42]和張雪紅等[43]研究結果基本一致,表明對原始光譜做變換并應用到作物生理參數反演方面是可行的。基于變換光譜構建的最優窄波段光譜指數與LAI的相關性顯著高于敏感波段,并且最優窄波段光譜指數與LAI的關系均表現為非線性,更適合用一元二次函數進行擬合,其中基于連續統去除光譜的 3類光譜指數所篩選的最佳波段組合均為 738和822 nm,且 NDSI(738,822)表現最優。基于多光譜指數構建的多元線性回歸模型相較于一元回歸模型精度并無明顯提高。但多元回歸模型和一元線性回歸模型均表明基于連續統變換光譜的模型精度較原始光譜和一階導數變換光譜下的模型精度有所提高,這與張雪紅等[43]、張金恒[44]的研究結果一致,表明以連續統去除法處理光譜反射率并用來建模是可行的,基于連續統去除光譜的 NDSI(738,822)指數具備較好的LAI的估算能力。以9個最優窄波段光譜指數為自變量所構建的BP神經網絡模型和隨機森林回歸模型相對于傳統的回歸模型精度明顯提高,其中隨機森林回歸模型精度最高,這是由于隨機森林算法建模能較好的容忍一些噪聲和異常值,只要調參精確,就不易出現過擬合,更適合解決一些非線性問題。基于該模型,本研究區拔節期冬小麥的LAI空間分布如圖5所示,基本與實際情況相符。但該模型對不同小麥品種、不同地理位置和生長環境下的LAI估算的適用性還有待深入研究。

圖5 基于LAI-RFR的研究區冬小麥葉面積指數反演圖Fig.5 Inversion map of winter wheat Leaf Area Index (LAI) in the study area based on LAI-RFR

4 結 論

無人機被廣泛應用在農業遙感方面,可以獲得作物各個生育期的信息,實現低空遙感作業與地面數據測量同步進行,為區域尺度高光譜遙感理論提供依據。本研究采用無人機低空遙感平臺獲取高光譜遙感影像,構建基于窄波段光譜指數的冬小麥葉面積指數(Leaf Area Index,LAI)估算模型,結果表明基于不同變換光譜的多個窄波段光譜指數構建隨機森林算法能夠獲得最優的LAI估算精度,其相對預測偏差為2.01,驗證集決定系數和均方根誤差分別為0.77和0.27。該模型可以作為關中地區拔節期冬小麥葉面積指數無人機高光譜遙感估算的基本模型,進而實現小區域尺度冬小麥LAI的快速準確獲取,為后期作物長勢以及估產等提供理論依據。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 黄色污网站在线观看| 黄片一区二区三区| 午夜无码一区二区三区| 亚洲成人播放| 久久99热这里只有精品免费看| 人妻无码中文字幕一区二区三区| 亚洲成人在线免费| 99视频国产精品| 国产一区在线视频观看| 国产美女在线免费观看| 日韩高清欧美| 东京热高清无码精品| 一级毛片网| 亚洲另类色| 日a本亚洲中文在线观看| 国产成人精品午夜视频'| 亚洲欧美另类日本| 国产手机在线观看| 麻豆国产在线观看一区二区| 国产福利小视频高清在线观看| 精品久久国产综合精麻豆| 伊在人亚洲香蕉精品播放| 国产精品久久久久婷婷五月| 国产精品成人啪精品视频| 大陆国产精品视频| 国产精品污污在线观看网站| 无码精品福利一区二区三区| 国产丝袜啪啪| 亚洲精品不卡午夜精品| 国产高潮流白浆视频| 免费无码在线观看| 国产亚洲视频播放9000| 亚洲综合九九| 欧美日韩中文国产va另类| 亚洲AV无码一二区三区在线播放| 久久精品无码一区二区国产区| 成人国产小视频| 美女免费精品高清毛片在线视| 91福利国产成人精品导航| 国内精品小视频在线| 超碰免费91| 久久久久中文字幕精品视频| 国产黄网永久免费| 一级一毛片a级毛片| 麻豆国产精品一二三在线观看| 久久国产热| a天堂视频| 久久无码免费束人妻| 九九久久精品免费观看| 波多野结衣中文字幕一区二区| 久久精品91麻豆| 国产va在线观看| 老司国产精品视频91| 色综合五月婷婷| 国产视频一二三区| 国产一区免费在线观看| 精品综合久久久久久97| 国产一区二区三区免费观看| 亚洲人人视频| 国模粉嫩小泬视频在线观看| 毛片视频网| 亚洲一区二区黄色| 国产一级小视频| 黄色网站不卡无码| 夜夜爽免费视频| 国产色伊人| 亚洲区一区| 欧美成人午夜视频| 成人字幕网视频在线观看| 国产精品对白刺激| 先锋资源久久| 亚洲人成网站观看在线观看| 亚洲成人精品| 国产门事件在线| 国产99视频免费精品是看6| 成人一级黄色毛片| 亚洲午夜福利精品无码不卡| 五月婷婷综合在线视频| 亚洲午夜福利精品无码不卡| 色天天综合久久久久综合片| 国产精品妖精视频| 免费人成视网站在线不卡|