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

丹江口水庫入庫流量平滑修正方法研究

2024-05-20 00:00:00曾凡林牛文靜張成孝嚴方家邢雯慧
人民長江 2024年2期

摘要:準確的入庫流量計算是科學開展水庫調度的重要基礎。針對反推入庫流量存在的劇烈跳動現象,開展了入庫流量推求方法研究,論述了常用推求方法的原理及適用性,聚焦丹江口水庫庫容特點,采用控制變量法,針對反推時段數目、滑動平均算法及壩前水位代表性等影響靜庫容反推計算結果的主要因素進行了定量對比分析,比選提出了適用于丹江口水庫的入庫流量平滑修正計算方法。研究結果表明:丹江口水庫采用庫容反推法,耦合6 h反推、6 h滑動、多站水位算數平均值作為壩前水位可獲得兼具準確性和平滑性的入庫流量計算結果。所提方法在保證洪峰流量和峰現時間更接近真實情況的同時,入庫過程方差更小、平滑度更高,可為丹江口水庫科學精細化調度管理提供重要的技術支撐。

關 鍵 詞:入庫流量; 水量平衡; 靜庫容反推; 滑動平均; 丹江口水庫

中圖法分類號: TV697 文獻標志碼: A DOI:10.16232/j.cnki.1001-4179.2024.02.012

0 引 言

入庫流量一般指單位時間內不同地點匯集到水庫的水量,即入庫斷面洪水與入庫區間洪水之和[1-3]。入庫流量計算是水庫調度最重要的基礎工作之一,準確的入庫流量可為水庫科學管理、精細調度提供關鍵數據基礎[4-6]。

通常情況下,由于觀測水位誤差的放大效應、水庫動庫容等因素的影響,計算所得入庫流量易出現鋸齒狀“波動”,甚至出現負值等不合理現象,極大干擾水庫調度工作的順利開展,影響合理決策[7-10]。為避免入庫流量過程出現劇烈跳動等現象,國內外諸多學者針對入庫流量反演和平滑問題開展了大量研究,取得了一定進展。王世策等[4]針對大型水庫計算入庫流量波動過大問題開展了研究,針對水位代表性不足等問題提出了解決辦法;唐海華等[5]開展了三峽水庫入庫流量計算方法研究,提出了分段河槽蓄水量計算水庫合成入庫流量的方法;張俊等[6]針對水庫入庫過程推算結果的“鋸齒”現象,提出了基于最小二乘曲線的入庫平滑處理方法;王俊莉[7]針對烏江梯級水庫對比分析了幾種入庫流量計算方法,并將七點滑動平均法引入計算,明顯提升了入庫流量過程的平穩性;Deng等[11]考慮水庫入庫流量過程連續性,通過優化反推水庫入庫流量的解析公式,保證水庫入庫流量過程光滑。Zhang等[12]基于多種模型開發出一種機器學習算法并用于水庫入庫流量預報,取得了較好效果。

現有研究表明,包括動庫容分析、泄流曲線校正、機組NHQ關系曲線修正、梯級水量平衡關系率定、卡爾曼濾波校正、七點滑動平均校正等在內的諸多方法,均可以用于入庫流量計算過程,但平滑效果“因庫而異”[13-14]。因此,在實際應用過程中,需要針對水庫特性及報汛能力擇優考慮。本文立足南水北調中線水源工程丹江口水庫的實際生產需求,針對其入庫流量推求方法開展研究,擬確定一種合理可行的計算方法以獲得準確、平滑的入庫過程,為丹江口水庫科學精細調度提供技術支撐。

1 研究區域概況

漢江是長江中游的重要支流,也是長江流域豐枯變化最大的支流之一[15]。漢江干流發源于秦嶺南麓,經漢中盆地與褒河匯合后始稱漢江,于武漢匯入長江,全流域集水面積15.9萬km2,干流全長1 577 km,大致呈東西向,兩岸坡陡谷深,河道灘多水急,平均比降在0.6‰以上;流域包括陜西省、河南省、湖北省、四川省及重慶市一部分,北以秦嶺及外方山與黃河為界,東北以伏牛山及桐柏山與淮河為界,西南以米倉、大巴山、荊山與嘉陵江、沮漳河相鄰,東南為廣闊平原;流域支流有褒河、湑水河、子午河、牧馬河、任河、嵐河、月河、旬河、夾河、堵河、丹江等。

