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

基于HSPF的東江分布式水文模型構建

2011-01-25 06:47:04董延軍鄧家泉鄭江麗馬志鵬
長江科學院院報 2011年9期
關鍵詞:模型

董延軍,鄧家泉,李 杰,鄭江麗,馬志鵬

(珠江水利委員會珠江水利科學研究院,廣州 510611)

基于HSPF的東江分布式水文模型構建

董延軍,鄧家泉,李 杰,鄭江麗,馬志鵬

(珠江水利委員會珠江水利科學研究院,廣州 510611)

針對東江流域,在基于HSPF水文模型的基礎上,構建了東江流域分布式水文模型。詳細介紹了HSPF分布式水文模型構建的原理、模型數據準備、模型構建以及模型驗證過程。從氣象分布變異的角度出發,將東江流域劃分成51個子流域,從物理分布變異的角度出發,將東江流域的土地利用類型與土壤類型和坡向重新疊置成29種特征類型的下墊面。通過模型的模擬效果來看,兩校準站(河源站和博羅站)的相對誤差均小于15%,NASH系數都在0.9以上,模擬效果良好。

HSPF模型;分布式水文模型;東江流域;水文響應

1 概述

東江是珠江流域3大水系之一,流域總面積為27 040 km2,其中廣東省境內23 540 km2,占流域總面積的87.06%。東江流域是廣州、惠州、東莞、深圳等地市的重要供水水源地,同時還擔負著對香港特別行政區的供水任務,總供水人口達3 000余萬人。東江水問題十分敏感,這是因為:①東江下游三角洲地區包括香港在內不但經濟地位十分重要,而且其政治地位也極為敏感。根據2008年統計資料,珠江三角洲4市(深圳、惠州、東莞和廣州)GDP總和為27 015.75億元,占廣東省GDP的58.87%,占全國GDP的11.37%;②隨著東江流域人口快速增長以及城市化進程速度的不斷加快,水資源的量和質與水資源快速增長形成的需求產生了尖銳的矛盾,這關系著東江流域特別是下游三角洲地區的經濟社會可持續發展。因此在可預見的未來,東江流域水資源問題研究是一個熱點問題。

目前在國內應用SWAT模型研究水文模擬和非點源模擬方面已有大量的應用成果問世[1-3],但是應用HSPF模型的成果卻并不多見,HSPF模型還不為大眾所熟知。根據國外特別是美國水文水質模擬的應用實踐來看,HSPF模型與SWAT模型一樣,也是應用十分廣泛的流域水文水質模擬工具,模擬精度絲毫不遜于SWAT模型,二者各有千秋,各有側重[4,5]。因此本文研究出于2方面的考慮,其一是引入HSPF模型,以東江流域為研究對象,建立基于HSPF的分布式水文模型;其二是為觀察分析東江水問題提供一個新的工具和視角。

2 基于HSPF的東江分布式水文模型的構建原理與方法

HSPF(Hydrology Simulation Program Fortran)模型是能長時間模擬水量與水質的半分布式水文模型。HSPF模型的介紹、結構、原理參見文獻[6,7]。HSPF模型起源于斯坦福模型Ⅳ,最早主要用于模擬水文過程,后來不斷衍生發展,增加了許多模擬水質的模塊,逐漸演變成了能夠模擬流域水文水質重要的非點源模擬工具。今天所看到的HSPF模型已經以WinHSPF的形式內嵌于BASINS系統。BASINS系統是一套基于GIS技術的整合式平臺,內嵌MapWindows GIS,WinHSPF,WDMUtil,GenScn,系統能夠實現地形、地貌、土地利用/覆被、土壤、流域等數據的自動生成和疊加處理,并擁有強大的前后處理功能,為快速、方便地建立分布式水文模型提供了條件。

下面重點介紹基于HSPF的東江分布式水文模型構建的原理與過程。

2.1 流域的空間異質性與空間離散

分布式模型相較集總式模型來說,一個顯著的優越性就是考慮流域空間的異質性。流域空間異質性是通過空間離散的方式來實現[8]。對于HSPF模型來說,空間離散是通過氣象分布變異和物理特性變異2種方式來進行劃分,形成不同類型的水文響應單元,如透水面積上的PERLND101、PERLND102、PERLND103,不透水面積上的IMPLND101,每個水文響應單元對應不同的水文響應過程。

