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

應(yīng)用改進復Morlet小波識別電力變壓器繞組模態(tài)參數(shù)

2014-09-05 08:39:58王豐華廖天明金之儉
振動與沖擊 2014年6期
關(guān)鍵詞:模態(tài)變壓器振動

駱 波, 王豐華, 廖天明, 金之儉

(1. 上海交通大學 電氣工程系,上海 200240; 2. 上海電力公司 運檢部,上海 200240)

作為電力系統(tǒng)中重要組成部分,電力變壓器的穩(wěn)定性、可靠性與整個電網(wǎng)安全運行密切相關(guān)。統(tǒng)計分析表明[1],多數(shù)變壓器故障來自繞組變形,而變壓器繞組主要由線餅、絕緣墊塊、壓緊裝置等構(gòu)成,結(jié)構(gòu)非常復雜,在變壓器運行過程中發(fā)生的各種繞組變形故障均會使其機械特性改變、模態(tài)參數(shù)變化,導致變壓器繞組振動特性及振動信號改變。因此,對繞組結(jié)構(gòu)進行模態(tài)分析對其振動特性、故障診斷、動力特性優(yōu)化設(shè)計及變形檢測等具有重要意義及工程應(yīng)用價值。

變壓器繞組模態(tài)參數(shù)主要包括固有頻率、振型、阻尼比等,能全面描述整個繞組機械系統(tǒng)固有動態(tài)特性,物理意義較明確。目前,獲取變壓器繞組類機械結(jié)構(gòu)模態(tài)參數(shù)途徑主要依賴繞組結(jié)構(gòu)的試驗測試,進而分析與參數(shù)識別。傳統(tǒng)模態(tài)參數(shù)識別方法主要包括最小二乘負指數(shù)法、時間序列分析法等時域辨識方法[2]及最小二乘圓擬合法、分區(qū)模態(tài)綜合法、頻域總體識別法[3]等兩類。但此方法大多限于時域或頻域中單獨對測試信號進行識別,參數(shù)辨識精度有限;而且對測試信號中噪聲較敏感,固有頻率識別與定階存在較大主觀性。因此,對電力變壓器繞組類結(jié)構(gòu)復雜強非線性機械結(jié)構(gòu)系統(tǒng)而言,需尋求新的模態(tài)參數(shù)識別方法以提高模態(tài)參數(shù)識別準確性,此亦為變壓器繞組結(jié)構(gòu)健康監(jiān)測的重要前提。

小波變換法因時頻局部化特性良好、可使多自由度系統(tǒng)模態(tài)自動解耦等特點被廣泛用于結(jié)構(gòu)模態(tài)參數(shù)識別領(lǐng)域[4]。但在電力變壓器繞組類機械結(jié)構(gòu)系統(tǒng)模態(tài)參數(shù)識別領(lǐng)域鮮有應(yīng)用。故本文以某10 kV電力變壓器繞組為研究對象,嘗試用復Morlet小波變換法對繞組軸向激振實驗振動響應(yīng)信號進行時頻分析,用改進Crazy Climber算法提取小波脊線,以期獲得較準確、高效的模態(tài)參數(shù)識別結(jié)果。為說明計算結(jié)果的正確性,本文給出使用多參考最小二乘復頻域法PolyMAX法獲所得繞組模態(tài)參數(shù)識別結(jié)果。

1 理論基礎(chǔ)

1.1 小波變換與小波脊

對任意信號x(t)∈L2(R),其連續(xù)小波變換為:

(1)

式中:a,b為尺度因子、平移變量;ψ*表示共軛。

本文采用復Morlet小波作為母小波,原因為[5]:① 該小波的實部、虛部幅值均為按指數(shù)衰減的簡諧振動信號,與動態(tài)系統(tǒng)自由響應(yīng)信號一致;② 該小波具有單值頻率,若分析信號與某尺度小波高度相關(guān),即在該尺度下小波系數(shù)較大,則小波頻率可表征產(chǎn)生該響應(yīng)信號的動態(tài)系統(tǒng)固有頻率。因此對密集模態(tài)的識別效果較好。Morlet小波時頻域定義為:

(2)

(3)

式中:ω0為Morlet小波中心角頻率。

