CN107492064A - 基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法 - Google Patents

Abstract

本发明属于图像加密领域，特别涉及一种基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，包含：首先，图像经过离散小波变换，得到稀疏系数矩阵；然后采用zigzag置乱方法对稀疏系数矩阵进行置乱，再利用初等元胞自动机进行置乱操作；最后用忆阻混沌系统产生的测量矩阵对置乱后的图像进行压缩感知，得到最终的密文图像；其中，明文图像通过作用于SHA‑512函数来产生混沌系统的初始值、元胞自动机初始构型，增强算法与明文图像的相关性。本发明采用初等元胞自动机和压缩感知相结合的图像加密技术，利用初等元胞自动机对图像进行置乱，通过压缩感知在实现图像压缩的同时对图像进行加密，减少传输的数据量，并防止图像信息泄露，具有较高的安全性能。

Description

基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密 方法

技术领域

本发明涉及图像加密领域，特别涉及一种基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密 方法。

背景技术

近年来，随着互联网时代的到来，我们生活中的绝大多数信息都离不开网络的支持，我们利用网络进 行视频会议、发送一些机密信息等。图像作为信息的载体被广泛使用，而数字图像以其能够长期保存、便 于传输等优点，在网络上的传输、处理也变得越来越频繁。但是，网络也存在安全隐患，尤其是图像在传 输过程中面临着许多问题：图像信息被不法分子窃取、篡改，传输过程中受到噪声的影响等，这些因素使 现代多媒体数据的安全性已成为当务之急。因此，为了使图像安全的在网络上传输，需要利用有效的加密 技术处理图像，对图像信息进行保护。此外，需要传输的图像的数据量一般都比较大，信息冗余度高，为 了减少网络传输的数据量，图像需要先经过压缩再进入信道传输。目前，采用压缩感知并结合其它加密方 法对图像进行加密已经成为信息安全研究的一个热点，具有巨大的应用潜力以及较高的实用价值。压缩感 知理论指出：通过开发信号的稀疏特性，在远小于Nyquist采样率的条件下，用随机采样获取信号的离散 样本，然后通过非线性重建算法完美的重建信号。Candes和Donoho于2006年正式提出了压缩感知的概 念，之后，许多基于压缩感知的加密算法也被提出。将压缩感知运用到图像加密处理中，在保证数据安全 性的同时减少网络传输的数据量和存储空间，具有广阔的应用空间。

目前，混沌系统因其对初始条件和控制参数的高敏感性等特点，被广泛应用于图像加密处理中。混沌 系统可以分为两类：一维和高维的混沌系统。其中，一维混沌系统的参数和变量较少，相空间轨迹简单， 因此可以通过混沌信号估计技术来预测。相比之下，高维混沌系统的变量和参数相对较多，且具有更复杂 的结构和更好的混沌性能。因此，高维混沌系统是图像加密的潜在理想模型。近年来，普通混沌系统得到 了较为全面的发展，除了Chen系统、Lorenz系统等经典的混沌系统外，还有一些改进的、具有更优越非 线性特性的系统。其中，闵富红等人通过分析一种基于双曲正弦函数的新型磁控忆阻器模型的电压和电流 的相轨迹关系，发现其拥有典型的忆阻特性，故提出了一种新型忆阻混沌系统。该忆阻混沌系统具有更强 的非周期性，对初始值敏感性强，密钥空间大，且在物理上可以实现，将其应用于图像加密，可以设计出 一些安全的加密算法。元胞自动机(CA)是时间和空间都离散的动力系统，可以由简单的局部规则获得复 杂的全局行为，基于这一特点，当前，元胞自动机广泛应用于各个领域。在元胞自动机的发展过程中，科 学家们构造了各种各样的元胞自动机模型，各种CA模型在维数、可能状态的数量、邻域关系和状态更新 规则方面不同。其中，初等元胞自动机(ECA)是最简单的一维元胞自动机，它的状态集S只有两个元素 {s1,s2}，即状态个数k＝2，邻居半径r＝1。尽管它们构造简单，但可以导致复杂的行为并产生有用的 操作。目前，Abdel Latif Abu Dalhoum等提出了分别使用初等元胞自动机和二维元胞自动机置乱的图像加 密算法，达到了较好的置乱效果。但是所提出的算法密钥与明文相关性低，难以抵抗选择明文攻击，此外 加密后的图像冗余的数据并没有减少。

发明内容

为克服现有技术中的不足，本发明提出一种基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加 密方法，在保证信息安全性的前提下，通过采用初等元胞自动机和压缩感知相结合的图像加密技术，利用 初等元胞自动机进行置乱，通过压缩感知减少数据量的传输；使得密钥与明文密切相关，密钥空间较大， 能够有效抵抗不同攻击，安全性高，信息安全得到进一步保障。

按照本发明所提供的设计方案，一种基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法， 包含：明文图像利用SHA-512函数获取忆阻混沌系统的初始值和元胞自动机初始构型；明文图像通过离散 小波变换，得到稀疏系数矩阵；通过zigzag置乱方法对稀疏系数矩阵进行置乱，再利用初等元胞自动机进 行置乱操作；将忆阻混沌系统的初始值带入忆阻混沌系统中产生测量矩阵，通过测量矩阵对置乱后的矩阵 进行压缩感知，得到最终密文图像。

上述的，具体包含如下步骤：

步骤1、采用离散小波变换DWT，对大小为N×N的明文图像I进行变换，得到大小为N×N的稀疏 系数矩阵I₁；

步骤2、利用SHA-512函数对明文图像I进行计算，得到一组512位的哈希值，作为图像密钥Key， 将图像密钥Key换算为64个十进制数k₁,k₂,...,k₆₄，计算忆阻混沌系统的初始值x₀、y₀、z₀、w₀；

步骤3、从图像密钥Key中选取2个长度为N的0、1序列串，作为初等元胞自动机(ECA)的初始 构型其中为初始行构型，为初始列构型；

步骤4、通过zigzag置乱方法将稀疏系数矩阵I₁转换为大小为1×N²的1维向量p₁；

步骤5、选择初等元胞自动机的演变规则，并采取周期型边界、半径为1的方式对ECA的初始行构型列构型分别演变e次，得到每次演变的行构型列构型 其中，d＝1,2,…,e，演变总次数e为密钥的一部分；

步骤6、构造一个N×N的空白矩阵I₂，对p₁执行w轮置乱操作，得到大小为N×N的置乱矩阵I₄， 其中，w是置乱总轮数且其为密钥的一部分；

步骤7、将步骤2得到的初始值x₀、y₀、z₀、w₀带入忆阻混沌系统，得到四个大小为1×N的混沌 序列X、Y、Z、W，且X＝[x₁,x₂,...,x_N]，Y＝[y₁,y₂,...,y_N]，Z＝[z₁,z₂,...,z_N]，W＝[w₁,w₂,...,w_N]； 修正混沌序列X、Y、Z，得到序列X_1＝[x′₁,x′₂,...,x′_N]，Y_1＝[y′₁,y′₂,...,y′_N]，Z_1＝[z′₁,z′₂,...,z′_N]， 然后利用序列X_1、Y_1、Z_1，得到三个大小为1×N的序列U₁、U₂、U₃，表示为： U₁(i)＝[X_1(i)+Y_1(i)-Z_1(i)]，U₂(i)＝[X_1(i)+Z_1(i)-Y_1(i)]， U₃＝[Y_1(i)+Z_1(i)-X_1(i)]，其中，X_1(i)、Y_1(i)、Z_1(i)、U₁(i)、U₂(i)、U₃(i)分别 是序列X_1、Y_1、Z_1、U₁、U₂、U₃的第i个元素，i＝1,2,...,N；

步骤8、计算明文图像I的方差值var并修正，得到var1；根据var1的值从U₁、U₂、U₃中选取一 个序列并修正，作为序列U，大小为1×N；利用序列U构造大小为m×N的循环矩阵Φ，作为测量矩阵， 其中，m＝CR×N，CR是压缩率；

步骤9、利用测量矩阵Φ对经过ECA置乱后的矩阵I₄进行压缩感知采样，得到最终大小为m×N的 密文图像C，表示为：C＝ΦI₄。

上述的，步骤2包含如下内容：

步骤2.1、将512位的密钥Key每8位为一组，转换为64个十进制数k₁,k₂,...,k₆₄，然后利用公 式得到四个值h₁、h₂、h₃、h₄：

