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

基于混合蛙跳算法的水文模型參數估計方法

2016-05-14 07:00:16火久元劉立群趙紅星
關鍵詞:優化模型

火久元,劉立群,趙紅星

(1.蘭州交通大學 電子與信息學院,蘭州 730070;

2.中國科學院 寒區旱區環境與工程研究所,蘭州 730000;

3.甘肅農業大學 信息科學技術學院,蘭州 730070)

?

基于混合蛙跳算法的水文模型參數估計方法

火久元1,2,劉立群3,趙紅星1

(1.蘭州交通大學 電子與信息學院,蘭州730070;

2.中國科學院 寒區旱區環境與工程研究所,蘭州730000;

3.甘肅農業大學 信息科學技術學院,蘭州730070)

摘要:傳統啟發式算法在水文模型參數估計中通常存在著易早熟和收斂速度慢等缺陷,為提高水文模型的參數優化精度和算法性能,引入混合蛙跳算法(SFLA),提出一種基于SFLA的水文模型參數估計方法,將該方法應用到新安江模型的參數估計中,并與基于遺傳算法(GA)的參數估計方法進行實驗對比分析。實驗結果表明:基于SFLA的參數優化方法在平均優化精度上相比遺傳算法提高了2.5%;在固定優化精度時,優化成功率相比遺傳算法提高了53.33%。證明了混合蛙跳算法應用于水文模型參數估計時,在收斂精度和收斂速度方面均有明顯優勢。

關鍵詞:混合蛙跳算法;新安江模型;參數估計;Nash-Sutcliffe效率系數

流域水文建模是洪水模擬、預報的重要工具之一,也是認識和管理水資源的重要技術手段[1]。水文過程非常復雜,需考慮降雨、土壤水分和蒸散量的時空變化,并涉及土地利用、植被覆蓋、土地坡度、土壤排水性能等相關因素。水文模型一般有大量的參數,包括模型輸入、中期狀態變量和模型輸出[2]。

由于流域水文模型的復雜性、模型結構的不確定性、模型參數較多且多數參數不可測以及觀測數據的誤差等原因,模擬結果往往有著大量的不確定性,導致水文模擬與預報相比存在較大的誤差和不確定性。因此,水文工作者首先要借助模型參數估計方法來確定模型參數值,然后利用水文模型模擬流域水文過程[3]。新安江模型由趙人俊等在1973年提出,并成功應用在中國的濕潤和半濕潤地區。新安江模型有多個參數,其中一些可以直接從流域的物理特性知識中獲取,但大部分參數不能直接通過觀察現場進行測量[4-5]。

目前,眾多學者傾向于研究各種參數優化方法來自動完成水文模型的參數優選過程,并取得了大量的研究成果,提出和引入了很多優化方法。其中被應用到水文模型的參數優化過程中的主要有粒子群算法(PSO)、遺傳算法(GA)、SCE-UA和模擬退火算法(SA)。但是,這些優化算法存在著一些缺陷,例如粒子群算法中的粒子之間協作不充分[6],遺傳算法中的遺傳算子很難確定[7],SCE-UA算法易陷入局部最優,且對模型結構和優化準則要求嚴格[8],模擬退火算法的收斂精度較低,且收斂速度較慢[9]。以上這些優化算法存在的不足之處嚴重影響了水文模型的參數優化結果和優化性能。

2003年,Eusuff 等提出一種基于群體的亞啟發式、協同搜索群智能算法(shuffled frog leaping algorithm,SFLA)[10]。該算法主要模擬青蛙在覓食過程中互相進行信息共享和交流的特征來完成基于群體的優化。混合蛙跳算法(SFLA)結合了基于遺傳基因的MA算法和基于群體覓食行為的粒子群算法二者的優點,因此具有參數少、概念簡單、全局尋優能力強、計算速度快、易于實現等特點[11-12],獲得了國內外學者的深入研究[13],已在無線電協作頻譜感知[14]、容量約束車輛路徑(CVRP)問題求解[15]、作業車間調度優化[16]、水資源網絡優化、PID 控制器參數優化等領域得到成功應用。

本文從提高水文模型參數的優化精度和優化算法的運行效率出發,利用混合蛙跳算法對水文模型參數進行優化,將混合蛙跳算法與新安江模型耦合,將算法與應用實例結合來分析和驗證算法,從而將混合蛙跳算法引入到水文模型參數估計中,為模型的參數估計提供新途徑,從而提高流域水文模型的模擬精度和運行效率。

