Multi-constraint inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraints

Abstract

The present invention discloses a multi-constraint inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraints, including: constructing a smooth constraint term for the inversion parameter model; presetting the clustering fuzziness value and the number of clustering centers according to prior physical property information to generate an inversion parameter clustering constraint term; presetting the initial clustering centers of the inversion parameters and presetting the initial clustering weights of each inversion parameter to 1; constructing a multi-constraint inversion objective function by combining the data fitting term, the parameter model constraint term and the parameter clustering constraint term; updating the parameter clustering membership, the clustering centers and the clustering weights of each inversion parameter in the inversion iteration and substituting them into the multi-constraint objective function. Taking the extreme value of the multi-constraint objective function, inversely iteratively updating each parameter model until the inversion meets the iteration termination condition. The present invention has the advantages of clear logic, robustness and reliability, and has strong practical application significance in the field of geophysical exploration technology.

Description

Technical Field

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

[0002] After decades of development, the transient electromagnetic method has become an effective and important method technology in field exploration tasks. Nowadays, it has important application values in many fields, especially in metal and water resource exploration and environmental monitoring. With the rapid development of transient electromagnetic measurement instruments, the acquisition time of signals gradually increases and the signal-to-noise ratio of late-stage signals gradually improves. It is found that when measuring in areas with polarizable bodies, obvious rapid attenuation anomalies often appear in the response curves and a large number of negative values appear in the late-stage data. After extensive research by many scholars, it is considered that this phenomenon is caused by the induced polarization effect of underground polarizable bodies.

[0003] A large number of previous studies have shown that when transient electromagnetic data is affected by the induced polarization effect, the conventional single resistivity inversion method cannot effectively fit the measured data, and the resistivity distribution information obtained from the inversion results has a large deviation from the actual underground distribution. Therefore, for the subsequent simulation explanation and information extraction of this phenomenon, a large number of studies have been invested in this field in recent years. In the prior art, such as the technologies of Li Lun's "Research on the Forward and Inverse Modeling of Long Wire Source Transient Electromagnetics with Induced Polarization Effect", Liao Wenpeng's "Research on the Forward and Inverse Modeling of Electric Source Transient Electromagnetics with Induced Polarization Effect Based on Genetic Algorithm", Zhang Feiliang's "Research on Two-Dimensional Forward Modeling of Transient Electromagnetics Considering Induced Polarization Effect", etc., through the multi-parameter Cole-Cole model complex resistivity that can quantitatively describe the induced polarization effect of underground rocks and ores, the transient electromagnetic response data affected by the induced polarization effect can be well simulated, and then the subsequent inversion extraction work can be realized. However, when performing inversion extraction, since multiple parameters need to be inverted simultaneously, the originally unstable inversion problem becomes more complex, greatly reducing the reliability of the inversion results.

[0004] Another example is the publicly disclosed technology of "Patent Publication No. CN115793064A, titled An Improved Method for Extracting Induced Polarization Information in Semi-Airborne Transient Electromagnetic Data", which includes: intercepting early-stage data for conventional resistivity inversion to obtain an initial resistivity model, and then jointly inverting and extracting multiple parameters, and introducing parameter range constraints and lateral smooth constraints of the parameter model to improve the stability of the inversion.

[0005] The above technologies have the following problems: Only the smooth constraint of the parameter model is used in the extraction method to improve the stability of the inversion, which is still insufficient for significantly improving the stability and reliability of multi-parameter inversion, and the weights of the adopted smooth constraints are the same, which will cause the inversion results to be overly smoothed and lose boundary information.

[0006] Therefore, there is an urgent need to propose a multi-constrained inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraints, which has simple logic, is robust and reliable. Summary of the Invention

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

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

[0009] Step S1, according to the number of measuring points and position information of the measuring lines in the detection area, preset the smooth constraint factor of the current measuring line; several measuring points are set on the measuring line; at the same time, preset the smooth constraint weight according to the distance between adjacent measuring points on the measuring line, generate several weighted parameter smooth constraint factors, and construct the smooth constraint term of the inversion parameter model.