其中，sum(k₄₉,k₅₀,...,k₆₄)表示对k₄₉,k₅₀,…,k₆₄进行求和，max(k₄₉,k₅₀,...,k₆₄)表示求 k₄₉,k₅₀,...,k₆₄的最大值，代表取x和y的异或，t₂、t₃、t₄均为密钥的一部分；

步骤2.2、将h₁、h₂、h₃、h₄带入如下公式计算混沌系统的初始值x₀、y₀、z₀、w₀：

，其中，abs(x)表示求x的绝对值，mod(a,b)表示a对b的取模运算。

上述的，步骤3具体包含如下内容：a)从明文图像的512位哈希值中选取初始行构型选取规 则为：当N≤512时，从512位哈希值中按逆序的方向依次选取N个0、1值作为初等元胞自动机的初始 构型当512<N≤1024时，先按逆序的方向依次选取512位，再按正序的方向依次选取，直到选 够N位0、1序列，得到当N>1024时，先逆序选取、再正序、再逆序；按照该选取规则直到选 够长度为N的0、1序列作为b)从明文图像的512位哈希值中选取初始列构型规则如下：当 N≤512时，从512位哈希值中按正序的方向依次选取N个0、1值作为初等元胞自动机的初始构型当512<N≤1024时，先按正序的方向依次选取512位，再按逆序依次选取，直到选够N位0、1序列， 得到当N>1024时，先正序选取、再逆序、再正序；按照该规则直到选够长度为N的0、1序列作 为

上述的，步骤4包含如下内容：

步骤4.1、将由步骤2.1得到的h₁、h₂、h₃、h₄带入如下公式计算zigzag置乱时的初始位置(x′₀,y′₀)：

，式中，t₅是密钥的一部分，N是图像的尺寸；并对初始位置(x′₀,y′₀)进行修正，得到最终的初始位 置(x′₀,y′₀)，所采用的修正公式如下：

，式中，表示取比x小的最大整数，N是图像的尺寸；

步骤4.2、根据修正后得到的初始位置(x′₀,y′₀)，将稀疏系数矩阵I₁采用zigzag置乱方法转换为大小 为1×N²的1维向量p₁。

上述的，步骤6中的w轮置乱操作内容如下：在每一轮置乱中，根据的值，将向 量p₁中的元素填入空白矩阵I₂中，得到一轮置乱后的大小为N×N的矩阵I₃；接着，将矩阵I₃按列优先 转换为大小为1×N²的1维向量p₁，继续执行下一轮操作，共进行w轮置乱，得到最终结果，其中， i＝1,2,...,N，j＝1,2,...,N，d＝1,2,…,e，代表元胞自动机第d次演变的结果，e是演变总次 数且其为密钥的一部分。

优选的，步骤6具体包含如下步骤：

步骤6.1、令置乱轮数s＝1；

步骤6.2、构造一个N×N的空白矩阵I₂，令演变次数d＝1；

步骤6.3、将p₁中元素依次插入到矩阵I₂中坐标为(1,1)的位置上，以行优先，顺序插入，直到I₂中 所有坐标为的位置上填满元素，其中，i＝1,2,...,N，j＝1,2,...,N，代表元胞自动机第d次演变的结果；

步骤6.4、令演变次数d＝d+1，循环执行步骤6.3到步骤6.4，直至d＝e，其中，e是演变总次数且其为 密钥的一部分；

步骤6.5、经过e次演变后，倘若p₁中还有元素没有插入矩阵I₂中，则将剩余元素按照行优先依次插 入矩阵I₂中空白位置即可，完成一轮置乱，得到置乱矩阵I₃(N×N)；

步骤6.6、令置乱轮数s＝s+1，将矩阵I₃按列优先转换为1维向量p₁，循环执行步骤6.2到步骤6.6， 直至s＝w，得到最终的置乱矩阵I₄，其中，w是置乱总轮数且其为密钥的一部分。

更进一步地，步骤6.3中的以行优先，顺序插入，内容如下：每一次插入时，将p₁中剩下的元素依次 插入到矩阵I₂中当前坐标为(1,1)且为空的位置上，同样，按行优先插入；如果I₂中某坐标为(1,1)的 位置上有元素，则不插入，且继续寻找下一个坐标为(1,1)的位置，直至矩阵I₂中所有(1,1)位置上填 满元素。

上述的，步骤7具体包含如下步骤：

步骤7.1、将初始值x₀、y₀、z₀、w₀带入忆阻混沌系统，迭代n₀+N次；舍弃前n₀个值，得到四 个大小为1×N的混沌序列X、Y、Z、W，且X＝[x₁,x₂,...,x_N]，Y＝[y₁,y₂,...,y_N]，Z＝[z₁,z₂,...,z_N]，W＝[w₁,w₂,...,w_N]，其中，忆阻混沌系统是基于双曲正弦函数的新型磁控忆阻器模型所构造的新型忆阻 混沌系统，表达式为：

，式中，x，y，z，w代表混沌系统的状态变量，a，b，c，d，e，r 为混沌系统的控制参数且为实常数；

步骤7.2、对混沌序列X、Y、Z进行修正，得到序列X_1、Y_1、Z_1：X_1＝[x′₁,x′₂,...,x′_N]， Y_1＝[y′₁,y′₂,...,y′_N]，Z_1＝[z′₁,z′₂,...,z′_N]，其中，修正公式如下：

，式中，floor(x)表示取比x小的最大整数，x_i、y_i、z_i分别是序列X、Y、Z的第i个元素，X_1(i)、 Y_1(i)、Z_1(i)分别是序列X_1、Y_1、Z_1的第i个元素，i＝1,2,...,N；

步骤7.3、利用序列X_1、Y_1、Z_1，得到三个大小为1×N的序列U₁、U₂、U₃： U₁(i)＝[X_1(i)+Y_1(i)-Z_1(i)]，U₂(i)＝[X_1(i)+Z_1(i)-Y_1(i)]， U₃＝[Y_1(i)+Z_1(i)-X_1(i)]，其中，X_1(i)、Y_1(i)、Z_1(i)、U₁(i)、U₂(i)、U₃(i)分别 是序列X_1、Y_1、Z_1、U₁、U₂、U₃的第i个元素，i＝1,2,...,N。

上述的，步骤8具体包含如下步骤：

步骤8.1、按照如下公式计算明文图像I的方差值var并修正，得到var1：

var1＝floor(mod(var×10³,3))+1

其中，I(i,j)是明文图像I的灰度值，N是图像的尺寸，var1＝1,2,3；

步骤8.2、根据var1的值从三个序列U₁、U₂、U₃中选取一个序列，并对选取的序列进行修正，得 到大小为1×N的序列U，选取序列及修正规则具体如下：

如果var1＝1，那么U(i)＝mod(U₁(i),1)；

如果var1＝2，那么U(i)＝mod(U₂(i),1)；

如果var1＝3，那么U(i)＝mod(U₃(i),1)；

，其中，U(i)、U₁(i)、U₂(i)、U₃(i)分别是序列U、U₁、U₂、U₃的第i个元素，i＝1,2,...,N；

步骤8.3、计算测量矩阵大小，包含：测量矩阵大小设定为m×N，其中，m的值通过如下式子计算得 到：

m＝CR×N

，式中，CR是压缩率，N为图像的尺寸；

步骤8.4、利用序列U通过如下公式构造大小为m×N的测量矩阵Φ，其中，Φ(1,N)＝U：

Φ(j,1)＝λ·Φ(j-1,N)

Φ(j,2:N)＝Φ(j-1,1:N-1)

，式中，2≤j≤m，λ>1，且λ为密钥的一部分。

本发明的有益效果：

1、本发明对明文图像经过离散小波变换后得到的稀疏系数矩阵进行置乱，置乱分为两部分，先采用 zigzag置乱方法对图像进行置乱，再用初等元胞自动机进行置乱，此过程提高了置乱程度。然后，采用混 沌系统产生的测量矩阵对图像进行压缩感知，得到最终的密文图像。其中，混沌系统是一种基于双曲正弦 函数的新型磁控忆阻混沌系统，对密钥及明文都非常敏感，密钥空间较大。

2、本发明中的混沌系统的初始值、元胞自动机的初始构型等均与原始图像的哈希值有关，不同的明 文图像产生的密钥不同。为了进一步增强算法与明文像素的相关性，有效抵抗明文攻击，在构建测量矩阵 时，使用原始图像的方差值来选取混沌序列组，进而产生循环矩阵作为测量矩阵。

3、本发明采用初等元胞自动机和压缩感知相结合的图像加密技术，利用初等元胞自动机对明文图像 进行置乱，通过压缩感知在实现图像压缩的同时对图像进行加密，能够有效减少数据量的传输，并防止图 像信息泄露，同时具有良好的鲁棒性以及较高的安全性能。

