发明名称 基于多核独立元分析的批量生产过程故障检测方法
摘要 一种基于多核独立元分析的批量生产过程故障检测方法,包括采集数据、数据处理、利用核主元分析对数据进行白化、利用修正ICA提取独立元及利用T2和SPE统计量进行故障检测步骤。本发明首次提出一种基于MKICA的批量生产过程监测方法,实现对复杂过程进行监测。能较早的检测到故障,从而减少工业生产过程中的损失。
申请公布号 CN101158693B 申请公布日期 2011.08.17
申请号 CN200710012956.3 申请日期 2007.09.26
申请人 东北大学 发明人 张颖伟;秦泗钊;王婷
分类号 G05B19/418(2006.01)I;G06F17/00(2006.01)I 主分类号 G05B19/418(2006.01)I
代理机构 沈阳东大专利代理有限公司 21109 代理人 梁焱
主权项 1.一种基于多核独立元分析的批量生产过程故障检测方法,其特征在于将该方法应用于青霉素发酵过程,包括以下步骤:步骤一、采集数据采集过程中相关变量的数据,对于每个故障,采集两组数据,即标准操作模式的训练数据和在线监测的实时工况数据;训练数据由400组数据构成,实时工况数据由800组数据构成;训练数据和实时工况数据都包含了9个监测变量,即:通风率、煽动功率、培养基供给率、培养基温度、溶解氧的浓度、培养容积、CO<sub>2</sub>浓度、PH值和产生的热量;其中,训练数据用于建立模型,实时工况数据用于故障检测;步骤二、数据处理用当前值补充丢失的数据,用均值和标准偏差规范化采集的观测数据,批量生产过程数据是关于批次、变量和时间的数据,首先把数据放入三维矩阵X(I×J×K)中,其中I是批次,J是变量的数目,K是每个批采样的次数,然后将矩阵X(I×J×K)变换成二维矩阵X(I×JK);步骤三、利用核主元分析对数据进行白化处理通过非线性映射将输入空间映射到一个特征空间,接着在此特征空间对观测数据进行白化处理,得到白化后的观测变量z;具体过程如下:首先进行非线性映射,设观测数据x<sub>k</sub>∈R<sup>JK</sup>,k=1,...,I,k是观察数据的个数,I是批次,J是变量的数目,K是每个批采样的次数,R<sup>JK</sup>为JK维空间,通过非线性映射Φ:R<sup>JK</sup>→F,F为特征空间,把原始空间的观测数据映射到高维特征空间中,Φ(x<sub>k</sub>)∈F,特征空间中协方差矩阵为<maths num="0001"><![CDATA[<math><mrow><msup><mi>C</mi><mi>F</mi></msup><mo>=</mo><mfrac><mn>1</mn><mi>I</mi></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover><mi>&Phi;</mi><mrow><mo>(</mo><msub><mi>x</mi><mi>k</mi></msub><mo>)</mo></mrow><mi>&Phi;</mi><msup><mrow><mo>(</mo><msub><mi>x</mi><mi>k</mi></msub><mo>)</mo></mrow><mi>T</mi></msup></mrow></math>]]></maths>其中Φ(x<sub>k</sub>),k=1,...,I,为映射到高维特征空间的数据,假定Φ(x<sub>k</sub>)为零均值和单位方差,令Θ=[Φ(x<sub>1</sub>),...,Φ(x<sub>I</sub>)],因而协方差矩阵C<sup>F</sup>能表示为<img file="FSB00000353017400012.GIF" wi="275" he="117" />定义一个I×IGram维的核矩阵K:[K]<sub>ij</sub>=K<sub>ij</sub>=〈Φ(x<sub>i</sub>),Φ(x<sub>j</sub>)〉=k(x<sub>i</sub>,x<sub>j</sub>)(1)有K=Θ<sup>T</sup>Θ,k(x<sub>i</sub>,x<sub>j</sub>)为核函数,x<sub>i</sub>和x<sub>j</sub>为观测数据,1≤i,j≤I,应用核函数k(x<sub>i</sub>,x<sub>j</sub>),可以不通过非线性映射就能在F中计算内积,这样就避免了在特征空间中执行非线性映射和内积计算,常用的核函数是径向基核函数<img file="FSB00000353017400021.GIF" wi="552" he="180" />多项式核函数k(x,y)=〈x,y〉<sup>r</sup>和S形核函数k(x,y)=tanh(β<sub>0</sub>〈x,y〉+β<sub>1</sub>),核函数的选择决定了映射Φ和特征空间F;由核矩阵K可知,高维空间中Φ(x<sub>k</sub>)能按以下方法进行中心化,即将Φ(x<sub>k</sub>)的中心化转化为对K的中心化处理,从下式可以得到中心化核矩阵<img file="FSB00000353017400022.GIF" wi="66" he="49" /><maths num="0002"><![CDATA[<math><mrow><mover><mi>K</mi><mo>~</mo></mover><mo>=</mo><mi>K</mi><mo>-</mo><msub><mn>1</mn><mi>I</mi></msub><mi>K</mi><mo>-</mo><msub><mrow><mi>K</mi><mn>1</mn></mrow><mi>I</mi></msub><mo>+</mo><msub><mn>1</mn><mi>I</mi></msub><msub><mrow><mi>K</mi><mn>1</mn></mrow><mi>I</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中,<img file="FSB00000353017400024.GIF" wi="552" he="246" />对<img file="FSB00000353017400025.GIF" wi="41" he="47" />进行特征值分解:<maths num="0003"><![CDATA[<math><mrow><mi>&lambda;&alpha;</mi><mo>=</mo><mover><mi>K</mi><mo>~</mo></mover><mi>&alpha;</mi><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow></math>]]></maths>可以得到<img file="FSB00000353017400027.GIF" wi="41" he="48" />的d个最大正特征值λ<sub>1</sub>≥λ<sub>2</sub>≥...≥λ<sub>d</sub>和相应的标准正交的特征向量α<sub>1</sub>,α<sub>2</sub>,...,α<sub>d</sub>,按经验选取满足<img file="FSB00000353017400028.GIF" wi="326" he="127" />的特征值,i=1,...,d,则协方差矩阵C<sup>F</sup>的d个最大正特征值是<img file="FSB00000353017400029.GIF" wi="294" he="108" />相应的标准正交特征向量v<sub>1</sub>,v<sub>2</sub>,...,v<sub>d</sub>可以表示为<maths num="0004"><![CDATA[<math><mrow><msub><mi>v</mi><mi>j</mi></msub><mo>=</mo><mfrac><mn>1</mn><msqrt><msub><mi>&lambda;</mi><mi>j</mi></msub></msqrt></mfrac><mi>&Theta;</mi><msub><mi>&alpha;</mi><mi>j</mi></msub><mo>,</mo><mi>j</mi><mo>=</mo><mn>1</mn><mo>,</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>,</mo><mi>d</mi><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow></math>]]></maths>C<sup>F</sup>的特征矩阵V=[v<sub>1</sub>,v<sub>2</sub>,...,v<sub>d</sub>]可以简单的表示如下:V=ΘHΛ<sup>-1/2</sup>        (5)其中Λ=diag(λ<sub>1</sub>,λ<sub>2</sub>,...,λ<sub>d</sub>)是由d个特征值构成的对角矩阵,H=[α<sub>1</sub>,α<sub>2</sub>,...,α<sub>d</sub>]是对应的特征向量构成的矩阵,用协方差矩阵C<sup>F</sup>的特征矩阵V把C<sup>F</sup>化为对角矩阵:<maths num="0005"><![CDATA[<math><mrow><msup><mi>C</mi><mi>F</mi></msup><mo>=</mo><mi>Vdiag</mi><mrow><mo>(</mo><mfrac><msub><mi>&lambda;</mi><mn>1</mn></msub><mi>I</mi></mfrac><mo>,</mo><mfrac><msub><mi>&lambda;</mi><mn>2</mn></msub><mi>I</mi></mfrac><mo>,</mo><mo>&CenterDot;</mo><mo>&CenterDot;</mo><mo>&CenterDot;</mo><mo>,</mo><mfrac><msub><mi>&lambda;</mi><mi>d</mi></msub><mi>I</mi></mfrac><mo>)</mo></mrow><msup><mi>V</mi><mi>T</mi></msup><mo>=</mo><mfrac><mn>1</mn><mi>I</mi></mfrac><mi>V&Lambda;</mi><msup><mi>V</mi><mi>T</mi></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>6</mn><mo>)</mo></mrow></mrow></math>]]></maths>令<maths num="0006"><![CDATA[<math><mrow><mi>P</mi><mo>=</mo><mi>V</mi><msup><mrow><mo>(</mo><mfrac><mn>1</mn><mi>I</mi></mfrac><mi>&Lambda;</mi><mo>)</mo></mrow><mrow><mo>-</mo><mn>1</mn><mo>/</mo><mn>2</mn></mrow></msup><mo>=</mo><msqrt><mi>I</mi></msqrt><mi>&Theta;H</mi><msup><mi>&Lambda;</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>7</mn><mo>)</mo></mrow></mrow></math>]]></maths>那么,P<sup>T</sup>C<sup>F</sup>P=E    (8)P为白化矩阵,E为单位矩阵,这样就得到了白化矩阵P;所述的求解白化后的观测变量z的具体过程如下:映射到特征空间中的数据Φ(x)经过如下转换可以被白化:z=P<sup>T</sup>Φ(x)(9)z为白化后的观测变量,<maths num="0007"><![CDATA[<math><mrow><mi>z</mi><mo>=</mo><msup><mi>P</mi><mi>T</mi></msup><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><msqrt><mi>I</mi></msqrt><msup><mi>&Lambda;</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><msup><mi>H</mi><mi>T</mi></msup><msup><mi>&Theta;</mi><mi>T</mi></msup><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><msqrt><mi>I</mi></msqrt><msup><mi>&Lambda;</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><msup><mi>H</mi><mi>T</mi></msup><mo>[</mo><mi>&Phi;</mi><mrow><mo>(</mo><msub><mi>x</mi><mn>1</mn></msub><mo>)</mo></mrow><mo>,</mo><mo>&CenterDot;</mo><mo>&CenterDot;</mo><mo>&CenterDot;</mo><mo>,</mo><mi>&Phi;</mi><mrow><mo>(</mo><msub><mi>x</mi><mi>I</mi></msub><mo>)</mo></mrow><msup><mo>]</mo><mi>T</mi></msup><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow></mrow></math>]]></maths><maths num="0008"><![CDATA[<math><mrow><mo>=</mo><msqrt><mi>I</mi></msqrt><msup><mi>&Lambda;</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><msup><mi>H</mi><mi>T</mi></msup><msup><mrow><mo>[</mo><mover><mi>k</mi><mo>~</mo></mover><mrow><mo>(</mo><msub><mi>x</mi><mn>1</mn></msub><mo>,</mo><mi>x</mi><mo>)</mo></mrow><mo>,</mo><mo>&CenterDot;</mo><mo>&CenterDot;</mo><mo>&CenterDot;</mo><mo>,</mo><mover><mi>k</mi><mo>~</mo></mover><mrow><mo>(</mo><msub><mi>x</mi><mn>1</mn></msub><mo>,</mo><mi>x</mi><mo>)</mo></mrow><mo>]</mo></mrow><mi>T</mi></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>10</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中,x是需要白化的数据,<img file="FSB00000353017400033.GIF" wi="186" he="68" />j=1,2,...,I,为中心化核函数;步骤四、利用修正ICA提取独立元利用修正ICA算法把白化后的观测变量z转换成独立元,并使独立元各变量之间尽可能相互统计独立;求解独立元y的具体过程如下:修正ICA方法可以找到p个独立元,y={y<sub>1</sub>,...,y<sub>p</sub>},为使元素彼此之间能尽可能的独立,y需满足E{yy<sup>T</sup>}=D=diag{λ<sub>1</sub>,...,λ<sub>p</sub>},λ<sub>1</sub>,...,λ<sub>p</sub>是p个独立元所对应的特征值,利用y=C<sup>T</sup>z    (11)其中C ∈R<sup>d×p</sup>为得分转换矩阵,C<sup>T</sup>C=D;定义归一化的独立元y<sub>n</sub>为:y<sub>n</sub>=D<sup>-1/2</sup>y=D<sup>-1/2</sup>C<sup>T</sup>z=C<sub>n</sub><sup>T</sup>z    (12)C<sub>n</sub>为标准得分转换矩阵,显然,D<sup>-1/2</sup>C<sup>T</sup>=C<sub>n</sub><sup>T</sup>,C<sub>n</sub><sup>T</sup>C<sub>n</sub>=E,E{y<sub>n</sub>y<sub>n</sub><sup>T</sup>}=E,由于z是观测数据经过白化后得到的,具有不相关性,所以可以选择z的第一个数据中p个独立元作为y<sub>n</sub>的初始值,因此可以把C<sub>n</sub><sup>T</sup>的初始矩阵设置为:<img file="FSB00000353017400034.GIF" wi="1556" he="105" />其中I<sub>p</sub>是p维单位矩阵,0是p ×(d-p)维零矩阵;利用修正ICA算法计算标准得分转换矩阵C<sub>n</sub>,由式(11)和(12)就可以得到独立元y,独立元可由下式获得:y=D<sup>1/2</sup>C<sub>n</sub><sup>T</sup>z=D<sup>1/2</sup>C<sub>n</sub><sup>T</sup>P<sup>T</sup>Φ(x)(14);步骤五、利用T<sup>2</sup>和SPE统计量进行故障检测采用T<sup>2</sup>和SPE统计量进行在线故障检测,当观测数据的统计量没有超出统计量规定的控制限时,则属于正常数据,反之属于异常数据,表明出现故障;批量生产过程统计变量的T<sup>2</sup>和SPE的计算:T<sup>2</sup>统计定义如下:T<sup>2</sup>=y<sup>T</sup>D<sup>-1</sup>y,<maths num="0009"><![CDATA[<math><mrow><mi>y</mi><mo>=</mo><msup><mi>D</mi><mrow><mn>1</mn><mo>/</mo><mn>2</mn></mrow></msup><msubsup><mi>C</mi><mi>n</mi><mi>T</mi></msubsup><msup><mi>P</mi><mi>T</mi></msup><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>15</mn><mo>)</mo></mrow></mrow></math>]]></maths>SPE统计定义如下:<maths num="0010"><![CDATA[<math><mrow><mi>SPE</mi><mo>=</mo><msup><mi>e</mi><mi>T</mi></msup><mi>e</mi><mo>=</mo><msup><mrow><mo>(</mo><mi>z</mi><mo>-</mo><mover><mi>z</mi><mo>^</mo></mover><mo>)</mo></mrow><mi>T</mi></msup><mrow><mo>(</mo><mi>z</mi><mo>-</mo><mover><mi>z</mi><mo>^</mo></mover><mo>)</mo></mrow><mo>=</mo><msup><mi>z</mi><mi>T</mi></msup><mrow><mo>(</mo><mi>E</mi><mo>-</mo><msub><mi>C</mi><mi>n</mi></msub><msup><msub><mi>C</mi><mi>n</mi></msub><mi>T</mi></msup><mo>)</mo></mrow><mi>z</mi><mo>,</mo><mi>z</mi><mo>=</mo><msup><mi>P</mi><mi>T</mi></msup><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>16</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中,<img file="FSB00000353017400043.GIF" wi="186" he="43" /><img file="FSB00000353017400044.GIF" wi="23" he="43" />可由下式得到:<maths num="0011"><![CDATA[<math><mrow><mover><mi>z</mi><mo>^</mo></mover><mo>=</mo><msub><mi>C</mi><mi>n</mi></msub><msup><mi>D</mi><mrow><mo>-</mo><mn>1</mn><mo>/</mo><mn>2</mn></mrow></msup><mi>y</mi><mo>=</mo><msubsup><mrow><msub><mi>C</mi><mi>n</mi></msub><mi>C</mi></mrow><mi>n</mi><mi>T</mi></msubsup><mi>z</mi><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>17</mn><mo>)</mo></mrow></mrow></math>]]></maths>因为y不服从高斯分布,所以T<sup>2</sup>的控制上界可由F分布确定;如果大部分非高斯性包含在提取的独立元中,剩余的子空间将包含大部分的高斯噪声,这个噪声可以认为是正态分布的,假设预测误差是正态分布的,SPE的控制界限可由下式加权χ<sup>2</sup>分布计算:<maths num="0012"><![CDATA[<math><mrow><mi>SPE</mi><mo>~</mo><msubsup><mi>&mu;&chi;</mi><mi>h</mi><mn>2</mn></msubsup></mrow></math>]]></maths>μ=b/2a,h=2a<sup>2</sup>/b(18)其中a和b分别是标准操作数据中SPE的估计均值和方差。
地址 110004 辽宁省沈阳市和平区文化路3号巷11号