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

基于數據驅動的電動汽車電池安全風險預測*

2023-05-29 10:00:04楊博聞程雅鈺
汽車工程 2023年5期
關鍵詞:特征模型

胡 杰,余 海,楊博聞,程雅鈺

(1.武漢理工大學,現代汽車零部件技術湖北省重點實驗室,武漢 430070;2.武漢理工大學,汽車零部件技術湖北省協同創新中心,武漢 430070;3.新能源與智能網聯車湖北工程技術研究中心,武漢 430070)

前言

隨著電動汽車的快速普及,消費者和汽車廠商對電池安全問題越來越重視,電池安全領域正在成為研究的熱點,以實現對電池安全風險的準確預測。Hong 等[1]通過對電池溫度實時監測并運用Z-score方法,提出了一種基于溫度異常的電池安全實時預測方法,實驗結果表明該方法能預測電池組發生溫度風險的位置和時間。Yang等[2]利用電壓信息估計短路電流的人工神經網絡(ANN),獲得所需電流,利用該電流預測ESC 單元的最高溫升及溫度分布,實驗結果表明該方案能對電池風險進行精度較高的診斷預測。詹華等[3]提出了一種基于LS-SVR 模型的電池風險預測方法,通過預測電壓不同狀態來評估動力電池的風險,實驗結果表明該方案能對動力電池風險進行有效的診斷預測。姚建平[4]分析了影響電池溫度上升的原因并基于EEMD-AR 模型提取特征,有效地提高了鋰電池風險診斷的準確度。張家昌[5]利用小波包分解的方法去除數據噪聲和提取特征,以此方法分解數據作為輸入的模型,對電池風險分類預測的準確率達94.13%。Feng 等[6]總結了處理數據和提取特征的方法,對電池不一致性診斷和評估的方法進行了分析歸納,并提出了自己的思路。Zhang 等[7]根據開路電壓和內阻分析了電池組的一致性,建立了內膜標準偏差與電池組OCV 標準偏差之間的函數關系,實驗結果表明該方案誤差在6%以內。Jiang等[8]設計了一種分析電池一致性的建模方法,該方法能夠描述電池一致性參數間的依賴結構和參數的統計特性,與現有的方法相比能提高描述準確度。王慶年等[9]通過挖掘數據構建健康因子,并基于ARIMA 模型對健康因子進行分析預測,結果表明該方案能夠對車輛進行精度較高的有效預警。

本文中基于大量電動汽車的實時數據和周圍環境數據來分析電池安全風險預測問題,并將預測問題解耦為機器學習預測和時間序列預測問題,最后通過Stacking 模型將二者進行集成改進。構建接近實際情況的多指標特征,模擬了電池復雜的使用環境;基于數據挖掘提取了行駛工況和駕駛風格特征以提高方案的泛化性和魯棒性;對基模型進行篩選優化后建立Stacking 集成模型以提高精度,最終實現了基于數據驅動的電動汽車電池安全風險的準確預測。

1 數據與問題分析

1.1 樣本數據

本文所有數據為北京地區車輛的運行數據,數據采集時間是2019年10月1日至2020年8月20日,車輛類型為純電動運營汽車,共計6 974 萬條數據,采集周期為10 s。共有21 種不同類型數據,其中部分數據如表1所示。

表1 部分車輛運行數據

1.2 問題分析

目前市面上的純電動汽車的電池絕大部分是鋰電池[10-11]。電池的內部材料和各種性能在車內外復雜的環境下和每次循環充放電的過程中不斷發生變化。

為模擬真實的電池使用場景,本文中基于電池充放電對電池安全的影響、使用環境對電池安全的影響、電池狀態即荷電狀態(state of charge,SOC)和健康狀態(state of health,SOH)對電池安全的影響,進行特征構造、模型篩選和Stacking 集成。本方案大致步驟為:對實車數據進行處理,將數據劃分為3個不同片段,并分為訓練集和驗證集;基于對數據的理解來進行特征提?。贿\用隨機森林和SHAP 組合模型篩選出最終特征;基于預測問題的解耦,分別進行模型的篩選優化;將最終篩選出的模型進行Stacking集成改進,并對模型進行數據實驗。具體的預測流程如圖1所示。

圖1 電池風險預測流程圖

2 數據清洗

