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

基于小波變換-SG濾波融合算法的平寨水庫場次洪水反推入庫流量去噪研究

2025-09-02 00:00:00高林宋金偉李作孝
人民珠江 2025年7期

中圖分類號:TV1 文獻(xiàn)標(biāo)識碼:A 文章編號:1001-9235(2025)07-0053-08

Abstract:Inresponsetoabnormalfuctuationsofwaterflowintheshapeofa\"sawtoth\"inthereverseinferenceofwaterbalance,a wavelettransform Savitzky-Golayfilterfusionalgorithmwasproposed.Thealgorithmfirstperformedwaveletdecompostionand thresholddenoisingpreprocessingontheresearchdatabysetingdiferentparametercombinations.Waveletbasisfunctionsand decompositionlevelswerethenselected.Basedonthis,theSavitzky-Golayfilterwasusedtoperformsecondarydenoisingonthe preprocesseddata,tofurtereliminateresidallow-leveloiseandimproveteoverallsmoingefect.inally,teinverseavelet transformwasemployedtoreconstructthecompletesignalprocess.Thewindowlengthandpolynomialorderwerefurtheroptimized throughdenoisingevaluationindicatorstoachevewavelettransformSavitzky-Golayflterfusiondenoising.Thismethodwasapplied tothehistoricalfloodeventsofngzhaiReservoirtoinfertheiflow,andtherelativeerroroffodpeak,peakimediercend smoothness were setasevaluation indicators.According totheresults,comparedwiththe fivepointthreeordersmothingmethod, slidingaveragemethod,andwavelettresholddenoisingmethod,thewavelettransformSGfilteringfusionalgorithmsignificantly improvesthesmoohnessanddenoisingefectwhileensuringthecharacteristicsofthereverseflowpeakinformationofteinflow.The researchresultsareof greatsignificancefortherefinedmanagementandoptimizedoperationofPingzhaiReservoirandserveasa

reference for similar projects.

Keywords:inboundtrafic;igzag;smothenoising;waveletthresholddenoisingmetod;waveletthreshlddenoisingmethod;fusior algorithm of wavelet transform and Savitzky-Golay filter

入庫流量是水文預(yù)報、水庫調(diào)度管理的重要參數(shù)之一。自前,大多數(shù)水庫的入庫流量難以直接通過監(jiān)測獲取,通常由水量平衡法反推獲得,由于水位波動誤差、水位庫容曲線精度、出庫流量測量誤差以及動庫容等各種誤差疊加因素的影響,反推獲得的入庫流量容易呈現(xiàn)鋸齒狀波動,甚者在枯水期還時常出現(xiàn)入庫流量為“負(fù)\"值現(xiàn)象,這種入庫流量的異常波動,容易引起調(diào)度人員判斷失誤。針對水庫人庫流量去噪問題,國內(nèi)外專家學(xué)者以往開展過大量的研究,常用方法有五點三次平滑法、滑動平均法及小波閾值去噪法等,例如武煒等[基于五點三次平滑法對廣西多個水電站反推人庫流量進(jìn)行了平滑改進(jìn)研究,金錚等2采用離散小波閾值去噪算法、五點三次平滑算法、Stein無偏風(fēng)險估計閾值法對三峽水庫反推入庫流量進(jìn)行了平滑消噪研究,曾凡林等[3從計算時段長、滑動平均算法及庫水位代表性等方面針對丹江口水庫入庫流量波動問題進(jìn)行了研究。

以往入庫流量去噪研究成果中,有的側(cè)重于溯因分析,從造成入庫流量“鋸齒狀”波動的水位觀測、水位-庫容曲線、泄洪設(shè)施泄流能力曲線、發(fā)電機(jī)組NHQ特性曲線等方面出發(fā),消除各類誤差疊加最終實現(xiàn)平滑去噪[3-5];有的側(cè)重于利用各種數(shù)值濾波算法對反推入庫流量進(jìn)行去噪平滑處理[1-2.6-8]其中,滑動平均法能一定程度消除入庫流量鋸齒狀波動,使入庫流量過程更平滑,但存在一定程度的“削峰\"現(xiàn)象,會造成入庫流量的過度“坦化”,且平滑處理后的入庫流量過程存在滯后現(xiàn)象[7]。五點三次法對于非平穩(wěn)的數(shù)據(jù),特別是對于陡增陡降數(shù)據(jù),其平滑效果較差,這也在一定程度上影響了其對密集鋸齒狀入庫流量過程的平滑效果;小波閥值去噪[9-12]因其在時域和頻域的多尺度分辨特性,局部優(yōu)化效果較好,但小波去噪中小波基函數(shù)、分解層數(shù)的選取和閾值的設(shè)定需要進(jìn)行優(yōu)選,且整體

