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

應用優化限制帶寬經驗模態分解法識別電力變壓器繞組模態參數

2014-09-07 05:52:56周求寬王豐華萬軍彪段若晨
振動與沖擊 2014年13期
關鍵詞:模態變壓器振動

周求寬,王豐華,萬軍彪,耿 超,段若晨

(1. 江西省電力科學研究院,南昌 330096;2. 上海交通大學 電氣工程系,上海 200240)

作為電力系統中的重要組成部分之一,電力變壓器的穩定性與可靠性與整個電網的安全運行密切相關。變壓器在正常運行時,繞組會受電動力激勵而產生振動,其激勵力的頻率為100 Hz。若變壓器繞組固有頻率與激勵力頻率相近,將會產生共振現象,使得變壓器繞組的振幅劇增,進而可能引起繞組嚴重變形、絕緣損壞甚至繞組失穩坍塌等嚴重故障。因此,使繞組固有頻率遠離激勵頻率對變壓器的優化設計與制造意義重大[1]。另一方面,變壓器繞組在長期運行過程中可能會遭受到短路沖擊的作用,短路沖擊中產生的強大電動力可能使繞組發生松動或局部變形,使繞組機械特性發生明顯改變,從而導致模態參數的變化,進而表現為變壓器繞組振動特性及振動信號的改變[2];因此,變壓器繞組模態參數的準確識別與變壓器的制造與運行密切相關,具有重要的理論研究和工程應用價值。

變壓器繞組的模態參數主要包括固有頻率、振型和阻尼比等,它們恰當全面的描述了整個繞組機械系統的整體固有動態特性。目前,對于變壓器繞組這類機械結構模態參數的獲取途徑主要是對繞組結構進行試驗測試和數據獲取,進而進行分析和參數識別。傳統的模態參數識別方法主要包括最小二乘負指數法、時間序列分析法等時域辨識方法[3]和最小二乘圓擬合法、分區模態綜合法、頻域總體識別法[4]等頻域識別方法兩大類,但這些方法大都僅限于在時域或者頻域中單獨對測試信號進行識別,參數辨識精度有限;同時,已有的識別方法大都對測試信號中的噪聲較為敏感,且固有頻率的識別與定階存在較大的主觀性;因此,對電力變壓器繞組這類結構復雜的強非線性機械結構系統來說,需要尋求新的模態參數識別方法提高模態參數識別的準確性,這是變壓器繞組結構健康監測的重要前提。

本文根據某10 kV油浸式電力變壓器繞組的軸向模態實驗的測試結果,嘗試引入經驗模態分解法(Empirical Mode Decomposition, EMD)獲取變壓器繞組的模態參數。為了抑制該方法本身固有的模態混疊問題,在原始信號中加入屏蔽信號的同時,提出使用粒子群優化算法(Particle Swarm Optimization, PSO)選取最優的屏蔽頻譜,藉此進一步提高EMD算法的分解能力。文中同時使用目前常用的頻域識別方法PolyMax法對實驗變壓器繞組的模態參數進行識別。

1 理論基礎

1.1 經驗模態分解

經驗模態分解法可將多成分的復雜信號分解成一系列的單頻率成分的本征模態函數(Intrinsic Mode Function, IMF),是對以傅里葉變換為主的穩態分解方法的重要突破,適合于分析非平穩、非線性信號,由Huang等[5]提出。其分解過程為[6]:

(1) 確定原始信號x(t)的極值點,并用三次樣條函數對極大值與極小值點進行差值,形成原始信號的上下包絡線。取上下包絡線的平均值,定義為m1(t)。用原始信號減去m1(t),得到一個新的數據序列h1(t),如式(1)所示。

h1(t)=x(t)-m1(t)

(1)

(2) 判斷h1(t)是否滿足IMF的兩個條件:① 在整個時間跨度內,極值點和過零點的個數相等或至多相差1個;② 在任意一點處,上包絡線和下包絡線的均值為零。若h1(t)滿足條件,則提取出第1個IMF分量,定義為c1(t);若不滿足條件則返回第(1)步繼續篩選,直至提取出IMF1。

