郭 慧，付永慶，蘇東林，劉 焱

（1.哈爾濱工程大學 信息與通信工程學院，哈爾濱150001；2.北京航空航天大學 電磁兼容技術研究所，北京100191）

0 引 言

電磁兼容測試的主要環(huán)節(jié)之一是電磁干擾測試［1］。由于電子設備可能由多個工作模塊構成或測試時有多個設備同時工作，因此，電磁干擾測試系統(tǒng)可視為一個多輸入單輸出系統(tǒng)或多輸入多輸出系統(tǒng)，輸入量指的是潛在干擾源發(fā)射的電磁波，輸出量指的是傳感器接收到的混合信號，其輸出數量由傳感器的數目而定。由于是多輸入系統(tǒng)，故電磁干擾測試方法的輸出是多個信號的混合，無法對干擾源進行獨立測試，為后續(xù)的干擾診斷和整改帶來很大困難。為此，需要一種可以在僅有混合信號的情況下對其分離的方法。

目前，盲源分離（Blind sources separation，BSS）能夠很好地分離混合信號，并形成很多成熟的算法，如過 完 備ICA 算 法［2］、FastICA 算 法［3］、JADE算法［4］等，且已被廣泛應用于圖像處理［5］、生物 醫(yī) 學［6］、雷 達 與 通 信 系 統(tǒng)［7］、地 質 信 號 處理［8］、數據挖掘［9］等領域。但是這些算法要求源信號的數目是已知的，且混合信號數目必須大于等于源信號數。對于干擾信源數未知，并且測試設備有限的情況，這種要求很難實現(xiàn)。針對這一問題，本文作者提出了一種新的欠定盲源分離（Underdetermined blind sources separation，UBSS）方法，用于分離源數未知，且接收設備較少（最少可達2）情況下的干擾混合信號。該方法基于源信號的稀疏特性，分離的干擾信號主要針對一些由晶振等引起的諧波干擾。分離結果即為各個干擾源的發(fā)射電磁波，據此，可以對電子設備進行干擾診斷和整改。

1 問題描述

電子設備的電磁干擾測試是一個平穩(wěn)線性瞬時混疊的過程，基于UBSS的電磁干擾分離模型如圖1所示。

圖1 基于UBSS的電磁干擾分離示意圖Fig.1 Schematic diagram of EMI separation based on UBSS

當源信號具有很好的稀疏特性時，任意時刻的源信號只有在一個采樣點的取值遠離零，在其他采樣點的取值為零或接近于零［10］，此時，混合信號的散點圖具有直線聚類特點。直線的數目就是源信號的數目，直線的斜率由混合矩陣的列向量決定。

2 源數和混疊矩陣估計

上述假設是某一時刻只有一個源信號存在，但在某些時刻不排除有多個信號同時存在的情況，這些時刻的觀測點會影響混合矩陣的估計精度。

通常情況下，某一信號單獨存在的時刻點是呈區(qū)域分布的，即相鄰的多個時刻均只存在同一信號，因此，采用“領域比值法”來獲取單源主導區(qū)間。在此基礎上，采用Hough加窗法分別估計源信號的數目和混疊矩陣。

2.1 單源主導區(qū)間提取

兩個混合信號中只有一個源信號存在時，有：

在該信號的主導區(qū)域內相鄰時刻的觀測點比值均為該信號對應的混合矩陣列向量元素的比值，該比值是個相等的常數。當有兩個源信號存在時，有：

t＝1，2，…，T；p ＝1，2，…，n；q＝1，2，…，n

在這兩個信號的主導區(qū)域內，只有當a2p／a1p＝a2q／a1q時，相鄰時刻的觀測點比值才會相等，因此，可以通過鄰域比值來提取單源主導區(qū)間。

對于m 維觀測向量，設Num 為搜索域的大小（通常取3～5），k＝T／Num 表示搜索域的數目，其中 · 表示比·小的最大整數，Bk＝［b1，b2，…，bNum］T表示搜索域中觀測點的鄰域比值向量，其中，，元素bij＝xi（tj）／x1（tj），i ＝1，2，…，m－1；j ＝1，2，…，Num；ti＝1 ＋（k－1）Num，…，Num ＋（k－1）Num。求Bk的方差，如式（3）所示：

考慮到噪聲的影響，設一個與噪聲相關的閾值ε，若Bk≤ε，則認為該搜索域是單源主導區(qū)間，保留；否則，舍棄該搜索域觀測數目。本文中m 取2，即考慮有兩個觀測信號的情況。

