Gm-APD 激光雷達(dá)的雙參量估計(jì)透煙霧成像算法

張銀波，李昊陽(yáng)，孫劍峰*，李思寧，姜 鵬，侯 躍，張海龍

（1. 哈爾濱工業(yè)大學(xué) 光電子技術(shù)研究所 可調(diào)諧（氣體）激光技術(shù)重點(diǎn)實(shí)驗(yàn)室，黑龍江 哈爾濱 150001；2. 復(fù)雜系統(tǒng)控制與智能協(xié)同技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100074）

1 引 言

Gm-APD 激光雷達(dá)在室外成像的過(guò)程中，天氣的改變嚴(yán)重影響其成像結(jié)果，特別是云霧、煙霧等散射介質(zhì)。顆粒對(duì)激光有較強(qiáng)的散射和吸收能力，導(dǎo)致回波信號(hào)信噪比下降，激光雷達(dá)探測(cè)距離縮短，成像能力嚴(yán)重受限［1］。因此，研究激光與煙霧相互作用的散射機(jī)理、分析煙霧后向散射回波特點(diǎn)以及建立煙霧散射回波模型，對(duì)于煙霧后向散射抑制、精確提取目標(biāo)回波信號(hào)、提高激光雷達(dá)的探測(cè)能力和環(huán)境感知能力尤其重要。

光子在煙霧中傳輸時(shí)，與粒子之間發(fā)生多次散射，接收到的回波光子在時(shí)域上滿足一定的分布規(guī)律。目前，對(duì)于單光子激光雷達(dá)煙霧回波信號(hào)的研究主要以分布模型擬合去霧法為主。該技術(shù)通過(guò)擬合接收光子在時(shí)域上的分布，利用目標(biāo)信號(hào)與煙霧信號(hào)分布的不同來(lái)實(shí)現(xiàn)去霧。其中，美國(guó)麻省理工學(xué)院（MIT）在分布模型擬合去霧中率先提出了Gamma 分布模型。2018 年，Satat 等［2］通過(guò)觀察發(fā)現(xiàn)，煙霧回波光子與目標(biāo)回波光子分別滿足Gamma 分布和高斯分布，通過(guò)極大似然估計(jì)得到Gamma 分布的兩個(gè)參數(shù)，隨后將目標(biāo)光子從煙霧散射光子中分離，最終實(shí)現(xiàn)煙霧的去除。在不同的煙霧光學(xué)厚度下，達(dá)到了較好的重建效果。不同于Satat 提出的Gamma 分布理論，2020 年Sang［3］等提出了基于激光雷達(dá)物理模型估計(jì)霧輪廓線的去霧算法，同時(shí)采用期望最大化（EM）算法進(jìn)行參數(shù)最優(yōu)化搜尋。同年，Mau［4］等 利 用EM 算 法 將SPAD 陣 列 中 每 個(gè) 像 素點(diǎn)的探測(cè)概率曲線擬合為對(duì)數(shù)正態(tài)分布與高斯分布的混合分布，并依據(jù)高斯分布參數(shù)得到目標(biāo)距離信息。哈爾濱工業(yè)大學(xué)［5-10］開展了基于Gm-APD 激光雷達(dá)的煙霧后向散射機(jī)理研究，建立了適用于Gm-APD 激光雷達(dá)的后向散射光子分布模型。根據(jù)具備遠(yuǎn)距離探測(cè)能力的Gm-APD 激光雷達(dá)，提出多種基于Gamma 模型參數(shù)估計(jì)的霧天激光雷達(dá)重構(gòu)算法，并在室內(nèi)、室外環(huán)境中驗(yàn)證了模型和算法的能力。

以上研究表明，目前的單光子激光雷達(dá)透煙霧成像主要以基于Gamma 模型參數(shù)估計(jì)的煙霧去除為主導(dǎo)。估計(jì)算法主要通過(guò)極大似然估計(jì)直接獲取Gamma 函數(shù)的兩個(gè)參數(shù)，并未對(duì)兩個(gè)參數(shù)的關(guān)系進(jìn)行討論。此外，目標(biāo)峰的存在影響煙霧回波信號(hào)的估計(jì)精度，降低圖像的重構(gòu)質(zhì)量。