(3) 將原始信號x(t)減去c1(t),得到新的數據序列定義為r1(t),將r1(t)作為原始信號重復步驟(1)和步驟(2),直至rn(t)小于預定值或為單調函數時,停止篩選。則原始信號x(t)分解為n個IMF分量及余量rn(t),如式(2)所示:

(2)

1.2 變壓器繞組模態參數識別

變壓器繞組為一典型多自由度系統,其軸向振動微分方程為[7]

(3)

式中:M為質量矩陣;K為剛度矩陣;C為阻尼矩陣;F(t)為激勵力矩陣。

根據模態疊加理論可知,繞組的自由振動響應為[8]

(4)

式中:ζr為阻尼比;ωr為無阻尼固有頻率;ωdr為有阻尼固有頻率,且在阻尼不大情況下ωr≈ωdr;Ar為響應幅值;φr為初始相位。

由式(4)可見,變壓器繞組的自由振動響應信號為其各個模態響應信號的疊加。而理想的EMD分解可將信號分解成含有單一頻率的IMF信號,則每個IMF包含有一階模態響應信號,其中第i階模態響應信號xi(t)可以表示為

xi(t)=Aie-ζiωitcos(ωdit+φi)

(5)

則可得xi(t)的解析信號為

zi(t)=xi(t)+jH[xi(t)]=Ai(t)e-jθi(t)

(6)

式中,H[xi(t)]為xi(t)的Hilbert變換。

由式(5)和式(6)可見,Ai(t)與θi(t)表達式分別式(7)和式(8)所示:

Ai(t)=Aie-ζiωit

(7)

θi(t)=ωdit+φi

(8)

對式(7)取對數運算有

lnAi(t)=-ζiωit+lnAi

(9)

由式(8)和式(9)可見,變壓器繞組固有頻率與阻尼比可以通過求取相位與幅值隨時間變化斜率得到。

由式(5)可知,應用Hilbert變換對固有頻率與阻尼比進行識別的前提條件為EMD算法能夠將各模態響應分解到不同的IMF中,每個IMF僅含有單一的頻率分量。然而,標準的EMD分解存在模式混疊現象,一個IMF分量中往往含有多個頻率成分,這將嚴重影響模態參數識別精度。因此本文提出使用優化限制帶寬EMD算法對變壓器繞組振動信號進行分析。

1.3 限制帶寬EMD算法及其改進

Deering等[5]提出在EMD分解中引入屏蔽信號來抑制其模式混疊現象,提高分接精度。何啟源[6]將該方法進行了進一步改進,并稱之為限制帶寬模態分解算法(Restricted Band EMD, RBEMD),其實現步驟為:

(1) 對于原始信號x(t),使用標準的EMD算法將信號分離出第1個IMF,定義為IMF1。由于模式混疊,其中可能含有多個頻率分量。

(2) 對IMF1進行Hilbert變換,并加入屏蔽信號頻率,如式(10)所示:

(10)

式中:a1(i),f1(i)分別為IMF1瞬時幅值與瞬時頻率;M為帶寬系數。

(3) 構造限制帶寬信號為:

s(t)=A0sin(2πft)

(11)

(4) 對信號y(t)=x(t)+s(t)進行EMD分解,得到其第一個IMF,定義為IMFs1。則原始信號的第一階本征模態函數為h1(t)=IMFs1-s(t)。

(5) 從原始信號中減去h1(t),將余量作為原始信號,重復上述過程直至分解出所有的IMF。

由前述分解步驟可見,步驟(2)中屏蔽信號的引入是該算法的關鍵。EMD具有二進濾波特性,在原始信號中加入屏蔽信號能夠調整EMD各階IMF帶寬,從而達到分離相鄰頻率分量的目的。其中屏蔽頻率的選取將對算法的模式混疊抑制能力產生較大影響,而屏蔽頻率主要由步驟(4)中帶寬系數M決定。目前,絕大多數文獻通過根據經驗確定M的取值,即給定M=1.42[10]。然而,該帶寬系數并非對所有信號具有普適性,不合適的M值將會顯著降低算法對模式混疊的抑制能力。針對這一問題,本研究在此提出使用基于粒子群的優化算法,對帶寬系數M值進行最優選取,并將其用在變壓器繞組的模態參數識別之中。

粒子群算法是由Kennedy等[11]提出的一種群體智能演化計算技術,其優化過程為:

隨機初始化粒子群N,第i個粒子用d維向量xi和vi分別表示其位置和速度,將其代入優化目標函數得出適應值,更新粒子速度和位置,通過迭代尋求最優解。粒子速度和位置的更新方程為:

(12)

(13)

該方法在迭代過程中每個粒子根據自身和群體發現的最優值修正自己的前進方向和速度,最終得到全局最優解,具有收斂速度快、局部搜索能力強等優點。

為了使選取的屏蔽信號頻率具有最佳的抑制模式混疊能力,本文在此定義模式混疊指數KM作為IMF函數單頻性是否良好的評價指標,它也是PSO優化算法的目標函數。其定義為

(14)

式中:A1為IMF中主頻率分量幅值,一般指研究中所關心的頻率或IMF中頻率較高且能量較大的頻率分量;A2至AN為IMF中除主頻率以外的頻率分量幅值。

若分解算法對模式混疊有良好抑制能力,所得IMF單頻性較好,則KM接近于1。

2 仿真計算

為了說明所提出的優化RBEMD法對變壓器繞組模態模態參數識別的有效性,本文在此通過構造雙自由度系統自由振動信號進行仿真計算。同時,為了體現優化RBEMD方法抑制模式混疊問題的優勢,仿真信號中給出的系統兩階固有頻率較為接近,如式(15)所示。圖1為仿真信號的時域波形。

x(t)=5e-64.18tsin(365×2πt)+

6e-88.61tsin(415×2πt)

(15)

圖1 自由振動信號時域波形

圖2為使用RBEMD算法對圖1所示的自由振動信號進行分解的結果,其中,帶寬系數M由文獻中提供的經驗值確定,即有M=1.42[10]。由圖可見,使用經驗的屏蔽信號頻率并不能將信號完全分解,IMF1中依然出現了較大能量的低頻信號,而IMF2中也有較多高頻分量出現。由前述模態參數識別過程可知,這將對模態參數識別結果產生影響。

圖3 改進RBEMD算法的分解結果

應用優化RBEMD算法對仿真信號進行分解的結果分別如圖3所示,其中,選取的最佳帶寬系數M為1.62。由圖可見,分解后的各個IMF中幾乎僅含有一個頻率分量,具有較好的單頻性,有效地抑制了模式混疊現象。

表1 仿真信號的模態參數識別結果

表1為根據前述分解結果對IMF進行Hilbert變換得到的仿真信號的固有頻率及阻尼比。改進前后的限制帶寬EMD方法識別結果與理論值對比如表1所示。由表可見,使用已有的限制帶寬EMD算法得到仿真信號的固有頻率與阻尼比與理論值相比存在較大誤差。由于模式混疊未能有效抑制,IMF1中存在低頻分量導致2階固有頻率低于理論值;IMF2中存在的高頻分量導致1階固有頻率高于理論值。而基于PSO優化的改進限制帶寬EMD算法識別得到的各階模態參數與理論值匹配良好,精度較高。

3 變壓器繞組模態實驗描述

實驗對象為一臺10 kV油浸式配電變壓器,使用三維振動加速度測試振動信號,圖4為變壓器繞組及實驗測點布置實物圖。試驗時,使用電磁激振器對繞組進行軸向激振所用激振信號為白噪聲信號,帶寬為20 kHz。白噪聲信號發生器經功率放大直接驅動電磁激振器,對變壓器繞組進行激振。使用DH5922振動信號采集與分析系統拾取各測點加速度傳感器的振動信號,采樣頻率為2.56 kHz,得到繞組的綜合頻響函數。

圖4 變壓器繞組及測點布置實物圖

圖5 變壓器繞組的頻響函數

圖6 變壓器繞組的自由振動時域信號

圖5為測試得到的變壓器繞組頻響函數, 使用半對數坐標表示。由圖可見,繞組頻響函數有4個主要頻率分量,表明繞組在1 000 Hz內存在4階固有頻率。對頻響函數進行傅里葉逆變換,即可得到繞組自由振動時域信號,如圖6所示。

4 結果分析

4.1 優化RBEMD的繞組模態參數識別結果