平滑效果有待提升。

目前,SG濾波被廣泛地運用于地理信息、光譜圖像數(shù)據(jù)平滑去噪[13-15],它最大的特點在于保留原始數(shù)據(jù)特征的同時整體平滑效果較好,算法實現(xiàn)不需要頻域轉(zhuǎn)換。針對上述情況,從小波變換與SG濾波算法的優(yōu)勢出發(fā),本文提出了基于小波變換-SG濾波融合的算法,并創(chuàng)新用于解決水庫入庫流量去噪問題;經(jīng)應(yīng)用于平寨水庫場次洪峰流量推算,表明該方法在保留入庫流量洪峰特征值的同時提高了整體平滑效果,研究成果對于平寨水庫精細(xì)化管理與優(yōu)化調(diào)度具有重要意義,供類似項目參考。

1研究方法

1.1 小波變換算法

小波變換理論由法國從事石油信號處理的工程師J.Morlet在1974年首次提出,是一種基于時間和頻率的局部變換的信號分析方法,能有效地從信號中提取信息,在時頻分析中可以表征信號局部特征,是一種窗口大小固定不變但形狀可以改變的時頻局部分析方法,可以處理非平穩(wěn)信號帶強(qiáng)噪聲信號[1];也是一種在信號處理中廣泛應(yīng)用的多分辨分析方法,主要目的是針對信號的不同頻率成分采用相應(yīng)的分辨率進(jìn)行有效分析[12]

采用小波變換去噪處理流程如下: ① 分析含噪信號特性,選擇適切的小波基函數(shù)和分解層數(shù),對含噪信號進(jìn)行小波分解; ② 選擇一個合適的閾值和閾值函數(shù)對小波分解產(chǎn)生的高頻系數(shù)進(jìn)行量化,抑制噪聲的小波系數(shù),小波系數(shù)幅值較大的被保留,幅值較小的被置為零; ③ 將經(jīng)過閾值處理的小波系數(shù)采取離散小波逆變換處理手段進(jìn)行重組,即去噪后的信號。離散小波變換定義為:

式中: Wf(j,k) 為離散小波變換系數(shù) σjf(σt) 為原始信號,本文中為入庫流量, m3/s;ψt 為小波基函數(shù) ψt

復(fù)共軛函數(shù);j為分解水平,也稱時間尺度水平;k為時間位置因子; a0 和 b0 均為常數(shù)。

為便于分析對比,本文分析計算時直接采用經(jīng)典的軟閾值處理變換公式:

式中: Wk 為經(jīng)過閾值處理的小波系數(shù); T 為閾值。小波逆變換重構(gòu)公式為:

式中: F(t) 為經(jīng)小波去噪后的入庫流量;其余變量意義同上。

1.2SG濾波算法

SG濾波由Savitzky和Golay在1964年首次提出,是一種用于時間序列數(shù)據(jù)平滑和降噪的數(shù)字信號處理技術(shù),它提供了一種通過將多項式函數(shù)擬合到數(shù)據(jù)的一小部分來估計噪聲數(shù)據(jù)中潛在趨勢或模式的方法。該方法很大程度上在保留原數(shù)據(jù)信息的同時,還能體現(xiàn)較好的平滑效果。

SG濾波本質(zhì)上是一種卷積滑動窗口加權(quán)平均算法,其加權(quán)系數(shù)取決于在一個濾波窗口內(nèi)給定高階多項式的最小二乘擬合次數(shù),見式(4):

式中: Yj 為擬合后的輸出數(shù)據(jù); Xj+i 為待擬合的數(shù)據(jù);ai 為擬合的平滑系數(shù),即在第 i 期原始數(shù)據(jù) Xj+i 在平滑窗口的權(quán)重; m 為窗口長度,奇數(shù); n 為滑動窗口的數(shù)據(jù)個數(shù), n=2m+1 。

1.3小波變換-SG濾波融合算法

