Chamber roof deformation prediction and analysis of underground mining using experimental design methodologies

The deformation prediction of the chamber roof is crucial in underground mining. Combined with Flac3D numerical simulation, the experimental design methodologies, including single-factor test (SFT), Plackett–Burman design (PBD), steepest ascent design, and response surface methodology, were used to evaluate the effect of multiple variables on the chamber roof deformation. Firstly, eight factors that affected the vertical displacement (Ds) of the chamber roof were selected, and the sensitive interval of every factor was obtained through SFT. Then, four factors that significantly affect the results were screened by PBD: cohesion (Co), stope length (Ls), stope width (Ws), and internal friction angle (fr). Twenty-nine groups of response surface schemes with 4 factors and 3 levels satisfying the Box–Behnken design (BBD) were simulated. Through the result analysis of variance (ANOVA) and sensitivity, the influence order of each factor on Ds can be determined: Ws >Ls >fr > Co > interaction between Co and fr > interaction between Ls and Ws. Finally, using the prediction model, the roof deformation of 0#N stope of a lead–zinc mine was predicted and the error was analyzed. The relative errors between the prediction value and the numerical simulation value, the measured value are 0.7% and 4.3%, respectively, which indicates that the prediction model is reasonable and has a certain reference value for mine safety.


Introduction
The chamber deformation of the underground stope is the focus of rock mechanics workers and mine managers, and it affects the production capacity of mines and the safety of personnel directly (Abdellah et al. 2014;Brown 2012). Although there are many factors to describe the chamber roof deformation, the displacement is most intuitive, so displacement prediction is one of the most critical challenges in mining, and it is of profound significance to mining safety. The mechanized panel slicing backfilling mining method has been widely used in China because of its high mining efficiency and exemplary safety. However, when this method is used in mine inclined orebody, pillar spalling, roof fall, and rockburst may occur in the chamber due to the unreasonable structure of the stope. As a result, it is necessary to establish a reasonable plan to predict stope deformation.
The factors affecting roof stability are very complex. The uncertainty of these factors will become stronger when they reach a certain underground depth, which seriously threatens safe mining and puts higher demands on the prediction and evaluation method (Qiu et al. 2018). Scholars have published many valuable research findings in this area in recent decades. In conclusion, the commonly used methods for stability analysis and deformation prediction of underground engineering are as follows: analytical method, empirical method, numerical modeling, in situ measurement, and model test method (Heidarzadeh et al. 2018;Zhang et al. 2022). The analytical method has good applicability to simple engineering problems. However, it has limited capacity to deal with sophisticated mining issues. Standard empirical methods (such as Barton Q classification method (Barton 2002), Mathews stability graph method (Mathews et al. 1980;Potvin 1988), extended Mathews stability graph (Mawdesley et al. 2001;Mawdesley 2004)) also have certain limitations, such as oversimplification of irregular geometries and failure to consider the effects of pillars in the chamber, thus the functionalities of these methods are restricted. In recent years, some scholars have applied artificial intelligence technology (AI) to study the pillars strength (Mohanto and Deb 2020;Wattimena et al. 2013;Zhou et al. 2015) and the stability of the stopes (Goh et al. 2017;Qi et al. 2018;Yu et al. 2020). This method is effective when analyzing a single object, but it has limited ability to study the chamber with pillars, and the collection of training samples is complicated. In addition, in situ measurements and indoor model experiments are time-consuming and costly. On the other hand, numerical methods can circumvent these limitations because they can adopt different failure criteria corresponding to various failure mechanisms and intricate mining problems can be solved by building accurate numerical models according to actual stope conditions Shnorhokian et al. 2015). Therefore, they have been widely employed in the research of underground excavation stability (Liu et al. 2018;Vallejos et al. 2018).
The chamber roof stability of the inclined orebody is affected by rock mass properties, in situ stress, blasting disturbance, and stope structure parameters. Among them, the stope structure parameters and rock mass properties are very essential (Heidarzadeh et al. 2018). In general, the influence of stope structure parameters and rock mass properties on the stope chamber stability with pillars is usually studied by empirical stability analyses and experiments. For instance, Milne et al. (2004) proposed that it was preferable to use the radius factor (RF) to determine the actual stope surface in the stability graph. Zhou et al. (2019) proposed a pillar risk assessment model combined with pillar stability, which quantified the impact of a single pillar failure on the danger of the whole stope. Li et al. (2020) used the model test method to simulate the pillar recovery, and they analyzed the instantaneous deformation responses and long-term settlements of the stope during the mining of the pillar. As a standard stability analysis method, numerical modeling also has been widely used to evaluate the effect of stope structure parameters and rock mass properties on the chamber roof stability. For instance, some scholars studied the effects of geometrical parameters (Heidarzadeh et al. 2019) and in situ stress uncertainties ) on the stope stability using Flac3D. Some others modified the rock stress factor in the stability graph method by combining empirical methods and numerical modeling (Jia et al. 2020), and investigated the damage behavior of rock mass using Flac3D (Yilmaz and Unlu 2013). However, the specific impacts of these factors on the chamber deformation have yet to be clarified. Moreover, while the methods discussed in this study are better than traditional research methods, they are incapable of solving the challenge of assessing multivariate systems.
The most effective methods to analyze the influence of stope structure characteristics and rock mass properties on stopes roof stability is probabilistic method and numerical modeling (Idris et al. 2011). Among the several probabilistic approaches used in stope stability concerns, experimental design methodologies like response surface methodology (RSM) have been widely employed to analyze the stability of underground excavations (Lü and Low 2011;Mollon et al. 2009;Pytel and Świtoń 2014). RSM can assess individual and interactive effects of multiple influencing factors on the stability evaluation index and establish a prediction and analysis model for the stope roof stability. The influence of each factor on stope stability and the interaction between each factor can be directly evaluated through the generated mathematical model. The current research aims to give a complete assessment of the influence of critical factors on the chamber roof stability of inclined orebody, and the main implementation steps of the experimental design methodologies are as follows: 1. The results of single-factor test (SFT) are used to select the suitable sensitive interval for each factor. 2. Plackett-Burman design (PBD) is used to screen out the significant factors. 3. The steepest ascent design (SAD) is used to determine the center point of the response surface, and RSM combined with Flac3D is employed to assess the individual and interactive effects of the significant factors. 4. A polynomial regression model of chamber roof displacement prediction is set up.
The prediction model, which can reflect the relationship between the vertical displacement (D s ) of the chamber roof and various factors, is then applied in the deformation prediction of underground stopes and achieved good results.

Site description
After more than 50 years of mining, a lead-zinc mine in Guangdong Province, China, has a large number of ore bodies with a dip angle between 35° and 55° currently, and the ore grade is high, a mechanized panel slicing backfilling mining method was adjusted to be used to mine such inclined orebody at present. The thickness of the inclined orebody in this mine is generally medium thick, so the stopes have changed from the vertical direction to along the direction of the orebody. The width of one slicing is typically similar to the stope thickness, up to 15-22 m, to control dilution and loss. The length of the stope becomes longer, up to 40 m approximately, to improve the mining efficiency (Zhang et al. 2022). Therefore, the exposed area of the chamber roof is comparatively large. A pillar is left in the chamber to improve mining safety, and the stope layout diagram is shown in Fig. 1 (taking the stope chamber of the 55° orebody as an example).
The entire zone of the stopes can be divided into the cutting groove area and lateral blasting area, and the lateral blasting area is blasted after the cutting groove area. During in situ blasting mining, one high layer (7-9 m) will be mined once, and the orebody will be recovered at the same time as the ore pillar. The occurrence state of the inclined orebody stope is complex, and the influence of various factors on the chamber stability is not apparent totally. Hence, it is urgent to establish a reasonable and practical stope stability prediction and evaluation method.

Numerical modeling
Selected the 0#N stope as the research object, and Flac3D was used to analyze the chamber roof stability with a pillar comprehensively. Figure 2 shows the built model of a 55° dip angle.  The foot wall and hanging wall of the orebody are all wall rock. According to the principle of Saint-Venant, the surrounding area that is significantly affected is 3-5 times the size of the excavation , and the size of the entire model is X × Y × Z = 100 m × 100 m × 70 m. Valid and stable numerical results were obtained by testing with different mesh sizes. A finer mesh is created around the stope chamber and pillar for higher accuracy, whereas the mesh becomes coarser for the wall rock. The mesh size of the numerical model is about 0.4-2.5 m, and the number of meshes is about 150,000. This study uses the Mohr-Coulomb constitutive model for Flac3D numerical simulation.
As shown in Fig. 3, we set four monitoring points: 1#, 2#, 3#, and 4#. The monitoring points utilized in this paper are also these four monitoring points, and the chamber displacement and the plastic zone were observed using these points.
Two types of typical rock samples were selected for the physical-mechanical parameters tests, and the parameters of rock mass required for numerical simulation are listed in Table 1, where ρ is the density, K is the bulk modulus, and G is the shear modulus.
The rock mass is already subject to in situ stress before an underground excavation. These stresses are due to the weight of the overburdened, tectonic activity or are related to the local topography. The underground stability will be affected by in situ stress. The stope studied in this paper is generally concentrated in the deep underground, and the actual measured gravity stress load is very close to the vertical component of in situ stress, which can be determined according to the occurrence depth of the orebody. The ratio of the maximum horizontal components to the vertical component is between 1.16 and 1.49. The direction of the maximum horizontal components is parallel to the strike of the orebody . Assuming that the rock mass is uniform and continuous, the initial stress fields can be calculated by the following formula:   where γ is the bulk density of the rock, and H is the depth of the rock mass. The boundary conditions of this numerical simulation adopt displacement constraints. As illustrated in Fig. 2, the four sides constrain the horizontal displacement, and the lower surface constrains the vertical displacement, the upper surface applies gravity stress load, and the four sides apply tectonic stress loads.

