Well-to-surface differential electromagnetic three-dimensional inversion method considering well casing effect and application

Abstract

The application provides a well-ground differential electromagnetic three-dimensional inversion method considering the effect of well casing, and relates to the geophysical logging technical field; an additional term is added to modify the average resistivity formula, so that numerical simulation containing the steel casing is realized without refining the grid, and the resistivity model containing the casing effect is inverted by the Gauss-Newton method. The resistivity model parameters containing the casing effect are inverted, the well-ground differential three-dimensional electromagnetic inversion considering the effect of well casing is realized, and the precision and accuracy of the well-ground electromagnetic inversion are further improved.

Description

Technical Field

[0001] This invention relates to the field of geophysical logging technology, specifically to a well-to-surface differential electromagnetic three-dimensional inversion method and its application that takes into account the effect of the well casing. Background Technology

[0002] Well-to-ground electromagnetic method refers to a type of electromagnetic depth sounding method that supplies high-power alternating current in the well and receives electromagnetic field signals on the surface. As a high-precision exploration method applied downhole, this method was introduced from Russia in the last century and has undergone a series of tests and applications in oil fields. It has achieved certain results in delineating oil and gas reservoirs, studying reservoir distribution, and dynamically monitoring water injection or grouting in oil wells, and has become a new technical method for studying complex underground geoelectric models.

[0003] Well-to-surface electromagnetic (WSMS) 3D inversion combines the advantages of well-based electromagnetic sounding and surface electromagnetic exploration. Through 3D inversion algorithms, it enables high-precision, high-resolution imaging of subsurface media. In well-to-surface 3D inversion, optimization algorithms such as finite difference, nonlinear conjugate gradient, or quasi-Newton methods are typically used to solve the objective function, and regularization techniques are employed to improve the stability and reliability of the inversion solution. Furthermore, key technologies such as sparsity regularization and data fitting terms are also widely applied in well-to-surface 3D inversion to enhance accuracy and efficiency. This technology not only has broad application prospects but also plays a crucial role in resource exploration, environmental monitoring, and engineering surveying.

[0004] However, in electromagnetic inversion, the steel sleeve can significantly affect the received electromagnetic signal, including changes in amplitude and phase difference. Therefore, the influence of the steel sleeve must be fully considered during electromagnetic inversion to improve its accuracy and precision.

[0005] Therefore, developing a novel well-to-surface differential electromagnetic three-dimensional inversion method that fully considers the influence of the steel casing on the received electromagnetic signals is an urgent problem that researchers in this field need to solve. Summary of the Invention

[0006] To address the aforementioned problems, this invention provides a well-to-surface differential electromagnetic three-dimensional inversion method that considers the effect of the well casing. It introduces a forward modeling method that considers casing resistivity into the Gauss-Newton method inversion, aiming to reduce unnecessary errors in the inversion caused by the influence of casing resistivity.

[0007] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0008] On the one hand, the present invention provides a well-to-surface differential electromagnetic three-dimensional inversion method that takes into account the effect of the well casing, comprising the following steps:

[0009] Step S1: Read the model file and data file of the inverted region to obtain the measured data;

[0010] Step S2: Based on Maxwell's equations in the frequency domain and Ohm's law, rewrite Maxwell's equations as Helmholtz equations, and use the pseudo-finite volume method to discretize the spatial domain, obtaining the discrete form:

[0011] Step S3: Reconstruct the average matrix based on whether the source device is completely located inside the sleeve.

[0012] Step S4: Use the average matrix A obtained in step S3 c2e Substituting into the discrete equation described in step S2, we obtain the simulation equation. The overall coefficient matrix A for each transmission frequency is obtained by analyzing the pseudo-finite volume method.

[0013] Step S5: Based on the measured data obtained in Step S1 and the simulation equation obtained in Step S4, construct the objective function and obtain the normal vector equation H. k δm k =-g k First-order partial derivative g k and Hessian matrix H k ;

[0014] Step S6: Obtain the sensitivity matrix using the quasi-forward modeling method, then solve for the Hessian matrix obtained in step S5, and finally solve the normal vector equation using the PCG method to calculate the update direction δm of the model vector in each iteration. k ;

