Abstract
The strength reduction method embedded in a distinct element code was used to analyse the stability of a slope in a coal mining area that had been reinforced twice, primarily with pile and retaining wall, followed by porous steel-tube bored grouting. For the primary reinforcement, the factor of safety was calculated, slip surface and failure mechanism were determined, and the damage phenomenon of primary reinforcement was analysed in detail. Failure time of slope without further strengthening was predicted by applying a new quantitative method based on monitoring displacement data. The slope instability at the primary reinforced stage was verified by these analyses. For the second reinforcement, the effect was evaluated by combining the new factor of safety and the final monitoring data, which validates the slope stability. Especially, variations of displacement and factor of safety due to water influence are analysed. Through this procedure, a systemic method for the slope safety evaluation and assurance is presented for engineering practice reference.
1. Introduction
The soil and rock material of slope is usually the porous unsaturated or saturated media [1], and the stability of the slope has a great relationship with the strength of the soil and rock material, the reinforcement measures, the influences of water, earthquake, and external load [2, 3]. Due to the disturbance caused during the process of infrastructure construction, rainfall in the monsoon season, the fragile geology, and the variation of topography, land sliding happens more often. In engineering practice, the landslides which are a potential threat to human life and transportation security are constantly monitored to ascertain that people could be safely evacuated and the damage to property could be minimized or eliminated. Therefore, slope stability is still an issue which needs to be analysed in detail. In the past decades, researchers developed plentiful methods for slope stability analysis. The authors recommend that the slope stability and safety assessment should include as many factors as possible.
Generally, the factor of safety (FOS) calculation is an appropriate and convenient way to evaluate slope stability. As a traditional and well-established method, the limit equilibrium method (LEM) is widely used by geotechnical engineers, not only for its simplicity in estimating FOS with a fewer input data but also for its ability to accommodate complex geometries and variable soil and water pressure conditions [4]. For example, Liu et al. [5] developed an efficient direct Monte Carlo simulation (MCS) method, called adaptive MCS, for slope system reliability analysis based on LEM; Wang et al. [6] proposed a 3D slope stability analysis method based on the Pan’s maximum principle and under the framework of LEM. Despite all of its benefits, LEM has some pivotal deficiencies, such as neglecting the stress-strain relation of geomaterials and assuming numerous arbitrary slip surfaces for FOS calculation, which cannot represent a realistic failure mechanism. The limit analysis method (LAM), including upper and lower boundary approaches, has been used and improved by many researchers for slope stability analysis since Drucker and Prager [7] applied plasticity limit theorems in soil mechanics. Chen [8] has systematically reviewed the theory. Zhao et al. [9] proposed a new parameter back analysis method by combining the 2D/3D upper bound LAM and reliability theory to accurately determine the shear strength parameters for a 3D slope with a single failure surface. Adopting finite elements and linear programming, Sloan [10, 11] and Kim et al. [12] conducted both lower bound and upper bound analyses. Although LAM is rigorous in the sense for statically admissible stress fields for lower bound analysis and kinematically admissible velocity fields for upper bound analysis, its application in complicated real problems is still limited and it is seldom used for routine designing. The strength reduction method (SRM) has been used in slope stability analysis as early as 1975 by Zienkiewicz et al. [13]. Afterwards, many researchers [14–18] have applied it. The major advantage of SRM is that the critical slip surface is determined by reducing selected strength parameters under the gravity load until failure occurs. Therefore, it is convenient to use SRM for slope stability analysis in engineering practice.
In addition, the failure time of slope (FTS) is another important quantitative method based on monitoring displacement data to analyse the slope stability. Considering complicated boundary conditions, doubtful triggering failure mechanisms, and the heterogeneity of the geomaterials, it is extraordinarily difficult to predict FTS by adopting the rheological theory and rock fracture mechanics. An empirical approach is diffusely preferred, which is derived from the accelerating creep stage and fundamentally using displacement or deformation rates for the failure indicators. Most of these empirical equations are based on power and exponential laws [19]. According to the measured slope displacements versus time before failure taking the form of the tertiary creep curve, Terzaghi [20] showed the presence of a connection between creep and landslides. Saito [21, 22] concluded a method for FTS prediction by laboratory measurements of the strain rate during secondary creep using load-controlled triaxial tests. Zavodni and Broadbent [23] introduced an equation to predict two stages of creep and the time of failure. Fukuzono [24] through an experimental study of small-scale slope models found the logarithm of acceleration displacement to be proportional to the logarithm of the velocity of surface displacement of the slope. Later, Fukuzono [24, 25] proposed an inverse velocity method to predict the time of failure. Based on the research of Fukuzono [24, 25] and through theoretical analysis on power law creep under different loading conditions, Voight [26] deduced the correlation of velocity, time at the prediction moment, and the time of failure. Rearranging the equation of Saito [21] and comparing it with Fukuzono’s [24] inverse velocity method, Mufundirwa [27] invented a new method termed the SLO method for prediction of FTS.
In this paper, the authors analysed the stability of a slope that was reinforced twice, primarily with pile and retaining wall and secondly with porous steel-tube bored grouting. Using the geological data, the FOS in primary reinforced condition was calculated via the Universal Distinct Element Code (UDEC) equipped with the Shear Strength Reduction Method (SRM), which is capable of automatically locating the critical failure surface. Subsequently, FTS was predicted with the SLO method of Mufundirwa [27]. After that, the second reinforcement effect was evaluated by combining the monitoring data and the numerical simulation results. In addition, variations of displacement and FOS due to water influence are analysed. Through this procedure, a systematic approach for the slope safety evaluation and assurance is presented for future engineering practice reference.
2. Engineering Background
The 55 m wide slope is located in a coal mining area, which is also near the mileage of K18+170 Beijing-Zhuhai Highway in the south of China. The slope, consisting of three steps, was primarily reinforced with pile and retaining wall. The first step was reinforced by piles and a 5 m high retaining wall (Figure 1). Covered with a protective frame, the second one was 10 m high. The third step of 10 m height had a protective frame and a row of piles on the top. The natural slope is located above these structures. The second and the third steps were backfilled with soil in order to lower the ramp rate from 1 : 0.75 to 1 : 1.
2.1. Engineering Geology
To determine the engineering geological characteristics and the slip surface, four exploration holes were drilled along the sliding direction of the slope. Detailed exploration results showed that the lithological layers of the slope from top to bottom consisted of the overlying Quaternary soil (I artificial fill soil, II Quaternary alluvium layer 0-4.8 m), highly weathered limestone (4.8-6.5 m), peat layer (6.5-10.0 m), moderately weathered limestone (10.0-13.75 m), carbonaceous shale (13.8-21.0 m), followed by thin and slightly weathered sandstone and alternating layers of carbonaceous shale (21.0-30.5 m).
On the basis of the geological mechanics analysis method, the engineers supposed the potential sliding surface to occur along the interface of the peat layer. An engineering geologic profile is shown in Figure 1.
2.2. Damage to Reinforced Structure
In October 2003, a curved crack of about 35 m was found behind the piles of the third step (Figure 2(a)). The crack after three years of its formation is shown in Figure 2(b). In addition to this crack, another 1.5 cm wide transverse crack appeared in the middle of the third protective frame, which grew to a width of 30 cm in January 2006 to such an extent that the protective surface of the mortar rubble loosened and water leaked through it.
(a)
(b)
2.3. Displacement of Slope
Being aware of the considerable change in appearance, the engineers deduced that the displacement rate of the slope was large. Due to the reason that this highway is a transportation artery in the south of China and a coal mining area is located behind the slope, in case of slope failure, there would have been a serious threat to human life and transportation security. A number of actions were taken to stabilize it, including safety monitoring, stability analysis, rereinforcement, and effect evaluation. The displacement versus time curves and displacement at depth as shown in Figure 3 were obtained through monitoring displacement at the surface and at depth. A tremendous increase in displacement occurred due to heavy rains during May 2006 to July 2006. Maximum surface displacement was noted to be about 80 mm, while the distinct displacement change at depth ranged from 6.5 m to 8.5 m along the peat layer. The potential slip surface is shown in Figure 1.
(a)
(b)
3. Stability Analysis Based on the Discrete Element Method with SRM
3.1. Constitutive Models of Material and Contacts
The Mohr-Coulomb failure criterion (equation (1)) is used to describe the failure of block materials. The constitutive model of contact is the Coulomb slip model with residual strength (equations (2) and (3), Figure 4) which is used to control the response of contact in the normal and shear direction [28]. where is the shear strength, is the cohesive strength, is the normal stress, and is the friction angle of the block. where and are the normal stress and shear stress, respectively, and are the normal stiffness and shear stiffness, respectively, and are the normal displacement and shear displacement, respectively, is the contact overlap tolerance, and are tensile strength and overlap tolerance strength, respectively, is the shear strength of contact, and are the cohesive strength and residual cohesive strength, respectively, and are the friction angle and residual friction angle, respectively, and is the incremental shear displacement.
(a)
(b)
3.2. Theory Background of SRM in UDEC
The UDEC is a two-dimensional numerical program based on the distinct element method for the discontinuum modelling [28]. In the strength reduction method, the selected strength properties are reduced until failure occurs and the FOS is calculated. The method is commonly applied with the Mohr-Coulomb failure criterion. In this case, the FOS is defined according to the following equations: where is the trial value of cohesion, is the trial value of the factor, is the cohesion, is the trial value of friction angle, and is the friction angle.
The original strength parameters of rock and soil (such as and ) are evaluated by tests during survey processes. When conducting the strength reduction method for the numerical model, the is chosen and the trail values of the strength parameters are obtained by equations (4 and 5). A series of simulations are conducted using the trial values of the factorto reduce the cohesionand friction angleuntil slope failure occurs (note that if the slope is initially unstable,andare increased until the limiting condition is found). One technique to find the strength values that correspond to the onset of failure is to monotonically reduce (or increase) the strength parameters in small increments until a failure state is found.
3.3. Basic Model
The lateral and vertical dimensions of the numerical model are , respectively. The mesh consists of deformable triangle zones. The lower boundary is assumed to be fixed, while the vertical boundaries at the left- and right-hand sides are assumed to be on rollers to allow movement of soil/rock layers.
The pile and retaining wall are represented by a block with the corresponding material strength, and the cap rigidity of the pile is free. The slope that was reinforced by two rows of piles and retaining wall is a 3D problem, as shown in Figure 5. If the 3D simulation method is used, the numerical model will be huge and the calculation time will increase several times than that of the 2D simulation method. However, the UDEC is a 2D software which simulates the unit thickness of the slope. For this reason, some equivalent conversion for piles needs to be done.
As shown in Figure 5, the sliding force in 2D became one unit of the total sliding force in 3D. Because the distribution of the piles is discontinuous, the bending stiffness of the pile should also be equivalently transformed, which means that the elastic modulus and moment of inertia should be changed from 3D to 2D conditions. Donovan et al. [29] put forward a simple and convenient method to convert material parameters of the same distance distribution structure obeying the equivalent conversion rule in the finite element analysis. Likewise, Zhang et al. [30] proposed the conversion method for pile, which is used in this paper.
The length of the pile (1) is 28 m in the upper part of the slope, and the length of the pile (2) is 12 m in the toe of the slope, respectively. The section of the pile (1) is with a 6.5 m distance between each of the two piles, and the section of the pile (2) is with a 6.0 m distance between the piles. As one part of the pile (2) was inserted into a continuous retaining wall, it was difficult to distinguish between them. So only the pile (1) and part of the pile (2) are converted, and the material parameters of the inserted part of the pile (2) are considered to be the same as that of the retaining wall for simplicity. Referring to the method of Zhang et al. [30], the width of the pile (1) and pile (2) are taken to be 3 m and 2 m, respectively, and the thickness of these two is taken as 1 m in the numerical model. The corresponding parameters of the pile were calculated as follows:
For pile (1)
For the part of the pile (2) which was not inserted into the retaining wall
Poisson’s ratio of the pile is 0.2, so the bulk and shear modulus can be obtained using the following equations. where is the equivalently transformed elastic modulus of the pile (1) in 2D condition, is the elastic modulus of the pile (1) and pile (2) in 3D condition, is the moment of inertia of the pile (1) in 3D condition, is the moment of inertia of the pile (1) in 2D condition, is the horizontal distance between two piles (1), is the bulk modulus, is the elastic modulus, is the Poisson’s ratio of the pile (1) and pile (2), and is the shear modulus.
Table 1 summarizes the material parameters used in the analyses. The Mohr–Coulomb constitutive model was used for each layer. The soil/rock properties used in this study were originally obtained from geological exploration.
In order to get the joint parameters between different layers. The in situ shear tests should be carried out or the intact samples with joint should be taken from the slope and be tested on the laboratory. However, the above tests were not conducted due to the difficulty of in situ shear tests and taking the intact sample with a joint between layers during that time. From a conservative point of view, the friction angle and cohesion of the joint between two layers are assumed to be the same as that of the weaker layer, and the joint tensile strength is not considered. The joint normal and shear stiffnesses were calculated using the following formula [28] and are equal to one another. where is the smallest width of an adjoining zone in the normal direction.
3.4. The Simulation Results
3.4.1. FOS Calculation
In order to reduce the left boundary effects on the inclined soil and rock layers, we have added the rectangle block on the left boundary. The rectangle block has several layers and the values of material parameters are the same with the inclined soil and rock layers. The failure state of velocity arrows identified using the SRM analysis shows that the minimum FOS is 1.06 (Figure 6). Obviously, the velocity at the outlines of the first layer and pile head is larger than that of the other layers. The slope surface movement is caused mainly by the first three layers because the modulus and strength parameters of these three highly weathered layers are small.
3.4.2. Displacement
Displacement in the lateral direction for the slope is depicted in Figure 7. It can be confirmed that there is large surface displacement in the second and third step, with maximum values of about 0.3 m and 0.35 m, respectively. The displacement is large in the cap of pile, in the upper part of the slope, which means that the bending stiffness of the pile was not adequate. Therefore, this phenomenon is logical in a way that the protective surface of mortar rubble loosens in these steps. Between the second step and the pile, there is a displacement difference which is shown behind the pile in Figure 7. This simulation result can explain why the crack behind the pile grew (shown in Figure 2).
3.4.3. Slip Surface
According to the maximum shear strain contours shown in Figure 8, the layer with a floating range of the shear strain from 0.04 to 0.1 represented by yellow and red is concluded to be the slip surface. It is the peat layer and is the weakest layer in the slope. This simulation result adheres to the deep displacement monitoring data shown in Figure 3. The failure mechanism of this slope is that the first two layers slid along the weakest peat layer whose strength was greatly influenced by water.
From the simulation results exhibiting small FOS and monitoring data confirming a large displacement rate, the authors deduce that this slope is not safe and the primary reinforcement process with pile and retaining wall has not reached the anticipant goal. Referring to Figures 7 and 8, it can be concluded that if the slope movement continues to increase, the pile will become unstable (bending stiffness does not satisfy the requirement) and the slope will fail.
4. Failure Time Prediction without the Second Reinforcement
In this section, the authors analyse the monitoring surface displacement data shown in Figure 3 and adopt the SLO method to predict the failure time of slope (FTS) since the FOS is very small.
4.1. Theory and Background of Predicting Method
As early as 1950, Terzaghi [20] concluded that there are some connections between the creep of rock or soil mass and landsliding. After that, numerous researches have been carried out to find out the slope failure characteristics corresponding to the standard creep curve. Some researchers have verified that a slope will become unstable if the slope’s sliding velocity accelerates or displacement versus time curve is similar to the tertiary creep curve. There are several laws describing the displacement versus time relationships during tertiary creep which can be used to forecast time of slope failure. These methods were proposed by Saito [21, 22], Fukuzono [24, 25], Voight [26], Fukui and Okubo [31], and Mufundirwa [27], as shown in Table 2.
In Table 2, where is the strain rate or velocity, is the strain or displacement, , , , and are constants, is the failure time, is the time, is the displacement, and is the velocity.
The SLO method used in this paper for failure time prediction was proposed by Mufundirwa [27]. After substituting displacement instead of strain and differentiating Fukui and Okubo’s equation [31], equation (9) is deduced.
By rearranging equation (9), equation (10) is obtained. where is the failure time evaluated as the slope of versus curve for equation (10); this new method is termed the SLO method.
Utilizing the monitoring data, the velocity is calculated using equation (11), which can filter the measured data in order to smoothen the short-term deformation deviations that can be insignificant or may cause “false” results [32]. where are the computed velocity points, is the monitoring displacement, is the corresponding time of monitoring displacement, is the time at the instant of prediction, and is the displacement at the instant of prediction. In this equation, the sampling value is selected to yield a positive velocity only.
4.2. Predicting Results
Based on the monitoring surface displacement of four sites mentioned in Figure 3, the prediction procedure is displayed as follows. Choosing equation (10) to calculate all data which are smoothened using equation (11), the predicted results are shown in Table 3 ( is the correlation ratio).
The monitoring data in Figure 3 include the recorded surface displacement from 1st of January 2006 to 16th of July 2006, a total of 196 days. However, the predicted results in Table 3 indicate the slope failure time to range from 200 to 212 days, which means that if this tertiary-like slope movement continues (without any further strengthening), the landsliding will occur between 22nd of July 2006 to 1st of August 2006. It can be concluded that the for sampling value is bigger than that for . However, there is a difference of just a few days.
5. Second Reinforcement
5.1. Reinforcement Scheme
In order to mitigate the landsliding, the second reinforcement project was launched soon after the cracking in the primary reinforced structure and after a large displacement rate of the slope was discovered. The second reinforcement project included porous steel-tubes bored grouting, with the grouting area covering the first three weak layers from the top (Figure 9).
The diameter and separation distance of steel-tubes are 90 mm and 1.5 m, respectively. The diameter of the drilling hole for installing the steel-tube is 250 mm, and the end of the drilling hole is in the bottom of moderately weathered limestone. The cement pastes are injected into Ø22 mm PVC pipe, first flow through interspace between the steel-tube and drilling hole, then entering into the soil and rock mass voids, with the grouting pressure controlled less than 2.5 MPa. After hardening of the cement pastes, the steel-tube and the surrounding geomaterials binded so strongly together that a small pile is produced.
5.2. Effectiveness Evaluation
5.2.1. Final Monitoring Results
Figure 10(a) depicts the final monitoring surface displacements. Because the soil pore pressure increased and the strength of the saturated peat layer decreased during the monsoon rains from the end of April 2006 to the end of July 2007, the surface displacement rate monitored for four sites is at a relatively high value with an average of 0.18-0.41 mm/d. As a result of the effect of the first porous steel-tube bored grouting, the tertiary-like slope surface movement displacement rate decreased during September 2006 and March 2007. Due to the monsoon in April 2007, the surface displacement rate increased again with an average of 0.12-0.22 mm/d; that is why the second porous steel-tube bored grouting was conducted soon after October 2007. Then onwards, the slope displacement grew slowly with an average rate of 0.02-0.03 mm/d. The displacement rate in 2007 was smaller than in 2006, which shows that the deformation of slope began to converge. While in the same year, the displacement rate in the rainy season was obviously bigger than in other seasons, which demonstrates that the water has a huge influence on slope stability.
(a)
(b)
Figure 10(b) displays the final displacement monitored at depth. It can be concluded that the slope’s movement occurs mainly in the above three layers at a depth of 0 m to 8.5 m (Quaternary soil, highly weathered limestone, and peat layer) and that the slope is stable at depth (displacement rate close to zero).
5.2.2. Displacement Simulation considering the Water Influence
In consideration of the quick growing displacement during the monsoon rains, the influence of water has been analysed. There are two competing effects on the slope deformation related to the rise of the water level: first, the increase in pore pressure will generate a decrease in effective stress leading to shear strength degradation; second, water flowing in contacts and entering into the material pores will increase in density and cause settlement. In the numerical model, we proceed to simulate water can flow in contacts, pore pressure in block material, and density increases with water table rising.
The flow rate in contact is given by [28, 33] where is the permeability factor in contact (whose theoretical value is ), is the dynamic viscosity of the fluid, is the contact hydraulic aperture, and is the length assigned to the contact. where and are the pressure in domain 1 and domain 2, respectively, is the fluid density, is the acceleration of gravity, and and are the -coordinates of the domain centers.
The contact hydraulic aperture is given by where is the contact aperture at zero normal stress and is the contact normal displacement.
For block material, shear strength is described as follows [34] where is the pore pressure in the block material.
After enforcing with grouting, the strength parameters of the first three layers increased. As a conservative analysis, only the steel-tube and cement composite piles (the diameter is 250 mm) are added in the numerical model and the parameters of the first three layers are invariable. The original parameters and equivalent conversion parameters which consider the 2D to 3D effect of the steel-tube and cement composite piles (refer to equations (6), (8), and (9)) are listed in Table 4.
The friction angle and cohesion of the joint between the two layers are assumed to be the same as that of the weaker layer. The joint normal and shear stiffnesses were calculated using equation (8). The failure of the block materials is controlled by equation (15). In the following simulations, we assign , , and [28].
The displacement simulation results taking account into water table variation influence are illustrated in Figure 11. Figures 11(a) and 11(b) show different water tables in the slope reinforced with pile and retaining wall and then by porous steel-tube bored grouting. Water tables 1-4 represent the water level in 0 m, in the dry season (2006/1/1), in medium (2006/5/16), and high level (2006/8/16), respectively. Especially, the increase of water level from water tables 3 to 4 presents the whole effects of the monsoon rains. Figure 11(c) shows the displacement in two conditions corresponding to Figure 11(a) and 11(b). The exact displacement is bigger than the monitoring data because the deformation of slope began long before 2006, but the monitoring data initiate in 2006. That is why our simulated displacement is bigger than the monitoring displacement. In addition, the displacement in simulation 1 is obviously larger than in simulation 2. Assumed displacement in water table 2 (2006/1/1) as the initial displacement (0 mm), the increment is shown in Figure 11(d). It can be deduced that the maximum displacement increment occurs when the slope is in simulation 1 (primary reinforcement), the minimum occurs when the slope is in simulation 2 (all porous steel-tubes bored grouting have been finished), and the monitoring data is between the maximum and the minimum.
(a)
(b)
(c)
(d)
5.2.3. New FOS Calculation
According to the different water levels in Figure 11(b), the new FOS values are calculated (Figure 12). As shown in Figure 12, the FOS is decreased with the increase of water level. The FOS is 1.17 during the highest water level, 1.39 during the dry season, and 1.42 without considering water influence, which is higher as compared to the FOS (1.06) grouting reinforced before. All the new FOS values are conservative calculation results, because only the steel-tube and cement composite piles are taken into account in the numerical model and the improvement for parameters of the first three layers by grouting is neglected.
The minimum new FOS and velocity vectors during the highest water level are shown in Figure 13.
Based on the comprehensive referencing, monitoring, and numerical simulation results, as shown in Figure 10, the increment speed of surface displacement and deep displacement of slope both became slow; as shown in Figure 12, all values of FOS for different water tables are bigger than 1.0. Therefore, the authors conclude that the slope became stable and safe after the second reinforcement project.
6. Conclusion
The stability of a slope was analysed at a primary (pile and retaining wall) and secondary (grouting) reinforcement stage. For the primary reinforced stage, the FOS is 1.06 calculated by the Shear Strength Reduction Method embedded in UDEC. In addition, the lateral displacement difference around the piles was chosen to explain the phenomenon that the crack emerged behind the piles of the third step. The shear strain contours demonstrate the sliding surface location which is the maximum shear strain area. In fact, it was the weakest layer of the slope. The failure mechanism of this slope is such that the first two layers slid along the weakest layer, i.e., the peat layer. Alongside, the failure time of slope was predicted applying the SLO method while using the monitoring data of the tertiary-like slope movement. The results show that without the second reinforcement scheme, the slope would have failed between 22nd July 2006 and 1st August 2006. The small FOS, failure time prediction results, monitoring data, and site damage phenomenon all show that the slope is unstable at this primary stage. For the second reinforced stage, the porous steel-tube bored grouting scheme was applied twice. When the parameters of the first three layers were increased by strengthening through grouting, the updated FOS increased to 1.42. Interpreting the final monitoring data and numerical simulations, two conclusions are deduced. One that the displacement rate in the rainy season was obviously higher than in other seasons and the water table has a great influence on the value of FOS. Secondly, the deformation of the slope began to converge and the slope is considered stable keeping in view the new FOS and the supporting monitoring data.
The SRM is a good option for obtaining FOS for this kind of slopes whose failure mechanism is such that several upper layers slide along the weakest layer. The SLO method is an innovative way of quantitative analysis of monitoring displacement data for slope failure time prediction. The authors highly recommend the system method considering FOS, failure time prediction results, monitoring data, and the damage phenomenon to conduct slope stability analysis and reinforcement effect evaluation. Future research is expected to compare more engineering practices for further validating this system method.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Authors’ Contributions
Wei Chen is responsible for the methodology, conceptualization, and writing of the original draft. Dongbai Li helped in data curation. Ting Ma contributed to the supervision and resources. Helin Fu contributed resources. Yanpeng Du wrote, reviewed, and edited the manuscript.
Acknowledgments
This work was supported by the Research Foundation of Education Bureau of Hunan Province (19K100).