Three-dimensional induced polarization transient electromagnetic forward method

Abstract

The application discloses a three-dimensional induced polarization transient electromagnetic forward method. The method takes the grounded long wire source transient electromagnetic detection as an object, establishes a calculation region containing a transmitting source, an observation point, a background medium and a three-dimensional abnormal body, inputs medium conductivity, magnetic permeability, polarization rate, time constant, frequency correlation coefficient, transmitting current waveform, observation point position and time sampling information, solves by spectral element space discretization, Caputo fractional order induced polarization constitutive relation and SOE memory variable recursion coupling, and obtains three-dimensional transient electromagnetic electric field response containing induced polarization effect. The application adopts non-uniform hexahedral spectral element discretization and high-order vector edge basis function, can better meet the electric field tangential continuity condition at the medium interface, and improves the spatial discretization precision of the electric field response in the complex three-dimensional model; the polarization memory effect is directly introduced in the time domain, which can reduce the influence of frequency sampling and transformation algorithm on the result; the SOE recursive memory variable can reduce the storage and operation burden of historical items in long time channel calculation.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration technology, and in particular to a three-dimensional induced transient electromagnetic forward modeling method. Background Technology

[0002] Transient electromagnetic methods have been widely used to detect underground electrical anomalies such as coal mine water hazards, goaf water accumulation, water-rich fractured zones, and concealed water-conducting structures. Among these methods, grounded long conductor sources can induce strong conduction currents, and the horizontal electric field component is highly sensitive to low-resistivity bodies, polarized bodies, and electrical interfaces. Therefore, three-dimensional transient electromagnetic electric field response forward modeling is an important foundation for data interpretation and inversion imaging. Existing three-dimensional transient electromagnetic forward modeling methods mostly employ finite difference methods, finite volume methods, or low-order finite element methods, which struggle to balance adaptability to complex media, spatial discretization accuracy, and degree-of-freedom scale. Furthermore, induced polarization effects are prevalent in underground rocks, minerals, clay minerals, and water-bearing fractured media, leading to late-stage tailing, slower decay, negative responses, or sign reversals in the response. Considering only resistivity differences can easily affect the determination of anomaly properties.

[0003] Therefore, there is an urgent need for a three-dimensional forward modeling method for induced transient electromagnetic field response that can take into account high-order spatial discretization accuracy, direct introduction of time-domain induced polarization effect, and efficient calculation of fractional-order history terms. Summary of the Invention

[0004] This invention solves the technical problems existing in the prior art by providing a three-dimensional induced transient electromagnetic forward modeling method.

