A preconditioned subset simulation method for efficient computation of groundwater pollution failure probability

Abstract

The application discloses a preconditioned subset simulation method for efficiently calculating groundwater pollution failure probability, and comprises the following steps: generating N groundwater aquifer Gaussian permeability coefficient fields through a sequential Gaussian method; substituting the N groundwater aquifer Gaussian permeability coefficient fields into a functional function of a groundwater pollutant transport model to calculate response values; arranging the response values in ascending order to obtain a seed sample and a threshold value; taking the seed sample as an initial sample of each Markov chain, and generating a candidate sample according to a preconditioned Crank-Nicolson technique; judging whether the candidate sample is within the failure threshold, accepting the candidate sample falling within the failure threshold, and rejecting the candidate sample falling outside the failure threshold; repeating the Markov chain sampling and sorting process until the threshold value is less than 0, and then terminating the sampling; and calculating the failure probability according to the sampled sample. When the method is used for processing a high-dimensional groundwater pollution small failure probability problem, the efficiency is higher than that of traditional Monte Carlo simulation and subset simulation, and the result is more stable.

Description

Technical Field

[0001] This invention belongs to the field of hydrological statistics technology, specifically relating to a precondition subset simulation method for efficiently calculating the probability of groundwater pollution failure. Background Technology

[0002] As groundwater pollution becomes increasingly serious, research on predicting groundwater solute transport is receiving more and more attention. However, the inherent heterogeneity of aquifers and the limitations of conceptual mathematical models make accurate prediction of pollutant transport very difficult.

[0003] Uncertainty is closely related to risk, and an important method for assessing the risk posed by uncertainty in a system is to calculate the system's failure probability. In the field of groundwater pollution, uncertainty may lead to an underestimation of pollutant concentrations. Therefore, the failure probability of groundwater pollution is defined as the probability that the pollutant concentration in a given area exceeds the allowable critical concentration.

[0004] Failure probability calculation methods mainly include the first-order reliability method (FORM), the response surface methodology (RSM), and numerical simulation methods. The first-order reliability method has high computational efficiency, and the response surface methodology can solve problems without a clearly defined explicit function. However, neither of these methods can handle high-dimensional, highly nonlinear groundwater flow and solute transport problems.

[0005] The most common numerical simulation method is Monte Carlo simulation, which offers high computational accuracy and can solve most failure probability problems. However, it incurs enormous computational costs when dealing with low failure probability problems. The basic idea of ​​subset simulation is to transform low failure probability events into several intermediate, higher failure probability events for calculation. By introducing a series of intermediate failure events, the low failure probability can be represented as a product of higher failure probabilities. Subset simulation is an efficient method for calculating failure probabilities, but traditional subset simulation methods use the Metropolis-Hastings algorithm for sampling. As the dimensionality of the problem increases, the acceptance rate of candidate samples obtained by the Metropolis-Hastings algorithm rapidly drops to zero. This results in low computational accuracy for high-dimensional problems using traditional subset simulation methods based on the Metropolis-Hastings algorithm. Summary of the Invention

[0006] Purpose of the invention: In order to overcome the shortcomings of the prior art, the present invention discloses a preconditioned subset simulation method for efficiently calculating the probability of groundwater pollution failure. This method, by coupling the preconditioned Crank-Nicholson technique and the subset simulation method, enables the subset simulation method to efficiently handle the problem of high-dimensional groundwater pollution failure probability.

[0007] Technical solution: The precondition subset simulation method for efficiently calculating the probability of groundwater pollution failure disclosed in this invention includes the following steps:

[0008] S1. Generate N Gaussian permeability coefficient fields for underground aquifers, set the conditional probability value p0, and assume a preconditional subset to simulate intermediate failure events F. j j = 1,...,M;

[0009] S2. Substitute the permeability coefficient field of the groundwater aquifer into the function of the groundwater pollutant transport model to calculate the corresponding response value. Arrange the response values ​​in ascending order to obtain the seed sample and threshold b1. Use the threshold to divide the failure domain.

