Weighted clustering constraint-based multi-constraint inversion method for transient electromagnetic data containing induced polarization effect

Abstract

The invention discloses a weighted clustering constraint-based multi-constraint inversion method for transient electromagnetic data containing an induced polarization effect. The method comprises the steps of constructing a smooth constraint term of an inversion parameter model; generating an inversion parameter clustering constraint term according to the prior physical property information, a preset clustering ambiguity value and the number of clustering centers; presetting an initial inversion parameter clustering center, and presetting an initial clustering weighted value of each inversion parameter as 1; forming a multi-constraint inversion objective function by the data fitting item, the parameter model constraint item and the parameter clustering constraint item; and updating a parameter clustering membership degree, a clustering center and a clustering weighted value of each inversion parameter in inversion iteration, and substituting into the multi-constraint objective function. And solving an extreme value of the multi-constraint objective function, and performing inversion iteration to update each parameter model until inversion meets an iteration termination condition. The method has the advantages of being clear in logic, stable, reliable and the like, and has high practical application significance in the technical field of geophysical exploration.

Description

Technical Field

[0001] The invention relates to the technical field of geophysical exploration, in particular to a multi-constraint inversion method of transient electromagnetic data containing induced polarization effect based on weighted clustering constraints. Background Art

[0002] After decades of development, transient electromagnetic method has become an effective and important method and technology in field exploration tasks. Now it has important application value in many fields, especially in metal and water resource exploration and environmental monitoring. With the rapid development of transient electromagnetic measuring instruments, the time of signal acquisition has gradually increased and the signal-to-noise ratio of late signals has gradually increased. It is found that when measuring in areas where polarized bodies exist, obvious rapid decay anomalies in the response curve and a large number of negative values ​​in the late data are often observed. After a large number of scholars' research, it is believed that this phenomenon is caused by the induced polarization effect of underground polarized bodies.

[0003] A large number of previous studies have shown that when transient electromagnetic data are affected by the induced polarization effect, the conventional single resistivity inversion method cannot effectively fit the measured data, and the resistivity distribution information reflected by the inversion results obtained is quite different from the actual underground distribution. Therefore, a large number of studies have been invested in this field in recent years for the subsequent simulation interpretation and information extraction of this phenomenon. In the existing technology, such as Li Lun's "Research on Transient Electromagnetic Forward and Inversion of Long Wire Source with Induced Polarization Effect", Liao Wenpeng's "Research on Transient Electromagnetic Forward and Inversion of Electrical Source with Induced Polarization Effect Based on Genetic Algorithm", Zhang Feiliang's "Research on Transient Electromagnetic Two-Dimensional Forward Modeling Considering Induced Polarization Effect", etc., they can well simulate the transient electromagnetic response data affected by the induced polarization effect by introducing the multi-parameter Cole-Cole model complex resistivity that can quantitatively describe the induced polarization effect of underground rocks and minerals, and thus realize the subsequent inversion extraction work. However, when performing inversion extraction, since multiple parameters need to be inverted simultaneously, the originally unstable inversion problem becomes more complicated, which greatly reduces the reliability of the inversion results.

[0004] Another example is the disclosed technology of "Patent Publication No. CN115793064A, entitled "An Improved Method for Extracting Polarized Information from Semi-Aerial Transient Electromagnetic Data", which includes: intercepting early data for conventional resistivity inversion to obtain an initial resistivity model, and then combining multiple parameters for inversion extraction, while introducing parameter range constraints and lateral smoothness constraints of the parameter model to improve the stability of the inversion.

[0005] The above technology has the following problems: the extraction method only uses parameter model smooth constraints to improve the stability of inversion, which is still insufficient to significantly improve the stability and reliability of multi-parameter inversion, and the weights of the smooth constraints used are consistent, which will cause the inversion results to be over-smoothed and lose boundary information.

[0006] Therefore, it is urgent to propose a multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints that is simple in logic, robust and reliable. Summary of the invention

