引　言

光学信息安全技术有着高维度、高并行处理速度以及能快速实现卷积和相关运算的特点^[1-3]，在信息安全领域，光学信息安全技术的优势越来越明显。其中，相位检索算法在光学安全领域引起了广泛的关注。但基于单图像的研究不能满足当下人们对信息处理效率的要求。与单幅图像的光学加密技术相比，多图像加密技术可以在一幅图像中包含多重图像的信息，提高加密图像的容量和效率，这是多用户认证和内容分散的基础。

多图像光学加密技术作为光学加密的一个重要分支，不仅提高了加密能力，而且减少了密文传输数据。以关联成像（ghost imaging, GI）为基础的多图像加密方法相继被提出。Li等^[4]提出了一种基于关联成像和坐标采样的多图像加密方法，将改进的Logistic 映射和坐标采样与关联成像相结合，减少了密文的传输量。Li等^[5]提出了基于计算关联成像（computational ghost imaging, CGI）和提升小波变换并结合异或操作的多图像加密方法，提高了加密系统的安全性。Wu等^[6]利用位置复用将测量的不同衍射距离的矢量强度进行迭加并与GI结合，提高了图像的传输效率。Zhang等^[7]将相位恢复与GI相结合，提出了一种基于相位恢复算法和GI的多图像全息加密技术。Sui等^[8]提出了一种基于强度方程传输的光学多图像认证方法，并应用强度方程传输技术实现了光学多图像认证。目前，这些多图像加密方法虽然提高了加密图像的数量，但也增加了系统的复杂性。同时，随着加密容量的增加，数据处理的时间和复杂度也随之增加。而且由于单一技术加密方法的局限性，这些多图像加密方法大多是基于多种技术手段的组合^[9-11]。

本文提出了一种级联相位检索算法（cascaded phase retrieval algorithm, CPRA）与CGI相结合的方法。它可以在没有串扰的情况下大大增强加密容量。首先，利用输入图像与目标图像之间的相位检索算法，对每个输入图像进行相位编码，得到目标图像与输入图像之间的相位掩膜对。然后，将得到的最终目标图像经一组调制散斑照明，并由桶探测器记录测量数据，获得最终的密文。该系统的密钥由CPRA加密阶段的两个相位掩膜、CPRA输出平面上的相位分布以及关联成像阶段的调制相位组成，大大提升了密钥空间的维度。本文通过仿真实验验证了加密系统的性能。实验结果表明，所提出的加密系统不仅具有较好的安全性，而且能够抵抗噪声、裁剪等攻击，同时也具有较好的加密容量。

1 基于CPRA与CGI的多图像加密原理 1.1 CPRA加密原理

基于单一的4f相关器的光学加密系统的体系结构如图1所示。输入输出平面和傅里叶平面的坐标分别用 $ \left( {x,y} \right) $ 和 $ \left( {u,v} \right) $ 表示； $ f\left( {x,y} \right) $ 表示待加密图像； $ g\left( {x,y} \right) $ 表示目标图像。随机相位掩码 $ RPM1 $ 和 $ RPM2 $ 分别写为 $\exp \left[ {{\rm{i}}\theta \left( {x,y} \right)} \right]$ 和 $\exp [ {\rm{i}}\varphi ( u, v ) ]$ 。因此，将 $ f\left( {x,y} \right) $ 加密到 $ g\left( {x,y} \right) $ 的过程实际上是一个寻找正确的相位掩膜对 $\exp \left[ {{\rm{i}}\theta \left( {x,y} \right)} \right]$ 和 $\exp \left[ {{\rm{i}}\varphi \left( {u,v} \right)} \right]$ 的问题。CPRA可以在待加密图像 $ f\left( {x,y} \right) $ 和目标图像 $ g\left( {x,y} \right) $ 之间的强度约束下进行搜索^[12]。CPRA由多个循环迭代组成，它受两个强度或振幅图像的约束，每个循环都涉及正向和反向传播。在第 $ k\left( {k = 1,2,3\cdots } \right) $ 步的前向迭代中，这两个相位分布函数可以描述为 $ \mathop \theta \nolimits^k \left( {x,y} \right) $ 和 $\varphi^k \left( {u,v} \right)$ 。因此，执行每次前向迭代后的输出图像可以写为

