High-boundary-sensitivity 3-d inversion method based on electromagnetic gradient constraint matrix

Abstract

The application discloses a kind of high boundary sensitivity three-dimensional inversion methods based on electromagnetic gradient constraint matrix, belong to electromagnetic signal inversion technical field, including: constructing initial inversion model and objective function, the model roughness term based on electromagnetic gradient constraint is introduced in the objective function;Electromagnetic gradient anomaly response under each measuring frequency at each measuring point is calculated, the electromagnetic gradient anomaly response vector of each measuring point is obtained, normalization processing is carried out, then combination is formed into matrix, obtains overall normalized observation gradient anomaly matrix, constructs electromagnetic gradient constraint matrix in data space;Electromagnetic gradient constraint matrix in data space is converted to model space by frequency-depth weighted mapping matrix;The objective function is solved using iterative optimization algorithm until model converges.This application realizes the dynamic adjustment of the variation amplitude weight of each grid element in model roughness term by introducing electromagnetic field gradient constraint matrix, and improves the target boundary identification precision.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of electromagnetic signal inversion, and relates to a high-resolution inversion method for three-dimensional electromagnetic data, in particular to a high-boundary-sensitivity three-dimensional inversion method based on an electromagnetic gradient constraint matrix. BACKGROUND

[0002] The underground resistivity structure imaging methods based on geophysical electromagnetic data mainly include three types: forward trial modeling method, apparent resistivity imaging method and inversion method. The three methods have their own advantages and are suitable for different geological scenes and resolution requirements. The forward trial modeling method generates a predicted response by constructing a hypothetical underground model, and continuously adjusts the model parameters to approach the observed data; when there is a lack of effective inversion program, this method has certain practical value in the interpretation of complex underground structures; but due to its subjective process and low efficiency, it has obvious limitations in practical application. The apparent resistivity imaging method realizes the preliminary characterization of the underground structure by directly mapping the observed electromagnetic response to the apparent resistivity; this method is fast in calculation and simple in operation, and is suitable for large-scale rapid evaluation, but it relies on the assumption of uniform half-space, and the depth resolution is limited, making it difficult to accurately reflect complex high-dimensional geological structures. In comparison, the electromagnetic inversion method realizes the quantitative and systematic reconstruction of the underground resistivity structure by minimizing the residual between the predicted data and the observed data; this method has high resolution and good generalization, and is suitable for fine interpretation of complex geological structures. In recent years, with the wide application of electromagnetic inversion methods in resource exploration, environmental investigation and geotechnical engineering, etc., its key role in the analysis of underground electrical structure has become increasingly prominent, and it has become one of the indispensable core technologies in geophysical exploration. At the same time, the continuous innovation of theoretical methods and the growing demand for high-precision and high-resolution imaging are driving the electromagnetic inversion technology to develop rapidly to a higher level.

[0003] Although geophysical electromagnetic inversion technology has made great progress, achieving accurate boundary determination and fine interpretation of geological information is still a continuous challenge. Recent advances in spatial adaptive regularization techniques can recover blocky and smooth features from complex geological models; including combining the minimum support gradient functional to enhance boundary clarity, applying spatially varying Lp norms according to local structural characteristics, adjusting roughness weights to strengthen structural similarity with the guide model, and using uncertainty information for local adjustment of regularization strength. The above strategies can improve the sensitivity to structural discontinuities while maintaining the structural stability of homogeneous regions, thereby improving the geological consistency of the inversion results.

[0004] However, these methods all rely to some extent on prior geological knowledge or initial models, which are difficult to obtain in practical applications. To address this limitation, another strategy is to design regularization schemes based on the characteristics of the observation data. Electromagnetic field gradient data is more sensitive to geological discontinuities and has better boundary identification capabilities; these characteristics can be used to improve the boundary identification capabilities of data interpretation methods. Summary of the Invention

[0005] To address the aforementioned technical problems, this invention provides a high-boundary-sensitivity 3D inversion method based on an electromagnetic gradient constraint matrix. This method directly incorporates electromagnetic gradient information into the inversion process, organically integrating its boundary enhancement characteristics into a robust inversion constraint framework. While ensuring model stability, it significantly improves the resolution of geological structure boundaries. This method effectively combines the advantages of gradient-based imaging for boundary identification with the quantitative rigor of inversion techniques. It provides a solid theoretical foundation for the quantitative interpretation of geophysical data and constructs a practically feasible solution. Compared to traditional apparent resistivity imaging methods and conventional inversion techniques relying on smoothing constraints, this invention, by introducing an electromagnetic field gradient constraint matrix, dynamically adjusts the weights of the variation amplitudes of each grid cell in the model's roughness term. This enhances boundary sensitivity while maintaining the continuity and stability of uniform regions. This mechanism effectively improves the sensitivity of the inversion objective equation to electrical boundary structures, significantly improves boundary identification accuracy, and effectively enhances the balance between data fitting consistency and model geological rationality.