丹江口水庫以上為漢江上游,位于秦嶺與大巴山之間,是南水北調中線工程水源區,長約925 km,集水面積9.52萬km2。漢江上游屬副熱帶季風區,流域多年平均年降水量約為700~1 100 mm,降水量年內分配不均勻,5~10月降水占全年的70%~80%,且夏、秋季洪水分期明顯,年最大洪水多發生在7月或9月。

丹江口水利樞紐位于湖北省丹江口市,漢江干流與支流丹江匯合處下游約800 m,是由漢江與支流丹江兩個庫區組成的并聯式水庫,控制流域面積95 200 km2。漢江上游干流及主要支流目前已建有石泉、安康、潘口、黃龍灘等多座大型水利樞紐,這些水電站的建成運行對丹江口水庫入庫徑流特性產生了一定的影響。丹江口水庫以上流域水系及主要水文站、水庫分布如圖1所示。

2 入庫流量推求原理及方法

對于已建水庫,入庫流量的計算方法可以分為正向計算的流量疊加法和反向計算的庫容反推法[3-6]。前者是一種順演方法,通常用于預報入庫流量,后者是一種反推方法,通常用于計算報汛入庫流量。本節分別針對兩類方法介紹其基本原理和計算流程,以確定適用于丹江口水庫的入庫流量計算方法。

2.1 流量疊加法

當水庫周邊有水文站,其控制流域面積占壩址以上流域面積比重較大且水文資料較完整可靠時,可由上游來水、區間徑流以及庫面產流疊加得到入庫流量,即流量疊加法。以丹江口水庫為例,丹江口水庫的入庫流量由3部分組成,分別是水庫回水末端附近漢江干流白河站及黃龍灘、荊紫關、西峽、西坪、賈家坊等水文站以上流域產生的洪水,干支流水文站以下到水庫周邊以上區間陸面產生的徑流以及水庫庫面產流。實際制作入庫流量預報時,為了準確將上游來水演算至入庫點,通常采用錯時疊加干支流入庫控制站實況及預報流量,并演算匯入區間來水,計算得到入庫流量的預報過程。具體計算公式如下:

Q丹江口入庫,t=Q白河,t-9+Q黃龍灘,t-3+Q賈家坊,t-6+Q荊紫關,t-3+Q西峽,t-3+Q西坪,t-6+Q區間,t(1)

式中:Q丹江口入庫,t為t時刻丹江口水庫的入庫流量計算值;Q區間,t為t時刻白河水文站至丹江口水庫壩址整個無控區間的流量過程;Q白河,t-9為t-9時刻白河水文站流量,其中白河水文站流量演進至丹江口水庫入庫點的平均傳播時間為9 h;其他符號依次類推。

計算流程如圖2所示。

2.2 庫容反推法

庫容反推法是實際水庫報汛工作中的常用方法,其基本原理是水量平衡方程,通過監測所得水庫下泄流量大小及同時段庫水位變化,結合水位-庫容關系曲線反推入庫流量。根據水庫庫容特性不同,可分為靜庫容反推法和槽蓄量反推法。

(1) 靜庫容反推法。靜庫容反推法是采用壩前某水位站水位值為代表,查算靜庫容水位-庫容曲線,得到計算時段初和時段末水位對應的水庫庫容,繼而得出該時段對應的水庫庫容變化量;計算該時段水庫總出庫水量,根據水量平衡原理,得到水庫在該時段內的總來水量;總來水量與時段相除,即得到計算時段內的平均入庫流量。具體計算公式如下:

式中:Q—入庫為時段平均入庫流量;Q—出庫為時段平均出庫流量;ΔV為時段始末庫容差;Δt為時段長;Z為時段內損失量,包括蒸發、滲漏等,一般較小時忽略不計。

靜庫容反推法是最常用的入庫流量反推計算方法,其所得入庫流量成果通常受水庫庫水位代表性及反推計算時段長的影響較大。

(2) 槽蓄量反推法。槽蓄量反推法與靜庫容反推法均基于水量平衡原理,不同在于其采用槽蓄量的變化量反映時段始末水庫水面線下的水庫蓄量差,以反映動庫容的影響。具體計算公式為

式中:Wt2為時段末時刻的水庫實際槽蓄量,Wt1為時段開始時刻的水庫實際槽蓄量。

河道槽蓄量是指河道上、下斷面間水面線以下的蓄水量。以丹江口水庫為例,若把丹江口水庫比作一個長河道,那么丹江口水庫槽蓄量就是從水庫回水末端斷面到壩前段斷面的水面線線下蓄水量(見圖3)。槽蓄量的計算一般采用斷面法,該方法使用固定的觀測斷面,計算精度取決于實測數據的精度和斷面密度,斷面越密,則計算精度越高,但計算工作越繁瑣。

