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

基于多尺度并行卡爾曼濾波算法的電池狀態參數估算

2022-04-12 01:42:51孔德昊劉勝永
廣西科技大學學報 2022年2期
關鍵詞:卡爾曼濾波模型

孔德昊 劉勝永

摘? ?要:準確估算并合理利用電池的荷電狀態(state of charge,SOC)與健康狀態(state of health,SOH)可以延長電池的使用壽命。為了實現準確的SOC-SOH在線估計,在擴展卡爾曼濾波的基礎上,采用多尺度并行擴展卡爾曼濾波估計算法(multi-scale double extended Kalman filter,MDEKF)提高估計精度。在建立電池2階RC等效電路模型上,利用最小二乘法對模型參數進行辨識,設計并行結構的濾波器進行電池SOC估計和參數修正,并以電池組容量值作為表征量對SOH進行估算。仿真實驗結果表明,SOC估計誤差由1.43%降低到1.10%,SOH估計結果穩定在0.5%以內,驗證了算法的快速收斂性和實時性。

關鍵詞:SOC-SOH估計;多尺度并行擴展卡爾曼濾波;2階RC等效模型

中圖分類號:TM912? ? ? ? ?DOI:10.16375/j.cnki.cn45-1395/t.2022.02.008

0? ? 引言

電池參數SOC和SOH是描述電池狀態的兩個關鍵參數,二者之間存在緊密聯系。在對SOC進行單獨估計時,若電池SOH已經產生偏移,會使得SOC估計產生較大誤差。所以,對二者進行準確的估計可以保證電池能夠合理的充放電,從而延長電池的壽命[1-2]。

國內外早期對SOC的估計方法有開路電壓法和安時積分法,前者需要測量靜止電壓,缺乏實時性;后者容易在積分過程中累計誤差,造成結果發散。文獻[3]用模糊規則優化算法改進神經網絡結構的方法,完成了對SOC的估計。但神經網絡算法離不開大量實驗數據的支持,并且訓練方法也會對運算速度和結果精度產生較大的影響,所以通過建立電池等效模型對SOC進行估算逐漸成為研究主流。文獻[4]通過建立SOC與開路電壓、電池模型參數的關系,提出自適應無跡卡爾曼濾波算法。文獻[5]引用噪聲自適應因子,利用平方根球形無跡卡爾曼濾波算法完成對SOC的估算。文獻[6]提出了自適應無跡卡爾曼濾波算法和H無窮濾波算法相結合的方法對電池狀態和參數聯合估計,實驗結果表明該方法有更好的魯棒性。文獻[7]采用自適應無跡卡爾曼濾波算法對列車蓄電池在實際工況下的SOC-SOH進行估算。但是,以上研究均忽略了當SOC初始值存在誤差時對估算結果的影響。文獻[8]證明了雙擴展卡爾曼濾波算法估計SOC的可行性,在SOC初值存在誤差的情況下,估計結果也可較好收斂于估計值,但論文未對SOH的估算進行討論。

針對以上問題,本文提出采用并行擴展卡爾曼算法進行SOC-SOH聯合估算,通過MDEKF算法對電池組SOC和模型參數同時進行估計,將電池組容量作為表征參數,通過實時的容量參數值求得相對應的SOH,在UDDS工況下驗證SOC-SOH的在線估計,達到了預期的誤差精度要求,解決了卡爾曼濾波無法同時估算SOC、SOH的問題,并提高了估算精度。

1? ? 電池模型

1.1? ?等效電路模型建立

文獻[9]表明,等效電路模型能夠準確地反映電池的物理變化和化學變化,且計算量低,可以成功應用于電池SOC與SOH的估算。故本文采用2階RC電路模型,如圖1所示。

圖1中:[Uocv、Uout]為開路電壓和端電壓,[I]為電池工作電流,[R0]為歐姆內阻,[C1、R1、U1]為電化學極化電容、電阻和電壓,[C2、R2、U2]為濃度極化電容、電阻和電壓。

根據基爾霍夫定律得到系統方程和觀測方程。定義放電電流方向為正,充電電流方向為負。

系統方程:

[dU1dt=-1R1C1U1+1C1I,dU2dt=-1R2C2U2+1C2I,dSSOCdt=-1QnI.]? ? ? ? ? ? (1)

觀測方程:

[Uout=Uocv(SSOC)-R0I-U1-U2]? ?,? ? ? ? ?(2)

式中:[Qn]為電池的標稱容量,[SSOC]為電池SOC? ? ?的值。

將式(1)進行在k時刻的離散化并化簡得到:

[xk=Axk-1+BIT, k,yk=Cxk+DIT, k+Uocv(SSOC).]? ? ? ? ?(3)

其中:[x=U1, U2, SSOCT],[y=Uout],

[A=1-TsR1C10001-TsR1C10001],[B=TsC1TsC2-TsQn],

[C=[-1, -1, 0]],[D=-R0],k表示時刻k中每個變量的值;[Ts]代表采樣周期,是從k時刻到k+1時刻的時間,是常量值;[IT, k]表示k時刻的電池工作電流。

1.2? ?模型參數辨識

采用最小二乘法[10]識別式(1)和式(2)中電池的未知參數[R0]、[R1]、[R2]、[C1]、[C2]和非線性函數[Uocv(SSOC)]。選取三星公司生產的INR18650-30Q鋰離子電池,10個電池并聯成電池組作為研究對象,實驗平臺為ITECH公司生產的ITS5300電池充放電測試系統(ITS 5300 Battery Test System)。電池基本參數:標稱容量為3 000 mA·h,標稱電壓為3.6 V,放電截止電壓為2.5 V,充電截止電壓為4.2 V。對電池進行恒流脈沖放電實驗,實驗中電池端電壓曲線圖如圖2(a)和圖2(b)所示,將安時積分法所得值設為端電壓真實值,相對比得到端電壓模型值與真實值誤差如圖2(c)所示。由圖2可以看出,兩者誤差值基本穩定在±0.05%以內,模型準確,可以進行電池模型參數辨識。

1.2.1? ? 電池歐姆內阻辨識

參數辨識的原理取決于電池的物理行為。根據圖2可知,一旦放電電流執行或停止,端電壓立即下降或上升,因為電容器[C1]、[C2]的電壓[U1]、[U2]在開始放電時不會突然改變,所以歐姆內阻[R0]可以通過開始/結束放電時端電壓的變化來確定。在圖2中選取一個放電周期,選取[ta]、[tb]、[tc]、[td]時間點的端電壓,可以計算歐姆內阻[R0]:

[R0=|VT(tb)-VT(ta)|+|VT(td)-VT(tc)|2|I|].? ? ? (4)

1.2.2? ? ?電池極化電阻和極化電容辨識

參數[R1]、[R2]、[C1]、[C2]的識別分為兩步:第一步確定時間常數[τ1]、[τ2];第二步識別參數[R1]、[R2]、[C1]、[C2],給出公式用以辨識極化電阻和極化電容:

[U(t)=U(t0)e-t-t0τ+IR(1-e-t-t0τ)]? ,? ? ? ? ?(5)

式中:[τ=RC],[t0]為松弛過程的初始時間。

通過使用MATLAB函數“自定義方程”在曲線擬合工具箱,確定了時間常數[τ];通過式(7)確定[U1(tc)]、[U2(tc)]。參數求解公式如下:

[R1=U1(tc)I(1-e-tc-taτ1),R2=U2(tc)I(1-e-tc-taτ2),C1=τ1R1,C2=τ2R2.]? ? ? ? ? ? ? ? ? ? ? ?(6)

電池參數辨識完成,結果如表1所示。

1.3? ? OCV-SOC曲線實驗標定

采用曲線擬合的方法對非線性函數OCV-SOC進行辨識[11],通過間歇恒流脈沖放電實驗并將電池在放電結束后靜置2 h,確保松弛過程結束。為了準確地擬合測量數據,采用6階多項式:

[Uovc(SSOC)=4.723 6S6SOC-10.556S5SOC+3.877 1S4SOC+7.131 5S3SOC-7.695 9S2SOC+]

[3.802 4S1SOC+2.869 1.]? ? ? ? ? ?   ? ? ? ? ? ? ? ? ? ? ? (7)

該多項式的擬合結果圖如圖3所示。

1.4? ?擴展卡爾曼濾波算法

