Abstract

In underground metal mines, goaf brings huge safety risks. It is an important part of mine to evaluate goaf stability and determine the best goaf handling measures. However, the evaluation and decision are often separated; they need to be unified. Meanwhile, stability evaluation and decision for goaf are a complex system engineering of rock masses in underground mining, and subjectivity exists in the evaluation and decision process. Under these conditions, it is necessary to minimize subjectivity, and the results of stability evaluation also need to be considered comprehensively in determining handling measures. In this paper, the fuzzy theory was adopted based on the fact that the goaf stability evaluation and handling measure decision were fuzziness. Firstly, the goaf stability model was established by a two-layer fuzzy comprehensive evaluation. It took into account 12-factor indexes of goaf with engineering empirical approaches and divided them into 3 categories according to their engineering categories. The model improved the applicability of the goaf stability evaluation results, and the results were the basis for goaf handling measures as well. Secondly, a decision model of the goaf handling measures was established by multiobjective fuzzy optimization. It consisted of five goaf handling measures and five evaluation indexes. The model provided a comprehensive decision and optimal scheme for goaf handling. Two models were also applied to the Paishanlou gold mine and achieved a good handling result. The practical application showed that the two models were feasible.

1. Introduction

In underground metal mines, a goaf is formed. When the ore is mined and the mined-out area is not filled, especially the room and pillar mining, it forms many goafs of different sizes, shapes, and buried depths. In 2015, according to the survey results of 457 large and medium-sized mines in 25 provinces and cities conducted by the State Administration of Work Safety of China, there are 432 million m3 goaf, 80% of the goaf is below 10000 m3, but the remaining 20% of the goaf is above 10000 m3 and the total volume accounts for more than 50%. In particular, the number of goafs above 30000 m3 is less than 6%, but the volume accounts for 30% of the total [1].

The harm of the underground goaf is significant to the mine. Firstly, during the long-term creep and groundwater action, the strengths of surrounding rock and pillar decrease, the pillars occur continuously unstable or the roof suddenly caving, and at the same time, shock gas and waves are caused. These bring huge safety risks to personnel, equipment, facilities, buildings, etc. On November 6, 2005, a large-scale collapse accident occurred in Xingtai gypsum mines in China, which caused a large number of casualties including 37 deaths and 38 injuries, and a total of 88 living rooms on the surface were destroyed. The direct property loss was up to 7.74 million RMB [2]. Secondly, during the process of mining, the surrounding rock in the goaf is affected by blasting vibration, which leads to the development of fractures in the rock mass. This may lead to the formation of fracture connection networks and water inrush accidents. Accidents may submerge the tunnel and cause losses. On the afternoon of March 28, 2009, a water burst accident occurred in the old goaf of the Xishimen iron mine, which caused 8 people to die. Therefore, timely and appropriate goaf handling is the most fundamental measure to eliminate the safety threat. And the stability evaluation of the goaf is the basis for determining the goaf handling measures.

The stability of goaf is restricted by many factors, such as rock strength, joint conditions, groundwater, goaf parameters, and engineering factors. These factors are interrelated with each other. Therefore, the stability evaluation of the goaf is an extremely complex system engineering. The evaluation of goaf stability is a comprehensive evaluation for the goaf risk degree. The evaluation results have important guiding significance for determining goaf handling measures and timing. Swift and Reddish [3] considered the probabilistic approach to analyze stability problems. Yavuz [4] proposed a method for estimating the distance to return of the cover pressure and the stress distribution in the goaf. Hu and Li [5] presented a Bayes discriminant analysis method to identify the risk of complicated goafs in mines. Wang et al. [6] analyzed the creep failure of a roof stratum seated on pillars in the mined-out area through a newly developed visco-elastic model. Zhou et al. [7] proposed an approach to forecasting large-scale goaf instability that combined particle swarm optimization and support a vector machine. Sun et al. [8] proposed a three-zone model to analyze and evaluate the stability of the goaf. Hu et al. [9] built the RS-TOPSIS model to predict the hazard degree of goafs based on the results of the expert investigation. Zhang et al. [10] developed a susceptibility assessment system to define the risk from mine collapse for coalfields across the mining area based on the principles of fuzzy mathematics and the analytical hierarchy process. Li et al. [11] proposed a method to calculate releasable space in strata based on the characteristics of pore distribution in the rock strata above the goaf. Xiao et al. [12] established the hazard evaluation model of goaf by using information entropy and unascertained measurement theory. Guo et al. [13] established a risk assessment model with seven main assessment factors for expressway construction site instability based on fuzzy theory. In these studies, the influencing factors are considered in detail. However, goaf stability is affected by many uncertain factors and the influence degree of each factor is also different. These factors and the evaluation results have certain fuzziness and subjectivity. Therefore, it is a good way to combine the fuzzy theory with the classical comprehensive evaluation method. This will make the evaluation results as objective as possible and obtain more appropriate evaluation results.

