基于新的五維多環多翼超混沌系統的圖像加密算法*

莊志本 李軍 劉靜漪 陳世強

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指數譜、平衡點、分岔圖、時間序列、相圖等, 從相圖中可以明顯看出, 該混沌系統能夠在多方向上產生多環多翼, 從分岔圖可以看出, 該系統參數的動態范圍很廣.同時, 將該五維多環多翼超混沌系統應用于物理混沌加密和代數加密的混合圖像加密算法, 并對混合加密系統進行了數值仿真實驗, 實驗結果驗證了該加密方法的正確性.因此, 本文提出的加密算法在數字圖像的保密通信中有很好的應用前景.