1 引 言

随着世界科技的快速发展，传统的卫星光学成像技术、星载SAR成像技术由于受到气候条件、光照条件影响及自身波段限制，分辨率和高程测量精度难以进一步提高，已不能满足全球精确测绘的需要(如绘制1∶10 000和1∶5000比例尺地图)^{[1, 2, 3, 4]}。激光技术与雷达技术结合产生的主动式传感器——激光雷达，具有不受国界限制、全天时全天候工作、高灵敏度、高分辨率、高测距精度的特点，因而受到了航空航天测绘领域的广泛研究和应用^{[5, 6, 7, 8, 9, 10]}，在地球极地冰盖探测、地表植被探测以及海洋探测中展示出了出色的高程测量能力，并且已经成功应用于月球、水星和火星的高程测量中^{[11, 12, 13]}。

星载激光雷达发散角较小(通常为几十微弧度)，但由于其轨道高度大(400～1000 km)，因此在地面形成的光斑一般为10～200 m。激光足印区域地形的起伏、地貌的分布对激光回波信号的时间、空间分布有很大影响，相当于对入射的激光脉冲进行了时间-空间调制^[14]，而回波信号的波形直接影响高程测量的精度。目前，关于激光雷达测绘性能的研究大多只针对机载试验与仿真^{[15, 16, 17]}，针对星载激光雷达三维测绘的研究还比较少，尤其是Geiger模式的单光子探测技术在航天测绘领域中的应用可行性及性能还未得到充分研究^[18]。本文首先介绍了激光三维测绘卫星高程测量的原理，并对不同地形、地貌建立物理模型，通过仿真手段得到了经过地表形貌调制后的脉冲波形的时间-空间分布，最后采用蒙特卡罗方法研究了不同地表形貌对Geiger探测模式的星载激光雷达测量精度的影响，并利用Matlab生成仿真DEM数据，对激光三维测绘卫星的测绘效果进行了仿真试验验证。

星载激光三维成像雷达探测流程如图 1所示。首先由激光发射器向探测目标地点发射激光脉冲，激光脉冲在经过大气衰减后到达地表反射区域并发生漫反射，被反射的激光回波再经过大气衰减后进入到激光探测器中。

图 1 星载激光三维成像雷达卫星测量示意图 Fig. 1 The surveying scheme of space borne 3D LiDAR

图选项

图 1中，T为大气单程透过率，n(x,y)为地面反射率分布函数，f(x,y)为地形分布函数，H为卫星轨道高度。在脉冲往返时间(几个毫秒)内，卫星会飞行一定距离(几十米)，实际激光脉冲经过的路程与垂直高度往返会有一定偏差，可以证明，卫星飞行在此段时间引起的高程测量误差很小，可以忽略。

在不考虑大气对脉冲的影响情况下，激光雷达回波信号可以看作源信号经过往返两次菲涅尔衍射过程和一次地表漫反射过程得到。激光脉冲的发射源信号一般表达式为

点击浏览公式

式中，A₀为振幅；σ_pulse为激光脉宽；z为传播距离；θ为激光发散角。垂直于传播方向应为二维(x,y)分布，且x、y方向完全等价，这里为了方便仅选取x方向的一维模式。发射脉冲的时间、空间分布及其联合分布如图 2所示。

图 2 激光脉冲的时间、空间分布特性 Fig. 2 Temporal and spatial distribution of laser pulse

图选项

可以看出，发射的激光脉冲在时间、空间维度上均为高斯分布模型。根据C.S.Gardner的反射回波理论^[19]，可以得到回波信号表达式

点击浏览公式

 式中，S为激光探测器光学口径；T为单程大气透过率；n(x)为地表反射率分布函数；A(x,z)、g(t)、f(x)分别为激光传播到z处足印振幅分布、激光脉冲时域分布和地表形状分布函数，表达式分别为

点击浏览公式

点击浏览公式

点击浏览公式

式中，E为单位脉冲能量；h(x)为地表相对起伏。

对于普通的激光测高仪，一般采用单点式、低重频测量模式，发射功率为几十毫焦，如ICESat卫星^{[11, 20]}，但对于激光三维测绘卫星，需要获得全球全覆盖的高程信息，因此要求足印较小，重复频率需达到千赫量级，且采用千元激光脉冲阵列发射模式。受限于卫星的电源输出功率，每个脉冲的能量为几百微焦量级甚至更小。理论估算可知，在1.5 m左右的探测口径条件下，回波能量约为10^-18 J，达到单光子量级，因此传统的线性探测模式已然不适用于激光三维测绘卫星，必须采用光子计数体制实现回波弱信号的探测。Geiger模式APD探测器可以达到单光子级的探测灵敏度，使得以Geiger模式探测器为基础的测距系统具有远距离测距的优势^{[21, 22]}，因此在激光三维测绘卫星中，选取Geiger模式的APD探测器实现单光子信号探测。

