Calibration and correction of afterpulse characteristics of a photomultiplier tube

近年来，单光子探测技术已经广泛地应用于通信、测距以及生物荧光检测等领域^[1-5]。光电倍增管^[6](photomultiplier tube，PMT)是一种常用的单光子探测器，它能提供快速响应和高增益，被广泛应用于机载海洋激光雷达系统中^[7-9]。在使用过程中光电倍增管不可避免地会受到海面强激光的照射，因此无法完全消除光电倍增管的后脉冲。后脉冲形成的主要原因是PMT内的残留气体与光阴极或电子倍增极产生的光电子之间相互碰撞而发生电离。后脉冲的存在会破坏高强度输入尾部低强度信号的测量，信号衰减的基线加到真实的激光雷达回波信号中，导致信号衰减比预期的更长，这种现象会使得测量误差增大。因此，在激光雷达系统中必须对其使用的PMT后脉冲特性进行检测，以便获得准确的回波信号。

国内外的研究人员对PMT后脉冲对激光雷达数据的影响以及消除办法做了一些研究^[10-13]。其中，CAMPBELL等人^[12]首次提出了一种校正PMT后脉冲效应的方法，即推导后脉冲概率法，并通过对闪烁星光的测量，验证了该校正方法的有效性。但由于实验中采用的两个用来对比的光电倍增管没有达到预期的效果，因此只能说明该方法在提供光电倍增管后脉冲消除办法方面有一定潜力，还需进一步研究。WILLIAMSON等人^[13]则提出在PMT的光电阴极外部放置金属环来降低PMT的信号噪声，但这种方法需要在不影响激光雷达衰减常数的情况下确定金属环在接收机中的最佳位置，因此在实际操作中有一定的困难。

本文中在研究光电倍增管的结构和标定原理的基础上，分别对波长为486nm和532nm的激光雷达系统进行了后脉冲的标定。通过实验得到了在强光入射下，后脉冲计数对后续入射光计数的影响，即后脉冲概率分布函数，实现了对机载海洋激光雷达飞行数据的校正，得到了较为准确的海洋探测深度。最后通过蒙特卡洛仿真模拟激光在水中的传输过程，将实验结果与模拟结果对比, 进一步说明了标定方法的可靠性和优势^[14]。

当激光照射目标时，激光雷达利用目标的反射、折射、散射和透射等产生的回波辐射进行探测^[15]。图 1所示是一个典型的激光雷达系统结构示意图。激光雷达主要由发射、接收和信号处理3个部分组成，其中接收部分主要包括光学接收天线和光电探测器，光电探测器的作用是将激光信号直接转变为电信号。

Figure 1.  Schematic diagram of lidar system structure

光电倍增管是一种真空玻璃管，其工作原理如图 2所示。它一般包括3个部分：光电阴极、一个阳极和若干个二次发射极。光电倍增管的工作原理是光子撞击光电阴极发生光电效应，从而产生光电子，产生的光电子在电压的偏转作用下被加速，聚焦到了二次发射极，遭到碰撞后的二次发射极会释放出更多的二次电子，以此类推，经过一系列的碰撞，使得电子数加倍，直到最后阳极收集电子，形成阳极电压或电流^[16]。

Figure 2.  Schematic diagram of the photomultiplier tube

光电倍增管在闪烁计数或激光脉冲检测中的脉冲检测模式下工作时，可以观察到一些后脉冲。后脉冲常常干扰大幅度脉冲后低电平信号的精确测量，降低了闪烁计数中的能量分辨率，并在脉冲计数应用中造成误差。因此，必须对激光雷达系统中使用的光电倍增管的后脉冲特性进行检测，以便校正其输出信号。

连续光源或者脉冲光源在经过适当的光衰减后进入光电倍增管，通过分析光电倍增管输出脉冲的统计分布来测算其后脉冲分布的概率。在分析时，需要对大量采样时间间隔内计数分布进行统计，从而计算出后脉冲的分布概率，因此处理时将数据进行了20000次累加。分析过程中遵循的基本理论原则如下：(1)光电倍增管的入射光源和噪声的分布均满足泊松分布统计特性；(2)每个事件均可以以一定概率引发一个或多个后脉冲，且后脉冲产生的概率与入射光强度相互独立。

