数控机床作为一种复杂的机电设备，其结构复杂，故障隐蔽性较高，难以发觉^[1]。数控机床状态预测有利于及时发现数控机床健康状况，进而对机床进行检修，有效地控制数控机床受损程度，降低生产成本。因此，对数控机床状态预测的研究具有重要的意义。

传统数控机床状态预测方法有神经网络、支持向量机、粒子群算法等。Zhang等^[2]提出串行灰色神经网络和并联灰色神经网络方法预测数控机床的热误差。Sheng等^[3]提出Back-Propagation神经网络模型预测数控机床刀具的状态监测和刀具磨损量的预测。杜柳青等^[4]通过构建一种小波神经网络模型，实现对数控机床运动精度的预测。石宁等^[5]结合遗传算法和支持向量机的优势，建立改进支持向量机预测模型对机械状态进行振动趋势预测。许志军^[6]提出采用粒子群算法优化支持向量机进行数控机床状态预测。针对传统的预测方法，神经网络存在容易陷入局部极小值、收敛速度太慢等问题，支持向量机中训练参数的选择对其预测性能有很大影响，且没有考虑机床状态时序性，限制传统预测方法的实用性和准确性。

时间序列预测是目前应用最广的状态预测方法^[7]，针对数控机床状态的时序性，张海波等^[8]使用随机事件序列法建立数控机床故障频率预测模型，以保证故障频率预测值的准确度。李永祥等^[9]提出了采用时间序列分析进行机床热误差建模方法，利用实测的热误差序列进行时序分析实时预测。杜柳青等^[10]提出混沌相空间重构理论的数控机床运动精度演化分析方法。为解决时间序列的相似性匹配问题，Agrawal等^[11]提出对时间序列在频域上进行处理的方法，达到了降低数据维数的目的，并且显著提高了相似匹配的效率；Lim等^[12]研究了相似性匹配中的多指标问题，并且详细讨论了滑动时间窗口维数与滑动步长对匹配精度的影响关系。虽然时间序列分析法广泛应用于状态预测，但大多是时间序列单维匹配，且没有考虑噪音数据，影响预测准确度。

针对以上算法中存在的问题，作者提出了一种基于多维时间序列的数控机床状态预测方法。利用多维时间序列表征机床未来状态，提出基于窗口滑动的多重匹配预测方法，为消除噪音数据的影响，提出了密度空间聚类算法的状态序列分析法。通过与现有算法进行实验对比，验证了本文方法的优势，为实现数控机床状态的准确预测提供了有价值的参考。

1 多维时间序列数据的分析处理方法 1.1 基于OPC技术的机床状态数据采集

OPC协议提供了3个主要的COM对象：OPC Server、OPC Group、OPC Items。OPC Server动态地创建或释放OPC Group，OPC Group提供包容OPC Items的机制，实现对OPC项的管理，OPC Items代表与OPC服务器到数据源（主轴电流、主轴负载、回转误差等）的连接，包括值、品质、时间戳3个基本属性。

SIEMENS 840D由人机界面，数控实时操作系统（NCK），数控单元和可编程逻辑控制器（PLC）组成。OPC数据采集客户端在人机界面上运行，使用组件对象模型（COM）访问OPC服务器，OPC服务器从NCU中的动态数据交换服务NCDDE请求相应的数据。结合多点接口（MPI）的特性，NCDDE服务器与PLC之间可进行数据交换，实现OPC客户端对数控机床数据的采集。

1.2 时间序列数据预处理 1.2.1 归一化处理

为消除时间序列相同量纲下大数值的影响，需要对各个参数进行归一化处理，针对机床参数数据的随机性和不稳定性，本文采用Min-max标准化法。

1.2.2 自回归移动平均模型

自回归移动平均模型（auto-regressive and moving average model，ARMA）是研究时间序列的重要方法，模型认为对于任意时间序列 ${X} = \left\{ {{x_1},{x_2},{x_3}, \cdots ,{x_t}} \right\}$ 在时刻 $t$ 的取值 ${x_t}$ 不仅与前 $n$ 个时刻的参数值 ${x_t}_{ - 1},{x_t}_{ - 2}, \cdots ,$ ${x_t}_{ - n}$ 有关，而且与之前 $m$ 个随机白噪声干扰项 ${a_t}_{ - 1},{a_t}_{ - 2}, \cdots ,$ ${a_t}_{ - m}$ 有关，则 ${x_t}$ 可表达为( $n$ ， $m$ )阶的ARMA模型^[13]：

