Transient electromagnetic forward modeling method and device with induced polarization effect and storage medium

Abstract

The invention provides a transient electromagnetic forward modeling method and device containing an induced polarization effect and a storage medium, and relates to the technical field of electrical digital data processing. The method comprises the following steps: constructing a geologic model; applying an excitation field to the geologic model; determining a cut-off time T1 and a cut-off frequency F1 within a time period from T0 to T2 after the excitation field disappears; obtaining time domain electromagnetic data in a time period from T0 to T1 and frequency domain electromagnetic data in a frequency section from Fmin to F1 based on the first mathematical model and a Maxwell equation set; based on the frequency domain electromagnetic data, the time domain electromagnetic data and a second mathematical model, obtaining response electromagnetic data containing an induced polarization effect at any position in the geologic model within a time period from T0 to T2; according to the transient electromagnetic forward modeling method containing the induced polarization effect, the problem that in the prior art, due to the fact that calculation is difficult in the frequency domain, transient electromagnetic three-dimensional rapid forward modeling containing the induced polarization effect is difficult to achieve can be solved.

Description

Technical Field

[0001] The present application relates to the field of electrical digital data processing, and in particular to a transient electromagnetic forward modeling method, device and storage medium including induced polarization effect. Background Art

[0002] Transient Electromagnetic Method (TEM) is a geophysical exploration method based on the principle of electromagnetic induction. When using TEM to explore a target area, the transmitter emits a primary electromagnetic field to the target area. The primary electromagnetic field excites the medium in the target area to generate a secondary electromagnetic field, which contains relevant geoelectric information of the medium. The secondary electromagnetic field is received by the receiver, and analyzed and processed to obtain the geoelectric information of the medium, and then predict whether the target area contains metal deposits and other issues of interest.

[0003] The research on transient electromagnetic method includes two aspects: forward modeling and inversion. Forward modeling refers to simulating the response electromagnetic data (secondary electromagnetic field) of a geological model under a certain excitation field (primary electromagnetic field) through a certain forward modeling method given the geological model of the target area. Inversion refers to obtaining the parameters of the geological model based on the response electromagnetic data.

[0004] At present, the 3D forward modeling of transient electromagnetics with induced polarization (IP) effect (also known as "IP effect") is generally completed in the frequency domain and then converted to the time domain through Fourier transform. The problem is that the frequency domain response needs to be calculated point by point within a wide frequency band, which takes a lot of time; and as the frequency increases, the size of the coefficient matrix increases, and the required computer memory and calculation time increase dramatically. This makes it difficult to achieve fast 3D forward modeling of transient electromagnetics with IP effect. Summary of the invention

[0005] The present application provides a method, device and storage medium for transient electromagnetic forward modeling including induced polarization effect, so as to solve the problem in the prior art that three-dimensional rapid forward modeling of transient electromagnetic including induced polarization effect is difficult to achieve due to difficulty in frequency domain calculation.

[0006] In order to achieve the above-mentioned object, the first aspect of the present application provides a transient electromagnetic forward modeling method including induced polarization effect, the method comprising the following steps:

[0007] Constructing a geological model, wherein the geological model includes resistivity parameters and polarizability parameters;

[0008] applying an excitation field to the geological model;

[0009] In the time period from T0 to T2 after the excitation field disappears, the cutoff time T1 and the cutoff frequency F1 are determined; wherein T0<T1<T2, F min <F1<F max , F min is the minimum frequency in the period from T0 to T2, F max is the maximum frequency in the time period from T0 to T2;

[0010] Based on the first mathematical model and Maxwell equations, the time domain electromagnetic data and F in the time period from T0 to T1 are obtained. min to frequency domain electromagnetic data within the F1 frequency range, wherein both the time domain electromagnetic data and the frequency domain electromagnetic data contain induced polarization effect; and

[0011] Based on the frequency domain electromagnetic data, the time domain electromagnetic data and the second mathematical model, response electromagnetic data containing induced polarization effect at any position in the geological model within the time period from T0 to T2 is obtained.

[0012] In a possible implementation manner, the second mathematical model uses orthogonal polynomials as basis functions.

[0013] In a possible implementation manner, the orthogonal polynomial is a Laguerre orthogonal polynomial, and the second mathematical model includes a time domain response function and a frequency domain response function;

[0014] The time domain response function is: Among them, t represents time, a n are the polynomial coefficients, is the normalized Laguerre orthogonal polynomial, l t is the time domain scale factor, l t For limiting the domain of the time domain response function to the Within the scope of action,

[0015] The frequency domain response function is: Where f represents frequency, a n are the polynomial coefficients, Φ n By the Perform Fourier transform to obtain l f is the frequency domain scale factor, l f The frequency domain response function is limited to the domain of the frequency domain response function in the Φ n Within the scope of action, j is the imaginary unit.

[0016] In a possible implementation, the polynomial coefficient a n Obtained by the following steps:

[0017] Acquire P time domain sampling points from the time domain electromagnetic data;

[0018] Based on the time domain sampling points and the time domain response function, a discrete time domain response polynomial is obtained.

[0019]

[0020] Acquire Q frequency domain sampling points from the frequency domain response data;

[0021] Based on the frequency domain sampling points and the frequency domain response function, a discrete frequency domain response polynomial is obtained.

[0022]

[0023] Based on the time domain response polynomial h(t p ) and the frequency domain response polynomial H(f q ) Calculate the polynomial coefficient a n .

[0024] In a possible implementation, the time domain response polynomial h(t p ) and the frequency domain response polynomial H(f q ) and optimizing the time domain scaling factor l t and the frequency domain scaling factor l f :

[0025] Step A: setting the initial expansion item number N0 and the initial scale factor, wherein the initial scale factor includes an initial time domain scale factor and an initial frequency domain scale factor;

[0026] Step B: Calculate min{e(l based on the initial expansion item number N0 and the initial scale factor t )}, obtain the optimal scale factor and the first polynomial coefficient, wherein the optimal scale factor includes the optimal time domain scale factor and the optimal frequency domain scale factor, wherein

[0027]

[0028] Step C: Based on the polynomial coefficients of M to N expansion terms, the expansion error and the matching error are calculated, where 0<M<N, and the expansion error is: The matching error is the e(l t );

[0029] Step D: Compare the expansion error with a preset expansion accuracy threshold, and compare the matching error with a preset matching accuracy threshold:

[0030] When the expansion error is greater than the expansion accuracy threshold, or the matching error is greater than the matching accuracy threshold, N0+I is used as the new initial number of expansion items, and the optimal scale factor is used as the new initial scale factor, and steps A to D are repeated, where I is a positive integer; when the expansion error is less than or equal to the expansion accuracy threshold, and the matching error is less than or equal to the matching accuracy threshold, the current number of expansion items and scale factor are the optimized number of expansion items and scale factor, and the optimization procedure ends.

[0031] In a possible implementation, the first mathematical model is a Cole-Cole model: Where j is the imaginary unit, ω is the angular frequency; σ * (ω) is the resistivity; σ0 is the zero-frequency conductivity, which corresponds to the conductivity of the earth when there is no induced polarization effect; m is the charging rate, which characterizes the strength of the induced polarization effect of the earth, 0≤m≤1; τ is the time constant, which characterizes the slowness of the induced polarization effect of the earth; c is the frequency correlation coefficient, which characterizes the change of the phase of the induced polarization effect with frequency, 0≤c≤1.

[0032] In a possible implementation manner, the time domain electromagnetic data and F in the time period from T0 to T1 are obtained based on the first mathematical model and the Maxwell equations. min Frequency domain electromagnetic data in the F1 frequency range, including:

[0033] Gridding the geological model to obtain grid units; and

[0034] Based on the first mathematical model and the Maxwell equations, the time domain electromagnetic data and the frequency domain electromagnetic data of each grid unit are obtained.

[0035] In a possible implementation, obtaining the time domain electromagnetic data and the frequency domain electromagnetic data of each grid unit based on the first mathematical model and the Maxwell equations includes:

[0036] Determine the resistivity of the nodes of each grid cell;

[0037] Calculating the resistivity of any point in each grid unit according to the resistivity of the node and the interpolation function to obtain a grid unit with parameter assignment; and

[0038] Based on the first mathematical model, the Maxwell equations and the grid units with parameter assignments, the time domain electromagnetic data and the frequency domain electromagnetic data of each grid unit are obtained.

[0039] A second aspect of the present application provides a transient electromagnetic forward modeling device including induced polarization effect, the device comprising:

[0040] A modeling module, which is used to construct a geological model, wherein the geological model includes resistivity parameters and polarizability parameters;

[0041] an excitation module, which is used to apply an excitation field to the geological model;

[0042] The determination module is used to determine the cutoff time T1 and the cutoff frequency F1 in the time period from T0 to T2 after the excitation field disappears; wherein T0<T1<T2, F min <F1<F max , F min is the minimum frequency in the period from T0 to T2, F max is the maximum frequency in the time period from T0 to T2;

[0043] The first acquisition module is used to obtain the time domain electromagnetic data and F in the time period from T0 to T1 based on the first mathematical model and the Maxwell equations. min to frequency domain electromagnetic data within the F1 frequency range, wherein both the time domain electromagnetic data and the frequency domain electromagnetic data contain induced polarization effect; and

[0044] The second acquisition module is used to obtain the response electromagnetic data containing induced polarization effect at any position in the geological model within the time period from T0 to T2 based on the frequency domain electromagnetic data, the time domain electromagnetic data and the second mathematical model.

[0045] A third aspect of the present application provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the transient electromagnetic forward modeling method including induced polarization effect provided in any of the above-mentioned embodiments.

[0046] Compared with the prior art, the beneficial effect of the transient electromagnetic forward modeling method with induced polarization effect provided by the present application is that the method provides a way of joint derivation in the time and frequency domains, which effectively avoids the shortcomings of separate derivation in the time and frequency domains. The method obtains time domain electromagnetic data in the high-frequency region and frequency domain electromagnetic data in the low-frequency region, avoiding the calculation of frequency domain data in the high-frequency region, thereby simplifying the calculation and greatly improving the calculation efficiency, solving the problem in the prior art that the three-dimensional rapid forward modeling of transient electromagnetic with induced polarization effect is difficult to achieve due to the difficulty in frequency domain calculation. The method performs a joint derivation in the time and frequency domains, obtains the response electromagnetic data of the entire time and frequency domains, and thus realizes the rapid forward modeling of transient electromagnetic with induced polarization effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0048] Figure 1 A flowchart of a transient electromagnetic forward modeling method including induced polarization effect according to an embodiment of the present application;

[0049] Figure 2 A flowchart for implementing a method for optimizing polynomial coefficients and scale factors in a transient electromagnetic forward modeling method including induced polarization effect according to an embodiment of the present application;

[0050] Figure 3 It is a schematic diagram of the structure of a transient electromagnetic forward modeling device including induced polarization effect according to an embodiment of the present application;

[0051] Figure 4 It is a schematic diagram of the structural composition of the forward modeling equipment according to the implementation mode of the present application. DETAILED DESCRIPTION

[0052] In the following description, specific details such as specific system structures, technologies, etc. are provided for the purpose of illustration rather than limitation, so as to provide a thorough understanding of the embodiments of the present application. However, it should be clear to those skilled in the art that the present application may also be implemented in other embodiments without these specific details. In other cases, detailed descriptions of well-known systems, devices, circuits, and methods are omitted to prevent unnecessary details from obstructing the description of the present application.

[0053] TEM has the advantages of being easy to implement and sensitive to good conductors, and is one of the important means of mineral resource exploration. TEM has been widely used in the fields of metal mineral exploration, geological surveys, etc. The presence of polarization characteristics in metal minerals causes the propagation of electromagnetic waves in the strata to produce dispersion effects (i.e., IP effect, induced polarization effect). The induced polarization effect can cause the transient electromagnetic response to be distorted or even cause a sign reversal, destroying the normal attenuation law of the transient electromagnetic response as the earth's electrical properties change, thereby bringing difficulties and interference to the processing of transient electromagnetic data, and even causing erroneous interpretations. However, the appearance of the induced polarization effect often means the existence of mineral resources, and it is necessary to study electromagnetic data containing the induced polarization effect.

[0054] For transient electromagnetic data containing induced polarization effect, complex resistivity model is often used for modeling. Complex resistivity is a variable related to time. When using the complex resistivity model for forward modeling, the complex resistivity and coefficient matrix both change with time, and the coefficient matrix equation needs to be repeatedly solved for each iteration, which cannot take advantage of the fast forward modeling speed in the time domain. Therefore, although TEM is a time domain method, at present, the three-dimensional forward modeling of transient electromagnetic containing induced polarization effect is generally performed in the frequency domain, and then converted to the time domain through Fourier transform. However, as mentioned in the background technology, in the frequency domain, due to computational difficulties, the three-dimensional rapid forward modeling of transient electromagnetic containing induced polarization effect is difficult to achieve.

