A method for identifying non-Gaussian aquifer parameters by fusing groundwater level and natural potential data based on convolutional neural network

Abstract

The present invention discloses a method for identifying non-Gaussian aquifer parameters based on a convolutional neural network and fusing groundwater level and natural potential data, comprising: 1: training a CVAE network using non-Gaussian aquifer parameter field samples, sampling from a standard normal distribution, and initializing an estimated set of latent vectors z; 2: inputting the estimated set of latent vectors z into a trained CVAE decoder to reconstruct a corresponding estimated set of non-Gaussian aquifer parameter fields; 3: during a harmonic pumping test, running a hydrogeophysical forward model based on the reconstructed estimated set of non-Gaussian aquifer parameter fields to obtain hydraulic head and natural potential simulation data; 4: combining hydraulic head and natural potential observation data, and iteratively updating the estimated set of latent vectors z using an ESMDA method; repeating steps 2 to 4 until the maximum number of iterations is reached; 5: for the posterior set of updated latent vectors z, reconstructing the estimated result of the non-Gaussian aquifer parameter field through a CVAE decoder.

Description

Technical Field

[0001] The present invention relates to the technical field of groundwater numerical simulation, and in particular to a method for identifying non-Gaussian aquifer parameters by fusing groundwater level and natural potential data based on a convolutional neural network. Background Art

[0002] Numerical simulation methods have been widely used in groundwater management and remediation of contaminated aquifers. Reliable groundwater models require accurate knowledge of the properties of underground aquifers, such as permeability. However, due to the limitations of the number and location of boreholes, it is usually difficult to obtain comprehensive information about aquifers by directly measuring aquifer-related parameters through boreholes. Therefore, the current main approach is to characterize aquifer structure through stochastic inversion methods based on indirect observation data (such as hydraulic head and solute concentration).

[0003] Ensemble-based data assimilation methods are a widely promoted class of stochastic inversion methods, such as Ensemble Smoother with Multiple Data Assimilation (ESMDA, Emerick and Reynolds, 2013), which have been applied in many fields (such as hydrogeology, petroleum engineering, and oceanography). However, data assimilation methods require that the unknown parameters obey a Gaussian distribution. This condition cannot be met in many practical problems. For example, in a channel-shaped aquifer with fluvial sediments, there are usually several highly contrasting facies coexisting, so that the permeability presents a non-Gaussian pattern (Zhou et al., 2014). Since the spatial distribution of non-Gaussian permeability is significantly different from the assumed Gaussian distribution, directly using traditional data assimilation methods to estimate the non-Gaussian permeability field cannot obtain the true distribution of aquifer permeability. Therefore, many methods have been developed to overcome this non-Gaussian problem. In recent years, deep learning methods have begun to be applied to subsurface inversion modeling. One strategy is to apply deep learning methods to parameterize the non-Gaussian permeability field into a latent variable that obeys a Gaussian distribution, thereby satisfying the Gaussian assumption of the data assimilation method. For example, Laloy et al. (2018) used a spatial generative adversarial neural (SGAN) network to parameterize complex non-Gaussian geological structures; Canchumuni et al. (2019) studied the performance of the convolutional variational autoencoder (CVAE) in parameterizing non-Gaussian lithofacies distribution. Kang et al. (2021a) coupled the CVAE network to ESMDA to invert the structure of underground pollution source areas with non-Gaussian distributions.

[0004] Although the above deep learning-based methods can effectively handle non-Gaussian problems, related studies only rely on hydrogeological observation data (such as hydraulic head and solute concentration) in invasive boreholes. Since obtaining these observation data is time-consuming, costly and labor-intensive, the amount of observation data available for assimilation is very limited. Under the condition of limited data, the estimated results of the non-Gaussian permeability field are difficult to characterize the fine structure of the true permeability field. Summary of the invention

[0005] Purpose of the invention: In order to overcome the problems of insufficient borehole observation data and limited estimation accuracy in the inversion of traditional non-Gaussian permeability coefficient fields, the present invention proposes a method for identifying non-Gaussian aquifer parameters by fusing groundwater level and natural potential data based on a convolutional neural network. By introducing the natural potential method in geophysical methods and fusing low-cost natural potential data, the natural potential data is used as a supplement under the premise that hydrogeological observation data are limited, so as to achieve fine identification of non-Gaussian permeability coefficient fields.

