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

基于小波包分解的AJS-GMDH月徑流時間序列預測研究

2022-07-08 08:36:34楊瓊波崔東文
水力發(fā)電 2022年6期
關鍵詞:模型

楊瓊波,崔東文

(1.云南省水文水資源局紅河分局,云南 紅河 661199;2.云南省文山州水務局,云南 文山 663000)

主要從事水文水資源分析研究和水資源管理保護等工作;崔東文(通訊作者).

徑流時間序列預測是水文預報研究領域的熱點和難點。提高徑流時間序列預測精度對區(qū)域水資源開發(fā)利用、防洪抗旱規(guī)劃、水資源管理保護、水庫優(yōu)化調(diào)度等具有重要意義。由于徑流時間序列人類活動、受氣候變化、土地利用及植被覆蓋等多重因素的影響,常表現(xiàn)出多尺度、非線性、非平穩(wěn)性等特征,基于“分解—預測—重構”模式的預測方法被廣泛應用于徑流時間序列預測[1- 4],并取得了較好的預測效果。目前,常見的時間序列分解方法有經(jīng)驗模態(tài)分解(Empirical mode decomposition,EMD)、完全集合經(jīng)驗模態(tài)分解(Complete ensemble empirical mode decomposition,CEEMD)、小波分解(Wavelet decomposition,WD)、奇異譜分解(Singular spectrum analysis,SSA)、變分模態(tài)分解(Variational mode decomposition,VMD)等;預測模型有BP神經(jīng)網(wǎng)絡、廣義回歸神經(jīng)網(wǎng)絡(GRNN)、支持向量機(SVM)、相關向量機(RVM)、極限學習機(ELM)、LSTM神經(jīng)網(wǎng)絡、秩次集對分析法等。

WD是一種被廣泛應用的時頻分析工具,主要通過伸縮平移運算對信號逐步進行多尺度細化,最終達到高頻處時間細分,低頻處頻率細分的目的,已在各領域得到廣泛應用。然而,WD不足之處在于其僅對信號低頻部分詳細分解而忽略了信號的高頻部分。小波包分解(Wavelet packet decomposition,WPD)源于WD,它在分解信號低頻子集同時,對高頻子集繼續(xù)分解,且能夠根據(jù)分解需要自行設定分解層數(shù)和分解使用的小波函數(shù)[5-7],已在各種時序數(shù)據(jù)分解中得到廣泛應用。人工水母搜索(Artificial Jellyfish Search,AJS)算法是Chou等人于2020年基于水母覓食行為而提出一種新型群體智能算法,已在工程結構優(yōu)化中得到應用[8]。數(shù)據(jù)分組處理方法(Group method of data handling,GMDH)是針對復雜系統(tǒng)自組織建模的一種前饋型多層神經(jīng)網(wǎng)絡,網(wǎng)絡僅通過遞階分層取舍神經(jīng)元和調(diào)整閾值,實現(xiàn)對復雜非線性規(guī)律進行擬合,目前已在各種預測研究中得到應用[9-11]。在實際應用中,GMDH網(wǎng)絡權值、多項式系數(shù)、網(wǎng)絡層最大神經(jīng)元數(shù)等參數(shù)對GMDH網(wǎng)絡預測性能有著重要影響,目前主要有GMDH網(wǎng)絡權值優(yōu)化[12]、多項式系數(shù)優(yōu)化[13]兩種途徑,而鮮見于針對GMDH網(wǎng)絡層最大神經(jīng)元數(shù)、最大網(wǎng)絡層數(shù)和學習率等關鍵參數(shù)的優(yōu)化。

本文以云南省龍?zhí)墩?952年1月~2016年12月月徑流時間序列預測為例,基于WPD、AJS與GMDH網(wǎng)絡構建WPD-AJS-GMDH徑流時間序列預測模型。內(nèi)容安排如下:①利用WPD及對比分解方法CEEMD、WD將月徑流時間序列數(shù)據(jù)分解為若干子序列分量;②介紹AJS算法,在不同維度條件下通過6個典型函數(shù)對AJS算法進行仿真驗證,采用AJS優(yōu)化GMDH網(wǎng)絡層最大神經(jīng)元個數(shù)、最大網(wǎng)絡層數(shù)、學習率、選擇壓力系數(shù),構建WPD-AJS-GMDH模型,并構建WPD-AJS-SVM、WPD-AJS-BP、WPD-GMDH、WPD-SVM、WPD-BP、CEEMD-AJS-GMDH、CEEMD-AJS-SVM、CEEMD-AJS-BP、CEEMD-GMDH、CEEMD-SVM、CEEMD-BP和WD-AJS-GMDH、WD-AJS-SVM、WD-AJS-BP、WD-GMDH、WD-SVM、WD-BP作對比分析模型;③利用所構建的18種模型對實例月徑流時間序列進行檢驗。

