Abstract

Roadway deformation and rock burst are the two key challenges faced by the safe operation of coal mines. Aiming at the issue of large deformation of the gob-side roadway under coal pillars in multiseam mining, this study has considered the case of the 8308 panel of Xinzhouyao coal mine in China. Based upon a combination of theoretical analysis, numerical simulations, and engineering practices, the mechanical model of “stress and deformation quantitative calculation of gob-side roadway under overlying coal pillars” was established in this study. The analytical solutions of the vertical stress distribution and the plastic zone of the gob-side roadway under overlying coal pillars were obtained. Finally, the accuracy of the mechanical model was verified using numerical simulations. The results showed that the coal pillar, upright above the gob-side roadway, and the cantilever roof around the gob-side roadway were the main factors leading to stress concentration and deformation around the gob-side roadway. For the particular cases considered in this study, the peak stress of the gob-side roadway could reach 1.8 times of the self-weight stress of overlying strata. The rates of the contribution of the gob-side roadway’s overlying pillar and the cantilever roof around the gob-side roadway to peak stress were 78.3% and 16%, respectively. The obtained results have an essential reference significance for stress calculations and rock burst prevention design of gob-side roadway under overlying coal pillars in multiseam mining.

1. Introduction

In the next few years, coal will still be the basic energy in China [13]. With the decrease in coal reserves in the first mining coal seam, the number of coal mines mining under the upper mined coal seam is rapidly increasing [46]. During mining in multiseam, the coal pillars’ leftover from multiseam mining affects the stress distribution and contributes to a more complicated condition, which makes it easy to cause large local deformation of the roadway or even rock burst in certain regions [711]. The existing theoretical models are unable to systematically calculate the stress distribution under the effect of overlying coal pillars in multiseam mining. Therefore, it has become an urgent task to establish a quantitative calculation method for determining the deformation and failure of roadways under the effect of the overlying coal pillars.

Large deformation of roadway is a feature of deep mining, and it will also be encountered in multiseam mining. There are many research results about deep mining, which can provide a reference for this study. Some scholars have studied the deformation mechanism of deep soft rock roadway using the angles of theory [1214], numerical simulation [1517], and field measurements [18, 19]. A previous study has also studied the measures to control the deformation of roadways [20, 21]. These theories and practices have sped up the research of roadway failure mechanisms. Generally, the first mining coal seam is considered relatively stable in terms of both vertical and horizontal stresses and ubiquitous deformation in the whole roadway. However, in the close multiseam mining, the residual coal pillar and goaf in the upper seam will cause the redistribution of stress in the lower seam, which will lead to higher stress concentration under the overlying coal pillars and lead to large deformation or rock burst in the roadway [6]. As a result, it is necessary to establish a quantitative calculation method for the stress distribution of multiseam mining and use the research results of deep mining for reference to calculate the failure of the multiseam mining roadway.

With regard to multiseam mining, a series of fruitful research has preliminarily revealed the stress distribution in multiseam mining. In theoretical research, a new approach was established using the pressure arch theory to calculate the pillar loads [22]. Sun et al. [23] studied the stability of overlying pillars using the theory of elastoplasticity and slip field. Xiao et al. [24] qualitatively explained the stress concentration caused by multiseam mining using the stress transfer theory. In terms of the numerical and physical simulations, Gao et al. [25], Bai et al. [26, 27], Kang et al. [28], and Ju et al. [29] applied numerical simulation and theoretical research to study the stress distribution of multiseam mining. Wang et al. [30] used UDEC simulation to obtain the stress release and stress concentration characteristics of multiseam mining. Ghabraie et al. [31], Tian et al. [32], and Huang et al. [33] used physical simulations to study the movement of roof strata and the instability of coal pillars during multiseam mining. Li et al. [34] studied the rock burst that was induced by driving an underprotective coal pillar using the seismic computed tomography technique and found that the concentrated stress generated by the overlying coal pillars shifted toward the depth and resulted in abnormal stress distribution. Shen et al. [35] used microseismic data to establish an energy density risk index to evaluate the risk of multiseam mining. Zhao et al. [11] analyzed several rock burst accidents and concluded that remaining overlying pillar induced rock burst by conveying high stress to the roadway below. These studies provide a reference for revealing the large deformation of the near-empty roadway under the overlying pillar in multiseam mining. However, the studies mainly lie within the regime of qualitative interpretation of stress concentration and instability in multiseam mining. The field practice shows that the deformation and damage to a roadway under the impact of overlying coal pillars depend on the coal seam spacing, goaf adjacent to the roadway, width of the overlying coal pillars, and certain other factors. The stress of roadway under goaf and coal pillar is very different. If the same support and pressure relief measures are adopted under the overlying coal pillars, it will often result in the wastage of human and material resources and would not allow effective control of the deformation and damage of roadway. Therefore, it is of significant theoretical and practical value to establish a quantitative calculation method to determine the stress concentration and roadway deformation and failure in multiseam mining.

