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

撫仙湖流域負荷削減的水質風險分析

2013-09-07 08:18:08張曉玲顏小品陽平堅TetraTechInc006EatonPlaceSte40FairfaxVA00USA北京大學環境科學與工程學院水沙科學教育部重點實驗室北京0087云南省環境科學研究院云南省高原湖泊流域污染過程與管理重點實驗室云南昆明65004
中國環境科學 2013年9期
關鍵詞:水質污染模型

鄒 銳 ,張曉玲 ,劉 永 ,趙 磊 ,朱 翔 *,顏小品 ,陽平堅 (.Tetra Tech, Inc. 006 Eaton Place, Ste 40, Fairfax, VA 00, USA.北京大學環境科學與工程學院,水沙科學教育部重點實驗室,北京 0087.云南省環境科學研究院,云南省高原湖泊流域污染過程與管理重點實驗室,云南 昆明 65004)

流域負荷削減,如最大日負荷總量(TMDL)計劃,被認為是改善水質最為直接和有效的手段.TMDL的核心思想是指在滿足水質標準的條件下,水體能接受的某種污染物的最大日負荷總量;經過不斷的改進和發展,TMDL計劃已逐步形成一套完整系統的總量控制策略和技術方法體系,成為水環境管理與決策的主要發展趨勢[1-3].在美國紐約州的奧農達加湖、北卡羅萊納州的紐斯河、佛羅里達州的奧基喬比湖、加利福尼亞南部的新港灣流域及我國的東湖、三峽庫區巫溪段等眾多案例的研究表明,推行TMDL計劃可有效改善水體的水質[4-12].

在制定TMDL的過程中,揭示并定量表征流域污染負荷與水質之間的響應關系是最為重要的步驟,其中選擇合適的模型來模擬負荷與水質響應是關鍵要素[13].復雜的三維水動力-水質數值模型近年來廣泛應用于TMDL研究的重要基礎[14-16].數值機理模型可實現對湖泊三維水流及相應的營養鹽遷移轉化的模擬,準確再現污染物在湖體的時空分布特征.同時,在實際的流域決策中,由于在數據、模型結果等方面存在不確定性,使得開展基于風險分析的TMDL研究成為必然,并藉此為決策者提供更為準確和多樣的支撐[16-19].所謂基于風險的污染負荷控制,是針對傳統的TMDL對環境標準的解譯方式提出的.在傳統的TMDL框架下,水質目標多是靜態解譯的,比如營養鹽,美國最常用的是生長季節平均值或者月均值等,一般不對允許的超標頻度做顯性的表達.但在實際操作中,流域管理者常常面臨費用與污染控制效益之間矛盾的問題.因此,決策者需要根據技術經濟可行性,考慮允許適當的超標風險,就需要研究相應于此風險的容量水平.此時,最可靠的方法就是用三維數值模型對不同達標頻度下的情景進行模擬分析,但由于數值模型中非線性過程的存在,直接應用三維數值模型存在計算量過大的瓶頸,因此無法在實際的TMDL中得到廣泛應用.人工神經網絡模型[20-21](NN)能夠模仿人腦的功能,在預先不知道具體函數形式的情況下,對模糊信息或復雜的非線性關系進行識別與映射,在環境[22-24]領域應用廣泛.本文以云南高原湖泊撫仙湖為例,基于前期研究建立的三維水動力-水質模型,在模型校正和驗證的基礎上,采用ANN耦合數值水質模型方法,以30d移動平均值為解譯方法,對撫仙湖流域污染負荷與湖泊達標風險之間的關系進行定量化,繼而核算不同風險下的 TMDL,為在撫仙湖實施智能流域管理和決策提供支持[25].

1 撫仙湖三維水動力-水質模型的構建

