统计过程控制图(statistical process control chart, SPCC)被广泛用来监控制造过程的平稳性。如果制造过程是平稳的，说明过程中只存在随机误差的影响，生产是正常(normal/natural)^[1]的；如果制造过程是非平稳的，说明除了随机误差外，还存在系统误差，即生产是不正常(unnormal /unnatural)的，存在影响生产的系统因素，比如工具或夹具出现缓慢磨损、生产环境温度和湿度变化、生产设备电压波动等。当生产过程不正常时，需要立即停机检查人、机、料、法、环、测等系统因素，以消除这些系统因素，使生产恢复正常。一般这些引起生产异常的系统因素对应着不同的控制图模式。按照统计质量控制手册^[2]统计，一共有15种类型的控制图模式(control chart patterns, CCP)。这15种控制图模式中有8种是最基本的，分别是正常模式(normal, NOR)、分层模式(stratification, STA)、系统模式(systematic, SYS)、循环模式(cyclic, CYC)、趋势上升(increasing trend, UT)、趋势下降(decreasing trend, DT)、阶跃上升(upward shift, US)和阶跃下降(downward shift, DS)，如图1所示。其他控制图模式一般是这8种基本模式的特殊形式或者是几种基本模式的组合。

控制图的模式对应生产过程的各种状态，如表1所示。因此，识别控制图的不同模式，可以诊断生产过程处于何种状态，并采取后续的修正措施。随着系统复杂程度、自动化程度越来越高，系统产生的数据量越来越大，产生速度越来越快，实时性要求越来越高，这种情况下只依靠人工来对数据进行处理和分析并进行控制图模式的识别非常困难，控制图模式识别非常必要。

控制图模式识别方法很多，主要方法有：1) 专家系统(expert system, ES)^[3-4]；2) 模糊聚类方法(fuzzy clustering, FC)^[5]；3) 支持向量机方法(support vector machine)^[6-9]；4) 人工神经网络方法(artificial neural network, ANN)^[10-12]；5)决策树方法(classification and regression tree, CART)^[13-14]等。专家系统一般需要人工制定规则或者依靠其他几类分类器产生规则；模糊聚类方法分类精度较低；人工神经网络存在过拟合、收敛慢、易陷入局部最小、结构不易确定、对训练样本数量要求高等问题；支持向量机则需要参数调优，一般依靠人工经验或者通过遗传算法、蚁群算法、网格法等确定最优参数；决策树方法不需要领域知识和参数假设，但是容易过拟合，且没有考虑属性之间的相关性。

随机森林是一种灵活且易于使用的机器学习算法，不需要参数调优，也可以在大多数情况下取得很好的预测效果，因此非常适合控制图模式的分类识别。

1 特征选择

进行控制图模式自动识别时，首先需要进行特征提取和选择。关于特征提取的相关研究较多，比如Pham等^[15]、刘玉敏等^[7]提取了统计特征。但是统计特征提取后没有有效利用原有监测窗口数据的顺序信息。Gauri等^[16]对基于特征的CCP识别方法进行了广泛的研究并提取了能够反映数据顺序信息的形状特征，基于CART算法选择了7个形状特征来进行分类识别。本文通过查阅相关资料，尽可能涵盖控制图的所有统计特征和形状特征，整理资料如下。

1.1 统计特征

主要的统计特征如下

均值 $\overline x =\displaystyle \sum\limits_{i = 1}^N {{x_i}} /N$ 。

离散度(标准差)： $ {s^2} = \displaystyle \sum\limits_{i = 1}^N {{{({x_i} - \overline x )}^2}} /(N - 1), \; $ 　　 $s = \sqrt {\displaystyle \sum\limits_{i = 1}^N {{{({x_i} - \overline x )}^2}} /(N - 1)} $ 。

最大值 $ {X_{\max }} \!=\! \max \;({x_1},{x_2},\cdots{x_i}){\text{；}}i \!=\! 1,2,\cdot \cdots N$ 。

最小值 $ {X_{\min }} \!=\! \min \;({x_1},{x_2},\cdots {x_i}){\text{；}}i \!=\! 1,2,\cdot \cdots N$ 。

极差 $R = {X_{\max }} - {X_{\min }} $ 。