附图说明：

图1为实施例一加密方法流程示意图；

图2为实施例二加密方法实现步骤示意图；

图3中：(a)为实施例二所用的元胞自动机初始构型选取顺序示意图，(b)为用元胞自动机表示矩阵 坐标的示意图；

图4为实施例二所用的元胞自动机初始构型产生方法的示意图；

图5中：(a)为实施例二中步骤6所用的矩阵A的示意图，(b)为zigzag置乱方法示意图，(c)为矩 阵A转换为1维向量P的示意图，(d)为元胞自动机170号演变规则示意图，(e)为初等元胞自动机的初 始构型R0、C0的示意图，(f)为周期型边界条件示意图；

图6中：(a)为实施例二中步骤6所用的空白矩阵B的示意图，(b)为元胞自动机初始构型及演变结 果示意图；

图7中：(a)为实施例二中步骤6所用的第一次演变之后矩阵B的示意图，(b)为第二次演变之后矩 阵B的示意图，(c)为第三次演变之后矩阵B的示意图，(d)为第一轮置乱结束后的矩阵B的示意图；

图8中：(a)为实施例三所用的512×512的Lena明文图像，(b)为压缩率为0.25的Lena加密图像， (c)为压缩率为0.25的解密图像，(d)为压缩率为0.5的Lena加密图像，(e)为压缩率为0.5的解密图 像，(f)为压缩率为0.75的Lena加密图像，(g)为压缩率为0.75的解密图像；

图9中：(a)为实施例三所用的使用错误密钥t₄＝2.670000000001时的解密图像，(b)为使用错误密 钥n₀＝801时的解密图像，(c)为使用错误密钥w＝2时的解密图像，(d)为使用错误密钥Key(1)＝6时 的解密图像；

图10中：(a)为实施例三所用的添加方差为0.0002的椒盐噪声后的加密图像，(b)为方差为0.0002 的椒盐噪声污染后的解密图像，(c)为添加方差为0.000001的高斯噪声后的加密图像，(d)为方差为0.000001 的高斯噪声污染后的解密图像；

图11中：(a)为实施例三所用的裁剪1/16后的加密图像，(b)为裁剪1/16后的解密图像，(c)为裁 剪1/32后的加密图像，(d)为裁剪1/32后的解密图像。

具体实施方式：

下面结合附图和技术方案对本发明作进一步详细的说明，并通过优选的实施例详细说明本发明的实施 方式，但本发明的实施方式并不限于此。

实施例一，针对现有图像加密过程中与明文相关性不高或者密钥空间不够大，参见图1所示，本实施 例提供一种基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，明文图像利用SHA-512函数 获取忆阻混沌系统的初始值和元胞自动机初始构型；明文图像通过离散小波变换，得到稀疏系数矩阵；通 过zigzag置乱方法对稀疏系数矩阵进行置乱，再利用初等元胞自动机进行置乱操作；将忆阻混沌系统的初 始值带入忆阻混沌系统产生测量矩阵，通过测量矩阵对置乱后的矩阵进行压缩感知，得到最终密文图像。

采用初等元胞自动机和压缩感知相结合的图像加密技术，利用初等元胞自动机进行置乱，通过压缩感 知减少数据量的传输；其中置乱分为两部分，先采用zigzag置乱方法对图像置乱，再用初等元胞自动机进 行置乱，此过程提高了置乱程度；最后采用忆阻混沌系统产生循环矩阵作为测量矩阵，对图像进行压缩感 知，得到最终的密文图像。其中，测量矩阵由一种基于双曲正弦函数的新型磁控忆阻混沌系统产生，且混 沌系统的初始值、元胞自动机初始构型都与原始图像的哈希值有关。该算法的密钥与明文密切相关，密钥 空间较大，置乱程度高，能够有效抵抗不同攻击，安全性高。

实施例二，参见图2所示，提供一种基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法， 其实现过程具体包含如下步骤：

步骤1、采用离散小波变换DWT，对大小为N×N的明文图像I进行变换，得到大小为N×N的稀疏 系数矩阵I₁。

步骤2、利用SHA-512函数对明文图像I进行计算，得到一组512位的哈希值并将它作为图像密钥Key， 然后将512位的图像密钥Key换算为64个十进制数k₁,k₂,...,k₆₄，计算忆阻混沌系统的初始值x₀、y₀、 z₀、w₀；具体通过如下步骤计算忆阻混沌系统的初始值：

步骤2.1、将512位的密钥Key每8位为一组，转换为64个十进制数k₁,k₂,...,k₆₄，然后利用公 式得到四个值h₁、h₂、h₃、h₄：

，其中，sum(k₄₉,k₅₀,...,k₆₄)表示对k₄₉,k₅₀,...,k₆₄进行求和，max(k₄₉,k₅₀,...,k₆₄)表示求 k₄₉,k₅₀,...,k₆₄的最大值，代表取x和y的异或，t₂、t₃、t₄是密钥的一部分。

步骤2.2、将h₁、h₂、h₃、h₄带入如下公式计算混沌系统的初始值x₀、y₀、z₀、w₀：

，其中，abs(x)表示求x的绝对值，mod(a,b)表示a对b的取模运算。

步骤3、从图像的512位哈希值中选取2个长度为N的0、1序列串，作为初等元胞自动机(ECA) 的初始构型具体选取方法及步骤如下：

步骤3.1、从图像的512位哈希值中选取的方法如下：当N≤512时，从512位哈希值中按图3 (a)中逆序的方向(即从后至前)依次选取N个0、1值作为初等元胞自动机的初始构型当 512<N≤1024时，先按逆序的方向依次选取512位，再按正序的方向(从前至后)依次选取，直到选 够N位0、1序列，得到当N>1024时，先逆序选取、再正序、再逆序，以此类推直到选够长度 为N的0、1序列作为如图4所示。

步骤3.2、从图像的512位哈希值中选取的方法如下：当N≤512时，从512位哈希值中按图3 (a)中正序的方向依次选取N个0、1值作为初等元胞自动机的初始构型当512<N≤1024时， 先按正序的方向依次选取512位，再按逆序的方向依次选取，直到选够N位0、1序列，得到当 N>1024时，先正序选取、再逆序、再正序，以此类推直到选够长度为N的0、1序列作为如图4 所示。

最终得到两个初始构型其中是初始行构型，是初始列构型，用的 每一次演变结果表示矩阵的行列坐标值，相当于一个N×N的矩阵，如图3(b)所示。

步骤4、采用zigzag置乱方法，将稀疏系数矩阵I₁(N×N)转换为1维向量p₁(1×N²)。其中zigzag 置乱时的初始位置(x′₀,y′₀)通过如下方法得到：

将由步骤2.1计算得到的h₁、h₂、h₃、h₄带入如下公式得到(x′₀,y′₀)：

式中，t₅是密钥的一部分，N是图像的尺寸。

接着，通过公式对(x′₀,y′₀)进行修正，得到最终的初始位置(x′₀,y′₀)：

，式中，表示取比x小的最大整数，N是图像的尺寸。

然后，根据初始位置(x′₀,y′₀)，将稀疏系数矩阵I₁(N×N)采用zigzag置乱方法转换为1维向量p₁(1 ×N²)。

步骤5、选择初等元胞自动机的某一规则，并采取周期型边界、半径为1的方式对ECA的初始行构型列构型分别演变e次，得到每次演变的行构型列构型 其中d＝1,2,…,e，演变总次数e是密钥的一部分。本实施例选取第170号规则， 如图5(d)所示。

初等元胞自动机的状态集S只有两个元素{s1,s2}，即状态个数k＝2，状态是离散时间步长中的元胞 的值，特别的，第一个状态称为初始状态(初始构型)。在状态集S中具体采用什么符号并不重要，它可 取{0，1}，{-l，1}等，重要的是S所含的符号个数，通常我们将其记为{0，1}。ECA的邻居半径r＝1， 即一个元胞由它的两个相邻元胞(左邻居和右邻居)的状态决定该元胞的下一个状态，如下：

state^r+1(i)＝f(state^r(i-1),state^r(i),state^r(i+1))

，其中，f(·)表示映射规则，state^r(i)表示第r次演变的第i个元胞的状态。规则定义了更新所有 元胞的同步状态的确定性方式。演变是根据规则将所有元胞的状态转移到下一个状态的过程。作为示例， 规则170(二进制：10101010)是state^r+1(i)＝0时 state^r(i-1),state^r(i),state^r(i+1)∈{110,100,010,000}和state^r+1(i)＝1时 state^r(i-1),state^r(i),state^r(i+1)∈{111,101,011,001}。周期型边界是指最左边(最右边)元胞的左(右)邻域是最右边(最左边)的元胞，如图5(f)所示。

