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

高速列車牽引傳動系統仿真可信度評估方法

2020-04-16 13:16:58彭小奇陽春華
鐵道學報 2020年1期
關鍵詞:故障信號

李 晟 彭小奇 彭 濤 陽春華

(1.中南大學 自動化學院,湖南 長沙 410083;2.江西理工大學 理學院,江西 贛州 341000;3.湖南第一師范學院 信息科學與工程學院,湖南 長沙 410205)

我國高速鐵路發展迅速,形成了具有中國特色的高鐵技術體系。牽引傳動系統是高速列車的關鍵系統[1],其能否穩定可靠工作是確保高速列車運行安全的關鍵。目前高速列車牽引傳動系統的應用驗證多采取虛擬仿真或半實物仿真方法,其結果能否正確反映實際系統的性能,取決于仿真系統的可信度,未經可信度評估的仿真系統沒有任何價值[2]。但長久以來,仿真可信度問題并沒有獲得足夠的關注[3-4]。

近年來,許多學者對高速列車牽引傳動系統仿真進行了研究。西南交通大學、北京交通大學的學者針對CRH 高速列車牽引傳動系統,利用MATLAB 搭建虛擬仿真系統或基于dSPACE 搭建半實物仿真平臺,通過仿真著重分析牽引時傳動電機電磁轉矩、三相電流、直流環節支撐電容電壓和列車行駛速度的變化情況[5-7]。文獻[8]以我國典型高速列車CRH2型動車組為例,通過仿真模型分析牽引傳動系統在不同線路條件下的運行情況。文獻[9]提出了基于HLARTI和反射內存網的混合網絡系統構架的CRH2型高速列車牽引傳動系統仿真平臺。但目前的文獻多數未對其搭建的仿真系統的可信度進行討論,或僅通過對仿真波形和實測波形進行目測對比來定性分析,而沒有給出客觀的評估結果。文獻[10]提出了一種改進的模糊層次分析法(Fuzzy AHP,FAHP),針對列車運行控制系統仿真進行可信度評估,該方法考慮了評估者對復雜系統判斷的模糊性,給出了可信度指標權重確定方法。FAHP能更好地克服人為主觀判別、偏頗對決策結果的不利影響,但并未給出可信度指標的具體量化方法且并不適用于牽引傳動系統仿真。由文獻檢索可知,在高速列車仿真領域目前還缺乏針對牽引傳動系統仿真的可信度評估理論體系。如何給出適用于高速列車牽引傳動系統仿真信號特點的可信度指標量化模型,如何利用評估專家的寶貴經驗確定可信度指標的權重并盡可能降低主觀偏好帶來的不利影響,是急待解決的關鍵問題。

本文針對牽引工況下,在理想情況即平直道路上運行時的高速列車牽引傳動仿真系統,提出一種新的仿真可信度評估方法。該方法首先根據信號的緩變特性,利用殘差相似度構造指標的整體可信度;再根據信號的瞬變特性,利用基于SVD-TLS算法改進的Prony分析采集信號的頻率、幅值和阻尼特征量,根據相似性原理利用各特征量構造指標的特征可信度;通過綜合整體可信度和特征可信度給出可信度指標的量化模型。其次,使用群組AHP 方法確定可信度指標的權重,有效降低了主觀偏好的不利影響。最終,給出仿真可信度的綜合量化模型,形成一套完整的仿真可信度理論體系。將此方法應用到CRH2型動車組牽引傳動系統的無故障仿真和故障注入仿真實驗的可信度評估過程中,準確客觀地評估了仿真結果的可信度;并將可信度評估結果用于仿真模型關鍵參數的修正,為參數修正提供了客觀的數據支持。實踐表明,本文所提方法是有效、合理的。

1 高速列車牽引傳動仿真系統簡介

高速列車牽引傳動系統一般由單相三電平脈沖整流器、直流環節支撐電容、三電平逆變器和牽引電機組成的牽引傳動系統主電路和牽引控制器組成。牽引控制器接收車載傳感器反饋信號,按給定的控制策略輸出控制指令,控制脈沖整流器、逆變器等,總體思路見圖1,仿真所涉及的電氣設備電氣參數均參考文獻[1]和文獻[6]中的數據設置。

圖1 基于Simulink牽引傳統系統虛擬離線故障注入仿真平臺總體思路