[0006] According to one aspect of the present invention, a high boundary sensitivity three-dimensional inversion method based on an electromagnetic gradient constraint matrix is ​​provided, comprising: constructing an initial inversion model and an objective function, wherein the initial inversion model is an initial resistivity model set based on experience, serving as the starting point for inversion iteration, and the objective function consists of three parts: a data residual term, a reference model constraint term, and a model roughness term based on electromagnetic gradient constraints; calculating the electromagnetic gradient anomaly response of a geodetic model containing anomalies at each measurement frequency at each measurement point; combining the anomaly responses at all frequencies at each measurement point to obtain the electromagnetic gradient anomaly response vector for each measurement point; and normalizing the electromagnetic gradient anomaly response vector for each measurement point to obtain the normalized observed gradient anomaly vector for each measurement point. The process involves: combining the normalized observation gradient anomaly vectors of all measurement points into a matrix to obtain the overall normalized observation gradient anomaly matrix; constructing an electromagnetic gradient constraint matrix in the data space based on this matrix; transforming the electromagnetic gradient constraint matrix in the data space to the model space using a frequency-depth weighted mapping matrix, where the frequency-depth weighted mapping matrix reflects the sensitivity of electromagnetic gradients at different frequencies to different depths; and solving the objective function using an iterative optimization algorithm based on the transformed electromagnetic gradient constraint matrix in the model space. In each iteration, the algorithm jointly optimizes and minimizes the data fitting error and the structure term, updating the model iteration step term from the initial inversion model until the model converges.

[0007] Optionally, the objective function is defined as:

[0008] ;

[0009] in, Let be the objective function. For inversion model, For data residuals, For model structure items, Parameters for controlling the balance of two items; ,in, For reference model constraint terms, This is the roughness term of the model based on electromagnetic gradient constraints;

[0010] The data residuals Defined as:

[0011] ;

[0012] in, This is a diagonal weighted matrix that includes the uncertainty of the measured electromagnetic data. To measure electromagnetic data, Forward electromagnetic response;

[0013] The reference model constraint terms Defined as:

[0014] ;

[0015] in, To control the weighting coefficients for the reference model, As the reference model weight matrix, For reference model;

[0016] The roughness term of the model based on electromagnetic gradient constraints Defined as:

[0017] ;

[0018] ;

[0019] in, These are the gradient roughness weight control coefficients; This is the electromagnetic gradient constraint regularization term; The electromagnetic gradient constraint matrix; This is the model gradient operator.

[0020] Optionally, the calculation of the electromagnetic gradient anomaly response of the geodetic model containing the anomaly at each measurement frequency at each measurement point includes:

[0021] For each measurement frequency at each measuring point, the electromagnetic gradient anomaly response at that frequency is represented by the electromagnetic field spatial derivative. This electromagnetic field spatial derivative is calculated using the Euclidean distance between adjacent measuring points and the difference in electromagnetic field amplitude, as shown in the following formula:

[0022] ;

[0023] in, For measuring points Frequency measurement Abnormal response to electromagnetic gradient; For measuring points Frequency measurement The amplitude of the secondary electromagnetic field under the condition; For measuring points +1 measurement frequency The amplitude of the secondary electromagnetic field under the condition; and They are measuring points +1 and measuring point The x-coordinate in a three-dimensional coordinate system; For measuring points The set of nearest measuring points; Let j be the neighboring measurement point number; For measuring points and The Euclidean distance between them; For measuring points and The absolute difference in electromagnetic field amplitude between them.

[0024] Optionally, the anomalous responses at all frequencies of each measuring point are combined to obtain the electromagnetic gradient anomalous response vector for each measuring point, expressed as:

[0025] ;

[0026] in, For measuring points The electromagnetic gradient anomaly response vector; For measuring points Anomaly response of electromagnetic gradient at the first measurement frequency; For measuring points Anomaly response of electromagnetic gradient at the second measurement frequency; For measuring points First Anomaly response of electromagnetic gradient at a measurement frequency; It is the set of real numbers; This represents the total number of measurement points; For measuring points The total number of corresponding electromagnetic gradient anomaly response data, i.e. the total number of measured frequencies.

[0027] Optionally, the electromagnetic gradient anomaly response vector at each measuring point is normalized to obtain the normalized observed gradient anomaly vector at each measuring point, including:

[0028] Find the maximum value of the electromagnetic gradient anomaly response at all measurement points and frequencies, and use a preset minimum value;

[0029] The larger of the maximum value and the preset minimum value is selected to normalize the electromagnetic gradient anomaly response vector of each measuring point:

[0030] ;

[0031] in, For measuring points The normalized observed gradient anomaly vector; This represents the maximum value of the electromagnetic gradient anomaly response at all measurement points and frequencies; This is set to a minimum value to avoid division by zero errors.