由式(2)、(3)可見,Morlet小波的時、頻域波形均為高斯窗函數(shù)形式,已具備分離出各階模態(tài)能力。對形如x(t)=A(t)cos(ωt)的漸進單頻信號進行Morlet小波變換有:

(4)

式中:Zx(t)為信x(t)的解析形式;ψa,b(a,b)為母小波;φ(a,b)(t)為相位。分別表示為:

Zx(t)=x(t)+jH[x(t)]=As(t)ejφ(t)

(5)

ψa,b(a,b)=Aψ(t)ejθ(t)

(6)

(7)

定義使φ′(x)=0且φ″(x)≠0的點x=x0稱為φ(x)=0的一階平穩(wěn)相位點[6]。在小波變換域平面內(nèi),若平穩(wěn)相位點t(a,b)滿足條件t(a,b)=b,則定義滿足條件的點(a,b)為小波變換脊點,脊點集合稱為小波變換脊線。

若t=t0為一階平穩(wěn)點,則對式(7)求導有:

(8)

令a=ar(b)為小波脊線,將t0=b代入上式有:

(9)

令φ′(b)=ωs(b),據(jù)式(9)得信號在時間b處頻率ωs(b)為[5]:

(10)

由上式可見,信號x(t)在時間b處頻率與小波在尺度a=ar(b)處值有對應(yīng)的定量關(guān)系。因此,只要獲得小波變換脊線所在尺度位置即可求出原信號在該處頻率特性。與小波脊對應(yīng)的小波系數(shù)連線WTs(ar(b),b)稱為小波骨架,其實部對應(yīng)信號本身[7]。

1.2 模態(tài)參數(shù)識別方法

單自由度線性系統(tǒng)自由振動運動微分方程解的形式為[6]:

x(t)=Ae-εωntsin(ωdt+η)

(11)

對式(11)進行小波變換:

(12)

將上式在t=b處進行Taylor級數(shù)展開:

(13)

忽略高階分量,上式可近似為:

(14)

選擇Morlet小波為母小波,將式(2)、(12)代入式(15)得自由衰減振動信號x(t)的Morlet小波變換:

(15)

得小波系數(shù)模值為:

(16)

當尺度a=ω0/ωd時,小波系數(shù)模取得極大值,幅值、相位分別為:

(17)

∠WT(a,b)=ωdb+η

(18)

a=ω0/ωd確定的尺度值為該小波變換的一階平穩(wěn)點,因此可由小波脊線提取方法獲得與該尺度位置對應(yīng)的固有頻率。此時尺度a對應(yīng)頻率即為系統(tǒng)固有頻率。而該尺度下小波系數(shù)模值對數(shù)對時間的斜率即為無阻尼固有頻率ωn與阻尼比ε的乘積,設(shè)斜率為k1,固有頻率為fk,則得系統(tǒng)在該階固有頻率下阻尼比:

ε=-k1/(2πfk)

(19)

由此可見,對系統(tǒng)自由振動信號進行小波變換后,提取小波脊線的尺度位置即可識別出系統(tǒng)該階固有頻率及阻尼比。對多自由度系統(tǒng),由于小波變換過程實為信號分解過程,只在子小波定義的窗內(nèi)對信號進行分解。在窗外,因小波具有快速衰減特性,幾乎不分解信號,此時小波變換過程即為時頻濾波器對信號帶通濾波。對特定頻率尺度a,只有與該參數(shù)對應(yīng)的模態(tài)頻率成分信號才能通過小波濾波器,其它頻率信號對此階信號影響可忽略不計。故多自由度系統(tǒng)可變?yōu)閱巫杂啥认到y(tǒng)分析[9]。

小波變換識別模態(tài)參數(shù)算法的輸入應(yīng)為系統(tǒng)按指數(shù)衰減的正弦自由振動信號。而變壓器繞組軸向激振實驗中采集的振動信號為白噪聲激勵下受迫振動,可通過計算系統(tǒng)頻響函數(shù),獲得系統(tǒng)脈沖響應(yīng)函數(shù):

(20)