隨著水文監測行業技術的不斷發展及斷面資料的日積月累,河道地形觀測數據精度日益提升,獲取難度下降,根據地形圖計算所得河道槽蓄量愈加精確。但實際中,因無法獲得連續的水面線過程,一般只能根據水庫的水位站布設將水庫的庫容進行離散分段,根據各分段的代表水位和分段水位庫容曲線計算槽蓄量,再根據上述公式計算獲得入庫流量,故通常也將此類入庫流量計算方法稱為分段庫容反推法。

丹江口水庫由漢江、丹江兩個并聯式水庫組成,兩庫庫容相差不大,入庫徑流量以漢江為主。丹江庫區為典型湖泊型水庫,漢江庫區在庫尾存在一定動庫容。相對丹江口水庫整體庫容而言,漢江庫區楔蓄庫容占比較小,動庫容影響極小。靜庫容曲線能夠準確反映水庫水位-蓄量關系,因此,本文采用靜庫容反推法計算丹江口水庫入庫流量。

3 靜庫容反推法影響因素分析

根據靜庫容反推法原理概述可知,該法計算所得入庫流量受計算時段長、滑動平均算法及庫水位代表性影響較為明顯。本節將從這3方面影響因素入手,通過對比分析,提出一種適用于丹江口水庫的入庫流量平滑修正計算方法,使得其計算成果平滑穩定,計算所得入庫洪峰、洪量等特征值具有更好的代表性,以更準確地描述丹江口水庫來水過程。

3.1 不同時段反推入庫流量

采用水庫運行管理單位漢江水利水電(集團)有限責任公司水調中心提供的基礎數據,以2017,2019,2021年秋汛期水庫出庫流量及庫水位資料系列(時段長為1 h)為代表,分別取1,2,3,6,12 h等不同時段長反推計算入庫流量,所得入庫過程對比如圖4所示。由圖4可見,不同時段長計算所得入庫流量過程總體呈鋸齒狀“跳動”,有時甚至出現負值,明顯與真實入庫過程不符。此外,不同反推計算時段長對反推過程的跳動性影響明顯,其中較小計算時段(1,2,3 h)反推流量過程跳動頻繁且劇烈,6 h時段長反推流量過程跳動幅度明顯小于3 h反推過程,12 h時段長反推的入庫流量過程較6 h反推過程又進一步平滑。隨著計算時段長度的增加,反推流量過程逐漸平滑,跳動幅度逐步減小,但其反推過程對入庫洪峰、洪量等特征值的捕捉效果逐漸下降。

3.2 不同滑動平均算法反推入庫流量

為了減弱或消除干擾因素的影響,進而提高入庫流量曲線的光滑度,需要對入庫流量計算結果進行平滑處理。平滑的原則是一方面要消除數據中存在的干擾成分,另一方面要保持原始數據的曲線特性不變。

一般常用且平滑效果較好的方法為五點三次平滑算法及滑動平均算法。五點三次平滑算法對流量修正具有較好的效果,計算入庫流量曲線通常既能消除數據中的干擾成分,又能保持原有曲線特性不變,但該方法需要未來兩個時段的流量數據,無法真正應用于實際報汛工作。因此,本節主要采用滑動平均算法進行計算?;瑒悠骄ㄓ嬎愎饺缦拢篞mt=W1Q1+W2Q2+…+WnQt-n+1(4)

式中:Qmt為t時段入庫流量修正值,m3/s;Qt為t時段入庫流量,Qt-1為t-1時段入庫流量;Wn為各時段入庫流量的權重,n=1,2,…,N,N為入庫流量數據的總個數。當W=1/N的時候,各個時段的權重相等,即為N個時段的算術平均,稱為簡單滑動平均算法;否則,可通過對不同的時段設置各不相同的權重,用以區分不同時段的重要性,稱為加權滑動平均算法。通常來說,距離t時段越近,權重越大。本節采用簡單滑動平均算法,以“20190916”洪水為例,計算時段步長取 1 h,分別選取不同的平滑步長(3,6,12 h)來分析對入庫流量的修正效果,并與報汛入庫流量過程進行對比分析(見圖5)。

由圖5分析可知,3 h的平滑步長計算入庫數據明顯波動較大,且洪峰較報汛明顯偏大;6 h的平滑步長計算入庫過程平滑度、洪峰及峰現時間與報汛入庫過程基本一致,且比3 h過程波動明顯減少;12 h平滑效果更好,但洪峰較報汛明顯偏小,峰現時間更提前。綜上,建議采用6 h進行時段平滑處理,能夠更有效地消除數據波動誤差,且保證入庫流量的合理性。

3.3 壩前水位代表性

壩前水位代表性是影響反推入庫流量過程形態的另一要素。目前丹江口水庫僅采用龍王廟站作為壩前水位代表站,本節分析在原有基礎上,增加考慮漢江庫區龍口站及丹江庫區窯溝站、巡路口站,共計4站的水位算術平均值作為丹江口水庫水位,以提高其代表性。根據水位-庫容曲線和出庫流量計算入庫流量,并與龍王廟單站反推結果進行對比分析。選取“20190916”和“20210929”兩場洪水過程,采用多站反推入庫流量與報汛入庫流量進行對比分析,如圖6~7所示。

由圖6~7可知,采用多站平均水位作為庫水位,反推計算所得入庫流量過程、整體峰型、峰值與報汛基本一致。由表1可知,多站推算較報汛入庫整體方差較小,所得入庫過程劇齒波動變幅明顯小于報汛數據,平滑性更好。

4 實例分析

4.1 動、靜庫容影響分析

對于湖泊型水庫,庫區水面近似為水平線,楔形庫容相對較小,壩前水位可以代表庫區水位,采用靜庫容可較好地代表水庫蓄水量。2013年丹江口大壩加高工程完成,正常蓄水位170.00 m,漢江庫區回水長約180 km,丹江庫區回水長約100 km,形成了漢江、丹江兩個并聯式水庫。兩庫庫容相差不大,但漢江來水量大于丹江,兩者水面線呈現不同的形式。漢江庫區屬山區河道型水庫,為典型的緩坡壅水曲線;丹江庫區屬湖泊型水庫,壅水曲線基本屬平直線狀。丹江口水庫來水以漢江為主,靜庫容以上相應的楔形庫容主要集中于漢江庫區。整體來講,丹江口水庫漢江庫區尾部回水形成的楔型水體相對其整個庫區巨大庫容來說較小,基于不同壩前水位和入庫流量工況下的漢江庫區動庫容大小如表2所列。

根據上述計算分析可知,丹江口水庫動庫容水量及所占靜庫容比例較小,故可以不考慮動庫容影響,采用靜庫容曲線、出庫流量及庫水位代表值查算蓄水體體積,進而反推入庫流量過程。根據前述分析,靜庫容反推入庫流量受計算時段長、滑動平均算法及庫水位代表性影響明顯,本節從這三方面入手,對比不同方法組合的優越性,提出一種較適用于丹江口水庫的入庫流量平滑修正計算方法。

4.2 兩種推求方法的比較分析

基于丹江口水庫2021年典型入庫過程對應的白河、黃龍灘等入庫控制站報汛流量及區間雨量資料系列,采用API-UH模型[16](即降雨產流采用降雨徑流相關圖,流域匯流采用單位線,兩者結合形成完整的降雨徑流模型)計算白河、黃龍灘等站至丹江口水庫壩址無控區間來水過程??紤]傳播時間,將上述控制站流量進行合成后疊加區間流量得到預報入庫流量,并與報汛入庫流量、平滑入庫流量進行比較,預報入庫流量與報汛入庫、平滑入庫的過程比較見圖8。由圖8可以看出,預報入庫流量與報汛入庫流量過程十分吻合,多數洪峰吻合較好,部分洪峰因區間來水較大而略有差別;高水時,報汛入庫流量的跳動性較小,低水時(流量量級在3 000 m3/s以下),報汛入庫流量較預報入庫流量跳動性較大。

統計分析2011~2021年丹江口水庫較大洪水過程對應預報入庫、平滑入庫流量(為靜庫容反推法計算入庫流量并經6 h平滑后的數據,以下簡稱“平滑入庫”)及報汛入庫的洪峰及過程洪量等特征(見表3),可見,三者的入庫洪峰相差不大,19場洪水平均相對誤差分別為2.9%和5.4%,最大相對誤差分別為5.8%和12.7%,其中90%的場次洪水,入庫流量洪峰相對誤差在5%以內。對于峰現時間,預報入庫與報汛入庫、平滑入庫的平均偏差為2.7 h和2.8 h,所有場次洪水誤差均在6 h以內。綜上可知,預報入庫流量過程和平滑入庫流量過程基本一致。