[0006] Technical solution: A method for identifying non-Gaussian aquifer parameters by fusing groundwater level and natural potential data based on convolutional neural network, comprising the following steps:

[0007] Step 1: Train the CVAE network using non-Gaussian aquifer parameter field samples and initialize the latent vector z estimation set by sampling from the standard normal distribution;

[0008] Step 2: Input the estimated set of latent vectors z into the trained CVAE decoder to reconstruct the corresponding estimated set of non-Gaussian aquifer parameter fields;

[0009] Step 3: During the harmonic pumping test, the hydrogeophysical forward model is run based on the reconstructed non-Gaussian aquifer parameter field estimation set to obtain hydraulic head simulation data and natural potential simulation data;

[0010] Step 4: Based on the hydraulic head simulation data and natural potential simulation data obtained in step 3, combined with the hydraulic head observation data and the natural potential observation data, the ESMDA method is used to iteratively update the latent vector z estimation set; repeat steps 2 to 4 until the maximum number of iterations is reached;

[0011] Step 5: For the updated latent vector z estimation set, reconstruct the estimation result of the non-Gaussian aquifer parameter field through the CVAE decoder;

[0012] The hydrogeophysical forward model is composed of a natural potential forward model coupled with a groundwater flow forward model through the Darcy velocity q and the intrinsic permeability k.

[0013] Furthermore, the non-Gaussian aquifer parameter fields to be estimated include permeability field, water storage rate field, porosity field and conductivity field.

[0014] Furthermore, the groundwater flow forward model is expressed as:

[0015]

[0016] Among them, S s is the water storage rate; H is the hydraulic head; t is the time; represents the divergence operator; represents the gradient operator; K is the permeability coefficient; W represents the source and sink per unit volume;

[0017] The groundwater flow forward model is subject to the following initial and boundary conditions:

[0018] H=H 0 at t=0 (2)

[0019] H=H D at Γ D (3)

[0020]

[0021] In the formula, H 0 is the initial water head; Γ D is the Dirichlet boundary; H D is Γ D Given water head on Γ N is the Neumann boundary; n is perpendicular to Γ N The exterior unit vector of .

[0022] Furthermore, the spontaneous potential forward model is expressed as:

[0023]

[0024] The boundary conditions are as follows:

[0025]

[0026]

[0027] Where σ is the conductivity of the porous medium, represents the electric potential, Γ D is the Dirichlet boundary, Γ N is the Neumann boundary; q is the Darcy velocity, Represents the effective charge density generated by groundwater flow.

[0028] Furthermore, the operation of the hydrogeophysical forward model includes first solving the groundwater flow forward model and then solving the natural potential forward model.

[0029] Furthermore, the CVAE network includes a CVAE encoder and a CVAE decoder; the encoder consists of six convolutional layers; and the CVAE decoder is a mirror-symmetric network of the encoder with a deconvolution layer.

[0030] Beneficial effects: Compared with the prior art, the present invention has the following advantages:

[0031] (1) The present invention overcomes the defect in the prior art that it is difficult to accurately describe the fine structure of the permeability coefficient of a complex non-Gaussian aquifer due to limited head observation data. By introducing the natural potential method, while conducting a pumping test, the natural potential change at the position of the non-polarized electrode caused by the groundwater flow is observed, thereby obtaining a large amount of geophysical data, which complements the traditional hydraulic head data and constrains the inversion of the non-Gaussian permeability coefficient field.

[0032] (2) Compared with the traditional hydraulic tomography method, the coupled hydrogeophysical inversion framework of the present invention couples the SP method to the inversion framework of the non-Gaussian permeability field, which can realize the fine identification of the non-Gaussian permeability field under the premise of limited head observation data, while ensuring the characterization accuracy, reducing the number of drilling holes and reducing the observation cost. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] Figure 1 is a flow chart of the present invention;

