1 前言

间歇生产方式在现代工业界所占比重越来越高，其过程监控和故障诊断也引起广泛关注^[1]。长期以来，过程监控的主流统计方法^[2-5]如主元分析(principal component analysis, PCA)、偏最小二乘(partial least squares, PLS)、独立元分析(independent component analysis, ICA)、典型相关分析(canonical correlation analysis, CCA)等在应用于间歇过程故障诊断时，都是采用将三维多批次历史数据展开成二维的多向(multi-way)处理方法，即使针对间歇过程的不同特点涌现出一系列改进算法^[6-14]，也大都对数据进行展开处理，如按批次展开、按变量展开成二维矩阵的方式等。此类多向展开的处理方法将三维多批次数据中每一批次的二维矩阵展开变成向量看待，势必会破坏原有数据的相关性，导致部分信息的丢失和特征不完全提取。

近年来，一些基于张量的数据降维方法，如张量PCA (tensor PCA, TPCA)^[15-16]、张量邻域保持嵌入算法(tensor neighborhood preserving embedding, TNPE)^[17]等，使用张量模式的特征提取方法对张量形式的数据进行特征提取，保证了张量形式的特征空间不被破坏，张量分析逐渐在计算机视觉、图像处理、模式识别等领域成为研究热点^[18]。本文借鉴以上方法中张量分析的思想，对间歇过程在线监控时，将三维数据的每一批次看作一个二阶张量，应用TPCA算法对三维数据直接提取特征信息，不破坏三维样本中每一批次的原始矩阵空间结构，并定义R²统计量和SPE统计量建立故障检测模型；另外，为能更及时检测出间歇过程中常出现的渐变故障(如间歇精馏过程中催化剂活性的缓慢变化、间歇反应釜换热器由于长时间使用结垢现象越来越严重而导致换热性能的逐渐下降)，本文利用累积和控制图(cumulative sum，CUSUM)分析思想，将其应用于TPCA算法中，提出累加和的张量主元分析(summed tensor principal component analysis, STPCA)，将三维历史数据按一定的滑动窗口累加求和，累积历史信息，再做TPCA算法处理，建立故障检测模型，从而提高对渐变故障在初始出现阶段检测的灵敏度。

2 TPCA方法

假设收集到的间歇过程历史数据矩阵为$\mathit{\boldsymbol{X}}(M \times K \times J) = \left\{ {{\mathit{\boldsymbol{X}}} \in {R} \otimes {R}, \;m = 1, 2, \ldots , M} \right\}$，$\mathit{\boldsymbol{X}}$中含有$M$个批次，每批次数据${\mathit{\boldsymbol{X}}_m}$分别有$J$个观测变量和$K$次采样。TPCA算法首先将每批次${\mathit{\boldsymbol{X}}_m}$看成一个${R^K} \otimes {R^J}$的二阶张量子空间，与PCA思想类似，TPCA目的是寻找两个低秩投影矩阵$\mathit{\boldsymbol{U}} \in {R^K} \otimes {R^{K'}}$和$\mathit{\boldsymbol{V}} \in {R^J} \otimes {R^{J'}}$，由于$K' \ll K$，$J' \ll J$，从而把原始数据X投影到最大方差的低维张量特征空间$\mathit{\boldsymbol{Y}}(M \times K' \times J') = \left\{ {{\mathit{\boldsymbol{Y}}_m} \in {R^{K'}} \otimes {R^{J'}}, \;{\mathit{\boldsymbol{Y}}_m} = {\mathit{\boldsymbol{U}}^{\rm{T}}}{\mathit{\boldsymbol{X}}_m}\mathit{\boldsymbol{V}}} \right\}$，以使重构误差减少到最小，其目标函数遵循式(1)^[16]。

$\arg \mathop {\min }\limits_{\mathit{\boldsymbol{U}}, \mathit{\boldsymbol{V}}} \sum\limits_{m = 1}^M {\left\| {\mathit{\boldsymbol{U}}{\mathit{\boldsymbol{Y}}_m}{\mathit{\boldsymbol{V}}^{\rm{T}}} - {\mathit{\boldsymbol{X}}_m}} \right\|_{\rm{F}}^2} \\ {\rm{s.t.}}\quad {\mathit{\boldsymbol{U}}^{\rm{T}}}\mathit{\boldsymbol{U}} = I\quad {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{V}} = I$

