激光雷达已在三维成像、高精度对地观测和深空探测等空间领域得到诸多应用，例如，NASA已经成功发射的火星轨道激光高度计（MOLA）、近地小行星激光高度计（NLR）、水星激光高度计（MLA）、地球科学激光高度计（GLAS）和月球轨道激光高度计（LOLA）等^[1-3]多个星载激光雷达系统，用于地球或行星表面的高程测量。我国也先后发射了嫦娥一号激光高度计、ZY-3系列的星载激光测距仪等星载激光雷达系统^[4-6]。上述激光雷达采用了大能量低重频激光发射技术和回波线性探测体制，这种探测方式需要发射大能量脉冲以获得足够的回波能量，保证高信噪比，进而确保其探测概率和测距精度。但是该类型的激光雷达难以实现对目标的多波束、高重频探测，空间分辨率较低。

基于上述问题，高重频（10 kHz上下）微脉冲（百μJ左右）光子计数激光雷达逐渐成为新型星载激光雷达的发展方向^[7-8]。高重频微脉冲激光发射技术显著降低了单个激光脉冲能量，从百mJ降低至百μJ水平，减小了对激光器本身和光学系统造成损坏的风险，降低了对卫星平台的资源需求，可实现多波束激光高重频发射，大大提高卫星对地采样密度和覆盖宽度。星载微脉冲光子计数激光雷达工作时，向地表发射高重频微脉冲激光，单光子探测器采用G-APD、PMT等特殊模式的高灵敏度探测器件和电路技术，实现对地表反射的单光子回波信号进行响应，产生回波光子计数点，显著降低了对回波脉冲能量的要求，同时单光子探测器也对各种背景噪声响应产生噪声光子计数点，所有的光子计数点累计形成点云。另外，时间相关光子计数（TCSPC）技术是一种基于统计学的微弱信号处理技术，可从时域上包含随机噪声的多脉冲累计的光子点云中提取与时间（距离）相关的有效信号，将真实测距数据从点云中提取出来^[9]。美国的先进地形激光测高系统（ATLAS）即采用了高重频微脉冲激光发射技术，激光器重频为10 kHz，波长为532 nm，同时采用具有单光子探测能力的多像元光电倍增管（Photomultiplier Tube, PMT）作为探测器以实现微脉冲光子计数体制^[10]。

微脉冲光子计数激光雷达系统的探测性能受到回波光子数、背景光强度、单光子探测器时间抖动及探测器死时间等多种因素的影响^[11-13]。与盖革模式雪崩光电二极管(APD)相比，光电倍增管（PMT）具有单光子探测能力，死时间更短，可以到3 ns，光敏面更大且暗噪声相对更小，有成熟的商用多阳极集成器件。但是，单像元光电倍增管仅能对死时间内到达的首个光子响应，而对死时间内的其他光子则无法响应，从而造成飞行时间测量偏短，这种现象称为首光子误差^[14]。此外，单像元光电倍增管受限于死时间，总的回波光子计数减少，测距精度不高。而对于多像元PMT探测器，每个像元可以独立计数，进一步减小等效死时间，降低了首光子误差，提高探测性能^[15]。

文中依据星载微脉冲光子计数激光雷达的工作原理，建立了基于多像元光电倍增管（PMT）的探测仿真分析模型，并对典型探测过程进行仿真分析，验证仿真模型的正确性。

激光雷达系统仿真模型主要包括发射机模型和接收机模型，地物模型采用简单的朗伯反射模型^[16]。暂不考虑大气透过率、光学系统效率、背景噪声（太阳背景噪声及器件暗电流噪声）、探测器死时间等因素的影响^[17]。

星载光子计数激光雷达采用高重频微脉冲激光器，激光脉冲存在着一定时间和空间的能量分布。文中建立激光发射仿真模型时，激光脉冲功率的时间分布和空间分布均采用高斯脉冲模型近似表示^[18]。激光发射模型主要参数为波长532 nm，发射重频10 kHz，脉冲能量E_t为40 μJ，激光脉冲半高宽τ为5 ns，发散角$\alpha $为0.1 mrad。

激光脉冲功率的时域分布曲线如图1所示，其数学模型如公式(1)所示：

式中：E_t为单脉冲能量；τ为激光脉冲的半高宽。

同时，激光发射脉冲能量的空间分布也是不均匀的，星载固体激光器的空间分布可近似为高斯脉冲，如图2所示，其数学表达式如公式(2)所示：

式中：x和y表示光斑内任意点相对光斑中心的位置；α_x和α_y为激光发散角；G_x和G_y为超高斯因子，用于描述能量分布形状。

图  1  激光脉冲的时域分布示意图