[0032] Optionally, the normalized observation gradient anomaly vectors of all measurement points are combined into a matrix to obtain the overall normalized observation gradient anomaly matrix. Based on this overall normalized observation gradient anomaly matrix, an electromagnetic gradient constraint matrix in the data space is constructed, including:

[0033] By combining the normalized observation gradient anomaly vectors of all measurement points, the overall normalized observation gradient anomaly matrix is ​​obtained:

[0034] ;

[0035] in, This is the overall normalized observation gradient anomaly matrix; This is the normalized observation gradient anomaly vector for measurement point 1; This is the normalized observation gradient anomaly vector for measurement point 2; For measuring points The normalized observed gradient anomaly vector;

[0036] The overall normalized observation gradient anomaly matrix is ​​then used. As an electromagnetic gradient constraint matrix in the data space:

[0037] ;

[0038] in, This represents the electromagnetic gradient constraint matrix in the data space.

[0039] Optionally, transforming the electromagnetic gradient constraint matrix in the data space to the model space using a frequency-depth weighted mapping matrix includes:

[0040] Define the frequency-depth weighted mapping matrix corresponding to each measurement point as follows:

[0041] ;

[0042] in, For measuring points Frequency-depth weighted mapping matrix at the location; For measuring points The weighting function for the first frequency at the first depth; For measuring points The weighting function for the first frequency at the second depth; For measuring points In the The weighting function of the first frequency at the layer depth; For measuring points The weighting function for the second frequency at the first depth; For measuring points The weighting function for the second frequency at the second depth; For measuring points In the The weighting function for the second frequency at layer depth; For measuring points At the depth of the first layer A weighted function of frequencies; For measuring points At the second layer depth A weighted function of frequencies; For measuring points In the Layer depth A weighted function of frequencies; Total number of floors;

[0043] Each element in the frequency-depth weighted mapping matrix is ​​defined as follows:

[0044] ;

[0045] In the formula, For measuring points exist Frequency Weighting function for layer depth; This is a constant used to control the severity of the weight decay with depth mismatch; For measuring points Model No. Layer center Coordinate values; Permeability in free space; Electrical conductivity; Number the floors; For measuring points Model No. Layer center Coordinate values;

[0046] The electromagnetic gradient constraints in the data space are transformed into the model space using the following formula:

[0047] ;

[0048] ;

[0049] in, The electromagnetic gradient constraint matrix is ​​normalized in the model space to achieve boundary-aware adjustment; It is the identity matrix; It is a frequency-depth weighted mapping matrix; This is the frequency-depth weighted mapping matrix at measurement point 1; This is the frequency-depth weighted mapping matrix at measurement point 2; For measuring points The frequency-depth weighted mapping matrix at that location.

[0050] Optionally, solving the objective function using an iterative optimization algorithm includes:

[0051] No. The step value for the next iteration is obtained by solving the following equation:

[0052] ;

[0053] ;

[0054] ;

[0055] In the formula, For the first The left-hand side of the next iteration inversion equation, i.e., the coefficient matrix; For the first The model obtained from the first step of calculation and the first step The difference between the models obtained from the first step; For the first The right-hand side of the inversion equation in the next iteration, i.e., the gradient direction; For orthogonal Jacobian matrix; For the first Model regularization parameters in step calculation; For the first The model predicts data in the next iteration; For the first The model obtained through step-by-step calculation.

[0056] The beneficial effects of this invention are:

[0057] Traditional inversion methods often use model constraints (such as smoothness constraints) that tend to generate smooth, blurred-boundary models, which can obscure the true sharp boundaries of geological bodies. This invention introduces a model roughness term based on electromagnetic gradient constraints. This electromagnetic gradient constraint matrix originates from the spatial gradient of actual measurement data. By integrating electromagnetic gradient data into the inversion process, regions with drastic field changes (high gradient) in the data space typically correspond to the physical property interfaces (i.e., boundaries) of the subsurface medium. After transforming these regions to the model space through frequency-depth mapping, this matrix can adaptively assign lower constraint weights to high-gradient regions (boundaries) and higher constraint weights to low-gradient regions (homogeneous medium) during the inversion process. This method effectively relaxes the smoothness constraints at the boundaries, allowing the inversion model to produce more drastic and sharper changes at these locations. Therefore, the final inversion results can more clearly and accurately characterize the geometry and spatial distribution of anomalies (such as ore bodies and tectonic interfaces), greatly improving the ability to identify boundaries and spatial resolution.

[0058] The constraint information in this invention comes directly from the spatial variation characteristics of the observation data itself, rather than from artificially imposed mathematical assumptions. This reduces subjective bias and improves work efficiency and standardization. The data gradient is insensitive to the overall magnitude but sensitive to relative changes, which makes it somewhat effective in suppressing uniform noise. In addition, normalization (dividing by the maximum value) and preset minimum values ​​avoid numerical instability. While pursuing high resolution, it maintains good noise resistance, making the inversion process more stable and less prone to producing false anomalies unrelated to the real geological structure. The results are more reliable and trustworthy.