1新安江模型

新安江模型是具有多參數,能適應我國的濕潤、半濕潤地區概念型流域的水文模型。該模型先將大的流域劃分成多個較小的單元流域,然后對單元流域進行產匯流計算,從而得到各個單元流域的出口流量數據;接下來進行出口以下的河道洪水演算,從而求出相應單元流域的出流過程;最后將各個單元流域的出流過程相加,得到該流域出口的總出流過程[17]。

本文主要基于二水源的新安江模型,該模型的結構如圖1所示。在模型中,為了方便地計算出流域的蒸散發量,將模型中的包氣帶土壤層分為3層,包括上層WU,下層WL和深層WD。模型利用流域土壤穩定下滲率FC將流域中的總徑流量劃分為地表產流RS和地下產流RG兩個部分。如表1所示,修改后的新安江流域水文模型共包含了10個參數[18]。表中列出了各個參數的名稱、意義以及在實際實驗流域內的參數取值范圍[18]。

圖1 修改后的新安江模型結構

參數意義單位取值范圍K流域蒸散發能力和實測水面蒸發比值-0.001~1.000IMP流域內不透水面積占全流域面積比值-0.001~0.500B蓄水容量曲線的方次-0.001~1.000WUM流域中的上層土壤蓄水容量mm5.000~30.000WLM流域中的下層土壤蓄水容量mm50.000~100.000WDM流域中的深層土壤蓄水容量mm50.000~200.000C流域中的下層和深層土壤蒸散發系數-0.001~0.300FC穩定入滲率mm/h0.001~50.000KKG基流匯流所用線性水庫的退水系數-0.001~0.990Kr表面產流的匯流參數-0.001~10.000

2混合蛙跳算法

Eusuff等在2003年基于粒子群算法(PSO)提出了混合蛙跳算法(SFLA)。混合蛙跳算法是一種新型的、基于群智能的啟發式優化方法。該方法主要結合了粒子群(PSO)優化算法和模因Memetic算法的優點,具有概念簡單、搜索能力較強、計算速度較快以及較易實現等優點。該算法將每只青蛙作為一個模因載體,然后通過不同的青蛙個體之間進行交流信息,從而促進和實現各個模因進化,解決最優化問題。全局信息交換和局部搜索的平衡策略使得混合蛙跳算法能較快地收斂,并獲得全局最優解[10]。SFLA優化算法的流程如圖2所示。

圖2 混合蛙跳算法流程

混合蛙跳算法主要分為以下步驟[10]:

步驟1初始化算法參數。青蛙種群大小設置為P,搜索空間維數設為S,模因組數設為M,每個模因組包含N只蛙,則種群大小P=M×N。算法局部搜索最大次數設為T1,全局最大混合迭代次數設為T2。

步驟2生成蛙群。隨機生成P=M×N只青蛙的初始群體,記為XP={xi|xi=(xi1,xi2,…,xis),i=1,2,…,P};然后將每只青蛙作為適應度函數的輸入變量,計算出適應度值f(xi);接下來依據適應度值對蛙群按降序方式進行排序,找到并記錄下蛙群中最優蛙的信息。

步驟3生成模因組。將青蛙按照逐一循環分配的方式分配至M個模因組。每個模因組中將包含N只蛙,然后比較排序和記錄每個模因組中的最優蛙和最差蛙的信息。

步驟4進化模因。對每個模因組中適應度值最差的個體x(k)w按公式(1)進行局部搜索和位置更新工作,并利用適應度函數計算其適應度值,然后對每一個模因組重新排隊尋找最差個體。如果前后兩次的進化并非是對同一青蛙的操作,那說明前一次的進化是有效的,從而使前一次的最差個體得到進化。如此進行多次迭代,便能使整個模因組朝著確切方向進化,從而使整個群體的青蛙所代表的解的質量得到改善;如果通過模因組的最優個體和種群的最優個體都不能使當前個體得到進化,則當前個體在參數的全局空間內隨機搜索,以拓展算法的全局搜索能力。重復上述局部搜索過程共T1次,得到進化后的模因組。

(1)