[0034] Figure 2 It is a non-Gaussian confined aquifer; (a) is a three-dimensional view of the hypothetical study area, the light and dark quadrilaterals represent the hydraulic tomography injection holes located at z = 3m and z = 6m, respectively, and the dark surface at x = 0m in the study area represents the tracer source in the solute migration scenario; (b) is a plan view of the aquifer;

[0035] Figure 3 Schematic diagram of aquifer; (a) original training image used to generate CVAE training samples; (b) and (c) distribution of lithofacies reference field; the area labeled 1 in the image represents the high permeability preferential flow channel, and the area labeled 0 represents the low permeability background lithofacies;

[0036] Figure 4Three lithofacies images randomly selected from the test set, the corresponding CVAE reconstructed images and histograms of latent variables, and three lithofacies images randomly generated by the CVAE decoder; where (a) are the three original input images; (b) are the corresponding three CVAE reconstructed images; (c) is the corresponding latent variable distribution; (d) are the three lithofacies images randomly generated by the trained CVAE network;

[0037] Figure 5 is the collective mean field, an individual field (selected randomly) and the collective variance field of ln K;

[0038] Figure 6 is the residual field of ln K;

[0039] Figure 7 is the water head cross-validation result, where (a) is the initial set; (b) is the validation result of Example 1; (c) is the validation result of Example 2; (d) is the validation result of Example 3;

[0040] Figure 8 : The mean field of the tracer concentration distribution at z = 5 m for the initial set and after 100, 150 and 200 days in the three exemplary cases. DETAILED DESCRIPTION

[0041] The technical solution of the present invention is now further described in conjunction with the accompanying drawings.

[0042] like Figure 1 As shown, the present invention proposes a non-Gaussian permeability field hydrological geophysical joint inversion framework that couples hydraulic tomography (HT) and natural potential method. In this framework, CVAE in deep learning and ESMDA in data assimilation method are used to deal with the inversion problem of non-Gaussian field, and the estimation effect of non-Gaussian permeability field is improved by integrating natural potential data and limited hydraulic head data.

[0043] The above-mentioned natural potential data are obtained by the Self-Potential method (SP), which is used as additional distribution information of the heterogeneous parameters of the underground aquifer. SP is a passive geophysical method that can be directly implemented along with the pumping test. The SP method involves placing a series of non-polarized electrodes on the surface or in an existing borehole, connecting these electrodes to a multi-channel voltmeter with high impedance (>10MOhm) and high sensitivity (0.1mV) to read the changes in the natural potential signal at the observation point. The SP electrical signal responds to the movement of underground fluids (the generated potential is called flow potential), and it is highly sensitive to groundwater flow (Revil and Jardani, 2013). During the pumping test, the potential difference between each SP electrode and the reference electrode is recorded to obtain the SP potential value. The obtained SP potential value is a response to the groundwater flow state, so the SP electrical signal can be used as an additional source of observation data for the inversion of underground aquifer parameters.

[0044] The joint inversion framework flow chart is as follows: Figure 1 As shown, Figure 1 In , i represents the number of data assimilation, NI represents the specified maximum number of data assimilation, K, Ss, θ, σ represent the permeability field, water storage field, porosity field and conductivity field respectively, and the hydraulic parameters to be estimated can be replaced according to the actual situation. Figure 1 , the method steps of the present invention are further explained.

[0045] Step 1: Before inversion, a set of non-Gaussian aquifer parameter field samples are used to train the CVAE network. After the trained CVAE network is established, the latent vector z that obeys the Gaussian distribution is used as the parameter to be estimated during the data assimilation process instead of the original non-Gaussian aquifer parameter field. In this step, the non-Gaussian aquifer parameter field samples can be generated by the following method but are not limited to this method. They are generated by multi-point geostatistical methods, and the pattern of the parameter field can be obtained by expert information.

[0046] Step 2: Initialize a set of estimates of the latent vector z by sampling from a standard normal distribution.

[0047] Step 3: Based on the estimated set of latent vectors z, the trained CVAE decoder is used to reconstruct the corresponding estimated set of non-Gaussian aquifer parameter fields.