如图 3所示，由于回波光子数少，回波信号成为不连续的光子分布，回波信号强的区域(波峰)对应的光子分布概率大(分布密集)，回波信号弱的区域(波腹)对应的光子分布的概率小(分布稀疏)。当APD处于Geiger探测模式下时，每当探测到一个光电子便能产生雪崩饱和，输出一次探测结果。与传统线性探测模式相同，Geiger探测模式同样采用记录脉冲发射与回波时间间隔来测量距离，测距公式为

点击浏览公式

图 3 单光子探测回波信号与连续波信号关系示意图 Fig. 3 Echo pulses from single-photon-detection and normal continuous LiDAR detection

图选项

不同的是，线性探测模式采用的是质心法或阈值法来确定回波脉冲的时刻，而Geiger模式在一次探测中只能采集到一个回波信号，这个时刻便认为是回波时刻。

由于回波光子数服从Poisson分布，Geiger模式下光电子探测的概率为^[8]

点击浏览公式

式中，N为回波光电子数，当信号和噪声同时存在时，N=Nsignal+N_noise。需要指出的是，Geiger探测模式下，无法区分探测到信号的是噪声还是回波信号，存在虚警问题，可以利用相关硬件、软件方法，如距离门控制、滤光衰减法、多脉冲探测等方法降低虚警概率^{[23, 24]}，这里不作过多说明。

采用蒙特卡罗方法仿真Geiger模式探测过程思路为：将选通距离门以硬件计时精度为间隔划分为一系列独立的时间区间；依次计算每个区间被触发的概率并与生成的随机数比较，直到探测器被触发输出触发时间。算法流程如图 4所示。多次重复以上过程，以输出时间的平均值作为回波时间。

图 4 单光子探测回波时刻仿真算法流程图 Fig. 4 The flow diagram of simulation algorithm of single-photon-detection

图选项

在仿真计算过程中，选取的相关参数如表 1所示。

通常来说，即使激光足印直径为几十米，地面形状起伏、反射率也具有复杂的分布。这里笔者将足印内的地表抽象成4种数学-物理模型：斜坡地形、台阶地形、植被地貌(反射率的“斜坡”分布)和分界地貌(反射率的“台阶”分布)，如图 5所示，其他复杂地形、地貌可以表示为简单模型的组合、叠加^[14]。

图 5 4种地形地貌物理模型 Fig. 5 4 models of geography and geomorphology

图选项

斜坡地形模型中，地表起伏分布为

点击浏览公式

定义斜坡地形高度真实值为斜坡的平均高度，即斜坡高度的一半。斜坡地形下，激光脉冲回波的时间、空间分布如图 6所示。

图 6 斜坡地形回波信号的时间、空间分布特性 Fig. 6 Temporal and spatial distribution of echo pulse from slope terrain

图选项

从图 6可以看到，斜坡地形不影响回波信号的空间分布，但信号回波时间随地表高度的变化有一定延时，延时的结果是脉冲被展宽。

图 7(a)给出了回波信号波形随斜坡高度变化的关系：随着斜坡高度的增加，脉冲宽度逐渐增大，且波峰位置也在逐渐右移。从图 7(b)可以看出，随着坡度变陡，波峰位置移动，仿真得到的探测高度值也逐渐增加，且满足关系

点击浏览公式

这正是斜坡中心高度的表达式，因此探测到的h值是平均高度的无偏测量。但由于脉宽的加大，探测值的标准差也与tan θ成正比增加。

图 7 斜坡回波信号波形及测量精度与斜坡倾角的关系 Fig. 7 Relationship between echo pulse,detection precision and dip angle of the slope

图选项

台阶地形模型中，地表起伏分布为

点击浏览公式

定义台阶顶和台阶底的平均值为台阶地形的高度真实值，即h₁/2。台阶地形下，激光脉冲回波的时间、空间分布如图 8所示。回波信号的空间分布依旧不受地形影响，但时间分布曲线中间出现一个凹陷，这是由于台阶的顶和底分别给予了入射的激光脉冲不同的固定延时。

图 8 台阶地形回波信号的时间、空间分布特性 Fig. 8 Temporal and spatial distribution of echo pulse from stepped terrain

