A method for ground penetrating radar parameter inversion based on non-monotonic line search

Abstract

The invention discloses a method for inverting parameters of a ground penetrating radar based on a non-monotonic line search, comprising establishing a real model and an initial model, constructing a target function with the real model and the initial model, obtaining the gradient of a dielectric constant and a gradient of conductivity in the function, and calculating the step length of the dielectric constant and the step length of the conductivity, updating the dielectric constant model and the conductivity model, and re-extracting data based on the new model as new theoretical data to recalculate the target function value; judging whether the new target function value satisfies a non-monotonic linear search condition and judging whether the inversion process meets a termination condition. The invention uses a non-monotonic decreasing model updating method, which only requires the target function value to decrease within a finite number of steps, thereby improving the convergence rate of the inversion. At the same time, the looser convergence condition in the method also reduces the possibility of the inversion falling into a local minimum point, thereby ensuring the stability and timeliness of the full waveform inversion of the ground penetrating radar.

Description

Technical Field

[0001] The invention relates to the technical field of ground penetrating radar imaging, and in particular to a ground penetrating radar parameter inversion method based on non-monotonic line search. Background Art

[0002] Ground penetrating radar is an efficient and high-resolution nondestructive testing method that can be applied to underground structure detection, geological disaster assessment, soil geological information assessment, hydrocarbon pollution detection and many other aspects. Ground penetrating radar can identify abnormal bodies by picking up typical features on the radar profile and assisting with manual judgment, and then determine the abnormal distribution of the detected body. However, this method is not suitable for complex underground structures. When there are multiple morphological targets underground, the profile features will overlap with each other, making it difficult to judge the abnormal bodies. Sometimes, it is even difficult to identify due to the similar features of different target bodies, so it is impossible to obtain the accurate distribution of underground target bodies. The parameter inversion method performs underground imaging based on the electrical differences of the medium, which can clearly reflect the spatial distribution of the detected body. Among them, the full waveform inversion method is a relatively popular parameter inversion method in recent years.

[0003] Ernst (2006) realized the full waveform inversion of ground penetrating radar in the time domain and accurately located the sub-wavelength anomaly; Meles (2010) used the full waveform inversion of borehole radar to obtain a clearer underground stratigraphic structure by the river; Hu Zhouwen (2018) used full waveform inversion to perform non-destructive testing on concrete, and Feng Deshan (2018) et al. used GPU parallelism to realize dual-parameter inversion of radar data. It can be seen that full waveform inversion of radar data is feasible.

[0004] Although full waveform inversion can improve the resolution of ground penetrating radar, the inversion itself has certain limitations. On the one hand, full waveform inversion uses a local optimization algorithm to optimize the model. Its disadvantage is that when the initial model of the inversion is not accurate enough, the full waveform inversion will fall into a local minimum due to the large difference between the initial waveform and the theoretical waveform, resulting in inversion failure. On the other hand, due to the large amount of data involved in full waveform inversion, its calculation itself takes a long time. In addition, the update of the model often requires a monotonic decrease, resulting in low computational efficiency of the inversion, which does not meet the timeliness requirements of the ground penetrating radar. Therefore, a ground penetrating radar parameter inversion method based on non-monotonic line search is urgently needed to solve the above problems. Summary of the invention

[0005] The present invention provides a method for ground penetrating radar parameter inversion based on non-monotonic line search, which can improve the calculation efficiency of inversion and avoid the inversion falling into a local minimum point. The method can effectively solve the problem that the inversion of the prior art is prone to fall into a local minimum point and cannot meet the timeliness of the radar as mentioned in the above background technology.

[0006] To achieve the above object, the present invention provides the following technical solution: a method for ground penetrating radar parameter inversion based on non-monotonic line search, comprising the following steps:

[0007] Step 1: Establish a real model of ground penetrating radar full waveform inversion, and perform forward simulation based on the real model and extract gathers to obtain observation data;

[0008] Step 2: Establish an initial model for full waveform inversion of ground penetrating radar and set the inversion termination accuracy. Perform forward simulation based on the initial model and extract gathers to obtain theoretical data.

[0009] Step 3: Construct the objective function and express it in the form of a two-norm:

[0010]