[0007] In view of the above problems, the purpose of the present invention is to provide a multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints. The technical solution adopted by the present invention is as follows:

[0008] The multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints is characterized by comprising the following steps:

[0009] Step S1, according to the number and position information of the measuring points of the measuring line in the detection area, preset the smooth constraint factor of the current measuring line; the measuring line is provided with a plurality of measuring points; at the same time, according to the distance between the adjacent measuring points on the measuring line, preset the smooth constraint weight, generate a plurality of weighted parameter smooth constraint factors, and construct the smooth constraint item of the inversion parameter model;

[0010] Step S2, combining the prior physical property parameter information of the detection area, counting the number of cluster centers, the number of inversion parameter types involved in clustering, and the number of parameters using weighted control, and constructing inversion parameter clustering constraint items;

[0011] Step S3, obtaining the observation data collected by the survey line, and constructing the survey line data fitting term; combining the survey line data fitting term with the inversion parameter model smooth constraint term and the inversion parameter clustering constraint term to generate a multi-constraint inversion objective function;

[0012] Step S4, preset the initial cluster center value, and set the initial cluster weight value of any parameter to 1;

[0013] Step S5, updating the membership value of any inversion parameter, and obtaining the membership value to update the cluster center; using the distance between the membership value of the inversion parameter and the cluster center updated by the membership value to update the cluster weighted value of any inversion parameter;

[0014] Step S6, combining the cluster weighted values ​​of the inversion parameters and using a damped least squares inversion method to solve the extreme value of the multi-constrained inversion objective function, and updating the model corresponding to any inversion parameter;

[0015] Step S7, obtain the fitting difference between the response data corresponding to the updated parameter model and the observation data collected by the survey line, and preset a fitting difference threshold; if the fitting difference between the response data corresponding to the updated parameter model and the observation data collected by the survey line is less than the fitting difference threshold, output the multi-parameter model corresponding to the detection area; otherwise, return to step S5 for updating.

[0016] Compared with the prior art, the present invention has the following beneficial effects:

[0017] (1) The present invention ingeniously uses the number and location information of the measuring points in the survey line and designs a parameter model smoothing factor weighted by the distance factor, so as to achieve the purpose of applying strong smoothing constraints to adjacent measuring points that are close to each other and applying weak smoothing constraints to adjacent measuring points that are far away, so that the smoothing constraints of the inversion results are more reasonable.

[0018] (2) The present invention cleverly uses a clustering algorithm to integrate prior physical property parameter information into the inversion, maximizing the auxiliary effect of the prior information on the inversion results and greatly improving the stability of multi-parameter inversion.

[0019] (3) The present invention cleverly updates the membership value of any inversion parameter, obtains the membership value to update the cluster center, and uses the distance between the membership value of the inversion parameter and the membership value to update the cluster center to update the cluster weight value of any inversion parameter. Then, the damped least squares inversion method is used to solve the extreme value of the multi-constrained inversion objective function, and the parameter model corresponding to any inversion parameter is updated. The present invention utilizes the difference in clustering effects between different parameters and true physical property values, and adaptively applies clustering weights through corresponding parameters, thereby further improving the effectiveness of clustering constraints and preventing excessive clustering of various parameters.

[0020] (4) The present invention adopts a multi-constrained inversion method, which effectively introduces the parameter model smoothness constraint term and the parameter clustering constraint term into the objective function, greatly improving the stability of the multi-parameter inversion and the reliability of the inversion results.

[0021] In summary, the present invention has the advantages of clear logic, robustness and reliability, and can realize reliable multi-parameter information extraction of transient electromagnetic data containing induced polarization effect, and has high practical value and promotion value in the field of geophysical exploration technology. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments are briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope of protection. For those skilled in the art, other related drawings can be obtained based on these drawings without creative work.

[0023] Figure 1 It is a schematic diagram of the process of the present invention.

