搜索

文章查询

x

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

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

庄志本 李军 刘静漪 陈世强

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

庄志本, 李军, 刘静漪, 陈世强
PDF
HTML
导出引用
导出核心图
  • 本文提出了一种基于新的五维多环多翼超混沌系统的数字图像加密方法. 首先, 将明文图像矩阵和五条混沌序列分别通过QR分解法分解成一个正交矩阵和一个上三角矩阵, 将混沌系统产生的五条混沌序列分别通过LU分解法分解成一个上三角矩阵和一个下三角矩阵, 分别将两个上三角矩阵和一个下三角矩阵相加, 得到五个离散后的混沌序列; 其次, 将明文图像矩阵分解出来的正交矩阵与五个混沌序列分解出来的五个正交矩阵相乘, 同时把明文图像矩阵分解出来的上三角矩阵中的元素通过混沌序列进行位置乱, 再将操作后的两个矩阵相乘; 最后, 将相乘后的矩阵通过混沌序列进行比特位位置乱, 再用混沌序列与其进行按位“异或”运算, 得到最终加密图像. 理论分析和仿真实验结果表明该算法的密钥空间远大于10200, 密钥敏感性强, 能够有效地抵御统计分析和灰度值分析的攻击, 对数字图像的加密具有很好的加密效果.
      通信作者: 陈世强, chensq8808@126.com
    • 基金项目: 其它-湖北民族学院博士启动基金项目(MY2018B014)
    [1]

    王平, 冯勇, 孙黎霞, 韩凤玲 2002 控制理论与应用 21 1

    Wang P, Feng Y, Sun L X, Han F L 2002 Control Theory & Appl. 21 1

    [2]

    禹思敏 2005 物理学报 54 1500

    Yu S M 2005 Acta Phys.Sin. 54 1500

    [3]

    Karthikeyan R, Serdar C, Peiman N, Abdul J M, Sajad J, Anitha K 2018 Eur. Phys. J. Plus 133 354

    [4]

    贾美美, 蒋浩刚, 李文静 2019 物理学报 68 130503

    Jia M M, Jiang H G, Li W J 2019 Acta Phys. Sin. 68 130503

    [5]

    Li Y X, Tang W K S, Chen G R 2005 Int. J. Bifurcation Chaos. 15 3367

    [6]

    彭再平, 王春华, 林愿, 骆小文 2014 物理学报 63 240506

    Peng Z P, Wang C H, Lin Y, Luo X W 2014 Acta Phys. Sin. 63 240506

    [7]

    刘杨 2015 博士学位论文 (哈尔滨: 哈尔滨工业大学)

    Liu Y 2015 Ph. D. Dissertation (Harbin: Harbin Institute of Technology) (in Chinese)

    [8]

    禹思敏 2018 新型混沌电路与系统的设计原理及其应用 (北京: 科学出版社) 第139—155页

    Yu S M 2018 Design Principles and Applications of New Chaotic Circuits and Systems (Beijing: Science Press) pp139– 155 (in Chinese)

    [9]

    Zhang L H, Liao X F, Wang X B 2005 Chaos, Solitons Fractals 24 759

    [10]

    Wong K, Kwor B, Law W 2008 Phys. Lett. A 372 2645

    [11]

    Zhang W, Yu H, Zhao Y L, Zhu Z L 2016 Signal Process. 118 36

    [12]

    Luo Y L, Zhou R L, Liu J X, Gao Y, Ding X M 2018 Nonlinear Dyn. 4 1

    [13]

    Ye G D, Pan C, Huang X L, Mei Q X 2018 Nonlinear Dyn. 20 18

    [14]

    Abanda Y, Tiedeu A 2016 IET Image Proc. 10 742

    [15]

    Zhang Y 2018 Inf. Sci. 255 31145

    [16]

    He Y, Zhang Y Q, Wang X Y 2018 Neural Comput. Appl. 10 1

    [17]

    Raza S F, Satpute V 2018 Nonlinear Dyn. 254 1

    [18]

    Ahmad J, Khan M A, Hwang S O, Khan J S 2017 Neural Comput. Appl. 28 953

    [19]

    Ahmad J, Khan M A, Ahmed F, Khan J S 2018 Neural Comput. Appl. 3 1

    [20]

    Sprott J C 1994 Phys. Rev. E 50 647

    [21]

    Enayatifar R, Abdullah A H, Isnin I F, Altameem A, Lee A 2017 Opt. Lasers Eng. 90 146

    [22]

    Liu H J, Wang X Y, Kadir A 2012 Appl. Soft Comput. 12 1457

  • 图 1  Lyapunov指数谱

    Fig. 1.  Lyapunov exponent diagram.

    图 2  系统(2)随$f$变化的分岔图

    Fig. 2.  Bifurcation diagram of system (2) variation with $f$.

    图 3  时间序列图 (a) $x\text- t$时间序列; (b) $y \text- t$时间序列; (c) $z \text- t$时间序列; (d) $w \text- t$时间序列; (e) $v \text- t$时间序列

    Fig. 3.  Time series diagram: (a) $x \text- t$ time series; (b) $y \text- t$ time series; (c) $z \text- t$ time series; (d) $w \text- t$ time series; (e) $v \text- t$time series

    图 4  三维相图 (a) $x \text- y \text- z$三维图; (b) $x \text- y \text- w$三维图; (c) $x \text- y \text- v$三维图; (d) $x \text- z \text- w$三维图; (e) $x \text- z \text- v$三维图; (f) $x \text- w \text- v$三维图; (g) $y \text- z \text- w$三维图; (h) $y \text- z \text- v$三维图; (i) $y \text- w \text- v$三维图; (j) $z \text- w \text- v$三维图

    Fig. 4.  Three-dimensional phase diagram: (a) $x \text- y \text- z$ Three-dimensional map; (b) $x \text- y \text- w$ Three-dimensional map; (c) $x \text- y \text- v$ Three-dimensional map; (d) $x \text- z \text- w$ Three-dimensional map; (e) $x \text- z \text- v$ Three-dimensional map; (f) $x \text- w \text- v$ Three-dimensional map; (g) $y \text- z \text- w$ Three-dimensional map; (h) $y \text- z \text- v$ Three-dimensional map; (i) $x \text- y \text- z$ Three-dimensional map; (j) $z \text- w \text- v$ Three-dimensional map.

    图 5  二维平面相图 (a) $x \text- y$平面; (b) $x \text- z$平面; (c) $x \text- v$平面; (d) $y \text- z$平面; (e) $y \text- w$平面; (f) $z \text- w$平面; (g) $z \text- v$平面

    Fig. 5.  Two-dimensional plane phase diagram: (a) $x \text- y$ flat; (b) $x \text- z$ flat; (c) $x \text- v$ flat; (d) $y \text- z$ flat; (e) $y \text- w$ flat; (f) $z \text- w$ flat; (g) $z \text- v$ flat.

    图 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.

    图 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.

    图 8  密钥敏感性测试图 (a)明文图像; (b)密文${{{Y}}_1}$(密钥为${y_0}$); (c)密文${{{Y}}_2}$(密钥为${y_1}$); (d) ${{{Y}}_1}$正确解密结果; (e) ${{{Y}}_1}$${y_1}$错误解密结果; (f) ${{{Y}}_2}$${y_0}$错误解密结果

    Fig. 8.  Key sensitivity tests: (a) Plain-image; (b) cipher ${{{Y}}_1}$ with key ${y_0}$; (c) cipher ${{{Y}}_2}$ with key ${y_1}$; (d) right decrypted ${{{Y}}_1}$; (e) decrypted ${{{Y}}_1}$ with ${y_1}$; (f) decrypted${{{Y}}_2}$ with${y_0}$.

    图 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.

    图 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.

    表 1  明文图像与加密图像的信息熵分析表

    Table 1.  Information entropy analysis table of plain text and encrypted image.

    图像Lena图像Baboon图像Boat图像
    原图像7.46447.37137.1267
    密文图像7.99947.99947.9993
    文献[11]7.99927.9993
    文献[12]7.99947.9993
    文献[13]7.99717.9993
    文献[14]7.9960
    文献[15]7.99937.9993
    文献[16]7.99937.9992
    文献[17]7.9974
    下载: 导出CSV

    表 2  加密图像不动点比分析表

    Table 2.  Encrypted image fixed point ratio analysis table.

    图像总像素数不动点数不动点比
    Lena图像26214410150.39%
    Baboon图像26214410140.39%
    Boat图像2621449990.38%
    下载: 导出CSV

    表 3  灰度平均变化值分析表

    Table 3.  Grayscale average change value analysis table.

    图像Lena图像Baboon图像Boat图像
    灰度平均变化值73.193770.858974.8383
    下载: 导出CSV

    表 4  密钥敏感性测试结果表

    Table 4.  Key sensitivity test result table.

    图像lena图像baboon图像boat图像
    指标NPCRUACINPCRUACINPCRUACI
    本文算法0.99640.33400.99620.33460.99580.3344
    文献[13]0.99620.33470.99600.3340
    文献[14]0.99570.3508
    文献[15]0.99610.33470.99610.3347
    文献[16]0.99620.33440.99610.3349
    下载: 导出CSV

    表 5  明文图像与密文图像相关系数测试结果表

    Table 5.  Plaintext image and ciphertext image correlation coefficient test result table.

    图像水平方向相关系数垂直方向相关系数对角线方向相关系数
    明文图像密文图像明文图像密文图像明文图像密文图像
    Lena0.9762–0.00840.96590.04610.94680.0131
    Baboon0.7204–0.00500.8264–0.00740.7046–0.0322
    boat0.96210.01060.82520.00870.8327–0.0423
    下载: 导出CSV
  • [1]

    王平, 冯勇, 孙黎霞, 韩凤玲 2002 控制理论与应用 21 1

    Wang P, Feng Y, Sun L X, Han F L 2002 Control Theory & Appl. 21 1

    [2]

    禹思敏 2005 物理学报 54 1500

    Yu S M 2005 Acta Phys.Sin. 54 1500

    [3]

    Karthikeyan R, Serdar C, Peiman N, Abdul J M, Sajad J, Anitha K 2018 Eur. Phys. J. Plus 133 354

    [4]

    贾美美, 蒋浩刚, 李文静 2019 物理学报 68 130503

    Jia M M, Jiang H G, Li W J 2019 Acta Phys. Sin. 68 130503

    [5]

    Li Y X, Tang W K S, Chen G R 2005 Int. J. Bifurcation Chaos. 15 3367

    [6]

    彭再平, 王春华, 林愿, 骆小文 2014 物理学报 63 240506

    Peng Z P, Wang C H, Lin Y, Luo X W 2014 Acta Phys. Sin. 63 240506

    [7]

    刘杨 2015 博士学位论文 (哈尔滨: 哈尔滨工业大学)

    Liu Y 2015 Ph. D. Dissertation (Harbin: Harbin Institute of Technology) (in Chinese)

    [8]

    禹思敏 2018 新型混沌电路与系统的设计原理及其应用 (北京: 科学出版社) 第139—155页

    Yu S M 2018 Design Principles and Applications of New Chaotic Circuits and Systems (Beijing: Science Press) pp139– 155 (in Chinese)

    [9]

    Zhang L H, Liao X F, Wang X B 2005 Chaos, Solitons Fractals 24 759

    [10]

    Wong K, Kwor B, Law W 2008 Phys. Lett. A 372 2645

    [11]

    Zhang W, Yu H, Zhao Y L, Zhu Z L 2016 Signal Process. 118 36

    [12]

    Luo Y L, Zhou R L, Liu J X, Gao Y, Ding X M 2018 Nonlinear Dyn. 4 1

    [13]

    Ye G D, Pan C, Huang X L, Mei Q X 2018 Nonlinear Dyn. 20 18

    [14]

    Abanda Y, Tiedeu A 2016 IET Image Proc. 10 742

    [15]

    Zhang Y 2018 Inf. Sci. 255 31145

    [16]

    He Y, Zhang Y Q, Wang X Y 2018 Neural Comput. Appl. 10 1

    [17]

    Raza S F, Satpute V 2018 Nonlinear Dyn. 254 1

    [18]

    Ahmad J, Khan M A, Hwang S O, Khan J S 2017 Neural Comput. Appl. 28 953

    [19]

    Ahmad J, Khan M A, Ahmed F, Khan J S 2018 Neural Comput. Appl. 3 1

    [20]

    Sprott J C 1994 Phys. Rev. E 50 647

    [21]

    Enayatifar R, Abdullah A H, Isnin I F, Altameem A, Lee A 2017 Opt. Lasers Eng. 90 146

    [22]

    Liu H J, Wang X Y, Kadir A 2012 Appl. Soft Comput. 12 1457

  • [1] 汪 敏, 胡小方, 伍小平. 物体内部三维位移场分析的数字图像相关方法. 物理学报, 2006, 55(10): 5135-5139. doi: 10.7498/aps.55.5135
    [2] 王兴元, 张继明. 一种基于混沌和汉明码的数字图像篡改检测及修复算法. 物理学报, 2014, 63(2): 020701. doi: 10.7498/aps.63.020701
    [3] 王兴元, 张继明. 一种基于抖动和混沌技术的数字图像篡改检测及修复算法. 物理学报, 2014, 63(21): 210701. doi: 10.7498/aps.63.210701
    [4] 林书庆, 江宁, 王超, 胡少华, 李桂兰, 薛琛鹏, 刘雨倩, 邱昆. 基于动态混沌映射的三维加密正交频分复用无源光网络. 物理学报, 2018, 67(2): 028401. doi: 10.7498/aps.67.20171246
    [5] 晋建秀, 丘水生. 基于物理混沌的混合图像加密系统研究. 物理学报, 2010, 59(2): 792-800. doi: 10.7498/aps.59.792
    [6] 彭再平, 王春华, 林愿, 骆小文. 一种新型的四维多翼超混沌吸引子及其在图像加密中的研究. 物理学报, 2014, 63(24): 240506. doi: 10.7498/aps.63.240506
    [7] 姚丽莉, 袁操今, 强俊杰, 冯少彤, 聂守平. 基于gyrator变换和矢量分解的非对称图像加密方法. 物理学报, 2016, 65(21): 214203. doi: 10.7498/aps.65.214203
    [8] 王静, 蒋国平. 一种超混沌图像加密算法的安全性分析及其改进. 物理学报, 2011, 60(6): 060503. doi: 10.7498/aps.60.060503
    [9] 朱从旭, 孙克辉. 对一类超混沌图像加密算法的密码分析与改进. 物理学报, 2012, 61(12): 120503. doi: 10.7498/aps.61.120503
    [10] 徐宁, 陈雪莲, 杨庚. 基于改进后多维数据加密系统的多图像光学加密算法的研究. 物理学报, 2013, 62(8): 084202. doi: 10.7498/aps.62.084202
    [11] 温贺平, 禹思敏, 吕金虎. 基于Hadoop大数据平台和无简并高维离散超混沌系统的加密算法. 物理学报, 2017, 66(23): 230503. doi: 10.7498/aps.66.230503
    [12] 石航, 王丽丹. 一种基于压缩感知和多维混沌系统的多过程图像加密方案. 物理学报, 2019, 68(20): 200501. doi: 10.7498/aps.68.20190553
    [13] 杨素丽, 符师桦, 蔡玉龙, 张迪, 张青川. 基于数字图像相关法的Mg含量对Al合金Protein-Le Chatelier效应影响的实验研究. 物理学报, 2017, 66(8): 086201. doi: 10.7498/aps.66.086201
    [14] 谢 鲲, 冯正进, 雷 敏. 一种超混沌系统的加密特性分析. 物理学报, 2005, 54(3): 1267-1272. doi: 10.7498/aps.54.1267
    [15] 匡锦瑜, 邓昆, 黄荣怀. 利用时空混沌同步进行数字加密通信. 物理学报, 2001, 50(10): 1856-1861. doi: 10.7498/aps.50.1856
    [16] 陈雪梅, 张静, 易兴文, 曾登科, 杨合明, 邱昆. 基于数字相干叠加的相干光正交频分复用系统中光纤非线性容忍性研究. 物理学报, 2015, 64(14): 144203. doi: 10.7498/aps.64.144203
    [17] 胡金秀, 高效伟. 变系数瞬态热传导问题边界元格式的特征正交分解降阶方法. 物理学报, 2016, 65(1): 014701. doi: 10.7498/aps.65.014701
    [18] 罗佳奇, 段焰辉, 夏振华. 基于自适应本征正交分解混合模型的跨音速流场分析. 物理学报, 2016, 65(12): 124702. doi: 10.7498/aps.65.124702
    [19] 段焰辉, 吴文华, 范召林, 罗佳奇. 基于本征正交分解的气动优化设计外形数据挖掘. 物理学报, 2017, 66(22): 220203. doi: 10.7498/aps.66.220203
    [20] 范虹, 韦文瑾, 朱艳春. 基于二维集合经验模式分解的距离正则化水平集磁共振图像分割. 物理学报, 2016, 65(16): 168701. doi: 10.7498/aps.65.168701
  • 引用本文:
    Citation:
计量
  • 文章访问数:  419
  • PDF下载量:  21
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-09-05
  • 修回日期:  2019-11-21
  • 刊出日期:  2020-02-01

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

  • 1. 湖北民族大学理学院, 恩施 445000
  • 2. 湖北民族大学信息工程学院, 恩施 445000
  • 3. 湖北民族大学新材料与机电工程学院, 恩施 445000
  • 通信作者: 陈世强, chensq8808@126.com
    基金项目: 其它-湖北民族学院博士启动基金项目(MY2018B014)

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

English Abstract

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

      目前, 构造新混沌系统的方法主要有三种, 一是通过引入一些非线性函数来产生多涡卷混沌吸引子[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 \in (0, 50)$时, 混沌系统都处于混沌状态; 将原文图像进行${\rm{QR}}$分解, 分别对正交矩阵${{Q}}$和上三角矩阵${{R}}$进行处理; 同时将混沌系统产生的5条混沌序列进行了${\rm{QR}}$分解和${\rm{LU}}$分解处理, 通过这些处理, 增大了穷举攻击的难度, 除穷举之外的一切基于明确的明文密文映射的攻击方法, 对该方法都将失效, 具有较高的安全性.

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

      i) ${{Q}}$的逆一定存在并且有${{{Q}}^{ - 1}} = {{{Q}}^{\rm{T}}}$;

      ii) ${{{Q}}^{ - 1}}{{Q}} = {{{Q}}^{\rm{T}}}{{Q}} = {{Q}}{{{Q}}^{ - 1}} = {{Q}}{{{Q}}^{\rm{T}}} = {{E}}$, (其中${{E}}$为单位矩阵);

      iii) 如果${{A}} = {{QP}}$, 那么${{P}} = {{{Q}}^{ - 1}}{{A}} = {{{Q}}^{\rm{T}}}A$;

      iv) 如果${{{Q}}_1}, {{{Q}}_2}, ..., {{{Q}}_n}$都是与${{Q}}$具有相同行列的正交矩阵, 且${{{A}}_1} = {{{Q}}_1}{{{Q}}_2}...{{{Q}}_n}{{P}}$, 则${{P}} = {{Q}}_n^{ - 1}{{Q}}_{n - 1}^{ - 1}...{{Q}}_1^{ - 1}{{{A}}_1}$.

    • 矩形矩阵的正交分解又称${\rm{QR}}$分解. ${\rm{QR}}$分解是把一个$m \times n$的矩阵${{A}}$分解成一个正交矩阵${{Q}}$和一个上三角矩阵${{R}}$. 即设矩阵

      ${{A}} = \left[ {\begin{array}{*{20}{c}} {{a_{11}}}&{{a_{12}}}&.&.&.&{{a_{1n - 1}}}&{{a_{1n}}} \\ {{a_{21}}}&{{a_{22}}}&.&.&.&{{a_{2n - 1}}}&{{a_{2n}}} \\ .&.&.&.&.&.&. \\ .&.&.&.&.&.&. \\ .&.&.&.&.&.&. \\ {{a_{m1}}}&{{a_{m2}}}&.&.&.&{{a_{mn - 1}}}&{{a_{mn}}} \end{array}} \right],$

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

      ${{Q}} = \left[ {\begin{array}{*{20}{c}} {{b_{11}}}&{{b_{12}}}&.&.&.&{{b_{1m}}} \\ {{b_{21}}}&{{b_{22}}}&.&.&.&{{b_{2m}}} \\ .&.&.&.&.&. \\ .&.&.&.&.&. \\ .&.&.&.&.&. \\ {{b_{m1}}}&{{b_{m2}}}&.&.&.&{{b_{mm}}} \end{array}} \right]$

      和一个上三角矩阵

      ${{R}} = \left[ {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}&{{c_{13}}}&.&.&.&{{c_{1n - 1}}}&{{c_{1n}}} \\ 0&{{c_{22}}}&{{c_{23}}}&.&.&.&{{c_{2n - 1}}}&{{c_{2n}}} \\ 0&0&{{c_{33}}}&.&.&.&{{c_{3n - 1}}}&{{c_{3n}}} \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ 0&0&0&.&.&.&0&{{c_{mn}}} \end{array}} \right].$

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

      $\begin{split} & {{A}} = {{Q}} \cdot {{R}} \\ =\; &\left[ {\begin{array}{*{20}{c}} {{a_{11}}}&{{a_{12}}}&.&.&.&{{a_{1n - 1}}}&{{a_{1n}}} \\ {{a_{21}}}&{{a_{22}}}&.&.&.&{{a_{2n - 1}}}&{{a_{2n}}} \\ .&.&.&.&.&.&. \\ .&.&.&.&.&.&. \\ .&.&.&.&.&.&. \\ {{a_{m1}}}&{{a_{m2}}}&.&.&.&{{a_{mn - 1}}}&{{a_{mn}}} \end{array}} \right]\\ =\; & \left[ {\begin{array}{*{20}{c}} {{b_{11}}}&{{b_{12}}}&.&.&.&{{b_{1m}}} \\ {{b_{21}}}&{{b_{22}}}&.&.&.&{{b_{2m}}} \\ .&.&.&.&.&. \\ .&.&.&.&.&. \\ .&.&.&.&.&. \\ {{b_{m1}}}&{{b_{m2}}}&.&.&.&{{b_{mm}}} \end{array}} \right] \\ & \times \left[ {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}&{{c_{13}}}&.&.&.&{{c_{1n - 1}}}&{{c_{1n}}} \\ 0&{{c_{22}}}&{{c_{23}}}&.&.&.&{{c_{2n - 1}}}&{{c_{2n}}} \\ 0&0&{{c_{33}}}&.&.&.&{{c_{3n - 1}}}&{{c_{3n}}} \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ 0&0&0&.&.&.&0&{{c_{mn}}} \end{array}} \right].\end{split}$

      $ $

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

      ${{X}} = \left[ {\begin{array}{*{20}{c}} {{a_{11}}}&{{a_{12}}}&.&.&.&{{a_{1n - 1}}}&{{a_{1n}}} \\ {{a_{21}}}&{{a_{22}}}&.&.&.&{{a_{2n - 1}}}&{{a_{2n}}} \\ .&.&.&.&.&.&. \\ .&.&.&.&.&.&. \\ .&.&.&.&.&.&. \\ {{a_{n1}}}&{{a_{n2}}}&.&.&.&{{a_{nn - 1}}}&{{a_{nn}}} \end{array}} \right],$

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

      ${{L}} = \left[ {\begin{array}{*{20}{c}} {{b_{11}}}&0&0&.&.&.&0&0 \\ {{b_{21}}}&{{b_{22}}}&0&.&.&.&0&0 \\ {{b_{31}}}&{{b_{32}}}&{{b_{33}}}&.&.&.&0&0 \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ {{b_{n1}}}&{{b_{n2}}}&{{b_{n3}}}&.&.&.&{{b_{nn - 1}}}&{{b_{nn}}} \end{array}} \right]$

      和一个上三角矩阵

      ${{U}} = \left[ {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}&{{c_{13}}}&.&.&.&{{c_{1n - 1}}}&{{c_{1n}}} \\ 0&{{c_{22}}}&{{c_{23}}}&.&.&.&{{c_{2n - 1}}}&{{c_{2n}}} \\ 0&0&{{c_{33}}}&.&.&.&{{c_{3n - 1}}}&{{c_{3n}}} \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ 0&0&0&.&.&.&0&{{c_{nn}}} \end{array}} \right].$

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

      $ \begin{split} &{{X}} = {{L}} \cdot {{U}} \\=\; & \left[ {\begin{array}{*{20}{c}} {{a_{11}}}&{{a_{12}}}&.&.&.&{{a_{1n - 1}}}&{{a_{1n}}} \\ {{a_{21}}}&{{a_{22}}}&.&.&.&{{a_{2n - 1}}}&{{a_{2n}}} \\ .&.&.&.&.&.&. \\ .&.&.&.&.&.&. \\ .&.&.&.&.&.&. \\ {{a_{n1}}}&{{a_{n2}}}&.&.&.&{{a_{nn - 1}}}&{{a_{nn}}} \end{array}} \right]\\ =\;& \left[ {\begin{array}{*{20}{c}} {{b_{11}}}&0&0&.&.&.&0&0 \\ {{b_{21}}}&{{b_{22}}}&0&.&.&.&0&0 \\ {{b_{31}}}&{{b_{32}}}&{{b_{33}}}&.&.&.&0&0 \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ {{b_{n1}}}&{{b_{n2}}}&{{b_{n3}}}&.&.&.&{{b_{nn - 1}}}&{{b_{nn}}} \end{array}} \right] \\ & \times \left[ {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}&{{c_{13}}}&.&.&.&{{c_{1n - 1}}}&{{c_{1n}}} \\ 0&{{c_{22}}}&{{c_{23}}}&.&.&.&{{c_{2n - 1}}}&{{c_{2n}}} \\ 0&0&{{c_{33}}}&.&.&.&{{c_{3n - 1}}}&{{c_{3n}}} \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ .&.&.&.&.&.&.&. \\ 0&0&0&.&.&.&0&{{c_{nn}}} \end{array}} \right].\end{split} $

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

      $\left\{ \begin{aligned} & \dot x = y, \\ & \dot y = - x + yz, \\ & \dot z = 1 - {y^2}. \end{aligned} \right.$

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

      $\left\{ \begin{aligned} & \dot x = ay, \\ & \dot y = - bx + cyz, \\ & \dot z = g - d{y^2}, \\ & \dot w = e{y^3}, \\ & \dot v = fz, \end{aligned} \right.$

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

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

      图  1  Lyapunov指数谱

      Figure 1.  Lyapunov exponent diagram.

    • $\dot x = \dot y = \dot z = \dot w = \dot v = 0$, 即

      $\left\{ \begin{aligned} & ay = 0,\\ & - bx + cyz = 0,\\ & g - d{y^2} = 0,\\ & e{y^3} = 0,\\ & fz = 0, \end{aligned} \right.$

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

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

      图  2  系统(2)随$f$变化的分岔图

      Figure 2.  Bifurcation diagram of system (2) variation with $f$.

    • 当系统参数$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  时间序列图 (a) $x\text- t$时间序列; (b) $y \text- t$时间序列; (c) $z \text- t$时间序列; (d) $w \text- t$时间序列; (e) $v \text- t$时间序列

      Figure 3.  Time series diagram: (a) $x \text- t$ time series; (b) $y \text- t$ time series; (c) $z \text- t$ time series; (d) $w \text- t$ time series; (e) $v \text- t$time series

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

      图  4  三维相图 (a) $x \text- y \text- z$三维图; (b) $x \text- y \text- w$三维图; (c) $x \text- y \text- v$三维图; (d) $x \text- z \text- w$三维图; (e) $x \text- z \text- v$三维图; (f) $x \text- w \text- v$三维图; (g) $y \text- z \text- w$三维图; (h) $y \text- z \text- v$三维图; (i) $y \text- w \text- v$三维图; (j) $z \text- w \text- v$三维图

      Figure 4.  Three-dimensional phase diagram: (a) $x \text- y \text- z$ Three-dimensional map; (b) $x \text- y \text- w$ Three-dimensional map; (c) $x \text- y \text- v$ Three-dimensional map; (d) $x \text- z \text- w$ Three-dimensional map; (e) $x \text- z \text- v$ Three-dimensional map; (f) $x \text- w \text- v$ Three-dimensional map; (g) $y \text- z \text- w$ Three-dimensional map; (h) $y \text- z \text- v$ Three-dimensional map; (i) $x \text- y \text- z$ Three-dimensional map; (j) $z \text- w \text- v$ Three-dimensional map.

      图  5  二维平面相图 (a) $x \text- y$平面; (b) $x \text- z$平面; (c) $x \text- v$平面; (d) $y \text- z$平面; (e) $y \text- w$平面; (f) $z \text- w$平面; (g) $z \text- v$平面

      Figure 5.  Two-dimensional plane phase diagram: (a) $x \text- y$ flat; (b) $x \text- z$ flat; (c) $x \text- v$ flat; (d) $y \text- z$ flat; (e) $y \text- w$ flat; (f) $z \text- w$ flat; (g) $z \text- v$ flat.

    • 给定一个$M \times N$的灰度图像${{P}}$, 加密步骤如下:

      Step1 输入灰度图像${{P}}$, 通过${\rm{QR}}$分解把${{P}}$分解成一个正交矩阵${{E}}$和一个$M \times N$的上三角矩阵${{R}}$;

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

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

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

      $\begin{split} & B = {\rm{reshape}}(y(201:(200 + M \times N),i),M,N)\\ & B1(:,:,i) = B,\quad {i = 1,2,3,4,5}\end{split}$

      其中$y(201:(200 + M \times N), i)$是把混沌序列$y(i)$的第201个值到第$200 + M \times N$ 个值取出来; ${\rm{reshape}}(y(201:(200 + M \times N), i), M, N)$是把取出来的$M \times N$个值转换成$M \times N$矩阵;

      Step5 第4步中的5个$M \times N$矩阵分别通过${\rm{QR}}$分解, 得到五个正交矩阵${{{{{E}}}}_1}, {{{E}}_2}, {{{E}}_3}, {{{E}}_4}, {{{E}}_5}$ 和五个上三角矩阵${{{R}}_1}, {{{R}}_2}, {{{R}}_3}, {{{R}}_4}, {{{R}}_5}$;

      Step6 改变初始值${y_0}$, 重复步骤二到步骤五$n(n \geqslant 2)$次, 得到$5 n$个正交矩阵和$5 n$个上三角矩阵. 分别记为${{{E}}_{11}}, {{{E}}_{12}}, {{{E}}_{13}}, {{{E}}_{14}}, {{{E}}_{15}}, {{{E}}_{21}}, \cdots, {{{E}}_{n1}}, $ ${{{E}}_{n2}}, {{{E}}_{n3}}, {{{E}}_{n4}},{{{E}}_{n5}}$;

      Step7 第一步中的正交矩阵${{E}}$与第六步中的正交矩阵${{{E}}_{11}}, {{{E}}_{12}}, {{{E}}_{13}}, {{{E}}_{14}}, {{{E}}_{15}}, {{{E}}_{21}}, \cdots, {{{E}}_{n1}},$$ {{{E}}_{n2}},{{{E}}_{n3}}, {{{E}}_{n4}}, {{{E}}_{n5}}$分别相乘, 得到矩阵${{E}}{{E}}$; 即${{E}}{{E}} = $$ {{{E}}_{11}} \times {{{E}}_{12}} \times {{{E}}_{13}}\, \times \,{{{E}}_{14}} \,\times\, {{{E}}_{15}} \,\times {{{E}}_{21}}\, \times \,\cdots \times {{{E}}_{n1}}\, \times$ ${{{E}}_{n2}} \times {{{E}}_{n3}} \times {{{E}}_{n4}} \times {{{E}}_{n5}} \times {{E}}$;

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

      $\begin{split} & B2 = {\rm{reshape}}(y(201:(200 + P \times P),i1),P,P)\\ & B3(:,:,i1) = B2,\quad {i1 = 1,2,3,4,5} , \end{split}$

      然后把这五个$P \times P$矩阵分别通过${\rm{LU}}$分解, 得到五个上三角矩阵${{{L}}_1}, {{{L}}_2}, {{{L}}_3}, {{{L}}_4}, {{{L}}_5}$ 和五个下三角矩阵${{{U}}_1}, {{{U}}_2}, {{{U}}_3}, {{{U}}_4}, {{{U}}_5}$;

      Step9 ${{{L}}_1}, {{{L}}_2}, {{{L}}_3}, {{{L}}_4}, {{{L}}_5}$${{{U}}_1}, {{{U}}_2}, {{{U}}_3}, {{{U}}_4}, {{{U}}_5}$分别对应相加, 即${{L}}{{{U}}_1} \!=\! {{{L}}_1} \!+\! {{{U}}_1}, {{L}}{{{U}}_2} \!=\! {{{L}}_2} \!+\! {{{U}}_2}, $${{L}}{{{U}}_3} = {{{L}}_3} + {{{U}}_3},\; {{L}}{{{U}}_4} = {{{L}}_4} + {{{U}}_4}, \;{{L}}{{{U}}_5} = {{{L}}_5} + {{{U}}_5}$. 然后分别取出${{L}}{{{U}}_i}(i$ $ = 1, 2, 3, 4, 5)$中的前$M$行与前$N$列位置上的元素, 从而得到五个$M \times N$矩阵. 分别记为${{M}}{{{N}}_1}, {{M}}{{{N}}_2}, {{M}}{{{N}}_3}, {{M}}{{{N}}_4}, {{M}}{{{N}}_5}$;

      Step10 ${{M}}{{{{{N}}}}_1}, {{M}}{{{{{N}}}}_2}, {{M}}{{{{{N}}}}_3}, {{M}}{{{{{N}}}}_4}, {{M}}{{{{{N}}}}_5}$${{{R}}_1},{{{R}}_2}, {{{R}}_3}, {{{R}}_4}, {{{R}}_5}$${{{R}}_2}, {{{R}}_3}, {{{R}}_4}, {{{R}}_5}$分别对应相加, 即${{{B}}_{41}} = {{M}}{{{{{N}}}}_1}\!+ {{{R}}_1}, $$ {{{B}}_{42}} \!=\! {{M}}{{{N}}_2} \!+\! {{{R}}_2}, {{{B}}_{43}} \!=\! {{M}}{{{N}}_3} $$+ {{{R}}_3}, {{{B}}_{44}} \!=\! {{M}}{{{N}}_4} \!+\! {{{R}}_4},{{{B}}_{45}} \!=\! {{M}}{{{N}}_5} \!+\! {{{R}}_5} $;

      Step11 ${{{B}}_{41}}, {{{B}}_{42}}, {{{B}}_{43}}, {{{B}}_{44}}, {{{B}}_{45}}$转换成1行$M \times N$列的矩阵, 同时把每个矩阵中的所有元素转换成0—255之间的整数, 分别得到矩阵${{{B}}_{51}}, {{{B}}_{52}}, {{{B}}_{53}}, {{{B}}_{54}}, {{{B}}_{55}}$;

      Step12 ${{{B}}_{51}}, {{{B}}_{52}}, {{{B}}_{53}}, {{{B}}_{54}}, {{{B}}_{55}}$作如下处理:

      $\begin{split} & {{{B}}_{61}}(i2) = {{{B}}_{51}}(i2) \oplus {{{B}}_{52}}(i2) \oplus {{{B}}_{53}}(i2) \oplus {{{B}}_{54}}(i2);\\ & {{{B}}_{62}}(i2) = {{{B}}_{51}}(i2) \oplus {{{B}}_{52}}(i2) \oplus {{{B}}_{53}}(i2) \oplus {{{B}}_{55}}(i2);\\ & {{{B}}_{63}}(i2) = {{{B}}_{51}}(i2) \oplus {{{B}}_{52}}(i2) \oplus {{{B}}_{54}}(i2) \oplus {{{B}}_{55}}(i2);\\ & {{{B}}_{64}}(i2) = {{{B}}_{51}}(i2) \oplus {{{B}}_{53}}(i2) \oplus {{{B}}_{54}}(i2) \oplus {{{B}}_{55}}(i2);\\ & {{{B}}_{65}}(i2) = {{{B}}_{52}}(i2) \oplus {{{B}}_{53}}(i2) \oplus {{{B}}_{54}}(i2) \oplus {{{B}}_{55}}(i2);\end{split}$

      其中$i2 = 1, 2, 3, \cdot \cdot \cdot, M \times N$;

      Step13 ${{{B}}_{61}}, {{{B}}_{62}}, {{{B}}_{63}}, {{{B}}_{64}}$作如下处理:

      $\begin{split} & {{{B}}_{71}} = {\rm{mod(round}}({{{B}}_{61}}(1:M) \times {10^5}),N - 1) + 1;\\ & {{{B}}_{72}} = {\rm{mod(round}}({{{B}}_{62}}(1:N) \times {10^5}),M - 1) + 1;\\ & {{{B}}_{73}} = {\rm{mod(round}}({{{B}}_{63}}(1:M \times N) \times {10^{10}}),7) + 1;\\ & {{{B}}_{74}} = {\rm{mod(round(}}{{{B}}_{64}}(({\rm{end}} - 7):{\rm{end}}) \times {10^{20}}),\\ & (M \times N - 1)) + 1;\\[-10pt]\end{split}$

      其中${\rm{mod}}$为求模运算, ${\rm{round}}(f)$是对$f$进行靠近取整;

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

      $\begin{split} & {{RR}}(i3,:) = {\rm{circshift}}({{R}}(i3,:),{{{B}}_{71}}(i3),2);\\ & i3 = 1,2,...,M,\end{split}$

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

      Step15 将${{RR}}$的所有列的元素进行扰乱, 扰乱公式如下:

      $\begin{split} &{{R}}{{{R}}_1}(:,i4) = {\rm{circshift}}({{RR}}(:,i4),{{{B}}_{72}}(i4),1);\\ & i4 = 1,2, \cdot \cdot \cdot,N,\end{split}$

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

      Step16 将${{EE}}$${{R}}{{{R}}_1}$相乘, 得到矩阵${{E}}{{{E}}_1}$, 即${{E}}{{{E}}_1} = {{EE}} \times {{R}}{{{R}}_1}$;

      Step17 将${{E}}{{{E}}_1}$转换成1行$M \times N$列的矩阵${{E}}{{{E}}_2}$, 再把${{E}}{{{E}}_2}$中的所有元素都映射成0—255之间的整数, 映射公式如下:

      $\begin{split} & {\rm{round}}\left(\frac{{{{E}}{{{E}}_2}(i5) - ({\rm{min}}({{E}}{{{E}}_2}) - 1)}}{{{\rm{max}}({{E}}{{{E}}_2}) - ({\rm{min}}({{E}}{{{E}}_2}) - 1)}} \times 255\right);\\ & i5 = 1,2, \cdot \cdot \cdot,M \times N;\\[-13pt]\end{split}$

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

      Step19 ${{E}}{{{E}}_2}$中的每一行二进制数分别进行扰乱, 得到矩阵${{E}}{{{E}}_3}$, 扰乱公式如下:

      $\begin{split} & {{E}}{{{E}}_3}(i6,:) = {\rm{circshif}}t({{E}}{{{E}}_2}(i6,:),{{{B}}_{73}}(i6),2);\\ & i6 = 1,2, \cdot \cdot \cdot,M \times N; \end{split}$

      Step20 ${{E}}{{{E}}_3}$中的每一列二进制数分别进行扰乱, 得到矩阵${{E}}{{{E}}_4}$, 扰乱公式如下:

      $\begin{split} & {{E}}{{{E}}_4}(:,i7) = {\rm{circshift}}({{E}}{{{E}}_3}(:,i7),{{{B}}_{74}}(i7),1);\\ & i7 = 1,2, \cdot \cdot \cdot,8; \end{split}$

      Step21 将${{E}}{{{E}}_4}$转换成十进制数, 再把已经转换成十进制数的${{E}}{{{E}}_4}$转换成1行$M \times N$列的矩阵${{E}}{{{E}}_5}$;

      Step22 ${{E}}{{{E}}_5}$${{{B}}_{65}}$进行按位“异或”运算, 得到一个新的序列${{E}}{{{E}}_6}$, 把${{E}}{{{E}}_6}$转换成$M$$ \times N$的矩阵${{E}}{{{E}}_7}$, ${{E}}{{{E}}_7}$为最终的加密图像.

    • Step1 输入图像矩阵${{E}}{{{E}}_7}$, 同时与${{{B}}_{65}}$进行按位“异或”运算, 得到一个新的序列${{E}}{{{E}}_8}$;

      Step2 ${{E}}{{{E}}_8}$转换成二进制数, 然后把行列进行还原, 得到矩阵${{E}}{{{E}}_9}$;

      Step3 ${{E}}{{{E}}_9}$转换成十进制数, 得到矩阵${{E}}{{{E}}_{10}}$, 再把${{E}}{{{E}}_{10}}$通过下列公式还原成${{E}}{{{E}}_1}$,

      $\begin{split} {{E}}{{{E}}_1}(j) =\; & ({{E}}{{{E}}_{10}}(j) \div 255) \times ({\rm{max}}({{E}}{{{E}}_2})\\ & - ({\rm{min}}({{E}}{{{E}}_2}) - 1)) \\ & + ({\rm{min}}({{E}}{{{E}}_2}) - 1),\end{split}$

      其中$j = 1, 2, \cdot \cdot \cdot, M \times N$;

      Step4 在${{E}}{{{E}}_1}$的左边依次乘上${{{E}}_{11}}, {{{E}}_{12}}, $$ {{{E}}_{13}}, {{{E}}_{14}}, {{{E}}_{15}}, {{{E}}_{21}}, \cdots, {{{E}}_{n1}}, {{{E}}_{n2}}, {{{E}}_{n3}}, {{{E}}_{n4}}, {{{E}}_{n5}}, E$的逆, 得到矩阵${{R}}{{{R}}_1}$;

      Step5 ${{R}}{{{R}}_1}$的行列进行还原, 得到矩阵${{R}}$, 再把矩阵${{E}}$${{R}}$相乘, 得到解密图像${{P}}$.

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

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

      图  6  数字图像加解密实验图 (a) lena原图; (b) lena加密图像; (c) lena解密图像; (d) baboon原图; (e) baboon加密图像; (f) baboon解密图像; (g) boat原图; (h) boat加密图像; (i) boat解密图像

      Figure 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.

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

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

      图  7  明文图像和密文图像直方图 (a) lena明文直方图; (b) lena密文直方图; (c) baboon明文直方图; (d) baboon密文直方图; (e) boat明文直方图; (f) boat密文直方图

      Figure 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.

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

      $H(m) = - \sum\limits_{i = 1}^L {p({m_i})} {\log _2}p({m_i}),$

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

      图像Lena图像Baboon图像Boat图像
      原图像7.46447.37137.1267
      密文图像7.99947.99947.9993
      文献[11]7.99927.9993
      文献[12]7.99947.9993
      文献[13]7.99717.9993
      文献[14]7.9960
      文献[15]7.99937.9993
      文献[16]7.99937.9992
      文献[17]7.9974

      表 1  明文图像与加密图像的信息熵分析表

      Table 1.  Information entropy analysis table of plain text and encrypted image.

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

      $\operatorname{BD} ({{G}},{{C}}) = \frac1{{MN}}{{\displaystyle\sum\limits_{i = 1}^M {\displaystyle\sum\limits_{j = 1}^N {f(i,j)} } }} \times 100\%,$

      其中$f(i, j) = \left\{ {aligned} & 1, \;\;{g_{ij}} = {c_{ij}} \\ & 0, \;\;{g_{ij}} \ne {c_{ij}}{aligned} \right.$. 由(15)式计算出本算法的不动点比如表2所列.

      图像总像素数不动点数不动点比
      Lena图像26214410150.39%
      Baboon图像26214410140.39%
      Boat图像2621449990.38%

      表 2  加密图像不动点比分析表

      Table 2.  Encrypted image fixed point ratio analysis table.

      ${\rm{GAVE}}({{C}},{{G}}) = \frac1{{MN}} {{\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {\left| {{c_{ij}} - } \right.\left. {{g_{ij}}} \right|} } }}, $

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

      图像Lena图像Baboon图像Boat图像
      灰度平均变化值73.193770.858974.8383

      表 3  灰度平均变化值分析表

      Table 3.  Grayscale average change value analysis table.

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

      图  8  密钥敏感性测试图 (a)明文图像; (b)密文${{{Y}}_1}$(密钥为${y_0}$); (c)密文${{{Y}}_2}$(密钥为${y_1}$); (d) ${{{Y}}_1}$正确解密结果; (e) ${{{Y}}_1}$${y_1}$错误解密结果; (f) ${{{Y}}_2}$${y_0}$错误解密结果

      Figure 8.  Key sensitivity tests: (a) Plain-image; (b) cipher ${{{Y}}_1}$ with key ${y_0}$; (c) cipher ${{{Y}}_2}$ with key ${y_1}$; (d) right decrypted ${{{Y}}_1}$; (e) decrypted ${{{Y}}_1}$ with ${y_1}$; (f) decrypted${{{Y}}_2}$ with${y_0}$.

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

      ${\rm{NPCR}} = \frac1{{M \times N}} {{\sum\limits_{i,j}^{M,N} {D\left( {i,j} \right)} }} \times 100\%,$

      ${\rm{UACI = }}\frac{{\displaystyle\sum\limits_{i,j}^{M,N} {\left| {{C_1}\left( {i,j} \right) - {C_2}\left( {i,j} \right)} \right|} }}{{M \times N \times 255}} \times 100\%,$

      其中$D\left( {i, j} \right) = \left\{ {aligned}& 1, \;\; {C_1}\left( {i, j} \right) \ne {C_2}\left( {i, j} \right) \\ & 0, \;\; {C_1}\left( {i, j} \right) = {C_2}\left( {i, j} \right) {aligned} \right.$; $M, N$为图像的大小, ${C_1}\left( {i, j} \right)$${C_2}\left( {i, j} \right)$分别是两幅图像在$\left( {i, j} \right)$位置上的像素值. NPCR和UACI的值越大, 说明两幅图像之间的差异就越大. 为了更好地评价算法的密钥敏感性, 我们测试了分别用密钥${y_0} \!=\! [1, 0.1, 2, 0.5, 0.4]$${y_1} \!= \![1, 0.1000000000000001,$ 2, 0.5, 0.4]加密的图像之间的NPCR和UACI的值, 测试值和文献[1316]的测试值如表4所列. 从表4可以看出, 测试结果良好, 因此提出的加密算法具有良好的密钥敏感性.

      图像lena图像baboon图像boat图像
      指标NPCRUACINPCRUACINPCRUACI
      本文算法0.99640.33400.99620.33460.99580.3344
      文献[13]0.99620.33470.99600.3340
      文献[14]0.99570.3508
      文献[15]0.99610.33470.99610.3347
      文献[16]0.99620.33440.99610.3349

      表 4  密钥敏感性测试结果表

      Table 4.  Key sensitivity test result table.

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

      图像水平方向相关系数垂直方向相关系数对角线方向相关系数
      明文图像密文图像明文图像密文图像明文图像密文图像
      Lena0.9762–0.00840.96590.04610.94680.0131
      Baboon0.7204–0.00500.8264–0.00740.7046–0.0322
      boat0.96210.01060.82520.00870.8327–0.0423

      表 5  明文图像与密文图像相关系数测试结果表

      Table 5.  Plaintext image and ciphertext image correlation coefficient test result table.

      图  9  baboon图像加密前后三个方向上的相关性分析图 (a), (b)对角相邻; (c), (d)水平相邻; (e), (f)垂直相邻;

      Figure 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.

      $E\left( u \right) = \frac{1}{N}\sum\limits_{i = 1}^N {{u_i}},$

      $D\left( u \right) = \frac{1}{N}\sum\limits_{i = 1}^N {{{({u_i} - E(u))}^2}},$

      ${\rm{Cov}}\left( {u,v} \right) = \frac{1}{N}\sum\limits_{i = 1}^N {({u_i} - E(u))({v_i} - E(v))},$

      ${r_{uv}} = \frac{{{\rm{Cov}}(u,v)}}{{\sqrt {D\left( u \right)} *\sqrt {D\left( v \right)} }}.$

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

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

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

      图  10  抗剪切攻击能力分析图 (a)剪切前密文; (b)剪切后密文; (c)剪切前解密; (d)剪切后解密

      Figure 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.

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

      图  11  抗噪声攻击能力分析图 (a)加噪前密文; (b)加噪后密文; (c)加噪前解密; (d)加噪后解密

      Figure 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.

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

参考文献 (22)

目录

    /

    返回文章
    返回