The Bainiuchang Mine is located 42.5 km at a direction of 75° from the downtown area of Mengzi City, Yunnan Province, China (coordinates: 103°45′ to 103°48′, 23°26′30″ to 23°29′30″). The study area (Figure 2), the Baiyang mining area, is located in the northwestern part of the mining area. Mining activities began in 1994, with the peak mining period occurring between 2000 and 2010.

Numerous factors contribute to the instability in mined-out areas. Through field surveys and literature review, 13 indicators, including rock mass structure, geological structure, rock quality, ore body dip angle, goaf volume, maximum exposed area, pillar area ratio, ratio of rise to span, embedding depth, exposure time of goaf, mined-out area water, situation of adjacent mined-out areas, and goaf management (Chen et al., 2024; He et al., 2022; Qin et al., 2019; Ren et al., 2022). An evaluation index system for the instability risk in mined-out areas was constructed. Based on the relevant literature and field surveys, the risk levels were identified as I, II, III, and IV, with the risk level classification detailed in Supplementary Table S1 (Guo et al., 2022; Jia et al., 2016; Lang and Chen, 2015). The corresponding standards for each evaluation indicator are listed in Supplementary Table S2 (Hu et al., 2017; Hu, 2022; Ma, 2015; Ren Honggang, 2020; Wang, 2020).

The investigation of the mined-out areas commenced in April 2021 and concluded in October 2022. The team conducted handheld 3D laser scanning surveys of accessible and safe mined-out areas (field surveys and scanning results are shown in Supplementary Figures S1A, C). For areas with high safety risks where direct human entry was not feasible, UAV photography equipment was used for surveying (field survey and scanning results are shown in Supplementary Figures S1B, D). A total of 195 mined-out areas (or groups of mined-out areas) were identified in the Baiyang mining area. Through field surveys and interviews with mining engineering technicians, data on 13 evaluation indicators were collected from 195 mined-out areas in the Baiyang area, as detailed in Table 1.

The integrity and distribution of data in Table 1 were assessed using violin plots. Owing to space constraints, violin plots for only four indicators (X3, X4, X5, and X6) are presented. Figure 3 shows that under natural conditions, the data patterns in Table 1 generally align with a normal distribution.

To reduce the subjectivity of the computational results and simplify the calculation process, this study employs the CRITIC method to comprehensively measure the contrast strength and conflict of the evaluation indicators. This method objectively quantifies the importance for ranking purposes. A preliminary ranking was established and then adjusted based on expert judgment and the initial ranking. Subsequently, a binary comparison method (Wang et al., 2014; Lu et al., 2014) was used to calculate the subjective weights of the indicators without consistency testing. This approach is relatively straightforward and can replace traditional methods such as the AHP to determine the subjective weights of indicators related to the instability risk of mined-out areas.

The Gini coefficient weighting method is an objective weighting method that calculates the Gini coefficient of the evaluation index and normalizes the Gini coefficient of each index to obtain the weight of the evaluation index. The calculation process is as follows (Dai et al., 2024; Zhao et al., 2024; Li, 2015; Li et al., 2014):

Calculate the Gini coefficient value of the evaluation index: G k = ∑ i = 1 n ∑ j = 1 n X k i − X k j / 2 n 2 μ k μ k ≠ 0 ∑ i = 1 n ∑ j = 1 n X k i − X k j / n 2 − n μ k = 0 (8) where G k is the Gini coefficient of the k-th index, n is the total number of index data points, and X k i is the i-th value of the k-th evaluation index. μ k is the expected value of the k-th indicator.

Evaluation index Gini coefficient weight calculation: The Gini coefficient weights of the different indices were obtained by normalizing the Gini coefficient values of each evaluation index. g k = G k / ∑ i = 1 m G i (9)

The problem of determining the weight of the index combination is solved using game-theory aggregation. Let the weight vector of m indicators calculated by the importance-binary comparative analysis method be W → 1 T , and the weight vector of m indicators obtained by the Gini coefficient objective weighting method be W → 2 T . Let W → T be the linear combination of the weight vectors W → 1 T and W → 2 T . By calculating, the deviation between W → T and W → 1 T and W → 2 T is minimized, α p can be obtained, and then normalized. The combined weight vector is obtained as follows (Dai et al., 2024; Zhao et al., 2024; Li, 2015; Li et al., 2014): W → = ∑ P = 1 2 α P W → P T (10) min ∑ P = 1 2 α P W ⇀ P T − W → l T (11)

