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

基于模態分解與特征匹配的串聯故障電弧識別方法研究

2021-12-01 07:41:16張桂青曹建榮張漢元
計算機測量與控制 2021年11期
關鍵詞:模態特征故障

陳 浩,閻 俏,2,張桂青,2,曹建榮,張漢元,2,莊 園,任 飛,田 豐

(1.山東建筑大學 信息與電氣工程學院,濟南 2501012.山東省智能建筑技術重點實驗室,濟南 250101)

0 引言

國家消防救援局發布的2020年火災出警數據顯示,電氣火災數量多年來依舊持續龐大,占各種火災數量的1/3以上,大型火災中有一半以上系電氣原因引發。常見的電氣火災主要由過載、過流、電弧、漏電等原因引起[1],相比其他火災誘因,串聯故障電弧隱蔽性強,不會引起斷路器或熔斷器動作,因此,高效可靠的實現對串聯故障電弧進行識別和監測至關重要。

近些年,學者們采用多種方法致力于串聯電弧故障的識別。文獻[2]以正交二次樣條小波為工具,通過檢測電流周期性奇異點間接檢測串聯電弧故障;文獻[3]基于小波變換和支持向量機,提取特征進行分類;文獻[4]利用SWT算法提取波形動態特征,輸入SVM進行決策判斷;文獻[5]利用馬爾拉特算法,將時頻信號分解為高低頻兩部分,將高頻系數均值及其方差作為判據。以上論文在一定程度上都實現了對串聯故障電弧的識別,但漏報誤報率較高,識別精確度有待提高。

互補集合經驗模態分解(CEEMD)方法可以使復雜信號自適應的分解為多個本征模態函數( IMF,intrinsic mode function)與余差的組合,其IMF分量包含了原信號的不同時間尺度的局部特征,基函數是由數據本身所分解得到,對非平穩、非線性數據的處理優勢明顯,同時又解決了經驗模態分解方法(EMD,empirical mode decomposition)的端點效應和模態混疊問題,因此近年來在不同工程領域得到了眾多應用。趙雪花等將CEEMD與廣義回歸神經網絡組合對汾河上游月徑流量進行了分析[6],深入分解挖掘了時間序列信號的模態分量特征,比較準確地獲得了月徑流預測值,應用于有限水資源的分配領域;秦慶山將模態分解方法應用于配電網單相接地故障選線[7],將各線路首端零序電流分解后獲取最高頻分量,解決了電力系統單相接地故障選線問題;韓博躍等提出了一種基于改進IWT-CEEMD-cICA的故障診斷方法[8],通過峭度值配合選取有效函數,與其他方法的配合,提取目標振動信號,應用于滾動軸承故障識別性。

為實現建筑內串聯故障電弧的快速可靠識別與監測,本文對電氣線路高頻電氣參數進行分析,采用基于模態分解與特征匹配的識別方法,在識別結果和識別速度上較常規方法均有所提高,同時分步式的識別方法適用于云邊協同式的串聯電弧監測系統的實際部署。

1 基于模態分解與特征匹配的串聯故障電弧識別算法

1.1 算法整體框架

串聯故障電弧是在電場作用下擊穿空氣介質激烈放電的燃燒過程,同時伴隨大量熱量的釋放。當配電線路上出現長期載流量超限、線路疏于檢修、絕緣老化、接線端子連接不良等情況時,極有可能出現持續性放弧[9],串聯電弧型電氣火災可等效為在原有線路中串聯一個不規則動態變化的電弧電阻,一般不會使原線路上的電流發生激烈變化,線路上的斷路器、熔斷器等傳統保護設備無法有效發揮作用[10]。長時間的小弧燃燒也會使局部接觸電阻增大,埋下電氣火災隱患。經模擬實驗探究,串聯電弧故障高頻周期電流受線路負載的影響極大,使得識別的難度增大,對監測方法和識別算法提出了較高要求。本文將電氣參數監測器采集到的高頻周期電流差值信號進行預處理后,輸入識別算法,算法主要分為以下4個部分。

1)數據預處理。根據電弧電流特征,對采集的數據進行巴氏濾波[11],后進行歸一化處理,能夠有效排除高頻干擾信號對電弧識別的干擾,避免濾波過程中特征量損失,同時避免因采集數據值過大而對細小電弧識別產生影響,提高識別的準確度[12]。