[0010] S3. Use the seed sample as the initial sample for each Markov chain, generate candidate samples according to the preconditioning Kranknick Nicholson technique, and determine whether the candidate samples are within the failure threshold.

[0011] Accept candidate samples within the failure threshold and reject candidate samples outside the failure threshold. Through sampling, there are a total of N groundwater aquifer permeability coefficient field samples within the failure threshold.

[0012] S4. Substitute the N groundwater aquifer permeability field samples within the failure threshold in S3 into the function function of the groundwater pollutant transport model to calculate the corresponding response values. Sort the response values ​​in ascending order to obtain the seed sample and threshold b. j The failure domain is divided using a threshold.

[0013] S5. Repeat processes S3 and S4 until the Np0th response value is less than 0, at which point the termination condition is met, sampling stops, and the failure probability is calculated.

[0014] Furthermore, in S1, N Gaussian permeability coefficient fields {x} of underground aquifers are generated through sequential Gaussian generation. i :i=1,...,N}, the preconditioned subset simulates intermediate failure events F j ={x:g(x)<b j}, j=1,...,M, j is the intermediate failure event number, x is the permeability coefficient field of the groundwater aquifer, g(x) is the function function of the groundwater pollutant transport model, b j For a specific threshold.

[0015] Furthermore, the function of the groundwater pollutant transport model in S2 is g(x) = C * -C, where C * C represents the allowable critical concentration value for pollutants, and C represents the pollutant concentration at a certain point in the sensitive area.

[0016] The response values ​​are sorted in ascending order, and the Np0th response value is used as the threshold b1. Groundwater aquifer permeability field samples with response values ​​less than the threshold fall into the failure region. Groundwater aquifer permeability field samples {x} within the failure threshold are then considered... i (1) :i=1,...,Np0} as seed samples.

[0017] Furthermore, S3 uses the preconditioning Crank-Nicholson technique based on the current sample x. t Generate candidate sample x t ′

[0018]

[0019] In the formula, β is the jump factor and C is the covariance matrix;

[0020] Each Markov chain is generated With 1 candidate sample and Np0 seed samples, there are a total of N groundwater aquifer permeability coefficient field samples in the failure domain.

[0021] Furthermore, S4 substitutes the N groundwater aquifer permeability coefficient field samples within the failure threshold in S3 into the functional function g(x)=C of the groundwater pollutant transport model. * -C calculates the corresponding response value, sorts the response values ​​in ascending order, and uses the Np0th response value as the threshold b. j Delineate the failure domain and sample the permeability field of the underground aquifer within the failure threshold {x}. i (j) :i=1,...,Np0} as seed samples.

[0022] Furthermore, after sampling stops when the termination condition is met in S5, the number N of groundwater aquifer permeability coefficient field samples within the failure domain is statistically analyzed. FM Calculate the failure probability of the Mth intermediate event.

[0023] Finally, the failure probability of the target failure event is calculated.

[0024] Beneficial effects: Compared with the prior art, the advantages of the present invention are:

[0025] 1. While ensuring the accuracy of calculating the probability of groundwater contamination failure, the computational efficiency of this invention is greatly improved compared with traditional Monte Carlo simulation;

[0026] 2. Under the same computation time, the results of this invention are more accurate and more stable than those of traditional subset simulation methods. Attached Figure Description

[0027] Figure 1This is a schematic diagram showing the locations of pollution sources and sensitive areas in the embodiment;

[0028] Figure 2 This is a flowchart of the method of the present invention;

[0029] Figure 3 The relationship curve between the number of samples and the relative error in the Monte Carlo simulation;

[0030] Figure 4 The relationship curves between the average failure probability and the number of samples in each layer were calculated for 10 iterations for the subset simulation method (SS) and the preconditional Crank-Nicholson coupled subset simulation method (pCN-SS).

