星载激光测高仪脚点坐标解算方法及机载验证*

韩晓爽，潘 超，刘宇哲，王丽东，戴学兵，赵一鸣

（北京遥测技术研究所 北京 100094）

引 言

激光高度计卫星的研制已有近30 年的发展历史。从1994 年美国NASA 首次将激光测高仪搭载在航天飞机进行对地观测探索开始，火星轨道激光高度计MOLA（Mars Orbiter Laser Altimeter）、搭载于ICESat-1（Ice，Cloud，and land Elevation Satellite）卫星的地球科学激光测高仪系统GLAS（Geoscience Laser Altimeter System）、水星激光高度计MLA（Mercury Laser Altimeter）及月球轨道激光高度计LOLA（Lunar Orbiter Laser Altimeter）分别在陆地高程、森林、冰盖等探测中取得了丰富的成果。2018 年发射的ICESat-2 卫星搭载更为先进的地形激光高度计载荷ATLAS（Advanced Topographic Laser Altimeter System），率先采用6 波束激光、微脉冲、单光子探测体制，凭借其精确的测高能力，在山岳、冰川及冰帽、地表地貌、河流及湖泊深度、近海水深、海洋潮汐、高纬地区大地水准面、海洋重力场等地球科学的多个领域均广泛应用。NASA 规划的最新一代星载全球高分辨率成像测绘激光雷达系统LIST（Lidar Surface Topography）采用线阵推扫模式，具备多达1000 束激光的探测能力，可提供全球三维地形数据，预计2025 年发射。

多年来，在国家各类项目的支持下，我国星载激光雷达技术也得到了长足发展，已经具有较深厚的技术储备。从“嫦娥”探月系列卫星上搭载的激光测高仪开始，搭载了激光雷达载荷的“资源三号”02星、“高分七号”已经分别于2016 年和2019 年成功发射并在轨应用，搭载了星载多波束激光雷达载荷的“陆地生态系统碳监测卫星”也将于2021 年发射，实现对全球植被高度与大气分布的协同探测。多家单位研制的机载多波束三维成像激光雷达类载荷也开展了飞行试验。其中，由北京遥测技术研究所自主研制的机载大光斑激光雷达于2017 年分别在石家庄、张家界开展了多架次飞行试验，实现了星载多波束激光雷达载荷的机载验证与反演算法优化。

在星载激光测高仪激光地面脚点坐标解算方面，文献[1-4]介绍了ICESat-1/GLAS 的脚点定位模型和误差。ICESat-1/GLAS 和ICESat-2/ATLAS 的ATBD（Algorithm Theoretical Basis Document）报告分别对脚点定位流程及误差有较为详细的介绍。文献[5-7]也较为系统地阐述了卫星激光测高严密几何模型及精度分析，同时，文献采用“资源三号”02 星的激光测高数据对其进行了处理和验证。

本文对单波束、全波形和多波束、单光子两种探测体制激光载荷的坐标解算方法进行分析，在上述方法的指引下，构建机载平台坐标转换模型及方法，并通过机载飞行试验进行了验证。

1 星载激光雷达坐标解算方法

星载激光雷达地面脚点坐标的解算方法可基本概括为：通过激光发射和接收的时间差t，以及激光传播速度c，求得测距值，再结合卫星搭载的GPS 和星敏感器获得卫星位置和姿态信息，求得激光地面点在目标坐标系下的三维坐标。

根据卫星实际情况，星载激光雷达地面脚点坐标解算涉及瞬时激光束坐标系、测量参考坐标系（光学平台坐标系）、卫星本体坐标系、国际天球参考框架（ICRF）下的J2000 坐标系和国际地球参考框架（ITRF）下的WGS84 坐标系五个坐标系[8]。

转换流程如图1 所示[9]，其中下标SC、ICRF 表示卫星本体坐标系，COM 为卫星质心，Ref 为激光参考点，Laser Spot 为激光地面脚点，Δref为激光参考点和卫星质心的位置偏移量，h为激光传输距离。由RSPOT=RSC+Δref+h可求得激光地面脚点在ICRF 框架下的坐标，再将其转换至ITRF 参考框架下的坐标。获得足印点的三维坐标（X,Y,Z）后，也可用大地椭球坐标系（B,L,H）来表示。

