过程监控中会用到很多中方法,如主成分分析(PCA)、慢特征分析(SFA)、概率MVA方法或独立成分分析(ICA)等为主流算法。
其中PCA主要多用于降维及特征提取,且只对正太分布(高斯分布)数据样本有效;SFA被用来学习过程监控的时间相关表示,SFA不仅可以通过监测稳态分布来检测与运行条件的偏差,还可以根据时间分布来识别过程的动态异常,多用于分类分析;概率MVA方法,多以解决动力学、时变、非线性等问题。
今天要介绍的是独立成分分析(ICA),由浅入深,细细道来。
此外文末还附有ICA可实现的代码哟~不要错过
独立成分分析(Independent Component Analysis,ICA)
基本原理
在信号处理中,独立成分分析(ICA)是一种用于将多元信号分离为加性子分量的计算方法。这是通过假设子分量是非高斯信号,并且在统计上彼此独立来完成的。ICA是盲源分离的特例。一个常见的示例应用程序是在嘈杂的房间中聆听一个人的语音的“ 鸡尾酒会问题 ”。
首先,引入一下经典的鸡尾酒宴会问题(Cocktail Party Problem)。
假设我们令一个未知的混合系数矩阵(mixing coefficient matrix)为A ,用来组合叠加信号S ,
ICA的不确定性(ICA ambiguities)
由于h 和s 都不确定,那么在没有先验知识的情况下,无法同时确定这两个相关参数。
ICA 算法
下面直接上ICA算法。
独立成分分析 ICA(Independent Component Correlation Algorithm)是一种函数,X为n维观测信号矢量,S为独立的m(m<=n)维未知源信号矢量,矩阵A被称为混合矩阵。ICA的目的就是寻找解混矩阵W(A的逆矩阵),然后对X进行线性变换,得到输出向量U。
此公式代表一个假设前提:每个人发出的声音信号各自独立。
求导可得,
举个Paper的栗子
下面为我们观测到的信号:
然后,再通过ICA还原后的信号为:
MATLAB代码实现
MATLAB代码:
Fast ICA
% Input:X 行变量维数,列采样个数;需要对原始矩阵转置 % Output:Sources重构的原信号, Q白化矩阵, P白化信号解混矩阵 function [Sources, Q, P] = FastICA(X, P) % 白化处理 [dim, numSample] = size(X); Xcov = cov(X'); [U, lambda] = eig(Xcov); Q = lambda^(-1/2)*U'; Z = Q*X; % FastICA maxiteration = 10000; %最大迭代次数 error = 1e-5; % 收敛误差 % P = randn(dim,dim); % 随机初始化P,并按照列更新 for k = 1:dim Pk = P(:,k); Pk = Pk./norm(Pk); % 向量归一化 lastPk = zeros(dim,1); % 0不需要再归一化 count = 0; while abs(Pk - lastPk)&abs(Pk + lastPk) > error count = count + 1; lastPk = Pk; g = tanh(lastPk'*Z); % g(y)函数 dg = 1 - g.^2; % g(y)的一阶导函数 %-------------------------------核心公式------------------------------------ Pk = mean(Z.*repmat(g,dim,1), 2) - repmat(mean(dg),dim,1).*lastPk; Pk = Pk - sum(repmat(Pk'*P(:,1:k-1),dim,1).*P(:,1:k-1),2); Pk = Pk./norm(Pk); %-------------------------------------------------------------------------- if count == maxiteration fprintf('第%d个分量在%d次迭代内不收敛!\n',k,maxiteration); break; end end P(:,k) = Pk; end Sources = P'*Z; % end
此外还有基于故障诊断的ICA算法代码实现->在这里
下面给出部分代码:
% 基于ICA的故障诊断 clear;clc;close all; load('MPD2000.mat'); Xnormal = MPD0'; % 数据归一化 [dim, numSample] = size(Xnormal); XnormalMean = mean(Xnormal, 2); XnormalStd = std(Xnormal, 0, 2); XnormalNorm = normalization(Xnormal, XnormalMean, XnormalStd); % 正常数据计算 解混矩阵W % P = rand(dim,dim)*100; load('P.mat'); [S, Q, P] = FastICA(XnormalNorm, P); W = P'*Q; % ------------------------利用2范数大小对W的行重新排列--------------------- Wnorm = zeros(dim,1); for k = 1:dim Wnorm(k) = norm(W(k,:)); end [Wnorm, indices] = sort(Wnorm, 'descend'); % -------------------------确定主导成分Sd与参与成分Se---------------------- threshold = 0.80; percentage = cumsum(Wnorm)./sum(Wnorm); for k = 1:dim if percentage(k) > threshold break; end end
效果如下: