A method and system for extracting induced polarization parameters from magnetotelluric observation data

Abstract

The present invention relates to an induced polarization parameter extraction method for magnetotelluric observation data, comprising: assigning initial values ​​to zero-frequency resistivity, polarizability, time constant and frequency correlation coefficient respectively to obtain an initial geoelectric model; establishing an objective function according to the standard error term between electromagnetic observation data and electromagnetic theoretical data and the focusing constraint condition of the geoelectric model; calculating the optimal damping factor according to the fitting difference between the electromagnetic observation data and the electromagnetic theoretical data; calculating the Hessian matrix according to the apparent resistivity and impedance phase in the electromagnetic theoretical data; calculating the correction amount according to the objective function, the optimal damping factor and the Hessian matrix; correcting the geoelectric model according to the correction amount to obtain an updated geoelectric model; judging whether the updated geoelectric model meets the preset iteration condition; if not, continuing to update the geoelectric model; if satisfied, obtaining induced polarization parameters extracted based on electromagnetic observation data. Extracting induced polarization parameters directly based on observation data makes the parameters more real and reliable.

Description

Technical Field

[0001] The present invention relates to the field of geophysical exploration technology, and in particular to a method and system for extracting induced polarization parameters from magnetotelluric observation data. Background Art

[0002] In geophysical exploration, a single exploration method can often only obtain one physical property description of the underground medium. This is also true in electrical exploration, but there are exceptions. The complex resistivity method (also called induced polarization method) uses the Cole-Cole model to describe the true resistivity of the medium, and can obtain four apparent parameters describing the physical properties of the underground medium. However, due to the mutual influence between the parameters and the weak forward response, they are generally obtained from active source exploration. The extraction of true induced polarization parameters from natural sources has progressed slowly. At present, it is still at the stage of obtaining the true induced polarization parameters from electric and magnetic fields rather than observational data, and is still far from practical application. In addition, due to the existence of multiple solutions in geophysical inversion methods, the extraction of multiple true parameters will reduce the reliability of the results. Therefore, it is very necessary to carry out research on stable extraction technology directly targeting observational data to accelerate the process of practical application. Summary of the invention

[0003] The purpose of the present invention is to provide a method and system for extracting induced polarization parameters from magnetotelluric observation data to solve the dilemma of directly fitting the electric and magnetic fields instead of directly fitting the first-hand field observation data, so as to make the extracted induced polarization parameters more realistic and more reliable.

[0004] To achieve the above object, the present invention provides the following solutions:

[0005] A method for extracting induced polarization parameters from magnetotelluric observation data, the method comprising:

[0006] Assigning initial values ​​to the zero-frequency resistivity, polarizability, time constant and frequency correlation coefficient respectively to obtain the initial geoelectric model; each of the initial values ​​belongs to a value within a preset range corresponding to each parameter; and making the initial geoelectric model the current geoelectric model;

[0007] Establishing a current objective function according to a standard error term between electromagnetic observation data and current electromagnetic theoretical data and a focusing constraint condition of a current geoelectric model; the current electromagnetic theoretical data is electromagnetic theoretical data corresponding to the current geoelectric model;

[0008] Calculating a current optimal damping factor according to a fitting difference between the electromagnetic observation data and the electromagnetic theoretical data;

[0009] Calculate the current Hessian matrix according to the apparent resistivity and impedance phase in the electromagnetic theory data;

[0010] Calculating a correction amount of the current geoelectric model according to the current objective function, the current optimal damping factor and the current Hessian matrix;

[0011] Correcting the current geoelectric model according to the correction amount to obtain a current updated geoelectric model;

[0012] Determining whether the currently updated geoelectric model meets a preset iteration condition;

[0013] If not, the updated geoelectric model is set as the current geoelectric model, and the process returns to the step of "establishing the current objective function according to the standard error term between the electromagnetic observation data and the electromagnetic theoretical data and the focusing constraint of the current geoelectric model";

[0014] If satisfied, the current updated geoelectric model is output to obtain induced polarization parameters extracted based on the electromagnetic observation data.

[0015] A system for extracting induced polarization parameters from magnetotelluric observation data, comprising:

[0016] An initial geoelectric model construction module is used to assign initial values ​​to zero-frequency resistivity, polarizability, time constant and frequency correlation coefficient to obtain the initial geoelectric model; each of the initial values ​​is within a preset range corresponding to each parameter; and the initial geoelectric model is set as the current geoelectric model;

[0017] An objective function building module is used to establish a current objective function according to a standard error term between electromagnetic observation data and current electromagnetic theoretical data and a focusing constraint condition of a current geoelectric model; the current electromagnetic theoretical data is electromagnetic theoretical data corresponding to the current geoelectric model;

[0018] An optimal damping factor calculation module, used for calculating the current optimal damping factor according to the fitting difference between the electromagnetic observation data and the electromagnetic theoretical data;

[0019] A Hessian matrix calculation module, used for calculating the current Hessian matrix according to the apparent resistivity and impedance phase in the electromagnetic theory data;

[0020] A model correction amount calculation module, used for calculating the correction amount of the current geoelectric model according to the current objective function, the current optimal damping factor and the current Hessian matrix;

[0021] A model updating module, used for correcting the current geoelectric model according to the correction amount to obtain a current updated geoelectric model;

[0022] A judgment module, used to judge whether the current updated geoelectric model meets the preset iteration condition;

[0023] If not, the current updated geoelectric model is set as the current geoelectric model, and the objective function construction module is executed;

[0024] If satisfied, the current updated geoelectric model is output to obtain induced polarization parameters extracted based on the electromagnetic observation data.

[0025] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0026] The present invention relates to an induced polarization parameter extraction method and system for magnetotelluric observation data, the method comprising: assigning initial values ​​to zero-frequency resistivity, polarizability, time constant and frequency correlation coefficient respectively to obtain an initial geoelectric model; establishing an objective function according to the standard error term between electromagnetic observation data and electromagnetic theoretical data and the focusing constraint condition of the geoelectric model; calculating the optimal damping factor according to the fitting difference between the electromagnetic observation data and the electromagnetic theoretical data; calculating the Hessian matrix according to the apparent resistivity and impedance phase in the electromagnetic theoretical data; calculating the correction amount according to the objective function, the optimal damping factor and the Hessian matrix; correcting the geoelectric model according to the correction amount to obtain an updated geoelectric model; judging whether the updated geoelectric model meets the preset iteration condition; if not, continuing to update the geoelectric model; if satisfied, obtaining the induced polarization parameters extracted based on the electromagnetic observation data. Extracting induced polarization parameters directly based on electromagnetic observation data makes the extracted induced polarization parameters more real and reliable. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0028] Figure 1 A flow chart of a method for extracting induced polarization parameters from magnetotelluric observation data provided in Example 1 of the present invention;

[0029] Figure 2 A schematic diagram of the induced polarization parameter results in the example of the verification extraction method provided in Example 1 of the present invention;

[0030] Figure 3 A block diagram of an induced polarization parameter extraction system for magnetotelluric observation data provided in Example 2 of the present invention. DETAILED DESCRIPTION

[0031] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0032] In geophysics, underground rocks and minerals are abstracted into various attribute bodies (density body, velocity body, resistivity body, magnetic body, etc.), and in electrical prospecting, underground media are considered to be resistivity bodies. Electromagnetic observation data can be collected in the field through electrical prospecting instruments. The mathematical calculation process of obtaining the underground electrical structure characteristics of the corresponding data collection area from the electromagnetic observation data is the inversion problem in electrical prospecting. How to achieve inversion specifically? The common practice is to first assume that the values ​​of the underground electrical body model (also known as the geoelectric model) are consistent. At this time, the initial geoelectric model is called the initial geoelectric model (for convenience, "geoelectric" is often omitted and referred to as the initial model, the same below). According to the characteristics of the electromagnetic field in geophysics, when the initial model is given, the forward response characteristics of the geoelectric model (that is, the theoretical data F[m0]) can be calculated, and the error (fitting error) between the theoretical data and the electromagnetic observation data is calculated to see if it meets the requirements. If it does not meet the requirements, the initial model is modified according to certain rules, and the error between the theoretical data of the modified geoelectric model and the electromagnetic observation data is recalculated, and the cycle is repeated until the error meets the requirements. At this point, the geoelectric model modified for the last time is considered to be the result model.

[0033] The purpose of the present invention is to provide a method and system for extracting induced polarization parameters from magnetotelluric observation data to solve the dilemma of directly fitting the electric and magnetic fields instead of directly fitting the first-hand field observation data, so as to make the extracted induced polarization parameters more realistic and more reliable.

[0034] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0035] Example 1

[0036] like Figure 1 As shown, this embodiment provides a method for extracting induced polarization parameters from magnetotelluric observation data, the method comprising:

[0037] Step S1: assigning initial values ​​to induced polarization parameters such as zero-frequency resistivity, polarizability, time constant and frequency correlation coefficient to obtain the initial geoelectric model; each of the initial values ​​belongs to a value within a preset range corresponding to each parameter; and the initial geoelectric model is set as the current geoelectric model.

[0038] Inversion is actually an iterative approximation mathematical process. To obtain the resulting geoelectric model, it is necessary to set the initial model in advance, and then modify the initial model again and again according to a certain rule until the preset conditions are met. The geoelectric model after the last modification is considered to be the resulting model. In this embodiment, the zero-frequency resistivity ρ0∈(10 -2 ~10 5 ), polarizability η∈(0~0.98), time constant τ∈(10 -3 ~5×10 3 ), frequency correlation coefficient c∈(0.1~0.6). Random values ​​are taken within the preset range of each parameter to form the initial geoelectric model, that is, the geoelectric model is recorded as a vector composed of four parameters: zero-frequency resistivity, polarizability, time constant and frequency correlation coefficient.

[0039] Step S2: Establish the current objective function according to the standard error term between the electromagnetic observation data and the current electromagnetic theoretical data and the focusing constraint of the current geoelectric model; the current electromagnetic theoretical data is the electromagnetic theoretical data corresponding to the current geoelectric model; the geoelectric model is often used in this field to calculate the electromagnetic theoretical data through forward modeling methods such as finite element method, finite difference method, integral equation method, etc.

[0040] In reality, the electrical characteristics of underground rocks and minerals rarely undergo large mutations, and the resistivity value changes relatively slowly and smoothly. Therefore, this gentle and smooth change can be used as a limiting means to constrain the geoelectric model, and it is believed that as long as the resistivity value changes suddenly, the resulting model is unreliable. How to abstract and theorize this limitation and constraint? The common practice is to introduce a new parameter (model roughness) that changes with the resistivity value of the model. In order to make the inversion method more applicable, a compromise parameter (called: damping factor) is added between the model roughness and the fitting error during the inversion process. The objective function is established under the idea of ​​"finding the solution with the minimum model roughness under a certain fitting error standard". The independent variable of the objective function is the geoelectric model parameter, and its general formula is as follows:

[0041] U[m]=a[m]+λ -1 b[m]

[0042] Among them: a[m] is the model roughness term, b[m] is the standard error term of electromagnetic observation data and electromagnetic theoretical data, λ is the damping factor, and m is the geoelectric model parameter. Substituting the specific data items into the general formula, the complete expression of the objective function can be obtained, as follows:

[0043]

[0044] ||▽m|| 2 is the model roughness term, ▽m is the gradient of the model parameters; is the standard error term, ||Wd-WF(m)|| is the data fitting error; X * is the expected value of the fitting error, d is the observed data (apparent resistivity and impedance phase), that is, b[m]; F[m] is the response data (electromagnetic theory data, also apparent resistivity and impedance phase), and W is the element of 1 / σ j A diagonal matrix, where σ j is the standard deviation of the jth observation data measured by the instrument.

[0045] In fact, model roughness is a constraint on the geoelectric model. This constraint is based on the known information that the electrical characteristics of underground rocks and minerals rarely undergo drastic mutations, and is therefore also called prior constraint information. This prior constraint information can play different roles, and the role of model roughness is only to reduce the number of unreliable geoelectric models. In view of the problem that the range of resistivity body anomalies in the inversion results is too large based on the electromagnetic field response law of electrical exploration, and the fact that the geoelectric model parameter in the conventional electrical exploration inversion problem is only one (resistivity), while the inversion induced polarization parameters include four (zero-frequency resistivity, polarizability, time constant, frequency correlation coefficient), and this multi-parameter inversion significantly increases the instability of the solution process. The present invention proposes to use a focused constraint method to reduce the number of unreliable solutions, narrow the range of anomalies, and improve the stability of the objective function solution process.

