High-resolution three-dimensional imaging based on all-fiber photon-counting Lidar system

光子计数激光雷达系统（Lidar）可以实现光子量级的高灵敏度探测、μrad级的高空间分辨率，近年来受到广泛关注^[1]。激光能量在远距离大气传输过程中衰减严重，回波信号甚至达到光子量级，传统探测方式难以识别。光子计数激光雷达发射高重频、窄脉宽的脉冲激光，利用单光子探测方式响应光子能量级别的微弱信号，通过计算激光飞行时间，获得几百甚至几十公里外的目标距离信息，进而重建三维图像^[2-3]。目前，其应用已覆盖远距离成像^[4-5]、大气探测^[6]、地形测绘^[7-8]等领域。

在光子计数激光雷达远距离成像领域，林肯实验室自2001年开展了大量研究工作，先后研制了多型机载激光雷达原理样机^[9]，在2~9 km空域完成飞行试验，验证了全天时三维成像能力^[10-11]。2015年，Harris公司推出了世界上第一款商用盖革模式激光雷达系统，探测距离达到6 km^[12]。这些研究均是基于盖革模式单光子雪崩二极管（GM-APD）阵列，受到阵列规模、暗计数、视场角等性能制约，难以进一步提升视场。2011年，NASA验证了基于4×4光纤阵列的机载测绘激光雷达，工作波长为1030 nm，探测距离达到10 km^[13]。国内在远距离成像方面的起步较晚，中国科学院上海技术物理研究所搭建的基于51波束光纤阵列的1064 nm光子计数激光雷达，可实现3 km以内的机载三维测绘^[14]。华东师范大学搭建的1550 nm收发一体式光子计数激光雷达，可以实现21 km的测距^[4]。中国科学技术大学搭建的1550 nm光子计数激光雷达系统，实现了45 km远距离成像^[15]。在这些已报道的光子计数激光雷达系统中，收发光学系统大多采用空间光传输方式，或空间光与光纤组合传输方式，增加了系统不稳定性。此外，扫描方式大多为摆镜扫描。为实现大范围扫描，摆镜口径通常较大，增大了系统尺寸；另一方面，摆镜旋转会导致反射图像发生几何畸变，影响成像质量。

文中设计了一种工作波长为1064 nm的全光纤光子计数激光成像雷达系统，通过整机扫描方式和亚像素扫描方法对3 km以上目标进行成像探测，空间分辨率超越衍射极限，并通过噪声阈值自适应的图像重建算法，重建多距离目标信息。该系统适用于远距离探测、地形测绘等领域。

文中设计的光子计数激光雷达系统，通过发射脉冲激光，记录脉冲飞行时间来探测目标距离，通过二维扫描来实现大视场成像。系统工作原理见图1，激光器在内部电路驱动下产生高重频、窄脉宽的脉冲激光，经准直扩束后出射。发射激光被目标反射后，进入接收望远镜。最终回波信号进入硅基单光子探测器，完成光电转换。时间相关单光子计数（TCSPC）模块记录激光发射时刻t₀和回波到达时刻t₁，则目标距离R为：

Figure 1.  Block diagram of photon-counting Lidar system

式中：c为光速。

利用TCSPC技术对单光子信号进行探测，并且记录光子响应的时间，具有极高的时间分辨率。该技术将时间轴分为离散的时间区间，时间区间的大小决定了系统的精度。当探测器探测到一个或多个光子，产生一个响应输出，记录下该次响应发生的时间，并对该时间区间内的光子计数值加1。当进行大量重复脉冲探测后，可以统计得到响应时间对应光子计数数目的统计直方图，直方图峰值表征回波幅度，峰值所在的时刻表征回波到达时间。直方图波形存在着一定的波形宽度，对应的是一定的时间跨度，即距离误差，所以需要选取一定的算法来得到最准确的距离信息。

（1）静态视场角

通常，静止状态下光纤光学系统的空间分辨率由发射视场角FOV决定，可近似表达为：

式中：MFD为光纤模场直径；f为望远镜焦距。

（2）亚像素扫描成像

在远距离成像条件下，即使是很小的发散角也会在目标位置产生较大的光斑。为提高空间分辨率，文中采用如图2所示的亚像素扫描方法^[16]。利用高精度转台，控制扫描角度步进小于FOV，对静态视场覆盖区域进行n（n>1）次采样。整个场景扫描完成后，通过图像重建算法实现三维成像，分辨率为静态视场角的1/n。

Figure 2.  Sub-pixel scanning method