步骤6、构造一个N×N的空白矩阵I₂，对p₁(1×N²)执行w轮置乱操作。在每一轮置乱中，根据 的值，将向量p₁中的元素填入空白矩阵I₂中，得到一轮置乱后的矩阵I₃(N×N)；接 着，将矩阵I₃按列优先转换为1维向量p₁(1×N²)，继续执行下一轮操作，共进行w轮置乱得到最终的 置乱矩阵I₄。其中i＝1,2,...,N，j＝1,2,...,N，d＝1,2,…,e，代表元胞自动机第d次演变的 结果，e是演变总次数、w是置乱总轮数且均为密钥的一部分。具体步骤如下：

步骤6.1、令置乱轮数s＝1。

步骤6.2、构造一个N×N的空白矩阵I₂，令演变次数d＝1。

步骤6.3、将p₁(1×N²)中元素依次插入到此时矩阵I₂中坐标为(1,1)(即) 的位置上，以行优先，顺序插入，直到此时I₂中所有坐标为的位置上填满元素。 每一次插入时将p₁中剩下的元素依次插入到矩阵I₂中当前坐标为(1,1)且为空的位置上，即矩阵I₂中坐 标为(1,1)且以前没有插入元素的位置，同样，按行优先插入。如果I₂中某坐标为(1,1)的位置上有元 素，则不插入，且继续寻找下一个坐标为(1,1)的位置，直至矩阵I₂中所有(1,1)位置上填满元素，其 中i＝1,2,...,N，j＝1,2,...,N，代表元胞自动机第d次演变的结果。

步骤6.4、令演变次数d＝d+1，循环执行步骤6.3到步骤6.4，共e次。其中，d＝1,2,…,e，e是演变总 次数且为密钥的一部分。

步骤6.5、经过e次演变后，倘若p₁中还有元素没有插入矩阵I₂中，则将剩余元素按照行优先依次插 入矩阵I₂中空白位置即可，此时，完成一轮置乱，得到置乱矩阵I₃(N×N)。

步骤6.6、令置乱轮数s＝s+1，将矩阵I₃按列优先转换为1维向量p₁(1×N²)，循环执行步骤6.2-6.6 共w次，即共进行w轮置乱，得到最终的置乱矩阵I₄。其中，s＝1,2,…,w，w是置乱总轮数且为密钥的一 部分。

为了便于描述本发明的置乱方法，举一个例子如下：

(1)假设存在一个4×4的矩阵A，如图5(a)。选取元素6为起点，利用图5(b)所示的zigzag置 乱方法，将其转换为1维向量P，结果如图5(c)所示。

(2)假设元胞自动机初始构型R0、C0如图5(e)所示，利用图5(d)所示170号规则演变3次， 得到Ri、Ci，其中i＝1,2,3，如图6(b)所示。在该例子中，R0、C0相当于步骤3中所述元胞自动机初 始行构型初始列构型Ri、Ci相当于步骤5中所述元胞自动机初始行列构型每一次演变的结果

(3)构造一个4×4的空白矩阵B，用序列Ri、Ci分别表示其行列坐标，其中i＝1,2,3，j＝1,2,3， 如图6(a)所示。

(4)利用步骤6所述置乱方法对一维向量P进行置乱。用序列R1代表矩阵B的行坐标，序列C1代 表矩阵B的列坐标。将P中元素，按行优先依次插入到矩阵B中坐标为(1,1)的位置，即 (R1(g),C1(h))＝(1,1)的位置上，得到图7(a)所示，其中g＝1,2,3,4，h＝1,2,3,4。

然后用R2、C2表示矩阵B的新的行列坐标值，将矩阵P中剩下元素依次插入到此时坐标为(1,1)且 为空的位置上，同样按行优先。这次演变之后，矩阵B中第三行、第二列坐标为(1,1)，但是该位置有元 素6，这是因为在上一次演变之后该位置坐标也为(1,1)且插入了元素6，此时跳过该位置不插入新元素， 继续寻找下一个坐标为(1,1)且为空的位置，直至矩阵B中所有(1,1)位置上填满元素，如图7(b)所 示。

接着用R3、C3表示矩阵B的新的行列坐标值，同上一步，将P中剩余元素依次插入到矩阵B中空白 位置。

最后，经过三次插入元素得到图7(c)所示的矩阵B，此时P中还有2个元素(5、9)，则将这2个 元素按行优先依次插入到矩阵B中空白位置，完成一轮置乱，得到此轮置乱后的矩阵B，如图7(d)所示。

步骤7、将步骤2得到的初始值x₀、y₀、z₀、w₀带入忆阻混沌系统，得到四个大小为1×N的混沌 序列X、Y、Z、W，且X＝[x₁,x₂,...,x_N]，Y＝[y₁,y₂,...,y_N]，Z＝[z₁,z₂,...,z_N]，W＝[w₁,w₂,...,w_N]。 修正混沌序列X、Y、Z，得到序列X_1＝[x′₁,x′₂,...,x′_N]，Y_1＝[y′₁,y′₂,...,y′_N]，Z_1＝[z′₁,z′₂,...,z′_N]， 然后利用序列X_1、Y_1、Z_1，得到三个大小为1×N的序列U₁、U₂、U₃： U₁(i)＝[X_1(i)+Y_1(i)-Z_1(i)]，U₂(i)＝[X_1(i)+Z_1(i)-Y_1(i)]， U₃＝[Y_1(i)+Z_1(i)-X_1(i)]，其中，X_1(i)、Y_1(i)、Z_1(i)、U₁(i)、U₂(i)、U₃(i)分别 是序列X_1、Y_1、Z_1、U₁、U₂、U₃的第i个元素，i＝1,2,...,N。具体步骤如下：

步骤7.1、将步骤2.2中计算得到的初始值x₀、y₀、z₀、w₀带入忆阻混沌系统，迭代n₀+N次(其 中n₀≥500且是密钥的一部分)。舍弃前n₀个值，得到四个大小为1×N的混沌序列X、Y、Z、W， 且X＝[x₁,x₂,...,x_N]，Y＝[y₁,y₂,...,y_N]，Z＝[z₁,z₂,...,z_N]，W＝[w₁,w₂,...,w_N]。

其中所涉及的混沌系统是一种基于双曲正弦函数的新型磁控忆阻器模型所构造的新型忆阻混沌系统， 它的表达式为：

，式中，x，y，z，w代表混沌系统的状态变量，a，b，c，d，e，r 为混沌系统的控制参数且为实常数。当a＝10，b＝8，c＝15，d＝5.2，e＝5且r＝1.5时，系统存在一个典 型的混沌吸引子，处于混沌状态。

步骤7.2、对混沌序列X、Y、Z进行修正，得到序列X_1、Y_1、Z_1：X_1＝[x′₁,x′₂,...,x′_N]， Y_1＝[y′₁,y′₂,...,y′_N]，Z_1＝[z′₁,z′₂,...,z′_N]，其中修正公式如下：

，式中，floor(x)表示取比x小的最大整数，x_i、y_i、z_i分别是序列X、Y、Z的第i个元素，X_1(i)、 Y_1(i)、Z_1(i)分别是序列X_1、Y_1、Z_1的第i个元素，i＝1,2,...,N。

步骤7.3、利用序列X_1、Y_1、Z_1，得到三个大小为1×N的序列U₁、U₂、U₃： U₁(i)＝[X_1(i)+Y_1(i)-Z_1(i)]，U₂(i)＝[X_1(i)+Z_1(i)-Y_1(i)]， U₃＝[Y_1(i)+Z_1(i)-X_1(i)]，其中，X_1(i)、Y_1(i)、Z_1(i)、U₁(i)、U₂(i)、U₃(i)分别 是序列X_1、Y_1、Z_1、U₁、U₂、U₃的第i个元素，i＝1,2,...,N。

步骤8、计算明文图像I的方差值var并修正，得到var1；根据var1的值从U₁、U₂、U₃中选取一 个序列并修正，作为序列U，大小为1×N；接着利用序列U构造大小为m×N的循环矩阵Φ作为测量矩 阵的具体步骤如下，其中m＝CR×N，CR是压缩率，var1为中间密钥。

步骤8.1、按照如下公式计算明文图像I的方差值var并修正，得到var1且var1＝1,2,3：