1 研究方法

1.1 小波包分解(WPD)

WPD是在小波分解(WD)的基礎上改進和發(fā)展而來的一種有效分解方法,能同時對高頻、低頻部分實施信號分解。理論上,3層小波包分解能夠提取絕大部分信號中的有效信息,逼近任意非線性函數(shù),從而達到解決實際問題的目的。小波包分解算法公式[14-16]為

(1)

重構算法為

(2)

1.2 人工水母搜索(AJS)算法

1.2.1 AJS算法數(shù)學描述

AJS算法是基于水母覓食行為而提出一種新型元啟發(fā)式算法。該算法通過模擬水母洋流移動、群內(nèi)移動(主動移動和被動移動)和移動切換策略來達到求解優(yōu)化問題的目的。AJS算法基于3個理想化規(guī)則:①水母要么跟隨洋流移動,要么在群內(nèi)移動,“時間控制策略”控制著兩種移動之間的切換;②水母在海洋中尋找食物,他們更容易吸引到食物豐富的區(qū)域;③找到的食物數(shù)量取決于水母位置及其相應的目標函數(shù)值[7]。

AJS算法數(shù)學描述簡述如下[7]:

(1)洋流移動。AJS通過計算海洋中每只水母到當前最佳水母位置的平均矢量值確定洋流方向(趨勢)。 洋流移動數(shù)學描述

(3)

(4)

表1 AJS算法標準測試函數(shù)尋優(yōu)結果

式中,Xi(t+1)為第i只水母第t+1次迭代位置;Xi(t)為第i只水母第t次迭代位置;其他參數(shù)意義同上。

(2)水母群移動。在水母群中,水母分為被動移動(A型)和主動移動(B型)兩種類型。當水母群剛形成時,大多數(shù)水母呈A型移動;隨著時間推移,大多數(shù)水母呈B型移動。A型移動表示水母在其自身位置周圍移動,其位置更新公式為

Xi(t+1)=Xi(t)+γ×rand(0,1)×(Ub-Lb)

(5)

式中,γ為運動系數(shù),本文取0.1;Ub、Lb分別為搜索空間上、下限值;其他參數(shù)意義同上。B型運動描述為:當水母Xj所處區(qū)域食物量超過水母Xi所處區(qū)域食物量時,水母Xi向水母Xj移動;反知,水母Xi遠離水母Xj。其位置更新公式描述如下

(6)

式中,Xj(t)為第j只水母的第t次迭代位置;f(Xi)、f(Xj)分別為水母Xi、水母Xj的適應度值;其他參數(shù)意義同上。

(3)時間控制切換策略。水母在水母群中的運動分為A型和B型,為更好地模擬水母在這兩種移動之間的切換,引入時間控制函數(shù)c(t)來描述這種情況,當c(t)≥c0時,水母跟隨洋流移動;當c(t)

c(t)=|(1-t/Tmax)×[2×rand(0,1)-1]|

(7)

式中,t為當前迭代次數(shù);Tmax為算法最大迭代次數(shù);c0為時間控制常數(shù),本文取0.5。

1.2.2 AJS算法仿真驗證

在不同維度條件下,選取Sphere等6個測試函數(shù)對AJS算法進行仿真驗證(見表1)。AJS算法通過Matlab 2018a M語言實現(xiàn),利用函數(shù)20次尋優(yōu)平均值評估AJS算法的性能。設置AJS算法最大迭代次數(shù)Tmax=3 000,種群規(guī)模npop=100。