[0010] Step S2, combine the prior physical property parameter information of the detection area, count the number of clustering centers, the number of types of inversion parameters participating in clustering, and the number of parameters controlled by weighting, and construct the clustering constraint term of the inversion parameters.

[0011] Step S3, obtain the observed data collected by the measuring line, and construct the fitting term of the measuring line data; combine the fitting term of the measuring line data with the smooth constraint term of the inversion parameter model and the clustering constraint term of the inversion parameters to generate a multi-constrained inversion objective function.

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

[0013] Step S5, update and obtain the membership value of any inversion parameter, and obtain the clustering center updated by the membership value; use the distance between the membership value of the inversion parameter and the clustering center updated by the membership value to update the clustering weight value of any inversion parameter.

[0014] Step S6, combine the clustering weight value of the inversion parameter and use the damped least squares inversion method to solve the extreme value of the multi-constrained inversion objective function, and update the corresponding model of any inversion parameter.

[0015] Step S7, obtain the fitting difference between the response data corresponding to the updated parameter model and the observed data collected by the measuring line, and preset the fitting difference threshold; if the fitting difference between the response data corresponding to the updated parameter model and the observed data collected by the measuring 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 adopts the number and position information of measurement points in the measuring line, and designs a parametric model smoothing factor weighted based on distance factors, so as to impose a strong smoothing constraint on adjacent measurement points that are relatively close, and impose a weak smoothing constraint on adjacent measurement points that are far apart, making the smoothing constraint of the inversion result more reasonable.

[0018] (2) The present invention ingeniously incorporates the prior physical property parameter information into the inversion by using the clustering algorithm, maximizing the auxiliary effect of the prior information on the inversion result, and greatly improving the stability of multi-parameter inversion.

[0019] (3) The present invention ingeniously updates the membership value of any inversion parameter, and obtains the updated clustering center of the membership value. The clustering weight value of any inversion parameter is updated by using the distance between the membership value of the inversion parameter and the updated clustering center of the membership value. Then, the extreme value of the multi-constrained inversion objective function is solved by using the damped least squares inversion method, and the parameter model corresponding to any inversion parameter is updated. By utilizing the difference in the clustering effect between different parameters and the true physical property values, the present invention adaptively applies clustering weights to the corresponding parameters, further improving the effectiveness of the clustering constraint and preventing the phenomenon of over-clustering of each parameter.

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

[0021] In summary, the present invention has the advantages of clear logic, robustness and reliability, etc., and realizes reliable multi-parameter information extraction for transient electromagnetic data containing induced polarization effects, having high practical value and popularization 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 will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention, and therefore should not be regarded as limiting the protection scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative efforts.

[0023] Figure 1 It is a schematic flow chart of the present invention.

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

[0025] Figure 3 It is the inversion result of Model 1. DETAILED DESCRIPTION OF THE INVENTION

[0026] To make the objectives, technical solutions, and advantages of this application clearer, the following further describes the present invention in conjunction with the accompanying drawings and embodiments. The implementation manners 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 those of ordinary skill in the art without creative efforts shall fall within the scope of protection of this application.

[0027] In this embodiment, the term "and / or" merely describes the association relationship of associated objects and indicates that three relationships may exist. For example, A and / or B may represent three situations: A exists alone, A and B exist simultaneously, and B exists alone.

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

[0029] In the embodiments of this application, words such as "exemplary" or "for example" are used to indicate examples, illustrations, or explanations. Any embodiment or design solution described as "exemplary" or "for example" in the embodiments of this application should not be construed as being more preferred or having more advantages than other embodiments or design solutions. Rather, the use of words such as "exemplary" or "for example" is intended to present relevant concepts in a specific manner.

[0030] In the description of the embodiments of this application, unless otherwise specified, the meaning of "a plurality of" refers to two or more. For example, a plurality of processing units refers to two or more processing units; a plurality of systems refers to two or more systems.

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

