随着北京冬奥会的成功举办，中国滑雪运动爱好者与日俱增^[1]。为保证滑雪体验，雪道质量的维护成为滑雪场的重要工作内容^[2]。压雪车作为一种广泛应用于滑雪场的特种作业车辆，其主要作业环境为夜间低温山地，能见度极低，且通常存在陡坡与悬崖等地势，对雪道质量维护作业甚至是压雪车行驶安全提出了极大挑战^[3-4]。

随着环境感知技术的发展，基于激光雷达的雪道分析受到了越来越多的关注^[5-8]，德国凯斯鲍尔已经研发出基于激光雷达的雪道环境感知系统^[9]。中国关于滑雪场雪道感知与重构的研究较少，目前已有研究手段多依赖于无人机机载激光雷达。鉴于人工造雪和日常运营等因素，滑雪场雪道各区域的状态随时会发生变化。为实时准确地识别雪道，提高雪道平整作业的效率和维护精度，提升行驶的安全性，基于压雪车车载激光雷达的实时雪道信息获取工作尤为重要。现有地面分割算法研究主要针对城市公路环境，例如基于高程地图、平面拟合等算法，其主要目的是去除地面点，以加快后续如聚类、跟踪和控制等任务^[10-12]。然而，由于室外大型滑雪场多依山而建地形复杂，同一雪道上坡度连续变化，城市公路环境下常用的地面分割方法难以直接用于高山滑雪场的雪道分割^[13-14]。同时雪道分割与提取旨在获得更丰富的地面信息，为压雪车标注可行驶区域，提高工作效率。因此有必要根据高山滑雪场的地形特点，提出针对性的雪道分割方法。

然而雪场作业时多存在降雪工况，会导致激光雷达前方存在飘雪噪声，激光脉冲在撞击雪花后的反射信号存在误差，发生误检。为消除飘雪噪声引起的不利影响，通常采用滤波算法对激光雷达数据进行预处理，达到滤除噪声的目的^[15-16]。目前，激光雷达飘雪噪声去除主要是基于距离和基于深度学习的滤波算法，这类方法有较高的滤波精度，但计算成本较高，不适合高山滑雪场点云的实时处理与分析。

因此，本文采用非重复性扫描的半固态式激光雷达，搭建了高山滑雪场夜间感知系统。提出了基于点云反射率分布模型的飘雪预处理方法，并基于扇形栅格对雪道坡度进行实时估计。根据激光点云在空间上的局部时域连续性特点，提出了基于时域波形阶跃值检测和包络方法，用于快速雪道分割计算，同时提出了适用于高山滑雪场的含雪道分割、特征分类匹配和特征约束的实时定位与建图(simultaneous localization and mapping，SLAM)方法。在张家口万龙滑雪场娇龙道与游龙道开展了测试，验证了所提出方法的有效性。

1 高山滑雪场夜间感知系统

车载激光雷达雪道环境感知系统框架如图 1(a)所示，载具为中国自主研发的国内首台室外大型压雪车SG400。平整车身前方雪道为压雪车夜间作业的主要任务工况之一，为最大限度获取车体前的雪道地形点云数据，车顶激光雷达与地面设计垂直距离h=2.3 m，向前倾斜角度α= 11°。压雪车在实际行驶与作业过程中振动剧烈，为减少车身振动对激光雷达硬件造成的影响，设计了悬臂减振支架如图 1(b)所示。山地滑雪场在雪季夜间温度最低可达-30 ℃，已超出激光雷达正常工作温度范围，设计了内置加热片的防寒外壳，以保证激光雷达在极寒环境下可正常工作。压雪车车载系统计算平台为NVIDIA Jetson AGX xavier(8核Carmel 64位ARM CPU、Volta GPU)。