中位数 ${X_{\rm median}} = \left\{ {\begin{align}{{x_{(n + 1)/2}},}\qquad&{\left( n {\text{为奇数}}\right);}\\{\dfrac{{{x_n} + {x_{(n + 1)/2}}}}{2},}&\;{\left( n{\text{为偶数}} \right) {\text{。}}}\end{align}} \right. $

离散系数(变异系数) ${V_s} = s/\overline x $ 。

偏态系数 ${\rm SK} = \dfrac{{N\displaystyle \sum\limits_{i = 1}^N {{{({x_i} - \overline x )}^3}} }}{{(N - 1)(N - 2){s^3}}} $ 。

峰态系数

$\qquad K = \dfrac{{N(N + 1)\displaystyle \sum\limits_{i = 1}^N {{{({x_i} - \overline x )}^4} - 3{{\left[\displaystyle \sum\limits_{i = 1}^N {{{({x_i} - \overline x )}^2}} \right]}^2}(N - 1)} }}{{(N - 1)(N - 2)(N - 3){s^4}}} {\text{。}}$

平均偏差 $ {X_{\rm d}} = \dfrac{1}{N}\displaystyle \sum\limits_{i = 1}^N {|{x_i} - \overline x |}$ 。

平均幅值 ${X_{\rm amp}} =\displaystyle \sum\limits_{i = 1}^N {|{x_i}|} /N $ 。

均方幅 $ {X_{\rm ms}}{\rm{ = }}{{X}_{{\max}}}/{\rm{SK}}$ 。

波形 $ {X_{\rm k}} = {X_{\rm msamp}}/{X_{\rm amp}}$ 。

脉冲 $ {X_{\rm I}}{\rm{ = }}{X_{\max }}/{X_{\rm amp}}$ 。

1.2 形状特征

1) 建立在整个区域的最小二乘线基础上。

SB为控制图最小二乘线斜率标识值。

$\begin{split} &\qquad B = \frac{{\displaystyle \sum\limits_{i = 1}^N {{y_i}(t_i - \overline t )} }}{{\displaystyle \sum\limits_{i = 1}^N {{{(t_i - \overline t )}^2}} }} {\text{，}}\\ &\qquad{\rm SB }= \left\{ {\begin{align} 0, \; \; \;&{\rm if} \; {B {\text{＜}} 0} ;\\ 1, \; \; \;&{\rm if} \; {B {\text{≥}} 0} {\text{。}} \\ \end{align}} \right. \end{split} $

其中，B是控制图的最小二乘线斜率， $t_i^{}$ 是第 $i$ 次时间检测点到原点的距离， ${y_i}$ 是质量特性的第 $i$ 次检测值。SB可以区分下降模式(DT/DS)和上升模式(UT/US)。

RVE为观测值方差( ${\rm SD}^2$ )与全模式的平均最小二乘误差平方和(MSE)之比，

$\qquad{\rm RVE }= {\rm SD}^2/{\rm MSE} {\text{，}}$

$ \qquad{\rm SD}^2 = \frac{1}{{N - 1}}\sum\limits_{i{\rm{ = }}1}^N {{{({y_i} - \overline y )}^2}} {\text{，}} $

$ \qquad{\rm MSE} = \frac{1}{{N - 2}}\left[{\sum\limits_{i = 1}^N {({y_i} - \overline y )} ^2} - \frac{{{{\left(\displaystyle\sum\limits_{i = 1}^N {{y_i}({t_i} - \overline t )} \right)^2}}}}{{{{\displaystyle\sum\limits_{i = 1}^N {({t_i} - \overline t )^2} }}}}\right] {\text{，}} $

其中 $N$ 是观测值的个数。 ${t_i} = ic(i = 1,2\cdots N)$ 是第 $i$ 个观测时间与原点的距离， $c$ 是观测间隔； ${y_i}$ 是第 $i$ 个观测点的质量特性值； $\overline {{y}} {\rm{ = }}\displaystyle \sum\limits_{i{\rm{ = }}1}^N {{y_i}} /N;\overline {{t}} {\rm{ = }}\displaystyle \sum\limits_{i{\rm{ = }}1}^N {{t_i}} /N$ 。

RVE的值对NOR、STA、SYS、CYC和MIX来说大约等于1，而对于UT、US、DT、DS来说大于1。

2) 建立在模式与某条直线的面积基础上。

