Fast evaluation method for phosphogypsum slope stability based on extreme learning machine

Abstract

The application discloses a kind of phosphogypsum slope stability rapid evaluation method based on extreme learning machine, belong to geotechnical engineering safety evaluation technical field.The method first constructs the t-Copula joint distribution model of phosphogypsum mechanical parameter, is updated by bayesian in situ adaptive calibration, and joint sampling generates synthetic training dataset;Then train extreme learning machine, and the output weight regularization increment migration correction is corrected with a small amount of field data;Finally, based on the joint distribution model, the parameters are transformed to the standard normal space, and the HL-RF iterative solution of the most likely instability point is driven by the analytical gradient of extreme learning machine, and the safety factor and other indicators are output simultaneously.The application solves the problem that phosphogypsum yard cannot quickly provide safety margin evaluation results considering joint degradation path when parameters are strongly correlated and data are scarce, and the whole-link calculation time is less than 20 milliseconds.

Description

Technical Field

[0001] This invention relates to the field of geotechnical engineering safety assessment technology, specifically to a rapid assessment method for the stability of phosphogypsum slopes based on limit learning machine. Background Technology

[0002] Phosphogypsum is a solid waste generated during the production of phosphoric acid. Approximately 4.5 to 5.5 tons of phosphogypsum are produced for every ton of phosphoric acid produced. Phosphogypsum is typically stored in dedicated stockpiles using a layered compaction method. As the stockpile volume continues to increase, the slope height of the stockpiles is constantly rising. Some stockpiles have developed steep slopes exceeding 60 meters in height, and their stability directly affects the safety of life and property of surrounding residents and the protection of the ecological environment.

[0003] Currently, the finite element method (FEM) or limit equilibrium method is commonly used for slope stability analysis of phosphogypsum stockpile sites. The FEM method determines the safety factor by progressively reducing soil strength parameters and solving the elastoplastic equations until the calculation fails to converge. This method can automatically search for the most dangerous slip surface without pre-setting failure mode assumptions and has been widely used in engineering practice. However, phosphogypsum is a low-strength loose material with a cohesion typically only 10 to 50 kPa. Approaching the limit state, the material enters a plastic flow stage, making it difficult for the iterative solution of the nonlinear equations to converge. A single calculation often takes several minutes or even hours, failing to meet the engineering requirements for real-time response to parameter changes. More importantly, the mechanical parameters of phosphogypsum exhibit strong correlations—cohesion is significantly negatively correlated with water content, and unit weight is positively correlated with water content. These parameters change simultaneously under external forces such as rainfall and surcharges, and their synergistic deterioration effect far outweighs the additive effect of independent parameter changes. Traditional deterministic analysis methods treat each parameter as an independent variable, which cannot accurately model the joint degradation behavior between parameters, nor can they answer the core engineering decision-making question of "how much more can the parameters deteriorate in which direction before instability", resulting in the inability to effectively transform safety assessment results into dynamic early warning thresholds. Summary of the Invention

[0004] Technical Objective: To address the shortcomings of existing technologies in phosphogypsum slope stability analysis, which cannot simultaneously consider computational efficiency, parameter correlation modeling, and joint quantitative assessment of safety margin, this invention discloses a rapid assessment method for phosphogypsum slope stability based on limit learning machine.

[0005] Technical solution: To achieve the above technical objectives, the present invention adopts the following technical solution: A rapid stability assessment method for phosphogypsum slopes based on extreme learning machine includes the following steps: Step S1: Construct a Copula joint distribution model of the mechanical parameters of phosphogypsum: Obtain statistical data on the cohesion, internal friction angle, bulk density and moisture content of phosphogypsum, fit the marginal probability distribution function of each parameter respectively, establish a joint distribution model of the four mechanical parameters using the t-Copula function, and estimate the correlation parameters of the joint distribution model using the maximum likelihood method.

[0006] There is a significant physical coupling relationship among the mechanical parameters of phosphogypsum. As the moisture content increases, water molecules enter the pore spaces between gypsum crystals, weakening the cementing force between crystals and leading to a decrease in cohesion. Simultaneously, the water filling the pores increases the mass per unit volume, resulting in increased bulk density. This moisture-driven change in cohesion and bulk density forms the physical basis for the correlation between phosphogypsum parameters. Modeling this joint correlation structure using the Copula function offers the advantage of separating the marginal distribution of each parameter from the correlation structure between them—the marginal distribution function independently characterizes the univariate statistical features of each parameter, while the Copula function specifically describes the dependencies between parameters. Both are unified into a complete joint distribution model through Sklar's theorem.

[0007] Regarding the choice of Copula function type, this invention adopts the t-Copula function. The technical reason is that the safety assessment of phosphogypsum slopes is most concerned with extreme conditions—namely, the parameter combination of a surge in water content and a simultaneous sharp drop in cohesion and internal friction angle after a heavy rain. This extreme combined deterioration corresponds to the tail region of the parameter joint distribution. The t-Copula function has symmetrical tail correlation, which can accurately characterize the joint probability of multiple parameters simultaneously taking extreme values ​​under extreme conditions. In contrast, the commonly used Gaussian Copula function is asymptotically independent at the tail, systematically underestimating the probability of extreme combined deterioration events, and thus underestimating the risk of instability.

[0008] The mathematical expression for the t-Copula function is:

[0009] Where R is a 4×4 correlation matrix, and ν is the degree of freedom parameter. Let be a multivariate t-distribution function with ν degrees of freedom and R correlation matrix. Let f be the inverse function of a univariate t-distribution with ν degrees of freedom. , , , These are the probability integral values ​​of cohesion c, internal friction angle φ, unit weight γ, and moisture content w after transformation by their respective marginal distribution functions. The degree of freedom parameter ν controls the strength of the tail correlation—the smaller ν is, the stronger the tail correlation; as ν approaches infinity, the t-Copula degenerates into a Gaussian Copula.

[0010] The marginal distribution of each parameter is selected based on the following criteria: cohesion c adopts a log-normal distribution. Because cohesion is always positive and its probability density is right-skewed; the internal friction angle φ follows a normal distribution. Because the variation range of the internal friction angle is relatively concentrated and its distribution is approximately symmetrical; the unit weight γ follows a normal distribution. Moisture content w is distributed using a Beta distribution. Because the moisture content is limited to the range of 0 to saturated moisture content, the Beta distribution can accurately describe the probabilistic characteristics of bounded variables.

[0011] Furthermore, after acquiring the field-measured data of the target stockpile, the Copula parameters are adaptively calibrated in the field through Bayesian updating. Specifically, the correlation matrix R and degrees of freedom ν estimated from literature statistics are used as the prior distribution, and the likelihood function of the field-measured data is used for updating to obtain the posterior parameters (R*, ν*). The technical significance of this Bayesian updating mechanism is that when the amount of field data is small, the posterior parameters are mainly dominated by prior information, ensuring the stability of the model; when the amount of field data increases, the posterior parameters gradually converge with the field data, improving the model's site specificity.

[0012] Step S2: Perform Copula joint sampling based on the joint distribution model to generate a synthetic parameter sample set that satisfies the physical correlation constraints between parameters. Calculate the corresponding safety factor label for each group of parameters in the synthetic parameter sample set using the finite element strength reduction method to form a synthetic training dataset.

[0013] Generated by sampling from the calibrated t-Copula model The specific process of sample grouping parameters is as follows: First, sample vectors are drawn from a multivariate t-distribution with degrees of freedom ν* and a correlation matrix R*. Then, it is transformed into a uniformly distributed sample using a univariate t-distribution function. Finally, the parameters are transformed into physical parameters through the inverse function of the marginal distribution of each parameter: , , , Each parameter sample generated thus satisfies the physical correlation constraints between phosphogypsum parameters.

[0014] For slope α, slope height H, and pore water pressure ratio These three parameters are slope geometry and hydrological parameters, not inherent material mechanical properties. Their values ​​are determined by engineering design and environmental conditions. There is no physical coupling between the parameters; therefore, they are sampled independently within their respective reasonable ranges for each project, evenly distributed. The slope angle α ranges from 15° to 45°, the slope height H ranges from 10m to 80m, and the pore water pressure ratio... The value range is from 0 to 0.5.

[0015] Substitute the parameters of each combination into the finite element strength reduction method to calculate the safety factor, forming a synthetic training dataset. ,in Let i be the parameter vector of the i-th group. This corresponds to the safety factor. Number of samples. The value is determined based on the parameter space dimension and accuracy requirements, and is usually between 2000 and 5000.

[0016] The technical advantage of Copula joint sampling in generating training sets compared to traditional Latin hypercube sampling lies in the following: Traditional methods treat each parameter as an independent variable and sample them separately, which may generate physically impossible parameter combinations, such as "high water content and high cohesion"—in actual phosphogypsum, high water content inevitably leads to weakened cementation and reduced cohesion, and the combination of both being high violates physical laws. Such unreasonable training samples occupy the network's fitting capacity, diluting the model's prediction accuracy within the physically reasonable parameter domain. Copula joint sampling fundamentally eliminates unreasonable parameter combinations, concentrating all training capacity within the physically possible parameter domain, thus improving the effective approximation accuracy of the Extreme Learning Machine.

[0017] Step S3: Construct and train the Extreme Learning Machine network: Set the input layer nodes, hidden layer nodes, and output layer nodes of the Extreme Learning Machine, randomly generate the hidden layer input weight matrix and bias vector and keep them unchanged, calculate the hidden layer output matrix based on the synthetic training dataset, and solve the output weight vector through the generalized inverse matrix.

[0018] Extreme Learning Machine (ELM) is a single-hidden-layer feedforward neural network. Its core characteristic lies in the fact that the hidden layer parameters (input weight matrix W and bias vector b) are randomly generated during network initialization and remain unchanged throughout the entire usage cycle. Only analytical computation is needed to solve for the output layer weights. This feature offers two technical advantages over traditional backpropagation neural networks: first, the training process eliminates the need for gradient iteration, requiring only a single matrix inversion operation, thus increasing training speed by several orders of magnitude; second, because the hidden layer parameters are fixed, the analytical solution for the output weights is the globally optimal solution, avoiding the local optima and gradient vanishing problems of the backpropagation algorithm. More importantly, the fixed hidden layer parameters allow for the analytical computation of the partial derivatives of the ELM's output function with respect to the input parameters. This property forms the mathematical basis for the joint safety margin back-analysis in step S5—only models with analytical gradients can support efficient first-order reliability analysis iterations.

[0019] The network structure design of the Extreme Learning Machine is as follows: The input layer includes 7 nodes, corresponding to cohesion c, internal friction angle φ, unit weight γ, water content w, slope α, slope height H, and pore water pressure ratio, respectively. The hidden layer consists of L nodes, with the value of L optimized within the range of 50 to 500 using 10-fold cross-validation to minimize the root mean square error of the cross-validation. The output layer consists of one node, corresponding to a safety factor Fs. The sigmoid function is used as the activation function for the hidden layers. .

[0020] The mathematical description of the training process is as follows: Randomly generate the input weight matrix. and bias vector ,in Represents the set of real numbers. For the synthetic training dataset... Given a set of input data, calculate the hidden layer output matrix. Its element in the l-th row and i-th column is ,in Let W be the l-th row vector. Output the weight vector. Solving using the Moore-Penrose generalized inverse:

[0021] in This is the safety factor label vector.

[0022] Step S4: Use the on-site measured data of the target storage yard to perform transfer correction on the output weight vector: Substitute the on-site measured data into the Extreme Learning Machine network to calculate the on-site hidden layer output matrix, and update the output weight vector through the regularized incremental correction formula to obtain the corrected Extreme Learning Machine model.

[0023] The technical motivation for transfer correction lies in the fact that although the synthetic training data in step S2 ensures the physical consistency of the parameter combination through Copula modeling, its statistical characteristics are derived from literature-synthesized data, and there are still site-specific biases between it and the specific target stockpile. For example, the phosphorus content of phosphogypsum in a certain stockpile may be particularly high, resulting in a lower mean cohesion than the literature average. If an extreme learning machine trained on synthetic data is directly used to predict the stockpile, systematic bias will occur.

