CN111896125B - 一种单光子计数成像的偏振去噪方法 - Google Patents

Abstract

本发明提供了一种单光子计数成像的偏振去噪方法，包括：通过斯托克斯(stokes)矢量和穆勒矩阵计算出两个方向的探测器探测得到的总强度，并定义一个新的反射率指标α_N；通过计算反射率α和线性去极化参数β的负对数似然函数，并添加非局部先验稀疏表示约束，构建单光子计数成像偏振去噪模型；通过交替最小化方法解决优化问题，得到并输出反射率和去极化参数图像。本发明方法可以在高背景噪声环境下，利用信号光与背景光之间的偏振差异，从根源上提高了信噪比，同时结合非局部先验稀疏表示约束，不仅能获得高质量的反映目标特征的反射率图像，还能得到目标本身的线性去极化参数图像，获取目标的偏振信息。

Description

一种单光子计数成像的偏振去噪方法

技术领域

本发明属于单光子计数成像领域，涉及一种单光子计数成像的去噪方法，适用于绝大多数主动成像系统如雷达等应用场景。

背景技术

单光子计数成像技术是一项通过对在目标探测所在场景中反射回来的光子进行计数，从而获得目标信息，进而重建目标图像的技术，在遥感以及空间探测等领域内都吸引了广泛的关注。于传统成像技术不同的是，单光子计数成像技术里使用的光电探测器是具有单光子探测能力的盖革模式雪崩光电二极管(GM-APD)。因此，它能在极低光照水平下探测到物体。然而，在低光照条件下，噪声会极大的影响图像重建质量。传统的去噪方法如非局部均值法[1]，虽然能直接的应用于噪声的去除，但它们忽略了噪声的特性。在低光照条件下，泊松分布往往被用来描述反射光子的计数过程。由此，2011年，Makitalo等人[2]提出采用基于方差稳定变换(VST)的高斯过程近似泊松分布，并利用块匹配和3D滤波(BM3D)降低噪声。然而，在低光照环境下，当返回的光子数量很少时，这种近似并不准确[3]。最近，最小化负泊松对数似然项[4]，并添加稀疏项作为先验约束[5]用以减少噪声成为单光子计数成像的一种常见手段。在单光子计数成像技术中，完整的噪声由暗计数和背景噪声组成。暗计数噪声是由探测器的暗计数以较小的计数率引起的，而背景噪声是由环境光产生的，在某些应用中，例如白天使用单光子激光，它可能是个相对值。在低信噪比(SNR)情况下，传统的去噪方法无法轻易的降低环境噪声的影响。

偏振是光波的一个重要特性。当光作用在任何物体上，发生反射、投射及散射，都有可能产生特定的偏振态。通过记录光与物质发生作用后的偏振特性图像，可以提供场景中的物体更多的特征，对于滤除背景或其它无效信息，增强图像质量都有很重要的意义。因此，偏振成像技术在目标探测与识别领域有着非常广泛的应用。我们注意到，大多数的主动成像技术如激光雷达等都带有偏振信息，而环境光是自然光不带有偏振信息。如何利用这种偏振差异，构建有效的重建算法，在高环境光情况下，实现单光子计数成像去噪，从而获得高质量的反映目标特征的图像，成为该领域研究的重点问题。

发明内容

本发明针对现有技术的不足，利用偏振信息，提供一种单光子计数成像的去噪方法，从而得到在高环境光状态下目标的发射率重建图像，同时获取目标本身的线性去极化参数。

本发明所采用的技术方案是一种单光子计数成像的偏振去噪方法，该方法包括以下步骤：

步骤1，通过斯托克斯(stokes)矢量和穆勒矩阵计算出两个方向的探测器探测得到的总强度，并定义一个新的反射率指标α_N；

步骤2，通过计算反射率α和线性去极化参数β的负对数似然函数，并添加非局部先验稀疏表示约束，构建单光子计数成像偏振去噪模型；

步骤3，通过交替最小化方法解决优化问题，得到并输出反射率和去极化参数图像。

进一步地，步骤1中，已知偏振态的用以照射目标的周期性脉冲激光表示为：

S_laser＝I₀[1 a₁ b₁ c₁]^T

其中，I₀为脉冲激光的光强度，a₁，b₁，c₁代表归一化后的stokes矢量，同时满足a₁ ²+b₁ ²+c₁ ²＝1。经目标后的反射光的偏振态可表示为：

S_re＝M_tar×S_laser

