A three-dimensional wavelet-Galerkin forward modeling method for DC resistivity

Abstract

The present invention belongs to the field of geophysical exploration, and is specifically a three-dimensional direct current resistivity wavelet Galerkin forward modeling method, which includes: performing non-uniform grid division on a geoelectric model calculation area into a plurality of grid units, obtaining grid unit numbers, node numbers and coordinates; reading designed geoelectric model parameters, including background layer parameters, grid unit numbers, node numbers and coordinates; calculating the primary field potential value of the underground background medium of the designed background layer; performing cyclic calculation of unit coefficient matrices for all grid units, and then synthesizing the coefficient matrices of all grid units as a whole; calculating the right-hand side term of the direct current resistivity method wavelet Galerkin linear equation group according to the primary field potential value; loading essential boundary conditions, solving the direct current resistivity method wavelet Galerkin linear equation group, obtaining the secondary field potential of each node, and then obtaining the total field potential value. The problem of slow calculation speed and low calculation accuracy existing in the existing method is solved.

Description

Technical Field

[0001] The invention belongs to the field of geophysical exploration, and specifically relates to a direct current resistivity wavelet Galerkin three-dimensional forward modeling method. Background Art

[0002] With the rapid and sustained development of economic construction, the shortage of mineral resources has become increasingly prominent. The current traditional mineral exploration and development technology is insufficient in theoretical research, and has limitations in detection depth and accuracy. As the target is more hidden and complex, the difficulty of oil and gas exploration and development is increasing accordingly. Therefore, conducting research on mineral resource exploration technology and improving exploration accuracy are of great significance to improving the benefits of mineral resource exploration and development.

[0003] As we all know, DC resistivity method, as one of the oldest exploration methods in geophysical exploration, has long been widely used in mineral resource exploration, hydrogeological survey, engineering survey and other fields. With the increasing complexity of exploration, how to improve the accuracy of exploration has become a hot research issue at home and abroad. Forward modeling is the basis of inversion. Therefore, the study of high-precision three-dimensional resistivity forward modeling is of great theoretical significance and practical value for achieving efficient three-dimensional inversion. So far, the commonly used DC resistivity numerical simulation methods mainly include integral equation method (IEM), finite difference method (FDM), finite element method (FEM), etc.

[0004] There are problems such as slow calculation speed and low calculation accuracy in the three-dimensional forward modeling of DC resistivity method. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a direct current resistivity wavelet Galerkin three-dimensional forward modeling method to solve the problems of slow calculation speed and low calculation accuracy in the existing method.

[0006] The present invention is achieved in this way.

[0007] A direct current resistivity wavelet Galerkin three-dimensional forward modeling method, the method comprising:

[0008] The geoelectric model calculation area is divided into multiple grid units by non-uniform grid division, and the grid unit number, node number and coordinate are obtained;

[0009] Read the designed geoelectric model parameters, including background layer parameters, grid unit numbers, node numbers and coordinates;

[0010] Calculate the primary field potential value of the underground background medium of the designed background layer;

[0011] Calculate the unit coefficient matrix for all grid cells in a loop, and then synthesize the coefficient matrix of all grid cells as a whole;

[0012] Calculate the right-hand side of the wavelet Galerkin linear equations of the DC resistivity method based on the primary field potential value;

[0013] The essential boundary conditions are loaded, and the DC resistivity method wavelet Galerkin linear equations are solved to obtain the secondary field potential of each node, and then the total field potential value is obtained.

[0014] Further:

[0015] The DC resistivity method wavelet Galerkin linear equations are:

[0016] K s u s =b

[0017] in:

[0018]

[0019]

[0020] K s is the coefficient matrix, u s To solve the field value of each node in the domain, b is the right-hand side term. Both are 2-term connection coefficients.

[0021] Further:

[0022] The DC resistivity method wavelet Galerkin linear equations are obtained by the following method:

[0023] The boundary value problem expression satisfied by the quadratic field potential using the Dirichlet boundary condition is:

[0024]

[0025] where u s is the secondary field potential to be solved, u p is the primary field potential, σ is the conductivity of the medium, σ p is the conductivity of the background medium, Δσ is the abnormal conductivity, Δσ=σ-σ p ;