压雪车的行驶速度通常在5~7 km/h之间，在作业时行驶速度会降低至4 km/h以下。雪道环境感知系统中的车载激光雷达Livox Mid-70的固定视角为70.4°，具有关注区域局部点云密度高，可靠性高，成本低等特点，可满足压雪车对前方环境的感知需求^[17]。Livox Mid-70激光雷达采用“拟人眼”的单线花瓣式扫描方式，激光点云密度在扫描中心处最高，并且随着点到扫描中心距离的增加而减少^[18]。通过两个棱镜的转速调整，Livox Mid-70实现了点云的非重复性扫描，可以在足够的累积时间情况下提供高分辨率扫描^[19]。在0.03、0.06、0.10 s的采样时间下，具体激光扫描累积路径如图 2所示，随着采样时间的增加累积扫描面积显著增大。

2 基于时域波形特征的雪道分割

Livox Mid-70的单束激光扫描所得花瓣状点云组的空间位置在时域上存在序列关系，即单束激光扫掠过单一平面时，所获取的花瓣状点组的坐标x，y，z随时间呈连续变化且较为平缓，不存在明显阶跃；扫掠过不同物体平面时，其在x，y，z坐标的一个或多个维度会出现明显的阶跃。基于这一特性，本文提出了一种基于时域波形阶跃值包络的雪道分割算法，算法流程如图 3所示。

2.1 基于点云反射率分布模型的飘雪预处理

激光雷达点云数据P包含一个物体在三维空间的位置x，y，z和目标物体的反射率R，P=[x, y, z, R]^[20-21]。经实测，空中飘浮雪颗粒对激光信号有很强的吸收作用，激光光强衰减导致反射率降低，依此特性将点云反射率作为划分飘雪和其他物体的重要依据^[22]。图 4展示了测试场地某段雪道的无降雪和降雪两种工况下连续200帧激光雷达数据的反射率分布，相比于无降雪工况，降雪工况下测量点数N在反射率R∈[0, 4]范围内大幅增加，导致雪道测量点R∈[150, 153]急剧减少。根据飘浮雪颗粒低反射率特性与实测对比数据，可知降雪工况下飘浮雪颗粒的反射率集中区间为R∈[0, 4]。由于飘雪集中在激光雷达前方，因此可将激光雷达前方2 m＜x＜6 m反射率R∈[0, 4]的点标记为飘雪噪声点。

图 5为激光点云扫描示意图，受飘浮雪颗粒遮挡的影响，无法获得真实的激光测距值。然而，直接过滤飘雪噪声，将导致扫描点空间位置在时域上的不连续性，影响后续基于时域波形阶跃特征的检测。因此，本文采用区间相邻点线性插值方法对过滤后的不连续雪点P_k的坐标进行插值拟合：

$ P_k=P_{m-1}+(k-m+1) \frac{P_{n+1}-P_{m-1}}{n-m+1} $

(1)

式中：P_m、P_n分别为飘浮雪颗粒遮挡的起始和终末点；k为遮挡点个数, k=1，2，3，…，取二者的相邻非雪点P_m-1和P_n+1相连直线的等间隔点的位置作为插值点的坐标。

2.2 基于扇形栅格的坡度估计

实验测试在张家口崇礼万龙雪场完成，其高级雪道设计的最大坡度θ_max =30°，平均坡度θ_mean=20°；中级道θ_max=25°，θ_mean=18°；初级道θ_max=10°，θ_mean=7°。因自然降雪和人工造雪导致雪道坡度变化较大，需要估计车前雪道坡度参与时域波形阶跃值检测计算。根据实际雪道情况，认为车前方10~40 m，与x轴正方向夹角-7°~7°的点云属于雪道点云，并以2 m为固定间隔，划分出雪道空间上共15个扇形栅格区间，如图 6所示。以各栅格区间内点云坐标平均值作为栅格代表点P，并计算相邻代表点之间的坡度角θ。

$ \theta=\frac{1}{n} \sum\limits_{n=1}^{14} \arctan \left(\frac{\bar{z}_{n+1}-\bar{z}_n}{\sqrt{\left(\bar{x}_{n+1}-\bar{x}_n\right)^2+\left(\bar{y}_{n+1}-\bar{y}_n\right)^2}}\right) $

(2)

2.3 点云时域波形阶跃值检测及区间包络

