GB/T 40859-2021 流式数据监测控制图

　　ICS 03 . 120 . 30 CCS A 4 1

　　中 华 人 民 共 和 国 国 家 标 准

　　GB/T 40859—2021

　　流式数据监测控制图

　　Controlchartforstreamingdatamonitoring

　　2021-10-1 1 发布 2022-05-01 实施

　　国家市场监督管理总局国家标准化管理委员会

　　发

　　布

　　GB/T 40859—202 1

　　GB/T 40859—202 1

　　前 言

　　本文件按照 GB/T 1 . 1—2020《标准化工作导则 第 1 部分：标准化文件的结构和起草规则》的规定起草。

　　请注意本文件的某些内容可能涉及专利。 本文件的发布机构不承担识别专利的责任。

　　本文件由全国统计方法应用标准化技术委员会(SAC/TC 21)提出并归口 。

　　本文件起草单位：清华大学、中国标准化研究院。

　　本文件主要起草人：王凯波、杨洛、丁文兴、张帆、赵静。

　　GB/T 40859—202 1

　　引

　　言

　　在现代复杂生产系统和装备运行过程中，由于自动化传感器的广泛应用，所获取的数据往往是流式数据的形式。 所谓流式数据，是一组数量庞大、顺次有序、快速高频到达的数据序列。 现有的控制图往往基于较大间隔进行离散采样获得的数据进行过程控制。 由于流式数据具有快速高频等特点，现有控制图通常无法适用于控制产生流式数据的过程。

　　传统的数据采集与分析往往是低频的，如每小时对一条生产水杯的生产线进行产品采样和关键质量指标测量，以判断生产过程的状态是否正常。 而在复杂系统的生产制造和运行过程中，通常会布局很多设备状态传感器和自动化检测仪器，以秒甚至更小的时间单位为间隔，对设备运行状态及产品质量指标进行实时、自动和持续的数据收集，由此产生流式数据。 这些流式数据包含比低频数据更加复杂、全面的信息，但其信息密度也相对较低，因此，需要合适的特征提取手段，挖掘流式数据中包含的信息。

　　20 世纪，有学者提出了多分辨率分析方法。 这是一种基于离散小波变换的数据分析方法，可以将数据分解到不同的频率尺度下进行分析。 此后，小波分析的理论方法不断完善。 多分辨率小波分析十分适用于处理高频的流式数据。 将高频的流式数据分解到不同的频率尺度后，可在不同的尺度下分别分析，提取流式数据中复杂的信息。

　　本文件可对高速铁路、装备制造、复杂产品制造等系统在运行过程中所产生的实时状态数据和产品质量数据进行监测，对关键信号的异常变化进行预警，有可能在早期发现危险故障，提升系统安全可靠性，提高产品质量。

　　GB/T 40859—202 1

　　流式数据监测控制图

　　1 范围

　　本文件给出了对单变量及多变量流式计量型数据进行监测的控制图方法。

　　本文件适用于对产生流式数据的一个或多个有一定关联关系的过程或质量变量，在均值和基本波动特征已知的条件下，对偏移和波动进行控制。

　　2 规范性引用文件

　　下列文件中的内容通过文中的规范性引用而构成本文件必不可少的条款。 其中，注 日期的引用文件，仅该日期对应的版本适用于本文件;不注日期的引用文件，其最新版本(包括所有的修改单)适用于本文件。

　　GB/T 3358 . 1 统计学词汇及符号 第 1 部分：一般统计术语与用于概率的术语

　　GB/T 3358 . 2 统计学词汇及符号 第 2 部分：应用统计

　　3 术语和定义

　　GB/T 3358 . 1 和 GB/T 3358 . 2 界定的以及下列术语和定义适用于本文件。

　　3.1

　　流式数据 streamingdata

　　一组数量庞大、顺次有序、快速高频到达的数据序列。

　　注 1 ：传感器监测数据、交通监测数据等往往都是以流式数据的形式被收集。

　　注 2：流式数据通常是由一个或多个数据源持续产生，随时间延续而增长的动态数据集合，具有采样间隔小、采样频率高、到达速度快、数据量大等特点。

　　3.2

　　EWMA控制图 EWMA chart

　　指数加权移动平均控制图 exponentiallyweightedmovingaveragecontrolchart

　　用指数平滑加权平均评估和监测过程水平的计量控制图。

　　[来源：GB/T 3358 . 2—2009,2 . 3 . 16,有修改]

　　3.3

　　EWMS控制图 EWMSchart

　　指数加权均方误差控制图 exponentiallyweightedmeansquarederrorcontrolchart

　　用指数平滑加权平均评估和监测过程偏移与变异的计量控制图。

　　3.4

　　MEWMA控制图 MEWMA chart

　　多变量指数加权移动平均控制图 multivariateexponentiallyweightedmovingaveragecontrolchart

　　将两个或两个以上相关的变量合成一个样本统计量，并用指数平滑加权平均评估和监测过程水平的计量控制图。

　　GB/T 40859—202 1

　　3.5

　　MEWMC控制图 MEWMCchart

　　多变量指数加权移动协方差矩阵控制图 multivariateexponentiallyweighted movingcovariance matrixcontrolchart

　　将两个或两个以上相关的变量合成一个样本统计量，并用指数平滑加权平均评估和监测过程变异的计量控制图。

　　3.6

　　小波变换 wavelettransform

　　用样本数据与小波函数的内积计算得到小波系数的过程。

　　3.7

　　离散小波变换 discretewavelettransform

　　用样本数据与正交小波基的内积计算得到小波系数的过程。

　　3.8

　　小波基函数 waveletbasisfunction

　　母小波 motherwavelet

　　均方可积且均值为 0 的函数。

　　注 1 ：常用小波基函数有 Haar小波、Meyer小波等。

　　注 2：小波基函数通常表示为 ψ(t)。

　　3.9

　　小波函数 wavelets

　　用小波基函数伸缩与平移得到的一组函数。

　　注：小波函数可表示为 ψj，k(t)=a-j/2 ψ(a-jt-kb)，其中 a 为伸缩参数，b为平移参数，j为尺度参数，k 为位置参数。

　　3 . 10

　　正交小波基 orthonormalwaveletbases

　　用某些小波基函数进行特定的离散化伸缩与平移，得到的可作为正交基的小波函数。

　　注 1 : Haar小波是可用于产生正交小波基的一种小波基函数。

　　注 2：通常选取伸缩参数 a=2,平移参数 b= 1。

　　3 . 1 1

　　小波系数 waveletcoefficient

　　用样本数据与小波函数内积的结果。

　　注 ：内积过程可表示为 dj，k= ∫f(t)ψj，k(t)e-t2 dt，其中 f(t)为样本数据。

　　3 . 12

　　父小波 fatherwavelet

　　尺度函数 scalingfunction

　　可生成正交基的范数为 1 的规范化函数，与母小波对应。

　　注 1 ：并不是所有的母小波都有对应的父小波。

　　注 2：可用于离散小波变换的母小波都有对应的父小波。

　　注 3：母小波用于表现样本波动信息，父小波用于表现样本近似信息。

　　3 . 13

　　多分辨率分析 multiresolutionanalysis

　　用离散小波变换得到的小波系数进行数据分析的方法。

　　注：离散小波变换得到的小波系数包含不同频率尺度下的信息。

　　GB/T 40859—202 1

　　3 . 14

　　近似层小波系数 approximatewaveletcoefficient

　　用样本数据与父小波内积的结果。

　　注 1 ：近似层小波系数包含样本趋势信息。

