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

卡洛特水電站入庫(kù)流量平滑修正方法研究

2024-12-31 00:00:00徐志劉天勇何松
水利水電快報(bào) 2024年11期

摘要:巴基斯坦卡洛特水電站采用靜庫(kù)容反推計(jì)算的入庫(kù)流量存在“鋸齒”問(wèn)題。為此,量化分析了卡洛特水庫(kù)水位觀測(cè)誤差對(duì)入庫(kù)流量的影響,并利用3點(diǎn)線(xiàn)性平滑法、5點(diǎn)線(xiàn)性平滑法、7點(diǎn)線(xiàn)性平滑法、五點(diǎn)三次函數(shù)擬合平滑法,對(duì)卡洛特水電站2023年8月和9月入庫(kù)流量中的“鋸齒”波動(dòng)進(jìn)行平滑。結(jié)果表明:卡洛特水電站正常蓄水位461 m時(shí),0.01 m的水位誤差造成的流量誤差為15 m3/s,超過(guò)允許入庫(kù)流量預(yù)報(bào)誤差(13 m3/s),庫(kù)水位觀測(cè)誤差對(duì)流量預(yù)報(bào)工作存在不利影響。五點(diǎn)三次平滑算法可以有效消除反推入庫(kù)流量數(shù)據(jù)中的“鋸齒”干擾成分,更全面地保留原入庫(kù)流量的水文特征,平滑結(jié)果較理想。研究成果可為調(diào)節(jié)庫(kù)容小、入庫(kù)流量預(yù)報(bào)精度要求高的中小水庫(kù)反推入庫(kù)流量提供參考。

關(guān)鍵詞:入庫(kù)流量; 靜庫(kù)容反推計(jì)算; 流量平滑修正; 五點(diǎn)三次平滑法; 卡洛特水庫(kù); 巴基斯坦

中圖法分類(lèi)號(hào):TV72

文獻(xiàn)標(biāo)志碼:A

DOI:10.15974/j.cnki.slsdkb.2024.11.003

文章編號(hào):1006-0081(2024)11-0021-05

0 引 言

巴基斯坦卡洛特電站工程設(shè)計(jì)初期考慮到水庫(kù)庫(kù)尾與上一梯級(jí)規(guī)劃的阿扎德帕坦(Azad Pattan)水電站壩址相接,采用水量平衡法計(jì)算入庫(kù)洪水,發(fā)現(xiàn)流量存在嚴(yán)重的“鋸齒”現(xiàn)象,增加了用其校核預(yù)報(bào)成果的難度。

入庫(kù)流量計(jì)算的準(zhǔn)確性問(wèn)題難以從測(cè)量手段解決,閔要武、張俊等[1-2]提出“動(dòng)庫(kù)容”計(jì)算的方式,在三峽水庫(kù)進(jìn)行了模擬并取得了一定的效果。然而,多數(shù)水電站不具備設(shè)立多個(gè)水庫(kù)水位站的條件。目前采用平滑算法對(duì)入庫(kù)流量結(jié)果進(jìn)行平滑修正是處理此類(lèi)問(wèn)題的常用手段[3-4],主要的平滑算法有 3 點(diǎn)、5 點(diǎn)、7 點(diǎn)線(xiàn)性平滑,五點(diǎn)三次函數(shù)擬合平滑,最小二乘曲線(xiàn)[5],樣條插值平滑[6],以及一些其他的濾波算法[7] 。本文以卡洛特水電站實(shí)際計(jì)算入庫(kù)流量為例,對(duì)比分析了 3 點(diǎn)、5 點(diǎn)、7點(diǎn)線(xiàn)性平滑法,以及五點(diǎn)三次函數(shù)擬合平滑法對(duì)入庫(kù)流量的修正效果。

1 工程概況

