CN111007529A - 一种全链路的光子计数激光测高仪点云的生成方法 - Google Patents

Abstract

本发明公开了一种全链路的光子计数激光测高仪点云的生成方法，建立了包含星载激光测高仪测高过程中的地面光斑坐标解算、大气折射偏移修正、目标几何和辐射特性、探测器响应等的全链路仿真模型，结合卫星平台参数，激光测高仪系统参数，环境参数，目标特性参数，仿真生成高精度的单光子点云。该激光测高仪点云的生成方法可以准确地调整各项系统参数，环境参数，具有快速，易于模拟各种情况，且扩展性好的优点。可以以ICESat‑2系统参数为输入，仿真探测随机生成的地形表面，生成了高质量光子点云。

Description

一种全链路的光子计数激光测高仪点云的生成方法

技术领域

本发明涉及数据仿真技术领域，具体涉及一种全链路的光子计数激光测高仪点云的生成方法。

背景技术

相较于线性体制的激光测高仪，光子计数激光测高仪使用单光子探测器作为接收器件，灵敏度提高了2到3个量级，很好的解决了单脉冲能量和重频之间的矛盾，带来了小体积、多波束、高重频等诸多优势，是未来激光测高仪的发展趋势。

本申请发明人在实施本发明的过程中，发现现有技术的方法，至少存在如下技术问题：

目前，国内光子计数体制的激光测高仪仿真模型，尤其是全链路仿真还鲜有报道。由于探测器工作在光子计数模式，不同于传统线性体制探测器，其输出是服从特定概率分布的离散的光子事件，因而最终的数据产品是由大量离散事件组成的光子点云。如果采用实际实验获取所需的不同条件的点云数据，则需要对多种条件下进行实验，不仅系统参数调整困难而且成本极高。

由此可知，现有技术中存在由于参数调整困难导致点云数据生成准确和效率不高的技术问题。

发明内容

有鉴于此，本发明提供了一种全链路的光子计数激光测高仪点云的生成方法，用以解决或者至少部分解决现有技术中存在由于参数调整困难导致的点云数据生成准确和效率不高的技术问题。

为了解决上述技术问题，本发明提供了一种全链路的光子计数激光测高仪点云的生成方法，包括：

步骤S1：建立激光发射坐标系、卫星本体坐标系和地面坐标系；

步骤S2：确定激光方向矢量、指向偏差和卫星在地面相对坐标系中的偏转，根据激光方向矢量、指向偏差和卫星在地面相对坐标系中的偏转，通过坐标转换关系得到地面坐标系下激光脉冲矢量；

步骤S3：根据激光脉冲矢量、激光发射点在相对地面坐标系下的坐标以及大地坐标系，计算出地面坐标系下所有光斑脚点的坐标和激光入射指向角；

步骤S4：获取光斑位置处的高程轮廓，对地面坐标系下所有光斑脚点的坐标和激光入射指向角进行修正，并根据修正后的光斑脚点的坐标、修正后的激光入射指向角以及光斑位置处的高程轮廓，计算光斑位置法向指向角；

步骤S5：获得光子计数激光测高仪的发射系统参数、激光雷达与目标之间的单程大气透过率以及激光脉冲在空域分布的初始分布，计算单脉冲到达地面光斑位置的能量分布；

步骤S6：获得地面目标反射率和激光雷达接收系统参数，根据单脉冲到达地面光斑位置的能量分布、目标反射率、激光雷达接收系统参数、光斑位置法向指向角与修正后的激光入射指向角以及目标高程轮廓，计算目标的单位冲击响应，再根据目标的单位冲击响应计算出光斑脚点的回波信号波形，目标高程轮廓为光斑位置处的高程轮廓；

步骤S7：根据激光雷达方程计算单位时间内探测器接收的由目标发射的背景光噪声；

步骤S8：根据光斑脚点的回波信号波形以及探测器接收的由目标发射的背景光噪，利用单光子探测器进行响应探测，生成全链路的光子计数激光测高仪点云。

在一种实施方式中，步骤S1具体包括：

以激光发射的参考起点为原点，激光出射方向为Z轴正方向，Y轴垂直于卫星轨道平面，右手准则确定X轴建立激光发射坐标系O_L-X_LY_LZ_L；以卫星质心为原点，X轴指向卫星飞行方向，Z轴指向天底方向，Y轴垂直于卫星轨道平面建立卫星本体坐标系O_B-X_BY_BZ_B；以WGS84坐标系为基础建立大地坐标系(O-XYZ)。

在一种实施方式中，步骤S2具体包括：

步骤S2.1：将激光方向矢量设置为绕卫星本体坐标系X_B轴逆时针旋转θ_lx、绕Y_B轴逆时针旋转θ_ly的单位矢量，指向偏差为Δθ_lx(t)和Δθ_ly(t)，根据公式(1)计算卫星本体坐标系激光脉冲方向矢量P_L(t)，

式(1)中，Δθ_lx(t)和Δθ_ly(t)是由各种因素引起的指向偏差，是随时间变化的量；

步骤S2.2：获取卫星在地面相对坐标系中的偏转，偏转包括偏航角ω、俯仰角

和侧滚角κ，根据公式(2)计算得到地面坐标系下激光脉冲矢量P(t)，

其中，ω、

κ分别为卫星定姿系统测定的卫星在地面相对坐标系中的偏航角、俯仰角和侧滚角，M(ω),

和M(κ)为坐标旋转矩阵。

在一种实施方式中，步骤S3具体包括：

步骤S3.1：获得卫星平台参数：卫星相对地面轨道高度H₀、卫星运行速度V、卫星轨道路径T_s，；

步骤S3.2：获得激光发射重频设定R_L，激光发射坐标系在卫星本体坐标系下的坐标；

