Wavelet Packet Method for Locating Critical Slip Surface Using the Strength Reduction Method

When the finite element-strength reduction method is used for two-dimensional slope stability analysis for elastic-perfectly plastic material, the failure criterion usually adopts the criterion of plastic zone penetration. That is, when the slope is in the limit equilibrium state, the plastic zone goes through the slope from the toe to the top. Meanwhile, the critical slip surface is composed of a series of points of maximum equivalent plastic strain along the depth direction. By deploying a set of parallel lines approximately perpendicular to the slope surface and picking out the points of these lines with the maximum equivalent strain points, we obtain a series of points taking on a wave shape, which constitutes a signal function. Subsequently, the wavelet packet analysis is used to smooth these points, i.e., locating the critical slip surface. The analysis of classic examples and comparison with Spencer’s method show that the proposed method in this paper is reasonable and effective.


Introduction
The finite element method is the most widely used method for slope stability analysis (SSA). Its strong ability in nonlinear analysis gives the finite element method (FEM) great potential in geotechnical engineering. It was first introduced to geotechnical engineering by Clough and Woodward [1]. In SSA, the factor of safety (FOS) is defined as the ratio of the actual soil shear strength to the minimum shear strength required to prevent failure [2]. As Duncan [3] pointed out, the FOS is the factor that continuously reduces the soil shear strength to bring the slope to the verge of failure. This method, which brings slope into a limit equilibrium state (LES) by reducing the shear strength of the soil, is called the strength reduction method (SRM), was first used in SSA by Zienkiewicz et al. [4], and has been studied by Naylor [5], Wong [6], Donald et al. [7], Matsui et al. [8], Duncan, and others. Griffiths et al. systematically summarized the advantages of the finite element method-strength reduction method compared with limit equilibrium methods (LEM) [9] in SSA; more detail can be found in [10].
In practice, there are two crucial issues in the SSA: calculation of the FOS and location of the critical slip surface (CSS). However, research into SRM is mainly concentrated on the computation of the FOS rather than locating the CSS. By bringing the slope to LES, the FOS agreeing with the LEM can be easily obtained. In contrast, the CSS of a slope is usually determined approximately through visualization techniques, such as the distribution of maximum shear-strain rates [4], shear-strain increment distributions [8], yielding regions of the slope [38], the velocity field at collapse [39], or other technical measures. However, to this day, there is no rigorous mathematical method to locate the CSS.
For two-dimensional SSA, Zheng et al. [40] developed a new approach to search the CSS for any yield criterion by using the equivalent plastic strain field in the LES, in which the least-square method is used to fit the CSS. However, for a complicated slope composed of multiple materials, its CSS shape may contain several turning points [41], and the leastsquare method will not be applicable in this case. For this reason, the wavelet method [42] was proposed and applied to the analysis of homogeneous and multi-layer slopes, in which the maximum equivalent plastic strain point is regarded as a signal function, and the high-frequency part of the signal is filtered to obtain a smooth CSS. The results show that the wavelet method has good applicability. Wavelet analysis is a rapidly developing new field in applied mathematics and engineering [43]. It can be used for damage detection [44], time-frequency analysis [45], signal-to-noise separation [46], multi-scale edge detection [47], and so on. Wavelet packet analysis (WPA) divides the time-frequency more carefully, and its resolution of the high-frequency part of the signal is higher than that of the dyadic wavelet. In other words, the dyadic wavelet only further decomposes the lowfrequency part into low-frequency and high-frequency parts, while wavelet packet decomposition further decomposes both low-frequency and high-frequency parts as shown in Figure 1. Hence, the wavelet packet has a wide range of application values [48]. The application of new mathematical methods for the solution of engineering problems can not only introduce the methods but also shorten the distance between the methods and applications. In this study, the SRM and the WPA are combined for the SSA. As in [42], the distribution of maximum equivalent plastic strain points in this study is regarded as a signal function, which is filtered by wavelet packet denoising technology to obtain a smooth CSS. The implementation step is as follows: First, the SRM is used to bring the slope into the LES and calculate the FOS. Second, deploying a series of parallel lines approximately perpendicular to the slope surface, and picking out the points of these lines with the maximum equivalent plastic strain points, is noted as y i (x i ). Third, the wavelet packet filter is used to smooth y i (x i ) and the location of the CSS is accordingly determined. Three numerical examples including a homogeneous slope, a slope containing a weak intercalated layer, and a soil-rock mixture slope are employed to validate the proposed method.