ACLPI为模式与中心线组成的面积与观测值标准差之比，

$\begin{split} & \qquad{\rm ACLPI }= [{\rm ACL}/(N - 1)]/{\rm SD} {\text{，}}\\ &\qquad{\rm SD }= \sqrt {\frac{1}{{N - 1}}\sum\limits_{i = 1}^N {({y_i} - \overline y } {)^2}} {\text{。}} \end{split} $

其中，ACL是模式与中心线组成的面积。

ACLMLC为每中心线交点的模式与中心线组成面积与观测标准差之比，

$\qquad{\rm ACLMLC }= \frac{{{\rm ACL/MLC}}}{\rm{SD}} {\text{。}}$

ALSPI为最小二乘线与观测值形成的平均间隔面积与观测值标准差之比，

$\qquad{\rm ALSPI }= [{\rm ALS}/(N - 1)]/{\rm SD}^2 {\text{，}}$

$\qquad{\rm SD}^2 = \frac{1}{{N - 1}}\sum\limits_{i{\rm{ = }}1}^N {{{({y_i} - \overline y )}^2}} {\text{。}}$

其中ALS是最小二乘线与观测值形成的总面积，可以由相邻两点与最小二乘线形成的梯形(如第5点和第6点)或三角形(如第10点和第11点)面积累积，如图2所示。

STA模式的ALSPI的值一般来说最大，SYS模式的值最低，其他模式中间。

ALSLSC为每最小二乘线交点的模式与最小二乘线组成面积与观测标准差之比，

$\qquad{\rm ALSLSC} = \frac{\rm{ALS/LSC}}{\rm{SD}} {\text{。}}$

RACLALS为ACL与ALS之比，

$\qquad{\rm RACLALS = {{ACL} / {ALS}}} {\text{。}}$

3) 建立在中心线及最小二乘线交点基础上。

PMLC为MLC与观测数量之比，

$\qquad{\rm PMLC = MLC}/N {\text{。}}$

PLSC为LSC与观测数量之比，

$\qquad{\rm PLSC = LSC}/N {\text{。}}$

PSMLSC为中心线交点与最小二乘线交点均值与观测数量之比，

$\qquad{\rm{PSMLSC}} = ({\rm{MLC}} + {\rm{LSC)}}/(2N) {\text{，}}$

其中，MLC是控制图与中心线的交点数量，LSC是最小二乘线与控制图的交点数量，N是控制图中点的数量。

SYS模式的PSMLSC值最大，NOR、STA、UT、DT和MIX模式的值中间，CYC、US和DS模式的PSMLSC值较小。

4) 建立在相邻连续点距离上。

RDIST为相邻两点平均距离与SD之比，

$\qquad{\rm RDIST }= \frac{1}{{N - 1}}\sum\limits_{i = 1}^N {\sqrt {{{({t_i} - \overline t )}^2} + {{({y_i} - \overline y )}^2}} } \Big/{{\rm SD} {\text{。}}}\Big.$

AASBP为所有通过连续相邻两点直线斜率的均值，

$\qquad{\rm AASBP }= \displaystyle \sum\limits_{i = 1}^{N - 1} {|({y_{i + 1}} - {y_i})/({t_{i + 1}} - {t_i})|/(N - 1)} {\text{。}}$

5) 建立在通过四均分区域中心点的直线上。

ASL为经过分块区域中间点组合的斜率平均值。其中分块区域指的是如图3所示将整个观测窗口分为4个区域，每个区域中间点的坐标计算方式为

$\left\{ \left[\displaystyle \sum\limits_{i = {n_1}}^{{n_1} + 7} {{t_i}} /8\right],\left[\displaystyle \sum\limits_{i = {n_1}}^{{n_1} + 7} {{y_i}} /8\right]\right\} $ ，　　其中， ${n_1}$ =1,9,17,25。任意两个区域的组合共有 ${C_4}^2 = 6$ 种，因此经过任意两个区域 $j$ 和 $k$ 中间点的直线有6条，其斜率为 ${S_{\!\!\!jk}}$ 。则