其中，M_tar为目标物体表面的穆勒矩阵。对于非双折射材料的光反射，穆勒矩阵可以近似为

其中，α为目标的反射率，β和χ分别为目标物体的线性去极化参数和圆形去极化参数。在接收系统中，存在一个相对于x方向成角度θ的偏振片，该偏振片的穆勒矩阵可表示为

在真实场景中，背景光的影响不可忽略，我们将其光强度表示为I_B。由于背景光为自然光，所以通过线性偏振片后，其光强度会变为原来的一半，因此，到达探测器的总光强I为

上式表明了，偏振片的角度对于信噪比(Signal Noise Ratio，SNR)的提升起到重要作用。使SNR最高的角度θ_P可以由下式求出

得到

这时在该角度的总光强I_P为

我们注意到，噪声的强度减少了一半，而信号光光强度的减少与目标的线性去极化参数β和激光的极化参数a₁和b₁相关。由此，我们定义一个新的反射率指标(NRI)α_N，拥有最高的SNR，表示为

α_N与反射率α以及去极化参数β高度相关联，因此，α_N有助于α和β的求解。

为求解出α和β，我们转动偏振片至角度θ_M，这时到达探测器的总光强I_M为

我们选择θ_M使a₁cos2θ_M+b₁sin2θ_M的值尽可能的减少，从而使得I_P和I_M之间的差异尽可能的大。为方便，这里我们确定θ_M＝θ_P+90°。

进一步的，步骤2的具体实现包括了以下子步骤：

步骤2.1，由单光子计数模块记录每个像素(i，j)上收集到的光子数X(i，j)。根据分别在θ_P和θ_M角度上得到的光子计数X_θP(i，j)和X_θM(i，j)，计算关于反射率α和线性去极化参数β的负对数似然函数：

其中

η为探测器的光电探测效率，D表示每个脉冲重复周期内探测器的暗计数，B和S分别代表每个脉冲重复周期内背景计数和总信号计数，N为激光照射到目标的每个像素的脉冲数。

m×n为单光子计数图像的分辨率，P_θP和P_θM分别为在角度为θ_P和θ_M的单个照明脉冲中检测到的光子数量X(i，j)的概率质量函数。

步骤2.2，由于NRI具有较高的SNR属性，同时与α以及β高度相关，由此，我们提出了一种非局部先验稀疏表示约束。首先，我们将图像分成重叠的小块进行处理，以反射率为例，给定一个噪声较大的色块α^*，基于块匹配算法从NRI中探究出K个相似的色块，生成一个具有较高信噪比的色块组，作为先验知识。

步骤2.3，结合上述约束，最终构建单光子计数成像偏振去噪模型：

其中后两项皆为先验约束，D为字典，利用K-SVD(一种经典的字典学习方法)从色块组α^*中学习获得字典，λ_α和λ_β为平衡不同项之间的参数，w_α为反射率图像α的稀疏表示，w_β为线性去极化参数图像β的稀疏表示。

进一步的，步骤3的具体实现包括了以下子步骤：

步骤3.1，我们通过交替最小化方法将该复杂的优化问题分解为以下子问题：

其中，上标t表示外部迭代次数。

步骤3.2，对于反射率α的子问题，其方程可以转化为以下形式：

其中

上标k表示内部迭代次数。该方程可以根据先验约束直接采用不同的现有方法来求解。

是

对于α_N的梯度，表示为

τ^k是一个参数，可以根据两点步长梯度法Barzilai-Borwein(BB)在[τ_min，τ_max]范围内进行选择，即

其中，[τ_min，τ_max]是τ^k的最小值和最大值，

是

对于α_N的二阶梯度，表示为

最后求得α_N。

步骤3.3，对于线性去极化参数β子问题，其方程可以转化成以下形式：

其中

上标k表示内部迭代次数。该方程可以根据先验约束直接采用不同的现有方法来求解。

是

对于β的梯度，表示为

τ^k是一个参数，可以根据改进的Barzilai-Borwein(BB)方法在[τ_min，τ_max]范围内进行选择，即

其中，

是

对于β的二阶梯度，表示为

最后求得β，最终根据α_N表达式求得α。

由上可见，本发明所提供的一种单光子计数成像的偏振去噪方法，能够在高环境噪声背景下实现对单光子计数图像的去噪，利用信号光与背景光之间的偏振差异，从根源上提高信噪比，同时结合非局部先验稀疏表示约束，不仅能获得高质量的反映目标特征的反射率图像，还能得到目标本身的线性去极化参数图像，获取目标的偏振信息。