[0024] The mathematical principle of migration correction is as follows: Suppose the target storage yard has Group of on-site measured data ( (The number of sets can be small, for example, 5 to 15 sets). Substitute these data into the Extreme Learning Machine network, and use the same W and b as in step S3 to calculate the on-site hidden layer output matrix. The output weight vector is updated using a regularized incremental correction formula:

[0025] in, Let λ be the on-site measured safety factor vector, λ be the regularization coefficient, and I be the... An identity matrix of order 1. In the formula ( ) represents the residual vector of the source domain model on the target domain data, and the correction amount. It is exactly equal to the minimum norm adjustment required to eliminate the residual. Regularization term. Its purpose is to prevent overfitting due to excessive adjustments under small sample conditions—the larger λ is, the more conservative the adjustments. The closer to the source domain training result The smaller λ is, The closer it gets to the on-site data. When Approaching zero time Degenerate into ,when When λ is sufficiently large and λ approaches zero It approaches the results of training in a pure target domain.

[0026] The regularization coefficient λ is obtained by... The results were determined using leave-one-out cross-validation on the field measured data. Specifically, one set of data was left out as the validation set each time, and the rest... The group is used as the calibration set. Candidate λ values ​​are iterated through, and those values ​​are selected such that... The value of λ that minimizes the sum of squared predicted residuals in leave-one-out verification.

[0027] The synergistic effect of migration correction and Copula joint distribution modeling is that: in step S1, the Bayesian update of Copula uses field data to adjust the joint distribution parameters, ensuring that the parameter distribution of the synthetic training set is close to the target stockpile; this makes the parameters trained on the synthetic data in step S3... The baseline deviation is already small, and only a small correction is needed in step S4 to eliminate the residual deviation. Therefore, even with very little field data (5 to 15 sets), the correction will not overfit – because the signal amplitude that needs to be corrected is limited, and the regularization constraint is sufficient to ensure stability.

[0028] When making a positive prediction of the current parameter state of the target stockpile, the current parameter vector is input. Calculate the predicted safety factor:

[0029] The single prediction time is less than 1 millisecond, of which This is the predicted value for the safety factor. For output weight vector The l-th component.

[0030] Step S5: Perform a joint safety margin back analysis on the current parameter status of the target stockpile.

[0031] Traditional safety assessments only output a single value for the safety factor, failing to answer core engineering decision-making questions such as "how much safety margin is left before instability" and "which parameters' joint degradation is most likely to trigger instability." This step upgrades the function from "deterministic safety factor assessment" to "probabilistic joint safety margin assessment" by performing first-order reliability analysis in the joint parameter space defined by Copula.

[0032] Step S5.1, Nataf Transformation. This transforms the mechanical parameters from physical space to standard normal space. Specifically, for each mechanical parameter... (i=1,2,3,4 correspond to c,φ,γ,w respectively), and are transformed into standard normal variables using the following formula. :

[0033] in Let i be the marginal distribution function of the i-th parameter calibrated in step S1. It is the inverse function of the standard normal distribution. The correlation matrix R* of Copula corresponds to the equivalent correlation matrix in the standard normal space under the Nataf transform. , The solution is obtained numerically by transforming the relationship with equal probability.

[0034] The physical significance of the Nataf transform lies in the fact that, in the transformed standard normal space, the joint distribution of parameters is a multivariate standard normal distribution. Its isodensity surface is an ellipsoid centered at the origin. The farther a point is from the origin, the more extreme the corresponding parameter combination (the lower the joint probability density). Therefore, the origin represents the most likely state of parameter values, while the Euclidean distance from the origin measures the degree to which the parameter deviates from the most likely state.

[0035] For slope α, slope height H, and pore water pressure ratio These parameters are treated as deterministic known quantities in the safety margin inverse analysis, do not participate in the Nataf transformation, and their values ​​are fixed as the current actual values.

[0036] Step S5.2, construct the limit state function. It is defined as follows:

[0037] in, Let be the prediction function of the corrected extreme learning machine model obtained in step S4. This is the inverse Nataf transform (transforming U in standard normal space back to X in physical space). The target safety factor threshold is set. The surface where the limit state function value is zero. The standard normal space is divided into a safe region (G>0, where the safety factor is greater than the target value) and an unstable region (G<0, where the safety factor is less than the target value).

[0038] Target safety factor threshold According to the current specifications, the value is set as follows: 1.30 for normal operating conditions, 1.20 for flood operating conditions, and 1.10 for earthquake operating conditions.

[0039] Step S5.3: Calculate the gradient of the limiting state function with respect to the standard normal space parameters. Gradient calculation is a core step in first-order reliability analysis. This invention utilizes the structural characteristic of the limited learning machine with fixed hidden layer parameters to directly calculate the gradient using analytical methods, avoiding the approximation errors and additional computational overhead of numerical differencing.

[0040] Extreme learning machine for the j-th physical space parameter The partial derivatives are calculated analytically as follows:

[0041] Where L is the number of hidden layer nodes. The corrected output weight vector The l-th component, The derivative of the Sigmoid function and , Let W be the l-th row vector of the hidden layer input weight matrix. Let be the l-th component of the bias vector b. for The j-th element, X, is the current physical space parameter vector. Since W and b are randomly fixed and known during training, The solution has been obtained in step S4, therefore the above partial derivatives are analytic functions of X.

[0042] Transform the gradient in physical space into its Jacobian matrix J using the Nataf transform, and then convert it into its gradient in standard normal space.

[0043] in represents the Jacobian matrix element of the inverse Nataf transform.

[0044] The key technical reason for choosing Extreme Learning Machine (ELM) as the surrogate model over other machine learning methods is revealed in this step: the hidden parameters of backpropagation neural networks are iteratively updated during training, and their gradient calculation requires numerical difference approximation, which not only introduces truncation error but also requires additional model calls for each gradient evaluation; models such as Support Vector Machines (SVM) and Random Forests lack continuous differentiability and cannot be directly used for gradient-driven iterative search. The fixed hidden parameters of ELM make it the only machine learning model that can simultaneously provide fast prediction and analytical gradients, enabling HL-RF (Hasofer-Lind and Rackwitz-Fiessler) iterations to converge in milliseconds.

[0045] Step S5.4: The HL-RF iterative algorithm is used to search for the most likely instability point in the standard normal space. The iterative formula of the HL-RF algorithm is:

[0046] in Let the coordinates be the standard normal space coordinates of the k-th iteration. For the limit state function in The gradient vector at that point, The initial point u0 is taken as the standard normal space coordinates of the current parameter state X0 after Nataf transformation. The iterative convergence criterion is that the change in the reliability index between two adjacent iterations is less than a preset threshold and the absolute value of the limit state function value is less than a preset threshold. Due to the smooth and continuous analytical gradient of the Extreme Learning Machine, the HL-RF iteration usually converges within 3 to 8 steps, with a total time of less than 10ms.

[0047] The physical meaning of the most likely instability point u* is: the point at which the safety factor is exactly equal to the target value. In the parameter combinations, u* is the most likely combination under the Copula joint distribution. In other words, the parameter degradation path corresponding to u* is the "shortest probability distance path" from the current safe state to the instability boundary.

[0048] Step S5.5: Extract the joint safety margin index system based on the most likely instability point.

[0049] Reliability index β: β is the Euclidean distance between the most likely point of instability and the origin of the standard normal space. The larger β is, the farther the current parameter state is from the instability boundary, and the more sufficient the safety margin.

[0050] Probability of destruction : ,in It is a standard normal distribution function. The probability of failure is a first-order approximation of the probability that the safety factor is lower than the target value.

[0051] Parameter sensitivity factor : . The square value Let represent the contribution rate of the i-th parameter to the instability risk, and the contribution rate of all parameters. The sum equals 1. Since the sensitivity factor is calculated in the standard normal space after Copula transformation, it automatically includes the effect of parameter correlation—when cohesion and moisture content are highly correlated, the sensitivity factor does not reflect the independent contributions of the two, but rather the net contribution of each to the joint instability risk under the condition of maintaining relevant constraints.

[0052] Step S5.6: Map the most likely instability point u* back to physical space using the inverse Nataf transform to obtain the joint critical values ​​of each parameter. The difference between the current parameter value and the joint critical value. This refers to the combined permissible variation of all parameters. The parameter with the largest squared value of the sensitivity factor is determined as the most sensitive parameter, and the combined permissible variation of this parameter is used as its dynamic early warning threshold.

[0053] The fundamental difference between the joint allowable change and the traditional independent analysis method lies in the fact that the traditional method analyzes "how much further the cohesion can decrease" or "how much further the moisture content can increase" independently along the coordinate axes of each parameter, and the allowable changes of each parameter are unrelated. In contrast, the joint allowable change of this invention reflects the joint degradation path under Copula-related constraints—when the cohesion degrades along Δc, the moisture content simultaneously changes along Δw. This joint degradation path is determined by the relevant structure of Copula and is the shortest path in a probabilistic sense to reach the instability boundary.

[0054] Beneficial Effects: The rapid assessment method for the stability of phosphogypsum slopes based on extreme learning machine provided by this invention has the following beneficial effects: First, by using Copula joint distribution modeling and joint sampling to generate training sets, unreasonable parameter combinations generated by traditional independent sampling are eliminated, allowing the entire fitting capacity of the extreme learning machine to be concentrated in the physically possible parameter domain, thus improving prediction accuracy with the same number of hidden layer nodes.

[0055] Second, by using output weight regularization incremental migration correction, the ability to adapt a general model to a specific stockpile is achieved with only 5 to 15 sets of field data, solving the problem of model availability under small sample conditions commonly faced by phosphogypsum stockpiles.

[0056] Third, by utilizing the structural characteristic of fixed hidden layer parameters of the extreme learning machine to realize analytical gradient calculation, the HL-RF iteration in the standard normal space is driven to converge in milliseconds. Simultaneously, five types of indicators are output: safety factor, reliability index, probability of failure, parameter sensitivity factor, and joint allowable variation. This realizes the functional upgrade from deterministic safety factor assessment to probabilistic joint safety margin assessment.

[0057] Fourth, the joint allowable variation of each parameter output by the joint safety margin back analysis can be directly converted into a dynamic early warning threshold that takes into account the correlation of parameters, which overcomes the shortcomings of traditional fixed thresholds that fail to report under conditions of coordinated deterioration of parameters or false alarms when the safety margin is sufficient. Attached Figure Description

[0058] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below.

[0059] Figure 1 This is a flowchart illustrating the overall process of a rapid assessment method for phosphogypsum slope stability based on extreme learning machine provided in this invention. The flowchart shows the complete technical chain from steps S1 to S5, including: S1 constructing a Copula joint distribution model (receiving 50-200 sets of literature statistical data as priors), S2 Copula joint sampling and batch FEM calculation, S3 extreme learning machine training, S4 field data transfer and correction (receiving 5-15 sets of field measured data), and the forward prediction stage (…). Single prediction time <1ms), and S5 joint safety margin back analysis (including six sub-steps S5.1~S5.6), the final output / / / Five categories of indicators (end-to-end <20ms).

