Transient electromagnetic two-parameter joint inversion method and device and storage medium

Abstract

The invention provides a transient electromagnetic two-parameter joint inversion method, a transient electromagnetic two-parameter joint inversion device and storage equipment, and relates to the technical field of electric digital data processing. The method comprises the following steps: acquiring detected actual electromagnetic data; establishing a geological model; performing grid division on the geological model; performing parameter assignment on each grid unit to form an initial parameter model; based on the geologic model, constructing a theoretical mathematical model by using a forward modeling method; constructing an objective function, wherein the objective function comprises an FCM clustering item; optimizing the target function to obtain a target geological model; according to the transient electromagnetic two-parameter joint inversion method provided by the invention, the problem that the two-parameter joint inversion is difficult to realize due to weak correlation between the resistivity and the polarizability in the prior art 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 dual-parameter joint inversion method, device and storage medium. 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] Transient electromagnetic data containing induced polarization (IP) effect (also known as "induced polarization effect") contain resistivity parameters and polarizability parameters. The correlation between resistivity and polarizability is weak, making it difficult to achieve joint inversion of these two parameters, and the inversion results cannot accurately reflect the geological conditions of the target area. Summary of the invention

[0005] The present application provides a transient electromagnetic dual-parameter joint inversion method, device and storage medium to solve the problem in the prior art that dual-parameter joint inversion is difficult to achieve due to weak correlation between resistivity and polarizability.

[0006] In order to achieve the above object, the first aspect of the present application provides a transient electromagnetic dual-parameter joint inversion method, which comprises the following steps:

[0007] Acquire the actual electromagnetic data detected;

[0008] Establish a geological model, the geological model includes a resistivity parameter m ρ and the polarizability parameter m η ;

[0009] Dividing the geological model into grids to obtain M grid units;

[0010] Assign parameters to each grid unit to form an initial parameter model mapr ,in,

[0011] is the initial resistivity parameter, is the initial polarizability parameter;

[0012] Based on the geological model, a theoretical mathematical model A(m) is constructed using a forward modeling method. The theoretical mathematical model includes the resistivity parameter m ρ and the polarizability parameter m η ;

[0013] Construct the objective function P(m),

[0014] in:

[0015] m is the parameter to be solved, m=(m ρ , m η ) T ;

[0016] W d is a covariance matrix constructed based on the actual electromagnetic data;

[0017] d obs is an actual data matrix constructed based on the actual electromagnetic data;

[0018] λ is the regularization factor, λ=(λ ρ ,λ η ) T , where λ ρ is the resistivity regularization factor, λ η is the polarizability regularization factor;

[0019] W m is the least squares model weighting matrix;

[0020] is the fuzzy C-means clustering term, β is the clustering term weight factor;

[0021] The objective function is optimized to obtain a target geological model.

[0022] In one possible implementation, in:

[0023] M is the number of grid cells in the geological model;

[0024] C is the number of cluster centers;

[0025] m j is the parameter to be solved for the jth grid unit, m j =(m ρj ,m ηj )T , where m ρj is the resistivity of the jth grid cell, m ηj is the polarizability of the jth grid cell;

[0026] v k is the kth cluster center, v k =(v ρk ,v ηk ) T , where v ρk is the cluster center of the kth resistivity, v ηk is the cluster center of the kth polarizability;

[0027] u jk is the degree of membership of the unsolved parameter of the jth grid unit to the kth cluster center;

[0028] q is the fuzzification parameter;

[0029] t k is the kth reference cluster center, t k =(t ρk ,t ηk ) T , where t ρk is the reference cluster center of the kth resistivity, t ηk is the reference cluster center of the kth polarizability;

[0030] ω k is the kth reference cluster center t k The weight factor of the kth reference cluster center t k confidence level.

[0031] In a possible implementation, before optimizing the objective function, the following steps are included:

[0032] Determine the resistivity parameter m of the node of each grid unit of the geological model ρ ;and

[0033] According to the resistivity parameter m of the node ρ and interpolation function to calculate the resistivity parameter m at any point in each grid cell ρ , get the grid cells assigned by the parameters.

[0034] In a possible implementation, the forward modeling method includes the following steps:

[0035] applying a forward excitation field to the geological model;

[0036] After the forward excitation field disappears, T 0 To T 2 Within the time period, determine the cutoff time T1 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;

[0037] Based on the first forward mathematical model and Maxwell equations, T 0 To T 1 Theoretical time domain electromagnetic data and F in the time period min To F 1 Theoretical frequency domain electromagnetic data within the frequency range, wherein both the theoretical time domain electromagnetic data and the theoretical frequency domain electromagnetic data contain induced polarization effect;

