一种计算相位成像方法及系统与流程

1.本发明涉及计算成像技术领域，更具体地，涉及一种计算相位成像方法及系统。


背景技术：

2.计算显微成像技术是先进制造、精密测量及生物医学等领域中最为重要的一种先进成像方法。计算成像技术充分挖掘光与物质作用时成像光学多维度信息，诸如通过光谱、相位、偏振、相干、衍射等多纬度光学信息，借助数值计算方法和计算机超快计算能力，实现照明和样品的快速精密计算成像。计算成像方法突破了传统光场显微镜数值孔径下的衍射极限，可实现高精密、高速度计算显微成像。
3.在实际应用场景下，由于受光学探测器采样噪声、样品台位置扰动、照明光源指向稳定误差等综合扰动噪声的影响，用于计算成像的原始采样数据夹杂了诸多扰动噪声，计算成像分辨率和成像效率严重受损。近年以来，许多先进的计算成像方法逐步提出，并取得了广泛的科研关注，诸如叠层衍射成像算法、x射线断层扫描算法等。此类计算成像方法对采样微弱噪声具备一定的鲁棒性，但是当复杂噪声场景下，此类方法通常情况下难以收敛至全局最优解。


技术实现要素：

4.本发明针对现有技术中存在的技术问题，提供一种计算相位成像方法及系统，通过对计算成像迭代更新过程中施加动态激活函数，提高算法在复杂噪声场景下的收敛鲁棒性，提升算法跳出局部最优陷阱的能力，加速收敛至全局最优，由此解决了现有计算成像算法噪声鲁棒性低、收敛精度差等技术问题。
5.根据本发明的第一方面，提供了一种计算相位成像方法，包括：
6.步骤1，基于照明光束及所述照明光束入射到样品的当前待测位置，获取所述样品的近场出射波分布和远场衍射光强分布；
7.步骤2，基于所述远场衍射光强分布，对所述近场出射波分布进行迭代求解，获取每一次迭代更新后当前待测位置的样品函数和照明函数；
8.步骤3，基于所述每一次迭代更新后当前待测位置的样品函数和照明函数，计算每一次迭代更新后的出射波增量梯度残差ε，当ε小于设定阈值时，执行步骤4；否则，执行步骤2；
9.步骤4，基于动态正则函数对步骤2中获取的迭代更新后当前待测位置的样品函数和照明函数进行优化，基于优化后的样品函数和照明函数，重复执行步骤2至步骤4，直至所述出射波增量梯度残差ε和迭代更新后的远场衍射光强分布与步骤1中所述远场衍射光强分布之间的均方差mse均小于设定阈值，迭代终止，求得当前待测位置的样品函数和照明函数。
10.根据本发明的第二方面，提供一种计算相位成像系统，包括：
11.获取模块，用于基于照明光束及所述照明光束入射到样品的当前待测位置，获取
所述样品的近场出射波分布和远场衍射光强分布；
12.迭代求解模块，用于基于所述远场衍射光强分布，对所述近场出射波分布进行迭代求解，获取每一次迭代更新后当前待测位置的样品函数和照明函数；
13.计算模块，用于基于所述每一次迭代更新后当前待测位置的样品函数和照明函数，计算每一次迭代更新后的出射波增量梯度残差ε；
14.判断模块，用于当ε小于设定阈值时，执行所述优化模块；否则，执行所述迭代求解模块；
15.优化模块，基于动态正则函数对所述迭代更新模块获取的迭代更新后当前待测位置的样品函数和照明函数进行优化，基于优化后的样品函数和照明函数，重复执行所述迭代求解模块、所述计算模块和所述判断模块，直至所述出射波增量梯度残差ε和迭代更新后的远场衍射光强分布与步骤1中所述远场衍射光强分布的均方差mse均小于设定阈值，迭代终止，求得当前待测位置的样品函数和照明函数。
16.根据本发明的第三方面，提供了一种电子设备，包括存储器、处理器，所述处理器用于执行存储器中存储的计算机管理类程序时实现计算相位成像方法的步骤。
17.根据本发明的第四方面，提供了一种计算机可读存储介质，其上存储有计算机管理类程序，所述计算机管理类程序被处理器执行时实现计算相位成像方法的步骤。
18.本发明提供的一种计算相位成像方法及系统，采集照明光束入射样品后的包含噪声的远场衍射光强，采用相干衍射成像方法对样品和照明光束复振幅信息进行迭代重构，对迭代过程施加动态激活函数，用于抑制噪声对计算收敛的扰动，具备跳出局部最优陷阱能力，最终获得全局最优解。本发明解决了现有相位计算成像方法在噪声扰动下成像分辨率损失严重、算法不收敛等技术问题。
附图说明
19.图1为本发明提供的一种计算相位成像方法流程图；
20.图2为照明光束入射采集装置示意图；
21.图3为样品函数和照明函数的求解结果示意图；
22.图4为本发明提供的一种计算相位成像系统的结构示意图；
23.图5为本发明提供的一种可能的电子设备的硬件结构示意图；
24.图6为本发明提供的一种可能的计算机可读存储介质的硬件结构示意图。
25.附图中，各标号所代表的名称如下：
26.201、照明光束，202、样品，203、近场出射波分布，204、夫琅禾费远场衍射，205、远场衍射光强分布，301、样品函数的振幅，302、样品函数的相位，303、照明函数的振幅，304、照明函数的相位。
具体实施方式
27.为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。另外，
本发明提供的各个实施例或单个实施例中的技术特征可以相互任意结合，以形成可行的技术方案，这种结合不受步骤先后次序和/或结构组成模式的约束，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时，应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
28.基于现有技术中计算成像算法噪声鲁棒性低、收敛精度差等技术问题，本发明提出了一种计算相位成像方法，参见图1，该计算相位成像方法主要包括以下步骤：
29.步骤1，照明光束入射到样品待测位置j，获取样品的近场出射波分布ψ
jr
和远场衍射光强分布i
ju
，j表示所述待测位置的序号，j＝1,2,...,j,j为自然数。
30.作为实施例，所述步骤1，照明光束入射到样品当前待测位置j，获取样品的近场出射波分布ψ
jr
和采集样品的远场衍射光强分布i
ju
，包括：基于样品当前待测位置j的样品函数o
jr
和照明函数p
jr
，获取当前待测位置样品的近场出射波分布ψ
jr
，ψ
jr
＝p
jr
·qjr
；基于夫琅禾费衍射法采集远场衍射光强分布i
ju
。
31.可参见图2，为照明光束入射及光强采集装置示意图，其中，当照明光束201入射样品202时，对于样品202被照射的当前待测位置j，样品的远场衍射光强分布205为采集而来，具体的，基于夫琅禾费衍射204采集样品的远场衍射光强分布205。样品的近场出射波分布203是根据样品在当前待测位置的样品函数和照明函数计算而来。
32.其中，基于夫琅禾费衍射法采集样品的远场衍射光强分布，样品的远场衍射光强分布可表示如下：
[0033][0034]
其中，λ为照明光束的波长，z为夫琅禾费衍射光路中从样品到衍射光强探测器之间的距离，j为照明光束入射到样品的当前待测位置j位置，u和r分别是空间域和频域上的位置坐标矩阵，k是对应照明波长λ下的波矢。
[0035]
需要说明的是，样品的远场衍射光强分布是包含噪声的衍射光强，本发明实施例中，含噪声的衍射光强的获取形式主要有两种：第一种是照明光束透射样品，衍射光路和衍射光强探测器在样品背离照明光束一侧；第二种是照明光束反射样品，衍射光路和衍射光强探测器在入射照明的反射光路上。衍射光强中的噪声包括探测器噪声、样品抖动噪声和光源退相干噪声中的一种或多种。
[0036]
夫琅禾费衍射光场模型需满足夫琅禾费衍射近似条件，一般用菲涅尔数f进行表征：
[0037][0038]
其中，a是样品微结构孔径的尺寸，如半径，l是微结构与衍射光强探测器之间的距离。
[0039]
样品的近场出射波分布为根据样品函数和照明函数计算得到，对于初始近场出射波分布，可初始化样品函数和照明函数，然后根据初始化样品函数和照明函数，计算样品的初始近场出射波分布。其中，样品函数和照明函数的初始化步骤一般有两种形式：其一是采用随机数矩阵赋值样品函数和照明函数；其二是给定样品函数或照明函数先验信息赋值给
样品函数和照明函数。
[0040]
步骤2，基于样品的所述远场衍射光强分布，对所述近场出射波分布进行迭代求解，获取每一次迭代更新后当前待测位置的样品函数和照明函数。
[0041]
作为实施例，所述步骤2，基于样品的所述远场衍射光强分布，对所述近场出射波分布进行迭代求解，获取每一次迭代更新后当前待测位置的样品函数和照明函数，包括：根据所述远场衍射光强分布，基于时空间和傅里叶空间之间的交替投影算法，对当前所述近场出射波分布进行迭代更新，以获取更新后的近场出射波分布；根据当前更新前的近场出射波分布和更新后的近场出射波分布，基于相干衍射成像算法对当前迭代前的样品函数和照明函数进行更新，获取当前迭代更新后当前待测位置的样品函数和照明函数。
[0042]
可以理解的是，上述步骤可基于初始化的样品函数和照明函数获取照明光束入射到样品当前待测位置的初始近场出射波分布和远场衍射光强分布，然后结合样品的远场衍射光强分布，对近场出射波分布进行迭代求解，以对近场出射波分布进行更新。具体的，根据远场衍射光强分布，基于傅里叶投影变换，对当前近场出射波分布ψ
jr
进行更新，获取更新后的近场出射波分布ψ
′
jr
，其具体更新公式如下：
[0043][0044]
其中，o
jr
和p
jr
为当前更新前的样品函数和照明函数，ψ
jr
为更新前的近场出射波分布，fft表示快速傅里叶变换，fft-1
表示快速傅里叶逆变换；j表示照明光束p
jr
入射样品第j处(当前待测位置)的交叠位置，ψ
ju
表示出射波ψ
jr
的频域表示，也可称为迭代更新的远场衍射光强分布，ψ
′
jr
为当前迭代更新后的近场出射波分布，i
ju
为远场衍射光强分布。
[0045]
作为实施例，所述根据当前更新前的所述近场出射波分布和更新后的第一近场出射波分布，基于相干衍射成像算法对当前更新前的样品函数和照明函数进行更新，获取当前迭代更新后当前待测位置的样品函数和照明函数，包括：
[0046][0047]
其中：
[0048][0049]
其中，o
′
jr
和p
′
jr
分别为每一次迭代更新后的样品函数和照明函数，wo和w
p
分别是o
′
jr
和p
′
jr
的更新权重因子，取值范围为(0,1)的区间常数，o
jr
和p
jr
为更新前的样品函数和照明函数，ψ
jr
为更新前的近场出射波分布，ψ
′
jr
为更新后的近场出射波分布；ξ和η是函数sg(n,flag)的更新步长，n为当前迭代次数，n为出射波增量梯度残差ε小于设定阈值时的迭代次数，sg(n,flag)是激活函数，flag是函数sg(n,flag)的选择开关，有0和1两种状态，其初
值设定为1，随着迭代n的进行，当ε趋近于零时，flag值取反，依次循环进行。
[0050]
可以理解的是，本发明实施例采用相干衍射成像算法对含噪声的远场衍射光强分布i
ju
进行迭代求解，输出迭代更新后的样品函数o
jr
和照明函数p
jr
。
[0051]
具体的，相干衍射成像算法，其核心原理是基于时空间和傅里叶空间之间交替投影迭代优化理论，迭代过程可参见公式(4)和公式(5)，获取当前迭代更新后当前待测位置的样品函数和照明函数。
[0052]
步骤3，基于每一次迭代更新后当前待测位置的样品函数和照明函数，计算每一次迭代更新后的出射波增量梯度残差ε，当ε小于设定阈值时，执行步骤4；否则，执行步骤2。
[0053]
作为实施例，所述步骤3中计算每一次迭代更新后的出射波增量梯度残差ε，包括：获取当前待测位置之前的每一个待测位置在每一次迭代更新后的样品函数和照明函数；基于所有完成迭代更新的待测位置在每一次迭代更新后的样品函数和照明函数，以及所有完成迭代更新的待测位置更新前的近场出射波分布，计算每一次迭代更新后的出射波增量梯度残差ε。
[0054]
可以理解的是，当获取当前迭代更新后当前待测位置的样品函数和照明函数，计算当前的出射波增量梯度残差ε。其中，对于每一次迭代，需要计算当前待测位置之前的每一个位置迭代更新后的样品函数和照明函数。基于当前包括当前待测位置在内的之前所有位置迭代更新后的样品函数和照明函数，以及所有位置更新前的近场出射波分布，计算当前的出射波增量梯度残差ε。
[0055]
其中，出射波增量梯度残差ε的计算公式可如下表示：
[0056][0057]
其中，∑为求和函数。
[0058]
根据式(6)计算出每次迭代更新后的出射波增量梯度残差ε，判断出射波增量梯度残差ε是否小于设定阈值，如果是，则执行步骤4，如果不是，则再次执行步骤2，直到迭代更新后出射波增量梯度残差ε小于设定阈值，执行步骤4，并记录下每次出射波增量梯度残差ε小于设定阈值时的迭代次数n。
[0059]
步骤4，基于动态正则函数对步骤2中获取的迭代更新后当前待测位置的样品函数和照明函数进行优化，基于优化后的样品函数和照明函数，重复执行步骤2至步骤4，直至所述出射波增量梯度残差ε和迭代更新后的远场衍射光强分布与步骤1中所述远场衍射光强分布之间的均方差mse均小于设定阈值，迭代终止，求得当前待测位置的样品函数和照明函数。
[0060]
作为实施例，所述步骤4中，基于动态正则函数对步骤2中获取的当前待测位置的样品函数和照明函数进行更新，包括：根据获取的当前的当前待测位置的样品函数和照明函数、当前更新后的近场出射波分布和更新前的近场出射波分布，基于所述动态正则函数，对所述样品函数和所述照明函数进行更新。
[0061]
可以理解的是，当出射波增量梯度残差ε小于设定阈值时，记录下此时的迭代次数n，对于迭代次数n迭代更新后获取的当前待测位置的样品函数和照明函数施加动态正则函数，对迭代次数n对应的样品函数和照明函数进行优化，基于优化后的样品函数和照明函数，然后继续执行步骤2。
[0062]
其中，动态正则函数可表示为：
[0063][0064]
其中：
[0065][0066]
其中，o
′
jr
和p
′
jr
分别为更新后的第一样品函数和第一照明函数，ρ是指数滑动平均系数，r(x)是重塑函数，random是生成随机数矩阵函数。
[0067]
可以理解的是，随着计算迭代的进行，出射波增量梯度残差ε开始下降收敛，当ε小于设定阈值时，此时迭代次数赋值到sg(n,flag)激活函数n，标志位flag求反，重塑函数r(x)取值随机矩阵，通过所述步骤4施加的动态正则函数，可以理解为，当迭代次数为n时，可以利用公式(7)和(8)替换公式(4)和(5)来求取迭代次数为n时的样品函数和照明函数。在其它的迭代次数时，利用公式(4)和(5)来求取样品函数和照明函数。
[0068]
对步骤2中迭代更新获取的样品当前待测位置的样品函数和照明函数，若当前迭代次数的出射波增量梯度残差ε小于设定阈值时，基于步骤4中的动态正则函数对当前的样品函数和照明函数进行正则化优化。
[0069]
基于正则化优化后的样品函数和照明函数，循环执行步骤2-步骤4进行循环迭代，得到最终样品当前待测位置的样品函数和照明函数。
[0070]
需要说明的是，步骤4的动态正则函数中引入重塑函数r(x)。当且仅当残差ε小于设定阈值时，表示迭代过程陷入局部最优，此时样品函数或照明函数引入重塑函数随机矩阵random(x)，将样品函数或照明函数进行随机矩阵赋值，强行打破局部最优陷阱，使其跳出局部最优，直至样品函数迭代更新寻找到全局最优过程。
[0071]
其中，作为实施例，所述重复执行步骤2至步骤4，最终求得当前待测位置的样品函数和照明函数，包括：基于优化后的样品函数和照明函数，重复执行步骤2至步骤4，直至所述近场出射波分布增量梯度残差ε和远场衍射光场与步骤1中所述远场衍射光强分布的均方差mse均小于设定阈值，迭代终止，求得当前待测位置的样品函数和照明函数，包括：基于优化后的样品函数和照明函数，重复执行步骤2至步骤4，当满足迭代终止条件时，获取当前待测位置的样品函数和照明函数；其中，当所述出射波增量梯度残差ε和迭代更新后的远场衍射光强分布之间的均方差mse均小于设定阈值时，判定满足迭代终止条件，其中：
[0072][0073]
其中，j表示所有待测位置数量；
[0074]
或者，若当前迭代次数达到设定的最大迭代次数，判定当前迭代满足迭代终止条件。
[0075]
可以理解的是，步骤4基于动态正则函数对样品函数和照明函数正则化优化后，基于优化后的样品函数和照明函数，然后循环执行步骤2、步骤3和步骤4，直到满足迭代终止条件，停止迭代。其中，本发明迭代终止条件的判定方式主要有两种：第一种是照明光束根据重构的近场出射波分布通过夫琅禾费衍射传播计算出的远场光强分布ψ
ju2
与步骤1输入的衍射光强i
ju
之间求均方差mse,用均方差mse结合当前的出射波增量梯度残差ε综合进行判定，比如，当均方差mse和出射波增量梯度残差ε均小于对应的设定阈值时，满足迭代终止条件；第二种是根据先验知识设定最大迭代循环次数进行判定，当迭代次数到达最大迭代循环次数时，停止迭代。当满足迭代终止条件时，输出样品当前待测位置的样品函数和照明函数，其中，图3为样品函数和照明函数的最终求解结果示意图，具体为输出样品函数的振幅301和相位302，以及输出照明函数的振幅303和相位304。其中，所求得的样品函数是一个复振幅矩阵，此样品函数通过计算相位成像的方法求得，是一种先进的计算成像方法，从方法论上讲完全区别于传统的几何光路追迹的光场成像方法。
[0076]
本发明采集照明光束入射样品后的包含噪声的远场衍射光强，采用相干衍射成像方法对样品和照明光束复振幅信息进行迭代重构，对迭代过程施加动态正则函数，用于抑制噪声对计算收敛的扰动，具备跳出局部最优陷阱能力，最终获得全局最优解，解决了现有相位计算成像方法在噪声扰动下成像分辨率损失严重、算法不收敛等技术问题。
[0077]
图4为本发明实施例提供的一种计算相位成像系统结构图，如图4所示，一种计算相位成像系统，包括获取模块401、迭代求解模块402、计算模块403、判断模块404和优化模块405，其中：
[0078]
获取模块401，用于基于照明光束及所述照明光束入射到样品的当前待测位置，获取所述样品的近场出射波分布和远场衍射光强分布；
[0079]
迭代求解模块402，用于基于所述远场衍射光强分布，对所述近场出射波分布进行迭代求解，获取每一次迭代更新后当前待测位置的样品函数和照明函数；
[0080]
计算模块403，用于基于所述每一次迭代更新后当前待测位置的样品函数和照明函数，计算每一次迭代更新后的出射波增量梯度残差ε；
[0081]
判断模块404，用于当ε小于设定阈值时，执行所述优化模块405；否则，执行所述迭代求解模块402；
[0082]
优化模块405，基于动态正则函数对所述迭代更新模块获取的迭代更新后当前待测位置的样品函数和照明函数进行优化，基于优化后的样品函数和照明函数，重复执行所述迭代求解模块402、所述计算模块403和所述判断模块404，直至所述出射波增量梯度残差ε和迭代更新后的远场衍射光强分布与步骤1中所述远场衍射光强分布之间的均方差mse均小于设定阈值，迭代终止，求得当前待测位置的样品函数和照明函数。
[0083]
可以理解的是，本发明提供的一种计算相位成像系统与前述各实施例提供的计算相位成像方法相对应，计算相位成像系统的相关技术特征可参考计算相位成像方法的相关技术特征，在此不再赘述。
[0084]
请参阅图5，图5为本发明实施例提供的电子设备的实施例示意图。如图5所示，本发明实施例提了一种电子设备500，包括存储器510、处理器520及存储在存储器510上并可
在处理器520上运行的计算机程序511，处理器520执行计算机程序511时实现所述计算相位成像方法的步骤。
[0085]
请参阅图6，图6为本发明提供的一种计算机可读存储介质的实施例示意图。如图6所示，本实施例提供了一种计算机可读存储介质600，其上存储有计算机程序611，该计算机程序611被处理器执行时实现计算相位成像方法的步骤。
[0086]
本发明实施例提供的一种计算相位成像方法及系统，主要具有以下有益效果：
[0087]
1、本发明通过对计算成像迭代更新过程中施加动态正则激活函数，提高算法在复杂噪声场景下的收敛鲁棒性，提升算法跳出局部最优陷阱的能力，加速收敛至全局最优。
[0088]
2、本发明所提供的计算相位成像方法具备系统精简、算法鲁棒、求解快速、成像分辨率高等优势，能够实现多场景、复杂噪声干扰下样品的快速精密计算成像求解。
[0089]
需要说明的是，在上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详细描述的部分，可以参见其它实施例的相关描述。
[0090]
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
[0091]
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式计算机或者其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0092]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0093]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0094]
尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。
[0095]
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包括这些改动和变型在内。