基于Simulink的牽引傳動系統故障注入/模擬Benchmark框架見圖2。由故障注入基準庫、故障注入模塊和故障注入控制器組成,采用修改信號和變量/參數的方式,完成牽引傳動控制系統各故障模式的注入、仿真/模擬。其中,故障注入基準庫提供所有故障的屬性信息(包括故障模式及故障數學描述等);故障注入模塊是根據已實現子系統中各故障類型的模擬方案,在Simulink環境下搭建出的故障注入模塊庫,包括變流器、控制器TCU、傳感器、牽引電機4個部分的故障注入模塊;故障注入控制器用于控制及實現具體故障信號的注入、仿真/模擬。

根據實際車載傳感器的分布,通常選擇電磁轉矩、三相電流、直流環節支撐電容電壓和列車行駛速度等作為反映仿真可信度的觀測信號[6]。這些信號各具特點,但均包含了緩變和瞬變信號的特征,尤其在故障注入時,瞬變信號的特征顯得尤為重要。因此在進行仿真可信度評估時,應同時考慮緩變信號和瞬變信號的特征提取。

2 仿真可信度評估方法

2.1 可信度指標的量化模型

2.1.1 基于殘差分析的整體可信度

殘差時序向量是將仿真信號時序向量與實測信號時序向量相減后得到的,它很好地反映了信號的整體誤差和緩變信號特征,據此可建立可信度指標的量化模型

圖2 基于Simulink的牽引傳動系統故障注入/模擬Benchmark框架

式中:xi、分別為實測信號、仿真信號幅值序列的第i個幅值;n是幅值序列的維數;σi為第i個點的可信度;αi為第i個點的重要程度;CE為仿真信號與實測信號的殘差可信度量化指標,稱為殘差相似度。易證明,σi關于單調遞減,因此CE也滿足單調性。殘差相似度作為2個時間序列間的距離,是一種靜態度量,體現了信號的整體誤差水平,因此定義CE為可信度指標的整體可信度。由于僅通過CE難以區分信號間不同瞬變性能的誤差,因此還需對信號誤差進行進一步分析,以提取信號的瞬變特征。

2.1.2 基于特征分析的特征可信度

特征值分析法的思路是比較仿真信號和實測信號的各特征值間的差異,通過這些差異來計算信號間的誤差。Prony分析可以提取信號的頻率、幅值、阻尼和能量等特征,它是用一組指數函數的線性組合來擬合等間隔采樣數據。Prony分析提取的特征量可以很好地表征信號的瞬變特性,故被廣泛應用于低頻振蕩分析、模態識別和雷達信號分析等領域[11-13]。

本文選取信號頻率、幅值和阻尼3個特征,構建頻率相似度、幅值相似度和阻尼相似度3個參量來表征信號瞬態特征的誤差。先給出信號的能量表達式,設n為采樣點數為采樣點幅值,則離散信號的能量公式為

通過Prony分析,提取實測信號和仿真信號的頻率、能量向量對分別為

式中:FM、FS分別為實測信號和仿真信號的頻率、能量向量對;fMi、fSi分別為實測信號和仿真信號頻率向量的第i個頻率分量;WMi、WSi分別為實測信號和仿真信號能量向量的第i個能量分量。通常n1≠n2,且n1和n2均較大,故采用以下方法降維并統一維數。方法步驟如下:

Step1將能量小于總能量0.1%的信號分量刪除,并按照能量從小到大的順序將剩余分量重新排序,獲得新的頻率、能量向量對

實驗表明:n3、n4遠小于n1、n2。f′Mi和f′Si分別為降維并統一維數后實測信號和仿真信號頻率向量的第i個頻率分量;W′Mi和W′Si分別為降維并統一維數后實測信號和仿真信號能量向量的第i個能量分量。

Step2生成新的頻率、能量向量對F″S。其生成規則是:F″S中第i個元素為F′S與F′M中第i個元素相距最近的元素。距離計算式為

式中:dij為F′S中第i個元素到F′M中第j個元素的距離;χ為可以根據具體信號調整的常數;a′Mi、a′Sj為對應頻率的幅值。

通過以上步驟,向量F″S與F′M維數一致,頻率相似度為式中:f′Mi、f″Si分別為實測信號與仿真信號中第i個點的頻率;σfi為第i個點的可信度;W′i表示降維并統一維數后第i個點的能量;βi為第i個點的權重;CF為頻率相似度。因非周期信號頻率為零,故在頻率相似度的計算中不予考慮。