var1＝floor(mod(var×10³,3))+1

，其中，I(i,j)是明文图像I的灰度值，N是图像的尺寸。

步骤8.2、根据var1的值从三个序列U₁、U₂、U₃中选取一个序列，并对选取的序列进行修正，得 到大小为1×N的序列U。选取序列及修正规则具体如下：

如果var1＝1，那么U(i)＝mod(U₁(i),1)；

如果var1＝2，那么U(i)＝mod(U₂(i),1)；

如果var1＝3，那么U(i)＝mod(U₃(i),1)；

，其中，U(i)、U₁(i)、U₂(i)、U₃(i)分别是序列U、U₁、U₂、U₃的第i个元素，i＝1,2,...,N。

步骤8.3、计算测量矩阵大小。测量矩阵大小为m×N，其中m的值通过如下式子计算得到：

m＝CR×N

，式中，CR是压缩率。

步骤8.4、利用序列U通过如下公式构造大小为m×N的测量矩阵Φ，其中Φ(1,N)＝U：

Φ(j,1)＝λ·Φ(j-1,N)

Φ(j,2:N)＝Φ(j-1,1:N-1)

，式中，2≤j≤m，λ>1，且λ为密钥的一部分。

步骤9、利用测量矩阵Φ(m×N)对经过ECA置乱后的矩阵I₄(N×N)进行压缩感知采样：C＝ΦI₄， 得到最终大小为m×N的密文图像C。

图像经过离散小波变换，得到稀疏系数矩阵；然后采用zigzag置乱方法对稀疏系数矩阵进行置乱，再 利用初等元胞自动机进行置乱；最后用忆阻混沌系统产生的测量矩阵对置乱后的图像进行压缩感知，得到 最终的密文图像。其中，明文图像通过作用于SHA-512函数来产生混沌系统的初始值、元胞自动机初始构 型，增强了算法与明文图像的相关性。采用初等元胞自动机和压缩感知相结合的图像加密技术，利用初等 元胞自动机对图像进行置乱，通过压缩感知在实现图像压缩的同时对图像进行加密，减少了传输的数据量， 并防止图像信息泄露，具有较高的安全性能。

实施例三，参见图8～11所示，在本实施例中，采用的编程软件为Matlab R2016a，选取大小为512× 512的Lena灰度图像为实验对象，具体加密过程如下：

步骤1：输入原始大小为512×512的Lena灰度图像，用I＝imread('Lena.bmp')读取图像信息，采用离 散小波变换(DWT)，对图像I进行变换，得到大小为512×512的稀疏系数矩阵I₁。

步骤2：利用SHA-512函数对明文图像I进行计算，得到一组512位的哈希值并将它作为图像密钥Key， 然后将512位的图像密钥Key换算为64个十进制数k₁,k₂,...,k₆₄，计算忆阻混沌系统的初始值。具 体步骤如下：

2.1)利用SHA-512函数对明文图像I进行计算，得到一组512位的哈希值(十六进制表示为【7C 2 3 4 4A C 5 6 8AE F B D B 8 7D A 3 3D C E 5E D 4C C 8 6 6 4D F D 03 7 3 0B 5 3F 8 5 4 5B D 3 1B F 9 9 0F 1C E A 7 3 2B 5D 7 4 0F 3C E 2 3 0 1 33 1 6 0F 7C 8 4 5 2 4F F 5 4 2 1 5 9 4 9 5 7 4 9 4 5 1 2D C 2F C D 6 6A D 4 09 1】)，并把它作为图像的密钥Key。然后将其每8位为一组，转换为64个十进制数(124 3568 172 86 138 239 189 184 125 63 61 206 94 212 204 134 100 223 208 55 48 18163 133 69 189 49 191 153 15 28 234 115 43 93 116 15 60 226 48 19 49 96 247200 69 36 255 84 33 89 73 87 73 69 18 220 47 205 102 173 64 145)，并将其定义为k₁,k₂,...,k₆₄，接着利用如下公式得到四个值h₁、h₂、h₃、h₄：

，其中，sum(k₄₉,k₅₀,...,k₆₄)表示对k₄₉,k₅₀,...,k₆₄进行求和，max(k₄₉,k₅₀,...,k₆₄)表示求 k₄₉,k₅₀,...,k₆₄的最大值，代表取x和y的异或，t₂＝33.2418，t₃＝3.5609，t₄＝2.67是密钥的一 部分。

通过计算得到的四个值为h₁＝34.1559、h₂＝0.3477、h₃＝0.3984、h₄＝18.1874。

2.2)将h₁、h₂、h₃、h₄带入如下公式计算混沌系统的初始值x₀、y₀、z₀、w₀：

，其中，abs(x)表示求x的绝对值，mod(a,b)表示a对b的取模运算。通过计算得到x₀＝0.1559， y₀＝0.3477，z₀＝0.5982，w₀＝0.5858。

步骤3：从明文图像I的512位哈希值中选取2个长度为512的0、1序列串，作为初等元胞自动机(ECA) 的初始行构型初始列构型具体选取方法及步骤如下：

从明文图像I的512位哈希值中按逆序的方向(即从后至前)依次选取512个0、1值作为初等元胞自 动机的初始行构型然后，从512位哈希值中按正序的方向(即从前至后)依次选取512个0、1值 作为初等元胞自动机的初始列构型其中的每一次演变结果用作置乱时行坐标，的每一次演 变结果用作列坐标。

步骤4：采用zigzag置乱方法，将大小为512×512的稀疏系数矩阵I₁转换为大小为1×262144的1 维向量p₁，其中zigzag置乱时的初始位置由计算得到，具体步骤如下：

4.1)将由2_.1)得到的h₁、h₂、h₃、h₄带入如下公式计算zigzag置乱时的初始位置(x′₀,y′₀)：

，式中，令密钥t₅＝1.0314，N是图像的尺寸：N＝512。

并对初始位置(x′₀,y′₀)进行修正，得到最终的初始位置(x′₀,y′₀)，所采用的修正公式如下：

，式中，表示取比x小的最大整数，N是图像的尺寸：N＝512。计算得到的初始位置 (x′₀,y′₀)＝(356,184)。

4.2)根据初始位置(x′₀,y′₀)＝(356,184)，将大小为512×512的稀疏系数矩阵I₁采用zigzag置乱方法 转换为大小为1×262144的1维向量p₁。

步骤5：选择初等元胞自动机的170号规则，并采取周期型边界、半径为1的方式对ECA的初始行构 型列构型分别演变e次，在本实施例中选取e＝9，得到每次演变的行构型 列构型其中d＝1,2,…,9，演变总次数e是密钥 的一部分。

步骤6：构造一个512×512的空白矩阵I₂，对p₁(1×262144)执行w轮置乱操作。在每一轮置乱 中，根据的值，将向量p₁中的元素填入空白矩阵I₂中，得到一轮置乱后的矩阵I₃(512 ×512)；接着，将矩阵I₃按列优先转换为1维向量p₁(1×262144)，继续执行下一轮操作，共进行w轮 置乱得到最终的置乱矩阵I₄。其中i＝1,2,...,512，j＝1,2,...,512，d＝1,2,…,9，代表元胞 自动机第d次演变的结果，w是置乱总轮数且为密钥的一部分，在本实施例中选取w＝3。具体步骤如下：

6.1)令置乱轮数s＝1。

6.2)构造一个512×512的空白矩阵I₂，令演变次数d＝1。

6.3)将p₁(1×262144)中元素依次插入到此时矩阵I₂中坐标为(1,1)的位置上，即：以行优先，顺序插入，直到此时I₂中所有坐标为的 位置上填满元素。每一次插入时将p₁中剩下的元素依次插入到矩阵I₂中当前坐标为(1,1)且为空的位置 上，即矩阵I₂中坐标为(1,1)且以前没有插入元素的位置，同样，按行优先插入。如果I₂中某坐标为(1,1) 的位置上有元素，则不插入，且继续寻找下一个坐标为(1,1)的位置，直至矩阵I₂中所有(1,1)位置上 填满元素，其中i＝1,2,...,512，j＝1,2,...,512，代表元胞自动机第d次演变的结果。

6.4)令演变次数d＝d+1，循环执行步骤6.3到步骤6.4共e次，在本实施例中选取e＝9，即d＝1,2,…,9。

6.5)经过9次演变后，倘若p₁中还有元素没有插入矩阵I₂中，则将剩余元素按照行优先依次插入矩 阵I₂中空白位置即可，此时，完成一轮置乱，得到置乱矩阵I₃(512×512)。