2)一階特征提取。一階提取為粗提取,以相對寬泛的閾值將可能的串聯故障電弧納入二階提取范圍,同時確定負載類型。

高頻電氣參數監測器實時采集被測線路上的電氣參數,基于事件觸發機制實現數據的上報處理,一階提取識別機制閾值1為電流有效值變化量,閾值2為電流諧波中各次諧波占比。閾值1觸發后,將上報3個周波中的第1周波數據按識別特征通過高斯函數進行數據擬合,擬合產生3種輸出結果:阻性、阻感、非線性,確定負載類型。閾值1和閾值2均通過后進行二階特征提取識別。

3)二階特征提取。一階特征識別通過后,持續捕捉上報此后連續多個周波電氣參數,進行二階特征提取識別。二階特征提取識別即對一階提取后的疑似異常數據的周期電流差值信號進行CEEMD分解,獲取IMF實際值,IMF作為二階特征用于算法下一步的特征匹配。

4)特征匹配。利用實驗獲取的正常數據與故障數據,構建差值信號送入CEEMD算法,獲取IMF理論值作為匹配模板,與步驟3)中IMF實際值進行特征匹配運算[15]。特征匹配選用相似距離函數動態時間規整算法(DTW,dynamic time warping),DTW基于動態規劃的思想,以兩時間序列相似度作為輸出結果,在本算法中代表了電弧識別的結果,即當DTW運算結果符合相似度判定條件時,認定當前線路存在串聯故障電弧。

實際監測中,一階特征識別提取通過后,持續向云平臺服務器上報數據,進行二階特征識別提取,直至一階特征提取不滿足,停止上報。即一階特征的提取識別在本地高頻電氣參數采集器中實現,二階特征的提取識別是在遠端服務器實現。由于串聯電弧故障持續時間存在較強隨機性,獲取的串聯電弧故障3周波包含的電弧信息長度同樣具有隨機性,所以在CEEMD分解的基礎上應用DTW方法,可以有效利用其動態延展特點,解決串聯故障電弧隨機出現在上半周波或下半周波對識別結果的影響,同時解決兩不同長度時間序列信號相似性的求解問題。識別算法整體框架流程圖如圖1所示。

圖1 算法整體框架流程圖

1.2 模態分解與特征匹配方法

串聯故障電弧可以等效為原有電氣設備與故障電弧特征的疊加,如果將串聯電弧故障發生時的周期電流和正常工作電流分別進行分解,提取串聯電弧故障特征,就能以此為依據,訓練電氣火災監測算法。在正常和故障兩個維度下分別輸入識別算法,構建故障向量模型,實際應用時,通過合理設置時間窗格,將窗格前后數據差向量輸入檢測算法,通過特征匹配方法實現串聯電弧故障識別。

常用的模態分解方法有傅里葉變換、小波變換以及希爾伯特-黃變換等[17]。傅里葉變換產生的分解結果會忽視原有時間序列中的先后順序,且其窗函數和時間窗尺度的選取會對分解結果產生較大影響。小波分解的基函數和分解層數對最終結果的影響較大,其選取具有很大的盲目性。希爾伯特-黃變換中EMD分解是一種自適應分解,分解產生的每一個模態分量都代表了一種特殊振蕩模式,其分解層數由信號細節特征決定[14],不需要人為選定基函數或給出過多參數設置,但經常存在端點效應和模態混疊現象。CEEMD分解則是對EMD分解的優化和改進,通過引入一定條件的白噪聲,通過多次集成平均的形式,使剩余噪聲維持在一定范圍內,同時通過端點延拓的方法,解決端點處EMD單向分解結果不理想和一種振蕩模式下混有一種以上振蕩信號的問題。

本文對CEEMD的完全自適應分解做了限制,加入迭代判斷后,自適應分解層數可控,將分解目標由最優解改變為快速確定解,取消余差判斷與運算輸出,將最后一層余差合并入分解結果,減少分解步驟,提高分解速度,同時改變了自適應分解的對象,對觸發機制前后的向量差進行分析,不依賴對原有信號或當前信號進行分析,使算法數據處理量更小,更適應快速反應的串聯電弧故障識別。

