基于分形混沌系统的多图像加密算法

Multi-Image Encryption Algorithm Based on Fractal Chaotic System

摘要: 构造了一个迭代函数系统，生成充满整个单位正方形的填充曲线。由迭代函数系统生成一个单位正方形上的分形混沌系统，并从混沌系统的轨迹、李亚普诺夫指数和0~1测试等方面验证了该系统具有优良的混沌性能。基于构造的二维分形混沌系统，提出了密钥相关的Z字形变换，并结合自适应Arnold变换和有限域运算设计了一个针对多个灰度图像的混沌图像加密算法，详细分析了该加密算法的安全性。

Abstract: An iterated function system is constructed to generate the filling curve on the unit square. One fractal chaotic system is derived by the iterated function system on the unit square. The chaotic natures of the fractal chaotic system are proved by the orbits, the Lyapunov exponents and 0~1 test, etc. Based on the constructed two-dimensional chaotic system, a chaotic image encryption algorithm for multiple gray images is designed by key-determined zigzag transformation, Arnold transformation and finite field operations. The security and performance of the proposed image encryption algorithm is analyzed in detail.

1. 引言

随着科学技术的不断进步，多媒体的信息安全已经成为人们越来越关心的问题，其中图像信息作为一种具有显著信息的载体，在军事、医学、教育和商业等方面具有广泛应用。所以图像加密已经成为信息安全领域的研究热点。由于图像信息的数据量大，冗余度高和相邻像素相关性强等特点，因此传统的加密方法如AES、DES等不适合数字图像的加密。因为混沌系统具有对初始值的敏感性、遍历性和随机性的特性，1989年R. Matthews提出了用混沌序列生成伪随机数用于数据加密。1998年，Fridrich提出了一种基于混沌的图像加密体系结构，利用混沌序列生成伪随机数，主要框架为置乱和扩散的数字图像加密方法被广泛利用在数字图像加密领域中 [1] 。

常用于图像加密的一维混沌系统有Logistic映射、Sine映射和Tent映射等，它们结构简单，易于实现 [2] 。但是由于计算设备精度有限，截断和舍入误差会使混沌系统陷入周期态，所以密码学家在一维混沌映射的基础上，提出了许多复合型的混沌系统。文献 [3] 提出了改进的Logistic-Tent混沌系统，文献 [4] 提出了复合的Logistic-Sine混沌系统，比原映射有更好的随机性和遍历性。文献 [5] 提出一个新的正弦–幂函数混沌系统，映射具有更长的周期窗口，更好的混沌特性。此外，高维Liu系统和Lorenz系统等也常常被用作图像加密 [6] [7] 。高维的混沌系统比起低维的混沌系统具有更大的密钥空间，用于加密时更不容易被暴力破解。但常常运算量大，耗时长，难以复现 [8] 。二维的混沌系统相较传统的低维系统安全性能更好，运算耗时和实现难度又优于高维系统，是密码学家用来图像加密的高效选择。文献 [9] 构造了一个具有大混沌范围的二维Logistic-Gaussian超混沌映射；文献 [10] 提出了一种新的二维正弦混沌系统，实验结果表明该系统比传统的系统具有更大的参数空间。

分形几何学用数学公式描述自然界中混沌、无序的复杂事物，使得分形与混沌之间有着紧密的天然联系。1890年科学家Peano发现可以使用迭代函数系统(IFS)生成能够充满平面的曲线(填充曲线) [11] ，这样的迭代函数系统在单位正方形内存在唯一的吸引子，由于IFS的压缩性，可以诱导得到分形吸引子上具有混沌特性的动力系统。本文对 [12] 的分形混沌系统进行改进，提出了一个二维三参数的分形混沌系统，给出了分形的迭代函数系统，并生成单位正方形的填充曲线。本文时序图，李雅普诺夫指数图和0~1测试等方面分析了该分形混沌系统的混沌特性。所提出的新分形混沌系统与低维系统相比具有更好的动力学特性，并且拥有运算速度快的优点；生成的混沌序列具有更好的遍历性和初值敏感性，适合生成用于图像加密的混沌序列。

置乱操作是通过打乱数字图像矩阵像素点的位置，实现图像的加密。矩阵置乱方法是直接用相应大小的矩阵，对图像数字矩阵位置进行变换。因为像素位置是整数，所以要求变换的矩阵的值是整数。著名的Arnold矩阵是整数矩阵且逆矩阵仍然是整数矩阵，一般应用在混沌加密中，相较于其他置乱方法来说，其置乱效果较好 [13] 。文献 [14] 将Arnold矩阵推广到高维情况；文献 [15] 分析了图像尺寸和Arnold变换周期的关系，为了改进周期性，提出了一种利用高维的Arnold变换加密彩色图像的置乱算法；文献 [16] 提出了一种基于Arnold的加密算法，在变换过程中加入密钥用来加密彩色图像；文献 [17] 运用三维Arnold映射结合混沌序列的算法对图像进行加密，算法取得了较好的加密效果。本文在 [18] 自适应图像加密算法的基础上，将分块之后的图像经过Arnold变换，结合与混沌序列的运算，以扩散的形式混淆每块图像的像素值，在有限次的变换后实现良好的加密效果。基于异或加法模运算和Arnold变换的扰乱和图像的明文相关，具有自适应性。

