远距离超分辨单光子成像重构方法[发明专利]

(12)发明专利申请

(10)申请公布号 CN 110187356 A(43)申请公布日 2019.08.30

(21)申请号 201910519460.8(22)申请日 2019.06.14

(71)申请人 中国科学技术大学

地址 230026 安徽省合肥市包河区金寨路

96号(72)发明人 徐飞虎　黄鑫　汪斌　张强　

潘建伟　(74)专利代理机构 中科专利商标代理有限责任

公司 11021

代理人 李坤(51)Int.Cl.

G01S 17/(2006.01)G01S 7/487(2006.01)

权利要求书2页 说明书10页 附图3页

()发明名称

远距离超分辨单光子成像重构方法(57)摘要

本公开提供一种远距离超分辨单光子成像重构方法，包括：步骤A：在探测信号中寻找回波信号并作门；步骤B：基于三维SPIRALTAP的凸优化算法进行计算，完成远距离超分辨单光子成像重构的数据模型；以及步骤C：将探测信号输入步骤B所建立的远距离超分辨单光子成像重构的数据模型中，完成远距离超分辨单光子成像的重构；利用寻找整体回波信号的门区间和基于三维SPIRALTAP的凸优化算法完成远距离超分辨单光子成像重构，缓解远距离回波信号的提取和存在光斑情况下提高分辨率的问题。

CN 110187356 ACN 110187356 A

权　利　要　求　书

1/2页

1.一种远距离超分辨单光子成像重构方法，包括：步骤A：在探测信号中寻找回波信号并作门；步骤B：基于三维SPIRALTAP的凸优化算法进行计算，完成远距离超分辨单光子成像重构的数据模型；以及

步骤C：将探测信号输入步骤B所建立的远距离超分辨单光子成像重构的数据模型中，完成远距离超分辨单光子成像的重构。

2.根据权利要求1所述的远距离超分辨单光子成像重构方法，所述步骤A，包括：步骤A1：将探测信号中所有像素上的所有光子计数统计到一起，得到分布特性；步骤A2：根据步骤A1所得分布特性对探测信号中噪声做拟合得到噪声计数；以及步骤A3：找到回波信号的门区间。

3.根据权利要求2所述的远距离超分辨单光子成像重构方法，所述步骤A3包括：步骤A31：找到回波信号的大门区间；以及步骤A32：在步骤A31找到的大门区间内进一步找到回波信号的小门区间。4.根据权利要求3所述的远距离超分辨单光子成像重构方法，所述大门区间为将探测信号实际计数率减去步骤A2拟合得到的噪声计数，用一个确定能包含所有探测信号计数的时间宽度对剩余计数进行统计，对应该时间宽度内计数最高的区间。

5.根据权利要求3所述的远距离超分辨单光子成像重构方法，所述小门区间，为在已知的所述大门区间内寻找满足条件的回波信号区间，得到最终的整体回波信号区间；所述满足的条件，为：

Cbin≥c*σ；其中，Cbin为一个时间精度中的剩余计数，σ为剩余计数的标准差；剩余计数Cbin＝Creal-(a*Cnoise_pseudo+b)，其中，Creal表示探测信号计数，Cnoise_pseudo是拟合出来的噪声计数。

6.根据权利要求1所述的远距离超分辨单光子成像重构方法，所述步骤B，包括：步骤B1：建立远距离情形下的正向探测模型；以及步骤B2：求解步骤B1所建立的探测模型。

7.根据权利要求1所述的远距离超分辨单光子成像重构方法，所述探测信号服从非均匀泊松分布。

8.根据权利要求6所述的远距离超分辨单光子成像重构方法，所述步骤B1中所述探测模型，表达如下：

