(高清版)GBT 40859-2021 流式数据监测控制图

流式数据监测控制图2021-10-11发布2022-05-01实施国家标准化管理委员会I前言 2规范性引用文件 13术语和定义 14符号 35应用条件 46基于离散小波变换的数据预处理 7单变量流式数据控制图的构建 67.1概述 67.2单变量流式数据EWMA控制图的构建 7.3单变量流式数据EWMS控制图的构建 67.4单变量流式数据控制图应用步骤 78多变量流式数据控制图的构建 78.1概述 8.2小波系数重组 8.3多变量流式数据MEWMA控制图的构建 88.4多变量流式数据MEWMC控制图的构建 98.5多变量流式数据控制图应用步骤 9附录A(资料性)离散小波变换及多分辨率分析原理 附录B(资料性)单变量流式数据应用实例 附录C(资料性)多变量流式数据应用实例 附录D(资料性)代码实现 参考文献 Ⅲ本文件按照GB/T1.1—2020《标准化工作导则第1部分：标准化文件的结构和起草规则》的规定起草。请注意本文件的某些内容可能涉及专利。本文件的发布机构不承担识别专利的责任。本文件由全国统计方法应用标准化技术委员会(SAC/TC21)提出并归口。本文件起草单位：清华大学、中国标准化研究院。在现代复杂生产系统和装备运行过程中，由于自动化传感器的广泛应用，所获取的数据往往是流式数据的形式。所谓流式数据，是一组数量庞大、顺次有序、快速高频到达的数据序列。现有的控制图往往基于较大间隔进行离散采样获得的数据进行过程控制。由于流式数据具有快速高频等特点，现有控制图通常无法适用于控制产生流式数据的过程。传统的数据采集与分析往往是低频的，如每小时对一条生产水杯的生产线进行产品采样和关键质量指标测量，以判断生产过程的状态是否正常。而在复杂系统的生产制造和运行过程中，通常会布局很多设备状态传感器和自动化检测仪器，以秒甚至更小的时间单位为间隔，对设备运行状态及产品质量指标进行实时、自动和持续的数据收集，由此产生流式数据。这些流式数据包含比低频数据更加复杂、全面的信息，但其信息密度也相对较低，因此，需要合适的特征提取手段，挖掘流式数据中包含的信息。20世纪，有学者提出了多分辨率分析方法。这是一种基于离散小波变换的数据分析方法，可以将数据分解到不同的频率尺度下进行分析。此后，小波分析的理论方法不断完善。多分辨率小波分析十分适用于处理高频的流式数据。将高频的流式数据分解到不同的频率尺度后，可在不同的尺度下分别分析，提取流式数据中复杂的信息。本文件可对高速铁路、装备制造、复杂产品制造等系统在运行过程中所产生的实时状态数据和产品质量数据进行监测，对关键信号的异常变化进行预警，有可能在早期发现危险故障，提升系统安全可靠1流式数据监测控制图本文件给出了对单变量及多变量流式计量型数据进行监测的控制图方法。本文件适用于对产生流式数据的一个或多个有一定关联关系的过程或质量变量，在均值和基本波动特征已知的条件下，对偏移和波动进行控制。2规范性引用文件下列文件中的内容通过文中的规范性引用而构成本文件必不可少的条款。其中，注日期的引用文件，仅该日期对应的版本适用于本文件；不注日期的引用文件，其最新版本(包括所有的修改单)适用于本文件。GB/T3358.1统计学词汇及符号第1部分：一般统计术语与用于概率的术语GB/T3358.2统计学词汇及符号第2部分：应用统计3术语和定义GB/T3358.1和GB/T3358.2界定的以及下列术语和定义适用于本文件。一组数量庞大、顺次有序、快速高频到达的数据序列。注1:传感器监测数据、交通监测数据等往往都是以流式数据的形式被收集。注2:流式数据通常是由一个或多个数据源持续产生，随时间延续而增长的动态数据集合，具有采样间隔小、米样频指数加权移动平均控制图exponentiallyweightedmovingaveragecontrolchart用指数平滑加权平均评估和监测过程水平的计量控制图。[来源：GB/T3358.2—2009,2.3.16,有修改]指数加权均方误差控制图exponentiallyweightedmeansquarederrorcontrolchart用指数平滑加权平均评估和监测过程偏移与变异的计量控制图。MEWMA控制图MEWMAchart多变量指数加权移动平均控制图multivariateexponentiallyweightedmovingaveragecontrolchart将两个或两个以上相关的变量合成一个样本统计量，并用指数平滑加权平均评估和监测过程水平的计量控制图。2多变量指数加权移动协方差矩阵控制图multivariateexponentiallyweightedmovingcovariance将两个或两个以上相关的变量合成一个样本统计量，并用指数平滑加权平均评估和监测过程变异的计量控制图。小波变换wavelettransform用样本数据与小波函数的内积计算得到小波系数的过程。离散小波变换discretewavelettransform用样本数据与正交小波基的内积计算得到小波系数的过程。小波基函数waveletbasisfunction均方可积且均值为0的函数。注2:小波基函数通常表示为ψ(t)。用小波基函数伸缩与平移得到的一组函数。注：小波函数可表示为ψ(t)=a-i2ψ(a-it-kb),其中a为伸缩参数，b为平移参数，j为尺度参数，k为位置参数。正交小波基orthonormalwaveletbases用某些小波基函数进行特定的离散化伸缩与平移，得到的可作为正交基的小波函数。注1:Haar小波是可用于产生正交小波基的一种小波基函数。注2:通常选取伸缩参数a=2,平移参数b=1。小波系数waveletcoefficient用样本数据与小波函数内积的结果。注：内积过程可表示为d=Jf(t)ψ;(t)e-²dt,其中f(t)为样本数据。可生成正交基的范数为1的规范化函数，与母小波对应。注1:并不是所有的母小波都有对应的父小波。注2:可用于离散小波变换的母小波都有对应的父小波。注3:母小波用于表现样本波动信息，父小波用于表现样本近似信息。多分辨率分析multiresolutionanalysis用离散小波变换得到的小波系数进行数据分析的方法。注：离散小波变换得到的小波系数包含不同频率尺度下的信息。3用样本数据与父小波内积的结果。注1:近似层小波系数包含样本趋势信息。注2:样本经离散小波变换通常仅保留一层近似层小波系数。细节层小波系数detailwaveletcoefficient用样本数据与母小波内积的结果。注1:细节层小波系数包含样本波动信息。注2:样本经离散小波变换会产生若干层近似层小波系数。4符号下列符号适用于本文件。ARL₀平均链长ARL过程受控状态下的平均链长CJ.k第J层第k个近似层小波系数CJ.k第J层第k个近似层小波系数向量d;.k第j层第k个细节层小波系数dj,k第j层第k个细节层小波系数向量hAMEWMA控制图控制限hcMEWMC控制图控制限I,p维的单位矩阵J离散小波变换分解层数LEWMA控制图控制限系数p多变量数据的维度Sj.i第j个EWMS控制图统计量Lc下控制限Uc上控制限X多变量流式数据x单变量流式数据第i维变量Y?MEWMA控制图统计量yj,i第j个MEWMC控制图统计量≥iEWMA控制图统计量λ₁EWMA控制图平滑系数λ₂EWMS控制图平滑系数λ₃MEWMA控制图平滑系数λ₄MEWMC控制图平滑系数μc近似层小波系数均值的目标值μc近似层小波系数向量均值的目标值μdj第j层细节层小波系数均值的目标值μdj第j层细节层小波系数向量均值的目标值4σ.近似层小波系数标准差的目标值Ga;第j层细节层小波系数标准差的目标值Ze近似层小波系数向量均值协方差矩阵的目标值第j层细节层小波系数向量均值协方差矩阵的目标值5应用条件现实的生产过程中，由于各种传感器的大量使用，过程样本数据往往以流式数据的形式被收集。流式数据具有采样频率高、采样间隔短、数据量大等特点。由于数据获取频率高，短时间内便能获得大量数据，提高了数据计算和存储的难度，故流式数据需采用新的监测方法。此外，生产过程中的过程信号往往十分复杂，可能同时存在着多种不同的特征信号。这些复杂的信号组合在一起会导致难以对故障的原因进行诊断，更严重的是有时多种特征信号的叠加可能会导致一些特征信号难以被监测。多种特征信号叠加的数据如图1所示。说明：xA——白噪声；xg——包含均值突变的数据；xc——包含波动突变的数据；xp——xA、xB、xc叠加后的数据。图1复杂的特征信号从图1中可看出，xg发生了均值的偏移，xc发生了波动幅度的变化。虽然这二者变化特征明显，但是对于由xA、xg、xc叠加得到的数据xp,其均值偏移与波动变化两种特征信号互相影响。若使用单一控制图监测xp,则对均值偏移与波动变化都难以进行有效识别。若可以将图1中叠加在一起的复杂信号分开进行监测，那么便会大大提高监测的效果。离散小波变换是一种对数据进行分解的方法，原始数据经过离散小波变换后，可在不同的频率尺度上监测不同的特征信号。此外，即使原始数据存在一定自相关性，经离散小波变换得到的小波系数的自相关性会显著5GB/T40859—2021降低或消失，亦即离散小波变换所得到的系数对原始数据中存在的自相关性不敏感，因此，离散小波变换的方法可适用于流式数据的监测中。在使用本控制图对过程进行监测时，应有一定量能够代表过程正常运行状态的历史数据，用于计算控制图参数的目标值。6基于离散小波变换的数据预处理为了实现对流式数据样本的监测，需要首先对所获取数据应用基于Haar小波的离散小波变换算法进行预处理，生成小波系数，然后将小波系数用于控制图监测。Haar小波是一种常用的正交小波，其定义为： (1)结合附录A中介绍的小波函数的公式(A.1),可知其在第j层的小波函数的定义域长度仅为2',因此仅使用有限的数据点即可计算一个小波系数。假定分解层数为J,流式数据序列为x={x₁,xz,…}。近似层小波系数和细节层小波系数的计算公式分别为：在监测过程中，需对新获取的数据点进行计数。假定计数为i,当i可以整除2′时，计算前j层小波系数。表1给出了应用小波变换得到细节层小波系数的示例。表中每一列表示一个新获取的数据，每一行表示小波分解的层，表内的d;..表示第j层的第k个小波系数。表中共展示了8个数据点，进行3层小波分解的系数。表中可看出，当计数为2时，可通过1、2两个数据点计算获得一个第一层的小波系数d.1;计数为8时，可通过1～8八个数据点计算第三层小波系数ds.1,同时通过5、6、7、8四个数据点计算第二层小波系数d₂.2,通过7、8两个数据点计算第一层小波系数d₁.4。当获取新的数据时，可重复该过程，不断获取新的小波系数。表1小波变换样例层数计数123456781d₁,i2d₂.i3通过该算法，可将实时获得的数据及时分解为各层小波系数，共可得到1层近似层系数{cj.,cj.z,…}以及J层细节层小波系数{d;.1,d;₂,…},j∈{1,…,J},用于对过程的监测。在构建控制图的过程中，分解层数J的选择取决于流式数据特征信号的频率范围，若要监测较低频的波动信号则选取较大的分解层数J。但分解层数J过大会导致计算复杂度增加，影响监测效率。通常分解层数J建议选为在应用控制图时，往往假设被监测特征服从正态分布。在流式数据服从正态分布的假设下，基于Haar小波的数据预处理得到的各层小波系数仍然满足正态分布的假设，可使用控制图进行监测。在过6GB/T40859—2021程受控状态下获取的数据经预处理得到的各层小波系数，可用于计算控制图参数的目标值。7单变量流式数据控制图的构建7.1概述流式数据经小波变换分解为小波系数后，使用EWMA控制图监测近似层小波系数，使用EWMS控制图监测各层细节层小波系数。附录B中给出了构建控制图的案例。7.2单变量流式数据EWMA控制图的构建EWMA控制图适用于监测数据的均值变化。应使用经离散小波变换得到的近似层小波系数{cj.1,cj.2,…}构建该控制图。EWMA统计量可表示为：将公式(4)展开，可以得到²;=λ₁cjx+λ₁(1-λ₁)cjx-i+λ₁(1-λ₁)²cjx-2+ (5)可以看出，EWMA控制图将所有以往的信息都包括进来计算统计量。距离当前时间越远，该数据点的信息在此时刻的统计量中的占比越小。EWMA控制图的控制限为： (7)其中，L为控制限系数。应用ARL分析的方法可以确定控制图的两个参数：控制限系数L和平滑系数λ₁。通常，λ₁的取值取决于偏移量的大小，一般推荐范围为0.05≤λ₁≤0.5。ARL=370时的参数选择如表2所示。表2ARL₀=370时EWMA控制图参数选择标准L7.3单变量流式数据EWMS控制图的构建EWMS控制图适用于监测数据的方差变化。应使用经离散小波变换得到的各层细节层小波系数{dj.i,d;₂,…},j∈{1,…,J},共构建J个EWMS控制图。其中第j个控制图的EWMS统计量可表示为：s;;=λ₂(d;;-pa,)²+(1-λ₂)s?.;- (8)其中λ₂∈(0,1)为平滑系数，s?.o=oa,。通常，绘制关于s;的控制图，而非关于s?.的控制图，即将迭代序列的平方根绘制入控制图。EWMS控制图的控制限为：7GB/T40859—2021…其中，α表示发生第一类错误的概率，Xi-₂(v)和X2(v)分别表示自由度为v的卡方分布的(1-α/2)×100%和(a/2)×100%分位数，v=(2-λ₂)/λ₂。λ₂的取值取决于偏移量的大小，一般推荐范围为0.01≤λ₂≤0.4。ARL₀=370时(α=0.0027),平滑系数λ₂与控制限的对应关系如表3所示。表3ARL₀=370时EWMS控制限选择标准7.4单变量流式数据控制图应用步骤单变量流式数据控制图构建的具体步骤如下。步骤1确定控制图参数：选取离散小波变换分解层数J、EWMA控制图的控制限L和平滑系数λ₁、EWMS控制图的第一类错误的概率α和平滑系数λ₂,并计算各控制图的控制限。步骤2计算小波系数：对新获取的数据x;进行计数，代入公式(2)、公式(3)获取新的近似层或细节层小波系数。步骤3计算统计量：使用新获得的小波系数，近似层小波系数代入公式(4)计算EWMA统计量，各层细节层小波系数分别代入公式(8)计算各层的EWMS统计量。步骤4绘制控制图：将各统计量及相应的控制限分别绘制在控制图上，共得到一张EWMA控制图及J张EWMS控制图。根据控制图判断受控状态。当任何一张控制图出现控制点超过控制限时，表明过程可能“失控”,需要对过程采取必要的措施，予以纠正。8多变量流式数据控制图的构建8.1概述假设某一时刻开始获得多变量流式数据序列X={x₁,x₂,…},变量的维度为p,即x₁=为达到监测目的，首先使用离散小波变换将每一个变量的数据分解为近似层小波系数和细节层小波系数，再将不同变量同一层的小波系数重组，构成多维的小波系数向量。使用MEWMA控制图监测近似层小波系数向量，使用MEWMC控制图监测各细节层小波系数向量。附录C中给出了构建控制图的案例。8.2小波系数重组使用离散小波变换将每一维数据都分解至J层。以2维数据为例：第1维数据{x(),x¹),…}可分解得到1层近似层系数{cy,c2,…}以及J层细节层小波系数{dh,d2,…},j∈{1,…,J};第2维数据{x{²),x₂),…}可分解得到1层近似层系数{cY:,c2,…}以及J层细节层小波系数{d},d}:2,…},j∈{1,…,J}。将2维数据的近似层系数组成向量，可以得到1层近似层系数向量{cj.1,cj.₂,…},其中cj.=(c},c{});将2维数据对应的细节层系数组成向量可得到J层细节层小波系数{d;.,d;.z,…},j∈{1,…,J},其中d;x=(d?,d{)。图2展示了变量维度为p时小波系数的重组过程。8GB/T40859—2021x(1)x(2}小波变换dg系数向量cyMEWMA控制图MEWMC控制图**MEWMC控制图细节层小波系数向量D;细节层小波系数向量D₁小波变换小波变换x(u)图2多变量流式数据控制图小波系数重组示意图8.3多变量流式数据MEWMA控制图的构建MEWMA控制图适用于监测多变量流式数据的均值变化。应使用8.2中经离散小波变换并重组得到的近似层小波系数向量{cj.,cj.z,…}构建控制图。定义向量序列Z;,Z₀=u.。每获取一个新的小波系数向量cj.,便可通过迭代计算得到Z;:Z;=λ₃cj.x+(1-λ₃)Z;-i……………(11)其中，常数λ₃∈(0,1)是平滑系数。判断数据失控的条件为：Y}=(Z;-μe)¹Ez¹(Z₁-μe)>hA…………(12)其中，Y?为该控制图的统计量，hA为MEWMA控制图的控制限；Z,为Z;的协方差矩阵，其可表示为：应用ARL分析的方法可以确定控制图的两个参数：控制限hA和平滑系数λ₃。λ₃的取值取决于偏移量的大小，一般推荐范围为0.05≤λ₃≤0.5。应选择合适的hʌ使得控制图可达到一定的受控ARL,ARL₀=370时hA选择的原则如表4所示。表4ARL₀=370时MEWMA控制图控制限选择标准p246898.4多变量流式数据MEWMC控制图的构建MEWMC控制图适用于监测多变量流式数据的方差及协方差变化。应使用8.2中经离散小波变换并重组得到的各层细节层小波系数向量{d,i,d;2,…},j∈{1,…,J}构建第j个MEWMC控制图时，定义向量序列S;i,S;.o=I,。构建控制图。每获取一个新的小波系数向量d,.,先进行一步转换： (14)S;x=λ₁U;U';x+(1-λ₄)S;x-1 (15)y,.=tr(S,)-ln|S;|-p>hc…………(16)迹的函数。构建该控制图同样需要确定两个参数：控制限hc和权重系数λ。λ₄的取值取决于偏移量的大小，一般推荐范围为0.01≤λ₄≤0.4。应选择合适的hc使得控制图可达到一定的受控ARL,ARLo=370时hc选择的原则如表5所示。表5ARL₀=370时MEWMC控制图控制限选择标准p24688.5多变量流式数据控制图应用步骤多变量流式数据控制图构建的具体步骤如下：步骤1确定控制图参数：选取离散小波变换分解层数J、MEWMA控制图的平滑系数λ₃、MEWMC控制图的平滑系数λ₄,并选择各控制图的控制限。步骤2计算小波系数向量：对新获取的数据x;进行计数，并分别代入公式(2)、公式(3)获取新的近似层或细节层小波系数，再将对应层小波系数重组为小波系数向量。步骤3计算统计量：使用新获得的小波系数向量，近似层小波系数向量代入公式(11)、公式(12)和公式(13)计算MEWMA统计量，各层细节层小波系数向量分别代入公式(14)、公式(15)和公式(16)计算各层的MEWMC统计量。步骤4绘制控制图：将各统计量及相应的控制限分别绘制在控制图上，共得到一张MEWMA控制图及J张MEWMC控制图。根据控制图判断受控状态。当任何一张控制图出现控制点超过控制限以上步骤往往需要使用计算机辅助实现，实现控制图的Python代码在附录D中列出。(资料性)离散小波变换及多分辨率分析原理多分辨率数据分析的原理是利用离散小波变换将数据分解到不同的尺度上进行分析，这对于流式数据的监测十分有效。离散小波变换是一种时频分析算法，与传统的傅里叶变换不同，它能保留原始数据在时域上的信息，即可将原始数据分解到不同的频率尺度上。用于过程监测的流式数据通常可通过小波函数分解成近似信号和细节信号的形式：为近似小波系数和细节小波系数。其表达式为：当k=m2i,m∈Z时，这一小波变换的方法被称作离散小波变换。在分解时，首先确定分解层数J,之后依次计算各层的小波系数，在每层中以k=i2',i∈Z逐步平移小波基函数。共可得到1层的近似小波系数和J层细节小波系数，如图A.1所示。细节小波系数细书小波系数原始数据第二层第一层近似小波系数近似小波系数第三层近似小波系数细节小波系数图A.1离散小波变换小波系数每一个小波系数都是由一段原始数据与该层小波基函数的内积得到，即该小波系数可以描述原始数据在一段时间内与该层小波基函数的相似程度。因此不同层的小波系数描述了原始信号在不同频域X2di21tlz2x——原始数据；d₁——第一细节层小波系数；d₂——第二细节层小波系数；c₂——近似层小波系数。从图A.2中可看出，原始信号x包含有均值偏移和方差变化两种不同的特征信号。利用离散小波变换将其分解到2层后，可看出均值偏移的信息表现在近似小波系数上，而波动变化的信息表现在细节系数上。因此，可对分解后的各层小波系数分别进行监测，这便是多分辨率分析的原理。(资料性)单变量流式数据应用实例在某直拉单晶硅生产过程中，需要对单晶的关键质量参数进行监测。记该质量参数为变量x。单晶拉制过程中，数据的获取间隔约1s。生产系统已经平稳运行了一段时间，获取了足够长的平稳数据X。本案例中相关数据已进行脱敏处理。在建立控制图前首先要基于历史平稳数据对控制图参数的目标值进行估算。当所获得历史数据存在异常值时，基于工程经验判断和统计分析，对异常值进行剔除，仅使用平稳数据进行参数的目标值估计。将过程平稳状态下获取的平稳数据X应用离散小波变换分解到2层，得到一层近似层小波系数{c₂.1,c2.2,…}和两层细节层小波系数{d₁1,d₁.2,…}和{d₂.,d₂.z,…}。根据这些小波系数可估计出各层小波系数的均值和方差的目标值，计算结果如表B.1所示。表B.1各层小波系数均值与方差小波系数均值方差0.000.14970.000.00380.000.0075现从该过程获取新的1000组数据x={x₁,…,x₁am},控制图监测。如图B.1所示。下面对这组新的数据实施时间/s流式数据x对数据x={x₁,…,x₁o}应用公式(2)、公式(3),并代入j=2,即可使用离散小波变换将数据分解到2层，同样得到一层近似层小波系数{c₂.1,c2.2,…,c2.250}和两层细节层小波系数{di.i,d₁.2,…,d1.sm}对于近似层小波系数，使用EWMA控制图监测变量的均值，选择平滑系数λ₁=0.1,控制限系数L=3。将{c₂.1,c₂.2,…,cz.25o}和μe=0.00代入公式(4),即可得到EWMA控制图的统计量²。将μ.=0.00和σ²=0.1497代入公式(6)、公式(7)可计算EWMA控制图的控制限。将统计量x;和控制限绘Z:时间/s图B.2近似层小波系数EWMA控制图从控制图中可看出，EWMA统计量在第580～904秒超出了控制限，可判断在此段时间内流式数据的均值发生了变化。对于两层细节层小波系数，分别使用EWMS控制图监测变量的方差，选择平滑系数λ₂=0.1,控制限系数α=0.0027。将{d₁.1,d₁.2,…,d₁.so}、pa=0.00、σ=0.0038和{d₂.1,d₂.z,…,d2.2so}、μa=0.00、σ²,=0.0075分别代入公式(8),即可得到EWMS控制图的统计量s₁.和s₂.。将λ₂=0.1代入公式(9)、公式(10)可计算EWMS控制图的控制限。将统计量s₁.、s₂.和控制限绘制成图，如图B.3和时间/s图B.3第1层细节层小波系数EWMS控制图时间/s图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秒内出现了超出控制限的情况。由于两个控制图分别包含原始数据在不同频域上的信息，应将两个控制图的结果结合在一起考虑，最终可判断流式数据在上述时间段内变量的方差发生了变化。(资料性)多变量流式数据应用实例在某直拉单晶硅生产过程中，涉及两个质量相关的变量A和B(该数据已进行数据脱敏),过程中数据的获取间隔是1s。生产系统已经平稳运行了一段时间，获取了足够长的平稳数据X。在建立控制图前首先要对控制图参数的目标值进行估算。当所获得历史数据存在异常值时，基于工程经验判断和统计分析，对异常值进行剔除，仅使用平稳数据进行参数的目标值估计。将过程平稳状态下获取的平稳数据X的每一个变量都用离散小波变换分解到2层，并重组得到一层近似层小波系数向量{cz.1,c₂.2,…}和两层细节层小波系数向量{d₁.i,d₁.₂,…}和{d₂.t,d₂.z,…}。根据这些小波系数向量可估计出各层小波系数的均值和方差的目标值，近似层小波系数向量均值为μ。=[0.00,0.00],协方差矩阵为=[0.1497,-0.0102;-0.0102,0.0579]。两层细节层小波系数向量均值为μa=μa₂=[0.00,0.00],协方差矩阵为a₁=[0.0038,0.0000;0.0000,0.0405],2a₂=[0.0075,-0.0017;—0.0017,0.0520]。现获取了1000组新的数据x={x₁,x₂,…,x₁oo},其中x₁=[x{A],x{B]]。两个变量如图C.1所示。时间/s图C.1变量示意图对于获取的数据X={x₁,x₂,…,xtooo}进行控制图监测。应用公式(2)、公式(3),并代入j=2,即可使用小波变换分别将两个变量的数据分解到2层，每个变量都可得到一层近似层小波系数{cz.1,c2.2,…,C2.250}和两层细节层小波系数{d₁.1,d₁.2,…,d1.so}和{d₂.1,d2.2,…,d₂.zso}。将两个变量对应层的小波系数重新组合，则每一层都包含两个变量的小波系数，得到小波系数向量{cz.1,c₂.2,…,c₂.250}、{d₁.i,d₁.2,…,d₁.soo}和{dz.1,dz.z,…,d₂.25o}。对于近似层小波系数向量，使用MEWMA控制图监测变量的均值，选择平滑系数λ₃=0.1。将{C₂.1,C₂.2,…,C₂.25o}和μ.代入公式(11),即可得到迭代向量Z;。再将Z;和2.代入公式(12)、公式(13),即可得到MEWMA控制图的统计量Y;。根据表4可知，该控制图的控制限hA=10.08。将统计量Y?和控制限hA绘制成图，如图C.2所示。Y²时间/s图C.2近似层小波系数MEWMA控制图从控制图中可看出，MEWMA统计量在第428～656秒、第684～806秒和第888～904秒内超出了控制限，可判断在此段时间内数据中的某一个或某几个变量发生了均值的变化。对于两层细节层小波系数向量，使用MEWMC控制图监测变量的方差和协方差，选择平滑系数λ₄=0.1。将{d₁.1,d₁.2,…,d₁.soo}、μa₁、Ea₁和{d₂.1,d₂.2,…,d₂.250}、ua₂、Za₂分别代入公式(14),即可得到转换后的向量序列U₁.;和U₂x。再将λ₄=0.1、U₁.;和U₂;代入公式(15)、公式(16),即可计算MEWMC控制图的统计量y和y₂。根据表5可知，该控制图的控制限hc=0.63。将统计量y.、y₂.;和控制限hc绘制成图，如图C.3和图C.4所示。时间/s图C.3第一层细节层小波系数MEWMC控制图GB/T40859—2021时问/s图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秒内出现了超出控制限的情况。由于两个控制图分别包含原始数据在不同频域上的信息，因此将两个控制图的结果结合在一起考虑，最终可判断原始数据在上述时间段内某一个或某几个变量的方差或几个变量间的协方差发生了变化。(资料性)代码实现D.1概述为了实现对生产过程的监测，需要使用计算机的辅助。且在本文件介绍的算法中，小波变换、矩阵计算等复杂的数学计算往往也可通过编程软件的内置函数进行。Python是目前应用最广泛的数据处理语言之一，以下是使用Python实现绘制两种控制图的主要代码。D.2单变量流式数据控制图构建算法classEWMA():H这是用于建立EWMA控制图的类Hdefinit(self,lambda_1=0.1,mu_c=0,slambda_1:平滑系数mu_c:均值目标值sigma_c:方差目标值level:小波系数的分解层数L:控制限系数self.value=[]#记录所有的统计量的值self.old_number=mu_c#用于表示z(i-1),其中z(0)=mu_cself.new_number=0#用于表示z(i)制限#上控#下控defupdate(self,data):获取新的小波系数后计算对应的控制图统计量data:新获取的小波系数n=len(data)#新获取的小波系数的总长度foriinrange(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_numberdefHplot(self,length):H将最新获取的部分数据绘制成控制图length:绘制控制图的长度(原始数据的长度)##Htotal_length=len(self.value)#目前已经计算出的统计量的长度adjusted_length小波系数的长度plot_length=int(math.floor(length*2*"(-self.level)))#原始数据的长度对应到该层=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-',linewidthplt.text(0,self.ucl,'s\mathregular{U_{CL}}S')plt.plot(t_label,self.lcl"np.ones(plot_length),'k--',linewidthplt.text(0,self.lcl,'s\mathregular{L_{CL}}s')plt.ylabel('z$_is',rotation=360,fontsize=16)plt.show()=0.6)#绘制上控制限=0.6)#绘制下控制限importmathimportnumpyasnpimportmatplotlib.pyplotaspltclassEWMS():这是用于建立EWMS控制图的类definit_(self,lambda_2=0.1,mu_d=0,sigma_d=1,level=1,ucl=1.32,lcl=0.67);lambda_2:平滑系数mu_c;均值目标值sigma_c:方差目标值level:小波系数的分解层数ucl:上控制限参数(标准差的倍数)lcl:下控制限参数(标准差的倍数)self.lambda_2=lambda_2self.mu_d=mu_dself.sigma_d=sigma_dself.level=levelself.value=[]#记录所有的统计量的值self.old_number=sigma_d#用于表示s^2(i-1),其中s^2(0)=sigma_dself.new_number=0#用于表示s^2(i)self.ucl=ucl*self.sigma_d**0.5#上控制限self.lcl=lcl*self.sigma_d**0.5井下控制限defupdate(self,data):获取新的小波系数后计算对应的控制图统计量data:新获取的小波系数n=len(data)#新获取的小波系数的总长度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_numberdefplot(self,length):将最新获取的部分数据绘制成控制图length:绘制控制图的长度(原始数据的长度)total_length=len(self.value)#目前已经计算出的统计量的长度小波系数的长度=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-',l=0.6)#绘制统计量plt.plot(t_label,self.ucl*np.ones(plot_length),'k-',linewidth=0.6)#绘制上控制限plt.text(0,self.ucl,'s\mathregular{U_{CL}}s')plt.plot(t_label,self.lcl*np.ones(plot_length),'k--',linewidthplt.ylabel('s$_is',rotation=360,fontsize=16)plt.show()构建控制图：假设构建控制图前获取了一组受控数据'data'。importnumpyasnpfromEWMAimportEWMA#引用前面定义的EWMA类fromEWMSimportEWMS#引用前面定义的EWMS类importpywt#用于小波变换的程序包#==========================================#这里应给定用于计算控制图参数目标值的平稳历史数据，赋值给data变量#也可以应用下行代码生成1000个仿真样本数据进行测试#通过受控数据计算控制图参数的目标值#对应7.4步骤1,将受控数据分解为各层小波系数，设置控制图参数并计算各参数的目标值J=2#分解层数为2wavelet_coefficient=pywt.wavedec(data,'haar',level=J)#使用离散小波变换将受控数据分解为各层小波系数#设置EWMA控制图的平滑系数，并计算近似层小波系数的均值、方差的目标值lamda_A2=0.1mu_A2=np.mean(wavelet_coefficient[0],axis=0)#设置第二层EWMS控制图的平滑系数，并计算第二层细节层小波系数的均值、方差的目标值lamda_D2=0.1mu_D2=np.mean(wavelet_coefficient[1],axis=0)sigma_D2=np.cov(wavelet_coefficient[1],rowvar=0)#设置第一层EWMS控制图的平滑系数，并计算第一层细节层小波系数的均值、方差的目标值lamda_D1=0.1GB/T40859—2021mu_D1=np.mean(wavelet_coefficient[2],axis=0)sigma_Dl=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([datal,data2,…])#也可以用下行代码生成1000个仿真样本数据进行测试test_data=np.random.normal(0,1,1000)N=len(test_data)#用于监测的样本数据的长度old_data=np.array([])#保存用于计算小波系数的数据i=0#计数为0,该计数用于判断是否更新小波系数whilen<N:#每当获得一组新的数据时，进行以下操作：new_data=test_data[forjinrange(1,J+1):#从最低层开始看，是否可以更新系数ifi%2*j==0:#如果i整除2^j,说明第j层可以更新一个系数=old_data[len(old_data)-2**j+1:len(old_data)]=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类ifj==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))#将新数据纳入老数据中，重新获取数据ifi==2*J:#如果计数达到了2^J,则归0na=np.array([])#清空用于计算小波系数的数据foriinrange(0,J+1):plt.rcParams['figure.dpi']=300control_chart[i].plot(N)#分别绘制各个控制图D.3多变量流式数据控制图构建算法importmathimportnumpyasnpimportmatplotlib.pyplotasplt这是用于建立MEWMA控制图的类definit(self,lambda_3=0.1,mu_c=0,sigma_c=1,level=1,ucl=8.64);HHHlambda_3:平滑系数mu_c:均值向量目标值sigma_c:协方差矩阵目标值level:小波系数的分解层数ucl:控制限HHself.lambda_3=lambda_3self.sigma_c=np.mat(self.value=[]#记录所有的统计量的值=mu_c#用于表示Z(i-1),=np.array([])#Z(i)其中Z(O)=mu_cself.ucl=ucl#控制限defupdate(self,data):#P1获取新的小波系数后计算对应的控制图统计量data:新获取的小波系数向量HH[n,p]=data.shape#n为小波系数的总长度，p为变量维数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_cself.value.append(np.float(np.dot(np.dot(trans,np.linalg.inv(sigma_z)),trans.T)))#对应公式(12)计算统计量并保存defself.old_number=self.new_numberdefplot(self,length):H将最新获取的部分数据绘制成控制图length:绘制控制图的长度(原始数据的长度)Htotal_length小波系数的长度plot_length=len(self.value)#目前已经计算出的统计量的长度=int(math.floor(length*2**(-self.level)))#原始数据的长度对应到该层=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-',linewidthplt.text(0,self.ucl,'$\mathregular{U_{CL}}$)plt.ylabel('$\mathregular{Y_i^2}$',rotation=360,fontsizeplt.show()=0.6)#绘制上控制限=16)importmathimportmatplotlib.pyplotaspltH这是用于建立MEWMC控制图的类definit(self,lambda_4=0.1,mu_d=0,sigma_d=1,level=1,ucl=1.169):H初始化MEWMC类lambda_4:平滑系数mu_d:均值向量目标值sigma_d:协方差矩阵目标值level:小波系数的分解层数ucl:控制限self.sigma_d=sigma_dself.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))foriinrange(0,len(eigVects[:,i]=eigVects[:,i]*(1/np.sqrt(eigVals[i]))self.Amatrix=eigVects.T#计算公式(14)中的A矩阵defupdate(self,data):获取新的小波系数后计算对应的控制图统计量data:新获取的小波系数向量[n,p]=data.shape#n为小波系数的总长度，p为变量维数trans=data[i,:]-self.mu_dU=np.dot(self.Amatrix,trans.T)UUT=np.dot(U,U.T)#对应公式(15)中的U(i)*U(i)^Tself.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_numberdefplot(self,length):def将最新获取的部分数据绘制成控制图length:绘制控制图的长度(原始数据的长度)total_length=len(self.value)#目前已经计算出的统计量的长度小波系数的长度=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-',linewidthplt.text(0,self.ucl,'s\mathregular{U_{CL}}$')plt.ylabel('$\mathregular{y_i}$',rotation=360,fontsizeplt.show()=0.6)#绘制上控制限=16)构建控制图：假设构建控制图前获取了一组p=2维的受控数据，分别记为a和b。importmatplotlib.pyplotaspltfromMEWMAimportMEWMA#引用前面定义的MEWMA类fromMEWMCimportMEWMC#引用前面定义的MEWMC类importpywt#用于小波变换的程序包#========================================#这里应给定用于计算控制图参数目标值的平稳历史数据，赋值给data变量#a=np.array([al,a2,…])#b=np.array([bl,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#分解层数为2wavelet_coefficient_a=pywt.wavedec(a,'haar',level=J)wavelet_coefficient_b=pywt.wavedec(b,'haar',level=J)#对应8.2,分别对各个变量进行小波变换，并将小波系数重组为系数向量wavelet_coefficient=[]range(0,len(wavelet_coefficient_a)):#将对应层系数组成矩阵=np.vstack((wavelet_coefficient_a[i],wavelet_coefficient_b[i]))=temp.Twavelet_coefficient.append(temp)#设置MEWMA控制图的平滑系数，并计算近似层小波系数向量的均值向量、协方差矩阵的目标值mu_A2=np.mean(wavelet_coefficient[0],axis=0)sigma_A2=np.cov(wavelet_coefficient[0],rowvar=0)#设置第二层MEWMC控制图的平滑系数，并计算第二层细节层小波系数向量的均值向量、协方差矩阵的目标值lamda_D2=0.1mu_D2=np.mean(wavelet_coefficient[1],axis=0)sigma_D2=np.cov(wavelet_coefficient[1],rowvar=0)#设置第一层MEWMC控制图的平滑系数，并计算第二层细节层小波系数向量的均值向量、协方差矩阵的目标值lamda_D1=0.1mu_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[o])#查表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])