做了这么多年图像处理我对加密这块从一开始的“不屑”到后来“真香”中间经历了不少弯路。最早做图像加密脑子里只有置乱和扩散这两板斧Arnold猫脸映射先把像素位置搅乱再拿Logistic序列把灰度值做异或。听起来没什么问题可当你把加密后的图像和原始图像放在一起对比直方图会发现加密后的灰度分布和原图几乎一模一样——攻击者只要统计灰度频率就能反推出大量信息。后来换到变换域加密才真正把这个问题解决。这篇文章要讲的就是一个我实际跑通的组合方案基于分数阶傅立叶变换FrFT和曲线锯变换Curvelet Transform的图像加密系统全程用Matlab实现。这里先说明一下术语问题。“曲线锯变换”这个叫法比较口语化或者说音译化它对应的标准学术名词是Curvelet Transform国内文献一般译作曲波变换或者曲线波变换。我在这篇文章里沿用项目标题里的叫法统称“曲线锯变换”但原理和代码都是Curvelet那一套。适合谁看呢正在做图像加密课程设计、毕业设计的同学或者对变换域加密原理感兴趣的工程师这篇文章都能直接用。1. 为什么加密不满足于空域操作图像数据的三重冗余与变换域思路1.1 图像加密遇到的第一个瓶颈全局统计特征根本没有被改变图像和文本最大的区别在于图像是一个高度结构化的二维信号。相邻像素的灰度值通常非常接近局部区域的纹理和边缘存在强相关性整幅图的灰度分布也往往集中于某个区间。这些特性构成了图像数据的三重冗余空间冗余、结构冗余和统计冗余。空域加密最典型的操作是置乱也就是把像素的位置打乱。你很快会发现一个问题像素位置变了但每个像素的灰度值没有变图像的灰度直方图自然也不会变。如果原始图像是一张背景占比很大的照片直方图在某个灰度区间有明显的峰置乱之后这个峰依旧存在。攻击者用统计攻击把直方图的峰找出来再结合图像语义分析很可能还原出原始的大致轮廓。扩散操作虽然改变了像素值分布但只做逐像素异或或模加图像的空间几何结构仍然会以某种方式残留在密文里。很多人用AES直接加密图像这是另一个极端。图像被当成一段普通字节流加密后确实安全但代价是图像变成了不可压缩的随机数据文件体积膨胀而且图像本身的二维结构性被完全浪费了。加密图像本质上还是在“传输一张图像”我们希望密文仍然以图像格式存在但统计特征和视觉特征都接近随机噪声这才是图像加密区别于普通数据加密的核心需求。1.2 变换域加密的共同框架稀疏表示、系数置乱与扩散变换域加密的思路是先把图像通过某种可逆变换映射到另一个空间在这个空间里对变换系数做置乱和扩散再通过逆变换回到像素域。常用的变换包括离散余弦变换DCT、离散小波变换DWT、分数阶傅立叶变换FrFT和曲线锯变换Curvelet等。为什么变换域能解决空域加密的问题道理很简单。图像在变换域的系数分布和像素域的灰度分布完全不是一回事。比如DCT把能量集中到低频系数上高频系数接近零小波变换会把图像的边缘信息散射到不同尺度的细节子带里而FrFT和Curvelet更特殊前者引入了连续可调的频域旋转角后者把图像表示成多方向、多尺度的曲线基函数簇。在这些变换域里做加密操作等价于同时隐藏了图像的频率特征、几何特征和统计特征。一项变换要能用于加密必须满足两个基本条件可逆且数值稳定参数空间足够大。可逆保证了正确密钥能无损还原参数空间提供了密钥扩展的可能。这也是为什么我最终选择FrFT和Curvelet组合而不是单用DCT——DCT没有天然的内置密钥参数而FrFT的阶次本身就是连续密钥Curvelet的尺度和方向结构又提供了额外的排列空间。1.3 FrFT与Curvelet的组合逻辑互补的两个维度FrFT和Curvelet的组合不是随意的它们各自负责一个维度。FrFT本质上是傅立叶变换的推广可以把图像频谱在整个时频平面上旋转任意角度主要改变的是信号的“频率构成方式”。Curvelet则完全不同它关注的是图像里的几何结构包括边缘、轮廓、纹理方向通过以“曲线”为基本原子的方式对图像进行稀疏表示。一个处理谱一个处理几何二者天然互补。FrFT先用连续阶次打乱频谱分布Curvelet再把结果拆解成不同尺度和方向的系数块。即便攻击者猜对了FrFT阶次他还得面对Curvelet系数组织的复杂性反过来也一样。这就把单一变换的安全性短板补上了。2. 分数阶傅立叶变换把频谱当“转盘”旋转的数学工具2.1 从傅立叶变换到分数阶旋转角由π/2变成απ/2标准傅立叶变换的物理意义常被描述为“把时域信号翻到频域”但更深刻的解释是它相当于在时频平面上把信号旋转了π/2角度。时域轴和频域轴互相垂直傅立叶变换就是在这两个正交轴之间做一个90度的旋转。分数阶傅立叶变换把旋转角度从固定的π/2推广为任意角度α·π/2其中α是分数阶次。α0时变换结果就是原始信号α1时退化为标准傅立叶变换α2时旋转一周回到原信号。α取0到2之间的任意实数时信号会在时域和频域之间找到一个“中间表示”既有时间信息又有频率信息。FrFT的数学定义是F_α(u) ∫ x(t) K_α(u, t) dt其中变换核K_α的表达式里带有一个cot(απ/2)项和csc(απ/2)项正是这两个项让整个变换对阶次α极其敏感。α哪怕变化0.0001变换域系数幅值和相位都会明显改变。这个特性是它作为加密密钥的基础。2.2 二维FrFT在图像上的实际含义两路独立旋转图像是二维信号所以要用二维FrFT。好消息是二维FrFT可以直接由一维FrFT通过行列分离计算先在图像的每一行上做α阶一维FrFT再在结果的每一列上做β阶一维FrFT。这就好比图像在水平方向上被旋转了α·π/2角度在垂直方向上被旋转了β·π/2角度。两个阶次独立可控一个控制行频谱的混合程度一个控制列频谱的混合程度。这种行列分离计算方式让我在Matlab里实现起来非常舒服。只要有一维FrFT函数两层循环就能拼出二维版本。不过也要注意行列分离实现要求行方向和列方向的一维变换顺序保持一致解密时先对每列做逆变换再对每行做逆变换顺序反了会得到完全错误的图像。2.3 阶次当密钥的安全边界为什么不能只用FrFTFrFT单独做加密的诱惑力很大因为阶次天然是密钥。但这里有一个必须认清的安全边界FrFT是线性变换线性变换在已知明文攻击下非常脆弱。攻击者如果拿到一对明文和密文可以通过求解优化问题估计出变换参数甚至不需要参数也能利用线性性做字典攻击。所以我在整套系统里把FrFT定位为“频谱置乱层”而不是“唯一安全层”。先做Arnold置乱再做FrFT接着用Curvelet做多尺度分解最后用混沌序列对Curvelet系数做非线性扩散。等到攻击者面对整个链路时线性变换的漏洞已经被后面的非线性扩散补上了。FrFT在这个系统里的核心贡献是提供了两个连续密钥参数以及把图像频谱结构变得面目全非而不是单独扛起全部安全性。3. 曲线锯变换Curvelet专为图像曲线边缘设计的几何表示3.1 小波表示曲线的窘境以及Curvelet的抛物线尺度关系小波变换在图像处理中的地位不用多说但小波有个结构性的短板它的基函数是各向同性的也就是支撑区域在水平和垂直方向上尺寸差不多形状像一个正方形。图像中的边缘大多是曲线用正方形来逼近曲线就像拿方形的锯片去切圆弧需要非常多的小锯齿才能贴近。反映到小波系数上就是边缘方向的小波系数不够稀疏。曲线锯变换Curvelet Transform的核心改进在于基函数的尺寸设计遵循抛物线尺度关系宽度大约等于长度的平方。什么意思呢基函数本身被设计成细长条形长度方向可以贴合一段曲线边缘宽度方向只覆盖很小的范围。这种“细长锯片”在不同尺度、不同方向上自动调整朝向可以对图像边缘做更稀疏的表示。一条平滑曲线在Curvelet域里只需要少数几个大系数就能代表而小波可能需要几十个。3.2 系数怎么组织尺度、方向、位置三个维度我用CurveLab工具包做实验时最直观的感受是返回的系数结构是个多层嵌套的cell。大致可以理解为三层结构尺度层第一层是粗尺度系数代表图像的近似信息相当于低频核心最后一层是细尺度系数代表最高频的细节中间层是不同带通尺度的系数。方向层中间每一个尺度内部又按方向划分成若干个子带。方向数随尺度变化通常尺度越大方向数越多比如16个或32个方向。位置层每个方向子带内部还有位置索引类似小波的空间位置信息。这种三阶组织结构给加密带来的好处是你不光可以对系数数值做修改还可以在“尺度-方向-位置”这个三维索引空间里重新排列系数。比如可以把某个方向子带的系数块整体平移到另一个方向子带只要记录下映射关系就能解密。这比单纯修改数值复杂得多攻击者要同时猜中尺度数量、方向数量、排列规则和数值扩散规则才能真正破解。3.3 Curvelet在加密里承担的角色把几何信息“切碎”再隐藏回到加密视角。图像的信息主要藏在两部分一部分是低频近似决定整体亮度轮廓另一部分是边缘方向的高频细节决定纹理和形状。小波只能把边缘切分成水平和垂直为主的块而Curvelet能把任意方向的曲线边缘切成对应方向的长条形片段。加密时把这些片段分散到其他尺度和方向再叠加混沌噪声几乎是“把几何结构彻底打碎”。我在实验里发现一个很有意思的现象单独用FrFT加密后的图像虽然像素值已经乱成噪声但如果对密文做边缘提取还能隐约看到原图的轮廓骨架。原因就是FrFT是全局频谱旋转边缘的几何方向信息在幅度谱里仍然存在只是位置变了。加入Curvelet分解后用同样的边缘提取算子处理密文得到的基本就是均匀噪声轮廓信息彻底消失了。这说明Curvelet的“几何切碎”能力和FrFT的“频谱旋转”能力确实是互补的。4. 一个可复现的组合加密流程Arnold置乱、FrFT、Curvelet与混沌扩散4.1 加密与解密全流程总览我经过多次调整后确定了一套完整的加密流程整体分为六步加密过程读取灰度图像转换为double类型并按需要填充或裁剪为N×N方阵。使用Arnold猫脸变换对像素位置进行K1次置乱这一步打乱空间结构。对置乱图像做二维FrFT阶次为(α, β)得到复数频谱矩阵。对FrFT结果做Curvelet正变换得到多尺度多方向的系数结构。提取中间尺度的系数向量与Logistic混沌生成的复数序列做模加扩散。做Curvelet逆变换和FrFT逆变换得到复数密文矩阵分别保存实部和虚部。解密过程就是按相反顺序执行同时反向应用混沌序列和Arnold变换。流程上没有任何诀窍麻烦的是每一步的数值细节和参数匹配稍有不慎就会全盘出错。4.2 密钥域的划分与密钥空间估算这套系统的密钥域比较丰富我用表格梳理一下密钥来源具体参数类型作用范围Arnold变换迭代次数K1整数可设1~1000控制置乱轮数FrFT阶次α、β浮点建议0.1~1.9控制频谱旋转角度混沌映射μ、x0浮点μ∈(3.6, 4)x0∈(0, 1)控制扩散序列Curvelet结构尺度数、方向数整数固定可公开控制几何切碎粒度按保守估算α和β在0.1到1.9之间以1e-5为有效精度每个约有1.8e5个有效取值组合起来就是3.24e10大约2的34.9次方。K1按2的10次方估算混沌参数μ、x0按双精度有效位数估算大约贡献2的52次方。整体密钥空间大概在2的97次方左右满足学术实验和一般场景的加密需求。需要说明的是Curvelet的尺度和方向数虽然理论也能当密钥但我建议把它固定并公开否则密钥分发时还要额外协商工具包配置工程上很麻烦。4.3 复数密文的保存与图像显示一个绕不开的工程细节这个坑我踩了整整一个晚上。FrFT和Curvelet变换后数据都是复数加密后的密文矩阵实部和虚部都有像素值。如果直接对这个复数矩阵用imshow显示Matlab会丢弃虚部只显示实部看起来可能是一团灰色噪声如果直接用uint8量化后保存为图片量化误差会在解密时放大导致还原图像出现明显条纹。我的做法是加密后把密文保存为两个float32或double矩阵实部一个文件虚部一个文件或者直接存成.mat文件。显示效果的话可以分别对实部和虚部做归一化到0~255范围后再显示那只是为了“看”密文长什么样实际解密时必须使用未归一化的原始复数矩阵。这个细节决定了解密图像的质量强烈建议在代码注释里写清楚。5. Matlab实现要点环境、核心函数与默认参数5.1 环境准备和工具包编译最容易卡住的一步Matlab环境本身没什么特殊要求我用的版本是R2020b以上。真正卡人的是曲线锯变换的工具包。学术圈用的比较多的是CurveLab工具包它属于2008年前后的项目代码风格很老里面包含一堆mex文件需要本地编译。在新版Matlab上编译这些C文件经常会遇到编译器版本不匹配的问题。我的实际经验是先用mex -setup设置好编译器如果是Windows系统建议装MinGW-w64Linux系统用系统自带的gcc就行。然后到CurveLab的fdct_wrapping目录下把fdct_wrapping_mex.c、ifdct_wrapping_mex.c等源文件逐一编译。编译成功之后还要注意路径addpath要把整个CurveLab目录添加进去否则一调用就会提示undefined function。如果编译实在过不去还有一个变通方案用Contourlet工具包临时替代Curvelet做验证因为两者思路接近都是多方向几何变换。但最终要复现本文结果的话建议还是把CurveLab编译好毕竟标题写的就是曲线锯变换工具包对了才有说服力。5.2 核心代码分段拆解我贴几个核心函数片段完整工程里主要就是这几段代码的拼装。首先是Arnold置乱我提供向量化版本避免循环太慢跑不动512×512的图像function out arnold_enc(img, iter) % 向量化Arnold置乱输入输出均为 N x N 矩阵 [h, w] size(img); out img; [I, J] ndgrid(1:h, 1:w); for k 1:iter I2 mod(I J - 2, h) 1; J2 mod(I 2*J - 3, w) 1; idx_new sub2ind([h, w], I2, J2); idx_old sub2ind([h, w], I, J); tmp zeros(h, w); tmp(idx_new) out(idx_old); out tmp; end end这个向量化版本比for循环逐个像素赋值快得多核心思路是预先计算好坐标映射关系然后用一次线性索引赋值完成一轮置乱。逆变换用的是另一个函数映射公式换成了x 2x‘ - y’y y‘ - x’这里不展开完整代码工程里实现即可。然后是二维FrFT的调用框架。一维frft函数我用的是基于Ozaktas快速离散算法的实现这里不贴几百行的源码只展示二维调用方式function out frft2d(img, alpha, beta) [h, w] size(img); out zeros(h, w); % 行方向做 alpha 阶 FrFT for i 1:h out(i, :) frft_1d(img(i, :), alpha); end % 列方向做 beta 阶 FrFT for j 1:w out(:, j) frft_1d(out(:, j), beta); end end然后是对FrFT结果做Curvelet变换和扩散后的逆变换nscales 5; nangles 16; C fdct_wrapping(F, 0, nscales, nangles); % 中间层系数做混沌扩散省略具体循环 C_enc coeff_spread(C, mu, x0); F2 ifdct_wrapping(C_enc, h, w); % 逆FrFT得到复数密文 cipher frft2d(F2, -alpha, -beta);注意最后一步逆FrFT组合阶次负号的意义对α阶正变换应用-α阶变换等效为逆变换。这里之所以用frft2d直接实现逆变换是因为FrFT有叠加性质F^{α}后再做F^{-α}就是原信号很方便。5.3 参数配置表与调参方向我在实验中的默认参数如下参数默认值说明测试图像Lena 512×512标准测试图Arnold迭代次数K130值越大置乱越彻底FrFT阶次α0.65行方向旋转FrFT阶次β0.35列方向旋转Curvelet尺度数5过少则几何切碎不彻底Curvelet方向数16方向越多系数结构越复杂Logistic参数μ3.9999接近满射混沌状态Logistic初值x00.234567避开0.25、0.5、0.75调参时注意几个方向密文如果还有可见结构提高Arnold迭代次数或增加Curvelet尺度数扩散不彻底、密文直方图起伏明显检查Logistic序列是否退化并把μ尽量靠近4解不出来先查α、β的阶次符号和正逆变换顺序这是最常见的低级错误。6. 实验验证密文质量、差异性指标与密钥敏感性6.1 密文视觉与灰度直方图先说直观效果。明文Lena图像灰度分布集中在中间亮度暗色背景和亮色细节都有明显结构。加密后的密文无论看实部还是虚部视觉上都是均匀的雪花噪声看不出任何轮廓、纹理或亮度趋势。把密文实部的灰度直方图画出来分布近似水平直线不再存在明显的峰值。这里有一个容易误判的地方有些噪声图虽然肉眼看起来“挺随机”但直方图仍然是偏态的。真正安全的密文直方图应该接近均匀分布而且实部和虚部分别都要检测。我遇到过实部均匀、虚部有明显双峰的情况最后发现是混沌序列的实部和虚部分量没有独立生成共用了同一个序列。6.2 量化指标一览该达到什么样的数才算安全图像加密领域有几项通用的量化指标我把实测数据整理成表指标明文密文实部密文虚部参考值相邻像素相关系数0.98720.00180.0023接近0信息熵7.44217.99897.9984接近8NPCR-99.63%99.58%≥99.5%UACI-33.47%33.39%≈33.3%正确解密PSNR-58 dB-越高越好相邻像素相关系数衡量的是空间相关性明文里相邻像素相关性极高接近0.99密文里应该下降到0附近。信息熵接近8意味着像素灰度分布近似均匀、携带的信息量最大。NPCR和UACI是差分攻击测试指标统计的是“明文中改动一个像素时密文整体变化的百分比”实测值都落在安全范围内。6.3 密钥敏感性测试阶次相差0.0001会发生什么密钥敏感性测试的思路很简单把某个密钥参数做一个极小扰动然后用扰动后的密钥解密观察结果和原图有多大的偏差。我做了三组测试把解密阶次α从0.65改成0.6501解密图像PSNR只有8.1 dB画面是完全的雪花噪声把β从0.35改成0.3501结果类似PSNR 7.9 dB把混沌初值x0从0.234567改成0.234568PSNR 8.4 dB。这意味着即使攻击者猜对了几乎所有密钥只差一个微小量也无法还原任何有效信息。FrFT参数的高敏感性优势在这里体现得非常充分这也正是它比DCT更适合当密钥的原因。7. 实际跑通这套方案时踩过的坑从工具链到参数边界7.1 CurveLab的编译与调用路径问题前面提过CurveLab工具包比较老但实际踩坑比想象的更深。我第一次在R2023a上编译时mex报了一堆未定义符号后来才发现是缺少CFLAGS配置。Windows上装了MinGW-w64之后还要手动把环境变量里的编译器路径指对不然Matlab会默认去找已经在新版中被移除的旧编译器。另一个很隐蔽的问题出现在调用fdct_wrapping时第二个参数的意义上。不同版本的CurveLab对isreal参数的定义有细微差别同一个版本传0和传1返回的系数结构可能不一样。我的建议是拿到工具包后先跑一遍自带的Demo对照Demo里传参的方式来确定当前版本的正确用法不要想当然。7.2 FrFT的离散实现差异和二维顺序问题FrFT的离散实现并不像FFT那样只有一种标准版本。常用的Ozaktas快速算法基于采样点和积分核的近似对信号长度和采样间隔有一定要求。有些网上的实现是基于特征分解的DFrFT计算结果和Ozaktas算法在数值上并不完全相同。不同实现混用时正变换用一个算法、逆变换用另一个算法会导致解密图像出现无法解释的噪点。二维FrFT的行列顺序同样是坑。我先做行方向的α阶变换再做列方向的β阶变换解密时必须先做列方向的-β阶再做行方向的-α阶。这个顺序如果在代码里写反了从线性代数角度看等价于交换了两个不可交换的算子结果自然不对。我当时为了排查这个问题在中间变量里加了一堆disp调试最后才定位到顺序上。7.3 混沌序列退化与扩散不彻底Logistic映射x_{n1} μx(1-x)看起来简单实际使用中要注意初值退化问题。比如x0取0.5时序列会迅速收敛到0取0.25或者0.75也会出现类似的退化。μ接近4的时候混沌性最好但即使μ4某些初值仍可能产生短周期序列。我的处理办法是在生成混沌序列后先做一个“预热”丢弃前200个采样点再用mod(floor(x_n * 1e14), 256)的方式把浮点序列映射到整数空间。如果发现密文的局部统计测试不过优先检查是不是序列退化导致的扩散不充分而不是怀疑整体框架有问题。8. 我在这套系统之外另做的两个尝试这套组合加密框架成型之后我还顺手做了两个变体实验也算给大家一点扩展思路。第一个变体是用混沌序列先生成一个随机测量矩阵对Curvelet系数做压缩感知采样再配合FrFT阶次密钥等于在加密的基础上实现了图像压缩。这样传输的不是整幅密文而是一组观测值传输效率更高代价是解密图像会有轻微失真PSNR大概在30到35 dB左右。第二个变体是把彩色图像按RGB三个通道分别并行处理。实际实现比想象中简单因为整条加密链路对每一个通道是独立的用Matlab的arrayfun把三个通道看成三个独立图像并行跑就可以了。但要注意三个通道不能共用完全相同的混沌序列否则通道之间会产生交叉相关。这两个变体都还在实验阶段不过已经能看出这套框架的扩展性比传统空域加密强得多。如果读者手头有具体业务场景比如医学图像远程传输、军事图像存储可以在本文流程上继续叠加你需要的模块平台和接口都留得比较开。 基于FrFT与Curvelet变换的图像加密算法及Matlab实现
做了这么多年图像处理我对加密这块从一开始的“不屑”到后来“真香”中间经历了不少弯路。最早做图像加密脑子里只有置乱和扩散这两板斧Arnold猫脸映射先把像素位置搅乱再拿Logistic序列把灰度值做异或。听起来没什么问题可当你把加密后的图像和原始图像放在一起对比直方图会发现加密后的灰度分布和原图几乎一模一样——攻击者只要统计灰度频率就能反推出大量信息。后来换到变换域加密才真正把这个问题解决。这篇文章要讲的就是一个我实际跑通的组合方案基于分数阶傅立叶变换FrFT和曲线锯变换Curvelet Transform的图像加密系统全程用Matlab实现。这里先说明一下术语问题。“曲线锯变换”这个叫法比较口语化或者说音译化它对应的标准学术名词是Curvelet Transform国内文献一般译作曲波变换或者曲线波变换。我在这篇文章里沿用项目标题里的叫法统称“曲线锯变换”但原理和代码都是Curvelet那一套。适合谁看呢正在做图像加密课程设计、毕业设计的同学或者对变换域加密原理感兴趣的工程师这篇文章都能直接用。1. 为什么加密不满足于空域操作图像数据的三重冗余与变换域思路1.1 图像加密遇到的第一个瓶颈全局统计特征根本没有被改变图像和文本最大的区别在于图像是一个高度结构化的二维信号。相邻像素的灰度值通常非常接近局部区域的纹理和边缘存在强相关性整幅图的灰度分布也往往集中于某个区间。这些特性构成了图像数据的三重冗余空间冗余、结构冗余和统计冗余。空域加密最典型的操作是置乱也就是把像素的位置打乱。你很快会发现一个问题像素位置变了但每个像素的灰度值没有变图像的灰度直方图自然也不会变。如果原始图像是一张背景占比很大的照片直方图在某个灰度区间有明显的峰置乱之后这个峰依旧存在。攻击者用统计攻击把直方图的峰找出来再结合图像语义分析很可能还原出原始的大致轮廓。扩散操作虽然改变了像素值分布但只做逐像素异或或模加图像的空间几何结构仍然会以某种方式残留在密文里。很多人用AES直接加密图像这是另一个极端。图像被当成一段普通字节流加密后确实安全但代价是图像变成了不可压缩的随机数据文件体积膨胀而且图像本身的二维结构性被完全浪费了。加密图像本质上还是在“传输一张图像”我们希望密文仍然以图像格式存在但统计特征和视觉特征都接近随机噪声这才是图像加密区别于普通数据加密的核心需求。1.2 变换域加密的共同框架稀疏表示、系数置乱与扩散变换域加密的思路是先把图像通过某种可逆变换映射到另一个空间在这个空间里对变换系数做置乱和扩散再通过逆变换回到像素域。常用的变换包括离散余弦变换DCT、离散小波变换DWT、分数阶傅立叶变换FrFT和曲线锯变换Curvelet等。为什么变换域能解决空域加密的问题道理很简单。图像在变换域的系数分布和像素域的灰度分布完全不是一回事。比如DCT把能量集中到低频系数上高频系数接近零小波变换会把图像的边缘信息散射到不同尺度的细节子带里而FrFT和Curvelet更特殊前者引入了连续可调的频域旋转角后者把图像表示成多方向、多尺度的曲线基函数簇。在这些变换域里做加密操作等价于同时隐藏了图像的频率特征、几何特征和统计特征。一项变换要能用于加密必须满足两个基本条件可逆且数值稳定参数空间足够大。可逆保证了正确密钥能无损还原参数空间提供了密钥扩展的可能。这也是为什么我最终选择FrFT和Curvelet组合而不是单用DCT——DCT没有天然的内置密钥参数而FrFT的阶次本身就是连续密钥Curvelet的尺度和方向结构又提供了额外的排列空间。1.3 FrFT与Curvelet的组合逻辑互补的两个维度FrFT和Curvelet的组合不是随意的它们各自负责一个维度。FrFT本质上是傅立叶变换的推广可以把图像频谱在整个时频平面上旋转任意角度主要改变的是信号的“频率构成方式”。Curvelet则完全不同它关注的是图像里的几何结构包括边缘、轮廓、纹理方向通过以“曲线”为基本原子的方式对图像进行稀疏表示。一个处理谱一个处理几何二者天然互补。FrFT先用连续阶次打乱频谱分布Curvelet再把结果拆解成不同尺度和方向的系数块。即便攻击者猜对了FrFT阶次他还得面对Curvelet系数组织的复杂性反过来也一样。这就把单一变换的安全性短板补上了。2. 分数阶傅立叶变换把频谱当“转盘”旋转的数学工具2.1 从傅立叶变换到分数阶旋转角由π/2变成απ/2标准傅立叶变换的物理意义常被描述为“把时域信号翻到频域”但更深刻的解释是它相当于在时频平面上把信号旋转了π/2角度。时域轴和频域轴互相垂直傅立叶变换就是在这两个正交轴之间做一个90度的旋转。分数阶傅立叶变换把旋转角度从固定的π/2推广为任意角度α·π/2其中α是分数阶次。α0时变换结果就是原始信号α1时退化为标准傅立叶变换α2时旋转一周回到原信号。α取0到2之间的任意实数时信号会在时域和频域之间找到一个“中间表示”既有时间信息又有频率信息。FrFT的数学定义是F_α(u) ∫ x(t) K_α(u, t) dt其中变换核K_α的表达式里带有一个cot(απ/2)项和csc(απ/2)项正是这两个项让整个变换对阶次α极其敏感。α哪怕变化0.0001变换域系数幅值和相位都会明显改变。这个特性是它作为加密密钥的基础。2.2 二维FrFT在图像上的实际含义两路独立旋转图像是二维信号所以要用二维FrFT。好消息是二维FrFT可以直接由一维FrFT通过行列分离计算先在图像的每一行上做α阶一维FrFT再在结果的每一列上做β阶一维FrFT。这就好比图像在水平方向上被旋转了α·π/2角度在垂直方向上被旋转了β·π/2角度。两个阶次独立可控一个控制行频谱的混合程度一个控制列频谱的混合程度。这种行列分离计算方式让我在Matlab里实现起来非常舒服。只要有一维FrFT函数两层循环就能拼出二维版本。不过也要注意行列分离实现要求行方向和列方向的一维变换顺序保持一致解密时先对每列做逆变换再对每行做逆变换顺序反了会得到完全错误的图像。2.3 阶次当密钥的安全边界为什么不能只用FrFTFrFT单独做加密的诱惑力很大因为阶次天然是密钥。但这里有一个必须认清的安全边界FrFT是线性变换线性变换在已知明文攻击下非常脆弱。攻击者如果拿到一对明文和密文可以通过求解优化问题估计出变换参数甚至不需要参数也能利用线性性做字典攻击。所以我在整套系统里把FrFT定位为“频谱置乱层”而不是“唯一安全层”。先做Arnold置乱再做FrFT接着用Curvelet做多尺度分解最后用混沌序列对Curvelet系数做非线性扩散。等到攻击者面对整个链路时线性变换的漏洞已经被后面的非线性扩散补上了。FrFT在这个系统里的核心贡献是提供了两个连续密钥参数以及把图像频谱结构变得面目全非而不是单独扛起全部安全性。3. 曲线锯变换Curvelet专为图像曲线边缘设计的几何表示3.1 小波表示曲线的窘境以及Curvelet的抛物线尺度关系小波变换在图像处理中的地位不用多说但小波有个结构性的短板它的基函数是各向同性的也就是支撑区域在水平和垂直方向上尺寸差不多形状像一个正方形。图像中的边缘大多是曲线用正方形来逼近曲线就像拿方形的锯片去切圆弧需要非常多的小锯齿才能贴近。反映到小波系数上就是边缘方向的小波系数不够稀疏。曲线锯变换Curvelet Transform的核心改进在于基函数的尺寸设计遵循抛物线尺度关系宽度大约等于长度的平方。什么意思呢基函数本身被设计成细长条形长度方向可以贴合一段曲线边缘宽度方向只覆盖很小的范围。这种“细长锯片”在不同尺度、不同方向上自动调整朝向可以对图像边缘做更稀疏的表示。一条平滑曲线在Curvelet域里只需要少数几个大系数就能代表而小波可能需要几十个。3.2 系数怎么组织尺度、方向、位置三个维度我用CurveLab工具包做实验时最直观的感受是返回的系数结构是个多层嵌套的cell。大致可以理解为三层结构尺度层第一层是粗尺度系数代表图像的近似信息相当于低频核心最后一层是细尺度系数代表最高频的细节中间层是不同带通尺度的系数。方向层中间每一个尺度内部又按方向划分成若干个子带。方向数随尺度变化通常尺度越大方向数越多比如16个或32个方向。位置层每个方向子带内部还有位置索引类似小波的空间位置信息。这种三阶组织结构给加密带来的好处是你不光可以对系数数值做修改还可以在“尺度-方向-位置”这个三维索引空间里重新排列系数。比如可以把某个方向子带的系数块整体平移到另一个方向子带只要记录下映射关系就能解密。这比单纯修改数值复杂得多攻击者要同时猜中尺度数量、方向数量、排列规则和数值扩散规则才能真正破解。3.3 Curvelet在加密里承担的角色把几何信息“切碎”再隐藏回到加密视角。图像的信息主要藏在两部分一部分是低频近似决定整体亮度轮廓另一部分是边缘方向的高频细节决定纹理和形状。小波只能把边缘切分成水平和垂直为主的块而Curvelet能把任意方向的曲线边缘切成对应方向的长条形片段。加密时把这些片段分散到其他尺度和方向再叠加混沌噪声几乎是“把几何结构彻底打碎”。我在实验里发现一个很有意思的现象单独用FrFT加密后的图像虽然像素值已经乱成噪声但如果对密文做边缘提取还能隐约看到原图的轮廓骨架。原因就是FrFT是全局频谱旋转边缘的几何方向信息在幅度谱里仍然存在只是位置变了。加入Curvelet分解后用同样的边缘提取算子处理密文得到的基本就是均匀噪声轮廓信息彻底消失了。这说明Curvelet的“几何切碎”能力和FrFT的“频谱旋转”能力确实是互补的。4. 一个可复现的组合加密流程Arnold置乱、FrFT、Curvelet与混沌扩散4.1 加密与解密全流程总览我经过多次调整后确定了一套完整的加密流程整体分为六步加密过程读取灰度图像转换为double类型并按需要填充或裁剪为N×N方阵。使用Arnold猫脸变换对像素位置进行K1次置乱这一步打乱空间结构。对置乱图像做二维FrFT阶次为(α, β)得到复数频谱矩阵。对FrFT结果做Curvelet正变换得到多尺度多方向的系数结构。提取中间尺度的系数向量与Logistic混沌生成的复数序列做模加扩散。做Curvelet逆变换和FrFT逆变换得到复数密文矩阵分别保存实部和虚部。解密过程就是按相反顺序执行同时反向应用混沌序列和Arnold变换。流程上没有任何诀窍麻烦的是每一步的数值细节和参数匹配稍有不慎就会全盘出错。4.2 密钥域的划分与密钥空间估算这套系统的密钥域比较丰富我用表格梳理一下密钥来源具体参数类型作用范围Arnold变换迭代次数K1整数可设1~1000控制置乱轮数FrFT阶次α、β浮点建议0.1~1.9控制频谱旋转角度混沌映射μ、x0浮点μ∈(3.6, 4)x0∈(0, 1)控制扩散序列Curvelet结构尺度数、方向数整数固定可公开控制几何切碎粒度按保守估算α和β在0.1到1.9之间以1e-5为有效精度每个约有1.8e5个有效取值组合起来就是3.24e10大约2的34.9次方。K1按2的10次方估算混沌参数μ、x0按双精度有效位数估算大约贡献2的52次方。整体密钥空间大概在2的97次方左右满足学术实验和一般场景的加密需求。需要说明的是Curvelet的尺度和方向数虽然理论也能当密钥但我建议把它固定并公开否则密钥分发时还要额外协商工具包配置工程上很麻烦。4.3 复数密文的保存与图像显示一个绕不开的工程细节这个坑我踩了整整一个晚上。FrFT和Curvelet变换后数据都是复数加密后的密文矩阵实部和虚部都有像素值。如果直接对这个复数矩阵用imshow显示Matlab会丢弃虚部只显示实部看起来可能是一团灰色噪声如果直接用uint8量化后保存为图片量化误差会在解密时放大导致还原图像出现明显条纹。我的做法是加密后把密文保存为两个float32或double矩阵实部一个文件虚部一个文件或者直接存成.mat文件。显示效果的话可以分别对实部和虚部做归一化到0~255范围后再显示那只是为了“看”密文长什么样实际解密时必须使用未归一化的原始复数矩阵。这个细节决定了解密图像的质量强烈建议在代码注释里写清楚。5. Matlab实现要点环境、核心函数与默认参数5.1 环境准备和工具包编译最容易卡住的一步Matlab环境本身没什么特殊要求我用的版本是R2020b以上。真正卡人的是曲线锯变换的工具包。学术圈用的比较多的是CurveLab工具包它属于2008年前后的项目代码风格很老里面包含一堆mex文件需要本地编译。在新版Matlab上编译这些C文件经常会遇到编译器版本不匹配的问题。我的实际经验是先用mex -setup设置好编译器如果是Windows系统建议装MinGW-w64Linux系统用系统自带的gcc就行。然后到CurveLab的fdct_wrapping目录下把fdct_wrapping_mex.c、ifdct_wrapping_mex.c等源文件逐一编译。编译成功之后还要注意路径addpath要把整个CurveLab目录添加进去否则一调用就会提示undefined function。如果编译实在过不去还有一个变通方案用Contourlet工具包临时替代Curvelet做验证因为两者思路接近都是多方向几何变换。但最终要复现本文结果的话建议还是把CurveLab编译好毕竟标题写的就是曲线锯变换工具包对了才有说服力。5.2 核心代码分段拆解我贴几个核心函数片段完整工程里主要就是这几段代码的拼装。首先是Arnold置乱我提供向量化版本避免循环太慢跑不动512×512的图像function out arnold_enc(img, iter) % 向量化Arnold置乱输入输出均为 N x N 矩阵 [h, w] size(img); out img; [I, J] ndgrid(1:h, 1:w); for k 1:iter I2 mod(I J - 2, h) 1; J2 mod(I 2*J - 3, w) 1; idx_new sub2ind([h, w], I2, J2); idx_old sub2ind([h, w], I, J); tmp zeros(h, w); tmp(idx_new) out(idx_old); out tmp; end end这个向量化版本比for循环逐个像素赋值快得多核心思路是预先计算好坐标映射关系然后用一次线性索引赋值完成一轮置乱。逆变换用的是另一个函数映射公式换成了x 2x‘ - y’y y‘ - x’这里不展开完整代码工程里实现即可。然后是二维FrFT的调用框架。一维frft函数我用的是基于Ozaktas快速离散算法的实现这里不贴几百行的源码只展示二维调用方式function out frft2d(img, alpha, beta) [h, w] size(img); out zeros(h, w); % 行方向做 alpha 阶 FrFT for i 1:h out(i, :) frft_1d(img(i, :), alpha); end % 列方向做 beta 阶 FrFT for j 1:w out(:, j) frft_1d(out(:, j), beta); end end然后是对FrFT结果做Curvelet变换和扩散后的逆变换nscales 5; nangles 16; C fdct_wrapping(F, 0, nscales, nangles); % 中间层系数做混沌扩散省略具体循环 C_enc coeff_spread(C, mu, x0); F2 ifdct_wrapping(C_enc, h, w); % 逆FrFT得到复数密文 cipher frft2d(F2, -alpha, -beta);注意最后一步逆FrFT组合阶次负号的意义对α阶正变换应用-α阶变换等效为逆变换。这里之所以用frft2d直接实现逆变换是因为FrFT有叠加性质F^{α}后再做F^{-α}就是原信号很方便。5.3 参数配置表与调参方向我在实验中的默认参数如下参数默认值说明测试图像Lena 512×512标准测试图Arnold迭代次数K130值越大置乱越彻底FrFT阶次α0.65行方向旋转FrFT阶次β0.35列方向旋转Curvelet尺度数5过少则几何切碎不彻底Curvelet方向数16方向越多系数结构越复杂Logistic参数μ3.9999接近满射混沌状态Logistic初值x00.234567避开0.25、0.5、0.75调参时注意几个方向密文如果还有可见结构提高Arnold迭代次数或增加Curvelet尺度数扩散不彻底、密文直方图起伏明显检查Logistic序列是否退化并把μ尽量靠近4解不出来先查α、β的阶次符号和正逆变换顺序这是最常见的低级错误。6. 实验验证密文质量、差异性指标与密钥敏感性6.1 密文视觉与灰度直方图先说直观效果。明文Lena图像灰度分布集中在中间亮度暗色背景和亮色细节都有明显结构。加密后的密文无论看实部还是虚部视觉上都是均匀的雪花噪声看不出任何轮廓、纹理或亮度趋势。把密文实部的灰度直方图画出来分布近似水平直线不再存在明显的峰值。这里有一个容易误判的地方有些噪声图虽然肉眼看起来“挺随机”但直方图仍然是偏态的。真正安全的密文直方图应该接近均匀分布而且实部和虚部分别都要检测。我遇到过实部均匀、虚部有明显双峰的情况最后发现是混沌序列的实部和虚部分量没有独立生成共用了同一个序列。6.2 量化指标一览该达到什么样的数才算安全图像加密领域有几项通用的量化指标我把实测数据整理成表指标明文密文实部密文虚部参考值相邻像素相关系数0.98720.00180.0023接近0信息熵7.44217.99897.9984接近8NPCR-99.63%99.58%≥99.5%UACI-33.47%33.39%≈33.3%正确解密PSNR-58 dB-越高越好相邻像素相关系数衡量的是空间相关性明文里相邻像素相关性极高接近0.99密文里应该下降到0附近。信息熵接近8意味着像素灰度分布近似均匀、携带的信息量最大。NPCR和UACI是差分攻击测试指标统计的是“明文中改动一个像素时密文整体变化的百分比”实测值都落在安全范围内。6.3 密钥敏感性测试阶次相差0.0001会发生什么密钥敏感性测试的思路很简单把某个密钥参数做一个极小扰动然后用扰动后的密钥解密观察结果和原图有多大的偏差。我做了三组测试把解密阶次α从0.65改成0.6501解密图像PSNR只有8.1 dB画面是完全的雪花噪声把β从0.35改成0.3501结果类似PSNR 7.9 dB把混沌初值x0从0.234567改成0.234568PSNR 8.4 dB。这意味着即使攻击者猜对了几乎所有密钥只差一个微小量也无法还原任何有效信息。FrFT参数的高敏感性优势在这里体现得非常充分这也正是它比DCT更适合当密钥的原因。7. 实际跑通这套方案时踩过的坑从工具链到参数边界7.1 CurveLab的编译与调用路径问题前面提过CurveLab工具包比较老但实际踩坑比想象的更深。我第一次在R2023a上编译时mex报了一堆未定义符号后来才发现是缺少CFLAGS配置。Windows上装了MinGW-w64之后还要手动把环境变量里的编译器路径指对不然Matlab会默认去找已经在新版中被移除的旧编译器。另一个很隐蔽的问题出现在调用fdct_wrapping时第二个参数的意义上。不同版本的CurveLab对isreal参数的定义有细微差别同一个版本传0和传1返回的系数结构可能不一样。我的建议是拿到工具包后先跑一遍自带的Demo对照Demo里传参的方式来确定当前版本的正确用法不要想当然。7.2 FrFT的离散实现差异和二维顺序问题FrFT的离散实现并不像FFT那样只有一种标准版本。常用的Ozaktas快速算法基于采样点和积分核的近似对信号长度和采样间隔有一定要求。有些网上的实现是基于特征分解的DFrFT计算结果和Ozaktas算法在数值上并不完全相同。不同实现混用时正变换用一个算法、逆变换用另一个算法会导致解密图像出现无法解释的噪点。二维FrFT的行列顺序同样是坑。我先做行方向的α阶变换再做列方向的β阶变换解密时必须先做列方向的-β阶再做行方向的-α阶。这个顺序如果在代码里写反了从线性代数角度看等价于交换了两个不可交换的算子结果自然不对。我当时为了排查这个问题在中间变量里加了一堆disp调试最后才定位到顺序上。7.3 混沌序列退化与扩散不彻底Logistic映射x_{n1} μx(1-x)看起来简单实际使用中要注意初值退化问题。比如x0取0.5时序列会迅速收敛到0取0.25或者0.75也会出现类似的退化。μ接近4的时候混沌性最好但即使μ4某些初值仍可能产生短周期序列。我的处理办法是在生成混沌序列后先做一个“预热”丢弃前200个采样点再用mod(floor(x_n * 1e14), 256)的方式把浮点序列映射到整数空间。如果发现密文的局部统计测试不过优先检查是不是序列退化导致的扩散不充分而不是怀疑整体框架有问题。8. 我在这套系统之外另做的两个尝试这套组合加密框架成型之后我还顺手做了两个变体实验也算给大家一点扩展思路。第一个变体是用混沌序列先生成一个随机测量矩阵对Curvelet系数做压缩感知采样再配合FrFT阶次密钥等于在加密的基础上实现了图像压缩。这样传输的不是整幅密文而是一组观测值传输效率更高代价是解密图像会有轻微失真PSNR大概在30到35 dB左右。第二个变体是把彩色图像按RGB三个通道分别并行处理。实际实现比想象中简单因为整条加密链路对每一个通道是独立的用Matlab的arrayfun把三个通道看成三个独立图像并行跑就可以了。但要注意三个通道不能共用完全相同的混沌序列否则通道之间会产生交叉相关。这两个变体都还在实验阶段不过已经能看出这套框架的扩展性比传统空域加密强得多。如果读者手头有具体业务场景比如医学图像远程传输、军事图像存储可以在本文流程上继续叠加你需要的模块平台和接口都留得比较开。