[0026] Secondary field potential u s and the primary field potential u p The wavelet scaling function is used to expand them as follows:

[0027]

[0028]

[0029] Substituting into the boundary value problem expression satisfied by the quadratic field potential, we get:

[0030]

[0031] make becomes:

[0032]

[0033] in

[0034]

[0035] Multiply both sides by ∫φ m,p‘ (x)dx∫φ m,q‘ (y)dy∫φ m,r‘ (z)dz gives:

[0036]

[0037] in

[0038]

[0039]

[0040]

[0041] in is the product integral of the wavelet scaling function and its derivative, called the 2-term connection coefficient, which is expressed as follows:

[0042]

[0043]

[0044]

[0045] According to the orthogonal property of wavelet scaling function, we can get:

[0046] ∫φ m,p‘ (x)φ m,p (x)dx=δ p,p′

[0047] ∫φ m,q‘ (y)φ m,q (y)dy=δ q,q′

[0048] ∫φ m,r‘ (z)φ m,r (z)dz=δ r,r′

[0049] Substituting into the equation have to

[0050]

[0051] in

[0052]

[0053] Further simplification yields:

[0054] K s u s =b.

[0055] Furthermore: the wavelet scaling function is a Daubechies wavelet scaling function.

[0056] Furthermore: a non-uniform geometric grid is used to discretize the study area.

[0057] Furthermore: a uniform half space or a horizontal layered medium is used as a background medium to calculate a primary field potential value.

[0058] Furthermore: the parallel direct solver MUMPS is used to solve the wavelet Galerkin linear equations of the DC resistivity method.

[0059] Further: The essential boundary condition is the Dirichlet boundary condition, and the potential value at the boundary is: u s | Γ0 =0.

[0060] Compared with the prior art, the present invention has the following beneficial effects:

[0061] The present invention proposes for the first time a three-dimensional numerical simulation algorithm of DC resistivity based on the wavelet Galerkin method of Daubechies wavelet function. The algorithm can effectively avoid the influence of singularity brought by the introduction of field sources, and discretize the research area with non-uniform grids, encrypt near the field sources and anomalies, and use sparse grids far away from the field sources and anomalies. This is conducive to the construction of geometrically complex geoelectric models and better simulation of real geological conditions. On the other hand, it can avoid the same degree of encryption in the entire area, effectively reducing the number of unknowns and thus reducing the amount of calculation. Therefore, the application of the above-mentioned technologies is combined to achieve efficient and fast numerical simulation of isotropic media by DC resistivity method.

[0062] The invention is applicable to numerical simulation research of geophysical exploration methods such as borehole electrical method, ocean direct current method, induced polarization, etc. with complex geometry and physical property distribution.

[0063] The present invention proposes an efficient numerical simulation technology to solve the three-dimensional forward modeling problem of DC resistivity, which can improve the calculation accuracy and reduce the time. First, it is assumed that the underground medium is an isotropic medium. And the DC resistivity wavelet Galerkin equation under the condition of isotropic conductivity based on the secondary field potential is derived for the first time. The algorithm can effectively avoid the influence of source point singularity and improve the accuracy of the numerical solution. In order to better simulate complex geological structures, the present invention uses non-uniform grids to discretize the study area, which can construct complex geological structures and physical property distributions, while avoiding the same degree of encryption of the entire area, reducing the amount of calculation. In view of the problem of slow convergence of conventional iterative algorithms for solving linear equations, it is proposed to use a parallel direct solver MUMPS to directly solve the linear equations, which effectively improves the calculation efficiency. The conventional iterative solution method further improves the accuracy of the numerical solution. Therefore, the present invention can realize high-precision and fast numerical simulation of DC resistivity under the condition of isotropic conductivity. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 1 is a scale function (a) and a wavelet function diagram (b) of the wavelet Galerkin method provided by an embodiment of the present invention;

[0065] Figure 2 is a flow chart of a method provided by an embodiment of the present invention;

[0066] Figure 3 is a horizontal layered isotropic conductivity medium model established in an embodiment of the present invention;

[0067] Figure 4 It is a comparison between the quasi-analytical solution of the isotropic layered model of the embodiment of the present invention and the wavelet Galerkin numerical solution and the finite element numerical solution;

[0068] Figure 5 is a comparison of relative errors of the wavelet Galerkin numerical solution and the finite element numerical solution and the quasi-analytical solution of the isotropic medium layered model according to the embodiment of the present invention;

[0069] Figure 6 The model of an abnormal body in an isotropic medium established in an embodiment of the present invention is as follows: (a) a vertical cross-sectional schematic diagram of the isotropic medium model at y=0; (b) a planar schematic diagram of the isotropic medium model at z=5m;

[0070] Figure 7 The apparent resistivity plane diagrams of the anomaly body model in the isotropic medium established in the embodiment of the present invention measured by the E-SCAN method are as follows: (a) the apparent resistivity plane diagram at z=0 with the power supply point coordinates of (0,0,0)m; (b) the apparent resistivity plane diagram at z=0 with the power supply point coordinates of (0,-12,0)m; (c) the apparent resistivity plane diagram at z=0 with the power supply point coordinates of (-12,-12,0)m. DETAILED DESCRIPTION

[0071] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0072] See also Figure 2 , the DC resistivity wavelet Galerkin three-dimensional forward modeling method includes the following steps:

[0073] 1) Starting from the Laplace equation, the wavelet Galerkin method is used to obtain the DC resistivity wavelet Galerkin linear equations under the condition of isotropic conductivity;

[0074] 2) Grid the geoelectric model study area to obtain parameters such as grid unit number, node number and coordinates;

[0075] 3) Read the designed geoelectric model parameters, including background layer parameters, conductivity model, grid unit number, node coordinates, etc.;

[0076] 4) Calculate the primary field potential corresponding to the background layer;

[0077] 5) Calculate the coefficient matrix for all grid cells in a loop, then synthesize the cell coefficient matrix of all cells as a whole. At the same time, calculate the right-hand side of the DC resistivity wavelet Galerkin linear equation group based on the primary potential calculated from the background layer.

[0078] 6) Load the essential boundary conditions, solve the linear equations, and obtain the quadratic field potential of each node;

[0079] 7) Add the secondary field potential to the primary field potential to obtain the total field potential value of all nodes.

[0080] Among them, starting from the boundary value problem of DC resistivity method for isotropic conductive media,

[0081] In order to avoid the influence of the source point singularity, the secondary field algorithm is introduced to carry out the three-dimensional forward modeling of DC resistivity. Here, the boundary value problem satisfied by the secondary field potential considering the Dirichlet boundary condition is directly given:

[0082]

[0083] where u s is the secondary field potential to be solved, u p is the primary field potential, which can be calculated by taking a uniform half space or full space or a layered medium as the background medium (step 4) (calculate the primary field potential). σ is the conductivity of the medium, σ pis the conductivity of the background medium, Δσ is the abnormal conductivity, that is, Δσ=σ-σ p .

[0084] Secondary field potential u s and the primary field potential u p Wavelet scaling functions are used for expansion respectively. The wavelet scaling function uses Daubechies wavelet scaling function.

[0085] The Daubechies wavelet scaling function is defined as:

[0086]

[0087] The corresponding wavelet function is defined as

[0088]

[0089] Where L is the support length of the wavelet, p j is the filter coefficient. For the scaling function of Daubechies wavelet, see Figure 1 , and by the multi-resolution property of wavelet, we have

[0090]

[0091] The wavelet basis function of the three-dimensional problem can be obtained by the product of the one-dimensional compact support wavelet basis functions in the x, y, and z directions as shown below:

[0092] φ r,p,q (x,y,z)=φ r (x)φ p (y)φ q (z) (5)

[0093] where φ r (x),φ p (y),φ q (z) is the wavelet scaling function in the x, y, and z directions.

[0094] In the present invention, the research area is divided into a finite number of hexahedral units e. In each unit, the conductivity is fixed and the secondary field potential u in the unit is constant. s and the primary field potential u p The wavelet scaling function is used to expand (corresponding to the grid division in step 2)

[0095]

[0096]

[0097] Substituting equations (6) and (7) into equation (1), we can obtain:

[0098]

[0099] make

[0100] Equation (8) can be transformed into

[0101]

[0102] in

[0103]

