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

融合GA與SVR算法的小麥條銹病特征優選與模型構建

2020-11-24 13:17:56白宗璠黃文江
農業機械學報 2020年11期
關鍵詞:特征模型

競 霞 張 騰 白宗璠 黃文江

(1.西安科技大學測繪科學與技術學院, 西安 710054; 2.中國科學院空天信息創新研究院, 北京 100094)

0 引言

小麥條銹病是小麥的主要病害之一,具有爆發性和流行性的特點,嚴重影響小麥的安全生產[1]。實時準確獲取小麥條銹病發病信息,對小麥病害的精準防治、減少農藥污染、提高小麥產量和質量具有重要意義。人工田間取樣的傳統病害調查方式費時、費力、時效性差,難以滿足多點同時大范圍監測的需求[2]。遙感技術具有快速、大面積、無損害等優點,已被廣泛應用于小麥病害監測中[3-5]。高光譜遙感數據能提供精細的地物光譜信息,定量識別作物受病蟲害的脅迫程度,受到了眾多學者的關注[6-9]。目前,小麥病害高光譜遙感監測主要利用光譜特征參量與病情指數之間的相關性作為優選敏感因子的標準,采用過濾法確定模型構建的特征參量[10-12],該方法從數據特征的結構出發,特征參量的選擇獨立于模型算法,算法簡單、通用性好[13],但忽略了各特征參量間的相關性,難以保證所選特征參量是適合于建模算法的最優因子,影響了模型構建的精度。封裝法將特征選擇算法和建模方法相結合,進行特征參量提取,能夠降低模型的復雜性,提高模型預測精度[14-16]。已有研究表明采用封裝法進行特征選擇和模型構建能夠有效剔除光譜冗余信息[17],以較少的敏感波段達到較高的反演精度[18],所建模型具有較強的預測能力[19]。在眾多的特征選擇和模型構建算法中,GA能夠模擬自然選擇和遺傳學機理的生物進化過程,具有原理和操作簡單、通用性強、能在短時間內達到全局最優的特點[20],在特征選擇中得到了廣泛應用[21-23]。而SVR算法具有適應性好、泛化能力強等特點,更適合小樣本數據[24]。

小麥受條銹病菌侵染后,其水分及葉綠素含量、光合速率和光能轉換率等生理生化指標均會發生變化[25]。反射率光譜數據能夠反映植被生化組分信息,但對作物光合活性不敏感,難以揭示植被光合生理狀態[26],且受陰影等非綠色景觀成分的背景噪聲影響較大[27]。自然光照條件下的葉綠素熒光受背景噪聲的影響較小[28],能夠直接探測到植被的光合作用狀態[29],遙感探測作物脅迫的精度較高。冠層SIF和反射率光譜在作物病害的遙感探測中各有優勢和局限性,已有研究表明,在反射率光譜數據中融合冠層SIF能夠提高小麥條銹病病情的估測精度[30-32]。基于此,本文將反射率吸收參數與冠層SIF數據作為初始特征集合,以SVR模型交叉驗證精度作為評價標準,利用封裝法將GA和SVR算法相融合,對特征子集與SVR參數進行聯合優選,建立基于遺傳算法優選特征參量的小麥條銹病嚴重度估測的GA-SVR模型,并將其與CC分析法提取特征參量所建模型的精度進行對比,探討GA-SVR模型在小麥條銹病遙感監測中的適用性,確定遙感監測小麥條銹病的敏感因子及模型。

1 材料和方法

1.1 試驗概況

試驗區位于河北省廊坊市中國農業科學院試驗站(39°30′40″N,116°36′20″E),小麥試驗品種為對條銹病比較敏感的銘賢169號,2018年4月9日(小麥起身期)采用濃度為0.09 mg/mL的孢子溶液利用噴霧法對小麥進行條銹病接種。試驗區小麥平均種植密度為113株/m2,分為健康組(編號為 A、D)和染病組(編號為B、C),健康組與染病組之間設置5 m隔離帶,并對健康組采用打藥處理防止病害發生。每個試驗組的面積為220 m2,每個組分為8個樣方(A1~A8、B1~B8、C1~C8、 D1~D8),即健康組和染病組各16個樣方。

