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

基于EEMD-AR模型的丹江口水庫年徑流隨機模擬與預(yù)報

2017-09-16 06:53:12練繼建孫蕭仲唐志波
水利水電科技進展 2017年5期
關(guān)鍵詞:模型

練繼建,孫蕭仲,馬 超,趙 明,唐志波

(天津大學(xué)水利工程仿真與安全國家重點實驗室,天津 300072)

基于EEMD-AR模型的丹江口水庫年徑流隨機模擬與預(yù)報

練繼建,孫蕭仲,馬 超,趙 明,唐志波

(天津大學(xué)水利工程仿真與安全國家重點實驗室,天津 300072)

基于水庫歷史年入庫徑流序列組分分析和識別,采用線性趨勢回歸檢驗法、有序聚類法、方差線譜法等方法,推求出序列趨勢項、跳躍項及周期項等確定性成分,提出基于集合經(jīng)驗?zāi)B(tài)分解法(EEMD方法)的水庫年徑流自回歸隨機模擬模型(EEMD-AR),并應(yīng)用于丹江口水庫的年徑流隨機模擬和預(yù)報中。通過EEMD分解,解決了當(dāng)?shù)そ谒畮鞖v史年徑流序列為非平穩(wěn)序列時不能直接應(yīng)用自回歸模型(AR) 進行隨機模擬和預(yù)報的問題。模擬結(jié)果表明,EEMD-AR模型能較好地模擬丹江口水庫年徑流序列并保持原歷史序列的統(tǒng)計特性,且模型預(yù)報精度符合要求。

徑流序列成分識別;EEMD-AR模型;徑流隨機模擬;丹江口水庫

隨著計算機技術(shù)的發(fā)展,水文隨機模擬技術(shù)除最初的線性平穩(wěn)隨機模型之外,人工神經(jīng)網(wǎng)絡(luò)[1]、小波分析[2]、混沌理論[3]、Copula方法[4]亦被引入到隨機模型中,并取得了較好的模擬效果。相較近來廣泛應(yīng)用于徑流模擬的分布式水文模型[5-6],傳統(tǒng)的線性平穩(wěn)隨機模型由于構(gòu)造簡單,在水文模擬中備受青睞,但一般僅適用于平穩(wěn)序列的隨機模擬,具有一定的局限性。集合經(jīng)驗?zāi)B(tài)分解法(EEMD方法)因其對非平穩(wěn)序列具有較強的處理能力,被不少學(xué)者引入到徑流隨機模擬和預(yù)報領(lǐng)域[7-10]。章國勇等[7]通過建立基于EEMD方法的LSSVM組合預(yù)測模型和基于動態(tài)逼近搜索粒子群的LSSVM參數(shù)尋優(yōu)方法,進行了江埡站的年徑流預(yù)測和模型驗證。Wang等[9]將EEMD方法結(jié)合人工神經(jīng)網(wǎng)絡(luò)模型進行年徑流序列的隨機模擬和預(yù)報研究,取得了較為理想的預(yù)報精度。

目前針對丹江口水庫入庫徑流方面的研究主要集中于徑流特征分析[11-12]和以月為時間尺度[13-15]的預(yù)報模型研究。然而,隨著丹江口水庫大壩加高和南水北調(diào)工程通水,針對丹江口水庫月尺度入庫徑流模擬和預(yù)報方面的已有研究不能完全滿足水庫多年調(diào)節(jié)的徑流輸入要求。此外,由于丹江口水庫歷史年徑流序列為非平穩(wěn)序列,不能直接采用傳統(tǒng)自回歸模型進行隨機模擬。基于此,本文考慮水庫年徑流年際間變化規(guī)律,在識別丹江口水庫年入庫徑流序列成分的基礎(chǔ)上,建立基于EEMD-AR的隨機模擬和預(yù)報模型,開展丹江口水庫非平穩(wěn)年徑流序列隨機模擬和預(yù)報研究,以期為丹江口水庫多年優(yōu)化調(diào)度提供可靠徑流數(shù)據(jù)。

1 研究思路和方法

1.1 研究思路

