作为混沌理论的重要课题，混沌相空间重构技术被广泛应用于多个领域^[1-5].当重构嵌入维数大于系统自由度时，重构轨迹与原系统高维流行空间微分同胚，可以恢复复杂系统的动力学特征，因此相空间重构技术有着十分重要的应用价值^[6-8].但是该技术存在多变量时滞参数的计算与选择这一重要问题.选择时滞过小，重构空间严重压缩扭曲; 选择过大, 则各变量坐标独立，很难提供系统吸引子的稳定投影^[9].所以对多变量最优时滞计算方法的研究具有重要理论意义.

在传统方法中，互信息法在多变量问题的应用上受到限制^{[1, 10]}.陶慧等^[2]使用遗传神经网络计算煤矿冲击地压的最佳重构参数; 张文专等^[3]将基于差分进化支持向量机与混沌理论相结合，根据模型预测效果进行时滞选择; 但受样本数量影响，易出现过拟合问题.张春涛等^[4]和张淑清等^[5]从信息熵的角度出发进行混沌时滞计算，利用了混沌系统的耗散性，物理意义明显，可解释性强，但计算复杂、误差较大，不易实现. Zhang等^[11]将最小微分熵原则与神经网络模型组合，选择结果更为精准，但是计算复杂和过于依赖模型的根本问题没有得到解决.上述方法难以推广用于多变量时滞的多参数同时选择. Palit等^[9、12-13]针对互信息法的局限性提出了自相关法和改进高维互信息法，但高维互信息只能近似估计，误差较大; 为此进一步研究了高维自相关算法，但可视化的定性评价方法仅适用于低维空间.

本文提出基于最大独立互相关算法进行多变量时滞的同步选择，并推广到井下瓦斯多变量时滞计算，解决多变量瓦斯预测参数选择困难的问题^[14].以联合递归图^[15-17]和共有近邻比值法^[18]评价结果的优劣.通过实验将本文算法与上述算法比较，并利用井下多变量瓦斯数据对算法进行真实数据验证.

1 最大独立多变量互相关算法

相空间重构实质上是对原始系统运动轨迹的还原过程，而响应变量与驱动变量的重构曲线是扭曲的系统轨迹.当时滞参数最优时, 相空间得到完美重构，此时的曲线扭曲程度最小，变量间的相关性也最小.根据Palit等^[13]的研究，相关系数仅能度量二元序列的线性关系，难以解决实际问题多变量非线性问题.因此，本文将相关系数改进为

$ \begin{array}{l} {R_{X, Y}}\left( \tau \right) = {\rm{ }}\\ \;\;\;\;\frac{{\sum\limits_{t = 1}^{N - \left( {m - 1} \right)\tau } {\left\{ {Y\left( t \right) - \overline {Y\left( t \right)} } \right\}} \left\{ {f\left( {X\left( t \right)} \right) - \overline {f\left( {X\left( t \right)} \right)} } \right\}}}{{\sqrt {\sum\limits_{t = 1}^N {{{\left\{ {Y\left( t \right) - \overline {Y\left( t \right)} } \right\}}^2}} } \sqrt {\sum\limits_{t = 1}^{N - \left( {m - 1} \right)\tau } {{{\left\{ {f\left( {X\left( t \right)} \right) - \overline {f\left( {X\left( t \right)} \right)} } \right\}}^2}} } }}. \end{array} $

式中：R()为互相关，m为嵌入维，τ为时滞参数，$X = \{ x\left( t \right)\} _{t = 1}^N$和$Y = \{ y\left( t \right)\} _{t = 1}^N$为两序列，N为序列长度，f()是从x(t)中对y(t)进行估计的函数，$\overline {Y\left( t \right)} $和$\overline {fX\left( t \right)} $分别为Y(t)和f(X(t))对应的平均值.

为将互相关系数用于解决多变量非线性问题，改进如下:

$ \begin{array}{l} R({\tau _1}, \ldots, {\tau _p}) = {\rm{ }}\\ \;\;\;\;\;\;\frac{{\sum\limits_{p, t = 1}^{N - \left( {m - 1} \right)\tau } {\left\{ {Q\left( t \right) - \overline {Q\left( t \right)} } \right\}\left\{ {W - \overline W } \right\}} }}{{\sqrt {\sum\limits_{p, t = 1}^N {{{\left\{ {Q\left( t \right) - \overline {Q\left( t \right)} } \right\}}^2}} } \sqrt {\sum\limits_{p, t = 1}^{N - \left( {m - 1} \right){\tau _{{\rm{max}}}}} {{{\left\{ {W - \overline W } \right\}}^2}} } }}. \end{array} $

式中：p为变量数，Q为含有p个变量的多元序列，W为对应的估值序列.

使用最小二乘曲面拟合作为优化方法，将p个变量拟合成一个曲面，变量间的相关性通过局部线性距离估计.

2 延迟时间的选择评价

本文以递归图和共有近邻比值法对实验结果进行评价.递归图是观测变量的时间序列在不同时刻访问同一区域二元图像^[14]，以黑色表示递归状态白色表示非递归状态，图 1为Lorenz系统x变量的递归图.

递归图的初衷是将高维动力系统的动态演变转化为一个N阶方阵并进行可视化展示^[15], 但可视化方法局限于三维空间展示，无法用于高维系统.本文使用联合递归密度来度量高维空间的扭曲程度，联合递归密度^[16]如下式:

$ \begin{array}{l} JR_{i, j}^M({\text{e}^{x(1)}}, \ldots, {\text{e}^{x(p)}}) = \mathop \prod \limits_{k = 1}^p R_{i, {\rm{ }}j}^{{m_k}}({\text{e}^{x(k)}}), {\rm{ }}\\ \;\;\;\;\;\;M = \sum\limits_{k = 1}^p {{m_k}}, {\rm{ }}i, {\rm{ }}j = 1, \ldots, N. \end{array} $

式中:JR()为联合递归密度函数，M为所有变量嵌入维之和，e为自然对数函数的底数.

为了对延迟时间选择进行定量评价，引入共有近邻比值法^[17]，度量相空间重构前后各变量轨迹共有近邻点的变化，比值越大则说明重构后共有近邻点减少的越多，相空间的扭曲程度越小重构效果越好，利用这一规律进行最优时滞选择的评价.

综上所述，本文提出的以联合递归密度和共有近邻比为评价方法的最大独立多变量互相关算法总结如下:

1) 初始化，设定嵌入维数、延迟时间和最大迭代次数，系统观测序列进行相空间延时重构，得到重构结果Y_n^m₁和X_n^m₂;

2) 对Y_n^m₁进行分段处理，得到p+1段子序列; $\{ {Y_{1, n}}, \ldots, {Y_{p + 1, n}}\} _{n = 1}^{^{150 \sim 500}}$;

3) 对p+1段序列进行曲面拟合，得到函数近似估计;

4) 将重构后的系统驱动变量代入拟合函数，计算得到估计值W，进行迭代运算;

5) 迭代结束，最大独立互相关对应的坐标值即为最优时滞.

3 实验验证与分析

为了验证本文算法相对其他算法的优越性，利用标准Lorenz混沌系统的仿真数据将本文算法与文献[2]的遗传神经网络预测法、文献[4]的极大联合熵法和传统的互信息法3种方法的计算结果对比，并采用递归图法和共有近邻比值法对计算结果进行评判. Lorenz系统的微分方程为

$ \begin{array}{l} \frac{{{\rm{d}}x}}{{{\rm{d}}t}} = s\left( {y - x} \right), \\ \frac{{{\rm{d}}y}}{{{\rm{d}}t}} = x\left( {r - z} \right) - y, , s = 10, r = 28, b = 8/3.\\ \frac{{{\rm{d}}z}}{{{\rm{d}}t}} = xy - bz{\rm{ }}. \end{array} $