雪道的障碍物主要包括：两侧树林、服务区建筑群、围栏、广告牌、照明灯等静态障碍以及同雪道内其他作业压雪车。当激光扫描雪道或障碍的过渡点时，三维坐标时域波形上的表现为阶跃特征。经实际测量可知，相邻两点间z轴高程差Δz及x轴距离差Δx是区分雪道面和非雪道障碍过渡点的明显特征。基于雪道与障碍点云特征，分别提出了z轴、x轴时域波形阶跃值检测判据，以获得雪道和障碍物之间的点云时域波形阶跃过渡点。

2.3.1 基于z轴的时域波形阶跃值检测

如图 7(a)所示，以建筑障碍为例，作为一种几何连续的障碍物，激光光束从花瓣中心P₀处以顺时针开始扫描，当从雪道表面扫描至建筑表面时，点P₃和P₄之间的高程差Δz_b剧增，即点云的z轴时域波形呈现在分界线处P₃出现阶跃现象(如图 7(c)所示)。同理，当激光从建筑表面扫描至雪道时，点P₈和P₉之间的高程差Δz_e骤减，点云z轴时域波形上P₉也会出现阶跃的现象。此条件下x轴的时域波形无明显阶跃现象(如图 7(b)所示)。

z轴方向，选择当前点P_i与前、后各6个点，计算相邻点的高程差。点P_i的z轴阶跃判据如下：

