時空域聯(lián)合編碼擴頻單光子計數(shù)成像方法*

沈姍姍 顧國華 陳錢 何睿清 曹青青

1) (南京工業(yè)職業(yè)技術大學航空工程學院,南京 210023)

2) (南京理工大學電子工程與光電技術學院,智能感知和微光成像實驗室,南京 210094)

3) (南京工程學院信息與通信工程學院,南京 211167)

本文結合空間編碼單像素成像技術和擴頻時間編碼掃描成像技術,提出一種時空域聯(lián)合編碼擴頻單光子計數(shù)成像方法.該方法具備可避免距離模糊、大時寬帶寬積的優(yōu)勢,并且在噪聲干擾下,能夠準確恢復距離像.本文推導了基于單光子探測的時空域聯(lián)合相關非線性探測模型、成像正向模型和信噪比模型,并通過凸優(yōu)化反演算法恢復深度圖像,理論模型和仿真實驗均證明,與傳統(tǒng)的基于空間編碼的單像素成像方法相比,本方法提高了重建的質量.其中,采用m 序列作為時間編碼,成像的噪聲魯棒性更高.和傳統(tǒng)的空間編碼單像素成像技術相比,本文提出方法的成像均方誤差降低了4/5,引入二次相關方法后,成像均方誤差降低了9/10.本文所提出的成像架構為非掃描激光雷達成像方法提供了新思路.

1 引言

單光子計數(shù)三維成像技術[1?3]廣泛應用于機器視覺、遙感和無人駕駛等領域.其中,利用目標場景的可壓縮性,從欠采樣線性投影序列中重建三維圖像是近年來廣泛研究的單光子計數(shù)三維成像方法.該方法在空間光調制器(digital micro mirror device,DMD)上加載一系列的隨機編碼.進而在空間域調制光信號,單光子探測器(single photon avalanche detector,SPAD)用于測量場景和隨機編碼的內積,從遠小于像素大小的投影中獲取目標信息,并在稀疏約束下根據(jù)欠采樣測量值重建圖像.Howland 等[4]從當前測量的信號中減去前一次測量的參考信號以削弱噪聲能量,實現(xiàn)32×32 分辨率高幀速動態(tài)目標跟蹤.Zhao 等[5]采用Hadamard基提出一種高效的單像素全彩成像方法.在微光環(huán)境下,Radwell 等[6]提出了一種基于深度學習的單像素成像優(yōu)化算法.李鳳強等[7]利用L2范數(shù)進行約束,并采用TwIST(two-step iterative shrinkage/thresholding)反演方法重建圖像,所提出的100 萬像素的壓縮感知(compressed sensing,CS)激光雷達的成像空間分辨率提高了4 倍.Liu 等[8]提出一種SP3I 技術以提高系統(tǒng)的接收效率和輸出信噪比(signal-to-noise ratio,SNR).Asmann 等[9]提 出一種陣列的壓縮感知激光雷達系統(tǒng),成像信噪比逼近全幀互相關方法的成像信噪比,并分析了全變分增強拉格朗日交替方向法(total variation augmented Lagrangian alternation direction algorithm,TVAL3)和梯度投影稀疏重建法(gradient projection for sparse reconstruction,GPSR)的重建結果.傳統(tǒng)的基于空間編碼的單像素激光雷達通常以固定頻率(10 kHz 至10 MHz)發(fā)射激光脈沖,限制了最大不模糊距離和信號探測速率.此外,分辨率和帶寬之間的矛盾難以調和,成像結果易受噪聲影響,使得此類方法成像信噪比低,僅限于室內和短距離成像的應用場合.

擴頻通信技術廣泛應用于無線傳感器網(wǎng)絡[10]、全球定位系統(tǒng)和掃描式激光雷達[11?18]等領域,其中,基于擴頻通信時間編碼的單光子計數(shù)掃描成像方法[11?18],通過接收端的匹配濾波器實現(xiàn)脈沖壓縮,高效提高輸出信噪比并減少信號衰減的概率,具有可避免距離模糊、大時寬帶寬積、低功耗和易于單片集成的優(yōu)勢.該方法通常采用偽隨機序列進行時間編碼,Ding 等[15]采用加權的多距離時間相關函數(shù)作為目標函數(shù),通過模擬退火方法優(yōu)化偽隨機序列結構,提高時間相關函數(shù)的距離分辨能力.Hiskett 等[16]設計“1”碼概率為0.1 的序列,并指出若“1”碼概率增加,噪聲相關水平也會增加.

