基于三階統(tǒng)計量的欠定盲源分離方法

鄒 亮 張 鵬 陳 勛

①(中國礦業(yè)大學信息與控制工程學院 徐州 221116)

②(中國科學技術大學信息科學技術學院 合肥 230026)

1 引言

盲源分離(Blind Source Separation, BSS)憑借著其強大的技術優(yōu)勢，在生物醫(yī)學信號處理、移動通信、地質(zhì)勘探、語音信號處理和故障檢測等領域得到了廣泛的運用[1–5]。它旨在將無法直接觀測到的源信號從被干擾和噪聲污染的混合信號中分離出來。所謂的“盲”有兩重含義：原始的源信號不能被直接觀測到；源信號如何混合也是未知的。以語音信號為例，當多個聲源同時存在時，如果忽略麥克風自身的非線性特性和信號的反射，麥克風的輸入信號為各個聲源的線性混合。以盲源分離為基礎的信號處理技術在大量實驗和應用中被證實可以有效分離出源信號，便于對信號的進一步分析。

根據(jù)源信號和觀測信號數(shù)目的相對多少，盲源分離分為超定、正定、欠定3種情況。超定和正定情況下，獨立成分分析(Independent Component Analysis, ICA)是最經(jīng)典的盲源分離算法之一[6]。其中，F(xiàn)astICA采用擬牛頓法實現(xiàn)參數(shù)尋優(yōu)，收斂速度較梯度下降算法明顯加快[7]。受信號采集條件和成本的限制，人們嘗試使用較少的傳感器實現(xiàn)對信號的檢測，欠定盲源分離(Underdetermined Blind Source Separation, UBSS)應運而生。大部分欠定盲源分離算法從信號的稀疏性出發(fā)，在時域或變換域(如短時傅里葉變換、小波變換等)中尋找單數(shù)據(jù)源點估計混合矩陣進而實現(xiàn)源信號的恢復[8–10]。較為經(jīng)典的算法包括基于時頻域比率的TIFROM(TIme Frequency Ratio Of Mixtures)和退化分離估計技術DUET(Degenerate Unmixing Estimation Technique)[11,12]。在這兩種算法的基礎上，衍生出了諸多欠定盲源分離算法，進一步放松了對源信號稀疏性的要求。文獻[13]基于對觀測信號傅里葉變換的實、虛部系數(shù)的比較，提出了一種新穎的單信號源點(Single Source Points, SSP)檢測方法，通過聚類分析估計混合矩陣。文獻[14]通過檢測在時頻域中對應點在復平面上是否共線實現(xiàn)SSP的篩選。文獻[15]通過對觀測信號的稀疏編碼，用L1正則化約束條件，提高了算法的魯棒性。文獻[16]利用特定源信號對應的單信號源點距離較小的特點，篩選出各個源信號對應的單信號源點，進而利用K-means實現(xiàn)混合矩陣的估計。然而，基于信號的稀疏性的欠定盲源分離算法往往對噪聲比較敏感。許多研究人員嘗試利用觀測信號的統(tǒng)計特性解決欠定盲源分離問題。文獻[17–19]通過對源信號的2階統(tǒng)計特性進行分析，利用張量的典范雙峰分解(canonical polyadic decomposition)和塔克分解(Tucker decomposition)，開發(fā)了一系列的欠定盲源分離算法。其中文獻[18]提出從多數(shù)據(jù)集角度解決欠定盲源分離問題，利用多數(shù)據(jù)集源信號的2階相關性構建3階張量，進而實現(xiàn)混合矩陣估計。然而，高斯噪聲等對稱分布噪聲的2階累積量不為零，影響了混合矩陣估計的準確率。

高斯信號所有關于分布的信息都包含在1階與2階矩以內(nèi)，2階以上的累積量恒等于0。為了更好地抑制對稱分布噪聲的影響，本文提出基于3階累積量和源信號自相關特性的欠定盲源分離方法，并將所提出的算法與經(jīng)典欠定盲源分離算法進行比較。在語音信號上的1000次蒙特卡羅實驗發(fā)現(xiàn)，本文提出的算法能夠更好地抑制噪聲對混合矩陣估計的影響。對本文感興趣的讀者可參考相關源代碼https://github.com/usefulbbs/thirdOrder。

2 信號模型

瞬時混合的盲源分離模型為

其中，X(t)=[x1(t),x2(t),...,xM(t)]T表 示M個傳感器接收的觀測信號，S(t)=[s1(t),s2(t),...,sN(t)]T表示N個潛在的源信號，A為M×N維的混合矩陣，E(t)表示信號采集過程中的噪聲。通常假設噪聲是均值為0且與源信號不相關的高斯信號。當M N時，式(1)表示欠定情形下的信號線性混合模型，其實質(zhì)是利用較少的觀測信號實現(xiàn)混合矩陣的估計和源信號的恢復。

在解決欠定盲源分離問題時，通常需要對問題的模型設定一定的假設條件，如信號的時頻域稀疏性、源信號的統(tǒng)計獨立性等，然后采用先估計混合矩陣后分離源信號的“兩步法”實現(xiàn)盲源分離。鑒于源信號稀疏性假設在現(xiàn)實情形中較難滿足，本文提出了利用混合信號3階統(tǒng)計量的研究思路，基于以下3個假設：

假設(1)：當i?=j時 ，si(t)和sj(t)不相關；

假設(2)：源信號為均值為零的非高斯信號；

假設(3)：混合矩陣A列向量中任意兩個列向量不共線。

近年來，世界主要國家及地區(qū)開始逐漸禁止塑料微珠在沖洗類化妝品等日化產(chǎn)品中的使用。究其原因，塑料微珠作為新型生態(tài)污染源對生態(tài)環(huán)境的影響及危害受到了業(yè)界的高度重視，也引起了熱議。為此，筆者收集并整理了塑料微珠的定義及其危害，以及近年來相關國家及地區(qū)對塑料微珠的一系列監(jiān)管動態(tài)，幫助企業(yè)了解及應對塑料微珠禁令；同時，也希望中國化妝品監(jiān)管機構積極應對并做出調(diào)整，順應塑料微珠禁用全球化的趨勢。

3 混合信號的3階累積量

對均值為零的觀測信號，其3階統(tǒng)計量可以表示為

其中，觀測信號的3階統(tǒng)計量中任何一個元素π(i,j,k)可以表示為

其中，⊙為哈達瑪積，Ai,:,Aj,:,Ak,:表示混合矩陣A的 第i,j,k行向量，其中i,j,k=1,2,...,M，3階張量π的尺寸為M×M×M。τ1,τ2分別表示對觀測信號的延時，×n表示張量的n-模乘積。為了簡便計算，文中令τ1=τ2。D=cum[S(t),S(t+τ),S(t+τ)]為源信號的3階累積量，其尺寸為N×N×N。由假設(1)和假設(2)可知D為超對角張量

通過改變延時τ的大小，構建了Q個不同的3階張量，對應著各個延時，具體可以表示為

其 中，τq為 觀 測 信 號 延 時 大 小，張 量D(τq)=cum(S(t),S(t+τq),S(t+τq))為超對角張量，其中q=1,2,...,Q。 本文將張量π(1),π(2),...,π(Q)堆疊在一起，組成一個4階張量W∈RM×M×M×Q，其公式滿足Wi,j,k,q=(π(q))i,j,k,i=1,2,...,M,j=1,2,...,M,k=1,2,...,M,q=1,2,...,Q。定義矩陣H∈RQ×N，其中的元素為Hq,n=，其中q=1,2,...,Q，n=1，2,...,N，A:,n表 示 為 混 合 矩 陣 的 第n列 列 向量，H:,n為 矩陣H的第n列 列向量。4階張量W可以表示為

從而混合矩陣的估計問題，可以轉(zhuǎn)化為求解一個凸優(yōu)化問題，目標函數(shù)為

4 信號恢復

對于超定或正定盲源分離問題，可以通過求混合矩陣的逆矩陣恢復源信號。然而，當觀測信號數(shù)目少于信號源數(shù)目時，混合矩陣是病態(tài)矩陣，不能通過直接求逆來恢復源信號。受文獻[14]啟發(fā)，本文采用廣義高斯模型對信號進行恢復。

對于任何欠定的非齊次線性方程，其解可以表示為通解與特解的和。對于本文的欠定盲源分離問題，方程AS=X的 解S可以表示為

其中，Sp表 示方程的特解，Sh表示方程的通解。特解Sp可以表示為

這里矩陣(A)?表 示矩陣A的偽逆。齊次方程的通解可以寫成

