基于新的五维多环多翼超混沌系统的图像加密算法*

庄志本 李军 刘静漪 陈世强

1) (湖北民族大学理学院, 恩施 445000)

2) (湖北民族大学信息工程学院, 恩施 445000)

3) (湖北民族大学新材料与机电工程学院, 恩施 445000)

本文提出了一种基于新的五维多环多翼超混沌系统的数字图像加密方法.首先, 将明文图像矩阵和五条混沌序列分别通过QR分解法分解成一个正交矩阵和一个上三角矩阵, 将混沌系统产生的五条混沌序列分别通过LU分解法分解成一个上三角矩阵和一个下三角矩阵, 分别将两个上三角矩阵和一个下三角矩阵相加,得到五个离散后的混沌序列; 其次, 将明文图像矩阵分解出来的正交矩阵与五个混沌序列分解出来的五个正交矩阵相乘, 同时把明文图像矩阵分解出来的上三角矩阵中的元素通过混沌序列进行位置乱, 再将操作后的两个矩阵相乘; 最后, 将相乘后的矩阵通过混沌序列进行比特位位置乱, 再用混沌序列与其进行按位“异或”运算, 得到最终加密图像.理论分析和仿真实验结果表明该算法的密钥空间远大于10200, 密钥敏感性强, 能够有效地抵御统计分析和灰度值分析的攻击, 对数字图像的加密具有很好的加密效果.

1 引 言

在网络高速发展的今天, 保密通信已成为当今备受关注的一个热门话题, 其中混沌系统和超混沌系统的复杂结构和动力学行为在通信保密中有着很好的应用前景.因此, 构造出能够产生具有复杂拓扑结构的超混沌系统是非常重要的.

目前, 构造新混沌系统的方法主要有三种, 一是通过引入一些非线性函数来产生多涡卷混沌吸引子[1−4]; 二是在三维混沌系统的基础上, 增加一个新的状态反馈控制器来构造出超混沌系统[5−7],这就使得到的新混沌系统满足超混沌系统的必要条件.在构造的过程中, 必须将新增加的控制器反馈到原始的控制器中, 同时将原始控制器反馈到新的控制器中, 这两个操作能够使此系统的四个控制器之间互相影响, 可使得关系更加复杂; 三是通过多个正李氏指数来构造超混沌系统[8].1979年,Rossler首次提出了超混沌系统, 从此, 专家和学者就相继提出了一些超混沌系统.2014年与2015年,彭再平等[6]和刘杨[7]分别提出的新的四维超混沌系统结构简单, 动力学特征复杂, 但不具有多环的特征, 且刘杨的系统方程处于混沌状态的参数不够大.本文提出了一种五维多环多翼超混沌系统, 能够产生明显的多环和多翼, 系统处于混沌状态的参数范围大, 且系统方程结构简单, 没有平衡点.

随着混沌理论及其应用的发展, 专家和学者提出了很多基于混沌的数字图像加密算法.相对于传统的加密算法, 基于混沌的图像加密算法在安全性、抗攻击能力、复杂度等方面都具有很好的特性[9,10].因此, 专家和学者提出了许多基于混沌的图像加密算法.为了提高后续密码系统的安全性, 2016年,Zhang等[11]提出基于3D比特位矩阵的混合置换架构; 2018年, Luo 等[12]根据密钥和明文图像信息, 设计了一种新颖的并行量化方法.而在图像的加密过程中, 主要是改变其像素的位置和像素值的大小, 目前主要有三种方法, 一是在比特位上用混沌序列进行行列扰乱, 二是直接对明文图像进行像素位置扰乱和改变像素值的大小[13−15], 三是用可逆矩阵去乘明文图像矩阵.针对方法一, 2018年,He等[16]提出对比特位高位进行扰乱, Raza和Satpute[17]提出先对高6位进行处理, 紧接着再对低6位进行扰乱, 但所用的混沌系统都不是超混沌系统, 也不能产生多环; 针对方法三, 2017年, Ahmad等[18,19]提出用正交矩阵去乘以明文图像矩阵, 以达到改变像素值的目的, 但不能很好地克服计算机精度对复杂的混沌动力学特性迅速退化问题.在本文提出的方案中, 当混沌系统的参数 f ∈(0,50) 时,混沌系统都处于混沌状态; 将原文图像进行 Q R 分解, 分别对正交矩阵 Q 和上三角矩阵 R 进行处理;同时将混沌系统产生的5条混沌序列进行了 Q R 分解和 L U 分解处理, 通过这些处理, 增大了穷举攻击的难度, 除穷举之外的一切基于明确的明文密文映射的攻击方法, 对该方法都将失效, 具有较高的安全性.