對周期電流CEEMD分解的流程圖如圖2所示[6]。

圖2 周期電流CEEMD分解流程圖

其中,IMF1的計算方法為:

(4)

R1(t)的計算方法為:

R1(t)=I(t)-IMF1

(5)

IMF2的計算方法為:

(6)

Rk(t)的計算方法為:

Rk(t)=Rk-1(t)-IMFk

(7)

IMFk+1的計算方法為:

(8)

Ek即為算子期望,ωi(t)為均值是零的典型高斯白噪聲信號,εk為系數允許在每個階段選擇信噪比,I(t)為原高頻周期電流序列。

動態時間規整(DTW)起初多被用于語音識別領域,DTW最顯著的優勢就是不受時間序列長度、頻率等條件的約束,通過動態匹配最優路徑,以距離作為衡量標準表征任意兩個時序信號的相似性。本文以IMF理論值與IMF實際值分別為橫縱坐標,構建n*m維矩陣網格,目的是為了把兩個時間序列進行延伸和縮短,產生部分理論對應點,使兩時間序列點點對應,通過依次計算數據點與理論對應點之間的距離,來尋找兩時間序列累加距離的最小值。

本文將IMF理論值作為S1:

S1=[x1,x2,x3,…,xm]

(9)

將IMF實際值作為S2:

S2=[y1,y2,y3,…,yn]

(10)

以上述兩時間序列對應元素組成距離代價矩陣D:

(11)

其中:d(xm,yn)為此兩序列的歐氏距離:

(12)

兩時間序列間的最優路徑ξ(S1,S2):

ξ(S1,S2)=D(xm,yn)+min[ξ(x(m-1),xn),

ξ(x(m-1),x(n-1) ),ξ(xm,x(n-1) ) ]

(13)

其中:

(14)

(15)

兩時間序列“相似度”越高,ξ(S1,S2)的數值越接近零,反之越大。

2 實驗驗證

2.1 總體方案

根據《電氣火災監控系統第4部分:故障電弧探測器》GB14287.4搭建交流串聯電弧故障實驗臺[7]。實驗臺由模擬電弧發生器、實驗電路和高頻電氣參數監測器組成,如圖3所示。模擬電弧發生器由銅電極與石墨電極組成,兩電極間隔可通過軌道螺旋調節,最小可使兩電極緊密接觸,保持穩定連接,最大可使電路完全斷開,處于絕對斷路狀態。

圖3 串聯電弧模擬實驗臺及實驗臺電路

2.2 數據采集與分析

采集正常運行數據時,將模擬電弧發生器兩電極緊密接觸,使之處于完全通路狀態,電氣參數監測器實時監測線路上的電氣參數,斷開電路,更換負載,重復實驗操作。

采集串聯電弧數據時,將模擬電弧發生器兩電極分開,保持較小空氣間隙,電氣參數監測器實時監測線路上的電氣參數,接通電路,可在兩電極間觀察到明顯的電弧聲光現象,斷開電路,更換負載,重復實驗操作。

考慮到實際建筑中電路性質以線性(阻性)和非線性(阻感性、阻感容性)狀態出現,數據采集過程中選用額定功率為800 W的電熱取暖器作為典型阻性負載、額定功率為220 W的風機作為典型阻感性負載,額定功率為60 W的LED白光燈(含芯片控制電路)作為非線性負載,分別模擬了正常工作和串聯電弧故障狀態下的運行狀況[19],圖4為經巴氏濾波處理后的高頻周期電流波形。

圖4 高頻周期電流波形圖

1)阻性負載:正常工作狀態下,由于阻性負載的線性特性,其高頻周期電流波形呈現220 V交流電源的標準正弦波。發生串聯電弧故障時,在正弦波近零點處當兩電極間電壓不足以擊穿中間空氣介質時,實驗電路呈現短期斷路狀態,高頻周期電流瞬間變零,隨著電壓的正弦變化,當兩電極間電壓達到臨界起弧值時(臨界起弧值與負載類型、電路構造、環境等因素相關),空氣介質再次被擊穿,實驗電路恢復通路狀態,高頻周期電流值按正弦波波動,在電弧發生時間段內周期性重復出現[21]。

