浅谈盲源分离算法在故障诊断中的应用（附matlab代码）
图2 信号的预处理和特征提取流程图
诊断方法的研究
特征提取方法虽然能够有效反映机械设备故障信息的特征，但是仍需要进一步实行有效的故障分类，即故障识别与预测。随着大数据、云计算、互联网、物联网等信息技术的发展，人工智能被应用到各个领域。在诊断方法层面，人工智能诊断方法通过对提取的特征信息的学习、 抽象、推理和决策，分析和判断故障的类型、故障的位置和故障的程度。
盲源分离
盲源分离（Blind Source Separation, BSS），即在混叠通道和源信号未知的情况下仅由观测信号（也称传感器采集的信号或混合信号）分离出源信号。盲源分离最早起源于鸡尾酒会问题，即在一个大型聚会场所上，几个固定的麦克风采集到在场所有说话人员的声音，但是机器无法通过麦克风采集到的混叠信号将对应的发生源辨识出来，盲源分离技术就是通过混叠信号得到源信号，其过程如图3所示。
图3 盲源分离简图
根据源信号的个数（n）与传感器数量（m）的关系，盲源分离模型被分成超定盲源分离（m>n）、正定盲源分离（m=n）与欠定盲源分离（m <n）。在实际情况中，人们希望利用尽可能少的传感器采集到的信号去分离出更多的源信号，故欠定盲源分离模型更符合实际情况，被更多的学者选择和研究。
源信号的混合过程可以分为线性和非线性两种情况，对于盲源分离模型而言，两种情况分别对应线性瞬时混合模型和卷积混合模型两种。由于受到相关理论和技术的限制，相比于卷积混合模型来说，线性瞬时混合模型更简单，可以认为是卷积混合模型一种特殊情况，长期以来，学者们研究工作主要集中在线性混合瞬时模型上。因此，本文以最具代表性的线性瞬时混合模型进行说明，BSS线性瞬时混合模型如下：
公式中，X(t)代表m个混合信号；A是一个的混合矩阵，表示源信号的线性混合过程；S(t)为n个源信号；N(t)为环境噪声，t为采样时间点。
在实际处理过程中，一般只知道混合信号这一先验条件，并不知道其它相关信息，盲源分离技术就是利用已知的m个混合信号这个唯一条件在误差容忍范围内分离出 n 个源信号。目前，盲源分离技术主要有三个研究方向：独立成分分析、稀疏成分分析以及盲解卷积。不同的盲源分离模型对应着不同的解决手段，在线性混合模型中，对于超定与正定模型，一般采用独立分量分析法对这两种模型进行求解；对于欠定分离模型，一般根据信号的稀疏性，利用稀疏成分分析法解决分离问题。为了更好的理解盲源分离，下面以一个欠定的线性混合仿真实验进行说明。
图4 源信号时域波形和频谱图
利用Matlab软件随机生成一个混合矩阵A和噪声矩阵N，用来模拟实际过程得到的两个混合信号，A的表达式如下:
将仿真数据带入前文BSS线性瞬时混合模型，可以得到混合信号X，其波形图和频谱图如图5示。
图5 混合信号时域波形和频谱图
通过分析图5中每个混合信号的频谱图无法较为完整的得到源信号的频率成分，即均存在频率成分被淹没、具体频率成分数量无法确定的问题，会严重影响对故障的准确判断。对此，采用盲源分离对两个混合信号进行分离，分离结果如图6所示。
图6 分离的源信号时域波形和频谱图
对比分析分离的源信号图6与仿真的源信号图4可知，采用盲源分离技术得到的3个源信号的频谱图与仿真的源信号基本一致，即仿真的源信号得到了较好的恢复，更能直观的分析得到混合故障中的每个单一故障。若对图6中的信号进行适当降噪处理，可得到如下波形和频谱图，如图7所示。
图7 降噪后的时域波形和频谱图
对比分析图7与仿真的源信号图4可知，降噪后的3个源信号的外形和频谱图均与仿真的源信号的一致，即每个源信号的特征频率都被准确得到分离。相比于单个分析每个混合信号，盲源分离能更准确、全面判断出三种故障。
上述仿真过程的Matlab代码如下所示：
clcclearclose all%% 仿真信号load A2（A2为文字内容中的混合矩阵A）N=1024*2; % 信号长度fs=1024; % 采样频率f0=50; % 基频t=0:1/fs:(N-1)/fs;s1=sin(2*pi*f0*t)+sin(2*pi*2*f0*t);s2=sin(2*pi*f0*t)+sin(2*pi*0.33*f0*t);s3=sin(2*pi*0.45*f0*t);S=[s1;s2;s3]; %源信号%仿真源信号时域、频域图像figure(1)subplot(321)plot(S(1,:))subplot(323)plot(S(2,:))subplot(325)plot(S(3,:))subplot(322)hua_fft(s1,fs,1,0,500)subplot(324)hua_fft(s2,fs,1,0,500)subplot(326)hua_fft(s3,fs,1,0,500)%% 混合信号及图像X=A*S+0.5*randn(2,N);%混合信号时域、频域图像figure(3)subplot(221)plot(X(1,:))subplot(223)plot(X(2,:))subplot(222)hua_fft(X(1,:),fs,1,0,500)subplot(224)hua_fft(X(2,:),fs,1,0,500)%% 盲源分离分解信号X_fft=X*dftmtx(N);%% 源信号信号在傅里叶基下的恢复(最短路径法)eS_fft1=SourceRecovery_ShortestPathMethon(X_fft,A);eS_1=real(eS_fft1/dftmtx(N));%% 恢复源信号时域图像figure(5)subplot(321)plot(eS_1(1,:))subplot(323)plot(eS_1(2,:))subplot(325)plot(eS_1(3,:))subplot(322)hua_fft(eS_1(1,:),fs,1,0,500)subplot(324)hua_fft(eS_1(2,:),fs,1,0,500)subplot(326)hua_fft(eS_1(3,:),fs,1,0,500)S_corr1=corr(S',eS_1');%% 适当降噪处理eS_fft2=abs(eS_fft1);c=max(eS_fft2,[],2);for j=1:3eS_fft1(j,find(eS_fft2(j,:)<(c(j,1))/3))=0;endeS_1=real(eS_fft1/dftmtx(N));%% 降噪后源信号时域图像figure(6)subplot(321)plot(eS_1(1,:))subplot(323)plot(eS_1(2,:))subplot(325)plot(eS_1(3,:))subplot(322)hua_fft(eS_1(1,:),fs,1,0,500)subplot(324)hua_fft(eS_1(2,:),fs,1,0,500)subplot(326)hua_fft(eS_1(3,:),fs,1,0,500)S_corr2=corr(S',eS_1')；%%%%用到的函数function eS = SourceRecovery_ShortestPathMethon(X,A)% 把矩阵A映射到C_N^(m)个子空间中来恢复源，使得重构误差最小、% X -- 观测信号（时频点复数）% A -- 混合矩阵% eS -- 估计出的源信号n = size(A,2);m = size(X,1);N = size(X,2);Combine = nchoosek(1:n, m); % 从1~n中，每次取m，所有组合为Combine混合矩阵A构成m*m型矩阵的几种可能组合eS = zeros(n, size(X,2));residual = zeros(N, size(Combine,1)); % 误差系数for j = 1:size(Combine,1)residual(:,j)=sum(abs(pinv(A(:,Combine(j,:)))* X).^2,1);end[~,u] = min(residual,[],2);for i = 1:size(Combine,1)idx = find(u==i); % 找出每种组合中的最小误差对应的eS(Combine(i,:),idx) = pinv(A(:,Combine(i,:)))* X(:,idx); % X(:,idx)表示误差最小的时频点endfunction hua_fft(y,fs,style,varargin)% 作用：画信号的幅值谱和功率谱% 输入：% y -- 输入信号（一个向量，行或列向量）% fs -- 采样频率% style - 当style=1,画幅值谱；当style=2,画功率谱;当style=其他的数字，那么画幅值谱和功率谱% - 当style=1时，还可以多输入2个可选参数，如：hua_fft(y,fs,1,0,500)表示查看[0,500]段的频率图像（第一个为设置的频率起点，第二个为终点）% 输出：% 直接绘制出对应的图像%%y=y-mean(y); % 去除信号中的直流分量，防止对频谱分析造成干扰nfft= 2^nextpow2(length(y)); % 自动计算最佳FFT步长nfft）y_ft=fft(y,nfft);% 计算功率谱（相当于y轴）y_p=y_ft.*conj(y_ft)/nfft; % conj()函数是求y函数的共轭复数，实数的共轭复数是他本身。如：conj(1+2i)=1-2i% 幅值谱和功率谱的横坐标（频率轴）y_f=fs*(0:nfft/2-1)/nfft; % DFT变换后对应的频率的序列if style==1 % 画幅值谱if nargin==3plot(y_f,2*abs(y_ft(1:nfft/2))/length(y));ylabel('幅值');xlabel('频率');title('信号幅值谱');else % 画出指定频率区间的幅值谱f1=varargin{1};fn=varargin{2};ni=round(f1 * nfft/fs+1);na=round(fn * nfft/fs+1);plot(y_f(ni:na),abs(y_ft(ni:na)*2/length(y)));endelseif style==2 % 画功率谱plot(y_f,y_p(1:nfft/2));ylabel('功率谱密度');xlabel('频率');title('信号功率谱');else % 画幅值谱和功率谱subplot(211);plot(y_f,2*abs(y_ft(1:nfft/2))/length(y));ylabel('幅值');xlabel('频率');title('信号幅值谱');subplot(212);plot(y_f,y_p(1:nfft/2));ylabel('功率谱密度');xlabel('频率');title('信号功率谱');endend
编辑：王畅
校核：李正平、张泽明、张勇、赵栓栓、陈凯歌
来源：故障诊断与python学习
