王 媛,程小勁
(上海工程技術大學 機械與汽車工程學院,上海 201620)
近年來,單光子探測技術已經廣泛地應用于通信、測距以及生物熒光檢測等領域[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個部分組成,其中接收部分主要包括光學接收天線和光電探測器,光電探測器的作用是將激光信號直接轉變為電信號。

Fig.1 Schematic diagram of lidar system structure
光電倍增管是一種真空玻璃管,其工作原理如圖2所示。它一般包括3個部分:光電陰極、一個陽極和若干個二次發射極。光電倍增管的工作原理是光子撞擊光電陰極發生光電效應,從而產生光電子,產生的光電子在電壓的偏轉作用下被加速,聚焦到了二次發射極,遭到碰撞后的二次發射極會釋放出更多的二次電子,以此類推,經過一系列的碰撞,使得電子數加倍,直到最后陽極收集電子,形成陽極電壓或電流[16]。

Fig.2 Schematic diagram of the photomultiplier tube
光電倍增管在閃爍計數或激光脈沖檢測中的脈沖檢測模式下工作時,可以觀察到一些后脈沖。后脈沖常常干擾大幅度脈沖后低電平信號的精確測量,降低了閃爍計數中的能量分辨率,并在脈沖計數應用中造成誤差。因此,必須對激光雷達系統中使用的光電倍增管的后脈沖特性進行檢測,以便校正其輸出信號。
連續光源或者脈沖光源在經過適當的光衰減后進入光電倍增管,通過分析光電倍增管輸出脈沖的統計分布來測算其后脈沖分布的概率。在分析時,需要對大量采樣時間間隔內計數分布進行統計,從而計算出后脈沖的分布概率,因此處理時將數據進行了20000次累加。分析過程中遵循的基本理論原則如下:(1)光電倍增管的入射光源和噪聲的分布均滿足泊松分布統計特性;(2)每個事件均可以以一定概率引發一個或多個后脈沖,且后脈沖產生的概率與入射光強度相互獨立。

在考慮后脈沖的情況下,由于每個后脈沖的出現需要一個初始事件,因此單位時間間隔Δt內零計數的概率p(0)不因后脈沖而改變,其概率函數可由下式得出[12]:
p(0)=e-rΔt
(1)
如果用pany表示特定事件后任意數量的后脈沖的概率,那么單位時間間隔內光子數為1的概率值p(1)為:
p(1)=e-rΔtrΔt(1-pany)=
p(0)rΔt(1-pany)
(2)
需要特別注意的是,在計算單位時間間隔內光子數為2的概率值p(2)時,需要在原始泊松分布概率值的基礎上減去任一計數產生的后脈沖概率值,同時增加光子計數為1時產生的后脈沖概率值,其表達式為:

(3)
式中,p1表示僅在事件后出現一個后脈沖的概率。
通過對(1)式~(3)式的處理,可以推導出光子的間隔速率r和特定數量的后脈沖的概率:

(4)
pany=1-p(1)/[p(0)rΔt]
(5)

(6)
假設p為單個初始后脈沖的概率,那么pany由下式給出:
pany=p+p2+p3…=p/(1-p)
(7)
如果一個初始后脈沖產生單個后續后脈沖的概率值為ps,那么它產生一個或多個后續后脈沖的概率值ps′為:
ps′=ps+ps2+ps3…=ps/(1-ps)
(8)
ps=ps′/(1+ps′)
(9)
將方程(3)改寫成:
p(0)rΔtp(1-ps′)
(10)
從而可得:

(11)
為了對光電倍增管的后脈沖效應做出標定,搭建了如圖3所示的測試系統。通過該系統對入射待測PMT的光源進行了必要的調制和適當的衰減,實驗結果用高速采集卡接收,其采集到的數據用計算機進行處理。測試系統主要包括激光光源、示波器、聲光調制器(acoustic optical modulator,AOM)、數字可變衰減器、倍頻晶體(periodically poled lithium niobate,PPLN)、待測光電倍增管等。波長為1064nm的連續光源通過光纖導入到額定電壓為12V的聲光調制器中,被調制為脈沖激光,為滿足實驗要求再通過倍頻晶體改變其波長大小到532nm,將出射532nm激光的光纖固定在支架上,并將其對準待測PMT進行測試,測試系統主要參量如表1所示。

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

Table 1 Main parameters of the test system
對PMT的后脈沖分布進行實驗室標定,并通過手持衰減器對光強進行不同程度的衰減,讀取高速采集卡采集到的數據,然后對這些數據進行處理,具體步驟見下。
(1)PMT本底噪聲修正。由于探測器的緩慢放電效應,探測器的基線存在隨采樣時間緩慢降低的情況。為了得到基線值的擬合函數,在探測器離散回波信號出現之后,讀取不同時刻下探測系統的基線值并作出散點圖,然后利用Curve Fitting工具箱對散點圖進行曲線擬合,如圖4所示。
根據曲線擬合結果,532nm和486nm基線幅值與采樣時間t的關系如下:
B532(t)=1.506e-0.00233(t-800)+
203.2e-3.869×10-6(t-800)
(12)
B486(t)=1.589e-0.03235(t-800)+
204.2e-2.651×10-6(t-800)
(13)
在進行數據處理之前,需要將實驗采集的數據減去基線幅值的變化,以便得到更為準確的接收光子數。
(2)背景光噪聲修正。在使用光電倍增管接收光子時,除了目標反射回來的激光回波信號,還摻雜著一定程度的背景光噪聲。為獲得背景光噪聲的產生概率,需要對所采集到的光子數據進行累加觀察其拖尾現象,累加結果如圖5所示。其中圖5a的結果表明,在光強衰減值為0dB時,光電倍增管的拖尾效應持續的時長約為1000ns,并且通過縱坐標可以大致看出經過20000次累加之后,系統接收到的背景光的數目為6個,因此相應采集到的背景光噪聲概率pn為:

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

Fig.5 Result of accumulating the photons collected by the photomultiplier tube was obtained under the noise
pn=6/20000=3/10000
(14)
為進一步確認背景光噪聲的概率值,在光強衰減值分別為40dB,50dB,60dB時對收集到的數據又進行了相同的處理,其中60dB累加結果如圖5b所示。對比兩個小圖可知,背景光噪聲概率為3/10000。
利用MATLAB編寫背景光噪聲過濾程序,具體程序流程如圖6所示。首先提取需處理的實驗室數據,并將數據的行數和列數分別命名為a,b;然后引入隨機變量i,j以及背景噪聲光子數Np,當滿足條件i≤a且j≤b時,系統借助隨機函數rand產生[0,1]區間內均勻分布的隨機數,并將隨機數大小與背景光噪聲概率3/10000相比較,統計所有小于3/10000的值的個數并賦值給Np;最后將原始數據在經過PMT本底噪聲修正的基礎上減掉背景光噪聲的光子數。

Fig.6 Program flow chart for the background noise of PMT
(3)累加降噪后的結果。對不同光衰減強度下的光電倍增管接收到的光子數進行20000次累加,得到結果如圖7所示(以532nm激光為例)。
由統計結果可知,隨著入射光強的不斷衰減,PMT接收到的光子數逐漸減少,同時后脈沖效應也逐漸減弱。在入射光強衰減值大于30dB后,后脈沖效應將不再明顯,并且入射光下的后脈沖分布已經足夠低。因此,對于入射光強低于30dB的情況,將按照30dB入射光下的后脈沖分布函數進行插值。

Fig.7 The number of photons received by PMT after noise reduction is accumulated
(4)后脈沖概率分布擬合函數。對光電倍增管的后脈沖部分進行擬合,得到探測器后脈沖分布概率隨采樣時間之間的變化關系式。單光子探測器的后脈沖概率分布函數可以用統一的雙指數函數進行描述[15]。
本文中借助的是MATLAB擬合工具箱Curve Fitting,擬合過程分為以下幾個步驟:首先導入不同光衰減下需要進行擬合的數據;然后選擇擬合方程為雙指數函數型;最后得到擬合方程曲線以及擬合函數,結果如圖8所示。

Fig.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]如下:
[Fc,a(i-j)-Fc,b(i-j)]+Fc,b(i-j)
(15)
式中,Fap(j;i)表示第ins所接收到的入射光強所產生的后脈沖在第jns的分布,Nap(i) 表示在第ins所接收到的后脈沖光子數,且Nap(i)∈[Nc,a,Nc,b],Nc,a,Nc,b為實驗室標定的在特定衰減值a,b下系統所采集到的入射光脈沖光子數,Fc,a,Fc,b為入射光子數是Nc,a,Nc,b時的后脈沖概率分布函數。
這一節中通過對實驗系統所得數據的處理,得到了光電倍增管后脈沖的分布情況隨采樣時間變化的函數關系式,并且利用插值公式求解了相鄰兩個入射光強度下的后脈沖分布概率,最后可通過線性插值法來分段統計其后脈沖光子數。
本次實驗中采用中國科學院上海精密機械研究所空間激光信息技術研究中心的研究團隊在我國南海水域測得的數據,測試中所用的激光雷達為自主研制的486nm和532nm藍綠波長多通道海洋激光雷達,飛行所得原始數據如圖9所示。