氣象分布變異的處理是根據流域內氣象站點的分布情況,按照Thiessen多邊形的原理為各個子流域分配氣象數據。物理特性分布變異主要是根據流域內土地利用分布情況,并參照土壤分布、坡向分布進行劃分,由于BASINS4.0系統只識別土地利用LUCC參數,因此在考慮流域空間分布變異方面,一般事先需要在ArcGIS或者ArcView工具中將土地利用分布、土壤分布、坡向分布進行疊置(UNION)分析處理,然后再以土地利用類型的方式導入疊置后的數據。

東江流域共有5個氣象站點,分別為連平站、尋烏站、河源站、惠州站、深圳站,Thiessen多邊形分配各子流域情況見圖1(a)。其中連平站控制11個子流域,編碼分別為6,12~15,17~19,21~23;尋烏站控制9個子流域,編碼分別為2~5,7~11;河源站控制17個子流域,編碼分別為16,20,24,25,27~38,40;惠州站控制13個子流域,編碼分別為1,39,41~50,52;深圳站控制1個子流域,編碼為51。

土利用類型與土壤類型和坡向疊置后重新分類(Reclassify),一共生成29種特征類型的下墊面,這29種新地土地利用類型就是用原土地利用類型、坡向圖和土壤類型圖在GIS平臺中疊加(Overlay)生成圖1(b)。

2.2 單元水文響應過程分析

相對于透水面積的水文響應過程來說,不透水面積的水文響應過程處理就非常容易,因此本部分重點介紹透水面積的水文響應過程原理。HSPF模型的水文過程是概念型的,水流的運動分成垂向運動和橫向運動,在垂向上模型分成5層,即植被截流層、上土壤層、下土壤層、淺層地下層、深層地下層;在水平方向上,根據模型設計的一系列規則,自上而下分別產生坡面漫流、壤中流、地下徑流3種徑流成份。

2.2.1 冠層截流模型

冠層截留模型在垂向上水量平衡關系表達為

式中:Wc2,Wc1分別表示計算時段初和時段末植被冠層截留蓄積量(mm);P(t)表示計算時段內降雨量或其它降雨補充量(mm)。

落地雨量計算公式為

圖1 東江流域空間離散Fig.1 Spatial discretization of Dongjiang watershed

式中:NP(t)表示準落地雨量(mm);IM表示冠層截流容量(mm);E截為截留層的蒸發量。

冠層部分沒有徑流發生,因此沒有水平方向上的水量分配。

2.2.2 上土壤層水量傳輸與分配

上土壤層垂向上水量平衡關系表達為

式中:UZS1(t),UZS2(t)為計算時段初與時段末的上土壤層蓄積量(mm);ΔD為地表滯蓄增量(mm);Pr(%)為ΔD進入上土壤層的百分比,具體計算公式見式(4);PERC為滯后下滲量;E上為上土壤層的蒸發量。

式中UZS,UZSN分別為上土壤層蓄積量和上土壤層額定容量。

由公式(3)可知一部分水量已進入上土壤層,那么另一部分則形成坡面滯蓄增量和壤中流滯蓄增量,進而產生水平方向上的坡面出流和壤中出流。水量平衡關系分別為式(5)和式(6):

式中:SURS1(t),SURS2(t)表示時段初和時段末的坡面滯蓄量;SURO表示坡面出流。

式中:SRGX1(t),SRGX2(t)表示時段初和時段末的壤中流滯蓄量;INTF表示壤中出流;ΔSRGX為壤中流滯蓄增量。

2.2.3 下土壤層水量傳輸與分配

下土壤層不產生出流,因此只有垂向方面的傳輸,沒有水平方向的分配。在垂向上下土壤層的水量平衡關系為

式中:LZS1(t),LZS2(t)為計算時段初與計算時段末的下土壤層蓄積量(mm);IND為直接下滲量;PERC含義同式(3);E下為下土壤層的蒸發量;Pg(%)為進入地下水蓄積的下滲百分量,計算公式見式(8)。

當LZS/LZSN≤1時,

式中LZS,LZSN分別為時段下土壤層蓄積量和下土壤層額定容量。

2.2.4 淺地下層水量傳輸與分配