由于預報入庫流量依賴于入庫站來水的資料及區間洪水的估算,當區間來水比例較大時,對區間洪水的估算直接影響入庫洪水的預報精度。“20190806”洪水主要受丹江口庫區強降雨影響,強降雨集中在丹江口庫區無控區域附近,區間匯流速度快,API模型對這種集中強降雨的估算明顯與實況有較大偏離,致使預報入庫流量與報汛入庫流量出現了明顯偏離(見圖9)。在峰后,對區間來水進行分割,分析區間洪峰約4 330 m3/s,而模型計算的區間洪峰為6 320 m3/s,加之白河等控制站合成演算后疊加錯時的選擇,導致該場洪水預報入庫流量洪峰與平滑入庫流量洪峰出現了明顯偏差。

根據上述分析,丹江口水庫的入庫流量可采用靜庫容反推法并經6 h平滑計算獲得。當區間強降雨占比較大、顯著影響入庫過程時,則需要更加關注區間產匯流計算模型的選擇與校正。

5 結 論

本文針對南水北調中線水源工程丹江口水庫開展了入庫流量推求方法研究,闡述了流量疊加法和庫容反推法兩類常用方法的基本原理和流程,并針對靜庫容反推法的3類主要影響因素進行了定量對比分析,確定了適用于丹江口水庫的入庫流量平滑修正計算方法,在保障入庫過程準確性的同時,獲得盡可能平滑的流量過程,并以丹江口水庫實例洪水過程進行了方法驗證。結果表明:耦合6 h反推、6 h滑動、多站水位算數平均值為壩前水位可獲得丹江口水庫更科學合理的入庫過程,在保證洪峰流量和峰現時間更接近真實情況的同時,減小過程方差、提高平滑度。所提方法能夠兼顧丹江口水庫入庫流量的準確性和平滑性,為水庫的科學精細調度提供重要技術支撐。此外,考慮丹江口水庫入庫流量在主汛期與非主汛期差異性較大,針對不同時期來水特點采用差異化反推方法可能會進一步提高入庫流量計算準確性,這也是進一步研究的重點方向。

參考文獻:

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

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

[3]趙凱華,武新宇,李保健,等.考慮梯級蓄泄的水庫預報入庫流量修正方法[J].水力發電學報,2014,33(2):73-81.

[4]王世策,胡曉勇.大型水庫計算入庫流量波動過大問題分析[J].安徽水利水電職業技術學院學報,2010,10(3):19-21.

[5]唐海華,陳森林,趙云發,等.三峽水庫入庫流量計算方法研究[J].中國農村水利水電,2008(4):26-28.

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

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

[8]董前進,傅建彬,陳森林.基于統計圖形的入庫流量預報誤差分布規律[J].水電能源科學,2011,29(4):5-7.

[9]閔要武,王俊,陳力.三峽水庫入庫流量計算及調洪演算方法探討[J].人民長江,2011,42(6):49-52.

[10]楊金標,舒凱,張后來,等.水庫水量平衡計算中消除水位跳變影響的過濾算法[J].人民長江,2019,50(8):98-102.

[11]DENG C,LIU P,GUO S L,et al.Estimation of nonfluctuating reservoir inflow from water level observations using methods based on flow continuity[J].Journal of Hydrology,2015,529:1198-1210.

[12]ZHANG Z L,SUN H,SU Y.Water use efficiency and its influencingfactors in arid areas of northwest China[J].Journal of Ecology and Rural Environment,2017,33(11):961-967.

[13]方崇惠,郭生練,段亞輝,等.還原入庫洪水的一種簡便新方法[J].巖土工程學報,2008,30(11):1743-1747.

[14]李榮容,吳國堯.水庫入庫流量抗差修正研究[J].中國農村水利水電,2008(11):12-14.

[15]張愛靜,姚文鋒,吳智健.丹江口水庫入庫徑流變化特征分析[J].人民長江,2020,51(3):81-86,93.

[16]鐘小燕,文磊,余鐘波.連續API模型在沂河臨沂站徑流預報中的應用[J].人民長江,2017,48(13):26-30.

(編輯:謝玲嫻)

Study on smooth correction method of inflow into Danjiangkou Reservoir

ZENG Fanlin1,NIU Wenjing2,ZHANG Chengxiao3,YAN Fangjia2,XING Wenhui2