卡洛特水庫(kù)總庫(kù)容1.88億m3,正常蓄水位以下庫(kù)容1.52億m3,死庫(kù)容1.03億m3,調(diào)節(jié)庫(kù)容0.49億m3,具有日調(diào)節(jié)性能。卡洛特水電站的商業(yè)模式為容量電價(jià)機(jī)制,根據(jù)巴基斯坦中央電力采購(gòu)局對(duì)卡洛特水電站的發(fā)電指標(biāo)及水資源利用考核要求,當(dāng)水電站入庫(kù)流量不大于水電站最大發(fā)電引用流量時(shí),電站每小時(shí)均需要充分利用水資源。根據(jù)水庫(kù)運(yùn)行模式,當(dāng)卡洛特水庫(kù)入庫(kù)流量為624.2~1 248.4 m3/s時(shí),水庫(kù)運(yùn)行處于“徑流模式”,購(gòu)電協(xié)議要求水位需維持在正常蓄水位461 m運(yùn)行,即保持出入庫(kù)平衡,這對(duì)水庫(kù)入庫(kù)流量預(yù)報(bào)精度提出了極高的挑戰(zhàn)。

2 入庫(kù)流量誤差分析

基于水量平衡的反推入庫(kù)流量誤差來(lái)源主要有庫(kù)容曲線(xiàn)誤差、水位觀測(cè)誤差和閘門(mén)啟閉引起的泄流量變化[8],以上誤差導(dǎo)致反推入庫(kù)流量偏離實(shí)際值,流量過(guò)程線(xiàn)呈現(xiàn)“鋸齒”狀。卡洛特水電站屬于河道型水庫(kù),于2021年11月下閘蓄水,蓄水時(shí)間短,水庫(kù)淤積對(duì)庫(kù)容曲線(xiàn)影響有限。本文定量分析水位觀測(cè)誤差對(duì)反推入庫(kù)流量的影響。

水庫(kù)反推入庫(kù)實(shí)測(cè)流量,都是根據(jù)庫(kù)水位和總出庫(kù)流量來(lái)計(jì)算。根據(jù)水庫(kù)水量平衡方程,計(jì)算公式為

Q1+Q22Δt-q1+q22Δt=V2-V1(1)

式中:Q1,q1為時(shí)段初入庫(kù)、出庫(kù)流量,m3/s;Q2,q2為時(shí)段末入庫(kù)、出庫(kù)流量,m3/s;V1,V2為時(shí)段初、末的水庫(kù)蓄水量,m3,由庫(kù)水位查水庫(kù)庫(kù)容曲線(xiàn)獲得;Δt為時(shí)間步長(zhǎng),s。

對(duì)式(1)取時(shí)段平均值,有:

QΔt-qΔt=V2-V1(2)

式中:Q為時(shí)段平均入庫(kù)流量,m3/s;q為時(shí)段平均總出庫(kù)流量,m3/s;V1=V(H1),V2=V(H2),H1和H2分別為時(shí)段初始和末了時(shí)刻的庫(kù)水位,m,V1和V2為相應(yīng)時(shí)刻的庫(kù)容,m3,由庫(kù)水位查水庫(kù)庫(kù)容曲線(xiàn)獲得。

根據(jù)庫(kù)水位查算庫(kù)容,再由式(2)推算時(shí)段平均入庫(kù)流量,誤差存在放大現(xiàn)象,其放大的倍數(shù)與水庫(kù)水面面積成正比,與時(shí)段長(zhǎng)成反比。設(shè)時(shí)段初始和末了時(shí)刻的庫(kù)水位差為ΔH=H2-H1,H1和H2間水庫(kù)橫斷面的平均面積為A(單位為km2),則式(2)可表示為

Q=ΔH×A3600×Δt+q(3)

當(dāng)時(shí)段始末庫(kù)水位存在誤差(單位為m)時(shí),就有相應(yīng)的流量誤差ΔQ:

ΔQ=1000××A3.6×Δt(4)