[0055] To solve this problem, the present application provides a forward modeling method to realize three-dimensional rapid forward modeling of transient electromagnetic including induced polarization effect. To make the purpose, technical solution and advantages of the present application clearer, the following will be described in conjunction with the accompanying drawings and through specific implementation methods.

[0056] In one embodiment, see Figure 1 , the forward modeling method provided in this application comprises the following steps:

[0057] S100: constructing a geological model, where the geological model includes resistivity parameters and polarizability parameters;

[0058] S200: applying an excitation field to the geological model;

[0059] S300: Determine the cutoff time T1 and the cutoff frequency F1 in the time period from T0 to T2 after the excitation field disappears; wherein T0<T1<T2, F min <F1<F max , F min is the minimum frequency in the period from T0 to T2, F max is the maximum frequency in the time period from T0 to T2;

[0060] S400: Based on the first mathematical model and Maxwell equations, obtain the time domain electromagnetic data and F in the time period from T0 to T1 min Frequency domain electromagnetic data within the F1 frequency range, wherein both the time domain electromagnetic data and the frequency domain electromagnetic data contain induced polarization effects; and

[0061] S500: Based on the frequency domain electromagnetic data, the time domain electromagnetic data and the second mathematical model, the response electromagnetic data including induced polarization effect at any position in the geological model within the time period from T0 to T2 is obtained.

[0062] In S100, the construction of the geological model can be realized with the help of computer software (such as COMSOL, GID). The geological model includes not only geometric information, but also resistivity and polarizability at any position. Resistivity corresponds to the normal response of transient electromagnetics, and polarizability corresponds to the induced electrical response.

[0063] In S200, the application of the excitation field can also be achieved with the help of computer software. In practical applications, the excitation field is generated by a transmitting device. There are many forms of transmitting devices, and the common ones are electric dipoles, magnetic dipoles, grounded power lines, and ungrounded return lines. Correspondingly, the transmitting device of the software model also has many forms, and this application does not limit this; that is, the forward modeling method provided in this application can be applied to the excitation field generated by any form of transmitting device.

[0064] In S300, the time period to be studied is determined, that is, T0 to T2, which can be called a sampling period. Forward modeling is to derive the response electromagnetic data of any point in the geological model during the sampling period, such as the magnitude and direction of the response magnetic field. Exemplarily, the timing can be started after the excitation field disappears (turns off), that is, T0=0, and T2 can be tens to thousands of milliseconds (ms). For example, the turn-off time can be 1μs, and the timing starts after the turn-off, T0=0, T2=200ms.

[0065] S300 also determines the cutoff time, i.e., T1. Determining T1 can divide the sampling period under study into a high-frequency period and a low-frequency period. The frequency of the response magnetic field decays over time. Before the cutoff time, i.e., the time period from T0 to T1, can be called a high-frequency period and an early period; after the cutoff time, i.e., the time period from T1 to T2, can be called a low-frequency period and a late period. In order to determine T1, in a relatively simple scenario without induced polarization effect, it can be studied which time period is suitable for early research and which time period is suitable for late research, so as to finally determine a reasonable early and late dividing moment, i.e., T1. The research on transient electromagnetic without induced polarization effect is already very mature, and T1 can be determined based on these previous experiences. For example, for different geological models, when the sampling period is from 10ms to 100ms, T1 can take values ​​in the range of 0.01ms to 1ms.

[0066] S300 also determines the cutoff frequency, i.e., F1. Before the cutoff frequency, i.e., the interval from F0 to F1, it can be called the high-frequency period and the early period; after the cutoff frequency, i.e., the interval from F1 to F2, it can be called the low-frequency period and the late period. In order to determine F1, in a relatively simple scenario without induced polarization effect, it is possible to study which frequencies are suitable for high-frequency research and which frequencies are suitable for low-frequency research, so as to finally determine a reasonable high- and low-frequency dividing point, i.e., F1. The research on transient electromagnetics without induced polarization effect is already very mature, and F1 can also be determined based on these previous experiences. For example, for different geological models, when the frequency of the secondary magnetic field decays in the interval [10Hz, 100kHz], the cutoff frequency can take values ​​in the interval [100Hz, 1000Hz].

[0067] In S400, the first mathematical model and the Maxwell equations are used to calculate the time domain electromagnetic data of the high frequency period (early period) and the frequency domain electromagnetic data of the low frequency period (late period). This avoids the calculation of frequency domain data in the high frequency area. In the high frequency area, because the frequency changes dramatically, if you want to describe the frequency change as accurately as possible, you need to calculate a large number of frequency points, which will cause the calculation complexity to increase sharply. In the low frequency area, the frequency changes slowly, and only a small number of frequency points need to be calculated, which simplifies the calculation and greatly improves the calculation efficiency.

[0068] It can be understood that the first mathematical model is generally a complex resistivity model. The complex resistivity model can characterize the induced polarization effect in transient electromagnetic data. In the process of studying the induced polarization effect, scholars have proposed a series of complex resistivity models (such as the Wait model, the Ward and Fraser model, the Madden and Cantwell model, etc.), which are not limited in this embodiment; that is, any complex resistivity model is applicable to this embodiment.

[0069] In S500, based on the second mathematical model and the frequency domain electromagnetic data and time domain electromagnetic data obtained in S400, the response electromagnetic data of the entire time and frequency domains are derived; that is, the response electromagnetic data of any point in the geological model at any time in the sampling period is obtained. Generally, the response electromagnetic data includes full-band frequency domain response electromagnetic data and full-time time domain response electromagnetic data.

[0070] Specifically, based on the second mathematical model and the frequency domain electromagnetic data and the time domain electromagnetic data obtained in S400, the response electromagnetic data in the entire time and frequency domains are derived, including:

[0071] Based on the second mathematical model and frequency domain electromagnetic data, the whole frequency band F min To F max Full-band frequency domain response electromagnetic data within; and

[0072] Based on the second mathematical model and the time-domain electromagnetic data, the full-time time-domain response electromagnetic data in the entire time period from T0 to T2 is derived.

[0073] Optionally, the second mathematical model is any mathematical model that can be used for numerical simulation. Specifically, the frequency domain electromagnetic data and the time domain electromagnetic data obtained in S400 are used to solve the undetermined parameters in the second mathematical model. After the undetermined parameters are determined, the second mathematical model can be used to calculate the response electromagnetic data at any position and any time in the geological model.

