新型多吸引子混沌系统及其应用

一种具有多吸引子的新型混沌系统：动态分析、电路实现与S盒设计

摘要

本文报道了一种具有三个非线性的新型三维混沌系统。该系统在不同参数下分别具有一个稳定平衡点、两个稳定平衡点和一个鞍结点、两个鞍焦点和一个鞍结点。此新型系统的一个显著特征是由不同初始值引起的多吸引子现象。随着参数的变化，系统表现出单稳态、双稳态、单周期性、双周期性、一个奇异吸引子以及两个共存的奇异吸引子。通过分析相应的平衡点并结合数值模拟方法，揭示了该系统的复杂动力学行为。此外，给出了用于实现该系统混沌吸引子的电子电路。利用该新型混沌系统设计了一个用于密码操作的S盒。并且，我们测试了所生成S盒的性能，并将其与现有的S盒研究进行了比较。

关键词 : 新型混沌系统；多吸引子；电子电路实现；S盒算法

引言

1963年著名的洛伦兹吸引子[1]的发现，掀起了混沌研究的热潮。在此后的几十年中，关于混沌控制、混沌化、同步和混沌应用的一系列有意义的成果不断涌现，人们对混沌的认识也发生了巨大变化。学者们开始更多地思考如何产生混沌，而非盲目地抑制混沌。三维自治常微分系统中混沌吸引子的产生尤其受到关注。众所周知，已发现大量具有混沌吸引子的典型系统，包括罗斯勒系统、陈系统、斯普罗特系统、吕系统等。[2–8]

随着对混沌的深入研究，科学家发现某些非线性动力系统不仅具有混沌吸引子，而且在一组固定参数值下还存在多个共存吸引子。这些共存吸引子可能是不动点、极限环、奇异吸引子等。吸引子的数量和类型通常与系统的参数和初始条件相关。每个吸引子都有其自身的吸引域，由导致长期行为趋于该吸引子的初始条件组成。多吸引子现象在许多生物系统和物理系统中均可观察到[9–11]。近年来，具有多吸引子的低维自治混沌系统激发了学者的研究热情。

李和斯普罗特通过数值分析在混沌系统中发现了多吸引子，并引入了偏移提升法和条件对称法来在微分系统中产生多吸引子[12–16]。肯涅等人分析了可由微分方程描述的简单混沌电路中的多吸引子[17,18]。包和徐提出了具有多个混沌吸引子的基于忆阻器的电路系统[19,20]。赖等人提出了一些具有多吸引子的三维和四维连续混沌系统[21–23]。魏等人试图通过分析系统的分岔来揭示多吸引子的内在机制[24]。混沌及多个共存吸引子的研究确实是学术界一个非常有趣的研究课题，有助于认识实际系统的动态演化，并推动复杂性科学的研究。

混沌系统已被发现应用于许多领域，其中最有价值的应用是密码学。由于混沌系统具有丰富的动态行为和初值敏感性，能够提供加密所必需的混淆和扩散特性[25,26]。S盒是分组加密算法中具有扰乱功能的最基本单元。一个良好的S盒可以使加密算法具有更高的安全性以及更强的抗攻击能力。尽管已有大量关于混沌S盒设计的研究，但根据一些独特的混沌系统生成S盒仍然具有重要意义。在将混沌系统应用于工程领域之前，有必要通过电子电路实现该系统，以证明其真实存在。基于电路理论和简单电路元件，可以在示波器上产生混沌信号。迄今为止，电子电路已成为分析混沌系统的重要工具[27–30]。

本文研究了一类具有以下特征的特殊多项式混沌系统：（i）系统包含三个非线性项xz、yz、xyz，且在变换(x,y, z) → (−x,−y,z)下具有不变性；（ii）系统可产生蝴蝶吸引子；（iii）在不同参数条件下，系统分别具有一个稳定平衡点、一个不稳定平衡点、两个稳定平衡点、三个不稳定平衡点，并表现出单稳态、双稳态、单周期性、双周期性、一个奇异吸引子以及两个共存的奇异吸引子。在分析该系统的动态行为之后，根据该系统设计了电子电路和S盒。

本文结构安排如下：第2节描述混沌系统并展示其蝴蝶吸引子。第3节分析平衡点的稳定性。第4节研究系统的动态行为。第5节考虑系统的电子电路实现。第6节基于该系统建立S盒，第7节总结本文结论。