[0048] Step 4: According to the reconstructed non-Gaussian aquifer parameter field estimation set, during the harmonic pumping test, the groundwater flow model and the natural potential forward model are run in sequence to obtain hydraulic head simulation data and natural potential simulation data. Among them, the SP signal is not only directly controlled by the permeability coefficient and conductivity, but also affected by the groundwater flow velocity, so other hydraulic parameters will also indirectly affect the SP data, and before solving the natural potential model, the groundwater flow model needs to be solved first.

[0049] The groundwater flow forward model is now further explained.

[0050] The governing equation for three-dimensional transient saturated groundwater flow with external sources / sinks is expressed as follows:

[0051]

[0052] Among them, S s is the water storage rate (m -1 ) ; H is hydraulic head (m); t is time (s); represents the divergence operator; represents the gradient operator; K is the permeability coefficient (m / s); W represents the source and sink per unit volume (s -1 ).

[0053] Formula (1) is subject to the following initial and boundary conditions:

[0054] H=H 0 at t=0 (2)

[0055] H=H D at Γ D (3)

[0056]

[0057] In the formula, H 0 is the initial water head (m); H D is the given water head on the Dirichlet boundary (m); Γ D is the Dirichlet boundary; n represents the external unit vector perpendicular to the Neumann boundary; Γ N is the Neumann boundary.

[0058] During the harmonic pumping test, a harmonic hydraulic tomography (HHT) method was used, and the harmonic pumping rate was expressed as follows:

[0059] Q=-q 0 cos(2πt / T) (5)

[0060] Among them, q0 Indicates the peak value of volume flow rate (kg / s); T is the period (s).

[0061] The natural potential forward model is now further explained.

[0062] For the spontaneous potential forward model, the total charge flux in the porous medium is given by the current density j (A / m 2 ) represents the sum of the conduction current density and the source current density:

[0063] j=j c +j s (6)

[0064]

[0065]

[0066] Among them, j c is the conduction current density (A / m 2 ), can be obtained by Ohm's law; j s is the source current density (A / m 2 ), σ represents the electrical conductivity of the porous medium (S / m); represents the electric potential (V); q is the Darcy velocity (m / s), and represents the effective charge density (C / m 3 );parameter The following empirical model can be used for calculation (related to the intrinsic permeability):

[0067]

[0068] Where k is the intrinsic permeability (m 2 ), which can be obtained by converting the permeability coefficient K, that is, K = ρgk / η (ρ, η represent the density and dynamic viscosity of pore water respectively).

[0069] Combining formula (6) to formula (8) and the continuity equation of charge (i.e. ), the natural potential can be obtained The governing equations of the field are:

[0070]

[0071] The boundary conditions are as follows:

[0072]

[0073]

[0074] The solution of the spontaneous potential forward model requires the provision of the conductivity (σ) field. Under general conditions, σ can be inferred by the Electrical Resistivity Tomography (ERT) method (see Ahmed et al., 2016; Kang et al., 2020 for details).

[0075] The natural potential forward model is coupled with the groundwater flow forward model through the Darcy velocity q and the intrinsic permeability k to form the hydrogeophysical forward model of the present invention. The hydrogeophysical forward model is combined with the CVAE-ESMDA inversion algorithm to obtain the complete non-Gaussian field hydrogeophysical joint inversion framework proposed by the present invention.

[0076] Step 5: Based on the hydraulic head simulation data and natural potential simulation data obtained in step 4, combined with the hydraulic head observation data and the natural potential observation data, the ESMDA method is used to iteratively update the latent vector z, and steps 3-5 are repeated until the maximum number of data assimilation times NI is reached.

[0077] Step 6: Based on the posterior set of the latent vector z, the CVAE decoder is used to reconstruct the posterior set of the non-Gaussian aquifer parameter field to be estimated, that is, the optimal estimation result.

[0078] The effectiveness of the method of the present invention is further illustrated with reference to examples.

[0079] Construct a hypothetical non-Gaussian confined aquifer (see Figure 2 ) as the research object, in Figure 2 In the figure, (a) is a three-dimensional view of the hypothetical study area. The light and dark quadrilaterals represent hydraulic tomography injection holes located at z = 3m and z = 6m, respectively. The dark surface at x = 0m in the study area represents the tracer source in the solute transport scenario. (b) is a plan view of the aquifer. The red cross marks represent boreholes, each of which contains eight sampling ports with a spacing of 1m (z = 1, 2, ..., 8m). The black dots represent non-polarized electrodes used to record natural potential signals. These 100 electrodes are installed at a depth of 0.4m below the ground. Electrodes are installed at the aquifer observation hole locations, but they are not marked in the figure for clarity. The hollow circles represent three monitoring wells, whose observation data are not used for inversion but only for cross-validation of the results.