對(duì)于現(xiàn)行的大多數(shù)水庫(kù)實(shí)時(shí)洪水預(yù)報(bào)系統(tǒng),時(shí)段間隔為1 h,則式(4)的誤差ΔQ與水面面積呈線(xiàn)性比例關(guān)系,卡洛特水電站流量誤差與水位誤差關(guān)系見(jiàn)表1。

卡洛特水電站正常蓄水位461 m,電站運(yùn)行水位長(zhǎng)期控制在460~461 m。庫(kù)水位為461 m時(shí),0.01 m的水位誤差,其流量差值為15 m3/s;0.02 m的水位誤差,其流量差值為31 m3/s;當(dāng)水位誤差為0.05 m時(shí),其流量差值就達(dá)到了77 m3/s。卡洛特水電站10月至次年2月多年平均流量約270 m3/s,按照5%的預(yù)報(bào)誤差,允許入庫(kù)流量預(yù)報(bào)誤差往往只有13 m3/s左右。因此,庫(kù)水位誤差變化0.01 m帶來(lái)的流量差值就足以影響預(yù)報(bào)精度評(píng)定結(jié)果。

3 研究方法

3.1 滑動(dòng)平均法

滑動(dòng)平均法在數(shù)學(xué)和統(tǒng)計(jì)學(xué)中是一種常用的數(shù)據(jù)平滑方法。通過(guò)對(duì)原始數(shù)據(jù)進(jìn)行平滑處理,可以減少數(shù)據(jù)中的噪聲和波動(dòng),使數(shù)據(jù)更加平穩(wěn)和可靠。基本原理是對(duì)時(shí)間序列資料進(jìn)行逐項(xiàng)推移計(jì)算,可按照3點(diǎn)、5點(diǎn)或7點(diǎn)等長(zhǎng)度逐步遞進(jìn),依次計(jì)算各長(zhǎng)度項(xiàng)數(shù)的時(shí)序平均值,以弱化隨機(jī)波動(dòng)的影響[9]。計(jì)算方法見(jiàn)式(5)。

Yk=∑pi=-qwiXk+ik=q+1,q+2,…,N-p(5)

式中:Y為平滑后的時(shí)間序列;X為原始時(shí)間序列;n為滑動(dòng)窗口長(zhǎng)度;N為時(shí)間序列的長(zhǎng)度。wi為加權(quán)系數(shù);p,q 為小于n 的任一正整數(shù),且p+q+1=n。

3.2 五點(diǎn)三次函數(shù)擬合平滑法

滑動(dòng)平均法[9]應(yīng)用領(lǐng)域比較廣泛,五點(diǎn)三次函數(shù)擬合平滑法在最小二乘法的基礎(chǔ)上加入了數(shù)據(jù)平滑技術(shù),使得該算法更加注重?cái)?shù)據(jù)的細(xì)節(jié)處理[10-11],且算法簡(jiǎn)單、運(yùn)算速度快。國(guó)內(nèi)水庫(kù)的相關(guān)研究[12-14]顯示,五點(diǎn)三次函數(shù)擬合平滑法在平滑入庫(kù)流量的問(wèn)題中取得了理想的結(jié)果。

設(shè)在2n+1個(gè)等距為h的節(jié)點(diǎn)x-n,x-n+1,…,x0, …,xn-1,xn上,數(shù)據(jù)分別為Y-n,Y-n+1,…,Y0,…,Yn-1,Yn,作交換t=(x-x0)/h,則原節(jié)點(diǎn)變?yōu)閠-n=-n, …,t0=0,…,tn=n,建立m次多項(xiàng)式,見(jiàn)式(6),系數(shù)由最小乘法確定,方差和見(jiàn)式(7)。

Yt=a0+a1t+a2t2+…+amtm(6)

F(a0,a1,…,am)=∑ni=-n∑mj=0ajtji-Yi2(7)

為使方差和最小,令:

Fak=0(k=0,1,…,m) (8)

可得方程組:

∑ni=-nYjtki=∑mj=0aj∑ni=-ntk+1i(9)

當(dāng)n取2(即5個(gè)節(jié)點(diǎn))、m取3時(shí),將計(jì)算所得的a0,a1,a2代入式(7),同時(shí)令t=-2,-1,0,1,2,得到五點(diǎn)三次平滑公式。設(shè)卡洛特水電站的入庫(kù)流量序列為Q1,Q2,…,Qt,…,Qn-1,Qn(ngt;5),將Qn代入五點(diǎn)三次平滑公式中的Yn,則得到五點(diǎn)三次平滑計(jì)算公式如下。當(dāng)節(jié)點(diǎn)個(gè)數(shù)超過(guò)5時(shí),除在序列兩端4個(gè)節(jié)點(diǎn)分別使用式(10)、式(11)、式(13)及式(14)進(jìn)行平滑,其余節(jié)點(diǎn)均用式(12)進(jìn)行平滑。

Q1=170(69Q1+4Q2+6Q3+4Q4-Q5)(10)

Q2=130(2Q1+27Q2+12Q3-8Q4+2Q5)(11)

Qt=135(-3Qt-2+12Qt-1+17Qt+12Qt+1-3Qt+2)(2lt;tlt;n-1)(12)

Qn-1=1352Qn-4-8Qn-3+12Qn-2+27Qn-1+2Qn(13)

Qn=170-Qn-4+4Qn-3-6Qn-2+4Qn-1+69Qn(14)

式中:Q1為t=1時(shí)經(jīng)五點(diǎn)三次平滑后的改進(jìn)流量,m3/s;Qt為當(dāng)2lt;tlt;n-1時(shí),平滑后的改進(jìn)流量,m3/s;Qn為t=n時(shí)平滑后的改進(jìn)流量,m3/s。

4 平滑效果評(píng)價(jià)指標(biāo)

本文引入平滑度和相關(guān)系數(shù)指標(biāo)評(píng)價(jià)平滑后的入庫(kù)流量。

(1) 平滑度由平滑后數(shù)據(jù)和原始數(shù)據(jù)差分平方和構(gòu)成。平滑度范圍為0~1,且越接近1曲線(xiàn)越光滑。

Esm=1-∑n-1i=1(Qde,i+1-Qde,i)2∑n-1i=1(Qori,i+1-Qori,i)2(15)

式中:Esm為平滑度,其值越大表示曲線(xiàn)越平滑;Qori,i,Qori,i+1 分別為第i,i+1時(shí)段平滑前的入庫(kù)流量,m3/s;Qde,i,Qde,i+1分別為第i,i+1時(shí)段平滑后的入庫(kù)流量,m3/s。

(2) 相關(guān)系數(shù)可以反映兩個(gè)隨機(jī)變量相關(guān)程度,記為P。P越接近于1,數(shù)據(jù)相關(guān)性越好。相關(guān)系數(shù)計(jì)算公式如式(16)所示。

P=1n-1∑ni=1[fi-uf][gi-ug]σfσg(16)

式中:uf,ug分別為f(i)和g(i)的均值;σf,σg分別為f(i)和g(i)的標(biāo)準(zhǔn)差;f(i),g(i)分別為i時(shí)段平滑前、后的入庫(kù)流量。

5 平滑結(jié)果對(duì)比分析

本文根據(jù)卡洛特水電站機(jī)組過(guò)流特性,結(jié)合電站“徑流模式”流量考核區(qū)間,選取2023年8月和9月的逐時(shí)入庫(kù)流量作為對(duì)象,對(duì)比目前常采用的 3 點(diǎn)、5 點(diǎn)、7 點(diǎn)線(xiàn)性平滑法,五點(diǎn)三次函數(shù)擬合平滑法的平滑度和相關(guān)系數(shù)(表2)。

