A small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization

Abstract

The present application belongs to the technical field of geophysical electromagnetic detection and inversion, and particularly is a small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization. The method comprises the following steps: S1: obtaining of multi-frequency and multi-offset secondary field data and establishment of observation error; S2: shallow encryption of non-uniform grid and parameterization of multi-unit logarithmic resistivity; S3: screening of multiple uniform initial models and pilot Ockham scanning; S4: normalized data Jacobian and sensitivity adaptive regularization; S5: sensitivity adaptive rescan and moderate IRLS block optimization; and S6: approximate model resolution and reliable interpretation depth. The present application can improve the discrete precision of near-surface structure while controlling the size of inversion parameters, reduce the influence of initial resistivity and regularization parameter selection on the inversion result, and obtain an Ockham model meeting the target fitting requirement without needing to accurately give the real stratum interface in advance.

Description

Technical Field

[0001] This invention relates to the field of geophysical electromagnetic detection and inversion technology, specifically a small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization. Background Technology

[0002] Frequency-domain electromagnetic methods transmit periodic electromagnetic fields and measure the secondary electromagnetic field response induced by the underground medium. Based on the variation of the response with frequency, transmit / receive offset, and underground electrical structure, the conductivity or resistivity distribution of the underground medium can be inferred. Small-loop frequency-domain electromagnetic systems are characterized by their small size, controllable frequency, lack of need for grounding electrodes, and suitability for urban roads and narrow spaces. They can be used to identify backfill disturbance layers, water-rich areas, seepage channels, clay-rich areas, and other low-resistivity anomalies.

[0003] In urban roads and shallow underground spaces, water-bearing backfill layers, seepage channels, clay layers, and adjacent artificial conductive facilities often create significant shallow low-resistivity responses. These responses are typically stronger than those of the underlying strata and alter the sensitivity of data at different frequencies and offsets to the depth direction. Multi-frequency and multi-offset observations can provide complementary information, but the constraint capacity of different observation data on each depth model unit is not uniform.

[0004] One-dimensional multi-element inversion does not require pre-determining the actual number of stratigraphic layers and interfaces, and can describe continuous electrical variations through a large number of thin-layer elements, making it suitable for scenarios where prior knowledge of underground structures is insufficient. However, as the number of model elements increases, the ill-conditioned nature and multiple solutions of the inversion problem intensify, necessitating a regularized stable solution. Traditional fixed regularization typically sets constraints based only on depth or empirical parameters, failing to fully utilize the differences in actual sensitivity reflected by the current data and the current model.

[0005] Existing technologies have the following main shortcomings: First, when using a uniform grid, if the grid is coarse, it is difficult to discretize the near-surface thin layer and shallow low-resistivity anomalies; if a fine grid is used throughout the entire depth, the number of model parameters becomes too large, the computational load increases, and the deep stability decreases. Second, when using uniform smoothing constraints, the shallow high-sensitivity region and the deep low-sensitivity region are subject to similar constraint strengths, making it difficult to simultaneously consider shallow resolution and deep stability. Third, fixed reference model constraints cannot change with data coverage, and when the deep sensitivity is low, non-physical oscillations deviating from the reference model are prone to occur. Fourth, ordinary L2 smoothing inversion will cause continuous tailing of electrical boundaries, while directly using strong block constraints may sacrifice data fitting and lead to over-sharpening. Fifth, some methods only evaluate the results based on data fitting error, lacking diagnostic information such as model element sensitivity, model resolution, and reliable interpretation depth. Therefore, this invention provides a small-loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization. Summary of the Invention

[0006] The purpose of this section is to outline some aspects of the embodiments of the present invention and to briefly describe some preferred embodiments. Simplifications or omissions may be made in this section, as well as in the abstract and title of this application, to avoid obscuring the purpose of these documents; however, such simplifications or omissions should not be construed as limiting the scope of the invention.

[0007] To address the aforementioned technical problems, according to one aspect of the present invention, the present invention provides the following technical solution:

[0008] A small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization includes the following steps:

[0009] S1: Use a single small loop transmitter and a small loop receiver to acquire the real part and imaginary part of the secondary magnetic field, expand the observations into an observation data vector, and construct the observation standard deviation vector;

[0010] S2: Based on the expected detection depth, establish a non-uniform one-dimensional grid with shallow densification and deep sparseness. Use the common logarithm of resistivity of each model unit as the inversion parameter to construct first-order and second-order depth difference matrices and reference model constraints.

[0011] S3: Set up multiple uniform initial resistivity models, select the initial model based on the fitting error of normalized data, and perform pilot scans on multiple smoothing regularization parameters.

[0012] S4: Calculate the Jacobian matrix of the predicted data after normalization of the observation standard deviation at the leading Occam's calculus model, and adjust the first-order and second-order depth difference constraints and the reference model constraints respectively.

[0013] S5: Scan the smoothing regularization parameter again to obtain the sensitivity adaptive smoothing model; update the iterative reweighting coefficients according to the first difference of the model, and gradually reduce the smoothing stability parameter to carry out iterative reweighted least squares block optimization; when the normalized fitting error of the block optimization model exceeds the allowable value, reduce the block regularization parameter and solve it again to obtain the block multi-unit resistivity model.

[0014] S6: Calculate the approximate model resolution, combine the relative sensitivity of the model elements and the resolution diagonal value to determine the reliable interpretation depth index, and output the blocky multi-element resistivity model, the relative sensitivity of the model elements, the resolution diagonal value, the data fitting and the normalized residual.