For the goaf left by underground mining, there are usually four goaf handling measures: closing, caving, reinforcing, and filling. However, the specific conditions of each mine are different, the positions and shape characteristics of each goaf are also different, and the goaf handling measures are often different. Therefore, it is an important part of the mine to determine the best goaf handling measures in terms of technology, economy, safety, and reliability according to various factors affecting the decision-making of handling measures.

In this paper, based on the goaf stability evaluation and decision-making of goaf handling measures that are fuzzy, the fuzzy theory is adopted to establish a stability evaluation model and a decision model for goaf in underground metal mines. In the establishment process, evaluation and decision are considered comprehensively. These models are applied to the Paishanlou gold mine.

2. Establishment of Goaf Stability Model

The two-layer fuzzy comprehensive evaluation mathematical model (in Figure 1) is adopted for the goaf stability evaluation, and the processes are as follows.

2.1. Determination of Factor

There were many factors affecting goaf stability, which could be summarized as the factors of rock mass quality, goaf parameters, and induced factors. The factors of rock mass quality could refer to rock mass classification standards, such as Q classification [14, 15], RMR classification [16, 17], MRMR classification [18, 19], and engineering rock mass classification standard [20]. These factors included rock uniaxial compressive strength, rock point load strength, rock quality index RQD, joint spacing, intactness index of the rock mass, volumetric joint count of the rock mass, and joint occurrence. The uniaxial compressive strength and the point load strength both represented the rock strength, and the point load test was a more practical, time-saving, and economical method compared to the uniaxial compressive strength [21]. Many studies [2225] also have shown that the point load strength has a good correlation with the uniaxial compressive strength. Therefore, the point load strength was chosen to represent the strength rock. The RQD, joint spacing, intactness index of the rock mass, and volumetric joint count of the rock mass [26] all represented the characteristics of the joint, and these factors were related or the same parameter was reflected in different aspects. Considering the convenience of field measurement, the intactness index of the rock mass was determined to represent the joint characteristics. Another factor was the joint occurrence. The study [20] has shown that the angle between the joint occurrence and the goaf direction affects the stability of the rock mass. Because the joint occurrence was complex in the actual rock mass, and it was difficult to study the influence of each joint occurrence, the dominant joint occurrence was selected.

The goaf parameters mainly included the goaf span, goaf area, goaf volume, goaf buried depth, goaf height, height span ratio, goaf shape, goaf engineering layout, and mine pillar and layout. Among them, the hydraulic radius could be used to express the goaf span, goaf area, goaf volume, goaf height, height span ratio, and goaf shape. The hydraulic radius refers to the ratio of area to the perimeter of the goaf.

The induced factors of the goaf could be summarized as groundwater, protective measures, mining disturbance, and influence of adjacent goaf.

Based on the above analysis, the actual application, and site conditions of the mine, the factor set U1 was established; U1 = {U11, U12, U13, U14, U15, U16, U17, U18, U19, U110, U111, U112} = {point load strength, intactness index of the rock mass, dominant joint occurrence, hydraulic radius, goaf buried depth, goaf height, engineering layout, mine pillar and layout, groundwater, protective measures, mining disturbance, adjacent goaf}. These factors are shown in Figure 1.

2.2. Factor Layer and Assessment Ranks

Due to the fact that many factors affected goaf stability, if a single-layer fuzzy comprehensive evaluation was adopted, the factor weight might not be determined, or the weight value was too small, which led to the distortion of the evaluation results. For this reason, 12 factors in the first layer were divided into 3 categories according to their engineering categories, which were integrated into 3 factors in the second layer; that is, U2 = {U21, U22, U23} = {rock mass quality factors, goaf parameters, induced factors}. The specific factor layering and meaning are shown in Figure 1.