[0060] Figure 2 This is a flowchart illustrating the Copula joint distribution modeling and training set generation process provided in this embodiment of the invention. The flowchart shows the detailed process of steps S1 and S2, including: literature statistics (50-200 samples) are processed through two parallel paths: marginal distribution fitting (c~LogN, φ~N, γ~N, w~Beta) and t-Copula parameter estimation (maximum likelihood method → ​​R,ν). These, along with field measured data (5-15 sets), are then updated using Bayesian methods to obtain the calibrated parameters θ*=(R*,ν*). Finally, Copula joint sampling (simultaneously receiving operating condition parameters α, H, ...) is performed. The uniform sampling is used to generate the synthetic training dataset Ds through batch computation of FEM (intensity reduction method → ​​Fs label).

[0061] Figure 3This is a schematic diagram of the network structure of the Extreme Learning Machine provided in an embodiment of the present invention; the diagram shows the three-layer network structure of the Extreme Learning Machine: the input layer contains 7 nodes (corresponding to cohesion c, internal friction angle φ, unit weight γ, water content w, slope α, slope height H and pore water pressure ratio, respectively). The hidden layer contains L nodes (g1, g2, g3, ..., g...). L The output layer contains one node (safety factor Fs). The diagram indicates that the input weight matrix W and bias vector b are randomized fixed parameters, and the output weights... These are trainable parameters. The activation function is labeled at the bottom. and the formula for calculating the analytical gradient .

