CN114779355A - Ground transient electromagnetic inversion method and device based on transmitting current full waveform - Google Patents

Abstract

The invention provides a ground transient electromagnetic inversion method and device based on the full waveform of the emission current. The method includes the following steps: constructing a one-dimensional inversion objective function based on a Tikhonov regularization method, and using the Gauss-Newton method for the inversion objective function. Perform optimization and output the final inversion model. In the process of optimizing the inversion objective function, the basic waveform of the emission current used for the one-dimensional inversion of the ground transient electromagnetic method is expanded from a sloped step waveform to include the on time, the settling time, the off time and the duty cycle. The full waveform of the emission current of the empty ratio, and the parameter of the number of cycles is included in the basic waveform of the emission current. This strategy not only helps to improve the accuracy of data processing and interpretation, but also helps to improve the design of instrument parameters.

Description

Ground Transient Electromagnetic Inversion Method and Device Based on Full Waveform of Transmitted Current

technical field

The invention relates to the field of geophysical detection, in particular to a ground transient electromagnetic method inversion method and device based on the full waveform of an emission current.

Background technique

The ground transient electromagnetic method mostly adopts the "bipolar" emission current waveform, that is, both positive and negative power supply in one emission cycle. As shown in Figure 1, both positive and negative supply half-cycles include current turn-on, settling, and turn-off phases, as well as a measurement phase. In the "bipolar" waveform, eddy currents are induced in the underground medium each time the emission current is turned on and turned off, and the eddy currents generated by the previous turn-on or turn-off of the emission current will inevitably partially offset the eddy currents generated by the next turn-off or turn-on. . Affected by the above phenomenon, there are differences in the magnitudes of the electromagnetic responses measured in the two half-cycles of the positive and negative power supply, respectively. In order to ensure that the absolute value of the electromagnetic response measured during two consecutive power supply half-cycles differs within a very small range, TEM usually does not start recording data until after several transmit cycles. In addition, in order to reduce the influence of random interference, the transient electromagnetic method data measurement usually adopts the multiple superposition technology, that is, the electromagnetic signal with the specified number of emission cycles needs to be continuously observed and superimposed. Therefore, the transient electromagnetic method observation data must be affected by the emission current waveform and the number of cycles.

At present, the ground transient electromagnetic method data processing software adopts the slope-step current waveform, that is, only considers the effect of the off time, while ignoring the effect of the settling time, the opening time and the number of cycles in the emission current waveform. However, under the background of high-conductivity exploration, the eddy current signal generated by the current turn-on is strong and decays slowly. Even if the stabilization time is as long as tens of milliseconds, it may not fully decay, and it will still partially offset the transient electromagnetic response after turn-off. Especially obvious. Therefore, the ground transient electromagnetic method data processing method that only considers the effect of the off time is not suitable for late data, and restricts the further improvement of data processing and interpretation accuracy.

SUMMARY OF THE INVENTION

The main technical problem solved by the present invention is that, on the basis of considering the influence of the turn-off time, the stabilization time and turn-on time in the emission current waveform, and the influence of the number of cycles are also considered, so as to improve the interpretation accuracy of the one-dimensional inversion of the ground transient electromagnetic method data. .

In order to achieve the above object, according to one aspect of the present invention, a ground transient electromagnetic method inversion method based on the full waveform of the emission current is provided, comprising the following steps:

设定为公式一：Step 1: Construct a one-dimensional inversion objective function based on the Tikhonov regularization method, the inversion objective function

Set to formula one:

where ||W _d (d ^obs -d ^pre )|| ² is the fitting difference function between the observed data and the forward model response, α||W _m (mm _* )|| ² is the model regularization function, | |·|| represents the l ₂ norm; among them, d ^obs is the observation data vector; d ^pre is the response vector of the forward model corresponding to the model parameter vector m, namely F(m), F represents the one-dimensional forward model of the transient electromagnetic method function; W _d is the data weighting matrix, here is a diagonal matrix, the elements of W _d are the reciprocal of the standard deviation of the observed data; α is the regularization factor; m _* is the reference model, which contains prior information; W _m is The model roughness operator, W _m is in the form:

Step 2: Use the Gauss-Newton method to optimize the inversion objective function, and output the final inversion model;

The process of optimizing the inversion objective function using the Gauss-Newton method includes:

Step S1: input the observation data and working parameters of the transient electromagnetic method;

Step S2: input the initial inversion model, reference model and control parameters;

Step S3: according to the working parameters, calculate the transient electromagnetic method forward response of the initial inversion model;

Step S3 includes:

S3.1: Calculate the magnetic induction intensity impulse response or step response excited by the lower step current waveform;

S3.2: According to the magnetic induction intensity impulse response or step response excited by the lower step current waveform, use the convolution algorithm to calculate the magnetic induction intensity impulse response or step response excited by the basic waveform of the emission current;

The basic waveform of the emission current is divided into two categories: the first type only includes the off time, that is, the ramp step waveform; the second type includes multiple cycles with different numbers and is determined by the on time, settling time, off time and duty cycle. The determined full waveform, the full waveform parameter depends on the actual current waveform, which is input by step S1;

Step S4: using the perturbation method to calculate the Jacobian matrix of the current inversion model (that is, the inversion model obtained by the k-th inversion iteration);

Step S5: generating a Gauss-Newton's normal equation according to the observation data, the reference model, the control parameters and the Jacobian matrix, and solving the Gauss-Newton's normal equation to obtain a model update amount;

Step S6: Update the current inversion model according to the model update amount, and calculate the forward response and root mean square fitting difference of the updated inversion model;

