潘新宇 毕筱雪 董政 耿直 徐晗 张一 董宇辉 张承龙†

1) (中国科学院高能物理研究所,北京同步辐射装置,北京 100049)

2) (中国科学院大学核科学与技术学院,北京 100049)

3) (中国科学院高能物理研究所,散裂中子源科学中心,东莞 523808)

随着同步辐射技术的发展和光源相干性的提升,叠层相干衍射成像(ptychography)得到快速发展.叠层相干衍射成像算法解决了传统相干衍射成像算法收敛速度较慢、容易陷入局部最优解和算法停滞等问题,具有成像视场大、算法鲁棒性强、对误差容忍性高、应用范围广等优点,正成为相干衍射成像领域的热点研究方向.本文首先介绍了叠层相干衍射成像算法提出的背景;然后详细总结了叠层相干衍射成像算法的发展脉络、主要的算法流程以及应用场景,并且介绍了叠层相干衍射成像与人工智能结合的新算法及应用潜力;最后介绍了叠层相干衍射成像算法具体的并行化实现及常用软件包.本文有助于建立叠层相干衍射成像领域算法本身、人工智能以及计算方法全局研究视角,对于促进叠层相干衍射成像方法学的系统发展具有重要的参考意义.

1 引言

传统的透射成像方法可以记录样品的吸收衬度图像,但是对于弱吸收样品仅重构样品的吸收衬度往往效果不佳,需要考虑相位衬度成像.相干衍射成像(coherent diffractive imaging,CDI)可以同时获得样品的吸收衬度和相位衬度信息,但其实验数据仅记录了衍射强度信息,而丢失了相位信息Φ(x,y).然而,只要样品尺寸足够小且入射光满足相干条件,那么在对衍射图案进行足够精细地采样后,原则上可以恢复出相位信息[1],相位恢复是衍射成像的核心问题.

早在1952 年Sayre[2]就提出CDI 理论最基础的部分,直至1982 年Bates[1,3]提出相位恢复的过采样理论,才使CDI 理论趋于成熟.但是CDI 发展初期由于硬件条件导致图像重构的效果很差,限制了该方法的推广与发展.1999 年,华人科学家缪健伟[4]在美国国家同步辐射光源(NSLS)上首次成功演示了CDI 成像实验,自此CDI 技术开始得到快速发展和广泛应用.1972 年,Gerchberg 和Saxton[5]提出第一个实用的迭代相位恢复算法Gerchberg-Saxton 算法,简称G-S 算法.但是该算法需要已知样品强度图像,在复空间满足的相位Φ(x,y)有无数个G-S 算法会出现重建结果不唯一的情况,限制了该算法能应用的场合.在G-S 算法的基础上,Fienup 于1978 年和1982 年分别提出了混合输入输出(hybrid input and output,HIO)算法[6]和误差递减 (errorreduce,ER)算法[7].然而,这些CDI 算法存在以下问题,仅利用物域约束和衍射强度的约束将会导致相位重建结果不唯一,即会出现孪生像问题[8];算法往往会因为发散或停滞而无法收敛到全局最小值;同时由于光学技术和工程限制,光斑的大小通常是微米量级,为了满足相干条件,通常要求样品尺度小于光斑大小,从而使得成像视场受限.

为了解决上述问题,2004 年Rodenburg[9]和Faulkner[10]提出一种基于全场扫描的无透镜成像方法叠层相干衍射成像 (ptychography) 方法以及叠层相干衍射成像迭代算法 (ptychographic iterative engine,PIE).它是基于经典的CDI 技术[5,11]发展而来,可以看作是扫描透射X 射线成像 (scanning transmission X-ray microscopy) 以及CDI 技术的结合.Ptychography 算法核心思想是要求相邻扫描区域有一部分重叠,重叠区域内物体分布函数的一致性构成了额外的约束条件,从而加速相位恢复算法的收敛,确保得到唯一解.运用迭代投影算法解决相位恢复问题理论上可以达到衍射极限精度[12].2007 年,实验成功验证了PIE 算法[13],2008 年ptychography 理论进一步成熟[14],实现了光斑和样品的同时重建功能.相比较传统的成像方法,ptychography 的优势在于: 1) 基于扫描的方式不受光斑有限视场的约束,能够对大尺度样品进行成像;2) 通过调节聚焦元件和样品平面位置得到纳米至微米量级的光斑,在分辨率范围上可动态调整并且可以接近衍射极限的分辨率;3) 可以同时重构出探针函数和样品分布函数,具有较强的数值反演能力和波前重建的功能;4) 对误差的容忍度较大.因为ptychography 具有这些优势,其在可见光、X 射线和电子束等不同领域的相位重建和波前检测[15]均有应用.同时ptychography 与其他光学技术结合,在生物医学[16,17]、化学[18]、计量学[19]以及超快激光[20]等领域也均有着广泛的应用.

本文首先介绍ptychography 领域中基本算法,然后总结ptychography 基本算法在实验效率提升、厚样品成像和误差校正等方面遇到的关键问题,并介绍优化这些问题过程中ptychography 算法的发展脉络.特别地,还将介绍ptychography与人工智能结合优化带来的新方法和新思路.接着总结ptychography 物理算法模型具体的计算机软件实现,以及海量计算需求带来的并行化方面的发展,最后展望在未来第四代同步辐射光源HEPS 上的发展前景.

2 Ptychography 算法

2.1 传统算法概述

Ptychography 最主流的算法是PIE 算法和DM 算法(DM 算法是基于差分图所设计的算法)[21].2007 年Rodenburg[13]成功完成第一个使用PIE算法的光学实验后,PIE 算法势如破竹地进入了成像领域.如图1 所示,已知入射光束,假设一个待测样品的初始分布函数为O(r),光束分布函数为P(r-Rj),不考虑样品厚度,出射光分布函数可以表示为

图1 (a) Ptychography 的基本光路图;(b) PIE 算法对应的算法流程和重建结果[12];(c) DM 算法对应的算法流程和重建结果[21]Fig.1.(a) Schematic diagram of the basic optical path of ptychography;(b) the corresponding algorithmic flow and reconstruction results of the PIE algorithms[12];(c) the corresponding algorithmic flow and reconstruction results of the DM algorithms[21].