[0104] To further simplify equation (9), multiply both sides of equation (9) by ∫φ m,p‘ (x)dx∫φ m,q‘ (y)dy∫φ m,r ‘ (z)dz,

[0105] Equation (10) is further transformed into

[0106]

[0107] in:

[0108]

[0109]

[0110]

[0111] in is the product integral of the wavelet scaling function and its derivative, generally called the 2-term connection coefficient. It can be expressed as follows:

[0112]

[0113]

[0114]

[0115] According to the orthogonal property of wavelet scaling function

[0116] ∫φ m,p‘ (x)φ m,p (x)dx=δ p,p′ (18)

[0117] ∫φ m,q‘ (y)φ m,q (y)dy=δ q,q′ (19)

[0118] ∫φ m,r‘(z)φ m,r (z)dz=δ r,r′ (20)

[0119] Substituting (18)-(20) into equation (9), we can obtain

[0120]

[0121] in

[0122]

[0123] Then equation (21) can be further simplified to

[0124] K s u s =b (23)

[0125] Where (step 5) (calculate the coefficient matrix and the right-hand side term)

[0126]

[0127]

[0128] By solving equation (23), the secondary field potential can be obtained, and then the total field potential can be calculated.

[0129] In order to solve equation (23), it is also necessary to add corresponding boundary conditions. The present invention adopts the Dirichlet boundary condition, that is, the potential value at the boundary is 0: (Step 6) (Adding essential boundary conditions)

[0130] u s | Γ0 =0 (24)

[0131] Solution of the linear equations in step 6:

[0132] Compared with the Krlov subspace iteration algorithm, the direct solution method is widely used to solve large sparse linear equations due to its advantages of high solution accuracy and good stability. In particular, it can solve multiple right-hand side terms at one time, which can effectively speed up the speed of solving forward problems. Therefore, the present invention adopts the parallel direct solver MUMPS to solve the linear equations.

[0133] In order to verify the correctness of the three-dimensional forward modeling of DC resistivity based on the wavelet Galerkin method (WGM) proposed in this paper, a two-layer horizontal layered medium model is designed. Figure 3As shown, the thickness of the first layer is 100m, the resistivity is 100 ohm-meter, and the resistivity of the second layer is 1000 ohm-meter. The origin of the coordinate is located at the center of the landmark, and a symmetrical four-pole device is used for measurement, with a current of 1A. The digital filtering algorithm is used to calculate the quasi-analytical solution, and the numerical solution calculated by the wavelet Galerkin method is compared with the analytical solution and the numerical solution calculated by the finite element method. The comparison results are shown in Figure 4 As shown. Figure 4 It can be seen that the numerical solution of the wavelet-Galerkin method is basically consistent with the analytical solution and the finite element solution, which verifies the correctness of the present invention. Figure 5 Table 1 is a comparison of the accuracy of the two numerical simulation algorithms, WGM and FEM. It can be clearly seen that the accuracy of the WGM solution is higher than that of the finite element method. This shows that the DC resistivity WGM forward algorithm proposed in the present invention is correct and feasible. Table 1 is a comparison of the time consumed by the finite element method and the WGM calculation. From Table 1, it can be seen that the wavelet Galerkin algorithm proposed in the present invention is faster than the finite element method. This further shows that the wavelet Galerkin method can be used as a new approach to numerical simulation of the DC resistivity method.

[0134] Table 1 Comparison of calculation time of different numerical simulation methods

[0135]

[0136] In order to further verify the effectiveness of the method of the present invention, a model with a low-resistance anomaly in a uniform half-space is designed, such as Figure 6 (a) and Figure 6 (b) As shown. The length, width and height of the low-resistance anomaly are all 6m, the burial depth is 3m, and the center position is (0, 0, 6). The resistivity is 50 ohm-meter, and the resistivity of the surrounding rock is 1000 ohm-meter. The coordinate origin is set at the center of the model. The power supply current is 1A, and three power supply points are designed for power supply respectively, with specific coordinates of (0, 0, 0), (0, -12, 0), and (-12, -12, 0). The E-SCAN method is used for measurement, and the calculation results are as follows Figure 5 shown.