6.6)令置乱轮数s＝s+1，将矩阵I₃(512×512)按列优先转换为1维向量p₁(1×262144)，循环执行 步骤6.2)-6.6)共w次，在本实施例中选取w＝3，即s＝1,2,3，得到最终的置乱矩阵I₄。

步骤7：将步骤2得到的初始值x₀、y₀、z₀、w₀带入忆阻混沌系统，得到四个大小为1×512的混 沌序列X、Y、Z、W，且X＝[x₁,x₂,...,x₅₁₂]，Y＝[y₁,y₂,...,y₅₁₂]，Z＝[z₁,z₂,...,z₅₁₂]， W＝[w₁,w₂,...,w₅₁₂]。修正混沌序列X、Y、Z，得到序列X_1＝[x′₁,x′₂,...,x′₅₁₂]，Y_1＝[y′₁,y′₂,...,y′₅₁₂]， Z_1＝[z′₁,z′₂,...,z′₅₁₂]，然后利用序列X_1、Y_1、Z_1，得到三个大小为1×512的序列U₁、U₂、U₃： U₁(i)＝[X_1(i)+Y_1(i)-Z_1(i)]，U₂(i)＝[X_1(i)+Z_1(i)-Y_1(i)]， U₃＝[Y_1(i)+Z_1(i)-X_1(i)]，其中，X_1(i)、Y_1(i)、Z_1(i)、U₁(i)、U₂(i)、U₃(i)分别 是序列X_1、Y_1、Z_1、U₁、U₂、U₃的第i个元素，i＝1,2,...,512。

7.1)将步骤2.2中计算得到的初始值x₀＝0.1559，y₀＝0.3477，z₀＝0.5982，w₀＝0.5858带入忆阻 混沌系统，迭代n₀+N次，在本实施例中选取n₀＝800，即800+512＝1312次。为了避免有害影响，舍弃 前800个值，得到四个大小为1×512的混沌序列X、Y、Z、W，且X＝[x₁,x₂,...,x₅₁₂]， Y＝[y₁,y₂,...,y₅₁₂]，Z＝[z₁,z₂,...,z₅₁₂]，W＝[w₁,w₂,...,w₅₁₂]。

其中所涉及的混沌系统是一种基于双曲正弦函数的新型磁控忆阻器模型所构造的新型忆阻混沌系统， 它的表达式为：

，式中，x，y，z，w代表混沌系统的状态变量，a，b，c，d，e，r 为混沌系统的控制参数且为实常数。当a＝10，b＝8，c＝15，d＝5.2，e＝5且r＝1.5时，系统存在一个典 型的混沌吸引子，处于混沌状态。

7.2)通过公式对混沌序列X、Y、Z进行修正，得到修正后的序列X_1、Y_1、Z_1： X_1＝[x′₁,x′₂,...,x′₅₁₂]，Y_1＝[y′₁,y′₂,...,y′₅₁₂]，Z_1＝[z′₁,z′₂,...,z′₅₁₂]：

，式中，floor(x)表示取比x小的最大整数，x_i、y_i、z_i分别是序列X、Y、Z的第i个元素，X_1(i)、 Y_1(i)、Z_1(i)分别是序列X_1、Y_1、Z_1的第i个元素，i＝1,2,...,512。

7.3)利用序列X_1、Y_1、Z_1，得到三个大小为1×512的序列U₁、U₂、U₃： U₁(i)＝[X_1(i)+Y_1(i)-Z_1(i)]，U₂(i)＝[X_1(i)+Z_1(i)-Y_1(i)]， U₃＝[Y_1(i)+Z_1(i)-X_1(i)]，其中，X_1(i)、Y_1(i)、Z_1(i)、U₁(i)、U₂(i)、U₃(i)分别 是序列X_1、Y_1、Z_1、U₁、U₂、U₃的第i个元素，i＝1,2,...,512。

其整个过程如下所示：

步骤8：计算明文图像I的方差值var并修正，得到var1，根据var1的值从U₁、U₂、U₃中选取一 个序列并修正，作为序列U，大小为1×512；接着利用序列U构造大小为256×512的循环矩阵Φ作为 测量矩阵。具体步骤如下：

8.1)首先，按照如下公式计算明文图像I的方差值var，在本实施例中，计算得到var＝2289：

式中，I(i,j)是明文图像I的灰度值，N是图像的尺寸：N＝512。

然后，通过如下公式修正var的值，得到var1：

var1＝floor(mod(var×10³,3))+1

本实施例中，计算得到var1＝1。

8.2)根据var1的值从U₁、U₂、U₃中选取一个序列并修正，作为序列U；选取序列及修正规则具 体如下：

如果var1＝1，那么U(i)＝mod(U₁(i),1)；

如果var1＝2，那么U(i)＝mod(U₂(i),1)；

如果var1＝3，那么U(i)＝mod(U₃(i),1)；

，其中U(i)、U₁(i)、U₂(i)、U₃(i)分别是序列U、U₁、U₂、U₃的第i个元素，i＝1,2,...,512。

在本实施例中，通过计算得到var1＝1，因此，选取序列U₁并修正，得到大小为1×512的序列U。

8.3)计算测量矩阵大小。测量矩阵大小为m×N，其中m的值通过如下式子计算得到：

m＝CR×N

，式中，CR是压缩率，N是图像的尺寸。本实施例中N＝512，选取CR＝0.5，则m＝0.5×512＝256， 故测量矩阵大小为256×512。

8.4)利用序列U通过如下公式构造大小为256×512的测量矩阵Φ，其中Φ(1,512)＝U：

Φ(j,1)＝λ×Φ(j-1,512)

Φ(j,2:512)＝Φ(j-1,1:511)

，式中，2≤j≤256，λ＝2.3且为密钥的一部分。

步骤9：利用测量矩阵Φ(256×512)对经过ECA置乱后的矩阵I₄(512×512)进行压缩感知采样： C＝ΦI₄，得到最终大小为256×512的密文图像C。

一个好的加密算法应该达到较好的加密效果和解密效果，且密钥空间足够大，密钥的敏感性足够高以 此来抵抗各种攻击。下面对本发明的图像加密算法进行安全性分析。

1、图像加密、解密效果好

好的加密算法不仅要能够达到较好的加密效果，即从密文中看不出明文图像的相关信息，同时接收方 根据相应的算法与密钥能恢复出较全的图像信息。为了评估加密质量和恢复效果，引入峰值信噪比PSNR 和平均结构相似性MSSIM。

峰值信噪比PSNR(Peak Signal to Noise Ratio)常常用来代表用户期望的重建图像的质量，PSNR的值 越大，原始图像和重建图像之间的相似度越高。

结构相似性SSIM(Structural Similarity)，是一种衡量两幅图像相似度的指标，分别从亮度、对比度、 结构三个方面度量图像相似性。SSIM取值范围为[0,1]，值越大，表示图像失真越小。在实际应用中，可 以利用滑动窗将图像分块，然后计算对应块的结构相似度SSIM，最后将平均值作为两图像的结构相似性 度量，即平均结构相似性MSSIM(MeanStructural Similarity)。

本文利用原始图像和重建图像之间的PSNR以及MSSIM来评估加密质量。首先，用如下公式定义PSNR：

，式中，X(i,j)和Y(i,j)分别代表明文图像X和解密图像Y的第(i,j)个像素值，MSE代表明文图像和 解密图像的均方差，N代表图像的尺寸，n为每像素的比特数，一般取8，即像素灰阶数为256。接着利用 如下公式定义SSIM和MSSIM：

，式中，μ_X、μ_Y分别表示明文图像X和解密图像Y的均值，σ_X、σ_Y分别表示明文图像和解密图 像的方差，σ_XY表示明文图像和解密图像的协方差，M表示图像分块的总数，在本实施例中，M＝64，C₁、 C₂、C₃为常数，为了避免分母为0的情况，取C₁＝(K₁×L)²，C₂＝(K₂×L)²，K₁＝0.01， K₂＝0.03，L＝255。

附图8是对512×512的Lena图像采用不同压缩率进行加密解密，得到的加密、解密结果。其中，图 8(a)是512×512的Lena明文图像，8(b)是压缩率CR＝0.25的Lena加密图像，8(c)是压缩率CR＝0.25 的解密图像，8(d)是压缩率CR＝0.5的Lena加密图像，8(e)是压缩率CR＝0.5的解密图像，8(f)是 压缩率CR＝0.75的Lena加密图像，8(g)是压缩率CR＝0.75的解密图像。

相应的，Lena图像在不同压缩率下的PSNR、MSSIM值如下：