Y～Poisson(h*X+B)        (3)；式中X表示场景，Y表示探测信号，X、Y均为三维矩阵(nx*ny*nt)，Yx，y)y，t表示像素点(x，处到达时间处于t这个时间精度里的光子的数量；关于场景X，对于任意像素点(x，y)处对应的t方向上的向量，仅有一处为非零值，该非零值的位置对应场景X的深度，该非零值的大小对应场景X的反射率信息。

9.根据权利要求6所述的远距离超分辨单光子成像重构方法，所述步骤B2中，定义L(X；Y，h，B)为X的负对数似然函数，则：

2

CN 110187356 A

权　利　要　求　书

2/2页

其中，X是一个用以描述场景信息的三维矩阵；Y是一个用以描述实测数据的三维矩阵；h是一个用以描述时空核函数的三维矩阵；B是一个用以描述背景噪声的三维矩阵；TV是惩罚约束项，nx、ny、nt是自变量矩阵X的三维尺寸；β是一个正实数，用以描述调节惩罚项的影响大小。

10.根据权利要求1至9任一项所述的远距离超分辨单光子成像重构方法，适用于待测目标的深度跨越范围为0～1500m。

3

CN 110187356 A

说　明　书

远距离超分辨单光子成像重构方法

1/10页

技术领域

[0001]本公开涉及激光雷达领域，尤其涉及一种远距离超分辨单光子成像重构方法。背景技术

[0002]随着单光子测量技术的快速发展，利用时间相关单光子测量技术(TCSPC)[1]对目标进行三维成像得到了很大发展[2-3]。TCSPC技术具有传统深度测量手段所无法达到的时间分辨能力(≤102ps)[1]。这一特性使得TCSPC技术在分子光谱，激光雷达和三维成像等多个方向得到大量的应用。在激光雷达领域，利用TCSPC技术进行三维成像已经有了长足发展[4-9]。但由于该技术单光子测量的特性，在进行三维成像时，每个像素点平均的信号光子计数只有十几甚至几个。这种情形下，Nyquist(奈奎斯特)采样定理已经难以使用，且在大量的噪声和暗计数背景中，找出这几个信号光子计数也存在着很大的困难。这些问题的存

8-12]