卡爾曼濾波[12]是一種高效率的遞歸濾波器,能夠從一系列不完全包含噪聲的測量中估計動態系統的狀態,其適用范圍為線性系統。在非線性系統中,通常采用擴展卡爾曼濾波(EKF)[13-14],EKF的基本思想是將非線性系統線性化。對非線性函數的Taylor展開式進行1階線性化階段,忽略其余高階項,從而將問題轉化為線性。

EKF算法在估算電池狀態時有很大的局限性。電池在使用的過程中,電池內部的鋰離子活躍性會逐漸變差,電池壽命會隨著參與化學反應的鋰離子數目減少而逐漸衰減。等效電路模型是根據電池內部的電化學反應而建立,當電池老化時,電池內部電化學反應程度發生變化,導致電池模型參數產生相應的改變。EKF算法在估算電池狀態時將電池內部參數當作定值,使得估計結果會有較大誤差。

隨著電池的使用,電池老化不可避免。由于電池老化產生的電池內部參數改變會導致估計結果產生較大的偏差,為了最大程度地減小此類偏差產生的影響,在EKF算法的基礎上進行改進,讓等效電路模型參數進行在線自我調節,使電池模型更加貼合真實情況下的電池,在保證模型的基礎上提高估計精度。

2? ? MDEKF算法的實現

提出多尺度并行擴展卡爾曼算法(MDEKF)對電池狀態和參數進行估計。該算法將設計2個并行的擴展卡爾曼濾波器(EKF)分別用于狀態濾波和參數濾波。考慮到電池內部參數隨時間改變程度較小,在多尺度的構想下分別對2個濾波器設定不同的采樣時間,以便達到減小計算量,使結果更加快速收斂的目的。

在電池模型建立中已經將電池模型進行了離散化處理,電池模型的狀態和參數的空間方程如? ? ?式(8)—式(12)所示:

[dθdt=θk+1-θkTs]? ,? ? ? ? ? ? ? ? ? ? ? ? (8)

[xk+1=F(xk, θk, Ik)=]

[1-TsR1C10001-TsR1C10001xx+TsC1TsC2-TsQnI+ωxk] ,? (9)

[θk+1=θk+ωθk]? ? ? ,? ? ? ? ? ? ? ? ? ? (10)

[Uout, k=G(xk, θk, Ik)=Uocv(SSOCk)-U1-U2-R0I+υk] ,? ? (11)

[xk=[U1, k, U2, k, SSOCk]T,θk=[R0, k, R1, k, R2, k, C1, k, C2, k, Qk].]? ? ? ? ? ? (12)

其中:下標[k]表示變量在步驟[k]中的值;[Ts]、[ts]為2個時間尺度下的采樣時間;[x]表示電池狀態,[θ]表示所有電池模型需要迭代的參數,[ωxk]和[ωθk]是電池狀態和參數在步驟k時的過程噪聲,與電池狀態和參數無關,它們的協方差矩陣分別為[Qx]和[Qθ],與式(18)中的[Qk]不為同一變量;[?xk]和[?θk]是電池狀態和參數在步驟[k]時的測量噪聲,它們的協方差矩陣為[Rk]。

多尺度并行擴展卡爾曼濾波算法(MDEKF)就是同時運行狀態濾波器和參數濾波器,在不同時間尺度下對電池狀態和參數進行運算,用電池容量表征SOH,從而完成該算法對SOC-SOH的聯合估算。算法實現結構如圖4所示。

狀態和參數的估計方程和具體實現步驟為:

Step 1? 算法初始化,對[θ0、x0、Pθ0、Px0]賦初始值。實現途徑見式(13):

[Pθ0=E[(θ-θ0)(θ-θ0)T], θ0=E[θ],Px0=E[(x-x0)(x-x0)T], x0=E[x].]? ?(13)

Step 2? 計算參數濾波器的時間更新方程。參數濾波器的采樣時間為[Ts=6 s],步長[k∈1, 2, …, ∞]。時間更新方程見式(14):

[P-θk=Pθk-1+Qθ, θ-k=θk-1].? ? ? ? ? ? ?(14)

Step 3? 計算參數濾波器的測量更新方程。測量更新方程見式(15),得到[θk]后將該值傳遞給狀態過濾器進行計算,在一個采樣周期后,狀態過濾器返回執行Step 2更新電池參數。

[Kθk=P-θk(Cθk)T[CθkP-θk(Cθk)T+Rk]-1,θk=θ-k+Kθk[Yk-G(xk, θ-k, Ik)],Pθk=P-θk-P-θkKθkCθk.]?   ? (15)