Figure 1.  Time domain distribution of laser pulse

图  2  激光脉冲的空间分布示意图

Figure 2.  Spatial distribution of laser pulse

仿真分析时，激光脉冲能量的空间分布采用$N$个正方形单元离散表示，如图3所示，图中外围的圆表示激光光斑，正方形为离散单元，每个离散单元的能量分布可根据离散单元几何中心坐标近似获得，如公式（3）和公式（4）所示：

式中：x_i和y_i表示光斑内离散单元几何中心点相对光斑中心的位置；N为离散单元数量；α为激光发散角，一般取α_x=α_y=α。

综上可得每个空间离散单元的激光发射模型如公式（5）所示：

图  3  激光脉冲的空间离散分布示意图

Figure 3.  Spatial discrete distribution of laser pulse

根据激光雷达作用距离方程可计算激光脉冲照射至朗伯体上后，某个空间离散单元的激光脉冲反射回的瞬时平均回波光子数为${S_{pe}}$，如公式（6）所示：

式中：${S_{pe}}$为瞬时平均回波光子数；${P_{lt}}$为发射激光脉冲的单个空间离散单元能量；$\lambda $为激光波长；$h$和$c$分别为普朗克常量和光速；$\rho $为反射率；$\theta $为发射光轴与地面目标法线的夹角；$FOV$为接收光学系统的视场角；$\alpha $为出射激光发散角；${A_r}$为接收面积；$R$为目标距离；${\eta _t}$为激光发射系统透过率；${\eta _r}$为激光接收系统透过率；${\eta _a}$为大气透过率。

将一个激光发射脉冲周期在时域上离散为n个大小相同的$\Delta t$时间间隔，则在单个时间区间内单个空间离散单元的平均回波光子总数为：

探测器件的探测过程主要取决于两个随机过程，一个是光电子的产生过程，一个是倍增激发的响应过程。若假设探测器量子效率为${\eta _q}$，不考虑噪声影响，则在单个时间区间内探测器产生的平均回波光子数${S_{pe}}$可近似采用公式（8）表示：

对于星载光子计数激光雷达，回波光子到达PMT探测器并产生光电子的概率可使用泊松分布近似描述^[19]。由以上分析可知，若在某个时间区间$\Delta t$内的平均回波光子总数为$ {S}_{p}(\Delta t)$，不考虑噪声影响，根据泊松分布可知，在该时间区间内产生的${S_{pk}}$个光电子的概率如下：

一般情况下，由于单个空间离散单元的回波光子数较少，笔者认为单个空间离散单元在时间区间$\Delta t$内，探测器只存在被倍增激发一次或者不激发的状态，不存在两次或多次响应的情况。则在单个时间区间内不产生光电子的概率为：

则在单个时间区间内至少产生1个光电子的概率为：

以上分析是探测过程中单个时间区间内回波光子转化为光电子的过程，下面分析在n个$\Delta t$时间区间中第j个时间区间（$1 \leqslant j \leqslant n$）发生倍增激发响应的概率。

若探测器在第j个时间区间内发生倍增激发响应，则在其之前的j−1个时间区间内必然无光电子激发且在第j个时间区间内至少产生一次倍增激发响应过程，由上可得，探测器在第j个时间区间内产生倍增激发响应的概率为：

式中：$\left[ {\displaystyle\prod\limits_{i = 1}^{j - 1} {{p_i}({S_{pki}} = 0)} } \right]$表示在j−1个时间区间内无光电子激发的概率；${p_j}({S_{pkj}} \geqslant 1)$表示在第j个时间区间内至少产生一次倍增激发响应的概率。

每产生一次倍增激发响应便获得一个光子计数，将单个离散单元在n个时间区间$\Delta t$内的探测器激发的次数依次累计叠加，即可获得该单个离散单元的光子计数点云。以此类推，将一个完整发射脉冲视为由N束小激光构成且每束小激光的功率不同。依次将N个空间离散单元的激光功率带入上述公式，重复计算即可累积获得一个完成脉冲的光子计数点云。

为了简化多像元光电倍增管的像元列阵问题，将每个激光脉冲离散单元的回波光子数均分到各个独立像元上，假设共有m个像元，则单个像元上的回波光子数如公式（13）所示。各像元分别对各自的回波光子产生光电子、发生倍增激发响应并产生光子计数，最后将各个像元的光子计数点云叠加后获得单个离散单元的完整的光子点云数据，将所有激光脉冲离散单元的所有回波光子点云数据累加，即为整个回波光子点云数据。

式中：m为像元数；$ {S}_{pm}(\Delta t)$为$\Delta t$时间区间内到达每个像元的平均光子总数。