The assessment ranks were a set reflecting goaf stability. Generally, a 4-rank classification was used; that is, V = {V1, V2, V3, V4} = {stable I, basically stable II, understable III, unstable IV}. Among them, the stable I indicated that the goaf did not need to be handled and monitored. The basically stable II indicated that the goaf needed to be handled and monitored to ensure safe production within the affected area. Understable III indicated that the goaf needed to be handled and monitored, and emergency plans needed to be formulated. The unstable IV indicated that handling measures needed to be taken immediately and strengthen monitoring, and personnel and equipment within its affected area must be evacuated immediately.

Based on the factor layer and assessment ranks, each factor index was divided into 4 levels according to the stability classification of the engineering rock mass. In addition, the results are shown in Table 1; there are both qualitative and quantitative factor indexes. The membership degree of each qualitative factor index could be obtained by counting the frequency by several measuring personnel, such as U13 dominant joint occurrence, U17 engineering layout, U18 mine pillar and layout, U110 protective measures, U111 mining disturbance, and U112 adjacent goaf. The membership degree of each quantitative factor index could be determined by measuring the value and membership function, such as the U11 point load strength, U12 intactness index of the rock mass, U14 hydraulic radius, U15 goaf buried depth, U16 goaf height, and U19 groundwater.

2.3. Membership Function

In order to determine the membership degree of each quantitative factor index and make the model more effective, the original measured value was normalized, and the input data were in the [0, 1] interval. The normalized formula iswhere f(Uij) is the normalized value; Uij is the original measured value; pimax and pimin are the maximum values and minimum of the original measured value, respectively; qimax and qimin are the maximum values and minimum of the quantitative range corresponding to the original measured value, respectively; Formula (1) shows that the larger Uij is, the more stable of the goaf. Formula (2) is that the larger Uij is, the more unstable it is.

After normalizing the original measured value, it was necessary to establish the membership function of each quantitative factor index, and the membership degree could be calculated. According to the conclusion that different membership functions were equivalent by Su et al. [27], this paper used the fuzzy reasoning method to establish membership functions based on the characteristics of fuzzy sets.

In rock engineering, the membership function usually adopted the intermediate type. That was, the membership degree was 0.5 at the endpoint of the interval range and the state was the fuzziest. The membership degree was 1 at the middle point and the neighborhood of the interval range and the state was the clearest. The membership function formula iswhere δ is the neighborhood value centered on the middle point, and the membership degrees are 1 in the neighborhood range.

2.4. Fuzzy Judgment of First Layer Factor

According to formulas (1) to (3), the membership degree of each quantitative factor index was calculated. Combining the membership degrees of the qualitative factor index, a fuzzy comprehensive evaluation matrix could be established. Based on the three categories, three fuzzy comprehensive evaluation matrixes could be established: where Ak(f(Uij)) is the membership degree; R21 is the evaluation matrix of the rock mass quality; R22 is the evaluation matrix of the goaf parameter; R23 is the evaluation matrix of the induced factor.

After establishing the fuzzy comprehensive evaluation matrixes, it needed to determine the weight of the first layer factor. In the evaluation of goaf stability, each factor index had a different influence on goaf stability, and it was necessary to determine the weight of each factor index. Weight was an important part of the comprehensive evaluation. Among the methods to determine the weight, the analytic hierarchy process (AHP) was an effective method to determine the weight coefficient.

The AHP [28, 29] divided a complex system into several levels and factors. By comparing the importance of the two factors, the weight coefficient was determined. It was an effective multiobjective planning method in system engineering. The essence of AHP was a kind of decision-making thinking mode, which combined qualitative analysis and quantitative analysis in the decision-making process. The steps of AHP were as follows: determination of judgment object and factors, establishment of a judgment matrix, calculation of the order of the relative importance, and consistency check.

The object was the stability of the goaf, and the factors had three sets included U21 = {U11, U12, U13}, U22 = {U14, U15, U16, U17, U18}, and U23 = {U19, U110, U111, U112}. In the judgment matrix, the value reflected the relative importance. Generally, the judgment matrix was obtained by the scale method of 1∼9 and its reciprocal. The judgment matrix can be established as follows: where P21 is the judgment matrix of rock mass quality, P22 is the judgment matrix of goaf parameters, and P23 is the judgment matrix of induced factors.