其中:[Cθk=?G(xk, θk, Ik)?θ|θ=θk].

Step 4? 計算狀態濾波器的時間更新方程。參數濾波器的采樣時間為[ts=0.1 s],步長[k∈1, 2, …, ∞]。時間更新方程見式(16):

[x-k=F(xk-1, θ-k-1, Ik-1)],[P-xk=Ak-1Pxk-1ATk-1+Qx]. (16)

其中:[Ak-1=?F(x, θk-1, Ik-1)?x|x=xk-1].

Step 5? 計算狀態濾波器的測量更新方程。測量更新方程見式(17),該式計算可得第k步的SOC,在一個采樣周期后,狀態過濾器返回執行Step 4。因為2個濾波器采樣周期不同,當2個濾波器采樣時間在宏觀時間尺度上相等時,狀態濾波器接受參數濾波器更新的參數[θk+1],繼續更新SOC。

[Kxk=P-xk(Cxk)T[CxkP-xk(Cxk)T+Rk]-1,xk=x-k+Kxk[Yk-G(xk, θ-k, Ik)],Pxk=P-xk-P-xkKxkCxk.]? ? ?(17)

其中:[Cxk=?G(x, θk, Ik)?x|x=xk].

在Step 3中更新了電池參數,提取電池容量參數Q,對SOH進行估計。在Step 3運行結束后返回Step 2更新參數,重復該過程,提取容量Q并更新SOH估計值([SSOH])。[SSOH]估計公式如式(18)所示:

[SSOHk=QkQ0]? .? ? ? ? ? ? ? ? ? ? ? ? (18)

上述步驟中,[θ0、x0、Pθ0、Px0]是算法的初始值,[θ0]可以由電池實驗數據計算得到,[x0]是電池狀態初值,[Pθ0、Px0]分別是參數濾波器和狀態濾波器的初始系統參數協方差矩陣,表示每次迭代的交付誤差。[Qx]和[Qθ]是過程噪聲的協方差矩陣,表示模型參數的初始誤差,[Rx]表示測量誤差。因此,通過降低[Q]值,使模型對電池模型的精度提高;通過降低[R]值,使算法實測值更加準確。[K]為卡爾曼增益,右上角的“-”表示算法中變量的預估值,頂部標簽“^”表示算法中變量的估計值。

如果狀態觀測器的[Qx]取值過小,算法將不再考慮端電壓的實測值,此時算法近似于庫倫計數法。當SOC初始值誤差較大時,得到的結果無法收斂;若[Rx]取值過小,算法將完全依賴于實測值,此時算法近似于開路電壓法,精度無法達到預期。對于參數過濾器來說,如果[Qθ]取值過小,算法將會認為初始電池參數的取值無限逼近真實值,將會使參數過濾器不再更新電池參數,使算法變為EKF;若[Rθ]取值過小,會使參數濾波更新范圍變大,導致算法發散。所以為保證算法性能,[Qx]、[Rx]、[Qθ]、[Rθ]應取合適的值,具體取值思路及極端取值造成的影響參考文獻[15]。本文取值為:

[Qx=10-600010-600010-6,Rx=1,Px=100010001,Qθ=diag(10-4, 10, 10-3, 10, 10-3, 10-3),Rθ=0.1,Pθ=diag(10-4, 10, 10-3, 10, 10-3, 10-2).]? ? (19)

3? ? 實驗驗證與分析

通過UDDS行駛工況對提出的MDEKF算法進行驗證。將離線估計的SOC作為真實值,并用EKF算法與本算法作比較,驗證算法精度與穩定性;將容量估算值作為表征量,估計出SOH,通過討論SOH的估計值再次反證SOC估計的準確性,從而完成雙擴展卡爾曼濾波的SOC-SOH估算。

圖5是MDEKF算法和EKF算法的SOC估計結果,MDEKF算法和EKF算法的SOC初始值設為0.8。從圖5中可以看出,2種算法都能在初始SOC不正確的情況下追蹤真實SOC,MDEKF算法的SOC估計精度的平均誤差控制在1.10%,比EKF算法的平均誤差1.43%降低了0.33%,而最大誤差也由EKF的2.68%降低至2.26%。由時間初始時間段可以看出,MDEKF算法和EKF算法都可以在SOC初始值錯誤的情況下,估計值收斂于真實值,而對算法的采樣時間進行多尺度思想優化后,MDEKF算法的收斂速度遠遠高于EKF算法。2種估計算法的數值結果如表2所示。