当光电倍增管检测以单位时间间隔速率为r出现的光子时，在没有后脉冲的情况下，从0开始每间隔Δt计数，单位时间内光子计数为0, 1, 2, 3…的概率分别为: $\mathrm{e}^{-r \Delta t}, r \Delta t \mathrm{e}^{-r \Delta t}, \frac{1}{2}(r \Delta t)^{2} \mathrm{e}^{-r \Delta t}, \frac{1}{6}(r \Delta t)^{3} \mathrm{e}^{-r \Delta t} \cdots$。

在考虑后脉冲的情况下，由于每个后脉冲的出现需要一个初始事件，因此单位时间间隔Δt内零计数的概率p(0)不因后脉冲而改变，其概率函数可由下式得出^[12]：

如果用p_any表示特定事件后任意数量的后脉冲的概率，那么单位时间间隔内光子数为1的概率值p(1)为：

需要特别注意的是，在计算单位时间间隔内光子数为2的概率值p(2)时，需要在原始泊松分布概率值的基础上减去任一计数产生的后脉冲概率值，同时增加光子计数为1时产生的后脉冲概率值，其表达式为：

式中，p₁表示仅在事件后出现一个后脉冲的概率。

通过对(1)式~(3)式的处理，可以推导出光子的间隔速率r和特定数量的后脉冲的概率：

假设p为单个初始后脉冲的概率，那么p_any由下式给出：

如果一个初始后脉冲产生单个后续后脉冲的概率值为p_s，那么它产生一个或多个后续后脉冲的概率值p_s′为：

将方程(3)改写成：

从而可得：

为了对光电倍增管的后脉冲效应做出标定，搭建了如图 3所示的测试系统。通过该系统对入射待测PMT的光源进行了必要的调制和适当的衰减，实验结果用高速采集卡接收，其采集到的数据用计算机进行处理。测试系统主要包括激光光源、示波器、声光调制器(acoustic optical modulator，AOM)、数字可变衰减器、倍频晶体(periodically poled lithium niobate，PPLN)、待测光电倍增管等。波长为1064nm的连续光源通过光纤导入到额定电压为12V的声光调制器中，被调制为脉冲激光，为满足实验要求再通过倍频晶体改变其波长大小到532nm，将出射532nm激光的光纤固定在支架上，并将其对准待测PMT进行测试，测试系统主要参量如表 1所示。

Figure 3.  Physical diagram and composition diagram of the test system

对PMT的后脉冲分布进行实验室标定，并通过手持衰减器对光强进行不同程度的衰减，读取高速采集卡采集到的数据，然后对这些数据进行处理，具体步骤见下。

(1) PMT本底噪声修正。由于探测器的缓慢放电效应，探测器的基线存在随采样时间缓慢降低的情况。为了得到基线值的拟合函数，在探测器离散回波信号出现之后，读取不同时刻下探测系统的基线值并作出散点图，然后利用Curve Fitting工具箱对散点图进行曲线拟合，如图 4所示。

Figure 4.  Curve-fit result for the baseline of 532nm and 486nm channels

根据曲线拟合结果，532nm和486nm基线幅值与采样时间t的关系如下:

在进行数据处理之前，需要将实验采集的数据减去基线幅值的变化，以便得到更为准确的接收光子数。

(2) 背景光噪声修正。在使用光电倍增管接收光子时，除了目标反射回来的激光回波信号，还掺杂着一定程度的背景光噪声。为获得背景光噪声的产生概率，需要对所采集到的光子数据进行累加观察其拖尾现象，累加结果如图 5所示。其中图 5a的结果表明，在光强衰减值为0dB时，光电倍增管的拖尾效应持续的时长约为1000ns，并且通过纵坐标可以大致看出经过20000次累加之后，系统接收到的背景光的数目为6个，因此相应采集到的背景光噪声概率p_n为：

Figure 5.  Result of accumulating the photons collected by the photomultiplier tube was obtained under the noise

为进一步确认背景光噪声的概率值，在光强衰减值分别为40dB，50dB，60dB时对收集到的数据又进行了相同的处理，其中60dB累加结果如图 5b所示。对比两个小图可知，背景光噪声概率为3/10000。

