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

        ?

        WASMOD水文模擬殘差統(tǒng)計特征檢驗

        2013-06-07 07:17:08李占玲徐宗學(xué)
        水利水電科技進(jìn)展 2013年1期
        關(guān)鍵詞:參數(shù)估計水文方差

        李占玲,徐宗學(xué),周 訓(xùn)

        (1.中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院,北京 100083;2.北京師范大學(xué)水科學(xué)研究院,北京 100875)

        WASMOD水文模擬殘差統(tǒng)計特征檢驗

        李占玲1,徐宗學(xué)2,周 訓(xùn)1

        (1.中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院,北京 100083;2.北京師范大學(xué)水科學(xué)研究院,北京 100875)

        針對采用普通最小二乘法對水文模型進(jìn)行參數(shù)估計時,模型殘差需滿足一定的內(nèi)在統(tǒng)計假定的要求,分別選用Levene檢驗、Kolmogorov-Smirnov檢驗和殘差自相關(guān)系數(shù)圖等方法對水量平衡模型WASMOD的模型殘差的同方差性、正態(tài)性以及相互獨立性等統(tǒng)計假定進(jìn)行了檢驗。結(jié)果表明:當(dāng)對原始數(shù)據(jù)未經(jīng)變換而直接采用普通最小二乘法進(jìn)行參數(shù)估計時,得到的模型殘差滿足相互獨立假定,但并不滿足同方差性和正態(tài)性假定;對原始流量數(shù)據(jù)進(jìn)行開根號變換,可以很好地解決模型殘差的異方差性和非正態(tài)性問題。在模型殘差統(tǒng)計假定得以滿足的條件下,在月時間尺度上WASMOD模型可以為鶯落峽流域的流量模擬與預(yù)報提供良好的工具。

        WASMOD;參數(shù)估計;普通最小二乘法;模型殘差;水量平衡模型;鶯落峽

        由于對時間和空間尺度的依賴,無論是何種類型的水文模型都有部分參數(shù)無法通過試驗直接確定,而必須借助參數(shù)估計方法來確定[1-2]。在進(jìn)行水文模型參數(shù)估計時,其中一種假定是假設(shè)最優(yōu)參數(shù)存在且唯一[2],此時模型參數(shù)估計問題就轉(zhuǎn)化成一個最優(yōu)化問題,通常采用優(yōu)化方法求解,而優(yōu)化方法大多都包含一些統(tǒng)計假定在內(nèi)[3-4]。普通最小二乘法(ordinary least square,OLS)是常用的一種水文模型參數(shù)估計方法,采用該方法進(jìn)行參數(shù)估計時,需假定模型殘差序列相互獨立且服從均值為零、方差穩(wěn)定的正態(tài)分布[4]。但很少有研究對這些內(nèi)在的統(tǒng)計假定進(jìn)行檢驗,且已有研究表明這些假定多數(shù)情況下都得不到滿足,這將意味著由此得到的最優(yōu)參數(shù)并不一定最優(yōu)[4-6];另外,一個模型對流域的降雨徑流關(guān)系描述的成功與否通常由模型殘差來表征,這進(jìn)一步說明了在建模及模型應(yīng)用過程中對殘差的有關(guān)統(tǒng)計假定進(jìn)行檢驗的必要性。本文以水量平衡模型WASMOD[7](thewater andsnowbalance modeling system,以下簡稱WASMOD模型)為例,選用OLS對模型進(jìn)行參數(shù)估計,在此基礎(chǔ)上探討模型殘差的一些統(tǒng)計特性。

        1 WASMOD模型

        WASMOD模型是由Xu等[7]開發(fā)的集總式概念性水量平衡模型,目前該模型已經(jīng)在不同氣候和土壤條件下的20多個國家的200多個流域進(jìn)行了應(yīng)用和驗證[7],在我國黑河流域[8-9]以及南部東江流域[6]也得到了很好的應(yīng)用。該模型既可以用于流量序列的外延和流域水量平衡計算,也可以用來模擬和預(yù)測氣候變化條件下的水文響應(yīng)及周、月尺度流量預(yù)測等。

        WASMOD模型結(jié)構(gòu)如圖1所示。根據(jù)氣溫的不同,降水分為降雨和降雪。降雨作為有效降水進(jìn)入到土壤含水層之前,有一部分被蒸發(fā)掉。降雪被加入到積雪層,其中有一部分降雪融化,進(jìn)入土壤含水層。土壤含水層中的水分一部分蒸發(fā),一部分形成地表徑流,還有一部分形成基流。以月時間尺度為例,WASMOD模型輸入數(shù)據(jù)包括研究區(qū)的月降水、月潛在蒸散發(fā)和月平均氣溫;輸出包括徑流和其他水量平衡要素,如實際蒸散發(fā)、地表徑流、基流等。模型主要包括積雪、融雪、實際蒸散發(fā)、基流、地表徑流、總徑流和流域水量平衡等計算模塊。有關(guān)該模型原理、結(jié)構(gòu)及應(yīng)用詳見文獻(xiàn)[7]。

        圖1 WASMOD模型結(jié)構(gòu)

        2 研究區(qū)概況及數(shù)據(jù)資料

        以黑河上游鶯落峽流域為研究區(qū)。該流域面積1.0萬km2,河道長303 km,河道兩岸山高谷深,河床陡峻,氣候陰濕寒冷,植被較好,多年平均氣溫不足2℃,年降水量為350 mm,是黑河流域的產(chǎn)流區(qū)。流域內(nèi)及周邊有6個雨量站(野牛溝、扎木什克、祁連、托勒、張掖、鶯落峽)和4個氣象站(野牛溝、祁連、托勒、張掖),站點位置如圖2所示。選擇流域控制站鶯落峽站的流量數(shù)據(jù)率定和驗證WASMOD模型,數(shù)據(jù)選用時段為1987年1月—2000年12月,其中,1987—1989年的數(shù)據(jù)作為預(yù)熱期,1990—1996年的作為率定期,1997—2000年的作為驗證期。數(shù)據(jù)取自中國西部環(huán)境與生態(tài)科學(xué)數(shù)據(jù)中心(http://westdc.westgis.ac.cn/)和數(shù)字黑河(http:// heihe.westgis.ac.cn/)。

        圖2 黑河流域及其上游鶯落峽位置示意圖

        為準(zhǔn)備WASMOD模型的輸入數(shù)據(jù),根據(jù)泰森多邊形法,由6個雨量站月降水?dāng)?shù)據(jù)計算得到流域面月降水?dāng)?shù)據(jù);由4個氣象站的月氣溫數(shù)據(jù)計算得到整個流域氣溫數(shù)據(jù);采用Penman-Monteith方法計算流域的潛在蒸散發(fā)[10-11]。將以上計算得到的流域面降水量、氣溫和潛在蒸散發(fā)數(shù)據(jù)作為WASMOD模型的輸入數(shù)據(jù)。

        3 模型參數(shù)估計及殘差統(tǒng)計特征分析方法

        采用OLS進(jìn)行參數(shù)估計時,模型殘差et需滿足以下條件[4,12]:①無偏差;②方差穩(wěn)定;③取值在時間上相互無關(guān);④服從正態(tài)分布。當(dāng)模型殘差存在異方差時,參數(shù)的OLS估計量不再是具有最小方差的估計量;參數(shù)估計量的方差是有偏的,參數(shù)的假設(shè)檢驗也并非最有效,這會降低模型估計和預(yù)測的精度。當(dāng)殘差項存在自相關(guān)時,會夸大估計參數(shù)的顯著性。因此,有必要對模型殘差的同方差性、相互獨立性和正態(tài)性等統(tǒng)計假定進(jìn)行檢驗。

        3.1 同方差性檢驗

        采用Levene檢驗方法檢驗殘差序列的同方差性[13-14]。計算Levene檢驗統(tǒng)計量F:

        式中:P為樣本組數(shù);mj為第j(j=1,2,…,P)組樣本數(shù);xij為第j組第i個原始數(shù)據(jù)經(jīng)數(shù)據(jù)轉(zhuǎn)換后的新的變量值。Levene檢驗統(tǒng)計量F服從自由度為P-1 和N-P-1(N為各組樣本數(shù)之和)的F分布。當(dāng)F≥F(α,P-1,N-P+1)時,則在α=0.05或α=0.01顯著性水平下拒絕H0(各組樣本之間具有同方差性),接受H1,可認(rèn)為各樣本方差不全相等;當(dāng)F<F(α,P-1,N-P+1)時,則不拒絕H0,可以認(rèn)為各樣本方差齊性。

        3.2 正態(tài)分布檢驗

        Kolmogorov-Smirnov(K-S)檢驗方法是可以用來檢驗單一樣本是否來自某一特定分布的非參數(shù)檢驗方法[15-16]。設(shè)理論累積頻率分布為F0(x)(這里的理論分布是正態(tài)分布),一組隨機(jī)樣本的經(jīng)驗累積頻率分布為FN(x),FN(x)=k/N,其中k為樣本出現(xiàn)的累積次數(shù)。把樣本觀測值從小到大排列,計算理論累積分布函數(shù)F0(x)和經(jīng)驗累積分布函數(shù)FN(x)。設(shè)D 為F0(x)和差距的最大值:

        當(dāng)實際觀測D>D(N,α)(D(N,α)為顯著性水平為α、樣本數(shù)為N時,D的拒絕臨界值),則拒絕H0(樣本來自的總體分布服從某特定分布);反之則接受H0假設(shè)。

        3.3 相互獨立性檢驗

        殘差的相互獨立性可以通過殘差自相關(guān)系數(shù)圖來判斷[4,8]。殘差序列et的k階自相關(guān)系數(shù)ρk=E[(et-ˉe)(et+k-ˉe)]/σ2,其中,ˉe和σ2為殘差項的均值和方差。ρk的估計量^ρk為

        當(dāng)N很大而k很小時,N/(N-k)→1,則k階自相關(guān)系數(shù)ρk的估計量為

        在α=0.05顯著性水平下,一個相互獨立序列的自相關(guān)系數(shù)的置信區(qū)間為

        4 結(jié)果分析與討論

        4.1 殘差序列et的統(tǒng)計假定檢驗

        定義模型殘差et=Yobs,t-^Yt,其中Yobs,t為流量觀測值,為流量模擬值,t=1,2,…,N。模型率定時目標(biāo)函數(shù)為殘差平方和最小?;诖四繕?biāo)函數(shù),求出最優(yōu)參數(shù)及對應(yīng)的殘差序列,并對此殘差序列的同方差性、相互獨立性和正態(tài)性進(jìn)行檢驗。

        首先將殘差序列et分成3組:第1組,Yobs,t<0.9對應(yīng)的殘差序列,序列長度m1;第2組,0.9≤Yobs,t≤1.1對應(yīng)的殘差序列,序列長度m2;第3 組,Yobs,t>1.1對應(yīng)的殘差序列,序列長度m3。通過計算可得令,其中eij為原始?xì)埐顢?shù)據(jù),為原始?xì)埐顢?shù)據(jù)中第j組(j=1,2,3)樣本的算術(shù)平均。由原始?xì)埐钚蛄猩尚碌男蛄衶xij}并計算Levenes檢驗統(tǒng)計量F。計算得到F=11.18,查F分布表得到其臨界值為F(0.01,2,60)=4.98,F(0.01,2, 120)=4.79。由于F=11.18>F(0.01,2,82),因此在α=0.01顯著性水平下,拒絕H0,認(rèn)為殘差序列et不具有同方差性。

        圖3(a)為殘差序列et的經(jīng)驗累積頻率曲線和理論累積頻率曲線以及二者差距的最大值D值,可以看出,兩條累積頻率曲線并不十分接近,計算得到的D=0.184>D(84,0.05)=0.148,所以拒絕H0假設(shè),即殘差序列et來自的總體分布不服從正態(tài)分布。

        圖4(a)為殘差序列et的自相關(guān)系數(shù)圖,可以看出,,即在α=0.05顯著性水平下,模型殘差序列et滿足相互獨立假設(shè)。

        上述研究表明,如果定義模型殘差為et=Yobs,t-,模型率定時目標(biāo)函數(shù)為殘差平方和最小,則得到的模型殘差序列et滿足相互獨立但不滿足同方差性和正態(tài)性假定。

        4.2 殘差序列ut的統(tǒng)計假定檢驗

        對于方差異性的模型殘差,常常需要對原始流量數(shù)據(jù)進(jìn)行某種變換,如開根號變換、Box-Cox變換、取對數(shù)變換等[4-5,8]。Xu[4]研究發(fā)現(xiàn),對模型驗證數(shù)據(jù)做開根號變換可以很好地解決WASMOD模型殘差方

        圖3 殘差序列et和ut的經(jīng)驗累積頻率曲線和理論累積頻率曲線

        圖4 殘差序列et和ut自相關(guān)系數(shù)

        同樣采用Levene檢驗方法對殘差ut的同方差性進(jìn)行檢驗,方法及步驟同前。計算得到對應(yīng)新的殘差序列ut統(tǒng)計量F=3.48,由于F<F(0.01,2,82),因此,可以認(rèn)為在α=0.01顯著性水平下,新的殘差序列ut具有同方差性。

        采用K-S方法檢驗殘差序列ut來自的總體分布是否服從正態(tài)分布。圖3(b)為殘差序列ut的經(jīng)驗累積頻率曲線和理論累積頻率曲線以及二者差距的最大值D值,可以看出,兩條累積頻率曲線很接近,計算得到的D=0.082<D(84,0.05)=0.148,接受H0假設(shè),即殘差序列ut來自的總體分布服從正態(tài)分布。

        通過殘差自相關(guān)系數(shù)圖來判斷新構(gòu)建的殘差序列是否相互獨立。圖4(b)為殘差序列ut的自相關(guān)系數(shù)圖,可以看出,,即在0.05顯著性水平下,模型殘差序列ut也相互獨立。

        以上檢驗說明,變換后的殘差序列滿足同方差性、正態(tài)性和相互獨立性等統(tǒng)計假定,模型參數(shù)的OLS估計量是具有最小方差的估計量,此時該模型對本研究區(qū)的流量具有良好的模擬和預(yù)測能力。

        4.3 流量模擬結(jié)果

        基于變換后的目標(biāo)函數(shù),得到模型參數(shù)的最優(yōu)值以及模型對流量模擬的最優(yōu)值。圖5為流域鶯落峽站實測與模擬月徑流過程的對比以及流域月降水量分布,可以看出,鶯落峽站的月模擬流量過程線與實測流量過程線在趨勢上吻合較好。在率定和驗證過程中,流量模擬的Nash-Sutcliffe有效性系數(shù)分別達(dá)到0.942和0.928,流量模擬相對誤差[8]分別為-0.4%和-4.0%,在±5%以內(nèi),說明模型模擬的總體水量平衡效果較好。其中,率定期1996年、1990年、1991等年份月流量過程模擬的最好,ENS分別達(dá)到0.984、0.973和0.958,相對誤差分別為-2.50%、-4.10%和3.10%;驗證期中1998年效果最好,ENS達(dá)到0.943,流量模擬相對誤差為-3.80%。比較未做變換時對應(yīng)的流量序列模擬值,在率定期和驗證期ENS分別為0.938和0.920,說明對流量數(shù)據(jù)進(jìn)行變換后模型的模擬效果并沒有降低??傮w而言,WASMOD模型在黑河鶯落峽流域的徑流模擬方面具有良好的適用性,在月時間尺度上可以為研究區(qū)域的徑流模擬與預(yù)報提供良好的工具。

        圖5 WASMOD模型對鶯落峽站月流量過程模擬結(jié)果

        5 結(jié) 語

        OLS是較為常用的水文模型參數(shù)估計方法之一,然而以往的研究中在采用OLS進(jìn)行參數(shù)估計時,往往缺乏或忽略了對模型殘差需滿足的一些內(nèi)在統(tǒng)計假定進(jìn)行檢驗。如果模型殘差的這些統(tǒng)計假定不能滿足,則由此得到的最優(yōu)參數(shù)并不一定最優(yōu),導(dǎo)致模型的模擬和預(yù)測出現(xiàn)偏差。本文采用OLS對水量平衡模型WASMOD進(jìn)行了參數(shù)估計,并選用Levene檢驗、Kolmogorov-Smirnov檢驗和殘差自相關(guān)系數(shù)圖等方法對模型殘差的同方差性、正態(tài)性以及相互獨立性等統(tǒng)計假定進(jìn)行了檢驗。結(jié)果表明,當(dāng)對原始數(shù)據(jù)未經(jīng)任何變換而直接采用OLS進(jìn)行參數(shù)估計時,得到的模型殘差滿足相互獨立假定,但并不滿足同方差性和正態(tài)性假定;當(dāng)把原始數(shù)據(jù)進(jìn)行開根號變換時,得到的殘差序列滿足同方差性和正態(tài)性假定。在此統(tǒng)計假定得以滿足的條件下,對黑河上游鶯落峽流域的流量進(jìn)行了模擬,結(jié)果表明,在月時間尺度上WASMOD模型可以很好地模擬鶯落峽流域的流量過程。

        [1]REFSGAARD J C,SETH S M,BATHURST J C,et al. Application of the SHE to catchments in India:partⅠ: general results[J].Journal of Hydrology,1992,140:1-23.

        [2]王書功.水文模型參數(shù)估計方法及參數(shù)估計不確定性研究[D].北京:中國科學(xué)院研究生院,2006.

        [3]WANG Y,DIETRICH J,VOSS F,et al.Identifying and reducing model structure uncertainty based on analysis of parameter interaction[J].Advances in Geosciences,2007, 11:117-122.

        [4]XU Chongyu.Statistical analysis of parameters and residuals of a conceptual water balance model-methodology and case study[J].Water Resources Management,2001,15:75-92.

        [5]YANG Jing,REICHERT P,ABBASPOUR K C,et al. Hydrological modelling of the Chaohe Basin in China: statistical model formulation and Bayesian inference[J]. Journal of Hydrology,2007,340:167-182.

        [6]LI Lu,XU Chongyu,XIA Jun,et al.Uncertainty estimates by Bayesian method with likelihood of AR(1)plus normal model and AR(1)plus multi-normal model in different time-scales hydrological models[J].Journal of Hydrology, 2011,406:54-65.

        [7]XU Chongyu.WASMOD:the water and snow balance modelling system[M]//SINGH V P,FREVERT D K. Mathematical Models of Small Watershed Hydrology and Applications:Chapter 17.Chelsea,Michigan,USA:Water Resources Publications,LLC,2002:555-590.

        [8]LI Zhanling,XU Zongxue,LI Zhanjie.Performance of WASMOD andSWATonhydrologicalsimulationin YingluoxiawatershedinnorthwestofChina[J]. Hydrological Processes,2011,25:2001-2008.

        [9]李占玲,徐宗學(xué).黑河流域上游山區(qū)徑流模擬及模型評估[J].北京師范大學(xué)學(xué)報:自然科學(xué)版,2010,46(3): 344-349.(LIZhanling,XUZongxue.Assessmenton hydrological models for runoff simulation in the upper reaches of the Heihe River Basin[J].Jonrnal of Beijing Normal University:Natural Science,2010,46(3):344-349.(in Chinese))

        [10]MONTEITH JL.Evaporation,theenvironment[C]// XIXth symposia of the society for Experimental Biology in the State and Movement of Water in Living Organisms. Cambridge:Cambridge University Press,1965:205-234.

        [11]左德鵬,徐宗學(xué),程磊,等.渭河流域潛在蒸散量時空變化及其突變特征[J].資源科學(xué),2011,33(5):975-982. (ZUO Depeng,XU Zongxue,CHENG Lei,et al.Spatialtemporalvariationsandmutationsofpotential evapotranspirationintheWeiheRiverBasin[J]. Resources Science,2011,33(5):975-982.(in Chinese))

        [12]單明,聶燕萍.最小二乘法直線擬合基本假定的幾點討論[J].大學(xué)物理實驗,2008,21(4):63-65.(SHAN Ming,NIE Yanping.Discussion on assumptions of least square linear fitting[J].Physical Experiment of College, 2008,21(4):63-65.(in Chinese))

        [13]龔秀芳,馮珍珍.幾種異方差檢驗方法的比較[J].菏澤師范??茖W(xué)校學(xué)報,2003,25(4):19-23.(GONG Xiufang,FENGZhenzhen.Comparingseveraltesting methodsofheteroscedasticity[J].JournalofHeze Teachers College,2003,25(4):19-23.(in Chinese))

        [14]郭顯光.方差齊次性檢驗的方法[J].統(tǒng)計教育,2005 (4):32-33.(GUOXianguang.Testingmethodsfor homoscedasticity[J].Statistical Education,2005(4):32-33.(in Chinese))

        [15]褚健婷,夏軍,許崇育,等.海河流域氣象和水文降水資料對比分析及時空變異[J].地理學(xué)報,2009,64(9): 1083-1092.(CHU Jianting,XIA Jun,XU Chongyu,et al. Comparisonandspatial-temporalvariabilityofdaily precipitation data of weather stations and rain gauges in Haihe River Basin[J].Acta Geographica Sinica,2009,64 (9):1083-1092.(in Chinese))

        [16]孫鵬,張強(qiáng),陳曉宏.鄱陽湖流域枯水徑流演變特征、成因與影響[J].地理研究,2011,30(9):1702-1712. (SUN Peng,ZHANG Qiang,CHEN Xiaohong.Changing propertiesoflowstreamflow:possiblecausesand implications[J].Geographical Research,2011,30(9): 1702-1712.(in Chinese))

        [17]NASH J E,SUTCLIFFE J V.River flow forecasting through conceptual models:partⅠ:a discussion of principles[J].Journal of Hydrology,1970,10,282-290.

        [18]張利平,曾思棟,王任超,等.氣候變化對灤河流域水文循環(huán)的影響及模擬[J].資源科學(xué),2011,33(5):966-974.(ZHANG Liping,ZENG Sidong,WANG Renchao,et al.Impacts of climate change on the hydrological cycle in the Luan River Basin[J].Resources Science,2011,33 (5):966-974.(in Chinese))

        Check of statistical features of model residuals of WASMOD

        //LI Zhanling1,XU Zongxue2,ZHOU Xun1(1.College of Water Resources and Environment,China University of Geosciences,Beijing100083;China;2.College of Water Sciences,Beijing Normal University,Beijing100875,China)

        When the parameters for hydrological models are estimated by using the ordinary least square method,the model residuals should satisfy certain statistical assumptions.The Levene test,Kolmogorov-Smirnov test and autocorrelation coefficient are employed to check the statistical assumptions of model residuals of the conceptual lumped Water And Snow balance MODeling system(WASMOD)such as homoscedasticity,normality and independence features.As the original discharge values without any transformation are calculated by using the ordinary least square method,the yielded model residuals are proved to be independent,which greatly deviate from the features of homoscedasticity and normality. However,a square root transformation for the discharge values is capable of solving the problems of heteroscedasticity and non-normality in model residuals.The proposed WASMOD model may be satisfactory tool for the simulation and predication of discharge in Yingluoxia watershed on monthly temporal scale when the statistical assumptions of the model residuals are satisfied.

        WASMOD;parameter estimation;ordinary least square method;model residual;water balance model; Yingluoxia watershed

        10.3880/j.issn.10067647.2013.01.003

        P333.6

        A

        10067647(2013)01001305

        2012-04-26 編輯:熊水斌)

        國家自然科學(xué)基金(41101038);中央高?;究蒲袠I(yè)務(wù)費專項(2011YXL038)

        李占玲(1980—),女,內(nèi)蒙古赤峰人,講師,博士,主要從事水文及水資源研究。E-mail:zhanling.li@cugb.edu.cn

        猜你喜歡
        參數(shù)估計水文方差
        2022年《中國水文年報》發(fā)布
        方差怎么算
        基于新型DFrFT的LFM信號參數(shù)估計算法
        概率與統(tǒng)計(2)——離散型隨機(jī)變量的期望與方差
        水文
        水文水資源管理
        計算方差用哪個公式
        方差生活秀
        水文
        Logistic回歸模型的幾乎無偏兩參數(shù)估計
        国产91成人自拍视频| 国产精品视频一区日韩丝袜| 日韩在线看片| 国产在线不卡免费播放| 热re99久久精品国产66热6| 少妇极品熟妇人妻高清| 亚洲一区二区岛国高清| 亚洲乱码中文字幕一线区| 亚洲综合久久久中文字幕| 亚洲av毛片在线播放| 91色区在线免费观看国产| 国产精品蝌蚪九色av综合网| 国产不卡视频一区二区三区| 色婷婷综合久久久久中文字幕| 玩弄放荡人妻少妇系列| 97人人超碰国产精品最新o| 天天综合久久| 国产精品一区二区三区黄片视频| 日本一区二区三区综合视频| 日本区一区二区三视频 | 亚洲精品久久久久久动漫| 中文字幕天堂网| 国产精品情侣露脸av在线播放| 扒下语文老师的丝袜美腿| 国产精品一区二区三区三| 人妻熟妇乱又伦精品hd| 免费毛片a线观看| 亚洲аv天堂无码| 亚洲va成无码人在线观看| 91久久精品一二三区色| 久久精品一区二区三区蜜桃| 国产老熟妇精品观看| 免费人妻无码不卡中文字幕18禁| 亚洲国产人在线播放首页| 无码免费人妻超级碰碰碰碰| 国产一区二区三区免费精品| 成人在线观看av毛片| 国产又a又黄又潮娇喘视频| 九九热在线视频观看这里只有精品| 永久免费看免费无码视频| 亚洲日产乱码在线中文字幕|