[0080] The aquifer is discretized into 40×40×10=16000 grid blocks, each of which is 1m×1m×1m in size. The aquifer consists of two lithofacies (see Figure 3 ) constitutes, in Figure 3 In the figure, (a) is the original training image used to generate CVAE training samples. (b) and (c) are the distributions of the lithofacies reference field, assuming that the permeability field (K) has the same distribution as the lithofacies. Figure 3 In the figure, the area labeled 1 and the surrounding background area labeled 0 represent the high permeability preferential flow channel and the low permeability background lithofacies, respectively. For the high permeability preferential flow channel, the K value is set to 10m / d, while the K value of the background medium is set to 0.4m / d. The reference (true) distribution of geological lithofacies is shown in Figure 3 As shown in b and c in Figure 1, porosity (θ), water storage rate (Ss) and conductivity (σ) are assumed to be unknown parameters and their distribution is consistent with the lithofacies. Therefore, in the inversion process, the distribution of lithofacies is first estimated, and then the lithofacies are filled with parameter values ​​to reconstruct the parameter field to be estimated. Table 1 summarizes the parameter values ​​used in the model.

[0081] Table 1 shows the parameter settings in the numerical experiments

[0082] parameter Value Regional scale (m) 40×40×10 Grid size (m) 1.0×1.0×1.0 Number of unknown variables 16000K Groundwater flow and spontaneous potential forward modeling Hydraulic gradient under steady flow conditions 0.006 Initial water head (m) 20 Permeability coefficient (m / d) High permeability lithofacies 10 Background lithofacies 0.4 <![CDATA[Water storage rate (m -1 )]]> High permeability lithofacies <![CDATA[2×10 -5 ]]> Background lithofacies <![CDATA[1×10 -4 ]]> Porosity High permeability lithofacies 0.3 Background lithofacies 0.4 Conductivity (S / m) High permeability lithofacies 0.001 Background lithofacies 0.025 Observation time 10,20,30,…,60s Number of sampling holes 64 Number of electrodes 164 Observation data volume Example 1 HHT 2688 Example 2 SP 7488 Example 3 Joint 10176 Solute transport model Simulation duration (d) 500 Tracer release time (d) 100 Observation time 10,20,30,…,500d Release source concentration (mg / L) 100 Longitudinal diffusion coefficient (m) 0.5 Horizontal diffusion coefficient (m) 0.05 Molecular diffusion coefficient 0

[0083] In order to weaken the influence of boundary conditions on hydraulic head and electrical signals, a large homogeneous buffer zone with constant hydraulic parameters and electrical conductivity is used to surround the study area.

[0084] For the forward modeling of groundwater flow, the top and bottom sides of the buffer zone are impermeable boundaries. The sides at x=-180m and x=220m are set as constant head boundaries of h=28m and h=20m, respectively. The remaining sides are zero flux boundary conditions. The initial head of the entire aquifer is set to 20m. Hydrogeophysical tests (i.e., HHT and SP) are carried out after the water flow reaches a steady state, and the hydraulic gradient in the steady state is about 0.006. For the forward modeling of natural potential, the top side of the buffer zone is an insulating boundary condition, and the remaining sides are set as zero potential boundaries.

[0085] Eleven boreholes were set up, each of which contained eight sampling ports with depths of z = 1, 2, ..., 8 m and a spacing of 1 m. Figure 2 The three boreholes marked with circles in b were not used in the inversion phase and were only used for cross-validation of the hydraulic head prediction. 100 non-polarized electrodes were placed 0.4 m below the ground surface and 64 non-polarized electrodes were placed at the sampling ports of the eight boreholes used for the inversion. Reference electrodes were placed in areas not affected by the pumping test.