[0024] Figure 2 It is a schematic diagram of parameters and data of Model 1 in the invention.

[0025] Figure 3 This is the inversion result of model 1. DETAILED DESCRIPTION

[0026] In order to make the purpose, technical scheme and advantages of this application clearer, the present invention is further described below in conjunction with the accompanying drawings and embodiments, and the embodiments of the present invention include but are not limited to the following embodiments. Based on the embodiments in this application, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of this application.

[0027] In this embodiment, the term "and / or" is merely a term used to describe the association relationship of associated objects, indicating that three relationships may exist. For example, A and / or B may represent three situations: A exists alone, A and B exist at the same time, and B exists alone.

[0028] The terms "first" and "second" in the description and claims of this embodiment are used to distinguish different objects rather than to describe a specific order of objects. For example, a first target object and a second target object are used to distinguish different target objects rather than to describe a specific order of target objects.

[0029] In the embodiments of the present application, words such as "exemplary" or "for example" are used to indicate examples, illustrations or descriptions. Any embodiment or design described as "exemplary" or "for example" in the embodiments of the present application should not be interpreted as being more preferred or more advantageous than other embodiments or designs. Specifically, the use of words such as "exemplary" or "for example" is intended to present related concepts in a specific way.

[0030] In the description of the embodiments of the present application, unless otherwise specified, the meaning of "multiple" refers to two or more than two. For example, multiple processing units refer to two or more processing units; multiple systems refer to two or more systems.

[0031] like Figure 1 As shown, this embodiment provides a multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints, which includes the following steps:

[0032] The first step is to consider transient electromagnetic data with induced polarization effect, and then use the multi-parameter Cole-Cole resistivity model to quantitatively describe the induced polarization effect of rocks and minerals for forward and inverse calculations, which is expressed as:

[0033]

[0034] in, represents the dispersive resistivity including induced polarization effect; ρ0 represents the zero-frequency resistivity; m0 represents the charging rate; τ represents the time constant; and c represents the frequency correlation coefficient.

[0035] Here, according to the number of measuring points and position information of the survey line in the detection area, the smooth constraint factor of the current survey line is preset; a plurality of measuring points are arranged on the survey line; at the same time, the smooth constraint weight is preset according to the distance between adjacent measuring points on the survey line, a plurality of weighted parameter smooth constraint factors are generated, and the smooth constraint item of the inversion parameter model is constructed.

[0036] Here, the expression of the weighted parameter smooth constraint factor is:

[0037]

[0038] in, Indicates that it acts on the zero-frequency resistivity parameter and the distance between the measuring points is D ref The weighted value of Indicates that it acts on the charging rate parameter and the distance between the measuring points is D ref The weighted value of Indicates that it acts on the time constant parameter and the distance between the measuring points is D ref The weighted value of It means that it acts on the frequency correlation coefficient parameter and the distance between the measuring points is D ref dis represents the distance matrix between each measuring point obtained according to the position information of the measuring point; k scale A constant that controls how much distance affects the weight.

[0039] Then, the smooth constraint term of the inversion parameter model is constructed, and its expression is:

[0040]

[0041] Among them, M represents the inversion parameter set, which consists of zero-frequency resistivity ρ0, charging rate m0, time constant τ and frequency correlation coefficient c; R represents the weighted parameter smoothness constraint factor.

[0042] In the second step, the number of cluster centers, the number of inversion parameter types involved in clustering, and the number of parameters using weighted control are counted based on the prior physical parameter information of the detection area, and the inversion parameter clustering constraint term is constructed. The expression is:

[0043]

[0044] Where M represents the inversion parameter set; N represents the number of inversion parameter types involved in clustering; t k Represents the true prior physical property parameter information under the control of the weight factor of the kth cluster center; represents the number of cluster centers, that is, the number of true prior physical property parameter information; u ijrepresents the membership matrix of the i-th type inversion parameter for the j-th cluster center; q represents the constant controlling the clustering ambiguity; ω k represents the weight factor controlling the kth cluster center on the inversion parameters; M ik represents the i-th type inversion parameter controlled by the weight factor of the k-th cluster center; v jk represents the jth cluster center obtained by inversion iteration update under the control of the weight factor of the kth cluster center; η j Represents the constant that controls the weight of the jth cluster center; d is the number of parameters involved in weighting, d≤N.