[0059] This invention introduces a frequency-depth weighted mapping matrix, which is based on the skin effect principle of electromagnetic field theory (the higher the frequency, the shallower the detection depth; the lower the frequency, the deeper the detection depth). By establishing a physical relationship between electromagnetic responses of different frequencies and model parameters of different depths, the abnormal responses observed in the data space are intelligently "assigned" to the model depth from which they most likely originate. This mapping makes the final gradient constraint matrix contain not only lateral structural information but also depth information, making the constraints more physically reasonable and improving the vertical resolution and depth positioning accuracy of the inversion results.

[0060] Furthermore, this method is an improvement within a mature smooth inversion framework, without completely overturning existing inversion processes and algorithms, making it easy to integrate and implement. Attached Figure Description

[0061] The accompanying drawings, which are included to provide a further understanding of the invention and form part of this invention, illustrate exemplary embodiments of the invention and are used to explain the invention, but do not constitute an undue limitation of the invention. In the drawings:

[0062] Figure 1 This is a flowchart of a high boundary sensitivity three-dimensional inversion method based on an electromagnetic gradient constraint matrix according to an embodiment of the present invention;

[0063] Figure 2 To simplify the schematic diagram of constructing the electromagnetic gradient constraint matrix under mesh construction, where, Figure 2 (a) is a conductivity model generated by a simple mesh. Figure 2 (b) is a schematic diagram of the constraint matrix corresponding to the model response, taken as a unit diagonal matrix. Figure 2 (c) is the electromagnetic gradient constraint matrix corresponding to the model response. Value selection diagram;

[0064] Figure 3 For four specific Weighted function under parameters The characteristics of variation with layer depth and skin depth, among which, Figure 3 (a) is The changing characteristics over time, Figure 3 (b) is The changing characteristics over time, Figure 3 (c) is The changing characteristics over time, Figure 3 (d) is The changing characteristics over time. Detailed Implementation

[0065] To enable those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present application, and not all of them. Based on the embodiments of the present application, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present application. It should be noted that, unless otherwise specified, the embodiments and features in the embodiments of the present application can be combined with each other.

[0066] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0067] The terms “comprising” and “having”, and any variations thereof, in the specification and claims of this application are intended to cover non-exclusive inclusion, for example, a process, method, product, or apparatus that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such process, method, product, or apparatus.

[0068] This invention provides an innovative electromagnetic field gradient-constrained inversion method. This method directly incorporates electromagnetic gradient information into the inversion process, organically integrating its boundary enhancement characteristics into a robust inversion constraint framework. While ensuring model stability, it significantly improves the resolution of geological structure boundaries. This method effectively combines the advantages of gradient imaging-based boundary identification with the quantitative rigor of inversion techniques. It provides a solid theoretical foundation for the quantitative interpretation of geophysical data and constructs a practically feasible solution. Compared to traditional apparent resistivity imaging methods and conventional inversion techniques relying on smoothing constraints, this method significantly improves boundary identification accuracy and effectively enhances the balance between data fitting consistency and model geological rationality.

[0069] Reference Figure 1 , Figure 1This is a flowchart of a high boundary sensitivity 3D inversion method based on an electromagnetic gradient constraint matrix according to an embodiment of the present invention, as shown below. Figure 1 As shown, the method includes:

[0070] S1, Construct the initial inversion model and objective function;

[0071] The initial inversion model is an empirically set initial resistivity model, which is a preliminary estimate of the underground geological structure and serves as the starting point for the inversion iteration, providing a starting point for subsequent inversion iterations.

[0072] The objective function is the function that needs to be minimized during the inversion process. In this invention, it consists of three parts: a data residual term, a reference model constraint term, and a model roughness term based on electromagnetic gradient constraints.

[0073] Data residuals: measure the difference between model-predicted data and observed data;

[0074] Reference model constraints: control the overall range of the model and prevent excessive deviation of model parameters;

[0075] Model roughness term based on electromagnetic gradient constraint: a key regularization component in the objective function. By dynamically adjusting the regularization intensity of each region in the model, it enhances the sensitivity to boundary regions, adjusts the constraint intensity of different regions, and improves the solvability and boundary recovery capability of the inversion problem.

[0076] In this embodiment of the invention, the objective function is defined as follows: ,in, This is an inversion model (usually the logarithm of conductivity or resistivity). For data residuals, For model structure items, To control the two balancing parameters, ,in, For reference model constraint terms, This is the roughness term of the model based on electromagnetic gradient constraints.

[0077] Data residuals Defined as:

[0078] ;

[0079] in, This is a diagonal weighted matrix that includes the uncertainty of the measured electromagnetic data. To measure electromagnetic data, This is the forward electromagnetic response.

[0080] Reference model constraint terms Defined as:

[0081] ;

[0082] in, To control the weighting coefficients for the reference model, As the reference model weight matrix, This is a reference model.