在數據的實際采集過程中,存在環境變化、設備不穩定、傳輸延遲等問題,故本文數據存在大量的問題數據。為了后續模型構建的需要、減小模型誤差和提高模型準確度,需要對數據進行清洗。

2.1 片段劃分

為了便于對問題數據進行處理,本文將數據劃分為不同的狀態片段:行駛片段、停車充電片段、停車靜置片段。具體的劃分規則如下:(1)車輛狀態保持不變;(2)充電狀態保持不變;(3)考慮到數據集的數據分布特點,為了更好地劃分充電、靜止和放電片段,數據間間隔不能小于240 s。劃分后的片段內數據少于100 條時,考慮其分析意義不大,且增加計算時長,故對其進行刪除處理。具體劃分流程如圖2所示。

圖2 片段劃分流程圖

2.2 數據處理

針對存在缺失值、數值為0、采集滯后、車輛狀態異常標記等不同問題數據,本文進行了差異化處理。

2.2.1 缺失數據處理

在數據中最高單體電壓、最低單體電壓和最高溫度存在大量缺失,缺失量分別為35%、47%和56%,針對這類數據有兩類及兩類以上的缺失進行刪除處理,然后計算缺失數據前后各5 個共10 個數據的標準差以判斷變化程度,當數據變化程度小時,取這10個值的平均值進行填充。此時這3類缺失量分別為8.5%、9%、11%,對這3 類缺失數據基于隨機森林算法進行回歸填補,具體的回歸填充流程圖如圖3所示。填補后的平均相對誤差如表2所示。

圖3 隨機森林填補缺失值

表2 平均相對誤差

2.2.2 異常數據處理

(1)在數據中存在大量數值為0 的異常值,具體如圖4 所示。從圖4 中可以看出,數據存在電壓、SOC、單體電池最小溫度和單體電池最大溫度為0的異常情況。針對此部分的異常數據,用上述缺失數據處理的方法進行填補。

圖4 數值為0的異常數據圖

(2)在數據中某些片段交界處會出現異常數據,例如前后兩個片段分別為行駛片段和停車充電片段,在停車充電片段開始處電流為-10 A,但車輛速度顯示大于0。對該部分數據分析后,認為是采集滯后的原因,利用速度數據來進行平移修正。

(3)在數據中存在車輛狀態標記異常的數據,例如車輛的充電狀態為停車充電,但運行狀態卻顯示行駛狀態(圖5)。考慮車輛運行狀態對于后續模型構建十分重要,并且該數據占總數據比重較大,因此不能刪除處理。分析后在已有的車輛運行狀態特征的基礎上利用速度和電流等特征進行修正處理。

圖5 車輛狀態標記異常圖

(4)將通用報警標志進行二進制轉換后,發現數據中只存在單體電池欠壓、SOC 低、一致性差報警、制動系統報警、SOC 跳變5 種通用報警。并發現放電過程中的數據,大量耦合了電池報警信息,且該部分數據的電壓、單體溫度等也出現了偏低或偏高,為了防止對后續模型預測造成影響,將放電過程數據進行刪除。

3 特征工程

3.1 多指標特征提取

為了更加真實地模擬電池實際工作場景,本文爬取了周圍環境數據作為補充參數,并在數據的理解上構建了多指標特征。

(1)時間信息特征提取

可根據解析出的時間特征得到是否為晚上、白天等類別特征,這些特征對電動汽車的安全風險存在一定影響。例如2019-11-01-06:31:46 特征轉換如表3所示。

表3 時間信息特征轉換

(2)電池健康狀態特征提取

本文根據電池健康狀態的標準定義進行健康狀態特征提取,其標準定義如下:

式中:CM為電池容量;CN為電池額定容量。

其中電池容量采用安時積分法進行計算:

由于該方法計算電容存在誤差,為了減少誤差提高模型精度,需要對電池容量進行修正。考慮到溫度對電池內部的化學反應以及能量消耗有很大影響,相比于常溫環境,低溫狀態的電池容量會有明顯減少,本文先基于不同溫度下電池容量值與溫度的變化關系,對電池容量進行修正,得到修正后電池容量和溫度的變化,如圖6所示。

圖6 單體電池平均溫度和修正后的容量值變化圖