類似可以給出幅值相似度的定義。在得到維數一致的2個幅值、能量向量對A″S、A′M后,幅值相似度為

式中:a′Mi、a″Si分別為實測信號與仿真信號中第i個點的幅值;σAi為第i個點的可信度;γi為第i個點的權重,計算式同式(10);CA為幅值相似度。

同理,得到維數一致的2個阻尼、能量向量對D″S、D′M,相應的阻尼相似度計算式為

式中:d′Mi、d″Si分別為實測信號與仿真信號中第i個點的阻尼因子;d′Mid″Si<0為信號衰減趨勢完全相反,故相似度為0;σdi為第i個點的相似度,其計算式同式(10);δi為第i個點的權重;CD為阻尼相似度。幅值相似度和阻尼相似度的計算均應考慮非周期分量信號。

由于頻率、幅值和阻尼相似度均源自于信號特征,可以認為3者具有同等重要性,即其權重應相同,故定義表征信號瞬變特性的特征可信度CT為

綜合考慮整體可信度和特征可信度所構建的仿真可信度評估指標的可信度量化模型為

式中:0≤λ≤1。λ越大,表示越重視信號緩變特征對可信度指標量化的影響;λ越小,表示越重視信號瞬變特征對可信度指標量化的影響。

2.2 基于群組AHP的指標權重確定方法

在仿真可信度評估過程中,僅依靠誤差分析,會丟失專家經驗中的有用信息,但評估專家的判斷又受其學術水平、工作經歷以及對待評對象的熟諳程度等主觀因素的影響。為此,本文引入群組AHP 方法來確定各可信度評估指標的權重,以便既吸收專家的寶貴經驗,又使評估結果更加客觀合理。

群組AHP 方法(Group AHP)[14]是由多位專家運用AHP法對待評對象進行分析,然后再對各專家的判斷矩陣進行聚合,得到綜合判斷矩陣,從而形成更加客觀合理的群組評估結果。其核心問題是專家權重的確定和綜合判斷矩陣的一致性分析。

邀請m位待評仿真系統相關領域的專家組成評估專家集合E={e1,e2,…,em},第k位專家的權重用ηk表示。針對待評仿真系統,每位專家按照AHP的1-9標度設置方法,給出相應的判斷矩陣Ak=(akij)l×l(k=1,2,…,m),其中l表示AHP方法中同層可信度指標數。對所有判斷矩陣進行一致性檢驗,使一致性比例CR小于0.1,經專家反復調整,直到所有判斷矩陣都通過一致性檢驗為止[15]。CR的表達式為

式中:n為判斷矩陣的階數;λmax為判斷矩陣的最大特征值;RI是可通過查表得到的平均隨機一致性指標。

判斷矩陣的一致性水平是專家評判可靠性和嚴謹程度的體現。CR越小,判斷矩陣的可靠程度越高,相應專家的判斷重要程度越高;CR越大,表明判斷矩陣的可靠性低,相應評估專家的判斷重要程度降低。因此評估專家權重為

由于CRk<0.1才能通過一致性檢驗,故μ取10較為合適。

群組AHP方法的另一個關鍵問題是將所有評估專家的判斷矩陣進行集結,得到包含所有專家判斷信息的綜合判斷矩陣,該矩陣也要通過一致性檢驗。文獻[16]已證明,若A1,A2,…,Am為m個專家給出的判斷矩陣且均為一致性可接受判斷矩陣,則其加權幾何平均綜合判斷矩陣也是一致性可接受矩陣。由此定義綜合判斷矩陣為

式中:?為加權幾何平均集結算子。

設l為同層可信度指標數,則上述矩陣元素的計算式為

可以求得歸一化的可信度評估指標的綜合權重向量ω=[ω1,ω2,…,ωl]。其中,ωi表示該層第i個可信度指標的綜合權重。

2.3 仿真系統可信度評估方法

本文提出的仿真系統可信度評估方法見圖3。

圖3 仿真系統可信度評估方法體系框

設Ci是根據式(18)獲取的該層第i個可信度量化值,則基于信號誤差分析和群組AHP 方法的仿真可信度評估量化模型為

3 仿真實驗與可信度評估算例