[0083] Model roughness term based on electromagnetic gradient constraint Defined as:

[0084] ;

[0085] In the formula, These are the gradient roughness weight control coefficients; The electromagnetic gradient constraint matrix; For model gradient operators; This refers to the electromagnetic gradient constraint regularization term.

[0086] By adjusting The value of the middle element can limit the drastic spatial changes of the model to achieve global smooth control. At the same time, by introducing the electromagnetic gradient constraint matrix, differential regulation can be implemented in local abnormal boundary regions. The regularization intensity can be appropriately relaxed in high gradient regions, thereby strengthening boundary features and enhancing structural analysis capabilities.

[0087] S2, construct the electromagnetic gradient constraint matrix in the data space;

[0088] The electromagnetic gradient constraint matrix is ​​calculated from observation data and is used to quantify the electromagnetic gradient characteristics at different locations and frequencies. Specifically, S2 includes:

[0089] S21, Calculate the electromagnetic gradient anomaly response of the geodetic model containing the anomaly at each measurement frequency at each measurement point;

[0090] The anomaly is reflected as a sudden change in the boundary, that is, a change in the electromagnetic response gradient.

[0091] measuring point Frequency measurement Electromagnetic gradient anomaly response From the measuring point Frequency measurement The spatial derivative of the electromagnetic field is obtained at the measuring point. Frequency measurement The spatial derivative of the electromagnetic field under these conditions is defined as:

[0092] ;

[0093] in, For measuring points Frequency measurement The amplitude of the secondary electromagnetic field under the condition; For measuring points +1 measurement frequency The amplitude of the secondary electromagnetic field under the condition; For measuring points +1 is the x-coordinate in a three-dimensional coordinate system (this three-dimensional coordinate system has the center of the traverse source as the origin, the direction of the survey line distribution as the x-axis, the vertical distance of the survey line from the center of the traverse source as the y-axis, and the vertical upward direction of the survey area as the z-axis. The traverse source is the long traverse emission source used for measurement. The survey line consists of measurement points distributed on the same straight line. The survey area is the measurement area that includes all survey lines). For measuring points The x-coordinate in a three-dimensional coordinate system; For measuring points The set of nearest measuring points; Let j be the neighboring measurement point number; For measuring points and The Euclidean distance between them; Two measuring points (measuring points) and The absolute difference in electromagnetic field amplitude between ).

[0094] The amplitude of the secondary electromagnetic field is the value of the secondary field calculated by the three-dimensional forward modeling method based on the total field and secondary field separation algorithm. It is the anomalous response caused by the anomalous body. The advantage of using the secondary field instead of the total field is that the secondary field mainly reflects the anomalous body itself. Therefore, using the secondary field can avoid the influence of the background field gradient, resulting in higher anomalous resolution.

[0095] S22, combine the abnormal responses at all frequencies at each measuring point to obtain the electromagnetic gradient abnormal response vector for each measuring point;

[0096] measuring point Electromagnetic gradient anomaly response vector at all frequencies Defined as:

[0097] ;

[0098] in, For measuring points Anomaly response of electromagnetic gradient at the first measurement frequency; For measuring points Anomaly response of electromagnetic gradient at the second measurement frequency; For measuring points First Anomaly response of electromagnetic gradient at a measurement frequency; It is the set of real numbers; This represents the total number of measurement points; Number the measurement points; For measuring points The total number of corresponding electromagnetic gradient anomaly response data, i.e. the total number of measured frequencies.

[0099] Furthermore, by combining the electromagnetic gradient anomaly response vectors of all measurement points, the overall electromagnetic gradient response in the data space is obtained:

[0100] ;

[0101] in, represents the overall electromagnetic gradient response in the data space, where each element represents the electromagnetic gradient anomaly response vector at all frequencies at the measurement point; This represents the electromagnetic gradient anomaly response vector at all frequencies at measurement point 1. This represents the electromagnetic gradient anomaly response vector at all frequencies at measurement point 2; For measuring points Electromagnetic gradient anomaly response vector at all frequencies.

[0102] S23, normalize the electromagnetic gradient anomaly response vector of each measuring point to obtain the normalized observed gradient anomaly vector of each measuring point;

[0103] Find the maximum value of the spatial derivative of the electromagnetic field (i.e., the electromagnetic gradient anomaly response) at all frequencies and all measurement points. To avoid division by zero errors, a preset minimum value is used. The electromagnetic gradient anomaly response vector at each measuring point is normalized to obtain the normalized observed gradient anomaly vector. The normalized observed gradient anomaly vector is defined as :

[0104] .

[0105] S24. Combine the normalized observation gradient anomaly vectors of all measurement points into a matrix to obtain the overall normalized observation gradient anomaly matrix, and construct the electromagnetic gradient constraint matrix in the data space based on the overall normalized observation gradient anomaly matrix.

[0106] Overall normalized observation gradient anomaly matrix Defined as:

[0107] ;