(1)式作为初值参与第一次迭代计算(r,R为空域坐标),迭代过程中首先将出射光变换到频域中:

然后将探测到的强度值Iobs(k,Rj) 替换之前的强度值,反变换回空域:

再更新样品和计算误差:

重复上述流程,直至完成所有位置的迭代.

PIE 算法的核心思想在于使用梯度下降方法(又称最速下降法)来找到全局最优解.对第n次迭代探测器平面的衍射光做频域约束后,二者的偏差可以表示为

误差看作关于不同迭代次数ψ(r) 的函数.对误差B关于ψ n(r) 进行一阶泰勒展开,再考虑理想情况下偏差为0,可得

可以看出(8)式正是梯度下降的形式,详细推导过程参见文献[7].

该算法依然存在一些问题,例如需要知道光斑探针分布,对物体的初始化分布函数敏感,不合适的初始化可能会使算法落入局部最优解或导致收敛速度慢、不稳定等.为了优化这一系列问题,Maiden和Rodenburg[22]根据实验以及数据的特性提出了ePIE 算法.与PIE 算法相比,ePIE 算法具有不同的待测样品分布函数和光斑波函数的重建更新形式,可以同时解出探针和样品分布.PIE 算法不仅在初始化上敏感,在具体的实验操作中,实验平台的移动难免会有抖动,光束分布与实际分布的重叠区域有一定的差距,重叠区域的多少会影响PIE算法的收敛性[23].2.2.3 节会详细介绍实验误差带来的影响以及在算法上如何优化来减小误差.

DM 算法在2009 年由Thibault 等[21]首次提出.DM 算法是一种非线性优化迭代算法,最早用于解析晶体学中的相位问题[24],如图1(c)所示.它有效地降低了图片伪影,提高了图像重建质量,并且可以同时重建探针和样品分布函数[10,21].虽然DM 算法不拘束于扫描顺序来计算样品和探针分布函数,但是需要在每一次迭代过程中,计算所有扫描位置的O(r) 和P(r-R) 来进行函数的更新,因此单次迭代的计算量较大.

2.2 传统算法优化

2.1节介绍的是ptychography 基本算法,但是在实验过程中依然存在许多问题和挑战.例如,周期性网格式扫描容易导致伪影,通过引入非周期或非对称的扫描路线[25](例如螺旋式,同心圆和随机扫描等) 可降低伪影的影响.此外,针对ptychography 的实验特点,研究人员发展了许多不同的实验方法来提高图像重建速度和成像分辨率等.针对传统ptychography 扫描完一个区域并且记录数据后再进行下一位置的扫描,导致实验效率低的问题,2.2.1 节介绍fly-PIE 算法和单次曝光ptychography算法(single-shot ptychography,SSP).针对传统算法仅适用于薄样品而实际样品厚度在实验中往往不能忽视的问题,2.2.2 节介绍厚样品成像算法3PIE.针对ptychography 实验结果对扫描位置的准确性和探针的相干模态敏感、实验误差等问题.2.2.3 节介绍面向位置和相干模态的优化算法.

2.2.1 效率优化算法

2.2.1.1 fly-PIE 算法

Ptychography 实验过程中需要对每一个离散扫描区域进行步进扫描,每扫描一块区域电机将会反复经历启动-停止过程,带来无意义的时间消耗.2014 年Clark 等[26]理论上证明了连续扫描和步进扫描采集到的图像重建后的效果相当,为提高扫描效率提供了理论基础,进而提出fly-PIE 算法.其算法思想是将探测器采集时间设为T,选择足够小的时间步长 Δt,将探测器采集时间分为n个区域,以恒定速度v扫过的小区域vΔt作为一个数据模态,不同数据模态求和后的值为探测器所采集的衍射强度值.第j个位置的衍射强度等于扫描时间T内强度的叠加:

其中,扫描时间T=t2-t1恒定不变,P z表示旁轴传播,P(r) 表 示探针函数,O(r+rj+vt) 表示第j个位置飞扫后的待测样品函数.如果选择一个小的时间区域则可以对(9)式做离散化处理:

每个模态代表一种完全相干的衍射模式,但不同模态之间非相干,其强度相加是一种非相干叠加,可较好地描述光或样品的部分相干状态.在部分相干叠加下的相干状态与完全相干叠加(传统ptychography)的效果持平,前者甚至更好.

2.2.1.2 SSP 算法

fly-PIE 算法虽然通过飞扫替换步进扫描缩短了电机移动的时间,提高了实验效率,但是依然存在以下问题.一方面,由于每个扫描位置需要足够的曝光时间获取衍射信号,电机移动的物理速度是有上限的,依然存在电机移动时间浪费.另一方面,电机快速移动中的抖动会产生光斑位置偏差,电机移动越快,出现光斑偏差的概率越大,进而导致重建的探针分布函数出现错误.2013 年中国科学院上海光学精密机械研究所潘兴臣等[27]提出了SSP算法.该算法避免了通过电机移动扫描样品,由于无需电机移动样品,也无需考虑光斑位置误差,可以一次性扫描大范围的样品面积使得实验效率大大提高.实验示意图如图2 所示,算法思想为一束相干光经过准直照射到网格光栅后发生衍射,网格光栅会使入射光发生分离.根据光栅的规格不同,照射在样品表面的光探针数量不同.通过设置光栅与样品的距离使相邻的光探针在样品表面有一定的重叠.在不放置测试样品的情况下,可以测得入射光通过网格光栅后的图样.通过相位恢复算法或传统的PIE 可以得到精确的入射光分布函数.根据光栅的间距以及物函数G(x,y) 可以给探测器端所记录的各个衍射图案进行标号 (m,n) .不同标号的光斑表示为

