沈姗姗 顾国华 陈钱 何睿清 曹青青

1) (南京工业职业技术大学航空工程学院,南京 210023)

2) (南京理工大学电子工程与光电技术学院,智能感知和微光成像实验室,南京 210094)

3) (南京工程学院信息与通信工程学院,南京 211167)

本文结合空间编码单像素成像技术和扩频时间编码扫描成像技术,提出一种时空域联合编码扩频单光子计数成像方法.该方法具备可避免距离模糊、大时宽带宽积的优势,并且在噪声干扰下,能够准确恢复距离像.本文推导了基于单光子探测的时空域联合相关非线性探测模型、成像正向模型和信噪比模型,并通过凸优化反演算法恢复深度图像,理论模型和仿真实验均证明,与传统的基于空间编码的单像素成像方法相比,本方法提高了重建的质量.其中,采用m 序列作为时间编码,成像的噪声鲁棒性更高.和传统的空间编码单像素成像技术相比,本文提出方法的成像均方误差降低了4/5,引入二次相关方法后,成像均方误差降低了9/10.本文所提出的成像架构为非扫描激光雷达成像方法提供了新思路.

1 引言

单光子计数三维成像技术[1−3]广泛应用于机器视觉、遥感和无人驾驶等领域.其中,利用目标场景的可压缩性,从欠采样线性投影序列中重建三维图像是近年来广泛研究的单光子计数三维成像方法.该方法在空间光调制器(digital micro mirror device,DMD)上加载一系列的随机编码.进而在空间域调制光信号,单光子探测器(single photon avalanche detector,SPAD)用于测量场景和随机编码的内积,从远小于像素大小的投影中获取目标信息,并在稀疏约束下根据欠采样测量值重建图像.Howland 等[4]从当前测量的信号中减去前一次测量的参考信号以削弱噪声能量,实现32×32 分辨率高帧速动态目标跟踪.Zhao 等[5]采用Hadamard基提出一种高效的单像素全彩成像方法.在微光环境下,Radwell 等[6]提出了一种基于深度学习的单像素成像优化算法.李凤强等[7]利用L2范数进行约束,并采用TwIST(two-step iterative shrinkage/thresholding)反演方法重建图像,所提出的100 万像素的压缩感知(compressed sensing,CS)激光雷达的成像空间分辨率提高了4 倍.Liu 等[8]提出一种SP3I 技术以提高系统的接收效率和输出信噪比(signal-to-noise ratio,SNR).Asmann 等[9]提 出一种阵列的压缩感知激光雷达系统,成像信噪比逼近全帧互相关方法的成像信噪比,并分析了全变分增强拉格朗日交替方向法(total variation augmented Lagrangian alternation direction algorithm,TVAL3)和梯度投影稀疏重建法(gradient projection for sparse reconstruction,GPSR)的重建结果.传统的基于空间编码的单像素激光雷达通常以固定频率(10 kHz 至10 MHz)发射激光脉冲,限制了最大不模糊距离和信号探测速率.此外,分辨率和带宽之间的矛盾难以调和,成像结果易受噪声影响,使得此类方法成像信噪比低,仅限于室内和短距离成像的应用场合.

扩频通信技术广泛应用于无线传感器网络[10]、全球定位系统和扫描式激光雷达[11−18]等领域,其中,基于扩频通信时间编码的单光子计数扫描成像方法[11−18],通过接收端的匹配滤波器实现脉冲压缩,高效提高输出信噪比并减少信号衰减的概率,具有可避免距离模糊、大时宽带宽积、低功耗和易于单片集成的优势.该方法通常采用伪随机序列进行时间编码,Ding 等[15]采用加权的多距离时间相关函数作为目标函数,通过模拟退火方法优化伪随机序列结构,提高时间相关函数的距离分辨能力.Hiskett 等[16]设计“1”码概率为0.1 的序列,并指出若“1”码概率增加,噪声相关水平也会增加.