本文結合空間編碼單像素成像技術和擴頻時間編碼掃描成像技術,提出了一種基于時空域聯(lián)合編碼的擴頻單光子計數(shù)三維成像架構.本方法采用偽隨機序列在時間域調制激光脈沖,替代激光脈沖的周期性調制方法,并將調制后的光脈沖投影到DMD,在空間域調制激光脈沖.SPAD 用于探測返回的光脈沖信號,最后將接收信號和參考信號互相關,從而建立時空域聯(lián)合相關非線性探測的數(shù)學模型,該模型將場景點多路復用到單光子探測器上,進而將壓縮感知凸優(yōu)化反演技術(如TwIST[19])引入到新的成像架構中,由此,建立成像正向模型和輸出信噪比模型,并研究死時間對成像質量的影響規(guī)律.仿真實驗和理論模型共同研究表明,通過優(yōu)化時域編碼的平衡度,可降低環(huán)境噪聲和暗計數(shù)的干擾,采用更長的序列可以獲得更高的信噪比和深度確度.此外,針對一般的時域編碼序列,理論模型和仿真實驗表明,引入的二次相關重建方法[20,21]可降低噪聲影響,提高成像質量.

2 成像理論

2.1 擴頻時間編碼單光子計數(shù)成像方法

擴頻時間編碼單光子計數(shù)成像系統(tǒng)中包括偽隨機序列發(fā)生器,其用于生成0,1 碼流,碼為1 時觸發(fā)激光器發(fā)射光脈沖,激光脈沖經(jīng)目標后漫反射得到回波信號,經(jīng)由光學系統(tǒng)接收,到達SPAD 感光面.SPAD 啟動一次探測并驅動光子到達時間記錄儀(time-to-digital converter,TDC)產(chǎn)生時間值.經(jīng)過上述發(fā)射接收過程,完成單個像素的有效探測.定義成像場景的反射率為r(x,y),激光脈沖由SPAD 和TDC 接收后,通過匹配濾波器進行脈沖壓縮,也即將接收序列s(t) 與發(fā)送副本序列f(t)互相關,則像素 (x,y)th的互相關信號為

其中T表示一個周期的序列長度;PR為信號光功率;?t為碼元寬度.選取適當?shù)木嚯x門信號可以消除后向散射峰值,(1)式中第一個最高的互相關峰值對應的橫坐標τ(x,y)為從場景目標反射的光子到達時間的估計值.因此,可采用最大似然估計法(maximum likelihood estimation,ML)[14]估 計τ(x,y)ML:

目標的距離值為

其中c為光速.

2.2 時空域聯(lián)合編碼擴頻單光子計數(shù)成像方法

本文提出的時空域聯(lián)合編碼擴頻單光子技術成像架構如圖1 所示.采用基于可編程邏輯器件的序列發(fā)生器發(fā)送擴頻偽隨機序列.并驅動激光器發(fā)射調制后的激光脈沖,在時間域完成激光脈沖的編碼;采用空間光調制器調制返回的激光脈沖,在空間域完成激光脈沖的編碼.采用偏振分束鏡將從空間光調制器反射的激光光束與從目標反射的激光光束分開,透鏡收集時空域聯(lián)合編碼的激光光束并導入SPAD,信號處理器記錄光子到達時間點,并且和模板序列做互相關運算,通過后續(xù)信號處理重建目標的三維圖像.設先后被第i次時域擴頻偽隨機序列和空域隨機矩陣調制后的光束,依次通過偏振分束鏡和透鏡,到達單光子探測器的總光子計數(shù)率記為Ri(t) :

圖1 時空域聯(lián)合編碼擴頻單光子計數(shù)成像架構Fig.1.Time-space united coding spread spectrum single photon counting imaging method architecture.

