Abstract
Open pit mines have a direct impact on the production efficiency and sustainable development of the enterprise by realm optimization. In the case of changing mining technology conditions and market demand, the original design realm of the mine needs to be optimized in time for the study. Therefore, based on the slope stability analysis, this paper changes the slope angle by changing the working platform size and step and section, and forms a total of four Programs I, II, III, and IV. Geostudio software was utilized to analyze the slope stability of four scenarios and obtain the slope safety coefficients under different scenarios, and the results show that Programs I, II, and III all meet the requirements of reasonable safety coefficients. Furthermore, the geological model of the slope in C profile in Program III was constructed, and the changes of slope stress and displacement during the step-by-step excavation process were analyzed, with a maximum of stress 3.5 MPa, total displacement 12.9 cm, shear should be incremental, plastic zone change is small, and the results showed that the slope structure parameters of Program III were reasonable. The geological drilling data of Jiguanshan molybdenum mine was processed with 3DMine software, the geological database was constructed, the ore body model was established, and the mine boundary optimization study was conducted, Program III was finally determined as the optimal option. The result of the boundary optimization improves the final slope angle by 1°, which can enable the mine to extract 13.334 million tons more ore, reduce the stripping ratio to 1.07m3/m3, and generate economic benefits of $116 million.
1. Introduction
Open pit mining is an economic, simple, and productive mining method, which can be effectively exploited at shallow depths, and there are a large number of open pit mines in China [1, 2]. With the development of mineral resources, most of the open pit mines are facing the problem of increasing the height of slope and mining to the design boundary year by year. In order to realize the safe and efficient exploitation of resources, it is necessary to carry out the slope stability analysis and boundary optimization in time [3, 4]. In open pit mining, the key factor affecting the stability of slopes and the efficiency of mining enterprises is the angle of the open slope size [5, 6]. A small slope angle will lead to an increase in the rock stripping ratio, which will increase the production cost and reduce the efficiency of the mine, and a large slope angle will increase the risk of slope instability and landslide. According to statistics [7, 8], in the process of large open pit mining, every 1° increase in final slope angle, hundreds of millions of dollars of rock stripping cost will be able to be saved. Therefore, on the basis of carrying out slope stability analysis research, reasonable selection of final slope angle of open pit mine is of great significance to improve the production efficiency of mining enterprises.
The stability of open slopes has always been a hot topic in geotechnical engineering research, and its research methods are generally divided into five categories: engineering geological analysis [9], rigid body limit equilibrium method [10], plastic limit analysis [11], numerical analysis [12], and reliability analysis, and the rigid body limit equilibrium method and numerical analysis have been applied by many scholars in recent years. Wu et al. [13] studied the overall stability of the slope in the process of increasing the slope width of the open pit coal mine slope based on the numerical analysis method of the finite difference software FLAC3D and gave an optimized final slope angle of 38° by combining the slope forces and displacements. Chen et al. [14] studied the effects of blast vibration, cumulative rock damage [15], blast vibration wave velocity variation [16], and other factors on slope stability. Shi et al. [17] analyzed the stability of saturated multilayered saturated soil slopes with rainfall infiltration based on numerical analysis method and limit equilibrium method. Based on the model study of instability of soil slope body, Ding et al. [18] explored the effect of seismic action on the stability of slopes. In summary, it can be seen that the study of slope stability not only requires the selection of reasonable calculation methods and analysis software but also requires the comprehensive consideration of a variety of factors, such as the influence of geological formations, hydrological conditions, blasting, and earthquake, in order to ensure the reliability of slope design.
Open pit mine boundary optimization is also an indispensable part of the mining process, which affects the sustainable development of the mine [19]. The calculation and analysis of the boundary optimization is a complex system engineering, which needs to consider various factors such as mine resources, technical aspects, and market aspects [20]. At present, there are floating cone method [21] and L-G graph theory method [22] for surface mine mining boundary optimization, among which L-G graph theory method has been widely used in surface mine boundary optimization, because it has a rigorous mathematical logic and can analyze the boundary model that can obtain the maximum economic value by creating a value model, adding economic constraints and side slope angle constraints.
The stability study of open-air slope and the realm optimization study complement each other and give feedback to each other, and the safe slope angle determined by the slope stability analysis can be used as the basis of realm optimization. At present, most scholars focus on the slope stability problem but rarely combine the slope stability analysis with the realm optimization, which makes the realm optimization lack a certain degree of rationality. Therefore, in this paper, based on the Jiguanshan molybdenum mine, four final slope angle solutions are formed by changing the working platform size and step and section, and the slope stability of the four solutions was analyzed by combining the limit equilibrium method and numerical analysis method. Further processing the drill hole data, constructing a geological database using 3DMine software, generating a three-dimensional model of the surface, the solid and block models of the ore body, and combining the L-G diagram theory by comparing the four design solutions, the optimal design solution was proposed and the final open pit boundary was circled, and the research results can provide reference for similar mines.
2. Project Overview
Jiguanshan molybdenum mine is a large molybdenum mine in Songshan District, Chifeng City, Inner Mongolia, China, with a simple topography and a single type of landform. The design mining depth is +814 m to +210 m, and the elevation of +814 m to +542 m is open pit mining, with a production scale of 1.8 million t/a. There are four groups of fracture structures in the mine area: NEE direction, NE direction, NW direction, and NNW direction, which affect the stability of the mine slope. There is some weathering in the upper part of the mining area, but other parts of the rock are intact. The mine rock is composed of quartz coarse-faced rock, coarse-faced tuff, rhyolite, etc. The rock is hard and dense, with good rock quality and good stability. The water-rich permeability of the aquifer in the mining area is very weak, and there is no surface water body.
The open pit of Jiguanshan molybdenum mine can be divided into 2 mining areas according to the geographical location, the mining area step height is 14 m, and the average side slope angle is 49°. The width of cleaning platform is 10 m, width of transportation platform is 15 m, and slope of transportation road is 8°. +542 m to +640 m elevation step slope angle is 70°, above 640 m elevation step slope angle is 65°, and the slope angle near the step surface is 50°. +542 m to +640 m elevation safety platform width is 3.0 m and above 640 m elevation safety platform width is 3.8 m. The lower plate mining area has been mined to the final boundary of the original design, and the upper plate mining area has been excavated to +654 m, the upper plate of the side slope began to be expanded, and the thickness of the ore body in the upper plate mining area is narrowed from the previous 190 m to 110 m. In general, the eastern part of the ore body has become thinner and the western part of the ore body has been re-detected, and the morphology of the ore body has changed from the previous one. Therefore, it is necessary to carry out timely research on slope angle optimization and realm circling.
3. Rock Mechanics Test and Rock Parameters Determination
As shown in Figure 1(a), representative rocks in the project area were sampled on site, and quartz rough surface rock, rough surface tuff, and rhyolite were processed into standard specimens to carry out uniaxial tensile and compressive tests, triaxial tests, and shear tests. Respectively, with uniaxial equipment using GAW-2000 Rock mechanics testing machine, the acoustic emission signal characteristics of rock loading process were monitored using PCI-2 acoustic emission (AE) system, and the equipment is shown in Figure 1(b): (① GAW-2000 Rock mechanics testing machine, ② LVDT system acoustic emission probe, ③ AE21C acoustic emission monitoring system).
The stress and strain curves of the whole process of compressional rupture of Trachyte tuff were obtained by uniaxial tests, and the relationship between the acoustic emission signal ringing counts and axial pressure is shown in Figures 2(a) and 2(b), respectively. According to the data in Figure 2, the acoustic emission activity of Trachyte tuff is less in the compression-density stage and elastic stage. Before the arrival of the peak strength, when the crack accelerates and the crack gradually penetrates, the acoustic emission activity is the most intense, and the ring count also reaches the peak. After the peak strength, the acoustic emission activity is no longer obvious, because of the serious damage to the internal structure of the rock and the existence of a penetrating damage surface. But after the peak intensity, the acoustic emission activity is no longer obvious due to the serious damage of the internal structure of the rock and the existence of the damage surface.
(a) Stress-strain curve
(b) Stress-acoustic emission curve
The rock mechanical parameters obtained from the test cannot be directly used for the evaluation of slope stability. By combining the joints and fractures in the rock specimen itself, engineering geological conditions, rock type, and other factors, the rock parameters were obtained using discounting methods such as RMR values for rock classification and GSI [23, 24] values for rock structure surface. Rock parameters are shown in Table 1.
4. Stability Analysis of Slopes in Open Pit Mines
Combined with the finite element software GeoStudio, Bishop method, Janbu method, and Morgenstern-price method are selected to calculate the slope safety coefficient. The results show that the Janbu method calculates to the smallest safety coefficient, and the reference Janbu safety coefficient can better ensure the stability of the slope; therefore, only the results of Janbu method are discussed in this paper. Because GeoStudio can only calculate the slope safety coefficient, in order to grasp the parameters of slope mechanics, it is further combined with FLAC3D to construct a three-dimensional geological model of the slope in the key area and carry out the analysis of stress and displacement of the whole process of slope mining.
4.1. Janbu Method Calculation Principle
The Janbu method [25] is a safety factor calculation method that considers all equilibrium conditions and is able to produce arbitrarily shaped sliding surfaces. Figure 3 shows the schematic diagram of the basic principle of the method.
(a) Schematic diagram of strip division
(b) Stress diagram of the th strip
The method first assumes the slide thrust line position, as shown in Figure 3(a), and then strip-separates the slide, where the force state of the th strip is shown in Figure 3(b), which establishes the moment balance equation and derives the slope safety factor by multiple iterations. Its simplified Janbu method is calculated by the formula. where is safety factor, is self-weight of the strip, is pore water pressure, is bottom inclination angle, is width of the strip, is cohesion, is friction angle, is inter-strip fast thrust, is shear force on the slip surface, and is combined force between the strips, .
4.2. Analysis of Factors Influencing Slope Stability
4.2.1. Consideration of Seismic Effects
The rock vibration of the open slope directly affects the stability of the slope and seriously threatens the mine safety production. There are two main causes of slope vibration: blasting vibration caused by frequent blasting in mines and natural earthquakes in the nearby mining areas. The area where the mine is located is seismically inactive and the peak acceleration is 0.15 g. The impact of blasting load is calculated by the following formula. where is blast vibration influence coefficient, is amplification factor spectrum, is maximum peak acceleration, and is gravitational acceleration.
According to Mao et al. [26] and others on the seismic response spectrum site method, coefficient =0.2 substituted into formula (2) can be obtained =0.0336.
4.2.2. Consideration of Groundwater Pressure
According to the hydrogeological survey results of Jiguanshan molybdenum mine, the water-rich permeability of the aquifer in the mining area is very weak, and there is no surface water body, and the slope is less affected by water. Therefore, this calculation is based on a semi-empirical and semi-theoretical analysis idea, favoring safety considerations, and the location of the stable infiltration surface inside the slope is determined using the secondary parabolic depth reduction curve [27] to calculate the hydrostatic pressure inside the rock body.
4.2.3. Determination of Reasonable Side Slope Safety Factor
The value of reasonable slope safety factor is an empirical value obtained based on the summary of a large number of engineering experiences. It is generally believed that the value of slope safety factor obtained through calculation is greater than the reasonable slope safety factor, which indicates that the slope stability is better, and vice versa, the slope has the risk of landslide instability. The determination of reasonable slope safety factor should combine the influence of slope rock parameters characteristics, slope height, groundwater, earthquake, blasting load, and other factors. Therefore, this paper calculates the safety factor of slope by the following three working conditions, respectively. The reasonable safety factor of working condition combination 1 (self-weight + groundwater) is 1.25~1.20, the reasonable safety factor of working condition combination 2 (self-weight + groundwater + earthquake) is 1.23~1.18, and the reasonable safety factor of working condition combination 3 (self-weight + groundwater + blasting power) is 1.20~1.15.
4.3. Slope Angle Design Program
As shown in Figure 4(a), starting from the west end of the open pit mining boundary, eight representative slope profiles are selected clockwise, among which profiles A, B, C, and D are selected for the upper plate and profiles E, F, G, and H are selected for the lower plate. The slope stability model is established according to the section position and stratum distribution. Limited to space, only the geological model map of Program III C is displayed, as shown in Figure 4(b).
In this calculation, the slope angle of the steps from +542 m to +626 m elevation is 70°, the slope angle of the steps above +626 m elevation is 65°, and the local part near the surface is 50°; in order to study the reasonable final slope angle, the platform width is calculated and analyzed according to the implementation of four excavation programs. Program I: safety platform is 12 m wide, cleaning platform is 13 m wide, and transportation platform is 15 m wide. Program II: safety platform is 11 m wide, cleaning platform is 12 m wide, and transportation platform is 15 m wide. Program III: safety platform is 10 m wide, cleaning platform is 10 m wide, and transportation platform is 15 m wide. Program IV: safety platform is 9 m wide, cleaning platform is 9 m wide, and transportation platform is 15 m wide. The final slope angle is shown in Table 2.
4.4. Analysis of Slope Safety Factor Calculation Results
From Figure 5(a), it can be seen that the minimum safety factor for each profile slope in condition 1 (self-weight + groundwater) is 1.38, which is greater than the reasonable safety factor of 1.2 for this case. The minimum safety factor of slope of each profile of Program I is 1.25 in the case of working condition 2 (self-weight + groundwater + earthquake), which is greater than the reasonable safety factor of 1.18 in this working condition. The minimum safety factor in the case of working condition 3 (self-weight + groundwater + blasting vibration) is 1.2, which is greater than the reasonable safety factor of 1.15 for this working condition. The overall safety coefficient values are all greater than the reasonable safety coefficients under the corresponding working conditions, indicating that the stability of the slope is better when this scheme is adopted. From Figure 5(b), it can be seen that with the increase of the slope angle, the slope safety coefficient of each profile tends to decrease, and the overall safety coefficient value is greater than the reasonable safety coefficient under the corresponding working condition, and the slope stability with Program II is better. From Figure 5(c), it can be seen that the safety factor of H-profile of Program III is 1.14 in working condition 3 (self-weight + groundwater + blasting vibration), which is less than the reasonable safety factor of 1.15, and the slope has the risk of instability. From Figure 5(d), it can be seen that the safety coefficient of each profile slope of Program IV is less than the reasonable safety coefficient at working condition 3. The overall safety coefficient of the upper and lower slopes is higher than that of the upper slope, and the safety coefficient is the lowest in the case of working condition (3) (self-weight + groundwater + blasting vibration), so the mine should reasonably arrange the shell holes during blasting to reduce the impact of blasting on the stability of the slope.
(a) Program I
(b) Program II
(c) Program III
(d) Program IV
The safety coefficient of the slope of each profile of Program IV is too small, and the slope is easy to destabilize and damage, and the realm optimization is not suitable. In Program III, profile H can be reduced from 50.87° to 50.37° by widening the first sweeping platform of the slope from 10 m to 12 m, and the safety factor is calculated to be 1.159, which is greater than the reasonable safety factor and can be optimized by the realm. The overall safety coefficient of each profile of Programs I and II is greater than the reasonable safety coefficient, indicating that the slope stability is better under the adoption of this scheme and the realm optimization can be adopted.
4.5. Numerical Analysis of Side Slope Excavation
4.5.1. Modeling
Since Geostudio software can only calculate the safety factor and cannot observe the slope displacement and stress distribution, FLAC3D software was used to carry out further analysis. Combined with the above analysis, the C profile of Program III was selected to establish a numerical model with a model size of 600 m ×50 m ×350 m and excavated to the final realm slope in three steps. Stage 1: from the top to +654 m elevation. Stage 2: from +654 m to +598 m elevation. Stage 3: rock body below +598 m elevation, as shown in Figure 6.
Slope step excavation model lithology and parameters are shown in Table 1, and combined with the rock type meter parameters of Jiguanshan molybdenum mine, this paper adopts the M-C principal structure; the expressions are as follows. where tensile stress is “positive,” compressive stress is “negative,” is maximum principal stress, is minimum principal stress, is cohesion, and is angle of internal friction. When >0, the material is damaged by shear.
The slope of the molybdenum mine in Jiguanshan is at a higher elevation than the surrounding ground surface, so the mechanical boundary is only set as the initial ground stress field under the self-weight of the rock body. The upper boundary is a free surface, and the lower boundary is set as a fixed constraint: considering the area of the actual landslide, the horizontal motion of the right, back, left, and front boundaries of the model is constrained, while allowing its motion in the vertical direction.
4.5.2. Analysis of Simulation Results of Step-by-Step Excavation
As can be seen from Figure 7, after the excavation of stage one is completed, the rock on the excavation surface is removed, the stress on the slope surface is released, and the stress relaxation zone is formed on the slope surface. The maximum principal stress at the bottom of the slope is 1.8 MPa, and no tensile stress appears. The location of the foot of the slope has the largest increment of shear strain, with an increment of about 5.7e10-4, indicating that the most easily damaged area is the foot of the slope. On the slope surface, due to the excavation interface of the slope body, the maximum displacement is about 12.8 cm, and the lower part of the quarry pit blocks the slope’s slide, which ensures the stability of the slope. The plastic zone of the slope is concentrated at the bottom of the slope, and shear damage is the main cause, but no shear damage occurs in the current state of the slope. The overall stability of the slope is good in the overall analysis.
(a) Maximum principal stress cloud
(b) Shear strain increment cloud
(c) Total displacement cloud map
(d) Plastic damage zone
From Figure 8, it can be seen that after the completion of excavation in stage two, the maximum principal stress at the bottom of the slope is 2.7 MPa, which increases by 0.9 MPa compared with the previous stage, and the overall stress is small. The incremental shear strain of the slope is about 7.3e10-4, and the location of the maximum incremental shear strain is also at the foot of the slope. The maximum displacement of the slope surface is about 12.8 cm. The plastic zone is concentrated at the bottom of the slope, and shear damage is the main factor. The analysis of the overall stability of the slope is good.
(a) Maximum principal stress cloud
(b) Shear strain increment cloud
(c) Total displacement cloud map
(d) Plastic damage zone
From Figure 9, it can be seen that after the completion of excavation in stage three, the maximum principal stress at the bottom of the slope is 3.5 MPa, which is 1.7 MPa higher compared with stage one. The increment of shear strain in the slope is about 9.6e10-4, and the location of the maximum increment of shear strain is at the location of the foot of the slope. The maximum displacement of the slope face is about 12.9 cm. The plastic zone is concentrated at the bottom of the slope and is dominated by shear damage, but no shear damage occurs in the final realm slope. Shear damage does not occur in the final realm slope. As the excavation proceeds, the overall slope angle steepens and the step height increases; the slope tensile stress increases and the displacement increases accordingly, indicating that when the general slope reaches instability, the damage of the rock body is mostly gradual, and the plastic shear damage area of the macroscopic homogeneous rock body usually starts from a point at the foot of the slope to expand gradually; the internal force of the rock body is adjusted continuously, and the plastic shear zone extends upward until the overall instability along the slip crack surface. From the numerical analysis results, it can be seen that under the given lithological parameters of Jiguanshan molybdenum mine, the numerical calculation results of excavation show that the slope stress is small, the extension of shear zone is limited, the overall stability of the slope is good, there is no danger of landslide, and the slope Program III is reasonable.
(a) Maximum principal stress cloud
(b) Shear strain increment cloud
(c) Total displacement cloud map
(d) Plastic damage zone
5. Realm Optimization Study
5.1. Realm Optimization Process
The process of open pit realm optimization is using 3DMine mining design software, constructing geological database through mine geological exploration drilling data, circling the solid model, block model, and open pit surface, combining geostatistical methods to assign grade to the block model, and adding economic parameter constraints, slope angle constraints, and software analysis of economic maximization realm. The flow chart is shown in Figure 10.
5.2. Geological Database Construction
The geological resource database was established by collecting and organizing the information from the borehole data of the Jiguanshan mine, combined with the numerical software 3Dmine. In this paper, the existing geological database of the mine is used. By supplementing the exploration drill holes, the geological database is updated, while the surface model uses the mining status quo surface provided by the mine to provide a constraint upper limit for the final realm optimization, as shown in Figure 11.
(a) Drill hole effect
(b) 3D model of the ground surface
5.3. Data Preprocessing
Once the geological database has been created, further preprocessing of the data is required, including generation of combined samples and ultra-high-grade processing.
The combined sample length of Jiguanshan molybdenum mine ranges from 8 m to 12m, with the maximum not exceeding 12 m, and the average of 10 m is taken as the length of the combined sample. The results of the processing of exceptionally high grades have a great influence on the calculation of the statistical parameters of the data, the properties of the experimental variation function, and the grade of the surrounding valuation section blocks, and are an important factor causing errors in the amount of metal and ore. Therefore, the upper grade limit of molybdenum was taken as the exceptionally high value in this paper. Further statistical analysis of the samples was performed and the results are shown in Table 3.
5.4. Variance Function Fitting
Variance function fitting [4] is a mathematical method has been used to characterize the spatial correlation of data values, which is characterized by the inverse of the correlation between data points and their distance in space. The anisotropy parameters of the variance function for the Cockscomb Mountain molybdenum mine are shown in Table 4.
5.5. Block Model Construction
The block property model contains properties such as grade, specific gravity, and lithology inside the block, and can better describe the distribution of properties within the block by preprocessing the grade information of known sample points in the geological database. The block size division of the block model needs to consider the impact on kriging valuation, and choosing a relatively small unit block size can more accurately circumscribe the ore body boundary. The block model parameters determined to be used in this study are shown in Table 5, and the model attributes are set in Table 6.
By analyzing the molybdenum sample values of the Jiguanshan mine, the results showed that the coefficient of variation of the molybdenum metal grade of the mine is large and the compliance is low when using the distance power inverse method, and the ordinary kriging valuation method is used for a more realistic attribute assignment. Using a series of parameters obtained by fitting the above variation function and running the ordinary kriging module in 3Dmine, the final valuation results were saved in the block property Mo, with a statistical average grade of 0.088% Mo. The valuation results and coloring according to the grade property are shown in Figure 12.
5.6. Economic Model and Definition of the Slope Angle
The economic model is the basis for the optimization of the open pit realm, combining the actual production cost, sales cost, ore recovery rate, and other factors to carry out the calculation and analysis of the optimal realm and find the optimal mining realm. The main economic parameters required for the optimization of the realm in this paper are shown in Table 7, and the slope angles are adopted from Program I to IV used in the slope stability study above.
5.7. Analysis of Realm Optimization Results
Combining the block model, economic parameters, and slope angle, the realm optimization analysis was carried out, and an open pit was formed after optimization, as shown in Figure 13. And the results of sales price, ore volume, stripping ratio, and net present value (NPV) within the realm for the four scenarios were further generated, as shown in Table 8.
According to Table 8, in the slope Programs I~IV, as the slope angle increases, the realm stripping gradually decreases and the NPV increases, but the difference between the NPV of the four programs is small. Combined with the previous slope stability analysis, it is known that the slope angle of each profile of Program IV is too large and the slope has the risk of destabilization and damage. The slope safety coefficient of each profile of Programs I~III is greater than the reasonable safety coefficient and the slope stability is better, and according to the NPV results of realm optimization, Program III is the optimal slope program, and Program III is determined as the final realm of optimization, which obtains a net present value of $116 million.
6. Discussion and Conclusion
In this paper, the safety coefficients of typical slope profiles of each partition and numerical simulation of mining step-by-step excavation were calculated by the limit equilibrium method, and the optimization of the slope stability and the overall edge angle of the Jiguanshan open pit was studied. The ore body was re-circulated and the realm was optimized according to the slope angle Programs I~IV, and the following conclusions were obtained through the above research calculations. (1)Using Janbu method, the stability coefficients of the slopes of each profile A~G of the four programs were calculated and analyzed, and the results showed that Programs I~III meet the requirements of reasonable safety coefficients, and the safety coefficients of the slopes of Program IV do not meet the requirements of reasonable safety coefficients. Further constructing the three-dimensional geological model of profile C of Program III, the maximum principal stresses increased continuously during the whole excavation process of the slopes of profile C, with a maximum of 3.5 MPa, displacement 12.9 cm, shear should be incremental, plastic zone change is small, and slope stability is good(2)Realm optimization takes into account the price of molybdenum, mining cost and beneficiation cost, recovery rate, and other factors. The four options designed for the slope angle open pit model were compared and analyzed, and Program III was finally determined to be the best, which can improve the overall slope angle of the mine by about 1°(3)The mine was able to extract 13.3340 million tons more ore after the boundary optimization, and the stripping ratio was reduced to 1.07 m3/m3, generating an additional economic benefit of $116 million
Data Availability
The experimental data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare no conflict of interest.
Authors’ Contributions
Yongyue Hu was responsible for the experiments and data analysis; Changhong Li was responsible for the methodology and conceptualization; Juzhou Li and Dayu Long were responsible for the visualization, data curation, and resources; Yu Wang was responsible for the supervision, funding acquisition, and project administration.
Acknowledgments
This study was supported by the National Natural Science Foundation of China (52174069), Beijing Natural Science Foundation (8202033), and National Key Technologies Research & Development Program (2018YFC0808402).