Experimental design
As mentioned in the last paragraph of the introduction, there are four basic steps in the experiment design. Firstly, the sensitive interval of each factor is obtained by SFT. Then, PBD is used to screen out the factors that have a significant impact on the results. Next, SAD is used to determine the center point of the response surface. Finally, a response surface scheme is created. Thus, we can determine the relationship between each influencing factor and the chamber roof displacement.

Single-factor test design (SFT)
The main parameters related to ore strength are: density (ρ), bulk modulus (K), cohesion (C o ), internal friction angle (f r ), shear modulus (G), and tension strength (T). In slope engineering, shear failure dominates, and the factors that have a significant effect on the safety are C o and G, while for underground stopes, C o , f r , and T are the main parameters affecting ore strength and stope safety (Chen et al. 2017;Shi et al. 2011), so other parameters remain unchanged and only change these three factors when we consider the effect of ore strength on stope stability.
In addition, we selected orebody dip angle (α) and four factors closely related to the stope structure parameters: pillar height (H p ), pillar width (W p ), stope length (L s ), and stope width (W s ). Therefore, a total of 8 factors that have a more significant effect on the chamber roof were considered in SFT. In actual production, the height required for rock drilling in the chamber and the working height of the workers determined the pillar height, so H p is generally 3.0-4.2 m. As mentioned earlier, W s can reach 15-22 m to improve mining efficiency, α of the inclined orebody ranges from 35° to 55°, the three factors mentioned above are greatly restricted by the site conditions, so their sensitive range is also set within this range. The following fixed α, H p , and W s at 55°, 4.2 m, and 15 m, respectively. SFT was performed on W p , L s , C o , f r , and T, and determined the sensitivity ranges of these factors. As illustrated in Table 2, fixed four factors, and the other factor is changed for the experiment. The influence factors (X) and 10%X were used to determine the error to prevent the MPa and monitored the vertical displacement of the roof by changing f r (the range of f r is 20° to 60°, and the simulation was performed once every 10° increments), the displacement of each monitoring point was collected by the levels f r and f r ± 1°. After that, change T (initial W p = 3 m, L s =45 m, C o = 5 MPa, f r = 40°) from 1 to 9 MPa by uniformly increasing the strength 2 MPa, and gather the displacement date of 1#, 2#,3#, 4# points by the levels T-0.2 MPa, T, T + 0.2 MPa. And then, initial L s =45 m, C o = 5 MPa, f r = 40°, T = 3 MPa, and monitored the vertical displacement of the roof by changing W p (the range of W p is 2 m to 6 m, and the simulation is performed once every 1 m increment), the displacement of each monitoring point was collected by the levels W p and W p ± 1 m. Finally, change L s (initial W p = 3 m, C o = 5 MPa, f r = 40°, T = 3 MPa) from 25 to 65 m by uniformly increasing the 10 m, and gather the displacement date of 1#, 2#,3#, 4# points by the levels L s -10 m, L s , L s + 10 m.