[0005] This invention provides a three-dimensional induced transient electromagnetic forward modeling method, comprising: Formula construction (1); In equation (1), For vector differential operators, Let μ be the electric field intensity at time t for position vector r = (x, y, z), and μ be the magnetic permeability. Let be the magnetic field strength at time t for position vector r=(x,y,z). Let r be the dielectric current density at time t for position vector r = (x, y, z). Let r be the applied source current density at time t for the position vector r=(x,y,z); Taking the curl of the first equation in equation (1) and substituting the second equation into it, we obtain the vector diffusion equation with the electric field as the fundamental unknown: (2); Construct higher-order interpolation basis functions, and further construct vector edge basis functions that satisfy the tangential continuity condition. The electric field in the e-th element is approximated as follows: (3); In equation (3), Let r be the electric field of the element at position r and time t within the e-th spectral element. This represents the total number of degrees of freedom within the element. For higher-order vector edge basis functions, For reference unit coordinates, F is the unknown electric field on the corresponding edge degree of freedom. e For the first Geometric mapping of spectral units, Let r be the dielectric current density at position r and time t within the e-th spectral element. For the unknown quantity of dielectric current density on the corresponding edge degree of freedom; Within each spectral element, equation (3) is used as a finite-dimensional expansion of the electric field, and the same set of vector boundary basis functions is used. By performing a weighted integration of equation (2) using the Galerkin weight function, and constructing the element weighting residual corresponding to the i-th weight function of the e-th element based on equation (2), it can be written as: (4); In equation (4), This is the weighted residual corresponding to the i-th Galerkin weight function in the e-th cell. Let e ​​be the integration region of the e-th spectral element in physical space. Let i be the basis function of the vector edge within the e-th unit. Let e ​​be the electric field within the e-th unit. Let e ​​be the current density within the e-th cell. For the applied source current density, This represents the volume integral region in physical space. Substituting equation (3) into equation (4), we perform integral by parts on the curl-curl term and, combined with the tangential continuity condition and boundary conditions of the vector edge basis functions, obtain the curl inner product form of the basis functions; at the same time, we transform the two main left-hand terms in equation (2) into: (5); Substituting equation (5) and the source term into the weighted residual (4), we obtain the element matrix equation: (6); In equation (6), This is the right-hand term of the unit formed by the change in the emission source current; Assemble equation (6) of each unit according to the common edge degrees of freedom, and use Delta line source terms distributed along the conductor path for the grounded long conductor source to obtain the semi-discrete spectral element matrix equation: (7); In equation (7), The global stiffness matrix is ​​formed by assembling the stiffness matrices of each element. Let $\mathbf{ ... The global mass matrix is ​​formed by assembling the mass matrices of each unit. For global dielectric current density, This is the equivalent right-hand term formed by the change in current from a grounded long conductor source; Ohm's law in the frequency domain Combining and rearranging with the Cole-Cole model, we get: (8); In equation (8), Zero-frequency conductivity, The frequency domain electric field intensity at angular frequency ω It is the product of the imaginary unit i and the angular frequency ω. The frequency correlation coefficient, It is a time constant. The frequency domain dielectric current density at angular frequency ω Polarizability; Taking the inverse Fourier transform of equation (8) yields the time-domain Caputo fractional induced polarization constitutive relation: (9); In equation (9), is a Caputo fractional-order time derivative operator of order c with a lower limit of integration of 0; Constructing the Caputo fractional derivative: (10); In equation (10), Let be a time function whose fractional derivative is to be found. For the Gamma function, For historical time variables, Let f be the first derivative of the function f with respect to time at the historical time η; In the time step At this point, the integration interval of equation (10) is split into local terms of the current step and historical terms of the past: (11); Construct the kernel function that appears in the historical integral within a given time interval: (12); In equation (12), The number of exponent terms. For the first The weights of each index term, For the first One index node; Based on equation (12), for any historical function Construct the first One memory variable: (13); In equation (13), For the l-th SOE memory variable corresponding to the history function f in t k The value at time, ; Substituting equations (10) to (13) into equation (9) and rearranging, the current density at the current moment can be written in a recursive form with respect to the current electric field, the electric field at the previous moment, the current density at the previous moment, and the SOE memory variable: (14); In equation (14), and They are respectively The current density vector and electric field vector at time t; , , , , They are respectively from , , , , The coefficients are jointly determined by the SOE node weights; and These are the SOE memory variables corresponding to the electric field and current density, respectively; For any time function ,exist Derivative of construction time: (15); In equation (15), Let f be a time function at time t k The discrete value at time is obtained by substituting equation (14) into equation (7) and then using equation (15) to discretize. You can get the first A system of global linear equations at time steps: (16); In equation (16), For the first The system matrix at each time step includes the spectral element curl stiffness matrix, mass matrix, induced electrostatic constitutive coefficients, and time discretization coefficients; e k Let $\mathbf{ ... The right-hand term includes the emission source term, the field value terms of the first two steps, the current density term of the previous step, and the SOE memory variable term. After solving equation (16), the global electric field edge degrees of freedom e at the current time step are obtained. k Then, the current density j is updated using equation (14). k And update the SOE memory variables corresponding to the electric field and current density according to equation (13). and The updated e k、 j k and memory variables and The historical information used for the next time step continues to participate in the recursion until all time sampling points have been calculated. Based on the spectral element where the observation point is located and its local coordinates, interpolation is performed using the same set of higher-order vector basis functions to obtain the electric field response at the observation point: (17); In equation (17), For the first Each observation point is at The three-dimensional electric field response at time t, For the spatial location of the observation point, This is the spectral element interpolation operator corresponding to the observation point.

[0006] Specifically, in equation (5), the first equation is obtained by integrating the curl-curl term by parts and combining the tangential continuity of the vector edge basis functions and the corresponding boundary conditions, while the second equation is obtained directly by expanding the projection of the current density in the same edge basis function space.

[0007] Specifically, it also includes: The response of the horizontal electric field component Ex along the direction of the grounded long conductor source at the observation point in the background model, which does not contain target anomalies and maintains consistency in emission source, observation point, mesh partitioning, and time sampling, as a function of time, is calculated. ; The response of the same electric field component Ex at the same observation point as the time-varying point was calculated in the model after incorporating the three-dimensional polarization anomaly. ; Subtracting the values ​​point by point from the same observation point and the same time channel yields the difference response. ; Through the difference response This provides a basis for anomaly localization, boundary identification, and subsequent multi-parameter inversion.

[0008] One or more technical solutions provided in this invention have at least the following technical effects or advantages: Taking transient electromagnetic detection of a grounded long-conducting source as the object, a computational domain is established, including the source, observation point, background medium, and three-dimensional anomaly. Inputs include the medium's conductivity, permeability, polarizability, time constant, frequency correlation coefficient, emission current waveform, observation point location, and time sampling information. Through spectral spatial discretization, Caputo fractional-order induced polarization constitutive relations, and SOE memory variable recursive coupling, the three-dimensional transient electromagnetic field response with induced polarization effect is obtained. This invention employs non-uniform hexahedral spectral element discretization and high-order vector edge basis functions, which better satisfies the tangential continuity condition of the electric field at the medium interface and improves the spatial discretization accuracy of the electric field response in complex three-dimensional models. This invention directly introduces the polarization memory effect in the time domain, reducing the impact of frequency sampling and transformation algorithms on the results. Furthermore, SOE recursive memory variables can reduce the storage and computational burden of historical terms in long-term calculations. More importantly, this invention highlights the response characteristics of the horizontal electric field component to polarizability, time constant, frequency correlation coefficient, resistivity difference, and the location of three-dimensional anomalies under the condition of grounded long conductor source. It can be used to analyze phenomena such as late tailing, slowed decay, negative response, and sign reversal caused by induced polarization, providing a high-precision forward modeling basis for the interpretation of transient electromagnetic electric field data containing induced polarization and subsequent three-dimensional multi-parameter inversion. Attached Figure Description

[0009] Figure 1 A flowchart of the three-dimensional induced transient electromagnetic forward modeling method provided in the embodiments of the present invention; Figure 2 This is a schematic diagram of the three-dimensional vector edge basis functions in an embodiment of the present invention, representing the vector edge basis functions in the x, y, and z directions, respectively, used to illustrate the construction of basis functions for discrete spectral space; Figure 3 This is a schematic diagram of a layered dielectric grounding long conductor source model and a source-receiver device in an embodiment of the present invention, used to illustrate the three-dimensional calculation area, the arrangement of the transmitting source and the observation point; Figure 4 This is a comparison diagram of the response and relative error of the layered medium model Ex in the embodiments of the present invention, used to illustrate the accuracy comparison between the spectral element method, the finite element method and the analytical solution; Figure 5 This is a comparison chart of the SOE approximation curve and the exact weakly singular kernel function in the embodiments of the present invention, used to illustrate the approximation effect of exponential and approximation on fractional kernel functions; Figure 6 This is a verification error diagram of SOE time input in an embodiment of the present invention, used to illustrate the error control of SOE approximation under a preset target error; Figure 7 This is a diagram of the SOE exponent nodes and weight distribution in an embodiment of the present invention, used to illustrate how the exponent terms cover changes in the early and late memory kernels on a logarithmic scale; Figure 8 This is a comparison diagram of the Ex response of the model containing the induced polarization effect and the one-dimensional semi-analytical solution in the embodiments of the present invention, used to illustrate the reliability of the forward modeling results containing the induced polarization effect; Figure 9 This is a schematic diagram of a three-dimensional dual-polarization anomaly model in an embodiment of the present invention, used to illustrate the spatial arrangement of the emission source, observation array, and three-dimensional anomaly. Figure 10 The graphs show the electric field difference response ΔEx curves under different time channels in this embodiment of the invention, used to illustrate the time-domain anomalous response changes caused by the three-dimensional polarization anomaly. Figure 11 These are cross-sectional and plan views of the electric field difference response ΔEx of the three-dimensional polarization anomaly in this embodiment of the invention, used to illustrate the influence of the spatial location, boundary and electrical differences of the anomaly on the electric field anomaly response. The meanings of the English and Chinese labels in the attached figure are as follows: Source indicates the source of transmission, Receiver indicates the receiving point, SEM indicates the spectral element method, FEM indicates the finite element method, Analytical indicates the analytical solution, SOE approximation indicates the exponential and approximate solution, Relative error indicates the relative error, and Nodes / Weights indicates the exponential nodes and weights. Detailed Implementation

[0010] like Figure 1 As shown, the technical solution in this embodiment of the invention addresses the technical problems existing in the prior art, and the overall approach is as follows: First, forward modeling parameters and model information are input, including the computational domain, dielectric parameters, Cole-Cole induced polarization parameters, grounded long conductor source, observation point location, and time sampling information. Then, based on the quasi-static Maxwell equations, a three-dimensional transient electromagnetic field control equation is established with the electric field as the fundamental unknown. On this basis, this invention performs calculations from two aspects: spatial discretization and induced polarization time-domain modeling.

[0011] In terms of spatial discretization, the computational domain is partitioned into a non-uniform hexahedral spectral element mesh. GLL nodes, higher-order vector edge basis functions, and corresponding spectral element discretization schemes are constructed and assembled to form the spectral element stiffness matrix, mass matrix, and source term vector. For induced polarization time-domain modeling, the time-domain fractional constitutive relation is derived from the Cole-Cole frequency-domain induced polarization model, and the Caputo fractional derivative is used to describe the time memory effect of the polarization medium. Furthermore, the fractional derivative is decomposed into local and historical terms, and the SOE exponent and approximation are used to accelerate the processing of the historical term. A finite number of memory variables are constructed and recursively updated.

[0012] Subsequently, the discrete matrix of the spectral space is coupled with the fractional-order induced polarization (SOE) memory term to construct a linear system of equations for the three-dimensional transient electromagnetic field response containing the induced polarization effect. During the calculation, the initial electric field or initial time-time response is first obtained, and then a variable-step-size second-order back-Eulerian scheme is used for time advancement, combined with a direct solver to solve the electric field degrees of freedom step-by-step. While waiting for the last time sampling point, the right-hand side, SOE memory variables, and electric field response are continuously updated; after reaching the last time sampling point, the calculated global electric field degrees of freedom are interpolated, and the Ex, Ey, and Ez components at the observation points are extracted, with a focus on analyzing the Ex component parallel to the direction of the grounded long conductor source. Finally, the three-dimensional forward modeling results are output, and the accuracy, stability, and applicability of the method can be evaluated by combining layered model verification, SOE error verification, and three-dimensional polarization anomaly response analysis.

[0013] To better understand the above technical solutions, the following will provide a detailed explanation of the technical solutions in conjunction with the accompanying drawings and specific implementation methods.

[0014] The three-dimensional induced transient electromagnetic forward modeling method provided in this embodiment of the invention includes: Formula construction (1); In equation (1), For vector differential operators, Let μ be the electric field intensity at time t for position vector r = (x, y, z), and μ be the magnetic permeability. Let be the magnetic field strength at time t for position vector r=(x,y,z). Let r be the dielectric current density at time t for position vector r = (x, y, z). Let r be the applied source current density at time t for the position vector r=(x,y,z); Taking the curl of the first equation in equation (1) and substituting the second equation into it, we obtain the vector diffusion equation with the electric field as the fundamental unknown: (2); Equation (2) shows that the response of the underground electric field is controlled by two parts: one is the spatial curl term, which reflects the diffusion of the electromagnetic field in the three-dimensional medium; the other is the current density time variation term, which reflects the influence of the medium's conductivity and polarization effect on the transient process. The subsequent spectral element discretization, induced constitutive model and time progression are all based on Equation (2).

[0015] Construct higher-order interpolation basis functions, and further construct vector edge basis functions that satisfy the tangential continuity condition. The electric field in the e-th element is approximated as follows: (3); In equation (3), Let r be the electric field of the element at position r and time t within the e-th spectral element. This represents the total number of degrees of freedom within the element. For higher-order vector edge basis functions, For reference unit coordinates, , F is the unknown electric field on the corresponding edge degree of freedom. e For the first The geometric mapping of spectral elements allows non-uniform hexahedral elements in complex three-dimensional regions to be uniformly transformed into a numerical integration problem over a reference domain. Let r be the dielectric current density at position r and time t within the e-th spectral element. For the unknown quantity of dielectric current density on the corresponding edge degree of freedom; Within each spectral element, equation (3) is used as a finite-dimensional expansion of the electric field, and the same set of vector boundary basis functions is used. By performing a weighted integration of equation (2) using the Galerkin weight function, and constructing the element weighting residual corresponding to the i-th weight function of the e-th element based on equation (2), it can be written as: (4); In equation (4), This is the weighted residual corresponding to the i-th Galerkin weight function in the e-th cell. Let e ​​be the integration region of the e-th spectral element in physical space. Let i be the basis function of the vector edge within the e-th unit. Let e ​​be the electric field within the e-th unit. Let e ​​be the current density within the e-th cell. For the applied source current density, This represents the volume integral region in physical space. Substituting equation (3) into equation (4), we perform integral by parts on the curl-curl term and, combined with the tangential continuity condition of the vector edge basis functions and the boundary conditions, obtain the curl inner product form of the basis functions. Simultaneously, since the current density and electric field use the same set of vector edge basis functions for expansion, the current density time derivative term can be written as the product of the basis function inner product and the time derivative of the edge degrees of freedom. Therefore, the two main left-hand terms in equation (2) are transformed into: (5); In equation (5), the first equation is obtained by integrating the curl-curl term by parts and combining the tangential continuity of the vector edge basis functions and the corresponding boundary conditions, while the second equation is obtained directly by expanding the projection of the current density in the same edge basis function space.

[0016] Substituting equation (5) and the source term into the weighted residual (4), we obtain the element matrix equation: (6); In equation (6), The element right-hand side is formed by the change in the source current, and the electric field expansion coefficients are used to form the element stiffness matrix through curl integration. The current density time term is multiplied by the basis functions to form the unit mass matrix. .

[0017] Assemble equation (6) of each unit according to the common edge degrees of freedom, and use Delta line source terms distributed along the conductor path for the grounded long conductor source to obtain the semi-discrete spectral element matrix equation: (7); In equation (7), The global stiffness matrix is ​​formed by assembling the stiffness matrices of each element. Let $\mathbf{ ... The global mass matrix is ​​formed by assembling the mass matrices of each unit. For global dielectric current density, This is the equivalent right-hand term formed by the change in current from a grounded long conductor source; Ohm's law in the frequency domain Combining and rearranging with the Cole-Cole model, we get: (8); In equation (8), Zero-frequency conductivity, The frequency domain electric field intensity at angular frequency ω It is the product of the imaginary unit i and the angular frequency ω. The frequency correlation coefficient, , It is a time constant. The frequency domain dielectric current density at angular frequency ω Polarizability; Taking the inverse Fourier transform of equation (8) yields the time-domain Caputo fractional induced polarization constitutive relation: (9); In equation (9), The Caputo fractional time derivative operator with order c and a lower limit of integration of 0 is given by equation (9). Equation (9) shows that the dielectric current density at the current moment depends not only on the current electric field, but also on the changes in the historical electric field and historical current density. This equation is the basis for the direct introduction of the induced polarization effect in the time domain in this invention.

[0018] Constructing the Caputo fractional derivative: (10); In equation (10), Let be a time function whose fractional derivative is to be found. For the Gamma function, For historical time variables, Let f be the first derivative of the function f with respect to time at the historical time η; In the time step At this point, the integration interval of equation (10) is split into local terms of the current step and historical terms of the past: (11); In equation (11), the first term is only related to the current time step. The first term relates to the field values ​​within the time frame and can be directly discretized using piecewise linear interpolation. The second term contains information from all past time steps, which is the main reason for the increased storage and computational complexity of the direct fractional-order algorithm. Therefore, this invention only uses SOE acceleration for historical terms, while local terms are still processed using direct discretization.

[0019] Construct the kernel function that appears in the historical integral within a given time interval: (12); In equation (12), The number of exponent terms. For the first The weights of each index term, For the first The significance of equation (12) lies in transforming the historical kernel of the power function, which is difficult to directly recursively derive, into a superposition of several exponential functions, while the exponential function has a natural time recursive property.

[0020] Based on equation (12), for any historical function Construct the first One memory variable: (13); In equation (13), For the l-th SOE memory variable corresponding to the history function f in t k The value at time l is used to characterize the historical contribution accumulated by the history function under the l-th exponential kernel. Equation (13) shows that the historical contribution of the current time step can be obtained by combining the memory variables of the previous time step and the local integral within the current time step, thus eliminating the need to store the field values ​​at all historical moments. In actual calculations, the contribution within the current time step can be obtained by combining the memory variables of the previous time step and the local integral within the current time step. By using linear interpolation, equation (13) is transformed into one containing only... and The explicit update form.

[0021] Substituting equations (10) to (13) into equation (9) and rearranging, the current density at the current moment can be written in a recursive form with respect to the current electric field, the electric field at the previous moment, the current density at the previous moment, and the SOE memory variable: (14); In equation (14), and They are respectively The current density vector and electric field vector at time t; , , , , They are respectively from , , , , The coefficients are jointly determined by the SOE node weights; and These are the SOE memory variables corresponding to the electric field and current density, respectively; this formula preserves the recursive structure necessary for algorithm implementation.

[0022] This invention employs variable-step long-time sampling and utilizes an implicit second-order back-out Euler scheme to approximate the time derivative. For any time function... ,exist Derivative of construction time: (15); In equation (15), Let f be a time function at time t k The discrete values ​​at time points, Equation (15) can adapt to non-uniform time steps, making the early response have denser time sampling and the late response have higher computational efficiency.

[0023] Substitute equation (14) into equation (7), and then use equation (15) to discretize. You can get the first A system of global linear equations at time steps: (16); In equation (16), For the first The system matrix at each time step includes the spectral element curl stiffness matrix, mass matrix, induced electrostatic constitutive coefficients, and time discretization coefficients; e k Let $\mathbf{ ... The right-hand term includes the emission source term, the field value terms of the first two steps, the current density term of the previous step, and the SOE memory variable term. After solving equation (16), the global electric field edge degrees of freedom e at the current time step are obtained. k Then, the current density j is updated using equation (14). k And update the SOE memory variables corresponding to the electric field and current density according to equation (13). and The updated e k、 j k and memory variables and The historical information used for the next time step continues to participate in the recursion until all time sampling points have been calculated. Based on the spectral element where the observation point is located and its local coordinates, interpolation is performed using the same set of higher-order vector basis functions to obtain the electric field response at the observation point: (17); In equation (17), For the first Each observation point is at The three-dimensional electric field response at time t, For the spatial location of the observation point, This is the spectral interpolation operator corresponding to the observation point. In interpreting data from grounded long conductor sources, the horizontal electric field component parallel to the transmitting conductor can be extracted as a key focus. .

[0024] To highlight the local response differences caused by the three-dimensional polarization anomaly, two forward modeling calculations were performed using the background model and the anomaly model. Specifically, this also included: The response of the horizontal electric field component Ex along the direction of the grounded long conductor source at the observation point in the background model, which does not contain target anomalies and maintains consistency in emission source, observation point, mesh partitioning, and time sampling, as a function of time, is calculated. ; The response of the same electric field component Ex at the same observation point as the time-varying point was calculated in the model after incorporating the three-dimensional polarization anomaly. ; Subtracting the values ​​point by point from the same observation point and the same time channel yields the difference response. The difference response can weaken the dominant components of the background field and the source field, and enhance the local disturbances caused by differences in polarizability, resistivity and the spatial location of the anomalous body. It can further form survey curves, time trace profiles or planar contour maps.

[0025] Through the difference response This provides a basis for anomaly localization, boundary identification, and subsequent multi-parameter inversion.

[0026] Example 1: Accuracy Verification of Spectral Spatial Discretization and Layered Medium In Example 1, a layered dielectric model of a grounded long conductor source is first established. The transmitting source and receiving point are arranged on the ground surface, and the underground dielectric is divided according to its layered resistivity structure. (Appendix) Figure 2 A schematic diagram of three-dimensional vector edge basis functions is given, representing the construction methods of edge basis functions in the x, y, and z directions, respectively, to illustrate the physical meaning of the higher-order vector edge basis functions in equation (3). This type of basis function uses the edge degrees of freedom as unknowns, making it suitable for describing the tangential continuity condition in a three-dimensional electromagnetic field.

[0027] Appendix Figure 3A model of a long conductor grounded in a layered medium and a schematic diagram of the source-receiver device are presented. This diagram corresponds to the source-transmitter settings in equations (2) and (5). The source is located on the Earth's surface, and the observation point is located at a specified location near the source. The underground contains a layered electrical interface. This model can be used to verify the spatial discretization accuracy of the spectral element method in regular layered media.

[0028] Appendix Figure 4 The results of comparing the response and relative error of the layered medium model Ex are presented. Comparing the third-order spectral element results of this invention with the low-order finite element results and analytical solutions, it can be seen that the spectral element results and analytical solutions agree well, and lower full-time relative errors can be obtained under similar or fewer degrees of freedom conditions. This embodiment illustrates that the spectral element spatial discretization described in equations (3) to (5) can improve the accuracy of three-dimensional transient electromagnetic field response calculations—degree-of-freedom utilization efficiency.

[0029] Example 2: SOE approximation and verification of response with induced polarization effect In Example 2, a uniformly polarized half-space or a polarized low-resistivity layered dielectric model was selected to verify the reliability of the Caputo fractional induced polarization term and SOE memory variable recursion described in equations (8) to (12). Example parameters include a background resistivity of 100 Ω·m, a low-resistivity layer resistivity of 20 Ω·m, a polarizability of m = 0.1, a time constant of τ = 1.0, and a frequency correlation coefficient of c = 0.5.

[0030] Appendix Figure 5 A comparison between the SOE approximation curve and the exact weakly singular kernel function is given. This figure corresponds to Equation (10), which illustrates that the power function kernel can be approximated by a finite number of exponential functions. The two curves highly overlap in a double logarithmic coordinate system, indicating that the SOE approximation can preserve the main variation characteristics of the fractional kernel function.

[0031] Appendix Figure 6 A time-input verification error plot for SOE is presented. This plot illustrates that, given the target error, the overall relative error of the SOE approximation is within a controllable range. Since SOE transforms historical convolution into a finite number of memory variables in Equation (13), the number of exponential terms is much smaller than the total number of time steps, thus reducing the storage and computational load for long-duration forward modeling.

[0032] Appendix Figure 7 A distribution diagram of the SOE exponent nodes and weights is provided to illustrate the coverage of the exponent nodes and weights on a logarithmic scale. This diagram corresponds to s in equation (12). l and ω l This indicates that the SOE term can take into account both early rapid changes and late long memory effects.

[0033] Appendix Figure 8The results comparing the Ex response of the model with the induced polarization effect with the one-dimensional semi-analytical solution are presented. The calculation results can characterize the mid-to-late-stage decay, zero-crossing, and negative response characteristics caused by the polarization effect. Except for the local amplification of the relative error near the zero-crossing point due to the reference response amplitude being close to zero, most of the time channel errors remain within a reasonable range, indicating that the direct time-domain induced polarization spectral element forward modeling method described in equations (9) to (16) can effectively simulate the time-domain electric field response with the induced polarization effect.

[0034] Example 3: Electric Field Response Analysis of a Three-Dimensional Polarization Anomaly In Example 3, a grounded long conductor source model containing a three-dimensional polarization anomaly is established. An array of observation points is set up on the ground surface. The Ex components of the background model and the model containing the anomaly are calculated respectively, and the difference response ΔEx is obtained. (Appendix) Figure 9 A schematic diagram of a three-dimensional dual-polarization anomaly model is provided. This diagram illustrates the spatial arrangement of the emission source, observation array, and three-dimensional anomaly, and also demonstrates that the present invention can handle problems in non-layered three-dimensional polarized media.

[0035] Appendix Figure 10 The response curves of ΔEx in the direction of the measured line under different time traces are presented. As time progresses, the amplitude and shape of the anomalous response change, with local negative responses and curve curvature gradually emerging. These results indicate that the horizontal electric field component of a grounded long conductor source is highly sensitive to the location of the three-dimensional polarimeter, polarization parameters, and resistivity differences.

[0036] Appendix Figure 11 A profile of the electric field difference response ΔEx of a three-dimensional polarization anomaly is presented. The figure shows significant local perturbations, negative responses, and boundary gradient changes near the anomaly's projection area, indicating that this invention can be used not only for accuracy verification of layered models but also for electric field response analysis of complex three-dimensional polarization anomalies. Identification of the profile response provides a forward modeling basis for subsequent three-dimensional polarization anomaly localization and parameter inversion.

[0037] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0038] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0039] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0040] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0041] Any aspects of this invention not described in detail in the embodiments are well-known techniques to those skilled in the art. Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of this invention and not to limit it. Although this invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of this invention without departing from the spirit and scope of this invention, and all such modifications and substitutions should be covered within the scope of the claims of this invention.

