基于DCT离散余弦变换与混沌置乱的图像加密算法及MATLAB实现

做图像处理的人多少都会碰到这类任务把一张图片加密到视觉上完全不可读又能通过密钥无损还原。我最近把基于DCT离散余弦变换的图像加密算法完整跑了一遍从原理推导到MATLAB代码再到直方图、相关系数、信息熵这些安全性指标全部打通了。这里把整个思路、可复现代的代码以及踩过的坑整理出来给正在做数字图像处理课设、信息安全方向入门或者单纯想了解变换域加密原理的朋友参考。这套方案最大的优点是逻辑清晰、复现成本低用MATLAB自带函数就能实现。1. 为什么选DCT做图像加密思路与场景分析1.1 图像加密到底要解决什么问题图像数据和普通文本不一样它有很强的空间冗余性相邻像素灰度值往往高度相关。如果直接套用DES、AES这类传统分组密码把图像的二进制字节流当成普通数据加密虽然也能得到噪声图但你会发现三个问题一是运算量大图像动辄几百万像素实时性很难保证二是不符合图像本身的结构特点密码学上的扩散和混淆效果在图像上未必直观三是很多场景只需要保证视觉不可读并不需要做到和金融数据同等级别过重的加密反而拖慢速度。所以图像加密领域一般走的是另一条路利用图像自身的空间结构做文章。最常见的框架是“变换域处理 混沌序列密钥”也就是先把图像从像素域变换到频率域再在频率域对系数进行置乱或扰动最后反变换回空间域得到密文图。DCT离散余弦变换就是频率域工具里最经典、最成熟的一个这也是为什么很多课程设计和论文都选择DCT作为切入点。1.2 DCT变换的核心把像素空间搬到频率空间离散余弦变换的本质是把一个信号用不同频率的余弦函数叠加来表示。对图像这样的二维信号二维离散余弦变换公式可以写成F(u,v) (2/N)·c(u)·c(v)·ΣΣ f(x,y)·cos((2x1)uπ/(2N))·cos((2y1)vπ/(2N))其中 c(0)1/√2c(u0)1。这个式子看着复杂实际用MATLAB一行dct2()就搞定了。关键要理解它的物理含义变换后得到的系数矩阵里左上角是低频分量右下角是高频分量能量高度集中在低频区域。也就是说图像的大部分视觉内容被压缩进了少量低频系数中而高频系数数值很小主要对应边缘和纹理细节。这里有个很直观的类比DCT相当于把一幅图像“拆解”成不同频率成分的叠加就像把一段音乐拆成不同音高的音符。JPEG压缩之所以选用DCT而不是傅里叶变换是因为DCT是实数运算、没有复数开销而且对自然图像的能量集中效果更好。既然JPEG能靠DCT把图像压缩到十分之一大小还能保持视觉质量那反过来在DCT系数上做加密操作自然也就能以很小的代价彻底破坏图像的视觉内容。1.3 方案选型为什么是分块DCT加混沌置乱确定了在DCT域做文章之后下一个问题是直接整图DCT还是分块DCT。整图DCT的理论能量集中度最高但问题也很明显整图变换后低频系数全部堆在左上角一个区域如果对系数做置乱低频区域系数会扎堆移动加密分布不均匀而且整图DCT的运算量是O(N²logN)图像稍大就慢。分块DCT则不同把图像切成一个个8×8的小块分别做变换这是JPEG验证过的成熟方案块与块之间独立运算系数分布更均匀置乱起来也更彻底。所以我的方案选择分块DCT块大小默认取8。另一个关键决策是密钥的生成方式。加密系统的安全性很大程度上取决于密钥传统伪随机数发生器生成的序列周期有限密钥空间小容易受统计攻击。我选用Logistic混沌映射作为密钥序列来源它可以说是混沌加密里的“标配”x(n1) μ·x(n)·(1 - x(n))当控制参数μ处在3.57到4之间时系统处于混沌状态初值x0的微小差异会让迭代序列完全分道扬镳这正是密码学最看重的“初值敏感性”。我把μ和x0共同作为密钥把生成的混沌序列排序得到一组乱序索引用这个索引去置乱DCT系数矩阵。加密后的图像经过逆DCT回到空间域视觉上完全是噪声而解密时只要密钥有一丁点偏差置乱恢复就会错位全图彻底无法还原。2. 加密算法设计与MATLAB代码实现2.1 算法整体框架和参数设计整个加密流程可以拆成五个环节预处理读入图像转灰度图转double类型并做边界填充使宽高能被8整除。分块DCT对每个8×8小块分别做离散余弦变换得到系数矩阵。生成混沌序列用Logistic映射迭代长度等于系数总个数。系数置乱把混沌序列排序得到索引用索引重排所有DCT系数。逆变换输出把置乱后的系数矩阵分块进行逆DCT裁掉填充部分得到密文图像。解密是严格的反向流程先对密文图做分块DCT再用同样的混沌序列索引把系数恢复到原始位置最后逆DCT还原明文图像。这里有一个容易忽略但非常重要的设计点明文图像是0到255的整数DCT变换后系数变成了浮点数置乱后可以做逆变换再回到空间域加密到解密的整体链路里没有做任何量化所以理论上可逆。只是浮点运算会带来极其微小的误差解密后取round()就能还原回原始像素值。如果中间加了量化或者取整操作那就是有损加密解密图像会失真。参数方面核心参数如下表参数取值说明分块大小8×8JPEG标准块大小能量集中效果好Logistic控制参数μ3.993.57到4之间混沌区间越接近4混沌性越强混沌初值x00.3取值范围(0,1)作为密钥迭代长度等于像素总数每生成一个值对应一个系数位置2.2 三个核心模块的代码拆解模块一Logistic混沌序列生成function seq logistic_map(mu, x0, n) % 生成Logistic混沌序列 % mu: 控制参数建议取值3.57~4 % x0: 初值范围(0,1) % n: 序列长度 seq zeros(1, n); seq(1) x0; for k 2:n seq(k) mu * seq(k-1) * (1 - seq(k-1)); end end这个函数很朴素就是按公式迭代。但有几点实操经验第一初值x0不能取0、0.5、1这类不动点否则序列卡死在固定值上第二μ不要取到4.0的极限值数值计算时序列很容易溢出到负值后续排序就没法用第三序列的前几项还没有充分进入混沌状态实际使用时我通常从第50项开始取或者干脆先迭代100次再开始采样。模块二分块DCT变换分块DCT有两种写法。教学演示可以用两层for循环对每个块调用dct2这样逻辑直观dctCoef zeros(padM, padN); for i 1:block_size:padM for j 1:block_size:padN block imgPad(i:iblock_size-1, j:jblock_size-1); dctCoef(i:iblock_size-1, j:jblock_size-1) dct2(block); end end如果想写得更简洁推荐用blockproc函数它专门用于分块处理图像一行的效果等同上面两个循环fun (block_struct) dct2(block_struct.data); dctCoef blockproc(imgPad, [block_size block_size], fun);blockproc还支持并行和更复杂的块操作对性能有要求的场景更合适。我习惯把blockproc封装成一个函数加密解密都能用。模块三排序置乱与逆置乱置乱是整个加密的核心环节。思路是把混沌序列排序得到从小到大的索引顺序然后按这个索引重新排列DCT系数向量。关键点在于置乱和恢复的写法必须配套% 置乱方向cipherVec(j) plainVec(idx(j)) [~, idx] sort(seq); cipherVec plainVec(idx); % 恢复方向plainVec(idx(j)) cipherVec(j) recoverVec zeros(size(cipherVec)); recoverVec(idx) cipherVec;这个写法经常有人写反。解密时recoverVec(idx) cipherVec的意思是“把密文向量里的第j个值放回到原始位置idx(j)上去”因为加密时第idx(j)个原始值被移动到了第j个位置。想清楚这一步解密代码就不会错。2.3 完整加密解密主函数下面是完整可运行的加密函数我加了比较详细的注释function cipher dct_image_encrypt(plain, mu, x0, block_size) % 基于DCT与Logistic混沌置乱的图像加密 % 输入 % plain - 灰度图像矩阵uint8或double范围0~255 % mu - Logistic控制参数建议3.99 % x0 - 混沌初值建议0.3 % block_size - 分块大小默认8 % 输出 % cipher - 密文图像矩阵double类型 if nargin 4, block_size 8; end plain im2double(plain); % 注意这里归一化到0~1 [M, N] size(plain); % 1. 填充至分块大小的整数倍 padM ceil(M / block_size) * block_size; padN ceil(N / block_size) * block_size; imgPad zeros(padM, padN); imgPad(1:M, 1:N) plain; % 2. 分块DCT fun (block_struct) dct2(block_struct.data); dctCoef blockproc(imgPad, [block_size block_size], fun); % 3. 生成混沌序列并排序得到置乱索引 total padM * padN; seq logistic_map(mu, x0, total 100); seq seq(101:end); % 跳过前100次迭代保证充分混沌 [~, idx] sort(seq); % 4. 对DCT系数做全局置乱 vec dctCoef(:); vecShuffled vec(idx); dctShuffled reshape(vecShuffled, padM, padN); % 5. 分块逆DCT得到密文 funInv (block_struct) idct2(block_struct.data); cipher blockproc(dctShuffled, [block_size block_size], funInv); cipher cipher(1:M, 1:N); % 裁掉填充区 end解密函数是逆向流程代码几乎镜像function recover dct_image_decrypt(cipher, mu, x0, block_size) % 基于DCT与Logistic混沌置乱的图像解密 % 输入输出说明同加密函数 if nargin 4, block_size 8; end [M, N] size(cipher); padM ceil(M / block_size) * block_size; padN ceil(N / block_size) * block_size; imgPad zeros(padM, padN); imgPad(1:M, 1:N) cipher; % 1. 分块DCT回到系数域 fun (block_struct) dct2(block_struct.data); dctCoef blockproc(imgPad, [block_size block_size], fun); % 2. 用相同密钥生成相同索引 total padM * padN; seq logistic_map(mu, x0, total 100); seq seq(101:end); [~, idx] sort(seq); % 3. 逆置乱恢复原始系数 vec dctCoef(:); recoverVec zeros(size(vec)); recoverVec(idx) vec; dctRecover reshape(recoverVec, padM, padN); % 4. 逆DCT回到像素域 funInv (block_struct) idct2(block_struct.data); recover blockproc(dctRecover, [block_size block_size], funInv); recover recover(1:M, 1:N); end测试脚本可以这样写img imread(cameraman.tif); mu 3.99; x0 0.3; cipher dct_image_encrypt(img, mu, x0, 8); recover dct_image_decrypt(cipher, mu, x0, 8); figure; subplot(1,3,1); imshow(img); title(原始图像); subplot(1,3,2); imshow(cipher, []); title(密文图像); subplot(1,3,3); imshow(recover, []); title(解密图像);这里有个显示细节我必须强调密文图像矩阵是double类型而且DCT逆变换后可能有负值直接imshow(cipher)会因为数据范围不是0~255而显示异常甚至全黑全白。要加上[]参数让MATLAB自动把最小值映射到黑色、最大值映射到白色也就是imshow(cipher, [])。解密图像同理因为浮点误差的存在直接用imshow显示可能发灰偏暗建议先对解密图像做round()再显示或者也加上[]。3. 实验效果与安全性指标验证3.1 主观效果加密后还能不能看跑完上面的测试代码视觉效果很明显原始cameraman图像是一张清晰的灰度照片密文图像则完全是雪花噪声看不出任何结构轮廓。这里面的原理在于DCT系数矩阵里的数值代表的是不同频率成分当我把低频系数和高频系数的位置完全打乱后图像的主要能量信息被分散到各个地方反变换回空间域后原本规则的像素排列被彻底破坏。解密图像如果密钥正确视觉上和原始图像几乎无差别。由于浮点运算误差解密结果一般是原始值的±0.01以内显示出来肉眼看不出任何区别。这里再强调一次如果要做像素级别的比对判定是否无损应该对解密图像做im2uint8(round(recover))后再比较通常能完全还原。3.2 客观指标直方图、相关系数与信息熵视觉效果之外学术界评估图像加密算法一般看三个客观指标像素直方图、相邻像素相关系数、信息熵。先看直方图。原始图像的直方图通常有非常明显的波峰波谷反映像素灰度分布的不均匀性。好的加密算法会让密文图像像素值均匀分布直方图趋近一根水平线。用MATLAB看直方图很简单figure; subplot(1,2,1); imhist(uint8(255*im2double(img))); title(原始图像直方图); subplot(1,2,2); imhist(uint8(255*mat2gray(cipher))); title(密文图像直方图);注意密文要先归一化到0~255再转uint8计算直方图。再看相邻像素相关系数。自然图像相邻像素之间相关性极强相关系数通常接近1。加密后这个值应该跌到接近0说明相邻像素已经没有统计关系。计算代码function r corr_coef(img, direction) % 计算图像相邻像素相关系数 % direction: h水平, v垂直, d对角 img double(img); [M, N] size(img); switch direction case h x img(1:M, 1:N-1); y img(1:M, 2:N); case v x img(1:M-1, 1:N); y img(2:M, 1:N); case d x img(1:M-1, 1:N-1); y img(2:M, 2:N); end x x(:); y y(:); r sum((x - mean(x)) .* (y - mean(y))) / ... sqrt(sum((x - mean(x)).^2) * sum((y - mean(y)).^2)); end实测典型结果原始图像水平方向相关系数在0.95以上密文图像三个方向的相关系数基本都落在-0.05到0.05之间说明置乱彻底切断了像素间关系。信息熵衡量数据的不确定性和随机性8位图像的理想最大熵是8。熵越大说明像素值分布越随机。计算公式和代码p imhist(uint8(img)) / numel(img); p p(p ~ 0); H -sum(p .* log2(p));原始图像信息熵通常在7.2到7.5之间密文图像信息熵能到7.99以上非常接近理论最大值说明密文像素取值几乎是等概率的统计攻击很难找到规律。3.3 抗差分攻击与密钥敏感性测试差分攻击是图像加密领域最常见的攻击方式攻击者不断修改明文图像的单个像素观察密文变化试图反推密钥。对应指标是NPCR和UACI。NPCR衡量两幅图像不同像素的比例UACI衡量像素值差异的平均程度。对一幅256级灰度图像两个随机噪声图像的NPCR理论上接近99.6%UACI理论上接近33.4%。实现代码function [npcr, uaci] calc_npcr_uaci(img1, img2) img1 double(img1); img2 double(img2); D img1 ~ img2; npcr sum(D(:)) / numel(D) * 100; uaci sum(abs(img1(:) - img2(:))) / (255 * numel(D)) * 100; end测试方法是取一张图像加密修改其中一个像素比如把像素(50,50)从原值改成原值1再加密计算两个密文的NPCR和UACI。实测这个方案的NPCR普遍在99.6%以上UACI在33%上下说明单个像素的变化会被DCT变换和全局系数置乱扩散到整幅图像抗差分攻击能力是合格的。密钥敏感性测试更加直观用x00.3加密再用x00.3000000001解密肉眼上看解密结果仍然是雪花噪声PSNR只有8dB左右完全无法辨认。这说明方案的密钥空间足够大初值在10的负10次方级别的差异就足以让混沌序列完全改变置乱索引瞬间失效。这也是混沌系统做图像加密的核心优势。4. 实操中的常见问题与避坑记录4.1 运行环境和工具箱依赖先说环境。整套代码依赖MATLAB的Image Processing Toolbox因为dct2、idct2、blockproc、im2double、imread这几个函数都在这个工具箱里。光装基础MATLAB而不装工具箱会报Undefined function dct2。安装工具箱后可以运行ver命令确认工具箱列表。另外新版MATLAB对某些函数有推荐替代但dct2目前仍然是二维DCT的通用入口兼容性没问题。说到MATLAB版本不同版本对oo架构的支持和对图片读取格式支持略有差异但对这套加密算法来说没有任何影响R2018b之后的版本都能直接跑。有一点要提醒部分旧版本blockproc函数性能较差如果在旧版遇到运行很慢换成双层for循环也可以逻辑等价。4.2 高频报错与处理办法我整理了几个复现时最容易踩的坑直接做成速查表报错或异常可能原因解决办法Error using zeros维度错误输入图像是彩色RGBsize返回3维先rgb2gray转灰度Undefined function dct2未安装Image Processing Toolbox安装工具箱或改用自写DCT矩阵Subscript indices must either be real positive integers混沌序列排序索引出现非整数/负数检查μ是否超出混沌区间或x0取值不合法解密图像完全是噪声解密密钥与加密密钥不一致检查μ、x0是否完全相等尤其注意浮点精度解密图像方向错乱或局部错位填充区域没有裁掉或置乱恢复写法反向确认recoverVec(idx) vec的写法imshow显示全黑全白double类型矩阵带负值默认显示范围失效统一用imshow(cipher, [])图像边缘出现黑边填充导致逆变换后边界波纹用相同参数填充0解密后裁掉代码里已处理这里重点展开讲一个我最开始踩过的坑加密好的图像直接imwrite存盘再用imread读回来解密结果乱成一团。原因很简单DCT系数置乱后逆变换得到的密文是double浮点数imwrite保存成PNG时会自动把像素值量化到0~255整数这个过程是有损的原始密文的浮点信息丢失解密时逆置乱到的系数已经不是原来的数值。正确做法是加密后如果要存盘传输应该把密文矩阵归一化后转成uint8保存解密前先转回double但这样做本质上已经把加密过程从无损变成了有损。在课程设计或论文里如果要求显示“无损恢复”就别走存盘再读回这条链路直接在内存里完成加密解密。4.3 加密强度与计算效率的平衡分块大小这个参数值得多说两句。默认8×8来自JPEG标准但不同块大小对加密质量有细微影响。块越小局部能量越集中置乱粒度越细但分块数量多循环次数多效率下降块越大比如16×16或整图变换单个块的能量集中度越高置乱后的扩散效果更难预测但运算量增长。我实测下来注意力放在8×8和16×16即可两者安全性指标差异很小8×8在视觉效果上噪声分布更均匀。另一个容易被忽视的点是混沌序列的迭代起点。我之前说过序列前几项可能还没进入混沌状态所以代码里跳过了前100项。这个处理虽然让密钥空间少了一小块但避免了排序索引在初始阶段出现规律性安全性反而更高。至于计算效率分块DCT加全局置乱的总耗时主要卡在DCT变换上一张512×512的灰度图在普通笔记本上跑一遍加密大概零点几秒满足实验和演示需求。如果要进一步提速可以对置乱做向量化比如直接用reshape加sort的索引一次性完成全局置乱而不是逐块循环。我上面的主函数已经用了向量化置乱块循环交给blockproc处理整体性能在同类实现里算是不错的。最后再提醒一个彩色图像的处理思路。题目如果要求加密彩色图像不建议直接把RGB三个通道分别加密因为通道间相关性会被破坏且密文色彩杂乱无规律。规范做法是转成YCbCr空间对Y亮度通道做整套DCT加密CbCr色度通道做轻量置乱这样能兼顾安全性和效率。当然如果是大作业演示直接对RGB三个通道分别执行同一套加密流程也完全可行只是密文图像会变成三通道的彩色噪声视觉冲击力反而更强。这类基于变换域的图像加密方案后续还能扩展的方向不少比如把DCT换成小波变换、在系数域加入混沌异或扰动、采用双混沌映射组合增加密钥复杂度、结合Arnold猫映射做空间域预置乱等。我个人建议刚接触的同学先把灰度图、分块DCT、系数置乱这条主线跑通理解透每一步为什么这么做再逐步加花。跑通一遍之后你收获的不只是代码还有对频率域操作和加密框架的整体认知这套思路迁移到其他变换域加密方案上也是通用的。