基于数据驱动的过程监控和软测量已经成为工业应用研究的重要领域^[1-3]. 当操作条件、外界环境、过程本身固有因素或者特定需求发生变化时，工业生产过程呈现多模态特征^[4]. 为了更好地描述多模态过程的特性，在过程建模前须对多模态数据进行模态识别，即对没有模态指示信息的数据进行类别划分. 模态识别是建立多模态过程监控和软测量的首要前提^[5].

现阶段，针对工业多模态过程数据不同的分布特性，许多聚类方法应用在模态识别中^[6-7]. Du等^[8]采用K-means聚类进行模态识别，提取所有模态的共同特征和特有特征后对多模态数据进行建模. K-means方法在对含有过渡过程数据进行模态识别时结果会变差^[9]. 模糊C均值聚类（fuzzy C-means，FCM）在K-means聚类基础上，引入隶属度进行软划分. Zhang等^[10]在离线建模过程中采用核模糊C均值聚类算法，将多模态数据分为C个模糊组，对各模态建立KPCA模型. 当数据呈现线性可分或超球面形状时，利用FCM进行模态识别可以取得较好的效果. 高斯混合模型（Gaussian mixture model，GMM）可以拟合任意分布的数据. Jiang等^[11]在多模态过程监控中，假设每一模态对应一个高斯分量，利用高斯混合模型进行模态识别. 但工业过程中过渡数据会影响高斯密度函数的参数，导致参数优化收敛的不稳定^[12].

基于窗口的模态识别方法被提出^[13-14]来获取局部数据变化信息，进而更好地表征模态变化关系. Zheng等^[15]计算滑动窗口数据的PLS相似因子对多模态过程稳定模态进行识别；针对过渡模态的识别问题，张淑美等^[16]对窗口均值向量进行聚类，实现稳定模态和过渡模态的识别，解决聚类算法参数选取不准确造成的稳定模态被埋没现象；Song等^[17]针对局部样本间关系，扩展原始样本后引入局部异常因子，提取局部样本密度信息来判断稳定模态和过渡模态；考虑到工业过程不同阶段样本间不同的依赖性，He等^[18]利用基于互信息的相似度识别稳定模态和过渡模态. 上述方法在一定程度实现了稳定模态和过渡模态的识别，但过渡模态是从一种稳定运行模式向另一种稳定运行模式渐进式量变^[19]，在复杂工业过程中过渡模态的过程特性在不同阶段会呈现分段性变化，在识别过渡模态后进行子模态的识别分析是有必要的.

Tan等^[20]认为多模态过程潜在变量相关关系跟随模态变化呈现分段性，利用窗口数据的负载矩阵相似度进行模态识别，在对稳定模态过渡模态进行粗划分后进一步对过渡模态细划分出过渡子模态，所提方法中的相似度基于欧氏距离定义，而PCA提取子空间即负载矩阵属于格拉斯曼流形，是特殊的黎曼空间，属于非欧式空间^[21]. 利用欧氏距离度量非欧式空间的相似度，可能导致时序下的模态识别不准确. 在多模态工业过程中，合适的度量指标对实现模态识别是至关重要的.

本研究提出基于最大均值差异的多模态过渡模态识别方法. 引入最大均值差异作为局部数据分布度量指标，刻画数据概率分布的均值和方差，有效地表征多模态过程数据特性，实现多模态过程下模态的准确识别. 分析局部数据分布差异变化趋势识别模态，削弱噪声对模态识别的影响. 通过仿真实验验证方法的有效性.

在多模态生产过程中，稳定工作模式在同一工作点附近波动，局部数据分布逼近同一概率分布. 当模态发生变化时，过程变量随工作点的改变发生渐变，过程数据的概率分布随之发生变化.

在多模态运行过程中，假设取同一模态下的局部数据 ${{X}} = \{ {{{x}}_1},{{{x}}_2}, \cdot \cdot \cdot, {{{x}}_r}\}$服从分布p，局部数据 ${{Y}} = \{ {{{y}}_1},{{{y}}_2}, \cdot \cdot \cdot, {{{y}}_s}\}$服从分布q，r和s为局部样本个数. 定义分布p和q的差异为 $d(p,q)$. 在同一稳定模态下的局部数据分布差异处于某一变化范围，即 $d(p,q) \leqslant {d_{{\rm{s}}}}$， ${d_{{\rm{s}}}}$为稳定模态局部数据分布最大差异. 在模态发生变化后，原模态局部数据服从分布p，与模态改变后的局部数据服从分布q之间的差异 $d(p,q) > {d_{{\rm{s}}}}$.