步骤S3.3：根据卫星平台参数、激光发射重频设定R_L，激光发射坐标系在卫星本体坐标系下的坐标，计算出地面坐标系激光发射坐标[X_L0(n)，Y_L0(n)，Z_L0(n)]；

步骤S3.4：根据地面坐标系激光发射坐标以及激光脉冲矢量P(t)，计算出地面坐标系的所有光斑脚点坐标[X₀(n),Y₀(n)]以及激光入射指向角θ₀(x,y)。

在一种实施方式中，步骤S4具体包括：

步骤S4.1：获取大气参数：不同大气压位置的大气温度T_a、湿度H_a、位势高度G_a；

步骤S4.2：根据大气参数，构建大气折射偏移修正模型，并利用大气折射偏移修正模型对地面坐标系下所有光斑脚点的坐标和激光入射指向角进行修正，得到相对地面坐标系下的修正后的光斑坐标[X(n),Y(n)]及入修正后的射指向角θ_L(n)，其中大气折射偏移修正模型中大气折射率如公式(3)：

其中，k₁(λ)和k₂(λ)是与波长λ有关的函数，T为温度，P_d和P_w、Z_d和Z_w分别为干空气和水汽的分压强与压缩率；

步骤S4.3：获取光斑位置处的高程轮廓ξ(x,y)，根据修正后的光斑脚点的坐标、修正后的激光入射指向角以及光斑位置处的高程轮廓，计算光斑位置法向指向角θ_n(x,y)。

在一种实施方式中，步骤S5具体包括：

步骤S5.1：获取光子计数激光测高仪发射系统参数：激光单脉冲发射能量E_t、激光脉冲均方根脉宽σ_s和发射光学系统效率η_t；

步骤S5.2：获得激光雷达与目标之间的单程大气透过率T_v、激光脉冲在空域分布的初始分布，根据公式(4)计算单脉冲到达地面光斑位置[X(n),Y(n)]的能量分布E_i(x,y)：

E_i(x,y)＝η_tT_vE₀(x,y) (4)

其中，激光脉冲在空域分布为高斯分布，E₀(x,y)为激光脉冲初始分布，η_r为光学系统的传输效率，T_v为星载激光测高仪与目标之间的大气透过率。

在一种实施方式中，步骤S6具体包括：

步骤S6.1：获得地面目标反射率ρ，激光雷达接收系统参数：探测器的量子效率η_q、接收光学系统效率η_r、接收光学望远镜的有效口径A_r；

步骤S6.2：根据公式(5)计算出信号回波信号光子数分布Ns(x,y)，

其中，η_r为光学系统的传输效率，hv为单个光子能量，A_r为接收光学系统的有效口径，T_v为星载激光测高仪与目标之间的大气透过率，ρ为随机地形表面反射率，H₀是卫星轨道高度，θ_g为激光入射指向角与光斑位置法向指向角之间的夹角；

步骤S6.3：根据信号回波信号光子数分布和目标高程信息计算出的激光往返传播时间Δt(x,y)，累加同一时刻的回波信号，得到针对该目标的单位冲击响应W_s0(t)；

步骤S6.4：获得激光单脉冲的时域分布为S(t)，将目标的单位冲击响应W_s0(t)与激光单脉冲的时域分布为S(t)进行卷积，得到脚点的回波信号波形W_s(t)。

在一种实施方式中，步骤S7具体包括：

获得太阳高度角θ_h和方位角θ_a，结合随机地形的坡度角、接收系统接收光学系统的滤波片带宽△λ、探测器的量子效率η_q以及大气层外太阳光谱辐照度N_λ ⁰，根据激光雷达方程，计算背景光噪声N_b，计算公式为(6)：

其中，N_λ ⁰为与波长相关的大气层外太阳光谱辐照度，△λ表示接收光学系统的滤波片带宽，R为地球半径，hv为单个光子能量，ρ为随机地形表面反射率，θ_r为接收单元的视场角，A_r为接收光学系统的有效口径，η_r为光学系统的传输效率，Ψ为太阳照射方向与目标表面法向方向的夹角，T_v是激光雷达与目标之间的双程大气透过率。

在一种实施方式中，步骤S8具体包括：

步骤S8.1：获得单光子探测器的死区时间t_d、时间分辨率τ；i

步骤S8.2：结合有效信号回波波形W_s(t)和背景光噪声N_b，计算出每个时间片的平均信号光子数；

步骤S8.3：根据探测器响应模型计算出每个时间片的探测概率P_sn(t,t+τ)，其中，探测器响应模型为公式(7)：

P_sn(t,t+τ,λ>0)＝1-exp(-K_sn(t,t+τ)) (7)

其中，K_sn(t,t+τ)为时间片τ内的平均信号光子数；

步骤S8.4：使用蒙特卡洛仿真方法模拟探测器的响应过程及死区效应，产生0-1均匀分布随机数，当小于探测概率P_sn(t,t+τ)，则生成光子事件，从前至后，每响应一个光子事件，后续t_d时间内不响应光子事件，生成最终光子点云。

本申请实施例中的上述一个或多个技术方案，至少具有如下一种或多种技术效果：

本发明提供的一种全链路的光子计数激光测高仪点云的生成方法，首先建立激光发射坐标系、卫星本体坐标系和地面坐标系；然后进行地面光斑坐标解算得到地面坐标系下所有光斑脚点的坐标和激光入射指向角；接着利用大气折射延迟修正模型对地面坐标系下所有光斑脚点的坐标和激光入射指向角进行修正，结合光斑位置处的高程轮廓，计算光斑位置法向指向角；接着计算单脉冲到达地面光斑位置的能量分布；接下来计算出光斑脚点的回波信号波形，然后根据激光雷达方程计算单位时间内探测器接收的由目标发射的背景光噪声；最后根据光斑脚点的回波信号波形以及探测器接收的由目标发射的背景光噪，利用单光子探测器进行响应探测，生成全链路的光子计数激光测高仪点云。