式中N表示x軸或y軸方向的像素個數(shù);Φi(x,y) 表示第i次調制的隨機矩陣;Bi1(x,y,t) 表示空域外部環(huán)境噪聲光子計數(shù)率;Bi2(t) 表示時域外部環(huán)境噪聲光子計數(shù)率.第i次時空域聯(lián)合調制后的光束被SPAD 探測后產(chǎn)生響應,該過程為非齊次泊松過程.假設探測器的量子效率為η,則探測到的光子計數(shù)率為ηRi(t) .探測器暗電流產(chǎn)生響應的過程為齊次泊松過程,d表示暗電流響應速率,可得探測到的光子計數(shù)率為

采用相關接收技術,將(5)式和模板序列f(t) 互相關,時空域聯(lián)合相關的非線性探測模型為

(6)式中,第一項為時空相關信號,外部環(huán)境噪聲Bi1(x,y,t) 和Bi2(t)與模板序列互相關后分別得到噪聲項暗計數(shù)噪聲d和模板序列互相關得到噪聲項Bd(τ) .

測量向量為C ∈RM×1,空間光調制器加載的測量矩陣為Φ ∈RM×Q,投影次數(shù)為M,噪聲隨機向量為b ∈RQ×1,像素總數(shù)為Q=N2.(6)式表明Ci(τ)包括信號光子Cis和噪聲光子Cin,本文通過給出Ci(τ)中變量的統(tǒng)計特征來建立信噪比模型.在微光環(huán)境下,單光子探測服從泊松過程[22],假設噪聲和信號是獨立的,且期望和方差相等.則第i次調制得到的Ci(τ) 中包含的信號光子Cis的數(shù)學期望為

(11)式表明輸出信噪比與序列長度成正比.分母的后三項均包括f(t+τ),如果序列中“1”碼的概率幾乎等于“–1”碼的概率,則噪聲能量可忽略不計.因此,通過合理配置時域編碼序列正負比特的平衡性,可降低噪聲能量.其次,考慮探測器的死時間,引入探測到的信號光子或噪聲光子的總概率,死時間增大,雪崩觸發(fā)概率降低,探測到信號光子的總概率降低,輸出信噪比降低.

2.3 二次相關重建方法

根據(jù)(6)式,互相關函數(shù)Ci(τ) 中包含信號和噪聲項,其中,第一項是信號項.只有當目標平面和調制平面之間在空間上滿足嚴格的一對一相關性,以及接收到的光子到達時間點和發(fā)送的副本序列之間在時間上滿足嚴格的一對一的相關性,才能生成無噪聲干擾的信號項.噪聲主要來自系統(tǒng)和環(huán)境,這兩類噪聲不相關,將兩者相乘不會生成直流分量.因此考慮通過二次相關來降低噪聲.假設矩陣Φ的每一列表示為φz ∈RM×1,M=1,2,???,Q,則加性噪聲b可以等效為隨機矩陣V與信號p的乘積,表達為b=V p,其中V與Φ非相關.對(7)式做二次相關:

其中?∈RQ×1,H=cov(Φ,Φ) .壓縮感知充分利用信號的可壓縮性,通過優(yōu)化的方法恢復稀疏表示的測量值.如果p在稀疏基Ψ的表示下是足夠稀疏的,則有p=ΨZ,其中Z是稀疏系數(shù).(18)式等效為?=HΨZ=AZ,其中A為感知矩陣且A=HΨ,采用以下凸優(yōu)化反演方法重建圖像

其中λ是正則化參數(shù);κ(Z) 是正則項.實驗中,采用TwIST[19]方法求解(19)式.為了平衡計算的魯棒性和復雜性,引入最大似然估計(ML)來優(yōu)化重建.

3 實驗結果與分析

