针对非均匀性光强下光强传输方程的无边界误差求解方法与流程

本发明属于光学测量中的相位恢复与定量相位成像技术，特别是一种非均匀性光强下光强传输方程的无边界误差求解方法。

背景技术：

相位恢复是光学测量与成像技术的一个重要课题，无论在生物医学还是工业检测领域，相位成像技术都在发挥着重要的作用。纵观光学测量近半个世纪的进展，最经典的相位测量方法是干涉测量法。然而，干涉测量法的缺点也十分明显：干涉测量一般需要高度相干性的光源(如激光)，从而需要较为复杂的干涉装置；引入额外的参考光路会使光学系统对测量环境的要求变得十分苛刻；高相干性的光源引入的散斑相干噪声会限制成像系统的空间分辨率与测量精度。

不同于干涉测量法，还有一类非常重要的相位测量技术并不需要借助干涉，它们统称为相位恢复。由于可以直接测量光波场的振幅/强度而无法直接测量光波场的相位分布，所以可以将由强度分布来恢复(估算)相位这一过程考虑为一个数学上的“逆问题”，即相位恢复问题。相位恢复方法还可细分为迭代法与直接法。光强传输方程法是相位恢复方法中的一种典型的直接法。光强传输方程是一个二阶椭圆偏微分方程，其阐明了沿着光轴方向上光强度的变化量与光轴垂直的平面上光波的相位的定量关系。在光强轴向微分以及光强分布已知的情况下，通过数值求解光强传输方程可直接获取相位信息。相比于干涉法，直接求解光强传输方程获取相位信息主要优点包括：(1)非干涉，仅仅通过测量物面光强直接求解相位信息，不需要引入额外参考光；(2)非迭代，通过直接求解微分方程获得相位；(3)可以很好的应用于白光照明，如传统明场显微镜中的科勒照明(illumination)；(4)无需相位解包裹，直接获取相位的绝对分布，不存在一般干涉测量中的2π相位包裹问题；(5)无须复杂的光学系统，对于实验环境没有苛刻的要求，振动不敏感。

针对光强传输方程的求解，目前已有许多方法提出：如格林函数法(m.reedteague，"deterministicphaseretrieval:agreen’sfunctionsolution,"j.opt.soc.am.73,1434-1441(1983).)、泽尼克多项式展开法(t.e.gureyevandk.a.nugent,"phaseretrievalwiththetransport-of-intensityequation.ⅱ.orthogonalseriessolutionfornonuniformillumination,"j.opt.soc.am.a13,1670-1682(1996).)、快速傅里叶变换法(l.j.allenandm.p.oxley,"phaseretrievalfromseriesofimagesobtainedbydefocusvariation,"optcommun199,65-75(2001).)。然而，在这些求解光强传输方程的算法中，往往需要引入teague辅助函数，从而将光强传输方程转换为两个标准的泊松方程简化求解。然而，teague辅助函数隐含一个较强的假设，即光强的横向能流场是个保守场。一般情况下该假设并不成立，所以teague辅助函数必定会引起求解误差(称为相位差异)，导致传统求解方法并不能给出光强传输方程的精确解(j.a.schmalz,t.e.gureyev,d.m.paganin,andk.m.pavlov,"phaseretrievalusingradiationandmatter-wavefields:validityofteague'smethodforsolutionofthetransport-of-intensityequation,"phys.rev.a84,023808(2011))，尤其是在边界处，解得的相位值误差会更加明显。所以需要针对此问题(光强传输方程解得的相位存在误差；光强非均匀时，光强传输方程可能无解)采取有效措施，将会大大影响光强传输方程法求解相位的准确度，使其难以应用在高精度相位测量以及定量相位成像领域。

技术实现要素：

本发明的目的在于提供一种针对非均匀性光强下光强传输方程的无边界误差求解方法，以解决在非均匀性光强下的精确的相位恢复问题。

实现本发明目的的技术解决方案为：一种针对非均匀性光强下光强传输方程的无边界误差求解方法，首先初始化光强的轴向微分以及相位值，然后计算得到的非精确的相位值以及其所对应的光强轴向微分值，获得当前迭代后的光强轴向微分值与上一轮的光强轴向微分值之间的残差，每次迭代完成后，判断光强轴向微分误差以及所对应的相位残差是否足够小，当满足停止迭代条件时，得到的相位值就是所求的精确相位值，能够准确的求解在非均匀光强下的光强传输方程，稳定并且精确地获得待测物体相位。