$\qquad\text{ASL}=\sum\limits_{\begin{smallmatrix} j,k , j{{<}}k \end{smallmatrix}}{{{S}\!\!\!_{jk}}/6,j=1,2,3;k=2,3,4}{\text{。}}$

SRANGE为如上所述6条直线斜率的变化范围，即最大值与最小值之差。

$ \begin{split} & \qquad{\rm SRANGE = max}\,({S\!\!\!_{jk}}) - {\rm min}\,({S\!\!\!_{jk}}), \\ &j = 1,2,3;k = 2,3,4;j {\text{＜}} k {\text{。}} \end{split} $

阶跃模式(US/DS)的SRANGE值比趋势模式(UT/DT)大；CYC模式的SRANGE值比NOR、STA、SYS和MIX模式的大，除非CYC模式的4个均分区域中每个区域都正好是一个循环，这种情况概率比较小。

6) 建立在四均分区域两两组合形成的6分段的最小二乘线基础上。

ABL为6分段的最小二乘线斜率均值，

$\qquad{\rm ABL} = \frac{1}{6}\sum\limits_ {j {\text{＜}} k} {{b\!\!_{jk}}}, j = 1,2,3;k = 2,3,4{\text{。}}$

BRANGE为6个区域最小二乘线斜率的极差，

$\qquad{\rm BRANGE} = \max \;({b\!\!_{jk}}) - \min \;({b\!\!_{jk}}) {\text{。}}$

REAE为上述的四均分区域可以两两组合形成6个组合，整个控制图区域的MSE与这6个区域的MSE均值之比，

$\qquad{\rm REAE }= \frac{\rm{MSE}}{{\dfrac{1}{6}\displaystyle\sum\limits_{i = 1}^6 {{\rm MSE}_i} }} {\text{，}}$

$\qquad{\rm MSE} = \frac{1}{{N - 2}}\left({\sum\limits_{i = 1}^N {({y_i} - \overline y )} ^2} - \frac{{{{\left(\displaystyle\sum\limits_{i = 1}^N {{y_i}({t_i} - \overline t )} \right)}^2}}}{{{{\displaystyle\sum\limits_{i = 1}^N {({t_i} - \overline t )} ^2}}}}\right) {\text{。}}$

CYC和阶跃模式(US/DS)的REAE值比1大，NOR、STA、SYS和趋势模式(UT/DT)的REAE值约为1，对于MIX模式来说，REAE值小于1，因此可以将MIX模式区分开。

RVAE为整个区域观测值的方差(SD²)与6个区域的MSE均值之比，

$\qquad{\rm RVAE }= \frac{{{{\rm SD}^2}}}{{\dfrac{1}{6}\displaystyle\sum\limits_{i = 1}^6 {{\rm MSE}_i} }} {\text{。}}$

DABL为不同组合区域的最小二乘线斜率，

$\qquad{\rm DABL} = \frac{1}{4}\sum\limits_{j = 1}^4 {{\rm db}_j}, j = 1,2,3,4 {\text{。}}$

其中， ${\rm db}_1$ 是第1和第2个区域的最小二乘线斜率； ${\rm db}_2$ 是第1、第2和第3个区域的最小二乘线斜率； ${\rm db}_3$ 是第3和第4个区域的最小二乘线斜率； ${\rm db}_4$ 是第2、第3个区域和第4个区域的最小二乘线斜率。

DBRANGE为上述4个区域斜率的极差，

$\qquad{\rm DBRANGE }= \max\; ({\rm db}) - \min \;({\rm db}) {\text{。}}$

7) 基于分界点的区域划分。

PMSE为两分段最小二乘误差平方和(SE₁和SE₂)的最小平均值(PMSE)。

两分段指的是将整个观测窗口区域分为两部分。对于阶跃模式来说，如果分成的两部分恰好是阶跃前后两部分，则两段分别做最小二乘回归线，这两条直线应该平行于X轴，分别计算两个分段的最小二乘误差平方(MSE1和MSE2)，其均值PMSE应该最小。因为阶跃位置事先未知，因此可以通过确定PMSE最小来确定阶跃位置。计算方法如下。