2. 混沌系统的描述

本文提出的混沌系统可以用以下微分方程组表示：
$$
\begin{cases}
\dot{x} = ax - yz, \
\dot{y} = -by + xz, \
\dot{z} = -cz + xyz + k,
\end{cases}
\tag{1}
$$
状态变量为 $(x, y, z) \in \mathbb{R}^3$，参数向量为 $(a, b, c, k) \in \mathbb{R}^4$。通过在MATLAB软件（MATLAB 8.0，MathWorks，美国马萨诸塞州纳蒂克）上进行数值仿真，可以观察到蝴蝶吸引子。在参数 $(a, b, c, k) = (4, 9, 4, 4)$ 和初始条件 $(1, 1, 1)$ 下，系统(1)的相图如图1所示。该图直观地表明，系统(1)呈现出一个吸引子，因为系统轨迹最终会进入一个有界区域。系统的李雅普诺夫指数计算结果为 $l_1 = 1.7729$、$l_2 = 0.0000$、$l_3 = -7.5949$。李雅普诺夫维数为 $D_l = 2 - l_1 / l_3 = 2.2334$，因此可以确定该吸引子为混沌吸引子。由两个非常接近的初始条件 $(1, 1, 1)$ 和 $(1, 1, 1.001)$ 生成的变量 $z$ 的时间序列如图2所示。起初它们几乎相同，但经过若干次迭代后差异逐渐增大。这说明系统(1)具有对初始条件的敏感依赖性，其长期行为无法预测。系统(1)的庞加莱映射是通过选取截面获得的
$$
\Delta_1 = {(x, y) \in \mathbb{R}^2 \mid z = 10} \quad \text{和} \quad \Delta_2 = {(y, z) \in \mathbb{R}^2 \mid x = 0}.
$$
如图3所示，庞加莱映射是一组点集，这与混沌的特性一致。

3. 平衡点的稳定性

假设参数a、b、c、k均为正实数。通过求解 $\dot{x} = \dot{y} = \dot{z} = 0$，可得到系统(1)的平衡点。当 $k \geq c\sqrt{ab}$ 时，系统(1)只有一个平衡点 $O(0, 0, k/c)$。当 $k < c\sqrt{ab}$ 时，系统(1)具有三个平衡点，如下所示：
$$
O(0, 0, k/c), \quad O_1(\sqrt{(c\sqrt{ab}-k)/a}, \sqrt{(c\sqrt{ab}-k)/b}, \sqrt{ab}), \quad O_2(-\sqrt{(c\sqrt{ab}-k)/a}, -\sqrt{(c\sqrt{ab}-k)/b}, \sqrt{ab}).
$$

命题1 。假设 $b > a > 0$，$k > 0$，且参数 $c$ 满足以下条件：
$$
\frac{k}{\sqrt{ab}} < c < \frac{2k[(a^2 + b^2)\sqrt{ab} + k(b-a)]}{\sqrt{ab}[(a + b)^2\sqrt{ab} + k(b-a)]},
\tag{2}
$$
则系统(1)的平衡点 $O_1$ 和 $O_2$ 是渐近稳定的。

证明 。通过对系统(1)在平衡点处进行线性化，得到雅可比矩阵为
$$
H =
\begin{pmatrix}
a & -z & -y \
z & -b & x \
yz & xz & xy - c
\end{pmatrix}.
\tag{3}
$$
通过使用 $|\lambda I - H| = 0$，可在 $O_1, O_2$ 处得到相应的特征方程
$$
\lambda^3 + w_1\lambda^2 + w_2\lambda + w_3 = 0,
\tag{4}
$$
其中
$$
w_1 = (b - a + \frac{k}{\sqrt{ab}}),
\quad
w_2 = (a - b)(c - \frac{2k}{\sqrt{ab}}),
\quad
w_3 = 4ab(c - \frac{k}{\sqrt{ab}}).
$$
根据劳斯–赫尔维茨判据，当方程(4)的所有根的实部均为负值时，平衡点 $O_1, O_2$ 是稳定的。这要求 $w_1 > 0$，$w_2 > 0$，$w_3 > 0$ 以及 $w_1w_2 > w_3$。容易验证，当 $w_1 > 0$，$w_2 > 0$，$w_3 > 0$ 且 $b > a > 0$，$k > 0$ 时，
$$
c > \frac{k}{\sqrt{ab}}.
\tag{5}
$$
为了使 $w_1w_2 > w_3$，参数 $c$ 应满足
$$
c < c_0 = \frac{2k[(a^2 + b^2)\sqrt{ab} + k(b-a)]}{\sqrt{ab}[(a + b)^2\sqrt{ab} + k(b-a)]}.
\tag{6}
$$
由于 $c_0 < 2k / \sqrt{ab}$，则当 $b > a > 0$，$k > 0$, $\frac{k}{\sqrt{ab}} < c < c_0$ 时，$O_1, O_2$ 是渐近稳定的。当参数 $c$ 经过临界值 $c_0$ 时，发生双Hopf分岔，从 $O_1, O_2$ 分支出两个极限环，且系统(1)失去稳定性。