Considering the width of the overlying coal seam goaf and coal pillar, the interval between coal seams, whether the adjacent working face of the roadway is mined out, the nature of coal and rock, and other factors, we established the model of “large deformation mechanism of the gob-side roadway under overlying coal pillars,” which can realize the quantitative calculation of the deformation and failure of the roadway in the near-empty area under the overlying pillar in multiseam mining. The research results will have important guiding significance for evaluating support effect and formulating reasonable support measures.

2. Engineering Background

2.1. Geologic Conditions of the Working Face

This study considers the 8308 panel in Xinzhouyao mine in Datong coalfield as the research object. Datong coalfield is a typical multiseam mining area with residual pillars. Figure 1 shows the location of the mine.

The 8308 panel is located in the east second panel of 14-3# coal seam and is 272–300 m deep. The roof and floor of the 8308 panel consist of hard sandstone. The distance between 14-3# and 11-2# coal seams is about 30 m. Figure 2(a) shows the synthesis columnar of 8308 panel.

The extraction method of 14–3# coal seam consists of fully mechanized longwall mining. The western side of the 8308 panel is 8310 gob, whereas the width of the interval coal pillar is 20 m. The eastern side consists of solid coal. The gob-side and nongob-side roadways contain retained bottom coal, which is 3.6 m wide and 3 m high. The overlying 11–2# coal seam, which is 30 m above the 14–3# coal seam, has been mined, and a large number of irregular pillars are left in this area. The gob-side roadway is adjacent to 8310 gob of 14–3# coal seam and passes through the pillar between 8210 gob and 8215 gob of 11–2# coal seam 337m-418 m area (R1 in Figure 2) away from the stopping line. It also passes through the pillar between 8215 gob and 8101 gob of 11–2# coal seam in the 503 m–678 m area (R2 in Figure 2(b)) away from the stopping line. The headgate passes through the pillar between 8208 gob and 8213 gob of 11–2# coal seam 161 m-175 m area away from the stopping line, whereas the pillar between 8213 gob and 8211 gob of 11–2# coal seam is in 266 m-400 m area away from the stopping line. Moreover, the pillar between 8215 gob and 8211 gob of 11–2# coal seam is in the 577 m-600 m area away from the stopping line. The spatial distribution relationship of the 8308 panel is shown in Figure 2(b).

2.2. Deformation Characteristics of the Roadway

The deformation of the gob-side roadway is greater than that of the nongob-side roadway in the roadway driving and mining of 8308 panel. The deformation and floor heave are particularly serious in R1 and R2 areas of the gob-side roadway. The two sides of the roadway move by more than 2.0 m, and the single props sink into the bottom coal up to 0.6 m. The deformation and failure of the gob-side roadway are shown in Figure 3. In contrast, although a certain degree of roadway deformation occurs in the area where the nongob-side passes over the overlying pillar, it is within the controllable range.

It is well known that, for the on-site roadway deformation condition, the deformation under the overlying coal pillars is much severer than that under the goaf, and the deformation in the gob-side roadway is severer than that in the nongob-side roadway. Roadway deformation under overlying coal pillars of multiseam mining is regionally different. The study of its deformation and realization of the quantitative calculations have certain guiding significances for preventing disasters under similar conditions. Therefore, considering the width of the overlying goaf and coal pillar, the interval between the coal seams, whether the adjacent working face of the roadway is mined out or not, the natures of coal and rock, and other factors, the model of “large deformation mechanism of the gob-side roadway under overlying coal pillars” can be established. The model can then be used to quantitatively calculate the loading and damage of the gob-side roadway under different types of overlying coal pillars.

3. Model of Stress and Deformation Calculations in the Gob-Side Roadway under Overlying Coal Pillars