淺層地下水蓄積量通過下土壤層的滲漏而得到補充,其值為Pg(%)·(IND+PERC),并通過地下水出流GWF、地下水蒸發E地和繼續向深層地下水部分滲漏的量而失掉水分。其垂向上水量平衡關系為

式中:SGW1(t),SGW2(t)表示時段初與時段末的淺層地下水蓄積量;K24L表示進入地下水蓄積的下滲量分配到深層地下水的份額(無量綱);E地為地下水蒸發量。

水平方向水量分配表現為產生地下出流:

式中:KGW表示地下水蓄泄系數,KVARY表示地下水退水率。

2.2.5 深地下層水量傳輸與分配

HSPF模型對深層地下水蓄積的設計是只有入流,沒有出流。深層的入流量為

式中K24L是進入地下水蓄積的下滲量分配到深層地下水的份額,無量綱,以百分數表示。

2.3 RCHRES匯流模型

早期的斯坦福模型流域匯流部分是運用單位線法和馬斯京根法原理進行匯流計算,HSPF模型對匯流部分進行了改進和調整,構建RCHRES(Reach&Reservoir)模塊。RCHRES模塊假定水體為單向流動,將地表水體運動視作線性波(Linear Wave),模擬河道中地表水體的非恒定流運動[9]。RCHRES模塊不僅實現了流域匯流計算,而且能夠實現水動力學計算、水體溫度擴散、污染物遷移轉化等功能,使得HSPF模型不但能夠模擬水量,而且也能夠模擬水質,大大拓寬了HSPF模型功能和使用范圍,使HSPF模型有一個飛躍式的發展。

2.3.1 RCHRES模型概化

RCHRES模型將河段或水庫概化成圖2所示。水體的入流成份是降雨,出流部分包括灌溉用水、發電用水、其它抽取的水量(Pumping)、溢流等。

圖2 RCHRES概化圖Fig.2 Generalization of RCHRES

HSPF模型約定所有的進入RCHRES的入流需要通過唯一的入口INFLO,從入口INFLO進入的水量稱之為IVOL。RCHRES的出流可以通過若干個出口,譬如第N個出流定義為OFLO(N),相應出口通過的水量稱之為OVOL(N)。出口通過的水量之和稱之為ROVOL。

2.3.2 RCHRES模型原理

徑流演算模型的基本方程為

式中:VOL為計算時段末的RCHRES水量;VOLS為時段初的RCHRES水量;PRSUPY為河段水體表面的降雨量;VOLEV河段水體表面的蒸發量。當水體量足夠時,采用線性關系的表達式,總出流量ROVOL可表達為

式中:KS是權重因子;ROS是計算起始時段的出流量;COKS=1.0-KS;ROD是計算末時段的待求出流量;dT指時間隔;KS的范圍一般在0≤KS≤1,一般不要超過0.5。聯立(12)和(13)有

式中:VOLT=IVOL+PRSUPY-VOLEV+VOLS。

式(14)中有2個未知量,一個是VOL,另一個是ROD,因此需要另外增加一個函數關系表達式來完成上述方程的求解。如果不考慮時間因素,每個出口的流量與體積的函數關系表示為

聯立式(14)與(15)、(16)進行求解,求解過程如圖3。

圖3 RCHRES模塊求解示意圖Fig.3 Solution of RCHRES module

2.4 模型數據準備

對于任何一種分布式水文模型來說,一般都需要龐大的流域數據量來支持,甚至還可能會涉及到數據的合成和分解,數據準備工作量是相當大的。BASINS系統提供了一個功能非常完備的前處理工具WDMUtil,為分布式水文模型構建提供了極大的方便。根據HSPF模型的要求,WDMUTIL需要導入降雨資料、蒸發資料,其它的資料如風速資料(Wind)、溫度(ATEM)、露點(DEWP)、太陽輻射(SOLR)等信息,雖然本次模型不需要,但從模型構建的角度和完整性來說,需要虛擬構建,即值設為“0”的時間序列。據此依次將東江流域內的5個氣象站點的數據導入WDMUtil。每個站點不同類型的數據對應不同的數據源(DSN),也就是FORTRAN語言對應的通道號。界面顯示如圖4。

圖4 水文氣象數據導入WDMUtil工具Fig.4 WDMUtil tool for importing hydrometeorological data