電池放電電流的變化對電池的容量也有較大影響,電池的放電電流會影響到電池的內部化學反應從而加速容量衰減。本文基于計算的電池容量值與電池容量真實值之間的比值,對電池容量進行修正,得到電池容量在修正前后的差值與放電電流的關系,如圖7所示。

圖7 電流修正后電池容量差值與放電電流關系圖

電池容量經過修正后仍存在部分離群點,針對離群點則采用箱型圖來進行分析檢測,檢測結果如圖8 所示。由箱型圖可知,將邊緣線之外的離群點予以剔除。

圖8 各月電池容量箱型圖

(3)能耗信息特征提取

電動汽車在實際運行過程中,電池容量會不斷衰減從而影響電池安全性能,本文基于式(3)來計算電池能耗并提取電池能耗特征[12]。

式中:SOC(t0)和SOC(t1)分別為在t0和t1時刻的SOC值;U為放電過程中的電壓值;Cn為電池容量值;Qd為制動回收部分貢獻的能量,計算公式見式(4)。采用放電片段之前的充電片段來計算放電片段的電池容量,計算公式見式(2)。整車能量消耗中的汽車制動能量主要是制動能量回收貢獻的能量[13]。本文根據汽車制動能量回收貢獻率和各部分的效率來計算制動回收貢獻的能量。

式中:Ui表示輸出電壓;IiiN表示制動電流的大??;ηes、ηmd、ηtd、ζ分別表示儲能裝置平均能量效率、電機驅動系統效率、傳動系統效率、制動能量貢獻率,這些系數分別為0.902、0.801 1、0.901、0.539。

(4)電池一致性特征提取

電池是由多個單體電池組成的整體而非單個單體電池,各單體電池本身的連接方式、性能、溫度、電壓等都存在著差異[14]。電池不一致現象會對電池產生很大的影響,由于數據維度的原因,本文只分析了電池的溫度一致性和電壓一致性。

本文用單體溫度極差和單體電壓極差來表示電池的溫度一致性和電壓一致性。單體電池溫度一致性的計算公式如下:

式中:Tempmax、Tempmin分別表示單體溫度最大值和單體溫度最小值;Umax、Umin分別表示單體電壓的最大值與最小值分別表示溫度一致性得分系數和電壓一致性得分系數,且都為0.2,一致性得分越接近1,說明電池一致性就越好。

(5)行駛工況特征提取

基于前文片段劃分的結果,本文選取片段內的平均車速、駐車時間比、勻速時間比、加速時間比4個特征用于工況片段聚類分析。其中平均車速vm、駐車時間比例Pi、勻速時間比例Pc、加速時間比例Pa的具體計算公式如下:

式中:n是片段長度;vi是第i時刻的速度;ns是速度為0且加速度絕對值小于0.1的長度;nc為加速度小于0.1的長度;na為加速度大于0.1的長度。

本文基于K-means 聚類對工況片段進行聚類,運用“肘部法”得到SSE 和聚類數的關系,如圖9 所示。從圖9 中可知,當聚類數為3 時,SSE 的下降坡度明顯變緩,后續變化不大,因此確定最佳類簇數目為3。通過分析并可視化片段中的平均速度與駐車比例之間的關系,得到的聚類結果如圖10 所示。在圖10 中可以看出行駛工況分為3 類,用1、2、3 表示:1 代表的工況駐車比例高,平均速度低,劃分為城市擁堵工況;2 代表的工況駐車比例適中,平均速度適中,劃分為市郊工況;3 代表的工況駐車比例較低,平均速度較高,劃分為高速工況。

圖9 肘部法確定類簇數圖

圖10 K-means聚類結果圖

(6)駕駛風格特征提取

在實際駕駛過程中,不同的駕駛員對應不同的駕駛風格,不同的駕駛風格會對行駛狀態產生不同的影響,并進一步影響電池狀態。由于行駛工況會對駕駛風格產生一定影響,故需要針對不同的工況來對駕駛員的風格進行分析。駕駛員主要是通過制動和加速踏板將行駛意圖傳遞給車輛,根據踏板的開合程度與開合度變化率可以評估其駕駛風格,踏板的開合度和開合度變化率可通過汽車的加速度與加速度變化率來表示[15]。