[0046] Therefore, the expression of the objective function of the geoelectric model obtained in the i-th iteration is:

[0047]

[0048] Among them, U[m i ] is the objective function of the geoelectric model obtained in the i-th iteration; m i is the geoelectric model obtained in the i-th iteration; It is the minimum support functional of the prior constraints, which is used to stabilize the entire target functional so that the volume of the abnormal body in the model parameters is minimized; ▽m i is the gradient of the geoelectric model obtained at the i-th iteration; β is a non-zero decimal, V is the volume domain that constitutes the geoelectric model parameters; λ i is the damping factor; is the standard error term, ||Wd-WF(m i )|| is the data fitting error; X * is the expected value of the fitting difference between electromagnetic observation data and electromagnetic theoretical data; d is the electromagnetic observation data, F[m i ] is the electromagnetic theoretical data corresponding to the geoelectric model obtained in the i-th iteration, and W is a matrix normalized by the standard deviation of the data, that is, W is a matrix with elements of 1 / σ j A diagonal matrix, where σ jis the standard deviation of the jth electromagnetic observation data measured by the instrument.

[0049] In this way, a focusing constraint is successfully added to the objective function, which can improve the resolution and stability of the results.

[0050] Step S3: Calculate the current optimal damping factor according to the fitting difference between the electromagnetic observation data and the electromagnetic theoretical data.

[0051] The damping factor is a compromise parameter between the model roughness term and the data fitting error term. The optimal damping factor needs to be determined by searching. The search process needs to meet two principles: one is to reduce the data fitting error to the expected value of the fitting error, and the other is to keep the target fitting error at the expected value while searching for the model that minimizes the objective function. At this time, the corresponding damping factor is the optimal damping factor λ. * ;

[0052] The expression for calculating the optimal damping factor during the i-th iteration is:

[0053] X(λ 2i )=min||Wd-WF(m i (λ 2i ))||

[0054]

[0055] Among them, λ i * is the optimal damping factor for the i-th iteration; λ 2i is the threshold of the optimal damping factor for the i-th iteration.

[0056] Step S4: Calculate the current Hessian matrix according to the apparent resistivity and impedance phase in the electromagnetic theory data; specifically, obtain the current Hessian matrix by taking partial derivatives of the current apparent resistivity and the current impedance phase respectively;

[0057] Forward response F[m i ] is the geoelectric model parameter m i A nonlinear function of the previous model parameter m i-1 The model parameter m can be obtained by performing a first-order Taylor series expansion linearization at i-1 The Hessian matrix A i .

[0058] F[m i ]=F[m i-1 ]+A(m i -m i-1 )+ε, ε is an infinitesimal quantity;

[0059]

[0060] Since the actual electromagnetic observation data and forward response data are both apparent resistivity ρ s and impedance phase φ, the model parameters are zero-frequency resistivity ρ0, polarizability η, time constant τ, and frequency correlation coefficient c. Therefore, the Hessian matrix obtained in the i-th iteration is:

[0061]

[0062]

[0063] In the above formula, The calculation formula is:

[0064]

[0065] Where ω and μ are the angular frequency and the dielectric constant in vacuum, respectively; Z i is the electromagnetic wave impedance at the i-th iteration, Re() and Im() are the real and imaginary parts respectively, ρ i is the complex resistivity calculated by the single Cole-Cole model, σ i is i The reciprocal of m i is the geoelectric model parameter vector [ρ 0i ,η i ,τ i ,c i ], since the geoelectric model is composed of four parameters, we can also directly understand m i is the geoelectric model under the i-th iteration; and The calculation of resistivity can refer to the method of extracting resistivity in the field.

[0066] Step S5: Calculate the correction amount of the current geoelectric model according to the current objective function, the current optimal damping factor and the current Hessian matrix.

[0067] Step S51: deforming the minimum support functional in the current objective function to obtain the current objective penalty function.

[0068] Specifically, the minimum support functional is deformed to obtain a first deformation expression; the first deformation expression is:

[0069] make The second deformation expression is obtained; the second deformation expression is: S βMSi =(w βMSi (m i -m i-1 ),w βMSi (mi -m i-1 ));

[0070] Introduce the integrated sensitivity diagonal matrix S and combine w βMSi Get S βMSi The third deformation expression of W; the third deformation expression is: βSi =w βMSi ·S; where w βMSi is a diagonal matrix.

[0071] The target penalty function is obtained according to the current target function, the first deformation expression, the second deformation expression and the third deformation expression; the target penalty function is:

[0072]

[0073] Step S52: Calculate the partial derivative of the objective penalty function and introduce the current optimal damping factor and the current Hessian matrix to calculate the correction amount of the current geoelectric model. By solving , we can get the model correction amount;

[0074] The correction amount of the geoelectric model obtained in the i-th iteration is expressed as:

[0075]

[0076] Among them, Δm i ' is the correction amount obtained in the i-th iteration; Δd i =dF(m i )+A i m i , A i is the Hessian matrix obtained at the i-th iteration; here is the weight coefficient, which is related to the focusing constraint; I is the unit matrix.

[0077] Step S6: correcting the current geoelectric model according to the correction amount to obtain a current updated geoelectric model;

[0078] That is m i =m i-1 +Δm′ i ; Where m0 is the initial geoelectric model;

[0079] Step S7: determining whether the currently updated geoelectric model meets a preset iteration condition;

[0080] If not, the currently updated geoelectric model is set as the current geoelectric model, and the process returns to step S2 "establishing the current objective function according to the standard error term between the electromagnetic observation data and the electromagnetic theoretical data and the focusing constraint condition of the current geoelectric model";

[0081] If satisfied, the current updated geoelectric model is output to obtain induced polarization parameters extracted based on the electromagnetic observation data.

[0082] The determining whether the currently updated geoelectric model meets the preset iteration condition specifically includes:

[0083] It is determined whether the updated values ​​of the zero-frequency resistivity, the polarizability, the time constant and the frequency-correlation coefficient in the currently updated geoelectric model are within the corresponding preset ranges to obtain a first judgment result; limiting the range to the preset range can effectively reduce results that have no physical meaning.

[0084] If the first judgment result is no, the zero-frequency resistivity, the polarizability, the time constant and the frequency correlation coefficient whose updated values ​​are not within the preset range are reassigned to obtain the current modified geoelectric model, and the current modified geoelectric model is set as the current geoelectric model and returns to the step of "establishing the current objective function according to the standard error term between the electromagnetic observation data and the electromagnetic theoretical data and the focusing constraint of the current geoelectric model". The reassignment here can be manually assigning random values ​​within the preset range, or taking the value obtained in the previous iteration process that is within the preset range.

[0085] If the first judgment result is yes, then judging whether the standard deviation (fitting error) of the electromagnetic theoretical data and the electromagnetic observation data corresponding to the current updated geoelectric model reaches a threshold value to obtain a second judgment result, or judging whether the current number of corrections reaches a set maximum number of corrections to obtain a third judgment result;

[0086] If both the second judgment result and the third judgment result are negative, the current updated geoelectric model is set as the current geoelectric model and the process returns to the step of "establishing the current objective function according to the standard error term between the electromagnetic observation data and the electromagnetic theoretical data and the focusing constraint condition of the current geoelectric model";

[0087] If at least one of the second judgment result and the third judgment result is yes, the current updated geoelectric model is output to obtain induced polarization parameters extracted based on the electromagnetic observation data.

[0088] In this embodiment, the induced polarization parameter focusing extraction technology method for magnetotelluric observation data provided by the present invention can directly use the magnetotelluric data collected outdoors. This method is no longer limited to the fitting of theoretical electric and magnetic fields, and has taken a very important step towards the practical application of multi-property parameter extraction in the field of electromagnetic methods for geophysical exploration. At the same time, a focusing constraint model is added to the objective function, so that problems of sharp boundaries can be better solved, and the inversion stability and the resolution of abnormal bodies can be effectively improved. This technical method comprehensively considers the electromagnetic coupling effect of electromagnetic waves, makes full use of electromagnetic data at different frequencies, and obtains the resistivity and electrochemical parameters of underground media. While searching for conventional geological targets, polymetallic ores associated with metal sulfides and other ore bodies with high polarizability can also be discovered, providing a favorable means for solving resource exploration.