Slope Stability Analysis by the SRM
For SRM, there may be a significant difference in the failure loads or FOS if different failure criteria are selected [49]. At present, there are several common failure criteria: (1) non-convergence of the numerical solution [10,38,50]; (2) a sudden increase in the deformation of a characteristic location [51]; (3) connectivity of the shear increment strain or shear-strain region [8]. Since many factors may lead to the non-convergence of this numerical calculation, non-convergence does not necessarily mean the failure of the structure, thus the first criterion is considered insufficient. For the second criterion, this method needs to constantly draw the displacement-reduction factor curve to judge the position of the turning point through observation, so this criterion is also considered not objective enough.
For an elastic-perfectly plastic model with an associated flow rule, when a point within the slope is in the plastic state, Zheng et al. [52] proved that the elastoplastic matrix D ep at that point must be a positive semidefinite matrix with a rank defect of one and found out the direction of the singularity of D ep . Based on this fact, Zheng suggested that the failure criterion of the slope should be considered the moment when the strip runs through from the toe to the top of the slope, and all the elements in the strip are in the plastic state. Therefore, the third criterion, i.e., plastic zone penetration, was chosen as the failure criterion in this study, and the assumptions are adopted too.
In practice, we take the reduction factor F i corresponding to the LES as the FOS. Once the failure criterion is selected, the FOS and the distribution of the equivalent plastic strain (EPS) in the LES can be easily obtained by RSM. The EPS is defined as follows: where ∆ε p is incremental plastic strain vector in a load; the sum is performed with all load steps.
In this paper, the failure criterion based on plastic zone penetration is adopted. Once a slope arrives at the LES, a contour line of ε p with a small value, it passes through the slope from the toe to the top. The larger ε p of a contour line, the larger the plastic deformation of the closed area of the contour line, and the more serious the slope damage is [41]. Therefore, it can be known that the CSS goes through the points at which the equivalent plastic strain obtains the maximum in the vertical direction. Due to the interpolation characters of the isoparametric elements, the maximizers will generally not be on a smooth curve but take the form of a wave. Based on this observation, a new method was proposed by Zheng et al. for searching the CSS.

Wavelet Transformation and Reconstruction
A wavelet is a special kind of function with oscillation and can quickly decay to zero, and tends to be irregular and asymmetric [53]. In general, the mother wavelet ψ(t) is restricted, and function ψ(t) should be governed as follows whereψ(ω) is Fourier transform of ψ(t). By the dilation and translation of ψ(t), the wavelet basis ψ a,b (t) can be defined as follows where the parameter a > 0 indicates the scale factor and the parameter b∈R represents the displacement factor.
For the signal function f (t) ∈ L 2 (R), after wavelet decomposition, the coefficients are where ψ * (t) is the complex conjugate function of ψ(t).
The inverse wavelet transform formula is Defining a = 2 −j , b = k2 −j , and j, k ∈ Z, the corresponding binary discrete wavelet function sequence ψ j,k can be formulated as and when f (t) ∈ L 2 (R), the discrete wavelet transform is

Multi-Resolution Analysis
Mallat first proposed a multi-resolution analysis (MRA) and presented a fast algorithm based on wavelet transformation and reconstruction, named the Mallat algorithm. According to the Mallat algorithm, if f k contains the sampling information of f (t), and a 0,k = f k , the orthogonal wavelet decomposition formula is a j,k = ∑ n a j−1,k h(n − 2k), k = 0, 1, 2, . . . , N − 1 d j,k = ∑ n a j−1,k g(n − 2k), k = 0, 1, 2, . . . , N − 1 (8) where c j,k is the scaling coefficients, d j,k is the wavelet coefficients, j is the layer of decomposition, and N is the sampling points. In this, h and g are orthogonal filters generated by the scale function φ(t) and wavelet function ψ(t), respectively. These are defined as The reconstruction procedure is an inverse transform of the decomposition. The reconstruction formula is a j−1,k = ∑ n a j,n h k−2n + ∑ n d j,n g k−2n (10) A wavelet has the property of multi-resolution, which can treat signal data in multiresolution analysis and decompose a signal interwoven with different frequencies into different frequency sub-signals.

