本發(fā)明涉及機(jī)械振動(dòng)信號(hào)和聲輻射信號(hào)處理領(lǐng)域，具體涉及一種基于同步壓縮短時(shí)傅里葉變換的盲分離混合矩陣估計(jì)方法。

背景技術(shù)：

機(jī)械系統(tǒng)的振動(dòng)和噪聲對(duì)系統(tǒng)的性能與安全有重要影響，如潛艇推進(jìn)系統(tǒng)振動(dòng)過(guò)大容易被敵方發(fā)現(xiàn)且干擾自身的聲吶探測(cè)，火箭發(fā)動(dòng)機(jī)振動(dòng)過(guò)大容易導(dǎo)致其攜帶的衛(wèi)星及儀器不能正常工作等。因此，尋找振動(dòng)噪聲源并采取措施減小其影響至關(guān)重要。傳感器檢測(cè)的觀察信號(hào)往往是多個(gè)振動(dòng)源疊加，而且鑒于機(jī)械系統(tǒng)的復(fù)雜性，通常缺乏振動(dòng)源混合特性的先驗(yàn)知識(shí)，給振動(dòng)噪聲源的尋找?guī)?lái)了挑戰(zhàn)。盲源分離能夠在未知源信號(hào)和混合方式的情況下，對(duì)源信號(hào)進(jìn)行估計(jì)，很好地解決了以上難題。

然而，在工程實(shí)際中，傳感器的數(shù)目經(jīng)常會(huì)少于振動(dòng)源的數(shù)目，如下面幾種情況：(1)傳感器安裝數(shù)目較少，某些裝備不能隨意安裝傳感器，如火箭和導(dǎo)彈等一旦型號(hào)確定，不能隨意更改結(jié)構(gòu)方案；(2)傳感器長(zhǎng)期工作在惡劣環(huán)境中引起故障，沒(méi)有采集到數(shù)據(jù)，在數(shù)據(jù)缺失的情況下，重新采集數(shù)據(jù)費(fèi)時(shí)費(fèi)力；(3)數(shù)據(jù)可用性差，傳感器安裝的位置不合理，對(duì)信號(hào)振動(dòng)不敏感。此時(shí)，傳感器采集的觀測(cè)信號(hào)數(shù)目少于源信號(hào)數(shù)目，屬于欠定情況，傳統(tǒng)的盲源分離方法(如快速獨(dú)立分量分析，自然梯度算法等)效果較差，因此研究在欠定情況下混合矩陣的估計(jì)方法具有重要的學(xué)術(shù)意義和工程價(jià)值。

目前，欠定情況下混合矩陣的估計(jì)方法主要利用信號(hào)在時(shí)域、頻域以及時(shí)頻域的稀疏性提取單源點(diǎn)，進(jìn)而對(duì)單源點(diǎn)進(jìn)行聚類(lèi)分析，可以得到對(duì)應(yīng)于每個(gè)源的單源點(diǎn)，即對(duì)應(yīng)于混合矩陣的每一列。然而對(duì)于機(jī)械系統(tǒng)而言，很難保證所有的振動(dòng)信號(hào)在時(shí)域和頻域上均具有稀疏性，而基于時(shí)頻域的混合矩陣估計(jì)方法又多為基于stft和小波變換等，這些傳統(tǒng)時(shí)頻處理方法得到的時(shí)頻變換能量聚集性低，造成復(fù)雜機(jī)械信號(hào)在時(shí)頻域上稀疏性不足，進(jìn)而導(dǎo)致單源點(diǎn)估計(jì)具有較大誤差，在低信噪比條件下，甚至不能正確估計(jì)混合矩陣。

技術(shù)實(shí)現(xiàn)要素：

為了解決現(xiàn)有技術(shù)中的問(wèn)題，本發(fā)明提出一種基于同步壓縮短時(shí)傅里葉變換的盲分離混合矩陣估計(jì)方法，能夠提高觀察信號(hào)時(shí)頻域的稀疏性，提高了有用信號(hào)的能量，進(jìn)而減弱噪聲的影響，正確估計(jì)混合矩陣。

為了實(shí)現(xiàn)以上目的，本發(fā)明所采用的技術(shù)方案為：包括以下步驟：