　　注 2：样本经离散小波变换通常仅保留一层近似层小波系数。

　　3 . 15

　　细节层小波系数 detailwaveletcoefficient

　　用样本数据与母小波内积的结果。

　　注 1 ：细节层小波系数包含样本波动信息。

　　注 2：样本经离散小波变换会产生若干层近似层小波系数。

　　4 符号

　　下列符号适用于本文件。

　　GB/T 40859—202 1

　　σc 近似层小波系数标准差的目标值

　　σdj 第 j 层细节层小波系数标准差的目标值

　　Σc 近似层小波系数向量均值协方差矩阵的目标值

　　Σdj 第 j 层细节层小波系数向量均值协方差矩阵的目标值

　　ψ 母小波

　　5 应用条件

　　现实的生产过程中，由于各种传感器的大量使用，过程样本数据往往以流式数据的形式被收集。 流式数据具有采样频率高、采样间隔短、数据量大等特点。 由于数据获取频率高，短时间内便能获得大量数据，提高了数据计算和存储的难度，故流式数据需采用新的监测方法。

　　此外，生产过程中的过程信号往往十分复杂，可能同时存在着多种不同的特征信号。 这些复杂的信号组合在一起会导致难以对故障的原因进行诊断，更严重的是有时多种特征信号的叠加可能会导致 一些特征信号难以被监测。 多种特征信号叠加的数据如图 1 所示。

　　说明：

　　xA — 白噪声;

　　xB —包含均值突变的数据;

　　xC —包含波动突变的数据;

　　xD —xA 、xB 、xC 叠加后的数据。

　　图 1 复杂的特征信号

　　从图 1 中可看出，xB 发生了均值的偏移，xC 发生了波动幅度的变化。 虽然这二者变化特征明显，但是对于由 xA 、xB 、xC 叠加得到的数据 xD,其均值偏移与波动变化两种特征信号互相影响。 若使用单一控制图监测 xD,则对均值偏移与波动变化都难以进行有效识别。

　　若可以将图 1 中叠加在一起的复杂信号分开进行监测，那么便会大大提高监测的效果。 离散小波变换是一种对数据进行分解的方法，原始数据经过离散小波变换后，可在不同的频率尺度上监测不同的特征信号。 此外，即使原始数据存在一定自相关性，经离散小波变换得到的小波系数的自相关性会显著

　　GB/T 40859—202 1

　　降低或消失，亦即离散小波变换所得到的系数对原始数据中存在的自相关性不敏感，因此，离散小波变换的方法可适用于流式数据的监测中。

　　在使用本控制图对过程进行监测时，应有一定量能够代表过程正常运行状态的历史数据，用于计算控制图参数的目标值。

　　6 基于离散小波变换的数据预处理

　　为了实现对流式数据样本的监测，需要首先对所获取数据应用基于 Haar 小波的离散小波变换算法进行预处理，生成小波系数，然后将小波系数用于控制图监测。

　　Haar小波是一种常用的正交小波，其定义为：

　　…………………………( 1 )

　　结合附录 A 中介绍的小波函数的公式(A. 1),可知其在第 j 层的小波函数的定义域长度仅为 2j，因此仅使用有限的数据点即可计算一个小波系数。 假定分解层数为 J，流式数据序列为 x={x1 , x2 ,…} 。近似层小波系数和细节层小波系数的计算公式分别为：

　　cJ，k 2J+i …………………………( 2 )

　　dj，k ………………

　　在监测过程中，需对新获取的数据点进行计数。 假定计数为 i，当 i 可以整除 2j 时，计算前 j 层小波系数。 表 1 给出了应用小波变换得到细节层小波系数的示例。 表中每一列表示一个新获取的数据，每一行表示小波分解的层，表内的 dj，k表示第j 层的第k 个小波系数。 表中共展示了 8 个数据点，进行3 层小波分解的系数。 表中可看出，当计数为 2 时，可通过 1 、2 两个数据点计算获得一个第一层的小波系数 d1 , 1 ;计数为 8 时，可通过 1~8 八个数据点计算第三层小波系数 d3 , 1 ,同时通过 5、6、7、8 四个数据点计算第二层小波系数 d2 , 2 ,通过 7、8 两个数据点计算第一层小波系数 d1 , 4 。当获取新的数据时，可重复该过程，不断获取新的小波系数。

　　表 1 小波变换样例

