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

基于SWMM的城市社區尺度管網排水模擬
——以福州市某排水小區為例

2022-11-15 11:11:40葉陳雷徐宗學雷曉輝陳陽李鵬王京晶
南水北調與水利科技 2022年2期
關鍵詞:模型

葉陳雷,徐宗學,雷曉輝,陳陽,李鵬,王京晶

(1.北京師范大學水科學研究院 城市水循環與海綿城市技術北京市重點實驗室,北京 100875;2.中國水利水電科學研究院水資源研究所,北京 100038;3.東北大學資源與土木工程學院,沈陽 110819)

近年來,我國經濟社會快速發展,城市面積不斷擴大,城市人口持續增長,城市水問題隨之而來。城市“看海”層出不窮,逢雨必澇逐漸演變為我國很多大中城市的痼疾。氣候變化與城市化對城市水文循環規律產生了深遠的影響[1],“城市雨島”與“城市熱島”效應被廣泛關注,城市內澇已成為危害居民財產安全的重要問題[2-5]。城市化使下墊面對降雨的響應機制發生改變,不透水面積增多,植被覆蓋顯著減少,徑流量增加[6-7]。

城市雨洪模擬方法主要包括水文學方法、水動力學方法以及水文水動力結合的方法[8-9]。自然情況下,城市中的降雨經過初期損失后形成凈雨,再經過坡面匯流進入城市管網中,最終由管網系統排出。在這一過程中,水流在城市地表的產流、匯流以及在管網中的匯流過程極為重要,是研究城市洪澇過程的重要工作。國內外相關研究人員[10-12]已針對城市雨洪模擬開展了大量研究。對于城市降雨徑流的計算,我國自主研發的雨洪模型較少,主要集中在對徑流量的計算上。岑國平[13]通過對城市雨水的地表產流、地表匯流、管網匯流過程進行研究,建立了城市雨水徑流計算模型。周玉文等[14]將徑流分為地表徑流和管網徑流兩個階段,構建了城市雨水徑流模型。目前應用廣泛的城市雨洪模擬軟件包括美國環境保護署開發的SWMM(storm water management model)[15]、丹麥DHI開發的MIKE URBAN[16]以及英國Wallingford開發的InfoWorks ICM[17]。相比于其他模型軟件,SWMM由于其開源的特性,被廣大從事科學研究的人員應用在城市雨洪模擬中。石小芳等[18]基于SWMM對寧波地區構建了臺風和非臺風暴雨情景內澇響應模型,對比分析了兩種不同暴雨情景下內澇的特征及差異。戎貴文等[19]基于SWMM中的低影響開發(low impact development,LID)模塊以淮南某一工業園區為研究對象,分析了多種LID組合措施對水質水量的影響。目前,雖然關于SWMM的應用研究已較為豐富,但對于城市社區尺度下管網排水全過程的分析尚待進一步完善。鑒于SWMM開源的特點,本文基于實測場次降雨水位數據,采用優化算法對模型參數進行自動率定,并對各水文參數進行局部敏感性分析。在此基礎上,選取不同重現期與不同降雨歷時的組合工況對該管網排水系統進行計算。通過數據準備、模型構建、參數自動率定及敏感性分析、組合情景模擬及管網排水分析等主要步驟,以我國沿海地區福州市晉安區某典型排水小區為例,提出一套較為完整的城市排水內澇分析體系。

1 研究區數據

福建省會福州是典型的沿海易澇城市,位于閩江下游地區,屬于亞熱帶季風氣候。近年來,福州頻繁受臺風暴雨影響,強降雨天氣嚴重影響人們的日常生活。如:2005年的“龍王”臺風對福州市造成了極大破壞,城區主要交通道路幾乎全部被淹;2016年的“鲇魚”臺風登陸福州時,12 h降雨量達到195 mm,經濟損失14.2億元[20-22]。選擇福州市晉安區某典型居民區為研究區域,小區西面為茶園路,北面為晉安河支流,東面和南面以居民區街道與附近其他居民區相隔開,總面積約11.48 hm2。小區中有社區、學校、廣場等典型居民區建筑設施。建模時,將房屋、廣場、道路等概化為不透水面,將綠地等概化為透水面。圖1(a)為該居民區位置示意圖。利用SWMM對其排水系統進行評估,使用數據主要包括小區排水管網設施、土地類型以及小區地形。研究小區建模數據包括小區內的管網數據、地形數據和土地類型數據。管網數據包括檢查井、管道、排口3個組成部分,其中:檢查井數據包括檢查井的位置、井底高程、井深;管道數據包括各管段的位置、形狀、尺寸;排口信息包括排口位置等。