由于無法獲取CRH2型動車組牽引傳動控制系統實際信號的錄波數據,本文利用dSPACE 仿真器和實物控制器搭建圖4所示的系統半實物實時仿真平臺進行半實物仿真,并以基于該平臺的仿真結果作為參照的實測數據。

圖4 基于dSPACE的高速列車組牽引傳統系統半實物實時仿真平臺實現總體思路

3.1 Prony分析階數確定

在利用Prony進行信號分析時,階數選擇對分析的結果有較大的影響,甚至導致所提取的特征信息嚴重失真,無法用于可信度評估。故提取特征信息前須確定合理的階數。

影響階數選擇的指標主要有:動態變化率DVR和信噪比SNR[17]。DVR反映了重構數據與實測數據的整體誤差大小,DVR越低,重構數據與實測數據越吻合。而SNR越大,表示Prony估計的精度越高。

文獻[18]指出,Prony分析在實際應用中取能使SNR達到40 dB且DVR小于1%的階數,即可滿足運行時間和精度要求。為得到滿足2個指標的階數,本文采用奇異值-總體最小二乘法(Singular Value Decomposed-Total Least Square,SVD-TLS)對原始信號進行去噪并確定Prony分析的階數[19],再使用Prony算法提取信號特征。

3.2 無故障仿真算例

分別在基于Simulink的虛擬離線仿真平臺和基于dSPACE的半實物實時仿真平臺上進行無故障仿真。采樣間隔3×10-3s,時長10 s。仿真結果見圖5。通過目測法定性對比逆變器a相電流、列車行駛速度、牽引電機電磁轉矩、直流環節上支撐電容電壓4個觀測點的信號可見,兩者的仿真結果很接近。但定性分析的結果并不客觀,還需要進行定量的客觀分析。

圖5 Simulink和dSPACE仿真結果對比

以逆變器a相電流信號的重構為例,由SVD-TLS算法確定的Prony分析階數為58時,原信號與重構信號波形見圖6,兩者幾乎重合,此時SNR=41.983 4,DVR=0.006 33,符合精度要求。可見,本文采用的基于SVD-TLS算法改進的Prony分析所得到的頻率、幅值、阻尼因子和能量向量是有效的。

采用本文所提出的相似度量化公式分別求各相似度的量化值,結果見表1。

圖6 原始信號與Prony重構信號比較

表1 無故障仿真各相似度量化值 %

由于是無故障仿真,故可以認為整體可信度和特征可信度的重要性相同。令式(18)中的λ=0.5,利用式(17)、式(18)分別計算各可信度指標的可信度量化值,可得:逆變器a相電流可信度C1=91.54%,列車行駛速度可信度C2=96.83%,牽引電機電磁轉矩可信度C3=94.34%,直流環節上支撐電容電壓可信度C4=95.29%。

邀請4位相關領域內專家組成仿真可信度評估專家集合E={e1,e2,e3,e4},每位專家按照AHP 方法將逆變器a相電流、列車行駛速度、牽引電機電磁轉矩、直流環節上支撐電容電壓4個可信度指標進行兩兩間的重要性比較,按照1-9標度法分別給出4個判斷矩陣

分別對4個判斷矩陣進行一致性檢驗,由式(19)、式(20)可以得到一致性比例CR,分別為CR1=0.005 4,CR2=0.027 3,CR3=0.044 0,CR4=0.017 8,均小于0.1,故4個判斷矩陣均為一致性可接受矩陣。取μ=10,利用式(21)可以確定專家權重η分別為η1=0.289 5,η2=0.239 7,η3=0.211 9,η4=0.259 0。

利用式(22)、式(23)計算可得綜合判斷矩陣A為

計算一致性比例CR=1.015 2×10-4,遠小于0.1,且有CR<CRi(i=1,2,3,4)。可見群組AHP方法所獲得的權重比AHP 方法獲得的權重更加客觀準確,有效降低了主觀意志的影響。

根據綜合判斷矩陣,采用特征值法可以求得可信度評估指標的綜合權重為

根據式(24)可得基于Simulink的CRH2高速列車組牽引傳動系統仿真綜合可信度為

仿真綜合可信度的分析計算結果非常接近100%,表明基于Simulink的仿真系統與作為參照的dSPACE仿真系統的仿真結果一致。

上述結果表明,本文提出的仿真可信度指標合理有效,通過群組AHP 方法確定權重能吸收評估專家的寶貴經驗,使可信度指標權重更加客觀合理;同時也證明,基于Simulink的虛擬離線仿真系統的仿真結果是可信的,該仿真系統模型是準確有效的。