Claims

1. A three-dimensional induced transient electromagnetic forward modeling method, characterized in that, include: Formula construction (1); In equation (1), For vector differential operators, Let μ be the electric field intensity at time t for position vector r = (x, y, z), and μ be the magnetic permeability. Let be the magnetic field strength at time t for position vector r=(x,y,z). Let r be the dielectric current density at time t for position vector r = (x, y, z). Let r be the applied source current density at time t for the position vector r=(x,y,z); Taking the curl of the first equation in equation (1) and substituting the second equation into it, we obtain the vector diffusion equation with the electric field as the fundamental unknown: （2）； Construct higher-order interpolation basis functions, and further construct vector edge basis functions that satisfy the tangential continuity condition. The electric field in the e-th element is approximated as follows: (3); In equation (3), Let r be the electric field of the element at position r and time t within the e-th spectral element. This represents the total number of degrees of freedom within the element. For higher-order vector edge basis functions, For reference unit coordinates, F is the unknown electric field on the corresponding edge degree of freedom. e For the first Geometric mapping of spectral units, Let r be the dielectric current density at position r and time t within the e-th spectral element. For the unknown quantity of dielectric current density on the corresponding edge degree of freedom; Within each spectral element, equation (3) is used as a finite-dimensional expansion of the electric field, and the same set of vector boundary basis functions is used. By performing a weighted integration of equation (2) using the Galerkin weight function, and constructing the element weighting residual corresponding to the i-th weight function of the e-th element based on equation (2), it can be written as: (4); In equation (4), This is the weighted residual corresponding to the i-th Galerkin weight function in the e-th cell. Let e ​​be the integration region of the e-th spectral element in physical space. Let i be the basis function of the vector edge within the e-th unit. Let e ​​be the electric field within the e-th unit. Let e ​​be the current density within the e-th cell. For the applied source current density, This represents the volume integral region in physical space. Substituting equation (3) into equation (4), we perform integral by parts on the curl-curl term and, combined with the tangential continuity condition and boundary conditions of the vector edge basis functions, obtain the curl inner product form of the basis functions; at the same time, we transform the two main left-hand terms in equation (2) into: （5）； Substituting equation (5) and the source term into the weighted residual (4), we obtain the element matrix equation: (6); In equation (6), This is the right-hand term of the unit formed by the change in the emission source current; Assemble equation (6) of each unit according to the common edge degrees of freedom, and use Delta line source terms distributed along the conductor path for the grounded long conductor source to obtain the semi-discrete spectral element matrix equation: (7); In equation (7), The global stiffness matrix is ​​formed by assembling the stiffness matrices of each element. Let $\mathbf{ ... The global mass matrix is ​​formed by assembling the mass matrices of each unit. For global dielectric current density, This is the equivalent right-hand term formed by the change in current from a grounded long conductor source; Ohm's law in the frequency domain Combining and rearranging with the Cole-Cole model, we get: (8); In equation (8), Zero-frequency conductivity, The frequency domain electric field intensity at angular frequency ω It is the product of the imaginary unit i and the angular frequency ω. The frequency correlation coefficient, It is a time constant. The frequency domain dielectric current density at angular frequency ω Polarizability; Taking the inverse Fourier transform of equation (8) yields the time-domain Caputo fractional induced polarization constitutive relation: (9); In equation (9), is a Caputo fractional-order time derivative operator of order c with a lower limit of integration of 0; Constructing the Caputo fractional derivative: (10); In equation (10), Let be a time function whose fractional derivative is to be found. For the Gamma function, For historical time variables, Let f be the first derivative of the function f with respect to time at the historical time η; In the time step At this point, the integration interval of equation (10) is split into local terms of the current step and historical terms of the past: （11）； Construct the kernel function that appears in the historical integral within a given time interval: (12); In equation (12), The number of exponent terms. For the first The weights of each index term, For the first One index node; Based on equation (12), for any historical function Construct the first One memory variable: (13); In equation (13), For the l-th SOE memory variable corresponding to the history function f in t k The value at time, ; Substituting equations (10) to (13) into equation (9) and rearranging, the current density at the current moment can be written in a recursive form with respect to the current electric field, the electric field at the previous moment, the current density at the previous moment, and the SOE memory variable: (14); In equation (14), and They are respectively The current density vector and electric field vector at time t; , , , , They are respectively from , , , , The coefficients are jointly determined by the SOE node weights; and These are the SOE memory variables corresponding to the electric field and current density, respectively; For any time function ,exist Derivative of construction time: (15); In equation (15), Let f be a time function at time t k The discrete value at time is obtained by substituting equation (14) into equation (7) and then using equation (15) to discretize. You can get the first A system of global linear equations at time steps: (16); In equation (16), For the first The system matrix at each time step includes the spectral element curl stiffness matrix, mass matrix, induced electrostatic constitutive coefficients, and time discretization coefficients; e k Let $\mathbf{ ... The right-hand term includes the emission source term, the field value terms of the first two steps, the current density term of the previous step, and the SOE memory variable term. After solving equation (16), the global electric field edge degrees of freedom e at the current time step are obtained. k Then, the current density j is updated using equation (14). k And update the SOE memory variables corresponding to the electric field and current density according to equation (13). and The updated e k、 j k and memory variables and The historical information used for the next time step continues to participate in the recursion until all time sampling points have been calculated. Based on the spectral element where the observation point is located and its local coordinates, interpolation is performed using the same set of higher-order vector basis functions to obtain the electric field response at the observation point: (17); In equation (17), For the first Each observation point is at The three-dimensional electric field response at time t, For the spatial location of the observation point, This is the spectral element interpolation operator corresponding to the observation point.

2. The three-dimensional induced transient electromagnetic forward modeling method as described in claim 1, characterized in that, In equation (5), the first equation is obtained by integrating the curl-curl term by parts and combining the tangential continuity of the vector edge basis functions and the corresponding boundary conditions. The second equation is obtained directly by expanding the projection of the current density in the same edge basis function space.

3. The three-dimensional induced transient electromagnetic forward modeling method as described in claim 1 or 2, characterized in that, Also includes: The response of the horizontal electric field component Ex along the direction of the grounded long conductor source at the observation point in the background model, which does not contain target anomalies and maintains consistency in emission source, observation point, mesh partitioning, and time sampling, as a function of time, is calculated. ; The response of the same electric field component Ex at the same observation point as the time-varying point was calculated in the model after incorporating the three-dimensional polarization anomaly. ; Subtracting the values ​​point by point from the same observation point and the same time channel yields the difference response. ; Through the difference response This provides a basis for anomaly localization, boundary identification, and subsequent multi-parameter inversion.

Images