基于机器学习（霍特林统计量，高斯混合模型，支持向量机）的工业数据异常检测（MATLAB R2021B）
近年来，隨着集散控制系统、工业物联网、智能仪表等信息技术在现代工业生产系统中的应用，生产过程的运行状态能够以大量数据的形式被感知和记录。基于数据的故障诊断方法以过程数据为基础，采用统计分析、统计学习、信号处理等方法，充分挖掘过程数据的隐含信息，建立数据驱动的模型实现故障诊断。事实上，基于数据的故障诊断研究起源于上世纪的控制图；后来，研究者们进一步提出了累积和控制图、指数加权滑动平均控制图等方法。然而，这些方法只关注单个变量的变化，忽略了变量之间的相互关系，并不适用于复杂的、大规模的现代工业生产过程。因此，研究者们开始利用面向多变量场景的统计分析方法进行故障诊断，主要包括主成分分析、独立成分分析、偏最小二乘、Fisher判别分析等。利用这些统计分析方法可对高维过程数据进行降维处理，通过提取过程数据的本质特征，实现故障诊断；进一步地，为了提高建模准确性，改善故障诊断效果，研究者们结合具体的过程数据特性（如非线性、非高斯、时变、多工况、动态、随机、离群点等），对上述统计分析方法进行了完善。统计学习方法主要包括支持向量机、人工神经网络、K近邻、决策树、随机森林、高斯混合模型、隐马尔科夫模型、贝叶斯网络等；信号处理方法主要包括小波变换、频谱分析等。基于数据的故障诊断方法不但能够规避对复杂过程机理的依赖，而且具有很强的实用性和通用性，已成为当今故障诊断领域的研究热点。值得说明的是，从过程数据挖掘的角度来看，统计分析、统计学习等方法均可归类为机器学习方法，二者的共同之处在于对高维过程数据进行深入挖掘，提取过程数据中隐含的高价值信息。这些高价值信息可作为反映过程状态的本质特征，利用这些高价值信息，可建立相应的故障检测、溯源、分类等模型，有效地完成各种故障诊断任务。
鉴于此，采用基于机器学习方法（霍特林统计量，高斯混合模型，支持向量机）对工业数据进行异常检测,运行环境为MATLAB R2021B，数据集为C-MAPSS涡扇发动机数据集。
目前大多数论文的实验中使用的数据集是美国国家航空航天局的 C-MAPSS 数据集。C-MAPSS 数据集是由模拟航空发动机的模拟软件生成。模拟发动机的结构图如下：
监控涡扇发动机状况为21个传感器。由于传感器的单位不同，传感器记录的数值的量级也有所差异，位于10的-2次方到10的3次方之间。例如，燃烧室油气比数值的量级是10的-2次方；低压涡轮冷气流量数值的量级是10的1次方。表2-2描述了NASA的C-MAPSS数据集。由于不同的操作条件和故障模式，数据集可以分为四个子数据集，依次是FD001、FD002、FD003和FD004。每个子数据分为训练集和测试集，记录了发动机的3种操作设置和21个传感器数据。每个子数据集通过.txt文件单独保存。在.txt文件中，每一行记录了一个引擎某个时间刻的3种操作设置和21个传感器数据。关于故障模式和操作条件方面，FD001和FD002子数据集包含一种故障模式（高压压气机退化），FD003和FD004包含两种故障模式（高压压气机退化和风扇退化）；FD001和FD003只有一种操作条件，FD002和FD004有六种操作条件。由于FD002和FD004子数据集引擎的操作环境复杂多变，FD002和FD004子数据集中RUL的预测更加困难。
%% Set Variable Names
varName = {'Unit', 'Time', 'Setting1', 'Setting2', 'Setting3', 'FanInletTemp',...
'LPCOutletTemp', 'HPCOutletTemp', 'LPTOutletTemp', 'FanInletPres', ...
'BypassDuctPres', 'TotalHPCOutletPres', 'PhysFanSpeed', 'PhysCoreSpeed', ...
'EnginePresRatio', 'StaticHPCOutletPres', 'FuelFlowRatio', 'CorrFanSpeed', ...
'CorrCoreSpeed', 'BypassRatio', 'BurnerFuelAirRatio', 'BleedEnthalpy', ...
'DemandFanSpeed', 'DemandCorrFanSpeed', 'HPTCoolantBleed', 'LPTCoolantBleed'};
dataIn.Properties.VariableNames = varName;
dataIn.Properties.VariableNames = varName;
%% Extract Effective Sensors
selectedVarName = varName([1:2, 7:9, 12:14, 16:20, 22, 25:26]);
NumOfUnits = 100;
dataIn = dataIn(:,selectedVarName);
for kk = 1:NumOfUnits
idx = dataIn.Unit == kk;
dataIn.Time(idx) = (dataIn.Time(idx) - max(dataIn.Time(idx)));
end
%% Labeling the condition into 4 categories based on the remaining cycles till failure occurs
% - 0~50 cycles : urgent
% - 51~125 cycles : short
% - 126~200 cycles: medium
% - 201~ : long
catname = {'urgent','short','medium','long'};
label = discretize(-dataIn.Time, [0 51 126 201 inf],'categorical',catname);
dataIn.Label = label;
dataOut = dataIn;完整数据和代码通过知乎学术咨询获得：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。