Zigzag变换方法是用Z字形的扫描办法扫描想要加密的图像矩阵，然后按照顺序放入新矩阵，即加密后的矩阵，是一种像素级的置乱方法。因为Z形扫描方法容易实现，加密效果好，所以被广泛使用。文献 [19] 使用随机置乱和Zigzag形置乱方法加密图像；文献 [20] 提出了四个不同方向的Zigzag扫描算法；文献 [21] 基于神经网络算法，提出任意起点的Z字形扫描方式；文献 [22] 将传统算法改进，提出了三维的Zigzag扫描方法，周期性更大，安全性更好。本文提出了适用于任何大小矩阵的Zigzag型扫描，在加密多个灰度图像的算法中，使用明文图像Hash函数生成的密钥决定了每幅图像Zigzag型变换的起点和终点，使算法具有更好的明文敏感性和随机性。

扩散算法是在不改变像素位置的基础上，改变像素值的大小，使用有限域实现图像加密算法的扩散，是安全又高效的选择。基于有限域上的多项式运算，其运算规律和结构，生成的有限域上的四则运算表在模运算的结果下，个数有限，大小有限制，所以其在计算机科学与信息理论方面都有较为广泛的应用 [23] 。在有限域理论中GF(257)的算术运算以257为模，此集合共包含257个元素。在乘法运算时，剔除元素0以减少因其造成的信息损失，这时集合中元素的个数为256，正好与图像像素点灰度值的取值范围对应 [24] 。文献 [25] 构造了GF(2⁸)的加法、减法表，扩散加密图像；文献 [26] 构造了GF(257)的乘法表，与混沌序列结合实现图像的扩散；本文构造了GF(257)的乘法、除法表和加法、减法表，在查表法的基础上，对多个灰度图像合成的大图像进行扩散。利用构造的分形迭代系统生成的混沌序列，将明文信息扩散加密，加密速度和效果良好。

本文算法包含三点创新性：(1) 构造了一个二维三参数的基于迭代函数系统的分形混沌系统，该系统具有优良的混沌性能；(2) 设计了Arnold变换相关的自适应性加密方法；(3) 提出了一个由密钥决定的Zigzag扫描方法，加密方法具有很好的明文相关性；(4) 算法使用三个阶段的混淆、置乱和扩散构成，适用于由任意大小组成的多个灰度图像的加密，加密效果和加密速度优良。本文的剩余部分安排如下：第二节构造了一个基于迭代函数系统的动力系统，并从数值上验证其混沌性能；第三节提出一个基于所构造的分形混沌系统的多图像的加密算法；第四节仿真实验验证了算法的性能；第五节为结语。

2. 单位正方型上的迭代函数系统及其分形混沌系统

2.1. 迭代函数系统

给定参数 $0<a,b,c<1$ ，考虑下面的迭代函数系统IFS $\left\{R^2;w_1,w_2,w_3,w_4\right\}$ ：

$w_1\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}0& a\\ b& 0\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)$ ； $w_2\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}c& 0\\ 0& 1-b\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}0\\ b\end{array}\right)$ ；

$w_3\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}1-c& 0\\ 0& 1-b\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}c\\ b\end{array}\right);w_4\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}0& -\left(1-a\right)\\ -b& 0\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}1\\ b\end{array}\right)$ (1)

由于 $0<a,b,c<1$ ，IFS(1)是一个压缩的迭代函数系统，因此在R²中有唯一的吸引子分形。由于单位正方形A = [0,1]²满足

$W\left(A\right)=w_1\left(A\right)\cup w_2\left(A\right)\cup w_3\left(A\right)\cup w_4\left(A\right)=A$

因此，IFS的分形为单位正方形A。用确定性算法，当选择不同的初始集合B₀进行迭代得到集合则有 $B_k=W^k\left(B_0\right),k=0,1,2,\cdots ,$ 分形几何理论知道。图1为参数取值为 $\left(0.5,0.5,0.5\right)$ 和 $\left(0.3,0.51,0.6\right)$ 迭代6次后得到的单位正方形内的填充曲线。可以看出曲线可以充满整个区域。

Figure 1. Filling curves generated by IFS (1)

图1. IFS (1)生成的填充曲线

该迭代函数系统的分形A上有一个伴随的位移动力系统 $\left\{A;S\right\}$ 。 $\left\{A;S\right\}$ 可以表示为公式(2)