In the order of the relative importance, the maximum eigenvalue of P21 with λ21max = 3.04, and the eigenvector X21 is as follows:

The fuzzy relative weight coefficient was obtained by normalizing the eigenvector X21. The weight vector C21 is

The maximum eigenvalue of P22 with λ22max = 5.13; the eigenvector X22 is as follows:

The weight vector C22 is

The maximum eigenvalue of P23 with λ23max = 4.01; the eigenvector X23 is as follows:

The weight vector C23 is

In the consistency check, it was necessary to check whether the weight was reasonable. The formula of the consistency index (CI) is as follows:where n is the number of factors.

The consistency ratio CR can be obtained as follows:where RI is the average random consistency index. n is 3, 4, and 5, and RI is 0.52, 0.89, and 1.12, respectively.

When the CR was less than 0.10, it showed that the judgment matrix satisfied the consistency test; that is, the distribution of the weight coefficient was reasonable. The results were CI21 = 0.02, CR21 = 0.04, CI22 = 0.033, CR22 = 0.029, CI23 = 0.003, and CR23 = 0.004. These results showed that the distribution of the weight coefficient was reasonable.

After the weight vector C2i and the fuzzy comprehensive evaluation matrix R2i were determined, fuzzy subset B2i could be obtained through a fuzzy linear change. where “◦” is a synthetic operator for weight vector C2i and fuzzy comprehensive evaluation matrix R2i.

2.5. Fuzzy Judgment of Second Layer Factor

Based on the fuzzy comprehensive evaluation matrix of the first layer, the fuzzy comprehensive evaluation matrix R of the second layer can be established:

The next was to determine the weight of the second layer factor. AHP was still used to determine the weight of the second layer factor. The judgment object was goaf stability. The set of judgment factors was U2 = {U21, U22, U23}. When determining the relative importance value in judgment matrix P, the main basis was the current research results or scoring value in the common evaluation methods. The judgment matrix P could be obtained.

In the order of the relative importance, the maximum eigenvalue of P with λmax = 3.02; the eigenvector X is as follows:

The fuzzy relative weight coefficient was obtained by normalizing the eigenvector X. The weight vector is

In the consistency check, the consistency index CI was 0.01, and the consistency ratio CR was 0.02 and less than 0.10, which showed that the distribution of the weight coefficient was reasonable.

After the weight vector C and the fuzzy comprehensive evaluation matrix R were determined, fuzzy subset B could be obtained through a fuzzy linear change. That is,

According to the principle of maximum membership, the maximum value bmax = max{bi} was obtained, and the value of i was 1 to 4. The rank of bmax = max{bi} was the rank of comprehensive evaluation.

2.6. Application in Paishanlou Gold Mine

The Paishanlou gold mine belongs to the China National Gold Group, which is the key enterprise in Northeast China. The Paishanlou gold deposit is a large-scale hydrothermal gold deposit. The surrounding rock is mylonite, and the remote rock is primary mylonite and mylonitic rock. The dip angle of the ore body is 35°∼55°. The ore body is located between the elevations of 465∼−30 m. After open pit mining to +300 m, the mine transferred to underground mining. Due to the slow dip angle of the ore body and the increase the production capacity, the open stope method is adopted in the two middle sections of 300 m and 225 m, forming large-scale discontinuous goafs (in Figure 2).

The horizontal area of the goafs has reached 11250 m2 above the 275 m level, and the volume is approximately 500000 m3. The roof of the goaf obliquely intersects with the surface. In the section of exploration line X (in Figure 2(b)), the dip angle of the goaf roof is 25°–35°, the surface boundary is the open-pit mining boundary, and the boundary slope angle is 29°–50°. The thinnest thickness of the goaf roof is only 36 m. The original open pit is used to store dry toxic tailings. The design storage elevation is +420 m. The locations of the open pit and goaf are shown in Figure 2(a).

Therefore, it was urgent to evaluate goaf stability in the Paishanlou gold mine. According to the results of the stability evaluation and the technical conditions of the mine, the optimal goaf handling measures were determined. The ultimate goals were to achieve safe mining of ore bodies and safe management of goafs. Then, stability evaluation of the upper goaf and lower goaf was carried out by the goaf stability model. The indexes and results are shown in Table 2.