[0015] Step S7: Obtain the model update vector: Use the inverse tracking method to perform an inaccurate line search on the update direction to determine the update step size α. k The model vector used in the next iteration is m. k+1 =m k +α k δm k ;

[0016] Step S8: Repeat steps S5, S6, and S7 until the requirements are met.

[0017] Preferably, in step S1, the model file and data file include the following: mesh division, steel sleeve length, resistivity, cross-sectional area, location, transmitting and receiving devices, transmission frequency, and measured data.

[0018] Preferably, step S2 is performed while neglecting displacement current.

[0019] Preferably, in step S2, the discrete expression is:

[0020] (Cf2e C e2f +iωμdiag(A c2e m))E=-iωμJ s (1);

[0021] Among them, C f2e With C e2f These are the discrete matrices of the twisting operator from the center of the mesh to the edge and the discrete matrices of the twisting operator from the edge to the center of the mesh, respectively.

[0022] ω is the angular frequency at the current desired frequency;

[0023] A c2e This is the average matrix from the center of the grid cell to the edge of the grid.

[0024] m is the conductivity vector discretized onto the grid cells, and its unit is S / m;

[0025] E is the electric field vector, and its unit is V / m;

[0026] J s For external field sources, the unit is A / m. 2 ;

[0027] μ is the magnetic permeability, with units of H / m;

[0028] i is the imaginary unit.

[0029] Preferably, in step S2, the Helmholtz equation is:

[0030]

[0031] in ρ is the Hamiltonian operator; ρ is the resistivity, with units of Ωm.

[0032] Preferably, in step S2, the discrete expression is obtained through the following steps:

[0033] S2-1. Based on the Helmholtz equation, the curl operation matrix C is approximately calculated using the integral of the finite volume method (FVM). e2f C f2e and average matrix A c2e ;

[0034] In the grid cells excluding the casing, the average matrix A c2e Used to average the following relationship:

[0035]

[0036] in,

[0037] ρ is the average resistivity that is kept constant within the control volume;

[0038] ρ i It is the resistivity of the i-th cell;

[0039] S i It is the cross-sectional area of ​​the i-th cell portion within the control volume;

[0040] S2-2. From equation (4), the average resistivity ρ defined on the grid edge can be obtained as the area-weighted average of four adjacent resistivities;

[0041] S2-3. Obtain the curl operation matrix C e2f C f2e and average matrix A c2e Then, the discrete form of the Helmholtz equation was further derived.

[0042] (C f2e C e2f +iωμdiag(A c2e m))E=-iωμJ s (1);

[0043] Where m is the vector of the resistivity model.

[0044] Preferably, in step S3, the situation is divided into two cases based on whether the source device is completely inside the sleeve: one is that the source is completely inside the sleeve, and the other is that the source is not inside or near the sleeve.

[0045] Preferably, if the emission source is entirely inside the sleeve, an additional term is added to modify the resistivity averaging process, where the average matrix A in the mesh cells containing the sleeve is... c2e This is used to express the following relationship, while the average discretization process of the remaining cell grids remains unchanged.

[0046]

[0047] Where S5 is the cross-sectional area of ​​the steel pipe, typically 10. -2 In the above formula, S does not need to be changed; the simulation of thin steel sleeves is carried out by adding the product of the sleeve cross-sectional area and its inherent resistivity.

[0048] Preferably, if the measurement is not performed inside or near the sleeve, the sleeve is considered equivalent to a long conductor with comprehensive conductivity properties, which are the product of the cross-sectional area and the inherent resistivity; that is, the resistivity of the unit grid is replaced by...

[0049]

[0050] For any well path that does not necessarily coincide with the grid edge, the conductivity of the casing segment is redistributed to the eight nearest edges using orthogonal decomposition and trilinear interpolation. The contributions of all adjacent casing segments are summed to form the total edge conductivity for each edge, and the average matrix A is reconstructed. c2e And rewritten as

[0051] Preferably, in step S4, the simulation equation is:

[0052]

[0053] Preferably, the specific steps of step S4 are as follows:

[0054] The average matrix A obtained in step S3 c2e Instead of finding the average matrix of the Helmholtz equation in step 2, each electric dipole source is orthogonally decomposed into x, y and z components, and then each component is trilinearly distributed to the eight nearest grid edges in the same direction;

