CN107025637B - 基于贝叶斯估计的光子计数集成成像迭代重构方法 - Google Patents

Abstract

本发明公开了一种基于贝叶斯估计的光子计数集成成像迭代重构方法，根据光子计数过程的泊松分布以及图像空间相关性，建立光子数多点估计的后验概率模型，在迭代过程中引入泊松去噪算法来更新后验概率模型的参数。本发明实现光子计数集成成像系统在极少量光子数探测情况下的成像，提高目标三维重构质量。

Description

基于贝叶斯估计的光子计数集成成像迭代重构方法

技术领域

本发明属于光子计数成像领域，具体涉及一种基于贝叶斯估计的光子计数集成成像迭代重构方法。

背景技术

集成成像由于能够记录三维目标不同视角的2D图像，捕捉光线的方向和强度信息，因而在三维成像领域得到了广泛的关注，同时在目标识别、深度估计以及遮拦物去除等多个应用方向体现其独特的优势。集成成像可以借助微透镜阵列或者相机阵列来实现，而合成孔径集成成像系统则是通过单个相机的平移对三维场景进行多通道观测，每个通道即代表每次相机对场景的图像采集，采集得到的二维图像则被称为元素图像。

光子计数成像系统摆脱了传统成像系统需要大量光子的局限性，能够在微弱光条件下实现少量光子成像。由于光子计数过程的泊松分布特性，修复光子计数图像的问题变成了对泊松分布期望值的有效估计，而集成成像对于场景数据采集的冗余性恰巧弥补了单张光子计数图像采样不足的劣势，因此借助集成成像的系统平台，光子计数技术能够更好地应用在微弱光成像领域。

在光子计数集成成像重构领域中，通常采用泊松分布来建立光子计数探测模型，传统方法是将极大似然估计(MLE)用于集成成像的深度重构，极大似然估计算法尽管能够产生比较好的结果，但缺乏目标物体的先验信息，在重建效果方面仍有不足，当环境光子低到一定程度时，其应用则会受到限制。因此，人们提出采用Gamma分布作为元素图像的先验分布，对场景进行贝叶斯估计重构，然而Gamma分布并不是描述目标特性最恰当的分布，它只能作为一种粗略描述。而另一种重构方向是将集成成像的重构问题转变成求逆问题，构造成像系统矩阵H，引入正则化因子进行罚似然最大期望值(PMLEM)估计重构，该方法由于达到最优结果的条件时迭代式中惩罚项的选择以及参数调节的不确定性，从而导致该方法场景自适应性不强。

发明内容

本发明的目的在于提供一种基于贝叶斯估计的光子计数集成成像迭代重构方法，实现极微弱光环境下目标高质量三维成像的方法。

实现本发明的技术解决方案为：一种基于贝叶斯估计的光子计数集成成像迭代重构方法，步骤如下：

第一步，通过光子计数泊松过程的仿真，得到集成成像的光子计数元素图像C_kl，建立光子数估计的贝叶斯后验概率模型；

第二步，根据贝叶斯理论，通过第一步贝叶斯后验概率模型的后验概率及其先验分布，计算贝叶斯后验概率模型的后验均值作为像素光子数的贝叶斯估计值

第三步，利用第二步元素图像像素光子数的贝叶斯估计值进行目标图像重构，即通过虚拟的小孔阵列将元素图像按比例放大，其放大比例为M＝z/g，由重构目标距离相机透镜的距离z以及相机透镜与传感器的距离g决定；放大翻转后的元素图像在同一平面上相互叠加，从而重构出了z深度的贝叶斯重构的目标图像而虚化了不属于该深度位置的目标；

第四步，对第三步得到的重构的目标图像进行泊松去噪处理，得到深度切片图像利用深度切片图像以及元素图像C_kl更新参数ζ和η，回到第二步代入计算，若深度切片图像与前一次该步骤的去噪结果均方误差小于设定阈值ε，则步骤终止，将第三步的作为重构结果；如果不小于设定阈值ε，利用以及元素图像C_kl更新参数ζ和η，并回到第二步代入计算，直到达到设置的迭代次数。

本发明与现有技术相比，其显著优点:(1)采用迭代的方式更新贝叶斯估计的参数，更接近光子数探测的真实值。(2)基于图像空间相关性，利用邻域像素建立多点贝叶斯估计模型。(3)引入泊松去噪算法，更好地估计了像素间的强度关系，优化了贝叶斯估计的参数。

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

附图说明

图1是合成孔径集成成像系统示意图。

图2是正常光照条件下集成成像7*8幅元素图像。

图3是正常光照下z＝97cm处目标重构图像。

图4(a)列是Np＝100、150、200、300、500五种情况下的部分光子计数元素图像；(b)列是Np＝100、150、200、300、500五种情况下的极大似然估计重构的目标图像；(c)列是Np＝100、150、200、300、500五种情况下基于贝叶斯估计的光子计数集成成像迭代重构结果。

图5是基于贝叶斯估计的光子计数集成成像迭代重构方法原理示意图。