撫仙湖三維水動力-水質模型的構建基于美國 EPA開發的環境流體動力學模型(EFDC)平臺[14-18],在前期的研究中[3],已經建立了整個湖體的三維動態模型,并基于 4個監測點位進行模擬校正,主要的模擬指標為COD、TN和TP;主要的模擬步驟包括:網格生成、初始條件配置、邊界條件設定、模型校準及應用.為精確描述撫仙湖的湖岸線,采用了曲線網格法;將撫仙湖在水平方向上劃分為 323個正交曲線網格,每個網格在垂直方向上按西格瑪坐標切割成為50層以代表撫仙湖的深水特征.設置2009年1月1日觀察到的水面高度1723.28m為初始高程,初始溫度以1月初的觀測值 13.5℃為基礎,3個速度向量按水動力學常規初始化為0.0m/s,挑選2009年1月6日獲得的水質數據TP 0.005mg/L、TN 0.171mg/L、COD 0.98mg/L為初始條件,湖流和營養物質的水平邊界條件是以2009年1~12月主要入湖河流的流域模型結果為基礎;驅動流體模型的大氣邊界數據來自澄江縣氣象站獲得的每小時的氣象數據.在模擬中,首先對水動力過程進行校準,計算步長為30s,模擬的時間涵蓋2009年全年,模擬變量為流場、水位和溫度,校準參數是湖水水位和水溫.在水動力校準的基礎上,對水質過程進行模擬.撫仙湖水質模型的模擬時段與水動力模型相同,關鍵校驗參數為各個水質狀態變量的一階衰減速率及沉降速率,模型進行了大約20次迭代校準.模型經過校正和驗證,模擬的流場、水位、溫度及水質與觀測到的時空分布狀況匹配良好,表明模型可用于進一步的決策分析.在此基礎上,應用該模型對2種TMDL情景進行了分析[3],發現在2種極端的風險水平上,核算的TMDL值差別很大,表明要實現對流域管理的有效支持,有必要探索風險與流域污染負荷之間的關系,而非僅依據2種極端情景來指導流域管理決策.

2 EFDC-NN耦合模型

在進行 TMDL核算之前,首先必須確定對環境質量標準的解譯方式.對于撫仙湖這類高原湖泊而言,雖然直觀上的水質目標是以營養鹽濃度描述的,但實際上,決策中真正的水質風險是富營養化和藻類爆發問題.考慮到富營養化及藻類的爆發不僅需要較高的營養鹽濃度,同時需要足夠的維持這種較高濃度的時間.基于這種考慮,在前期的研究中[3],以瞬時極大值或全年平均值來解譯水質目標的情景,但這種方法只能提供上、下限范圍,而不具實際決策意義.因此,本文以月為時間尺度來探究營養鹽達標的問題.不同于一般的按月進行評價的方法,采用30d移動平均的方法來表征水質特征,即:在開展TMDL的風險研究時,水質達標的評估方式是以全湖表層30d移動平均濃度滿足GB3838-2002的I類標準[26]為依據.除了時間尺度,在實際管理中另一個重要的決策維度是風險,即湖體TP、TN、COD濃度達標的頻度與流域污染負荷削減(或增加)間的不確定性響應關系.在明確了不同的污染負荷水平對應的湖泊水質風險后,決策者才能根據實際情況做出判斷,制定合理的適應性流域管理方案[25].由于撫仙湖三維水動力-水質模型每次情景運行時間需要 2d以上,而要獲得完整的流域污染負荷和風險響應關系則有可能需要上百次的模擬,顯然這在計算上是不可行的.為了高效實現對撫仙湖污染負荷與水質風險響應關系的定量化,本研究采用EFDC-NN耦合模型將水質風險與流域污染負荷關聯起來.

2.1 耦合模型構建

首先用已開發和校正的 EFDC模型進行了一系列的情景分析,并在模型開發和校驗年的水文條件下,按照月均值的解譯方式要求對 TP、TN、COD結果進行了30d移動平均,以便對神經網絡(NN)進行訓練和驗證;而后以流域污染負荷削減(或增加)與達標頻度為輸入變量,以相應的TP、TN和COD濃度為輸出變量構建神經網絡模型,選用BP網絡架構.一個隱層含有2m+1個神經元的3層神經網絡模型,可以實現一個m維的實向量與一個n維實向量之間的函數映射[20].在撫仙湖BP模型中,由于m=2,因此5個或5個以下隱層神經元就應足以完成對染負荷削減(或增加)和達標頻度與TP、TN和COD之間的函數映射.但由于神經網絡這類人工智能模型,其收斂性和泛化度確認都沒有明確的指標,所以為了表達模型分析可能存在的不確定性,在本研究中構造了6個BP模型,其隱層神經元數目分別為2、3、4、5、6、7.

2.2模型訓練與校正