針對反推入庫流量呈現(xiàn)異常波動的現(xiàn)象,本文提出利用小波變換-SG濾波融合去噪的數(shù)值濾波方法,該方法意在融合小波閾值去噪和SG濾波在去噪方面的各自優(yōu)勢,摒除單一濾波算法的局限性。具體實現(xiàn)步驟如下: ① 小波變換預(yù)處理,選定適宜的小波基函數(shù),利用小波變換對反推入庫流量系列進(jìn)行多尺度分解預(yù)處理,通過去噪效果量化指標(biāo)展開分析評價進(jìn)行小波變換參數(shù)優(yōu)選; ② 應(yīng)用SG濾波,以第 ① 步首次優(yōu)選參數(shù)為基礎(chǔ)將小波閾值去噪?yún)?shù)統(tǒng)一,做到控制唯一變量,利用SG濾波對前述預(yù)處理的數(shù)據(jù)進(jìn)行二次去噪,根據(jù)去噪需求調(diào)整濾波器窗口長度及多項式階數(shù); ③ 信號重構(gòu),利用逆小波變換將經(jīng)過SG濾波處理后的入庫流量系列重構(gòu)成一個完整的信號,得到不同參數(shù)組合去噪處理后的入庫流量系列; ④ 融合參數(shù)優(yōu)選,通過評價指標(biāo)對不同SG濾波的窗口長度及多項式階數(shù)組合方案的去噪效果進(jìn)行二次優(yōu)選,選定最優(yōu)去噪?yún)?shù)組合,實現(xiàn)小波變換-SG濾波融合去噪。

1.4效果評價指標(biāo)

針對汛期場次洪水反推人庫流量過程平滑去噪效果,本文選取洪峰信息刻畫指標(biāo)和平滑度指標(biāo)共2個指標(biāo)進(jìn)行評價。

a)洪峰信息刻畫指標(biāo)。洪峰信息刻畫指標(biāo)反映洪峰信息保留程度,包括洪峰相對誤差 δQ 和峰現(xiàn)時差 Δt 。洪峰相對誤差 δQ 可明確修正前后洪峰的誤差變化規(guī)律,其值越接近0表示洪峰信息保留程度越好;峰現(xiàn)時差 Δt 表示平滑后洪峰出現(xiàn)的時間點與原數(shù)據(jù)系列的時間點差值,其值越接近0表示洪峰信息保留程度越好。其中洪峰相對誤差 計算見式(5):

式中: QFi 為時段 i 水庫入庫流量修正值, m3/s : Qi 為時段i水庫入庫流量反推值, m3/s 。

b)平滑度指標(biāo)。平滑程度指標(biāo)由平滑后流量與原始流量的一階差分的標(biāo)準(zhǔn)差比值構(gòu)成,范圍在0~1 ,且越接近1表明平滑后流量序列相比原反推流量的平滑程度更高。

式中: Esm 表示平滑度,其值范圍為 0~1 ,值越大表示曲線越平滑,修正效果越好。

2 應(yīng)用對象概況

2.1平寨水庫基本情況

黔中工程是貴州省首個大型跨地區(qū)、跨流域的長距離調(diào)水工程,地處長江和珠江兩大流域分水嶺地帶,受水區(qū)涉及貴州省貴陽市、安順市、六盤水市、畢節(jié)、黔南自治州等4市1州的10個縣(區(qū))和貴陽、安順市區(qū)。工程以灌溉、城鎮(zhèn)供水為主,兼顧發(fā)電等綜合利用,主要解決貴陽市、安順市城市供水和沿線4個地州市7個縣42個鄉(xiāng)鎮(zhèn)51.17萬畝(1畝約等于 667m2 )耕地的灌溉,及5個縣城、28個鄉(xiāng)鎮(zhèn)供水,解決35萬人、31.5萬頭大牲畜的飲水。黔中工程水源樞紐為平寨水庫,位于六枝特區(qū)與織金縣交界的三岔河中游木底河段平寨附近,是烏江流域上游三岔河的龍頭水庫,距貴陽 200km ,距六枝特區(qū) 45km ,壩址以上集水面積 3492km2 ,工程規(guī)模為I等大(1)型,大壩為混凝土面板堆石壩,壩頂長355m ,最大壩高 157.5m ,水庫正常蓄水位 1331.00m ,死水位 1305.00m ,調(diào)節(jié)庫容4.48億 m3 ,總庫容10.89億 m3

2.2數(shù)據(jù)來源及選取