[0015] As a preferred embodiment of the small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization described in this invention, the specific method of S1 is as follows: a small loop transmitter and a small loop receiver are arranged in the detection area. The transmitting loop is sequentially supplied with periodic currents at multiple transmission frequencies. The receiving loop obtains the secondary magnetic field response generated by the underground medium at multiple transmit-receive offset positions. The received signal is decomposed into real response and imaginary response with the phase of the transmitting current as a reference.

[0016] Let the set of transmission frequencies be . In the formula, For the first One transmission frequency, This represents the total number of transmission frequencies.

[0017] Let the set of transmit and receive offsets be: In the formula, Let q be the transmit / receive offset. This represents the total transmit and receive offset.

[0018] When the emission ray is the radius When the small circular loop appears, the first... The equivalent magnetic dipole moment at each transmission frequency is: In the formula, For the equivalent magnetic dipole moment, For the number of turns of the transmit loop, For the first... One transmission frequency and the Transmit / receive offset Let the complex response of the second magnetic field in the sensitive direction of the receiving loop be... The direct first-order magnetic field complex response is ;

[0019] When expressed in parts per million, the real and imaginary observations are respectively: , ;

[0020] Following the order of real parts first, then imaginary parts within each frequency-offset combination, all observations are expanded as follows: The total number of individual observations is: ;set up for The first in One single observation, of which ;

[0021] Set the real part fixed absolute error base respectively and the imaginary part with a fixed absolute error base The standard deviations of the real and imaginary observations corresponding to the i-th transmit frequency and the q-th transmit / receive offset are respectively: , In the formula, The relative error coefficient;

[0022] Arrange the standard deviations in the same order according to the observed data vectors. And construct the observation error weighting matrix: .

[0023] As a preferred embodiment of the small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization described in this invention, the standard deviation is used to normalize the data residuals, normalize the data Jacobian, evaluate the model fit, and calculate the model resolution. It does not require adding random disturbances to the observation data, nor does it change the original observation values. The fixed absolute error floor is determined based on the instrument background, field test, repeated observation error, and field noise level.

[0024] As a preferred embodiment of the small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization described in this invention, the specific method of S2 is that the non-uniform one-dimensional grid adopts a smaller model unit thickness near the ground surface, and the model unit thickness is gradually increased with the depth.

[0025] A uniform half-space is set at the bottom, and the number of model elements with finite thickness is [number missing]. If a half-space is set at the bottom, the total number of model parameters is: ;No. The resistivity of each model element is denoted as The inversion parameters are defined as follows: The model vector is: ;

[0026] A first-order fundamental difference matrix is ​​established based on the center-to-center distance of the model units and the depth-based fundamental weights. Second-order fundamental difference matrix and the basic reference weight matrix The first-order difference is used to suppress non-physical jumps between adjacent model units, the second-order difference is used to suppress high-frequency oscillations formed by three consecutive model units, and the reference model constraint is used to limit the low-sensitivity parameters from deviating excessively from the initial background.

[0027] As a preferred embodiment of the small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization described in this invention, the specific method of S3 is as follows: setting Uniform initial resistivity ,in For each uniform model, the predicted data are calculated using a one-dimensional small loop frequency-domain electromagnetic forward modeling operator. And calculate the normalized root mean square error: ;

[0028] The uniform model with the smallest normalized root mean square error was selected as the initial model. and reference model For candidate smoothing regularization parameters The joint residual vector for smooth inversion is constructed as follows: It is formed by sequentially concatenating the normalized data residuals, first-order difference residuals, second-order difference residuals, and reference model constraint residuals. After expansion... ,in, The first in the joint residual vector Each scalar residual component , The total number of elements in the joint residual vector is given by the formula. The second-order difference relative weights are used for the smoothing phase; These are the constraint coefficients for the reference model during the smoothing phase.

[0029] A nonlinear least squares method with parameterized upper and lower bounds is used to solve the problem, and a soft L1 robust loss is employed to reduce the dominance of a few large residuals on model updates: the soft L1 robust loss introduces a degree parameter. If the value is greater than 0, the smoothing inversion objective function is uniformly written as: In the formula, Let be the total number of elements in the joint residual vector, where The scale controlling the transition of residuals from approximately quadratic penalty to approximately linear penalty; when the absolute value of the residuals is much smaller than... When the residual value is close to the quadratic loss, the dominance of a single large residual on the model update is reduced when the absolute value of the residual is large.

[0030] Sensitivity adaptive rescan and IRLS trial calculation use the same soft L1 scale parameter, and scan multiple times sequentially. Candidate values; when multiple candidate models satisfy When, choose Larger candidate models are used as leading Occam models Prioritize obtaining the smoothest model that meets the data fitting requirements; otherwise, select the candidate model with the smallest normalized root mean square error; if no candidate model meets the target error, select the candidate model with the smallest normalized root mean square error.

[0031] As a preferred embodiment of the small-loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization described in this invention, the specific method of S4 is as follows:

[0032] In the leading Occam model Calculate the normalized data Jacobian matrix. , its first line, number The column represents the rate of change of the k-th observation in the observation data vector relative to the j-th logarithmic resistivity parameter: ;

[0033] When using finite difference: In the formula, For the first Finite difference perturbation of each model parameter; For the first A unit direction vector with one element being 1 and the rest being 0. If the positive perturbation exceeds the upper limit of the parameters, then the negative perturbation is used.