Under the specific condition in multiseam mining, the stress concentration areas that are affected by the overlying pillar convey stress to the gob-side roadway according to the variation pattern of stress transfer. Meanwhile, the cantilever roof is formed in the adjacent goaf, which forms the lateral abutment pressure on the surrounding rock of the gob-side roadway. According to the theory of elastic mechanics, the coal rock mass under the overlying coal pillar is regarded as an ideal elastic body, whereas the stress distribution on the overlying coal pillar is simplified as a uniform load. With this backdrop, the model of “large deformation mechanism of the gob-side roadway under overlying coal pillars” is established, as shown in Figure 4. Moreover, , , and are the stress values of the overlying coal pillars, whereas , , and are the stress values of the overlying goafs.

3.1. Stress Calculation and Stress Transfer in the Overlying Coal Seam

The loads on the pillars that are leftover from the overlying coal seam mining are the superposition of the self-weight of overlying strata on the pillars and the lateral abutment pressure produced by the overhanging roof in the goaf area of one or both sides of the pillars. After the overlying coal seam is mined, the load on goaf and pillars tends to evenly distribute after a long rebalancing of coal and rock structures. The unit length load of the coal pillar can be estimated according to the load estimation model of the coal pillar, as shown in Figure 5 [36].

Assuming that the width of goaf on two sides of the coal pillar is different, the average load on coal pillars can be given by

Here, D1 and D2 are the widths of goaf on two sides of coal pillar; is average stress of coal pillar; B1 and B2 are the widths of two different coal pillars; H is the depth of coal seam; and is roof caving angle.

The difference between the self-weight stress of overlying rock in coal seam and the load of remaining coal pillars is the load of goafs. Then, the average load of goafs can be written as

After obtaining the stress distribution of coal pillar and goaf in the overlying coal seam, the stress transfer model of overlying coal pillars and goafs is established. As shown in Figure 6, the X-axis is in the vertical direction, while the Y-axis is downwards.

After the mining of the lower coal seam, the roof will produce a caving zone, fissure development zone, and bending subsidence zone. Through calculation, the height of the caving zone and fissure development zone is greater than the interval between the two coal seams, so the overlying coal pillar 1 loses the ability of stress transfer to point M. Therefore, the stress transmitted by is not calculated here. Taking the AB segment pillar as an example, the stress transferred from the overlying coal pillars and the goafs to point M is calculated. The coordinates of points A, B, and M are (xA, yA), (xB, yB), and (x, y), respectively. The infinitesimal length is taken at the distance from the origin of the coordinate in the segment AB. According to the theory of elasticity, the stress caused by the element at point M is given by

By integrating each infinitesimal stress in the segment AB using equation (3), and based upon to , the vertical stress transferred from segment AB to point M is given by

Equation (4) is integrated to obtain

The stress of the overlying residual coal pillars and goafs transferring to point M within the influencing range of the lower section of the coal pillar is calculated using equation (5). Moreover, the stress value of the overlying residual coal pillars and goafs transferred to the coal pillar section can be determined.

3.2. Calculation of the Lateral Abutment Pressure

After mining out the adjacent working face, the roof forms a cantilevered block on the district sublevel pillar (coal pillar between the next working face and goaf of the last working face), causing lateral abutment pressure. As shown in Figure 7, the angle between the roof fracture line and the horizontal line after the collapse of the roof is called the roof caving angle. The angle between the initiation line of the roof’s rotational deformation and the horizontal line is called the bedrock movement angle.

The self-weight of the rock is considered in the calculation of stress transfer of overlying coal pillars and goafs. Therefore, when calculating the stress of the cantilever roof in the adjacent goaf, only the lateral abutment pressure between the two layers of coal is considered.

The lateral abutment pressure () consists of the gravity stress () and the additional stress increment () transferred from the cantilever roof around the gob-side roadway and is given by

, , and are the i-th increment of additional stress transferred from overburden rock block along goaf to lateral coal body, and m is the number of roof cantilever blocks between the lower coal seam and overlying coal seam.

The cantilevered block is approximately calculated as an elastic body. The additional loads are distributed in the form of uniform loads on the surface of the bearing rock strata. Since only the strata between the two coal seams are considered, and the interval between them is small, the force of the cantilevered block is calculated according to the structure without any articulation. The additional load generated by the cantilevered block is , where is the sum of the weight of the cantilevered block in the i-th layer and the weight of the cantilevered block between the roof of the i-th layer and the overlying coal seam. is the vertical distance from the initiation point of the cantilevered block’s rotational deformation in layer i to the coal wall.

