Displacement Prediction of the Muyubao Landslide Based on a GPS Time-Series Analysis and Temporal Convolutional Network Model

Landslide displacement prediction is an essential base of landslide hazard prevention, which often needs to establish an accurate prediction model. To achieve accuracy prediction of landslide displacement, a displacement prediction model based on a salp-swarm-algorithm-optimized temporal convolutional network (SSA-TCN) is proposed. The TCN model, consisting of a causal dilation convolution layer residual block, can flexibly increase the receptive fields and capture the global information in a deeper layer. SSA can solve the hyperparameter problem well for TCN model. The Muyubao landslide displacement collected from a professional GPS monitoring system implemented in 2006 is used to analyze the displacement features of the slope and evaluate the performance of the SSA-TCN model. The cumulative displacement time series is decomposed into trend displacement (linear part) and periodic displacement (nonlinear part) by the variational modal decomposition (VMD) method. Then, a polynomial function is used to predict the trend displacement, and the SSA-TCN model is used to predict the periodic displacement of the landslide based on considering the response relationship between periodic displacement, rainfall, and reservoir water. This research also compares the proposed approach results with the other popular machine learning and deep learning models. The results demonstrate that the proposed hybrid model is superior to and more effective and accurate than the others at predicting the landslide displacement.

Introduction
Landslides are common natural geological disasters that cause a lot of property loss and personal injury every year [1,2]. The landslide early warning system (LEWS), including real-time monitoring methods such as remote sensing (RS), the global position system (GPS), and the geographic information system (GIS), is an effective approach to identify the early displacement and development trend of landslides [3][4][5]. The landslide surface displacement prediction based on real-time monitoring data is an important part for LEWS [6,7]. It can provide advanced landslide displacement information to decision makers, which reduces the damage caused by landslides. Therefore, the accurate prediction of landslide displacement can be an effective method to reduce landslide risks.
Numerous approaches have been developed for the purpose of landslide displacement prediction. They can be divided into two main types: physically-based approaches and data-driven approaches [8][9][10]. Among these models, the physically-based models mainly take geotechnical characteristics into analysis and quantify the slope displacement by Convolutional neural networks are widely applied in extracting features from multiple dimensions data with Euclidean space [33,34]. A typical CNN model consists of a convolution layer, an activation layer, a pooling layer, and a fully connected layer. On the convolution layer, a set of trainable filters slide over the input data to perform a convolution operation and extract features; they are called feature maps. The individual units in feature maps correspond to the local positions of original input that the filter currently covers [35] ( Figure 1). In addition, at every local region of input data, the convolution filter uses the same weights. This is the weight-sharing strategy, which can reduce the complexity of the convolution layer and the number of parameters that are required to be trained. The pooling layers usually follow the convolution layer, and adjacent information in the feature maps is merged by max pooling, average pooling, or other operations. Finally, the fully connected layer is used to transform the features processed by convolution and pooling layers into one-dimensional vectors for conveniently processing and outputting results. The activation layer puts nonlinear functions into the neural network, which can improve the nonlinear expression ability of the model. The commonly used activation equations are the Sigmiod function, the ReLu function, and the Tanh function, etc. Actually, according to different tasks and dimensions of input, data that needs to be processed by CNN can be divided into three different dimensions network structures. The one-dimensional CNN (1D-CNN) model is used for sequential data and is applicable to signal processing problems. The time-series data were processed as different channels data that not only can build nonlinear relationships between different input features but also merge nearby time-step data by the convolution layer and pooling layer in the 1D-CNN model. So, the 1D-CNN model can extract temporal features from original sequence data. The two-dimensional CNN (2D-CNN) model is mainly used for image and text recognition, and the three-dimensional CNN (3D-CNN) model is used for video data. Thus, in this study, the 1D-CNN model was also applied in landslide displacement prediction for competition with the TCN model. The main contributions of the present study are summarized as follows: (i) a hybrid model based on TCN model was proposed for displacement prediction. (ii) The linear and the nonlinear parts of the total displacement were extracted by a time-series decomposition method that can reduce the complexity of the time series. (iii) A large-scale landslide with step-like displacement behavior in the TGRA of China was used to validate the feasibility and effectiveness of the proposed model.

Convolutional Neural Network (CNN)
Convolutional neural networks are widely applied in extracting features from multiple dimensions data with Euclidean space [33,34]. A typical CNN model consists of a convolution layer, an activation layer, a pooling layer, and a fully connected layer. On the convolution layer, a set of trainable filters slide over the input data to perform a convolution operation and extract features; they are called feature maps. The individual units in feature maps correspond to the local positions of original input that the filter currently covers [35] (Figure 1). In addition, at every local region of input data, the convolution filter uses the same weights. This is the weight-sharing strategy, which can reduce the complexity of the convolution layer and the number of parameters that are required to be trained. The pooling layers usually follow the convolution layer, and adjacent information in the feature maps is merged by max pooling, average pooling, or other operations. Finally, the fully connected layer is used to transform the features processed by convolution and pooling layers into one-dimensional vectors for conveniently processing and outputting results. The activation layer puts nonlinear functions into the neural network, which can improve the nonlinear expression ability of the model. The commonly used activation equations are the Sigmiod function, the ReLu function, and the Tanh function, etc. Actually, according to different tasks and dimensions of input, data that needs to be processed by CNN can be divided into three different dimensions network structures. The one-dimensional CNN (1D-CNN) model is used for sequential data and is applicable to signal processing problems. The time-series data were processed as different channels data that not only can build nonlinear relationships between different input features but also merge nearby time-step data by the convolution layer and pooling layer in the 1D-CNN model. So, the 1D-CNN model can extract temporal features from original sequence data. The two-dimensional CNN (2D-CNN) model is mainly used for image and text recognition, and the three-dimensional CNN (3D-CNN) model is used for video data. Thus, in this study, the 1D-CNN model was also applied in landslide displacement prediction for competition with the TCN model.

Temporal Convolutional Network (TCN)
The TCN module is a new and unique CNN model that combines the advantages of the time-series modeling capabilities of RNN and CNN's parallel processing. In particular, the core components of TCN contain two causal dilated convolution layers and a residual block. The causal convolution layer stresses the order for sequence processing, and the dilated convolution layer can flexibly adjust the size of the receptive filed in TCN. The residual block is an effective method that can avoid the degradation of the deep network [36]. Causal convolutions are convolutions where the output at time t can only depend on the inputs that are no later than time t, as shown in Figure 2a. However, the insufficient receptive of the causal convolution layer has a limited ability to capture global information. In order to capture much longer sequence information, the dilation strategy through interval sampling is applied to obtain larger perception. For a given input 1-D dimension sequence x, the output at time t of the dilation convolution layer can be defined as: where the s is the output and f is a 1-D dimension filter with size k. * is the convolution operation, d is the dilation factor, and t − d · i indicates the cells covered in the perspective field. The dilation factor d and the filter size k in the TCN can flexibly adjust the size of the receptive filed to obtain features of the input sequence with different time scales. Figure 2b shows an example of the dilation strategy, where the filter size is two and the interval sampling period d = 1, 2.

Temporal Convolutional Network (TCN)
The TCN module is a new and unique CNN model that combines the advantages of the time-series modeling capabilities of RNN and CNN's parallel processing. In particular, the core components of TCN contain two causal dilated convolution layers and a residual block. The causal convolution layer stresses the order for sequence processing, and the dilated convolution layer can flexibly adjust the size of the receptive filed in TCN. The residual block is an effective method that can avoid the degradation of the deep network [36].

Causal Dilation Convolution Layer
Causal convolutions are convolutions where the output at time t can only depend on the inputs that are no later than time t, as shown in Figure 2a. However, the insufficient receptive of the causal convolution layer has a limited ability to capture global information. In order to capture much longer sequence information, the dilation strategy through interval sampling is applied to obtain larger perception. For a given input 1-D dimension sequence x, the output at time t of the dilation convolution layer can be defined as: where the s is the output and f is a 1-D dimension filter with size k. * is the convolution operation, d is the dilation factor, and t − d·i indicates the cells covered in the perspective field. The dilation factor d and the filter size k in the TCN can flexibly adjust the size of the receptive filed to obtain features of the input sequence with different time scales. Figure  2b shows an example of the dilation strategy, where the filter size is two and the interval sampling period d = 1, 2.

Residual Block
The residual block can be expressed as that the original input of the TCN module is added to the output of the causal dilation convolution layer, and then the addition is output by an activation function (Figure 2c). The residual block can be calculated as where x and o are the input and output of the TCN module; F represents a series of operations, including causal dilated convolution, and activation function. In order to make the

Residual Block
The residual block can be expressed as that the original input of the TCN module is added to the output of the causal dilation convolution layer, and then the addition is output by an activation function (Figure 2c). The residual block can be calculated as where x and o are the input and output of the TCN module; F represents a series of operations, including causal dilated convolution, and activation function. In order to make the input x with the same size of F(x), the 1 × 1 convolution layer is added to express the input. So, the Equation (2) can summarized as follows:

Total Displacement Decomposition Based on Variational Mode
Previous works [20] indicated that the cumulative displacement of the landslide was affected by the internal geological conditions (lithology, structural conditions, topography, etc.) and external environment factors (rainfall, reservoir water level). The step-like displacement of landslide is a non-stationary time-series varying with time, which can be divided into trend displacement and periodic displacement. Geological conditions are the main factor affecting the trend displacement, and the curve of trend displacement shows an approximate monotonous increasing trend for a long period. The displacement fluctuation in a short period is affected by external environmental factors. Therefore, the accumulated displacement can be decomposed as follows: where t is the time, x t is the accumulated displacement, s t is the trend displacement, and p t is the periodic displacement. Variational mode decomposition (VMD) is a signal decomposition technology using a non-recursive decomposition mode [37]. Based on convergence conditions, it can control the number of modes obtained by original data. Assuming that each mode has a limited bandwidth and is compacted around their respective center frequencies, the goal is to minimize the sum of the bandwidth of each mode. The constraint is that the sum of the mode functions v k is equal to the original signal. The mathematical description of the variational problem is as follows: where e −jω k t is the vector of the center frequency on the complex plane, v k represents the mode function, ω k is the center frequency of each component, and f is the original data.
In order to solve the constrained variational problem, a quadratic penalty term and lagrange multiplier are used to transform the problem unconstrained. The augmented lagrangian, L, is depicted as follows: where α is the penalty factor and λ is the lagrange multiplier.
In the process of solving the constrained variational problems, the alternating direction multipliers method (ADMM) is used to constantly update v k , ω k, and λ for finding the saddle point of the lagrange expression, that is, the solution of the above equation. The adaptive decomposition of the signal is completed according to the frequency domain characteristics of the original signal [38].

Salp Swarm Algorithm
Salp swarm algorithm (SSA) [32] is a global optimization algorithm based on swarm biological intelligence proposed in recent years. The SSA simulates the movement and foraging behavior of salp chain in nature [39]. The salps are divided into two subgroups equally, named leaders and followers, respectively. Assuming that there is a food source obtained is also constantly updated. For the leaders, the movement method is carried out according to the following formula: where x j i shows the position of the i-th salp (leader) in the j-th dimension, i ≤ m/2; m is the number of salps; F j is the position of the food source in the j-th dimension; and ub j and lb j are the upper and lower bounds in the same dimension, respectively. c 2 and c 3 are generated randomly in range (0, 1), and c 1 is calculated as follows: where t is the current iteration step and t max is the max number of iteration step. c 1 is decreased during the iteration process for strengthening the different capabilities of the algorithm in different search stages.
The position of the followers is updated according to the following equation: where x j i shows the position of i-th follower in j-th dimension, i > m 2 .

Hybrid TCN Forecasting Model and Implementation Procedure
In this study, we present a hybrid model that uses the swarm intelligence algorithm as a parameter optimizer to automatically determine the hyperparameters for the TCN model. Compared with the machine learning model, the complicated model of deep learning requires extensive trial and error when determining the structure of the model. The SSA has an efficient search strategy and strong spatial search ability, which are likely to provide greater possibilities to obtain global optimal solutions. Any application of the swarm intelligence algorithm should meet two basic requirements: (i) a defined function, which can converse input to output; (ii) an evaluation indexed for an individual solution [40]. The optimizing deep learning model with the SSA meets the above requirements in the following two respects: (1) The deep learning model is often regarded as a 'black box' model. Although it is difficult to give specific mathematical expressions, the model is still a one-to-one correspondence between input and output when the weights and biases are determined. In this study, each individual in the SSA will be a potential solution for setting the hyperparameters of the TCN model. An individual in SSA can be expressed as a real vector S. For example, in TCN, the hyperparameters such as the kernel size, the number of kernels, and the number of neural units should be contained in S.
(2) The evaluation scheme is used to calculate the fitness value of each potential solution to determine the performance of the solution. In this study, the root mean squared error (RMSE) is chosen to evaluate the solution. Moreover, the mean absolute percentage error (MAPE) and correlation coefficient (R) are also adopted to verify the performance of the predictive model. These error criteria methods are defined as: Remote Sens. 2022, 14, 2656 where s i andŝ i represent the measured and predicted values; s andŝ represent the mean of observed and predicted values; and n is the number of predicted values. The implementation procedure of the hybrid model based on TCN for "step-like" displacement prediction is summarized as follows.
Step 1: Randomly generate m solutions in the k-dimensional search space.
Step 2: Introduce the solution S as hyperparameters of TCN model and train the TCN model on the training dataset.
Step 3: Evaluate the prediction results by calculating the fitness value of term RMSE on the validation dataset and sort all the solutions based on the fitness value. The solution corresponding to the optimal fitness value is assigned to F.
Step 4: Update the position of the leading and follower salps separately by the update formulas in Section 2.4.
Step 6: Return the TCN predictive model. Based on the time-series decomposition method mentioned above, we proposed a temporal convolutional network (TCN) optimized by SSA as shown in Figure 3. Firstly, the total displacement of the landslide is decomposed into two displacement sub-sequences named periodic term and trend term. The trend term is fitted and predicted by a polynomial function. Through the analysis of the monitoring data, the external candidate factors that trigger the periodic displacement are preliminarily selected, and the maximal information coefficient (MIC) [41] method is used to quantitatively determine the correlation value between the candidate factors and periodic displacement. The SSA is applied to optimize the hyperparameters of the TCN model. Then, the trained model is used to predict the periodic component displacement. The predicted result of the cumulative displacement is the sum of the predicted periodic displacement and trend displacement.

Topography and Geology of the Muyubao Landslide
The Muyubao landslide on the slope of the south bank of the Yangtze River is located in the county of Zigui country of Hubei province, 56 km away from the Three Gorges Dam (Figure 4a). The average annual rainfall in this study region is 1493.2 mm. The rainfall is concentrated from July to September, counting 70% of the cumulative rainfall each year.

Topography and Geology of the Muyubao Landslide
The Muyubao landslide on the slope of the south bank of the Yangtze River is located in the county of Zigui country of Hubei province, 56 km away from the Three Gorges Dam (Figure 4a). The average annual rainfall in this study region is 1493.2 mm. The rainfall is concentrated from July to September, counting 70% of the cumulative rainfall each year. The Muyubao landslide is a large-scale unstable slop that includes a volume of 96 million m 3 sliding mass, an area of 1.8 km 2 . The landslide has a maximum north-to-south length of about 1500 m and an average width of about 1200 m on the different direction (Figure 4b).
The upper boundary of the landslide is at an elevation of 540 m above sea level (a.s.l.), and the toe of the landslide dips into the Yangtze River at 100 m (a.s.l.). The upper portion of the landslide shows a steep slope gradient of nearly 25 • and gently lowers with the elevation decrease. The sliding mass is mainly composed of a surficial quaternary deposit and fractured quartz sandstone, as shown in Figure 4c. The surficial deposit is composed of silty clay and rubble, which makes up 40-60% of the total volume. The thickness of sliding mass ranges from 60 m to 120 m, and the front part is an uplift platform formed by shear sliding with a thickness of 80~120 m. The main sliding direction of the landslide is 20 • , which faces the Yangtze River. The bedrock mainly consists of Jurassic siltstone and Triassic quartz sandstone with dip directions of 27 • and dip angles of 25 • . The sliding zone is a thickness of 0.1 to 0.3 m layer of dark-gray to black coal and shale that extends along the interface between the bedrock and the sliding mass [42,43].

Displacement Characteristics Analysis
The landslide displays slow creep character before the initial TGR impoundment (June 2003), and activity signs including the ground fracture and infrastructure damage occurred rarely. However, by the filling-drawdown cycle of the reservoir level after 2003, the landslide deformation is gradually obvious, especially after the first trial impound-(c)

Displacement Characteristics Analysis
The landslide displays slow creep character before the initial TGR impoundment (June 2003), and activity signs including the ground fracture and infrastructure damage occurred rarely. However, by the filling-drawdown cycle of the reservoir level after 2003, the landslide deformation is gradually obvious, especially after the first trial impoundment to 156 m a.s.l. (November 2006). Since then, a total of 12 GPS points is set up on the Muyubao landslide in a regular form of three rows and four columns for surface displacement monitoring. The distribution of monitoring points is depicted in Figure 4b, and the other two reference stations (ZG12 and ZG13) are located on stable rock, which is out of the sliding area. The GPS is a kind of GNSS that can provide millimeter accuracy for landslide displacement monitoring [44]. The long-term GPS displacement monitoring results in Figure 5 indicated that most of the area of Muyubao exists at a similar displacement rate (except for ZG01, which is affected by the topography conditions). The displacement around section I-I' is the most significant, with a total displacement of 2447.9 mm (ZG01), 1693.0 mm (ZG02), 1678.9 mm (ZG03), and 1756 mm (ZG04). The displacement recorded by ZG05, ZG06, ZG07, and ZG08 around the central section of II-II' of the landslide are similar at cumulative displacement between 1560 mm with 1630 mm and an average velocity of about 10.96 mm/month. Compared with the other points, the displacement data measured by the four GPS station sites around profile III-III' on the west of the landslide show the lowest displacement, with an average cumulative displacement of 1521 mm. During the whole 12-year monitoring period, the horizontal displacement of the Muyubao landslide displays a linear creep character, but several significant displacement events were observed in the overall monitoring period. There are four evident acceleration events (2007, 2008, 2012, and 2017); each rapidly increased event was usually during a periodic of approximately 5 months (from November to April of the following year). During the four periods, the displacement increased approximately 30% of the total displacement observed in the December 2018. Moreover, the landslide displacement rapidly increased corresponding to the periods of relatively higher reservoir water level every year, while during the rest of the year, the landslide displacement experienced a slower increase, which corresponds to the relatively lower reservoir level period.

