基于增量主成份分析的牛奶凝结过程实时监测方法

摘要

本发明公开了一种基于增量主成份分析的牛奶凝结过程实时监测方法，它是在牛奶凝结过程中的近红外光谱数据基础上，利用聚类方法将数据聚类到多个窗口中，然后使用增量主成份分析方法对窗口内数据进行降维，而且当后续的批次数据加入窗口内时，可以通过新增加的样本数据来修正原来的数据，从而可以对窗口内的特征值和特征向量进行持续的更新，最终动态窗口中的特征数据将精确反映出牛奶凝结过程。本发明的方法可以对近红外光谱数据进行更为可靠的数据分析，最终实现了对牛奶凝结过程的精准监测。通过反复试验验证，本发明的方法在对牛奶凝结过程进行检测时，准确率高达90％以上，从而可以用于更高质量奶酪产品的生产。

说明书

技术领域

本发明涉及一种牛奶凝结过程实时监测方法，特别涉及一种基于增量主成份分析的牛奶 凝结过程实时监测方法。

背景技术

牛奶凝结用于制作奶酪的步骤是通过切割形成乳凝胶，所切割的时间最终影响奶酪的产 量和品质。如果切割的时间太早导致产量的损失，切割的时间太晚导致得到高水分的奶酪。 现在的奶酪乳制品制作中，工厂操作员主要是根据他们的工作经验手工切割牛奶凝结形成的 乳凝胶，不仅成本高还带有主观因素，这阻碍了流程的自动化和影响了工作效率。因此，我 们有必要优化一种确定切割时间的自动化方法。

牛奶凝乳酶诱导凝胶是需要经历三个基本阶段：

(1)初始酶水解，J-酪蛋白将改变其酪蛋白胶束，最终形成副酪蛋白。

(2)副酪蛋白的聚合，其中聚合速率取决于自由副酪蛋白的浓度，意味着这一阶段是依 赖于速率和第J酪蛋白水解的程度；

(3)凝胶过程，形成聚合物网络聚合胶束链，也称为凝胶固化。过度之间的阶段，不易 察觉，因为头部和尾部连续的阶段在一定程度上是重叠的。

发明内容

本发明的目的在于，提供一种基于增量主成份分析的牛奶凝结过程实时监测方法。它可 以精确检测牛奶凝结过程中主要成分的变化，从而实现牛奶凝结过程的实时监测，以用于作 为切割乳凝胶的时间标准，进而提高奶酪的产量和品质。

本发明的技术方案：基于增量主成份分析的牛奶凝结过程实时监测方法，包括以下步 骤：

S1.采集多个批次的牛奶凝结过程中各个时间点的近红外光谱数据作为原始数据；

S2.对原始数据进行标准化的预处理，得到多个批次的标准化数据；

S3.通过聚类算法将第一个批次的标准化数据进行聚类操作；，生成多个聚类，并将这些 聚类标识为不同的窗口；

S4.对每个窗口中的标准化原始数据进行主成份分析(常规主成份分析方法)，得到每 个窗口内的数据载荷和数据对应的主成份得分；

S5.根据时间点的对应关系，将对下一个批次的标准化数据增加到各个窗口中，并使用 增量主成份分析方法对窗口中的新增数据进行修正以及预测它们的主成份得分，并形成当前 窗口的数据特征；

S6.对剩余批次标准化数据逐个重复步骤S5；

S7.采集待监测的牛奶凝结过程中的近红外光谱实时数据，并以步骤S5相同的方法将近 红外光谱实时数据新增到窗口中，通过观察每个窗口内近红外光谱数据的主成份得分来判断 牛奶凝结处于哪个阶段(通过窗口去实现观察)。

上述的基于增量主成份分析的牛奶凝结过程实时监测方法中，所述步骤S2的预处理使 用标准化函数，即$X_0=\frac{X-\bar{X}}{\sqrt{\frac{{(X-\bar{X})}^2}{n}}}.$

前述的基于增量主成份分析的牛奶凝结过程实时监测方法中，所述步骤S3中采用欧式 距离作为聚类的指标。多维欧式距离d计算公式为

$d=\sqrt{\underset{i=1}{\overset{m}{\Sigma }}{(x_{i1}-x_{i2})}^2}$

