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

基于物候特征的農田土壤表層含水率遙感反演

2021-11-09 07:44:06張曉春劉海若嚴憶輝
灌溉排水學報 2021年10期
關鍵詞:模型

張曉春,劉海若,嚴憶輝,唐 蓉,張 煜

(1.武漢大學 水資源與水電工程科學國家重點實驗室,武漢430072;2.長江水利委員會長江科學院 空間信息技術應用研究所,武漢430010)

0 引 言

【研究意義】傳統的土壤含水率測定常以點代面進行定點測量,但取樣點的選取對結果有較大影響,能代表的實際范圍也比較有限[1]。而基于遙感影像監測土壤含水率具有監測時空范圍大、分析速度快和經濟高效等優點[2]。遙感監測土壤含水率方法主要分為光學遙感和微波遙感2 種監測方法。光學遙感監測有基于指數的方法,例如距平植被指數、植被狀態指數(VCI)、溫度植被干旱指數(TVDI)和垂直干旱指數(PDI)等,以及基于土壤熱慣量的方法[3]。光學遙感利用的波段主要為熱紅外和可見光—近紅外波段,對土層的穿透能力有限,一般只能監測表層幾微米的土壤含水率,且光學影像受云雨天氣的影響較大。微波遙感中的微波波段可以監測到土壤表層幾厘米到十幾厘米的土壤含水率,特別是星載合成孔徑雷達(SAR)的工作波段波長比一般微波波段波長長,且不受云雨遮蓋的影響,具有全天候全天時優點,但會受到地表植被覆蓋的影響[4-5]。可基于光學數據消除植被覆蓋的影響,并基于雷達數據反演得到土壤含水率分布圖,因此與使用單一數據相比,融合多種數據源能夠提高土壤含水率反演精度[6]。

【研究進展】利用光學數據消除植被覆蓋影響的模型主要是半經驗模型,其中水云模型應用最為廣泛。Attema 等[7]以農作物為研究對象,通過研究農作物覆蓋地表植被后向散射特性,提出估算農作物覆蓋地表土壤水分的“水云”模型(Water-Cloud Model)。水云模型將輻射傳輸模型中的植被冠層概化成水平均勻的云層,不考慮植被和土壤表層之間的多次散射,重要的變量僅為冠層高度和云密度,因而可假設與冠層體積含水率成比例,能夠有效反演耕地土壤含水率。Baghdadi 等[8]基于哨兵一號衛星影像和光學影像,在水云模型中添加NDVI數據作為輸入數據,研究不同極化條件下雷達信號數據的差異性。曾旭婧等[9]基于Sentine-1A(哨兵一號)雙極化雷達影像和Landsat 8光學影像,對不同植被覆蓋條件下的土壤含水率進行研究,發現VV 極化圖像和雙極化組合在低植被覆蓋條件下更具有優勢。胡佩敏等[10]利用ASAM GM 和MODIS 數據監測湖北省四湖地區的土壤表層含水率的時空分布,結果表明該方法實際可行。張友靜等[11]利用ASAM APP 和TM 數據,基于水云模型下反演小麥覆蓋條件下的土壤含水率,與實測數據相比,其相關性可達90%以上。綜上所述,基于水云模型監測農作物覆蓋區的土壤含水率是可行的。

【切入點】Zribi 等[12]研究表明土壤相對體積含水率和遙感反演的后向散射系數具有良好的線性相關性。目前研究大多數是基于某一年的點測土壤含水率,通過線性回歸法建立后向散射系數與土壤含水率的關系,并基于遙感反演的后向散射系數獲取土壤含水率分布圖[13-14]。然而不同的作物在不同時期具有不同的需水量,不同物候時期的降水量也不同,且作物需水量和區域降水量直接關系到土壤含水率,因此需進一步開展農田土壤水分研究[15],考慮作物物候特征,將土壤后向散射系數轉換為土壤含水率。

【擬解決的關鍵問題】選用Sentine-1 雷達影像和環境衛星CCD 光學數據,對研究區表層土壤含水率進行反演監測,并輔以landsat 8 衛星對環境衛星影像進行幾何校正。本研究首先對遙感影像進行預處理,從Sentine-1 SAR 影像提取總體后向散射系數,并從環境衛星影像提取NDVI植被指數,基于水云模型消除植被覆蓋的影響,推求出土壤的后向散射系數。然后根據作物的物候特征,利用土壤后向散射系數與土壤含水率之間的線性關系分生育階段推求土壤水分,最后用點測土壤含水率對遙感反演值進行驗證。