根据局部数据服从分布的变化趋势，将多模态过程划分成若干子模态，针对不同子模态的分布特性建立子模型，能更好揭示多模态过程的运行状态和分布变化规律. 基于此，引入衡量p、q分布差异的指标进行模态识别. Gretton等^[22]定义最大均值差异（maximum mean discrepancy，MMD）：

(1) ${\rm{MMD}} [F,p,q] = \mathop {\sup \;} \limits_{f \in F} ({{{E}}_{{{x}} \sim p}}[f({{x}})] - {{{E}}_{{{y}} \sim q}}[f({{y}})]).$

其中：x、y为局部数据样本；F为特征空间映射到实数集的所有函数 $f$集合，MMD的取值取决于使差异最大的投影 $f$. 当F是再生核希尔伯特空间 ${\cal{H}}$的单位球时， ${{\rm{MMD}} ^2}[F,p,q]$表示为

(2) $ \begin{split} {{\rm{MMD}} ^2}[F,p,q] = & {\left[ {\mathop {\sup }\limits_{||f|{|_{\cal{H}}} \leqslant 1} ({{{E}}_p}[f({{x}})] - {{{E}}_q}[f({{y}})])} \right]^2} = \\ & {\left[ {\mathop {\sup }\limits_{||f|{|_{\cal{H}}} \leqslant 1} < {{{\mu }}_p} - {{{\mu }}_q},f{ > _{\cal{H}}}} \right]^2} = \\ & ||{{{\mu }}_p} - {{{\mu }}_q}||_{\cal{H}}^2. \end{split} $

由式(2)可知，MMD刻画了高维特征空间下的任意阶统计量. 针对样本集 ${{X}}$和 ${{Y}}$之间的 ${{\rm{MMD}} ^2}$的经验估计表示为

(3) $ \begin{split} {{\rm{MMD}} ^2}[F,{{X}},{{Y}}] = & ||{{{\mu }}_X} - {{{\mu }}_Y}||_{\cal{H}}^2 = {({{{\mu }}_X} - {{{\mu }}_Y})^{\rm{T}}}({{{\mu }}_X} - {{{\mu }}_Y}) = \\ & {{\mu }}_X^{\rm{T}}{{{\mu }}_X} - {{\mu }}_X^{\rm{T}}{{{\mu }}_Y} - {{\mu }}_Y^{\rm{T}}{{{\mu }}_X} + {{\mu }}_Y^{\rm{T}}{{{\mu }}_Y} = \\ & < {{{\mu }}_X},{{{\mu }}_X} > + < {{{\mu }}_Y},{{{\mu }}_Y} > - \\ & 2 < {{{\mu }}_X},{{{\mu }}_Y} > .\\[-10pt] \end{split} $

将 ${{{\mu }}_X} = \dfrac{1}{r}\displaystyle\sum\nolimits_{i = 1}^r {\phi ({{{x}}_i})}$， ${{{\mu }}_Y} = \dfrac{1}{s}\displaystyle\sum\nolimits_{i = 1}^s {\phi ({{{y}}_i})}$代入，其中 $\phi ( \cdot )$为映射函数，满足 $< \phi ({{x}}),\phi ({{y}}){ > _{\cal{H}}} = k({{{x}},{{y}}})$， $k({{{x}},{{y}}})$为正定核函数. 可得

(4) $\begin{split} {{\rm{MMD}} ^2}[F,{{X}},{{Y}}] = & \frac{1}{{{r^2}}}\sum\limits_{i,j = 1}^r {k({{{x}}_i},{{{x}}_j})} - \frac{2}{{rs}}\sum\limits_{i,j = 1}^{r,s} {k({{{x}}_i},{{{y}}_j})} + \\ & \frac{1}{{{s^2}}}\sum\limits_{i,j = 1}^s {k({{{y}}_i},{{{y}}_j})} . \\[-15pt] \end{split} $

通过引入最大均值差异指标计算相邻局部数据分布间的差异，将分布相近的局部数据识别为同一模态，分布差异较大的局部数据识别为不同模态.

引入稳定模态阈值 $\alpha $和过渡模态阈值 $\beta $. 通过度量稳定模态和过渡模态的数据分布差异，识别过渡模态起始点.