$S\left(x,y\right)=\{\begin{array}{c}{w}_1{}^{-1}\left(x,y\right),\left(x,y\right)\in w_1\left(G\right)\left\{\left(a,y\right):0\le y\le b\right\},\\ w_2{}^{-1}\left(x,y\right),\left(x,y\right)\in w_2\left(G\right)\left\{\left(x,b\right):0\le x\le c\right\},\\ w_3{}^{-1}\left(x,y\right),\left(x,y\right)\in w_3\left(G\right)\left\{\left(c,y\right):b\le y\le 1\right\},\\ w_4{}^{-1}\left(x,y\right),\left(x,y\right)\in w_4\left(G\right)\left\{\left(x,b\right):a\le x\le 1\right\}.\end{array}$ (2)

分形几何理论表明S(A) = A，A是S的排斥吸引子，分形A上的动力系统 $\left\{G;S\right\}$ 具有混沌特性。推导得到定义在A上的动力系统 $S\left(x,y\right)$ 的具体形式如(3)所示：

$S\left(x,y\right)\{\begin{array}{l}\left(\frac{x}{b},\frac{x}{a}\right),\left(x,y\right)\in \left[0,a\right)\times \left[0,b\right],\\ \left(\frac{x}{c},\frac{y-b}{1-b}\right),\left(x,y\right)\in \left[0,c\right]\times \left(b,1\right],\\ \left(\frac{x-c}{1-c},\frac{y-b}{1-b}\right),\left(x,y\right)\in \left(c,1\right]\times \left[b,1\right],\\ \left(\frac{b-y}{b},\frac{1-x}{1-a}\right),\left(x,y\right)\in \left\{\left[a,1\right]\times \left[0,b\right),\left(a,b\right)\right\}.\end{array}$ (3)

2.2. 性能分析

下面选定初始值 $x_0=0.313$ ， $y_0=0.519$ ，系统参数 $a=0.313$ ， $c=0.633$ ， $b=0.718$ ，迭代10,000次，产生序列X，Y分析分形混沌系统(3)在单位正方形的混沌特性。

2.2.1. 轨迹分析

在给定的初始值和参数下，分形动力系统迭代产生的轨迹图和序列X，Y的时序图如图2所示。可以看出系统迭代生成的点，随机且均匀的遍布在整个单位正方形内，因此提出的动力系统(3)在相空间内，具有良好的遍历性。

Figure 2. Chaotic orbit and time series diagram

图2. 轨迹图和时序图

2.2.2. 分叉图

动力系统的分叉图是描述迭代序列随着参数选择的变化而变化的图像，分叉图可以直观地观察出生成序列的周期点或周期窗口，混沌序列的分叉图应该具有很好的随机性和遍历性，即具有较大的混沌参数范围。分别固定 $b,c$ 的值，参数取值范围 $a\in \left(0,1\right)$ ，刻画序列 $X,Y$ 随a的变换而变化的图像；固定 $a,c$ 的值，参数取值范围 $b\in \left(0,1\right)$ ，刻画序列 $X,Y$ 随b的变换而变化的图像；固定 $a,b$ 的值，参数取值范围 $c\in \left(0,1\right)$ ，刻画序列 $X,Y$ 随c的变换而变化的图像。所得结果如图3所示，可以看出分形混沌系统的混沌参数范围均是 $\left(0,1\right)$ 。

2.2.3. Lyapunov指数

混沌的动力系统应该具有对初始值的敏感性，Lyapunov指数可以衡量动力系统的初值敏感性，从而判断系统是否混沌。若一维的系统具有正的李雅普诺夫值，则认为系统是混沌的；n维的系统具有n个Lyapunov指数值，若n维系统的正的最大李亚普诺夫指数值大于0，则认为系统具有混沌特性，高维度的混沌系统具有超混沌性，一般来说，超混沌系统比一般的混沌系统具有更加复杂的动力学行为，在实际应用中更能满足加密信息的需求。下面探究当参数变换时，二维分形动力系统的Lyapunov指数变换，从而确定系统是否混沌，计算混沌范围。

对于二维混沌系统，

$\left(x_i,y_i\right)=f^i\left(x_0,y_0\right)$ (4)

$\left(x_0,y_0\right)$ 是系统的初始值，i表示迭代次数， $i=1,2,\cdots ,n$ 。将i时刻系统 $f\left(x,y\right)$ 的Jacobi矩阵记为J_i，迭代产生的Jacobi矩阵为 $J_n,J_{n-1},\cdots ,J_1$ ， $Q_0=I$ 。对矩阵 $J_n,J_{n-1},\cdots ,J_1$ ，进行QR分解：

$\begin{array}{c}qr\left[J_n,J_{n-1},\cdots ,J_1\right]=qr\left[J_n,J_{n-1},\cdots ,\left(J_1{Q}_0\right)\right]=qr\left[J_n,J_{n-1},\cdots ,\left(J_2{Q}_1\right)\right]\left[R_1\right]\\ =qr\left[J_n,J_{n-1},\cdots ,J_4\left(J_3{Q}_2\right)\right]\left[R_2{R}_1\right]\\ =\cdots =qr\left[J_n,J_{n-1},\cdots ,J_{i+1}\left(J_i{Q}_{i-1}\right)\right]\left[R_{i-1}{R}_{i-2}\cdots R_2{R}_1\right]\\ =\cdots =Q_n\left[R_n\cdots R_2{R}_1\right]=Q_{n}R\end{array}$

其中 $qr[\cdot ]$ 表示QR分解。从 $J_1$ 开始，在第i步中，得到用 $Q_{i-1}$ 右乘J_i得到 $B_i=J_i{Q}_{i-1}$ ，再做QR分解 $B_i=Q_i{R}_i$ ， $i=1,2,\cdots ,n$ 。得到阶数为2 × 2的对角矩阵 $R=R_n\cdots R_2{R}_1$ ，其对角线上的每一个元素是所有二阶矩阵R_i对角线上对应元素的乘积。

Figure 3. Bifurcation diagram

图3. 分叉图

混沌系统 $f\left(x,y\right)$ 的两个Lyapunov指数可以近似计算得到：

${\lambda }_1=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}|R\left(1,1\right)|},{\lambda }_2=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}|R\left(2,2\right)|}$ (5)

分别固定 $b,c$ 的值，参数取值范围 $a\in \left(0,1\right)$ ，刻画Lyapunov指数值随a的变换而变化的图像；固定 $a,c$ 的值，参数取值范围 $b\in \left(0,1\right)$ ，刻画Lyapunov指数值随b的变换而变化的图像；固定 $a,b$ 的值，参数取值范围 $c\in \left(0,1\right)$ ，刻画Lyapunov指数值随c的变换而变化的图像。如图4所示，加入参数c前后分形动力系统Lyapunov指数值的图像，可以明显的看出，当三个参数变化时，产生序列的Lyapunov指数值都大于0，所以分形混沌系统具有超混沌性。相比较本文的分形混沌系统具有更大的Lyapunov指数值，一些序列在整个区间包括端点的值都很稳定，具有更好的初值敏感性。