針對(duì)上述問(wèn)題，本文從Gamma 模型參數(shù)關(guān)系推導(dǎo)和重構(gòu)算法兩個(gè)方面展開研究。首先，介紹了基于Gm-APD 激光雷達(dá)的觸發(fā)模型以及根據(jù)探測(cè)概率求解實(shí)際接收的回波信號(hào)原理，基于光子與煙霧粒子的碰撞理論，以Mie 散射理論為基礎(chǔ)詳細(xì)推導(dǎo)Gamma 模型的兩個(gè)參數(shù)，并且建立衰減系數(shù)μ與散射次數(shù)k之間的物理關(guān)系。根據(jù)推導(dǎo)的μ，k關(guān)系式提出一種雙參量估計(jì)算法，考慮了如何精確估計(jì)μ，k兩個(gè)參數(shù)。最后，利用仿真實(shí)驗(yàn)檢驗(yàn)μ，k關(guān)系式的正確性，利用室內(nèi)實(shí)驗(yàn)驗(yàn)證了提出算法的透煙霧成像能力。

2 透霧理論模型

2.1 Gm-APD 激光雷達(dá)觸發(fā)模型

Gm-APD 激光雷達(dá)在探測(cè)的過(guò)程中，目標(biāo)光子和噪聲光子組成的回波信號(hào)入射到探測(cè)器光敏面上實(shí)現(xiàn)光電子轉(zhuǎn)換。目標(biāo)反射光子和噪聲光子均滿足Poisson 分布，因此，在t1和t2探測(cè)時(shí)間間隔內(nèi)產(chǎn)生k個(gè)光電事件的概率［9，11］為：

式中：M（t1，t2）為在時(shí)間區(qū)間（t1，t2）內(nèi)探測(cè)器探測(cè)到的回波光子數(shù)。

Gm-APD 在時(shí)間t1和t2之間不發(fā)生觸發(fā)，即產(chǎn)生0 個(gè)光電子事件的概率為exp［-M（t1，t2）］。產(chǎn)生光電事件的概率為1-exp［-M（t1，t2）］。然而，由于探測(cè)器死時(shí)間的存在［16］，第i個(gè)時(shí)間bin的探測(cè)概率與之前時(shí)間bin 的探測(cè)概率有關(guān)。若門控時(shí)間小于探測(cè)器死時(shí)間，在距離門內(nèi)Gm-APD 探測(cè)器只能響應(yīng)一次光子計(jì)數(shù)信號(hào)，即單次觸發(fā)模式。因此，第i個(gè)時(shí)間bin 中的探測(cè)概率為［17］：

根據(jù)探測(cè)得到的探測(cè)概率曲線，結(jié)合式（2）可以得到第i個(gè)bin 內(nèi)的回波光子數(shù)為：

通過(guò)求解距離門內(nèi)每個(gè)時(shí)間bin 內(nèi)的光子數(shù)，最終得到實(shí)際接收的回波信號(hào)。本文的主要工作就是基于實(shí)際接收的回波信號(hào)進(jìn)行目標(biāo)信號(hào)的提取與圖像重建。

2.2 基于Gamma 模型的參數(shù)推導(dǎo)

由于煙霧顆粒對(duì)激光具有強(qiáng)散射和吸收作用，激光在煙霧中傳播時(shí)的回波能量隨距離的增長(zhǎng)呈e 指數(shù)衰減，并且激光傳輸時(shí)光子存在連續(xù)散射事件，每個(gè)散射事件相互獨(dú)立，多個(gè)獨(dú)立指數(shù)隨機(jī)變量的總和服從Gamma 分布。其表達(dá)式如下［2］：

式中：μ代表衰減系數(shù)，k代表碰撞次數(shù)。

Gamma 模型涉及的兩個(gè)參數(shù)，衰減系數(shù)μ和碰撞次數(shù)k，均與激光在煙霧中的傳輸過(guò)程密切相關(guān)。本文從光子與煙霧顆粒的碰撞過(guò)程這個(gè)微觀角度出發(fā)，推導(dǎo)這兩個(gè)參數(shù)的關(guān)系。

