A method for simultaneous inversion of dual parameters of three-dimensional frequency-domain ground penetrating radar

Abstract

The present invention discloses a three-dimensional frequency domain ground penetrating radar dual-parameter synchronous inversion method, including inputting observation data, setting a three-dimensional initial model, and initializing parameters; forward modeling the established three-dimensional initial model; constructing an inversion objective function and decomposing it into two alternating sub-problems; initializing dielectric constant and conductivity to obtain an initial pilot model vector; solving sub-problem 1 according to the pilot model vector and updating the model medium parameter vector; solving sub-problem 2 according to the obtained model medium parameter vector and updating the pilot model vector; synchronously and alternately updating until the inversion termination condition is met and the inversion result is output. The present invention can synchronously update dielectric constant and conductivity values ​​to improve the efficiency of inversion; improve the total variation model constraint method to the full waveform inversion process, increase the robustness of the algorithm, combine the frequency weight factor, and use a parallel method to achieve 3D GPR quantitative imaging, so that the imaging process is more efficient and accurate.

Description

Technical Field

[0001] The invention belongs to the field of geophysical detection, and in particular relates to a three-dimensional frequency domain ground penetrating radar dual-parameter synchronous inversion method. Background Art

[0002] Ground penetrating radar (GPR) is an electromagnetic nondestructive detection technology for locating underground anomalies, structures and characteristics or invisible targets inside objects. In actual GPR detection, the position, size and dimensions of the detection target can be roughly inferred by relying solely on the acquired radar profile. However, bottom penetrating radar is mostly limited to rough estimation, and cannot provide parameters such as the dielectric constant and conductivity of the anomaly to achieve accurate definition of the target attributes. Full waveform inversion (FWI) is an important emerging imaging method. GPR full waveform inversion can obtain physical parameters such as dielectric constant and conductivity of the medium, which is helpful for accurate description of lithology and prediction of anomalies, and is a powerful tool for reconstructing complex geological structures near the surface.

[0003] At present, FWI has received extensive attention and research in GPR imaging, and can be divided into time domain FWI and frequency domain FWI. In the frequency domain, by inverting the data from low frequency to high frequency, the nonlinearity of the problem can be effectively reduced during the calculation. By reasonably selecting a few frequency component data for inversion, the same effect as the time domain inversion can be achieved.

[0004] Given that the 3D scattering of electromagnetic transverse waves is more complex than the scattering of 2D longitudinal waves, the nonlinearity between the simulated data and the model parameters is stronger, and there are more underground parameters in the 3D model to fit the observed data (more multi-solution). In addition, solving the 3D inversion problem requires greater storage costs, and the multiple calculation costs of the forward modeling problem involved are much higher than 2D. Therefore, most of the existing technologies are about the research of 2D ground penetrating radar data FWI. However, for objects with 3D changes, the 2D imaging method limits the applicability and accuracy of FWI.

[0005] Considering that the full waveform inversion of ground penetrating radar is a highly nonlinear and ill-posed problem. Lavoué et al. introduced Tikhonov regularization in the 2D frequency domain GPR full waveform inversion to improve the high amplitude oscillation phenomenon in the conductivity model. However, the inversion of Tikhonov regularization tends to produce a smooth model. Watson uses an approximate total variation regularization to constrain the dielectric constant model, so that its inversion results have a sharp interface. This approximate total variation method may be unstable and the convergence of the inversion cannot be guaranteed. Therefore, it is necessary to find new methods for full waveform inversion of GPR data. Summary of the invention

[0006] The object of the present invention is to provide a three-dimensional frequency domain ground penetrating radar dual-parameter synchronous inversion method, which can improve the inversion accuracy of the ground penetrating radar.

[0007] The three-dimensional frequency domain ground penetrating radar dual-parameter synchronous inversion method provided by the present invention comprises the following steps:

[0008] S1. Input observation data, set up the 3D initial model, and initialize parameters;

[0009] S2. Forward modeling of the established three-dimensional initial model;

[0010] S3. construct an inversion objective function and decompose it into two alternating sub-problems, including sub-problem 1 and sub-problem 2, sub-problem 1 is to solve the model medium parameter vector, and sub-problem 2 is to solve the pilot model vector;

[0011] S4. Initialize the dielectric constant and conductivity to obtain an initial pilot model vector;

[0012] S5. Solve subproblem 1 according to the pilot model vector and update the model medium parameter vector;

[0013] S6. Solve subproblem 2 according to the obtained model medium parameter vector and update the leading model vector;

[0014] S7. Repeat steps S5-S6, and update synchronously and alternately until the inversion termination condition is met, then end the inversion and output the inversion result.

[0015] The step S2 includes forward modeling the established three-dimensional initial model using the finite element method, specifically using the linear vector basis function finite element method to discretize the simulation area, dividing the grid into rectangular block units, and using anisotropic perfectly matched layers at the truncation boundary; the frequency domain Maxwell equations satisfied by the three-dimensional GPR are expressed as:

[0016]

[0017] in, represents the electric field strength; ε represents the current dielectric constant; μ represents the magnetic permeability; σ represents the current conductivity; ω represents the angular frequency, ω=2πf, f represents the frequency; j represents the imaginary unit, is the current density.

[0018] Find the linear equations for a multi-source system:

[0019]

[0020] in, represents the coefficient matrix common to all frequency-dependent sources; represents the discrete electric field value; represents the discretized excitation source term; ω represents the angular frequency; ε represents the dielectric constant; σ represents the conductivity; s represents the excitation source at the current position; the coefficient matrix common to all sources related to frequency By all unit matrices synthesis:

[0021]

[0022] in, Represents the volume integral within the current unit e; and is a vector basis function; when the calculation area is discretized using hexahedral units: l,n=1,2,…,12; for the matrix Perform an LU decomposition; the forward operator for a single source and a single frequency is described as:

[0023]

[0024] in, Represents the simulation data under the angular frequency ω and the current excitation source s; The corresponding position matrix of the receiving point representing the current source; represents the discrete electric field value; represents the coefficient matrix common to all frequency-dependent sources; represents the discretized excitation source term.

[0025] Step S3, using known measured data to reconstruct the model parameters in the underground medium, including the spatial distribution of the dielectric constant ε and the conductivity σ; and introducing the total variation regularization method to invert the objective function Defined as:

[0026]