本发明的方法包含星载激光测高仪测高过程中的地面光斑坐标解算、大气折射偏移修正、目标几何和辐射特性、探测器响应等的全链路仿真模型，并结合卫星平台参数，激光测高仪系统参数，环境参数，目标特性参数，仿真生成高精度的单光子点云。该方法可以准确地调整各项系统参数，环境参数，具有快速，易于模拟各种情况，且扩展性好的优点，提高了点云的生成效率和准确性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例中一种全链路的光子计数激光测高仪点云的生成方法的流程图；

图2是星载激光测高仪运行时各坐标关系示意图。

图3是星载激光测高仪光束偏转示意图，其中，H₀为星载激光测高仪轨道高度以，

为天顶角，R_e表示地球半径，X₀表示激光脚点理论位置，X'表示激光脚点的实际位置，则Δs表示激光脚点的横向偏移量；θ为光束与地表水平面之间的夹角,θ_b为激光脚点-激光测高仪连线与地面之间的夹角，则Δθ为由于大气折射引起的角度偏离。

图4是仿真1km×1km随机地形表面等高图，卫星平台发射激光数量为3波束。其中圆点为不考虑指向偏差在天底方向下接的光斑位置，其中黑色实心圆为不考虑大气折射引起角度偏离时的光斑位置，X型为考虑了各项输入参数后的实际光斑位置。

图5是测高仪非天底方向入射地表斜面示意图，其中

为激光与天底方向的夹角，ξ(x,y)为目标地表表面高程轮廓，S_∥和S_⊥分别表示平行和垂直系统飞行方向的目标地表斜坡角，Δξ表示地表粗糙度引起的高程变化。

图6是仿真模型以ICESat-2系统参数为输入，仿真探测随机生成的地形表面，最后产生的点云图。其中，虚线为沿轨方向地形表面的实际高程轮廓。

具体实施方式

本发明提出一种全链路的光子计数激光测高仪点云的生成方法，可以对激光测高过程中激光脉冲发射，经大气传输、被随机地形反射最后再经过大气传输，最后被探测器接收并经过量化采集得到点云的全过程进行仿真。

为了达到上述目的，本发明的主要构思如下：

建立了包含星载激光测高仪测高过程中的地面光斑坐标解算、大气折射偏移修正、目标几何和辐射特性、探测器响应等的全链路仿真模型，结合卫星平台参数，激光测高仪系统参数，环境参数，目标特性参数，仿真生成高精度的单光子点云。该仿真方法可以准确地调整各项系统参数，环境参数，具有快速，易于模拟各种情况，且扩展性好的优点。仿真方法以ICESat-2系统参数为输入，仿真探测随机生成的地形表面，生成了高质量光子点云。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例一

本实施例提供了一种全链路的光子计数激光测高仪点云的生成方法，请参见图1，该方法包括：

步骤S1：建立激光发射坐标系、卫星本体坐标系和地面坐标系；

步骤S2：确定激光方向矢量、指向偏差和卫星在地面相对坐标系中的偏转，根据激光方向矢量、指向偏差和卫星在地面相对坐标系中的偏转，通过坐标转换关系得到地面坐标系下激光脉冲矢量；

步骤S3：根据激光脉冲矢量、激光发射点在相对地面坐标系下的坐标以及大地坐标系，计算出地面坐标系下所有光斑脚点的坐标和激光入射指向角。

步骤S4：获取光斑位置处的高程轮廓，对地面坐标系下所有光斑脚点的坐标和激光入射指向角进行修正，并根据修正后的光斑脚点的坐标、修正后的激光入射指向角以及光斑位置处的高程轮廓，计算光斑位置法向指向角。

步骤S5：获得光子计数激光测高仪的发射系统参数、激光雷达与目标之间的单程大气透过率以及激光脉冲在空域分布的初始分布，计算单脉冲到达地面光斑位置的能量分布；

步骤S6：获得地面目标反射率和激光雷达接收系统参数，根据单脉冲到达地面光斑位置的能量分布、目标反射率、激光雷达接收系统参数、光斑位置法向指向角与修正后的激光入射指向角以及目标高程轮廓，计算目标的单位冲击响应，再根据目标的单位冲击响应计算出光斑脚点的回波信号波形，目标高程轮廓为光斑位置处的高程轮廓；

步骤S7：根据激光雷达方程计算单位时间内探测器接收的由目标发射的背景光噪声；

步骤S8：根据光斑脚点的回波信号波形以及探测器接收的由目标发射的背景光噪，利用单光子探测器进行响应探测，生成全链路的光子计数激光测高仪点云。

在一种实施方式中，步骤S1具体包括：

以激光发射的参考起点为原点，激光出射方向为Z轴正方向，Y轴垂直于卫星轨道平面，右手准则确定X轴建立激光发射坐标系O_L-X_LY_LZ_L；以卫星质心为原点，X轴指向卫星飞行方向，Z轴指向天底方向，Y轴垂直于卫星轨道平面建立卫星本体坐标系O_B-X_BY_BZ_B；以WGS84坐标系为基础建立大地坐标系(O-XYZ)。

其中，激光测高仪运行时各个坐标系时间的转换关系如图2所示，在仿真中，考虑到运算时间和效率，通常只能对数km范围内的光斑脚点进行仿真，卫星与地面的关系可以进行一定简化，模型假设姿态测量坐标系、轨道坐标系和卫星本体坐标系相互重合。其中(O_B-X_BY_BZ_B)为卫星本体坐标系，(O_L-X_LY_LZ_L)为激光测量坐标系，(O_G-X_GY_GZ_G)为GPS天线测量坐标系，(O-XYZ)为地面坐标系，圆圈为地面光斑。

在一种实施方式中，步骤S2具体包括：

