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

集總式喀斯特水文模型構建及其應用

2017-03-27 02:35:27許波劉董增川
水資源保護 2017年2期
關鍵詞:模型

許波劉,董增川,洪 嫻

(河海大學水文水資源學院,江蘇 南京 210098)

集總式喀斯特水文模型構建及其應用

許波劉,董增川,洪 嫻

(河海大學水文水資源學院,江蘇 南京 210098)

針對大部分喀斯特流域屬于無資料地區的現狀,構建資料要求較低的集總式喀斯特水文模型(lumped karst hydrological model,LKHM)。LKHM以三水源新安江模型為基礎,針對喀斯特流域二元三維的流域性質,在匯流計算層做了改進。將巖溶地貌概化為3類:溶溝、溶隙、管道,采用滯后演算法模擬溶溝對地表徑流的調蓄作用,采用若干線性水庫串聯模擬不同大小口徑的溶隙對地下徑流的調蓄作用,將所有管道概化為一個大管道模擬地表徑流與地下徑流的交換;最后在貴州平湖流域進行模型檢驗和應用。結果表明:LKHM模擬精度較高,說明模型的結構、參數合理,具備一定的實用價值。

集總式喀斯特水文模型;巖溶地貌;溶溝;管道;溶隙;線性水庫

喀斯特流域是指存在地下水和地表水對可溶性巖石進行破壞和改造形成的地貌的流域。與非喀斯特流域相比,喀斯特流域有著以下獨特的地理特征和水文特征[1]:①流域不閉合,相鄰流域存在相當數量的地表水與地下水交換;②流域內地下水面以上多垂直節理,孔道、溶隙、落水洞、漏斗等地貌發育較好;③流域植被以及地表土壤發育較差,裸露的巖石多,覆蓋層薄;④流域內存在較大的地下水蓄水庫,蓄水容量與垂直方向節理發育的大小和數量有關;⑤流域存在非巖溶地貌區。

喀斯特流域的上述特征,一方面使大量地表水匯入地下,導致流域內嚴重的缺水狀況;另一方面又由于流域內排水不暢,引發不少洪澇災害。研究應用在喀斯特流域的水文模型是解決這些問題行之有效的辦法。

國內外學者對此進行了廣泛的研究。國外學者首先相繼提出了巖溶水雙重、三重和四重介質模型。由于喀斯特流域地下水系統介質結構及水動力機理的復雜性,從20世紀70年代開始,一些學者以系統論為指導展開了集總式水文模型的研究,其中以嚴啟坤[2]、袁道先等[3]為代表。隨著計算機性能的提高與3S技術的發展,很多學者開始將分布式、半分布式水文模型引入到喀斯特流域,常用的包括SWAT模型、TOPMODEL與SWMM。Schomberg等[4-6]在改進的SWAT模型基礎上建立了巖溶流域分布式水文模型。SWAT模型具有輸入數據易獲取、計算效率高的優點,但是主要用于長時段的模擬,難以用于喀斯特流域短時段的洪水過程模擬。索立濤等[7]在子流域中采用TOPMODEL作產流計算,解決了在巖溶地區無法使用的問題。潘歡迎[8]提出了改進的KARST-TOPMODEL模型。TOPMODEL的優點在于所需要的參數少且容易獲得,缺點是對地下水文過程描述過于簡單,且不能用DEM來生成喀斯特流域的天坑和暗河。Campbell等[9-11]對SWMM進行了相關研究。SWMM的優點在于能夠較真實地反映洪水波在管道中的運移過程,缺點是對管道的幾何特征較敏感,僅適合于一些具有詳細勘測資料的管道系統的水文模擬。

總體而言,分布式水文模型考慮了氣象、下墊面等影響因子的差異性,能夠較為準確地描述喀斯特流域內的水文過程,但要求研究區提供大量的數據,這制約了其在喀斯特流域的應用。集總式水文模型以系統為研究對象,雖然部分參數物理意義不太明確,且往往具有多解性,但由于其結構簡單,對數據要求相對較低,因此,在大部分屬于無資料地區的喀斯特流域水文模擬中具有不可忽視的優勢。本文針對巖溶地貌的特點,建立了精度較高、資料要求較低的集總式水文模型(lumped karst hydrological model,LKHM)。

1 模型結構