[0038] Based on the theoretical frequency domain electromagnetic data, the theoretical time domain electromagnetic data and the second forward mathematical model, T 0 To T 2 Theoretical response electromagnetic data including induced polarization effect at any position in the geological model within a time period; and

[0039] The theoretical mathematical model A(m) is constructed according to the theoretical response electromagnetic data.

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

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

[0042] 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,

[0043] The frequency domain response function is: Where f represents frequency, an 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.

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

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

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

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

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

[0049] 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 .

[0050] 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 :

[0051] Step A: Set the initial number of expansion items N 0 , an initial scale factor, wherein the initial scale factor includes an initial time domain scale factor and an initial frequency domain scale factor;

[0052] 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

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

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

[0055] 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.

[0056] The second aspect of the present application provides a transient electromagnetic dual-parameter joint inversion device, the device comprising:

[0057] An acquisition module, used for acquiring the actual electromagnetic data detected;

[0058] The first modeling module is used to establish a geological model, wherein the geological model includes a resistivity parameter m ρ and the polarizability parameter m η ;

[0059] A grid division module, used for performing grid division on the geological model to obtain M grid units;

[0060] The assignment module is used to assign parameters to each grid unit to form an initial parameter model m apr ,in, is the initial resistivity parameter, is the initial polarizability parameter;

[0061] The second modeling module is used to construct a theoretical mathematical model A(m) based on the geological model using a forward modeling method. The theoretical mathematical model includes the resistivity parameter m ρ and the polarizability parameter m η ;

[0062] The third modeling module is used to construct the objective function P(m).

[0063]

[0064] m is the parameter to be solved, m=(m ρ , mη ) T ;

[0065] W d is a covariance matrix constructed based on the actual electromagnetic data;

[0066] d obs is an actual data matrix constructed based on the actual electromagnetic data;

[0067] λ is the regularization factor, λ=(λ ρ ,λ η ) T , where λ ρ is the resistivity regularization factor, λ η is the polarizability regularization factor;

[0068] W m is the least squares model weighting matrix;

[0069] is the fuzzy C-means clustering term, β is the clustering term weight factor;

[0070] The optimization module is used to optimize the objective function to obtain a target geological model.

[0071] The third aspect of the present application provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements any transient electromagnetic dual-parameter joint inversion method provided in the first aspect of the present application.

[0072] Compared with the prior art, the transient electromagnetic dual-parameter joint inversion method provided by the present application has the following beneficial effects: based on the fact that resistivity and polarizability have good clustering characteristics, the FCM clustering term is introduced. This realizes the transient electromagnetic multi-parameter joint inversion, solving the problem in the prior art that the two-parameter joint inversion is difficult to achieve due to the weak correlation between resistivity and polarizability. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] 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.

[0074] Figure 1 A flowchart for implementing the transient electromagnetic dual-parameter joint inversion method provided in the embodiment of the present application;

[0075] Figure 2A flowchart for implementing the forward modeling method in the transient electromagnetic dual-parameter joint inversion method provided in the embodiment of the present application;

[0076] Figure 3 A flowchart for implementing a method for optimizing polynomial coefficients and scale factors in a forward modeling method in a transient electromagnetic dual-parameter joint inversion method provided in an embodiment of the present application;

[0077] Figure 4 A schematic diagram of the structure of a transient electromagnetic dual-parameter joint inversion device provided in an embodiment of the present application;

[0078] Figure 5 A schematic diagram of the structural composition of the inversion device provided in an embodiment of the present application. DETAILED DESCRIPTION

[0079] 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.

[0080] 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.

[0081] At present, transient electromagnetic inversion is mainly based on single-parameter inversion of resistivity. Multi-parameter inversion is different from single-parameter inversion. It not only needs to consider the data constraints of inversion, but also needs to consider the relationship between physical parameters, and thus construct parameter constraints.

[0082] In multi-parameter joint inversion, the main methods for constructing parameter constraints are physical property constraints and structural constraints. When there is a clear rock physics coupling relationship between different physical parameters (such as resistivity and polarizability), the coupling relationship can be used to construct parameter constraints. However, this coupling relationship is usually obtained by an empirical formula formed by multiple experiments, and is only applicable to typical local areas. Only when there is an accurate expression relationship between different physical parameters, the coupling relationship between them has universal significance. The above parameter constraints can also be constructed through structural constraints. Structural constraints mean that when the type of rock or ore changes, the physical parameters will undergo similar changes. Based on structural constraints, the cross-gradient method is produced in the prior art, that is, the structural coupling relationship between different physical parameters is constructed through a cross-gradient operator, thereby realizing the joint inversion of the spatial structural similarity of different physical parameters.

[0083] In order to make the purpose, technical solutions and advantages of the present application clearer, specific embodiments will be described below in conjunction with the accompanying drawings.

