基于FastICA 的工频干扰消除算法

张君生，刘瑞安，王 磊

（天津师范大学物理与电子信息学院，天津 300387）

独立成分分析以非高斯源信号为研究对象，在统计独立的假设下，对多路观测到的混合信号进行盲分离，从而较完好地分离出隐含在混合信号中的源信号.近年来，通过ICA 实现盲信源分离已经逐渐成为信号处理领域的热点问题之一.从1990 年代以来，研究人员从不同角度提出多种ICA 算法.其中，Nadal 等[1]提出了信息最大化（Infomax）基本原理，Bell 等[2]给出该方法的随机梯度学习规则，并将其分别运用于预测、时间序列分析和盲源分离等研究中.Hyvärinen 等[3]提出并详细论述了最大似然估计法（MLE）.

应用过程中，ICA 算法所假设的条件，如信号之间的独立性、线性、瞬时性以及不考虑混合时产生的噪声等条件不一定能够得到满足，导致算法的应用效果并不理想.为解决这些问题，很多改进的ICA 方法应运而生，包括非线性ICA 法、约束ICA法和多分辨率子波段分解ICA 法等[4-9].除了对基本ICA 方法的改进外，基于信号时间结构和稀疏性的ICA 新算法不断被提出，奇异混合[10-11]、时延混合[12]和噪声ICA[13]等问题也吸引了众多学者的关注.

实际应用中，有用信号与噪声信号往往是独立的，本研究将较弱的有效地震信号和较强的工频干扰看成2 个独立的信源，接收到的信号是两者混合，采用FastICA 算法消除地震信号中的工频干扰.

1 ICA 算法

ICA 的基本模型[14]可以具体描述为：

式（1）中：x=（x1，x2，…，xn）T为n 维随机观测向量；s=（s1，s2，…，sm）T是m 维未知源信号，si称为独立分量；A 为满秩的混合矩阵（为了简化模型，假设A 为m=n 的方阵）；ai是混合矩阵的基向量.

为了寻找变换矩阵或称解混矩阵W，对x 进行线性变换，得到n 维输出向量

在比例允许和顺序不确定的前提下，y=sˆ是对si的一个估计量.以2 个源信号为例，ICA 线性组合模型的基本运算流程如图1 所示.

2 FastICA 算法

FastICA 算法又称Fixed-Point 算法，是一种通过牛顿迭代法使目标函数最优化的迭代算法，具有收敛速度快、稳健性高等优点.与梯度算法相比，FastICA 算法不要求学习步长，可通过确定数据块的学习方式调整分离矩阵系数.为减少算法需要估计的参数，简化算法的计算需要对数据进行预处理.FastICA 算法的预处理步骤是对观测信号进行白化处理，使白化后的分量x 为非相关的，且满足E{xxT}=I.

在非高斯性最大化方法中，FastICA 算法基于负熵来进行估计.为求出y=wTx 的投影方向，首先利用式（3）进行最大化

式（3）中：w 是m 维变量.在预处理后，式（3）的最大化问题可以转化为E[G（wTx）]的优化问题.在约束条件E[（wTx）2]的限制下，E[G（wTx）]的优化可以通过式（4）获得，

式（4）中：β 是常量，可通过β=E[w0Txg（w0Tx）]获得；w0是w 的初始值，函数g（u）为

式（5）中：1≤a1≤2；a2≈1.利用牛顿迭代法解式（4），式（4）左边部分记为F（w），其雅可比（Jacobian）矩阵

为了简化求矩阵的逆，利用式（7）对式（6）中的E[xxTg′（wTx）]进行近似，

自媒体又称“个人媒体”或“公民媒体”，是普通大众以现代化和电子化的手段，传递规范性及非规范性信息，提供与分享他们自身的事实、新闻的途径。本研究通过微信和新浪微博等自媒体平台，在搜索栏输入“黄山民宿”，剔除纯图片，纯营销，纯介绍性文字，得到最终有效游记和评价共198条。

亦可以得到式（7）的牛顿迭代法的近似解

为了提高算法的稳定性，在迭代后利用wk+1=wk+1/‖wk+1‖归一化w.在式（8）两边分别乘以β=E[g′（wTx）]，得到固定点算法的迭代公式

