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

基于替代模型的地下水DNAPLs污染源反演識別

2019-02-13 06:09:52侯澤宇盧文喜
中國環境科學 2019年1期
關鍵詞:優化模型

侯澤宇,盧文喜*,王 宇

?

基于替代模型的地下水DNAPLs污染源反演識別

侯澤宇1,2,盧文喜1,2*,王 宇1,2

(1.吉林大學,地下水與資源環境教育部重點實驗室,吉林 長春 130021;2.吉林大學環境與資源學院,吉林 長春 130021)

應用基于核極限學習機替代模型的模擬-優化理論和方法研究解決了地下水DNAPLs污染源及含水層參數的同步反演識別問題.結果表明:1)核極限學習機替代模型對模擬模型有較高的逼近精度,能夠識別并模仿模擬模型的輸入-輸出關系,絕大部分相對誤差小于5%,平均相對誤差僅有2.98%;2)以替代模型代替模擬模型,大幅度地減小了模擬-優化過程的計算負荷,將反演識別時間由傳統方法的83天減少到3小時,并能夠保持較高的計算精度;3)應用基于模擬退火的粒子群優化算法求解優化模型,能夠以較快的速度搜尋到全局最優,同時避免搜索過程陷于局部極小解.

DNAPLs;污染源反演識別;模擬-優化;多相流模擬;核極限學習機替代模型

石油類有機污染物在水中的溶解度一般很小,進入含水層后通常以非水相流體(Non-Aqueous Phase Liquids, NAPLs)的形式存在.在與水接觸的過程中,NAPLs會不斷向水中溶解釋放.密度大于水的重非水相流體(Dense NAPLs, DNAPLs)具有高密度、低水溶性和高界面張力的特性,實際修復過程中,常常面臨污染物去除率低、修復過程時耗長和修復費用昂貴的困難[1].因此,制定合理高效的修復方案就顯得尤為重要.而合理高效修復方案的制定則要以對于含水層中DNAPLs污染源狀況的識別和掌握為前提.

由于地下水埋藏于地下巖土介質之中,地下水污染通常具有存在的隱蔽性和發現的滯后性特點,致使人們對于地下水污染源的狀況都缺乏了解和掌握.這給地下水污染修復方案設計、風險評估和污染責任認定都帶來了很大的困難.因此,有關地下水污染源反演識別的研究就顯得格外重要[2-3].

地下水污染源反演識別是指根據有限的地下水污染監測數據,反演求解描述地下水污染的數學模擬模型,從而識別確定地下水污染源的個數、空間位置和釋放歷史等相關信息[4-5],在數學上屬于數理方程反問題(或稱數理方程反演問題).相對而言,數理方程的正問題(或稱正演問題)具有較長的研究和應用歷史.正演問題一般都屬適定性問題,是發展較為完善的問題.而反演問題的發展歷史相對較短,興起于20世紀60年代,其理論和應用尚處于發展階段.反演問題大都具有非線性和病態性(不適定性)的特點,原因在于求解反演問題時,已知的信息量遠遠少于待求的信息量,這也正是反演問題的難點之所在[6].

目前地下水污染源反演識別的研究現狀和發展趨勢具有如下特點:1.研究對象多為假想例子[7-10],實際例子十分少見.對于假想例子,可以假定含水層參數是已知的,只有地下水污染源的相關信息是未知待求的.而對于實際例子,不僅地下水污染源的相關信息是未知的,含水層參數也是未知的,待求變量增加.而且含水層參數與污染源之間還是相互影響的,如果含水層參數取值不準確,必然導致污染源的識別結果不準確.因此,今后需要加強對于實際例子的研究,并在反演過程中對污染源和含水層參數都要進行識別.2.模擬-優化方法目前已成為地下水污染源反演識別的主要方法之一[11].然而,若在優化模型的迭代求解過程反復調用模擬模型,會帶來龐大的計算負荷,嚴重制約了模擬-優化方法在反演識別實際應用中的可行性,因此,研究建立模擬模型的替代模型成為近年來研究進展中的前沿問題之一[12-13].3.近年來,遺傳算法、模擬退火算法、蟻群算法等啟發式算法越來越多地被使用到污染源反演識別問題中[14-17].啟發式算法具有全局搜索能力,但對于初始點的依賴仍較為嚴重[18-19].因此,尚有待于尋求不依賴于初值選擇的非線性優化模型的有效解法.本文旨在運用模擬-優化反演的理論和方法,針對這些問題開展進一步的研究,解決地下水重非水相流體污染的污染源反演識別問題.

