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

電池等效電路模型極化參數估計的收斂性分析

2025-03-07 00:00:00張維剛曾嘉博周維陳正吳頌潘文軍
湖南大學學報·自然科學版 2025年2期

摘要:為探究動力電池極化參數估計在不同工況下的誤差收斂特性,采用李雅普諾夫第二法對估計器的穩定性和收斂性進行分析. 基于等效電路模型和遺忘因子遞歸最小二乘法推導得到極化參數估計誤差的狀態方程. 采用李雅普諾夫第二法對狀態方程進行分析,得到估計器漸進收斂的必要條件,即持續變化的電流輸入,并提出一種圖解法用于分析誤差的動態收斂特性和論證此必要條件的合理性. 采用實驗校準的電池模型產生的數據對理論分析過程和結果進行驗證. 結果表明,估計器在持續變化的電流輸入下能逐漸收斂到真值附近,且在變化劇烈的正負交變工況下具備更好的收斂性.

關鍵詞:電池管理系統;等效電路模型;參數辨識;漸進收斂性;李雅普諾夫第二法

中圖分類號:TM912.8 文獻標志碼:A

在電動汽車、智能電網等許多應用中,動力電池內部狀態和參數估計是能源管理系統的關鍵. 準確的電池參數估計依賴于精確且計算高效的電池模型及算法[1]. 近年來隨著研究的深入,參數估計算法層出不窮,電池模型也朝著高精度、高效率的方向發展. 然而,模型參數的在線辨識仍是該領域的一大難題,研究者提出了多種在線估計等效電路模型(equivalent circuit model,ECM)參數的算法,包括最小二乘法及其變體(如遞歸最小二乘法和加權最小二乘濾波器)[2]、雙擴展卡爾曼濾波器[ 3]、自適應無跡卡爾曼濾波器[ 4]、粒子濾波器[5]等. 但這些研究成果并未在工業界得到廣泛應用,最重要的原因是參數在線辨識的魯棒性欠佳,辨識結果在某些情況下攜帶較大的誤差甚至發散. 鑒于此,有學者對在線辨識時參數的可辨識性及其與輸入數據特性的關系展開了研究.

Sharma 等[6]采用Fisher 信息矩陣分析了一階ECM中內阻和容量參數的可辨識性,并通過克拉默-拉奧(Cramer-Rao, CR)邊界量化了聯合估計誤差的下界與傳感器誤差、電流大小和電池開路電壓(opencircuit voltage, OCV)? 荷電狀態 (state of charge,SoC) 曲線斜率的關系. Rothenberger等[7]使用同樣的方法分析了參數CR邊界與電流軌跡的幅值和頻率的關系,并據此設計了一條可提高辨識精度的電流軌跡. 在上述研究的基礎上,Lin等[8]分析了SoC、容量和內阻在典型電流輸入下單獨估計和聯合估計的CR邊界變化,并得到了聯合估計相比單獨估計時放大系數的表達式. Lin[9]推導了極化參數CR邊界關于電流的頻率響應表達式. Song等[10]在此基礎上得到了OCV、歐姆內阻和極化參數的CR邊界關于電流的頻率響應表達式,并測試了不同電流輸入以及不同參數聯合估計時的估計精度. 研究者意識到輸入數據特征對參數辨識精度的重要性,Song等[11]提出了使用多頻率濾波器分離電流以提高辨識精度的方法,并據此設計了一種順序估計算法[12]. 各種噪聲和參數初值誤差也會影響估計精度. Lin[13]分析了存在測量偏差時SoC的穩態收斂情況. Shen等[14]分析了容量和SoC的初值存在誤差時SoC誤差的穩態收斂值與修正增益之間的關系. Zhou等[15]基于參數的置信區間評估了二階ECM中各參數的可辨識性,并探究了數據量、采樣頻率和傳感器噪聲等對可辨識性的影響. Andersson等[16]指出,并非所有模型參數都能被準確辨識,也并非所有參數都與模型性能相關,參數辨識不僅需要了解其所針對的物理系統,還應理解參數辨識背后的理論和概念.

上述研究雖量化了參數可辨識性與輸入數據之間的關系,但并未解決參數估計器在不同工況下的收斂性這一重要問題. 為填補該項研究空白,本文基于李雅普諾夫第二法推導得到了估計器漸進收斂的必要條件,并創新性地提出了一種圖解法,詳細分析估計器的收斂機制. 本文的研究成果對電池模型參數在線估計算法的設計具有一定指導意義.