[0084] In the process of realizing the present invention, the inventors found that although the spatial correlation between resistivity and polarizability is weak, typical ore bodies all have classification characteristics such as low resistance and high polarizability or high resistance and high polarizability, and the intersection of parameters has good clustering characteristics. Based on this, the inventors proposed a multi-parameter joint inversion method based on fuzzy clustering constraints.

[0085] See also Figure 1 In one embodiment, the dual-parameter joint inversion method provided in the present application comprises the following steps:

[0086] S100: Acquire the actual electromagnetic data detected;

[0087] S200: Establishing a geological model, wherein the geological model includes a resistivity parameter m ρ and the polarizability parameter m η ;

[0088] S300: Gridding the geological model to obtain M grid units;

[0089] S400: Assign parameters to each grid unit to form an initial parameter model m apr , in is the initial resistivity parameter, is the initial polarizability parameter;

[0090] S500: Based on the geological model, a theoretical mathematical model A(m) is constructed using a forward modeling method. The theoretical mathematical model includes the resistivity parameter m ρ and the polarizability parameter m η ;

[0091] S600: constructing an objective function P(m),

[0092]

[0093] m is the parameter to be solved, m=(m ρ , m η ) T ;

[0094] W d is a covariance matrix constructed based on the actual electromagnetic data;

[0095] d obs is an actual data matrix constructed based on the actual electromagnetic data;

[0096] λ is the regularization factor, λ=(λ ρ ,λ η ), where λ ρ is the resistivity regularization factor, λ η is the polarizability regularization factor;

[0097] W m is the least squares model weighting matrix;

[0098] is the fuzzy C-means clustering term, β is the clustering term weight factor;

[0099] S700: Optimizing the objective function to obtain a target geological model.

[0100] When performing transient electromagnetic inversion, a field experiment is conducted; that is, a transmitting device is used to transmit a primary electromagnetic field to a target area, and then a receiving device is used to receive a secondary electromagnetic field (ie, the actual electromagnetic data in S100).

[0101] Before inversion, it is necessary to establish a geological model corresponding to the above-mentioned target area (i.e., the geological model in S200). The geological model is generally established using computer software (such as COMSOL, GID) and is gridded, and each grid unit has its own resistivity and polarizability. The geological model can represent the underground resistivity and polarizability distribution of the target area.

[0102] The resistivity and polarizability distribution represented by the initial geological model is generally inconsistent with the actual situation. The inversion process is the process of correcting the initial geological model using the actual electromagnetic data measured by field experiments. The inverted geological model (i.e., the target geological model obtained in S700) can represent the actual resistivity and polarizability distribution of the underground of the target area, and the distribution can be used to predict whether there are mineral deposits in the underground of the target area and other issues of interest.

[0103] It can be seen from the objective function P(m) constructed in S600 that the inversion will use the theoretical mathematical model A(m) established by the forward modeling. Through the theoretical mathematical model A(m), the electromagnetic response data of the geological model in theory at any target time can be calculated. In this embodiment, the existing forward modeling method can be used.

[0104] The objective function P(m) is a function formed by introducing the Fuzzy C-means Algorithm (FCM) to perform joint constraints on the basis of the regularization constraints of the classic minimum structure model. In the objective function P(m), It is used to characterize the difference between actual electromagnetic data and theoretical electromagnetic data; where W d is a covariance matrix constructed based on actual electromagnetic data, which is used to prevent the error in the actual electromagnetic data from having too much impact on the optimization of the objective function P(m). Specifically, the actual electromagnetic data is generally a waveform that changes with time; accordingly, d obs Generally, the matrix is ​​obtained by discretizing the waveform based on time. For example, d obs It can be dB / dt corresponding to a series of moments, where B represents magnetic induction intensity, t represents time, and d is the differential sign. Accordingly, A(m) can represent the theoretical electromagnetic data at the corresponding moments.

[0105] The Gauss-Newton iteration method can be used to optimize the objective function P(m); that is, apr is the initial value, based on a certain step size or model change (can be positive or negative, Δm ρ and Δm η ), so that the parameter m to be solved of each grid unit changes until the preset iteration stop condition is reached (for example, the preset number of iteration steps or the preset P(m) accuracy is reached), and the target geological model is obtained. At this time, it can be considered that the resistivity and polarizability of each grid unit are consistent with the actual situation, so as to obtain the actual resistivity and polarizability distribution of the underground in the target area.

[0106] In the objective function P(m), is a regular constraint term, which is used to speed up the convergence of the function. apr It is constantly changing, starting from the initial value and being continuously given a better solution.

[0107] In the objective function P(m), It is an FCM clustering item; it can establish the connection between resistivity and polarizability based on the classification characteristics such as low resistance and high polarization or high resistance and high polarization discovered by the inventor.