圖1 研究區及SWMM概化圖

本文使用2 m分辨率DEM數據計算子匯水區平均坡度,管網數據及地形數據來自于福州市城區水系聯排聯調中心,土地類型數據通過目視解譯得到。概化的小區SWMM見圖1,共計排口6個,檢查井78個,管段78根。圖1(b)為SWMM的概化示意圖,本文使用的部分建模數據見表1。

表1 SWMM需要的基礎管網數據

2 主要方法

2.1 SWMM

SWMM是一個分布式水文模型,將子匯水區作為基本水文響應單元,產匯流過程在每個基本單元上進行計算。根據不同的城市地表特性,每個基本單元被概化為透水區和不透水區兩部分。在不透水區中,根據地表是否能夠洼蓄劃分為有洼蓄不透水區和無洼蓄不透水區兩部分。基于研究區地形數據,利用ArcGIS表面分析計算百分比坡度,通過分區統計可得各子匯水區的平均坡度為2.09%。該模型[23]提供了下滲、坡面匯流及管網匯流等主要的水文水動力計算。

SWMM的下滲模塊提供了Horton模型、Green-Ampt、Curve Number等不同的下滲模型,本文選擇廣泛使用的Horton模型進行計算。Horton模型可以較好地描述下滲率隨降雨動態變化過程的規律,其計算公式為

ft=fc+(f0-fc)e-kt

(1)

式中:ft為土壤t時刻下滲率,mm/h;f0為初期下滲率,mm/h;fc為穩定下滲率,mm/h;t為時間,h;k為衰減指數,1/h。

在計算坡面匯流時,SWMM將每一個子匯水區概化為一個非線性水庫,將最大洼地蓄水量作為水庫的容量,根據水量平衡方程計算每一個水庫中的水深,通過曼寧公式計算得到地表出流量為

(2)

式中:Q為出流量,m3/s;W為子流域寬度,m;S為平均坡度;n為地表曼寧系數;d為水深,m;dp為洼蓄深度,m。

城市雨洪過程中,經下滲損失后形成的凈雨由坡面匯流進入管網系統。在管網匯流階段,SWMM提供了恒定流、運動波、動力波3種計算方式。其中:恒定流假設管網中任意時刻水流運動要素保持不變,是對實際水流運動最為簡化的處理方式;運動波通過聯立求解簡化的連續方程和動量方程計算水量,但是其無法有效模擬復雜逆坡、回水、環狀管網等復雜的實際情形;動力波通過求解完全的圣維南方程組,可以對復雜流態進行模擬。因此,本文采用動力波對管網水流進行模擬。

(3)

(4)

式中:Q為流量,m3/s;H為水深,m;A為過水斷面面積,m2;g為重力加速度,9.8 m/s2;Sf為摩阻比降;t為時間,s;x為距離,m。

2.2 遺傳算法

遺傳算法(genetic algorithm,GA)是一種通過模擬自然選擇和進化過程來搜索最優解集的算法,常被使用在最優化問題的研究中。由于SWMM開源的特性,可以在遺傳算法優化參數過程中,將SWMM的模擬值與實測數據計算得到的評價指標作為遺傳算法種群個體中的適應度函數,進而完成遺傳算法的更新迭代。遺傳算法主要包括以下3個基本的遺傳算子,分別為交叉操作、變異操作和選擇操作。基于遺傳算法優化SWMM參數的主要步驟見圖2。

圖2 遺傳算法優化SWMM參數步驟

初始種群的初始化。SWMM需率定的參數構成一個解空間,設置進化代數計數器以及最大進化代數,并隨機生成M個個體作為初始群體P(0)。

個體適應度計算。計算群體P(t)中各個個體的適應度。

將交叉算子作用于群體,交叉算子是遺傳算法中起核心作用的算子。

將變異算子作用于群體。即是對群體中的個體串的某些基因座上的基因值作變動。群體P(t)經過選擇、交叉、變異運算之后得到下一代群體P(t+1)。