Figure 4. Lyapunov exponent

图4. Lyapunov指数

2.2.4. 自相关分析

时间序列 $X=\left\{x_1,\cdots ,x_n\right\}$ 与其k期延迟序列 $Y=\left\{x_{k+1},\cdots ,x_{n+k}\right\}$ 的自相关系数 ${\rho }_k$ 的定义如下：

${\rho }_k=\frac{{\displaystyle \underset{i=k+1}{\overset{n}{\sum }}\left(x_i-\bar{x}\right)\left(x_{i-k}-\bar{x}\right)}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(x_i-\bar{x}\right)}}$ (6)

其中n表示序列的长度，k是延迟期数，x_i表示时间序列在时刻i的值， $\bar{x}$ 表示序X的平均值。自相关系数是描述一个序列X与其延迟了k期的序列之间的相关程度，可以用来刻画序列的噪声性质。具有伪随机性的序列类似白噪声，序列的其延迟了 $k\left(k\ne 0\right)$ 期的自相关系数应该接近0，也就是序列与自身的时间序列的相关系数图与δ函数一致 [12] 。如图5所示，分形混沌系统产生的两个序列 $X,Y$ 与自身的时间延迟序列的自相关系数图与δ函数极其近似，说明产生的序列具有很好的伪随机性、混沌性。

Figure 5. Autocorrelation coefficient diagram

图5. 自相关系数图

2.2.5. 信息熵

信息熵(香农熵)反应了信息的不确定性，一般来说，熵值越大，不确定性越大。可以用来描述混沌映射生成的序列的随机性 [27] ，序列X的信息熵可以定义为：

$H=-{\displaystyle \underset{i=1}{\overset{n}{\sum }}p\left(x_i\right)}{\log }_{2}p\left(x_i\right)$ (7)

其中， $p\left(x_i\right)$ 表示x_i出现的概率，n表示将序列分割的区间个数，在实验中，设置n = 256。计算理论值时，假设序列中的随机数落入区间内的概率都相等，则信息熵的理论值 $H={\log }_{2}256=8$ 。如图6所示，当三个参数值分别变化时，信息熵值都接近序列的最大值，所以序列的随机性良好 [4] 。

Figure 6. Information entropy diagram

图6. 信息熵图

2.2.6. 0~1测试

0~1测试可以用来检验非线性系统是否处于混沌状态 [28] 。考虑动力系统生成的序列 $\varphi \left(n\right),n=1,2,\cdots ,N$ ，计算相关系数：

$K=\text{corr}\left(\alpha ,\beta \right)=\frac{\mathrm{cov}\left(\alpha ,\beta \right)}{\sqrt{\mathrm{var}\left(\alpha \right)\mathrm{var}\left(\beta \right)}}\in \left[-1,1\right]$ (8)

其中 $\alpha =\left(1,2,\cdots ,n\right)$ ， $\beta =\left(D\left(1\right),D\left(2\right),\cdots ,D\left(n\right)\right)$ ，

$\mathrm{cov}\left(\alpha ,\beta \right)=\frac{1}{n}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\left(\alpha \left(j\right)-\bar{\alpha }\right)}\left(\beta \left(j\right)-\bar{\beta }\right)$ ， $\bar{\alpha }=\frac{1}{n}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\alpha \left(j\right)}$ ， $\mathrm{var}\left(\alpha \right)=\mathrm{cov}\left(\alpha ,\alpha \right)$ ，

$D\left(n\right)=M\left(n\right)-\left(\underset{N\to \infty }{\lim }\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\varphi {\left(j\right)}^2}\right)\frac{1-\cos \left(nc\right)}{1-\cos \left(c\right)}$ ， $c\in \left(0,\pi \right)$ ，

$M\left(n\right)=\underset{N\to \infty }{\lim }\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\left({\left[p\left(j+n\right)-p\left(j\right)\right]}^2+{\left[q\left(j+n\right)-q\left(j\right)\right]}^2\right)}$ ，

$p\left(n\right)=\frac{1}{n}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\varphi \left(j\right)\cos \left(jc\right)}$ ， $q\left(n\right)=\frac{1}{n}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\varphi \left(j\right)\sin \left(jc\right)}$ 。

只有当 $K\approx 1$ 时，认为非线性系统是混沌的，迭代次数 $N=1000$ ， $c\in \left(\pi /5,4\pi /5\right)$ 。对分形混沌系统，可以从图7看到，当控制参数在 $\left(0,1\right)$ 区间变化时，相关系数 $K\approx 1$ ，说明序列的混沌性能良好。

Figure 7. 0~1 test

图7. 0~1测试

2.2.7. 排列熵

排列熵 [29] 是一种动力系统的分析标准，能够准确地揭示混沌映射迭代序列的伪随机性。设迭代序列为 ${\left\{x_t\right\}}_{t=1}^T$ 加入嵌入维度n和时间延l = 1，将序列{x_t}按递增排序： $x_t=\left[x_t,x_{t+l},\cdots ,x_{t+\left(n-1\right)l}\right]$ 。对于n维向量可能出现的排序方式有n!种，π表示排列方式， $f\left(\pi \right)$ 表示它在时间序列中的概率，每种排序出现的概率