Step S7: Determine whether the current number of iterations or the root mean square fitting difference satisfies the iteration termination condition in the control parameter, if so, stop the iteration and output the final inversion model, otherwise, return to step S4 to start the next iteration.

Preferably, in step S1, the observation data includes: the magnetic induction intensity impulse response or step response normalized by the equivalent area of the receiving probe and the emission current intensity corresponding to the observation time series and the corresponding measurement standard deviation;

The working parameters include: working device, size and position of the emission loop, emission current waveform parameters, measurement point position and observation time sequence;

Wherein, the transmit current waveform parameters include: fundamental frequency, duty cycle, turn-on time, stabilization time and turn-off time.

Preferably, in step S2, the initial inversion model is a one-dimensional layered medium, and the input initial inversion model includes: inputting the number of layers of the layered medium, the thickness and conductivity of each layer of the medium;

Input prior information in the form of input reference model;

The control parameters include: target root mean square fitting error, maximum number of iterations, initial value of regularization factor, and update parameters.

Preferably, step S31 includes:

The rectangular emission loop source laid on the surface of the horizontal layered medium carries a current of a certain frequency f, and the component B _z (f) of the magnetic induction in the frequency domain excited at the measuring point (x, y, 0) along the z-axis direction is expressed as Formula two:

和

表达式依次为公式三至公式六：in,

and

The expressions are in order from Formula 3 to Formula 6:

其中k_x和k_y为傅里叶变换中的波数；J₁为1阶的第一类贝塞尔函数；r_L、r_-L、r_W和r_-W表示测点(x,y,0)到矩形发射回线四条边的距离，r_L、r_-L、r_W和r_-W的表达式分别为：r_L＝[(x′-x)²+(L-y)²]^1/2、r_-L＝[(x′-x)²+(-L-y)²]^1/2、r_W＝[(W-x)²+(y′-y)²]^1/2和r_-W＝[(-W-x)²+(y′-y)²]^1/2，其中四条边坐标依次表示为(x′,L,0)、(x′,-L,0)、(W,y′,0)和(-W,y′,0)；In the formula, E _x and E _y represent the electric field components along the x and y axes, respectively, I represents the current intensity, W and L represent the half-length and half-width of the rectangular emission loop source; z ₀ is the resistivity in vacuum, satisfying z ₀ =iωμ ₀ , i represents a complex unit, ω is the angular frequency, ω=2πf, μ ₀ is the magnetic permeability in vacuum; r _TE is the reflection coefficient, and the conductivity and thickness information of each layer are included in the reflection coefficient; the λ form for

where k _x and k _y are the wave numbers in the Fourier transform; J ₁ is the first-order Bessel function of the first order; r _L , r _-L , r _W and r _-W represent the measuring points (x, y, 0) The distance to the four sides of the rectangular emission return line, the expressions of r _L , r _-L , r _W and r _-W are respectively: r _L =[(x′-x) ² +(Ly) ² ] ^{1/ 2} , r _-L = [(x'-x) ² +(-Ly) ² ] ^1/2 , r _W =[(Wx) ² +(y'-y) ² ] ^1/2 and r _-W = [(-Wx) ² +(y′-y) ² ] ^1/2 , where the coordinates of the four sides are expressed as (x′,L,0), (x′,-L,0), (W,y′ in turn ,0) and (-W,y′,0);

The 12th-order Gauss-Legendre integral and the Hankel transform are used to solve the outer and inner integrals in Equation 3 to Equation 6 respectively, and obtain the magnetic induction intensity in the frequency domain;

或阶跃响应b_z(t)，表达式为公式七和八：Perform sine or cosine transformation on the imaginary part of the magnetic induction intensity in the frequency domain respectively to obtain the magnetic induction intensity impulse response excited by the lower step current waveform

or the step response b _z (t), expressed as Equations seven and eight:

Preferably, step S32 includes:

Convolve the calculated magnetic induction intensity impulse response or step response excited by the lower step current waveform with the current versus time differential in any basic waveform, and this process is expressed as formula 9:

表示电流对时间t的微分，τ表示积分变量；In the formula, f _waveform (t) represents the magnetic induction intensity impulse response or step response excited by the basic waveform, f _stepoff (t) represents the magnetic induction intensity impulse response or step response excited by the lower step current waveform,

represents the differential of the current with respect to time t, and τ represents the integral variable;

The trapezoidal integration method is used to solve Formula 9, and the magnetic induction intensity impulse response or step response excited by the basic waveform is obtained.

Preferably, in step S4, the disturbance amount is determined according to the parameters in the current inversion model m _k .

Preferably, in step S5, for the inversion objective function, the normal equation of the k+1 iteration of the Gauss-Newton method is Formula 10:

表示第k次反演模型m_k的正演响应，即F(m_k)；In the formula, δm _k+1 is the amount to be calculated, that is, the model update amount, and J _k is the first derivative of the forward function F at the k-th inversion model m _k , that is, the first derivative of the current inversion model m _k . Order derivative, namely Jacobian matrix, superscript T means transposition,

represents the forward response of the k-th inversion model m _k , namely F(m _k );

The preconditioned conjugate gradient method is used to solve Equation 10 to obtain δm _k+1 ; in the inversion, the regularization factor decreases proportionally with the number of inversion iterations.

Preferably, step S6 includes:

According to the model update amount, the current inversion model is updated by formula 11:

m _k+1 =m _k +aδm _k+1

In the formula, a is the model update step size, δm _k+1 is the model update amount, m _k is the inversion model obtained by the k-th inversion iteration, and m _k+1 is the inversion model obtained by the k+1-th inversion iteration. Inversion model, k represents the number of inversion iterations;