$ D_z=\left\{\begin{array}{l} +1, \quad\left|\Delta z_{i-k}\right|>\frac{z_{\mathrm{th}}}{\cos \theta}, k=1, \cdots, 6 \\ 0 \end{array}\right. $

(3)

$ U_z=\left\{\begin{array}{l} +1, \quad\left|\Delta z_{i+k}\right|>\frac{z_{\text {th }}}{\cos \theta}, k=1, \cdots, 6 \\ 0 \end{array}\right. $

(4)

$ L_{z i}= \begin{cases}1, & U_z>4, D_z>4 \\ 0\end{cases} $

(5)

式中：z_th = 0.05为高程差判据阈值，θ为坡度估计值，D_z为当前点P_i时域前序点集权重统计，U_z为点P_i的后序点集权重统计，若D_z>4或U_z>4，则当前点P_i被标记为z轴时域上的阶跃点。

2.3.2 基于x轴的时域波形阶跃值检测

如图 8(a)所示，以远距离山体或地形凹坑的几何间断障碍为例，激光光束从花瓣中心P₀处以顺时针开始扫描，受坡度β和激光扫描角度的影响，花瓣从雪道P₂直接扫掠至障碍P₃处，随后花瓣扫描回转后从障碍P₇扫描回到了雪道P₈。图 8(b)、8(c)分别展示了点云的x轴和z轴时域波形，在时域波形上雪道点P₂到障碍点P₃处Δx_b剧增，在高障碍点的P₇到雪道点P₈处则发生了Δx_e骤减。z轴的时域波形上未发生较大阶跃。

如图 9(a)所示，以包含树木和地面点的树林为例，作为一种几何间断且不规则的障碍，当激光从雪道开始扫描至树林区域时，地面点P′₅和P′₇被树木遮挡，得到了表示树木的点P₅、P₇。在x轴时域波形的点P₃和P₉之间发生无序的阶跃变化，z轴波形存在类似微小幅度变化。

基于时域的x轴坐标阶跃值检测的相关判据如下。判断当前点P_i与其在时序上的前、后各4个点的x向距离阈值计算权重，并根据P_i前、后两点来判断当前点是否为跳跃阶跃点。对于当前点P_i，点云的x轴距离阶跃的相关判据如下：

$ S_x=\left\{\begin{array}{l} +1, \quad x_i-x_{i+2 j-1}>x_{\mathrm{th}} \cos \theta, j=1,2,3,4 \\ 0 \end{array}\right. $

(6)

$ \operatorname{Inf}_i=\left\{\begin{array}{l} \text { Ture }, \quad x_i-x_{i-1} \geqslant 0, x_i-x_{i+1} \geqslant 0 \\ \text { False } \end{array}\right. $

(7)

$ L_{x i}= \begin{cases}1, & S_x>6, \operatorname{Inf}_i=\text { True } \\ 0\end{cases} $

(8)

式中：x_th=0.05为x轴上的距离阶跃阈值，θ为坡度估计值，S_x为时序上x距离的权重统计，Inf为拐点判据，若S_x>6或当前点为拐点时，则该点被标记为x轴时域波形上的阶跃点。

基于时域波形阶跃特征，寻找雪道-障碍交界点并完成点云的雪道与非雪道分割。对单个花瓣区间内出现的波形阶跃值点云进行包络，识别非雪道障碍。通过计算激光扫描特征角φ确定单个花瓣区间(如图 5所示)，以初始阶跃点作为障碍起始点、末尾阶跃点作为障碍结束点进行区间包络，识别当前花瓣上的障碍。激光扫描特征角φ计算公式为

$ \varphi=\arccos \left(\frac{z}{\sqrt{x^2+y^2+z^2}}\right) $

(9)

3 含雪道分割的滑雪场实时定位与建图

为实现对压雪车前方环境的实时感知并将雪道分割结果可视化表达，基于Loam Livox算法的基本结构，提出了一种适用于高山滑雪场的激光雷达实时定位与建图算法，如图 10所示。基于点云反射率分布模型的飘雪噪声预处理，利用点云时域波形阶跃值包络方法，完成雪道分割。根据所得雪道与障碍标签，分别在对应的点集中进行特征提取：根据曲率的大小划分平面与边缘两类特征点。曲率c的计算公式为^[23-24]

$ c=\frac{1}{|S| \cdot\left\|X_{(k, i)}^L\right\|} \| \sum\limits_{j \in S, j \neq i}\left(X_{(k, i)}^L-X_{(k, j)}^L\right)|| $

(10)

式中：X_{(k, i)}^L为当前帧某点的坐标，S为当前点i相邻的前后n个点组成的集合，X_{(k, j)}^L为邻域集合S中的点j的坐标。将曲率c>0.100的点作为边缘点，c＜0.005的点作为平面点。通过特征分类，最终获得4类特征点：雪道平面点、雪道边缘点、障碍平面点及障碍边缘点。

在迭代位姿优化过程中，雪道平面点和障碍平面点在其局部平面特征地图中利用KD树寻找最近邻点并计算特征值，分析特征值判断是否为面特征，并根据式(11)构建点到面的残差r_p2p；同理，雪道边缘点和障碍边缘点根据式(12)构建点到线的残差r_e2e。

$ r_{\mathrm{p} 2 \mathrm{p}}=\frac{\left(P_{\mathrm{p}}-P_1\right)^{\mathrm{T}}\left(P_3-P_5\right) \times\left(P_3-P_1\right)}{\left|\left(P_3-P_5\right) \times\left(P_3-P_1\right)\right|} $

(11)

$ r_{\mathrm{e} 2 \mathrm{e}}=\frac{\left|\left(P_{\mathrm{e}}-P_5\right) \times\left(P_{\mathrm{e}}-P_1\right)\right|}{\left|P_5-P_1\right|} $

(12)

式中：P_p为平面点，P_e为边缘点。最后将残差r_p2p和r_e2e添加到姿态优化中，经迭代优化后得到点云最优位姿。最优位姿的连续积累可以用于更新压雪车的位置，跟踪其行驶轨迹。

为提升雪道SLAM的实时性，在建图过程中对特征点的数量进行约束。实验中发现，雪道点云以平面特征点为主，但平面特征点占比过高易引起位姿估计误差，并且容易对车载计算平台造成算力负担。因此根据雪场环境特点，将雪道平面特征点限制为曲率最小的i个点，在保留特征的基础上，降低了雪道大平面对点云位姿估计和建图带来的影响。

4 结果与分析

高山滑雪场具有场景大、特征少、坡度变化大等特点，为验证基于激光反射率的飘雪噪声滤波、基于时域波形阶跃值包络的雪道分割和基于雪道分割、特征分类与特征约束的SLAM方法的有效性，本文选取了张家口万龙滑雪场的两条典型赛道：游龙道与娇龙道开展实验，如图 11所示。

4.1 飘雪预处理结果

图 12展示了压雪车在滑雪场上坡无降雪工况和降雪工况行驶时激光雷达采集的某单帧点云数据。通过对比点云特征，相比于无降雪工况，图 12(d)中由于飘雪颗粒的影响，在雷达前方2~ 6 m处聚集了大量低反射率点。利用飘雪颗粒的反射率分布模型识别出激光雷达前方2~6 m的低反射率点云为飘雪颗粒噪声。如图 13所示，单帧点云中雪颗粒可以准确地被标记出来，基于反射率分布的雪颗粒标记方法效果良好。

由于激光雷达的花瓣式扫描特性，时域上的点云在各坐标上应该具有连续的周期性变化特征。图 14(a)为降雪工况下激光扫描的单帧点云x坐标的原始时域波形图，雷达前方的雪颗粒噪声导致测量点偏离坐标值。同时，由于雪道坡度变化较大，雷达随车身俯仰角度变化将导致点云数据中存在一定数量无回波信号点，也会导致坐标值偏离。利用相邻点线性插值解决时域波形的不连续问题，降噪插值后的时域波形如图 14(b)所示。

4.2 雪道分割结果

在完成点云预处理后，可以得到如图 15(a)、15(c) 所示的坐标x、z在时域0 s＜t＜0.1 s上的连续波形。依据x和z坐标的时域波形阶跃值判据，分别对其时域上的阶跃值位置进行识别标记(L_x，L_z)，如图 15(b)、15(d)所示。在x方向上，每一个峰值所在的邻近区间内都存在若干阶跃值，表明车前距离较远位置处存在障碍。在0.02 s＜t＜0.04 s区间内，存在明显的障碍开始点和结束点，符合几何连续障碍的特点。在t=0.08 s附近区间，连续出现的阶跃值表明存在几何不连续的障碍。结合z坐标相关的时域波形特征，可知在0 s＜t＜0.03 s，两轴阶跃值判据具有相似检测结果，均能在相近的位置识别出障碍物所产生的点云时域波形阶跃。而对于t>0.04 s的区间，二者对阶跃值的检测并不相同，可以实现对障碍物识别的互相补充。

对x轴和z轴波形阶跃值结果取并集L后，得到了如图 16(a)所示的阶跃标记点集。并通过单帧点云的特征角变化(如图 16(b)所示)，得到了各单个花瓣的扫描时间范围。基于时域波形阶跃值检测结果，以单个花瓣周期内的初始阶跃点作为包络起点、末尾阶跃点作为包络终点进行包络，得到如图 16(c)所示的阶跃值包络图，包络内的点均为障碍点。

实际雪道由于坡度和周围特征复杂多变，为验证雪道分割算法的有效性，分别对双侧障碍、转向工况和不同坡度分岔雪道3个典型雪道地形进行雪道分割检验。各个雪道基于点云时域波形阶跃值包络方法的分割效果分别如图 17~19所示。图 17为压雪车行驶在双侧均存在障碍物的雪道分割效果，此时雪道坡度在θ=17°左右波动。由于雪道坡度变化的影响，部分雪道被误判为障碍(如图 17(a)所示)。图 17(b)为实时坡度估计两种情况下的雪道分割结果，基于坡度估计的分割效果体现出了更好的地形适应性。图 17(b)中雪道两侧障碍具有两种不同特征：其中，y>0区间内的障碍点分布较为规律，符合侧方建筑物障碍的特点；y＜0区间内的障碍点的分布更趋离散，符合非连续障碍的特点。在两种障碍特征同时存在情况下，提出的雪道分割算法具有较好的准确性。

高山滑雪场雪道具有转弯角度大的特点，压雪车向y>0方向转向工况下的雪道点云分割结果如图 18所示。从分割结果可知，在60 m＜x＜80 m区间存在前方离散障碍，符合实际转弯工况下的障碍特点。压雪车行驶作业时，存在雪道分岔与交汇工况，因此分岔路段的雪道可行驶区域识别就尤为重要。如图 19(a)所示为雪道分岔工况的雪道分割结果。图 19(b)、19(c)分别为点云在xoz、xoy两个平面上的投影，雪道1坡度约为θ=34°，雪道2坡度约为θ=16°，在两条雪道坡度差Δθ=18°的地形条件下，雪道分割算法仍能准确分割两条雪道点云。综上所述，在3种高山滑雪场的典型地形条件下，所提出的雪道分割算法可以实现标记障碍物和分割压雪车可行驶雪道点云的功能。

并从准确率A和运行时间t两方面对算法进行评估:

$ A=\frac{N_{\mathrm{TP}}+N_{\mathrm{PP}}}{N_{\mathrm{TP}}+F_{\mathrm{TP}}+N_{\mathrm{PP}}+F_{\mathrm{PP}}} $

(13)

式中: N_TP为提取正确的雪道点云数量，N_PP为提取正确的非雪道障碍点云数量，F_TP为提取错误的雪道点云数量，F_PP为提取错误的非雪道障碍点云数量。通过对在万龙滑雪场多条雪道上采集的点云数据进行标记，得到雪道分割的平均准确率A为98.54%，平均运行时间t为2.36 ms。

4.3 高山滑雪场SLAM建图结果与分析

图 20为降雪条件下，压雪车在娇龙道行驶作业时的SLAM结果。图 20(a)、20(b)分别为无飘雪噪声处理和基于点云反射率分布模型去除飘雪噪声后的雪道建图结果。图 20(c)、20(d)分别为图 20(a)、20(b)区域1的局部放大，经对比可以明显看出，飘雪噪声处理方法在SLAM建图过程中能够有效过滤激光雷达前方聚集性飘雪颗粒。建筑障碍的区域内建筑点云受飘雪颗粒的影响，配准误差增大，导致模糊和错位的现象(如图 20(a)，区域2所示)，而过滤飘雪噪声后建筑障碍点云构图的错位问题得到明显改善(如图 20(b)，区域2所示)。

图 21、22分别为针对娇龙道和游龙道的基于时域波形阶跃值包络算法的含特征分类和约束的雪道分割SLAM效果。图 21(a)可知，雪道分割算法在坡度变化稳定的娇龙道上可以准确划分出雪道及护栏、建筑、树木等障碍物，实现雪道可行驶区域的分割。图 21(b)为建筑障碍的局部放大图，在本文提出的SLAM算法下，建筑物可以准确被标注为障碍点，且细节构建良好。图 22展示了未进行特征分类匹配和进行特征分类匹配的游龙道雪道分割效果。由图 22(a)、22(b)可知，雪道分割算法在雪道坡度大、存在连续大转弯以及非连续性障碍特征少的游龙道上仍能准确划分出雪道及树木等障碍。如图 22(c)、22(d)所示，在具有转角变化大的转弯处，相比于未进行特征分类匹配的建图算法，特征分类匹配的建图效果更加精确，模糊程度明显下降。

图 23为针对游龙道的原方法(Loam Livox)、基于雪道分割特征分类和基于特征分类与约束的建图运行时间对比。本文提出的基于雪道分割特征约束的分类匹配方法在运行时间上明显低于原方法和直接特征分类匹配方法。在高山滑雪场的定位建图中，特征约束的匹配方法可以有效地提高计算效率，且运行时间均低于20 ms，有效地增强了建图的实时性。

5 结论

1) 利用飘雪噪声低反射率和集中分布于雷达前方的特点对其准确标记，并通过线性插值方法去除噪声对时域波形连续性的干扰，有效解决了由于飘雪噪声引起的地图重复构建和错位构建问题。

2) 通过娇龙道和游龙道两条雪道进行实验验证，测试结果表明含坡度实时估计的时域波形阶跃值检测包络的雪道分割算法能够消除高山滑雪场的坡度影响，准确分割出不同坡度变化下的雪道，对单帧点云数据处理的平均耗时为2.36 ms，平均分割准确率可达98.54%。

3) 通过实验验证，改进的SLAM方法计算效率对比原方法有明显提升，运行时间均降至20 ms以内，有效地增强了建图的实时性。