飞行时间成像系统的数字仿真技术*

卢纯青,刘兴潭,武延鹏,杨孟飞

0 引 言

基于连续调制光信号的飞行时间技术是近年来三维成像感知技术的重要发展方向之一，在自动驾驶、SLAM等场景得到逐步应用[1-2].该技术根据幅度调制光信号从光源传播到场景，经过场景表面反射返回传感器的时间差实现对照射物面的三维数据感知.飞行时间成像数字仿真技术是研究将上述物理过程建模，对预设的场景条件和传感器参数建立起可形成实时深度数据集成的数字运行环境，对研究飞行时间传感器在多种光照条件下的成像性能评估和基于三维数据的机器视觉算法的仿真验证具有重要的意义.

目前飞行时间成像仿真的实现方法分为基于特定噪声模型的三维数据仿真[2-3]和基于光信号传播物理过程的三维数据仿真[4-5]两大类，其中基于光信号传播物理过程的三维数据仿真以其数据-参数耦合程度高的特点而具有较大的优势[7].通过对上述两类仿真方法的研究，两类仿真方法为研究传感器参数集与三维成像之间的耦合关系提供了特定的平台，但并未考虑环境光与测量数据间的非定常耦合关系，而在本文研究的有关小行星附着降落的应用场景中，当前的环境光照条件是影响着陆敏感器设计中需要重点考虑的关键因素之一.因此本文以小行星附着着陆任务需求为背景，研究并建立了带有环境光照模型的仿真环境，验证了包括环境光模型的飞行时间成像仿真方法，对于飞行时间成像技术在光照条件下的成像性能的评估具有重要意义和应用前景.

1 飞行时间成像原理

飞行时间传感器基于光信号从光源传播到场景，并返回到传感器的时间t测量距离.设光源是与传感器位置相同的点光源，光在场景和传感器之间传播，则测量距离d为[7]

(1)

幅度调制的飞行时间成像系统的光源发射波长为850 nm±30 nm的光信号，传感器的像元以电信号作为测量信号s和参考信号g的相关函数[3-4]

(2)

假设光源的驱动信号为调制周期为fmod的正弦函数，则相关函数呈现周期函数的形式，有[3-4]

(3)

a为函数幅值，b为函数偏置量，φ为距离相关的相移量，有[3-4]

(4)

一般使用四象限采样法对相关函数进行采样，设采样值为Ii,则有[3-4]

(5)

通过采样值计算相移量φ[3-4]

(6)

相关函数的采样值Di具有二维矩阵的形式，又称为相位图.每个像元的采样值Di由像元电路的两个读出信号NA,i和NB,i计算得到[4]

Ii=NA,i-NB,i

(7)

像元接收光信号产生光生电荷后，由参考信号g控制的电场梯度，导致两个势阱A和B中的光生电荷分布变化.NA,i和NB,i分别为两个势阱读出的电信号.

2 飞行时间成像物理过程模型

假设调制光源L和传感器C位于空间中的相同位置，该位置为坐标原点，场景表面的反射模型为红外波段的朗伯反射，像元的读出电路模型为双读出电路型.

2.1 光源模型

设飞行时间成像的光源为点光源，具有以下3个参数：

1)光源功率PL,物理单位为W；

(8)

对于各向同性光源，立体角ωL为4π.对于具有孔径角ψL的光源，立体角由球冠的面积得到，在单位球面上的高度h为[8]

(9)

则有[3-4]

(10)

2.2 反射模型

(11)

入射角[8]

(12)

(13)

设A为场景表面无限小的面元，则有[8]

(14)

(15)

场景表面P点的辐照度EP为[8]

(16)

设场景表面为反照率系数ρ∈[0,1]的朗伯体，则P点的辐射出射度为[8]

(17)

由于朗伯辐射体在各方向上的辐射相等，因此[8]

(18)

因此P点到传感器的辐射LP为[8-9]

(19)

图1 场景表面反射过程示意图Fig.1 Schematic diagram of the surface reflection process of the scene

2.3 透镜模型