首先,對水庫歷史年徑流序列進行成分分析和識別。其中,非周期成分中的趨勢項采用線性趨勢回歸檢驗法、Spearman秩次相關(guān)檢驗法和Mann-Kendall檢驗法;跳躍項采用有序聚類法、秩和檢驗法和 Mann-Kendall 檢驗法;周期成分采用方差線譜法和累積解釋方差圖法進行檢驗。其次,對剔除趨勢、跳躍及周期等確定性成分后的隨機成分進行EEMD分解,并對分解后的各階固有模態(tài)函數(shù)IMF(intrinsic mode functions)序列進行AR模型建模,殘差序列進行多項式模擬,形成基于EEMD-AR的丹江口水庫年徑流隨機模擬模型和預(yù)報模型,研究技術(shù)路線見圖1。

圖1 研究技術(shù)路線

1.2 研究方法

EEMD方法是在傳統(tǒng)經(jīng)驗?zāi)B(tài)分解法(EMD方法)的基礎(chǔ)上,通過添加有限的且非無窮小振幅的白噪聲序列,從而改善傳統(tǒng)經(jīng)驗?zāi)B(tài)分解法中的模態(tài)混疊[16],其分解步驟如下:

(1)

步驟2 將新合成序列yi(t)進行EMD分解并得到對應(yīng)的各IMF分量cij和殘差ri:

(2)

式中:cij為第i次EMD分解后的第j階IMF分量;n為第i次EMD分解對應(yīng)的IMF分量數(shù);ri為第i次EMD分解得到的殘差。

步驟3 重復(fù)前2個步驟N次(N為總體平均次數(shù)),得到N組新的合成序列及相應(yīng)的各IMF分量和殘差。

步驟4 將N組新合成序列分解后的各IMF分量取總體平均,得到最終的各IMF分量:

(3)

2 基于EEMD-AR的丹江口水庫年徑流隨機模擬和預(yù)報模型

將1956—2009年丹江口水庫年徑流序列(圖2)分成1956—2003年和2004—2009年兩組,采用1956—2003年徑流序列建立基于EEMD-AR的徑流隨機模擬和預(yù)報模型,采用2004—2009年徑流序列進行模型預(yù)報精度檢驗。

作品改編權(quán)保護的歷史之維 ..................................李 楊 06.32

圖2 丹江口水庫年徑流序列(1956—2009年)

2.1 丹江口水庫年徑流序列成分分析

2.1.1 非周期成分識別

a. 趨勢項診斷。采用Mann-Kendall檢驗法、Spearman秩次相關(guān)檢驗法及線性趨勢回歸檢驗法對丹江口水庫年徑流序列(1956—2003)進行趨勢項檢驗,結(jié)果表明:Mann-Kendall檢驗法下丹江口水庫年徑流序列標準統(tǒng)計變量Z<0(Z=-1.73)且其絕對值小于臨界值Z1-a/2=1.96。因此,丹江口水庫年徑流序列有一定的減小趨勢,但是該趨勢不顯著。此外,Spearman秩次相關(guān)檢驗法計算得統(tǒng)計量T=1.75,線性趨勢回歸檢驗法計算得到統(tǒng)計量d=-0.42,其絕對值均小于自由度為n-2的t分布在5%顯著水平下的臨界值t1-a/2,n-2=2.01,表明Spearman秩次相關(guān)檢驗法及線性趨勢回歸檢驗法的檢驗結(jié)果與Mann-Kendall檢驗法一致(標準統(tǒng)計變量Z、統(tǒng)計量T和統(tǒng)計量d的計算可參考王文圣等的《隨機水文學(xué)》[18])。

b. 變異項診斷。采用有序聚類法、滑動秩和檢驗法及Mann-Kendall檢驗法對丹江口水庫年徑流序列(1956—2003)進行跳躍診斷,各檢驗結(jié)果見圖3~4。其中,有序聚類法選取檢驗結(jié)果中統(tǒng)計量S值的最小值點作為最可能變異點;滑動秩和檢驗法選取統(tǒng)計量U值超過臨界值最多的點作為最可能變異點;Mann-Kendall檢驗法則選取統(tǒng)計曲線UB和UF在臨界區(qū)間內(nèi)的交點作為最可能變異點。由圖3~4可知,3種不同檢驗方法識別的變異結(jié)果一致,丹江口水庫年徑流序列的最可能變異點為1990年。

圖3 有序聚類法和滑動秩和檢驗法檢驗結(jié)果

圖4 Mann-Kendall檢驗法檢驗結(jié)果

根據(jù)變異項診斷結(jié)果,對存在跳躍變異的序列進行預(yù)處理,即將變異點之后的序列加上變異前后兩序列均值之差,從而剔除確定性成分中的跳躍項部分,使序列具有一致性[18]。

2.1.2 周期成分識別