[0137] Figure 7 This is the apparent resistivity plane diagram when measured using the E-SCAN method. Figure 7 (a)-(c) are plane diagrams of apparent resistivity anomalies with power supply points of (0,0,0), (0,-12,0), and (-12,-12,0), respectively. Figure 7 The apparent resistivity plane diagram of the isotropic medium containing anomaly model established in the embodiment of the present invention when measured by the E-SCAN method: (a) the apparent resistivity plane diagram at z=0 with the power supply point coordinates (0,0,0)m; (b) the apparent resistivity plane diagram at z=0 with the power supply point coordinates (0,-12,0)m; (c) the apparent resistivity plane diagram at z=0 with the power supply point coordinates (-12,-12,0)m. Figure 7 (a) It can be seen that when the power supply point is located directly above the low-resistance body, the apparent resistivity plane diagram presents a regular symmetrical circular anomaly, and the size of the anomaly corresponds to the projection position of the anomaly body on the surface. Figure 7 (b)-(c) show that high and low resistance anomalies appear in pairs, and the centers of high and low resistance anomalies are on the same straight line as the power supply point. Table 2 shows the CPU time comparison between the finite element method and the wavelet-galerkin method. It can be found from Table 2 that under isotropic conditions, the time consumed by WGM is less than that of FEM.

[0138] Table 2 Comparison of calculation time of different numerical simulation methods under different models

[0139]

[0140] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A DC resistivity wavelet Galerkin three-dimensional forward modeling method, characterized in that: The method includes: The geoelectric model calculation area is divided into multiple grid units by non-uniform grid division, and the grid unit number, node number and coordinate are obtained; Read the designed geoelectric model parameters, including background layer parameters, grid unit numbers, node numbers and coordinates; Calculate the primary field potential value of the underground background medium of the designed background layer; Calculate the unit coefficient matrix for all grid cells in a loop, and then synthesize the coefficient matrix of all grid cells as a whole; Calculate the right-hand side of the wavelet Galerkin linear equations of the DC resistivity method based on the primary field potential value; Load the essential boundary conditions, solve the DC resistivity method wavelet Galerkin linear equations, obtain the secondary field potential of each node, and then obtain the total field potential value; The DC resistivity method wavelet Galerkin linear equations are: K s u s =b in: K s is the coefficient matrix, u s To solve the field value of each node in the domain, b is the right-hand side term. are 2-term connection coefficients; the DC resistivity method wavelet Galerkin linear equations are obtained by the following method: The boundary value problem expression satisfied by the quadratic field potential using the Dirichlet boundary condition is: where u s is the secondary field potential to be solved, u p is the primary field potential, σ is the conductivity of the medium, σ p is the conductivity of the background medium, Δσ is the abnormal conductivity, Δσ=σ-σ p ; Secondary field potential u s and the primary field potential u p The wavelet scaling function is used to expand them as follows: Substituting into the boundary value problem expression satisfied by the quadratic field potential, we get: make becomes: in Multiply both sides by ∫φ m,p‘ (x)dx∫φ m,q‘ (y)dy∫φ m,r‘ (z)dz gives: in in is the product integral of the wavelet scaling function and its derivative, called the 2-term connection coefficient, which is expressed as follows: According to the orthogonal property of wavelet scaling function, we can get: ∫φ m,p‘ (x)φ m,p (x)dx=δ p,p′ ∫φ m,q‘ (y)ϕ m,q (y)to=δ q,q′ ∫φ m,r‘ (of)φ m,r (z)dz=δ r,r′ Substituting into the equation have to in Further simplification yields: K s u s =b。 2. The method according to claim 1, characterized in that: The wavelet scaling function is the Daubechies wavelet scaling function.

3. The method according to claim 1, characterized in that: The study area is discretized using non-uniform geometric grids.

4. The method according to claim 1, characterized in that: A uniform half space or a horizontal layered medium is used as the background medium to calculate the primary field potential value.

5. The method according to claim 1, characterized in that: The parallel direct solver MUMPS is used to solve the wavelet Galerkin linear equations of the DC resistivity method.

6. The method according to claim 1, characterized in that: The essential boundary condition is the Dirichlet boundary condition, and the potential value at the boundary is: u s | Γ0 =0.

Images