2.2 源數和混疊矩陣估計

經過單源主導區(qū)間提取后的數據具有較好的稀疏特性，呈現(xiàn)直線聚類特點，由此引入Hough變換［11］，將觀測信號域內的散點變換到變換空間中的某些位置，變換式為：

式中：ρ為觀測點到原點的距離；θ為觀測點與原點形成的直線與x 軸的夾角。

在圖2 中，原空間中的5 個點經式（4）變換后，屬于同一條直線上的點在變換空間中變成一簇交于點（0，3）的正弦曲線，如實線所示。而不屬于該直線上的點在變換空間中是一條不經過交點的正弦曲線，如虛線所示。

圖2 Hough變換檢測直線Fig.2 Line detection by Hough transform

為了減小算法的復雜度，在進行Hough變換前，先對混合信號進行單位化處理，即：

式中：‖·‖2表示l2范數。

經過單位化處理后，所有觀測點到原點的距離ρ＝1。所有觀測點變換后，對變換域中的變換量θ（t）進行統(tǒng)計，引入函數：

式中：k＝1，2，…，π／h，表示直方圖的分區(qū)數目，其中h是Hough 變換的量化步長，k 越大，即h越小，變換量的分類精度越高 ；［· ］ 表示比·小的最大整數。則落入直方圖某一區(qū)域的變換量數目為：

得到直方圖，若直接搜索其峰值點，容易陷入局部最大，造成誤估計。為此，對直方圖進行加窗處理，處理后重新計算位置點的統(tǒng)計數目為：

對新的直方圖搜索峰值數即為干擾源數，峰值對應的角度估計值與混合矩陣元素有式（9）的關系，由此可得到混合矩陣的估計。

式中：a1i表示混和矩陣中第1行第i列的元素。

3 混合信號分離

分離混合信號滿足如下的約束條件：

針對兩個混合信號的情況，滿足任意時刻起主導作用的源信號數小于等于2的條件的混合信號可以分離為：

式中aj和ak構成了時刻t的分離矩陣。因此，混合信號的分離關鍵是尋找使分離信號總和最小時的分離矩陣。針對混合信號數為2的情況，提出一種“夾角差排序法”。由Hough加窗法獲得的混合矩陣列向量在散點圖中是一組單位化過原點的向量，如圖3所示。

圖3 混合矩陣列向量的散點圖Fig.3 Scatter plot column vector of mixing matrix

在散點圖中，混合矩陣列向量是單位圓上的一點，如a1、a2和a3，x是時刻t的混合信號，計算每個列向量及混合信號與x1軸的夾角的公式為：

式中：P 是列向量或混合信號；sign（P（2））表示取P向量中第二項的符號；dot表示內積；x1是橫軸，計算時代入（1，0）即可。

得到所有的夾角后，用列向量的夾角減去混合信號的夾角，將所有的夾角差從小到大排序，選擇零度角前后的夾角差對應的列向量，如果夾角差均大于或小于零度角，則選擇最小和最大夾角差對應的列向量，所選列向量即構成該時刻的分離矩陣，按照式（11）遍歷所有時刻，便可得到分離信號。

4 實驗分析

4.1 算法性能評價標準

采 用 “泛 化 交 擾 誤 差”（Generalized crosstalking error，GCE）和“互相關系數”（Cross correlation coefficient，CCC）分別評估混合矩陣的估計精度和混合信號的分離精度?；旌暇仃嘇和估計矩陣A＿est的“泛化交擾誤差”定義為：

式中：Π 表示所有n×n維可逆矩陣組成的集合，而且這些矩陣的每一列只有一個非零元素；A＿estB 表示A＿est與一個尺度矩陣和置換矩陣的乘積，以消除A＿est 的幅度不確定性和排序不確定性，當且僅當A＿est 和A 完全等價時，Err（A，A＿est）＝0，而Err（A，A＿est）的值越小，則表示兩矩陣越接近。

盲源分離的互相關系數度量定義如下：

式中：sj為源信號，為還原信號，i＝1，2，…，n。在該度量下，0≤ξij≤1，ξij 的值越接近1，說明與sj越相似；反之，ξij 越接近0，說明與sj越不相似。

4.2 數據預處理

首先，去除具有諧波干擾的測試信號的幅度影響，目的是使信號具有稀疏性。具體方法為：逐次判斷目標域（一般為3～5個）中采樣點的最小值，并將該值賦給目標域中其他采樣點，全部完成后得到的曲線即為測試信號的幅度包絡，認為測試信號減去幅度包絡得到的信號x～（t）具有一定的稀疏特性。