2 相关知识

2.1 矩阵的一些基本性质及运算规律

设 Q 是一个正交矩阵, P 为与 Q 满足乘法条件的一般矩阵, 那么:

i) Q 的逆一定存在并且有 Q−1=QT;

ii) Q−1Q=QTQ=QQ−1=QQT=E, (其中 E 为单位矩阵);

iii) 如果 A =QP , 那么 P =Q−1A=QTA ;

iv) 如 果 Q1,Q2,...,Qn都 是 与 Q 具 有 相 同 行列 的 正 交 矩 阵 , 且 A1=Q1Q2...QnP , 则

2.2 矩形矩阵的正交分解与一般方阵的高斯消去法分解

矩形矩阵的正交分解又称 Q R 分解.Q R 分解是把一个 m ×n 的矩阵 A 分解成一个正交矩阵Q和一个上三角矩阵 R .即设矩阵

通过调用MATLAB程序 [ Q ,R]=qr(A) , 就可以把矩阵 A 分解成一个正交矩阵

和一个上三角矩阵

很显然矩阵 R 是与矩阵 A 具有相同大小的上三角矩阵, 矩阵 Q 是 m ×m 矩阵, 并且还要满足

高斯消去法分解又称 L U 分解.L U 是把任意一个方阵 X 分解成一个下三角矩阵 L 和一个上三角矩阵 U .即设矩阵

通过调用MATLAB程序 [ L ,U]=lu(X) ,就可以把矩阵 X 分解成一个下三角矩阵

和一个上三角矩阵

很显然矩阵 L 和 U 都与方阵 X 具有相同的大小.并且满足

3 新的五维多环多翼超混沌系统

1994年,Sprott[20]总结了许多三维混沌系统,其A系统方程如下:

在系统(1)的基础上,我们引入了两个控制器 w ,v ,并将原始控制器 y 反馈到新的控制器 w 中,将原始控制器 z 反馈到新的控制器 v 中,这两个操作能够使此系统的五个控制器之间互相影响,可使得关系更加复杂.新构造出来的五维多环多翼超混沌方程如下:

方程中的 a,b,c,d,e,f,g 是系统参数，且均取大于0的数.

3.1 Lyapunov指数谱分析

混沌方程的主要动力学特性可以通过它的Lyapunov指数谱来进行分析.当系统参数 a =1.5 ,b=c=d=e=f=g=1时,该系统的 Lyapunov指数谱如图1所示.从图1中可以明显地看出,该系统有两个Lyapunov指数为正,一个Lyapunov指数为零,还有两个Lyapunov指数为负,因此系统(2)处于超混沌状态.

图1 Lyapunov 指数谱Fig.1.Lyapunov exponent diagram.

3.2 平衡点分析

很显然,无论系统参数 a ,b,c,d,e,f,g 取大于0的任何值,方程(3)都没有解,因此该系统方程无平衡点.

3.3 分岔图分析

方程参数 a =1.5 , b =c=d=e=g=1 ,参数f 在变化时,关于 f 的分岔图如图2所示.从图2中可以看出,当 f ∈(0,50) 时,系统 (2) 都是处于混沌状态.因此,该系统处于混沌状态的参数动态范围很大.

图2 系统 (2) 随 f 变化的分岔图Fig.2.Bifurcation diagram of system (2)variation with f .

3.4 时间序列图

当系统参数 a =1.5 ,b=c=d=e=f=g=1时, x ,y,z,w,v 的值随时间 t 的变化序列图如图3所示.从图3中也可以看出,在系统参数 a =1.5 ,b=c=d= e =f=g=1 时,系统 (2)处于混沌状态.