步骤S2.1：将激光方向矢量设置为绕卫星本体坐标系X_B轴逆时针旋转θ_lx、绕Y_B轴逆时针旋转θ_ly的单位矢量，指向偏差为Δθ_lx(t)和Δθ_ly(t)，根据公式(1)计算卫星本体坐标系激光脉冲方向矢量P_L(t)，

式(1)中，Δθ_lx(t)和Δθ_ly(t)是由各种因素引起的指向偏差，是随时间变化的量；

步骤S2.2：获取卫星在地面相对坐标系中的偏转，偏转包括偏航角ω、俯仰角

和侧滚角κ，根据公式(2)计算得到地面坐标系下激光脉冲矢量P(t)，

其中，ω、

κ分别为卫星定姿系统测定的卫星在地面相对坐标系中的偏航角、俯仰角和侧滚角，M(ω),

和M(κ)为坐标旋转矩阵。

在一种实施方式中，步骤S3具体包括：

步骤S3.1：获得卫星平台参数：卫星相对地面轨道高度H₀、卫星运行速度V、卫星轨道路径T_s，；

步骤S3.2：获得激光发射重频设定R_L，激光发射坐标系在卫星本体坐标系下的坐标；

步骤S3.3：根据卫星平台参数、激光发射重频设定R_L，激光发射坐标系在卫星本体坐标系下的坐标，计算出地面坐标系激光发射坐标[X_L0(n)，Y_L0(n)，Z_L0(n)]；

步骤S3.4：根据地面坐标系激光发射坐标以及激光脉冲矢量P(t)，计算出地面坐标系的所有光斑脚点坐标[X₀(n),Y₀(n)]以及激光入射指向角θ₀(x,y)。

具体来说，步骤S2～步骤S3可以通过地面光斑坐标解算模型来实现，其中地面光斑坐标解算模型的相关说明如下：

假设激光发射坐标系在卫星本体坐标系下的坐标为[X_L,Y_L,Z_L]，激光方向矢量是绕卫星本体坐标系X_B轴逆时针旋转θ_lx，绕Y_B轴逆时针旋转θ_ly的单位矢量。则在卫星本体坐标系下，激光脉冲的方向矢量可以表示为：

当考虑平台受热等因素引起的指向随机误差时，激光脉冲的方向矢量可以表示为公式(1)的形式，式中，Δθ_xl(t)和Δθ_xl(t)是由各种因素引起的指向偏差，是随时间变化的量。由于激光脉冲的方向矢量是单位矢量，因此坐标系的平移对其没有影响，因此在地面相对坐标系下激光脉冲的方向矢量表示公式(2)的形式。

其中，公式(2)中坐标旋转矩阵M(ω),

和M(κ)可以分别表示为

此时，在不考虑大气折射引起光斑偏移的情况下，认为激光脉冲沿发射的方向直线传播，依据激光脉冲方向矢量P＝[m,n,p]^T，结合激光发射点在相对地面坐标系下的坐标(x_l0,y_l0,z_l0)，可以给出激光脉冲路径的直线参数方程

与地面相对坐标系的XOY平面，即z＝0联立，即可解出激光光斑在地面相对标系下的坐标为(x_l0-z₁₀n/p,y_l0-z₁₀m/p,0)。

在一种实施方式中，步骤S4具体包括：

步骤S4.1：获取大气参数：不同大气压位置的大气温度T_a、湿度H_a、位势高度G_a；

步骤S4.2：根据大气参数，构建大气折射偏移修正模型，并利用大气折射偏移修正模型对地面坐标系下所有光斑脚点的坐标和激光入射指向角进行修正，得到相对地面坐标系下的修正后的光斑坐标[X(n),Y(n)]及入修正后的射指向角θ_L(n)，其中大气折射偏移修正模型大气折射率如公式(3)：

其中，k₁(λ)和k₂(λ)是与波长λ有关的函数，T为温度，P_d和P_w、Z_d和Z_w分别为干空气和水汽的分压强与压缩率；

步骤S4.3：获取光斑位置处的高程轮廓ξ(x,y)，根据修正后的光斑脚点的坐标、修正后的激光入射指向角以及光斑位置处的高程轮廓，计算光斑位置法向指向角θ_n(x,y)。

具体来说，步骤S4可以通过构建大气折射延迟修正模型来实现，具体如下：

大气折射率受湿度、温度、压强等因素影响，在激光传输路径中往往呈现非均匀分布。星载激光测高仪发出的光束在进入大气层后，由于大气折射率不均匀，将发生折射现象，激光光程增加，激光脚点发生水平偏移，这为激光脚点精确定位引入系统误差。大气延迟修正模型计算激光传输路径中每个位置的折射率，通过光线追迹的方法解算系统误差，精确定位激光脚点。

对于测高仪的激光光波段，折射率n可以用公式(3)表示，公式(3)中，k₁(λ)和k₂(λ)是与波长λ有关的函数，T为温度，P_d和P_w、Z_d和Z_w分别为干空气和水汽的分压强与压缩率。由于地球引力作用，温度，压强与湿度主要体现在垂直方向上的不均匀，水平方向上的相对均匀。在进行光线追迹时，利用26个不同位势高度的温度、压强与湿度计算出该位置的折射率，结合折射率在不同位势高度满足的统计规律，对30km的大气层每10m进行一次分层，由三次样条插值估计该层的折射率，认为同一层内的折射率保持不变。

如图3所示，不考虑大气折射时的光路，此时光斑位于X₀处；而根据实际情况，光线在经过大气时，会不断发生折射，而最终处于地表的X′处，故X₀与X′之间的距离Δs，就是大气带来的横向漂移。θ为光束与地表水平面之间的夹角,θ_b为激光脚点-激光测高仪连线与地面之间的夹角，则由于大气折射引起的角度偏离满足公式Δθ＝θ_b-θ。每一层的折射率已经由上述插值方法获得，依据折射定律对其进行光线追迹，计算得到地面的横向漂移修正量Δs和角度偏离修正项Δθ，对地面光斑坐标及入射角进行修正。