[0055] The additional current density from the source is expressed as J s The non-zero elements in the source term J s Together with iωμ, they form the right-hand side term (rhs), which is represented in symbolic form as follows:

[0056] AE = rhs (7);

[0057] The direct solver PARDISO is used to solve for E in the above equation; finally, the required data at any position is calculated as follows:

[0058] d = QE (8);

[0059] Where Q is a sparse matrix used to perform trilinear interpolation from the grid edges to the observation location.

[0060] Preferably, in step S5, the step of constructing the objective function is as follows:

[0061] S5-1: Logarithmically process the inverted resistivity ρ, i.e.

[0062] m = -ln(ρ) (9);

[0063] S5-2: Use the Tikhonov regularization method, i.e.

[0064]

[0065] in,

[0066] m is an unknown model vector;

[0067] and These are the data fitting term and the model normalization term, respectively. Together, they ensure that a suitable optimal solution is found.

[0068] λ is the Tikhonov regularization parameter, used to balance the weights of the data fitting term and the model regularization term in the objective function.

[0069] S5-3: Use the normal equation of the Gauss-Newton method to obtain the minimum value of the objective function, that is, the update vector δm from the k-th to the (k+1)-th iteration. k The solution formula is as follows, which gives the equation of the normal vector.

[0070] H k δm k =-g k (11);

[0071] Where g k and H k The distribution of the first-order partial derivatives of the objective function and the Hessian matrix is ​​obtained from the following equation.

[0072]

[0073] In the formula The sensitivity matrix is ​​represented by F, which is the first-order partial derivative of the forward operator with respect to the model parameters. The first-order partial derivative of the objective function is g. k The Hessian matrix H can be obtained directly. k Further calculations are needed after solving the sensitivity matrix.

[0074] Preferably, in step S6, the steps for solving the Hessian matrix in step S5 are as follows:

[0075] S6-1: Using the quasi-forward modeling method, for the case of a single frequency and field source, the electromagnetic field at the measurement point can be expressed as:

[0076] d(m)=ψ(E,m)(14);

[0077] Where ψ is the difference function between the electric field and the signal at the measuring point, the above equation is solved and transposed;

[0078] S6-2: The sensitivity matrix J can be solved by the following formula.

[0079] J T =G T A(m) -1 F T +Q T (15);

[0080] In the formula

[0081]

[0082] Substitute the obtained sensitivity matrix into the formula for solving the Hessian matrix in step 5 to solve for the Hessian matrix.

[0083] Preferably, the specific steps of step S7 are as follows:

[0084] S7-1: Update step size α k Start with a unit step size and test whether the condition is met using the following formula. If the condition is not met, gradually decrease the step size until the condition is met:

[0085]

[0086] S7-2: In each Gauss-Newton iteration, the model vector update direction and update step size are calculated using the normal equation vector and the above formula, respectively. After obtaining the model update vector, the model vector used in the next iteration is...

[0087] m k+1 =m k +α k δm k (18);

[0088] Preferably, step S8, achieving the inversion requirement is mainly based on three aspects: reaching the maximum number of iterations, the data fitting difference reaching the target value, and the model update vector being less than a predetermined value.

[0089] Furthermore, this invention provides the application of the well-to-surface differential electromagnetic three-dimensional inversion method considering the role of the well casing in the oilfield exploration and development process.

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

[0091] This invention provides a multi-directional, multi-attribute, five-dimensional seismic fracture prediction method based on ensemble learning. It fuses extracted multi-directional and multi-angle seismic attributes to form a fused attribute that can reflect the characteristics of the fracture to the greatest extent. It uses density-based noise spatial clustering (DBSCAN) technology to fuse the fracture volume, which can more comprehensively characterize reservoir fractures. Attached Figure Description

[0092] Figure 1 The flowchart of a multi-directional, multi-attribute, five-dimensional earthquake crack prediction method based on ensemble learning provided by this invention.

[0093] Figure 2 This is an improved finite volume method in a Cartesian coordinate grid. The left side shows a schematic diagram of the control volume defined around the edge in the y direction; the right side shows a schematic diagram of the calculation of average conductivity.

[0094] Figure 3This is a diagram showing the location of the measuring points in Example 1.

[0095] Figure 4 This is a schematic diagram of the test model for Example 1.