1 研究方法

地下水污染源反演識別是指根據有限的地下水污染監測數據,對描述地下水污染的數學模擬模型進行反演求解,從而識別確定含水層中地下水污染源的相關信息,包括污染源的個數、空間位置和釋放歷史.其中釋放歷史是指污染物釋放強度隨時間的變化過程,也就是污染物在各時段的釋放強度.

針對前述地下水污染源反演識別研究現狀和發展趨勢中存在的問題,本文基于模擬-優化理論與方法,首先建立污染場地的多相流模擬模型,然后建立模擬模型的核極限學習機替代模型,并通過應用基于模擬退火的粒子群優化算法求解替代模型與優化模型的耦合模型,實現地下水DNAPLs污染源及含水層參數的同步識別.

1.1 模擬-優化方法

模擬-優化方法中的模擬指的是地下水溶質運移模擬模型(簡稱模擬模型),用來描述地下水污染源特征與監測點處污染物濃度的激勵-響應關系.優化指的是優化模型,其目的是解決優選問題.對于污染源反演識別問題,模擬-優化方法是將模擬模型與優化模型耦合,即將模擬模型嵌入優化模型中,使模擬模型成為優化模型的一部分.運用啟發式算法等方法對優化模型進行求解,通過不斷迭代搜索污染源的位置、釋放強度等變量或它們的組合,使得由模擬模型獲得的污染物計算濃度與實際監測濃度最為接近.

1.2 替代模型

優化模型的迭代求解過程需反復調用模擬模型,會帶來龐大的計算負荷,嚴重制約了模擬-優化方法在反演識別實際應用中的可行性[20],因此,研究建立模擬模型的替代模型成為近年來研究進展中的前沿問題之一.替代模型在功能上逼近模擬模型,能夠以很小的計算負荷逼近模擬模型的輸入-輸出響應關系[21-22].本文應用核極限學習機法建立多相流模擬模型的替代模型.

核極限學習機(Kernel Extreme learning machines, KELM)將核函數引入到極限學習機(ELM)的訓練中,采用核映射替代傳統ELM中的隨機映射,能夠產生穩定的輸出結果,其分類和擬合能力均優于非核的ELM方法,泛化能力和學習速度大幅提高.目前尚未見有運用核極限學習機替代模型進行地下水DNAPLs污染源反演識別的研究報道,可通過實例應用分析核極限學習機替代模型的可靠性和優缺點.

ELM的核矩陣為

式中:是訓練樣本輸入在特征空間的映射矩陣.

訓練后的KELM輸出函數表達如下:

1.3 優化模型

優化模型的一般形式為:

粒子群優化算法(Particle swarm optimization, PSO)是一種基于種群協作的隨機尋優算法.模擬鳥群覓食行為的粒子群優化算法與遺傳算法相似,通過迭代尋找最優解.

模擬退火法(Simulated annealing, SA)是一種眾所周知的元啟發式優化算法,模擬高溫固體的退火過程.由于在搜索過程中具有概率突跳能力,模擬退火法可以有效地避免搜索過程陷入局部極小解[19].

將模擬退火法與粒子群優化算法結合,能有效吸取兩種方法的優點,在不依賴于初值選擇的條件下,以較快的收斂速度獲得最優解,避免搜索過程陷于局部極小解.

基于模擬退火的粒子群優化算法的基本運算過程可以概括為[23-24]:

3)確定初始溫度:

4)計算當前溫度下各粒子的適配值:

5)計算慣性權重:

式中:max和min為的最大值和最小值,f為第個粒子當前的目標函數值,f和min為當前所有目標函數值的平均值和最小值.

式中:1和2為0到1之間均勻分布的隨機數.

7)重新計算粒子的適應度,并與之前作對比,更新pbest,gbest.進行退火操作:

式中:為退火常數.=+1,如果達到最大迭代次數,令gbest為最優解;否則返回步驟3.

2 實例研究

2.1 問題描述

以吉林省吉林市某化工污染場地為研究實例,污染很可能是由一個爆炸事故引起的.根據地下水水樣測試分析結果可知,研究區孔隙潛水含水層中的地下水受到了有機污染物——氯苯的污染.收集研究區先前的相關水文地質研究項目報告,根據野外現場調查工作,了解掌握研究區的地質和水文地質條件以及地下水污染源可能的分布狀況.在對已有資料充分分析的基礎上,利用場地內的3口水位水質監測井開展水位、水質同步動態監測工作,每月進行一次水位、水質監測.

2.2 模擬模型的建立

計算目的層為松散巖類孔隙潛水含水層,研究區可概化為非均質、各向同性,并含有若干夾層的三維多相流模型.污染場地附近無天然邊界,可在污染物遷移影響可忽略不計的地段劃定邊界.其中,東西邊界概化為一類邊界,地下水流向由東向西;南北邊界由流面組成,概化為零通量邊界;計算模擬區的上部為潛水面,是水量交換邊界,下部為隔水層,可概化為零通量邊界(圖1).水和氯苯的物理化學參數詳見表1.

圖1 計算區及邊界條件

表1 水和氯苯的物理化學參數

根據污染場地的水文地質概念模型,初步建立如下多相流數學模型:

模型中的源匯項(污染源信息)是未知的,無法進行模型的校正及檢驗,使得模型中的含水層參數取值也無法確定.因此,該模型并不完整,是為進行地下水污染源反演識別而初步建立的.

圖2 計算區剖分示意

在后文的研究中,涉及到對污染源信息及含水層參數賦值并求解多相流模擬模型.進行多相流模擬模型的求解時,運用由美國德克薩斯大學研制的UTCHEM軟件[25].研究區剖分圖見圖2.水質水位監測井的相對位置如圖3所示.

2.3 替代模型的建立

通過敏感性分析(通過給定一個假想的污染源,運行多相流數值模擬模型實現)可知,孔隙度、滲透率、縱向水相彌散度、橫向水相彌散度這4個參數的敏感度較高,而縱向油相彌散度、橫向油相彌散度的敏感度幾乎為零.因此,將孔隙度、滲透率、縱向水相彌散度、橫向水相彌散度作為待識別參數,而將其他參數作為已知常數,根據經驗在其取值范圍內取值.

最終,本文污染源反演識別問題的待識別變量為:污染源的縱向坐標、污染源的橫向坐標、污染物泄漏量、污染物遷移轉化時長、孔隙度、滲透率、縱向水相彌散度、橫向水相彌散度,共8個變量,各參數的初始估計范圍是通過野外現場調研查和專家經驗估計所得.以待識別變量作為優化的可控輸入變量建立替代模型,替代模型的輸出為末時刻三口監測井位置含水層底部的污染物濃度,即本文建立的核極限學習機替代模型有8個輸入變量和3個輸出變量.由于地下水污染物濃度的影響因素眾多,因此,污染物濃度的監測值在一年內有劇烈的變化.本文將三口監測井中的含水層底部污染物濃度監測序列分別取平均值,近似表示各監測井位置末時刻含水層底部污染物濃度實際監測值.

本文設計的核極限學習機替代模型需要100個訓練樣品,20個檢驗樣品.檢驗樣本由不同于訓練樣本的樣品點組成,用來檢驗替代模型對模擬模型的逼近精度.故需要進行2次拉丁超立方抽樣,分別抽取訓練樣本和檢驗樣本.抽樣結束后,將樣本方案多相流數值模擬模型,獲得相應的輸出(當前三口監測井位置含水層底部的污染物濃度).

