光学相位恢复算法

前言

经典的光学相位恢复算法有 1971 年提出的 GS 算法和 1980 年提出的 YG 算法，分别适用于特定光学衍射过程和一般光学衍射过程

Gerchberg-Saxton Algorithm

基本定义

二维光场在物平面的入射光场可表示为
$$E_1(x_1,y_1)={\rho }_1(x_1,y_1)\cdot e^{ı{\varphi }_1(x_1,y_1)}$$
同理，像平面的出射光场可表示为
$$E_2(x_2,y_2)={\rho }_2(x_2,y_2)\cdot e^{ı{\varphi }_2(x_2,y_2)}$$
其中 ${\rho }_1(x_1,y_1),{\rho }_2(x_2,y_2)$ 为物平面和像平面光场的振幅分布，${\varphi }_1(x_1,y_1),{\varphi }_2(x_2,y_2)$ 为物平面和像平面光场的相位分布

在夫琅和费衍射或透镜衍射中，入射光场与出射光场之间存在傅立叶变换关系

正过程可以近似表示为
$$E_2(x_2,y_2)=\mathcal{F}[E_1(x_1,y_1)]$$
其逆过程为
$$E_1(x_1,y_1)={\mathcal{F}}^{-1}[E_2(x_2,y_2)]$$
在 GS 算法中振幅分布 ${\rho }_1(x_1,y_1),{\rho }_2(x_2,y_2)$ 已知，通过多次迭代的方法求解未知相位分布 ${\varphi }_1(x_1,y_1),{\varphi }_2(x_2,y_2)$，初值敏感，易陷入局部最优，后续有诸多改进方法如加权 GS 算法等

算法逻辑

1. 随机初始化入射光场的相位 ${\varphi }_1^{\prime }(x_1,y_1)$ ，与已知的入射光场振幅分布 ${\rho }_1(x_1,y_1)$ 合成入射光场 $E_1^{\prime }(x_1,y_1)$

2. 入射光场经过夫琅和费衍射或透镜汇聚过程等价通过傅立叶变换得射光场 $E_2(x_2,y_2)$ 的相位 ${\varphi }_2(x_2,y_2)$

3. 由上步的新相位和已知出射光场振幅分布 ${\rho }_2(x_2,y_2)$ 合成出射光场 $E_2^{\prime }(x_2,y_2)$

4. 出射光经过夫琅和费衍射或透镜汇聚的逆过程等价通过傅立叶反变换得入射光场 $E_1(x_1,y_1)$ 的相位 ${\varphi }_1(x_1,y_1)$

5. 由上步的新相位和已知入射光场振幅分布 ${\rho }_1(x_1,y_1)$ 合成入射光场 $E_1^{\prime }(x_1,y_1)$

6. 回到步骤 2 进行下一次循环，不断迭代直至光场趋于稳定 $min\frac{1}{k}{\sum }_{x_1,x_2,y_1,y_2=1}^{x_1,x_2,y_1,y_2=k}[\mid E_1(x_1,y_1){\mid }^2-\mid E_1^{\prime }(x_1,y_1){\mid }^2{]}^2$，$k$ 为变量的上限值

杨顾算法

基本定义

YG 算法扩展了 GS 算法，解决入射光场与出射光场之间非幺正变换的情况，可普遍应用于一般衍射器件，后续有诸多改进方法如加权 YG 算法等

高维光场进行 $m$ 次均匀采样，因此物平面的入射光场可表示为
$$E_{1i}={\rho }_{1i}\cdot e^{ı{\varphi }_{1i}},i=1,2,\cdots ,m$$
在传播轴上进行 $n$ 次均匀采样，因此像平面的出射光场可表示为
$$E_{2j}={\rho }_{2j}\cdot e^{ı{\varphi }_{2j}},j=1,2,\cdots ,n$$