計算目標函數。將更新得到的種群個體作為需要率定的模型參數,SWMM通過讀取其輸入文件(.inp)進行計算,(.inp)是一個結構化的文本格式,模型參數在文本中位于對應的模塊下,通過Python語言對其進行參數的讀取和寫入。在參數替換為GA更新后的個體數值后,即產生新的(.inp)文件,即更新后的模型驅動數據,利用第三方庫Pyswmm可以方便地調用模型計算引擎,以獲取SWMM計算結果,提取目標點位水位或流量過程作為模擬值,通過模擬值與實測值計算目標函數。該步驟是通過遺傳算法率定SWMM參數的核心。

選擇操作。將選擇算子作用于群體。選擇的目的是把優化的個體直接遺傳到下一代或通過配對交叉產生新的個體再遺傳到下一代。選擇操作是建立在群體中個體的適應度評估基礎上的。

終止條件判斷。在每輪迭代更新種群后,計算納什效率系數或判斷是否已滿足迭代次數要求,如滿足條件則停止遺傳算法迭代過程,以算法進化過程中的最大適應度個體作為最優解輸出,并終止循環。

2.3 Morris敏感性分析

本文選用Morris篩選法分析參數的敏感性[24]。Morris篩選法是從模型參數中選出某一個變量,在其取值范圍內進行變動,使自變量以固定步長百分率改變,并以此驅動模型,從而得到不同xi對應的目標函數y的結果。最終敏感性判別因子取多個Morris系數的平均值,其計算公式為

(5)

式中:y*為改變參數后的模型輸出;y為改變參數前的模型輸出;Δxi為參數的變化。本文針對通過遺傳算法率定的參數,選取一定范圍改變參數值,并固定其他參數,以研究區主要排水口的洪峰流量為衡量指標,分析各參數的靈敏性。敏感性系數的計算公式為

(6)

式中:Y0為基準輸出結果;Yi為參數第i次變化的輸出;Pi為參數相對于基準輸入的百分比變化量;n為參數的改變次數;S為參數的敏感性判別因子。當|S|≥1時,表明參數高度敏感;當0.2≤|S|<1時,表明參數中高等敏感;當0.05≤|S|<0.2時,表明參數中低等敏感;當|S|<0.05時,表明參數不敏感。

3 結果分析

3.1 參數優化及其敏感性

在構建研究小區SWMM的過程中,綜合利用泰森多邊形以及研究區的衛星影像圖對其進行子流域劃分。計算時每個子流域可以看作一個集總式水文模型,且包括下滲、地表匯流以及管網匯流3個主要的物理過程。因此,SWMM降雨徑流模擬中主要有以下3個模塊的參數,包括:子匯水區面積(hm2)、不透水百分比(%)、特征寬度(m)、平均坡度百分比(%)、不透水地表曼寧系數、透水地表曼寧系數、不透水地表洼地蓄水深度(mm)、透水地表洼地蓄水深度(mm)等子匯水區模塊參數;最大入滲速率(mm/h)、最小入滲速率(mm/h)、衰減指數(1/h)等下滲模塊參數;糙率等管道匯流模塊參數。其中:子匯水區面積(hm2)、特征寬度(m)等參數由實際地形數據提取得到,不需要進行人為率定;而其他參數如糙率、最大入滲速率(mm/h)、最小入滲速率(mm/h)等雖然具有明確的物理意義,但是一般需要通過實測數據與模擬數據的對比來進行率定。由于其物理意義明確,這類參數通常具有一定的取值范圍,而其中一些參數對模型模擬結果影響比較大,稱之為敏感參數。準確率定模型參數,使模型能夠更好地體現研究區實際特征,是建模過程中十分重要的步驟。在實際工作中,往往有大量的監測數據,手動率定往往帶有極大的主觀性,且工作量很大。與手動調參相比,遺傳算法能夠實現模型參數的自動尋優,為模型率定工作提供了便利。本文在遺傳算法中,設置種群參數為30,迭代次數為200,其余超參數的突變概率設置為0.1,交叉概率為0.5,適應度函數使用納什效率系數(Nash-Sutcliffe coefficient of efficiency,ENS),其計算表達式為

(7)