圖6是在UDDS工況下采用MDEKF算法進行的電池實際容量[Q]估計和SOH估計。由于該工況時間周期較短,電池充放電次數還不足以使電池老化、出現額定容量下降的現象,所以默認在該工況運行中電池SOH真實值為1。由[Q]估計結果得知,電池在工況實驗內估計容量值較好收斂于電池額定容量;SOH估計值相較于真實值,溢出偏移誤差最大值為0.58%,虧損偏移誤差為0.41%,估計結果總體誤差不超過0.99%,由此證明了電池處于健康狀態,驗證了算法的預測性能準確。估計誤差數值如表3所示。

4? ? 結論

本文在2階RC等效電路的基礎上提出了多尺度并行擴展卡爾曼濾波算法(MDEKF),通過仿真測試驗證算法的估計性能。結果表明,MDEKF算法可以實現SOC-SOH的聯合估算,保證了估計精度,同時降低了運算量,使得SOC估計可以在初始值錯誤的情況下,能夠快速作出反應并收斂于真實值,避免了因初始值錯誤而導致無法對電池狀態進行估計。估計結果收斂后與EKF估計結果相比,平均誤差由1.43%降低到1.10%。

參考文獻

[1]? ? ?馬華,從長杰,王馳偉.儲能用鋰離子動力電池研究進展[J].化學工業與工程,2014,31(3):26-33.

[2]? ? ?孫道明.動力鋰離子電池SOC和容量估計方法研究[D].杭州:浙江大學,2021.

[3]? ? ?錢建文,杜翀,田欣,等.一種改進T-S模糊神經網絡估計鋰電池SOC的方法[J].電源技術,2020,44(9):1270-1273.

[4]? ? ?宮兵,凌六一,何業梁,等.基于自適應無跡卡爾曼濾波的鋰電池SOC估計[J].電源技術,2020,44(11):1594-1599.

[5]? ? ?劉楠,葉望博,李明,等.自適應平方根球型卡爾曼濾波下的SOC估計[J].電源技術,2021,45(2):173-176,239.

[6]? ? ?王志福,李仁杰,李霞.基于混合AUKF和HIFF的鋰離子電池SOC估計[J].電池,2021,51(4):380-384.

[7]? ? ?郭佑民,戴銀娟,付石磊.城軌車用復合動力儲能系統蓄電池SOC和SOH估計[J].鐵道科學與工程學報,2020,17(11):2920-2928.

[8]? ? ?唐傳雨,韓華春,史明明,等.基于DEKF的儲能電池系統SOC估計方法研究[J].電力工程技術,2021,40(3):7-14.

[9]? ? ?吳小慧,張興敢.鋰電池二階RC等效電路模型參數辨識[J].南京大學學報(自然科學版),2020,56(5):754-761.

[10]? ?陳息坤,孫冬.鋰離子電池建模及其參數辨識方法研究[J].中國電機工程學報,2016,36(22):6254-6261.

[11]? ?ZHENG F D,XING Y J,JIANG J C,et al.Influence of different open circuit voltage tests on state of charge online estimation for lithium-ion batteries[J].Applied Energy,2016,183:513-525.

[12]? ?彭丁聰.卡爾曼濾波的基本原理及應用[J].軟件導刊,2009,8(11):32-34.

[13]? ?萬亞坤,李陽春,馬浩天,等.基于擴展卡爾曼濾波算法的鋰電池SOC估計[J].蓄電池,2020,57(5):243-246,250.

[14]? ?王黨樹,王新霞.基于擴展卡爾曼濾波的鋰電池SOC估算[J].電源技術,2019,43(9):1458-1460.

[15]? ?ZHU Q,XIONG N,YANG M L,et al.State of charge estimation for lithium-ion battery based on nonlinear observer:an H∞ method[J].Energies,2017,10(5):1-19.

Estimation of battery state parameters based on a multi-scale double extended Kalman filter algorithm

KONG Dehao1,LIU Shengyong*1,2