渐晕是指由于离轴光线的有效孔径光阑减小，导致在边缘视场的光强度衰减的现象.对于飞行时间成像系统，光信号的强度与测量误差之间呈反比关系，即光强度衰减使距离测量的误差增大.为了仿真渐晕现象，设传感器S的光学系统F数为NS，则由辐射LP在传感器上的辐照度ES由下式表达[14]:

(20)

其中θL的变化区间为[0,90]，在该区间上单调递减，表明像平面的辐照度ES随像方视场的增大而衰减.

2.4 传感器模型

飞行时间成像传感器可根据系统设计指标进行参数化，形成仿真参数集[13,18-19].

传感器为W×H的像元阵列，像元面积为AS.假设每个像元接收场景表面P点面元的辐射，则该像元接收的光功率PS为[4]

PS=ES·AS

(21)

其中PS是调制光源的峰值功率，调制信号占空比为50%，因此在积分时间T内像元采集的能量为[4]

(22)

该能量在像元中转换为电子空穴对，产生光生电荷的数量由量子效率νq和波长λ决定[4]:

(23)

其中ћ为普朗克常量，c为光速，q为元电荷量.像元产生的光电荷Ntotal由势阱A和势阱B收集，电荷数分别为NA和NB,有[4]

Ntotal=NA+NB

(24)

相关函数采样使用四相采样法，有[4]

(25)

(26)

其中D为解调对比度，相移φ通过下式计算[4]:

(27)

上述模型通过像元累积的总电荷量计算相移.

2.5 空域超采样模型

上述物理过程中仅将每个像元和场景表面的某个点对应，实际的物理过程中，每个像元通过光学系统与场景表面的一个区域相对应，该区域的面积与光学系统的焦距相关.当该区域的深度信息不均匀时，在像元接收信号时产生相位混叠，产生飞行像元误差.

为了对飞行像元误差进行仿真，在像元阵列的模型中引入空域超采样[10-11]，将每个像元分为WS×HS个子像元，子像元与场景表面呈理想的小孔成像模型，每个像元的深度计算值为子像元深度计算值的矢量和.

另一方面，由于像元的读出电路和控制电路包括像元的一部分面积，导致感光面积小于其物理尺寸，即存在填充因子o∈[0,1],因此每个子像元的入射光能量需要考虑填充因子对入射光能量的约束[6].

2.6 时域超采样模型

当使用四相采样法时，采样值Ii组成的4个相位图均是在同一场景中获取的，即场景为静态场景，相对位置和姿态没有变化.当场景的深度信息在积分时间内产生变化时，将产生运动伪影.在仿真系统中模拟运动伪影需要时域超采样，将积分时间T分为t个间隔，最终的相位图是每个子积分时间的相位图之和.

2.7 灰度-距离误差模型

灰度误差是反射光信号幅值的函数，取决于传感器接收的反射光强度.表现形式为当物体的反射率越低时，传感器接收到的光信号幅值越小，距离的测量结果误差越大.因此场景表面的反射性质对测量结果具有较大的影响[15-16].

传感器接收到的采样数据I可以计算得到反射信号的幅值A、相位差φd和偏置c.这3个物理量与距离测量的误差关系为[16]

(28)

(29)

(30)

因此深度测量数据的误差与信号幅值A有直接的关系，当信号幅值A越大，即反射光信号的信噪比越高，测量误差越小.当光源的调制频率确定时，相位差φd和测量距离d呈线性关系，因此测量误差和信号幅值之间的关系可简化表达为

(31)

由于信号幅值以模拟量的形式输出转换为灰度值数据，在模数转换的动态范围内亮度与灰度值呈线性关系，因此式(31)可进一步写作

(32)

2.8 光信号饱和模型

