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

國Ⅵ重型柴油機NOx排放故障診斷算法研究

2023-05-26 08:06:28王志堅王秀雷趙康荏史彥曉帥石金李國祥
內燃機學報 2023年3期
關鍵詞:故障診斷模型

王志堅 ,王秀雷 ,趙康荏,史彥曉,帥石金,李國祥

(1. 濰柴動力股份有限公司,山東 濰坊 261061;2. 清華大學 車輛與運載學院,北京 100084;3. 山東大學 能源與動力工程學院,山東 濟南 250061)

柴油機匹配選擇性催化還原(SCR)技術可滿足國Ⅴ、國Ⅵ排放法規,是目前主流的NOx排放控制技術路線.NOx排放故障診斷是排放法規最重要的車載診斷系統(OBD)功能要求之一,對于國Ⅵ重型柴油機,根據GB17691—2018重型柴油車污染物排放限值及測量方法(中國第六階段),當NOx排放超過限值1.2g/(kW·h)時,在臺架的OBD認證過程中,一個冷熱態全球瞬態試驗循環(WHTC)結束前,OBD系統報NOx排放超標故障,激活轉矩限制器.國Ⅵ重型柴油機使用兩只NOx傳感器分別直接測量發動機的NOx原始排放和SCR下游NOx排放,NOx排放故障診斷受到NOx傳感器精度和信號漂移的影響,且車載NOx傳感器對NH3有交叉敏感性[1].為此國內外公司和學者做了大量研究,提出了基于MAP映射和基于物理模型的診斷方法,前者需大量的標定試驗,成本高且靜態條件下標定的MAP無法準確表達瞬態條件下的NOx排放狀態[2],基于物理模型的方法需要電控系統很強的計算能力,且部分模型參數難以從實車上獲取[3-4],從而限制了其在實際產品上的應用.

由于現有NOx排放故障診斷算法存在準確率偏低的不足,有必要應用智能算法提高準確率.數據驅動模型的方法因在開發成本和時間上的優越性,被廣泛應用于柴油機故障診斷模型的建立.戴金池等[5]基于神經網絡實現柴油機的瞬態NOx排放數據預測,Pu等[6]基于機器學習實現納米級的顆粒物分布規律預測,Najafi等[7]基于支持向量機和自適應神經模糊系統實現油耗、煙度、爆發壓力和增壓器轉速等的數據預測.孫立民[8]基于概率神經網絡(PNN),對船舶用柴油機高壓系統節流閥狀態進行故障診斷研究,診斷準確率為95%.文獻[9—10]基于神經網絡結合缸蓋振動和曲軸角加速度信號對柴油機失火進行故障診斷,基于分類的思想診斷準確率高達100%.張永祥等[11]基于加窗和卷積神經網絡算法,可以有效地對柴油機拉缸進行故障診斷,Ding等[12]基于小波包能量和深度卷積神經網絡實現對機床主軸軸承的智能診斷,Ince等[13]提出一種將數據生成和故障診斷結合起來的算法,該算法能夠自動提取故障特征數據,且具有優越的故障診斷能力,但是以上都是將穩態、瞬態分開建模,并應用不同的模型參數進行故障診斷,沒有考慮實際道路條件下的模型泛化性能.因此,有必要建立一種模型,能夠對瞬態和實際道路過程進行準確的NOx排放故障診斷,更好地適應產品的應用需求.

柴油機的運行是一個連續的過程,之前時刻的NOx排放狀態將會影響下一個時刻柴油機的NOx排放狀態,NOx排放故障診斷是一個典型的基于時間序列的診斷問題.對NOx排放進行故障診斷是典型的分類問題,受多種因素的影響,魯棒性要求也較高,支持向量機(SVM)作為一種典型的二分類機制,少數的支持向量將決定最終的分類結果,適合于魯棒性要求較高的故障診斷問題.基于以上分析,筆者針對實際情況,將NOx排放故障診斷抽象為時間序列的分類問題,提出一種基于時間窗口與SVM的NOx排放故障診斷算法,針對WHTC和實際道路(PEMS)試驗數據建立統一的故障診斷模型,結合診斷放行、延遲時間條件,使用時間窗口將特征輸入離散化,使用SVM算法將時間窗口分類,分別將WHTC循環作為訓練集和測試集;WHTC循環作為訓練集,PEMS數據作為測試集;WHTC循環與PEMS數據的集合作為訓練集和測試集,模型均具有很高的診斷準確率,說明模型具備良好的泛化性能.同時對比了PNN及長短期記憶神經網絡(LSTM)等其他神經網絡模型,基于時間窗口與SVM的NOx排放故障診斷模型具有更高的診斷準確率和泛化性能.

1 發動機參數與試驗設計

1.1 主要參數和試驗設備

表1~表3為柴油發動機、SCR后處理及整車主要技術參數.

表1 發動機主要技術參數Tab.1 Engine main technical parameters

表2 SCR后處理主要技術參數Tab.2 Main technical parameters of SCR aftertreatment

表3 整車主要技術參數Tab.3 Vehicle main technical parameters

發動機試驗臺架總體布置示意如圖1所示.試驗采用湘儀動力測試儀器有限公司CAC265交流電力測功機測量轉速和轉矩,采用AVL公司483微煙度計測量碳煙(soot)排放,且不受顆粒物中的可溶性有機物(SOF)、硫酸鹽和礦物質等成分的影響,HORIBA公司MEXA-7200D氣體分析儀測量排氣中NOx、氣態HC等.試驗使用國Ⅵ標準柴油.

圖1 發動機試驗臺架示意Fig.1 Schematic of experimental set-up

1.2 試驗設計

基于時間窗口和SVM的NOx排放診斷策略,需要兼容WHTC(OBD認證)和整車路譜工況(診斷策略的整車適應性驗證),并且區分正常SCR催化消聲器和劣化SCR催化消聲器(使用熱老化方法降低SCR轉化效率,使得NOx排放結果控制在1.0~1.2g/(kW·h)).因此,基于臺架WHTC和整車路譜工況及正常SCR催化消聲器和劣化SCR催化消聲器,分別采集試驗數據用于模型的訓練和驗證.

針對臺架,運行熱態WHTC,采集正常SCR催化消聲器的試驗數據(圖2a),NOx排放為0.35g/(kW·h),滿足法規要求,SCR上游NOx排放平均值為453.6×10-6、下游排放平均值為9.2×10-6,SCR上游排氣溫度平均值為275.0℃.圖2b為采集劣化SCR催化消聲器的試驗數據,NOx排放為1.18g/(kW·h),SCR上游NOx排放平均值為486.0×10-6、下游排放平均值為116.0×10-6,SCR上游排氣溫度平均值為271.1℃.

圖2 熱態WHTC循環下的NOx排放Fig.2 NOx emissions during hot WHTC cycle

針對整車,根據GB 17691—2018實際道路行駛測量方法(PEMS),進行整車路譜工況數據采集,發動機臺架進行路譜模擬,包括10%負荷和100%負荷兩種應用場景.圖3為路譜工況,10%負荷的路譜平均轉速為1133r/min,平均轉矩為416N·m;100%負荷的路譜平均轉速為1131r/min,平均轉矩為877N·m.整車PEMS雖然只認證整車排放,并不認證NOx排放故障診斷功能,但因為PEMS認證是基于整車的實際道路工況,適合驗證NOx排放故障診斷策略的整車適應性.因此,用PEMS認證的整車工況驗證NOx排放故障診斷策略的整車適應性.

圖3 路譜工況Fig.3 Working condition of road spectrum

采集正常SCR催化消聲器的整車路譜試驗數據,10%負荷和100%負荷的NOx排放分別為0.34g/(kW·h)、0.59g/(kW·h),滿足PEMS法規要求,10%負荷、100%負荷路譜的SCR上游NOx排放平均值分別為501.7×10-6、687.8×10-6,下游NOx排放平均值分別為3.6×10-6、25.7×10-6,SCR上游排氣溫度平均值分別為274.3℃、314.0℃.

采集劣化SCR催化消聲器的整車路譜試驗數據,10%負荷和100%負荷的NOx排放分別為1.71g/(kW·h)、1.66g/(kW·h);10%負荷和100%負荷的路譜的SCR上游NOx排放平均值分別為484.5×10-6、755.3×10-6,下游NOx排放平均值分別為100.4×10-6、99.8×10-6;SCR上游排氣溫度平均值分別為265.2℃、345.1℃.

2 試驗數據處理

2.1 相關度分析及輸入特征提取

灰色關聯度分析(GRA)是一種用于多因素分析的統計方法,在一個灰色系統中,評價結果受其他誘發因素影響的相對強弱,通過灰色關聯度來衡量相似或相異程度[13].相比其他分析方法,灰色關聯度分析在樣本多少及樣本有無規律的情況下同樣適用[5],數據處理的流程如下.

(1) 確定灰色關聯度分析法的比較矩陣(轉速、噴油量、冷卻液溫度、SCR上/下游NOx排放、SCR上/下游排氣溫度、廢氣流量、尿素噴射量、機油壓力、中冷后溫度、中冷后壓力、排氣背壓和煙度,共14個變量)和參考矩陣(NOx排放,將NOx體積流量轉化為NOx排放),將各參數與NOx排放進行關聯,來判斷不同參數變化對NOx排放的影響,即

通過除以矩陣每列的最大值將量綱歸一化為

(2) 計算灰色關聯系數為

(3) 計算灰色關聯度,比較矩陣和參考矩陣的灰色關聯度為

式中:Xi、Y0分別為比較矩陣和參考矩陣;n、i分別為輸入特征變量的個數和數據采集的條數;ξ0i(k)為第k個特征輸入的關聯系數;γ0i[Y0, Xi′(k)]為第k個特征輸入的關聯度;ρ為分辨系數,設定值為0.5.

根據式(5)可以得到灰色關聯度向量,通過大小次序可以判斷出各變量對NOx排放的影響程度,圖4所示灰色關聯度的前8位,都大于0.85,根據關聯度分析和工程經驗,將轉速、噴油量、SCR上/下游NOx排放、SCR上/下游排氣溫度、廢氣流量及尿素噴射量共8個變量作為最終的輸入特征,特征輸入能夠很好地描述柴油機的工作狀況,且符合比較容易準確獲取的原則[7-8].

圖4 灰色關聯度分析Fig.4 Grey correlation analysis

NOx排放的故障診斷結果為診斷模型的輸出結果,NOx排放為參考矩陣,最終確定的8個變量作為診斷模型的特征輸入.

2.2 時間窗口計算

圖5 為歸一化數據處理后的時間窗口持續時間以及間隔時間,為清晰地表達時間動態變化過程,將輸入特征進行預處理,取時間窗口的特征平均值、方差值分別代表發動機當前的運行狀態與穩定性,作為診斷模型的特征輸入.

圖5 特征輸入的時間窗口Fig.5 Time window for feature input

采用時間窗口法進行NOx排放診斷,測試數據從0s開始,隨著測試的進行,累積時間逐漸增加,當累積時間達到系統設定值時停止計算,間隔時間達到設定值后,重新進行時間窗口計算,如圖6所示.特征輸入計算式為

圖6 時間窗口特征輸入計算示意Fig.6 Schematic diagram of feature input calculation for time window

2.3 排放診斷放行條件

發動機在低負荷區,后處理的溫度低,轉化效率降低,易產生故障診斷的誤報錯;在高負荷區,后處理溫度高,且廢氣流量較大,后處理的化學反應處在極限條件下,也易產生故障診斷的誤報錯.基于上述問題,限定發動機NOx排放故障診斷的放行范圍如表4所示.其中,SCR排氣溫度是SCR上/下游排氣溫度的算術平均值.

表4 診斷放行范圍Tab.4 Diagnose release range

滿足診斷放行邊界后,發動機及SCR系統需要一定的延遲時間進入穩定工作狀態,狀態穩定后放行故障診斷,提高故障診斷穩定性,避免誤診斷,參考文獻[2,4—5]放行延遲時間,取值為10s.

3 NOx排放診斷模型建立

3.1 分類模型對比分析

對于如何區分正常與劣化的SCR催化消聲器,可以抽象為時間序列分類問題,用于數據分類的模型有SVM、PNN和LSTM等.其中,SVM是建立在統計學習理論上,SVM模型將向量映射到一個更高維的空間里,在這個空間里建立一個最大間隔超平面H,在超平面H的兩邊建有兩個互相平行的超平面H1、H2,平行超平面間的距離越大,分類器的總誤差越小,SVM遵循結構風險最小化原理,泛化能力強[14-15];PNN是用線性學習算法來完成非線性模式分類工作,能夠達到快速分類處理的效果,PNN模型具備結構簡潔、訓練難度小并且收斂速度快等優點[16-18];LSTM是為了解決長期以來的問題而專門設計出來的,LSTM能夠存儲并傳遞長期記憶和短期記憶,并將其用于當前的數據預測中,時間越臨近,對當前數據預測的影響權重越大,訓練工作量大,學習時間長,特別適合預測復雜的連續時間數據[5,19].

對于NOx排放故障診斷這種復雜的時間序列分類問題,將連續時間數據處理為時間窗口,弱化了窗口之間的時間先后影響,理論上并不能發揮LSTM的優勢.考慮到ECU計算能力,不可能做到實時的機器學習和參數優化,只能離線學習,在線使用,因而在考慮準確率的基礎上,還要求高泛化性能,即模型要適應WHTC循環和不同的路譜工況,而PNN通常需要更多的數據對網絡進行訓練才能進一步提高網絡模型的故障診斷準確率[16-18],PNN相對不適用于此應用場景.因此,最終選擇SVM算法研究NOx排放故障診斷.

3.2 SVM模型建立

3.2.1 SVM模型

使用SVM算法進行分類處理,SVM算法核心思想是在特征空間中找到最優分類面,使帶有不同標簽的特征向量到最優分類面的距離最大化.圖7為SVM模型架構示意,最優分類超平面H的公式為

圖7 SVM模型架構示意Fig.7 Architecture diagram of SVM model

式中:m為訓練集樣本數量;i為樣本空間中第i個樣本;w為權值矩陣;b為偏置;C為懲罰因子;a為拉格朗日系數,取值為1;iξ為松弛變量.

對樣本的輸入特征值進行標準化,消除量綱的影響,即

式中:x為標準化前的樣本點;x′為標準化后的樣本點;μ為樣本的均值;σ為樣本的標準差.將發動機NOx排放超1.2的狀態記作故障狀態,輸出標定為1;不超1.2的狀態記作正常狀態,輸出標定為0.

3.2.2 核函數選擇

常見的核函數主要包括線性核函數、多項式核函數和徑向基核函數(RBF)等,采用徑向基核函數完成低維特征到高維特征的映射[6,20],徑向基核函數通常根據樣本數據ix與樣本數據中心點xj之間的歐氏距離來定義其在高維特征中的值.通常,樣本數據距離樣本中心越近將會有更大的特征映射值,距離樣本中心越遠將會有更小的特征映射值,其表達式為

3.2.3 學習算法

由于SVM是借助二次規劃來求解支持向量,當樣本數量較大時,特征輸入矩陣計算的復雜度較大,針對此問題,采用序列最小優化(SMO)的啟發式算法[8-10],ξi的定義是降低了異常點對分類模型的影響,在實際的推導過程中轉變為0 ≤ ai≤ C 的約束條件,通過拉格朗日函數,SVM最優超平面的目標函數可優化為

式中:iα為拉格朗日因子;iy為樣本觀測值.SMO算法通過每次只優化兩個變量,將其他的變量都視為常數,通過反復迭代求得權重w與偏置b,使樣本點最終滿足Karush-Kuhn-Tucker(KKT)條件,KKT的條件為

式中:yi為第i個樣本的觀測值;wxi+b為第i個樣本的預測值;α*i為拉格朗日因子的解.

3.3 模型參數確認

3.3.1 貝葉斯優化算法

貝葉斯優化算法(Bayesian optimization,BO)是一種基于后驗概率的尋優方式[5],其核心思想在于根據高斯先驗擴充目標函數的統計學模型,利用統計模型尋找目標函數的最優位置.其優化過程如下.

(1)通過當前樣本建立高斯過程回歸模型(Gaussian process,GP),即

式中:X、y分別為樣本的屬性和觀測值;X*為待預測的樣本集合;,若考慮到輸入數據的噪聲方差 nσ,則有,I為單位矩陣,K由平方指數協方差函數k(x,x′|θ)所得,一般描述的公式為

(2) 選取目標函數所對應的最優解作為下一個迭代更新的指標,通過擴充樣本數據來完成迭代尋優的算法,當前最優解的均值與方差將會決定下次尋優.貝葉斯尋優過程如式(19)所示.

式中:μ、σ分別是樣本空間D中點x點的均值與方差;β用于調節平均值和標準差的重要程度,取值為1.為了避免陷入尋優局部最小值,在迭代尋優過程中若當前最優解的方差σ<tσ(設定tσ為0.5),將判定當前尋優陷入局部最優解,此時將修改高斯過程模型中協方差函數 k(x, x′|θ)的核參數θ,將核參數θ乘以當前迭代次數,并更新尋優模型.通過反復迭代直到滿足迭代終止條件,獲得目標函數的最優解.

3.3.2 參數優化

SVM算法的權重矩陣w與偏置b在訓練的過程中能夠自動地迭代更新,而模型的徑向基因子g和懲罰因子C是超參數,合理地選擇超參數能夠有效地提升分類模型的準確性[5-6].利用貝葉斯優化算法優化分類模型超參數,設在D維參數空間中,通過擴充滿足邊界條件內的目標函數最優值,最終到達模型最優位置,如圖8所示.貝葉斯優化SVM模型主要包括初始化、目標函數計算、更新尋優模型和迭代尋優等操作[6-7].

圖8 故障診斷模型的建模過程示意Fig.8 Modeling process of fault diagnosis model

懲罰因子C與徑向基因子g構造尋優參數Xi=(CXi, gXi),CXi的變化范圍為[0.0001,100000],gXi的變化范圍為[0.0001,100000],最大進化次數為50,損失函數定義為

式中:n為驗證集樣本數量;y?j為預測值;yj為實值;I{}為指示函數;wj為樣本觀測值權重.

利用5折交叉驗證進行完成模型尋優[5],損失函數如圖9a所示.當時間窗口為105s時,超參數優化過程如圖9b所示,經貝葉斯優化后,懲罰因子C為2.94,徑向基因子g為1.96,WHTC訓練集上診斷準確率從86.63%提高到99.61%.

圖9 貝葉斯優化過程Fig.9 Bayesian optimization process

4 NOx排放故障診斷結果及分析

4.1 診斷結果

根據式(21)計算NOx排放故障診斷準確率.

式中:η為準確率;rw為診斷正確的窗口數;sw為有效窗口總數.

