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

        ?

        一元及多元信號分解發(fā)展歷程與展望

        2024-02-05 12:38:08陳啟明文青松蘇宏業(yè)
        自動化學報 2024年1期
        關鍵詞:時頻時變模態(tài)

        陳啟明 文青松 郎 恂 謝 磊 蘇宏業(yè)

        在一個飛速發(fā)展的信息社會中,信號是無處不在的,例如人們隨時可以聽到的語音信號、隨時可以看到的視頻信號、伴隨生命始終存在的生理信號、新冠肺炎感染患者的數量、工業(yè)生產中的控制監(jiān)測信號、黑洞碰撞的引力波信號等[1].這些信號是變化的,變化的信號構成了五彩斑斕的世界,如何描述這些變化的信號以及揭示這些信號中潛在的信息,是信號處理研究的重要任務[1].一般而言,在傳統(tǒng)的信號處理中,人們會從時域或頻域兩個角度去描述信號.但單純從時域或頻域角度出發(fā)描述信號較為片面,無法刻畫非平穩(wěn)和非線性信號的重要性質.在現實世界中,所獲得的信號大部分都是非平穩(wěn)和非線性的,為了滿足現實信號處理需求,時頻聯(lián)合分析方法相關研究應運而生[2].

        1998 年,Huang 等[3]提出一種自適應的非平穩(wěn)非線性信號分解方法,成為信號處理發(fā)展的一個重要里程碑,隨后很快發(fā)展為“后小波”時代時頻分析領域的研究熱點[4].信號分解假設復雜的非平穩(wěn)非線性信號由多個簡單的子信號組成,通過分析這些子信號的特征,可以間接或直接揭示原始復雜信號的時頻信息,進而為各個領域的信號處理任務提供有力工具.近年來,多元/多變量/多通道信號分解理論方興未艾,在諸多領域得到了成功應用,但目前尚未見到相關綜述報道.為了填補這個空缺,本文從單變量和多變量兩個方面對主流信號分解方法及其主要改進進行綜述,分析比較這些方法的原理和優(yōu)缺點,并在最后進行了研究展望.

        1 傳統(tǒng)信號分解與時頻表示

        信號處理的一個重要任務就是分析信號中的頻譜成分,同時還希望知道不同頻率成分所出現的時間.傳統(tǒng)的信號處理方法主要集中于信號變換,例如通過傅里葉變換(Fourier transform,FT)及其反變換,人們構建起了信號時域與頻域之間轉換的橋梁.以傅里葉變換為基礎的分析方法已經統(tǒng)治了線性時不變與平穩(wěn)信號處理領域近200 年,尤其是Cooley 等[5]利用傅里葉算子的周期性和對稱性,提出了快速傅里葉變換算法,將N點傅里葉變換的乘法計算量從N2次降為 (N/2×log2N) 次,這成為數字信號處理發(fā)展史上的轉折點和里程碑.以此為契機,伴隨著超大規(guī)模集成電路和計算機領域的迅猛發(fā)展,傅里葉變換不但已成為一個重要的數學分支,而且也成為信號分析和處理的重要工具,并在眾多領域得到了廣泛應用.但是,在科學研究與工程應用中,研究人員發(fā)現,傅里葉變換具有較多不足,主要體現在三個方面[1]: 1)傅里葉變換缺乏時間和頻率的定位功能.傅里葉變換得到的結果是信號在整個積分區(qū)間的時間范圍內所具有的頻率特征的平均表示,無法通過傅里葉變換知道在某一個特定時刻或較短時間范圍內的信號頻率信息.2)傅里葉變換對于非平穩(wěn)信號的局限性.只有時不變的信號才能夠展開為無窮多個復正弦函數的和,而且這無窮多個復正弦信號的幅度、頻率和相位都不隨時間變化,即取某一特定常數.因此,傅里葉變換只適合于處理平穩(wěn)的時不變信號,但從實際過程中采集到的數據往往具有時變特性[1],屬于非平穩(wěn)信號[5-6].3)傅里葉變換在分辨率上的局限性.傅里葉變換受到不定原理的制約[1],無法根據信號的特點來自適應調整時域分辨率和頻域分辨率.但是,正是傅里葉變換的這些不足成為了幾十年來推動人們尋找新的信號分析與處理方法的動力.

        Gabor[7]提出短時傅里葉變換(Short-time Fourier transform,STFT)進行時域和頻域的聯(lián)合分析.這種方法簡單易懂,但是不能自動調節(jié)時域窗口和頻域窗口.Cohen[8]給出了各種時頻分布的統(tǒng)一表示形式,稱為Cohen 類時頻分布,其中最核心的就是Wigner-Ville 分布(Wigner-Ville distribution,WVD).Wigner 分布是由Wigner[9]在1932年提出的概念,但直到1948 年,Ville 才將其應用于信號分析,因此命名為Wigner-Ville 分布.實際上,WVD 可理解為在一個特定區(qū)間上的傅里葉變換,因此它仍然受到不定原理的制約.WVD 的時頻分析性能也易受交叉項的影響.此外,WVD 的瞬時頻率是時間的單值函數,在處理多分量信號時,只能給出多個頻率在一個時間點上的均值,這樣就無法刻畫多分量信號的頻率分布.小波變換(Wavelet transform,WT)是過去20 多年信號處理領域最重要的進展之一[5],它最大的特點是在基本小波中引入了尺度因子,使得小波具有自動調節(jié)時域和頻域分辨率的能力.然而,在對信號進行小波變換前如何選擇小波基函數,仍是一個未解決的現實問題,一般需要使用者不斷試用各類小波.

        2 一元信號分解

        上述傳統(tǒng)信號處理方法都是基于基函數展開的思路,具有簡單、唯一和對稱等優(yōu)點[10],但它們的基函數都是預先定義好且固定不變的,缺乏靈活性,而且受到不定原理的制約,其時頻分析結果也比較模糊[11].在許多涉及非平穩(wěn)信號的實際應用中,這些缺陷都是亟待解決的.因為基于數據驅動的信號分解與時頻分析方法對輸入數據很少有或幾乎沒有先驗假設,所以人們對這些方法產生了極大的興趣.

        這一趨勢始于90 年代末,當時Huang 等[3]提出一種遞歸算法,稱為經驗模態(tài)分解(Empirical mode decomposition,EMD).EMD 通過利用信號極值的遞歸篩選過程,將輸入信號分解為固有的振蕩模式,稱為本征模態(tài)函數(Intrinsic mode function,IMF).一個振蕩信號能夠被稱為IMF 需滿足兩個條件: 1)在其時間區(qū)間內,模態(tài)的極值點數目和過零點的數目應當相等或最多相差一個;2)在其時間區(qū)間內,分別由信號的局部極大值和極小值確定的上包絡和下包絡的均值為零.其中,第1 個條件保證了IMF 是一個窄帶信號;第2 個條件則從信號局部特征的角度考慮,避免了由于信號波形不對稱而引起的瞬時頻率波動.圖1 是一個典型的IMF 波形示意圖,其極值點數目和過零點數目相同,且上包絡和下包絡的均值為零[1,5].

        圖1 一個IMF 的波形示意圖Fig.1 Waveform diagram of an IMF

        由上述討論和圖1 可以看出,IMF 類似于傅里葉變換得到的一個具有固定幅度和固定頻率的分量,但是IMF 包含了幅度調制和頻率調制的特性,因此更具一般性.

        EMD 的篩分過程如圖2 所示[1].首先,對信號極值進行插值,并對上包絡和下包絡進行平均,從而獲得信號的局部均值,這些局部均值可視為信號中的低頻成分估計;然后,將低頻成分從輸入信號中迭代分離出來,得到高頻(快速振蕩)成分.這樣就完成了一次篩分[12].重復篩分過程,直到輸入信號中的所有主要振蕩模態(tài)都被提取出來.由于EMD完全是由數據驅動的,避免了線性時頻變換方法的很多局限性,如受Heisenberg 不定原理限制而導致的有限時頻分辨率和由使用固定基函數而獲得的無意義的模態(tài)等.因此,EMD 自1998 年問世,便在科學界產生了重大影響,揭開了后小波時代時頻分析領域的新篇章,并被應用到了眾多工程領域,如語音增強[13]、圖像處理[14]、設備診斷學[15]、生物醫(yī)學[16]、氣候學[17]、地球物理學[18]等.

        圖2 EMD 的篩分過程示意圖Fig.2 Schematic diagram of sifting process of EMD

        雖然EMD 在處理非線性非平穩(wěn)信號方面具有較多優(yōu)點,但是仍然有一些待解決的問題[18-20],主要包括迭代停止準則與樣條函數的選擇問題、端點效應與模態(tài)混疊問題、Hilbert 變換與分量正交性問題等.EMD 的一些改進版本,例如集成EMD (Ensemble empirical mode decomposition,EEMD)[21]、互補EMD (Complementary ensemble empirical mode decomposition,MCEEMD)[22]、中值EEMD(Median ensemble empirical mode decomposition,MEEMD)[23]和中值互補EEMD[24]等,在一定程度上修復了這些問題,但 EMD (包括分解算法、IMF的定義以及上述的EMD 算法)的各種改進版本都是建立在經驗基礎之上的,目前仍然缺乏堅實的數學理論分析做支撐.對此,Huang 等[18,20]指出,目前EMD 的理論水平猶如小波變換在20 世紀80 年代初的水平,他特別期望能有如Daubechies I.那樣,為小波變換奠定堅實理論基礎的學者出現,能將基于EMD 的時頻分析方法置于堅實的數學理論之上[5].為了克服EMD 理論分析上的困難,部分學者嘗試用基于魯棒約束優(yōu)化的方法取代EMD 中包絡和局部均值估計[25-26],從理論上保證了算法能收斂到全局最優(yōu),但這些方法會遇到沒有可行解的問題[27].Lin 等[28]提出迭代濾波分解(Iterative filtering decomposition,IFD),試圖彌補EMD 數學理論上的不足.IFD 使用特定的移動均值替代EMD篩分過程中的包絡均值,實現了自適應濾波,提升了分解穩(wěn)定性,并嚴格證明了在一定條件下的算法收斂性[29].Cicone 等[30]進一步為IFD 設計了濾波器長度自適應更新策略,使得濾波器滿足迭代濾波收斂的充分條件,為非線性和非平穩(wěn)信號處理提供了一個完整的局部分析工具箱.

        還有一些其他相關方法,例如匹配追蹤[31]這種經典的原子分解算法.與常見的正交基函數相比,原子字典具有過完備性,因此可以更加靈活地表征復雜信號.匹配追蹤算法將原子庫中與當前信號最相關的原子作為當前最優(yōu)原子,經過多次迭代,可將信號表示成多個最優(yōu)原子的線性疊加模型.匹配追蹤算法的時頻分析效果與原子字典復雜程度密切相關,當分析復雜信號時,往往需要復雜原子模型,使得分解結果不稀疏,降低了算法的效率,也會導致時頻特征間斷問題[32].文獻[33-36]受EMD 和壓縮感知理論的啟發(fā),在本征模態(tài)函數組成的字典中,尋找多尺度數據的稀疏表示,將信號分解問題轉化為非線性L1優(yōu)化問題,并提出一種迭代算法遞歸求解該非線性稀疏優(yōu)化問題,實現了稀疏時頻分析.這些稀疏時頻分析方法的分解結果與EMD 較類似,但其性能不受EMD 中停止準則的影響,抗噪聲和抑制端點效應的能力也優(yōu)于EMD.此外,Peng等[37]和Guo 等[38]提出基于算子的信號分解方法,該方法采用零空間追蹤的方式,自適應估計算子和參數,將信號中的成分分離.這些方法的性能取決于所構造的算子,常見的算子有奇異局部線性算子[37]和復數微分算子[38]等.

        除了EMD 的改進版本以外,學術界還提出一些類似于EMD 迭代篩分過程的新分解算法,來嘗試解決這些問題.例如 Smith[39]通過分離調頻和調幅成分,提出局部均值分解(Local mean decomposition,LMD)算法.Frei 等[40]通過分段計算旋轉分量在每一個局部時間區(qū)間上的瞬時頻率和瞬時幅值,提出計算復雜度較低的本征時間尺度分解(Intrnsic time-scale decomposition,ITD)算法,可以實現數據的在線分解與實時處理.ITD 的分解過程[1]如圖3 所示,其中Xk表示輸入信號x(t) 的極值點;虛線L(t) 和虛線H(t) 分別表示基線和旋轉分量.需要指出的是,上述分解算法仍然都是基于經驗的分解算法,分解結果也會受到端點效應和模態(tài)混疊的影響.另外,這些方法都需要準確估計信號的局部極值或均值,但實際數據往往會被噪聲干擾,不易準確估計這些極值或均值,這使得這些方法的抗噪聲能力較差.

        圖3 ITD 分解過程示意圖Fig.3 Schematic diagram of ITD decomposition process

        另一類類似于EMD 的數據驅動信號處理方法包括同步壓縮變換(Synchro squeezed transform,SST)[41]和經驗小波變換(Empirical wavelet transform,EWT)[12].與EMD 不同,SST 是通過一套合理且方便的數學框架,來提取信號中的模態(tài),它既可以在STFT 域,也可以在小波域中工作[42].SST首先通過頻率重整算子[43]銳化信號的STFT 譜圖或小波譜圖;然后,在模態(tài)總數已知條件下,采用脊線提取技術,來估計瞬時頻率;最后,通過在相應的脊線附近對重整STFT 或重整小波變換積分,實現模態(tài)重構.需要注意的是,SST 要求模態(tài)在時頻面上各自獨立,無法處理時頻信息交叉的情況.

        EWT 是加州大學洛杉磯分校的Gilles[12]提出的,這種方法融合了經驗模態(tài)分解和小波變換的優(yōu)點,通過峰值檢測機制對頻譜進行分割,并在每個分割區(qū)間構建小波濾波器,實現將一個復雜信號分解為一系列具有緊支撐頻譜的調頻調幅信號.圖4是EWT 的模態(tài)頻譜分割示意圖[32],其中豎虛線表示檢測到的模態(tài)邊界頻率.EWT 的效果與基于頻譜分割所構造自適應的小波濾波器是否精確有很大關系,但對實際信號進行準確的頻譜分割是一項頗有挑戰(zhàn)性的任務[44].

        圖4 EWT 的模態(tài)頻譜分割示意圖Fig.4 Schematic diagram of modal spectrum division of EWT

        同樣來自加州大學洛杉磯分校的Dragomiretskiy等[45]在2013 年提出變分模態(tài)分解算法(Variational mode decomposition,VMD).VMD 的目標是將輸入信號分解為一系列具有稀疏特性的模態(tài),這里的稀疏特性指的是所有模態(tài)都是集中在各自中心頻率附近的窄帶信號.為了實現這個目標,VMD 構建了一個帶約束的變分優(yōu)化問題,其中目標函數是最小化所有模態(tài)的帶寬,約束條件是分解得到的模態(tài)能夠完全重構輸入信號.VMD 目標函數的構造分為3 個步驟: 1)對每個模態(tài)進行Hilbert 變換,得到其解析信號;2)將解析信號的頻譜平移到零中頻,得到基帶信號;3)利用H1高斯平滑度,估算每個模態(tài)的基帶信號帶寬,并將這些帶寬的和最小化作為目標函數.VMD 所建立的優(yōu)化問題可以通過交替方向乘子法在頻域求解,最終得到輸入信號中包含的模態(tài)和對應的中心頻率.

        與之前基于經驗的信號分解方法不同,VMD的目標函數具有較完備的數學理論基礎,它的求解公式表明,VMD 在本質上是一個自適應最優(yōu)Wiener濾波器組,其中心頻率如圖5 中豎虛線所示[32].由于其優(yōu)良的特性,VMD 一經提出,就受到了極大的關注,目前已成功在機械故障診斷[46]、風速預測[47]、疾病診斷[48]、金融數據分析[49]、新型冠狀病毒肺炎預測[50]等領域得到應用.

        圖5 VMD 原理示意圖Fig.5 Schematic diagram of VMD principle

        VMD 的出現是信號分解領域研究的一個分水嶺.在此之前,以EMD 為代表的基于經驗的分解方法雖然得到了廣泛的關注與應用,但在數學基礎理論方面的進展艱難而緩慢,以至于EMD 最開始被提出的時候頗有爭議,甚至受到了質疑[51].VMD的出現改變了信號分解方法研究全憑經驗的現狀,提供了具有堅實數學優(yōu)化理論的研究視角.但是VMD 也有一些不足,例如VMD 的分解性能與模態(tài)數量K和懲罰系數α這兩個參數密切相關,并且這兩個參數無法自適應獲得,需要用戶事先指定.

        目前,對于VMD 的主要改進工作也集中在如何確定這兩個參數的問題上.這些改進VMD 參數依賴性的方法可分為兩類: 1)固定懲罰系數α,通過窮舉法迭代優(yōu)化模態(tài)數量K.例如Li 等[52]使用近似完全重構的判據,來確定合適的模態(tài)數量;Lian等[53]通過判斷提取出來的模態(tài)特性,來選取合適的模態(tài)數量;Cai 等[54]利用傅里葉變換得到的頻譜,作為判斷模態(tài)數量是否合適的標準.上述這類改進方法雖然簡單,但由于VMD 的性能受到模態(tài)數量K和懲罰系數α的共同調節(jié)[55],只考慮模態(tài)數量的作用存在一定的風險.2)同時對模態(tài)數量K和懲罰系數α進行優(yōu)化.這類方法的思路是,首先,針對特定的應用場景,構造適應度函數;然后,使用一些啟發(fā)式智能算法,對參數進行尋優(yōu).如人工魚群算法[56]和蝗蟲算法[57]等.盡管這類基于智能尋優(yōu)的改進方法考慮了模態(tài)數量K和懲罰系數α這兩個參數,但這些改進方法沒有觸及VMD 的本質,且有效性受限于特定場景下的適應度函數和優(yōu)化算法,難以推廣到其他場景.為了解決VMD 的參數選取問題,Chen等[58]提出一種自整定變分模態(tài)分解算法(Self-tuning variational mode decomposition,SVMD).SVMD 將VMD 的目標函數由同時提取K個模態(tài)改進為每次只提取一個模態(tài),且在提取模態(tài)的過程中,懲罰系數α不再采用固定值,而是會根據信號特點自適應地對α進行更新.

        VMD 實際上只是對窄帶信號的分解效果較好,這意味著基于VMD 的方法處理帶寬較大時變信號的能力非常有限.為了解決這個問題,2017 年,Chen等[59]提出非線性調頻模態(tài)分解算法(Nonlinear chirp mode decomposition,NCMD).NCMD 通過解調技術[60],將帶寬較寬的時變信號進行頻率解調,可以使之轉換為窄帶信號.NCMD 頻率解調過程如圖6 所示[59],圖中曲線、上方橫線和下方橫線分別表示原信號、解調信號和基帶信號的時頻曲線,BW和 B Wmin分別表示原始信號和解調信號的帶寬.利用這一特性,NCMD 將信號分解問題轉化為解調問題,并建立了帶約束的非線性調頻模態(tài)變分優(yōu)化模型.

        圖6 NCMD 頻率解調過程示意圖Fig.6 Schematic diagram of NCMD frequency demodulation process

        與VMD 類似,NCMD 亦采用了交替方向乘子法求解所提出的變分優(yōu)化模型.求解所得更新公式表明,NCMD 可以被視為一個時頻濾波器組,該濾波器組的中心頻率就是估計得到的瞬時頻率.NCMD 的計算過程可以概括為: 首先,用當前估計得到的瞬時頻率對解調的正交信號進行迭代更新;然后,通過反正切解調技術[61],利用正交信號的相位信息進一步更新瞬時頻率;重復上述兩個步驟,直到解調后的信號具有最窄的頻帶.與VMD 相比,NCMD 的最大優(yōu)勢在于能夠處理非平穩(wěn)信號,并提供準確的時頻信息.雖然NCMD 提出的時間很短,但是由于其在非平穩(wěn)信號時頻分析上的出色性能,目前已被應用于碰撞摩擦故障檢測[62]、生理信號監(jiān)測[63]、高速鐵路故障診斷[64]等領域.值得指出的是,雖然NCMD 相比于VMD 有較大優(yōu)勢,但NCMD也具有和VMD 類似的不足,即分解性能依賴先驗知識,要求使用者預先指定分解模態(tài)數量和懲罰系數.

        為了解決這個問題,Chen 等[65]提出NCMD 的改進版本,稱為自適應調頻模態(tài)分解算法(Adaptive chirp mode decomposition,ACMD),并將其成功應用到旋轉機械的碰撞摩擦故障檢測中[62].ACMD采用一種類似于匹配追蹤[31,34]的貪婪算法,對信號模態(tài)逐一進行遞歸估計,直到提取出信號中的所有的模態(tài).雖然ACMD 在自適應時頻分析上顯示出重大優(yōu)勢,但Chen 等[66-68]發(fā)現,ACMD 算法的停止條件對旋轉機械故障信號比較有效,但在分解工業(yè)過程振蕩信號時,頻繁出現難以收斂導致過分解問題,因此改進了ACMD 算法的收斂判據,提出快速自適應調頻模態(tài)分解算法[67],緩解了ACMD 過分解問題,并減少了分解用時.

        與NCMD 同年,Chen 等[69-70]提出本征調頻模態(tài)分解算法(Intrinsic chirp mode decomposition,ICMD).考慮到很多應用場景采集的信號模態(tài)的瞬時頻率和瞬時幅值具有連續(xù)光滑的特性,ICMD 采用傅里葉級數或多項式函數逼近模態(tài)的瞬時頻率和瞬時幅值,從而構造了一種通用的非線性調頻模態(tài)參數化模型,能夠有效刻畫不同調制程度的信號分量.在此模型基礎上,ICMD 利用廣義參數化時頻變換方法[71]估計模態(tài)的瞬時頻率,進而通過正則化最小二乘法估計模態(tài)的瞬時幅值,最終實現信號模態(tài)的重構.ICMD 在非平穩(wěn)信號的分解效果上,與NCMD 難分伯仲,但由于ICMD 并未采用NCMD所使用的交替方向乘子法這種迭代式的優(yōu)化方法,而是采用簡單的最小二乘法,即可得到結果,因此在算法復雜度上,ICMD 具有顯著優(yōu)勢.目前,ICMD已被成功應用于信號消噪[72]、雷達信號處理[73]、機械故障診斷[74]等領域.

        本節(jié)對眾多單變量信號分解方法進行了介紹和分析,選擇EMD、VMD 和NCMD 作為代表方法進行案例驗證,其中EMD 是最經典的信號分解方法,已得到廣泛關注和應用,因此選之作為經驗性信號分解方法的代表;VMD 是信號分解研究由經驗性轉向數學優(yōu)化建模的標志性方法,因此選之作為基于優(yōu)化的信號分解方法的代表;NCMD 是最近提出的能處理時變信號的分解方法,推動了基于優(yōu)化的信號分解方法向縱深發(fā)展,因此選之作為近年提出的有較大影響力的信號分解方法的代表.本文以單變量信號(1)為例,它包含3 個模態(tài),其中2 個模態(tài)是固定頻率的正弦函數,頻率分別是10 Hz 和30 Hz,最后一個模態(tài)具有時變瞬時頻率(50+20t)Hz,信號中加入了一定程度的噪聲:

        EMD、VMD 和NCMD 的分解結果如圖7、圖8和圖9 所示,圖中x為原始信號,IMF 為分解所得模態(tài).可以看出,EMD 的第1 個模態(tài)提取得較好,但第2 個模態(tài)和第3 個模態(tài)之間出現了明顯的模態(tài)混疊效應,這是EMD 的典型問題;VMD 由于其原理所限,難以處理帶寬較大的時變信號,因此其分解結果受時變模態(tài)的影響,出現了很大誤差;而NCMD 由于其算法中采用了解調算子處理時變瞬時頻率,能夠做到不受時變瞬時頻率的影響,正確提取出了3 個模態(tài).這些實驗結果與本文中對各個方法的分析是一致的.

        圖7 EMD 分解結果Fig.7 The decomposition results of EMD

        圖8 VMD 分解結果Fig.8 The decomposition results of VMD

        圖9 NCMD 分解結果Fig.9 The decomposition results of NCMD

        綜上所述,歷經20 余年的發(fā)展,單變量信號分解研究取得了豐碩成果,已成為現代信號處理領域的研究熱點.本文將常見方法根據其作用域分為時域、頻域和時頻域等多個類別,并在表1 總結了它們的優(yōu)點和局限性.

        表1 常見單變量信號分解方法歸類總結Table 1 Classification and summary of common univariate signal decomposition methods

        3 多元信號分解

        3.1 多元信號分解方法

        雖然基于數據驅動的非平穩(wěn)信號分解和時頻分析工作在蓬勃發(fā)展,但拓展現有的數據驅動方法來處理非平穩(wěn)多變量/多通道/多元信號,也引起了學術界和工程界的極大興趣[75-76].由于傳感器和計算機技術的進步,在現代科學和工程應用中,廣泛存在對多變量/多通道/多元類型的數據處理方法的需求[77],例如基于多通道腦電[78]或心電信號[79]的分類、多變量信號消噪[80-81]、圖像融合[14]等.針對非平穩(wěn)數據的多變量信號分解與時頻分析技術的主要挑戰(zhàn)有以下兩點: 1)具有模態(tài)齊整特性.即不同通道中具有相同或相似的頻率的模態(tài)出現在同一尺度[82].2)提取出多通道之間的相關信息[83].在信號分解領域的研究中,多元、多通道和多變量表示相同的意思.圖10 給出了多元/多通道/多變量信號分解領域的一些術語解釋,以便理解本文后續(xù)內容.

        圖10 多元/多通道/多變量信號分解領域術語的圖形化解釋Fig.10 Graphical interpretation of terms in multivariate signal decomposition

        最直接的多變量信號處理方法是使用單變量的方法,逐一、單獨分析多變量信號的每一個通道.但文獻[75,83-84]指出,對于多變量信號,如果采用單變量信號分解技術逐一處理每個變量,會導致信息泄露和變量之間相關信息缺失等問題,且不滿足模態(tài)齊整[75]的要求.因此,需要為單變量信號分解方法應用至多變量情形開發(fā)特定的拓展方法,以便直接在多變量信號所在的多維空間里處理多變量信號.這也是研究者們開展多元信號分解研究的目的.參照文獻[85]的案例,本文以單變量ICMD 和多變量ICMD 處理多變量信號(2)為例:

        單變量和多變量信號分解方法處理多變量信號(2)結果見圖11 和圖12.可以看出,用單變量方法分解多變量信號時,不僅無法滿足模態(tài)齊整要求,而且分解誤差也較大;反之,多變量信號分解方法的性能令人滿意,分解所得模態(tài)的誤差較小,不同通道中具有相似頻率的成分也位于同一個尺度,滿足模態(tài)齊整性質要求.因此,研究多變量信號分解方法處理多變量信號是十分必要和有意義的.

        圖11 單變量ICMD 分解多變量信號的結果Fig.11 The decomposition results of multivariate signals by the univariate ICMD

        圖12 多變量ICMD 分解多變量信號的結果Fig.12 The decomposition results of multivariate signals by the multivariate ICMD

        多變量信號分解工作是從EMD 應用于復數信號處理開始的.2007 年,Tanaka 等[86]基于復數域的性質,巧妙地利用了原始的EMD,來分解二元時間序列,并把這種方法命名為CEMD (Complex empirical mode decomposition),揭開了數據驅動信號分解與時頻分析方法向多變量拓展的序幕.但Tanaka 等[86]的CEMD 僅局限于處理復數信號,不適用于三變量以上情況.此外,CEMD 無法保證復數信號的實部和虛部分解所得的IMF 數量一致,因此這種思路的適用性非常有限.隨后,Altaf 等[87]將EMD 中的單變量的“振蕩”與雙變量的“旋轉”概念相對應,認為雙變量信號是快速旋轉和慢速旋轉的疊加,并給出了復數空間中極值的定義,在此基礎上,提出旋轉復數EMD (Rotation complex empirical mode decomposition,RCEMD).

        遺憾的是,Altaf 等[87]的工作只將這個思路應用到了復數信號上,并未充分挖掘其潛能.很快,Rilling等[84]受RCEMD 的啟發(fā),將雙變量輸入信號通過投影映射到單位圓上,并在此基礎上,清晰明確地給出了在雙變量空間中信號極值、均值和包絡的定義,正式提出了雙變量EMD (Bivariate empirical mode decomposition,BEMD).BEMD 的分解原理見圖13,兩種二維包絡的均值計算示意圖見圖14[84],圖中曲線上藍點表示極值點,中心紅點表示均值點.圖15 進一步對二維局部極值點進行了說明,根據微積分中局部極值點的相關定義,當選擇Y軸方向為投影方向時,A點將被視為局部極大值點;同樣,B點可視為相應的局部極小值點.因此,以投影方向為基準,可將多變量數據的極值點與投影方向的極值點一一對應,構成了后來基于高維空間投影的多變量信號分解技術的基礎.

        圖15 二維局部極值點示例Fig.15 Example of two-dimensional local extreme points

        Ur Rehman 等[88]意識到Rilling 等[84]的雙變量拓展思路可以被推廣到更多變量的情形.Ur Rehman等[88]通過在高維空間中,建立均勻分布的投影向量集,分別計算輸入信號在各個方向上的投影包絡線,然后通過計算包絡線的均值,定義多變量信號的局部均值函數,并在此基礎上,提出三變量EMD (Trivariate empirical mode decomposition,TEMD)[88]、四變量EMD (Quadrivariate empirical mode decomposition,QEMD)[89]和多變量EMD (Multivariate empirical mode decomposition,MEMD)算法[75].Ur Rehman 等[75]提出的MEMD 算法正式開啟了多變量信號分解與時頻分析的新時代,迅速受到了來自各個領域研究人員的關注.目前,已被應用于腦機接口[90]、設備診斷[91]、因果分析[92]、地球物理[93]、生物醫(yī)學[78]等諸多領域.然而,MEMD 也繼承了原始EMD 的所有局限性,例如對采樣頻率敏感、噪聲魯棒性差以及EMD 算法的經驗性.此外,MEMD隨著輸入信號數目的增加,計算量會出現指數型增長.為了提高MEMD 的計算效率,Lang 等[94]提出快速MEMD 算法(Fast multivariate empirical mode decomposition,FMEMD),重新定義了多變量IMF 的基本概念,并在此基礎上,架構了FMEMD與EMD 方法之間的映射關系,使得FMEMD 計算量與輸入信號的數量無關.FMEMD 中一個雙變量信號及其投影信號、局部均值和多變量IMF 的示意圖見圖16、圖17、圖18 和圖19.

        圖16 雙變量信號Fig.16 Bivariate signal

        圖17 雙變量信號的投影信號Fig.17 Projection signal of bivariate signal

        圖19 多變量IMFFig.19 Multivariate IMF

        Lang 等[95-96]創(chuàng)造性地提出了兩種多變量形式的ITD 算法,分別稱為間接多變量ITD (Indirect multivariate intrinsic time-scale decomposition,IMITD)和直接多變量ITD (Direct multi-variate intrinsic time-scale decomposition,DMITD).IMITD 與FMEMD 的思路相似,即通過定義多變量與單變量運算空間映射關系,架構超定線性方程組,以求解多變量分解結果.需要指出的是,IMITD 使用Halton-Hammersley 采樣技術獲得高維空間的采樣點.與使用等角度采樣技術獲得的采樣點相比,Halton-Hammersley 采樣可以使得采樣點對高維空間的覆蓋更加均勻.

        等角度采樣和Halton-Hammersley 序列采樣在三維球體上的采樣點見圖20 和圖21.可以看出,等角度采樣結果在球體兩極更為集中,不利于后續(xù)的均勻投影操作.IMITD 在局部特征處理上,比MEMD 的效果要好,計算效率也較低,但如果投影方向選擇不恰當,會導致IMITD 不能正確提取出期望的基線,從而引起IMITD 分解產物沒有物理意義.

        圖20 等角度采樣Fig.20 Uniform angle sampling

        圖21 Halton-Hammersley 序列采樣Fig.21 Halton-Hammersley sequences based sampling

        DMITD 是通過恰當定義多變量極值點、多變量基線節(jié)點和多變量基線算子等概念,直接拓展ITD 算法至多變量情形,實現ITD 算法在多變量空間的運算操作.與IMITD 相比,DMITD 對投影方向的魯棒性要好一些,但運算效率低.上述多變量信號分解方法都繼承了與之相對應的單變量信號分解的缺陷,例如模態(tài)混疊與端點效應.此外,由于這些方法的思路是通過將多變量輸入信號投影映射到高維的空間中再分解,因此效果均與投影向量的數目和方向有關.目前還沒有明確的依據來確定投影方案,相關研究尚待完善.

        多變量SST (Multivariate synchrosqueezed transform,MSST)[97]和多變量EWT (Multivariate empirical wavelet transform,MEWT)[98]放棄了將輸入多變量信號投影映射到高維空間的做法,開辟了基于小波變換的多變量信號處理新思路.具體地,MSST 首先對每個信號通道分別應用標準SST 算法;然后,對時頻域進行自適應劃分,以分離輸入數據中的單模態(tài)多變量振蕩;最后,估計出多變量瞬時頻率和幅值,并在此基礎上,計算出一個多變量同步壓縮變換算子.MSST 可以得到多變量信號清晰的時頻譜,對于探索性的數據分析非常有用.但MSST 只給出了時頻譜表達,不能重構模態(tài),因此MSST 的應用范圍是有限的[97].MEWT 首先采用模態(tài)估計過程來獲取多變量數據中的最優(yōu)信號;然后,對其相應的頻譜進行分割,以恢復出所有輸入信號通道中的所有模態(tài).這種方法存在EWT 固有的缺陷,即需要基于有效的頻譜分割來顯式構造自適應小波濾波器組[98],這個缺陷對實際物理系統(tǒng)中的信號,難以做到頻譜的有效分割.

        2019 年,受單變量調制振蕩信號推廣至多變量形式的啟發(fā)[99-100],Ur Rehman 等[83]提出多變量VMD算法(Multivariate variational mode decomposition,MVMD).雖然2017 年Wang 等[101]就提出了復數VMD (Complex variational mode decomposition,CVMD),但CVMD 遭遇了與CEMD 一樣的問題,即僅局限于處理復數信號,不適用于三變量以上的情況,且無法保證復數信號的實部和虛部分解所得的IMF 數量一致,因此這種思路的適用性非常有限.MVMD 從輸入多變量信號中,尋找一組共同的多變量調制振蕩,這些振蕩在完全重構輸入數據所有通道的同時,具有最小的帶寬和.MVMD是極具潛力的多變量信號分解方法,這種拓展形式不僅繼承了標準VMD 的許多理想性質,還展現出優(yōu)良的模態(tài)齊整特性.MVMD 現已在風機故障診斷[102]、腦電信號檢測[103]、基因工程[104]和信號消噪[105]等領域得到成功應用.但 MVMD 直接將模態(tài)頻譜的重心的估計作為中心頻率,因此不適宜處理帶寬較大的時變信號.此外,MVMD 也不能直觀地提供時頻信息.

        目前,大多數多變量信號處理方法局限于處理窄帶信號,對時變信號的分解能力有限,因此Chen 等[77]受MVMD 的啟發(fā),提出多變量非線性調頻模態(tài)分解算法(Multivariate nonlinear chirp mode decomposition,MNCMD),解決了連續(xù)時變條件下多變量信號分解與時頻分析問題.MNCMD 現已被成功應用于過程控制系統(tǒng)中的復雜多重廠級振蕩根因分析[106-107].MNCMD 雖然在時變多元信號分解與時頻表示任務上性能突出,但是它的計算復雜度較高,達到隨后,Chen等[85,108]進一步結合傅里葉級數建模,提出多變量本征調頻模態(tài)分解算法(Multivariate intrinsic chirp mode decomposition,MICMD),能夠以 O (N) 復雜度達到與MNCMD相似的時變多元信號分解和時頻分析效果.

        本節(jié)選用MEMD、MVMD 和MNCMD 作為代表性方法,進行案例驗證實驗.其中,MEMD 使用了高維空間投影的方式進行多變量拓展,也是第1 個多變量信號分解方法,因此選之作為以高維空間投影進行多變量拓展的代表方法;MVMD 使用了多元調制振蕩的方式進行多變量拓展,因此選之作為這種多變量拓展方式的代表;MNCMD 是最近提出的能處理時變多元信號的分解方法,因此選之以展示對時變多元信號的處理能力.

        本文以多變量信號(3)為例,該信號由四個通道組成:x1通道包含一個時變模態(tài)和兩個時不變模態(tài),x2通道只包含一個時變模態(tài),x3通道包含兩個時不變模態(tài),x4通道包含一個時變模態(tài)和一個時不變模態(tài).信號中加入了噪聲:

        MEMD、MVMD 和MNCMD 的分解結果見圖22、圖23 和圖24,圖中第1 行為原始信號,IMF1~IMF5為分解所得模態(tài).可以看出,在MEMD的分解結果中,產生了很多冗余模態(tài),雖然在該結果中,時不變模態(tài)提取的效果不錯,但時變模態(tài)出現了嚴重的分裂現象;MVMD 雖然避免了MEMD產生冗余模態(tài)過多問題,但由于其原理所限,分解結果中的時變模態(tài)誤差很大,尤其是在時變模態(tài)的高頻部分,有一部分直接泄露到了第2 個模態(tài)中;而MNCMD 能夠很好地提取出這些時變和時不變模態(tài),并且表現出模態(tài)齊整性.這些實驗結果與本文中對各方法優(yōu)缺點的分析是一致的.

        圖22 MEMD 的分解結果Fig.22 The decomposition results of MEMD

        圖23 MVMD 的分解結果Fig.23 The decomposition results of MVMD

        圖24 MNCMD 的分解結果Fig.24 The decomposition results of MNCMD

        綜上所述,多元信號分解歷經10 余年的發(fā)展,涌現了諸多多元信號分解拓展方式,例如高維空間投影和多元調制振蕩等,提出了很多各具特點的多元信號分解方法.隨著多元信號分解研究的興起與蓬勃發(fā)展,研究者們在此過程中發(fā)現了新的問題,給出了新的定義,發(fā)展了新的概念,揭示了新的性質,極大豐富了信號分解的研究和應用范圍.表2總結了主要多元信號分解方法的優(yōu)點和局限性.

        表2 多元信號分解方法歸類總結Table 2 Classification and summary of multivariate signal decomposition methods

        3.2 多元信號分解應用

        多元信號分解方法已經在很多領域得到了廣泛應用.本文以生物醫(yī)學工程、工業(yè)控制系統(tǒng)、機械故障監(jiān)測、時間序列預測等領域的研究為例,進行說明.

        1)在生物醫(yī)學工程領域,多元信號分解方法的典型應用是心電信號ECG (Electroen cephalo gram)和腦電信號EEG 的處理.ECG 信號是反映心臟基本功能和病理信息的重要參考值.ECG 信號的采集需要在胸部的右上、左上、左下和右下側連接多個導聯(lián),以觀察心電圖的變化,這就構成了多個信號通道,形成了多變量信號.例如文獻[109-110]使用MEMD 處理ECG 信號,移除了基線漂移的影響,提高了心血管疾病診斷的準確性.在進行腦電數據采集時,通常會根據需要選取不同數量的導聯(lián)或電極點,在不同的腦區(qū)采集EEG 信號,這就形成了多變量EEG 信號.在腦機接口研究中,各種復雜腦電信號的非穩(wěn)定性和多通道性一直阻礙著常規(guī)基于基函數信號處理方法的應用及多通道同步分析.文獻[78,111]利用MEMD 提取多通道腦電信號中的特征,提升了癲癇疾病診斷效果.文獻[112]利用MVMD 分析多通道腦電信號的時頻特征,并在人類情感識別任務中取得了優(yōu)異的性能.

        2)在工業(yè)控制系統(tǒng)中,多變量信號分解方法在性能評估上也得到了廣泛應用.由于控制系統(tǒng)的大規(guī)模、高集成、強耦合的特點[113-114],不同設備或單元之間會相互影響,極易在系統(tǒng)的多個部位表現出相似的故障特征.典型的例子是過程控制系統(tǒng)中的廠級振蕩.控制系統(tǒng)規(guī)模龐大、機理復雜,使得廠級振蕩表現出非線性、非平穩(wěn)、多模態(tài)、強噪聲等特征,這影響了廠級振蕩的檢測和診斷效果,嚴重時甚至威脅整個系統(tǒng)的穩(wěn)定性和安全性[115].文獻[95,116-117]分別使用了MEMD、MITD 和MNCMD,對不同的工業(yè)控制系統(tǒng)廠級振蕩數據進行分析,一致認為,多元信號分解方法是目前廠級振蕩檢測最主流和最有效的手段之一,有利于提升控制系統(tǒng)性能水平.

        3)近年來,多變量信號分解方法在機械故障監(jiān)測領域得到了較多關注.例如文獻[118]將MEMD與Teager 能量譜結合,獲取信號的故障特征信息,實現微小故障特征的提取,并在核主泵軸承外圈早期故障檢測實際任務中的表現優(yōu)于基于EMD 的方法.文獻[119]在改進MVMD 自適應性的基礎上,實現了一種具有突出優(yōu)勢的軸承故障診斷方法.文獻[120]利用MNCMD 能處理時變多元信號的特點,實現了對轉子軸承系統(tǒng)在非平穩(wěn)過程(如啟動和關閉)中,瞬時振動狀態(tài)的實時分析與監(jiān)測.

        4)多變量信號分解可以用于提升多變量時間序列預測性能.因為現實時間序列數據往往具有高度非線性和非平穩(wěn)特征,多變量信號分解方法不僅能有效處理非線性和非平穩(wěn)因素,還能充分保留變量之間的相關信息,以提升多變量時間序列預測性能.文獻[121]提出一種基于MEMD 和支持向量回歸的混合預測模型,它的新穎之處主要在于MEMD的應用,使得多元數據分解能夠有效提取不同相關變量之間的固有信息.在多個數據集上的實驗結果表明,基于MEMD 的混合模型是一種很有前景的電力峰值負荷預測方法.除了電力負荷預測之外,多元信號分解還廣泛應用于各種行業(yè)的時間序列預測,例如經濟領域的股價指數預測[122]和原油價格預測[123]、物理學中的太陽輻射預測[124-125]、生態(tài)環(huán)境中的空氣質量預測[126]等.

        綜上所述,多元信號分解方法并不局限應用于某一類場景,可以與多個學科交叉,為各行各業(yè)賦能,從而產生研究價值和應用價值.目前,多元信號分解方法的應用仍然在蓬勃發(fā)展中,表3 歸納總結了幾種常見多元信號分解方法的適用場景.

        表3 常見多元信號分解方法的適用場景Table 3 Applicable scenarios of common multivariate signal decomposition methods

        4 總結與展望

        非線性和非平穩(wěn)是現實信號的普遍規(guī)律,時頻聯(lián)合分析技術正是應現實科學研究和工程應用需求而產生和發(fā)展的.瞬時頻率是研究非平穩(wěn)和非線性信號最有力的工具[127],但對于實際的非線性和非平穩(wěn)信號,由于Bedrosian 定理[128]的限制,難以直接求解瞬時頻率[129].一個自然的做法是將所研究的信號分解為一個個單分量的信號,每個單分量信號只包含一種振蕩模態(tài),這就是基于信號分解的時頻分析研究最初的想法.本文從單變量和多變量兩個方面,梳理了信號分解領域的發(fā)展歷程與研究現狀,比較分析了主流信號分解方法的優(yōu)缺點,在此基礎上,本文認為未來可以從以下五個方向進一步探索:

        1)間歇信號分解理論與技術.雖然Xie 等[130]結合K均值聚類方法,提升了ITD 處理間歇信號的能力,但對于大多數方法,特別是VMD、NCMD、ICMD 等基于優(yōu)化建模的信號分解方法,均要求瞬時頻率和瞬時幅值都是連續(xù)光滑函數,這意味著間歇信號并不滿足這些信號分解方法的假設條件,因此現有的信號分解方法處理間歇信號的能力極為有限.進一步地,目前多變量間歇信號分解理論與技術是一個開放性問題.對于MEMD、MITD、MVMD、MNCMD 和MICMD 等多元信號分解方法中出現的模態(tài)齊整特性,在多變量間歇信號中該如何定義和體現,仍然是一個尚未討論的問題.因此,研究復雜間歇信號分解理論和開發(fā)復雜間歇信號時頻分析技術,對豐富信號分解與時頻分析領域研究具有重要意義,也將極大拓展信號分解理論方法應用范圍[131].

        2)自適應或無參化信號分解理論與技術.雖然SVMD 等方法在一定程度上為特定的信號分解技術提供了參數整定或自適應更新方案,但是這些方案依然會有一些超參數或閾值需要提前指定.需要注意的是,很多改進方案(如EEMD、CEEMD (Complementary ensemble empirical mode decomposition)等),雖然提升了EMD 抗模態(tài)混疊和抗端點效應的能力,但也引入了新的超參數;另外,還有一些多變量拓展方式也引入了超參數,例如基于高維空間投影思路的多元信號分解方法就會引入投影向量方向和數量這兩個超參數.這些超參數會影響到信號分解方法在處理實際信號時的實用性和便捷性.本文認為至少可以從以下兩個角度來緩解這個問題: a)將已有分解方法與參數尋優(yōu)技術結合,為特定場景下的應用需求提供定制化的參數整定方案[53,132-134];b)研究新的參數自適應調整的或無參化的信號分解方法,這具有相當的挑戰(zhàn)性.

        3)新型多變量拓展技術.多變量/多元/多通道信號分解與時頻分析研究方興未艾.對于現有的多變量信號分解與時頻分析技術,目前主要是基于高維空間投影[75,94-95]和多變量調制振蕩[97,100]兩種多元拓展思路.基于高維空間投影方法會帶來投影向量方向和投影向量數目敏感性的問題,基于多變量調制振蕩方法則有時會出現零頻分量波動的問題.如何解決現有多變量拓展技術的不足和開發(fā)新的多變量拓展方案,是一件富有挑戰(zhàn)而又令人期待的任務.

        4)信號分解方法的理論完備性.對于經典的基于經驗性的EMD、LMD、ITD 等方法,雖然應用領域廣泛、使用效果較好,但是其缺乏嚴格的數學理論支持,亟待能有如Daubechies I.那樣能為小波變換奠定堅實理論基礎的學者出現,能將基于EMD的時頻分析方法置于堅實的數學理論之上[1].對于新興的基于優(yōu)化的VMD、NCMD 等方法,雖然具備優(yōu)化模型的理論支持,但其算法收斂性和結果唯一性的證明仍然是一個懸而未決的問題,這使得在使用這類方法時仍然把握不足.因此,進一步完備信號分解理論基礎,明確算法收斂條件,對該領域的長期化與實用化發(fā)展具有重大意義.

        5)新領域和交叉領域的探索融合.信號分解屬于基礎的信號處理技術,是工程學科的基礎學科,是一門“使能”學科,它既可以探索新領域,發(fā)現新需求,又可以與傳統(tǒng)學科領域結合起來,為其他學科的傳統(tǒng)需求賦能.例如引力波的時頻分析研究[135],機械、工業(yè)控制等領域的故障檢測與診斷研究[136-137];復雜系統(tǒng)因果分析研究[106,117]、腦機接口與生理疾病特征識別[110]、新型冠狀病毒肺炎相關問題分析與預測[50,138]、云計算集群中的復雜時序信號周期檢測與分解應用[139-141].信號分解與Transformer[142]結合提出的新型神經網絡,使單變量時序預測和多變量時序預測誤差分別降低了14.8%和22.6%[143-144].諸如此類的舊問題與新需求廣泛存在,蓬勃發(fā)展的信號分解理論方法為這些問題和需求的解決提供了有效的新途徑,這些新領域和交叉領域的發(fā)展,也啟發(fā)著信號分解技術新的突破方向.

        猜你喜歡
        時頻時變模態(tài)
        基于時變Copula的股票市場相關性分析
        智富時代(2017年4期)2017-04-27 17:08:47
        煙氣輪機復合故障時變退化特征提取
        國內多模態(tài)教學研究回顧與展望
        基于MEP法的在役橋梁時變可靠度研究
        基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
        基于時頻分析的逆合成孔徑雷達成像技術
        對采樣數據序列進行時頻分解法的改進
        由單個模態(tài)構造對稱簡支梁的抗彎剛度
        計算物理(2014年2期)2014-03-11 17:01:39
        雙線性時頻分布交叉項提取及損傷識別應用
        淺析《守望燈塔》中的時頻
        国产一区二区三区日韩在线观看| 国精品人妻无码一区免费视频电影| 国产激情久久久久影院老熟女免费 | 欧美日本精品一区二区三区| 精品久久欧美熟妇www| www插插插无码视频网站| 91蜜桃国产成人精品区在线| 亚洲人成网站色在线入口口| 国内成+人 亚洲+欧美+综合在线| 久久99精品久久久久久野外| 网友自拍人妻一区二区三区三州| 91精品国产91综合久久蜜臀| 日本一本之道高清不卡免费| 伊人久久成人成综合网222| 特黄三级一区二区三区| 亚洲精品中文字幕视频色| 中文人妻熟妇乱又伦精品| 亚洲欧美日韩国产色另类| 亚洲成人黄色av在线观看| 国产毛片黄片一区二区三区| 成人免费一区二区三区| 亚洲—本道中文字幕久久66| 亚洲女同性恋激情网站| 免费国产在线精品一区| 欧美成人一区二区三区在线观看| 中文亚洲成a人片在线观看| 亚洲精品在线免费视频| 日韩中文字幕免费视频| 中日韩欧美在线观看| 亚洲国产人成自精在线尤物| 国产精品亚洲а∨无码播放| 欧美人与动牲交a欧美精品| 91精品国产91久久综合桃花| 国产黄久色一区2区三区| 日产学生妹在线观看| 日本欧美在线播放| 色综合久久人妻精品日韩| 亚洲色图片区| 无码国产精品一区二区免费16| 亚洲妇女av一区二区| 情爱偷拍视频一区二区|