[0034] No. The original sensitivity of each model unit is: A moving average operator along the depth direction is used. Smooth the original sensitivity and normalize the maximum value: ;

[0035] Set the relative sensitivity for truncation: ,in, As the lower limit of relative sensitivity, the first-order difference constraint is adjusted according to the truncated relative sensitivity of adjacent model units, the second-order difference constraint is adjusted according to the truncated relative sensitivity of three consecutive model units, and the reference model constraint is adjusted according to the truncated relative sensitivity of a single model unit.

[0036] The first-order difference constraint adjustment factor is: The second-order difference constraint adjustment factor is: The reference model constraint adjustment factor is: In the formula, To adjust the intensity of differential constraint sensitivity, The sensitivity adjustment intensity is constrained by the reference model;

[0037] This yields the first-order difference matrix, the second-order difference matrix, and the reference model weight matrix for sensitivity adaptation. , , The lower the relative sensitivity, the stronger the corresponding constraint; the higher the relative sensitivity, the closer the corresponding constraint is to the basic constraint.

[0038] use , and The regularization parameter is scanned again near the leading Occam's parameter, and the sensitivity adaptive smoothing model is obtained according to the selection rule of S3. and its regularization parameters .

[0039] As a preferred embodiment of the small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization described in this invention, the specific method of S5 is as follows:

[0040] by As the initial model for block optimization, the first In the next reweighted iteration, the first-order difference gradient is calculated: , No. The reweighting coefficients for each first-order difference term are: In the formula, For the first The smoothing and stabilizing parameters of the next iteration; and These are the lower and upper limits of the reweighted coefficients, respectively.

[0041] The smoothing stability parameter decreases successively according to the following formula: The joint residual vector for block optimization is: In the formula, For block regularization parameters; The relative weights are the second-order difference of the block stage; The constraint coefficients of the block stage reference model are given in the formula. The decreasing coefficient, To provide a lower limit for smooth and stable parameters;

[0042] Let the normalized root mean square error of the sensitivity adaptive smoothing model be... The target fitting error is The allowable error for the block model is: ,in, The fitting error protection coefficient is set when the normalized root mean square error obtained in a trial calculation is greater than 1. At that time, Multiply by an attenuation coefficient less than 1 The solution was then recalculated, and after a preset number of iterations, a block-shaped multi-element model was obtained. .

[0043] As a preferred embodiment of the small-loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization described in this invention, the specific method of step S6 is as follows:

[0044] Normalized data Jacobian matrix at the leading Occam's model And the approximate model resolution is calculated using the regularization matrix of the selected sensitivity adaptive smoothing model. The regularization matrix is: The approximate model resolution matrix is: In the formula, For numerical stability damping, for The first-order identity matrix, the second-order identity matrix The resolution diagonal value of each model parameter is: Let the first The bottom boundary depth of a finite thickness model element is The relative sensitivity threshold is The resolution diagonal threshold is The reliable explanation depth metric is defined as: If there is no finite thickness model element that simultaneously satisfies the threshold, then Take zero.

[0045] As a preferred embodiment of the small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization described in this invention, the blocky multi-unit resistivity model is used to identify underground electrical changes and low-resistivity anomaly ranges, and the reliable interpretation depth index is only used for result diagnosis and interpretation marking, and is not used as a pre-given true layer thickness, true interface or deep target location.

[0046] Compared with the prior art, the beneficial effects of the present invention are:

[0047] (1) Using a shallow-densified and deep-seated non-uniform grid, and combining multiple uniform initial models and multiple regularization parameters to carry out pilot scanning, the discretization accuracy of near-surface structures can be improved while controlling the scale of inversion parameters, reducing the influence of initial resistivity and regularization parameter selection on inversion results, and obtaining an Occam model that meets the target fitting requirements without needing to accurately give the real stratigraphic interface in advance.

[0048] (2) The relative sensitivity of each model unit is calculated based on the normalized data Jacobi matrix, so that the sensitivity evaluation comprehensively considers the effects of observation error, transmission frequency, transmission and reception offset and response components; the first-order difference, second-order difference and reference model constraint are adaptively adjusted according to the sensitivity, which can enhance the stability of the deep low-sensitivity region, suppress non-physical oscillations and thin-layer artifacts, while retaining the data-driven update capability of the shallow high-sensitivity region.

[0049] (3) The mild iterative reweighted least squares block optimization with gradually decreasing smoothing stability parameters and normalized data fitting error protection can reduce the smoothing tail of the electrical interface without significantly reducing the degree of data fitting. At the same time, the relative sensitivity of the model unit, the diagonal value of the approximate model resolution and the reliable interpretation depth index are output, thereby avoiding the evaluation of the credibility of the inversion results based solely on the data fitting error. Attached Figure Description

[0050] To more clearly illustrate the technical solutions of the embodiments of the present invention, the present invention will be described in detail below with reference to the accompanying drawings and detailed embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. Wherein:

[0051] Figure 1 This is a flowchart of a small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization according to the present invention;

[0052] Figure 2 This is a diagram showing the shallow densification, deep gradual sparse non-uniform inversion grid and cell thickness variation in Embodiment 2 of the small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization of the present invention. (a) is a diagram of the non-uniform inversion network, and (b) is a diagram of cell thickness variation with depth.

[0053] Figure 3 This is a comparison diagram of the synthetic reference model and the first-stage IRLS multi-unit inversion model in Embodiment 3 of the present invention, which is a small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization.