Taking the second-layer cantilever block as an example, the stress of the additional load of the single-layer cantilever block transferred to point M is calculated. The coordinates of point C at the intersection of the coal wall extension line and the upper boundary of the second-layer cantilever block are (xC, yC), whereas the coordinates of point D at the intersection of the bedrock moving line and the upper boundary of the second-layer cantilever block are (xD, yD). If the length of the infinitesimal is at the distance from the origin of the coordinates of the segment CD, the stress caused by the infinitesimal at point M is given by

Here, q2 is the additional load generated by the second cantilevered block.

Integrating the stresses of each infinitesimal in the segment CD from to (equation (7)), the stress value of the additional load of the second-layer cantilever block transferred to point M is given by

Equation (7) is integrated to obtain

The surrounding rock transferring stress around the gob-side roadway, transmitted by the cantilever roof around the gob-side roadway, is given by

According to equation (9), the stress of the cantilevered block in each layer at point M can be calculated. The lateral abutment pressure at point M is obtained by combining the calculation results with those of equation (10).

3.3. Calculation of the Plastic Zone Range

The large deformation of the roadway is induced by the formation and development of the plastic zone of coal and rock around the roadway. The scope and distribution of plastic deformation directly represent the size and shape of the deformation of the roadway. Therefore, it is necessary to further study the size and distribution of plastic deformation after obtaining the superposition value of the stress transmitted by the overlying coal pillars and the cantilever roof around the gob-side roadway. Cheng [37] discussed the adaptability of the correlation of the plastic zone to the surrounding rock of the roadway. In their study, the range of the plastic zone is calculated using the Kastner correlation. In the calculations, the roadway is simplified as a circular one. According to the stress calculation model of the gob-side roadway under the overlying coal pillars, the stress value is calculated and combined with the results of the Kastner correlation to calculate the plastic zone range of the roadway surrounding the rockwhere is the radius of the outer circle of the roadway, is the depth of the plastic zone, is the coefficient of lateral pressure, is the cohesion of the surrounding rock, is the angle of internal friction, is the support resistance, and is the superposition value of the transfer stress of the overlying coal pillars and the cantilever roof around the gob-side roadway.

The above process constitutes a model for quantitatively calculating the stress and deformation of the gob-side roadway under overlying coal pillars. The model fully considers the role of overlying coal pillars, goaf, and the cantilever roof around the gob-side roadway and calculates the range of the plastic zone of gob-side roadway based on multistress superposition results.

4. Stress and Plastic Zone Analysis Based on the Model

4.1. Calculation of Stress and Plastic Zones
4.1.1. Stress Transfer Calculation of 11–2# Coal Seam

From the derived stress transferring model and specific calculation results, the stress value of the remaining pillars and goafs of 11-2# coal seam transferred to lower coal seam is calculated. The width of the coal pillars and goafs in Section I-I, as shown in Figure 1, is measured at 37 m, 51 m, and 44 m (in width) from left to right. The width values of goaf successively are 58 m and 100 m. The buried depth of 11–2# seam is 272 m. The roof caving angle is given by , and the bedrock movement angle is given by . The distance between the 11–2# seam and 14–3# seam is 30 m.

According to equation (5), the stress values of the residual pillars and goafs in 11–2# coal seam near the gob-side roadway are calculated. Because the 8310 panel has been mined out, the stresses of the residual pillars in the upper part of the 8310 gob are not calculated. Then, the stress transfer superposition values of a certain point near the gob-side roadway are obtained, which are given by

4.1.2. Calculation of the Lateral Abutment Pressure

According to the lateral abutment pressure calculation model, the pressure produced by the cantilevered block in 8310 gob is calculated to be kN/m3, m.

When calculating the lateral pressure, only the one produced by the roof between the mining seam and its upper seam is considered. According to the borehole histogram of the 8308 panel, the roof is divided into four layers, and the thickness values are 3 m, 4 m, 6 m, and 17 m from bottom to top. The lateral abutment pressure of 8310 gob can be obtained by incorporating the parameters into equation (10), which results in

4.1.3. Calculation of the Total Stress