1 參數估計誤差的狀態方程

1.1 基于一階ECM 的參數估計

本節基于一階ECM 設計電池極化參數的在線估計器[17-18]. 一階ECM結構如圖1所示.

根據基爾霍夫定律和電池SoC定義,可得一階ECM的狀態方程和輸出方程如下:

Vp,k = Vp,k - 1e-Δt/R1C1 + Ik - 1 R1(1 - e-Δt/R1C1 ) (1)

ξk = ξk - 1 + ηc Ik - 1Δt (3 600C)-1 (2)

Uk = Voc(ξk ) + Vp,k + Ik R0 (3)

式中:ξ 為電池的SoC; Voc(ξk )為電池的開路電壓,是關于SoC的單調遞增函數;Vp,k 和Uk 分別為極化電壓和端電壓;R1 和C1 表示極化電阻和極化電容;R0 為歐姆內阻;C 為電池容量;ηc 為庫倫效率;Δt 為采樣周期;I 為電流,規定充電方向為電流正方向;k 為采樣時刻.

由式( 1) 可知,R1、C1與極化電壓之間的關系為非線性,直接辨識這兩個參數十分困難[19-20]. 因此,本文采用文獻[21]中的方法,通過辨識兩個與R1、C1相關的參數來完成對R1、C1 的辨識. 令

令 yk = Vp,k,將式(1)轉化為如下辨識方程:

yk = yk - 1 p1,k + Ik - 1 p2,k (5)

很明顯,式(5)是關于待辨識參數的線性方程.為了便于閱讀,將式(5)表示成矩陣形式:

yk = φk θk (6)

式中:φk 和θk 分別為辨識方程的輸入向量和待辨識參數向量.

參數修正過程可寫成:

θ?k = θ?k - 1 + Kk (y mk - y?k ) (8)

式中:Kk=[ K1,k, K2,k ]T 為估計器的更新增益;y mk 為yk的偽測量值,由電壓、電流的測量值以及SoC、R0 的估計值計算得到;右上標m 表示變量的測量值;上標 ^ 表示變量的估計值.

y mk = Uk - Voc(ξk ) - Ik R0 (9)

不失一般性,本文選擇最基本的遺忘因子遞歸最小二乘法(forgetting factor recursive least squares,FFRLS)作為辨識算法. 算法的增益向量Kk 和協方差矩陣Pf,k 的更新方程如下:

式中:λ f 為辨識過程的遺忘因子,是一個常數.

1.2 誤差狀態方程

為分析估計誤差的動態收斂特性,需先得到其狀態方程,再根據李雅普諾夫第二法[22-23]對系統的穩定性進行分析. 首先,記各變量的誤差如下:

式中:右上標*表示變量的真值;Δ表示變量的誤差;eI,k 和eU,k 分別為電流和電壓的測量噪聲.

聯立式(1)、式(5)~式(7)、式( 11),經簡化整理后可得:

式中:χ 為OCV - SoC 曲線在相應SoC 區間內的斜率.

將式(12)代入式(8),可得參數估計誤差的狀態方程,并整理成如下矩陣形式:

估計器的穩定性由系統矩陣決定,估計器穩定,誤差可收斂. 然而,在誤差收斂過程中,會受到各種外部誤差的干擾,影響估計誤差最終的收斂程度,也有可能影響估計器的穩定性. 從式(14)中可以看到,外部誤差既會影響系統矩陣,又會引入額外估計誤差,最終引起穩態誤差. 雖然外部干擾會影響系統矩陣,但系統矩陣結構并未改變,系統穩定性是否受影響還需進一步分析.

2 誤差動態特性的理論分析

2.1 估計器的漸近穩定性條件分析

得到估計器的誤差狀態方程的系統矩陣后,便可根據李雅普諾夫第二法分析系統的穩定性. 估計器本質上是一種線性時變的離散系統,該系統的能量函數可設為V (Δθk ):