3.5 五维多环多翼混沌相图

当系统参数 a =1.5 ,b=c=d=e=f=g=1时,系统(2)产生的多环多翼的三维相图如图4所示.方程(2)产生的多环多翼的平面相图如图5所示.从图4和图5可以明显地看出,该混沌系统能在多个方向上产生多环多翼.

4 算法描述

4.1 加密算法描述

给定一个 M ×N 的灰度图像 P ,加密步骤如下:

Step1输入灰度图像 P ,通过 Q R 分解把P分解成一个正交矩阵 E 和一个 M ×N 的上三角矩阵 R ;

Step2输入超混沌系统的初始值y0=[1,0.1,2,0.5,0.4]和步长 l =0.02 ,求出迭代总时间 T =l×(P×P+ 250).其中 P =max(M,N) ;

Step3调用 ode45 函数,迭代系统 (3),产生5条混沌序列;

Step4分别把5条混沌序列作如下处理:

其中y(201:(200+M×N),i)是把混沌序列y(i)的第201个值到第 2 00+M×N 个值取出来;reshape(y(201:(200+M×N),i),M,N)是把取出来的 M ×N 个值转换成 M ×N 矩阵;

Step5第4步中的5个 M ×N 矩阵分别通过 Q R 分解,得到五个正交矩阵E1,E2,E3,E4,E5和五个上三角矩阵 R1,R2,R3,R4,R5;

Step6改变初始值 y0,重复步骤二到步骤五n(n≥2)次,得到 5n 个正交矩阵和 5n 个上三角矩阵.分别记为E11,E12,E13,E14,E15,E21,···,En1,En2,En3,En4,En5;

Step7第一步中的正交矩阵 E 与第六步中的正 交 矩 阵E11,E12,E13,E14,E15,E21,···,En1,En2,En3,En4,En5分别相乘,得到矩阵 E E ;即EE=E11×E12×E13× E14× E15×E21× ···×En1×En2×En3×En4×En5×E;

Step8第三步中的五条混沌序列分别作如下处理:

然后把这五个 P ×P 矩阵分别通过 L U 分解,得到五个上三角矩阵 L1,L2,L3,L4,L5和五个下三角矩阵 U1,U2,U3,U4,U5;

图3 时间序列图 (a) x -t 时间序列;(b) y -t 时间序列;(c) z -t 时间序列;(d) w -t 时间序列;(e) v -t 时间序列Fig.3.Time series diagram:(a) x -t time series;(b) y -t time series;(c) z -t time series;(d) w -t time series;(e) v -t time series.

Step9L1,L2,L3,L4,L5与U1,U2,U3,U4,U5分别对应相加,即LU1=L1+U1,LU2=L2+U2,LU3=L3+U3,LU4=L4+U4,LU5=L5+U5.然后分别取出 L Ui(i = 1,2,3,4,5) 中的前 M 行与前 N 列位置上的元素,从而得到五个 M ×N 矩阵.分别记为 M N1,MN2,MN3,MN4,MN5;

Step10MN1,MN2,MN3,MN4,MN5与R1,R2,R3,R4,R5R2,R3,R4,R5分别对应相加,即B41=MN1+R1,B42=MN2+R2,B43=MN3+R3,B44=MN4+R4,B45=MN5+R5;

Step11B41,B42,B43,B44,B45转换成 1行M×N列的矩阵,同时把每个矩阵中的所有元素转换成 0—255之间的整数,分别得到矩阵B51,B52,B53,B54,B55;

Step12B51,B52,B53,B54,B55作如下处理:

其中i2=1,2,3,···,M × N;

图4 三维相图 (a) x -y-z 三维图; (b) x -y-w 三维图; (c) x -y-v 三维图; (d) x -z-w 三维图; (e) x -z-v 三维图; (f) x -w-v 三维图;(g) y -z-w 三维图; (h) y -z-v 三维图; (i) y -w-v 三维图; (j) z -w-v 三维图Fig.4.Three-dimensional phase diagram: (a) x -y-z Three-dimensional map; (b) x -y-w Three-dimensional map; (c) x -y-v Threedimensional map; (d) x -z-w Three-dimensional map; (e) x -z-v Three-dimensional map; (f) x -w-v Three-dimensional map; (g)y-z-w Three-dimensional map; (h) y -z-v Three-dimensional map; (i) x -y-z Three-dimensional map; (j) z -w-v Three-dimensional map.

