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

        ?

        基于SARIMA-GS-SVR組合模型的短期電力需求預(yù)測(cè)

        2022-07-17 06:04:34王萬(wàn)雄
        電子科技 2022年8期
        關(guān)鍵詞:殘差預(yù)測(cè)函數(shù)

        劉 晗,王萬(wàn)雄

        (甘肅農(nóng)業(yè)大學(xué) 理學(xué)院,甘肅 蘭州 730070)

        目前,工廠、企業(yè)規(guī)模的逐步擴(kuò)大和生活水平的提高使得電力消費(fèi)直線增長(zhǎng)。電力生產(chǎn)的控制和最佳規(guī)劃需要精確的電力需求預(yù)測(cè)。電力預(yù)測(cè)時(shí)期包括短期、中期和長(zhǎng)期,其中短期預(yù)測(cè)顯示了天、周和季節(jié)的系統(tǒng)變化,是該領(lǐng)域的研究熱點(diǎn)[1-3]。此外,組合模型利用了單一模型的優(yōu)勢(shì),提供了比單一模型更穩(wěn)定、可靠的預(yù)測(cè)結(jié)果,具有更重要的研究?jī)r(jià)值[4],已在經(jīng)濟(jì)、工業(yè)等領(lǐng)域有了廣泛的應(yīng)用[5-8]。文獻(xiàn)[9]應(yīng)用混合深度神經(jīng)網(wǎng)絡(luò)(Convolutional Long Short-Term Memory,CLSTM)模型預(yù)測(cè)短期電力需求。文獻(xiàn)[10]將小波變換(Wavelet Transform,WT)和人工神經(jīng)網(wǎng)絡(luò)(Artificial Neural Network,ANN)混合使用,提高了預(yù)測(cè)精度。文獻(xiàn)[11]利用前饋(Back Propagation,BP)神經(jīng)網(wǎng)絡(luò)和模糊系統(tǒng)的方法進(jìn)行了電力負(fù)荷預(yù)測(cè)。相比其它算法,神經(jīng)網(wǎng)絡(luò)算法在大樣本條件下性能優(yōu)良,隨著數(shù)據(jù)量的增多,其計(jì)算能力也更強(qiáng)。但是,神經(jīng)網(wǎng)絡(luò)算法可解釋性低,計(jì)算成本昂貴,易陷入局部最優(yōu),因此有學(xué)者提出了支持向量機(jī)(Support Vector Machine,SVM)預(yù)測(cè)模型。

        與神經(jīng)網(wǎng)絡(luò)算法相比,SVM算法的泛化能力強(qiáng),更容易收斂到全局最優(yōu)[12]。本文應(yīng)用傳統(tǒng)算法與機(jī)器學(xué)習(xí)算法相結(jié)合的方法來(lái)預(yù)測(cè)短期電力需求時(shí)間序列。由于受溫度、工作日等因素影響,短期電力需求時(shí)間序列呈現(xiàn)明顯的周期性特點(diǎn),而季節(jié)差分自回歸移動(dòng)平均(Seasonal Auto Regressive Integrated Moving Average,SARIMA)模型在處理周期性數(shù)據(jù)上表現(xiàn)良好[13]。用電量與人口、電價(jià)等諸多因素有關(guān),這些都會(huì)導(dǎo)致電力需求時(shí)間序列呈現(xiàn)高度的非線性特征,而支持向量回歸(Support Vector Regression,SVR)模型則在處理非線性數(shù)據(jù)上有較大的優(yōu)勢(shì)[12]。綜上所述,本文充分考慮了這兩個(gè)模型的優(yōu)點(diǎn)。此外,鑒于SVR模型參數(shù)的選擇將直接影響其預(yù)測(cè)精度,為解決傳統(tǒng)SVR因人工經(jīng)驗(yàn)設(shè)置參數(shù)產(chǎn)生的局部最優(yōu)問(wèn)題,本文選用網(wǎng)格搜索(Grid Search,GS)算法對(duì)SVR的懲罰因子C和核函數(shù)參數(shù)γ進(jìn)行全局尋優(yōu)。綜合考慮以上因素,本文最終建立了GS算法優(yōu)化的SARIMA-SVR模型,即SARIMA-GS-SVR組合模型,并建立了SARIMA、SVR、GS-SVR與指數(shù)平滑模型共同預(yù)測(cè)短期電力需求時(shí)間序列。

        1 基本理論

        1.1 SARIMA的基本理論

        自回歸移動(dòng)平均模型ARMA(p,q)描述平穩(wěn)的時(shí)間序列為

        xt=θ1xt-1+θ2xt-2+…+θpxt-p+εt+

        ω1εt-1+ω2εt-2+…+ωqεt-q

        (1)

        式中,θ1,θ2,…,θp和ω1,ω2,…,ωq為模型參數(shù);εt為誤差項(xiàng)。式(1)由兩部分組成

        xt1=θ1xt1-1+θ2xt1-2+…+θpxt1-p+εt1

        (2)

        xt2=εt2+ω1εt2-1+ω2εt2-2+…+ωqεt2-q

        (3)

        式(2)為p階自回歸模型AR(p),θ1,θ2,…,θp為自回歸系數(shù),簡(jiǎn)記為

        (4)

        式中,B為時(shí)間后移算子。

        式(3)為q階移動(dòng)平均模型MA(q),ω1,ω2,…,ωq為移動(dòng)平均系數(shù),簡(jiǎn)記為式(5)。

        xt2=Wq(B)εt2,

        (5)

        不平穩(wěn)的序列需要通過(guò)差分處理,平穩(wěn)后使用ARMA模型建模。令φ(B)=Bp(B)(1-B)d,得到

        Wq(B)εt2=θp(B)?dxt1

        (6)

        式中,d為差分次數(shù),且?d=(1-B)d。此時(shí)模型為積累式自回歸移動(dòng)平均模型,記為ARIMA (p,d,q)。

        具有周期性的時(shí)間序列需要先處理其周期性,之后用ARIMA模型建模,模型簡(jiǎn)寫(xiě)為

        θp(B)ψp(Bs)?d(1-Bs)Dxt1=Wq(B)φQ(Bs)εt2

        (7)

        其中

        ψp(Bs)=1-ψ1Bs-ψ2B2s-…-ψpBps

        (8)

        φQ(Bs)=1-φ1Bs-φ2B2s-…-φQBqs

        (9)

        式中,s為周期長(zhǎng)度;D為季節(jié)差分的次數(shù),模型記為ARIMA (p,d,q)(P,D,Q)s。

        1.2 SVR的基本理論

        SVM是以統(tǒng)計(jì)理論為基礎(chǔ)的機(jī)器學(xué)習(xí)算法,包括支持向量分類和SVR[14],是一種基于結(jié)構(gòu)風(fēng)險(xiǎn)最小化的算法。SVM有良好的泛化能力[15],在解決非線性問(wèn)題上效果良好,因此在分類和回歸問(wèn)題中得到了廣泛的運(yùn)用[16]。SVR是基于SVM的回歸算法,它的基本理論是通過(guò)核函數(shù)把數(shù)據(jù)從樣本空間映射到高維特征空間,然后在特征空間中進(jìn)行線性回歸,尋找最優(yōu)的回歸超平面,進(jìn)而實(shí)現(xiàn)數(shù)據(jù)預(yù)測(cè)的目的。

        假設(shè)訓(xùn)練樣本集S={(si,yi)},i=1,2,…,n,xi∈Rd(R為實(shí)數(shù)域,d為維數(shù)),xi為輸入,yi為對(duì)應(yīng)輸出。SVR的目標(biāo)為搜尋回歸函數(shù)f(x),使f(xi)與yi的偏差盡可能小,回歸函數(shù)為

        f(x)=ωTφ(x)+b

        (10)

        式中,ω為權(quán)值系數(shù);φ(x)為非線性變換函數(shù);b為偏執(zhí)項(xiàng)?;貧w預(yù)測(cè)的偏差不可避免,因此需要引入不敏感損失函數(shù)

        L[y,f(x)]=max{0,|y-f(x)|-ε}

        (11)

        式中,L[y,f(x)]為ε損失函數(shù),若預(yù)測(cè)偏差小于ε,則損失為0,否則將偏差減去ε。

        (12)

        (13)

        (14)

        1.3 指數(shù)平滑法的基本理論

        Holt-Winters模型是一種三參數(shù)指數(shù)平滑法,用來(lái)預(yù)測(cè)有水平項(xiàng)、趨勢(shì)項(xiàng)和季節(jié)性波動(dòng)的時(shí)間序列[18]。該方法提供了λ、β和γ3個(gè)參數(shù),分別對(duì)應(yīng)當(dāng)前點(diǎn)的水平、趨勢(shì)和季節(jié)部分。平滑參數(shù)λ控制水平的指數(shù)型下降,β控制斜率的指數(shù)型下降,γ控制季節(jié)指數(shù)下降。參數(shù)的取值范圍均為[0~1],取值越大說(shuō)明越近的觀測(cè)值占有的權(quán)重越大。t+1時(shí)刻的預(yù)測(cè)值可表示為

        (15)

        式中,a為t+1時(shí)刻序列的水平項(xiàng)系數(shù);b為斜率;xt為當(dāng)前時(shí)刻值;st+1為t+1時(shí)刻序列的季節(jié)效應(yīng)。

        2 數(shù)據(jù)采集及分析

        本文實(shí)驗(yàn)使用的數(shù)據(jù)來(lái)自美國(guó)國(guó)家航空和宇宙航空局(NASA,http://www.eia.doe.gov/)發(fā)布的加利福尼亞州(加州)電力需求歷史數(shù)據(jù),樣本數(shù)據(jù)時(shí)間范圍為2020年2月16日~2020年2月26日,共計(jì)264組,監(jiān)測(cè)間隔為1小時(shí)。圖 1為該時(shí)間段內(nèi)加州小時(shí)電力需求數(shù)據(jù)圖。由圖1可知,加州電力需求序列呈現(xiàn)周期性特征,在每天約7時(shí)和約19時(shí)達(dá)到峰值,在約13時(shí)達(dá)到谷值。原始序列分解為趨勢(shì)項(xiàng)、周期項(xiàng)和殘差項(xiàng),具體如圖2所示。圖2表明序列有很強(qiáng)的周期性特征和一定的趨勢(shì)特征。2020年2月22日與23日是周末,趨勢(shì)圖和殘差圖在這兩天出現(xiàn)下跌情況,原因可能為周末休息等因素。

        圖1 加州實(shí)際電力需求時(shí)序圖Figure 1. Sequence diagram of California actual electricity demand

        圖2 電力需求因素分解綜合圖Figure 2. Comprehensive decomposition diagram of power demand factors

        3 SARIMA-GS-SVR組合預(yù)測(cè)模型

        3.1 GS算法

        GS算法是通過(guò)遍歷給定的參數(shù)組合來(lái)優(yōu)化模型表現(xiàn)的方法[19]。本文設(shè)定懲罰因子C的變化區(qū)間為C∈[C1,C2],核函數(shù)參數(shù)γ的變化區(qū)間為γ∈[γ1,γ2]。實(shí)驗(yàn)使用GS算法對(duì)區(qū)間里的每對(duì)參數(shù)(C′,γ′)進(jìn)行模型訓(xùn)練,并計(jì)算均方誤差(Mean Square Error,MSE),最后采用MSE最小的參數(shù)作為最佳模型參數(shù)。算法步驟如下:

        步驟1初始化SVR參數(shù)C和γ,令C=C0,γ=γ0;

        步驟2設(shè)定C的步長(zhǎng)為Cs,更新式為C=C+Cs。γ的步長(zhǎng)為γs,更新式為γ=γ+γs;

        步驟3將參數(shù)(C,γ)帶入SVR構(gòu)建SVR模型,采用5折交叉驗(yàn)證法計(jì)算訓(xùn)練樣本5次交叉驗(yàn)證MSE的均值;

        圖3 GS算法流程圖Figure 3. Flow chart of GS algorithm

        步驟4將得到的MSE與上一步進(jìn)行對(duì)比,若MSE減小,則(C,γ)替換上一步的(C,γ),反之不操作;

        步驟5判斷C是否達(dá)到最大值C2,若是,進(jìn)行步驟6;反之,返回步驟2;

        步驟6判斷γ是否達(dá)到最大值γ2,若是,進(jìn)行步驟7;反之,返回步驟2;

        步驟7當(dāng)前SVR的參數(shù)(C,γ)達(dá)到最優(yōu),輸出結(jié)果。

        3.2 模型原理

        由前面數(shù)據(jù)分析可知原始序列包括線性時(shí)間序列分量、周期性時(shí)間序列分量和非線性時(shí)間序列分量,因此單一模型對(duì)電力需求預(yù)測(cè)有較大的局限。鑒于SARIMA模型對(duì)線性數(shù)據(jù)的擬合優(yōu)勢(shì)和SVR模型對(duì)非線性數(shù)據(jù)的擬合優(yōu)勢(shì),可利用GS算法對(duì)SVR模型的C和γ參數(shù)進(jìn)行尋優(yōu)。本文使用SARIMA-GS-SVR模型來(lái)預(yù)測(cè)短期電力需求。

        假設(shè)電力需求序列Yt由線性部分Lc和非線性部分Nt組成,即

        Yt=Lt+Nt

        (16)

        SARIMA-GS-SVR組合模型預(yù)測(cè)步驟如下:

        步驟1將原始序列因素分解,觀測(cè)序列周期項(xiàng)、趨勢(shì)項(xiàng)和殘差項(xiàng)特征;

        步驟2序列分為訓(xùn)練子集和測(cè)試子集,對(duì)前者差分,觀測(cè)序列平穩(wěn)性;

        步驟3查看序列自相關(guān)圖和偏自相關(guān)圖,根據(jù)特征擬合SARIMA模型;

        步驟5將SARIMA模型預(yù)測(cè)的殘差,以3個(gè)歷史殘差作為輸入,下一個(gè)殘差作為輸出,轉(zhuǎn)化為有監(jiān)督學(xué)習(xí)數(shù)據(jù)集;

        步驟6歸一化數(shù)據(jù)集并將其劃分為訓(xùn)練集和測(cè)試集;

        步驟7利用GS算法尋優(yōu)到最佳SVR參數(shù)組合(C,γ);

        步驟8將訓(xùn)練集數(shù)據(jù)帶入SVR模型,利用最優(yōu)參數(shù)進(jìn)行模型訓(xùn)練;

        (17)

        3.3 模型設(shè)置

        3.3.1 核函數(shù)的選取

        構(gòu)建SVR模型包括參數(shù)確定和核函數(shù)選取。需要確定的參數(shù)有正則化參數(shù)C、核函數(shù)參數(shù)γ,這兩個(gè)非負(fù)參數(shù)對(duì)于模型準(zhǔn)確預(yù)測(cè)具有重要作用[20]。RBF核函數(shù)可直觀反映兩個(gè)數(shù)據(jù)的距離,比其它核函數(shù)效果更好,應(yīng)用更廣泛[21],因此本文選用如下RBF核函數(shù)

        (18)

        式中,γ為控制半徑,且γ=1/2σ2;σ2為核函數(shù)的方差。

        3.3.2 數(shù)據(jù)預(yù)處理

        為了使數(shù)據(jù)具有相同尺度的量綱以減少誤差,在進(jìn)行SVR預(yù)測(cè)之前,對(duì)數(shù)據(jù)進(jìn)行歸一化處理

        (19)

        式中,xstd為歸一化后的值;xmin和xmax分別為數(shù)據(jù)的最小值和最大值。

        3.4 模型預(yù)測(cè)流程

        將SARIMA模型預(yù)測(cè)的殘差基于GS算法尋優(yōu)到的最佳參數(shù)(C,γ)輸入到SVR中,建立SARIMA-GS-SVR組合模型進(jìn)行電力需求預(yù)測(cè)。具體流程如圖4所示。

        圖4 SARIMA-GS-SVR模型流程圖Figure 4. Flow chart of SARIMA-GS-SVR model

        4 實(shí)驗(yàn)設(shè)置及結(jié)果分析

        4.1 實(shí)驗(yàn)參數(shù)設(shè)置

        基于Python 3仿真環(huán)境,以加州電力需求數(shù)據(jù)為例,建立SARIMA、SVR、GS-SVR、SARIMA-GS-SVR和指數(shù)平滑5種預(yù)測(cè)模型。設(shè)置GS算法參數(shù):模型為SVR,交叉驗(yàn)證折數(shù)CV=5,訓(xùn)練核數(shù)n_jobs=4,評(píng)分目標(biāo)函數(shù)為訓(xùn)練樣本5折交叉驗(yàn)證的均方誤差MSE,設(shè)置核函數(shù)參數(shù)γ和懲罰因子C的搜索上限均為28,下限為2-8,搜索步長(zhǎng)s=2.5,不敏感損失因子ε在尋優(yōu)過(guò)程中保持0.1不變。設(shè)置SVR模型參數(shù):選擇RBF徑向基核函數(shù),核函數(shù)參數(shù)γ和懲罰因子C的取值為GS算法尋優(yōu)到的最佳參數(shù)組合。設(shè)置指數(shù)平滑模型參數(shù):λ=0.951 2,β=0.935 0,γ=1,a=26 482.129 4,b=-712.614 9。

        4.2 實(shí)驗(yàn)結(jié)果及分析

        選取2020年2月16日~2020年2月25日范圍內(nèi)的240組數(shù)據(jù)作為預(yù)測(cè)模型的訓(xùn)練子集,并選取2020年2月26日內(nèi)的24組數(shù)據(jù)作為測(cè)試子集,進(jìn)行模型預(yù)測(cè)。由數(shù)據(jù)分析圖2看出,原始序列存在多種復(fù)雜特征,序列以天為周期,對(duì)訓(xùn)練子集進(jìn)行一階24步差分后所得序列如圖5所示。由圖5可知,去除個(gè)別極值點(diǎn)影響之外,序列均值和標(biāo)準(zhǔn)差在0值附近小范圍波動(dòng)。自相關(guān)圖與偏自相關(guān)圖如圖6所示,兩圖在一階延遲后都基本落在兩倍標(biāo)準(zhǔn)差范圍之內(nèi)。根據(jù)圖6,擬合模型SARIMA(1,1,1) (0,1, 0)24。該模型對(duì)測(cè)試子集的預(yù)測(cè)結(jié)果如圖8所示。

        圖5 電力需求差分序列時(shí)序圖Figure 5. Sequence diagram of power demand differential sequence

        圖6 訓(xùn)練樣本自相關(guān)與偏自相關(guān)圖Figure 6. Autocorrelation and partial autocorrelation graph of training sample

        表1列出了GS-SVR與SARIMA-GS-SVR模型利用GS算法尋優(yōu)到的最優(yōu)參數(shù)以及對(duì)應(yīng)的最小目標(biāo)函數(shù)MSE。將SARIMA模型預(yù)測(cè)的殘差利用GS算法尋優(yōu)到最佳參數(shù)(C=10.045 1,γ=3.590 8),帶入SVR模型得到2020年2月26日殘差預(yù)測(cè)結(jié)果,如圖7所示。

        表1 GS算法尋優(yōu)的最佳參數(shù)Table 1. Optimal parameters searched by GS algorithm

        由圖7可以看出,GS-SVR模型預(yù)測(cè)的殘差與實(shí)際殘差走勢(shì)一致,表明該模型具有良好的預(yù)測(cè)效果。SARIMA、SVR、GS-SVR、SARIMA-GS-SVR和指數(shù)平滑5種預(yù)測(cè)模型對(duì)測(cè)試子集的預(yù)測(cè)結(jié)果如圖8所示。從圖8 可以看出,相比其它4種模型的預(yù)測(cè)值,SARIMA-GS-SVR的預(yù)測(cè)值更接近實(shí)際值。SVR單一模型的預(yù)測(cè)效果差于GS-SVR組合模型,說(shuō)明了利用GS算法進(jìn)行SVR模型參數(shù)尋優(yōu)的必要性。將5種模型的預(yù)測(cè)值與實(shí)際值作比較,結(jié)果如表2所示。

        圖7 實(shí)際殘差值與GS-SVR預(yù)測(cè)殘差值對(duì)比Figure 7. Comparison of actual residual and GS-SVR forecasting residual

        圖8 電力需求實(shí)際值與5種模型預(yù)測(cè)值對(duì)比Figure 8. Comparison of actual power demand values and predicted values of five models

        表2 模型預(yù)測(cè)結(jié)果Table 2. Prediction results of models

        續(xù)表

        4.3 實(shí)驗(yàn)誤差對(duì)比

        本文應(yīng)用平均絕對(duì)誤差(Mean Absolute Error,MAE)、平均絕對(duì)百分比誤差(Mean Absolute Percent Error,MAPE)、均方根誤差(Root Mean Square Error,RMSE)和組合預(yù)測(cè)相對(duì)單一預(yù)測(cè)模型提高的預(yù)測(cè)精度(Accuracy Improvement,AI)來(lái)評(píng)價(jià)模型效果。各指標(biāo)如式(20)~式(23)所示,其中N為樣本個(gè)數(shù),xt為真實(shí)值,t為預(yù)測(cè)值,s為單一預(yù)測(cè)模型絕對(duì)誤差之和,sc為組合預(yù)測(cè)模型絕對(duì)誤差之和。當(dāng)AI>0時(shí),組合模型的預(yù)測(cè)效果優(yōu)于單一模型;反之,組合模型的預(yù)測(cè)效果較差。

        (20)

        (21)

        (22)

        (23)

        表3列出了5種模型的預(yù)測(cè)誤差。相比SVR,GS-SVR的預(yù)測(cè)精度提高了35.622 7%。相比SARIMA,SARIMA-GS-SVR的預(yù)測(cè)精度提高了29.181 2%。SARIMA-GS-SVR的MAE、MAPE和RMSE的誤差指標(biāo)評(píng)價(jià)值均小于其他4種模型。誤差評(píng)價(jià)結(jié)果表明,組合模型SARIMA-GS-SVR模型的預(yù)測(cè)效果最好。

        表3 模型預(yù)測(cè)誤差Table 3. Forecasting error of models

        5 結(jié)束語(yǔ)

        本文基于SARIMA模型預(yù)測(cè)的殘差,提出用GS算法對(duì)SVR的參數(shù)進(jìn)行優(yōu)化,并建立了SARIMA-GS-SVR組合預(yù)測(cè)模型,以美國(guó)加利福尼亞州2020年2月26日實(shí)際歷史數(shù)據(jù)做預(yù)測(cè)分析。 結(jié)果表明,SARIMA-GS-SVR組合預(yù)測(cè)模型預(yù)測(cè)效果優(yōu)于GS-SVR、SARIMA、SVR和指數(shù)平滑4種模型。本文提出的SARIMA-GS-SVR組合預(yù)測(cè)模型在一定程度上提高了短期電力需求的預(yù)測(cè)精度,但在SVR參數(shù)優(yōu)化算法的選擇上仍有一定的提升空間。在未來(lái)的研究中,計(jì)劃選擇其他算法來(lái)優(yōu)化SVR參數(shù),并進(jìn)行組合預(yù)測(cè)模型分析來(lái)繼續(xù)優(yōu)化預(yù)測(cè)精度。

        猜你喜歡
        殘差預(yù)測(cè)函數(shù)
        無(wú)可預(yù)測(cè)
        黃河之聲(2022年10期)2022-09-27 13:59:46
        選修2-2期中考試預(yù)測(cè)卷(A卷)
        選修2-2期中考試預(yù)測(cè)卷(B卷)
        基于雙向GRU與殘差擬合的車輛跟馳建模
        二次函數(shù)
        第3講 “函數(shù)”復(fù)習(xí)精講
        二次函數(shù)
        基于殘差學(xué)習(xí)的自適應(yīng)無(wú)人機(jī)目標(biāo)跟蹤算法
        函數(shù)備考精講
        基于遞歸殘差網(wǎng)絡(luò)的圖像超分辨率重建
        无码三级在线看中文字幕完整版| 精品无码久久久久成人漫画| 亚洲女同免费在线观看| 中文字幕国产精品中文字幕| 日日噜狠狠噜天天噜av| 国产精品久久国产精品99 gif| 国产主播一区二区三区在线观看| 亚洲国产成人手机在线电影| 内射少妇36p九色| 成人av一区二区三区四区 | 波多野结衣中文字幕一区二区三区| 人妻在线有码中文字幕 | 中文字幕精品人妻在线| 日韩亚洲一区二区三区在线| 国内免费AV网站在线观看| 大地资源网高清在线播放| 丰满少妇在线播放bd| 在线视频自拍视频激情| 亚洲最新中文字幕一区| 色伊人国产高清在线| 国产精品第一二三区久久蜜芽 | 伊人不卡中文字幕在线一区二区| 啪啪无码人妻丰满熟妇| 毛茸茸性xxxx毛茸茸毛茸茸| 精品淑女少妇av久久免费| 97人妻人人做人碰人人爽| 精品人妻一区二区三区久久| 日本一区二区三级在线| 国产精品日本一区二区三区在线| 国产区高清在线一区二区三区| 不卡a v无码在线| 国内无遮码无码| 日本免费人成视频播放| 国产成人亚洲综合无码| 中国丰满熟妇av| 国产精品久久久久久久妇| 国产激情艳情在线看视频| 国产精品亚洲av无人区一区香蕉| 国产成人av三级在线观看韩国 | 无码中文字幕人妻在线一区二区三区| 久久99久久99精品中文字幕 |