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

南方丘陵稻田土堿解氮高光譜特征及反演模型研究

2015-01-04 06:19:26葉英聰謝碧裕匡麗花
自然資源遙感 2015年2期

郭 熙,葉英聰,謝碧裕,匡麗花,謝 文

(1.江西省農業科學院土壤肥料與資源環境研究所,南昌 330200;2.中國農業科學院農業資源與農業區劃所,北京 100081;3.江西農業大學江西省鄱陽湖流域農業資源與生態重點實驗室,南昌 330045)

0 引言

土壤光譜反射特性是土壤理化特征和內在結構的光譜行為的綜合,其研究是土壤遙感的物理基礎,也為土壤屬性研究提供了一個新的途徑和指標[1]。基于土壤反射光譜特性的高光譜遙感技術的廣泛應用將為農業生產管理、土壤定量遙感研究、土壤水土流失防治、土壤資源保護、陸地生態系統相關研究等提供工作方法[2-3]。土壤光譜是各種土壤理化性質的綜合反映,可以通過光譜特征來估計土壤重金屬、有機質等物質的含量[4-5]。土壤反射光譜特性是土壤的基本性質之一,是土壤理化參數的綜合反映,土壤礦物、有機質、水分、土壤質地、粗糙度和秸稈殘留物等均影響土壤的光譜反射率[6]。

土壤堿解氮包括無機態氮(銨態氮、硝態氮)及易水解的有機態氮(氨基酸、酰銨和易分解的蛋白質),易被植物吸收,是反映土壤供氮能力的指標之一。作為一項主要的土壤氮素(N)分析項目,土壤堿解氮質量分數的高低,能大致反映近期土壤氮素的供應情況,與作物生長和產量有一定的相關性,可作為土壤有效氮的指標。實驗室內土壤堿解氮含量的測定有堿解擴散法、堿解蒸餾法和儀器法等3種。堿解擴散法中堿解擴散時間為24 h,測定時間較長[7];堿解蒸餾法雖蒸餾時間為8 min,但加熱蒸餾裝置落后,需人為經驗控制,操作繁瑣[8]。

大量研究表明,350~2 500 nm波段的高光譜反射數據能夠反映土壤理化參數的細微差異,可以用于土壤理化參數的反演。其中土壤水分、有機質、N素及鐵的氧化物與土壤反射率的關系得到了大量的研究,并建立了定量估測模型。劉煥軍分析了松嫩平原主要土壤的反射光譜特征及其與土壤理化參數的關系,提出了基于光譜反射特性的土壤有機質含量高光譜模型[6];何挺研究了土壤氧化鐵與土壤光譜反射特性的關系,構建了土壤氧化鐵指數,并應用逐步回歸分析法構建了反演模型[9];徐馳研究了內蒙古河套灌區土壤含鹽量、pH含量與土壤光譜反射特性的關系,并構建了反演模型[10]。

國內外學者利用遙感技術研究土壤性質,多在可見光—近紅外光譜條件下,應用多元線性、偏最小二乘回歸等方法對影響土壤光譜反射特性的主要因素(如有機質、水分、鹽漬度、總氮及重金屬)含量進行預測[11]。應用高光譜數據預測土壤堿解氮含量相對較少,陳紅艷以山東潮土為例,對土壤堿解氮含量構建了高光譜估測模型[12];徐麗華研究了三峽庫區王家溝小流域,選擇了特征光譜和特征波段,建立了紫色土和水稻土的堿解氮養分預測模型[13]。總體來說,基于南方丘陵稻田土高光譜特性與土壤堿解氮含量關系的研究較少。

本研究以江西省興國縣稻田土高光譜反射率為研究對象,分析南方丘陵稻田土堿解氮的光譜響應波段,提取定量光譜指標,建立基于反射光譜特征的南方丘陵稻田土高光譜反演模型,為南方丘陵稻田土堿解氮的快速測定提供支持。

1 材料與方法

1.1 研究區概況

以江西省興國縣為研究區,包括永豐鄉、高行鎮、社富鄉、隆坪鄉、鼎龍鄉和長崗鄉6個鄉鎮。興國縣地處我國中亞熱帶南部,E115°01'~115°51',N26°03'~26°41'之間,主要地貌類型有河谷沖積平原、崗地、丘陵山地等,母巖主要有第四紀紅色粘土、砂頁巖、花崗巖、千枚巖等,主要的土壤類型為紅壤、黃壤、紫色土和水稻土。研究區域水稻土又分為表潛青潮沙泥田、漚水紫泥田、淹育沙泥田和潴育潮沙泥田、潴育麻沙泥田共5個土屬。興國縣屬雙季稻區,是江西省糧食主產縣之一。常年水稻種植面積在49 667 hm2以上,總產量約 26.5 ×104t,年人均占有量約為340 kg。