步驟5混合模因組。將進化后的模因組內的蛙重新混合,計算每個青蛙的適應度值,然后按照適應度值遞增的方式進行青蛙排序,并根據排序結果更新蛙群中的最優蛙。

步驟6判斷算法終止條件。判斷青蛙族群中全局最優蛙是否達到收斂精度或者全局進化次數是否達到T2。若未達到終止條件,則轉至步驟3,青蛙群體再次進行全局信息交換;否則,停止混合蛙跳算法的運行,并輸出全局最優蛙的信息,即獲得的優化問題的最優解。

3SFLA與新安江模型應用模型

基于混合蛙跳優化算法實現對新安江模型參數的優化應用,模型如圖3所示。

圖3 新安江模型參數估計的應用模型

該應用模型實現了混合蛙跳優化算法與新安江水文模型的耦合,通過混合蛙跳算法產生優化后的模型參數,然后在模型中利用目標函數來調用新安江水文模型進行模擬計算,得出該優化參數在模型模擬中的目標函數值;然后混合蛙跳算法將該反饋值作為算法優化改進的信息,經過大量的重復迭代工作使得算法進化,從而找到符合精度的參數優化結果(全局最優解)或達到最大進化次數。圖3中的破折號框包括了從數據庫或數據文件中獲取到的基礎數據。

在應用混合蛙跳算法優化新安江模型參數時,利用目標函數判斷新個體的優良程度。首先以新的個體作為新安江模型的參數,新安江模型利用此參數和驅動數據,模擬產生模擬數據,將模擬數據與真實的觀測數據對比,統計NSE值反映模擬量與觀測量之間的逼近程度,當NSE的值越大,說明逼近程度越好,意味著個體所反映的參數越優秀;反之,參數越差。通過算法的不斷迭代,找到能夠使模型的模擬量與觀測值逼近程度更高的參數。NSE表達式如下:

(2)

4實例應用及分析

4.1研究區域

本文的研究區域選用黑河流域中的扎馬什克水文站以上流域。扎馬什克水文站海拔高度為2 635 m,流域控制面積約為4 589 km2。在本文的模型參數優化方法的實例應用實驗中,共收集了扎馬什克水文站1990—2000年間逐日流量資料以及野牛溝氣象站在1990—2000年(共11年間)的每6 h/d的氣象數據(降水、氣溫、氣壓、風速和水面蒸發等)[18]。

4.2實驗說明和參數設置

在內存為4G,處理器為Intel(R)Core i5-3750 3.40 GHz的計算機上,采用VC++6.0進行實驗。在實驗過程中,為了比較混合蛙跳算法在優化新安江模型時的算法性能,與應用于新安江模型的參數優化的遺傳算法(GA)[18]進行比較。在此,簡稱混合蛙跳算法和遺傳算法優化新安江模型分別為SFLAXAJ和GAXAJ。對SFLAXAJ和GAXAJ可以從以下兩個方面進行實驗比較:① 固定算法迭代次數,比較2種模型在達到一定迭代次數后的目標函數值的大小;② 固定優化精度,即給定目標函數值,比較2種算法收斂至特定目標函數值所需要的迭代次數。通過這2個不同角度的對比,實現對2種優化模型的收斂能力和速度的比較。下面對SFLAXAJ和GAXAJ2種模型的參數進行設置。

SFLAXAJ:群體大小SN=100,模因組個數M=10,每個模因組包含的個體數N=10,局部迭代次數I=10;

GAXAJ:群體大小SN=100,變異概率Pb=0.1,雜交概率Pz=0.6,復制概率Pf=0.3。

4.3實驗結果與分析

4.3.1固定模型迭代次數

設置SFLAXAJ和GAXAJ兩種模型的全局迭代次數G=1 000。當2種模型全局迭代次數達到1 000次時,統計2種模型的實驗結果,主要包括10個參數結果以及目標函數NSE值,NSE值越大,說明獲得的結果越好,反之越差。為了更加客觀地說明2種模型的實驗結果,使每個模型分別獨立運行30次,分析2種模型的工作性能,統計2種模型在循環30次過程中所獲得的目標函數最優值、最差值、平均值以及30次結果的方差,統計結果如表2所示。在30次運行過程中,2種模型獲得的最優參數結果和目標函數值如表3所示。圖4顯示了2種模型運行30次的平均收斂過程。