命题2 。假设 $b > a > 0$，$c > 0$，$k > 0$，则：（i）当 $c > k / \sqrt{ab}$ 时，平衡点 $O$ 不稳定；（ii）当 $c \leq k / \sqrt{ab}$ 时，平衡点 $O$ 渐近稳定。

证明 。在 $O$ 处计算的特征方程为
$$
(\lambda + c)[c^2\lambda^2 + (b - a)c^2\lambda + k^2 - abc^2] = 0.
\tag{7}
$$
如果 $c > k / \sqrt{ab}$，那么方程(7)有一个实部为正值的根。因此，$O$ 是不稳定的。如果 $c < k / \sqrt{ab}$，方程(7)的所有根都具有负实部，这意味着 $O$ 是稳定的。当 $c = k / \sqrt{ab}$ 时，方程(7)有三个根 $\lambda_1 = 0$, $\lambda_2 = a - b$, $\lambda_3 = -c$。因此，$O$ 是非双曲平衡点。可以通过应用中心流形定理验证 $O$ 是渐近稳定的。

4. 多吸引子的演化

本节详细研究了系统(1)的复杂动力学行为。通过包括分岔图、相图、李雅普诺夫指数和庞加莱映射在内的仿真实验，对系统(1)进行了深入且直观的观察。系统(1)中的吸引子具有分形性质，表现出丰富的混沌动力学特性。随着参数的变化，系统(1)经历了稳定状态、周期状态和混沌状态。对于不同的初值，系统(1)呈现出具有独立吸引域的不同类型的吸引子。

4.1. 参数c的动力学演化

考虑在给定参数条件 $a = 2$, $b = 8$, $k = 4$ 下，系统(1)随参数 $c$ 的动态演化。图4a 显示了系统(1)随 $c \in (0, 6)$ 变化的分岔图，其中红色分支和蓝色分支分别由初值 $x_{01} = (1, 1, 1)$ 和 $x_{02} = (-1, -1, 1)$ 得到。红色分支与蓝色分支的重叠区域表明 $x_{01}, x_{02}$ 的轨迹最终趋于同一吸引子，而分离区域则表明 $x_{01}, x_{02}$ 的轨迹趋于不同的吸引子。图4b 为初始值 $x_{01}$ 下系统(1)的李雅普诺夫指数。结果显示，随着 $c$ 的变化，系统(1)经历了稳定状态、周期状态和混沌状态。当 $c \in (0, 1)$ 时，系统(1)为单稳态，因其仅有一个稳定平衡点。当 $c \in (1, 1.396)$ 时，由于存在两个稳定平衡点，系统(1)呈现双稳态。当 $c$ 增大并越过临界值 $c_0 = 1.396$ 时，系统(1)在平衡点处发生双Hopf分岔。当 $c \in (1.396, 1.516)$ 时，系统(1)表现为双周期性。当 $c \in (1.516, 2.257)$ 时，系统(1)转变为单周期状态。当 $c \in (2.821, 3.096)$ 时，系统(1)从初值 $x_{01}, x_{01}$, 产生两个奇异吸引子。当 $c \in (3.310, 4.167) \cup (4.370, 5.324) \cup (5.480, 6)$ 时，系统(1)仅有一个混沌吸引子。表1 描述了系统(1)在不同c值下的吸引子。图5中的相图说明了系统(1)中存在不同类型的吸引子。