本文结合空间编码单像素成像技术和扩频时间编码扫描成像技术,提出了一种基于时空域联合编码的扩频单光子计数三维成像架构.本方法采用伪随机序列在时间域调制激光脉冲,替代激光脉冲的周期性调制方法,并将调制后的光脉冲投影到DMD,在空间域调制激光脉冲.SPAD 用于探测返回的光脉冲信号,最后将接收信号和参考信号互相关,从而建立时空域联合相关非线性探测的数学模型,该模型将场景点多路复用到单光子探测器上,进而将压缩感知凸优化反演技术(如TwIST[19])引入到新的成像架构中,由此,建立成像正向模型和输出信噪比模型,并研究死时间对成像质量的影响规律.仿真实验和理论模型共同研究表明,通过优化时域编码的平衡度,可降低环境噪声和暗计数的干扰,采用更长的序列可以获得更高的信噪比和深度确度.此外,针对一般的时域编码序列,理论模型和仿真实验表明,引入的二次相关重建方法[20,21]可降低噪声影响,提高成像质量.

2 成像理论

2.1 扩频时间编码单光子计数成像方法

扩频时间编码单光子计数成像系统中包括伪随机序列发生器,其用于生成0,1 码流,码为1 时触发激光器发射光脉冲,激光脉冲经目标后漫反射得到回波信号,经由光学系统接收,到达SPAD 感光面.SPAD 启动一次探测并驱动光子到达时间记录仪(time-to-digital converter,TDC)产生时间值.经过上述发射接收过程,完成单个像素的有效探测.定义成像场景的反射率为r(x,y),激光脉冲由SPAD 和TDC 接收后,通过匹配滤波器进行脉冲压缩,也即将接收序列s(t) 与发送副本序列f(t)互相关,则像素 (x,y)th的互相关信号为

其中T表示一个周期的序列长度;PR为信号光功率;∆t为码元宽度.选取适当的距离门信号可以消除后向散射峰值,(1)式中第一个最高的互相关峰值对应的横坐标τ(x,y)为从场景目标反射的光子到达时间的估计值.因此,可采用最大似然估计法(maximum likelihood estimation,ML)[14]估 计τ(x,y)ML:

目标的距离值为

其中c为光速.

2.2 时空域联合编码扩频单光子计数成像方法

本文提出的时空域联合编码扩频单光子技术成像架构如图1 所示.采用基于可编程逻辑器件的序列发生器发送扩频伪随机序列.并驱动激光器发射调制后的激光脉冲,在时间域完成激光脉冲的编码;采用空间光调制器调制返回的激光脉冲,在空间域完成激光脉冲的编码.采用偏振分束镜将从空间光调制器反射的激光光束与从目标反射的激光光束分开,透镜收集时空域联合编码的激光光束并导入SPAD,信号处理器记录光子到达时间点,并且和模板序列做互相关运算,通过后续信号处理重建目标的三维图像.设先后被第i次时域扩频伪随机序列和空域随机矩阵调制后的光束,依次通过偏振分束镜和透镜,到达单光子探测器的总光子计数率记为Ri(t) :

图1 时空域联合编码扩频单光子计数成像架构Fig.1.Time-space united coding spread spectrum single photon counting imaging method architecture.

式中N表示x轴或y轴方向的像素个数;Φi(x,y) 表示第i次调制的随机矩阵;Bi1(x,y,t) 表示空域外部环境噪声光子计数率;Bi2(t) 表示时域外部环境噪声光子计数率.第i次时空域联合调制后的光束被SPAD 探测后产生响应,该过程为非齐次泊松过程.假设探测器的量子效率为η,则探测到的光子计数率为ηRi(t) .探测器暗电流产生响应的过程为齐次泊松过程,d表示暗电流响应速率,可得探测到的光子计数率为

采用相关接收技术,将(5)式和模板序列f(t) 互相关,时空域联合相关的非线性探测模型为