After the model is updated, the method in step S3 is used to calculate the TEM forward response of the updated inversion model, and the root mean square fitting error x _rms is calculated according to formula 12, and the root mean square fitting The difference is a measure of the fit of the data:

对应的标准差，

为反演模型m_k+1相应正演响应中的第m个数据。In the formula, M is the number of data, s _m is the mth observation data

The corresponding standard deviation,

is the mth data in the corresponding forward response of the inversion model m _k+1 .

Preferably, step S7 includes:

It is judged whether the root mean square fitting difference is less than or equal to the target fitting difference, or the set number of iterations is reached. If so, the iteration is stopped and the final inversion model is output; otherwise, the next iteration is started.

According to another aspect of the present invention, a ground transient electromagnetic method inversion device based on the full waveform of the emission current is provided, comprising:

设定为公式一：The inversion objective function building module is used to construct a one-dimensional inversion objective function based on the Tikhonov regularization method, and the inversion objective function

Set to formula one:

where ||W _d (d ^obs -d ^pre )|| ² is the fitting difference function between the observed data and the forward model response, α||W _m (mm _* )|| ² is the model regularization function, | |·|| represents the l ₂ norm; among them, d ^obs is the observation data vector; d ^pre is the response vector of the forward model corresponding to the model parameter vector m, namely F(m), F represents the one-dimensional forward model of the transient electromagnetic method function; W _d is the data weighting matrix, here is a diagonal matrix, the elements of W _d are the reciprocal of the standard deviation of the observed data; α is the regularization factor; m _* is the reference model, which contains prior information; W _m is The model roughness operator, W _m is in the form:

The inversion objective function optimization module is used to optimize the inversion objective function by using the Gauss-Newton method, and output the final inversion model;

The inversion objective function optimization module includes:

a first input unit, used for inputting the observation data and working parameters of the transient electromagnetic method;

a second input unit for inputting the initial inversion model, the reference model and the control parameters;

a forward response calculation unit, configured to calculate the transient electromagnetic method forward response of the initial inversion model according to the working parameters;

Wherein, the forward response calculation unit includes:

a first calculation unit, configured to calculate the magnetic induction intensity impulse response or step response excited by the lower step current waveform;

a second calculation unit, configured to calculate the magnetic induction intensity impulse response or step response excited by the basic waveform of the emission current by using a convolution algorithm according to the magnetic induction intensity impulse response or step response excited by the lower step current waveform;

The basic waveform of the emission current is divided into two categories: the first type only includes the off time, that is, the ramp step waveform; the second type includes multiple cycles with different numbers and is determined by the on time, settling time, off time and duty cycle. The determined full waveform, the full waveform parameter depends on the actual current waveform, and is input by the first input unit;

The Jacobian matrix calculation unit is used to calculate the Jacobian matrix of the current inversion model using the perturbation method;

A normal equation generating and solving unit, used for generating a Gauss-Newton normal equation according to the observation data, the reference model, the control parameters and the Jacobian matrix, and solving the Gauss-Newton normal equation to obtain a model update quantity;

an iterative update unit, configured to update the current inversion model according to the model update amount, and calculate the forward response and root mean square fitting difference of the updated inversion model;

The iteration termination unit is used to judge whether the current iteration number or the root mean square fitting difference satisfies the iteration termination condition in the control parameter, if so, stop the iteration and output the final inversion model, otherwise, return to the Jacobian The matrix calculation unit starts the next iteration.

The technical scheme provided by the invention has the following beneficial effects:

A ground transient electromagnetic method inversion method and device based on the full waveform of the emission current are provided, and the basic waveform of the emission current used for the ground transient electromagnetic method data processing is expanded from a sloped step waveform to include the opening time and the stabilization time. , off-time and duty cycle emission current full waveform, and the number of cycles parameter is included in the emission current basic waveform. This strategy not only helps to improve the accuracy of data processing and interpretation, but also helps to improve the design of instrument parameters.

Description of drawings

The present invention will be further described below in conjunction with the accompanying drawings and embodiments, in which:

Figure 1 is a schematic diagram of a single cycle "bipolar" waveform;

2 is a flowchart of a ground transient electromagnetic method inversion method based on the full waveform of the emission current in Embodiment 1 and Embodiment 2;

Fig. 3 is a schematic diagram of the basic waveform of the emission current in the first and second embodiments, wherein Fig. 3(a) waveform 1 (existing method), Fig. 3(b) waveform 2 (the present invention), Fig. 3(c) waveform 3 (the present invention), Figure 3 (d) waveform 4 (the present invention);

4 is a schematic diagram of a layered medium model in Embodiment 1 and Embodiment 2;

Fig. 5 is a geoelectric model obtained by one-dimensional inversion based on the four basic waveforms of Fig. 3 in Example 1, wherein Fig. 5(a) waveform 1, Fig. 5(b) waveform 2, Fig. 5(c) waveform 3, Figure 5(d) Waveform 4;

Fig. 6 is the field observation data (magnetic induction intensity step response) curve in embodiment 2;

Fig. 7 is the geoelectric model obtained based on the one-dimensional inversion of the four basic waveforms of Fig. 3 for the observation data of Fig. 6;

FIG. 8 is a structural diagram of a ground transient electromagnetic method inversion device based on the full waveform of the emission current provided by the present invention.

Detailed ways

In order to have a clearer understanding of the technical features, objects and effects of the present invention, the specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

Example 1: Theoretical Model

设定为：The one-dimensional inversion objective function is constructed based on the Tikhonov regularization method, and the Gauss-Newton method is used to optimize it. Inversion objective function