[0089] The following is a verification of the IP parameter focusing extraction technique for magnetotelluric observation data provided in this embodiment.

[0090] Design a double prism low-resistance high-polarization model. All parameters of the anomalies on both sides are consistent, and the specific induced polarization parameters are zero-frequency resistivity of 10Ω·m, polarizability of 0.5, time constant of 50, and frequency correlation coefficient of 0.5. The distance between the two prisms is 2000m, the burial depth is 1000m, and the length and width are both 1000m×300m; buried in a uniform half-space with zero-frequency resistivity of 100Ω·m, polarizability of 0.5, time constant of 5, and frequency correlation coefficient of 0.2. The magnetotelluric data includes the apparent resistivity and impedance phase of the two polarization modes of TM and TE. The data contains 21 observation points with a spacing of 50m between the observation points; there are 21 observation frequencies in total, ranging from 10000Hz to 0.1Hz. The number of grid divisions is 50×80, the air layer expansion factor is 1.5, the expansion layer is 10 layers, and the thickness of the first layer of grid is the initial thickness; the number of left and right expansion layers is 8 layers, the expansion factor is 0.8, and the maximum vertical depth is 80km.

[0091] Figure 2The results of induced polarization parameter extraction (inversion results). From the zero-frequency resistivity results, it can be seen that the geometric shapes of the two low-resistance and high-polarization anomalies have been well restored; but due to the TE mode effect, the longitudinal elongation of the anomalies is large, making the overall result resemble a "heart". In contrast, the inversion results of the polarizability show that the geometric shape of the anomalies is more consistent with the real model, and there is no obvious longitudinal stretching effect. The overall inversion effect of the time constant is better than that of the zero-frequency resistivity, and slightly worse than that of the polarizability; but there are problems with the inversion values. The real time constant of the anomaly is 50, and the values ​​of the entire result are too small. At the same time, the frequency correlation coefficient does not have a situation where the value is overall too small; the spatial distribution of the anomaly is also basically restored. In summary, it can be found that the focus constraint extraction technology for magnetotelluric observation data can effectively obtain the induced polarization parameters of the underground medium, and the multi-attribute physical property parameters have restored the geometric shape of the anomaly.

[0092] Example 2

[0093] like Figure 3 As shown, in this embodiment, a system for extracting induced polarization parameters for magnetotelluric observation data is provided, comprising:

[0094] The initial geoelectric model construction module M1 is used to assign initial values ​​to the zero-frequency resistivity, polarizability, time constant and frequency-related coefficient to obtain the initial geoelectric model; each of the initial values ​​is within the preset range corresponding to each parameter; and the initial geoelectric model is set as the current geoelectric model;

[0095] The objective function building module M2 is used to establish the current objective function according to the standard error term between the electromagnetic observation data and the current electromagnetic theoretical data and the focusing constraint condition of the current geoelectric model; the current electromagnetic theoretical data is the electromagnetic theoretical data corresponding to the current geoelectric model;

[0096] An optimal damping factor calculation module M3, used to calculate the current optimal damping factor according to the fitting difference between the electromagnetic observation data and the electromagnetic theoretical data;

[0097] A Hessian matrix calculation module M4, used for calculating the current Hessian matrix according to the apparent resistivity and impedance phase in the electromagnetic theory data;

[0098] A model correction amount calculation module M5 is used to calculate the correction amount of the current geoelectric model according to the current objective function, the current optimal damping factor and the current Hessian matrix;

[0099] A model updating module M6, used for correcting the current geoelectric model according to the correction amount to obtain a current updated geoelectric model;

[0100] A judgment module M7 is used to judge whether the current updated geoelectric model meets the preset iteration condition;

[0101] If not, the current updated geoelectric model is set as the current geoelectric model, and the objective function construction module is executed;

[0102] If satisfied, the current updated geoelectric model is output to obtain induced polarization parameters extracted based on the electromagnetic observation data.

[0103] In this specification, each embodiment is described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the embodiments can be referred to each other. For the system disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant parts can be referred to the method part.