[0108] Compared with the prior art, the inversion method provided in this embodiment introduces the FCM clustering term based on the fact that resistivity and polarizability have good clustering characteristics. This realizes the transient electromagnetic multi-parameter joint inversion, solving the problem in the prior art that the two-parameter joint inversion is difficult to achieve due to the weak correlation between resistivity and polarizability.

[0109] For example, in S400, the initial parameter model m apr It can be obtained by using existing methods, for example, it can be obtained by using a one-dimensional inversion method based on a priori constraint objective function. For example, the initial parameter model m can be determined first apr , and then determine the optimal grid division based on it.

[0110] In one embodiment, the FCM clustering term middle,

[0111] in:

[0112] M is the number of grid cells in the geological model;

[0113] C is the number of cluster centers;

[0114] m j is the parameter to be solved for the jth grid unit, m j =(m ρj ,m ηj ) T , where m ρj is the resistivity of the jth grid cell, m ηj is the polarizability of the jth grid cell;

[0115] v k is the kth cluster center, v k =(v ρk ,v ηk ) T , where v ρk is the cluster center of the kth resistivity, v ηk is the cluster center of the kth polarizability;

[0116] u jk is the degree of membership of the unsolved parameter of the jth grid unit to the kth cluster center;

[0117] q is a fuzzification parameter (exemplarily, q=2);

[0118] t k is the kth reference cluster center, t k =(t ρk ,t ηk ) T , where t ρk is the reference cluster center of the kth resistivity, t ηk is the reference cluster center of the kth polarizability;

[0119] ω k is the kth reference cluster center t k The weight factor of the kth reference cluster center t k confidence level.

[0120] This embodiment provides a The specific form of u jk Non-negative and satisfies:

[0121] The process of implementing the FCM algorithm is to minimize the clustering constraint function in an iterative convergence process. When the minimum value of For the membership value u jk and cluster center v k Take the derivative and take zero, and use the Lagrange multiplier method to get u jk 、v k Finally, the Gauss-Newton algorithm can be used to optimize and solve the objective function, and the belonging of the object to be classified can be determined according to the value of the elements in the solved membership matrix.

[0122] Specifically, u is obtained by the Lagrange multiplier method. jk 、v k The expressions are:

[0123]

[0124] FCM clustering terms The first term in can be expressed as:

[0125]

[0126] In this implementation, the specific process of optimizing the parameter m to be solved of the objective function using the Gauss-Newton method is as follows:

[0127] Taylor expansion is performed on the theoretical mathematical model A(m) constructed by the forward modeling method, and the result of the n+1th iteration of the objective function is:

[0128] in:

[0129] m (n) is the nth resistivity and polarizability model vector;

[0130] J (n) is the Jacobian matrix, in is the sensitivity matrix of the nth resistivity, is the sensitivity matrix of the nth polarizability;

[0131] Δm (n) is the model change of the n+1th resistivity and polarizability,

[0132]

[0133]

[0134] The objective function is iterated for the n+1th time. ( n+ 1) (m) vs. Δm (n) Taking the derivative and taking 0, we get the following Gauss-Newton equation:

[0135] in:

[0136] Where I is the identity matrix, W x , W z are the roughness matrices in the horizontal and vertical directions, respectively.

[0137] Cluster centers of the n+1th iteration and They are:

[0138]

[0139] The objective function is iterated for the n+1th time. ( n+ 1) (m) vs. Δm (n) Taking the derivative and taking 0, the Gauss-Newton equation above can be expressed as:

[0140] Ex = f, where:

[0141] x=Δm (n) ;

[0142]

[0143] The stable biconjugate gradient method can be used to optimize the Gauss-Newton equations. Exemplarily, the iterative algorithm is as follows:

[0144] (1) Let x 0 =Δm 0 =0, let the initial residual be x 0 =Δm 0 =0;

[0145] (2) Initialize the number of iterations i to 0 and set the initial search direction to p 0 =r 0 ;

[0146] (3) Calculate step size α i=(r i ,r 0 ) / (Er i ,r 0 ), calculation parameters:

[0147]

[0148] (4) Update the model change x i+1 =x i +α i p i +σ i s i , update the residual r i+1 =s i -σ i Es i ;

[0149] (5) Calculation Update search direction p i+1 =r i+1 +β i (p i -σ i Ep i ), and let k = k + 1.

[0150] Repeat the above steps (1) to (5) until the preset iteration stop condition is reached, the iteration is stopped, and the optimal model change parameter is obtained. Then, the final parameter m to be solved is obtained according to the search step size, and each element in the final parameter m to be solved is assigned to each grid unit of the geological model, and the target geological model in S700 can be obtained. The final parameter m to be solved is: m (n+1) =m (n) +μΔm (n) , where: μ is the variable Δm along the model (n) The search step length, using linear search, is expressed as:

[0151] Where ξ is a very small positive real number.

[0152] In one embodiment, parameter assignment may be performed on each grid cell of the geological model by the following steps:

[0153] S801: Determine the resistivity parameter m of the node of each grid unit of the geological model ρ ;and

[0154] S802: According to the resistivity parameter m of the node ρ and interpolation function to calculate the resistivity parameter m at any point in each grid cell ρ , get the grid cells assigned by the parameters.

[0155] In the prior art, when building the relationship between parameters based on structural constraints and realizing multi-parameter joint inversion, physical property parameters are often assigned to each grid unit. During the iterative optimization process, when the parameter model m undergoes a significant change, it means that the structural characteristics of the geological model medium have undergone a significant change, so the geological model needs to be re-divided and optimized, and the grid optimization requires a lot of time. This will reduce the forward and inversion speed.

[0156] In this embodiment, the resistivity at the same position in the grid unit is represented by the node resistivity through a scalar interpolation function. Compared with the traditional vector finite element, the resistivity is assigned to the node instead of the unit, and the resistivity in the unit is represented by interpolation. In this way, the geological model does not need to consider the structural characteristics of the medium when meshing, thereby avoiding re-meshing, greatly reducing the workload of finite element pre-processing. The method provided in this embodiment has the advantages of the meshless method, which can reduce the difficulty of complex medium segmentation and is suitable for the forward modeling of resistivity gradient geological models such as deep mineral resources. At the same time, it can be unified with the inversion grid to avoid repeated updating and optimization of the grid during inversion.

[0157] In one embodiment, see Figure 2 , the forward modeling method used in S500 includes the following steps:

[0158] S501: applying a forward excitation field to the geological model;

[0159] S502: T after the forward 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;

[0160] S503: Based on the first forward mathematical model and Maxwell equations, obtain T 0 To T 1 Theoretical time domain electromagnetic data and F in the time period min To F 1 Theoretical frequency domain electromagnetic data within the frequency range, the theoretical time domain electromagnetic data and the theoretical frequency domain electromagnetic data both contain induced polarization effect;

[0161] S504: Based on the above theoretical frequency domain electromagnetic data, theoretical time domain electromagnetic data and the second forward mathematical model, obtain T 0 To T 2 Theoretical response electromagnetic data including induced polarization effect at any location in the above geological model within the time period; and

[0162] S505: Construct the theoretical mathematical model A(m) according to the theoretical response electromagnetic data.

[0163] The forward modeling excitation field in S501 is generally an excitation field simulated by a computer, and the forward modeling excitation field is consistent with the excitation field that generates the actual electromagnetic data in S100. In practice, the excitation field is generally generated by a transmitting device; the transmitting device has various forms, and the common ones are electric dipoles, magnetic dipoles, grounded power lines, and ungrounded return lines.

[0164] In S502, the time period studied by the forward modeling is determined, namely, T 0 To T 2 , which can be called the sampling period. Forward modeling is to derive the theoretical response electromagnetic data of any point in the geological model during the sampling period, such as the magnitude and direction of the theoretical response magnetic field. For example, the timing can be started from the time when the forward modeling excitation field disappears (is turned off), that is, T 0 =0, T 2 It can be tens to thousands of milliseconds (ms). For example, the turn-off time can be 1μs, starting from the time of turn-off, T 0 =0, T 2 =200ms.

[0165] In S502, a cut-off time, T 1 . Determine T 1 The sampling period under study can be divided into a high frequency period and a low frequency period. The frequency of the response magnetic field decays with time. Before the cutoff time, i.e., T 0 To T 1 The time period can be called the high-frequency period and the early period; after the cutoff time, that is, T 1 To T 2 time periods, which can be called low-frequency periods and late periods. In order to determine T 1 , in a relatively simple scenario without induced polarization effect, we can study 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, that is, T 1 The research on transient electromagnetic without induced polarization effect is very mature. 1 It can be determined based on these previous experiences. For example, for different geological models, when the sampling period is 10ms to 100ms, T 1The value can be in the range of 0.01ms to 1ms.

[0166] In S502, the cutoff frequency, F 1 Before the cutoff frequency, that is, F 0 To F 1 The interval can be called the high-frequency period and the early period; after the truncation frequency, that is, F 1 To F 2 The interval can be called the low-frequency period and the late period. 1 , in a relatively simple scenario without induced polarization effect, we can 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-low frequency dividing point, that is, F 1 The research on transient electromagnetic without induced polarization effect is very mature. 1 It 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 [10 Hz, 100 kHz], the cutoff frequency can be taken in the interval [100 Hz, 1000 Hz].

[0167] In S503, the theoretical time domain electromagnetic data of the high frequency period (early period) and the theoretical frequency domain electromagnetic data of the low frequency period (late period) are calculated respectively using the first forward mathematical model and the Maxwell equations. This avoids the calculation of frequency domain data in the high frequency region. In the high frequency region, 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 region, 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.

[0168] It can be understood that the first forward 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.

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

[0170] Specifically, based on the second forward mathematical model and the theoretical frequency domain electromagnetic data and the theoretical time domain electromagnetic data obtained in S503, the theoretical response electromagnetic data in the entire time and frequency domains are deduced, which may include:

[0171] Based on the second forward mathematical model and theoretical frequency domain electromagnetic data, the entire frequency band F min To F max The full-band theoretical frequency domain response electromagnetic data within; and

[0172] Based on the second forward mathematical model and theoretical time-domain electromagnetic data, the entire period T 0 To T 2 The full-time theoretical time domain response electromagnetic data within.

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

[0174] 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 S503, theoretical time domain electromagnetic data are obtained in the high-frequency region, and theoretical frequency domain electromagnetic data are obtained 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, so that the three-dimensional rapid forward modeling of transient electromagnetics containing induced polarization effect can be realized. In S504, a joint derivation of the time and frequency domains is performed to obtain the theoretical response electromagnetic data of the entire time and frequency domains, thereby realizing the rapid forward modeling of transient electromagnetics containing induced polarization effect.

[0175] For a target area, if the theoretical response electromagnetic data obtained by forward modeling is consistent with the actual response electromagnetic data obtained in the inversion, it indicates that the geological model initially established and iterated 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 the inversion and even the TEM. The forward modeling method provided in this embodiment has fast calculation speed and high efficiency, thereby effectively ensuring that the transient electromagnetic dual-parameter joint inversion is successfully implemented and smoothly implemented.

[0176] 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.

[0177] In one implementation, the first forward mathematical model in S503 is a 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 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.

[0178] The left side of the Cole-Cole model equal sign is resistivity, and the right side is the parameter that characterizes the induced polarization effect. 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 established in S200 and the forward excitation field given in S501, the Cole-Cole model combined with the Maxwell equations can calculate the theoretical electromagnetic data containing the induced polarization effect at any position in the geological model.

[0179] 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 the time period when the early complex resistivity changes little and will not have a significant impact on the transient electromagnetic response, so as to reasonably determine the cutoff time T 1 .

[0180] In one embodiment, in S503, based on the first forward mathematical model and the Maxwell equations, T is obtained. 0 To T 1 Theoretical time domain electromagnetic data and F in the time period min To F 1 Theoretical frequency domain electromagnetic data within the frequency range may include the following steps:

[0181] S5031: Gridding the geological model in S200 to obtain grid units; and

[0182] S5032: Based on the first forward mathematical model and the Maxwell equations, obtain theoretical time domain electromagnetic data and theoretical frequency domain electromagnetic data of each grid cell.

[0183] 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 theoretical electromagnetic data at any position of the geological model in the time domain. As for the theoretical frequency domain electromagnetic data, the time domain can be converted to the frequency domain by Fourier transform.

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

[0185] S5033: Determine the theoretical resistivity of the nodes of each grid unit;

[0186] S5034: Calculate the theoretical resistivity of any point in each grid unit according to the theoretical resistivity of the node and the interpolation function to obtain a grid unit with parameter assignment; and

[0187] S5035: Based on the first forward mathematical model, the Maxwell equations and the grid units with parameter assignments, obtain the theoretical time domain electromagnetic data and the theoretical frequency domain electromagnetic data of each grid unit.

[0188] Compared with the traditional vector finite element method, in this embodiment, the theoretical 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 in 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 in 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.

[0189] 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.

[0190] Exemplarily, in this embodiment, when the theoretical time-domain electromagnetic data is calculated using the first forward 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.

[0191] In this example, in terms of loading the forward modeling 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.

[0192] Exemplarily, in this embodiment, when the theoretical frequency domain electromagnetic data is calculated using the first forward 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.

[0193] 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 forward 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.

[0194] In this example, the abnormal field method can be used for loading the forward 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.

[0195] In one implementation, the second forward mathematical model in S504 uses orthogonal polynomials as basis functions.

[0196] 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.

[0197] 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 theoretical time domain electromagnetic data and the late theoretical frequency domain electromagnetic data obtained in S503, thereby obtaining a complete theoretical response in the entire time domain and frequency domain.

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

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

[0200] 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,

[0201] 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.