對比了時間窗口時長對診斷準確率的影響,如圖10所示.其中,WHTC-PEMS表示“-”前面數據作為訓練集,“-”后面數據作為測試集.隨著窗口時間的延長,模型預測準確率逐漸上升并最終趨于穩定,繼續延長時間窗口意義不大,甚至準確率反而下降,拐點大約在105s.結合工程應用,推薦時間窗口取值為105s,此時利用WHTC作為模型訓練集,在PEMS數據上的故障診斷準確率達到96.99%.其他測試集與訓練集的建模組合,其模型故障診斷準確率均可達到99%以上,說明SVM模型具備良好的泛化性能.

圖10 時間窗口時長對診斷準確率的影響Fig.10 Influence of time window duration on diagnostic accuracy

4.2 WHTC與PEMS工況診斷準確率對比分析

WHTC與PEMS工況的數據分布對比如圖11所示.對于轉速分布規律,WHTC與PEMS工況相比,WHTC、10%負荷PEMS工況和100%負荷PEMS工況在低轉速工況區間(1000r/min以下)的數據占比分別為16.7%、6.4%和14.3%;在中間轉速工況區間(1000~1300r/min)的數據占比分別為74.4%、92.3%和 65.5%;在高轉速工況區間(1300r/min以上)的數據占比分別為8.9%、1.2%和20.2%.WHTC的低轉速工況區間數據比較完整,中間、高轉速工況區間數據相對不完整.對于噴油量分布(與負荷分布相當),WHTC與PEMS工況相比,WHTC、10%負荷PEMS工況和100%負荷PEMS工況在低噴油量工況區間(75mg/ cyc以下)其數據占比分別為40.5%、57.5%和24.0%;在中間噴油量工況區間(75~150mg/cyc)的數據占比分別為44.2%、42.5%和51.3%;在高噴油量工況區間(150mg/cyc以上)的數據占比分別為15.3%、0.0%和24.7%.WHTC在低、中間噴油量工況區間數據比較完整,高噴油量工況區間數據相對不完整.

圖11 WHTC與PEMS工況數據分布對比Fig.11 Comparison of distribution of working condition data between WHTC and PEMS

只利用WHTC數據訓練故障診斷模型,在PEMS數據進行驗證,由于訓練數據在中間、高轉速工況區間與高噴油量工況區間的數據相對不完整,會影響模型在對應工況區間的故障診斷準確率.WHTC數據作為訓練集,PEMS數據作為測試集,故障診斷模型在10%負荷、100%負荷PEMS工況的故障診斷準確率分別為98.67%和94.81%,如表5所示,高負荷的診斷準確率相對較低,與分析結果相符.不同SCR狀態下,故障診斷模型在柴油機PEMS工況數據的診斷準確率分別為97.29%、96.69%,正常與裂化SCR催化消聲器的診斷準確率差異不大,見表6.

表5 不同負荷下柴油機PEMS工況數據的故障診斷準確率Tab.5 Fault diagnostic accuracy of diesel engine PEMS operating data under different loads

表6 不同狀態下柴油機PEMS工況數據的故障診斷準確率Tab.6 Fault diagnostic accuracy of diesel engine PEMS operating data under different conditions

4.3 不同診斷模型對比分析

在Matlab中建立PNN模型[11-13]與LSTM模型[5],PNN拓撲結構為:1層輸入層,5層隱藏層,1層分類層;LSTM的拓撲結構為:1層輸入層,1層LSTM層,2層全連接層,1層分類輸出層.采用相同的方法訓練并優化超參數,PNN的超參數主要為隱藏層神經元個數與徑向基因子,LSTM的超參數主要為LSTM層神經元個數、批處理尺寸及學習率.表7為不同模型對PEMS試驗數據的診斷準確率,說明SVM模型在面對未知輸入時,依然能夠利用訓練集數據,結合時間窗口計算和放行條件,使得SVM模型相對PNN和LSTM具有更高的泛化能力.

表7 不同模型的故障診斷準確率Tab.7 Fault diagnosis accuracy of different models

5 結論