從式(24)中可以看到,參數誤差在修正過程中表現為:估計值的誤差Δθ?k - 1 向目標值的誤差Δθk 所在的直線φkΔθk = 0(后文簡稱目標直線)不斷靠近(圖2). 經過k 時刻的修正,能使Δθ?k 到達目標直線附近. 若下一時刻的輸入φk + 1 = φk,目標直線不變,那么目標值的誤差也不變,即Δθk + 1 = Δθk. 由式(17)可知,此時V?(Δθk + 1) = V?(Δθk ) = 0,誤差能量保持不變. 在電池參數辨識中,也有類似的結論[10-11]——恒流輸入時估計器無法進一步收斂;相反,若φk + 1 ≠ φk,目標直線變化,Δθ?k 向目標值的誤差Δθk + 1所在的直線φk + 1Δθk + 1 = 0 不斷靠近,此時V?(Δθk + 1) ≠ V?(Δθk ). 由此可知,在不同輸入下,V?(Δθk )和V?(Δθk + 1)不會同時等于0.

從上述證明過程可知,在單步修正過程中會存在一種P 矩陣使得Q 矩陣半正定,即使輸入φ 攜帶誤差也不會對Q 矩陣的正定性造成影響. 由此可知,雖然外部干擾影響系統矩陣,但系統穩定性不受其影響. 在連續兩步修正過程中,不斷變化的輸入φ 是系統漸近穩定的必要條件.

2.3 基于圖解法的漸近穩定性分析

為了進一步分析參數估計器的收斂過程,本節利用圖解法來描述能量函數、誤差修正方向和目標直線的幾何特征. 由式(16)可知,能量函數的幾何形狀由Pk 決定,為一系列同心橢圓,且誤差Δθ1,k、Δθ2,k越大,能量橢圓面積越大,誤差能量也就越大,反之亦然. 由式(14)可知,誤差修正方向為兩個時刻的誤差點形成的向量,其斜率km 為:

當輸入一定時,整個平面的誤差修正方向都是相同的,且誤差都沿著修正方向往目標直線φkΔθk = 0靠近.

由(20)可知,系統在單步修正過程中的V?(Δθk )不恒小于0,即誤差能量的變化是不確定的.因此,本節將結合上述幾何特征來定性分析誤差能量的變化情況. 當φk 和Pk 一定時,能量函數、誤差修正方向和目標直線也隨之被確定. 結合圖3(a)進行分析,在k 時刻,每一條誤差修正方向總會與能量橢圓相切并形成一個切點,多個切點可以連成一條直線. 由切點連成的直線與目標直線將誤差平面劃分為3個區域(A、B、C). 若初始誤差Δθk - 1 位于A、C 兩部分時,經過k 時刻的修正,Δθk - 1 將沿著誤差的修正方向抵達Δθk. 顯然,由Δθk 形成的能量橢圓面積小于Δθk - 1,誤差能量減小,誤差漸近收斂;反之,誤差漸近發散. 因此,根據構造的Pk 和Δθk - 1 所在位置的不同,誤差有可能是漸近收斂的,也有可能是漸近發散的.

若改變Pk 矩陣使得由切點連成的直線與目標直線所夾的B 區域面積越小,誤差就越有可能收斂,如圖3(b)所示. 若在k 時刻合理構造一種Pk矩陣,使得由切點連成的直線與目標直線重合,如圖3(c)所示,便可使除了目標直線上誤差能量不變,其余位置的誤差能量都減小,誤差將最大可能地漸近收斂. 因此,通過構造Pk 矩陣使單步修正過程中V?(Δθk )不大于0的方法是可行的,這與前文的分析是一致的.

結合誤差能量的變化,誤差收斂原理圖解如圖4所示. 當輸入不斷發生變化時,目標直線持續旋轉,誤差沿著修正方向朝目標直線靠近,誤差能量逐漸減小,誤差在趨近目標直線的過程中逐漸接近0.

3 動態工況下的仿真與實驗驗證

通過對估計器的誤差收斂特性進行分析可知,變化的輸入是估計器漸近收斂的必要條件,且誤差能量在收斂的過程中會跟隨修正方向不斷減小. 本節選用2.6 Ah/3.6 V的三元鋰電池設計實驗與仿真,對上述理論分析進行驗證. 實驗采用新威CT-4004-5V100A型高性能電池檢測系統作為電池監測設備.在真實工況中無法實時獲取極化參數的真值,為了更好地評價估計器的估計效果,需構建一個一階ECM,用于生成仿真工況數據.

3.1 電池模型參數獲取與實驗驗證