(1.School of Electrical, Electronic and Computer Science, Guangxi University of Science and Technology,

Liuzhou 545616, China; 2. Guangxi Key Laboratory of Automobile Components and Vehicle Technology(Guangxi University of Science and Technology), Liuzhou 545006, China)

Abstract: Accurate estimation of state of charge (SOC) and state of health (SOH) of the battery and? ? rational utilization can extend its life. Multi-scale double extended Kalman filter (MDEKF) algorithm was adopted based on extended Kalman filter algorithm to improve the accuracy of SOC-SOH online estimation. Firstly, 2-RC equivalent circuit model was established. Secondly, the least square method was used to identify the model parameters. Then, a filter with double structures was designed to? ? ? ? ? ?estimate the battery SOC and modify the parameters, and battery SOH was estimated with battery pack capacity as representation. The simulation results verify the fast convergence and real-time performance of the algorithm. The SOC estimation error decreases from 1.43% to 1.10%, and the SOH estimation? ?result is stable within 0.5%.

Key words: SOC-SOH estimation; multi-scale double extended Kalman filter; 2-RC equivalent model

(責任編輯:黎? 婭)

猜你喜歡
卡爾曼濾波模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
改進的擴展卡爾曼濾波算法研究
測控技術(2018年12期)2018-11-25 09:37:34
基于遞推更新卡爾曼濾波的磁偶極子目標跟蹤
3D打印中的模型分割與打包
基于模糊卡爾曼濾波算法的動力電池SOC估計
電源技術(2016年9期)2016-02-27 09:05:39
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
基于擴展卡爾曼濾波的PMSM無位置傳感器控制
電源技術(2015年1期)2015-08-22 11:16:28
基于EMD和卡爾曼濾波的振蕩信號檢測
主站蜘蛛池模板: 精品国产一区91在线| 亚洲免费福利视频| 亚洲三级视频在线观看| 国产亚洲欧美在线专区| 国产主播福利在线观看| 精品人妻一区二区三区蜜桃AⅤ| 少妇精品在线| 不卡无码网| 久久黄色影院| 精品在线免费播放| 香蕉视频在线观看www| 欧美亚洲国产一区| 国产小视频免费| 欧美成人影院亚洲综合图| 日韩中文无码av超清| 99免费在线观看视频| 自拍亚洲欧美精品| 黄色网页在线播放| 精品无码一区二区三区电影| 高清欧美性猛交XXXX黑人猛交| 国产系列在线| 日韩精品无码一级毛片免费| 国产色偷丝袜婷婷无码麻豆制服| 91小视频在线| 国产精品污视频| 欧美午夜视频在线| 色偷偷av男人的天堂不卡| 午夜a级毛片| 亚洲bt欧美bt精品| 国产毛片片精品天天看视频| 国产精品久久久久无码网站| 亚洲最新网址| 日韩麻豆小视频| 国产青榴视频| 无码国内精品人妻少妇蜜桃视频| 色首页AV在线| 97精品伊人久久大香线蕉| 女人18毛片水真多国产| 人妻精品全国免费视频| 九九九精品成人免费视频7| 欧美日韩国产系列在线观看| 色屁屁一区二区三区视频国产| 一级全黄毛片| 久草网视频在线| 88av在线看| 免费一级全黄少妇性色生活片| 国产成人免费| 无码高潮喷水在线观看| 综合色天天| 亚洲第一成年网| 精品视频一区在线观看| 中美日韩在线网免费毛片视频| 国产最爽的乱婬视频国语对白| 第九色区aⅴ天堂久久香| 亚洲欧美日韩中文字幕在线一区| 国产小视频在线高清播放| 色欲不卡无码一区二区| 欧美午夜在线视频| 欧美97色| 中国一级特黄大片在线观看| 在线99视频| 国产成年女人特黄特色毛片免| 日本精品视频一区二区| 国产激爽爽爽大片在线观看| 污污网站在线观看| 欧美在线天堂| 婷婷六月综合网| 狠狠五月天中文字幕| 天天干伊人| 欧美成人第一页| 性激烈欧美三级在线播放| 在线免费看片a| 亚洲成aⅴ人片在线影院八| 扒开粉嫩的小缝隙喷白浆视频| 一本大道东京热无码av| 福利在线不卡一区| 性69交片免费看| 亚洲天堂视频网站| 久久网综合| 欧美中文一区| 欧美在线中文字幕| 国产午夜人做人免费视频中文 |