為了驗證本文提出成像方法的有效性,采用MATLAB 進行成像仿真實驗.由MATLAB 生成10 GHz 偽隨機序列驅動激光器發(fā)射激光脈沖,偽隨機序列中1 碼的概率由隨機數(shù)閾值確定,設置為0.1,長度為2048.灰度圖像“飛機”的深度范圍為20—20.8 m,如圖2(a)所示.將生成的M個64×64 大小的高斯矩陣依次加載至DMD 上,接收端數(shù)字化時間編碼和空間編碼的信號后,結合TwIST 和ML 法重建圖像.圖2(b)—(g)給出采用不同壓縮比(5%—85%)進行重建的結果,壓縮比表示壓縮的程度,定義為,壓縮程度越高,測量次數(shù)越少,重建速度越快,但重建質量越差.采用不同壓縮比探測和重建圖像的總時間在8—78 s 之間.圖2(e)—(g)中45%或更大的壓縮比可實現(xiàn)高分辨率和高對比度的重建.由于ML 的穩(wěn)定性,本方法的重建結果與參考文獻[7,24]中強度圖像重建工作所得出的結論一致,即25%的壓縮比也可清晰呈現(xiàn)目標場景,如圖2(d)所示.M序列和Gold 序列是CDMA 系統(tǒng)上行鏈路中最常用的偽隨機序列[25,26].為了比較不同序列的抗噪聲能力,在仿真中引入泊松噪聲.平均噪聲光子計數(shù)之和定義為三種噪聲光子計數(shù)值之和,記為mn.假設平均信號光子計數(shù)ms不變且ms=5.噪聲干擾由泊松隨機數(shù)發(fā)生器生成,將mn=1的噪聲引入成像仿真,由圖3(a)和圖3(b),1 碼概率為0.3 的偽隨機序列成像效果優(yōu)于1 碼概率為0.1 的序列成像效果.由于m序列中1 碼的概率和–1 碼的概率基本相同,圖3(d)中所示的m序列的成像效果優(yōu)于其他序列的成像效果,與(11)式序列平衡度和成像信噪比的關系相互驗證.考慮單光子探測器死時間,根據(jù)(11)式,數(shù)值模擬死時間和成像信噪比的制約關系,如圖3(a)所示.隨著死時間的增大,信號光子探測概率降低,信噪比降低;相同死時間下,碼流速度越大,在死時間內到達探測器的信號光脈沖數(shù)量越多,探測到的信號光子數(shù)目減少,信噪比降低.16384 碼長下1 碼概率為0.3,引入不同死時間后的成像仿真結果如圖3(e)—(h)所示,實驗結果表明死時間減小,信噪比增大,成像質量提高,變化規(guī)律和理論模型數(shù)值模擬結果一致.定義深度圖像重建質量的均方誤差(mean squared error,MSE)為

圖2 不同壓縮比的時空域聯(lián)合編碼擴頻單光子計數(shù)成像 (a) “飛機”深度真值;(b)?(g)不同壓縮比的深度成像 (b) 5%,(c)15%,(d) 25%,(e) 45%,(f) 65%,(g) 85%.顏色條表示深度數(shù)值,單位為mFig.2.Time-space united coding spread spectrum single photon counting image with different compression proportions: (a) Depth ground truth of ‘a(chǎn)irplane’;(b)?(g) depth maps with different compression proportions of (b) 5%,(c)15%,(d) 25%,(e) 45%,(f) 65%and (g) 85%.Colorbars show the depth with unit of meter.

圖3 時空域聯(lián)合編碼擴頻單光子計數(shù)成像性能模擬 (a) 1 碼概率為0.1 的偽隨機序列深度圖像;(b) 1 碼概率為0.3 的偽隨機序列深度圖像;(c) Gold 序列的深度圖像;(d) m 序列的深度圖像;(e) 死時間對成像性能的影響理論模擬;(f) 死時間為200 ns;(g) 死時間為20 ns;(h) 死時間為1ns.顏色條表示深度數(shù)值,單位為mFig.3.Time-space united coding image spread spectrum single photon counting imaging performance simulation: (a) Depth maps by pseudo-random sequences with ‘1’ bit of 0.1;(b) depth maps by pseudo-random sequences with ‘1’ bit of 0.3;(c) depth maps by Gold sequences;(d) depth maps by m-sequences;(e) theoretical simulation of dead time influence on imaging performance;(f) dead time is 200 ns;(g) dead time is 20 ns;(d) dead time is 1 ns.Colorbars show the depth information with unit of meter.

設第n個像素重建的深度為 Zn;Gn為第n個像素的深度真值;N2 為像素數(shù).該項指標用于衡量成像準確程度,誤差越小,確度越高.

當噪聲mn=10 時,圖4(a)—(d)給出m序列長度分別為2048,4096,8192 和16384 重建的深度圖.隨著序列長度的增加,圖像的信噪比增強,與(11)式一致.圖4(j)給出碼長1024到32768 重建圖像的MSE,32768 碼長的MSE為0.033 m,1024 碼長的MSE 為0.1 m.由于m序列良好的自相關性和ML 的穩(wěn)定性[18],在環(huán)境噪聲和暗計數(shù)噪聲的干擾下,仿真設定基于空間編碼的單像素成像方法的積分時間與本文所提出的16384 碼長的積分時間相同,圖4(d)所示本文提出的成像方法優(yōu)于圖4(i)所示傳統(tǒng)的成像方法,也表明了大多數(shù)單像素激光雷達在遠距離和戶外成像質量不高的原因.