步驟1，將所有待分析的觀察信號(hào)進(jìn)行短時(shí)傅里葉變換得到相應(yīng)的時(shí)頻域復(fù)矩陣；

步驟2，將獲取的時(shí)頻域復(fù)矩陣對(duì)時(shí)移因子求偏導(dǎo)，并與時(shí)頻域復(fù)矩陣取比值，得到瞬時(shí)頻率估計(jì)算子矩陣，利用瞬時(shí)頻率估計(jì)算子矩陣對(duì)時(shí)頻域復(fù)矩陣在頻率方向上進(jìn)行積分重排，實(shí)現(xiàn)時(shí)頻域能量重排至能量重心，得到時(shí)頻域復(fù)矩陣的同步壓縮時(shí)頻矩陣；

步驟3，利用稀疏編碼技術(shù)在同步壓縮時(shí)頻矩陣中尋找來(lái)自同一個(gè)1-d子空間的時(shí)頻點(diǎn)，構(gòu)造一個(gè)基于l1范數(shù)的目標(biāo)函數(shù)，實(shí)現(xiàn)單源點(diǎn)快速提取；

步驟4，將已提取的單源點(diǎn)進(jìn)行層次聚類(lèi)，得到所有類(lèi)的聚類(lèi)中心，每個(gè)聚類(lèi)中心對(duì)應(yīng)混合矩陣的一列，進(jìn)而利用聚類(lèi)中心的方向?qū)崿F(xiàn)對(duì)混合矩陣的估計(jì)。

所述步驟1中對(duì)每個(gè)待分析的觀察信號(hào)進(jìn)行短時(shí)傅里葉變換，對(duì)第i個(gè)待分析的觀察信號(hào)xi＝[xi1xi2…xin]進(jìn)行短時(shí)傅里葉變換得到對(duì)應(yīng)的時(shí)頻域復(fù)矩陣為：

其中i＝1,2,…,m，n表示采樣點(diǎn)個(gè)數(shù)，在時(shí)頻域復(fù)矩陣中，行表示頻率，列表示時(shí)間，表示的第k列，k＝1,2,…,k。

所述步驟2中采用同步壓縮變換對(duì)時(shí)頻域復(fù)矩陣進(jìn)行后處理，首先將時(shí)頻域復(fù)矩陣對(duì)時(shí)移因子求偏導(dǎo)；

然后得到的偏導(dǎo)與時(shí)頻域復(fù)矩陣本身取比值，得到瞬時(shí)頻率估計(jì)算子矩陣其中，i＝1,2,…,m，表示將時(shí)頻域復(fù)矩陣對(duì)其時(shí)移因子u求偏導(dǎo)，為虛數(shù)單位，

最后利用瞬時(shí)頻率估計(jì)算子矩陣ωi對(duì)時(shí)頻域復(fù)矩陣在頻率方向上進(jìn)行積分，得到時(shí)頻域復(fù)矩陣的同步壓縮時(shí)頻矩陣：yi＝[yi,1yi,2…yi,k]，其中i＝1,2,…,m，yi,k表示yi的第k列。

所述步驟2中對(duì)同步壓縮時(shí)頻矩陣yi進(jìn)行重新排列，使得每個(gè)yi轉(zhuǎn)換成一個(gè)行向量，即其中i＝1,2,…,m，t表示轉(zhuǎn)置，則所有觀察信號(hào)的同步壓縮時(shí)頻矩陣組合成矩陣：

所述步驟3中同步壓縮時(shí)頻矩陣yi的每一列是由不同的觀察信號(hào)在同一時(shí)頻點(diǎn)的復(fù)數(shù)值組成的列向量，對(duì)應(yīng)于觀察信號(hào)在此時(shí)頻點(diǎn)的方向，令其中為的第q列，q＝1,2,…,q，則利用中的其他列進(jìn)行稀疏編碼，即s.t.cii＝0，其中為編碼系數(shù)。