Plackett-Burman Design (PBD)
It is necessary to screen all the factors and identify the significant factors to reduce the number of tests after obtaining the sensitive interval of each factors. PBD, which has been used in this paper, can also be called screening experimental design, mainly used to screen important factors, reduce workload, and save costs (Vanaja and Rani 2007;Zhao et al. 2010). This design does not take into consider the interaction between factors but only determines the factors that significantly impact the results to screen and optimize the experimental factors. Here, we created a PBD with 12 tests, and 2 levels were taken for the 8 factors of the design. The selected factor's high level is set to be 1.0-2.0 times the low level, and other test conditions such as in situ stress and model mesh sizes are fixed.

Steepest ascent design (SAD)
After founding the significant factors through the PBD, it is necessary to change the significant factors toward the maximum direction of the response value at the same time to approach the maximum response interval and use SAD to obtain the center point of the response surface, so that an effective response surface equation can be built. The factors that have no significant effect on roof displacement are fixed in SAD. In contrast, the significant influencing factors take the change direction of the test results as the raised path so that the optimal value area can be approached as quickly as possible ).

Box-Behnken Design (BBD)
RSM is used to evaluate the effects of selected parameters associated with orebody properties parameters, orebody dip angle, and stope structure parameters on the obtained chamber roof displacement. It can be used to design experiments, optimize mathematical modeling, and get the optimal level value of each factor (Mason et al. 2003). The fundamental thought 1 3 of RSM is to obtain a large amount of data through tests, estimate the implicit statistical relationship between the input variables and the response factor. Box-Behnken Design (BBD) is a commonly used RSM, which can be used to evaluate the nonlinear relationship between indicators and factors. Compared with the RSM, although the orthogonal method is simpler to calculate, it cannot find the exact optimal combination and the optimal value of the response variable in a given area. The best combination obtained by the orthogonal experiment can only be limited to the level determined by oneself before rather than the best combination necessarily. The RSM considers the random error of the experiment, and the obtained prediction model is continuous, which can optimize and evaluate the level and interaction of each factor.
Two main steps are taken when using RSM: create a mathematical equation to identify the input parameters and interactions that have a substantial impact on the result; optimize each critical parameter to obtain the desired response value (Ahmadi et al. 2014). A polynomial regression model with quadratic coefficients as formula (4) can be obtained by fitting the experimental data: where D p stands for the predicted displacement, β 0 , β i , β ij , β ii are constant coefficient, linear coefficients, interaction coefficients, and quadratic coefficients, respectively, and ε is considered as the random error. The significant input variables are denoted by x i and x j , and they can be converted by formula (5): where X i is coded of the input variable, x 0 represents the actual value of the input variable at the center point, and Δx is the steps between independent variables. The entire experimental design process is shown in Fig. 4.