图2 (a) fly-PIE 的原理示意图[26];(b) SSP 算法原理示意图[27];(c) fly-PIE 算法单模态和多模态的重建图像对比;(d) SSP 实验对强散射样品和蜜蜂翅膀的光斑、模量和相位结果Fig.2.(a) Diagram of fly-PIE[26];(b) diagram of SSP[27];(c) comparison of reconstructed images of fly-PIE algorithm using single and multi-modal;(d) probe,modulus,and phase results of the SSP experiment for stronger scattered samples and the bee wings.

(11)式对应于照射至样品表面的光斑探针,同时可作为图像重建的入射光探针函数.

2016 年Sidorenko 和Cohen[28]进一步改进了SSP 算法.他们基于4f成像系统(由两块透镜同轴共焦组成,物平面与像平面距离大于等于4 倍焦距)将网格光栅换为微孔阵列来进行分光,有效地提升了成像效果.由于物体同时被多个重叠探针照射,CCD 端的强度表示为

其中,角标m是微孔阵列标号,k m为相应微孔分光的波矢.由于(12)式求和号可以放在绝对值外,因此可以对第m个探针的单个模式应用传统的PIE 算法或是 ePIE 算法进行重建.但是因为总的强度是求和的,所有光束的多个衍射图案实际上会相互干扰.因此,当使用区块中的衍射图案进行重建时,来自其他针孔的光束对衍射图案的贡献被视为噪音.为此,Sidorenko 和Cohen[28]计算了微孔数造成的噪声极限,以此更改焦距、物距等光学系统参数,减小噪声的影响.同时他们还改进了 SSP实验不同孔径构型,可能会提高算法的收敛速度.对于噪声减小的优化方案,尝试不同的微孔形状或许会减少噪声的影响从而提高鲁棒性.2018 年,Chen 等[29]提出了多路复用的SSP 实验技术,不仅能重建探测光束的振幅和相位,还能恢复其偏振态.目前,SSP 算法已经运用到了多个领域,例如光学编码[30]和加密系统[31,32]等.

虽然fly-PIE 和SSP 两种算法分别从时间和空间上提高了算法效率,但这两种实验算法本身仍然存在缺陷.fly-PIE 的问题在于如何确定部分相干模态的数量,而且使用部分相干模态相加会导致强度降低.SSP 技术的局限性在于: 1) 分束的衍射信号均由一个探测器同时采集,衍射信号的空间频率范围将非常有限,不同的高频衍射信号会有一定弱干涉,无法获得高分辨率图像信息;2) 实际应用中分光光栅的制作成本相对较高.改进后的SSP虽然使用了不同的光栅构型,但是噪声问题变得更加严重.如何降低SSP 噪声以及不同光栅构型对噪声的影响是后续发展的优化方向.

提高成像效率的算法还有傅里叶叠层显微镜成像(Fourier ptychographic microscopy,FPM).FPM 是一种合成孔径技术,能够在非相干光实验条件下重建图像.由于篇幅限制不在此详述,读者可阅读相关文献[33-36].

2.2.2 厚样品成像算法

Ptychography 实验技术一般适用于薄样品衍射成像,厚样品容易出现探测器信号缺失,重建图像不准确等问题.但这样就限制了ptychography的应用范围.Maiden 等[37]提出了一种新的算法,可以重建厚样品,称为3PIE 算法.他们在图像重建中将样品切片,对每一片进行重建,然后再按照传播顺序一层一层地恢复,如图3 所示.下面给出3PIE 算法的具体流程.

图3 (a) 3PIE 算法的示意图[37];(b) 使用 3D-ptychography 技术对部分纸张组织的图像采集结果[37] (图(b)前两幅是上层组织的模量和相位,后两幅是下层组织的模量和相位)Fig.3.(a) Schematic diagram of the 3PIE algorithm[37];(b) the result of image acquisition of some paper tissues using the 3Dptychography technique[37] (The first two pictures shows the modulus and phase of the upper tissue,and the last two shows the modulus and phase of the lower tissue in Fig.(b)).

1) 确定所需重建的图像数据,标记其强度为Ic(u),初始化探针函数和样品中每层切片的分布函数.

2)计算第一层的出射波前

其中,Rc是扫描路径.

3)将第一层出射波前传播至下一层,作为下一层的入射波前

4)计算第二层的出射波前

同样再将其传播至第三层作为入射波前

5)如此计算直至最后一层出射波前ψe,N(r) .

6)使用傅里叶变换或菲涅耳衍射传播至探测器,得到频域空间的理论值

7)用实际探测得到的强度值替代理论计算得到的强度,

8)再将其进行反向传播,得到第N层切片的修正估计值:

设置更新函数:

9)更新第N层切片入射波函数和样品分布函数,

10) 再反向传播至上一层切片,得到上一层修正的出射波:

11) 重复步骤9)和步骤10),计算每一层切片,直到第一层切片的波前.

12) 更新波前探针函数和样品分布函数,

13) 将每层的样品分布和探针更新:

并将其作为下一次迭代的初值.

14) 重复上述步骤1)—13),直到满足终止条件或达到固定迭代次数.

此外,3PIE 算法要求在扫描期间待测样品不能发生变化,入射光探针的状态稳定,收集的强度数据集输入模型时与相应切片传播的测量衍射图案相匹配.3PIE 算法继承了PIE 算法的易理解性.但是由于对样品进行切片操作,同时还有扫描曝光操作,无疑会对计算量和计算速度提出挑战.切片数量的多少会直接影响重建的效果以及计算效率.但是在噪声方面,切片操作会让噪声的权重降低从而提高图片重建的质量.

2020 年,Barutcu 等[38]提出了一种与梯度下降算法结合的Ptycho-Tomography 算法,避免了厚样品重建中不同层之间的伪影问题.在材料领域[39-41],3D ptychography 算法可以重建出更多的细节而被广泛应用.3D-ptychography 算法还与SSP 算法相结合[42],在近轴向对厚样品照明具有明显的抗干扰性和良好的图像重建质量.

2.2.3 误差校正算法