1.2 數據獲取

1.2.1冠層光譜測量

試驗分別于2018年5月18日、5月24日和5月30日3個時期測量小麥條銹病不同病情下的冠層光譜數據,測量儀器為美國ASD公司生產的ASD Field Spec 4型地物光譜儀,其光譜分辨率為3 nm,采樣波長范圍為300~2 500 nm。為了減弱觀測角度和太陽天頂角對測量結果的影響,每次光譜測量時間均為北京時間11:00—12:30,且在數據測試前通過標準BaSO4參考板對冠層輻亮度數據進行校正。測量時探頭高度距離小麥冠層1.3 m,每個區域重復測量10次,取其平均值作為本次小麥冠層輻亮度。然后計算反射率,計算式為

(1)

式中ρ——冠層反射率

Ltarget——目標輻亮度,μW/(cm2·nm·sr)

Lboard——參考板輻亮度,μW/(cm2·nm·sr)

ρboard——參考板反射率

1.2.2病情指數調查

冠層光譜測量的同時采用5點取樣法調查小麥條銹病病情嚴重度,在每個樣方內選取對稱的5點,每點調查30株小麥的單葉病情嚴重度,并將其分為9個級別,即:0、1%、10%、20%、30%、45%、60%、80%和100%,通過記錄各嚴重度小麥葉片數計算病情指數[33],計算式為

(2)

式中DI——病情指數

t——各梯度級別

m——最高梯度

f——各梯度葉片數

1.2.3冠層SIF提取

冠層SIF提取方法主要包括基于輻亮度數據的直接提取算法和基于反射率數據的間接提取算法兩種。基于輻亮度的冠層SIF直接提取算法是利用夫瑯和費暗線原理,通過FLD(Fraunhofer line discrimination)、3FLD(3 bands FLD)以及iFLD(Improved FLD)等算法[34]估測SIF強度。已有研究表明,FLD算法在一定程度上會高估熒光值,iFLD和3FLD方法則能夠得到較為可靠的熒光估測結果[35],而3FLD算法較iFLD算法更穩健[36-37]。因此,本文利用3FLD方法估算O2-A波段(760 nm)冠層SIF的絕對強度,計算式為

(3)

其中

ωleft=(λright-λin)/(λright-λleft)

(4)

ωright=(λin-λleft)/(λright-λleft)

(5)

λin、λleft、λright——吸收線內、左、右波段波長

ωleft、ωright——吸收線左、右兩參考波段權重

Iin、Ileft、Iright——吸收線內、左、右波段的太陽輻照度,μW/(cm2·nm)

Lin、Lleft、Lright——吸收線內、左、右波段的植被冠層反射輻亮度,μW/(cm2·nm·sr)

除了利用夫瑯和費吸收線直接提取冠層SIF外,基于反射率比值的方法也可獲取650~800 nm波段的熒光信息,是提取冠層SIF強度的間接方法[38]。基于反射率比值提取冠層SIF主要是利用熒光對650~800 nm波段范圍內反射率的影響,獲取與葉綠素熒光相關的光學指數,包括反射率比值熒光指數、反射率導數熒光指數和填充指數3種[39]。反射率比值熒光指數通過對680 nm或740 nm附近受熒光影響強的波段和一個受熒光影響較少或不受影響波段的比值運算,消除反射率信息,獲取熒光強度。反射率導數熒光指數利用反射率一階導數的比值反映熒光信息,能夠增強熒光對反射率的影響,突出由熒光引起的細微變化。填充指數通過兩個波段反射率的差間接反映熒光信息,但該指數受夫瑯和費暗線的深度和熒光大小的共同影響,由于大氣和太陽觀測幾何的變化都會影響夫瑯和費暗線的深度,所以該指數僅適用于在相同時間和觀測條件下的數據對比[40]。