對經(jīng)過預(yù)處理后的丹江口水庫年徑流序列(1956—2003)采用方差線譜法和累積解釋方差圖法進行周期分析,見圖5和圖6。圖5中占明顯優(yōu)勢的振幅Aj所對應(yīng)頻率wj的倒數(shù)即為周期,圖6中轉(zhuǎn)折點對應(yīng)的諧波個數(shù)即為周期數(shù)(橫坐標標目m為諧波個數(shù),縱坐標標目Bm為累積解釋方差)。確定序列周期后,可根據(jù)式(4)確定序列的傅里葉系數(shù)(a0、aj和bj),如表1所示,并將其代入式(5)中以確定序列中的周期成分。

(4)

(5)

式中:a0、aj和bj均為序列傅里葉系數(shù);wj為頻率;l為諧波總個數(shù);n′為序列長度。

圖5 方差線譜

圖6 累積解釋方差

表1 剔除跳躍項后的序列周期計算結(jié)果及相應(yīng)傅里葉系數(shù)

2.2 丹江口水庫年徑流序列EEMD分解

根據(jù)確定性成分識別成果,將原始丹江口水庫年徑流序列先進行跳躍項剔除預(yù)處理,即將跳躍變異點后序列值加上前后序列均值之差;在此基礎(chǔ)上減去由式(5)得出的周期成分,從而剔除周期項,最終可得序列隨機成分(趨勢項不顯著,不作處理)。對剔除確定性成分后的序列隨機項進行EEMD分解和ADF平穩(wěn)性檢驗,結(jié)果如下:

a. EEMD分解結(jié)果。采用EEMD法對處理后的丹江口水庫年徑流序列(1956—2003)隨機項進行分解,總體平均次數(shù)N取為1 000次,噪聲標準差為原始序列標準差的20%。經(jīng)EEMD分解后的各階固有模態(tài)函數(shù)IMF1~IMF4序列及殘差如圖7所示。

圖7 丹江口水庫年徑流序列隨機項EEMD分解結(jié)果

b. ADF平穩(wěn)性檢驗。將IMF1~IMF4序列采用Eviews軟件進行ADF單位根法檢驗,檢驗結(jié)果顯示經(jīng)成分剔除后的丹江口水庫年徑流序列的t值(-0.42)大于顯著性水平為5%時對應(yīng)的臨界值(-1.95),且p值(0.53)大于0.05,因此其為非平穩(wěn)序列,不能直接采用AR模型進行模擬。而經(jīng)EEMD分解后的IMF1~IMF4序列的t值(分別為-7.06、-2.54、-3.26和-3.41)均小于顯著性水平為5%時對應(yīng)的臨界值,且p值均小于0.05,即IMF1~IMF4均為平穩(wěn)序列,可采用AR模型進行隨機模擬。

2.3 基于EEMD-AR的丹江口水庫年徑流隨機模擬

2.3.1 殘差序列的多項式模擬

采用多項式對EEMD分解后的殘差序列進行模擬,模擬結(jié)果如下:

Rt=6×10-5t3-0.012 7t2+1.721 2t+358.87

(6)

式中:Rt為殘差項模擬值;t為年份序號。

2.3.2 IMF1~IMF4序列AR模型建立

在建立AR模型前,對IMF1~IMF4序列進行獨立性檢驗,檢驗結(jié)果顯示分解后的IMF序列均為相依序列,可采用AR模型進行隨機模擬。

a. AR模型階數(shù)確定。采用AIC準則對IMF1~IMF4序列進行模型階數(shù)選擇,并根據(jù)AIC值最小原則及相應(yīng)的偏態(tài)系數(shù)Cs,最終選定IMF1序列采用正態(tài)AR(3)模型,IMF2序列采用偏態(tài)AR(4)模型,IMF3序列采用偏態(tài)AR(3)模型,IMF4序列采用偏態(tài)AR(4)模型進行隨機模擬。

b. IMF1~IMF4序列AR模型建立。建立IMF1~IMF4序列對應(yīng)的AR模型,其中AR(p)模型的數(shù)學(xué)表達式見式(7)[14],各IMF序列AR模型參數(shù)見表2,表中σ為序列標準差。

(7)

式中:μ為序列均值;φ1,φ2,…,φp為自回歸系數(shù);εt為純隨機序列。

表2 IMF1~IMF4序列AR模型參數(shù)