式中:X(jω)為信號x(t)傅氏變換;Y(jω)為實驗中所采信號傅氏變換。由于系統(tǒng)脈沖響應(yīng)函數(shù)h(t)即為頻響函數(shù)H(jω)的傅里葉逆變換,而脈沖響應(yīng)函數(shù)為系統(tǒng)自由振動信號,滿足小波算法的輸入要求。

1.3 基于改進Crazy Climber 算法的小波脊線提取

小波變換模態(tài)參數(shù)識別方法關(guān)鍵為小波脊線的提取。目前該提取方法主要有按定義計算的相位法及模值法、改進的局部模極大值法等[10]。但均存在對信號噪聲較敏感、不能有效提取多條脊線等問題。因此,為能在噪聲環(huán)境下準確進行多階模態(tài)參數(shù)識別,本文引入自適應(yīng)優(yōu)化思想到Crazy Climber算法[11]中,完成復Morlet小波變換系數(shù)矩陣中小波脊線的提取。設(shè)信號經(jīng)小波變換后時頻矩陣為B×A矩陣,矩陣中每點(i,j)的小波變換系數(shù)值為M(i,j)。則據(jù)Crazy Climber算法提取小波脊線的基本過程為:

(1) 初始化N個可移動爬升點、小波變換視頻矩陣尺寸相同爬升點觀測矩陣B及度量矩陣D,并將N個可移動爬升點均勻分布于時頻空間,度量矩陣D及觀測矩陣B初始化為零矩陣;

(2) 定義系統(tǒng)循環(huán)參量T0=max(M)-min(M),并令系統(tǒng)當前循環(huán)參量為Tt=T0,記錄可移動爬升點初始位置為Bk(t)=(i,j),k∈[1,2,…,N];

(3) 在t時刻,設(shè)可移動爬升點對應(yīng)位置為Bk(t)=(i,j)。t+1時刻,在不考慮邊界點時,可移動爬升點對應(yīng)位置Bk(t+1)=(i′,j′)確定規(guī)則為:①t時刻i值,按50%概率左移或右移一格,即i′=i-1或i′=i+1;②t時刻j值,按50%概率上移或下移一格,即j′=j-1或j′=j+1;若小波變換系數(shù)矩陣M(i′,j′)>M(i,j),則該點垂直移動,有Bk(t+1)=(i′,j′),否則該點按概率p沿垂直方向移動,即Bk(t+1)=(i′,j′),按概率1-po沿垂直方向不移動,即j′=j,Bk(t+1)=(i′,j′),概率p的計算公式為:

(21)

(4) 移動結(jié)束后,在度量矩陣D的相應(yīng)位置(i′,j′)上增加度量值M(i′,j′);

(5) 重復步驟(3)、(4),直至系統(tǒng)當前循環(huán)變量低于設(shè)定閾值Tf;

(6) 設(shè)定合適密度閾值Td,對最終所得度量矩陣D進行篩選低于該閾值的矩陣元素置零,此為降低信號中噪聲對計算結(jié)果影響、提高模態(tài)參數(shù)識別精度之關(guān)鍵;

(7) 據(jù)所得度量矩陣D在時間方向(j向)進行遞歸,給定任意點(i,j),在(i-Δi,j)、(i-Δi,j±Δj)中尋找最優(yōu)相鄰點,并形成一條脊線;

(8) 重復步驟(7),直至所有滿足要求的點均在脊線中,形成整個時頻平面的脊線。

(1)~(6)步為隨機移動階段,(7)~(8)步為鏈接遞歸階段。在脊線鏈接遞歸階段,遍歷網(wǎng)格大小選取直接影響小波脊線提取后形態(tài)。網(wǎng)格越小,提取脊線越平滑,模態(tài)參數(shù)辨識結(jié)果亦更準確。而網(wǎng)格大小亦直接影響密度閾值Td的選取,為保證脊線提取效果及模態(tài)參數(shù)識別精度,本文引入自適應(yīng)優(yōu)化算法確定密度閾值及網(wǎng)格尺寸的選取。過程為:① 在鏈接遞歸初始階段,引入脊線參量s為一設(shè)定常數(shù),其含義為所需提取的小波脊線數(shù)量。② 在鏈接遞歸過程中,將密度閾值Td設(shè)為相對較小數(shù)值,以確保真實模態(tài)能被識別。當一次脊線提取完后對模態(tài)數(shù)量進行統(tǒng)計,若超過脊線參量s,則增加密度閾值Td的數(shù)值,繼續(xù)循環(huán)該過程直至提取的脊線數(shù)量與模態(tài)參量一致;否則減小密度閾值Td的數(shù)值,繼續(xù)循環(huán)該過程直至提取的脊線數(shù)量與模態(tài)參量一致。在算法隨機移動過程較充分情況下,該自適應(yīng)過程可取得較理想結(jié)果。