As shown in Table 2, the stability rank of the upper goaf is unstable IV based on the principle of maximum membership. However, the membership in unstable IV is close to membership in understable III, and the stability ranks of the upper goaf are unstable IV and understable III. The stability rank of the lower goaf is basically II. The results showed that the goafs need to be handled to ensure safe mining, and these results were the basis for goaf handling measures as well.

3. Establishment of the Goaf Decision Model

In the process of the goaf decision model, many uncertain factors were involved, such as the geological conditions of the mine, goaf stability, goaf location, goaf shape, and filling of the goaf. Different goaf handling measures had different handling costs, handling effects, and time required for handling. They also depend on the specific production situations of the mine. The feasibility and other indexes were also different for the goaf with different stability ranks. Therefore, the advantages and disadvantages of goaf handling measures were fuzzy. And multiobjective fuzzy optimization theory could make a good comprehensive evaluation of the goaf handling effect and obtain the optimal handling scheme.

3.1. Multiobjective Fuzzy Optimization

For the goaf handling measures, there were usually four categories: closing, caving, reinforcing, and filling. Each category also included a series of subhandling measures, such as filling could be divided into waste rock dry filling, tailings filling, and cemented filling, and caving could be divided into induced natural caving and blasting caving. Considering the mine situation, the goaf handling measures were preliminarily determined, which were Ds = {ds1, ds2, ds3, ds4, ds5 } = {blasting and caving, induced natural caving, waste rock dry filling, tailings filling, cemented filling}.

The comprehensive indexes of goaf handling mainly included necessity, feasibility, requirements of laws and regulations, treatment effect, handling cost and work efficiency, and handling time. The necessity of handling was determined according to the goaf stability evaluation results, the results in Table 2, and mining status showed that goaf handling was necessary. The feasibility mainly considered safety (based on the goaf stability evaluation) and construction conditions. Therefore, the evaluation index set Ps was determined, Ps = {ps1, ps2, ps3, ps4, ps5} = {feasibility, requirements of laws and regulations, handling effect, handling cost, handling time}. According to the five goaf handling measures and five evaluation indexes, the index characteristic matrix is expressed as where xij is the characteristic value of index i and measure j.

In the process of fuzzy optimization or decision-making, in order to solve the possible objectivity of membership and the subjective in the process of determination, the relative membership matrix was established by using the index relative membership. For larger and better indexes (handling effect, feasibility, requirements of laws and regulations), the relative membership formula is

For smaller and better indexes (necessity of handling, handling cost, handling time), the relative membership formula iswhere rij is the relative membership degree of index i and measure j; ∨ and ∧ are the characteristics of taking large and taking small, respectively.

The index characteristic matrix X is transformed into the index relative membership matrix by formulas (21) and (22); that is,

The relative membership degree of the superior measure is

The relative membership degree of the inferior measure is

In the goaf handling system, the weight vector Ws is

The weight vector Ws was still determined by AHP, and the judgment matrix is

Calculating the maximum eigenvalue λmax = 5.11 of the judgment matrix Ps, the consistency index CI was 0.028, and the consistency ratio CR was 0.025 and less than 0.10. That is, the weight vector Ws is

In fuzzy theory, the weighted optimal distance Djg is

The weighted inferior distance Djb iswhere p is the distance parameter and p = 1 is the Hamming distance.

The optimization criterion was that the sum of the weighted optimal distance square and the weighted inferior distance square was the smallest. That is,

Finding the derivative of formula (31), and making it equal to zero, the solution is

Formula (32) was a fuzzy optimization theory model of the goaf handling measures. In the set of measures that meet the index constraints, the measure with the highest relative membership degree uj was the satisfactory measure, and the sequence of uj from the largest to the smallest was the satisfactory sequence.

3.2. Application in Paishanlou Gold Mine

To determine the goaf handling measures, a group of nine experts scored and calculated the feasibility, requirements of laws and regulations, handling effect, handling cost, and handling time for five goaf handling measures. The experts were composed of mining technicians, mine leaders, material planners, and university researchers. Then, the highest and the lowest scores were removed, respectively, and the mean value of the remaining results was taken. Among them, the feasibility, treatment effect, and requirements of laws and regulations were all scored on the 10-point system. The final score and calculation results are shown in Table 3. Then, the relative membership degree uj is calculated by the multiobjective fuzzy optimization model (in Table 3).

