双随机相位编码结合压缩感知的图像加密Matlab实现解析

做图像加密项目时，我最早接触双随机相位编码（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级，(P=46) 时，两个掩码的组合空间也已经到 (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加密”。具体流程：

1. 读入灰度图像，统一尺寸，转成double。
2. 对图像做二维DCT，得到稀疏系数矩阵。
3. 用高斯测量矩阵对系数矩阵的向量形式做投影，得到测量值。
4. 为了让测量值能进入二维DRPE链路，把测量值折叠成一个方阵。
5. 生成两块随机相位掩码，对折叠后的方阵做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 解密端的流程与重构细节

解密端流程：

1. 对收到的复数密文C做逆DRPE，得到含数值误差的测量值矩阵。
2. 拉直并取实部，因为DCT系数和测量值本来就是实数，DRPE逆运算后的虚部只是浮点计算误差。
3. 用OMP从测量值恢复DCT系数向量。
4. 把系数向量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迭代次数优化后的典型结果：

数值只是趋势参考，不同图像的绝对数值会有差异。但如果测量率已经到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}))，再把实际测量矩阵的行数定为 (m=P^2)。第二种的好处是没有任何补零操作，测量率和预期值之间只差一个很小的取整误差。我建议优先用第二种。

6.3 逆DRPE的conj与取实部

逆DRPE里使用conj(R2)和conj(R1)是最常见的坑。单位复数相位掩码的共轭就是倒数，这是由相位掩码的构造决定的。如果你生成掩码时用了exp(1j * rand * 2 * pi)，那逆操作必须用conj。如果用的是任意复数矩阵，那就不能简单共轭，而是要取逆矩阵或者求伪逆，但那样通常不叫DRPE。

解密后Y_hat会包含一层数值量级的虚部，这是因为FFT/IFFT的浮点运算不是完全精确的。直接把虚部丢掉不会丢失有效信息，因为明文侧测量值本身就是实数。我在测试中发现，如果不加real，OMP会因为复数虚部干扰导致恢复质量下降，所以解密代码里一定记得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，对亮度分量保持较高测量率，对两个色度分量使用较低测量率，利用人眼对色度细节不敏感的特性做视觉优化。

视频处理可以逐帧执行这套流程，但帧间解相关的效率不如专门的关键帧加密方案。如果你真的要把这个方案落到视频系统里，建议先对视频做运动估计，只对关键帧执行DRPE+CS加密，非关键帧用轻量级置乱即可。

我个人在实际调试中的体会是，这类算法最花时间的不是加密环节，而是解密端的重构参数调优。相位掩码和测量矩阵的随机性保证了安全性，但性能上限往往由稀疏基和OMP迭代次数决定。你先在一张小图上把链路完整跑通，再逐步放大尺寸、压低测量率，会比一上来就追求大图和极低压缩率舒服得多。