假设多模态过程离线建模数据表示为 ${{{{X}}}} \in {{\rm{{{R}}}}^{n \times m}}$，n为建模样本数，m为过程变量数.

借鉴Tan等^[20]的滑动窗口思想，如图1所示，引入窗宽为H的滑动窗，滑动步长为H，第K个窗口内数据为 ${{X}}_K^{\rm{T}} = \left[ {{{{x}}_{K1}},{{{x}}_{K2}},{{{x}}_{K3}} \cdots, {{{x}}_{KH}}} \right]$，x为 ${m \times 1}$的列向量. 引入大滑动窗H遍历数据，假设第1个窗口为稳定模态，从第2个窗口起，依次计算当前窗口K与前一窗口K−1的MMD指标. 当 ${\rm{MMD}} $ $ ({{X}}_K^{\rm{T}},{{X}}_{K - 1}^{\rm{T}}) > \alpha$时，表明窗口数据分布发生变化；否则K=K+1，继续遍历窗口.

当2个窗口数据分布发生变化时，存在2种情况：1）窗口K数据开始包含过渡模态，K−1窗口依旧处于稳定模态，计算MMD指标大于给定稳定模态阈值. 2）窗口K−1数据后半段进入过渡模态，由于过渡数据较少整个窗口依旧呈现稳定模态分布特性，窗口K数据属于过渡模态，此时计算K窗口与K−1窗口的MMD指标大于稳定模态阈值.

当 ${\rm{MMD}} ({{X}}_K^{\rm{T}},{{X}}_{K - 1}^{\rm{T}}) > \alpha$时，过渡模态可能处于K窗口或K−1窗口. 为了准确判断过渡起始点，从K−2窗口开始分析，建模数据表示为 ${{{X}}^{\rm{T}}}$，引入小滑动窗L遍历数据，小窗口用 ${\tilde {{X}}}_1^{\rm{T}},{\tilde {{X}}}_2^{\rm{T}},{\tilde {{X}}}_3^{\rm{T}} ,\cdots$ 表示，第k个小窗口内数据表示为 ${\tilde {{X}}}_k^{\rm{T}} = \left[ {{{{x}}_{k1}},{{{x}}_{k2}},{{{x}}_{k3}} \cdots, {{{x}}_{kL}}} \right]$. 取第1个小窗口为稳定模态参考窗口，从第2个小窗口起，依次计算当前窗口i与参考窗口的MMD指标. ${\rm{MMD}} ({\tilde {{X}}}_i^{\rm{T}},{\tilde {{X}}}_1^{\rm{T}}) > \beta $表明窗口数据与稳定模态数据分布不一致，该窗口为过渡模态起始窗口.

将经过识别的过渡模态起始窗口作为过渡子模态参考窗口，计算第j个窗口数据与参考窗口数据分布的最大均值差异 ${\rm{MMD}} ({\tilde {{X}}}_j^{\rm{T}},{\tilde {{X}}}_{k1}^{\rm{T}})$， ${\rm{MMD}} $ $ ({\tilde {{X}}}_j^{\rm{T}},{\tilde {{X}}}_{k1}^{\rm{T}}) < \beta $，表明当前窗口与参考窗口属于同一过渡子模态；反之，当前窗口与参考窗口属于不同过渡子模态，第k1个窗口到第j−1个窗口数据为第1个过渡子模态. 将第j个窗口作为新的过渡子模态参考窗口 ${\tilde {{X}}}_{k2}^{\rm{T}}$，按上述步骤继续进行分析，依次完成对过渡子模态的识别. 过程一直持续到属于某过渡子模态的连续窗口样本数大于最短稳定模态长度D（D根据经验知识和实际稳定模态运行最短时间选定），过渡模态结束，连续窗口为下一个稳定模态.根据本节策略，所有过渡子模态构成过渡模态，方法完成对过渡模态及过渡子模态的识别.

本方法中引入2个阈值，分别为稳定模态阈值 $\alpha $和过渡模态阈值 $\beta $，阈值的选取会影响模态识别结果. 稳定模态阈值 $\alpha $选取偏小，同一稳定模态中数据分布相近的窗口会被识别为过渡模态窗口；稳定模态阈值 $\alpha $选取偏大，过渡模态窗口未被识别. 过渡模态阈值 $\beta $的选取偏大，过渡模态中不同分布特性的数据会被识别为同一过渡子模态，后续对不同分布的数据建立软测量模型会降低模型预测精度；过渡模态阈值 $\beta $选取偏小，对过渡模态数据分布差异识别更敏感，导致识别出更多的过渡子模态，建模样本变少，模型更易受噪声干扰. 在实际应用中，通过试验选取恰当的阈值.