$P\left(\pi \right)=\frac{f\left(\pi \right)}{T-\left(n-1\right)l}$ ，排列熵的定义为： $H_p\left(\pi \right)=-{\displaystyle \underset{n=1}{\overset{n!}{\sum }}P\left(\pi \right)}\ln P\left(\pi \right)$ ，归一化为：

$H_p=\frac{H_p\left(n\right)}{\ln \left(n!\right)}\in \left[0,1\right]$ (9)

排列熵越小，输出值越有规律；反之，排列熵越大，输出值的随机性越大。当排列熵接近1时，混沌性能最好。测试序列长度N = 2000，嵌入维度n = 3，时间延迟l = 1。结果如图8所示，当参数在区间 $\left(0,1\right)$ 变化时，排列熵稳定接近1，因此分形系统生成的序列伪随机性较好，表现出较强的混沌行为。

2.2.8. 统计性能分析

先对分形位移动力系统迭代10⁶次生成的混沌序列X和Y进行预处理：

$X=\text{mod}\left(\text{floor}\left(X,{10}^9\right),2\right),Y=\text{mod}\left(\text{floor}\left(Y,{10}^{11}\right),2\right),$ (10)

生成0~1比特流，再对比特流序列进行NIST SP800-22测试，测试这段数据的统计特性，以此反应数据的随机性。生成P值期望值为[0.1,1]。实验结果如表1所示，分形系统生成的序列满足所有测试的标准。

Figure 8. Permutation entropy diagram

图8. 排列熵图

Table 1. NIST SP800-22 test

表1. NIST SP800-22测试结果

3. 图像加密算法

本文提出的加密算法包括三个阶段，下面介绍加密时要用到的算法。

3.1. Arnold变换

传统的Arnold变换，是将需要加密的图像和参数矩阵通过矩阵的乘法运算，进行加密。改进的Arnold映射是将需要置乱的图像矩阵展开为一个一维向量，其交换位置只需要考虑列标的变动，不需要考虑行标的变动，因为此时所有像素点位置的行标皆为1。对向量中的任意一点坐标 $\left(1,j\right)$ 进行Arnold映射置换后得到新的坐标位置 $\left(p,q\right)$ ，其中 $a,b$ 是参数，公式如下：

$\left[\begin{array}{c}p\\ q\end{array}\right]=\left[\begin{array}{cc}1& a\\ b& ab+1\end{array}\right]\left[\begin{array}{c}1\\ j\end{array}\right]$ (11)

3.2. Zigzag扫描

Zigzag扫描方法是利用Z字形的扫描顺序扫描矩阵，将得到的一维向量重新放入新矩阵的置乱方法。 本文提出了一种由密钥决定扫描起点和长度的Zigzag扫描方法。密钥由多个灰度图像的Hash函数决定，生成每个图像不同的扫描始终点。如图9所示，选择需要加密矩阵的节点(又密钥量化之后的值决定)，将元素选择之后进行Z字形变换，再按照节点的顺序放入原矩阵。最后加密完成的各个矩阵合成大矩阵，长度不足的地方再后方补0，方便后续的计算。这种扫描方式比起传统的变换具有很好的明文敏感性。

Figure 9. Flow diagram of Zigzag transformation determined by secret key

图9. 密钥决定的Z字形变换流程图

3.3. 有限域扩散

基于有限域上的多项式运算生成的数域GF(Pⁿ)称为Pⁿ阶的伽罗瓦域，其中P是素数。伽罗瓦域中有四则运算的定义，其中GF(257)域的算术运算以257为模，剔除会造成信息损失的元素0之后，集合包含256个元素，与数字图像像素点的灰度值的取值范围 $\left[0,255\right]$ 刚好对应。所以GF(257)适合用于图像加密的扩散算法。将加法和乘法作为加密表，减法和除法作为解密表，使用查表法进行加密。运算速度快且加密效果优良。

3.4. 图像加密算法

本文提出了一个基于分形混沌映射图像加密算法，算法的流程如图10所示。使用分形混沌系统生成的序列决定各个环节的加密，经过Arnold变换，Zigzag扫描和有限域扩散后得到加密后的图像。

Figure 10. Flow chart of encryption algorithm

图10. 加密算法流程

1) 图像预处理

输入多个灰度图 $P_1,P_2,\cdots ,P_k$ ，先将组合成为大矩阵：

$P=\left(\begin{array}{cccc}{P}_1& P_2& \cdots & P_{k_2}\\ P_{k_2+1}& P_{k_2+2}& & P_{2k_2}\\ & ⋮& \ddots & ⋮\\ P_{\left(k_1-1\right)k_2+1}& P_{\left(k_1-1\right)k_2+2}& \cdots & P_{k_1\times k_2}\end{array}\right)$ (12)

其中 $k=k_1\times k_2$ ，P_i图像的大小是M_i × N_i，P的大小为M × N。

2) 密钥和混沌序列生成

利用大矩阵P和Hash函数生成256位的密钥 $k=k_1{k}_2\cdots k_{32}$ ，选择动力系统的初值和参数：初始值 $x_0=0.313$ ， $y_0=0.519$ ，系统参数 $a=0.313$ ， $c=0.633$ ， $b=0.718$ ，迭代 $M\times N+t$ 次，过渡态长度t。得到两个长度为M × N的混沌序列X和Y。

3) 自适应性的Arnold扰乱

