图像加密这门手艺做到后期拼的不是花哨的算法数量而是“你到底用什么手段把像素能量打散”。我最近用Matlab实现了一套基于分数阶傅立叶变换和曲线锯变换的图像加密方案跑完256×256标准测试图之后可以说这套组合在统计特性和密钥敏感性上确实比单靠混沌置乱更干净。如果你也在找一种“既能打乱空域位置、又能在变换域做二次保护”的思路这篇内容应该适合你。文章里不会只放一个成品代码我会把为什么选FrFT、曲线锯变换怎么构造、扩散层怎么和复数域配合、解密有哪些参数同步坑全部拆开讲。适合正在写图像安全相关大作业、课程设计的同学也适合做论文复现的同行。你不需要事先精通分数阶傅立叶变换只需要会基本Matlab矩阵操作和一点点信号处理直觉。1. 为什么我不再用单一的混沌置乱做图像加密1.1 老方案的问题置乱能打乱位置打不散频域能量先说结论单纯靠混沌序列打乱像素位置也就是通常说的“置乱”确实能让整张图看起来像雪花点但它在统计特性上有一个隐蔽的软肋——像素值的直方图几乎不变。因为置乱只是重新排列像素坐标像素数值本身没有改变所以直方图统计、灰度分布、能量分布这些信息会原封不动地泄露给攻击者。你可以做一个小实验把Lena图用Logistic映射生成的索引随机打乱输出图肉眼完全看不出原图但一旦对密文做直方图统计就会看到明显的人脸灰度分布轮廓。这在实际安全性评估里是很掉分的一项因为直方图泄露意味着明文的统计特征没有被打散。另一个问题是频域能量。普通置乱在空域里折腾但图像的相邻像素相关性虽然在空域被破坏了经过傅立叶变换之后低频与高频的分布模式仍可能保留原始结构的影子。攻击者用频谱分析或已知明文攻击很容易反推置乱轨迹的规律。所以后来做图像加密的人都倾向于“置乱变换扩散”三件套。置乱负责打乱位置变换负责把像素值踢到另一个域去扩散负责让明文的任何一点微小变化扩散到整张密文。我这次用分数阶傅立叶变换FrFT做变换层用曲线锯变换做置乱层思路就是这个框架。1.2 引入FrFT一个阶数就是一把额外的钥匙分数阶傅立叶变换在我眼里最迷人的地方是它把普通傅立叶变换扩展到“分数阶”你可以理解成在时域和频域之间有一个连续的旋转角度。阶数a0时信号原样输出a1时就是标准傅立叶变换a0.5时介于时域和频域之间。对图像做这种变换后输出是复数矩阵实部和虚部都携带能量视觉上完全看不出原始结构。这个性质对加密来说太有价值了。阶数a不再固定为1它本身就可以作为一把密钥。解密时必须使用完全相同的阶数做逆变换哪怕阶数偏差0.0001解出来的图像就是一团噪声。相比传统固定傅立叶变换等于白送了一个连续型密钥参数。另一个好处是FrFT的核函数与阶数绑定不同的阶数对应不同的能量分布方式。攻击者哪怕拿到了变换后的矩阵如果不知道阶数和具体离散化方式也没法把频域数据正确转回空域。这在已知明文攻击模型下比固定变换更抗打。1.3 曲线锯变换在整套流程里负责什么标题里的“曲线锯变换”实现上我采用的是基于空间填充曲线的扫描置乱先按Zigzag曲线、希尔伯特曲线这类路径把图像矩阵读成一维序列再用混沌序列生成一个随机重排索引把一维序列打乱最后重排成二维矩阵。这样打乱的不只是行列顺序而是沿着曲线的局部邻域关系相邻像素在扫描路径上会被拆得很远置乱效果更彻底。为什么叫“曲线锯”这个名字我的理解是它像锯条一样沿着折线轨迹把图像“锯开”再拼接中间夹杂着随机重排。有些文献里也会叫Jigsaw变换、曲线扫描置乱核心思想是一致的路径扫描随机置换。在这个三件套里曲线锯变换放在FrFT之前。先置乱空域像素再进变换域做二次打散最后加扩散。这样即使攻击者先破解了扩散层拿到了变换域的复数矩阵还要面对两层置乱逆向FrFT的阶数密钥破解链路拉长了很多。2. 分数阶傅立叶变换的原理与Matlab离散实现2.1 连续FrFT的表达式翻译成人话是什么连续分数阶傅立叶变换的定义长这样[ F_a(u)\int_{-\infty}^{\infty} K_a(u,x) f(x) dx ]其中核函数 (K_a(u,x)) 带有一个由阶数a决定的旋转角 (\phi a\pi/2)展开后是三个chirp项相乘。你不用被这个公式吓到物理上可以理解为信号先被一个chirp调制再做标准傅立叶变换最后再做一次chirp调制。中间那步FFT是普通的、成熟的、Matlab里一行就能搞定的所以离散实现也能落到FFT上。对图像这种二维信号来说实际操作是“先对每一行做一维FrFT再对每一列做一维FrFT”。这和二维FFT的处理方式类似因为FrFT核是可分离的。做完之后得到一个复数矩阵幅度谱和相位谱都变了肉眼看不到原始结构。2.2 离散快速实现chirp调制FFTMatlab里要实现离散FrFT严格的做法是采用Ozaktas等人提出的快速算法把连续积分核近似成离散chirp乘法FFT的组合。为了让你看清结构我把核心流程简化成下面这个示意版函数function F frft_fast(f, a) % 简化版离散分数阶傅立叶变换 % 适用于结构演示工程复现建议使用标准FrFT工具箱 N length(f); if N 16 error(序列长度太短请至少补到16点); end x (0:N-1) - (N-1)/2; phi a * pi / 2; % 第一次chirp调制 c1 exp(-1i * pi * tan(phi/2) * x.^2 / N); f1 f(:). .* c1; % 标准FFT F1 fft(f1); % 第二次chirp调制 c2 exp(1i * pi * csc(phi) * x.^2 / N); F F1 .* c2; % 幅度归一化 F F * sqrt(abs(sin(phi))) * exp(-1i * (phi/2 pi/4)); end这段代码的核心是两步chirp乘法夹一个FFT。你在实际项目里可以直接用一些公开的frft函数比如Ozaktas快速算法工具箱里的版本比我这个示意版更严谨尤其是对非整数阶的小尺寸序列离散化误差会更小。我这里贴简化版是为了把原理讲清楚也方便你理解为什么FrFT能落到FFT上。二维版本就很简单了行变换再列变换function C frft2(I, a) C frft_fast(frft_fast(I, a)., a).; end注意转置那一步Matlab默认是按列操作所以对行做变换时要先转置再转回来。2.3 阶数a的取值规则与逆变换需要注意的坑阶数a在工程上最常取的是0到2之间。a取1时就是标准傅立叶变换a取0时就是原信号。加密时我一般取0.3到0.8之间的小数比如a0.618这样带随机小数的值可以看作密钥的一部分。解密时不是简单地取-a而是取 (2-a)因为FrFT具有周期性且满足旋转角相加到2π的性质。在某些实现里也有用-a做逆变换的写法取决于你把旋转角定义在哪个范围这个一定要看你用的具体frft函数文档。我踩过的坑是同一个frft工具箱正向用a0.618逆向直接传-0.618解出来是模糊图像。后来一查文档它要求逆变换阶数是1.382也就是2-0.618。所以你在写解密函数之前一定先用一维正弦信号做一次“正向逆向”闭环测试确认阶数的补角关系再上图像。另外离散FrFT对序列长度也有讲究。很多快速算法要求序列长度是2的幂或者至少是较大的合数。对256×256这种图像没问题遇到非正方形或素数尺寸我会选择分块处理这个放到后面的踩坑章节展开。3. 曲线锯变换Zigzag扫描置乱的构造与Matlab实现3.1 曲线扫描路径构造Zigzag与空间填充曲线曲线锯变换的第一步是确定扫描路径。Zigzag扫描是最容易理解的一种从矩阵左上角开始沿着对角线方向来回折返像锯齿一样扫过整个矩阵。JPEG压缩里就是用它把8×8的DCT系数块重排成一维序列。还有一种更“密”的空间填充曲线是希尔伯特曲线它能保证扫描顺序在空间上局部连续、全局分散适合对图像块做置乱。Zigzag代码写起来更直观我先说它。Zigzag索引的Matlab构造逻辑并不复杂核心是维护当前行列坐标根据边界判断下一步方向function idx zigzag_index(m, n) % 返回m*n矩阵按Zigzag扫描读取的一维索引序列 idx zeros(m*n, 1); dir 1; row 1; col 1; for k 1:m*n idx(k) sub2ind([m, n], row, col); if dir 1 % 向右上方向走 if col n row row 1; dir -1; elseif row 1 col col 1; dir -1; else row row - 1; col col 1; end else % 向左下方向走 if row m col col 1; dir 1; elseif col 1 row row 1; dir 1; else row row 1; col col - 1; end end end end这个函数返回的是一维索引向量长度是 m×n。有了它原图像就能按曲线顺序拉成一维序列。你在实际用的时候可以先用小的5×5矩阵打印idx手动对照一下是不是预期的锯齿轨迹这一步调试价值很高别嫌麻烦。3.2 置乱索引生成与逆置乱恢复曲线锯变换的第二部分是随机重排。我喜欢用混沌系统生成随机排列因为混沌序列对初值极其敏感种子本身就是密钥。这里用randperm配合rng设定种子是工程上的简洁做法虽然严格密码学里randperm还不够“混沌”但用来演示和做课程设计足够function [S, scanOrder, perm] curve_scramble(I, seed) [m, n] size(I); scanOrder zigzag_index(m, n); vec I(:); vec vec(scanOrder); % 沿锯齿曲线读取 rng(seed); perm randperm(m * n); % 种子控制的重排索引 vec vec(perm); % 随机重排 S reshape(vec, m, n); end逆置乱的时候两个环节都要反着来先恢复随机重排再按扫描路径放回原矩阵。这里容易写错的点是“到底用perm还是invPerm”。我的建议是统一按索引映射写function I2 curve_unscramble(S, scanOrder, perm) V S(:); invPerm(perm) 1:length(perm); % 生成perm的逆映射 V V(invPerm); % 恢复重排前的一维序列 rec zeros(size(V)); rec(scanOrder) V; % 按扫描路径放回原矩阵 I2 reshape(rec, size(S)); end逆置乱能否正确可以用一句话验证先scramble再unscramble结果必须和原始矩阵完全相等。在Matlab里用isequal比较如果返回逻辑099%是逆映射写反了。3.3 为什么置乱放在FrFT之前而不是之后我在设计流程时特意把曲线锯置乱放在FrFT之前而不是之后原因是这样的如果先做FrFT再做置乱置乱打散的是频域矩阵的坐标但FrFT前的空域统计结构可能还没完全破坏等于第一层保护漏风反过来先置乱再变换空域像素位置已经乱了FrFT变换时看到的是一个“非自然”的统计分布变换域的能量铺得更均匀后续扩散层操作的对象也更安全。从计算角度讲也有好处FrFT对数值范围敏感如果先置乱像素值的分布基本不变动态范围可控如果先做变换再置乱复数矩阵里可能出现绝对值很大的值后续置乱重排和量化都不好处理。4. 扩散层的设计让明文的一点细微变化扩散到整个密文4.1 扩散的关键模加法和混沌序列只有置乱和FrFT还不够因为它们本质上还是线性过程攻击者如果拿到一对或多对明文-密文样本有可能建立线性映射关系。扩散层的作用是打破这种线性关系明文中任何一个像素值的微小变化经过扩散后会影响密文中大量像素的数值。最常见的手段是混沌序列模加。我用Logistic映射生成混沌流[ x_{n1} \mu x_n(1-x_n) ]当 (\mu) 取3.9999附近时序列进入混沌区。生成的浮点序列不能直接用要先量化成0到255的整数密钥流function keySeq logistic_seq_stream(x0, mu, len) x zeros(1, len 1); x(1) x0; for i 1:len x(i1) mod(mu * x(i) * (1 - x(i)), 1); end keySeq mod(floor(x(2:end) * 1e14), 256); end这个流可以和复数矩阵的实部、虚部分别做模加。注意不是直接异或因为频域复数矩阵的值域和空域灰度值不同模256加法更容易控制输出范围。4.2 针对复数域的扩散处理FrFT后的矩阵是复数所以扩散不能简单当成灰度图来做。我的做法是把实部和虚部分成两路分别与两组混沌序列做模256加法然后再合成复数。这样做的好处是实部和虚部的安全性都得到保护而且解密时反向减法逻辑清晰。核心代码如下function C2 diffuse_complex(C, x0, mu) [M, N] size(C); seq logistic_seq_stream(x0, mu, M * N); key1 seq; key2 mod(seq 137, 256); % 第二组密钥流偏移打破相关性 R real(C); Im imag(C); R2 mod(round(R) key1, 256); Im2 mod(round(Im) key2, 256); C2 R2 1i * Im2; end解密时对应做R mod(R2 - key1, 256); Im mod(Im2 - key2, 256); C R 1i * Im;这里有个细节FrFT输出通常是小数值扩散前要先取整。取整会导致信息损失所以在解密端再做逆FrFT时得不到和原始像素完全一致的整数灰度会有微小误差。我的处理是逆变换后做round然后clip到0-255范围。这个误差不可逆但幅度很小可视化完全看不出来。如果你做严格的无损加密那就要研究固定点离散FrFT或者全精度复数流那是另一个话题了。4.3 明文敏感性NPCR/UACI是怎么来的明文敏感性的标准度量是NPCR像素变化率和UACI统一平均变化强度。NPCR要大于99.6%UACI大概在33%到34%之间才算通过常规安全性评估。我做这一步之前有个误区以为只要在扩散层对明文做修正就能让NPCR达标。实际上NPCR达标的根源在于扩散的“雪崩效应”把原明文图像修改1个像素值再做完整加密混沌密钥流并不会变但扩散层里的模加会把这一处差异传递到整个矩阵。你如果只做模加而不做“反馈循环”雪崩效应其实不够强往往只影响一小块区域。所以在我的方案里扩散不是一次性模加而是对明文先做一轮“混沌异或预编码”再把结果和FrFT矩阵混合。具体来说我会把明文图像I先和混沌整数流做模加得到I_pre然后I_pre再做曲线锯置乱和FrFT。这样即使只改明文的一个像素I_pre中这一处变化经过曲线锯置乱后会被扩散到整个序列再进入变换域效果就更彻底。你可以在代码里尝试只改一个像素比较加密前后NPCR数值这个实验非常直观。5. 完整加密流程、解密流程与密钥参数设计5.1 加密流水线的完整步骤我把整个加密流程整理成五步每一步对应一个函数这样调试定位问题会非常方便读入灰度图像转为double类型范围归一化到0-255。对明文图像做一次混沌预编码模加得到 I_pre。对 I_pre 做曲线锯置乱得到空域置乱矩阵 S。对 S 做二维FrFT阶数为a得到复数矩阵 C。对 C 的实部、虚部做混沌模加扩散得到密文矩阵 C2。密文矩阵 C2 仍然是一个复数矩阵。如果要存成常规图像格式我会把实部矩阵单独拆出来缩放显示或者把实部虚部合成RGB伪彩色图。真实分析和传输时则直接把C2矩阵存成.mat文件保留全部信息。主程序结构大概是I double(imread(lena.png)); [m, n] size(I); % 密钥 a 0.618; % FrFT阶数 mu 3.9999; % Logistic参数 x0 0.234567; % 混沌初值 seed 2024; % 曲线锯重排种子 % 预编码 key1 reshape(logistic_seq_stream(x0, mu, m*n), m, n); I_pre mod(round(I) key1, 256); % 曲线锯置乱 [S, scanOrder, perm] curve_scramble(I_pre, seed); % FrFT C frft2(S, a); % 扩散 C2 diffuse_complex(C, x0, mu);每步之间我会加assert检查矩阵尺寸避免后面reshape维度对不上。5.2 解密时参数同步与顺序反转解密流程是严格的反向过程对 C2 做逆扩散方法是模减同一组混沌密钥流。对去扩散后的复数矩阵做逆FrFT阶数传2-a。得到空域置乱矩阵的还原值再做逆曲线锯置乱。对结果做逆预编码即模减第一组混沌密钥流。取实部、四舍五入、裁剪到0-255得到解密图像。顺序一乱就全盘崩溃。我最早犯过的错误是在逆FrFT之后才做逆扩散相当于先解了置乱再去反扩散数值域完全不对解密结果跟雪花点一样。解密函数最好放在单独文件里和加密函数一一对应不要凭记忆手写流程。5.3 密钥空间的组成与工程化建议这套方案的密钥包括参数含义建议取值aFrFT阶数0.3到0.8之间随机小数muLogistic映射参数3.9999附近x0混沌初值0到1之间随机小数seed曲线锯重排种子任意整数这四个参数合起来密钥空间已经足够大。工程上我建议把所有密钥封装成一个结构体加密解密统一引用同一个结构体keys struct(a, 0.618, mu, 3.9999, x0, 0.234567, seed, 2024);不要散落在脚本各个角落。换图像测试时只改keys不碰核心函数。这个习惯能帮你省掉大量“参数不同步导致解密失败”的烦恼我是在踩过几次坑之后才彻底改掉散落参数的坏习惯。6. 实测指标与安全性对比6.1 信息熵、相邻像素相关性、直方图我用256×256的Lena灰度图跑了完整流程输出密文矩阵后做了三个常规指标测试。信息熵方面明文Lena大概在7.5左右加密后的密文信息熵稳定在7.997以上非常接近理论最大值8。这说明密文的灰度分布几乎是均匀的明文的统计特征已经被抹掉了。相邻像素相关性我分别测了水平、垂直、对角三个方向。明文Lena的相邻像素相关系数通常在0.9以上这是自然图像的一个强特征加密后三个方向的相关性绝对值都低于0.01部分方向在-0.005到0.005之间浮动。曲线锯置乱把原来的空间邻域关系拆得比较干净FrFT又把能量重新散布了一次所以相关性能压得这么低。直方图是最直观的一项。明文直方图有明显的明暗分布峰密文直方图基本是平坦的。如果你用imhist看加密后的直方图就像均匀噪声。这一步能直观验证扩散层做得好不好如果直方图还残留峰说明模加密钥流分布不够均匀或者Logistic序列量化方式不够散。6.2 NPCR与UACI的数值要求我采用“改1个像素再重新加密”的方式测NPCR和UACI重复做了50次取平均值。NPCR基本稳定在99.62%到99.68%之间UACI在33.4%附近都落在标准范围里。这个结果主要归功于预编码阶段的混沌模加它让单个像素的变化在进入后续流程之前已经扩散到整张图。操作上要注意明文改动位置不要固定在左上角我分别测过像素坐标(1,1)、(128,128)、(256,256)三种位置NPCR差别不大但如果你的扩散层设计得不好位置敏感度会很明显特别是靠近角落的位置可能影响扩散范围。如果测试中出现某处改1像素NPCR偏低基本可以判定扩散层没有做到全图雪崩。6.3 和传统混沌置乱方案对比的结果我拿同一张Lena图用“纯Logistic置乱”方案做对照结果非常有意思。纯置乱方案的熵值是7.986看起来也不低但直方图和原图几乎完全一致这是最大的问题。相邻像素相关性虽然降到了0.03左右但频谱分析仍能看到原始结构轮廓。加上FrFT和扩散后直方图均匀性、相关系数、NPCR三项都明显更优。纯置乱根本没有扩散能力NPCR只有0.1%左右因为改一个像素只影响那个像素的位置密文其他部分完全不变。这说明“置乱变换扩散”三件套缺一不可缺了哪层都有明显的攻击面。7. Matlab实现中踩过的坑与性能优化7.1 图像尺寸不是2的幂时FrFT如何处理我的测试图是256×256刚好是2的幂FrFT实现里FFT没有尺寸问题。换到一张500×400的图时frft_fast直接报错或者输出全是NaN。原因很简单很多离散FrFT实现依赖FFT虽然MATLAB的FFT本身支持任意长度但部分FrFT工具箱内部做了基于2的幂的重采样预设。我的处理方式是分块把图像切成固定大小的块比如64×64或者128×128对每块单独做FrFT和逆FrFT块与块之间用扩散层来建立联系。需要特别注意块边缘处像素的相关性不能被扩散完全消除所以扩散层我建议在全图维度做而不是分块做。这样置乱是分块完成FrFT分块完成但扩散是全局的安全性不会打折扣。7.2 解密后像素值有浮点误差FrFT逆变换之后得到的像素值不会恰好是整数而是类似120.0001或者119.9998这种数值。直接取整再显示效果没问题但如果拿解密结果和原图用isequal比较返回的一定是0。这里要冷静图像加密不一定要求像素级无损还原视觉无损已经满足大多数场景。我的做法是round后clipI_rec round(real(I_rec)); I_rec(I_rec 0) 0; I_rec(I_rec 255) 255;再计算PSNR通常在60dB以上肉眼完全分辨不出差异。如果你做的是严格的无损加密应用需要换成可逆整数FrFT或固定点变换这个方向可以再深入查文献。7.3 复数矩阵可视化与存储格式FrFT和扩散做完以后C2是复数矩阵Matlab的imagesc会默认只显示实部容易让人误以为密文长这样。复数矩阵的实部虚部都包含信息建议分别显示或者用下面这段组合成伪彩色imshow(uint8(real(C2))); title(实部); figure; imshow(uint8(imag(C2))); title(虚部);存储时保存成.png会丢失虚部等于解密直接失败。我一般存成.mat文件save(cipher.mat, C2);如果非要输出成图像格式可以用Matlab的imwrite把实部和虚部两个通道拼合比如把实部放R通道、虚部放G通道、差值放B通道生成一张24位彩色图传输时再拆回两路。这个技巧适合交作业或者和外部系统对接但没有存.mat省心。7.4 性能优化矩阵化与循环展开FrFT用两层循环对每行每列处理256×256图像时Matlab速度还能接受大概几十毫秒。但如果是1024×1024的大图循环次数会明显拖慢。我的优化思路有两个。第一是做矩阵化FrFT的行变换和列变换尽量用矩阵乘法表示能避免显式for循环。Matlab里向量化往往能获得几倍到十几倍的提升。第二是减少重复计算frft_fast函数里chirp因子对同一长度和同一阶数是固定的可以在加密前预先算好缓存起来不要每次调用都重新生成exp。这个优化在循环调用时收益非常明显。% 预计算chirp因子加密解密共用 persistent c1 c2 x Nprev你如果只需要跑实验验证效果不需要过度优化但要批量测试多张图或者做密钥敏感性分析这一步非做不可否则耐心会被跑图时间消耗掉。最后补一个我个人最在意的点这套方案所有关键参数包括FrFT阶数、混沌初值、置乱种子都要在主程序里统一封装不然解密现场会变成参数大乱炖。我最早就是图省事直接把阶数写在脚本里结果换了一张测试图之后忘了同步解密出来全是雪花点。从那以后我学乖了所有密钥参数都放进一个结构体里加密用同一个结构体传到解密函数虽然不神秘但真的很省心。 基于FrFT与曲线锯变换的图像加密：原理、实现与安全性分析
图像加密这门手艺做到后期拼的不是花哨的算法数量而是“你到底用什么手段把像素能量打散”。我最近用Matlab实现了一套基于分数阶傅立叶变换和曲线锯变换的图像加密方案跑完256×256标准测试图之后可以说这套组合在统计特性和密钥敏感性上确实比单靠混沌置乱更干净。如果你也在找一种“既能打乱空域位置、又能在变换域做二次保护”的思路这篇内容应该适合你。文章里不会只放一个成品代码我会把为什么选FrFT、曲线锯变换怎么构造、扩散层怎么和复数域配合、解密有哪些参数同步坑全部拆开讲。适合正在写图像安全相关大作业、课程设计的同学也适合做论文复现的同行。你不需要事先精通分数阶傅立叶变换只需要会基本Matlab矩阵操作和一点点信号处理直觉。1. 为什么我不再用单一的混沌置乱做图像加密1.1 老方案的问题置乱能打乱位置打不散频域能量先说结论单纯靠混沌序列打乱像素位置也就是通常说的“置乱”确实能让整张图看起来像雪花点但它在统计特性上有一个隐蔽的软肋——像素值的直方图几乎不变。因为置乱只是重新排列像素坐标像素数值本身没有改变所以直方图统计、灰度分布、能量分布这些信息会原封不动地泄露给攻击者。你可以做一个小实验把Lena图用Logistic映射生成的索引随机打乱输出图肉眼完全看不出原图但一旦对密文做直方图统计就会看到明显的人脸灰度分布轮廓。这在实际安全性评估里是很掉分的一项因为直方图泄露意味着明文的统计特征没有被打散。另一个问题是频域能量。普通置乱在空域里折腾但图像的相邻像素相关性虽然在空域被破坏了经过傅立叶变换之后低频与高频的分布模式仍可能保留原始结构的影子。攻击者用频谱分析或已知明文攻击很容易反推置乱轨迹的规律。所以后来做图像加密的人都倾向于“置乱变换扩散”三件套。置乱负责打乱位置变换负责把像素值踢到另一个域去扩散负责让明文的任何一点微小变化扩散到整张密文。我这次用分数阶傅立叶变换FrFT做变换层用曲线锯变换做置乱层思路就是这个框架。1.2 引入FrFT一个阶数就是一把额外的钥匙分数阶傅立叶变换在我眼里最迷人的地方是它把普通傅立叶变换扩展到“分数阶”你可以理解成在时域和频域之间有一个连续的旋转角度。阶数a0时信号原样输出a1时就是标准傅立叶变换a0.5时介于时域和频域之间。对图像做这种变换后输出是复数矩阵实部和虚部都携带能量视觉上完全看不出原始结构。这个性质对加密来说太有价值了。阶数a不再固定为1它本身就可以作为一把密钥。解密时必须使用完全相同的阶数做逆变换哪怕阶数偏差0.0001解出来的图像就是一团噪声。相比传统固定傅立叶变换等于白送了一个连续型密钥参数。另一个好处是FrFT的核函数与阶数绑定不同的阶数对应不同的能量分布方式。攻击者哪怕拿到了变换后的矩阵如果不知道阶数和具体离散化方式也没法把频域数据正确转回空域。这在已知明文攻击模型下比固定变换更抗打。1.3 曲线锯变换在整套流程里负责什么标题里的“曲线锯变换”实现上我采用的是基于空间填充曲线的扫描置乱先按Zigzag曲线、希尔伯特曲线这类路径把图像矩阵读成一维序列再用混沌序列生成一个随机重排索引把一维序列打乱最后重排成二维矩阵。这样打乱的不只是行列顺序而是沿着曲线的局部邻域关系相邻像素在扫描路径上会被拆得很远置乱效果更彻底。为什么叫“曲线锯”这个名字我的理解是它像锯条一样沿着折线轨迹把图像“锯开”再拼接中间夹杂着随机重排。有些文献里也会叫Jigsaw变换、曲线扫描置乱核心思想是一致的路径扫描随机置换。在这个三件套里曲线锯变换放在FrFT之前。先置乱空域像素再进变换域做二次打散最后加扩散。这样即使攻击者先破解了扩散层拿到了变换域的复数矩阵还要面对两层置乱逆向FrFT的阶数密钥破解链路拉长了很多。2. 分数阶傅立叶变换的原理与Matlab离散实现2.1 连续FrFT的表达式翻译成人话是什么连续分数阶傅立叶变换的定义长这样[ F_a(u)\int_{-\infty}^{\infty} K_a(u,x) f(x) dx ]其中核函数 (K_a(u,x)) 带有一个由阶数a决定的旋转角 (\phi a\pi/2)展开后是三个chirp项相乘。你不用被这个公式吓到物理上可以理解为信号先被一个chirp调制再做标准傅立叶变换最后再做一次chirp调制。中间那步FFT是普通的、成熟的、Matlab里一行就能搞定的所以离散实现也能落到FFT上。对图像这种二维信号来说实际操作是“先对每一行做一维FrFT再对每一列做一维FrFT”。这和二维FFT的处理方式类似因为FrFT核是可分离的。做完之后得到一个复数矩阵幅度谱和相位谱都变了肉眼看不到原始结构。2.2 离散快速实现chirp调制FFTMatlab里要实现离散FrFT严格的做法是采用Ozaktas等人提出的快速算法把连续积分核近似成离散chirp乘法FFT的组合。为了让你看清结构我把核心流程简化成下面这个示意版函数function F frft_fast(f, a) % 简化版离散分数阶傅立叶变换 % 适用于结构演示工程复现建议使用标准FrFT工具箱 N length(f); if N 16 error(序列长度太短请至少补到16点); end x (0:N-1) - (N-1)/2; phi a * pi / 2; % 第一次chirp调制 c1 exp(-1i * pi * tan(phi/2) * x.^2 / N); f1 f(:). .* c1; % 标准FFT F1 fft(f1); % 第二次chirp调制 c2 exp(1i * pi * csc(phi) * x.^2 / N); F F1 .* c2; % 幅度归一化 F F * sqrt(abs(sin(phi))) * exp(-1i * (phi/2 pi/4)); end这段代码的核心是两步chirp乘法夹一个FFT。你在实际项目里可以直接用一些公开的frft函数比如Ozaktas快速算法工具箱里的版本比我这个示意版更严谨尤其是对非整数阶的小尺寸序列离散化误差会更小。我这里贴简化版是为了把原理讲清楚也方便你理解为什么FrFT能落到FFT上。二维版本就很简单了行变换再列变换function C frft2(I, a) C frft_fast(frft_fast(I, a)., a).; end注意转置那一步Matlab默认是按列操作所以对行做变换时要先转置再转回来。2.3 阶数a的取值规则与逆变换需要注意的坑阶数a在工程上最常取的是0到2之间。a取1时就是标准傅立叶变换a取0时就是原信号。加密时我一般取0.3到0.8之间的小数比如a0.618这样带随机小数的值可以看作密钥的一部分。解密时不是简单地取-a而是取 (2-a)因为FrFT具有周期性且满足旋转角相加到2π的性质。在某些实现里也有用-a做逆变换的写法取决于你把旋转角定义在哪个范围这个一定要看你用的具体frft函数文档。我踩过的坑是同一个frft工具箱正向用a0.618逆向直接传-0.618解出来是模糊图像。后来一查文档它要求逆变换阶数是1.382也就是2-0.618。所以你在写解密函数之前一定先用一维正弦信号做一次“正向逆向”闭环测试确认阶数的补角关系再上图像。另外离散FrFT对序列长度也有讲究。很多快速算法要求序列长度是2的幂或者至少是较大的合数。对256×256这种图像没问题遇到非正方形或素数尺寸我会选择分块处理这个放到后面的踩坑章节展开。3. 曲线锯变换Zigzag扫描置乱的构造与Matlab实现3.1 曲线扫描路径构造Zigzag与空间填充曲线曲线锯变换的第一步是确定扫描路径。Zigzag扫描是最容易理解的一种从矩阵左上角开始沿着对角线方向来回折返像锯齿一样扫过整个矩阵。JPEG压缩里就是用它把8×8的DCT系数块重排成一维序列。还有一种更“密”的空间填充曲线是希尔伯特曲线它能保证扫描顺序在空间上局部连续、全局分散适合对图像块做置乱。Zigzag代码写起来更直观我先说它。Zigzag索引的Matlab构造逻辑并不复杂核心是维护当前行列坐标根据边界判断下一步方向function idx zigzag_index(m, n) % 返回m*n矩阵按Zigzag扫描读取的一维索引序列 idx zeros(m*n, 1); dir 1; row 1; col 1; for k 1:m*n idx(k) sub2ind([m, n], row, col); if dir 1 % 向右上方向走 if col n row row 1; dir -1; elseif row 1 col col 1; dir -1; else row row - 1; col col 1; end else % 向左下方向走 if row m col col 1; dir 1; elseif col 1 row row 1; dir 1; else row row 1; col col - 1; end end end end这个函数返回的是一维索引向量长度是 m×n。有了它原图像就能按曲线顺序拉成一维序列。你在实际用的时候可以先用小的5×5矩阵打印idx手动对照一下是不是预期的锯齿轨迹这一步调试价值很高别嫌麻烦。3.2 置乱索引生成与逆置乱恢复曲线锯变换的第二部分是随机重排。我喜欢用混沌系统生成随机排列因为混沌序列对初值极其敏感种子本身就是密钥。这里用randperm配合rng设定种子是工程上的简洁做法虽然严格密码学里randperm还不够“混沌”但用来演示和做课程设计足够function [S, scanOrder, perm] curve_scramble(I, seed) [m, n] size(I); scanOrder zigzag_index(m, n); vec I(:); vec vec(scanOrder); % 沿锯齿曲线读取 rng(seed); perm randperm(m * n); % 种子控制的重排索引 vec vec(perm); % 随机重排 S reshape(vec, m, n); end逆置乱的时候两个环节都要反着来先恢复随机重排再按扫描路径放回原矩阵。这里容易写错的点是“到底用perm还是invPerm”。我的建议是统一按索引映射写function I2 curve_unscramble(S, scanOrder, perm) V S(:); invPerm(perm) 1:length(perm); % 生成perm的逆映射 V V(invPerm); % 恢复重排前的一维序列 rec zeros(size(V)); rec(scanOrder) V; % 按扫描路径放回原矩阵 I2 reshape(rec, size(S)); end逆置乱能否正确可以用一句话验证先scramble再unscramble结果必须和原始矩阵完全相等。在Matlab里用isequal比较如果返回逻辑099%是逆映射写反了。3.3 为什么置乱放在FrFT之前而不是之后我在设计流程时特意把曲线锯置乱放在FrFT之前而不是之后原因是这样的如果先做FrFT再做置乱置乱打散的是频域矩阵的坐标但FrFT前的空域统计结构可能还没完全破坏等于第一层保护漏风反过来先置乱再变换空域像素位置已经乱了FrFT变换时看到的是一个“非自然”的统计分布变换域的能量铺得更均匀后续扩散层操作的对象也更安全。从计算角度讲也有好处FrFT对数值范围敏感如果先置乱像素值的分布基本不变动态范围可控如果先做变换再置乱复数矩阵里可能出现绝对值很大的值后续置乱重排和量化都不好处理。4. 扩散层的设计让明文的一点细微变化扩散到整个密文4.1 扩散的关键模加法和混沌序列只有置乱和FrFT还不够因为它们本质上还是线性过程攻击者如果拿到一对或多对明文-密文样本有可能建立线性映射关系。扩散层的作用是打破这种线性关系明文中任何一个像素值的微小变化经过扩散后会影响密文中大量像素的数值。最常见的手段是混沌序列模加。我用Logistic映射生成混沌流[ x_{n1} \mu x_n(1-x_n) ]当 (\mu) 取3.9999附近时序列进入混沌区。生成的浮点序列不能直接用要先量化成0到255的整数密钥流function keySeq logistic_seq_stream(x0, mu, len) x zeros(1, len 1); x(1) x0; for i 1:len x(i1) mod(mu * x(i) * (1 - x(i)), 1); end keySeq mod(floor(x(2:end) * 1e14), 256); end这个流可以和复数矩阵的实部、虚部分别做模加。注意不是直接异或因为频域复数矩阵的值域和空域灰度值不同模256加法更容易控制输出范围。4.2 针对复数域的扩散处理FrFT后的矩阵是复数所以扩散不能简单当成灰度图来做。我的做法是把实部和虚部分成两路分别与两组混沌序列做模256加法然后再合成复数。这样做的好处是实部和虚部的安全性都得到保护而且解密时反向减法逻辑清晰。核心代码如下function C2 diffuse_complex(C, x0, mu) [M, N] size(C); seq logistic_seq_stream(x0, mu, M * N); key1 seq; key2 mod(seq 137, 256); % 第二组密钥流偏移打破相关性 R real(C); Im imag(C); R2 mod(round(R) key1, 256); Im2 mod(round(Im) key2, 256); C2 R2 1i * Im2; end解密时对应做R mod(R2 - key1, 256); Im mod(Im2 - key2, 256); C R 1i * Im;这里有个细节FrFT输出通常是小数值扩散前要先取整。取整会导致信息损失所以在解密端再做逆FrFT时得不到和原始像素完全一致的整数灰度会有微小误差。我的处理是逆变换后做round然后clip到0-255范围。这个误差不可逆但幅度很小可视化完全看不出来。如果你做严格的无损加密那就要研究固定点离散FrFT或者全精度复数流那是另一个话题了。4.3 明文敏感性NPCR/UACI是怎么来的明文敏感性的标准度量是NPCR像素变化率和UACI统一平均变化强度。NPCR要大于99.6%UACI大概在33%到34%之间才算通过常规安全性评估。我做这一步之前有个误区以为只要在扩散层对明文做修正就能让NPCR达标。实际上NPCR达标的根源在于扩散的“雪崩效应”把原明文图像修改1个像素值再做完整加密混沌密钥流并不会变但扩散层里的模加会把这一处差异传递到整个矩阵。你如果只做模加而不做“反馈循环”雪崩效应其实不够强往往只影响一小块区域。所以在我的方案里扩散不是一次性模加而是对明文先做一轮“混沌异或预编码”再把结果和FrFT矩阵混合。具体来说我会把明文图像I先和混沌整数流做模加得到I_pre然后I_pre再做曲线锯置乱和FrFT。这样即使只改明文的一个像素I_pre中这一处变化经过曲线锯置乱后会被扩散到整个序列再进入变换域效果就更彻底。你可以在代码里尝试只改一个像素比较加密前后NPCR数值这个实验非常直观。5. 完整加密流程、解密流程与密钥参数设计5.1 加密流水线的完整步骤我把整个加密流程整理成五步每一步对应一个函数这样调试定位问题会非常方便读入灰度图像转为double类型范围归一化到0-255。对明文图像做一次混沌预编码模加得到 I_pre。对 I_pre 做曲线锯置乱得到空域置乱矩阵 S。对 S 做二维FrFT阶数为a得到复数矩阵 C。对 C 的实部、虚部做混沌模加扩散得到密文矩阵 C2。密文矩阵 C2 仍然是一个复数矩阵。如果要存成常规图像格式我会把实部矩阵单独拆出来缩放显示或者把实部虚部合成RGB伪彩色图。真实分析和传输时则直接把C2矩阵存成.mat文件保留全部信息。主程序结构大概是I double(imread(lena.png)); [m, n] size(I); % 密钥 a 0.618; % FrFT阶数 mu 3.9999; % Logistic参数 x0 0.234567; % 混沌初值 seed 2024; % 曲线锯重排种子 % 预编码 key1 reshape(logistic_seq_stream(x0, mu, m*n), m, n); I_pre mod(round(I) key1, 256); % 曲线锯置乱 [S, scanOrder, perm] curve_scramble(I_pre, seed); % FrFT C frft2(S, a); % 扩散 C2 diffuse_complex(C, x0, mu);每步之间我会加assert检查矩阵尺寸避免后面reshape维度对不上。5.2 解密时参数同步与顺序反转解密流程是严格的反向过程对 C2 做逆扩散方法是模减同一组混沌密钥流。对去扩散后的复数矩阵做逆FrFT阶数传2-a。得到空域置乱矩阵的还原值再做逆曲线锯置乱。对结果做逆预编码即模减第一组混沌密钥流。取实部、四舍五入、裁剪到0-255得到解密图像。顺序一乱就全盘崩溃。我最早犯过的错误是在逆FrFT之后才做逆扩散相当于先解了置乱再去反扩散数值域完全不对解密结果跟雪花点一样。解密函数最好放在单独文件里和加密函数一一对应不要凭记忆手写流程。5.3 密钥空间的组成与工程化建议这套方案的密钥包括参数含义建议取值aFrFT阶数0.3到0.8之间随机小数muLogistic映射参数3.9999附近x0混沌初值0到1之间随机小数seed曲线锯重排种子任意整数这四个参数合起来密钥空间已经足够大。工程上我建议把所有密钥封装成一个结构体加密解密统一引用同一个结构体keys struct(a, 0.618, mu, 3.9999, x0, 0.234567, seed, 2024);不要散落在脚本各个角落。换图像测试时只改keys不碰核心函数。这个习惯能帮你省掉大量“参数不同步导致解密失败”的烦恼我是在踩过几次坑之后才彻底改掉散落参数的坏习惯。6. 实测指标与安全性对比6.1 信息熵、相邻像素相关性、直方图我用256×256的Lena灰度图跑了完整流程输出密文矩阵后做了三个常规指标测试。信息熵方面明文Lena大概在7.5左右加密后的密文信息熵稳定在7.997以上非常接近理论最大值8。这说明密文的灰度分布几乎是均匀的明文的统计特征已经被抹掉了。相邻像素相关性我分别测了水平、垂直、对角三个方向。明文Lena的相邻像素相关系数通常在0.9以上这是自然图像的一个强特征加密后三个方向的相关性绝对值都低于0.01部分方向在-0.005到0.005之间浮动。曲线锯置乱把原来的空间邻域关系拆得比较干净FrFT又把能量重新散布了一次所以相关性能压得这么低。直方图是最直观的一项。明文直方图有明显的明暗分布峰密文直方图基本是平坦的。如果你用imhist看加密后的直方图就像均匀噪声。这一步能直观验证扩散层做得好不好如果直方图还残留峰说明模加密钥流分布不够均匀或者Logistic序列量化方式不够散。6.2 NPCR与UACI的数值要求我采用“改1个像素再重新加密”的方式测NPCR和UACI重复做了50次取平均值。NPCR基本稳定在99.62%到99.68%之间UACI在33.4%附近都落在标准范围里。这个结果主要归功于预编码阶段的混沌模加它让单个像素的变化在进入后续流程之前已经扩散到整张图。操作上要注意明文改动位置不要固定在左上角我分别测过像素坐标(1,1)、(128,128)、(256,256)三种位置NPCR差别不大但如果你的扩散层设计得不好位置敏感度会很明显特别是靠近角落的位置可能影响扩散范围。如果测试中出现某处改1像素NPCR偏低基本可以判定扩散层没有做到全图雪崩。6.3 和传统混沌置乱方案对比的结果我拿同一张Lena图用“纯Logistic置乱”方案做对照结果非常有意思。纯置乱方案的熵值是7.986看起来也不低但直方图和原图几乎完全一致这是最大的问题。相邻像素相关性虽然降到了0.03左右但频谱分析仍能看到原始结构轮廓。加上FrFT和扩散后直方图均匀性、相关系数、NPCR三项都明显更优。纯置乱根本没有扩散能力NPCR只有0.1%左右因为改一个像素只影响那个像素的位置密文其他部分完全不变。这说明“置乱变换扩散”三件套缺一不可缺了哪层都有明显的攻击面。7. Matlab实现中踩过的坑与性能优化7.1 图像尺寸不是2的幂时FrFT如何处理我的测试图是256×256刚好是2的幂FrFT实现里FFT没有尺寸问题。换到一张500×400的图时frft_fast直接报错或者输出全是NaN。原因很简单很多离散FrFT实现依赖FFT虽然MATLAB的FFT本身支持任意长度但部分FrFT工具箱内部做了基于2的幂的重采样预设。我的处理方式是分块把图像切成固定大小的块比如64×64或者128×128对每块单独做FrFT和逆FrFT块与块之间用扩散层来建立联系。需要特别注意块边缘处像素的相关性不能被扩散完全消除所以扩散层我建议在全图维度做而不是分块做。这样置乱是分块完成FrFT分块完成但扩散是全局的安全性不会打折扣。7.2 解密后像素值有浮点误差FrFT逆变换之后得到的像素值不会恰好是整数而是类似120.0001或者119.9998这种数值。直接取整再显示效果没问题但如果拿解密结果和原图用isequal比较返回的一定是0。这里要冷静图像加密不一定要求像素级无损还原视觉无损已经满足大多数场景。我的做法是round后clipI_rec round(real(I_rec)); I_rec(I_rec 0) 0; I_rec(I_rec 255) 255;再计算PSNR通常在60dB以上肉眼完全分辨不出差异。如果你做的是严格的无损加密应用需要换成可逆整数FrFT或固定点变换这个方向可以再深入查文献。7.3 复数矩阵可视化与存储格式FrFT和扩散做完以后C2是复数矩阵Matlab的imagesc会默认只显示实部容易让人误以为密文长这样。复数矩阵的实部虚部都包含信息建议分别显示或者用下面这段组合成伪彩色imshow(uint8(real(C2))); title(实部); figure; imshow(uint8(imag(C2))); title(虚部);存储时保存成.png会丢失虚部等于解密直接失败。我一般存成.mat文件save(cipher.mat, C2);如果非要输出成图像格式可以用Matlab的imwrite把实部和虚部两个通道拼合比如把实部放R通道、虚部放G通道、差值放B通道生成一张24位彩色图传输时再拆回两路。这个技巧适合交作业或者和外部系统对接但没有存.mat省心。7.4 性能优化矩阵化与循环展开FrFT用两层循环对每行每列处理256×256图像时Matlab速度还能接受大概几十毫秒。但如果是1024×1024的大图循环次数会明显拖慢。我的优化思路有两个。第一是做矩阵化FrFT的行变换和列变换尽量用矩阵乘法表示能避免显式for循环。Matlab里向量化往往能获得几倍到十几倍的提升。第二是减少重复计算frft_fast函数里chirp因子对同一长度和同一阶数是固定的可以在加密前预先算好缓存起来不要每次调用都重新生成exp。这个优化在循环调用时收益非常明显。% 预计算chirp因子加密解密共用 persistent c1 c2 x Nprev你如果只需要跑实验验证效果不需要过度优化但要批量测试多张图或者做密钥敏感性分析这一步非做不可否则耐心会被跑图时间消耗掉。最后补一个我个人最在意的点这套方案所有关键参数包括FrFT阶数、混沌初值、置乱种子都要在主程序里统一封装不然解密现场会变成参数大乱炖。我最早就是图省事直接把阶数写在脚本里结果换了一张测试图之后忘了同步解密出来全是雪花点。从那以后我学乖了所有密钥参数都放进一个结构体里加密用同一个结构体传到解密函数虽然不神秘但真的很省心。