對于單峰函數(shù),AJS算法在不同維度條件下20次尋優(yōu)精度均在7.20×10-83以上,表現(xiàn)出較好的尋優(yōu)精度;對于多峰函數(shù),AJS算法在不同維度條件下20次尋優(yōu)精度均在1.22×10-12以上,具有較好全局尋優(yōu)能力??梢?,對于表1中的函數(shù),在不同維度條件下AJS算法均表現(xiàn)出較好尋優(yōu)效果。

1.3 GMDH網(wǎng)絡

GMDH網(wǎng)絡是針對復雜系統(tǒng)自組織建模的一種前饋型多層神經(jīng)網(wǎng)絡,也稱多項式網(wǎng)絡,目前已在電主軸熱位移[17]、經(jīng)濟短期預測[18]、碳價格預測[19]、潮汐預報[20]、用戶流失預測[21]等方面得到應用。GMDH網(wǎng)絡由系統(tǒng)各輸入單元交叉組合產(chǎn)生活動神經(jīng)元,從已產(chǎn)生的一代神經(jīng)元中選擇若干與目標變量最接近的神經(jīng)元,通過強強結合再次產(chǎn)生新的神經(jīng)元,……。通過重復這一優(yōu)勝劣汰和逐漸進化的過程,直到獲得外準則值最低的最優(yōu)神經(jīng)元,即完成GMDH網(wǎng)絡模型訓練過程[9-11]。

yi=f(xi1,xi2,…,xin),i=1,2,…,N

(8)

式中,yi為實際輸出;xi1,xi2,…,xin為輸入;N為輸入變量的數(shù)量。

(9)

(10)

在GMDH預測過程中,除網(wǎng)絡權值、多項式系數(shù)對GMDH預測效果產(chǎn)生直接影響外,GMDH網(wǎng)絡層最大神經(jīng)元數(shù)、最大網(wǎng)絡層數(shù)、學習率等關鍵參數(shù)對GMDH網(wǎng)絡預測性能同樣具有重要影響。為進一步提升GMDH網(wǎng)絡預測性能,拓展GMDH網(wǎng)絡優(yōu)化方向,本文利用AJS算法對GMDH網(wǎng)絡層最大神經(jīng)元數(shù)、最大網(wǎng)絡層數(shù)、學習率、選擇壓力系數(shù)進行優(yōu)化,以期進一步提高GMDH網(wǎng)絡預測性能。

1.4 建模流程

WPD-AJS-GMDH模型預測實現(xiàn)如下:

步驟一,利用WPD將原月徑流時間序列分解為8子序列分量[3,0]~[3,7]。利用自相關函數(shù)法(AFM)確定各子序列分量的輸入、輸出向量,合理劃分訓練、預測樣本[22],見圖1。

步驟二,構建訓練樣本均方誤差作為AJS-GMDH模型適應度函數(shù)

(11)

步驟三,設置AJS算法水母種群規(guī)模npop、最大迭代次數(shù)Tmax和算法終止條件。利用式(12)初始化水母空間位置Xi(i=1,2,…,npop),則

Xi+1=ηXi(1-Xi),0≤X0≤1

(12)

式中,Xi為第i只水母位置的logistic混沌值;X0用于生成水母的初始種群,X0∈(0,1),且X0?{0.00,0.25,0.50,0.75,1.00 };參數(shù)η設置為4.0。

步驟四,計算每只水母Xi適應度值,找到并保存最佳水母適應度fbest和最佳水母個體位置X*。

步驟五,利用式(7)計算時間控制函數(shù)c(t)。若c(t)≥0.5,利用式(4)更新水母位置;否則,利用式(5)、(6)更新水母位置。

步驟六,計算當前每只水母Xi適應度值;比較并保存當前最佳水母適應度值fbest和最佳水母個體位置X*。

步驟七,判斷算法是否達到終止條件,若是,則輸出最佳水母個體位置X*,算法結束。

步驟八,輸出最佳水母個體位置X*。即,GMDH網(wǎng)絡層最大神經(jīng)元數(shù)、最大網(wǎng)絡層數(shù)、學習率、選擇壓力系數(shù);利用AJS-GMDH模型對分量[3,0]~[3,7]進行預測,疊加重構后即得到實例月徑流最終預測結果。

步驟九,利用平均絕對百分比誤差(MAPE)、平均絕對誤差(MAE)、納什系數(shù)(NSE)和均方根誤差(RMSE)對各模型有效性進行評估。即