[0045] The third step is to obtain the observation data collected by the survey line and construct the survey line data fitting term; the survey line data fitting term is combined with the inversion parameter model smooth constraint term and the inversion parameter clustering constraint term to generate a multi-constraint inversion objective function. The expression is:

[0046]

[0047] Among them, α represents the weight constant of the smooth term of the control parameter model, and its value is 1; β represents the clustering constraint term of the control parameter, and its value is 2; Represents the line data fitting term.

[0048] Then, the line data fitting term The expression is:

[0049]

[0050] Among them, W d represents the observation data weighting matrix; d obs represents the observed data; F(M) represents the forward operator.

[0051] The fourth step is to preset the initial cluster center value and set the initial cluster weight value of any parameter to 1, that is:

[0052] ω k =1, k=1,2…d;

[0053] The fifth step is to update the membership value of any inversion parameter and obtain the membership value to update the cluster center; the distance between the membership value of the inversion parameter and the cluster center updated by the membership value is used to update the cluster weight value of any inversion parameter.

[0054] Here, the membership value u of the inversion parameter ij The expression is:

[0055]

[0056] Among them, ω krepresents the weight factor controlling the kth cluster center on the inversion parameters; M ik represents the i-th type inversion parameter controlled by the weight factor of the k-th cluster center; v jk It represents the jth cluster center obtained by inversion iterative update under the control of the weight factor of the kth cluster center.

[0057] In addition, the membership value updates the cluster center v j The expression is:

[0058]

[0059] Among them, M i represents the i-th type inversion parameter; ω k represents the weight factor controlling the kth cluster center on the inversion parameters; t kj Represents the jth true prior physical property parameter information under the control of the weight factor of the kth cluster center.

[0060] Here, the distance between the membership value of the inversion parameter and the membership value to update the cluster center is used to update the cluster weight value of any inversion parameter. Its expression is:

[0061]

[0062] Among them, ω k Represents the weight factor that controls the k-th cluster center on the inversion parameters.

[0063] In the sixth step, the extreme value of the multi-constrained inversion objective function is solved by combining the cluster weighted values ​​of the inversion parameters and using the damped least squares inversion method, and the model corresponding to any inversion parameter is updated, which is expressed as:

[0064] ΔM=A -1 B

[0065]

[0066] Where J represents the Jacobian matrix; λ represents the damping factor; I represents the identity matrix; u ij Represents the membership matrix of the i-th type inversion parameter for the j-th cluster center.

[0067] In the seventh step, the fitting difference between the response data corresponding to the updated parameter model and the observation data collected by the survey line is obtained, and a fitting difference threshold is preset; if the fitting difference between the response data corresponding to the updated parameter model and the observation data collected by the survey line is less than the fitting difference threshold, the multi-parameter model corresponding to the detection area is output; otherwise, return to step 5 for updating.

[0068] Case 1

[0069] like Figure 2 As shown, this embodiment provides a multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints, which designs a typical Kimberlite area model. Figure 2 (a) is the model distribution of a typical Kimberlite area model. Figure 2 (b) in the figure is the parameter information of the model of a typical Kimberlite area. Here, the resistivity of the 20-meter-thick surface layer is 100Ωm, there is a crater-shaped polarization body in the middle (zero-frequency resistivity 10Ωm, charging rate 0.5, time constant 3ms, frequency correlation coefficient 0.5), and the bottom is a volcanic channel with a resistivity of 50Ωm, and the background rock resistivity is 500Ωm. A total of 51 transient electromagnetic measurement points were designed horizontally, with a measurement point interval of 20 meters. A rectangular transmitting coil (100m×100m) was used, and a transient electromagnetic device was set up for measurement at the center point. A one-dimensional forward simulation was performed on each measurement point and 3% random noise was applied to generate simulation data.

