A New Method of Ski Tracks Extraction Based on Laser Intensity Information

Abstract

At present, the mainstream laser point cloud classification algorithms are mainly based on the geometric information of the target. Nevertheless, if there is occlusion between the targets, the classification effect will be negatively affected. Compared with the above methods, a new method of ski tracks extraction using laser intensity information based on target reflection is presented in this paper. The method can complete the downsampling of the point cloud datasets of ski tracks under the condition that the information of the target edge is complete. Then, the clustering and extraction of ski tracks are effectively accomplished based on the smoothing threshold and curvature between adjacent point clouds. The experimental results show that, different from the traditional methods, the composite classification method based on the intensity information proposed in this paper can effectively extract ski tracks from the complex background. By comparing the proposed method to the Euclidean distance method, the clustering segmentation method, and the RANSAC method, the average extraction accuracy is increased by 16.9%, while the over extraction rate is reduced by 8.4% and the under extraction rate is reduced by 8.6%, allowing us to accurately extract the ski track point cloud of a ski resort.

1. Introduction

With the successful holding of the 24th Beijing Winter Olympic Games in 2022, people are progressively interested in participating in winter sports, especially skiing [1,2]. In the meantime, researchers and scientific research are becoming increasingly interested in skiing. In recent years, a brand-new high-tech method to collect high-precision 3D point cloud data and extract the needed information of ski tracks by LiDAR has been developed [3,4,5,6]. For the managers of ski resorts, the relevant measurement and statistics (such as thickness, slope et al.) is an important part of daily maintenance. Therefore, a key step is to effectively separate the ski tracks from the background by semantic segmentation of 3D point cloud datasets.

In addition, as the key technology of point cloud processing, the semantic segmentation of laser point cloud is widely used in artificial intelligence, reverse engineering, target reconstruction, robot, medical image analysis, and other related fields [7,8]. The point cloud semantic segmentation classifies three-dimensional points according to their local features, classifies points with the same or similar attributes into a set of disjoint points, and points in the same set that have similar features [9]. The semantic segmentation, which is helpful in scene analysis, such as location and object recognition, classification and feature extraction, and object recognition, is one of the key links in point cloud data processing.

The current recognition method of the target classification for 3D datasets research mainly uses the coordinate information of point clouds [10,11,12]. The geometric dimension of the target is calculated by the coordinate information of the point cloud; thus the target classification and recognition can be realized. Guo M et al. divided the point clouds into independent targets and classified the targets by extracting the volume characteristics of each target [13]. Xu X et al. divided the three-dimensional datasets into slices according to height and extracted the length and width features of each slice to detect pedestrian targets [14]. The classification method based on the coordinate information of point clouds is easily affected by the occlusion of target, and if the targets have similar geometric characteristics, they cannot be distinguished effectively. This situation makes these methods unable to effectively separate the ski tracks from the complex background. The LiDAR can accurately record coordinate information and reflection intensity of an object, which can efficiently reflect the reflection ability of the object, and has great application prospects in semantic segmentation for multiscale targets. Generally, the reflection intensity of the target is affected by many factors such as atmospheric conditions, incident azimuth, and emission distance, resulting in different reflection intensities of objects [15,16,17,18,19]. However, due to artificial snow and its high hardness at ski resorts, the laser reflection is similar. The research on how to extract the parameters of the object characteristics from the intensity information is of great significance in strengthening the utilization of the intensity information and broadening target classification methods.

In this paper, a composite ski tracks extraction algorithm based on laser point cloud intensity information, which can effectively extract the ski tracks under the complex background and provide convenience for the daily operation and maintenance of ski resorts, is proposed.

2. Study Area and Data

The LiDAR datasets used in this paper were acquired through the unmanned aircraft LiDAR System (RIEGL VUX-1 UAV System). The collection site is located in Zhangjiakou Wanlong ski resort, the most famous ski resort in China, near Yanqing District, Beijing (249 km away from Beijing, 2110.3 m above sea level, more than 200,000 square meters).

The collection time was 12:00 to 14:00, and the temperature of the ski resort was −18 °C. Figure 1a shows the environmental conditions of the Wanlong ski resort at the time of collection, Figure 1b is an aerial picture taken during acquisition, and Figure 2 is an introduction to the system architecture.