[0086] like Figure 3 At the water injection hole shown in a, T = 60s and q 0=5kg / s were used to conduct eight harmonic pumping tests in sequence. During each pumping test, for simplicity, only the aquifer response of 60s period was simulated, and the head changes at the aquifer observation holes (except the boreholes currently used for pumping and injection) were observed every 10 seconds. In addition, 100 surface electrodes and 64 borehole electrodes were used to simultaneously record the SP potential anomalies caused by the harmonic pumping test. In this way, a total of 2688 heads and 7488 SP data were obtained from 8 sets of HHT and SP tests. Gaussian white noise with standard deviations of 0.1m and 0.1mV, respectively, was added to the heads and SP data calculated by the reference model (reference K field forward model) as hypothetical observation values.

[0087] Three examples are designed to verify the advantages of fusing hydraulic head and SP data. In Example 1 (Case 1), only hydraulic head is used to estimate the non-Gaussian K field. In Example 2 (Case 2), only inversion SP data is considered. In Example 3 (Case 3), hydraulic head and SP observation data are used simultaneously to characterize the non-Gaussian permeability field.

[0088] The architecture of the CVAE network is shown in Table 2. The input data of CVAE is a lithofacies image of a non-Gaussian permeability field, where high permeability lithofacies and background lithofacies correspond to values ​​1 and 0, respectively. The encoder consists of six convolutional layers. The decoder is a mirror-symmetric network of the encoder with a deconvolution layer.

[0089] The training set is from Figure 3 The images are generated by cropping 40×40×10 images from the 340×200×80 training images shown in a. By using a cropping step of (1,1,1), a total of 53,856 petrographic field images are obtained. Among them, 48,380 samples are used for training, and the remaining samples are used as test samples to verify the performance of CVAE. The reconstruction accuracy of CVAE is evaluated using the following expression:

[0090]

[0091] Where y represents the reconstruction accuracy of CVAE; G is the total number of grid blocks (16,000 in this case); S n is the size of the training set or test set; and Represent the lithofacies value (0 or 1) of the jth pixel of the ith sample before and after CVAE reconstruction, respectively.

[0092] The dimension of the latent vector is 500. The reference field ( Figure 3b and c) are randomly selected from the test samples. Based on the Adam optimizer (Kingma and Ba, 2014), 120 epochs of training are performed on an NVIDIA Tesla V100 GPU with a learning rate of 0.0003 and a batch size of 32. After training, the CVAE decoder is combined with ESMDA for non-Gaussian data assimilation. In ESMDA, the ensemble size is 500. According to the suggestion of Emerick and Reynolds (2013), the number of iterations is set to four, and the expansion factors of each iteration are 9.33, 7, 4, and 2, respectively.

[0093] The reconstruction accuracy of the CVAE training set and test set is 0.9086 and 0.9108 respectively. Figure 4 Three lithofacies images randomly selected from the test set are shown, along with the corresponding CVAE reconstructed images and histograms of the latent variables. The results show that the trained CVAE can successfully reproduce the lithofacies spatial distribution pattern of the input image. In addition, the generated latent variables roughly follow a standard normal distribution. Using these normally distributed latent variables as the parameters to be estimated, data assimilation can be performed without violating the Gaussian assumption of ESMDA. Figure 4 Three lithofacies images randomly generated by the CVAE decoder are also shown (in the fourth column of the figure). The results show that the trained CVAE decoder is able to generate lithofacies distribution images that conform to the training sample pattern. In summary, CVAE can be used as an effective parameterization tool to handle non-Gaussian inversion problems. In addition, it reduces the dimensionality of the variables to be solved from 16,000 to 500, which greatly reduces the computational cost of the inversion process.

[0094] Table 2 shows the architecture of the CVAE network.

[0095]

[0096]

[0097] Among them, BN stands for Batch Normalization (Ioffe and Szegedy, 2015), and ReLU stands for Rectified Linear Unit function.