其中x_i1代表第一个点的第i维坐标，x_i2代表第二个点的的第i维坐标。的第1维坐 标，x_i2代表x_i2的第2维坐标，d代表多维空间的点点之间的几何距离。前述的基于增量主 成份分析的牛奶凝结过程实时监测方法中，所述步骤7中每个窗口内的数据特征包括窗口内 数据的个数、窗口内数据的主成份载荷、主成份的得分及得分所对应的时间点。

前述的基于增量主成份分析的牛奶凝结过程实时监测方法中，所述步骤S5中，当下一 个批次的标准化原始数据中有不符合加入已有窗口的条件的数据，则建立新的窗口，并将这 些数据放入新建的窗口，并使用主成份分析方法进行载荷和得分的计算；如果窗口数量达到 上限，则合并最近的两个窗口，否则直接建立一个新的窗口。

所述步骤S5中的增量主成份分析方法(IPCA)，它不必假设输入样本为零均值并可以 动态地更新均值，其描述如下：

给定一个d*n的近红外光谱数据矩阵X_n＝{X₁,X₂,...X_n}，其中，X_i表示第i个批次 的数据(d表示d维数向量)。

针对第一批次的近红外光谱的数据,利用PCA算法可以求得PCA得分矩阵F₀：

X₀Q₀＝F₀

其中P₀表示第0批次数据的载荷，即单位特征向量。

$Q_0=X_0^{-1}{F}_0={\left(X_0^T{X}_0\right)}^{-1}{X}_0^T{F}_0=K_0^{-1}{X}_0^T{F}_0$

令K₀是一个常数。

当给定新输入样本X₁，有如下

$\left(\begin{array}{c}{X}_0\\ X_1\end{array}\right)F_1=\left(\begin{array}{c}{X}_0\\ X_1\end{array}\right)$

$Q_1={\left(\begin{array}{c}{X}_0\\ X_1\end{array}\right)}^{-1}\left(\begin{array}{c}{F}_0\\ F_1\end{array}\right)={\left({\left(\begin{array}{c}{X}_0\\ X_1\end{array}\right)}^T\left(\begin{array}{c}{X}_0\\ X_1\end{array}\right)\right)}^{-1}{\left(\begin{array}{c}{X}_0\\ X_1\end{array}\right)}^T\left(\begin{array}{c}{F}_0\\ F_1\end{array}\right)=K_1^{-1}{\left(\begin{array}{c}{X}_0\\ X_1\end{array}\right)}^T\left(\begin{array}{c}{F}_0\\ F_1\end{array}\right)$

其中Q₁表示第1批次数据的载荷。

其中$K_1={\left(\begin{array}{c}{X}_0\\ X_1\end{array}\right)}^T\left(\begin{array}{c}{X}_0\\ X_1\end{array}\right)=\left[X_0^T{X}_1^T\right]\left(\begin{array}{c}{X}_0\\ X_1\end{array}\right)=X_0^T{X}_0+X_1^T{X}_1=K_0+X_1^T{X}_1$

$K_1=K_0+X_1^T{X}_1=X_0^T{X}_0+X_1^T{X}_1$

$\left(\begin{array}{c}{\left(\begin{array}{c}{X}_0\\ X_1\end{array}\right)}^T\left(\begin{array}{c}{F}_0\\ F_1\end{array}\right)=X_0^T{F}_0+X_1^T{F}_1=K_0{K}_0^{-1}{X}_0^T{F}_0+X_1^T{F}_1=K_0{P}_0+X_1^T{F}_1\\ =\left(K_1-X_1^T{T}_1\right)Q_0+X_1^T{F}_1\end{array}\right)$

$Q_1=K_1^{-1}\left(\right(K_1-X_1^T{F}_1)Q_0+X_1^T{F}_1)=Q_0+K_1^{-1}{X}_1^T(F_1-X_1{Q}_0)$

综上所述，对于每一个新输入样本X_n，都可以根据之前的样本对新输入的样本进行修 正，并预测它在整个集合中的特征值和特征向量。

$Q_{k+1}=Q_k+K_{k+1}^{-1}{X}_{k+1}^T(F_{k+1}-X_{k+1}{Q}_k)$

$K_{k+1}=K_k+X_k^T{X}_k$