表2 兩種模型循環30次獲得的目標

表3 兩種模型循環30次獲得的最優參數結果和目標函數值

圖4 兩種模型的平均收斂過程

從實驗統計的結果可以看出:SFLAXAJ模型比GAXAJ模型具有更優的結果,在30次運行結果中,SFLAXAJ模型比GAXAJ模型的平均優化精度提高了2.5%。同時,SFLAXAJ的方差值比GAXAJ要小很多,表明SFLAXAJ模型運行更加穩定。通過結果對比,說明SFLAXAJ模型具有比GAXAJ模型更加優秀的全局收斂能力和更高的收斂精度。這是由于在GAXAJ模型中,部分較為優秀的個體會以較大的概率作為父代個體進化計算產生子代個體,在產生新群體過程中,缺乏對整個舊群體信息的充分利用,經過多次迭代,將會使子代群體快速地向同一區域聚集。在GAXAJ模型中,變異算子的加入在一定程度上能增強算法的全局收斂能力,但由于新安江水文模型較為復雜,變異算子使群體跳出局部極值區域的能力相對較弱,導致GAXAJ模型出現“早熟”現象。在SFLAXAJ模型中,每個模因組中的最壞個體以較小的步長逐步向該模因組中的最優個體或者全局最優個體迭代靠攏,或者是在最壞個體的解得不到改善的情況下,將會在給定取值范圍內隨機生成新的個體。這種進化機制不會使所有粒子迅速地向同一極值區域迅速收斂,從而增強了每個粒子的搜索范圍。同時在所有模因組完成局部迭代后,將會對所有模因組混合按照所有個體的目標函數值降序排序重新劃分模因組,使得每個模因組中的所有個體的覆蓋空間有所增大,增強了全局收斂能力。通過多個模因組的共同組合,大大提升了SFLAXAJ模型的全局收斂能力和收斂精度。

4.3.2固定優化精度

設定目標函數值NSE=0.75,比較SFLAXAJ模型和GAXAJ模型收斂至該優化精度值時所需要的迭代次數。每個模型分別獨立運行30次,比較兩種模型優化結果的成功率(當迭代次數大于 1 000 仍未收斂至特定目標函數值時,認為模型運行失敗)、最優值、最差值、平均值以及30次結果的方差,統計結果如表4所示。

表4 兩種模型循環30次獲得的目標

從實驗結果可以看出:在固定優化精度時,SFLAXAJ模型在30次運行過程中的成功率為100%,而GAXAJ的成功率不到50%,SFLAXAJ模型比GAXAJ模型的平均成功率提高了53.33%。在平均迭代次數方面,SFLAXAJ模型比GAXAJ模型低很多。這是由于在GAXAJ模型中,群體的進化缺乏對最優秀的個體的充分利用,同時,在父代個體進化計算產生子代個體時,缺乏對父代個體和子代個體優劣程度的對比,很難保證新的群體整體解的質量一定優于舊的群體解的質量,增加了群體進化的不確定性,導致模型的收斂精度和收斂速度不理想。而在SFLAXAJ模型中,群體中最為優秀的個體往往攜帶了重要的群體進化的啟發信息,每個模因組中的最壞個體進化都是向模因組中的最優個體或者全局最優個體方向的進化搜索,充分利用了精英個體所攜帶的啟發信息。即便是在最壞個體解的質量得不到改善的前提下,仍將會在給定空間內隨機生成新的個體,擴展了群體的搜索空間,使模型的局部搜索能力和全局搜索能力得到了很好的結合。同時,在SFLAXAJ模型中,除了每個模因組最壞個體的搜索,其他個體并不參與進化,這樣能保證除群體部分最壞個體外,其他個體并不會出現退化現象,使群體的進化朝著確定方向收斂,提高了模型的收斂精度和收斂速度。

4.3.3實驗總結

對SFLAXAJ模型進行實驗仿真,同時與文獻[18]中的GAXAJ模型進行實驗對比。通過對實驗結果的定量分析,得出SFLAXAJ模型無論是在全局收斂能力,還是收斂精度與速度方面,都比GAXAJ模型更具優勢。

5結束語