图选项

图 9给出了台阶地形对回波信号及探测精度的影响。从图 9(a)可以看出，随着台阶高度的增加，回波波形中间凹陷逐渐加大，直至分裂成两个高斯脉冲波形。图 9(b)给出了测量值、标准差与台阶高度的关系，可以看到，高度测量值与台阶高度成线性正比关系，且比例系数为1/2，是对台阶真实值的无偏测量。

图 9 台阶回波信号波形及测量精度与台阶高度的关系 Fig. 9 Relationship between echo pulse,detection precision and height of the step

图选项

形成这种现象的原因是由于回波信号的两个波形相对于中间点成对称分布，因此在大量重复情况下，探测值在平均后趋于台阶高度的一半，即趋于真实高度值。随着台阶高度的增加，波峰间距增大，测量的标准差也逐渐增大。

植被地貌情况下，由于物种的类型和稀疏程度不同，反射率呈现特定的分布，这里抽象为反射率的“斜坡”模型

点击浏览公式

反射率“斜坡”地貌下，激光脉冲回波的时间、空间分布如图 10所示。

点击浏览公式

图 10 植被地貌回波信号的时间、空间分布特性 Fig. 10 Temporal and spatial distribution of echo pulse from vegetation terrain

图选项

从图 10可以看出，由于地表反射率呈“斜坡”分布，回波信号中心向反射率大的方向偏移，但不影响信号的时间分布，因此高程测量误差来源于回波光子地面分布点的偏移。假设地面探测点偏移量为100 m，卫星轨道高度为400 km，则造成的高程测量误差为

误差在可以忽略范围内。实际上，地面反射率的不均匀分布造成激光光斑中心的偏移量很小，且在足印范围内，而激光三维测绘卫星的足印不大于10 m，因此植被地貌造成的高程测量误差可以忽略不计。

在水域与陆地、植被与沙地等分界面上，反射率会出现“台阶”式的跳跃，用以下模型表示

点击浏览公式

反射率“台阶”地貌下，激光脉冲回波的时间、空间分布如图 11所示。

图 11 分界地貌回波信号的时间、空间分布特性 Fig. 11 Temporal and spatial distribution of echo pulse from dividing terrain

图选项

从图 11可以看出，分界地貌不影响回波信号的时间分布，但由于反射率的“台阶”式分布，回波信号在分界线的两边强度会有不均等的衰减，结果同样是回波信号强度重心的空间偏移。根据3.3小节的分析可知，该偏移造成的高程测量误差很小，可以忽略不计。

采用Matlab建立模拟的DEM数据模型，如图 12(a)所示。仿真DEM数据面积为500 m×500 m，网格大小为1 m×1 m，包括斜坡、台阶地形和分界、植被地貌，其中，斜坡、台阶处地表反射率取0.3，植被地貌反射率呈正弦函数变化，在DEM图上以正弦起伏表示。

图 12 激光三维测绘仿真试验 Fig. 12 Simulation experiment of 3D LiDAR surveying

图选项

激光三维测绘取网格分辨率为4 m×4 m。图 12中(b)、(c)、(d)分别为采用激光雷达探测1次、3次、8次时的结果，对应的测量误差平均值分别为：0.5 m、0.3 m、0.05 m。可以看出，在斜坡、台阶地形处，激光雷达的探测效果与原DEM数据很好地吻合，数据有一定随机误差，随着探测次数的累积误差逐渐减小；在分界、植被地貌处，探测结果也与原DEM数据具有良好的一致性。

本文建立了不同地形、表貌物理模型，分析了不同地表空间起伏、反射率分布情况下的激光回波特性，并采用蒙特卡罗方法研究了不同地表形貌对Geiger探测模式的激光三维测绘卫星测量精度的影响，通过仿真试验验证，得到如下结论：

(1) 不同地形对回波脉冲的影响主要体现在脉冲的时间分布上，表现为回波信号的展宽、分裂等，但大量重复得到的测量平均值是平均高度值的无偏测量；随着地表坡度的增大、台阶高度的增大，测量标准差均增大。

(2) 不同地貌对回波信号的影响主要体现在脉冲的空间分布上，表现为回波信号强度中心的偏移，但偏移量较小，且造成的测量误差很小，可以忽略不计。

研究结果表明：激光三维测绘卫星对垂直陡变地形(如城市建筑)、斜坡地形(如山坡)、分界地貌(如水陆分界)和折射率起伏地貌(如植被地表等)具有良好的探测能力。