1 研究區與數據來源

1.1 研究區概況

研究區域為安徽省蚌埠市固鎮縣和亳州市利辛縣(圖1)。固鎮縣氣候屬亞熱帶和暖溫帶過渡帶,氣候兼有南北之長,四季分明,光照充足,年平均氣溫14.9 ℃,降水量871 mm,日照時間2 170 h。利辛縣區域總面積約2 005 km2。利辛縣屬暖溫帶半濕潤季風氣候區,且有明顯的過渡性特征,氣候溫和,光照充足,雨量適中,但集中于夏季,冬季時間長并且較為干燥。安徽省區位優勢明顯,常年農作物種植面積超過860 萬hm2,為全國糧食主產省。其中,固鎮縣與利辛縣均為傳統農業大縣,種植業發達。2 個縣夏收和秋收主要作物為冬小麥和夏玉米,本研究在固鎮縣實測了冬小麥土壤含水率,在利辛縣實測了夏玉米土壤含水率,具有較長時間序列的實測數據,具有一定的代表性。

圖1 研究區地理位置Fig.1 Location of the study area

1.2 數據來源

1.2.1 遙感數據及預處理

1)哨兵一號數據

哨兵一號(Sentinel-1)衛星由哨兵1-A 和哨兵1-B組成,均載有C 波段合成孔徑雷達。Sentinel-1 單顆衛星重訪周期為12 d,2 顆衛星組成系統的重訪周期為6 d。哨兵一號衛星擁有干涉寬幅模式(IW)和波模式(WV)2 種主要工作模式,還有條帶模式(SM)和超寬幅模式(EW)2 種附加模式。哨兵數據來源于NASA EarthData 數據網站(https://earthdata.nasa.gov/)。研究選取2016年到2018年IW 模式和VV 極化條件下的110 幅景象。Sentinel-1 原始影像帶有因噪聲、輻射、地形起伏等因素造成的干擾,需對其進行預處理。研究使用SNAP Desktop 軟件對110 幅雷達影像進行處理,主要包括輻射校正、多視處理、斑點噪聲抑制、幾何校正、拼接和裁剪等步驟。

2)環境衛星數據