Wavelet Packet Transform and Reconstruct
In the MRA of signal function, the original signal is first decomposed into two parts, low frequency and high frequency, and then the low-frequency part is further decomposed into low frequency and high frequency, and then it is decomposed gradually. However, in wavelet packet decomposition (WPD), each scale is broken down into high frequency and low frequency simultaneously, so the resolution is higher.
The wavelet package {u n } n∈Z + can be created by scale function φ(t) and wavelet function ψ(t). The definition is as follows where n is a wavelet packet parameter, h and g are orthogonal filters [54], u 0 = φ, and u 1 = ψ. Binary scaling and translation can also be performed for u n u j,n,k (t) = 2 −j/2 u n 2 −j t − k For orthogonal wavelet packet analysis, the algorithm is very similar to the Mallat algorithm of orthogonal wavelet analysis. The following recursive formula can achieve the wavelet packet transform where d j,n (k) is the wavelet packet coefficient, defined by The reconstruction wavelet packet procedure is an inverse transform of the decomposition. The reconstruct recursive equation is The essentials of wavelet packet decomposition are to divide a signal into different frequency spans layer-by-layer. If the original signal contains enough information, the frequency span can be divided very subtly.
The equation deduces that if the decomposition layer is chosen appropriately, we can obtain the beginning and end of expected frequency spans-and with this method, the original signal when every frequency is wider enough [55]. Thus, there is no need to locate the frequency spans accurately for the ingredients which we needed so that it can adjust frequency drift.
The advantage of reconstruction is that we can choose all the frequency bins or some parts of them and set the others (noise or random interruption) as zero. Thus, if we can decompose a signal into different frequency bins, it is much easier to extract noise during the reconstruction. Figure 1 illustrates the three-layer wavelet and wavelet packet decomposition, where a represents the low-frequency part and d represents the high-frequency part.
According to the entropy standard, a set of orthogonal bases selected from the wavelet packet is called a wavelet packet basis of L 2 (R). Different wavelet packet bases result in different signal decomposition. The optimal wavelet tree can be selected according to the entropy standards [56]. Matlab software provides Shannon entropy, entropy threshold (threshold), norm entropy (norm), Sure entropy, log energy entropy (log energy), and userdefined entropy. Among these standards, Shannon entropy is used to select the optimal wavelet tree in this article.

Threshold Processing of Wavelet Packet Decomposition Coefficients
After selecting the optimal wavelet packet tree, we can use a threshold to remove the noise in wavelet packet nodes. There are many threshold denoising methods in practice. The commonly used threshold functions are soft threshold and hard threshold, and the calculation formulas are as follows: Soft threshold:ŵ Hard threshold:ŵ Based on Matlab software, we can use the "ddencmp.m" function to obtain the default threshold value of the signal data. The methods to obtain these values include the soft and hard threshold methods. We can also use the "wpthcoef.m" function to deal with the threshold values of wavelet packet coefficients [57].