電池模型參數可通過連續充放電脈沖實驗進行離線辨識獲取初值,再由動態工況進行優化調整獲取最終的阻抗參數,如圖5所示. 其中,單個脈沖曲線包含10 s的1 C放電、40 s的靜置以及10 s的1 C充電,通過對放電脈沖回彈的電壓曲線進行擬合可得阻抗參數的初值. 間隔10% SoC重復上述流程,可得電池模型阻抗參數的初值. 隨后采用US06工況進行動態測試,對阻抗參數進行優化調整獲取最終電池阻抗參數及OCV-SoC曲線,分別如圖6(a)和圖6(b)所示. 調整的目的是使離線辨識的阻抗參數在動態工況中更為適用,進一步減小模型端電壓誤差.

調整的步驟為:

1)繪制US06工況下模型的端電壓誤差曲線.

2)設置電壓誤差閾值的絕對值為20 mV.

3)對于超過閾值的SoC 區間,固定歐姆內阻、OCV-SoC曲線不變,通過調整對應SoC區間的極化電容與極化電阻,從而減小端電壓誤差.

經優化調整后,構建的一階電池模型產生的仿真電壓與實驗電壓的最大電壓誤差小于20 mV,如圖6(c)所示. 為進一步驗證該模型在其他工況下的適用性,采用動態壓力測試工況(dynamic stress test,DST)、中國輕型車行駛工況(China light-duty vehicletest cycle,CLTC)對電池模型產生的端電壓進行驗證,如圖7所示. 由圖7可知,該電池模型在不同工況下產生的最大電壓誤差均小于30 mV,平均絕對誤差小于6 mV,證明該模型對動態工況具有良好的適應性,可用于后續仿真驗證.

3.2 估計器漸近收斂性的仿真驗證

由電池模型生成動態工況下的電壓數據作為電壓真值. 為了模擬外部干擾,在電壓真值中添加均值為2 mV、方差為3 mV 的測量噪聲[25]作為電壓測量值,并輸入參數估計器中進行參數辨識. 同時,為了更好地對比不同工況下估計器的收斂程度,生成了兩種典型動態工況,如圖8所示.

圖9展示了在不同工況中參數估計器參數誤差變化軌跡. 由圖9可知,在動態工況1中,經過一段時間的修正,估計器誤差能收斂到0附近,表明該估計器在動態輸入下具有良好的收斂性. 最后誤差在0附近震蕩,表明外部干擾的存在最終只會影響估計器的穩態特性. 由于工況2中的電流變化一開始較為平緩,其估計誤差先是增加,而后才逐漸減小. 很明顯,工況2中誤差最終的收斂程度不如工況1,且p1的誤差顯著小于0. 該結果驗證了前文理論分析的正確性,即估計器在變化的輸入下能漸近收斂. 同時,該結果進一步表明,工況的變化越劇烈,估計器的收斂性會越好.

以工況1為例對估計器的動態收斂過程進行詳細分析,根據參數誤差的動態變化軌跡,可繪制誤差修正方向和誤差能量橢圓圖,如圖10所示. 由圖10可知,估計器的誤差修正方向與理論分析的結果是一致的,整個誤差平面的修正方向以目標直線為界同側一致,異側相反. 由于外部干擾的存在,修正方向無法完全平行于目標直線和能量橢圓等高線交點的切線,只能近似平行該切線,但外部干擾的存在并沒有影響估計器的收斂性. 當輸入不斷變化時,誤差沿著修正方向不斷向目標直線靠近,直至誤差減小到零附近. 工況2 的收斂原理類似,故不再重復展示.

4 結 論

本文基于遞歸最小二乘法推導了鋰電池一階等效電路模型極化參數估計誤差的狀態方程,明確了系統矩陣是估計器穩定性的決定因素. 采用李雅普諾夫第二法對估計器的系統矩陣進行分析,得出持續變化的電流輸入是估計器漸近收斂的必要條件.提出了一種圖解法對誤差的修正方向和收斂軌跡進行了分析,分析表明誤差在沿著修正方向朝目標直線靠近的過程中逐漸收斂,誤差能量得到減小. 采用基于實驗數據校準的電池模型對理論分析進行了仿真驗證,結果表明估計器在持續變化的輸入下能漸近收斂,且誤差會跟隨其修正方向不斷收斂到零附近.

本文研究表明,估計器在持續變化的工況下,可展現出更好的收斂性. 后續的研究可基于這一結論開發自適應的工況篩選策略和算法,以提升估計器的收斂性能.

參考文獻