其中，V是一個大小為N×(N ?M)矩陣，該矩陣的列向量為矩陣A零空間的基。矩陣Z則是大小為(N ?M)×T的任意矩陣(T表示每個觀測信號樣本點的個數(shù))?；仃嘨可以根據(jù)混合矩陣A計算得到，因此估計N維的源信號的問題轉(zhuǎn)化為計算N ?M維潛在變量Z。本文采用廣義高斯模型來對源信號的概率密度函數(shù)進行建模，廣義高斯模型可以表示為

其中，Γ(·) 為 伽馬函數(shù)，u和σ分別為單個信號的均值和方差。通常人們關注的主要是源信號的變化而非信號的幅度，此處假設各個源信號的均值u為0。β為源信號峰態(tài)的控制參數(shù)，σ為源信號的方差。本文采用期望最大算法(Expectation Maximization,EM)估計模型的隱變量。本文重點關注欠定情形下混合矩陣的估計，源信號恢復的具體細節(jié)請參考文獻[14]。

基于3階統(tǒng)計量的欠定盲源分離算法的主要步驟如表1所示。本文默認選擇最小延時為0.0625 ms，延時間隔為0.125 ms，延時數(shù)目Q為10。

表1 本文算法步驟

5 仿真實驗和結果分析

5.1 性能評價指標

本文采用歸一化均方誤差(Normalized Mean Square Error, NMSE)和平均絕對皮爾遜相關系數(shù)(Mean Absolute Pearson Correlation Coefficient,MAPCC)來衡量混合矩陣的估計精度和源信號的估計精度[13,18]，其數(shù)學表達式分別如式(15)

其中，sn表示第n個通道源信號的真實值，s?n表示恢復出的sn，σ表示信號的標準差，c ov()表示信號之間的協(xié)方差。源信號真實值與恢復結果之間的相關系數(shù)越高表示信號恢復越準確。MAPCC為源信號真實值之間PCC絕對值的平均數(shù)，可以表示為

其中，MAPCC的取值范圍為[0,1]。由于噪聲等因素的影響，很難完全恢復源信號。通常認為，MAPCC大于0.8時表明信號恢復的性能較好。文中所有實驗的計算環(huán)境配置如下：處理器Intel(R)Core(TM) i5-10200H CPU @2.40 GHz, 16 GB內(nèi)存，計算軟件MATLAB 2019b。

5.2 仿真實驗

本文選擇4通道語音信號作為源信號(如圖1)，采樣頻率為16000 Hz，采樣點長度為160000。

圖1 4個語音源信號

經(jīng)過3× 4的混合矩陣混合成3路觀測信號(如圖2)，其中隨機產(chǎn)生的混合矩陣為

圖2 3路觀測信號

使用本文提出的基于3階統(tǒng)計量的算法對上述欠定盲源分離問題進行求解，求得混合矩陣

評估發(fā)現(xiàn)，估計誤差為–27.4030 dB。恢復的源信號如圖3所示，源信號的估計值和真實值之間的皮爾遜相關系數(shù)的絕對值為0.9757，恢復的源信號和真實源信號十分接近。

圖3 恢復的源信號

5.3 性能比較

本文在觀測信號的信噪比取不同值時，與常見的欠定盲源分離算法的性能進行了比較，進一步驗證了本文提出的算法的有效性。為克服混合矩陣的隨機性對實驗結果的影響，通過1000次蒙特卡羅實驗，產(chǎn)生1000個3×4的混合矩陣，并將源信號混合成1000組觀測信號，計算相關算法性能的平均值。其中，與原文中參數(shù)設置一致，文獻[13]中參數(shù)?θ=0.8?，文獻[17]和文獻[18]中延時步長為0.125 ms，延時數(shù)量為10。1000次蒙特卡羅對比實驗結果如圖4所示，從圖中可以看出，隨著信噪比的降低，6種算法對混合矩陣估計的準確率逐漸提高?；?階或高階統(tǒng)計量的估計方法，在信噪比大于15 dB后性能基本穩(wěn)定，例如本文所提方法在信噪比為15 dB和30 dB時，對混合矩陣的估計誤差分別為–20.3548 dB和–20.0518 dB?；谛盘栂∈栊缘那范ぴ捶蛛x算法在信噪比大于10 dB后，隨著信噪比增加性能迅速提升，例如文獻[14]算法在10 dB, 20 dB和30 dB時，對混合矩陣的估計誤差分別是–5.3294 dB, –11.6122 dB和–17.7624 dB。主要原因是在當混合信號受到噪聲干擾時，信號的稀疏性受到較大影響，信號稀疏性的假設不再成立。然而，在現(xiàn)實的信號采集過程中，信號極易受到噪聲影響。以腦電信號為例，其強度只有微伏級，很容易受到肌電、眼電等非腦信號的干擾。本文提出的算法具有較強的抗干擾能力，能夠準確地估計混合矩陣。本文的重點為混合矩陣的估計，而信號恢復可參考期望最大算法、稀疏編碼等方法[14,16]。鑒于此，基于上述混合矩陣的估計，本文使用廣義高斯模型和期望最大算法進行信號恢復，結果如表2所示。當信噪比較低時，平均皮爾遜相關系數(shù)較小。例如，當信噪比為–5 dB時，提取出的源信號與真實值有較大差別，MAPCC在0.3附近，所恢復出的信號無使用價值。隨著信噪比的升高，信號恢復效果逐漸變好。本文所提算法在信噪比為30 dB時，MAPCC可達0.88，與其他算法相比，取得了最優(yōu)性能。

