亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        面向水文模擬的大規(guī)模多級并行參數(shù)率定框架

        2021-09-13 02:13:12全婷李強聶寧明田在榮
        青島大學學報(自然科學版) 2021年3期
        關鍵詞:優(yōu)化算法

        全婷 李強 聶寧明 田在榮

        摘要:為了實現(xiàn)大尺度水文模擬中的參數(shù)率定,提出一種基于優(yōu)化算法的大規(guī)模多級并行參數(shù)率定框架。首先利用MPI劃分子通信域的技術,實現(xiàn)了多級并行處理框架,其次設計了基于對等模式的整體架構,以充分利用處理器資源,最后使用大量非阻塞式通信的方式優(yōu)化了計算效率,減少了進程間等待。將該框架應用于HIMS水文模型的參數(shù)率定,試驗結果表明,對等多級并行框架相對于主從并行框架具有更好的尋優(yōu)效果,利用非阻塞式通信,在尋優(yōu)效率上有所提升。該框架能夠高效地利用大規(guī)模處理器且有效地縮短運行時間,提升了參數(shù)優(yōu)化的整體效率,具有良好的擴展性。

        關鍵詞:水文模擬;大規(guī)模參數(shù)率定;對等式多級并行框架;非阻塞式通信;優(yōu)化算法

        中圖分類號:TP338.6 文獻標志碼:A

        文章編號:1006-1037(2021)03-0014-08

        隨著GIS/RS等高新技術的發(fā)展,分布式水文模型大量出現(xiàn),如英國、丹麥和法國的水文工作者在1986年合作研發(fā)的MIKE-SHE模型[1]。此后,SWAT、IHDM等水文模型的出現(xiàn),也驗證了分布式大規(guī)模水文模擬實現(xiàn)的可行性。雖然國內水文模擬的發(fā)展起步較晚,但也涌現(xiàn)了HIMS(Hydroinformatic Modeling System)[2]、GBHM(Geomorphology-Based Hydrological Mode)[3]、新安江模型[4]以及EasyDHM(Easy Distributed Hydrological Model)[5]等分布式水文模型,在各流域模擬上展現(xiàn)了出色的模擬效果。分布式水文模型一般包含多個參數(shù),其中大部分參數(shù)有明確的物理含義,能夠根據(jù)實測的數(shù)據(jù)進行確定,但是在實際中仍然有很多參數(shù)需要通過率定獲得。由于分布式水文模型需要率定的參數(shù)空間維度通常十幾、二十維,高維參數(shù)空間搜索計算量巨大,傳統(tǒng)搜索方式已經難以實施。20世紀80年代初興起的啟發(fā)式算法,包括遺傳算法、模擬退火算法、粒子群算法和人工神經網(wǎng)絡等,極大地推動了分布式水文模型中參數(shù)優(yōu)化的發(fā)展。Wang[6]較早的將GA方法應用于徑流模型的產流參數(shù)率定;Cheng等[7]將GA與模糊優(yōu)選原理相結合,將其用于新安江三水源模型的參數(shù)率定;馬海波等[8]將SCE-UA算法應用到TOPMODEL模型的參數(shù)率定;Thiemann等[9]以美國密西西比河Leaf流域(19 944 km2)為例,將貝葉斯回歸參數(shù)估計算法應用于Nash-Cascade模型和Sacramento模型,驗證了貝葉斯算法在參數(shù)優(yōu)化方面有很好的效果。但是由于傳統(tǒng)參數(shù)優(yōu)化過程都是串行執(zhí)行,運行時間較長。特別在大尺度水文模擬的參數(shù)率定中,運行時間有可能長達幾個小時到幾個月不等,率定優(yōu)化效率問題成為制約水文模型發(fā)展的重要問題。近些年,隨著高性能計算技術的飛速發(fā)展,在科學研究的諸多領域,并行計算都顯著地提高了大規(guī)模數(shù)據(jù)處理的效率,并且可以通過算法的設計來提高最終結果的可行性。因此,越來越多的科研工作者將此技術應用于水文模擬中。Kollet等[10]證明了將并行計算用于水文模擬正確性與可行性;Yalew等[11]還探索了將并行計算用于SWAT中的水文模擬過程;Li等[12]運用分解法和MasterSlave范式開發(fā)了動態(tài)并行水文模擬模型。通過水文模擬和高性能計算的結合,水文模型的計算效率有了相應的提高。但是這些研究都是針對水文模擬過程的加速,而花費時間最長的參數(shù)率定過程的求解加速并未直接涉及。當前針對參數(shù)率定優(yōu)化加速方法的研究主要以優(yōu)化算法的改進為主,如Huo等[13]將改進的多核并行人工蜂群優(yōu)化算法應用于水文模擬參數(shù)優(yōu)化過程中;闞光遠等[14]對SCE-UA算法并行化;李強等[15]設計的基于克里金插值理論的對等式并行框架。而以參數(shù)率定并行框架為主的研究并不多見,如李強等[16]設計的應用于SWAT模型的參數(shù)敏感度分析的主從模式并行框架;申蒙蒙等[17]通過改進優(yōu)化算法中的相關算子實現(xiàn)的雙層主-從并行框架等。這些并行框架的加入雖然使得優(yōu)化算法的搜索效率和運行效率相對串行有了很大的提升,但這些率定框架的并行擴展性還有待提高。在現(xiàn)有水文模擬參數(shù)率定并行實現(xiàn)中,參數(shù)率定主從模式并行框架可擴展規(guī)模的的數(shù)量級皆在百核數(shù)左右。這是由于主從模式中的主進程和從進程的頻繁數(shù)據(jù)交換和大量從進程的請求會導致從進程等待,從而影響最終的運行效率和并行規(guī)模的擴大。如果想要實現(xiàn)萬核數(shù)甚至數(shù)十萬核數(shù)的大規(guī)模水文模擬參數(shù)率定,需要突破當前已有的率定框架設計思路,從應用本身所具有的特點出發(fā),設計真正適用于水文模擬參數(shù)率定的大規(guī)模率定方法。針對以上問題,本文提出一套針對大規(guī)模水文模擬參數(shù)率定的基于非阻塞通信的對等式多級并行框架。通過劃分子通信域的方法,將框架分成多級并行結構。由于框架的特殊性,在進行諸如粒子群優(yōu)化的過程中,采用異步非阻塞通信技術,使得群體內部的信息能夠得到分享,從而進一步提高框架的搜索效率。

        1 參數(shù)率定

        在采用水文模型進行徑流模擬和預報時,關鍵在于尋求水文模型的最優(yōu)參數(shù)值。水文模型的參數(shù)優(yōu)化過程主要分為兩部分:模型率定和模型校驗,本文主要研究模型率定。在模型率定時,將實測降雨、蒸發(fā)、徑流等資料輸入到水文模型中,以模型的模擬流量和實測流量的誤差最小為目標函數(shù),采用優(yōu)化算法來搜尋模型參數(shù)的最優(yōu)值。選擇NS效率系數(shù)[18]來驗證水文模型模擬徑流的精確度

        其中,NS為納什效率系數(shù),yc為預報值,m3/s;yo為實測值,m3/s;o為實測值的平均值,m3/s;n為系列長度。NS越接近1,該水文模型的質量越好,可信度更高;NS接近0,該水文模型的模擬結果接近觀測值的平均值水平,即總體結果可信,但是誤差還是較大;NS遠遠小于0,則該水文模型不可信。

        參數(shù)率定的主要計算過程是根據(jù)率定優(yōu)化算法選取步長和方向,得到新的粒子點,每個粒子代表一組參數(shù);本文中的水文模型寫成了一個函數(shù)接口,通過一組參數(shù)來調用水文模型函數(shù)接口進行水循環(huán)模擬,得到徑流量模擬值,最后使用式(1)得到NS效率系數(shù)。整個率定過程可以作為一個子任務,分配到多個進程中,利用多個進程并行計算,最終得到最優(yōu)的NS效率系數(shù)?;趨?shù)優(yōu)化算法特點,如若想要更加均衡的使用更多的處理器來進一步提高搜索效率,則可以通過MPI劃分子通信域的方式來加大計算規(guī)模。

        2 多級并行參數(shù)率定框架

        2.1 基于子通信域的多級并行框架

        為了更好地提高計算效率,本文選擇以進程級并行作為主要并行策略。MPI(Message Passing Interface)[19]是一個性能很好的并行程序接口。具有兩種最基本的并行程序設計模式:對等模式和主從模式[20]。由于對結果的收集需要進行大量的規(guī)約和廣播操作,同時參數(shù)率定過程可以獨立為一個模塊,本文并行設計模式采用的是對等模式,更好的適應大規(guī)模測試。對等模式設計如圖1所示。本文程序的基礎并行框架采用對等模式,但為了率定優(yōu)化需要,邏輯上分為三級:主管理進程、子管理進程和計算進程。為了提高框架的并行效率,同時為了適應后續(xù)工作的異構結構特征設計,本文考慮劃分子通信域。

        算法設計如下:將整個通信域MPI_COMM_WORLD中所有進程創(chuàng)建進程組MPI_GROUP_WORLD;排除主進程獲取新的進程組newgrp;將新進程組newgrp創(chuàng)建為一個通信域newcomm;獲取通信域newcomm中進程的進程號newrank和進程數(shù)newnum_procs;將通信域newcomm劃分為subSpace個子通信域subcomm;將通信域MPI_COMM_WORLD的前subSpace+1個進程劃為組group,將組group創(chuàng)建為通信域comm。其中,通信域comm中的0號進程為主管理進程,子通信域subcomm中的0號進程為子管理進程,子通信域subcomm中的其他進程為計算進程,至此框架的三級劃分完成,如圖2所示。

        2.2 多級并行框架流程

        采用MPI劃分子通信域實現(xiàn)的三級并行框架,其中的主管理進程、子管理進程和計算進程之間的關系如圖3所示。主管理進程和子管理進程劃分在一個通信域comm中,子管理進程和計算進程劃分在多個子通信域subcomm中。首先主管理進程進行一個非阻塞式接收,用于接收來自子管理進程的收斂信號,子管理進程進行一個非阻塞式廣播和非阻塞式接收,用于接收來自主管理進程的退出信號和計算進程的收斂信號,計算進程進行一個非阻塞式廣播,用于接收來自子管理進程的退出信號。然后進入?yún)?shù)率定,獲得各自最優(yōu)解,同時對得到的最優(yōu)解進行收斂性分析。主管理進程在進行收斂性分析時,若自身收斂或者子管理進程也收斂,則向所有其他子管理進程發(fā)送結束信號;若二者都未收斂,則回到參數(shù)率定過程,以此迭代,直至找到最優(yōu)解。子管理進程在進行收斂性分析時,若自身收斂,則同時向主管理進行非阻塞式發(fā)送以及向計算進程進行非阻塞式廣播,通知它們退出;考慮同時收到主管理進程和計算進程的收斂信號或者只收到其中一個收斂信號的可能性,進行最優(yōu)解比較和非阻塞式通信,通知對應的進程退出;若本進程未收斂,同時也未收到收斂信號,則回到參數(shù)率定過程,以此迭代,直至找到最優(yōu)解。計算進程在進行收斂性分析時,若自身收斂,檢查是否收到子管理的收斂信號,若收到,則比較最優(yōu)解,否則直接向子管理進程進行發(fā)送結束信號;若本進程未收斂,同時也未收到收斂信號,則回到參數(shù)率定過程,以此迭代,直至找到最優(yōu)解。

        2.3 面向并行框架的非阻塞通信優(yōu)化

        本文并行框架中采用自定義的規(guī)約操作,調用MPI_Reduce( )函數(shù)接口對主管理進程和子管理進程中的最優(yōu)解進行收集和比較,使得框架整體所有進程搜索得到的最優(yōu)解能夠被迅速規(guī)約并將當前全局最優(yōu)解分享到每一個進程。該框架可與率定算法緊密結合,框架中計算進程可以映射成率定算法的搜索粒子點,框架中子通信域可以映射成率定算法的種群,可以與粒子群算法、遺傳算法等無縫結合。本文程序中主管理進程采用的率定算法為遺傳算法[21],子管理進程和計算進程采用的率定算法為粒子群算法[22]。結合該特性,本文在每個子通信域中的子管理進程和計算進程之間加入了非阻塞式通信,計算進程異步地發(fā)送當前最優(yōu)值和對應參數(shù)組到子管理進程,同時異步的接收子管理進程發(fā)送過來的整個子通信域的最優(yōu)值和對應參數(shù)組,子管理進程異步的循環(huán)接收來自計算進程的最優(yōu)值和對應參數(shù),比較后,將該子通信域的最優(yōu)值和對應參數(shù)發(fā)送給計算進程。由于非阻塞式通信主要用于計算和通信的重疊,通信操作全部后臺運行,從而提高整個程序的執(zhí)行效率。

        3 數(shù)值試驗

        將本文的框架應用在HIMS水文模型參數(shù)率定模塊的數(shù)值實驗中,實驗主要采用拉薩河流域數(shù)據(jù)針對HIMS模型的12個參數(shù)進行參數(shù)率定。通過比較和分析數(shù)值實驗結果,來探討框架效率以及并行性能。測試計算環(huán)境為中國科學院計算機網(wǎng)絡信息中心超級計算機“元”和國產先進計算系統(tǒng)。

        3.1 參數(shù)尋優(yōu)效果的驗證

        首先利用HIMS水文模型對拉薩河流域進行參數(shù)率定優(yōu)化的計算結果分析該并行框架的尋優(yōu)效果和率定算法的正確性,尋優(yōu)效果將從以下四個方面分析。

        (1) 普通優(yōu)化算法和并行混合算法框架的尋優(yōu)效果對比。本次實驗將HIMS水文模型中參數(shù)環(huán)境設置為粒子數(shù)20個,模擬500次,結果如圖4所示。并行混合算法優(yōu)化的NS效率系數(shù)在同樣的迭代次數(shù)下大于其他單獨的優(yōu)化算法的計算結果,說明參數(shù)率定并行框架有相對突出的優(yōu)化效果,達到了水文模型對流域進行參數(shù)率定優(yōu)化的所需達到的基本要求。

        (2) 主從模式和對等模式的尋優(yōu)效果對比。在兩種模式框架對比中,采用三種規(guī)模進行數(shù)值測試,分別為:小規(guī)模21核,粒子數(shù)200個;中規(guī)模201核,粒子數(shù)2 010個;大規(guī)模801核,粒子數(shù)8 010個。通過NS效率系數(shù)和運行時間對比兩者的尋優(yōu)效果,結果如圖5所示。每種情況測試10次取均值。可以看出,在率定精確性方面,對等模式相對主從模式略占優(yōu)勢,對等模式比主從模式模擬的NS效率系數(shù)高0.01左右。但從計算效率來看,對等模式在大規(guī)模率定過程中的運行時間相對主從模式的運行時間要短很多,對等模式在大規(guī)模參數(shù)率定中發(fā)揮較好的效果。

        (3) 加非阻塞通信和不加非阻塞通信的尋優(yōu)效果對比。對比加入非阻塞通信和不加入非阻塞通信的兩種情況下的NS效率系數(shù)(圖6),同時測試了迭代50次(圖6(a))和迭代500次(圖6(b))的情況。由圖6(a)可知,在固定迭代50次后,加入非阻塞通信的NS效率系數(shù)比未加入非阻塞通信的NS效率系數(shù)高一些,同時隨著進程數(shù)的增加,NS效率系數(shù)變高;由圖6(b)可知,在固定迭代500次后,加入非阻塞通信的NS效率系數(shù)比未加入非阻塞通信的NS效率系數(shù)稍高一點,但總體上相差不大。以上分析說明在迭代次數(shù)較少的情況下,加入非阻塞通信的框架尋優(yōu)效果較好,迭代次數(shù)增加后,兩者之間差距縮小,說明加入非阻塞通信能夠加快算法收斂,從而整體更快的貼近最優(yōu)值。

        (4) 算法正確性驗證。通過對比率定過程中的徑流量和實測徑流量來驗證算法正確性,結果如圖7所示。伴隨著迭代次數(shù)的增加,水文模擬出來的徑流量總體趨勢上越來越貼近實際測量值,說明通過率定算法的參數(shù)集越來越靠近最優(yōu)參數(shù)集,通過這組最優(yōu)參數(shù)集水文模擬得到的徑流量也會與實測值更加貼近吻合。

        3.2 多級并行參數(shù)率定框架的擴展性測試

        3.2.1 主從模式和對等模式可擴展性對比 針對主從模式框架和對等模式框架在強擴展性方面進行測試。模擬總粒子數(shù)為9 000粒子,隨著進程的增加,每個進程中的粒子數(shù)減少。其中,由于主從框架的特性,在數(shù)據(jù)整理時,總進程中去除了前四個只進行優(yōu)化算法的固定進程。兩者對應的加速比對比如圖8所示??梢钥闯?,隨著進程數(shù)的增加,對等模式框架的加速比呈指數(shù)上漲,呈現(xiàn)較好的擴展性,而主從模式框架的加速比在90個進程之后幾乎沒有變化,甚至在900進程處呈下降趨勢??梢?,對等模式框架在大規(guī)模率定中的擴展性遠遠好于主從框架。

        3.2.2 對等模式在大規(guī)模方面的擴展性測試 (1) 弱擴展性測試。在多級并行參數(shù)率定框架的并行弱擴展測試中,在保持每個進程所分到的計算量基本一致的情況下,通過增加進程數(shù)來測試程序的并行性能。測試以6 251核為基準,每種情況測試10次取均值,迭代次數(shù)為500次。所有算例均設置每個子通信域50個進程,每個進程1 600個粒子,保證每個進程分到的任務量基本一致。測試的運行時間和并行效率如圖9所示。從測試結果可以看到,以6 251核為基準,到十萬核時并行效率還保持在86%,整個框架始終運行在較高的效率之上,擴展性較好。

        (2) 強擴展性測試。參數(shù)率定框架的強擴展測試實驗中,在保證問題計算規(guī)模一定的情況下,測試程序隨著核數(shù)增多而導致程序性能的變化情況。本文測試以6 251核為基準,將核數(shù)翻倍,總粒子數(shù)在1.6億左右,迭代次數(shù)500次。測試的運行時間和加速比如圖10所示,測試的并行效率如圖11所示??梢钥闯觯瑥姅U展性測試具有較為理想的加速比與并行效率,在十萬核之內保持85%的并行效率,說明劃分平衡度較好,且非阻塞式通信使得通信與計算重疊,由于通信產生的核間等待很少??傮w擴展性良好。

        4 結論

        本文采用了MPI中的非阻塞式通信實現(xiàn)了一個多級并行框架,同時將多種優(yōu)化算法相結合,用于水文模型中的參數(shù)優(yōu)化模塊。分析HIMS模型的參數(shù)優(yōu)化實驗結果可知,此框架求出的NS效率系數(shù)有明顯的提高,并較大幅度減少了運行時間。并行框架具有良好的可擴展性,在大規(guī)模水文模擬中有較好的發(fā)揮。下一步的工作主要是考慮并行框架在異構眾核架構上的加速,同時探討其他的優(yōu)化算法與該并行框架的結合,使得在參數(shù)率定過程中進一步提高并行度。

        參考文獻

        [1]ABBOTT M B, BATHURST J C, CUNGE J A, et al. An introduction to the European Hydrological System-Systeme Hydrologique Europeen, “SHE”, 1: History and philosophy of a physically-based, distributed modelling system[J]. Journal of Hydrology, 1986, 87(1): 45-59.

        [2]LIU C M, WANG Z G, ZHENG H X, et al. Development of hydro-informatic modelling system and its application[J]. Science in China Series E: Technological Sciences, 2008, 51(4): 456-466

        [3]ALAM Z, RAHMAN M, ISLAM A. Assessment of climate change impact on the Meghna River Basin using geomorphology based hydrological model (GBHM)[C]// 3rd International Conference on Water and Flood Management (ICWFM-2011). Dhaka, 2011.

        [4]趙人俊, 王佩蘭. 新安江模型參數(shù)的分析[J]. 水文, 1988(6): 2-9.

        [5]雷曉輝, 廖衛(wèi)紅, 蔣云鐘, 等. 分布式水文模型EasyDHM模型[J]. 水利信息化, 2010, (3): 31-37.

        [6]WANG Q. The genetic algorithm and its application to calibrating conceptual rainfall-runoff models[J]. Water Resources Research, 1991, 27(9): 2467-2471.

        [7]CHENG C, OU C, CHAU K. 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]馬海波, 董增川, 張文明, 等. SCE-UA算法在TOPMODEL參數(shù)優(yōu)化中的應用[J]. 河海大學學報,2006, 34(4): 361-365.

        [9]THIEMANN M, TROSSET M, GUPTA H, et al. Bayesian recursive parameter estimation for hydrologic models[J]. Water Resources Research, 2001, 37(10): 2521-2535.

        [10] KOLLET S J, MAXWELL R M, WOODWARD C S, et al. Proof of concept of regional scale hydrologic simulations at hydrologic resolution utilizing massively parallel computer resources[J]. Water Resources Research, 2010, 46(4): W04201.1-W04201.7.

        [11] YALEW S G, VAN GRIENSVEN A, KOKOSZKIEWICZ L. Parallel computing of a large scale spatially distributed model using the Soil and Water Assessment Tool (SWAT) [C]// 2010 International Congress on Environmental Modelling and Software Modelling for Environment′s Sake, Fifth Biennial Meeting. Ottawa, 2010:1182-1189.

        [12] LI T J, WANG G Q, CHEN J, et al. Dynamic parallelization of hydrological model simulations[J]. Environmental Modelling & Software, 2011, 26(12): 1736-1746.

        [13] HUO J Y, LIU L Q, ZHANG Y N. An improved multi-cores parallel artificial Bee colony optimization algorithm for parameters calibration of hydrological model[J]. Future Generation Computer Systems, 2018, 81: 492-504.

        [14] 闞光遠, 洪陽, 梁珂, 等. 基于GPU加速的水文模型參數(shù)率定[J]. 人民長江, 2019, 50(5): 65-69,75.

        [15] 李強, 陸忠華, 王彥棡, 等. 基于克里金插值的SWAT參數(shù)率定大規(guī)模并行方法[J]. 計算機應用研究, 2016,33(1): 60-63.

        [16] 李強, 陸忠華, 王彥棡, 等. 基于高性能計算的SWAT參數(shù)敏感度分析并行框架[J]. 計算機應用研究, 2015, 32(1): 41-44+70

        [17] 申蒙蒙, 陸忠華, 王彥棡. 水文模擬中并行參數(shù)優(yōu)化算法[J]. 計算機工程與設計, 2017, 38(4): 1002-1007.

        [18] NASH J E, SUTCLIFFE J V. River flow forecasting through conceptual models part I—A discussion of principles[J]. Journal of Hydrology, 1970, 10(3): 282-290.

        [19] GROPP W, LUSK E, DOSS N, et al. A high-performance, portable implementation of the MPI message passing interface standard[J]. Parallel Computing, 1996, 22(6): 789-828

        [20] 王亞茹, 王鵬, 王德志. 基于MPI的多核并行模式的性能測試與分析[J]. 成都信息工程大學學報, 2018, 33(6): 617-623.

        [21] MAN K F, TANG K S, KWONG S. Genetic algorithms: concepts and designs[M]. Berlin: Springer Science & Business Media, 2012.

        [22] 王蓉, 劉遵仁, 高陽. 基于粒子群優(yōu)化和鄰域粗糙集的快速約簡算法[J]. 青島大學學報(自然科學版), 2017, 30(3):51-54.

        猜你喜歡
        優(yōu)化算法
        淺議小學數(shù)學口算教學的有效策略
        云計算平臺聯(lián)合資源調度優(yōu)化算法研究
        PLC故障檢測優(yōu)化算法
        原子干涉磁力儀信號鑒頻優(yōu)化算法設計
        故障樹計算機輔助分析優(yōu)化算法研究與應用
        混沌優(yōu)化算法在TSP問題的應用
        基于混沌初始化和高斯擾動的煙花算法
        計算機時代(2016年7期)2016-07-15 16:12:30
        再制造閉環(huán)供應鏈研究現(xiàn)狀分析
        二進制數(shù)轉十進制優(yōu)化算法探討
        故障樹計算機輔助分析優(yōu)化算法的實踐應用
        科技傳播(2016年3期)2016-03-25 00:23:31
        日韩一区二区三区无码影院| 国产亚洲精品综合99久久| 在线观看免费视频发布白白色| 亚洲视频一区二区三区视频| 中文字幕欧美人妻精品一区| 最新国产三级| 国产精品亚洲一区二区三区正片 | 亚洲色欲色欲www| 无码中文字幕加勒比一本二本 | 精品久久久中文字幕人妻| 精品国精品国产自在久国产应用| 无遮高潮国产免费观看韩国| 最新69国产精品视频| 中文字幕乱码高清完整版| 欧美俄罗斯乱妇| 无码三级国产三级在线电影| 国产的自拍av免费的在线观看 | 日韩一二三四精品免费| 大香蕉视频在线青青草| 国产精品视频一区二区三区不卡| 精品少妇ay一区二区三区| 无码AV无码免费一区二区| 大香蕉国产av一区二区三区 | 91精品国产福利在线观看麻豆| 色avav色av爱avav亚洲色拍| 成人欧美在线视频| 国产av精品久久一区二区| 成年av动漫网站18禁| 中国xxx农村性视频| 国产精品国三级国产av| 顶级高清嫩模一区二区| 网禁拗女稀缺资源在线观看| 亚洲人成综合网站在线| 国产日韩精品视频一区二区三区 | 亚洲欧洲高潮| 在线观看日韩精品视频网站| 亚洲综合av一区二区三区蜜桃| 国产99久久精品一区二区| 丰满熟妇人妻无码区| 少妇被爽到高潮喷水免费福利 | 国产亚洲一区二区三区综合片|