(13)

圖1 實例月徑流預測流程

2 實例應用

龍?zhí)墩疚挥谠颇鲜∥纳绞信手ㄦ?zhèn)龍?zhí)墩?,系紅河流域瀘江水系盤龍河干流控制站,控制徑流面積3 128 km2,為我國重要的水文站和中央報汛站。

2.1 時序數(shù)據(jù)分解

2.1.1 WPD分解

基于demy小波包基,利用WPD將龍?zhí)墩?952年1月~2016年12月逐月徑流數(shù)據(jù)進行3層小波包分解,得到8個子序列分量數(shù)據(jù),見圖2。

從圖2來看,[3,0]~[3,3]為低頻空間數(shù)據(jù)波形,振幅較大、波長較短,大致反映月徑流時間序列的變化趨勢;[3,4]、[3,5]、[3,6]、[3,7]為高頻空間數(shù)據(jù)波形,從[3,4]分量到[3,7]分量,其振幅逐漸減小、波長逐漸變長,大致反映月徑流時間序列的波動情況[4]。

2.1.2 CEEMD分解

利用CEEMD方法對龍?zhí)墩緦嵗聫搅鲾?shù)據(jù)進行分解,分解結果見圖3。從圖3可以看出,高頻子序列RSE的頻率最高、幅度最大、波長最短;中頻子序列IMF4~IMF8和低頻子序列IMF1~IMF3的振幅逐漸減小,頻率逐漸降低,波長逐漸增大,周期性逐漸增強,符合經(jīng)驗模態(tài)分解規(guī)律[18]。

2.1.3 WD分解

利用db5小波將龍?zhí)墩緦嵗聫搅鲾?shù)據(jù)進行分解,得到1個低頻率分量rA5和5個高頻率分量rD1、rD2、rD3、rD4、rD5,見圖4。從圖4可以看出,高頻分量反映原始數(shù)據(jù)的變化特征,低頻分量反映原始數(shù)據(jù)的變化趨勢。

圖2 實例月徑流時間序列WPD分解3D效果

圖3 實例月徑流時間序列CEEMD分解3D效果

圖4 實例月徑流時間序列WD分解3D效果

2.2 時間序列建模

本文采用自相關函數(shù)法(AFM)確定各分量的輸入、輸出向量,即通過計算各子序列分量的自相關系數(shù),在滯后數(shù)≥3情況下,選取自相關系數(shù)最大時所對應的滯后數(shù)作為各子序列分量最優(yōu)輸入維數(shù),通過最優(yōu)輸入維數(shù)構建預測模型輸入、輸出向量,見表2。本文利用實例前649~657個月徑流數(shù)據(jù)作為訓練樣本,后120個月作為預測樣本[4]。

表2 各分量自相關系數(shù)、輸入維數(shù)及序列長度

從表2可以看出,經(jīng)WPD分解的各子序列分量自相關系數(shù)絕對值均在0.591以上;經(jīng)CEEMD分解的IMF1分量、經(jīng)WD分解的rD1分量自相關系數(shù)絕對值小于0.25,這在很大程度上制約分解效果和預測精度。

2.3 參數(shù)設置及預測分析

2.3.1 參數(shù)設置

(1)WPD-AJS-GMDH、WPD-GMDH、CEEMD-AJS-GMDH、CEEMD-GMDH、WD-AJS-GMDH、WD-GMDH模型:設置WPD-GMDH、CEEMD-GMDH、WD-GMDH模型網(wǎng)絡層最大神經(jīng)元數(shù)m=20;最大網(wǎng)絡層數(shù)r=5;選擇壓力系數(shù)α=0.6;學習率p=0.85;WPD-AJS-GMDH、CEEMD-AJS-GMDH、WD-AJS-GMDH模型AJS算法最大迭代次數(shù)Tmax=50,種群規(guī)模npop=30;網(wǎng)絡層最大神經(jīng)元數(shù)m搜索空間[5,50],最大網(wǎng)絡層數(shù)r搜索空間[2,10],選擇壓力系數(shù)α搜索空間[0.01,0.9],學習率p搜索空間[0.01,0.9],其他采用默認值。