環境衛星系列衛星由2 顆光學衛星HJ-1A 和HJ-1B,以及1 顆雷達衛星HJ-1C 組成。其中A 星和B 星均搭載CCD 相機。環境衛星數據來源于中國科學院的地理空間數據云網站(http://www.gscloud.cn/)。研究選擇固鎮縣2015年到2018年期間的55 幅CCD影像和利辛縣2020年的11 幅CCD 影像作為NDVI數據源。預處理步驟包括輻射定標、大氣校正、幾何校正、圖像裁剪等。以Landsat 8 影像作為基準影像,對環境衛星影像進行幾何校正,其中Landsat 8 數據來自USGS 網站(https://earthexplorer.usgs.gov/)。以上遙感影像均采用UTM/WGS84 投影系統。

1.2.2 地面點測數據

在固鎮縣和利辛縣均采集了地面點土壤含水率數據,其中固鎮縣僅有冬小麥種植區田間土壤含水率數據,利辛縣僅有夏玉米種植區土壤含水率數據。

固鎮縣土壤含水率實測數據在新馬橋灌排試驗站(位置如圖2(a)所示)采集,采集方法是在試驗田選取3 個點,在地表下0~10 cm 處用環刀取土,用鋁盒封存。在實驗室用稱質量法測得土壤質量含水率。取3 次測量的平均值,利用環刀法測定土壤體積質量為1.36 g/cm3,并換算得到土壤體積含水率,測定的土壤田間持水率為40.6%。固鎮縣站點數據時間如表1 所示。

利辛縣實測數據的設備是智墑土壤監測設備,設備基于頻域反射法(Frequency Domain Reflectometry,FDR),通過測量電解質常量的變化測量土壤的體積含水率。研究區監測點一和監測點二分別設有智墑一號和智墑二號儀器,分別于2020年6月21日和8月24日安裝,按照1 h 間隔監測土壤表層10 cm 體積含水率,保證數據測量時間與衛星過境時間一致,監測點經緯度使用手持GPS 采集。利辛縣監測點數據時間如表1 所示。

表1 2016—2020年試驗區實測數據日期Table 1 Date of measured data in the test area from 2016 to 2020

以固鎮縣2017年和2018年實測土壤含水率作為擬合數據,2016年的實測數據作為遙感反演結果的驗證數據。利辛縣共有2 個監測點(位置如圖2(b)所示),采用較為新型的土壤監測設備測得土壤含水率,以利辛縣監測點一的數據作為擬合數據源,監測點二的數據作為遙感反演結果的驗證數據。

圖2 地面實測數據位置Fig.2 Location of ground measured data

1.2.3 種植結構數據

首先基于2017年與2020年的landsat 8 影像對固鎮和利辛區域進行監督分類,分出城鎮和水體部分,并通過掩膜提取后得到作物種植區。利用HJ 衛星提取2017年固鎮縣地區長時間序列影像以及2020年利辛縣地區長時間序列影像的NDVI值,分別結合冬小麥與夏玉米的NDVI特征曲線進行決策樹分類,篩選出符合特征的農作物種植區。通過查詢安徽省農業統計年鑒以及政府統計公告,得到2017年固鎮縣種植冬小麥的總面積以及2020年利辛縣種植夏玉米的總面積,對分類結果進行修正,最終得到固鎮縣冬小麥種植區域和利辛縣冬小麥種植區域分布圖。農業統計年鑒在中國知網的統計數據平臺下載(https://data.cnki.net/),政府統計公告從利辛縣人民政府官網得到(http://lxxxgk.bozhou.gov.cn/)。具體的區域提取精度如表2 所示。

2 研究方法

2.1 基于水云模型計算土壤后向散射系數

水云模型即植被散射理論模型,是Attema 等[7]基于輻射傳輸方程提出的經典方程。該模型認為雷達傳感器所接受的散射回波是由于植被和土壤之間散射而成,并假設植被中的植被層水分是均勻分布的。模型忽略植被與土壤表面的多次相互散射,模型中的變量僅為植被高度和植被含水率[16]。水云模型用于去除植被對后向散射系數的影響,計算式為:

式中:σ0(θ)為總的后向散射系數;為土壤表層后向散射系數;為植被散射系數;γ2(θ)為植被雙層衰減因子(透過率);θ為入射角(rad);A、B值取決于植被類型及電磁波頻率的經驗系數;mv為植被含水率(kg/m2)。

參數A、B值的大小與植被類型密切相關,不同植被覆蓋條件下的A、B值有著較大的差異。研究區植被覆蓋條件相對較復雜,主要包括小麥、玉米和油菜等,且各自分布不集中。參考Bindlish 等[17]總結所得到的在不同地表覆蓋條件下的水云模型參數,選取其中所有植被的綜合方式能更好地反映研究區的實際情況,即A=0.001 2,B=0.091。

選取Jackson[18]的經驗模型計算植被含水率:

式中:VMC為植被含水率;NDVI為歸一化植被指數。根據式(1)—式(4)和水云模型參數A、B值的選取,即可提取出在VV 極化條件下考慮作物覆蓋下的土壤后向散射系數:

計算土壤體積相對含水率和遙感反演的后向散射系數(分貝)的線性相關性時,需要將無單位的后向散射系數轉換成單位為分貝的后向散射系數:

式中:為土壤表層后向散射系數;為單位為分貝的土壤表層后向散射系數。

2.2 結合物候特征的土壤含水率遙感反演模型

2.2.1 基于后向散射系數的土壤含水率計算方程

本文土壤含水率(Soil Moisture,SM)統一用土壤相對體積含水率表示。土壤體積相對含水率和遙感反演的后向散射系數具有線性相關性[12]:

式中:C和D均為無單位的經驗系數;σ0soil-db為單位為分貝的土壤表層后向散射系數。

2.2.2 結合物候特征的時間分段方法

本文研究對象為農作物區,研究區域作物類型以及不同生長期需水的不同會對表層土壤含水率產生不同程度的影響。作物需水量和區域降雨量也直接關系到土壤含水率的值,因此需要根據不同的時間段并結合作物物候特征來分別確定式(7)中系數C和D值。研究區域種植作物主要為夏玉米和冬小麥,以下將結合2 種作物的物候特征介紹分段方法。

作物系數Kc是指作物不同發育期中需水量與可能蒸散量的比值。冬小麥不同生育期具有不同的作物系數Kc,作物系數變化能反映出作物需水量的變化,而土壤含水率與作物需水量具有相關性。因缺少淮北平原生育期作物系數的研究資料,且淮北平原與黃淮海平原地理位置接近,因此,用黃淮平原冬小麥生育期作物系數[19]近似代替淮北平原冬小麥作物系數。如表3 所示,研究區域冬小麥種植期一般從10月初開始一直到次年6月,生育期較長。在冬小麥的生長發育過程中,從播種期到分蘗期作物系數Kc值逐漸上升,需水量逐漸增多。從越冬期開始,冬小麥停止生長,作物系數Kc值降低,需水量達到低谷。返青期后氣溫逐漸恢復正常,冬小麥開始重新生長,作物需水量逐漸上升。根據冬小麥生育期作物需水特征的變化,將全生育期分為3 個階段:播種—分蘗期(10—12月)、越冬期(1—3月)、返青期—成熟期(4—6月)。在研究區實地田間調查也發現,一般6月開始種植夏玉米,到10月收割,接著11月種植冬小麥,第二年4月左右收割。研究區域夏玉米種植期一般從6月開始一直到9月底,生育期較短,將夏玉米全生育期劃分為一段。

表3 冬小麥全生育期作物系數(Kc)的變化[19]Table 3 Variation of crop coefficient(Kc)during the whole growth period of winter wheat

根據往年實測點土壤含水率,發現以上劃分的冬小麥3 個時間段和夏玉米全生育期內相應的土壤含水率變化規律較為明顯。冬小麥從播種開始生長發育,此時冬季降水少且作物需水量也小,土壤含水率維持在較低水平,隨后進入越冬期。越冬期內作物生長緩慢,需水量低,土壤含水率變化平穩。從返青期開始作物生長加快,作物需水量變大,土壤含水率迅速降低并持續最低。而夏玉米全生育期內的土壤含水率基本維持在比較均衡的水平。同時,在氣候變化背景下,作物的生育期也存在年際變化。就華北平原地區冬小麥而言,在1961—2017年中各個生育期均有推遲變化的趨勢,平均推遲3~7 d,但各生育期所在月份幾乎不會產生變化[20]。故在以月為單位的函數中,作物物候特征及生育期變化可以忽略。

因此,本文將全年共劃分為4 個時間段來分別通過土壤后向散射系數計算土壤含水率。4 個時間段為:冬小麥的3 個生育期播種—分蘗期(10—12月)、越冬期(1—3月)、返青期—成熟期(4—6月),和夏玉米的全生育期(6—9月)。

2.2.3 結合物候特征的土壤含水率分生育階段估算模型

土壤含水率遙感反演模型采用分生育階段函數方法來轉換土壤后向散射系數,分別為冬小麥越冬期(overwintering period,OW)、冬小麥返青—成熟期(growing stage from reviving to maturity,RM)、夏玉米全生育期(whole grow phase of summer-sown corn,SC)、冬小麥播種—分蘗期(growing stage from sowing to tillering,ST),分生育階段函數為:

式中:Ci、Di分別為對應生育期的經驗系數,其中i=1、2、3、4;SMST為冬小麥播種—分蘗期土壤含水率;SMOW為冬小麥越冬期土壤含水率;SMRM為冬小麥返青—成熟期土壤含水率;SMSC為夏玉米全生育期土壤含水率;σ0soil-db是單位為分貝的土壤表層后向散射系數。

3 結果與分析

3.1 去除植被前后的后向散射系數對比

通過水云模型去除地表植被覆蓋影響后,后向散射系數會衰減,并得到更為準確的土壤后向散射系數。以固鎮縣2017年數據為例,去除植被影響前后的后向散射系數如圖3 所示。圖中實線為去除前的后向散射系數,虛線為去除后的后向散射系數,去除后的后向散射系數整體比去除前的值小。在利用水云模型分離出植被層散射和吸收的貢獻后,地表后向散射系數均有衰減,如圖3 所示,基于水云模型去除后的后向散射系數比去除前的減小5 db 左右,去除效果較好。

圖3 去除植被影響前后效果Fig.3 Effect before and after removing vegetation

3.2 結合物候特征的土壤含水率遙感反演經驗系數

基于2017年和2018年固鎮縣地區土壤含水率實測值和遙感計算的后向散射系數值擬合冬小麥3 個生育期的經驗系數值,基于2020年利辛縣監測點一的土壤含水率實測值和遙感計算的土壤后向散射系數值擬合夏玉米生育期的經驗系數值,并對農作物各生育期后向散射系數與土壤相對體積含水率數據進行擬合后做統計分析,得到復相關系數R2、平均絕對誤差MAE和均方根誤差RMSE,最后得到遙感反演模型的系數如表4 所示。

表4 研究區不同作物生育期的系數值Table 4 Coefficient values of different crop growth periods in the study area

如表4 所示,固鎮縣冬小麥播種期—分蘗期時間段內模擬精度不高,原因在于此時間段內地表幾乎無植被覆蓋,冬小麥根系開始生長但此時根長密度較低,導致土壤含水率受實時天氣等因素影響較大,模擬精度降低。具體的分生育階段函數為:

3.3 模型驗證與精度評價

基于實測土壤含水率數據,對土壤含水率遙感反演模型進行驗證和精度評價,如圖4 所示。

冬小麥分生育階段土壤含水率模型以2015年11月至2016年5月固鎮縣新馬橋試驗點的數據進行驗證。如圖4(a)所示,計算得2 組數據復相關系數R2為0.73,平均絕對誤差MAE=0.051,最大誤差MaxE=0.072,最小誤差MinE=0.025,均方根誤差RMSE=0.056。夏玉米土壤含水率模型以2020年利辛縣監測點二采集的數據進行驗證。如圖4(c)所示,計算得2 組數據復相關系數R2達到0.82,平均絕對誤差MAE=0.069,最大誤差MaxE=0.112,最小誤差MinE=0.027,均方根誤差RMSE=0.078。2 種作物土壤含水率模擬精度滿足要求。

如圖4(b)所示,實測土壤含水率與反演土壤含水率總體趨勢幾乎相同,但存在差距。根據相關新聞(http://www.weather.com.cn/anhui/tqyw/01/2642620.s html)與研究資料[21],安徽省2016年降水量為55 a來最多,4月起降水明顯增多,4—7月安徽省平均降雨較往年偏多近6 成,且6月左右暴雨洪澇災害嚴重,因此實際土壤含水率較往年偏高,反演值較實測值偏低。

圖4 遙感反演精度驗證Fig.4 Verification of remote sensing inversion accuracy

3.4 土壤含水率反演結果

根據上述分生育階段函數反演了固鎮縣2016—2018年農作區的土壤含水率空間分布,以及利辛縣2020年農作區的土壤含水率分布(圖5),在2017年固鎮縣冬小麥3 個生育期挑選3 幅土壤相對體積含水率圖,以及2020年利辛縣夏玉米生育期中的一幅土壤相對體積含水率圖,作為示例展示。

圖5 農作物區土壤相對體積含水率Fig.5 Relative volumetric water content of each growth period

如圖5(a)—圖5(c)所示,從時間上冬小麥越冬期土壤含水率較低,其他生育期土壤含水率明顯高于越冬期土壤含水率,與冬小麥生育期需水特性符合。從空間上,固鎮縣東南方的土壤含水率略高于西北方,與固鎮縣自西北向東南傾斜的地勢吻合。如圖5(d)所示,從空間上,利辛縣東南方的土壤含水率略高于西北方,與利辛縣西北略高,東南略低的地勢吻合。

4 討論

傳統的烘干法監測土壤含水率是各國公認的最為準確且經典的方法,但監測費時費力且具有很大的空間局限性。利用已建立的模型,可直接得到農田區域的土壤含水率,在時間和空間尺度上都有更好的發展。在未來的研究中,利用大面積反演的土壤含水率,結合區域水文模型反演出農田中水分變化,進而求得農田澇漬災害,對農業生產和災害評估進行指導。

基于水云模型結合物候特征反演土壤表層含水率,在地表植被覆蓋的條件下擁有較好精度。本研究中,安徽省固鎮縣冬小麥各生育期土壤含水率變化不一致,其中返青—拔節期土壤含水率變化最大,越冬期土壤含水率變化最小,同時冬小麥生育期出現日期與持續時間的年際波動大,年內變化趨勢不明顯的結果,與張佩[22]等針對淮北地區冬小麥發育期特征的研究結論相一致。這也是用不同年同測站的數據進行驗證精度不高,用同年不同測站的數據進行驗證精度較好的原因。由圖5 可知,反演的土壤含水率空間分布結果與區域地勢和地區氣候有較好的一致性,在地勢較低的地方含水率高,地勢較高的地方含水率相對較低,這進一步表明,基于遙感數據反演的土壤含水率空間分布結果可靠。

固鎮縣冬小麥播種時間大多在10月初,在剛播種的一段時間內,農作物區土壤表面很少或無植被覆蓋而形成裸地,且通過遙感反演的土壤表層含水率隨天氣波動大,進一步削弱遙感數據反演的精度。相反,在較多植被覆蓋的生育期內土壤含水率保持較為平穩,反演精度高。實際應用中,針對冬小麥播種期—分蘗期的生育階段,在接下來的工作中還需進一步單獨進行研究,選擇不同的土壤水分反演方式提高反演精度。從驗證結果也可看出,結合物候特征的土壤表層含水率反演模型具有較好的精度,盡管用于模型反演的地面實測數據與區域較為有限,且影響反演精度的因素較多,但結果仍具有一定代表性,能為相關研究提供一些思路與參考。

5 結論

1)通過遙感數據計算2017年和2018年固鎮縣冬小麥種植區以及2020年利辛縣夏玉米種植區的后向散射系數,結合作物的不同物候特征分生育階段得到后向散射系數與土壤相對體積含水率之間的關系,相關系數分別達到0.80、0.91、0.79 和0.40,反演效果較好。

2)利用實測土壤相對體積含水率對模型進行驗證,結果表明該方法反演的土壤含水率在一般降雨年份內具有較好的模擬精度,相關系數達到0.73 和0.82,可靠性較高。

3)本研究根據不同作物的物候特征進行分生育階段求解可以提高模擬的精確度,為在農作物覆蓋條件下進行土壤含水率的反演提供了一定的依據,在農業生產指導以及澇漬災害預警方面具有一定的應用前景。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 婷婷成人综合| 自拍偷拍欧美日韩| 天天综合色天天综合网| 久久精品亚洲热综合一区二区| 最新日韩AV网址在线观看| 四虎亚洲精品| 精品国产一区91在线| 亚洲最大情网站在线观看| 朝桐光一区二区| 欧美一级在线| 精品国产一二三区| 狠狠色综合网| 日本黄色不卡视频| 精品乱码久久久久久久| 久久中文电影| 亚洲无线一二三四区男男| 91久久夜色精品| 亚洲色图欧美视频| 午夜国产小视频| 午夜一区二区三区| 亚洲人成高清| 国产精品真实对白精彩久久| 国产网站免费观看| 高清欧美性猛交XXXX黑人猛交| 亚洲视频无码| 中日韩一区二区三区中文免费视频| 四虎国产精品永久在线网址| 国产男女XX00免费观看| 熟妇丰满人妻| 91成人在线免费视频| 国产免费羞羞视频| 日本在线免费网站| 一本大道无码日韩精品影视| 精品日韩亚洲欧美高清a| 国产自在线播放| 日韩欧美国产成人| 夜夜操天天摸| 激情亚洲天堂| 91精品啪在线观看国产60岁 | 黄网站欧美内射| 91在线丝袜| 久久久噜噜噜| 超清无码一区二区三区| 2020精品极品国产色在线观看| 国产精品白浆在线播放| 国产精品第一区| 日本三区视频| 666精品国产精品亚洲| 欧美一级专区免费大片| 国内丰满少妇猛烈精品播| 国产精品lululu在线观看| 亚洲永久视频| 亚洲看片网| 色窝窝免费一区二区三区| 国产日韩久久久久无码精品| 91国内外精品自在线播放| 欧美在线三级| 国产aⅴ无码专区亚洲av综合网| 一级毛片在线直接观看| 久久久亚洲国产美女国产盗摄| 国产爽爽视频| 全部毛片免费看| 久久精品国产91久久综合麻豆自制| 新SSS无码手机在线观看| 日韩欧美国产成人| 91久久国产综合精品| 国产无套粉嫩白浆| 日韩精品久久久久久久电影蜜臀| 一级毛片a女人刺激视频免费| 久久无码免费束人妻| 欧美色图第一页| 欧美综合区自拍亚洲综合绿色| 亚洲欧美日韩成人在线| www亚洲天堂| 国产va在线观看免费| 国产成人区在线观看视频| 一本色道久久88综合日韩精品| 一本色道久久88| 亚洲欧州色色免费AV| 欧美午夜理伦三级在线观看| 精品无码国产一区二区三区AV| 日韩人妻无码制服丝袜视频|