(6)式中,第一项为时空相关信号,外部环境噪声Bi1(x,y,t) 和Bi2(t)与模板序列互相关后分别得到噪声项暗计数噪声d和模板序列互相关得到噪声项Bd(τ) .

测量向量为C ∈RM×1,空间光调制器加载的测量矩阵为Φ ∈RM×Q,投影次数为M,噪声随机向量为b ∈RQ×1,像素总数为Q=N2.(6)式表明Ci(τ)包括信号光子Cis和噪声光子Cin,本文通过给出Ci(τ)中变量的统计特征来建立信噪比模型.在微光环境下,单光子探测服从泊松过程[22],假设噪声和信号是独立的,且期望和方差相等.则第i次调制得到的Ci(τ) 中包含的信号光子Cis的数学期望为

(11)式表明输出信噪比与序列长度成正比.分母的后三项均包括f(t+τ),如果序列中“1”码的概率几乎等于“–1”码的概率,则噪声能量可忽略不计.因此,通过合理配置时域编码序列正负比特的平衡性,可降低噪声能量.其次,考虑探测器的死时间,引入探测到的信号光子或噪声光子的总概率,死时间增大,雪崩触发概率降低,探测到信号光子的总概率降低,输出信噪比降低.

2.3 二次相关重建方法

根据(6)式,互相关函数Ci(τ) 中包含信号和噪声项,其中,第一项是信号项.只有当目标平面和调制平面之间在空间上满足严格的一对一相关性,以及接收到的光子到达时间点和发送的副本序列之间在时间上满足严格的一对一的相关性,才能生成无噪声干扰的信号项.噪声主要来自系统和环境,这两类噪声不相关,将两者相乘不会生成直流分量.因此考虑通过二次相关来降低噪声.假设矩阵Φ的每一列表示为φz ∈RM×1,M=1,2,⋅⋅⋅,Q,则加性噪声b可以等效为随机矩阵V与信号p的乘积,表达为b=V p,其中V与Φ非相关.对(7)式做二次相关:

其中∆∈RQ×1,H=cov(Φ,Φ) .压缩感知充分利用信号的可压缩性,通过优化的方法恢复稀疏表示的测量值.如果p在稀疏基Ψ的表示下是足够稀疏的,则有p=ΨZ,其中Z是稀疏系数.(18)式等效为∆=HΨZ=AZ,其中A为感知矩阵且A=HΨ,采用以下凸优化反演方法重建图像

其中λ是正则化参数;κ(Z) 是正则项.实验中,采用TwIST[19]方法求解(19)式.为了平衡计算的鲁棒性和复杂性,引入最大似然估计(ML)来优化重建.

3 实验结果与分析