附图说明

图1是本发明总体流程图。

图2是本发明非局部先验稀疏表示约束示意图。

图3是本发明的算法流程图。

图4是本发明实施例中仿真得到的传统的不带偏振信息的单光子计数图像与带偏振信息的两个角度的单光子计数图像的视觉比较。

图5是本发明提出的方法与其他对比方法重建的反射率图像的视觉比较。

图6是本发明提出的方法得到的线性去极化参数图像。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下参照附图并举实施例，对本发明进一步详细说明。

如图1所示，本实施例提供了一种在高背景噪声状态下实现单光子计数成像的偏振去噪方法，具体包括如下几个步骤：

步骤1：通过斯托克斯(stokes)矢量和穆勒矩阵计算出两个方向的探测器探测得到的总强度，并定义一个新的反射率指标α_N；

传统的单光子计数成像系统，一般由激光光源、光学收发扫描系统、单光子探测器、时间相关单光子技术模块以及数据处理与系统控制模块这几个部分组成，通过光栅扫描的方式使用周期性脉冲激光照亮场景，由单光子探测器探测反射信号，通过计数器进入数据处理模块，重建目标图像。而我们的单光子计数成像系统仅仅在传统的系统中添加一个偏振片(探测器前)，用以获得偏振信息，其他原理相同。

已知偏振态的用以照射目标的周期性脉冲激光表示为：

S_laser＝I₀[1 a b c]^T

其中，I₀为脉冲激光的光强度，a，b，c代表归一化后的stokes矢量，同时满足a²+b²+c²＝1。经目标后的反射光的偏振态可表示为：

S_re＝M_tar×S_laser

其中，M_tar为目标物体表面的穆勒矩阵。如[6，7]所述，对于非双折射材料的光反射，穆勒矩阵可以近似为

其中，α为目标的反射率，β和χ分别为目标物体的线性去极化参数和圆形去极化参数。在接收系统中，存在一个相对于x方向成角度θ的偏振片，该偏振片的穆勒矩阵可表示为

在真实场景中，背景光的影响不可忽略，我们将其光强度表示为I_B。由于背景光为自然光，所以通过线性偏振片后，其光强度会变为原来的一半，因此，到达探测器的总光强I为

上式表明了，偏振片的角度对于SNR的提升起到重要作用。使SNR最高的角度θ_P可以由下式求出

得到

这时在该角度的总光强I_P为

我们注意到，噪声的强度减少了一半，而信号光光强度的减少与目标的线性去极化参数β和激光的极化参数a和b相关。由此，我们定义一个新的反射率指标(NRI)α_N，拥有最高的SNR，表示为

α_N与反射率α以及去极化参数β高度相关联，因此，α_N有助于α和β的求解。

为求解出α和β，我们转动偏振片至角度θ_M，这时到达探测器的总光强I_M为

我们选择θ_M使acos2θ_M+bsin2θ_M的值尽可能的减少，从而使得I_P和I_M之间的差异尽可能的大。为方便，这里我们确定θ_M＝θ_P+90°。

步骤2：通过计算反射率α和线性去极化参数β的负对数似然函数，并添加非局部先验稀疏表示约束，构建单光子计数成像偏振去噪模型；

步骤2.1：通过计算反射率α和线性去极化参数β的负对数似然函数。由单光子计数模块记录每个像素(i，j)上收集到的光子数X(i，j)。非均匀泊松过程可用于描述由单光子计数模块计数的光子的精确模型，于是，零检测的概率P₀可以表示为

其中

η为探测器的光电探测效率，D表示每个脉冲重复周期内探测器的暗计数，B和S分别代表每个脉冲重复周期内背景计数和总信号计数。D，B，S皆为恒定值。我们用一系列N个脉冲照射每个像素(i，j)。对于单光子探测器来说，只有零探测或一个探测这两种可能。因此，在角度为θ的单个照明脉冲中检测到的光子数量X(i，j)的概率质量函数可以表示为

由此，根据分别在θ_P和θ_M角度上得到的光子计数X_θP(i，j)和X_θM(i，j)，计算关于反射率α和线性去极化参数β的负对数似然函数：

其中

m×n为单光子计数图像的分辨率。

步骤2.2：添加非局部先验稀疏表示约束。由于NRI具有较高的SNR属性，同时与α以及β高度相关，由此，我们提出了一种非局部先验稀疏表示约束。如图2所示，首先，我们将图像分成重叠的小块进行处理，以反射率为例，给定一个噪声较大的色块α，基于块匹配算法从NRI中探究出K个相似的色块，生成一个具有较高信噪比的色块组，作为先验知识。在该实施例中β的先验约束与α的先验约束相同。

