亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        局部特征尺度分解結合局部均值解調的齒輪故障診斷

        2017-01-09 02:38:50崔偉成許愛強孟凡磊
        中國機械工程 2016年24期
        關鍵詞:調頻均值齒輪

        崔偉成 許愛強 李 偉 孟凡磊

        海軍航空工程學院,煙臺,264001

        局部特征尺度分解結合局部均值解調的齒輪故障診斷

        崔偉成 許愛強 李 偉 孟凡磊

        海軍航空工程學院,煙臺,264001

        為了準確地進行故障診斷,根據(jù)齒輪故障振動信號的多分量調幅-調頻特征,提出了一種新的解調方法——局部均值解調法,將之與局部特征尺度分解相結合進行齒輪故障診斷。該診斷方法首先對齒輪振動信號運用局部特征尺度分解,得到若干個內稟尺度分量,然后應用局部均值解調法求取每個分量的調頻分量,最后根據(jù)瞬時頻率的頻譜進行故障診斷。采用仿真信號將局部均值解調法與Hilbert解調法、經(jīng)驗調幅調頻分解法進行了對比,結果表明,局部均值解調法的精確性更好。通過齒輪故障振動數(shù)據(jù)的分析,驗證了局部特征尺度分解結合局部均值解調的故障診斷方法的有效性。

        局部特征尺度分解;局部均值解調;瞬時頻率譜;齒輪故障診斷

        0 引言

        齒輪傳動是一種常見的運動和動力傳遞方式,廣泛地應用于機械設備中。齒輪故障會導致機械設備的整體性能下降,甚至引起嚴重設備事故從而造成重大經(jīng)濟損失[1-2],因此,研究齒輪故障診斷方法具有實際意義。而齒輪的振動信號蘊含了齒輪的運行狀態(tài)信息,所以,通常通過振動信號分析進行齒輪故障診斷。

        齒輪故障的振動信號一般具有非線性、非平穩(wěn)的特性,為了準確地將信號分解,進而得到分量的局部特征,許多學者開展了時頻分析技術的研究[3-5]。其中,經(jīng)驗模態(tài)分解(empirical mode decomposition,EMD)是一個研究的熱點。作為一種有效的自適應時頻分析方法,EMD與Hilbert變換(Hilbert transform, HT)相結合,被廣泛地應用于齒輪故障診斷領域[6-8]。但EMD在使用過程中存在一些問題,如過包絡、欠包絡、端點效應、模態(tài)混疊等[6-8]。因此,國內外學者一直在尋求EMD的改進、替代方法。

        程軍圣等[9]在本征時間尺度分解(intrinsic time-scale decomposition, ITD)[10]的基礎上創(chuàng)造性地提出了局部特征尺度分解(local characteristic-scale decomposition,LCD)方法,該方法在保證分解所得分量物理意義清晰的基礎上,具有計算時間短、端點效應不明顯等優(yōu)勢,開拓了自適應時頻分析方法的新思路[11-12]。

        齒輪故障振動信號可以看成多個調幅-調頻信號的疊加,通過LCD可將信號分解成單個調幅-調頻分量之和,對每個分量進行解調再進行故障診斷是一種有效的方法。常用的解調方法有Hilbert解調和經(jīng)驗調幅調頻分解(empirical AM-FM decomposition, EAD)。Hilbert法的不足是:包絡本身光滑性較差,且端點效應明顯,求取信號瞬時頻率時會出現(xiàn)負頻率現(xiàn)象[13]。EAD法可避免Hilbert變換引起的端點效應及負頻率,但EAD法只是一種經(jīng)驗方法,缺乏嚴密的理論推導[14]。因此,研究新的解調方法具有一定的工程意義。

        為了準確分析齒輪故障振動信號蘊含的故障特征,本文提出了一種新的解調方法——局部均值解調法,并結合LCD進行故障診斷。

        1 基本理論

        1.1 LCD方法

        LCD方法將一個信號x(t)分解成若干個內稟尺度分量(intrinsic scale components,ISC)之和[10-12]。每個ISC分量必須滿足2個條件:①任意2個相鄰的極值點符號互異;②考慮所有的極值點Xk及對應的時刻τk(k=1,2,…,M,其中M為極值點的個數(shù))。

        取兩個相距最近且符號互異的極值點(τk,Xk)和(τk+2,Xk+2),按照下式定義一個τk+1時刻的函數(shù)值:

        (1)

        則Ak+1與極值點Xk+1的比值關系不變,即

        αAk+1+(1-α)Xk+1=0

        (2)

        式中,α∈(0,1)為一常量,典型地,α=0.5。

        為了篩選出ISC分量,LCD方法構造了信號的均值曲線,將均值曲線不斷地從原信號中分離,其迭代過程與EMD方法類似,不同之處在于均值曲線的構造方法。LCD方法用均值點構造均值曲線,均值點的定義為

        Lk+1=αAk+1+(1-α)Xk+1

        (3)

        式中,Lk+1為均值點,在迭代結束時數(shù)值為0。

        標準LCD算法采用分段線性方法由均值點Lk計算均值曲線。為了得到更光滑的ISC分量,本文采用三次樣條方法替代分段線性方法。

        1.2 局部均值解調法

        文獻[14]提出了EAD法,指出該方法大多數(shù)情況下較Hilbert法有優(yōu)勢,筆者通過理論分析及數(shù)值仿真發(fā)現(xiàn),在ISC分量的解調中,有必要對EAD法進行改進。

        EAD法的基本思路是:將信號的所有極值點取絕對值,用三次樣條函數(shù)插值形成包絡函數(shù),用原始信號除以包絡函數(shù)得到標準化信號;若標準化信號為調頻分量,則分解結束,否則將標準化信號作為原始信號,重復分解過程。通常情況下,迭代2~3次,標準化分解就會結束,信號可以表示為調幅分量(包絡)和調頻分量的乘積。

        在上述過程中,有兩點需要注意:

        (1)EAD法認為信號是零均值對稱的,所以將所有極值點取絕對值,進而求包絡函數(shù)。但在實際應用中,LCD分解得到的分量往往不能滿足EAD法的假設,只能“局部”滿足“零均值對稱”,這一點從均值點的定義及均值曲線的構造過程可以看出。為了從數(shù)值上說明這個問題,圖1給出了調幅-調頻信號x(t)的均值曲線。x(t)為

        x(t)=0.3(1+0.5 sin5πt)sin(120πt+20πt2)

        圖1 調幅調頻信號的均值曲線

        由圖1可以看出,LCD方法定義的均值曲線是零均值的近似,在局部范圍內與零值存在較大的偏差,按照均值曲線迭代得到的ISC分量并非是零均值對稱的,而EAD法按照全局零均值處理求包絡理論上存在誤差,利用局部信息則能減小其影響。因此,在構造包絡函數(shù)時,考慮用“局部均值曲線”代替“全局均值曲線”。

        (2)EAD法采用三次樣條插值方法求包絡函數(shù),存在“過包絡”、“欠包絡”的問題,從而導致解調精度下降。為了避免這種情況,考慮用滑動平均法代替三次樣條插值。

        基于以上兩點考慮,參考局域均值分解(local mean decomposition,LMD)[15]的思想,提出了一種局部均值解調(local mean envelope,LME)法,將信號的調幅部分和調頻部分分離。與EAD法相比,該方法在包絡函數(shù)的構造方法上作了改進,其余步驟相同。包絡函數(shù)的構造方法為:計算信號的局部均值點,采用滑動平均法求局部均值曲線,用原信號減去局部均值曲線,得到零均值對稱曲線,用滑動平均法計算包絡函數(shù)。計算過程如下。

        (1)確定第j個ISC分量ISCj(t)的所有局部極值點ni及其時刻tni,計算均值點mi、局部幅值ai:

        (4)

        (5)

        (2)用滑動平均法得到均值函數(shù)m11(t)、包絡估計函數(shù)a11(t)。

        (3)將局部均值函數(shù)m11(t)從原始信號ISCj(t)中分離出來,得到

        h11(t)=ISCj(t)-m11(t)

        (6)

        (4)用h11(t)除以包絡估計函數(shù)對h11(t)進行解調,得到

        s11(t)=h11(t)/a11(t)

        (7)

        (5)重復步驟(1)~(4)k次,直至得到一個純調頻信號s1k(t)(瞬時幅值均小于1)及k個包絡信號a1k(t),按照下式計算包絡信號:

        (8)

        從LME解調法的流程可以看出,其實質就是將一個ISC分量視為一個乘積函數(shù)(product function,PF)分量和一個剩余信號之和,然后按LMD方法迭代的第一步求出一個純調頻信號和一個包絡信號,剩余信號視為計算誤差舍去。ISC分量與PF分量的定義內涵類似,剩余信號很小,可作為誤差舍去。

        1.3 瞬時頻率計算

        對純調頻信號s1k(t)利用反正切函數(shù)計算相位:

        (9)

        將φ(t)展開并求導,可得到瞬時頻率。

        求導計算會導致局部極值點附近的瞬時頻率出現(xiàn)畸點,可對瞬時頻率作平滑或濾波處理。本文采用滑動平均法,對任何一個瞬時頻率的計算結果取連續(xù)3個采樣點的平均值。

        1.4 基于瞬時頻率頻譜的齒輪故障診斷

        當齒輪存在故障時,其振動信號會出現(xiàn)調幅、調頻現(xiàn)象。忽略傳遞函數(shù)的影響,齒輪故障振動信號可表示為[8]:

        (10)

        am(t)=Xm(1+dm(t))

        (11)

        Φm(t)=2πmzfrt+φm+bm(t)

        (12)

        式中,fr為軸的轉頻;z為齒輪的齒數(shù);Xm為第m階嚙合頻率諧波分量的幅值;φm為第m階嚙合頻率諧波分量的初相位;dm和bm(t)為第m階嚙合頻率諧波分量的幅值和相位調制函數(shù),兩者都是以fr及其倍頻為重復頻率的周期函數(shù)。

        根據(jù)瞬時頻率fm(t)的定義:

        (13)

        可以看出fm(t)可反映bm(t)的信息。

        應用LCD方法對y(t)進行分解,對ISC分量求出瞬時頻率fm(t)。fm(t)是一個以齒輪轉軸轉頻及其倍頻為中心的頻率分量。對fm(t)進行頻譜分析進而得到瞬時頻率譜,由瞬時頻率譜可以直觀地判斷fr及其倍頻是否存在,從而進行故障診斷。

        2 仿真數(shù)據(jù)分析

        考察下式所示的仿真信號:

        (14)

        其中,x(t)為仿真信號,由調幅調頻信號x1(t)、正弦信號x2(t)合成;采樣頻率為fs=1024 Hz,仿真時間t∈[0,1]。對仿真信號進行LCD,得到2個ISC分量和一個剩余信號r(t),仿真信號及分解結果見圖2。

        圖2 仿真信號及LCD分解結果

        由圖2可以看出,LCD能將調頻調幅信號、正弦信號分解出來,剩余信號幅值很小,具有良好的分解能力。

        為了比較Hilbert法、EAD法和LME法的效果,分別采用三種方法對兩個ISC分量進行包絡解調,將瞬時幅值直接給出。對EAD法和LME法采用反正切法計算瞬時頻率,三種方法的瞬時頻率計算結果均采用三點平滑處理。

        圖3、圖4所示是兩個ISC分量瞬時幅值的計算結果,可以看出三種方法解調結果均能反映原始信號的變化趨勢,其中,Hilbert法計算結果波動最大且端點處的計算值嚴重偏離理論值;EAD法效果優(yōu)于Hilbert法;LME法求得的瞬時幅值最貼近理論值,端點效應最小,更符合原始信號的實際特征。

        圖3 第1個ISC分量ISC1(t)的瞬時幅值

        圖4 第2個ISC分量ISC2(t)的瞬時幅值

        圖5 第1個ISC分量ISC1(t)的瞬時頻率

        圖6 第2個ISC分量ISC2(t)的瞬時頻率

        圖5、圖6所示是兩個ISC分量瞬時頻率的計算結果,可以看出不同方法計算結果差別較大:①在ISC1(t)的計算結果中,Hilbert法端點效應最明顯,兩端的計算誤差已超過實際頻率;EAD法的端點效應較Hilbert法有明顯改善,但中間段的精確度有所下降;LME法計算結果最好;②在ISC2(t)的計算結果中,Hilbert法端點效應最小,EAD法端點效應最大;而中間段的計算效果,LME法最好,Hilbert法與EAD法相當。

        總的來說,在瞬時幅值的計算方面,三種方法計算誤差均不大,LME法優(yōu)于EAD法,Hilbert法效果最差;在瞬時頻率的計算方面,計算誤差均較大,從端點和中間段的綜合效果來看,LME法最優(yōu),Hilbert法最差。

        3 試驗數(shù)據(jù)分析

        圖7 齒輪箱傳動結構

        齒輪故障試驗數(shù)據(jù)源于QPZZ-Ⅱ旋轉機械振動分析及故障診斷試驗平臺系統(tǒng),其齒輪箱傳動結構如圖7所示。變頻調速電機通過聯(lián)軸節(jié)驅動小齒輪,大齒輪與小齒輪直接嚙合。大小齒輪均為圓柱齒輪,大齒輪齒數(shù)為75,小齒輪齒數(shù)為55。

        人為地將小齒輪的一個齒尖切割約5 mm,模擬小齒輪斷齒中度故障。試驗中,大齒輪輸出軸負載為零。設置電機軸轉速為880 r/min,實測轉速878 r/min,則小齒輪的轉頻f1=14.6 Hz,大齒輪的轉頻f2=10.7 Hz,齒輪嚙合頻率為804.8 Hz。采用加速度傳感器采集振動信號,傳感器在輸出軸電機側軸承處垂直于齒輪箱上表面安裝,信號采樣頻率為fs=5120 Hz,計算數(shù)據(jù)點N=1024。

        (a)時域波形

        圖8給出了原始信號的時域波形圖及包絡譜。從時域波形圖上可以看出原始信號具有調幅-調頻特性,并且存在周期性沖擊。為了驗證LCD本身的能力,本文沒有采用中值濾波、SVD等降噪技術,包絡譜直接采用Hilbert法求出,為了觀察特征頻段的特征,只給出了低頻段(0~1000 Hz)的波形。從包絡譜上可以看出20 Hz、45 Hz處有明顯譜峰。20 Hz可近似認為是大齒輪轉頻的2倍頻,45 Hz可近似認為是小齒輪轉頻的3倍頻,由此可以判定齒輪箱出現(xiàn)了故障,但不能判定哪個齒輪出現(xiàn)了故障。

        (b)包絡譜

        對原始采樣信號進行LCD,得到4個ISC分量和1個剩余信號。圖9是4個ISC分量的時域波形??梢钥闯觯孩貺CD類似于自適應濾波器,4個ISC分量所包含的頻率段逐漸降低;②每個ISC分量的幅值依次減?。虎勖總€ISC分量均可近似認為是調幅-調頻分量。因此,LCD是一種有效的信號自適應分解方法。

        圖9 原始信號LCD的ISC分量時域波形

        圖10 應用Hilbert法解調的ISC分量瞬時頻率譜

        圖11 應用EAD法解調的ISC分量瞬時頻率譜

        圖12 應用LME法解調的ISC分量瞬時頻率譜

        對4個ISC分量分別應用三種方法解調,并求瞬時頻率譜,將低頻段(0~1000 Hz)的結果列于圖10~圖12??梢钥闯觯孩僭诜至?的瞬時頻率譜中,小齒輪轉頻3倍頻處存在相對明顯的譜線;②在分量2、3的瞬時頻率譜中,小齒輪轉頻2倍頻處存在相對明顯的譜線;③分量4的瞬時頻率譜能清晰反映小齒輪轉頻及其2倍頻。因此, LCD結合三種解調方法求出的瞬時頻率譜均能得出小齒輪存在故障的正確結論。

        但是Hilbert法求得的瞬時頻率譜中存在明顯的高頻能量泄漏現(xiàn)象,其他頻段的譜線雜亂,干擾較大;在EAD法求得的瞬時頻率譜中,高頻段的干擾得到明顯的抑制;而LME法求得的瞬時頻率譜高頻干擾最小,信噪比最高。因此,LME法在齒輪故障診斷中具有一定的優(yōu)越性。

        4 結論

        本文根據(jù)齒輪故障振動信號的多分量調幅-調頻特征,提出了局部特征尺度分解(LCD)和局部均值解調(LME)相結合的故障診斷方法。通過仿真信號驗證了LME解調的精確性;通過齒輪斷齒故障振動數(shù)據(jù)的分析,驗證了LCD結合LME的故障診斷方法的有效性,該方法不僅能夠診斷出齒輪箱故障與否,還能準確地定位發(fā)生故障的齒輪,具有一定的工程應用價值。

        [1] 康海英, 祁彥潔, 閻文, 等. 運用階次跟蹤和奇異譜降噪診斷齒輪早期故障[J]. 振動測試與診斷, 2010,30(6):662-665. Kang Haiying, Qi Yanjie, Yan Wen, et al. Early Fault Diagnosis and Singular of Gear Spectrum by Order Tracking Denoising[J]. Journal of Vibration Measurement & Diagnosis, 2010,30(6): 662-665.

        [2] 趙軍, 崔穎, 賴欣歡, 等. 隨機共振降噪下的齒輪微弱故障特征提取[J]. 中國機械工程, 2014, 25(4): 539-546. Zhao Jun, Cui Ying, Lai Xinhuan, et al. Weak Texture Extraction of Gear Faults Based on Stochastic Resonance Denoising[J]. China Mechanical Engineering, 2014, 25(4): 539-546.

        [3] Huang N E, Shen Z, Long R S, et al. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-stationary Time Series Analysis[J]. Proc. Roy. Soc. London, 1998, 454: 903-995.

        [4] Huang N E, Shen Z, Long S R. A New View of Nonlinear Water Waves: the Hilbert Spectrum[J]. Annual Review of Fluid Mechanics, 1999, 31(1): 417-457.

        [5] Huang N E, Chiang W, Busalacchi A J. A Bridge Health Monitoring Method Using the Hilbert-Huang Transform: a Case Study [R]. Washington: Nasa, 2000.

        [6] Huang N E,Wu Z H. A Review on Hilbert-Huang Transform: Method and Its Applications to Geophysical Studies[J]. Adv. Adapt. Data Anal., 2009, 1: 1-23.

        [7] 于德介, 程軍圣, 楊宇. 機械故障診斷的Hilbert-Huang變換方法[M]. 北京: 科學出版社. 2006.

        [8] 于德介, 程軍圣, 楊宇. Hilbert-Huang變換在齒輪故障診斷中的應用[J]. 機械工程學報, 2005, 41(6): 102-107. Yi Dejie, Cheng Junsheng, Yang Yu. Application of Hilbert-Huang Transform Method to Gear Fault Diagnosis[J]. Chinese Journal of Mechanical Engineering,2005,41(6):102-107.

        [9] 程軍圣, 鄭近德, 楊宇. 一種新的非平穩(wěn)信號分析方法—局部特征尺度分解法[J]. 振動工程學報, 2012, 25(2): 215-220. Cheng Junsheng,Zheng Jinde,Yang Yu. A New Non-Stationary Signal Analysis Approach—The Local Characteristic-Scale Decomposition Method[J]. Journal of Vibration Engineering,2012,25(2):215-220.

        [10] Frei M G,Osorio I. Intrinsic Time-scale Decomposition: Time-frequency-energy Analysis and Real-Time Filtering of Non-stationary Signals[J]. Proc. Royal Soc. A, 2007, 463: 321-342.

        [11] 楊宇, 曾鳴, 程軍圣. 局部特征尺度分解方法及其分解能力研究[J]. 振動工程學報, 2012, 25(5): 602 -609. Yang Yu, Zeng Ming, Cheng Junsheng. Research on Local Characteristic-scale Decomposition and Its Capacities[J].Journal of Vibration Engineering,2012,25(2):602-609.

        [12] 楊宇, 曾鳴, 程軍圣. 局部均值尺度分解方法及其分量判據(jù)研究[J]. 中國機械工程, 2013, 24(2): 195-208. Yang Yu, Zeng Ming, Cheng Junsheng. Research on Local Characteristic-scale Decomposition and Its Stopping Criteria[J]. China Mechanical Engineering, 2013, 24(2): 195-208.

        [13] 戴豪民, 許愛強. 瞬時頻率計算方法的比較研究和改進[J]. 四川大學學報(自然科學版), 2014, 51 (6): 1197-1203. Dai Haomin, Xu Aiqiang. Comparative Research and Improvement on the Calculation Method of Instantaneous Frequency[J]. Journal of Sichuan University(Natural Science Edition), 2014,51(6):1197-1203.

        [14] Huang N E. Computing Instantaneous Frequency by Normalizing Hilbert Transform: U.S., 6901353[P]. 2005-05-31.

        [15] Smith J S. The Local Mean Decomposition and Its Application to EEG Perception Data[J]. Journal of the Royal Society Interface, 2005, 2(5): 443-454.

        (編輯 蘇衛(wèi)國)

        Gear Fault Diagnosis Based on LCD and LME Demodulation Approach

        Cui Weicheng Xu Aiqiang Li Wei Meng Fanlei

        Naval Aeronautical and Astronautical University,Yantai,Shandong,264001

        In order to diagnose the gear fault accurately, according to the amplitude-modulated and frequency-modulated characteristics of the gear fault vibration signals a modle was proposed based on LCD and LME. Firstly,the original vibration signals were disposed by LCD.Secondly, the intrinsic mode components of LCD were envelope demodulated by LME. Finally, the gear faults were diagnosed by the instantaneous frequency spectrum. The analysis of the simulated data shows that the LME mode may get more accurate results of instantaneous amplitude and instantaneous frequency than that of Hilbert mode and empirical AM-FM decomposition mode. The analysis of gear fault data shows that the method based on LCD and LME may realize the fault diagnosis effectively.

        local characteristic-scale decomposition(LCD); local mean envelope(LME) demodulation; instantaneous frequency spectrum; gear fault diagnosis

        2015-12-18

        國家部委預研基金資助項目(9140A27020214JB1446)

        TN911.23;TP206.3

        10.3969/j.issn.1004-132X.2016.24.012

        崔偉成,男,1981年生。海軍航空工程學院飛行器工程系講師、博士研究生。研究方向為裝備智能故障診斷。發(fā)表論文13篇。許愛強,男,1963年生。海軍航空工程學院飛行器檢測與應用研究所教授。李 偉,男,1980年生。海軍航空工程學院飛行器工程系講師。孟凡磊,男,1978年生。海軍航空工程學院飛行器工程系講師。

        猜你喜歡
        調頻均值齒輪
        東升齒輪
        內燃機工程(2021年6期)2021-12-10 08:07:46
        考慮頻率二次跌落抑制的風火聯(lián)合一次調頻控制
        能源工程(2021年5期)2021-11-20 05:50:42
        你找到齒輪了嗎?
        異性齒輪大賞
        齒輪傳動
        均值不等式失效時的解決方法
        均值與方差在生活中的應用
        調頻發(fā)射機技術改造
        調頻激勵器干擾的排除方法
        關于均值有界變差函數(shù)的重要不等式
        一区二区视频网站在线观看| 亚洲av无码av制服丝袜在线| 91日韩高清在线观看播放| 亚洲人成无码网站十八禁| 福利视频偷拍一区二区| 99麻豆久久久国产精品免费| 亚洲精品国产福利一二区| 中字亚洲国产精品一区二区| 中文字幕人妻激情在线视频| 国内永久福利在线视频图片| 欧洲熟妇色xxxx欧美老妇多毛图片| 国产小屁孩cao大人| 国产91大片在线观看| 久久国语露脸国产精品电影| 国产熟人av一二三区| 99在线视频精品费观看视| 亚洲精品中文字幕码专区| 亚洲国产精品成人久久久| 欧美亚洲日本国产综合在线| 色系免费一区二区三区| va精品人妻一区二区三区| 高潮抽搐潮喷毛片在线播放| 亚洲av有码在线天堂| 国产日韩午夜视频在线观看| 白白色发布免费手机在线视频观看| 日韩一区国产二区欧美三区 | 欧美日韩中文国产一区发布| 一区二区在线亚洲av蜜桃| 日韩精品国产精品亚洲毛片| 无码人妻一区二区三区兔费 | 国产一区二区三区尤物| 99久久99久久久精品齐齐| 天天躁人人躁人人躁狂躁 | 国产韩国一区二区三区| 少妇粉嫩小泬喷水视频| 伊人婷婷在线| 一区二区亚洲精美视频| 丰满熟妇乱又伦精品| 亚洲欧洲日产国码无码久久99| 亚洲福利第一页在线观看| 国产一区二区三区毛片|