特定流域的集總式模型結構決定于該流域的性質。喀斯特流域的典型流域性質為由不同類型的含水介質構成的二元三維空間系統。二元指流域有地表、地下2個水系,它們的邊界一般不重合,通過水力聯系構成一個緊密的整體。地表水系主要由溶洼、溶溝、漏斗、石芽、峰叢及溶水洞等組成,地表河流一般作為洪水時期暴雨的轉移通道;地下水系主要由溶洞、管道、溶隙、伏流及地下河等組成,地下河流是一般時期中小徑流的轉移通道。地表水系是徑流的形成場和分配場,地下水系是徑流的調蓄場和輸移場,2個水系存在比較大的徑流交換[12]。三維指喀斯特流域的能量、物質遷移和轉化是一個從上游到下游、從地表到地下的三維立體過程,各處地貌性質不一,對徑流的調蓄作用相差很大,這是區別于非喀斯特流域的平面特征。因此,對集總式模型結構的要求有:①必須反映二元水系的徑流場性質對徑流最終形成的影響;②必須反映三維立體空間內差異很大的各類地貌調蓄作用。

LKHM基于三水源新安江模型結構,前3層與新安江模型一致[13]。產流計算層采用蓄滿產流方法,是因為在喀斯特流域地表滯水作用比較弱,降雨直接進入地下水系中,控制產流量的是降雨總量而不是降雨強度。匯流計算是模型的核心,基于流域的二元三維性質和巖溶地貌的特點,LKHM重新構建了匯流計算層。

喀斯特流域巖溶地貌廣為發育,調蓄作用非常強,模擬出這些巖溶地貌的調蓄作用并將其反映在匯流計算層是本模型的出發點。作為集總式水文模型,將整個流域看成一個整體,采用概化研究的方法,將流域整體內的多種巖溶地貌概化為3類。

a. 溶溝。包括石芽,峰叢等地貌,發育在地面,主要影響地表徑流的匯流,調蓄作用較弱,采用滯后演算法進行地表水匯流計算,即可模擬這一類地貌對地表徑流的推移坦化作用,計算公式為

(1)

(2)

(3)

式中:Q(t)為t時刻河段出流量,m3/s;I(t)為t時刻河段入流量,m3/s;τ為滯后時間,h;Δt為計算時段長,h;K為出流系數;C0、C1為系數。

b. 管道。由發育最為強烈的石灰巖形成,包括落水洞、巖溶漏斗等地貌,位于地表與地下之間。管道是地表水系與地下水系進行水量交換的主要途徑,是地下河形成的主要介質。地表徑流直接通過管道匯入地下徑流接受地下的調蓄,基本不經過地表的調蓄作用。模型中將所有管道的過流能力概化為一個大管道的過流能力,即單位時間所有管道通過的水量等于該大管道通過的水量。大管道過流能力用參數cf表示,cf用來分離出匯入地下水的地表徑流,其作用類似于二水源劃分結構中的穩定下滲率fc,即時段內地表水徑流量小于cf的,地表水全部匯入地下水中,作為地下水進行匯流計算;時段內地表水徑流量大于cf的,匯入地下水的徑流量為cf,其余部分作為地表水進行匯流計算,計算公式為

當Rs0

(4)

當Rs0>cf時,

(5)

式中:Rs0為二水源劃分的地表水徑流量,mm;Rg0為二水源劃分的地下水徑流量,mm;Rs為經過管道過流后剩余的地表水徑流量,mm;Rg為經過管道過流后的地下水徑流量,mm;cf為大管道過流能力,mm。

c. 溶隙。由發育較差的白云巖組成,位于地面以下。溶隙發育非常密集,在整個流域地下立體空間內都存在著,壤中流和地下徑流在匯流的過程中,流經溶隙時會受到很強的調蓄作用,且不同大小的溶隙調蓄作用不同。繼續采用概化分析的方法,將溶隙的大小進行概化,按地面以下接近地面位置的溶隙口徑的大小概化為5種:極大型、大型、中型、小型、極小型。而對于某一種特定的溶隙,由于應力的作用,它的口徑從地面往地下是越來越小的,到地下最底層位置時溶隙都變成極小型。向下入滲的徑流,由于溶隙口徑的減小,在不同大小的溶隙交界面上,一部分徑流將排出地下立體空間,一部分將繼續下滲。繼續下滲的徑流,在下一個不同大小的溶隙交界面上,繼續重復一部分徑流排出、一部分徑流下滲的過程。因壤中流與地下徑流運動過程相似,故將壤中流并入地下徑流處理。地下徑流排出前所經過溶隙的大小、數量不同,受到的調蓄作用也不同。