具体实施方式

结合图5，本发明基于贝叶斯估计的光子计数集成成像迭代重构方法，步骤如下：

第一步，通过光子计数泊松过程的仿真，得到集成成像的光子计数元素图像C_kl，建立光子数估计的贝叶斯后验概率模型；

第二步，根据贝叶斯理论，通过第一步贝叶斯后验概率模型的后验概率及其先验分布，计算贝叶斯后验概率模型的后验均值作为像素光子数的贝叶斯估计值

第三步，利用第二步元素图像像素光子数的贝叶斯估计值进行目标图像重构，即通过虚拟的小孔阵列将元素图像按比例放大，其放大比例为M＝z/g，由重构目标距离相机透镜的距离z以及相机透镜与传感器的距离g决定；放大翻转后的元素图像在同一平面上相互叠加，从而重构出了z深度的贝叶斯重构的目标图像而虚化了不属于该深度位置的目标；

第四步，对第三步得到的重构的目标图像进行泊松去噪处理，得到深度切片图像利用深度切片图像以及元素图像C_kl更新参数ζ和η，回到第二步代入计算，若深度切片图像与前一次该步骤的去噪结果均方误差小于设定阈值ε(ε根据场景的差异需要人为设置)，则步骤终止，将第三步的作为重构结果；如果不小于设定阈值ε，利用以及元素图像C_kl更新参数ζ和η，并回到第二步代入计算，直到达到设置的迭代次数，迭代次数一般为2-5次。

在第一步中，建立光子数估计的贝叶斯后验概率模型为：

其中符号的上标p、p’代表像素在元素图像上的位置；下标k、l代表元素图像在x、y轴方向上的位置，取值范围分别为0～K-1、0～L-1；为光子计数元素图像上p点的探测值，为p点相邻像素p’点的探测值，则是p点光子数的期望值，其先验分布服从形状参数α和尺度参数β的Gamma分布。

在第二步中，计算后验均值作为像素光子数的贝叶斯估计值如下式：

定义形状参数尺度参数η＝(2+β)^-1作为迭代算法参数，迭代初始值令参数ζ＝0，η＝1，整个迭代过程中ζ和η进行更新。

在第四步中，利用以及元素图像C_kl更新参数ζ和η的步骤如下：

(1)由于元素图像的稀疏性，认为以p点为中心的邻域空间U(p)(U(p)默认表示元素图像上以像素位置p为中心的3*3邻域空间)里至多一个像素值不为0；因此，对于不为0的元素图像像素点保留其原值，即ζ＝0，η＝1；

(2)对于为0的像素点利用其邻域像素点进行估计，更新公式如下：

其中指去噪后的深度切片图像上与元素图像上相对应的目标像素点，同理，是与对应的目标像素点。

实施例

结合图5，本发明基于贝叶斯估计的光子计数集成成像迭代重构方法，步骤如下:

第一步，在正常光照条件下，对目标场景进行2-D图像采集，得到一系列元素图像如图2。如图1所示，在正常光照条件下，以S_x＝S_y＝5mm的间隔沿着x、y轴移动相机，对目标场景进行2D横向采样，获得7*8幅、尺寸为200*153的强度元素图像，如图2。

第二步，引入可调节的单幅元素图像光子总数Np，由强度元素图像仿真光子计数元素图像，计算公式如下：

C_x～poisson(N_pλ_x)

其中I_x为强度元素图像，下标x代表元素图像上的像素位置，N为元素图像总像素数，C_x为光子计数元素图像上x位置的泊松分布采样，该位置泊松分布的期望值为N_pλ_x。当Np＝100、150、200、300、500时，部分光子计数元素图像如图4(a)列所示。

第三步，根据第二步得到的光子计数元素图像C_kl，建立光子数估计的贝叶斯后验概率模型：

其中符号的上标p代表像素在元素图像上的位置；下标k、l代表元素图像在x、y轴方向上的位置，取值范围分别为0～K-1、0～L-1；f为后验概率。为光子计数元素图像上p点的探测值，为p点相邻像素p’点的探测值。则是p点光子数的期望值，其先验分布服从形状参数为α和尺度参数为β的Gamma分布。

第四步，计算第三步贝叶斯后验均值作为像素光子数的贝叶斯估计如下式：

定义形状参数尺度参数η＝(2+β)^-1作为本发明的迭代算法参数，迭代初始值令参数ζ＝0，η＝1，迭代过程中ζ和η进行更新。

第五步，利用第四步元素图像像素光子数的估计值进行目标图像重构。通过虚拟的小孔阵列将元素图像以一定比例放大，其放大因子为M＝z/g，由重构目标距离相机透镜的距离z以及相机透镜与传感器的距离g决定；放大翻转后的元素图像在同一平面上相互叠加，从而重构出了z深度的目标图像而虚化了不属于该深度位置的目标，计算公式如下：

已知相机透镜与传感器的间隔g＝4cm，场景中的目标离相机透镜的距离z＝97cm，这个深度的目标重构图像如图3。在元素图像光子总数Np＝100、150、200、300、500的条件下分别进行光子计数集成成像重构。