所述步驟3中若某一時(shí)頻點(diǎn)為單源點(diǎn)，則此時(shí)頻點(diǎn)的編碼系數(shù)中只有一個(gè)元素非零，即同一個(gè)源對(duì)應(yīng)的單源點(diǎn)所對(duì)應(yīng)的時(shí)頻列向量位于同一個(gè)1-d子空間中，則對(duì)l0范數(shù)進(jìn)行優(yōu)化尋找單源點(diǎn)，在ci足夠稀疏時(shí)，l0范數(shù)優(yōu)化與l1范數(shù)優(yōu)化等價(jià)，l0范數(shù)為：min||ci||0,s.t.cii＝0，l1范數(shù)為：min||ci||1,s.t.cii＝0。

所述步驟3基于l1范數(shù)建立如下目標(biāo)函數(shù)：

s.t.cii＝0

利用該目標(biāo)函數(shù)尋找與編碼系數(shù)ci只有一個(gè)元素為非零值相對(duì)應(yīng)的矩陣中的列向量，則該列向量對(duì)應(yīng)的時(shí)頻點(diǎn)為單源點(diǎn)，其中λ是正則化參數(shù)。

所述尋找與編碼系數(shù)ci只有一個(gè)元素為非零值優(yōu)化為：|cij|＞ε1，cip＜ε2且cii＝0，其中ε1和ε2表示設(shè)定的閾值，i≠j，i≠p且j≠p。

所述步驟4中采用下列公式評(píng)價(jià)估計(jì)的混合矩陣的準(zhǔn)確性：

其中，ai為原始混合矩陣a的第i列，為估計(jì)的混合矩陣對(duì)應(yīng)的第i列，e為準(zhǔn)確性指標(biāo)，代表估計(jì)的混合矩陣和原始混合矩陣的相近程度。

與現(xiàn)有技術(shù)相比，本發(fā)明利用stft同步壓縮方法將待分析信號(hào)從時(shí)域變換到時(shí)頻域，極大地提高了待分析信號(hào)在時(shí)頻平面的能量聚集性，從而提高其在時(shí)頻域的稀疏性。作為本方法的前處理，同步壓縮短時(shí)傅里葉變換具有較強(qiáng)的魯棒性。因此本方法相比將傳統(tǒng)時(shí)頻分析作為前處理的方法，能夠得到具有更強(qiáng)的魯棒性的結(jié)果。本發(fā)明在提取單源點(diǎn)的過(guò)程中，利用稀疏編碼技術(shù)，尋找編碼系數(shù)中只有一個(gè)元素非零對(duì)應(yīng)的時(shí)頻點(diǎn)，實(shí)現(xiàn)位于同一個(gè)1-d子空間的時(shí)頻點(diǎn)的搜尋。通過(guò)最小化編碼系數(shù)的l1范數(shù)，構(gòu)造一個(gè)誤差最小目標(biāo)函數(shù)實(shí)現(xiàn)單源點(diǎn)的提取，能夠進(jìn)一步增強(qiáng)對(duì)噪聲的抑制作用，最后通過(guò)層次聚類(lèi)法實(shí)現(xiàn)混合矩陣的估計(jì)。本發(fā)明僅需觀察信號(hào)大于等于2個(gè)，即可實(shí)現(xiàn)任意源數(shù)目混合矩陣的估計(jì)，并且當(dāng)信噪比較低時(shí)，仍然保持良好的估計(jì)性能，具有較好的魯棒性。

附圖說(shuō)明

圖1為實(shí)施例1的仿真信號(hào)波形圖；

圖2為本實(shí)施例1的仿真信號(hào)的stft時(shí)頻圖；

圖3為實(shí)施例1的仿真信號(hào)的stft同步壓縮時(shí)頻圖；

圖4為實(shí)施例2的六個(gè)源信號(hào)波形圖；

圖5為實(shí)施例2的六個(gè)源信號(hào)經(jīng)過(guò)混合矩陣后得到的兩個(gè)觀察信號(hào)的混合信號(hào)波形圖；

圖6為實(shí)施例2不同信噪比下混合矩陣估計(jì)誤差對(duì)比圖；

圖7為實(shí)施例2不同信噪比下時(shí)間消耗對(duì)比圖；

圖8為實(shí)施例2不同單源點(diǎn)數(shù)量不同信噪比下混合矩陣估計(jì)誤差對(duì)比圖；

圖9為實(shí)施例2不同單源點(diǎn)數(shù)量不同信噪比下時(shí)間消耗對(duì)比圖。

具體實(shí)施方式