Fig.9 Flight data of airborne marine lidar
對機載海洋激光雷達飛行數據的校正流程可以概括為以下3個步驟:(1)去除背景光噪聲;(2)結合第3.1節中所得PMT后脈沖概率分布函數與 (15) 式求出位于不同衰減值區間的后脈沖光子數;(3)將原始飛行數據減去背景光噪聲和后脈沖光子數,即可得到探測器接收到的真實光子數。
通過實驗室標定系統得到了PMT的后脈沖概率分布函數,利用后脈沖的概率分布函數對機載海洋激光雷達的機載數據進行后脈沖校正。圖10所示為機載海洋激光雷達測得的數據經過后脈沖校正后的結果。
由圖10可知,對比后脈沖校正前后的數據,校正前的數據顯示后脈沖拖尾效應較為明顯,532nm,486nm激光通道的最大測深深度分別約為110m,200m;經過后脈沖校正后的波形無明顯后脈沖拖尾,且532nm,486nm激光通道所探測到的最大水深為77m,102m左右。由此可得,機載雷達數據經過后脈沖校正后去除了假信號,使得測深距離更為準確,并且由圖可以看出,由于散射層的存在,兩個激光通道的衰減率在40m水深處均有所增加,且隨著海水深度的增加,渾濁程度也隨之增加,486nm激光的優勢減弱,因此水深大于40m后,486nm通道與532nm通道的衰減速率更加接近。

