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

        ?

        抗環(huán)境音干擾的設(shè)備聲音故障監(jiān)測方法①

        2019-07-23 02:07:52李蘭村廉東本毛立爽
        計算機系統(tǒng)應(yīng)用 2019年6期
        關(guān)鍵詞:降維聚類維度

        李蘭村,廉東本,毛立爽

        1(中國科學院大學,北京 100049)

        2(中國科學院 沈陽計算技術(shù)研究所,沈陽 110168)

        1 引言

        變壓器等大型設(shè)備在運行時的聲音具有一定的辨識性,目前普遍通過擁有多年設(shè)備檢修經(jīng)驗的技術(shù)人員對變壓器運行聲音進行監(jiān)聽,并判斷設(shè)備是否正常運行.但該方式有兩個不足之處,一是需要技術(shù)人員時刻在設(shè)備前監(jiān)測設(shè)備運行聲音,若在惡劣環(huán)境條件下安全難以保障;二是人耳的聽覺具有特定的聽覺感知范圍,只能感知特定頻率范圍和特定響度范圍內(nèi)的聲音[1].近年來,多個領(lǐng)域借鑒語音識別技術(shù)取得了豐富的研究成果.通過在變壓器旁安裝拾音器采集聲音并對其進行實時監(jiān)測,不僅不會干擾變壓器的正常運行,而且裝置安裝簡單,采集信號方便.但在實際環(huán)境中存在復(fù)雜環(huán)境噪聲,可能會對采集的聲音信號造成一定影響,從而導(dǎo)致聲音識別的準確率下降.

        現(xiàn)有的研究中通常采用分類算法建立聲音模型,例如文獻[2]中作者提出了一種基于改進的支持向量機分類算法建立聲音模型的解決方案.但分類算法只能在分類個數(shù)已知的條件下建立聲音模型,對于由復(fù)雜環(huán)境音而導(dǎo)致分類個數(shù)不確定的條件下建立聲音模型的情形并不適用.此外,現(xiàn)有研究對建立的聲音模型采用剝離環(huán)境音的方式檢測設(shè)備是否發(fā)生故障,該種檢測方式在環(huán)境音影響較小的條件下對聲音模型的去噪效果理想,但在多種復(fù)雜環(huán)境音的影響下,對聲音模型的去噪將變得困難且效果不理想,從而導(dǎo)致較高的誤警率,并降低整體的檢測效率.本文提出的方法是對保留環(huán)境音的設(shè)備工作聲音數(shù)據(jù)進行聚類,并對得到的訓練模型進行迭代優(yōu)化.

        2 模型構(gòu)建

        本文采用的模型主要由數(shù)據(jù)預(yù)處理模塊、特征提取模塊、特征降維模塊、特征訓練模塊、標準集增強模塊構(gòu)成.模型數(shù)據(jù)流程圖如圖1所示.

        原始聲音數(shù)據(jù)經(jīng)過預(yù)處理、特征提取和降維后分成兩部分,一部分用于訓練標準集,另一部分用于匹配增強標準集.

        2.1 數(shù)據(jù)預(yù)處理模塊

        采集運行中的設(shè)備聲音時,由于拾音器參數(shù)不同或拾音器所放置的方向和距離的不同會導(dǎo)致采集的聲音品質(zhì)不同,因此需要對采集到的聲音進行歸一化處理,降低不同樣本間因采集因素造成的差異.

        圖1 模型數(shù)據(jù)流程圖

        處理聲音信號時,需使用快速傅里葉變換(FFT)明確聲音中各個頻率成分的分布.傅里葉變換要求輸入的信號是平穩(wěn)的,雖然聲音信號從時間軸上看具有非平穩(wěn)特性,但在極短的時間內(nèi),則呈現(xiàn)平穩(wěn)信號所具有的特性.因此,將一段聲音信號切分成若干小段,就能將其近似地認為是平穩(wěn)的信號,且該小段被稱為幀長.但聲音信號具有連續(xù)性和相關(guān)性,因此不能講聲音直接均分.且為了提高傅里葉變換的分辨率,在變換前需要對每一幀信號進行”加窗”,即將其與一個窗函數(shù)相乘,使得該幀信號的兩端逐漸趨于零.常見的窗函數(shù)有漢寧窗和矩形窗,但漢寧窗的旁瓣衰減較大,具有更平滑的低通特性,能夠在較高程度上反應(yīng)短時信號的頻率特性.而矩形窗的旁瓣峰值較大,導(dǎo)致頻譜泄露較為嚴重[2],因此本文采用漢寧窗函數(shù).由于加窗操作會將一幀信號的兩端趨于0,所以需要將截取的后一幀向前一幀移動一段距離,稱為幀移,一般重疊距離為幀長的1/3 至1/2.漢寧窗函數(shù)的公式如下:

        式(1)中,a一般取0.46.

        2.2 特征提取模塊

        通過對聲音信號特征的提取,分析和處理,最終為之后的處理得到主要的聲音特征數(shù)據(jù).選取最優(yōu)的聲音特征參數(shù)對于聲音特性的識別監(jiān)測至關(guān)重要,只有根據(jù)盡量降低計算復(fù)雜度和降低特征維數(shù)的要求來選擇不同的聲音特征參數(shù),才能保證故障監(jiān)測系統(tǒng)的識別正確率和效率.在非語音類聲音的模式識別領(lǐng)域,因為聲音信號在時域上變化快速且不穩(wěn)定,所以通常將其轉(zhuǎn)換到頻域上分析特征參數(shù).倒頻譜特征參數(shù)因為具有分離頻譜上的高低頻的優(yōu)點而被廣泛使用[1],常用的倒頻譜特征參數(shù)有梅爾頻率倒譜系數(shù)(MFCC)、LPC 倒譜系數(shù)(LPCC)、線性預(yù)測倒譜系數(shù)(LPC)等.由于梅爾頻率倒譜系數(shù)是根據(jù)人耳頻率的非線性特性而提取出來的倒譜參數(shù),梅爾頻率濾波器組模擬人的耳蝸聲道模型,實現(xiàn)了模擬人類特定的聽覺掩蔽和聽覺修復(fù)等功能,所以更加適用于需要人耳判別聲音的場景,本文選用提取MFCC 來轉(zhuǎn)化聲音數(shù)據(jù),梅爾標度與頻率的關(guān)系如式(2)所示:

        式中f代表頻率.

        提取MFCC 系數(shù)首先要將預(yù)處理后的聲音信號進行FFT,計算出每幀數(shù)據(jù)的頻譜參數(shù),以得到在頻譜上的能量分布.其次,設(shè)置M(24)個帶通濾波器,每個濾波器具有三角形濾波器的特性,可以對頻譜起到平滑化且消除諧波的作用.然后,計算每個濾波器組輸出的對數(shù)能量.最后,將對數(shù)能量經(jīng)過離散余弦變換,求出1×12 的MFCC 特征系數(shù)矢量.根據(jù)MFCC 提取聲音特征矩陣流程如圖2所示.

        圖2 提取MFCC 特征矩陣流程

        2.3 特征降維模塊

        提取出的MFCC 聲音特征參數(shù)存在數(shù)據(jù)維度高、數(shù)據(jù)有冗余等缺陷從而導(dǎo)致存儲空間占用過多和計算量過大等問題.通過降維可以舍棄貢獻度不高的維度,從而達到減少計算量和計算時間的目的.此外,提取到的聲音中包含復(fù)雜的環(huán)境音,當數(shù)據(jù)受到噪聲影響時,低貢獻度維度數(shù)據(jù)往往與噪聲有關(guān),因此舍棄這部分數(shù)據(jù)也能減少噪聲對數(shù)據(jù)本身的影響.由于提取到的聲音數(shù)據(jù)特征在不同維度上的方差分布不均勻,因此更適合采用PCA 降維算法.PCA 通過計算聲音數(shù)據(jù)的協(xié)方差矩陣及其特征值分解得到特征值和相應(yīng)的特征向量劃分數(shù)據(jù)的主成分.PCA 降維算法的目標是去掉不同數(shù)據(jù)維度之間的相關(guān)性,并盡量使降維后的數(shù)據(jù)維度能夠很好地保持原始數(shù)據(jù)的特征.

        PCA 降維首先要對聲音數(shù)據(jù)標準化,使得每個維度的數(shù)據(jù)均落入相同的區(qū)間;其次,獲取聲音特征矩陣的協(xié)方差矩陣;然后對協(xié)方差矩陣進行特征值分解,獲取各維度的特征值和對應(yīng)的特征向量;最后,將特征向值降序排列后通過維度選取算法獲取降維幅度.

        在降維的過程中,對降維幅度的選取至關(guān)重要.在現(xiàn)有的研究中,普遍采用交叉驗證或逆向重構(gòu)兩種方法來獲取最優(yōu)的降維幅度.為避免額外開銷(如學習器等),本文采用逆向重構(gòu)的方法獲取最優(yōu)降維幅度.根據(jù)實際生產(chǎn)經(jīng)驗及研究,本文設(shè)定重構(gòu)閾值t的值為95%[3],然后選取使得式(3)成立的最小維度.

        其中,d為原始數(shù)據(jù)的維度,d’為降維后的維度,λi為該維度特征向量的特征值.

        2.4 特征訓練模塊

        標準集的設(shè)計需要能使訓練數(shù)據(jù)獲取最優(yōu)分類,由于設(shè)備聲音在運行過程中會與環(huán)境音混合,設(shè)備正常工作聲音的種類數(shù)量及聲音類別標號均未知,因此需要提前確定類別標號的有監(jiān)督分類學習算法是不可行的[4-6],綜上,本文采用無監(jiān)督的聚類學習算法.

        聚類處理將數(shù)據(jù)對象集劃分成多個簇,使得簇內(nèi)的對象具有較高的相似性,且與其他簇中對象的相似性較低.通過聚類算法能發(fā)現(xiàn)任意形狀的簇,且具有處理噪聲點數(shù)據(jù)的能力,這些優(yōu)點在變壓器聲音數(shù)據(jù)進行模式匹配時具有輔助作用.在基于密度的聚類算法中比較常見的有DBSCAN(Density-Based Spatial Clustering of Application with Noise)和OPTICS(Odering Point To Identify the Cluster Structure),兩者較為相似,均通過核心對象及其鄰域發(fā)現(xiàn)稠密區(qū)域,該稠密區(qū)域被稱為簇.

        DNSCAN 算法中有兩個重要的參數(shù):分別為Eps(定義密度時的鄰域半徑)和MinPts(定義核心點時的閾值),本文將其簡化為ε和M.對于數(shù)據(jù)集合X={x(1),x(2),…,x(N)},引入以下概念:

        (1)鄰域

        設(shè)x∈X,則

        為x的ε鄰域.

        (2)密度

        設(shè)x∈X,則

        為x的密度.

        (3)核心點

        設(shè)x∈X,若ρ(x)≥Μ,則x稱為X的核心點,X中所有核心點組成的集合記為Xc,X中所有非核心點組成的集合為Xnc.

        (4)邊界點

        若x∈Xnc,且x落在某個核心點的ε鄰域內(nèi),則稱x為X的一個邊界點,X中所有邊界點組成的集合記為Xbd.

        (5)噪音點

        設(shè)Xnoi=X(Xc∪Xbd),若x∈Xnoi,則稱x為噪音點.

        (6)直接密度可達

        設(shè)x,y∈X,若x∈Xc,y∈Nε(x),則稱y由x直接密度可達.

        (7)密度可達

        設(shè)P(1),P(2),…,P(n)∈X(n≥2),若滿足P(i+1)是從P(i)直接密度可達的,則稱P(n)是從P(1)密度可達的.

        (8)密度相連

        設(shè)x,y,z∈X,若y,z都可由x密度可達,則稱y,z密度相連.

        DBSCAN 算法的核心思想是:從某個隨機選定的核心點開始,不斷向密度可達的區(qū)域擴張,最終得到一個包含核心點和邊界點的簇,使得簇中任意兩點密度相連.由于DBSCAN 算法初始時設(shè)定了兩個固定的輸入?yún)?shù)ε和Μ,這需要提前根據(jù)經(jīng)驗來設(shè)定,設(shè)置的細微不同可能造成聚類結(jié)果的巨大差異.Μ的設(shè)定可以根據(jù)一個原則,即Μ≥m+1(m表示數(shù)據(jù)維度),但ε的設(shè)定是比較困難的,因為對密度差異較大的數(shù)據(jù)很難界定一個密度,使得邊界點不被遺落且有差異的簇不被歸到一個簇中.

        OPTICS 算法不需要設(shè)置全局參數(shù),而是為聚類分析生成一個增廣簇排序,它代表了各樣本點基于密度的聚類結(jié)構(gòu),從該排序中可以獲得基于任何半徑和任何閾值的聚類[7].OPTICS 算法是對DBSCAN 算法的一種改進,因此上文提到的關(guān)于DBSCAN 算法中的概念在OPTICS 算法中依然適用,但增加了以下兩個概念.

        (9)核心距離

        設(shè)x∈X,對于給定的參數(shù)ε和Μ,稱使得x成為核心點的最小鄰域半徑為x的核心距離,記為cd(x).

        (10)可達距離

        設(shè)x,y∈X,若Nε(x)<M,則

        OPTICS 算法過程如下:

        算法1.OPTICS 算法輸入:數(shù)據(jù)樣本X,初始化所有點的可達距離與核心距離均為MAX,并初始化鄰域半徑ε和閾值Μ.1)建立兩個隊列,有序隊列包含核心點及該核心點的直接密度可達點,結(jié)果隊列包含存儲樣本輸出及處理次序;2)若X中數(shù)據(jù)全部處理完畢,則轉(zhuǎn)到步驟7),否則從X中選擇一個為核心對象且未處理過的點,將此核心點放入結(jié)果隊列,將此核心點的直接密度可達點按可達距離升序排列放入有序隊列;3)若有序序列為空,回到步驟2),否則從有序隊列中取出第一個點;4)若該點不是核心點則回到步驟3),若該點是核心點是且該點不在結(jié)果隊列中則將該點放入結(jié)果隊列;5)若該點是核心點,找到該點所有直接密度可達點,并將這些點放入有序隊列中,且將有序隊列中的點按照可達距離重新排序,若該點已經(jīng)在有序隊列中且新的可達距離更小,則更新該點的可達距離.6)重復(fù)步驟3)、4)、5),直至有序隊列為空.7)算法結(jié)束.

        至此,通過算法獲取到一個有序數(shù)組和每個元素的核心距離和可達距離,然后利用這些數(shù)據(jù)生成一個聚類結(jié)果,加入一個鄰域半徑參數(shù)ξ,ξ<ε.由于需要將數(shù)據(jù)集合分成k個簇及噪音點集合,因此引入簇標記數(shù)組如下:

        2.5 標準集增強模塊

        將訓練樣本數(shù)據(jù)通過特征訓練模塊訓練后得到了具有不同密度的簇,每一簇代表一種噪聲混合變壓器正常工作聲音或僅為變壓器正常工作聲音.這些簇均屬于變壓器在正常工作環(huán)境下會形成的簇,但簇卻是不完整的.由于環(huán)境干擾音種類繁多,在前期無法把各種環(huán)境干擾音全部采集到,若在對聲音特征進行匹配識別時,出現(xiàn)了沒有經(jīng)過訓練的環(huán)境干擾音,則會出現(xiàn)匹配失敗的情況,但該聲音的確是正常的聲音,因此為了避免多次出現(xiàn)類似的情況,本文引入標準集增強模塊,標準集增強模塊能將新出現(xiàn)的正常聲音加入到標準集當中.

        標準集增強模塊的流程如下:

        (1)在已經(jīng)存在標準集的情況下,假設(shè)已知參考集為C={C1,C2,…,Cn},其中每一類代表一個簇,假設(shè)增加的測試樣本數(shù)據(jù)為T={t(1),t(2),…,t(s)},其中s表示樣本數(shù)據(jù)的幀數(shù).當前的測試樣本數(shù)據(jù)中既存在非正常數(shù)據(jù)幀也存在提前標記完成的正常數(shù)據(jù)幀;

        (2)將T中的每幀與C中每個簇按照式(8)進行歐式距離的運算.

        式(8)中,i表示樣本幀號,j表示參考集簇的序號,dij表示第i幀樣本數(shù)據(jù)與第j個簇之間的距離,得到每幀與各個簇之間的距離為:di1,di2,…,din,其中i為幀號,n為簇的個數(shù),取

        若dmin小于提前設(shè)定好的閾值,則標記為正常數(shù)據(jù),否則,標記為異常數(shù)據(jù).

        (3)比較數(shù)據(jù)匹配前后的標記,若正常數(shù)據(jù)幀被標記為非正常數(shù)據(jù),說明該正常數(shù)據(jù)在初始的訓練集中并沒有出現(xiàn)過,便將其進行聚類,形成新的一個簇,加入到C中,此時簇C={C1,C2,…,Cn,Cn+1}.

        (4)至此一次標準集增強處理完畢,若再有新數(shù)據(jù)的加入,則進行增強處理,訓練集就會一直優(yōu)化迭代.

        數(shù)據(jù)經(jīng)過標準集增強模塊后,在標準集中沒有出現(xiàn)過且的確屬于設(shè)備正常工作的聲音種類會被加入到標準集中,進而在下次與樣本聲音數(shù)據(jù)進行匹配時候?qū)⒃擃惵曇糇R別為正常聲音,而不會將其識別為異常聲音進而誤警.本文通過實驗,驗證經(jīng)過標準集增強模塊后的訓練集在與測試數(shù)據(jù)匹配時的誤警率隨著迭代次數(shù)增加而逐漸降低.

        3 實驗分析

        3.1 實驗環(huán)境

        實驗運行在Windows 10 操作系統(tǒng)上,實驗PC 機主頻為3.6 GHz,內(nèi)存4 GB.由于變壓器故障樣本獲取非常困難,為了測試算法對故障情況的識別能力,本文合成了一批故障樣本.本文對樣本進行篩選并選取了有代表性的聲音種類,包含正常聲音種類20 種,非正常的故障聲音3 種,每種聲音截取10 分鐘.訓練樣本從正常聲音數(shù)據(jù)中隨機選擇15 種,測試樣本個數(shù)為40 個,且每一個測試樣本均隨機從所有聲音數(shù)據(jù)中選擇5 種.

        3.2 實驗設(shè)計和結(jié)果

        首先將訓練樣本和參考樣本進行預(yù)處理和特征提取.在分幀階段,設(shè)置幀長為40 ms,幀移為20 ms,故訓練樣本共有4.5×105幀,每個測試樣本有1.5×105幀.在加窗階段對聲音的每一幀乘以漢寧窗函數(shù).對預(yù)處理后的數(shù)據(jù)提取MFCC 特征參數(shù)并將其進行標準化處理后的特征矩陣如圖3所示.

        圖3 MFCC 特征矩陣

        其中特征矩陣共12 列,每一列代表每一幀的一個屬性值,特征矩陣共640000 行,每一行代表初始聲音數(shù)據(jù)中的一幀.

        然后計算該MFCC 特征矩陣的協(xié)方差矩陣,并將該協(xié)方差矩陣進行特征值分解,得到如圖4所示的降序排列的12 個特征值及每個特征值對應(yīng)的特征向量.

        圖4 MFCC 特征矩陣的特征值集合

        根據(jù)式(3)可知,僅用前三個特征值對應(yīng)的特征向量便可使重構(gòu)率達到97.15%,滿足不小于重構(gòu)閾值95%的要求,因此可以使用PCA 算法將原矩陣降維至3 維矩陣,降維后的特征矩陣如圖5所示.

        圖5 降維后的特征矩陣

        對降維后的訓練樣本聲音特征矩陣采用OPTICS算法進行聚類,根據(jù)實際開發(fā)經(jīng)驗及文獻研究,本實驗采用的OPTICS 算法中鄰域半徑ε 的值為0.24,閾值Μ的值為16 時,聚類結(jié)果最好,然后將測試樣本依次與標準集進行匹配.由于變壓器正常工作聲音數(shù)據(jù)與標準集偏差小于0.48,而變壓器故障聲音數(shù)據(jù)與標準集偏差大于1.5,因此為減少模型對噪音點的敏感度,本實驗將匹配閾值設(shè)置為0.5,并記錄下每次匹配增強后的標準集的結(jié)果.分別對原始標準集和經(jīng)過40 次測試樣本迭代增強后的標準集計算正常率和誤警率,得到的結(jié)果如表1所示.

        表1 模型識別正常率和誤警率(%)

        圖6為模型的準確率和誤警率隨著迭代次數(shù)的增加而改變的折線圖.

        3.3 實驗結(jié)果分析

        從表1與圖6可以看出,模型的正常率隨著迭代次數(shù)的增加不斷升高,因為設(shè)備故障聲音數(shù)據(jù)占總樣本數(shù)據(jù)的13%,且每次迭代選擇的樣本數(shù)據(jù)比較大,所以正常率最終只會逼近但實際生產(chǎn)中變壓器發(fā)生故障的情況較少,所占總樣本數(shù)據(jù)的比例不高,因此實際生產(chǎn)中正常率會遠高于當前數(shù)值.誤警率是一個更能反映模型優(yōu)勢的指標,可以看出,迭代10 次后的誤警率比迭代5 次后要高,因為在迭代前期,新出現(xiàn)的正常聲音數(shù)據(jù)較多,出現(xiàn)誤報的次數(shù)就相對較多.當新出現(xiàn)的正常聲音數(shù)據(jù)均存入標準集后,隨著迭代次數(shù)的增加,誤警率將不斷降低.實驗結(jié)果說明,本文提出的模型能準確識別非正常聲音,且隨著迭代次數(shù)的增加,模型的誤警率將持續(xù)降低,可信賴程度將不斷升高.

        圖6 模型識別正常率和誤警率折線圖

        4 結(jié)論與展望

        本文基于變壓器等大型設(shè)備在運行過程中發(fā)聲具有辨識性和平穩(wěn)性,且容易受各種環(huán)境音干擾的特點,提出了一種抗多種環(huán)境音干擾的設(shè)備聲音故障監(jiān)測方法.一方面,將樣本數(shù)據(jù)經(jīng)過特征提取與OPTICS 聚類等模塊得到初始標準集,另一方面,將測試樣本與已有標準集不斷匹配和增強獲得更優(yōu)標準集.通過實驗,驗證了本文提出的方法具有可行性和有效性,可以提高設(shè)備狀態(tài)監(jiān)測數(shù)據(jù)的可靠性,進而有助于更快推進變壓器等大型設(shè)備自動監(jiān)測的進程.

        猜你喜歡
        降維聚類維度
        Three-Body’s epic scale and fiercely guarded fanbase present challenges to adaptations
        降維打擊
        海峽姐妹(2019年12期)2020-01-14 03:24:40
        淺論詩中“史”識的四個維度
        中華詩詞(2019年7期)2019-11-25 01:43:00
        基于DBSACN聚類算法的XML文檔聚類
        電子測試(2017年15期)2017-12-18 07:19:27
        光的維度
        燈與照明(2016年4期)2016-06-05 09:01:45
        “五個維度”解有機化學推斷題
        基于改進的遺傳算法的模糊聚類算法
        一種層次初始的聚類個數(shù)自適應(yīng)的聚類方法研究
        拋物化Navier-Stokes方程的降維仿真模型
        計算物理(2014年1期)2014-03-11 17:00:18
        基于特征聯(lián)合和偏最小二乘降維的手勢識別
        97无码免费人妻超级碰碰夜夜| 国产一区二区不卡av| 一本大道在线一久道一区二区| 99国产精品欲av麻豆在线观看| 高清日韩av在线免费观看| 国产精品久久久久乳精品爆| 免费无码肉片在线观看| 国产精品久久一区性色a| 亚洲av第一区国产精品| 18禁黄网站禁片免费观看女女| 人人妻人人玩人人澡人人爽| 偷拍熟女亚洲另类| 亚洲乱妇熟女爽到高潮视频高清| 中文字幕网伦射乱中文| 中文字幕+乱码+中文字幕无忧| 久久久久无码精品国| 久久伊人亚洲精品视频 | 国产a级毛片久久久精品毛片| 亚洲日韩中文字幕一区| 国产亚洲第一精品| 丝袜美腿一区在线观看| 亚洲va韩国va欧美va| 国产精品久久久av久久久| 一区二区三区国产97| 91精品久久久老熟女91精品| 久久精品成人一区二区三区| 久久精品国产亚洲av蜜臀| 精品亚洲少妇一区二区三区| 亚洲午夜精品第一区二区| 成人亚洲一区二区三区在线| 精品欧美一区二区在线观看 | 国产精品久久久久久妇女6080| 四虎成人精品国产永久免费| 日韩精品视频高清在线| 琪琪的色原网站| 国产九色AV刺激露脸对白| 中文字幕熟女激情50路| 777米奇色狠狠俺去啦| 91视频88av| 国产精品一区又黄又粗又猛又爽 | 久久精品这里只有精品|