[0054] Figure 4 This is a diagnostic diagram of the relative sensitivity of the model unit and the approximate model resolution in the fourth embodiment of the small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization of the present invention, wherein (a) is the relative sensitivity of the model unit and (b) is the approximate model resolution.

[0055] Figure 5 The image shows the fitting results of the observation data and the predicted data for each transmit / receive offset in the small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization of the present invention. In this image, (a) is the fitting result of the real part data and (b) is the fitting result of the imaginary part data.

[0056] Figure 6 This is a normalized residual diagram of the observation data and prediction data of each transmit / receive offset in a small loop frequency domain electromagnetic inversion method based on sensitivity adaptive regularization according to the present invention, wherein (a) is the normalized residual of the real part and (b) is the normalized residual of the imaginary part. Detailed Implementation

[0057] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0058] Secondly, the present invention is described in detail with reference to the schematic diagrams. When detailing the embodiments of the present invention, for ease of explanation, the cross-sectional views illustrating the device structure may be partially enlarged, not according to the usual scale. Furthermore, the schematic diagrams are merely examples and should not limit the scope of protection of the present invention. In addition, actual fabrication should include three-dimensional spatial dimensions of length, width, and depth.

[0059] To make the objectives, technical solutions, and advantages of the present invention clearer, the embodiments of the present invention will be described in further detail below with reference to the accompanying drawings.

[0060] This invention proposes a small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization. The overall process is as follows: Figure 1 As shown, the process includes six steps: establishing multi-frequency, multi-offset data and observation errors; non-uniform multi-element parameterization; uniform initial model selection and leading Occam's sweep; normalized data Jacobian and sensitivity adaptive regularization; sensitivity adaptive resweep and mild iterative reweighted least squares (IRLS) block optimization; reliability diagnosis; and result output. Specifically:

[0061] S1: Acquisition of Multi-Frequency, Multi-Offset Secondary Field Data and Establishment of Observation Errors: A small-loop transmitter and a small-loop receiver are deployed in the detection area. Periodic currents are sequentially passed through the transmitter at multiple transmission frequencies, while the receiver acquires the secondary magnetic field response generated by the underground medium at multiple transmit / receive offset positions. The received signal is decomposed into real and imaginary responses with reference to the phase of the transmitted current.

[0062] Let the set of transmission frequencies be . In the formula, For the first One transmission frequency, Let be the total number of transmission frequencies. Let the set of transmit / receive offsets be: In the formula, Let q be the transmit / receive offset. This represents the total transmit and receive offset. When the transmit loop is the radius... When the small circular loop appears, the first... The equivalent magnetic dipole moment at each transmission frequency is: In the formula, Equivalent magnetic dipole moment, unit: ; This refers to the number of turns of the transmit loop; This represents the amplitude of the transmitting current. For the [missing value]... One transmission frequency and the Transmit / receive offset Let the complex response of the second magnetic field in the sensitive direction of the receiving loop be... The direct first-order magnetic field complex response is When expressed in parts per million (ppm), the real and imaginary observations are respectively: , Following the order of real parts first, then imaginary parts within each frequency-offset combination, all observations are expanded as follows: The total number of individual observations is: ;set up for The first in One single observation, of which To comprehensively consider the differences in observed amplitude, instrument background noise, and the differences in measurement accuracy between the real and imaginary responses, a fixed absolute error floor for the real part is set separately. and the imaginary part with a fixed absolute error base The fixed absolute error floor is determined based on instrument air-field testing, repeatability measurement error, receiver system background noise, or ambient noise level. The standard deviations of the real and imaginary observations corresponding to the i-th transmit frequency and the q-th transmit / receive offset are respectively: , In the formula, The standard deviations are the relative error coefficients. The standard deviations are arranged in the same order according to the observed data vectors. And construct the observation error weighting matrix: The standard deviation is used to normalize data residuals, normalize data Jacobian, evaluate model fit, and calculate model resolution. It does not require adding random disturbances to the observed data, nor does it change the original observed values. The fixed absolute error floor can be determined based on the instrument background, airfield testing, repeated observation errors, and field noise levels.

[0063] S2: Shallow Refined Non-Uniform Mesh and Multi-Element Logarithmic Resistivity Parameterization: A shallowly refined, gradually sparser non-uniform one-dimensional mesh is established based on the expected detection depth. A smaller model element thickness is used near the surface to improve the discretization accuracy of the road structure layer and shallow low-resistivity anomalies; the model element thickness increases piecewise with depth to control the number of parameters and improve the stability of deep inversion. A uniform half-space is set at the bottom. The number of model elements with a finite thickness is set to... If a half-space is set at the bottom, the total number of model parameters is: . No. The resistivity of each model element is denoted as The inversion parameters are defined as follows: The model vector is: A first-order fundamental difference matrix is ​​established based on the center spacing and depth-based fundamental weights of the model elements. Second-order fundamental difference matrix and the basic reference weight matrix First-order differences are used to suppress non-physical jumps between adjacent model elements, second-order differences are used to suppress high-frequency oscillations formed by three consecutive model elements, and reference model constraints are used to limit the excessive deviation of low-sensitivity parameters from the initial background. The difference matrix is ​​scaled according to the center-to-center spacing of the model elements to make the constraint strength of elements with different thicknesses comparable.