2 仿真分析

為說明所提結(jié)合改進Crazy Climber算法復Morlet小波變換進行信號處理及模態(tài)識別的優(yōu)越性,構(gòu)造自由衰減仿真信號為:

x(t)=eε1ωn1πtsin(ωd1t)+

eε2ωn2πtsin(ωd2t+θ)

(22)

式中:ωd1≈ωn1=160π;ωd2≈ωn2=220π;ε1=-0.13;ε2=-0.08;t為時間;采樣頻率256 Hz,采樣時間1 s。為說明信號中噪聲對用Crazy Climber算法提取小波脊線影響,在式(22)中加入信噪比SNR=-3 dB的白噪聲。

圖1為信號x(t)時域波形。圖2為采用復Morlet小波與原Crazy Climber算法所得小波脊線。由兩圖看出,信號x(t)的兩個頻率分量分別在80 Hz及110 Hz附近,即式(22)仿真信號中的頻率分量。說明用復Morlet小波與Crazy Climber算法計算結(jié)果的正確性。但在原始信號中加入信號比為SNR=-3 dB的白噪聲后,提取的小波脊線尤其110 Hz頻率分量對應(yīng)的小波脊線出現(xiàn)彎曲現(xiàn)象,說明信號中的噪聲對原Crazy Climber算法的小波脊線提取結(jié)果有影響。

圖1 信號x(t)的時域波形

用改進Crazy Climber算法對加入信號比SNR= -3 dB白噪聲的仿真信號進行小波脊線提取后計算結(jié)果見圖3。由3圖看出,小波脊線彎曲現(xiàn)象已消失,表明本文改進Crazy Climber算法的有效性。

圖2 用Crazy Climber算法所得小波脊線圖

3 變壓器繞組軸向激振實驗分析

3.1 變壓器繞組軸向激振實驗描述

實驗對象為一臺餅式繞組結(jié)構(gòu)10 kV變壓器。用三維振動加速度測試振動信號,測點布置實物見圖4。試驗時用電磁激振器對繞組進行軸向激振,激振位置見圖5。所用激振信號為白噪聲信號,帶寬20 kHz。白噪聲信號發(fā)生器經(jīng)功率放大直接驅(qū)動電磁激振器對變壓器繞組進行激振。用DH5922振動信號采集與分析系統(tǒng)拾取各測點加速度傳感器振動信號,獲得繞組綜合頻響函數(shù)。

圖4 繞組測點布置圖

3.2 結(jié)果與分析

圖6為根據(jù)變壓器繞組軸向激振實驗所得變壓器繞組振動頻響函數(shù)(Frequency Response Function, FRF),用半對數(shù)坐標表示。由圖6看出,頻響函數(shù)在229 Hz,322 Hz,425 Hz,724 Hz左右有明顯疊加峰值結(jié)果。對其進行傅里葉逆變換,獲得自由振動信號見圖7。

圖6 振動頻響函數(shù)

圖7 自由振動時域信號

選復Morlet小波為小波母函數(shù),對圖7的自由振動信號進行連續(xù)小波變換,獲得小波變換時頻見圖8。由圖8看出,小波變換已將各階固有頻率所在頻帶清晰分離。但僅由時頻圖較難獲得精確固有頻率,需進一步提取小波脊。圖9為用Crazy Climber算法據(jù)圖8小波變換結(jié)果中提取的小波脊線圖。由圖9看出,繞組前四階固有頻率分別為227 Hz,330 Hz,427 Hz,732 Hz。分別提取四頻率處小波系數(shù),據(jù)式(19)得阻尼比分別為4.71%,4.47%,4.97%,3.05%。