$ \begin{split} {{\displaystyle g}}^{k}\left(x,y\right) &\mathrm{exp}\left[{\rm{i}}{{\displaystyle \phi }}^{k}\left(x,y\right)\right]=\\ &{{\displaystyle IFT}}\left(FT\left\{f\left(x,y\right)\mathrm{exp}\left[{\rm{i}}{{\displaystyle \theta }}^{k}\left(x,y\right)\right]\right\}\times \right.\\ & \left.\mathrm{exp}\left[{\rm{i}}{{\displaystyle \varphi }}^{k}\left(u,v\right)\right]\right) \end{split} $

(1)

式中： $ FT $ 和 $\mathop {IFT}\nolimits$ 分别表示傅里叶变换和傅里叶反变换； $ \mathop g\nolimits^k \left( {x,y} \right) $ 是第 $ k $ 次迭代近似振幅图像； $ \mathop \phi \nolimits^k \left( {x,y} \right) $ 是输出平面上的第 $ k $ 次迭代相位； $ f\left( {x,y} \right) $ 是每个前向迭代中的约束条件，并保持不变。反向迭代中的图像也可以用同样的方式来处理。在这种情况下，在每次反向迭代中， $ \mathop g\nolimits^k \left( {x,y} \right) $ 被 $ g\left( {x,y} \right) $ 作为加密的目标图像。而这两个相位分布函数在每个循环中都可以被更新为

$ \begin{split} \varphi ^{k + 1} \left( {u,v} \right) =& angle\left\{ {\dfrac{{FT\left\{ {g\left( {x,y} \right)\exp \left[ {{\rm{i}} \phi ^k \left( {x,y} \right)} \right]} \right\}}}{{FT\left\{ {f\left( {x,y} \right)\exp \left[ {{\rm{i}}\theta ^k \left( {x,y} \right)} \right]} \right\}}}} \right\} \\ \theta ^{k + 1} \left( {x,y} \right) =& angle\left\{ { {IFT} \left( {FT\left\{ {g\left( {x,y} \right)\exp \left[ {{\rm{i}}\phi ^k \left( {x,y} \right)} \right]} \right\} \times}\right.}\right.\\ &\left.{\left.{\exp \left[ { - {\rm{i}} \varphi ^{k + 1} \left( {u,v} \right)} \right]} \right)} \right\} \\[-13pt] \end{split} $

(2)

式中， $ angle\left\{ \cdot \right\} $ 表示相位提取操作。在每个循环中，两个相位分布同时被修改。初始阶段 $ \theta^1 \left( {x,y} \right) $ 和 $ \varphi ^1 \left( {u,v} \right) $ 随机分布在区间 $\left[ {0,2\text{π} } \right]$ 上。对于不同的初始相分布，获得的相位对如式（2）所示，是不同的。该算法的迭代次数由 $ g\left( {x,y} \right) $ 和 $ g^k \left( {x,y} \right) $ 之间的相关系数（correlation coefficient, CC）来控制。CC可以表示为

$ CC = \dfrac{{{cov} \left[ {g\left( {x,y} \right),\mathop g\nolimits^k \left( {x,y} \right)} \right]}}{{\mathop \sigma \nolimits_g \mathop \sigma \nolimits_{\mathop g\nolimits^k } }} $

(3)

式中： $ {cov} \left[ {g\left( {x,y} \right),\mathop g\nolimits^k \left( {x,y} \right)} \right] $ 是 $ g\left( {x,y} \right) $ 和 $ \mathop g\nolimits^k \left( {x,y} \right) $ 之间的协方差； $ \mathop \sigma \nolimits_g $ 和 $ \mathop \sigma \nolimits_{\mathop g\nolimits^k } $ 是 $ g\left( {x,y} \right) $ 和 $ \mathop g\nolimits^k \left( {x,y} \right) $ 的标准差。当 $ g\left( {x,y} \right) $ 与 $ \mathop g\nolimits^k \left( {x,y} \right) $ 之间的CC值大于设定的阈值，则迭代过程将停止。为了得到加密的目标图像，对输入平面上的 $f\left( {x,y} \right)\exp \left[ {{\rm{i}}\mathop \theta \nolimits^k \left( {x,y} \right)} \right]$ 进行傅里叶变换，并乘以相位项 $\exp \left[ {{\rm{i}} \varphi ^k \left( {u,v} \right)} \right]$ ，然后对它们进行傅里叶反变换，最后进行绝对值运算，得到加密的目标图像。最终得到的振幅图像 $ \mathop g\nolimits^k \left( {x,y} \right) $ 近似于目标图像 $ g\left( {x,y} \right) $ 。根据参考文献 [13] 可知，在迭代次数达到10次后，便可以获得较为不错的效果。