2.5 基于HSPF的東江分布式模型構建過程

2.5.1 流域網絡拓撲空間建立

目前流域網絡拓撲空間建立普遍采用D8算法。D8算法是利用DEM網格,按最陡坡度方向決定地表匯流方向、生成河道網的方法。這些成果被ARC/INFO及ARCVIEW等GIS軟件采用,也被BASINS系統采用,為流域水文模擬帶來了極大方便。在BASINS4.0里,這部分功能集成在流域自動劃分工具(Automatic Delineation)當中。流域網絡拓撲空間將各子流域與東江水系按照流域上下游屬性有機地統一起來,見圖5。

圖5 東江流域拓撲空間建立Fig.5 Topological space of Dongjiang watershed

圖6 東江流域HSPF模型拓撲結構圖Fig.6 Topological structure of HSPF Model of Dongjiang watershed

2.5.2 流域屬性構建

流域屬性主要是通過BASINS GIS的環境,提取出流域水系的河道水深、寬度、坡度等地形基礎資料,然后通過曼寧公式計算FTABLEs表。這就是所謂標準的計算方法(Standard Method)。標準方法認定河道糙率n的值為0.05,并認定河道的斷面形狀是復合梯形斷面形式。梯形斷面的主槽部分邊坡為1∶1,河漫灘部分邊坡為0.5∶1。如果有具體的河道地形資料,可以用標準方法對其進行修正。

2.5.3 BASINS向WinHSPF跳轉

經過基礎數據準備,土地利用類型數據、WDM水文氣象時間序列和BASINS系統空間屬性前處理工作,HSPF模型運行前提條件已經齊備。跳轉后生成的HSPF模型再經流域氣象站點分隔處理后,就建立起一個基于HSPF的東江流域分布式水文模型,見圖6。圖6中RCHRES1為東江流域出口河段,即博羅段,RCHRES32為河源站控制段。這兩個控制段作為本次研究的水文校準站。

3 模型驗證

3.1 模型的校準

HSPF模型率定的基本步驟是:首先運行UCI文件,然后觀察模型輸出的模擬結果,最后根據實際觀測結果,對模型的相關參數反復進行調整,使得模型模擬結果滿足精度要求。參數調整過程采用以下4個步驟進行調整:①根據模型手冊以及美國應用HSPF模型的40個流域使用的參數數據庫HSPFPARM軟件給定所選參數的初始值;②調整參數值,實現年水量總平衡,其關鍵參數主要包括LZSN和UZSN;③徑流成份劃分調整(地表徑流、壤中流、基流)。其關鍵性參數主要包括INFILT,LZETP,KVARY和AGWRC。INFILT控制著有降雨轉化成地表徑流和地下徑流的比例。KVARY和AGWRC控制著地下退水過程;④水位(流量)曲線形狀調整。其關鍵參數包括INTFW和IRC。水位(流量)形狀受地表徑流和壤中流量雙重影響,INTFW和IRC控制著進入河道的產流量和時間。詳細的模型校準方法和參數校準見文獻[6]。參數校準見表1。

3.2 模擬效果評價

選用相對誤差(RE)和NASH-Sut-cliff標準來評估模型在校準和驗證過程中的模擬效果。NASH公式計算如下:

式中:Qi為第i時刻觀測流量(m3/s);Q'i為第i時刻模擬流量(m3/s);n為時段總數;ˉQ為平均觀測流量。NASH系數是描述計算值對目標值擬合精度的無量綱統計參數,一般取值在0~1之間,并認為0.9以上為甲等,0.7~0.9為乙等,0.5~0.7為丙等,低于0.5則認為是偏差過大,結果不適宜。

表1 不同性質地面參數取值Table 1 Parameters of different land types

模擬效果要求模擬值與實測值年均誤差RE小于實測值的15%,NASH系數為乙等以上。

研究運用東江河源站和博羅站觀測的天然徑流量與模擬值進行對比。設置校準期為1956年1月1日至1965年12月31日,驗證期1為1966年1月1日至1969年12月31日,驗證期2為1970年1月1至1973年12月31日。