一般而言,汛期與枯水期入庫流量均存在鋸齒狀波動,但汛期更關(guān)注洪峰特征信息,枯水期更關(guān)注負(fù)值率情況,導(dǎo)致2個不同時期入庫流量去噪研究方式方法有所區(qū)別。本文從平寨水庫汛期運行遭遇過的不同年份場次洪水出發(fā),選取其中\(zhòng)"20210917\"\"20220526\"\"20220619\"等8場洪水反推入庫流量系列作為研究對象,研究洪水序列包含了平寨水庫不同年份小、中、大等洪峰量級的完整起漲、退水過程。選取的各場次洪水過程見圖1,由圖可見各場次洪水過程呈鋸齒狀變化,導(dǎo)致部分?jǐn)?shù)據(jù)失真,不利于汛期水庫精細(xì)化調(diào)度。

圖1選取的平寨水庫各場次洪水反推入庫流量過程系列

Fig.1Series of selectedPingzhaiReservoirflood dischargeprocessforeachevent

3小波變換-SG濾波算法融合研究

3.1小波變換參數(shù)取值優(yōu)選

小波變換去噪主要在于選擇適合的小波基函數(shù)及分解級數(shù),基于平寨水庫入庫流量系列形態(tài)及去噪目的,本文選擇適用于非線性(非對稱性)特點的Daubechies(db)系列小波基函數(shù),初選有db2、db3、db4、db5等4種小波基函數(shù)進(jìn)行小波分解;考慮到一般分解級數(shù)達(dá)到3已經(jīng)足夠平滑大部分噪聲信號,因此本文小波分解級數(shù)擬為1、2、3。應(yīng)用上述小波基函數(shù)及小波分解級數(shù)進(jìn)行參數(shù)組合,并利用小波閾值去噪進(jìn)行預(yù)處理,以洪峰信息刻畫指標(biāo)及平滑度指標(biāo)量化去噪效果,以“20220526\"場次洪水為例,各參數(shù)組合方案去噪效果對比見圖2。

3 原始反推流量 400 原始反推流量 400 原始反推流量小波去噪后流量 200 小波去噪后流量 300 小波去噪后流量100 1000 010 20 30 40 50 60 0 10 20 30 40 50 60 0 10 20 30 40 50 60時段/h 時段/h 時段/ha)小波基函數(shù)db2,分解級數(shù)1 b)小波基函數(shù)db2,分解級數(shù)2 c)小波基函數(shù)db2,分解級數(shù)3原始反推流量 _S 400 原始反推流量 (s 400 原始反推流量

3.2SG濾波參數(shù)取值優(yōu)選

根據(jù)公式設(shè)定,SG濾波窗口長度必須大于多項式階數(shù)且為奇數(shù),較長的窗口可以更好地去除高頻噪聲但保留數(shù)據(jù)特征細(xì)節(jié)越差;較高的多項式階數(shù)可以更好地保留數(shù)據(jù)的高頻特征,然而可能會導(dǎo)致過擬合,一般情況下,多項式階數(shù)設(shè)置為3已經(jīng)足夠平滑大部分噪聲信號。綜合分析,本文設(shè)定窗口長度3、5、7與多項式階數(shù)1、2、3的組合進(jìn)行參數(shù)優(yōu)選,以洪峰信息刻畫指標(biāo)及平滑度指標(biāo)量化去噪效果,以“20220526\"場次洪水為例,各參數(shù)組合方案去噪效果對比見圖3。

3.3小波變換-SG濾波算法融合

融合應(yīng)用時反推入庫流量系列經(jīng)不同參數(shù)的小波變換參數(shù)組合預(yù)處理后,對去噪效果進(jìn)行首次篩選,首次篩選小波基函數(shù)db3,分解級數(shù)1效果最優(yōu),以此為基礎(chǔ)將小波閾值去噪?yún)?shù)統(tǒng)一,做到控制唯一變量,再進(jìn)一步經(jīng)SG濾波進(jìn)行二次去噪,提高整體平滑效果,各參數(shù)組合方案融合法針對研究數(shù)據(jù)所列8場洪水反推流量序列去噪后效果對比見表1,經(jīng)優(yōu)選小波變換-SG濾波算法融合效果以參數(shù)小波基函數(shù)db3,分解級數(shù)1與SG濾波窗口長度5,多項式階數(shù)2組合最優(yōu)。為便于直觀對比,截取單一算法及本文融合算法針對“20220526\"場次洪水去噪的效果對比見圖4。

原始反推流量 : 400- 原始反推流量 S. 400- 原始反推流量 SG去噪后流量 SG去噪后流量 SG去噪后流量

表1不同參數(shù)組合的小波變換-SG濾波融合算法去噪效果對比

differentparametercombinations

