Abstract
Support crushing and roof fall accidents frequently occur in the large mining height stope covered by the thin bedrock and large-thick alluvium (TBLA). A comprehensive method including field monitoring, numerical modelling, and theory analysis was performed to investigate the strata movement and its induced ground behaviors of a deep-buried stope under TBLA. The filed measurement results indicated that strong roof periodic weighting phenomenon occurred regularly at an interval of 18–23 m during panel 3301 retreating. During weighting period, dynamic characteristics of hydraulic supports are prominent and characterized by high safety valve opening rate and poor working conditions including tilt down operation and excessive end distance. Additionally, numerical modelling results revealed that the thickness of bedrock affects the movement behaviors of the roof strata and its generated spatial structure significantly. Based on the above results, an integrated control measure was proposed and successfully applied to the test site.
1. Introduction
Coal resources are rich in reserves and have obvious regional distribution characteristic in China. East China is an important coal energy-producing area in China, in which the development and utilization of coal resources play a decisive role in promoting the economic rapid development [1–4]. Currently, the coal resources covered by the shallow alluvium are nearly exhausted in the East China. As a result, many mines are exploiting the coal resources covered by the large-thick alluvium of the Quaternary system. Compared with the coal mines in Western China, coal mining in East China encounters new challenges, such as large burial depth, thin bedrock, and large-thick alluvium [5–8]. Consequently, many disasters, such as water inrush, support crushing, and coal mass leakage at the end face, occurred in the field. Take Qidong Coal Mine in Anhui Province as an example; a total of 17 times water inrush and support crushing accidents occurred in the panels such as 3222, 7114, and 6130. Therefore, water inrush occurred in the panel 3222; it even flooded the mine, resulting in a restoring fee of nearly 1 billion.
To date, a number of studies have been conducted to investigate the strata movement and its induced ground behaviors of stopes covered by TBLA in China. For example, Xu et al. [9] studied the load transfer mechanism of loose aquifers and its influence on the composite failure of key stratum by a physical simulation experiment. Fang et al. [10] investigated the influence of the thickness and mechanical properties of the loose layer on the ground response of the stope. Zhang et al. [11] revealed the correlation between the spatial structure and ground behaviors of the roof strata in a deep-buried stope covered by TBLA in Gubei Coal Mine. Du et al. [12] explored the quantitative relationship between the working resistance of the supports and the roof subsidence in a fully mechanized caving mining stope covered by thin bedrock. Considering the influence of arch structure formed in alluvium on the bedrock, Wang et al. [13] presented a calculation formula for the working resistance of the hydraulic supports in the stope covered by alluvium. Zhang et al. [14] established the strata structure model and deduced the height calculation formula of the fracture zone based on geological conditions of Xieqiao coal mine. The above studies provided a better understanding of the roof strata movement principle and corresponding ground control in the stope covered by TBLA. However, according to the literature review of these studies, it can be found that the studies have two major limitations. The first is that the above studies are mostly focused on Huainan, Huaibei, and Yongcheng mining areas, where the thickness of loose layer, coal seam, and bedrock are 200–400 m, 3.0–6.5 m, and 18–43 m, respectively. However, there are few case studies on the ground behaviors of a large mining height stope, which has the characteristics of large buried depth, thin bedrock, and thick alluvium in Juye mining area. The second is that the previous studies are performed mainly by a simple simulation method or small-scale short-term field monitoring. Due to the complexity of the geological conditions, the studies on the ground behaviors of a specific stope should be implemented by an integrated method of all-around field monitoring, numerical simulation, theoretical calculation, and other methods.
Juye mining area is located in the southwest of Shandong province, China, with proven geological reserves of 5.57 billion tons. A total of seven pairs of coal mines in this mining area generate 18 million tons coal every year. Different from other mining areas such as Huainan, Huaibei, and Yongcheng, the coal seam in Juye mining area is 5.8–9.9 m, which is covered by more than 600 m thick alluvium and 40–150 m thick bedrock. The above complex geological condition resulted in many difficulties in the coal mining activities. Various accidents such as support crushing and end face leaking roofs occurred frequently in the field. According to incomplete statistics, more than ten times large-scale roof falling and support crushing accidents occurred in the Juye mining area in 2018 and 2019. Therefore, it is urgent to conduct in-depth and systematic research to investigate the strata movement and its induced ground behaviors of the stopes in the Juye mining area.
The objective of this paper is to investigate the roof strata movement process and its induced ground behaviors of a stope covered by TBLA in the Juye mining area. This paper is organized as follows. In Section 2, the detailed geological and production conditions are introduced. In Section 3, the ground response of the stope, including the working resistance and geometric position of the hydraulic supports, roof leakage, and the coal wall slip, was monitored and analyzed systematically. In Section 4, the distribution and evolution characteristics of roof strata convergence and stress changes during coal extraction stages were presented through numerical modelling. In Section 5, the movement process of the roof strata movement associated with its induced ground behaviors of the stope was revealed to determine a reasonable working resistance of the hydraulic supports. In Section 6, the validity of the newly control measures is verified by a field application. The roof strata movement process and its induced ground behaviors in this study can provide some reference for other similar mining TBLA projects.
2. Engineering Background
Guotun Coal Mine is a representative mine that has the most significant change in bedrock thickness in Juye mining area, and large-scale roof falling and support crushing accidents occurred many times. Guotun coal mine has certain representative, so it was selected for this case study. Its geographical coordinates are 115°50′-116°0′0″ east longitude and 35°27′00″-35°34′30″ north latitude. The mine field is about 16 km long from north to south and 14 km wide from east to west, with an area of 222.1 km2, and the occurrence area of the main mining #3 coal seam is 89.8 km2. The mine covered by the Tertiary and Quaternary of the Cenozoic is a completely concealed North China type Carboniferous Permian coalfield. Generally, the thickness of the Cenozoic is 530–650 m, with an average of 590 m. The main seam (#3 coal seam), which has an average thickness of 7.02 m and an average dip angle of 5.5°, has stable occurrence and simple structure. The roof with great thickness variation, which has a false roof of 0–3.13 m in some parts, is dominated by fine and siltstone, and the floor is dominated by mudstone and siltstone.
The study site is panel 3301, which is the first mining face of the three mining areas with a dip length of 188 m and a strike length of 1212 m. A total of 111 hydraulic supports of model ZY13000/28/63D are applied to basic support, in which the pressure of the relief valve and the pump station is 46.9 MPa and 31.5 MPa, respectively. The section map of panel and geological column chart of panel 3301 is illustrated in Figure 1. According to the investigation results of multiple geological boreholes, it can be revealed that the thickness of bedrock of panel 3301 begins to increase gradually since the position of open-off cut has a range of 35–110 m.
Panel 3301 adopts the comprehensive mechanized coal mining technology of full height at one time by the caving method to control roof. Additionally, it is the first mining face of the three mining areas with a tendency length of 188 m and a strike length of 1212 m. A total of 111 hydraulic supports of model ZY13000/28/63D are applied to basic support, in which the pressure of the relief valve and the pump station is 46.9 MPa and 31.5 MPa, respectively.
3. Field Investigation
Ground pressure can reflect the overlying strata movement of the stope and provide basic data for subsequent roof strata spatial structure analysis. In this section, field measurements including the working resistance and geometric position of the support, roof leakage, and the rib slip are conducted to obtain detailed ground behaviors law during panel 3301 retreating. The monitoring lasted from June 2, 2019, to December 23, 2019. Twenty YHY.60 stress meters were installed to monitor the working resistance of the support [15], the monitoring stations layout, and the hydraulic support resistance monitoring equipment as shown in Figure 2. Because the monitoring results in close region station have similar tendencies, only the measurement data from the station (#95, #55, #15) are used for the analysis.
3.1. Analysis of Roof Weighting Law
It can be observed in Figure 3 that #15 support suffers three times roof weighting events with an average weighting step of 18.5 m. The average working resistance of the hydraulic support during weighting period is about 11035 KN, while that in the normal is about 7621 KN. Thus, the corresponding dynamic load coefficient of the support is 1.44. For the #55 support, two times weighting events occurred with an average weighting step of 23 m. The working resistance of the hydraulic support in the weighting period and the normal is 11565 KN and 7021 KN, respectively, with a dynamic load coefficient of 1.64. For #95 support, three times roof weighting events occurred with an average weighting step of 17.5 m. The working resistance of the hydraulic support in the weighting period and the normal is 12278 KN and 7786 KN, respectively, with a dynamic load coefficient of 1.57.
(a)
(b)
(c)
The working resistance was divided into low resistance interval (0∼5200 KN), normal resistance interval (5200∼10400 KN), and high resistance interval (10400∼13000 KN), which are 0∼40%, 40%∼80%, and 80%∼100% of the rated working resistance, respectively. The distribution characteristics of the working resistance of #15, #55, and #95 can be obtained. For the #15 support, the proportion of the working resistance in the low, normal, and high resistance interval is 17.12%, 56.12%, and 26.74%, respectively. For the #55 support, the proportion of the working resistance in the low, normal, and high resistance interval is 40.62%, 45.08%, and 14.28%, respectively. For the #95 support, the proportion of the working resistance in the low, normal, and high resistance interval is 43.33%, 46.89%, and 9.76%, respectively.
It can be observed clearly in Figure 4 that the supports in the lower-middle area were in a low resistance working state for a longer time than that in the upper area, which may be attributed to the diversity of the roof strata conditions. The roof strata in the stope upper area were strong siltstone with a large thickness, which generated more vertical loads to the hydraulic support, resulting in the supports in a high working resistance state. However, roof strata in the stope lower-middle area were a multilayer roof, which is composed by weak mudstone or siltstone with a small thickness. Due to the lower antideforming capability of the multilayer roof, the roof separation was prone to occurrence and further caused the vertical loads transferred to the nearby supports. As a result, the phenomenon of “one support in low working resistance and its adjacent supports in high working resistance” occurred frequently in the stope lower-middle area [16].
(a)
(b)
(c)
3.2. Support Position Information
A comprehensive information of the geometric position of the hydraulic supports, including pitch angle of a canopy and the support end distance, was performed during panel 3301 retreating. The pitch angle refers to the angle canopy deviating from the horizontal plane, and the end distance refers to the distance between the front end of the canopy and the coal rib.
Figure 5 depicts the statistics of pitching angle of top beam and end face distance during hydraulic supports service period. Based on the pitch angle characteristics of the canopy, the working condition of the hydraulic support can be divided into tilt-up and tilt-down working condition, which is with pitch angles of 0–5° and 0–10° from the horizontal plane, respectively. It can be observed in Figure 5(a) that the hydraulic supports in high-efficiency and low-efficiency working condition account for 59% and 41%, respectively. When the supports are in low-efficiency working condition, the coal and rock mass in the end face are unstable and prone to falling. This is because the supports in the low-efficiency condition reduced the contact area between the canopy and immediate roof, obstructing the working resistance transferring from the supports to the immediate roof. As a result, the immediate roof showed subsidence and further generated more loads to the hydraulic support. In turns, the supports working condition deteriorates further.
(a)
(b)
The statistics of support end distance in panel 3301 are presented in Figure 5(b). About 61% of the supports are working good; therefore, about 57% of the supports have an end distance less than 0.3 m, and 39% of the supports have an end distance of 0.3–0.8 m. It should be noted that about 4% of the supports are working poor with an end distance greater than 0.8 m, which is caused by the tilt down operation of some hydraulic supports and the supports following-up untimely after coal mass mined.
3.3. Analysis of Rock Mass Caving Accidents in the End Face
The unreasonable position of the support will reduce the stability of the coal and rock mass in the end face, which may lead to rib spalling and roof falling accidents. In the field, many rib spalling and roof falling accidents occurred during the panel 3301 retreating. Figure 6 shows a rib and roof falling accidents in panel 3301; it can be found that the rib and roof falling accidents are frequent, and the breaking degree is small.
For the rib spalling, the depth of the rib spalling in 0–0.25 m and 0.25–0.5 m accounts for 74% and 12%, respectively, while the depth of the rib spalling in 0.5–1.0 m accounts for 14%. For the roof falling, the depth of the roof falling was concentrated in the range of 0.25–0.5 m, accounting for 55%, while the roof falling depth within 0.5–1.0 m accounts for 20%.
It can be found from the above accidents that the roof falling accidents mainly occurred at the false roof of mudstone or siltstone; see Figure 7. The roof falling process in panel 3301 is illustrated in Figure 8. Under the effect of the front abutment pressure, massive cracks developed in the roof mudstone or siltstone immediately above the virgin coal mass. As the panel retreated and supports advanced, the immediate roof was above the hydraulic supports. Due to the repeated support effect of the hydraulic supports, the immediate roof turned to be damaged seriously and even falling. In this situation, the supports cannot achieve their full capacity, which in turn resulted in rib coal mass unstable and even spalled. Meanwhile, the abnormal support geometric position, including the supports following-up untimely and tilt down operation, was also an indispensable reason for the roof falling and rib spalling [17].
4. Numerical Model Generation
4.1. Numerical Simulation Model and Simulation Plan
The characteristics of stress and displacement of stope surrounding rock were investigated during the panel retreating by a numerical simulation model. The dimensions of the model were 200 m × 100 m × 120 m (see Figure 9), which were determined based on the geological and production conditions of the panel 3301. Based on filed measurement results, a vertical stress of 15 MPa was applied to the model upper boundary to simulate an overburden of weight, while the horizontal-to-vertical stress ratio is set to 1.3 in the x-and y-directions. Note that there are many bedding, discontinuities, and joints in coal or rock mass, which affect their failure characteristics. Although finite difference method using FLAC3D is incapable of truly capturing their characteristics, coal or rock mass properties that are upscaled from intact cores properties using the Hoek-Brown criterion can be performed in the model. In other words, the coal and rock properties in our model are effective properties that accounted for the rock discontinuities. The coal and rock mass parameters were modified by criterion of Hoek-Brown [18–21], as shown in Table 1.
(a)
(b)
In this simulation, the Mohr–Coulomb model was used for the rock strata, and the numerical model was solved using the following steps: (1) establishment of numerical model, (2) calculation of original rock stress balance, and (3) mining in panel 3301 (the panel with mining height of 6 m has been excavated for 100 m in total and ten times). Meanwhile, due to the prominent variation of bedrock thickness in panel 3301 (35 m–110 m), the influence of bedrock thickness on the overlying strata movement rule of stope was studied. And combined with filed investigation results at panel 3301, the characteristics of roof stress and displacement filed evolution under the bedrock thickness (43.5 m, 48.5 m, 53.5 m, 68.5 m, and 94 m) were analyzed.
4.2. Analyses of Model Results
4.2.1. Vertical Displacement Analysis
Investigations have shown that the exposed area of the roof continues to expand as panel retreating. In the action of underground pressure, immediate roof was bent, showed subsidence, and collapsed successively, and the main roof subsequently showed periodic breaking. The displacement distribution of immediate and main roof is shown in Figures 10 and 11. ① When the panel retreated from 10 m to 20 m, the subsidence of immediate and main roof is 187–453 mm and 494∼822 mm, respectively, and the immediate roof and main roof have showed separation; ② when the panel retreated further, from 30 m to 40 m, the subsidence of main roof rapidly increased to 758–1086 mm, and the main roof was first breaking during first weighting. It can be inferred that the main roof is instable due to the sliding of structure, which results in bench convergence, roof cutting, and support crushing; ③ with further retreating of the panel, from 50 m to 60 m, the subsidence of main roof gradually increased to 1370–1597 mm, and main roof was breaking during periodic weighting. Compared to first breaking, the increased rate of main roof subsidence is small, and ground behaviors are relatively mild, which is almost consistent with the field measured results.
(a)
(b)
Based on the above results, it can be inferred that the subsidence of the immediate roof is significantly larger than that of the main roof during panel retreating. It should be noted that the main roof will show instability and sliding when the separation deformation is large [22], which generated more dynamic loads to the hydraulic supports and led to the pressure of the supports suddenly increasing or even support crushing.
4.2.2. Vertical Stress Analysis
Figure 12 shows the distribution characteristics of roof stress during panel retreating. The roof stress distribution has obvious zoning characteristics (caving zone and fracture zone). For rock mass in the caving zone, it has tensile failure and tensile stress relief zone, whereas, for the rock mass in fracture zone, it is in compression and forms a compressive stress zone. Further, the interface between two stress zones is approximately arched, and the arch foot is located in the rib and the coal of the open-off cut. The high stress areas on both sides of the goaf generated 70° with the horizontal direction, which extends outward to the top of the model; that is, the displacement angle of stratum is about 70°.
Figure 13 shows the vertical stress curves of immediate roof and main roof, and Figure 14 shows vertical stress contours and vectors of main roof. The abutment pressure of deep-buried stope covered by TBLA is characterized by high pressure, large dynamic load coefficient, and small influence range. The peak value of the vertical stress of the immediate roof is 22.5 MPa, and stress concentration factor (kp) is 1.5 during the first weighting of the panel, which is at the position 6.7 m away from the panel. During the normal retreating period of the panel, the peak value of the vertical stress of the immediate roof is 23.61 MPa, and stress concentration factor (kp) is 1.57, which is at the position 7.5 m away from the panel. It is noted that the above stress distribution characteristics are closely related to the sliding instability of strata and small size of broken rock block.
(a)
(b)
4.3. Strata Movement Law at Various Bedrock Thickness
Figure 15 illustrates the changes of stope roof subsidence at various bedrock thickness. It can be observed that the roof subsidence decreases continuously and presents periodic change regularities with the increase of the thickness of bedrock. When the thickness of bedrock is less than 53 m, the roof subsidence is significant and changes at a large rate. As the thickness of bedrock decreases from 53.5 m to 43.5 m, the maximum displacement is 2.46 m, 3.50 m, and 5.40 m, respectively. When the thickness of bedrock is greater than 53.5 m, the roof subsidence tends to be stable gradually, and the largest subsidence had a little fluctuation, in the range of 2.27∼2.46 m. It is well known that the roof displacement is related to the overlying strata movement. Based on ground pressure theory and the field pressure result, the strata movement law of roof could be summarized as follows: ① for the bedrock with thickness greater than 53.5 m, it has the characteristics of high bearing performance capacity and can form a relatively stable balance of roof structure; thus, the roof subsidence is small and stable. ② For the bedrock with thickness less than 53.5 m, due to the small thickness and low bearing performance of bedrock, the bedrock will slip and become unstable [23], leading to the roof producing large vertical displacement in a short time, further resulting in strong ground pressure such as roof cutting and support crushing on-site [24–27].
5. Rock Stability Control of the Thin Bedrock
5.1. Mechanical Principles of the Thin Bedrock
According to the analysis of geological conditions, ground response behavior law, and roof strata spatial structure during 3301 panel retreating, the mechanical principles of stability of the thin bedrock could be summarized as follows: ① when stope is covered by the thin bedrock and large-thick alluvium, the main roof mainly showed sliding instability, and the broken rock block further generated large dynamic loads to immediate roof and hydraulic support, which will lead to strongly ground pressure phenomena including roof falling and support crushing. ② Panel 3301 adopts the comprehensive mechanized coal mining technology of full height at one time, and the height is 6 m. Due to large mining height, weak strata, and waste rock in seam of coal mass, the coal wall is unstable and prone to falling. As a result, support end distance increases rapidly and further leads to the area of roof falling. In turn, the rib slip deteriorates further. ③ When support end distance is too large, which is caused by immediate roof cutting, inaccurate operation, and tilt down of hydraulic supports, only improving hydraulic support horizontal force cannot ensure the stability of coal and rock mass in end face. Once the rock mass in the end face showed leakage, the immediate roof turned into seriously falling, resulting in hydraulic supports being unable to provide greater support resistance, and the relationship between surrounding rock and hydraulic support deteriorated.
5.2. Control Strategy of the Thin Bedrock
Based on the above analysis, it can be inferred that, in the large mining height stope covered by TBLA, the key technology to control the stability of surrounding rock adjusts the operating mode of hydraulic supports to reduce dynamic load influence caused by sliding and instability of roof strata. The main measures on adjusting the operating mode of hydraulic support are as follows:(1)Reasonable support working resistance: according to analytical and numerical modelling results, it can be inferred that roof strata spatial structure type depends on the thickness of bedrock, and the hydraulic support needs to bear more loads when bedrock thickness is small. During panel 3301 retreating period, the thickness of bedrock begins to increase gradually since the position of open-off cut is with a range of 35–110 m, and the rock strata combination tends to be complicated. Correspondingly, the roof strata spatial structure gradually transits from single key layer structure to multikey layer structure, and the stability of roof strata structure is improved gradually. Therefore, in actual engineering practice, for the purpose of achieving a coal mine safety goal, the hydraulic supports with abundant working resistance should be selected.(2)Ensure sufficient initial support force: ① based on the analysis of the field results of ground pressure, it can be obtained that when the supports have sufficient initial support force, the coal and rock mass in the end face are stable. This is because the supports owned sufficient initial support force increased the contact area between the canopy and immediate roof, avoiding the dangerous of supports failure. ② In order to avoid roof cutting and support crushing induced by sliding instability of roof strata, the height of hydraulic support legs should be increased to maximum before the rock strata breaks.(3)Improving the management level of hydraulic supports and strengthening the fault diagnosis capability of hydraulic supports: ① during panel 3301 retreating period, through improving the proficiency of workers in support operation, ensure the support end distance and pitch angle of a canopy should be ensured in 0.7 m and 5°, respectively. ② Based on a field observation of panel 3301, it can be obtained that the support working resistance reduces rapidly when the hydraulic system is damaged. Therefore, it is necessary to carry out timely fault diagnosis and ensure panel retreating rapidly in weighting period, further avoiding excessive working resistance of the hydraulic supports due to stagnant of the panel.
5.3. Field Application and Implications
5.3.1. Field Application
To evaluate the change law of the ground pressure with the thickness of bedrock, the working status information of the hydraulic supports was recorded. Measurement stations were arranged within 200 m from the stopping line in panel 3301. The working status information of the hydraulic supports is shown in Table 2. For roof with a thick bedrock, the working resistance of the hydraulic supports in the weighting period and the normal is 8891 KN and 6744 KN, respectively, with a dynamic load coefficient of 1.57. Compared to the roof with thin bedrock, the working resistance of the hydraulic support is significantly reduced, and the dynamic load characteristic is weakened. This is because the thick bedrock has higher load-bearing capacity, forming a double key layer structure. Therefore, the upper key layer can form triple articulated arch and obstruct the alluvium loads transferring to the hydraulic supports. As a result, the hydraulic supports only need to bear the weight of the rock strata below the upper key layer and further reduce loads to the hydraulic supports. Based on the above analysis, it can be obtained that the hydraulic supports have the characteristics of larger loads and obvious dynamic loads when the bedrock thickness is small. However, when the bedrock thickness is larger, the action of dynamic loads is significantly reduced, which confirms the correctness of the single-double key bed theory of thin bedrock stope.
In actual engineering practice, many measures including improving initial support force, ensuring sufficient retractability of supports during weighting period, and removing the crushing supports timely were applied during panel retreating. In total, 96% of hydraulic supports are in state of tilt-up working condition, end distance is less than 0.8 m, and the rate of roof falling height greater than 0.5 m has been reduced from 20% to 5%. Meanwhile, rib sliding accidents decreased significantly, and the depth is less than 0.3 m. The photograph in Figure 16 shows the support-rock mass control effective in panel 3301. Engineering practice has proven that the above measures can effectively reduce accidents in the large mining height stope covered by the TBLA.
5.3.2. Implications
It is extremely significant to clarify the law of stope ground pressure for coal mining. However, according to the literature review in the introduction, it can be found that the approaches for the large mining height stope covered by the TBLA have some limitations. Therefore, it is important to propose an integrated method of all-around field monitoring, numerical simulation, and theoretical calculation.
The findings in this paper are beneficial for a better understanding of the ground response of the large mining height stope covered by the TBLA. It should be noted that the study is based on a specific geological condition. And more case studies are needed to prove the accuracy of the law of ground response of the large mining height stope covered by the TBLA. However, the modelling procedure and evaluation method presented in this paper can be used for other similar projects.
6. Conclusion
(1)During panel 3301 retreating period, strong roof periodic weighting phenomenon occurred regularly, and the weighting step is small. Due to high with long time of mine pressure, dynamic characteristics of hydraulic supports are prominent and characterized by high safety valve opening rate and poor working conditions including tilt down operation and excessive end distance.(2)The numerical results show that separation deformation is significant; that is, the subsidence of the immediate roof is significantly greater than the subsidence of the main roof. When the separation deformation is large, and the working resistance of the support is insufficient, sliding instability of the main roof will generate more loads to the immediate roof, further resulting in support crushing.(3)Multiple measures, including improving sufficient initial support force and the worker operation level and ensuring sufficient retractability of supports, successfully controlled the stability of support-surrounding rock in the large mining height stope covered by TBLA. The ground pressure law of stope and measures presented in this paper can provide references for other similar geological conditions.Data Availability
The data used to support the findings of this study are included in the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this study.
Authors’ Contributions
Xiaoping Li conceptualized the study; Chao Wang and Huaixuan Cao were responsible for data curation; Chao Wang and Xianyang Yan performed formal analysis; Xipo Zhao and Guanglei Zhou investigated the data; Shibao Shen was responsible for software; Xiaoping Li and Guangchao Zhang were responsible for methodology; Guangzhe Tao was responsible for software; Guangchao Zhang performed study supervision; Xiaoping Li prepared the original draft; Guangchao Zhang and Guangzhe Tao reviewed and edited the manuscript. Guanglei Zhou has made an important contribution to the field investigation and data analysis. However, his contributions are mixed with others due to the longtime interval. Therefore, the authors sincerely request to add his name in terms of their actual contribution to this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (no. 51904164) and the Natural Science Foundation of Shandong Province (ZR2018QEE001).