请参见图4和图5，为仿真1km×1km随机地形表面等高图以及测高仪非天底方向入射地表斜面示意图。

在一种实施方式中，步骤S5具体包括：

步骤S5.1：获取光子计数激光测高仪发射系统参数：激光单脉冲发射能量E_t、激光脉冲均方根脉宽σ_s和发射光学系统效率η_t；

步骤S5.2：获得激光雷达与目标之间的单程大气透过率T_v、激光脉冲在空域分布的初始分布，根据公式(4)计算单脉冲到达地面光斑位置[X(n),Y(n)]的能量分布E_i(x,y)：

E_i(x,y)＝η_tT_vE₀(x,y) (4)

其中，激光脉冲在空域分布为高斯分布，E₀(x,y)为激光脉冲初始分布，η_r为光学系统的传输效率，T_v为星载激光测高仪与目标之间的大气透过率。

在一种实施方式中，步骤S6具体包括：

步骤S6.1：获得地面目标反射率ρ，激光雷达接收系统参数：探测器的量子效率η_q、接收光学系统效率η_r、接收光学望远镜的有效口径A_r；

步骤S6.2：根据公式(5)计算出信号回波信号光子数分布Ns(x,y)，

其中，η_r为光学系统的传输效率，hv为单个光子能量，A_r为接收光学系统的有效口径，T_v为星载激光测高仪与目标之间的大气透过率，ρ为随机地形表面反射率，H₀是卫星轨道高度，θ_g为激光入射指向角与光斑位置法向指向角之间的夹角；

步骤S6.3：根据信号回波信号光子数分布和目标高程信息计算出的激光往返传播时间Δt(x,y)，累加同一时刻的回波信号，得到针对该目标的单位冲击响应W_s0(t)；

步骤S6.4：获得激光单脉冲的时域分布为S(t)，将目标的单位冲击响应W_s0(t)与激光单脉冲的时域分布为S(t)进行卷积，得到脚点的回波信号波形W_s(t)。

具体来说，步骤S5～步骤S6可以通过构建有效信号回波模型来实现，相关理论与实现过程如下：

光子计数激光测高仪激光发射脉冲的初始分布一般可用高斯函数来描述。依据Gardner理论，在只考虑大气衰减的条件下，把随机地形表面当作漫反射目标物，可以把激光脉冲到达地球表面目标物的过程看作夫琅禾费衍射，而高斯分布的激光脉冲经过夫琅禾费衍射仍然服从高斯分布。所以可以用高斯函数描述光斑到达随机地形时的归一化的空域和时域分布：

其中，T₀为回波信号的时间重心，x₀，y₀分别为光斑中心横纵坐标，σ_t，σ_s为激光脉冲时域空域的脉宽。一般将随机地形表面近似为朗伯反射体，则由激光雷达方程(公式5)可以计算单光子探测器件接收到的单脉冲平均回波信号光子数N_s(x,y)，其中E_i(x,y)是激光到达地面的单脉冲能量，可以表示为公式(4)，

其中，E₀(x,y)为激光脉冲初始分布，η_r是光学系统的传输效率，hv是单个光子能量，A_r是接收光学系统的有效口径，T_v是星载激光测高仪与目标之间的大气透过率，ρ是随机地形表面反射率，H₀是卫星轨道高度，θ_g是激光测高仪的指向角与目标表面法向方向的夹角，其中目标法向方向可用目标高程轮廓ξ(x,y)求出。

不同面元的高程不同，回波信号光子的返回的时间也不同，其返回时间T₀可以近似表示为T₀(x,y)≈2(H₀-(x,y))/c。激光雷达探测器只分辨回波信号的时间信息，所以累加平均信号回波光子数N_s(x,y)在同一时间的能量，得到回波信号波形W_s0(t)。而脉冲时域分布为式(9)所描述的高斯分布，所以最终回波波形W_s(t)由W_s0(t)与脉冲时域分布S(t)的卷积得到，卷积公式为：

在一种实施方式中，步骤S7具体包括：

获得太阳高度角θ_h和方位角θ_a，结合随机地形的坡度角、接收系统接收光学系统的滤波片带宽△λ、探测器的量子效率η_q以及大气层外太阳光谱辐照度N_λ ⁰，根据激光雷达方程，计算背景光噪声N_b，计算公式为(6)：

其中，N_λ ⁰为与波长相关的大气层外太阳光谱辐照度，△λ表示接收光学系统的滤波片带宽，R为地球半径，hv为单个光子能量，ρ为随机地形表面反射率，θ_r为接收单元的视场角，A_r为接收光学系统的有效口径，η_r为光学系统的传输效率，Ψ为太阳照射方向与目标表面法向方向的夹角，T_v是激光雷达与目标之间的双程大气透过率。

具体来说，步骤S7可以通过构建噪声模型来实现，相关理论与实现过程如下：

光子计数激光测高仪在测距过程中的噪声主要为背景光噪声，同样将随机地形表面近似为朗伯反射体，由激光雷达方程可知单位时间内探测器接收到由目标反射的背景光噪声Nb由公式(6)确定。

式中，N_λ ⁰是与波长相关的大气层外太阳光谱辐照度，对于532nm波长，日照条件下N_λ ⁰的值约为1.83w/m²·nm。△λ表示接收光学系统的滤波片带宽，θ_r是接收单元的视场角(全角)，A_r为接收光学系统的有效口径，η_r是光学系统的传输效率，Ψ是太阳照射方向与目标表面法向方向的夹角，T_v是激光雷达与目标之间的双程大气透过率,用该式计算出。仿真过程中输入特定条件下ICESat-2搭载的ATLAS工作环境和系统参数(目标坡度为0度、正入射、望远镜直径为0.8m、大气单程透过率为0.8、工作波长为532nm、日照条件下N_λ ⁰的值约为1.83w/m²·nm、接收光学系统的滤波片带宽为38pm、光学传输效率为0.504、探测器的量子效率为0.15，接收视场角为85μrad)，在随机地形中计算的噪声率为MHz量级。