水文模型參數優化對水文模型的成功應用至關重要,但水文模型本身較為復雜,造成了水文模型參數估計較為困難,用傳統的辦法很難獲得較優秀的參數結果。本文提出了基于混合蛙跳算法的水文模型參數方法,并將該算法與具體的新安江模型進行耦合,選特定區域進行實驗仿真。通過與基于遺傳算法的水文模型參數估計方法定量對比分析可知:基于混合蛙跳算法的水文模型參數估計更具優勢,全局收斂能力更好,收斂速度更快,證明了混合蛙跳算法可以應用于水文模型的參數估計。

參考文獻:

[1]FRANCHINI M,GALEATI G.Comparing several genetic algorithm schemes for the calibration of conceptual rainfall-runoff models[J].Hydrological Sciences Journal,1997,42(3):357-379.

[2]VRUGT J A,GUPTA H V,BOUTEN W.Real-time data assimilation for operational ensemble streamflow forecasting[J].J.Hydrometeorol, 2006(7):548-564.

[3]HAISHEN LU,TING HOU,HORTON R,et al.The streamflow estimation using the Xinanjiang rainfall runoff model and dual state-parameter estimation method[J].Journal of Hydrology,2013,480:102-114.

[4]ZHAO R J,ZHANG Y L,FANG L R,et al.The Xinanjiang model[C]//Hydrological Forecasting Symposium:Proceedings of the Oxford Symposium.Wallingford:IAHS Press,1980.

[5]ZHAO R J.The Xinanjiang model applied in China[J].J.Hydrol.1992,135:371-381.

[6]江燕,劉昌明,胡軼松,等.新安江模型參數優選的改進粒子群算法[J].水利學報,2007,39(10):1200-1206.

[7]CHENG C T,OU C P,CHAU K W.Combining a fuzzy optimal model with a genetic algorithm to solve multi-objective rainfall-runoff model calibration[J].Journal of Hydrology,2002,268(1):72-86.

[8]李致家,王壽輝,HAPUARACHCHI H A P.SCE-UA方法在新安江模型參數優化中的應用[J].湖泊科學,2001 (4):304-314.

[9]孟新華,涂啟玉,周年華,等.基于遺傳模擬退火算法的新安江模型參數優選[J].水電自動化與大壩監測,2009 (3):64-67.

[10]EUSUFF M M,LANSEY K E.Optimization of water distribution network design using the shuffled frog leaping algorithm[J].J of Water Resources Planning and Management,2003,129(3):210-225.

[11]崔文華,劉曉冰,王偉,等.混合蛙跳算法研究綜述[J].控制與決策,2012,27(4):481- 486.

[12]駱劍平,李霞,陳泯融.混合蛙跳算法的Markov 模型及其收斂性分析[J].電子學報,2010,38(12):2875- 2880.

[13]劉立群,王聯國,韓俊英,等.基于全局共享因子的混合蛙跳算法[J].計算機工程,2013,30(10):162-155.

[14]鄭仕鏈,樓才義,楊小牛.基于改進混合蛙跳算法的認知無線電協作頻譜感知[J].物理學報,2010,59(5):3611-3617.

[15]駱劍平,李霞,陳泯融.基于改進混合蛙跳算法的CVRP求解[J].電子與信息學報,2011,33(2):429- 434.

[16]蔡良偉,李霞.基于混合蛙跳算法的作業車間調度優化[J].深圳大學學報:理工版,2010,27(4):391-395.

[17]趙人俊.流域水文模擬—新安江模型和陜北模型[M].北京:水利水電出版社,1984:98-142.

[18]王書功.水文模型參數估計方法及參數估計不確定性研究[M].鄭州:黃河水利出版社,2010:5-142.

(責任編輯楊黎麗)

Research of Parameters Estimation Method for Hydrological Model Based on Shuffled Frog Leaping Algorithm

HUO Jiu-yuan1, 2, LIU Li-qun3, ZHAO Hong-xing1

(1.School of Electronic and Information Engineering, Lanzhou Jiaotong University,Lanzhou 730070, China;2.Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences, Lanzhou 730000, China;3.Information Institute of Science and Technology,Gansu Agricultural University, Lanzhou 730070, China)