[0062] Figure 4 The flowchart illustrates the joint safety margin inverse analysis provided in this embodiment of the invention; the flowchart shows the detailed process of step S5, including six sub-steps: based on the current parameter state Starting from the S5.1Nataf transform ( Mapped to the standard normal space, the limit state function is constructed via S5.2. Analytical gradient calculation using S5.3ELM ( After Jacobian matrix transformation, and through S5.4HL-RF iterative search ( →Converging to u*, including iterative loops between S5.3 and S5.4), and then extracting the index via S5.5 ( , , Finally, the combined allowable change ΔX is calculated by S5.6 and the dynamic early warning threshold is output.

[0063] Figure 5 This diagram illustrates the search for the most probable instability point in the standard normal space provided in this embodiment of the invention. Using a two-dimensional standard normal space (u1, u2) as an example, it shows that the limiting state surface G(U) = 0 divides the space into a safe region (G(U) > 0, located on one side of the origin) and an instability region (G(U) < 0). The current parameter state u0 is located within the safe region. The HL-RF iteration path starts from u0 and converges iteratively to the most probable instability point u* on the limiting state surface. The reliability index β is equal to the Euclidean distance ||u|| from u to the origin O, and the gradient vector... The direction points towards the instability region. The figure also marks the isoprobability density ellipsoid (dashed arc). Detailed Implementation

[0064] The present invention will now be described more clearly and completely by way of a preferred embodiment in conjunction with the accompanying drawings, but this does not limit the invention to the scope of the described embodiment.

[0065] Example 1 This example focuses on a phosphogypsum stockpile of a phosphate chemical company. The stockpile has been in operation for eight years, currently has a slope height of approximately 45 meters, a gradient of approximately 28 degrees, and is equipped with a three-level horse trail. The area receives an average annual rainfall of approximately 1200 mm, with the flood season concentrated between June and September. The stockpile has been equipped with 12 pore water pressure monitoring sensors and 8 sets of moisture content monitoring probes, providing 12 sets of valid mechanical test data available on-site.

[0066] The implementation process of step S1 (e.g.) Figure 1 and Figure 2 (as shown) A total of 156 sets of physical and mechanical test data for phosphogypsum were collected from publicly available literature. The statistical characteristics of each parameter are as follows: the mean cohesion c is 28.5 kPa, and the standard deviation is 11.2 kPa; the mean internal friction angle φ is 29.8°, and the standard deviation is 3.6°; the mean bulk density γ is 15.8 kN / m³. 3 The standard deviation is 1.4 kN / m. 3 The mean moisture content w is 0.32, and the standard deviation is 0.08.

[0067] The marginal distribution fitting results are as follows: c follows a LogNormal (3.28, 0.38) distribution, φ follows a Normal (29.8, 3.6) distribution, γ follows a Normal (15.8, 1.4) distribution, and w follows a Beta (9.6, 20.4) distribution. The goodness of fit of each distribution was confirmed by the KS test, and the significance level was greater than 0.05.

[0068] The t-Copula parameters are estimated using the maximum likelihood method, and the key elements in the correlation matrix R are: (A strong negative correlation between cohesion and moisture content). (Positive correlation between bulk density and moisture content) (A weak positive correlation between cohesion and the angle of internal friction), degrees of freedom .

[0069] Bayesian updates were performed using 12 sets of field measurement data from the stockyard to obtain the calibrated parameters: , , After the update The absolute value of the value decreased slightly, indicating that the negative correlation between the cohesion and moisture content of the phosphogypsum in this stockpile was slightly weaker than the average level in the literature. This is consistent with the engineering practice of improving the density of the stockpile due to the use of higher compaction work.

[0070] The implementation process of step S2 (e.g.) Figure 2 (as shown) 3000 sets of mechanical parameter samples were generated from the calibrated t-Copula model. Slope α was sampled uniformly within the range of 15° to 45°, slope height H was sampled uniformly within the range of 10m to 80m, and pore water pressure ratio was also sampled. Uniform sampling was performed within the range of 0 to 0.5. Each set of parameters was then substituted into the SIGMA / W+ module of GeoStudio software, and the safety factor was calculated using the finite element strength reduction method, resulting in a synthetic training dataset containing 3000 records.

[0071] The parameters of the synthetic training set were checked by scatter plot. It was confirmed that there were no unreasonable parameter combinations such as "high moisture content - high cohesion". Cohesion and moisture content showed a clear negative correlation trend, which was consistent with the physical expectation.

[0072] The implementation process of step S3 (e.g.) Figure 3 (as shown) An extreme learning machine network was constructed, with 7 nodes in the input layer, 280 nodes in the hidden layer determined by 10-fold cross-validation (searching stepwise with a step size of 50 within the range of 50 to 500 nodes; the root mean square error of cross-validation was minimized to 0.042 with 280 nodes), and 1 node in the output layer. The input weight matrix was randomly generated. and bias vector The hidden layer output matrix is ​​generated using a uniform distribution within the range [-1, 1]. The weight vector is output by solving the generalized inverse. .

[0073] Coefficient of determination on the training set The root mean square error (RMSE) was 0.038; on the independently plotted test set of 600 lines, =0.961, RMSE=0.052, indicating good generalization performance of the model.

[0074] The implementation process of step S4 (e.g.) Figure 1 (as shown) Using 12 sets of field measurement data from the stockpile, the field hidden layer output matrix was calculated. The regularization coefficient λ = 0.015 was determined through leave-one-out cross-validation. Substituting this into the incremental correction formula yields... .

[0075] Before calibration, the root mean square error of the source domain model on 12 sets of field data was 0.098; after calibration, it decreased to 0.031, with a relative improvement rate of 68.4%. The relative error of the calibration model in predicting the safety factor of the stockyard is less than 3%, which meets the engineering accuracy requirements.

[0076] The implementation process of step S5 (e.g.) Figure 1 , Figure 4 and Figure 5(as shown) The current measured parameters of the stockpile are: c0 = 24.3 kPa, φ0 = 28.5°, γ0 = 16.2 kN / m. 3 , w0=0.35, α0=28°, H0=45m, r u0 =0.18. Set the target safety factor F. target =1.30.

[0077] Positive prediction result: Fs predict =1.42, indicating that the current state is safe.

[0078] like Figure 4 and Figure 5 As shown, a joint safety margin inverse analysis is performed: (c0,φ0,γ0,w0) is mapped to the standard normal space via Nataf transformation to obtain u0=(-0.38,-0.36,0.29,0.38). HL-RF iteration is performed with u0 as the initial point. After 6 iterations, convergence is achieved, yielding the most likely instability point u*=(-1.52,-0.98,0.83,1.41), with a convergence accuracy ||G(u*)||<10. -6 .

[0079] The reliability index β = ||u*|| = 2.47, corresponding to the probability of failure P. f =Φ(-2.47)=6.76×10 -3 That is, the probability that the safety factor is lower than 1.30 under the current state is about 0.68%.

[0080] Sensitivity factors for each parameter: (Cohesion) =0.33 (moisture content) =0.16 (angle of internal friction). =0.13 (bulk density). The combined contribution of cohesion and moisture content reaches 71%, making them the parameter pair that most needs to be monitored.

[0081] Transforming u* back to physical space: c* = 16.8 kPa, φ* = 25.7°, γ* = 17.1 kN / m 3 w*=0.42. The combined allowable variation of all parameters is: =7.5kPa, =2.8°, =0.9kN / m 3 , =0.07. That is, when the moisture content increases from the current 0.35 to 0.42, and the cohesion decreases from 24.3 kPa to 16.8 kPa, the safety factor will drop to the critical value of 1.30. The moisture content monitoring and early warning threshold for this stockpile should be dynamically set to 0.42, and the cohesion monitoring and early warning threshold should be dynamically set to 16.8 kPa.

[0082] The entire calculation time for step S5 is 14 milliseconds, including 2 milliseconds for Nataf transformation, 11 milliseconds for HL-RF iteration (including 6 ELM positive evaluations and gradient calculations), and 1 millisecond for index extraction.

[0083] Comparative verification: with 10 6 Using the Monte Carlo simulation as a baseline, 10 samples were randomly drawn from the Copula model. 6 The safety factor was calculated by substituting each set of parameter samples into the corrected extreme learning machine model. The proportion of samples with Fs < 1.30 was 7.02 × 10⁻⁶. -3 The failure probability P obtained from the first-order reliability analysis of this invention f =6.76×10 -3 The relative deviation from the Monte Carlo results was 3.7%, validating the accuracy of the method. The Monte Carlo simulation took approximately 15 minutes (10... 6 The sublimity learning machine evaluation, while the joint safety margin back analysis of this invention only requires 14 milliseconds, improving computational efficiency by approximately 6.4 × 10⁻⁶. 4 times.

[0084] Example 2 This embodiment verifies the performance of migration correction under different amounts of field data. Using the same source domain model as in Embodiment 1 (steps S1 to S3), n... t = Transfer correction was performed on 5, 8, 10, 12, and 15 sets of field data to evaluate the prediction accuracy of the corrected model on the independent validation set.

[0085] The coefficient of determination R of the calibrated model on the validation set under various data sizes 2 They are respectively: n t =5 when R 2 =0.891, n t =8 times R 2 =0.923, n t =10 times R 2 =0.948, n t =12 times R 2 =0.961, n t =15 hours R 2 =0.972. As a control, the source domain model trained directly using LHS sampling without Copula quality assurance achieved the same result under the same conditions. t =5 when R 2=0.782, n t =12 times R 2 =0.918. The comparison shows that Copula's quality assurance improves the accuracy of transfer correction by approximately 10-14% under small sample conditions.

[0086] Example 3 This embodiment verifies the advantages of joint safety margin back analysis compared to independent analysis. For the stockpile in Embodiment 1, the permissible variations of each parameter are calculated using both the joint back analysis of this invention and the traditional independent back analysis.

[0087] Independent back analysis: With the remaining parameters fixed at their current values, independently solve for the critical value that makes Fs = 1.30 for each parameter. The results are: =14.2 kPa, =10.1 kPa; =0.46, =0.11.

[0088] Joint back analysis (this invention): =16.8kPa, =7.5kPa; =0.42, =0.07.

[0089] The permissible variation given by the joint back analysis is significantly smaller than that of the independent analysis. This is because the independent analysis assumes that only the moisture content increases while the cohesion remains constant. However, in phosphogypsum, this assumption violates physical laws—an increase in moisture content inevitably leads to a simultaneous decrease in cohesion. The joint back analysis considers this synergistic deterioration effect, providing a more conservative and reasonable warning threshold. If the lenient threshold of the independent analysis (moisture content 0.46) is used, when the actual combined degradation occurs—with the moisture content rising to 0.43 and the cohesion simultaneously decreasing to 17.5 kPa—the safety factor has already dropped to 1.28, below the target value, but the moisture content has not yet reached the independent threshold of 0.46, resulting in a false negative. The joint threshold of this invention (moisture content 0.42) can trigger a timely warning under this condition.

[0090] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A rapid assessment method for the stability of phosphogypsum slopes based on extreme learning machine, characterized in that, Specifically, the following steps are included: S1. Obtain statistical data on the cohesion, internal friction angle, bulk density and moisture content of phosphogypsum, fit the marginal probability distribution function of each parameter respectively, establish a joint distribution model of the four mechanical parameters using the t-Copula function, and estimate the correlation parameters of the joint distribution model by the maximum likelihood method. S2. Based on the joint distribution model, Copula joint sampling is performed to generate a synthetic parameter sample set that satisfies the physical correlation constraints between parameters. For each group of parameters in the synthetic parameter sample set, the corresponding safety factor label is calculated using the finite element strength reduction method to form a synthetic training dataset. S3. Set the input layer nodes, hidden layer nodes, and output layer nodes of the extreme learning machine. Randomly generate the hidden layer input weight matrix and bias vector and keep them fixed. Calculate the hidden layer output matrix based on the synthetic training dataset and solve for the output weight vector through the generalized inverse matrix. S4. Use the on-site measured data of the target stockpile to perform transfer correction on the output weight vector: Substitute the on-site measured data into the extreme learning machine network to calculate the on-site hidden layer output matrix, and update the output weight vector through the regularized incremental correction formula to obtain the corrected extreme learning machine model. S5. Perform joint safety margin back analysis on the current parameter state of the target storage yard: Based on the joint distribution model, transform the current parameters from the physical space to the standard normal space, construct the limit state function with the corrected limit learning machine model, use the analytical gradient driven iterative algorithm of the limit learning machine to solve the most likely instability point in the standard normal space, calculate the reliability index, failure probability and sensitivity factor of each parameter based on the most likely instability point, and output the joint safety margin assessment result.

2. The rapid assessment method for the stability of phosphogypsum slopes based on extreme learning machine according to claim 1, characterized in that, In step S1, the t-Copula function is expressed as: ， Where R is the four-dimensional correlation matrix, and ν is the degree of freedom parameter. Let be the multivariate t-distribution function with ν degrees of freedom. Let f be the inverse function of a univariate t-distribution with ν degrees of freedom. , , , These are the probability integral values ​​of cohesion, internal friction angle, unit weight, and moisture content after transformation by their respective marginal distribution functions.

3. The rapid assessment method for the stability of phosphogypsum slopes based on extreme learning machine according to claim 2, characterized in that, In step S1, the marginal distribution of cohesion adopts a log-normal distribution, the marginal distribution of internal friction angle adopts a normal distribution, the marginal distribution of bulk density adopts a normal distribution, and the marginal distribution of moisture content adopts a Beta distribution.

4. The rapid assessment method for the stability of phosphogypsum slopes based on extreme learning machine according to claim 1, characterized in that, Step S1 also includes a field adaptive calibration sub-step for Copula parameters: after obtaining the field measured data of the target stockpile, the correlation parameters estimated by the statistical data are used as the prior distribution, and Bayesian updates are performed using the likelihood function of the field measured data to obtain the calibrated correlation parameters.

5. The rapid assessment method for the stability of phosphogypsum slopes based on extreme learning machine according to claim 1, characterized in that, In step S2, the method for generating the synthetic parameter sample set is as follows: the probability integral value vector is obtained by sampling from the joint distribution model, and then the probability integral value vector is transformed into a physical parameter vector by the inverse function of the marginal distribution of each parameter; for slope, slope height and pore water pressure ratio, they are sampled independently within their respective reasonable engineering ranges according to uniform distribution.

6. The rapid assessment method for the stability of phosphogypsum slopes based on extreme learning machine according to claim 1, characterized in that, In step S3, the input layer of the extreme learning machine includes seven input nodes, corresponding to cohesion, internal friction angle, bulk density, water content, slope, slope height, and pore water pressure ratio, respectively; the activation function of the hidden layer is the Sigmoid function; the formula for solving the output weight vector is: ， in, The output weight vector is L×1. This is the hidden layer output matrix. Let L be the safety factor label vector, and L be the number of hidden layer nodes. This represents the number of training samples used in the synthesis.

7. The rapid assessment method for the stability of phosphogypsum slopes based on extreme learning machine according to claim 1, characterized in that, In step S4, the regularization increment correction formula is: ， in, Let L×1 be the corrected output weight vector. The output weight vector is L×1. This is the hidden layer output matrix in the field. This is the vector of the safety factor measured on-site. Let I be the regularization coefficient. An identity matrix of order 1. This represents the number of sets of data measured on-site.

8. The rapid assessment method for the stability of phosphogypsum slopes based on extreme learning machine according to claim 7, characterized in that, The regularization coefficient λ is determined by leave-one-out cross-validation on field measured data, wherein the objective function of the leave-one-out cross-validation is to minimize the sum of squared predicted residuals.

9. A rapid assessment system for the stability of phosphogypsum slopes based on extreme learning machines, characterized in that, include: The joint distribution modeling unit is used to obtain statistical data on the cohesion, internal friction angle, bulk density and moisture content of phosphogypsum. The marginal probability distribution functions of each parameter are fitted respectively. The joint distribution model of the four mechanical parameters is established by using the t-Copula function, and the correlation parameters of the joint distribution model are estimated by the maximum likelihood method. The training data generation unit is used to perform Copula joint sampling based on the joint distribution model, generate a synthetic parameter sample set that satisfies the physical correlation constraints between parameters, and calculate the corresponding safety factor label for each group of parameters in the synthetic parameter sample set using the finite element strength reduction method to form a synthetic training dataset. The Extreme Learning Machine Training Unit is used to construct the Extreme Learning Machine network. It randomly generates the hidden layer input weight matrix and bias vector and keeps them fixed. It calculates the hidden layer output matrix based on the synthetic training dataset and solves the output weight vector through the generalized inverse matrix. The migration correction unit is used to substitute the on-site measured data of the target stockpile into the Extreme Learning Machine network to calculate the on-site hidden layer output matrix, and update the output weight vector through the regularized incremental correction formula to obtain the corrected Extreme Learning Machine model. The joint safety margin inverse analysis unit is used to transform the current parameters from the physical space to the standard normal space based on the joint distribution model, construct the limit state function with the corrected limit learning machine model, and use the analytical gradient-driven iterative algorithm of the limit learning machine to solve for the most likely instability point in the standard normal space. Based on the most likely instability point, the reliability index, failure probability and sensitivity factors of each parameter are calculated, and the joint safety margin assessment results are output.

10. A rapid assessment system for the stability of phosphogypsum slopes based on extreme learning machine as described in claim 9, characterized in that, The joint safety margin inverse analysis unit is also used to map the most likely instability point back to the physical space through the Nataf inverse transformation to obtain the joint critical value of each parameter, calculate the difference between the current parameter value and the joint critical value as the joint allowable change of each parameter, and determine the most sensitive parameter based on the square value of the sensitivity factor of each parameter. The joint allowable change corresponding to the most sensitive parameter is used as the dynamic early warning threshold of that parameter.

Images