目前多基于邏輯回歸、支持向量機和神經網絡算法等對駕駛風格進行識別,但本文數據采樣間隔為10 s,采樣頻率太低,用多種模型識別后發現效果不佳,故本文選擇片段內加速度標準差和加速度變化率標準差作為評判指標。為避免在不同工況下出現識別偏向性,對不同工況加速度標準差和加速度變化率標準差取平均值作為修正系數SDJJ和SDaa。將片段內標準差除以所屬工況的標準差平均值作為評判標準:

式中:n表示加速度絕對值大于0.1 的片段長度;Ji表示第i個加速度絕對值大于0.1 的加速度變化率大?。籥i表示第i個加速度絕對值大于0.1 的加速度大??;m表示所屬工況的片段長度。

用加速度標準差與不同工況下的加速度標準差平均值的比值來反映踏板開合度情況,用加速度變化率的標準差與不同工況下加速度變化率標準差的平均值的比值來反映踏板開合度的變化程度。根據踏板開合度、踏板開合度變化率和實際用車情況,劃分的駕駛風格如表4所示。

表4 駕駛風格劃分

3.2 多指標特征篩選

為了方便對后續模型的統一處理和提高模型的穩定性,本文對特征的連續型變量采用歸一化處理,特征的離散型變量采用獨熱編碼。

為了精簡特征,提高模型的泛化能力,需要剔除無關特征和冗余特征。本文先運用隨機森林模型對特征進行初篩選,將重要度得分靠后的特征剔除;利用SHAP 計算特征對模型結果的貢獻度;最后用隨機森林和SHAP 模型相組合的方法來對多指標特征進行篩選,具體特征篩選流程如圖11 所示。特征最終的篩選結果如圖12所示。

圖11 RF-SHAP特征篩選流程

圖12 RF-SHAP特征篩選結果

在篩選非線性關系的特征中,樹模型的篩選效果較好,但單棵決策樹則容易出現過擬合的現象,故用隨機森林模型對特征進行初步篩選,以避免出現過擬合現象??紤]到隨機森林模型的可解釋性不強,而SHAP 的方法能夠進行很好解釋,故用隨機森林模型對多指標特征進行初步篩選后,再用SHAP方法計算特征對模型結果的具體貢獻值,具體公式如下:

式中:g為解釋模型;z′表示相應特征是否能夠被觀察到,其中結構化數據的每個特征都能被觀察到;M為輸入特征的數目;?0是用于解釋模型的固定常數。由于隨機森林模型輸入的是結構化數據,因此z′應該是為1的向量,則

對于特征初篩中的隨機森林模型來說,需要計算出每個特征的Shapley 值即該特征的貢獻值后對所有可能的特征組合進行加權求和:

式中:S為模型中特征的子集;M為特征數量;fx(S∪i)和fx(S)分別表示特征子集S∪i和特征子集S的預測;fx(S∪i) -fx(S)則表示特征的邊際貢獻值為特征子集S的權重,其與順序存在一定的關系。

根據圖12 中的特征篩選結果,提取出重要程度最高的13 個特征作為最終輸入特征,具體如表5 所示(表5 還包括其他特征)。在最終特征中電池信息類的特征重要度得分最高,對模型效果影響最大,而能耗特征、行駛工況特征、駕駛風格特征和天氣季節特征等對模型效果也有一定影響。

表5 篩選后的最終特征

4 基于多指標電池安全風險預測模型

4.1 預測方法

考慮到電池安全狀況是在復雜車載環境下,在使用過程中隨著電池各項性能指標變化而變化,也在使用的過程中隨著使用時間的增加而不斷變化,即電池安全風險預測與電池關鍵性能指標有關,也與時間有關,本文中將電池安全風險預測問題解耦為傳統機器學習預測問題和時間序列預測問題兩種。對于前者采用基于多指標的機器學習預測方法,后者是在機器學習模型的基礎上加入時間特征進行處理。并對兩個模型建立新的Stacking 模型進行預測,評價指標選擇受試者工作特征(receiver operating characteristic,ROC)曲線下的面積(area under curve,AUC)值以保證模型預測精度以及魯棒性。

表5 特征中max_alarm_lvl 即安全風險等級分為0、1、2、3 的4 種不同風險等級,其中:0 級是安全狀態;1 級是一種相對安全的狀態,并不影響行駛;2 級風險表示電池已經處于相對危險的狀態,需要駕駛員停車仔細檢查;3 級風險是最高級別風險,需要駕駛員立即停車并下車,車輛很可能出現熱失效和失火現象。以篩選出的13 個特征作為模型輸入,以安全風險等級作為模型預測輸出,對車輛風險進行預測。

分類模型輸出一般為概率值,設置的閾值不同,準確度也會不同,模型預測輸出分為4 類:正類樣本被預測成正類,即真正類(True Positive,TP);正類樣本被預測成負類,即假負類(False Negative,FN);負類樣本被預測正類,即假正類(False Positive,FP);負類樣本被預測負類,即真負類(true negative,TN)。具體關系如表6 所示,AUC 為ROC 曲線下的面積,ROC 是以False Positive Rate(FPR)為橫坐標、True Negative Rate(TNR)為縱坐標的曲線,其定義如式(18)和式(19)所示,故AUC 的值越趨近于1,模型預測效果越好。

表6 模型輸出值之間的關系

4.2 模型篩選

模型篩選是為了能快速簡潔地對多個模型進行對比擇優,從而免去單個模型的繁瑣。對于機器學習模型的篩選,本文中基于劃分好的訓練數據集來進行交叉驗證,其計算公式如式(20)所示;對于時間序列模型的篩選,則是基于自定義時間的交叉驗證法來對模型進行驗證。

(1)機器學習模型篩選

對機器學習模型的訓練、篩選和驗證不包括time 和max_alarm_lvl_shift 的時間特征。具體的交叉驗證結果如圖13 所示。由圖可知XGBoost 和隨機森林模型預測精度靠前,雖然XGBoost 的模型運行時間稍長但精度高,且運行時間在可接受范圍內,故本文選擇XGBoost 作為機器學習模型的最終 模型。

圖13 機器學習模型篩選結果

(2)時間序列模型篩選

對時間序列模型的篩選即在機器學習模型特征的基礎上加入time 和max_alarm_lvl_shift 的時間特征,采用自定義時間的交叉驗證方法,根據月份以4∶1 的比例劃分好數據進行交叉驗證。驗證結果如圖14所示。從圖中可看出隨機森林和model-cb模型精度靠前,但model-cb模型的運行時間過長,故本文選擇隨機森林模型作為時間序列模型的最終模型。

圖14 時間序列模型篩選結果圖

4.3 模型集成

利用貝葉斯優化法來對模型進行參數優化,貝葉斯優化法是通過對種群進行宏觀迭代的方法來實現參數尋優,可以轉換為未知函數在已知空間的最優解問題,對于集成模型和樹模型,貝葉斯參數優化法是一種很有效的方法[16]。優化后的結果如表7 所示。從表7 中可以看出兩個模型效果都有了一定提升,XGBoost 模型AUC 得分提升了1.100%,RF 模型AUC得分提升了1.122%。

表7 模型前后優化對比表

RE 模型為時間序列模型其輸入特征多了time和max_alarm_lvl_shift 兩個時間特征,XGBoost 模型為機器學習模型其輸入與時間序列無關,本文Stacking集成模型的方法為:為充分利用數據擴大訓練樣本提高預測精度,將訓練數據劃分為4 個月、5個月和6 個月的數據來訓練時間序列RF 模型,將得到的RE 模型預測后1 個月的數據;將全部訓練數據作為輸入變量來訓練XGBoost 模型,得到機器學習模型XGBoost 對數據進行預測,最后將3 個RF 模型與XGBoost 模型一起作為基模型輸入到簡單的邏輯回歸LR 模型中,集成得到最終的預測結果。具體Stacking集成模型的結構如圖15所示。

圖15 Stacking集成模型結構圖

圖15 中:X為XGBoost 模型的全部數據輸入變量,X4、X5、X6為RF 模型的4 個月、5 個月、6 個月的輸入變量,其中不同輸入變量結構如圖16所示。

圖16 不同輸入變量的結構圖

4.4 數據實驗

為了驗證Stacking 集成模型的優越性,采用圖17 中的方法對其進行數據實驗。首先根據月份將訓練數據劃分為4 個月、5 個月和6 個月的數據,將其作為不同輸入變量來訓練3 個不同的RF 模型;將全部訓練數據作為輸入變量來訓練XGBoost 模型,最后將3個不同RF模型與XGBoost模型的結果進行集成,得到最終預測結果。