[0108] in, This is the normalized observation gradient anomaly vector for measurement point 1; This is the normalized observation gradient anomaly vector for measurement point 2; For measuring points The normalized observed gradient anomaly vector;

[0109] Based on the overall normalized observation gradient anomaly matrix Construct the electromagnetic gradient constraint matrix in the data space:

[0110] ;

[0111] in, This represents the electromagnetic gradient constraint matrix in the data space.

[0112] S3, which correlates frequency and layer depth through skin depth;

[0113] measuring point exist Frequency Layer depth weighting function Defined as:

[0114] ;

[0115] in, For measuring points Model No. Layer center Coordinate values; It is a constant, controlling the degree to which the weights decay with depth mismatch; Permeability in free space; Electrical conductivity; Number the floors; For measuring points Model No. Layer center Coordinate values.

[0116] Weighting function It can achieve the control of the induction response intensity of different depth layers, so that the electromagnetic gradient constraint can be accurately mapped from the observation data space to the model space, reflecting the sensitivity of electromagnetic gradients of different frequencies to different depths, and improving the depth consistency and physical rationality of the inversion.

[0117] S4, Establish the electromagnetic gradient constraint matrix in the model space;

[0118] The electromagnetic gradient constraint matrix in the data space is transformed into the model space using a frequency-depth weighted mapping function, and a normalized electromagnetic gradient constraint matrix in the model space is constructed:

[0119] ;

[0120] ;

[0121] in, The electromagnetic gradient constraint matrix is ​​normalized for the model space, which assigns low weights to high gradient regions and high weights to low gradient regions to achieve boundary-aware adjustment. It is the identity matrix; It is a frequency-depth weighted mapping matrix; This is the frequency-depth weighted mapping matrix at measurement point 1; This is the frequency-depth weighted mapping matrix at measurement point 2; For measuring points Frequency-depth weighted mapping matrix at the location;

[0122] Among them, measuring points Frequency-depth weighted mapping matrix at [location] for:

[0123] ;

[0124] For measuring points The weighting function for the first frequency at the first depth; For measuring points The weighting function for the first frequency at the second depth; For measuring points In the The weighting function of the first frequency at the layer depth; For measuring points The weighting function for the second frequency at the first depth; For measuring points The weighting function for the second frequency at the second depth; For measuring points In the The weighting function for the second frequency at layer depth; For measuring points At the depth of the first layer A weighted function of frequencies; For measuring points At the second layer depth A weighted function of frequencies; For measuring points In the Layer depth A weighted function of frequencies; Total number of floors;

[0125] S5, iteratively update the model;

[0126] In each iteration, the model parameters are updated by optimizing the objective function. The step value for the next iteration is obtained by solving the following equation:

[0127] ;

[0128] ;

[0129] ;