[0202] In this embodiment, Laguerre orthogonal polynomials are used as the basis functions of the time domain response h(t). Compared with Hermite orthogonal polynomials, Laguerre orthogonal polynomials are more suitable as the expansion basis functions of the electromagnetic response. Because the domain of Hermite orthogonal polynomials is [-∞, +∞], it is necessary to perform t 0 The delay factor can match the time domain response. If the delay factor is not selected properly, it will have a great impact on the result. The domain of the Laguerre orthogonal polynomial is [0, +∞], and it can match the time domain response without any delay; and the lack of delay variables can more effectively speed up parameter optimization.

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

[0204] The normalized Laguerre orthogonal polynomial is:

[0205] 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:

[0206] 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:

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

[0208] S5041: Acquire P time domain sampling points from the theoretical time domain electromagnetic data;

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

[0210] S5043: Obtain Q frequency domain sampling points from the theoretical frequency domain response data;

[0211] S5044: Based on the frequency domain sampling points and the frequency domain response function, obtain a discrete frequency domain response polynomial

[0212] S5045: 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 .

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

[0214] 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.

[0215] In one embodiment, Figure 3 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 :

[0216] S5046: Set the initial expansion item number N 0, an initial scale factor, wherein the initial scale factor includes an initial time domain scale factor and an initial frequency domain scale factor;

