亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        DEM分辨率對(duì)分布式水文模擬的影響

        2013-06-07 07:17:17李致家章玉霞
        水利水電科技進(jìn)展 2013年5期
        關(guān)鍵詞:模型

        姚 成,李致家,章玉霞

        (河海大學(xué)水文水資源學(xué)院,江蘇南京 210098)

        DEM分辨率對(duì)分布式水文模擬的影響

        姚 成,李致家,章玉霞

        (河海大學(xué)水文水資源學(xué)院,江蘇南京 210098)

        以數(shù)字高程模型為基礎(chǔ),采用柵格型新安江模型進(jìn)行不同分辨率情況下的洪水模擬,分析了分辨率對(duì)分布式水文模擬的影響。研究結(jié)果表明:數(shù)字高程模型分辨率越高,流域數(shù)字化的精度相對(duì)越高,但模型的模擬結(jié)果并不一定越好;針對(duì)不同分辨率的數(shù)字高程模型,在重新率定匯流參數(shù)的情況下,柵格型新安江模型的模擬精度相當(dāng),均能取得較好的應(yīng)用效果;對(duì)于中尺度流域而言,采用30″分辨率的數(shù)字高程模型不僅能夠提高模型的運(yùn)行時(shí)效,也能取得較高的模擬精度。

        數(shù)字高程模型;柵格型新安江模型;水文模擬;分辨率;參數(shù)估計(jì)

        流域水文模型作為一種模擬流域內(nèi)復(fù)雜水文現(xiàn)象的重要工具,一直是水文學(xué)研究的熱點(diǎn)課題之一。就目前水文模型發(fā)展趨勢(shì)而言,分布式水文模型由于其能夠更好地考慮氣候和下墊面因子的空間變異性,能夠更好地利用地理信息系統(tǒng)(GIS)與遙感(RS)技術(shù)等空間信息模擬流域的降雨徑流響應(yīng),現(xiàn)已成為水文模型發(fā)展的重要方向[1-2]。

        在應(yīng)用分布式水文模型時(shí),往往需要流域內(nèi)部詳細(xì)的水流流向、河網(wǎng)水系、流域形狀等地形地貌特征,需要對(duì)研究流域進(jìn)行數(shù)字化。隨著GIS的應(yīng)用與計(jì)算機(jī)能力的提高,利用數(shù)字高程模型(DEM)數(shù)字化流域和提取流域地貌特征的技術(shù)已得到了廣泛研究,DEM的引入大幅度促進(jìn)了分布式水文模型的進(jìn)一步發(fā)展[3]。然而,在DEM平臺(tái)上構(gòu)建分布式水文模型,DEM分辨率對(duì)流域地貌特征的提取結(jié)果以及產(chǎn)匯流模擬精度有著重要的影響[4-5]。理論上, DEM分辨率越高,所提取的流域地貌特征越接近于真實(shí)情況,但對(duì)基于DEM的分布式水文模型而言,高分辨率的DEM會(huì)使得模型的計(jì)算負(fù)擔(dān)過重,在很大程度上制約了模型的運(yùn)行時(shí)效,這必將導(dǎo)致分布式水文模型在實(shí)際應(yīng)用中受到限制。為了更真實(shí)高效地描述現(xiàn)實(shí)世界中的降雨徑流形成機(jī)理,有必要進(jìn)一步開展DEM分辨率對(duì)分布式水文模擬影響的研究。

        本文以皖南山區(qū)呈村流域?yàn)槔?在3″、12″、21″與30″共4種DEM分辨率情況下,對(duì)研究流域進(jìn)行了數(shù)字化,分析了分辨率對(duì)地貌特征提取結(jié)果的影響,并在此基礎(chǔ)上,采用基于DEM的柵格型新安江模型[6-8]進(jìn)行呈村流域的洪水模擬,探討了不同DEM分辨率對(duì)模型參數(shù)以及模擬結(jié)果的影響。

        1 柵格型新安江模型

        1.1 模型原理

        柵格型新安江模型以流域內(nèi)每個(gè)DEM柵格作為計(jì)算單元,并假設(shè)降雨和下墊面特征在每個(gè)單元內(nèi)空間分布均勻,模型只考慮各要素在不同單元之間的變異性。模型構(gòu)建以蓄滿產(chǎn)流原理為基礎(chǔ),其蒸散發(fā)、產(chǎn)流與分水源計(jì)算方法和新安江模型一致,采用概念性方法。因此,柵格型新安江模型屬于松散結(jié)構(gòu)分布式水文模型[1]。

        在模擬洪水時(shí),對(duì)于任意柵格單元,在考慮柵格間水量交換與河道排水網(wǎng)絡(luò)影響的基礎(chǔ)上,分別采用3層蒸散發(fā)模型與蓄滿產(chǎn)流原理對(duì)單元內(nèi)的蒸散發(fā)E與產(chǎn)流量R進(jìn)行計(jì)算,并采用自由水蓄水庫(kù)結(jié)構(gòu)對(duì)水源進(jìn)行劃分,包括地表徑流Rs、壤中流Ri以及地下徑流Rg3種水源,再根據(jù)柵格間的匯流演算次序,依次將各種水源演算至流域出口。對(duì)于單元間的匯流過程模擬,采用的是逐柵格馬斯京根匯流演算方法[6-9]。圖1為柵格單元內(nèi)產(chǎn)流與分水源示意圖,圖中P為降雨量;WM、SM分別為張力水、自由水蓄水容量;WUM、WLM分別為上層、下層張力水蓄水容量;S為實(shí)際的自由水含量;Kg、Ki分別為地下水徑流、壤中水徑流的出流系數(shù);θwp、θfc、θs分別為凋萎含水量、田間持水量與飽和含水量。

        圖1 柵格單元內(nèi)產(chǎn)流與分水源示意圖

        1.2 模型參數(shù)

        柵格型新安江模型以新安江模型為基礎(chǔ),其參數(shù)與新安江模型參數(shù)基本一致,如圖1中的WM、SM、Kg與Ki等。不同的是,當(dāng)利用該模型進(jìn)行分布式水文模擬時(shí),不能完全以集總的方式對(duì)模型參數(shù)賦值,而是要獲得參數(shù)的空間分布。因此,不能單一地根據(jù)觀測(cè)資料來率定參數(shù),而應(yīng)該在率定參數(shù)之前,先根據(jù)參數(shù)的物理意義,基于地形地貌、土壤、植被等下墊面條件,結(jié)合理論推導(dǎo),合理估計(jì)出參數(shù)的空間分布[10-12]。客觀的估計(jì)方法不僅有利于甄選出敏感參數(shù),降低參數(shù)之間的相依性,簡(jiǎn)化參數(shù)的率定工作,還可以保證模型參數(shù)在整個(gè)流域上的相容性,以更真實(shí)地模擬流域的水文響應(yīng)[13]。

        以WM、SM、Kg與Ki4個(gè)參數(shù)為例,WM表示的是張力水蓄水容量,對(duì)應(yīng)的土層為“全土層”,即整個(gè)包氣帶[14-15],而張力水可以認(rèn)為介于凋萎含水量與田間持水量之間,則有:

        式中:La為包氣帶土層厚度。SM表示的是自由水蓄水容量,對(duì)應(yīng)的土層為“表土層”,即腐殖質(zhì)土層[14-15],而自由水可以認(rèn)為介于田間持水量與飽和含水量之間,則有:

        式中:Lh為腐殖質(zhì)土層厚度。

        利用式(1)和式(2)計(jì)算每個(gè)柵格單元蓄水容量時(shí),土壤水分常數(shù)θwp、θfc、θs可以根據(jù)每個(gè)單元對(duì)應(yīng)的土壤類型估算得到,土層厚度La、Lh可以根據(jù)每個(gè)單元對(duì)應(yīng)的地形指數(shù)及植被覆蓋估算得到。利用式(1)和式(2)不僅可以獲得WM與SM的空間分布,也可以統(tǒng)計(jì)得到研究流域張力水與自由水容量曲線。

        兩個(gè)出流系數(shù)Kg與Ki屬于并聯(lián)參數(shù),其和代表自由水出流的快慢,與單元內(nèi)的土壤類型有關(guān),可以將θs與θfc作為衡量自由水出流快慢的指標(biāo);Kg/Ki代表地下水與壤中流的比,該比值可以通過θwp來反映[16],則有:

        式中:moc與mr為兩個(gè)校正系數(shù),根據(jù)已有的研究成果,在估計(jì)參數(shù)空間分布時(shí),mr≈1[13-16],而moc可以通過分水源中的結(jié)構(gòu)性約束(即Kg+Ki=0.7)試算得到。

        當(dāng)采用逐柵格馬斯京根法進(jìn)行單元間匯流演算時(shí),對(duì)于其參數(shù)流量比重因子xe與蓄量常數(shù)ke的確定,模型暫沒有考慮其空間分布不均的問題,而是利用流域內(nèi)柵格單元統(tǒng)一賦值的方法,需要根據(jù)流域的觀測(cè)資料進(jìn)行率定和檢驗(yàn)。ke屬于敏感參數(shù),反映的是河道水流、地表徑流、壤中流與地下徑流在柵格單元內(nèi)的匯流時(shí)間,由于不同的徑流成分在匯流速度上存在著差異,因此在率定ke時(shí)應(yīng)予以區(qū)分。

        2 研究流域概況

        呈村流域位于皖南山區(qū),鄰近東南沿海,四季分明,氣候溫和,屬亞熱帶季風(fēng)氣候區(qū),流域控制面積290km2。流域內(nèi)植被良好,雨量充沛,多年平均降雨量約為2 100 mm,流域降水在年內(nèi)年際分配極不均勻,為典型的濕潤(rùn)流域。呈村流域的植被類型主要包括常綠針葉林、落葉闊葉林、混合林、森林地、林地草原、牧草地與作物地,土壤類型主要為黏壤土。流域內(nèi)有10個(gè)雨量站,具有1986—1999年32場(chǎng)洪水過程的時(shí)段資料,其中,前22場(chǎng)洪水過程用于模型參數(shù)的率定,后10場(chǎng)洪水過程用于模型參數(shù)的檢驗(yàn)。

        3 結(jié)果分析

        3.1 地貌特征提取結(jié)果

        根據(jù)3″分辨率的DEM數(shù)據(jù)[17],利用最近鄰算法生成12″、21″與30″的DEM數(shù)據(jù),在這4種分辨率情況下,對(duì)呈村流域進(jìn)行了數(shù)字化,并提取了流域地貌特征。圖2為根據(jù)上述4種分辨率的DEM生成的流域邊界及水系,生成水系時(shí),對(duì)應(yīng)圖2(a)~(d)所用的閾值(網(wǎng)格數(shù))分別為100、25、15和10。

        圖2 根據(jù)不同分辨率DEM生成的呈村流域邊界與水系

        由圖2可以看出,DEM分辨率越高,所生成的流域邊界及水系越精細(xì)。分析地貌特征結(jié)果后發(fā)現(xiàn),分辨率的變化對(duì)所提取的呈村流域面積影響相對(duì)較小,變化幅度在1 km2以內(nèi)。此外,隨著DEM分辨率的降低,流域平均高程及河道平均坡度有增大的趨勢(shì),而最大匯流路徑長(zhǎng)度有減小的趨勢(shì),三者變化幅度分別在13 m、0.6%、12 km以內(nèi)。

        3.2 洪水模擬結(jié)果

        以流域數(shù)字化與地貌特征提取結(jié)果為基礎(chǔ),將柵格型新安江模型用于呈村流域的洪水模擬??紤]到分辨率對(duì)模型運(yùn)算效率有較大影響,先針對(duì)30″DEM,結(jié)合土壤、植被等下墊面資料,對(duì)模型產(chǎn)流參數(shù)的空間分布進(jìn)行了估計(jì),并在此基礎(chǔ)上,對(duì)模型匯流參數(shù)進(jìn)行了率定,其中,河道水流、地表徑流、壤中流和地下徑流對(duì)應(yīng)的ke率定結(jié)果分別為0.06、0.10、0.60和1.00,xe的率定結(jié)果均為0.3。然后在不改變模型參數(shù)的情況下對(duì)其他幾個(gè)不同分辨率DEM進(jìn)行了應(yīng)用,結(jié)果表明,完全根據(jù)30″DEM對(duì)應(yīng)的參數(shù)模擬呈村流域的洪水,對(duì)于3″、12″與21″DEM而言,峰值會(huì)減小,峰現(xiàn)時(shí)間也會(huì)滯后。分辨率越高,即柵格單元越小,洪水過程的推移坦化現(xiàn)象越明顯。分析其原因,可能主要是由于匯流參數(shù)ke取值未作改變所造成的。ke反映了不同水源在柵格單元上的匯流時(shí)間,隨著DEM分辨率的提高,計(jì)算單元尺度減小,相應(yīng)的匯流時(shí)間應(yīng)有所減小,因此需要對(duì)ke進(jìn)行重新率定。通過對(duì)ke重新率定后發(fā)現(xiàn),若以30″DEM對(duì)應(yīng)的ke率定結(jié)果為準(zhǔn),不同分辨率對(duì)應(yīng)的k′e≈ke/n(n為DEM分辨率變化的倍數(shù)),即可以根據(jù)計(jì)算單元的尺度對(duì)該參數(shù)進(jìn)行校正。表1為只率定ke而不改變其他參數(shù)情況下,不同分辨率DEM所對(duì)應(yīng)的柵格型新安江模型模擬結(jié)果的合格率、平均相對(duì)誤差水平以及確定性系數(shù)統(tǒng)計(jì),各精度評(píng)定指標(biāo)均根據(jù)GB/T 22482—2008《水文情報(bào)預(yù)報(bào)規(guī)范》計(jì)算得到。由表1可以看出,通過對(duì)匯流參數(shù)ke進(jìn)行尺度校正,模型均能取得較高的模擬精度,表明流域匯流速度與總的匯流時(shí)間不受計(jì)算單元尺度的影響。

        表1 不同分辨率DEM對(duì)應(yīng)的柵格型新安江模型模擬結(jié)果誤差統(tǒng)計(jì)

        表1和圖2的應(yīng)用結(jié)果表明,雖然利用高分辨率DEM能夠獲得高精度的數(shù)字流域結(jié)果,但并不一定能夠獲得更高精度的分布式水文模擬結(jié)果。由表1可以看出,在3″、12″、21″與30″這4種DEM分辨率情況下,柵格型新安江模型的模擬精度相當(dāng),采用30″DEM不僅能夠提高模型的運(yùn)行時(shí)效,也能取得較好的應(yīng)用效果。此外,模型在進(jìn)行柵格單元產(chǎn)流計(jì)算時(shí),其參數(shù)可以根據(jù)每個(gè)單元內(nèi)的土壤類型、植被覆蓋等下墊面特征進(jìn)行提取,在分析DEM分辨率對(duì)模型影響時(shí),所用的土壤及植被類型數(shù)據(jù)并沒有改變,因此提取的參數(shù)分布變化不大。以WM與SM兩個(gè)參數(shù)為例,圖3為不同分辨率DEM對(duì)應(yīng)的蓄水容量分布曲線(橫坐標(biāo)為不同的WM與SM值所占的流域面積比例),該曲線是利用式(1)和式(2)計(jì)算每個(gè)單元對(duì)應(yīng)的WM與SM后統(tǒng)計(jì)出的兩者在流域上總的分布情況。圖3表明不同分辨率DEM對(duì)應(yīng)的張力水和自由水蓄水容量曲線相似度較高,總體趨勢(shì)基本一致。由此可以看出,采用客觀合理的參數(shù)估計(jì)方法能夠在一定程度上降低參數(shù)對(duì)計(jì)算單元尺度的依賴性。

        圖3 不同分辨率DEM對(duì)應(yīng)的蓄水容量分布曲線

        4 結(jié) 語(yǔ)

        隨著計(jì)算機(jī)、地理信息系統(tǒng)、遙感等新興技術(shù)以及水文學(xué)理論的進(jìn)一步發(fā)展,基于DEM的分布式水文模型成為水文模型新的發(fā)展方向。本文采用的柵格型新安江模型,即是以DEM為基礎(chǔ)構(gòu)建的松散結(jié)構(gòu)分布式模型,該模型可以較好地考慮降雨和地形地貌、土壤、植被等下墊面條件的空間變異性。

        考慮到在實(shí)際應(yīng)用與研究中,DEM分辨率對(duì)分布式模型有著重要的影響,以皖南山區(qū)呈村流域?yàn)槔?針對(duì)3″、12″、21″與30″這4種分辨率的DEM,探討了分辨率對(duì)柵格型新安江模型的影響。研究結(jié)果表明,在不同分辨率情況下進(jìn)行呈村流域的洪水模擬,通過對(duì)匯流參數(shù)ke進(jìn)行尺度校正,模型均可以取得較高的模擬精度,利用土壤、植被等下墊面特征估算模型參數(shù)有利于降低參數(shù)的尺度依賴性;此外, DEM分辨率越高,流域數(shù)字化的精度相對(duì)越高,但模型的應(yīng)用效果并不一定越好。一般而言,分辨率過高,會(huì)導(dǎo)致分布式模型的計(jì)算負(fù)擔(dān)過重,在一定程度上會(huì)降低其實(shí)用性。對(duì)于中尺度流域而言,采用30″DEM不僅能夠提高模型的運(yùn)行時(shí)效,也能取得較好的應(yīng)用效果。

        柵格型新安江模型在實(shí)際應(yīng)用時(shí),只考慮了降雨與下墊面特征在不同計(jì)算單元之間的空間變異性,而在每個(gè)單元內(nèi),則認(rèn)為各個(gè)要素空間分布均勻,如何更好地考慮其在計(jì)算單元內(nèi)的空間變異性,還有待更深入的研究。

        [1]芮孝芳,黃國(guó)如.分布式水文模型的現(xiàn)狀與未來[J].水利水電科技進(jìn)展,2004,24(2):55-58.(RUI Xiaofang, HUANG Guoru.The present and future of the distributed hydrological models[J].Advances in Science and Technology of Water Resources,2004,24(2):55-58.(in Chinese))

        [2]熊立華,郭生練.分布式流域水文模型[M].北京:中國(guó)水利水電出版社,2004.

        [3]芮孝芳,石朋.數(shù)字水文學(xué)的萌芽及前景[J].水利水電科技進(jìn)展,2004,24(6):55-58.(RUI Xiaofang,SHI Peng.Seeds and prospects of digital hydrology[J]. Advances in Science and Technology of Water Resources, 2004,24(6):55-58.(in Chinese))

        [4]VIEUX B.Distributed hydrologic modeling using GIS [M].Dordrecht:Kluwer Academic Publishers,2001.

        [5]BEVEN K J.Rainfall-runoff modelling:the primer[M]. West Sussex,UK:John Wiley&Sons,2012.

        [6]李致家,姚成,汪中華.基于柵格的新安江模型的構(gòu)建和應(yīng)用[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,2007,35(2): 131-134.(LI Zhijia,YAO Cheng,WANG Zhonghua. Development and application of grid-based Xinanjiang model[J].Journal of Hohai University:Nature Sciences, 2007,35(2):131-134.(in Chinese))

        [7]李致家,姚成,章玉霞,等.柵格型新安江模型的研究[J].水力發(fā)電學(xué)報(bào),2009,28(2):25-33.(LI Zhijia, YAO Cheng,ZHANG Yuxia,et al.Study on gird-based Xinanjiangmodel[J].JournalofHydroelectric Engineering,2009,28(2):25-33.(in Chinese))

        [8]YAO Cheng,LI Zhijia,BAO Hongjun,et al.Application ofadevelopedGrid-XinanjiangmodeltoChinese watersheds for flood forecasting purpose[J].Journal of Hydrologic Engineering,2009,14(9):923-93.

        [9]姚成,章玉霞,李致家.擴(kuò)散波與馬斯京根法在柵格匯流演算中的應(yīng)用比較[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版, 2013,41(1):6-10.(YAO Cheng,ZHANG Yuxia,LI Zhijia.Application and comparison of cell-to-cell diffusion wave and Muskingum routing methods[J].Journal of Hohai University:Nature Sciences,2013,41(1):6-10.(in Chinese))

        [10]芮孝芳,凌哲,劉寧寧,等.新安江模型的起源及對(duì)其進(jìn)一步發(fā)展的建議[J].水利水電科技進(jìn)展,2012,32 (4):1-5.(RUI Xiaofang,LIN Zhe,LIU Ningning,et al. Origin of Xinanjiang model and its further development [J].Advances in Science and Technology of Water Resources,2012,32(4):1-5.(in Chinese))[11]姚成,紀(jì)益秋,李致家,等.柵格型新安江模型的參數(shù)估計(jì)及應(yīng)用[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,2012,40 (1):42-47.(YAO Cheng,JI Yiqiu,LI Zhijia,et al. Parameter estimation and application of grid-based Xin' anjiang model[J].Journal of Hohai University:Natural Sciences,2012,40(1):42-47.(in Chinese))

        [12]YAO Cheng,LI Zhijia,YU Zhongbo,et al.A priori parameterestimatesforadistributed,grid-based Xinanjiang model using geographically based information [J].Journal of Hydrology,2012,468/469:47-62.

        [13]ANDERSON R,KOREN V,REED S.Using SSURGO data to improve Sacramento Model a priori parameter estimates [J].Journal of Hydrology,2006,320:103-106.

        [14]趙人俊.流域水文模擬[M].北京:水利電力出版社, 1984.

        [15]趙人俊,王佩蘭.新安江模型參數(shù)的分析[J].水文, 1988(6):2-9.(ZHAO Renjun,WANG Peilan.Analysis of the parameters of the Xinanjiang model[J].Journal of China Hydrology,1988(6):2-9.(in Chinese))

        [16]KOREN V,SMITH M,WANG D,et al.Use of soil property data in the derivation of conceptual rainfall-runoff model parameters[C]//Proceedings of the 15th Conference on Hydrology.Long Beach:AMS,2000:103-106.

        [17]中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心國(guó)際科學(xué)數(shù)據(jù)服務(wù)平臺(tái).SRTM 90米分辨率數(shù)字高程原始數(shù)據(jù)[EB/OL]. [2012-09-28].http://datamirror.csdb.cn.

        Effects of DEM resolution on distributed hydrologic simulation

        //YAO Cheng,LI Zhijia,ZHANG Yuxia(College of Hydrology and Water Resources,Hohai University,Nanjing210098,China)

        Based on the digital elevation model(DEM),the Grid-Xin'anjiang model was employed to simulate flood events at different DEM resolutions to assess the effects of resolution on distributed hydrologic simulations.The results indicate that the simulated results may not be improved by using higher resolution DEM data while the accuracy of the watershed delineation can be enhanced.The Grid-Xin'anjiang model is able to produce comparable and promising results for different DEM resolutions with a recalibration of the routing parameters.For medium-sized watersheds,improved computational efficiency and simulation accuracy can be achieved by using the 30 arc-second DEM.

        digital elevation model;grid-Xin'anjiang model;hydrologic simulation;resolution;parameter estimation

        10.3880/j.issn.10067647.2013.05.003

        P333.2

        A

        10067647(2013)05001104

        20121128 編輯:熊水斌)

        國(guó)家自然科學(xué)基金(41101017,41130639,51179045);公益性行業(yè)(氣象)科研專項(xiàng)(GYHY201006037)

        姚成(1982—),男,江蘇揚(yáng)州人,講師,博士,主要從事水文預(yù)報(bào)研究。E-mail:yaocheng@hhu.edu.cn

        猜你喜歡
        模型
        一半模型
        一種去中心化的域名服務(wù)本地化模型
        適用于BDS-3 PPP的隨機(jī)模型
        提煉模型 突破難點(diǎn)
        函數(shù)模型及應(yīng)用
        p150Glued在帕金森病模型中的表達(dá)及分布
        函數(shù)模型及應(yīng)用
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
        3D打印中的模型分割與打包
        好吊妞视频这里有精品| 久久aⅴ无码av高潮AV喷| 亚洲国产成人精品一区刚刚| 99国产精品久久一区二区三区| 热久久美女精品天天吊色| 伊伊人成亚洲综合人网7777| 久久久久无码精品国| 中文字幕一区二区三区乱码人妻| 国语自产偷拍在线观看| 国产一区二区不卡老阿姨 | 久久99精品国产麻豆| 久久久国产精品免费无卡顿| 日本一区人妻蜜桃臀中文字幕 | 亚洲av无码成人网站在线观看| 蜜臀av免费一区二区三区| 国产美女三级视频网站| 在线国人免费视频播放| 久久久久久国产精品免费免费| 97se在线观看| 国产特黄1区2区3区4区| 蜜桃视频在线免费观看| 怡红院免费的全部视频| 波多野结衣一区二区三区视频| 青青草免费视频一区二区| 亚洲一区二区三区四区五区黄| 一本色道av久久精品+网站 | 偷拍夫妻视频一区二区| 丰满人妻被黑人猛烈进入| 久久一区二区三区四区| 在线视频自拍视频激情| 人人妻人人澡人人爽欧美一区双| 久久国产色av| 日本一区二区在线播放观看| 亚洲免费国产中文字幕久久久| 亚洲第一av导航av尤物| 一区二区久久不射av| 国产亚洲综合另类色专区| 成人毛片av免费| 国产精品亚洲A∨天堂不卡| 在线小黄片视频免费播放| 国产尤物精品视频|