By superimposing the stress transfer value of 11–2# coal seam and the lateral abutment pressure, the vertical stress distribution curve of 14–3# coal seam is obtained, as shown in Figure 8.

4.1.4. Calculation of the Plastic Zone in the Gob-Side Roadway

Because of the influence of the overlying coal pillars and the lateral abutment pressure produced by the cantilever roof around the gob-side roadway, the stress around the gob-side roadway is not uniformly distributed. According to the results of the mechanical model, the stress values in the 20 m range on the left and right sides of the central line of the gob-side roadway are integrated. Then, the stress values on the left and right sides of the roadway are obtained as follows:where is the stress value on the left side of the roadway midline, ΦL is the integral value of stress within 20 m on the left side of the roadway midline, is the stress value on the right side of the roadway midline, and ΦR is the integral value of stress within 20 m on the right side of the roadway midline.

After obtaining the vertical stress on both sides of the gob-side roadway, the plastic zone range of the roadway is calculated using equation (11). For this purpose, various parameters have the following values: radius of the outer circle of the roadway , cohesion  MPa of the surrounding rock of the coal body, cohesion MPa of the surrounding rock of roof and floor sandstone, friction angle of coal body, and internal friction angle of the roof and floor sandstone. The role of the supporting resistance is not considered for the time being.

The vertical stress values of the roadway’s two sides are incorporated into the plastic zone to calculate the radii of the roof and floor and the coal body. The left and right sides of the roof and floor are 1.8 m and 1.72 m, respectively, while the two sides of the coal body reach the values of 4.73 m and 4.29 m.

The calculation boundary is shown in Figure 9.

4.2. Comparative Analysis with Numerical Simulation

Based on the actual geological and mining conditions of the 8308 panel, FLAC3D is used to simulate the surrounding rock stress of the gob-side roadway in the 8308 panel and compared with the theoretically calculated results. The rock strata were modeled using the Mohr–Coulomb model. The top boundary condition of the model is stress constraint, while the load equivalent to the self-weight of the overlying strata (6 MPa) is applied on the top. According to the in situ stress test results, the lateral pressure coefficient is set to be 1.2. The bottom surface is constrained by normal and tangential displacements. The thickness values of 14–3# coal seam, floor, and roof are 4.5 m, 20 m, and 30 m, respectively. The thickness values of 11–2# coal seam and the roof are 5 m and 15 m, respectively. In order to avoid the boundary effect, the geological conditions of the 8308 panel are extended around the model. The left and right sides extend by 100 m and 50 m, respectively. The front and back extend by 50 m each. The mechanical parameters of 14–3# coal seam and bottom roof and 11–2# coal seam and roof are presented in Table 1. After mining the working faces of 11–2 # and 14–3 # coal seams, the constitutive model of the goaf immediate roof is changed to the strain-softening model. When the goaf immediate roof does not yield, the tensile strength is 0, and the internal friction angle and cohesion remain unchanged. When the plastic yield degree of direct jacking is 1%, the cohesion decreases to 14.11% and the internal friction angle decreases to 50.87%. When the plastic zone is completely formed in the immediate top, the cohesion is reduced to 4.70% and the internal friction angle is reduced to 16.14%. The numerical simulation model is shown in Figure 10.

According to the model, the plastic zone and vertical stress of 14–3# coal seam at Section I-I are obtained, as shown in Figure 11.

Comparing the results of the mechanical model and the numerical simulation at Section I-I, the numerical simulation results show that the depth of the plastic zone on both sides of the gob-side roadway is about 4.2 m each, whereas the corresponding theoretically calculated value is 4.29–4.73 m. The correlation coefficient of stress calculation results is 0.89. The numerical simulation results show that the roof caving angle and the bedrock movement angle are about 75°. This is consistent with the field measurement and the assignment of the theoretical model. The results of the numerical simulation for the plastic zone, roof collapse, and stress are highly correlated with the results of the mechanical model. Therefore, both the theoretical model and numerical simulation can be used to calculate the plastic zone of the gob-side roadway under the overlying coal pillar.

5. Discussion on the Deformation Mechanism of the Gob-Side Roadway

By analyzing the mechanical model of “stress and deformation quantitative calculation of gob-side roadway under overlying coal pillars” and the calculation results of the 8308 panel as an example, the main influencing factors of roadway’s large deformation can be obtained, and the degree of the influence of each main influencing factor can be quantitatively analyzed. Combined with the existing support measures, it is concluded that the main influencing factors of large deformation in the gob-side roadway under overlying coal pillars are as follows.