替代模型訓練樣本數是根據替代模型輸入輸出的維度和經驗確定.通過之前的研究經驗發現,當替代模型的訓練樣本數量較少時,替代模型的精度會隨著訓練樣本數量的增多而顯著提高,而當訓練樣本數較多時,樣本數量對替代模型精度的影響并不明顯.訓練樣本的獲取需要花費比較大的工作量,因此,通過權衡確定訓練樣本數為100.

根據核極限學習機方法的原理在MATLAB平臺上編寫程序,建立核極限學習機模型.運用100組訓練樣本對核極限學習機模型進行訓練,將20組檢驗樣本的輸入數據輸入訓練好的核極限學習機模型,獲得輸出響應(污染物濃度數據),并將其與多相流模擬模型的輸出響應對比.替代模型的輸出變量為3個,由20組檢驗樣本共可獲得60個輸出數據.在與與多相流模擬模型的輸出響應對比后,可得到60個相對誤差數據,相對誤差(絕對值)的分布情況如圖4、圖5所示.結果表明:對于檢驗樣本,替代模型的相對誤差水平很低,絕大部分相對誤差小于5%.

圖4 替代模型相對誤差分布

替代模型的各項精度評估指標如表2所示.核極限學習機模型的各項指標均比較優秀,即核極限學習機模型對模擬模型有較高的逼近精度,能夠識別并取代模擬模型的輸入-輸出關系.綜上,核極限學習機模型可以代替模擬模型與優化模型耦合進行污染源反演識別.

圖5 替代模型相對誤差頻次累積曲線

表2 替代模型精度分析

2.4 優化模型的建立及求解

本文用以識別污染源和含水層參數的優化模型的目標函數可表示為:

將目標函數與污染源相關變量約束以及含水層參數約束等約束條件組合,就構成了污染源反演識別的優化模型.在MATLAB平臺上編寫程序,應用基于模擬退火的粒子群優化算法對建立的非線性規劃優化模型進行求解.通過不斷更新粒子的飛行速度和位置,使粒子的適應度不斷提高,最終得到滿足所有約束條件的最優解.目標函數值收斂曲線如圖6所示.優化求解過程中,在迭代進行到20次時,目標函數的種群最優值即達到了收斂,且目標函數值的下降幅度非常大.說明應用基于模擬退火的粒子群優化算法求解優化模型,能夠以較快的速度搜尋到全局最優,且不依賴于初始值的選擇.

圖6 優化識別過程中目標函數值收斂曲線

求解非線性優化模型獲得的最優解即為污染源與含水層參數的識別結果,如表3所示,識別出的污染源位置見圖7.根據污染源和含水層參數的識別結果重新建立研究區的多相流模擬模型,可以獲得當前含水層中DNAPLs的分布狀況,如圖8所示.

圖中污染物分布中心不在“S”點附近,是由于污染很可能是由一個爆炸事故引起的,污染物屬于瞬時排放,而且根據識別結果,事故發生時間距今超過了十年,油相污染物全部溶于水中,污染中心會隨著水流移動偏離污染源.識別出的污染源位置可以為污染責任認定提供可靠依據,而識別出的污染物分布情況對于修復方案設計和污染風險評價具有重要參考價值.

表3 污染源與含水層參數識別結果

模擬-優化方法的主要計算負荷來自于反復多次調用模擬模型.該硝基苯污染含水層多相流數值模擬模型在CPU為Intel core i5 3.0GHz內存為8GB的計算機上運行一次平均需要大約400s.如果用傳統的模擬-優化方法來解決此問題,模擬模型需要運行18000(遺傳算法初始種群數量為60,遺傳代數為300)次,花費約83d時間.

圖7 污染源位置識別結果示意

圖8 污染物分布情況示意