1.2 CGI加密原理

GI作为一种非局域成像方式，在成像时，经过空间光调制器调制的光束先被分成参考光路和物体所在光路，然后经光学器件收集光学信息。参考光路和物体光路的光学信息经过二阶关联计算，可以实现对信息的重构^[14]。而CGI方法则简化了装置，省略了参考光路，见图2。

在CGI过程中，只需要一个桶探测器来获取密文 $ \mathop B\nolimits_i $ 。 $ \mathop B\nolimits_i $ 可表示为

$ \mathop B\nolimits_i = \displaystyle\int {T\left( {x,y} \right)} \mathop I\nolimits_i \left( {x_p,y_q} \right){\rm{d}}x{\rm{d}}y $

(4)

式中， $ T\left( {x,y} \right) $ 为要加密的图像。在加密过程中，将大小为 $ n \times n $ 的二维图像 $ T\left( {x,y} \right) $ 转化为一个一维列向量（ $ \mathop n\nolimits^2 \times 1 $ ）。由式（4）可以得到第 $ i $ 次的光强分布函数为 $ \mathop I\nolimits_i \left( {\mathop x\nolimits_p ,\mathop y\nolimits_q } \right) $ ，其中 $i = 1,2, \cdots, N$ ； $p,q = 1, 2, \cdots ,n$ 。光强分布函数的矩阵表达式为

$ \mathop {\boldsymbol{I}}\nolimits_i = \left | {\begin{array}{*{20}{c}} {\mathop I\nolimits_{11}^i }& \cdots &{\mathop I\nolimits_{1n}^i } \\ \vdots & & \vdots \\ {\mathop I\nolimits_{n1}^i }& \cdots &{\mathop I\nolimits_{nn}^i } \end{array}} \right| $

(5)

式中， $ \mathop I\nolimits_{nn}^i $ 是第 $ i $ 次测量矩阵中，第 $ n $ 行第 $ n $ 列的元素。光强分布矩阵的大小为 $ n \times n $ ，将矩阵拉伸成一个一维行向量，其大小为（ $ 1 \times \mathop n\nolimits^2 $ ），该向量可表示为

$ \begin{split} \mathop {\boldsymbol{I}}\nolimits_i =& \left[ \begin{array}{*{20}{c}} {\mathop I\nolimits_{11}^i } & {\mathop I\nolimits_{12}^i } & \cdots & {\mathop I\nolimits_{1n}^i } & {\mathop I\nolimits_{21}^i } & {\mathop I\nolimits_{22}^i } & \cdots \end{array} \right. \\ & \left. \begin{array}{*{20}{c}} & {\mathop I\nolimits_{nn - 1}^i } & {\mathop I\nolimits_{nn}^i} \end{array} \right] \end{split}$

(6)

但CGI的检测和重构时间较长。为解决这一问题，有学者提出一种压缩感知（compressed sensing, CS）的CGI算法。它利用压缩感知的特性，以较少的测量次数恢复出高质量的图像。当测量次数为 $ M $ （ $ M < \mathop n\nolimits^2 $ ）时，将生成 $ M $ 个大小为 $ n \times n $ 的一维行向量。由此得到一个大小为 $ M < \mathop n\nolimits^2 $ 测量矩阵，以此用作CGI的测量矩阵。桶探测器接收到的总强度 $ \mathop {\left\{ {\mathop B\nolimits_i } \right\}}\nolimits_{i = 1}^M $ 可以表示为：

$ \begin{split} \left[ {\begin{array}{*{20}{c}} {\mathop B\nolimits_1 } \\ {\begin{array}{*{20}{c}} {\mathop B\nolimits_2 } \\ {\begin{array}{*{20}{c}} \vdots \\ {\mathop B\nolimits_M } \end{array}} \end{array}} \end{array}} \right] =& \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\mathop I\nolimits_{11}^1 } \\ {\mathop I\nolimits_{11}^2 } \\ \vdots \\ {\mathop I\nolimits_{11}^M } \end{array}}&{\begin{array}{*{20}{c}} \cdots \\ \cdots \\ \\ \cdots \end{array}} \end{array}}&{\begin{array}{*{20}{c}} {\mathop I\nolimits_{1n}^1 } \\ {\mathop I\nolimits_{1n}^2 } \\ \vdots \\ {\mathop I\nolimits_{1n}^M } \end{array}}&{\begin{array}{*{20}{c}} \cdots \\ \cdots \\ \\ \cdots \end{array}}&{\begin{array}{*{20}{c}} {\mathop I\nolimits_{nn}^1 } \\ {\mathop I\nolimits_{nn}^2 } \\ \vdots \\ {\mathop I\nolimits_{nn}^M } \end{array}} \end{array}} \right] \cdot \\ &\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\mathop T\nolimits_{11} } \\ \vdots \end{array}} \\ {\mathop T\nolimits_{1n} } \\ \vdots \\ {\mathop T\nolimits_{nn} } \end{array}} \right]\\[-45pt] \end{split} $

(7)

在解密过程中，光强分布函数 $ {I_i}\left( {x,y} \right) $ 和探测得到的光强分布数据 $ \mathop B\nolimits_i $ 联合计算可解出秘密图像信息，可以表示为：

$ {T_{{\rm{CS}}}}\left( {x,y} \right) = \left\langle {\mathop B\nolimits_i \mathop I\nolimits_i \left( {x,y} \right)} \right\rangle - \left\langle {\mathop B\nolimits_i } \right\rangle \left\langle {\mathop I\nolimits_i \left( {x,y} \right)} \right\rangle $

(8)

式中， $ \left\langle \cdot \right\rangle $ 表示平均操作。

根据CS重建原始信号的原理，可以通过凸优化算法重建原始图像。

$ \mathop T\nolimits_{{\rm{CS}}} = \mathop T\nolimits' ,\arg \min \left\| {{\mathit{\Psi}} \left\{ {\mathop T\nolimits' \left( {x,y} \right)} \right\}} \right\| $

(9)

2 基于CPRA与CGI的多图像加解密流程

加密流程：

结合基于 4f 的 N个单图像加密方法对多图像进行编码，如图3所示。在该方法中，N个目标图像对应 N个输入图像，利用输入图像 $ \mathop g\nolimits_{0n}^k \left( {x,y} \right) $ 与目标图像 $ \mathop g\nolimits_{0n + 1} \left( {x,y} \right) $ 之间的相位检索算法，对每个输入图像编码不同的 $\exp \left( {{\rm{j}} \varphi_k } \right)$ 和 $\exp \left( {{\rm{j}} \theta _k } \right)$ 。通过使用这些相位掩码，可以对图像 $ \mathop g\nolimits_{0n}^k \left( {x,y} \right) $ 进行如下编码

$ g_{0n + 1}^k \exp \left( {{\rm{i}} \phi _n } \right) = IFT\left\{ {FT\left[ { g_{0n}^k \exp \left( {{\rm{i}} \theta _n } \right)} \right]\exp( {\rm{i }}\varphi _n) } \right\} $

(10)

式中：下标 $ n $ 表示图像的数量。由此输入图像 $ \mathop g\nolimits_{0n}^k $ 则是由 $\mathop g\nolimits_{0n}^k \exp \left( {{\rm{i}}\mathop \phi \nolimits_{n - 1} } \right)$ 经过光学加密系统加密所得， $\mathop g\nolimits_{0n + 1}^k \exp \left( {{\rm{i}}\mathop \phi \nolimits_n } \right)$ 则可以表示为

$ \begin{split} g_{0n + 1}^k \exp \left( {{\rm{i}} \phi _n } \right) =& IFT\left\{ {FT\left[ { g_{0n}^k \exp \left( {{\rm{i}} \phi _{n - 1} } \right)\exp \left( { - {\rm{i}} \phi _{n - 1} } \right) \times } \right.} \right.\\ & \left.{ \left.{ \exp \left( {{\rm{i}} \theta _n } \right)} \right]\exp \left( {{\rm{i}} \varphi _n } \right)} \right\}\\[-10pt] \end{split} $

(11)

进一步整理多图像加密的等式为

$ \begin{split} g_{0n + 1}^k \exp \left( {{\rm{i}} \phi _n } \right) =& IFT\left[ {FT\left\{ { g_{0n}^k \exp \left( {{\rm{i}} \phi _{n - 1} } \right) } \right.} \right.\times\\ &\left.{\left.{\exp \left[ {{\rm{i}}\left( { \theta _n - \phi _{n - 1} } \right)} \right]} \right\}} \right]\exp \left( {{\rm{i}} \varphi _n } \right) \end{split} $

(12)