[0011] Among them, D cal represents theoretical data, D obs represents observation data;

[0012] Step 4: Calculate the first-order derivative of the objective function, obtain the gradient of the dielectric constant and the gradient of the conductivity, and calculate the step length of the dielectric constant and the step length of the conductivity;

[0013] Step 5: Update the dielectric constant model and the conductivity model, re-extract data based on the new model as new theoretical data, and recalculate the objective function value;

[0014] Step 6: Determine whether the new objective function value satisfies the non-monotonic linear search condition. If it does not satisfy the search condition, recalculate the step size until it satisfies the search condition.

[0015] Step 7: Determine whether the inversion process meets the inversion termination accuracy condition set in step 2. If the termination condition is not met, update the gradient of the dielectric constant and conductivity, repeat steps 4 to 6 until the termination condition is met, and then output the final inversion result.

[0016] According to the above technical solution, the real model of the full waveform inversion of the ground penetrating radar is a partial marmousi model, which is assigned values ​​according to the numerical variation range of the dielectric constant and the conductivity. When the real model is forward simulated, the frequency domain finite difference method based on the staggered grid is used for numerical simulation, and the boundary adopts the fully matched layer boundary condition.

[0017] According to the above technical solution, in step 2, the initial model is obtained after smoothing.

[0018] According to the above technical solution, in step 4, the gradient expressions of dielectric constant and conductivity are as follows:

[0019] The gradient expression of the dielectric constant is,

[0020] The gradient expression of conductivity is:

[0021] Where T represents transpose, A is the sparse operator, and Δd is the difference between simulated data and observed data.

[0022] According to the above technical solution, sequential inversion is adopted for the two parameters of dielectric constant and conductivity, that is, when inverting one of them, the other parameter is taken as a constant.

[0023] According to the above technical solution, in step 4, when calculating the step length, the initial parameters α, δ, ρ and ε are set, where α>0, δ>0, ε>0 and 0<ρ<1, and the step length α of the dielectric constant is k and the conductivity step size α k Satisfy the formula

[0024] According to the above technical solution, in step 5, according to the formula Update the dielectric constant model and conductivity model, where x k represents the dielectric constant and conductivity of the kth iteration.

[0025] According to the above technical solution, in step 6, judging whether the non-monotonic linear search condition is met is specifically to judge whether the objective function satisfies the inequality If this inequality does not hold, then by adjusting h k Until the step size α k Satisfy the above formula to determine the final step size, where h k is to make α k The smallest non-negative integer that satisfies the above inequality.

[0026] In the inversion process, the purpose of inversion is to make the objective function drop to the allowable error range, but it does not require the objective function to be strictly monotonically decreased. Blindly requiring the objective function to be monotonically decreased will reduce the convergence speed of the objective function and increase the amount of calculation;

[0027] Compared with the prior art, the present invention uses a non-monotonic decreasing model updating method, which only requires the objective function value to decrease within a finite number of steps, thereby improving the convergence rate of the inversion. At the same time, the looser convergence conditions in this method also reduce the possibility of the inversion falling into a local minimum, thereby ensuring the stability and timeliness of the full waveform inversion of the ground penetrating radar. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] The accompanying drawings are used to provide further understanding of the present invention and constitute a part of the specification. They are used to explain the present invention together with the embodiments of the present invention and do not constitute a limitation of the present invention.

[0029] In the attached picture:

[0030] Figure 1 It is a schematic diagram of the process of the present invention;

[0031] Figure 2 The real model diagram, the initial model diagram and the inversion model diagram of the dielectric constant of the present invention are;

[0032] Figure 3 The real model diagram, the initial model diagram and the inversion model diagram of the electrical conductivity of the present invention are;

[0033] Figure 4 It is a schematic diagram of the non-monotonic linear search method of the present invention. DETAILED DESCRIPTION

[0034] The preferred embodiments of the present invention are described below in conjunction with the accompanying drawings. It should be understood that the preferred embodiments described herein are only used to illustrate and explain the present invention, and are not used to limit the present invention.

[0035] Example: Figure 1 As shown, a method for ground penetrating radar parameter inversion based on non-monotonic line search includes the following steps:

[0036] Step 1: Establish a realistic model of full waveform inversion of ground penetrating radar to simulate actual data collection;