信号饱和是由于场景同时被飞行时间成像系统的光源和其他光源照明(如太阳、室内照明)，进入像元的光功率较大，导致测量数据的误差增大甚至错误的现象，其来源可分为主动光源饱和、环境光饱和两大类.出现该现象的机理是由于传感器像元以势阱收集光生电荷，可收集的电荷量存在上限值，称为满阱电荷量，数学上是光信号与电荷量之间呈非线性关系.当对光信号进行相关采集的势阱均达到满阱时，像元电路无法通过参考电压的值与势阱间的电荷量差异体现二者的相关性，导致测量结果的错误.图3为像元接收到的包含有背景光的光信号示意图，背景光信号的数学模型为不随时间变化的定制，调制光信号随时间周期性变化，二者线性叠加，当由于背景光信号产生的光生电荷在进行相关信号采集的势阱中均达到满阱值时，不存在势阱之间的信号差，此时入射光信号和参考信号没有相关性，输出的距离测量结果等于0.

图2 环境光的干扰示意图Fig.2 Schematic diagram of ambient light interference

图3为同一像元在无环境光干扰和有环境光干扰时的输出信号对比结果，其中图3(a)为无环境光照明时像元的输出信号，图3(b)为存在环境光照明时像元的输出信号.图中的Ua和Ub分别为每个像元内部两个势阱的输出电压，表征了像元采集的光生电荷量；ΔUab为像元输出的差分信号，即Ua-Ub.此信号与光源驱动信号进行相关运算，完成相位信息的解算.从图中可以看出，环境光干扰对输出信号的干扰体现为每个势阱输出信号幅值增大，即共模干扰.由于飞行时间像元的输出信号为差分形式，因此抑制了由于环境光产生的共模干扰，在势阱的动态范围内ΔUab保持不变，不影响像元的测量结果.但是环境光对飞行时间成像的影响体现为像元的动态范围的衰减，当环境光产生的光生电荷量达到势阱的满阱电荷量时

Ua=Ub=Umax

(33)

ΔUab=0

(34)

此时将ΔUab与光源驱动信号进行相关运算，解算的相移量φ等于0，因此输出的距离信息d也为0.

根据飞行时间测距原理，其测距性能与光信号的强度相关.一般使用调制对比度Km定义为[26]

(35)

其中Imax和Imin分别为传感器接收的光强度的最大值和最小值.传感器像元电路将采集的光信号与转换为模拟信号，光强信号由电信号表征.Imax和Imin转换为模拟信号电压和基准电压之差ΔUmax和ΔUmin.在无环境光时，信号的调制对比度为

(36)

当存在环境光信号时，调制对比度为

(37)

其中Pbackground为环境光在场景表面照明的光功率，Popt为飞行时间成像系统调制光源的光功率.当照明光源光功率为定值时，如果无环境光存在，则

Pbackground=0

(38)

K=KPMD

(39)

当环境光存在时，有

Pbackground>0

(40)

K

(41)

即环境光对场景表面的光强度增大导致光信号的调制对比度衰减.

图3 存在环境光干扰的像元输出信号对比图Fig.3 Pixel output signal comparison chart with ambient light interference

3 仿真算法组成

根据上述飞行时间成像的物理过程，仿真算法可分为成像系统参数模块、环境光参数模块、场景参数模块、场景数据生成模块和场景数据输出模块五个部分.各模块之间的关系如图4所示.其中，成像系统参数模块用于输入传感器参数，包括光学系统焦距、像元数、帧率、调制光源波长、调制频率等；环境光参数模块用于输入环境光参数，包括环境光波长、环境光强度、光传输过程的衰减系数等；场景参数模块用于设置场景的三维顶点数据、表面反射模型及反射率；场景数据生成模块通过渲染器生成传感器光源照明和环境光共同作用下的场景深度数据；数据输出模块将场景深度数据存储为特定格式，是仿真算法与外部程序的交互接口.

图4 仿真算法组成框图Fig.4 Block diagram of simulation algorithm composition

4 仿真参数设计

传感器参数集以TI公司的飞行时间传感器为参考[13]，主要参数如表1所示.系统光源和传感器中心位于空间中的同一位置，并将该点设置为成像系统的坐标系原点，光源的强度分布设定为理想的高斯光源,生成的深度数据均定义在相机坐标系中.

表1 飞行时间成像系统仿真参数集Tab.1 Time-of-flight imaging system simulation parameters set

对于近地轨道的小行星表面，环境光照具有以下4个特点：