在星载激光雷达地面坐标点的解算过程中，基本沿用以上原理进行坐标转换。根据测高体制、波束分布的不同，在具体的转换流程中也存在区别。下文分别对ICESat-1、ICESat-2的坐标解算方法进行分析，探寻不同载荷体制下坐标解算过程中的异同。

1.1 ICESat-1/GLAS

ICESat-1 于2003 年1 月发射成功，其搭载的GLAS 是目前较为成功的对地星载激光测高系统。GLAS为单波束、全波形测高系统，脉冲频率40Hz，沿轨方向足印间隔170m，地面足印大小55m～110m，其科学探测目标为南北极冰盖监测、陆地及水文地貌监测、植被高度监测和云-气溶胶高度分布监测。公开发布的数据产品包括GLA0-GLA15 在内的0 级、1 级和2 级数据产品。

GLAS 的激光脚点解算方法符合星载激光雷达坐标解算流程，但过程有所简化。采用ANC04 惯性坐标系旋转矩阵、ANC08 卫星精密定轨（POD）数据、ANC09 激光指向和精密定姿（PAD）数据、ANC25时间转换数据，与计算得到的激光测距值联合解算，得到激光地面点在ITRF 框架下的坐标。其计算流程如图2 所示。

图1 星载激光雷达地面脚点解算流程Fig.1 Spaceborne laser altimeter footprint coordinate calculation process

图2 ICESat-1/GLAS 地面坐标解算流程Fig.2 Footprint coordinate calculation process of ICESat-1/GLAS

由于GLAS 发射波形可近似为高斯脉冲，激光经地表反射后的回波波形可视作一个或多个高斯脉冲的叠加，采用波形分解的方法对回波波形进行拟合，从而求得激光单向传输的距离。

1.2 ICESat-2/ATLAS

已于2018 年9 月15 日发射的ICESat-2 卫星是ICESat-1 卫星的延续，搭载了更先进的地形激光测高系统载荷ATLAS。ATLAS 为非扫描推帚式、高重频、微脉冲、光子计数式的多光束激光测高载荷，轨道高度500km，轨道倾角92°，重访周期91 天。ATLAS 激光脉冲频率达10kHz，激光脚点沿轨方向间距0.7m，地面足印直径为17m；采用6 波束探测，其中3 波束能量较强，3 波束能量较弱，能量比为4:1；每组强激光脚点之间的垂轨距离3.3km，沿轨距离2.5km，强-弱激光脚点之间的垂轨距离为90m，如图3 所示[10,11]。相较GLAS 载荷，ATLAS 载荷可提供空间分辨率和平面定位精度更高的数据产品。

ATLAS 的激光脚点解算方法总体上与GLAS 相似，采用ATL02 数据、ANC03 惯性坐标系旋转矩阵、ANC04 卫星精密定轨（POD）数据、ANC05 激光指向和精密定姿（PAD）数据联合解算，得到激光地面点在ECF 坐标系下的坐标，其实现流程如图4 所示。

图3 ICESat-2 激光地面足印示意图Fig.3 Laser footprints of ICESat-2

图4 ICESat-2/ATLAS 地面坐标解算流程Fig.4 Footprint coordinate calculation process of ICESat-2/ATLAS

基于脉冲时间飞行法TOF（time of flight）的光子计数激光雷达通过记录起始脉冲信号与目标返回信号的时间差来测量距离。设激光光源发射脉冲信号的时刻为tT，脉冲信号传播至目标表面经漫反射、接收光学系统收集至探测器的时刻为tR，介质中的光速为c，则理想情况下探测距离为

针对光子计数体制的特点，载荷ATLAS 的地面脚点坐标解算方法也进行了相应的变化。在ATL02数据中，将N个光子的数据合并为一组进行处理，每组光子对应时间约为5ms，对应沿轨距离约为40m，通过解算获得所有N个光子的地面脚点坐标。后续在进行大气参数和大气延迟校正时，为节省存储与运算资源，从N个光子中选取第j个光子，其指向向量作为后续校正所用向量[12,13]。

具体解算步骤如下：