Searching the CSS Based on the Wavelet Packet
In finite element analysis, due to the interpolation characters of the isoparametric elements, y i (x i ) is not on a smooth curve but in the form of a wave. Each line has only one maximum point, and these points constitute a set of good signal data y i (x i ). Our task is to use a wavelet packet to smooth these points and locate the CSS.
The algorithm steps of the proposed method are as follows: (1) Combined with FEM, the SRM is used to bring the slope into the LES and obtain the FOS of the slope. The obtained data points are decomposed by wavelet packet. According to the given signal data, a suitable wavelet basis function ψ a,b and the decomposition level is determined. Although there exist a huge number of possible wavelet bases, this paper is limited to the Daubechies orthogonal wavelet family. After the test, the Daubechies six-coefficient wavelets were used in all work. (4) The wavelet packet coefficients are processed by threshold quantization. According to the distribution of the data and its high-frequency part, each frequency band's wavelet packet decomposition coefficients are processed by a threshold. (5) Wavelet packet reconstruction is performed. The processed wavelet packet decomposition coefficients are reconstructed by wavelet packet, and the location of the CSS is obtained by connecting the corresponding data points.
The implementation of the proposed method is summarized in the flow chart as in Figure 2.

Numerical Examples
Although more complicated constitutive models can be used in finite elemen sis, the most commonly used model in engineering analysis is the elastic-perfectl model based on Mohr-Coulomb's yield criterion [41]. Therefore, in the following ical example analysis, the Mohr-Coulomb's yield surface and the associated fl were used in the material model. Pore pressure was not considered. Based on th strain assumption, the ANSYS software platform was used to simulate the slop four-node isoparametric element.

An Idealized Homogeneous Slope
To verify the proposed method, a homogeneous slope studied in [58] was a which was composed of homogeneous, isotropic soil. The dimensions and the parameters of the slope are shown in Figure 3. Four discrete models were used to the convergence of the proposed method, as shown in Figure 4, in which the bo conditions of the computational models were as follows: the left and right bou were fixed in the horizontal direction, and the bottom boundary was fixed in the h tal and vertical direction.  20m

Numerical Examples
Although more complicated constitutive models can be used in finite element analysis, the most commonly used model in engineering analysis is the elastic-perfectly plastic model based on Mohr-Coulomb's yield criterion [41]. Therefore, in the following numerical example analysis, the Mohr-Coulomb's yield surface and the associated flow rule were used in the material model. Pore pressure was not considered. Based on the plane strain assumption, the ANSYS software platform was used to simulate the slope with a four-node isoparametric element.

An Idealized Homogeneous Slope
To verify the proposed method, a homogeneous slope studied in [58] was adopted, which was composed of homogeneous, isotropic soil. The dimensions and the material parameters of the slope are shown in Figure 3. Four discrete models were used to analyze the convergence of the proposed method, as shown in Figure 4, in which the boundary conditions of the computational models were as follows: the left and right boundaries were fixed in the horizontal direction, and the bottom boundary was fixed in the horizontal and vertical direction.

Numerical Examples
Although more complicated constitutive models can be used in finite elemen sis, the most commonly used model in engineering analysis is the elastic-perfectly model based on Mohr-Coulomb's yield criterion [41]. Therefore, in the following ical example analysis, the Mohr-Coulomb's yield surface and the associated flo were used in the material model. Pore pressure was not considered. Based on th strain assumption, the ANSYS software platform was used to simulate the slope four-node isoparametric element.