[0027] in, represents the model medium parameter vector, Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the prior model vector, Represents the prior model vector corresponding to the current dielectric constant, Represents the prior model vector corresponding to the current conductivity; represents the data objective function; ε represents the current dielectric constant; σ represents the current conductivity; α ε , α σ , β ε and β σ represents a regularization parameter with a positive value; ||·|| TV represents the total variation operator; the inversion objective function It can be further decomposed into subproblem 1 and subproblem 2. Subproblem 1 is to solve the model medium parameter vector, and subproblem 2 is to solve the prior model vector.

[0028] The sub-question 1 is specifically:

[0029]

[0030] in, represents the model medium parameter vector, Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the model medium parameter vector after the kth iteration; represents the inversion objective function with respect to the model medium parameter vector; represents the data objective function; ε represents the current dielectric constant; σ represents the current conductivity; and represents the regularization parameter with a positive value at the kth iteration; represents the leading model vector corresponding to the current dielectric constant after the k-1th iteration, represents the leading model vector corresponding to the current conductivity after the k-1th iteration; where the regularization parameter of the kth iteration and They are:

[0031] and

[0032] in, and is a dimensionless number; represents the data objective function; Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the leading model vector corresponding to the current dielectric constant after the k-1th iteration, Represents the leader model vector corresponding to the current conductivity after the k-1th iteration.

[0033] The sub-question 2 includes:

[0034]

[0035] Among them, ε represents the current dielectric constant; σ represents the current conductivity; β σ and β ε represents a regularization parameter with a positive value; ||·|| TVrepresents the total variation operator; represents the leading model vector corresponding to the current dielectric constant after the kth iteration; Represents the leading model vector corresponding to the current dielectric constant; represents the objective function of the leading model vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current dielectric constant after the kth iteration; represents the leading model vector corresponding to the current dielectric constant after the kth iteration; Represents the leading model vector corresponding to the current dielectric constant; represents the objective function of the leader model vector corresponding to the current conductivity; represents the model medium parameter vector corresponding to the current dielectric constant after the kth iteration; among them, the regularization parameter of the kth iteration and They are:

[0036] and in, and is a dimensionless number; represents the data objective function; Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the leading model vector corresponding to the current dielectric constant after the k-1th iteration, Represents the leader model vector corresponding to the current conductivity after the k-1th iteration.

[0037] The step S5 specifically includes: solving sub-problem 1 using the limited memory BFGS method, and for the model medium parameter vector after the k+1th iteration

[0038]

[0039] in, represents the model medium parameter vector after the kth iteration; a (k ) is the search step size of the kth iteration; represents the search direction of the kth iteration, represents an approximation of the inverse Hessian matrix; represents the objective function of the model medium parameter vector for the kth iteration The gradient of

[0040] The objective function of the model medium parameter vector for the kth iteration Gradient Specifically:

[0041]

[0042] in, represents the gradient of the objective function with respect to the model medium parameter vector; represents the data objective function; and Indicates that the value of the k-th iteration is positive regularization parameter; Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; Represents the prior model vector corresponding to the current dielectric constant after the k-1th iteration; Represents the prior model vector corresponding to the current conductivity after the k-1th iteration; Represents the gradient of the data objective function with respect to the model medium parameter vector corresponding to the current dielectric constant; represents the gradient of the data objective function with respect to the model medium parameter vector corresponding to the current conductivity;

[0043] For the underground medium unit e, the gradient of the data objective function of the model medium parameter vector corresponding to the current dielectric constant Adjoint state method is used to obtain:

[0044]

[0045] in, is the model medium parameter vector; N ω Indicates the maximum number of angular frequencies; N s Indicates the maximum number of excitation sources; i indicates the frequency number; r indicates the excitation source number; is the difference between the observed data and the simulated data, specifically represents observation data; represents simulated data; T represents transpose, and * represents conjugation; represents the real part operator; represents discrete electric field values; represents the impedance matrix; The corresponding position matrix of the receiving point representing the current source; is a sparse matrix, Impedance matrix The model medium parameter m corresponding to the current unit e e sensitivity, is a diagonal matrix, the coefficient corresponding to the current unit e is a non-zero element, through Calculated; its model medium parameter vector m e Including dielectric constant and conductivity;

[0046] Solving the cell matrix For the model medium parameter corresponding to the current dielectric constant ε in the corresponding current unit e Sensitivity Solving the cell matrix For the model medium parameter corresponding to the current conductivity σ in the corresponding current unit e Sensitivity

[0047] and

[0048] in, represents the model medium parameter corresponding to the current dielectric constant ε in the current unit e; ω represents the angular frequency; j represents the imaginary unit, ε0 is the dielectric constant of free space; For vector basis functions, when the computational domain is discretized using hexahedral units: l,n=1,2,…,12; Represents the volume integral within the current unit e; represents the model medium parameter corresponding to the current conductivity σ in the current unit e; β represents a regularization parameter with a positive value; σ0 is the model background conductivity.

[0049] The step S6 specifically includes solving subproblem 2 using the Split-Bregman iterative method and updating:

[0050] Subproblem 2 is expressed using the unified TV model:

[0051]

[0052] in, represents the prior model vector, Represents the prior model vector corresponding to the current dielectric constant, Represents the prior model vector corresponding to the current conductivity; represents the model medium parameter vector after the kth iteration; Represents the model medium parameter vector corresponding to the current dielectric constant after the kth iteration; represents the model medium parameter vector corresponding to the current conductivity after the kth iteration; ||·|| TV represents the total variation operator; μ TV =μ TV,ε ,μ TV,σ , while μ TV , μ TV,ε and μ TV,σ represents the regularization parameter;

[0053] In the total variation operator ||·|| TV In the example, for a 3D isotropic regularization operator It is expressed as:

[0054]

[0055] Among them, e represents the current unit; Represents the gradient of the prior model vector on the x-axis; Represents the gradient of the prior model vector on the y-axis; Represents the gradient of the prior model vector on the z-axis; sets auxiliary variables and The L1 constraint and L2 penalty components of the segmentation problem give the optimality condition of the Bregman formula for the unified TV model splitting:

[0056]