$\qquad{\rm PMSE} = ({\rm SE}_1 + {\rm SE}_2)/(N - 4) {\text{，}}$

$ \qquad{\rm SE}_1 = {\sum\limits_{i = 1}^m {({y_i} - \overline y )} ^2} - \frac{{{{\left(\displaystyle\sum\limits_{i = 1}^m {{y_i}({t_i} - \overline {{t_1}} )} \right)^2}}}}{{{{\displaystyle\sum\limits_{i = 1}^m {({t_i} - \overline {{t_1}} )} ^2}}}} {\text{，}} $

$ \qquad{\rm SE}_2 = {\sum\limits_{i = m + 1}^N {({y_i} - \overline y )} ^2} - \frac{{{{\left(\displaystyle\sum\limits_{i = m + 1}^N {{y_i}({t_i} - \overline {{t_2}} )} \right)^2}}}}{{{{\displaystyle\sum\limits_{i = m + 1}^N {({t_i} - \overline {{t_2}} )} ^2}}}} {\text{，}} $

其中，m是阶跃位置，一般仿真时可以设定 $m \in [8,24] $ ； $ \overline {{t_1}} = \displaystyle\sum\limits_{i = 1}^m {{t_i}} /m, \;\overline {{t_2}} =\displaystyle \sum\limits_{i = m + 1}^N {{t_i}} /m $ 。

AABPE为两分段区域最小二乘回归线斜率均值( $\overline {{B}} $ )的绝对值，

$\qquad{\rm AABPE }= \left|\frac{{{B_1} + {B_2}}}{2}\right| {\text{。}}$

ABDPE为整体最小二乘回归线斜率(B)与两分段区域最小二乘回归线斜率均值( $\overline {{B}} $ )之差的绝对值。

${\rm ABDPE }= |B - \overline B | $ ，

其中， $ \qquad \qquad B = \displaystyle\sum\limits_{i = 1}^N {{y_i}({t_i} - \overline t )} \Big/\displaystyle\sum\limits_{i = 1}^N {{{({t_i} - \overline t )}^2}} {\text{，}}\Big.$

$ \begin{array}{l} \qquad \overline B = \dfrac{{{B_1} + {B_2}}}{2} =\\ \dfrac{{\displaystyle\sum\limits_{i = 1}^m {{y_i}({t_i} - \overline {{t_1}} )} /\displaystyle\sum\limits_{i = 1}^m {{{({t_i} - \overline {{t_1}} )}^2}} + \displaystyle\sum\limits_{i = m + 1}^N {{y_i}({t_i} - \overline {{t_2}} )} /\displaystyle\sum\limits_{i = m + 1}^N {{{({t_i} - \overline {{t_2}} )}^2}} }}{2} {\text{。}} \end{array} $

SASDPE表示将整个控制图分为两部分，假如分割点为m，使两分段的最小二乘误差平方和PMSE最小，那么整体最小二乘回归线斜率(B)与两分段区域最小二乘回归线斜率之差绝对值的和就是SASDPE。

$\qquad{B_1} = \sum\limits_{i = 1}^m {{y_i}({t_i} - \overline {{t_1}} )}\Big /\sum\limits_{i = 1}^m {{{({t_i} - \overline {{t_1}} )}^2}} {\text{，}} \Big.$

$\qquad{B_2} = \sum\limits_{i = m + 1}^N {{y_i}({t_i} - \overline {{t_2}} )}\Big /\sum\limits_{i = m + 1}^N {{{({t_i} - \overline {{t_2}} )}^2}}\Big.{\text{，}}$

$\qquad{\rm{SASDPE}} = {{|B - }}{{{B}}_{\rm{1}}}{\rm{|}} + {{|B - }}{{{B}}_{\rm{2}}}{{|}} {\text{。}}$

对于阶跃模式(US/DS)来说，其SASDPE的值比趋势模式(UT/DT)高；MIX、CYC、SYS模式的SASDPE比NOR和STA模式高。

SASPE为两分段区域最小二乘回归线斜率绝对值之和， ${\rm SASPE} = |{B_1}| + |{B_2}| {\text{。}} $