1.2.4吸收特征參數計算

吸收特征是植物組織結構、色素含量、水分和蛋白質中各種基團對反射光譜響應的重要特征[33]。反射率光譜中的細微吸收特征能夠在連續統去除的歸一化過程中得到增強[41]。連續統是指將光譜曲線上凸出的“峰”值點用直線連接,并使所形成的折線在峰值點上的外角大于180°[42],連續去除后的光譜反射率計算式為

(6)

式中ρc——連續統線

原始光譜經過連續統去除后,波長起點與終點相對反射率等于1,處于起點與終點間的波長,其反射率都小于1。本文在小麥冠層反射率光譜連續統去除基礎上參照文獻[43-45]提取了吸收峰深度(H)、吸收波長位置(P)、吸收峰面積(A)、吸收峰左面積(A1)、吸收峰右面積(A2)、對稱度(S)和面積歸一化吸收深度(NMAD)等特征參量。

吸收峰深度H計算式為

H=1-ρ′min

(7)

式中ρ′min——連續統去除后的最小反射率

吸收波長位置P計算式為

P=λH

(8)

式中λH——吸收峰深度對應的波長

吸收峰面積A計算式為

(9)

式中a——吸收帶起始波長

b——吸收帶終止波長

di——吸收深度

Δλ——波長增量

吸收峰左面積A1計算式為

(10)

吸收峰右面積A2計算式為

(11)

對稱度S計算式為

S=A1/A

(12)

面積歸一化吸收深度NMAD計算式為

NMAD=H/A

(13)

1.3 GA-SVR算法

1.3.1SVR算法

SVR算法的基本思想是利用訓練樣本建立一個回歸超平面,將樣本逼近超平面以使樣本點到該平面的總偏差達到最小[46],其計算式為

(14)

β——偏置φ——映射函數

為避免高維內積運算,引入核函數k(xi,x)代替(φ(xi),φ(x)),則式(14)轉換為

(15)

SVR算法常用的核函數有線性核函數、徑向基(Radial basic function,RBF)核函數、多項式核函數以及Sigmoid核函數等[47]。相較于線性核函數僅能處理線性問題,RBF核函數可處理自變量與因變量之間復雜的非線性問題。此外,應用RBF核函數的效果通常優于多項式核函數與Sigmoid核函數,并且需要設置的參數更少[48]。因此,本文在基于GA-SVR算法構建小麥條銹病嚴重度估測模型時采用了RBF核,其計算式為

(16)

式中g——核參數

1.3.2GA

GA是一種模擬自然選擇和自然遺傳的隨機搜索與最優化求解方法,它將具體的求最優解問題進行編碼轉換為特定數量并具有基因序列的種群,其中種群中的每一個體代表每一個解,以適應度函數評價種群中個體的求解結果,并根據個體的適應度,進行選擇、交叉和變異3種遺傳操作,以保證在進化過程中種群的適應性不斷增強,進而實現個體的優勝劣汰,獲取問題的最優解[49](圖1)。

圖1 遺傳算法流程圖Fig.1 Flow chart of genetic algorithm

1.3.3基于GA-SVR算法的特征優選

利用GA-SVR算法優選模型特征參量是基于二進制編碼方法對初始群體中的個體染色體進行編碼(圖2),圖中每個染色體由3部分構成,第1部分C1-Cx與第2部分g1-gy分別為懲罰參數C與核參數g的二進制編碼,其中下標x和y表示懲罰參數C與核參數g的二進制位數,決定解碼為十進制數的精度。第3部分f1-fz為所有特征集合的編碼,下標z為特征總個數。當編碼中的二進制數為“1”時表示該特征被選中,為“0”時表示該特征未被選中。

圖2 編碼示意圖Fig.2 Coding diagram