[0064] S3: Multiple Uniform Initial Model Selection and Leading Occam Scan: Settings Uniform initial resistivity ,in For each uniform model, the predicted data are calculated using a one-dimensional small loop frequency-domain electromagnetic forward modeling operator. And calculate the normalized root mean square error: The uniform model with the smallest normalized root mean square error is selected as the initial model. and reference model For candidate smoothing regularization parameters The joint residual vector for smooth inversion is constructed as follows: It is formed by sequentially concatenating the normalized data residuals, first-order difference residuals, second-order difference residuals, and reference model constraint residuals. After expansion... ,in, The first in the joint residual vector Each scalar residual component , The total number of elements in the joint residual vector is given by the formula. The second-order difference relative weights are used for the smoothing phase; These are the constraint coefficients for the reference model during the smoothing phase. A nonlinear least squares method with parameterized upper and lower bounds is used to solve them, and a soft L1 robust loss is employed to reduce the dominance of a few large residuals on model updates: the soft L1 robust loss introduces a degree parameter. >0. To ensure that the scale parameters in the technical solution and the embodiments completely correspond, the smoothing inversion objective function is uniformly written as: In the formula, Let be the total number of elements in the joint residual vector, where The scale controlling the transition of residuals from approximately quadratic penalty to approximately linear penalty; when the absolute value of the residuals is much smaller than... When the time approaches quadratic loss, the dominance of a single large residual on model updates is reduced when the absolute value of the residual is large. Sensitivity adaptive rescanning and IRLS trial calculations use the same soft L1 scaling parameter. Multiple scans are performed sequentially. Candidate values; when multiple candidate models satisfy When, choose Larger candidate models are used as leading Occam models Prioritize obtaining the smoothest model that meets the data fitting requirements; otherwise, select the candidate model with the smallest normalized root mean square error; if no candidate model meets the target error, select the candidate model with the smallest normalized root mean square error.

[0065] S4: Normalized Data Jacobian and Sensitivity Adaptive Regularization: In the Preceding Occam Model Calculate the normalized data Jacobian matrix. , its first line, number The column represents the rate of change of the k-th observation in the observation data vector relative to the j-th logarithmic resistivity parameter: When using finite difference: In the formula, For the first Finite difference perturbation of each model parameter; For the first A unit direction vector with one element set to 1 and the rest to 0. If the positive perturbation exceeds the upper limit of the parameters, a negative perturbation is used. The original sensitivity of each model unit is: A moving average operator along the depth direction is used. Smooth the original sensitivity and normalize the maximum value: To avoid an excessively large constraint factor due to extremely low sensitivity, a cutoff relative sensitivity is set: .in, This is the lower limit of relative sensitivity. The first-order difference constraint is adjusted based on the truncated relative sensitivity of adjacent model elements, the second-order difference constraint is adjusted based on the truncated relative sensitivity of three consecutive model elements, and the reference model constraint is adjusted based on the truncated relative sensitivity of a single model element. The first-order difference constraint adjustment factor is: The second-order difference constraint adjustment factor is: The reference model constraint adjustment factor is: In the formula, To adjust the intensity of differential constraint sensitivity, The sensitivity adjustment intensity is constrained by the reference model. This yields the first-order difference matrix, second-order difference matrix, and reference model weight matrix for sensitivity adaptation. , . The lower the relative sensitivity, the stronger the corresponding constraint; the higher the relative sensitivity, the closer the corresponding constraint is to the basic constraint. (Using...) , and The regularization parameter is scanned again near the leading Occam's parameter, and the sensitivity adaptive smoothing model is obtained according to the selection rule in step S3. and its regularization parameters .

[0066] S5: Sensitivity Adaptive Rescan and Mild IRLS Block Optimization: with As the initial model for block optimization, the first In the next reweighted iteration, the first-order difference gradient is calculated: , No. The reweighting coefficients for each first-order difference term are: In the formula, For the first The smoothing and stabilizing parameters of the next iteration; and These are the lower and upper limits of the reweighted coefficients, respectively. The smoothing stability parameter decreases successively according to the following formula: The joint residual vector for block optimization is: In the formula, For block regularization parameters; The relative weights are the second-order difference of the block stage; The constraint coefficients of the block stage reference model are given in the formula. The decreasing coefficient, To define the lower bound of the smoothing stability parameters, let the normalized root mean square error of the sensitivity adaptive smoothing model be... The target fitting error is The allowable error for the block model is: .in, The fitting error protection coefficient is set when the normalized root mean square error obtained in a trial calculation is greater than 1. At that time, Multiply by an attenuation coefficient less than 1 The solution was then recalculated, and after a preset number of iterations, a block-shaped multi-element model was obtained. .

[0067] Step S6: Approximate Model Resolution and Reliable Interpretation Depth: Using the normalized data Jacobian matrix from the leading Occam's Matrix. The approximate model resolution is calculated using the regularization matrix of the selected sensitivity adaptive smoothing model. The regularization matrix is: The approximate model resolution matrix is: In the formula, For numerical stability damping, for The first-order identity matrix, the second-order identity matrix The resolution diagonal value of each model parameter is: Let the first... The bottom boundary depth of a finite thickness model element is The relative sensitivity threshold is The resolution diagonal threshold is The reliable explanation depth metric is defined as: If no finite-thickness model element simultaneously satisfies the threshold, then Set to zero. The stated reliable interpretation depth is only used as a diagnostic indicator under the current data, grid, and regularization conditions, and is not considered as the true stratigraphic interface.

[0068] Example:

[0069] Example 1: Synthetic data and error parameters of multi-frequency and multi-offset small loop lines.