圖17 Stacking集成模型實驗流程圖

為了進一步驗證模型的性能,本文基于測試數據對XGBoost、RF 和Stacking 集成模型進行模型測試,其模型測試結果如表8 所示。由表8 可知XGBoost、RF 模型在測試數據集上的AUC 得分與訓練集上的結果類似,表明各個模型沒有出現過擬合,模型結果符合預期并具有良好的泛化性和魯棒性。

表8 各個模型測試結果表

其中XGBBoost、RF 和Stacking 集成模型的ROC曲線如圖18 所示。由圖18 可知Stacking 集成模型相比前兩者單模型能明顯提高效果,可以對電動汽車電池風險進行更準確的預測。

圖18 各個模型ROC曲線圖

5 結論

本文中通過提取汽車內外多維度信息,用隨機森林和SHAP 組合模型的方式對特征進行篩選,模擬了復雜的電池應用場景??紤]了汽車實際工況和駕駛員的駕駛習慣,使方案有較好的泛化性和魯棒性。將電池安全風險預測問題解耦為傳統機器學習預測問題和時間序列預測問題,分別采用XGBoost模型和RF 模型進行預測,并將兩者進行Stacking 集成改進后明顯提高了預測精度,根據集成模型的預測結果AUC值為0.991 7,表明該方案能對電動汽車的電池安全風險做出準確的預測,可以為安全化、智能化的電池管理系統提供輔助決策信息。

猜你喜歡
特征模型
一半模型
抓住特征巧觀察
重要模型『一線三等角』
新型冠狀病毒及其流行病學特征認識
重尾非線性自回歸模型自加權M-估計的漸近分布
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 亚洲—日韩aV在线| 久久国产精品影院| 久久精品国产精品青草app| 久久亚洲中文字幕精品一区| 日韩高清中文字幕| 91精品国产一区自在线拍| 2021国产精品自产拍在线| 五月激情综合网| 东京热一区二区三区无码视频| 尤物亚洲最大AV无码网站| 手机永久AV在线播放| 真实国产乱子伦高清| 伊人激情综合网| …亚洲 欧洲 另类 春色| 日韩午夜福利在线观看| 她的性爱视频| 无码aaa视频| 国产在线自乱拍播放| 在线欧美日韩| 欧美精品v| 91网址在线播放| 亚洲va在线∨a天堂va欧美va| 香蕉综合在线视频91| 日韩 欧美 国产 精品 综合| 免费看的一级毛片| 一本大道香蕉中文日本不卡高清二区| 丁香婷婷激情网| 久99久热只有精品国产15| 欧美精品1区| 91精品视频播放| 99尹人香蕉国产免费天天拍| 欧美亚洲一二三区| 久草国产在线观看| a级高清毛片| 亚洲黄色网站视频| 狠狠做深爱婷婷久久一区| 91视频99| 久久婷婷综合色一区二区| 亚洲视频影院| 国产拍在线| 国产精品极品美女自在线| 精品久久国产综合精麻豆| 男人天堂亚洲天堂| 福利在线不卡一区| 久久久受www免费人成| 日本午夜精品一本在线观看| 国产亚洲欧美日韩在线观看一区二区| 毛片三级在线观看| www.youjizz.com久久| 日韩国产黄色网站| 五月天福利视频| 伊人欧美在线| 无码精品福利一区二区三区| 国产成人综合网| 日韩午夜片| 无码啪啪精品天堂浪潮av| 真实国产乱子伦视频| 自拍欧美亚洲| 精品欧美视频| 99视频在线免费| 色综合成人| 九色视频线上播放| 波多野结衣视频一区二区| 中文字幕永久视频| 免费在线不卡视频| 国产成人久久综合777777麻豆| 色婷婷亚洲综合五月| 国产一级毛片高清完整视频版| 99re在线观看视频| 国产真实自在自线免费精品| 国产在线日本| 先锋资源久久| 欧美色伊人| 国产丰满大乳无码免费播放| 亚洲第一黄色网| 国产特一级毛片| 久久99国产精品成人欧美| 亚洲免费三区| 一本大道香蕉高清久久| 国产a v无码专区亚洲av| 欧美不卡二区| 久久国产精品国产自线拍|