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

產(chǎn)流誤差動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法的應(yīng)用

2016-08-05 07:41:14束慧連包為民楊?yuàn)檴?/span>劉可新趙麗平楊小強(qiáng)
水力發(fā)電 2016年4期

束慧連,包為民,2,楊?yuàn)檴?劉可新,趙麗平,楊小強(qiáng)

(1.河海大學(xué)水文水資源學(xué)院,江蘇南京210098;2.河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,江蘇南京210098)

?

產(chǎn)流誤差動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法的應(yīng)用

束慧連1,包為民1,2,楊?yuàn)檴?,劉可新1,趙麗平1,楊小強(qiáng)1

(1.河海大學(xué)水文水資源學(xué)院,江蘇南京210098;2.河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,江蘇南京210098)

摘要:動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法(Dynamic System Response Curve,DSRC)是一種基于微分響應(yīng)的有明確物理意義并且有效的實(shí)時(shí)修正新方法,屬于向誤差源頭進(jìn)行修正的方法。通過計(jì)算產(chǎn)流的系統(tǒng)響應(yīng)曲線對(duì)產(chǎn)流系列進(jìn)行修正,實(shí)現(xiàn)流量過程的修正。將該方法與新安江模型相結(jié)合應(yīng)用于閩江上游的水吉流域,并將修正效果與遞推最小二乘法進(jìn)行了比較。利用水吉流域29場(chǎng)洪水資料進(jìn)行模擬和修正,結(jié)果表明該方法可以有效提高洪水預(yù)報(bào)的精度,且其效果優(yōu)于遞推最小二乘法。

關(guān)鍵詞:產(chǎn)流誤差;動(dòng)態(tài)系統(tǒng)響應(yīng)曲線;水吉流域

0引言

實(shí)時(shí)洪水預(yù)報(bào)作為防洪減災(zāi)的非工程措施,有其重要的經(jīng)濟(jì)和民生意義。然而,水文模型預(yù)報(bào)洪水習(xí)慣上是采用觀測(cè)的歷史水文資料來確定模型參數(shù),然后用于未來的洪水預(yù)報(bào)中。但這樣的預(yù)報(bào)方案往往在實(shí)時(shí)洪水預(yù)報(bào)中得不到滿意的結(jié)果。因此,為了提高實(shí)時(shí)洪水預(yù)報(bào)的精度,有必要進(jìn)行實(shí)時(shí)誤差修正。所謂實(shí)時(shí)洪水預(yù)報(bào)誤差修正就是對(duì)在水文模型中沒有考慮或是考慮不得當(dāng)卻對(duì)預(yù)報(bào)結(jié)果有影響的因素,利用實(shí)測(cè)信息和其他可利用信息進(jìn)行實(shí)時(shí)校正,達(dá)到彌補(bǔ)水文模型不足的目的[1]。按修正

內(nèi)容劃分,實(shí)時(shí)修正方法可以分為:模型誤差修正,如自回歸模型[2];模型參數(shù)修正,如神經(jīng)網(wǎng)絡(luò)算法[3];模型輸入修正,如抗差估計(jì)理論[4];模型狀態(tài)修正,如卡爾曼濾波方法[5- 6];綜合修正,就是前四種方法的結(jié)合[7]。

現(xiàn)有的誤差修正方法種類很多[8-12],源于海洋與氣象學(xué)的數(shù)據(jù)同化方法也被應(yīng)用于水文預(yù)報(bào)的誤差修正中[13]。然而,常見的誤差修正方法存在各種各樣的問題,如沒有明確物理意義、預(yù)見期的損失、資料不能合理利用等,再加上很多修正方法存在假設(shè)條件,當(dāng)實(shí)際情況不能滿足假設(shè)條件的時(shí)候修正效果顯然不會(huì)理想。例如:卡爾曼濾波方法雖然在修正過程中明確地考慮了模型和資料的不確定性,但是假設(shè)條件不易得到滿足,應(yīng)用于分布式水文模型、多參數(shù)模型等時(shí)計(jì)算量大,結(jié)構(gòu)復(fù)雜,實(shí)時(shí)洪水預(yù)報(bào)難以提供誤差修正所需的信息量;自回歸模型假設(shè)誤差序列按時(shí)間排序是相關(guān)的,但是當(dāng)預(yù)報(bào)量發(fā)生突變時(shí)破壞了相關(guān)關(guān)系,因此不再適用,修正能力較弱,而且隨著預(yù)見期的增加修正效果迅速下降;神經(jīng)網(wǎng)絡(luò)算法雖然能模擬變量間的非線性關(guān)系,但其結(jié)構(gòu)不唯一,且收斂極為緩慢。