[1] HAVANGI R. Adaptive robust unscented Kalman filter withrecursive least square for state of charge estimation of batteries[J].Electrical Engineering, 2022, 104(2): 1001-1017.

[2] ZHENG Y J,CUI Y F,HAN X B,et al. Lithium-ion batterycapacity estimation based on open circuit voltage identificationusing the iteratively reweighted least squares at different aginglevels[J].Journal of Energy Storage,2021,44:103487.

[3] LIU Z,CHEN S H,JING B Q,et al. Fractional variable-ordercalculus based state of charge estimation of Li-ion battery usingdual fractional order Kalman filter[J].Journal of Energy Storage,2022,52:104685.

[4] HOU J,LIU J W,CHEN F W,et al.Robust lithium-ion state-ofchargeand battery parameters joint estimation based on anenhanced adaptive unscented Kalman filter[J]. Energy,2023,271:126998.

[5] XING L K,LING L Y,GONG B,et al. State-of-chargeestimation for lithium-ion batteries using Kalman filters based onfractional-order models[J]. Connection Science,2022,34(1):162-184.

[6] SHARMA A, FATHY H K. Fisher identifiability analysisfor a periodically-excited equivalent-circuit lithium-ion batterymodel[C]//2014 American Control Conference, June 4-6,2014.Portland, OR, USA: IEEE, 2014: 274-280.

[7] ROTHENBERGER M J,ANSTROM J,BRENNAN S,et al.Maximizing parameter identifiability of an equivalent-circuitbattery model using optimal periodic input shaping[C]//ASME2014 Dynamic Systems and Control Conference,October 22-24,2014. San Antonio,Texas,USA: American Society of MechanicalEngineers,2014.

[8] LIN X F,STEFANOPOULOU A G. Analytic bound on accuracyof battery state and parameter estimation[J]. Journal of theElectrochemical Society,2015,162(9):A1879-A1891.

[9] LIN X F. Analytic analysis of the data-dependent estimationaccuracy of battery equivalent circuit dynamics[J].IEEE ControlSystems Letters,2017,1(2):304-309.

[10] SONG Z Y,HOFMANN H,LIN X F,et al. Parameteridentification of lithium-ion battery pack for different applicationsbased on Cramer-Rao bound analysis and experimental study[J].Applied Energy,2018,231:1307-1318.

[11] SONG Z Y,WANG H,HOU J,et al. Combined state andparameter estimation of lithium-ion battery with active currentinjection[J]. IEEE Transactions on Power Electronics, 2020,35(4): 4439-4447.

[12] SONG Z Y,HOU J,LI X F,et al. The sequential algorithm forcombined state of charge and state of health estimation of lithiumionbattery based on active current injection[J]. Energy,2020,193:116732.

[13] LIN X F. Theoretical analysis of battery SOC estimation errorsunder sensor bias and variance[J]. IEEE Transactions onIndustrial Electronics,2018,65(9):7138-7148.

[14] SHEN P, OUYANG M G, HAN X B,et al. Error analysis ofthe model-based state-of-charge observer for lithium-ionbatteries[J].IEEE Transactions on Vehicular Technology,2018,67(9): 8055-8064.

[15] ZHOU W,HUANG R J,LIU K,et al. A novel interval-basedapproach for quantifying practical parameter identifiability of alithium-ion battery model[J]. International Journal of EnergyResearch, 2020, 44(5): 3558-3573.

[16] ANDERSSON M,STREB M,KO J Y,et al. Parametrization ofphysics-based battery models from input-output data:a review ofmethodology and current research[J].Journal of Power Sources,2022,521:230859.

[17] NASERI F,SCHALTZ E,STROE D I,et al. An enhancedequivalent circuit model with real-time parameter identificationfor battery state-of-charge estimation[J].IEEE Transactions onIndustrial Electronics, 2022, 69(4): 3743-3751.

[18] HUA X,ZHANG C,OFFER G.Finding a better fit for lithium ionbatteries:a simple,novel,load dependent,modified equivalentcircuit model and parameterization method[J].Journal of PowerSources,2021,484:229117.

[19] SHU X,CHEN Z,SHEN J W,et al.State of charge estimation forlithium-ion battery based on hybrid compensation modeling andadaptive H-infinity filter[J]. IEEE Transactions on TransportationElectrification, 2023, 9(1): 945-957.