1.2 研究方法

1.2.1土壤樣品采集

2013年4月,在興國縣境內采集0~20 cm耕層土樣,采樣同時利用GPS記錄樣點的經緯度。根據鄉鎮到村依次排列編號:1—2號在永豐鄉豪溪村,3—8號在永豐鄉馬良村,9—14號在永豐鄉馬良村,15—20號在隆坪鄉蘭溪村,21—26號在鼎龍鄉楊村,27—32號在長崗鄉榔木村,33—38號在社富鄉東韶村,39—43號在高行鎮高興村,總計43個樣點。1—2號樣點土屬為淹育沙泥土,3—14號樣點土屬為漚水紫泥田,15—20,33—38號樣點土屬為潴育麻沙泥田,21—26號樣點土屬為表潛青潮沙泥田,27—32及39—43號樣點土屬為潴育潮沙泥田。

1.2.2 土壤高光譜測定

首先在室內對43個土樣進行研磨和風干,過1mm篩。光譜測試采用美國Spectra Vista公司生產的SVC HR-768地物波譜儀,測定43個稻田土樣的反射光譜,其波長范圍為350~2 500 nm,采樣間隔為1.5 nm(350 ~1 000 nm 區間)、7.5 nm(1 000 ~1 850 nm區間)和5 nm(1 850~2 500nm區間),輸出波段數為768個。光源是功率為1 000 W的鹵素燈,距土壤樣品表面100 cm,天頂角為30°,采用25°視場角探頭,探頭位于土壤樣本垂直上方15 cm處。測試之前先除去輻射強度中暗電流的影響,以白板進行定標,每個土樣掃描時間為5 s。為保證試驗數據的準確度,每個土樣采集3條光譜曲線,經算術平均后得到該土樣的實際光譜反射數據。

1.2.3 土壤堿解氮含量測定

堿解氮含量采用堿解擴散法測定[7]。精確稱土樣2 g均勻鋪在擴散皿外室中,在擴散皿內室中加入 2 mL H3BO3(ρ(H3BO3)=20 g·L-1)指示劑溶液;然后在擴散皿的外室邊緣涂上堿性膠液,迅速加入 10 mL NaOH 溶液(c(NaOH)=1.0 mol·L-1)于皿的外室中;隨后放入40℃恒溫箱中,24 h后取出;取出擴散皿,用 HCl(c(HCl)=0.01 mol·L-1)標準溶液,半微量滴定管滴定內室硼酸中所吸收的NH3,由藍色到微紅色為終點,記下HCl用量;最后,通過公式計算土樣中堿解氮的含量。堿解氮含量描述性統計量見表1。

表1 43個土樣堿解氮含量描述性統計表Tab.1 Descriptive statistics of paddy soil available nitrogen content

1.2.4 光譜數據處理與分析

光譜數據處理方法主要涉及重采樣、連續統去除、光譜微分處理及光譜數據變換方法。

1)重采樣處理方法。利用SVC HR-768地物波譜儀自帶的光譜數據處理軟件對原始光譜數據進行重采樣處理,重采樣間隔為1 nm,再導出重采樣數據。SVC HR-768地物波譜儀測定土壤樣點的反射光譜數據在1 000 nm,1 900 nm附近的接縫處以及前后邊緣波段344~399 nm,2 451~2 494 nm處噪聲較大,使得重采樣所得的原始光譜曲線的相鄰波段之間存在信息重合,導致整個光譜數據信息冗余,而其他波段的信噪比高,約為1 000∶1。因此,對光譜數據以10 nm為間隔進行算術平均運算,處理后的光譜曲線更加平滑,同時仍然維持了原光譜的主要特征,故對導出的每個土樣的光譜曲線去除前后噪聲較大的344~380 nm和2 451~2 494 nm邊緣波段。

2)連續統去除法。連續統去除法也叫去包絡線法,是一種典型的光譜分析法,定義為逐點直線連接隨波長變化的吸收或反射凸出的“峰”值點,并使折線在“峰”值點上的外角大于180°。它可以有效地突出光譜曲線吸收和反射特征,并將其歸一到一致的光譜背景上,有利于和其他光譜曲線進行特征數值比較,從而提取出特征波段進行分類識別,去掉包絡線后變為光譜波段深度曲線。光譜波段深度曲線計算公式為