模型選取EFDC情景分析得到并經過月均處理的700組數據作為訓練數據,另外400組數據作為驗證數據.網絡的訓練采用多次重啟的方法,以避免陷入局部最優網絡,而訓練的收斂條件定為當模型泛化度開始下降為止,這樣做有助于避免過度訓練問題.由于研究對象存在大量的數據,可以有效地解決神經網絡的泛化度問題.在神經網絡對每種污染負荷的預測值和EFDC模擬值間進行線性擬合,線性系數越接近于1,表示二者越近似相等;R2值越接近于 1,表示線性擬合的可信度越高;因此可以用線性系數和R2值表示預測精度.表1展示了不同網絡結構的模擬結果,6個神經網絡對每種污染負荷的響應值近似等于EFDC的模擬值,說明6個神經網絡模型雖然隱含神經元數目不同,但均能較精確地描述真實系統中的響應關系.

表1 六個網絡結構對污染負荷的模擬精度Table 1 Accuracy of the six network structures

2.3 撫仙湖人工神經網絡的響應分析

撫仙湖神經網絡的響應是通過分析 TP、TN和 COD濃度在一系列不同的風險水平(達標頻度)、不同污染負荷削減率/增加比例下的響應來實現的.神經網絡模型本質上是一類統計模型,因此在應用該模型進行響應分析時不能無限度地改變輸入以獲取輸出.鑒于此,在模擬和情景分析中將污染負荷削減率/增加比例的上限設定為100%(污染負荷最多可增加 100%),下限定為20%(污染負荷最多可削減20%);將達標頻度的降低比例上限設定為 25%,即:達標頻度不能低于75%,以確保撫仙湖的水質風險不會太高.然后分別用 6個神經網絡模型模擬污染負荷削減率/增加比例在不同的達標頻度下的 TP、TN和 COD響應.污染負荷削減率/增加比例和達標頻度均以1%為離散步長,對TN單因子、TP單因子和COD單因子進行一系列的模擬.模擬過程中,先設定達標頻度值,分析污染負荷削減率(或增加比例)每改變1%對TP、TN和COD的響應,再依次改變達標頻度數值.圖 1展示了幾個典型的達標頻度下削減率/增加比例與TP、TN和COD的響應曲線.由圖1可見,TP、TN和COD的濃度都隨著污染負荷削減率的減少和增加比例的增加而升高,隨達標頻度的降低而降低.

圖1 不同達標頻度下削減率/增加比例與TP、TN和COD響應曲線Fig.1 TP, TN and COD responses to watershed loading variations under different risks

3 撫仙湖流域TMDL的風險分析

通過響應分析可知,6個網絡結構的模擬結果基本一致.當達標頻度一定時,隨著污染負荷削減率的減少,TP、TN和COD濃度逐漸升高,經過削減率為0這一特殊點后,隨著污染負荷增加比例的升高,污染物濃度繼續升高.綜合分析基于不同風險的情況,在污染負荷削減率/增加比例一定時,隨著達標頻度的增加,TP、TN和COD的響應值升高,因為達標頻度的增加意味著基于風險的TMDL達標頻率升高、超標頻率降低,會產生濃度更低的污染負荷,從而允許有更高的環境容量.

具體而言,當達標頻度為100%時,為了達到I類水質標準,TP、COD負荷還有繼續上升的空間,TP允許增加比例為14%~18%,COD允許增加比例為 9%~11%.這里的增加水平是區間值,反映了由于采用不同的神經網絡架構產生的不確定性.但是由于 6個神經網絡的模擬結果一致性很高,均高于 99%.在 100%達標頻度水平下,TN 超過 I類水標準,因此需要在現狀負荷的基礎上削減 13%~14%.95%的達標頻度情況與 100%達標頻度相似,但是隨著達標頻度的降低,對污染負荷的要求嚴格程度降低,相應地污染負荷削減率降低、增加比例升高,因而TP、COD允許增加比例分別為 36%~38%、20%~21%,TN 削減比例為4%~6%.對于 90%的達標頻度而言,由于要求更加寬松,TN不再需要進行削減,相反也有增加的余地,TN、TP、COD允許增加比例分別為2%~4%、56%~60%、30%~32%.對于85%的達標頻度而言,TN、TP、COD允許增加比例進一步增加,分別為 9%~12%、76%~82%、39%~41%.對于80%的達標頻度,TP負荷的增加比例達到預設上限(100%)時,仍可滿足 I類水質標準,TN、COD允許增加比例分別為16%~19%、48%~50%.當75%的達標頻度時,TP增加比例仍可達到預設上限,TN、COD允許增加比例分別為22%~24%、56%~58%.