(1)

可通过以下迭代方法实现式(1)中U、V的求解：

Step1：初始化${\mathit{\boldsymbol{U}}_1} = I \in {R^{K \times K}}$、${\mathit{\boldsymbol{V}}_1} = I \in {R^{J \times J}}$、最大迭代次数${N_{\max }}$。

Step2：设置循环变量$n$初始值为1，进入循环迭代运算。

Step3：令${\mathit{\boldsymbol{Y}}_{1m}} = {\mathit{\boldsymbol{U}}_n}^{\rm{T}}{\mathit{\boldsymbol{X}}_m}$，求${Y_1}$的协方差阵${\mathit{\boldsymbol{C}}_{{{\rm{Y}}_1}}} = \frac{1}{M}\sum\limits_{m = 1}^M {{\mathit{\boldsymbol{Y}}_{1m}}^{\rm{T}}{\mathit{\boldsymbol{Y}}_{1m}}} $。

Step4：对${\mathit{\boldsymbol{C}}_{{Y_1}}}$进行特征值分解${\mathit{\boldsymbol{C}}_{{Y_1}}} = {\mathit{\boldsymbol{V}}_{n + 1}}\mathit{\boldsymbol{ \boldsymbol{\varLambda}}} {\mathit{\boldsymbol{V}}_{n + 1}}^{\rm{T}}$，得到特征值矩阵$\mathit{\boldsymbol{ \boldsymbol{\varLambda}}} $和特征向量矩阵${\mathit{\boldsymbol{V}}_{n + 1}}$，$\Lambda = {\rm{diag}}({\lambda _1}, {\lambda _2}, \ldots , {\lambda _{J'}})$只选取前$J'$个最大的特征值组成，${\mathit{\boldsymbol{V}}_{n + 1}}$则是由前$J'$个特征值对应的特征向量矩阵。

Step5：令${\mathit{\boldsymbol{Y}}_2}_m = {\mathit{\boldsymbol{X}}_m}{\mathit{\boldsymbol{V}}_n}$，求${\mathit{\boldsymbol{Y}}_2}$的协方差阵${\mathit{\boldsymbol{C}}_{{{\rm{Y}}_2}}} = \frac{1}{M}\sum\limits_{m = 1}^M {{\mathit{\boldsymbol{Y}}_{2m}}^{\rm{T}}{\mathit{\boldsymbol{Y}}_{2m}}} $。

Step6：特征值分解${\mathit{\boldsymbol{C}}_{{{\rm{Y}}_2}}} = {\mathit{\boldsymbol{U}}_{n + 1}}\mathit{\boldsymbol{ \boldsymbol{\varLambda}}} {\mathit{\boldsymbol{U}}_{n + 1}}^{\rm{T}}$，得到特征值矩阵$\mathit{\boldsymbol{ \boldsymbol{\varLambda}}} $和特征向量矩阵${\mathit{\boldsymbol{U}}_{n + 1}}$，$\Lambda = {\rm{diag}}({\lambda _1}, {\lambda _2}, \ldots , {\lambda _{K'}})$保留前$K'$个最大的特征值，${\mathit{\boldsymbol{U}}_{n + 1}}$是由前$K'$个特征值对应的特征向量矩阵。

Step7：若$\left\| {{\mathit{\boldsymbol{U}}_{n + 1}} - {\mathit{\boldsymbol{U}}_n}} \right\|_F^2 < \mathit{\boldsymbol{\varepsilon}} $和$\left\| {{\mathit{\boldsymbol{V}}_{n + 1}} - {\mathit{\boldsymbol{V}}_n}} \right\|_{\rm{F}}^2 < \mathit{\boldsymbol{\varepsilon}} $成立，则结束，否则$n = n + 1$，并跳转到Step3。

Step8：输出最终的$\mathit{\boldsymbol{U}} \in {R^{K \times K'}}$、$\mathit{\boldsymbol{V}} \in {R^{J \times J'}}$。

注意，TPCA变换处理的数据X是经过标准化处理的，其处理过程如图 1所示，第1步将初始数据$\underline{\mathit{\boldsymbol{X}}}(M \times K \times J)$展开成大小为$M \times KJ$的二维矩阵，沿着每列进行均值为0方差为1的归一化处理，第2步再恢复成三维结构数据$\mathit{\boldsymbol{X}}(M \times K \times J)$。此处的展开仅用于数据的预处理，应用TPCA变换时使用的仍是三维矩阵$\mathit{\boldsymbol{X}}(M \times K \times J)$。因此，预处理后的三维数据$\mathit{\boldsymbol{X}}(M \times K \times J)$通过以上Step1~8的迭代算法获得U和V后，经过映射公式${\mathit{\boldsymbol{Y}}_m} = {\mathit{\boldsymbol{U}}^{\rm{T}}}{\mathit{\boldsymbol{X}}_m}\mathit{\boldsymbol{V}}$，则可被投影到低维张量特征空间得到$\mathit{\boldsymbol{Y}}(M \times K' \times J') = \left\{ {{\mathit{\boldsymbol{Y}}_m} \in K' \times J', \;\;m = 1, 2, \ldots , M} \right\}$。

3 基于STPCA方法的间歇过程故障诊断 3.1 STPCA算法

各类间歇过程中的渐变故障很常见，渐变故障是某变量操作值缓慢偏离正常变化区间、不断累积的过程，故障初期变化量较小，很多故障检测方法对该类故障存在预报大大延迟的问题。为提高检测故障灵敏度，本文利用PAGE等^[19]提出的累积和(CUSUM)控制图思想，提出一种基于累加和的张量主元分析算法(STPCA)，将经过预处理后的正常批次数据$X(M \times K \times J)$按照每批次内滑动窗口大小进行历史信息的叠加，数据处理公式如下：

${\rm{summed}}\_{\mathit{\boldsymbol{x}}_m}(k) = \left\{ \begin{gathered} \sum\limits_{{k_1} = 1}^k {{\mathit{\boldsymbol{x}}_m}({k_1})} \quad \quad k \leqslant h \\ \sum\limits_{{k_1} = k - h + 1}^k {{\mathit{\boldsymbol{x}}_m}({k_1})} \quad \quad k > h \\ \end{gathered} \right.$

(2)

${\rm{summed}}\_{\mathit{\boldsymbol{x}}_m}(k)$表示第$m$批次中按照式(2)累加和后的第$k$采样时刻观测值，$h$为累加求和的宽度，构成移动的滑动窗口，当$k > h$，每一时刻的数据值都由前$h$个多变量样本求和得到，从而得到求和后的${\rm{summed}}\_\mathit{\boldsymbol{X}}(M \times K \times J) = \{ {\rm{summed}}\_{\mathit{\boldsymbol{X}}_m} \in {R^K} \otimes {R^J}, {{\kern 1pt} _{}}m = 1, 2, \ldots , M\} $，然后用第1节中TPCA算法对其进行特征提取，求得变换矩阵U和V，为下一步故障诊断打下基础。

3.2 基于STPCA算法的在线监控

在线监控过程主要包含2个阶段。

第1阶段：离线模型的建立

① 将三维训练样本$\underline{\mathit{\boldsymbol{X}}} (M \times K \times J)$按照图 1进行标准化处理，得到$X(M \times K \times J)$。

② 对$\mathit{\boldsymbol{X}}(M \times K \times J)$按式(2)进行累加和处理，获得${\rm{summed}}\_\mathit{\boldsymbol{X}}(M \times K \times J)$。

③ 运用TPCA算法，确定投影矩阵$\mathit{\boldsymbol{U}} \in {R^{K \times K'}}$和$\mathit{\boldsymbol{V}} \in {R^{J \times J'}}$，计算经过投影后的低维张量特征空间数据$\mathit{\boldsymbol{Y}} = \left\{ {{\mathit{\boldsymbol{Y}}_m} \in {R^{K'}} \otimes {R^{J'}}, \;m = 1, 2, \ldots , M} \right\}$，${\mathit{\boldsymbol{Y}}_m}$是累加和处理后${\rm{summed}}\_{\mathit{\boldsymbol{X}}_m}$的低维投影，变换公式为：${\mathit{\boldsymbol{Y}}_m} = {\mathit{\boldsymbol{U}}^T} * {\rm{summed}}\_{\mathit{\boldsymbol{X}}_m} * \mathit{\boldsymbol{V}}\;$。

④ 建立故障检测模型。

对基于TPCA和STPCA 2种算法的间歇过程故障检测，本文都采用R²和SPE统计量。

R²是基于支持向量数据域描述(support vector domain description，SVDD)^[20]方法建立的统计量，表示每批次${\mathit{\boldsymbol{X}}_m}$的投影${\mathit{\boldsymbol{Y}}_m}$到特征空间中心的距离，其计算式如下^[21]：

${{\rm{R}}^2}(m) = {\rm{tr}}({\mathit{\boldsymbol{Y}}_m}{\mathit{\boldsymbol{Y}}_m}^{\rm{T}}) - 2\sum\limits_{i = 1}^M {{\alpha _i}{\rm{tr}}({\mathit{\boldsymbol{Y}}_m}{\mathit{\boldsymbol{Y}}_i}^{\rm{T}})} + \sum\limits_{i = 1}^M {\sum\limits_{j = 1}^M {{\alpha _i}{\alpha _j}{\rm{tr}}({\mathit{\boldsymbol{Y}}_i}{\mathit{\boldsymbol{Y}}_j}^{\rm{T}})} } $

(3)

拉格朗日乘子${\alpha _i}$通过下式的最大化问题求解^[21]：

(4)

$C$为惩罚因子，可取1。

SPE统计量用来测量残差空间的变化。先通过U、V重构${\rm{summed}}\_{\mathit{\boldsymbol{X}}_m}$的估计值，得到残差矩阵${E_m} = {\rm{summed}}\_{\mathit{\boldsymbol{X}}_m} - \mathit{\boldsymbol{U}}{\mathit{\boldsymbol{Y}}_m}{\mathit{\boldsymbol{V}}^{\rm{T}}}$，从而得SPE统计量的计算式：

${\rm{SPE}}{_m}(k) = {\mathit{\boldsymbol{e}}_m}(k){\mathit{\boldsymbol{e}}_m}{(k)^{\rm{T}}}$

(5)

R²和SPE统计量皆可通过核密度估计法^[22]确定各自95%或99%的置信限(也称控制限、阈值)，在每个采样时刻$k$分别利用核密度估计法计算对应采样点所有批次统计量的阈值，因此可得到动态变化的控制限，当统计量超出各自控制限所判定的正常变化范围，则意味着故障的出现。

第2阶段：在线故障检测

① 采集新批次第$k$采样时刻数据${\mathit{\boldsymbol{x}}_{{\rm{new}}, k}} \in {R^{1 \times J}}$，根据离线建模标准化处理时得到的均值和方差进行预处理。

② 对预处理后的数据${\mathit{\boldsymbol{x}}_{{\rm{new}}, k}}$按照式(2)进行累加求和处理，得到${\rm{summed}}\_{\mathit{\boldsymbol{x}}_{{\rm{new}}, k}}$。

③ 计算新的投影矩阵和残差：

${\mathit{\boldsymbol{Y}}_{{\rm{new}}, k}} = {\mathit{\boldsymbol{u}}_k}^{\rm{T}}*{\rm{summed}}\_{\mathit{\boldsymbol{x}}_{{\rm{new}}, k}}*\mathit{\boldsymbol{V}}$

(6)

${\mathit{\boldsymbol{e}}_{{\rm{new}}, k}} = {\rm{summed}}\_{\mathit{\boldsymbol{x}}_{{\rm{new}}, k}} - {\mathit{\boldsymbol{u}}_k} * {\mathit{\boldsymbol{Y}}_{{\rm{new}}, k}} * {\mathit{\boldsymbol{V}}^T}$

(7)

其中${\mathit{\boldsymbol{u}}_k}$表示矩阵U中的第$k$行向量。

④ 参考式(3)和(5)分别计算新的在线R²和SPE统计量，得到

${{\rm{R}}_{{\rm{new}}}}^2(k) = {\rm{tr}}({Y_{{\rm{new}}, k}}{\mathit{\boldsymbol{Y}}_{{\rm{new}}, k}}^{\rm{T}}) - 2\sum\limits_{i = 1}^M {{\alpha _i}{\rm{tr}}({Y_{{\rm{new}}, k}}{\mathit{\boldsymbol{Y}}_i}^{\rm{T}})} + \sum\limits_{i = 1}^M {\sum\limits_{j = 1}^M {{\alpha _i}{\alpha _j}{\rm{tr}}({\mathit{\boldsymbol{Y}}_i}{\mathit{\boldsymbol{Y}}_j}^{\rm{T}})} } $

(8)

${\rm{SPE}}{_{{\rm{new}}}}(k) = {\mathit{\boldsymbol{e}}_{{\rm{new}}, k}}{\mathit{\boldsymbol{e}}_{{\rm{new}}, k}}^{\rm{T}}$

(9)

并判断统计量是否偏离各自阈值，以确定间歇过程运行状况。

4 仿真结果分析

分批补料的盘尼西林发酵过程是典型的间歇式生产过程，本文基于BIROL等^[23]开发的标准仿真软件Pensim 2.0产生的盘尼西林发酵数据进行在线监控算法仿真研究。本仿真由该软件在正常操作范围内对初始条件加入一定的微小偏移和适当噪声生成30批次数据作为训练样本，每批次含11个观测变量，反应时间相等，均为400 h，采样间隔为0.5 h，从而得到$\underline{\mathit{\boldsymbol{X}}} (30 \times 800 \times 11)$，用于建立MPCA、TPCA、STPCA监控模型。另外，利用Pensim2.0产生6个故障批次，表 1是3种跳变故障，表 2是渐变类型故障。

建模时，3种方法都需要对数据进行标准化，不同的是，MPCA方法如图 1中Step 1按批次展开成二维，按列作归一化处理，得到二维矩阵$\mathit{\boldsymbol{X}}(30 \times {8_{}}800)$，即$\mathit{\boldsymbol{X}}(M \times KJ)$，在此基础上建立MPCA模型。TPCA和STPCA算法则需将标准化后的二维矩阵$\mathit{\boldsymbol{X}}(30 \times 8800)$再按图 1中Step 2还原成$\mathit{\boldsymbol{X}}(30 \times 800 \times 11)$，在保持和不破坏原来每批次数据结构基础上建立各自的监控模型。因此，传统PCA算法会把每批次的数据${\mathit{\boldsymbol{X}}_m} \in {R^K} \otimes {R^J}$转换为向量${\mathit{\boldsymbol{x}}_m} \in {R^{KJ}}$看待，特征值分解的矩阵大小是$KJ \times KJ$，基于张量的TPCA方法其特征值分解的矩阵大小仅为$K \times K$或$J \times J$，在特征值分解计算复杂度上要更低。另外，STPCA相比TPCA方法在数据处理上多一步累加和处理，以增大对渐变故障检测的灵敏度。最后的监控模型所用到统计量也有差异：MPCA使用${T^2}$和SPE，TPCA和STPCA则运用本文定义的R²和SPE统计量。

当新批次中出现故障1、2、3跳变故障时，基于MPCA的在线故障检测结果如图 2、4、6所示，MPCA的SPE统计量能准确预报故障的发生(140~300采样点)，但在非故障时段有较高误报，从图中可明显看到有多点超出控制限，特别是对故障3达到了5.12%的误报率(见表 3)。相对而言，图 3、5、7所示TPCA故障检测结果中，SPE统计图既能准确检出故障，又具有较低误报率，由表 3可见，对前两个故障的误报率为0，仅对故障3有约1.63%的误报，但仍低于MPCA方法的结果。与SPE统计量相比，MPCA方法的${T^2}$统计图检测效果较差，只能对故障1有效，对故障2和3完全不敏感，无任何超限警报；反观TPCA方法的R²统计图，在采样点140~300仅无法检测出故障3，能完全准确检测到故障2，对故障1仅有少量点的漏报。从误报率情况分析，TPCA的R²统计量也优于MPCA的${T^2}$统计量。总体分析，TPCA两种统计图对3个故障的检出率高于MPCA方法，且具有较低的误报率。可见，将利用张量分解直接处理三维数据的TPCA方法引入到间歇过程故障检测中，构造R²和SPE统计量进行实时监控是可行并具有一定效果的。

为提高检测渐变故障的灵敏度，本文在TPCA基础上提出基于STPCA的间歇过程故障诊断方法，同样使用R²和SPE统计图进行故障分析。仿真实验设置了表 2所述3种故障，分别是通风速率、搅拌功率、底物流加速率从第200采样点以一定速度呈现斜坡变化，模拟间歇过程中可能出现的渐变类故障。图 8到图 13反映了MPCA、STPCA方法对3个故障批次的检测效果。

从图 8、10和12所示结果观察，MPCA的${T^2}$和SPE统计量均在一定时刻越过阈值不断增加，表明渐变型故障的发生；同样，图 9、11、13中STPCA方法的统计结果图也给出了故障警报。但明显的是，两种方法的SPE统计图比基于得分的${T^2}$和基于特征空间的R²统计量趋于更灵敏，给出故障预报要早。为直观对比不同方法的检测效果优劣，本文通过表 4列出两种方法各统计量所确定的故障出现时刻，并假设统计量连续8个采样点超过控制限则认定该采样时刻为故障检测点。对比SPE统计量，针对渐变故障4、5、6，STPCA分别领先MPCA方法9个、28个、30个采样点提前给出故障指示；STPCA的R²统计图检测故障的能力虽逊于SPE统计量，却仍比MPCA的${T^2}$统计量分别早117、336、56个采样时刻超出控制限。综合看来，两种方法检测渐变故障时，相对故障实际发生时刻都有或多或少的滞后，原因在于故障变量刚开始是缓慢变化的，偏离正常范围需要一定的时间间隔，只有当变量值渐渐漂移累积超出一定量后，监测系统才能判断过程出现了非正常的变化。但STPCA方法将历史批次数据进行叠加处理，放大故障变量的初始偏移，并在张量空间进行特征值分解、不破坏原有数据结构基础上得到监控统计量，结合以上结果图和表 4，不难看出该方法能更早检测出渐变故障，体现该方法的优越性。

需要说明的是，STPCA方法更适用于渐变故障存在的场合，尤其是故障出现早期，通过CUSUM累积过程中的微小变化，提高对渐进、缓慢变化类故障初期较小故障偏移的检测。同时，该方法中累加求和宽度$h$的取值与检测效果关系密切，若$h$取值太小，利用CUSUM累积的历史信息不够丰富，对渐变故障出现初期检测准确率会较低，反之一味追求快速检测出故障而将$h$取得过大，会因带入部分噪声信息而增加误报率，导致故障检测效果下降。目前总体上，$h$大小的选择没有统一标准，多是根据经验或交叉验证方法确定。本文通过反复试验，选择适当的$h$值，既能提高对渐变故障的检出率，较早发出故障报警，又防止出现较高误报的情况。另外，STPCA算法在故障检测时，对新批次数据引入CUSUM的处理，不可避免会产生故障检测延时，因此，本方法并不适用于故障检测实时性要求较高的场合。

5 结论

传统将三维数据变成二维的间歇过程故障诊断方法，破坏数据原有的空间结构和相关特性，必然会对过程监控效果产生不利影响。本文将基于张量分解的TPCA算法引入到间歇过程监控中，直接将三维数据作为张量目标进行运算，保留了数据完整性，定义R²和SPE统计量建立故障检测模型。经仿真验证，基于TPCA的间歇过程故障诊断方法适合检测突变故障，检测效果较好。针对间歇过程中常发生的渐变故障，本文提出STPCA方法，通过对三维数据多个连续时刻的采样值求和，不断放大历史信息，同时又保证不破坏三维原始数据结构，因此在检测渐变类型故障的仿真实验中，STPCA方法有更好效果，提高故障检测灵敏度，相应降低了漏检率。

符号说明：

C  —惩罚因子    summed_x_m(k)  —样本训练集中，累加求和处理后的第m批次第k次观测向量

e_m(k)  —第m批次的第k个残差向量                

e_{new, k}  —新的残差向量  summed_x_new, _k  —新批次数据中，累加求和后的第k次观测向量

h  —累加求和宽度  U  —X_m的左酉投影矩阵

J'  —求右酉矩阵V对应保留的最大特征值个数  V  —X_m的右酉投影矩阵

K'  —求左酉矩阵U对应保留的最大特征值个数  X_m  —第m批次二阶张量历史数据

R²(m)  — R²统计量：Y_m到特征空间中心的距离  Y_m  —X_m的低秩投影特征空间

SPE_m}(k)  —第m批次的第k个SPE统计量  Y_{mew, k}  —summed_x_new, _k到特征空间的投影

                a_i  —拉格朗日乘子