[0057] Where λ = 2μ TV ;μ TV =μ TV,ε ,μ TV,σ , while μ TV , μ TV,ε and μ TV,σ represents the regularization parameter; represents the identity matrix; represents the model medium parameter vector after the kth iteration; Indicates the three directions of x, y, and z in three dimensions; represents the kth iteration, Auxiliary variables in direction; represents the kth iteration, The Bregman coefficient in the direction, its iterative formula is represents the k+1th iteration, Bregman coefficient in direction; represents the k+1th iteration, Auxiliary variables in direction; express and The vector difference is The transposed matrix of the gradient in the direction; Δ is the Laplace operator; represents the prior model vector after the k+1th iteration;

[0058] set up Auxiliary variables Solve explicitly using the generalized contraction formula:

[0059]

[0060] Where λ = 2μ TV ;μ TV =μ TV,ε ,μ TV,σ , while μ TV , μ TV,ε and μ TV,σ represents the regularization parameter; represents the shrinkage coefficient, Represents the prior model vector after the kth iteration Directional gradient; Represents the gradient of the prior model vector in the x direction after the kth iteration; Represents the gradient of the prior model vector in the y direction after the kth iteration; Represents the gradient of the prior model vector in the z direction after the kth iteration; represents the kth iteration, Bregman coefficient in direction; represents the Bregman coefficient in the x direction of the kth iteration; represents the Bregman coefficient in the y direction of the kth iteration; represents the Bregman coefficient in the z direction at the kth iteration.

[0061] The step S7 is specifically to repeat steps S5-S6, update subproblem 1 and subproblem 2 synchronously and alternately, until the set inversion termination condition is met, end the inversion and output the inversion result: the set inversion termination condition includes reaching the set iteration accuracy or reaching the maximum number of iterations; the modified total variation problem alternately solves subproblem 1 and subproblem 2; in the kth iteration, according to step S4, use the variable and Solve the model medium parameter vector after the kth iteration and represents the leading model vector corresponding to the current dielectric constant after the k-1th iteration; represents the pilot model vector corresponding to the current conductivity after the k-1th iteration; then according to step S5, the model medium parameter vector after the kth iteration is used and Solving for variables and represents the leading model vector corresponding to the current dielectric constant after the kth iteration; represents the leading model vector corresponding to the current conductivity after the kth iteration; the obtained variable and Used in the next iteration step; for the first iteration step, the pilot model vector corresponding to the initial current dielectric constant Equal to the model medium parameter vector corresponding to the initial current dielectric constant The leading model vector corresponding to the initial current conductivity Equal to the model medium parameter vector corresponding to the initial current conductivity Get the inversion sequence:

[0062]

[0063]

[0064] The three-dimensional frequency domain ground penetrating radar dual-parameter synchronous inversion method provided by the present invention adopts a finite element solver to simulate a GPR forward modeling problem; uses the L-BFGS method in the quasi-Newton method to solve the inversion problem, thereby avoiding the calculation of the Hessian matrix; through parameter transformation and parameter adjustment factors, it is convenient to synchronously update the dielectric constant and the conductivity value within the same descent step to improve the inversion efficiency; introduces an improved total variation model constraint method into the GPR full waveform inversion process to increase the robustness of the algorithm, combines with the frequency weight factor, and adopts a parallel method to realize the quantitative imaging of 3D GPR data, thereby making the imaging process more efficient and accurate. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0066] Figure 2 a and Figure 2 b is a schematic diagram of a three-dimensional model of the model 1 according to an embodiment of the present invention.

[0067] Figure 3 a and Figure 3 b is a schematic diagram of the full waveform inversion results of model 1 of an embodiment of the present invention.

[0068] Figure 4 a and Figure 4 b is a schematic diagram of relative data residual and reconstruction error of model 1 of an embodiment of the present invention.

[0069] Figure 5 a and Figure 5 b is a schematic diagram of a three-dimensional model of model 2 of an embodiment of the present invention.

[0070] Figure 6 a. Figure 6 b. Figure 6 c and Figure 6 d is a schematic diagram of the full waveform inversion results under different regularization strategies of Model 2 of an embodiment of the present invention. DETAILED DESCRIPTION

[0071] like Figure 1The flowchart of the method of the present invention is as follows: The three-dimensional frequency domain ground penetrating radar dual-parameter synchronous inversion method provided by the present invention comprises the following steps:

[0072] S1. Input observation data, set up the 3D initial model, and initialize parameters;

[0073] S2. Forward modeling of the established three-dimensional initial model;

[0074] S3. construct an inversion objective function and decompose it into two alternating sub-problems, including sub-problem 1 and sub-problem 2, sub-problem 1 is to solve the model medium parameter vector, and sub-problem 2 is to solve the pilot model vector;

[0075] S4. Initialize the dielectric constant and conductivity to obtain an initial pilot model vector;

[0076] S5. Solve subproblem 1 according to the pilot model vector and update the model medium parameter vector;

[0077] S6. Solve subproblem 2 according to the obtained model medium parameter vector and update the leading model vector;

[0078] S7. Repeat steps S5-S6, and update synchronously and alternately until the inversion termination condition is met, then end the inversion and output the inversion result.

[0079] The step S2 includes forward modeling the established three-dimensional initial model using the finite element method, specifically using the linear vector basis function finite element method to discretize the simulation area, dividing the grid into rectangular block units, and using anisotropic perfectly matched layers at the truncation boundary: The frequency domain Maxwell equations satisfied by the three-dimensional GPR are expressed as:

[0080]

[0081] in, Indicates the electric field strength; represents the magnetic field strength; ε represents the current dielectric constant; μ represents the magnetic permeability; σ represents the current conductivity; ω represents the angular frequency, ω=2πf, f represents the frequency; j represents the imaginary unit, represents the current density; the time harmonic factor of this equation is e jωt ;

[0082] The finite element method of linear vector basis function is used for discretization, and anisotropic perfectly matched layers are used at the truncation boundary. After overall merging, the linear equations of the multi-source system are obtained:

[0083]

[0084] in, represents the coefficient matrix common to all frequency-dependent sources; represents the discrete electric field value; represents the discretized excitation source term; ω represents the angular frequency; ε represents the dielectric constant; σ represents the conductivity; s represents the excitation source at the current position; the coefficient matrix common to all sources related to frequency By all unit matrices synthesis:

[0085]

[0086] in, Represents the volume integral within the current unit e; and is a vector basis function; when the computational domain is discretized using hexahedral units: l,n=1,2,…,12; here, rectangular block units are used, which are easy to implement. The linear system is solved using the direct solver MUMPS. For a given medium and frequency, only the matrix Perform an LU decomposition to effectively solve multiple right-hand side problems.