[0031] Figure 5 The coefficient of variation is the result calculated by the two methods. Detailed Implementation

[0032] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0033] The invention was verified in a saturated steady-state confined aquifer as a method for efficiently calculating the failure probability of high-dimensional groundwater pollution. The aquifer space was set as 2000[L]*2000[L], with a thickness of 20[L], and discretized into 20*20*1 cells, each cell being 100[L]*100[L]*20[L]. The east and west boundaries of the model were constant head boundaries with given heads of 300[L] and 50[L] respectively, and the north and south boundaries were impermeable boundaries.

[0034] In this embodiment, the permeability coefficient field is treated as the only uncertainty variable, and it follows a log-multiple Gaussian distribution. The mean value is generated by a sequential Gaussian procedure, which is ln[LT]. -1 The standard deviation is 1ln[LT]. -1 The maximum and minimum correlation lengths are 800[L] and 600[L], respectively. The aquifer is a log-permeability random field with an exponential variogram of 135 degrees and a dimension of 400. The initial concentration of the aquifer is 0[ML]. -3 Other groundwater flow and solute transport model parameters were set as follows: porosity 0.3 [-], longitudinal dispersion 3.0 [L], and the ratio of lateral dispersion to longitudinal dispersion 0.5. It was assumed that solute transport of pollutants only involved advection and dispersion processes, and that the entire solute transport process was unsteady. The total simulation time was 8000 [T].

[0035] like Figure 1 The watershed shown has a pollution source (e.g., a landfill) upstream, with a constant mass loading rate of 1000 [ML]. -3The continuous discharge of pollutants into aquifers threatens the groundwater environment in sensitive downstream areas (such as residential areas).

[0036] Calculation of pollutant transport equations using the MT3DMS groundwater solute transport model:

[0037]

[0038] In the formula: C is the pollutant concentration [ML] -3 ]; t is time [T]; θ is the divergence operator; θ is the porosity; D m Molecular dispersion coefficient [L] 2 T -1 ]; α is the dispersion tensor [L]; For the Laplace operator; q s The volumetric flow velocity per unit volume of the aquifer [T] -1 ];C s Source flux concentration; v is the velocity vector [LT] -1 ].

[0039] The steady-state head H is obtained by calculating the steady-state equation of groundwater flow using the MODFLOW groundwater flow model:

[0040]

[0041] Where: K is the permeability coefficient [LT] -1 W represents the source and sink volume per unit volume of the aquifer [LT]. -1 ].

[0042] pCN-SS: A Preconditional Crank Nicholson Coupled Subset Simulation Method

[0043] SS: Subset Simulation Method

[0044] MCS: Monte Carlo Simulation Method

[0045] For this embodiment, a failure event is defined as the pollutant concentration C at any point within a sensitive area exceeding the allowable critical concentration value C within a certain period of time. * According to relevant regulations, different pollutants and regions have different permissible critical concentration values. This embodiment is a simplification of the actual situation, using the permissible critical concentration value C. * Set to 3.5 [ML] -3 The pollutant concentration C in the sensitive area was calculated using a groundwater flow and solute transport model.

[0046] C = F(x)

[0047] In the formula: x is the permeability coefficient of the random field.

[0048] The system's function is:

[0049] g(x) = C * -C

[0050] like Figure 2 As shown, the precondition subset simulation method for efficiently calculating the probability of groundwater pollution failure according to the present invention includes the following steps:

[0051] S1. Generate N Gaussian permeability coefficient fields for underground aquifers, set the conditional probability value p0, and assume a preconditional subset to simulate intermediate failure events F. j j = 1,...,M;

[0052] S2. Substitute the permeability coefficient field of the groundwater aquifer into the function of the groundwater pollutant transport model to calculate the corresponding response value. Arrange the response values ​​in ascending order to obtain the seed sample and threshold b1. Use the threshold to divide the failure domain.