從表2可以看出:4種方法得到的8月和9月入庫(kù)流量平滑度從高到低依次是7點(diǎn)線(xiàn)性平滑法、5點(diǎn)線(xiàn)性平滑法、3點(diǎn)線(xiàn)性平滑法、五點(diǎn)三次函數(shù)擬合平滑法。在相關(guān)系數(shù)指標(biāo)上,4種方法平滑后的相關(guān)系數(shù)均在0.89以上;平滑度最高的7點(diǎn)線(xiàn)性平滑法在9月的相關(guān)系數(shù)反而是最低的,說(shuō)明存在過(guò)度平滑現(xiàn)象;五點(diǎn)三次函數(shù)平滑在研究時(shí)段內(nèi)相關(guān)系數(shù)指標(biāo)最高,均在0.95以上,說(shuō)明較好地保留了原始數(shù)據(jù)的特征。

選擇卡洛特水電站“徑流模式”上限流量的時(shí)間段8月5~9日、下限流量時(shí)間段8月15~19日和枯水期流量時(shí)間段9月22~25日的流量平滑過(guò)程放大觀察(圖1~3)。

可以看出,所有方法都有平滑效果,基本能準(zhǔn)確反映對(duì)應(yīng)洪峰時(shí)間,各方法的洪量相等,但過(guò)度平滑帶來(lái)的“削峰填谷”效應(yīng)在電站運(yùn)行中會(huì)帶來(lái)一定風(fēng)險(xiǎn)。在洪峰流量過(guò)度削減情況下,電站調(diào)度若不能及時(shí)響應(yīng)可能會(huì)發(fā)生庫(kù)水位超限的情況。而低流量的“虛高”效應(yīng)也會(huì)導(dǎo)致庫(kù)水位急劇下降,導(dǎo)致水頭減小,影響發(fā)電效益。

卡洛特水電站在“徑流模式”下水庫(kù)需保持出入庫(kù)平衡,水位控制在461 m,過(guò)度平滑會(huì)使得電站運(yùn)行違反規(guī)定,引起購(gòu)電方罰款。 五點(diǎn)三次平滑算法可以將反推入庫(kù)流量數(shù)據(jù)中“鋸齒”現(xiàn)象的干擾成分消除,保持?jǐn)?shù)據(jù)的變化趨勢(shì),可以認(rèn)為該方法的平滑結(jié)果較理想。

6 結(jié) 論

本文根據(jù)卡洛特水電站特性,分析了庫(kù)水位的觀測(cè)誤差引起的反推入庫(kù)流量誤差,并基于3點(diǎn)線(xiàn)性平滑法、5點(diǎn)線(xiàn)性平滑法、7點(diǎn)線(xiàn)性平滑法、五點(diǎn)三次函數(shù)擬合平滑法,對(duì)卡洛特水電站2023年8月和9月入庫(kù)流量中的“鋸齒”問(wèn)題進(jìn)行平滑,同時(shí)采用平滑度、相關(guān)系數(shù)指標(biāo)對(duì)其平滑效果進(jìn)行評(píng)價(jià),得到結(jié)論如下。

(1) 卡洛特水電站正常蓄水位461 m時(shí),0.01 m的水位誤差,其流量差值為15 m3/s;流域10月至次年2月的多年平均流量約270 m3/s,按照5%的預(yù)報(bào)誤差,允許入庫(kù)流量預(yù)報(bào)誤差往往只有13 m3/s左右,庫(kù)水位觀測(cè)誤差對(duì)預(yù)報(bào)工作極為不利。