核極限學習機模型運行一次需要0.8s,如果將核極限學習機模型引入模擬-優化模型中來代替模擬模型,優化求解過程只需要3h.因此,核極限學習機模型的引入在保證精度的同時大大減小了整個計算過程的計算負荷.

3 結論

3.1 核極限學習機替代模型對模擬模型有較高的逼近精度,能夠識別并模仿模擬模型的輸入-輸出關系,絕大部分相對誤差小于5%,平均相對誤差僅有2.98%,可以代替模擬模型與優化模型耦合進行污染源反演識別.

3.2 應用基于模擬退火的粒子群優化算法可以高效地求解優化模型.優化求解過程中,在迭代進行到20次時,目標函數的種群最優值即達到了收斂,而且目標函數值的下降幅度非常大.說明應用基于模擬退火的粒子群優化算法求解優化模型,能夠以較快的速度搜尋到全局最優,且不依賴于初始值的選擇.

3.3 以替代模型代替模擬模型,大幅度地減小了模擬-優化過程的計算負荷,將反演識別時間由傳統方法的83d減少到3h,并能夠保持較高的計算精度.

3.4 應用基于替代模型的模擬-優化方法可以有效地完成地下水DNAPLs污染源與含水層參數的同步識別工作,為修復方案設計、污染風險評估和污染責任定提供可靠依據.

[1] Qin X S, Huang G H, He L. Simulation and optimization technologies for petroleum waste management and remediation process control [J]. Journal of Environmental Management, 2009,90(1): 54-76.

[2] Ayvaz M T. A linked simulation–optimization model for solving the unknown groundwater pollution source identification problems [J]. Journal of Contaminant Hydrology, 2010,117(1-4):46-59.

[3] 趙 瑩.基于模擬-優化方法的地下水污染源反演識別[D]. 長春:吉林大學, 2017.Zhao Y. Inversion identification of groundwater pollution source based on simulation-optimization method [D]. Changchun: Jilin University, 2017.

[4] Sun N Z. Inverse problems in groundwater modeling [M]. Springer Netherlands, 2009.

[5] 張 宇.基于隨機統計方法的地下水污染源反演識別研究 [D]. 長春:吉林大學, 2016.Zhang Y. Inversion recognition of groundwater pollution source based on stochastic statistical method [D]. Changchun: Jilin University, 2016.

[6] Mahinthakumar G K, Sayeed M. Hybrid genetic algorithm – Local search methods for solving groundwater source identification inverse problems [J]. Journal of Water Resources Planning and Management, 2005,131(1):45–57.

[7] Liu X, Cardiff M A, Kitanidis P K. Parameter estimation in nonlinear environmental problems [J]. Stochastic Environmental Research and Risk Assessment, 2010,24(7):1003–1022.

[8] Wang H, Jin X. Characterization of groundwater contaminant source using Bayesian method [J]. Stochastic Environmental Research and Risk Assessment, 2013,27(4):867-876.

[9] Gurarslan G, Karahan H. Solving inverse problems of groundwater- pollution-source identification using a differential evolution algorithm [J]. Hydrogeology Journal, 2015,23(6):1109-1119.

[10] 顧文龍.基于MCMC方法的地下水污染源反演識別[D]. 長春:吉林大學, 2017.Gu W L. Inversion identification of groundwater pollution source based on MCMC method [D]. Changchun: Jilin University, 2017.

[11] Srivastava D, Singh R M. Breakthrough curves characterization and identification of an unknown pollution source in groundwater system using an artificial neural network (ANN) [J]. Environmental Forensics, 2014,15(2):175-189.

[12] Zhao Y, Lu W X, Xiao C N. A Kriging surrogate model coupled in simulation-optimization approach for identifying release history of groundwater sources [J]. Journal of Contaminant Hydrology, 2016,185: 51-60.

[13] 肖傳寧.地下水污染源反演識別研究[D]. 長春:吉林大學, 2017.Xiao C N. Inversion identification of groundwater pollution sources [D]. Changchun: Jilin University, 2017.