式中:Rs,R,Rc分別為光譜波段深度、原始光譜和光譜包絡線;λ為波長。

3)光譜一階微分方法。光譜微分方法(導數算法)是常用的光譜增強方法,該方法對光譜信噪比非常敏感。研究表明,光譜的低階微分處理對噪聲影響敏感性較低,對不同的背景和噪聲有去除作用[14]。微分光譜可以消除基線漂移或平緩背景干擾的影響,并可以提供比原始光譜更高的分辨率和更清晰的光譜輪廓變換。其基本處理方法是:先確定導數窗口寬度Δλ,根據導數的定義計算波長λ的導數,再逐步移動依次計算所有波長的導數,由此得出導數光譜。土壤高光譜在波長λ處的一階微分光譜計算公式為

式中:λ 為波長,λ =351,352,…,2 499 nm;R(λ)為波長λ的光譜反射系數;R'(λ)為波長λ的光譜反射系數的一階微分。在實際計算中,一般用光譜的差分作為微分的有限近似。土壤原始反射率光譜經過微分變換后,曲線隨波長變化更明顯,更能凸顯細微信息差異引起的變化,也更易找出曲線拐點位置。

1.2.5 模型建立與驗證

土壤堿解氮反演模型的預測精度采用預測值和實測值的決定系數R2,均方根誤差RMSE和相對分析誤差RPD(檢驗樣本標準SD與預測均方根誤差RMSE的比值)來衡量。R2越大,RMSE越小,說明模型的精度越高。另外,當RPD>2時,表明模型具有極好的預測能力;當1.4≤RPD<2時,表明模型可對樣品作粗略估測;而RPD<1.4時則表明模型無法對樣品進行預測[15]。

2 結果與分析

2.1 南方丘陵稻田土室內反射光譜特征

對稻田土光譜曲線進行去包絡線處理,能有效突出光譜曲線吸收和反射特征,并將其歸一到一致的光譜背景上,有利于和其他光譜曲線進行特征數值比較,從而提取出特征波段進行分類識別。

本研究選取8號、23號和40號樣本進行研究分析,3個樣本土屬分別為漚水紫泥田、潴育麻沙泥田和潴育潮沙泥田,堿解氮含量分別為152 mg/kg,154 mg/kg和155 mg/kg。對比3個樣本光譜重采樣反射率來探討同一堿解氮含量下土屬與樣本光譜反射率的關系(由于淹育沙泥土和表潛青潮泥田2個土屬樣本數較少,因而沒有在這次研究中探討)。同一堿解氮含量不同土屬光譜重采樣曲線見圖1。

圖1 同一堿解氮含量不同土屬光譜重采樣曲線Fig.1 Spectrum resample curves of the same available nitrogen content,different soil types

從圖1可以看出,同一堿解氮含量水平3個不同土屬樣本光譜曲線反射率差異不大。結合周學文等[16]研究南方丘陵地區水田土壤所采用的土壤養分豐缺度與分級標準,根據堿解氮含量將土樣分為“極缺”、“缺”、“中”、“豐富”和“偏高”5個等級,相應的堿解氮含量區間分別為<50 mg/kg,[50,100)mg/kg,[100,150)mg/kg,[150,200)mg/kg 和≥200 mg/kg。將43個土樣根據堿解氮含量進行分區,對屬于同一區間的樣點光譜重采樣數據取均值,然后進行去包絡線處理。

圖2為不同堿解氮區間含量土壤反射率曲線。

圖2 不同堿解氮含量區間光譜重采樣曲線Fig.2 Spectrum resample curves of different available nitrogen content range

由圖2可知,受水分吸收的影響,1 350~1 420 nm和1 844~1 944 nm這2個波段出現水分吸收谷,2 444~2 494 nm波譜范圍內反射率數據噪聲較大。總體來說,稻田土光譜反射率曲線較平滑;在344~1 294 nm波譜范圍內,光譜曲線呈現單調上升的趨勢,曲線上凸;1 350~1 420 nm波譜范圍內,出現水分吸收谷;1 420~1 800 nm波譜范圍內,光譜曲線呈現單調上升趨勢,出現明顯的波峰;1 800~1 920 nm波譜范圍內,出現水分吸收谷,且受噪聲影響較大;1 920~2 360 nm波譜范圍內,在2 194 nm波段出現明顯的光譜吸收谷;2 360~2 494 nm波段內,光譜曲線呈現單調下降趨勢。總體而言,對比堿解氮含量不同的4條曲線,堿解氮含量越高,光譜反射率越低。

