做图像加密方向的复现实验时最让人头疼的不是算法本身有多难而是那些看似不起眼的细节相位掩模用rand还是randn生成、DRPE加密后的复振幅要不要拆成实虚两路、压缩感知的测量矩阵到底该乘在密文上还是明文上。这些坑我踩了个遍今天把一套能跑的“DRPE压缩感知”安全图像加密方案完整拆开连同MATLAB代码思路和调参经验一起分享出来。这套方案适合这几类人正在做光学加密或信息安全方向的研究生需要在论文里复现DRPE与压缩感知组合算法的同学以及想快速验证“加密同时压缩”思路的MATLAB开发者。看完之后你不仅能跑通代码还能理解每一步设计背后的逻辑包括怎么构造矩阵、怎么设置稀疏度、怎么设计安全分析实验。1. 为什么非要DRPE和CS绑在一起1.1 DRPE的天然短板加密很强数据量翻倍双随机相位编码Double Random Phase EncodingDRPE是1995年由Refregier和Javidi提出的经典光学加密方案。它的核心思想很直白在4f光学系统中输入平面放置一个随机相位掩模傅里叶频谱面放置另一个随机相位掩模两个掩模共同决定密钥。明文经过这两次随机相位调制后输出就是一幅统计特性接近白噪声的复振幅分布图。这套机制的优势是密钥空间极大且光学实现简单但数字仿真中它有一个很扎眼的短板加密结果是复数值既有实部又有虚部。一张256×256的灰度图加密后如果要完整保存密文需要分别存储实部矩阵和虚部矩阵数据量直接翻倍。这在带宽受限的传输场景下非常不友好。更麻烦的是复振幅在直方图、像素相关性等统计指标上虽然是“噪声样”但密文尺寸和明文尺寸完全一样这给了攻击者额外的侧信道信息——知道密文大小就能推断明文大小。1.2 压缩感知能解决什么压缩感知Compressive SensingCS的出现正好弥补这个短板。它的核心结论是如果信号在某个变换域稀疏就能用远低于奈奎斯特采样率的测量数量去采集它并通过优化算法精确重建。对应到加密场景就是在加密端同时对数据进行降维把原本尺寸翻倍的复数密文压缩到一个更小的测量向量。但必须说清楚压缩感知不是简单的“压缩”它是测量。理论上用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图像MN1024H矩阵只有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Ψs3.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(k1, idx1) sqrt(1/n); else D(k1, idx1) sqrt(2/n) * cos(pi * (2*idx1) * k / (2*n)); end end end end有了Ψ之后明文和稀疏系数之间的关系是x_vec Ψ·ss Ψ·x_vec。这里Ψ正交所以正变换和逆变换都是矩阵乘法不需要额外工具箱。3.2 从yΦx到yΦHΨs的逻辑演进标准CS模型是y Φxx是目标信号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图像MN1024这个指数是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这个参数需要谨慎设置。实测下来CR0.5时测量值mMN/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和Ay拆成三个级联操作例如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实测经验CR0.5、K120时对标准cameraman图重建PSNR大约在30~33dBSSIM在0.95左右。CR降到0.25时PSNR会掉到25dB以下图像边缘开始发糊。这时候不要盲目调大KK过大反而会把重建野值引进来。5.4 重建失败时的排查顺序如果解密图像一团乱按照这个顺序排查第一步确认加密端的A和解密端的A是否一致种子有没有传错第二步打印s_hat的长度和支撑集数量确认OMP确实迭代了K次第三步看r的残差能量是否在下降如果残差大且不下降多半是CR太低或者K太小第四步检查f_rec是否有严重的负值分布如果有说明最小二乘解不稳定考虑用QR解替代左除。这些排查步骤比直接重新写代码要高效得多。6. 安全性测试怎么设计密钥、统计与鲁棒性6.1 密钥敏感性测试密码算法的安全性首先看密钥敏感性。测试方法很简单分别用正确的种子和改掉一位的种子解密同一份密文观察解密图的PSNR。正确密码重建PSNR应在30dB左右错误密码重建PSNR应低于10dB肉眼看起来完全是噪声。我常用一个技巧修改种子后得到的解不仅视觉是噪声其分布还接近均匀随机这样就能体现DRPECS组合对密钥的完全依赖。注意错误密钥解密时解密算法本身不会报错它照样完成OMP迭代只是A的列不再包含正确的DRPE信息重建的系数会在DCT域四处发散。6.2 统计攻击分析统计攻击主要看密文的两个指标直方图和相邻像素相关性。DRPE加密后的密文测量值y直方图应当接近高斯或均匀分布没有明显的峰谷结构相邻元素间的皮尔逊相关系数应该接近0。对复数测量值y计算相关系数时可以先拆成实虚部分别统计实部相邻元素相关性和虚部相邻元素相关性。实测中CR0.5时实部和虚部的相邻相关系数都能控制在0.05以下肉眼和计算都看不出明文结构。另一个常用指标是信息熵。密文信息熵接近理论最大值时说明密文不确定性高。对量化到8bit的密文理论最大熵是8bit实际中DRPE密文熵能到7.9以上。这套指标组合写进论文“抗统计攻击分析”一节是够用的。6.3 抗裁剪与抗噪声鲁棒性光学加密系统经常遇到信道噪声和数据裁切的问题所以还要测试鲁棒性。压缩感知本身自带一定容错能力测量值y在传输中遭到部分损坏只要损坏比例不太高OMP重建依然能恢复主体信息。测试方法是对y随机加高斯白噪声或者随机置零一部分元素幅度从5%到20%递增观察PSNR变化曲线。实测下来y被随机置零30%时重建PSNR仍能维持在20dB左右但会明显出现块状伪影。这个鲁棒性来自CS测量的全局性质每个测量值都包含整个明文的信息局部丢失不会立刻摧毁全部内容。6.4 关于已知明文攻击的简单讨论严格的安全性还需要考虑已知明文攻击。所谓已知明文攻击就是攻击者拿到一组明文和对应密文想反推密钥。在DRPECS框架里如果攻击者知道明密文对理论上可以构造出关于Φ和H的部分约束方程但H中包含R1和R2两个随机矩阵每个像素相位连续变化约束方程的非线性很强目前公开文献里还没有看到能在密钥空间不缩减前提下有效求解的方法。这也是这类光学加密方案长期活跃的学术原因。不过做工程落地的话我的建议是不要把DRPE的随机种子当作固定长期密钥最好采用“一次一密”的会话密钥即每次加密都重新生成种子用密码学手段先安全协商种子再用于图像加密。这样可以规避很多实际攻击场景。落到实处的经验汇总最后聊几句贴近操作的心得。第一次复刻这套方案时我一度把CR设置为0.75理由是“压缩率越高越好”结果解密图像全是噪点。后来才意识到CR是测量次数占比CR0.75意味着测量矩阵没有把数据压下来反而因为感知矩阵更庞大数值条件数变差重建更容易病态。所以CR并不是越高越好也不是越低越好0.3到0.5是常见甜点区间。另外一定要记住加密对象是二维图像但在模型里它被列向量化了。很多人在调试时分不清vecot和矩阵的reshape导致解密图像横竖方向对不上。我的建议是所有中间变量一律按MN×1的列向量处理只有最终恢复明文时才reshape成M×N。这套方案如果要扩展到实际系统下一步就是我把显式矩阵H替换成函数句柄形式把DCT基换成小波基再用真实光学系统采集的数据测试。只要理解本文的yΦHΨs这个统一模型各种扩展都只是工程细节核心逻辑不会变。