表1。不同c值下系统(1)的吸引子。

参数 $c$ 的值	平衡点	吸引子类型	图
$c = 0.8$	稳定点：$(0, 0, 5)$	一个点吸引子	图5a
$c = 1.2$	鞍结点: $(0, 0, 3.3333)$ 稳定点：$(\pm0.6325,\pm0.3162, 4)$	一对点吸引子	图5b
$c = 1.4$	鞍结点：$(0, 0, 2.8571)$ 鞍焦点: $(\pm0.8944,\pm0.4472, 4)$	一对极限环	图5c
$c = 1.6$	鞍结点：$(0, 0, 2.5)$ 鞍焦点：$(\pm1.0954,\pm0.5477, 4)$	对称极限环	图5d
$c = 2.9$	鞍结点：$(0, 0, 1.3793)$ 鞍焦点：$(\pm1.9494,\pm0.9748, 4)$	一对奇怪吸引子	图5e
$c = 3.1$	鞍结点: $(0, 0, 1.2903)$ 鞍焦点：$(\pm2.0494,\pm1.0247, 4)$	一对极限环	图5f
$c = 3.2$	鞍结点：$(0, 0, 1.25)$ 鞍焦点：$(\pm2.0976,\pm1.0488, 4)$	对称极限环	图5g
$c = 3.6$	鞍结点：$(0, 0, 1.1111)$ 鞍焦点：$(\pm2.2804,\pm1.1402, 4)$	蝴蝶型奇异吸引子	图5h

4.2. 参数k的动态演化

系统(1)在参数 $(a, b, c) = (4, 9, 4)$，$k \in (5, 25)$ 下的分岔图如图6a所示，其中红色分支和蓝色分支分别由初值 $x_{01}, x_{02}$ 得到。显然，当参数 $k$ 从5增加到25时，系统(1)的状态从混沌转变为周期，再变为稳定。这一点也可以通过图6b中的李雅普诺夫指数来说明。最大李雅普诺夫指数在 $k \in (5, 13.6) \cup (13.9, 14.8) \cup (15.4, 15.9)$ 时为正值，在 $k \in (19.8, 25)$ 时为负值，在 $k \in (13.7, 13.8) \cup (14.9, 15.3) \cup (16, 19.7)$ 时等于零。当 $k = 5$ 分别为15、18、25时，可以观察到系统(1)的奇异吸引子、极限环和稳定点，其相图如图7所示。当 $k = 19$ 为20时，可以观察到系统(1)的两个共存周期吸引子和两个共存点吸引子，如图8所示。

5. 电子电路实现

文献中有许多与基于混沌的应用相关的工作[31–36]。这里，我们将展示系统(1)的电路实现，以实际获得其混沌吸引子。图5e中的数值仿真显示了在参数 $(a, b, c, k) = (2, 8, 2.9, 4)$ 和初始条件 $(\pm1, \pm1,1)$ 下，系统(1)中存在两个共存的奇异吸引子。为了实现该状态下系统(1)的电路实现，我们需要避免电路元件饱和，而实现这一目标的有效方法是通过对系统(1)的变量进行缩放来降低电路电压值。在缩放过程中，假设 $(X,Y, Z) = (x, y, z/2)$，然后得到缩放后的系统为
$$
\begin{cases}
\dot{X} = aX - 2YZ, \
\dot{Y} = -bY + 2XZ, \
\dot{Z} = -cZ + XYZ + \frac{k}{2}.
\end{cases}
\tag{8}
$$
图9给出了缩放系统(8)在参数 $(a, b, c, k) = (2,8, 2.9,4)$ 下的新相图。显然，缩放过程并未对系统(1)造成根本性改变，而只是将变量限制在更小的区域 $\Omega = {(x, y, z) \mid x, y \in (-5, 5), z \in (0, 10)}$ 内。