[0053] S3. Use the seed sample as the initial sample for each Markov chain, generate candidate samples according to the preconditioning Kranknick Nicholson technique, and determine whether the candidate samples are within the failure threshold.

[0054] Accept candidate samples within the failure threshold and reject candidate samples outside the failure threshold. Through sampling, there are a total of N groundwater aquifer permeability coefficient field samples within the failure threshold.

[0055] S4. Substitute the N groundwater aquifer permeability field samples within the failure threshold in S3 into the function function of the groundwater pollutant transport model to calculate the corresponding response values. Sort the response values ​​in ascending order to obtain the seed sample and threshold b. j The failure domain is divided using a threshold.

[0056] S5. Repeat processes S3 and S4 until the Np0th response value is less than 0, at which point the termination condition is met, sampling stops, and the failure probability is calculated.

[0057] In practice:

[0058] (1) Generate N Gaussian permeability coefficient fields {x} of underground aquifers using sequential Gaussian generation. i Given the subset i = 1, ..., N, and a conditional probability value p0 = 0.1, simulate intermediate failure events F using a preconditional subset. j ={x:g(x)<b j}, j=1,...,M, j is the intermediate failure event number, x is the permeability coefficient field of the groundwater aquifer, g(x) is the function function of the groundwater pollutant transport model, b j For a specific threshold.

[0059] (2) Substitute the permeability field of the aquifer into the function of the groundwater pollutant transport model to calculate the corresponding response value. The function is g(x) = C * -C, where C * The allowable critical concentration value for pollutants is taken as 3.5 [ML]. -3 [C] represents the pollutant concentration at a point in the sensitive area. The response values ​​are sorted in ascending order, and the Np0th response value is used as the threshold b1. Groundwater aquifer permeability field samples with response values ​​less than the threshold fall into the failure region. Groundwater aquifer permeability field samples {x} within the failure threshold are considered... i (1) :i=1,...,Np0} as seed samples.

[0060] (3) Using the seed sample from the previous layer as the initial sample for each Markov chain, the preconditioning Crank-Nicholson technique is applied based on the current sample x. t Generate candidate sample x t ′

[0061]

[0062] In the formula, β is the jump factor and C is the covariance matrix.

[0063] Determine whether a candidate sample is within the failure threshold. Accept candidate samples within the failure threshold and reject candidate samples outside the failure threshold.

[0064] Each Markov chain is generated With 1 candidate sample and Np0 seed samples, there are a total of N groundwater aquifer permeability coefficient field samples in the failure domain.

[0065] (4) Substitute the permeability field samples of N groundwater aquifers within the failure threshold into the function g(x)=C of the groundwater pollutant transport model. * -C calculates the corresponding response value, sorts the response values ​​in ascending order, uses the Np0th response value as a threshold to divide the failure domain, and classifies the samples {x} within the failure threshold. i (j) :i=1,...,Np0} as seed samples.

[0066] (5) After the sampling is stopped when the termination condition is met, count the number of samples of the permeability coefficient field of the underground aquifer within the failure domain. Calculate the failure probability of the Mth intermediate event.

[0067] Finally, the failure probability of the target failure event is calculated.

[0068] like Figure 3The figure shows the relationship between the number of samples and the relative error in the Monte Carlo simulation. When the number of samples reaches 1,000,000, the relative error is 3.9%, indicating that the Monte Carlo simulation results are reliable. At this point, the failure probability calculated by the Monte Carlo simulation is 6.56 × 10⁻⁶. -4 This value will be used as a reference value for the accurate results of this embodiment.

[0069] like Figure 4 The figure shows the relationship between the average failure probability and the number of samples per layer obtained from 10 calculations using the subset simulation method (SS) and the preconditioned Crank-Nicholson coupled subset simulation method (pCN-SS). It can be seen that both methods require a much smaller number of samples than Monte Carlo simulation, with pCN-SS showing significantly better results than SS. When the sample size per layer is set to 1000, the pCN-SS results converge, while the SS results remain unstable.