[0032] In the first step, when considering transient electromagnetic data with induced polarization effect, it is necessary to use the multi-parameter Cole-Cole resistivity model that quantitatively describes the induced polarization effect of rock and ore for forward and inverse calculations, and its expression is:

[0033]

[0034] where represents the dispersive resistivity with induced polarization effect; ρ0 represents the zero-frequency resistivity; m0 represents the chargeability; τ represents the time constant; c represents the frequency correlation coefficient.

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

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

[0037]

[0038] Among them, represents the weighted value acting on the zero-frequency resistivity parameter when the measuring point spacing is D ref ; represents the weighted value acting on the chargeability parameter when the measuring point spacing is D ref ; represents the weighted value acting on the time constant parameter when the measuring point spacing is D ref ; represents the weighted value acting on the frequency correlation coefficient parameter when the measuring point spacing is D ref ; dis represents the distance matrix between each measuring point obtained according to the position information of the measuring points; k scale represents a constant that controls the influence degree of the distance on the weight.

[0039] Subsequently, 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 is composed of zero-frequency resistivity ρ0, chargeability m0, time constant τ, and frequency correlation coefficient c; R represents the weighted parameter smooth constraint factor.

[0042] In the second step, combining the prior physical property parameter information of the detection area, the number of clustering centers, the number of types of inversion parameters participating in clustering, and the number of parameters controlled by weighted are statistically analyzed, and the clustering constraint term of the inversion parameters is constructed. The expression of this clustering constraint term of the inversion parameters is:

[0043]

[0044] Among them, M represents the inversion parameter set; N represents the number of types of inversion parameters participating in clustering; t k represents the true prior physical property parameter information under the control of the weight factor at the kth clustering center; represents the number of clustering centers, that is, the number of true prior physical property parameter information; u ijrepresents the membership degree matrix of the i-th type of inversion parameter with respect to the j-th clustering center; q represents a constant controlling the clustering fuzziness; ω k represents the weight factor of the k-th clustering center for the inversion parameter; M ik represents the i-th type of inversion parameter under the control of the weight factor of the k-th clustering center; v jk represents the j-th clustering center obtained by inversion iteration update under the control of the weight factor of the k-th clustering center; η j represents a constant controlling the weight of the j-th clustering center; d is the number of parameters participating in weighting, and d ≤ N.

[0045] Thirdly, obtain the observed data collected by the survey line and construct the survey line data fitting term; combine 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. This multi-constraint inversion objective function is expressed as:

[0046]

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

[0048] Then, the survey line data fitting term is expressed as:

[0049]

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

[0051] Fourthly, preset the initial clustering center value and set the initial clustering weight value of any parameter to 1, that is:

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

[0053] Fifthly, update and obtain the membership degree value of any inversion parameter, and obtain the clustering center updated by the membership degree value; use the distance between the membership degree value of the inversion parameter and the clustering center updated by the membership degree value to update the clustering weight value of any inversion parameter.

[0054] Here, the membership degree value u of the inversion parameter ij is expressed as:

[0055]

[0056] where ω kDenote the weight factor of the k-th clustering center for the inversion parameters; M ik Denote the i-th type of inversion parameters under the control of the weight factor of the k-th clustering center; v jk Denote the j-th clustering center obtained by inversion iteration update under the control of the weight factor of the k-th clustering center.

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

[0058]

[0059] Among them, M i Denote the i-th type of inversion parameters; ω k Denote the weight factor of the k-th clustering center for the inversion parameters; t kj Denote the j-th true prior physical property parameter information under the control of the weight factor of the k-th clustering center.

[0060] Here, the clustering weighted value of any inversion parameter is updated by using the distance between the membership degree value of the inversion parameter and the clustering center updated by the membership degree value, and its expression is

[0061]

[0062] Among them, ω k Denote the weight factor of the k-th clustering center for the inversion parameters.

[0063] Step 6: Combine the clustering weighted value of the inversion parameters and use the damped least squares inversion method to solve the extreme value of the multi-constrained inversion objective function, and update the model corresponding to any inversion parameter. Its expression is:

[0064] ΔM = A -1 B

[0065]

[0066] Among them, J represents the Jacobian matrix; λ represents the damping factor; I represents the identity matrix; u ij Denote the membership degree matrix of the i-th type of inversion parameters for the j-th clustering center.

[0067] Step 7: Obtain the fitting difference between the response data corresponding to the updated parameter model and the observed data collected by the survey line, and preset the fitting difference threshold; if the fitting difference between the response data corresponding to the updated parameter model and the observed 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 5 for update.

[0068] Case 1

[0069] Such asFigure 2 As shown in the figure, this embodiment provides a multi-constraint inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraints, and designs a typical kimberlite area model. Among them, Figure 2 (a) in the figure shows the model distribution of the typical kimberlite area model, Figure 2 (b) in the figure shows the parameter information of the typical kimberlite area model. Here, the resistivity of the surface layer with a thickness of 20 meters is 100 Ωm. There is a volcanic crater-shaped polarizable body in the middle (zero-frequency resistivity 10 Ωm, chargeability 0.5, time constant 3 ms, frequency correlation coefficient 0.5). The bottom is a volcanic channel with a resistivity of 50 Ωm, and the resistivity of the background rock mass is 500 Ωm. A total of 51 transient electromagnetic measurement points are designed horizontally, and the interval between the measurement points is 20 meters. A rectangular transmitting coil (100 m × 100 m) is used, and the transient electromagnetic device is set to measure at the center point. One-dimensional forward simulation is carried out for each measurement point and 3% random noise is applied to generate simulation data.

[0070] Case 2

[0071] As Figure 3 shown in the figure, in this case, the single-point one-dimensional inversion based on the Cole-Cole model and the inversion method in this embodiment are used to perform inversion calculations on the simulation data. The initial models are all set as 10-layer homogeneous layered half-spaces, and the thickness of each layer is 10 meters. The parameters of each layer are zero-frequency resistivity 100 Ωm, chargeability 0.05, time constant 1 ms, and frequency correlation coefficient 0.5. Among them, when performing the inversion method calculation of this application, the true physical property parameter information of the surface layer and the polarizable body is used as prior information and put into the clustering inversion. The inversion results are as Figure 3 shown in the figure, where Figure 3 (a) in the figure represents the result of single-point multi-parameter inversion. In the absence of constraint conditions, the continuity of the inversion result interface is poor, and it is impossible to distinguish the distribution and boundary information of underground geology well. At the same time, there is a large deviation between the induced polarization parameter information of the middle polarizable body and the actual situation. Figure 3 (b) in the figure represents the result of weighted clustering multi-constraint inversion. 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 position boundary and physical property distribution information of the underground polarizable body, greatly improving the stability and reliability of the inversion, and providing reliable multi-parameter information for subsequent interpretation.

[0072] The above embodiments are only the preferred embodiments of the present invention, and do not limit the protection scope of the present invention. Any changes made by using the design principle of the present invention and non-creative labor on this basis shall fall within the protection scope of the present invention.

Claims

1. A multi-constraint inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraints, characterized in that It includes the following steps: Step S1: Preset the smooth constraint factor of the current survey line according to the number of measurement points and position information of the measurement points on the survey line in the detection area; several measurement points are set on the survey line; meanwhile, preset the smooth constraint weight according to the distance between adjacent measurement points on the survey line, generate several weighted parameter smooth constraint factors, and construct the smooth constraint term of the inversion parameter model; Step S2: Combine the prior physical property parameter information of the detection area, count the number of clustering centers, the number of types of inversion parameters participating in clustering, and the number of parameters using weighted control, and construct the clustering constraint term of the inversion parameters; Step S3: Obtain the observed data collected by the survey line and construct the survey line data fitting term; combine the survey line data fitting term with the smooth constraint term of the inversion parameter model and the clustering constraint term of the inversion parameters to generate a multi-constraint inversion objective function; Step S4: Preset the initial clustering center value and set the initial clustering weight value of any parameter to 1; Step S5: Update to obtain the membership value of any inversion parameter and obtain the clustering center updated by the membership value; Update the clustering weight value of any inversion parameter by using the distance between the membership value of the inversion parameter and the clustering center updated by the membership value; Step S6: Combine the clustering weight value of the inversion parameter and use the damped least squares inversion method to solve the extreme value of the multi-constraint inversion objective function, and update 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 observed data collected by the survey line, and preset the fitting difference threshold; if the fitting difference between the response data corresponding to the updated parameter model and the observed 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 update.

