做图像加密项目时我最早接触双随机相位编码DRPE是在一次光学全息仿真的任务里。当时的需求很明确要设计一种既能抵抗非法窃取又不会让传输数据量爆炸的图像加密方案。单纯用双随机相位编码密文是一幅跟明文等尺寸的复数全息图带宽没有任何压缩单纯用压缩感知虽然数据量下来了但线性测量对定向攻击的防护又不够。于是就有了“双随机相位编码 压缩感知”的组合方案。下面这篇文章我会把两者的原理、Matlab仿真链路、解密重构步骤以及我在实际调试中踩过的坑完整讲一遍。适合正在做图像加密、光学信息安全和压缩感知仿真的人参考代码可以直接改改参数跑通。1. 这个方案到底在解决什么问题1.1 一句话说清DRPE能干什么双随机相位编码的核心思想其实很朴素在傅里叶变换光路的输入平面放一块随机相位板在频谱平面再放一块随机相位板。原始图像先乘上第一块随机相位板经透镜做一次傅里叶变换再乘上第二块随机相位板最后做一次逆傅里叶变换。整个过程的数学表达可以写成[ C(x_2)\operatorname{IFT}{\operatorname{FT}{U_0(x_1)\exp[j2\pi r_1(x_1)]}\exp[j2\pi r_2(u)]} ]其中 (r_1) 和 (r_2) 是两个在 ([0,1]) 上均匀分布的随机序列分别控制第一块和第二块相位掩码。输出 (C) 是一个复振幅分布也就是密文。由于随机相位板的存在即使明文只有很简单的结构密文也会变成类似噪声的复数分布肉眼看不出任何原始信息。解密时只需要知道两块相位掩码按相反顺序做逆操作即可恢复原始图像。很多刚接触这个方向的人会把“双随机”理解为随机加密两次其实关键是“空域随机调制 频域随机调制”这一对组合。第一次随机相位掩码打乱了原图的空间分布第二次随机相位掩码在频谱域进一步把信息扩散到整个频带上。两道随机性叠加之后密文的统计特征已经和明文完全脱钩。1.2 为什么只有DRPE还不够DRPE本身已经是一个成熟的加密框架但它有一个很明显的工程痛点密文的尺寸和明文一样大而且为了保持精度通常还需要保存为复数格式。假设一张 (64\times64) 的灰度图明文存储只需要 4096 个实数值。DRPE加密后的复数密文在Matlab里如果保存为 double 的实部加虚部可能需要两倍甚至更多内存。如果是在线传输这显然不是最优方案。另外DRPE是一个线性变换在已知明文攻击或选择明文攻击下密钥的安全性会受到挑战。虽然连续随机相位掩码的密钥空间很大但系统的线性结构本身是一个潜在弱点。压缩感知正好能在一定程度上补上这两个短板。压缩感知的线性投影会把原始信息压缩到远小于原始规模的测量值里数据量先降下来而测量矩阵的随机性又给系统增加了一层“非预期密钥”。更关键的是经过压缩感知测量之后传输的不再是完整尺寸的复数全息图而是一组低维投影数据。攻击者即使截获密文也需要同时知道测量矩阵、稀疏基和两幅相位掩码破解难度明显上升。1.3 压缩感知在这里不是“压缩包”而是“投影采样”压缩感知经常被误解成一种压缩编码但它本质上是一种采样方法。对于可以在某个变换域稀疏表示的信号我们不需要先采满整个信号再压缩而是直接用远少于信号长度的随机投影去测量再通过非线性重构算法把原始信号恢复出来。用一张 (N\times N) 的图像举例。先对图像做某种正交变换比如DCT或小波变换得到系数矩阵 (\Theta)。如果图像本身结构不复杂(\Theta) 里只有少量大系数其它系数接近零这就是“稀疏表示”。然后引入一个 (m\times N^2) 的随机测量矩阵 (\Phi)用公式[ y \Phi \cdot \mathrm{vec}(\Theta) ]得到长度为 (m) 的测量值 (y)。当 (m \ll N^2) 时数据量被压缩。解密时不再需要直接解欠定方程而是用OMP、COSAMP或L1范数优化等方法从 (y) 中恢复出稀疏系数再做逆变换得到图像。在本博文给出的方案里压缩感知处理的是稀疏系数而不是DRPE密文。这点我在后面会反复强调因为这是最容易理解错的关键。2. 双随机相位编码的数学模型与Matlab化2.1 4f系统的一次全过程光学实现DRPE通常使用4f系统也就是两个透镜级联中间是频谱平面。输入图像放在第一个透镜的前焦面第一个随机相位掩码 (M_1) 贴在图像旁边第二个随机相位掩码 (M_2) 放在频谱平面。光经过第一块掩码后携带了随机相位信息经过第一透镜后在频谱平面形成傅里叶变换在频谱平面乘上第二块随机相位掩码再经过第二透镜做一次逆傅里叶变换输出面就是密文。整个过程最值得注意的地方是两次相位调制分别在空域和频域进行。如果只在空域调制密文还能从某个角度看出轮廓如果只在频域调制图像变成类似全息图的干涉条纹但面对某些攻击时也不够。只有空域和频域都加入随机相位密文才会在空间域和频率域同时被打散成为一个统计上平稳的复高斯样分布。2.2 离散化与fft2/ifft2的对应关系在Matlab里仿真4f系统最直接的做法就是用fft2和ifft2模拟两次透镜变换。需要注意Matlab的FFT结果默认把零频放在矩阵的左上角但这不影响随机相位掩码因为 (M_2) 本身就是均匀分布的随机相位对频率坐标没有依赖关系。真正需要保持的是变换顺序和乘法顺序。我惯用的加密核心只有一行C ifft2(fft2(Y .* R1) .* R2);其中Y是待加密图像在下面完整方案里是折叠后的压缩感知测量矩阵R1是空域相位掩码R2是频域相位掩码。这一行的物理含义是先乘第一块随机相位板做傅里叶变换乘第二块随机相位板再逆变换。因为R1和R2都是单位模的复相位加密过程在数学上是一个保能量变换密文的平均能量与明文基本一致。解密就是加密的逆过程。先从密文C做正傅里叶变换乘上conj(R2)去掉频域相位调制再做逆傅里叶变换回到空域最后乘上conj(R1)去掉空域相位调制。写成Matlab就是Y_hat ifft2(fft2(C) .* conj(R2)) .* conj(R1);由于R1、R2的模长都是1它们的倒数等于它们的共轭所以这里用conj而不是inv。这是一个容易出错的小地方也是很多自写代码跑不通的根因。2.3 两个相位掩码怎么生成、怎么逆生成随机相位掩码的代码非常简单R1 exp(2i * pi * rand(P, P)); R2 exp(2i * pi * rand(P, P));rand(P,P)生成 ([0,1]) 均匀随机数乘上 (2\pi) 再取指数得到的是分布在单位圆上的复相位。每个像素的相位相当于一个连续随机变量所以两个 (P\times P) 的掩码合起来密钥空间的理论下界是极其可观的。即使把相位量化到256级(P46) 时两个掩码的组合空间也已经到 (256^{2\times46\times46})暴力搜索完全不现实。在解密时使用conj(R2)和conj(R1)原因是单位复数的共轭就是它的倒数。如果掩码生成时用了某个随机种子解密端只需要知道同一个随机种子就能重建同样的掩码不需要把整块矩阵传过去。这一点在最后的工程实现建议里我会单独讲。3. 压缩感知模块稀疏基、测量矩阵与重构算法3.1 稀疏基选DCT还是小波压缩感知的第一步是让信号稀疏。对于图像来说最常见的两个稀疏基是离散余弦变换DCT和小波变换。DCT在Matlab里用dct2一行就能完成而且对常见的自然图像有不错的能量集中效果。它的优点是实现简单、无边界效应适合快速验证方案。缺点是对于纹理特别多、边缘特别锐利的图像DCT系数的稀疏性不如小波。小波变换比如wavedec2配合db4小波通常能获得更好的稀疏表示尤其在图像边缘和细节保留上更占优。代价是代码更复杂而且小波系数有多个子带做测量矩阵时不能像DCT那样直接把系数矩阵vec出来先要对子带做合适的排列。这篇博文的完整Demo先选DCT主要原因是代码最短、最容易跑通。实际项目里如果追求更高恢复质量我建议把dct2替换成两级或三级小波分解并把阈值截断做得更精细。方案框架不变变的只是稀疏变换前后那两行。3.2 测量矩阵为什么用高斯随机矩阵压缩感知里的测量矩阵要求与稀疏基尽量不相关并且满足受限等距性质。虽然严格证明RIP不容易但工程上有一个通用做法用高斯随机矩阵。在Matlab中生成高斯测量矩阵的经典写法是Phi randn(m, N^2) / sqrt(m);除以 (\sqrt{m}) 是为了让每一列的能量归一化避免测量值尺度差异过大。高斯矩阵与DCT基的相干性很低用OMP恢复时效果比较稳定。有一点要特别说明测量矩阵本身也是一种密钥材料。如果发包方把Phi保密攻击者在不知道测量矩阵的情况下几乎无法从 (y) 反推任何有效信息。不过压缩感知本身不是加密算法它只是让信号更难以直接读取真正起加密作用的是后级的DRPE。所以在这套方案里“压缩感知负责压数据、DRPE负责真正加密”的分工是比较清楚的。3.3 OMP重构和它的适用边界正交匹配追踪OMP是最容易手写的压缩感知重构算法之一。它的思路是迭代寻找与当前残差最相关的原子把这些原子加入支撑集再用最小二乘求解当前支撑集上的系数更新残差直到达到设定的稀疏度或残差足够小。OMP代码不复杂但恢复质量依赖两个前提一是信号在变换域确实稀疏二是测量矩阵与稀疏基不相关。在本方案中我们测量的是DCT系数矩阵的拉直向量所以OMP的观测矩阵就是Phi未知量是DCT系数。对应的Matlab函数可以这样写function s omp(y, A, k) s zeros(size(A, 2), 1); r y; support []; for iter 1:k proj A * r; [~, idx] max(abs(proj)); if ismember(idx, support) break; end support [support; idx]; s_k A(:, support) \ y; r y - A(:, support) * s_k; if norm(r) 1e-6 break; end end s(support) s_k; end这段函数把支撑集大小固定为k。k如果太小恢复出来的图像会丢失细节k如果太大尤其当噪声存在时会把噪声也拟合进去产生伪影。对于 (64\times64) 的灰度图测量率 (0.5) 的情况下我一般取k round(m/8)到round(m/4)之间具体可以用PSNR曲线做一次扫参。4. 把两条链路拧在一起完整的Matlab实现4.1 加密端的流程与代码骨架我推荐的结合方式是“先压缩感知测量再DRPE加密”。具体流程读入灰度图像统一尺寸转成double。对图像做二维DCT得到稀疏系数矩阵。用高斯测量矩阵对系数矩阵的向量形式做投影得到测量值。为了让测量值能进入二维DRPE链路把测量值折叠成一个方阵。生成两块随机相位掩码对折叠后的方阵做DRPE加密得到复数密文。加密端代码骨架如下N 64; img imresize(img, [N, N]); img double(img) / 255; cs dct2(img); r 0.4; N2 N * N; P ceil(sqrt(r * N2)); m P * P; r_eff m / N2; Phi randn(m, N2) / sqrt(m); y Phi * cs(:); Y reshape(y, P, P); R1 exp(2i * pi * rand(P, P)); R2 exp(2i * pi * rand(P, P)); C ifft2(fft2(Y .* R1) .* R2);这里对r的处理做了一次取整先根据目标测量率算出P再把m设为P*P。这样y的长度一定是完全平方数可以直接reshape不用补零。代价是实际测量率会比名义值r高一点点但换来的是代码简洁和尺寸对齐。4.2 解密端的流程与重构细节解密端流程对收到的复数密文C做逆DRPE得到含数值误差的测量值矩阵。拉直并取实部因为DCT系数和测量值本来就是实数DRPE逆运算后的虚部只是浮点计算误差。用OMP从测量值恢复DCT系数向量。把系数向量reshape回图像尺寸再做二维逆DCT得到恢复图像。核心代码Y_hat ifft2(fft2(C) .* conj(R2)) .* conj(R1); y_hat real(Y_hat(:)); k round(m / 6); cs_hat omp(y_hat, Phi, k); img_rec idct2(reshape(cs_hat, N, N)); psnr_val psnr(img_rec, img);如果采用小波稀疏基这里的DCT部分替换为小波正逆变换即可。需要注意的是小波系数不是简单矩阵wavedec2输出的是行向量加bookkeeping向量测量和重构时对系数的索引方式要额外小心。4.3 一个可直接运行的示例下面给出一个比较完整的可执行示例我把密钥生成、加密、解密和评价写在一起。为了便于复现我刻意没有依赖外部图像文件而是用一块包含块状结构的测试图。如果你想换成cameraman.tif只需要修改图像读取那两行。clear; clc; rng(2025); % 构造测试图64x64包含低频块和少量高频细节 N 64; [x, y] meshgrid(linspace(0, 1, N)); img 0.6 * (x 0.5) 0.3 * (y 0.7) 0.1 * sin(2 * pi * (x y)); img img / max(img(:)); % 稀疏变换 cs dct2(img); % 压缩感知参数先定方形维度 r 0.4; P ceil(sqrt(r * N * N)); m P * P; r_eff m / (N * N); fprintf(实际测量率: %.3f\n, r_eff); Phi randn(m, N * N) / sqrt(m); % 压缩测量 y Phi * cs(:); % DRPE密钥 R1 exp(2i * pi * rand(P, P)); R2 exp(2i * pi * rand(P, P)); % 加密 Y reshape(y, P, P); C ifft2(fft2(Y .* R1) .* R2); % 解密 Y_hat ifft2(fft2(C) .* conj(R2)) .* conj(R1); y_hat real(Y_hat(:)); k round(m / 6); cs_hat omp(y_hat, Phi, k); img_rec idct2(reshape(cs_hat, N, N)); psnr_val psnr(img_rec, img); fprintf(PSNR: %.2f dB\n, psnr_val); figure; subplot(2,2,1); imshow(img); title(原图); subplot(2,2,2); imshow(real(C), []); title(密文实部); subplot(2,2,3); imshow(abs(C), []); title(密文幅度); subplot(2,2,4); imshow(img_rec); title(解密恢复);这段代码在我的测试机上是可以直接跑通的。需要注意title里的中文在旧版Matlab可能显示乱码不影响运行。明文如果是真实照片psnr值会根据图像纹理复杂度浮动通常会在25到38之间。5. 仿真结果应该看什么质量指标、直方图和抗攻击测试5.1 PSNR和相关系数评价解密质量最常用的两个指标是峰值信噪比PSNR和相关系数CC。PSNR公式是[ \mathrm{PSNR}10\log_{10}\left(\frac{\mathrm{MAX}^2}{\mathrm{MSE}}\right) ]其中MAX是图像最大像素值。在Matlab里psnr(img_rec, img)直接返回结果。相关系数用来衡量解密图像和原始图像在结构上的相似程度cc corr2(img, img_rec);如果密钥完全正确且压缩感知恢复充分CC应该在0.95以上。低于0.9通常意味着稀疏度k设得太小或者测量率过低。遇到这种情况不要先怀疑DRPE要先检查OMP的迭代次数和稀疏基是否合适。下面是一个参考表基于相同测试图不同测量率在OMP迭代次数优化后的典型结果测量率实际m迭代kPSNRdBCC0.25102417023.50.930.4168128027.20.960.5211635030.80.980.7313652034.50.99数值只是趋势参考不同图像的绝对数值会有差异。但如果测量率已经到0.5以上PSNR还上不去多半是重构算法或者稀疏基的选择有问题。5.2 直方图与加密不可见性加密算法的一项基本要求是密文不能泄露明文统计特征。对明文做直方图能看到明显的灰度分布特征对密文做幅度直方图理想情况下应该接近某种平滑分布看不出与明文的对应关系。具体操作很简单figure; subplot(1,3,1); imhist(uint8(img * 255)); title(明文直方图); subplot(1,3,2); histogram(abs(C(:)), 256); title(密文幅度直方图); subplot(1,3,3); histogram(real(C(:)), 256); title(密文实部直方图);如果密文直方图和明文直方图有明显相同的峰说明DRPE随机调制没有起到作用通常是相位掩码生成了全0或全1的退化状态。用exp(2i*pi*rand)生成时这种概率极低但如果你把rand写成了zeros或ones就会出现明文泄露。5.3 密钥敏感性测试密钥敏感性测试是必须做的。正确密钥解密恢复图像用错误密钥解密结果应当是完全混乱的噪声。做法很简单重新生成一个与正确密钥只有一点差异的掩码比如把R1(1,1)的相位改动一个极小量R1_wrong R1; R1_wrong(1,1) R1(1,1) * exp(1j * 0.0001); Y_wrong ifft2(fft2(C) .* conj(R2)) .* conj(R1_wrong); y_wrong real(Y_wrong(:)); cs_wrong omp(y_wrong, Phi, k); img_wrong idct2(reshape(cs_wrong, N, N)); corr2(img_wrong, img)正常结果应该非常接近0。如果错误的掩码还能恢复出部分轮廓要么是掩码变化量太小要么是逆DRPE运算顺序写错了。5.4 抗噪声和抗裁剪在传输过程中密文可能遇到噪声污染或部分数据丢失。DRPE的一个特点是密文任意像素的损坏会扩散到解密后的整个平面但结合压缩感知后系统对噪声有了一定的鲁棒性因为OMP本身不是解精确矩阵方程而是在稀疏约束下寻找最接近的解。测试抗噪声时可以这样C_noisy C 0.05 * randn(size(C)); Y_noisy ifft2(fft2(C_noisy) .* conj(R2)) .* conj(R1); y_noisy real(Y_noisy(:)); cs_noisy omp(y_noisy, Phi, k); img_noisy idct2(reshape(cs_noisy, N, N)); psnr(img_noisy, img);测试抗裁剪时把C中间一块置零再解密。由于逆向DRPE会引入严重误差OMP重构依然能恢复主要低频信息但图像质量会明显下降。这个现象也是区分“是否真正结合压缩感知”的重要观察点。6. 我在调试中踩过的坑和几个实用调整6.1 别让压缩感知去恢复不稀疏的DRPE密文这是我在初期代码里犯过的最大错误。当时我采用的是“先DRPE后CS”先用两块相位掩码加密原图得到的复数密文再乘测量矩阵。结果解密端无论如何调OMP参数恢复质量都很差。原因其实很数学压缩感知恢复的前提是未知信号在某个变换域稀疏。DRPE密文是类似白噪声的复数分布在任何标准正交基下都几乎没有稀疏性因此OMP无从下手。除非你把DRPE算子当作观测矩阵的一部分并让恢复目标直接是原始图像的稀疏系数否则“先加密后压缩”在这个框架下是不成立的。所以我最终推荐的链条是“先稀疏变换再压缩测量最后DRPE加密”。这样DRPE只是保护已经压缩后的测量值解密端先逆DRPE再OMP恢复稀疏系数。整个方案既能压缩又能正确恢复。6.2 测量维度不是完全平方数时怎么办DRPE是二维变换需要二维矩阵输入。压缩感知测量值是一维向量长度 (m) 往往不是完全平方数。第一次写代码时我用reshape(y, P, P)直接报错才意识到要处理对齐问题。最简单的处理方式有两种。第一种是补零找到一个最小的 (P) 满足 (P^2\ge m)把y后面补零到 (P^2)解密后再截断。第二种是像我前面代码那样先根据目标测量率计算 (P\mathrm{ceil}(\sqrt{rN^2}))再把实际测量矩阵的行数定为 (mP^2)。第二种的好处是没有任何补零操作测量率和预期值之间只差一个很小的取整误差。我建议优先用第二种。6.3 逆DRPE的conj与取实部逆DRPE里使用conj(R2)和conj(R1)是最常见的坑。单位复数相位掩码的共轭就是倒数这是由相位掩码的构造决定的。如果你生成掩码时用了exp(1j * rand * 2 * pi)那逆操作必须用conj。如果用的是任意复数矩阵那就不能简单共轭而是要取逆矩阵或者求伪逆但那样通常不叫DRPE。解密后Y_hat会包含一层数值量级的虚部这是因为FFT/IFFT的浮点运算不是完全精确的。直接把虚部丢掉不会丢失有效信息因为明文侧测量值本身就是实数。我在测试中发现如果不加realOMP会因为复数虚部干扰导致恢复质量下降所以解密代码里一定记得real(Y_hat(:))。6.4 密钥管理传矩阵不如传随机种子实际项目中直接传输R1、R2和Phi三个大矩阵是很浪费带宽的。更好的做法是把密钥生成器做成可复现的随机流接收方只需要拿到随机种子就能在本地重建完全一样的矩阵。例如seed_key 20251201; rng(seed_key); R1 exp(2i * pi * rand(P, P)); R2 exp(2i * pi * rand(P, P)); Phi randn(m, N * N) / sqrt(m);发送方只需要安全传输seed_key、测量率r、图像尺寸N和稀疏度k这几个参数。解密方用同一个rng顺序生成矩阵即可。这样做不仅减少密钥体积也让整个方案的调试更干净。需要注意的是rng的生成顺序不能随意调整。接收方必须和发送方使用完全相同的随机数生成顺序否则重建出的矩阵是不同的。如果中间插入了其它随机数调用后面的矩阵全都会变。6.5 彩色图像和视频怎么扩展上面所有的讨论都基于灰度图。如果是彩色图像最省事的做法是分别处理R、G、B三个通道每个通道独立压缩、加密、解密最后合并。代价是计算时间增加三倍压缩率不变。更高效的做法是先把彩色图像从RGB转到YCbCr对亮度分量保持较高测量率对两个色度分量使用较低测量率利用人眼对色度细节不敏感的特性做视觉优化。视频处理可以逐帧执行这套流程但帧间解相关的效率不如专门的关键帧加密方案。如果你真的要把这个方案落到视频系统里建议先对视频做运动估计只对关键帧执行DRPECS加密非关键帧用轻量级置乱即可。我个人在实际调试中的体会是这类算法最花时间的不是加密环节而是解密端的重构参数调优。相位掩码和测量矩阵的随机性保证了安全性但性能上限往往由稀疏基和OMP迭代次数决定。你先在一张小图上把链路完整跑通再逐步放大尺寸、压低测量率会比一上来就追求大图和极低压缩率舒服得多。