将最后一张目标图像经一组调制散斑照明，并由桶探测器记录测量数据，过程见公式（7）。将获取的密文 $ \mathop B\nolimits_i $ 、每对输入图像与目标图像之间相位掩码 $\exp \left( {{\rm{j}} \varphi _k } \right)$ 和 $\exp \left( {{\rm{j}}\mathop \theta \nolimits_k } \right)$ 以及光强分布函数 $ {I_i}\left( {x,y} \right) $ 传输给接收方，加密过程就完成了。

解密流程：

在解密过程中，光强分布函数 $ {I_i}\left( {x,y} \right) $ 和探测得到的光强分布数据 $ \mathop B\nolimits_i $ 联合计算可解出秘密图像信息，具体过程可以表示为公式（8），得到初级加密后的最后一张目标图像。多图像编码方法的解密阶段计算为

$ \begin{split} g_{0n}^k \exp \left( {{\rm{i}} \phi _{n - 1} } \right) =& IFT\left[ {FT\left\{ { g_{0n + 1}^k \exp \left( {{\rm{i}} \phi _n } \right) \times} \right.} \right.\\ & \left.{ \left.{\exp \left( { - {\rm{i}} \varphi _n } \right)} \right\}} \right] \times \exp \left[ { - {\rm{i}}\left( { \theta _n - \theta _{n - 1} } \right)} \right] \end{split} $

(13)

由此可以依次解密出第 $N,N - 1,N - 2, \cdots ,1$ 张明文图像。

3 仿真与分析

利用CPRA与CGI相结合的方法，将 4 张明文图像加密成最终的桶探测器值，使用正确密钥解密结果如图4所示。图4（a）为输入图像，图4（b）~（d）为目标图像，图4（e）~（h）为密钥完全正确情况下解密结果。解密图像与原输入图像间 CC值分别为0.9975，0.9932，0.9951，0.9967。结果表明，该系统解密图像的质量非常好。

为了验证该加密系统的级联特性，在图5中给出不同阶次相位掩码错误时的解密情况。图5（a）~（d）显示了CPRA解密阶段最后一个编码图像的相位密钥错误时的解码图像；图5（e）~（h）显示了二阶编码图像的相位掩码错误时对应的解密图像；图5（i）~（l）显示了在第一阶编码的图像的相位掩码错误时对应的解密图像。从图中可以看出，如果在解密的第一阶段使用了错误的密钥，则无法访问图像。这是由于其级联特性，各用户在接收端相互依赖。在解密前一个图像之前，当前的图像将不会被解密，因为输入的图像是连续加密的。

在信息的加密和传输过程中，噪声的影响是不可避免的。因此，一个合格的加密系统需要具有一定的鲁棒性。引入了标准差（ $ \sigma $ ）分别为0.1，0.2，0.3时的高斯噪声来模拟由光学器件引起的热噪声。

从表 1可以看出：（1）对于高斯噪声，随着标准差的增加，重构效果呈下降趋势；（2）当噪声标准差达到0.3时，解密图像与原始图像之间的 CC保持在0.48以上，仍能得到图像的部分信息。结果表明，系统中光学器件产生的热噪声对图像的解密质量有较大的影响，但该系统具有一定的抗噪声攻击能力。

对密文在不同位置进行不同程度的裁剪，选用 CC作为评价指标来验证加密系统的抗裁剪攻击性能。由表 2和图6可知：（1）随着密钥裁剪面积的增加，重构图像与原始图像之间的 CC 值逐渐减小，重构图像的清晰度越低，越难分辨出原始图像的信息；（2）当裁剪比例为30%时，重构出的图像与原始图像的CC值均在0.78以上，能够大致地分辨出原始图像的细节信息，重构的效果较好；（3）当密钥的裁剪比例为50%时，重构出的图像与原始图像的 CC值均大于0.50，仍可分辨出加密图像的部分信息，说明该方法能够较好地抵抗裁剪攻击。综上可得：本文所提出的加密系统不仅对噪声攻击表现出良好的鲁棒性，而且能够抵抗裁剪攻击，进一步表明该加密系统具有较好的鲁棒性。

4 结 论

本文提出了一种基于CPRA和CGI的多图像加密方法，该方法可以通过纯光学方法实现。与传统的CPRA方法相比，系统中增加了CGI加密过程，不仅增加了密钥空间，提高了整个系统的安全性，而且密文的形式为一系列桶探测器值，利于传输与存储，由此提高了整个加密系统的抗裁剪攻击能力。该方法拓展了CGI加密的容量，具有更好的实用性和应用前景。