与现有技术相比，本发明在牛奶凝结过程中的近红外光谱数据基础上，利用聚类方法将 数据聚类到多个窗口中，然后使用增量主成份分析方法对窗口内数据进行降维，而且当后续 的批次数据加入窗口内时，可以通过新增加的样本数据来修正原来的数据，从而可以对窗口 内的特征值和特征向量进行持续的更新，最终动态窗口中的特征数据将精确反映出牛奶凝结 过程。本发明的方法结合了基于窗口的聚类方法和增量主成份分析方法，可以对近红外光谱 数据进行更为可靠的数据分析，最终实现了对牛奶凝结过程的精准监测。通过反复试验验 证，本发明的方法在对牛奶凝结过程进行检测时，准确率高达90％以上，可以用于生产更高 质量奶酪。

附图说明

图1是本发明的流程示意图；

图2是本发明的实施例的近红外光谱数据维数图；

图3是本发明的实施例的第一批次PC1得分百分比图；

图4是第二批次数据的修正并预测PC1得分和未经修正PC1关系图；

图5是本发明实施例的误差率示意图。

具体实施方式

下面结合实施例对本发明作进一步的说明，但并不作为对本发明限制的依据。对本发明 所做的任何形式上的变通和/或改变都将落入本发明保护范围。

本发明的实施例：基于增量主成份分析的牛奶凝结过程实时监测方法，如图1所示，包 括以下步骤：

S1.采集多个批次的牛奶凝结过程中各个时间点的近红外光谱数据作为原始数据；

S2.对原始数据进行标准化的预处理，得到多个批次的标准化数据；

S3.通过聚类算法将第一个批次的标准化数据进行聚类操作，生成多个聚类，并将这些 聚类标识为不同的窗口；

S4.对每个窗口中的标准化原始数据进行主成份分析，得到每个窗口内的数据载荷，并 计算出主成份得分；

S5.根据时间点的对应关系，将对下一个批次的标准化原始数据增加到各个窗口中，并 使用增量主成份分析方法对窗口中的新增数据进行修正以及预测它们的载荷和主成份得分， 并形成当前窗口的数据特征；

S6.对剩余批次的标准化数据逐个重复步骤S5；

S7.采集待监测的牛奶凝结过程中的近红外光谱实时数据，并以步骤S5相同的方法将近 红外光谱实时数据归类到相应窗口中，依据窗口内的近红外光谱实时数据的主成分得分来判 断当前牛奶凝结状态，并进行相应处理。e

其具体实施时可采用如下步骤：

本实施例使用的是牛奶凝结过程的十二批次的近红外光谱数据。

步骤1：数据的提取，每个批次中包含的数据结构是一个x对应一个时间点t。其中t 为牛奶凝结过程的时间，x为相对应t时间点的1557维数的近红外光谱数据。

batch₁：(x₁₁,t₁₁),(x₁₂,t₁₂)...(x_1n,t_1n)

......

batch_m：(x_m1,t_m1),(x_m2,t_m2)...(x_mn,t_mn)

步骤2：对数据进行预处理，对batch₁进行预处理得到X₀，预处理使用标准化函数。即

$X_0=\frac{X-\bar{X}}{\sqrt{\frac{{(X-\bar{X})}^2}{n}}}$

步骤3：计算X₀的多维欧式距离d

$d=\sqrt{\underset{i-1}{\overset{m}{\Sigma }}{(x_{i1}-x_{i2})}^2}$

其中x_i1代表第一个点的第i维坐标，x_i2代表第二个点的的第i维坐标。

转成关于多维欧式距离的方阵A₁

以下步骤4-7，通过聚类算法建立窗口，并计算数据的主成分得分。

步骤4：设定初始值i＝1,j＝1,n＝1(1≤n≤k).其中k为A₁的总列数。若a_nj≤σ(σ为 设定的一个阈值)，σ＞0，则原(x_1n,t_1n)放入一个窗口window_i中，若a_n"j＜σ则进入步骤 5，若a_n"j＞σ，则j＝n",n＝n",i＝i+1重复本步骤。步骤4中i代表窗口计数器，表示第i 个窗口，而n和j分别表示矩阵A₁中行和列的坐标，决定a的位置，取行位置变为n"，列位 置也变为n"的a。

步骤5：若n<k,则n＝n+1,并返回步骤4，若n＝k则进入步骤6。

步骤6：窗口内的数据进行主成分分析，计算得到载荷矩阵和主成分得分。