SFT results
The difference between the vertical displacement of the four monitoring points 1#, 2#, 3#, and 4# on the chamber roof is almost nonexistent, so the average value of the vertical displacement Ds is taken to replace the displacement of those four points. The fluctuation of the vertical displacement of each scheme at 10%X is investigated using error bars, and Fig. 5 depicts the results. Figure 5a illustrates a negative connection between C o and D s . The impact process can be broken down into 3 stages. The shear failure region (C o < 5 MPa) is the first stage, which can be produced by low ore strength or external stresses loads, the vertical displacement of the roof is close to or even more than 30 mm, and the possibility of failure is very high (Chen et al. 2016), as shown in Figs. 6a and 7a (C o = 2 MPa), the entire pillar is in shear failure, the roof and both sides of the chamber are subjected to a large shear force, the pillar undergoes compression deformation under the external loads, and the side walls also tend to collapse to the interior of the chamber. The chamber roof has a downward vertical displacement of about 33 mm, and the bottom plate has an upward bulge of about 25 mm. Stage II is a sensitive area (5 MPa < C o < 11 MPa), in which the vertical displacement decreases significantly with C o increasing. Stage III is the stable area (C o > 11 MPa), in which the vertical displacement does not change dramatically with the increase of C o , the error bars show the variations in displacement.
A declining and nonlinear relation between f r and D s is displayed in Fig. 5b, and it can be separated into 3 stages. Stage I is the failure area (f r < 35°), which may be caused by the low ore strength and occurs on the pillar and both sides of the chamber, the entire pillar and the chamber roof are all in shear failure, as shown in Fig. 6b (f r = 20°). As shown in Fig. 7b (f r = 20°), the chamber roof has a downward vertical displacement of more than 30 mm, and the safety condition of the chamber is poor. As seen by the error bars, Stage II is a sensitive area (35° < f r < 55°) that fluctuates with variations in f r . Stage III is the stable area (f r > 55°), in which the vertical displacement does not change significantly with the increase of f r .
As illustrated in Fig. 5c, D s satisfies the nonlinear declining relation for different T, 3 stages exist: Stage I is the failure area (T < 3 MPa), in which the low ore strength may cause. Stage II is a sensitive area (3 MPa < T < 5.5 MPa), the error bars show that the change in D s is noticeable in this location. Stage III is the stable area (T > 5.5 MPa), in which the vertical displacement does not change significantly with the increase of T. Figure 5d illustrates a negative connection between W p and D s , 3 stages exist: Stage I is the failure area (W p < 3 m), which may be caused by the too narrow pillar and the large load in the stud. Stage II is a sensitive area (3 m < W p < 4.5 m). Stage III is the stable area (W p > 4.5 m), in which the vertical displacement does not change substantially with the increase of W p . Figure 5e illustrates a positive connection between L s and D s , 3 stages exist: Stage I is the stable area (L s < 35 m), in which the vertical displacement is small, and the stope safety is good. Stage II is a sensitive area (35 m < L s < 50 m). Stage III is the failure area (L s > 50 m), which the too-long stope may cause. In this paper, the mathematical function fitting (Table 3) is used to study the nonlinear relationship between D s and C o , f r , T, W p , L s , and it can be seen from the table below that the fitting coefficients are all above 0.9, and the fitting degree is reasonable.

PBD simulation results
Twelve groups of specific tests were implemented according to the scheme of PBD, and the vertical displacements of the 1-4# monitoring points D s were averaged and listed in Table 4. The significance of the factors was analyzed.
The findings in Table 4 showed a broad range of D s from 18.21 to 41.09 mm in the 12 tests, and this variety revealed that optimizing influencing factors were essential for enhancing chamber roof stability. The above test results are analyzed; the estimated effect, regression coefficient, and corresponding t and p values are given in Table 5. The calculated standard deviation of the calculation is 1.16029. The sum of prediction errors squares   (Table 5) show that C o , L s , W s , f r were the significant parameter, at the same time, H p , W p , α, T were considered nonsignificant, they were fixed as 3.6 m, 4.5 m, 55°, 1.2 MPa, respectively, in the following path of SAD and BBD tests.