基于GA-SVR算法進行特征優選時,需要對種群數量、最大迭代次數、交叉概率和遺傳概率等參數進行設置。為了保證遺傳算法具有全局搜索能力,在上述遺傳參數設置的基礎上按 “單點交叉”方式對種群中的個體進行兩兩交叉,同時以基本位變異的方式依據較小的變異概率對交叉完成后的個體進行變異操作,以防止算法陷入局部最優。經過選擇-交叉-變異等遺傳操作形成子代種群,重復進行解碼、適應度評價、種群更新等步驟,在滿足迭代終止條件時算法停止。

1.4 模型構建與評價指標

為了保證評價結果的穩定性和可靠性,提高模型的泛化能力,減弱樣本數據對評價結果的影響,本文將52個原始樣本重復進行3次隨機分組(記為a、b、c 3組),每組按照3∶1的比例隨機劃分為訓練集和驗證集,其中39個訓練樣本用于特征選擇與模型構建,13個驗證樣本用于模型評價。

選用預測值和實測值之間的決定系數R2與均方根誤差(RMSE)2個指標進行模型精度評價,其中決定系數R2越大,均方根誤差(RMSE)越小,預測樣本與實測樣本偏差越小,模型精度越高。

2 結果與分析

2.1 特征參量優選

特征參量很大程度上影響了模型的性能,為了提高小麥條銹病遙感監測精度,利用反射率光譜數據在作物生化參數探測方面的優勢和SIF在光合生理診斷方面的優勢,計算了吸收深度等7個反射率光譜吸收特征及15個冠層SIF數據,利用封裝法將GA和SVR算法相結合進行小麥條銹病遙感監測特征參量優選,并將其與CC特征優選方法進行對比,分析兩種方法提取特征參量所建模型精度,以確定遙感探測小麥條銹病的敏感因子及模型。

2.1.1CC特征優選

CC是一種過濾性的特征選擇方法,具有計算量小及所選特征參量通用性較好的優點。本文在計算冠層SIF及反射率吸收特征參數與病情嚴重度間相關系數的基礎上,挑選達到0.001極顯著水平的特征參量作為遙感監測小麥條銹病嚴重度的敏感因子。

(1)冠層SIF強度特征參量

分別利用輻亮度和反射率計算了小麥條銹病不同病情下的冠層SIF信息。利用輻亮度數據計算冠層SIF時,由于O2-A (760 nm)波段氧氣吸收形成的夫瑯和費暗線特征明顯[50],熒光估測精度高[51],因此本文利用輻亮度基于3FLD方法估測了760 nm處的SIF強度,并將計算得到的SIF絕對強度除以夫瑯和費吸收線內的太陽入射輻照度,獲取O2-A波段的SIF相對強度,以消除太陽光照強度等外界因素對葉綠素熒光提取結果的影響[52-53]。

(17)

式中F′——SIF相對強度

利用反射率間接估測冠層SIF時,由于填充指數法僅適用于在相同時間和觀測條件下的數據對比,因此本文僅計算了反射率比值熒光指數和反射率一階導數熒光指數。目前常用的反射率比值熒光指數主要有ZARCO-TEJADA等[54-55]提出的

表1 冠層SIF與病情指數相關系數(n=52)Tab.1 Correlation coefficient between canopy SIF and disease index (n=52)

(2)反射率吸收特征優選

在對400~550 nm和550~770 nm兩個吸收波段范圍內的反射率光譜數據進行連續統去除的基礎上提取了吸收位置等光譜吸收特征參數,并計算其與小麥條銹病嚴重度的相關性(表2)。由表2可以看出,在400~550 nm和550~770 nm 2個吸收波段范圍內共有吸收位置(P)、最大吸收深度(H)以及面積歸一化最大吸收深度(NMAD)等9個特征參量與小麥條銹病嚴重度達到了0.001水平的極顯著相關,可以作為遙感監測小麥條銹病嚴重度的模型輸入變量。

表2 吸收特征參數與病情指數相關系數(n=52)Tab.2 Correlation coefficient between absorption characteristics and disease index (n=52)