为了验证本文提出成像方法的有效性,采用MATLAB 进行成像仿真实验.由MATLAB 生成10 GHz 伪随机序列驱动激光器发射激光脉冲,伪随机序列中1 码的概率由随机数阈值确定,设置为0.1,长度为2048.灰度图像“飞机”的深度范围为20—20.8 m,如图2(a)所示.将生成的M个64×64 大小的高斯矩阵依次加载至DMD 上,接收端数字化时间编码和空间编码的信号后,结合TwIST 和ML 法重建图像.图2(b)—(g)给出采用不同压缩比(5%—85%)进行重建的结果,压缩比表示压缩的程度,定义为,压缩程度越高,测量次数越少,重建速度越快,但重建质量越差.采用不同压缩比探测和重建图像的总时间在8—78 s 之间.图2(e)—(g)中45%或更大的压缩比可实现高分辨率和高对比度的重建.由于ML 的稳定性,本方法的重建结果与参考文献[7,24]中强度图像重建工作所得出的结论一致,即25%的压缩比也可清晰呈现目标场景,如图2(d)所示.M序列和Gold 序列是CDMA 系统上行链路中最常用的伪随机序列[25,26].为了比较不同序列的抗噪声能力,在仿真中引入泊松噪声.平均噪声光子计数之和定义为三种噪声光子计数值之和,记为mn.假设平均信号光子计数ms不变且ms=5.噪声干扰由泊松随机数发生器生成,将mn=1的噪声引入成像仿真,由图3(a)和图3(b),1 码概率为0.3 的伪随机序列成像效果优于1 码概率为0.1 的序列成像效果.由于m序列中1 码的概率和–1 码的概率基本相同,图3(d)中所示的m序列的成像效果优于其他序列的成像效果,与(11)式序列平衡度和成像信噪比的关系相互验证.考虑单光子探测器死时间,根据(11)式,数值模拟死时间和成像信噪比的制约关系,如图3(a)所示.随着死时间的增大,信号光子探测概率降低,信噪比降低;相同死时间下,码流速度越大,在死时间内到达探测器的信号光脉冲数量越多,探测到的信号光子数目减少,信噪比降低.16384 码长下1 码概率为0.3,引入不同死时间后的成像仿真结果如图3(e)—(h)所示,实验结果表明死时间减小,信噪比增大,成像质量提高,变化规律和理论模型数值模拟结果一致.定义深度图像重建质量的均方误差(mean squared error,MSE)为