The results of fuzzy optimization for the goaf decision are shown in Table 3, the handling measure of the upper goaf is waste rock dry filling, and the handling measure of the lower goaf is induced natural caving. However, considering the relevance of goafs (in Figure 2) and the urgency of handling, a comprehensive handling measure is proposed as shown in Figure 3. Firstly, the filling shaft is drilled in the open pit slope to fill the upper goaf with waste rock (in Figure 3(a)). The mining truck is used to fill the upper goaf quickly (in Figure 4(a)), and the total amount of waste rock filling is 372000 m3. After the upper goaf is filled, the cemented filling is carried out at the top of the upper goaf, to better protect the stability of the slope.

Secondly, the tunnel is excavated in the pillar of the upper and lower goaf, and blasting medium-length holes are drilled so that the upper and lower goaf is connected by blasting and this engineering induces natural caving of the rock mass above the goaf (in Figure 3(b)). Two days after blasting, a surface subsidence pit is formed, and it has a depth of 10 to 43 m and a total area of 10451 m2 (as shown in Figure 5(a)).

Finally, the surface subsidence pit is filled with waste rock and covered with soil (in Figure 5). When the subsidence pit is stable, the waste rock in the waste rock field is transported by mining trucks, and the subsidence pit is filled from north to south (in Figure 3(c)). By September 2012, the landform was basically restored, and the total amount of waste rock filled is 333566 m3. After the completion of waste rock filling, the engineering of covering with new soil is carried out in layers (in Figure 5). The first layer is covered with 0.5 m loess (in Figure 4(b)), after leveling, the protection against osmosis fabric is covered, and then the new soil is covered with a 0.5 m loess layer to protect against osmosis fabric. In the spring of 2013, soybean, corn, and other crops were planted on the new soil, to prevent soil erosion and achieve economic benefits.

4. Conclusions

To comprehensively solve the safety risks from goafs in underground metal mines, minimize subjectivity in the evaluation and decision process, and comprehensively consider different handling measures, an evaluation model of goaf stability and a decision model of handling measures were provided based on fuzzy theory.

The goaf stability model was established by the two-layer fuzzy comprehensive evaluation. It took into account 12-factor indexes of goafs with engineering empirical approaches and divided them into 3 categories according to their engineering categories. The membership degree of indexes was determined based on the normalized formula and membership function, and weights were determined by the analytic hierarchy process. This model improved the applicability of the goaf stability evaluation results, and the results were the basis for goaf handling measures. The model was applied to the Paishanlou gold mine. The results showed that the stability ranks of the upper goaf were unstable IV and understable III, and the stability rank of the lower goaf was basically II. In addition, the goaf needed to be handled.

Based on the results of the goaf stability evaluation, the decision model of the goaf handling measures was established by multiobjective fuzzy optimization. It consisted of five goaf handling measures and five evaluation indexes. The model provided a good comprehensive decision and optimal scheme for goaf handling. This model was also applied to the Paishanlou gold mine. The results showed that the handling measure of the upper goaf was waste rock dry filling, and the handling measure of the lower goaf was induced natural caving. These results were integrated with the actual situation of the mine, and a comprehensive handling measure was proposed. That was, the first step was to fill the upper goaf with waste rock, and cemented filling was also carried out at the top of the upper goaf. The second step was to connect the upper and lower goaf and induce natural caving of the rock mass above the goaf. The third step was to fill the surface subsidence pit with waste rock and cover it with soil.

By 2020, the change in the covered soil in the surface subsidence pit was tracked and observed, as shown in Figure 6. There was small subsidence and some cracks only at the corners, and the integrity of the covered soil was not destroyed. This showed that the comprehensive handling measure was successful. It also showed that the results of goaf stability evaluation and handling measure decisions were correct, and the stability model and decision model for goafs were feasible.

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.

Acknowledgments

This study was jointly supported by grants from the National Key Research and Development Program (Grant no. 2020YFC1909801), Fundamental Research Funds for the Central Universities (Grant no. N2101024), Scientific research project of Inner Mongolia University of Technology (Grant no. ZZ202114), and Natural Science Foundation of Liaoning Province (Grant no. 2021-MS-096). The authors are grateful for the support.