在一种实施方式中，步骤S8具体包括：

步骤S8.1：获得单光子探测器的死区时间t_d、时间分辨率τ；

步骤S8.2：结合有效信号回波波形W_s(t)和背景光噪声N_b，计算出每个时间片的平均信号光子数；

步骤S8.3：根据探测器响应模型计算出每个时间片的探测概率P_sn(t,t+τ)，其中，探测器响应模型为公式(7)：

P_sn(t,t+τ,λ>0)＝1-exp(-K_sn(t,t+τ)) (7)

其中，K_sn(t,t+τ)为时间片τ内的平均信号光子数；

步骤S8.4：使用蒙特卡洛仿真方法模拟探测器的响应过程及死区效应，产生0-1均匀分布随机数，当小于探测概率P_sn(t,t+τ)，则生成光子事件，从前至后，每响应一个光子事件，后续t_d时间内不响应光子事件，生成最终光子点云。

具体来说，步骤S8可以通过构建探测器响应模型来实现，相关理论与实现过程如下：

由于单光子探测器响应信号的有无，不能反映信号的大小，探测的光子事件呈现概率分布。所以用概率统计的方法，分析探测器每一微小时间段τ内发生光子事件的概率(该值取探测器响应的最小时间分辨率)。对于目前星载激光测高仪使用的主流单光子探测器GM-APD和PMT，其(t,t+τ)时间段产生的总的平均光电子数K_sn(t,t+τ)可以描述为平均信号光电子数K_s(t,t+τ)、平均背景光噪声光电子数K_b(t,t+τ)和平均暗噪声光电子数N_d(t,t+τ)之，计算公式为公式(15)

其中，N_s(t,t+τ)、N_b(t,t+τ)分别表示(t,t+τ)时间段平均回波信号光子数与平均背景光噪声光子数。依据噪声模型，背景光噪声率一般可达MHz量级。而随着近年探测器技术的不断发展，暗噪声率已经降至KHz量级或更低，远小于背景光噪声，在探测器响应时可以忽略不记，所以式(15)可以简化为公式(16)：

K_sn(t,t+τ)≈η_qN_s(t,t+τ)+τη_qN_b(t,t+τ) (16)

结合探测器探测理论与统计光学理论，可知单光子探测器在(t,t+τ)时间段产生光电子的概率P_sn(t,t+τ)为公式(7).

设定探测门范围为R_g，则探测过程中共分为N_g＝R_g/τ个时间区间，探测器死区时间为t_d时，则探测器每探测一个光子，后续的N_d＝t_d/τ个时间区间无法响应光子事件。根据公式(7)可以计算出门范围内的每一个区间的光子探测概率P_sn(t,t+τ)，本实施方式中利用蒙特卡洛仿真的方法来模拟探测器的响应及死区过程，首先生成长度为N_g的0-1均匀分布随机数序列，然后将每个时间片的探测概率与对应的随机数进行比较，小于探测概率则认为被探测到，模拟光电子的随机产生过程，最后由前至后响应光子事件，每响应一个光子事件，后面的N_d个区间均不再响应光子事件，最终生成单光子点云。请参见图6，为仿真模型以ICESat-2系统参数为输入，仿真探测随机生成的地形表面，最后产生的点云图。其中虚线为沿轨方向地形表面的实际高程轮廓。

综上所述，本发明主要有以下优点：

1)依据激光测高仪系统参数，环境参数和目标高程轮廓与反射系数，结合卫星平台参数，可以快速仿真高精度单光子点云。且各项参数调节十分方便，易于实现，仿真结果与单光子探测器实际结果非常接近，可以用于后续数据处理与精度分析。

2)根据地面光斑位置解算模型、大气折射延迟修正模型、有效信号回波模型、噪声模型和探测器响应模型，描述了光子计数激光测高仪的完整探测过程。该点云仿真方法每一模块耦合度较低，当更换激光发射系统，传输环境，目标物或者探测器时，只需更换相应的仿真模块，有着良好的可扩展性。

下面通过具体示例对本发明提供的方法进行说明。

对于一个参数已知的单光子激光雷达系统，基于随机地形的单光子激光雷达点云生成方法的步骤如下：

a.如图2所示建立坐标系：以激光发射的参考起点为原点，激光出射方向为Z轴正方向，Y轴垂直于卫星轨道平面，右手准则确定X轴建立激光发射坐标系(O_L-X_LY_LZ_L)；以卫星质心为原点，X轴指向卫星飞行方向，Z轴指向天底方向，Y轴垂直于卫星轨道平面建立卫星本体坐标系(O_B-X_BY_BZ_B)；以WGS84坐标系为基础建立大地坐标系(O-XYZ)。

b.将激光方向矢量绕卫星本体坐标系X_B轴逆时针旋转角度θ_lx，绕Y_B轴逆时针旋转角度θ_ly，各种因素引起的随时间变化的指向偏差Δθ_xl(t)和Δθ_xl(t)作为参数，利用公式(1)计算卫星本体坐标系激光脉冲方向矢量P_L(t)。将卫星在地面相对坐标系中的偏航角ω、俯仰角

和侧滚角κ作为参数，依据公式(2)计算得到地面坐标系下激光脉冲矢量P(t)。