動(dòng)態(tài)響應(yīng)曲線修正方法(Dynamic System Response Curve,DSRC)是一種追溯至誤差源頭的修正方法,沒有采用傳統(tǒng)的誤差修正方法的修正模式——在預(yù)報(bào)結(jié)果中進(jìn)行誤差修正;而是通過修正用于進(jìn)行洪水預(yù)報(bào)的資料達(dá)到誤差修正的目的。新安江模型分為四個(gè)部分,產(chǎn)流誤差動(dòng)態(tài)響應(yīng)曲線修正方法將產(chǎn)流以下的部分作為響應(yīng)系統(tǒng),并用差分近似代替微分,得到時(shí)段產(chǎn)流量所對(duì)應(yīng)的系統(tǒng)響應(yīng)曲線,將實(shí)測(cè)流量和計(jì)算流量的差值作為修正信息,用最小二乘估計(jì)修正產(chǎn)流過程,達(dá)到修正計(jì)算流量過程的目的[14]。本文在閩江上游的水吉流域?qū)⒃撜`差修正方法與新安江模型相結(jié)合,比較分析了修正前后的結(jié)果,驗(yàn)證了該方法的可行性和有效性。

1修正方法介紹

圖1 系統(tǒng)示意

新安江模型將模型結(jié)構(gòu)分散為蒸散發(fā)計(jì)算、產(chǎn)流計(jì)算、分水源計(jì)算和匯流計(jì)算四個(gè)層次,匯流計(jì)算又包含坡地匯流和河道匯流兩部分,坡地匯流計(jì)算采用線性水庫,河道匯流采用馬斯京根分段連續(xù)演算。文獻(xiàn)[14]提出將產(chǎn)流誤差修正的動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法是將新安江模型的分水源計(jì)算和匯流計(jì)算這兩個(gè)層次作為一個(gè)系統(tǒng),如圖1所示,系統(tǒng)響應(yīng)曲線由線性差分近似代替微分得到,修正信息為實(shí)測(cè)流量系列,采用最小二乘估計(jì)原理,對(duì)產(chǎn)流量R進(jìn)行修正,最終實(shí)現(xiàn)計(jì)算流量QC的修正。

1.1動(dòng)態(tài)系統(tǒng)響應(yīng)曲線

由文獻(xiàn)[14]可知,所謂動(dòng)態(tài)系統(tǒng)響應(yīng)曲線,就是給某時(shí)段的產(chǎn)流量ri增加一個(gè)單位的值,用變化后的產(chǎn)流系列帶入模型計(jì)算得到新的流量過程,這個(gè)新的流量過程與改變之前得到的模型計(jì)算流量過程Qc相減得到的曲線就是對(duì)應(yīng)ri的系統(tǒng)響應(yīng)曲線。由于該曲線受到其他時(shí)段產(chǎn)流量的影響實(shí)時(shí)變化,因此稱作動(dòng)態(tài)系統(tǒng)響應(yīng)曲線。文獻(xiàn)[14]對(duì)其計(jì)算過程進(jìn)行了推導(dǎo)。

方法的關(guān)鍵是U矩陣的求解,即

(1)

式中,RC為模型初始計(jì)算徑流系列[rc1,…,rcn];t為時(shí)間;θ為模型參數(shù);L為樣本長(zhǎng)度;Q(RC,θ,t)為模型初始計(jì)算流量過程QC。

U矩陣的每一項(xiàng)都可以通過線性差分來代替微分近似求得

(t

(2)

得到U矩陣的具體計(jì)算步驟為:

(1)首先直接通過新安江模型進(jìn)行計(jì)算,得到初始的計(jì)算徑流過程RC和流量過程QC。

(2)在以計(jì)算得到的徑流系列的基礎(chǔ)上給i時(shí)段的徑流量增加1個(gè)單位值,將變化后的徑流系列帶入新安江模型的分水源計(jì)算和匯流計(jì)算中,求得對(duì)應(yīng)的流量過程。

(3)用步驟(2)中得到的流量過程與初始計(jì)算流量過程相減,得到對(duì)應(yīng)時(shí)段徑流的響應(yīng)曲線,即U矩陣中的第i列。

(4)完成所有列的計(jì)算則得到U矩陣。

計(jì)算得到U矩陣后,產(chǎn)流量的修正量ΔR通過最小二乘法得到,即

ΔR=(UTU)-1UT(Q(R,θ,L)-Q(RC,θ,t))

(3)

將修正后的產(chǎn)流量R′=RC+ΔR帶入新安江模型,計(jì)算得到修正后的計(jì)算流量過程,可以實(shí)現(xiàn)誤差修正。

1.2遞推最小二乘法

描述時(shí)間系列的線性回歸模型[7],其觀測(cè)系列為(x11,x12,…,x1p;y1),(x21,x22,…,x2p;y2),…,(xt1,xt2,…,xtp;yt)。其中,x為自變量;y為因變量;t為時(shí)刻;p為模型回歸階數(shù)。若已知t時(shí)刻之前的自變量和因變量的觀測(cè)值,要估計(jì)t+1時(shí)刻的因變量,可根據(jù)這些觀測(cè)信息先用最小二乘法估計(jì)出參數(shù),然后再估計(jì)y值。遞推最小二乘法可以解決當(dāng)中計(jì)算量的問題和隨著信息量的增大計(jì)算機(jī)容量的問題。

將每個(gè)觀測(cè)值帶入回歸方程

(4)

寫成矩陣形式為

Yt=XtCt+Et

(5)

表1水吉流域使用新安江模型的參數(shù)

WMUMLMKCBIMCSMEXKGKICGCICSKEXE16320801.350.7900.1613.751.50.450.40.9940.920.7451.99770.2175

注:表中WM為流域平均張力水容量,mm;UM為上土層張力水容量,mm;LM為下土層張力水容量,mm;KC為蒸散發(fā)折算系數(shù);B為蓄水容量面積分布曲線指數(shù);IM為不透水面積占全流域面積比例;C為深層蒸散發(fā)擴(kuò)散系數(shù);SM為自由水蓄水容量,mm;EX為自由水蓄水容量面積分布曲線指數(shù);KG和KI分別為自由水蓄水庫對(duì)地下水和壤中流的日出流系數(shù);CG、CI、CS分別為地下水、壤中流、地面徑流;KE和XE為馬斯京根法參數(shù),其中KE單位為h。

在t-1時(shí)刻和t時(shí)刻可得參數(shù)的最小二乘估計(jì)

(6)

(7)

經(jīng)過一系列的簡(jiǎn)化和計(jì)算可以得到參數(shù)的遞推最小二乘估計(jì)

(8)

2應(yīng)用與檢驗(yàn)

2.1流域簡(jiǎn)介

閩江流域處于亞熱帶季風(fēng)氣候區(qū),雨量豐沛,氣候溫和。流域多年平均降水量1 724 mm,年內(nèi)分布不均,年際變化大,其中60%左右的年降水量主要集中在4月~6月,豐水年降雨量是枯水年的2~3倍。受氣候的影響,流域內(nèi)多年平均蒸發(fā)量為915.0 mm[15],夏季蒸發(fā)量較大,冬季蒸發(fā)量較小。

水吉流域地處閩江上游地區(qū),南平以上閩江支流南浦溪流經(jīng)該流域,流域面積3 767 km2,多年平均氣溫15~20℃。本文選用水吉流域1988年~1999年之間的資料進(jìn)行分析,流域年平均降水量為1 820 mm,年平均徑流深972.0 mm[16],徑流系數(shù)大于0.4,屬于典型的濕潤(rùn)地區(qū),使用三水源新安江模型進(jìn)行洪水模擬預(yù)報(bào)是較為合理的。

水吉流域境內(nèi)有11個(gè)雨量站,分別為:游楓、東坑、浦城、仙陽、高坊、山下、洋源、舊館、管坦、水吉、外屯,有一個(gè)蒸發(fā)站,以水吉站作為預(yù)報(bào)站。

2.2應(yīng)用檢驗(yàn)

本文選用水吉流域1988年~1999年間29場(chǎng)洪水資料,使用新安江模型進(jìn)行模擬,模型使用的參數(shù)采用日模參數(shù)率定獲得的結(jié)果,具體見表1。

驗(yàn)證產(chǎn)流誤差的動(dòng)態(tài)響應(yīng)曲線修正方法的修正效果,具體通過以下3項(xiàng)指標(biāo)來衡量修正效果:

(1)確定性系數(shù)DC

(9)

(2)洪峰相對(duì)誤差ΔQP

ΔQP=((QCP-QOP)/QOP)×100%

(10)

(3)洪量相對(duì)誤差ΔQT

ΔQT=((RC-RO)/RO)×100%

(11)

29場(chǎng)洪水修正前和修正后的結(jié)果比較及修正效果見表2、表3。43880620和43960317場(chǎng)次洪水的修正結(jié)果如圖2和圖3所示。

表2水吉流域29場(chǎng)洪水實(shí)時(shí)修正前后計(jì)算結(jié)果

洪號(hào)實(shí)測(cè)值ROP/mmQOP/m3·s-1新安江模型ΔQP/%ΔRT/%DC動(dòng)態(tài)系統(tǒng)響應(yīng)曲線法ΔQP/%ΔRT/%DC遞推最小二乘法ΔQP/%ΔRT/%DC4388022876.61810-15.91-4.230.914-12.451.50.974-5.841.440.9404388050934.98313.981.180.7231.411.430.8285.02-0.320.7904388062094.82680-32.57-4.270.6787.112.080.9145.65-4.70.8264389061846.52020-26.5-20.710.8382.652.740.993-9.55-6.980.94543890629119.614109.8319.940.761-0.2-1.660.9402.763.480.8894390060123.4636-32.360.450.663-0.10.310.9770.83-1.550.83543900608114.12085-19.63-14.940.864-17.83-1.270.924-18.53-0.260.8694391032651.97442.0829.860.225-1.04-3.310.897-0.951.90.8264391042634.5540.0618.8311.80.5037.12-4.370.65221.8614.890.50943920321130.3101012.8819.70.6717.659.030.75411.99.390.8304392050125.458030.374.610.424-4.18-0.430.9920.79-3.60.92043920616542010-0.32.530.953-0.46-1.80.9916.152.930.93343920704158.434702.168.80.9540.2-1.040.982-6.281.780.9874393050280.2112026.21-1.820.80824.921.660.87618.333.20.90643930615223.923603.873.950.9463.663.070.9524.812.70.9574393063094.42380-27.66-13.420.8222.690.440.99312.7-0.770.9834394042562.21410-31.52-15.430.748-0.53-0.660.94117.320.480.94943940614127.73250-1.89-15.970.865-0.011.750.9814.051.630.9444395060346.31896.675.760.640.9140.620.250.99512.311.080.9804395061485.4175010.79-11.430.897-0.67-2.520.9889.28-0.090.9784395081312.622524.7924.16-0.1311.320.740.88230.121.130.93043960317791190-8.1915.070.7030.515.540.94514.582.630.9234397060554172023.0919.760.7614.79-0.750.9606.111.520.94543980301114.61196.019.314.860.9507.46-1.380.96411.12-0.330.9744398050967.7201012.4919.440.9317.22-0.050.9627.413.310.9664399041554.41840-15.444.310.930-14.143.430.9452.261.120.94043990515661177.14-0.2719.990.849-8.7-1.160.966-1.671.030.9494399071549.7940-12.63-19.890.558-2.27-0.740.90412.26-1.210.9134399082591.51390-0.0820.150.820-5.970.010.9071.115.770.904

表3水吉流域29場(chǎng)洪水實(shí)時(shí)修正平均修正效果

評(píng)價(jià)模型評(píng)價(jià)指標(biāo)ΔQp/%ΔR/%DC新安江模型14.5312.180.743動(dòng)態(tài)系統(tǒng)響應(yīng)曲線法5.101.900.930遞推最小二乘法9.022.800.905

圖2 43880620場(chǎng)次洪水的誤差修正結(jié)果

圖3 43960317場(chǎng)次洪水的誤差修正結(jié)果

從圖2、3和表2、3可以看出:采用基于動(dòng)態(tài)系統(tǒng)響應(yīng)曲線的產(chǎn)流誤差修正方法之后,洪峰相對(duì)誤差、洪量相對(duì)誤差均有所減小,確定性系數(shù)均有所提高,平均修正效果也好于遞推最小二乘法。計(jì)算過程中發(fā)現(xiàn),采用遞推最小二乘法修正時(shí)會(huì)出現(xiàn)時(shí)刻流量值突漲的現(xiàn)象。如43960317場(chǎng)次的洪水,修正后雖然各項(xiàng)指標(biāo)有所提高但洪峰偏大且峰現(xiàn)時(shí)間提前,實(shí)際上修正后并不符合要求。而動(dòng)態(tài)系統(tǒng)響應(yīng)曲線方法則沒有這一缺點(diǎn)。從后者模擬效果得到了明顯的提升可以認(rèn)為,產(chǎn)流誤差動(dòng)態(tài)響應(yīng)曲線修正方法可以有效應(yīng)用于水吉流域。

3結(jié)論

產(chǎn)流誤差動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法是一種新的實(shí)時(shí)校正方法,其物理意義明確、結(jié)構(gòu)簡(jiǎn)明,既不損失預(yù)見期也不增加模型的參數(shù)個(gè)數(shù),且能有效保證預(yù)報(bào)精度。在水吉流域的應(yīng)用進(jìn)一步驗(yàn)證了該方法的合理性和實(shí)用性,將29場(chǎng)洪水的洪峰平均相對(duì)誤差從14.53%降低到5.10%,洪量平均相對(duì)誤差從12.18%降低到1.90%,平均確定性系數(shù)從0.743提高到0.930。與遞推最小二乘法比較也體現(xiàn)了其優(yōu)越性,整體修正效果相對(duì)較好。雖然該方法也存在問題,如U矩陣的穩(wěn)定性問題等,但值得在更多流域進(jìn)行應(yīng)用驗(yàn)證推廣和改進(jìn),以推進(jìn)實(shí)時(shí)洪水預(yù)報(bào)誤差修正方法的發(fā)展。

參考文獻(xiàn):

[1]包為民. 水文預(yù)報(bào)[M]. 北京: 中國(guó)水利水電出版社, 2006.

[2]周軼, 李致家. 改進(jìn)最小二乘遞推算法的洪水預(yù)報(bào)應(yīng)用研究[J]. 水力發(fā)電, 2006, 32(8): 14- 16.

[3]覃光華, 丁晶. 基于人工神經(jīng)網(wǎng)絡(luò)的卡爾曼濾波實(shí)時(shí)校正技術(shù)[J]. 水力發(fā)電, 2002(11): 9- 12.

[4]包為民, 嵇海祥, 胡其美, 等. 抗差理論及在水文學(xué)中的應(yīng)用[J]. 水科學(xué)進(jìn)展, 2003, 14(4): 528- 532.

[5]WU Xiaoling, WANG Chuanhai, CHEN Xi, et al. Kalman Filtering Correction in real-time forecasting with hydrodynamic model[J]. Journal of Hydrodynamics, Ser. B, 2008, 20(3): 391- 397.

[6]KOMMA J, BL?SCHL G, RESZLER C. Soil moisture updating by Ensemble Kalman Filtering in real-time flood forecasting[J]. Journal of Hydrology, 2008, 357(3): 228- 242.

[7]瞿思敏, 包為民. 實(shí)時(shí)洪水預(yù)報(bào)綜合修正方法初探[J]. 水科學(xué)進(jìn)展, 2003, 14(2): 167- 171.

[8]HINO M. Runoff forecasts by linear predictive filter[J]. Journal of the Hydraulics Division, 1970, 96(3): 681- 702.

[9]BIONDI D, De LUCA D L. A Bayesian approach for real-time flood forecasting[J]. Physics and Chemistry of the Earth, Parts A/B/C, 2012, 42/44: 91- 97.

[10]HSU M H, FU J C, LIU W C. Flood routing with real-time stage correction method for flash flood forecasting in the Tanshui River, Taiwan[J]. Journal of Hydrology, 2003, 283(1): 267- 280.

[11]MADSEN H, SKOTNER C. Adaptive state updating in real-time river flow forecasting—A combined filtering and error forecasting procedure[J]. Journal of Hydrology, 2005, 308(1): 302- 312.

[12]PAGANO T C, WANG Q J, HAPUARACHCHI P, et al. A dual-pass error-correction technique for forecasting streamflow[J]. Journal of Hydrology, 2011, 405(3): 367- 381.

[13]YU Z, LIU D, Lü H, et al. A multi-layer soil moisture data assimilation using support vector machines and ensemble particle filter[J]. Journal of Hydrology, 2012, 475: 53- 64.

[14]司偉, 包為民, 瞿思敏. 洪水預(yù)報(bào)產(chǎn)流誤差的動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法[J]. 水科學(xué)進(jìn)展, 2013, 24(004): 497- 503.

[15]林虹. 閩江流域?qū)崟r(shí)洪水預(yù)報(bào)系統(tǒng)模型研究[D]. 南京: 河海大學(xué), 2005.

[16]張章新. 閩江流域水文特性分析[J]. 水文, 2000, 20(6): 55- 58.

(責(zé)任編輯陳萍)

收稿日期:2015- 06- 23

基金項(xiàng)目:國(guó)家自然科學(xué)基金面上項(xiàng)目(51279057/41371048/51479062);國(guó)家自然科學(xué)基金重大項(xiàng)目(51190090);國(guó)家自然科學(xué)基金重大項(xiàng)目子項(xiàng)目(51190091);中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金資助(2014B35314)

作者簡(jiǎn)介:束慧連(1989—),女,江蘇鹽城人,碩士研究生,主要從事水文物理規(guī)律模擬、水文預(yù)報(bào)方面研究.

中圖分類號(hào):P338(257)

文獻(xiàn)標(biāo)識(shí)碼:A

文章編號(hào):0559- 9342(2016)04- 0022- 05

Application of Runoff Error Correction Based on Dynamic System Response Curve

SHU Huilian1, BAO Weimin1,2, YANG Shanshan1, LIU Kexin1, ZHAO Liping1, YANG Xiaoqiang1

(1. College of Hydrology and Water Resources, Hohai University, Nanjing 210098, Jiangsu, China; 2. State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, Jiangsu, China)

Abstract:The correction method of Dynamic System Response Curve (DSRC) is a new real-time correction method that is based on differentiation and is with clear physical meaning. This method can trace back to the source of error. The dynamic system response curve of runoff is calculated and then it is used to correct the runoff, and eventually, the correction of flood process is realized. To verify the correction effectiveness of DSRC method, it is combined with the Xin’anjiang model and is tested in Shuiji Basin of Minjiang River, and the method is also compared to the recursive least square method. Total 29 floods in Shuiji Basin are calculated and corrected, and the results of correction show that the correction method of dynamic system response curve has a remarkable effect in improving the accuracy of flood forecasting and it is more effective than the recursive least square method.

Key Words:runoff error; dynamic system response curve; Shuiji Basin

主站蜘蛛池模板: 国产一级妓女av网站| 国产精品性| 超碰色了色| 免费在线播放毛片| 亚洲综合天堂网| 91久久偷偷做嫩草影院电| 69综合网| 色偷偷一区二区三区| 久久人人97超碰人人澡爱香蕉| www.av男人.com| 在线日韩一区二区| 欧美综合成人| 人人91人人澡人人妻人人爽| 久久网欧美| 中文字幕永久在线观看| 亚洲人成高清| 国产精品吹潮在线观看中文| 亚洲精品第五页| 国产91成人| 无码久看视频| 中国毛片网| 精品国产亚洲人成在线| 无码精品国产dvd在线观看9久| 国产精品片在线观看手机版| 久久精品国产91久久综合麻豆自制| 亚洲日韩精品伊甸| 国产交换配偶在线视频| 在线视频精品一区| 精品伊人久久久香线蕉| 四虎成人精品在永久免费| 日韩美一区二区| 亚洲二区视频| 国产精品尹人在线观看| 国产亚洲视频在线观看| 国产日韩欧美成人| 日本高清成本人视频一区| 成人小视频网| 国产日本欧美亚洲精品视| 免费在线观看av| 国产色婷婷| 日本免费精品| 视频二区欧美| 高潮毛片免费观看| 在线精品视频成人网| 国产精品观看视频免费完整版| 日韩精品毛片| 久久综合丝袜日本网| 国产在线第二页| 尤物精品国产福利网站| 直接黄91麻豆网站| 在线免费无码视频| 国产激情在线视频| 无码人中文字幕| 日本亚洲国产一区二区三区| 亚洲日韩精品无码专区97| 久久国产精品国产自线拍| 国产成人av一区二区三区| 国产精品视屏| 亚洲第一黄片大全| 亚洲伊人电影| 女人18毛片水真多国产| 成人在线欧美| 一本色道久久88综合日韩精品| 亚洲精品欧美日本中文字幕| 亚洲欧州色色免费AV| 精品国产污污免费网站| 五月激情婷婷综合| 成人福利在线视频| 精品撒尿视频一区二区三区| 在线观看的黄网| 国产精品女主播| 免费在线国产一区二区三区精品| 中文字幕日韩久久综合影院| 一级毛片在线播放免费| 黄片一区二区三区| 中文字幕资源站| 亚洲成人黄色网址| 国产在线视频福利资源站| 国产精品13页| 免费啪啪网址| 国产精品综合色区在线观看| 噜噜噜久久|