设P₁图像的大小是M_i × N_i选择序列 $S=X\left(1:M_1\times N_1\right)$ ，

$S_1=\mathrm{mod}\left(\text{floor}\left(S\times {10}^{13}\right),256\right)$ (13)

重新组合成为大小为M₁ × N₁的矩阵S₁，对P₁进行一次Arnold变换后，组合扩散如下：

$R_1=\mathrm{mod}\left(\text{Arnold}\left(P_1,1\right)+S_1,256\right)$ (14)

再用R₁和P₂异或运算扰乱下一个图像P₂， $P_2=\text{bitxor}\left(R_1,P_2\right)$ 用扩散之后的P₂加密 $P_3=\text{bitxor}\left(P_2,P_3\right)$ 以此类推，直到

$P_k=\text{bitxor}\left(P_{k-1},P_k\right)$ (15)

最后再用P_k改变P₁的像素值： $P_1=\text{bitxor}\left(P_1,P_k\right)$ 。混淆每块图像的像素值，在有限次的变换后实现良好的加密效果。基于异或加运算和Arnold变换的扰乱和图像的明文相关，具有明文敏感性。

4) Zigzag型扫描

使用明文相关的密钥决定每个图像的扫描起点和长度，以k = 4为例

$D=k_{21}+k_{22}+k_{23}+k_{24}+k_{25}+k_{26};$

$d_1=\text{mod}\left(D\times 2^8+k_{21}+k_{24},M_1\times N_1\right);$

$d_2=\mathrm{mod}\left(d_1+D\times 2^7+k_{22}+k_{25},M_1\times N_1\right);$

$d_3=\mathrm{mod}\left(d_2+D\times 2^6+k_{23}+k_{26},M_1\times N_1\right);$

得到三个扫描起点 $d_1,d_2,d_3$ 。按从小到大排序，分别作为各个图像的扫描长度。若 $d_1<d_2<d_3$ ，则图像P₁从 $\left(1,1\right)$ 开始扫描，扫描长度为d₁，图像P₂从d₁开始扫描，扫描长度为 $d_2-d_1,\cdots ,$ 以此类推，直到图像P₄被扫描完成。再将所有的图片按照顺序组合为密文图像C₁使用明文图像Hash函数生成的密钥决定了每幅图像Zigzag型变换的起点和终点，使算法具有更好的明文敏感性和密钥敏感性。

5) 有限域扩散

整数集合 $\left\{0,1,\cdots ,P-1\right\}$ 按通常的加法(减法)、乘法(除法)代数运算在模P意义下构成一个有限域GF(P)。即：对于加法，如果 $a,b\in \text{GF}\left(P\right)$ ，则 $a+b\equiv r\left(\mathrm{mod}P\right)$ ， $r\in \text{GF}\left(P\right)$ ；对于乘法，如果 $a,b\in \text{GF}\left(P\right)$ ，则 $a\cdot b\equiv s\left(\mathrm{mod}P\right)$ ， $s\in \text{GF}\left(P\right)$ 。本文构造了GF(257)的乘法表和除法表，剔除元素0，用乘法表加密，除法表解密。能够达到无损解密。选择序列 $T=Y\left(M\times N+1:2\times M\times N\right)$ ，

$S_2=\mathrm{mod}\left(\text{floor}\left(T\times {10}^{10}\right),256\right)$ (16)

将C₁展开成为一维向量，利用混沌序列S₂，进行正向和逆向的扩散。将明文信息扩散加密，加密速度和效果良好。

算法的解密过程为加密的逆过程。

4. 安全性分析

4.1. 仿真实验结果

(a) (b) (c)

Figure 11. Encryption and decryption test results. (a) Plaintext image; (b) cipher text image; (c) decrypt image

图11. 加密及解密测试结果。(a) 明文图像；(b) 密文图像；(c) 解密图像

本文使用MATLAB2017b进行实验仿真，选择大小为512 × 512的Lena、Peppers、Baboon和Barbara进行加密和解密，结果如图11所示。加密图像呈现随机均匀分布的噪声，利用算法图像能够无损解密。

4.2. 密钥敏感性分析

密钥敏感性分析是分析同一明文图像在密钥发生微小变化时，得到的两个密文图像之间的差别。如果两个密文图像之间的差别显著，则表示图像密码系统的密钥敏感性强；若两个密文图像之间差别较小，则称密钥敏感性差。良好的图像密码系统应该具有颇强的密钥敏感性。

衡量两幅相同大小图像 $P_1,P_2$ 的差别有两种方法，分别是计算图像的NPCR和UACI值。NPCR是计算两幅图像之间不同像素点值所占的比例：

$\mathrm{NPCR}\left(P_1,P_2\right)=\frac{1}{M\times N}{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\left|\text{Sign}\left(P_1\left(i,j\right)-P_2\left(i,j\right)\right)\right|}}\times 100\%$ (17)

其中Sign是符号函数， $\text{Sign}\left(x\right)=\{\begin{array}{l}1,x>0\\ 0,x=0\\ -1,x<0\end{array}$ 。UACI是计算两幅图像不同像素点值的差别程度：

$\text{UACI}\left(P_1,P_2\right)=\frac{1}{M\times N}{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\left|\frac{P_1\left(i,j\right)-P_2\left(i,j\right)}{255}\right|}}\times 100\%$ (18)

理论值分别为NPCR = 99.6094%，UACI = 33.4635%。