[0070] like Figure 5 As shown, the coefficients of variation for the results calculated by the two methods are displayed. It can be seen that the coefficient of variation for the pCN-SS result is smaller, which means that pCN-SS is more stable than SS.

[0071] The above analysis shows that the coupled preconditioning Crank Nicholson subset simulation method used in this invention can efficiently calculate the probability of high-dimensional groundwater contamination failure.

Claims

1. A precondition subset simulation method for efficiently calculating the probability of groundwater pollution failure, characterized in that, Includes the following steps: S1. Generate a Gaussian permeability coefficient field for an underground aquifer, set a conditional probability value p0, and assume a preconditional subset to simulate intermediate failure events F. j j=1,…,M; S2. Substitute the permeability coefficient field of the groundwater aquifer into the function of the groundwater pollutant transport model to calculate the corresponding response value. Arrange the response values ​​in ascending order to obtain the seed sample and threshold b1. Use the threshold to divide the failure domain. S3. Use the seed sample as the initial sample for each Markov chain, generate candidate samples according to the preconditioning Kranknick Nicholson technique, and determine whether the candidate samples are within the failure threshold. Accept candidate samples within the failure threshold and reject candidate samples outside the failure threshold. Through sampling, there are a total of N groundwater aquifer permeability coefficient field samples within the failure threshold. S4. Substitute the N groundwater aquifer permeability field samples within the failure threshold in S3 into the function function of the groundwater pollutant transport model to calculate the corresponding response values. Sort the response values ​​in ascending order to obtain the seed sample and threshold b. j The failure domain is divided using a threshold. S5. Repeat processes S3 and S4 until the Np0th response value is less than 0, at which point the termination condition is met, sampling stops, and the failure probability is calculated. In S1, N Gaussian permeability coefficient fields {x} of underground aquifers are generated using sequential Gaussian generation. i : i=1,…, N}, the preconditioned subset simulates intermediate failure events F j ={x:g(x)<b j }, j=1,…M, j is the intermediate failure event number, x is the permeability coefficient field of the groundwater aquifer, g(x) is the function function of the groundwater pollutant transport system, b j For the threshold; The function of the groundwater pollutant transport model in S2 is g(x) = C * -C, where C * C represents the allowable critical concentration value for pollutants, and C represents the pollutant concentration at a certain point in the sensitive area. The response values ​​are sorted in ascending order, and the Np0th response value is used as the threshold b1. Groundwater aquifer permeability field samples with response values ​​less than the threshold fall into the failure region. Groundwater aquifer permeability field samples {x} within the failure threshold are then considered... i （1） : i=1,…,Np0} as seed samples; S3 uses the preconditioning Crank Nicholson technique based on the current sample x. t Generate candidate sample x t ’ ; ； In the formula, β is the jump factor and C is the covariance matrix; Each Markov chain generates (1 / p0)-1 candidate samples, plus Np0 seed samples, so that there are a total of N A sample of the permeability coefficient field of an underground aquifer; S4 substitutes the N aquifer permeability coefficient field samples within the failure threshold in S3 into the function function g(x) = C of the groundwater pollutant transport model. * -C calculates the corresponding response value, sorts the response values ​​in ascending order, and uses the Np0th response value as the threshold b. j Delineate the failure domain and sample the permeability field of the underground aquifer within the failure threshold {x}. i （j） : i=1,…, Np0} as seed samples.

2. The precondition subset simulation method for efficiently calculating the probability of groundwater pollution failure according to claim 1, characterized in that, After sampling stops when the termination condition is met in S5, the number of samples of the permeability coefficient field of the underground aquifer within the failure domain is counted. Calculate the failure probability of the Mth intermediate event. ; Finally, the failure probability of the target failure event is calculated. .

Images