通过前文对PIE 算法的简要概述可以看出迭代算法通常需要探针位置足够精确,否则需要能够设计出可以自动修正探针位置的算法.上述算法需要光斑及相关参数严格满足两点: 一是扫描过程中位置足够精确;二是函数参数设置足够确定或位置参数个数不能超过一定数量[43].电机的移动、辐射的热效应、光源的变化、地面的微振动等因素影响会对高分辨纳米成像实验产生严重影响.例如,光斑抖动会造成样品曝光时间内,光斑呈现模态不唯一的情况(光源的辐射场发生改变,相应的相干模式也会发生变化);辐射的热效应也会使样品内部发生改变,不再是静态混合;此外,同步辐射装置中的振动、镜子的热变形或是缺陷、样品的偏移等都会使理论扫描位置和实际扫描位置出现误差.下文主要介绍多模态条件下优化位置误差和修正参数不确定性的算法.

2.2.3.1 多模态优化算法

多模态可以理解为一束光中不相同的相位、模量信息的叠加,导致光是部分相干.早在1998 年Nugent 研究小组[44]明确定义了部分相干光场.他们用平均坡印廷矢量定义光场,并且结合连续性方程推导出部分相干光的相位表达和恢复方程.而后在2009 年Nugent 研究小组[45,46]基于迭代算法将部分相干光进行模式分解,从而达到信息的均匀与恢复.如图4(a)所示,多模态产生原因主要有三点[47]:1) 光源本身的辐射,内部有一定的退相干以及光源的振动导致的部分相干;2) 样品内部发生量子态的混合和平稳随机过程;3) 探测器的点扩散特性所导致的衍射信号模糊,形成另一种退相干机制.

图4 (a) 多模态产生原因的示意图和多模态重建图[47] (其中模糊的图片为单模态重建，下方为探针的 5 个主要模态);(b) 单层WS2 样品在不同电子剂量条件下单模态和多模态的重建图像[49];(c),(d) 分别使用OPRP 算法和ePIE 算法的重建结果和探针的4 个主要模态占比[48]Fig.4.(a) Schematic diagram of the causes of multimodal generation and the multimodal reconstruction[47] (The blurred image is the unimodal reconstruction,and the 5 main modes of the probe are shown below);(b) the reconstructed images of unimodal and multimodal reconstructions of the monolayer WS2 sample under different electron dose conditions[49];(c),(d) the reconstruction results using OPRP algorithm and ePIE algorithm respectively and the four major mode occupancies of the probe[48].

因为在实验中存在扫描过程,扫描过程的各种振动使得光线的模态发生改变(位置和波前等因素)从而导致退相干效应.并且因为需要移动平台且每一步的移动可能会有些许误差,导致重叠区域发生变化.而且大型同步辐射装置使用的 X 射线辐射有些许不同,在考虑辐射效应下,有些样品的结构或许会发生改变.以上所有因素叠加会导致整个实验出现多模态情况.如果将每一次扫描位置中的模态视为一个系统整体,不同模态之间是非相干叠加.这样最后的图像重建便是所有不同模态重建出的图像的非相干叠加.

假设有m=1,2,···,M个模态产生,不同模态集合可以记作

与 PIE 算法基本相同,其算法如下.

1) 计算第j个扫描位置M个模态的传播函数

其中下标k代表迭代次数.

2)傅里叶变换后传播至探测器,计算误差函数

3) 对不同模态传播函数的振幅进行加权修正:

4) 进行傅里叶反演至空域.

5) 修正物体和探针:

6) 顺次改变扫描位置,重复上述步骤1)—5),完成第k次迭代.得到M个探针和样品函数,并计算每个像素的平均误差:

7)其中,X和Y代表探测面积.

8)上述计算值作为下一次迭代的初值,迭代计算直至误差最小.

多模态的产生是不可避免的,但总体可将其分解成多个部分相干光的叠加.2016 年Odstrcil 等[48]提出了OPRP 算法(orthogonal probe relaxation ptychography).他们利用奇异值分解,将探针分解成多个模式,得到了不错的图像重建结果.Chen[49]在2020 年使用多模态探针重建了纳米结构材料的原子结构,并且与单一模态实验结果对比发现在追求高分辨率、高精度和大视场成像时多模态重建的重要性.同时多模态相比传统的算法会产生更大的数据量.为了提高图像重建的效率,需要采用GPU 并行加速.PyNX[50],Ptypy[51]和Ptychopy[52]等软件均支持GPU 并行处理图像数据,提高重建效率.

2.2.3.2 位置精度误差修正算法

上文提到ptychography 实验中光斑探针位置会对图像重建有明显的影响.Rodenburg 等[43]指出探针倾斜、漂移和位置随机误差都会对图像重建质量有影响.在TEM-ptychography 上进行了一系列的模拟实验,结果如图5 所示.探针位置的变化不仅仅会使图像重建的质量下降,而且如果数据集整体发生位置偏移,会使得其自身是自洽的,但是重建出的图像整体会有一定的错误,从而误导我们对图像的分析.因此非常有必要修正光斑位置误差.未来第四代同步辐射光源HEPS(high energy photon source)上,高相干性光源会降低多模态的影响,但由于各种振动,工程设计和光学元件产生的位置误差是无法避免的.如何设计出更高效,更准确的位置误差算法是未来提高图像重建质量的重要方向之一.

图5 位置误差对图像重建质量影响示意图[43] (a) 对应图(b)—(e)的探针位置;(b) 随机位置误差重建的结果;(c) 探针逆时针倾斜重建的结果;(d) 探针漂移重建的结果;(e) 探针位置在垂直方向上有扩展的重建结果Fig.5.Diagram of the effect of position error on image reconstruction quality[43]: (a) Probe positions in Fig.(b)—(e);(b) result of the random position error reconstruction;(c) result of counterclockwise tilt of the probe reconstruction;(d) result of the probe drift reconstruction;(e) result of the reconstruction of the probe position with extension in the vertical direction.

退火算法[53]的具体流程和ePIE 算法基本类似,除了在样品和探针迭代修正方程中添加了对位置误差的修正.退火算法主要目的是修正探针位置的平移、偏移、倾斜、漂移以及随机位置误差的产生.退火算法的示意图如图6 所示,蓝色小圆圈代表应测位置,蓝色矩形代表实测位置,虚线圆圈代表产生随机位置误差.