使用本文提出的優化RBEMD算法對圖5所示的變壓器繞組自由振動信號分解時得到的前4階IMF及其頻譜如圖7所示。其中使用粒子群算法中得到的最優帶寬系數為M=1.55,模式混疊指數KM=0.84。對圖7所示的各個IMF進行Hilbert變換,求取其瞬時幅值與相位。進一步對瞬時幅值與相位進行最小二乘擬合,則可根據式(8)和式(9)得到變壓器繞組的固有頻率及阻尼比。

圖7 優化RBEMD的分解結果

圖8和圖9分別為變壓器繞組前4階IMF幅值和相位的計算結果及其擬合曲線,其中,圖8中的縱坐標為IMF幅值的自然對數值。根據各個擬合曲線的斜率,可得到變壓器繞組的前4階固有頻率及阻尼比如表2所示。

表2 變壓器繞組模態參數識別結果

圖8 IMF的幅值及其擬合曲線

圖9 IMF的相位及其擬合曲線

4.2 PolyMax法的繞組模態參數識別結果

為了進一步說明前述計算結果的正確性,本文在此使用多參考點最小二乘復頻域法PolyMax法根據模態實驗中測得的振動頻響函數對變壓器繞組模態參數進行識別。該方法由Peeters等[12]提出,它以頻響函數為基礎,采用右矩陣分數模型,并對非線性參數做線性化處理獲得誤差方程,再用最小二乘法求出分子、分母多項式系數。再將多項式系數代入誤差方程,可獲得縮減正則方程,求解該方程即可求得模態參數[12]。計算結果如表3所示。

表3 PolyMax法模態參數識別結果

由表2和表3可見,文中提出的方法與PolyMax法在識別繞組的各階固有頻率時吻合程度良好。在阻尼比識別方面,這兩種方法對第1階、第3階和第4階的阻尼比識別結果較為接近,誤差均在0.5 %以內,而第2階阻尼比的識別結果有較大誤差。為了進一步對第2階阻尼比結果進行說明,文中在此使用工程中廣泛應用的半功率帶寬法對2階固有頻率對應的阻尼比進行估計[13],其計算公式為:

(16)

由變壓器繞組頻響函數及式(16)計算可得,二階固有頻率對應阻尼比約為2.94 %,與優化RBEMD算法識別結果接近,而PolyMax算法識別結果存在較大誤差。究其原因在于,PolyMax法需要選擇合適的擬合頻段進行計算。但對于變壓器繞組頻響函數,其一、二階固有頻率相距緊密,且均有較大阻尼,頻響函數呈現“寬胖”形態。頻響函數在前兩階固有頻率對應峰值附近出現疊加,兩者相互影響。因此,PolyMax法在二階固有頻率所在頻段進行最小二乘擬合計算阻尼比時會出現較大誤差。

4.3 RBEMD的繞組模態參數識別結果

為了說明EMD固有的模態混疊對變壓器繞組模態參數識別結果的影響,本文在此給出了使用RBEMD法的模態參數識別結果。

圖10 IMF的頻譜分布

圖9為使用RBEMD法對變壓器繞組自由振動信號分解時的得到的各個IMF的頻譜分布。與圖6(b)的分解結果比較可見,各個IMF均有不同程度的模態混疊現象出現。另外,雖然720 Hz與490 Hz模態響應分量被分到兩個IMF中,但除主頻率之外仍存在較大能量的其它頻率分量;而230 Hz與320 Hz模態響應則出現了較嚴重的模態混疊現象。IMF3和IMF4中混疊了多個較大能量的頻率分量,這必將嚴重影響繞組模態參數的識別結果。

表4 RBEMD法模態參數識別結果

類似地,根據圖9的計算結果及前述計算方法,可得到對應的繞組模態參數識別結果如表4所示。對比表2、表3和表4結果可見,受模態混疊的影響,RBEMD法對繞組的模態參數識別結果也存在一定的誤差,特別是模式混疊嚴重的第1階和第2階固有頻率處的識別誤差高于20 %。這是由于傳統RBEMD未對屏蔽信號頻率進行優化選擇,由經驗公式(10)計算得到屏蔽頻率并不一定適合各IMF的分解。相應地,IMF3與IMF4中混入的高頻分量在對相位進行最小二乘擬合求固有頻率過程中產生了較大誤差,進而進一步影響了通過幅值擬合計算得到的阻尼比識別結果。由此可見,相比傳統的RBEMD算法,本文所提出的優化RBEMD算法可以顯著提高變壓器繞組模態參數識別的準確率。

