異常檢測通常是根據已有的觀測數據建立正常行為模型，從而將不同機制下產生的遠離正常行為的數據劃分為異常類，進而實現對異常狀態的檢測。常用的異常檢測方法主要有：統計方法、信息度量方法、譜映射方法、聚類方法、近鄰方法和分類方法等。

早期的異常檢測方法主要是基于統計理論，通過對觀測數據進行統計建模或密度估計來評價異常，如果數據的分布假設成立，統計方法能夠得到一個滿意的可證明的解，且為無監督過程，但該方法對數據的統計假設較為依賴，特別是對于高維復雜數據，統計分布往往不能準確假設。而信息度量方法雖不需對數據分布進行假設，但其檢測性能取決于信息度量指標，若信息度量不能明顯反映異常數據與正常數據的差異，則該方法將失效。基于主成分分析或矩陣分解的空間映射方法，通過將正常和異常數據嵌入到一個低維子空間進行劃分，并可自動執行維度約簡，適合于高維復雜數據，但在空間映射過程中易丟失部分有用的數據信息。

隨著人工智能理論的產生和發展，許多機器學習和數據挖掘的方法被應用于異常檢測，主要有基于聚類或近鄰的無監督學習方法和基于分類的有監督學習方法。對于聚類方法，不同的計算正常類數據間距離或相似度的策略對聚類結果有很大影響，如果形成的聚類簇太大，易將異常數據包含進去而導致誤判。對于近鄰方法，通過計算正常類數據的k近鄰或局部異常因子來進行異常檢測，體現了數據的局部分布特性，檢測性能得到了提高，但近鄰方法需要存儲所有的正常數據，來計算與測試數據間的距離，復雜度較高，而且如果正常數據的分布比較稀疏（即沒有足夠的近鄰點），易導致較大的檢測誤差。對于分類方法，通過使用可獲得的標記數據作為訓練樣本來構造分類模型，對測試數據進行分類，能夠實現一類或多類分類異常檢測。

鑒于此，采用支持向量機、孤立森林和LSTM自編碼器方法對機械狀態進行異常檢測，運行環境為MATLAB R2021B。數據集包含來自工業機器的三軸振動測量值， 在計劃維護之前和維護之后采集數據。假定在定期維護后采集的數據代表機器的正常運行狀況，維護前的數據代表正常或異常情況。每軸的數據存儲在單獨的列中，每個文件包含 7000 個測量值。

?

function[net]=AE(X,Xtest,Labels,options)
% S-AEEL:Stacked AutoEncoder with Embedded Labels 
% it allows labels embedding inside the hidden layer
% Inputs
% number_neurons:number of neurons in hidden layer.
% X: the training inputs.
% Xtest: the testing inputs
% prefomance: RMSE of training.
% options :training options (3 options -->)
% 
% get options
number_neurons=options.number_neurons; % number of neurons
Layers=options.Layers;% number of layers
ActivF=options.ActivF;% activation function
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
X = scaledata(X,0,1);% data Normalization
x=X;% save a copy from original input;
alpha=size(X);% get information for weights generations
% 1:generate a random input weights
input_weights{1}=rand(number_neurons,alpha(2))*2-1;for i=1:Layers
% 2:calculating the hidden layer
tempH=input_weights{i}*X';
CAT=repmat(Labels',size(tempH,1),1);
tempH=[tempH+CAT];
tempH=[tempH;Labels'];
% activation function
switch lower(ActivF)case {'sig','sigmoid'}%%%%%%%% Sigmoid H = 1 ./ (1 + exp(-tempH));case {'sin','sine'}%%%%%%%% SineH = sin(tempH);    case {'hardlim'}%%%%%%%% Hard LimitH = double(hardlim(tempH));case {'tribas'}%%%%%%%% Triangular basis functionH = tribas(tempH);case {'radbas'}%%%%%%%% Radial basis functionH = radbas(tempH);%%%%%%%% More activation functions can be added here                
end
% 3: calculate the output weights beta
B{i}=pinv(H') * X ; %Moore-Penrose pseudoinverse of matrix
% calculate the output : Unlike other networks the AEs uses the same weight
% beta as an input weigth for coding and output weights for decoding
% we will no longer use the old input weights:input_weights. 
Hnew=X*B{i}';
output=Hnew*pinv(B{i}');
input_weights{i+1}=B{i};
end
% 4:calculate the prefomance
labels_hat=scaledata(Hnew(:,end-(size(Labels,2))+1:end),min(Labels),max(Labels));% recovered labels
prefomance=sqrt(mse(x-output));% RMSE of training samples recontruction
rec_acc=sqrt(mse(labels_hat-Labels));% % RMSE of training labels recovering
%% mappings for test set
H_test=Xtest*B{i}';
%% save results
net.xtr_mapp=Hnew; % feature mappings of training inputs
net.xts_mapp=H_test; % feature mappings of testing inputs
net.recovered_labels=labels_hat;% 
net.learning_loss=prefomance;
net.Learning_weights=B;
net.xtr_hat=output;
net.rec_acc=rec_acc;
end

工學博士，擔任《Mechanical System and Signal Processing》《中國電機工程學報》《控制與決策》等期刊審稿專家，擅長領域：現代信號處理，機器學習，深度學習，數字孿生，時間序列分析，設備缺陷檢測、設備異常檢測、設備智能故障診斷與健康管理PHM等。完整數據和代碼通過知乎學術咨詢獲得：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1