其中 ${\rho }_1,{\rho }_2$ 和 ${\varphi }_1,{\varphi }_2$ 分别是物平面和像平面光场的振幅和相位

物平面的入射光场和像平面的出射光场之间满足一般光学变换关系，可以是非幺正变换即变换矩阵可以是非酉矩阵

酉矩阵定义，设存在一个 $n$ 阶矩复数矩阵 $U$ 满足
$$U^{H}U=U{U}^H=E_n$$
则 $U$ 为酉矩阵，其中 $U^H$ 是 $U$ 的共轭转置矩阵且 $U^H=U^{-1}$，$E_n$ 是 $n$ 阶单位矩阵，酉矩阵 $U$ 行与列的二阶范数都是 $1$ 且 $\parallel Ux{\parallel }_2=\parallel x{\parallel }_2$，当元素都为实数是就是正交矩阵

正过程为
$$E_{2j}=\sum _{i=1}^m{G}_{ij}\cdot E_{1i}$$

其中 $G$ 是核函数根据描述光场传播的衍射公式以及物平面与像平面之间相互关系决定

定义所得光场和所需光场之间的二阶范数来描述两者的距离
$$D^2=\parallel E_2-G\cdot E_1{\parallel }_2$$

当 $D^2$ 取极小值时，物平面的入射光场的相位 ${\varphi }_1$ 与像平面出射光场的相位 ${\varphi }_2$ 可表示为
$$\begin{array}{rl}{\varphi }_{1k}& =arg\lfloor \sum _{j=1}^n{G}_{jk}^H{E}_{2j}-\sum _{j\ne k}{A}_{kj}{E}_{1k}\rfloor \end{array}\text{\hspace{1em}(1)}$$

$$\begin{array}{rl}{\varphi }_{2k}& =arg\lfloor \sum _{k=1}^m{G}_{kj}{E}_{1k}\rfloor \end{array}\text{\hspace{1em}(2)}$$

其中 $A=G^{+}G,A^H=A$ 而 $G^{+}$ 是 $G$ 的广义逆矩阵 Moore-Penrose 伪逆，$\lfloor f(x)\rfloor$ 是对 $f(x)$ 向下取整

Moore-Penrose 伪逆定义，设存在一个矩阵 $A$ 有
$$A^{+}=\underset{x\to 0}{\lim }(A^{T}A+xI{)}^{-1}{A}^T$$
计算式
$$A^{+}=V^{+}\Sigma U^H$$
先对 $A$ 采用奇异值分解得到
$$\mathrm{S}\mathrm{V}\mathrm{D}(A)=U\Sigma V^H$$
$U,V^H$ 是酉矩阵，$\Sigma$ 是一个除了主对角线上的元素以外全为 $0$ 的矩阵，主对角线上的每个元素都是奇异值，若 $A$ 为方阵则 $\Sigma$ 即为特征值矩阵 $D$

对于 $Ax=y$ 如果无解我们使用 $\parallel y{A}^{+}-x{\parallel }_2$ 最小的 Moore-Penrose 伪逆，如果多解我们使用 $\parallel y{A}^{+}{\parallel }_2$ 最小的 Moore-Penrose 伪逆

算法逻辑

1. 随机初始化入射光场的相位 ${\varphi }_1$ ，与已知的入射光场振幅分布 ${\rho }_1$ 合成初始化的入射光场 $E_1$

2. 入射光场经过正过程通过 $(2)$ 得出射光场像平面 $E_2$ 的相位 ${\varphi }_2$

3. 由上步的新相位和已知出射光场振幅分布 ${\rho }_2$ 合成出射光场 $E_2$

4. 出射光场通过 $(1)$ 得出入射光场 $E_1$ 更为准确的相位 ${\varphi }_1$

5. 由上步的新相位和已知入射光场振幅分布 ${\rho }_1$ 合成入射光场 $E_1$

6. 回到步骤 2 进行下一次循环，不断迭代直至光场趋于稳定是时终止