(1)来自太阳的能量占环境光的绝大部分，强度为1个太阳常数即1 367 W/m2；

(2)光照形式为平行光；

(3)无大气衰减；

(4)光线矢量与场景表面法线空间角的变化范围为0°～180°.

根据上述4个特点，环境光参数的设置如表2所示，分别对环境光的强度、波长、矢量方向和衰减进行了配置.环境光强度的配置范围为0.0f至10000.0f，对应环境光强度可从0个太阳常数和1个太阳常数之间可调；环境光波长通过光谱矢量表达，表征环境光的光谱分量；光照矢量定义在成像系统坐标系中，为一个单位矢量，表征环境光传播方向；传输衰减为光能量到达场景表面时的损耗比例，设置为0.

5 数字仿真结果

5.1 仿真精度评估

在上述仿真参数集的基础上，设定场景为余弦辐射体平面，该平面到图像中心的距离分别为0.5 m、1 m、2.5 m和5 m.图5～8分别为图像中心到场景平面的距离分别为上述距离时的三维坐标值和真值.从仿真结果中可以看出，飞行时间成像的测量值在场景深度真值的基础上叠加噪声，并且由于照明光源的高斯分布和光学系统渐晕的影响，二维图像的灰度随着像方视场的增大而衰减，导致深度信息的随机误差增大.在各个测量距离上计算测量距离的均值和均方差如表3所示.

表2 环境光仿真参数集Tab.2 Ambient light parameters set

图5 0.5 m距离的余弦辐射体平面仿真结果和真值Fig.5 Cosine radiator plane simulation results and ground truth at a distance of 0.5 m

图6 1 m距离的余弦辐射体平面测量值和真值Fig.6 Cosine radiator plane simulation results and ground truth at a distance of 1 m

图7 2.5 m距离的余弦辐射体平面测量值和真值Fig.7 Cosine radiator plane simulation results and ground truth at a distance of 2.5 m

根据表3的计算结果，在典型测量距离(0.5、1、2.5、5 (m))上，飞行时间成像的距离测量的均值与理论值较为符合，测量值的均方差不大于0.055，表明在目标为余弦辐射体时，当成像灰度值相同时，飞行时间成像测量值的随机噪声与测量距离无关.

表3 典型测量距离余弦辐射平面仿真值的统计结果Tab.3 Statistical results of typical measured distance cosine radiation plane simulation results

5.2 灰度值与距离测量误差的关系

仿真的二维图像灰度值和测量精度之间的关系如图9～12所示.从仿真结果中可以看到，当二维图像中的像元灰度范围在210至255时，测量值的方差较为稳定，约为0.1 m，当灰度值低于210时，测量误差与灰度值之间呈现负相关，当灰度值为70时，测量结果的最大误差为0.25 m，测量值方差的变化趋势和理论计算一致.根据上述计算结果，该仿真系统在明确了飞行时间成像系统参数集的前提下，可对成像系统的测量能力进行评估.

图9 幅度值和测量误差仿真统计结果(0.5 m)Fig.9 Amplitude value and measurement error simulation statistics (0.5 m)

图10 幅度值和测量误差仿真统计结果(1 m)Fig.10 Amplitude value and measurement error simulation statistics (1 m)

图11 幅度值和测量误差仿真统计结果(2.5 m)Fig.11 Amplitude values and measurement errors simulation statistics (2.5 m)

图12 幅度值和测量误差仿真统计结果(5 m)Fig.12 Amplitude value and measurement error simulation statistics (5 m)

5.3 环境光干扰

包括环境光干扰的飞行时间成像仿真的参数设置如第5节中所述，环境光与成像系统光轴呈45°空间角入射到场景表面.图13～18表示了环境光对飞行时间成像的干扰情况.当没有环境光的情况下，成像系统的测量结果仅有随机噪声.当环境光功率与成像系统的调制光源功率相等时，场景中受到环境光照射的区域由于多路径效应导致测量结果的误差增大.当环境光功率增大至调制光源功率的10倍时，由于光信号饱和导致照射区域中部分像元的测量结果等于0，之后随着环境光功率的进一步增大，饱和像元数量越多，测量结果中的0值越多.该变化趋势和理论分析的结果一致.环境光功率和调制光源功率之比与饱和像元数量之间的变化曲线如图19所示.