[0070] A vertical magnetic dipole equivalent circular small loop source is used, with the transmitter located at coordinates (0, 0, 0.25 m). The receivers are sequentially arranged along the x-direction at offsets of 1 m, 2 m, 4 m, 6 m, and 8 m, with both the transmission and reception directions in the z-direction. The secondary magnetic field is converted to its real and imaginary parts in ppm using the primary field as a reference. To ensure a one-to-one correspondence between the implementation parameters and the symbols in step S1, the observation system parameters for this embodiment are shown in Table 1.

[0071] Table 1: Observation System Parameters in this Embodiment

[0072]

[0073] In this embodiment, each offset uses the same real and imaginary parts to fix the absolute error base, that is, for any q, the base is taken as... =0.01 ppm =0.05 ppm. The transmission frequency uses denser logarithmic sampling in the 50 Hz to 5 kHz range and relatively sparser logarithmic sampling in the 5 kHz to 100 kHz range, with specified frequency points added; all frequencies are aligned to a 25 Hz frequency grid. Each offset contains 75 frequency points, and each frequency-offset combination contains both real and imaginary observations, resulting in 375 frequency-offset combinations and 750 individual observations. The finite layer thickness of the synthetic reference model is [0.50, 1.50, 3.00, 3.00, 4.00] m, with a resistivity of [500, 18, 160, 90, 28, 120] Ω·m; where 0.50–2.00 m is a shallow low-resistivity layer, 8.00–12.00 m is a deep low-resistivity layer, and below 12.00 m is a 120 Ω·m half-space.

[0074] Example 2: Non-uniform grid, parameter range and pilot scan.

[0075] like Figure 2 As shown, this embodiment employs a non-uniform one-dimensional inversion mesh that is shallowly refined and gradually thins out at deeper layers. Wherein, Figure 2 (a) shows the non-uniform inversion mesh generation. Figure 2(b) shows the relationship between model element thickness and depth. The model element thicknesses in the ranges of 0–1 m, 1–3 m, 3–6 m, 6–10 m, 10–20 m, and 20–30 m are 0.04 m, 0.10 m, 0.20 m, 0.50 m, 1.25 m, and 2.00 m, respectively, corresponding to 25, 20, 15, 8, 8, and 5 model elements of finite thickness. As can be seen from the figure, the near-surface grid is denser, and the model element thickness increases segment by segment with depth, thereby improving the discretization accuracy of shallow structures while controlling the scale of inversion parameters.

[0076] The non-uniform one-dimensional mesh is segmented according to Table 2. Number of elements in the finite thickness model. The total number of model parameters after adding the bottom half space The model contains 45 finite-thickness model elements in the range of 0 to 3 m. The upper and lower bounds of the model resistivity are as follows: The candidate set for uniform initial resistivity is: The candidate set of leading smoothing regularization parameters is as follows: Target normalized root mean square error The maximum number of function evaluations for each candidate model in the pilot scan is set to... The trust region reflection least squares algorithm with upper and lower bounds is adopted, and the soft L1 loss scale parameter is used. The depth-based regularization weights and reference weights in the program are based on the depth of the model unit centers. Set to: , The weights mentioned are only used as a basis for construction in this embodiment. , and The preferred implementation is not limited to other equivalent weighting forms that continuously vary with depth. In the pilot scan, the second-order difference relative weights are used during the smoothing phase. The reference model constraint coefficient is 0.025. The value is 0.010. The sensitivity adaptive smoothing model is obtained using the same Occam's selection rule as the pilot scan. and its regularization parameters .

[0077] Table 2: Piecewise parameters of non-uniform one-dimensional inversion grid

[0078]

[0079] Example 3: Sensitivity Adaptive Regularization and Mild IRLS Block Optimization.

[0080] In this embodiment, the iterative reweighted block optimization parameters are shown in Table 3, and the results are as follows: Figure 3 IRLS block multi-unit model Full data normalized root mean square error The value is 0.02950, ​​which is less than the target value. Within the shallow low-resistivity region of 0.50–2.00 m, the thickness-weighted geometric mean resistivity of the first-stage multi-element model is approximately 22.58 Ω·m; within the deep low-resistivity region of 8.00–12.00 m, it is approximately 35.29 Ω·m. Although the first-stage multi-element results still exhibit smoothing tails and are not used to directly determine precise formation interfaces, they can simultaneously reflect the low-resistivity variations in both shallow and deep layers, indicating that sensitivity adaptive regularization and mild IRLS optimization are feasible.

[0081] Table 3: Iterative Reweighted Block Optimization Parameters

[0082]

[0083] Example 4: Model resolution and reliable interpretation depth.

[0084] In this embodiment, the resolution and depth diagnostic parameters are shown in Table 4, and the data and final results are as follows: Figures 4-6 Relative sensitivity of model units The range is approximately 0.0185 to 1.0000, normalized resolution diagonal value. The range is approximately 0.0477 to 1.0000. According to... and The combined criteria and the reliable interpretation depth obtained in this embodiment The value is 30 m, which represents the bottom boundary of the finite grid. This value reaches the bottom boundary of the finite grid in this embodiment, indicating that under the conditions of synthetic data without random noise, the current error model, and relative resolution normalization, there are still model elements at a depth of 30 m that satisfy the joint threshold. See the results below. Figure 4 .

[0085] Table 4: Resolution and Depth Diagnostic Parameters

[0086]

[0087] Example 5: Practical data application and parameter setting principles.