為說明用復小波變換識別模態(tài)參數(shù)結(jié)果的正確性,本文亦用PolyMAX方法對測試所得振動數(shù)據(jù)模態(tài)參數(shù)進行識別[12]。表1為兩種方法計算結(jié)果。由表1看出,兩種方法對繞組固有頻率識別一致性良好,說明用復小波變換法識別的繞組模態(tài)參數(shù)的正確性。此外,據(jù)PolyMAX方法計算所得各階固有頻率阻尼比差別較大,尤其二階固有頻率對應(yīng)阻尼明顯小于其它各階,此由于PolyMax法對頻響函數(shù)擬合時噪聲影響所致。可見PolyMax在識別過程中受振動噪聲影響較大。當頻響函數(shù)存在較大噪聲時,PolyMax計算結(jié)果對頻率擬合范圍較敏感,需根據(jù)頻響函數(shù)峰值預(yù)先確定擬合區(qū)間而存在較大主觀性。因此,基于小波變換的模態(tài)識別方法在變壓器繞組模態(tài)識別中抗干擾性更好、識別精度更高。

4 結(jié) 論

(1) 本文根據(jù)實體變壓器繞組軸向激振實驗測試結(jié)果,提出用復Morlet小波變換及改進Crazy Climber算法對變壓器繞組模態(tài)參數(shù)進行識別。計算結(jié)果與用PolyMAX方法識別結(jié)果的良好吻合說明該方法的正確性。變壓器繞組前四階固有頻率均遠離100 Hz電動力激勵頻率,說明變壓器繞組結(jié)構(gòu)設(shè)計較合理。

(2) 本文所提復Morlet小波變換及改進Crazy Climber算法在變壓器繞組類復雜結(jié)構(gòu)模態(tài)參數(shù)識別中優(yōu)勢明顯、識別精度高、抗干擾性能強,且能清晰刻畫信號能量隨時間頻率分布。研究結(jié)果可為變壓器繞組結(jié)構(gòu)振動特性分析及基于振動分析法的變壓器繞組狀態(tài)監(jiān)測提供依據(jù)。

參 考 文 獻

[1]金文龍,陳建華,李光范,等.全國110 kV及以上等級電力變壓器短路損壞事故統(tǒng)計分析[J]. 電網(wǎng)技術(shù),1999,23(6):21-25.

JIN Wen-long, CHEN Jian-hua, LI Guang-fan,et al. Statistics and analysis on power transformer damage caused by short circuit fault in 110kV and higher voltage classes[J]. Power System Technology, 1999,23(6):21-25.

[2]林循泓,潘得引,臧朝平,等,振動模態(tài)參數(shù)識別及其應(yīng)用[M].南京:東南大學出版社,1994.

[3]周 云,易偉建.用PolyMAX方法進行彈性地基板的實驗?zāi)B(tài)分析[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.

[4]何啟源,湯寶平,程發(fā)斌. 基于修正Morlet小波的自適應(yīng)模態(tài)參數(shù)識別[J]. 中國機械工程,2007,18(20): 2476-2480.

HE Qi-yuan,TANG Bao-ping,CHENG Fa-bin. Modified Morlet wavelet-based adaptive modal parameter identification[J]. China Mechanical Engineering, 2007,18(20): 2476-2480.

[5]孫 鵬,丁幼亮,張勁泉,等. 基于Morlet小波變換的結(jié)構(gòu)密集模態(tài)參數(shù)識別[J]. 東南大學學報,2012,42(2):339-345.

SUN Peng, DING You-liang, ZHANG Jin-quan,et al. Modal identification of closely spaced modes based on Morlet wavelet transform[J]. Journal of Southeast University,2012,42(2):339-345.

[6]劉 寧. 系統(tǒng)模態(tài)參數(shù)小波辨識方法研究[D].天津:天津大學,2006.

[7]伊廷華,李宏男,王國新.基于小波變換的結(jié)構(gòu)模態(tài)參數(shù)識別[J].振動工程學報, 2006,19(1):51-56.

YI Ting-hua,LI Hong-nan,WANG Guo-xin. Structural modal parameter identification based on wavelet transform[J].Journal of Vibration Engineering, 2006,19(1):51-56.

[8]Lardies J, Ta M N, Berthiller M. Modal parameter estimation based on the wavelet transform of output data[J]. Archive of Applied Mechanics, 2004, 73(9-10): 718-733.

[9]羅光坤,張令彌. 基于Morlet小波變換的模態(tài)參數(shù)識別研究[J].振動與沖擊,2007,26(7):135-139.

LUO Guang-kun, ZHANG Ling-mi.Study on identification of modal parameters based on Morlet wavelet transformation[J]. Journal of Vibration and Shock, 2007, 26(7):135-139.

[10]Carmona R, Hwang W L. Characterization of signals by the ridges of their wavelet transforms[J].IEEE Transactions on Signal Proeessing,1997,45(10):2586-2590.

[11]Carmona R, Hwang W L, Torresani B. Multi-ridge detection and time-frequency reconstruction[J].IEEE Transactions on Signal Processing,1999, 47:480-492.

[12]謝小平,韓 旭,吳長德,等,基于PolyMAX方法的某轎車車身實驗?zāi)B(tài)分析[J]. 汽車工程,2009,31(5):440-444.

