Effect of amplitude freedom on imaging of diffractive optical elements

衍射光学元件(diffractive optical elements, DOE)^[1]可对光波进行相位和振幅调制，在光束整形、衍射成像等领域有广泛的应用前景。DOE的设计与光学变换中的相位恢复相似，最早的相位恢复算法是1972年由GERCHBERG与SAXTON提出的GS算法^[2]，它本质上是一种最速下降算法，容易陷入局部最优解，且再现像光强分布均匀性较低。目前主要从两方面对GS算法进行改进，一是改进算法本质，如杨-顾算法^[3]、桑涛算法^[4]和遗传算法^[5]等。二是在迭代过程中优化振幅和相位。常用的振幅优化方法有：增加振幅自由度^[6]、自适应加法(adaptive-additive，AA)^[7]、利用3个面(输入面和两个输出面)的强度信息恢复相位^[8]、用实际振幅与迭代后振幅加权替换原振幅^[9]，以及引入梯度下降和加权^[10]等。相位优化方法有：利用迭代中某两次相位值加权得到下次计算的相位值^[11]、微分中值定理替换相位^[12]、二系数误差扩散滤波^[13]，还有限制相位分布抑制散斑^[14]等。

作者从光波的传播特性入手：由于衍射过程一定存在泄漏，则目标图像的背景区域(光场振幅一般设为零)存在一定量的光强分布，最简单的方法是为目标图像增加一个均匀振幅自由度，从而使设计过程更符合实际衍射情况。再考虑到衍射实际上是低通滤波的过程，图像衍射后振幅分布不再均匀，背景区域振幅实际呈减幅振荡分布，所以在为目标图像增加均匀振幅(uniform amplitude)分布的基础上，应该再构建减幅振荡的振幅分布，即为其叠加混合振幅(mixed amplitude)才更符合实际光场分布。

GS算法和模拟退火算法(simulated annealing algorithm，SA)^[15]可用于光束整形，将高斯光束整形为光强均匀分布的平顶光束和环形光束等^[16-19]。用这两种算法设计了傅里叶变换型DOE，发现添加振幅自由度后取得了较好效果，而且在背景振幅相同的情况下，叠加混合振幅得到的再现像光强均匀性更高。

分别用GS算法和模拟退火算法设计傅里叶型DOE，如图 1所示。已知输出的衍射像面Σ₂上目标图像的振幅分布为B_uv，(u, v)为输出面点的坐标，且DOE为相位型元件，即输入的物面Σ₁上振幅分布A_xy=1，(x, y)为输入面点的坐标，FT为傅里叶变换(Fourier transform)，IFT为傅里叶逆变换(inverse Fourier transform)，设计结果为Σ₁面上的相位分布。

Figure 1.  DOE's design schematic

GS算法迭代过程如下：(1)生成随机相位分布φ₀，与已知振幅A_xy(A_xy=1)组成初始物面分布，即f=A_xyexp(jφ₀)；(2)将复振幅分布f代入(1)式进行傅里叶变换，得到衍射像分布F；(3)用目标图像振幅B_uv代替F中的振幅，与F中的相位组成新的衍射像分布F′，即F′=B_uvexp[jangle(F)]，angle()代表计算相位；(4)将F′代入(2)式进行逆傅里叶变换，得到物面分布f′，用振幅1代替f′中振幅，与其相位组成新的物面分布f，即f=1exp[jangle(F')]；(5)回到第(2)步，进行下一次循环，直到误差函数达到给定精度ε，或循环次数达到阈值停止迭代，即可得到所需的DOE相位分布φ=angle(f′) ^[20]:

模拟退火算法作为局部搜索算法的扩展，其核心思想是在每一次修改模型的过程中，在最优状态下的邻域内随机产生1个新的状态模型，然后以一定的概率选择邻域中评价函数更好的状态，这种接受新模型的方式使其成为一种全局最优算法。

模拟退火算法的迭代过程如下：(1)给定初始相位φ₀，初始温度t₀，最低温度t_end；(2)对当前相位φ₀做邻域扰动，生成新相位解：φ₁=φ₀+S(x, y)，S(x, y)为邻域扰动因子；(3)初始相位φ₀和新相位解φ₁分别与给定振幅A_xy(A_xy=1)组合，计算衍射像分布F₀和F₁；(4)计算F₀和F₁与目标图像振幅B_uv间的均方(mean square, MS)误差E_{MS, 0}和E_{MS, 1}；(5)判断两个误差的大小，若E_{MS, 1} < E_{MS, 0}，用衍射像分布F₁逆运算得到DOE分布f₁，否则重复步骤(2)~步骤(5)，直到达到当前温度下的最大循环；(6)降低温度进行下一次迭代，直到达到最低温度t_end，输出优化后的相位，即所求相位^[15]。

模拟中使用均方误差E_MS判断输出面振幅与目标图像振幅间的误差，它常用来判断设计结果是否满足要求，E_MS定义如下：

DOE衍射像期望区域的光强不均匀性I_RMS(均方根root mean square, RMS)与振幅不均匀性A_RMS是DOE设计的另一重要评价标准，它们的定义如下：

(3) 式~(5)式中，(u, v)为输出面Σ₂上采样点的坐标; U, V为行与列的采样点数，计算中均取512;I_uv和A_uv分别为Σ₂上期望区域点(u, v)处的光强和振幅; I和A为期望区域的平均光强与平均振幅。

用标量衍射理论来阐明改进方法的原理。1678年, 惠更斯为了描述波的传播过程提出关于子波的设想; 1818年, 菲涅耳引入干涉概念补充了惠更斯原理，提出了计算衍射光场分布U(P)的方法，如图 2所示。其表达式为：

Figure 2.  Diffraction of an aperture

式中，S为积分曲面，C为常数，U(P₀)为衍射面上任一点P₀处光场的复振幅，θ为PP₀与过P₀点波面的法线n之间的夹角，K(θ)表示子波源P₀对P点光场的贡献，它与角度θ有关，k是波数，r是P与P₀之间的距离。但此表达式仅基于子波源和相干叠加的假说，缺乏严格的波动理论基础支撑，C与K(θ)的具体形式也无法确定。

1882年, 基尔霍夫从亥姆霍兹方程出发，利用格林定理求解波动方程，导出了基尔霍夫衍射公式：

式中, 格林函数G=exp(jkr)/r，推导中用到两点假设：(1)在衍射孔径Σ上(如图 2所示)场分布U(P₀)及其偏导数$\partial$U(P₀)/$\partial$n与没有衍射屏时完全相同；(2)在位于衍射屏几何阴影区S₁的那部分上，场分布U及其偏导数$\partial$U(P₀)/$\partial$n恒为0。以此为基础推导出了著名的菲涅耳-基尔霍夫衍射公式：

式中, λ是光波波长，λ=2π/k。正是这两点假设备受质疑^[21]，因为对于波动方程，如果它的一个解在任意非无限小的面上，光场的复振幅和它的法向导数都为零，则这个解在整个空间都为零。直接与这个基本假设相矛盾。为此索末菲选择了与基尔霍夫不同的格林函数G，从而在边界条件中不必规定屏后光场及其法向导数均为零，克服了基尔霍夫边界条件的矛盾。

实际情况是，衍射屏会在一定程度上影响孔径附近的场分布，屏后的阴影区光场分布不会完全为零，光场会泄漏到屏后孔径之外^[21]。比如常用的直边菲涅耳衍射，其衍射图样的强度呈减幅振荡分布，不存在光强为零的区域。所以在空间中任意处有光源存在，便不会有光强分布为0的区域，反之亦然。于是，在衍射光学元件设计中必须考虑到这一特性，继而作相应的改进：目标图像的背景不应该为零，必须增加一个振幅自由度。最简单的方法是将目标图像的背景振幅设置为不为零的均匀振幅分布。但注意到衍射是一个低通滤波的过程，通常会引起如直边衍射的振铃现象，所以增加的振幅分布取成减幅振荡与均匀振幅的混合会更合理。

模拟设计参量如下:DOE与目标图像大小均设为512pixel×512pixel，衍射过程借助快速傅里叶变换完成。所用计算机为安装有Windows7的惠普电脑，处理器3.2GHz，内存4GB，MATLAB为R2014b版64位软件。设计DOE使用的目标图像可分为期望的信号窗口区域与背景部分，通常信号窗口振幅值设为1，背景区域振幅分布为0，图 3分别是目标像光强的2-D与3-D分布。先用GS算法和模拟退火算法对背景振幅分布为0的目标图像完成DOE设计，结果如图 4所示，分别是GS算法和SA算法的结果。模拟计算得到的GS算法I_RMS=13.4%，E_MS=9.8%，SA算法I_RMS=14.4%，E_MS=9.8%，信号窗口的光强分布不均匀，设计结果不理想。

Figure 3.  Target image intensity distribution

Figure 4.  Reconstructed images intensity distribution by GS algorithm and SA algorithm

接着为目标图像的背景分别叠加均匀振幅和混合振幅，再用GS算法和SA算法进行验证。获得减幅振荡背景的过程为：对目标图像“光”字作傅里叶变换得到频谱，接着对频谱进行理想低通滤波(滤波窗为80pixel×80pixel的矩函数)，做傅里叶逆变换得到背景区域有振铃现象的“光”字，再将信号窗的振幅替换为原本的振幅值1，得到目标图像，这时的背景区域呈减幅振荡分布。这样做的原因为：衍射其实是一个低通滤波的过程，在频域中使用矩形理想低通滤波器，相当于在空域中将原始图像与sinc函数作卷积，由于sinc函数是振荡变化的，卷积计算后原始图像上会叠加上振荡信号，即产生振铃效应。经验证，滤波窗的大小对设计结果没有影响。实际计算中，将背景振幅分别设为信号窗振幅0.03, 0.05, 0.08, 0.1, 0.2和0.3倍的均匀值，再在此振幅上叠加减幅振荡背景分布。

背景振幅选取为0~0.3基于如下考虑：衍射造成光泄露，其光强肯定弱于照射光强。以典型的直边衍射为例，如果按几何光学而不考虑衍射，观察屏上与直边边缘对应的位置处光强为零；但实际有衍射发生，同一位置处光强为照射光强的0.25倍(且为泄露光强的最大值)，换算成振幅就是0.5倍，这个比例关系与直边衍射的距离和照明光波波长均无关。当然，衍射屏以及衍射类型不同时，泄露光强的最大值与照明光强的比例关系会有所不同。但以直边衍射为参考，均匀背景取值最大值可以达到0.5。实际计算发现，当背景振幅最大值超过0.3时均匀性不再继续下降，反而开始变差，所以将背景振幅值选取为0~0.3。

设计中所用背景增加均匀振幅的目标图像如图 5a所示，背景增加混合振幅的目标图像如图 6a所示，此时的背景振幅为0.08，图 5b和图 6b分别为叠加均匀背景和混合背景时目标像254行归一化振幅分布。

Figure 5.  Target image with uniform amplitude 0.08 and its amplitude distribution of 254 rows

Figure 6.  Target image with mixed amplitude 0.08 and its amplitude distribution of 254 rows

首先在GS算法中进行设计，计算中有一点需要注意，相同振幅下每次计算使用的初始随机相位是固定的，选用当前振幅情况下成像效果最好时的相位，单次设计过程所花时长为41s。设计得到DOE，再现像第254行(见图 5、图 6中剖线位置)归一化光强分布如图 7所示。图 7a是叠加均匀背景的7组结果，图 7b是叠加混合背景的7组结果。选取背景振幅值为0.08时的再现像光强做3-D分布对比，如图 8所示。其中图 8a是叠加均匀背景的结果，图 8b是叠加混合背景的结果。增加均匀振幅和混合振幅情况下，光强不均匀性对比见图 9a, 振幅不均匀性对比见图 9b。

Figure 7.  Contrast of light intensity distribution of reconstructed image in GS algorithm

Figure 8.  Comparisons of reconstructed images in GS algorithm

Figure 9.  I_RMS and A_RMS the uniform amplitude and the mixed amplitude in GS algorithm

从图 8可以看到, 此时的信号窗内，光强分布较均匀，噪声较小，成像质量得到了改善，且叠加混合背景振幅的方法，比只加均匀背景振幅得到的再现像光强分布更均匀。从图 7、图 9可以看到，叠加均匀振幅和混合振幅时的再现像光强不均匀性都从14%左右降低到了1%以下，振幅不均匀性均从7%左右降到了0.5%以下，为目标图像增加均匀振幅和混合振幅的方法，确实有效降低了再现像光强不均匀性I_RMS和振幅不均匀性A_RMS。

接下来验证在模拟退火算法中使用两种方法的效果。目标图像选取同前。单次设计过程所用时长约82s。设计得到DOE的再现像第254行归一化光强分布, 如图 10所示。图 10a是叠加均匀背景的7组结果，图 10b是叠加混合背景的7组结果。选取背景振幅值为0.08时的再现像光强3维分布对比，如图 11所示，其中图 11a是叠加均匀背景的结果，图 11b是叠加混合背景的结果。增加均匀振幅和混合振幅情况下，光强不均匀性对比如图 12a所示，振幅不均匀性对比如图 12b所示。

Figure 10.  Contrast of light intensity distribution of reconstructed image in SA algorithm

Figure 11.  Comparisons of reconstructed images in SA algorithm

Figure 12.  I_RMS and A_RMS comparison of the uniform amplitude and the mixed amplitude in SA algorithm

从图 11可以看到，DOE再现像的光强分布均匀了很多，且在背景振幅相同情况下，混合振幅比均匀振幅效果更好。从图 10和图 12中可以看到，在模拟退火算法中，叠加均匀振幅和混合振幅时的再现像光强不均匀性也均从14%左右降低到了1%以下，振幅不均匀性均从7%左右降到了0.5%以下，通过改变目标图像的背景振幅值，也可大幅度降低再现像光强不均匀性和振幅不均匀性，且也是混合背景振幅情况下降低的速度更快。

至此验证了添加均匀背景振幅与混合背景在GS算法和SA算法中的有效性，下面将GS算法和SA算法，分别加均匀背景振幅与加混合背景的均方误差E_MS做个对比，结果如图 13所示。

Figure 13.  E_MS comparison of the uniform amplitude and the mixed amplitude

由图 13可以看到，GS算法(见图 13a)和SA算法(见图 13b)中均方误差都从10%左右降到了0.1%以下，叠加均匀振幅和混合振幅的两种方法对于降噪均有明显效果，但为目标图像叠加混合背景的方法效果更明显，误差E_MS降低速度更快。

图 14分别是GS算法中为目标图像增加均匀振幅和混合振幅后，再现像第254行归一化光强分布。实线为背景振幅为0时的光强，虚线为背景振幅为0.2时的光强。从图中可以看出，虽中间部分的信号窗光强均匀性大幅提高，但两边部分即背景区域的噪声也变大很多，且背景振幅为0.3时噪声变大会更明显，SA算法中与GS算法基本相同。这是因为叠加上的均匀振幅以及减幅振荡背景严格意义上都是不希望存在的噪声，背景越大，偏离目标图像越大，因而增大了误差。GS算法和SA算法中叠加混合背景的方法在背景振幅为0.05时就有了明显降噪效果，且未增大再现像背景噪声，显然增加混合背景的方法在降低不均匀性、提高信噪比方面更有优势。

Figure 14.  Intensity comparison of the uniform amplitude and the mixed amplitude in GS algorithm

由于光波在传播过程中有一定程度的泄漏，不存在振幅和法向导数完全为零的区域，且衍射后的图像边缘光强为减幅振荡分布，综合这两点，为目标图像增加了混合型振幅自由度，用GS算法和模拟退火算法设计傅里叶变换型DOE, 并进行了验证。模拟仿真结果表明：增加振幅自由度后，元件再现像光强不均匀性从不增加时的14%左右降低到1%以下，均方误差从10%左右降到0.1%以下。该方法对两种算法均适用。