对于window_i内的数据X，进行如下处理

X＝PΔQ^T

F＝PΔ＝PΔQ^T＝XQ

得到window_i的载荷Q_i和初始时间t_i1和截止时间t_i2、(x_1n,t_1n)的PCA得分F。

步骤7：对第二批次的数据，利用增量PCA算法对载荷进行预测。

对于batch₂:(x₂₁,t₂₁),(x₂₂,t₂₂)...(x_2n,t_2n)，若t_i1≤t_2n≤t_i2，则将x_2n加入到窗口 window_i，并得到的x_2n修正后的PCA得分为利用增量PCA算法对载荷进行预测：

${\hat{F}}_{2n}=x_{2n}{Q}_i$

${\hat{F}}_{2n}=x_{2n}(Q_i+K_{2n}^{-1}{x}_{2n}^T(F_{2n}-x_{2n}{Q}_i\left)\right)$

其中$K_{2n}=K_i+x_{2n}^T{x}_{2n}$

若t_i1＞t_2n或者t_2n>t_i2，重新创建一个新窗口window(可定义一个不同于已有窗口的下 标)，将数据x_2n放入该窗口内，进行如下处理：

X＝PΔQ^T

F＝PΔ＝PΔQ^T＝XQ

得到该窗口的载荷Q和数据对应的PCA得分F,这里的Q和F与窗口的小标相同。

(步骤7中，t_2n表示时间，如果t_2n在时间t_i1和t_i2之间就加入窗口如果不在t_i1和t_i2之 间就重新创建一个窗口)。

步骤8：针对多批次后续进入的数据，对于batch_n:(x_n1,t_n1),(x_n2,t_n2)...(x_nn,t_nn)重复步 骤7的方法。

步骤9：当所有批次的数据处理完成后，对新采集的近红线光谱实时数据采用与步骤7 相同的方法进行归类到已有窗口中，并预测窗口中的近红线光谱实时数据的主成份得分，具 体运算公式如下

${\hat{F}}_{new}=x_{new}{Q}_i$

${\hat{F}}_{new}=x_{new}\left(Q_i+K_{new}^{-1}{x}_{new}^T\left(F_{new}-x_{new}{Q}_i\right)\right)$

其中$K_{new}=K_i+x_{new}^T{x}_{new}.$

最终可以通过上述算法对新增加的数据进行实时的监控，利用多个窗口观测每个窗口内 数据的一些特征，如特征向量、PCA载荷以及PCA得分，并用增量PCA算法对每个进入窗口 内的新增数据进行修正以及预测它们的载荷和PCA得分。

近红外光谱数据是通过近红外光谱分析仪与SMA光纤反射探头测得，本实施例中采用 12批数据，每批数据的格式是一个时间点对应1557维近红外光谱数据，下表为第一批次奶 酪凝结过程的数据样本。

如果针对第一批次的t＝0,t＝4.8,t＝6.7,t＝12.2分钟的近红外光谱数据作出二维数据图 示化，可以形成如图2所示的二维图像,以便更加直观的分析其特征。

从图中可以清晰的看出牛奶凝结的全过程。横坐标是0-1557维数，纵坐标是近红外光 谱分析仪所测得的近红外光谱数据。图中显示牛奶凝结的过程中，经历了三个阶段的变化, 导致主要的成分在不断的变化。

将牛奶凝结过程的十二批次数据样本约2/3作为训练数据，建立主成分PCA模型，1/3 数据作为检验数据。

首先对1557维数的第一批数据X进行分解，表示为m个向量的外积之和,表示为

$X=f_1{q}_1^T+f_2{q}_2^T+...+f_m{q}_m^T$

式中满足f_i＝Xq_i,i＝1,2,...,m。当X存在一定程度的线性相关时，则X的信息量主要体现在 前q个载荷向量上，因此有

$X=f_1{q}_1^T+f_2{q}_2^T+...+f_q{q}_q^T+E$

前面q个主元代表了数据的主要信息，E表示误差矩阵，可以忽略不计，将数据X用 [f₁,f₂,...f_q]近似表示，即

$X=f_1{q}_1^T+f_2{q}_2^T+...+f_q{q}_q^T$

q＜＜m，从而实现了数据的压缩。PCA算法的关键是求解载荷向量，常用的方法主要是奇异 值分解和特征值分解等。本文所使用的方法主要是基于奇异值分解的方法。