⑥利用ANC03 数据得到光子到达地面点时刻tB(i)下ECI 坐标系向ECF 坐标系的旋转矩阵，并计算激光地面点在ECF 坐标系下的坐标向量。

⑦将激光地面点在ECF 坐标系空间直角坐标系下的坐标转化为大地坐标系下的坐标形式。

2 机载激光雷达坐标解算方法

机载激光雷达坐标解算模型与星载模型较为相似，通过发射点到地面反射点距离、惯性导航系统（INS）测定的飞机姿态参数及GPS 测定的飞机精确位置联合解算获取地面点在目标坐标系下的三维坐标。星载平台通过星敏感器或陀螺仪定姿，机载激光雷达通过惯导定姿，在坐标转换过程中不涉及天球坐标系，不需考虑岁差、章动、极移等现象。参考激光测高卫星定位模型和ICESat-1/GLAS 公布的处理流程，可构建机载平台坐标转换模型，进行激光地面脚点坐标解算，以下分别对机载激光雷达涉及的坐标系统和解算方法进行介绍。

2.1 坐标系统

坐标系统包括瞬时激光束坐标系、激光扫描参考坐标系、惯性平台参考坐标系、导航坐标系和大地坐标系。其定义分别如下[14]：

①瞬时激光束坐标系：定义激光发射参考点为原点O，X轴为飞行方向，Z轴为瞬时光束方向，O-XYZ构成右手系。

②激光扫描参考坐标系（L 系）：定义激光发射参考点为原点O，XL轴为飞行方向，ZL轴为扫描系统的零点方向，O-XLYLZL构成右手系。

③惯性平台参考坐标系（I 系）：定义惯性平台中心为原点O，XI轴、YI轴ZI轴方向按照IMU 参考标架定义，O-XIYIZI构成右手系。

④导航坐标系（N 系）：又称为当地切平面坐标系，GPS 天线相位中心为原点O，XN轴指向真北，YN轴指向东，ZN轴为铅直向下方向，O-XNYNZN构成右手系。

⑤大地坐标系（E 系）：采用WGS84 坐标系，地球质心为原点O，Z轴指向国际时间局（BIH）1984.0定义的协议地级（CTP）方向，XE轴指向BIH1984.0 的协议子午面和CTP 赤道的交点，O-XEYEZE构成右手系。

2.2 坐标转换过程

机载激光雷达对地定位中的坐标转换顺序为瞬时激光束坐标系→激光扫描参考坐标系（L 系）→惯性平台参考坐标系（I 系）→导航坐标系（N 系）→大地坐标系（E 系）[15,16]，各坐标系的定义2.1 节已进行阐述，具体转换过程如下。

①瞬时激光束坐标系转L 系

在激光扫描参考坐标系下定义扫描角θ为激光束与ZL的夹角，逆时针转动方向为正，如图5 所示。

当激光测距值为ρ时，地面点在L 系下的坐标可表示为

图5 激光扫描参考坐标系Fig.5 Laser scanning reference system

② L系转I 系

可通过测量得到L 系和I 系坐标原点间的距离，即惯性平台中心与激光发射参考点的偏移量为[ΔXLΔYLΔZL]T；坐标轴之间的夹角为偏心角α、β、γ，即惯导与激光扫描仪在安装过程中三个坐标轴产生的固定轴线偏差转角。地面点在I 系下的坐标可表示为

③I 系转N 系

N 系和I 系坐标原点间的距离为偏移量[ΔXNΔYNΔZN]T，坐标轴之间的夹角可由IMU 记录的三个姿态角航向角H(heading)、俯仰角P(pitch)和侧滚角R(roll)表示，地面点在N 系下的坐标可表示为

④N 系转E 系

GPS 系统可以实时获得天线相位中心即导航坐标系原点的参心大地坐标（B,L,H）（纬度，经度和椭球高），[XgpsYgpsZgps]T为GPS 在WGS84 坐标系下的坐标，可利用下面的公式进行计算：

其中，e 为第一偏心率，N为卯酉圈曲率半径，，a 为WGS84 椭球的长半轴。

地面点在WGS84 下的坐标可表示为

根据上述的坐标转换公式，最终可求得激光地面点在WGS84 坐标系下的坐标表示为