圖3為不同堿解氮含量區間土壤去包絡線。

圖3 不同堿解氮含量區間光譜去包絡線Fig.3 Continuum removed curve of different available nitrogen content range

比較圖3可以發現,堿解氮含量不同,稻田土去包絡線曲線在≥700 nm的光譜范圍內差異較小;而<700 nm的光譜范圍內稻田土去包絡線曲線呈現顯著差異,在394~444 nm波段,堿解氮含量區間為[50,100)mg/kg 和[100,150)mg/kg,曲線出現密集且連續的波峰波谷,堿解氮含量為[150,200)mg/kg和≥200 mg/kg處則呈現單調下降趨勢。在494 nm波段4條曲線均出現波谷,但波谷深度差異明顯,[50,100)mg/kg 吸收深度為 0.114,[100,150)mg/kg 吸收深度為 0.119,[150,200)mg/kg 吸收深度為0.257,≥200mg/kg吸收深度為0.336。以上分析說明,堿解氮含量與稻田土反射率光譜特征具有一定的相關性,稻田土光譜曲線在小于700 nm波段內存在敏感波段。

2.2 堿解氮含量與光譜參數的相關分析

為了分析43個樣本土壤光譜反射系數及其各種變換形式與土壤堿解氮含量在不同波段的相關性,將原始數據在ENVI軟件中進行了重采樣處理,對重采樣的光譜反射率R進行重采樣去包絡線處理Rc1,一階微分R',倒數,對數lnR,平方根R,立方根,倒數的一階微分()',對數的一階微分(lnR)',平方根的一階微分)',立方根的一階微分)',平方根的倒數,倒數的對數ln(),對數的倒數,對數倒數的一階微分()',倒數對數的一階微分(ln)',平方根倒數的一階微分)'等16種數學變換。表2為各種處理中相關系數絕對值最大的波段表。稻田土光譜曲線變換后相關系數見圖4。

表2 各種變換相關系數絕對值最大的波段匯總Tab.2 Sum of maximum band correlation absolute values of each transformation

圖4 稻田土光譜曲線變換后相關系數圖Fig.4 Correlation coefficient diagram of paddy soil spectrum curves after transformation

如圖4所示,堿解氮與稻田土反射率的一階微分的相關系數最高,最高值在2 194 nm,達0.646 7;對數倒數的一階微分在694 nm處相關系數也高達-0.632 4。結合表2和圖4可以看出,堿解氮與稻田土反射率相關性高的波段為644~694 nm和2 044~2 194 nm。這是因為這2個波段處在稻田土去包絡線吸收谷的邊緣,堿解氮含量高低不同,此波段內堿解氮反射率的變化率不同,即堿解氮含量與該波段反射率曲線的性狀特征相關性顯著。

綜合以上研究表明,稻田土光譜反射率受堿解氮含量影響的敏感波段為694 nm,2 058 nm和2 189 nm。

2.3 土壤堿解氮含量高光譜反演模型的構建

結合前文研究的敏感波段694 nm,2 058 nm和2 189 nm,選取光譜反射率重采樣對數倒數的一階微分作為變換形式用于構建反演模型。本研究中采用了偏最小二乘回歸(partial least squares regression,PLSR)計算,能有效地提取對系統解釋能力最強的綜合變量,排除無解釋作用的信息,使之對因變量有最強的解釋能力[17]。回歸模型的建模精度采用決定系數R2,標準偏差SD,均方根誤差RMSE以及標準偏差和均方根誤差的比值RPD來進行評價。R2越大,RPD值越大,說明反演模型的預測精度越好。PLSR計算公式為

式中:AN為土壤堿解氮含量;R694為694 nm波段稻田土反射率重采樣一階微分;R2058為2 058 nm波段稻田土反射率重采樣一階微分;R2189為2 189 nm波段稻田土反射率重采樣一階微分。稻田土堿解氮預測模型參數見表3。

表3 稻田土堿解氮預測模型參數表Tab.3 Index of paddy soil available nitrogen content prediction model

由表3可知,模型決定系數R2=0.56,RPD=1.51,說明模型穩定性較強,初步具備了預測能力。

3 結論