基于式（9），FastICA 算法即固定点迭代算法的基本结构为：

（1）权向量w 的初始化选择，k=0.

（2）对wk+1进行调整.

3 实验结果

地震记录中的有效信号可以看作是地震子波与地层反射系数序列的卷积，假设地层反射系数是非高斯分布的，即地震信号也非高斯分布，由于有效地震信号和50 Hz 工频干扰分别由不同的信源产生，彼此相互独立，因此地震信号和工频干扰信号满足ICA 的应用限制条件[16].

为了满足FastICA 的应用条件，构造观测信号x2（t）=a21s（t）+a22n（t）.实际情况下，s（t）未知，工频干扰信号满足n（t）=A0sin（2πf0t+θ0）.在不考虑幅值的情况下，构造观测信号x2（t）=0×s（t）+a22（sin 2πf0t+θ0）以符合ICA 的要求.观测信号可进一步表示为

式（10）中：b=a12A0cos θ0，c=a12A0sin θ0.可以通过构造2 组参考信号s2=sin（2πf0t）和s3=cos（2πf0t）与含有工频干扰和地震信号的混合信号共同作为ICA 的输入信号，分离出有效地震信号s（t），这样就避免了对θ0的直接估计[17]，其数学模型表示为

利用FastICA 算法对式（11）中的混合矩阵A 和分离矩阵W=A-1进行估计，并利用式（12）求出不含工频干扰的地震信号s（t）.

式（12）中：X 和S 分别为观测信号向量和源信号向量；A 为混合矩阵.

在地震信号中选取3 000 个采模点进行降噪实验，FastICA 算法对混有工频干扰的地震信号的干扰消除结果如图2 所示.

图2c 为混入50 Hz 工频干扰的地震源信号，图2d 和图2e 为依据前文介绍算法构造出的两路参考信号.图2f～图2h 为经过FastICA 算法处理后三路信号（图2c～图2e）的恢复信号，可以看出，图2g 所示的第2 路恢复信号为所需要的信号.

信号的相关系数描述了信号之间的线性相关性，表1 给出经过20 次实验所得的分离信号和源信号的相关系数绝对值.由表1 中数据可以看出，分离信号与源信号的相关系数绝对值均高于0.998 8，说明分离效果较为理想.

表1 分离信号与源信号的相关系数绝对值Tab.1 Absolute values of correlation of separative signal and source signal

4 对比实验

4.1 不同方法的比较

采用有限长单位冲激响应（Finite Impulse Response，FIR）陷波滤波器和最小均方误差（Least Mean Square，LMS）自适应滤波器2 种传统消除工频干扰的方法与FastICA 算法进行比较.实验中，陷波器的采样频率为fs=1 000 Hz，通过将边界频率转换为归一化频率可减少其带宽；LMS 自适应滤波器的步长因子为0.001，自适应滤波器可以自动调整滤波器系数，以跟踪输入过程的变化，并实现工频干扰的自适应抵消.

FIR 陷波滤波器对混有工频干扰的地震信号滤波的结果如图3 所示.由图3 可以看出，输出信号（图3c）较输入信号（图3a）存在明显延迟，严重影响了地震信号分离的时效性.

利用LMS 自适应滤波器消除地震信号中工频干扰的结果如图4 所示.加入正弦干扰（即工频干扰）的原始信号（图4b）在经过LMS 自适应滤波器的降噪处理后，得到的输出信号如图4c 所示.图4d为误差信号，对于地震信号而言，此结果的信号误差较大.

表2 给出FIR 陷波滤波器、LMS 自适应滤波器和FastICA 算法输出结果相关系数绝对值的比较结果.由表2可以明显看出，采用FastICA 算法得到的恢复信号与源信号的相关系数要高于其他2 种典型的传统滤波方法，这说明FastICA 算法消除地震信号中的工频干扰能力优于传统方法.

表2 3 种方法相关系数绝对值的对比Tab.2 Comparison of absolute values of correlation of 3 methods

4.2 采样点对结果的影响

考虑到ICA 本身的统计特性，减少各信号的采样点至1 000 个时，实验结果如图5 所示.