[0070] Case 2

[0071] like Figure 3 As shown, this case is a single-point one-dimensional inversion based on the Cole-Cole model for simulated data and an inversion trial calculation of the data using the inversion method in this embodiment. The initial model is set to 10 layers of uniform layered half-space, each layer is 10 meters thick, and the parameters of each layer are zero-frequency resistivity 100Ωm, charging rate 0.05, time constant 1ms, and frequency correlation coefficient 0.5. When performing the trial calculation of the inversion method of this application, the real physical parameter information of the surface layer and the polarization body is put into the cluster inversion as prior information. The inversion results are shown in Figure 3 As shown, Figure 3 (a) in the figure shows the result of single-point multi-parameter inversion. In the absence of constraints, the continuity of the inversion result interface is poor, and it is impossible to distinguish the distribution and boundary information of underground geology. At the same time, the induced polarization parameter information of the intermediate polarization body deviates greatly from the actual situation. Figure 3 (b) in the figure shows the weighted clustering multi-constraint inversion result. The distribution and interface of each geological body have been significantly improved. At the same time, the results of the induced polarization parameters can accurately reflect the location boundary and physical property distribution information of the underground polarization body, which greatly improves the stability and reliability of the inversion and provides reliable multi-parameter information for subsequent interpretation.

[0072] The above embodiments are only preferred embodiments of the present invention and are not intended to limit the protection scope of the present invention. Any changes that adopt the design principles of the present invention and are made through non-creative work on this basis should fall within the protection scope of the present invention.

Claims

1. A multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints, characterized in that: The following steps are involved: Step S1, according to the number and position information of the measuring points of the measuring line in the detection area, preset the smooth constraint factor of the current measuring line; the measuring line is provided with a plurality of measuring points; at the same time, according to the distance between the adjacent measuring points on the measuring line, preset the smooth constraint weight, generate a plurality of weighted parameter smooth constraint factors, and construct the smooth constraint item of the inversion parameter model; Step S2, combining the prior physical property parameter information of the detection area, counting the number of cluster centers, the number of inversion parameter types involved in clustering, and the number of parameters using weighted control, and constructing inversion parameter clustering constraint items; Step S3, obtaining the observation data collected by the survey line, and constructing the survey line data fitting term; combining the survey line data fitting term with the inversion parameter model smooth constraint term and the inversion parameter clustering constraint term to generate a multi-constraint inversion objective function; Step S4, preset the initial cluster center value, and set the initial cluster weight value of any parameter to 1; Step S5, updating the membership value of any inversion parameter, and obtaining the membership value to update the cluster center; The cluster weight value of any inversion parameter is updated by using the distance between the membership value of the inversion parameter and the cluster center updated by the membership value; Step S6, combining the cluster weighted values ​​of the inversion parameters and using a damped least squares inversion method to solve the extreme value of the multi-constrained inversion objective function, and updating the model corresponding to any inversion parameter; Step S7, obtain the fitting difference between the response data corresponding to the updated parameter model and the observation data collected by the survey line, and preset a fitting difference threshold; if the fitting difference between the response data corresponding to the updated parameter model and the observation data collected by the survey line is less than the fitting difference threshold, output the multi-parameter model corresponding to the detection area; otherwise, return to step S5 for updating.

2. The multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints according to claim 1 is characterized in that: The step S1 also includes: when the transient electromagnetic data is affected by the induced polarization effect, a multi-parameter Cole-Cole resistivity model that quantitatively describes the induced polarization effect of rocks and minerals is used to perform forward and inverse calculations, and the expression is: in, represents the dispersive resistivity including induced polarization effect; ρ0 represents the zero-frequency resistivity; m0 represents the charging rate; τ represents the time constant; and c represents the frequency correlation coefficient.