其中，[XoffsetYoffsetZoffset]T为激光参考点到GPS 天线中心的偏移量。

通过上述坐标系统的转换，可以得到每个激光地面点在地心坐标系WGS84 下的三维坐标值。可以根据需要将激光地面点的空间直角坐标转换为大地坐标。

3 机载验证系统——大光斑激光雷达

机载大光斑激光雷达是为满足碳监测卫星多波束激光雷达载荷的算法验证需求而研制的机载验证载荷，采用单波束、全波形体制，通过将机载大光斑激光雷达与航空相机、航空惯导系统、高稳平台、控制及存储系统等设备的集成和优化，形成机载大光斑激光雷达系统，以40Hz 的重频发射激光脉冲进行地表探测[17,18]。系统参数指标如表1 所示，系统组成如图6 所示。

表1 机载大光斑激光雷达系统参数指标Table 1 Parameters of airborne large spot LiDAR system

4 机载飞行验证试验数据结果

机载大光斑激光雷达系统于2017年4 月上旬，在石家庄栾城机场附近区域开展了首次机载飞行试验，并开展了点云激光雷达同步试验比对、地面外业试验交叉比对与精度评估。此次飞行试验包括两个架次的校飞试验，获取了检校区与测区共计10 条航线，约300GB 的有效数据。系统再次于2017 年12 月下旬，在张家界荷花机场附近区域开展了两个架次的飞行试验。此次飞行试验航线总长超过980km，航高超过3000m，航时总计11 小时，获取了检校区与测区共计四个测区的有效数据。

图6 机载大光斑激光雷达系统构成Fig.6 System configuration of airborne large spot LiDAR system

对机载大光斑激光雷达系统于2017 年12 月25 日在张家界附近采集到的一段数据进行解算，可得到地面激光点的轨迹如图7 所示，部分数据如表2 所示。

图7 部分地面激光点轨迹Fig.7 Partial laser footprint track

表2 经解算得到的部分激光地面点坐标Table 2 Partial calculated laser footprint coordinates

获得回波信号数据后，通过高斯滤波方法对信号进行去噪处理，再用最小二乘曲线拟合法或最大值检测法可以得到该点位测距值[19,20]。本论文通过计算回波信号的二阶差分获取峰值，采用公式如下：

其中，w为波形数据，有：

利用find 函数获得满足式的元素的序号，k值所在位置即为回波峰值所在的点位，当回波信号有多个峰值时，通常选取回波强度最大的位置作为峰值所在点。利用峰值位置计算出激光测距值后，根据第

2.2 节介绍的坐标转换方法可解算出对应的激光地面点坐标。

以表1 第8 行数据为例，此时UTC 为12987.5，采集到的回波信号波形如图8 所示，通过地图定位可知该位置为树木，坐标解算结果与信号符合。表1 第25 行数据对应UTC 为13017.85，采集到的回波信号波形如图9 所示，通过地图定位可知该位置为裸地，信号强度过大，存在饱和现象，坐标解算结果与信号符合。

图8 激光回波信号数据与对应地面脚点坐标（UTC=12987.5）Fig.8 Reflected laser echo signal data and footprint coordinate（UTC=12987.5）

图9 激光回波信号数据与对应地面脚点坐标（UTC=13017.85）Fig.9 Reflected laser echo signal data and footprint coordinate（UTC=13017.85）

5 结束语

本文通过分析ICESat-1 和ICESat-2 星载激光测高载荷地面脚点坐标解算方法，比较了多波束、单光子体制对地激光三维测绘雷达与传统的单波束、全波形体制激光测高雷达在地面坐标解算过程的异同，并建立了机载激光雷达坐标解算方法与流程，通过自主研制的机载大光斑激光雷达与飞行试验，对坐标解算方法进行验证与分析。随着单光子探测技术的发展，光子计数激光雷达正在逐步应用于机载、星载对地观测领域，北京遥测技术研究所也正在进行机载光子计数激光雷达研制。通过对ICESat-2 地面坐标解算流程作进一步理论分析，改进并构建适用于机载光子计数激光雷达的坐标解算方法，并应用于机载飞行验证试验，是下一步工作的重点。