[0096] Figure 5 The image shows a comparison of the forward modeling response at the center of the ground projection of the anomalous body, considering and not considering the resistivity of the casing in Example 1.

[0097] Figure 6 This is a slice of the inversion results considering the resistivity of the bushing in Example 1. Detailed Implementation

[0098] To make the technical means, creative features, achieved objectives, and effects of this invention readily understandable, the invention is further illustrated below with specific embodiments. However, these embodiments are merely preferred embodiments and not all embodiments. Other embodiments obtained by those skilled in the art based on the embodiments described herein without creative effort are all within the scope of protection of this invention. It is worth noting that the raw materials used in this invention are all common commercially available products, and their sources are not specifically limited. The technical and scientific terms used in the embodiments have the meanings commonly understood by those skilled in the art to which this invention pertains.

[0099] Example 1

[0100] A well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of well casing is illustrated in the flowchart below. Figure 1 As shown, the specific steps include:

[0101] Step S1: Read the model file and data file of the inversion region. The model file and data file should at least include: mesh division, steel sleeve length, resistivity, cross-sectional area, location, transmitting and receiving devices, transmission frequency, and measured data.

[0102] Step S2: Neglecting displacement current, we rewrite the equations as Helmholtz equations based on Maxwell's equations in the frequency domain and Ohm's law:

[0103]

[0104] The curl operation matrix C is then approximated by integration using the finite volume method (FVM). e2f C f2e and average matrix A c2e In the grid cells excluding the casing, the average matrix A c2e Used to average the following relationship:

[0105]

[0106] in,

[0107] ρ is the average resistivity that is kept constant within the control volume;

[0108] ρ i It is the resistivity of the i-th cell;

[0109] S i It is the cross-sectional area of ​​the i-th cell portion within the control volume.

[0110] From the above equation, we can derive the average resistivity ρ defined at the grid edge, which is the area-weighted average of four adjacent resistivities. See the schematic diagram of the pseudo-finite volume method discretization. Figure 2 .

[0111] From equation (4), we can derive the average resistivity ρ defined on the grid edge, which is the area-weighted average of four adjacent resistivities;

[0112] Obtain the curl operation matrix C e2f C f2e and average matrix A c2e Then, the discrete form of the Helmholtz equation was further derived.

[0113] (C f2e C e2f +iωμdiag(A c2e m))E=-iωμJ s (1);

[0114] Where m is the vector of the resistivity model.

[0115] Step S3: The discussion is based on whether the emission source is entirely inside the sleeve. If the emission source is completely inside the sleeve, an additional term is added to modify the resistivity averaging process. In the mesh cells containing the sleeve, the average matrix A... c2e This is used to express the following relationship, while the average discretization process of the remaining cell grids remains unchanged.

[0116]

[0117] Where S5 is the cross-sectional area of ​​the steel pipe, typically 10. -2 In the above formula, S does not need to be changed (square meters).

[0118] The simulation of thin steel sleeves is performed by adding the product of the sleeve's cross-sectional area and its inherent resistivity.

[0119] If the emission source is not measured inside or near the sleeve, the sleeve can be considered equivalent to a long conductor with comprehensive conductivity properties, which are the product of the cross-sectional area and the inherent resistivity. That is, the resistivity of the cell grid is replaced by...

[0120]

[0121] For any well path that does not necessarily coincide with the grid edges, the conductivity of the casing segment is redistributed to the eight nearest edges using orthogonal decomposition and trilinear interpolation. The contributions of all adjacent casing segments are summed to form the total edge conductivity for each edge, and the average matrix A is reconstructed. c2e And rewritten as

[0122] Step S4: Use the average matrix A obtained in step 3 c2e Instead of obtaining the average matrix of the Helmholtz equation in step 2, the simulation equation is:

[0123]

[0124] Each electric dipole source is orthogonally decomposed into x, y, and z components, and each component is then trilinearly distributed to the eight nearest grid edges in the same direction. The additional current density from the source is expressed as J s The non-zero elements in the source term J s Together with iωμ, they form the right-hand side term (rhs), which is represented in symbolic form as follows:

[0125] AE = rhs (7)

[0126] The direct solver PARDISO is used to solve for E in the above equation. Finally, the required data at any location is calculated as follows:

[0127] d = QE (8)