[0037] Among them, the true model of the full waveform inversion of the ground penetrating radar is a partial marmousi model, which is assigned according to the numerical variation range of the dielectric constant and the conductivity. The true model of the dielectric constant and the numerical true model of the conductivity are as follows: Figure 2 (a) and Figure 3 (b)

[0038] in, Figure 2 and Figure 3 The actual depth of the model in is 1.5m×2m, which is discretized into 30×40 grids with a grid spacing of 0.05m.

[0039] Based on the real model, forward simulation is carried out and gathers are extracted to obtain observation data, which is equivalent to the radar data collected in the actual project. When forward modeling the real model, the frequency domain finite difference method based on staggered grids is used for numerical simulation, and the boundary adopts the fully matched layer boundary condition;

[0040] Analogous to the equation form of the forward part in seismic full waveform inversion, the wave equation of ground penetrating radar full waveform inversion can be written as: A(m1 m2)U=S, where A is the forward operator, U is the wave field, S is the field source, m1 and m2 are the dielectric constant and conductivity parameters respectively; when m1 and m2 both represent the real model, the gather extracted from the wave field U is the observed radar data.

[0041] Step 2: Establish the initial model of the ground penetrating radar full waveform inversion and set the inversion termination accuracy. The initial model is obtained after smoothing. The initial model of dielectric constant and the initial model of conductivity are as follows: Figure 2 (b) and Figure 3 (b)

[0042] Based on the initial model, forward simulation is performed and gathers are extracted to obtain theoretical data. The principle is similar to that of step 1. The forward equation is still A(m1 m2)U=S. When m1 and m2 both represent the initial model, the theoretical radar data is obtained by extracting gathers from the wave field U.

[0043] Step 3: Construct the objective function and express it in the form of a two-norm:

[0044]

[0045] Among them, D cal represents theoretical data, D obs represents observation data;

[0046] Step 4: Calculate the first-order derivative of the objective function to obtain the gradient of the dielectric constant and the gradient of the conductivity;

[0047] The gradient expression of the dielectric constant is,

[0048] The gradient expression of conductivity is:

[0049] Where T represents transpose, A is the sparse operator, and Δd is the difference between simulated data and observed data;

[0050] When calculating the dielectric constant gradient, the conductivity is fixed as a constant, and when calculating the conductivity gradient, the dielectric constant is considered as a constant, e.g., when calculating the dielectric constant gradient, the conductivity is fixed as a constant, and when calculating the conductivity gradient, the dielectric constant is considered as a constant;

[0051] Set the initial parameters α, δ, ρ and ε, where α>0, δ>0, ε>0 and 0<ρ<1, and calculate the step size α of the dielectric constant k and the conductivity step size α k , satisfying the formula

[0052] Step 5: According to the formula Update the dielectric constant model and conductivity model, where x k The dielectric constant and conductivity of the k-th iteration are expressed, and the data are re-extracted based on the new model as new theoretical data to recalculate the objective function value;

[0053] Step 6: Determine whether the new objective function value satisfies the non-monotonic linear search condition. If it does not meet the search condition, recalculate the step size until it meets the search condition. Specifically: Determine whether the objective function satisfies the inequality If this inequality does not hold, then by adjusting h k Until the step size α k Satisfy the above formula to determine the final step size, where h k is to make α k The smallest non-negative integer that satisfies the above inequality, and proceed to step 7 after the conditions are met;

[0054] Step 7: Determine whether the inversion process meets the inversion termination accuracy condition set in step 2. If the termination condition is not met, set k = k + 1, update the gradients of the dielectric constant and conductivity, repeat steps 4 to 6 until the termination condition is met, output the final inversion result, and the inversion ends. Figure 2 (c) and Figure 3 (c) shows the inversion model of dielectric constant and conductivity respectively.