表2 各算法恢復出的源信號與真實源信號間的平均絕對皮爾遜相關系數(shù)

圖4 欠定盲源分離算法混合矩陣估計的歸一化均方誤差比較

為了探究本文算法的性能隨觀測信號數(shù)目和源信號數(shù)目的變化規(guī)律，本文比較了當觀測信號數(shù)目為3時，在源信號數(shù)目不同時算法的性能及當源信號數(shù)目為7時，在觀測信號數(shù)目不同時算法的性能，結果分別如圖5(a)和圖5(b)所示。當觀測信號數(shù)目一定時，源信號數(shù)目越多，混合矩陣中待估計參數(shù)越多，估計難度越大。當源信號數(shù)目一定時，隨著觀測信號數(shù)目增多，用于估計混合矩陣的統(tǒng)計信息越充足，估計越準確。

圖5 源信號和觀測信號數(shù)目變化時，算法的性能

為了探究算法的時效性，本文進行了1000次蒙特卡羅實驗，隨機產(chǎn)生1000個混合矩陣，并選擇在信噪比為30 dB時計算各算法平均運行1次所需要的時間，其結果如表3所示。本文算法平均運行1次僅需要1.2 ms，顯著少于文獻[13–16]中的算法，與文獻[17,18]所需時間相似。文獻[13–16]基于信號的稀疏性實現(xiàn)混合矩陣的估計，對應算法的耗時主要取決于篩選單數(shù)據(jù)源點所需的時間。其中，文獻[13,14]根據(jù)信號在時頻域中對應時頻點(time-frequency point)的實部和虛部向量是否共線，篩選單信號源點，該類算法所消耗的時間和信號長度成正比；文獻[15,16]計算信號在時頻域中每個點與其他點的距離，通過合理設置閾值，篩選出單信號源點，該類算法所消耗的時間和信號長度的平方成正比。文獻[17,18]和本文算法則通過觀測信號的統(tǒng)計量構建張量，并通過張量分解求解出混合矩陣，該類算法中主要是張量分解比較耗時，受信號長度影響不大。

表3 各種算法平均運行1次所需要的時間(s)

6 結束語

針對欠定情形下的混合矩陣估計問題，研究人員多通過稀疏理論來展開，而忽略了噪聲自身的統(tǒng)計特性及其對信號稀疏性的影響。本文圍繞欠定盲源分離中的混合矩陣估計，提出了一種基于信號3階統(tǒng)計特性的欠定盲源分離算法。利用源信號的自相關特性，構建3階張量，并將不同時延對應的3階堆砌成4階張量，進而通過典范雙峰分解和奇異值分解實現(xiàn)混合矩陣的估計。1000次蒙特卡羅實驗表明，本文所提算法可以有效抑制噪聲的影響、克服對信號稀疏性的約束，提升混合矩陣的估計精度，并在此基礎上實現(xiàn)源信號恢復。然而，需要指出的是，源信號數(shù)目和混合矩陣的估計是當前欠定盲源分離研究的兩大挑戰(zhàn)?，F(xiàn)有的欠定盲源分離方法多假設源信號數(shù)目已知，而忽視了源信號數(shù)目對混合矩陣估計結果的影響，這無疑會造成實際運用時的難度。后續(xù)的研究中，作者擬引入貝葉斯網(wǎng)絡和期望最大算法等技術手段，以進一步提高源信號數(shù)目和混合矩陣估計的準確性。