(2) 卡洛特水電站8月和9月入庫(kù)流量平滑度從高到低依次是7點(diǎn)線(xiàn)性平滑法、5點(diǎn)線(xiàn)性平滑法、3點(diǎn)線(xiàn)性平滑法、五點(diǎn)三次函數(shù)擬合平滑法;五點(diǎn)三次平滑法相關(guān)系數(shù)在0.95以上,表現(xiàn)最佳。綜合分析認(rèn)為,五點(diǎn)三次平滑算法不僅可以消除反推入庫(kù)流量數(shù)據(jù)中“鋸齒”現(xiàn)象的干擾成分,還可以更全面地保證原入庫(kù)流量的水文特征,其平滑結(jié)果較理想,在調(diào)節(jié)庫(kù)容小、入庫(kù)流量預(yù)報(bào)精度要求高的中小水庫(kù)具有適用價(jià)值。

參考文獻(xiàn):

[1] 閔要武,王俊,陳力.三峽水庫(kù)入庫(kù)流量計(jì)算及調(diào)洪演算方法探討[J].人民長(zhǎng)江,2011,42(6):49-52.

[2] 張俊,閔要武,陳力,等.基于動(dòng)庫(kù)容曲線(xiàn)的水庫(kù)調(diào)洪高水位查算方法研究[J].水文,2016,36(5):63-68.

[3] 唐海華,陳森林,趙云發(fā),等.三峽水庫(kù)入庫(kù)流量計(jì)算方法研究[J].中國(guó)農(nóng)村水利水電,2008(4):26-28.

[4] 陳建,李允軍,王建平,等.基于梯度自適應(yīng)優(yōu)化的入庫(kù)流量數(shù)據(jù)平滑算法研究[J].人民長(zhǎng)江,2021,53(3):92-97.

[5] 張俊,廖勝利,程春田,等.基于最小二乘曲線(xiàn)的水電站入流平滑處理[J].水電能源科學(xué),2011,29(4):55-56,154.

[6] 楊海星,黎浩,曹凈.基于樣條函數(shù)的監(jiān)測(cè)數(shù)據(jù)平滑算法與Matlab實(shí)現(xiàn)[J].長(zhǎng)春工程學(xué)院學(xué)報(bào)(自然科學(xué)版),2018(3):92-97,120.

[7] 楊金標(biāo),舒凱,張后來(lái),等.水庫(kù)水量平衡計(jì)算中消除水位跳變影響的過(guò)濾算法[J].人民長(zhǎng)江,2019,50(8):98-102.

[8] 徐楊,劉亞新,汪濤,等.基于LSTM的三峽水庫(kù)短期上游水位預(yù)測(cè)方法研究[J].水利水電快報(bào),2022,43(10):13-18.

[9] 李征,房宏才,柯熙政,等.滑動(dòng)平均法在MEMS陀螺信號(hào)趨勢(shì)項(xiàng)提取中的應(yīng)用[J].電子測(cè)量與儀器學(xué)報(bào),2019,33(7):43-49.

[10] 蔡天凈,唐瀚.Savitzky-Golay平滑濾波器的最小二乘擬合原理綜述[J].數(shù)字通信,2011(1):63-68,82.

[11] 李匡,柳俊,張子平,等.基于自適應(yīng)權(quán)重的滑動(dòng)平均入庫(kù)流量修正方法[J].中國(guó)防汛抗旱,2021,31(5):20-24.

[12] 武煒,陳標(biāo),吳劍鋒,等.基于五點(diǎn)三次平滑算法的入庫(kù)流量反推研究[J].水利水電技術(shù),2013,44(12):100-102.

[13] 金錚,張行南,夏達(dá)忠,等.基于多種平滑算法的三峽水庫(kù)反推入庫(kù)流量消噪[J].中國(guó)農(nóng)村水利水電,2021(3):48-53.

[14] 李凱凱.張峰水庫(kù)入庫(kù)流量推算方法與應(yīng)用研究[J].山西水土保持科技,2021(2):38-40,48.

(編輯:江 文)

Smoothing correction method for inflow into Karot Hydropower Station

XU Zhi,LIU Tianyong,HE Song

(Jingjiang Bureau of Hydrology and Water Resources Survey,Bureau of Hydrology of Changjiang Water Resources Commission,Jingzhou 434000,China)

