叠层成像相位重建的原理分析与流程可视化

作为一种利用衍射强度的测量数据重建样品空间结构的无透镜成像技术, CDI 技术综合了衍射手段、采样理论和相位恢复算法，基本概念于1969年提出[1-3]。Ptychography这个单词是由 Hoppe [1-3]提出来的，以表示通过卷积定理来解决相位恢复问题的方法，最初用于解决 X 射线和电子束成像中的晶体测量问题。Ptycho来自希腊语“πτυξ”，意思是“重叠”，因此 Ptychography算法的特点之一就是相邻两次扫描位置的扫描区域部分重叠，能够实现全场成像和避免使用支撑域。CDI实验于1999年在同步辐射 (synchrotron radiation source, SRS)X 射线上得以首次实现[4]。

叠层扫描光路模型

2004 年, Rodenburg 和 Faulkner[5] 追溯了 Hoppe 早期的工作, 提出叠层扫描成像迭代引擎 (ptychographic iterative engine, PIE) 算法, 有效地克服了孤立样本限制和像模糊问题。

求解样品相位的过程就是利用非凸交替投影在样品面重叠约束集S和观察面振幅约束集M之间反复进行投影，找到两复希尔伯特空间的交集作为最终的解。关于复希尔伯特空间的知识可以看连接“https://www.zhihu.com/question/19967778/answer/2572975608

”。

L是复希尔伯特空间，I_\bar{D} 是D补集的特征函数，S是衍射出射面重叠约束，M是观察面振幅约束，通过非凸交替投影迭代求得交集x ，更详细的数学细节请看文献[6]和书籍“Topics in Fixed Point Theory”链接“https://link.springer.com/book/10.1007/978-3-319-01586-6”。

自由空间中的平面波

α,  β ,γ 分别是k 矢量与轴 x , y , z的夹角，角谱传播满足亥姆赫兹方程

解得

利用空间频率

角谱传播方程写为

简写为

其中定义传递函数为

衍射平面的透过率角谱函数

入射光的角谱函数

出射面光波角谱

观察面角谱函数

如果入射光的角谱是单位δ(fx, fy )函数，则

取傍轴近似

传递函数为

观察面光场的傅里叶积分

对于观察面角谱进行傅里叶逆变换

对应有空间域上的积分

其中脉冲响应函数

出射面光场复振幅

再次观察积分式子

计算脉冲响应函数

此处我们看到清晰的对应关系。衍射积分还可以表示为另一种傅里叶变换表达式

利用傅里叶移位定理，可以将图像傅里叶变换后的频域乘上相位再经傅里叶逆变换实现空域上的位移。

首先做变量代换

代入傅里叶变换

即

来看几个简单例子，假设有N1xN2像素的图像，在频率域上选定（s1,s2）的坐标点集作为平移的位置索引。

首先利用matlab生成坐标网格，左上角是负的右下角是正的

[s1,s2]=meshgrid(linspace(-100,100,6),linspace(-100, 100,6))

如下图所示将光斑图像进行二维傅里叶变换后乘以带有位置索引的相位，若坐标R=（100，100），分别选取红点坐标（-100，-100）和蓝点坐标（100，100）进行平移。

指数相位中的参数（b1,b2）会带着负号直接加到空域坐标上去，所以对于红点是O(r-R-（-R）)=O(r)，图像向左上角平移刚好回到中心，对于蓝点是O(r-R)，图像向右下角平移。再如下图，取(-R)=（-100，-100）坐标，将图像平移到左上角。

就如第三节的扫描位置数量，探针光束会依次在样品表面的6x6=36个位置进行照射，透过（或是反射）的光束，携带有光束和样品表面的振幅以及相位信息，经过z2距离的菲涅尔衍射传播后被CCD探测器所接收。注意探测器只能接收到光场的模值平方数据。由样品面的重叠区域给出相位约束，由观察面的测量值给出振幅约束，两约束集都是复希尔伯特空间，还原样品表面的相位信息则需要通过非凸交替投影在两个集合之间进行迭代来求解。每个位置都要进行下图中的一次循环，36个位置全部循环过后称为一次完整迭代，一般要经过多次完整迭代才能让图像变得更清晰。

移位函数定义

重建的迭代流程图

在上图中，(1)(2)为初始的样品和探针场函数，初始样品函数可以设置为0矩阵，由于探针处在中心位置，需基于傅里叶移位定理用 imshift() 函数将图像需要重建的区域移动到中心，并与探针相乘作为出射面的场函数U0。

正向传播过程，场函数U0经过fft2变换后得到角谱A0，利用角谱传播方程A0H=A将A0乘以传递函数H得到观察面角谱，最后进行ifft2就得到了观察面场函数分布U。这个过程可以定义为传播函数 propagate，即输入U0输出U。

在观察面上获取U的相位即可，振幅用测量值替换，这是重叠约束带给我们的支撑信息，由于要求两个圆形区域有60%以上的重叠率，故每次扫描到的新区域都带有部分上一次位置的相位信息，再通过振幅的约束，更新了采集区域，不断的强化这一过程。

逆向传播过程，右乘传递函数的共轭得到反向角谱传播方程，同样的经过变换后得到新的出射面场函数U0。

利用小框中的更新函数可以分解出样品和探针的函数。

定义正向移动光斑的函数imshift(）为f+，反向平移光斑时需要在坐标上面加负号，所以记为f-。由于第一次分解前是空白的，便于描述我们从第二次开始描述。如下图所示，第二次循环重建开始，红圈处是需要重建的区域，注意黑色光斑的中心区域坐标是Rj，红斑虚线框区域的中心坐标则是第二个需要重建的位置（你可暂且认为它是R2）。将红色虚线框用f+移动到中心后与同样处在中心探针相乘得到出射面场函数，再经过正propagate振幅替换逆propagate得到新出射面场函数，做差后的图像再反向平移到原来红色虚线框的位置，根据更新函数的形式还要将差值图像乘以同样平移到红色虚线框位置的探针，得到新增加的重建图像，最后叠加到上一次得到的重建图像上作为下一次重建区域选取的对象，这样虚线框的位置就重建好了。继续进行剩余的34个位置的重建，待36个位置全部重建完成后，完成一次完整循环，接着重复进行完整循环，位置不必按顺序来进行选取，可以在36个位置上随机进行挑选。

第二次迭代时光斑随流程的变化图

下图中左边是原始样品，右边是经过10次完整循环重建得到的重建图。源程序代码可参考清华大学高云辉的github主页“https://github.com/Yunhui-Gao/ptychographic-phase-retrieva”。

原图                         重建图

随着投影迭代的次数的增加，图像上的每个像素点的复值都逐渐收敛到了两个约束集合的交集上，使得图像达到整体结构收敛。下面的动图给出了图像上像素点在复平面上的迭代收敛过程，能够大致估计出结构收敛的程度大小。

图像复值像素点迭代动图