[0128] Where Q is a sparse matrix used to perform trilinear interpolation from the grid edges to the observation location.

[0129] Step S5: Based on the measured data obtained in Step 1 and the simulation equations in Step 4, construct the objective function. First, perform logarithmic processing on the retrieved resistivity ρ, i.e.

[0130] m=-ln(ρ) (9)

[0131] Then, the Tikhonov regularization method was used, namely...

[0132]

[0133] Where m is an unknown model vector. and These are the data fitting term and the model regularization term, respectively, which work together to ensure a suitable optimal solution is found. λ is the Tikhonov regularization parameter, used to balance the weights of the data fitting term and the model regularization term in the objective function.

[0134] Subsequently, the minimum value of the objective function is obtained using the normal equations of the Gauss-Newton method, which is the update vector δm from the k-th to the (k+1-th)-th iterations. k The solution formula is as follows, which gives the equation of the normal vector.

[0135] H k δm k =-g k (11)

[0136] Where g k and H k The distribution of the first-order partial derivatives of the objective function and the Hessian matrix is ​​obtained from the following equation.

[0137]

[0138] In the formula The sensitivity matrix is ​​represented by F, which is the first-order partial derivative of the forward operator with respect to the model parameters. The first-order partial derivative of the objective function is g. k The Hessian matrix H can be obtained directly. k Further calculations are needed after solving the sensitivity matrix.

[0139] Step S6: Differentiate the simulation equations from Step 4 using the quasi-forward modeling method to obtain the sensitivity matrix. Using the quasi-forward modeling method, for a single frequency and field source, the electromagnetic field at the measurement point can be expressed as:

[0140] d(m)=ψ(E,m) (14)

[0141] Where ψ is the difference function between the electric field and the signal at the measuring point, the sensitivity matrix J can be obtained by solving the above equation and transposing it.

[0142] J T =G T A(m) -1 F T +Q T (15)

[0143] In the formula

[0144]

[0145] Substitute the obtained sensitivity matrix into the formula for solving the Hessian matrix in step 5 to solve for the Hessian matrix. Then, use the PCG method and the normal equations in step 5 to obtain the model vector update direction for each Gauss-Newton iteration.

[0146] Step S7: Obtain the model update vector. Use the inverse tracking method to perform an inaccurate line search on the update direction to determine the update step size α. k Update step size α kStart with a unit step size and test whether the condition is met using the following formula. If the condition is not met, gradually decrease the step size until the condition is met:

[0147]

[0148] In each Gauss-Newton iteration, the model vector update direction and update step size are calculated using the normal equation vector and the above formula, respectively. After obtaining the model update vector, the model vector used in the next iteration is...

[0149] m k+1 =m k +α k δm k (18)

[0150] Step S8: Repeat steps 5, 6, and 7 until the requirements are met. The process of exiting the Gauss-Newton iteration process upon meeting the inversion requirements is primarily based on three criteria: reaching the maximum number of iterations, the data fit difference reaching the target value, and the model update vector being less than a predetermined value.

[0151] To test the application effect of the present invention, such as Figure 3 A three-dimensional model was designed, featuring a 5Ω·m low-resistivity body within a uniform 100Ω·m half-space. The total mesh size is 44×44×68. The emission source is a 300m differential source, centered at (0,0,-6350)m, as shown in the diagram. Figure 4 As shown.

[0152] Two forward modeling methods were used respectively: (1) considering the effect of bushing resistivity; (2) not considering the effect of bushing resistivity. Figure 5 The forward modeling responses of the anomaly within the ground projection (-175, -175) m using the two forward modeling methods described above were compared. It can be seen that the casing resistivity has a significant impact on the forward modeling results. Therefore, the influence of casing resistivity should be considered when performing well-to-ground differential electromagnetic three-dimensional forward modeling.

[0153] Furthermore, the response obtained by forward modeling considering casing resistivity is used as a data file for well-to-surface differential electromagnetic three-dimensional inversion considering the effect of the well casing. Figure 6 The inversion results are presented in a format that allows for visualization.

[0154] Finally, it should be noted that the above content is only used to illustrate the technical solution of the present invention, and is not intended to limit the scope of protection of the present invention. Simple modifications or equivalent substitutions made by those skilled in the art to the technical solution of the present invention do not depart from the essence and scope of the technical solution of the present invention.