set as:

In the formula, the first term is the fitting difference function between the observed data and the forward model response, the second term is the model regularization function, and ||·|| represents the l ₂ norm. Among them, d ^obs is the observation data vector; d ^pre is the response vector of the forward model corresponding to the model parameter vector m, namely F(m), F represents the one-dimensional forward model function of the transient electromagnetic method; W _d is the data weighting matrix, this where is a diagonal matrix, the elements of W _d are the reciprocal of the standard deviation of the observed data; α is the regularization factor; m _* is the reference model, and the prior information is included. W _m is the model roughness operator, and its form is:

As shown in FIG. 3 , in the inversion of this embodiment, four basic waveforms of the transmit current are successively used, and 10 inversion calculations are performed for each waveform.

As shown in Figure 2, the process of each inversion calculation includes the following steps:

Step S1: input the observation data and working parameters of the transient electromagnetic method;

When the theoretical model is adopted, the observation data includes two parts: the forward model response of the specific geoelectric model and the working parameters (generated by the method in step S3) and random noise.

The geoelectric model and working parameters are as follows: the geoelectric model is a 5Ωm uniform half space; the emission source is a square loop with a side length of 100m, which is laid on the surface of the uniform half space; the device type is a center loop device; the emission current is 1A; the emission current The basic waveform is a bipolar current waveform with a fundamental frequency of 25 Hz (50% duty cycle) adopted by the Canadian PROTEM instrument. The corresponding time series contains 30 moments in the range of 6.8–6978 μs. Among them, the corresponding period of the 25Hz transmission fundamental frequency is 0.04s, that is, the period of 0.25 is 10ms, and the turn-on time, stabilization time and turn-off time are set to 1ms, 8.99ms and 0.01ms in turn. The forward model responds to the induced electromotive force (magnetic induction intensity impulse response) excited by the fundamental waveform of 8.25 cycles.

的0.5–3％(详见表1)。Random noise is generated as follows: Gaussian noise is first generated and then multiplied by the standard deviation as random noise. The standard deviation is set to the forward response at a certain moment

0.5–3% (see Table 1 for details).

The forward response is summed with random noise as observational data for inversion.

Table 1 Theoretical model forward modeling response standard deviation setting table

Step S2: input the initial inversion model, reference model and control parameters;

The initial inversion model is 39 layers of layered dielectrics, each layer of dielectric resistivity is 20Ωm, and the thickness is stepped from 5m to 40m (the 1st to 10th layers are 5m, the 11th to 25th layers are 10m, and the 26th to 35th layers are 20m, 40m for layers 36 to 38, and a uniform half space for layer 39), and the layer thickness is fixed in the inversion. No reference model. The target rms fit difference is 1, and the maximum number of iterations is 30. The initial value of the regularization factor α ₀ is 360000, and its update scheme is α _k+1 =α _k /1.5, where k is the number of inversion iterations.

Step S3: According to the transient electromagnetic method working parameters input in step S1, use the method described in step S3 to calculate the magnetic induction intensity impulse response under the excitation of the basic waveform

再采用卷积算法计算发射电流基本波形(见图3)激发的磁感应强度脉冲响应

图3中，图3(a)对应波形1为斜阶跃波形，被现有资料处理方法所采用，图3(b)、3(c)和3(d)为本发明使用的发射电流基本波形，分别对应波形2–4。图3(b)中，基本波形为正向供电发射电流波形，即0.25个周期；图3(c)中的基本波形在图3(b)所示波形的基础上，增加了一个周期，即1.25个周期；图3(d)中的基本波形在图3(b)所示波形的基础上，增加了两个周期，即2.25个周期。其发射电流波形参数(基频、占空比、开启时间、稳定时间和关断时间)由步骤S1中输入值确定。在一些实施例中，本发明的发射电流基本波形数量以及周期数目可以根据实际需要进行调整。Step S3 This step is the key to distinguish it from other inversion methods, that is, the influence of the full waveform of the emission current is included in the forward response calculation. Step S3 can be subdivided into two steps, first calculate the magnetic induction intensity impulse response excited by the next step current waveform

Then the convolution algorithm is used to calculate the magnetic induction intensity impulse response excited by the basic waveform of the emission current (see Figure 3).

In Fig. 3, the waveform 1 corresponding to Fig. 3(a) is a sloped step waveform, which is adopted by the existing data processing method. Figs. 3(b), 3(c) and 3(d) are the basic emission currents used in the present invention. waveforms, corresponding to waveforms 2–4, respectively. In Fig. 3(b), the basic waveform is the forward power emission current waveform, that is, 0.25 cycles; the basic waveform in Fig. 3(c) is based on the waveform shown in Fig. 3(b), and one cycle is added, that is, 1.25 cycles; the basic waveform in Fig. 3(d) adds two cycles to the waveform shown in Fig. 3(b), that is, 2.25 cycles. Its emission current waveform parameters (fundamental frequency, duty cycle, turn-on time, settling time and turn-off time) are determined by the input values in step S1. In some embodiments, the number of basic waveforms of the emission current and the number of cycles of the present invention can be adjusted according to actual needs.

The rectangular emission loop source laid on the surface of the horizontal layered medium carries a current of a certain frequency f, and the frequency domain magnetic induction intensity z component B _z (f) excited at the measuring point (x, y, 0) is expressed as:

和

表达式依次为：in,

and

The expressions are in turn:

其中k_x和k_y为傅里叶变换中的波数；J₁为1阶的第一类贝塞尔函数；r_L、r_-L、r_W和r_-W表示测点(x,y,0)到矩形发射回线四条边的距离，r_L、r_-L、r_W和r_-W的表达式为：In the formula, I represents the current intensity, and W and L represent the half-length and half-width of the rectangular emission loop. z ₀ =iωμ ₀ , i represents a complex unit, ω is the angular frequency, ω=2πf, μ ₀ is the magnetic permeability in vacuum. r _TE is the reflection coefficient, and the conductivity and thickness information of each layer are included in the reflection coefficient. The λ form is

where k _x and k _y are the wave numbers in the Fourier transform; J ₁ is the first-order Bessel function of the first order; r _L , r _-L , r _W and r _-W represent the measuring points (x, y, 0) The distance to the four sides of the rectangular emission return line, the expressions of r _L , r _-L , r _W and r _-W are:

r _L =[(x'-x) ² +(Ly) ² ] ^1/2 , r _-L =[(x'-x) ² +(-Ly) ² ] ^1/2 , r _W =[(Wx ) ² +(y′-y) ² ] ^1/2 and r _-W =[(-Wx) ² +(y′-y) ² ] ^1/2 , where the coordinates of the four sides are expressed as (x′,L ,0), (x',-L,0), (W,y',0) and (-W,y',0).

The 12th-order Gauss-Legendre integral and the Hankel transform are used to solve the outer and inner integrals in equations (4) to (7), respectively, and obtain the magnetic induction in the frequency domain.

表达式为：The sine transformation is carried out on the imaginary part of the magnetic induction intensity in the frequency domain respectively, and the magnetic induction intensity impulse response excited by the lower step current waveform is obtained.

The expression is:

Then the convolution algorithm is used to solve the electromagnetic field excited by the basic waveform of the emission current. In the convolution algorithm, first calculate the magnetic induction intensity impulse response (f _stepoff ) excited by the lower step current waveform according to formula (3), and then calculate the impulse response with the current versus time differential (dI/dt) in a certain basic waveform. Convolution operation, the magnetic induction intensity impulse response (f _waveform ) excited by the basic waveform can be obtained. This process can be represented by the following equation:

Equation (8) is solved by the trapezoidal integration method.

Step S4: Calculate the Jacobian matrix of the current inversion model by using the perturbation method. In this embodiment, there are 39 layers of dielectric. Since the layer thickness is fixed, the resistivity of each layer of dielectric needs to be perturbed one by one, and a complete Jacobian matrix can be obtained after perturbation 39 times. The disturbance amount is 10% of a certain parameter of the current inversion model. In some embodiments, the number of layers of the layered medium can be selected according to actual needs, and the corresponding disturbance times and disturbance amounts can also be adjusted accordingly.

Step S5: Substitute into steps S3 and S4 to generate the forward response and Jacobian matrix respectively, and then substitute the current regularization factor to generate the Gauss-Newton normal equation. The normal equation is solved by the preconditioned conjugate gradient method, and the termination condition is that the number of iterations reaches 200 or the relative residual is less than or equal to 10 ^-6 .

Step S6: Calculate the model update amount and the root mean square fitting difference, where the model update step size is set to 1.

Step S7: judging the termination condition. In this theoretical model study, the RMS fitting errors were all less than the target value within 10 inversion iterations.

The above inversion procedure was run 10 times. The observational data used for the inversion contains random noise, so the geoelectric model obtained from each inversion is slightly different. The inversion results are shown in Figure 5, where the solid line represents the theoretical model and the dashed line represents the inversion model. When waveform 1 (the existing method) is used, most of the inversion geoelectric models are significantly different from the theoretical model (see Figure 5 ( a)); when waveforms 2, 3 and 4 (the present invention) are used, most of the inversion geoelectric models agree with the theoretical models (see Figures 5(b), 5(c) and 5(d)). Therefore, compared with the existing data processing methods based on the sloped step current waveform, the ground transient electromagnetic method inversion method based on the full waveform of the emission current proposed by the present invention has obvious advantages and achieves the expected goal.

Example 2: Field observation data

Step S1: input the observation data and working parameters of the transient electromagnetic method;

The observation data collection site is a tidal flat in the East China Sea. This area is in an alternating land and sea environment. The shallow strata are generally saturated with seawater, and the resistivity is only 0.1 to 0.8Ωm, and the resistivity of the deep strata increases slightly. The observed data is the step response of the magnetic induction intensity of a measuring point, which is collected by a superconducting quantum interferometer (SQUID) produced by Crone Instruments in Canada. The emission source is a rectangular loop with a side length of 400m×200m, the loop center coordinate is (0,0,0)m, and the measurement point coordinate is (-100,0,0)m. The emission current is 25A; the basic waveform of the emission current is the bipolar current waveform (50% duty cycle) of the emission fundamental frequency of 0.833333Hz adopted by the Canadian Crone Pulse EM instrument, and the corresponding time series includes 45 moments with a range of -0.15ms –222.7ms. Among them, the corresponding period of the 0.833333Hz transmission fundamental frequency is 1.2s, that is, the period of 0.25 is 300ms, and the turn-on time, stabilization time and turn-off time are set to 0.5ms, 299.5ms and 1ms in turn. In the inversion, only the observation data from time 2 to time 45 are used, see Figure 6. Since the observation data collected by the Crone Pulse EM instrument does not contain standard deviation, the present invention matches the standard deviation for the measured data according to the characteristics of the transient electromagnetic method with strong early signal and high signal-to-noise ratio, and weak late signal and low signal-to-noise ratio. See Table 2.

Table 2 Standard deviation setting table of field observation data

Step S2: input the initial inversion model, reference model and control parameters;

The initial inversion model is 46 layers of layered dielectrics, each layer of dielectric resistivity is 1Ωm, and the thickness is stepped from 10m to 40m (the first to 20th layers are 10m, the 21st to 30th layers are 20m, and the 31st to 45th layers are 40m, the 46th layer medium is a uniform half space), and the layer thickness is fixed in the inversion. No reference model. The target rms fit difference is 1, and the maximum number of iterations is 30. The initial value of the regularization factor α ₀ is 360000, and its update scheme is α _k+1 =α _k /1.5, where k is the number of inversion iterations.

Step S3: According to the transient electromagnetic method working parameters input in step S1, the method described in step S3 is used to calculate the magnetic induction intensity step response b _z (t) under the excitation of the basic waveform.

Step S4: Calculate the Jacobian matrix of the current inversion model by using the perturbation method. In this embodiment, there are 46 layers of medium. Since the layer thickness is fixed, the resistivity of each layer of medium needs to be perturbed one by one, and a complete Jacobian matrix can be obtained after perturbation 46 times.

Step S5: Substitute the forward response and Jacobian matrix generated in step S3 and step S4 respectively, and then substitute the current regularization factor to generate the Gauss-Newton normal equation. The normal equation is solved by the preconditioned conjugate gradient method, and the termination condition is that the number of iterations reaches 200 or the relative residual is less than or equal to 10 ^-6 .

Step S6: Calculate the model update amount and the root mean square fitting difference, where the model update step size is set to 1.

Step S7: judging the termination condition. In this theoretical model study, the RMS fitting errors were all less than the target value within 10 inversion iterations.

The inversion results are shown in Fig. 7. The inversion results obtained by waveform 1 (existing method) and the inversion results obtained by waveforms 2, 3 and 4 (the present invention) are significantly different from 100 to 1000 m underground, and the The inversion results obtained by waveforms 2, 3 and 4 (the present invention) are more in line with the real geological conditions. Therefore, the ground transient electromagnetic method inversion method based on the full waveform of the emission current proposed by the present invention achieves the expected goal.

As shown in FIG. 8 , this embodiment also provides a ground transient electromagnetic method inversion device based on the full waveform of the emission current, including:

设定为公式一：The inversion objective function building module 01 is used to construct a one-dimensional inversion objective function based on the Tikhonov regularization method, and the inversion objective function

Set to formula one:

where ||W _d (d ^obs -d ^pre )|| ² is the fitting difference function between the observed data and the forward model response, α||W _m (mm _* )|| ² is the model regularization function, | |·|| represents the l ₂ norm; among them, d ^obs is the observation data vector; d ^pre is the response vector of the forward model corresponding to the model parameter vector m, namely F(m), F represents the one-dimensional forward model of the transient electromagnetic method function; W _d is the data weighting matrix, here is a diagonal matrix, the elements of W _d are the reciprocal of the standard deviation of the observed data; α is the regularization factor; m _* is the reference model, which contains prior information; W _m is The model roughness operator, W _m is in the form:

The inversion objective function optimization module 02 is used to optimize the inversion objective function by using the Gauss-Newton method, and output the final inversion model;

The inversion objective function optimization module 02 includes:

The first input unit 021, for inputting the observation data and working parameters of the transient electromagnetic method;

The second input unit 022 is used for inputting the initial inversion model, the reference model and the control parameters;

The forward response calculation unit 023 is used to calculate the transient electromagnetic method forward response of the initial inversion model according to the input working parameters;

The forward response calculation unit 023 includes:

The first calculation unit 0231 is used to calculate the magnetic induction intensity impulse response or step response excited by the lower step current waveform;

The second calculation unit 0232 is configured to use a convolution algorithm to calculate the magnetic induction intensity impulse response or step response excited by the basic waveform of the emission current according to the magnetic induction intensity impulse response or step response excited by the lower step current waveform;

The basic waveform of the emission current is divided into two categories: the first type only includes the off time, that is, the ramp step waveform; the second type includes multiple cycles with different numbers and determined by the on time, settling time, off time and duty cycle. Full waveform, the full waveform parameter depends on the actual current waveform, which is input by the first input unit 021;

Jacobian matrix calculation unit 024, used for calculating the Jacobian matrix of the current inversion model by the perturbation method;

The normal equation generation and solving unit 025 is used to generate the Gauss-Newton's normal equation according to the input observation data, the reference model, the control parameter and the calculated Jacobian matrix, and solve the Gauss-Newton's normal equation to obtain the model update amount;

an iterative update unit 026, configured to update the current inversion model according to the model update amount, and calculate the forward response and the root mean square fitting difference of the updated inversion model;

The iteration termination unit 027 is used to judge whether the current iteration number or the root mean square fitting difference satisfies the iteration termination condition in the control parameter, if so, stop the iteration and output the final inversion model, otherwise, return the Jacobian matrix The computing unit 024 starts the next iteration.

It should be noted that, herein, the terms "comprising", "comprising" or any other variation thereof are intended to encompass non-exclusive inclusion, such that a process, method, article or system comprising a series of elements includes not only those elements, It also includes other elements not expressly listed or inherent to such a process, method, article or system. Without further limitation, an element qualified by the phrase "comprising a..." does not preclude the presence of additional identical elements in the process, method, article or system that includes the element.

The above-mentioned serial numbers of the embodiments of the present invention are only for description, and do not represent the advantages or disadvantages of the embodiments. In a unit claim enumerating several means, several of these means may be embodied by one and the same item of hardware. The use of the words first, second, and third, etc. do not denote any order, and these words may be construed as identifications.

The above are only preferred embodiments of the present invention, and are not intended to limit the scope of the present invention. Any equivalent structure or equivalent process transformation made by using the contents of the description and drawings of the present invention, or directly or indirectly applied in other related technical fields , are similarly included in the scope of patent protection of the present invention.

Claims (10)

1. a ground transient electromagnetic method inversion method based on the full waveform of emission current, is characterized in that, comprises the following steps: Step 1: Construct a one-dimensional inversion objective function based on the Tikhonov regularization method, the inversion objective function

Set to formula one:

where ||W _d (d ^obs -d ^pre )|| ² is the fitting difference function between the observed data and the forward model response, α||W _m (mm _* )|| ² is the model regularization function, | |·|| represents the l ₂ norm; among them, d ^obs is the observation data vector; d ^pre is the response vector of the forward model corresponding to the model parameter vector m, namely F(m), F represents the one-dimensional forward model of the transient electromagnetic method function; W _d is the data weighting matrix, here is a diagonal matrix, the elements of W _d are the reciprocal of the standard deviation of the observed data; α is the regularization factor; m _* is the reference model, which contains prior information; W _m is The model roughness operator, W _m is in the form:

Step 2: using the Gauss-Newton method to optimize the inversion objective function, and output the final inversion model; The process of optimizing the inversion objective function using the Gauss-Newton method includes: Step S1: input the observation data and working parameters of the transient electromagnetic method; Step S2: input the initial inversion model, reference model and control parameters; Step S3: according to the working parameters, calculate the transient electromagnetic method forward response of the initial inversion model; Step S3 includes: S3.1: Calculate the magnetic induction intensity impulse response or step response excited by the lower step current waveform; S3.2: According to the magnetic induction intensity impulse response or step response excited by the lower step current waveform, use the convolution algorithm to calculate the magnetic induction intensity impulse response or step response excited by the basic waveform of the emission current; The basic waveform of the emission current is divided into two categories: the first type only includes the off time, that is, the ramp step waveform; the second type includes multiple cycles with different numbers and is determined by the on time, settling time, off time and duty cycle. The determined full waveform, the full waveform parameter depends on the actual current waveform, which is input by step S1; Step S4: using the perturbation method to calculate the Jacobian matrix of the current inversion model; Step S5: generating a Gauss-Newton's normal equation according to the observation data, the reference model, the control parameters and the Jacobian matrix, and solving the Gauss-Newton's normal equation to obtain a model update amount; Step S6: Update the current inversion model according to the model update amount, and calculate the forward response and root mean square fitting difference of the updated inversion model; Step S7: Determine whether the current number of iterations or the root mean square fitting difference satisfies the iteration termination condition in the control parameter, if so, stop the iteration and output the final inversion model, otherwise, return to step S4 to start the next iteration. 2. The ground transient electromagnetic method inversion method based on the full waveform of emission current as claimed in claim 1, is characterized in that, in step S1, described observation data comprises: the equivalent area of receiving probe corresponding to observation time series and and Emission current intensity normalized magnetic induction intensity impulse response or step response and corresponding measurement standard deviation; The working parameters include: working device, size and position of the emission loop, emission current waveform parameters, measurement point position and observation time sequence; Wherein, the transmit current waveform parameters include: fundamental frequency, duty cycle, turn-on time, stabilization time and turn-off time. 3. The ground transient electromagnetic method inversion method based on the full waveform of the emission current according to claim 1, wherein in step S2, the initial inversion model is a one-dimensional layered medium, and the input initial inversion model is The modeling model includes: the number of input layered dielectric layers, the thickness and conductivity of each layer of dielectric; Input prior information in the form of input reference model; The control parameters include: target root mean square fitting error, maximum number of iterations, initial value of regularization factor, and update parameters. 4. The ground transient electromagnetic method inversion method based on the full waveform of emission current as claimed in claim 1, is characterized in that, step S31 comprises: The rectangular emission loop source laid on the surface of the horizontal layered medium carries a current of a certain frequency f, and the component B _z (f) of the magnetic induction in the frequency domain excited at the measuring point (x, y, 0) along the z-axis direction is expressed as Formula two:

in,

and

The expressions are in order from Formula 3 to Formula 6:

In the formula, E _x and E _y represent the electric field components along the x and y axes, respectively, I represents the current intensity, W and L represent the half-length and half-width of the rectangular emission loop source; z ₀ is the resistivity in vacuum, satisfying z ₀ =iωμ ₀ , i represents a complex unit, ω is the angular frequency, ω=2πf, μ ₀ is the magnetic permeability in vacuum; r _TE is the reflection coefficient, and the conductivity and thickness information of each layer are included in the reflection coefficient; the λ form for

where k _x and k _y are the wave numbers in the Fourier transform; J ₁ is the first-order Bessel function of the first order; r _L , r _-L , r _W and r _-W represent the measuring points (x, y, 0) The distance to the four sides of the rectangular emission return line, the expressions of r _L , r _-L , r _W and r _-W are: r _L =[(x'-x) ² +(Ly) ² ] ^1/2 , r _-L =[(x'-x) ² +(-Ly) ² ] ^1/2 , r _W =[(Wx ) ² +(y′-y) ² ] ^1/2 and r _-W =[(-Wx) ² +(y′-y) ² ] ^1/2 , where the coordinates of the four sides are expressed as (x′,L ,0), (x',-L,0), (W,y',0) and (-W,y',0); The 12th-order Gauss-Legendre integral and the Hankel transform are used to solve the outer and inner integrals in Equation 3 to Equation 6 respectively, and obtain the magnetic induction intensity in the frequency domain; Perform sine or cosine transformation on the imaginary part of the magnetic induction intensity in the frequency domain respectively to obtain the magnetic induction intensity impulse response excited by the lower step current waveform

or the step response b _z (t), expressed as Equations seven and eight:

5. The ground transient electromagnetic method inversion method based on the full waveform of emission current as claimed in claim 1, is characterized in that, step S32 comprises: Convolve the calculated magnetic induction intensity impulse response or step response excited by the lower step current waveform with the current versus time differential in any basic waveform, and this process is expressed as formula 9:

In the formula, f _waveform (t) represents the magnetic induction intensity impulse response or step response excited by the basic waveform, f _stepoff (t) represents the magnetic induction intensity impulse response or step response excited by the lower step current waveform,

represents the differential of the current with respect to time t, and τ represents the integral variable; The trapezoidal integration method is used to solve Formula 9, and the magnetic induction intensity impulse response or step response excited by the basic waveform is obtained. 6 . The ground transient electromagnetic method inversion method based on the full waveform of the emission current according to claim 1 , wherein, in step S4 , the disturbance amount is determined according to the parameters in the current inversion model m _k . 7 . 7. The ground transient electromagnetic method inversion method based on the full waveform of emission current as claimed in claim 1, is characterized in that, in step S5, for inversion objective function, the normal equation of Gauss-Newton method k+1 iteration For formula ten:

In the formula, δm _k+1 is the amount to be calculated, that is, the model update amount, and J _k is the first derivative of the forward function F at the k-th inversion model m _k , that is, the first derivative of the current inversion model m _k . Order derivative, namely Jacobian matrix, superscript T means transposition,

represents the forward response of the k-th inversion model m _k , namely F(m _k ); The preconditioned conjugate gradient method is used to solve Equation 10 to obtain δm _k+1 ; in the inversion, the regularization factor decreases proportionally with the number of inversion iterations. 8. The ground transient electromagnetic method inversion method based on the full waveform of emission current as claimed in claim 1, is characterized in that, step S6 comprises: According to the model update amount, the current inversion model is updated by formula 11: m _k+1 =m _k +aδm _k+1 In the formula, a is the model update step size, δm _k+1 is the model update amount, m _k is the inversion model obtained by the k-th inversion iteration, and m _k+1 is the inversion model obtained by the k+1-th inversion iteration. Inversion model, k represents the number of inversion iterations; After the model is updated, the method in step S3 is used to calculate the TEM forward response of the updated inversion model, and the root mean square fitting error x _rms is calculated according to formula 12, and the root mean square fitting The difference is a measure of the fit of the data:

In the formula, M is the number of data, s _m is the mth observation data

The corresponding standard deviation,

is the mth data in the corresponding forward response of the inversion model m _k+1 . 9. The ground transient electromagnetic method inversion method based on the full waveform of emission current as claimed in claim 1, is characterized in that, step S7 comprises: Determine whether the root mean square fitting difference is less than or equal to the target fitting difference, or reach the set number of iterations, if so, stop the iteration and output the final inversion model; otherwise, start the next iteration. 10. A ground transient electromagnetic method inversion device based on the full waveform of emission current, characterized in that, comprising: The inversion objective function building module is used to construct a one-dimensional inversion objective function based on the Tikhonov regularization method, and the inversion objective function

Set to formula one:

where ||W _d (d ^obs -d ^pre )|| ² is the fitting difference function between the observed data and the forward model response, α||W _m (mm _* )|| ² is the model regularization function, | |·|| represents the l ₂ norm; among them, d ^obs is the observation data vector; d ^pre is the response vector of the forward model corresponding to the model parameter vector m, namely F(m), F represents the one-dimensional forward model of the transient electromagnetic method function; W _d is the data weighting matrix, here is a diagonal matrix, the elements of W _d are the reciprocal of the standard deviation of the observed data; α is the regularization factor; m _* is the reference model, which contains prior information; W _m is The model roughness operator, W _m is in the form:

The inversion objective function optimization module is used to optimize the inversion objective function by using the Gauss-Newton method, and output the final inversion model; The inversion objective function optimization module includes: a first input unit, used for inputting the observation data and working parameters of the transient electromagnetic method; a second input unit for inputting the initial inversion model, the reference model and the control parameters; a forward response calculation unit, configured to calculate the transient electromagnetic method forward response of the initial inversion model according to the working parameters; Wherein, the forward response calculation unit includes: a first calculation unit, configured to calculate the magnetic induction intensity impulse response or step response excited by the lower step current waveform; a second calculation unit, configured to calculate the magnetic induction intensity impulse response or step response excited by the basic waveform of the emission current by using a convolution algorithm according to the magnetic induction intensity impulse response or step response excited by the lower step current waveform; The basic waveform of the emission current is divided into two categories: the first type only includes the off time, that is, the ramp step waveform; the second type includes multiple cycles with different numbers and is determined by the on time, settling time, off time and duty cycle. The determined full waveform, the full waveform parameter depends on the actual current waveform, and is input by the first input unit; The Jacobian matrix calculation unit is used to calculate the Jacobian matrix of the current inversion model using the perturbation method; A normal equation generating and solving unit, used for generating a Gauss-Newton normal equation according to the observation data, the reference model, the control parameters and the Jacobian matrix, and solving the Gauss-Newton normal equation to obtain a model update quantity; an iterative update unit, configured to update the current inversion model according to the model update amount, and calculate the forward response and root mean square fitting difference of the updated inversion model; The iteration termination unit is used to judge whether the current iteration number or the root mean square fitting difference satisfies the iteration termination condition in the control parameter, if so, stop the iteration and output the final inversion model, otherwise, return to the Jacobian The matrix calculation unit starts the next iteration.

Images