2)阻感負載:正常工作狀態下,由于受電感器件性能的影響,其電流相位滯后于電壓,伴隨有一定的畸變,由于不同電路的電阻和電感值比例存在差異,其阻抗角略有不同,所以其波形變化與阻感負載的具體數值有關,但波形趨勢一致。發生串聯電弧故障時,阻感負載呈現出與阻性負載類似的電流平肩部特征,其發生原因一致,此處不再贅述,但由于正常工作狀態時存在的畸變現象,發生串聯電弧時的畸變仍然存在。

3)非線性負載:正常工作狀態下,由于內部電子器件阻性、感性、容性的復雜組合,整體電路對外表現出典型的非線性特征,周期性明顯。發生串聯電弧故障時,由于電容的充放電效應,在平肩部階段的首和尾會出現較為明顯的反向放電過程,周期電流波形圖中表現為反向尖峰電流。

2.3 算法驗證

一階提取為雙閾值判定,由于一般建筑常規用電器功率為幾十到幾百瓦不等,正常工作電流有效值一般在0.2 A以上,故閾值1取為0.15 A,經先期實驗分析對比,閾值2諧波占比數值由模擬實驗得出。

n次諧波占比計算方法:

(1)

其中,Pn為n次諧波占比,Nn為n次諧波數值,∑13i=2Ni為除基波外2~13次諧波數值的總和(①排除基波:基波數值較高,影響高次諧波占比;②13次以上諧波基本以0形式存在,與整體識別關聯性不大)。

綜上,本實驗中閾值1為預處理后電流有效值變化量≥0.15 A,閾值2為電流諧波中各次諧波占比,詳見表1。閾值1和閾值2均通過后進行二階特征提取識別。

表1 一階特征提取諧波占比閾值表

二階提取是對一階提取通過后的高頻周期電流數據進行CEEMD分解,獲取IMF實際值,將其與IMF理論值一并送入DTW進行運算,通過合理設置DTW結果的判定條件實現特征識別。表2為二階特征提取判別條件表。

表2 電弧識別結果相似度判別條件

本文1.1節中介紹了本文采用的事件觸發機制,3個周波共上報256點數據。為方便計算,IMF理論值與IMF實際值組成的DTW計算矩陣以100為矩陣分度值。

1)阻性負載:由于阻性電路的線性關系,根據圖4可知,其上報的3周波數據具有高度的對稱性,所示實際參與特征匹配的數據點僅為3周波中任一周波的1/2半波,根據1.2節DTW輸出結果的計算辦法,完全匹配下的最優路徑即方陣對角線路徑,按照1.2節計算辦法約為9 051,按照表2轉換為百分制為0.905 1%,為保證一定的容錯機制調整為1%。

2)阻感負載:由于阻感電路的滯后特性,根據圖4可知,其上報的3周波數據均實際參與特征匹配,按照上述阻性負載的計算辦法約為5.430 5%,為保證一定的容錯機制調整為6%。

3)非線性負載:由于非線性電路的特性,按照上述阻感負載的計算辦法約為5.430 5%,因為電弧燃燒段的首尾均隨機發生反向尖峰電流,即在特征匹配過程中,最優路徑的極限不會是多個方陣對角線的疊加,經實驗測定,取其最優路徑為方陣對角線疊加的1.3~1.4倍,即7.059 65%~7.602 7%,為保證一定的容錯機制調整為8%。

但需要特別說明,表2中判別條件數值為本實驗條件下兼顧準確率和漏報率的參考值,同時為便于描述,將輸出結果轉換為百分數表示,實際應用時應在此基礎上綜合考量線路負載復雜程度做出調整。如:實際線路雖為阻感負載,但電感值較小,此時判別條件應取[1%,6%]。

1)模態分解:

分別將三種負載情況下的正常、電弧差值數據進行模態分解,分解如圖5~7所示。上述分解中,高斯白噪聲個數p取值為4,自上而下首行是高頻周期電流信號,2~5行波形對應本征模態函數IMF1~4。

圖5 阻性負載高頻周期電流分解圖

圖6 阻感負載高頻周期電流分解圖

圖7 非線性負載高頻周期電流分解圖

2)特征匹配:

為驗證本文方法的有效性,在2.3的基礎上,將IMF理論值與IMF實際值共同送入DTW運算,計算后對輸出結果進行了數值統計,散點分布圖如圖8所示。