首先引入修正向量C j,以及没有随机位置误差的传播函数

然后设置m=1,2,···,M个随机位置误差r0(k)Δm,其中r0(k) 的取值随迭代的次数线性减小.Δm=[δx,m,δy,m],随机位置误差向量的取值在±1 之间.含有随机位置误差的传播函数为

与ePIE 不同的是,需要选择出与真实情况最为接近的传播函数作为算法的迭代函数,做法如下.

1)将ψ0和ψ m做傅里叶变换Ψ m(u)=Fψm得到探测器上的被估计值Ij(u) (u为倒易空间坐标).

2)计算实际探测与估计值的误差:

3)选取误差结果最小的传播函数作为迭代算法中的反向传播波函数,记作Ψ n(u) .

4)将探测的振幅值替代理论估计值,得到反向传播后的出射波函数:

5)更新位置修正向量

6)更新样品分布函数:

7)更新探针函数

重复上述步骤直到所有扫描位置均完成迭代.

之后,Maiden 等[53]提出可以修正全局位置误差的算法,称为e-pcPIE.其算法思想是额外考虑以下误差: 探测器平面相对样品平面的旋转、样品和探测器距离测量的不准确造成缩放比例的误差,以及数据获取过程中的线性漂移误差.与上述方法一致,只是将修正向量加入关于旋转和距离的误差向量和标量,每次迭代更新时用随机小量更新迭代,重建后的图像质量显著提升.

2.2.3.3 相关匹配算法

相关匹配算法(serial cross-correlation)[54]在交叉相关算法[55](cross-correlation)基础上演化而来.由于前后两个扫描位置的重叠具有一定的关联性,相关匹配算法的核心思想是运用互相关函数求极值,得到两次迭代过程中样品分布函数修正的位置误差,从而对每一次迭代中的位置进行修正，如图7 所示.与ePIE 算法的具体迭代过程类似,相关匹配算法在更新样品分布函数后对其位置进行更新修正.样品分布函数宗量为实域r和位置误差修正R j,m,角标j代表扫描位置,m代表迭代次数.根据相关函数

图7 (a),(b) 在每一个扫描位置x 和y 方向的位置误差[54];(c),(d) 使用可见光的X 射线重建样品图像结果[54];(e)—(h) 使用调制器的X 射线重建样品图像结果[54];(i)，(j) 使用KB 镜的X 射线重建样品图像结果[54].其中图(c),(e),(g),(i)是使用相关匹配算法优化的结果，图(d),(f),(h),(j)是未使用相关匹配算法优化的结果Fig.7.(a),(b) Position errors in the x and y directions at each scan position[54];(c),(d) X-ray reconstructed sample image results using visible light[54];(e)—(h) X-ray reconstructed sample image results using modulator[54];(i),(j) X-ray reconstructed sample image using KB mirror[54].Fig.(c),(e),(g),(i) and Fig.(d),(f),(h),(j) are the results optimized with and without using the corss-correlation algorithm,respectively.

对其求极值可以得到位置误差修正值C j,m.将下一次迭代的位置误差修正为

参数κ是一个常数或是随迭代次数减小的线性量.函数Λ m(r) 代表光斑函数,选择对样品探测贡献最大的区域.将修正后的宗量带入样品或探针函数后进行重构即可.

位置误差对相位项的影响较大,相关匹配算法解决的是平移误差所带来的相位影响.相关匹配算法是运用图像处理领域中自相关系数的思想实现的一种算法.将前后两次的迭代过程整体视为一张具有相互关联性的图像,从而找到其中最大关联,进而修正位置误差.

2.3 Ptychography 与人工智能

近年来人工智能在图像处理领域取得了很多不错的成绩,例如图像识别、阈值分割、目标检测、去噪和超分辨恢复等.国际上各个国家也在大力发展人工智能领域在基础科学等实际应用方面的价值.美国在过去5 年资助资金达14.24 亿美元多[56].图像处理与光学有密不可分的关系,科学家们希望利用人工智能在图像处理上的准确性、鲁棒性、泛化性以及提升效率等优势,对光学领域的发展做出推动.并且人工智能的深度学习领域适合解决多种任务模式相混合、任务复杂等问题.复杂问题的解往往是在高维空间，而传统的计算方式几乎不可能完成.由此催生了AI for Science 这一科学范式.人工智能可以赋予科学家更高效、更强大、更准确的计算推理能力,将人工智能作为工具在科学领域使用将成为主流.Ptychography 作为光学领域分支,不仅有物理建模过程还有后续的图像处理需求.尽管传统ptychography 及其优化算法在不同应用场景下均取得了不错的图像重建效果,但是对于多种实验误差相互耦合的情形,传统算法将束手无策.机器学习、深度学习等人工智能算法,可一次性完成多种任务的处理,适合处理复杂任务相互耦合的情形.基于人工智能的ptychography 图像重建方法逐渐成为研究热点.

2017 年Maiden 等[57]基于机器学习中的动量梯度下降算法思想提出了mPIE,经过一定次数迭代后在待测物体的分布函数更新公式中添加动量项使得算法的迭代次数大量降低,收敛速度显著加快.同年,Kappeler 等[58]首次提出基于卷积神经网络 (CNN)对傅里叶叠层成像 (Fourier ptychography,FP)的图像进行重建.Kappeler 搭建了PtychNet,其架构如图8 所示,包含3 个卷积层,其中第1 层是64 核9×9 的卷积层,第2 层是32 核5×5 的卷积层,最后通过ReLU(linear rectification function,ReLU) 激活函数全连接.在扫描位置无重叠和有重叠两种情形下,与迭代误差下降（iterative error reduction algorithm,IERA）算法[59]重建进行对比,PtychNet 在没有重叠的条件下与IERA 算法持平,但是重叠率大于50%的情况下图像重建的信噪比略低.不过PtychNet 的计算速度却快于IERA 算法,并且在低迭代次数的情况下,PtychNet 重建的图像的信噪比高于IERA 算法.而后FP 与机器学习方法结合,也提出了其他不同的网络构造[60,61],同样有很好的图像重建效果.Ptychography 领域的工作逐渐发展起来.