2.1.2GA-SVR特征優選

利用GA-SVR算法優選遙感監測小麥條銹病敏感因子時,GA及SVR參數設置對迭代尋優的速率與準確率有較大的影響,基于此本文通過對訓練樣本的多次仿真確定SVR算法中懲罰參數C與核參數g的尋優范圍及GA各參數設置(表3)。然后在上述參數設置的基礎上將GA-SVR特征選擇算法獨立運行20次,單次運行最大迭代次數設定為100,以訓練集10折交叉驗證均方誤差(Mean square error,MSE)作為評價個體優劣的適應度,其中適應度最小的個體所對應的特征參量以及懲罰參數C和核參數g即為GA-SVR模型構建時的最優特征集合以及SVR算法的最佳參數組合。表4為a組數據運行20次得到的不同適應度及其對應的最優特征集合。

表3 GA-SVR各參數設置Tab.3 GA-SVR parameter settings

表4 a組數據單次適應度最小值及其對應特征與參數組合Tab.4 Minimum single fitness of group a of data and their corresponding characteristics and parameter combinations

由表4可知,第8次運行時的適應度為GA-SVR算法20次運行中的最佳適應度,對應的懲罰參數C(94.937 9)與核參數g(0.967 0)為SVR算法的最優參數組合,特征集合{4,5,17,22,24,25,27}為遙感監測小麥條銹病嚴重度的最佳特征參量。同理得到b組和c組的最優特征集合及SVR算法的最優參數組合(表5)。

2.2 模型構建和精度評價

2.2.1模型構建

基于GA-SVR算法優選的最優特征組合建立遙感監測小麥條銹病嚴重度的GA-SVR模型(圖3),并將其與CC法算法優選的特征參量集合所建的CC-SVR模型(圖4)進行對比分析。

表5 3組數據GA-SVR優選特征組合及參數組合Tab.5 Three sets of data GA-SVR preferred feature combination and parameter combination

圖3 GA-SVR模型訓練集預測結果(n=39)Fig.3 GA-SVR model prediction result on training set(n=39)

圖4 CC-SVR模型訓練集預測結果(n=39)Fig.4 CC-SVR model prediction result on training set(n=39)

由圖3、4可以看出,3個樣本組中GA-SVR模型預測DI值和實測DI值之間的R2均在0.9以上,相對于同組CC-SVR模型分別提高了9.9%、10.0%、5.6%,RMSE分別降低了26.3%、26.5%、19.3%,GA-SVR模型對小麥條銹病嚴重度的估測精度優于CC-SVR模型。說明將GA-SVR模型優選的特征作為模型輸入能夠明顯地提高模型對訓練集病情指數的預測能力,而CC法優選的特征雖然對病情指數敏感,但從CC-SVR模型對訓練集病情指數的預測精度來看,其并非為SVR模型的最佳輸入特征。

2.2.2模型評價

為了在樣本數量有限的情況下更客觀地評價模型精度,盡可能減少“過擬合”問題對模型評價結果的影響,采用能夠充分利用樣本中所有數據的保留樣本交叉檢驗方式對CC分析法以及GA-SVR方法所選特征構建的小麥條銹病嚴重度預測模型進行評價(表6)。由表6可以看出,3組驗證樣本中GA-SVR模型的預測精度較同組CC-SVR模型均有不同程度的提高,其中GA-SVR模型估測DI值和實測DI值間的R2較同組CC-SVR模型至少提高了2.7%,平均提高了17.8%,RMSE較同組CC-SVR模型至少減少了10.1%,平均減少了32.1%。綜合分析圖3、圖4和表6可知,無論是訓練集樣本還是驗證集樣本,GA-SVR模型對DI的預測能力均優于CC-SVR模型,將特征優選和模型構建相融合的GA-SVR模型能夠提高小麥條銹病嚴重度的預測精度,且模型穩定性和泛化能力均優于CC-SVR模型,更適合于小麥條銹病嚴重度的敏感因子優選與模型構建。這是因為CC算法是以顯著檢驗水平作為特征優選的條件,僅考慮了特征參量與病情嚴重度之間的相關性,沒有考慮不同特征參量之間是否存在數據的冗余以及所選的特征參量是否最適合于模型構建算法,影響了模型的預測精度。而GA-SVR算法綜合利用了GA的全局搜索特點和SVR算法較優的學習能力,考慮了SVR參數對特征參量優選的影響,將所有特征與SVR參數進行聯合編碼,以10折交叉驗證最小MSE作為特征與參數的優選依據,對所有可能選擇到的特征進行有向搜索,保證了所選特征的敏感性,提升了模型的預測精度。