样本数据矩阵X的奇异值分解课表示为

X＝UΕV^T

式中U＝[u₁u₂…u_n]∈R^n×n,V＝[v₁v₂…v_m]∈R^m×m

其中U^TU＝V^TV＝I,π₁＞π₂＞…＞π_m为矩阵X的奇异值。进而：

其中是Σ的特征值，而v_i是Σ的特征向量，因此有Xv_i＝p_iu_i令p_iu_i＝f_i,q_i＝v_i，所以有

$X={\mathrm{UEV}}^T=\underset{i=1}{\overset{m}{\Sigma }}{\pi }_i{u}_i{v}_i^T=\underset{i=1}{\overset{m}{\Sigma }}{f}_i{q}_i^T$

通过上述奇异值分解方法对十二批数据进行PCA降维数，并算出累积变量的解释程度。 样本数据所包含的信息总和为λ₁+λ₂+…+λ_m,第i个主元中包含的信息为λ_i,其对数据的解 释率为λ_i/(λ₁+λ₂+…+λ_m)。一般主元的解释率一般选为85％，这种方法的优点是简单、计 算量较小、可靠。4组批次的数据在主成分1和主成分2中累积变量的解释程度都占90％以 上，主成分1远远大于主成分2。本文使用主成分1作为原始数据的代替，下文中将不再阐 述。

通过奇异值分解的方法对十二批次的近红外光谱数据进行PCA降维处理，最终得到每个 批次各个时间点所对应的PC1得分。从图3中可以看出，牛奶凝结过程主要变化趋势，我们 不能准确的估算出这些波段是属于牛奶凝结过程中的哪一个阶段，即初始酶水解,副酪蛋白 的聚合,凝胶过程。但是我们可以作出PC1得分作为时间的函数，从而进一步的观测数据的 变化。

我们只能从图3中粗略的观测到在7～9分钟时初始酶开始水解，10～45分钟时副酪蛋 白的聚合和凝胶过程。

经过一系列的数据处理，我们将一批次的数据聚类到了17个窗口内，通过每个窗口， 可以观察到每个窗口一些参数，例如每个窗口的载荷向量，每个窗口所包含的数据的数量以 及每个窗口的初始时间点和结束时间点。

窗口中的一些参数如下表。

该显示了第一批次的近红外光谱数据经过窗口后，分成了17类具有相同特性的数据， 在窗口1、窗口2、窗口3内的数据最多，由此可以判断，窗口1为初始酶开始水解阶段， 窗口2为时副酪蛋白的聚合阶段，窗口3为凝胶过程阶段。上述过程和结果只是针对单独批 次进行的处理，可以对单独批次的近红外光谱数据进行观测，并实时监控数据的变化，但是 存在着一些误差。

将第二批次的数据引入第一批次数据的窗口后，经过增量PCA算法的处理，得到第二批 次数据的修正后的PC1得分结果。图4显示了第二批次数据的修正并预测PC1得分和未经修 正PC1的关系，从图中可以看出预测的PC1得分与原来PC1得分相差不大，并有细微的变 化。出现上述结果的原因是第一批的数据对第二批次的数据进行了细微的修正，并且对原来 数据结果并没有造成巨大的偏差。

本文通过R²对上述结果进行检验，这里， $SSE=\underset{i=1}{\overset{n}{\Sigma }}{w}_i{(y_i-{\hat{y}}_i)}^2,SST=\underset{i=1}{\overset{n}{\Sigma }}{w}_i{(y_i-\bar{y})}^2.$

其正常取值范围为0至1，越接近1，表明方程的变量对y的解释能力越强，这个模型 对光谱数据状态判断的越准确，本实验中计算得到R²的值为0.98，说明本发明可以通过对 红外光谱数据的聚类和增量主成分分析，能够很好的判断牛奶在凝结成奶酪过程中的具体凝 聚状态。

从图4中可以看出28分钟和66分钟2个时间点存在着巨大偏差，检验时应当剔除异常 数据,使得检验结果更加准确。从图5中0～60分钟可以看出误差率在5％内,个别数据在10％ 左右。说明该算法能够预测新增数据所表述的凝结状态可行，本发明的方法可以利用已有的 数据对新增的数据进行修正，误差能够控制在10％以内。