文中设计的光子计数激光雷达系统主要由发射光学系统、接收光学系统、时间相关单光子计数模块、扫描系统和上位机组成。如图1所示，发射和接收光学系统为全光纤结构。发射光学系统由光纤激光器、准直镜构成。接收光学系统由滤波器、接收望远镜和单光子探测器组成。伺服系统由二维转台、伺服控制系统、角度传感器、支架等构成。上位机负责控制激光器的脉宽、重频、泵浦功率、TCSPC采集时间和采集模式，以及二维转台的转速和步进，并记录角度传感器的实时角度，完成数据采集和图像重建。

光子计数激光成像雷达的典型工作波长包括532 nm、1064 nm、1030 nm和1550 nm等。波长对系统性能起到关键作用，例如太阳在532 nm的辐射量是1064 nm的三倍以上，因此大多数光子计数激光成像雷达工作于近红外波段。综合考虑太阳背景辐射、单光子探测器的光谱响应效率和技术成熟度，文中系统采用工作波长1064.4 nm的脉冲光纤激光器，重频为40 kHz。激光器采用主振荡器加功率放大器（MOPA）结构设计，经过优化设计的多级光放大，实现峰值功率6 μJ，脉宽5 ns的激光输出，线宽约为0.04 nm。激光器输出尾纤为芯径9/125 μm的SMF-28e光纤。发射准直镜焦距80 mm。根据公式（2）计算可得，系统的静态视场角FOV约为131 μrad。

接收光学系统中，设计接收视场角为发射光束发散角的两倍。接收望远镜焦距为40 mm，孔径为20.5 mm。为了降低背景辐射噪声，接收准直镜前安装了带宽1 nm的滤光片，通过螺纹固紧于接收望远镜的前端。接收信号直接耦合进SMF-28e光纤，进入探测器，完成光电转换。

目前，能实现光子量级响应的探测器主要包括光电倍增管（PMT）、盖革模式单光子雪崩二极管（GM-APD）和超导探测器。其中，PMT对近红外响应较低，体积较大，不适用于光子计数激光雷达。盖革模式单光子雪崩二极管（SPADs）是一种概率器件，可以一定的概率对单个光子发生雪崩。现有SPAD单光子探测器的货架产品主要采用Si基或InGaAs探测器^[17]。InGaAs探测器的死时间和暗计数都较大。Si基探测器较为成熟，响应波段主要为可见光，对1064 nm波长的探测效率为2%，对1550 nm波长几乎无响应。其他报道的提高探测效率的方法包括超导探测器和上转换技术等^[18]，但仍处于实验室研究阶段，尚不成熟。因此，文中采用Si基单光子探测器，暗计数为100 cps，死时间为22 ns。

光学系统主要参数见表1。

为扩大成像视场，现有光子计数激光雷达系统基本采用摆镜扫描方法：在出射/接收光路中加入反射镜，通过旋转反射镜，实现水平和俯仰方向的光束扫描。该方法存在以下两个问题：一是反射镜口径较大，二是反射图像存在几何畸变。反射镜示意图见图3，对于收发同轴光学系统，设收发望远镜直径为w₁，反射镜与激光发射光路的夹角为θ；对于收发分置系统，w₁为发射和接收望远镜的最大包络距离，其他与收发同轴光学系统相同。则反射镜直径w₂应大于w₁/sinθ，以避免边缘能量损失。因此，望远镜口径越大，扫描角度越大，需要的反射镜口径越大。

Figure 3.  Beam reflection by a mirror

此外，摆镜引起的几何畸变如图4所示，摆镜的两个相互垂直的旋转轴为r轴和p轴，p轴在xy平面内，目标在xz平面内。如图4（a）所示，当r轴与z轴平行，p轴与x轴呈45°时，目标经摆镜反射后，反射像上下不变，左右翻转，与目标镜像对称。当r轴保持与z轴平行，摆镜围绕r轴转动时，如图4（b）所示，与z轴平行的目标的反射像在z轴上的投影长度不变；但与x轴平行的目标反射后长度发生拉伸或收缩。当p轴保持与x轴呈45°时，摆镜围绕p轴转动，如图4（c）所示，与x轴平行的目标的反射像在y轴上的投影长度不变；与z轴平行的目标经反射镜反射后，在yz平面内发生旋转。

Figure 4.  Geometry distortion caused by a swing mirror

为了避免反射像的几何畸变，系统采用整机扫描方案，即光学系统发射的激光直接照射目标，通过将光学系统整体旋转，实现水平和俯仰方向的扫描。如图5所示，将一个一维转台水平放置，以实现水平方向扫描；将另一个一维转台垂直安装于水平转台上，以实现俯仰方向扫描。光学系统中，发射准直镜、接收望远镜、滤波片及其调整、夹持结构件安装于俯仰转台转动面的转接板上，激光器和单光子探测器安装于俯仰转台底座的转接板上。转台水平扫描范围为360°，俯仰扫描范围为±30°。俯仰和水平方向的角度分辨率均为0.002°（约等于35 μrad）。因此，扫描模式的空间分辨率理论上小于凝视模式的1/3。

Figure 5.  Design of the rotation platform

伺服控制过程如下：第一步，开启激光器和TCSPC，将激光器内部时钟作为TCSPC同步时钟；第二步，转台开始转动，当按照设定角度步进值转动一定角度后，转台停止转动；第三步，TCSPC按照设定的采集时间采集回波光子数据；第四步，判断转动角度是否达到设定转动范围，若达到，则扫描结束，关闭激光器、TCSPC和转台，若未达到，则跳转至第二步，继续扫描。

为了精确提取距离信息，系统采用多距离图像重建算法^[19]。对于某一像素，在连续时间域上的光子探测率函数λ(t)可以表达为：

式中：*表示卷积；h(t)为目标的响应函数；s(t)为探测器对发射脉冲的归一化响应函数；η为量子效率；b为噪声常数，是背景辐射噪声和暗计数的总和。

对于真实场景来说，若在一个脉冲周期T内接收到K个不同距离的回波信号，则在光子计数直方图中将出现多个峰值^[20]。该像素点接收的目标的脉冲响应函数h(t)可以表达为：

式中：δ(·)为冲激函数；r_i和d_i分别为第i个反射体的反射率和距离；t为时间。

当探测器接收信号之后，以一定的时间分辨率对接收到的信号离散采样，经过N个脉冲周期的积累，光子计数直方图中第k个时间单元的光子数y_k服从泊松分布，即

式中：x表示连续函数h(t)采样后的离散形式，是一个n×1向量；Sx表示s(t)离散采样后与x的卷积，S为一个n×n矩阵；b为常数；I_n为常数1组成的n×1向量。

目标图像重建问题即从光子直方图数据y求解目标响应x的问题。为了方便计算，可以将y的概率分布函数写成关于x的负对数似然函数L(x:y,S, b)。因此，多距离图像重建问题可以转化成在x的非零元素等于K的条件下，最小化负对数似然函数的问题：

式中：τ > 0；L(x:y,S,b)和l₁范数为关于x的凸函数，所以，最小化负对数似然函数也是一个凸优化问题。公式（6）中，各输入变量如下：

(1) y为测量目标物获得的光子直方图数据。

(2) S为先通过测量获得探测器对发射脉冲的响应直方图，再根据直方图拟合成归一化响应函数s(t)，最终建立的S矩阵。文中拟合的归一化响应函数为：

式中：μ为发射激光脉宽。

(3) b为根据y计算的噪声常数。对于全天时、大幅宽的单光子成像场景，背景噪声会随着扫描时长的增加而变化。定义系统盲区距离L_b内的光子直方图数据矩阵y_b为：

在评价噪声水平时，需要计算系统盲区接收光子数的平均值E_yb，用于表征噪声均值；计算系统盲区接收光子数方差σ_yb，用于表征噪声波动情况。实际上，由于接收噪声光子数也服从泊松分布，E_yb近似等于σ_yb，为了节约计算资源，可以只计算E_yb。根据噪声均值的高低，文中对噪声常数b采取分段赋值的方法：

(4) τ = 3b；

利用增加非零限制的迭代收缩阈值算法^[19]，可以求解该问题。文中对单个像素最多提取2个距离信息，即K = 2。对所有M×N个像素运行迭代收缩阈值算法，若存在目标，则将最大响应幅度存入原始幅度像矩阵E₁，对应的目标距离存入原始距离像矩阵L₁；若存在2个目标，则将次级响应幅度存入原始幅度像矩阵E₂，对应的目标距离存入原始距离像矩阵L₂。

在获得原始距离像L₁、L₂和原始幅度像E₁、E₂后，文中根据相邻像素的关联性进一步剔除图像噪声，算法流程如下：

(1) 将所有E₁中的非零像素按照幅值从小到大排列，定义排在前1/q的像素的最大幅值为噪声幅度阈值；将E₁中所有小于噪声幅度阈值的像素赋值为0；

(2) 以E₁中任一非零像素为中心，若至少有7个相邻像素为0，则将该像素赋值为0，

(3) 若E₁中某个像素的值为0，则将E₂、L₁、L₂中对应像素赋值为0；

(4) 输出优化后的L_opt1、L_opt2、E_opt1和E_opt2。