c. 最終模擬序列生成。將各階IMF1~IMF4序列的AR模型采用Monte Carlo法進行隨機模擬,并將模擬結(jié)果與殘差項模擬結(jié)果求和,并在此基礎(chǔ)上還原原序列中的確定性成分,即可得到大量丹江口水庫年徑流隨機模擬序列,如式(8)和式(9)所示:

(8)

(9)

式中:x為最終模擬序列;xqd為確定性成分;yIMFi為第i階IMF模擬序列;R為殘差模擬序列。

2.3.3 模擬結(jié)果檢驗

表3 模擬結(jié)果檢驗

2.4 基于EEMD-AR的丹江口水庫年徑流預(yù)報

基于EEMD-AR的丹江口水庫年徑流預(yù)報模型與隨機模擬模型表達形式大致相同,預(yù)報值為確定性成分(跳躍成分、周期成分)、IMF1~IMF4序列和殘差R預(yù)報值的疊加之和,唯一區(qū)別是IMF1~IMF4序列的AR預(yù)報模型不包含隨機項εt。采取預(yù)留的2004—2009年丹江口水庫實測年徑流來檢驗預(yù)報模型。丹江口水庫2004—2009年的年徑流預(yù)報結(jié)果如表4所示,預(yù)報精度均滿足GB/T 22482—2008《水文情報預(yù)報規(guī)范》的精度要求,平均誤差為9.54%,最大誤差不超過18%。預(yù)報結(jié)果表明基于EEMD-AR的丹江口水庫年徑流預(yù)報模型可較好地預(yù)報未來年徑流情況。

表4 年徑流量預(yù)報結(jié)果檢驗

3 結(jié) 論

本文基于歷史年徑流序列進行序列成分識別,開展基于EEMD-AR的丹江口水庫年徑流隨機模擬和預(yù)報。丹江口水庫歷史年徑流序列(1956—2003)成分識別結(jié)果顯示:①丹江口水庫年徑流有一定的減小趨勢,但不顯著;②水庫年徑流序列跳躍顯著且最可能變異點為1990年;③水庫年徑流序列存在3.2 a和9.6 a的周期變化。

剔除趨勢、跳躍及周期等確定性成分后的隨機序列通過EEMD分解,較好地解決了非平穩(wěn)序列不能直接進行AR模型隨機模擬的問題。通過對EEMD分解后的IMF1~IMF4序列進行AR模型建模以及殘差序列進行多項式擬合,建立基于EEMD-AR的丹江口水庫年徑流隨機模擬模型和預(yù)報模型。結(jié)果表明,模擬序列各統(tǒng)計指標相對誤差絕對值最大值不超過9%,能較好地保持原序列的統(tǒng)計特性;預(yù)報序列平均誤差9.54%,最大誤差不超過18%。因此,所提出的模型能較好地適應(yīng)丹江口水庫年入庫徑流序列的模擬和預(yù)報。

[ 1 ] 閻俊愛,鐘登華.基于遺傳算法的神經(jīng)網(wǎng)絡(luò)優(yōu)化預(yù)測模型及其在年徑流預(yù)報中的應(yīng)用[J].水利水電技術(shù),2003,34(6):1-4.(YAN Junai,ZHONG Denghua.Forecasting method of neural network optimization based on genetic algorithm and its application in annual runoff forecast[J].Water Resources and Hydropower Engineering,2003,34(6): 1-4.(in Chinese))

[ 2 ] 杜富慧.基于小波分析-模糊神經(jīng)網(wǎng)絡(luò)的徑流預(yù)報模型[J].水利水電技術(shù),2013,44(2): 5-8.(DU Fuli.Wavelet analysis-fuzzy neural network based runoff forecasting mode[J].Water Resources and Hydropower Engineering,2013,44(2):5-8.(in Chinese))

[ 3 ] 計亞麗.基于混沌理論的呼倫湖流域降雨徑流時間序列預(yù)測分析[D].呼和浩特:內(nèi)蒙古農(nóng)業(yè)大學(xué),2012.

[ 4 ] 周研來,梅亞東,張代青,等.一種新的徑流過程隨機模擬方法[J].水利水電科技進展,2011,31(3):9-12.(ZHOU Yanlai,MEI Yadong,ZHANG Daiqing,et al.A aew stochastic method for streamflow hydrographs[J].Advances in Science and Technology of Water Resources,2011,31(3):9-12.(in Chinese))

[ 5 ] 劉柏君,權(quán)錦,楊超慧,等.基于SWAT模型的青海省主要流域徑流模擬[J].水資源保護,2016,32(6):39-44.(LIU Bojun,QUAN Jin,YANG Chaohui,et al.Simulation of runoff of major basins in Qinghai Province based on SWAT model[J].Water Resources Protection,2016,32(6):39-44.(in Chinese))