表6 小區試驗數據驗證結果(m=13)Tab.6 Verification results of plot experiment (m=13)

2.2.3大田調查數據驗證

2018年5月12日于陜西省漢中市寧強縣大田小麥種植區域調查條銹病的發病狀況,利用ASD Field Spec 4型地物光譜儀同步獲取42個對應樣本點的冠層光譜數據,并基于此數據對本文所建立的模型進行驗證(表7),以評價模型的可靠性與普適性。

表7 大田調查數據驗證結果(m=42)Tab.7 Verification results of field survey data (m=42)

由表7可知,3個樣本組建立的GA-SVR模型預測DI值與實測DI值間的R2均高于0.7,比CC-SVR模型至少提高了66.2%,RMSE均不高于0.150,比CC-SVR模型至少減少了31.5%,聯合GA與SVR算法的特征優選方法能夠提高小麥條銹病的遙感探測精度,GA-SVR模型具有較好的可靠性和普適性。

3 討論

融合GA和SVR算法優選小麥條銹病遙感監測的敏感因子時,3個樣本組中優選的特征參量的一致性不高,說明GA-SVR方法進行特征選擇時,受到原始樣本的影響較大。這是因為封裝法是以模型的精度作為特征優選的標準,樣本的差異影響交叉驗證的精度,小樣本訓練集的隨機劃分難以代表所有樣本間的差異,進而導致3個樣本組所選的特征參量不一致,能否利用大量的訓練樣本解決這一問題還有待繼續研究。

本文協同冠層SIF數據和反射率光譜進行小麥條銹病遙感監測時,反射率數據僅利用了對植被葉綠素含量敏感的可見光波段的吸收特征,而小麥在受到條銹病菌侵染時,水分、花青素等生化組分含量也會發生相應的變化[58]。在構建小麥條銹病嚴重度反演模型時,若增加能夠反映各種生化組分的光譜指數,分析不同光譜指數對條銹病害脅迫的響應特性,能否提高小麥條銹病的遙感監測精度還有待進一步研究。

綜合利用反射率光譜在作物生化參數探測方面的優勢及SIF在光合生理診斷方面的優勢構建了小麥條銹病嚴重度估測模型,并沒有考慮反射率吸收特征和SIF數據對模型的貢獻率問題,如何量化反射率光譜數據和SIF數據對遙感監測小麥條銹病嚴重度的貢獻率,確定各個輸入特征在模型中的權重,進一步提高小麥條銹病的遙感監測精度,則是一個值得研究的問題。

在利用SVR算法構建小麥條銹病嚴重度估測模型時,僅考慮了高斯核函數,并沒有分析所選特征參量與小麥條銹病嚴重度之間的映射關系,確定不同特征參量與病情嚴重度之間的最優映射函數進而提高小麥條銹病嚴重度的遙感估測精度則是下步要研究的問題。

本文在構建小麥條銹病嚴重度反演模型時,僅選用了對小樣本適用性較好的SVR算法,并未與其它機器學習算法進行對比分析。已有研究將GA與隨機森林算法結合進行敏感波段的選擇,并取得了較好的研究成果[18],但是將GA與其他機器學習算法結合進行小麥條銹病敏感特征選擇能否改善小麥條銹病嚴重度反演精度還需進一步研究。