[0098] Figure 5 The ensemble mean field, a randomly chosen posterior field, and the ensemble variance field of ln K estimated using different types of datasets (Cases 1-3) are shown. Figure 6The residual field between the posterior mean field and the reference field is given. The second and third columns are the residual fields at depths of z = 2m and z = 9m, respectively; for a more intuitive evaluation, the NRMSE at z = 2m and z = 9m is given according to formula (14). For ease of comparison, the estimated results of the prior set are also shown. Since the probability of randomly generating high-permeability lithofacies and background lithofacies at the same pixel is equal, the ln K mean field estimated by the prior set tends to a constant field. In addition, since no data assimilation is performed, the uncertainty (variance) of the prior field is large.

[0099] In Example 1, the estimated ln K field roughly reproduces the main features of the reference field. However, due to the limited amount of head data, the estimated field fails to capture the small-scale features of the reference field. For example, Figure 5 The ln K structure in the black dot box in Example 1 only provides a smooth and fuzzy estimate, which fails to reflect the fine structure of the reference field. In addition, the variance field of Example 1 shows that the near-surface region with x>20m and y<20m has a high uncertainty, and the residual field map of Example 1 ( Figure 6 The residual at z = 9 m is relatively larger than that at z = 2 m, indicating that the HHT method cannot reproduce the lnK field near the surface with a higher resolution due to the limited number of head piezometers installed close to the surface.

[0100] Similar to Example 1, the lnK field estimated in Example 2 reproduces the main structure of the reference field, but fails to characterize some fine inhomogeneous features. Figure 6 It can be seen that, probably due to the low accuracy of SP data, the residual in the lower part of the aquifer (z = 2 m) in Example 2 is larger than that in Example 1. However, for the upper part of the aquifer (z = 9 m), due to the high density of non-polarized electrodes installed on the surface, SP can better estimate the non-Gaussian ln K field, and the uncertainty near the surface is low.

[0101] In Example 3, considering the importance of prior estimation to the data assimilation method, 20 latent vector posterior samples are taken from Example 1 and Example 2 to replace the same number of prior samples in Example 3, thereby accelerating the convergence of the joint inversion and improving the estimation accuracy. This strategy is very effective in highly nonlinear inversion problems (e.g., Wen and Chen, 2005; Chang et al., 2021). Compared with Examples 1 and 2, the lnK field estimated in Example 3 reproduces the reference field with higher resolution. For example, Figure 5 In the figure, the lnK field in the black dot box of Example 3 well describes the characteristics of the reference field. Figure 6It is shown that for Example 3, the residuals in both the shallow and deep parts of the aquifer are significantly reduced, indicating that the effects of hydraulic head and SP data are complementary. In addition, compared with using a single type of data set, the uncertainty (i.e., ensemble variance) of Example 3 is significantly reduced. This is because (1) the information contained in the two data sets is complementary, which can better constrain the inversion and reduce the uncertainty in the shallow and deep areas; (2) a better prior set is used in Example 3.

[0102] In order to further compare the lnK estimation results, a quantitative indicator, normalized root mean square error (NRMSE), is introduced to evaluate the accuracy of the estimated ln K set mean field. The formula of NRMSE is as follows:

[0103]

[0104] Among them, N u is the number of grids (16000 in this case); K e,i and K r,i are the estimated mean and reference values ​​of the permeability coefficient in the i-th grid respectively; K max and K min are the permeability values ​​corresponding to high permeability lithofacies and background rocks, respectively; the NRMSE of each example is shown in Table 3. Consistent with the previous results, the NRMSE of Example 1 and Example 2 are close, while the NRMSE of Example 3 (joint inversion) is the lowest, which shows that the non-Gaussian permeability field can be better characterized by jointly inverting hydraulic head and SP data.

[0105] Using three validation wells (see Figure 2 ) head data, we tested the head prediction ability of the mean field estimated by the permeability coefficient. The cross-validation results are shown in Figure 7 Among them, the hydraulic head predicted by Example 3 fits the reference hydraulic head best, which shows that the improved permeability field estimation result in Example 3 can better predict the movement of underground fluids.

[0106] The tracer migration test was performed to further evaluate the predictive ability of the permeability estimation field. For the solute transport model, the tracer was uniformly released on the x = 0 m side of the study area for 100 days (see Figure 2 ), the source concentration is 100 mg / L. The simulation time is set to 500 days. The parameters used in the tracer transport model are shown in Table 1. For each K field sample in the initial and posterior sets, the solute transport forward model is run to obtain the tracer concentration field corresponding to each sample.

[0107] Figure 8The predicted mean fields of the tracer concentration distribution at z = 5 m are shown for the initial ensemble and after 100, 150 and 200 days for the three example cases. Figure 8 The fourth column gives the residual field between the predicted concentration field and the reference concentration field at t=150d. The tracer plume predicted according to the initial set migrates in a uniform form and does not show any characteristics of the reference concentration field. In contrast, the tracer plumes predicted by Examples 1-3 reproduce the characteristics of the reference tracer plume to varying degrees. The prediction results using only one data set (Example 1 and Example 2) can roughly reflect the main structure of the reference concentration field, but cannot characterize the fine structure of the tracer plume. For example, in some areas of the two cases, large residuals are still shown. However, when the hydraulic head and SP data are jointly inverted (Example 3), the distribution of the tracer plume is better predicted, and the residual between the predicted concentration field and the reference concentration field is very small. In general, the tracer migration prediction results further show that the SP method is indeed helpful in characterizing complex non-Gaussian permeability fields and predicting groundwater flow and solute transport processes.

[0108] Table 3 shows the NRMSE of the ln K mean field

[0109]

[0110]

Claims

1. A method for identifying non-Gaussian aquifer parameters by fusing groundwater level and natural potential data based on convolutional neural network, characterized in that: The following steps are involved: Step 1: Train the CVAE network using non-Gaussian aquifer parameter field samples and initialize the latent vector z estimation set by sampling from the standard normal distribution; Step 2: Input the estimated set of latent vectors z into the trained CVAE decoder to reconstruct the corresponding estimated set of non-Gaussian aquifer parameter fields; Step 3: During the harmonic pumping test, the hydrogeophysical forward model is run based on the reconstructed non-Gaussian aquifer parameter field estimation set to obtain hydraulic head simulation data and natural potential simulation data; Step 4: Based on the hydraulic head simulation data and natural potential simulation data obtained in step 3, combined with the hydraulic head observation data and the natural potential observation data, the ESMDA method is used to iteratively update the latent vector z estimation set; repeat steps 2 to 4 until the maximum number of iterations is reached; Step 5: For the updated latent vector z estimation set, reconstruct the estimation result of the non-Gaussian aquifer parameter field through the CVAE decoder; The hydrogeophysical forward model is composed of a natural potential forward model coupled with a groundwater flow forward model through Darcy velocity q and intrinsic permeability k; The groundwater flow forward model is expressed as: Among them, S s is the water storage rate; H is the hydraulic head; t is the time; represents the divergence operator; represents the gradient operator; K is the permeability coefficient; W represents the source and sink per unit volume; The groundwater flow forward model is subject to the following initial and boundary conditions: H=H0 at t=0 (2) H=H D at Γ D (3) Where H0 is the initial water head; Γ D is the Dirichlet boundary; H D is Γ D Given water head on Γ N is the Neumann boundary; n is perpendicular to Γ N The exterior unit vector of The spontaneous potential forward model is expressed as: The boundary conditions are as follows: Where σ is the conductivity of the porous medium, represents the electric potential, Γ D is the Dirichlet boundary, Γ N is the Neumann boundary; q is the Darcy velocity, Represents the effective charge density generated by groundwater flow.

2. According to claim 1, a method for identifying non-Gaussian aquifer parameters based on convolutional neural network fusion of groundwater level and natural potential data, characterized in that: The non-Gaussian aquifer parameter field includes a permeability field, a water storage rate field, a porosity field and a conductivity field.

3. The method for identifying non-Gaussian aquifer parameters based on convolutional neural network fusion of groundwater level and natural potential data according to claim 1 is characterized by: The operation of the hydrogeophysical forward model includes first solving the groundwater flow forward model and then solving the natural potential forward model.

4. The method for identifying non-Gaussian aquifer parameters based on convolutional neural network fusion of groundwater level and natural potential data according to claim 1 is characterized by: The CVAE network includes a CVAE encoder and a CVAE decoder; the encoder consists of six convolutional layers; the CVAE decoder is a mirror-symmetric network of the encoder with a deconvolution layer.

Images