步骤2.3：构建单光子计数成像偏振去噪模型。结合上述约束，最终构建单光子计数成像偏振去噪模型：

其中后两项皆为先验约束，D为字典，λ_α和λ_β为平衡不同项之间的参数。

步骤3：通过交替最小化方法解决优化问题，得到并输出反射率和去极化参数图像。

步骤3.1：分解优化问题。我们通过交替最小化方法将该复杂的优化问题分解为以下子问题：

其中，上标t表示外部迭代次数。

步骤3.2：交替求解。以反射率α的子问题为例，其方程可以转化为以下形式：

其中

上标k表示内部迭代次数。λ_αN为正则化参数。该方程可以根据先验约束直接采用不同的现有方法来求解。

是

对于α_N的梯度，表示为

τ^k是一个参数，可以根据改进的Barzilai-Borwein(BB)方法在[τ_min，τ_max]范围内进行选择，即

其中，

是

对于α_N的二阶梯度，表示为

以此类推，β子问题求解方式与上述相同，最后求得α_N与β，最终根据α_N表达式求得α，具体算法流程如图3所示。

具体的，这里使用的测试数据是根据成像过程仿真得出。在仿真过程中，设置照射激光的偏振态参数a＝1，b＝0，c＝0，每个像素照射的脉冲数N为40，每个脉冲重复周期内的背景计数B为0.1，光子计数图像的像素值为200×200。仿真结果如图4，得到传统的不带偏振信息的单光子计数图像X₀与带偏振信息的两个角度的单光子计数图像X_θP和X_θM。如图4所示，图像X_θP的图像质量明显优于图像X_θM。这时因为在背景光(自然光)通过偏振片之后，光强度变为原始强度的一半。由于信号光的偏振特性，光强度的减小要小得多，因此SNR相应增加。由于θ_M＝θ_P+90°，以使两个角度方向之间的光强度差异尽可能大，因此在图像X_θM中几乎没有物体信息。使用主成分分析与局部像素分组(LPG-PCA)和BM3D+VST方法和我们的方法进行比较。在该实施例中，设置参数内外迭代次数皆为5，正则化参数λ_αN＝0.01，λ_β＝0.01。有关各算法重建出的反射率图像视觉质量对比如图5。为了对重建结果进行定量的评价，我们选用峰值信噪比(PSNR)以及结构相似性(SSIM)作为评价指标，并分别改变变量N，B以及激光的偏振态参数a，b(c²＝1-a²-b²)，重建反射率结果如表1所示。

表1不同重建方法的客观比较结果

可以看到，我们的方法利用信号光与背景光之间的偏振差异，从根源上提高了信噪比，同时结合非局部先验稀疏表示约束，可以很好的重建目标的反射率图像，无论在视觉质量还是评价指标上都占据优势，表现更好。与此同时，如图6所示，还能得到目标本身的线性去极化参数图像，获取目标的偏振信息。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，本说明书中所描述的具体实施例仅仅是对本发明作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

参考文献

[1]L.Lu,W.Jin,and X.Wang,“Non-local means image denoising with a softthreshold,”IEEE Signal Process.Lett.22,833–837(2015).

[2]M.Makitalo and A.Foi,“Optimal inversion of the anscombetransformation in low-count poisson image denoising,”IEEE Transactions onImage Process.20,99–109(2011).

[3]E.D.Kolaczyk,“Wavelet shrinkage estimation of certain poissonintensity signals using corrected thresholds,”Stat.Sinica9,119–135(1999).

[4]D.Shin,F.Xu,D.V enkatraman,R.Lussana,F.Villa,F.Zappa,V.Goyal,F.Wong,and J.Shapiro,“Photon-efficient imaging with a single-photon camera,”Nat.Commun.7,12046(2016).

[5]D.Shin,A.Kirmani,V.Goyal,and J.Shapiro,“Photon-efficientcomputational 3-d and reflectivity imaging with single-photon detectors,”IEEETransactions on Comput.Imaging 1(2015).

Claims (5)

1.一种单光子计数成像的偏振去噪方法，其特征在于，包含以下步骤：

步骤1，通过斯托克斯stokes矢量和穆勒矩阵计算出两个方向的探测器探测得到的总强度，并定义一个新的反射率NRI指标α_N；