然后，對上述信號的各列數據求l2范數，去除‖X（t）‖2≤0.001的點，減小噪聲的影響，對剩余信號按式（5）進行單位化處理。

最后，對數據進行對稱化。為了保證聚類方向唯一，簡化計算，對數據進行如下處理：

式中：sign（·）表示取·的符號。

4.3 仿真數據分析

通過仿真得到信噪比為35dB下的兩個晶振波形圖，分別為6 MHz和10 MHz，如圖4（a）所示，可以看出在較寬的頻帶內源信號的采樣值大多數接近于10dBuV／m，幅度很平坦。隨機選取式（16）中的2×2的混合矩陣，得到混合信號如圖4（b）所示。

混合信號經預處理后，采用本文算法得到混合信號的分離結果，如圖4（c）所示。

對比分離信號和源信號，除了信號的幅度有區(qū)別外，諧波信號的波形信息基本恢復，證明本文算法對這類諧波信號具有很好的分離效果，且適用于正定盲源分離問題。同時，分離信號與源信號的互相關系數平均值可達0.9936，表明算法具有很高的準確性。

圖4 仿真數據分離Fig.4 Simulation data separation

為了表明本文算法在欠定情況下的不同處理階段均具有較好的準確性和抗噪聲性能，在2 ～40dB的噪聲范圍內，對2m3s模型（2 個混合信號，3個源信號）的晶振混合信號分別采用Hough加窗法、K 均值法［12］、模糊C 均 值法［13］、勢函數法［14］進行50次Monte Carlo運算，得到不同算法對混合矩陣估計精度的GCE－SNR 關系圖如圖5所示。當混合矩陣已知時，采用夾角排序法、l1范數法［15］、子空間投影法［16］對混合信號進行50次Monte Carlo運算，得到幾種算法對混合信號分離精度的CCC－SNR 關系圖如圖6所示。

圖5 幾種算法對混合矩陣估計精度的性能比較圖Fig.5 Performance comparison of estimating mixing matrix among proposed algorithm and other methods

從圖5中可以看出，采用本文提出的Hough加窗法得到的GCE值始終小于0.2，比K 均值法和模糊C 均值法得到的結果小得多，尤其在信噪比低的情況下，優(yōu)勢更突出。勢函數法性能較接近本文算法，但由于Hough加窗法是基于提取出的單源主導區(qū)間，因此具有更高估計精度。仿真結果說明本文算法還具有較好的抗噪聲性能和較高的估計精度。

圖6 幾種算法對混合信號分離精度的性能比較圖Fig.6 Performance comparison of separating mixing signals among proposed algorithm and other methods

圖6 中由于子空間投影法只對任意時刻同時存在的信號源數目少于觀測信號數目的情況效果較好，而l1范數法選擇分離矩陣的標準具有一定的局限性，因此，這兩種算法對混合信號的分離效果不如本文提出的夾角差排序法。采用夾角差排序法在較低信噪比下的互相關系數仍可達到0.9以上，表明本文算法具有較好的抗噪聲性能。

4.4 實測數據分析

以機載電臺的電磁干擾測試為例，為了突出算法可以分離欠定情況較嚴重的混合信號，在機載電臺3m 的范圍內放置3個頻率不同且正在工作的晶振模塊來模擬測試環(huán)境中存在的其他設備，在暗室中對其進行測試得到機載電臺的電磁干擾測試頻譜如圖7（a）所示。從圖7（a）中可以看出，干擾測試頻譜并不像仿真數據那樣幅度平坦，測試數據中伴隨著被試品的頻率響應，造成數據不滿足稀疏假設條件，因此，若采用本文算法需首先去除幅度影響，圖7（a）中紅色曲線即為兩次測試數據的幅度包絡，去除幅度包絡后的測試數據如圖7（b）所示，可以看出處理后的數據與圖4（b）中的數據特性類似，可以對其進行分離，分離結果如圖7（c）所示。

圖7 實測數據分離Fig.7 Measured data separation

從分離結果可知，測試數據中有3種諧波成分，分別為7.2、4.5和5.5 MHz，與實驗中放置的3個晶振模塊頻率一致，若屏蔽其中一個晶振模塊，如5.5 MHz的晶振，對機載電臺再次進行電磁干擾測試，得到整改后的測試頻譜，其與整改前的測試頻譜對比圖如圖8所示。

通過對比整改前后的波形圖，可以清楚地看出5.5 MHz的諧波成分已被消除，表明本文算法對設備的諧波電磁干擾診斷具有一定的可行性，可以為后續(xù)的電磁干擾診斷及整改提供依據。