可以看出，在不同的压缩率下，Lena图像可以被有效加密且从密文中得不到明文图像的相关信息； PSNR>30dB，MSSIM>0.99，原始图像和重构图像之间的相似度较高，图像失真小，因此该算法可以达到 很好的加密效果与恢复效果。

2、密钥空间足够大，足以抵抗穷举攻击

一般来说，密钥空间越大，算法抵抗穷举攻击的能力也就越强。一个安全的图像密码系统应该具有大 于2¹⁰⁰的密钥空间，以抵抗各种攻击。本文密钥具体包括：1)由SHA-512函数产生的512位哈希值；2) 舍弃的忆阻混沌序列的个数n₀；3)初等元胞自动机的演变次数e；4)置乱轮数w；5)求循环矩阵的参数 λ；6)计算初始值的四个参数t₂、t₃、t₄、t₅；7)中间密钥var1。

若设置精度为10^-14，则该算法的密钥空间约为(10¹⁴)⁵＝10⁷⁰>2²³²，如果再加上512位的哈希值，整 体密钥空间远远大于2¹⁰⁰，可见本发明的算法密钥空间足够大，能够有效的抵抗穷举攻击。

3、密钥敏感性高

密钥敏感性是密码学的基本特性之一，一个好的密码算法应当对密钥极其敏感。当攻击者用一个跟密 钥很相近的数据进行图像破解时，不能恢复出原始的图像，则说明该加密系统对密钥敏感。混沌系统对初 始条件和控制参数极其敏感，任何微小的初始偏差都会被指数式放大，因此混沌加密算法的安全性与密钥 敏感性有很大关系。本发明所使用的忆阻混沌系统，对初始值敏感性强，具有很强的抗破译功能。附图9 为密钥敏感性实验，其中忆阻混沌系统的控制参数是a＝10，b＝8，c＝15，d＝5.2，e＝5且r＝1.5，附图 8(e)是加密所选取的置乱轮数w＝3，演变次数e＝9，λ＝2.3，t₂＝33.2418，t₃＝3.5609，t₄＝2.67， t₅＝1.0314，n₀＝800，CR＝0.5，512位哈希值密钥Key为(十六进制表示为【7 C 2 3 4 4A C 56 8 A E F B D B 8 7 D A 3 3 D C E 5E D 4C C 8 6 6 4D F D 0 3 7 3 0B 5 3F 8 54 5B D 3 1B F 9 9 0F 1C E A 7 3 2 B 5D 7 4 0F 3C E 2 3 0 1 3 3 1 6 0F 7 C 8 45 2 4F F 5 4 2 1 5 9 4 9 5 7 4 9 4 5 1 2D C 2F C D 6 6A D 4 0 9 1】)时的正确解密图像。

附图9是本发明的加密算法密钥敏感性测试结果。其中，图9(a)是当其他密钥不变，参数 t₄＝2.670000000001时的解密图像；图9(b)是当其他密钥不变，参数n₀＝801时的解密图像；图9(c) 是当其他密钥不变，置乱轮数w＝2时的解密图像；图9(d)是当其他密钥不变，密钥Key(十六进制表 示)第一个值由7变为6时的解密图像。

4、可以抵抗一定的噪声、剪切攻击，鲁棒性强

(1)抗噪声攻击

图像在网络传输过程中或者在剪切、复制、移动的过程中，可能会受到机器自身的因素或者其它外界 因素而使获得的图像中含有噪声。本文通过产生不同程度的各类噪声然后添加到图8(d)所示的Lena密 文图像上，然后对其进行解密处理来测试本发明所提出的加密算法对噪声的抵抗能力。附图10是本发明 的抗噪声攻击能力测试结果，其中，附图10(a)是添加了方差为0.0002的椒盐噪声后的加密图像，图10 (b)是方差为0.0002的椒盐噪声污染后的解密图像，图10(c)是添加了方差为0.000001的高斯噪声后 的加密图像，图10(d)是方差为0.000001的高斯噪声污染后的解密图像。

从图10可以看出，即使密文图像受到一定的噪声攻击，也可以恢复出较清晰的图像，保留了明文图 像所包含的重要信息，这证明了本文所提出的算法具有一定的抗噪声攻击能力。

(2)抗剪切攻击

图像在处理或者传输过程中，由于设备的不健全或是出现误差等原因，有时我们只能获得部分加密图 像的信息。因此，研究部分信息丢失对解密图像的影响也十分重要。为了测试本发明所提出的加密算法的 抗剪切能力，对图8(d)所示的Lena加密图像进行部分裁剪，然后利用正确的密钥得出恢复的图像，测 试该算法对噪声的抵抗能力。附图11是本发明的抗剪切攻击能力测试结果。其中，图11(a)和(c)分 别显示了被遮挡1/16和1/32内容的加密图像，相应的解密结果如图11(b)和(d)所示。

从图11中可以看出，密文图像被裁剪掉一部分后，恢复的图像仍然保留了明文图像所包含的重要信 息，这证明了本文所提出的算法具有一定的抗剪切攻击能力。

通过测试可以看出，本算法可以抵抗一定的噪声攻击、剪切攻击，具有良好的鲁棒性，算法安全性高。

由上述实施例可知，本发明提供的加密算法可对N×N(N＝2^t，t为正整数)的灰度图像进行高安全的 加密，在信息加密领域具有广阔的应用前景。

本发明并不局限于上述具体实施方式，本领域技术人员还可据此做出多种变化，但任何与本发明等同 或者类似的变化都应涵盖在本发明权利要求的范围内。

Claims (10)

1.一种基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，其特征在于，包含：明文图像利用SHA-512函数获取忆阻混沌系统的初始值和元胞自动机初始构型；明文图像通过离散小波变换，得到稀疏系数矩阵；通过zigzag置乱方法对稀疏系数矩阵进行置乱，再利用初等元胞自动机进行置乱操作；将忆阻混沌系统的初始值带入忆阻混沌系统中产生测量矩阵，通过测量矩阵对置乱后的矩阵进行压缩感知，得到最终密文图像。

2.根据权利要求1所述的基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，其特征在于，具体包含如下步骤：

步骤1、采用离散小波变换DWT，对大小为N×N的明文图像I进行变换，得到大小为N×N的稀疏系数矩阵I₁；

步骤2、利用SHA-512函数对明文图像I进行计算，得到一组512位的哈希值，作为图像密钥Key，将图像密钥Key换算为64个十进制数k₁,k₂,...,k₆₄，计算忆阻混沌系统的初始值x₀、y₀、z₀、w₀；

步骤3、从图像密钥Key中选取2个长度为N的0、1序列串，作为初等元胞自动机ECA的初始构型其中为初始行构型，为初始列构型；

步骤4、通过zigzag置乱方法将稀疏系数矩阵I₁转换为大小为1×N²的1维向量p₁；

步骤5、选择初等元胞自动机的演变规则，并采取周期型边界、半径为1的方式对ECA的初始行构型列构型分别演变e次，得到每次演变的行构型列构型其中，d＝1,2,…,e，演变总次数e为密钥的一部分；

步骤6、构造一个N×N的空白矩阵I₂，对p₁执行w轮置乱操作，得到大小为N×N的置乱矩阵I₄，其中，w是置乱总轮数且其为密钥的一部分；

步骤7、将步骤2得到的初始值x₀、y₀、z₀、w₀带入忆阻混沌系统，得到四个大小为1×N的混沌序列X、Y、Z、W，且X＝[x₁,x₂,...,x_N]，Y＝[y₁,y₂,...,y_N]，Z＝[z₁,z₂,...,z_N]，W＝[w₁,w₂,...,w_N]；修正混沌序列X、Y、Z，得到序列X_1＝[x′₁,x′₂,...,x′_N]，Y_1＝[y′₁,y′₂,...,y′_N]，Z_1＝[z′₁,z′₂,...,z′_N]，然后利用序列X_1、Y_1、Z_1，得到三个大小为1×N的序列U₁、U₂、U₃，表示为：U₁(i)＝[X_1(i)+Y_1(i)-Z_1(i)]，U₂(i)＝[X_1(i)+Z_1(i)-Y_1(i)]，U₃＝[Y_1(i)+Z_1(i)-X_1(i)]，其中，X_1(i)、Y_1(i)、Z_1(i)、U₁(i)、U₂(i)、U₃(i)分别是序列X_1、Y_1、Z_1、U₁、U₂、U₃的第i个元素，i＝1,2,...,N；