表2和圖2展示了不同風險水平下撫仙湖TN、TP和COD的核算結果(削減率/增加比例取區間中值).由于在75%和80%達標頻度下TP增加比例達到了預設上限,因此本研究中沒有再核算相應的環境容量值.在管理上的含義就是,如果決策者對 TP的超標風險容忍度比較高,那么將不會顯性考慮TP的負荷削減問題.

圖2 TP、TN和COD在不同達標頻度下的TMDLFig.2 TMDL of Lake Fuxian at different risk levels

表2 在不同風險水平下撫仙湖的水環境容量Table 2 Environment capacity of Lake Fuxian at different risk levels

4 討論

從TN、TP、COD的響應可知,撫仙湖水質應以TN為治理重點,其次是COD.雖然TP負荷可增加的空間較大,但作為富營養化的關鍵因子之一也不容忽視.基于風險的TMDL研究為決策管理部門提供了強有力的科學依據,具有非常強的現實意義.但達標頻度具體如何選取,應綜合考慮撫仙湖現實情況及技術經濟可行性,由決策管理部門根據需求來決定.需要指出的是,撫仙湖是深水貧營養湖,內源或藻類動力對營養鹽與流域負荷的響應關系沒有淺水湖泊(如:滇池和異龍湖)明顯;在淺水湖泊中,流域負荷和湖泊水質之間的關系就顯示很強的非線性,因此在淺水湖泊應用本方法時需進行對應的改進.

本研究采用的 30d移動平均值的解譯方式不像以瞬時濃度為解譯方案(即要求在湖體內表層任何地方、任何時段污染物濃度都滿足I類水標準)會導致管理目標過于嚴格,可能面臨經濟和技術上不可行局面,也不像年平均的解譯方案管理目標過于寬松、容易引起水質惡化,是嚴格程度、可實施性及風險性均介于二者之間的解譯方案.

本文的解譯方式是按照一般湖泊標準以月為時間尺度來進行的,但是必須考慮到,30d移動取平均值雖然在大多數時候都是比較可靠的,但是還是存在下面的風險,即以 30d為時間尺度往往會忽略掉這期間的濃度峰值情況;而撫仙湖的敏感水層(表層水)短期營養鹽濃度過高就可能導致藻類的大量爆發,給湖泊生態系統帶來危害.因此,究竟選取何種時間尺度作為解譯方式才能使管理決策更具有保障,還需要根據撫仙湖的實際情況在今后進行深入研究,并顯性地將營養鹽與藻類的動力學關系表達出來.

水環境容量和TMDL并不是新概念,但本文的研究對二者是個有益的擴充(表2).在傳統的環境容量或TMDL計算中,一般是以一種典型情景作為依據來核算環境容量.但如本文所述,所謂水環境容量或者TMDL本身并不是一個孤立的值,而是和管理決策風險緊密相關的變量.因此,在核算TMDL時,如果能夠同時考慮到風險這一決策維度,獲得的結果對實際管理將更為有益.這種基于模型并包括不確定性特征的流域決策模式,充分利用信息和定量決策手段,將為流域決策提供更為直接、多樣和科學的決策集,減少決策失誤,因此也是目前廣受關注的智能流域管理的重要途徑.

5 結論

5.1 構建的三維水動力-水質模型可以有效地模擬撫仙湖流場、水位、溫度與水質的實際情況;6個神經網絡模型對TN、TP和COD的擬合度均大于 0.97,能較精確地反映真實系統的響應關系.

5.2 通過EFDC-NN耦合模型進行風險分析可知,當達標頻度一定時,TN、TP和COD的濃度隨污染負荷削減率的減少而逐漸升高.在污染負荷削減率/增加比例一定時,TN、TP和COD的環境容量隨達標頻度的增加而增加.

5.3 當達標頻度為 100%時,為了達到 I類水質,TP允許增加 14%~18%,COD允許增加9%~11%,TN 需要削減 13%~14%;達標頻度為90%時,TN、TP、COD允許增加比例分別為2%~4%、56%~60%、30%~32%;達標頻度為80%時,TP增加比例達到預設上限,TN、COD允許增加比例分別為16%~19%、48%~50%;達標頻度為75%時,TP增加比例仍可達到預設上限,TN、COD允許增加比例分別為22%~24%、56%~58%.流域管理者可依據不同的風險與管理費用偏好實施流域污染負荷削減.