在，使得适用于单光子计数三维成像的算法应运而生[2，。[0003]2014年，来自MITs Research Laboratory of Electronics的Ahmed Kirmani等人研发了一套通过光栅式扫描的首光子成像(First-photon imaging，FPI)系统[9]。该套系统可以仅仅依靠每个像素点处探测器接收到的第一个光子来进行高质量的三维图像重建，开创了少光子成像的先河。但该系统在实际应用上存在较大的局限性，逐点扫描的成像方式使得成像速度受到。[0004]此后，Shin等人将首光子成像的框架由逐点扫描拓展到了阵列的情形[8][10][11]，使得少光子成像技术向实用化迈出了坚实的一步。Shin等人基于光子的量子性和阵列的情形提出了相应的测量模型，将重构问题描述为一个求解极大似然估计的凸优化问题[11]。在此类凸优化问题中，光子的量子性通常表现为泊松分布特性。在成像领域，早在2012年便已提出成熟的的泊松分布凸优化求解算法，即SPIRALTAP[13]。Shin等人利用SPIRALTAP包，有效地实现了少光子条件下的重构问题。[0005]由于光的量子性，探测端在探测少光子时会受到泊松噪声的影响。而对于自然场景中的目标物体，其反射率和空间深度分布，通常具有很强的空间关联性，即相近的空间点，反射率和空间深度通常相近。利用该特点，可以通过TV(total-variation)[14]约束方式来优化成像结果。经验标明，TV的约束方式对于抑制泊松噪声有很好的效果。Shin等人通过添加约束项TV(total-variation)[14]来抑制泊松噪声，从而获得了较好的重构效果。[0006]更进一步，Rapp等人考虑环境噪声很大的情况，提出了相应的滤噪手段[12]。当信噪比很低时，对采集到的信号先进行滤噪处理，再进行重构，以此提高重构精度。[0007]随着这些算法的发展，信号光子计数可以做到～1ppp(光子每像素)，信噪比也可以做到极低(SBR～0.04)。这样的条件，使得将TCSPC技术应用于远距离三维成像变得很有潜力。因为在几十公里的距离上，激光传播路径上的衰减极大，回波的信号光子在单光子量级，且掩盖在大量噪声当中。

[0008]已有的工作将远距离单光子激光雷达做到的10km[5]。但在10km这个距离上，其成像结果较为简单，只是个斜坡。我们考虑在更远的距离上，对目标物体进行完整的三维成

4

CN 110187356 A

说　明　书

2/10页

像。在这样的问题当中，存在以下问题：很大比重的背景噪声；由于光斑的存在，每个像素中将存在多个不同的深度信息。现有的算法大多注重于如何从大量背景噪声中提取出极少的信号光子信息，但对于光斑的问题考虑却极少。在几十公里的距离上，光斑带来的问题将是不可回避的问题。

[0009]激光雷达采用脉冲激光主动照明，该脉冲激光的脉冲周期通常会远远大于目标物体的深度跨越范围，保证目标的回波信号能完成的分布在一个周期内。单光子探测器与激光同步计数，故总体的计数区间为脉冲激光的脉冲周期。但考虑到回波信号的范围远远小于脉冲周期，故在信号时间范围外的单光子计数均来源于环境的背景光和单光子探测器的暗计数与后脉冲。在信噪比较高的情形下，这些计数不会给成像带来太多影响，但在信噪比极低的情形下，整个周期上，噪声光子与信号光子的计数比率可以达到102∶1，甚至更低。此时噪声计数将对图像重构产生极大的影响。很直接的解决方案便是将信号范围外的所有计数抛弃掉，用余下的计数来做图像重构，即作门的方法。在已有的工作中[8]，已将作门的方法引入。但在该工作中，场景为几米的距离范围，存在大量场景带来的回波，整体回波信号较为复杂，如图1所示，其寻峰方案适用场景比较有限。而在远距离情形下，回波信号具有一些特征，已发表工作[8]中的方案在该情形下并不适用。我们根据回波信号的特征，提出了自己的寻峰和作门方案。[0010]另一方面，激光在经过远距离传播后会出现发散的情况，导致最后照射到物体上的不再是一个一个理想的激光点而是一个一个的圆形光斑，不同于激光点只对应单一深度，一个光斑可以包含多个深度。

[0011]文献[8]-[12]等一系列工作所考量的都是理想的激光点情形，所提出的一系列测量模型和相匹配算法都是基于单一深度去考虑的，以致于在远距离成像中这一系列的成像算法表现欠佳。[0012][1].O′Connor D.Time-correlated single photon counting[M].Academic Press，2012.

[0013][2].Buller G，Wallace A.Ranging and three-dimensional imaging using time-correlated single-photon counting and point-by-point acquisition[J].IEEE Journal of selected topics in quantum electronics，2007，13(4)：1006-1015.[0014][3].Altmann Y，McLaughlin S，Padgett M J，et al.Quantum-inspired computational imaging[J].Science，2018，361(03)：eaat2298.[0015][4].A.McCarthy，X.Ren，A.Della Frera，N.R.Gemmell，N.J.Krichel，C.Scarcella，A.Ruggeri，A.Tosi，and G.S.Buller，“Kilometer-range depth imaging at 1,550nm wavelength using an InGaAs/InP single-photon avalanche diode detector，”Opt.Express 21(19)，22098-22113(2013).[0016][5].Pawlikowska A M，Halimi A，Lamb R A，et al.Single-photon three-dimensional imaging at up to 10kilometers range[J].Optics express，2017，25(10)：11919-11931.[0017][6].Kang Y，Li L，Liu D，et al.Fast long-range photon counting depth imaging with sparse single-photon data[J].IEEE Photonics Journal，2018.[0018][7].Tobin R，Halimi A，McCarthy A，et al.Long-range depth profiling of 

5

CN 110187356 A

说　明　书

3/10页

camouflaged targets using single-photon detection[J].Optical Engineering，2017，57(3)：031303.[0019][8].Shin D，Xu F，Venkatraman D，et al.Photon-efficient imaging with a single-photon camera[J].Nature communications，2016，7：12046.[0020][9].Kirmani A，Venkatraman D，Shin D，et al.First-photon imaging[J].Science，2014，343(6166)：58-61.[0021][10].Shin D，Kirmani A，Goyal V K，et al.Computational 3D and reflectivity imaging with high photon efficiency[C]//Image Processing (ICIP)，2014IEEE International Conference on.IEEE，2014：46-50.[0022][11].Shin D，Kirmani A，Goyal V K，et al.Photon-efficient computational3-D and reflectivity imaging with single-photon detectors[J].IEEE Transactions on Computational Imaging，2015，1(2)：112-125.[0023][12].Rapp J，Goyal V K.A few photons among many：Unmixing signal and noise for photon-efficient active imaging[J].IEEE Transactions on Computational Imaging，2017，3(3)：445-459.[0024][13].Harmany Z T，Marcia R F，Willett R M.This is SPIRAL-TAP：Sparse Poisson intensity reconstruction algorithms-theory and practice[J].IEEE Transactions on Image Processing，2012，21(3)：1084-1096.[0025][14].Rudin L I，Osher S，Fatemi E.Nonlinear total variation based noise removal algorithms[J].Physica D：nonlinear phenomena，1992，60(1-4)：259-268.[0026][15].Shin D，Shapiro J H，Goyal V K.Photon-efficient super-resolution laser radar[C]//Wavelets and Sparsity XVII.International Society for Optics and Photonics，2017，10394：1039409.[0027]Snyder D L，Miller M I.Random point processes.j[J].Willey&Sons，New York，1975.

[0028]公开内容[0029](一)要解决的技术问题[0030]基于上述问题，本公开提供了一种远距离超分辨单光子成像重构方法，以缓解现有技术中远距离成像时背景噪声占比大，光斑影响大导致的远距离成像方法效果差等技术问题。[0031](二)技术方案[0032]在本公开中，提供一种远距离超分辨单光子成像重构方法，包括：[0033]步骤A：在探测信号中寻找回波信号并作门；[0034]步骤B：基于三维SPIRALTAP的凸优化算法进行计算，完成远距离超分辨单光子成像重构的数据模型；以及[0035]步骤C：将探测信号输入步骤B所建立的远距离超分辨单光子成像重构的数据模型中，完成远距离超分辨单光子成像的重构。[0036]在本公开实施例中，所述步骤A，包括：[0037]步骤A1：将探测信号中所有像素上的所有光子计数统计到一起，得到分布特性；

6

CN 110187356 A[0038]

说　明　书

4/10页

步骤A2：根据步骤A1所得分布特性对探测信号中噪声做拟合得到噪声计数；以及

[0039]步骤A3：找到回波信号的门区间。[0040]在本公开实施例中，所述步骤A3包括：[0041]步骤A31：找到回波信号的大门区间；以及[0042]步骤A32：在步骤A31找到的大门区间内进一步找到回波信号的小门区间。[0043]在本公开实施例中，所述大门区间为将探测信号实际计数率减去步骤A2拟合得到的噪声计数，用一个确定能包含所有探测信号计数的时间宽度对剩余计数进行统计，对应该时间宽度内计数最高的区间。[0044]在本公开实施例中，所述小门区间，为在已知的所述大门区间内寻找满足条件的回波信号区间，得到最终的整体回波信号区间；所述满足的条件，为：Cbin≥c*σ；其中，Cbin为一个时间精度中的剩余计数，σ为剩余计数的标准差；剩余计数Cbin＝Creal-(a*Cnoise_pseudo+b)，其中，Creal表示探测信号计数，Cnoise_pseudo是拟合出来的噪声计数。[0045]在本公开实施例中，所述步骤B，包括：[0046]步骤B1：建立远距离情形下的正向探测模型；以及[0047]步骤B2：求解步骤B1所建立的探测模型。[0048]在本公开实施例中，所述探测信号服从非均匀泊松分布。[0049]在本公开实施例中，所述步骤B1中所述探测模型，表达如下：[0050]Y～Poisson(h*X+B)  (3)；[0051]式中X表示场景，Y表示探测信号，X、Y均为三维矩阵(nx*ny*nt)，Yx，y，t表示像素点(x，y)处到达时间处于t这个时间精度里的光子的数量；关于场景X，对于任意像素点(x，y)处对应的t方向上的向量，仅有一处为非零值，该非零值的位置对应场景X的深度，该非零值的大小对应场景X的反射率信息。[0052]在本公开实施例中，所述步骤B2中，定义L(X；Y，h，B)为X的负对数似然函数，则：

[0053]

[00]

其中，X是一个用以描述场景信息的三维矩阵；Y是一个用以描述实测数据的三维矩阵；h是一个用以描述时空核函数的三维矩阵；B是一个用以描述背景噪声的三维矩阵；TV是惩罚约束项，nx、ny、nt是自变量矩阵X的三维尺寸；β是一个正实数，用以描述调节惩罚项的影响大小。

[0056]在本公开实施例中，所述的远距离超分辨单光子成像重构方法，适用于待测目标的深度跨越范围为0～1500m。[0057](三)有益效果

[0058]从上述技术方案可以看出，本公开远距离超分辨单光子成像重构方法至少具有以下有益效果其中之一或其中一部分：

7

[0055]

CN 110187356 A[0059][0060][0061]

说　明　书

5/10页

(1)可以应对背景噪声十分强烈的情形，具有全天时的工作能力；(2)具有一定的透雾、透雨的能力，可以工作在较为恶劣的天气条件下；(3)配合硬件的精细采集，可以实现超越光学口径极限的分辨率。

附图说明

[0062]图1为现有技术中探测信号数据的寻峰方案示意图。

[0063]图2为本公开实施例远距离超分辨单光子成像重构方法中探测信号的数据分布情况示意图。

[00]图3为本公开实施例对图2所示的探测信号数据的处理过程示意图；其中，图3(a)为对探测信号进行噪声拟合的结果；图3(b)为探测信号扣除拟合噪声的结果；图3(c)为根据(b)中结果取出的包含回波信号的大门区间；图3(d)为在大门区间内进一步精细寻找包含回波信号的小门区间；

[0065]图4为采用不同的成像方法对3.75km外电线杆的成像结果对比示意图；其中，图4(a)显示了用被动天文相机拍到的画面，图4(b)是直接对探测信号做平均得到的深度图，图4(c)是用本公开的方法对探测信号进行重构得到的深度图，图4(d)为普通激光雷达对分辨板进行成像的模拟结果；图4(e)为远距离超分辨单光子成像重构方法配合视场交叠约为1/2时扫描得到的模拟结果；图4(f)为远距离超分辨单光子成像重构方法配合视场交叠约为3/4时扫描得到的模拟结果。

[0066]图5为本公开实施例对3km外的自然景观和21.5km外人造景观的成像结果示意图；其中，图5(a)为自然景观的辅助相机被动成像结果示意图，图5(b)为采用本公开的远距离超分辨单光子成像重构方法得到的自然景观的深度图数据重现的三维场景；图5(c)为采用本公开的远距离超分辨单光子成像重构方法得到的自然景观的深度图；图5(d)为人造景观的辅助相机被动成像结果示意图，图5(e)为采用本公开的远距离超分辨单光子成像重构方法得到的人造景观的深度图数据重现的三维场景；图5(f)为采用本公开的远距离超分辨单光子成像重构方法得到的人造景观的深度图。

[0067]图6为本公开实施例远距离超分辨单光子成像重构方法的流程图。

具体实施方式

[0068]本公开提供了一种远距离超分辨单光子成像重构方法，利用寻找整体回波信号的门区间和基于三维SPIRALTAP的凸优化算法完成远距离超分辨单光子成像重构，缓解远距离回波信号的提取和存在光斑情况下提高分辨率的问题。[0069]为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开进一步详细说明。[0070]在本公开实施例中，提供一种远距离超分辨单光子成像重构方法，结合图2至图6所示，所述远距离超分辨单光子成像重构方法，包括：[0071]步骤A：在探测信号中寻找回波信号并作门；[0072]在本公开实施例中，所述步骤A包括：[0073]步骤A1：将探测信号中所有像素上的所有光子计数统计到一起，得到分布特性；[0074]远距离情形下，由于目标特性和系统的特点，所有像素累加后的探测信号具有很

8

CN 110187356 A

说　明　书

6/10页

多独特的特点。我们先将探测信号的所有像素上的所有光子计数统计到一起，对总体的计数的分布特性做分析。这样做是因为单个像素上的积分时间过短，信号光子数过少，也不包含目标物体的整体深度区间。[0075]远距离情形下，目标物体通常是大型的建筑目标，目标的深度跨越范围比较大(～102m)，其对应的回波信号通常也在102ns的量级；由于目标物体的深度分布的不规则性，回波信号的形状也不再是简单的高斯型，无法通过特性高斯拟合寻找；远距离情形下，通常需要大口径的望远镜，除目标物体上的一次漫反射回波外，基本不会有其他回波信号被收集到，故整体的回波信号很明显，不会存在信号外的峰；由于远距离的衰减，回波信号在单光子量级，故整体的回波信号光子计数相对于所有的噪声计数是很少的；另外，在收发同轴的系统中，存在系统噪声，由于需要高功率激光器，激光器的底噪在后光路中的后向散射会带来大量的背景噪声计数，这些计数不携带任何目标物体的信息，且占比极大，甚至会高于环境光和暗计数。

[0076]这些特点导致远距离情形下的单光子成像计数分布不同于常见的回波信号，如图2所示，是将所有像素的计数统计到一起后作出的柱状图。该图中呈现逐渐上升的底噪便是因为采用了收发同轴的方式，激光器中用于主放大的泵浦源在一个周期中不断积蓄能量，工作物质的自发辐射光随着能量的积蓄不断增强，带来越来越多的噪声计数。[0077]步骤A2：根据步骤A1所得分布特性对探测信号中噪声做拟合得到噪声计数；[0078]以一定的bin(实验时间精度，比如1ns的时间精度，体现在坐标轴中，即横轴最小一格为1ns)的宽度来统计所有的探测信号计数分布特性，根据噪声的分布特性对噪声做拟合。

[0079]由于整体的回波信号计数占总探测信号计数的比例很小，故通过多项式拟合，便可以使得拟合得到的函数接近噪声的分布，如图3-(a)所示。[0080]步骤A3：找到回波信号的门区间；包括：[0081]A31：找到回波信号的大门区间；

[0082]将探测信号实际计数率减去步骤A2拟合得到的噪声计数，用一个确定能包含所有探测信号计数的时间宽度对剩余计数进行统计，对应该时间宽度内计数最高的区间，我们称该区间为大门区间；

[0083]值得注意的是大门区间的数目不限于1，对于多个比较的目标，会存在多个完全的信号峰，便可以采用多个大门区间，例如a，b，c三个参数，是考虑到不同的噪声情形下，可能需要适当调节对噪声的扣除，对于积分时间长，噪声计数抖动小的，a，b，c可以依次设为1，0，1，而对于积分时间短噪声抖动比较大的情形，可适当调大这些参数。[0084]步骤A32：在步骤A31找到的大门区间内进一步找到回波信号的小门区间；[0085]进一步的，然后在这(几个)大门区间里，去掉大门区间里计数很小的bin。，此时提供可调节的参数a，b，c，则剩余计数为Cbin＝C-real-(a*Cnoise_pseudo+b)，其中Creal表示实际测出来的一个bin里的计数，即探测信号计数，Cnoise_pseudo是拟合出来的噪声计数；计算剩余计数的标准差σ，对于剩余计数，满足Cbin≥c*σ的bin记为回波信号区间，其中Cbin为一个bin中的剩余计数。在已知的大门区间内寻找满足条件的回波信号区间，得到最终的整体回波信号区间，我们称该区间为小门区间。[0086]对于图2中所示的计数情形，可以通过该寻峰作门的方式，得到如图3所示的门区

9

CN 110187356 A

说　明　书

7/10页

间。如图3所示，由于信号相对于噪声的计数是很少的，可以用Matlab内置的拟合方式，二次型拟合，便可简单实现对噪声的拟合，并很快找出信号峰所在时间段。[0087]步骤B：基于三维SPIRALTAP的凸优化算法进行计算，完成远距离超分辨单光子成像重构的数据模型；包括：[0088]步骤B1：建立远距离情形下的正向探测模型；[00]通常，在不考虑激光发散和视场角发散的情况下，我们认为接收端与场景中的像素点一一对应，对场景中任一像素点(x，y)，以h(t)形式的脉冲光照亮一次该像素点，与该像素点相对应的探测器所能接收到的光通量r具有如下形式[11]：[0090]r(t；x，y)＝ref(x，y)·h(t-2dep(x，y)/c)+b  (1)[0091]其中，ref表示反射率，dep表示深度，b是背景光强。[0092]然而，在远距离的情形下，由于视场角的发散，一个像素点会对应相对较大的视场(50公里外，10urad的发散角，一个像素会对应50cm的视场)，探测端会接收到从该视场所包含的多个不同的深度所反射回来的光子。记探测端视场中心对应的位置为(x，y)，该视场为FOV，(x′，y′)为FOV中的任一位置，则探测端接收到的光通量可以描述为[13]：[0093]R(t；x，y)＝∫，y-y′)ref(x′，y′)ht(t-2dep(x′，y′)/c)dx′dy′+b  (2)；x′，y′∈FOVhxy(x-x′[0094]上式中，R(t；x，y)表示像素点(x，y)处收集到的到达时间为t的光子数密度，hxy表示光斑在视场内的强度分布，ht表示脉冲光的脉宽和探测器抖动带来的到达时间上的分布(一般简单视作高斯分布)，ref是反射率，b是背景光强。[0095]实际应用中，时间和空间上都存在一个最小精度，而非连续值，这种情况下公式(2)存在离散的表达形式，光通量R可以描述为场景的反射率和深度X与时空核函数h的卷积R＝h*X+B，其中h是hxy和ht的直积，表示时空核函数对X的卷积操作，即光强分布hxy和脉冲时间分布ht对测量结果带来的影响；也是同样规格的一个三维矩阵，表示噪声的影响。同时我们根据光子探测的相关理论[16]，知道探测信号Y服从一个非均匀泊松分布，从而我们可以得到：

[0096]Y～Poisson(h*X+B)  (3)；[0097]式中式中X表示场景，Y表示探测信号，X、Y均为三维矩阵(nx*ny*nt)，Yx，y，t表示像素点(x，y)处到达时间处于t这个时间精度里的光子的数量；每个矩阵元的值表示像素点(x，y)处到达时间处于nt这个bin里的光子的数量；X也是一个三维矩阵(nx*ny*nt)，用来描述场景，这是一个稀疏矩阵，对于任意像素点(x，y)处对应的t方向上的向量，仅有一处为非零值，该非零值的位置对应场景X的深度，该非零值的大小对应场景X的反射率信息。[0098]步骤B2：求解步骤B1所建立远距离探测模型；

[0099]我们所需要做的事情是根据探测信号结果Y估计场景反射率和深度，也就是X。我们采用极大似然估计的思想，同时根据式(3)，我们定义关于X的负对数似然函数L(X；Y，h，B)，上述估计问题可以描述为一个凸优化问题：

[0100]

为了求解上述凸优化问题，我们在[13]中的SPIRATAP包的基础上做了一定的修改，将二维情形拓展到了三维情形。求解极值的思路依旧沿用[13]中添加约束项迭代求解的思路，同时TV项的选取也是基于[13]中的TV表示的拓展：

10

[0101]

CN 110187356 A

说　明　书

8/10页

[0102]

上式中，X是一个用以描述场景信息的三维矩阵；Y是一个用以描述实测数据的三

维矩阵；h是一个用以描述时空核函数的三维矩阵；B是一个用以描述背景噪声的三维矩阵；TV是惩罚约束项，具体表达式为公式(5)的第二行，nx、ny、nt是自变量矩阵X的三维尺寸；β是一个正实数，用以描述调节惩罚项的影响大小。[0104]步骤C：将探测信号输入步骤B所建立的远距离超分辨单光子成像重构的数据模型中，完成远距离超分辨单光子成像的重构。[0105]如图4所示，我们对约4公里(3.75km)外的一个电线杆进行了主动成像，图4(a)显示了当时用被动天文相机拍到的画面，图4(b)是直接对探测信号做平均得到的深度图，图4(c)是用本公开的方法对探测信号进行重构得到的深度图，可以看到，本公开的方法可以恢复出一个较为清晰的电线杆的深度图。[0106]如图4(d)至图4(f)所示，我们对远距离超分辨单光子成像重构方法在超分辨上的表现进行了研究，我们借助一个分辨板进行成像模拟，图4(d)给出了普通激光雷达对分辨板进行成像的模拟结果；图4(e)、图4(f)给出了远距离超分辨单光子成像重构方法配合交叠扫描得到的模拟结果，其中图4(e)对应的视场交叠约为1/2，图4(f)对应的视场交叠约为3/4，可以看到，随着交叠程度的越来越高，恢复结果的分辨程度也越高。

[0107]图5展示了本公开所述远距离超分辨单光子成像重构方法在现实场景当中的具体运用，将探测信号导入编写好的相关重构程序，即可获得重构结果。为进一步说明该算法的实用性，图5分别展示远距离超分辨单光子成像重构方法对自然景观和人造景观的重构结果以显示该重构算法对不同场景较好的适应性。如图5所示，分别是对3km外的自然景观进行成像和对21.5km外人造景观进行成像。图5(a)和图5(d)为我们展示了辅助相机被动成像的结果；图5(c)和图5(f)为我们展示了该算法重构出来的深度图；图5(b)和图5(e)显示了利用重构得到的深度图数据重现出的三维场景，可以看到远距离超分辨单光子成像重构方法对自然景观和人造景观都有较好的重构效果。[0108]至此，已经结合附图对本公开实施例进行了详细描述。需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换。[0109]依据以上描述，本领域技术人员应当对本公开远距离超分辨单光子成像重构方法有了清楚的认识。[0110]综上所述，本公开提供了一种远距离超分辨单光子成像重构方法，利用寻找整体回波信号的门区间和基于三维SPIRALTAP的凸优化算法完成远距离超分辨单光子成像重构，缓解远距离回波信号的提取和存在光斑情况下提高分辨率的问题。[0111]还需要说明的是，实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向，并非用来本公开的保护范围。贯穿附图，相同的元素由

[0103]

11

CN 110187356 A

说　明　书

9/10页

相同或相近的附图标记来表示。在可能导致对本公开的理解造成混淆时，将省略常规结构或构造。

[0112]并且图中各部件的形状和尺寸不反映真实大小和比例，而仅示意本公开实施例的内容。另外，在权利要求中，不应将位于括号之间的任何参考符号构造成对权利要求的。

[0113]除非有所知名为相反之意，本说明书及所附权利要求中的数值参数是近似值，能够根据通过本公开的内容所得的所需特性改变。具体而言，所有使用于说明书及权利要求中表示组成的含量、反应条件等等的数字，应理解为在所有情况中是受到「约」的用语所修饰。一般情况下，其表达的含义是指包含由特定数量在一些实施例中±10％的变化、在一些实施例中±5％的变化、在一些实施例中±1％的变化、在一些实施例中±0.5％的变化。[0114]再者，单词“包含”不排除存在未列在权利要求中的元件或步骤。位于元件之前的单词“一”或“一个”不排除存在多个这样的元件。[0115]说明书与权利要求中所使用的序数例如“第一”、“第二”、“第三”等的用词，以修饰相应的元件，其本身并不意味着该元件有任何的序数，也不代表某一元件与另一元件的顺序、或是制造方法上的顺序，该些序数的使用仅用来使具有某命名的一元件得以和另一具有相同命名的元件能做出清楚区分。[0116]此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无于以上所列，且可根据所需设计而变化或重新安排。并且上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。

[0117]本领域那些技术人员可以理解，可以对实施例中的设备中的模块进行自适应性地改变并且把它们设置在与该实施例不同的一个或多个设备中。可以把实施例中的模块或单元或组件组合成一个模块或单元或组件，以及此外可以把它们分成多个子模块或子单元或子组件。除了这样的特征和/或过程或者单元中的至少一些是相互排斥之外，可以采用任何组合对本说明书(包括伴随的权利要求、摘要和附图)中公开的所有特征以及如此公开的任何方法或者设备的所有过程或单元进行组合。除非另外明确陈述，本说明书(包括伴随的权利要求、摘要和附图)中公开的每个特征可以由提供相同、等同或相似目的的替代特征来代替。并且，在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。[0118]类似地，应当理解，为了精简本公开并帮助理解各个公开方面中的一个或多个，在上面对本公开的示例性实施例的描述中，本公开的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而，并不应将该公开的方法解释成反映如下意图：即所要求保护的本公开要求比在每个权利要求中所明确记载的特征更多的特征。更确切地说，如下面的权利要求书所反映的那样，公开方面在于少于前面公开的单个实施例的所有特征。因此，遵循具体实施方式的权利要求书由此明确地并入该具体实施方式，其中每个权利要求本身都作为本公开的单独实施例。

[0119]以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本公开的具体实施例而已，并不用于本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保

12

CN 110187356 A

说　明　书

10/10页

护范围之内。

13

CN 110187356 A

说　明　书　附　图

1/3页

图1

图2

14

CN 110187356 A

说　明　书　附　图

2/3页

图3

图4

15

CN 110187356 A

说　明　书　附　图

3/3页

图5

图6

16