[0074] Compared with the prior art, the beneficial effect of the forward modeling method provided in this embodiment is that it provides a method of joint derivation in the time and frequency domains, which effectively avoids the shortcomings of separate derivation in the time and frequency domains. In S400, time domain electromagnetic data is obtained in the high-frequency region, and frequency domain electromagnetic data is obtained in the low-frequency region, avoiding the calculation of frequency domain data in the high-frequency region, thereby simplifying the calculation, greatly improving the calculation efficiency, and solving the problem in the prior art that the three-dimensional rapid forward modeling of transient electromagnetics containing induced polarization effects is difficult to achieve due to the difficulty in frequency domain calculation. In S500, a joint derivation of the time and frequency domains is performed to obtain the response electromagnetic data of the entire time and frequency domains, thereby realizing the rapid forward modeling of transient electromagnetics containing induced polarization effects.

[0075] For a target area, if the response electromagnetic data obtained by forward modeling is consistent with the inversion electromagnetic data, it indicates that the geological model set in the forward modeling is consistent with the actual geological conditions. Then, the target area can be explored and developed according to the geological model. However, in practical applications, the electromagnetic data of a forward modeling generally cannot conform to the inversion electromagnetic data. Therefore, in the inversion process, it is necessary to continuously perform forward modeling (generally, each inversion iteration requires a forward modeling) to make the forward and inversion electromagnetic data fit and conform. It can be seen that forward modeling is the basis of inversion. If you want to obtain the geological conditions of the target area through inversion, you generally need to perform forward modeling many times. If the speed of forward modeling is slow, it will hinder the smooth implementation of inversion and even TEM. The forward modeling method provided by the present application realizes rapid forward modeling, thereby effectively solving this problem.

[0076] In addition to serving the inversion, another significance of forward modeling is that it can theoretically (with the help of computer software) discuss how a geological model responds to the stimulation of a set primary field and what form the electromagnetic data containing induced polarization effect is generated; without the need to make an actual geological model and use an actual transmitting device to generate a primary field. Forward modeling provides convenience for theoretically grasping the response law of transient electromagnetics.

[0077] In one implementation, the first mathematical model in S400 is a Cole-Cole model: Where j is the imaginary unit, ω is the angular frequency; σ * (ω) is the resistivity; σ0 is the zero-frequency conductivity, which corresponds to the conductivity of the earth when there is no dispersion effect (induced polarization effect); m is the charging rate, which characterizes the strength of the earth's dispersion effect (induced polarization effect), 0≤m≤1; τ is the time constant, which characterizes the slowness of the earth's induced polarization process (induced polarization effect); c is the frequency-dependent coefficient, which characterizes the change of the dispersion effect (induced polarization effect) phase with frequency, 0≤c≤1.

[0078] The Cole-Cole model has resistivity on the left side of the equal sign and parameters characterizing the induced polarization effect on the right side. That is to say, the Cole-Cole model unifies the polarizability of the induced polarization effect into the resistivity. The resistivity calculated by the Cole-Cole model can characterize both the normal response of transient electromagnetics and the response containing the induced polarization effect. Based on the geological model given in S100 and the excitation field given in S200, the Cole-Cole model combined with the Maxwell equations can calculate the electromagnetic data containing the induced polarization effect at any position in the geological model.

[0079] In addition, based on the Cole-Cole model, when the polarizability is determined, the influence of the time constant τ and the frequency correlation coefficient c on the morphological characteristics of the time-varying complex resistivity can be analyzed to find out the time period in which the early complex resistivity changes little and will not have a significant impact on the transient electromagnetic response, thereby reasonably determining the cutoff time T1.

[0080] In one embodiment, in S400, based on the first mathematical model and the Maxwell equations, the time domain electromagnetic data and F in the time period from T0 to T1 are obtained. min The frequency domain electromagnetic data in the F1 frequency range may include the following steps:

[0081] S401: Gridding the geological model in S100 to obtain grid cells; and

[0082] S402: Based on the first mathematical model and the Maxwell equations, obtain time domain electromagnetic data and frequency domain electromagnetic data of each grid unit.

[0083] Specifically, the method provided in this embodiment can be implemented by the finite element method. In addition to the finite element method, in other embodiments, the integral equation method and the finite difference method can also be used to calculate the electromagnetic data of any position of the geological model in the time domain. As for the frequency domain electromagnetic data, the time domain can be converted to the frequency domain by Fourier transform.

[0084] In one implementation, the above S402 obtains the time domain electromagnetic data and the frequency domain electromagnetic data of each grid unit based on the first mathematical model and the Maxwell equations, which may include the following steps:

[0085] S402-1: Determine the resistivity of the nodes of each grid unit;

[0086] S402-2: Calculate the resistivity of any point in each grid unit according to the resistivity of the node and the interpolation function to obtain a grid unit with parameter assignment; and

[0087] S402-3: Based on the first mathematical model, the Maxwell equations and the grid units with parameter assignments, obtain the time domain electromagnetic data and the frequency domain electromagnetic data of each grid unit.

[0088] Compared with the traditional vector finite element method, in this embodiment, the resistivity is assigned to the node rather than the unit, and the interpolation function is used to represent the resistivity in the unit. In this way, the structural characteristics of the geological model do not need to be considered when meshing, and the grid units do not need to be re-divided as the geological model changes, as long as the nodes are in the geological model. In the traditional forward modeling method, when the geological model is continuously optimized and updated during the inversion process, if the geological model changes greatly, it is necessary to re-mesh; and the method provided by this embodiment does not need to re-mesh, thereby avoiding the complex model pre-processing work of the traditional method and improving the forward modeling speed. The method provided by this embodiment has the advantages of the meshless method (however, unlike the meshless method, the method of this embodiment needs to consider the constituent parameters of the unit (for example, the range of the unit, the location of the node, etc.) during the interpolation calculation), reduces the difficulty of meshing, and is suitable for the forward modeling of the resistivity gradient geological model.

[0089] Specifically, the above interpolation function may be a scalar interpolation function. Interpolation may be performed using the following model: in is a scalar interpolation function. In addition, the field at any position in the grid cell can be represented by edges and vector interpolation functions.

[0090] Exemplarily, in this embodiment, when the time-domain electromagnetic data is calculated using the first mathematical model and the Maxwell equations, the variation of the electric field control equation can be derived directly from the time-domain Maxwell equations. The time term can be discretized by backward difference, which is unconditionally stable and has a loose restriction on the time step. Optionally, a larger time step can be selected while ensuring the calculation accuracy. For the entire time zone, a segmented processing strategy can be adopted, with a uniform iterative step size within the segment, and only one stiffness matrix decomposition is required, thereby improving efficiency and shortening calculation time.

[0091] In this example, in terms of loading the excitation field, for each grid unit, the current density can be applied to the edge that coincides with the horizontal component of the electric field. This direct source addition method has better applicability. The waveform of the excitation field can be set arbitrarily in theory, and the method provided in this embodiment can perform full waveform three-dimensional forward modeling of a transient electromagnetic field of arbitrary waveform.