The seven ski tracks of Wanlong ski resort are all made of extremely hard artificial snowPlease verify meaning is retained. The environment of the snow field is relatively complex, including not only five main snow tracks, but also two training tracks, as well as trees, riprap, etc. A DJI M600 PRO UAV was used as the flight platform with LiDAR because of its excellent stability and flight performance, and the fact that its battery life can reach 30 min. In the acquisition process, the scanning frequency was 550 KHz, the overlap rate of sidebands was 20 percent, the flight altitude was 110 m, the flight speed was 5 m/s, and we conducted four flight scanning missions in total.

3. Methodology

3.1. Downsampling with Edge-Preserving Information

The downsampling of the point cloud datasets can reduce the point cloud density, reduce the requirements for computing power, and improve the operation speed. The paper uses a point cloud downsampling framework based on an edge detection algorithm. Its technical roadmap is shown in Figure 3.

First, the two-dimensional images from multiple angles were processed to calculate the two-dimensional feature points, such as Holistically-Nested Edge Detection (HED), binarization, contour thinning, etc., and then the three-dimensional features were calculated by using the corresponding relationship between the two-dimensional image and the three-dimensional model. Finally, the non-feature points were sampled to complete the downsampling of the model. In brief, the point cloud downsampling network based on edge detection algorithm mainly uses two-dimensional images with several different angles.

3.2. The First Segmentation Based on Intensity Information

The first segmentation and first classification of targets using LiDAR intensity information need to analyze the influence of each factor on intensity. The intensity returns a value of 3D LiDAR and can be regarded as the discrete integer value of the received power of LiDAR after certain signal processing. The echo intensity of LiDAR is mainly characterized by the LiDAR equation [20,21,22]. According to the LiDAR equation, the laser echo light power is determined by the laser transmitting power, the transmitting optical system parameters, the atmospheric laser transmission characteristics, the receiving optical system parameters, and the target reflection characteristics [23,24]. The expression is as follows:

$$P_r=\frac{4P_t}{\pi {\theta }_t^2{R}_t^2}·{\tau }_t{\eta }_t{\tau }_r{\eta }_r·\rho cos{\theta }_{i}cos{\theta }_r{A}_t\frac{D^2}{4R_r^2}$$

(1)

where $P_r$ represents the received echo power of the target, $P_t$ represents the power of the laser emission, ${\theta }_t$ represents the laser reflection angle, $R_t$ is the distance from the laser light source to the target, ${\tau }_t$ represents the atmospheric transmittance between the laser light source and the target, ${\eta }_t$ represents the optical efficiency of the transmitting optical system, ${\tau }_r$ s the atmospheric transmittance between the target and the receiver, ${\eta }_r$ represents the optical efficiency of the receiving optical system, $\rho$ represents the reflectivity of the target, ${\theta }_i$ and ${\theta }_r$ are the incident and reflection angles of the target relative to the LiDAR, $A_t$ represents the effective area of the laser beam on the target, $R_r$ represents the distance from the receiver to the target, and D represents the aperture of the receiving detector.

The common three-dimensional scanning LiDAR system is generally a single station structure integrating transmitting and receiving [25,26,27,28,29]. In this case, there are the following equation relations: $R_t=R_r=R$, ${\theta }_i={\theta }_r=\theta$, and ${\tau }_t={\tau }_r=e^{-\alpha R}$, where $\alpha$ is the atmospheric extinction coefficient. In summary, the LiDAR equation can be derived in the following form:

$$P_r=P_t{\eta }_t{\eta }_r\frac{D^2}{4}·\frac{e^{-2\alpha R}}{R^2}·\rho co{s}^2\theta$$

(2)

Then, the downsampled point cloud datasets were attached with different colors according to different intensities. On this basis, according to the intensity distribution of point cloud, the statistical histogram of the intensity information of laser point cloud was drawn, the suitable segmentation range was selected, the starting and ending points of intensity segmentation value were set in the program, and the initial segmented point cloud was obtained.

3.3. The Second Segmentation with The Cloth Method

A cloth simulation filtering (CSF) algorithm is based on a simple physical process simulation [30,31,32,33]. The algorithm assumes that a virtual cloth falls by gravity on the terrain surface. If the cloth is soft enough, it will be attached to the terrain, and the shape of the fabric is the $DSM$. When the terrain is turned over, the fabric shape that falls on the surface is $DEM$, as shown in Figure 4a. The algorithm uses a cloth simulation to extract the ground points of laser point cloud data.