REPEPE为整体平均最小二乘误差平方(MSE)与两分段最小二乘误差平方和(SE₁和SE₂)的最小平均值(PMSE)。

$\qquad{\rm REPEPE }={\rm MSE/PMSE} {\text{。}}$

RVPEPE为 ${\rm SD}^2$ 与PMSE的比值，

$\qquad{\rm RVPEPE} = \frac{{{{\rm SD}^2}}}{{{\rm PM}S{\rm E}}} {\text{。}}$

1.3 特征组合

上文共总结了14个统计特征和27个形状特征，不同的特征组合其分类效果不同。文献[15]从14个统计特征提取了3个统计特征。在1.1节中，离散度 $s$ 、均方幅值 ${X_{\rm ms}}$ 、脉冲 ${X_{\rm I}}$ ，证明其优于原始数据的分类识别；文献[16]提取了7个形状特征。在1.2节中，最小二乘线斜率标识值SB，模式与中心线组成的面积与观测值标准差之比ACLPI，中心线交点与最小二乘线交点均值与观测数量之比PSMLSC，经过分块区域中间点组合的斜率平均值ASL，6条直线斜率变化范围SRANGE，全控制图区域MSE与6个组合区域的MSE均值之比REAE，整体最小二乘回归线斜率与两分段区域最小二乘回归线斜率均值之差的绝对值ABDPE)，证明其分类识别效果较好。本文对这些特征进行了组合，形成了如下5种特征组合方案。方案1：7个形状特征；方案2：7个形状特征+3个统计特征；方案3：7个形状特征+14个统计特征；方案4：27个形状特征；方案5：27个形状特征+14个统计特征。

2 随机森林

随机森林(random forrest, RF)是利用多棵树对样本进行训练并预测的一种分类器，对于含有噪声及含缺失值的数据具有很好的预测精度，并且可以处理大量的输入变量，具有较快的训练速度，近些年来已经被广泛应用于分类、特征选择等诸多领域^[17]。

RF对于每棵树使用的训练集是从总训练集中有放回采样出来的，这意味着总训练集中有些样本可能多次出现在一棵树训练集中，也可能从未出现在一棵树训练集中。

随机森林和决策树的区别主要是：一方面随机森林中的决策树在分裂过程选择最佳特征进行分裂时，是从无放回抽取的特征中选择最优的特征进行分裂。这样能够使得随机森林中的决策树彼此之间差异性增大，提升系统的多样性，从而提升分类性能；另一方面是随机森林没有剪枝步骤，这样并不会导致过度拟合的问题。

随机森林的构建过程如下。

1) 从原始训练集中使用Bootstraping方法随机有放回采样选出m个样本，共进行n_tree次采样，生成n_tree个训练集。

2) 对于n_tree个训练集，分别训练n_tree个决策树模型。

3) 对于单个决策树模型，训练样本特征的个数为m_try，每次分裂时根据信息增益/信息增益比/基尼指数选择最好的特征进行分裂；每棵树都一直这样分裂下去，直到该节点的所有训练样例都属于同一类。

4) 将生成的多棵决策树组成随机森林,按多棵树分类器投票决定最终分类结果。

3 仿真验证 3.1 仿真数据产生

控制图中以 $\overline x $ 控制图最为常见，因此以 $\overline x $ 控制图为例进行仿真验证。采用蒙特卡罗方法仿真生产中实时测量值，公式如下^[18]：

$\qquad S(t) = D + r(t)\sigma {\rm{ + }}d(t) {\text{。}}$

$S(t)$ 是实时测量值； $D$ 是设计目标值； $d(t)$ 为过程异常因素导致的波动； $r(t)$ 为过程随机因素导致的波动； $\sigma $ 是过程标准差。

每种质量模式的仿真公式如下。

1) 正常模式(normal)。

$ \qquad S(t) = D + r(t)\sigma ,r(t)\sim N(0,1) {\text{。}} $

2) 分层模式(stratification)。

$\qquad S(t) = D + r(t)\sigma ',\;\;\sigma ' = 0.2\sigma {\rm{\sim }}0.4\sigma {\text{。}}$

3) 系统模式。

$ \qquad S(t) = D + r(t)\sigma + d{( - 1)^i},\;d = 1\sigma \sim 3\sigma {\text{。}} $