激光在煙霧中傳播時(shí)其衰減滿足朗伯比爾定律：

式中：Iin，Iout分別為激光透過(guò)煙霧前后的光強(qiáng)，μ為介質(zhì)的衰減系數(shù)，l為激光在介質(zhì)中傳播的距離。

將式（5）變形得到衰減系數(shù)μ的表達(dá)式為：

假設(shè)一次成像過(guò)程結(jié)束后，全部光子的碰撞次數(shù)的均值為k，由米氏理論可知光子與散射介質(zhì)顆粒發(fā)生單次碰撞會(huì)損耗能量，損耗的比例為散射效率因子與消光效率因子的比值，即Qsca/Qext。

散射效率因子與消光效率因子的求解公式如下［12］：

式中：α=πd/λ，為煙霧顆粒的尺寸因子；d為煙霧平均粒徑；λ為激光波長(zhǎng)；an，bn為Mie 散射系數(shù)，他們的值與粒子尺寸因子α和介質(zhì)顆粒復(fù)折射率m密切相關(guān)［18］。

當(dāng)激光在衰減系數(shù)為μ的散射介質(zhì)中傳輸l路徑時(shí)，與煙霧顆粒發(fā)生k次碰撞后的光強(qiáng)近似滿足以下關(guān)系：

將式（8）帶入式（6）中得到Gamma 函數(shù)兩個(gè)參量μ和k之間的關(guān)系，即：

可以看出，衰減系數(shù)μ與散射次數(shù)k滿足線性關(guān)系，即隨著散射次數(shù)k的增加，衰減系數(shù)μ呈線性增大的趨勢(shì)。然而，增大的斜率與煙霧自身的消光系數(shù)、散射系數(shù)和激光在煙霧中傳輸?shù)拈L(zhǎng)度等參量有關(guān)。相較于文獻(xiàn)［2］中模型參量的定義，本文從Mie 散射理論出發(fā)，建立了衰減系數(shù)μ與散射次數(shù)k之間的物理關(guān)系。

2.3 透霧重建算法

Gm-APD 激光雷達(dá)透霧成像的關(guān)鍵是能夠精確地估計(jì)煙霧的散射信號(hào)，進(jìn)而實(shí)現(xiàn)煙霧散射光子的抑制?；贕amma 分布模型的煙霧回波信號(hào)估計(jì)，其研究核心是實(shí)現(xiàn)兩個(gè)參量μ與k的精確估計(jì)。不同于采用極大似然估計(jì)算法直接求解μ與k［2］，本文的參 量估計(jì)算 法 是基于2.2 節(jié)中μ與k的線性關(guān)系估計(jì)這兩個(gè)參數(shù)的。

參數(shù)估計(jì)主要采用兩步策略。第一步是利用Gamma 模型的數(shù)學(xué)性質(zhì)求解衰減系數(shù)μ。由式（4）的Gamma（μ，k）可知，其均值E和方差D分別 滿 足E=k/μ，D=k/μ2。衰 減 系 數(shù)μ為 均 值E和方差D的比值，即μ=E/D。根據(jù)式（3）計(jì)算得到的回波信號(hào)，可以得到每個(gè)像素接收回波信號(hào)的均值E和方差D。第二步是利用本文推導(dǎo)得到的μ，k關(guān)系求解k。其中，傳輸距離l為實(shí)際目標(biāo)距離，散射效率因子與消光效率因子可根據(jù)實(shí)際的煙霧類型確定。

由于回波信號(hào)中包含目標(biāo)反射光子與煙霧散射光子，這嚴(yán)重限制了均值E和方差D的求解。因此，這里采用DOG 小波函數(shù)對(duì)回波信號(hào)進(jìn)行連續(xù)小波變換［19］，通過(guò)峰值選取確定目標(biāo)回波峰并進(jìn)行去除。

本文提出的基于μ，k線性關(guān)系的雙參量估計(jì)算法，通過(guò)估計(jì)面陣Gm-APD 激光雷達(dá)中每個(gè)像素的回波信號(hào)來(lái)實(shí)現(xiàn)煙霧干擾抑制。算法的具體流程如圖1 所示。