河源水文站月徑流過程模擬效果見表2,河源水文站月徑流過程模擬結果與實測結果的對比見圖7;博羅水文站日徑流過程模擬效果見表3,博羅水文站日徑流過程模擬結果與實測結果的對比見圖8。從表1和表2的結果來看,2個校準站的NASH系數都在0.9以上,屬于甲等,本研究所采用的HSPF模型在月徑流過程模擬中取得了較高的模擬精度,在校準期前后的2個驗證期均得到良好驗證。

表2 河源站月徑流過程模擬效果Table 2 Simulation results of monthly runoff at Heyuan station

圖7 河源站月徑流過程模擬結果與實測結果的對比Fig.7 Comparison of simulated and measured monthly runoff curves at Heyuan station for two calibration periods

表3 博羅站月徑流過程模擬效果Table 3 Simulation results of monthly runoff at Boluo station

圖8 博羅站月徑流過程模擬結果與實測結果的對比Fig.8 Comparison of simulated and measured monthly runoff curves at Boluo station for two calibration periods

4 結語

本文根據東江流域的特點,在基于美國環保署(EPA)開發的HSPF模型基礎上,構建了東江流域分布式水文模型,對HSPF模型的原理、分布式水文模型的數據準備、構建以及驗證等關鍵性環節進行了有效地探討。從氣象分布變異的角度出發,將東江流域劃分成51個子流域,從物理分布變異的角度出發,將東江流域的土地利用類型與土壤類型和坡向重新疊置成29種特征類型的下墊面。通過模型的模擬效果來看,兩校準站(河源站和博羅站)的相對誤差均小于15%,NASH系數都在0.9以上,模擬效果良好,這充分證明基于HSPF的分布式水文模型可以適用于東江流域的水文模擬,為東江流域水循環和水資源問題研究提供了一個新的研究方法。

[1] 代俊峰,崔遠來.基于SWAT的灌區分布式水文模型—Ⅰ.模型構建的原理與方法[J].水利學報,2009,2 (40):145-152.(DAI Jun-feng,CUI Yuan-lai.Distributed Hydrological Model for Irrigation Area Based on SWAT—Ⅰ.Principle and Method[J].Journal of Hydraulic Engineering,2009,2(40):145-152.(in Chinese))

[2] 黃清華,張萬昌.SWAT分布式水文模型在黑河干流山區流域的改進與應用[J].南京林業大學學報,2004,3 (28):22-26.(HUANG Qing-hua,ZHANG Wan-chang.Improvement and Application of GIS-Based Distributed SWAT Hydrological Modeling on High Altitude,Cold,Semi-arid Catchment of Heihe River Basin[J].Journal of Nanjing Forestry University,2004,3(28):22-26.(in Chinese))

[3] 張雪松,郝芳華,楊志峰,等.基于SWAT模型的中尺度流域產流產沙模擬研究[J].水土保持研究,2003,12 (10):38-42.(ZHANG Xue-song,HAO Fang-hua,YANG Zhi-feng,et al.Runoff and Sediment Yield Modeling in Meso-Scale Watershed Based on SWAT Model[J].Research of Soil and Water Conservation,2003,12(10): 38-42.(in Chinese))

[4] JASWINDER S,KNAPP H V,DEMISSIE M.Hydrologic Modeling of the Iroquois River Watershed Using HSPF and SWAT[J].Journal of the American Water Resources Association,2005,41(2):343-360.

[5] IM S,BRANNAN K,MOSTAGHIMI S,et al.A Comparison of SWAT and HSPF Models for Simulating Hydrologic and Water Quality Responses from an Urbanizing Watershed[C]∥The Society for Engineering in Agricultural,food,and Biological Systems.Presentation at the ASAE Annual International Meeting.Nevada,July 27-30,2003.

[6] Hydrological Simulation Program-Fortran(HSPF)USER’s Manual[K].California:Aqua Terra Consultants,2001.

[7] 董延軍,李 杰,鄭江麗,等.流域水文水質模擬軟件(HSPF)應用指南[M].鄭州:黃河水利出版社,2009.(DONG Yan-jun,LI Jie,ZHENG Jiang-li,et al.The HSPF Software Manual for Watershed Simulation on Hydrology and Water Quality[M].Zhengzhou:Yellow River Water Conservancy Press,2009.(in Chinese))