1)不同堿解氮含量的南方丘陵稻田土光譜曲線在波長<700 nm波譜范圍內隨著堿解氮含量的增高,光譜反射率降低,吸收深度加大的趨勢。

2)通過分析南方丘陵稻田土堿解氮含量與光譜反射率16種數學變換的相關系數,提取敏感波段為694 nm,2 058 nm和2 189 nm。

3)基于反射光譜特征的南方丘陵稻田土堿解氮高光譜反演模型穩定性強(R2=0.56),具備一定預測能力(RPD=1.51),能夠用于南方丘陵稻田土含量速測。

4)堿解氮含量與南方丘陵稻田土反射光譜特征呈現較好的相關性,基于反射光譜特征的南方丘陵稻田土堿解氮高光譜反演模型具備一定預測能力。本次研究樣點包含5種土屬,但個別土屬采樣點較少且采樣區域較集中,因而未從土屬對土壤光譜的影響的角度進行深入研究。影響土壤反射光譜特征的影響因素較多,如鐵、全氮、磷、有機質等化學指標。僅僅考慮堿解氮含量影響所建立的反演模型預測能力還受到一定的限制。因此,在下一步工作中,需要充分考慮其他因素的綜合影響,建立引入其他理化參數的南方丘陵稻田土養分反演模型。

[1] 黃應豐,劉騰輝.土壤光譜反射特性與土壤屬性的關系:以南方主要土壤為例[J].土壤通報,1989,20(4):158-160,176.Huang Y F,Liu T H.The relationship between soil properties and soil spectral reflectance characteristics:A case study of south main soil[J].Chinese Journal of Soil Science,1989,20(4):158-160,176.

[2] Fox G A,Sabbagh G J.Estimation of soil organic matter from red and near-infrared remotely sensed data using a soil line Euclidean distance technique[J].Soil Science Society of America Journal,2002,66(6):1922-1929.

[3] He Y,Song H Y,Pereira A G,etal.A new approach to predict N,P,K and OM content in a loamy mixed soil by using near infrared reflectance spectroscopy[C]//International Conference on Intelligent Computing,Lecture Notes in Computer Science Volume3644.Hefei,China:Springer-Verlag,2005:859-867.

[4] Kemper T,Sommer S.Estimate of heavy metal contamination in soils after a mining accident using reflectance spectroscopy[J].Environmental Science and Technology,2002,36(12):2742-2747.

[5] Deaton BC,Balsam W L.Visible spectroscopy-a rapid method for determining hematite and goethite concentration in geological materials[J].Journal of Sedimentary Research,1991,61(4):628-632.

[6] 劉煥軍,張 柏,劉殿偉,等.松嫩平原典型土壤高光譜定量遙感研究[J].遙感學報,2008,12(4):647-654.Liu H J,Zhang B,Liu D W,et al.Study on quantitatively remote sensing typical soils in Songnen plain,northeast China[J].Journal of Remote Sensing,2008,12(4):647-654.

[7] 魯如坤.土壤農業化學分析方法[M].北京:中國農業科技出版社,2000.Lu R S.Chemical Analysis of Soil Agriculture[M].Beijing:China Agricultural Science Press,2000.

[8] 王曉嵐,卡麗畢努爾,楊文念.土壤堿解氮測定方法比較[J].北京師范大學學報:自然科學版,2010,46(1):76-78.Wang X L,Kalibinuer,Yang W N.Comparison of methods for de-termining alkali-hydrolyzed nitrogen in soil[J].Journal of Beijing Normal University:Natural Science,2010,46(1):76-78.

[9] 何 挺,王 靜,程 燁,等.土壤氧化鐵光譜特征研究[J].地理與地理信息科學,2006,22(2):30-34.He T,Wang J,Chen Y,et al.Study on spectral features of soil Fe2O3[J].Geography and Geo-Information Science,2006,22(2):30-34.

[10] 徐 馳,曾文治,伍靖偉,等.內蒙古河套灌區土壤含鹽量和pH高光譜反演研究[J].灌溉排水學報,2013,32(3):39-43.Xu C,ZengW Z,Wu JW,et al.Retrieval of soil salt content and pH using hyper spectral remote sensing in Inner Mongolia[J].Journal of Irrigation and Drainage,2013,32(3):39-43.

[11] 胡 芳,藺啟忠,王欽軍,等.土壤鉀含量高光譜定量反演研究[J].國土資源遙感,2012(4):157-162.doi:10.6046/gtzyyg.2012.04.26.Hu F,Lin Q Z,Wang Q J,et al.Quantitative inversion of soil potassium content by using hyper spectral reflectance[J].Remote Sensing for Land and Resources,2012(4):157-162.doi:10.6046/gtzyyg.2012.04.26.