圖1 雙參量估計(jì)算法處理流程Fig.1 Processing flow chart of dual-parameter estimation algorithm

3 實(shí) 驗(yàn)

為了檢驗(yàn)本文提出算法對(duì)煙霧干擾的抑制能力，開展了室內(nèi)高濃度煙霧條件下的Gm-APD激光雷達(dá)透霧成像實(shí)驗(yàn)，如圖2 所示。Gm-APD激光雷達(dá)采用激光脈沖飛行時(shí)間測(cè)量（TOF）方式獲得目標(biāo)的三維距離信息。探測(cè)器的觸發(fā)信號(hào)為激光脈沖同步信號(hào)。采用FPGA 實(shí)現(xiàn)控制命令的發(fā)送以及圖像數(shù)據(jù)的傳輸。目標(biāo)的三維圖像在計(jì)算機(jī)上顯示。雷達(dá)的具體參數(shù)見文獻(xiàn)［9］。

圖2 Gm-APD激光雷達(dá)透霧成像實(shí)驗(yàn)示意圖及場(chǎng)景Fig.2 Schematic diagram and experimental scene for Gm-APD lidar imaging through fog

為了減小煙霧的擴(kuò)散范圍、提高其分布的均勻性，本實(shí)驗(yàn)采用長(zhǎng)度為1 m，直徑為450 mm 的塑料箱作為煙霧箱，其兩端安裝具有90%透過(guò)率的光學(xué)透鏡。煙霧箱前端與雷達(dá)的距離為5 m。目標(biāo)在煙霧箱后端。本文的煙霧發(fā)生器可產(chǎn)生粒徑分布為1～2 μm 的顆粒，1 064 nm 激光對(duì)應(yīng)的復(fù)折射率與煤煙型顆粒折射率接近。

為避免背景光子對(duì)實(shí)驗(yàn)結(jié)果的影響，本文的實(shí)驗(yàn)在暗室中進(jìn)行；當(dāng)煙霧箱中不充入煙霧時(shí)，對(duì)目標(biāo)進(jìn)行探測(cè)，利用傳統(tǒng)的峰值算法進(jìn)行圖像重建，并將這些成像結(jié)果作為無(wú)煙霧時(shí)的標(biāo)準(zhǔn)圖像；利用煙霧發(fā)生器向霧箱內(nèi)充入煙霧，待煙霧在霧箱中擴(kuò)散均勻后，利用激光雷達(dá)透過(guò)煙霧對(duì)目標(biāo)進(jìn)行數(shù)據(jù)采集；利用不同的算法對(duì)有煙霧的圖像進(jìn)行去霧成像，并與標(biāo)準(zhǔn)圖像進(jìn)行對(duì)比，檢驗(yàn)算法的去霧成像能力。

4 結(jié)果與分析

4.1 蒙特卡洛仿真驗(yàn)證

為了驗(yàn)證式（9）的正確性，本文首先利用蒙特卡洛算法實(shí)現(xiàn)煙霧環(huán)境下的回波信號(hào)仿真。在仿真中，選用煤煙型顆粒的復(fù)折射率，m=1.75i～0.44i［13］，煙霧顆粒平均粒徑d為1 μm［14］，激光波長(zhǎng)λ為1.064 μm。通過(guò)調(diào)節(jié)能見度參數(shù)，實(shí)現(xiàn)不同煙霧濃度條件下的激光后向散射回波光子分布的仿真。 蒙特卡洛仿真算法見文獻(xiàn)［7］。

當(dāng)能見度V為0.12 km 時(shí)，回波光子的時(shí)域分布結(jié)果如圖3（a）所示。可以看出，回波信號(hào)主要為雙峰分布，第一個(gè)峰為煙霧散射光子，第二個(gè)峰為目標(biāo)峰。為了提高Gamma 擬合參數(shù)估計(jì)結(jié)果的準(zhǔn)確性，估計(jì)算法只需要對(duì)煙霧散射光子進(jìn)行討論。因此，需要對(duì)回波信號(hào)中的目標(biāo)回波信號(hào)峰進(jìn)行剔除。此處選用一次函數(shù)代替目標(biāo)回波峰，對(duì)處理后的數(shù)據(jù)進(jìn)行Gamma 擬合得到標(biāo)準(zhǔn)的模型參數(shù)。去除目標(biāo)峰后的擬合結(jié)果如圖3（b）所示，Gamma 模型擬合結(jié)果與原始數(shù)據(jù)更貼切。