[0092] Exemplarily, in this embodiment, when the frequency domain electromagnetic data is calculated using the first mathematical model and the Maxwell equations, the frequency domain vector finite element model can be derived from the Maxwell equations, the first type of boundary conditions can be added, and the variation of the frequency domain electric field control equation can be established. The grid division, the field at any position in the grid unit, and the resistivity assignment method in the grid unit can all adopt the same strategy as the time domain vector finite element method in the previous example, thereby realizing the three-dimensional forward modeling of the frequency domain vector finite element transient electromagnetic with continuous resistivity variation.

[0093] In this example, the electric field of a horizontal electric dipole in a uniform half-space at any position in space can be derived first, and then the response of the excitation field in the uniform half-space can be loaded into the abnormal field control equation of the electric field as the background field using the abnormal field method. The solution region can be discretized and interpolated using the Witney interpolation basis function to achieve three-dimensional forward modeling of the frequency domain rate.

[0094] In this example, the abnormal field method can be used to load the excitation field. Since the power supply term is not included in the functional of the abnormal electric field, the background field is usually the main part of the total electric field. The background field can be accurately obtained by analytical solution or digital filtering solution; the abnormal field accounts for a small proportion of the total field and can be obtained by solving the finite element equations.

[0095] In one implementation, the second mathematical model in S500 uses orthogonal polynomials as basis functions.

[0096] Under the excitation of electromagnetic pulse signals with limited bandwidth, the pulse width of the time domain electromagnetic response and the bandwidth of the frequency domain electromagnetic response are limited. In other words, the energy of the time domain response and the frequency domain response are always concentrated in a certain time domain and frequency domain range respectively, so the superposition of orthogonal basis functions with relatively concentrated energy can be used to simultaneously represent the time domain response and the frequency domain response.

[0097] That is to say, the transient electromagnetic energy is almost entirely concentrated in the limited range of the time domain and the frequency domain, so that the time domain response can be expanded into an orthogonal polynomial with undetermined coefficients. The orthogonal polynomial is subjected to Fourier transform to obtain a frequency domain response model similar to the time domain response model. Then, a joint extrapolation calculation in the time and frequency domains is performed to determine the undetermined coefficients of the time domain response model and the frequency domain response model from the early time domain electromagnetic data and the late frequency domain electromagnetic data obtained in S400, thereby obtaining the complete response of the entire time domain and frequency domain.

[0098] Optionally, the above orthogonal polynomial may be a Hermite orthogonal polynomial.

[0099] In one implementation, the orthogonal polynomial is a Laguerre orthogonal polynomial, and the second mathematical model in S500 includes a time domain response function and a frequency domain response function;

[0100] The time domain response function is: Among them, t represents time, a n are the polynomial coefficients, is the normalized Laguerre orthogonal polynomial, l t is the time domain scale factor, l t For limiting the definition domain of the time domain response function to the Within the scope of action,

[0101] The frequency domain response function is: Where f represents frequency, a n are the polynomial coefficients, Φ n By the Perform Fourier transform to obtain l f is the frequency domain scale factor, l f The frequency domain response function is limited to the domain of the frequency domain response function in the Φ n Within the scope of action, j is the imaginary unit.

[0102] In this implementation, Laguerre orthogonal polynomials are used as basis functions of the time domain response h(t). Compared with Hermite orthogonal polynomials, Laguerre orthogonal polynomials are more suitable as expanded basis functions of electromagnetic responses. Because the domain of Hermite orthogonal polynomials is [-∞,+∞], a time delay of t0 is required to match the time domain response. If the delay factor is not selected appropriately, it will have a great impact on the results. The domain of Laguerre orthogonal polynomials is [0,+∞], and no time delay is required to match the time domain response; and the lack of delay variables can more effectively speed up parameter optimization.

[0103] The n-th order Laguerre polynomial is:

[0104] The normalized Laguerre orthogonal polynomial is:

[0105] Select the normalized Laguerre orthogonal polynomial As the basis function of the time domain response h(t), then select the appropriate scale factor l t , limiting the definition domain of the time domain response to the scope of the basis function can obtain the above time domain response function:

[0106] Set the basis function The Fourier transform of is Φ(f), and the frequency domain response function can be obtained by performing Fourier transform on both sides of the above time domain response function:

[0107] In one embodiment, the polynomial coefficient a n It can be obtained by following the steps below:

[0108] S501: Acquire P time domain sampling points from the time domain electromagnetic data;

[0109] S502: Based on the time domain sampling points and the time domain response function, obtain a discrete time domain response polynomial

[0110] S503: Acquire Q frequency domain sampling points from the frequency domain response data;

[0111] S504: Obtaining a discrete frequency domain response polynomial based on the frequency domain sampling points and the frequency domain response function

[0112] S505: Based on the time domain response polynomial h(t p ) and the frequency domain response polynomial H(f q ) Calculate the polynomial coefficient a n .

[0113] This embodiment provides a method for calculating the undetermined coefficients in the above-mentioned time-domain and frequency-domain response functions using the time-domain electromagnetic data and the frequency-domain electromagnetic data obtained in S400.

[0114] In this implementation, the time and frequency domain response polynomials can be expressed as a matrix equation: λv=γ. The complex matrix equation is a linear algebraic equation system (which contains P+Q equations), which can be converted into a real matrix equation for solution. Moreover, if the number of equations is greater than the number of coefficients to be determined, the least squares method can be used for solution.

[0115] In one embodiment, Figure 2 As shown, the time domain response polynomial h(t p ) and the frequency domain response polynomial H(f q ) and optimizing the time domain scaling factor l t and the frequency domain scaling factor l f :

[0116] S600: Setting an initial expansion item number N0 and an initial scale factor, wherein the initial scale factor includes an initial time domain scale factor and an initial frequency domain scale factor;

[0117] S601: Calculate min{e(l based on the initial expansion item number N0 and the initial scale factor t)}, obtaining an optimal scale factor and a first polynomial coefficient, wherein the optimal scale factor includes an optimal time domain scale factor and an optimal frequency domain scale factor, wherein:

[0118]

[0119] S602: Calculate the expansion error and the matching error based on the polynomial coefficients of M to N expansion items, where 0<M<N, and the expansion error is: The matching error is the e(l t );

[0120] S603: Compare the expansion error with a preset expansion accuracy threshold, and compare the matching error with a preset matching accuracy threshold:

[0121] When the expansion error is greater than the expansion accuracy threshold, or the matching error is greater than the matching accuracy threshold, N0+I is used as a new initial expansion item number, the optimal scale factor is used as a new initial scale factor, and steps S601, S602, and S603 are repeated, where I is a positive integer (exemplarily, I=1);

[0122] When the expansion error is less than or equal to the expansion accuracy threshold, and the matching error is less than or equal to the matching accuracy threshold, the current number of expansion items and the scale factor are the optimized number of expansion items and the scale factor, and the optimization procedure ends.

[0123] This embodiment provides a method for solving the polynomial coefficient a n In the process of , a method is used to optimize the number of expansion items and scale factors. For different models, the optimal number of expansion items varies greatly. If the number of expansion items is too small, the interpolation and extrapolation results will be inaccurate; if the number of expansion items is too large, it will cause oscillation in the interpolation area or the extrapolation area, which will also cause distortion of the calculation results, and the scale factor will also change with the number of expansion items. For this situation where multiple parameters change at the same time, adopting a reasonable optimization algorithm to optimize the parameters can effectively improve the accuracy of the calculation results.

[0124] According to the convergence condition, the N-order expansion error can be obtained as:

[0125] The accumulated matching error of all early responses and low-frequency responses can be expressed as:

[0126]

[0127] Obviously, for a given number of expansion items, choosing different scale factors will result in different expansion coefficients and corresponding expansion errors and matching errors. Therefore, when the number of expansion items is given, the best choice is to minimize the matching error to obtain the optimal scale factor. This requires solving the nonlinear problem about the scale factor, i.e., min{e(l t )}. The number of expansion items N is determined by the matching error and the expansion error. To obtain the optimal expansion so that the calculation result meets the accuracy requirement, the above S600 to S603 can be used to optimize the number of expansion items and the scale factor.

[0128] In S601, the nonlinear problem about the scale factor is solved, that is, min{e(l t )} can be solved by using existing technology, for example, the gradient method (also called iterative method) can be used for solving.

[0129] In S602, M to N expansion items are the last expansion items, and the number of the last expansion items can be determined based on actual conditions or experience; for example, it can be any value in the range of 5 to 10.

[0130] It should be understood that the size of the serial number of each step in the above implementation does not mean the order of execution. The execution order of each step should be determined by its function and internal logic, and the serial number of each step should not constitute any limitation on the implementation process of the implementation method of the present invention.

[0131] In one embodiment, the present application provides a transient electromagnetic forward modeling device including induced polarization effect; Figure 3 As shown, the transient electromagnetic forward modeling device 30 including induced polarization effect may include a modeling module 31, an excitation module 32, a determination module 33, a first acquisition module 34 and a second acquisition module 35, which are respectively described as follows:

[0132] The modeling module 31 is used to construct a geological model, wherein the geological model includes resistivity parameters and polarizability parameters;

[0133] The excitation module 32 is used to apply an excitation field to the geological model;

[0134] The determination module 33 is used to determine the cutoff time T1 and the cutoff frequency F1 in the time period from T0 to T2 after the excitation field disappears; wherein T0<T1<T2, F min <F1<F max , F min is the minimum frequency in the period from T0 to T2, F max is the maximum frequency in the time period from T0 to T2;

[0135] The first acquisition module 34 is used to obtain the time domain electromagnetic data and F in the time period from T0 to T1 based on the first mathematical model and the Maxwell equations.min Frequency domain electromagnetic data within the F1 frequency range, the time domain electromagnetic data and the frequency domain electromagnetic data both contain induced polarization effect;

[0136] The second acquisition module 35 is used to obtain the response electromagnetic data containing induced polarization effect at any position in the geological model within the time period from T0 to T2 based on the frequency domain electromagnetic data, the time domain electromagnetic data and the second mathematical model.

[0137] In one embodiment, the second mathematical model uses orthogonal polynomials as basis functions.

[0138] In one embodiment, the orthogonal polynomial is a Laguerre orthogonal polynomial, and the second mathematical model includes a time domain response function and a frequency domain response function;

[0139] The time domain response function is: Among them, t represents time, a n are the polynomial coefficients, is the normalized Laguerre orthogonal polynomial, l t is the time domain scale factor, l t For limiting the definition domain of the time domain response function to the Within the scope of action,

[0140] The frequency domain response function is: Where f represents frequency, a n are the polynomial coefficients, Φ n By the Perform Fourier transform to obtain l f is the frequency domain scale factor, l f The frequency domain response function is limited to the domain of the frequency domain response function in the Φ n Within the scope of action, j is the imaginary unit.

[0141] In one embodiment, to obtain the polynomial coefficients a n , the second acquisition module 35 is used for:

[0142] Acquire P time domain sampling points from the time domain electromagnetic data;

[0143] Based on the time domain sampling points and the time domain response function, a discrete time domain response polynomial is obtained.

[0144]

[0145] Acquire Q frequency domain sampling points from the frequency domain response data;

[0146] Based on the frequency domain sampling points and the frequency domain response function, a discrete frequency domain response polynomial is obtained.

[0147]

[0148] Based on the time domain response polynomial h(t p ) and the frequency domain response polynomial H(f q ) Calculate the polynomial coefficient a n .

[0149] In one embodiment, in order to optimize the time domain response polynomial h(t p ) and the frequency domain response polynomial H(f q ) and optimizing the time domain scaling factor l t and the frequency domain scale factor lf, the second acquisition module 35 is used for:

[0150] Step 3510: Setting the initial expansion item number N0 and the initial scale factor, wherein the initial scale factor includes an initial time domain scale factor and an initial frequency domain scale factor;

[0151] Step 3520: Calculate min{e(l t )}, obtain the optimal scale factor and the first polynomial coefficient, wherein the optimal scale factor includes the optimal time domain scale factor and the optimal frequency domain scale factor, wherein

[0152]

[0153] Step 3530: Based on the polynomial coefficients of M to N expansion terms, calculate the expansion error and the matching error, where 0<M<N, and the expansion error is: The matching error is the e(l t );

[0154] Step 3540: Compare the expansion error with a preset expansion accuracy threshold, and compare the matching error with a preset matching accuracy threshold:

[0155] When the expansion error is greater than the expansion accuracy threshold, or the matching error is greater than the matching accuracy threshold, N0+I is used as the new initial number of expansion items, and the optimal scale factor is used as the new initial scale factor, and steps A to D are repeated, where I is a positive integer; when the expansion error is less than or equal to the expansion accuracy threshold, and the matching error is less than or equal to the matching accuracy threshold, the current number of expansion items and scale factor are the optimized number of expansion items and scale factor, and the optimization procedure ends.