大滑动窗口切割数据执行 $n/H$次，计算MMD指标的时间复杂度为 $O(4{H^2})$^[22]，当比较计算两两窗口数据MMD时，执行 $(n/H) {(H + H)^2}$次，根据阈值确定过渡模态起始窗口执行 $n/H$次，对剩余数据分析，小滑动窗口切割数据执行 $\tilde n/L$次，其中 $\tilde n$为 ${\tilde {{X}}}$的样本数量. 当确定过渡子模态计算MMD时，执行 $(\tilde n/L) {(L + L)^2}$次. 方法复杂度为 $O(4Hn + n/H + $ $ 4L\tilde n + \tilde n/L)$. 相较于利用负载矩阵相似度（loading matrix similarity，LMS）的模态识别方法，由于求解H阶矩阵的负载矩阵的时间复杂度为 $O({H^3})$^[23]，大于计算MMD指标的时间复杂度，当面对高维数据时本研究方法更高效. 基于聚类的模态识别方法由于实际迭代次数不确定，时间复杂度变化不稳定.

仿真包含5个监控变量 ${{{y}}_1}{\text{、}}{{{y}}_2}{\text{、}}{{{y}}_3}{\text{、}}{{{y}}_4}{\text{、}}{{{y}}_5}$， ${{{x}}_1}{\text{、}}{{{x}}_2}$作为输入变量. 包含渐变过渡趋势的多模态过程数据方程为

(5) ${{{y}}_1} = 0.576\;8{{{x}}_1} + 0.376\;6{{{x}}_2} + {{{e}}_1},$

(6) ${{{y}}_2} = 0.738\;2{{x}}_1^2 + 0.056\;6{{{x}}_2} + {{{e}}_2},$

(7) ${{{y}}_3} = 0.829\;1{{{x}}_1} + 0.400\;9{{x}}_2^2 + {{{e}}_3},$

(8) ${{{y}}_4} = 0.651\;9{{{x}}_1}{{{x}}_2} + 0.207\;0{{{x}}_2} + {{{e}}_4},$

(9) ${{{y}}_5} = 0.397\;2{{{x}}_1} + 0.804\;5{{{x}}_2} + {{{e}}_5}.$

式中： $ {{e}}_{1}$、 ${{e}}_{2}$、 ${{e}}_{3}$、 ${{e}}_{4}$、 ${{e}}_{5}$为均值为0，标准差为0.01的白噪声. 首先令输入变量 ${{{x}}_1}$服从均值为5，方差为0.36的正态分布 $N(5,0.36)$，输入变量 ${{{x}}_2}$服从均值为20，方差为0.49的正态分布 $N(20,0.49)$，按表1所示定义稳定模态P、Q. 依次在第401、501、601、701个采样点处对源 ${{{x}}_1}$引入期望值的阶跃，产生4种过渡子模态为PQ_1、PQ_2、PQ_3、PQ_4. 仿真共采集1 200个样本，其中稳定模态P、过渡模态PQ、稳定模态Q各采集400个样本. 如图2所示为5个输出变量的变化曲线.

利用MMD度量识别过渡模态时，选取 $H = 100$， $L = 25$， $\alpha = 0.3$， $\beta = 0.2$. 引入K-means、FCM、GMM算法作为对比方法，模态识别结果如图3所示。图中，M为识别模态标签. 由图3可知，基于聚类的模态识别方法将同一模态中波动的部分识别为2类不同模态，这是数据波动性和噪声对计算均值产生的影响. 基于负载矩阵相似度和最大均值差异的模态识别结果如图4（a）、（b）所示，图4（b）清晰地识别出400到800样本点间的4个过渡子模态，由结果可知基于MMD的模态识别结果更准确. 分析原因，以仿真数据中每100个样本构成的窗口数据为分析单元，第1个窗口（即1~100样本）为参考窗口，计算各窗口与参考窗口的负载矩阵相似度和MMD指标， 如图4（c）、（d）所示。图中， ${\varphi _{\rm{LMS} }}$为负载矩阵相似度，W为窗口序号， ${\varphi _{\rm{MMD}}}$为归一化后的MMD指标. 由图4（c）、（d）可知，MMD指标随着期望的增大逐渐增大，负载矩阵相似度随着期望的增大未能呈现逐渐减小的趋势，说明负载矩阵相似度只在一定程度上度量了数据模态变化趋势，未能准确划分出模态；MMD指标可以更有效捕捉到时序上的模态变化趋势，基于MMD的模态识别方法能更准确地识别出过渡模态.