圖8 整改前、后的測試信號波形圖Fig.8 Test signals before and after rectification

5 結束語

機載設備具有高度集成、數目繁多等特點，集成了很多種類的時鐘晶振、電源開關等易造成諧波干擾的器件，為此，基于欠定盲源分離技術提出適用于具有諧波特性的機載設備電磁干擾分離方法。首先，采用鄰域比值法提取信號單源主導區(qū)間；然后，在單源主導區(qū)間上采用Hough加窗法對電磁干擾源的數目和混合信道進行估計；最后，采用夾角差排序法對電磁干擾混合觀測信號進行分離，得到的分離信號即可作為后續(xù)干擾源分析診斷和故障整改的對象。該方法解決了傳統(tǒng)測試方法無法對多個同時工作的設備進行單獨測試和測試資源有限導致的混合信號無法分離兩大問題。仿真實驗和實測數據表明本文算法具有較好的分離效果和一定的可行性。

［1］Prasad Kodali V.工程電磁兼容［M］.陳淑鳳，高攸綱，蘇東林，周璧華譯.北京：人民郵電出版社，2006：5－8.

［2］Pham D T，Cardoso J F.Blind separation of instantaneous mixtures of non stationary sources［J］.IEEE Transactions on Signal Processing，2000，49（9）：1837－1848.

［3］Reyhani N，Ylipaaavalniemi J，Vigario R，et al.Consistency and asymptotic normality of FastICA and bootstrap FastICA［J］.Signal and Processing，2012，92（8）：1767－1778.

［4］Cardoso J F，Souloumiac A.Blind beamforming for non－Gaussian signals［J］.IEEE Proceedings－F，1993，140（6）：362－370.

［5］Badawi W K M，Chibelushi C C，Patwary M N.et al.Specular－based illumination estimation using blind signal separation techniques［J］.IET Image Processing，2012，6（8）：1181－1191.

［6］Mammone N，La Foresta F，Morabito F C.Automatic artifact rejection from multichannel scalp EEG by wavelet ICA［J］.IEEE Sensors Journal，2012，12（3）：533－542.

［7］陳曉軍，成昊，唐斌.基于ICA 的雷達信號欠定盲分離算法［J］.電子與信息學報，2010，32（4）：919－924.Chen Xiao－jun，Chen Hao，Tang Bin.Underdetermined blind radar signal separation based on ICA［J］.Journal of Electronics ＆Information Technology，2010，32（4）：919－924.

［8］Ilin A，Valpola H，Oja E.Semiblind source separation of climate data detects E1 Nino as the component with the highest variability［C］∥Proceedings of International Joint Conference on Neural Networks，Montreal，Canada，2005：1722－1727.

［9］Takahata A K，Nadalin E Z，F(xiàn)errari R，et al.Unsupervised processing of geophysical signals：a review of some key aspects of blind deconvolution and blind source separation［J］.IEEE Signal Processing Magazine，2012，29（4）：27－35.

［10］Shindo H，Hirai Y.Blind source separation by ageometrical method［C］∥Proceeding of the International Joint Conference on Neural Networks，Honolulu，2002：1108－1114.

［11］郭斯羽，孔亞廣，張熙芳.基于Hough變換的角點檢測算法［J］.儀器儀表學報，2008，29（11）：2424－2429.Guo Si－yu，Kong Ya－guang，Zhang Xi－fang.Corner detection algorithm based on Hough transform［J］.Chinese Journal of Scientific Instrument，2008，29（11）：2424－2429.

［12］Onoda T，Sakai M，Yamada S.Careful seeding method based on independent components analysis for k－means clustering［J］.Journal of Emerging Technologies in Web Intelligence，2012，4（1）：51－59.

［13］Sun H J，Wang S R，Jiang Q S.FCM－based model selection algorithm for determining the number of cluster［J］.Pattern Recognition，2004，37（10）：2027－2037.

［14］Bofill P，Zibulevsky M.Underdetermined blind source separation using sparse representations［J］.Signal Processing，2001，81（11）：2353－2362.

［15］Li Y，Amari S I，Cichocki A.Underdetermined blind source separation based on sparse representation［J］.IEEE Transactions on Signal Processing，2006，54（2）：423－437.

［16］Aissa－El－Bey A，Linh－Trung N.Underdetermined blind separation of nondisioint sources in the timefrequency domain［J］.IEEE Transaction on Signal Processing，2007，55（3）：897－907.