CN102759408B - 一种单光子计数成像系统及其方法 - Google Patents

Abstract

本发明涉及一种单光子计数成像系统及其方法，该系统包括：滤光片、第一透镜(1)、DMD控制系统、第二透镜(2)、单光子计数器和数据处理单元；其中，DMD结合第一透镜(1)和第二透镜(2)，用于将二维图像数据转化为一维数据序列完成被测信号的压缩采样，极弱光通过滤光片滤除杂光，经第一透镜(1)在DMD控制系统处成像，并由DMD控制系统控制光子被反射到第二透镜(2)的概率，经过第二透镜(2)控制光子聚焦；数据处理单元结合单光子计数器完成稀疏重建，数据处理单元根据单光子计数器在一定时间内对光子进行计数折算成探测到光子数的概率作为测量值，和DMD控制系统上的测量矩阵经过最优化算法重建光子密度图像，解算出二维图像。

Description

一种单光子计数成像系统及其方法

技术领域

本发明涉及极弱光探测的技术领域，特别涉及一种单光子计数成像系统及其方法，它采用压缩传感理论和DLP技术，利用点探测器可以实现极弱光对象的高质量二维成像。

背景技术

普通成像设备通过记录观察对象上某点的光强和位置获得图像。当观察对象的光强衰减到一定程度，达到单光子水平，就变成了离散脉冲信号。单光子是一种极弱光，被认为是光不可分割的最小能量单位，是可以探测的极限。单光子探测技术应用于生物自发光、医疗诊断、非破坏性物质分析、天文观测、光谱测量、量子光学等领域，并在其中扮演着重要角色。研究极弱光成像探测技术在这些领域应用的发展非常有意义。

光子计数成像就是一种极弱光探测技术，通常它通过记录成像位置的光子计数以及探测到光子的概率，在数据处理端进行累计和融合获得一幅图像，其核心是面元探测器，面元探测器规模(阵列大小)、灵敏度范围、以及响应波段直接影响能否获得单光子水平的图象获取质量。但是，用于单光子探测水平的面元探测器不但价格昂贵，只能够在少数波段可以实现，且面元探测器灵敏度低，存在技术不甚成熟与极弱光对象二维成像的强烈需求之间的矛盾。

压缩传感理论(CS理论)是由E.J.Candes等人提出的，它打破了传统的线性采样模式，表明可压缩信号的少量线性随机投影中包含足够的信息来重建原信号。

CS理论包括两部分：压缩采样和稀疏重建。

压缩采样是被测信号由高维向低维映射的过程。假设x∈Rⁿ是被测数据，y∈R^k是观测数据，Φ∈R^k×n是随机投影矩阵(k＜＜n)，e∈R^k是测量噪声，那么，压缩采样过程可以描述为(1)式：

y＝Φx+e              (1)

如果x是变换域稀疏的，即θ＝Ψx，Ψ是稀疏变换矩阵，那么(1)式变化为(2)式：

y＝ΦΨ^-1θ+e         (2)

随机投影矩阵Φ，也叫测量矩阵，需要满足RIP(Restricted Isometry Property)：

$(1-{\mathrm{\δ}}_s){\left|\right|x\left|\right|}_2^2\≤{\left|\right|\mathrm{Ax}\left|\right|}_2^2\≤(1+{\mathrm{\δ}}_s){\left|\right|x\left|\right|}_2^2$

其中，δ_s定义为使所有s稀疏的向量x都满足此不等式的最小常数，并且δ_s＜1。

另外，Φ与Ψ越不相关，采样所需的测量数k越小，所以，一般情况下Φ设计为随机矩阵。

稀疏重建实际上是在已知观测数据y和测量矩阵Φ的条件下求解(1)式中的x，这是一个ill-posed问题，一般用最优化方法求解，可描述为(3)式：

$\underset{x\∈R^n}{\min }(\frac{1}{2}{\left|\right|y-\mathrm{\Φx}\left|\right|}_2^2+\mathrm{\τ}{\left|\right|x\left|\right|}_1)---\left(3\right)$

如果x是变换域稀疏的，对应于(2)式的重建问题可以描述为(4)式：

$\underset{x\∈R^n}{\min }(\frac{1}{2}{\left|\right|y-\mathrm{\Φx}\left|\right|}_2^2+\mathrm{\τ}{\left|\right|\mathrm{\Ψx}\left|\right|}_1)---\left(4\right)$

(3)式和(4)式中，第一项是最小二乘约束，记为f(x)；第二项是对x稀疏度的一种约束；两项之和是目标函数，记为

DLP技术是美国德州仪器公司(TI)提出的一项技术，它与数字视频或图形信号结合起来，其微镜和透镜系统可以将数字图像反射到屏幕或其它表面，其核心是DLP芯片——数字微镜器件(Digital Micro-mirror Device，简称DMD控制系统)，这目前是世界上最精密的光开关。它包含一个多达200万个安装在铰链上的微镜的矩阵，每个微镜的大小小于人的头发丝的五分之一，每一个微镜都可以在一定角度范围内(通常-12°和+12°)摆动，如果把这两种状态记为0和1，那么，用脉宽调制波(PWM)来驱动微镜，使其在0和1之间高速抖动，就可以实现中间状态。DMD控制系统及其相关的精密的电子元件就是所谓的DLP技术，该技术有成熟产品，广泛用于投影仪器等产品中。

这种“先采样，后重建”的思想使得将二维信号转换为随时间分布的一维信号，并且由单个探测器采样成为可能。而点探测器无论在探测灵敏度、波长范围具有更宽的选择范围，成本优势明显，利用点探测器实现单光子记数成像成为未来单光子水平成像的重要发展趋势。

发明内容

本发明的目的在于，为解决目前面元探测器灵敏度低、技术不甚成熟与极弱光对象二维成像的强烈需求之间的矛盾，从而提供一种单光子计数成像系统及其方法，以压缩传感(Compressive Sensing，CS)理论为基础，采用DLP技术将图象信号随机变化变成随机的光强信号，再以单光子计数器为探测元件获得记数信号，实现了用点探测器对极弱光对象进行二维成像。

为实现上述目的，本发明提供了一种单光子计数成像系统，其特征在于，该单光子计数成像系统采用压缩传感理论和DLP技术，并以单光子计数器为探测元件，实现了单光子级别的极弱光对象的二维成像，所述的单光子计数成像系统包括：滤光片、第一透镜1、DMD控制系统、第二透镜2、单光子计数器和数据处理单元；其中，

DMD结合第一透镜1和第二透镜2，用于将二维图像数据转化为一维数据序列完成被测信号的压缩采样，极弱光通过滤光片滤除杂光，经第一透镜1在DMD控制系统处成像，并由DMD控制系统控制光子被反射到第二透镜2的概率，经过第二透镜2控制光子聚焦；

数据处理单元结合单光子计数器完成稀疏重建，数据处理单元根据单光子计数器在一定时间内对光子进行计数折算成探测到光子数的概率作为测量值，和DMD控制系统上的测量矩阵经过最优化算法重建光子密度图像，解算出二维图像。

作为上述技术方案的一种改进，所述的第二透镜2聚焦后至单光子计数器的光路上还设有光衰减器，用于将光衰减到单光子探测器的工作范围。该光衰减器的设计是为了防止被测光子密度过大和单光子计数器的门控时间过长引起的饱和。

为实现上述的另一发明目的，本发明还提供了一种单光子计数成像方法，该方法采用了压缩传感理论和DLP技术，并以单光子计数器为探测元件，实现了单光子级别的极弱光对象的二维成像，其步骤包括：

1)压缩采样的步骤；

所述的压缩采样是由DMD控制系统结合第一透镜1和第二透镜2，用于将二维图像数据转化为一维数据序列完成被测信号的压缩采样，

极弱光通过滤光片滤除杂光，经第一透镜1在DMD控制系统处成像，并由DMD控制系统控制光子被反射到第二透镜2的概率，经过第二透镜2控制光子聚焦；

2)稀疏重建的步骤；

所述的稀疏重建是由数据处理单元根据单光子计数器在一定时间内对光子进行计数折算成探测到光子数的概率作为测量值，和DMD控制系统上的测量矩阵经过最优化算法重建光子密度图像，解算出二维图像。