[0104] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A method for extracting induced polarization parameters from magnetotelluric observation data, characterized in that: The method comprises: Assigning initial values ​​to zero-frequency resistivity, polarizability, time constant and frequency correlation coefficient respectively to obtain an initial geoelectric model; each of the initial values ​​belongs to a value within a preset range corresponding to each parameter; and making the initial geoelectric model the current geoelectric model; Establishing a current objective function according to a standard error term between electromagnetic observation data and current electromagnetic theoretical data and a focusing constraint condition of a current geoelectric model; the current electromagnetic theoretical data is electromagnetic theoretical data corresponding to the current geoelectric model; Calculating a current optimal damping factor according to a fitting difference between the electromagnetic observation data and the electromagnetic theoretical data; Calculate the current Hessian matrix according to the apparent resistivity and impedance phase in the electromagnetic theory data; Calculating a correction amount of the current geoelectric model according to the current objective function, the current optimal damping factor and the current Hessian matrix; Correcting the current geoelectric model according to the correction amount to obtain a current updated geoelectric model; Determining whether the currently updated geoelectric model meets a preset iteration condition; If not satisfied, the current updated geoelectric model is set as the current geoelectric model, and the process returns to step "establishing the current objective function according to the standard error term between the electromagnetic observation data and the current electromagnetic theoretical data and the focusing constraint condition of the current geoelectric model"; If satisfied, the current updated geoelectric model is output to obtain induced polarization parameters extracted based on the electromagnetic observation data.

2. The method according to claim 1, characterized in that The expression of the current objective function is established based on the standard error term between the electromagnetic observation data and the current electromagnetic theoretical data and the focusing constraints of the current geoelectric model: Among them, U[m i ] is the objective function of the geoelectric model obtained in the i-th iteration; m i is the geoelectric model obtained in the i-th iteration; is the minimum support functional of the prior constraints; is the gradient of the geoelectric model obtained at the i-th iteration; β is a non-zero decimal, V is the volume domain that constitutes the geoelectric model parameters; λ i is the damping factor; is the standard error term, ||Wd-WF(m i )|| is the data fitting error; X * is the expected value of the fitting difference between electromagnetic observation data and electromagnetic theoretical data; d is the electromagnetic observation data; F(m i ) is the electromagnetic theoretical data corresponding to the geoelectric model obtained in the ith iteration; W is the matrix normalized by the standard deviation of the data.

3. The method according to claim 2, characterized in that The expression for calculating the optimal damping factor is: X(λ i )=min||Wd-WF(m i (l 2i ))|| l i * =max{λ i |X(λ i )=X * ,l i >l 2i } Among them, λ i * is the optimal damping factor for the i-th iteration; λ 2i is the threshold of the optimal damping factor for the i-th iteration.

4. The method according to claim 2, characterized in that: The calculating of the current Hessian matrix according to the apparent resistivity and impedance phase in the electromagnetic theory data specifically includes: The partial derivatives of the current apparent resistivity and the current impedance phase are respectively obtained to obtain the current Hessian matrix; Among them, the expressions of the apparent resistivity and the impedance phase obtained in the i-th iteration are: The expression of the Hessian matrix obtained in the i-th iteration is: In the formula, ρ si and φ i Respectively represent the apparent resistivity and impedance phase; ω, μ are the angular frequency and the dielectric constant in vacuum; Z i is the electromagnetic wave impedance at the i-th iteration, Re() and Im() are the real and imaginary parts respectively.

5. The method according to claim 2, characterized in that: The calculating the correction amount of the current geoelectric model according to the current objective function, the current optimal damping factor and the current Hessian matrix specifically includes: Transforming the minimum support functional in the current objective function to obtain the current objective penalty function; The partial derivative of the aforementioned objective penalty function is obtained and the current optimal damping factor and the current Hessian matrix are introduced to calculate the correction amount of the current geoelectric model.

