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

        ?

        基于深度時間卷積神經(jīng)網(wǎng)絡的風電功率預測

        2022-03-14 11:43:38王碩禾張嘉姍常宇健張國駒
        濟南大學學報(自然科學版) 2022年2期
        關鍵詞:電功率殘差分量

        劉 晗,王碩禾, 2,張嘉姍,常宇健, 2,張國駒

        (1.石家莊鐵道大學 電氣與電子工程學院,河北 石家莊 050043;2.河北省分布式能源應用技術創(chuàng)新中心,河北 石家莊 050043;3.中國科學院電工研究所,北京 100190)

        風能作為一種可再生能源,是新能源系統(tǒng)的重要組成部分,具有豐富性與環(huán)境友好性[1]。雖然我國近海地區(qū)風能豐富易得;但由于風能在氣象條件上的隨機性、波動性以及空間和時間上的不確定性,因此濱海地區(qū)風力發(fā)電的預測與調(diào)度相對于傳統(tǒng)能源變得較為困難[2]。

        近年來,國內(nèi)外學者針對風力發(fā)電預測的研究很多,取得了豐碩的成果[3-6]。風力發(fā)電的預測是利用歷史功率數(shù)據(jù)或者氣象數(shù)據(jù)對未來功率的產(chǎn)生進行建模的過程,主要建模方法分為物理、統(tǒng)計和組合預測的方法[7]。物理、統(tǒng)計方法分別在長期、短期預測中具有突出的優(yōu)勢[8],組合預測方法則將不同算法組合在一起,利用每種算法的優(yōu)勢來提高預測精度。

        模態(tài)分解技術可提取原始的風電功率序列所包含的特征,近年來被廣泛利用。經(jīng)驗模態(tài)分解(empirical mode decomposition,EMD)[9]算法將原始序列平穩(wěn)化處理,將其分解成有限個代表不同時間尺度上樣本特征的本征模態(tài)函數(shù)(intrinsic mode function,IMF),但是存在模態(tài)混疊的缺點。自適應集成經(jīng)驗模態(tài)分解(complete ensemble empirical mode decomposition with adaptive noise,CEEMDAN)[10]算法在前人的基礎上對相應的模型進行了改進,通過在原始信號中添加自適應的白噪聲,既能解決EMD算法的模式混淆問題,又能抵消白噪聲殘留。

        卷積神經(jīng)網(wǎng)絡是深度學習中常見的網(wǎng)絡結(jié)構(gòu),利用卷積核與輸入的卷積運算來提取特征信息,近年來已成功應用于模式識別[11]、圖像分類[12]、特征提取[13]、故障診斷[14]等領域。因果膨脹卷積依照輸入的時間順序?qū)v史數(shù)據(jù)進行建模,在不影響運算復雜度的情況下使卷積核的運算范圍擴大,使得卷積核的接受視野拓寬。深度時間卷積網(wǎng)絡(deep temporal convolutional networks, DTCN)[15]將因果膨脹卷積用于時間序列建模,經(jīng)過實驗驗證,時間卷積結(jié)構(gòu)在各種任務和數(shù)據(jù)集上的表現(xiàn)優(yōu)于常規(guī)循環(huán)神經(jīng)網(wǎng)絡,如長短時記憶網(wǎng)絡(long short-term memory,LSTM)。而在DTCN的訓練過程中,學習率動態(tài)調(diào)整是值得關注的問題,合適的學習率有助于DTCN的快速收斂。余弦退火[16](cosine annealing,CA)是學習率優(yōu)化策略的一種,與熱重啟技術相結(jié)合,在訓練后期可改善模型陷入局部最優(yōu)的狀態(tài)。

        本文中將CEEMDAN與DTCN引入風力發(fā)電功率預測中,CEEMDAN將風電功率分解成若干個模態(tài)函數(shù)用于特征分析,DTCN中的因果膨脹卷積對分解后的功率進行信息提取,避免了信息泄露。通過改進的余弦退火(improved cosine annealing,ICA)算法優(yōu)化DTCN學習率,使模型易于收斂,利用深度殘差結(jié)構(gòu)降低梯度彌散和梯度爆炸。本文中依據(jù)環(huán)渤海某地風電場數(shù)據(jù),建立8個對照預測模型,通過實驗驗證了本文中提出的模型的預測精度與準確性,為濱海地區(qū)風電功率預測提供參考。

        1 基于CEEMDAN的風電功率特征提取方法

        1.1 CEEMDAN原理

        風力發(fā)電的預測可以看作是時間序列建模的過程。以歷史序列作為輸入,t+1時刻的風電功率作為輸出進行預測;但這種建模過程往往學習到的是輸入到輸出的相關性,忽視了風電功率內(nèi)在特征對預測的作用,因此對歷史風電功率序列進行分解,有利于從信號特征上對風電功率進行分析。

        為了克服EMD過程中有用信號消失和模態(tài)混疊的問題,CEEMDAN可用于信號分析,以平滑高頻本征模態(tài)函數(shù)中的噪聲。假設x(n)為原始風電功率序列,其中n為采樣點,fk(n)為EMD得到的第k個模態(tài)分量,w(n)為滿足N(0,1)分布的白噪聲序列,第k次分解時設置其偏差為αk。定義Ei(·)通過EMD得到第i個模態(tài)分量的算子。CEEMDAN算法的流程如下。

        1)構(gòu)造信號序列,在原始風電功率序列中加入白噪聲。

        x′(n)=x(n)+αkw(n),

        (1)

        式中x′(n)為構(gòu)造的信號序列。

        2)計算第1個模態(tài)分量f1(n),并得到余項r1(n)。

        (2)

        r1(n)=x′(n)-f1(n),

        (3)

        3)利用r1(n)計算f2(n),則f2(n)和余項r2(n)分別為

        (4)

        r2(n)=r1(n)-f2(n)。

        (5)

        4)計算第k個模態(tài)函數(shù)分量fk(n)和第k個余項rk(n),分別為

        (6)

        rk(n)=rk-1(n)-fk(n)。

        (7)

        5)重復步驟4),得到K個IMF,最終的余項為

        (8)

        當殘差小于某個閾值后,迭代終止,最終得到若干個模態(tài)函數(shù)以及余項序列。不同的模態(tài)函數(shù)代表著原始風電功率信號在不同頻率下的特征。利用CEEMDAN進行信號分析,以提取原始風電功率序列內(nèi)在特征。

        1.2 基于排列熵的IMF重構(gòu)

        CEEMDAN產(chǎn)生的模態(tài)函數(shù)數(shù)量較多,如果全部輸入至預測網(wǎng)絡中,會帶來模型龐大、訓練時間過長的問題。文獻[17]中定義了一種衡量復雜度的定量計算方式——排列熵(permutation entropy,PE)。為了簡化問題,本文中將排列熵引入分析序列的復雜度,用于重構(gòu)IMF分量。排列熵的具體計算方式如下。

        1)設長度為T的序列為{y(t),t=1,2,…,T},對序列y(t)進行相空間重構(gòu),得到矩陣Y,

        (9)

        式中:l為重構(gòu)分量數(shù);m為嵌入維數(shù);τ為延遲因子。

        3)計算向量Y′(r)出現(xiàn)的概率P1,P2,…,Pl,排列熵Hp(m)可用香農(nóng)熵的形式定義為

        (10)

        4)為了方便表示,通常將排列熵Hp(m)歸一化表示為

        (11)

        由式(11)可得出,當Pr=1/m!時,Hp(m)有最大值ln(m!)。歸一化后的Hp在[0,1]范圍內(nèi),即

        0≤H(m)/ln(m!)≤1。

        (12)

        Hp反映了時間序列的復雜程度,其值越小說明時間序列越規(guī)則,反之則隨機性越大。將CEEMDAN產(chǎn)生的IMF分量利用PE算法計算出各個分量的熵值,可定量判斷各個分量的復雜程度,并將此作為依據(jù)對各個分量進行分類,對復雜度較低的分量進行合并,從而達到減少輸入維度、提高計算效率的目的。

        2 DTCN序列建模

        風電功率是基于時間變化的一維序列,而DTCN解決的是序列建模的問題,因此利用DTCN本身對于一維數(shù)據(jù)強大的處理能力,可以有效地對風電功率進行精準預測。

        2.1 風電功率序列建模

        時間卷積網(wǎng)絡對時間序列的建模方式如下:假定有一輸入序列x=(x0,x1,x2, …,xT),對應的期望輸出為y=(y0,y1,y2, …,yT),其中序列x、y要滿足因果關系,即預測時刻t的輸出yt,必須使用先前觀察得到輸入x0,x1,x2, …,xt-1。DTCN的建??梢允钱a(chǎn)生映射的任何函數(shù)f,表示為

        (13)

        風電功率序列的預測問題可表示為

        (14)

        構(gòu)建風電功率序列預測值與觀測值的損失函數(shù),本文中采用log-cosh函數(shù),定義為

        (15)

        2.2 DTCN的結(jié)構(gòu)

        2.2.1 因果膨脹卷積原理

        一維卷積可對時間序列進行運算并提取各種特征,但隨著時間序列長度的增長,規(guī)則的卷積網(wǎng)絡則需要更多的卷積層以接收更長的序列。膨脹卷積在傳統(tǒng)卷積上進行了改進,通過在每2個相鄰的卷積核之間引入一個固定的跳躍以擴大卷積核的運算范圍,從而降低運算的復雜度,提升運算的效率。層數(shù)為L且卷積核大小為k的接收域r的表達式為

        r=2L-1k。

        (16)

        因果卷積在時間維度上僅依賴于歷史數(shù)據(jù),保持因果性,且達到提取歷史數(shù)據(jù)特征的作用,實現(xiàn)對未來數(shù)據(jù)預測的能力。時間序列s上的因果膨脹卷積運算F定義為

        (17)

        式中:p=(p1,p2,…,pt)為輸入向量;*為膨脹卷積運算符;f為卷積核;k為卷積核大?。籶s-d·i為參與膨脹卷積運算的輸入向量元素,其中d為擴張因子,當d=1時,膨脹卷積轉(zhuǎn)換為規(guī)則卷積。

        因果膨脹卷積的結(jié)構(gòu)如圖1所示。由圖可以看出,具有3個擴張卷積層、膨脹率為4且卷積核大小為2的網(wǎng)絡接收域為8,即8個輸入值決定一個輸出。理論上,因果膨脹卷積可以接受任意時間步長的數(shù)據(jù)作為輸入。

        圖1 因果膨脹卷積結(jié)構(gòu)

        2.2.2 DTCN殘差模塊

        DTCN殘差模塊中采用殘差連接結(jié)構(gòu),即卷積運算后的網(wǎng)絡輸出與殘差通道進行元素相加運算,減少了反向傳播的過程中因網(wǎng)絡層數(shù)過多而導致的梯度彌散或者梯度爆炸。殘差結(jié)構(gòu)可繼續(xù)加深神經(jīng)網(wǎng)絡使模型具有更好的性能。殘差輸出o為

        o=φ(F(X)+X),

        (18)

        式中:F(X)為卷積運算模塊;X為殘差通道;φ為激活函數(shù),本文中選取的為修正線性單元(rectified linear units,Relu)函數(shù)。

        DTCN的殘差模塊結(jié)構(gòu)如圖2所示。DTCN的殘差模塊由殘差函數(shù)F(X)和X組成。F(X)通道中輸入經(jīng)過2次運算單元,每個單元首先經(jīng)過膨脹卷積層,后進入權值正則化層。權值正則化層目的是控制過擬合現(xiàn)象。然后將Relu函數(shù)作為卷積層的激活函數(shù),最后添加失活層。右側(cè)為殘差通道,由輸入經(jīng)過一個卷積核大小為1的卷積層,目的是使F(X)和X兩路輸出張量的維度保持一致。最終2個通道經(jīng)過一個逐元素相加網(wǎng)絡層運算,作為殘差單元的輸出。

        F(X)—卷積運算模塊;X—殘差通道;?—逐元素相加。

        DTCN的總體結(jié)構(gòu)如圖3所示,該結(jié)構(gòu)由DTCN殘差模塊堆疊和全連接層組成。DTCN的視野域取決于網(wǎng)絡深度n以及卷積核大小k和擴張因子d。DTCN殘差模塊解決了網(wǎng)絡深度發(fā)展的問題,因此DTCN通過堆疊殘差模塊,可以使模型具有更好的泛化能力。

        圖3 深度時間卷積網(wǎng)絡模型結(jié)構(gòu)

        2.3 基于改進CA的學習率熱重啟

        在DTCN的訓練中,學習率是一個非常重要的參數(shù)。學習率過大易使模型發(fā)散,過小則會增加模型訓練時間,合適的學習率調(diào)整策略有助于模型擺脫陷入局部最優(yōu)的狀態(tài)。CA優(yōu)化學習率采用循環(huán)變化的策略,在每一個訓練批次中動態(tài)調(diào)整學習率,并保存最優(yōu)模型,避免陷入局部最優(yōu)。CA基本公式為

        (19)

        式中:ηt為第t次迭代的學習率;ηmax、ηmin為學習率的最大值與最小值;Tcur為當前運行步;Ti為重啟周期。

        熱重啟的方法為模型訓練經(jīng)過Ti步之后,對模型進行重啟。重啟時改變學習率的值,并以之前的權重值為初始值重新訓練。通常Ti的值呈線性增加;但由于每次重啟時學習率的初值均相同,不利于收斂,因此需要對重啟后的初值進行改進,即

        ηmax(Ti+1)=ρηmax(Ti),

        (20)

        式中:ηmax(Ti)、ηmax(Ti+1)分別為第Ti、Ti+1階段熱重啟的最大(初始)學習率,數(shù)值范圍為ηmin<ηmax(Ti+1)<ηmax(Ti);ρ為減小因子,0<ρ<1。

        通過ρ每次重啟時減小學習率初值。這種方法消除了每個周期后具有相似模型權重的問題,在模型訓練后期易于收斂。

        3 實驗設置與結(jié)果分析

        實驗環(huán)境為Windows10操作系統(tǒng),以深度學習框架Keras為前端,TensorFlow為后端進行系統(tǒng)模型的搭建。硬件環(huán)境為Intel Core i5 7500處理器,內(nèi)存16 GB;圖形處理器(GPU)為NVIDIA?GeForce GTX 2070,顯存為8 GB。實驗程序在GPU上運行。

        3.1 實驗總體模型與參數(shù)設置

        本文中提出的風電功率預測方法的總體思路是:首先采集風電功率原始數(shù)據(jù),利用CEEMDAN提取其內(nèi)在特征,得到若干個IMF分量;計算各個IMF的排列熵值,根據(jù)熵值的計算結(jié)果對各個IMF進行重構(gòu),合并復雜度較低的IMF分量;由ICA優(yōu)化的DTCN(ICA-DTCN)對重構(gòu)完成的向量分別進行預測,將反歸一化的預測值相疊加,通過0~t時間點的數(shù)據(jù)最終得到t+1時刻,即未來1 h風電功率的預測值。實驗的總體模型如圖4所示。

        圖4 風電功率預測實驗總體模型

        該實驗利用環(huán)渤海地區(qū)某15 MW風電場2012年5月1日—7月30日的1 h風電功率數(shù)據(jù),個數(shù)為2 252,90%的數(shù)據(jù)為訓練集,其余10%為測試集,并在訓練集中隨機選擇10%的數(shù)據(jù)作為驗證集,作用是在訓練時調(diào)整參數(shù),提高模型的泛化能力。

        在輸入至DTCN預測之前,首先對重構(gòu)后的分量進行標準化處理,將發(fā)電功率映射到[-1,1],即

        (21)

        預測的結(jié)果反歸一化得

        (22)

        本實驗參數(shù)配置如下:圖2中5個殘差塊中膨脹卷積A的膨脹率均設置為1,膨脹卷積B膨脹率設置為1、2、4、8、16。構(gòu)建的損失函數(shù)為log-cosh函數(shù),優(yōu)化器的選擇為Adam。訓練時增加梯度裁剪方法,將梯度值限制在區(qū)間[-c,c],以防止梯度彌散或者爆炸。經(jīng)多次試驗,得到較優(yōu)超參數(shù),見表1。

        表1 實驗超參數(shù)設置

        3.2 模型評價指標

        (23)

        (24)

        (25)

        (26)

        3.3 實驗結(jié)果與分析

        3.3.1 CEEMDAN結(jié)果

        風電功率原始數(shù)據(jù)如圖5(a)所示,可以看出風電功率序列具有較大的隨機性和波動性。為了提取分析歷史功率特征,對原始風電功率數(shù)據(jù)進行CEEMDAN。分解過程中,添加白噪聲組個數(shù)為300,標準差為0.5,迭代終止閾值為0.001。

        通過CEEMDAN算法,原始風電功率序列共提取出為8個特征互異的模態(tài)分量和1個余項信號,分解得到的IMF分量1—8以及余項如圖5(b)—(j)所示。由圖可以看出,分解后的功率頻率依次減小,復雜度也相應降低,但各個IMF頻率和幅值范圍各不相同,需要對分解后的向量進行約簡處理,重構(gòu)輸入向量。

        圖5 原始風電功率序列及模態(tài)分解后本征模態(tài)函數(shù)(IMF)分量

        3.3.2 基于排列熵的風電功率序列重構(gòu)

        重構(gòu)時,PE重構(gòu)維數(shù)m一般為3~7,延遲因子τ通常為1。根據(jù)文獻[18]中的偽近鄰算法,本文中選擇重構(gòu)維數(shù)m=3,延遲因子τ=1,計算各個模態(tài)函數(shù)的樣本熵,結(jié)果如圖6所示。由圖可以看出,前4個模態(tài)分量的PE值相對較大,可單獨作為重構(gòu)后的向量,f5以及之后的模態(tài)函數(shù)PE值都較小,且數(shù)值接近,為避免計算復雜度,將PE值相似的f5、f6相加作為第5個分量,f7之后的模態(tài)分量相加作為重構(gòu)后的第6個分量。重組后的分量如表2所示,將重組后的特征序列分別進行預測,反歸一化后的預測結(jié)果相加作為最終結(jié)果。

        f1, f2, …, f8—IMF分量;R—IMF余項。

        表2 重構(gòu)序列組成成分

        3.3.3 ICA優(yōu)化的學習率熱重啟結(jié)果分析

        選取Ti=[4,8,16,32],ρ=0.9,經(jīng)過ICA優(yōu)化的學習率變化曲線如圖7所示。由圖可以看出,學習率在經(jīng)過步長Ti之后重啟,且重啟后學習率的初值減小,最終趨于平緩。

        圖7 改進余弦退火算法優(yōu)化的學習率

        在同等條件下比較CA優(yōu)化和ICA優(yōu)化的DTCN模型,兩者驗證集中的損失函數(shù)與迭代次數(shù)的關系如圖8所示。由圖可見,在訓練第1—8輪迭代,兩者的損失函數(shù)快速減小,且ICA-DTCN模型的損失函數(shù)值均較小,因此該模型收斂更快。

        圖8 不同學習率優(yōu)化模型的驗證集損失函數(shù)對比

        3.3.4 功率預測結(jié)果分析

        為了驗證提出的CEEMDAN-ICA-DTCN結(jié)果的一般性,選取BP、DTCN、LSTM、ICA-DTCN、EMD-DTCN、CEEMDAN-BP、CEEMDAN-LSTM和CEEMDAN-DTCN等網(wǎng)絡作為對照模型,為了保證各個模型的相對公平,模型訓練參數(shù)的設置基本相同。各個模型訓練參數(shù)如表3所示

        表3 不同模型的訓練參數(shù)個數(shù)

        CEEMDAN-ICA-DTCN模型以及其他對照模型在測試集上的風電功率預測值與真實值對比如圖9所示。由圖可以直觀地看出,在測試集上,本文中所提出的CEEMDAN-ICA-DTCN模型的預測結(jié)果可以較為準確地跟蹤風電功率觀測值。

        BP、DTCN、LSTM、ICA、EMD和CEEMDAN分別為反向傳播神經(jīng)網(wǎng)絡、深度時間卷積網(wǎng)絡、長短時記憶網(wǎng)絡、改進余弦退火、經(jīng)驗模態(tài)分解和自適應集成經(jīng)驗模態(tài)分解。

        根據(jù)模型評價指標,分別計算各個模型的相關誤差,結(jié)果如表4所示。由表可知:在單一模型中,DTCN模型的σMAPE、σRMSE、σMAE分別為24.8%、851 kW、556 kW,均為單一模型中最小,r2為0.948,為單一模型中最大。而ICA-DTCN模型的σMAPE、σRMSE、σMAE與r2的數(shù)值分別為23.7%、803 kW、518 kW、0.950,比未經(jīng)優(yōu)化的DTCN在原始數(shù)據(jù)曲線擬合程度上有所提升。在組合模型中,CEEMDAN-ICA-DTCN模型的σMAPE、σRMSE、σMAE分別為12.8%、292 kW、200 kW,均為最小,相比較于ICA-DTCN、EMD-DTCN、CEEMDAN-BP、CEEMDAN-LSTM和CEEMDAN-DTCN等模型,r2分別提升4.5%、2.3%、4.4%、1.2%和0.1%,綜合表現(xiàn)是9種模型中最好的。

        表4 實驗結(jié)果評價指標

        由表4還可以看出,單一模型的預測精度低于分解與預測的組合模型的,原因是原始風電功率序列經(jīng)過分解后,提取出其內(nèi)部特征,降低復雜度,因此組合模型的預測精度要優(yōu)于單一模型的。同時,改進余弦退火優(yōu)化的學習率也提高了訓練模型的精度。由此可知,本文中所提出的CEEMDAN-ICA-DTCN模型相較于對照模型組的風電功率預測準確度最高。針對濱海地區(qū)風速波動大的特點,該方法對近海地區(qū)的超短期風電功率預測具有較好的適應性,預測結(jié)果較為準確。

        4 結(jié)論

        針對濱海地區(qū)風電功率受離岸風影響、波動大、隨機性強等問題,本文中建立了基于CEEMDAN-ICA-DTCN的超短期風電功率預測模型。根據(jù)算例結(jié)果得出以下結(jié)論:

        1)利用CEEMDAN算法處理復雜功率序列,克服了EMD常規(guī)分解方法存在的問題,并利用排列熵相關理論,對分解后的序列進行約簡處理,重構(gòu)后的序列可以直接輸入DTCN進行預測。

        2)利用DTCN對風電功率序列進行建模,因果膨脹卷積對功率序列進行特征提取,堆疊的殘差單元可以有效地學習序列歷史特征,最終通過全連接網(wǎng)絡,對未來1 h的風電功率進行預測。

        3)經(jīng)過改進的余弦退火優(yōu)化學習率后的DTCN提高了收斂速度,在預測精度上也比未經(jīng)優(yōu)化的一般模型有所提升。

        由實驗結(jié)果對比可以看出,CEEMDAN-ICA-DTCN模型能對濱海地區(qū)的風電功率進行有效的預測,解決了短時間內(nèi)風電功率波動大、預測難等問題,為濱海地區(qū)海上風力發(fā)電并網(wǎng)規(guī)劃以及能源調(diào)度提供參考,助力我國新能源產(chǎn)業(yè)的發(fā)展。

        猜你喜歡
        電功率殘差分量
        基于PCC-CNN-GRU的短期風電功率預測
        基于雙向GRU與殘差擬合的車輛跟馳建模
        帽子的分量
        輕松上手電功率
        你會計算電功率嗎
        基于殘差學習的自適應無人機目標跟蹤算法
        一物千斤
        智族GQ(2019年9期)2019-10-28 08:16:21
        基于遞歸殘差網(wǎng)絡的圖像超分辨率重建
        自動化學報(2019年6期)2019-07-23 01:18:32
        解讀電功率
        論《哈姆雷特》中良心的分量
        一区二区三区国产黄色| 视频一区精品自拍| 资源在线观看视频一区二区| 日本一区二区三区爱爱视频| 久久人人爽av亚洲精品| 99精品热这里只有精品| 免费看欧美日韩一区二区三区| 一区二区三区手机看片日本韩国| 尤物在线观看一区蜜桃| 人人爽久久涩噜噜噜av| 亚洲最新版无码AV| 精品午夜中文字幕熟女| 夜夜爽夜夜叫夜夜高潮| 亚洲av综合久久九九| 久久国产成人午夜av影院| 中文字幕人妻激情在线视频| 国产一二三四2021精字窝| 1区2区3区高清视频| 国产真实乱XXXⅩ视频| 蜜臀av一区二区三区| 国产日韩精品欧美一区喷水| 亚洲日韩精品国产一区二区三区| 国产精品天干天干在线观蜜臀| 久久免费亚洲免费视频| 免费无码毛片一区二区app| 国产一国产一级新婚之夜| 中文字幕中文字幕人妻黑丝| 激情人妻另类人妻伦| 国产成人精品日本亚洲11| 无码免费午夜福利片在线| 深夜黄色刺激影片在线免费观看 | 全免费a级毛片免费看视频| 一片内射视频在线观看| 久久国产黄色片太色帅| 最近中文字幕视频完整版在线看 | 日韩精品在线免费视频| 国产涩涩视频在线观看| 韩国日本亚洲精品视频 | 亚洲一区二区三区精品久久| 日韩有码中文字幕在线观看| 日本一卡2卡3卡四卡精品网站|