本发明与现有技术相比，其显著优点：(1)能够高效准确地求解出相位值，降低了传统方法采用teague辅助函数所引起的求解误差(称为相位差异)，尤其是在边界处求得的相位误差。(2)在每次求解光强轴向微分误差项所对应的相位差异时，能够避免光强传输方程无解的情况，这样可以提高算法的鲁棒性，降低非均匀光强对相位求解过程的影响。(3)在非精确求解与光强轴向微分误差相对应的相位差异时，利用常数替代非均匀光强分布，减少以及的计算次数，降低了算法的复杂度，从而减少单次迭代循环运算量。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是本发明的无边界误差求解方法步骤流程示意图。

图2(a)是仿真的非均匀的光强图。

图2(b)是仿真的理论的相位图。

图2(c)是仿真的光阑图。

图3(a)是快速傅里叶变换算法求解得到的原始相位图。

图3(b)是本方法求解得到的原始相位图。

图3(c)是快速傅里叶变换算法求解得到的原始相位误差图。

图3(d)是本方法求解得到的原始相位误差图。

具体实施方式

本发明针对非均匀性光强下光强传输方程的无边界误差求解方法，首先初始化光强的轴向微分以及相位值，然后计算得到的非精确的相位值以及其所对应的光强轴向微分值，同时可以获得本轮光强轴向微分值与上一轮的光强轴向微分值之间的残差。每次迭代完成后，判断光强轴向微分误差以及所对应的相位残差是否足够小(达到可以忽略不计)，当满足停止迭代条件时，得到的相位值就是所求的精确相位值，最终能够准确的求解在非均匀光强下的光强传输方程，稳定并且精确地获得待测物体相位。

如图1所示，本发明针对非均匀性光强下光强传输方程的无边界误差求解方法，具体实现步骤如下：

步骤一，初始化光强的轴向微分以及相位值。

采集欠焦、聚焦、过焦三幅光强分布，分别记作i+、i0、i-，其中i+与i-的离焦距离δz相等且方向相反，光强轴向微分信号由中心有限差分法估计得到：

至此完成了非均匀性光强下光强传输方程的高效准确求解的数据准备工作，即聚焦光强分布i0，首次迭代的光强轴向微分估计δj0，初始化相位信息φ0＝0，迭代次数n＝1。

步骤二，非精确的相位δφn求解，利用离散余弦变换或者快速傅里叶变换求解光强传输方程，得到一个非精确的相位δφn，关键在于求解光强传输方程非精确解时，利用一个常数值代替非均匀光强分布矩阵，从而有效地解决了非均匀光强下光强传输方程可能无解的问题，而且极大地简化了非精确相位δφn的求解过程：

首先，求解光强传输方程利用常数imax取代i，imax为与拍摄到的聚焦位置光强图即i0中的最大值，这样可以有效地避免非均匀光强情况下i0存在0值的情况，解决方程可能无解的问题；而且由于imax为常数值(非精确估计)，则光强传输方程可以化简为所以其中为逆拉普拉斯运算符号，k为波数(λ为照明光的波长)。在光强传输方程中传统方法i为拍摄到的聚焦位置光强图即i0，引入teague辅助函数()，首先求解得到ψ，即然后求解整个过程需要分别进行两次以及运算。本方法中，i取值为一个常数，该常数为i0的最大值，记作imax，既可以有效地避免非均匀光强情况下i0存在0值的情况，导致无法求解从而导致φ无解，又可以将光强传输方程化简为从而快速求解得到整个过程只需要进行一次运算。(为哈密顿算符，为逆拉普拉斯运算符号，k为波数)。

基于傅里叶变换或者离散余弦变化，求解出非精确相位φ即δφn轴向微分为步骤四中解得的轴向微分误差δjn-1。上述非精确解φ的求解转变为光强轴向微分的估计，以及的计算。先进行光强轴向微分估计：当迭代次数n＝1时，轴向微分值取值当n>1时，轴向微分取值为δjn＝jn-1-jn，即相邻两次迭代的轴向微分的残差，jn的取值取决于步骤三。然后快速求解通过快速傅里叶变换进行实现：

其中，代表傅里叶变换，代表傅里叶逆变换，(u,v)是与空间坐标相对应的频域坐标。

通过离散余弦变换进行实现：

其中，dct代表离散余弦变换，dct^-1代表离散逆余弦变换，(u,v)是与空间坐标相对应的频域坐标。对于二维数据矩阵f(x,y),x＝0,1,...,n-1,y＝0,1,...,m-1，其中，(x,y)代表二维空间坐标，n,m为图像的二维尺寸，其离散余弦正弦逆变换可以定义为：

式中常系数au与av分别为

步骤三，求解δφn对应的轴向微分值jn，将非精确解δφn重新代入光强传输方程的右端，得到人为计算的光强轴向微分

