1 引言

激光反射断层成像(LRT)是近年来发展起来的一种新型激光成像方式。传统光学镜头存在衍射极限,不适用于远距离目标成像,而激光成像雷达能克服衍射极限,且受外界环境的影响较小^[1]。开展LRT回波信号的建模仿真与分析研究对LRT目标重构技术的发展具有重要意义。Parker等^[2-3]在计算机断层成像(CT)中引入激光反射层析,并进行了实验验证。Matson等^[4-5]研究了LRT的理论和应用。林放等^[6-8]在特征点跟踪、二维成像、图像配准、脉冲压缩解卷积等方面进行了实验研究。

国内外学者在利用激光束进行目标探测方面也进行了大量的研究,其中三维(3D)激光雷达仿真系统主要包括瑞典**研究局研发的3D成像仿真模型^[9-10]和美国林肯实验室为激光成像雷达“竖锯”开发的雷达仿真系统^[11]。目前,已开发的完整可实际应用的仿真软件有美国洛克希德·马丁公司的ALASS^[12]和美国犹他州立大学研发的USU Ladar SIM^[13]。国内学者也在一定程度上开展了激光成像雷达系统的建模仿真研究^[14-16],金晓峰等^[17]提出了反射层析激光雷达小系统的成像模拟,但其在复杂轮廓目标的应用上有待进行进一步深入研究。

根据LRT的基本原理,要实现目标二维轮廓图像重构的建模仿真,需要着重解决激光子光束与探测目标表面的交点获取问题及激光子光束与不同特性的目标表面作用过程的建模问题。Carlsson等^[9-12]采用从目标向成像平面作垂直映射的方法,这种方法无法精细观察单个激光束的作用过程,且与激光束和探测目标之间的实际作用过程不符。当目标处于运动状态时,Graham等^[18-19]主要采用光线追踪算法获取了激光子光束与探测目标的交点,其目标场景模型采用的是体素表示法,这种方法没有对扩束激光进行细分。因此,如何快速有效地获取激光子光束与探测目标表面的交点是迫切需要解决的一个问题。

本文结合激光雷达探测的基本原理,展开了对基于LRT的目标重构建模仿真。研究了激光子光束与探测目标表面的交点获取问题,准确获取了交点处的距离、反射率、入射角等相关信息;研究了激光子光束与不同特性的目标表面作用过程,得到了激光照射目标一周的回波数据。运用滤波反投影算法完成了目标二维轮廓的图像重构,重构结果验证了仿真模型的正确性,为进一步开展复杂场景下LRT的目标重构和空间目标探测提供了参考。

2 LRT原理

LRT的基本重建算法是滤波反投影(FBP)算法,其数学原理是傅里叶切片定理。如图1所示,平行激光光束照射一个二维目标,照射角度记为ϕ,目标在角度ϕ的投影定义为

p(r,ϕ)=∫f(x,y)ds,(1)

图 1. LRT的示意图

Fig. 1. Schematic of LRT

下载图片 查看所有图片

式中(x,y)为目标在Oxy坐标系中的位置;积分路径L_r,ϕ为垂直于光照方向的直线;极径r=xcosϕ+ysin ϕ;f(x,y)为目标的反射系数分布;s为积分路径。平行箭头的展开区域为激光照射区域,投影的距离范围取决于被激光照射区域的距离范围。激光无法穿透目标,实际上反射系数只在表面不为零,即

f(x,y)=0, (x,y)∉D,(2)

式中D为目标表面点集。因此,实际投影可表示为

p(r,ϕ)=∫f(x,y)ds。(3)

目标重构过程与投影数据的获取是逆向的过程,需要对采集的投影数据进行预处理以去除掉信号采集过程中的噪声。投影数据的配准是将采集的回波信号按照统一的时间基准进行重新排列。接受回波信号是目标表面的反射率分布信息与激光脉冲的卷积,故在目标重构前需要对接受信号进行解卷积运算以得到目标的实际反射率分布投影数据。将全部探测角度下的解卷积后投影数据进行直接反投影运算,可得到物体重建的二维轮廓图像。