$\begin{aligned}[b]{x_t} = & {\phi _1}{x_{t - 1}} + {\phi _2}{x_{t - 2}} + \cdots + {\phi _n}{x_{t - n}} + {a_t} - \\ & {\theta _1}{a_{t - 1}} - {\theta _2}{a_{t - 2}} -\cdots - {\theta _m}{a_{t - m}}\end{aligned}$

(1)

式中，误差干扰项 $a_t \sim NID(0,\sigma _a^2)$ ， $\phi_i(0 \le i \le n)$ 称为自回归系数， ${\theta _i}(0 \le i \le m)$ 表示移动平均系数。当 ${\theta _i} = 0$ 时，称为自回归模型记为AR（ $n$ ）（auto-regressive）；当 ${\phi _i} = 0$ 时，称为移动平均模型记为MA （ $m$ ）（moving average）。

考虑机床参数的实际特征，通过对序列 ${X}$ 建立自回归方程（AR）。若序列 ${X} = \left\{ {{x_1},{x_2}, \cdots ,{x_n}, \cdots {x_{2n}}} \right\}$ ，则序列元素满足方程：

$\left\{ \begin{aligned}& {x_{n + 1}} = {\phi _1}{x_n} + {\phi _2}{x_{n - 1}} +\cdots{\rm{ + }}{\phi _{n - 1}}{x_2}{\rm{ + }}{\phi _n}x{}_1{\rm{ + }}{{{a}}_{\rm{1}}},\\& {x_{n + 2}} = {\phi _1}{x_{n + 1}} + {\phi _2}{x_n} +\cdots{\rm{ + }}{\phi _{n - 1}}{x_3}{\rm{ + }}{\phi _n}x{}_2{\rm{ + }}{{{a}}_2},\\& \;\;\quad\quad\quad\quad\quad\quad\quad \cdots\;\cdots\\& {x_{2n}} = {\phi _1}{x_{2n - 1}} + {\phi _2}{x_{2n - 2}} + \cdots{\rm{ + }}{\phi _{n - 1}}{x_{n + 1}}{\rm{ + }}{\phi _n}x{}_n{\rm{ + }}{{{a}}_n}\end{aligned} \right.$

(2)

1.3 多维时间序列状态模型及度量模型

多维时间序列是指观测对象的一组指标在不同时期上获得的一系列观测值，并按时间顺序排列而构成的数字集合，如图1所示。若基于OPC状态数据采集获得的参数种类有 $m$ 个，且以采样频率100 ms获得 $n$ 个数值点，则机床过程数据形成 $m$ 行 $n$ 列的多维时间序列矩阵 ${{T}_{m \times n}}$ 。

参数集合根据相应分类规则分为单维时间序列向量 ${{{x}_1}}$ 和多维时间序列矩阵 ${X} = \left\{ {{x_1},{x_2},{x_3}, \cdots , {x_m}} \right\}$ ，并随着时间 $t$ 的变化不断延伸。预测机床未来的状态需要利用机床历史的状态数据，应建立数控机床的状态模型和研究 ${t_1}$ 时刻之后时间序列的预测方法。

在机床参数时间序列矩阵 ${{T}_{{{m}} \times n}}$ 中，令采样时间跨度为 $N$ 的序列矩阵集合 ${{T}_{{{m}} \times n}}$ 为机床的状态模型，表征机床设备状态在时段 $N$ 内的变化规律。图1中时间窗口 $N$ 内的矩阵集合 ${{T}} = \left\{ {{{T}_1},{{T}_2},{{T}_3}, \cdots ,{{T}_{n/N}}} \right\}$ 均为机床的历史状态，状态模型 ${{T}_{\left( {n/N} \right) + 1}}$ 为需要预测的机床未来状态。针对机床的状态模型集合 ${{T}_{{{m}} \times n}}$ ，建立统一的状态度量标准和评价函数是对模型进行定量处理的基础，采用特征向量、特征趋势距离标示状态模型，利用差异度对状态模型进行对比衡量。

对于任意的状态模型 ${T}_k \in {{T}_{m \times N}}$ ，特征向量 ${\mathit{\boldsymbol{F}}_k}$ 为（ ${f_{kx}}$ ， ${f_{ky}}$ ），分别表示机床状态矩阵中某时刻点向量与平均值向量的最大、最小偏差值。

$\begin{aligned}[b]& {{F}_k}({{T}_{m \times N}}) = ({f_{kx}},{f_{ky}}) =\\ & \left({{Max}_j}^N = \left| {\sqrt {\sum\limits_{i = 1}^m {{{({T_i}({t_j}))}^2}} } - \frac{{\sum\limits_{j = 1}^N {\sqrt {\sum\limits_{i = 1}^m {({T_i}({t_j}))^2} } } }}{N}} \right|, \right. \\ & \left. Mi{n_j}^N{\rm{ = }}\left| {\sqrt {\sum\limits_{i = 1}^m {{{({T_i}({t_j}))}^2}} } - \frac{{\sum\limits_{j = 1}^N {\sqrt {\sum\limits_{i = 1}^m {({T_i}({t_j}))^2} } } }}{N}} \right|\right)\end{aligned}$

(3)

式中， ${T_i}\left( {{t_j}} \right)$ 表示参数 ${m_i}$ 在一个状态段N内时刻点 ${t_j}$ 的参数值。

对于任意的状态模型 ${{T}_{{k}}} \in {{T}_{{{m}} \times n}}$ ，二范数 ${D_k}$ 为模型的特征趋势距离，衡量机床状态模型的变化规律和平稳性特征，其表达式为：

$D_k({{T}_{m \times N}}) = {\left\| {{T_i}({t_j})} \right\|_2} = \sqrt[2]{{\sum\limits_{j = 1}^N {\sum\limits_{i = 1}^m {{{({T_i}({t_j}))}^2}} } }}$

(4)

对于任意的状态模型 ${{T}_{k}} \in {T}_{{m} \times n}$ ，特征向量和特征趋势距离组成了二元组 ${{G}_k} = \left( {{F_k},{D_k}} \right)$ ，则任意两个状态模型 ${{T}_i}$ 和 ${{T}_j}$ 的差异度 ${\delta _{ij}}$ 为：

$\begin{aligned}[b]& {\delta _{ij}} = \left| {{G_i} - {G_j}} \right| = \sqrt[2]{{{{({F_i} - {F_j})}^2} + {{({D_i} - {D_j})}^2}}} =\\ & \sqrt[2]{{{{({{({f_{kx}})}_i} - {{({f_{kx}})}_j})}^2} + {{({{({f_{ky}})}_i} - {{({f_{ky}})}_j})}^2} + {{({D_i} - {D_j})}^2}}}\end{aligned}$

(5)

状态模型的差异度可衡量机床两个状态之间的综合差别，差异度越大则说明状态矩阵吻合度越低，相似性最小，反之相似性最大。

2 基于多重匹配的状态预测模型研究 2.1 多重匹配预测建模理论

图1中，时间序列矩阵 ${\mathit{\boldsymbol{T}}_{n/N + 1}}$ 为需要预测的机床未来状态，将机床状态对应参数的序列值进行独立分析，滑动时间窗口建模过程如图2所示。

图2中，时间窗口长度为 $l$ 的矩阵 ${({{X}_{n + 1}})_{m \times l}}$ 为需要预测的机床状态矩阵， ${({{X}_n})_{m \times l}}$ 为机床当前状态矩阵。从矩阵 ${{X}_n}$ 的当前采样时刻元素 ${x_n}$ 开始，以滑动单位 $w$ ，时间窗口长度 $l$ 进行逆向移动，可得第 $k$ 个时间窗口代表的状态矩阵通项为：

$\begin{aligned}[b]{{X}_k} = & [x(n - wk - l + 1),x(n - wk - l + 2),\\ & x(n - wk - l + 3),\cdots,x(n - wk)]\end{aligned}$

(6)

式中， $k$ =1， $2,\cdots,$ 代表第 $k$ 个历史状态， $w$ =1， $2,\cdots,$ 为滑动时间窗口的滑动间隔长度。一般地，时间窗口长度 $l$ 大于滑动时长 $w$ ，且矩阵中的任意元素 ${x_i}$ （ $0 \le i \le l$ ）均为 $m \times 1$ 的列向量。结合滑动单位为 $w$ 的时间窗口，离散时间序列被分割成了元素长度均为 $l$ 的独立历史状态集合 $\tilde {X} = \{ {X}_1,{X}_2,\cdots,{X}_k,\cdots\} $ 。

基于窗口滑动的多重匹配预测是在历史状态集合 $\tilde {X}$ 中寻找与当前状态矩阵 ${\mathit{\boldsymbol{X}}_n}$ 相关性最大，即寻找相似度最高的某一历史状态 ${\mathit{\boldsymbol{X}}_k}$ ，且状态 ${\mathit{\boldsymbol{X}}_k}$ 的跟随状态 $\mathit{\boldsymbol{X}}_k^*$ 具有与 ${\mathit{\boldsymbol{X}}_{n + 1}}$ 相同的变化趋势和规律，视为 ${\mathit{\boldsymbol{X}}_{n + 1}}$ 的预测值，如图3所示。矩阵 ${\mathit{\boldsymbol{X}}_k}$ 和 $\mathit{\boldsymbol{X}}_k^*$ ，以及矩阵 ${\mathit{\boldsymbol{X}}_n}$ 和 ${\mathit{\boldsymbol{X}}_{n + 1}}$ 分别构成图3中时间跨度为N的一个完整机床状态。多重匹配即在寻找历史状态 ${\mathit{\boldsymbol{X}}_n}$ 的过程中，利用相似性分析和聚类分析分别完成初始寻找和精确匹配。

2.2 状态序列相似性分析

机床状态的相似性分析是指在机床海量历史状态数据中筛选与给定序列具有相似行为规律（变化趋势、数值大小）的状态矩阵集合，主要包括序列分割和相似性度量。滑动时间窗口序列分割方法有效避免了参数在时间轴的扭曲，保证了历史状态向量具有相同的长度。结合机床参数的实际特征，提出基于矩阵相关性的 $\beta $ -耦合度相似度量标准。

对于任意的机床状态矩阵 ${{X}_{m \times n}} \!\!=\!\! {[{x_1},{x_2},\cdots\!,{x_i},\cdots\!,{x_m}]^{\rm T}}$ 和 ${{Y}_{m \times n}} = {[{y_1},{y_2},\cdots,{y_i},\cdots,{y_m}]^{\rm T}}$ ，其 $\beta $ –耦合相似度定义为：

${\rm cos}\; \theta (X,Y) = \sum\limits_{i = 1}^m {\left| {\frac{{\sum\limits_{j = 1}^n {({x_{ij}} - \lambda ({x_i})) \times ({y_{ij}} - \lambda ({y_i}))} }}{{\sum\limits_{j = 1}^n {{{({x_{ij}} - \lambda ({x_i}))}^2} \times {{({y_{ij}} - \lambda ({y_i}))}^2}} }}} \right|} /m$

(7)

式中， $\lambda ({x_i}) = \displaystyle\frac{1}{n}\sum\limits_{j = 1}^n {x_{ij}} $ ， $\lambda ({y_i}) = \displaystyle\frac{1}{n}\sum\limits_{i = 1}^n {{y_{ij}}} $ ， $\lambda ({x_i})$ 表示状态矩阵 ${{X}_{m \times n}}$ 中第 $i$ 行元素的平均值， $m$ 为状态矩阵的行数即元素种类个数， $n$ 为状态矩阵的列数即每个参数的长度。

$\beta $ –耦合相似度可衡量机床两个状态矩阵中参数变化趋势的平均相似性和趋近性，应满足以下性质：

1）对称性且非负： $C{o_\beta }({X_1},{Y_1}) = C{o_\beta }({Y_1},{X_1}) \ge 0$ ；

2）三角不等式： $0 \le C{o_\beta }({X_1},{Y_1}) \le C{o_\beta }({X_1},{\textit{Z}_1}) + C{o_\beta }({\textit{Z}_1},{Y_1})$ 。

根据式（7）的相似性度量模型，对于给定相似性阈值 ${\varphi _{\rm d}} \in (0,1)$ ，当2个序列的相似度满足 $C{o_\beta }({X_1},{Y_1}) \ge {\varphi _{\rm d}}$ 时，称序列互为 $\beta $ –耦合相似序列。相似度 $C{o_\beta }$ 越大，说明两个状态模型中各行元素的平均相关性和相似度越大。相似度量模型不仅用于判断两个状态矩阵的相似性，还用于确定最优滑动单位 $\mathord{\buildrel{\hbox{$\scriptscriptstyle\frown$}} \over w} $ 和预测时长 $\mathord{\buildrel{\hbox{$\scriptscriptstyle\frown$}} \over l} $ ，最优参数应满足：

$\begin{aligned}[b]& Number\{ C{o_\beta }(T_{\rm history}^*(\mathord{\buildrel{\hbox{$\scriptscriptstyle\frown$}} \over w} ,\mathord{\buildrel{\hbox{$\scriptscriptstyle\frown$}} \over l} ),{{X}_n}) \ge {\varphi _{\rm d}}\} \ge \\ & \quad\quad\quad Number\{ C{o_\beta }(T_{\rm history}^*(w,l),{{X}_n}) \ge {\varphi _{\rm d}}\} \end{aligned}$

(8)

式中， $T_{\rm history}^*(\mathord{\buildrel{\hbox{$\scriptscriptstyle\frown$}} \over w} ,\mathord{\buildrel{\hbox{$\scriptscriptstyle\frown$}} \over l} )$ 代表以参数 $\mathord{\buildrel{\hbox{$\scriptscriptstyle\frown$}} \over w} $ 和 $\mathord{\buildrel{\hbox{$\scriptscriptstyle\frown$}} \over l} $ 构造的历史状态集合， $T_{\rm history}^*(w,l)$ 为任意参数下的历史状态集合， $Number$ 为历史集合中满足相似性要求的状态矩阵个数。

2.3 状态序列聚类分析 2.3.1 聚类度量函数

度量函数是聚类分析的核心，是分析评价聚类元素共同点和差异度的关键因素。对于任意的两个状态序列 ${{X}_{m \times n}} = {[{x_1}\!,\cdots\!,{x_i}\!,\cdots\!,{x_m}]^{\rm T}}$ 和 ${{Y}_{m \times n}} = \left[ {{y_1}\!, \cdots\!,{y_i} \!,\cdots\!,} \right.$ $ {\left. {{y_m}} \right]^{\rm{T}}}$ 其一般形式为：

$d({X},{Y}) = \sqrt[\kappa ]{{{{\left| {{Y} - w \times {X}} \right|}^\kappa }}} = \sqrt[\kappa ]{{\sum\limits_{i = 1}^m {{{\sum\limits_{j = 1}^n {\left| {{y_{ij}} - w{x_{ij}}} \right|^\kappa} } }} }}(\kappa \ge 1)$

(9)

式中： $w$ 为关联权重；当 $\kappa $ 取值为1、2、∞时，度量函数 $d({X},{Y})$ 分别为曼哈顿距离、欧几里得距离、切比雪夫距离。

2.3.2 密度空间聚类算法

密度空间聚类算法（density-based spatial clustering of applications with noise，DBSCAN）是一种基于高密度联通区域的聚类算法，将类簇定义为高密度相连点的最大结合，可有效排除噪声数据的干扰，且能发现任意形状的类簇，如图4所示。图4中Q₁、Q₂、Q₃为核心对象，p₁、p₂、p₃、p₄、p₅、p₆为直接密度可达，q₁、q₂、q₃、q₄、q₅、q₆为密度可达，则称对象q₁到对象p₂、对象q₂到对象p₃、对象q₃到对象p₅密度相连。

密度空间聚类算法实际上是将数据样本聚类为同簇样本之间均为密度相连点，而不同簇之间不相连，不在所有簇中的样本点作为噪声数据而丢弃。为保证各维度指标在相同的刻度级别，本文选用欧几里得距离作为聚类度量指标。

密度空间聚类算法的难点主要在于核心半径 $r$ 、最小数目 $MinPts$ 的选取，为提高算法鲁棒性，通过下式确定参数取值：

$\begin{aligned}& r = w \times ({d_{\max }} - {d_{\min }}),\\ & MinPts \in (0,Num/3)\end{aligned}$

(10)

式中： ${d_{\max }}$ 和 ${d_{\min }}$ 分别为聚类元素状态中欧几里得距离的最大值和最小值； $w$ 为缩小权重参数，取值范围为（0，1）； $Num$ 为聚类元素的总个数。

若机床历史状态矩阵聚类中心簇个数为 $k$ ，则聚类簇中心点 ${C_1},{C_2}, \cdots\!,\;{C_k}$ 应满足：

$D({C_k}) = Min\sum\limits_{p = 1}^{M_k} {d({C_k}} ,{{U}_p})$

(11)

式中， $M_k$ 为第 $k$ 类中所有历史矩阵状态的样本数， ${\mathit{\boldsymbol{U}}_p}$ 为历史状态矩阵。

针对聚类元素为状态矩阵 ${s}$ 时，密度空间聚类算法首先计算集合 ${S} = \left\{ {{{s}_1},{{s}_2},{{s}_3}, \cdots ,{{s}_k}, \cdots ,{{s}_f}} \right\}$ 中任意两个状态矩阵的欧几里得距离 ${d_{i,j}}$ （ ${{{s}_i}}$ ， ${{{s}_j}}$ ， ${i \ne j}$ ），根据式（10）确定核心半径 $r$ ；然后，根据初始缩小权重 $w$ ，最小半径数目 $MinPts$ ，任意选择一个未处理状态矩阵 ${{{s}_k}}$ ，寻找所有与其密度相连的状态矩阵 $\mathit{\boldsymbol{S}}^*= \left\{ {{s}_1^*,{s}_2^*,{s}_3^*, \cdots ,} \right.$ $\left. {{s}_k^*, \cdots\!,\; {s}_f^*} \right\}$ ，再从 $\mathit{\boldsymbol{S}}^*$ 中任选一个状态矩阵 ${{s}_k^*}$ 继续寻找与其密度相连状态矩阵，直到集合 $\mathit{\boldsymbol{S}}$ 中的所有状态矩阵被遍历，利用式（11）得到聚类簇 $\mathit{\boldsymbol{C}} = \left\{ {{C_1},{C_2}, \cdots\!,\;{C_k}} \right\}$ ；最后，利用式（5）分别计算当前状态 ${\mathit{\boldsymbol{X}}_{{n}}}$ 与聚类中心 ${C_1},{C_2}, \cdots\!,\;{C_k}$ 的差异度 $\delta _{ij}$ 值可完成 ${\mathit{\boldsymbol{X}}_n}$ 的精度匹配。

3 数控机床状态的预测实验

本次实验的环境为宁江机床厂的SIEMENS 840D系列机床和“S”试件，CPU为Inter Core i5-6200U 2.30 GHz，4 G内存，测试软件为MATLAB 2014a。为验证本文算法的有效性，分别进行了基于多重匹配的向量形式、矩阵形式和AR模型的预测对比实验。

3.1 时间序列数据预处理

利用OPC数据采集客户端采集宁江机床SIEMENS 840D的运行数据，且运行数据在时间上具有一定的重复性。采集的主要参数有：主轴电机温度、主轴电机负载、主轴电流、主轴回转误差4项参数，采样周期为100 ms，机床切削S试件的累计有效加工时长（去除主轴转速为0时刻的数据）为12.8 h，产生的数据量总数为：4 $ \times $ 10 $ \times $ 60 $ \times $ 60 $ \times $ 12.8=1 843 200。

图5是4项参数在120 800～140 800 ms（20 s）内经过归一化处理后的图像，由图5可知，机床的回转误差、电机负载和主轴电流具有无序性和无规律性，而主轴电机温度取值比较均匀，主要在几个参数值中变化。

3.2 状态序列相似性分析

采用矩阵形式、向量形式对主轴4项参数进行相似性分析，取机床状态 $\left[ {{{T}},{{ T}} + l} \right]$ 为机床当前状态 ${\mathit{\boldsymbol{X}}_n}$ ，T代表采样时间，并确定最优滑动单位 $\mathord{\buildrel{\hbox{$\scriptscriptstyle\frown$}} \over w} $ 和预测时长 $\mathord{\buildrel{\hbox{$\scriptscriptstyle\frown$}} \over l} $ 。令 $w$ 以1个采样时间为公差在1～100个采样时间内依次递增，同时 $l$ 以10个采样时间为公差在10～50个采样时间内依次递增。在给定相似度阈值区间和滑动窗口范围内，按式（7）和（8）求出最大相似集，其结果分别为图6和7（以T=362 500，主轴回转误差为例）。

图6中，对机床主轴的4项参数以矩阵形式回滑1 000个时间窗口，相似度阈值取0.2、0.3和0.5。由图6可以看出，当预测时长为1 s时，图像峰值最高说明相似集合中元素最多，结合实际情况，预测1 s的数据意义不大。因此，最终的最大相似集合选用预测时长为24 s，滑动单位为2。同理，由图7可知，若主轴回转误差以单维向量的形式进行相似性分析，则应采用预测时长为24 s，滑动时长为2.3 s。

针对机床主轴的海量数据，以单维向量的形式将耗费大量计算时间，因此，对4种参数的状态预测采用矩阵形式，采用图6确定的最优预测时长和滑动单位，并回滑1 000个时间窗口，最大相似集合的求解结果如表1所示。

3.3 状态序列聚类分析

根据状态序列相似性分析得到在1 000个历史滑动窗口中满足 $\beta $ –相似耦合度 ${\varphi _{\rm d}} \ge 0.2$ 的481个历史状态矩阵 ${{\mathit{\boldsymbol{T}}}_{4 \times 240}}$ 进行密度空间聚类分析，利用式（10）确定参数为 $w$ =0.3， $MinPts$ =6。由于聚类元素为4行240列的矩阵元素难以在平面图中展示，因此利用状态特征值 $D\left( {{\mathit{\boldsymbol{T}}_K}} \right)$ 和特征向量模 $\left\| {\mathit{\boldsymbol{F}}\left( {{\mathit{\boldsymbol{T}}_K}} \right)} \right\|$ 分别作为横纵坐标画出其聚类结果如图8（a）所示。对聚类后的簇中心和当前状态 ${\mathit{\boldsymbol{X}}_n}$ 进行差异度分析，其结果如图8（b）所示。

图8（a）中，聚类矩阵元素的不同类簇之间沿着45°角排布，同一类簇元素大多沿着特征距离的 $x$ 轴方向排布，这是由于聚类的指标为欧几里得距离近似于矩阵的特征距离。同时由于部分聚类元素可能具有相同的特征向量和特征距离，因此聚点结果存在叠合的现象。由图8（b）可知，4个聚类簇的中心元素，其中聚类簇1的中心状态矩阵和当前机床状态矩阵的差异度最小，因此，簇1中心点代表的状态 ${C_1}$ 为多重匹配的结果。通过计算发现聚类簇 ${C_1}$ 是滑动过程中的第769个窗口，利用 ${C_1}$ 之后的240个采样点作为预测值。

3.4 机床状态预测结果

在完成时间序列数据预处理、状态序列相似性分析、状态序列聚类分析之后，对数控机床主轴4种参数使用基于多重匹配的矩阵形式、向量形式和传统的AR模型机床状态预测，其参数预测归一化结果如图9所示。为了便于区分，图中曲线分别做了不同程度的平移。图9（a） 中矩阵预测与向量预测分别向上平移了0.1与0.25个单位，图9（b）中向量形式预测向上平移了0.05个单位，图9（c）中向量形式预测向下平移了0.001个单位，图9（d）中矩阵预测与向量预测分别向上平移了0.15与0.3个单位。

对预测结果的评价，除差异度分析外，还可通过曲线的拟合量进行分析，采用最大误差 ${E_{\max }}$ 、平均误差 ${E_a}$ 、均方根误差MSE、相对误差率 $ee$ 进行评价。矩阵形式、向量形式以及AR模型的预测结果及分析如表2所示。

从表2的预测结果可以看出，基于多重匹配理论的矩阵形式、向量形式状态预测方法的各项误差均比传统的AR模型低，特别是基于多重匹配的向量形式状态预测方法可以保证较高的预测置信度，其相对误差最大值仅为10%，某些参数的相对误差低于1%，且差异度较小。相比于上述两种预测方法，基于传统的自回归移动平均模型精度较低，最大误差、平均误差、均方误差以及相对误差均较大，主轴回转误差相对误差最大值超过90%，失去了参数的真实变化规律，用于机床状态预测效果较差。因此，针对数控机床的多维时间序列，本文提出基于多重匹配的状态预测方法具有较高的预测置信度，可用于数控机床状态预测。

4 结　论

1）研究了多维时间序列数据的分析处理方法。利用OPC技术进行机床状态数据采集，对状态数据进行预处理，建立多维时间序列状态模型及其度量模型。

2）建立了基于多重匹配的机床状态预测模型。利用多重匹配预测建模理论，结合 $\beta $ –耦合相似度和空间密度算法开展了状态序列相似性和聚类分析。

3）通过实验验证了方法的可行性。实验结果表明，当预测时长为24 s，滑动单位为2 s时的相似度最高，预测实验结果表明，基于多重匹配的机床状态预测方法比传统预测方法置信度更高，预测误差更低。

4）研究中忽略了多维时间序列的相似性分析时计算量的问题，下一步将针对这一情况对算法作出改进。