[14] 李守巨,劉迎曦.基于蟻群算法的地下水污染源辨識方法[J]. 計算機科學, 2005,32(8):289-291.Li S J, Liu Y X. Groundwater pollution source identification method based on ant colony algorithm [J]. Computer Science, 2005,32(8):289- 291.

[15] 江思珉,張亞力,蔡 奕,等.單純形模擬退火算法反演地下水污染源強度[J]. 同濟大學學報(自然科學版), 2013,41(2):253-257.Jiang S M, Zhang Y L, Xai Y. Inversion of groundwater pollution source intensity by simplex simulated annealing algorithm [J]. Journal of Tongji University (Natural Science Edition), 2013,41(2):253-257.

[16] Jha M, Datta B. Application of dedicated monitoring-network design for unknown pollutant-source identification based on dynamic time warping [J]. Journal of Water Resources Planning and Management, 2015,141(11):04015022.

[17] 肖傳寧,盧文喜,安永凱,等.基于兩種耦合方法的模擬-優化模型在地下水污染源識別中的對比[J]. 中國環境科學, 2015,35(8):2393- 2399.Xiao C N, Lu W X, An Y K, et al. Comparison of simulation- optimization models based on two coupling methods in the identification of groundwater pollution sources [J]. China Environmental Science, 2015,35(8):2393-2399.

[18] 李敏強,寇紀淞,林 丹,等.遺傳算法的基本理論與應用[M]. 北京:科學出版社, 2002:63-65.Li M Q, Kou J S, Lin D. The basic theory and application of genetic algorithm [M]. Beijing: Science Press, 2002:63-65.

[19] 余 麗.地理網絡中求解TSP的并行混合算法研究[D]. 北京:中國地質大學, 2013. Yu L. Research on Parallel Hybrid Algorithm for Solving TSP in Geographical Network [D]. Beijing: China Geosciences University, 2013.

[20] Sreekanth J, Datta B. Multi-objective management of saltwater intrusion in coastal aquifers using genetic programming and modular neural network based surrogate models [J]. Journal of Hydrology, 2010,393(3/4):245–256.

[21] Yao W, Chen X, Huang Y, et al. A surrogate-based optimization method with RBF neural network enhanced by linear interpolation and hybrid infill strategy [J]. Optimization Methods and Software, 2014, 29(2):406-429.

[22] 姜 雪.DNAPLs污染含水層修復方案優選及不確定性分析[D]. 長春:吉林大學, 2016.Jiang X. Optimization and uncertainty analysis of DNAPLs contaminated aquifer repair scheme [D]. Changchun: Jilin University, 2016.

[23] Manoochehri M, Kolahan F. Integration of artificial neural network and simulated annealing algorithm to optimize deep drawing process [J]. International Journal of Advanced Manufacturing Technology, 2014,73(1-4):241-249.

[24] 侯澤宇.基于替代模型的DNAPLs污染含水層修復方案優選過程的不確定性分析[D]. 長春:吉林大學, 2015. Hou Z Y. Uncertainty analysis of optimization process of DNAPLs contaminated aquifer restoration scheme based on surrogate model [D]. Changchun: Jilin University, 2015.

[25] Delshad M, Pope G A, Sepehrnoori K. A compositional simulator for modeling surfactant enhanced aquifer remediation, 1. formulation [J]. Journal of Contaminant Hydrology, 1996,23(4):303–327.

Surrogate-based source identification of DNAPLs-contaminated groundwater.

HOU Ze-yu1,2, LU Wen-xi1,2*, WANG Yu1,2

(1.Key Laboratory of Groundwater Resources and Environment Ministry of Education, Jilin University, Changchun 130021, China;2.College of Environment and Resources, Jilin University, Changchun 130021, China)., 2019,39(1):188~195