将非精确解φ(当代入光强传输方程左侧的光强轴向微分值记为δjn时，φ记作δφn)，即步骤二中φ的求解可以记作的步骤二求得的相位值记作δφn，代入光强传输方程的右端，得到人为计算的光强轴向微分

由于上述δφn的求解是非精确的(i近似估计为imax，而且对于利用辅助函数求解光强传输方程，本身就隐含一个较强的假设，即光强的横向能流场是个保守场。一般情况下该假设并不成立，所以辅助函数必定会引起求解误差)，其对应的实际轴向微分值为并不等于求解δφn时所代入方程中的轴向微分值δjn-1。

通过快速傅里叶变换进行实现：

其中，代表傅里叶变换，代表傅里叶逆变换，(u,v)是与空间坐标相对应的频域坐标。

通过离散余弦变换进行实现：

其中，dct代表离散余弦变换，dct^-1代表离散逆余弦变换(dct与dct^-1的具体求解过程见步骤二)，(u,v)是与空间坐标相对应的频域坐标。

步骤四，获得上一次求解得到的光强轴向微分jn-1与本轮光强轴向微分jn之间的误差函数δjn：将步骤二中得到的非精确的相位δφn(即相位误差补偿项)与上一轮迭代中获得的非精确解φn-1相加，得到本轮非精确相位值φn。

具体地，求解相邻两次的光强轴向微分的误差函数δjn＝jn-1-jn时，n代表目前迭代的次数即n＝1,2....，(j0为初始化的光强轴向微分，)；步骤二中求得的对应于光强轴向微分值δjn-1时的相位解即为相位误差补偿项其与上一轮迭代得到的非精确解φn-1相加得到新一轮的非精确解φn＝φn-1+δφn(φ0＝0为初始化相位值)。

步骤五，判断迭代终止条件。判断循环条件，不满足则一直重复步骤二到步骤四，一直到满足终止条件为止，最终获得物体的相位图。

光强轴向微分的误差函数δjn小于一个给定的阈值(阈值εj<0.1)，相位误差补偿项δφn小于一个给定的阈值(阈值εφ<0.1)，或者超过最大迭代次数(限定的最大迭代次数nmax≥2)，三者都不满足时，迭代次数n←n+1(n的取值从1开始取)，三者任意一者满足时可终止迭代，最终得到了光强传输方程的精确解，也就是待求的相位分布φn，即物体的相位图。

通过上述步骤可以看出，本发明采用迭代补偿的方法，能够高效准确地求解出相位值，降低了传统方法采用teague辅助函数所引起的求解误差(称为相位差异)，尤其是在边界处求得的相位误差。并且在每次求解光强轴向微分误差项所对应的相位差异时，该方法能够避免光强传输方程无解的情况，这样可以提高算法的鲁棒性，降低非均匀光强对相位求解过程的影响；此外，在非精确求解与光强轴向微分误差相对应的相位差异时，利用常数代替聚焦面非均匀的光强信息，减少了单次循环中的以及的运算次数，降低了算法的复杂度。

为了测试一种针对非均匀性光强下光强传输方程的无边界误差求解方法，采用一组仿真实例来验证本方法的有效性，在该仿真中，待测光波场的振幅被定义在256×256的像素网格上，像元尺寸为2μm×2μm，波长为633nm。如图2(a)、图2(b)所示，相位与光强图被定义为i(x,y)＝exp(-a0x²-b0y²)/(2σ²)，φ(x,y)＝a1x³+b1y³+a2x⁶+b2y⁶+0.82，其中，a0＝1,b0＝1,σ＝1.8×10^-4，a1＝-0.7,b1＝2,a2＝10,b2＝-10，且光场中存在一个不规则的光阑，光阑形状如图2(c)。利用恢复出的相位与理论相位之间的误差值去客观地衡量相位恢复的精度。图3(a)-图3(d)显示了该算法的仿真结果。

图3(a)显示了快速傅里叶变换算法求解得到的原始相位(其中为了避免光强传输方程无解的情况，将聚焦图i0中光强值为0的值近似为一个趋于0的数)，对应的图3(c)中显示了快速傅里叶变换算法求解得到的原始相位的误差，由于可获得精确解的充要条件在大多数位置都无法满足，所以可以明显看出很大幅度的相位差异，尤其是光阑边界处的较大相位差异，这个误差水平对于定量相位恢复来说是不可接受的。经过提出的算法求解得到的相位如图3(b)所示，对应求解得到的原始相位的误差如图3(d)，其残余误差明显下降，残余误差已经减小到了一个可以忽略的水平。从图3(b)和图3(d)中可以看出，使用本发明能够稳定地求解得到相位值，且与实际真实的相位值基本吻合，误差值基本可以忽略，说明本发明能够实现非均匀光强下光强传输方程的无边界误差求解。