2. The multi-constraint inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraints according to claim 1, characterized in that In the above Step S1, it also includes: when the transient electromagnetic data is affected by the induced polarization effect, the multi-parameter Cole-Cole resistivity model for quantitatively describing the induced polarization effect of rock and ore is used for forward and inverse calculations, and its expression is: Among them, represents the dispersion resistivity with induced polarization effect; ρ0 represents the zero-frequency resistivity; m0 represents the chargeability; τ represents the time constant; c represents the frequency-dependent coefficient.

3. The multi-constraint inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraints according to claim 1 or 2, characterized in that, In the above Step S1, preset the smooth constraint factor of the current survey line according to the number of measurement points and position information of the measurement points on the survey line in the detection area; preset the smooth constraint weight according to the distance between adjacent measurement points on the survey line, and generate several weighted parameter smooth constraint factors, and its expression is: Among them, represents the weighted value acting on the zero-frequency resistivity parameter when the measuring point spacing is D ref ; represents the weighted value acting on the chargeability parameter when the measuring point spacing is D ref ; represents the weighted value acting on the time constant parameter when the measuring point spacing is D ref ; represents the weighted value acting on the frequency correlation coefficient parameter when the measuring point spacing is D ref ; dis represents the distance matrix between each measuring point obtained according to the position information of the measuring points; k scale represents a constant that controls the influence degree of the distance on the weight.

4. The multi-constraint inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraint according to claim 2, characterized in that, The smooth constraint term of the inversion parameter model has the following expression: Where, M represents the inversion parameter set, which consists of zero-frequency resistivity ρ0, chargeability m0, time constant τ, and frequency correlation coefficient c; R represents the weighted parameter smooth constraint factor.

5. The multi-constraint inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraint according to claim 4, wherein In the step S2, the inversion parameter clustering constraint term is expressed as: Among them, M represents the set of inversion parameters; N represents the number of types of inversion parameters participating in clustering; t k represents the true prior physical property parameter information under the control of the weight factor of the k-th 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 of inversion parameter with respect to the j-th cluster center; q represents a constant controlling the clustering fuzziness; ω k represents the weight factor of the k-th cluster center for the inversion parameters; M ik represents the i-th type of inversion parameter under the control of the weight factor of the k-th cluster center; v jk represents the j-th cluster center obtained by inversion iteration update under the control of the weight factor of the k-th cluster center; η j represents a constant controlling the weight of the j-th cluster center; d is the number of parameters participating in weighting, d ≤ N.

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

7. The multi-constraint inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraints according to claim 5, wherein, In the step S5, the membership value u of the inversion parameter ij is expressed as: where, ω k represents the weight factor of the k-th clustering center for the inversion parameters; M ik represents the i-th type of inversion parameter under the control of the weight factor of the k-th clustering center; v jk represents the j-th clustering center obtained by inversion iteration update under the control of the weight factor of the k-th clustering center.

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

9. The multi-constraint inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraints according to claim 8, characterized in that Update the clustering weight value of any inversion parameter by using the distance between the membership value of the inversion parameter and the clustering center updated by the membership value, and its expression is where ω k represents the weight factor for the k-th cluster center to the inversion parameter.

10. The multi-constraint inversion method for transient electromagnetic data with induced polarization effect based on weighted clustering constraint according to claim 9, characterized in that, In the above Step S6, use the damped least squares inversion method to solve the extreme value of the multi-constraint inversion objective function, and update the model corresponding to any inversion parameter, 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 inversion parameter with respect to the j-th clustering center.

Images