直接反投影重构算法的主要思想是将每一次测得的投影数据按角度为ϕ的扫描路径“原路”回抹到路径经过的像素上。即是将所有经过这个坐标点的值进行累加或累加后取均值,再将该值作为重建的二维轮廓图像的该像素点的值。但这样会产生伪影,影响图像分辨,因此通常使用滤波反投影算法,重建的目标反射系数分布的图像表示为

gFB(x,y)=∑i=1mF1-1ωF1[p(r,ϕ)],(4)

式中m为投影数量;F₁、 F1-1分别为一维傅里叶变换和傅里叶逆变换的算符; ω为由直角坐标系转到极坐标系产生的雅可比因子,起斜坡滤波器的效果。

3 仿真系统建模

3.1 目标的投影角度旋转

在LRT的目标轮廓重构过程中,需要采集环绕目标360°的全角度信息,根据目标与激光成像雷达的相对运动关系可确定探测雷达的位置,通过将目标旋转不同的角度来采集不同角度的目标回波信息。在LRT过程中,根据环绕目标旋转一周需要采集的回波信号的组数调整目标的单次旋转角度。

3D空间中的旋转变换比二维空间中的旋转变换复杂,除了需要指定旋转角外,还需指定旋转轴。若以坐标系的三个坐标轴x、y、z分别作为旋转轴,则点实际上只在垂直坐标轴的平面上作二维旋转,此时用二维旋转公式即可直接推出3D旋转变换矩阵。规定在右手坐标系中,物体旋转的正方向是右手螺旋方向,即从该轴正半轴向原点看为逆时针方向(俯视图角度)。绕z轴正向旋转γ角后的坐标为

(x'y'z')=(xyz)cosγsinγ0-sinγcosγ0001。(5)

旋转后点的z坐标值保持不变,x、y坐标的变化相当于在xOy平面内作正γ角的旋转。

物体绕平行于某一坐标轴的旋转变换的基本步骤为:1)平移物体使旋转轴与所平行的坐标轴重合;2)沿该坐标轴进行指定角度的旋转;3)平移物体使旋转轴移回到原位置。

3.2 扩束激光在目标场景中的投影点获取

在LRT目标重构的过程中,激光源发射的激光束为一束扩束光束,一定距离处的光斑直径明显大于目标的尺寸。当扩束激光照射目标表面时,激光光斑往往能完全覆盖目标,在每个激光束中,由于目标表面的空间和反射率分布的特性,不同部分对激光束的作用过程存在差异。为了准确反映激光扩束光束与目标表面不同部分不同特性的作用过程,需要运用微分的思想将一束宽束激光分解为一定数量的子光线来分别求取每条子光线与目标的交点,这样就能得到扩束激光与目标不同部位的作用过程,最终不同部位的子光束的回波在时间域上的线性叠加即为扩束激光的回波信号。将一束激光细分为N×N条子光线,如图2所示,通过计算每条子光线与目标场景模型最近的交点,即能通过交点处的目标表面信息得到每个微小面元的回波信号。

图 2. 扩束激光的子光线点阵

Fig. 2. Dot matrix diagram of sub-rays of expanded laser

下载图片 查看所有图片

3.2.1 需求分析

通过3DSMAX软件绘制探测目标模型的3D形状,并设置其在目标场景坐标系中的方位、缩放比例、反射率等,再将绘制完成的模型以3DS的文件形式导出,通过Visual Studio软件编写的相关程序读取目标场景模型的三角形面元各顶点的坐标等相关信息。故激光束与探测目标的求交问题实际上相当于每个子光线与三角形面元求交点问题。直线与多边形求交点是计算机图形学中的一个基础问题,一般通过联立方程组求解^[20-21]。在实际建模过程中,为了能够尽可能贴近实际的探测目标,模型的面元数量往往十分庞大,同时,为了能在建模过程获得准确的回波数据,激光束分解的子光线数量也较多,每条子光线都需要与目标面元进行求交运算。假如目标场景模型由1000个面元组成,激光束细分为30×30个子光线,获取激光束的投影点就需进行9×10⁵次求交运算,计算量巨大,因此必须解决光线求交的计算效率问题。

3.2.2 基于激光子光线面阵的快速求交算法

扩束激光的子光线集只与探测目标中的很少一部分面元相交,若能在每条子光线与目标场景求交点前先初步判断哪些面元可能与当前激光束内的子光线相交,就能在一定程度上提高光线求交的效率。

图 3. 光束平面与面元位置关系示意图

Fig. 3. Positional relationship between beam plane and element

下载图片 查看所有图片

如图3所示,将一束激光细分为N×N条子光线点阵,取点阵每一列中的最边缘的两个边界点P₁、P₂和激光源所在位置O点构成一个三角形OP₁P₂。探测目标中的三角形面元与激光束的一列点阵所在的平面可能存在以下几种位置关系:

1) 三角形面元的三个顶点全部处在平面的左侧,如图3中的A面元;

2) 三角形面元的三个顶点全部处在平面的右侧,如图3中的C面元;

3) 三角形面元与平面相交,这时面元的顶点分布在子光线一列点阵形成的平面的两侧,如图3中的B面元。

在第一种和第二种情况下,三角形面元完全处于激光束所在平面的范围外,与激光子光线不存在交点,可以将其剔除。光线求交的具体步骤如下。

1) 激光源所在的坐标O点和P₁、P₂所产生的两条边界光线的方向矢量F₁、F₂确定的激光一列子光线所在的平面方程为

a1x+b1y+c1z=d1,(6)n1=F1×F2=[a1b1c1]T,(7)

式中a₁、b₁、c₁为系数;T代表求转置;d₁为平面方程常数项;n₁为平面方程的法向量。

2) 确定点( xiyizi)^T(i为)与平面a₁x+b₁y+c₁z=d₁的位置关系。当点位于包含法线向量的半空间(正半空间)时,存在

a1xi+b1yi+c1zi-d1>0。(8)

当点位于另一侧时,存在

a1xi+b1yi+c1zi-d1<0。(9)

三角形面元的三个顶点的坐标矩阵为

M1=x1y1z1x2y2z2x3y3z3。(10)

故判定矩阵M为

M=M1·n1-D=x1y1z1x2y2z2x3y3z3·a1b1c1-d1d1d1,(11)

式中D为平面方程常数项列矩阵。

根据(11)式可以判断该三角形面元的三个顶点与激光束一列点阵所在平面的位置关系。

3) 判断三角形面元与激光束所在平面的位置关系。M为3×1矩阵,若矩阵M中的元素全部大于零,该面元位于平面的正半空间;若矩阵M中的元素全部小于零,该面元位于平面的负半空间。这两种情况下面元与激光的子光线不存在交点。当矩阵M中的元素即存在正数又存在负数时,该面元与激光束所在的平面相交,此时三角形面元保留。

4) 对于保留下来的面元,按照经典光线求交算法求得每条子光线与面元的交点。

5) 分别取每一列点阵进行N次循环,直至完成所有子光线与探测目标的求交运算。

3.3 激光束与目标场景的作用过程建模

3.3.1 激光发射脉冲的时空分布

在激光脉冲与目标场景作用过程的建模仿真中,影响回波信号的时间延迟和强度的因素主要包括交点距激光源或光电探测器的距离、探测目标的表面反射率及子光线相对于探测目标表面法线的入射角。激光雷达信号的时空分布函数^[22]为

U(x,y,z)=I0·p(t)·I(x,y),(12)

式中p(t)为发射激光脉冲的时间分布函数;t为时间;I₀为单个激光发射脉冲的能量;I(x,y)为扩束激光的空间分布函数。因此,可对扩束激光脉冲与探测目标的作用过程分别从时间域和空间域进行考虑。

p(t)一般采用较为简单的高斯分布模型^[22],即

p(t)=2T12ln2πexp·-4ln2×t-3/2×T12T122,(13)

式中T_1/2为激光脉冲的半峰全宽。

I(x,y)为光束瞬时功率沿激光照射的射线传播方向的空间分布,通常情况下可表示为径向高斯分布,即