Step13B61,B62,B63,B64作如下处理:

其 中 m od 为 求 模 运 算 , r ound(f) 是 对 f 进 行 靠 近取整;

Step14上三角矩阵 R 的所有行的元素进行扰乱,扰乱公式如下:

其 中 Z (i3,:) 是 矩 阵 Z 的 第 i3 行 的 所 有 列 ,circshift(A,k,2)是把行向量 A 的所有元素按顺时针方向移动 k 个单位;

Step15将 R R 的所有列的元素进行扰乱,扰乱公式如下:

图5 二维平面相图 (a) x -y 平面; (b) x -z 平面; (c) x -v 平面; (d) y -z 平面; (e) y -w 平面; (f) z -w 平面; (g) z -v 平面Fig.5.Two-dimensional plane phase diagram: (a) x -y flat; (b) x -z flat; (c) x -v flat; (d) y -z flat; (e) y -w flat; (f) z -w flat; (g)z-v flat.

其 中 Z (:,i4) 是 矩 阵 Z 的 第 i4 列 的 所 有 行 ,circshift(B,k,1)是把列向量 B 的所有元素按顺时针方向移动 k 个单位;

Step16将 E E 与 R R1相乘, 得到矩阵 E E1,即 E E1=EE×RR1;

Step17将 E E1转换成 1 行 M ×N 列的矩阵EE2, 再把 E E2中的所有元素都映射成0—255之间的整数, 映射公式如下:

Step18第十七步中经过映射公式处理了的EE1转换成二进制数, 得到矩阵 E E2;

Step19E E2中的每一行二进制数分别进行扰乱, 得到矩阵 E E3, 扰乱公式如下:

Step20E E3中的每一列二进制数分别进行扰乱, 得到矩阵 E E4, 扰乱公式如下:

Step21将 E E4转换成十进制数, 再把已经转换成十进制数的 E E4转换成1行 M ×N 列的矩阵 E E5;

Step22E E5与 B65进行按位“异或”运算, 得到一个新的序列 E E6, 把 E E6转换成 M × N 的矩阵 E E7, E E7为最终的加密图像.

4.2 解密算法描述

Step1输入图像矩阵 E E7, 同时与 B65进行按位“异或”运算, 得到一个新的序列 E E8;

Step2E E8转换成二进制数, 然后把行列进行还原, 得到矩阵 E E9;

Step3E E9转 换 成 十 进 制 数 , 得 到 矩 阵EE10, 再把 E E10通过下列公式还原成 E E1,

其中 j =1,2,···,M × N ;

Step4在 E E1的左边依次乘上E11,E12,E13,E14,E15,E21,···,En1,En2,En3,En4,En5,E的逆, 得到矩阵 R R1;

Step5R R1的行列进行还原, 得到矩阵 R ,再把矩阵 E 与 R 相乘, 得到解密图像 P .

5 实验结果

5.1 实验平台

PC 机配置: Intel(R) Core(TM) i3-4170 CPU@ 3.70 GHz, 内存 4 GB, Windows7 32 位操作系统.通过 Matlab R2014 a 编写程序实现上述加密算法.

5.2 实验结果

实验选取了经典的 lena, baboon, boat灰度图像, 其大小均为512 × 512.明文图像、加密图像和解密图像如图6所示.

6 安全性分析

6.1 密钥空间分析

图6 数 字图像 加解密 实验图 (a) lena 原图; (b) lena 加 密图像; (c) lena 解 密图像; (d) baboon 原图; (e) baboon 加密图 像;(f) baboon 解密图像; (g) boat 原图; (h) boat 加密图像; (i) boat解密图像Fig.6.Digital image encryption and decryption experiment: (a) Original Lena image; (b) encrypted Lena image; (c) decrypted Lena image; (d) original baboon image; (e) encrypted baboon image; (f) decrypted baboon image; (g) original boat image.; (h) encrypted boat image; (i)decrypted boat image.