4 結束語

基于小麥條銹病不同病情嚴重度下的冠層SIF數據與反射率吸收特征參量,分別利用GA-SVR算法與CC算法優選了小麥條銹病遙感監測的敏感因子,在此基礎上構建了小麥條銹病的遙感監測模型。結果表明,GA-SVR算法優選的小麥條銹病遙感監測的敏感因子為7~8個,而CC分析法則優選出了22個因子,GA-SVR算法極大地降低了特征空間的維度,減少了特征參量間的冗余,降低了模型的復雜性。從模型的預測性能來看,無論是小區試驗數據還是野外大田數據,GA-SVR模型的預測精度均高于CC-SVR模型。以小區試驗數據為驗證樣本時,GA-SVR模型預測DI值與實測DI值之間R2均在0.9以上,RMSE均小于0.1,比CC-SVR模型的R2至少提高了2.7%,RMSE至少降低了10.1%。以大田數據為驗證樣本時,GA-SVR模型預測DI值與實測DI值間的R2比CC-SVR模型至少提高了66.2%,RMSE至少降低了31.5%。

猜你喜歡
特征模型
一半模型
抓住特征巧觀察
重要模型『一線三等角』
新型冠狀病毒及其流行病學特征認識
重尾非線性自回歸模型自加權M-估計的漸近分布
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 免费 国产 无码久久久| 99久久精品免费看国产免费软件| 男女猛烈无遮挡午夜视频| 茄子视频毛片免费观看| 久久综合久久鬼| a毛片在线免费观看| 怡红院美国分院一区二区| 久久人人爽人人爽人人片aV东京热| jijzzizz老师出水喷水喷出| 色综合激情网| 国产一区自拍视频| www.精品国产| 亚洲香蕉在线| 久久综合伊人77777| 国产99精品久久| 久草视频中文| 99精品高清在线播放| 老色鬼久久亚洲AV综合| 亚洲无码高清视频在线观看| 一级毛片在线播放| 任我操在线视频| 国产精品yjizz视频网一二区| 在线观看av永久| 亚洲精品图区| 国产成人欧美| 亚洲欧美另类中文字幕| 91色国产在线| 精品国产aⅴ一区二区三区| 亚洲人成网7777777国产| 91在线丝袜| 深爱婷婷激情网| 亚洲人成网站日本片| 国产真实乱人视频| 欧美日一级片| 网友自拍视频精品区| 素人激情视频福利| 在线精品亚洲国产| 2022国产91精品久久久久久| 老司机久久精品视频| 国产h视频在线观看视频| 国产女人在线| 全部免费毛片免费播放| 在线播放真实国产乱子伦| 91麻豆精品国产高清在线| 亚洲AⅤ无码日韩AV无码网站| 午夜a视频| 97人人做人人爽香蕉精品| 成人国产精品一级毛片天堂| 四虎国产精品永久一区| 国产精品三级av及在线观看| 久久国产精品国产自线拍| 香蕉视频在线精品| 一区二区自拍| 国产美女精品人人做人人爽| 亚洲一级毛片免费看| 日韩专区第一页| 欧美视频在线第一页| 91在线中文| 亚洲浓毛av| 蜜臀av性久久久久蜜臀aⅴ麻豆| 国产视频入口| 久久国产高潮流白浆免费观看| 亚洲成肉网| 欧美a在线视频| 精品一区二区三区自慰喷水| 97青草最新免费精品视频| 亚洲日产2021三区在线| 国产成人高清精品免费软件 | 伊人网址在线| 国产成人综合网在线观看| 在线观看国产一区二区三区99| 美女一区二区在线观看| 丰满人妻久久中文字幕| 在线国产毛片| 在线观看国产小视频| 蜜臀AVWWW国产天堂| 三上悠亚在线精品二区| 午夜高清国产拍精品| 国产精品爆乳99久久| 成人福利在线视频免费观看| 99er这里只有精品| 欧美亚洲综合免费精品高清在线观看|