图2 不同压缩比的时空域联合编码扩频单光子计数成像 (a) “飞机”深度真值;(b)−(g)不同压缩比的深度成像 (b) 5%,(c)15%,(d) 25%,(e) 45%,(f) 65%,(g) 85%.颜色条表示深度数值,单位为mFig.2.Time-space united coding spread spectrum single photon counting image with different compression proportions: (a) Depth ground truth of ‘airplane’;(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 时空域联合编码扩频单光子计数成像性能模拟 (a) 1 码概率为0.1 的伪随机序列深度图像;(b) 1 码概率为0.3 的伪随机序列深度图像;(c) Gold 序列的深度图像;(d) m 序列的深度图像;(e) 死时间对成像性能的影响理论模拟;(f) 死时间为200 ns;(g) 死时间为20 ns;(h) 死时间为1ns.颜色条表示深度数值,单位为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 为像素数.该项指标用于衡量成像准确程度,误差越小,确度越高.

当噪声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 的稳定性[18],在环境噪声和暗计数噪声的干扰下,仿真设定基于空间编码的单像素成像方法的积分时间与本文所提出的16384 码长的积分时间相同,图4(d)所示本文提出的成像方法优于图4(i)所示传统的成像方法,也表明了大多数单像素激光雷达在远距离和户外成像质量不高的原因.

其次,在本文提出的成像方法基础上,引入二次相关重建方法,重建的深度图4(e)—(h) MSE 比不采用二次相关重建的深度图4(a)—(d) MSE 小,准确度高,和图4(j)一致.传统的基于空间编码的单像素成像的MSE 为0.201 m,引入二次相关法后成像的MSE 为0.02 m.本文提出的成像方法的MSE 是传统成像方法的1/10,成像误差降低了9/10.最后,采用Sobel算子对图4(d)、图4(h)和图4(i)进行边缘检测,结果分别如图5(a)、图5(b)和图5(c)所示.引入二次相关方法后,边缘检测的噪声鲁棒性更高,边界更加平滑和清晰.

图4 不引入二次相关法的时空域联合编码扩频单光子计数深度成像,码长为 (a) 2048,(b) 4096,(c) 8192,(d) 16384.引入二次相关法的时空域联合编码扩频单光子计数深度图像,码长为(e) 2048,(f) 4096,(g) 8192,(h)16384;(i) 传统的基于空间编码的单像素深度图像;(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) 不引入二次相关法的时空域联合编码扩频单光子成像边缘检测;(b) 引入二次相关法的时空域联合编码扩频单光子计数成像边缘检测;(c) 传统的基于空间编码的单像素成像边缘检测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.

最后,采用本课题组单光子计数扫描成像系统进行实际场景的模拟测试[27,28].噪声光子计数值为100 (counts/s),单个像素积分时间为1 s,成像系统的时间分辨率为8 ps,图6(a)为采用该系统测量的高精度深度图,将其作为真实深度图的模拟.采用本文提出的方法进行目标探测和重建,结果如图6(b)—(h)所示.由于所采用的系统深度分辨率比图像“飞机”高,因此,在相同条件下,图6(b)、图6(e)、图6(h)中的MSE 分别比图4(i)、图4(d)、图4(h)中的MSE 小.m序列长度为16384,当mn=10时,传统的基于空间编码的单像素深度成像如图6(b)所示.在不引入二次相关法的情况下,当mn=5 到10 时,时空域联合编码扩频单光子计数成像分别如图6(c)、图6(d)和图6(e)所示.时空域联合编码扩频单光子计数成像性能优于基于空间编码的单像素成像性能.基于空间编码的单像素成像的MSE 为0.176 m,而提出的时空域联合编码扩频单光子计数成像的MSE 为0.035 m.因此,本文提出方法的MSE 降低了4/5.当mn=5—10 时,采用二次相关法重建的深度图分别如图6(f)、图6(g)和图6(h)所示,与不采用二次相关法的成像性能相比,二次相关法重建的深度图MSE 平均减少量约为0.018 m.与图6(b)传统空间编码单像素成像方法相比,引入二次相关法后MSE 降低了9/10,图7(b)的成像边缘更清晰,对噪声的鲁棒性更强,相比其他边缘检测结果更加平滑.针对(19)式采用TVAL3,GPSR,FISTA 重建的深度图如图8所示,重建结果表明TVAL3 的性能最好.

图6 单光子计数扫描成像测试图像深度重建仿真 (a) 玩偶深度真值;(b) mn=10,MSE=0.176 m,传统的基于空间编码的单像素成 像深度 图;在(c) mn=5,MSE=0.019 m,(d) mn=8,MSE=0.03 m,(e) mn=10,MSE=0.035 m 条件下,不引入二次相关法的时空域联合编码扩频单光子计数深度图;在(f) mn=5,MSE=0.005m,(g) mn=8,MSE=0.011 m,(h) mn=10,MSE=0.017 m条件下,引入二次相关法的时空域编码扩频单光子计数深度图,深度单位为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) 不引入二次相关法的时空域联合编码扩频单光子成像边缘检测;(b) 引入二次相关法的时空域联合编码扩频单光子成像边缘检测;(c) 传统的基于空间编码的单像素成像边缘检测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 不同的凸优化反演方法的图像重建 (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 的均方误差和重建时间,以充分研究凸优化反演方法的性能.TwIST方法的MSE 如图6所示.噪声分别为mn=5 和mn=10,TwIST 方法的运行时间分别为20 s 和28 s.FISTA 和GPSR反演的图像质量几乎与TwIST 相同,FISTA 和GPSR 运行时间更短.TVAL3 需要对参数进行更多调整,虽然重建质量高,但耗费时间长.

表1 不同重建方法的性能统计Table 1.Performance statistical record of different reconstruction methods.

4 结论

本文提出了一种结构紧凑、非扫描的时空域联合编码成像架构.基于单光子探测理论,建立了时空域联合相关的探测模型、成像正向模型和输出信噪比模型.理论模型和实验结果证明,通过配置时域编码的平衡度可优化成像质量.将二次相关法引入到成像模拟中,实验结果与理论模型相符,均证明该技术具有探测距离远、抗干扰能力强、准确度高的优势.今后的工作中,课题组将构建时空域联合编码扩频单光子计数无扫描成像系统.希望该项工作能够提高机器视觉、远程探测和无人驾驶等应用的成像性能,为传统的单光子计数单像素激光雷达的研究打开新的局面.