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

        ?

        感潮河段單站等時(shí)距水位非線性優(yōu)化預(yù)報(bào)模型

        2017-05-15 09:49:49王如云
        水道港口 2017年2期
        關(guān)鍵詞:分潮潮位潮汐

        王如云,李 丹,周 鈞,汪 天

        (1.河海大學(xué)海洋學(xué)院,南京210098;2.河海大學(xué)港口海岸與近海工程學(xué)院,南京210098;3.河海大學(xué)水文水資源學(xué)院,南京210098)

        感潮河段單站等時(shí)距水位非線性優(yōu)化預(yù)報(bào)模型

        王如云1,李 丹2,周 鈞3,汪 天2

        (1.河海大學(xué)海洋學(xué)院,南京210098;2.河海大學(xué)港口海岸與近海工程學(xué)院,南京210098;3.河海大學(xué)水文水資源學(xué)院,南京210098)

        目前感潮河段水位預(yù)報(bào)模型基于當(dāng)前站水位、上游徑流水位和下游海洋潮位數(shù)據(jù)進(jìn)行預(yù)報(bào),導(dǎo)致預(yù)見(jiàn)期受上游來(lái)流時(shí)間制約、上下游站數(shù)據(jù)缺失情況下無(wú)法進(jìn)行預(yù)報(bào)、數(shù)據(jù)采集成本高等問(wèn)題。文章假設(shè)感潮河段水位由徑流因素和潮汐因素共同線性作用,基于單站等時(shí)距水位數(shù)據(jù),建立了包含一個(gè)徑流動(dòng)力修正系數(shù)和若干個(gè)潮汐調(diào)和分潮系數(shù)的非線性優(yōu)化預(yù)報(bào)模型,并結(jié)合曲線擬合最小二乘法、潮汐調(diào)和分析法、梯度下降法,給出一個(gè)尋優(yōu)的迭代算法,對(duì)模型參數(shù)進(jìn)行優(yōu)化辨識(shí)。模型應(yīng)用于長(zhǎng)江感潮河段部分觀測(cè)站,水位預(yù)報(bào)精度較好,具有可行性和實(shí)用價(jià)值。

        等時(shí)距;感潮河段;潮汐調(diào)和分析;非線性優(yōu)化預(yù)報(bào)模型

        感潮河段是指河流下游潮區(qū)界至河口之間的河段。與上游河段不同,感潮河段受到上游徑流和下游潮汐共同作用,水位預(yù)報(bào)工作較為困難。

        目前感潮河段水位預(yù)報(bào)的方法有很多,根據(jù)預(yù)報(bào)前對(duì)水位數(shù)據(jù)處理方法的不同,可將其分為兩種類型,第一類是直接對(duì)水位進(jìn)行預(yù)報(bào),主要利用水動(dòng)力模型[1-2]、概率統(tǒng)計(jì)方法[3]及神經(jīng)網(wǎng)絡(luò)模型[4-5]等完成預(yù)報(bào)。第二類則是將感潮河段水位進(jìn)行潮、洪分離后分別對(duì)徑流水位和天文潮水位進(jìn)行預(yù)報(bào)校準(zhǔn)。

        針對(duì)感潮河段水位的復(fù)雜性,將感潮河段水位進(jìn)行潮、洪分離后處理不失為一個(gè)好方法。早在1989年,陳尚渭等人[6]就提出感潮河段水位是上游徑流和海洋潮波的復(fù)合。2003年,黃國(guó)如、芮孝芳[7]將徑流和天文潮水位分離研究,提出了建立感潮河段水位與上游徑流、下游天文潮因素之間的線性關(guān)系,采用頻率組合法預(yù)報(bào)感潮河段的設(shè)計(jì)洪水位。2013年,Smith等[8]將蘇格蘭迪河感潮河段水位進(jìn)行潮、洪分離,對(duì)迪河兩年的天文潮水位數(shù)據(jù)進(jìn)行潮汐調(diào)和分析,并與非線性預(yù)報(bào)的徑流水位進(jìn)行耦合,利用若干站點(diǎn)數(shù)據(jù)進(jìn)行逐次修正完成預(yù)報(bào)。

        此外,潮汐調(diào)和分析方法在感潮河段水位預(yù)報(bào)方面的應(yīng)用也十分廣泛。20世紀(jì)80年代,陳尚渭等[6]根據(jù)長(zhǎng)江上游水位及流量數(shù)據(jù)、下游天文潮特征,提出時(shí)間序列分析與潮汐調(diào)和分析相結(jié)合的方法進(jìn)行潮位預(yù)報(bào)。1985年,Speer等[9]指出了感潮河段中,天文潮與徑流的相互作用、河床的摩擦、河道不斷變化的幾何形態(tài)存在非線性關(guān)系,直接影響了感潮河段的水位。1999年,Godin[10]證明在Saint Lawrence河段可以用河流流量及潮汐分潮的振幅、潮差來(lái)預(yù)報(bào)上游站的水位,并對(duì)潮汐調(diào)和常數(shù)的變化進(jìn)行校準(zhǔn),進(jìn)一步指出徑流流量與潮汐調(diào)和常數(shù)之間的變化關(guān)系。2010年,Nidzieko[11]發(fā)現(xiàn)了感潮河段在漲潮和退潮過(guò)程中存在不對(duì)稱現(xiàn)象,直接導(dǎo)致了感潮河段高低潮(極值水位)出現(xiàn)的時(shí)間和潮位不規(guī)律。2013年,Li Guofang等[12]利用調(diào)和分析方法,構(gòu)建了長(zhǎng)江感潮河段潮位預(yù)報(bào)模型,同時(shí)通過(guò)對(duì)實(shí)測(cè)潮位的預(yù)報(bào),得到了預(yù)報(bào)誤差與徑流之間存在的關(guān)系,建立了以大通水文站實(shí)測(cè)流量為輸入數(shù)據(jù)的實(shí)時(shí)校正模型。同年,Moftakhari等[13]通過(guò)連續(xù)32 d的潮汐數(shù)據(jù)進(jìn)行調(diào)和分析,通過(guò)得到的參數(shù)對(duì)舊金山潮汐數(shù)據(jù)以及薩克拉門托河的水位流量數(shù)據(jù)進(jìn)行校正,較為準(zhǔn)確地反應(yīng)了該感潮河段的水位及流量狀況。

        目前的感潮河段站位的水位預(yù)報(bào)方法基本都基于對(duì)上游徑流水位觀測(cè)站、下游海洋潮位站,以及當(dāng)前感潮河段站位的水位觀測(cè)資料的分析,建立預(yù)報(bào)模型。然而該方法存在如下問(wèn)題:(1)預(yù)見(jiàn)期受上游來(lái)流時(shí)間制約;(2)如果出現(xiàn)上下游站位數(shù)據(jù)缺失,預(yù)報(bào)將無(wú)法進(jìn)行;(3)為了對(duì)感潮河段的水位進(jìn)行預(yù)報(bào),除了建立當(dāng)前觀測(cè)站外,還必須在其上游和下游分別建立徑流水位和海洋潮位觀測(cè)站,勢(shì)必增加很大費(fèi)用。為了克服這些問(wèn)題,建立了只需要當(dāng)前站水位的數(shù)據(jù)進(jìn)行水位預(yù)報(bào)的非線性數(shù)學(xué)模型,其基本思想是,在假設(shè)感潮河段水位由徑流因素和潮汐因素共同線性作用引起的情況下,通過(guò)考慮對(duì)徑流采用一個(gè)帶動(dòng)力修正項(xiàng)的水位預(yù)報(bào)模型,對(duì)潮汐采用調(diào)和分潮和函數(shù)的潮位預(yù)報(bào)模型情況下,基于曲線擬合最小二乘法,建立了關(guān)于一個(gè)動(dòng)力修正系數(shù)和若干個(gè)分潮調(diào)和系數(shù)的非線性目標(biāo)函數(shù),并結(jié)合調(diào)和分析法、梯度下降法,給出了一個(gè)尋優(yōu)的迭代算法,對(duì)模型參數(shù)進(jìn)行優(yōu)化辨識(shí)。經(jīng)過(guò)對(duì)長(zhǎng)江感潮河段一些觀測(cè)站水位進(jìn)行預(yù)報(bào)檢驗(yàn),結(jié)果表明預(yù)報(bào)精度較好,具有可行性和實(shí)用價(jià)值。

        1 感潮河段單站等時(shí)距水位非線性優(yōu)化預(yù)報(bào)模型

        假設(shè)感潮河段某水位觀測(cè)站等間距時(shí)刻ti的水位為ηi

        式中:η(t)為感潮河段水位。

        由于感潮河段水位受上游徑流和下游天文潮共同作用,所以假設(shè)感潮河段水位可分離為徑流水位和天文潮水位之和,即有

        這里ηˉ(t)為徑流水位,η~(t)為天文潮水位。

        假設(shè)徑流水位可以用3個(gè)連續(xù)等距時(shí)刻的水位,采用一個(gè)帶動(dòng)力修正項(xiàng)的水位預(yù)報(bào)模型進(jìn)行預(yù)報(bào),形式如下

        這里α[ηˉ(ti-l-h)-ηˉ(ti-l-2h)]稱為動(dòng)力修正項(xiàng),α稱為動(dòng)力修正系數(shù),l為預(yù)見(jiàn)期,h為3個(gè)連續(xù)時(shí)刻的時(shí)間間距。

        假設(shè)天文潮水位可以用含如下若干個(gè)調(diào)和分潮和函數(shù)的形式表示

        這里A0、Ak、Bk為調(diào)和分潮系數(shù),m為分潮總數(shù),ωk為第k個(gè)分潮的角速度。

        分別利用式(3)和(4)對(duì)徑流水位ηˉ(ti)和天文潮水位η~(t)進(jìn)行預(yù)報(bào),為求出動(dòng)力修正系數(shù)和調(diào)和分潮系數(shù)的值,根據(jù)最小二乘法建立水位均方差目標(biāo)函數(shù)如下

        這里預(yù)見(jiàn)期l和時(shí)間間距h為整數(shù)。

        據(jù)式(2)有

        這里

        據(jù)式(4)有

        這里

        故有

        由于Gc(k,i)和Gs(k,i)為 α的函數(shù),為使目標(biāo)函數(shù)達(dá)到最小值,由 φ(α,A1,…,Am,B1,…,Bm)分別對(duì)α,A1,…,Am,B1,…,Bm的導(dǎo)數(shù)為零,形成了非線性的代數(shù)方程組。利用調(diào)和分析法和梯度下降法給出了辨識(shí)模型參數(shù)α,A1,…,Am,B1,…,Bm的迭代算法如下:

        (1)給定辨識(shí)精度指標(biāo)ε>0,α的初始值α0。

        (2)當(dāng) α=α0時(shí),為得到 φ(α0,A1,…,Am,B1,…,Bm)的最小值,需滿足 φ(α0,A1,…,Am,B1,…,Bm)對(duì) Al、 Bl(l=1,2,…,m)的導(dǎo)數(shù)等于0。又因?yàn)?0,其中x表示A1,…,Am,B1,…,Bm,則有

        經(jīng)化簡(jiǎn)可得2m元一次方程組

        解方程組(13)可以得到A1,…,Am,B1,…,Bm的值。

        (3)給定學(xué)習(xí)速率初值λ=λ0,學(xué)習(xí)次數(shù)限定值Ktotal,學(xué)習(xí)次數(shù)計(jì)數(shù)器knum=1。

        (4)將α作為未知量,對(duì)式(11)求導(dǎo),可以得到

        表1 基本調(diào)和常數(shù)Tab.1 Basic harmonic constant

        表2 南京2011~2012年水位預(yù)報(bào)模型后報(bào)誤差Tab.2 Hindcast error of water level of Nanjing from 2011 to 2012

        表3 南京2011~2012年水位后報(bào)參數(shù)預(yù)報(bào)2013~2015年水位誤差Tab.3 Forecast error of the water level from 2013 to 2015 based on the hindcast parameters of the water level of Nanjing from 2011 to 2012

        表4 南京2013~2015年預(yù)報(bào)誤差分布(預(yù)見(jiàn)期6 h)Tab.4 Distribution of forecast error of Nanjing from 2013 to 2015(forecast period 6 h)

        表5 南京2013~2015年預(yù)報(bào)誤差分布(預(yù)見(jiàn)期24 h)Tab.5 Distribution of forecast error of Nanjing from 2013 to 2015(forecast period 24 h)

        利用迭代得到的模型參數(shù)α,A1,…,Am,B1,…,Bm,可以完成當(dāng)前站位的水位預(yù)報(bào)工作。對(duì)于任意時(shí)刻t,由式(2)和(3)得

        這 里η(t)、η(t-h)和η(t-2h)為 已 知 數(shù) 據(jù) ,將α,A1,…,Am,B1,…,Bm帶入式(15),結(jié)合式(4),可以得到當(dāng)前站位t+l時(shí)刻的預(yù)報(bào)水位η(t+l)。

        圖1 南京2013年預(yù)報(bào)水位誤差與實(shí)際水位對(duì)比分析(局部)Fig.1 Contrast between the error of forecast water level and the actual water level of Nanjing in 2013

        2 水位預(yù)報(bào)模型的檢驗(yàn)

        在模型的計(jì)算檢驗(yàn)過(guò)程中,選擇潮汐的11個(gè)分潮進(jìn)行調(diào)和分析計(jì)算,包括4個(gè)全日分潮、4個(gè)半日分潮和3個(gè)主要的淺水分潮,其分潮角速度見(jiàn)表1。

        將感潮河段單站等時(shí)距水位優(yōu)化預(yù)報(bào)模型應(yīng)用于長(zhǎng)江下游感潮河段。限于已有數(shù)據(jù),利用南京2011~2012年全年等時(shí)距水位數(shù)據(jù),對(duì)模型參數(shù)α和調(diào)和常數(shù)進(jìn)行辨識(shí),并作后報(bào)檢驗(yàn)(表2)。再對(duì)2013~2015年的南京站水位進(jìn)行預(yù)報(bào),預(yù)報(bào)的均方差見(jiàn)表3。

        后報(bào)及預(yù)報(bào)的預(yù)見(jiàn)期可根據(jù)需要進(jìn)行調(diào)整,這里選擇6~72 h間共12個(gè)時(shí)間長(zhǎng)度作為預(yù)見(jiàn)期。由表2、表3數(shù)據(jù)可以看出,整體水位后報(bào)及預(yù)報(bào)精確度良好,尤其當(dāng)預(yù)見(jiàn)期為12、24、36、48、60、72 h,誤差較小。同時(shí),針對(duì)預(yù)報(bào)的誤差情況進(jìn)行了分析,表4、表5主要反映了南京站預(yù)報(bào)水位在預(yù)見(jiàn)期為6 h和24 h的誤差分布,可以看到當(dāng)預(yù)見(jiàn)期為6 h、限定誤差精度為0.3 m時(shí),合格率在97%~98%。為保證合格率在85%以上,預(yù)見(jiàn)期可以選取36 h以內(nèi)時(shí)長(zhǎng)。而圖1可以看到,水位預(yù)報(bào)誤差整體較小,部分誤差較大的預(yù)報(bào)時(shí)刻主要集中在徑流變化趨勢(shì)較大的時(shí)間段。

        模型的驗(yàn)證可以看到,整體預(yù)報(bào)誤差和傳統(tǒng)的結(jié)合上下游數(shù)據(jù)的水位預(yù)報(bào)模型結(jié)果[14]相比不差,部分地方甚至較好,且該模型可以針對(duì)不同預(yù)見(jiàn)期進(jìn)行預(yù)報(bào)。由此可見(jiàn),該水位優(yōu)化預(yù)報(bào)模型具有一定的可行性和實(shí)用性。

        3 總結(jié)與展望

        限于目前感潮河段的水位預(yù)報(bào)基本都需要上游徑流水文觀測(cè)站和下游海洋潮位觀測(cè)站的數(shù)據(jù),若出現(xiàn)上下游站位數(shù)據(jù)缺失等突發(fā)情況時(shí),感潮河段單站等時(shí)距水位優(yōu)化預(yù)報(bào)模型可以作為很好的替代模型使用,其獨(dú)立性保證了模型不會(huì)受到上下游數(shù)據(jù)缺失的影響而進(jìn)行較為準(zhǔn)確的水位預(yù)報(bào)。同時(shí),與現(xiàn)有模型相比,該水位預(yù)報(bào)模型更好地挖掘了現(xiàn)有的水文數(shù)據(jù)信息,數(shù)據(jù)利用率較高。且該模型整體預(yù)報(bào)精度較好,預(yù)見(jiàn)期可根據(jù)需要靈活調(diào)整,在保證精度的前提下,最長(zhǎng)可以預(yù)報(bào)36 h的水位,適用性更強(qiáng)。因此,在只關(guān)注感潮河段水位變化的情況下,便不需要花費(fèi)更多的人力和物力設(shè)立上下游的水文站。后期,可以將該模型改造成實(shí)時(shí)校準(zhǔn)并連續(xù)預(yù)報(bào)的模式,加入洪水期和枯水期因素,增加其預(yù)報(bào)精度。還可以將模型與現(xiàn)有感潮河段水位模型進(jìn)行結(jié)合,進(jìn)一步改進(jìn)現(xiàn)有水位預(yù)報(bào)方法,提高水文數(shù)據(jù)的利用率,改善預(yù)報(bào)結(jié)果。

        [1]劉曉波,彭文啟,劉靜玲,等.基于POM模式的感潮河段三維水動(dòng)力模擬[J].水力發(fā)電學(xué)報(bào),2009,28(1):89-94. LIU X B,PENG W Q,LIU J L,et al.3?D numerical simulation of tidal flow based on POM model[J].Journal of Hydroelectric Engi?neering,2009,28(1):89-94.

        [2]張小琴,包為民,馬德蓮.雙向波水位演算模型參數(shù)動(dòng)態(tài)修正[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,2010,21(5):71-76. ZHANG X Q,BAO W M,MA D L.Dynamic correction of parameters for bi?directional stage routing model[J].Journal of Hohai University:Natural Sciences,2010,21(5):71-76.

        [3]林勛勵(lì).長(zhǎng)江口天文潮預(yù)報(bào)修正的一種方法[J].海洋預(yù)報(bào),1985(1):34-42. LIN X L.An adjusting method of astronomical tide forecast in the Changjiang estuary[J].Marine Forecast Service,1985(1):34-42.

        [4]Supharatid S.Tidal?level forecasting and filtering by neural network model[J].Coastal Engineering Journal,2003,45(1):119-137.

        [5]Hidayat H,Hoitink A J F,Sassi M G,et al.Prediction of Discharge in a Tidal River Using Artificial Neural Networks[J].JOUR?NAL OF HYDROLOGIC ENGINEERING,2014,19(8):DOI:10.1061/(ASCE)HE.1943-5584.0000970.

        [6]陳尚渭,金兆森.長(zhǎng)江潮位預(yù)報(bào)方法的研究[J].水利學(xué)報(bào),1989,6:41-47. CHEN S W,JIN Z S.Study of Method of Predicting Tidal Level at Lower Reach of the Yangtze River[J].Journal of Hydraulic Engi?neering,1989,6:41-47.

        [7]黃國(guó)如,芮孝芳.感潮河段設(shè)計(jì)洪水位計(jì)算的頻率組合法[J].水電能源科學(xué),2003,21(2):72-74. HUANG G R,RUI X F.Design Flood Water Level for Tidal Reach with Frequency Combination Method[J].Water Resources and Power,2003,21(2):72-74.

        [8]Smith P J,Beven K J,Horsburgh K.Data?based mechanistic modelling of tidally affected river reaches for flood warning purposes-An example on the River Dee,UK[J].Quarterly Journal of the Royal Meteorological Society,2013,671(139):340-349.

        [9]Speer P E,Aubrey D G.A study of non?linear tidal propagation in shallow inlet/estuarine systems Part II:Theory[J].Estuarine,Coastal and Shelf Science,1985,21(2):207-224.

        [10]Godin G.The propagation of tides up rivers with special considerations on the upper Saint Lawrence River[J].Estuarine,Coastal and Shelf Science,1999,48(3):307-324.

        [11]Nidzieko N J.Tidal asymmetry in estuaries with mixed semidiurnal/diurnal tides[J].Journal of Geophysical Research:Oceans(1978-2012),2010,115(C8).

        [12]Li G F,Xiang X Y,Wu J,et al.Long?Term Water?Level Forecasting and Real?Time Correction Models in the Tidal Reach of the Yangtze River[J].JOURNAL OF HYDROLOGIC ENGINEERING,2013,18(11):1 437-1 442.

        [13]Moftakhari H R,Jay D A,Talke S A,et al.A novel approach to flow estimation in tidal rivers[J].Water Rrsources Research,2013,49:4 817-4 832.

        [14]方新.長(zhǎng)江感潮河段水位過(guò)程預(yù)報(bào)模型研究[D].南京:南京師范大學(xué),2014.

        Nonlinear optimal forecast model of water level at time knots with constant time interval from a single station of tidal reach

        WANG Ru?yun1,LI Dan2,ZHOU Jun3,WANG Tian2
        (1.College of Oceanography,Hohai University,Nanjing 210098,China;2.College of Harbor,Coastal and Offshore Engineering,Hohai University,Nanjing 210098,China;3.College of Hydrology and Water Resources,Hohai University,Nanjing 210098,China)

        The current water level forecast models of tidal reach are subject to the observation data of water lev?el from the current tidal reach station,the runoff stations upstream and the tide stations downstream.It leads to the problems as follows:the leading time is subject to the time of the flow from the upstream;if the data at the upstream or the downstream are missed,the forecast would not be able to proceed;the acquisition of the data is high in cost. Assuming that the water level of tidal reach was linearly controlled by the runoff and tide,the paper built a nonlin?ear optimal forecast model with respect to one dynamic correction coefficient of runoff and a set of tidal harmonic co?efficients,which was based on the water level at time knots with constant time interval from a single station.Then, combining the curve fitting least square method,the harmonic analysis method and the gradient descent method,an iterative algorithm was proposed for solving the forecast model.The model was applied to some stations of the tidal reach of Yangtze River for test.Numerical experiences show the feasible and applicable of the model.

        time knots with constant time interval;tidal reach;tidal harmonic analysis;nonlinear optimal fore?cast model

        TV 133;P 731.23

        A

        1005-8443(2017)02-0132-05

        2016-12-01;

        2017-03-28

        中國(guó)江蘇省水利科技重點(diǎn)項(xiàng)目(2010500312);中央高?;究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(2014B06314)

        王如云(1963-),男,安徽蕪湖人,教授,主要從事物理海洋學(xué)和港口海岸與近海工程研究。

        Biography:WANG Ru?yun(1963-),male,professor.

        猜你喜歡
        分潮潮位潮汐
        潮汐與戰(zhàn)爭(zhēng)(上)
        基于距離倒數(shù)加權(quán)的多站潮位改正方法可行性分析
        大亞灣雙峰水位的形成條件及準(zhǔn)調(diào)和分量應(yīng)用的分析
        唐山市警戒潮位標(biāo)志物維護(hù)研究
        山東鄰海長(zhǎng)周期分潮對(duì)深度基準(zhǔn)面的影響分析
        絕美海灘
        多潮位站海道地形測(cè)量潮位控制方法研究
        基于改進(jìn)的OLS-RBF模型的感潮河段潮位預(yù)測(cè)研究
        潮汐式灌溉控制系統(tǒng)的設(shè)計(jì)及應(yīng)用
        電子制作(2017年9期)2017-04-17 03:00:56
        干法紙的潮汐
        生活用紙(2016年6期)2017-01-19 07:36:25
        国产成人亚洲合色婷婷| 久久天天躁狠狠躁夜夜躁2014| 国产综合久久久久久鬼色| 毛片内射久久久一区| 欧美黑人疯狂性受xxxxx喷水| 久久亚洲Av无码专区| 国内精品久久久久久久亚洲| 国产精品一区二区久久精品蜜臀| 好看的中文字幕中文在线| 人妖一区二区三区视频| 人妻少妇被粗大爽.9797pw| 午夜福利av无码一区二区| 中国丰满大乳乳液| 国产精品国产午夜免费福利看| 日韩精品一级在线视频| 亚洲av日韩专区在线观看| 人妻丰满熟妇aⅴ无码| 少妇仑乱a毛片| 日韩久久一级毛片| 亚洲春色视频在线观看| 国产片在线一区二区三区| 国产乱子伦精品无码专区| 免费人成网ww555kkk在线| 97人人超碰国产精品最新o| 免费中文熟妇在线影片| 久久精品国产亚洲AⅤ无码剧情| 亚洲av中文字字幕乱码软件 | 久久老熟女一区二区三区福利| 午夜福利理论片在线观看播放 | 无码av一区二区大桥久未| 少妇被躁爽到高潮无码文| 久久久精品国产亚洲AV蜜| 国产精品三级av一区二区| av免费一区二区久久| 亚洲国产精品成人久久| 国产真实老熟女无套内射| 一区一级三级在线观看| 伊人影院在线观看不卡| 国产黄色一区二区在线看| 亚洲av高清在线观看一区二区 | 日本护士一区二区三区高清热线|