c.根据卫星平台参数(卫星相对地面轨道高度H₀、卫星运行速度V、卫星轨道路径T_s)，激光发射重频设定R_L，激光发射坐标系在卫星本体坐标系下的坐标[X_L,Y_L,Z_L]，得到相对地面坐标系下的激光发射坐标[X_L0(n)，Y_L0(n)，Z_L0(n)]，然后结合步骤a得到的激光脉冲矢量P(t)，联立公式(11)，并结合z＝0计算出不考虑大气偏折时，地面坐标系的所有光斑脚点坐标[X(n),Y(n)]及激光入射指向角θ₀(x,y)。

d.根据大气环境参数(不同大气压位置的大气温度T_a，大气湿度H_a，位势高度G_a)，利用公式(8)得到不同高度的折射率N，依据折射率分布规律进行三次样条插值后，依据光线追迹方法得到横向漂移修正量Δs和角度偏离修正项Δθ，得到修正后相对地面坐标系下的光斑坐标[X(n),Y(n)]及入射指向角θ(n)。

e.根据光子计数激光测高仪发射系统参数(激光单脉冲发射能量E_t，激光脉冲均方根脉宽σ_s，发射光学系统效率η_t)，激光雷达与目标之间的单程大气透过率T_v，设定激光脉冲在空域分布为高斯分布，激光脉冲初始分布E₀(x,y)，利用式(4)计算单脉冲到达地面光斑位置[X(n),Y(n)]的能量分布E_i(x,y)。

f.设定地面目标反射率ρ，设定激光雷达接收系统参数，探测器的量子效率η_q，接收光学系统效率η_r，接收光学望远镜的有效口径A_r。通过式(5)计算出信号回波信号光子数分布N_s(x,y)，结合高程信息ξ(x,y)计算出的激光往返传播时间Δt(x,y)，累加同一时刻的回波信号，得到针对该目标的单位冲击响应W_s0(t)。设定激光单脉冲的时域分布为S(t)，则脚点的回波信号波形W_s(t)可通过W_s0(t)与S(t)卷积得到。

g.设定太阳高度角θ_h和方位角θ_a，结合随机地形的坡度角，接收系统接收光学系统的滤波片带宽△λ，探测器的量子效率为η_q，结合大气层外太阳光谱辐照度N_λ ⁰，依据式(6)，计算背景光噪声N_b。

h.设定单光子探测器的死区时间为t_d，时间分辨率为τ。结合步骤f得到的有效信号回波波形W_s(t)和步骤g得到背景光噪声N_b，计算出每个时间片的平均信号光子数K_sn(t,t+τ)，依据探测器响应模型计算出每个时间片的探测概率P_sn(t,t+τ)。使用蒙特卡洛仿真方法模拟探测器的响应过程及死区效应，产生0-1均匀分布随机数，若小于探测概率P_sn(t,t+τ)，则生成光子事件，从前至后，每响应一个光子事件，后续t_d时间内不响应光子事件，生成最终光子点云。

本发明中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (9)

1.一种全链路的光子计数激光测高仪点云的生成方法，其特征在于，包括：

步骤S1：建立激光发射坐标系、卫星本体坐标系和地面坐标系；

步骤S2：确定激光方向矢量、指向偏差和卫星在地面相对坐标系中的偏转，根据激光方向矢量、指向偏差和卫星在地面相对坐标系中的偏转，通过坐标转换关系得到地面坐标系下激光脉冲矢量；

步骤S3：根据激光脉冲矢量、激光发射点在相对地面坐标系下的坐标以及大地坐标系，计算出地面坐标系下所有光斑脚点的坐标和激光入射指向角；

步骤S4：获取光斑位置处的高程轮廓，对地面坐标系下所有光斑脚点的坐标和激光入射指向角进行修正，并根据修正后的光斑脚点的坐标、修正后的激光入射指向角以及光斑位置处的高程轮廓，计算光斑位置法向指向角；

步骤S5：获得光子计数激光测高仪的发射系统参数、激光雷达与目标之间的单程大气透过率以及激光脉冲在空域分布的初始分布，计算单脉冲到达地面光斑位置的能量分布；

步骤S6：获得地面目标反射率和激光雷达接收系统参数，根据单脉冲到达地面光斑位置的能量分布、目标反射率、激光雷达接收系统参数、光斑位置法向指向角与修正后的激光入射指向角以及目标高程轮廓，计算目标的单位冲击响应，再根据目标的单位冲击响应计算出光斑脚点的回波信号波形，目标高程轮廓为光斑位置处的高程轮廓；

步骤S7：根据激光雷达方程计算单位时间内探测器接收的由目标发射的背景光噪声；

步骤S8：根据光斑脚点的回波信号波形以及探测器接收的由目标发射的背景光噪，利用单光子探测器进行响应探测，生成全链路的光子计数激光测高仪点云。

2.如权利要求1所述的方法，其特征在于，步骤S1具体包括：

以激光发射的参考起点为原点，激光出射方向为Z轴正方向，Y轴垂直于卫星轨道平面，右手准则确定X轴建立激光发射坐标系O_L-X_LY_LZ_L；以卫星质心为原点，X轴指向卫星飞行方向，Z轴指向天底方向，Y轴垂直于卫星轨道平面建立卫星本体坐标系O_B-X_BY_BZ_B；以WGS84坐标系为基础建立大地坐标系(O-XYZ)。

3.如权利要求1所述的方法，其特征在于，步骤S2具体包括：

步骤S2.1：将激光方向矢量设置为绕卫星本体坐标系X_B轴逆时针旋转θ_lx、绕Y_B轴逆时针旋转θ_ly的单位矢量，指向偏差为Δθ_lx(t)和Δθ_ly(t)，根据公式(1)计算卫星本体坐标系激光脉冲方向矢量P_L(t)，

式(1)中，Δθ_lx(t)和Δθ_ly(t)是由各种因素引起的指向偏差，是随时间变化的量；

步骤S2.2：获取卫星在地面相对坐标系中的偏转，偏转包括偏航角ω、俯仰角