圖3 仿真數(shù)據(jù)及擬合結(jié)果Fig.3 Simulation data and fitting results

為了提升數(shù)據(jù)分析的合理性與精度，在每個(gè)能見度下進(jìn)行多次仿真，并且對(duì)去除信號(hào)峰后的數(shù)據(jù)進(jìn)行擬合。表1 所示為能見度V從0.12 km增加至0.16 km 時(shí)的擬合結(jié)果及誤差分布。標(biāo)準(zhǔn)參數(shù)μ，k是由十組數(shù)據(jù)擬合的參數(shù)求均值得到。 當(dāng)V=0.12 km 時(shí)，k的 標(biāo) 準(zhǔn) 值 均 值 為1.944，方差為0.013。由表1 可知，利用μ，k關(guān)系公式計(jì)算的參數(shù)與極大似然估計(jì)的參數(shù)誤差均在7%以下，證明了式（9）的有效性。

表1 不同能見度條件下公式求解參數(shù)的誤差分布Tab.1 Error distribution of formula solution parameters under different visibility conditions

4.2 室內(nèi)實(shí)驗(yàn)數(shù)據(jù)處理與分析

4.2.1 單個(gè)像素參數(shù)估計(jì)

為了檢驗(yàn)本文提出算法對(duì)煙霧回波信號(hào)的估計(jì)能力，對(duì)室內(nèi)透霧成像數(shù)據(jù)的單個(gè)像素回波進(jìn)行擬合。選擇像素點(diǎn)（20，45）的信號(hào)進(jìn)行分析，結(jié)果如圖4（a）所示?；夭ㄖ邪瑑蓚€(gè)信號(hào)峰（煙霧散射信號(hào)和目標(biāo)信號(hào)），煙霧散射信號(hào)的強(qiáng)度大于目標(biāo)的回波信號(hào)，并且目標(biāo)峰在煙霧回波信號(hào)的下降沿內(nèi)。

圖4 單個(gè)像素回波信號(hào)的擬合結(jié)果Fig.4 Fitting results of single pixel echo signal

分別利用文獻(xiàn)［2］中MIT 的擬合算法和本文提出的算法對(duì)這個(gè)回波信號(hào)進(jìn)行估計(jì)。值得注意的是，這兩種算法都是直接估計(jì)計(jì)算得到回波光子信號(hào)。我們研究的重點(diǎn)是對(duì)包含煙霧散射光子和目標(biāo)反射光子的回波信號(hào)進(jìn)行估計(jì)。然而，兩種算法的擬合結(jié)果差別較大。從圖4（a）中可見，MIT 算法的擬合效果較差。這是由于目標(biāo)回波峰的存在影響了Gamma 擬合的精度。Gamma 擬合得到的曲線將目標(biāo)回波覆蓋，導(dǎo)致在對(duì)消過(guò)程中目標(biāo)信息的丟失。從圖4（c）中可見，雙參量估計(jì)法的擬合效果更優(yōu)。算法通過(guò)兩步策略的估計(jì)方式，實(shí)現(xiàn)了只對(duì)煙霧回波峰數(shù)據(jù)的擬合，抑制了目標(biāo)信號(hào)對(duì)估計(jì)結(jié)果的影響。

分析圖4（b）和4（d）兩算法的殘差分布可知，MIT 算法殘差分布中的目標(biāo)信息完全丟失，錯(cuò)誤地將煙霧的回波信息當(dāng)作目標(biāo)回波信息。雙參量估計(jì)算法殘差分布圖中目標(biāo)信息保留完整，克服了煙霧回波殘差對(duì)目標(biāo)信息提取精度的影響，提高了對(duì)煙霧背后目標(biāo)信息的感知能力。

4.2.2 室內(nèi)透霧圖像重構(gòu)

