key: cord-0074358-ebbcyjn9 authors: Jin, Ai-Bing; Liu, Bing; Gao, Yong-Tao; Sun, Hao title: Evaluation of safety and deformation characteristics of the secondary stope sandwiched between backfills in underground iron mines date: 2022-01-31 journal: Bull Eng Geol Environ DOI: 10.1007/s10064-022-02581-7 sha: b7dc068d18584de66cb215b328ac7d12f18db76e doc_id: 74358 cord_uid: ebbcyjn9 In a large underground iron mine (over 2000, 000 t/a) in China, a large volume of mined void could not be filled in time, which posed a threat to the working face. Physical and numerical models were used to investigate the stability of the mined void sandwiched between the backfills. Digital image correlation (DIC) method was employed in the physical model to calculate the deformation field of the roof and backfills. Stress meters and multi-point borehole extensometers were installed in the mine to observe the change of vertical stress and horizontal displacement in real backfills. Principal component analysis (PCA) and cluster analysis were applied to classify the deformation state around mined void. The deformation field was exhibited by the contour and curves of displacement, which was consistent between the physical model and numerical simulation. The contour indicated that highly deformed zones were ellipse-shaped and wedge-shaped in the backfills and roof respectively. The deformation process was classified into three states (stable state, transition state, and unstable state) by using the density-based clustering. The data obtained from the in situ monitoring showed that the vertical stress increased by 60–80 kPa and the horizontal deformation was small in the backfills close to mined void. Classification of deformation states indicated that the mined void (temporary opening) remained stable during production plan. The stability of mine-out void in underground mines is an important issue in mining activities and surface subsidence (Li et al. 2004; Bakhtavar and Yousefi 2019; Xu et al. 2015) . In early 2020, the pandemic COVID-19 has negatively influenced the importation of iron ore rocks into China. Consequently, this limitation has stimulated the price of ore rocks and caused enlarged production in mines. In some underground mines, the filling up of mined voids could be delayed due to the limited filling capacity and production plan (Brady and Brown 2006) . From the early 1990s, backfilling has been developed at underground mines, including Mount Isa Mine in Australia (Villaescusa 1996) and Neves Corvo Mine in Portugal (Been et al. 2002) . This study is based on an underground iron mine located in China in which the filling of the mined void has been delayed. Figure 1 shows the unmined pillar and mined-out void in the orebody. Through long-term mining activities, a large volume of the mined void was produced. While part of the void was filled with cemented tailings backfills (CTB), the remainder was left unfilled. For the pillar sandwiched between the backfills, the pillar recovery (secondary stopes) generally began after the adjacent void (primary stopes) had been filled for more than 1 month and the desired load-bearing capacity of the CTB was actualized. For one stope with dimensions of 80 m × 60 m × 18 m in the mine, mining activities could last for 2 months. During the pillar recovery, the collapse of the roof and backfills may occur due to the change in the ground pressure and blast vibration. The development of mathematical methods provided a basis for assessment of underground opening and strength of surrounding rock (Jiang et al. 2020; Yang et al. 2019; Zhang et al. 2020; Zhou et al. 2020; Mazraehli and Zare 2020) . But more studies were required for assessing the stability of secondary stopes during service life of underground mines. The physical model test is a traditional method used to analyze the stability of underground openings, such as tunnels and underground mines (Yang et al. 2018; Kang et al. 2018; Xie et al. 2020) . Deformation measurement and stress monitoring are the two main techniques applied in physical model tests. Digital image correlation (DIC) method was developed about 30 years ago (Keating et al. 1975) . As an optical technique, DIC method has been previously employed in tests of small-scale specimens and physical model tests to obtain the change in the deformation field under loading (Zhang and Zhao 2013; Liu et al. 2020b; Xiang et al. 2018; Huang et al. 2020b ). Numerical simulation plays an important role in the mechanical analysis of geotechnical engineering (Cundall and Strack 1979) . For an underground space, Particle Flow Code (PFC) could be used to model the deformation of the opening and fracture of the surrounding rock, such as tunneling and mined void (Xiang et al. 2018; Wang et al. 2017) . In this study, the physical model test was conducted for the mined void sandwiched between the backfills. The deformation field around the void was calculated using DIC and the inner stress was recorded by the aid of stress sensors and a continuous acquisition system. As a comparison, PFC confirmed the deformation field of physical model. Based on the curves of displacement field obtained by DIC, multivariate statistical analysis was applied in the assessment of deformation field. Cluster analysis and principal component analysis (PCA) were used to classify the deformation state. The results of the in situ monitoring in the underground mine were also listed to provide reference data for similar underground mines. This study focused on the stability of the roof and backfills around the mined void in the underground iron mine. Figure 1 shows the mined void in the orebody, which lies about 300-600 m below the ground. For the ore rock, the main metallic mineral was magnetite, while the gangue minerals were mainly lizardite, clinochlore, and calcite. For a complete orebody, the average ore grade was about 54% and the surrounding rock contained skarn and diorite. Due to the high grade, the average density of the ore rock exceeded 3600 kg/m 3 . The strike direction of the orebody and the dip angle was recorded as 78° and 10-30°, respectively. Figure 2a shows the physical model and the location of the stress sensors. The model dimensions were 1 m × 1 m × 0.4 m. The ore rock, surrounding rock, secondary stope, and backfills were denoted as OR, SR, SS, and BF, respectively, in the different figures and tables. The geometric similarity ratio, volume-weight similarity ratio, and stress similarity ratios in this physical model were determined as C l = 200, C γ = 1.48, and C σ = 296, respectively. The height and width of the stopes and mined voids in the mine were mainly 40 and 20 m, respectively. Thus, the height and width of the excavated stopes were determined as 200 and 100 mm in the physical model, respectively. The thickness of the roof in the mine was about 8 m, and similarly, the roof was determined as 40 mm in the physical model. Table 1 highlights the basic mechanical parameters and density values for the ore rock, surrounding rock, and CTB based on uniaxial compressive tests. Table 2 shows the mix proportions and parameters for different materials in the physical model. In Table 2 , the alternative addition, denoted as AA, was barite powder for OR and SR, while that of the BF was sawdust to adjust the density. The stress sensors were installed on the backfills near the void. Four measurement cross-sections were made in two backfills close to the void. While the upper two crosssections, denoted as B and D, were located 10 mm away from the roof, the other two cross-sections were located 60 mm away from the roof. The distance between each sensor in a cross-section is displayed in Fig. 2a . In the following sections, the sensors in each cross-section were denoted as A-1 and so on. As shown in Fig. 2b , per layer mixture were added to the loading frame and the model was compacted. Next, a wooden beam was fixed on the loading frame for subsequent layer mixture (Fig. 2c ). Before excavation, irregular-shaped speckles were prepared on the surface of the model. Figure 2d shows the distribution of the speckles for the DIC measurement. An acrylic plate was fixed in front of the model. During the preparation of the physical model, the secondary stope (unmined ore reserve) was initially sandwiched between two backfills and then initial stress was applied to the whole model. A board was fixed in front of the secondary stope to mark the position of the stope and make the excavation easy. Moreover, speckles should not be prepared on the surface of the secondary stope and this board could avoid creating speckles on the surface of the stope accidentally. After the initial stress reached the required stress, the board was removed and the stope was excavated. The images of speckles and the Figure 3 shows the size distribution of the particles in the PFC model. The region marked by a white rectangle corresponded to the secondary stope (unmined ore reserve). The particles of the roof (ore rock) were surrounded by a yellow rectangle with a thickness of 4 mm. The regions marked by three black rectangles were designated as backfills. One region was located above the roof while the secondary stope was sandwiched between the others. Except for the regions of backfills and ore rock, the remaining part of the PFC model corresponded to surrounding rock. The dimensions of secondary stopes (unmined ore reserves) and backfills (mined primary stopes) were 100 × 200 mm (width × height). The geometry of the PFC model was consistent with the physical model in Fig. 2a . In total, regions close to the secondary stope were filled up with small particles to model the failure process, while the particles in other regions were relatively large in size, to reduce the number of particles. The region of surrounding rock was majorly made up of particles having a radius of 0.9 ~ 1.2 mm, while that of the particles in the region of the roof was less than 0.6 mm. For the region of backfills close to the secondary stope, the radius of particles was also below 0.6 mm. In the other region of BF, the radius of particles was within the range of 0.6-0.9 mm. A PFC model consists of particles and contacts. In PFC simulations, the contact models include linear contact, contact bond, parallel bond, flat joint, and smooth joint. Each contact model can simulate mechanical behavior of different materials. For example, linear contact only provides contact force (including compressive contact force and frictional behavior) and it can be used to reproduce the deformation of soil generally. Contact bond model and parallel bond model (PBM) can be subjected to tension and they have been widely used in rock mechanics and fracture mechanics (Potyondy and Cundall 2004; Itasca 2015) . Compared with contact bonds, parallel bonds can also transmit bending moments between particles. In this study, PBM was applied between the particles. When the shear/tensile stress exceeds the shear/tensile strength, the bond breaks and a micro-crack forms. After a number of micro-cracks initiate and propagate, the whole PFC model fails. This criterion possesses the advantage that the crack initiation and propagation of brittle materials can be reproduced spontaneously under a wide variety of loading conditions. For underground engineering, PBM has been used to analyze long-term stability and failure process of tunnels and mined void under static loading Huang et al. 2020a) . Compared with previous studies (Liu et al. 2017; Xu et al. 2020; Huang et al. 2020a) , the mesoscale parameters of weak rock-like materials can be adjusted to reproduce the mechanical behavior of similar materials. The procedures of calibration process include the genesis of rectangular models (100 mm × 50 mm), the simulated uniaxial compression tests, and comparison of mechanical parameters. The main mesoscale parameters of the PBM are the elastic modulus ( E ), bond normal-to-shear stiffness ratio ( k n ∕k s ), bond tensile strength ( ̃c ), bond cohesion strength ( c ), and friction angle ( ). The elastic modulus was matched by setting ̃c and c to a large value and varying E . The peak strength was matched by varying ̃c , c , and . After a few iterations, when the uniaxial compressive strength (UCS) and elastic modulus of a rectangular model were close to the results in Table 2 , the micro-parameters of the rectangular model were determined as input micro-parameters in simulation of the physical model. Table 3 lists the density (m d ) of the balls and mesoscale parameters of the PBM for OR, SR, and BF in the PFC model. This study mainly focused on the stability of the roof of the backfills lying at a distance of 400 m below the surface level. The in situ stress was estimated to be about 15 MPa; thus, the mined void was excavated after the Y-and X-stress attained a value of 50 kPa. Figure 4a and b show the X-and Y-displacement near the void after the excavation. In the X-displacement field, an ellipse-shaped zone formed in the left backfill, while in that of the Y-displacement field, a wedge-shaped zone formed in the roof. In these zones, the displacement was relatively larger than the other region, which indicated that the adjacent excavation led to the sudden deformation of the roof and backfills. The polygonal zone marked in Fig. 4b could be regarded as the region affected by the excavation process, which would be also observed in the numerical simulation. With the monitoring time increasing, the displacement of the roof and backfills became large. In Fig. 4c and e, the X-displacement increased within the ellipse-shaped zone, partially reaching up to 2 mm. The ellipse-shaped zones in the X-displacement field indicated the large and small amounts of the X-displacement in the middle and on the top of the backfills, respectively. From Fig. 4b -d, the Y-displacement of the roof increased partially and locally at the early stage of the excavation process. From Fig. 4d -f, the Y-displacement of the whole roof increased totally and this phenomenon was caused by the settlement of the upper materials and the self-deformation of the roof. To measure the displacement quantitatively, the average displacement was calculated as shown in Fig. 5 . Small measurement zones were selected in the adjacent backfills and roof. For the measurement zones in the roof, the results of four measurement zones in one column were plotted in one figure and the measurement zones were denoted as U1-U4 from the bottom to the top. For the measurement zones in the backfills, the X-displacement of the five measurement zones in one row was plotted in one figure. For the right backfills, the five measurement zones in one row were denoted as R1-R5 from the left to the right, while for the left backfills, they were denoted as L5-L1 from the left to the right. Based on the magnitude of the displacement, the curves of the average displacement contained 88 data points and they could be divided into three stages as shown in Fig. 6 . The early stage of the excavation was denoted as S-1 (image 1 ~ image 37). When the displacement around the void became unstable, this stage was denoted as S-3 (image 81 ~ image 88). The stage between S-1 and S-3 was denoted as S-2 (image 38 ~ image 81), which could be regarded as the transition stage from stability to instability. The initial increase of the X-displacement of the backfills was observed at the beginning of S-1, with little difference between each other. After S-1, the difference in the X-displacement (selfdeformation) increased. The starting points of S-2 and S-3 were determined by the magnitude of self-deformation. At the starting points of S-2, self-deformation began to increase and the magnitude was small. The X-displacement of L1, L2, and R1-R3 remained steady or increased, while the X-displacement of the other measurement zones decreased slightly. The curves of the average displacement revealed that parts of two backfills near the void showed a tendency to move toward the void at S-2. At the starting points of S-3, self-deformation became large and unstable. During this stage, the X-displacement of L1-L5 increased significantly, which indicated a tendency of the left backfill to be unstable. In Fig. 7 , the Y-displacement of the roof induced by excavation was about 0.2 mm in S-1. This result was close to the deformation observed in other physical models. For example, the displacement of the roof in a physical model (C l = 170) based on Sijiaying Iron Mine was in the range of 0.05-0.4 mm (Li et al. 2015) . After an initial increase in S-1, the curves of the Y-displacement almost became flat. At the beginning of S-2, the difference in the Y-displacement was induced between the curve of U1 and the other curves as shown in Fig. 7c and d. This difference was also observed in the middle of S-2, as shown in Fig. 7a . This phenomenon indicated that the self-deformation of the roof began increasing. At S-3, the difference in the Y-displacement became notable among the four curves in each plot (Fig. 7a-d) , which indicated that the self-deformation of the roof became unstable. Due to large deformation, the DIC calculation became invalid after S-3. Figure 8 shows the change of the vertical stress in the backfills near the void. After excavation, the Y-stress for cross-section A-C increased by about 1-2 kPa in S-1. For cross-section D, the Y-stress of D-2 and D-3 increased by about 7 kPa, while the increase in Y-stress of D-1 and D-4 was smaller than 1 kPa. This discrepancy could be related to the uneven compaction of materials during the completion of the physical model. At S-2, the increase of the Y-stress in cross-sections C and D was more significant than the other two, which was an indication that the inner stress increased more notably for the left backfills. At S-3, Y-stress in the backfills changed a lot. In one backfill, the stress of some sensors increased while the others As shown in Fig. 9 , the numerical simulation also displayed a similar deformation field in comparison with the physical model. Figure 9a and b show that parts of backfills and the roof moved toward the void at the early stage of opening. Figure 9a shows the measurement circles installed in the model. The X-displacement of the backfills was represented with m1-m18, while that of the Y-displacement was denoted with m19-m27. As shown in Fig. 9b , a polygonal zone, which was like the phenomenon depicted in Fig. 4b , also formed around the void. Figure 9c and d show the evolution of the displacement field and the crack initiation at S-3. The X-displacement in the ellipse-shaped zones increased and the settlement also formed in the roof. Figure 9e and f show the crack coalescence and interaction after S-3. As a comparison, fracture process of the physical model was also displayed in Fig. 10 . In Fig. 10a , the breakage could be observed in backfills. From Fig. 10b to c, wedge-shaped breakage occurred. From the view of deformation field and fracture process, the similarities of between the PFC model and the physical model provided confidence that the two models behave realistically. Figure 11 shows the curves of the Y-displacement recorded represented by m19-m27. After excavation, the Y-displacement was about 0.2 mm and a flat region was observed in the curves, which was close to the results obtained through the DIC. When the time step of the calculation was less than 1.8 × 10 5 , the Y-displacement was almost the same for m19-m27. When the time step exceeded 1.8 × 10 5 , the difference in the Y-displacement increased for all the curves, which is an indication of the increment of the self-deformation of the roof. In Fig. 12 , the horizontal displacement for the two backfills significantly increased to about 0.25 mm during the early stage of opening, and then, the curves of horizontal displacement changed slightly. After the calculation of 1.8 × 10 5 time steps, the X-displacement of the measurement circles close to the void, such as m9 and m16, increased sharply, which revealed that the parts of the backfills close to the void showed a tendency to collapse. This phenomenon showed a good agreement between the experimental and numerical results. PCA is a dimensionality reduction method to explain the variability of high-dimensional data by a few linear combinations of original variables, which are defined as the principal components. In most multivariate cases, two or three principal components, fewer than the number of original variables, can almost express the variability of the original sample with little loss of information. The main steps of PCA are briefly described below (Anderson 2003; Härdle and Simar 2007) . is defined as: Denote the eigenvalues of this matrix as 1 , 2 , ⋯ , p and the corresponding eigenvectors as = 1 , 2 , ⋯ , p , where i = u 1i , u 2i , ⋯ , u pi � , i = 1, 2, ⋯ , p , then the principal components = 1 , 2 , ⋯ , p can be written as: Cluster analysis has been widely used in image processing, pattern recognition, and machine learning to classify elements into clusters on the basis of their similarity (Bewley and Upcroft 2013; Rodriguez and Laio 2014; Dong et al. 2018; Sepúlveda et al. 2018) . Typical clustering methods included K-means clustering, grid-based clustering, and distributionbased clustering. However, some classical clustering methods, such as K-means clustering, can be not suitable for nonconvex-shaped data points. As a comparison, density-based clustering (Kriegel et al. 2011) behaves well for data points with complex geometries. The results of DIC analysis shown Figs. 6 and 7 were analyzed by the combination of density-based clustering with PCA. Density-based clustering might be not suitable for highdimensional data because the number of variables affected (2) = � (3) y 1 = u 11 X 1 + ⋯ + u p1 X p the measure of the similarity between data points. PCA could reduce the number of variables for high-dimensional data before density-based clustering. The number of curves in Figs. 6 and 7 was 56 and the number of images was 88. In this assessment, each image was regarded as an observation and curves corresponded to variables. The sample size was 88 and the number of variables was up to 56, correspondingly. Four indexes were used to assess the stability of the left backfill ( I if ), the right backfill ( I rg ), the roof ( I rf ), and the whole mined void ( I mv ). Based on the number of curves in Figs. 6 and 7, I if , I rg , I rf , and I mv included 20 variables, 20 variables, 16 variables, and 56 variables, respectively. The values in input matrix equaled the ratio of original observations to the peak value. By applying PCA, q 1 for y lf , y rg , and y rf were 0.985, 0.978, and 0.998, respectively, which indicated that the first principal component could explain more than 95% of total variability. The first principal component of I if , I rg , and I rf , denoted as y lf , y rg , and y rf , was as follows: (4) y lf = 0.240X 1 + 0.225X 2 + 0.202X 3 + 0.194X 4 + 0.186X 5 +0.297X 6 + 0.247X 7 + 0.242X 8 + 0.232X 9 + 0.217X 10 +0.264X 11 + 0.234X 12 + 0. y rf = 0.271X 1 + 0.237X 2 + 0.234X 3 + 0.234X 4 + 0.258X 5 +0.252X 6 + 0.240X 7 + 0.238X 8 + 0.285X 9 + 0.257X 10 +0.242X 11 + 0.240X 12 + 0.271X 13 + 0.258X 14 + 0.237X 15 +0.238X 16 Fig. 12 (a, b) Curves of the horizontal displacement with respect to time in PFC (a) Horizontal displacement field of the left backfills (b) Horizontal displacement field of the right backfills The density-based clustering was applied to classify the state of deformation field in a 3D coordinate system ( y lf , y rg , and y rf ) based on the observations of the first principal component. During the whole process of deformation, the observations of deformation field changed markedly from the stable state to the unstable state. The criterion for classifying the state was that the observations in one state were close to each other while the dissimilarity was relatively great among observations of different states. A new densitybased clustering method, called method I, was applied in this section (Rodriguez and Laio 2014) and it was based on two assumptions: (1) cluster centers are surrounded by neighbors with lower local density (ρ); (2) The centers are at a relatively large distance (δ) from any points with a higher local density. The assumptions were in agreement with the criterion. Figure 13 shows the identification of cluster centers in method I. Figure 14a shows the sequence of data points in the 3D coordinate system. In Fig. 14b , four clusters were called db-1 ~ db-4 for method I and the classification of deformation process included three states. Stable state corresponded to db-1 (the deformation process before image 21) and unstable state corresponded to db-4 (the deformation process after image 81). The process from image 21 and image 81 could be regarded as transition state. In the Appendix, the comparison between two clustering methods also indicated that the deformation process could be classified into three states. Compared with Figs. 6 and 7, the stable state, transition state, and unstable state approximately corresponded to S-1, S-2, and S-3, respectively. The starting point of transition state was earlier than S-2 and this difference indicated that the classification of stable state was relatively conservative by multivariate statistical analysis. This classification could provide an early warning for increase of self-deformation. Method I provided a method to measure the similarity of deformation field in one state. The local density ρ of db-1, in the top left corner of Fig. 14b , was highest among four clusters as shown in Fig. 13 , which indicated that the variation of y lf , y rg , and y rf was small in the stable state. For the unstable state, the local density ρ of db-4 was lowest and the data points were far apart in the bottom right corner of Fig. 14b . This phenomenon indicated that the deformation field changed markedly during this stage. In situ monitoring Figure 15 shows a backfill-related monitoring program in the underground mine. The monitoring devices included stress meters (GK-4300BX) and multi-point borehole extensometers (BGK-A3), provided by China Geokon Instruments Co., Fig. 13 Identification of cluster centers in method I Fig. 14 (a, b) The state of deformation field classified by the clustering method (a) The index of images (b) Result of cluster analysis Ltd. The monitoring devices were installed on the top of the backfills before mining activities progressed in stopes 1 and 2. Two measurement points were selected in the backfills and they were denoted as CS-a and CS-b. The boreholes for monitoring devices were drilled from stope 2 because drill headings have been completed on the top of stope 2. The length of the connecting rods between the anchors and displacement transducers was 20, 16, 12, and 8 m in one borehole. CS-a and CS-b were denoted as CS-a-1 to CS-a-4 and CS-b-1 to CS-b-4, respectively. The length of boreholes was about 5 m for stress meters. After the installation of devices, the boreholes were filled with cement grout. The mining direction was also marked in Fig. 15 and the ore rock near CS-b would be mined before CS-a. The acquisition lasted for 42 days and the interval was about 7 days. Figure 16 shows the lateral displacement of the backfills in CS-a and CS-b. The displacement in CS-a was relatively small, which indicated that the deformation of adjacent backfills was small. In Fig. 16b , the third acquisition indicated that the displacement in CS-b-1 to CS-b-3 increased by about 4 mm, whereas that in CS-b-4 increased by about 3 mm. The sudden increase could be attributed to the excavation of adjacent ore rock and the fracture induced by a nearby blasting process. The curves of the displacement Fig. 15 Installment of the displacement transducers and the borehole stress meters Fig. 16 (a, b) Curves of the horizontal displacement of the backfills with respect to time (a) Deformation of CS-a (b) Deformation of CS-b remained steady in the following acquisitions, which were consistent with the curves of S-1 in the physical and numerical investigations as shown in Figs. 6 and 12. Figure 17 shows the stress redistribution on the top of the backfills caused by the excavation of stope 2. The incremental vertical stress increased to 60-80 kPa. Considering that the UCS of backfills was about 3 MPa in indoor tests (Liu et al. 2020a; Wang et al. 2020; Xue et al. 2020 ), the small change of the vertical stress might not have a significant influence on the backfills. The increase in the vertical stress in CS-b occurred mainly during the first to the third acquisition. For CS-a, the increase mainly occurred during the fourth to the last acquisition. The significant increase in the vertical stress of CS-b occurred earlier than that of CS-a and the lateral displacement in CS-b was also relatively large. This phenomenon might be attributed to the earlier excavation of the ore rock near CS-b. The stability of the backfills could be assessed considering the displacement and stress. In the physical model, the maximum displacements of the left and the right backfills were about 0.5 and 0.25 mm during stable state, respectively. In the numerical simulation, the maximum displacement was about 0.25 mm after 5 × 10 4 time steps. The geometric similarity ratio C l was determined as 200 and real lateral deformation was in the range of 40-100 mm during stable state. The results of deformation monitoring in the mine showed that the backfills were at the early stage of stable state. In the physical model, the increase in the vertical stress was generally close to 1 kPa in the backfills at the beginning of S-1. The stress similarity ratio C σ = 296 and the real vertical stress was nearly 300 kPa, larger than the stress of in situ monitoring, which also indicated that the backfills were stable after excavation. As shown in Table 4 , the vertical displacement of the roof was also recorded by field measurements, such as multipoint borehole extensometers. These field measurements in underground mines showed that the displacement of the roof was mainly in the range of 3-40 mm. In the physical and numerical model, the deformation of the roof induced by excavation was about 0.2 mm at an early stage. C l was 200, which suggested that the in situ deformation could reach 40 mm. Based on the curves of the vertical displacement in the DIC results and the numerical model, the vertical displacement of the roof consisted of two parts. While a part was related to the compaction of upper materials, the other part was the self-deformation of the roof. The self-deformation of the roof could be estimated by the difference between the curves of vertical displacement, as shown in Figs. 7 and 11. In the physical test, the difference was small at S-1 but got notable at S-2 and S-3. This phenomenon could be also observed in the numerical model. This result indicated that the roof was stable at the early stage of the excavation; however, a large self-deformation of the roof could occur after a long-time opening. (Yu et al. 2018) Diorite 15 m −171 m 3 mm Iron mine (Yu et al. 2014) Marble 15 m −230~−290 m 5-25 mm Gold mine (Xie 2015) Syenite 20 m −1264~−1224 m 12~40 mm Chongli Zijin Mine (Sheng et al. 2016) Syenite 20 m −1344 m 6~20 mm Due to the high demand for iron ore rock in China, a large volume of mined voids are left unfilled. The stability of the roof and backfills near the mined void was investigated through a physical model, numerical simulation, and field measurement. The main conclusions herein are listed as follows: 1. The physical and numerical model indicated that the excavation of the stope sandwiched between the backfills led to a limited deformation of the roof and backfills at the early stage of the excavation process. After a longtime opening, stress redistribution could be observed in the backfills and the displacement of the backfills and the roof increased significantly, which poses a threat to the stability of the mined void. 2. In the displacement field obtained by DIC and numerical simulation, the highly deformed zones, which were ellipse-shaped and wedge-shaped, formed in backfills and roof, respectively. The curves of displacement field indicated that the self-deformation of the model changed significantly after the early stage of excavation. 3. Multivariate statistical analysis was applied to assess the deformation field of backfills and roof in the physical model. The density-based clustering could classify three states of deformation field, including stable state, transition state, and unstable state in the physical model. The characteristics of deformation field could also be measured by local density of data points in the densitybased clustering method. 4. The results of in situ monitoring indicated that the mined void was at the early stage of stable state. After the extraction of ore rock, mined void should be filled immediately to avoid the long-period opening. Analysis of ground vibration risk on mine infrastructures: integrating fuzzy slack-based measure model and failure effects analysis Liquefaction potential of paste fill at Neves Corvo mine Portugal Advantages of exploiting projection structure for segmenting dense 3D point clouds Rock mechanics for underground mining A discrete numerical model for granular assemblies Clustering based on grid and local density with priority-based expansion for multi-density data Applied multivariate statistical analysis Instability mechanism of shallow tunnel in soft rock subjected to surcharge loads Experimental analysis of progressive failure behavior of rock tunnel with a fault zone using noncontact DIC technique PFC-Particle Flow Code in 2 and 3 dimensions Version 5.0 Documentation Set of version 5.00.21 Evaluation of safety and deformation characteristics of cemented tailings backfill mining disturbed area near shafts: a case study in China A physical and numerical investigation of sudden massive roof collapse during longwall coal retreat mining An improved method of digital image correlation Density-based clustering Engineering geology, ground surface movement and fissures induced by underground mining in the Jinchuan Nickel Mine Similar material model research on the stope stability at cut and fill phase Influence of water loss on mechanical properties of superfine tailing-blast-furnace slag backfill Numerical simulation of stressstrain behaviour of cemented paste backfill in triaxial compression Experimental and numerical study on crack development characteristics between two cavities in rock-like material under uniaxial compression An application of uncertainty analysis to rock mass properties characterization at porphyry copper mines A bonded-particle model for rock Clustering by fast search and find of density peaks Fuzzy clustering with spatial correction and its application to geometallurgical domaining Practice on the stability monitoring technique of laminated goaf during pillar recovery Excavation design for bench stoping at Mount Isa mine A tutorial on spectral clustering Height of the mining-induced fractured zone above a coal face Mechanical behavior, acoustic emission properties and damage evolution of cemented paste backfill considering structural feature Application of transparent soil model test and DEM simulation in study of tunnel failure mechanism Study on earth pressure of stope and sedimentation transformation law of a certain underground mine Roof deformation associated with mining of two panels in steeply dipping coal seam using subsurface subsidence prediction model and physical simulation experiment Experimental and numerical studies on the reinforcing mechanisms of geosyntheticreinforced granular soil under a plane strain condition A computational model of safe thickness of roof under filling body based on cusp catastrophe theory Mining induced strata movement and roof behavior in underground coal mine Fiber length effect on strength properties of polypropylene fiber reinforced cemented tailings backfill specimens with different sizes Physical experiment and numerical modelling of tunnel excavation in slanted upper-soft and lower-hard strata Large-scale model experiment and numerical simulation on convergence deformation of tunnel excavating in composite strata Underground pressure law and stability monitoring in caving to backfilling method overlaying open-pit Space-time rule of the control action of filling body for the movement of surrounding rock in method of the delayed filling open stopping. Caikuang Yu Anquan Gongcheng Xuebao Fracture pattern of overlying strata in multiple coal seam mining in a physical model vis-àvis MATLAB analysis and geological radar Determination of mechanical properties and full-field strain measurements of rock material under dynamic loads Prediction of rockburst risk in underground projects developing a neuro-bee intelligent system As a comparison, spectral clustering (Von Luxburg 2007), called method II, was also used to classify the state of deformation field. Method II included dimensionality reduction before clustering, which made it suitable for high-dimensional data. The input matrix in spectral clustering was I mv (88 × 56). The results of method II were displayed in the same 3D coordinate system.The result of method I is shown in Fig. 18a . In Fig. 18b , four clusters were called sp-1 ~ sp-4 for method II. The comparison between two methods indicated that the classification of stable state and unstable state was consistent with method I. The main difference between two methods was the separation of transition state and this discrepancy might be caused by two factors. One factor was the different theory of two clustering methods. The other one was the different input data in two methods. In method I, the input data were y lf , y rg , and y rf after PCA. In method II, the input data were the original data (88 × 56). Although q 1 of the first principal component was larger than 95%, a little unexplained variation might still affect the results of clustering. Availability of data and material Not applicable. Fig. 18 (a, b) The comparison between two clustering methods (a) Method I (b) Method II Code availability Not applicable. The authors declare no competing interests.