假設5種溶隙占巖溶地區的面積比例分別為X1、X2、X3、X4、X5,地下徑流從接近地面的整個平面均勻下滲,則每種大小的溶隙分得的地下徑流比例也為X1、X2、X3、X4、X5。從X1面積上向下入滲的地下徑流,首先經極大型溶隙的調蓄后,在極大型溶隙與大型溶隙的交界面上,一部分徑流將排出地下立體空間,假設這部分比例(以下簡稱排出比例)為Y1,剩下的(1-Y1)比例的地下徑流繼續下滲,經大型溶隙的調蓄后,在大型溶隙與中型溶隙的交界面上,一部分徑流排出,假設這部分占繼續下滲的徑流比例為Y2,以此類推,假設在中型溶隙與小型溶隙交界面上的排出比例為Y3,小型溶隙與極小型溶隙交界面上排出比例為Y4,最后剩下的地下徑流經極小型溶隙調蓄后排出整個地下溶隙調蓄系統。從X2面積上向下入滲的地下徑流,由于接近地面處為大型溶隙,所以首先被大型溶隙調蓄,在大型溶隙與中型溶隙的交界面上,排出比例也為Y2,一部分繼續下滲,在中型溶隙與小型溶隙交界面上排出比例為Y3,小型溶隙與極小型溶隙交界面上排出比例為Y4,最后剩下的地下徑流經極小型溶隙調蓄后排出整個地下溶隙調蓄系統。以此類推,最后從X5面積上向下入滲的地下徑流,由于接近地面處為極小型溶隙,只經過極小型溶隙的調蓄作用,隨即排出整個地下溶隙調蓄系統。

圖1 從極大型溶隙下滲的地下徑流的調蓄-排出過程

以X1面積上下滲的地下徑流為例,推導各類線性水庫串聯的瞬時單位線。設總的地下徑流為I,則X1面積上分得的徑流為X1I。對于只經過極大型線性水庫調蓄的比例為Y1的出流,水量平衡方程和蓄量方程可寫為

(6)

S=K1Q

(7)

式中:I為地下總徑流;Q為K1線性水庫調蓄后的出流;S為流域蓄量。

求得K1線性水庫單獨調蓄瞬時單位線為

(8)

對于經過K1、K2線性水庫串聯調蓄的比例為Y2的出流,根據線性水庫串聯的瞬時單位線公式,K1、K2線性水庫串聯的瞬時單位線為

(9)

對于經過2個以上線性水庫串聯調蓄的出流,n個線性水庫串聯的瞬時單位線公式為

(10)

式中:Ki(i=1,2,…,n)為線性水庫調蓄系數;A為比例系數,計算時根據串聯水庫的具體類型、數量而定。

然此本后來湮沒無聞,不知何故。魯迅即言“至清,皆稱“《稽中散集》”,仍十卷。其稱“《嵇康文集》者,無聞。”[10]70。

從X2、X3、X4、X5面積上向下入滲的地下徑流經過線性水庫串聯的瞬時單位線可以此類推。地下徑流匯流是線性系統,符合疊加原理,可求得線性水庫串聯的總瞬時單位線為

u(0,t)=∑un(0,t)

(11)

由總瞬時單位線可推得總時段單位線為

(12)

即可求得總的地下徑流經整個地下溶隙系統調蓄后的匯流值。

以上是對巖溶地貌區匯流的處理。對于非巖溶地貌區,地表徑流匯流仍采用滯后演算法,并入巖溶地貌區的地表徑流處理。下滲的地表徑流首先在淺層土壤相對不透水層形成壤中流,運動特性與經過一次小型溶隙調蓄的地下水相似,所以并入這部分處理;最后在深層土壤形成地下徑流匯出流域,運動特性與經過一次極小型溶隙調蓄的地下水相似,并入這部分處理。

2 模型應用

2.1 流域概況

將LKHM應用于云貴高原東南部的平湖流域進行水文模擬。平湖流域屬珠江流域西江水系,位于107°03′E~107°37′E、 25°51′N~26°06′N,行政區劃屬貴州省平塘縣,平湖流域高程見圖2。平湖水文站以上流域面積為1 418 km2,依次為中山、低山地形,地勢北高南低,地層主要有石炭系、二疊系、三疊系。平湖流域主要由新橋河水系和京舟河水系匯合形成,巖溶地貌廣泛發育,主要有巖溶洼地、溶隙、落水洞、漏斗等。河流水系明暗相間,具有典型的喀斯特流域二元三維特點。消水洞洞口發育較小,導致在每年汛期洪水短時間難以排泄,造成不同程度的洪澇災害。