[20] WU J J,FANG C,JIN Z Y,et al.A multi-scale fractional-orderdual unscented Kalman filter based parameter and state of chargejoint estimation method of lithium-ion battery[J]. Journal ofEnergy Storage,2022,50:104666.

[21] CHEN Z G,ZHOU J X,ZHOU F,et al. State-of-chargeestimation of lithium-ion batteries based on improved H infinityfilter algorithm and its novel equalization method[J].Journal ofCleaner Production,2021,290:125180.

[22] 王長江,姜濤,陳厚合,等.基于相位校正李雅普諾夫指數的電力系統暫態電壓穩定評估[J].電工技術學報,2021,36(15):3221-3236.

WANG C J,JIANG T,CHEN H H,et al. Transient voltagestability assessment of power systems based on phase correctionmaximum Lyapunov exponent [J] . Transactions of ChinaElectrotechnical Society,2021,36(15):3221-3236.(in Chinese)

[23] ZHOU P,HU X K,ZHU Z G,et al. What is the most suitableLyapunov function?[J]. Chaos,Solitons amp; Fractals,2021,150:111154.

[24] WANG R N,WANG Q W,LIU L S.Solving a system of sylvesterlikequaternion matrix equations[J]. Symmetry,2022,14(5):1056.

[25] CUI Z R,CUI N X,WANG C Y,et al.A robust online parameteridentification method for lithium-ion battery model underasynchronous sampling and noise interference [J]. IEEETransactions on Industrial Electronics,2021,68(10):9550-9560.

基金項目:湖南省自然科學基金資助項目(2022JJ30158),Natural Science Foundation of Hunan Province(2022JJ30158);國家自然科學基金資助項目(51705139),National Natural Science Foundation of China(51705139)

主站蜘蛛池模板: a在线亚洲男人的天堂试看| 国产精品自拍合集| www.av男人.com| 亚洲A∨无码精品午夜在线观看| 国产精品密蕾丝视频| 性喷潮久久久久久久久| 亚洲福利视频一区二区| 91精品网站| 丁香亚洲综合五月天婷婷| 91亚洲精选| 久草热视频在线| 国产精品福利在线观看无码卡| 91美女在线| 98精品全国免费观看视频| 亚洲综合网在线观看| 99一级毛片| 国产一级视频在线观看网站| 国产网友愉拍精品视频| 色偷偷av男人的天堂不卡| 亚洲视频免| 青青青伊人色综合久久| 色哟哟国产精品| 免费一级成人毛片| 国产亚洲欧美日韩在线一区| 亚洲欧美综合在线观看| 国产一区二区三区在线观看视频| 玖玖精品在线| 欧美成人精品高清在线下载| 国产精品亚欧美一区二区三区| 久久一本精品久久久ー99| 99re在线视频观看| 日韩成人免费网站| 亚洲色图欧美激情| 人妻无码中文字幕一区二区三区| av一区二区三区高清久久| 91小视频在线观看| 欧美另类图片视频无弹跳第一页| 国产网友愉拍精品视频| 亚洲欧美国产高清va在线播放| 91成人在线观看| 亚洲黄色高清| 亚洲天堂成人在线观看| 曰AV在线无码| 六月婷婷激情综合| 麻豆精品国产自产在线| 成人国产精品一级毛片天堂| 国产真实乱子伦视频播放| lhav亚洲精品| 99热这里只有成人精品国产| 青草91视频免费观看| 97免费在线观看视频| 97精品久久久大香线焦| 色屁屁一区二区三区视频国产| 青青网在线国产| 韩日午夜在线资源一区二区| 国产青青草视频| 国产第八页| 日韩精品少妇无码受不了| 午夜福利视频一区| 午夜国产不卡在线观看视频| 无码福利日韩神码福利片| 亚洲人成在线精品| 国产99热| 熟妇人妻无乱码中文字幕真矢织江 | 国产乱子伦一区二区=| 亚洲Va中文字幕久久一区| 亚洲日韩AV无码一区二区三区人| 国产亚洲精品97在线观看| 午夜一级做a爰片久久毛片| 久久99精品久久久大学生| 小说 亚洲 无码 精品| av在线5g无码天天| 亚洲日韩高清在线亚洲专区| 欧美伊人色综合久久天天| 国产成人一二三| 国产高潮流白浆视频| 色有码无码视频| 日韩精品一区二区三区大桥未久| 不卡无码网| 久久黄色一级视频| 国产视频一区二区在线观看| 亚洲精品日产精品乱码不卡|