(1.Hanjiang Water Conservancy and Hydropower (Group) Limited Liability Company,Wuhan 430048,China; 2.Bureau of Hydrology,Changjiang Water Resources Commission,Wuhan 430010,China; 3.Hanjiang Hydrology and Water Resources Survey Bureau,Hydrology Bureau of Changjiang Water Resources Commission,Xiangyang 441022,China)

Abstract: Accurate inflow into reservoir calculation is an important basis for scientific reservoir operation.However,the severe jumping phenomenon of reverse deduced inflow remains a difficult problem.In this study,the principles and applicability of commonly used methods for calculating inflow are discussed.Focusing on the characteristics of Danjiangkou Reservoir capacity,we quantitatively compared and analyzed the main factors that affect the inverted calculation of static storage capacity by using the single variable methods,such as the number of reverse calculation periods,the sliding average algorithm,and representativeness of water level in front of the dam.Then,a smooth correction calculation method for the inflow into Danjiangkou Reservoir was proposed.The results showed that for Danjiangkou Reservoir,the inverted deduced reservoir capacity method can obtain both quasi-deterministic and smooth inflow calculation results,coupled with 6 h inversion,6 h sliding,and the multi-station average levels as the water level in front of the dam.The proposed method ensured that the peak flow and peak time are closer to the real situation,while the variance of the inflow flow was smaller and the smoothness was higher.The research results can provide important technical support for the scientific and refined operation of Danjiangkou Reservoir.

Key words: inflow into reservoir;hydrologic balance;reverse deduction by static storage capacity;moving average;Danjiangkou Reservoir

收稿日期:2023-03-22;接受日期:2023-11-13

基金項目:國家重點研發計劃青年科學家項目(2023YFC3210500);湖北省自然科學基金聯合基金項目(2023AFD203);漢江集團科研項目“現狀人類活動條件下的漢江流域暴雨洪水特征及規律研究”(H202112)

作者簡介:曾凡林,男,正高級工程師,主要從事水庫調度工作。E-mail:11059994@qq.com

通信作者:牛文靜,女,高級工程師,博士,主要從事水文水資源預報與水庫調度工作。E-mail:dgniuwenjing@163.com

主站蜘蛛池模板: 欧美在线国产| 亚洲一级毛片| 91成人在线观看视频| 老司国产精品视频91| 2021国产精品自拍| A级毛片高清免费视频就| 欧美不卡在线视频| 国产熟女一级毛片| 成人国产免费| 亚洲综合九九| 天堂成人av| 免费看久久精品99| 国产亚洲精品资源在线26u| 国产男人天堂| 免费黄色国产视频| 波多野结衣视频一区二区| 91黄视频在线观看| 国产精品综合久久久| 国产精品不卡永久免费| 亚洲日韩精品伊甸| 91亚洲影院| 999国产精品| 91伊人国产| 免费jjzz在在线播放国产| 国产精品久久精品| 欧美在线网| 97se亚洲综合| 国产一区免费在线观看| 97久久精品人人| 免费A∨中文乱码专区| 中文字幕在线欧美| 亚洲无码电影| 热伊人99re久久精品最新地| 亚洲愉拍一区二区精品| 日本不卡免费高清视频| 一本久道热中字伊人| 天天综合亚洲| 欧美三级日韩三级| 国产欧美又粗又猛又爽老| 福利在线不卡| 国产丝袜一区二区三区视频免下载| 女人天堂av免费| 亚洲a级在线观看| 国产精品久久久久久久久久98| 美美女高清毛片视频免费观看| 老司机精品一区在线视频| 夜精品a一区二区三区| 国产微拍一区二区三区四区| 伊人91在线| 欧美精品影院| 亚洲伦理一区二区| 91啦中文字幕| 日韩美一区二区| 激情六月丁香婷婷| 午夜激情婷婷| 丰满少妇αⅴ无码区| 国产精品成人一区二区| 国产尤物在线播放| 久草视频精品| 99久久精品国产麻豆婷婷| 九九热这里只有国产精品| 欧美日韩激情在线| 亚洲乱码在线播放| 丁香亚洲综合五月天婷婷| 亚洲无线一二三四区男男| 尤物成AV人片在线观看| 久久这里只有精品23| 日韩欧美视频第一区在线观看| 亚洲自拍另类| 亚洲成A人V欧美综合天堂| 国产第一色| 99精品影院| 久久国产精品嫖妓| 99免费视频观看| 国产成人高清亚洲一区久久| 一级毛片免费高清视频| 亚洲第一国产综合| 精品久久人人爽人人玩人人妻| 欧美色视频网站| 婷婷六月综合网| 无码在线激情片| 国产精品对白刺激|