其次,在本文提出的成像方法基礎上,引入二次相關重建方法,重建的深度圖4(e)—(h) MSE 比不采用二次相關重建的深度圖4(a)—(d) MSE 小,準確度高,和圖4(j)一致.傳統(tǒng)的基于空間編碼的單像素成像的MSE 為0.201 m,引入二次相關法后成像的MSE 為0.02 m.本文提出的成像方法的MSE 是傳統(tǒng)成像方法的1/10,成像誤差降低了9/10.最后,采用Sobel算子對圖4(d)、圖4(h)和圖4(i)進行邊緣檢測,結果分別如圖5(a)、圖5(b)和圖5(c)所示.引入二次相關方法后,邊緣檢測的噪聲魯棒性更高,邊界更加平滑和清晰.

圖4 不引入二次相關法的時空域聯(lián)合編碼擴頻單光子計數(shù)深度成像,碼長為 (a) 2048,(b) 4096,(c) 8192,(d) 16384.引入二次相關法的時空域聯(lián)合編碼擴頻單光子計數(shù)深度圖像,碼長為(e) 2048,(f) 4096,(g) 8192,(h)16384;(i) 傳統(tǒng)的基于空間編碼的單像素深度圖像;(j) 引入二次相關法(點虛線)和不引入二次相關法(虛線)的深度MSEFig.4.Time-space united coding image spread spectrum single photon counting imaging without second correlated method by code length of (a) 2048,(b) 4096,(c) 8192,(d)16384.Depth maps simulations with second correlated method by code length of (e) 2048,(f) 4096,(g) 8192 and (h) 16384.(i) Traditional single pixel imaging method based on space coding.(j) The depth MSE with second correlated method (dot dashed line) and without second correlated method (dashed line).

圖5 成像邊緣檢測 (a) 不引入二次相關法的時空域聯(lián)合編碼擴頻單光子成像邊緣檢測;(b) 引入二次相關法的時空域聯(lián)合編碼擴頻單光子計數(shù)成像邊緣檢測;(c) 傳統(tǒng)的基于空間編碼的單像素成像邊緣檢測Fig.5.Image edge detection: (a) Time-space united coding spread spectrum single photon counting imaging without second correlated method image edge detection;(b) time-space united coding spread spectrum single photon counting imaging with second correlated method image edge detection;(c) traditional single pixel imaging method based on space coding image edge detection.

最后,采用本課題組單光子計數(shù)掃描成像系統(tǒng)進行實際場景的模擬測試[27,28].噪聲光子計數(shù)值為100 (counts/s),單個像素積分時間為1 s,成像系統(tǒng)的時間分辨率為8 ps,圖6(a)為采用該系統(tǒng)測量的高精度深度圖,將其作為真實深度圖的模擬.采用本文提出的方法進行目標探測和重建,結果如圖6(b)—(h)所示.由于所采用的系統(tǒng)深度分辨率比圖像“飛機”高,因此,在相同條件下,圖6(b)、圖6(e)、圖6(h)中的MSE 分別比圖4(i)、圖4(d)、圖4(h)中的MSE 小.m序列長度為16384,當mn=10時,傳統(tǒng)的基于空間編碼的單像素深度成像如圖6(b)所示.在不引入二次相關法的情況下,當mn=5 到10 時,時空域聯(lián)合編碼擴頻單光子計數(shù)成像分別如圖6(c)、圖6(d)和圖6(e)所示.時空域聯(lián)合編碼擴頻單光子計數(shù)成像性能優(yōu)于基于空間編碼的單像素成像性能.基于空間編碼的單像素成像的MSE 為0.176 m,而提出的時空域聯(lián)合編碼擴頻單光子計數(shù)成像的MSE 為0.035 m.因此,本文提出方法的MSE 降低了4/5.當mn=5—10 時,采用二次相關法重建的深度圖分別如圖6(f)、圖6(g)和圖6(h)所示,與不采用二次相關法的成像性能相比,二次相關法重建的深度圖MSE 平均減少量約為0.018 m.與圖6(b)傳統(tǒng)空間編碼單像素成像方法相比,引入二次相關法后MSE 降低了9/10,圖7(b)的成像邊緣更清晰,對噪聲的魯棒性更強,相比其他邊緣檢測結果更加平滑.針對(19)式采用TVAL3,GPSR,FISTA 重建的深度圖如圖8所示,重建結果表明TVAL3 的性能最好.