由OrCAD-PSpice程序（OrCAD 16.6，OrCAD公司，美国俄勒冈州希尔斯伯勒）生成的系统(8)的电路图如图10所示。该电路具有对应于变量X、Y、Z的三个输入（或输出）信号，信号间的运算通过基本电子元件实现，包括电阻、电容、TL081运算放大器(运放)和AD633乘法器。通过固定 $R_1 = R_3 = 20K\Omega$、$R_2 = 200 K\Omega$、$R_4 = 50 K\Omega$、$R_5 = R_6 = 100 K\Omega$、$R_7 = 138K\Omega$、$R_8 = 3000 K\Omega$、$R_9 = 4 K\Omega$、$C_1 = C_2 = C_3 = 1 nF$、$V_n = -15 V$、$V_p = 15 V$，并在图11所示的电子电路板上运行该电路，可在示波器上获得电路输出。图12中的示波器图形与数值模拟结果具有良好的一致性。

6. S盒设计及其性能分析

本节旨在通过应用系统(1)提出一种新的混沌S盒算法。在算法设计中，首先进行随机数生成，然后生成S盒。S盒生成算法的伪代码如算法1所示。为了建立S盒生成算法，我们首先输入系统(1)的参数 $(a, b, c, k) = (2, 8, 2.9, 4)$ 和初始值 $(x_0, y_0, z_0) = (-1,-1, 1)$，然后生成浮点数输出。为了在混沌系统的分析中产生更加随机的输出，我们选择一个合适的步长间隔 $\Delta h$ 并将其作为采样值使用。通过设置合适的步长间隔 $\Delta h = 0.000001$，可以获得更多的随机序列。利用RK4算法结合初始条件和指定采样值求解系统(1)，得到时间序列。在我们设计的混沌S盒算法中，使用了系统(1)的y、z相位的输出。从这两个相位获得的浮点数值（32位）被转换为二进制系统。通过对两个相位生成的32位数序列的最低有效位(LSB)各取8位，对这些值进行异或运算。所得到的新8位值被转换为十进制数。如果该十进制数此前已生成并包含在S盒中，则丢弃该值；否则将其添加到S盒中。以此方式继续该过程，直到在S盒中获得256个不同的值（介于0到255之间）。生成的S盒如表2所示。

算法1 S盒生成算法的伪代码。

1: 开始;
2：输入系统的参数和初始值；
3：以步长间隔 $\Delta h$ 进行采样
4: $i = 1$，S盒 $=[]$；
5: 当 $(i < 257)$ 时执行
6: 使用RK4算法求解系统并获得时间序列 $(y，z)$；
7: 将浮点数转换为二进制数；
8: 从RNG $y \oplus z$ 相位中提取最低有效位‐8的值；
9: 将二进制数转换为十进制数（8位）
10: 如果(S‐盒 =中是否存在小数值 yes) 那么
11: 转到第6步。
12: 否则(S‐盒 =中是否存在小数值 no)
13: S盒[i] ← 十进制值；
14: $i++$;
15: end
16： 结束
17: S‐盒 ← 重塑(S盒,16,16);
18: 准备使用基于混沌的S盒 16 × 16；
19: 结束。

表2。系统(1)的混沌S盒。

199	30	5	41	38	140	230	139	66	0	11	195	76	204	54	23
254	198	50	108	231	92	87	182	217	28	56	253	219	232	215	49
102	151	68	86	176	248	12	32	126	249	141	154	82	138	174	165
145	62	115	150	201	104	170	148	78	97	192	247	252	96	211	153
45	98	40	91	109	113	196	107	209	83	144	120	191	75	242	208
175	246	100	181	85	70	197	136	235	210	93	216	71	105	162	149
88	240	31	238	42	171	90	73	112	243	255	128	239	121	26	34
25	226	59	244	135	142	53	36	146	157	117	124	116	10	205	60
173	29	2	72	203	3	214	224	127	241	143	74	6	156	122	61
110	8	1	233	79	51	77	47	236	222	185	152	180	15	103	234
206	227	169	202	137	221	177	179	163	52	245	67	89	80	220	7
237	183	17	4	101	37	39	57	178	194	58	69	213	147	18	228
46	35	225	84	14	125	95	134	129	63	99	55	106	161	218	27
250	21	13	24	207	193	48	184	189	114	111	167	16	160	188	123
155	132	158	130	118	166	164	168	33	159	223	64	44	81	190	172
212	20	229	186	65	251	133	22	131	43	119	94	19	9	187	200