Abstract:

The \"jagged\" issue exists in the inflow into the Karot Hydropower Station calculated by reversing the static reservoir capacity.We quantitatively analyzed the impact of the Karot Reservoir′s water level measurement error on the inflow and smoothed the jagged fluctuations in the inflow for August and September in 2023 by using the 3-point linear smoothing method,5-point linear smoothing method,7-point linear smoothing method,and the cubical smoothing algorithm with 5-point approximation.The results showed that when the normal storage level of Karot Hydropower Station was 461 meters,a water level error of 0.01 meters caused a flow error of 15 m3/s,exceeding the allowable inflow forecast error (13 m3/s).The reservoir water level measurement error had an adverse effect on the flow forecasting work.The quintic smoothing algorithm can effectively eliminate the \"jagged\" interference components in the reverse-calculated inflow data,comprehensively retain the hydrological characteristics of the original inflow,and the smoothing results were relatively more ideal.The research results can provide a reference for the small and medium reservoirs with small adjustable capacities and high requirements for inflow forecasting accuracy in reversing the inflow calculation.

Key words:

inflow into reservoir; calculation by reversing static reservoir capacity; flow smoothing correction; cubical smoothing algorithm with 5-point approximation; Karot Reservoir; Pakistan

主站蜘蛛池模板: 毛片网站观看| 成人国产精品网站在线看 | 国产精品思思热在线| 亚洲成人黄色在线| 无码粉嫩虎白一线天在线观看| 欧美在线视频a| 国产精品成人AⅤ在线一二三四| 日韩欧美视频第一区在线观看| 国产美女91视频| 国产欧美在线| 成人福利在线观看| 免费看黄片一区二区三区| 婷五月综合| 色婷婷电影网| 久久国产av麻豆| 亚洲男人的天堂网| 9久久伊人精品综合| 欧美中文字幕在线视频| 中文字幕久久波多野结衣| 日本色综合网| 国产精品亚洲一区二区三区z| 久久精品波多野结衣| 国产成人精彩在线视频50| 国产在线精品人成导航| 免费人成视网站在线不卡| 香蕉久久国产超碰青草| 九九九九热精品视频| 狠狠操夜夜爽| 激情综合激情| 最新无码专区超级碰碰碰| 亚洲swag精品自拍一区| 亚洲网综合| 中文字幕一区二区人妻电影| 国产成人精品视频一区二区电影 | 日本精品αv中文字幕| 999精品色在线观看| 亚洲视频在线观看免费视频| 久久国产精品影院| 高清乱码精品福利在线视频| 久久女人网| a亚洲视频| 久久综合九色综合97婷婷| 欧美yw精品日本国产精品| 欧美特黄一免在线观看| 国产网站一区二区三区| 国产夜色视频| 在线国产欧美| 丁香五月婷婷激情基地| 欧类av怡春院| 中文字幕在线播放不卡| 国产成在线观看免费视频| 毛片一级在线| 黄色污网站在线观看| 青青极品在线| 青草国产在线视频| 国产亚洲视频在线观看| 91九色最新地址| 国产亚洲视频在线观看| 欧美视频在线不卡| 久久久久久久久久国产精品| 国产在线自乱拍播放| 亚洲香蕉伊综合在人在线| 成人一区专区在线观看| 亚洲成在线观看| 国产精品亚洲精品爽爽| 欧美午夜在线视频| 广东一级毛片| 亚洲精品成人片在线观看| 精品视频在线一区| 国产激情影院| 亚洲色图另类| 亚洲成人动漫在线观看| 人妻丝袜无码视频| 999福利激情视频| 久久久久青草大香线综合精品| 亚洲综合狠狠| 亚洲第一av网站| 亚洲乱码在线视频| 中文字幕无线码一区| 国产色网站| 91小视频在线播放| 国产人成在线视频|