Fig.10 Marine flight data processing results of airborne marine lidar
為驗證后脈沖分布函數法校正機載數據的正確性,利用蒙特卡洛仿真模擬激光水下傳輸過程,并將其結果與實驗結果進行對比。蒙特卡洛方法是在對大量光子隨機運動軌跡計算的基礎上,模擬光在渾濁介質傳播的整個物理過程[18-20]。仿真參量與實驗情況保持一致,具體參量設置如下:模擬的光子包數為10000個,飛機高度H=300m,接收光學天線面積為0.2m2,海水衰減系數C=0.05m-1,最大散射次數設為50,接收機視場角為20mrad,光電倍增管相關系數不變,如表1所示。圖11中顯示了實驗結果與仿真結果的對比。

Fig.11 Comparison of experimental and simulated waveforms
為了更加直觀地反映實驗波形和模擬波形之間的相關關系,本文中將通過計算兩條曲線的相關系數來定量表示其相關程度。具體做法為:分別提取同一波長下的實驗波形和模擬波形數據,將同一橫坐標下對應的兩個點劃分為一組,共提取了663組這樣的點,將實驗波形上提取到的663個點的縱坐標值構成向量y1,同樣的,從模擬波形上提取到的663個點的縱坐標值構成了向量y2,計算y1和y2的相關系數R,最后通過計算得到激光波長為486nm,532nm的實驗波形曲線和模擬波形曲線的相關系數分別為0.9689,0.8648。根據參考文獻[21]可知,相關系數的絕對值大于0.8表示兩個變量高度線性相關,說明實驗方法正確地獲得了激光雷達的真實信號,即對其可行性進行了驗證。
由圖11可知,相比實驗波形,模擬波形看起來更加平滑,這是因為模擬狀態下噪聲信號幾乎為零;并且對比同一波長下的實驗波形與模擬波形可知,隨著探測水域深度的增加,實驗波形的衰減速度明顯要比模擬波形快一些,且實驗測得水深比模擬結果小10m左右,這是由于真實狀態下,海水中的懸浮物和不溶物對光信號產生了漫反射,使得激光雷達未能在有效接收面積內接收更多的光子,而模擬中未考慮這些情況;另外,由于蒙特卡洛仿真在模擬激光水下傳輸時能夠到達接收機的光子十分有限,因此需要耗費大量時間來得到一個可靠的統計結果,而后脈沖概率函數法是在原有機載信號的基礎上進行校正,故可以很好地避免這個缺點。
針對機載海洋激光雷達在使用中普遍存在的光電倍增管后脈沖的問題,本文中通過對光電倍增管后脈沖的標定,獲得了光電倍增管后脈沖概率分布函數,在不搭建仿真模型的前提下正確地校正了機載海洋激光雷達實測數據中的后脈沖區域,并通過與蒙特卡洛仿真結果的對照分析,驗證了光電倍增管標定方法的準確性和可靠性。因此,利用后脈沖概率分布函數能夠更加快捷、簡便地校正光電倍增管的后脈沖效應,節約了大量時間,提高了以光電倍增管為核心器件的機載海洋激光雷達測量數據的精確度,為進一步研究海洋特性奠定了良好的基礎。
感謝中國科學院上海光學精密機械研究所賀巖、羅遠等人提供的實驗設備及機載數據。