3.3 實際故障仿真算例

3.3.1 故障注入仿真實驗

牽引變流器是牽引傳動系統的動力源和能量轉換核心設備,一旦發生故障,將使列車運行在不穩定狀態,嚴重時將導致列車喪失動力,被迫停駛。2007年5月28日CRH2型012A 動車組7車主斷路器斷開,復位后故障未消除,經查是牽引變流器發生故障,因此采取切除8車牽引變流器、保留75%動力維持運行。下載故障數據,并通過廠家專用分析軟件對其進行分析,確定逆變單元IGBT 模塊內部U 相上橋臂2號功率器件短路是導致該故障的真正原因。采用基于Simulink的虛擬仿真系統對這次故障進行仿真。設定故障嚴重程度因子為0.95,故障注入時間3 s,選擇中間電路上下支撐電容電壓,以及列車行駛速度和電磁磁矩的仿真結果與實際錄波進行對比(見圖7),仿真波形與實錄波形變化趨勢基本一致,局部存在差別,表明通過故障注入仿真實驗可以較為成功地復現實際故障。

圖7 Simulink故障仿真結果和實錄波形結果對比

3.3.2 仿真結果可信度分析

為客觀評估仿真結果與實際故障間的差異,并借助可信度來修正仿真模型的關鍵參數,以提高仿真的準確度,利用本文提出方法定量分析仿真結果與實際故障間的差異。當故障嚴重程度因子設置為0.95時,各可信度指標的量化值見表2。

表2 IGBT 短路故障各相似度量化值(故障嚴重程度因子0.95) %

4位評估專家給出的判斷矩陣如下

取μ=10,利用式(21)計算得到專家權重為

由式(22)、式(23)得到綜合判斷矩陣為

用特征值法求得基于群組AHP的可信度評估指標的綜合權重分別為:直流環節上支撐電容電壓ω1=0.238 5,下支撐電容電壓ω2=0.253 8,列車行駛速度ω3=0.314 9,牽引電機電磁轉矩ω4=0.192 8。由式(17)、式(18)分別計算各可信度指標的可信度綜合量化值,可得:直流環節上支撐電容電壓可信度C1=85.92%,下支撐電容電壓可信度C2=84.24%,列車行駛速度可信度C3=99.72%,牽引電機電磁轉矩可信度C4=85.57%。

因故障狀態下,瞬變特征比緩變特征更為重要,故取式(18)中的λ=0.3,最后由式(24)計算得到基于Simulink 虛擬故障注入仿真系統的綜合可信度為R0.95=89.77%。

由計算結果可見,本次故障注入仿真的可信度較高,仿真波形與實際波形吻合較好,定性分析和定量分析結果一致。可見,利用本文提出的方法得到的仿真可信度評估結果是正確有效的。

3.3.3 基于可信度的仿真模型關鍵參數修正

仿真模型關鍵參數的設置不同,會導致仿真結果間存在差異。關鍵參數設置的不合理,甚至會導致仿真結果的嚴重失真。通過定性比較來修正模型參數的難度和主觀性較大,而利用本文所提出的仿真可信度計算方法,可通過不同的關鍵參數設置計算出對應的仿真可信度量化值,將其作為修改關鍵參數的客觀依據。以故障注入仿真實驗時,故障模型中的故障嚴重程度因子這一關鍵參數的選擇為例,驗證本文方法能夠對模型關鍵參數的修正提供指導。

在IGBT 短路或開路故障仿真中,故障嚴重程度因子多采用經驗設定,主觀性很大,也無法判斷該參數設定值的優劣。為此,設置故障嚴重程度因子分別為0.95、0.90、0.85,并進行故障注入仿真。IGBT 短路故障會影響列車的平穩運行,對列車速度的影響較為明顯。通過與CRH2型動車組上速度傳感器采集的故障數據進行對比,得到列車速度的仿真波形與實際波形,見圖8。由圖8可以看到,故障嚴重程度因子分別為0.95、0.90、0.85時,仿真波形與實錄波形變化趨勢基本一致,通過定性分析無法得知哪個參數更為合理。

圖8 不同參數時Simulink故障仿真結果和實錄波形結果對比

