CN106646625A - Transient electromagnetic inversion method for sharp boundary model - Google Patents

Abstract

The invention discloses a transient electromagnetic inversion method based on a sharp boundary model, comprising the following steps: inputting secondary field data observed by a transient electromagnetic method; inputting inversion parameters and constructing an initial inversion model; Mismatch functional between forward modeling theoretical data; construct minimum gradient support functional that constrains resistivity change; construct inversion objective function from mismatch functional and minimum gradient support functional; take initial model as initial condition, adopt Gaussian The Newton method solves the corresponding resistivity when the objective function is minimized, and obtains the inversion result. The invention can improve the resolution of the abrupt change surface of underground resistivity, and makes up for the defect that Occam inversion does not reflect clearly on the electrical interface. It is of great significance for the performance and later data interpretation.

Description

A Transient Electromagnetic Inversion Method for Sharp Boundary Models

technical field

The invention relates to the technical field of geophysical detection, in particular to a transient electromagnetic inversion method of a sharp boundary model.

Background technique

As a powerful means of geophysical exploration, transient electromagnetic method can solve various geophysical problems such as mineral exploration, energy, engineering, hydrology, environmental geological survey, and archaeological exploration. It uses step waves or other pulse currents to emit a pulsed magnetic field underground, and the underground conductive geological body generates an induced eddy current under the excitation of the primary field, and then stimulates a secondary magnetic field. When the emission current is turned off, the primary field disappears, and the subsurface eddy current decays with time, and its decay time is related to the electrical parameters (volume, structure, resistivity, buried depth) of the conductive geological body. By observing the secondary field value, the electrical parameters of the underground heterogeneous body can be reconstructed. This process is called inversion. On the contrary, the process of calculating the change of the secondary field using the underground structure with known electrical parameters is called forward modeling.

The commonly used TEM one-dimensional inversion methods include Marquardt method and Occam method. The Marquardt method directly fits the forward modeling data of the hypothetical model with the measured data, and iteratively calculates the thickness and resistivity of each layer of the earth model. This method has fast calculation speed and simple principle, and the result can reflect the abrupt boundary of earth resistivity, but the inversion process is unstable, and it is easy to fall into local optimal solution. The Occam method imposes a longitudinal smoothing constraint on the basis of the Marquardt method constraint function to minimize the two-norm of the longitudinal resistivity gradient, adopts a fixed thickness of the earth model, and uses the Gauss-Newton method to iteratively solve the resistivity of each layer of the earth. This method is stable in the inversion process and can find a global optimal solution that meets the longitudinal smoothness constraints. However, the resolution of some geoelectric interfaces with sudden changes in resistivity is not good, and the layered structure of the earth cannot be accurately described.

Contents of the invention

(1) Technical problems to be solved

The invention provides a transient electromagnetic inversion method of a sharp boundary model, which is used to solve the unstable inversion of the Marquardt method, the Occam method to smooth the change of the earth resistivity as much as possible, and the layered boundary resolution of the sudden change of the resistivity. However, it is impossible to accurately describe the layered structure of the earth and the combination of the above two methods is not convenient enough in practical engineering applications.

(2) Technical solutions

The present invention provides a transient electromagnetic inversion method for a sharp boundary model. The inversion is to reconstruct the electrical parameters of the inhomogeneous body through the measured secondary field data. The transient electromagnetic inversion method includes constructing the inversion initial The model is characterized in that it also includes the steps of:

Constructing a minimum gradient support functional according to the electrical parameters to be inverted of the inversion initial model; and

Construct the inversion objective function according to the minimum gradient support functional.

The above scheme also includes:

Forward modeling the inversion initial model to obtain a theoretical response value;

A mismatch functional is constructed according to the measured secondary field data and the theoretical response value.

In the above scheme, the construction of the inversion objective function also includes:

An inversion objective function is constructed according to the mismatched functional.

In the above scheme, the minimum gradient support functional is constructed using the following formula:

In the formula, P _s is the minimum gradient support functional, β is called the focusing factor, which is a positive number much smaller than 1, m is the electrical parameter to be inverted, and R is the matrix describing the gradient change: N is the number of inversion layers.

In the above scheme, the inversion objective function is constructed using the following formula:

P=P _t +P _s