[0130] In the above formula, For the first The left-hand side of the next iteration inversion equation, i.e., the coefficient matrix; For the first The model obtained from the first step of calculation and the first step The difference between the models obtained from the first step; For the first The right-hand side of the inversion equation in the next iteration, i.e., the gradient direction; The orthogonal Jacobian matrix (the orthogonal Jacobian matrix describes the degree of influence of changes in model parameters on the model's predicted data; it is a matrix connecting the model space and the data space); For the first Model regularization parameters in step calculation; For the first The model predicts data in the next iteration; For the first The model obtained by step calculation ( , , , , , (The explanation can be found in the steps above.)

[0131] The model change was calculated iteratively. The final model can be obtained when the iteration stops. To achieve electromagnetic inversion.

[0132] Reference Figure 2 , Figure 2 The electromagnetic gradient constraint matrix is ​​constructed. A diagram illustrating the generation of a simple mesh: Figure 2 (a) is a conductivity model generated by a simple mesh. Figure 2 (b) is a schematic diagram of the constraint matrix corresponding to the model response, taken as a unit diagonal matrix. Figure 2 (c) is the electromagnetic gradient constraint matrix corresponding to the model response. Value illustration.

[0133] Electromagnetic gradient constraint matrix In the medium, the value of the element at the boundary of the medium is smaller than that of other diagonal elements (the value is 1 in the homogeneous region).

[0134] By adjusting the electromagnetic gradient constraint matrix This allows the mesh values ​​at the corresponding medium boundary to be smaller than the corresponding mesh values ​​in the identity matrix, thereby reducing the electromagnetic gradient constraint matrix. Used as a modifier, it increases the contribution of model changes in high gradient regions, highlights the boundary regions with drastic changes in electromagnetic fields, and improves the boundary sensitivity in the roughness vector; in weak gradient regions, the contribution of the constraint matrix is ​​small, thus preserving the natural smoothness of homogeneous regions.

[0135] Figure 3 Four different Weight function under different values Characteristics of variations in numbers at different layer depths and skin depths:

[0136] Assume the background resistivity is ; Take 40 logarithmically uniform frequencies from 1Hz to 1000Hz, corresponding to depths from -1500m to -50m, with layer depths from 0m to -1500m, each layer spaced 50m apart; when Values Weight function at time like Figure 3 As shown in (a)~(d);

[0137] The results show that each frequency has a strong influence on the skin depth corresponding layer;

[0138] The smaller the value, the smoother the correlation between skin depth (frequency) and layer depth.

[0139] This invention proposes an innovative electromagnetic field gradient-constrained inversion method that integrates electromagnetic gradient data into the inversion process and embeds its boundary enhancement properties into a robust inversion regularization framework. This method combines the qualitative advantages of gradient imaging with the quantitative rigor of inversion methods, providing theoretical support and a feasible solution for the refined interpretation of geophysical electromagnetic data.

[0140] By introducing an electromagnetic field gradient constraint matrix, dynamic adjustment of the weights of each grid cell variation in the model roughness term is achieved, enhancing boundary sensitivity while maintaining the continuity and stability of the uniform region. This mechanism effectively improves the sensitivity of the inverted objective equation to the electrical boundary structure, significantly improves the resolution of the target boundary, and overcomes the limitations of traditional apparent resistivity imaging methods and conventional smoothing constraint inversion in terms of boundary resolution capabilities.

[0141] By incorporating the horizontal variation of the electromagnetic field into the constraint system, the response characteristics of the structural change region are enhanced, thereby improving the model resolution while increasing the identification accuracy of complex geological targets.

[0142] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A high boundary sensitivity three-dimensional inversion method based on an electromagnetic gradient constraint matrix, characterized in that, include: An initial inversion model and an objective function are constructed. The initial inversion model is an initial resistivity model set based on experience, which is the starting point for the inversion iteration. The objective function consists of three parts: a data residual term, a reference model constraint term, and a model roughness term based on electromagnetic gradient constraints. Calculate the electromagnetic gradient anomaly response of the geodetic model containing the anomaly at each measurement frequency at each measurement point; By combining the abnormal responses at all frequencies at each measuring point, the electromagnetic gradient abnormal response vector of each measuring point is obtained; The electromagnetic gradient anomaly response vector at each measuring point is normalized to obtain the normalized observed gradient anomaly vector at each measuring point. The normalized observation gradient anomaly vectors of all measurement points are combined into a matrix to obtain the overall normalized observation gradient anomaly matrix, and the electromagnetic gradient constraint matrix in the data space is constructed based on the overall normalized observation gradient anomaly matrix. The electromagnetic gradient constraint matrix in the data space is transformed into the model space through a frequency-depth weighted mapping matrix, wherein the frequency-depth weighted mapping matrix is ​​used to reflect the sensitivity of each frequency electromagnetic gradient to different depths. Based on the electromagnetic gradient constraint matrix after transformation to the model space, an iterative optimization algorithm is used to solve the objective function. In each iteration, the data fitting error and structure term are jointly optimized and minimized. The model iteration step term is updated from the initial inversion model until the model converges.

2. The high boundary sensitivity three-dimensional inversion method based on the electromagnetic gradient constraint matrix according to claim 1, characterized in that, The objective function is defined as follows: ； in, Let be the objective function. For inversion model, For data residuals, For model structure items, Parameters for controlling the balance of two items; ,in, For reference model constraint terms, This is the roughness term of the model based on electromagnetic gradient constraints; The data residuals Defined as: ； in, This is a diagonal weighted matrix that includes the uncertainty of the measured electromagnetic data. To measure electromagnetic data, Forward electromagnetic response; The reference model constraint terms Defined as: ； in, To control the weighting coefficients for the reference model, As the reference model weight matrix, For reference model; The roughness term of the model based on electromagnetic gradient constraints Defined as: ； ； in, These are the gradient roughness weight control coefficients; This is the electromagnetic gradient constraint regularization term; The electromagnetic gradient constraint matrix; This is the model gradient operator.

3. The high boundary sensitivity three-dimensional inversion method based on the electromagnetic gradient constraint matrix according to claim 2, characterized in that, The calculation of the electromagnetic gradient anomaly response of the geodetic model containing the anomaly at each measurement frequency at each measurement point includes: For each measurement frequency at each measuring point, the electromagnetic gradient anomaly response at that frequency is represented by the electromagnetic field spatial derivative. This electromagnetic field spatial derivative is calculated using the Euclidean distance between adjacent measuring points and the difference in electromagnetic field amplitude, as shown in the following formula: ； in, For measuring points Frequency measurement Abnormal response to electromagnetic gradient; For measuring points Frequency measurement The amplitude of the secondary electromagnetic field under the condition; For measuring points +1 measurement frequency The amplitude of the secondary electromagnetic field under the condition; and They are measuring points +1 and measuring point The x-coordinate in a three-dimensional coordinate system; For measuring points The set of nearest measuring points; Let j be the neighboring measurement point number; For measuring points and The Euclidean distance between them; For measuring points and The absolute difference in electromagnetic field amplitude between them.

4. The high boundary sensitivity three-dimensional inversion method based on the electromagnetic gradient constraint matrix according to claim 3, characterized in that, By combining the anomalous responses at all frequencies at each measuring point, the electromagnetic gradient anomalous response vector for each measuring point is obtained, expressed as: ； in, For measuring points The electromagnetic gradient anomaly response vector; For measuring points Anomaly response of electromagnetic gradient at the first measurement frequency; For measuring points Anomaly response of electromagnetic gradient at the second measurement frequency; For measuring points First Anomaly response of electromagnetic gradient at a measurement frequency; It is the set of real numbers; This represents the total number of measurement points; For measuring points The total number of corresponding electromagnetic gradient anomaly response data, i.e. the total number of measured frequencies.

5. The high boundary sensitivity three-dimensional inversion method based on the electromagnetic gradient constraint matrix according to claim 4, characterized in that, The electromagnetic gradient anomaly response vector at each measuring point is normalized to obtain the normalized observed gradient anomaly vector at each measuring point, which includes: Find the maximum value of the electromagnetic gradient anomaly response at all measurement points and frequencies, and use a preset minimum value; The larger of the maximum value and the preset minimum value is selected to normalize the electromagnetic gradient anomaly response vector of each measuring point: ； in, For measuring points The normalized observed gradient anomaly vector; This represents the maximum value of the electromagnetic gradient anomaly response at all measurement points and frequencies; This is a preset minimum value to avoid division by zero errors.

6. The high boundary sensitivity three-dimensional inversion method based on the electromagnetic gradient constraint matrix according to claim 5, characterized in that, The normalized observation gradient anomaly vectors of all measurement points are combined into a matrix to obtain the overall normalized observation gradient anomaly matrix. Based on this overall normalized observation gradient anomaly matrix, an electromagnetic gradient constraint matrix in the data space is constructed, including: By combining the normalized observation gradient anomaly vectors of all measurement points, the overall normalized observation gradient anomaly matrix is ​​obtained: ； in, This is the overall normalized observation gradient anomaly matrix; This is the normalized observation gradient anomaly vector for measurement point 1; This is the normalized observation gradient anomaly vector for measurement point 2; For measuring points The normalized observed gradient anomaly vector; The overall normalized observation gradient anomaly matrix is ​​then used. As an electromagnetic gradient constraint matrix in the data space: ； in, This represents the electromagnetic gradient constraint matrix in the data space.

7. The high boundary sensitivity three-dimensional inversion method based on the electromagnetic gradient constraint matrix according to claim 6, characterized in that, Transforming the electromagnetic gradient constraint matrix in the data space to the model space using a frequency-depth weighted mapping matrix includes: Define the frequency-depth weighted mapping matrix corresponding to each measurement point as follows: ； in, For measuring points Frequency-depth weighted mapping matrix at the location; For measuring points The weighting function for the first frequency at the first depth; For measuring points The weighting function for the first frequency at the second depth; For measuring points In the The weighting function of the first frequency at the layer depth; For measuring points The weighting function for the second frequency at the first depth; For measuring points The weighting function for the second frequency at the second depth; For measuring points In the The weighting function for the second frequency at layer depth; For measuring points At the depth of the first layer A weighted function of frequencies; For measuring points At the second layer depth A weighted function of frequencies; For measuring points In the Layer depth A weighted function of frequencies; This represents the total number of floors. Each element in the frequency-depth weighted mapping matrix is ​​defined as follows: ； In the formula, For measuring points exist Frequency Weighting function for layer depth; This is a constant used to control the severity of the weight decay with depth mismatch; For measuring points Model No. Layer center Coordinate values; Permeability in free space; Electrical conductivity; Number the floors; For measuring points Model No. Layer center Coordinate values; The electromagnetic gradient constraints in the data space are transformed into the model space using the following formula: ； ； in, The electromagnetic gradient constraint matrix is ​​normalized in the model space to achieve boundary-aware adjustment; It is the identity matrix; It is a frequency-depth weighted mapping matrix; This is the frequency-depth weighted mapping matrix at measurement point 1; This is the frequency-depth weighted mapping matrix at measurement point 2; For measuring points The frequency-depth weighted mapping matrix at that location.

8. The high boundary sensitivity three-dimensional inversion method based on the electromagnetic gradient constraint matrix according to claim 7, characterized in that, Solving the objective function using an iterative optimization algorithm includes: No. The step value for the next iteration is obtained by solving the following equation: ； ； ； In the formula, For the first The left-hand side of the next iteration inversion equation, i.e., the coefficient matrix; For the first The model obtained from the first step of calculation and the first step The difference between the models obtained from the first step; For the first The right-hand side of the inversion equation in the next iteration, i.e., the gradient direction; For orthogonal Jacobian matrix; For the first Model regularization parameters in step calculation; For the first The model predicts data in the next iteration; For the first The model obtained through step-by-step calculation.

Images