[ 6 ] 李致家,屈晨陽,黃鵬年,等.CASC2D模型與GSSHA模型在欒川流域的徑流模擬[J].河海大學(xué)學(xué)報(自然科學(xué)版),2017,45(1):1-6.(LI Zhijia,QU Chenyang,HUANG Pengnian,et al.Simulation of runoff in Luanchuan Basin using CASC2D and GSSHA models[J].Journal of Hohai University(Natural Sciences),2017,45(1):1-6.(in Chinese))

[ 7 ] 章國勇,伍永剛,楊林明,等.基于參數(shù)優(yōu)化的EEMD-LSSVM年徑流組合預(yù)測模型[J].水資源與水工程學(xué)報,2013,24(6):1-5.(ZHANG Guoyong,WU Yonggang,YANG Linming,et al.EEMD-LSSVM model of combined forecast of annual runoff based on parameter optimization[J].Journal of Water Resources & Water Engineering,2013,24(6):1-5.(in Chinese))

[ 8 ] ZHAO X,CHEN X.Auto regressive and ensemble empirical mode decomposition hybrid model for annual runoff forecasting[J].Water Resources Management,2015,29(8): 2913-2926.

[ 9 ] WANG W,CHAU K,QIU L,et al.Improving forecasting accuracy of medium and long-term runoff using artificial neural network based on EEMD decomposition[J].Environmental Research,2015,139: 46-54.

[10] 馬超,姜璇.基于EEMD-ANN的水庫年徑流預(yù)測[J].水電能源科學(xué),2016,34(8):32-35.(MA Chao,JIANG Xuan.Annual reservoir runoff forecasting using ensemble empirical mode decomposition and artificial neuron network[J].Water Resources and Power,2016,34(8):32-35.(in Chinese))

[11] 鄭曉東,魯帆,嚴登華,等.基于小波分析的丹江口水庫水文年度入庫徑流變化特征[J].水科學(xué)與工程技術(shù),2011(2):1-3.(ZHENG Xiaodong,LU Fan,YAN Denghua,et al.Variation characteristics of the Danjiangkou Reservoir hydrological annual runoff based on the morlet wavelet[J].Water Sciences and Engineering Technology,2011(2):1-3.(in Chinese))

[12] 李響,鄒寧,呂孫云.丹江口水庫天然入庫徑流演變規(guī)律研究[J].南水北調(diào)與水利科技,2015,13(2): 206-210.(LI Xiang,ZOU Ning,Lü Sunyun.Variation of natural inflow runoff in the Danjiangkou Reservoir[J].South-to-North Transfers and Water Science & Technology,2015,13(2): 206-210.(in Chinese))

[13] 冉篤奎,李敏,武晟,等.丹江口水庫中長期徑流量的多模型預(yù)報結(jié)果分析及綜合研究[J].水利學(xué)報,2010,41(9): 1069-1073.(RAN Dukui,LI Min,WU sheng,et al.Research on multi-model forecasts in mid-long term runoff in Danjiangkou Reservoir[J].Journal of Hydraulic Engineering,2010,41(9): 1069-1073.(in Chinese))

[14] 王元超.丹江口水庫中長期徑流預(yù)報及水質(zhì)水量聯(lián)合模擬技術(shù)[D].北京:中國水利水電科學(xué)研究院,2015.

[15] 馮小沖,王銀堂,劉勇,等.基于物理統(tǒng)計方法的丹江口水庫月入庫徑流預(yù)報[J].河海大學(xué)學(xué)報(自然科學(xué)版),2011,39(3):242-247.(FENG Xiaochong,WANG Yintang,LIU Yong,et al.Monthly runoff forecast for Danjiangkou Reservoir based on physical statistical methods[J].Journal of Hohai University(Natural Sciences),2011,39(3): 242-247.(in Chinese))

[16] WU Z,HUANG N E.Ensemble empirical mode decomposition: a noise-assisted data analysis method[J].Advance in Adaptive Data Analysis,2009,1(1): 1-41.

[17] HUANG N E,SHEN Z,LONG S R,et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings of the Royal Society A:Mathematical Physical and Engineering Sciences,1998,454: 903-995.

[18] 王文圣,丁晶,金菊良,等.隨機水文學(xué)[M].北京:中國水利水電出版社,2008:41,165-167.