I(x,y)=I(r)=ω0ω2exp-2r2ω2=ω0ω2exp-2x2+y2ω2,(14)

式中ω₀为激光源处的光束半径;ω为径向距离R处的光束半径。

3.3.2 目标对激光雷达信号的响应函数

激光雷达的距离方程通常用来计算特定条件下激光雷达接受的回波能量,其公式为

Pr=4IPtTa1ηtπθt2R12ΓTa24πR22πDr2ηr4,(15)

式中P_r为回波信号的功率(接受功率);P_t为激光器的发射功率;I为光束分布函数;T_a1为发射机到目标的大气透过率;η_t为发射光学系统的透过率;θ_t为激光束的发散角;R₁为发射机到目标的距离;Γ为目标的激光雷达截面;T_a2为目标到接收机的大气透过率;R₂为目标到接收机的距离;D_r为接受空间(光学镜头的直径);η_r为接受光学系统的透过率。

对于稳定的收发合置或当发射系统与接受系统十分靠近时,有R₁=R₂=R;若探测环境为几乎真空的太空环境,T_a1和T_a2的数值为1。故(15)式可以简化为

Pr=KPtTa2ΓDr2ηtηr4πθt2R4。(16)

雷达截面的定义^[23]为

Γ=4πR2IrIi,(17)

式中I_r为反射回波在探测器处的能流密度;I_i为入射到目标上的能流密度。粗糙表面单元照射面积dS的截面dΓ和双向反射分布函数f_r之间的关系为

dΓ=4πfrcosθrcosθidS,(18)

式中θ_r为反射角;θ_i为入射角。对于朗伯体来讲,f_r=ρ/π为常量,其中ρ为半球反射率(通常简称为反射率)。若激光雷达系统采用收发合置的探测系统,则θ_i=θ_r=θ;故(18)式可替代为

dΓ=4ρcos2θdS。(19)

考虑当前激光束中的一条子光线r_ray(x,y)(x=1,…,N;y=1,…,N),定义整个探测目标对应的观测区域的每个子光线交点到探测器的距离为R(x,y),每个子光线交点处的目标反射率分布为I_inten(x,y),子光线的入射角为θ_angle(x,y)。根据激光雷达距离的方程式,对应于每个子面元的观测距离和反射率,子面元对扩束激光脉冲的响应函数为

hp(x,y,t)=Ta2Dr2ηtηrI(x,y)Iinten(x,y)πθt2R4cos2θ·δt-2R(x,y)c,(20)

式中δ(t)冲击函数。

3.3.3 激光雷达的全波形回波信号

探测目标上每一个面元到探测器的距离不同,故探测器接受到每个部分的回波信号时间也不同。根据实验常用的雪崩光电二极管(APD)的响应特性可知,到光电传感器距离相同的子面元的响应时间是同一时刻,这个过程对于探测器而言是一个回波信号强度的叠加过程,如图4所示,由此可得到一个时间序列信号,即整个目标的响应函数h_p(t)为

hp(t)=∑Sobject(x,y)hp(x,y,tj),(21)

式中S_object(x,y)为激光束所有子面元的集合;t_j(j=1,2,…)为探测器接受子面元反射回波信号的时间。

图 4. 子光线回波信号叠加的示意图

Fig. 4. Schematic of superposition of sub-ray echo signals

下载图片 查看所有图片

对当前激光束对应的探测目标的响应函数h_p(t)与激光雷达的响应函数及信号的时间分布函数p(t)进行卷积运算,即可得到激光光束与探测目标作用后返回激光雷达接收器的全波形的回波信号 Prreceiver(t)为

Prreceiver(t)=hp(t)*I0p(t)=I0∫-∞+∞hp(τ)p(t-τ)dτ,(22)

式中*表示卷积运算;τ为积分变量。

图5所示为激光信号与目标作用前后的脉冲波形。

图 5. 激光信号与目标作用前后的脉冲波形。(a)发射脉冲信号;(b)回波信号

Fig. 5. Pulse shapes before and after interaction between laser signal and target. (a) Emitted pulse signal; (b) echo signal