步骤8、计算明文图像I的方差值var并修正，得到var1；根据var1的值从U₁、U₂、U₃中选取一个序列并修正，作为序列U，大小为1×N；利用序列U构造大小为m×N的循环矩阵Φ，作为测量矩阵，其中，m＝CR×N，CR是压缩率；

步骤9、利用测量矩阵Φ对经过ECA置乱后的矩阵I₄进行压缩感知采样，得到最终大小为m×N的密文图像C，表示为：C＝ΦI₄。

3.根据权利要求2所述的基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，其特征在于，步骤2包含如下内容：

步骤2.1、将512位的密钥Key每8位为一组，转换为64个十进制数k₁,k₂,...,k₆₄，然后利用公式得到四个值h₁、h₂、h₃、h₄：

其中，sum(k₄₉,k₅₀,...,k₆₄)表示对k₄₉,k₅₀,…,k₆₄进行求和，max(k₄₉,k₅₀,...,k₆₄)表示求k₄₉,k₅₀,...,k₆₄的最大值，代表取x和y的异或，t₂、t₃、t₄均为密钥的一部分；

步骤2.2、将h₁、h₂、h₃、h₄带入如下公式计算混沌系统的初始值x₀、y₀、z₀、w₀：

，

其中，abs(x)表示求x的绝对值，mod(a,b)表示a对b的取模运算。

4.根据权利要求2所述的基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，其特征在于：步骤3具体包含如下内容：a)从明文图像的512位哈希值中选取初始行构型选取规则为：当N≤512时，从512位哈希值中按逆序的方向依次选取N个0、1值作为初等元胞自动机的初始构型当512<N≤1024时，先按逆序的方向依次选取512位，再按正序的方向依次选取，直到选够N位0、1序列，得到当N>1024时，先逆序选取、再正序、再逆序；按照该选取规则直到选够长度为N的0、1序列作为b)从明文图像的512位哈希值中选取初始列构型规则如下：当N≤512时，从512位哈希值中按正序的方向依次选取N个0、1值作为初等元胞自动机的初始构型当512<N≤1024时，先按正序的方向依次选取512位，再按逆序依次选取，直到选够N位0、1序列，得到当N>1024时，先正序选取、再逆序、再正序；按照该规则直到选够长度为N的0、1序列作为

5.根据权利要求3所述的基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，其特征在于，步骤4包含如下内容：

步骤4.1、将由步骤2.1得到的h₁、h₂、h₃、h₄带入如下公式计算zigzag置乱时的初始位置(x′₀,y′₀)：

，

式中，t₅为密钥的一部分，N是图像的尺寸；并对初始位置(x′₀,y′₀)进行修正，得到最终的初始位置(x′₀,y′₀)，所采用的修正公式如下：

，

式中，表示取比x小的最大整数，N是图像的尺寸；

步骤4.2、根据修正后得到的初始位置(x′₀,y′₀)，将稀疏系数矩阵I₁采用zigzag置乱方法转换为大小为1×N²的1维向量p₁。

6.根据权利要求2所述的基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，其特征在于，步骤6中的w轮置乱操作内容如下：在每一轮置乱中，根据的值，将向量p₁中的元素填入空白矩阵I₂中，得到一轮置乱后的大小为N×N的矩阵I₃；接着，将矩阵I₃按列优先转换为大小为1×N²的1维向量p₁，继续执行下一轮操作，共进行w轮置乱，得到最终结果，其中，i＝1,2,...,N，j＝1,2,...,N，d＝1,2,…,e，代表元胞自动机第d次演变的结果，e是演变总次数且其为密钥的一部分。

7.根据权利要求2或6所述的基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，其特征在于，步骤6具体包含如下步骤：

步骤6.1、令置乱轮数s＝1；

步骤6.2、构造一个N×N的空白矩阵I₂，令演变次数d＝1；

步骤6.3、将p₁中元素依次插入到矩阵I₂中坐标为(1,1)的位置上，以行优先，顺序插入，直到I₂中所有坐标为的位置上填满元素，其中，i＝1,2,...,N，j＝1,2,...,N，代表元胞自动机第d次演变的结果；

步骤6.4、令演变次数d＝d+1，循环执行步骤6.3到步骤6.4，直至d＝e，其中，e是演变总次数且其为密钥的一部分；

步骤6.5、经过e次演变后，倘若p₁中还有元素没有插入矩阵I₂中，则将剩余元素按照行优先依次插入矩阵I₂中空白位置即可，完成一轮置乱，得到置乱矩阵I₃(N×N)；

步骤6.6、令置乱轮数s＝s+1，将矩阵I₃按列优先转换为1维向量p₁，循环执行步骤6.2到步骤6.6，直至s＝w，得到最终的置乱矩阵I₄，其中，w是置乱总轮数且其为密钥的一部分。

8.根据权利要求7所述的基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，其特征在于，步骤6.3中的以行优先，顺序插入，内容如下：每一次插入时，将p₁中剩下的元素依次插入到矩阵I₂中当前坐标为(1,1)且为空的位置上，同样，按行优先插入；如果I₂中某坐标为(1,1)的位置上有元素，则不插入，且继续寻找下一个坐标为(1,1)的位置，直至矩阵I₂中所有(1,1)位置上填满元素。

9.根据权利要求2所述的基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，其特征在于，步骤7具体包含如下步骤：

步骤7.1、将初始值x₀、y₀、z₀、w₀带入忆阻混沌系统，迭代n₀+N次；舍弃前n₀个值，得到四个大小为1×N的混沌序列X、Y、Z、W，且X＝[x₁,x₂,...,x_N]，Y＝[y₁,y₂,...,y_N]，Z＝[z₁,z₂,...,z_N]，W＝[w₁,w₂,...,w_N]，其中忆阻混沌系统是基于双曲正弦函数的新型磁控忆阻器模型所构造的新型忆阻混沌系统，表达式为：

，

式中，x，y，z，w代表混沌系统的状态变量，a，b，c，d，e，r为混沌系统的控制参数且为实常数；

步骤7.2、对混沌序列X、Y、Z进行修正，得到序列X_1、Y_1、Z_1：X_1＝[x′₁,x′₂,...,x′_N]，Y_1＝[y′₁,y′₂,...,y′_N]，Z_1＝[z′₁,z′₂,...,z′_N]，其中，修正公式如下：

，

式中，floor(x)表示取比x小的最大整数，x_i、y_i、z_i分别是序列X、Y、Z的第i个元素，X_1(i)、Y_1(i)、Z_1(i)分别是序列X_1、Y_1、Z_1的第i个元素，i＝1,2,...,N；

步骤7.3、利用序列X_1、Y_1、Z_1，得到三个大小为1×N的序列U₁、U₂、U₃：U₁(i)＝[X_1(i)+Y_1(i)-Z_1(i)]，U₂(i)＝[X_1(i)+Z_1(i)-Y_1(i)]，U₃＝[Y_1(i)+Z_1(i)-X_1(i)]，其中，X_1(i)、Y_1(i)、Z_1(i)、U₁(i)、U₂(i)、U₃(i)分别是序列X_1、Y_1、Z_1、U₁、U₂、U₃的第i个元素，i＝1,2,...,N。

10.根据权利要求2所述的基于忆阻混沌系统、初等元胞自动机和压缩感知的图像加密方法，其特征在于，步骤8具体包含如下步骤：

步骤8.1、按照如下公式计算明文图像I的方差值var并修正，得到var1：

var1＝floor(mod(var×10³,3))+1

其中，I(i,j)是明文图像I的灰度值，N是图像的尺寸，var1＝1,2,3；

步骤8.2、根据var1的值从三个序列U₁、U₂、U₃中选取一个序列，并对选取的序列进行修正，得到大小为1×N的序列U，选取序列及修正规则具体如下：

如果var1＝1，那么U(i)＝mod(U₁(i),1)；

如果var1＝2，那么U(i)＝mod(U₂(i),1)；

如果var1＝3，那么U(i)＝mod(U₃(i),1)；

其中，U(i)、U₁(i)、U₂(i)、U₃(i)分别是序列U、U₁、U₂、U₃的第i个元素，i＝1,2,...,N；

步骤8.3、计算测量矩阵大小，包含：测量矩阵大小设定为m×N，其中，m的值通过如下式子计算得到：

m＝CR×N ，式中，CR是压缩率，N为图像的尺寸；

步骤8.4、利用序列U通过如下公式构造大小为m×N的测量矩阵Φ，其中，Φ(1,N)＝U：

Φ(j,1)＝λ·Φ(j-1,N)

Φ(j,2:N)＝Φ(j-1,1:N-1) ，式中，2≤j≤m，λ>1，且λ为密钥的一部分。