[0087] The forward operator for a single source and a single frequency can be described as:

[0088]

[0089] in, Represents the simulation data under the angular frequency ω and the current excitation source s; The corresponding position matrix of the receiving point representing the current source; represents the discrete electric field value; represents the coefficient matrix common to all frequency-dependent sources; represents the discretized excitation source term.

[0090] The step S3 specifically includes that the GPR full waveform inversion is essentially to reconstruct the model parameters in the underground medium using known measured data, including the spatial distribution of the dielectric constant ε and the conductivity σ; according to the optimal fit between the forward simulation data and the measured data, the data objective function Defined as:

[0091]

[0092] in, is the model medium parameter vector; N ω Indicates the number of angular frequencies; N s represents the number of excitation sources; η i Represents the weighting factor of multi-frequency data; ω i represents the i-th frequency; represents the jth excitation source; is the difference between the observed data and the simulated data, specifically represents observation data; represents simulated data; T represents transposition, * represents conjugation; η i Represents the weighting factor for multi-frequency data:

[0093]

[0094] Since the dielectric constant and conductivity differ greatly in magnitude, the inversion calculation is inconvenient. In order to update the dielectric constant and conductivity model at the same time in each iteration, it is necessary to consider how to balance the influence of different parameters. After the change, the model medium parameter vector for:

[0095]

[0096] in, Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; ε r represents the relative dielectric constant, ε r =ε / ε0, ε0 is the free space dielectric constant, ε is the current dielectric constant; σ r Relative conductivity, σ r =σ / σ0, σ0 is the model background conductivity, the model background conductivity is the conductivity of the set reference medium, σ is the current conductivity; к represents the parameter adjustment factor, which is used to control the relative conductivity σ in the optimization process r Relative dielectric constant ε r The weights of are used to avoid the uncertainty in the inversion process caused by inaccurate definitions of relative conductivity and relative permittivity.

[0097] Since the FEI of GPR is an ill-posed problem, in order to improve the robustness of FWI, an improved total variation regularization method is often introduced in the inversion process. Defined as:

[0098]

[0099] in, represents the model medium parameter vector; Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the prior model vector, Represents the prior model vector corresponding to the current dielectric constant, Represents the prior model vector corresponding to the current conductivity; represents the data objective function; represents the model medium parameter vector corresponding to the current dielectric constant; ε represents the current dielectric constant; σ represents the current conductivity; α ε , α σ , β ε and β σ represents a regularization parameter with a positive value; ||·|| TV represents the total variation operator; the inversion objective function It is further decomposed into two alternating minimization subproblems:

[0100] Sub-question 1:

[0101]

[0102] Sub-question 2:

[0103]

[0104] in, represents the model medium parameter vector; Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the model medium parameter vector after the kth iteration; represents the objective function with respect to the model medium parameter vector; represents the data objective function; ε represents the current dielectric constant; σ represents the current conductivity; α ε , α σ , α ε,σ , β σ and β ε represents a regularization parameter with a positive value; represents the prior model vector corresponding to the current dielectric constant after the k-1th iteration, represents the prior model vector corresponding to the current conductivity after the k-1th iteration; ||·|| TV represents the total variation operator; represents the leading model vector corresponding to the current dielectric constant after the kth iteration; Represents the leading model vector corresponding to the current dielectric constant; represents the objective function of the leading model vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current dielectric constant after the kth iteration; represents the leading model vector corresponding to the current dielectric constant after the kth iteration; Represents the leading model vector corresponding to the current dielectric constant; represents the objective function of the leader model vector corresponding to the current conductivity; Represents the model medium parameter vector corresponding to the current dielectric constant after the kth iteration.

[0105] Subproblem 1 is to use a prior model based on Tikhonov regularization and multiple parameters, that is, the prior model vector corresponding to the current dielectric constant and the current conductivity after the k-1th iteration Solve the model medium parameter vector after the kth iteration For the full waveform inversion problem, Tikhonov regularized inversion can ensure the stability of the inversion;

[0106] Subproblem 2 is to perform two standard L2-total variation minimization methods to solve the model medium parameter vector corresponding to the current dielectric constant after the kth iteration with different parameters: The prior model vector corresponding to the current conductivity To maintain the inversion result, the inversion result is the model medium parameter vector corresponding to the current dielectric constant after the kth iteration The model medium parameter vector corresponding to the current conductivity The interface of the waveform is clear; solving the equations by alternating iteration can not only enhance the clarity of the interface, but also ensure that the full waveform inversion is more stable.

[0107] Among them, the regularization parameter of the kth iteration is and They are:

[0108] and

[0109] in, and is a dimensionless number. In this embodiment, represents the data objective function; Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the leading model vector corresponding to the current dielectric constant after the k-1th iteration, represents the leading model vector corresponding to the current conductivity after the k-1th iteration; the regularization parameter of the kth iteration and The iterative equation can realize the adaptive selection of regularization factors.

[0110] The step S5 specifically includes: solving sub-problem 1 by using the limited memory BFGS method, and synchronously updating the dielectric constant model and the conductivity model;

[0111] The L-BFGS algorithm is used to solve the model medium parameter vector The inverse problem is that for the model medium parameter vector after the k+1th iteration

[0112]

[0113] in, represents the model medium parameter vector after the kth iteration; a (k) is the search step size for the kth iteration, calculated using the inexact line search method with the strong Wolfe criterion; represents the search direction of the kth iteration, It represents an approximation of the inverse Hessian matrix, which can be calculated using two loop recursive algorithms; represents the objective function of the model medium parameter vector for the kth iteration The gradient is:

[0114]

[0115] in, represents the gradient of the objective function with respect to the model medium parameter vector; represents the data objective function; and Indicates that the value of the k-th iteration is positive regularization parameter; Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the leading model vector corresponding to the current dielectric constant after the k-1th iteration; represents the leading model vector corresponding to the current conductivity after the k-1th iteration; Represents the gradient of the data objective function with respect to the model medium parameter vector corresponding to the current dielectric constant; Represents the gradient of the data objective function with respect to the model medium parameter vector corresponding to the current conductivity.

[0116] For the underground medium unit e, the gradient of the data objective function of the model medium parameter vector corresponding to the current dielectric constant Adjoint state method is used to obtain:

[0117]