根据步骤一集成成像系统的采集参数，元素图像的放大因子M＝z/g＝24.25，元素图像的排布间隔S_x＝S_y＝5mm，元素图像在x、y轴方向的采样数K＝8、L＝7。

第六步，对第五步得到的贝叶斯重构图像进行泊松去噪处理，这里的泊松去噪算法采用基于方差稳定变换的三维块匹配(BM3D)算法。若去噪后的深度切片图像与前一次该步骤的去噪结果均方误差小于设定阈值ε＝0.01，则步骤终止，将第五步骤的作为本发明的重构结果；否则，利用以及元素图像C_kl更新参数ζ和η，并回到第四步。参数更新方法如下：

1.由于元素图像的稀疏性，认为以p点为中心的邻域空间U(p)(U(p)默认表示元素图像上以像素位置p为中心的3*3邻域空间)里至多一个像素值不为0。因此，对于不为0的元素图像像素点保留其原值，即ζ＝0，η＝1；

2.对于为0的像素点利用其邻域像素点进行估计，更新公式如下：

其中指去噪后的深度切片图像上与元素图像上相对应的目标像素点，同理，是与对应的目标像素点。

图4(b)为极大似然估计在z深度重构的目标图像。采用提出的基于贝叶斯估计的光子计数集成成像迭代重构方法对z深度的目标重构结果如图4(c)所示。将图4(b)、(c)的重构图像结合表1列出的不同Np条件下重构图像峰值信噪比数据，可以看出随着Np的变化，本发明重构图像的峰值信噪比呈现递增性；与MLE进行对比，当Np较小时，迭代重构算法优于MLE的差距并不大，而随着Np的增加，这种差距逐渐增大，本发明体现出显著的优势性。当Np＝500时，本发明提出的算法高出MLE达到3.3dB，图像目标的轮廓以及细节更加清晰。表1不同Np条件下的光子计数重构图像峰值信噪比

Np	MLE	本发明
			100	13.1435dB	13.9796dB
150	13.2019dB	14.5406dB
			200	13.2541dB	14.9824dB
300	13.3567dB	15.7778dB
			500	13.5626dB	16.8864dB

Claims (3)

1.一种基于贝叶斯估计的光子计数集成成像迭代重构方法，其特征在于步骤如下：

第一步，通过光子计数泊松过程的仿真，得到集成成像的光子计数元素图像C_kl，建立光子数估计的贝叶斯后验概率模型；

第二步，根据贝叶斯理论，通过第一步贝叶斯后验概率模型的后验概率及其先验分布，计算贝叶斯后验概率模型的后验均值作为像素光子数的贝叶斯估计值

第三步，利用第二步元素图像像素光子数的贝叶斯估计值进行目标图像重构，即通过虚拟的小孔阵列将元素图像C_kl按比例放大，其放大比例为M＝z/g，由重构目标距离相机透镜的距离z以及相机透镜与传感器的距离g决定；放大翻转后的元素图像C_kl在同一平面上相互叠加，从而重构出了z深度的贝叶斯重构的目标图像而虚化了不属于该深度位置的目标；

第四步，对第三步得到的重构的目标图像进行泊松去噪处理，得到深度切片图像利用深度切片图像以及元素图像C_kl更新参数ζ和η，回到第二步代入计算，若深度切片图像与前一次该步骤的去噪结果均方误差小于设定阈值ε，则步骤终止，将第三步的作为重构结果；如果不小于设定阈值ε，利用以及元素图像C_kl更新参数ζ和η，并回到第二步代入计算，直到达到设置的迭代次数；

利用以及元素图像C_kl更新参数ζ和η的步骤如下：

(1)由于元素图像的稀疏性，认为以p点为中心的邻域空间U(p)里至多一个像素值不为0；因此，对于不为0的元素图像像素点保留其原值，即ζ＝0，η＝1；

(2)对于为0的像素点利用其邻域像素点进行估计，更新公式如下：

其中指去噪后的深度切片图像上与元素图像上相对应的目标像素点，同理，是与对应的目标像素点。

2.根据权利要求1所述的基于贝叶斯估计的光子计数集成成像迭代重构方法，其特征在于第一步中，建立光子数估计的贝叶斯后验概率模型为：

其中符号的上标p、p′代表像素在元素图像上的位置；下标k、l代表元素图像在x、y轴方向上的位置，取值范围分别为0～K-1、0～L-1；为光子计数元素图像上p点的探测值，为p点相邻像素p’点的探测值，则是p点光子数的期望值，其先验分布服从形状参数α和尺度参数β的Gamma分布。

3.根据权利要求1所述的基于贝叶斯估计的光子计数集成成像迭代重构方法，其特征在于第二步中，计算后验均值作为像素光子数的贝叶斯估计值如下式：

定义参数参数η＝(2+β)^-1作为迭代算法参数，迭代初始值令参数ζ＝0，η＝1，迭代过程中ζ和η进行更新。