为了确定生成的S盒能够有效抵抗攻击，对其进行了若干性能测试。我们主要关注以下几项测试：非线性、输出比特独立准则(BIC)、严格雪崩准则(SAC)以及差分逼近概率(DP)。此外，还将本文提出的新S盒与Chen[37], Khan[38], Wang[39], Ozkaynak[40], Jakimoski[41], Hussain[42], Tang[43]等人提出的现有混沌S盒的性能对比结果列于表3中。

非线性被认为是所有性能测试中最核心的部分。由系统(1)生成的S盒的非线性值分别为104、106、104、104、108、104、110和104。因此，其平均值、最小值和最大值分别计算为105、104和110。通过比较表3中其他S盒的非线性，我们可以说新的S盒在某种程度上优于其他S盒。

表3. 不同混沌 S盒 的比较（BIC：比特独立性准则；SAC：严格雪崩准则；DP：差分概率）。

S盒	非线性 Min	非线性 Avg	非线性 Max	BIC-SAC Min	BIC-SAC Avg	BIC-SAC Max	BIC 非线性 Min	BIC 非线性 Avg	BIC 非线性 Max	SAC Min	SAC Avg	SAC Max	DP Min	DP Max
提出的S盒	104	105	110	0.5028	102.75	—	0.3906	0.5014	0.5937	10	—	—	—	—
Chen [37]	100	103	106	0.5024	103.1	—	0.4218	0.5000	0.6093	14	—	—	—	—
Khan [38]	96	103	106	0.5010	100.3	—	0.3906	0.5039	0.6250	12	—	—	—	—
Wang [39]	102	104	106	0.5070	103.8	—	0.4850	0.5072	0.5150	12	—	—	—	—
Ozkaynak [40]	100	103.2	106	0.5009	103.7	—	0.4218	0.5048	0.5938	10	—	—	—	—
Jakimoski [41]	98	103.2	108	0.5031	104.2	—	0.3761	0.5058	0.5975	12	—	—	—	—
Hussain [42]	102	105.2	108	0.5053	104.2	—	0.4080	0.5050	0.5894	12	—	—	—	—
Tang [43]	99	103.4	106	0.4995	103.3	—	0.4140	0.4987	0.6015	10	—	—	—	—

SAC由Webster等人提出[44]。一般来说，SAC的建立意味着当单个比特发生变化时，每个输出比特中大约有一半会发生变化。表3给出了新S盒的SAC平均值、最小值和最大值，分别为0.5014、0.3906和0.5937。显然，新S盒的平均值接近理想值0.5。BIC也是Webster等人提出的一个重要准则[44]，可用于部分衡量密码系统的安全性。通过翻转明文的一位所生成的向量集，被用来测试其与所有雪崩变量对的独立性。在测量雪崩效应之间的关系时，需要变量对来计算相关值[45]。在计算BIC值的同时会计算BIC‐SAC和BIC‐非线性值。当查看表3中的数值时，BIC‐SAC值计算如下：平均值为0.5028，最小值为0.4394，最大值为0.5312。该平均值几乎等于最优值0.5。

DP是比汉等人建立的用于测试S盒的另一个性能指标[46]。在此分析中，确定了S盒输入和输出位之间异或分布的平衡性。输入和输出位之间异或分布概率越接近，通常表示S盒抵抗差分攻击的能力越强。较低的DP值表明S盒更具抗攻击能力。新S盒的最小和最大DP值分别为4.0和10。从表3可知，新S盒的DP值与唐和奥兹卡亚克提出的S盒相同。

在通过一些重要指标测试新S盒的性能并与其他S盒进行比较后，可以确定由系统(1)生成的新S盒比其他S盒具有更好的性能。因此，它将更适用于抗攻击和强加密。

7. 结论

本文研究了一种具有多吸引子的特殊混沌系统。该系统的复杂动力学行为主要通过数值模拟进行展示。分岔图和相图表明，随着系统参数的变化，该系统呈现出一对点吸引子、一对周期吸引子以及一对奇异吸引子。此外，设计了用于实现该系统混沌吸引子的电子电路。同时，利用该混沌系统生成了一个新的S盒，并对该S盒的性能进行了评估与比较。结果表明，新S盒的性能优于一些现有的S盒。

事实上，具有多吸引子的混沌系统研究是当前的热点问题。与此相关的更多重要课题将在我们未来的论文中加以探讨。