[0156] In one embodiment, the first mathematical model is a Cole-Cole model: Where j is the imaginary unit, ω is the angular frequency; σ * (ω) is the resistivity; σ0 is the zero-frequency conductivity, which corresponds to the conductivity of the earth when there is no induced polarization effect; m is the charging rate, which characterizes the strength of the induced polarization effect of the earth, 0≤m≤1; τ is the time constant, which characterizes the slowness of the induced polarization effect of the earth; c is the frequency correlation coefficient, which characterizes the change of the phase of the induced polarization effect with frequency, 0≤c≤1.

[0157] In one embodiment, the first acquisition module 34 is further used for:

[0158] Gridding the geological model to obtain grid units; and

[0159] Based on the first mathematical model and the Maxwell equations, the time domain electromagnetic data and the frequency domain electromagnetic data of each grid unit are obtained.

[0160] In one embodiment, the first acquisition module 34 is further used for:

[0161] Determine the resistivity of the nodes of each grid cell;

[0162] Calculating the resistivity of any point in each grid unit according to the resistivity of the node and the interpolation function to obtain a grid unit with parameter assignment; and

[0163] Based on the first mathematical model, the Maxwell equations and the grid units with parameter assignments, the time domain electromagnetic data and the frequency domain electromagnetic data of each grid unit are obtained.

[0164] In one embodiment, the present application provides a forward modeling device; Figure 4 As shown, the forward modeling device 4 may include a processor 40 and a memory 41. The memory 41 is used to store a computer program 42, and the processor 40 is used to call and run the computer program 42 stored in the memory 41 to perform the steps in the above-mentioned embodiments of the transient electromagnetic forward modeling method containing induced polarization effect, for example Figure 1 Alternatively, the processor 40 is used to call and run the computer program 42 stored in the memory 41 to implement the functions of each module / unit in the above-mentioned device implementation, for example Figure 3 The functions of the modules 31 to 35 are shown.

[0165] Exemplarily, the computer program 42 may be divided into one or more modules / units, which are stored in the memory 41 and executed by the processor 40 to implement the technical solution provided by the present application. The one or more modules / units may be a series of computer program instruction segments capable of completing specific functions, which are used to describe the execution process of the computer program 42 in the controller 4. For example, the computer program 42 may be divided into Figure 3 Modules 31 to 35 are shown.

[0166] Those skilled in the art will understand that Figure 4 This is only an example of a forward modeling device, and does not constitute a limitation on the forward modeling device. The forward modeling device may include, but is not limited to, a processor 40 and a memory 41. The forward modeling device may include Figure 4 More or fewer components, or a combination of certain components, or different components are shown; for example, the forward modeling device may also include input and output devices, network access devices, buses, etc.

[0167] The processor 40 may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or any other conventional processor, etc.

[0168] The memory 41 may be an internal storage unit of the forward modeling device 4, such as a hard disk or memory of the forward modeling device 4. The memory 41 may also be an external storage device of the forward modeling device 4, such as a plug-in hard disk, a smart media card (SMC), a secure digital (SD) card, a flash card, etc. equipped on the forward modeling device 4. Further, the memory 41 may also include both an internal storage unit and an external storage device of the forward modeling device 4. The memory 41 is used to store the computer program 42 and other programs and data required by the forward modeling device 4. The memory 41 may also be used to temporarily store data that has been output or is to be output.

[0169] Those skilled in the art can clearly understand that the division of the above-mentioned functional units and modules is only for the convenience and simplicity of description, and does not constitute a limitation of the present application. In practical applications, the target function can be assigned to different functional units and modules as needed; that is, the internal structure of the device is divided into different functional units or modules to complete all or part of the target functions described above. The functional units and modules in the implementation method can be integrated into a processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit; the above-mentioned integrated unit can be implemented in the form of hardware or in the form of software functional units. In addition, the specific names of the functional units and modules are only for the convenience of distinguishing each other, and are not used to limit the scope of protection of the present application. The specific working process of the units and modules in the above-mentioned devices and equipment can refer to the corresponding description in the aforementioned method implementation method, which will not be repeated here.

[0170] For each of the above embodiments, the description has different emphases, either focusing on method steps or focusing on devices and equipment; for parts that are not described or recorded in detail in a certain embodiment, reference can be made to the relevant descriptions of other embodiments.

[0171] Those skilled in the art will appreciate that the functional units and algorithm steps involved in each embodiment of the present application can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether they are executed in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, and these implementations should not be considered to exceed the scope of the present invention.

[0172] It is understood that the devices, equipment and methods provided in the embodiments of the present application are merely illustrative, and they can also be implemented in other ways. For example, the division of the modules or units is only a logical function division, and there may be other division methods in actual implementation; for example, multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the coupling, communication and connection between the modules / units can be implemented in electrical, mechanical or other forms through some interfaces.

[0173] The units provided as separate components may or may not be physically separated, and the components provided as modules / units may or may not be physical units; that is, they may be located in one place or distributed on multiple network units. Those skilled in the art may select some or all of the modules / units according to actual needs to achieve the purpose of the embodiments of the present application.

[0174] In one embodiment, the present application further provides a computer-readable storage medium, in which a computer program is stored. When the computer program is executed by a processor, the technical solution in the embodiment of the transient electromagnetic forward modeling method containing induced polarization effect provided in the present application is implemented.

[0175] The computer-readable storage medium can be sold or used as an independent product. Based on this understanding, the various method steps provided in this application can be completed by controlling the relevant hardware through a computer program; wherein the computer program includes a computer program code, and the computer program code can be in source code form, object code form, executable file or some intermediate form, etc. The computer-readable medium can be: any entity or device that can carry the computer program code, a recording medium, a USB flash drive, a mobile hard disk, a magnetic disk, an optical disk, a computer memory, a read-only memory (ROM), a random access memory (RAM), an electrical carrier signal, a telecommunication signal, and a software distribution medium, etc.

[0176] The above-described embodiments are only used to illustrate the technical solutions of the present application, rather than to limit them. Although the present application has been described in detail with reference to the aforementioned embodiments, a person skilled in the art should understand that it is still possible to make simple modifications to the technical solutions described in the above embodiments, or to make equivalent replacements for some of the technical features therein. These simple modifications or replacements cannot cause the modified technical solutions to essentially deviate from the inventive concept / spirit and scope of the present application, and they should all be deemed to be included in the protection scope of the present application.

Claims

