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

        ?

        基于MEEMD-QUATRE-BILSTM的短期光伏出力區(qū)間預測

        2023-06-13 00:00:00張程林谷青匡宇
        太陽能學報 2023年11期
        關鍵詞:數(shù)據(jù)挖掘

        收稿日期:2022-07-19

        基金項目:國家自然科學基金(51677059);福建省財政廳專項(GY-Z220230);福建省自然科學基金(2023J01951)

        通信作者:張 程(1982—),男,博士、副教授,主要從事電力系統(tǒng)穩(wěn)定性分析、廣域監(jiān)測方面的研究。zhangcheng@fjut.edu.cn

        DOI:10.19912/j.0254-0096.tynxb.2022-1067 文章編號:0254-0096(2023)11-0040-15

        摘 要:提出一種基于改進集成經(jīng)驗模態(tài)分解(MEEMD)和擬仿射變換(QUATRE)優(yōu)化雙向長短期記憶神經(jīng)網(wǎng)絡(BILSTM)的光伏出力區(qū)間預測模型。通過主成分分析法(PCA)對時間序列進行降維處理,利用K-均值算法將降維數(shù)據(jù)分成3種類型氣象數(shù)據(jù);然后采用MEEMD對每類光伏出力序列進行分解,將其輸入QUATRE優(yōu)化BILSTM神經(jīng)網(wǎng)絡和核密度估計算法(KDE)聯(lián)合構建的短期光伏出力區(qū)間預測模型。最后基于寧夏光伏電站實例仿真評估模型區(qū)間預測性能,實驗結果表明該模型可生成高水平光伏預測區(qū)間,能夠為電力系統(tǒng)經(jīng)濟穩(wěn)定運行提供可靠的決策保障。

        關鍵詞:光伏發(fā)電;數(shù)據(jù)挖掘;預測;改進的集成經(jīng)驗模態(tài)分解;擬仿射變換進化算法;雙向長短期記憶網(wǎng)絡

        中圖分類號:TM615 """" 文獻標志碼:A

        0 引 言

        隨著全球常規(guī)能源局勢空前危急以及氣候變暖問題日益嚴峻,可再生清潔能源發(fā)電比例迅速提高[1]。在碳達峰、碳中和的“雙碳”背景下,中國光伏能源正在大規(guī)模接入電網(wǎng)[2]。光伏發(fā)電受天氣因素復雜性和無序性影響,光伏能源直接并入電網(wǎng)系統(tǒng)將會威脅配網(wǎng)調度運行的安全穩(wěn)定[3]。對光伏出力的準確預測能有效降低輸入電網(wǎng)能量的隨機性,提高光伏系統(tǒng)并網(wǎng)可靠性。

        光伏出力預測方法主要包括確定性點預測法[4]和基于概率學的區(qū)間預測法[5]。傳統(tǒng)確定性點預測已經(jīng)無法對光伏的不確定性因素作出有效判斷,難以滿足電網(wǎng)調度日益復雜的風險決策與評估需求,而概率區(qū)間預測能通過不同置信區(qū)間對未來光伏出力的概率信息進行分析判斷,量化不確定性影響因子對光伏系統(tǒng)的干擾,具有更重要的實際工程意義[6]。

        針對光伏時間序列的多維性與復雜性,文獻[7]利用主成分分析法(principal component analysis,PCA)對光伏數(shù)據(jù)影響因子進行篩選,降低原始時間序列的維度,消除模型輸入數(shù)據(jù)的相關性與冗余性,然后輸入長短期記憶神經(jīng)網(wǎng)絡(long short-term memory,LSTM)實現(xiàn)光伏功率預測;文獻[8]描述了一種結合集成經(jīng)驗模態(tài)分解(ensemble empirical mode decomposition,EEMD)和相關支持向量機(support vector machine,SVM)的光伏功率區(qū)間預測方法,但EEMD分解方法存在耗時過長和重構后噪聲過多等問題;為有效去除噪聲問題,文獻[9]提出基于改進的集成經(jīng)驗模式分解(modified ensemble empirical mode decomposition,MEEMD)結合最小二乘支持向量機(least squares support vector machine,LSSVM)的風功率預測模型,在分解功率序列后通過樣本熵理論合并相似模態(tài)分量,以減少模型輸入量,去除噪聲的同時,有效提升模型運行效率。

        由于光伏時間序列具有前后時序信息相關聯(lián)等特性,文獻[7]采用經(jīng)驗模態(tài)分解-主成分分析-長短期記憶神經(jīng)網(wǎng)絡(empirical mode decomposition-principal component analysis-long short-term memory,EMD-PCA-LSTM)模型對多變量時間序列實現(xiàn)動態(tài)時間建模,實驗結果證明該模型具有較高的預測精度。但LSTM在較長時序預測過程中存在極易丟失重要信息的問題,文獻[10]提出一種結合注意力機制的小波分解-雙向長短時記憶神經(jīng)網(wǎng)絡(wavelet-bidirectional long short-term memory,W-BILSTM)模型,實現(xiàn)有選擇性地選取有效信息進行發(fā)電功率預測,同時通過注意力機制的映射加權和學習參數(shù)矩陣賦予BILSTM隱含層不同的權重,仿真結果證明了該模型性能的優(yōu)越性。

        在光伏系統(tǒng)實際出力預測情形中,光伏時間序列受天氣因素隨機特性影響,易出現(xiàn)功率曲線波動性強、不同天氣類型功率出力差別大、光伏時序包含噪聲多且具有冗余性等問題。綜上所述,為了有效解決上述問題以及神經(jīng)網(wǎng)絡初始超參數(shù)具有隨機性等缺點,本文提出基于改進的集成經(jīng)驗模態(tài)分解(MEEMD)結合擬仿射變換算法(quasi-affine transformation evolutionary, QUATRE)優(yōu)化雙向長短期記憶神經(jīng)網(wǎng)絡(BILSTM)的短期光伏出力區(qū)間預測模型。首先,采用主成分分析法(priacipal component analysis,PCA)篩選光伏出力多維時間序列的關鍵影響因子,降低光伏時序維度;其次,使用K-均值算法將歷史降維序列分為晴、多云、雷雨3類數(shù)據(jù)集;利用MEEMD將每類數(shù)據(jù)集分解成不同頻率的模態(tài)分量(intrinsic mode functions,IMF),并分別將其輸入QUATRE-BILSTM神經(jīng)網(wǎng)絡進行功率預測,將各IMF預測結果疊加得到最終點預測值;最后,利用核密度估計法(kernel density estimatio,KDE)[11]得到滿足各置信水平的預測區(qū)間。在實例仿真分析中將上述模型與其他幾種模型相比較,結果表明所提模型具有更高的準確性和更強的魯棒性,驗證了該模型的有效性。

        1 研究理論

        1.1 PCA降維

        PCA通過原始數(shù)據(jù)構造一組全新的空間變量,既能夠代表原始數(shù)據(jù)的絕大多數(shù)信息又互不相關,從而實現(xiàn)數(shù)據(jù)的空間降維。

        主成分實現(xiàn)步驟如下[12]:

        1)對原始數(shù)據(jù)進行標準化處理。

        原始數(shù)據(jù)矩陣[X]為:

        [X=x11x12…x1px21x22…x2p????xn1xn2…xnp]""" (1)

        式中:[p]——指標數(shù)量;[n]——樣本數(shù)量。

        對原始數(shù)據(jù)作標準化處理:

        [x?ij=xij-xjsj] (2)

        式中:[xij]——第[j]維變量數(shù)據(jù);[xj]——第[j]維數(shù)據(jù)算術平方根;[sj]——第[j]維數(shù)據(jù)標準差。

        2)計算樣本相關系數(shù)矩陣。

        [R=r11r12…r1pr21r22…r2p????rp1rp2…rpp]" (3)

        [rij=k=1n(xki-xi)(xkj-xj)k=1n(xki-xi)2k=1n(xkj-xj)2]"""""" (4)

        式中:[xki、xkj]——[xi、xj]第[k]行元素;[xi、xj]——[xi、xj]平均值。

        3)計算相關系數(shù)矩陣[R]的特征值[(λ1,λ2,…,λp)]及相應特征向量[αi]。特征值[λp]特征方程為:

        [R-λpI=0]"""" (5)

        式中:[I]——單位向量。

        特征向量[αi]為:

        [αi=(αi1,αi2,…,αip)],[i=1,2,…,p]" (6)

        4)選擇主成分數(shù)量。由主成分分析所得[p]個主成分,計算各部分累計貢獻率[φt],選取前[t]個達到85%以上的確定為主成分數(shù)量。

        [φt=i=1tλii=1pλi, i=1,2,…,p;t=1,2,…,p]" (7)

        式中:[λi]——相關系數(shù)矩陣[R]第i個特征值;[i=1pλi]——矩陣[R]特征值之和。

        1.2 K-均值聚類法

        K-均值是最普及的聚類算法,它將數(shù)據(jù)對象分割成[K]組聚類簇。為了進一步提升聚類效果,使其在天氣因素與輸出功率等方面具有相似性。本文采用歷史輸出功率和氣象影響因素各6項統(tǒng)計指標(均值[P]、標準差[σ]、偏態(tài)系數(shù)[sk]、峰值系數(shù)[kur]、變異系數(shù)[cv、]總功率[Psum])共12個特征作為樣本特征值[13],計算式如下:

        [P=1Ni=1NPi]"""""" (8)

        [σ=1Ni=1NPi-P2] (9)

        [sk=Ni=1N(Pi-P)3(N-1)(N-2)σ]" (10)

        [kur=i=1N(Pi-P)4(N-1)σ]""" (11)

        [cv=σP] (12)

        [Psum=i=1NPi]""""" (13)

        式中:[N]——樣本點數(shù)量;[P]——平均光伏功率,W;[Pi]——各時刻光伏功率,W。

        K-均值具體方法:

        輸入:聚類數(shù)量[K],數(shù)據(jù)樣本[N=x1,x2,…,xn]。

        輸出:滿足誤差平方和最小標準的[K]個聚類。

        1)將輸入樣本[N]隨機分成[K]簇初始值,在每簇樣本中選擇一個對象作為初始聚類中心。

        2)確定每個簇的中心對象[mj],由歐氏距離計算得到每個樣本點與相對應[mj]的距離,根據(jù)中心距離最小原則再次對相應目標進行分配。

        3)重新計算每個有變化簇的聚類中心對象。

        4)循環(huán)2)、3),直至每組簇的數(shù)值不再發(fā)生改變,此時為最優(yōu)聚類結果。

        [mj=1ωjxi∈ωjxi]"""""" (14)

        [E=j=1kxi∈ωjxi-mj2]" (15)

        式中:[ωj]——第[j]組聚類樣本點數(shù)量;[ωj]——第[j]類數(shù)據(jù)樣本;[E]——數(shù)據(jù)樣本聚類誤差;[k]——設定聚類數(shù)量;[xi]——[ωj]的樣本點;[mj]——[ωj]的聚類中心對象。

        算法中,當[K]值增大至真實聚類數(shù)時,[E]值的減小速度將迅速趨于平穩(wěn)。因此,綜合評價[E]值的總體下降趨勢以及聚類指標作用等重要因素后,最終通過肘方法(Elbow)確定本文真實聚類數(shù)[K]值的大小。

        1.3 改進集成經(jīng)驗模態(tài)分解

        MEEMD通過添加白噪聲于信號波形中,再對信號作經(jīng)驗模態(tài)分解(empirical mode decomposition,EMD),減輕EMD模態(tài)混疊現(xiàn)象,有效避免噪聲等異常信號對于原始信號的干擾。光伏功率序列是一種非線性、非平穩(wěn)性信號,本文采用MEEMD [14]進行分解,具體實施過程如下:

        1)在光伏功率原始信號[st]中加入絕對值相等且方向相反的兩組白噪聲[nit]、[-nit],得到噪聲信號為:

        [s+i(t)=st+ainits-i(t)=st-ainit],[i=1, 2,…,Ne]""""" (16)

        式中:[ai]——白噪聲幅值系數(shù);[Ne]——白噪聲對數(shù)。

        2)利用EMD分別對[si+t]和[si-t]進行分解,[Iij+t]和[Iij-t]表示所得本征模態(tài)分量(IMF),分解結果表示為:

        [s+itEMDI+i1ts-i(t)EMDI-i1t]""""" (17)

        式中:[Ii1t]——第[i]個信號一階本征模態(tài)函數(shù)分量。

        3)將[I+i1t]、[I-i1t]進行平均處理,消除白噪聲殘余:

        [I1t=12Nei=1NeI+i1t+I-i1t]" (18)

        4)對[I1t]進行排列熵(PE)計算:

        [HP(m)=-g=1kPglnPg]" (19)

        式中:[m]——重構信號嵌入維數(shù);[HP(m)]——時間序列隨機程度,其值越高,隨機性越強。

        當[Pg=1/m]時,[HP(m)≈ln(m!)]達到極大值,故將排列熵[HP(m)]標準化處理:

        [HP=-g=1kPglnPg/ln(m!)]""" (20)

        此時[0≤HP≤1],當分解所得序列分量為非平穩(wěn)信號,排列熵值較大。因此,通過設置閾值[θ0]用于篩選信號熵值中非平穩(wěn)干擾信號。若[I1t]熵值大于閾值[θ0],則信號異常;反之,信號正常。經(jīng)反復實驗,本文閾值取0.5~0.6,故設置[θ0]為0.6[15]。

        5)若[I1t]異常,則返回步驟1),直至IMF分量[IPt]為正常序列信號,執(zhí)行步驟6)。

        6)去除原始光伏功率信號[st]中異常信號后,再對其余信號作EMD分解。對所得IMF分量由高頻至低頻排序,最后MEEMD分解過程簡化為:

        [stMEEMDj=1NeIjt+rt] (21)

        式中:[rt]——殘余分量。

        1.4 樣本熵

        樣本熵(sample entropy, SE)是一種改進的復雜性測試方法[16]。SE值越大,光伏功率分解所得新分量序列越復雜;SE值越小,新分量復雜程度越低,自身相似度越高。

        樣本熵可由[QSampEn(m,r,N)]表示,其表達式為:

        [QSampEn(m,r,N)=-lnBm+1(r)Bm(r)] (22)

        式中:[r]——相似閾值;[Bm(r)]——兩序列在[r]下滿足[m]個點的概率。

        樣本熵值大小與參數(shù)[m]和[r]的取值密切相關,根據(jù)研究結果反復實驗,設置參數(shù)[m]=1或2,[r=0.1~0.25Sd]([Sd]為時間序列的標準差)所得樣本熵值具備合理的理論依據(jù)。本文中,[m]=2,[r]=0.2Sd。

        2 QUATRE優(yōu)化BILSTM神經(jīng)網(wǎng)絡

        2.1 QUATRE算法基本原理

        QUATRE是文獻[17]提出的一種迭代效率高、收斂性優(yōu)異的新型算法[18]。

        本文采用QUATRE作為進化方法,在優(yōu)化過程中有效利用粒子協(xié)作,演變過程如下:

        [B=Xgbest+c?(Xr1-Xr2)X?M?X+M?B]""""" (23)

        [X=[X1,X2,…,Xps]T]"nbsp;"" (24)

        [Xgbest=[Xbesti,Xbesti,…,Xbesti]T]""""" (25)

        式中:[B]——進化指導矩陣;[Xgbest]——[ps]個粒子全局最優(yōu)坐標矩陣;[c]——差分矩陣系數(shù)因子;[Xr1-Xr2]——差分矩陣;[X]——粒子的坐標矩陣;[M]——協(xié)作搜索矩陣,其值均為0或1;[M]——[M]的關聯(lián)矩陣;[?]——矩陣元素按位相乘;[ps]——種群規(guī)模;[Xbesti]——第[i]個粒子取得最佳適應度。

        初始矩陣[Minit]為元素值為1的下三角矩陣,通過兩個連續(xù)操作可實現(xiàn)矩陣[Minit]到矩陣[M]的轉化:

        1)對矩陣[Minit]的行向量元素進行任意排列;

        2)將步驟1)矩陣的行向量進行任意排列。

        [Minit=111…11…1~11…11…11=M]"" (26)

        QUATRE通過目標函數(shù)評價個體適應度,并比較個體適應度值判斷其優(yōu)劣。相比其他算法,QUATRE克服了粒子群優(yōu)化算法收斂速度過慢等問題;同時解決了差分算法和悟空進化算法(monkey king evolution,MKE)于高維度時的偏置問題[19]。

        2.2 BILSTM神經(jīng)網(wǎng)絡

        2.2.1 LSTM

        LSTM [20]是一種深度神經(jīng)網(wǎng)絡,通過引入門控機制,使模型有選擇性地保留傳輸長期時序數(shù)據(jù)信息的功能,能夠準確有效地學習到長期依賴信息[21]。

        LSTM具有循環(huán)重復單元結構,如圖1所示,由3個門及1個核心計算節(jié)點(記憶存儲單元)構成。輸入門[dt]:控制信息輸入;輸出門[ot]:控制信息輸出;遺忘門[ft]:調節(jié)細胞信息是否更新并保持;記憶存儲單元狀態(tài)[ct]:記錄細胞實時狀態(tài)信息。遺忘門、輸入門以及輸出門聯(lián)合實現(xiàn)向單元狀態(tài)的把控,選擇性地對單元狀態(tài)增加或移除信息。

        LSTM訓練過程如下:

        1)遺忘門執(zhí)行遺忘選擇部分:

        [ft=σ(Wf[ht-1,it]+bf)]" (27)

        [σ(x)=11+exp(-x)]"""""" (28)

        式中:[ft]——遺忘門輸入;[σ]——激活函數(shù);[Wf]——遺忘門隱藏單元輸入權值矩陣;[ht-1]——上一時刻功率預測值;[it]——輸入特征向量;[bf]——遺忘門偏置矩陣。

        2)輸入門及長短期記憶神經(jīng)網(wǎng)絡部分:

        [dt=σ(Wi[ht-1,it]+bd)] (29)

        [gt=tanh(Wc[ht-1,it]+bc)]"""" (30)

        [ct=ft⊙ct-1+dt⊙gt]""""" (31)

        [tanh(x)=exp(x)-exp(-x)exp(x)+exp(-x)]""""" (32)

        式中:[dt]——輸入門輸入;[Wi]——輸入門隱藏單元權值矩陣;[bd]——輸入門偏置矩陣;[gt]——候選單元狀態(tài);[tanh(·)]——雙曲正切激活函數(shù);[Wc]——記憶單元權值矩陣;[bc]——記憶單元偏置矩陣;[ct]——記憶單元輸出;[ct-1]——上一時刻狀態(tài)(長期記憶狀態(tài));[⊙]——向量中元素按位相乘。

        3)輸出門確定粒子狀態(tài)并傳送至此刻隱含層狀態(tài)中:

        [ot=σ(Wo[ht-1,it]+bo)]"""""" (33)

        [ht=ot⊙tanh(ct)]""" (34)

        式中:[ot]——輸出門狀態(tài);[Wo]——輸出門隱藏單元權值矩陣;[bo]——輸出門偏置矩陣;[ht]——單元輸出光伏功率預測值。

        2.2.2 BILSTM

        由于LSTM為單向提取序列信息,而光伏功率序列具有信息前后關聯(lián)、隨機性強等特性。為了能夠合理利用氣象因素和歷史光伏功率數(shù)據(jù),提升光伏發(fā)電功率預測能力,本文選用BILSTM構建預測模型。

        BILSTM網(wǎng)絡[22]是由正向和反向LSTM神經(jīng)網(wǎng)絡組成,可對時間序列實現(xiàn)正反向兩次LSTM訓練,有效提升特征選取的全面性和完整性。BILSTM結構[23]如圖2所示。

        圖2中前向LSTM層輸出[ht]與后向LSTM層輸出[ht]相連接,經(jīng)過加權融合得到最終功率輸出值[Ot]。BILSTM更新方程為:

        [ht=LLSTM(ht-1,it)] [t=1,2,…,n]"" (35)

        [ht=LLSTM(ht+1,it)] [t=1,2,…,n]" (36)

        [Ot=f(Whht+Whht-1+bt)]"""" (37)

        式中:[it]——輸入特征向量;[ht]、[ht]——前向和后向功率預測值;[LLSTM(·)]、[LLSTM(·)]——LSTM網(wǎng)絡雙向計算過程;[Wh]、[Wh]——雙向輸出連接權值矩陣;[bt]——輸出層偏置;[Ot]——網(wǎng)絡最終輸出功率預測值。

        3 基于MEEMD-QUATRE-BILSTM預測模型的建立

        3.1 預測步驟

        由于光伏輸出功率具備強烈的波動性與隨機性,傳統(tǒng)光伏功率預測流程未對歷史數(shù)據(jù)進行降噪處理,大量噪聲可能導致光伏出力無法進行較準確地預測分析。因此,本文提出基于MEEMD-QUATRE-BILSTM的組合預測模型,該模型整體預測流程如圖3所示,具體步驟如下:

        1)初始化BILSTM網(wǎng)絡結構,對訓練集和測試集進行歸一化處理。

        2)將數(shù)據(jù)進行PCA降維,篩選出光伏功率出力關鍵影響因子,消除數(shù)據(jù)不同時間序列的冗余性與相關性。

        3)結合歷史功率數(shù)據(jù)和6個影響因子統(tǒng)計指標,利用K-均值算法將歷史光伏出力數(shù)據(jù)聚類,按照天氣類型劃分為晴、多云及雷雨天氣。

        4)通過MEEMD將3種天氣類型原始功率序列分解為不同頻率IMF分量和殘余分量(residual error,RE),合并樣本熵值相近的子序列。

        5)初始化父代種群參數(shù),包括種群個體數(shù)目[ps、]搜索空間[[Xmin,Xmax]]、空間維度[D]、差分矩陣系數(shù)因子[c]、最大迭代次數(shù)[Jmax]、種群位置[X=[X1,X2,…,Xps]T],設當前位置為種群個體歷史最優(yōu)解[Xpbest]和全局最優(yōu)解[Xgbest]。

        6)根據(jù)式(5)進化父代種群,產生子代種群并更新坐標位置[X]。

        7)確定粒子目標函數(shù)值[fitness(Xi)]。

        8)將子代目標函數(shù)值與歷史最優(yōu)函數(shù)值以及全局最優(yōu)函數(shù)值相比較,若[fitness(Xi)gt;fitness(Xpbest)],則局部最優(yōu)位置[Xpbest]由當前位置替代,并更新[fitness(Xpbest)];若[fitness(Xi)gt;fitness(Xgbest)],則更新全局最優(yōu)位置[Xgbest]其函數(shù)值[fitness(Xgbest)]。

        9)判斷迭代是否終止,若迭代達到最大次數(shù)或滿足收斂精度,則迭代終止,輸出最終優(yōu)化結果;否則回轉至步驟2),繼續(xù)搜尋最優(yōu)結果。

        10)將QUATRE全局最優(yōu)位置各維數(shù)值設為BILSTM網(wǎng)絡最優(yōu)超參數(shù)組合。

        11)分別將重構后IMF分量和殘差分量輸入QUATRE- BILSTM模型中,得到各分量預測結果,將其合并形成功率點預測結果。

        12)通過KDE對點預測絕對誤差直方圖近似擬合出對應的概率密度函數(shù)曲線最終得到預測區(qū)間。MEEMD-QUATRE-BILSTM光伏功率模型建模完成。

        3.2 模型參數(shù)設定

        經(jīng)多次實驗調優(yōu),BILSTM超參數(shù)[24]確定如下:輸入層維數(shù)為3;輸入層時間步數(shù)為5;隱含層數(shù)量為2;每層隱含層單元數(shù)均為50;輸出層維數(shù)為1;初始學習率為0.001。

        QUATRE參數(shù)設定:子代種群數(shù)量[ps]為100;搜索最小范圍[Xmin]取[(0,0,0,0,0)T];最大范圍[Xmax]取[(1,0.1,1,1,5)T];空間維度D為5;差分矩陣系數(shù)因子[c]取0.7;最大迭代次數(shù)[Jmax]取500。

        3.3 模型評價指標

        本文模型精度評價分析采用平均絕對誤差(mean absolute error,MAE)、均方誤差(mean squared error, MSE)、確定系數(shù)(coefficient of determination,R2)和均方根誤差(root mean square error,RMSE);區(qū)間預測評價指標采用KL散度(Kullback-Leibler divergence, KLD)、誤差平方和(sum of squares due to error,SSE)、預測區(qū)間覆蓋率(prediction interval coverage probability, PICP)、預測區(qū)間平均寬度(prediction interval normalized average,PINAW),以證明所提方法的有效性。以上指標計算公式[25]不再展開介紹。

        4 算例分析

        4.1 數(shù)據(jù)來源

        為驗證MEEMD分解的有效性以及PCA降維、K-均值聚類的必要性,證明基于MEEMD-QUATRE-BILSTM光伏出力區(qū)間預測組合模型的優(yōu)越性。本文選取寧夏某光伏發(fā)電站2021年5月1日—10月7日,每天07:00—19:00采集間隔15 min的7680個數(shù)據(jù)采樣點作為實驗樣本,僅分析白天功率有效時段。

        將樣本數(shù)據(jù)集中的7680個時間數(shù)據(jù)轉化為適于BILSTM網(wǎng)絡數(shù)據(jù)集,經(jīng)K-均值聚類后分別將每種天氣類型數(shù)據(jù)的前80%設為訓練集,后20%設為測試集。

        4.2 數(shù)據(jù)處理

        4.2.1 數(shù)據(jù)降維

        考慮光伏出力影響因素眾多,各因素間存在數(shù)據(jù)的冗余性,導致預測模型輸入維度過于龐大,計算耗時時間長。本文通過PCA降維,篩選出對原始特征序列有強代表性的綜合變量,實驗結果如圖4及表1所示。

        由圖4可見,當主成分數(shù)量大于3時,其特征值減小緩慢,出現(xiàn)明顯拐點,說明此刻主成分所對應貢獻率不再大幅增加,故主成分數(shù)量確定為3。由表1可知,前3部分的累計貢獻率達到91.962%,所包含原始信息量超過85%,因此選擇前3個特征序列代替原變量作為預測模型輸入。

        4.2.2 數(shù)據(jù)聚類選取K值

        在K-均值中,通常依靠肘方法確定最佳聚類數(shù)[K]。肘方法是對[K]值逐一進行聚類并記錄對應的誤差平方和([VSSE]),以此畫出[K]值與[VSSE]的關系圖,如圖5所示。

        根據(jù)手肘圖的拐點確定最佳[K]的取值,由圖5確定當[K=3]時聚類效果最佳,不同天氣類型功率曲線聚類結果如圖6所示,將天氣類型劃分為晴、多云、雷雨天氣。同種天氣類型下光伏出力曲線相似,由圖6a可看出,晴天功率曲線總體較平穩(wěn);由圖6b可看出,多云天氣由于云層遮擋陽光,曲線出現(xiàn)間歇性波動;由圖6c可看出,雷雨天氣曲線呈現(xiàn)不規(guī)則強烈波動,其相應功率值也較低。

        4.2.3 數(shù)據(jù)分解

        光伏實驗樣本數(shù)據(jù)為非平穩(wěn)時序信號,受天氣變化影響,具有一定隨機性和突變性。本文應用MEEMD對光伏出力原始時間序列進行分解處理,為了驗證MEEMD方法的有效性與合理性,圖7為3種天氣類型下MEEMD分解得到不同尺度的IMF模態(tài)分量和RS殘余分量,突出顯示了原始光伏出力序列局部特性。

        由圖7可知,MEEMD將晴、多云、雷雨天氣原始光伏出力時間序列分別分解為8個IMF分量與1個殘余分量RS。其中,晴與多云天氣原始數(shù)據(jù)出力平緩,僅有IMF1、IMF2頻率變化較大;而雷雨天氣原始數(shù)據(jù)由于出力差別大,IMF1、IMF2、IMF3頻率波動變化異常強烈,每組IMF均一定程度上反映光伏功率的局部特征。由此可見,光伏數(shù)據(jù)受天氣隨機性影響越大、功率波動越強,異常數(shù)據(jù)分量含量越多,相應IMF頻率波動越厲害。MEEMD將原始強波動性數(shù)據(jù)分解篩選為波動性相對較弱的數(shù)據(jù),降低待預測光伏數(shù)據(jù)的復雜程度,從而提升了預測精準度。

        由于IMF分量過多,若針對每個IMF分量逐一輸入QUATRE-BILSTM模型進行預測分析,將會擴大計算規(guī)模。本文應用復雜系統(tǒng)理論的樣本熵法,對分解后的9個分量進行復雜程度分析,減小光伏出力預測計算量,3種天氣類型下各IMF分量樣本熵值分析結果如圖8所示。

        圖8 各IMF分量樣本熵值

        由圖8可知,3種類型天氣樣本熵值整體逐級下降,復雜程度隨IMF頻率降低而逐漸減少,符合樣本熵理論自身特性,證明該理論適用于短期光伏出力預測。

        經(jīng)過MEEMD分解后,將樣本熵值相鄰部分的IMF分量進行疊加處理,達到減少模態(tài)分量的目的。合并重構后的新子序列結果如圖9所示,圖9a和圖9b合并為5個IMF分量與1個殘余分量RS;圖9c雷雨天氣由于復雜天氣因素影響,重構為6個IMF分量與1個殘余分量RS。

        4.3 模型適用性評估

        4.3.1 點預測分析

        為驗證預測模型的可行性,選取不同天氣類型下各模型預測效果進行對比分析。利用寧夏歷史光伏數(shù)據(jù)集,分別通過卷積神經(jīng)優(yōu)化BILSTM網(wǎng)絡組合模型(CNN-BILSTM)、麻雀算法優(yōu)化BILSTM神經(jīng)網(wǎng)絡組合模型(SSA-BILSTM)、擬仿射算法優(yōu)化LSTM神經(jīng)網(wǎng)絡組合模型(QUATRE-LSTM)、單一

        BILSTM神經(jīng)網(wǎng)絡模型與所提QUATRE-BILSTM預測模型對3種典型天氣(晴天、多云、雷雨天氣)光伏功率出力情形進行預測對比,結果如圖10、圖11及表2所示。

        由圖10a可知,在晴天天氣情況下,光伏功率出力總體趨勢平穩(wěn)。正午輻照度在13:00左右達到峰值,各模型在峰值附近預測效果存在一定差距,MEEMD-QUATRE-BILSTM模型追峰效果最佳。結合圖11a可知,晴天預測功率絕對誤差均在3 kW以內,低于其他預測模型。

        由圖10b可知,對多云天氣而言,由于云層對陽光不同程度遮擋,各模型預測精度整體出現(xiàn)偏差。單一BILSTM神經(jīng)網(wǎng)絡模型因未對噪聲及超參數(shù)進行優(yōu)化處理,導致其波動時刻反應遲鈍,與實際功率曲線存在較大差異;MEEMD-CNN-BILSTM與MEEMD-SSA-BILSTM模型雖然對噪聲及超參數(shù)作優(yōu)化處理,但其優(yōu)化效果顯然不如QUATRE算法顯著。結合圖11b可知,雖然多云天氣實際功率出現(xiàn)多處波動,但所提預測模型絕對誤差仍在5 kW以內,擁有較高的預測精度。

        由圖10c可知,在雷雨天氣情況下,光伏出力曲線由于日照強度不足以及復雜氣象因素引起的不確定性變化,使得實際輸出功率曲線具有頻繁無規(guī)律的波動,故各模型大部分時刻預測值與真實值存在一定幅度偏離。結合圖11c可得,MEEMD-QUATRE-BILSTM預測誤差較其余對比模型更小,即便在極端天氣情況下絕對誤差仍能控制在8 kW以內,可以證明該模型在光伏功率預測中更具適用性。

        表2為不同天氣類型下模型預測評價指標,相比單一BILSTM神經(jīng)網(wǎng)絡模型,經(jīng)降噪、超參數(shù)優(yōu)化處理后的MEEMD-QUATRE-BILSTM模型,[TMAE]提升65.2%~80.9%、[TMAE]提升64.5%~79.1%、[TR-square]提升70.9%~3.1%、[TRMSE]提升64.6%~78.8%;MEEMD-QUATRE-BILSTM模型相比MEEMD-QUATRE-LSTM模型,[TMAE]提升65.0%~72.3%、[TMAE]提升59.7%~71.6%、[TR-square]提升2.1%~44.6%、[TRMSE]提升59.9%~71.1%,表明BILSTM相比LSTM加強了數(shù)據(jù)有效特征時序的前后關聯(lián),提高預測精準度;MEEMD-QUATRE-BILSTM模型相比MEEMD-CNN-BILSTM和MEEMD-SSA-BILSTM模型,[TMAE]進一步提升42.4%~63.6%、[TMAE]進一步提升37.9%~59.7%、[TR-square]進一步提升1.0%~38.2%、[TRMSE]進一步提升37.9%~59.7%,實驗結果表明QUATRE優(yōu)化算法的優(yōu)越性,可以明顯改善模型的結構參數(shù)。

        4.3.2 區(qū)間預測分析

        預測誤差分布特性是建立光伏出力預測區(qū)間的基礎,誤差分布擬合情況的準確性直接決定區(qū)間預測的精準程度。

        為研究非參數(shù)KDE方法對于擬合預測誤差分布的優(yōu)越性,采用正態(tài)分布、邏輯分布、t Location-Scale分布等與非參數(shù)KDE模型分別對3種類型天氣誤差樣本進行擬合對比,通過擬合效果比較分析各方法的準確性,用于確定預測區(qū)間上下限。

        針對3種不同天氣類型預測誤差分布擬合情況進行對比分析,各誤差模型概率密度分布如圖12所示。MEEMD-QUATRE-BILSTM預測模型所得誤差分布具有偏度特性,出現(xiàn)多峰值情況,圖12a晴天預測誤差分布中出現(xiàn)明顯雙峰值;在圖12b多云預測誤差分布中出現(xiàn)明顯三峰值;而在圖12c雷雨天預測誤差分布中出現(xiàn)明顯四峰值,因此正態(tài)分布、邏輯分布、t Location-Scale分布、極值分布等參數(shù)估計方法無法對此類預測誤差分布形成多峰擬合曲線,只有非參數(shù)KDE擬合曲線滿足預測誤差分布的多峰值特性,其余分布方案皆無法對實際預測誤差分布實現(xiàn)有效擬合。圖13中3種天氣類型誤差累積概率分布結果表明,非參數(shù)KDE累積概率分布與實際累積概率分布曲線變化趨勢最為接近,擬合精度較高;其余累積概率分布則與實際累積概率分布偏差較大,擬合效果不太理想。結合表3可知,在3種天氣類型下,非參數(shù)KDE分布的[TKLD]和[TSSE]擬合評價指標均優(yōu)于其余分布方案。通過擬合效果對比分析可知,非參數(shù)KDE具有適應性靈活、擬合預測誤差準確度高等優(yōu)點。

        為驗證本文MEEMD-QUATRE-BILSTM模型在不同置信水平下的區(qū)間預測效果,選取80%、85%、90%、95%這4個置信水平為例。預測日選取與點預測數(shù)據(jù)一致,結果如圖14所示。由圖14可知,本文模型實際值絕大部分均落于置信水平為95%的置信區(qū)間以內,少部分落在80%置信水平預測區(qū)間之外。在不同天氣類型下,區(qū)間預測的平均寬度相差較大。由圖14a和圖14b可知,天氣因素較為穩(wěn)定,所提方法預測區(qū)間寬度較窄;由圖14c可知,雷雨天受復雜氣象因素影響,不確定性因素增多,使得預測區(qū)間寬度明顯增大,預測效果不佳。但預測區(qū)間總體效果較好,其區(qū)間上下限、功率變化趨勢均能有效包裹實際功率曲線,區(qū)間覆蓋率基本符合預先設定的置信度,滿足實際情況需求,體現(xiàn)了本文模型的有效性。

        評估指標如表4所示,光伏功率區(qū)間預測平均寬度隨置信度的增加而增加,雷雨天受復雜氣象因素影響,平均寬度和覆蓋率指標相比晴天和多云天氣均較差。上述結果說明,非參數(shù)KDE能較好地擬合出預測誤差分布,提供穩(wěn)定的上下限偏差值。同時,MEEMD-QUATRE-BILSTM模型能夠保證在追蹤光伏出力時序波動較劇烈情形的同時,具有較窄的區(qū)間上下限。

        總體而言,當光伏功率出力波動性強、受不確定氣象因素影響時,預測區(qū)間范圍也會發(fā)生波動,但區(qū)間預測結果基本符合置信水平,算例結果驗證了本文模型的適用性和優(yōu)越性。

        5 結 論

        針對光伏出力時間序列的間歇性與不穩(wěn)定性,本文提出基于MEEMD-QUATRE-BILSTM的短期光伏出力區(qū)間預測組合模型。通過寧夏光伏電站實例仿真評估該模型,得到如下主要結論:

        1)利用PCA對影響光伏出力的氣象因子進行降維,提取包含原始信息85%以上的代表變量,在保證精度的情況下簡化模型輸入向量維數(shù)。

        2)使用K-均值算法將歷史數(shù)據(jù)集分為晴、多云、雷雨天3類,相同天氣類型可減輕光伏功率數(shù)據(jù)的波動性,提高模型預測精度。

        3)采用MEEMD有效解決了EMD噪聲處理不當?shù)葐栴};樣本熵通過對相似分量進行疊加,有效提高了模型運行效率。

        4)BILSTM網(wǎng)絡超參數(shù)的選取對模型預測精度至關重要,利用QUATRE迭代速度快、控制參數(shù)少等特點對BILSTM超參數(shù)進行快速尋優(yōu),通過實驗對比證明其優(yōu)化效果的優(yōu)越性。

        5)經(jīng)過數(shù)據(jù)集降維聚類去噪處理,提高了QUATRE-BILSTM模型點預測精度;通過非參數(shù)KDE擬合光伏預測誤差分布情況,得到高質量區(qū)間誤差上下限,從而提高區(qū)間預測相關評估指標值?;诠夥娬静煌鞖忸愋?,在80%~95%的置信區(qū)間中,預測區(qū)間覆蓋率基本與設定置信水平相符,滿足實際情況對于預測結果的要求,充分驗證了所提模型的有效性和廣泛適用性,能夠為光伏并網(wǎng)系統(tǒng)的經(jīng)濟穩(wěn)定運行提供可靠且具備參考意義的預測信息。

        [參考文獻]

        [1]"""" PILLOT B, MUSELLI M, POGGI P, et al. Historical trends in global energy policy and renewable power system issues in Sub-Saharan Africa:the case of solar PV[J]. Energy policy,2019, 127: 113-124.

        [2]"""" 李升, 衛(wèi)志農, 孫國強, 等. 大規(guī)模光伏發(fā)電并網(wǎng)系統(tǒng)電壓穩(wěn)定分岔研究[J]. 電力自動化設備, 2016, 36(1): 17-23.

        LI S, WEI Z N, SUN G Q, et al. Voltage stability bifurcation of large-scale grid-connected PV system[J]. Electric power automation equipment, 2016, 36(1): 17-23.

        [3]"""" 崔楊, 張匯泉, 仲悟之, 等. 考慮需求響應的含光熱電站可再生能源高滲透率電力系統(tǒng)多源優(yōu)化調度[J]. 高電壓技術, 2020,46(5): 1486-1496.

        CUI Y, ZHANG H Q, ZHONG W Z, et al. Multi-source optimal scheduling of renewable energy high-permeability power system with CSP plants considering demand response[J]. High voltage engineering, 2020, 46(5): 1486-1496.

        [4]"""" 賴昌偉, 黎靜華, 陳博, 等. 光伏發(fā)電出力預測技術研究綜述[J]. 電工技術學報, 2019, 34(6): 1201-1217.

        LAI C W, LI J H, CHEN B, et al. Review of photovoltaic power output prediction technology[J]. Transactions of China Electrotechnical Society, 2019, 34(6): 1201-1217.

        [5]"""" 楊錫運, 張璜, 關文淵, 等. 基于滑動分塊百分位數(shù)Bootstrap法的風電功率概率區(qū)間預測[J]. 太陽能學報, 2019, 40(2): 430-437.

        YANG X Y, ZHANG H, GUAN W Y, et al. Probabilistic intervals forecasting of wind power based on moving block percentile bootstrap method[J]. Acta energiae solaris sinica, 2019, 40(2): 430-437.

        [6]"""" 李燕青, 杜瑩瑩. 基于雙維度順序填補框架與改進Kohonen天氣聚類的光伏發(fā)電短期預測[J]. 電力自動化設備, 2019, 39(1): 60-65.

        LI Y Q, DU Y Y. Short-term photovoltaic power forecasting"" based"" on"" double-dimensional"" sequential imputation framework and improved Kohonen clustering[J]. Electric power automation equipment, 2019, 39(1): 60-65.

        [7]"""" 張雲(yún)欽, 程起澤, 蔣文杰, 等. 基于EMD-PCA-LSTM的光伏功率預測模型[J]. 太陽能學報, 2021, 42(9): 62-69.

        ZHANG Y Q, CHENG Q Z, JIANG W J, et al. Photovoltaic power prediction model based on EMD-PCA-LSTM[J]. Acta energiae solaris sinica, 2021, 42(9): 62-69.

        [8]"""" BEHERA M K, NAYAK N. A comparative study on short-term PV power forecasting using decomposition based optimized""" extreme""" learning""" machine""" algorithm[J]. Engineering science and technology, 2020, 23(1): 156-167.

        [9]"""" 魏樂, 李思瑩. 基于MEEMD-LSSVM的風電功率超短期預測研究[J]. 智慧電力, 2020, 48(5): 21-26.

        WEI L, LI S Y. Ultra-short-term forecast for wind power based on MEEMD-LSSVM[J]. Smart power, 2020, 48(5): 21-26.

        [10]""" 謝小瑜, 周俊煌, 張勇軍, 等. 基于W-BiLSTM的可再生能源超短期發(fā)電功率預測方法[J]. 電力系統(tǒng)自動化, 2021, 45(8): 175-184.

        XIE X Y, ZHOU J H, ZHANG Y J, et al. W-BiLSTM based ultra-short-term generation power prediction method of renewable energy[J]. Automation of electric power systems, 2021, 45(8): 175-184.

        [11]""" 楊茂, 張強. 風電功率超短期預測誤差的非參數(shù)估計分布研究[J]. 東北電力大學學報, 2018, 38(1): 15-20.

        YANG M, ZHANG Q. The research of ultra short-term wind power prediction error distribution based on nonparametric estimation[J]. Journal of Northeast Electric Power University, 2018, 38(1): 15-20.

        [12]""" BIANCHI F M, DE SANTIS E, RIZZI A, et al. Short-term electric load forecasting using echo state networks and PCA decomposition[J]. IEEE access, 2015, 3: 1931-1943.

        [13]""" LIN P J, PENG Z N, LAI Y F, et al. Short-term power prediction for photovoltaic power plants using a hybrid improved Kmeans-GRA-Elman model based on multivariate meteorological factors and historical power datasets[J]. Energy conversion and management, 2018, 177: 704-717.

        [14]""" GU X H, YANG S P, LIU Y Q, et al. Rolling element bearing faults diagnosis based on kurtogram and frequency domain correlated kurtosis[J]. Measurement science and technology, 2016, 27(12): 125019.

        [15]""" 鄭近德, 程軍圣, 楊宇. 改進的EEMD算法及其應用研究[J]. 振動與沖擊, 2013, 32(21): 21-26, 46.

        ZHENG J D, CHENG J S, YANG Y. Modified EEMD algorithm and its applications[J]. Journal of vibration and shock, 2013, 32(21): 21-26, 46.

        [16]""" CHEN J Y, ZHOU D, LYU C, et al. An integrated method based on CEEMD-SampEn and the correlation analysis algorithm for the fault diagnosis of a gearbox under different working conditions[J]. Mechanical systems and signal processing, 2018, 113: 102-111.

        [17]""" MENG Z Y, PAN J S, XU H R. QUasi-affine TRansformation Evolutionary (QUATRE) algorithm: a cooperative swarm based algorithm for global optimization[J]. Knowledge-based systems, 2016, 109: 104-121.

        [18]""" MENG Z Y, PAN J S. QUasi-Affine TRansformation Evolution with External ARchive (QUATRE-EAR): an enhanced"""" structure"""" for"""" Differential"""" Evolution[J]. Knowledge-based systems, 2018, 155: 35-53.

        [19]""" MENG Z Y, PAN J S. Monkey King Evolution: a new memetic evolutionary algorithm and its application in vehicle"" fuel"" consumption"" optimization[J]." Knowledge-based systems, 2016, 97: 144-157.

        [20]""" SONG X Y, LIU Y T, XUE L, et al. Time-series well performance prediction based on Long Short-Term Memory (LSTM) neural network model[J]. Journal of petroleum science and engineering, 2020, 186: 106682.

        [21]""" KARIM F, MAJUMDAR S, DARABI H, et al. LSTM fully convolutional networks for time series classification[J]. IEEE access, 2017, 6: 1662-1669.

        [22]""" XIA T B, SONG Y, ZHENG Y, et al. An ensemble framework based on convolutional bi-directional LSTM with multiple time windows for remaining useful life estimation[J]." Computers" in" industry," 2020," 115: 103182.

        [23]""" WANG S X, WANG X, WANG S M, et al. Bi-directional long short-term memory method based on attention mechanism and rolling update for short-term load forecasting[J]. International journal of electrical power amp; energy systems, 2019, 109: 470-479.

        [24]""" LI W D, NG W W Y, WANG T, et al. HELP: an LSTM-based approach to hyperparameter exploration in neural network learning[J]. Neurocomputing, 2021, 442: 161-172.

        [25]""" 龔鶯飛, 魯宗相, 喬穎, 等. 光伏功率預測技術[J]. 電力系統(tǒng)自動化, 2016, 40(4): 140-151.

        GONG Y F, LU Z X, QIAO Y, et al. An overview of photovoltaic energy system output forecasting technology[J]. Automation of electric power systems, 2016, 40(4): 140-151.

        SHORT-TERM PV OUTPUT INTERVAL PREDICTION BASED ON

        MEEMD-QUATRE-BILSTM

        Zhang Cheng1,2,Lin Guqing1,Kuang Yu1

        (1. School of Electronic Electrical and Physics, Fujian University of Technology, Fuzhou 350118, China;

        2. Fujian Provincial University Engineering Research Center for Simulation Analysis and Integrated Control of Smart Grid,

        Fujian University of Technology, Fuzhou 350118, China)

        Abstract:In this paper, a prediction model of PV output interval based on improved ensemble empirical mode decomposition (MEEMD) and quasi affine transformation (QUATRE) optimized BILSTM was proposed. Principal component analysis(PCA) was used to reduce dimension of time series, and K-means clustering algorithm was used to divide dimension reduction data into three types of meteorological data. Then MEEMD was used to decompose the output sequence of each type of PV and input it into the short-term PV output interval prediction model jointly constructed by QUATRE optimized BILSTM neural network and kernel density estimation(KDE) algorithm. Finally, the interval prediction performance of the model is evaluated based on the example of Ningxia photovoltaic power station. The experimental results show that the model can generate a high level of photovoltaic prediction interval and provide reliable decision-making guarantee for the economic and stable operation of the power system.

        Keywords:photovoltaic power; data mining; forecasting; modified ensemble empirical mode decomposition; quasi-affine transformation evolutionary algorithm; bi-directional long-short-term memory (BILSTM) network

        猜你喜歡
        數(shù)據(jù)挖掘
        基于數(shù)據(jù)挖掘的船舶通信網(wǎng)絡流量異常識別方法
        探討人工智能與數(shù)據(jù)挖掘發(fā)展趨勢
        數(shù)據(jù)挖掘技術在打擊倒賣OBU逃費中的應用淺析
        基于并行計算的大數(shù)據(jù)挖掘在電網(wǎng)中的應用
        電力與能源(2017年6期)2017-05-14 06:19:37
        數(shù)據(jù)挖掘技術在中醫(yī)診療數(shù)據(jù)分析中的應用
        一種基于Hadoop的大數(shù)據(jù)挖掘云服務及應用
        數(shù)據(jù)挖掘在高校圖書館中的應用
        數(shù)據(jù)挖掘的分析與探索
        河南科技(2014年23期)2014-02-27 14:18:43
        基于GPGPU的離散數(shù)據(jù)挖掘研究
        利用數(shù)據(jù)挖掘技術實現(xiàn)LIS數(shù)據(jù)共享的開發(fā)實踐
        欧美在线不卡视频| 亚洲午夜成人精品无码色欲| 在线 | 一区二区三区四区| 国产精品嫩草影院AV| 亚洲av中文无码乱人伦在线咪咕| 国产又黄又湿又爽的免费视频| 新婚人妻不戴套国产精品| 欧美国产精品久久久乱码| 美女裸体无遮挡免费视频的网站| 99久久免费中文字幕精品| 国内自拍色第一页第二页| 精品丰满人妻无套内射| 日日人人爽人人爽人人片av| 亚洲人成伊人成综合网中文| 手机av在线中文字幕| 久久久久久人妻一区精品| 青青在线精品2022国产| 久久久成人av毛片免费观看| 中文字幕一区二区三区视频| 小宝极品内射国产在线| 中文 国产 无码免费| 狠狠久久av一区二区三区| 尤物yw午夜国产精品视频| 国产精品igao视频| 日本精品免费一区二区三区| 国产自拍一区二区三区| 久久国产成人精品国产成人亚洲 | 天天干夜夜操| 美女极度色诱视频国产免费| 成人自拍三级在线观看| 国自产拍偷拍精品啪啪一区二区| 少妇被粗大的猛进69视频| 91精品国产综合久久青草| 国产一区av男人天堂| 真人新婚之夜破苞第一次视频| 99精品热6080yy久久| 国产精品美女主播在线| 国产大片黄在线观看| 无码人妻一区二区三区在线视频 | 国产精品自产拍在线观看免费| 久久久免费精品国产色夜|