圖8 DTW結果散點分布圖

2.4 結果分析

2.3節中,我們按照本文方法對3類典型負載正常工作與串聯故障電弧分別進行一階、二階特征提取識別,并給出了此過程的關鍵分解圖與最終識別結果的散點分布圖。

通過最終的識別結果散點分布圖可以發現,在某一類負載識別區內,本文方法的DTW輸出結果具有明顯的分層現象,上部矩形識別框和下層圓形識別框分別為非電弧狀態與電弧狀態,未出現明顯混疊現象。

1)縱向對比:

每一類負載識別區內,識別結果中以表2中設定的閾值為界,均出現兩層分層現象。阻性非電弧、阻感非電弧、非線性非電弧(3類正常負載)DTW判別條件聚堆識別結果最小值均分別大于1%、6%、8%,阻性電弧、阻感電弧、非線性電弧(3類故障負載)DTW輸出值聚堆識別結果最大值均分別小于1%、6%、8%。

2)橫向對比:

每一類負載識別區內,串聯故障電弧DTW判別條件分層聚堆最靠下,即識別結果值最小,根據DTW的計算原理,即相似度最高。非電弧數據由于運行狀態的不同,呈現各自分散聚堆現象,但均未發生串聯電弧故障,DTW輸出結果均大于電弧數據。因此本文使用的方法及設置的閾值對串聯故障電弧的識別真實有效。

本文方法閾值分級設置且具有容錯機制,一定程度上限制了散點的分散程度,但仍有部分散點具有離散趨勢,主要與放弧強度有關,模擬電弧的釋放過程伴隨著眾多外界因素的干擾,中途可能出現短暫的熄弧與復燃過程,此部分散點是影響識別結果準確率的主要因素。

2.5 與其他識別方法的比較

為驗證提出方法對串聯電弧故障識別的有效性,本文選取了GLGCO-SVM[22]與TDV-CNN[23]進行了對比分析,兩篇文獻中使用的典型負載類型均與本文相似。GLGCO-SVM的卷積神經網絡相對復雜,尤其是在特征提取方面,它以稀疏編碼的方式獲取關鍵信息增加了算法的復雜度,運行時間延長;TDV-CNN是基于灰度圖像處理的方法,當發生串聯電弧故障時線路上會產生大量高頻信號,對信號進行合理濾波后,只保留電弧特征信號,但其濾波僅保留電弧特征信號的準確度難以保證[24]。

本文從時效性和準確性指標進行方法對比。時效性指標采用算法運行時間作為評價標準,準確性指標采用準確率作為評價標準。本文使用的計算機設備Intel(R) Core(TM) i5-6200U CPU @ 2.30 GHz 2.40 GHz,機帶RAM4.00 GB,64位操作系統,Windows 10 家庭版,版本號20H2。

1)時效性:

220 V、2~10 A交流電下,串聯故障電弧短時間內就能產生2 000~4 000 ℃的高溫,超過0.5 A就能夠拉弧[25]。因此,串聯故障電弧的識別對時效性要求很高。本文將10組串聯故障電弧數據送入算法,算法運行時間統計表如表3所示。

表3 算法運行時間統計表

根據表3的時間統計,識別方法可以在0.3 s左右的時間內實現故障電弧的識別,滿足《電氣火災監控系統第4部分:故障電弧探測器》GB 14287.4的規定要求。

2)準確性:

串聯故障電弧的識別不僅要求反應快速,更要識別準確。將實驗獲取的60組串聯電弧故障數據隨機混入600組正常數據中,組成測試集,分別送入識別算法。串聯電弧故障識別準確率如表4所示。

表4 串聯電弧故障識別結果準確率

根據表4的統計,本文從三種負載類型識別結果準確率及均值四個維度進行對比,針對各類型串弧故障的識別準確率代表了單體識別的有效性,識別結果準確率均值代表了對不同負載類型普適性。

CEEMD-DTW方法對隨機產生的故障數據均表現出較好的適應性,對阻性串弧、阻感串弧、非線性串弧都具有91%以上的識別準確度。