3. The multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints according to claim 1 or 2, characterized in that: In step S1, the smooth constraint factor of the current survey line is preset according to the number of survey points and position information of the survey line in the detection area; the smooth constraint weight is preset according to the distance between adjacent survey points on the survey line, and several weighted parameter smooth constraint factors are generated, and the expression thereof is: in, Indicates that it acts on the zero-frequency resistivity parameter and the distance between the measuring points is D ref The weighted value of Indicates that it acts on the charging rate parameter and the distance between the measuring points is D ref The weighted value of Indicates that it acts on the time constant parameter and the distance between the measuring points is D ref The weighted value of It means that it acts on the frequency correlation coefficient parameter and the distance between the measuring points is D ref dis represents the distance matrix between each measuring point obtained according to the position information of the measuring point; k scale A constant that controls how much distance affects the weight.

4. The multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints according to claim 2 is characterized in that: The inversion parameter model smoothness constraint term The expression is: Among them, M represents the inversion parameter set, which consists of zero-frequency resistivity ρ0, charging rate m0, time constant τ and frequency correlation coefficient c; R represents the weighted parameter smoothness constraint factor.

5. The multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints according to claim 4 is characterized in that: In step S2, the inversion parameter clustering constraint item The expression is: Where M represents the inversion parameter set; N represents the number of inversion parameter types involved in clustering; t k Represents the true prior physical property parameter information under the control of the weight factor of the kth cluster center; represents the number of cluster centers, that is, the number of true prior physical property parameter information; u ij represents the membership matrix of the i-th type inversion parameter for the j-th cluster center; q represents the constant controlling the clustering ambiguity; ω k represents the weight factor controlling the kth cluster center on the inversion parameters; M ik represents the i-th type inversion parameter controlled by the weight factor of the k-th cluster center; v jk represents the jth cluster center obtained by inversion iteration update under the control of the weight factor of the kth cluster center; η j Represents the constant that controls the weight of the jth cluster center; d is the number of parameters involved in weighting, d≤N.

6. The multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints according to claim 5 is characterized in that: In step S3, the line data fitting term The expression is: Among them, W d represents the observation data weighting matrix; d obs represents the observed data; F(M) represents the forward operator; The multi-constraint inversion objective function The expression is: Among them, α represents the weight constant of the smooth term of the control parameter model; β represents the clustering constraint term of the control parameter.

7. The multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints according to claim 5, characterized in that: In step S5, the membership value u of the inversion parameter ij The expression is: Among them, ω k represents the weight factor controlling the kth cluster center on the inversion parameters; M ik represents the i-th type inversion parameter controlled by the weight factor of the k-th cluster center; v jk It represents the jth cluster center obtained by inversion iterative update under the control of the weight factor of the kth cluster center.

8. The multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints according to claim 7 is characterized in that: The membership value updates the cluster center v j The expression is: Among them, M i represents the i-th type inversion parameter; ω k represents the weight factor controlling the kth cluster center on the inversion parameters; t kj Represents the jth true prior physical property parameter information under the control of the weight factor of the kth cluster center.

9. The multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints according to claim 8, characterized in that: The cluster weight value of any inversion parameter is updated by using the distance between the membership value of the inversion parameter and the membership value to update the cluster center. The expression is: Among them, ω k Represents the weight factor that controls the k-th cluster center on the inversion parameters.

10. The multi-constraint inversion method for transient electromagnetic data containing induced polarization effect based on weighted clustering constraints according to claim 9, characterized in that: In step S6, the extreme value of the multi-constrained inversion objective function is solved based on the damped least squares inversion method, and the model corresponding to any inversion parameter is updated, and its expression is: ΔM=A- 1 B Where J represents the Jacobian matrix; λ represents the damping factor; I represents the identity matrix; u ij Represents the membership matrix of the i-th type inversion parameter for the j-th cluster center.

Images