An Idealized Homogeneous Slope
To verify the proposed method, a homogeneous slope studied in [58] was a which was composed of homogeneous, isotropic soil. The dimensions and the m parameters of the slope are shown in Figure 3. Four discrete models were used to the convergence of the proposed method, as shown in Figure 4, in which the bo conditions of the computational models were as follows: the left and right bou were fixed in the horizontal direction, and the bottom boundary was fixed in the h tal and vertical direction.   The EPS contours by the SRM using different models are plotted in Figure 5. The results indicate that the denser the mesh, the smaller the areas of the plastic zone. This is because of the mesh dependence of the FEM, which can be further proved by the distribution of maximum equivalent plastic strain, that is, the calculation model with denser mesh has lower data fluctuation. The traditional SRM can only give the contours of equivalent plastic strain as shown in Figure 5, and the specific position of the CSS cannot be given. The EPS contours by the SRM using different models are plotted in Figure 5. The results indicate that the denser the mesh, the smaller the areas of the plastic zone. This is because of the mesh dependence of the FEM, which can be further proved by the distribution of maximum equivalent plastic strain, that is, the calculation model with denser mesh has lower data fluctuation. The traditional SRM can only give the contours of equivalent plastic strain as shown in Figure 5, and the specific position of the CSS cannot be given.
Taking the slope vertex node as a characteristic point, the reduction factor Z i versus horizontal displacement u x curve is drawn as shown in Figure 6. As can be seen from the figure, the horizontal displacement of the slope vertex node has an obvious inflection point. Taking the slope vertex node as a characteristic point, the reduction factor Zi versus horizontal displacement ux curve is drawn as shown in Figure 6. As can be seen from the figure, the horizontal displacement of the slope vertex node has an obvious inflection point.        The slope toe was defined as the coordinate origin, X-axis along the slope surface, and Y-axis perpendicular to the slope surface. The distribution of the maximum equivalent plastic strains at the LES is illustrated in Figure 8, along with a direction perpendicular to the slope surface with an interval of 0.5 m. After the wavelet packet decomposition of the signal function, it could be observed that the high-frequency part of the signal was mainly concentrated in the first three orders. The CSS is obtained through a wavelet packet filter with db6 at the third scale level and the CSS in the LES are given in Figure 9. As shown in Figure 9, using different mesh density calculation models, the CSSs obtained by the proposed method are in good agreement with the Spencer method. The results show that the proposed method can overcome the mesh dependence of FEM, and it is reasonable to take the point with the maximum equivalent plastic strain point as the location of the CSS. lent plastic strains at the LES is illustrated in Figure 8, along with a direction perpendicular to the slope surface with an interval of 0.5 m. After the wavelet packet decomposition of the signal function, it could be observed that the high-frequency part of the signal was mainly concentrated in the first three orders. The CSS is obtained through a wavelet packet filter with db6 at the third scale level and the CSS in the LES are given in Figure 9. As shown in Figure 9, using different mesh density calculation models, the CSSs obtained by the proposed method are in good agreement with the Spencer method. The results show that the proposed method can overcome the mesh dependence of FEM, and it is reasonable to take the point with the maximum equivalent plastic strain point as the location of the CSS.

Slope Containing Weak Intercalated Layer
This example is Problem 3a issued by Australia Computer Association and Design Society (ACADS) [58][59][60][61] and is selected to verify that the proposed method can also accommodate inhomogeneous slopes. The slope geometry and the soil properties are illustrated in Figure 10. The mesh generation of the slope model is shown in Figure 11, and its boundary displacement constraints are the same as in Example 1. The contours of EPS by the SRM are plotted in Figure 12. As shown in Figure 12, the plastic zone first appeared in the weak interlayer, and the areas of the plastic zone increased as the reduction factor Fi increased. As Figure 13 shows, if a group of parallel lines are arranged with an interval of 0.6 m, the distribution of the maximum equivalent plastic strains is irregular. The analysis results are shown in Figures 13 and 14 with db6 at the third scale level and are consistent with the Spencer method. According to the proposed method, the FOS is 1.31, which is close to the reference results of 1.24-1.27 given by ACADS [58].
16.5m lar to the slope surface with an interval of 0.5 m. After the wavelet packet decomposition of the signal function, it could be observed that the high-frequency part of the signal was mainly concentrated in the first three orders. The CSS is obtained through a wavelet packet filter with db6 at the third scale level and the CSS in the LES are given in Figure 9. As shown in Figure 9, using different mesh density calculation models, the CSSs obtained by the proposed method are in good agreement with the Spencer method. The results show that the proposed method can overcome the mesh dependence of FEM, and it is reasonable to take the point with the maximum equivalent plastic strain point as the location of the CSS.