决定图像加密算法强度的最重要因素之一是密钥空间的大小.本文的初始密钥由 y0中的五个初始值、系统参数 a ,b,c,d,e,f,g 、正交矩阵 E 和原图像的最大值最小值组成, 以计算机精度为 1 0−15计算的话, 本算法的密钥空间远大于 1 0200> 2 100 .如果一种图像加密算法的密钥空间大于 2100, 则它就是安全的[21,22].因此本算法是足够安全的.

6.2 直方图分析

直方图可以很好地反映图像像素值的分布情况, 直方图越平坦则像素值分布就越均匀.图7分别是lena, baboon和boat的原图像直方图和加密后图像的直方图.

6.3 信息熵分析

信息熵是最重要的随意因素之一.计算公式如下:

这里的 p (mi) 是 mi的机率, L 是 mi的总数量.对于灰度图像来说, 信息熵的最大值为 8.Lena,baboon和boat加密前后的信息熵与文献[11−17]的测试值如表1所列.仿真实验结果表明本算法具有良好的加密效果.

图7 明文图像和密文图像直方图 (a) lena 明文直方图; (b) lena 密文直方图; (c) baboon 明文直方图; (d) baboon 密文直方图;(e) boat 明文直方图; (f) boat 密文直方图Fig.7.Histogram of plaintext and ciphertext images (a) Plaintext Lena image histogram; (b) ciphertext Lena image histogram; (c)plaintext baboon image histogram; (d) ciphertext baboon image histogram; (e) plaintext boat image histogram; (f) ciphertext boat image histogram.

表1 明文图像与加密图像的信息熵分析表Table 1.Information entropy analysis table of plain text and encrypted image.

6.4 不动点比和灰度平均变化值分析

不动点比为图像加密后灰度值未发生变化的像素点占所有像素点的百分比, 计算公式如(15)式所示; 而灰度平均变化值能更好地评价加密图像灰度变化的程度, 计算公式如(16)式所示.

其中 G 为明文图像, C 为密文图像.根据 (16) 式计算出本算法的灰度平均变化值如表3所列.

表2 加密图像不动点比分析表Table 2.Encrypted image fixed point ratio analysis table.

表3 灰度平均变化值分析表Table 3.Grayscale average change value analysis table.

6.5 密钥敏感性分析

密钥的微小改变就会导致密文的极大改变, 这就是密钥敏感性.实验以经典的lena图像为例.图8展示了本算法对初始密钥的敏感性.图8(a)是lena的明文图像; 图8(b)和图8(c)分别是用密钥y0=[1,0.1,2,0.5,0.4]和y1=[1,0.1000000000000001,2,0.5,0.4]加密的密文图像 Y1和 Y2; 图8(d)是Y1用 y0正确解密的结果图; 图8(e) 和图8(f) 是 Y1和Y2分别用错误的密钥 y1和 y0解密的结果.图8说明尽管密钥 y0和 y1之间仅有微小的变化, 密文图像 Y1和 Y2却不能相应的用密钥 y1和 y0正确解密.

两幅图像之间的差异还可以用像素变化率(NPCR)和归一化平均变化强度(UACI)来度量.计算公式分别如下:

6.6 相邻像素相关性分析

相邻像素的相关性用来评价图像加密算法在消除明文图像相邻像素相关性方面的效果.本文分别在lena, baboon和boat的明文图像和密文图像中随机选取了5000个相邻像素点, 按(19)—(22)式计算明文图像与密文图像在水平方向的相邻像素的相关性系数、垂直方向的相邻像素的相关性系数和对角线方向的相邻像素的相关性系数.测试结果如表5所列.Baboon明文图像和密文图像在水平方向、垂直方向和对角线方向的像素相关性图如图9所示.