本文通過參考文獻中給出的模型參數范圍,以匯水區及管道參數的上、下閾值作為遺傳算法種群個體的上下邊界。參考文獻[25]與文獻[26],分別設置SWMM中子匯水區模塊、下滲模塊、管道模塊中的參數,見表2。經過200次迭代,遺傳算法計算得到的參數值已趨于收斂,最終得到的參數值如下所述。子匯水區模塊參數中,不透水面曼寧系數為0.022,透水面曼寧系數為0.226,不透水區洼地蓄水深度為0.659 mm,透水區洼地蓄水深度為5.998 mm;下滲模塊參數中,最大下滲率為83.62 mm/h,最小下滲率為1.534 mm/h,衰減系數為3.502,干燥時間為5 d;管道模塊參數中,糙率為0.01。

表2 SWMM關鍵參數的取值范圍

率定期選取20180531場次降雨,該場次總降雨量為43.8 mm,降雨歷時為200 min,計算得到納什效率系數為0.82。實測最大水位為4.64 m,模型計算得到的水位為4.66 m。實測峰現時間為當日20:40,模型計算得到峰現時間為20:35。將模型率定得到的參數代入SWMM,并用驗證場次的降雨數據驅動模型,將模擬結果與驗證場次的實測數據進行比較。驗證期選取20180620場次降雨,該場次總降雨量為53.1 mm,降雨歷時為240 min,計算得到納什效率系數為0.81。實測最大水位為4.750 m,模型計算得到的水位為4.758 m。實測峰現時間為當日18:30,模型計算得到峰現時間為18:25。因此,模型模擬的效果較好,可以用于實際洪澇情景計算。兩次峰值的模擬水位和實測水位誤差均在2%以內,模型的擬合度較高。率定參數后的SWMM在兩場降雨實測水位與模擬水位的對比見圖3。

圖3 模擬水位與實測水位

由于小區實測資料比較缺乏,無法使用更多的場次降水資料進行模型率定和驗證。較少的實測數據較易導致模型的適用性降低。因此,為進一步深化對模型參數的認識,有必要對模型主要參數的靈敏性進行一定的分析,找到敏感性強的參數,在模型的后續應用中重點率定這些參數。對于不敏感的參數,結合參數的物理意義及其取值區間進行適當調整。利用修正的Morris法,將率定后的參數作為基準參數值。為分析主要參數的靈敏性,在原參數基礎上設置-30%、-20%、-10%、10%、20%、30%的變動,驅動模型計算得到多組輸出結果。選擇排水口的洪峰流量和小區平均徑流系數為衡量指標,分析和討論各參數的敏感性。

圖4為采用Morris篩選法得到的洪峰流量對模型參數敏感性系數S分布,在計算過程中發現,下滲模塊的參數干燥時間在變動中對結果的影響為0,這里不在圖中作出。從圖4可以看出,隨著重現期的變化,各參數敏感性有所變化,其中:在重現期為1 a時,各參數敏感性系數分布較為接近,均位于[0.05,0.20)內,屬于中等敏感參數;在重現期為20 a 時,最大下滲率位于[0.2,1.0)內,屬于敏感參數,其他參數仍為中等敏感,但判別因子S有所增大;在重現期為50 a時,衰減系數與最小下滲率位于[0.2,1.0)內,屬于敏感參數,其他參數為中等敏感。在該小區模型中,最大下滲率、最小下滲率以及衰減系數的敏感性系數隨雨強的變化更為明顯。

圖4 Morris篩選法洪峰流量對參數的敏感性系數S分布

3.2 SWMM情景模擬

以3種設計降雨驅動模型,設計雨型來自福建省《城市及部分縣城暴雨強度公式》[27],共設置16種降雨方案。其中:設置4種降雨重現期,分別為1 a一遇、10 a一遇、20 a一遇、50 a一遇;同時設置4套降雨歷時方案,分別為2、6、12、24 h。所有設計暴雨方案見表3。20 a一遇的3種降雨歷時下的設計降雨過程線見圖5。

表3 設計降雨方案

圖5 20 a一遇的3種設計降雨過程線

設置4種降雨重現期以及4種降雨歷時,通過SWMM計算得到在16種情景下小區排口的洪峰流量。小區內共有2個主要的排水管網系統,因此分別統計了2個排口處的流量過程線,見圖6。圖6中的C17為連接排口A的管道,C74為連接排口B的管道,分別用以統計模型計算的過流量。