圖2 平湖流域高程

收集平湖流域基礎數據資料,需要流域內5個雨量站多年逐日降雨資料、平湖站E601B蒸發皿多年逐日實測水面蒸發資料與多年逐日流量數據。選取2006—2012年(缺2011年)共6年降雨、蒸發、流量數據進行模型的率定與檢驗,其中2006—2009年數據用于模型的參數率定,2010與2012年數據用于模型的檢驗。

2.2 參數率定

采用考慮擁擠度的多目標粒子群優化算法[14](MOPSO-CD)進行LKHM參數率定。選用了兩個目標函數,一個是基于確定性系數的目標函數F1:

(13)

另一個是總徑流量的相對誤差與洪峰流量的相對誤差各取0.5作為權重的目標函數F2:

(14)

總徑流量與洪峰流量分別從整體和局部控制流量過程,為防止程序自動優選參數時只考慮一方最優而放大另一方的相對誤差,本研究將兩者綜合為一個目標函數,并認為重要性相等,權重均取為0.5。F1與F2兩者相互制約[14-15],即一方最優時會犧牲另一方的取值,因此選用這兩個目標函數組成Pareto前沿來評價參數值的率定。

LKHM中一共包括27個參數,在模型中輸入參數的上下限[16],確定參數的尋優空間。上下限可以根據參數的物理意義來確定,也可以根據流域特性經驗確定。27個參數的物理意義以及上下限見表1。參數個數雖多,但相同種類的參數具有一定的聯系,在非劣解集評價時可人為設置約束條件進行非劣解的選取。滯后演算法中滯時τ未設為參數,通過比較流域內多場降雨峰值與洪水峰值滯后情況,將τ值定為1 d。

表1 參數物理意義及上下限

用2006—2009年的數據率定出模型的參數,再用檢驗期2010年與2012年降雨量、蒸發量數據計算出2010年與2012年模擬徑流深。與實測徑流深比較,計算2006—2010年、2012年年徑流深模擬值與實測值的絕對誤差、相對誤差與確定性系數(表2)。2006—2010年、2012年模擬與實測流量過程見圖3~8。

表2 率定期與檢驗期年徑流深模擬值與實測值的相對誤差、確定性系數

圖3 2006年模擬與實測流量過程線

圖5 2008年模擬與實測流量過程線

圖6 2009年模擬與實測流量過程線

圖7 2010年模擬與實測流量過程線

圖8 2012年模擬與實測流量過程線

由表2可知,年徑流深相對誤差絕對值最大為7.90%,平均相對誤差為4.14%,表明年徑流深基本是平衡的,LKHM對年徑流深模擬具有很高的精度。結合表2分析圖3~8可知,圖4中2007年模擬精度最高,確定性系數高達0.81,整個流量過程都基本重合;圖3中2006年、圖7中2010年、圖8中2012年確定性系數相似,總體模擬精度較高,2006、2010年洪峰模擬精度稍差,2012年除洪峰外的流量過程模擬精度稍差;圖6中2009年大洪峰和峰現時間能夠模擬出來,若干小洪峰和洪峰過后的退水過程模擬較差;圖5中的2008年精度最差,確定性系數僅為0.57,次大一級的洪峰誤差很大,流量過程誤差也較大。初步分析原因是該年度汛期內降雨量大,洪水過程較為頻繁,洪峰較多,導致模擬有困難,可取之處是峰現時間基本無偏差。6年確定性系數平均值為0.74,除了2008年以外各年模擬與實測流量過程基本吻合,峰現時間幾乎毫無偏差。總體而言LKHM精度較高,在無資料喀斯特地區具備一定的應用價值。

3 結 論

a. 構建了應用在喀斯特流域的LKHM。LKHM中多處采用概化研究的方法,將復雜的巖溶地貌概化為3類:溶溝、溶隙、管道。溶溝的調蓄作用采用滯后演算法模擬,溶隙的調蓄作用采用線性水庫串聯模擬,管道的調蓄作用采用概化為一個大管道的方法模擬。概化過程具有很強的經驗性,對研究區資料要求較低。

