基于欠定盲源分离的电磁干扰分离方法

郭 慧，付永庆，苏东林，刘 焱

（1.哈尔滨工程大学 信息与通信工程学院，哈尔滨150001；2.北京航空航天大学 电磁兼容技术研究所，北京100191）

0 引 言

电磁兼容测试的主要环节之一是电磁干扰测试［1］。由于电子设备可能由多个工作模块构成或测试时有多个设备同时工作，因此，电磁干扰测试系统可视为一个多输入单输出系统或多输入多输出系统，输入量指的是潜在干扰源发射的电磁波，输出量指的是传感器接收到的混合信号，其输出数量由传感器的数目而定。由于是多输入系统，故电磁干扰测试方法的输出是多个信号的混合，无法对干扰源进行独立测试，为后续的干扰诊断和整改带来很大困难。为此，需要一种可以在仅有混合信号的情况下对其分离的方法。

目前，盲源分离（Blind sources separation，BSS）能够很好地分离混合信号，并形成很多成熟的算法，如过 完 备ICA 算 法［2］、FastICA 算 法［3］、JADE算法［4］等，且已被广泛应用于图像处理［5］、生物 医 学［6］、雷 达 与 通 信 系 统［7］、地 质 信 号 处理［8］、数据挖掘［9］等领域。但是这些算法要求源信号的数目是已知的，且混合信号数目必须大于等于源信号数。对于干扰信源数未知，并且测试设备有限的情况，这种要求很难实现。针对这一问题，本文作者提出了一种新的欠定盲源分离（Underdetermined blind sources separation，UBSS）方法，用于分离源数未知，且接收设备较少（最少可达2）情况下的干扰混合信号。该方法基于源信号的稀疏特性，分离的干扰信号主要针对一些由晶振等引起的谐波干扰。分离结果即为各个干扰源的发射电磁波，据此，可以对电子设备进行干扰诊断和整改。

1 问题描述

电子设备的电磁干扰测试是一个平稳线性瞬时混叠的过程，基于UBSS的电磁干扰分离模型如图1所示。

图1 基于UBSS的电磁干扰分离示意图Fig.1 Schematic diagram of EMI separation based on UBSS

当源信号具有很好的稀疏特性时，任意时刻的源信号只有在一个采样点的取值远离零，在其他采样点的取值为零或接近于零［10］，此时，混合信号的散点图具有直线聚类特点。直线的数目就是源信号的数目，直线的斜率由混合矩阵的列向量决定。

2 源数和混叠矩阵估计

上述假设是某一时刻只有一个源信号存在，但在某些时刻不排除有多个信号同时存在的情况，这些时刻的观测点会影响混合矩阵的估计精度。

通常情况下，某一信号单独存在的时刻点是呈区域分布的，即相邻的多个时刻均只存在同一信号，因此，采用“领域比值法”来获取单源主导区间。在此基础上，采用Hough加窗法分别估计源信号的数目和混叠矩阵。

2.1 单源主导区间提取

两个混合信号中只有一个源信号存在时，有：

在该信号的主导区域内相邻时刻的观测点比值均为该信号对应的混合矩阵列向量元素的比值，该比值是个相等的常数。当有两个源信号存在时，有：

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

在这两个信号的主导区域内，只有当a2p／a1p＝a2q／a1q时，相邻时刻的观测点比值才会相等，因此，可以通过邻域比值来提取单源主导区间。

对于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≤ε，则认为该搜索域是单源主导区间，保留；否则，舍弃该搜索域观测数目。本文中m 取2，即考虑有两个观测信号的情况。

2.2 源数和混叠矩阵估计

经过单源主导区间提取后的数据具有较好的稀疏特性，呈现直线聚类特点，由此引入Hough变换［11］，将观测信号域内的散点变换到变换空间中的某些位置，变换式为：

式中：ρ为观测点到原点的距离；θ为观测点与原点形成的直线与x 轴的夹角。

在图2 中，原空间中的5 个点经式（4）变换后，属于同一条直线上的点在变换空间中变成一簇交于点（0，3）的正弦曲线，如实线所示。而不属于该直线上的点在变换空间中是一条不经过交点的正弦曲线，如虚线所示。

图2 Hough变换检测直线Fig.2 Line detection by Hough transform

为了减小算法的复杂度，在进行Hough变换前，先对混合信号进行单位化处理，即：

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

经过单位化处理后，所有观测点到原点的距离ρ＝1。所有观测点变换后，对变换域中的变换量θ（t）进行统计，引入函数：

式中：k＝1，2，…，π／h，表示直方图的分区数目，其中h是Hough 变换的量化步长，k 越大，即h越小，变换量的分类精度越高 ；［· ］ 表示比·小的最大整数。则落入直方图某一区域的变换量数目为：

得到直方图，若直接搜索其峰值点，容易陷入局部最大，造成误估计。为此，对直方图进行加窗处理，处理后重新计算位置点的统计数目为：

对新的直方图搜索峰值数即为干扰源数，峰值对应的角度估计值与混合矩阵元素有式（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和估计矩阵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）所示。

对比分离信号和源信号，除了信号的幅度有区别外，谐波信号的波形信息基本恢复，证明本文算法对这类谐波信号具有很好的分离效果，且适用于正定盲源分离问题。同时，分离信号与源信号的互相关系数平均值可达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 均值法得到的结果小得多，尤其在信噪比低的情况下，优势更突出。势函数法性能较接近本文算法，但由于Hough加窗法是基于提取出的单源主导区间，因此具有更高估计精度。仿真结果说明本文算法还具有较好的抗噪声性能和较高的估计精度。

图6 几种算法对混合信号分离精度的性能比较图Fig.6 Performance comparison of separating mixing signals among proposed algorithm and other methods

图6 中由于子空间投影法只对任意时刻同时存在的信号源数目少于观测信号数目的情况效果较好，而l1范数法选择分离矩阵的标准具有一定的局限性，因此，这两种算法对混合信号的分离效果不如本文提出的夹角差排序法。采用夹角差排序法在较低信噪比下的互相关系数仍可达到0.9以上，表明本文算法具有较好的抗噪声性能。

4.4 实测数据分析

以机载电台的电磁干扰测试为例，为了突出算法可以分离欠定情况较严重的混合信号，在机载电台3m 的范围内放置3个频率不同且正在工作的晶振模块来模拟测试环境中存在的其他设备，在暗室中对其进行测试得到机载电台的电磁干扰测试频谱如图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的谐波成分已被消除，表明本文算法对设备的谐波电磁干扰诊断具有一定的可行性，可以为后续的电磁干扰诊断及整改提供依据。

图8 整改前、后的测试信号波形图Fig.8 Test signals before and after rectification

5 结束语

机载设备具有高度集成、数目繁多等特点，集成了很多种类的时钟晶振、电源开关等易造成谐波干扰的器件，为此，基于欠定盲源分离技术提出适用于具有谐波特性的机载设备电磁干扰分离方法。首先，采用邻域比值法提取信号单源主导区间；然后，在单源主导区间上采用Hough加窗法对电磁干扰源的数目和混合信道进行估计；最后，采用夹角差排序法对电磁干扰混合观测信号进行分离，得到的分离信号即可作为后续干扰源分析诊断和故障整改的对象。该方法解决了传统测试方法无法对多个同时工作的设备进行单独测试和测试资源有限导致的混合信号无法分离两大问题。仿真实验和实测数据表明本文算法具有较好的分离效果和一定的可行性。

［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，Ferrari 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.