其中，α为反射率，β为线性去极化参数，a₁，b₁代表归一化后的stokes矢量中的两个参数；

步骤2，通过计算反射率α和线性去极化参数β的负对数似然函数，并添加非局部先验稀疏表示约束，构建单光子计数成像偏振去噪模型；

步骤3，通过交替最小化方法解决优化问题，得到并输出反射率和去极化参数图像。

2.如权利要求1所述的一种单光子计数成像的偏振去噪方法，其特征在：步骤1中两个角度满足，θ_P为任意角度，θ_M＝θ_P+90°。

3.如权利要求I所述的一种单光子计数成像的偏振去噪方法，其特征在：步骤1中计算两个方向的探测器探测得到的总强度的具体实现方式如下，

已知偏振态的用以照射目标的周期性脉冲激光表示为：

S_laser＝I₀[1 a₁ b₁ c₁]^T

其中，I₀为脉冲激光的光强度，a₁，b₁，c₁代表归一化后的stokes矢量，同时满足a₁ ²+b₁ ²+c₁ ²＝1，经目标后的反射光的偏振态表示为：

S_re＝M_tar×S_laser

其中，M_tar为目标物体表面的穆勒矩阵；对于非双折射材料的光反射，穆勒矩阵近似为

其中，α为目标的反射率，β和χ分别为目标物体的线性去极化参数和圆形去极化参数；在接收系统中，存在一个相对于x方向成角度θ的偏振片，该偏振片的穆勒矩阵表示为

在真实场景中，将背景光的光强度表示为I_B，由于背景光为自然光，所以通过线性偏振片后，其光强度会变为原来的一半，因此，到达探测器的总光强I为

上式表明了，偏振片的角度对于SNR的提升起到重要作用，使SNR最高的角度θ_P由下式求出

s.t.a₁ ²+b₁ ²+c₁ ²＝1.

得到

这时在该角度的总光强I_P为

为求解出α和β，转动偏振片至角度θ_M，θ_M＝θ_P+90°，这时到达探测器的总光强I_M为，

4.根据权利要求2或3所述的一种单光子计数成像的偏振去噪方法，其特征在于：步骤2的具体实现包括以下子步骤：

步骤2.1，由单光子计数模块记录每个像素(i，j)上收集到的光子数X(i，j)，根据分别在θ_P和θ_M角度上得到的光子计数X_θP(i，j)和X_θM(i，j)，计算关于反射率α和线性去极化参数β的负对数似然函数：

其中

η为探测器的光电探测效率，D表示每个脉冲重复周期内探测器的暗计数，B和S分别代表每个脉冲重复周期内背景计数和总信号计数，N为激光照射到目标的每个像素的脉冲数，

m×n为单光子计数图像的分辨率，P_θP和P_θM分别为在角度为θ_P和θ_M的单个照明脉冲中检测到的光子数量X(i，j)的概率质量函数；

步骤2.2，将图像分成重叠的小块进行处理，给定一个色块α^*，基于块匹配算法从NRI中探究出K个相似的色块，生成一个具有一定信噪比的色块组，作为先验知识；

步骤2.3，结合上述约束，最终构建单光子计数成像偏振去噪模型：

s.t.α＝D×w_α，β＝D×w_β

其中后两项皆为先验约束，D为字典，λ_α和λ_β为平衡不同项之间的参数，w_α为反射率图像α的稀疏表示，w_β为线性去极化参数图像β的稀疏表示。

5.根据权利要求4所述的一种单光子计数成像的偏振去噪方法，其特征在于：步骤3的具体实现包括以下子步骤：

步骤3.1，通过交替最小化方法将复杂的优化问题分解为以下子问题：

其中，上标t表示外部迭代次数；

步骤3.2，对于反射率α的子问题，其方程可以转化为以下形式：

其中

上标k表示内部迭代次数；

是

对于α_N的梯度，表示为

τ^k是一个参数，根据Barzilai-Borwein方法在[τ_min，τ_max]范围内进行选择，即

其中，[τ_min，τ_max]是τ^k的最小值和最大值，

是

对于α_N的二阶梯度，表示为

最后求得α_N；

步骤3.3，对于线性去极化参数β子问题，其方程可以转化成以下形式：

其中

上标k表示内部迭代次数，

是

对于β的梯度，表示为

τ^k是一个参数，根据Barzilai-Borwein方法在[τ_mmin，τ_mmax]范围内进行选择，即

其中，

是

对于β的二阶梯度，表示为

最后求得β，最终根据α_N表达式求得α。

Images