(2)WPD-AJS-SVM、WPD-SVM、CEEMD-AJS-SVM、CEEMD-SVM、WD-AJS-SVM、WD-SVM模型:WPD-SVM、CEEMD-SVM、WD-SVM模型懲罰因子、核函數(shù)參數(shù)和不敏感損失系數(shù)通過人工試湊法確定;WPD-AJS-SVM、CEEMD-AJS-SVM、WD-AJS-SVM模型AJS算法最大迭代次數(shù)Tmax=50,種群規(guī)模npop=30;不敏感損失系數(shù)、懲罰因子、核函數(shù)參數(shù)搜索空間分別設置為[0.001,1]、[0.01,100]、[0.01,100],輸入數(shù)據(jù)采用[0,1]進行歸一化處理。

(3)WPD-AJS-BP、WPD-BP、CEEMD-AJS-BP、CEEMD-BP、WD-AJS-BP、WD-BP模型:設置WPD-BP、CEEMD-BP、WD-BP模型隱層傳遞函數(shù)、輸出層傳遞函數(shù)、訓練函數(shù)分別為tansig、purelin、trainlm,設定期望誤差0.000 001,最大訓練輪回1 000次,隱層節(jié)點數(shù)為輸入維數(shù)的2倍減1;WPD-AJS-BP、CEEMD-AJS-BP、WD-AJS-BP模型AJS算法最大迭代次數(shù)Tmax=50,種群規(guī)模npop=30,權值閾值搜索空間∈[-1,1],BP部分設置同WPD-BP、CEEMD-BP、WD-BP模型。輸入數(shù)據(jù)采用[0,1]進行歸一化處理。

2.3.2 實例預測結果分析

基于表2構建各分量的輸入、輸出向量,并利用所建立的WPD-AJS-GMDH等18種模型對各分量進行預測,將預測結果疊加重構后即得到實例月徑流最終預測結果。各模型采用上述MAPE、MAE、RMSE、NSE進行評估,結果見表3,預測效果見圖5、6。

由表3及圖5、6可以得出以下結論:

(1)WPD-AJS-GMDH、WPD-AJS-SVM、WPD-AJS-BP模型的預測效果優(yōu)于WPD-GMDH、WPD-SVM、WPD-BP模型,遠優(yōu)于其他模型。其中,尤以WPD-AJS-GMDH模型預測的精度最高、效果最好,其平均絕對百分比誤差2.14%,平均絕對誤差0.25 m3/s,納什系數(shù)0.999 4,均方根誤差0.331 m3/s。將WPD-AJS-GMDH模型用于月徑流時間序列預測是可行和可靠的,模型具有更高的預測和更好的泛化能力。

(2)從不同分解方法預測結果對比來看,基于WPD分解的同種模型的預測效果要遠優(yōu)于基于CEEMD、WD分解的同種模型,表明WPD對月徑流時間序列的分解效果要優(yōu)于CEEMD、WD方法。

表3 實例月徑流時間序列預測結果對比

圖5 實例月徑流時間序列預測相對誤差3D

圖6 實例月徑流時間序列預測絕對誤差3D

(3)基于AJS算法優(yōu)化的同種模型的預測效果要優(yōu)于未經(jīng)AJS算法優(yōu)化的同種模型,表明AJS算法能有效優(yōu)化GMDH網(wǎng)絡層最大神經(jīng)元數(shù)、最大網(wǎng)絡層數(shù)、學習率和選擇壓力系數(shù),以及SVM懲罰因子、核函數(shù)參數(shù)、不敏感損失系數(shù)和BP網(wǎng)絡權閾值,有效提高GMDH網(wǎng)絡、SVM、BP網(wǎng)絡的預測性能。

(4)從圖5、6可看出,WPD-AJS-GMDH模型預測的絕大多數(shù)樣本絕對百分比誤差、平均絕對誤差均在0值附近波動,具有較優(yōu)的預測效果;次之為WPD-AJS-SVM、WPD-AJS-BP模型;再次為WPD-GMDH、WPD-SVM、WPD-BP模型,其他12種模型則均較差。

3 結 論

本研究將WPD、AJS與GMDH網(wǎng)絡進行結合,構建WPD-AJS-GMDH月徑流時間序列預測模型,利用云南省龍?zhí)墩?952年1月~2016年12月共780組月徑流時序數(shù)據(jù)進行實證分析,并建立基于CEEMD、WD分解的GMDH、SVM、BP模型做對比分析。結果如下:

(1)在不同維度條件下,AJS算法均具有較好的尋優(yōu)效果,將AJS算法用于GMDH網(wǎng)絡層最大神經(jīng)元數(shù)、最大網(wǎng)絡層數(shù)、學習率和選擇壓力系數(shù)尋優(yōu)是可靠的。

(2)WPD-AJS-GMDH模型對實例后120個月月徑流時間序列預測的平均絕對百分比誤差為2.14%,平均絕對誤差為0.25 m3/s,納什系數(shù)0.999 4,均方根誤差0.331 m3/s,預測效果略優(yōu)于WPD-AJS-SVM、WPD-AJS-BP模型,優(yōu)于WPD-GMDH、WPD-SVM、WPD-BP模型,遠優(yōu)于其他模型。將WPD-AJS-GMDH模型用于月徑流時間序列預測是可行和可靠的,模型具有更高的預測和更好的泛化能力。

(3)實例驗證表明,WPD能將徑流時序數(shù)據(jù)分解為更具規(guī)律的子序列分量,有效弱化復雜環(huán)境對徑流時間序列的影響,降低預測復雜度,其分解效果明顯優(yōu)于CEEMD、WD方法。

(4)實例驗證表明,AJS算法能有效優(yōu)化GMDH網(wǎng)絡層最大神經(jīng)元數(shù)、最大網(wǎng)絡層數(shù)、學習率和選擇壓力系數(shù),不但克服了GMDH網(wǎng)絡人工選擇關鍵參數(shù)的不足,而且大大提高了GMDH網(wǎng)絡的預測精度和智能化水平,具有較好的實用價值。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 久久国产成人精品国产成人亚洲 | 性69交片免费看| 国产簧片免费在线播放| 久久久久亚洲精品成人网| 永久在线播放| 啪啪免费视频一区二区| 日本人妻丰满熟妇区| 精品超清无码视频在线观看| 国产肉感大码AV无码| 国产成a人片在线播放| 久久无码av三级| 在线精品自拍| 91九色视频网| 欧美综合区自拍亚洲综合天堂| 激情国产精品一区| 97综合久久| 在线精品自拍| 人妻无码一区二区视频| a级毛片免费看| 91视频日本| 国产在线麻豆波多野结衣| 天天综合网站| 久久国产毛片| 999国内精品久久免费视频| 欧美色伊人| 中文字幕免费播放| 欧美成人午夜影院| 亚洲男人的天堂在线| 国产成人av一区二区三区| 国产精品无码在线看| 国产99视频精品免费观看9e| 福利一区三区| 波多野结衣一级毛片| 97无码免费人妻超级碰碰碰| 午夜日b视频| 欧美19综合中文字幕| 欧美a在线看| 免费无码AV片在线观看中文| 亚洲欧洲自拍拍偷午夜色无码| 日韩av高清无码一区二区三区| 国产亚卅精品无码| 国产精品专区第一页在线观看| 欧美成人精品在线| 97se亚洲综合在线韩国专区福利| 欧美激情伊人| 国产杨幂丝袜av在线播放| 亚洲日本在线免费观看| 亚洲精品视频免费观看| 欧美人与牲动交a欧美精品| 国产成a人片在线播放| 毛片一区二区在线看| 无码高清专区| www.精品国产| 亚州AV秘 一区二区三区| 国产欧美日韩在线一区| 欧美一区精品| 亚洲美女操| 亚洲中文久久精品无玛| 久久99国产精品成人欧美| 国产久草视频| 无码一区二区三区视频在线播放| 69综合网| 久久无码免费束人妻| 综合亚洲网| 在线无码av一区二区三区| 在线视频精品一区| 国产在线精品人成导航| 亚洲无码熟妇人妻AV在线| 91小视频在线播放| 欧美国产精品不卡在线观看| 日韩区欧美区| 91九色国产porny| 日韩毛片免费视频| 欧美在线中文字幕| 中国国产A一级毛片| 亚洲人成亚洲精品| 亚洲综合色婷婷| 日本a级免费| 国产成人区在线观看视频| 精品無碼一區在線觀看 | 99精品伊人久久久大香线蕉| 国产亚洲男人的天堂在线观看|