The variable weight theory, introduced by Wang (1985) in the 1980s, posits that the weight of an evaluation indicator dynamically adjusts based on changes in its corresponding state indicator value. Subsequently, Li (1995); Li (1996) provided axiomatic definitions and calculation formulas for three types of variable weights: penalty, incentive, and neutral.

Let the state vector of the evaluation index x j be X = x 1 , x 2 , ⋯ , x m , and let the constant weight vector W 0 = w 1 0 , w 2 0 , ⋯ , w m 0 be obtained using the above combination weighting method. Variable Weight Vector W x = w 1 x , w 2 x , ⋯ , w m x . The variable weight calculation formula (Equation 12), S X = S 1 X , S 2 X , ⋯ , S m X is the state-variable weight function. W X = W 0 • S X ∑ i = 1 m w i 0 • S i X = w 1 0 • S 1 X ∑ i = 1 m w i 0 • S i X , w 2 0 • S 2 X ∑ i = 1 m w i 0 • S i X , ⋯ , w m 0 • S m X ∑ i = 1 m w i 0 • S i X (12)

Duan (2003) proposed a penalty-incentive segmented coordination theorem for state variable weight functions, dividing the penalty phase into initial penalty, strong penalty, and veto stages. Li and Wu (2015); Wu and Li (2016) constructed a four-segment (Equation 13) local variable weight state function based on the variable weight theorems proposed by Li Hongxing et al., further dividing the incentive interval into initial and strong incentive segments. S j x = e a 1 d j 1 − x + c − 1 c e a 2 x − d j 2 + c − 1 e a 3 x − d j 3 + e a 3 d j 3 − d j 2 + c − 2   0 , d j 1 ) [ d j 1 , d j 2 ) [ d j 2 , d j 3 ) [ d j 3 , 1 ) (13)

Among them, c, a ₁, a ₂, and a ₃ are the weighting adjustment parameters, and d _j1, d _j2, and d _j3 are the threshold values for the variable weight intervals of the j-th factor. The variable weight vector is calculated according to Equations 12, 13. Currently, there is no unified method for determining variable weight intervals and variable weight parameters. Li and Wu (2015); Wu and Li (2016) proposed using the K-means method to classify the normalized data of the primary control factors. Then, the threshold values for the variable weight intervals of each factor are determined based on the classification-critical values of various indicators and the calculation formulas. By sequentially constructing evaluation units that meet the constraints and linking them with the constant weight weights to form equations, four weighting adjustment parameters are solved.

However, this method has the following shortcomings (Li and Wu, 2015; Li et al., 2023; Wu and Li, 2016):

1. Clustering analysis results can only indicate data similarity and aggregation degree but cannot reflect the magnitude relationship of the factor state values. The basis for determining interval thresholds through clustering analysis is insufficient.

Based on the characteristics of the variable weight model for evaluating the instability risk of mined-out areas, improvements were made to address its shortcomings.

1. For the local state function, the penalty and incentive intervals in Equation 13 were adjusted. The state function is shown in Equation 14 and illustrated in Supplementary Figure S3.

According to Equation 1, the indicators are processed as dimensionless. Using Equations 2, 3, the standard deviation of the evaluation indicators for the goaf instability risk is calculated. Using Equations 4, 5, the correlation coefficients of the evaluation indicators are calculated, and a heat map (Figure 5) is drawn. Using Equation 6, the importance scores of the evaluation indicators are calculated and ranked preliminarily. The results are summarized in Table 3.

After consulting relevant experts, the preliminary ranking results of the 13 evaluation indicators were adjusted. The adjusted ranking is: X1 ≻ X12 ≻ X13 ≻ X2 ≻ X7 ≻ X8 ≻ X5 = X6 ≻ X3 ≻ X11 ≻ X9 ≻ X4 ≻ X10. The membership degrees of the factors relative to each other were determined based on the ranking results, as shown in Table 2 and Supplementary Figure S2. The weights of the evaluation indicators were normalized, and the calculation results are listed in Table 4.

The Gini coefficients of the 13 evaluation indices were calculated using Equation 8 and normalized using Equation 9 to obtain the objective weights of the indices: The calculation results are listed in Table 5.

According to Equations 10, 11, the combination coefficients of the subjective and objective weights were determined as α 1 = 0.33 , α 2 = 0.79 . After normalization, the weight coefficients were reallocated, resulting in α 1 * = 0.29 , α 2 * = 0.71 . The combined weights of the 13 evaluation indices that influence the instability risk of the mined-out areas were calculated using Equation 10: The results are as follows: W = (0.0666, 0.0920, 0.0286, 0.0233, 0.0919, 0.1099, 0.0977, 0.0731, 0.1315, 0.1064, 0.1078, 0.0345, 0.0367).