Claims

1. A well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of well casing, characterized in that, Includes the following steps: Step S1: Read the model file and data file of the inverted region to obtain the measured data; Step S2: Based on the Maxwell equations in the frequency domain and Ohm's law, rewrite the Maxwell equations into Helmholtz equations, and use the pseudo-finite volume method to discretize the spatial domain to obtain the discrete form; Step S3: Reconstruct the average matrix based on whether the source device is completely located inside the sleeve. Step S4: Use the average matrix obtained in step S3 Substituting into the discrete equation described in step S2, we obtain the simulation equation, and use the mimicry finite volume method to analyze and obtain the overall coefficient matrix A for each transmission frequency; Step S5: Based on the measured data obtained in Step S1 and the simulation equation obtained in Step S4, construct the objective function and obtain the normal vector equation H. k δm k =-g k First-order partial derivative g k and Hessian matrix H k ; Step S6: Obtain the sensitivity matrix using the quasi-forward modeling method, then solve for the Hessian matrix obtained in step S5, and finally solve the normal vector equation using the PCG method to calculate the update direction δm of the model vector in each iteration. k ; Step S7: Obtain the model update vector: Use the inverse tracking method to perform an inaccurate line search on the update direction to determine the update step size α. k The model vector used in the next iteration is m. k+1 =m k +α k δm k ; Step S8: Repeat steps S5, S6, and S7 until the requirements are met.

2. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of well casing as described in claim 1, characterized in that, In step S1, the model file and data file contain the following information: mesh generation, steel sleeve length, resistivity, cross-sectional area, location, transmitting and receiving devices, transmission frequency, and measured data.

3. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of well casing as described in claim 1, characterized in that, In step S2, the discrete expression is: (C f2e C e2f +iωμdiag(A c2e m))E=-iωμJ s (1); Among them, C f2e With C e2f These are the discrete matrices of the twisting operator from the center of the mesh to the edge and the discrete matrices of the twisting operator from the edge to the center of the mesh, respectively. ω is the angular frequency at the current desired frequency; A c2e This is the average matrix from the center of the grid cell to the edge of the grid. m is the conductivity vector discretized onto the grid cells, and its unit is S / m; E is the electric field vector, and its unit is V / m; J s For external field sources, the unit is A / m. 2 ; μ is the magnetic permeability, with units of H / m; i is the imaginary unit.

4. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of the well casing as described in claim 1, characterized in that, In step S2, the Helmholtz equation is: in ρ is the Hamiltonian operator; ρ is the resistivity, with units of Ωm.

5. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of the well casing according to claim 1, characterized in that, In step S2, the discrete expression is obtained through the following steps: S2-1. Based on the Helmholtz equation, the curl operation matrix C is approximately calculated using the integral of the finite volume method (FVM). e2f C f2e and average matrix A c2e ; In the grid cells excluding the casing, the average matrix A c2e Used to average the following relationship: in, ρ is the average resistivity that is kept constant within the control volume; ρ i It is the resistivity of the i-th cell; S i It is the cross-sectional area of ​​the i-th cell portion within the control volume; S2-2. From equation (4), the average resistivity ρ defined on the grid edge can be obtained as the area-weighted average of four adjacent resistivities; S2-3. Obtain the curl operation matrix C e2f C f2e and average matrix A c2e Then, the discrete form of the Helmholtz equation was further derived. (C f2e C e2f +iωμdiag(A c2e m))E=-iωμJ s (1); Where m is the vector of the resistivity model.

6. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of the well casing according to claim 1, characterized in that, Step S2 is performed while neglecting displacement current.

7. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of the well casing according to claim 1, characterized in that, In step S3, the situation is divided into two cases based on whether the source device is completely inside the sleeve: one is that the source is completely inside the sleeve, and the other is that the source is not inside or near the sleeve.

8. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of the well casing according to claim 7, characterized in that, If the emission source is entirely inside the sleeve, an additional term is added to modify the resistivity averaging process. In the mesh cells containing the sleeve, the average matrix A is... c2e This is used to express the following relationship, while the average discretization process of the remaining cell grids remains unchanged. Where S5 is the cross-sectional area of ​​the steel pipe, typically 10. -2 In the above formula, S does not need to be changed; the simulation of thin steel sleeves is carried out by adding the product of the sleeve cross-sectional area and its inherent resistivity.