The cloth is a grid made up of a large number of interrelated nodes, as shown in Figure 4b, which is called the particle spring model [34]. In this model, the cloth points are connected by a virtual spring, and the interaction between the points follows the law of elasticity. In other words, the cloth point will produce displacement after being stressed. To simulate the shape of the cloth at a certain time, the position of the cloth points in the $3D$ space needs to be calculated. Combined with Newton’s second law, the relationship between the position of the distribution point and the interaction force follows the following formula:

$$X(t+\Delta t)=2X\left(t\right)-X(t-\Delta t)+\frac{G}{m}\Delta t^2$$

(3)

where m represents the mass of the distribution points and is set as constant 1; X represents the node position at a certain time; $\Delta t$ represents the time step; G represents the gravitational constant. Given the time step and the initial position of the node, the position of the current point will be obtained.

In the process of cloth simulation, the position of the cloth point after the displacement caused by gravity is calculated by Formula (3). In addition, the blank area of the overturned surface is usually micro-terrain such as buildings or deep pits. Therefore, to limit the movement of distribution points, it is necessary to calculate the elevation difference between adjacent distribution points to correct the position of distribution points moved by the force between adjacent nodes. If two adjacent nodes are movable points and have different elevation values, they move the same distance in the opposite direction in the vertical direction. If there is only one immovable point between them, the other point moves. If both are at the same elevation, neither is moved. The corrected displacement of each cloth point is calculated as follows:

$$d=\frac{1}{2}b(p_i-p_0)·n$$

(4)

where d is the movement vector of the node; $p_0$ is the current position of the node to be moved; $p_i$ is the adjacent node position of $p_0$; n is the standard vector in the vertical direction, $n={(0,0,1)}^T$; b is the parameter used to determine whether the node moves. When the node is a movable point, set $b=1$, otherwise $b=0$.

3.4. Ski Tracks Extraction

Before the semantic segmentation of the ski tracks and the completion of the extraction work, it is necessary to establish a topological relationship for the three-dimensional datasets of point cloud, determine the seed nodes in datasets to be segmented, and according to the growth criteria, judge whether the neighborhood points and seed nodes in the seed stage belong to the same surface. If the conditions are met, the point continues to grow as a new seed point, and so on. The datasets of ski-track point cloud with the same attributes are divided into the same region until the growth ceases. In the paper, the k-d tree method is used to establish the topological relationship between the ski-track datasets, because it has the advantages of fast neighborhood point search speed and high segmentation efficiency.

3.4.1. K-Neighborhood Acquisition of k-d Tree

The topological relationship between the scattered point cloud datasets of Wanlong ski resort must be established before point clouds can be grown and segmented. With the establishment of a topological relationship for the point cloud, the processing range of data can be effectively reduced, and the algorithm will perform more efficiently. Currently, octree and k-d trees are commonly used point cloud data organizing algorithms. In spite of the fact that the time complexity of searching K adjacent points by both of them is $O\left(n\right)$, the octree method does not provide a satisfactory location for a single point. Therefore, this paper uses k-d trees to establish the topological relationship of point clouds, and the FLANN library is used to implement k-d tree search in PCL. K-dimensional trees are used to organize the set of points in k-dimensional space. It is similar to a binary search tree but with additional constraints. The three-dimensional point cloud is analyzed using a three-dimensional k-d tree. The k-d tree places the value A of the specified dimension on the root, and then $B,D$, and E with smaller values are found on the left subtree, and the larger values $C,F$, and G are found on the right subtree. Once this process is repeated for both the left and right subtrees, only one element will remain in the last tree. The specific workflow is shown in Figure 5.

3.4.2. Establishment of Criteria for Regional Growth and Selection of Initial Seed Points

In region growing segmentation, the choice of seed points determines the accuracy of the segmentation. We will use the point with the smallest curvature as the seed point in this paper and grow from the point with the smallest curvature. Starting to grow at the point with the lowest curvature is the smoothest area, which can reduce the total number of areas and increase the efficiency. To determine the curvature of the point cloud in the first-segmented points, the following formula should be used.

(1) Calculate the average curvature $K_h$ of a point P in the surface. The average curvature $K_h$ satisfies the following relationship:

$$2K_h\overrightarrow{n}=lim\left[diam\left(A\right)\right]⟶0\frac{▿A}{A}$$