利用MATLAB编写背景光噪声过滤程序，具体程序流程如图 6所示。首先提取需处理的实验室数据，并将数据的行数和列数分别命名为a，b；然后引入随机变量i，j以及背景噪声光子数N_p, 当满足条件i≤a且j≤b时，系统借助随机函数rand产生[0, 1]区间内均匀分布的随机数，并将随机数大小与背景光噪声概率3/10000相比较，统计所有小于3/10000的值的个数并赋值给N_p；最后将原始数据在经过PMT本底噪声修正的基础上减掉背景光噪声的光子数。

Figure 6.  Program flow chart for the background noise of PMT

(3) 累加降噪后的结果。对不同光衰减强度下的光电倍增管接收到的光子数进行20000次累加，得到结果如图 7所示(以532nm激光为例)。

Figure 7.  The number of photons received by PMT after noise reduction is accumulated

由统计结果可知，随着入射光强的不断衰减，PMT接收到的光子数逐渐减少，同时后脉冲效应也逐渐减弱。在入射光强衰减值大于30dB后，后脉冲效应将不再明显，并且入射光下的后脉冲分布已经足够低。因此，对于入射光强低于30dB的情况，将按照30dB入射光下的后脉冲分布函数进行插值。

(4) 后脉冲概率分布拟合函数。对光电倍增管的后脉冲部分进行拟合，得到探测器后脉冲分布概率随采样时间之间的变化关系式。单光子探测器的后脉冲概率分布函数可以用统一的双指数函数进行描述^[15]。

本文中借助的是MATLAB拟合工具箱Curve Fitting，拟合过程分为以下几个步骤：首先导入不同光衰减下需要进行拟合的数据；然后选择拟合方程为双指数函数型；最后得到拟合方程曲线以及拟合函数，结果如图 8所示。

Figure 8.  The fitting image of PMT and its probability distribution function under different optical attenuation values

图 8所得拟合曲线对应的探测器后脉冲函数分别为：F(x)=1597e^-0.03378x+82.01e^-0.003129x(见图 8a)；F(x)=1597e^-0.03697x+82.01e^-0.003427x(见图 8b)；F(x)=310.2e^-0.02931x+11.28e^-0.001428x(见图 8c)；F(x)=101.2e^-0.02491x+82.01e^-0.00255x(见图 8d)；F(x)=44.33e^-0.07589x+82.01e^-0.004463x(见图 8e)；F(x)=10.31e^-0.04199x+82.01e^-0.000504x(见图 8f)；F(x)=0.494e^-0.005483x+3.958e^-0.02201x(见图 8g)。由图 8可知，随着光衰减值越来越大，探测器在后脉冲区间内接收到的光子数从1600减少到8，这说明探测器后脉冲效应与激光强度成正相关。

为了获得任意光强度下探测器的后脉冲分布，需要对位于相邻两个入射光强度下的待校正的后脉冲分布函数采用线性插值的方法求解，插值公式^[17]如下：

式中，F_ap(j; i)表示第i ns所接收到的入射光强所产生的后脉冲在第j ns的分布，N_ap(i) 表示在第i ns所接收到的后脉冲光子数，且N_ap(i)∈[N_c，a, N_c，b], N_c，a，N_c，b为实验室标定的在特定衰减值a，b下系统所采集到的入射光脉冲光子数，F_c，a，F_c，b为入射光子数是N_c，a，N_c，b时的后脉冲概率分布函数。

这一节中通过对实验系统所得数据的处理，得到了光电倍增管后脉冲的分布情况随采样时间变化的函数关系式，并且利用插值公式求解了相邻两个入射光强度下的后脉冲分布概率，最后可通过线性插值法来分段统计其后脉冲光子数。

本次实验中采用中国科学院上海精密机械研究所空间激光信息技术研究中心的研究团队在我国南海水域测得的数据，测试中所用的激光雷达为自主研制的486nm和532nm蓝绿波长多通道海洋激光雷达，飞行所得原始数据如图 9所示。

Figure 9.  Flight data of airborne marine lidar

对机载海洋激光雷达飞行数据的校正流程可以概括为以下3个步骤：(1)去除背景光噪声；(2)结合第3.1节中所得PMT后脉冲概率分布函数与(15) 式求出位于不同衰减值区间的后脉冲光子数；(3)将原始飞行数据减去背景光噪声和后脉冲光子数，即可得到探测器接收到的真实光子数。