StochasticsimulationandpredictionofannualrunoffintheDanjiangkouReservoirbasedonEEMD-ARmodel

LIAN Jijian, SUN Xiaozhong, MA Chao, ZHAO Ming, TANG Zhibo

(StateKeyLaboratoryofHydraulicEngineeringSimulationandSafety,TianjinUniversity,Tianjin300072,China)

Based on the analysis and identification of the annual runoff sequence components of the Danjiangkou Reservoir, deterministic components such as the trend term, the jumping term and the periodic term were derived by using linear trend regression analysis method, sequential cluster method and variance spectrum method, etc. A stochastic auto-regression model of annual runoff based on Ensemble Empirical Mode Decomposition (EEMD) was proposed (EEMD-AR) and it was applied to the stochastic simulation and prediction of the annual runoff in the Danjiangkou Reservoir. Through the EEMD decomposition, the problem that stochastic simulation and prediction by auto-regression (AR) model cannot be directly applied due to the non-stationary historical runoff sequence of the Danjiangkou Reservoir has been solved. The simulation results show that EEMD-AR model can simulate and predict the annual runoff sequence of the Danjiangkou Reservoir in a good forecast accuracy and it maintain the statistical characteristics of the original historical sequence.

Runoff sequence components identification; EEMD-AR model; runoff stochastic simulation; Danjiangkou Reservoir

國家重點研發(fā)計劃水資源高效利用專項(2016YFC0402203)

練繼建(1965—),男,教授,主要從事水利水電工程安全研究。E-mail: jjlian@tju.edu.cn

馬超(1981—),男,副教授,主要從事水文水資源研究。E-mail: mac_tju@126.com

10.3880/j.issn.1006-7647.2017.05.003

TV124

:A

:1006-7647(2017)05-0016-06

2016-10-19 編輯:駱 超)

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 欧美97色| 这里只有精品免费视频| 亚洲日韩AV无码一区二区三区人| 国产欧美日韩va另类在线播放| 国产欧美在线视频免费| 亚洲欧美另类日本| 亚洲国产欧美国产综合久久| 亚洲男人天堂2020| 国产内射一区亚洲| 午夜不卡视频| 亚洲中文字幕久久精品无码一区| 欧美精品亚洲精品日韩专区| 亚洲高清日韩heyzo| 久久毛片网| 欧美视频免费一区二区三区| 午夜欧美理论2019理论| 啪啪免费视频一区二区| 国产在线八区| 狠狠色综合网| 一区二区在线视频免费观看| 91福利在线观看视频| 人妻精品久久无码区| 一区二区无码在线视频| 精品视频免费在线| 日韩 欧美 国产 精品 综合| 国产波多野结衣中文在线播放| 成人国产小视频| 亚洲资源站av无码网址| 午夜毛片免费看| 不卡午夜视频| 成年人午夜免费视频| 超清无码熟妇人妻AV在线绿巨人| 91视频日本| 成年午夜精品久久精品| 国产一区二区三区免费| 成人在线观看不卡| 在线亚洲精品福利网址导航| 性做久久久久久久免费看| 99久久精彩视频| 中文字幕在线看| 国产精品无码一二三视频| 2021最新国产精品网站| 黄色免费在线网址| 免费aa毛片| 国产chinese男男gay视频网| 亚洲国产精品美女| 国产精品青青| 久久综合国产乱子免费| 午夜免费视频网站| 久久国产拍爱| 欧美爱爱网| 久久久久久午夜精品| 欧美不卡视频在线观看| 真实国产精品vr专区| 免费无码AV片在线观看国产| 456亚洲人成高清在线| 国产最爽的乱婬视频国语对白| 亚洲国内精品自在自线官| 日日拍夜夜操| 91www在线观看| 亚洲一区二区三区中文字幕5566| 亚洲色精品国产一区二区三区| 国产一在线| 91精品国产情侣高潮露脸| 91原创视频在线| 97色婷婷成人综合在线观看| 国产www网站| 日a本亚洲中文在线观看| 亚洲色欲色欲www在线观看| 在线观看亚洲人成网站| 国产乱子伦一区二区=| 国产91无码福利在线| 亚洲一区精品视频在线| 久爱午夜精品免费视频| 国产一级毛片yw| 在线播放精品一区二区啪视频| 亚洲欧美另类日本| 国模极品一区二区三区| 国产呦精品一区二区三区下载 | 国产精品成人免费视频99| 狠狠操夜夜爽| 激情亚洲天堂|