SAD simulation results
By analyzing the PBD results of Table 5, it is found that L s and W s positively affect the D s . In contrast, C o , f r have a negative effect on the D s , the path of steepest ascent was moved along the direction in which C o , f r decreased and L s , W s increased. Five tests were designed as presented in Table 6.
It is noted that the D s of the five SAD tests has been decreasing from Table 6. The analysis shows that increasing the C o , f r , and decreasing the L s , W s of the stope will lead to smaller D s and better chamber safety. However, mining efficiency will be significantly reduced when the stope length and stope width are reduced to a certain extent. The rock drilling, charging, mining, and other works will also be significantly hindered. So the goal of optimization is to increase the stope structure parameters as much as possible under the premise of safety. Therefore, the displacement failure criterion of 30 mm (Chen et al. 2016) is taken as the "top of the ascent" in SAD. The highest response was 29.31 mm with C o

BBD results
A 3-level, 4-factor BBD model that requires 29 simulations was employed to assess the individual or combined effects between the four significant factors and the response variable D s . The four factors, namely C o , L s , W s , and f r , were considered independent input variables. We can predict the D s using comprehensive factors of the sensitive areas (  Table 7, the simulation order 1 to 24 are factorial points, the order 25 to 29 correspond to the zero factor points, and the zero factor points are used to determine the pure error and the variance. The design matrix of tested variables and the results are represented, where D s is the average displacement of the monitoring points (1#, 2#, 3#, 4#) obtained by using Flac3D, and the predictive values of average displacement D p are calculated by the BBD analysis. The adequacy and the significance level for the regression model were checked using ANOVA. The regression relationship between D s and the independent variables A(C o ), B(L s ), C(W s ), D(f r ) is significant if p value ≤ 0.05 of the model item, and when p value ≤ 0.01, it means that the regression relationship is highly significant. A more minor "Lack of fit" represents a more significant p value, which means the model is more significant with linear and quadratic model terms. The closer the coefficient of determination (R 2 ) is to 1, the better the correlation of the model. When R 2 is greater than 80%, the correlation of the regression equation is thought to be good. Table 8 shows the results, the "F value" of the model was 316.28 and the "p value" < 0.0001, indicating that the model was highly significant, the value of 1.46 showed that the "Lack of fit" was not the significant compared to the pure error, so the model was suitable for making a prediction. Linear terms of A(C o ), B(L s ), C(W s ), D(f r ), quadratic terms of A 2 , B 2 , C 2 , D 2 and interactive terms of AD, BC were significant for D S , and the other terms are insignificant. A larger F value indicates that the factor is more critical in the model. From this, the sensitivity of the four factors is W s >L s >f r > C o .
Based on the data acquired from the Flac3D, a statistical analysis via RSM was used to construct the best-fit regression model for D s as follows: For the programming regression equation: At the same time, the actual regression equation obtained from the programming regression equation according to the actual testing data is as follows: The coefficient of determination R 2 and the adj-R 2 were calculated as 0.9968, 0.9937, respectively, indicating that the testing and the predicted values are in perfect accordance and that RSM may be represented by a quadratic function. Table 7 lists the predicted value D p of the chamber roof displacement obtained by using Eq. (7), and the prediction error is also counted in Table 7 to validate the prediction model, the prediction error ε can be calculated by Eq. (8) as follows: The model has favorable reliability because the statistics in Table 7 shows that the ε of the 16th group is 2.12%, and the rest errors are less than 2%, the prediction error is tiny, and the model prediction accuracy is satisfactory.

Comprehensive study of BBD
It can be seen from the above analysis that D s is not only related to the single factor of A(C o ), B(L s ), C(W s ), D(f r ), but also to the interaction between the four factors. The order of influence of each factor on stope displacement D s can be determined: W s >L s >f r > C o > interaction between C o and f r > interaction between L s and W s .
The three-dimensional response surface (Fig. 8) can be plotted using the regression equation by selecting any two parameters from the four factors above, and the correlation of each variable on D s can be reflected by the horizontal plane of the three-dimensional surface.
The effect of C o and L s on D s is given in Fig. 8a. D s increases when C o decreases and L s increases; the growing trend indicates that L s has a more significant effect on D s than C o . As shown in Fig. 8b, D s increases when C o decreases and W s increases, the curvature of the graph indicates that W s has a more significant effect on D s than C o . The adverse effects of f r and C o on D s are depicted in Fig. 8c, which shows that D s increases when f r and C o decrease and the effect of f r is almost equivalent C o , there is a zero factor point (when C o (0) = 8 MPa, L S (0) = 43 m, W S (0) = 19 m, f r (0) = 45°), which has a mutation phenomenon, that is, the D s increases significantly after the zero factor point. As indicated in Fig. 8d, W s and L s positively affect D s , and the influence of W s on D s is more significant from the response surface increasing trend. The adverse effects of f r and positive effects of L s on D s are depicted in Fig. 8e, the effect of L s and f r on D s is similar to the growing trend of the response surface. Figure 8f shows the influence of W s and f r on D s , displacement 1 3

Engineering applications
In the practical mine production process, obtaining optimal points from the RSM should not be the only goal. Although the most intuitive judgment for promoting the security of the chamber roof is the smallest vertical displacement, we need to increase the structural parameters as much as possible based on security to improve mining efficiency and optimize the benefits of production. The roof deformation of the 0#N stope was predicted using Eq. (7) obtained from the study. For the 0# N stope, the orebody dip angle is 55°, the main parameters of ore properties are: T is 1.2 MPa, C o is 5.3 MPa, and f r is 35°, W p is 4.5 m, H p is 3.6 m, L s is 40 m, and W s is 18 m, as shown in Fig. 10a.The stope deformation under this structural parameter was calculated using the prediction model and verified by numerical simulation. The stope contours of the vertical displacement and the plastic zone distribution are given in Fig. 9.
The roof displacement calculated by the prediction model is 30.18 mm, and the displacement obtained by the numerical simulation is 30.39 mm, they only differ by 0.21 mm, the error between them calculated by formula (8) is only 0.7%, indicating that the prediction model satisfies the high accuracy.
To further validate the accuracy of the model, a monitoring point on the roof center was set on the right side of the pillar in 0#N stope (as shown in Fig. 10b), the measured displacement of the point is 28.89 mm, and the error is 4.3% compared to the prediction value. Field support may explain the more meaningful error in 0#N stope, and arranging multiple monitoring points to take the average value is also conducive to reducing the error.
It can be observed in Fig. 10b that the stope chamber is stable, and the safety of the chamber is guaranteed after the bolt and mesh support. Therefore, the research conclusions in this paper are consistent with the actual situation at the chamber site, and the prediction model has definite guiding significance for the stope safety. It can be seen from the contour of the plastic zone distribution that the pillar is in a shear failure state, and the pillar on the site has a small amount of ore spalling, which is also consistent with the simulation results.

Conclusions
Combined with Flac3D, the experimental design methodologies, including SFT, PBD, SAD, and RSM, were used to evaluate the individual and interactive effects of eight factors related to chamber roof deformation of an underground mine in China. The eight factors are related to orebody properties, orebody dip angle, and stope structure parameters.
(1) The sensitivity interval of multiple factors was obtained through a mix of SFT and numerical simulation. Reliable curve fitting and correlation coefficients were calculated by fitting multiple functions, and the fitting coefficients are all above 0.9. (2) Through PBD, the difference significance parameters of each factor were obtained, and four factors that had a significant impact on the results were screened out: C o , L s , W s , and f r . The center point of the response surface was determined by SAD: C o (8 MPa), L s (43 m), W s (19 m), f r (45°). (3) The numerical simulation was organized using a 3-level, 4-factor BBD model that required 29 simulations, linear terms of A(C o ), B(L s ), C(W s ), D(f r ), quadratic terms of A 2 , B 2 , C 2 , D 2 and interactive terms of AD, BC were significant according to ANOVA and the sensitivity analysis. The polynomial equation accurately described the mathematical relationship between the vertical displacement (D s ) and these significant factors. The order of influence of each factor on D s can be determined: W s >L s >f r > C o > interaction between C o and f r > interaction between L s and W s . (4) The prediction model was applied to predict the roof deformation of 0#N stope in a lead-zinc mine. The deformation relative errors between the prediction model and the numerical simulation, the field measured data are 0.7%, 4.3%, respectively, which indicates that the prediction model is effective and has some reference value for mining safety.