根据优化后的距离像L_opt1、L_opt2和转台扫描的相对角度信息，即可重建激光三维点云。

实验采用逐行扫描模式，某一行水平扫描结束后，俯仰方向步进0.002°，水平方向回到初始角度，开始下一行扫描。扫描按照从左到右，从上到下的顺序进行，每个像素的采集时间为100 ms。整个过程由角度传感器记录实时角度信息。

实验被测目标为2.1 km以外某山顶天线塔和3.1 km以外某大楼。扫描天线塔当日天气为晴天，能见度8 km；扫描大楼当日天气为晴天，能见度14 km。实验装置安装在实验室内。

图6给出了天线塔塔底部分的幅度像成像结果。图6（a）为可见光相机拍摄的目标照片。作为对比，首先将噪声常数b设为噪声均值的固定倍数，即b = 2.5E_yb，原始幅度像重建结果见图6（b）。分析每个像素的光子直方图数据可得，从第一行像素至最后一行，噪声均值从0.28逐渐降低至0.09。这是因为随着扫描时间的增加，太阳逐渐落山，导致背景辐射强度逐渐降低。图6（b）中，在绝大部分没有目标存在的天空区域，噪声依然作为目标被提取了出来，在空间上呈连续分布，塔身的爬架几乎湮没在噪声中，图像质量较差。图6（c）为采用文中的噪声自适应方法获得的原始幅度像，天空区域的噪声大幅下降，目标和噪声的幅度对比度增大。进一步剔除图像噪声，设置幅度阈值为1，大部分噪点被剔除，但尚有少数离散的噪点存在。根据相邻像素的关联性，可以将离散分布的噪点剔除，最终获得的幅度像如图6（d）所示，可以较为清晰地辨认包括攀爬架在内的天线塔结构，噪点抑制效果明显。

Figure 6.  Intensity image of the bottom of the poly. (a) Photo of the target taken by a visible-band camera; (b) Noise threshold b = 2.5 E_yb; (c) Self-adaptive noise threshold; (d) Noise censoring results

消除图像噪点后的天线塔重建图像见图7，天线塔框架结构清晰，背景噪声抑制效果较好。该实验证明了文中提出的噪声常数自适应的图像重建算法可以有效提取目标特征，抑制背景噪声。

Figure 7.  Image reconstruction of the poly. (a) Intensity image; (b) Range image; (c) 3D point cloud reconstruction

为了进一步验证多距离图像的重建效果，笔者对存在前后遮挡关系的目标进行了成像实验。如图8所示，视场范围内1.1 km处的电线对3.1 km处的某建筑存在部分遮挡。

Figure 8.  Photo of the target scene taken by a visible-band camera

重建结果见图9，E₁为提取最大幅值信号重建的图像，对应的目标像为L₁；E₂为提取次级幅值信号重建的图像，对应的距离像为L₂。点云重建结果见图10。从重建结果可以看出，位于前景的电线和被电线遮挡的建筑物能够同时被提取出来，实验结果验证了文中的图像重建算法可以重建多距离目标。

Figure 9.  Multi-range image reconstruction results. (a) Intensity image E₁; (b) Intensity image E₂; (c) Depth image L₁; (d) Depth Estimation L₂

进一步对建筑物进行5次重复测距，每次累计时间均为100 ms。5次测距结果为3 135.05 m、3 134.82 m、 3134.90 m、3 134.98 m和3 135.12 m。根据均方根误差公式

式中：n为测量次数；R_i为每次距离测量值；$\overline R$为测量平均值。计算可得5次测量的均方根误差约为0.11 m。因此，文中系统的测距精度为0.11 m。

Figure 10.  3D point cloud reconstruction. (a) Multi-range object reconstruction; (b) Enlarged image of the building 3.1 km away; (c) Enlarged image of the cables 1.1 km away

文中提出了一种1064 nm全光纤光子计数激光雷达系统，介绍了系统工作原理和硬件构成。然后研究了光子计数图像重建算法，给出了自适应噪声阈值的计算方法。实验结果证明：系统在白天成功重建了3.1 km以外的目标，空间分辨率为0.11 m，超越衍射极限，成功提取多距离目标信息，测距精度为0.11 m，最终重建的图像特征清晰，噪声抑制效果较好。该系统适用于地面固定式或车载式激光雷达，用于远距离目标探测。

但是，文中搭建的系统尚存在扫描成像时间较长、视场较小、探测距离较短等问题。在接下来的工作中，笔者将继续研究基于光纤阵列的光子计数激光雷达系统，进一步缩短成像时间、扩大静态视场范围和提高探测距离，推动该项技术在机载测绘领域的发展。