[0118] in, is the model medium parameter vector; N ω Indicates the maximum number of angular frequencies; N s Indicates the maximum number of excitation sources; i indicates the frequency number; r indicates the excitation source number; is the difference between the observed data and the simulated data, specifically represents observation data; represents simulated data; T represents transpose, and * represents conjugation; represents the real part operator; represents discrete electric field values; represents the impedance matrix; The corresponding position matrix of the receiving point representing the current source; is a sparse matrix, Impedance matrix The model medium parameter m corresponding to the current unit e e sensitivity, is a diagonal matrix, and only the coefficient corresponding to the current unit e is a non-zero element, which can be obtained by The model medium parameter vector m is calculated. e Includes dielectric constant and conductivity.

[0119] Solving the cell matrix For the model medium parameter corresponding to the current dielectric constant ε in the corresponding current unit e Sensitivity Solving the cell matrix For the model medium parameter corresponding to the current conductivity σ in the corresponding current unit e Sensitivity

[0120]

[0121] in, represents the model medium parameter corresponding to the current dielectric constant ε in the current unit e; ω represents the angular frequency; j represents the imaginary unit, ε0 is the dielectric constant of free space; For vector basis functions, when the computational domain is discretized using hexahedral units: l,n=1,2,…,12; Represents the volume integral within the current unit e; represents the model medium parameter corresponding to the current conductivity σ in the current unit e; β represents a regularization parameter with a positive value; σ0 is the model background conductivity.

[0122] The step S6 is specifically to solve subproblem 2, that is, the auxiliary variables of the two parameters, using the Split-Bregman iterative method and update them, which can be done by the following steps:

[0123] Subproblem 2 is represented by the following unified TV model:

[0124]

[0125] in, represents the prior model vector, Represents the prior model vector corresponding to the current dielectric constant, Represents the prior model vector corresponding to the current conductivity; represents the model medium parameter vector after the kth iteration; Represents the model medium parameter vector corresponding to the current dielectric constant after the kth iteration; represents the model medium parameter vector corresponding to the current conductivity after the kth iteration; ||·|| TV represents the total variation operator; μ TV =μ TV,ε ,μ TV,σ , μ TV,ε and μ TV,σ represents the regularization parameter;

[0126] For a 3D isotropic regularization operator It is expressed as:

[0127]

[0128] Among them, e represents the current unit; Represents the gradient of the prior model vector on the x-axis; Represents the gradient of the prior model vector on the y-axis; Represents the gradient of the prior model vector on the z-axis; Subproblem 2 cannot be solved directly Therefore, you need to set the auxiliary variable and The L1 constraint and L2 penalty components of the segmentation problem give the optimality condition of the Bregman formula for the unified TV model splitting:

[0129]

[0130] Where λ = 2μ TV ;μ TV =μ TV,ε ,μ TV,σ , while μ TV,ε and μ TV,σ represents the regularization parameter; represents the identity matrix; represents the model medium parameter vector after the kth iteration; Indicates the three directions of x, y, and z in three dimensions; represents the kth iteration, Auxiliary variables in direction; represents the kth iteration, The Bregman coefficient in the direction, its iterative formula is represents the k+1th iteration, Bregman coefficient in direction; represents the k+1th iteration, Auxiliary variables in direction; express and The vector difference is The transposed matrix of the gradient in the direction; Δ is the Laplace operator; represents the prior model vector after the k+1th iteration;

[0131] set up Auxiliary variables Solve explicitly using the generalized contraction formula:

[0132]

[0133] Where λ = 2μ TV ;μ TV =μ TV,ε ,μ TV,σ , while μ TV , μ TV,ε and μ TV,σ represents the regularization parameter; represents the shrinkage coefficient, Represents the prior model vector after the kth iteration That is, the gradient in the x, y or z direction; Represents the gradient of the prior model vector in the x direction after the kth iteration; Represents the gradient of the prior model vector in the y direction after the kth iteration; Represents the gradient of the prior model vector in the z direction after the kth iteration; represents the kth iteration, Bregman coefficient in direction; represents the Bregman coefficient in the x direction of the kth iteration; represents the Bregman coefficient in the y direction of the kth iteration; represents the Bregman coefficient in the z direction for the kth iteration.

[0134] The step S7 is specifically to repeat steps S5-S6, and synchronously and alternately update subproblem 1 and subproblem 2 until the set inversion termination condition is met, and then the inversion is terminated and the inversion result is output. The inversion termination condition is set to include satisfying the iteration accuracy or reaching the maximum number of iterations. The modified total variation problem alternately solves subproblem 1 and subproblem 2; in the kth iteration, according to step S4, the variable and Solve the model medium parameter vector after the kth iteration and represents the leading model vector corresponding to the current dielectric constant after the k-1th iteration; represents the pilot model vector corresponding to the current conductivity after the k-1th iteration; then according to step S5, the model medium parameter vector after the kth iteration is used and Solving for variables and represents the leading model vector corresponding to the current dielectric constant after the kth iteration; represents the leading model vector corresponding to the current conductivity after the kth iteration; the obtained variable and Used in the next iteration step; for the first iteration step, the pilot model vector corresponding to the initial current dielectric constant Equal to the model medium parameter vector corresponding to the initial current dielectric constant The leading model vector corresponding to the initial current conductivity Equal to the model medium parameter vector corresponding to the initial current conductivity Get the inversion sequence:

[0135]

[0136]

[0137] The method of the present invention is further described below in conjunction with an implementation case:

[0138] like Figure 2 a and Figure 2 b is a schematic diagram of a three-dimensional model of the model 1 according to an embodiment of the present invention. Figure 2 a represents the dielectric constant ε of model 1 of the embodiment of the present invention r Schematic diagram of the three-dimensional model; Figure 2 b is a schematic diagram of a three-dimensional model of the conductivity of the model 1 of the embodiment of the present invention. There is a 0.48m thick air layer above the surface, and the dielectric constant of the underground background medium is ε r =4 and σ = 3mS / m. There are 4 anomalies at a depth of 0.72m underground, each with a size of 0.72m×0.72m×0.72m. The electrical parameters of the two blue anomalies are ε r =1 and σ=0; the red anomaly is ε r=8 and σ = 10mS / m, the computational domain consists of 50×50×50 grids, the grid size is 0.12m×0.12m×0.12m, and the PML thickness is 0.48m. The excitation source of the forward simulation uses the Ricker wavelet with a main frequency of 200MHz. 121 sources and 420 receivers (represented by black and white “·” respectively) are evenly distributed in the air layer 0.12m below the surface, with spacing of 0.48m and 0.24m respectively. The signal of each source will be recorded by all receivers.