(5)

In Formula (5), $\overrightarrow{n}$ stands for the normal vector, A stands for an infinitesimal region around P, $diam\left(A\right)$ represents its diameter, and ∇ is the gradient operator about point P.

(2) Discretize Equation (1) to obtain the average curvature of $P_i$:

$$K_h\left(P_i\right)=\frac{n}{4A_{min}}\gamma$$

(6)

$$\gamma =\sum _{j\in N\left(i\right)}(cot{\alpha }_{ij}+cot{\beta }_{ij})(P_i-P_j)$$

(7)

In Formulas (6) and (7), ${\alpha }_{ij}$ and ${\beta }_{ij}$ are included angles between the line connecting point $P_i,P_j$ and their normal vectors, respectively.

This algorithm calculates the angle between the normal of the neighborhood point and that of the current seed point and adds the neighborhood points that are less than the smooth threshold of the current area. The smooth threshold is implemented in the $PCL$ library by calling the setSmoothnessThreshold member function and then checking the curvature of each neighborhood point and adding them to the current seed point sequence. To estimate the normal vector of the point cloud, the $PCL$ library uses the midpoint and the neighborhood fitting plane to determine the normal vector of the point. The problem can be rephrased into a least squares plane fitting estimation problem. In other words, estimating the surface normals is equivalent to analyzing the eigenvectors and eigenvalues of a covariance matrix, which is created from the nearest neighbors of the query point. Specifically, for each point $P_i$, the corresponding point covariance matrix C is as follows:

$$C=\frac{1}{k}\sum _{i=1}^k(P_i-\overline{P}){(P_i-\overline{P})}^T$$

(8)

$$C\left({\overrightarrow{n}}_j\right)={\lambda }_j{\overrightarrow{n}}_j$$

(9)

In Formulas (8) and (9), the number of adjacent points of $P_i$ is given by k, $\overline{P}$ is the three-dimensional centroid of the closest element, ${\lambda }_j$ is the jth eigenvalue of the covariance matrix, and ${\overrightarrow{n}}_j$ is the jth eigenvector. The eigenvector associated with the smallest eigenvalue is the normal vector of the plane, i.e., the normal of the point cloud. By calling the third Eigen library in $PCL$, it is possible to estimate the covariance matrix of a point set.

3.4.3. Setting the Threshold for Curvature

The setting of the curvature threshold in segmenting regions directly impacts the quality of segmentation. By selecting the curvature threshold incorrectly, one can cause oversegmentation or undersegmentation. The curvature threshold value is directly related to the complexity of the actual object; for example, the curvature threshold value of the gentler object is larger, whereas the curvature threshold value of the sharper object is smaller. The curvature threshold can be set in the $PCL$ library by using the setCurvatureThreshold member function. Based on the points data of the Wanlong ski resort used in this article, the curvature threshold is set between $0.0016$ and $0.0023$.

The growth radius and vertical distance between the neighboring points and the seed surface were judged by setting the range of spatial threshold. If the angle is less than the threshold, the current point is added to the current area. The normal vector of the fitting surface is calculated and compared with the normal vector of the seed surface to judge whether the two included angles belong to the threshold range. If the curvature is less than the critical value, the point is added to the seed point. In summary, the growth segmentation process is shown in Figure 6, and the point cloud extraction algorithm for ski tracks is in Algorithm 1.