此外,對實驗變壓器繞組的模態參數識別結果表明,該變壓器繞組各階固有頻率均遠離其正常運行時的激勵頻率100 Hz,不易出現共振。且繞組阻尼較大,在突然遭遇較強電動力后可以很快恢復穩定,繞組在正常工況下整體穩定性較高。

5 結 論

本文在對一臺10 kV變壓器繞組進行模態實驗的基礎上,提出一種優化RBEMD算法對繞組的模態參數進行識別。研究結果表明:

(1) 所提出的優化RBEMD法較好地識別出了變壓器繞組的前四階固有頻率和阻尼比,計算結果與PolyMax法的識別結果吻合良好,且計算精度更高,說明了該方法的有效性;

(2) 所提出的優化RBEMD在準確識別變壓器繞組模態參數的同時,通過選取最優屏蔽頻率,較好地抑制了EMD分解過程中的模態混疊現象;

(3) 試驗用變壓器繞組的各階固有頻率均遠離其正常運行時的激勵頻率100 Hz,且阻尼比較大,說明該變壓器繞組的結構設計較為合理。

文中模態試驗測試及模態參數識別方法雖以一臺10 kV變壓器為研究對象,但繞組模態參數測試方法和所提出的模態參數識別方法均可用于更高電壓等級的電力變壓器繞組,因此,本文的研究結果可為變壓器繞組的優化設計、制造及振動監測提供理論依據,具有較大的工程應用價值。

[1] Beltle M, Tenbohlen S. Usability of vibration measurement for power transformer diagnosis and monitoring [C]. Proceedings of 2012 International Conference on Condition Monitoring and Diagnosis, Bali, Indonesia, 2012: 281-284.

[2] 郭健,林鶴云,徐子宏,等. 用有限元方法分析電力變壓器繞組軸向穩定性[J],高電壓技術,2007, 33(11): 209-212.

GUO Jian, LIN He-yun, XU Zi-hong, et al. Analysis of axial stability of power transformer windings using finite element[J]. High Voltage Engineering, 2007, 33(11): 209-212.

[3] 林循泓,潘得引,臧朝平,等. 振動模態參數識別及其應用[M].南京:東南大學出版社,1994.

[4] 周云,易偉建. 用PolyMAX方法進行彈性地基板的實驗模態分析[J].振動與沖擊,2007,26(7):139-145.

ZHOU Yun, YI Wei-jian, Experimental modal analysis of a slab on elastic foundation by PolyMAX method[J]. Journal of Vibration and Shock, 2007, 26(7): 139-145.

[5] Deering R, Kaiser J F. The use of a masking signal to improve empirical mode decomposition [C]. IEEE International Conference on Acoustic, Speech and Signal Processing, USA, 2005 : 485-488.

[6] 何啟源. 基于現代時頻分析的環境激勵模態參數識別方法研究[D].重慶: 重慶大學, 2009.

[7] Yang J N, Lei Y, Pan S, et al. System identification of linear structures based on hilbert-huang spectral analysis Part1: normal modes [J]. Earthquake Engineering and Structural Dynamics, 2003, 32: 1443-1467.

[8] 秦世強,蒲黔輝,施洲. 環境激勵下大型橋梁模態參數識別的一種方法[J].振動與沖擊,2012, 31(2): 95-100.

QIN Shi-qiang, PU Qian-hui, SHI Zhou. A method of modal parameter identification using ambient vibration measurements for a large-scale bridge[J]. Journal of Vibration and Shock, 2012, 31(2): 95-100.

[9] 司馬文霞,王荊,楊慶,等. Hilbert-Huang變換在電力系統過電壓識別中的應用[J].高電壓技術,2010,36(6):1480-1486.

SIMA Wen-xia, WANG Jin, YANG Qing, et al. Application of Hilbert-Huang transform to power system over-voltage recognition[J]. High Voltage Engineering, 2010, 36(6): 1480-1486.