　　通过该算法，可将实时获得的数据及时分解为各层小波系数，共可得到 1 层近似层系数{cJ，1 , cJ，2 , …}以及 J 层细节层小波系数{dj，1 ,dj，2 ,…} ,j∈ { 1 , … , J｝，用于对过程的监测。 在构建控制图的过程中，分解层数 J 的选择取决于流式数据特征信号的频率范围，若要监测较低频的波动信号则选取较大的分解层数 J。但分解层数 J 过大会导致计算复杂度增加，影响监测效率。 通常分解层数 J 建议选为2 或 3 即可。

　　在应用控制图时，往往假设被监测特征服从正态分布。 在流式数据服从正态分布的假设下，基于Haar小波的数据预处理得到的各层小波系数仍然满足正态分布的假设，可使用控制图进行监测。 在过

　　GB/T 40859—202 1

　　程受控状态下获取的数据经预处理得到的各层小波系数，可用于计算控制图参数的目标值。

　　7 单变量流式数据控制图的构建

　　7 . 1 概述

　　流式数据经小波变换分解为小波系数后，使用 EWMA 控制图监测近似层小波系数，使用 EWMS控制图监测各层细节层小波系数。 附录 B 中给出了构建控制图的案例。

　　7 . 2 单变量流式数据 EWMA控制图的构建

　　EWMA控制图适用于监测数据的均值变化。 应使用经离散小波变换得到的近似层小波系数{ cJ，1 , cJ，2 ,…}构建该控制图 。

　　EWMA统计量可表示为：

　　EWMA控制图的控制限为：

　　UCL=μc+L …………………………( 6 )

　　LCL =μc -L …………………………( 7 )

　　其中，L 为控制限系数。 应用 ARL分析的方法可以确定控制图的两个参数：控制限系数 L 和平滑系数λ1 。通常，λ1 的取值取决于偏移量的大小，一般推荐范围为 0 . 05≤λ1≤0 . 5 。ARL0 = 370 时的参数选择如表 2 所示。

　　表 2 ARL0 = 370 时 EWMA控制图参数选择标准

　　7 . 3 单变量流式数据 EWMS控制图的构建

　　EWMS控制图适用于监测数据的方差变化。 应使用经离散小波变换得到的各层细节层小波系数

　　{dj，1 ,dj，2 ,…} ,j∈ { 1 , … ,J｝，共构建 J个 EWMS控制图。

　　其中第 j 个控制图的 EWMS统计量可表示为：

　　s，i= λ2(dj，i-μdj)2 +(1-λ2)s，i- 1 …………………………( 8 )

　　其中 λ2 ∈ (0,1)为平滑系数，s，0 = σj 。

　　通常，绘制关于 sj，i的控制图，而非关于 s，i的控制图，即将迭代序列的平方根绘制入控制图。

　　EWMS控制图的控制限为：

　　UCL = σdj …………………………( 9 )

　　GB/T 40859—202 1

　　LCL = σdj槡 …………………………( 10 )

　　其中，α 表示发生第一类错误的概率，χ- α/2) ×100%和(α/2) ×100%分位数，ν=(2-λ2)/λ2 。λ2 的取值取决于偏移量的大小，一般推荐范围为 0.01≤λ2≤0.4。ARL0 =370 时 (α=0.002 7),平滑系数 λ2 与控制限的对应关系如表 3 所示。

　　表 3 ARL0 = 370 时 EWMS控制限选择标准

　　7 . 4 单变量流式数据控制图应用步骤

　　单变量流式数据控制图构建的具体步骤如下。

　　步骤 1 确定控制图参数：选取离散小波变换分解层数 J、EWMA 控制图的控制限 L 和平滑系数λ1、EWMS控制图的第一类错误的概率 α 和平滑系数λ2 ,并计算各控制图的控制限。

　　步骤 2 计算小波系数：对新获取的数据 xi 进行计数，代入公式(2)、公式(3) 获取新的近似层或细节层小波系数。

　　步骤 3 计算统计量：使用新获得的小波系数，近似层小波系数代入公式(4) 计算 EWMA 统计量，各层细节层小波系数分别代入公式(8)计算各层的 EWMS统计量。

　　步骤 4 绘制控制图：将各统计量及相应的控制限分别绘制在控制图上，共得到一张 EWMA 控制图及 J 张 EWMS控制图。 根据控制图判断受控状态。 当任何一张控制图出现控制点超过控制限时，表明过程可能“失控”，需要对过程采取必要的措施，予以纠正。

　　8 多变量流式数据控制图的构建

　　8 . 1 概述

　　假设某 一 时 刻 开 始 获 得 多 变 量 流 式 数 据 序 列 X ={ x1 , x2 , … } , 变 量 的 维 度 为 p，即 x1 = (x1(1) ,…,x1(p))T 。为达到监测 目的，首先使用离散小波变换将每一个变量的数据分解为近似层小波系数和细节层 小 波 系 数，再 将 不 同 变 量 同 一 层 的 小 波 系 数 重 组，构 成 多 维 的 小 波 系 数 向 量。 使 用MEWMA控制图监测近似层小波系数向量，使用 MEWMC控制图监测各细节层小波系数向量。 附录 C 中给出了构建控制图的案例。

　　8 . 2 小波系数重组

　　使用离散小波变换将每一维数据都分解至 J 层 。 以 2 维数据为例：第 1 维数据{x1(1) , x2(1) ,…}可分解得到 1 层近似层系数{ cJ(1,)1 , cJ(1,)2 , …} 以及 J 层细节层小波系数{ dj(1, , dj(1,2) , …} , j ∈ { 1 , … , J｝;第 2 维数据{ x1(2) , x2(2) , …} 可 分 解 得 到 1 层 近 似 层 系 数{ cJ(2,)1 , cJ(2,)2 , … } 以 及 J 层 细 节 层 小 波 系 数{dj(2, ,dj(2,2) ,…} ,j∈ { 1 , … ,J｝。

　　将 2 维数据的近似层系数组成向量，可以得到 1 层近似层系数向量{CJ，1 , CJ，2 , …},其中 CJ，1 = (cJ(1,)1 , cJ(2,)1 );将 2 维数据对应的细节层系数组成向量可得到 J 层细节 层 小 波 系 数{dj，1 , dj，2 , … } , j∈ { 1 , … ,J｝，其中 dj，1 =(dj(1, ,dj(2, )。

　　图 2 展示了变量维度为 p 时小波系数的重组过程。

　　GB/T 40859—202 1

　　图 2 多变量流式数据控制图小波系数重组示意图

　　8 . 3 多变量流式数据 MEWMA控制图的构建

　　MEWMA控制图适用于监测多变量流式数据的均值变化。 应使用 8 . 2 中经离散小波变换并重组得到的近似层小波系数向量{CJ，1 , CJ，2 ,…}构建控制图。

　　定义向量序列 zi，z0 =uc。每获取一个新的小波系数向量 CJ，i，便可通过迭代计算得到 zi：

　　c …………………………( 13 )

　　应用 ARL分析的方法可以确定控制图的两个参数：控制限hA 和平滑系数λ3 。λ3 的取值取决于偏移量的大小，一般推荐范围为 0 . 05≤λ3≤0 . 5 。应选择合适的 hA 使得控制图可达到一定的受控 ARL, ARL0 = 370 时 hA 选择的原则如表 4 所示。

　　表 4 ARL0 = 370 时 MEWMA控制图控制限选择标准

　　GB/T 40859—202 1

　　8 . 4 多变量流式数据 MEWMC控制图的构建

　　MEWMC控制图适用于监测多变量流式数据的方差及协方差变化。 应使用 8 . 2 中经离散小波变换并重组得到的各层细节层小波系数向量{dj，1 ,dj，2 ,…} ,j∈ { 1 , … , J｝构建控制图。

　　构建第 j 个 MEWMC 控制图时，定义向量序列 sj，i，sj，0 = Ip 。每获取一个新的小波系数向量dj，i，先进行一步转换：

　　其中，yj，i为该控制图的统计量，p为变量的维度，hc 为 MEWMC控制图的控制限，tr( )为求矩阵迹的函数。 构建该控制图同样需要确定两个参数：控制限 hC 和权重系数 λ4 。λ4 的取值取决于偏移量的大小，一般推荐范围为 0 . 01≤λ4≤0 . 4 。应选择合适的 hC 使得控制图可达到一定的受控 ARL, ARL0 = 370 时 hC 选择的原则如表 5 所示。

　　表 5 ARL0 = 370 时 MEWMC控制图控制限选择标准

　　8 . 5 多变量流式数据控制图应用步骤

　　多变量流式数据控制图构建的具体步骤如下：

　　步骤 1 确 定 控 制 图 参 数：选 取 离 散 小 波 变 换 分 解 层 数 J、MEWMA 控 制 图 的 平 滑 系 数 λ3 、 MEWMC控制图的平滑系数 λ4,并选择各控制图的控制限。

　　步骤 2 计算小波系数向量：对新获取的数据 xi 进行计数，并分别代入公式(2)、公式(3) 获取新的近似层或细节层小波系数，再将对应层小波系数重组为小波系数向量。

　　步骤 3 计算统计量：使用新获得的小波系数向量，近似层小波系数向量代入公式(11)、公式(12)和公式(13)计算 MEWMA统计量，各层细节层小波系数向量分别代入公式(14)、公式(15) 和公式(16)计算各层的 MEWMC统计量。

　　步骤 4 绘制控制图：将各统计量及相应的控制限分别绘制在控制图上，共得到一张 MEWMA 控制图及 J 张 MEWMC控制图。 根据控制图判断受控状态。 当任何一张控制图出现控制点超过控制限时，表明过程可能“失控”，需要对过程采取必要的措施，予以纠正。

　　以上步骤往往需要使用计算机辅助实现，实现控制图的 Python代码在附录 D 中列出。

　　GB/T 40859—202 1

　　附 录 A

　　(资料性)

　　离散小波变换及多分辨率分析原理

　　多分辨率数据分析的原理是利用离散小波变换将数据分解到不同的尺度上进行分析，这对于流式数据的监测十分有效。 离散小波变换是一种时频分析算法，与传统的傅里叶变换不同，它能保留原始数据在时域上的信息，即可将原始数据分解到不同的频率尺度上。

　　用于过程监测的流式数据通常可通过小波函数分解成近似信号和细节信号的形式：

　　J

　　fcJ，kφJ，k dj，kψj，k ……………………

　　其中，j和 k分别为尺度参数和位置参数，φ 是低频的小波基父函数，ψ是高频的母小波;cJ，k和dj，k分别为近似小波系数和细节小波系数。 其表达式为：

　　φJ，k ，cJ，k = φJ，k e-t2 dt ………… ψj，k(t)=t-k)，dj，k=〈f，ψj，k〉= ∫f(t)ψj，k(t)e-t2 dt …………( A.3 )

　　当 k=m2j，m ∈ Z 时，这一小波变换的方法被称作离散小波变换。 在分解时，首先确定分解层数J，之后依次计算各层的小波系数，在每层中以 k=i2j，i∈Z逐步平移小波基函数。 共可得到 1 层的近似小波系数和 J 层细节小波系数，如图 A. 1 所示。

　　图 A.1 离散小波变换小波系数

　　每一个小波系数都是由一段原始数据与该层小波基函数的内积得到，即该小波系数可以描述原始数据在一段时间内与该层小波基函数的相似程度。 因此不同层的小波系数描述了原始信号在不同频域上的信息，如图 A. 2 所示。

　　GB/T 40859—202 1

　　说明：

　　狓 —原始数据;

　　d1 —第一细节层小波系数;

　　d2 —第二细节层小波系数;

　　c2 — 近似层小波系数。

　　图 A.2 小波变换效果

　　从图 A. 2 中可看出，原始信号 狓 包含有均值偏移和方差变化两种不同的特征信号。 利用离散小波变换将其分解到 2 层后，可看出均值偏移的信息表现在近似小波系数上，而波动变化的信息表现在细节系数上。 因此，可对分解后的各层小波系数分别进行监测，这便是多分辨率分析的原理。

　　GB/T 40859—202 1

　　附 录 B

　　(资料性)

　　单变量流式数据应用实例

　　在某直拉单晶硅生产过程中，需要对单晶的关键质量参数进行监测。 记该质量参数为变量 x。单晶拉制过程中，数据的获取间隔约 1 s。 生产系统已经平稳运行了一段时间，获取了足够长的平稳数据X。本案例中相关数据已进行脱敏处理。

　　在建立控制图前首先要基于历史平稳数据对控制图参数的目标值进行估算。 当所获得历史数据存在异常值时，基于工程经验判断和统计分析，对异常值进行剔除，仅使用平稳数据进行参数的 目标值估计 。将过程平稳状态下获取的平稳数据 X 应用离散小波变换分解到 2 层，得到一层近似层小波系数{ c2 , 1 , c2 , 2 ,…} 和两层细节层小波系数{d1 , 1 ,d1 , 2 ,…}和{d2 , 1 , d2 , 2 ,…} 。根据这些小波系数可估计出各层小波系数的均值和方差的目标值，计算结果如表 B. 1 所示。

　　表 B.1 各层小波系数均值与方差

　　现从该过程获取新的 1 000 组数据 x={x1 ,…,x1000} ,如图 B. 1 所示。 下面对这组新的数据实施控制图监测。

　　o200400600 8001000

　　图 B.1 流式数据 x

　　GB/T 40859—202 1

　　对数据 x={x1 ,…,x1000}应用公式(2)、公式(3),并代入 j=2, 即可使用离散小波变换将数据分解到 2 层，同样得到一层近似层小波系数{c2 , 1 ,c2,2,…,c2,250} 和两层细节层小波系数{d1 , 1 ,d1 , 2 ,…,d1,500}和{d2 , 1 , d2 , 2 ,…, d2,250} 。

　　对于近似层小波系数，使用 EWMA 控制图监测变量的均值，选择平滑系数 λ1 =0 . 1, 控制限系数L= 3 。将{c2 , 1 , c2 , 2 ,…,c2,250} 和 μc =0 . 00 代入公式(4) ,即可得到 EWMA控制图的统计量 义i。将 μc =

　　0 . 00 和 σ =0 . 149 7 代入公式(6)、公式(7) 可计算 EWMA 控制图的控制限。 将统计量 义i 和控制限绘

　　制成图，如图 B. 2 所示。

　　02004006008001000时间/s

　　图 B.2 近似层小波系数 EWMA控制图

　　从控制图中可看出，EWMA统计量在第 580~904 秒超出了控制限，可判断在此段时间内流式数据的均值发生了变化。

　　对于两层细节层小波系数，分别使用 EWMS控制图监测变量的方差，选择平滑系数 λ2=0 . 1,控制

　　限系数 α=0.002 7。将{d1 , 1 , d1 , 2 ,…, d1,500} 、μd1 = 0 . 00 、σ1 = 0 . 003 8 和{d2 , 1 , d2 , 2 ,…, d2,250} 、μd2 =

　　0 . 00 、σ2 =0 . 007 5 分别代入公式(8) ,即可得到 EWMS控制图的统计量 S1,i和 S2 ,i。将 λ2 = 0 . 1 代入公

　　式(9)、公式(10)可计算 EWMS 控制图的控制限。 将统计量 S1,i、S2,i和控制限绘制成图，如图 B. 3 和图 B. 4 所示。

　　GB/T 40859—202 1

　　图 B.3 第 1 层细节层小波系数 EWMS控制图

　　图 B.4 第 2 层细节层小波系数 EWMS控制图

　　从图 B. 3 中可看出，第 1 层细节层小波系数的 EWMS 统计量在第 158 ~ 240 秒、第 278~298 秒 、第 490~678 秒、第 738~786 秒和第 892~962 秒内出现了超出控制限的情况;从图 B. 4 中可看出，第 2层细节层小波系数的 EWMS统计量在第 468~500 秒、第 560~600 秒和第 676~708 秒内出现了超出控制限的情况。 由于两个控制图分别包含原始数据在不同频域上的信息，应将两个控制图的结果结合在一起考虑，最终可判断流式数据在上述时间段内变量的方差发生了变化。

　　GB/T 40859—202 1

　　附 录 C

　　(资料性)

　　多变量流式数据应用实例

　　在某直拉单晶硅生产过程中，涉及两个质量相关的变量 A 和 B(该数据已进行数据脱敏)，过程中数据的获取间隔是 1 s。 生产系统已经平稳运行了一段时间，获取了足够长的平稳数据 X。

　　在建立控制图前首先要对控制图参数的目标值进行估算。 当所获得历史数据存在异常值时，基于工程经验判断和统计分析，对异常值进行剔除，仅使用平稳数据进行参数的目标值估计。 将过程平稳状态下获取的平稳数据 X 的每一个变量都用离散小波变换分解到 2 层，并重组得到一层近似层小波系数向量{c2, 1 , c2,2 ,…} 和两层细节层小波系数向量{d1 , 1 ,d1 ,2 ,…} 和{d2 , 1 ,d2 ,2 ,…} 。根据这些小波系数向量可估计出各层小波系数的均值和方差的目标值，近似层小波系数向量均值为 μc=[0.00,0.00],协方

　　差矩阵为 Σc= [0.149 7 , -0.010 2 ; -0.010 2 , 0.057 9]。两层细节层小波系数向量均值为 μd1 = μd2 = [0.00,0.00], 协方差矩阵为 Σd1 = [0. 003 8 , 0. 000 0 ; 0. 000 0 , 0. 040 5] , Σd2 = [0. 007 5 , - 0. 001 7 ;

　　-0 . 001 7 ,0 . 052 0] 。

　　现获取了 1 000 组新的数据 x={x1 ,x2,…,x1000},其中 x1 =[x1(A)，x1(B)]。两个变量如图 C. 1 所示。

　　图 C.1 变量示意图

　　对于获取的数据 X ={x1 , x2 ,…,x1000}进行控制图监测。 应用公式(2)、公式(3),并代入 j=2, 即可使用小波变换分别将两个变量的数据分解到 2 层，每个变量都可得到一层近似层小波系数{c2 , 1 , c2 , 2 , …, c2,250} 和两层细节层小波系数{d1 , 1 , d1 , 2 ,…, d1,500} 和{d2 , 1 , d2 , 2 ,…, d2,250} 。将两个变量对应层的小波系数重新组合，则每一层都包含两个变量的小波系数，得到小波系数向量{c2, 1 , c2,2 , … , c2,250 } 、 {d1 , 1 , d1,2,…, d1,500} 和{d2 , 1 , d2,2,…, d2,250} 。

　　对于近似层小波系数向量，使用 MEWMA 控制图监测变量的均值，选择平滑系数 λ3 = 0 . 1 。将{c2, 1 , c2,2 , … , c2,250 } 和 μc 代入公式 (11) , 即可得到迭代向量 zi。再将 zi 和 Σc 代入公式(12)、公

　　式(13) ,即可得到 MEWMA控制图的统计量 y。 根据表 4 可知，该控制图的控制限 hA = 10 . 08 。将统

　　计量 y 和控制限hA 绘制成图，如图 C. 2 所示。

　　GB/T 40859—202 1

　　图 C.2 近似层小波系数 MEWMA控制图

　　从控制图中可看出，MEWMA统计量在第 428~656 秒、第 684~806 秒和第 888~904 秒内超出了

　　控制限，可判断在此段时间内数据中的某一个或某几个变量发生了均值的变化。

　　对于两层细节层小波系数向量，使用 MEWMC 控制图监测变量的方差和协方差，选择平滑系数 λ4 =0 . 1 。将{d1 , 1 ,d1,2,…,d1,500} 、μd1 、Σd1 和{d2 , 1 , d2,2,…, d2,250} 、ud2 、Σd2 分别代入公式 ( 14 ) ,即可得到转换后的向 量 序 列 U1 ,i 和 U2 ,i。再 将 λ4 = 0 . 1 、U1 ,i 和 U2,i代 入 公 式 (15)、公 式 (16) , 即 可 计 算MEWMC控制图的统计量 y1,i和 y2 ,i。根据表 5 可知，该控制图的控制限 hC = 0 . 63 。将统计量 y1 ,i 、 y2,i和控制限hC 绘制成图，如图 C. 3 和图 C. 4 所示。

　　图 C.3 第一层细节层小波系数 MEWMC控制图

　　GB/T 40859—202 1

　　图 C.4 第二层细节层小波系数 MEWMC控制图

　　从图 C. 3、图 C. 4 中可看出，第一层细节层小波系数的 MEWMC统计量在第 246~250 秒、第 340~ 358 秒、第 478~588 秒、第 610~654 秒、第 774~812 秒和第 956~988 秒内出现了超出控制限的情况;

　　第二层细节层小波系数的 MEWMC统计量在第 560~592 秒、第 624~680 秒和第 964~976 秒内出现了超出控制限的情况。 由于两个控制图分别包含原始数据在不同频域上的信息，因此将两个控制图的结果结合在一起考虑，最终可判断原始数据在上述时间段内某一个或某几个变量的方差或几个变量间的协方差发生了变化。

　　GB/T 40859—202 1

　　附 录 D

　　(资料性)代码实现

　　D.1 概述

　　为了实现对生产过程的监测，需要使用计算机的辅助。 且在本文件介绍的算法中，小波变换、矩阵计算等复杂的数学计算往往也可通过编程软件的内置函数进行。 Python 是目前应用最广泛的数据处理语言之一，以下是使用 Python实现绘制两种控制图的主要代码。

　　D.2 单变量流式数据控制图构建算法

　　建立 EWMA类：

　　import math

　　import numpy as np

　　import matplotlib.pyplot as plt class EWMA() :

　　" " "

　　这是用于建立 EWMA控制图的类

　　" " "

　　def _init_(self, lambda_1 = 0.1 , mu_c = 0 , sigma_c = 1 , level = 1 , L = 2.715) :

　　" " "

　　初始化 EWMA类

　　参数：

　　lambda_1：平滑系数

　　mu_c： 均值目标值

　　sigma_c：方差目标值

　　level： 小波系数的分解层数

　　L： 控制限系数

　　" " "

　　self.lambda_ 1 = lambda_ 1

　　self. mu_c = mu_c

　　self.sigma_c = sigma_c

　　self.level = level

　　self.value = [] #记录所有的统计量的值

　　self.old_number = mu_c #用于表示 z(i-1), 其中 z(0) = mu_c

　　self.new_number = 0 #用于表示 z(i)

　　self.ucl = self.mu_c + L* (self.sigma_c* self.lambda_1/(2-self.lambda_1)) ** 0.5 #上控

　　制限

　　self.lcl = self.mu_c - L* (self.sigma_c* self.lambda_1/(2-self.lambda_1)) ** 0.5 #下控

　　制限

　　def update(self, data) :

　　GB/T 40859—202 1

　　" " "

　　获取新的小波系数后计算对应的控制图统计量

　　参数：

　　data： 新获取的小波系数

　　" " "

　　n = len(data) #新获取的小波系数的总长度

　　for i in range(0, n) :

　　self.new_number = self.lambda_1 * data[i] + (1-self.lambda_1) * self.old_number #

　　对应公式(4)

　　self.value.append(self.new_number)

　　self.old_number = self.new_number

　　def plot(self, length) :

　　" " "

　　将最新获取的部分数据绘制成控制图

　　参数：

　　length：绘制控制图的长度(原始数据的长度)

　　" " "

　　total_length = len(self.value) # 目前已经计算出的统计量的长度

　　adjusted_length = int(math.floor(length* 2 **(-self.level))) #原始数据的长度对应到该层

　　小波系数的长度

　　plot_length = min(total_length, adjusted_length) #若目前计算出的统计量不足则绘制 已

　　计算的统计量

　　t_label = np.arange(2 ** self.level, plot_length* 2 ** self.level+1 , 2 ** self.level) #该层小波

　　系数对应到原始数据的时间戳

　　plt.plot(t_label, np.asarray(self.value[-plot_length:]) .reshape(plot_length,) , k-，linewidth = 0.6) #绘制统计量

　　plt.plot(t_label , self.ucl* np.ones(plot_length) , k--，linewidth = 0.6) #绘制上控制限

　　p lt.text(0 , self.ucl , $ /mathregular{U_{CL}} $ )

　　plt.plot(t_label , self.lcl* np.ones(plot_length) , k--，linewidth = 0.6) #绘制下控制限

　　p lt.text(0 , self.l cl , $ /mathregular{L_{CL}} $ )

　　plt.xlabel(时间(秒)，fontproperties="Sim Hei")

　　p lt.ylabel( z $ _i $ , rotation = 360 , fontsize = 16)

　　plt.show()

　　建立 EWMS类：

　　import math

　　import numpy as np

　　import matplotlib.pyplot as plt class EWMS() :

　　" " "

　　这是用于建立 EWMS控制图的类

　　" " "

　　GB/T 40859—202 1

　　def _init_(self, lambda_2 = 0.1 , mu_d = 0 , sigma_d = 1 , level = 1 , ucl = 1.32 , lcl = 0.67) :

　　" " "

　　初始化 EWMS类

　　参数：

　　lambda_2：平滑系数

　　mu_c： 均值目标值

　　sigma_c：方差目标值

　　level： 小波系数的分解层数

　　ucl： 上控制限参数(标准差的倍数)

　　lcl： 下控制限参数(标准差的倍数)

　　" " "

　　self.lambda_2 = lambda_2

　　self. mu_d = mu_d

　　self.sigma_d = sigma_d

　　self.level = level

　　self.value = [] #记录所有的统计量的值

　　self.old_number = sigma_d #用于表示 s Λ 2(i-1), 其中 s Λ 2(0) =sigma_d

　　self.new_number = 0 #用于表示 s Λ 2(i)

　　self.ucl = ucl* self.sigma_d** 0.5

　　self.lcl = lcl* self.sigma_d** 0.5

　　#上控制限

　　#下控制限

　　def update(self, data) :

　　" " "

　　获取新的小波系数后计算对应的控制图统计量

　　参数：

　　data： 新获取的小波系数

　　" " "

　　n = len(data) #新获取的小波系数的总长度

　　for i in range(0, n) :

　　self.new_number = self.lambda_2 * ( data[ i ]-self. mu_d ) ** 2 + ( 1-self. lambda_2 ) * self.

　　old_number #对应公式(8)

　　self.value.append(math.sqrt(self.new_number)) #计算统计量 s(i)

　　self.old_number = self.new_number

　　def plot(self, length) :

　　" " "

　　将最新获取的部分数据绘制成控制图

　　参数：

　　length：绘制控制图的长度(原始数据的长度)

　　" " "

　　total_length = len(self.value) # 目前已经计算出的统计量的长度

　　adjusted_length = int(math.floor(length* 2 **(-self.level))) #原始数据的长度对应到该层

　　小波系数的长度

　　GB/T 40859—202 1

　　plot_length = min(total_length, adjusted_length) #若目前计算出的统计量不足则绘制 已

　　计算的统计量

　　t_label = np.arange(2 ** self.level, plot_length* 2 ** self.level+1 , 2 ** self.level) #该层小波

　　系数对应到原始数据的时间戳

　　plt.plot(t_label, np.asarray(self.value[-plot_length:]) .reshape(plot_length,) , k-，linewidth = 0.6) #绘制统计量

　　plt.plot(t_label , self.ucl* np.ones(plot_length) , k--，linewidth = 0.6) #绘制上控制限

　　p lt.text(0 , self.ucl , $ /mathregular{U_{CL}} $ )

　　plt.plot(t_label , self.lcl* np.ones(plot_length) , k--，linewidth = 0.6) #绘制下控制限

　　p lt.text(0 , self.l cl , $ /mathregular{L_{CL}} $ )

　　plt.xlabel(时间(秒)，fontproperties="Sim Hei")

　　p lt.ylabel( s $ _i $ , rotation = 360 , fontsize = 16)

　　plt.show()

　　构建控制图：

　　假设构建控制图前获取了一组受控数据 data 。

　　import numpy as np

　　import matplotlib.pyplot as plt

　　from EWMA import EWMA #引用前面定义的 EWMA类

　　from EWMS import EWMS #引用前面定义的 EWMS类

　　import pywt #用于小波变换的程序包

　　# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = #这里应给定用于计算控制图参数目标值的平稳历史数据，赋值给 data变量

　　# data = np.array([data1, data2,…])

　　#也可以应用下行代码生成 1000 个仿真样本数据进行测试

　　data = np.random.normal(0,1,1000)

　　#通过受控数据计算控制图参数的目标值

　　sample_size = len(data)

　　#对应 7 . 4 步骤 1,将受控数据分解为各层小波系数，设置控制图参数并计算各参数的目标值J = 2 #分解层数为 2

　　wavelet_coefficient = pywt.wavedec(data, haar , level=J) #使用离散小波变换将受控数据分解为

　　各层小波系数

　　#设置 EWMA控制图的平滑系数，并计算近似层小波系数的均值、方差的目标值lamda_A2 = 0 . 1

　　mu_A2 = np.mean(wavelet_coefficient [0] , axis = 0)

　　sigma_A2 = np.cov(wavelet_coefficient [0] , rowvar = 0)

　　#设置第二层 EWMS控制图的平滑系数，并计算第二层细节层小波系数的均值、方差的目标值lamda_D2 = 0 . 1

　　mu_D2 = np.mean(wavelet_coefficient [1] , axis = 0)

　　sigma_D2 = np.cov(wavelet_coefficient[1] , rowvar = 0)

　　#设置第一层 EWMS控制图的平滑系数，并计算第一层细节层小波系数的均值、方差的目标值lamda_D1 = 0 . 1

　　GB/T 40859—202 1

　　mu_D1 = np.mean(wavelet_coefficient [2] , axis = 0)

　　sigma_D1 = np.cov(wavelet_coefficient [2] , rowvar = 0)

　　#查表 2 设置 EWMA控制限系数 L=2.715,创建 EWMA类

　　control_chart_A2 = EWMA(lamda_A2, mu_A2, sigma_A2,2,2.715)

　　control_chart_A2.update(wavelet_coefficient[0])

　　#查表 3 设置上控制限为 1 . 50 倍标准差、下控制限为 0 . 54 倍标准差，创建 EWMS类control_chart_D2 = EWMS(lamda_D2, mu_D2, sigma_D2,2,1.50,0.54)

　　control_chart_D2.update(wavelet_coefficient[1])

　　#查表 3 设置上控制限为 1 . 50 倍标准差、下控制限为 0 . 54 倍标准差，创建 EWMS类control_chart_D1 = EWMS(lamda_D1, mu_D1, sigma_D1,1,1.50,0.54)

　　control_chart_D1.update(wavelet_coefficient[2])

　　control_chart = [ control_chart_A2 , control_chart_D2 , control_chart_D1 ]

　　# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = #数据监测

　　#这里应给定用于监测的样本数据，赋值给 test_data变量# test_data = np.array([data1, data2,…])

　　#也可以用下行代码生成 1000 个仿真样本数据进行测试

　　test_data = np.random.normal(0,1,1000)

　　N = len(test_data) #用于监测的样本数据的长度

　　old_data = np.array([]) #保存用于计算小波系数的数据

　　i = 0 #计数为 0,该计数用于判断是否更新小波系数n = 0

　　while n < N: #每当获得一组新的数据时，进行以下操作：

　　new_data = test_data[ n]

　　i = i + 1

　　for j in range(1,J+1) : #从最低层开始看，是否可以更新系数

　　if i % 2 ** j = = 0 : #如果 i 整除 2 Λ j,说明第 j 层可以更新一个系数

　　temp = old_data[len(old_data) -2**j+1: len(old_data)]

　　temp = np.insert(temp , 2 ** j - 1 , values= new_data, axis = 0)

　　#这个系数由新数据以及前面的 2 Λ j- 1 个数据共同产生

　　new_coef = pywt.wavedec(temp , haar , level=j) #获得新的系数

　　detail_new_coef = new_coef [1] #只取第 j 层细节层系数

　　wavelet_coefficient [J-j+1] = np.r_[wavelet_coefficient[J-j+1], detail_new_coef]

　　#更新该层系数

　　control_chart[J-j+1]. update(detail_new_coef) #将新的系数纳入对应的 EWMS 类

　　中，更新统计量

　　if j = = J: #如果为最高分解层数，那么也要更新近似层

　　# 以下步骤与上面类似

　　approximate_new_coef = new_coef [0]

　　wavelet_coefficient [0] = np.r_[wavelet_coefficient [0], approximate_new_coef] control_chart[0].update(approximate_new_coef)

　　old_data = np.hstack((old_data, new_data)) #将新数据纳入老数据中，重新获取数据if i = = 2 ** J : #如果计数达到了 2 Λ J,则归 0

　　GB/T 40859—202 1

　　i = 0

　　old_data = np.array([]) #清空用于计算小波系数的数据

　　n = n + 1

　　for i in range(0,J+1) :

　　plt.rcParams[ figure.dpi ] = 300

　　control_chart[i].plot(N) #分别绘制各个控制图

　　D.3 多变量流式数据控制图构建算法

　　建立 MEWMA类：

　　import math

　　import numpy as np

　　import matplotlib.pyplot as plt class MEWMA() :

　　" " "

　　这是用于建立 MEWMA控制图的类

　　" " "

　　def _init_(self, lambda_3 = 0.1 , mu_c = 0 , sigma_c = 1 , level = 1 , ucl = 8.64) :

　　" " "

　　初始化 MEWMA类

　　参数：

　　lambda_3：平滑系数

　　mu_c： 均值向量目标值

　　sigma_c：协方差矩阵目标值

　　level： 小波系数的分解层数

　　ucl： 控制限

　　" " "

　　self.lambda_3 = lambda_3

　　self. mu_c = mu_c

　　self.sigma_c = np.mat(sigma_c)

　　self.level = level

　　self.value = [] #记录所有的统计量的值

　　self.old_number = mu_c #用于表示 Z(i-1), 其中 Z(0) = mu_c

　　self.new_number = np.array([]) #Z(i)

　　self.ucl = ucl #控制限

　　def update(self, data) :

　　" " "

　　获取新的小波系数后计算对应的控制图统计量

　　参数：

　　data：新获取的小波系数向量

　　" " "

　　[n, p] = data.shape # n 为小波系数的总长度，p 为变量维数

　　GB/T 40859—202 1

　　for i in range(0, n) :

　　self.new_number = self.lambda_3 * data[i , :] + (1-self.lambda_3) * self.old_number #

　　对应公式(11)

　　sigma_z = self.lambda_3/(2-self.lambda_3) * self.sigma_c #对应公式 (13)

　　trans = self.new_number-self. mu_c

　　self.value.append(np.float(np.dot(np.dot(trans, np.linalg.inv(sigma_z)) , trans.T))) #

　　对应公式(12)计算统计量并保存

　　self.old_number = self.new_number

　　def plot(self, length) :

　　" " "

　　将最新获取的部分数据绘制成控制图

　　参数：

　　length：绘制控制图的长度(原始数据的长度)

　　" " "

　　total_length = len(self.value) # 目前已经计算出的统计量的长度

　　adjusted_length = int(math.floor(length* 2 **(-self.level))) #原始数据的长度对应到该层

　　小波系数的长度

　　plot_length = min(total_length, adjusted_length) #若目前计算出的统计量不足则绘制 已

　　计算的统计量

　　t_label = np.arange(2 ** self.level, plot_length* 2 ** self.level+1 , 2 ** self.level) #该层小波

　　系数对应到原始数据的时间戳

　　plt.plot(t_label, np.asarray(self.value[-plot_length:]) .reshape(plot_length,) , k-，linewidth = 0.6) #绘制统计量

　　plt.plot(t_label , self.ucl* np.ones(plot_length) , k--，linewidth = 0.6) #绘制上控制限

　　p lt.text(0 , self.ucl , $ /mathregular{U_{CL}} $ )

　　plt.xlabel(时间(秒)，fontproperties="Sim Hei")

　　p lt.ylabel( $ /mathregular{Y_i Λ 2} $ , rotation = 360 , fontsize = 16)

　　plt.show()

　　建立 MEWMC类：

　　import math

　　import numpy as np

　　import matplotlib.pyplot as plt class MEWMC() :

　　" " "

　　这是用于建立 MEWMC控制图的类

　　" " "

　　def _init_(self, lambda_4 = 0.1 , mu_d = 0 , sigma_d = 1 , level = 1 , ucl = 1.169) :

　　" " "

　　初始化 MEWMC类

　　参数：

　　lambda_4：平滑系数

　　mu_d： 均值向量目标值

　　GB/T 40859—202 1

　　sigma_d：协方差矩阵目标值

　　level： 小波系数的分解层数

　　ucl： 控制限

　　" " "

　　self.lambda_4 = lambda_4

　　self.mu_d = np.mat(mu_d)

　　self.sigma_d = sigma_d

　　self.level = level

　　self.value = [] #记录所有的统计量的值

　　self.old_number = np.eye(len(mu_d)) #用于表示 S(i-1), 其中 S(0) = I(p)

　　self.new_number = np.array([]) #用于表示 S(i)

　　self.ucl = ucl #控制限

　　eigVals, eigVects=np.linalg.eig(np.mat(sigma_d))

　　for i in range(0, len(eigVals)) :

　　eigVects[ : , i] = eigVects[ : , i] * (1/np.sqrt(eigVals[i]))

　　self.Amatrix = eigVects.T #计算公式(14) 中的 A 矩阵

　　def update(self, data) :

　　" " "

　　获取新的小波系数后计算对应的控制图统计量

　　参数：

　　data：新获取的小波系数向量

　　" " "

　　[n,p] = data.shape # n 为小波系数的总长度，p 为变量维数

　　for i in range(0, n) :

　　trans = data[ i , :] -self. mu_d

　　U = np.dot(self.Amatrix, trans.T)

　　UUT = np.dot(U, U.T) #对应公式(15) 中的 U(i) * U(i) Λ T

　　self.new_number = self.lambda_4 * UUT + ( 1-self. lambda_4) * self.old_number #对

　　应公式(15)

　　self.value.append(np.trace(self.new_number)-np.log(np.linalg.det(self.new_number))- p) #对应公式(16)计算统计量并保存

　　self.old_number = self.new_number

　　def plot(self, length) :

　　" " "

　　将最新获取的部分数据绘制成控制图

　　参数：

　　length：绘制控制图的长度(原始数据的长度)

　　" " "

　　total_length = len(self.value) # 目前已经计算出的统计量的长度

　　adjusted_length = int(math.floor(length* 2 **(-self.level))) #原始数据的长度对应到该层

　　小波系数的长度

　　GB/T 40859—202 1

　　plot_length = min(total_length, adjusted_length) #若目前计算出的统计量不足则绘制 已

　　计算的统计量

　　t_label = np.arange(2 ** self.level, plot_length* 2 ** self.level+ 1 , 2 ** self.level) #该层小波

　　系数对应到原始数据的时间戳

　　plt.plot(t_label, np.asarray(self.value[-plot_length:]) .reshape(plot_length,) , k-，linewidth = 0.6) #绘制统计量

　　plt.plot(t_label , self.ucl* np.ones(plot_length) , k--，linewidth = 0.6) #绘制上控制限

　　p lt.text(0 , self.ucl , $ /mathregular{U_{CL}} $ )

　　plt.xlabel(时间(秒)，fontproperties="Sim Hei")

　　p lt.ylabel( $ /mathregular{y_i} $ , rotation = 360 , fontsize = 16)

　　plt.show()

　　构建控制图：

　　假设构建控制图前获取了一组 p= 2 维的受控数据，分别记为 a 和 b。

　　import numpy as np

　　import matplotlib.pyplot as plt

　　from MEWMA import MEWMA #引用前面定义的 MEWMA类

　　from MEWMC import MEWMC #引用前面定义的 MEWMC类

　　import pywt #用于小波变换的程序包

　　# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = #这里应给定用于计算控制图参数目标值的平稳历史数据，赋值给 data变量

　　# a = np.array([a1, a2,…])

　　# b = np.array([b1, b2,…]) # data = np.vstack((a, b)) .T

　　#也可以应用下行代码生成 1000 个仿真样本数据进行测试

　　# a = np.random.normal(0,1,1000)

　　# b = np.random.normal(0,1,1000) # data = np.vstack((a, b)) .T

　　#通过受控数据计算控制图参数的目标值

　　sample_size = data.shape[0] #受控数据长度p = data.shape[1] #变量维度

　　#对应 8 . 5 步骤 1,设置控制图参数并计算各参数的目标值J = 2 #分解层数为 2

　　wavelet_coefficient_a = pywt.wavedec(a, haar , level=J)

　　wavelet_coefficient_b = pywt.wavedec(b , haar , level=J)

　　#对应 8 . 2,分别对各个变量进行小波变换，并将小波系数重组为系数向量wavelet_coefficient = []

　　for i in range(0, len(wavelet_coefficient_a)) : #将对应层系数组成矩阵

　　temp = np.vstack((wavelet_coefficient_a[i], wavelet_coefficient_b[i]))

　　temp = temp.T

　　wavelet_coefficient.append(temp)

　　#设置 MEWMA控制图的平滑系数，并计算近似层小波系数向量的均值向量、协方差矩阵的目标值lamda_A2 = 0 . 1

　　GB/T 40859—202 1

　　mu_A2 = np.mean(wavelet_coefficient[0] , axis = 0)

　　sigma_A2 = np.cov(wavelet_coefficient[0] , rowvar = 0)

　　#设置第二层 MEWMC控制图的平滑系数，并计算第二层细节层小波系数向量的均值向量、协方差矩阵的目标值

　　lamda_D2 = 0 . 1

　　mu_D2 = np.mean(wavelet_coefficient[1] , axis = 0)

　　sigma_D2 = np.cov(wavelet_coefficient[1] , rowvar = 0)

　　#设置第一层 MEWMC控制图的平滑系数，并计算第二层细节层小波系数向量的均值向量、协方差矩阵的目标值

　　lamda_D1 = 0 . 1

　　mu_D1 = np.mean(wavelet_coefficient[2] , axis = 0)

　　sigma_D1 = np.cov(wavelet_coefficient[2] , rowvar = 0)

　　#查表 4 设置 MEWMA控制限控制限 h_A=10.08,创建 MEWMA类control_chart_A2 = MEWMA(lamda_A2, mu_A2, sigma_A2,2,10.08)

　　control_chart_A2.update(wavelet_coefficient[0])

　　#查表 5 设置 MEWMC控制限控制限 h_C=0.63,创建 MEWMC类control_chart_D2 = MEWMC(lamda_D2, mu_D2, sigma_D2,2,0.63)

　　control_chart_D2.update(wavelet_coefficient[1])

　　#查表 5 设置 MEWMC控制限控制限 h_C=0.63,创建 MEWMC类control_chart_D1 = MEWMC(lamda_D1, mu_D1, sigma_D1,1,0.63)

　　control_chart_D1.update(wavelet_coefficient[2])

　　control_chart = [ control_chart_A2 , control_chart_D2 , control_chart_D1 ]

　　# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = #数据监测

　　#这里应给定用于控制图监测的样本数据，赋值给 test_data变量

　　# a = np.array([a1, a2,…])

　　# b = np.array([b1, b2,…])

　　# test_data = np.vstack((a, b)) .T

　　#也可以应用下面三行代码仿真生成 1000 个二维样本数据进行测试

　　# a = np.random.normal(0,1,1000)

　　# b = np.random.normal(0,1,1000) # test_data = np.vstack((a, b)) .T

　　N = test_data.shape[0] #用于监测的样本数据的长度

　　old_data = np.array([]) #保存用于计算小波系数的数据i = 0 #计数为 0,该计数用于判断是否更新小波系数

　　n = 0

　　while n < N: #每当获得一组新的数据时，进行以下操作： new_data = test_data[ n]

　　i = i + 1

　　for j in range(1,J+1) : #从最低层开始看，是否可以更新系数

　　if i % 2 ** j = = 0 : #如果 i 整除 2 Λ j,说明第 j 层可以更新一个系数

　　temp = old_data[len(old_data) -2** j+1: len(old_data) , :]

　　temp = np.insert(temp , 2 ** j - 1 , values= new_data, axis = 0)

　　GB/T 40859—202 1

　　#这个系数由新数据以及前面的 2 Λ j- 1 个数据共同产生

　　new_coefficient_a = pywt.wavedec(temp[ : ,0] , haar , level=j)

　　new_coefficient_b = pywt.wavedec(temp[ : ,1] , haar , level=j)

　　#两个变量分别计算这层小波系数，每个变量都得到一个新系数

　　new_coefficient = []

　　for k in range(0, len(new_coefficient_a)) : #对应层重组系数矩阵

　　temp1 = np.vstack((new_coefficient_a[k], new_coefficient_b[k]))

　　temp1 = temp1 . T

　　new_coefficient.append(temp1)

　　detail_new_coefficient = new_coefficient[1] #只取第 j 层细节层系数

　　wavelet_coefficient[J-j+1] = np.r_[wavelet_coefficient[J-j+1] , detail_new_coeffi- cient] #更新该层系数

　　control_chart[J-j+1]. update(detail_new_coefficient) # 将新的系数纳 入 对 应 的

　　MEWMC类中，更新模型

　　if j = = J: #如果达到了最高分解层数，那么也要更新近似层

　　# 以下步骤与上面类似

　　approximate_new_coefficient = new_coefficient[0]

　　wavelet_coefficient[0] = np.r_[wavelet_coefficient[0],approximate_new_coefficient] control_chart[0].update(approximate_new_coefficient)

　　old_data = np.vstack((old_data.reshape( -1,2) , new_data)) #将新数据纳入老数据中，重新获

　　取数据

　　if i = = 2 ** J : #如果计数达到了 2 Λ J,则归 0

　　i = 0

　　old_data = np.array([]) #清空用于计算小波系数的数据

　　n = n + 1

　　for i in range(0,J+1) :

　　plt.rcParams[ figure.dpi ] = 300

　　control_chart[i].plot(N) #分别绘制各个控制图

　　GB/T 40859—202 1

　　参 考 文 献

　　[1] GB/T 17989 . 1—2020 控制图 第 1 部分：通用指南

　　[2] ISO 7870-6 Control charts—Part 6 : EWMA control charts

　　[3] Montgomery D.C. Introduction to Statistical Quality Control. John Wiley & Sons , 2012.

　　[4] Daubechies I.Ten lectures on wavelets. Society for Industrial and Applied Mathematics, 1992.

　　[5] Mallat S.G. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1989 , 11(7) pp.674-693.

　　[6] Ganesan R. , Das T.K. , Venkataraman V. Wavelet-based multiscale statistical process mo- nitoring : A literature review.IIE Transactions. 2004,36(9) pp.787-806.

　　[7] MacGregor J.F. , Harris T.J. The exponentially weighted moving variance. Journal of Quali- ty Technology. 1993 , 25 pp.106-118.

　　[8] Lowry C.A. , Woodall W. H. , Champ C.W. , and Rigon S.E. A Multivariate Exponentially Weighted Moving Average Control Chart. Technometrics. 1992,34(1) pp.46-53.

　　[9] Hawkins , D.M. , Maboudou-Tchao , E.M. Multivariate Exponentially Weighted Moving Co- variance Matrix. Technometrics. 2008,50(2) pp.155-166.