表3、表4列出了故障嚴重程度因子為0.90和0.85時對應的相似度量化值,計算可得綜合可信度量化值分別為R0.90=90.94%,R0.85=86.90%。可見,當故障嚴重程度因子設定為0.90時,仿真綜合可信度最高,且各相似度量化值也不同程度的高于0.95和0.85。因此,選擇故障嚴重程度因子為0.90更加準確合理。

表3 IGBT 短路故障各相似度量化值(故障嚴重程度因子0.90) %

表4 IGBT 短路故障各相似度量化值(故障嚴重程度因子0.85) %

4 結論

(1)提出了一套完整的仿真可信度評估方法。提出了仿真的整體可信度的概念,通過殘差分析來構造整體可信度量化模型;提出了仿真的特征可信度的概念,通過Prony分析提取信號特征量,并通過提取的特征量構造特征可信度量化模型;將整體可信度和特征可信度相結合,得到各可信度指標的可信度量化模型;基于群組AHP 方法確定各可信度指標的權重,有效降低了確定權重過程中的主觀性,在吸收專家經驗的同時,使評估結果更加科學合理;通過SVD-TLS算法確定Prony分析的階數,使特征可信度的計算準確可靠。

(2)將虛擬離線仿真系統的仿真結果分別與基于dSPACE的半實物實時仿真系統的數據和實際發生的某次故障的數據進行對比。實例表明,利用本文所提出的方法可以計算仿真系統的綜合可信度量化值,從而客觀的對仿真結果的準確性進行評估。

(3)依據綜合可信度量化值在不同仿真模型關鍵參數設置下的計算結果,可對仿真模型的關鍵參數進行修正,進一步提高仿真的準確度。

猜你喜歡
故障信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
孩子停止長個的信號
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
故障一點通
故障一點通
故障一點通
主站蜘蛛池模板: 久久精品亚洲热综合一区二区| 99re66精品视频在线观看| 欧美成人午夜影院| 丰满的熟女一区二区三区l| 免费一极毛片| 色偷偷一区| 久久伊人色| 成人精品亚洲| 小说 亚洲 无码 精品| 最新国产麻豆aⅴ精品无| jizz国产视频| 久久一日本道色综合久久| 久久视精品| 毛片免费网址| 亚洲五月激情网| 国产在线视频二区| 国产美女精品一区二区| 国产真实二区一区在线亚洲| 欧美人与牲动交a欧美精品| 一级看片免费视频| 免费在线a视频| 久久精品国产999大香线焦| 免费国产在线精品一区| 久久国产精品无码hdav| 亚洲无码不卡网| 亚洲乱伦视频| 国产成人久视频免费| 亚洲 成人国产| 中文无码精品A∨在线观看不卡 | 亚洲女同一区二区| 99热这里只有精品国产99| 五月婷婷综合网| 经典三级久久| 欧美区一区二区三| 国产微拍精品| 亚洲h视频在线| а∨天堂一区中文字幕| 欧美精品在线免费| 亚洲一区二区三区国产精华液| 久久精品一品道久久精品| 精品国产网| 亚洲高清无在码在线无弹窗| 正在播放久久| 青青草国产精品久久久久| 日韩福利在线观看| 911亚洲精品| 色婷婷天天综合在线| 88国产经典欧美一区二区三区| 国产在线一二三区| 亚洲AV一二三区无码AV蜜桃| 欧美日本在线播放| 凹凸国产熟女精品视频| 99热这里只有精品在线观看| 欧美在线精品一区二区三区| 91在线视频福利| 国产又爽又黄无遮挡免费观看| www中文字幕在线观看| 国产日韩精品一区在线不卡| 欧美a级在线| 亚洲国产天堂久久综合| 九色在线视频导航91| 9啪在线视频| 欧美日韩成人| 国产三区二区| 中文无码毛片又爽又刺激| 91原创视频在线| 久久香蕉国产线| 国产鲁鲁视频在线观看| 91午夜福利在线观看| 亚洲色欲色欲www在线观看| 欧美性色综合网| 天天综合色天天综合网| 亚洲欧美成人在线视频| 国产传媒一区二区三区四区五区| 久久精品一品道久久精品| 午夜精品一区二区蜜桃| 亚洲欧洲日产国产无码AV| 国产一级在线观看www色 | 国产精品欧美亚洲韩国日本不卡| 成年片色大黄全免费网站久久| 黄色污网站在线观看| 在线中文字幕日韩|