对组合之后的大图像进行密钥敏感性测试。本文的密钥由外部参数 $x_0,y_0,a,b,c$ 和加密算法的密钥 $d_1,d_2,d_3$ 组成。对密钥分别改变10⁻¹⁵单位，分别计算同一明文图像加密之后的密文图像的NPCR和UACI，实验结果如表2所示。可以看出NPCR和UACI的测试值均接近理论期望值，表明本文提出的算法具有良好的密钥敏感性。

Table 2. Key sensitivity test

表2. 密钥敏感性测试

4.3. 密钥空间

Table 3. Comparison of key space

表3. 密钥空间大小比较

密钥空间是指所有合法密钥可构成的空间，图像密码系统的密钥空间应该足够大，从而可以有效地对抗穷举攻击。本文混沌系统的控制参数为 $K=\left\{x_0,y_0,a,b,c\right\}$ ，初值和参数的取值范围为 $\left(0,1\right)$ 。测试中计算精度为10⁻¹⁵，所以算法的密钥空间大小为 ${10}^{15\times 5}={10}^{75}\gg {10}^{30}$ 其中 $2^{100}\approx {10}^{30}$ 是能够抵抗穷举攻击的最小值。所以本文的算法可以抵抗暴力攻击。与其他算法的密钥空间比较如表3所示。

4.4. 直方图统计分析

方差和标准差是度量数据的离散度，支撑直方图用来衡量图像数据，图形直方图越均匀，直方图的方差越低。将区间分为256份，方差计算公式为：

$\alpha =\frac{1}{256}{\displaystyle \underset{i=1}{\overset{256}{\sum }}\left(x_i-\bar{x}\right)}$ (19)

其中，均值 $\bar{x}=\frac{M\times N}{256}$ ，图像大小为M × N，x_i为直方图中各值出现的频率。 图像的标准差为 $\beta =\sqrt{\alpha }$ 。

对明文图像和密文图像分别画出图像直方图，从图12可以看出，加密后的图像变得非常均匀，破坏了图像的统计特征。

4.5. 相关性分析

数字图像的特点就是相邻像素之间有较大的相关性，加密算法的效果希望能够消除这种相关性。从需要考察的图像中任取N = 3000对相邻的像素点。记它们的灰度值为 $\left(u_i,v_i\right)$ ， $i=1,2,\cdots ,N$ ，则向量u = {u_i}和v = {v_i}间的相关系数计算如下：

$r_{uv}=\frac{\mathrm{cov}\left(u,v\right)}{\sqrt{D\left(u\right)}\sqrt{D\left(v\right)}},\mathrm{cov}\left(u,v\right)=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left(u_i-E\left(u\right)\right)\left(v_i-E\left(v\right)\right)}$ (20)

$D\left(u\right)=\frac{1}{N}{\left(u_i-E\left(u\right)\right)}^2,D\left(v\right)=\frac{1}{N}{\left(v_i-E\left(v\right)\right)}^2$ (21)

Figure 12. Histogram of plaintext image and cipher text image

图12. 明文和密文图像的直方图

分析明文图像(Lena, Peppers, Baboon, Barbara)和密文图像水平、垂直和对角方向上的像素间的相关性，加密前后图像的相关系数的计算结果如表4示。实验结果表明，明文图像水平、垂直和对角方向的相关系数接近1，而密文图像在各个方向的相关系数接近0，本文的加密算法可以有效地破坏相邻像素之间地相关性，保护明文图像的信息。以Lena图像为例，密文图像相关系数的比较如表5示。

Table 4. Correlation coefficient of adjacent pixels

表4. 相邻像素的相关系数

Table 5. Correlation coefficient of Lena cipher text image

表5. Lena图像密文图像的相关系数比较

4.6. 信息熵分析

信息熵可以反应图像信息的不确定性。信息熵越大，不确定性越大，可视信息越少，信息熵的计算公式为：

$H=-{\displaystyle \underset{i=1}{\overset{L}{\sum }}p\left(i\right)}{\log }_{2}p\left(i\right)$ (22)

其中L为图像的灰度等级数，p(i)表示灰度值i出现的概率。在实验中，对于L = 256的灰度随机图像，信息熵的理论值H = log₂256 = 8。明文图像和加密图像的信息熵如表6所示，算法的信息熵接近8。因此，通过分析加密图像的信息来获取明文信息近乎不可能。

Table 6. Information entropy of plaintext image and cipher text image

表6. 明文图像和密文图像的信息熵

4.7. 差分攻击分析

差分攻击是在明文信息发生细微变化时，使用同一密钥得到的两幅密文图像之间的差别，如果两个密文图像的差别迥异，则称该图像密码系统具有优良的明文敏感性；如果这两个密文图像的差别不大，则称该图像密码系统具有较弱的明文敏感性。使用NPCR和UACI的值来衡量图像之间的差异。

本文对明文图像随机选取200个像素，在每一次实验中只改变一个像素位置灰度值一个单位，用相同的密钥进行加密，得到密文图像。计算原始密文图像和新的密文图像的NPCR和UACI值，计算其平均值。测试结果NPCR = 99.6031%，UACI = 33.4500%。NPCR和UACI的测试值均接近理论期望值，表明本文提出的算法具有良好的抵抗差分攻击的能力。

4.8. 鲁棒性分析

在真实应用环境中，图像信息传输时，会受到噪声信息干扰或者损失部分图像信息，一个好的算法要能够拥有抵抗这些攻击的能力。

4.8.1. 噪声攻击

