基于密钥控制测量矩阵的图像压缩加密混合算法Matlab实现
最近在做一个图像安全传输的小项目需要在压缩数据的同时保证内容不可见翻了几篇文献发现大家都绕不开一个思路把压缩感知里的测量矩阵当成密钥载体也就是俗称的“密钥控制测量矩阵”。这个思路很讨巧——测量矩阵本来就必须由收发双方共享才能重建那干脆让密钥来决定测量矩阵长什么样于是压缩过程本身就变成了加密过程一套流程同时完成两件事。这篇文章把我用Matlab复现这类“图像压缩加密混合算法”的完整过程记录下来包括密钥怎么一步步变成测量矩阵、压缩端和解压端怎么组织代码、OMP重建要做到什么程度、以及我踩过的几个版本兼容和参数调优的坑。适合正在做图像加密或压缩感知课题的本科生、研究生也适合想快速上手Matlab复现论文算法的工程师。1. 为什么说“测量矩阵”才是混合算法的命门1.1 压缩感知的采样模型天然就是一种加密投影先回顾一下压缩感知Compressed Sensing, CS的核心模型。假设一张图像对应的列向量是 \(x \in R^N\)但 \(x\) 本身往往不是稀疏的需要在某个变换域里才能稀疏表示比如小波域或DCT域[ x \Psi s ]其中 \(\Psi\) 是稀疏基矩阵\(s\) 是稀疏系数向量大部分分量接近0。然后我们用一个测量矩阵 \(\Phi \in R^{M \times N}\) 去线性投影[ y \Phi x \Phi \Psi s ]这里 \(y\) 就是测量值维度从 \(N\) 降到了 \(M\)且 \(M \ll N\)。重建的时候要从 \(y\) 反推出 \(s\)理论上这是个欠定问题但因为 \(s\) 稀疏可以通过OMP、CoSaMP或L1范数优化等算法恢复。关键点来了整个环节里测量矩阵 \(\Phi\) 的地位极其特殊。接收方如果不知道 \(\Phi\) 具体是什么那么拿到 \(y\) 也几乎不可能还原图像。这不就是一个现成的加密要素吗打个比方压缩感知相当于把一张清晰的图投影成几张混叠的剪影不知道投影的方式测量矩阵你就没法从剪影还原原图。于是让密钥控制测量矩阵的生成就是把这套线性投影机制直接变成密码学意义上的加密变换。当然这里有一个必须说清楚的前提经典的高斯随机测量矩阵本身并不是为加密而设计的它只关心满足RIP性质、保证重建质量。当它承担加密职责时需要考虑密钥空间够不够大、测量矩阵是否会被逆推出来。这就是为什么纯CS加密往往强度有限需要通过混合算法在测量矩阵之外再加置乱、扩散等操作。1.2 密钥控制测量矩阵的两种主流构造流派我复现过的文献里密钥控制测量矩阵的构造基本可以分成两派。第一派伪随机数生成器种子控制。用密钥数字作为随机数种子直接生成高斯随机测量矩阵。比如Matlab里用 rng(key) 或 RandStream 生成 \(\Phi randn(M,N)/\sqrt{M}\)。这个方案实现成本最低重建质量也最接近理论值因为高斯随机矩阵的RIP性质非常好。缺点是安全性偏弱随机数生成算法是公开的攻击者一旦猜到种子格式就可能暴力搜索而且测量矩阵的非零元素数量是M×N密钥空间和存储开销都比较可观。第二派混沌映射控制。用Logistic映射这样的混沌系统生成确定性序列再reshape成测量矩阵。比如[ x_{n1} \mu \cdot x_n \cdot (1-x_n), \quad \mu \in (3.57, 4] ]初值 \(x_0\) 和参数 \(\mu\) 就构成密钥混沌系统对初值极度敏感差这么 \(10^{-15}\)迭代出来的序列就完全不同。这一派的优点是密钥形式灵活可以做到“一次一密”而且确定性的混沌序列在收发两端都能很方便地复现。缺点是需要小心处理暂态效应生成序列时要先丢弃前面若干点否则测绘矩阵与初值的关联太明显容易泄漏信息。不管是哪一派落到这篇论文标题里的“混合算法”通常都是在密钥控制测量矩阵的基础上再对测量值或系数做一次置乱/扩散相当于给压缩感知加了第二道锁。测量矩阵承担“压缩采样加密”后处理承担“混淆扩散”两者结合才是完整的新型图像压缩加密混合算法。2. Matlab实现的核心骨架密钥如何一步步变成测量矩阵再回到图像2.1 整体流程怎么组织代码层面我习惯把整个算法拆成发送端和接收端两个子流程。发送端压缩加密读入图像灰度化转成double类型。选定稀疏基 \(\Psi\)我用的是DWT小波基选sym4。计算稀疏系数 \(s \Psi^T x\)。用密钥生成测量矩阵 \(\Phi\)。测量得到 \(y \Phi x\)。对 \(y\) 做量化必要时再做置乱得到最终密文。接收端解密重建用同一个密钥恢复出 \(\Phi\)。对密文做逆置乱、反量化得到 \(y\)。用OMP算法从 \(y\) 和 \(\Phi\Psi\) 中恢复稀疏系数 \(\hat{s}\)。反变换 \(\hat{x} \Psi \hat{s}\)得到重建图像。主脚本结构大概长这样%% 参数区 imgPath lena256.png; blockSize 64; % 分块尺寸 Mratio 0.5; % 测量率 M/N sparseK round(blockSize^2 * 0.1); % 稀疏度估计 %% 读取与转换 img imread(imgPath); if size(img,3) 3 img rgb2gray(img); end imgD double(img); %% 压缩加密 [compressed, Phi, keyInfo] cs_encrypt_compress(imgD, blockSize, Mratio, sparseK); %% 解密重建 imgRec cs_decrypt_decompress(compressed, Phi, keyInfo, size(imgD));2.2 密钥控制的高斯测量矩阵生成代码这里我给出我实际用过的方案先用 RandStream 保证不同Matlab版本行为一致同时把任意长度的密钥字符串映射成随机数种子。function Phi keyed_gaussian_measurement_matrix(key, M, N) % key: 字符串或者数字作为随机种子 if ischar(key) key string2hash(key); % 用djb2之类的hash把字符串转成uint32 end s RandStream(mt19937ar, Seed, key); Phi randn(s, M, N) / sqrt(M); end为什么用 RandStream 而不是老的 rand(seed, key)后一种方式在MATLAB R2011a之后就推荐不要再用了到了较新的版本里 rand(state) 这种老语法已经完全不工作。我这里用 RandStream(mt19937ar, Seed, key) 是mt19937这个梅森旋转算法的稳定入口只要版本支持生成出来的序列完全一致。密钥如果是个字符串我通常会先做个哈希将任意长度输入变为uint32值因为RandStream的Seed参数接受的范围有限直接塞一个很长的字符串会报错。2.3 用Logistic混沌序列构造测量矩阵的细节混沌序列构造测量矩阵我常用的做法是function Phi logistic_measurement_matrix(x0, mu, M, N) burnIn 500; % 丢弃前500个迭代点避免暂态效应 totalLen burnIn M * N; x zeros(totalLen, 1); x(1) x0; for i 2:totalLen x(i) mu * x(i-1) * (1 - x(i-1)); end seq x(burnIn1:end); Phi reshape(seq, M, N); Phi Phi * sqrt(12 / M); % 归一化让方差落在合理范围 end这里有两个工程细节值得展开。第一为什么丢弃前500个点Logistic映射从初值开始迭代的头几十个点与初值高度相关如果不丢弃攻击者可能通过分析测量矩阵的反推来逼近密钥。丢弃足够长的迭代序列后混沌轨道已经“满血”进入混合状态这时候提取的序列统计特性更接近随机。论文里也有类似操作只是代码注释经常省略自己实现时千万别省。第二归一化的问题。直接reshape出的混沌序列均值不为0方差也不固定如果直接当测量矩阵用测量值的动态范围会漂移影响重建。我习惯把 \(\Phi\) 缩放成类似高斯随机矩阵的量纲即每个元素方差接近 \(1/M\)。上面的 sqrt(12/M) 是因为Logistic混沌序列在均匀分布假设下方差约1/12乘以sqrt(12/M)后整体方差近似1/M跟高斯随机矩阵处于同一水平。2.4 重建端OMP核心流程与Matlab函数要点测量值怎么重建是整个算法能不能用的关键。OMP正交匹配追踪是最容易理解和实现的算法之一核心思路就是贪心每次从字典里找出与当前残差最相关的原子用最小二乘更新系数再更新残差。function s_hat OMP(y, A, K) % y: 测量值 (M x 1) % A: 感知矩阵 Phi*Psi (M x N) % K: 稀疏度 [M, N] size(A); s_hat zeros(N, 1); r y; idx zeros(K, 1); for t 1:K corr A * r; [~, i] max(abs(corr)); idx(t) i; A_t A(:, idx(1:t)); s_t pinv(A_t) * y; %% 用pinv而不是A_t\y数值更稳 r y - A_t * s_t; end s_hat(idx(1:t)) s_t; end这段代码里有一个实际使用中很重要的细节pinv(A_t) * y 而不是 A_t \ y。当A_t的病态程度较高时左除可能给出很大的系数震荡pinv通过SVD求最小二乘解会更稳定。代价是计算量稍大但对于图像块尺寸在32~64范围内的应用场景这点开销完全可以接受。OMP的K值也就是稀疏度不能拍脑袋乱给。我一般会在压缩前先对稀疏系数 \(s\) 做一次能量分析比如用 wenergy 看小波系数每一层的能量占比然后选择一个能保留95%能量的系数个数作为K。如果K给得太小重建图像会明显模糊K给得太大又容易把噪声当成有效原子重建反而出现颗粒状伪影。2.5 大图的块效应与分块策略这个算法直接拿整张图建测量矩阵会遇到一个现实问题测量矩阵是 \(M \times N\) 的稠密矩阵。假设图是512×512N是262144就算M/N取0.3\(\Phi\) 也有几十万行乘二十多万列Matlab内存直接爆掉。所以我在工程实现里都是分块处理常见做法是把图像切成64×64的小块每个块单独做CS测量和重建。块尺寸测量率M/N单块测量矩阵大小重建质量适用场景32×320.5512×1024中等快速演示64×640.52048×4096较好默认推荐128×1280.58192×16384内存压力大追求细节时使用分块造成的副作用是块与块之间的边界可能出现接缝尤其重建质量差的时候特别明显。缓解办法有两个一是对图像先做重叠分块块与块之间重叠几个像素重建后只取中间区域二是重建后加一步去块效应滤波例如用 imgaussfilt 做轻微平滑。当然如果分块尺寸足够大64以上接缝通常不明显可以不用额外处理。3. 仿真结果怎么看PSNR、SSIM、相关系数与安全性测试3.1 图像质量的基本指标复现混合算法的第一步是要确认压缩重建质量没有因为加了密钥控制而明显变差。我常用的指标是PSNR和SSIM。PSNR的计算很简单mse mean((imgD(:) - imgRec(:)).^2); psnr 10 * log10(255^2 / mse);SSIM用一个Matlab内置函数 ssim(imgRec, imgD) 就能算。一般来说PSNR高于30dB、SSIM高于0.9肉眼看起来就相当不错了。我在256×256测试图上实测下来测量率M/N0.5时采用keyed高斯测量矩阵配合DWT稀疏基、OMP重建PSNR大致在33~36dB之间M/N降到0.25PSNR大约在28dB左右。混沌Logistic测量矩阵的重建质量略低于高斯随机矩阵差1~2dB但安全性特征更明显密钥空间更大。3.2 密钥敏感性测试错一个比特全图崩溃加密算法最基础的检验就是密钥敏感性。做这个测试的方式很简单用原始密钥重建一幅图再用一个仅相差 \(10^{-15}\) 的密钥重建另一幅图然后对比两幅图。我实测的结果是正确密钥重建时PSNR约34dB错误密钥重建时PSNR直接掉到9dB以下肉眼看起来就是纯噪声原图信息几乎完全泄漏不出来。相关系数 corr2(imgD, imgRecWrong) 也基本接近0。这个结果说明两点混沌映射对初值的敏感性确实完整体现在了图像重建上测量矩阵作为密钥载体的设计方向是成立的。3.3 对噪声和量化误差的鲁棒性混合算法不能只在无噪理想环境里跑实际信道多少会有干扰。我做过的鲁棒性测试大致分两类。一类是测量值经过信道时叠加高斯噪声比如测量值SNR在20dB、15dB、10dB时重建质量的变化。压缩感知本身对噪声有一定的抑制能力测量值稍微带噪重建图像PSNR的下降幅度通常小于噪声幅度的下降幅度。另一类是量化误差的影响。测量值 \(y\) 是浮点数传输出去必须量化成有限bit。我对比过8bit量化和4bit量化对重建质量的影响8bit基本无损4bit时PSNR会下降2~4dB图像仍有轮廓但细节丢失严重。如果论文里说自己的算法支持量化后重建最好在4bit到8bit之间都测一组数据这样审稿人或者读者会比较信服。4. 我在复现中踩过的坑参数、工具箱与版本兼容4.1 rand(state) 的版本兼容坑这个坑我印象太深了。网上能找到的很多老代码用 rand(seed, key) 或者 rand(state, key) 来控制随机测量矩阵但这套语法在新版MATLAB里要么报错要么行为跟旧版完全不一样。我在R2020a上跑一个老脚本直接提示 state is an invalid argument。解决方法是用 rng(key) 或者 RandStream。如果你要复现一篇论文里面说“用高斯随机测量矩阵”的实验我建议直接用rng(key); Phi randn(M, N) / sqrt(M);而如果是做算法比较多组实验之间需要完全相同的测量矩阵那就应该用 RandStream 显式保存流对象而不是依赖全局随机数状态。4.2 OMP迭代次数和系数个数不匹配图像会“发灰”有段时间我的重建图像总是灰蒙蒙的像蒙了一层雾。排查下来发现是稀疏度K设得偏大OMP强行选了过多原子去拟合导致引入了一些本不该出现的重建噪声。这个问题的解决思路是不要静态设置K而是根据块内能量动态调整。我写了一个简单函数计算小波系数按降序累加的能量占比保留95%能量需要多少个系数就取这个数量作为K。这样每次重建都会自适应图像质量更稳定。4.3 uint8和double类型问题Matlab里图像处理最常见的类型坑在CS混合算法里同样躲不开。imread 读进来的是uint8直接做矩阵运算时乘法溢出、减法负数截断都会出问题。我的习惯是第一步就 double(img)而且乘除法都放在double域里做最后显示或保存再转回uint8。DWT变换也跟数据类型有关系。wavedec2 对double类型的输入处理起来最直接uint8输入会被内部转换但显示出来的量化特性不一样。源代码里如果混用了uint8和double经常会出现PSNR反复横跳的情况这种情况多半不是算法问题是类型转换位置的问题。4.4 密钥管理怎么做才不容易翻车还有一个容易被忽视的问题密钥本身怎么从字符串变成测量矩阵参数。我遇到过这样一个场景——密钥写成字符串直接拿来当Logistic映射的初值或者当随机种子结果不同长度的密钥互相之间没有形成足够大的密钥空间暴力搜索难度不够。稳妥的做法是把密钥先做一次哈希比如djb2或SHA-1摘要把hash结果映射到一个浮点初值或一组整数种子再进入测量矩阵生成流程。策略上属于简单实用但能有效拉开不同密钥之间的距离。5. 一些个人体会这类论文复现下来我的体会是学术代码容易跑通但算得漂亮需要耐心做参数微调。密钥控制测量矩阵是整个混合算法里“压缩和加密握手”的地方设计得好密钥空间、重建质量和计算效率三者能兼顾设计得粗糙加密强度上去了重建图像却糊成一片或者重建质量挺好但密钥空间太小经不起推敲。如果时间有限我建议先从经典高斯随机测量矩阵加OMP重建做起跑通全流程再换成混沌密钥控制安全性测试务必做一组错误密钥重建对比这是整个算法成立的最直接证据。重建端如果追求更高图像质量可以在输出后加一个BM3D或维纳滤波做后处理PSNR往往能再往上提0.5~1dB又不需要改动加密端。这些小技巧都是自己反复调代码换来的。