[10] 馬燕峰,趙書強. 用改進的Hilbert-Huang變換辨識電力系統低頻振蕩[J]. 高電壓技術, 2012, 38(6): 1492-1499.

MA Yan-feng, ZHAO Shu-qiang. Identification of low-frequency oscillations in power system based on improved Hilbert-Huang transformer[J]. High Voltage Engineering, 2012, 38(6): 1492-1499.

[11] 朱顯輝,崔淑梅,師楠,等. 電動汽車電機故障時間的粒子群優化灰色預測[J],高電壓技術,2012,38(6):1391-1396.

ZHU Xian-hui, CUI Shu-mei, SHI Nan, et al. Grey prediction model of electric vehicle motor based on particle swarm optimization[J]. High Voltage Engineering, 2012, 38(6): 1391-1396.

[12] Peeters B, Van der Auweraer H, et al. The polyMAX frequency domain method: A new standard for modal parameter estimation [J]. Shock and Vibration, 2004, 11: 395-409.

[13] 應懷樵,劉進明,沈松. 半功率帶寬法與INV阻尼計法求阻尼比的研究[J]. 噪聲與振動控制, 2006,26(2): 4-6.

YING Huai-qiao, LIU Jin-ming, SHEN Song. Half-power bandwidth method and INV damping ratio solver study [J]. Noise and Vibration Control, 2006,26(2): 4-6.

猜你喜歡
模態變壓器振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
理想變壓器的“三個不變”與“三個變”
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
開關電源中高頻變壓器的設計
中立型Emden-Fowler微分方程的振動性
一種不停電更換變壓器的帶電作業法
變壓器免維護吸濕器的開發與應用
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 日韩 欧美 国产 精品 综合| 亚洲国产成熟视频在线多多| 国产精品天干天干在线观看| 亚洲美女AV免费一区| 激情综合网址| 国产精品主播| 国产成人精品亚洲日本对白优播| 无码人妻免费| 美女国内精品自产拍在线播放| 人妻21p大胆| 免费看美女毛片| 亚洲国产日韩在线观看| 国产精品19p| 亚洲天堂视频网| 狠狠色香婷婷久久亚洲精品| 狠狠操夜夜爽| 国产swag在线观看| 天天色天天综合| 色综合五月婷婷| 国产欧美视频在线| 首页亚洲国产丝袜长腿综合| 亚洲中文字幕久久精品无码一区| 无套av在线| 亚洲码一区二区三区| 综合色亚洲| 亚洲精品视频免费观看| 久久这里只有精品2| 91区国产福利在线观看午夜 | 九九这里只有精品视频| 国产精品黑色丝袜的老师| 一级毛片免费不卡在线视频| 精品国产自在在线在线观看| 久热re国产手机在线观看| 亚洲精品成人片在线观看| 国产男女XX00免费观看| 国产成人毛片| 精品夜恋影院亚洲欧洲| 精品国产中文一级毛片在线看| 国产成人精品无码一区二| 欧美一区二区三区国产精品| 狠狠色丁香婷婷| 免费女人18毛片a级毛片视频| 亚洲成人在线播放 | 婷婷色中文网| 欧美成人日韩| 干中文字幕| 99在线观看精品视频| 高潮爽到爆的喷水女主播视频| 夜色爽爽影院18禁妓女影院| 欧美不卡二区| 成人午夜精品一级毛片| 欧美成人第一页| 亚洲va视频| 99久视频| 日本不卡在线视频| 毛片大全免费观看| 国产主播福利在线观看| 国产丝袜无码精品| 国产精品香蕉在线观看不卡| 波多野结衣第一页| 特黄日韩免费一区二区三区| 国产在线无码一区二区三区| 九九热免费在线视频| 免费高清毛片| h网址在线观看| 成年人久久黄色网站| 日韩一区二区三免费高清| 中文字幕亚洲综久久2021| 亚洲日本在线免费观看| 亚洲欧美日韩中文字幕一区二区三区| 欧美福利在线| 色婷婷色丁香| 99福利视频导航| 欧美综合在线观看| 国产精品第一区| 日本不卡在线播放| 国产三级韩国三级理| 日韩不卡高清视频| 亚洲成人77777| 国产日韩久久久久无码精品| 狠狠色狠狠色综合久久第一次| a级高清毛片|