Tab.1Comparisonofdenoisingeffectsoffusionalgorithm ofwavelettransformandSavitzky-Golayfilterwith

圖4單一算法及本文融合算法針對\"20220526\"場次洪水去噪效果對比

Fig.4Comparisonofdenoising effectsofasinglealgorithm and the fusionalgorithm proposedin thisarticle for the \"20220526\" flood event

4結(jié)果討論

為了評估本文提出的小波變換-SG濾波融合算法去噪效果,進(jìn)一步對比分析了五點三次平滑法、滑動平均法、小波閾值去噪法對3.2節(jié)選取的所有場次洪水反推入庫流量系列進(jìn)行平滑去噪,利用洪峰相對誤差、峰現(xiàn)時差及平滑度作為去噪效果量化指標(biāo),結(jié)果見表2。

表2不同方法針對平寨水庫各場次洪水反推入庫流量系列平滑去噪成果對比

Tab.2Comparisonofsmoothdenoisingresultsofdifferent methodsforpredictingtheinletflowseriesofvarious floodsinPingzhaiReservoir

由表2可以看出,各種去噪方法均能對場次洪水反推入庫流量進(jìn)行平滑去噪。在洪峰信息刻畫方面,五點三次平滑法對洪峰流量保留程度最好,小波變換-SG濾波融合法和小波閾值去噪法依次次之,滑動平均法過度坦化了洪峰,去噪之后造成的洪峰絕對偏差分別是其余3種方法的2.07倍、1.88倍、1.97倍。在峰現(xiàn)時差方面,滑動平均法去噪后造成的偏差最大,平均偏差達(dá)到了 1.25h ,五點三次平滑法、小波閾值去噪法、小波變換-SG濾波融合法造成的偏差均較小,分別為 0.25,0.38,0.38h ,時間偏差在 0.5h ,小于本文設(shè)定偏差閾值,對洪水實時調(diào)度也不會產(chǎn)生明顯影響。在平滑度方面,滑動平均法最優(yōu)但過度坦化了場次洪水過程,五點三次平滑法和小波閾值去噪法平滑效果接近,而在應(yīng)用小波變換-SG濾波融合優(yōu)化去噪后,在很好的保留洪峰信息特征的同時,相比單一小波變換去噪直接將平滑度提升了0.09,明顯改善了去噪效果。

另外,五點三次平滑法、滑動平均法、小波閾值去噪法、SG濾波-小波變換融合法對研究數(shù)據(jù)去噪后造成的流量偏差絕對均值分別為 0.05%.0.18% !0.03%.0.02% ,4種方法造成的水量平衡偏差均較小,尤其以SG濾波-小波變換融合法的偏差基本可忽略不計,最大程度保留了原始場次洪水入庫流量過程。

為便于直觀對比,本文截取不同方法針對其中“20220526\"典型場次洪水反推人庫流量過程去噪效果對比見圖5。

圖5不同方法對\"20220526\"場次洪水反推入庫流量平滑去噪效果對比

Fig.5Comparisonofdenoisingoptimizationeffectsofdifferentmethodsonthereverseinferenceofinflowflowforflood event“20220526\"

5 結(jié)論

a)提高水庫入庫流量精度是一個系統(tǒng)工程,目前各類方法互有利弊,本文從數(shù)值濾波融合去噪的角度率先提出了小波變換-SG濾波融合算法,該方法吸收了小波變換和SG濾波在數(shù)據(jù)去噪方面的各自優(yōu)勢,在保證了反推入庫流量洪峰信息特征的同時,顯著提升了平滑度,融合算法具有較好的應(yīng)用推廣價值。

b)相比于五點三次平滑法、滑動平均法、小波閾值去噪法,小波變換-SG濾波融合法在洪峰信息刻畫、峰現(xiàn)時差、平滑度等方面綜合效果最佳。

c)應(yīng)用小波變換-SG濾波融合法對不同量級場次洪水反推入庫流量進(jìn)行平滑去噪后,對場次洪水總體水量平衡的影響基本可忽略不計,應(yīng)用該方法幾乎不會帶來“后遺癥”。

隨著各種數(shù)值濾波新算法的不斷出現(xiàn),針對人庫流量異常波動的融合去噪新機(jī)制有待進(jìn)一步研究。

參考文獻(xiàn):

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

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

[3]曾凡林,牛文靜,張成孝.丹江口水庫入庫流量平滑修正方法研究[J].人民長江,2024,55(2):93-100.