圖6 單光子計數(shù)掃描成像測試圖像深度重建仿真 (a) 玩偶深度真值;(b) mn=10,MSE=0.176 m,傳統(tǒng)的基于空間編碼的單像素成 像深度 圖;在(c) mn=5,MSE=0.019 m,(d) mn=8,MSE=0.03 m,(e) mn=10,MSE=0.035 m 條件下,不引入二次相關法的時空域聯(lián)合編碼擴頻單光子計數(shù)深度圖;在(f) mn=5,MSE=0.005m,(g) mn=8,MSE=0.011 m,(h) mn=10,MSE=0.017 m條件下,引入二次相關法的時空域編碼擴頻單光子計數(shù)深度圖,深度單位為mFig.6.Depth reconstruction simulation by using single photon counting scanning imaging test image: (a) Ground truth of ‘doll’;(b) the depth map by traditional single pixel imaging method based on space coding image when mn=10,MSE=0.176 m .Depth maps by time-space united coding spread spectrum single photon counting imaging without second correlated method when (c) mn=5,MSE=0.019 m,(d) mn=8,MSE=0.03 m and (e) mn=10,MSE=0.035 m .Depth maps by time-space united coding spread spectrum single photon counting imaging with second correlated method when (f) mn=5,MSE=0.005 m,(g) mn=8,MSE=0.011 m,and (h) mn=10,MSE=0.017 m .The depth unit is m.

圖7 成像邊緣檢測 (a) 不引入二次相關法的時空域聯(lián)合編碼擴頻單光子成像邊緣檢測;(b) 引入二次相關法的時空域聯(lián)合編碼擴頻單光子成像邊緣檢測;(c) 傳統(tǒng)的基于空間編碼的單像素成像邊緣檢測Fig.7.Image edge detection: (a) Time-space united coding spread spectrum single photon counting imaging without second correlated method image edge detection;(b) time-space united coding spread spectrum single photon counting imaging with second correlated method image edge detection;(c) traditional single pixel imaging method based on space coding image edge detection.

圖8 噪聲 mn=10 不同的凸優(yōu)化反演方法的圖像重建 (a) TVAL3;(b) GPSR;(c) FISTAFig.8.Different convex optimization inversion methods imaging reconstruction when mn=10: (a) TVAL3;(b) GPSR;(c) FISTA.

表1 中計算了不同噪聲水平下TVAL3,GPSR和FISTA 的均方誤差和重建時間,以充分研究凸優(yōu)化反演方法的性能.TwIST方法的MSE 如圖6所示.噪聲分別為mn=5 和mn=10,TwIST 方法的運行時間分別為20 s 和28 s.FISTA 和GPSR反演的圖像質量幾乎與TwIST 相同,FISTA 和GPSR 運行時間更短.TVAL3 需要對參數(shù)進行更多調整,雖然重建質量高,但耗費時間長.

表1 不同重建方法的性能統(tǒng)計Table 1.Performance statistical record of different reconstruction methods.

4 結論

本文提出了一種結構緊湊、非掃描的時空域聯(lián)合編碼成像架構.基于單光子探測理論,建立了時空域聯(lián)合相關的探測模型、成像正向模型和輸出信噪比模型.理論模型和實驗結果證明,通過配置時域編碼的平衡度可優(yōu)化成像質量.將二次相關法引入到成像模擬中,實驗結果與理論模型相符,均證明該技術具有探測距離遠、抗干擾能力強、準確度高的優(yōu)勢.今后的工作中,課題組將構建時空域聯(lián)合編碼擴頻單光子計數(shù)無掃描成像系統(tǒng).希望該項工作能夠提高機器視覺、遠程探測和無人駕駛等應用的成像性能,為傳統(tǒng)的單光子計數(shù)單像素激光雷達的研究打開新的局面.