Abstract:Traditional heuristic algorithms usually have defects such as prematurity and slow convergence for parameter estimation in Hydrologic model. In order to improve the accuracy of parameter optimization and performance of optimization algorithm for hydrological model, Shuffled Frog Leaping Algorithm (SFLA) was introduced to propose a new method of Hydrological model parameters estimation based on SFLA. Then the method was applied to estimate the parameters of Xin’anjiang model, and compared with the parameter estimation method based on Genetic algorithm (GA) through experimental analysis. Experimental results show that the parameter optimization method based on SFLA could improve 2.5% of the average optimize accuracy and to GA, and 53.33% of success optimize rate under the fixed precision condition than the GA method. It is proved that SFLA could be applied for parameter estimation in Hydrologic model and it has more obvious advantages in both the speed and the accuracy of convergence.

Key words:shuffled frog leaping algorithm; Xin’anjiang model; parameter estimation; Nash-Sutcliffe efficiency coefficient

中圖分類號:TP3

文獻標識碼:A 1674-8425(2016)03-0080-07

doi:10.3969/j.issn.1674-8425(z).2016.03.014

作者簡介:火久元(1978—),男,甘肅永登人,博士,副教授,主要從事智能計算、模型參數優化研究。

基金項目:國家自然科學基金資助項目(61462058);蘭州交通大學青年科學基金資助項目(2013032)

收稿日期:2015-12-10

引用格式:火久元,劉立群,趙紅星.基于混合蛙跳算法的水文模型參數估計方法[J].重慶理工大學學報(自然科學),2016(3):80-86.

Citation format:HUO Jiu-yuan, LIU Li-qun, ZHAO Hong-xing.Research of Parameters Estimation Method for Hydrological Model Based on Shuffled Frog Leaping Algorithm[J].Journal of Chongqing University of Technology(Natural Science),2016(3):80-86.

猜你喜歡
優化模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: hezyo加勒比一区二区三区| 久热99这里只有精品视频6| 亚洲天堂区| 国产va在线观看| 久久天天躁狠狠躁夜夜躁| 97se亚洲综合不卡 | 少妇被粗大的猛烈进出免费视频| 热re99久久精品国99热| 国产杨幂丝袜av在线播放| 亚洲色图欧美激情| 日韩精品一区二区深田咏美| 91丝袜乱伦| 国产在线自乱拍播放| 国内精品免费| 国内丰满少妇猛烈精品播| 国产日韩欧美一区二区三区在线| 欧美成人影院亚洲综合图| 久久国产乱子伦视频无卡顿| 青青青视频免费一区二区| 亚洲丝袜第一页| 亚洲欧美日本国产专区一区| 亚洲欧美自拍一区| 午夜精品影院| 999福利激情视频| 国内精自线i品一区202| 伊人久久大香线蕉aⅴ色| 欧美成人午夜影院| 亚洲精品无码抽插日韩| 国产精品乱偷免费视频| 成人在线亚洲| 国产日产欧美精品| 免费一极毛片| 91精品在线视频观看| 一本久道久久综合多人| 成色7777精品在线| 伊人色综合久久天天| 人人看人人鲁狠狠高清| 毛片大全免费观看| 国产麻豆精品手机在线观看| 亚洲综合色吧| 亚洲一区色| 日韩福利在线视频| 日韩 欧美 小说 综合网 另类 | 国产亚洲精品在天天在线麻豆| 99精品视频九九精品| 成人午夜久久| 久久久久亚洲精品成人网 | 日韩精品中文字幕一区三区| 精品无码日韩国产不卡av| 亚洲人成高清| 欧洲亚洲欧美国产日本高清| 黄色国产在线| 1级黄色毛片| yy6080理论大片一级久久| 91色在线视频| 亚洲人成日本在线观看| 色综合久久88| 一级做a爰片久久毛片毛片| 东京热高清无码精品| 成人福利在线免费观看| 亚洲精品第1页| 国产哺乳奶水91在线播放| 亚洲成人精品| 天天躁夜夜躁狠狠躁躁88| 国产不卡网| 日韩色图在线观看| 青青国产视频| 伊人色在线视频| 久久综合结合久久狠狠狠97色 | 91欧美在线| 国产午夜无码专区喷水| 国产精品视频免费网站| 试看120秒男女啪啪免费| 亚洲精品不卡午夜精品| 欧美精品v| 日韩高清一区 | 久久久久国产精品嫩草影院| 免费国产高清视频| 欧美中文字幕第一页线路一| 国内精品手机在线观看视频| 国产乱肥老妇精品视频| 91啪在线|