在密文图像中加入不同强度的椒盐噪声，并用正确的密钥解密，结果如图13所示。结果表明，解密后的图像仍然可以被识别，说明算法能够抵抗噪声攻击。

4.8.2. 裁剪攻击

(a) 噪声强度0.01的密文图像 (b) 噪声强度0.05的密文图像(c) 噪声强度0.1的密文图像(d) 噪声强度0.01的解密图像 (e) 噪声强度0.05的解密图像 (f) 噪声强度0.1的解密图像

Figure 13. Cipher text images and decrypted images after pepper and salt noise attack

图13. 椒盐噪声攻击后的密文图像和解密后图像

对密文图像进行不同程度的裁剪攻击，用正确的密钥解密，实验结果如图14所示。明文图像能够得到较大程度的还原，说明本文的算法能够抵抗裁剪攻击。

(a) 裁剪1/16的密文图像(b) 裁剪1/8的密文图像 (c) 裁剪1/4的密文图像(d) 裁剪1/16的解密图像 (e) 裁剪1/8的解密图像 (f)裁剪1/4的解密图像

Figure 14. Cipher text images and decrypted images after clipping attack

图14. 裁剪攻击后的密文图像和解密后图像

4.8.3. PSNR和SSIM

均方误差(MSE)和峰值信噪比(PSNR)可以用来评价图像的质量，PSNR表示一副图像加入噪声后和原图像的差距程度。PSNR越高，图像与原始图像差别越小，效果越好。均方差用MSE表示，是各数据偏离真实值的距离平方和的平均数，均方差值越小越好。对大小为M × N，数据类型为uint8的原始图像X和处理后的图像Y，灰度值最大值为255，均方误差定义为：

$\text{MSE}\left(X,Y\right)=\frac{1}{M\times N}{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left|X\left(i,j\right)-Y\left(i,j\right)\right|}^2}}$ (23)

峰值信噪比定义为：

$\text{PSNR}=10\times {\log }_{10}\left(\frac{255\times 255}{\text{MSE}\left(X,Y\right)}\right)$ (24)

与MSE和PSNR衡量两幅数字图像之间的绝对差异不同，SSIM模拟人眼的视觉感受，考虑图像的亮度、对比度和结构来评价图像数据 [34] ：

$\text{SSIM}\left(X,Y\right)=\frac{\left(2{\mu }_X{\mu }_Y+C_1\right)\left(2{\sigma }_{XY}+C_2\right)}{\left({\mu }_X^2+{\mu }_Y^2+C_1\right)\left({\sigma }_X^2+{\sigma }_Y^2+C_2\right)}$ (25)

其中μ_X和μ_Y分别表示图像X和Y的样本均值，σ_X和σ_Y表示图像X和Y的标准差，σ_XY为X和Y之间的协方差，C₁和C₂是分母趋于零时稳定除法的常数，一般取C₁ = 0.01，C₂ = 0.03。SSIM的取值范围是 $\left[0,1\right]$ ，当X = Y时SSEM = 1。

将测试的明文图像加密之后加入噪声或剪切攻击，再用同一密钥解密，计算原明文图像和差分攻击后的明文图像之间的PSNR和SSIM值，用来衡量收到攻击之后明文图像信息的损失程度。测试结果如表7所示，本文提出的算法得到的PSNR值较大，高于30 dB；在裁剪攻击下的SSIM值都较大，甚至接近1；说明明文图像失真程度较小，两个图像在人眼观察下区别较小，本文提出的算法能够抵抗噪声和裁剪攻击，其中抗裁剪攻击能力较好。以Lena图像为列，与其他算法的PSNR和SSIM值比较为表8所示。

Table 7. PSNR (dB) and SSIM of the image

表7. 图像的PSNR (dB)和SSIM

Table 8. Comparison of PSNR and SSIM values (512 × 512 Lena image)

表8. PSNR和SSIM值比较(512 × 512的Lena图像为例)

4.9. 加密及解密速度

实验使用一台配置AMD R7-5800H、16 G内存的笔记本电脑，算法基于MATLAB2017b语言实现。加密及解密速度测试使用尺寸为1024 × 1024的图像。加密速度测试中，对同一图像加密100次，记录每次加密耗时，取计算平均值，解密速度测试也采取相同的方式。测试结果算法加密平均速度为0.2081773 s，解密平均速度为0.2347666 s。这说明算法具有有效的加密及解密速度，可适用于实际通信环境。与其他算法的运算时间比较如表9，本文的算法具有优势。

Table 9. Comparison of encryption and decryption time

表9. 运算时间比较

5. 结语

本文构造了一个新的分形混沌系统，并且从数值上的分析验证了该动力系统的混沌特性，改进之后的动力系统具有更好的混沌特性。利用该分形位移动力系统所生成的混沌序列，设计了一个针对多个灰度图像的扩散–置乱–扩散的加密算法。加密算法包括三个阶段，第一阶段利用改进的Arnold变换将首图像加密，再利用前一个图像加密下一个图像，以此类推，将混沌序列的信息扩散；第二阶段用明文图像的密钥决定扫描起点和长度的Zigzag变换置乱图像；第三阶段使用Galois有限域扩散多个图像组成的大图像。实验结果表明，该算法能有效抵抗蛮力攻击，统计攻击和差分攻击，也能抵抗噪声和裁剪攻击，具有优良的安全性能，而且计算速度较快。

致谢

本文为广东省基础与应用基础研究基金资助项目(No. 2020B1515310018, 2023A1515030199)。

NOTES

^*通讯作者。