由图5 可知，ICA 算法的准确性取决于采样点的数量.采样点越少，即先验知识越少，恢复信号与原信号的相关系数就越小；采样点越多，得到的相关系数越大，但随之而来的算法的时间度亦将变大，从而给算法的时效性带来了影响.

表3 给出1 000 个采样点时FastICA 算法消除工频干扰后分离信号与源信号相关系数的绝对值.同样进行20 次实验，1 000 个采样点相关系数绝对值的均值为0.971 8，显然小于3 000 个采样点的均值0.998 8，说明采样点减少后，恢复出来的信号失真率较高.进一步说明采样点对FastICA 算法的影响，也说明ICA 算法这种基于信号统计特性的算法存在的不足之处.

表3 选取1 000 个采样点时，分离信号与源信号的相关系数绝对值Tab.3 Results of absolute values of correlation of separative signal and source signal in 1 000 sampling points

5 结论

根据ICA 的数学模型，分析了不同ICA 算法以及它们之间的联系.将FastICA 算法应用于消除地震信号中的工频干扰，并对输出信号的相关系数绝对值进行对比.结果表明：经FastICA 算法处理后，信号的相关系数绝对值优于传统方法，说明在消除微信号工频干扰的问题中，与传统工频干扰消除技术相比，FastICA 算法可以更准确地从混合信号中分离出实际所需要的源信号，证明FastICA 算法在消除地震信号中的工频干扰方面具有良好的应用前景.

[1]NADAL J，PARGA N.Nonlinear neurons in low-noise limit：A factorial code maximizes information transfer[J].Network，1994，5（3）：565—581.

[2]BELL A，SEJNOWSKI T.An information maximization approach to blind separation and blind de-convolution[J].Neural Computation，1995，7（6）：1129—1159.

[3]HYVÄRINEN A，OJA E.Independent component analysis：Algorithms and applications[J].Neural Networks，2000，13（4/5）：411—430.

[4]HYVÄRINEN A，PAJUNEN P.Nonlinear independent component analysis：Existence and uniqueness results[J].Neural Networks，1999，12（3）：429—439.

[5]JUTTEN C，ZADEH M，HOSSEINI S.Three easy ways for separating nonlinear mixtures[J].Signal Processing，2004，84（2）：217—229.

[6]LIU W，RAJAPAKSE J.Approach and applications of constrained ICA[J].Trans Neural Networks，2005，16（1）：203—212.

[7]PLUMBLEY M.Algorithms for nonnegative component independent analysis[J].Trans Neural Networks，2003，14（3）：534—543.

[8]王世海，陈向东，毕雪，等.多分辨率子带分解的独立分量分析算法在红外图像去噪上的应用[J].电子技术应用，2007（6）：66—68.

[9]LEE T，LEWICKI M，GIROLAMI M，et al.Blind source separation of more sources than mixtures using over-complete representations[J].Signal Processing Letters，1999，4（5）：205—208.

[10]白琳，陈豪.一种奇异混合阵的盲信号提取算法[J].空间电子技术，2012（1）：6—10，22.

[11]李远清，张丽清，刘永清.一种奇异混合信号盲分离的神经网络模型[J].控制理论与应用，2001，18（1）：105—108.

[12]肖俊，何为伟.一种基于信息最大化的时延混合盲分离算法[J].井冈山师范学院学报：自然科学版，2005，26（3）：12—14.

[13]HYVÄRINEN A，KARHUNEN J，OJA C.Independent Component Analysis[M].New York：Wiley，2001：293—304.

[14]杨福生，洪波，唐庆玉.独立分量分析及其在生物医学工程中的应用[J].国外医学：生物医学工程分册，2000，23（3）：129—134.

[15]郭松.独立分量分析及其应用研究[D].武汉：武汉大学，2004：13—20.

[16]刘喜武，刘洪，李幼铭.独立分量分析及其在地震信息处理中的应用初探[J].地球物理学进展，2003，18（1）：90—96.

[17]吴小培，詹长安，周荷琴，等.采用独立分量分析的方法消除信号中的工频干扰[J].中国科技大学学报，2000，30（6）：671—676.