1. A transient electromagnetic forward modeling method including induced polarization effect, It is characterized in that include: Constructing a geological model, wherein the geological model includes resistivity parameters and polarizability parameters; applying an excitation field to the geological model; After the excitation field disappears, T 0 To T 2 Within the time period, determine the cutoff time T 1 and the cutoff frequency F 1 ; where T 0 <T 1 <T 2 , F min <F 1 <F max , F min T 0 To T 2 Minimum frequency in the time period, F max T 0 To T 2 Maximum frequency in the time period; Based on the first mathematical model and Maxwell equations, we obtain T 0 To T 1 Time domain electromagnetic data and F min To F 1 Frequency domain electromagnetic data within a frequency range, both the time domain electromagnetic data and the frequency domain electromagnetic data contain induced polarization effects; and Based on the frequency domain electromagnetic data, the time domain electromagnetic data and the second mathematical model, T is obtained. 0 To T 2 Response electromagnetic data including induced polarization effect at any position in the geological model within a time period.

2. The transient electromagnetic forward modeling method according to claim 1, It is characterized in that The second mathematical model uses orthogonal polynomials as basis functions.

3. The transient electromagnetic forward modeling method according to claim 2, It is characterized in that The orthogonal polynomial is a Laguerre orthogonal polynomial, and the second mathematical model includes a time domain response function and a frequency domain response function; The time domain response function is: Among them, t represents time, a n are the polynomial coefficients, is the normalized Laguerre orthogonal polynomial, l t is the time domain scale factor, l t For limiting the definition domain of the time domain response function to the Within the scope of action, The frequency domain response function is: Where f represents frequency, a n are the polynomial coefficients, Φ n By the Perform Fourier transform to obtain l f is the frequency domain scale factor, l f The frequency domain response function is limited to the domain of the frequency domain response function in the Φ n Within the scope of action, j is the imaginary unit.

4. The transient electromagnetic forward modeling method according to claim 3, It is characterized in that The polynomial coefficient a n Obtained by the following steps: Acquire P time domain sampling points from the time domain electromagnetic data; Based on the time domain sampling points and the time domain response function, a discrete time domain response polynomial is obtained. Acquire Q frequency domain sampling points from the frequency domain response data; Based on the frequency domain sampling points and the frequency domain response function, a discrete frequency domain response polynomial is obtained. Based on the time domain response polynomial h(t p ) and the frequency domain response polynomial H(f q ) Calculate the polynomial coefficient a n .

5. The transient electromagnetic forward modeling method according to claim 4, It is characterized in that The time domain response polynomial h(t p ) and the frequency domain response polynomial H(f q ) and optimizing the time domain scaling factor l t and the frequency domain scaling factor l f : Step A: Set the initial number of expansion items N 0 and an initial scale factor, wherein the initial scale factor comprises an initial time domain scale factor and an initial frequency domain scale factor; Step B: Based on the initial number of expansion items N 0 and the initial scale factor to calculate min{e(l t )}, obtain the optimal scale factor and the first polynomial coefficient, wherein the optimal scale factor includes the optimal time domain scale factor and the optimal frequency domain scale factor, wherein Step C: Based on the polynomial coefficients of M to N expansion terms, the expansion error and the matching error are calculated, where 0<M<N, and the expansion error is: The matching error is the e(l t ); Step D: Compare the expansion error with a preset expansion accuracy threshold, and compare the matching error with a preset matching accuracy threshold: When the expansion error is greater than the expansion accuracy threshold, or the matching error is greater than the matching accuracy threshold, N 0 +I is used as the new initial number of expanded items, the optimal scale factor is used as the new initial scale factor, and steps A to D are repeated, where I is a positive integer; when the expansion error is less than or equal to the expansion accuracy threshold, and the matching error is less than or equal to the matching accuracy threshold, the current number of expanded items and scale factor are the optimized number of expanded items and scale factor, and the optimization procedure ends.

6. The transient electromagnetic forward modeling method according to claim 1, It is characterized in that The first mathematical model is the Cole-Cole model: Where j is the imaginary unit, ω is the angular frequency; σ * (ω) is the resistivity; σ 0 is the zero-frequency conductivity, corresponding to the conductivity of the earth when there is no induced polarization effect; m is the charging rate, which characterizes the strength of the induced polarization effect of the earth, 0≤m≤1; τ is the time constant, which characterizes the slowness of the induced polarization effect of the earth; c is the frequency correlation coefficient, which characterizes the change of the induced polarization effect phase with frequency, 0≤c≤1.

7. The transient electromagnetic forward modeling method according to claim 1, It is characterized in that Based on the first mathematical model and Maxwell equations, T is obtained 0 To T 1 Time domain electromagnetic data and F min To F 1 Frequency domain electromagnetic data in the frequency range, including: Gridding the geological model to obtain grid units; and Based on the first mathematical model and the Maxwell equations, the time domain electromagnetic data and the frequency domain electromagnetic data of each grid unit are obtained.

8. The transient electromagnetic forward modeling method according to claim 7, It is characterized in that The step of obtaining the time domain electromagnetic data and the frequency domain electromagnetic data of each grid unit based on the first mathematical model and the Maxwell equations includes: Determine the resistivity of the nodes of each grid cell; Calculating the resistivity of any point in each grid unit according to the resistivity of the node and the interpolation function to obtain a grid unit with parameter assignment; and Based on the first mathematical model, the Maxwell equations and the grid units with parameter assignments, the time domain electromagnetic data and the frequency domain electromagnetic data of each grid unit are obtained.

9. A transient electromagnetic forward modeling device with induced polarization effect, It is characterized in that include: A modeling module, which is used to construct a geological model, wherein the geological model includes resistivity parameters and polarizability parameters; an excitation module, which is used to apply an excitation field to the geological model; A determination module is used to determine T after the excitation field disappears. 0 To T 2 Within the time period, determine the cutoff time T 1 and the cutoff frequency F 1 ; where T 0 <T 1 <T 2 , F min <F 1 <F max , F min T 0 To T 2 Minimum frequency in the time period, F max T 0 To T 2 Maximum frequency in the time period; The first acquisition module is used to obtain T based on the first mathematical model and the Maxwell equations. 0 To T 1 Time domain electromagnetic data and F min To F 1 Frequency domain electromagnetic data within a frequency range, the time domain electromagnetic data and the frequency domain electromagnetic data both contain induced polarization effect; as well as The second acquisition module is used to obtain T based on the frequency domain electromagnetic data, the time domain electromagnetic data and the second mathematical model. 0 To T 2 Response electromagnetic data including induced polarization effect at any position in the geological model within a time period.

10. A computer-readable storage medium storing a computer program. It is characterized in that When the computer program is executed by a processor, the transient electromagnetic forward modeling method including induced polarization effect according to any one of claims 1 to 8 is implemented.

Images