9. A well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of well casing, as described in claim 7, is characterized in that, If the emission source is not measured inside or near the sleeve, the sleeve can be considered equivalent to a long conductor with comprehensive conductivity properties, which are the product of the cross-sectional area and the inherent resistivity. That is, the resistivity of the cell grid is replaced by... For any well path that does not necessarily coincide with the grid edge, the conductivity of the casing segment is redistributed to the eight nearest edges using orthogonal decomposition and trilinear interpolation. The contributions of all adjacent casing segments are summed to form the total edge conductivity for each edge, and the average matrix A is reconstructed. c2e And rewritten as 10. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of the well casing according to claim 1, characterized in that, In step S4, the simulation equation is:

11. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of well casing according to claim 1, characterized in that, The specific steps of step S4 are as follows: The average matrix A obtained in step S3 c2e Instead of finding the average matrix of the Helmholtz equation in step 2, each electric dipole source is orthogonally decomposed into x, y and z components, and then each component is trilinearly distributed to the eight nearest grid edges in the same direction; The additional current density from the source is expressed as J s The non-zero elements in the source term J s Together with iωμ, they form the right-hand side term (rhs), which is represented in symbolic form as follows: AE = rhs (7); The direct solver PARDISO is used to solve for E in the above equation; finally, the required data at any position is calculated as follows: d = QE (8); Where Q is a sparse matrix used to perform trilinear interpolation from the grid edges to the observation location.

12. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of well casing according to claim 1, characterized in that, In step S5, the steps for constructing the objective function are as follows: S5-1: Logarithmically process the inverted resistivity ρ, i.e. m = -ln(ρ) (9); S5-2: Employ the Tikhonov regularization method, i.e. in, m is an unknown model vector; and These are the data fitting term and the model normalization term, which work together to ensure that a suitable optimal solution is found. λ is the Tikhonov regularization parameter, used to balance the weights of the data fitting term and the model regularization term in the objective function; S5-3: Use the normal equation of the Gauss-Newton method to obtain the minimum value of the objective function, that is, the update vector δm from the k-th to the (k+1)-th iteration. k The solution formula is as follows, which gives the equation of the normal vector. H k nm k = -g k (11); Where g k and H k The distribution of the first-order partial derivatives of the objective function and the Hessian matrix is ​​obtained from the following equation. In the formula F represents the sensitivity matrix, which is the first-order partial derivative of the forward operator with respect to the model parameters; where the first-order partial derivative of the objective function g k The Hessian matrix H can be obtained directly. k Further calculations are needed after solving the sensitivity matrix.

13. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of the well casing according to claim 1, characterized in that, In step S6, the steps for solving the Hessian matrix from step S5 are as follows: S6-1: Using the quasi-forward modeling method, for the case of a single frequency and field source, the electromagnetic field at the measurement point can be expressed as: d(m)=ψ(E,m)(13); Where ψ is the difference function between the electric field and the signal at the measuring point, the above equation is solved and transposed; S6-2: The sensitivity matrix J can be solved by the following formula. J T =G T A(m) -1 F T +Q T (14); In the formula Substitute the obtained sensitivity matrix into the formula for solving the Hessian matrix in step 5 to solve for the Hessian matrix.

14. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of well casing according to claim 1, characterized in that, The specific steps of step S7 are as follows: S7-1: Update step size α k Start with a unit step size and test whether the condition is met using the following formula. If the condition is not met, gradually decrease the step size until the condition is met: S7-2: In each Gauss-Newton iteration, the model vector update direction and update step size are calculated using the normal equation vector and the above formula, respectively. After obtaining the model update vector, the model vector used in the next iteration is... m k+1 =m k +a k δm k (17)。 15. The well-to-surface differential electromagnetic three-dimensional inversion method considering the effect of the well casing according to claim 1, characterized in that, Step S8, achieving the inversion requirement, is mainly based on three aspects: reaching the maximum number of iterations, the data fitting difference reaching the target value, and the model update vector being less than a predetermined value.

16. The application of the well-to-surface differential electromagnetic three-dimensional inversion method considering the role of well casing as described in any one of claims 1-15 in the oilfield exploration and development process.

Images