通過計算可知,模型較為準確地模擬了小區的降雨徑流過程。降雨歷時為2 h時:當重現期為1 a時,排口A的洪峰流量為0.39 m3/s,峰現時間為1:00;當重現期為10 a時,排口A的洪峰流量為0.59 m3/s,峰現時間為0:55;當重現期為20 a,排口A的洪峰流量為0.67 m3/s,峰現時間為0:55;當重現期為50 a時,排口A的洪峰流量為0.75 m3/s,峰現時間為0:55。對于排口B以及其他重現期的降雨規律類似,隨著降雨重現期的增加,該小區的洪峰流量逐漸增加,峰現時間較為接近。從圖6還可以看出:當降雨重現期為20 a時,2 h降雨歷時下排口B的洪峰流量為0.32 m3/s,峰現時間為0:55,而此時的降雨峰值出現時間為0:50;6 h降雨歷時下排口B的洪峰流量為0.32 m3/s,峰現時間為2:30,而此時的降雨峰值出現時間為2:25;12 h降雨歷時下排口B的洪峰流量為0.34 m3/s,峰現時間為4:55,而此時的降雨峰值出現時間為4:50;24 h降雨歷時下排口B的洪峰流量為0.32 m3/s,峰現時間為9:40,而此時的降雨峰值出現時間為9:35。可以看到在不同的降雨歷時下,小區排口峰值流量比較接近,而峰現時間直接受雨型的影響,一般出現在降雨峰值后的5 min。為進一步分析小區管網排水特征,統計各情景方案下小區排口的洪峰流量,見表4。相比于降雨歷時,小區洪峰流量受降雨重現期或雨強的影響較大。

表4 不同內澇情景下的洪峰流量(排口B)

圖6 不同情景下排口流量過程模擬

3.3 管網排水能力分析

為進一步分析小區管網排水特征,在4種降雨重現期與4種降雨歷時的組合下驅動SWMM,并將模型計算時間統一設置為48 h,統計各情景方案下管網排水特征,見表5。在不同重現期下,均有部分節點和管道發生超載。節點超載個數和管道超載個數達52個以上。對于同一降雨歷時的降雨,隨著降雨重現期的增大,由于降雨強度變大,同一時間內城市下墊面收集的雨水更多,匯入排水管網的水量更大,這使得節點超載數、節點溢流數、滿流管道數以及滿管時間均呈現增加的趨勢。但節點平均超載時長卻呈現減少的趨勢,是由于在降雨重現期增大導致超載的節點隨之增多所致。當降雨重現期不變,降雨歷時變大時,節點平均超載時長、節點超載數、滿流管道數基本保持不變,而節點溢流數減少,滿管時間增加。其中,在更短的降雨歷時下,同一重現期降雨在時程上的分配更為密集,當管道排水能力不足時,水從管道兩側進入檢查井,此時更容易造成檢查井溢流,實際中對于短歷時強降雨要引起足夠重視。總的來說,各項指標受降雨重現期的影響較大,而節點超載數、滿流管道數受降雨歷時影響相對較小。

表5 不同內澇情景下管網排水特征統計

4 結 論

以福州市某典型居民小區為例,基于實測的管道液位數據,利用遺傳算法對SWMM參數進行優化率定,并利用Morris法對影響產匯流計算的主要參數敏感性進行了分析。在此基礎上,利用不同重現期的降雨驅動模型,比較分析各種情景下小區排水管網的排水能力。主要結論如下:

基于遺傳算法率定研究區SWMM參數,20180531場次降雨計算得到納什效率系數為0.82;實測最大水位為4.64 m,模型計算得到的水位為4.66 m;實測峰現時間為當日20:40,模型計算得到峰現時間為20:35。20180620場次降雨計算得到納什效率系數為0.81;實測最大水位為4.750 m,模型計算得到的水位為4.758 m;實測峰現時間為當日18:30,模型計算得到峰現時間為18:25。基于Morris篩選法進行SWMM參數敏感性分析,得到洪峰流量對模型參數敏感性系數S分布圖。結果顯示,降雨重現期較小時,參數敏感性系數分布相對均勻。隨著降雨重現期增大,參數敏感性分布相對集中。衰減系數、最大下滲率、最小下滲率以及管道曼寧系數比其他參數更為敏感。