(1) 建立了基于時間窗口和SVM的NOx排放故障診斷模型,以轉速、噴油量、SCR系統上/下游NOx排放、排氣流量、SCR系統上/下游排氣溫度和尿素噴射量作為模型輸入特征,利用WHTC循環和PEMS試驗數據進行訓練和驗證,能夠取得較高的診斷準確率.

(2) WHTC循環試驗數據作為訓練集,貝葉斯參數優化后,模型對PEMS試驗數據的故障診斷準確率達到96.99%,同時在不同的訓練集、驗證集組合中都能夠保持99%以上的故障診斷準確率,說明模型具有較高的泛化性能.

(3) 對比PNN及LSTM神經網絡等其他類型模型,基于時間窗口與SVM的NOx排放故障診斷模型具有更高的診斷準確率和泛化性能.

猜你喜歡
故障診斷模型
一半模型
凍干機常見故障診斷與維修
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
基于量子萬有引力搜索的SVM自駕故障診斷
3D打印中的模型分割與打包
因果圖定性分析法及其在故障診斷中的應用
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
基于LCD和排列熵的滾動軸承故障診斷
基于WPD-HHT的滾動軸承故障診斷
機械與電子(2014年1期)2014-02-28 02:07:31
主站蜘蛛池模板: 欧美色香蕉| 久久超级碰| 日韩第一页在线| 2020国产精品视频| 97人人模人人爽人人喊小说| 亚洲美女操| 亚洲欧美天堂网| 亚洲v日韩v欧美在线观看| 亚洲bt欧美bt精品| 欧美日韩亚洲国产| 国产jizzjizz视频| 欧美黑人欧美精品刺激| 国产老女人精品免费视频| 国产麻豆精品手机在线观看| 日韩欧美中文字幕在线精品| 国产视频只有无码精品| 朝桐光一区二区| 国产在线观看91精品| 国产一区二区三区夜色| 国产视频大全| 大香网伊人久久综合网2020| 亚洲不卡网| 亚洲国产成人在线| 欧美中文字幕无线码视频| 国产精品吹潮在线观看中文| 精品国产中文一级毛片在线看| 亚洲av无码牛牛影视在线二区| 国产精品综合久久久| 91黄视频在线观看| 久久亚洲国产最新网站| 日韩A级毛片一区二区三区| 欧美不卡视频一区发布| 伊大人香蕉久久网欧美| 中日韩欧亚无码视频| 国产一区亚洲一区| 欧美日韩精品综合在线一区| 精品小视频在线观看| 欧美成人国产| 国产精品香蕉| 欧美日韩第三页| 欧美一区二区三区欧美日韩亚洲| 亚洲日本中文字幕天堂网| 亚洲成人一区二区三区| 日韩午夜片| 亚洲日本www| 国产97视频在线| 午夜福利网址| 97国产精品视频自在拍| 国产丝袜精品| 91久久国产成人免费观看| 97国产一区二区精品久久呦| 毛片一区二区在线看| 日本亚洲欧美在线| 99尹人香蕉国产免费天天拍| 精品久久综合1区2区3区激情| 久草国产在线观看| 青青草国产免费国产| 国内毛片视频| 日韩大片免费观看视频播放| 亚洲福利一区二区三区| 久久亚洲高清国产| 亚洲精品成人福利在线电影| 狠狠ⅴ日韩v欧美v天堂| 无码一区二区三区视频在线播放| 九色视频一区| 精品成人免费自拍视频| 亚洲天堂网视频| 热这里只有精品国产热门精品| 日韩不卡免费视频| 色婷婷成人网| 欧美亚洲国产视频| 亚洲中文字幕在线精品一区| 国产精品高清国产三级囯产AV| 久久99国产综合精品1| 国产欧美视频在线| 国产真实乱了在线播放| 亚洲最黄视频| 波多野结衣一区二区三区四区视频 | 国产91蝌蚪窝| 久久性妇女精品免费| 亚洲成人高清在线观看| 国产成人精品男人的天堂下载|