图8 密钥敏感性测试图 (a) 明文图像; (b) 密文 Y1 (密钥为 y0); (c) 密文 Y2 (密钥为 y1); (d) Y1 正确解密结果; (e) Y1 用 y1 错误解密结果; (f) Y2 用 y0 错误解密结果Fig.8.Key sensitivity tests: (a) Plain-image; (b) cipher Y1 with key y0 ; (c) cipher Y2 with key y1 ; (d) right decrypted Y1 ;(e) decrypted Y1 with y1 ; (f) decrypted Y2 with y0 .

表4 密钥敏感性测试结果表Table 4.Key sensitivity test result table.

表5 明文图像与密文图像相关系数测试结果表Table 5.Plaintext image and ciphertext image correlation coefficient test result table.

从表5中可以看出, 明文图像在三个方向上的相关系数都大于0.6, 说明明文图像的相邻像素之间的相关性非常强; 而密文图像在三个方向上的相关系数都接近于0, 说明密文图像的相邻像素之间的相关性已被打破.

从图9可以看出, 明文图像的像素点集中分布在对角线附近, 说明明文图像的像素点之间的相关性很强; 而密文图像的像素点分布得比较均匀、散乱, 说明密文图像的像素之间的相关性已被破坏.因此提出的算法可以很好地降低图像相邻像素之间的强相关性, 对图像加密具有良好的效果.

图9 baboon 图像加密前后三个方向上的相关性分析图 (a), (b)对角相邻; (c), (d)水平相邻; (e), (f)垂直相邻;Fig.9.Correlation analysis chart in three directions before and after baboon image encryption: (a), (b) Diagonally adjacent; (c), (d)horizontally adjacent; (e), (f) vertically adjacent.

6.7 抗剪切能力分析

为测试该算法的抗剪切能力, 我们剪掉Lena加密图像中间60 × 60大小的图像, 如图10(b)所示, 再对剪切过的加密图像进行解密, 解密图像如图10(d)所示.图10(a)为原加密图像, 图10(c)为原加密图像的解密图像.对比图10(c)和图10(d)可以发现, 图10(d)中有些点的像素值发生了改变, 但仍然可以显示出明文图像的大致信息.因此,加密图像在遭受剪切攻击后仍然具有一定的解密效果.

6.8 抗噪声能力分析

图10 抗剪切攻击能力分析图 (a) 剪切前密文; (b)剪切后密文; (c)剪切前解密; (d) 剪切后解密Fig.10.Anti-shear attack capability analysis chart: (a) Ciphertext before cutting; (b) ciphertext after cutting; (c) decrypted image before cutting; (d) decrypted image after cutting.

图11 抗噪声攻击能力分析图 (a)加噪前密文; (b)加噪后密文; (c)加噪前解密; (d)加噪后解密Fig.11.Anti-noise attack capability analysis chart: (a) Ciphertext before adding noise; (b) ciphertext after adding noise; (c) decrypted image before adding noise; (d) decrypted image after adding noise.

椒盐噪声是图像中常见到的一种噪声, 当影像讯号受到突如其来的强烈干扰、类比数位转换器或位元传输错误等都可能产生椒盐噪声.为了测试本算法的抗椒盐噪声攻击能力, 我们对Lena加密图像加入1%的椒盐噪声, 加噪后的加密图像如图11(b)所示, 再对加噪后的加密图像进行解密,解密图像如图11(d)所示.图11(a)为原加密图像,图11(c)为原解密图像.对比图11(c)和图11(d)可以发现, 图11(d)中有些点的像素值发生了改变, 但仍然可以显示出原始明文图像的大致信息.这说明加密图像在遭受到椒盐噪声攻击后仍然具有一定的解密效果.

7 结 论

本文提出了一种新的五维多环多翼超混沌系统, 该系统结构简单.对该混沌系统的基本动力学特征进行了理论分析和数值仿真实验, 包括Lyapunov指数谱、平衡点、分岔图、时间序列、相图等, 从相图中可以明显看出, 该混沌系统能够在多方向上产生多环多翼, 从分岔图可以看出, 该系统参数的动态范围很广.同时, 将该五维多环多翼超混沌系统应用于物理混沌加密和代数加密的混合图像加密算法, 并对混合加密系统进行了数值仿真实验, 实验结果验证了该加密方法的正确性.因此, 本文提出的加密算法在数字图像的保密通信中有很好的应用前景.