b. LKHM應用在云貴高原東南部的平湖流域,利用2006—2012年(缺2011年)共6年的降雨、蒸發、流量數據進行模擬。模擬結果年徑流深相對誤差絕對值最大為7.90%,平均相對誤差為4.14%,確定性系數平均值為0.74,表明LKHM精度較高,適用于資料缺乏的喀斯特流域。

[ 1 ] 雷曉輝,蔣云鐘,王浩,等. 分布式水文模型Easy DHM[M]. 北京: 中國水利水電出版社, 2010: 265-266.

[ 2 ] 嚴啟坤. 一個巖溶地下河流域模型及其應用[J]. 中國巖溶,1988, 7(2): 25-34. (YAN Qikun. A model of karstic underground drainage system and its application [J]. Carsologica Sinica, 1988, 7(2): 25-34. (in Chinese))

[ 3 ] 袁道先,戴愛國,蔡五田,等. 中國南方裸露型巖溶峰叢山區巖溶水系統及其數學模型的研究:以桂林丫吉村為例[M]. 桂林: 廣西師范大學出版社, 1996: 88-118.[ 4 ] SCHOMBERG J D, HOST G, JOHNSON L B. Evaluating the influence of landform, surficial geology, and land use on streams using hydrologic simulation modeling[J]. Aquatic Sciences, 2005, 67(4): 528-540.

[ 5 ] 任啟偉. 基于改進SWAT模型的西南巖溶流域水量評價方法研究[D]. 北京:中國地質大學, 2006.

[ 6 ] 代俊峰,郭純青,方榮杰. 西南巖溶灌區水文特性及其模擬模型的構建[J]. 水資源與水工程學報,2011, 22(4): 11-15. (DAI Junfeng, GUO Chunqing, FANG Rongjie. Hydrologic characteristics in southwest karst irrigated area and construction of simulation model[J]. Journal of Water Resources and Water Engineering, 2011, 22(4): 11-15. (in Chinese))

[ 7 ] 索立濤,萬軍偉,盧學偉. TOPMODEL模型在巖溶地區的改進與應用[J]. 中國巖溶,2007, 26(1): 67-70. (SUO Litao, WAN Junwei, LU Xuewei. Improvement and application of TOPMODEL in karst region[J]. Carsologica Sinica, 2007, 26(1): 67-70. (in Chinese))[ 8 ] 潘歡迎. 巖溶流域水文模型及應用研究[D]. 武漢:中國地質大學, 2014.

[ 9 ] CAMPBELL C W, SULLIVAN S M. Simulating time-varying cave flow and water levels using the storm water management model[J]. Engineering Geology, 2002, 65(2/3): 133-139.

[10] PETERSON E W, WICKS C M. Assessing the importance of conduit geometry and physical parameters in karst systems using the storm water management model (SWMM)[J]. Journal of Hydrology, 2006, 329(1/2): 294-305.

[11] 章程,蔣勇軍,袁道先,等. 利用SWMM模型模擬巖溶峰叢洼地系統降雨徑流過程:以桂林丫吉試驗場為例[J]. 水文地質工程地質,2007, 34(3): 10-14. (ZHANG Cheng, JIANG Yongjun, YUAN Daoxian, et al. Rainfall-runoff simulation of a typical karst fengcong depression system using SWMM model:a case study of the Yaji experimental site in Guilin[J]. Hydrogeology & Engineering Geology, 2007, 34(3): 10-14. (in Chinese))

[12] 王臘春,史運良. 西南喀斯特山區三水轉化與水資源過程及合理利用[J]. 地理科學,2006, 26(2): 173-178. (WANG Lachun, SHI Yunliang. Formation process and rational use of water resources and transform of rainfall, surface water and underground water in karst mountainous area in southwest China[J]. Scientia Geographica Sinica, 2006, 26(2): 173-178. (in Chinese))

[13] 趙人俊. 流域水文模擬:新安江模型與陜北模型[M]. 北京: 水利電力出版社, 1984.

[14] 宋萬禎,雷曉輝,黃曉敏,等. 考慮擁擠度的多目標粒子群優化算法在馬斯京根參數估計中的應用[J]. 水電能源科學,2013, 31(1): 38-41. (SONG Wanzhen, LEI Xiaohui, HUANG Xiaomin,et al. Application of multi-objective particle swarm optimization in muskingum parameter estimation considering crowding distance[J]. Water Resources and Power, 2013, 31(1): 38-41. (in Chinese))

[15] 王宇暉. 基于分布式水文模型的流域水循環及面源伴生過程模擬研究與應用[D].上海:東華大學, 2012.