5.1. Overlying Coal Pillars

As shown in Figure 2(b), the gob-side roadway of the 8308 panel is located directly below the overlying coal pillar at Section 1. After the 11–2# coal seam has been mined out, the stress transfer results of the overlying coal pillars and goafs are obtained from the model, as shown in Figure 12. The biggest contribution to the stress comes from the coal pillar above the roadway. The peak stress reaches the value of 12.06 MPa, and the stress transmitted by the coal pillar upright above is 78.3%. The peak stress value is 1.8 times of the self-weight stress of overlying strata. The average stress transmitted by the pillar above the roadway is 81.5% of the total average stress.

5.2. Goaf Adjacent to the Roadway

Figure 13 shows the lateral abutment pressure of Section 1 obtained by the model. The peak stress value of the lateral bracing pressure is 1.93 MPa. The peak stress position is shown in Figure 13, which accounts for 16% of the total vertical stress.

Due to the 8310 gob on the left of the 8308 panel, the horizontal stress of the coal pillar in the section decreases. Zhao [38] showed that, under certain other conditions, the smaller the lateral pressure coefficient, the greater the depth of the plastic zone on both sides of the roadway is. Therefore, the roadway adjacent to goaf is more prone to large deformation.

5.3. Bottom Coal Retained in the Roadway and Insufficient Support

The roadway of the Xinzhouyao coal mine adopts a rectangular layout, as shown in Figure 14. The roof of the roadway is supported by four rows of bolts, which are 20 mm in diameter, 2 m in length, and 0.8 m × 1.04 m in row spacing. The construction of the two rows of roof anchor cables is 1.82 m away from both sides of the roadway, while the row spacing is 3 m. The two sides of the roadway support the three rows of bolts, while the row spacing is 1 m × 1.2 m. The bottom coal is retained in the roadway. The single prop support is adopted within 50 m of the advanced working face.

Comparing the plastic zone calculated using the theoretical model (Figure 9) with the actual support parameters on-site, the length of the roadway bolt is found to be only 1.7 m, whereas the bolts are all in the plastic zone. It is difficult to achieve a substantial roadway support effect.

The established model can quantitatively calculate the stress and the degree of deformation under the comprehensive effects of gob and overlying coal pillars. Furthermore, the study also analyzed the degree of influence that various factors had on the stress and degree of deformation. The research results will have important guiding significance for evaluating the support and the formulation of related measures.

It is worth noting that the model and numerical simulation will simplify the engineering geological conditions in the calculations. Therefore, there will be some errors between the theoretical and simulation results and the field situations. Due to this reason, it is necessary to combine the results with the on-site data to determine the damage assessment and support scheme formulation of the gob-side roadway under the overlying coal pillar.

6. Conclusions

(1)We established the mechanical model of “stress and deformation quantitative calculation of the gob-side roadway under overlying coal pillars.” According to the established model, we studied the influence degree of the main control factors on the failure of the gob-side roadway in Section I-I. The peak value of the stress around the gob-side roadway reaches 12.06 MPa. The gob-side roadway’s overlying pillar transmits 78.3% stress. The lateral abutment pressure transmits 16% stress. The peak value of the stress is 1.8 times of the self-weight stress of overlying strata.(2)Through the quantitative calculation of the stress and plastic zone of the gob-side roadway under the overlying coal pillars, the main influencing factors of the large deformation of the gob-side roadway are obtained: the stress concentration caused by the overlying coal pillars, the lateral abutment pressure produced by the cantilever roof, the bottom coal retained in the roadway, and insufficient support.(3)The established model can quantitatively calculate the stress and plastic zone distribution of the roadway under the action of multiple force sources and obtain the influence degree of the main control factors. The obtained results have an essential reference significance for evaluating support effect and formulating reasonable support measures of the gob-side roadway under overlying coal pillars in the multiseam mining.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest in this work.

Acknowledgments

This work was financially supported by the National Natural Science Foundation of China (grant nos. 51634001 and 51774023), the State Key Research Development Program of China (grant no. 2016YFC0801408), and the Major Science and Technology Innovation Project of Shandong Province (2019SDZY01). The authors would like to express their gratitude to all the agencies for funding this research.