[0088] Actual detection data according to transmission frequency Transmit / receive offset Secondary field actual part and the imaginary part of the second field Four columns arranged. Program check. , The observed values ​​are finite, and each frequency-offset combination has no duplicate records. In actual data processing, the instrument parameters should be set accordingly. , , , return line height Transmit and receive directions and response units. If repeated observation statistics are missing, [the following can be specified]: and As an initial setting, and adjusted based on open-field testing, repeatability error, and on-site noise levels; during a complete inversion process, the two error parameters remain unchanged. The model parameter range can be initially set to... , Uniform initial resistivity set Preceding regularization parameter set IRLS iteration count The diagnostic threshold can be set with reference to Examples 2 to 5, and then adjusted according to the instrument's dynamic range, the expected resistivity of the site, and the actual fitting situation.

[0089] Although the present invention has been described above with reference to embodiments, various modifications can be made and components can be replaced with equivalents without departing from the scope of the invention. In particular, as long as there is no structural conflict, the features in the disclosed embodiments can be combined with each other in any manner. The lack of an exhaustive description of these combinations in this specification is merely for the sake of brevity and resource conservation. Therefore, the present invention is not limited to the specific embodiments disclosed herein, but includes all technical solutions falling within the scope of the claims.

Claims

1. A small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization, characterized in that, Includes the following steps: S1: Use a single small loop transmitter and a small loop receiver to acquire the real part and imaginary part of the secondary magnetic field, expand the observations into an observation data vector, and construct the observation standard deviation vector; S2: Based on the expected detection depth, establish a non-uniform one-dimensional grid with shallow densification and deep sparseness. Use the common logarithm of resistivity of each model unit as the inversion parameter to construct first-order and second-order depth difference matrices and reference model constraints. S3: Set up multiple uniform initial resistivity models, select the initial model based on the fitting error of normalized data, and perform pilot scans on multiple smoothing regularization parameters. S4: Calculate the Jacobian matrix of the predicted data after normalization of the observation standard deviation at the leading Occam's calculus model, and adjust the first-order and second-order depth difference constraints and the reference model constraints respectively. S5: Scan the smoothing regularization parameter again to obtain the sensitivity adaptive smoothing model; update the iterative reweighting coefficients according to the first difference of the model, and gradually reduce the smoothing stability parameter to carry out iterative reweighted least squares block optimization; when the normalized fitting error of the block optimization model exceeds the allowable value, reduce the block regularization parameter and solve it again to obtain the block multi-unit resistivity model. S6: Calculate the approximate model resolution, combine the relative sensitivity of the model elements and the resolution diagonal value to determine the reliable interpretation depth index, and output the blocky multi-element resistivity model, the relative sensitivity of the model elements, the resolution diagonal value, the data fitting and the normalized residual.

2. The small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization according to claim 1, characterized in that, The specific method of S1 is as follows: a small loop transmitter and a small loop receiver are arranged in the detection area. The transmitting loop is sequentially supplied with periodic current at multiple transmission frequencies. The receiving loop obtains the secondary magnetic field response generated by the underground medium at multiple transmit-receive offset positions. The received signal is decomposed into real part response and imaginary part response with the phase of the transmitting current as a reference. Let the set of transmission frequencies be . In the formula, For the first One transmission frequency, This represents the total number of transmission frequencies. Let the set of transmit and receive offsets be: In the formula, Let q be the transmit / receive offset. This represents the total transmit and receive offset. When the emission ray is the radius When the small circular loop appears, the first... The equivalent magnetic dipole moment at each transmission frequency is: In the formula, For the equivalent magnetic dipole moment, For the number of turns of the transmit loop, For the first... One transmission frequency and the Transmit / receive offset Let the complex response of the second magnetic field in the sensitive direction of the receiving loop be... The direct first-order magnetic field complex response is ; When expressed in parts per million, the real and imaginary observations are respectively: , ; Following the order of real parts first, then imaginary parts within each frequency-offset combination, all observations are expanded as follows: The total number of individual observations is: ;set up for The first in One single observation, of which ; Set the real part fixed absolute error base respectively and the imaginary part with a fixed absolute error base The standard deviations of the real and imaginary observations corresponding to the i-th transmit frequency and the q-th transmit / receive offset are respectively: , In the formula, The relative error coefficient; Arrange the standard deviations in the same order according to the observed data vectors. And construct the observation error weighting matrix: .

3. The small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization according to claim 2, characterized in that, The standard deviation is used to normalize the data residuals, normalize the data Jacobian, evaluate the model fit, and calculate the model resolution. It does not require adding random disturbances to the observed data, nor does it change the original observed values. The fixed absolute error floor is determined based on the instrument background, field testing, repeated observation errors, and the field noise level.

4. The small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization according to claim 1, characterized in that, The specific method of S2 is that the non-uniform one-dimensional grid adopts a smaller model unit thickness near the surface, and the model unit thickness is gradually increased with the depth. A uniform half-space is set at the bottom, and the number of model elements with finite thickness is [number missing]. If a half-space is set at the bottom, the total number of model parameters is: ;No. The resistivity of each model element is denoted as The inversion parameters are defined as follows: The model vector is: ; A first-order fundamental difference matrix is ​​established based on the center-to-center distance of the model units and the depth-based fundamental weights. Second-order fundamental difference matrix and the basic reference weight matrix The first-order difference is used to suppress non-physical jumps between adjacent model units, the second-order difference is used to suppress high-frequency oscillations formed by three consecutive model units, and the reference model constraint is used to limit the low-sensitivity parameters from deviating excessively from the initial background.