作为上述技术方案的一种改进，所述的方法包括步骤如下：

所述的压缩采样，是被测信号由高维向低维映射的过程：

11)假设x∈Rⁿ是被测数据，y∈R^k是观测数据，Φ∈R^k×n是随机投影矩阵(k＜＜n)，e∈R^k是测量噪声，那么，压缩采样的过程可以描述为(1)式：

y＝Φx+e             (1)

如果x是变换域稀疏的，即θ＝Ψx，Ψ是稀疏变换矩阵，那么(1)式变化为(2)式：

y＝ΦΨ^-1θ+e        (2)

式中，Ψ为小波变换矩阵，Φ为Gaussian随机矩阵；

12)假设测量数为k，二维图像的像素个数为n，则(1)式中的测量矩阵则为Φ＝{Φ₁，…，Φ_i，…，Φ_k}，Φ_i是Φ的第i行，把的二维图像的列首尾相连，化成n×1的一维列向量，对应(1)式中的x，其中的每一个元素代表相应位置处的光子密度；DMD控制系统具有同样的分辨率，它的列首尾相连，化成1×n的一维行向量，对应测量矩阵Φ中的一行，其中的每个元素代表相应位置处光子透射到第二透镜2的概率；

13)假设测量周期为T，在这段时间内，DMD控制系统保持不变，单光子计数器探测到的光子数为N，那么N/T就相当于光子密度图像与DMD控制系统上的随机数阵列的内积值，对应于(1)式中观察向量y的一个元素式中，Φ_i，j、x_j分别是Φ_i和x的第j个元素；根据测量矩阵改变DMD控制系统，重复k次测量，就可以得到整个观测数据y；

14)所述的稀疏重建是在已知观测数据y和测量矩阵Φ的条件下求解(1)式中的x，一般用最优化方法求解，可描述为(3)式：

$\underset{x\∈R^n}{\min }(\frac{1}{2}{\left|\right|y-\mathrm{\Φx}\left|\right|}_2^2+\mathrm{\τ}{\left|\right|x\left|\right|}_1)---\left(3\right)$

如果x是变换域稀疏的，对应于(2)式的重建问题可以描述为(4)式：

$\underset{x\∈R^n}{\min }(\frac{1}{2}{\left|\right|y-\mathrm{\Φx}\left|\right|}_2^2+\mathrm{\τ}{\left|\right|\mathrm{\Ψx}\left|\right|}_1)---\left(4\right)$

(3)式和(4)式中，第一项是最小二乘约束，记为f(x)；第二项是对x稀疏度的一种约束；两项之和是目标函数，记为

作为上述技术方案的进一步的改进，所述的最优化方法是采用IWT(IterativeWavelets Thresholding)算法，把本次迭代的估计值做DWT变换，对变换系数做阈值处理，再做DWT反变换得到下一次迭代的估计值；

如果阈值处理函数描述为S(u，v)＝sign(u)max{|u|-v，0}，那么该算法可以描述为：

$x^{t+1}={\mathrm{\Ψ}}^{-1}\·S(\mathrm{\Ψ}(x^t-\frac{1}{\mathrm{\α}}\▿f\left(x^t\right)),\frac{\mathrm{\τ}}{\mathrm{\α}})$

其中，α＝||Φ^TΦ||₂。

作为上述技术方案的另一种改进，所述的方法在所述的第二透镜2聚焦后至单光子计数器的光路上，还通过光衰减器将光衰减到单光子探测器的工作范围。

本发明的优点在于：本发明采用压缩传感理论(Compressive Sensing，简称CS)结合数字光处理(Digital Light Processing，简称DLP)技术的方案来解决利用点探测器实现高探测灵敏度的成像问题，其灵敏度可以达到单光子水平，分辨率与DMD控制系统直接相关，而DMD控制系统目前可以达到很高的分辨率。本发明可广泛应用于生物自发光检测、医疗诊断、非破坏性物质分析、天文观测、国防军事、光谱测量、量子电子学等领域。