[8] 劉昌明,鄭紅星,王中根.流域水循環分布式研究[M].鄭州:黃河水利出版社,2006.(LIU Chang-ming,ZHENG Hong-xing, WANG Zhong-gen. Distributed Hydrological Cycle of River Basins[M].Zhengzhou: Yellow River Conservancy Press,2006.(in Chinese))

[9] 芮孝芳.水文學原理[M].北京:中國水利水電出版社,2004.(RUI Xiao-fang.Hydrology Theory[M].Beijing:China Water Power Press,2004.(in Chinese))

Distributed Hydrological Model for Dongjiang Watershed Based on HSPF Model

DONG Yan-jun,DENG Jia-quan,LI Jie,ZHENG Jiang-li,MA Zhi-peng
(Pearl River Hydraulic Research Institute,Guangzhou 510611,China)

A distributed hydrological model is established for Dongjiang watershed based on Hydrological Simulation Program-Fortran(HSPF)Model.The theory of the construction,the data preparation,the model construction and model calibration are presented in details.The watershed is divided into 51 subfields from the perspective of spatial heterogeneity;while in terms of physical heterogeneity,the land use type,soil type and slope aspect are overlaid into 29 types of underlying surfaces.The relative errors at two calibrated stations(Heyuan station and Boluo station)are less than 15%and the NASH coefficients are larger than 0.9,which indicates a favorable simulation effect.

HSPF model;distributed hydrological model;Dongjiang watershed;hydrologic response

P334

A

1001-5485(2011)09-0057-07

2010-12-03

國家東江水專項(2008ZX07211-10-02)

董延軍(1972-),男,山西原平人,高級工程師,主要從事水資源規劃與管理的研究,(電話)13719292058(電子信箱)d_yj@163.com。

(編輯:王 慰)

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 日韩欧美国产精品| 亚洲一区二区三区中文字幕5566| 欧美日韩国产在线人| 亚洲第一色网站| 久久一级电影| 国产打屁股免费区网站| av一区二区三区在线观看| 波多野结衣无码中文字幕在线观看一区二区 | 伊人久久福利中文字幕| 中文字幕佐山爱一区二区免费| 久久www视频| 亚洲成肉网| 色婷婷综合激情视频免费看| 黄色网页在线播放| 日韩高清在线观看不卡一区二区| 色老二精品视频在线观看| 在线观看国产黄色| 美女免费精品高清毛片在线视| 欧美日韩午夜| 亚洲看片网| 国产激情影院| 国语少妇高潮| 高清不卡毛片| 国产97色在线| 久青草免费在线视频| 日韩东京热无码人妻| 精品综合久久久久久97超人该| 国产主播在线一区| 国产手机在线观看| 亚洲福利片无码最新在线播放 | 国产玖玖玖精品视频| 2021国产乱人伦在线播放| 欧美亚洲一区二区三区导航| 久久久91人妻无码精品蜜桃HD| 国产丝袜啪啪| 中文字幕亚洲综久久2021| 一区二区三区在线不卡免费| 中文字幕在线永久在线视频2020| 99热最新在线| 真人高潮娇喘嗯啊在线观看| 国内精品免费| 国产精品无码影视久久久久久久| 国产高清色视频免费看的网址| 欧美一级爱操视频| 国产精品久久久免费视频| 国产精品lululu在线观看| yjizz视频最新网站在线| 老色鬼久久亚洲AV综合| 熟女视频91| 先锋资源久久| 国产精品99久久久| 91成人精品视频| 久久无码免费束人妻| 国产凹凸一区在线观看视频| 精品国产www| 国产国模一区二区三区四区| 国产中文一区a级毛片视频| 欧类av怡春院| 极品尤物av美乳在线观看| 九九热精品视频在线| 精品一区二区三区四区五区| 国产欧美综合在线观看第七页| 国产熟女一级毛片| 欧美成人aⅴ| 国产福利在线观看精品| 午夜无码一区二区三区| 国模视频一区二区| 色综合综合网| 91青青草视频| 久久这里只有精品免费| 亚洲成在线观看 | 伊人久久影视| 久久久久久国产精品mv| 一级成人a毛片免费播放| 99这里精品| 亚洲日韩精品伊甸| 天堂久久久久久中文字幕| 五月激情综合网| 国产日本一区二区三区| av色爱 天堂网| 91在线国内在线播放老师| 久久99久久无码毛片一区二区|