CN114707395A - Groundwater organic pollution source inversion method based on heuristic homotopy algorithm - Google Patents
Groundwater organic pollution source inversion method based on heuristic homotopy algorithm Download PDFInfo
- Publication number
- CN114707395A CN114707395A CN202111602792.6A CN202111602792A CN114707395A CN 114707395 A CN114707395 A CN 114707395A CN 202111602792 A CN202111602792 A CN 202111602792A CN 114707395 A CN114707395 A CN 114707395A
- Authority
- CN
- China
- Prior art keywords
- model
- pollution source
- homotopy
- optimization
- algorithm
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/06—Multi-objective optimisation, e.g. Pareto optimisation using simulated annealing [SA], ant colony algorithms or genetic algorithms [GA]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Artificial Intelligence (AREA)
- Computer Hardware Design (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention relates to a groundwater organic pollution source inversion method based on a hyperheuristic homotopy algorithm, which comprises the following steps: establishing a numerical model of underground water organic pollution multiphase flow migration according to observation data, and determining pollution source characteristics to be inverted, aquifer parameters and value ranges of variables in the model; preparing a training sample set; training a single machine learning substitution model of the multiphase flow numerical model; establishing a combined intelligent substitution model; establishing a nonlinear programming optimization model for inversion identification of pollution source characteristics and aquifer parameters, and coupling an integrated learning intelligent surrogate model into the nonlinear programming optimization model; solving the optimization model by a heuristic-homotopy algorithm, and inversely identifying the characteristic value of the pollution source and the parameter value of the aquifer to be solved; the method can solve the problem of inversion of the underground water organic pollution source and multiphase flow migration parameters, and provides technical support for accurate simulation and prediction of the spatial and temporal distribution of the organic pollutants in the underground water.
Description
Technical Field
The invention relates to the technical field of underground water numerical simulation inverse problems, in particular to an underground water organic pollution source inversion method based on a heuristic homotopy algorithm.
Background
Organic pollutants have the characteristics of low water solubility, high toxicity, high interfacial tension and the like, are gathered and retained at the top or the bottom of an aquifer after entering a groundwater system (the gathering position depends on the density of the pollutants is less than that of water or more than that of the water), and are continuously dissolved and released into the water in the process of contacting with the water, so that serious and durable pollution is caused. In the actual remediation process of organic pollution, the difficulties of low pollutant removal rate, long time consumption in the remediation process and high remediation cost are often faced. Therefore, it is very important to make a reasonable and efficient repair scheme. The establishment of a reasonable and efficient remediation scheme is premised on the identification and mastery of the pollution source conditions in the aquifer.
However, underground water is buried in underground rock-soil media, and underground water pollution usually has the characteristics of concealment and discovery hysteresis, so that people are lack of understanding and mastering the conditions of underground water pollution sources, and a plurality of key parameters related to the migration process of organic pollutants in an underground water system are difficult to directly obtain by the conventional measurement means. Aiming at the problem, the inversion calculation by utilizing actual observation data and combining a data assimilation method is a main solution at present.
Disclosure of Invention
The invention aims to provide a ground water organic pollution source inversion method based on a heuristic homotopic algorithm, which can be used for identifying the characteristics of the ground water organic pollution source and estimating key parameters in a multiphase flow migration model, so that the accurate simulation and prediction of the ground water organic pollution multiphase flow migration process are realized, and important precondition and basic conditions are provided for the reasonable formulation of a ground water pollution remediation scheme, pollution risk evaluation and pollution responsibility determination.
In order to solve the technical problems, the technical scheme of the invention is as follows: 1. an underground water organic pollution source inversion method based on a heuristic homotopy algorithm is characterized by comprising the following steps: the method comprises the following steps:
s1: mastering hydrogeological conditions of a polluted site by means of field on-site investigation and dynamic monitoring to obtain groundwater water quality dynamic monitoring data;
s2: carrying out generalization treatment on the hydrogeological conditions of the site, and primarily establishing a numerical simulation model of underground water organic pollution multiphase flow for describing the transport mechanism of the organic pollutants in the underground water;
s3: determining pollution source characteristics to be identified and aquifer parameters in a pollution source inversion identification problem; according to the value range of the variable to be identified in the model, a plurality of groups of samples are randomly sampled and are substituted into the multiphase flow numerical model established in S2 one by one to obtain model responses corresponding to each group of samples, and a training sample set and an inspection sample set which are formed by model input-model response sample pairs are formed;
s4: establishing a single substitution model of the multiphase flow numerical model by adopting different machine learning methods according to the input-output training sample set obtained in the step S3;
s5: determining the weight and the core parameter of each single surrogate model established in S4 according to the test sample set obtained in S3 by applying a group intelligent optimization algorithm, and establishing a group intelligent integrated learning surrogate model;
s6: establishing a nonlinear programming optimization model for inversion identification of pollution source characteristics and aquifer parameters, and embedding the group intelligent ensemble learning substitution model established in S5 into the optimization model as one of constraint conditions;
s7: rewriting the optimization model established in S6 based on homotopy theory to obtain a series of inversion homotopy optimization models;
s8: and (4) solving the homotopy optimization model series established by S7 in sequence by adopting a hyper-heuristic algorithm, wherein the solution of the last optimization model is the inversion identification result of the underground water pollution source.
Further, a UTCHEM program is used for constructing a groundwater organic pollution multiphase flow migration model in S2.
Further, the characteristics of the pollution source to be identified in S3 include: longitudinal coordinates of the pollution source, transverse coordinates of the pollution source, migration and conversion time length of the pollutants and leakage amount of the pollutants; the aquifer parameters to be identified include: porosity, permeability, longitudinal aqueous phase dispersion, transverse aqueous phase dispersion; the random sampling is realized by a Latin hypercube sampling method.
Further, the single surrogate model modeling method in S4 includes: kriging, support vector regression, wavelet kernel extreme learning machine; the method for modeling the Kriging comprises the following steps: the regression equation for Kriging (Kriging) can be expressed as:
wherein f (x) is [ < f >1(x),f2(x),K,fk(x)]TBeing a known basis function, β ═ β (β)1,β2,K,βk)TThe regression parameters corresponding to the basis functions are obtained by estimating training samples; z (x) is a local deviation term, the mean is 0, and the variance isIts covariance can be expressed as:
r (u, v) is a correlation function of n-dimensional vectors u and v:
in the formula, alphaiFor the parameter to be determined, uiAnd viThe ith element being u and v;
for a given sample P ═ P1,p2,K,pm]TAnd corresponding output response Q ═ Q1,q2,K,qm]TThe predicted output response y (x) for any input vector x can be written as:
where r is the correlation vector between x and sample P:
r(x)=[R(x,p1),R(x,p2),K,R(x,pm)]T (5)
f is the response column vector for sample P:
r is a correlation matrix among sample points of the sample P:
the Support Vector Regression (SVR) modeling method is as follows:
SVR maps input data to a high-dimensional space through a nonlinear mapping function to perform linear regression, and balances training fitting precision and prediction precision; for a given training input X ═ X1,x2,K,xm]TEach element represents an N-ary input:
xi=(xi,1,xi,2,K,xi,N) I ═ 1,2Km and output Y ═ Y (Y)1,y2,K,ym)T(ii) a The nonlinear regression equation can be expressed as:
f(x)=<w,Φ(x)>+b (9)
wherein w is (w)1,w2,K,wN) As weight vectors, b as fitting errors,<w,x>represents the inner product of w and Φ (x); phi (x) is a nonlinear mapping function which can map the input vector from the input space to the high-dimensional space; this problem can be expressed as an optimization problem as follows:
Wherein | | w | | non-conducting phosphor2Is the norm of w, the constant C and the variable ξiAndε is the acceptable deviation for penalty factors and relaxation variables;
constructing a Lagrange function:
wherein L is Lagrangian, ηi、αiAndis a lagrange multiplier of not less than 0; with respect to the optimal solution, it is preferred that,andare all 0:
combining the above properties with equation 11, the optimization problem in equation 10 can be rewritten as a dual form:
From equation 16, it can be seen that:
thus, a regression function can be obtained:
the most common use of kernel functions instead of inner products < Φ (x), Φ (x') >, is gaussian:
And finally obtaining:
the fitting error b can be calculated by using the Karush-Kuhn-Tucker (KKT) condition;
the Wavelet Kernel Extreme Learning Machine (WKELM) modeling method is as follows:
for training sample (x)j,tj) J is 1, K, N, output y of the extreme learning machinejCan be expressed as:
q(xj)=[p(ω1xj+b1),p(ω2xj+b2),K,p(ωLxj+bL)]T (23)
wherein p (-) is an excitation function, and the input node passes through a weight vector omegaiConnected to the ith hidden neuron, i ═ 1, K, L, biThreshold for the ith hidden neuron, p (ω)ixj+bi) As an output function of the ith hidden neuron, betaiA weight vector connecting the ith hidden neuron with the output neuron;
Qβ=Y (24)
if an extreme learning machine model with L hidden nodes is able to learn N training samples unbiased, then the following formula exists:
in the formula, tjRepresents a target value;
equation 26 can be abbreviated as:
Qβ=T (27)
least square solution 27
β=Q+T (28)
In the formula, Q+Moore-Penrose generalized inverse of Q; q+Calculated from the following equation
Q+=QT(QQT)-1 (29)
For training sample (x)j,tj) The original optimization problem of j ═ 1, K, N, KELM is expressed as
s.t.q(xj)T·β=tj-ξj (30)
In the formula, C represents a regularization parameter capable of balancing training errors and algorithm complexity, ξjRepresents an error;
the optimization problem can be converted into a lagrange dual form solution:
in the formula, thetaiRepresents a lagrange operator; the dual problem can be solved by applying the condition of Karush-Kuhn-Tucker (KKT):
a least squares solution for the weight vector β can be calculated from equation 32:
according to the kernel function theory, a wavelet implicit mapping function K (x) is constructedi,xj) Instead of a random mapping function q (x)j):
KELM=QQT (34)
KELM(i,j)=q(xi)T·q(xj)=K(xi,xj) (35)
The trained KELM output function is expressed as follows:
wherein T is ═ T1,K,tN]T;
The expression of the wavelet kernel function is:
in the formula of alphaw,γ1,w,γ2,wIs an adjustable parameter.
Further, the ensemble learning surrogate model in S5 is a weighted linear superposition of the single surrogate model in S4, and the output expression thereof is as follows:
in the formula (I), the compound is shown in the specification,in order to integrate the predictive values of the surrogate model,the predicted values, w, of the Kriging model, the support vector regression model and the wavelet kernel limit learning model respectively1,w2,w3The sum of the weights is 1;
a key link for constructing the ensemble learning substitution model is determination of an integrated weight value and core parameters of each single substitution model, and the higher the prediction precision of the substitution model is, the higher the weight of the substitution model is; the poorer the accuracy of the surrogate model, the smaller its weight should be; the following method is adopted as a determination method for integrating the weight and the model parameter;
firstly, establishing an optimization model, wherein decision variables of the optimization model are weight values in the ensemble learning substitution model and variable parameter values of each single substitution model, and the optimal solution of the optimization model enables the root mean square error between the output of the ensemble learning substitution model obtained by using a test sample and the output of a simulation model to be minimum;
in the formula, yi(xk,pi) Using parameter vector p for ith substitution modeliAnd the kth test sample input xkFor inputting the corresponding output response, i is 1,2, K, N, K is 1,2, K, M, yactual(xk) Inputting x for the k test samplekOutputting response by the corresponding simulation model, wherein M is the number of the test samples; solving the optimization model by using a group intelligent optimization algorithm to obtain integrated weight and eachAnd constructing a group intelligent integrated learning substitution model by using the optimal value of the variable parameter of the single substitution model.
Further, the nonlinear optimization model in S6 takes the sum of squares of absolute deviations of the simulated calculated concentration values and the measured concentration values of the monitoring well as the minimum as an objective function; using variables to be identified in the underground water pollution source inversion identification problem as decision variables, wherein the decision variables comprise: longitudinal coordinates of the pollution source, transverse coordinates of the pollution source, migration and conversion time of the pollutants, leakage amount of the pollutants, porosity, permeability, longitudinal aqueous phase dispersity and transverse aqueous phase dispersity; taking a swarm intelligence integrated learning intelligent substitution model as equality constraint of a pollutant migration and transformation rule, and taking a pollution source characteristic and an aquifer parameter value range as inequality constraint conditions; the expression is as follows:
Ckthe actual monitoring value of the pollutant concentration in the monitoring well of the kth, K is 1,2, K, n,calculating a corresponding simulation calculation value;the method comprises the steps of calculating a value vector for a simulation calculation value of pollutant concentration, namely the output of an intelligent substitution model of a swarm intelligent integrated learning machine, wherein s is a value vector of a relevant variable of a pollution source, and p is a value vector of an aquifer parameter; the constraints respectively represent the longitudinal coordinates (X) of the pollution sourcei) Transverse coordinate (Y) of the source of pollutioni) Duration of contaminant migration and conversion (T)i) Amount of leakage of contaminants (V)i) Porosity (theta), permeability (K), longitudinal oil phase dispersion (D)oil,l) Transverse oil phase dispersivity (D)oil,t) Upper and lower bounds of the variation range.
Further, the homotopy optimization model series construction in the step 7 firstly needs to introduce homotopy parameters λ by means of homotopy algorithm idea, and construct homotopy functions H (s, p, λ) so that when λ is 0, the solution of the equation set H (s, p, λ) is 0 is the assumed groundwater pollution source and aquifer parameters, and when λ is 1, the solution of the equation set H (s, p, λ) is 0 is the true groundwater pollution source and aquifer parameters to be solved; starting from arbitrarily assumed pollution source information and aquifer parameter values, sequentially solving optimization problems corresponding to a series of homotopic equations by path tracking, gradually approaching and finally obtaining real pollution source information and aquifer parameter values to be solved;
the homotopy function takes a linear homotopy as follows:
H(s,p,λ)=λF(s,p)+(1-λ)G(s,p) (40)
wherein F (s, p) ═ F (s, p) -Cobs,G(s,p)=f(s,p)-C0;
In the formula, s is a value vector of a pollution source characteristic related variable, and p is a value vector of an aquifer parameter; λ is homotopy parameter, and the value range is [0,1 ]](ii) a f () represents a multiphase flow simulation model or a surrogate model thereof; cobsActually monitoring concentration vectors for monitoring points; c0The method comprises the steps of representing calculation of concentration vectors of monitoring points obtained by substituting any assumed underground water pollution source characteristics into a simulation model;
the homotopy function H continuously depends on homotopy parameter lambda, and the value range of lambda is within the range of [0,1]Is divided into0=0<λ1<…<λ N1, a series of equations is obtained:
H(s,p,λi)=f(s,p)-(λi·Cobs+(1-λi)·C0)=0,i=1,2,K,N (41)
if λi+1-λiSufficiently small, the solution of the adjacent equation (s, p)i+1And (s, p)iAre very close; according to the homotopy algorithm idea, the optimization model in the S6 is rewritten to obtain an optimization model series corresponding to a homotopy equation set series of the underground water pollution source inversion identification problem, as shown in a formula (42);the calculated concentration of the pollutants in the kth monitoring well obtained by substituting the characteristics of the underground water pollution source and the aquifer parameter values which are randomly assumed into the simulation model is represented;
further, the meta-heuristic algorithm in S7 adopts random selection as a high-level strategy, that is, a low-level heuristic algorithm is randomly selected from a given set in each solving, where the set includes a genetic algorithm, a particle swarm optimization algorithm, and a simulated annealing algorithm;
and (3) solving the homotopy optimization model series in sequence by using a hyper-heuristic algorithm (HHA), wherein in the solving process, because the path tracking process is gradually evolved, the solutions of two adjacent optimization problems are also very close to each other:
in the formula (I), the compound is shown in the specification,optimizing the initial population of the model for the ith homotopy,(s)i,opt,pi,opt) The optimal solution of the ith homotopy optimization model obtained by applying a hyper-heuristic algorithm can be used as a basis for generating the initial population of the next optimization model:
in the formula, vs,low,vs,up,vp,low,vp,upRepresenting the upper limit and the lower limit of an adjusting displacement vector corresponding to the pollution source characteristic and the aquifer parameter in the hyperheuristic algorithm;
the solution of the ith optimization problem is used as a generation basis for solving the initial value of the (i + 1) th optimization problem, so that the optimization process is ensured to be gradually evolved, and the premature convergence problem caused by the fact that the initial value is far away from the optimal solution is avoided; and the solution of the last optimization model is the inversion identification result of the underground water pollution source.
The invention has the advantages that:
1. according to the invention, a group intelligent integrated learning substitution model of the underground water organic pollution multiphase flow numerical model is established, and the inversion calculation efficiency is improved.
2. The invention constructs a hyper-heuristic homotopy algorithm for inversion optimization solution and reduces the calculation error of inversion.
3. The method solves the problems of identification of the characteristics of the underground water organic pollution source and correction of key parameters in the multiphase flow migration model, thereby realizing accurate simulation and prediction of the space-time distribution of the underground water organic pollution and providing important precondition basic conditions for reasonable formulation of an underground water pollution remediation scheme, pollution risk assessment and pollution responsibility confirmation.
Drawings
The invention is described in further detail below with reference to the drawings and the detailed description.
Fig. 1 is a flow architecture of the superheuristic homotopy algorithm solution optimization model in S7 and S8.
FIG. 2 is a diagram showing the relative positions of the site generalization and water quality monitoring wells of the embodiment.
FIG. 3 is a fitting scatter plot of the inspection output and the numerical simulation model output of each of the single surrogate models and the ensemble learning intelligent surrogate model established at S4 and S5.
FIG. 4 is a graph of the convergence curves of different variable identification values in the process of solving the optimization of the hyper-heuristic homotopy algorithm.
FIG. 5 is a comparison graph of an actual pollution plume and an identified pollution plume for an example: (a) actual pollution plume; (b) and identifying pollution plumes.
Detailed Description
In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some, but not all, embodiments of the present invention. The components of embodiments of the present invention generally described and illustrated in the figures herein may be arranged and designed in a wide variety of different configurations.
Thus, the following detailed description of the embodiments of the present invention, presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
It should be noted that: like reference numbers and letters refer to like items in the following figures, and thus, once an item is defined in one figure, it need not be further defined and explained in subsequent figures.
In the description of the present invention, it should be noted that the terms "center", "upper", "lower", "left", "right", "vertical", "horizontal", "inner", "outer", etc. indicate orientations or positional relationships based on the orientations or positional relationships shown in the drawings or the orientations or positional relationships that the products of the present invention are conventionally placed in use, and are only used for convenience in describing the present invention and simplifying the description, but do not indicate or imply that the devices or elements referred to must have a specific orientation, be constructed and operated in a specific orientation, and thus, should not be construed as limiting the present invention. Furthermore, the terms "first," "second," "third," and the like are used solely to distinguish one from another and are not to be construed as indicating or implying relative importance.
Furthermore, the terms "horizontal", "vertical" and the like do not imply that the components are absolutely horizontal or hanging, but may be slightly inclined. For example, "horizontal" merely means that the direction is more horizontal than "vertical" and does not mean that the structure must be perfectly horizontal, but may be slightly inclined.
In the description of the present invention, it should also be noted that, unless otherwise explicitly specified or limited, the terms "disposed," "mounted," "connected," and "connected" are to be construed broadly and may, for example, be fixedly connected, detachably connected, or integrally connected; can be mechanically or electrically connected; they may be connected directly or indirectly through intervening media, or they may be interconnected between two elements. The specific meanings of the above terms in the present invention can be understood in specific cases to those skilled in the art.
Example 1: an organically polluted diving aquifer which can be generalized into a homogeneous and isotropic three-dimensional multiphase flow model; no natural boundary is arranged near the polluted site, and the boundary is defined at a position where the influence of the migration of pollutants can be ignored; wherein, the northeast boundary and the southwest boundary are generalized into a type of boundary; the southeast boundary and the northwest boundary are composed of flow surfaces and generalized to zero-flux boundaries; the lower part of the calculation simulation area is a water-resisting layer which can be generalized to a zero-flux boundary, the upper part of the calculation simulation area is a diving surface which is a water exchange boundary, and the thickness of the water-bearing layer is slowly changed along the direction of underground water flow, so that the calculation simulation area is generalized to an equal-thickness water-bearing layer; the physicochemical parameters of water and chlorobenzene, organic contaminants, are detailed in table 1;
TABLE 1 physicochemical parameters of Water and organic contaminants
The organic pollution in the field is caused by a punctiform pollution source, the pollutants enter an aquifer in a short time, and most of the following time is a natural dissolution and diffusion stage of the pollutants in the aquifer; therefore, the variables to be identified in the simulation model (the variables to be identified in the pollution source inversion identification problem) can be finally determined as: 1. the lateral, longitudinal coordinates of the source of contamination (contamination source location); the pollution source is located at the top of the aquifer by default, and the vertical coordinate does not need to be identified; 2. the length of time for the transfer of the contaminant (simulation period); 3. the leakage amount of the pollutants (the migration and conversion time of the pollutants and the leakage amount of the pollutants are the pollution source release history in the pollution source inversion identification problem); 4. aquifer parameters (including porosity, permeability, longitudinal aqueous phase dispersion, transverse aqueous phase dispersion); acquiring pollutant concentration monitoring data by using 5 water level water quality monitoring wells in a field as known information for inversion solution; the site generalization and monitoring well location are shown in figure 2.
The actual values of the pollution source information and aquifer parameters given in the example are shown in table 2;
table 2 pollution source information and aquifer parameters in the examples
Based on the information, using the UTCHEM to construct a model, establishing a simulation model of the groundwater DNAPLS polluted multiphase flow related to the embodiment, and using the established simulation model to perform forward forecasting calculation to obtain the space-time distribution condition of the pollutant concentration in the seepage field, namely the pollutant concentration at the monitoring well is used as actual monitoring data, and the pollutant concentration at the bottom of the aquifer at the five monitoring wells at the end is shown in a table 3;
TABLE 3 example forward model calculation output (actual monitoring data)
By using the method, the characteristics of each pollution source and the parameters of the aquifer are subjected to inversion calculation according to the concentration monitoring data in the table 3;
in this embodiment:
the prior intervals of all parameters in the step 3 are shown in a table 4;
TABLE 4 Parametric prior distribution characterization
The number of training sample samples in S3 is 100, and 20 test samples are additionally adopted;
s4, the fitting condition of the inspection output and the numerical model output of each single surrogate model and the integrated learning surrogate model established in the step 5 is shown in figure 3;
when the optimal model is solved by applying the hyper-heuristic homotopy algorithm in S7, the interval of homotopy parameters is 0.1, an equation series containing 10 homotopy equations is constructed, and the optimal model is rewritten according to the homotopy equation series, so that corresponding 10 optimal models can be obtained.
When solving the homotopy optimization model series in the S8, the maximum evolutionary algebra solved by the first nine optimization models is 60, and the maximum evolutionary algebra solved by the last optimization model is 120; the inversion identification result and the relative error are obtained and shown in the table 5;
TABLE 5 identification of results and relative errors
According to the results in table 5, the inversion values are very close to the true values, and the errors are less than 5%;
the progressive optimization ability of the hyper-heuristic-homotopy algorithm can be clearly shown by the convergence curves identified by different pollution source characteristics and aquifer parameters, see fig. 4; in the course of iterative optimization, the identification values of the different variables, which are accompanied by a change in the homotopy parameters, gradually approach the true values, in particular the longitudinal and transverse coordinates of the contamination source, and the contamination leakage. The initial large-scale search problem is converted into a plurality of local search stages, and the optimal value of each stage is easy to solve; although an undesirable solution phase that falls short of premature convergence may occur, a homotopy algorithm-based progressive solution mechanism may reduce and disperse the risk, and subsequent solution phases may again bring the convergence trend toward a global optimum.
Substituting the identification value into the multiphase flow numerical model established in S2, calculating the distribution condition of the pollutants in the aquifer at the end moment, and comparing the distribution condition with the actual distribution of the pollutants at the end moment of the embodiment, as shown in FIG. 5; the shapes of the pollution plumes are very close, and the current or forecast distribution condition of pollutants in a future aquifer can be accurately calculated by applying a simulation model established by the recognition result, so that a reliable basis is provided for design of a groundwater pollution remediation scheme and risk assessment.
If the pollution source inversion identification problem in the embodiment is solved by using a traditional simulation-optimization method, a multiphase flow simulation model 39600 (the number of initial populations of the hyper-heuristic method is 60, and the maximum evolutionary algebra is 660 times) needs to be called; the polluted aquifer multiphase flow simulation model in the example is operated on a computer with the CPU of Intel core i53.0 GHz and the memory of 8GB, the operation time is about 600 seconds on average, and the whole optimization solving process takes about 275 days. And the integrated learning substitution model only needs 1.2 seconds to operate once, and if the integrated learning substitution model is used for substituting for the simulation model, the optimization solving process is shortened to 13.2 hours.
Therefore, the method can effectively solve the problems of parameter identification and model correction, not only improves the inversion efficiency, but also reduces the inversion error.
It will be understood by those skilled in the art that the present invention is not limited to the embodiments described above, which are described in the specification and illustrated only to illustrate the principle of the present invention, but that various changes and modifications may be made therein without departing from the spirit and scope of the present invention, which fall within the scope of the invention as claimed. The scope of the invention is defined by the appended claims and equivalents thereof.
Claims (8)
1. An underground water organic pollution source inversion method based on a heuristic homotopy algorithm is characterized by comprising the following steps: the method comprises the following steps:
s1: mastering hydrogeological conditions of a polluted site by means of field on-site investigation and dynamic monitoring to obtain groundwater water quality dynamic monitoring data;
s2: carrying out generalization treatment on site hydrogeological conditions, and primarily establishing a numerical simulation model of underground water organic pollution multiphase flow for describing a transport mechanism of organic pollutants in underground water;
s3: determining pollution source characteristics to be identified and aquifer parameters in a pollution source inversion identification problem; according to the value range of the variable to be identified in the model, a plurality of groups of samples are randomly sampled and are substituted into the multiphase flow numerical model established in S2 one by one to obtain model responses corresponding to each group of samples, and a training sample set and an inspection sample set which are formed by model input-model response sample pairs are formed;
s4: according to the input-output training sample set obtained in the S3, establishing a single substitution model of the multiphase flow numerical model by adopting different machine learning methods;
s5: determining the weight and the core parameter of each single surrogate model established in S4 according to the test sample set obtained in S3 by applying a group intelligent optimization algorithm, and establishing a group intelligent integrated learning surrogate model;
s6: establishing a nonlinear programming optimization model for inversion identification of pollution source characteristics and aquifer parameters, and embedding the group intelligent integrated learning substitution model established in S5 into the optimization model as one of constraint conditions;
s7: rewriting the optimization model established in S6 based on homotopy theory to obtain a series of inversion homotopy optimization models;
s8: and (4) solving the homotopy optimization model series established by S7 in sequence by adopting a hyper-heuristic algorithm, wherein the solution of the last optimization model is the inversion identification result of the underground water pollution source.
2. The groundwater organic pollution source inversion method based on the heuristic-homotopy algorithm of claim 1, wherein: and in the S2, a UTCHEM program is used for constructing a groundwater organic pollution multiphase flow migration model.
3. The groundwater organic pollution source inversion method based on the heuristic-homotopy algorithm of claim 1, wherein: the pollution source characteristics to be identified in S3 include: longitudinal coordinates of the pollution source, transverse coordinates of the pollution source, migration and conversion time length of the pollutants and leakage amount of the pollutants; the aquifer parameters to be identified include: porosity, permeability, longitudinal aqueous phase dispersion, transverse aqueous phase dispersion; the random sampling is realized by a Latin hypercube sampling method.
4. The groundwater organic pollution source inversion method based on the heuristic-homotopy algorithm as claimed in claim 1, characterized in that: the single surrogate model modeling method in S4 includes: kriging, support vector regression, wavelet kernel extreme learning machine; the method for modeling the Kriging comprises the following steps: the regression equation for Kriging (Kriging) can be expressed as:
wherein f (x) is [ < f >1(x),f2(x),K,fk(x)]TIs a known basis function, β ═ β (β)1,β2,K,βk)TThe regression parameters corresponding to the basis functions are obtained by estimating training samples; z (x) is a local deviation term, the mean is 0, and the variance isIts covariance can be expressed as:
r (u, v) is a correlation function of n-dimensional vectors u and v:
in the formula, alphaiFor the parameter to be determined, uiAnd viThe ith element being u and v;
for a given sample P ═ P1,p2,K,pm]TAnd corresponding output response Q ═ Q1,q2,K,qm]TThe predicted output response y (x) for any input vector x can be written as:
where r is the correlation vector between x and sample P:
r(x)=[R(x,p1),R(x,p2),K,R(x,pm)]T (5)
f is the response column vector for sample P:
r is a correlation matrix among sample points of the sample P:
the Support Vector Regression (SVR) modeling method is as follows:
SVR is a linear regression that maps input data to a high-dimensional space by a non-linear mapping function whereThe fitting precision and the prediction precision are trained and balanced; for a given training input X ═ X1,x2,K,xm]TEach element represents an N-ary input:
xi=(xi,1,xi,2,K,xi,N) I ═ 1,2Km and output Y ═ Y (Y)1,y2,K,ym)T(ii) a The nonlinear regression equation can be expressed as:
f(x)=<w,Φ(x)>+b (9)
wherein w ═ w1,w2,K,wN) As weight vectors, b as fitting errors,<w,x>represents the inner product of w and Φ (x); phi (x) is a nonlinear mapping function which can map the input vector from the input space to the high-dimensional space; this problem can be expressed as an optimization problem as follows:
Wherein | | w | | non-conducting phosphor2Is the norm of w, the constant C and the variable ξiAndε is the acceptable deviation for penalty factors and relaxation variables;
constructing a Lagrange function:
wherein L is Lagrangian, ηi、αiAndis a lagrange multiplier of not less than 0; with respect to the optimal solution, it is preferred that,andare all 0:
combining the above properties with equation 11, the optimization problem in equation 10 can be rewritten as a dual form:
From equation 16, it can be seen that:
thus, a regression function can be obtained:
the most common use of kernel functions instead of inner products < Φ (x), Φ (x') >, is gaussian:
formula 16 can be rewritten as:
And finally obtaining:
the fitting error b can be calculated by using the Karush-Kuhn-Tucker (KKT) condition;
the Wavelet Kernel Extreme Learning Machine (WKELM) modeling method is as follows:
for training sample (x)j,tj) J is 1, K, N, output y of the extreme learning machinejCan be expressed as:
q(xj)=[p(ω1xj+b1),p(ω2xj+b2),K,p(ωLxj+bL)]T (23)
wherein p (-) is an excitation function, and the input node passes through a weight vector omegaiConnected to the ith hidden neuron, i ═ 1, K, L, biThreshold for the ith hidden neuron, p (ω)ixj+bi) As an output function of the ith hidden neuron, betaiA weight vector connecting the ith hidden neuron with the output neuron;
equation 22 can be expressed as follows:
Qβ=Y (24)
wherein β ═ β1,K,βL]T,Y=[y1,K,yN]TQ is the hidden layer output matrix of ELM:
if an extreme learning machine model with L hidden nodes is able to learn N training samples unbiased, then the following formula exists:
in the formula, tjRepresents a target value;
equation 26 can be abbreviated as:
Qβ=T (27)
least square solution 27
β=Q+T (28)
In the formula, Q+Moore-Penrose generalized inverse of Q; q+Is calculated by the following formulaCalculating out
Q+=QT(QQT)-1 (29)
For training sample (x)j,tj) The original optimization problem of j ═ 1, K, N, KELM is expressed as
s.t.q(xj)T·β=tj-ξj (30)
In the formula, C represents a regularization parameter capable of balancing training errors and algorithm complexity, ξjRepresents an error;
the optimization problem can be converted into a lagrange dual form solution:
in the formula, thetaiRepresents a lagrange operator; the dual problem can be solved by applying the condition of Karush-Kuhn-Tucker (KKT):
a least squares solution for the weight vector β can be calculated from equation 32:
according to the kernel function theory, a wavelet implicit mapping function K (x) is constructedi,xj) Instead of a random mapping function q (x)j):
KELM=QQT (34)
KELM(i,j)=q(xi)T·q(xj)=K(xi,xj) (35)
The trained KELM output function is expressed as follows:
wherein T is ═ T1,K,tN]T;
The expression of the wavelet kernel function is:
in the formula of alphaw,γ1,w,γ2,wIs an adjustable parameter.
5. The underground water organic pollution source inversion method based on the meta-heuristic-homotopy algorithm as claimed in claim 1 is characterized in that: the ensemble learning surrogate model in S5 is formed by weighted linear superposition of the single surrogate model in S4, and the expression of the output is as follows:
in the formula (I), the compound is shown in the specification,in order to integrate the predictive values of the surrogate model,the predicted values, w, of the Kriging model, the support vector regression model and the wavelet kernel limit learning model respectively1,w2,w3The sum of the weights is 1;
a key link for constructing the ensemble learning substitution model is determination of an integrated weight value and core parameters of each single substitution model, and the higher the prediction precision of the substitution model is, the higher the weight of the substitution model is; the poorer the accuracy of the surrogate model, the smaller its weight should be; the following method is adopted as a determination method for integrating the weight and the model parameter;
firstly, establishing an optimization model, wherein decision variables of the optimization model are weight values in the ensemble learning substitution model and variable parameter values of each single substitution model, and the optimal solution of the optimization model enables the root mean square error between the output of the ensemble learning substitution model obtained by using a test sample and the output of a simulation model to be minimum;
in the formula, yi(xk,pi) Using parameter vector p for ith substitution modeliAnd the kth test sample input xkFor inputting the corresponding output response, i is 1,2, K, N, K is 1,2, K, M, yactual(xk) Inputting x for the k test samplekOutputting response by the corresponding simulation model, wherein M is the number of the test samples; and solving the optimization model by applying a group intelligent optimization algorithm to obtain the integrated weight and the optimal value of the variable parameter of each single substitution model, and constructing a group intelligent integrated learning substitution model.
6. The underground water organic pollution source inversion method based on the meta-heuristic-homotopy algorithm as claimed in claim 1 is characterized in that: the nonlinear optimization model in the S6 takes the minimum sum of squares of absolute deviations of the simulated calculation concentration values and the measured concentration values of the monitoring well as an objective function; the method takes variables to be identified in the underground water pollution source inversion identification problem as decision variables, and comprises the following steps: longitudinal coordinates of the pollution source, transverse coordinates of the pollution source, migration and conversion time of the pollutants, leakage amount of the pollutants, porosity, permeability, longitudinal aqueous phase dispersity and transverse aqueous phase dispersity; taking a swarm intelligence integrated learning intelligent substitution model as equality constraint of a pollutant migration and transformation rule, and taking a pollution source characteristic and an aquifer parameter value range as inequality constraint conditions; the expression is as follows:
Ckthe actual monitoring value of the pollutant concentration in the monitoring well of the kth, K is 1,2, K, n,calculating a corresponding simulation calculation value;the method comprises the steps of calculating a value vector for a simulation calculation value of pollutant concentration, namely the output of an intelligent substitution model of a swarm intelligent integrated learning machine, wherein s is a value vector of a relevant variable of a pollution source, and p is a value vector of an aquifer parameter; the constraints respectively represent the longitudinal coordinates (X) of the pollution sourcei) Transverse coordinate (Y) of the source of pollutioni) Duration of pollutant migration conversion (T)i) Amount of leakage of contaminants (V)i) Porosity (theta), permeability (K), longitudinal oil phase dispersion (D)oil,l) Transverse oil phase dispersivity (D)oil,t) Upper and lower bounds of the variation range.
7. The groundwater organic pollution source inversion method based on the heuristic-homotopy algorithm as claimed in claim 1, characterized in that: the homotopy optimization model series in step 7 is constructed by introducing homotopy parameters λ and constructing homotopy functions H (s, p, λ) by means of homotopy algorithm thought, so that when λ is 0, the solution of the equation set H (s, p, λ) is 0 is the assumed groundwater pollution source and aquifer parameters, and when λ is 1, the solution of the equation set H (s, p, λ) is 0 is the real groundwater pollution source and aquifer parameters to be solved; starting from arbitrarily assumed pollution source information and aquifer parameter values, sequentially solving optimization problems corresponding to a series of homotopic equations by path tracking, gradually approaching and finally obtaining real pollution source information and aquifer parameter values to be solved;
the homotopy function takes a linear homotopy as follows:
H(s,p,λ)=λF(s,p)+(1-λ)G(s,p) (40)
wherein F (s, p) ═ F (s, p) -Cobs,G(s,p)=f(s,p)-C0;
In the formula, s is a value vector of a pollution source characteristic related variable, and p is a value vector of an aquifer parameter; λ is homotopy parameter, and the value range is [0,1 ]](ii) a f () represents a multiphase flow simulation model or a surrogate model thereof; cobsActually monitoring concentration vectors for monitoring points; c0The method comprises the steps of representing calculation of concentration vectors of monitoring points obtained by substituting any assumed underground water pollution source characteristics into a simulation model;
the homotopy function H continuously depends on homotopy parameter lambda, and the value range of lambda is within the range of [0,1]Is divided into0=0<λ1<…<λN1, a series of equations is obtained:
H(s,p,λi)=f(s,p)-(λi·Cobs+(1-λi)·C0)=0,i=1,2,K,N (41)
if λi+1-λiSufficiently small, the solution of the adjacent equation (s, p)i+1And (s, p)iAre very close;
the optimization model in S6 is rewritten according to homotopy algorithm ideaObtaining an optimization model series corresponding to the homotopic equation set series of the underground water pollution source inversion identification problem, as shown in a formula (42);the calculated concentration of the pollutants in the kth monitoring well obtained by substituting the characteristics of the underground water pollution source and the aquifer parameter values which are randomly assumed into the simulation model is represented;
8. the underground water organic pollution source inversion method based on the meta-heuristic-homotopy algorithm as claimed in claim 1 is characterized in that: the hyper-heuristic algorithm in the S7 adopts random selection as a high-level strategy, namely, a low-level heuristic algorithm is randomly selected from a given set in each solving, and the set comprises a genetic algorithm, a particle swarm optimization algorithm and a simulated annealing algorithm;
and (3) solving the homotopy optimization model series in sequence by using a hyper-heuristic algorithm (HHA), wherein in the solving process, because the path tracking process is gradually evolved, the solutions of two adjacent optimization problems are also very close to each other:
in the formula (I), the compound is shown in the specification,for the initial population of the ith homotopy optimization model,(s)i,opt,pi,opt) For ith homotopy mode of optimization derived by applying hyper-heuristic algorithmThe optimal solution of the model can be used as a basis for generating the initial population of the next optimization model:
in the formula, vs,low,vs,up,vp,low,vp,upRepresenting the upper limit and the lower limit of an adjustment displacement vector corresponding to the pollution source characteristic and the aquifer parameter in the hyperheuristic algorithm;
the solution of the ith optimization problem is used as a generation basis for solving the initial value of the (i + 1) th optimization problem, so that the optimization process is ensured to be gradually evolved, and the premature convergence problem caused by the fact that the initial value is far from the optimal solution is avoided; and the solution of the last optimization model is the inversion identification result of the underground water pollution source.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111602792.6A CN114707395A (en) | 2021-12-24 | 2021-12-24 | Groundwater organic pollution source inversion method based on heuristic homotopy algorithm |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111602792.6A CN114707395A (en) | 2021-12-24 | 2021-12-24 | Groundwater organic pollution source inversion method based on heuristic homotopy algorithm |
Publications (1)
Publication Number | Publication Date |
---|---|
CN114707395A true CN114707395A (en) | 2022-07-05 |
Family
ID=82167223
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202111602792.6A Pending CN114707395A (en) | 2021-12-24 | 2021-12-24 | Groundwater organic pollution source inversion method based on heuristic homotopy algorithm |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114707395A (en) |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115330132A (en) * | 2022-07-20 | 2022-11-11 | 中交上海航道局有限公司 | Method for water quality distribution reverse-time inversion of wide and shallow river in sudden pollution accident |
CN115329607A (en) * | 2022-10-14 | 2022-11-11 | 山东省鲁南地质工程勘察院(山东省地质矿产勘查开发局第二地质大队) | System and method for evaluating underground water pollution |
CN115510977A (en) * | 2022-09-28 | 2022-12-23 | 吉林建筑大学 | Geological statistical pattern recognition method based on Bayes ensemble learning machine |
CN115935626A (en) * | 2022-11-25 | 2023-04-07 | 河南大学 | Inversion method of river water-underground water vertical transient interaction water flux based on VAE-PINN algorithm |
CN116312824A (en) * | 2023-03-20 | 2023-06-23 | 北京市生态环境保护科学研究院 | Water-bearing layer heterogeneity recognition method aiming at non-Fick dispersion of groundwater pollutants |
-
2021
- 2021-12-24 CN CN202111602792.6A patent/CN114707395A/en active Pending
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115330132A (en) * | 2022-07-20 | 2022-11-11 | 中交上海航道局有限公司 | Method for water quality distribution reverse-time inversion of wide and shallow river in sudden pollution accident |
CN115510977A (en) * | 2022-09-28 | 2022-12-23 | 吉林建筑大学 | Geological statistical pattern recognition method based on Bayes ensemble learning machine |
CN115329607A (en) * | 2022-10-14 | 2022-11-11 | 山东省鲁南地质工程勘察院(山东省地质矿产勘查开发局第二地质大队) | System and method for evaluating underground water pollution |
CN115329607B (en) * | 2022-10-14 | 2023-02-03 | 山东省鲁南地质工程勘察院(山东省地质矿产勘查开发局第二地质大队) | System and method for evaluating underground water pollution |
CN115935626A (en) * | 2022-11-25 | 2023-04-07 | 河南大学 | Inversion method of river water-underground water vertical transient interaction water flux based on VAE-PINN algorithm |
CN115935626B (en) * | 2022-11-25 | 2023-09-08 | 河南大学 | Inversion method of river water-underground water vertical transient interaction water flow |
CN116312824A (en) * | 2023-03-20 | 2023-06-23 | 北京市生态环境保护科学研究院 | Water-bearing layer heterogeneity recognition method aiming at non-Fick dispersion of groundwater pollutants |
CN116312824B (en) * | 2023-03-20 | 2023-08-22 | 北京市生态环境保护科学研究院 | Water-bearing layer heterogeneity recognition method aiming at non-Fick dispersion of groundwater pollutants |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN114707395A (en) | Groundwater organic pollution source inversion method based on heuristic homotopy algorithm | |
Hai et al. | Global solar radiation estimation and climatic variability analysis using extreme learning machine based predictive model | |
Shiri et al. | Comparison of genetic programming with neuro-fuzzy systems for predicting short-term water table depth fluctuations | |
Shiri et al. | Predicting groundwater level fluctuations with meteorological effect implications—A comparative study among soft computing techniques | |
CN105243435B (en) | A kind of soil moisture content prediction technique based on deep learning cellular Automation Model | |
Liu et al. | Application of artificial neural networks in global climate change and ecological research: An overview | |
CN112905560A (en) | Air pollution prediction method based on multi-source time-space big data deep fusion | |
CN114492164A (en) | Organic pollutant migration numerical model substitution method based on multi-core extreme learning machine | |
Phillips | Stability implications of the state factor model of soils as a nonlinear dynamical system | |
CN109001136A (en) | A kind of COD on-line monitoring method based on ultraviolet visible light absorption spectrum | |
CN108843312B (en) | Oilfield reservoir in-layer heterogeneity integrated evaluating method | |
Kişi | Evolutionary fuzzy models for river suspended sediment concentration estimation | |
Duka et al. | Processing, neural network-based modeling of biomonitoring studies data and validation on Republic of Moldova data | |
CN115952685B (en) | Sewage treatment process soft measurement modeling method based on integrated deep learning | |
CN112668606B (en) | Step type landslide displacement prediction method based on gradient elevator and quadratic programming | |
CN114399073A (en) | Ocean surface temperature field prediction method based on deep learning | |
CN115906626A (en) | Underground water organic pollution mutual-feedback inversion tracing method based on artificial intelligence | |
Xiaoyi et al. | Research on hybrid mechanism modeling of algal bloom formation in urban lakes and reservoirs | |
MacBean et al. | Quantifying and reducing uncertainty in global carbon cycle predictions: Lessons and perspectives from 15 years of data assimilation studies with the ORCHIDEE terrestrial biosphere model | |
CN116383589A (en) | Spatial interpolation prediction method for soil nutrients | |
Jiang et al. | Simultaneous identification of contaminant sources and hydraulic conductivity field by combining geostatistics method with self-organizing maps algorithm | |
Bandara et al. | Towards soil property retrieval from space: Proof of concept using in situ observations | |
Saeed et al. | Water temperature prediction in a subtropical subalpine lake using soft computing techniques | |
Bian et al. | Bayesian ensemble machine learning-assisted deterministic and stochastic groundwater DNAPL source inversion with a homotopy-based progressive search mechanism | |
CN114034334A (en) | Method for identifying pollution source and flow of karst pipeline |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination |