做图像加密方向的复现实验时,最让人头疼的不是算法本身有多难,而是那些看似不起眼的细节:相位掩模用rand还是randn生成、DRPE加密后的复振幅要不要拆成实虚两路、压缩感知的测量矩阵到底该乘在密文上还是明文上。这些坑我踩了个遍,今天把一套能跑的“DRPE+压缩感知”安全图像加密方案完整拆开,连同MATLAB代码思路和调参经验一起分享出来。
这套方案适合这几类人:正在做光学加密或信息安全方向的研究生,需要在论文里复现DRPE与压缩感知组合算法的同学,以及想快速验证“加密同时压缩”思路的MATLAB开发者。看完之后你不仅能跑通代码,还能理解每一步设计背后的逻辑,包括怎么构造矩阵、怎么设置稀疏度、怎么设计安全分析实验。
1. 为什么非要DRPE和CS绑在一起?
1.1 DRPE的天然短板:加密很强,数据量翻倍
双随机相位编码(Double Random Phase Encoding,DRPE)是1995年由Refregier和Javidi提出的经典光学加密方案。它的核心思想很直白:在4f光学系统中,输入平面放置一个随机相位掩模,傅里叶频谱面放置另一个随机相位掩模,两个掩模共同决定密钥。明文经过这两次随机相位调制后,输出就是一幅统计特性接近白噪声的复振幅分布图。
这套机制的优势是密钥空间极大且光学实现简单,但数字仿真中它有一个很扎眼的短板:加密结果是复数值,既有实部又有虚部。一张256×256的灰度图,加密后如果要完整保存密文,需要分别存储实部矩阵和虚部矩阵,数据量直接翻倍。这在带宽受限的传输场景下非常不友好。更麻烦的是,复振幅在直方图、像素相关性等统计指标上虽然是“噪声样”,但密文尺寸和明文尺寸完全一样,这给了攻击者额外的侧信道信息——知道密文大小就能推断明文大小。
1.2 压缩感知能解决什么
压缩感知(Compressive Sensing,CS)的出现正好弥补这个短板。它的核心结论是:如果信号在某个变换域稀疏,就能用远低于奈奎斯特采样率的测量数量去采集它,并通过优化算法精确重建。对应到加密场景,就是在加密端同时对数据进行降维,把原本尺寸翻倍的复数密文压缩到一个更小的测量向量。
但必须说清楚,压缩感知不是简单的“压缩”,它是测量。理论上用m行测量矩阵去采样一个MN维信号,输出是一个m维向量,且m远小于MN。这个测量过程本身带有随机性,测量矩阵可以当作额外密钥,因此CS层等于二次加密。
1.3 常见复现方案里的逻辑坑
很多期刊论文和开源代码在融合DRPE和CS时,采用两步走的流程:先对明文做普通DRPE加密,得到复数密文;再把密文的实虚部分拆开当成普通图像,去套传统的CS框架。我第一次按这个思路复现时,解密结果惨不忍睹。原因在于:压缩感知重建的前提是信号在某个变换域稀疏,而DRPE加密后的密文分布接近随机噪声,在DCT或小波域根本不稀疏。你用OMP去重建一个不稀疏的密文向量,恢复精度会崩得很厉害。
正确的做法是建立一个统一感知模型:y = Φ·H·f,其中f是明文,H是DRPE加密对应的线性酉变换,Φ是CS测量矩阵。也就是说,测量对象虽然是加密后的数据,但重建模型直接利用明文在某个稀疏基Ψ下的稀疏性。这样OMP恢复的是明文在Ψ下的稀疏系数,而不是密文本身。后面我会详细讲这个矩阵怎么构造。
2. DRPE数字化的关键:相位掩模与傅里叶算子的矩阵化
2.1 DRPE的数学结构
设明文图像为f(x,y),两个随机相位掩模分别为R1(x,y)=exp[j2πp(x,y)]和R2(u,v)=exp[j2πq(u,v)],其中p和q是[0,1]均匀分布的随机矩阵。标准DRPE加密过程可以写成:
E(x,y) = IFT{ FT{ f(x,y)·R1(x,y) }·R2(u,v) }
也就是说:明文先乘一个空间域随机相位掩模,做傅里叶变换后在频域乘第二个随机相位掩模,再逆傅里叶变换回空间域。解密过程是逆运算:
f(x,y) = IFT{ FT{ E(x,y) }·conj(R2) }·conj(R1)
这里用到了随机相位掩模的模值为1的性质,所以除法等价于乘共轭。这是DRPE数字实现中一个非常实用的简化:不要直接写复数除法,而是用conj(R)去乘,既稳定又省计算量。
2.2 相位掩模生成:rand和randn差别很大
这个坑我不得不先说:生成相位掩模必须用rand产生均匀分布,而不是randn。原因很直观,相位p(x,y)应当是[0,1]上的均匀随机数,乘上2π后均匀覆盖整个相位圆。如果用randn生成,相位分布会集中在中心区域,相当于某些相位值出现的概率远高于其他值,密文的统计噪声特性会被破坏。
MATLAB里的标准生成方式是:
rng(2024, 'twister'); [M, N] = size(f); R1 = exp(1j * 2 * pi * rand(M, N)); rng(2025, 'twister'); R2 = exp(1j * 2 * pi * rand(M, N));这里rng的种子就是密钥的一部分。你甚至可以把两个种子的组合当作主密钥来设计,后面安全性分析时会算密钥空间。
2.3 把DRPE写成显式矩阵算子H
要让DRPE和压缩感知在同一个数学模型里工作,不能只写fft2和ifft2的流程式代码,最好把整个DRPE加密过程构建成一个大矩阵H,让密文向量等于H乘明文向量。这里需要用到DFT矩阵和Kronecker积。
一维DFT矩阵D可以用MATLAB的dftmtx(N)获得,二维FFT算子的向量化形式是F2 = kron(D, D)。fft2的向量化等价于F2乘以列向量化的图像,ifft2等价于F2'/N^2乘该向量。于是DRPE正向加密算子可以构造为:
D = dftmtx(N); F2 = kron(D, D); invF2 = F2' / (N*N); H = invF2 * diag(R2(:)) * F2 * diag(R1(:));这里的尺寸是MN×MN。对32×32图像,MN=1024,H矩阵只有1024×1024复数元素,内存占用约16MB,完全可控。如果你换128×128的图像,直接构造H就是16384×16384复数矩阵,大概4GB内存,个人电脑基本扛不住。所以教学演示建议用32×32或64×64的小图,后面我会讲怎么用函数句柄绕开大矩阵问题。
2.4 复数矩阵的存储处理
DRPE输出E是复数矩阵。直接存储复数需要两个double数组,这正好和CS测量输出叠加。实际中用MATLAB保存时,我会把复数测量值拆成实部和虚部两列,如下:
cipher_real = real(y); cipher_imag = imag(y); cipher_pack = [cipher_real, cipher_imag];这样做的好处是方便后续量化和信道传输,坏处是密文存储量的实际占用增加。不过CS已经把维度从MN降到m了,拆成实虚后总量是2m,只要满足2m < MN就还有压缩收益。换句话说,压缩率CR需要按这个规则重新核算,通常加密场景下CR取0.25到0.5之间比较合理。
3. 把压缩感知的测量方程改写成 y=ΦHΨs
3.1 可显式构造的正交稀疏基:2D-DCT
压缩感知要求信号在某个变换域是稀疏的。自然图像在DCT或小波域能量高度集中,这是JPEG和JPEG2000能压缩的原因。为了教学清晰地显式构造基矩阵,我选择2D-DCT系数作为稀疏域,主要因为它可以写成两个一维DCT矩阵的Kronecker积:
Ψ = kron(DCT_M, DCT_N)
这样Ψ矩阵是MN×MN的正交矩阵,列向量化处理很方便。不用依赖额外的图像处理工具箱,用循环自己生成一维DCT矩阵即可:
function D = mydctmtx(n) D = zeros(n, n); for k = 0:n-1 for idx = 0:n-1 if k == 0 D(k+1, idx+1) = sqrt(1/n); else D(k+1, idx+1) = sqrt(2/n) * cos(pi * (2*idx+1) * k / (2*n)); end end end end有了Ψ之后,明文和稀疏系数之间的关系是x_vec = Ψ·s,s = Ψ'·x_vec。这里Ψ正交,所以正变换和逆变换都是矩阵乘法,不需要额外工具箱。
3.2 从y=Φx到y=ΦHΨs的逻辑演进
标准CS模型是y = Φx,x是目标信号,y是测量值。现在x其实是明文f,但我们在测量前对明文施加了DRPE变换H,所以测量过程实际是y = Φ(Hf)。由于f = Ψs,最终得到:
y = Φ·H·Ψ·s = A·s
其中A = Φ·H·Ψ就是感知矩阵,s是明文在DCT域的稀疏系数。OMP重建的目标是s,得到s之后再用f = Ψ·s恢复明文。
这个改写是整个方案里最关键的一步。它避免了一个经典误区:不是先重建密文再做逆DRPE,而是把DRPE变换H作为感知矩阵的一部分,让重建算法直接利用明文在DCT域的稀疏先验。这样即使在CR较低的时候,重建质量也远好于“先重建密文再解密”的分步方案。
3.3 测量矩阵与感知矩阵的性质
测量矩阵Φ我选用高斯随机矩阵,原因是它满足受限等距性质(RIP)的概率很高。具体构造如下:
rng(2026, 'twister'); Phi = randn(m, MN) / sqrt(m);这里的除以sqrt(m)很重要,它保证Φ的列具有近似单位范数,不会因为列范数差异过大而让OMP的选择偏向某一列。m = round(CR * MN)是测量行数,CR为压缩率。
构造完A之后要注意,A本身是复数矩阵,因为H是复数的。OMP算法处理复数矩阵完全没问题,只需要把内积运算自然扩展为共轭内积。MATLAB里A' * r会自动做共轭转置,所以复数OMP的代码和实数版本长得一样,唯一要注意的是相关性计算后要取abs找最大值。
3.4 密钥体系的扩展
传统的DRPE密钥只有R1和R2两个随机相位掩模。加入CS层后,测量矩阵Φ的随机种子也可以作为一个密钥分量参与加密。哪怕攻击者拿到了密文y和感知矩阵A的一部分,缺少Φ的种子就无法重建正确的s。密钥空间可以这样估算:如果相位掩模每个像素量化到256个相位等级(8 bit),那R1和R2组合起来就有2^(8×MN×2)种可能。对32×32图像,MN=1024,这个指数是16384 bit,远超AES-256。虽然这种直接对比在严格密码学里不太严谨,但用来回答审稿人关于密钥空间的问题已经足够了。
4. 加密端MATLAB实现:从明文到复数测量值
4.1 主函数完整流程
我把整个加密过程封装成一个函数cs_drpe_encrypt,输入是明文图像f、两个随机种子、压缩率CR和稀疏度K,输出是测量值y和感知矩阵A。主流程分四步:构造DCT基、构造DRPE算子H、构造测量矩阵Φ、计算测量值y。
function [y, A, Psi, H] = cs_drpe_encrypt(f, seed1, seed2, seedPhi, CR) [M, N] = size(f); MN = M * N; x = f(:); % 构造2D-DCT稀疏基 DM = mydctmtx(M); DN = mydctmtx(N); Psi = kron(DM, DN); % 生成两个随机相位掩模 rng(seed1, 'twister'); R1 = exp(1j * 2 * pi * rand(M, N)); rng(seed2, 'twister'); R2 = exp(1j * 2 * pi * rand(M, N)); % 构造DRPE算子H D = dftmtx(N); if M == N F2 = kron(D, D); else DM_fft = dftmtx(M); F2 = kron(D, DM_fft); end invF2 = F2' / (M*N); H = invF2 * diag(sparse(R2(:))) * F2 * diag(sparse(R1(:))); % 构造高斯测量矩阵 m = max(round(CR * MN), 1); rng(seedPhi, 'twister'); Phi = randn(m, MN) / sqrt(m); % 统一感知矩阵与测量 A = Phi * H * Psi; y = A * (Psi' * x); end这里我特意用了sparse构造对角阵,因为R1和R2的对角阵虽然大,但只有对角线非零,用sparse能节省大量内存和乘法计算。
4.2 为什么感知矩阵A要在加密端和解密端同时持有
注意我在函数里返回了A,这意味着接收方必须持有与发送方相同的A才能重建。A里面包含R1、R2和Φ的全部信息,也就是完整密钥。在实际传输中,发送方没有必要把整个A发给接收方,只需要共享三个种子和CR、K这些公共参数,接收方本地重新生成A即可。种子就是密钥,这个设计思路和对称加密体系很像。
4.3 密文量的平衡与参数选择
CR这个参数需要谨慎设置。实测下来,CR=0.5时测量值m=MN/2,拆成实虚部后总存储量恰好等于一个MN大小的实数矩阵,也就是和原始明文尺寸持平,这是“加密但不膨胀”的临界点。CR低于0.5时,端口传输的数据量会小于明文;CR高于0.5时,数据量超过明文,压缩感知的意义就变小了。所以我一般推荐CR取0.3到0.5之间。
稀疏度K则是OMP迭代次数的上界,需要和图像内容匹配。32×32的测试图,DCT系数能量集中在前100~200个系数,K取120比较稳妥。K太小会丢失高频细节,图像变模糊;K太大会引入噪声,重建结果出现振铃。这个值可以在加密端先算一下能量累积比例再定。
4.4 小图测试与真实图像生效的差异
用32×32小图验证算法时,所有矩阵都能显式构造,调试起来很舒服。但实验报告里总不能只放32×32的模糊图。要处理128×128以上的图像,H和Ψ的显式矩阵会非常巨大,需要在迭代算法里引入函数句柄,避免存储整个A矩阵。思路是把Ax和A'y拆成三个级联操作,例如Ax = Phi(H*(Psi*s)),先用函数形式实现Phi、H、Psi各自的作用,再用随机Kaczmarz或共轭梯度法做重建。篇幅关系这里不展开完整代码,但只要理解了本文的矩阵模型,改造方向是明确的。
5. 解密端OMP重建与图像恢复的完整链路
5.1 复数OMP的实现细节
正交匹配追踪(OMP)是最直观的稀疏重建算法:迭代地找出与残差最相关的原子,把该原子加入支撑集,用最小二乘更新系数,再更新残差。复数场景下唯一需要留意的就是相关性的度量,用A'*r然后取模值找最大值,而不是直接取内积实部。一个健壮的OMP实现如下:
function s_hat = omp_cs(y, A, K) [m, n] = size(A); r = y; supp = []; x_hat = zeros(n, 1); for iter = 1:K corr = A' * r; [~, pos] = max(abs(corr)); if ismember(pos, supp) corr(pos) = 0; [~, pos] = max(abs(corr)); end supp = sort([supp, pos]); A_s = A(:, supp); coef = A_s \ y; x_hat = zeros(n, 1); x_hat(supp) = coef; r = y - A * x_hat; end s_hat = x_hat; end这里最小的防御逻辑是处理重复原子。虽然理论情况下OMP不会重复选择同一个原子,但由于浮点误差或者感知矩阵列之间的微小相关性,重复选择偶有发生。如果重复了也不做处理,支撑集反复加入同一列,最小二乘矩阵会奇异,重建直接失败。
5.2 解密主流程
解密端的完整步骤是:用OMP从y和A恢复DCT稀疏系数s_hat,然后乘以Ψ得到明文向量,再reshape减去虚部残余。注意我保留了少量虚部残余处理,因为复数最小二乘可能引入微小虚部。
function f_rec = cs_drpe_decrypt(y, A, Psi, K, M, N) s_hat = omp_cs(y, A, K); x_hat = Psi * s_hat; f_rec = reshape(real(x_hat), M, N); f_rec = max(min(f_rec, 1), 0); end最后的max/min裁剪是把重建值约束到[0,1]区间。这一步很多人漏掉,结果PSNR算出来很高,但图像上有负值或超过255的异常点,显示时直接白花花的。
5.3 质量评估:PSNR和SSIM
重建质量评估用PSNR和SSIM两个指标就够了。PSNR关注像素级误差,SSIM关注结构相似度。计算代码如下:
function [psnr_val, ssim_val] = evaluate(f_orig, f_rec) mse_val = mean((f_orig(:) - f_rec(:)).^2); psnr_val = 10 * log10(1 / mse_val); % 图像归一化到[0,1] ssim_val = ssim(f_orig, f_rec); % 需要Image Processing Toolbox end实测经验:CR=0.5、K=120时,对标准cameraman图重建PSNR大约在30~33dB,SSIM在0.95左右。CR降到0.25时PSNR会掉到25dB以下,图像边缘开始发糊。这时候不要盲目调大K,K过大反而会把重建野值引进来。
5.4 重建失败时的排查顺序
如果解密图像一团乱,按照这个顺序排查:第一步确认加密端的A和解密端的A是否一致,种子有没有传错;第二步打印s_hat的长度和支撑集数量,确认OMP确实迭代了K次;第三步看r的残差能量是否在下降,如果残差大且不下降,多半是CR太低或者K太小;第四步检查f_rec是否有严重的负值分布,如果有,说明最小二乘解不稳定,考虑用QR解替代左除。这些排查步骤比直接重新写代码要高效得多。
6. 安全性测试怎么设计:密钥、统计与鲁棒性
6.1 密钥敏感性测试
密码算法的安全性,首先看密钥敏感性。测试方法很简单:分别用正确的种子和改掉一位的种子解密同一份密文,观察解密图的PSNR。正确密码重建PSNR应在30dB左右,错误密码重建PSNR应低于10dB,肉眼看起来完全是噪声。我常用一个技巧:修改种子后,得到的解不仅视觉是噪声,其分布还接近均匀随机,这样就能体现DRPE+CS组合对密钥的完全依赖。注意错误密钥解密时,解密算法本身不会报错,它照样完成OMP迭代,只是A的列不再包含正确的DRPE信息,重建的系数会在DCT域四处发散。
6.2 统计攻击分析
统计攻击主要看密文的两个指标:直方图和相邻像素相关性。DRPE加密后的密文测量值y,直方图应当接近高斯或均匀分布,没有明显的峰谷结构;相邻元素间的皮尔逊相关系数应该接近0。对复数测量值y计算相关系数时,可以先拆成实虚部,分别统计实部相邻元素相关性和虚部相邻元素相关性。实测中,CR=0.5时实部和虚部的相邻相关系数都能控制在0.05以下,肉眼和计算都看不出明文结构。
另一个常用指标是信息熵。密文信息熵接近理论最大值时,说明密文不确定性高。对量化到8bit的密文,理论最大熵是8bit;实际中DRPE密文熵能到7.9以上。这套指标组合写进论文“抗统计攻击分析”一节是够用的。
6.3 抗裁剪与抗噪声鲁棒性
光学加密系统经常遇到信道噪声和数据裁切的问题,所以还要测试鲁棒性。压缩感知本身自带一定容错能力:测量值y在传输中遭到部分损坏,只要损坏比例不太高,OMP重建依然能恢复主体信息。测试方法是对y随机加高斯白噪声,或者随机置零一部分元素,幅度从5%到20%递增,观察PSNR变化曲线。
实测下来,y被随机置零30%时,重建PSNR仍能维持在20dB左右,但会明显出现块状伪影。这个鲁棒性来自CS测量的全局性质:每个测量值都包含整个明文的信息,局部丢失不会立刻摧毁全部内容。
6.4 关于已知明文攻击的简单讨论
严格的安全性还需要考虑已知明文攻击。所谓已知明文攻击,就是攻击者拿到一组明文和对应密文,想反推密钥。在DRPE+CS框架里,如果攻击者知道明密文对,理论上可以构造出关于Φ和H的部分约束方程,但H中包含R1和R2两个随机矩阵,每个像素相位连续变化,约束方程的非线性很强,目前公开文献里还没有看到能在密钥空间不缩减前提下有效求解的方法。这也是这类光学加密方案长期活跃的学术原因。不过做工程落地的话,我的建议是不要把DRPE的随机种子当作固定长期密钥,最好采用“一次一密”的会话密钥,即每次加密都重新生成种子,用密码学手段先安全协商种子,再用于图像加密。这样可以规避很多实际攻击场景。
落到实处的经验汇总
最后聊几句贴近操作的心得。第一次复刻这套方案时,我一度把CR设置为0.75,理由是“压缩率越高越好”,结果解密图像全是噪点。后来才意识到CR是测量次数占比,CR=0.75意味着测量矩阵没有把数据压下来,反而因为感知矩阵更庞大,数值条件数变差,重建更容易病态。所以CR并不是越高越好,也不是越低越好,0.3到0.5是常见甜点区间。
另外一定要记住,加密对象是二维图像,但在模型里它被列向量化了。很多人在调试时分不清vecot和矩阵的reshape,导致解密图像横竖方向对不上。我的建议是所有中间变量一律按MN×1的列向量处理,只有最终恢复明文时才reshape成M×N。
这套方案如果要扩展到实际系统,下一步就是我把显式矩阵H替换成函数句柄形式,把DCT基换成小波基,再用真实光学系统采集的数据测试。只要理解本文的y=ΦHΨs这个统一模型,各种扩展都只是工程细节,核心逻辑不会变。