下面結(jié)合具體的實(shí)施例和說(shuō)明書(shū)附圖對(duì)本發(fā)明作進(jìn)一步的解釋說(shuō)明。

本發(fā)明包括以下步驟：

步驟1，將待分析的觀測(cè)信號(hào)x進(jìn)行短時(shí)傅里葉變換得到相應(yīng)的時(shí)頻域復(fù)矩陣

對(duì)每個(gè)待分析觀察信號(hào)xi(t)依次進(jìn)行短時(shí)傅里葉變換操作，其中x(t)＝[x1(t)x2(t)…xm(t)]^t為所有觀察信號(hào)，xi(t)表示x(t)的第i個(gè)觀察信號(hào)，m表示觀察信號(hào)的個(gè)數(shù)；對(duì)第i個(gè)待分析的觀察信號(hào)xi＝[xi1xi2…xin]進(jìn)行短時(shí)傅里葉變換得到對(duì)應(yīng)的時(shí)頻域復(fù)矩陣其中i＝1,2,…,m，n表示采樣點(diǎn)個(gè)數(shù)，在中，行表示頻率，列表示時(shí)間，表示的第k列，k＝1,2,…,k；

步驟2，將獲取的時(shí)頻域復(fù)矩陣對(duì)時(shí)移因子求偏導(dǎo)，并與取比值，得到瞬時(shí)頻率估計(jì)算子矩陣ωx，利用ωx對(duì)在頻率方向上進(jìn)行重排，實(shí)現(xiàn)時(shí)頻域能量重排至能量重心，得到的同步壓縮時(shí)頻矩陣y，提高x在時(shí)頻域上的稀疏性；

采用同步壓縮變換對(duì)步驟1中得到的時(shí)頻域復(fù)矩陣進(jìn)行后處理，將對(duì)時(shí)移因子求偏導(dǎo)，再與本身取比值，得到瞬時(shí)頻率估計(jì)算子其中i＝1,2,…,m，表示將對(duì)其時(shí)移因子求偏導(dǎo)，為虛數(shù)單位，

利用得到的瞬時(shí)頻率估計(jì)算子ωi對(duì)在頻率上進(jìn)行積分，得到同步壓縮后的復(fù)矩陣，即為的同步壓縮變換yi＝[yi,1yi,2…yi,k]，其中i＝1,2,…,m，yi,k表示yi的第k列；

將得到的yi進(jìn)行重新排列，使得每個(gè)yi轉(zhuǎn)換成一個(gè)行向量，即其中i＝1,2,…,m，t表示轉(zhuǎn)置，則所有觀察信號(hào)的同步壓縮變換組合成一個(gè)矩陣

步驟3，利用稀疏編碼技術(shù)，將單源點(diǎn)的尋找轉(zhuǎn)化為尋找來(lái)自同一1-d子空間的時(shí)頻點(diǎn)的問(wèn)題，構(gòu)造一個(gè)基于l1范數(shù)的目標(biāo)函數(shù)，實(shí)現(xiàn)單源點(diǎn)快速提取；

同步壓縮變換矩陣的每一列是由不同的觀察信號(hào)在同一時(shí)頻點(diǎn)的復(fù)數(shù)值組成的列向量，對(duì)應(yīng)于觀察信號(hào)在此時(shí)頻點(diǎn)的方向，令其中為的第q列，q＝1,2,…,q，則可以被中的其他列進(jìn)行稀疏編碼，即s.t.cii＝0，其中為編碼系數(shù)；

若某一時(shí)頻點(diǎn)為單源點(diǎn)(在此時(shí)頻點(diǎn)上僅存在一個(gè)源)，則此時(shí)頻點(diǎn)的編碼系數(shù)中只有一個(gè)元素非零，即同一個(gè)源對(duì)應(yīng)的單源點(diǎn)所對(duì)應(yīng)的時(shí)頻列向量位于同一個(gè)1-d子空間中，則尋找單源點(diǎn)的問(wèn)題轉(zhuǎn)化為以下l0范數(shù)優(yōu)化問(wèn)題：min||ci||0,s.t.cii＝0，此優(yōu)化問(wèn)題在ci足夠稀疏時(shí)，與下述的l1范數(shù)優(yōu)化問(wèn)題等價(jià)：min||ci||1,s.t.cii＝0；