XIE Xiao-ping, HAN Xu, WU Chang-de,et al.Experimental modal analysis for a car Body-in-white based on PolyMAX method [J].Automotive Engineering, 2009, 31(5): 440-444.

猜你喜歡
模態(tài)變壓器振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
理想變壓器的“三個不變”與“三個變”
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
開關(guān)電源中高頻變壓器的設(shè)計
中立型Emden-Fowler微分方程的振動性
一種不停電更換變壓器的帶電作業(yè)法
變壓器免維護吸濕器的開發(fā)與應(yīng)用
國內(nèi)多模態(tài)教學研究回顧與展望
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
UF6振動激發(fā)態(tài)分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 成人韩免费网站| 免费a在线观看播放| 国产成人精品18| 成人午夜福利视频| 亚洲精品福利视频| 日韩乱码免费一区二区三区| 毛片一区二区在线看| 国产美女视频黄a视频全免费网站| 国产精品福利尤物youwu| 国产亚洲精| 国产精品视频3p| 欧美不卡视频一区发布| 免费一看一级毛片| 久久免费视频6| 久久a毛片| 亚洲综合18p| 久久综合结合久久狠狠狠97色| 成人免费一区二区三区| 国产白浆一区二区三区视频在线| 免费精品一区二区h| 九色在线观看视频| 午夜精品久久久久久久无码软件 | 毛片免费高清免费| 97国内精品久久久久不卡| 国产成人91精品| 黄色网页在线观看| 国产免费黄| 国产拍在线| 97在线观看视频免费| 久久久久久久久亚洲精品| 日韩美毛片| 国产精品自拍合集| 国产成人精品优优av| 狠狠色综合网| 国产va在线观看免费| 人妻精品全国免费视频| 大乳丰满人妻中文字幕日本| 久久国产高潮流白浆免费观看| 香蕉综合在线视频91| 五月六月伊人狠狠丁香网| 在线一级毛片| 国产高清在线观看| 91色综合综合热五月激情| 国产剧情一区二区| 四虎综合网| 免费啪啪网址| 97一区二区在线播放| 国产精品网曝门免费视频| 国产精品蜜臀| 国产成人精品视频一区视频二区| 亚洲一区第一页| 亚洲v日韩v欧美在线观看| 久久香蕉国产线| 国产福利一区二区在线观看| 日韩免费中文字幕| 97人人做人人爽香蕉精品| 国产91无毒不卡在线观看| 五月婷婷综合在线视频| 亚洲伊人久久精品影院| 亚洲最新网址| 98超碰在线观看| 在线观看无码a∨| 在线a视频免费观看| 婷婷色一区二区三区| 国产一国产一有一级毛片视频| 99久久精品免费看国产免费软件| 91久久国产热精品免费| 日本高清免费不卡视频| 不卡无码网| 日韩午夜伦| 香港一级毛片免费看| 狠狠色综合网| 91视频免费观看网站| 成人在线视频一区| 亚洲一区二区三区香蕉| 女人毛片a级大学毛片免费| 黄色在线不卡| 亚洲国产日韩一区| 人妻一区二区三区无码精品一区| 国产主播喷水| 波多野结衣一区二区三视频| 十八禁美女裸体网站|