Groundwater contamination source identification (GCSI) is critical for taking effective actions in designing remediation strategies, estimating risks, and confirming responsibility. Surrogate-based simulation-optimization technique was applied to source identification and parameter estimation of DNAPLs-contaminated aquifer in this article. The results showed that: 1) kernel extreme learning machines (KELM) surrogate model approximated the simulation model accurately. It could simulate the input/output relationship of the simulation model with most of the relative errors less than 5%, and the mean relative error was only 2.98%; 2) Replacing the simulation model with a KELM model considerably reduced the computational burden of the simulation-optimization process and maintained high computation accuracy, the identification time was reduced to 3hours from 83days; 3)Simulated annealing-based particle swarm optimization algorithm is efficient in searching the global optimal solution of the nonlinear programming optimization model, and avoiding the optimization process trapping into local optimum.

Dense non-aqueous phase liquids (DNAPLs);contamination source identification;simulation-optimization;multi-phase flow simulation;KELM surrogate model

X523

A

1000-6923(2019)01-0188-08

侯澤宇(1989-),男,河北唐山人,博士后,主要從事地下水污染控制與修復研究.發表論文20余篇.

2018-06-04

國家自然科學基金項目(41672232);吉林省科技發展計劃項目(20170101066JC)

* 責任作者, 教授, luwx999@163.com

猜你喜歡
優化模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 福利在线免费视频| 日韩高清在线观看不卡一区二区| 国产亚洲精| 午夜少妇精品视频小电影| 久久精品亚洲中文字幕乱码| 亚洲AV成人一区二区三区AV| 国产特级毛片| 国产青榴视频在线观看网站| 熟女日韩精品2区| 亚洲成综合人影院在院播放| 狠狠做深爱婷婷久久一区| 国产成人精品午夜视频'| 国产福利免费视频| 乱人伦99久久| 精品无码人妻一区二区| 欧美午夜网| 中文字幕人妻av一区二区| 国产成年女人特黄特色大片免费| 亚洲三级视频在线观看| 日韩欧美在线观看| 人人妻人人澡人人爽欧美一区| 国产成人一区二区| 一级毛片免费的| 精品久久久久无码| 热伊人99re久久精品最新地| 亚洲开心婷婷中文字幕| 精品一区二区三区自慰喷水| 手机在线国产精品| 国产成人亚洲无吗淙合青草| 久久这里只有精品66| 欧美精品1区| 幺女国产一级毛片| 日本精品一在线观看视频| 中美日韩在线网免费毛片视频 | 熟女成人国产精品视频| 2048国产精品原创综合在线| 久久夜色精品国产嚕嚕亚洲av| 一本久道久久综合多人| 欧美日韩一区二区在线播放 | 精品国产99久久| 久久黄色视频影| 欧美成人精品在线| 国产成人成人一区二区| 亚洲精品欧美日本中文字幕| 热re99久久精品国99热| 免费中文字幕在在线不卡 | 久久99热这里只有精品免费看| 国内嫩模私拍精品视频| 欧美成人亚洲综合精品欧美激情| 亚洲二三区| 久久黄色毛片| 一本大道香蕉久中文在线播放| 国产色伊人| 国产成人高清亚洲一区久久| 少妇极品熟妇人妻专区视频| 超薄丝袜足j国产在线视频| 日本在线欧美在线| 小13箩利洗澡无码视频免费网站| 国产专区综合另类日韩一区| 亚洲视频色图| www欧美在线观看| a色毛片免费视频| 欧美一区二区丝袜高跟鞋| 欧美国产日本高清不卡| 日韩区欧美区| 久久国产亚洲偷自| 欧美一级在线| 日日拍夜夜嗷嗷叫国产| 日韩成人免费网站| 亚洲黄色激情网站| 日本久久网站| 亚洲精品无码不卡在线播放| 福利视频久久| 毛片免费在线视频| 无遮挡国产高潮视频免费观看| 曰韩人妻一区二区三区| 国内精品自在自线视频香蕉| 中文无码精品A∨在线观看不卡 | 国产亚洲视频免费播放| 欧洲一区二区三区无码| 99在线观看国产| 国产在线98福利播放视频免费|