為了提高本方法的穩(wěn)定性和魯棒性，構(gòu)造如下目標(biāo)函數(shù)：s.t.cii＝0，尋找與編碼系數(shù)ci只有一個(gè)元素為非零值相對(duì)應(yīng)的中的列向量，其中λ是正則化參數(shù)，則該列向量對(duì)應(yīng)的時(shí)頻點(diǎn)為單源點(diǎn)，考慮到噪聲的影響，尋找編碼系數(shù)ci只有一個(gè)元素為非零值優(yōu)化為|cij|＞ε1，cip＜ε2且cii＝0，其中ε1和ε2表示設(shè)定的閾值，i≠j，i≠p且j≠p；

步驟4，將已提取的單源點(diǎn)進(jìn)行層次聚類(lèi)，得到估計(jì)的聚類(lèi)中心，每個(gè)類(lèi)的中心對(duì)應(yīng)混合矩陣的一列，進(jìn)而利用聚類(lèi)中心的方向?qū)崿F(xiàn)混合矩陣的估計(jì)；

利用步驟3的單源點(diǎn)進(jìn)行層次聚類(lèi)，可獲得所有類(lèi)的聚類(lèi)中心，每個(gè)聚類(lèi)中心與混合矩陣的某一列相對(duì)應(yīng)，通過(guò)聚類(lèi)中心的方向可得到混合矩陣每列的估計(jì)，從而得到估計(jì)的混合矩陣

為了評(píng)價(jià)本方法估計(jì)混合矩陣的準(zhǔn)確性，利用下面的性能指標(biāo)：

其中ai為原始混合矩陣a的第i列，為估計(jì)的混合矩陣對(duì)應(yīng)的列，e為準(zhǔn)確性指標(biāo)，代表估計(jì)的混合矩陣和原始混合矩陣的相近程度，即估計(jì)的準(zhǔn)確性。

實(shí)施例1，為了說(shuō)明同步壓縮變換相對(duì)于傳統(tǒng)時(shí)頻變換的優(yōu)越性，選取一個(gè)仿真信號(hào)，其波形圖如圖1所示。仿真信號(hào)的表達(dá)式為：f(t)＝cos{2π[0.1t²+3sin(2t)+10t]}+e^-0.2tcos[2πt(40+t^1.3)]，采樣時(shí)間為10s，采樣頻率為200hz。從表達(dá)式中看出其由兩部分組成，一部分的頻率為0.1t+[3sin(2t)]/t+10，另一部分的頻率為40+t^1.3，兩部分的頻率都隨著時(shí)間的變化而變化。

將該仿真信號(hào)進(jìn)行stft變換，得到的時(shí)頻圖如圖2所示，其頻率為兩部分，但是由于時(shí)頻聚集性不高，可讀性不強(qiáng)。

圖3所示為該仿真信號(hào)經(jīng)過(guò)stft同步壓縮變換得到的時(shí)頻圖，與圖2相比，時(shí)頻聚集性得到明顯提高，而且由于頻率上的能量重排，增加了有用信號(hào)的能量。與圖2相比，圖3中的最大能量由0.7增加到了4。與單獨(dú)的stft相比，stft同步壓縮變換極大地提高了時(shí)頻平面的稀疏性。

實(shí)施例2，選擇6個(gè)語(yǔ)音信號(hào)作為源信號(hào)，其波形圖如圖4所示，雖然時(shí)域上有一定的稀疏性，但是6個(gè)語(yǔ)音信號(hào)在大部分時(shí)間內(nèi)重疊，即在大多數(shù)時(shí)間點(diǎn)存在多個(gè)語(yǔ)音信號(hào)。選擇原始混合矩陣如下：a＝[0.2588,0.7071,0.9659,0.9659,0.7071,0.2588；-0.9659，-0.7071，-0.2588，0.2588，0.7071,0.9659]。

圖5所示為實(shí)施例2中6個(gè)語(yǔ)音信號(hào)經(jīng)過(guò)原始混合矩陣a后得到的兩個(gè)觀察信號(hào)，即，此時(shí)相當(dāng)于共有2個(gè)傳感器，而源信號(hào)的個(gè)數(shù)為6，為欠定盲源分離情況，一般的超定或者適定盲源分離分離方法失效。