采用四阶Runge-Kutta法求解方程，设定相关参数分别与陶慧^[2]、张春涛^[4]一致，进行互信息和本文算法实验.陶慧设定步长为0.04，x₀=15.34, y₀=13.68, z₀=37.91, m₁=m₂=3，对x、y各采样6 000点去除前5 500点的暂态过程; 张春涛设定步长为0.02，x₀=8, y₀=9, z₀=2, m₁=m₂=m₃=2，对x、y各采样8 000点去除前5 000点的暂态过程.直接采用文献[2, 4]中的结果分别为τ₁=3, τ₂=4和τ₁=9, τ₂=9, τ₃=9，二变量和三变量的互信息法和本文算法结果如表 1与图 2~5(时滞为量纲一变量).

对于二变量问题，本文算法选择结果共有近邻比最大(1.08)，表明相空间重构后Lorenz系统的二元变量轨迹联合空间轨迹展开程度最大，所以共有近邻点数减少最多，并且联合递归密度最大，说明变量内部的相似度最高，能够更好地恢复原Lorenz系统的动态特性，即选择结果最优; 其次是互信息法分别确定的结果，再次是遗传神经网络法，文献[2]中实验结果证明的遗传神经网络法优于互信息法, 是因为该方法过于依靠神经网络的非线性拟合能力，并且实验只选取了Lorenz系统的500个数据，神经网络精度依赖于带有标签信息的大数据量，迭代次数多、计算量大、普适性差，而对于相空间重构的时滞参数选择问题很难获取带有标记的大数据量训练样本，人工标注成本较高，只能对少量样本使用.而本文算法对数据量没有具体要求, 且能够对多变量同时确定延迟时间，优于互信息法.对于三变量问题，本文算法的共有近邻点变化最大，三变量联合递归状态呈现最多，相空间重构效果最好，选择最优，其次是极大联合熵法，再次互信息法.极大联合熵法克服了只能进行二元选择的局限性，但是多变量熵值的计算复杂，很难推广到更多变量的情况，而本文算法避免了该问题.但是在二元和三元变量小数据量情况下，算法的优越性没有得到最大体现，所以实验结果中的评价参数与其他算法相差较小.

为了说明本文算法的实际应用价值，对某煤矿井下巷道的瓦斯检测数据进行相空间重构，通过真实数据对算法进行验证.将整个回采工作面视为一个局部动态复杂系统，设工作面瓦斯为该复杂系统的响应变量，进风巷、上隅角和回风巷为系统的驱动变量，选择2015年1月1日至1月2日303回采工作面的每分钟历史检测瓦斯浓度数据，受篇幅所限在这里不具体列出数据.设定重构嵌入维为3，对巷道中4个位置的瓦斯浓度数据进行归一化处理，如表 2所示.

该巷道煤岩瓦斯动态系统受进风瓦斯、上隅角瓦斯、采面瓦斯和回风巷瓦斯共同影响，互相关最小为2.27×10^-8，4个位置对应最优时滞的归一化结果为16-3-10-11.进风巷瓦斯延迟时间最大，其次为采面和回风巷瓦斯时滞，上隅角瓦斯的延迟时间最小且与其他3个变量相差较多.由于巷道上隅角区域不易通风，且温度湿度都较高，容易积聚瓦斯且与进风巷、采面和回风巷瓦斯的交互流通较少.而进风巷-采面-回风巷3个位置通风顺畅，对应的3个瓦斯变量间相互作用、相互影响，所以四变量间最大独立时，这3个位置瓦斯变量的延迟时间较为接近. 0.58的共有近邻比和0.34%的联合递归密度说明受4个变量作用的局部巷道煤岩瓦斯动态系统较为复杂，导致4个巷道位置对应的瓦斯变量重构后的高维空间轨迹依旧严重的压缩扭曲，从而造成四变量较少的联合递归状态.

4 结论

1) 本文针对相空间重构中多变量时滞选择计算的问题，提出了最大独立互相关算法, 并引入最大联合递归密度与最大共有近邻比的方法作为实验评价标准.

2) 通过仿真实验证明了本文算法在面对二变量与三变量的时滞计算问题时优于传统方法.对井下进风巷、上隅角、采面和回风巷的瓦斯多变量进行实验，结果优于传统方法，具有现实意义.

3) 未来将研究多变量非线性混沌系统相空间重构时滞与嵌入维参数的耦合计算.