In the formula, P is the inversion objective function, P _t is the mismatch functional, α is a regularization factor, which is used to adjust the proportion between the mismatch functional and the minimum gradient support functional during the inversion process, and P _s is the minimum gradient support functional.

The above scheme also includes:

Taking the initial model as the initial condition, the electrical parameters corresponding to the minimum value of the inversion objective function are solved to obtain the inversion result.

(3) Beneficial effects

The method provided by the invention can obtain the following beneficial effects:

1. The method of the present invention inherits the advantages of stable Occam inversion by using the regularized objective function, and the inversion result does not depend on the initial model, and the global optimal solution can be found.

2. The method of the present invention selects the minimum gradient support functional as the resistivity constraint condition, produces a result similar to Marquardt inversion, can highlight the stratum boundary with sudden change in resistivity, and makes up for the defect that Occam inversion does not reflect clearly on the electrical interface.

3. The method of the present invention provides specific formulas for each step in the process of use, is easy to operate, and can realize the inversion results of mixed use of Occam inversion and Marquardt inversion in one step.

Description of drawings

Fig. 1 schematically shows a flowchart of a transient electromagnetic inversion method for a sharp boundary model according to an embodiment of the present invention.

Fig. 2 schematically shows corresponding data of a transient electromagnetic secondary field according to an embodiment of the present invention.

Fig. 3 schematically shows the inversion result of a sharp boundary model according to an embodiment of the present invention.

detailed description

In order to make the object, technical solution and advantages of the present invention clearer, the present invention will be further described in detail below in conjunction with specific embodiments and with reference to the accompanying drawings.

Fig. 1 schematically shows a flowchart of a transient electromagnetic inversion method for a sharp boundary model according to an embodiment of the present invention. as the picture shows,

In step S1, the inversion parameters are input, and the inversion parameters include transmitting current I, radius R of the transmitting coil, number of turns N _t of the transmitting coil, observation offset r, and effective area S of the receiving coil. According to the embodiment of the present invention, the number of turns of the transmitting coil is N _t =1, the transmitting current is I=5A, the transmitting radius R=50m, the offset distance of the observation point is r=0m, and the effective area of the receiving coil is 2000m ² .

In step S2, the inversion initial model is constructed. The inversion initial model includes the inversion layer number N, the thickness of each layer h=[h ₁ , h ₂ ,..., h _N-1 ], and the inversion initial resistivity m= [m ₁ , m ₂ , . . . , m _N ]. Since the inversion resolution decreases with depth, h _i+1 /h _i >1 is usually set. According to the embodiment of the present invention, the inversion model is 30 layers, the thickness increase ratio h _i+1 /h _i =1.12, the thickness of the first layer h ₁ =5m, and the initial resistivity of each layer is 100Ωm.

In step S3, the initial model is forward modeled to obtain theoretical data. According to an embodiment of the present invention, the initial inversion model is forward modeled to obtain a theoretical response F(m).

In step S4, input the measured secondary field data, input the secondary field data observed by the transient electromagnetic method, including the acquisition time T and its corresponding field value d _obs . Among them, the field value d _obs needs to satisfy the attenuation characteristics, and the part containing the primary field in the early stage and the part oscillating due to noise in the late stage need to be eliminated. According to the embodiment of the present invention, as shown in FIG. 2 , the data is obtained by forward modeling of a three-layer earth model, the thicknesses of which are 100 m and 200 m respectively, and a uniform half-space. The corresponding resistivities are 300Ωm, 90Ωm, 300Ωm, respectively.

In step S5, a mismatch functional between the measured secondary field data and the theoretical response value is constructed. According to the embodiment of the present invention, the mismatch functional can be written as:

P _t ＝||W _d (F(m)-d _obs )|| ²

In the formula, W _d is the weight matrix describing the proportion of data, F(m) is the theoretical response value, and W _d can be selected as:

In the formula, σ _i is the noise of each time measurement channel, and M is the total number of time measurement channels.

According to the embodiment of the present invention, it is assumed that the proportions of trace data at each time are consistent, that is, the weight matrix is a unit diagonal matrix.

In step S6, the minimum gradient support functional is constructed. According to the embodiment of the present invention, constructing the minimum gradient support functional that constrains the resistivity change can be written as:

In the formula, β is called the focusing factor, which is a positive number much smaller than 1. R is a matrix describing the gradient change, which can be written as:

According to the embodiment of the present invention, in order to highlight the abrupt boundary of the earth resistivity, β should be as small as possible, and β=5×10 ⁻¹¹ is set.

In step S7, the inversion objective function is constructed from the mismatch functional and the minimum gradient support functional. According to the embodiment of the present invention, the inversion objective function can be written as:

P=P _t +αP _s

In the formula, α is called a regularization factor, which is used to adjust the proportion between the mismatch functional and the minimum gradient support functional during the inversion process. According to the embodiment of the present invention, the search range of α is set to 10 ^-4 -10 ⁰ .

In step S8, taking the initial model as the initial condition, the Gauss-Newton method is used to solve the corresponding resistivity when the objective function is minimized, and the inversion result is obtained. The theoretical response value F(m) is approximately written as:

F(m)＝F(m ₀ )+J ₀ (m ₁ -m ₀ )

In the formula, m ₀ is the initial condition of a certain iteration, m ₁ is the resistivity to be obtained in this iteration, J ₀ is the derivative matrix of the forward function to the initial condition in a certain iteration, each element of which is:

Substitute F(m)=F(m ₀ )+J ₀ (m ₁ -m ₀ ) into the objective function P=P _t +αP _s , let Get the iterative formula:

After each iteration, update the initial conditions with the newly obtained results, when the fit is poor The iteration terminates when less than the expected fit difference. Each iterative process needs to dynamically select the α value within a certain interval to minimize the objective function. The embodiments of the present invention select but are not limited to the Gauss-Newton method.

In step S9, the final inversion result is obtained, that is, the electrical parameters of the inhomogeneous body. Mutation boundaries have higher resolution.

The specific embodiments described above have further described the purpose, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above descriptions are only specific embodiments of the present invention, and are not intended to limit the present invention. Within the spirit and principles of the present invention, any modifications, equivalent replacements, improvements, etc., shall be included in the protection scope of the present invention.

Claims (6)

1. A transient electromagnetic inversion method of a sharp boundary model, said inversion is to reconstruct the electrical parameters of an inhomogeneous body by measuring secondary field data, and the transient electromagnetic inversion method comprises constructing an inversion initial model, It is characterized in that it also includes the steps of: Constructing a minimum gradient support functional according to the electrical parameters to be inverted of the inversion initial model; and Construct the inversion objective function according to the minimum gradient support functional. 2. the transient electromagnetic inversion method of the sharp boundary model according to claim 1, is characterized in that, also comprises: Forward modeling the inversion initial model to obtain a theoretical response value; A mismatch functional is constructed according to the measured secondary field data and the theoretical response value. 3. the transient electromagnetic inversion method of sharp boundary model according to claim 1, is characterized in that, described construction inversion objective function also comprises: An inversion objective function is constructed according to the mismatched functional. 4. the transient electromagnetic inversion method of the sharp boundary model according to claim 1, is characterized in that, use following formula to construct described minimum gradient support functional as: ${\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\frac{\mathrm{<span\; class="notranslate">\; R</span>}\mathrm{<span\; class="notranslate">\; m</span>}}{{<span\; class="notranslate">\; \&lsqb;</span>\mathrm{<span\; class="notranslate">\; R</span>}\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; \&Center\; Dot;</span>\mathrm{<span\; class="notranslate">\; R</span>}\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&beta;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&rsqb;</span>}^{\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; |</span>{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$ In the formula, P _s is the minimum gradient support functional, β is called the focusing factor, which is a positive number much smaller than 1, m is the electrical parameter to be inverted, and R is the matrix describing the gradient change: N is the number of inversion layers. 5. the transient electromagnetic inversion method of sharp boundary model according to claim 1, is characterized in that, use following formula to construct described inversion objective function as: P=P _t +αP _s In the formula, P is the inversion objective function, P _t is the mismatch functional, α is a regularization factor, which is used to adjust the proportion between the mismatch functional and the minimum gradient support functional during the inversion process, and P _s is the minimum gradient support functional. 6. the transient electromagnetic inversion method of the sharp boundary model according to claim 1, is characterized in that, also comprises: Taking the initial model as the initial condition, the electrical parameters corresponding to the minimum value of the inversion objective function are solved to obtain the inversion result.