在相同條件下,將本文方法與GLGCO-SVM方法對比,二者識別準確率均值類似,但前者針對阻性負載串聯故障電弧的識別優勢明顯,且算法結構清晰簡單,可在嵌入式設備中實現,可有效實現云邊協同識別;與TDV-CNN方法進行對比,本文方法針對非線性負載串聯故障電弧的識別優勢明顯,且識別結果準確率均值也高于對比方法。由于TDV-CNN方法在池化階段會丟失大量的數據特征,對隨機性較強、周期性較差的數據進行特征匹配時識別結果并不理想,因此其針對非線性串弧表現出的隨機性和非周期性特征識別效果較差,導致識別準確率均值較低。

3 結束語

1)本文依托監測線路高頻電氣參數提供了一種串聯電弧故障識別方案,不易受到外界溫度、煙霧濃度等環境因素的影響,滿足時效性要求,可信度較高。

2)本文針對串聯電弧故障的識別要求,采用事件觸發機制,將事件前后向量差作為識別對象,通過高斯擬合方法確定線路負載類型,巧妙設置了一階、二階特征提取識別閾值,利用CEEMD-DTW算法進行串聯電弧故障識別,普適性較強。

3)本文主要以某一類型負載的串聯電弧故障為研究對象。但實際中,除大型或重要單一設備外,線路監測對象復雜組合情況下,其非線性特征可能會有變化,實現準確的串聯電弧故障識別仍然是這一領域的難點。本文方法應用于復雜組合負載情況時有待進一步分析和研究。

猜你喜歡
模態特征故障
故障一點通
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
奔馳R320車ABS、ESP故障燈異常點亮
國內多模態教學研究回顧與展望
故障一點通
基于HHT和Prony算法的電力系統低頻振蕩模態識別
江淮車故障3例
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 尤物国产在线| 国产精品天干天干在线观看| 波多野结衣一区二区三区四区| 国产精品私拍在线爆乳| 91精品免费高清在线| 一级毛片免费不卡在线| 精品成人免费自拍视频| 在线欧美国产| 亚洲综合第一页| 成人一级黄色毛片| 91久久精品国产| 国产另类视频| 国产精品视频999| 99久视频| 国产成年女人特黄特色大片免费| 亚洲精品视频免费看| 久久黄色影院| 试看120秒男女啪啪免费| 女人av社区男人的天堂| 久久国产精品电影| 国产亚洲日韩av在线| 另类综合视频| 国产精品自拍露脸视频| 91精品免费高清在线| 日本手机在线视频| 日韩在线播放中文字幕| 成人国产精品一级毛片天堂| 色综合天天娱乐综合网| 无码在线激情片| 欧美日韩成人在线观看| 国产在线自揄拍揄视频网站| 亚洲精品无码在线播放网站| 亚洲第一视频免费在线| 日韩精品久久久久久久电影蜜臀| 亚洲伊人久久精品影院| 精品国产污污免费网站| 久久夜色精品| 欧美日韩免费在线视频| 亚洲一区无码在线| 亚洲激情99| 玖玖精品视频在线观看| 视频二区亚洲精品| 亚洲色图欧美激情| 国产无码网站在线观看| 综合色天天| 色综合成人| 精品视频免费在线| 午夜福利在线观看成人| 19国产精品麻豆免费观看| 亚洲资源站av无码网址| 99免费视频观看| 在线另类稀缺国产呦| 成人午夜亚洲影视在线观看| 国产黄在线观看| 亚洲天堂2014| 亚洲日韩AV无码一区二区三区人| 精品视频福利| 中文字幕有乳无码| 免费a级毛片视频| 国产成人精品免费视频大全五级| 日韩 欧美 小说 综合网 另类| 国产无码制服丝袜| 青青热久免费精品视频6| 亚洲香蕉久久| 成人一级黄色毛片| 日韩人妻无码制服丝袜视频| 亚洲天堂精品视频| 精品国产欧美精品v| 波多野结衣爽到高潮漏水大喷| av在线手机播放| 黄色网址免费在线| 国产91全国探花系列在线播放| 国产成人亚洲欧美激情| 四虎成人在线视频| 国产一区二区三区精品久久呦| 欧美亚洲第一页| 亚洲中文无码h在线观看 | 欧美精品啪啪| 天天综合网站| 91精品国产情侣高潮露脸| 制服丝袜一区二区三区在线| 国产高清在线观看91精品|