图8 PtychNet 网络架构示意图和重建图像结果[58] (a),(b) 扫描无重叠的图像重建;(c),(d) 扫描有重叠的图像重建Fig.8.Schematic of the PtychNet network architecture and reconstructed image results[58]: (a),(b) The image reconstruction with no overlap of scans;(c),(d) the image reconstruction with overlapping scans.

探测器的本底噪声和光源所带来的泊松噪声在实验中一般不可忽略.噪声导致迭代算法重构出的图像质量会变差.Metzler 等[62]在2018 年结合正则化去噪（regularization by denoising）[63]以及DnCNN 网络[64]构造出prDeep 网络.因为其构造形式使得prDeep 可以应用于相位恢复和去噪等各种问题领域.图像伪影问题在重建算法里也很常见,2019 年Işıl 等[65]通过结合深度神经网络（DNN）和 HIO 算法构建了新的相位恢复网络.他们将DNN 网络嵌入HIO 的迭代过程中，每一次迭代过程中用DNN 网络去除图像伪影,DNN 网络中的卷积自带部分去噪和分辨率恢复的能力,最后再用一次DNN 网络进一步恢复高分辨图像,效果显著.

Cherukara 等[66]构造了网络PtychoNN,可以基于光源数据同时重建出振幅和相位信息.如图9所示,与传统的ePIE 算法进行对比,发现稀疏采样条件下,PtychoNN 表现出比ePIE 更好的图像重建效果.这意味着重建图像时不一定需要完全满足重叠面积50%以上,可以适当放宽条件.重叠条件的放宽意味着光源振动、探测器振动等因素导致的位置误差和多模态也可以忽略.他们还指出稀疏采样使得采样点密度减少,使得特定分辨率成像的辐射剂量降低.

图9 (a)—(d) PtychoNN 与ePIE 重建的结果对比[66] (a),(c) 使用ePIE 算法重建的结果;(b),(d) 使用PtychoNN 重建的结果.(e),(f) PtychoNN 在少量数据集上的重建质量的对比,发现在至少800 个数据集上可以得到合理结果;(g) PtychoNN 网络的结构Fig.9.The left figure above shows and the right figure shows the comparison of the results of the reconstruction of PtychoNN and ePIE[66]: (a),(c) The results of the reconstruction using the ePIE algorithm;(b),(d) the results of the reconstruction using PtychoNN.(e),(f) Comparison the reconstruction quality of PtychoNN on a small number of datasets,and find that reasonable results can be obtained on at least 800 datasets;(g) structure of the PtychoNN network.

语音信号处理领域同样存在相位恢复问题.对ptychography 成像技术与语音信号处理中短时傅里叶变换(short-time Fourier transform,STFT)的相似性进行分析后,Welker 等[67]在2022 年提出了深度迭代投影(deep iterative projections,DIP)神经网络算法(图10).对音频信号做STFT 后，使用幅度谱还原信号但丢失相位信息.他们在语音相位检索DeGLI (deep Griffin-Lim iteration)架构的基础上优化ptychography 算法,将交替投影(alternative projection,AP)算法的输出结果作为DNN 网络的输入,DNN 旁支对其进行结果的预测并且与AP 算法的结果相结合后作为下一轮迭代的输入.他们还通过模拟比较了将DIP 网络作为初始迭代步骤,然后转用DM 算法进行迭代和只使用DM 算法两种实验结果,发现DIP 网络作为初始重建步骤在低迭代轮次下的收敛速度优于仅使用DM 算法.同时他们还比较了相同的误差阈值下,不同算法达到阈值的迭代次数和时间.结果也表明DIP 算法要优于传统的相位恢复算法.神经网络的优势在于更快的图像重建速度和低迭代次数图像的精度比传统算法收敛结果更好.

图10 (a) DIP 网络的结构示意图;(b) 使用DIP 网络、DM 算法与AP 算法重建图像结果的对比[67],其中DIPcg 和DIPvm 指两种不同模式的DIP 网络.右上方表格为误差阈值为0.1 的情况下各个算法需要的迭代次数和时间,可以看到DIP 网络表现很好.(c) DIP 网络运用到MINST 库上的测试结果Fig.10.(a) Schematic diagram of the structure of DIP network;(b) the comparison of the reconstructed image results using DIP network,DM algorithm and AP algorithm in reconstructed image results[67],where DIPcg and DIPvm refer to two different modes of DIP networks.The table on the top right shows the number of iterations and time required by each algorithm for an error threshold of 0.1,and we can see that the DIP network performs well.(c) Test results of the DIP network applied to the MINST library.

此外,其他传统ptychography 算法也在结合人工智能算法对相位恢复算法进行优化.例如针对SSP 的神经网络算法[68],具有浴幕效应的神经网络算法[69],在无标注数据情况下自主学习相位问题的神经网络 AutoPhaseNN[70]等.

神经网络拓展了CDI 的方法学,但是重建策略上仍然存在部分问题.首先,对于监督学习算法,数据集的制作至关重要.现阶段大部分监督学习任务大都使用模拟的ptychography 数据,而未提及真实实验上的重建.模拟数据几乎未考虑位置误差、多模态、光源的部分相干,与实际情况出入较大.最后重建效果难以保证.对于非监督学习算法,自动编码和解码器虽然有一定优势.但所需要的数据量大,训练时间长,并且还需考虑具体实验中的噪声影响.针对不同线站所收集的数据形式,需要精心设计相应的数据处理方式来适应神经网络的训练.

3 Ptychography 算法并行实例

随着 ptychography 广泛的应用,ptychography实验过程的自动化、智能化计算需求增多.一方面迫切需要集成数据采集、数据处理、图像重建等的软件体系.另一方面随着同步辐射光源性能的提升,将会采集到海量的科学大数据,迫切需要能够应对海量数据传输和处理的高性能软硬件体系.为了方便读者能够及时根据自身的硬件平台情况测试与应用相关ptychography 软件,本文按照软件依赖的硬件平台划分为GPU 并行软件和CPU 并行软件,分别在3.1 节和3.2 节介绍.