Displacement Decomposition Result
Through a detailed observation of the displacement time series, several displacement curves can be found with similar displacement trend such as ZG02 and ZG03, ZG05 and ZG06, ZG07 and ZG08, ZG10 and ZG11, and ZG04 and ZG12. Hence, the other two points

Displacement Decomposition Result
Through a detailed observation of the displacement time series, several displacement curves can be found with similar displacement trend such as ZG02 and ZG03, ZG05 and ZG06, ZG07 and ZG08, ZG10 and ZG11, and ZG04 and ZG12. Hence, the other two points (ZG01, ZG09) and the five typical points (ZG03, ZG06, ZG08, ZG11, and ZG12) from each pair of GPS points with similar displacement were selected to evaluate the performance of the proposed method. These typical points were located around the three geological profiles (I-I', II-II' and III-III') from the trailing edge to the former edge which can respect to the displacement status of the whole slope. For each GPS point, a total of 146 displacement samples were obtained by monthly recording from November 2006 to December 2018. In order to clearly provide the results of the proposed method in this study, in the following parts, the calculation results visualized in figures were based on the three points (ZG01, ZG06, and ZG12), and all seven points of the calculation results were specifically presented in tables. Figure 6 shows the trend and periodic displacement time-series that were extracted by VMD. The curves with approximately linear growth correspond to the trend displacements. The curves with high fluctuation frequency correspond to the periodic displacement components and show a similar one-year cycle of displacement characterized by acceleration from October to March or April of the following year and descending during April to September. During the four acceleration events (2007, 2008, 2012, and 2017), the average increments of periodic displacements were 50.72 mm, 72.32 mm, 81.27 mm, and 71.27 mm, and the average monthly trend displacements were 24.63 mm, 19.96 mm, 14.87 mm, and 15.13 mm, respectively. In the rest of these years, the average monthly trend displacements were 15.58 mm, 19.08 mm, 16.49 mm, and 12.53 mm, respectively. The trend displacement increase rates were similar even during the period with a displacement rapid increase, so the trend displacements can show an approximate line increase.

Data Preparation
The monitoring data collected from November 2006 to December 2017 of the selected seven points are chosen as the training sample, and the test samples contain the data from January to December 2018 test samples. The standardization function is used to standardize the training sample and test sample, respectively. The standardization function is calculated as where x is the standardization result, x is the original value, and u and σ are the sample mean and sample standard deviation of the selected dataset.

Analysis the Influencing Factors of Periodic Displacement
As the input of the periodic displacement prediction model, the trigger factor plays a key role in influencing the output of the final result. Hence, it is necessary to analyze the relationship between the influence factors and the periodic displacement. According to the aforementioned analysis, the landslide displacement shows a certain correlation with the rainfall and reservoir water level. The effect of reservoir operations on slope stability can be divided into a seepage effect or buoyancy effect [42]. The former can be regarded as the dynamic water seepage effect caused by the imbalance of the water levels in the changing process of the reservoir water level. The latter can be considered as the buoyancy effect by submerging different volumes of the slope to reduce effective stress. So, the average reservoir level in the previous month (a1) and the change of reservoir level during the last 1 month (a2) are considered as influencing factors on the reservoir water. The rainfall effect contributes to the destabilization of the soil by replenishing the groundwater and reducing effective stresses. The rainfall effect can be divided into "lateral flow" and "vertical flow" [45]. The lateral flow is considered as affecting the stability of the slope in the mid to long term timescale, and the vertical flow works in the short term [46]. Therefore, the candidate factors affecting rainfall on different timescales include the monthly rainfall (b1) and the cumulative rainfall of the last 2 months (b2). On the other hand, the different evolution states of the landslide make the displacement of the landslide respond differently to the external triggering factors [47,48]. The evolution status of the landslide can be obtained through the history movement information of the landslide. Moreover, the valuable information may be mined from the extracted displacements and velocities to predict the periodic displacement. Therefore, the velocity of cumulative displacement (c1) and the velocity of the periodic displacement (c2) are also selected as candidate factors.

Data Preparation
The monitoring data collected from November 2006 to December 2017 of the selected seven points are chosen as the training sample, and the test samples contain the data from January to December 2018 test samples. The standardization function is used to standardize the training sample and test sample, respectively. The standardization function is calculated as The displacement of slope is usually considered as a complex multi-dimension nonlinear system at spatiotemporal scale [49]. The complexity and uncertainty of the geological conditions and the randomness and suddenness of external factors are the two main aspects that impart the nonlinear characteristics to the displacement system. Correlation coefficients, such as Pearson's correlation coefficient, are usually used as feature selection algorithms due to their simple but practical character. However, these methods cannot usually provide a nonlinear relationship between inducing factors and displacement [50]. The maximal information coefficient (MIC) based on the statistical method is used for capturing the correlation between two variables. MIC not only calculates linear and non-linear correlations but also can identify important and difficult-to-recognize correlations [41,51].
Hence, the MIC was applied to identify and screen the environmental factors that induce landslide displacement in this research. The MIC value is between 0 and 1. If two variables show strong dependence, the MIC value is close to 1. Otherwise, it will tend to 0. Figure 7 presents linear and nonlinear relationships (MIC values) among the selected factors and the displacement at selected points of the Muyubao landslide. The MIC values of all candidate factors are equal to or greater than 0.2, revealing the strong correlation relationship between candidate factors and periodic displacement. So, the six candidate factors are selected as the input factors to predict the periodic displacements. In the following, the six factors data were constructed into a six-channel data, and each channel was used to input one factor time-series data, which is similar to the input data shown in Figure 1a, which can be processed by the TCN and 1D-CNN models.

Hyper-Parameters and Models Training
In this research, a series experiments were coded by python 3.6 and implemented on a personal computer (CPU: AMD Ryzen 5 2400G 3.6 GHz, 8 GB), including the proposed model and the other comparison models in the later section. In the TCN module, the receptive fields were directly affected by the dilation factors d and the kernel size (k). Insufficient reception hinders the extraction of long-term history information, while the excessively large receptive files increase the calculation burden. Hence, the SSA optimizer were used to determine the hyper-parameters and the model training parameters of the TCN model. In this study, there are six initial input channels of TCN, which corresponds to the number of input factors. Moreover, three causal dilation layers were constructed for the TCN module, and the dilation factors were set to {1, 2, 4}. The kernel size (K), the number of filters (Nf), and the initial learning rate (lr) were optimized by SSA. The specific search space of the hyperparameters was defined as a space R 7 ∈ [1, 6] 3 × [1, 64] 3 × [0.0001, 0.05]. For the SSA optimizer, the population size and iteration epochs are 20 and 200, respectively. The loss function was RMSE, and the iteration process was stopped when the loss of validation dataset was minimized. In comparison, the other two models that has performed well in the landslide displacement prediction (i.e., SVM, LSTM) [22,52] and the 1D-CNN model were applied in competition with the proposed model. The SSA was used to determine the hyperparameters of the three comparison models, and the specific parameters are shown in the Table 1. For the LSTM model, parameter Nu represents the LSTM unit number, and Nl represents the number of layers of the LSTM. In the SVM model, the kernel function was set to "rbf", with the parameters C and γ representing the regularization parameter and the kernel parameter, respectively. In the 1D-CNN model, the hyperparameters Nk, Kc represent the number of kernels and the kernel size in the convolution layer, respectively.

Hyper-Parameters and Models Training
In this research, a series experiments were coded by python 3.6 and implemented on a personal computer (CPU: AMD Ryzen 5 2400G 3.6 GHz, 8 GB), including the proposed model and the other comparison models in the later section. In the TCN module, the receptive fields were directly affected by the dilation factors d and the kernel size (k). Insufficient reception hinders the extraction of long-term history information, while the excessively large receptive files increase the calculation burden. Hence, the SSA optimizer were used to determine the hyper-parameters and the model training parameters of the TCN model. In this study, there are six initial input channels of TCN, which corresponds to the number of input factors. Moreover, three causal dilation layers were constructed for the TCN module, and the dilation factors were set to {1, 2, 4}. The kernel size (K), the number of filters (N f ), and the initial learning rate (lr) were optimized by SSA. The specific search space of the hyperparameters was defined as a space R 7 ∈ [1, 6] 3 × [1, 64] 3 × [0.0001, 0.05]. For the SSA optimizer, the population size and iteration epochs are 20 and 200, respectively. The loss function was RMSE, and the iteration process was stopped when the loss of validation dataset was minimized. In comparison, the other two models that has performed well in the landslide displacement prediction (i.e., SVM, LSTM) [22,52] and the 1D-CNN model were applied in competition with the proposed model. The SSA was used to determine the hyperparameters of the three comparison models, and the specific parameters are shown in the Table 1. For the LSTM model, parameter N u represents the LSTM unit number, and N l represents the number of layers of the LSTM. In the SVM model, the kernel function was set to "rbf", with the parameters C and γ representing the regularization parameter and the kernel parameter, respectively. In the 1D-CNN model, the hyperparameters N k , K c represent the number of kernels and the kernel size in the convolution layer, respectively. Table 1. The hyperparameters of TCN and the three comparison models.

Model
Hyperparameter Set

Trend Displacement Prediction
The trend displacement is a smooth monotone increasing sequence (Figure 6a). A polynomial with five orders is used to fit the trend displacement in the training sample. The least squares method is applied to find the best fitting values for the parameters of the polynomial function. The polynomial function is as follows: where y t is the fitted values of the trend term displacement at the time t, and a, b, c, d, e, and f are the coefficients, where a cannot be zero. The fitting polynomial of degree five and the prediction results of the trend displacements are shown in Figure 8 and Table 2. As exhibited in Figure 8, all the polynomials of degree five can accurately fit the increasing trend of landslide displacements and predict the following development in the test datasets. The R values of ZG01, ZG03, ZG06, ZG08, ZG09, ZG11, and ZG12 on the test samples are 0.991, 0.991, 0.982, 0.989, 0.989, 0.984, and 0.990, respectively. The fitting results of the trend displacements are close to the monitored displacement, mainly because the VMD can reduce the complexity of the original data and obtain relatively stable displacement components.

Trend Displacement Prediction
The trend displacement is a smooth monotone increasing sequence (Figur polynomial with five orders is used to fit the trend displacement in the training The least squares method is applied to find the best fitting values for the param the polynomial function. The polynomial function is as follows: The fitting polynomial of degree five and the prediction results of the trend d ments are shown in Figure 8 and Table 2. As exhibited in Figure 8, all the polyno degree five can accurately fit the increasing trend of landslide displacements and the following development in the test datasets. The R values of ZG01, ZG03, ZG0 ZG09, ZG11, and ZG12 on the test samples are 0.991, 0.991, 0.982, 0.989, 0.989, 0.9 0.990, respectively. The fitting results of the trend displacements are close to th tored displacement, mainly because the VMD can reduce the complexity of the data and obtain relatively stable displacement components.   (14)).

Periodic Displacement Prediction
According to the parameters setting the result in Section 4.4, the TCN prediction model was established. Figure 9 and Table 3 exhibit the prediction results produced by the trained models and compared with the measured values. It clearly shows that the prediction results of the TCN model have a similar fluctuation to the measured displacements. Most of the prediction values are close to the measured displacements. For the seven selected sites, the average values of RMSE, MAPE, and R on the test samples are 6.37 mm, 1.487%, and 0.911, respectively. The prediction results compared with the monitoring displacements indicate that the TCN model runs well on the periodic displacement prediction. In short, the TCN is a suitable tool for building nonlinear relationships between the external factors and the periodic displacement.      (Table 4). Figure 10 shows that the predicted cumulative displacements based on the polynomial of degree five and the TCN model are in good agreement with the measured displacement of the monitoring points. This is due to the hybrid models providing an integrated learning strategy, which can reduce the complexity of the time-series displacement, and relatively simple nonlinear and nonlinear models can realize high accuracy prediction. In the early monitoring period, the step-like displacement is obvious, but the step-like displacement feature is not significant at the end of the monitoring period. Compared with the different displacement status, we think the effect of the fluctuation of the reservoir water level on the slope will be weakened in the future. Other slopes with similar monitoring data in the Three Gorges Reservoir area also illustrate this point, such as the Baishuihe landslide [53] and the Qianjiangping landslide [54]. The surface displacement prediction is an import part of LEWS; the slope failure warning criterion is usually based on the displacement value and corresponding indices including the displacement rate, the displacement tangent angle and the inverse velocity method [55,56]. The displacement value cannot consistently produce an accurate slope failure warning for a certain landslide at present, due to the complexity and uncertainty of the geological conditions and external factors. However, these critical criteria of landslide failure can be obtained by accurate displacement prediction results, which can contribute to evacuating people from risk slop in advance and reduce landslide risk. So, in the future we will concentrate on the critical criterion of landslide failure, based on this method. displacement tangent angle and the inverse velocity method [55,56]. The displacement value cannot consistently produce an accurate slope failure warning for a certain landslide at present, due to the complexity and uncertainty of the geological conditions and external factors. However, these critical criteria of landslide failure can be obtained by accurate displacement prediction results, which can contribute to evacuating people from risk slop in advance and reduce landslide risk. So, in the future we will concentrate on the critical criterion of landslide failure, based on this method.

Comparing the 1D-CNN Model with the TCN Model
Although few studies introduce the CNN algorithm to perform landslide displacement prediction, it has useful applications in many other time-series prediction fields. So, in this competition, we constructed a 1D-CNN model for periodic displacement prediction, and three statistical metrics-RMSE, MAPE, and R-are used to evaluate the performance of the 1D-CNN and TCN models. Figure 11 and Table 5 show the specific comparison results. In Table 5, both the 1D-CNN and TCN models can provide relatively good accuracy of the periodic displacements, indicating that the 1D-CNN model also has potential to simulate landslide periodic displacement for LEWS. However, among 21 evaluation indices of the 7 points, 15

Comparison among TCN, LSTM, and SVM
To further evaluate the performance of the proposed model, the TCN model was compared with the SVM and LSTM models, which have performed well for displacement prediction too. The specific comparison results of the prediction performance indices are shown in Table 6, and the prediction values for the testing sample are shown in Figure 12. Generally, the three proposed models can predict the trend of the monitoring displacements. However, for the SVM model, the error and fluctuation compared with the monitoring displacement are larger than the other two models. From Table 6, compared with the static SVM model, the time-series modeling models (LSTM and TCN) are more competitive and reach a relatively high prediction accuracy. This result is consistent with the previous findings that the dynamic modeling method performed better than the static model [16,23,24]. Moreover, on most of the testing sites, the TCN model was outperformed by the other models. The performance of the LSTM model was lower than the TCN model in this competition. Specifically, compared with the SVM model, the average improvement ratios of RMSE, MAPE, and R of TCN were 39.58%, 66.69%, and 8.99%, respectively. For the LSTM model, the average improvement ratios of RMSE, MAPE, and R of TCN were 8.35%, 17.25%, and 0.56%, respectively.
In this competition, we also computed the time consumption for training the SVM, LSTM, and TCN models. The suitable SVM parameters were searched for using a metaheuristic SSA, and the parameters in the TCN and LSTM models were optimized based on the gradient descent algorithm (Adam). The training time for SVM (10.26 s) is lower than the LSTM (15.15 s) and TCN (11.57 s) models because there were only two parameters optimized by SSA, which was much less than the number of parameters in the two time-  To further evaluate the performance of the proposed model, the TCN model was compared with the SVM and LSTM models, which have performed well for displacement prediction too. The specific comparison results of the prediction performance indices are shown in Table 6, and the prediction values for the testing sample are shown in Figure 12. Generally, the three proposed models can predict the trend of the monitoring displacements. However, for the SVM model, the error and fluctuation compared with the monitoring displacement are larger than the other two models. From Table 6, compared with the static SVM model, the time-series modeling models (LSTM and TCN) are more competitive and reach a relatively high prediction accuracy. This result is consistent with the previous findings that the dynamic modeling method performed better than the static model [16,23,24]. Moreover, on most of the testing sites, the TCN model was outperformed by the other models. The performance of the LSTM model was lower than the TCN model in this competition. Specifically, compared with the SVM model, the average improvement ratios of RMSE, MAPE, and R of TCN were 39.58%, 66.69%, and 8.99%, respectively. For the LSTM model, the average improvement ratios of RMSE, MAPE, and R of TCN were 8.35%, 17.25%, and 0.56%, respectively. In this competition, we also computed the time consumption for training the SVM, LSTM, and TCN models. The suitable SVM parameters were searched for using a metaheuristic SSA, and the parameters in the TCN and LSTM models were optimized based on the gradient descent algorithm (Adam). The training time for SVM (10.26 s) is lower than the LSTM (15.15 s) and TCN (11.57 s) models because there were only two parameters optimized by SSA, which was much less than the number of parameters in the two time-series models. Compared with the LSTM model, the training consumption of TCN decreased by 23.63%. The main reason for this is that the parameters in the filters of the TCN model are shared across a convolution layer and a pooling layer, which can reduce the number of parameters. In particular, the LSTM is a sequence model, a lot of gate units were used to memorize the history information, and the calculation of the later state needs to wait for the former state calculation to be completed [25,57]. These results support the effectiveness of the proposed TCN model in processing the time series and improving the prediction accuracy, especially for periodic displacement. creased by 23.63%. The main reason for this is that the parameters in the filters of the TCN model are shared across a convolution layer and a pooling layer, which can reduce the number of parameters. In particular, the LSTM is a sequence model, a lot of gate units were used to memorize the history information, and the calculation of the later state needs to wait for the former state calculation to be completed [25,57]. These results support the effectiveness of the proposed TCN model in processing the time series and improving the prediction accuracy, especially for periodic displacement.

Conclusions
Based on the long-term professional GPS monitoring in the Muyubao landslide, a large-scale landslide in the TGRA of the chain, this study proposed a hybrid method of step-like landslide displacement prediction based on the TCN model, which consists of causal dilation layers and a residual block. The time-series analysis and the VMD theory were used to decompose the cumulative displacement, reduce the complexity of the data, and obtain relatively clear displacement patterns. For cumulative displacement prediction, the average values of RMSE, MAPE, and R of the proposed method are 15.89 mm, 0.85%, and 0.922, respectively. The results indicate that the prediction values of the proposed model are in good agreement with the actual values.
Moreover, this study focuses on the periodic displacement prediction, which is difficult to use to predict final prediction results. The periodic displacement prediction results indicated that the average RMSE (6.37 mm), MAPE (1.49%), and R (0.911) of the TCN model were the best, compared with the two commonly used machine learning and deep learning methods (SVM and LSTM). Compared with the static SVM model, the average improvement ratios range from 8.99% to 66.69%; compared with the sequence model LSTM, the accuracy of the proposed TCN model improved from 0.56% to 17.25%. Moreover, the training time of the TCN decreases by 23.63% compared with LSTM model. This is due to the TCN combined with the time-series modeling ability of RNN and CNN's global feature extraction. Convolution and pooling layers are widely applied in the TCN model, which can extract global features and reduce the number of parameters by using a weight-sharing strategy. The causal dilation layer in TCN can flexibly increase the receptive fields in the deeper layer, which contributes to the longer memory time series. Overall, the feasibility and superiority of the proposed hybrid method are demonstrated by successful implication in the seven landslides displacement time-series prediction. The proposed method can decompose cumulative displacement into relatively simple displacement components that are also suitable for other landslide displacement predictions characterized by slow-moving deformation and step-like displacements effected by periodical or seasonal factors.


extfieldsofstudy:
Computer Science
Mathematics

text:
SHREC 2022: Fitting and recognition of simple geometric primitives on point clouds

This paper presents the methods that have participated in the SHREC 2022 track on the fitting and recognition of simple geometric primitives on point clouds. As simple primitives we mean the classical surface primitives derived from constructive solid geometry, i.e., planes, spheres, cylinders, cones and tori. The aim of the track is to evaluate the quality of automatic algorithms for fitting and recognising geometric primitives on point clouds. Specifically, the goal is to identify, for each point cloud, its primitive type and some geometric descriptors. For this purpose, we created a synthetic dataset, divided into a training set and a test set, containing segments perturbed with different kinds of point cloud artifacts. Among the six participants to this track, two are based on direct methods, while four are either fully based on deep learning or combine direct and neural approaches. The performance of the methods is evaluated using various classification and approximation measures.

Introduction
Fitting and recognition are of paramount importance in multiple application domains, such as reverse engineering and data compression. In most real-world scenarios, acquisition methods tend to produce data affected by noise and potentially by other imperfections (e.g., non-uniform or low sampling density, misalignment and missing data). The literature on methods for fitting and recognition of simple primitives is quite consolidated [1] and recently we have seen the growth of data-driven methods, see for instance [2,3,4,5,6] resulting in a growing demand for appropriate datasets for training and their evaluation. It is therefore worth analyzing and comparing the outcome of simple primitive fitting methods on a controlled dataset with different kinds of perturbations and primitive variations and in terms of different evaluation measures.

Related datasets
Given the strong interest in the problem of fitting geometric primitives, several datasets are becoming available but, generally, they are not specifically designed for simple geometric primitives, possibly affected by perturbations, missing points and local variations. A first collection of one million CAD models targeted to data-driven problems was provided by [7], in which each object is associated to a ground truth for different quantities, as for example the indices of the patch segmentation and the analytic representation of surfaces and curves. A lot of surface types are considered and the models are represented as triangle meshes. Then, it is worth mentioning [8], a benchmark of point clouds representing CAD objects aimed at evaluating methods for detecting simple geometric primitives. Despite having some features similar to our purposes, such as sub-sampled clouds and missing parts, it is not big enough and was designed to test model segmentation methods. A dataset of mechanical components is proposed by [2], which generated the point samples from some models of the American National Standards Institute (ANSI) provided by TraceParts. In this collection, four types of primitives are considered (plane, sphere, cylinder, cone). Another synthetic dataset is presented in [9], but it proposes only cylinders, spheres and cubes as surface types. These existing datasets are not suitable for our purposes since our intent is to analyze the different methods for fitting geometric primitives on point clouds that present different types of specific artifacts. Indeed, in our dataset all point clouds have been artificially generated from exact surfaces and then optionally perturbed by simulating uniform or Gaussian noise of different intensity, undersampling, missing parts, small deformations, or a combination of such point cloud artifacts. In addition, we are not interested in the whole segmentation process, but we are focused on a sub-issue of it, that is the problem of classification and recognition of segments.

Contribution
The aim of this SHREC track is to evaluate the quality of automatic algorithms for fitting and recognising geometric primitives in point clouds. The goal is to identify, for each point cloud, its primitive type (i.e., planes, spheres, cylinders, cones and tori) and some geometric descriptors.
To this aim, a dataset of 3D segments, corresponding to parts of primitives and represented as point clouds, divided in training set and test set was generated. For each segment in the training set, its primitive type and an analytical representation of it were made available to the track participants. For each point cloud in the test set, the goal of the contest is twofold: to recognise the primitive type from the data and to provide an analytical representation of the primitive recognised. The relevance or non-relevance of the simple primitives has been evaluated on the basis of a ground truth that derived from the knowledge of the primitive that originated that segment. Finally, the performance of the methods is discussed using classification, recognition and approximation measures that quantify the adherence to the original primitives and the goodness of fit.

Organization
The remainder of this paper is organized as follows. Section 2 describes the characteristics of the proposed benchmark: the dataset, the ground truth, and the quality measures chosen to evaluate classification, approximation and recognition of primitives. Section 3 presents a description of the six automatic algorithms that participate to this SHREC track: two of them are fully-based on neural networks, two are direct methods, while the remaining ones use deep learning in an initial classification step and proceed by applying direct methods. Finally, Section 4 provides an accurate comparative analysis of the results obtained by the six methods.

Dataset and ground truth
The dataset is composed of 46, 925 three-dimensional segments represented as point clouds; it is divided into a training set and a test set that contain, respectively, 46, 000 and 925 point clouds. Each point cloud is provided in a TXT file listing one triplet per line.
The point clouds (see Figure 1) are generated by sampling classical surface primitives derived from constructive solid geometry (i.e., planes, cylinders, spheres, cones and tori), by means of the following procedure. First, point clouds are sampled by using parametric equations in their canonical form, i.e., centered at the origin of the coordinate axes and with rotational axis aligned with the z−axis; the geometric quantities that define each primitive (i.e., amplitudes and radii, if any) are assigned randomly. Second, to obtain segments of different shapes, several cuts are enforced by exploiting random planes. Third, we randomly apply translations and/or rotations to recover primitives in their general position. Lastly, each point cloud is (potentially) processed so that it can be: • A0 -Clean. No perturbation is applied.
• A1 -Perturbed by uniform noise of different intensities. The noise is obtained by sampling uniform distributions of the form U(− 1 n , 1 n ), being 3 ≤ n ≤ 20 random, and adding such perturbations to a random percentage of the points.
• A2 -Perturbed by Gaussian noise of different intensities. The noise is obtained by sampling normal distributions among N (− 1 n , 4 n 2 ), being 10 ≤ n ≤ 30 random, and adding such perturbations to a random percentage of the points.
• A3 -Clean but affected by undersampling. We randomly select a percentage of points to be removed.
• A4 -Clean but affected by missing parts. We randomly choose a point of the point cloud; then, we remove all points contained inside the sphere having the selected point as center and radius assigned randomly.
• A5 -Perturbed by uniform noise of different intensities and undersampled.
• A6 -Perturbed by Gaussian noise of different intensities and undersampled.
• A7 -Perturbed by uniform noise of different intensities and affected by missing parts.
• A8 -Perturbed by Gaussian noise of different intensities and affected by missing parts.
• A9 -Clean but with local deformations. We randomly select a point of the point cloud and we apply a bivariate Gaussian centered at the selected point with a random covariance matrix.
For the sake of brevity, we will often refer to these point cloud artifacts as perturbation types. For each perturbation type A0, . . . , A8, training set and test set contain, respectively, 5, 000 and 100 point clouds; for perturbation type A9, training set and test set count, respectively, 1, 000 and 25 point clouds.
• Cylinder. We list the primitive type (for cylinders, codified as 2), the radius r, the unit vector determining the rotational axis a, and a point P sampled on the axis. In this case, the file contains the vector v = [2, r, a t , P t ] t .
• Sphere. We list the primitive type (for spheres, codified as 3), the radius r and the center C. The ground truth corresponds to the vector v = [3, r, C t ] t .
• Cone. We list the primitive type (for cones, codified as 4), half the aperture 2 α, the unit vector determining the rotational axis a and the vertex C. Then, the file provides the vector v = [4, α, a t , C t ] t .
• Torus. We list the primitive type (for tori, codified as 5), the major and minor radii, respectively r and R, the unit vector determining the rotational axis a and the center C. In this case, the file contains the vector v = [5, r, R, a t , C t ] t .
For each point cloud, the track participants were required to return a TXT file per point cloud (of the test set), with the same syntax employed for the ground truth files of the training set.

Evaluation measures
Three tasks are considered in the evaluation phase: • Classification task. Methods are evaluated on the basis of their capability for identifying, for each point cloud of the test set, its primitive type.
• Recognition task. Methods are studied in terms of their ability to recover the geometric descriptors (e.g, radii of a torus, vertex of a cone) behind each segment.
• Fitting task. Regardless of the primitive type, methods' performance is studied with respect to the approximation of each point cloud.
To quantify the performance of automatic algorithms in classifying, recognizing and fitting simple geometric primitives on point clouds under different conditions, we therefore consider a variety of classification and approximation measures. This analysis will be further specified according to the type of geometric primitive and point cloud artifacts.

Classification measures
The confusion matrix [10] is a popular way to visualize the classification performance of a method. A confusion matrix CM is a square matrix whose order equals the number of classes in the dataset under study. Diagonal elements count true positives, i.e., all those items which have been correctly labeled as members of their ground truth classes: precisely, CM(i, i) is the number of items that have been correctly predicted as elements of class i. Off-diagonal elements give the numbers of items mislabeled by the classifier; precisely, CM(i, j), with j = i, is the number of elements wrongly labeled as belonging to class j rather than to class i. Ideals classifiers have diagonal classification matrices.
True Positive and Negative Rates True Positive Rate (TPR) refers to the method's ability to correctly identify positives (e.g., the percentage of cats correctly classified as cats). Likewise, True Negative Rate (TNR) measures the classifier's ability to identify negatives (e.g., the percentage of non-cats correctly classified as non-cats). In statistics, TPR and TNR are also called sensitivity and specificity. A perfect classifier is 100% sensitive and 100% specific.
Positive and negative predicted values Besides TPR and TNR, to quantify the likelihood that a method returns a true positive (or true negative) rather than a false-positive (or a false-negative) we consider the Positive Predictive Value (PPV) and the Negative Predictive Value (NPV). PPV is the ratio of true positives to the number of items labeled as positives by the classifier. Similarly, NPV is the ratio of true negatives to the number of items classified as negatives.
Accuracy Accuracy describes how often a classifier is correct in its predictions: precisely, it is defined as the proportion of predictions that the classifier got right.
Macro-average Macro-average permits to evaluate the overall performance of the classifier treating all classes equally, it is defined as the arithmetic mean of one of the metrics mentioned above, computed independently for each class.

Fitting and recognition measures
To measure the approximation accuracy of a specific primitive, we exploit the geometric descriptors predicted by each method and the corresponding parametric representation to sample densely the recognised surface primitive, say S P . Let us consider a point cloud P to be evaluated; we use the following two measures to evaluate the approximation accuracy: • Mean Fitting Error (MFE): where d is the minimum Euclidean distance between x and every point in S P , and l is the diagonal of the minimum bounding box containing P.
• Directed Hausdorff distance: with d the Euclidean distance.
Finally, the recognition accuracy is evaluated in terms of L 2 norm between the column vectors of the ground truth and that obtained by each method. Specifically, following the notation introduced in Section 2.1, we compare the geometric descriptors recognised with those of the ground truth considering the L 2 distance: where v G and v P here denote the ground truth and the prediction vectors, respectively, where we remove the first entry and, for planes and cylinders, the point (since they are not comparable).

Description of the methods
Ten groups initially subscribed to this track, with a final number of six methods submitted. In the remainder of this section we describe the methods, here simply referred to as M1, M2, . . . , M6. The six methods can be grouped as: • Fully direct methods: accelerated Hough Transform (M1), and Histograms of local features + parameter fitting (M2).

Overall strategy
The proposed approach is based on an accelerated version of the Hough transform, a well known technique in the pattern recognition community, which have been recently extended to algebraic objects, see [12]. The general (algebraic) HT-framework deals with the problem of finding a (hyper)surface, within a family of (hyper)surfaces dependent on a set of unknown parameters, that best approximates an input shape. More precisely, given an input point cloud and a family of parametrized primitives with some unknown parameters, the core idea is to consider the parameter space corresponding to the selected family of primitives, to discretize it into cells, and to vote a cell every time the Hough transform of an input point intersects it; the most voted cell uniquely identifies the optimal solution of the recognition problem. Despite the robustness to various point cloud artifacts (e.g., noise, outliers, missing data and uneven sampling), a severe limitation of this paradigm is due to the curse of dimensionality; indeed, the dimension of the parameter space corresponds to the number of parameters that are meant to be estimated through the voting procedure. Only a few families of primitives can be considered in their general position, i.e., planes and spheres, but even for ellipsoids it is required that the cloud is centred in the origin of the Cartesian axes [13] to make the problem computationally affordable.
The proposed method takes in input a point cloud P and, for each family of primitives F (among planes, cylinders, spheres, cones and tori), it works in three main steps.
• Pre-processing. Independently from the primitive type, the normal vectors of the points are computed using the Hough transform technique and the family of planes in the Hesse normal form, as described in [14]. The adoption of a voting strategy makes this computation robust to most point cloud artifacts. Once the candidate tangent planes have been estimated, the process is specialized for each type of primitive. Specifically, the rotations and/or translations necessary to bring the point cloud in standard form are based on the geometric properties of the primitive type in question. The result of this process is a new point cloud P and an initial estimates of the parameters of the primitive to be recognised. These parameters differ according to the type of family selected; for example in case of cylinders, tori or spheres, they correspond to the radii, while in case of cones they coincide with half the aperture.
• HT-based surface recognition. The new set of points P is the input of the classical HT-based recognition algorithm. The estimates obtained in the pre-processing step allow to localize the search in a dimensionally-reduced parameter space. The method returns the parameters of the surface that best fits P with respect to the given family of primitives. From these parameters, the geometric descriptors required by the track are obtained.
• Evaluation of the approximation accuracy. To evaluate the recognition accuracy of a specific primitive, the recognised surface is sampled from the parametric equation and the Mean Fitting Error as defined in Equation 1 is calculated between P and the sampled surface. Finally, the initial roto-translation is applied backwards in order to obtain the geometric descriptors for the point cloud P in its original position.
The three steps here described are repeated for each family of surfaces (that is, planes, cylinders, spheres, cones and tori). Then, the primitive type and the geometric descriptors of the surface with the lowest MFE is returned by the method. Figure 2 provides a graphical abstract of method M1.

Computational aspects
The tests were performed on a desktop equipped with an Intel Core i9 processor (at 3.6 GHz), in a Windows system 64 bits. The method was implemented in MATLAB. The execution time for each point cloud of the test set is, on average, 123 seconds.

Overall strategy
The proposed method pre-processes the input point cloud normalizing it in a unit cube. Then, the task is addressed with two main modules: a point cloud classification with k−nearest neighbours (k-NN) using a histogram of local features (Module 1), divided into several steps, and a shape parameter fitting (Module 2).
In the first step of Module 1, all objects are classified into two types: planar and non-planar segments. The latter are then classified into subcategories following two strategies: one considers cone and cylinder as two categories (Strategy 1), while the other on considers them as one category (Strategy 2). Both of them consider sphere and torus as one category. Finally, the primitive types are further specialized (to distinguish: tori from spheres for Strategy 1; cones, cylinders, tori and spheres for Strategy 2).
• Module 1: object shape classification. In the first step, both strategies search for all planes.
The Histogram of Oriented Normal Vectors (HONV, see [15]) is used to represent a point cloud, specifically, a 2D histogram (25 x 25). Each cell with a value greater than or equal to 80% of the maximum value in the histogram is shifted to the center of the 2D histogram. In this way, a multiple 2D histogram representation of a single object is created and the k−NN, with k = 3, is used to classify planar and non-planar objects. If a point cloud is classified as plane, the PCA is applied to further verify its planar property and to remove false-positive cases.
In the second step, the aim is to classify non-planar instances. Once the categories are defined accordingly to the strategy mentioned before, the HONV feature is used to represent the point cloud and the k−NN method is applied with k = 5. To classify sphere and torus, the Surflet-pair-relation histograms [16] are used in place of HONVs: in this way, each object can be represented with only one histogram. Finally, the k−NN with k = 15 is used to distinguish spheres from tori.
In Strategy 2, an additional step is required to classify cylinder and cone. In particular, the HONV feature is used to represent the point cloud and the k−NN with k = 5 is employed to differentiate between cylinder and cone.
• Module 2: finding shape parameters. For each primitive type, different approaches are considered to estimate and fine tune the shape parameters.
-Plane. The average value of all the points considered to be on the surface is calculated. Then, PCA is used to find the eigenvector corresponding to the smallest eigenvalue, which is taken as approximation of the normal vector.
-Sphere. A least squares sphere fit algorithm 3 is employed to find the center and radius of a sphere. In this case, only 30% of points are used and the process loops 100 times to find the parameters. Using the estimated parameter values, the number of points belonging to the surface of a sphere is computed, through the difference between the distance from the point to the center and the radius. If the difference is larger than 5% of the radius, that point is considered as an outlier. Finally, the set of parameter values with the most number of inliers is chosen.
-Cylinder. A cylinder fitting approach 4 is used to find the rotational axis and the radius by minimizing the common least squares loss.
-Cone. The RANSAC method proposed in [17] is followed to estimate the parameters for conic segments.
-Torus. The parameter values are initialised: the center in the origin, the rotation axis equivalent to the z−axis, the main radius 1 and the minor radius 0.1. First, the points are fitted like a 3D circle (find center, radius, axis of rotation) by optimizing the L 2 loss function. Then, the minor radius is calculated by averaging the distance of that circle to the set of points. Figure 3 provides a graphical abstract of method M2.

Computational aspects
The classification method needs at most 20 seconds on AMD Ryzen 5 3550H, when applying k−NN approach to classify all objects in each step of Strategy 1 and Strategy 2; this corresponds to an average of 0.02 seconds per model. Each point cloud might need multiple steps, and it takes at most 1 second to get the final primitive type prediction for an object. The shape parameter fitting methods take approximately 60 minutes on an Intel i5 10510U for all the 925 objects, which means about 3.89 seconds per model, on average.

Overall strategy
As for the previous method, M3 is based on a point cloud pre-processing, consisting of a centering and a scaling, and two consecutive steps: shape classification and parameter regression. This method considers the PointNet [18] architecture as backbone for all networks. A first neural network is used to return the class probabilities for the five classes: plane, cylinder, sphere, cone, and torus.
Depending on the predicted class, each point cloud is fed into one of five specific networks to obtain the primitive parameters. Networks differ in the final layer, whose outputs must depend on the primitive type; the last layer of the backend consists of 256 linear neurons, which means that the input for the last layer of the regression is a 256-dimensional vector. Each primitive parameter is thus generated by a different MLP head which consists of a linear layer. For example, the plane network has two heads in the final layer: a 256 × 3 linear layer for the normal and a 256 × 3 linear layer for the point in plane parameter.
When it comes to the classification task, the loss function used to train the network is the crossentropy loss between the one-hot activation functions from ground-truth and the network output. For the regression step, each primitive type has its own loss, defined as the sum of simpler functions: • The scalar loss L x,x compares a ground truth scalar parameters x with its predictionx, and is defined as the mean squared error. It can be used for radii and angles.
• The point loss L p,p is defined as the Euclidean distance between ground-truth point p and output pointp. It is used for points, e.g., for the sphere centers and the cone vertices.
• The vector loss controls the direction between a ground-truth vector v G and the output vector v O , and is defined as follow: It is used for normal vectors and rotational axes. The L 2 regularization between v G and v O enforces the componentwise similarity between vectors. The exponent of the cosine function is meant to improve the learning of normal directions.
The total loss of each primitive type is defined by summing the above-defined losses: Finally, the parameters are transformed with respect to the inverse translations and scaling applied in the preprocessing. The training is performed after splitting the original training set is split into two parts: a (sub-)training set (80%) and a validation set (20%). Figure 4 provides a graphical abstract of method M3.

Overall strategy
An initial preprocessing is applied to the input point clouds. First, they are translated to the origin, by subtracting the respective centroids. Following that, a normalization is applied, so that all points lie within the unit sphere. Random rotations about the x-, y-and z-axes are applied, for data augmentation purposes. The point clouds are finally converted into volumetric grids.
The backbone neural network for all models was originally proposed in the Minkowski Engine [19] as a candidate for the ModelNet classification challenge [20]. It is based on sparse voxel convolution, which minimizes the memory requirements associated with voxel based methods.
The network comprises of a voxel-wise input MLP and an input convolution. This is followed by three strided convolutions that progressively reduce the resolution of the voxel grid. After that, the output of the convolutional layers are projected on the initial grid, concatenated and finally passed through three strided convolutional layers. To extract the feature vector, Max and Mean pooling is applied to the resulting grid and the results are concatenated. All layers are followed by Batch Normalization and LeakyReLU activation function. For this track a 3-layer MLP head was used, with the number of output parameters changing according to the task.
M2's strategy involves two consequent steps: • Classification step. An instance of the model described above, trained for classification, is tasked to find which primitive the input points where sampled from.
• Fitting step. Five more copies of the same model, each one trained on a specific type of primitive, are considered for the fitting task. During this step, the primitive type is known, Figure 5: A graphical representation of the strategy adopted in method M4.
as it was predicted in the previous step. The model corresponding to that type of primitive is activated to regress the primitive's parameters. Although the network's performance is generally quite good, it can be increased even further by performing least-squares fitting. This is accomplished by using the network's prediction as an excellent initial guess, thereby making the fitting process very quick. For the least squares fitting, an existing python library geomfitty 5 is used and extended to include cones.
For the classification model, M2 minimizes the cross entropy loss. For the regression models, M4 uses a weighted combination of four basic losses, depending on the type of parameters to be predicted for a given primitive type: • Axis to axis loss. This loss is applied to predict direction vectors. Then, given two such vectors v 1 and v 2 the loss is defined as: Vectors are forced to have the same direction by minimizing the angle between them.
• Point to axis loss. This loss is used for anchor points that help to define the axes of various shapes. It is defined as the distance between the query point and the line that passes through the shape's axis.
• Point to point loss. This loss is used for the centers of spheres and tori. For a given point and a ground truth center it is defined as the Euclidean distance between them.
• Scalar to scalar loss. The remaining parameters such as radii and angles are simply scalars, so this loss is calculated via the mean squared error.
It is important to note that neural networks perform very poorly when tasked to predict numbers that are distributed over a wide range of values. The unit sphere normalization applied to the input helps with this problem, as it restricts said range to roughly (0, 1) for most types of parameters. In order to map these values back to the correct range we apply the inverse transformations that we applied to the input points.
The network was trained and evaluated on a 85% − 15% training-validation split of the given dataset and was found to perform very well, reaching over 95% validation accuracy. To further improve the performance, a backtracking strategy is implemented to deal with heavily perturbed primitives. After the first shape's parameters are regressed, the distance between the shape and the input points is measured. If the distance exceeds a pre-defined threshold, then it is assumed that Figure 6: A graphical representation of the strategy adopted in method M5.
our initial classification guess was wrong, so the network's guess with the second highest probability is taken and the process is repeated. For peculiar cases where the least-squares fitting diverges, the regressor network's predicted parameters is used as the final estimation. Figure 5 provides a graphical abstract of method M4.

Computational aspects
After experimenting on the training set, the voxel size was set to 0.05. For the minimization problem, M4 consider the optimization algorithm Adam [21].

Overall strategy
The input point clouds are pre-processed; specifically, each segment is re-sampled in order to obtain the same number of points (8100 points), as the number of points in each object is not constant. This strategy augments the data by holding two of the three (x, y, z) position and adds a small shift in the remainder coordinate to create new data points. Finally, normalization and random rotations about the z−axis are applied to increase the effectiveness of the training process. Similarly to M3, the core component for point cloud classification is the PointNet [18] architecture. This architecture has three main modules: the max pooling layer to aggregate information symmetrically from all the points symmetrically, a global and local feature combination structure and two joint alignment networks. The symmetrical max pooling layer helps to make model invariant to the orders of the input data. Input points are fed in a collection of multi-layer perceptron networks to extract global features. After computing the global point cloud feature vector, data are fed back to concatenate the global feature with each of the point features. M5's strategy works in two main steps: • Classification step. The data is trained by the above model with different parameters and the top three models are obtained with the highest accuracy in the validation test. The final label is returned after considering the majority votes from the models.
• Shape parameter fitting. For each primitive type, different ways are considered to estimate the shape parameters of each point cloud.
-Plane. The regression with least mean square loss is used to fit a 3D plane on the input point cloud.
-Sphere. The least mean square method is used to calculate the most appropriate sphere to the point cloud.
-Cylinder. The least square fitting method is used to fit the point clouds to an infinitely long cylinder, parameterized by an anchor point, the direction and the radius of the cylinder.
-Cone. Tangent planes are computed to approximate the vertex of the cone; then the circle that fits the set of points having the same distance from the vertex is used to find the rotational axis and half the aperture.
-Torus. The SciPy library is used to determine the torus that best fits the point cloud.
For the least squares fitting, the python library geomfitty is used (see method M4).
The networks are trained and evaluated by splitting the initial dataset in two parts, training set (80%) and validation set (20%); the validation accuracy is above 80%. Figure 6 provides a graphical abstract of method M5.

Computational aspects
In the experiments, a learning rate of 2.5 × 10 −4 is used, and Adam is used for optimization. The training is run on Google Colab Pro -GPU P100, for a total of 10 epochs. The average training times, per epoch, are: 50 minutes for the PointNet network and 8 min for all 3D geometry regressors. The execution time is, on average, 36 seconds for the classification of the 925 test segments, while it is 10 minutes to estimate their shape parameters.

M6: AlexNet and ISO reconstruction standards
This method is proposed by Yifan Qie and Nabil Anwer. The method implementation can be found at https://github.com/YifanQie/SHREC22_fitting_LURPA.git.

Overall strategy
Input point clouds are first preprocessed by centering them to the origin and scaling them into a unit sphere. Then, they are rotated so that the first principal axis of each point cloud lies along z−axis.
The primitive type identification of this method is based on AlexNet [22] which is a convolutional neural network architecture with 60 million parameters and 650,000 neurons. Segments are turned into images as input for training the neural network; images are collected from the principal directions for feature extraction and recognition. All the layers, except for the last three layers, are extracted form AlexNet. The updated last three layers consists of a fully connected layer, a softmax layer, and a classification output layer, where the fully connected layer is set to return as many numbers as the number of primitive types defined in contest (i.e., 5). Regarding the specificity of the input point clouds, transfer learning is used in the approach since the first layers in the pretrained networks have proved to be able to address simple pattern recognition in the images.
This strategy involves two main steps: • Classification step. The pre-trained neural network is used for the test set to provide the primitive type that identify the input point cloud.
• Fitting step. The surface reconstruction is tested using the methods considered in the ISO standards on geometrical product specifications and verification (GPS), detailed in [23]. Once the primitive type of a point cloud, that is the invariance class in ISO GPS terminology, is identified by the fine-tuned AlexNet, an ideal feature defined by a parametrized equation can be guaranteed using least square fitting to obtain a reconstructed surface, see [24].
The network was trained and evaluated on a 95% − 5% training-validation split of the given dataset, reaching over 87% validation accuracy. Figure 7 provides a graphical abstract of method M6.

Comparative analysis
The performance of each method presented in Section 3 is here quantitatively evaluated on the basis of the measures described in Section 2.2. For the sake of conciseness (and of the readers), we here report only tables referring to the whole test set; additional tables obtained by focusing on specific point cloud artifacts are provided as additional material in A and B. Figure 8 shows the confusion matrices over the whole test set. All the six methods have a high number of true positives, being cones and tori the most difficult primitive types to be recognized. A more thorough analysis is provided by Table 1, which contains the following information: Positive Predicted Value (PPV), Negative Predicted Value (NPV), True Positive Rate (TPR), True Negative Rate (TNR), Accuracy (ACC) and their macro-averages. We can notice that: • Methods M2, M3, M4 and M5 have NPV (slightly) higher than PPV: it is more likely for these methods to be right when reporting a negative rather than a positive. Methods M1 and M6 have mixed outcomes.

Classification task
• Methods M3 and M5 have TNR higher than TPR, making them more reliable in correctly finding true negatives rather than true positives. Methods M1, M2, M4 and M6 presents, for some primitive types, a TPR higher than the TNR: this means that, in some cases, the classifier is more likely to identify correctly true positives rather than true negatives.
• All methods have the accuracy over 90%, with the slightly lower scores for cylinders, cones and tori. Macro-averaged values of accuracy are always above 95%, with two of the three best scores reached by neural methods. Figure 9 clarifies the (global) impact of different point cloud artifacts upon the six methods, with respect to the classification measures. We observe that: • M1, M2, M4 and M6 reach their best performance on clean data. Surprisingly enough, M2, M4 and M6 perform best on data with small deformations: these three approaches appear to be able to extract some additional hidden information to improve their prediction capability.
• M1 suffers the most from undersampled data. On the other hand, coherently with the Hough paradigm on which the method relies on, it is not much affected by missing data.
• M3 and M4 have their lowest performance, respectively, in case of uniform and Gaussian noise. However, the grouped bar charts suggest these two methods are the least prone to misclassification. This result confirms the power of neural methods in handling classification tasks.  Table 1: Classification metrics for the whole test set, with respect to each primitive type: T1=plane, T2=cylinder, T3=sphere, T4=cone, T5=torus. Here, gold, silver and bronze identify the first, second and third best performance. The last column contains the macro averages.
• In spite of having, in most cases, the lowest classification performance, M2 has still no measure that goes below 80%. Compared to the other method using PointNet (M3), here the implementation seems to overfit way more the training data.
• M5 has a relatively high performance on undersampled data, when compared to other point cloud artifacts.
• The combination of noise and undersampling has proved to be the most challenging perturbation for M6 which, however, exhibits one of the best overall performance on small deformations.
By checking which segments are misclassified, it can be additionally observed that the segments which results being more problematic are, in many cases, the small ones. Tables 3-12 in A provide a complete overview of classification measures, specified with respect to the primitive type. Table 2 provides various statistics of the measures introduced in Section 2.2.2, when the whole test set is considered. It contains the following information: first, second and third quartiles, mean value and standard deviation. More specifically, these quartiles split a set of sorted real numbers into four parts of approximately equal cardinality: the first (Q1) and the third (Q3) quartiles are defined as the values such that, respectively, 25% and 75% of the numbers lie below them; the second quartile (Q2) is the median of the set. Quartiles' significance is due to their ability to identify possible outliers. Based on quartiles, we can notice that: Figure 9: Bar graphs of the macro averages for the classification measures, grouped by method. The legend reflects the color encoding of the perturbation types.

Recognition and fitting tasks
• The L 2 distance provides low values for all methods with regard to the three quartiles. The lowest Q1 has an order of magnitude of −14 and is obtained by method M6; the highest order of magnitude is that of method M3, and corresponds to −1. Q2 and Q3 have order of magnitudes between −3 and −1.
• With the sole exception of the third quartile of M5, quartiles of the MFEs lie under the order of magnitude of −2. The lowest values are obtained by M2.
• When it comes to quartiles, the Hausdorff distance the order of magnitude is generally (nonstrictly) lower than −1; the only exceptions fare the third quartiles of methods M3 and M5.
Based on the arithmetic mean and standard deviation, we can conclude that: • Methods M2, M5 and M6 are affected by numerical instability issues, which result in the presence of outliers; it is worth noting, however, this problem is limited to the top 25% of the errors.
• M2 presents 11 outliers corresponding to the recognition of cones in which the estimated halfangle estimate is close to zero and the estimated vertex is far away from the origin of the coordinate axes.
• M5 provides different outliers depending on the type of task. Specifically, it has 61 outliers in the parameter recognition of various cylinders, cones and tori; however, for these clouds, fitting errors are low. We therefore conclude that, in all these cases, the approximation is successful but at the cost of a poor recognition. From the fitting task point of view, MFE and directed Hausdorff distance have high values in correspondence of 18 outliers correctly identified as planes: here, the normal is correctly identified (thus, the low value of L 2 distance); on the other hand, the method fails in returning a point passing through each returned plane.
• M6 presents 30 outliers, which correspond to point clouds classified as cones; the problem these estimates suffer from are analogue to that of M2.
• M1 and M4 alternate the first and second places in terms of means and standard deviations, suggesting that direct and neural methods can both provide competitive solutions to the problem under study.
By checking the bar graphs of the log mean errors with respect to the perturbation type, see Figure 10, it is possible to note a generally low variation of the (average) performance of each method, with the exception of M2 and M6; for these results, the logarithmic scale was preferred because of the presence of outliers. It is indicative that, when the methods (i.e., all but M2 and M6) are right in predicting the primitive type, their estimates are not much influenced by the perturbation type. Moreover: • Methods M2 and M6 have better recognition and fitting performance when dealing with clean data or with point clouds having just small deformations. However, it is worth noticing once more that the poor performance mostly depends on the presence of outliers (see Table 2).
• M5 is on average better in fitting than in recognition: despite being the mean L 2 errors rather high, mean MFEs and directed Hausdorff distances are lower.
• M1 and M4 behave similarly in both recognition and fitting tasks.
As for classification measures, tables that specify recognition and fitting measures to the single point cloud artifacts are left to B (see Tables 13-22).

Concluding remarks
Six methods were submitted for evaluation to this SHREC track on the fitting and recognition of simple geometric primitives on point clouds. Two approaches -M1 and M2 -base their whole pipeline on direct methods; M3 and M4 are fully based in neural networks; M5 and M6 exploits both methodologies, that is, neural architectures for the classification task and direct methods for the fitting and recognition tasks. The dataset was specifically designed to allow the training of data-driven methods, with a training set of 46000 primitives and a test set of 925 models. Nine types of perturbations were randomly applied to each type of primitive in the ratio of 5000 training examples and 100 test samples, with the exception of local deformations that count 1000 training examples and 25 test samples. The outcome of the methods shows that, on simple primitives, data-driven methods are reaching a satisfactory level of maturity, obtaining the best scores both in primitive classification and in the parameter estimation. Direct methods remain fiercely competitive, and seem to suffer more in classification than in fitting and recognition. We believe that a similar comparative analysis will be possible in the future on models with multiple primitives.
As a possible future development, to better assess the potential of these methods in a reverse engineering context, it would be interesting to create a dataset with real data, e.g., captured with a high quality industrial 3D scanner, and derive the ground truth coming from the specifications used in the fabrication of the targets. To determine which approach (or which paradigm, between direct and neural methods) has the best generalization capability, it would be certainly worthwhile to keep the training set introduced in this track and use real data just in the test set. This direction is however rather challenging, as it requires having sufficient diversity, in terms of models to be scanned, between the primitives and parameters used for training and testing sets. Method  T1  T2  T3  T4  T5  Avg  M1 1