[0217] S5047: Based on the initial expansion item number N 0 、The initial scale factor calculation min{e(l 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:

[0218]

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

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

[0221] 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 as a new initial expansion item number, taking the optimal scale factor as a new initial scale factor, and repeating steps S5046, S5047 and S5048, where I is a positive integer (exemplarily, I=1);

[0222] 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.

[0223] 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.

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

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

[0226] 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.

[0227] In S5047, 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.

[0228] In S5048, 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.

[0229] 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.

[0230] In one embodiment, the present application provides a transient electromagnetic dual-parameter joint inversion device; Figure 4 As shown, the transient electromagnetic dual-parameter joint inversion device 30 may include an acquisition module 31 , a first modeling module 32 , a grid division module 33 , an assignment module 34 , a second modeling module 35 , a third modeling module 36 and an optimization module 37 .

[0231] The acquisition module 31 is used to acquire the actual electromagnetic data detected.

[0232] The first modeling module 32 is used to establish a geological model, wherein the geological model includes a resistivity parameter m ρ and the polarizability parameter m η .

[0233] The grid division module 33 is used to perform grid division on the geological model to obtain M grid units;

[0234] The assignment module 34 is used to assign parameters to each grid unit to form an initial parameter model m ap r, where is the initial resistivity parameter, is the initial polarizability parameter;

[0235] The second modeling module 35 is used to construct a theoretical mathematical model A(m) based on the geological model using a forward modeling method. The theoretical mathematical model includes the resistivity parameter m ρ and the polarizability parameter m η ;

[0236] The third modeling module 36 is used to construct the objective function P(m),

[0237] in:

[0238] m is the parameter to be solved, m=(m ρ , m η ) T ;

[0239] W d is a covariance matrix constructed based on the actual electromagnetic data;

[0240] d obs is an actual data matrix constructed based on the actual electromagnetic data;

[0241] λ is the regularization factor, λ=(λ ρ ,λ η ) T , where λ ρ is the resistivity regularization factor, λ η is the polarizability regularization factor;

[0242] W m is the least squares model weighting matrix;

[0243] is the fuzzy C-means clustering term, β is the clustering term weight factor;

[0244] The optimization module 37 is used to optimize the objective function to obtain a target geological model.

[0245] In one embodiment, the present application provides an inversion device; Figure 5 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 dual-parameter joint inversion method, 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 modes, such as Figure 4 The functions of modules 31 to 37 are shown.

[0246] 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 transient electromagnetic dual-parameter joint inversion method embodiment provided in the present application is implemented.

[0247] The embodiments described above 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 the technical solutions described in the aforementioned embodiments may still be modified, or some of the technical features may be replaced by equivalents. Such modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present application, and should all be included in the protection scope of the present application.

Claims

1. A transient electromagnetic dual-parameter joint inversion method, characterized in that: The following steps are involved: Acquire the actual electromagnetic data detected; Establish a geological model, the geological model includes a resistivity parameter m ρ and the polarizability parameter m η ; Dividing the geological model into grids to obtain M grid units; Assign parameters to each grid unit to form an initial parameter model m apr ,in, is the initial resistivity parameter, is the initial polarizability parameter; Based on the geological model, a theoretical mathematical model A(m) is constructed using a forward modeling method. The theoretical mathematical model includes the resistivity parameter m ρ and the polarizability parameter m η ; Construct the objective function P(m), in: m is the parameter to be solved, m=(m ρ , m η ) T ; W d is a covariance matrix constructed based on the actual electromagnetic data; d obs is an actual data matrix constructed based on the actual electromagnetic data; λ is the regularization factor, λ=(λ ρ ,λ η ) T , where λ ρ is the resistivity regularization factor, λ η is the polarizability regularization factor; W m is the least squares model weighting matrix; is the fuzzy C-means clustering term, β is the clustering term weight factor; The objective function is optimized to obtain a target geological model.

2. The transient electromagnetic dual-parameter joint inversion method according to claim 1 is characterized in that: in: M is the number of grid cells in the geological model; C is the number of cluster centers; m j is the parameter to be solved for the jth grid unit, m j =(m ρj ,m ηj ) T , where m ρj is the resistivity of the jth grid cell, m ηj is the polarizability of the jth grid cell; v k is the kth cluster center, v k =(v ρk ,v ηk ) T , where v ρk is the cluster center of the kth resistivity, v ηk is the cluster center of the kth polarizability; u jk is the degree of membership of the unsolved parameter of the jth grid unit to the kth cluster center; q is the fuzzification parameter; t k is the kth reference cluster center, t k =(t ρk ,t ηk ) T , where t ρk is the reference cluster center of the kth resistivity, t ηk is the reference cluster center of the kth polarizability; ω k is the kth reference cluster center t k The weight factor of the kth reference cluster center t k confidence level.

3. The transient electromagnetic dual-parameter joint inversion method according to claim 1 is characterized in that: Before optimizing the objective function, the following steps are included: Determine the resistivity parameter m of the node of each grid unit of the geological model ρ ;and According to the resistivity parameter m of the node ρ and interpolation function to calculate the resistivity parameter m at any point in each grid cell ρ , get the grid cells assigned by the parameters.

4. The transient electromagnetic dual-parameter joint inversion method according to any one of claims 1 to 3, characterized in that: The forward modeling method comprises: applying a forward excitation field to the geological model; In the time period from T0 to T2 after the forward 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; Based on the first forward mathematical model and Maxwell equations, the theoretical time-domain electromagnetic data and F in the time period from T0 to T1 are obtained. min Theoretical frequency domain electromagnetic data within the F1 frequency range, wherein both the theoretical time domain electromagnetic data and the theoretical frequency domain electromagnetic data contain induced polarization effect; Based on the theoretical frequency domain electromagnetic data, the theoretical time domain electromagnetic data and the second forward mathematical model, obtaining theoretical response electromagnetic data including induced polarization effect at any position in the geological model within a time period from T0 to T2; and The theoretical mathematical model A(m) is constructed according to the theoretical response electromagnetic data.

5. The transient electromagnetic dual-parameter joint inversion method according to claim 4 is characterized in that: The second forward mathematical model uses orthogonal polynomials as basis functions.

6. The transient electromagnetic dual-parameter joint inversion method according to claim 5 is characterized in that: The orthogonal polynomial is a Laguerre orthogonal polynomial, and the second forward modeling 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.

7. The transient electromagnetic dual-parameter joint inversion method according to claim 6 is characterized in that: The polynomial coefficient a n Obtained by the following steps: Acquire P time domain sampling points from the theoretical 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 theoretical 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 .

8. The transient electromagnetic dual-parameter joint inversion method according to claim 7 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: 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; 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 Step C: Calculate the expansion error and matching error based on the polynomial coefficients of M to N expansion terms, 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, 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.

9. A transient electromagnetic dual-parameter joint inversion device, characterized in that: include: An acquisition module, used for acquiring the actual electromagnetic data detected; The first modeling module is used to establish a geological model, wherein the geological model includes a resistivity parameter m ρ and the polarizability parameter m η ; A grid division module, used for performing grid division on the geological model to obtain M grid units; The assignment module is used to assign parameters to each grid unit to form an initial parameter model m apr ,in, is the initial resistivity parameter, is the initial polarizability parameter; The second modeling module is used to construct a theoretical mathematical model A(m) based on the geological model using a forward modeling method. The theoretical mathematical model includes the resistivity parameter m ρ and the polarizability parameter m η ; The third modeling module is used to construct the objective function P(m). in: m is the parameter to be solved, m=(m ρ , m η ) T ; W d is a covariance matrix constructed based on the actual electromagnetic data; d obs is an actual data matrix constructed based on the actual electromagnetic data; λ is the regularization factor, λ=(λ ρ ,λ η ) T , where λ ρ is the resistivity regularization factor, λ η is the polarizability regularization factor; W m is the least squares model weighting matrix; is the fuzzy C-means clustering term, β is the clustering term weight factor; The optimization module is used to optimize the objective function to obtain a target geological model.

10. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the transient electromagnetic dual-parameter joint inversion method described in any one of claims 1 to 8 is implemented.

Images