[1]Elshorbagy A, Teegavarapu R, Ormsbee L. Total maximum daily load (TMDL) approach to surface water quality management:concepts, issues, and applications [J]. Canadian Journal of Civil Engineering, 2005,32(2):442-448.

[2]Leclair V. Courts push states, EPA to create TMDL water programs [J]. Environmental Science and Technology, 1997,31(4):178-179.

[3]Zhao L, Zhang X L, Liu Y, et al. Three-dimensional hydrodynamic and water quality model for TMDL development of Lake Fuxian, China [J]. Journal of Environmental Science,2012,24(8):1355-1363.

[4]王彩艷,彭 虹,張萬順,等.TMDL技術在東湖水污染控制中的應用 [J]. 武漢大學學報(工學版), 2009,(05):665-668.

[5]Effler S W, O'Donnell S M, Matthews D A, et al. Limnologicaland loading information and a phosphorus total maximum daily load (TMDL) analysis for Onondaga Lake [J]. Lake and Reservoir Management, 2002,18(2):87-108.

[6]Stow C A, Borsuk M E. Assessing TMDL effectiveness using flow-adjusted concentrations: A case study of the Neuse River,North Carolina [J]. Environmental Science and Technology, 2003,37(10):2043-2050.

[7]Havens K E, Walker W W. Development of a total phosphorus concentration goal in the TMDL process for Lake Okeechobee,Florida (USA) [J]. Lake and Reservoir Management, 2002,18(3):227-238.

[8]Zheng Y, Keller A A. Stochastic watershed water quality simulation for TMDL development-A case study in the Newport Bay Watershed [J]. Journal of the American Water Resources Association, 2008,44(6):1397-1410.

[9]Santhi C, Williams J R, Dugas W A, et al. Water quality modeling of Bosque River Watershed to support TMDL analysis[J]. Total Maximum Daily Load (TMDL): Environmental Regulations, Proceedings. 2002:33-43.

[10]Bittencourt S, Gobbi E F. Maximum allowable phosphorus load in the Piraquara II Reservoir, a TMDL process application [J].Revista Brasileira de Ciencia Do Solo, 2006,30(3):595-603.

[11]Boyacioglu H, Alpaslan M N. Total maximum daily load (TMDL)based sustainable basin growth and management strategy [J].Environmental Monitoring and Assessment, 2008,146(1-3):411-421.

[12]丁京濤.大寧河巫溪段水體總磷 TMDL估算及分配研究 [D].北京:北京化工大學, 2009.

[13]Munoz-Carpena R, Vellidis G, Shirmohammadi A, et al.Evaluation of modeling tools for TMDL development and implementation [J]. Transactions of the ASABE, 2006,49(4):961-965.

[14]Liu Z J, Kingery W L, Huddleston D H, et al. Modeling nutrient dynamics under critical flow conditions in three tributaries of St.Louis Bay [J]. Journal of Environmental Science and Health Part A-Toxic/Hazardous Substances and Environmental Engineering,2008,43(6):633-645.

[15]Rodriguez H N, Cope B, Peene S J. Hydrodynamic and water quality modeling of ward cove, Alaska [J]. Estuarine and Coastal Modeling, Proceedings, 2004:628-645.

[16]Khangaonkar T, Yang Z Q. Coastal circulation and effluent transport modeling at Cherry Point, Washington [J]. Estuarine and Coastal Modeling, Proceedings, 2004:475-491.

[17]Wool T A, Davie S R, Rodriguez H N. Development of three-dimensional hydrodynamic and water quality models to support total maximum daily load decision process for the Neuse River Estuary, North Carolina [J]. Journal of Water Resources Planning and Management-ASCE, 2003,129(4):295-306.

[18]He G J, Fang H W, Bai S, et al. Application of a three-dimensional eutrophication model for the Beijing Guanting Reservoir, China [J]. Ecological Modeling, 2011,222(8):1491-1501.

[19]Zou R, Carter S, Shoemaker L, et al. An integrated hydrodynamic and water quality modeling system to support nutrient TMDL development for Wissahickon Creek [J]. Journal of Environmental Engineering, 2006,132(4):555-566.

[20]鄒 銳,張禎禎,劉 永,等.神經網絡模型用于數值水質模型逼近的適用性及非敏感參數的欺騙效應 [J]. 環境科學學報,2010,(10):1964-1970.