為了檢驗(yàn)提出算法對(duì)煙霧背后目標(biāo)的成像能力，本文對(duì)面陣Gm-APD 激光雷達(dá)64×64 個(gè)像素的回波信號(hào)進(jìn)行目標(biāo)光子提取，并且與傳統(tǒng)算法進(jìn)行對(duì)比。重建圖像如圖5 所示，可以看出，本文提出算法重構(gòu)的距離像與強(qiáng)度像均優(yōu)于峰值法和MIT 算法。

圖5 不同算法的圖像重建結(jié)果Fig.5 Image reconstruction results of different algorithms

由距離像可知，峰值法與MIT 算法重建的圖像不能識(shí)別和保留目標(biāo)像素點(diǎn)的信息，無(wú)法重構(gòu)出濃霧背后的目標(biāo)輪廓信息和深度信息，這兩種算法的抗煙霧干擾能力較弱。而雙參量估計(jì)算法很好地恢復(fù)了被煙霧遮擋目標(biāo)的輪廓信息，較好地重構(gòu)出目標(biāo)的深度圖像。

由強(qiáng)度像可知，峰值法和MIT 算法對(duì)強(qiáng)度像的恢復(fù)程度較低，兩個(gè)目標(biāo)反射光的強(qiáng)度與背景融為一體，很難區(qū)分目標(biāo)的輪廓。雙參量估計(jì)算法重構(gòu)的強(qiáng)度像清晰反映了目標(biāo)表面的強(qiáng)度分布情況，與背景明顯地分離開。

為了對(duì)比3 種算法的重建能力，本文選用結(jié)構(gòu)相似性（Structural Similarity，SSIM）［20］和目標(biāo)復(fù)原度（Target Recovery，TR）［9］兩個(gè)指標(biāo)來(lái)評(píng)價(jià)算法的圖像重建效果。SSIM 取值為［0，1］，SSIM 越大，圖像失真越小、圖像效果越好；TR 值越大，算法對(duì)場(chǎng)景的復(fù)原能力越強(qiáng)。

3 種算法重建圖像計(jì)算的評(píng)價(jià)指標(biāo)都是基于標(biāo)準(zhǔn)圖像得到的，結(jié)果如表2 所示。相對(duì)于峰值法和MIT 算法，本文算法具有最大的SSIM 和TR。與MIT 算法相比，本文算法重建強(qiáng)度像的SSIM 提高了0.228 9，距離像的TR 提高了73%。這表明本文提出的雙參量估計(jì)算法在面陣Gm-APD 激光雷達(dá)透煙霧成像中具有更強(qiáng)的煙霧抑制能力和目標(biāo)感知能力。

表2 不同算法處理圖像的評(píng)價(jià)指標(biāo)Tab.2 Evaluation index of image processed by different algorithms

5 結(jié) 論

為了提高Gm-APD 激光雷達(dá)對(duì)濃密煙霧背后的目標(biāo)的成像質(zhì)量，本文對(duì)MIT 提出的Gamma 分布模型進(jìn)行了完善，提出了一種基于雙參量估計(jì)的Gm-APD 激光雷達(dá)透煙霧成像算法，主要包括Gamma 模型參數(shù)關(guān)系推導(dǎo)和重構(gòu)算法兩個(gè)方面。首先，介紹了基于Gm-APD 激光雷達(dá)的觸發(fā)模型以及根據(jù)探測(cè)概率求解實(shí)際接收的回波信號(hào)原理，并且基于光子與煙霧粒子的碰撞理論和Mie 散射理論詳細(xì)推導(dǎo)了Gamma 模型兩個(gè)參數(shù)的物理關(guān)系。根據(jù)推導(dǎo)的關(guān)系式提出一種雙參量估計(jì)算法，并開展了仿真和室內(nèi)透霧的實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明，雙參量估計(jì)算法的重構(gòu)效果優(yōu)于MIT 算法，重構(gòu)的強(qiáng)度像的結(jié)構(gòu)相似性提高了0.228 9，重構(gòu)的距離像的目標(biāo)復(fù)原度提高了73%。該研究有效提升了Gm-APD 激光雷達(dá)在煙霧環(huán)境中的目標(biāo)感知能力，拓展了復(fù)雜環(huán)境下Gm-APD 激光雷達(dá)的應(yīng)用范圍。