[0139] The anomaly has a large contrast, and the multiple scattering between the four anomalies and the ground observation method increase the complexity of the inversion. 5% Gaussian noise is added to the observation data. The initial inversion model is a background homogeneous medium. The parameters of the air layer in the model are limited during the inversion process, while avoiding the singularity at the source and receiver positions. The parameters of the underground dielectric constant and conductivity to be determined are 210,000.

[0140] The full waveform inversion test was performed on a Dell Precision T7920 workstation with dual Intel(R) Xeon(R) Platinum 8168 CPUs @ 2.70GHz, 512GB of physical memory, and Windows 10 operating system. The forward simulation and gradient calculation modules in the inversion process were accelerated by CPU parallelism in frequency. Ten frequency data (i.e., 50, 60, 70, 80, 100, 120, 150, 180, 200, and 220MHz) were selected for synchronous multi-frequency inversion. The model results of relative permittivity and conductivity after 100 iterations are shown in Figure 2. Figure 3 As shown in the figure, TV regularization is used in the inversion process, and the parameter adjustment factor κ is set to 1.2 to ensure that the dielectric constant plays a major role in the inversion process. The inversion takes about 16.87 hours, and the maximum memory used in the entire algorithm is 300GB.

[0141] like Figure 3 a and Figure 3 b is a schematic diagram of the full waveform inversion results of model 1 of an embodiment of the present invention. Figure 3 a represents the full waveform inversion result of the dielectric constant of model 1 of the embodiment of the present invention; Figure 3 b represents the full waveform inversion result of the conductivity of model 1 of the embodiment of the present invention. Figure 3 a and Figure 3 As shown in b, the approximate shapes and physical parameters of the four anomalies were successfully reconstructed. Figure 2 a and Figure 3 a: The interface of the dielectric constant inversion result is relatively clear. The inversion result of the low dielectric constant anomaly is basically consistent with the true value, while the high dielectric constant anomaly is somewhat different from the true value result; opposite artifacts appear above the four anomalies. Figure 2 b and Figure 3 b: The volume of the conductivity anomaly has changed to a certain extent, and the center position has moved upward. Relatively speaking, the reconstruction of the dielectric constant is more accurate than the reconstruction of the conductivity, especially the interface of the anomaly is clearer.

[0142] like Figure 4 a and Figure 4 b is a schematic diagram of relative data residual and reconstruction error of model 1 of an embodiment of the present invention. Figure 4 a represents the relative data residual of model 1 of the embodiment of the present invention; Figure 4 b represents the reconstruction error of model 1 of the embodiment of the present invention. Figure 4 In a, the relative data residual decreases rapidly in the early stage of iteration, decreases toward the minimum value in the middle and late stages of iteration, and the convergence speed is more stable. Figure 4 The blue solid line in (b) shows a smaller dielectric constant reconstruction error curve than the red dotted line. According to the trend of the two convergence curves, more accurate reconstruction can be achieved after more iterations. However, since the physical parameters of the four anomalies have been well reconstructed, more iterations will take more time. The inversion results show that the synchronous multi-frequency strategy of the frequency domain weighted method can effectively reconstruct the three-dimensional distribution of the underground medium.

[0143] The advantages of the method of the present invention are as follows:

[0144] The inversion results of the above case show that the method of the present invention has good adaptability to ground multi-offset noise data. Not only can the dielectric constant and conductivity be inverted synchronously to improve the inversion accuracy, but also the clarity of the abnormal body interface is improved. The following is a further explanation combined with a case:

[0145] like Figure 5 a and Figure 5 b is a schematic diagram of a three-dimensional model of model 2 of an embodiment of the present invention. Figure 5 a is a schematic diagram of a three-dimensional model of the dielectric constant of model 2 according to an embodiment of the present invention; Figure 5 b is a schematic diagram of a three-dimensional conductivity model of Model 2 of the embodiment of the present invention. In order to observe the inversion effect of this algorithm for abnormal bodies at different depths, Figure 2 a and Figure 2 Based on the model 1 shown in b, the burial depth of two anomalies was changed to 1.08m, the other two remained unchanged, and the distance between the four anomalies was further reduced. The schematic diagram of 3D model 2 is shown in Figure 5As shown. Other segmentation and measurement parameters are consistent with 3D model 1. Gaussian noise with SNR = 30dB is added to the total field synthetic data as observation data, and other inversion parameters and strategies are consistent with 3D model 1. This algorithm uses the MTV algorithm to improve inversion, so on this basis, standard Tikhonov regularization and MTV regularization inversion experiments are carried out to compare the effects of different regularization methods on noisy data.

[0146] like Figure 6 a. Figure 6 b. Figure 6 c and Figure 6 d is a schematic diagram of the full waveform inversion results under different regularization strategies of Model 2 of an embodiment of the present invention. Figure 6 a is a schematic diagram showing a three-dimensional model of a dielectric constant using MTV regularization in Model 2 of an embodiment of the present invention; Figure 6 b is a schematic diagram of a three-dimensional model of a dielectric constant using standard Tikhonov regularization of Model 2 of an embodiment of the present invention; Figure 6 c is a schematic diagram of a three-dimensional model of conductivity using MTV regularization in Model 2 of an embodiment of the present invention; Figure 6 d is a schematic diagram of a three-dimensional model of conductivity using standard Tikhonov regularization for model 2 of the embodiment of the present invention. The inversion time of MTV regularization inversion is about 15.91 hours, and the inversion time of standard Tikhonov regularization inversion is about 21.57 hours (36% longer than the MTV regularization method). The maximum memory used in the entire algorithm is 300GB.

[0147] observe Figure 6 It can be found that Figure 6 a- Figure 6 The abnormal body interface in d is not as good as Figure 5 The interface of the real model in is clear. Figure 6 a and Figure 6 The interface in c is obviously Figure 6 b and Figure 6 The interface in d is clear, which shows that FWI and MTV regularization are more advantageous than standard Tikhonov regularization in maintaining sharp interfaces. Figure 6 a. Figure 6 b and Figure 6 c. Figure 6 The results in d show that, relatively speaking, the reconstruction of the dielectric constant is more accurate than that of the conductivity. This is because the data is more sensitive to the dielectric constant than to the conductivity. Figure 6 c. Figure 6 d and Figure 5 Compared with the real conductivity model of b, it is found that the conductivity anomaly has a certain change in volume and the center position has moved upward. In addition, the inversion parameters of the two shallow anomalies are closer to the real values.