本研究的实验数据来源于田纳西-伊斯曼（Tennessee Eastman, TE）过程的仿真. 以反应器压力设定值变化的过程为对象，选择TE过程中15个过程连续变量作为被监控变量^[20].

模拟40 h的仿真过程，压力初始值为2 800 kPa，过程运行至 $T = 20.824$ h，将反应器压力设定值减少为2 705 kPa，采集新模态过程数据. 仿真多模态过程数据共4 000个样本，包含2种稳定模态和1种过渡模态. 仿真过程变量V₁~V₁₅变化如图5所示.

取 $H = 100$， $L = 30$，经试验确定 $\alpha $和 $\beta $，取 $\alpha = $ $ 0.3$， $\beta = 0.06$. 最终模态识别结果如下：P（1~2 080样本点），PQ_1（2 081~2 110样本点），PQ_2（2 111~2 140样本点），PQ_3（2 141~2 170样本点），PQ_4（2 171~2 260样本点），Q（2 261~4 000样本点）. 结合图5分析可知，本方法识别出第2 081~2 110样本点作为第1个过渡子模态，在该过程下反应器温度由于压力设定值突然降低而降低. 在第2111~2170样本点上，很多过程变量正处在一定的趋势变化中，例如反应器、分离器和汽提塔温度处于持续升高的过程中，而分离器和汽提塔压力处于持续减小的过程. 基于MMD的模态识别方法在捕捉到变量分布变化趋势基础上，将该过程划分为分布差异较大的2个过渡子模态. 在第2 170~2 260样本点上，各变量在偏离原稳态工作点后逐渐达到新稳态. 第2 261样本点之后，局部数据分布之间的差异较小，步入下一个稳定模态.

为了验证方法的有效性，利用不同模态识别方法进行模态识别，建立多模型软测量实验，具体步骤如下。1）选择15 个过程变量作为软测量模型输入变量，将TE过程流9中的3个成分浓度（变量A、B、C）作为输出变量.2）将TE仿真实验中获取的4 000个离线样本进行标准化处理，选取包含完整过渡过程的中间2 800 个样本作为训练样本集. 在其他仿真条件不变的情况下，改变压力设定值至2700 kPa的过程作为工况1，改变压力设定值至2710 kPa的过程作为工况2，从2种工况中获取包含过渡过程的仿真数据作为测试样本集.3）在训练集上，将K-means、FCM、GMM、基于LMS模态识别方法以及本文基于MMD的模态识别方法划分后的数据子集，分别建立局部PLSR模型.4）针对测试样本，选择最相近的数据子集建立的软测量模型，对变量A、B、C进行预测.

利用均方根误差（root mean square error，RMSE）指标定量分析不同模态识别方法下的多模型软测量预测结果，如表 2所示.

由表2可知，基于MMD的模态识别方法下的多模型软测量预测效果更好. 基于聚类的模态识别方法可能将过渡过程中的分布接近但标签不一致的不同子阶段识别为同一子模态，而基于MMD的模态识别方法是根据时序划分过渡子模态，使得样本的标签分布基本一致，以提高后续建模精度. 相较于负载矩阵相似度划分模态，通过MMD直接度量分布差异，识别出的子模态样本间分布差异更明显，待测样本更易选择最合适的软测量模型进行预测，提高预测精度.

针对现有模态识别方法在多模态过渡数据上难以准确划分的问题，提出基于最大均值差异的多模态过程中的过渡模态识别方法. 在数值仿真中，所提方法准确识别出过渡子模态，实验结果表明，最大均值差异相较负载矩阵相似度，更能有效地刻画局部数据分布差异，提高模态识别精度，相比聚类模态识别方法，基于最大均值差异的模态识别方法可以更准确地识别模态. 在TE仿真实验中，在实现更准确识别过渡模态的基础上，将复杂多模态数据按时序识别为分布差异明显的子模态后建模，克服了过渡过程中具有相似数据分布的样本中标签分布不一致的问题，在过渡阶段的建模预测中具有较好的优势，证明本研究方法的有效性.本研究方法通过对数据的反复试验确定滑动窗的窗宽取值，在多模态过程中能否针对数据分布的变化自适应确定窗宽，需要进一步研究。