下载图片 查看所有图片

4 空间目标重构的仿真实现与分析

基于LRT的回波信号仿真系统,运用3DSMAX软件生成探测目标飞机的1∶1的3D立体模型,其中飞机翼展约为13 m,长约为19 m,机身高度约为2.8 m。设置激光雷达距离目标坐标原点约为10 km,激光束中心相对于目标场景的俯仰角为0°,激光束的发散角为1 mrad,激光束分解为30×30的子光束线阵。通过计算模型上28245个三角形面元和子光线的作用过程,生成全波形检测条件下激光雷达探测的回波数据。

图6(a)、(b)、(c)分别为全波形检测条件下目标的0°、60°、90°三个角度下获得的激光雷达探测目标的回波点阵图。全波形检测条件下每个角度下均可能出现多个探测点,不方便用角度-角度-距离像表述,这里均是经过坐标转换后的3D点云图。实际数据采样需要采集探测目标360°的反射回波信号,这里仅展示部分。

图 6. 不同角度目标的投影点信息。(a) 0°;(b) 60°;(c) 90°

Fig. 6. Projection point information of target under different angles. (a) 0°; (b) 60°; (c) 90°

下载图片 查看所有图片

在仿真实验中共进行了两组对比实验,分别为360°采样角度不同采样间隔条件下的回波数据以及相同采样间隔不同采样角度下的回波数据。在360°采样角度下,以1°、2°、4°间隔采集的360、180、90组数据运用滤波反投影算法重构目标的二维轮廓图像,结果如图8(a)、(b)、(c)所示,探测目标为F15战斗机模型,成像分辨率为2500 pixel×2500 pixel。由成像结果可知,在1°、2°采样间隔条件下目标的外围轮廓像比较清晰,当目标的采样间隔为4°时会有模糊。通过不同采样间隔的仿真实验能够得到LRT不同采样间隔对成像的影响,随着采样间隔的增大,图像会出现模糊,这也在一定程度上为实际应用提供了参考。在1°采样间隔条件下,进行了360°、180°、120°、90°不同采样角度的目标重构实验,遮挡区域的目标特性无法被脉冲捕捉,在相应角度会出现部分信息缺失。目标的细节存在一个可探测角度范围,因此遮挡区域对区域细节探测结果的影响取决于可探测角度范围的大小。图9所示为不同探测角度范围对目标重建图像的影响,可以看出,在大于180°的采样角度下可以收集到较完整的投影信息,当探测角度减小到120°时,目标出现模糊。战机的机翼部分存在重建模糊不全现象,这是由于仿真数据采集的角度为激光束垂直照射战机的侧面,而相对机身部分从激光探测角度来看,机翼的厚度十分纤薄,这会导致数据采集稀少,重构的目标二维轮廓在机翼部分出现模糊。

图 7. 3DSMAX软件绘制的目标原型。(a)俯视正面;(b)仰视背面

Fig. 7. Target prototype drawn by 3DSMAX software. (a) Top view of front side; (b) down view of back side

下载图片 查看所有图片

图 8. 不同角度间隔下采样的数据重构目标图像。(a) 1°;(b) 2°;(c) 4°

Fig. 8. Reconstructed target images based on data sampled under different angle intervals. (a) 1°; (b) 2°; (c) 4°

下载图片 查看所有图片

图 9. 不同探测角度下的目标重构图像。(a) 180°;(b) 120°;(c) 90°

Fig. 9. Reconstructed target images under different probing angles. (a) 180°; (b) 120°; (c) 90°

下载图片 查看所有图片

5 结论

针对LRT仿真系统方面研究存在的空缺问题,从激光子光线与探测目标的投影点相交、激光雷达信号与目标的表面作用过程及采样数据的滤波反投影重构三个方面系统研究了LRT过程中的激光发射、接受及数据处理的主要流程。提出了子光线与探测目标的快速求交算法;基于激光雷达的距离方程建立了激光束与目标场景的响应函数模型。在仿真实验中进行了不同采样间隔和不同探测角度的目标重构实验,初步得到了回波数据采集条件对目标重构的影响规律。