和侧滚角κ，根据公式(2)计算得到地面坐标系下激光脉冲矢量P(t)，

其中，ω、

κ分别为卫星定姿系统测定的卫星在地面相对坐标系中的偏航角、俯仰角和侧滚角，M(ω),

和M(κ)为坐标旋转矩阵。

4.如权利要求1所述的方法，其特征在于，步骤S3具体包括：

步骤S3.1：获得卫星平台参数：卫星相对地面轨道高度H₀、卫星运行速度V、卫星轨道路径T_s，；

步骤S3.2：获得激光发射重频设定R_L，激光发射坐标系在卫星本体坐标系下的坐标；

步骤S3.3：根据卫星平台参数、激光发射重频设定R_L，激光发射坐标系在卫星本体坐标系下的坐标，计算出地面坐标系激光发射坐标[X_L0(n)，Y_L0(n)，Z_L0(n)]；

步骤S3.4：根据地面坐标系激光发射坐标以及激光脉冲矢量P(t)，计算出地面坐标系的所有光斑脚点坐标[X₀(n),Y₀(n)]以及激光入射指向角θ₀(x,y)。

5.如权利要求1所述的方法，其特征在于，步骤S4具体包括：

步骤S4.1：获取大气参数：不同大气压位置的大气温度T_a、湿度H_a、位势高度G_a；

步骤S4.2：根据大气参数，构建大气折射偏移修正模型，并利用大气折射偏移修正模型对地面坐标系下所有光斑脚点的坐标和激光入射指向角进行修正，得到相对地面坐标系下的修正后的光斑坐标[X(n),Y(n)]及入修正后的射指向角θ_L(n)，其中大气折射偏移修正模型中大气折射率如公式(3)：

其中，k₁(λ)和k₂(λ)是与波长λ有关的函数，T为温度，P_d和P_w、Z_d和Z_w分别为干空气和水汽的分压强与压缩率；

步骤S4.3：获取光斑位置处的高程轮廓ξ(x,y)，根据修正后的光斑脚点的坐标、修正后的激光入射指向角以及光斑位置处的高程轮廓，计算光斑位置法向指向角θ_n(x,y)。

6.如权利要求1所述的方法，其特征在于，步骤S5具体包括：

步骤S5.1：获取光子计数激光测高仪发射系统参数：激光单脉冲发射能量E_t、激光脉冲均方根脉宽σ_s和发射光学系统效率η_t；

步骤S5.2：获得激光雷达与目标之间的单程大气透过率T_v、激光脉冲在空域分布的初始分布，根据公式(4)计算单脉冲到达地面光斑位置[X(n),Y(n)]的能量分布E_i(x,y)：

E_i(x,y)＝η_tT_vE₀(x,y) (4)

其中，激光脉冲在空域分布为高斯分布，E₀(x,y)为激光脉冲初始分布，η_r为光学系统的传输效率，T_v为星载激光测高仪与目标之间的大气透过率。

7.如权利要求1所述的方法，其特征在于，步骤S6具体包括：

步骤S6.1：获得地面目标反射率ρ，激光雷达接收系统参数：探测器的量子效率η_q、接收光学系统效率η_r、接收光学望远镜的有效口径A_r；

步骤S6.2：根据公式(5)计算出信号回波信号光子数分布Ns(x,y)，

其中，η_r为光学系统的传输效率，hv为单个光子能量，A_r为接收光学系统的有效口径，T_v为星载激光测高仪与目标之间的大气透过率，ρ为随机地形表面反射率，H₀是卫星轨道高度，θ_g为激光入射指向角与光斑位置法向指向角之间的夹角；

步骤S6.3：根据信号回波信号光子数分布和目标高程信息计算出的激光往返传播时间Δt(x,y)，累加同一时刻的回波信号，得到针对该目标的单位冲击响应W_s0(t)；

步骤S6.4：获得激光单脉冲的时域分布为S(t)，将目标的单位冲击响应W_s0(t)与激光单脉冲的时域分布为S(t)进行卷积，得到脚点的回波信号波形W_s(t)。

8.如权利要求1所述的方法，其特征在于，步骤S7具体包括：

获得太阳高度角θ_h和方位角θ_a，结合随机地形的坡度角、接收系统接收光学系统的滤波片带宽△λ、探测器的量子效率η_q以及大气层外太阳光谱辐照度N_λ ⁰，根据激光雷达方程，计算背景光噪声N_b，计算公式为(6)：

其中，N_λ ⁰为与波长相关的大气层外太阳光谱辐照度，△λ表示接收光学系统的滤波片带宽，R为地球半径，hv为单个光子能量，ρ为随机地形表面反射率，θ_r为接收单元的视场角，A_r为接收光学系统的有效口径，η_r为光学系统的传输效率，Ψ为太阳照射方向与目标表面法向方向的夹角，T_v是激光雷达与目标之间的双程大气透过率。

9.如权利要求1所述的方法，其特征在于，步骤S8具体包括：

步骤S8.1：获得单光子探测器的死区时间t_d、时间分辨率τ；i

步骤S8.2：结合有效信号回波波形W_s(t)和背景光噪声N_b，计算出每个时间片的平均信号光子数；

步骤S8.3：根据探测器响应模型计算出每个时间片的探测概率P_sn(t,t+τ)，其中，探测器响应模型为公式(7)：

P_sn(t,t+τ,λ>0)＝1-exp(-K_sn(t,t+τ)) (7)

其中，K_sn(t,t+τ)为时间片τ内的平均信号光子数；

步骤S8.4：使用蒙特卡洛仿真方法模拟探测器的响应过程及死区效应，产生0-1均匀分布随机数，当小于探测概率P_sn(t,t+τ)，则生成光子事件，从前至后，每响应一个光子事件，后续t_d时间内不响应光子事件，生成最终光子点云。

Images