此外，PMT探测器的工作原理导致其存在时间抖动特性，即探测器的输出电流脉冲与实际响应时间存在一定误差，误差可近似为泊松分布^[15]，文中仿真时，考虑了PMT探测器输出计数脉冲的时间抖动，每当获得一个光子计数并输出响应时，就增加一次探测器时间抖动，并获得最终目标时间。

首先根据激光发射模型计算获得单个激光发射脉冲的每个空间离散单元的脉冲功率${P_{lt}}$，根据公式（6）、公式（7）和公式（13）获得单个时间区间$\Delta t$内到达每个像元的平均光子总数$ {S}_{pm}(\Delta t)$。带入量子效率${\eta _q}$计算获得在单个时间区间内探测器单像元产生的平均光电子数${S_{pe}}$。

根据探测概率模型计算得到的每个时间区间内的光电子产生概率$p({S_{pk}} \geqslant 1)$及其在第j个时间区间内产生倍增激发响应的概率${p_j}$，建立上述两个随机过程的泊松分布发生器产生随机数来决定是否在第j个时间区间内完成光电探测。若随机数大于0，则表示有光子触发，计算光子飞行时间，并在该次飞行时间上附加PMT探测器随机的时间抖动，即获取一个光子计数。

在同一时间周期内，首先根据像元数累积不同像元之间的光子计数，然后再累积不同离散单元的光子计数，所有的单元累积完成后，单个脉冲的的光子计数点云计算完成。重复多个单脉冲的探测仿真过程，即可得到多脉冲累计的光子计数点云的统计结果。上述仿真过程均基于MATLAB软件编程实现（图4）。

为了验证多像元光电倍增管的探测优势，针对平坦地形，分别仿真分析了不同像元数下50000个脉冲的回波统计特性，累计光子计数曲线如图5所示。对每个脉冲，统一设置其平均回波光子数均为16个，探测器时间抖动均为100 ps，死时间为5 ns。从图5中可以明显看出，随着PMT像元数的增多，系统获取的光子计数统计波形越接近理想回波，飞行时间越靠近真实值。对不同像元数下50000个脉冲的回波光子计数结果进行统计分析，获得的光子飞行时间均值如图6所示。由图6可知，随着像元数的增加，光子飞行时间均值逐渐趋于平稳，更加迫近真实值，首光子误差明显减小，相对于单像元探测，16像元探测系统的首光子误差减小0.55 ns。

图  4  仿真算法流程示意图

Figure 4.  Flow diagram of the simulation algorithm

图  5  不同PMT像元数下的回波光子统计特性

Figure 5.  Statistical characteristic of echo photons under different numbers of PMT pixel

图  6  不同PMT像元数下的多脉冲飞行时间均值 (参考飞行时间为100 ns)

Figure 6.  Mean of multi-pulse TOF under different numbers of PMT pixel (The reference TOF is 100 ns)

文中仿真分析了不同地表坡度下10000个脉冲的回波统计特性。累计光子计数曲线如图7所示。对每个脉冲，统一设置其平均回波光子数为16，死时间为5 ns。从图7可知，随着地表坡度增大，回波激光脉宽明显展宽，系统所获取的光子计数统计波形越平坦，导致光子事件的不确定时域越大，不利于精确提取准确的光子飞行时间。对不同地表坡度下光子飞行时间标准差进行统计分析。光子飞行时间标准差定义如公式（14）所示，其代表了光子计数点云数据的离散程度或抖动程度，直接影响激光雷达的测量精度。由图8可知，随着地表坡度增大，光子飞行时间的标准差明显增大，当地表坡度增大到10°，光子飞行时间的标准差增大到5.33 ns，严重影响了激光雷达的测量精度。

图  7  不同坡度角的回波光子统计特性

Figure 7.  Statistical characteristic of echo photons with different slope angles

图  8  不同坡度角的飞行时间标准差

Figure 8.  Standard deviation of TOF with different slope angles

式中：$\sigma $为飞行时间标准差；$N$为计数光子总数；${t_i}$第i个计数光子的飞行时间；$\bar t$为N个计数光子飞行时间的平均值。

上述仿真分析中，均针对单一地形，采用脉冲多次累计的方法分析回波光子计数的统计特性。而卫星在轨飞行时，激光光斑随着地面起伏波动形成一定间距的足印轨迹，无法实现对地面单一目标的大数量级累计脉冲测量。因此，有必要仿真分析在复杂地形条件下的回波光子计数点云的统计特性。