3.1 GPU 算法并行实例

2014 年,Nashed 等[71]实现了一种基于GPU并行计算的ptychography 软件平台(图11).为了提高图像重建计算效率,将图像分割后用多个GPU 设备进行并行重建,提出了异步并行法和同步并行法.异步并行法通过MPI 对多个GPU 进行管理,将扫描区域分成多个子区域后分配给不同的GPU 进行图像重建,然后根据分配位置水平和垂直方向的梯度关系进行图像拼接.每个GPU 重建的图像信息是独立的,各个 GPU 之间不需要共享数据.与异步并行法不同,同步并行法是各个GPU 在处理不同扫描位置的数据时,相邻区域重叠处的数据在不同GPU 之间会有共享,确保重建图像的准确性.同步并行的数据共享类似于实验中的区域重叠,数据的共享也会加快算法的收敛,同时也符合直观的实验过程.各类ptychography 图像重建软件的GPU并行加速大多按照同步并行的方式.重建区域分割为多个不同子区域会导致重建未重叠区域边界信息时可能会产生伪影,上述两方法都有一定程度的克服.

图11 (a),(b) 分别为GPU 异步并行和同步并行的数据处理示意图[71];(c),(d) 分别为Ptychopy 软件的工作流程图和GUI 界面[52]Fig.11.(a),(b) Schematic diagrams of data processing with GPU asynchronous parallelism and synchronous parallelism,respectively[71];(c),(d) the workflow diagram and GUI interface of the Ptychopy software,respectively[52].

2016 年Mandula 团队[50]搭建了计算散射、ptychography 的PyNX 软件工作平台.最开始利用DM 算法和ePIE 算法进行图像的重建,并且讨论了高斯噪声条件下,实验数据的最小化和所有像素强度的平方差的等价性,再利用理论计算和观测数据的结果计算泊松似然估计得到最好的重建结果.但是PyNX 只有计算散射的模组可以进行GPU 并行加速.软件在2016 年成功运用欧洲同步辐射光源(European synchrotron radiation facility,ESRF)的数据重建出图像,并且运用最大似然法来优化更新探针和样品分布函数.2020 年,PyNX软件进一步拓展[72],可计算波前传播、分层成像以及所有的成像操作等.并且对ptychography 模块进行了更新重写,增加了计算CDI 和波前的模块.采用GPU 同步并行加速框架,部分底层采用C++语言使整体运算效率提高.同时,PyNX 支持多种算法混合使用,使重建过程不再拘泥于单一的算法,可以应用不同算法组成重建链来探索不同算法组合的可能.PyNX 集成了多种功能,并且高度集成的工具包可以使新算法的开发更加便捷,但没有GUI 界面会使其入门门槛较高.本文搭建了PyNX 运行环境,重建了不同大小的数据,比较了DM 算法和ML(max-likehood)算法[73]的重建时间和单次重建时间,如图12 所示.

图12 (a) 左图为使用Ptychopy 在同一数据集上测试不同的迭代次数的软件速度,右图为使用PyNX 测试模拟数据,分别处理200,500,1000,3000 张探针大小为256×256 的数据集.横坐标代表数据的大小,纵坐标代表算法实现的时间 (每一步DM 算法会分批次计算LLK(log-likehood)).(b) SHARP 的GPU 并行示意图和重建流程图[74].(c) 左图为PyNX 重建样品的部分相位图,右图为分配给不同 GPU 的扫描位置,部分位置的信息将共享[72]Fig.12.(a) The left figure in (a) shows the software speedup using Ptychopy to test different number of iterations on the same dataset,and the right figure shows the simulated data using PyNX to test and process 200,500,1000,and 3000 datasets with probe size of 256×256,respectively.The horizontal coordinate represents the size of the data,and the vertical coordinate represents the implementation time of the algorithm (LLK (log-likehood) is computed in batches for each DM algorithm step).(b) GPU parallelism of SHARP and the reconstruction flow chart[74].(c) The left figure in (c) shows the partial phase map of the PyNX reconstructed sample,and the right figure shows the scan positions assigned to different GPUs,and the information of some positions will be shared[72].

同年,Marchesini 团队[74]研制了一个高性能软件—SHARP,并应用在高级光源(advanced light source,ALS)中.Marchesini 等[75]设计了一套高性能、高效率的算法,通过对算子的操作实现RAAR 算法(relaxed averaged alternating reflections algorithm)来重建实验过程和实验数据结果,并且实现了高通量的流分析.GPU 并行计算也遵从同步并行计算的原则,将扫描区域划分子区域再由不同的GPU 重建.还设计了即时快速反馈系统,可以将实验结果快速反馈给用户,这种软件结构允许用户对他们的实验进行直观、灵活和反应灵敏的监测和控制.数据采集和分析之间需要紧密地整合,以使用户从中获得他们所期望的反馈.

2021 年Yue 等[52]针对计算量更大的数据集,实现了Ptychopy 计算软件.该软件采用混合式并行计算实现相位恢复.后端基于CPython,C++和CUDA,用户交互界面和系统控制基于Python库.它集成了GPU 图像重建模块和光束线软件模块(如光束线控制、数据收集和存储模块),向用户提供了用户交互式界面和脚本两种运行模式.Ptychopy 创建了两种工作模式,第一种是创建与GPU 数量相等的工作流.控制器会为指定的GPU创建相应的工作流,将图片重建任务分配给每一个GPU.GPU 采用先入先出的顺序处理所分配的作业.这种一对一的GPU 分配适用于作业数多但单个作业数据量小的场景.第二种是创建只有一条工作流,使用所有的GPU 并行计算.适用于单个作业的重建数据超过单个GPU 内存的场景.

3.2 CPU 算法并行实例