Variable weight intervals were determined using the cumulative frequency of the normalized values of the influencing factors. Owing to the absence of predefined thresholds for categorizing data frequencies as “extreme” or “unique,” this study temporarily defined the cumulative frequency thresholds for state values at 15%, 33%, and 67%. Specifically, the normalized value corresponding to a cumulative frequency of 15% was defined as a strong penalty interval, 33% as an initial penalty interval, and 67% as a neutral interval. The remaining interval was considered the incentive interval. The threshold values for the variable weight intervals were determined by plotting the cumulative frequency intervals of the state values of the main controlling factors (Supplementary Figure S5). Thresholds for factors heavily influenced by human intervention can be manually defined. The final interval thresholds are listed in Table 6.

Based on Equations 14–17, the variable weight intervals of the 13 evaluation indices listed in Table 6 were multiplied by constant weights to obtain the comprehensive variable weight interval matrix. The final divisions of the variable weight intervals are listed in Table 7.

Assuming that the accuracy rate of the instability risk evaluation results for the goaf has a function f (a1, a2, a3, c) related to the adjustment parameters, the Monte Carlo method is used to determine the combination of a1, a2, a3, and c within their ranges to maximize f (a1, a2, a3, c). The specific steps are illustrated in Figure 6.

With the aid of MATLAB software, 1,000 assignment tests were conducted to obtain the accuracy rate f (a1, a2, a3, c) of the evaluation results as a function of the adjustment parameters a1, a2, a3, c. The resulting variation curve is shown in Figure 7.

a3 adjusts the amplitude of the weight changes in the incentive interval, a1 and a2 adjust the amplitude of the weight changes in the punishment interval, and c controls the weight adjustments of the entire variable weight model. From Figure 7, it can be concluded that the adjustment parameters a1, a2, and c are negatively correlated with the accuracy of the evaluation results, whereas a3 is positively correlated with the evaluation results. In addition, a1 and a2 have weaker adjustment strengths, whereas a3 and ccc have stronger adjustment strengths. The final adjustment parameters of the variable weight model are listed in Table 8.

Based on the constructed state variable weight vector formula and the determined variable weight interval and adjustment parameter values, the state variable weight vector Equation 19 for the risk assessment of goaf instability in the Baiyang mine area was established according to Equation 14. S j x = e 0.29 * 0.25 − x + e 0.66 * 0.42 − 0.25 + 0.01 − 2 [ 0 , 0.25 ) e 0.66 * 0.42 − x + 0.01 − 1 [ 0.25 , 0.42 ) 0.01 [ 0.42 , 0.64 ) e 0.67 * x − 0.64 + 0.01 − 1 [ 0.64 , 1 ) (19)

According to Equation 12, the variable weight value of each evaluation index is calculated, and the relationship between the constant weight and variable weight was drawn (Figure 8). The variable weight data of the 13 evaluation indices were statistically analyzed using Origin software. The results are shown in Table 9.

From the relationship diagram of the constant and variable weight values of the evaluation indicators (Figure 8), it is evident that the constructed variable weight model effectively adjusts the weights of the evaluation indicators during the process of assessing the goaf instability risk in the Baiyang mining area. Specifically, the variable weight model includes the weights of six indicators: rock mass structure (X1), geological structure (X2), rock quality (X3), ore body dip angle (X4), condition of adjacent goafs (X12), and management measures (X13), which significantly increase their weight values. Conversely, the model has “penalized” the weights of seven indicators: goaf volume (X5), maximum exposed area (X6), pillar area ratio (X7), height-span ratio (X8), burial depth (X9), exposure time (X10), and water accumulation in the goaf (X11), reducing their weight values in some goafs.

The average weight of the “neutral” interval (neither penalized nor incentivized) is smaller than that of the other intervals. This is because, under the weight adjustments of the variable weight model, the weights of the factors in the “neutral” interval are compressed. As shown in Table 9, the difference between the average variable and constant weight was not significant. This indicates that the variable weight model not only meets the decision-makers’ preferences for the primary control factors under various combinations of states but also considers the relative importance of the indicator weights. Therefore, in terms of the weight adjustment effects, the constructed variable weight model effectively adjusted the weights of the 13 evaluation indicators according to the predetermined variable weight approach, thus satisfying the expected outcomes.