Slope Containing Weak Intercalated Layer
This example is Problem 3a issued by Australia Computer Association and Design Society (ACADS) [58][59][60][61] and is selected to verify that the proposed method can also accommodate inhomogeneous slopes. The slope geometry and the soil properties are illustrated in Figure 10. The mesh generation of the slope model is shown in Figure 11, and its boundary displacement constraints are the same as in Example 1. The contours of EPS by the SRM are plotted in Figure 12. As shown in Figure 12, the plastic zone first appeared in the weak interlayer, and the areas of the plastic zone increased as the reduction factor Fi increased. As Figure 13 shows, if a group of parallel lines are arranged with an interval of 0.6 m, the distribution of the maximum equivalent plastic strains is irregular. The analysis results are shown in Figures 13 and 14 with db6 at the third scale level and are consistent with the Spencer method. According to the proposed method, the FOS is 1.31, which is close to the reference results of 1.24-1.27 given by ACADS [58].
16.5m Figure 9. Critical slip surfaces based on wavelet packet and Spencer for Example 1.

Slope Containing Weak Intercalated Layer
This example is Problem 3a issued by Australia Computer Association and Design Society (ACADS) [58][59][60][61] and is selected to verify that the proposed method can also accommodate inhomogeneous slopes. The slope geometry and the soil properties are illustrated in Figure 10. The mesh generation of the slope model is shown in Figure 11, and its boundary displacement constraints are the same as in Example 1. The contours of EPS by the SRM are plotted in Figure 12. As shown in Figure 12, the plastic zone first appeared in the weak interlayer, and the areas of the plastic zone increased as the reduction factor Fi increased. As Figure 13 shows, if a group of parallel lines are arranged with an interval of 0.6 m, the distribution of the maximum equivalent plastic strains is irregular. The analysis results are shown in Figures 13 and 14 with db6 at the third scale level and are consistent with the Spencer method. According to the proposed method, the FOS is 1.31, which is close to the reference results of 1.24-1.27 given by ACADS [58]. lar to the slope surface with an interval of 0.5 m. After the wavelet packet decomposition of the signal function, it could be observed that the high-frequency part of the signal was mainly concentrated in the first three orders. The CSS is obtained through a wavelet packet filter with db6 at the third scale level and the CSS in the LES are given in Figure 9. As shown in Figure 9, using different mesh density calculation models, the CSSs obtained by the proposed method are in good agreement with the Spencer method. The results show that the proposed method can overcome the mesh dependence of FEM, and it is reasonable to take the point with the maximum equivalent plastic strain point as the location of the CSS.

Slope Containing Weak Intercalated Layer
This example is Problem 3a issued by Australia Computer Association and Design Society (ACADS) [58][59][60][61] and is selected to verify that the proposed method can also accommodate inhomogeneous slopes. The slope geometry and the soil properties are illustrated in Figure 10. The mesh generation of the slope model is shown in Figure 11, and its boundary displacement constraints are the same as in Example 1. The contours of EPS by the SRM are plotted in Figure 12. As shown in Figure 12, the plastic zone first appeared in the weak interlayer, and the areas of the plastic zone increased as the reduction factor Fi increased. As Figure 13 shows, if a group of parallel lines are arranged with an interval of 0.6 m, the distribution of the maximum equivalent plastic strains is irregular. The analysis results are shown in Figures 13 and 14 with db6 at the third scale level and are consistent with the Spencer method. According to the proposed method, the FOS is 1.31, which is close to the reference results of 1.24-1.27 given by ACADS [58].              Figure 15 illustrates the geometric model size of the soil-rock mixture slope [12] and the material parameters are given in Table 1. The rock gradation of the soil-rock mixture slope is shown in Table 2. Figure 15 illustrates the geometric model size of the soil-rock mixture slope [12] and the material parameters are given in Table 1. The rock gradation of the soil-rock mixture slope is shown in Table 2.   Figure 16 illustrates the mesh of the soil-rock mixture slope in ANSYS. By encrypting the mesh near the slope, more data points are available to be searched by the algorithm to improve the accuracy of the calculation and reduce the fluctuation of the equivalent plastic strain maximum distribution line, thus improving the smoothness of the slip surface by filtering with wavelets and reflecting more accurately the position of the CSS in the slope.    Figure 16 illustrates the mesh of the soil-rock mixture slope in ANSYS. By encrypting the mesh near the slope, more data points are available to be searched by the algorithm to improve the accuracy of the calculation and reduce the fluctuation of the equivalent plastic strain maximum distribution line, thus improving the smoothness of the slip surface by filtering with wavelets and reflecting more accurately the position of the CSS in the slope. Figure 15 illustrates the geometric model size of the soil-rock mixture slope [12] and the material parameters are given in Table 1. The rock gradation of the soil-rock mixture slope is shown in Table 2.   Figure 16 illustrates the mesh of the soil-rock mixture slope in ANSYS. By encrypting the mesh near the slope, more data points are available to be searched by the algorithm to improve the accuracy of the calculation and reduce the fluctuation of the equivalent plastic strain maximum distribution line, thus improving the smoothness of the slip surface by filtering with wavelets and reflecting more accurately the position of the CSS in the slope.  For this numerical example, the traditional LEM was no longer used. Therefore, the results of the proposed method were compared with the virtual element method [12]. The EPS contours by SRM are plotted in Figure 17. The plastic zone first appeared at the toe of the slop. With the increase in the reduction factor F i , the stress concentration gradually appeared at the toe of the slope and around the stone, and finally, the plastic zone bypassed the stone. The failure mechanism was different from that of the soil slope. With an interval of 0.5 m, Figure 18 displays the distribution of the maximum equivalent plastic strains and CSS. Due to the traditional SRM, the specific location for CSS was also not given in [12]. A wavelet packet decomposition with a db4 scale of a fourth was performed on the data for this example to calculate the curve of the CSS. An FOS of 1.09 was generated, which is close to the 1.07 evaluated by Sun et al. [12]. The filtered fitted curve reflects the influence of the stones on the position of the CSS in stony soil slopes.