Ptypy 软件[51]可用于比较不同参数设置后的重建结果,作为今后的标准模型框架来降低用户的入门门槛.Ptypy 通过消息传递接口(message passing interface,MPI)将任务分配至每个CPU 节点.同时它还分配了内存数据缓存区,当CPU 核需要处理衍射数据时,缓冲区的数据会并行传输至CPU端,如图13 所示.算法设计上讨论了同一次扫描中,探针的形态或是漂移以及样品的旋转都会使对应的函数发生变化导致结果不准确.但是在不同的扫描过程中会包含正确的结果,所以不同扫描数据之间可以实现数据共享来提高重建准确性.表1 汇总了当前主流软件所使用的GPU 型号、使用的算法、单次迭代计算时间以及估算的A100 GPU 计算时间.

图13 Ptypy 软件的API 控制与实验对应示意图[51],其中Pod 是封装好得到可调用包Fig.13.Diagram of the API control of the Ptypy corresponding to the experiment step[51],where Pod is wrapped to get the callable package.

表1 部分软件计算的硬件信息、使用软件的算法对单次迭代计算时间的估计值[76],以及使用A100 GPU 的性能速度预测值Table 1.Hardware information for part of the software as well as an estimate of the time for a single iterative computation using the various algorithm[76],and the performance is predicted using the A100 GPU.

3.3 其他算法实例

2020 年瑞士保罗·谢尔研究所(Paul Scherrer Institute,PSI)的Wakonig 团队[76]基于Matlab 研制了PtychoShelves 软件.PtychoShelves 的模块化框架并不依赖于面向对象的实现.新的引擎可以很容易地作为一个单一的函数被添加到引擎目录中,可以在原有基础上做任意修改.PtychoShelves使用MPI,OpenMP,Matlab MG 引擎实现了CPU和GPU 的并行加速[77],同时也实现了DM 算法、ePIE 算法和LSQML (least-squares maximumlikelihood)算法[78]的集成和使用.

2021 年Dieter 等[79]开发了Ptychography4.0软件.该软件致力于将叠层扫描显微成像技术借鉴到其他显微成像领域,同时建立一套广泛的显微成像模式标准.Ptychography4.0 软件包含使用单边带方法[80](single side band method)重建扫描电子显微镜实验图像的振幅和相位等实验数据以及传统的ptychography 实验图像.Ptychography4.0软件加入了探测器位置与探针位置的对应关系模块,可以修正由于位置未对齐而导致的图像模糊.并且加入了滤波后的移频操作,使傅里叶变换后的图像高频信号部分移至中心,效果明显.相较于其他软件,Ptychography4.0 并没有集成类似PIE 和DM 算法的模块,而是利用Python 软件包编写了一系列与实验流程相互对应的迭代步骤函数,利用这些函数来实现迭代算法.最后,表2 总结了常用ptychography 的软件特点及标准化的重建时间[76].

表2 不同软件的特点以及标准化重建时间[76]Table 2.Characteristics of different software and standardized reconstruction time[76].

4 结论与展望

随着同步辐射技术的发展,ptychography 面临的大部分问题逐渐得到解决,但依然存在许多问题.在算法上,高斯和泊松等多种噪声对迭代算法的影响,实域的物理约束如何在算法中体现等问题都有待优化.在工程上,ptychography 实验对探针位置精度较为敏感,如何进一步提高实验系统的稳定性还需进一步优化.

当前主流的ptychography 算法大多都是基于PIE 算法或DM 算法.PIE 算法的易理解性以及与实验投影过程的一致性,成为大部分相位恢复的主流选择.但在不同领域PIE 算法面临的问题也不同,例如在可见光领域信噪比的提升,电子束成像中的辐射损伤和滞后效应[81],X 射线领域光源的部分相干性和位置误差的影响等.DM 算法来自于晶体学中的差分法理论.但是对于DM 算法的优化较少.Ptychography 作为高分辨成像主流技术之一,在实验技术上也有一定的发展,例如新的实验仪器[82].在算法上,有EE-PIE 算法[83],ptychography 结合断层扫描的新算法[84]应用在生物医疗领域.Ptychography 实验技术现在已经广泛运用于材料[40,85],生物[86-88]以及量子态的观察[89]等领域,具有广泛的应用前景.

除此之外,ptychography 与人工智能相结合可以使得图像重建、相位恢复算法实现更高的效率和更佳的质量,在高能同步辐射光源上得到广泛的应用成为必然的趋势.然而,一方面随着四代同步辐射光源技术的发展,光源相干性的提高,纳米探针和相干散射等线站采集的数据量将会达到新的高度.且未来束线站会结合多种实验方法,如Ptychographic Tomography[90,91],Resonant Ptychography[92,93],InSitu Ptychography[94]等.Ptychography 与人工智能两者都属于经典的计算密集型任务,对计算能力的要求非常高,现有的计算方法面临严峻挑战.在未来,光源需要研究更高性能的数据采集和图像重构软硬件技术体系,满足ptychography 算法的海量计算需求.另一方面,随着高通量、多模态、原位探测器装置的升级,每秒将会采集海量的科学大数据,高通量大数据的涌入使得传统方法和机器学习方法都面临更加严重的挑战.两者都无法在短时间内进行图像处理.当前围绕高能同步辐射光源的数据处理模式主要采用云计算的架构模式,线站将采集的数据通过高带宽网络传输至云端,在云端进行计算后再传回本地.这种方式没有考虑到高能同步辐射光源数据处理的特殊性,无疑会显著增加传输和计算成本.给线站用户带来不必要的等待时间.随着边缘计算和云计算技术[95]的协同发展,未来光源的发展趋势是将探测器采集的数据进行线站端和云端协同处理.将一部分计算、存储资源部署在线站(边缘端),机器学习训练等海量计算任务依然发往云端计算集群进行处理.线站端每隔一段时间从云端获取关键的节点模型参数进行本地更新,基于端云协同技术完成本地轻量且重要的计算过程,从而达到快速实时 ptychography 图像处理和高效实验处理数据流.最终,通过ptychography 算法本身、人工智能、计算方法交叉研究促进ptychography 方法学的系统发展.