通过实验室标定系统得到了PMT的后脉冲概率分布函数，利用后脉冲的概率分布函数对机载海洋激光雷达的机载数据进行后脉冲校正。图 10所示为机载海洋激光雷达测得的数据经过后脉冲校正后的结果。

Figure 10.  Marine flight data processing results of airborne marine lidar

由图 10可知，对比后脉冲校正前后的数据，校正前的数据显示后脉冲拖尾效应较为明显，532nm，486nm激光通道的最大测深深度分别约为110m，200m；经过后脉冲校正后的波形无明显后脉冲拖尾，且532nm，486nm激光通道所探测到的最大水深为77m，102m左右。由此可得，机载雷达数据经过后脉冲校正后去除了假信号，使得测深距离更为准确，并且由图可以看出，由于散射层的存在，两个激光通道的衰减率在40m水深处均有所增加，且随着海水深度的增加，浑浊程度也随之增加，486nm激光的优势减弱，因此水深大于40m后，486nm通道与532nm通道的衰减速率更加接近。

为验证后脉冲分布函数法校正机载数据的正确性，利用蒙特卡洛仿真模拟激光水下传输过程，并将其结果与实验结果进行对比。蒙特卡洛方法是在对大量光子随机运动轨迹计算的基础上，模拟光在浑浊介质传播的整个物理过程^[18-20]。仿真参量与实验情况保持一致，具体参量设置如下：模拟的光子包数为10000个，飞机高度H=300m，接收光学天线面积为0.2m²，海水衰减系数C=0.05m^-1，最大散射次数设为50，接收机视场角为20mrad，光电倍增管相关系数不变, 如表 1所示。图 11中显示了实验结果与仿真结果的对比。

Figure 11.  Comparison of experimental and simulated waveforms

为了更加直观地反映实验波形和模拟波形之间的相关关系，本文中将通过计算两条曲线的相关系数来定量表示其相关程度。具体做法为：分别提取同一波长下的实验波形和模拟波形数据，将同一横坐标下对应的两个点划分为一组，共提取了663组这样的点，将实验波形上提取到的663个点的纵坐标值构成向量y₁，同样的，从模拟波形上提取到的663个点的纵坐标值构成了向量y₂，计算y₁和y₂的相关系数R，最后通过计算得到激光波长为486nm，532nm的实验波形曲线和模拟波形曲线的相关系数分别为0.9689，0.8648。根据参考文献[21]可知，相关系数的绝对值大于0.8表示两个变量高度线性相关，说明实验方法正确地获得了激光雷达的真实信号，即对其可行性进行了验证。

由图 11可知，相比实验波形，模拟波形看起来更加平滑，这是因为模拟状态下噪声信号几乎为零；并且对比同一波长下的实验波形与模拟波形可知，随着探测水域深度的增加，实验波形的衰减速度明显要比模拟波形快一些，且实验测得水深比模拟结果小10m左右，这是由于真实状态下，海水中的悬浮物和不溶物对光信号产生了漫反射，使得激光雷达未能在有效接收面积内接收更多的光子，而模拟中未考虑这些情况；另外，由于蒙特卡洛仿真在模拟激光水下传输时能够到达接收机的光子十分有限，因此需要耗费大量时间来得到一个可靠的统计结果，而后脉冲概率函数法是在原有机载信号的基础上进行校正，故可以很好地避免这个缺点。

针对机载海洋激光雷达在使用中普遍存在的光电倍增管后脉冲的问题，本文中通过对光电倍增管后脉冲的标定，获得了光电倍增管后脉冲概率分布函数，在不搭建仿真模型的前提下正确地校正了机载海洋激光雷达实测数据中的后脉冲区域，并通过与蒙特卡洛仿真结果的对照分析，验证了光电倍增管标定方法的准确性和可靠性。因此，利用后脉冲概率分布函数能够更加快捷、简便地校正光电倍增管的后脉冲效应，节约了大量时间，提高了以光电倍增管为核心器件的机载海洋激光雷达测量数据的精确度，为进一步研究海洋特性奠定了良好的基础。

感谢中国科学院上海光学精密机械研究所贺岩、罗远等人提供的实验设备及机载数据。