[16] 趙人俊,王佩蘭. 新安江模型參數的分析[J]. 水文,1988(6): 2-9. (ZHAO Renjun, WANG Peilan. Analysis of Xinanjiang model parameter[J]. Journal of China Hydrology, 1988(6): 2-9. (in Chinese))

Lumped karst hydrological model and its application

XU Boliu, DONG Zengchuan, HONG Xian

(CollegeofHydrologyandWaterResources,HohaiUniversity,Nanjing210098,China)

Based on the fact that most karst basins are ungauged basins, a lumped karst hydrological model called the LKHM, which requires fewer data than typical models, was established in this study. According to the two-element and three-dimensional properties of karst basins, the LKHM was established based on the three-water sources of the Xin’anjiang model and improved in its calculation of confluence. The karst landform was classified into three types through a generalizability study: grooves, fractures, and pipes. The lag-and-route method was used to simulate the storage effects of grooves on surface runoff. Several linear reservoirs in series were used to simulate the storage effects of fractures with different diameters on subsurface runoff. All pipes were generalized into one large pipe in order to simulate the exchange of surface water and subsurface runoff. The LKHM was applied to the Pinghu Basin in Guizhou Province. The results show that the model has high accuracy, indicating that its parameters and structure are reasonable and the model is applicable.

lumped karst hydrological model; karst landform; grooves; pipe; fracture; linear reservoir

10.3880/j.issn.1004-6933.2017.02.007

國家自然科學基金(51409091)

許波劉(1991—),男,碩士研究生,研究方向為水資源規劃與管理。E-mail:982395181@qq.com

董增川,教授,博士生導師。E-mail:dongzengchuan@163.com

P334.92

A

1004-6933(2017)02-0037-06

2016-08-04 編輯:徐 娟)

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 玩两个丰满老熟女久久网| 3344在线观看无码| 久久精品中文字幕免费| 免费国产高清视频| 午夜国产不卡在线观看视频| 美女啪啪无遮挡| 日韩精品专区免费无码aⅴ | 国产美女主播一级成人毛片| 免费国产无遮挡又黄又爽| 二级毛片免费观看全程| 中文字幕乱码二三区免费| 国产又粗又爽视频| 国产色爱av资源综合区| 色欲不卡无码一区二区| 六月婷婷综合| 无码一区中文字幕| 狠狠做深爱婷婷久久一区| 麻豆国产原创视频在线播放| 亚洲品质国产精品无码| 四虎成人免费毛片| 欧美激情视频二区| 18禁黄无遮挡网站| 亚洲欧美精品一中文字幕| 精品一区二区三区无码视频无码| 夜夜操天天摸| 国产精品一区不卡| 波多野结衣在线se| 又粗又大又爽又紧免费视频| 国产成人免费观看在线视频| 亚洲精品在线91| 国产后式a一视频| 日本午夜视频在线观看| 国产本道久久一区二区三区| 永久成人无码激情视频免费| a级免费视频| 欧美精品不卡| 一级看片免费视频| 国产欧美日本在线观看| 永久成人无码激情视频免费| 国产一区二区三区在线无码| 国产精品密蕾丝视频| 欧美不卡视频在线| 国产精品99r8在线观看| 国产乱子伦无码精品小说| 大陆精大陆国产国语精品1024| JIZZ亚洲国产| 午夜天堂视频| 奇米影视狠狠精品7777| 亚洲手机在线| 欧美在线三级| 亚洲视频免费播放| 国内精品视频在线| 97久久免费视频| 亚州AV秘 一区二区三区| 婷婷色中文网| 国产精品3p视频| 在线亚洲天堂| 亚洲无码高清一区二区| 好久久免费视频高清| 色偷偷男人的天堂亚洲av| 国产福利在线免费| 精品一区二区无码av| 四虎永久在线精品影院| 日本午夜影院| 欧美日韩精品综合在线一区| 2021无码专区人妻系列日韩| 国产区91| 人妻中文字幕无码久久一区| 中文字幕不卡免费高清视频| 亚洲成年网站在线观看| 日本人又色又爽的视频| 人妻一本久道久久综合久久鬼色| 国产精彩视频在线观看| 亚洲狼网站狼狼鲁亚洲下载| 亚洲第一福利视频导航| 2020国产免费久久精品99| 国产第二十一页| 国产屁屁影院| 亚洲色图欧美视频| 伊人久综合| 原味小视频在线www国产| 19国产精品麻豆免费观看|