图13 无环境光干扰的场景成像结果Fig.13 Imaging results without ambient light interference

图14 环境光功率等于1倍调制光源功率的场景成像结果Fig.14 Imaging results with ambient power equal to 1X power of the modulated source

图20和图21分别为无环境光干扰的场景和1个太阳常数照射漫反射表面的实际场景成像结果.经过对比可以看出，图20中的图像中心区域的漫反射表面没有环境光时，成像系统可对该区域生成深度数据，当漫反射表面受到1个太阳常数的环境光照射时，由于图像中心区域的光信号饱和，测量结果为0，即图21中的黑色区域，与仿真结果的变化情况相同.

图15 环境光功率等于2倍调制光源功率的场景成像结果Fig.15 Imaging results with ambient power equal to 2X power of the modulated source

图16 环境光功率等于10倍调制光源功率的场景成像结果Fig.16 Imaging results with ambient power equal to 10X power of the modulated source

图17 环境光功率等于20倍调制光源功率的场景成像结果Fig.17 Imaging results with ambient power equal to 20X power of the modulated source

图18 环境光功率等于40倍调制光源功率的场景成像结果Fig.18 Imaging results with ambient power equal to 40X power of the modulated source

图19 环境光功率比与饱和像元数关系曲线Fig.19 Relationship between ambient light power ratio and saturated pixel number

图20 无环境光照射的场景成像结果Fig.20 Scene imaging results without ambient light illumination

图21 单位太阳常数环境光照射的场景成像结果Fig.21 Scene imaging results of 1 solar constant ambient light illumination

5.4 以67P彗星为例的仿真验证试验

图22 67P/Churyumov-Gerasimenko彗星导航相机拍摄图像Fig.22 67P/Churyumov-Gerasimenko comet image

图23 灰度数据仿真结果Fig.23 Grayscale simulation results

图24 距离分布图的场景三维成像结果Fig.24 3D imaging results (top)and ground truth (bottom)of the range map

对小天体目标进行三维形貌数据遥感的形貌重构进行仿真，生成飞行时间成像系统的仿真测量数据，针对获取遥感的小天体目标为67P/Churyumov-Gerasimenko彗星，该彗星由欧空局发射的“罗塞塔”号探测器进行探测任务.图22为探测器搭载的宽视场遥感相机拍摄的灰度图像.图23为根据三维激光成像结果生成三维重构模型，使用飞行时间成像系统的仿真参数集对三维重构模型生成深度的测量结果，其高程分布结果如图24所示.从图中可以看出，飞行时间三维成像的仿真结果符合被测天体表面的高程分布，由于在仿真过程中融合了飞行时间成像系统的参数，仿真结果同时给出了地貌起伏和阴影等原因产生的测量数据误差.

6 结 论

本文针对小行星着陆应用需求，研究了飞行时间成像的光传输模型，建立了基于环境光照条件的飞行时间成像的数字仿真模型，该模型集成了完善的光源模型、光传输模型、场景表面反射模型、传感器电路模型和传感器噪声模型，可根据飞行时间传感器参数集对预设三维场景的光照条件生成飞行时间成像的数字仿真计算结果，给出包括二维灰度图像数据和三维深度的高程图像数据，并具备了完善的环境光条件影响机理模型，进行联合仿真验证试验，获得了明确的稳定的验模数据和验证结果.同时扩充了数字仿真系统的仿真对象的参数集范围，实现了对传感器的光源峰值功率、空间分辨率、积分时间、信号读出时间等主要参数的仿真优化，实现了对环境光影响相关参数的实时动态调整.该系统建立了基于小行星着陆应用数字仿真的研究验证平台，对于具有光照条件作用的飞行时间成像系统三维信息获取和以飞行时间成像系统作为数据源的三维信息算法以及方案的优化研究验证和评估具有重要的意义.