[21]王 愷,趙 宏,劉愛霞,等.基于風險神經網絡的大氣能見度預測 [J]. 中國環境科學, 2009,29(10):1029-1033.

[22]胡 康,萬金泉,馬邕文,等.基于模糊神經網絡的 A~2/O 工藝出水氨氮在線預測模型 [J]. 中國環境科學, 2012,32(2):260-267.

[23]王 祎,李靜文,邵 雪,等.基于計算智能的流域污染排放優化模式研究 [J]. 中國環境科學, 2012,32(1):173-180.

[24]劉 罡,李 昕,胡 非.大氣污染物濃度的神經網絡預報 [J].中國環境科學, 2000,20(5):429-431.

[25]劉 永,鄒 銳,郭懷成.智能流域管理研究 [M]. 北京:科學出版社, 2012.

[26]GB3838-2002 地表水環境質量標準 [S].

猜你喜歡
水質污染模型
一半模型
水質抽檢豈容造假
環境(2023年5期)2023-06-30 01:20:01
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
一月冬棚養蝦常見水質渾濁,要如何解決?這9大原因及處理方法你要知曉
當代水產(2019年1期)2019-05-16 02:42:04
堅決打好污染防治攻堅戰
當代陜西(2019年7期)2019-04-25 00:22:18
堅決打好污染防治攻堅戰
3D打印中的模型分割與打包
對抗塵污染,遠離“霾”伏
都市麗人(2015年5期)2015-03-20 13:33:49
水質總磷測定存在的問題初探
河南科技(2014年23期)2014-02-27 14:19:07
主站蜘蛛池模板: 国产欧美日韩在线一区| 国产成熟女人性满足视频| 91亚洲视频下载| 亚洲第一极品精品无码| 国产亚洲成AⅤ人片在线观看| 国产资源免费观看| 国产一区二区丝袜高跟鞋| 日韩免费中文字幕| 亚洲欧美成人在线视频| 在线观看视频一区二区| 麻豆精选在线| 亚洲日韩日本中文在线| 欧美不卡二区| 国产亚洲精久久久久久无码AV| 四虎亚洲精品| 色婷婷啪啪| 毛片一级在线| 国产区91| 波多野结衣AV无码久久一区| 97视频在线观看免费视频| 99视频精品全国免费品| 欧美亚洲国产一区| 操国产美女| 国产你懂得| 丰满人妻中出白浆| 91成人免费观看在线观看| 国产一区二区三区夜色| 亚洲无码精彩视频在线观看| 五月激激激综合网色播免费| 成人综合在线观看| 国产h视频免费观看| 亚洲av中文无码乱人伦在线r| 一区二区三区成人| 日本在线视频免费| 日本道中文字幕久久一区| 免费一级毛片完整版在线看| 国产午夜看片| 亚洲综合激情另类专区| 亚洲 成人国产| 中文字幕人成乱码熟女免费| 毛片免费观看视频| 婷婷午夜影院| 亚洲另类第一页| 免费99精品国产自在现线| 曰AV在线无码| 国产自在线播放| 国产正在播放| 久久精品亚洲专区| 丰满人妻久久中文字幕| 看你懂的巨臀中文字幕一区二区| 91久久天天躁狠狠躁夜夜| 在线播放精品一区二区啪视频| 天天色天天综合| 蝴蝶伊人久久中文娱乐网| 亚洲人精品亚洲人成在线| 女人18毛片水真多国产| 国产精品99一区不卡| 亚洲精品日产AⅤ| 一级全黄毛片| 亚洲国产理论片在线播放| 精品国产免费观看一区| 91年精品国产福利线观看久久 | 久久这里只有精品2| 亚洲中文在线看视频一区| 久久黄色一级视频| 无码AV动漫| 91啪在线| 日韩毛片在线视频| 天堂av综合网| 好吊妞欧美视频免费| 亚国产欧美在线人成| 在线免费看黄的网站| 午夜无码一区二区三区在线app| 香蕉久久国产精品免| 青草娱乐极品免费视频| 久久青草精品一区二区三区| 呦女亚洲一区精品| 狠狠色狠狠色综合久久第一次| 9966国产精品视频| 人妻丰满熟妇啪啪| 亚洲成aⅴ人在线观看| 日韩高清欧美|