4) 循环模式(周期模式)。

$\begin{split} & \qquad S(t) = D + r(t)\sigma {\rm{ + }}A\sin \left(\frac{{2\pi t}}{T}\right) {\text{，}}\\ & \qquad {A = 1.5\sigma \sim 2.5\sigma ,\;T = 8,16} {\text{。}} \end{split} $

5) 上升趋势模式。

$\qquad S(t) = D + r(t)\sigma {{ + g}}i,\;g = 0.05\sigma \sim 0.1\sigma {\text{。}}$

6) 下降趋势模式。

$\qquad S(t) = D + r(t)\sigma {{ - g}}i,g = 0.05\sigma \sim 0.1\sigma {\text{。}}$

7) 上升阶跃模式。

$\qquad S(t) = D + r(t)\sigma {{ + }}ks,$

$ s = 1.5\sigma \sim 2.5\sigma $ ，阶跃位置 $ P = 9,17,25 $ ，

$ \qquad {{\rm if}\;i {\text{＜}} P,\;k = 0;\;{\rm else}}\;{k = 1} {\text{。}} $

8)下降阶跃模式。

$ \qquad S(t) = D + r(t)\sigma {{ - }}ks {\text{，}} $

$s = 1.5\sigma \sim 2.5\sigma $ ，阶跃位置 $P = 9,17,25$ ,

$ \qquad {{\rm if}\;i {\text{＜}}P,\;k = 0;\;{\rm else}}\;{k = 1}{\text{。}} $

式中， $g$ 为上升或下降趋势的变化斜率；P为阶跃位置，用 ${t_0}$ 表示阶跃时刻，则阶跃前 $k = 0$ ，阶跃后 $k = 1$ ； $s$ 为阶跃幅度； $T$ 为周期模式的变化周期； $A$ 为周期变化幅度； $\sigma $ 取实际加工过程稳态时的标准差。

实验数据生成、特征提取和多分类支持向量机的训练和测试均在Matlab R2016a软件中编程实现，程序运行的硬件环境为2.3 GHz的CPU和8 GB内存。

3.2 仿真过程及结果

为了便于比较，仿真验证过程参考文献[19]，设置均值D=80，标准差 $\sigma $ =5，训练样本和测试样本分别为100×8个，即共8种控制图模式，每种类别中包含100个训练样本和100个测试样本，共测试10次。分类器选择支持向量机(SVM)、人工神经网络(ANN)、决策树(CART)和随机森林(RF)。特征选择上文所述5种组合方案。统计4种分类器在5种不同特征组合方案下得到的分类准确率和耗费CPU时间.其中支持向量机采用网格法参数寻优；人工神经网络采用误差反向传播算法；决策树采用基尼指数选择最优特征，采用交叉验证误差确定最小叶片尺寸；随机森林中n_tree=20，其他参数默认。测试结果如表2和图4所示。

从上面的图表可以看出，在任何一类特征下，随机森林和支持向量机的分类准确率都高于人工神经网络和决策树，而支持向量机前两种特征下分类准确率高于随机森林；后3种特征下随机森林略高于支持向量机。但从耗费时间来看，随机森林远远低于其他3种分类器，而且随着特征增多时间增长比较缓慢，见图5。因此综合分类准确率和消耗时间，随机森林对于控制图模式识别来说都是性能优异的分类器，其运算速度快，更适合实时在线模式识别。

4 总结

控制图模式的正确识别对于及时分辨生产中的异常状态、保持产品的高质量产出非常重要。随机森林是一种灵活且易于使用的机器学习算法，不需要参数调优，也可以在大多数情况下取得很好的预测效果。本文将其应用于控制图模式识别中，并与其他3种常用的控制图模式识别方法(支持向量机、人工神经网络、决策树)进行比较；通过蒙特卡罗仿真实验方法验证，随机森林方法在5种特征组合情况下，相比于其他3种分类方法在分类准确率和消耗时间方面都取得了较优异的结果，证明该方法适用于一般统计过程控制图的模式识别，不仅可以得到较高的分类准确率，而且还能够缩短识别时间，在智能制造的背景下更适合控制图模式的快速识别。