针对同一复杂地形，统一设置平均回波光子数为8，死时间为5 ns，不同像元数对回波光子计数点云的影响如图9所示。图9中横坐标为沿轨飞行距离，纵坐标为地表高程，红色曲线为沿飞行轨迹的真实地表高程曲线，蓝色点云为仿真获取的光子计数点云。由图9可知，提高像元数能够明显增加回波光子计数点的数量，大量光子计数点分布在实际高程曲线附近，更能反映真实的地表高程，提高针对复杂地形的测距精度。

图  9  复杂地形下的不同PMT像元数对平均回波光子点云的影响

Figure 9.  Effect of different numbers of PMT pixel on the average echo photons point cloud under the complicated terrain

图10为针对同一复杂地形，统一设置为4 pixelPMT探测、死时间均为5 ns时，不同平均回波光子数下获得的回波光子计数点云。从图中明显看出，随着平均回波光子数的增加，回波光子计数明显增多，并且分布更集中，与实际地表高程曲线吻合度越高。因此，提高平均回波光子数能够明显提高微脉冲光子计数激光雷达的测距精度。

图  10  复杂地形下不同平均回波光子数对回波光子点云的影响

Figure 10.  Effect of different number of average echo photons on the echo photons point cloud under the complicated terrain

由于星载微脉冲光子计数激光雷达仍处于前期论证研究阶段，难以开展试验验证，笔者通过光子计数激光雷达的机载飞行试验初步验证微脉冲光子计数激光雷达的探测性能，并完成光子计数激光雷达空地全链路的技术验证，为后续星载微脉冲光子计数激光雷达的仿真设计优化积累一定的试验基础。

机载飞行试验的主要技术指标如下：测距范围为(3.85±0.5) km，激光中心波长532.3 nm，脉宽1 ns，发射重频10 kHz，激光发散角0.1 mrad，接收视场角0.2 mrad。机载试验采用运五飞机，试验地点为东营市黄河口镇，图11为机载飞行试验的部分测试路径，箭头方向为飞行方向，由图可知，该路径上的典型特征即为厂房与油罐。

图  11  机载飞行试验的某段飞行路径

Figure 11.  Flight path of the airborne flight test

通过Matlab程序解算机载飞行试验路径上油罐和厂房这些典型地面特征的回波光子点云数据，具体如图12所示。由图可知，在相同的机载试验条件下，随着PMT像元数的增加，有效回波光子计数点明显增加，油罐顶部的轮廓特征更加清晰准确。这与复杂地形下的光子计数点云仿真分析结果一致，验证了仿真分析结果的正确性。但是在机载飞行试验中，随着PMT像元数的增加，背景噪声明显增大，这可以通过噪声消除算法消除噪声的影响，而在仿真分析中尚未考虑背景噪声的影响。仿真分析和试验结果说明，微脉冲光子技术激光雷达能够实现复杂地形下的高精度光子计数测距。

图  12  典型地貌下不同PMT像元数对回波光子点云的影响

Figure 12.  Effects of different numbers of PMT pixel on the echo photons point cloud under the typical landforms

在实验室条件下，通过固定光纤长度对机载样机进行测距精度标定，不同信噪比下其测距精度为0.12~0.20 m。对机载飞行试验数据，首先采用基于密度的聚类算法对点云数据进行自适应去噪，然后以较为平坦的屋脊线作为参照物，通过最小二乘法拟合获取9段理论屋脊线，并分别计算每段屋脊的数据点云的均方根误差为0.13~0.27 m，所有数据点总误差在0.18 m，与试验室标定精度接近。

微脉冲光子计数激光雷达已成为星载激光雷达的重要发展方向。文中依据微脉冲光子计数激光雷达的工作原理，建立了仿真模型并对其探测过程进行仿真，分析了特定条件下平坦地形、不同坡度地形和复杂地形的回波光子计数点云的统计特性，验证了探测模型的正确性。仿真结果表明，增加PMT探测器的像元数能够有效降低激光雷达系统的首光子误差，16 pixel PMT探测系统的首光子误差能够减小0.55 ns；地形坡度增大会导致激光雷达系统的测距标准差明显增大，故为保证测距精度，光子计数激光雷达测量的地表坡度最好不超过5°；对于复杂地形，增加PMT探测器像元数以及平均回波光子数均能够获取更接近地表轮廓的光子计数点云，提高对复杂地形的测距精度。另外，机载飞行试验验证了多像元PMT能够显著增加有效回波光子计数点，准确有效地反映了被测地形的轮廓特征，实现复杂地形下的高精度光子计数测距，这验证了仿真分析结果的正确性。下一步将继续优化仿真参数，对大气环境、复杂目标特性、背景光等其他各种因素对星载微脉冲光子计数激光雷达性能的影响进行深入研究分析。