[0055] In the inversion process, the purpose of inversion is to make the objective function drop to the allowable error range, but it does not require the objective function to be strictly monotonically decreased. Blindly requiring the monotonically decreased objective function will reduce the convergence speed of the objective function and increase the amount of calculation. Therefore, the present invention uses a non-monotonic linear search method. In addition to improving the inversion convergence rate, this method can also reduce the possibility of the inversion falling into a local minimum point, such as Figure 4 As shown in the figure, it is a schematic diagram of the non-monotonic linear search method. In the figure, X m1 and X m2 They represent the local minimum and the global minimum respectively, x1 represents the initial model, f(x) represents the objective function value, and for the monotonic line search method, since it requires the objective function value to decrease monotonically, that is, to satisfy f(x1)>f(x2)>f(x3)>f(x4), when the model is updated to X m1 When , the inversion will terminate because its objective function value will not continue to decrease, but the corresponding point at this time is not the global minimum point; while for the non-monotone line search method, its objective function only needs to satisfy max{f(x1),f(x2),f(x3)}>f(x4), so that the objective function has the opportunity to invert to the global minimum point, thereby ensuring the stability and timeliness of the full waveform inversion of the ground penetrating radar.

[0056] Finally, it should be noted that the above description is only a preferred example of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions described in the aforementioned embodiments or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A method for ground penetrating radar parameter inversion based on non-monotonic line search, characterized in that: The steps include: Step 1: Establish a real model of ground penetrating radar full waveform inversion, and perform forward simulation based on the real model and extract gathers to obtain observation data; Step 2: Establish an initial model for full waveform inversion of ground penetrating radar and set the inversion termination accuracy. Perform forward simulation based on the initial model and extract gathers to obtain theoretical data. Step 3: Construct the objective function and express it in the form of a two-norm: Among them, D cal represents theoretical data, D obs represents observation data; Step 4: Calculate the first-order derivative of the objective function, and then calculate the gradient of the dielectric constant and the gradient of the conductivity, as well as the step length of the dielectric constant and the step length of the conductivity; Step 5: Update the dielectric constant model and the conductivity model, re-extract data based on the new model as new theoretical data, and recalculate the objective function value; Step 6: Determine whether the new objective function value satisfies the non-monotonic linear search condition. If it does not satisfy the search condition, recalculate the step size until it satisfies the search condition. To determine whether the non-monotonic linear search condition is met, specifically, to determine whether the objective function satisfies the inequality If this inequality does not hold, then by adjusting h k Until the step size α k Satisfy the above formula to determine the final step size, where h k is to make α k The smallest non-negative integer that satisfies the above inequality; Step 7: Determine whether the inversion process meets the inversion termination accuracy condition set in step 2. If the termination condition is not met, update the gradient of the dielectric constant and conductivity, repeat steps 4 to 6 until the termination condition is met, and then output the final inversion result.

2. The method for ground penetrating radar parameter inversion based on non-monotonic line search according to claim 1, characterized in that: The real model of the ground penetrating radar full waveform inversion is a partial marmousi model, which is assigned according to the numerical variation range of the dielectric constant and conductivity. When the real model is forward simulated, the frequency domain finite difference method based on the staggered grid is used for numerical simulation, and the boundary adopts the perfectly matched layer boundary condition.

3. The method for ground penetrating radar parameter inversion based on non-monotonic line search according to claim 1, characterized in that: In step 2, the initial model is obtained after smoothing.

4. The method for ground penetrating radar parameter inversion based on non-monotonic line search according to claim 1, characterized in that: In step 4, the gradient expressions of dielectric constant and conductivity are as follows: The gradient expression of the dielectric constant is, The gradient expression of conductivity is: Where T represents transpose, A is the sparse operator, and Δd is the difference between simulated data and observed data.

5. The method for ground penetrating radar parameter inversion based on non-monotonic line search according to claim 4, characterized in that: Sequential inversion is used for the two parameters of dielectric constant and conductivity, that is, when inverting one of them, the other parameter is taken as a constant.

6. The method for ground penetrating radar parameter inversion based on non-monotonic line search according to claim 1, characterized in that: In step 4, when calculating the step size, the initial parameters α, δ, ρ and ε are set, where α>0, δ>0, ε>0 and 0<ρ<1, and the step size α of the dielectric constant is k and the conductivity step size α k Satisfy the formula 7. The method for ground penetrating radar parameter inversion based on non-monotonic line search according to claim 1, characterized in that: In step 5, according to the formula Update the dielectric constant model and conductivity model, where x k represents the dielectric constant and conductivity of the kth iteration.

Images