[0148] It is worth pointing out that the background value of the inversion result of MTV regularization is closer to the true model compared with the standard Tikhonov regularization. In general, there is no obvious artificial noise in the inversion results, which verifies the robustness of the method. In order to further quantitatively evaluate the accuracy of the inversion results, two different reconstruction error indicators are used:

[0149]

[0150]

[0151] in, represents the standard reconstruction error; represents the relative reconstruction error; represents the medium parameter vector in the real model; Represents the medium parameter vector obtained by the final inversion.

[0152] Table 1. The number of iterations and model error parameters of FWI under different regularization strategies

[0153] method Iterations <![CDATA[MSE(ε r )]]> MSE(σ) <![CDATA[Err m (e r )]]> <![CDATA[Err m (s)]]> MTV 72 0.058 0.149 0.117 0.240 Tikhonov 71 0.061 0.172 0.252 0.595

[0154] The iteration times and model error parameters of the two regularization schemes are listed in Table 1. The results show that the dielectric constant reconstruction error of the two regularization schemes is smaller than the conductivity reconstruction error; compared with the standard Tikhonov regularization, the MTV regularized FWI has a lower reconstruction error.

[0155] The inversion results under the two regularization strategies show that the MTV regularization-based method can improve the accuracy of dielectric constant and conductivity inversion, and is superior to the classic Tikhonov regularization method in terms of the clarity of the anomalous body interface.

Claims

1. A three-dimensional frequency domain ground penetrating radar dual-parameter synchronous inversion method, characterized in that The steps include: S1. Input observation data, set up the 3D initial model, and initialize parameters; S2. Forward modeling of the established three-dimensional initial model; S3. construct an inversion objective function and decompose it into two alternating sub-problems, including sub-problem 1 and sub-problem 2, sub-problem 1 is to solve the model medium parameter vector, and sub-problem 2 is to solve the pilot model vector; The sub-question 1 is specifically: in, represents the model medium parameter vector, Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the model medium parameter vector after the kth iteration; represents the inversion objective function with respect to the model medium parameter vector; represents the data objective function; ε represents the current dielectric constant; σ represents the current conductivity; and represents the regularization parameter with a positive value at the kth iteration; represents the leading model vector corresponding to the current dielectric constant after the k-1th iteration, represents the leading model vector corresponding to the current conductivity after the k-1th iteration; where the regularization parameter of the kth iteration and They are: and in, and is a dimensionless number; represents the data objective function; Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the leading model vector corresponding to the current dielectric constant after the k-1th iteration, represents the leading model vector corresponding to the current conductivity after the k-1th iteration; The sub-question 2 includes: Among them, ε represents the current dielectric constant; σ represents the current conductivity; β σ and β ε represents a regularization parameter with a positive value; TV represents the total variation operator; represents the leading model vector corresponding to the current dielectric constant after the kth iteration; Represents the leading model vector corresponding to the current dielectric constant; represents the objective function of the leading model vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current dielectric constant after the kth iteration; represents the leading model vector corresponding to the current dielectric constant after the kth iteration; Represents the leading model vector corresponding to the current dielectric constant; represents the objective function of the leader model vector corresponding to the current conductivity; represents the model medium parameter vector corresponding to the current dielectric constant after the kth iteration; among them, the regularization parameter of the kth iteration and They are: in, and is a dimensionless number; represents the data objective function; Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the leading model vector corresponding to the current dielectric constant after the k-1th iteration, represents the leading model vector corresponding to the current conductivity after the k-1th iteration; S4. Initialize the dielectric constant and conductivity to obtain an initial pilot model vector; S5. Solve subproblem 1 according to the pilot model vector and update the model medium parameter vector; solve subproblem 1 using the limited memory BFGS method, and for the model medium parameter vector after the k+1th iteration in, represents the model medium parameter vector after the kth iteration; a (k ) is the search step size of the kth iteration; represents the search direction of the kth iteration, represents an approximation of the inverse Hessian matrix; represents the objective function of the model medium parameter vector for the kth iteration The gradient of The objective function of the model medium parameter vector for the kth iteration Gradient Specifically: in, represents the gradient of the objective function with respect to the model medium parameter vector; represents the data objective function; and Indicates that the value of the k-th iteration is positive regularization parameter; Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; Represents the prior model vector corresponding to the current dielectric constant after the k-1th iteration; Represents the prior model vector corresponding to the current conductivity after the k-1th iteration; Represents the gradient of the data objective function with respect to the model medium parameter vector corresponding to the current dielectric constant; represents the gradient of the data objective function with respect to the model medium parameter vector corresponding to the current conductivity; For the underground medium unit e, the gradient of the data objective function of the model medium parameter vector corresponding to the current dielectric constant Adjoint state method is used to obtain: in, is the model medium parameter vector; N ω Indicates the maximum number of angular frequencies; N s Indicates the maximum number of excitation sources; i indicates the frequency number; r indicates the excitation source number; is the difference between the observed data and the simulated data, specifically represents observation data; represents simulated data; T represents transpose, and * represents conjugation; represents the real part operator; represents discrete electric field values; represents the impedance matrix; The corresponding position matrix of the receiving point representing the current source; is a sparse matrix, Impedance matrix The model medium parameter m corresponding to the current unit e e sensitivity, is a diagonal matrix, the coefficient corresponding to the current unit e is a non-zero element, through Calculated; its model medium parameter vector m e Including dielectric constant and conductivity; Solving the cell matrix For the model medium parameter corresponding to the current dielectric constant ε in the corresponding current unit e Sensitivity Solving the cell matrix For the model medium parameter corresponding to the current conductivity σ in the corresponding current unit e Sensitivity and in, represents the model medium parameter corresponding to the current dielectric constant ε in the current unit e; ω represents the angular frequency; j represents the imaginary unit, ε0 is the dielectric constant of free space; For vector basis functions, when the computational domain is discretized using hexahedral units: l,n=1,2,…,12; Represents the volume integral within the current unit e; represents the model medium parameter corresponding to the current conductivity σ in the current unit e; β represents a regularization parameter with a positive value; σ0 is the model background conductivity; S6. Solve subproblem 2 according to the obtained model medium parameter vector and update the pilot model vector; specifically, use the Split-Bregman iterative method to solve subproblem 2 and update it: Subproblem 2 is expressed using the unified TV model: in, represents the prior model vector, Represents the prior model vector corresponding to the current dielectric constant, Represents the prior model vector corresponding to the current conductivity; represents the model medium parameter vector after the kth iteration; Represents the model medium parameter vector corresponding to the current dielectric constant after the kth iteration; Represents the model medium parameter vector corresponding to the current conductivity after the kth iteration; TV represents the total variation operator; μ TV =μ TV,ε ,μ TV,σ , while μ TV , μ TV,ε and μ TV,σ represents the regularization parameter; In the total variation operator TV In the example, for a 3D isotropic regularization operator It is expressed as: Among them, e represents the current unit; Represents the gradient of the prior model vector on the x-axis; Represents the gradient of the prior model vector on the y-axis; Represents the gradient of the prior model vector on the z-axis; sets auxiliary variables and The L1 constraint and L2 penalty components of the segmentation problem give the optimality condition of the Bregman formula for the unified TV model splitting: Where λ = 2μ TV ;μ TV =μ TV,ε ,μ TV,σ , while μ TV , μ TV,ε and μ TV,σ represents the regularization parameter; represents the identity matrix; represents the model medium parameter vector after the kth iteration; Indicates the three directions of x, y, and z in three dimensions; represents the kth iteration, Auxiliary variables in direction; represents the kth iteration, The Bregman coefficient in the direction, its iterative formula is represents the k+1th iteration, Bregman coefficient in direction; represents the k+1th iteration, Auxiliary variables in direction; express and The vector difference is The transposed matrix of the gradient in the direction; Δ is the Laplace operator; represents the prior model vector after the k+1th iteration; set up Auxiliary variables Solve explicitly using the generalized contraction formula: Where λ = 2μ TV ;μ TV =μ TV,ε ,μ TV,σ , while μ TV , μ TV,ε and μ TV,σ represents the regularization parameter; represents the shrinkage coefficient, Represents the prior model vector after the kth iteration Directional gradient; Represents the gradient of the prior model vector in the x direction after the kth iteration; Represents the gradient of the prior model vector in the y direction after the kth iteration; Represents the gradient of the prior model vector in the z direction after the kth iteration; represents the kth iteration, Bregman coefficient in direction; represents the Bregman coefficient in the x direction of the kth iteration; represents the Bregman coefficient in the y direction of the kth iteration; represents the Bregman coefficient in the z direction at the kth iteration; S7. Repeat steps S5-S6, and update synchronously and alternately until the inversion termination condition is met, then end the inversion and output the inversion result.