本发明以压缩传感(Compressive Sensing，CS)理论为基础，以单光子点探测器为探测元件，用单光子计数器实现极弱光的二维成像，结构简单，灵敏度可以达到单光子水平，分辨率与DMD控制系统直接相关，而DMD控制系统目前可以达到很高的分辨率，解决了目前该领域中焦平面传感器灵敏度低、阵列规模小、探测波长范围相对单一与极弱光对象二维成像的强烈需求之间的矛盾。

附图说明

图1是本发明的单光子计数成像系统的结构示意图。

图2是对本发明的模拟实验结果，图2(a)是原始光子密度图像；图2(b)是一次测量中DMD控制系统上的随机矩阵，黑点代表0，白点代表1，灰点代表中间值；图2(c)是IWT算法的重建图像，与原始图像的相关系数Cov＝0.9783，信噪比PSNR＝23.95dB；图2(d)是IWT算法的残差图像。

具体实施方式

以下结合附图对本发明作进一步的详细说明。

如图1所示，观察对象发出的极弱光被滤光片滤除杂光，经第一透镜1在DMD控制系统处成像，由DMD控制系统控制光子被反射到第二透镜2的概率，经过第二透镜2的控制光子会聚点，光衰减器作用是当光过强时将光衰减到单光子探测器的工作范围，由单光子计数器在一定时间内对光子进行计数，该数值可以折算成探测到光子数的概率作为测量值。最后由数据处理单元根据测量值和DMD控制系统上的测量矩阵经过最优化算法重建光子密度图像。其中，光衰减器的设计是为了防止被测光子密度过大和单光子计数器的门控时间过长引起的饱和。

为了便于数学上的理解，假设测量数为k，(1)式中的测量矩阵写为Φ＝{Φ₁，…，Φ_i，…，Φ_k}，Φ_i是Φ的第i行。我们把的二维图像的列首尾相连，化成n×1的一维列向量，对应(1)式中的x，其中的每一个元素代表相应位置处的光子密度；DMD控制系统具有同样的分辨率，它的列首尾相连，化成1×n的一维行向量，对应测量矩阵Φ中的一行，其中的每个元素代表相应位置处光子透射到第二透镜2的概率。假设测量周期为T，在这段时间内，DMD控制系统保持不变，单光子计数器探测到的光子数为N，那么N/T就相当于光子密度图像与DMD控制系统上的随机数阵列的内积值，对应于(1)式中观察向量y的一个元素(Φ_i，j、x_j分别是Φ_i和x的第j个元素)。根据测量矩阵改变DMD控制系统，重复k次测量，就可以得到整个观测数据y，在物理上实现(1)式的过程。

根据光子学的知识，在一个元面积dA内，任意时刻在r点观察到一个光子的概率p(r)dA正比于该处光强。我们用生物芯片的灰度图像模拟光子密度图像。

生物芯片是典型的极弱光源，目前主要通过荧光标记的方法使其便于观察。实际上生物都有自发光的特性，并且自发光光谱包含很多重要的信息。采用光子计数成像技术就可以直接观测。实验中选择分辨率为64×64的生物芯片图像，灰度为256级，最高灰度级对应光子数为4.0×10²s^-1。在假设不知道原图像的情况下，采用Gaussian矩阵进行压缩采样，IWT稀疏重建算法进行图像重建，得到图2所示的结果，其中，(a)是原始光子密度图像；(b)是一次测量中DMD控制系统上的随机矩阵，黑点代表0，白点代表1，灰点代表中间值；(c)是IWT算法的重建图像；(d)是IWT算法的残差图像。

图2(c)所示的重建图像与原始图像图2(a)的相关系数Cov＝0.9783，信噪比PSNR＝23.95dB，满足对生物芯片图像的观测需求。实验结果表明，本发明硬件可以实现压缩采样，软件算法可以实现图像重建。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (2)

1.一种单光子计数成像方法，该方法采用了压缩传感理论和DLP技术，并以单光子计数器为探测元件，实现了单光子级别的极弱光对象的二维成像，其步骤包括：1)压缩采样的步骤；

所述的压缩采样是由DMD控制系统结合第一透镜(1)和第二透镜(2)，用于将二维图像数据转化为一维数据序列完成被测信号的压缩采样，极弱光通过滤光片滤除杂光，经第一透镜(1)在DMD控制系统处成像，并由DMD控制系统控制光子被反射到第二透镜(2)的概率，经过第二透镜(2)控制光子聚焦；

2)稀疏重建的步骤；

所述的稀疏重建是由数据处理单元根据单光子计数器在一定时间内对光子进行计数折算成探测到光子数的概率作为测量值，和DMD控制系统上的测量矩阵经过最优化算法重建光子密度图像，解算出二维图像；所述的方法包括步骤如下：

所述的压缩采样，是被测信号由高维向低维映射的过程：

11)假设x∈Rⁿ是被测数据，y∈R^k是观测数据，Φ∈R^k×n是随机投影矩阵(k＜＜n)，e∈R^k是测量噪声，那么，压缩采样的过程可以描述为(1)式：

y＝Φx+e          (1)

如果x是变换域Ψ下展开的系数是稀疏的，即x＝Ψθ，Ψ是稀疏变换矩阵，θ为稀疏系数，则(1)式变化为(2)式：

y＝ΦΨθ+e             (2)

式中，Ψ为小波变换矩阵，Φ为Gaussian随机矩阵；

12)假设测量数为k，二维图像的像素个数为n，则(1)式中的测量矩阵则为Φ＝{Φ₁,…,Φ_i,…,Φ_k}，Φ_i是Φ的第i行的行向量，Φ₁为Φ的第1行的行向量，Φ_k为Φ的第k行的行向量；把的二维图像的列首尾相连，化成n×1的一维列向量，对应(1)式中的x，其中的每一个元素代表相应位置处的光子密度；DMD控制系统具有同样的分辨率，它的列首尾相连，化成1×n的一维行向量，对应测量矩阵Φ中的一行，其中的每个元素代表相应位置处光子透射到第二透镜(2)的概率；

13)假设测量周期为T，在这段时间内，DMD控制系统保持不变，单光子计数器探测到的光子数为N，那么N/T就相当于光子密度图像与DMD控制系统上的随机数阵列的内积值，对应于(1)式中观察向量y的一个元素式中，Φ_i,j、x_j分别是Φ_i和x的第j个元素；根据测量矩阵改变DMD控制系统，重复k次测量，就可以得到整个观测数据y；

14)所述的稀疏重建是在已知观测数据y和测量矩阵Φ的条件下求解(1)式中的x，一般用最优化方法求解，可描述为(3)式：

$\underset{x\∈R^n}{\min }(\frac{1}{2}{\left|\right|y-\mathrm{\Φx}\left|\right|}_2^2+\mathrm{\τ}{\left|\right|x\left|\right|}_1)---\left(3\right)$

如果x是在变换域Ψ下是稀疏的，对应于(2)式的重建问题可以描述为(4)式：

$\underset{x\∈R^n}{\min }(\frac{1}{2}{\left|\right|y-\mathrm{\Φ\Ψ\θ}\left|\right|}_2^2+\mathrm{\τ}{\left|\right|\mathrm{\θ}\left|\right|}_1)---\left(4\right)$

(3)式和(4)式中，第一项是最小二乘约束，记为f(x)；第二项是对x稀疏度的一种约束；两项之和是目标函数，记为

所述的最优化方法是采用IWT算法，把本次迭代的估计值做DWT变换，对变换系数做阈值处理，再做DWT反变换得到下一次迭代的估计值；

如果阈值处理函数描述为S(u,v)＝sign(u)max{|u|-v,0}，那么，该算法可以描述为：

$x^{t+1}={\mathrm{\Ψ}}^{-1}\·S(\mathrm{\Ψ}(x^t-\frac{1}{\mathrm{\α}}\▿f\left(x^t\right)),\frac{\mathrm{\τ}}{\mathrm{\α}})$

其中，α＝||Φ^TΦ||₂。

2.根据权利要求1所述的单光子计数成像方法，其特征在于，所述的方法在所述的第二透镜(2)聚焦后至单光子计数器的光路上，还通过光衰减器将光衰减到单光子探测器的工作范围。