5. The small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization according to claim 1, characterized in that, The specific method of S3 is as follows: setting Uniform initial resistivity ,in For each uniform model, the predicted data are calculated using a one-dimensional small loop frequency-domain electromagnetic forward modeling operator. And calculate the normalized root mean square error: ; The uniform model with the smallest normalized root mean square error was selected as the initial model. and reference model For candidate smoothing regularization parameters The joint residual vector for smooth inversion is constructed as follows: It is formed by sequentially concatenating the normalized data residuals, first-order difference residuals, second-order difference residuals, and reference model constraint residuals. After expansion... ,in, The first in the joint residual vector Each scalar residual component , The total number of elements in the joint residual vector is given by the formula. The second-order difference relative weights are used for the smoothing phase; These are the constraint coefficients for the reference model during the smoothing phase. A nonlinear least squares method with parameterized upper and lower bounds is used to solve the problem, and a soft L1 robust loss is employed to reduce the dominance of a few large residuals on model updates: the soft L1 robust loss introduces a degree parameter. If the value is greater than 0, the smoothing inversion objective function is uniformly written as: In the formula, Let be the total number of elements in the joint residual vector, where The scale controlling the transition of residuals from approximately quadratic penalty to approximately linear penalty; when the absolute value of the residuals is much smaller than... When the residual value is close to the quadratic loss, the dominance of a single large residual on the model update is reduced when the absolute value of the residual is large. Sensitivity adaptive rescan and IRLS trial use the same soft L1 scale parameter, scanning multiple times sequentially. Candidate values; when multiple candidate models satisfy When, choose Larger candidate models are used as leading Occam models Prioritize obtaining the smoothest model that meets the data fitting requirements; otherwise, select the candidate model with the smallest normalized root mean square error; if no candidate model meets the target error, select the candidate model with the smallest normalized root mean square error.

6. The small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization according to claim 1, characterized in that, The specific method of S4 is as follows: In the leading Occam model Calculate the normalized data Jacobian matrix , its first line, number The column represents the rate of change of the k-th observation in the observation data vector relative to the j-th logarithmic resistivity parameter: ; When using finite difference: In the formula, For the first Finite difference perturbation of each model parameter; For the first A unit direction vector with one element being 1 and the rest being 0. If the positive perturbation exceeds the upper limit of the parameters, then the negative perturbation is used. No. The original sensitivity of each model unit is: A moving average operator along the depth direction is used. Smooth the original sensitivity and normalize the maximum value: ; Set the relative sensitivity for truncation: ,in, As the lower limit of relative sensitivity, the first-order difference constraint is adjusted according to the truncated relative sensitivity of adjacent model units, the second-order difference constraint is adjusted according to the truncated relative sensitivity of three consecutive model units, and the reference model constraint is adjusted according to the truncated relative sensitivity of a single model unit. The first-order difference constraint adjustment factor is: The second-order difference constraint adjustment factor is: The reference model constraint adjustment factor is: In the formula, To adjust the intensity of differential constraint sensitivity, The sensitivity adjustment intensity is constrained by the reference model; This yields the first-order difference matrix, the second-order difference matrix, and the reference model weight matrix for sensitivity adaptation. , , The lower the relative sensitivity, the stronger the corresponding constraint; the higher the relative sensitivity, the closer the corresponding constraint is to the basic constraint. use , and The regularization parameter is scanned again near the leading Occam's parameter, and the sensitivity adaptive smoothing model is obtained according to the selection rule of S3. and its regularization parameters .

7. The small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization according to claim 1, characterized in that, The specific method of S5 is as follows: by As the initial model for block optimization, the first In the next reweighted iteration, the first-order difference gradient is calculated: , No. The reweighting coefficients for each first-order difference term are: In the formula, For the first The smoothing and stabilizing parameters of the next iteration; and These are the lower and upper limits of the reweighted coefficients, respectively. The smoothing stability parameter decreases successively according to the following formula: The joint residual vector for block optimization is: In the formula, For block regularization parameters; The relative weights are the second-order difference of the block stage; The constraint coefficients of the block stage reference model are given in the formula. The decreasing coefficient, To provide a lower limit for smooth and stable parameters; Let the normalized root mean square error of the sensitivity adaptive smoothing model be... The target fitting error is The allowable error for the block model is: ,in, The fitting error protection coefficient is set when the normalized root mean square error obtained in a trial calculation is greater than 1. At that time, Multiply by an attenuation coefficient less than 1 The solution was then recalculated, and after a preset number of iterations, a block-shaped multi-element model was obtained. .

8. The small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization according to claim 1, characterized in that, The specific method of S6 is as follows: Normalized data Jacobian matrix at the leading Occam's model And the approximate model resolution is calculated using the regularization matrix of the selected sensitivity adaptive smoothing model. The regularization matrix is: The approximate model resolution matrix is: In the formula, For numerical stability damping, for The first-order identity matrix, the second-order identity matrix The resolution diagonal value of each model parameter is: Let the first The bottom boundary depth of a finite thickness model element is The relative sensitivity threshold is The resolution diagonal threshold is The reliable explanation depth metric is defined as: If there is no finite thickness model element that simultaneously satisfies the threshold, then Take zero.

9. The small-loop frequency-domain electromagnetic inversion method based on sensitivity adaptive regularization according to claim 8, characterized in that, The blocky multi-unit resistivity model is used to identify underground electrical variations and low-resistivity anomaly ranges. The reliable interpretation depth index is only used for result diagnosis and interpretation marking, and is not used as a pre-given true layer thickness, true interface, or deep target location.

Images