An Example of Soil-Rock Mixture Slope
results of the proposed method were compared with the virtual element method [12]. The EPS contours by SRM are plotted in Figure 17. The plastic zone first appeared at the toe of the slop. With the increase in the reduction factor Fi, the stress concentration gradually appeared at the toe of the slope and around the stone, and finally, the plastic zone bypassed the stone. The failure mechanism was different from that of the soil slope. With an interval of 0.5 m, Figure 18 displays the distribution of the maximum equivalent plastic strains and CSS. Due to the traditional SRM, the specific location for CSS was also not given in [12]. A wavelet packet decomposition with a db4 scale of a fourth was performed on the data for this example to calculate the curve of the CSS. An FOS of 1.09 was generated, which is close to the 1.07 evaluated by Sun et al. [12]. The filtered fitted curve reflects the influence of the stones on the position of the CSS in stony soil slopes.

Conclusions
To accurately locate the position of the CSS of a complicated slope in LES, a new method is proposed by combining the WPA and the SRM. The algorithm proposed in this paper includes two approaches. The first approach is an algorithm to bring the slope into the LES. Based on the SRM, the strain field of LES and FOS can be found. The second approach is the wavelet packet filter algorithm for smoothing the CSS. This algorithm applies wavelet packets to SSA. Through the analysis of three different examples and comparison with Spencer's method, the results show that: (1) the proposed method is reasonable and efficient and can be used to locate the CSS; (2) it is reasonable to take the maximum equivalent plastic strain point as the position of CSS. It should be pointed out that this study is limited to two-dimensional SSA with elastic-perfectly plastic material. The validation of the proposed method for the three-dimensional problems and the non-elastic-perfectly plastic material needs further study.

Conclusions
To accurately locate the position of the CSS of a complicated slope in LES, a new method is proposed by combining the WPA and the SRM. The algorithm proposed in this paper includes two approaches. The first approach is an algorithm to bring the slope into the LES. Based on the SRM, the strain field of LES and FOS can be found. The second approach is the wavelet packet filter algorithm for smoothing the CSS. This algorithm applies wavelet packets to SSA. Through the analysis of three different examples and comparison with Spencer's method, the results show that: (1) the proposed method is reasonable and efficient and can be used to locate the CSS; (2) it is reasonable to take the maximum equivalent plastic strain point as the position of CSS. It should be pointed out that this study is limited to two-dimensional SSA with elastic-perfectly plastic material. The validation of the proposed method for the three-dimensional problems and the nonelastic-perfectly plastic material needs further study.