在snr＝40db，窗長(zhǎng)為1024，窗每次移動(dòng)512采樣點(diǎn)條件下，利用本方法估計(jì)的混合矩陣為：和原始混合矩陣相比，估計(jì)混合矩陣與其絕對(duì)差值為：每個(gè)元素的絕對(duì)誤差均小于0.05，表明本方法具有良好的估計(jì)性能。

為了評(píng)價(jià)本發(fā)明方法估計(jì)混合矩陣的性能，利用下述指標(biāo)：

其中ai為原混合矩陣a的第i列，為本方法估計(jì)混合矩陣對(duì)應(yīng)的列。

圖6所示為實(shí)施例2中本發(fā)明方法在不同信噪比條件下的誤差對(duì)比圖，其參數(shù)如下：窗長(zhǎng)為1024，窗每次移動(dòng)采樣點(diǎn)數(shù)為512，每個(gè)數(shù)據(jù)均為相同參數(shù)下重復(fù)10次得到的平均值。圖6表明在不同的信噪比下，算法的性能誤差均小于0.03，表明本方法對(duì)噪聲具有較好的魯棒性。

圖7為實(shí)施例2中本方法在不同信噪比條件下所消耗的時(shí)間對(duì)比，參數(shù)和圖6相同，顯示每次消耗的時(shí)間在4s～6s，表明了本方法的高效性。

圖8所示為實(shí)施例2中提取不同的單源點(diǎn)數(shù)量在不同的信噪比下誤差對(duì)比圖，參數(shù)和圖6相同，表明隨著單源點(diǎn)數(shù)量的增加，本方法的誤差減小，而單源點(diǎn)數(shù)量增加到一定數(shù)量后，單源點(diǎn)數(shù)量對(duì)誤差的影響減小，但是在低信噪比條件下，單源點(diǎn)數(shù)量的增加，誤差反而有略微的增大。

圖9所示為實(shí)施例2中提取不同的單源點(diǎn)數(shù)量在不同的信噪比下消耗時(shí)間對(duì)比圖，參數(shù)和圖6相同，表明隨著單源點(diǎn)數(shù)量的增加，本方法消耗的時(shí)間逐漸增加，一方面是由于提取單源點(diǎn)數(shù)量增加導(dǎo)致的，另一方面是由于單源點(diǎn)數(shù)量增加，在聚類(lèi)時(shí)也需要消耗更多的時(shí)間。通過(guò)大量的試驗(yàn)和仿真驗(yàn)證，單源點(diǎn)的數(shù)量在300～1000時(shí)，本方法即可具有較高的運(yùn)算效率，并且具有較小的估計(jì)誤差。

本發(fā)明通過(guò)對(duì)觀察信號(hào)進(jìn)行stft同步壓縮變換，使之由時(shí)域變換到時(shí)頻域，通過(guò)在頻率上重排各時(shí)頻點(diǎn)的能量，極大地提高了待分析信號(hào)在時(shí)頻平面的能量聚集性，從而提高時(shí)頻域的稀疏性并減弱噪聲的影響。在提取單源點(diǎn)的過(guò)程中，利用稀疏編碼方法，尋找編碼系數(shù)中只有一個(gè)元素非零對(duì)應(yīng)的時(shí)頻點(diǎn)，實(shí)現(xiàn)位于同一個(gè)1-d子空間的時(shí)頻點(diǎn)的搜尋。通過(guò)最小化編碼系數(shù)的l1范數(shù)，構(gòu)造一個(gè)誤差最小目標(biāo)函數(shù)實(shí)現(xiàn)單源點(diǎn)的提取，進(jìn)一步增強(qiáng)了對(duì)噪聲的抑制作用，最后通過(guò)層次聚類(lèi)法實(shí)現(xiàn)混合矩陣的估計(jì)。本方法是一種新的欠定盲源分離混合矩陣估計(jì)方法，僅需觀察信號(hào)大于等于2個(gè)，即可實(shí)現(xiàn)任意源數(shù)目混合矩陣的估計(jì)，并且當(dāng)信噪比較低時(shí)，仍然保持良好的估計(jì)性能，具有較好的魯棒性，為實(shí)際振動(dòng)噪聲源分離提供有效支持。