[4]王俊莉.貴州烏江梯級水庫入庫流量計算分析及改進(jìn)方法[J].水利水電快報,2015,36(4):54-56.

[5]康傳雄,劉軒箕,王永文,等.反推水庫入庫流量的優(yōu)化修正方法研究[J].水電能源科學(xué),2024,42(6):158-161.

[6]張冀,張毅,周馳,等.入庫流量及其反演推算研究[J].云南水力發(fā)電,2022,38(8):29-33.

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

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

[9]彭楊,陳凱,紀(jì)昌明.向家壩水庫動庫容特性及影響分析[J].應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報,2018,26(2):239-248.

[10]周棟.水庫入庫流量鋸齒狀波動問題探討[J].陜西水利,2019,12(12):42-47.

[11]董博.小波變換技術(shù)在變形監(jiān)測中的應(yīng)用[J].測繪通報,2016(S2):150-152.

[12]于濤,劉國棟,李金耀,等.基于小波分析的壓力機(jī)振動信號處理與分析[J].鍛壓技術(shù),2022,47(2):152-157.

[13]吳迪,陳健,石滿,等.基于Savitzky-Golay濾波算法的FY-2F地表溫度產(chǎn)品時間序列重建[J].國土資源遙感,2019,31(2):59-65.

[14]魯一冰,劉文清,張玉鈞,等.一種自適應(yīng)層進(jìn)式Savitzky-Golay光譜濾波算法及其應(yīng)用[J].光譜學(xué)與光譜分析,2019,39(9):2657-2663.

[15]陳志剛,束炯.基于Savitzky-Golay方法的遙感影像融合[J].地理與地理信息科學(xué),2011,27(2):26-29.

(責(zé)任編輯:向飛)

主站蜘蛛池模板: yjizz国产在线视频网| 日韩成人高清无码| 色综合网址| 天天做天天爱天天爽综合区| 国产精品99一区不卡| 白丝美女办公室高潮喷水视频| 五月婷婷综合色| 亚洲AV无码精品无码久久蜜桃| 亚洲国产理论片在线播放| 亚洲天堂免费在线视频| 精品久久久久成人码免费动漫| 亚洲无码高清一区| 丁香亚洲综合五月天婷婷| 免费看一级毛片波多结衣| 国产杨幂丝袜av在线播放| 亚洲男人在线| 欧美性猛交xxxx乱大交极品| 久久人妻系列无码一区| 91精品网站| 免费不卡视频| 国产区福利小视频在线观看尤物| 2022国产91精品久久久久久| a级毛片网| 久久永久精品免费视频| 久久成人国产精品免费软件| 免费国产小视频在线观看| 亚洲人成网站观看在线观看| 国产精品lululu在线观看 | 91青草视频| 亚洲伦理一区二区| 国产偷国产偷在线高清| 免费高清a毛片| 日韩天堂视频| 欧美激情第一区| 国产另类视频| 国产女人18毛片水真多1| 国产欧美在线观看一区| 亚洲第一黄片大全| 日韩少妇激情一区二区| 97精品久久久大香线焦| 亚洲一级毛片在线观| 中文无码日韩精品| 亚洲色图欧美| 精品人妻一区二区三区蜜桃AⅤ| 五月婷婷丁香色| 国产视频一区二区在线观看| www.99精品视频在线播放| 青青草原国产免费av观看| 男女男免费视频网站国产| 无码免费视频| 无码福利日韩神码福利片| 日韩国产无码一区| 国产成人在线无码免费视频| 国产精品亚洲精品爽爽| 黄色不卡视频| 91小视频版在线观看www| 国产综合另类小说色区色噜噜| 色一情一乱一伦一区二区三区小说| 日韩经典精品无码一区二区| 欧美日韩成人| 久久综合五月| 国产女人在线| 亚洲欧美日韩成人高清在线一区| 97超爽成人免费视频在线播放| 国产精品私拍在线爆乳| 成人福利视频网| 国产91丝袜| 亚洲最猛黑人xxxx黑人猛交| 人禽伦免费交视频网页播放| av尤物免费在线观看| 美女亚洲一区| 国产亚洲视频免费播放| 午夜精品一区二区蜜桃| 国产爽妇精品| 91精品专区国产盗摄| 国产日产欧美精品| aaa国产一级毛片| 欧美一区二区人人喊爽| 国产成人综合日韩精品无码首页| 中文字幕伦视频| 国产精品毛片在线直播完整版| 一级毛片免费观看久|