6. The method according to claim 5, characterized in that The deformation of the minimum support functional in the current objective function to obtain the current objective penalty function specifically includes: The minimum support functional is deformed to obtain the first deformation expression; the first deformation expression is: m i-1 represents the geoelectric model obtained at the i-1th iteration; make The second deformation expression is obtained; the second deformation expression is: S βMSi =(w βMSi (m i -m i-1 ),w βMSi (m i -m i-1 )); Introduce the integrated sensitivity diagonal matrix S and combine w βMSi Get S βMSi The third deformation expression of W; the third deformation expression is: βSi =w βMSi ·S; The target penalty function is obtained according to the current target function, the first deformation expression, the second deformation expression and the third deformation expression; the target penalty function is: represents the target penalty function; λ i * is the optimal damping factor for the i-th iteration.

7. The method according to claim 6, characterized in that The expression for calculating the correction amount of the current geoelectric model is: Among them, Δm i ' is the correction amount obtained in the i-th iteration; Δd=dF(m0)+A0m0, A0 is the initial Hessian matrix; A i is the Hessian matrix obtained at the i-th iteration; is the weight coefficient, I is the unit matrix; W is the matrix normalized by the standard deviation of the data; d is the electromagnetic observation data; m0 is the initial geoelectric model; F(m0) is the electromagnetic theoretical data corresponding to the initial geoelectric model.

8. The method according to claim 1, characterized in that The determining whether the currently updated geoelectric model meets the preset iteration condition specifically includes: Determine whether the updated values ​​of the zero-frequency resistivity, the polarizability, the time constant and the frequency correlation coefficient in the currently updated geoelectric model are within the corresponding preset ranges to obtain a first determination result; If the first judgment result is no, the zero-frequency resistivity, the polarizability, the time constant and the frequency correlation coefficient whose updated values ​​are not within the preset range are reassigned to obtain the current modified geoelectric model, the current modified geoelectric model is set as the current geoelectric model and the step of "establishing the current objective function according to the standard error term between the electromagnetic observation data and the electromagnetic theoretical data and the focusing constraint condition of the current geoelectric model" is returned to; If the first judgment result is yes, then judging whether the standard deviation of the electromagnetic theoretical data and the electromagnetic observation data corresponding to the current updated geoelectric model reaches a threshold value to obtain a second judgment result, or judging whether the current number of corrections reaches a set maximum number of corrections to obtain a third judgment result; If both the second judgment result and the third judgment result are negative, the current updated geoelectric model is set as the current geoelectric model and the process returns to step "establishing the current objective function according to the standard error term between the electromagnetic observation data and the current electromagnetic theoretical data and the focusing constraint condition of the current geoelectric model"; If at least one of the second judgment result and the third judgment result is yes, the current updated geoelectric model is output to obtain induced polarization parameters extracted based on the electromagnetic observation data.

9. An induced polarization parameter extraction system for magnetotelluric observation data based on the method according to any one of claims 1 to 8, characterized in that: include: An initial geoelectric model construction module is used to assign initial values ​​to zero-frequency resistivity, polarizability, time constant and frequency correlation coefficient to obtain an initial geoelectric model; each of the initial values ​​is within a preset range corresponding to each parameter; and the initial geoelectric model is set as the current geoelectric model; An objective function building module is used to establish a current objective function according to a standard error term between electromagnetic observation data and current electromagnetic theoretical data and a focusing constraint condition of a current geoelectric model; the current electromagnetic theoretical data is electromagnetic theoretical data corresponding to the current geoelectric model; An optimal damping factor calculation module, used for calculating the current optimal damping factor according to the fitting difference between the electromagnetic observation data and the electromagnetic theoretical data; A Hessian matrix calculation module, used for calculating the current Hessian matrix according to the apparent resistivity and impedance phase in the electromagnetic theory data; A model correction amount calculation module, used for calculating the correction amount of the current geoelectric model according to the current objective function, the current optimal damping factor and the current Hessian matrix; A model updating module, used for correcting the current geoelectric model according to the correction amount to obtain a current updated geoelectric model; A judgment module, used to judge whether the current updated geoelectric model meets the preset iteration condition; If not, the current updated geoelectric model is set as the current geoelectric model, and the step of "establishing the current objective function according to the standard error term between the electromagnetic observation data and the current electromagnetic theoretical data and the focusing constraint condition of the current geoelectric model" in the objective function construction module is returned to be executed; If satisfied, the current updated geoelectric model is output to obtain induced polarization parameters extracted based on the electromagnetic observation data.

Images