2. The three-dimensional frequency domain ground penetrating radar dual parameter synchronous inversion method according to claim 1 is characterized in that The step S2 includes forward modeling the established three-dimensional initial model using the finite element method, specifically using the linear vector basis function finite element method to discretize the simulation area, dividing the grid into rectangular block units, and using anisotropic perfectly matched layers at the truncation boundary; the frequency domain Maxwell equations satisfied by the three-dimensional GPR are expressed as: in, represents the electric field strength; ε represents the current dielectric constant; μ represents the magnetic permeability; σ represents the current conductivity; ω represents the angular frequency, ω=2πf, f represents the frequency; j represents the imaginary unit, is the current density; Find the linear equations for a multi-source system: in, represents the coefficient matrix common to all frequency-dependent sources; represents the discrete electric field value; represents the discretized excitation source term; ω represents the angular frequency; ε represents the dielectric constant; σ represents the conductivity; s represents the excitation source at the current position; the coefficient matrix common to all sources related to frequency By all unit matrices synthesis: in, Represents the volume integral within the current unit e; and is a vector basis function; when the calculation area is discretized using hexahedral units: l,n=1,2,…,12; for the matrix Perform an LU decomposition; the forward operator for a single source and a single frequency is described as: in, Represents the simulation data under the angular frequency ω and the current excitation source s; The corresponding position matrix of the receiving point representing the current source; represents the discrete electric field value; represents the coefficient matrix common to all frequency-dependent sources; represents the discretized excitation source term.

3. The three-dimensional frequency domain ground penetrating radar dual-parameter synchronous inversion method according to claim 2 is characterized in that Step S3, using known measured data to reconstruct the model parameters in the underground medium, including the spatial distribution of the dielectric constant ε and the conductivity σ; and introducing the total variation regularization method to invert the objective function Defined as: in, represents the model medium parameter vector, Represents the model medium parameter vector corresponding to the current dielectric constant; Represents the model medium parameter vector corresponding to the current conductivity; represents the prior model vector, Represents the prior model vector corresponding to the current dielectric constant, Represents the prior model vector corresponding to the current conductivity; represents the data objective function; ε represents the current dielectric constant; σ represents the current conductivity; α ε , α σ , β ε and β σ represents a regularization parameter with a positive value; TV represents the total variation operator; the inversion objective function It can be further decomposed into subproblem 1 and subproblem 2. Subproblem 1 is to solve the model medium parameter vector, and subproblem 2 is to solve the prior model vector.

4. The three-dimensional frequency domain ground penetrating radar dual parameter synchronous inversion method according to claim 1 is characterized in that The step S7 is specifically to repeat steps S5-S6, synchronously and alternately update sub-problem 1 and sub-problem 2, until the set inversion termination condition is met, end the inversion and output the inversion result: the set inversion termination condition includes reaching the set iteration accuracy or reaching the maximum number of iterations; The modified total variation problem solves subproblems 1 and 2 alternately; in the kth iteration, according to step S4, the variable and Solve the model medium parameter vector after the kth iteration and represents the leading model vector corresponding to the current dielectric constant after the k-1th iteration; represents the pilot model vector corresponding to the current conductivity after the k-1th iteration; then according to step S5, the model medium parameter vector after the kth iteration is used and Solving for variables and represents the leading model vector corresponding to the current dielectric constant after the kth iteration; represents the leading model vector corresponding to the current conductivity after the kth iteration; the obtained variable and Used in the next iteration step; for the first iteration step, the pilot model vector corresponding to the initial current dielectric constant Equal to the model medium parameter vector corresponding to the initial current dielectric constant The leading model vector corresponding to the initial current conductivity Equal to the model medium parameter vector corresponding to the initial current conductivity Get the inversion sequence:

Images