Algorithm 1: The Point Cloud Extraction Algorithm for Ski Tracks
Require: Point cloud = {P}; Point normals = {N};
1:	Points curvatures ={K };
2:	Neighbor finding function = { $\Omega$$(.)$ };
3:	Curvature threshold = { $K_{th}$ };
4:	Angle threshold = { ${\theta }_{th}$ };
5:	Region list $R\leftarrow 0$;
6:	Available points list $\left\{A\right\}\leftarrow \{1,\dots ,|P\left|\right\}$.
Ensure: Clustering points set, M
7:	$\mathbf{while}\left\{A\right\}$ is not empty do
8:	Current region $\left\{R_K\right\}\leftarrow 0$
9:	Current seeds $\left\{S_K\right\}\leftarrow 0$
10:	Point with minimum curvature in $\left\{A\right\}\to 0$
11:	$\left\{S_K\right\}\leftarrow \left\{S_K\right\}\cup P_{min}$
12:	$\left\{R_K\right\}\leftarrow \left\{R_K\right\}\cup P_{min}$
13:	$\left\{A\right\}\leftarrow \left\{A\right\}\setminus P_{min}$
14:	for$i=0$ to size $<\left\{S_K\right\}>$ do
15:	Find nearest neighbours of current seed point $\left\{B_K\right\}\to \Omega \left(S_K\left\{i\right\}\right)$
16:	Current neighbour point $\left\{P_j\right\}\to B_K\left\{j\right\}$
17:	if $\left\{A\right\}$ contains $P_j$ and $co{s}^{-1}\left(\right|\left(N\{S_K\left\{i\right\},N\left\{S_K\left\{j\right\}\right\})\right|)<{\theta }_{th}$ then
18:	$\left\{R_K\right\}\leftarrow \left\{R_K\right\}\cup P_j$
19:	$\left\{A\right\}\leftarrow P_j$
20:	if $K\left\{P_j\right\}<K_{th}$ then
21:	$\left\{S_K\right\}\leftarrow \left\{S_K\right\}\cup P_j$
22:	end if
23:	end if
24:	end for
25:	Add current region to global segment list $\left\{R\right\}\leftarrow \left\{R\right\}\cup \left\{R_K\right\}$
26:	$\mathbf{end}\mathbf{while}$
27:	return M

4. Results and Discussion

In this paper, a total of four samples from two regions were selected for the downsampling experiment. The downsampling method of the point cloud datasets of the Wanlong ski resort based on edge detection algorithm can obtain the following four groups of results:

Figure 7 shows the downsampling results of the first laser point cloud datasets in region A; the number of raw data is 8,283,717 points, and the number after downsampling is 398,776 points. Figure 8 shows the downsampling results of the second laser point cloud datasets in region A; the number of raw data is 5,318,641 points, and the number after downsampling is 5,318,641 points.

Figure 9 shows the downsampling results of the first laser point cloud datasets in region B; the number of raw data is 2,648,473 points, and the number after downsampling is 293,198 points. Figure 10 shows the downsampling results of the second laser point cloud datasets in region B; the number of raw data is 497,563 points, and the number after downsampling is 242,744 points.

After downsampling, the edges of trees and ski tracks in point cloud datasets are obvious, and the information is clearly visible. Through analyzing the above four sets of results, we found that this method can effectively preserve the edge information of snow tracks, gentle slopes, and other places. The successful downsampling of laser point cloud datasets lays a foundation for the subsequent segmentation and classification of ground and ground objects, and the extraction of ski tracks.

Next, the target was preliminarily segmented and classified according to the intensity information of laser point cloud datasets. First, the points were attached with different $RGB$ values according to the echo intensity to display the segmentation and classification effect directly. Then, the statistical histogram of the laser point cloud was drawn according to the echo intensity. Finally, the point cloud datasets wase segmented according to the intensity information distribution, and the segmentation results are as follows:

After analyzing the intensity statistical histogram, in Figure 11 and Figure 12, it was determined that the intensity segmentation range of the second laser point cloud datasets in region A is 45,688 to 53,588. The result of point cloud sets segmentation is shown in Figure 13, where (a) is the non-ground points and (b) is the ground points.

The intensity statistical histogram based on the Figure 14 and Figure 15 indicates that the intensity segmentation range of the second laser point cloud datasets in region A is 47,587 to 53,488. The result of point cloud sets segmentation is shown in Figure 16, where (a) is the non-ground points and (b) is the ground points.

After analyzing the intensity statistical histogram, in Figure 17 and Figure 18, it is determined that the intensity segmentation range of the first laser point cloud datasets in region B is 46,073 to 52,903. The result of point cloud sets segmentation is shown in Figure 19, where (a) is the non-ground points and (b) is the ground points.

After analyzing the intensity statistical histogram, in the Figure 20 and Figure 21, it is determined that the intensity segmentation range of the second laser point cloud datasets in region B is 46,032 to 52,918. The result of point cloud sets segmentation is shown in Figure 22, where (a) is the non-ground points and (b) is the ground points.

Analyzing the above results of laser point cloud segmentation based on the intensity distribution of LiDAR point cloud can easily distinguish the snow-covered ground from the non-ground points, but it is still unable to distinguish the ski tracks from the snow-covered ground, so further segmentation operation is needed.

Therefore, to further separate the surface, the user-defined parameters of the cloth algorithm were set separately. After several adjustments, the best results can be obtained when cloth resolution is set to $0.2$, max iterations is set to 500, and the classification threshold is set to $0.5$. The point cloud segmentation results are shown in the Figure 23, Figure 24, Figure 25 and Figure 26.

The included angle between the normal of the seed point field and the normal of the current seed point is calculated and compared with the smoothing threshold. After several attempts, when the smoothing threshold is set to $2.4$ and the curvature is set to $0.2$, the ski tracks can be separated from other ground points successfully. The results are shown in the Figure 27, Figure 28, Figure 29 and Figure 30.

For further verification of the effectiveness, this paper compares the method presented in this paper with the Euclidean distance method, the cluster segmentation method, and the RANSAC method on three aspects: correct rate, oversegmentation rate, and undersegmentation rate. From the above four groups of test data, snow tracks were extracted using different methods, and the specific comparison results are shown in

Table 1. Comparison of ski-track extraction effects with different extraction algorithms.

Name of Points	Extraction Algorithm	Correct Rate (%)	Over-Extraction Rate (%)	Under-Extraction Rate (%)
First Data in Region A	Euclidean Distance	68.3	20.4	11.3
	Clustering Segmentation	76.9	6.2	16.9
	RANSAC Method	82.6	7.3	10.1
	Proposed Method	92.2	4.7	3.1
Second Data in Region A	Euclidean Distance	75.6	10.4	14.0
	Clustering Segmentation	74.3	12.6	13.1
	RANSAC Method	85.4	8.9	5.7
	Proposed Method	94.5	2.8	2.7
First Data in Region B	Euclidean Distance	59.8	26.3	13.9
	Clustering Segmentation	72.5	14.3	13.2
	RANSAC Method	80.1	10.7	9.2
	Proposed Method	89.4	6.8	3.8
Second Data in Region B	Euclidean Distance	72.2	15.4	12.4
	Clustering Segmentation	70.8	18.5	10.7
	RANSAC Method	85.6	7.8	6.6
	Proposed Method	92.9	5.1	2.0

.

Taking into account the analysis of Table 1, the extraction accuracy rate of the method proposed in this paper is the highest (92.2%), while the under-extraction rate is the lowest (3.1%), whereas the correct extraction rate of the European distance method is the lowest (68.3%), and the over-extraction rate is the highest (20.4%). As for the second group of experimental data, the extraction accuracy of the proposed method is highest (94.5%), and the under-extraction rate is lowest (2.7%), whereas the clustering segmentation method has the lowest correct extraction rate (74.3%), and the over-extraction rate is highest (12.6%). With the third group of experimental data, the extraction accuracy of most methods has decreased due to the complex test data environment; however, the extraction accuracy of the methods presented in this work is still the highest (89.4%), and the under extraction rate is the lowest (3.8%), whereas the European clustering method has the lowest correct extraction rate (59.8%) and the highest over extraction rate (26.3%). As evidenced by the fourth set of experimental data, the proposed method had the highest extraction accuracy (92.9%) and the lowest under-extraction rate (2.0%), while the clustering segmentation method had the lowest extraction accuracy (70.8%) and the highest over-extraction rate (18.5%).

Analyzing the results in Table 1, it can be seen that the method proposed in this article achieves the highest ski-track point cloud extraction accuracy, the lowest over-extraction rate, and the lowest under-extraction rate when compared with the European distance method, the clustering segmentation method, and the RANSAC method. Accordingly, compared with other methods, the average extraction accuracy of the ski-track point cloud is improved by 16.9%, the over-extraction rate is reduced by 8.4%, and the under-extraction rate is reduced by 8.6%, which enables effective extraction of the ski track points cloud of a ski resort.

5. Conclusions and Future Work

Compared with the current mainstream classification algorithms, this paper used a composite detection method based on the laser intensity information reflected by the target. The method in this paper could complete the downsampling of point cloud datasets under the condition of ensuring the integrity of target edge information and separating the ground points and non-ground points according to the intensity distribution histogram and cloth method. Finally, the clustering and extraction of ski tracks were effectively completed based on the smoothing threshold and curvature between adjacent point clouds. The future work will focus on the study of laser point cloud intensity distribution and simplify the extraction process. In addition, the depth and slope of the ski tracks will be calculated according to the snow-free and snow-covered scenes, respectively.