對研究小區進行管網排水計算,得到不同降雨情景下的洪峰流量、節點平均超載時長、節點超載數、節點溢流數、滿流管道數、滿管時間等指標。結果表明,在降雨歷時相同時,隨著降雨重現期的增加,該小區的洪峰流量逐漸增加,峰現時間較為接近。在不同的降雨歷時下,小區排口峰值流量比較接近,而峰現時間直接受雨鋒的影響,一般出現在降雨峰值后的5 min。在降雨歷時增加時,滿管時間相應增加。各項指標受降雨重現期的影響較大,其中,節點超載數、滿流管道數受降雨歷時影響相對較小。

本文在城市社區尺度下完成排水模型構建、參數率定及敏感性分析、組合工況模擬及排水能力量化分析等工作,較為完整地分析了城市排水內澇情勢,對相關防洪排澇工作有一定的借鑒意義。論文工作本身也存在以下不足:實測數據的缺乏使本文僅使用2場實測雨洪過程進行參數的率定和驗證,使模型參數在其他場次降雨中的泛化能力有所欠缺,未來應搜集更多的實測資料,不斷地對模型進行優化,增加模型的魯棒性和泛化能力。此外,更多的實測資料也有利于進一步對遺傳算法的超參數進行校核;本文僅對社區尺度洪澇過程中地表產流、地表匯流及管網匯流3個過程做了模擬,并未涉及河網匯流及地表漫流等重要過程,如直接導致城市內澇節點的溢水過程。在后續工作中需要將SWMM與其他模型工具相結合,綜合利用水文與水動力學方法,模擬城市尺度洪澇全過程。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲综合在线最大成人| 91精品啪在线观看国产91| 一区二区午夜| 免费Aⅴ片在线观看蜜芽Tⅴ | 亚瑟天堂久久一区二区影院| 国产精品短篇二区| 成年片色大黄全免费网站久久| 欧美精品v欧洲精品| 亚欧成人无码AV在线播放| 国产成人久久777777| 日本a级免费| 人妻精品全国免费视频| 精品少妇人妻av无码久久| 国产电话自拍伊人| 亚洲av成人无码网站在线观看| 免费视频在线2021入口| 性视频一区| 毛片网站在线看| 国产成人一区免费观看| 欧美激情视频一区二区三区免费| 欧美高清国产| 亚洲欧美人成电影在线观看| 国产成人免费| 亚洲成在线观看 | 久久综合丝袜日本网| 国产靠逼视频| 欧美亚洲一二三区| 亚洲综合18p| 亚洲精品成人福利在线电影| 小蝌蚪亚洲精品国产| 91尤物国产尤物福利在线| 综合亚洲网| 亚洲日韩AV无码一区二区三区人| 全部毛片免费看| 性做久久久久久久免费看| 97人妻精品专区久久久久| 色天天综合久久久久综合片| 伊人色天堂| 国产精品偷伦在线观看| 日韩视频免费| 国产拍在线| 亚洲伊人电影| 午夜激情福利视频| 国产精品一老牛影视频| 亚洲成人www| 久久久久久久久18禁秘| 亚洲精品视频网| 一级福利视频| 91探花国产综合在线精品| 国产精品女在线观看| 综合天天色| 欧美激情第一欧美在线| 欧美精品在线视频观看| 又黄又湿又爽的视频| 老司机久久99久久精品播放| 一区二区三区成人| 天天综合网色| 久久熟女AV| 99re经典视频在线| 国产精品网址你懂的| 免费在线不卡视频| 九色综合视频网| 亚洲最大福利视频网| 91福利片| 日韩精品无码不卡无码| 亚洲中文字幕无码爆乳| 丁香婷婷在线视频| 亚洲综合狠狠| 在线国产91| 国产一区二区免费播放| 无码不卡的中文字幕视频| 成人免费一区二区三区| 欧美中文字幕在线视频| 国产丝袜第一页| 99视频在线观看免费| 久热中文字幕在线| 日韩精品无码免费专网站| 亚洲综合欧美在线一区在线播放| 国产一级裸网站| 欧美日韩久久综合| 亚洲天堂久久新| 国产精品久久久久久久久|