[12] 陳紅艷.土壤主要養分含量的高光譜估測研究[D].泰安:山東農業大學,2012.Chen H Y.Hyper spectral Estimation of Major Soil Nutrient Content[D].Taian:Shandong Agricultural University,2012.

[13] 徐麗華.土壤養分預測方法的比較研究——以三峽庫區王家溝小流域為例[D].重慶:西南大學,2012.Xu L H.Comparison of Prediction Methods on Soil Nutrients:A Case Study of Wangjiagou Watershed in the Three Gorges Reservoir Area[D].Chongqing:Southwest University,2012.

[14] Cloutis E A.Hyper spectral geological remote sensing:Evaluation of analytical techniques[J].International Journal of Remote Sensing,1996,17(12):2215-2242.

[15] Chang CW,Laird D A.Near-infrared reflectance spectroscopic analysis of soil C and N[J].Soil Science,2002,167(2):110-116.

[16] 周學文,趙小敏,胡國瑞,等.南方丘陵地區水田土壤養分變異分析[J].江西農業大學學報,2009,31(5):919-926.Zhou XW,Zhao X M,Hu G R,et al.The analysis of soil nutrient variability in hilly paddy field of southern China[J].Acta Agriculturae Universitatis Jiangxiensis,2009,31(5):919-926.

[17] 王惠文.偏最小二乘回歸方法及其應用[M].北京:國防工業出版社,1999:119.Wang HW.Partial Least Square Regression Method and Its Application[M].Beijing:National Defense Industry Press,1999:119.

主站蜘蛛池模板: 久久99精品久久久大学生| 丁香六月激情综合| 天天综合网亚洲网站| 米奇精品一区二区三区| 日韩精品一区二区三区中文无码| 国产精选小视频在线观看| 日韩少妇激情一区二区| 五月激情婷婷综合| 日本成人精品视频| 香蕉网久久| 黄色一及毛片| 免费又黄又爽又猛大片午夜| 精品一区二区久久久久网站| 99国产精品国产| 伊人久久综在合线亚洲91| 亚洲福利视频一区二区| 不卡国产视频第一页| 国产最爽的乱婬视频国语对白| 91福利一区二区三区| 婷婷激情亚洲| 亚洲最大看欧美片网站地址| 一级毛片无毒不卡直接观看| 久久久久夜色精品波多野结衣| 高清色本在线www| 一级成人a毛片免费播放| 亚洲va在线观看| 综合人妻久久一区二区精品| 久久综合结合久久狠狠狠97色| 在线免费a视频| 成人年鲁鲁在线观看视频| 青青草欧美| 国产老女人精品免费视频| 一级毛片免费的| 91成人免费观看在线观看| 亚洲经典在线中文字幕| 欧美国产在线一区| 凹凸精品免费精品视频| 中文字幕自拍偷拍| 爽爽影院十八禁在线观看| 国产成人1024精品| 性欧美精品xxxx| 国产无码网站在线观看| 免费无码AV片在线观看国产| 欧美精品成人一区二区视频一| 久草热视频在线| 色综合久久综合网| 日日拍夜夜操| 少妇精品久久久一区二区三区| 亚洲精品无码AⅤ片青青在线观看| 日本五区在线不卡精品| 日本欧美成人免费| 国产在线日本| 国产精品久久久久久久久久98| 国产精品自在拍首页视频8| 国产美女久久久久不卡| 熟女成人国产精品视频| 日本欧美中文字幕精品亚洲| 一级成人a毛片免费播放| 伊人久久大香线蕉影院| 在线观看av永久| 亚洲一区二区三区在线视频| 91久久国产综合精品| 日韩在线成年视频人网站观看| 国产在线精彩视频论坛| 欧美日韩亚洲国产| 无码内射中文字幕岛国片| 国产成人精品综合| 欧美黑人欧美精品刺激| 亚洲欧洲自拍拍偷午夜色无码| 色有码无码视频| 国产女人喷水视频| 最新国产精品第1页| 免费人成又黄又爽的视频网站| 无码aaa视频| 婷婷综合色| 国产成人三级| 中文字幕 欧美日韩| 99成人在线观看| A级毛片高清免费视频就| 日本人妻丰满熟妇区| 国产精品第一区| 久久精品欧美一区二区|