Ground penetrating radar reverse time migration imaging method, system and device based on local Nyquist attenuation compensation and medium

Abstract

The invention discloses a ground penetrating radar reverse time migration imaging method, system and device based on local Nyquist attenuation compensation and a medium. The method comprises the steps that ground penetrating radar observation data are collected, and a speed model is constructed; calculating to obtain a local time window; executing forward propagation wave field simulation to obtain a sparse forward propagation wave field; executing back propagation wave field simulation according to the wave equation, and performing symbol inversion on a conductivity item in the wave equation to obtain a compensation back propagation wave field; performing cross-correlation calculation to obtain an imaging value; and generating a reverse time migration imaging profile. According to the invention, through cooperation of triple mechanisms of a local time window, Nyquist sparse sampling and attenuation compensation operation, the wave field storage demand and the calculation burden are reduced while the imaging precision is maintained; and the energy attenuation of electromagnetic waves in the conducting medium is compensated, the imaging energy of deep reflection signals is enhanced, the technical problems of high storage cost and insufficient deep imaging capability of a traditional method are solved, and a solution is provided for high-precision ground penetrating radar imaging under complex geological conditions.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration technology, and in particular to a ground-penetrating radar reverse time migration imaging method, system, device and medium based on local Nyquist attenuation compensation. Background Technology

[0002] Ground penetrating radar (GPR) is a geophysical technique that uses high-frequency electromagnetic waves to non-destructively probe underground structures. It is widely used in geological exploration, civil engineering inspection, and underground target identification. Reverse time migration (RTM), a two-way wave imaging method based on the full wave equation, can accurately describe complex wave phenomena such as electromagnetic wave propagation paths and waveform interference. Therefore, it has become a key technology for achieving high-precision, high-resolution GPR imaging, showing significant advantages, especially in the detection of underground cavities, buried objects, and complex layered structures.

[0003] However, applying reverse time migration (RTM) technology to ground-penetrating radar (GPR) data processing currently faces two major technical bottlenecks. First, imaging accuracy heavily relies on the accuracy of the velocity model; even small errors in the velocity model can lead to shifts in imaging position and energy distortion. Constructing a high-precision velocity model in complex media remains a challenge. Second, traditional cross-correlation imaging conditions require preserving the entire source wave field during forward propagation, resulting in significant memory overhead and input / output burdens. This makes computationally inefficient when processing large-scale two-dimensional or even three-dimensional data, hindering practical application. Although existing research has proposed sparse sampling strategies such as excitation amplitude imaging conditions to reduce storage costs, these methods struggle to effectively recover weak deep reflection signals, especially in environments with severe signal attenuation or strong radar noise, particularly in high-conductivity strata, leading to decreased imaging stability. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides a ground-penetrating radar reverse time migration imaging method, system, device, and medium based on local Nyquist attenuation compensation to address the problems of current ground-penetrating radar reverse time migration technology, such as the imaging accuracy being extremely sensitive to velocity model errors and the challenges in constructing the velocity model itself; and the need to store the entire wave field for traditional cross-correlation imaging, resulting in huge storage overhead and input / output burden, making it inefficient and impractical for large-scale data processing.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution: In a first aspect, the present invention provides a ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation, comprising: Collect ground-penetrating radar observation data and construct a velocity model; The imaging region is obtained by combining the velocity model with the observation geometry of the ground-penetrating radar observation data, and a corresponding local time window is calculated for each imaging point within the imaging region. A forward propagation wavefield simulation is performed based on the wave equation to obtain a full forward propagation wavefield. The full forward propagation wavefield and the local time window are then sparsely sampled and stored according to the Nyquist sampling law to obtain a sparse forward propagation wavefield. Based on the wave equation, a reverse propagation wave field simulation is performed, and the conductivity term in the wave equation is reversed to achieve attenuation compensation, thus obtaining a compensated reverse propagation wave field. The sparse forward propagation wavefield and the compensated reverse propagation wavefield are cross-correlated within the local time window to obtain the imaging values ​​of each imaging point. The imaging values ​​of all imaging points are then integrated to generate a reverse time-shifted imaging profile.

[0007] As a preferred embodiment of the ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation described in this invention, the steps of acquiring ground-penetrating radar observation data and constructing a velocity model include: Collect ground-penetrating radar observation data; Based on the ground-penetrating radar observation data, a velocity model of dielectric constant and conductivity is constructed.

[0008] The beneficial effects of this preferred technical solution are as follows: it constructs a velocity model that includes dielectric constant and conductivity, providing a physical parameter basis for subsequent wave field simulation and attenuation compensation calculation, supporting the implementation of subsequent local time window calculation and attenuation compensation calculation, improving imaging accuracy, and laying the foundation for accurate imaging of deep targets in highly conductive media.

[0009] As a preferred embodiment of the ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation described in this invention, the step of calculating the corresponding local time window includes: The principal reflection time corresponding to each imaging point within the imaging area is calculated using the velocity model. Based on the principal reflection time, the range of the local time window corresponding to each imaging point is determined, and the corresponding local time window is obtained.

[0010] The beneficial effects of this preferred technical solution are as follows: by adaptively determining the local time window for each imaging point, the interception of wavefield data is realized, which not only avoids the storage of invalid wavefield information and reduces the storage requirements of the entire wavefield, but also suppresses the interference of radar noise outside the local time window by selectively extracting effective signals.

[0011] As a preferred embodiment of the ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation described in this invention, the step of obtaining the sparse forward propagation wavefield includes: The full forward propagation wavefield is obtained by performing forward propagation wavefield simulation, and the wavefield data of each imaging point in the corresponding local time window is extracted through the full forward propagation wavefield. The wavefield data within the local time window is sparsely sampled using the Nyquist sampling theorem to obtain a sparse forward propagating wavefield.

[0012] The beneficial effects of this preferred technical solution are as follows: by combining the local time window and the Nyquist sampling law, the forward propagation wavefield is structurally compressed, which reduces the amount of data stored in the wavefield while retaining the sampling points of the effective imaging signal, thereby significantly reducing storage overhead and input / output burden, and providing technical support for achieving efficient reverse time migration imaging.

[0013] As a preferred embodiment of the ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation described in this invention, the step of obtaining the compensated backpropagating wavefield includes: Based on the wave equation controlling the propagation of electromagnetic waves, perform a reverse propagation wave field simulation; In the reverse propagation wave field simulation process, the conductivity term in the wave equation is inverted, and attenuation compensation calculation is performed. The compensated reverse propagation wavefield is obtained through the attenuation compensation calculation.

[0014] The beneficial effects of this preferred technical solution are as follows: by reversing the sign of the conductivity term in the wave equation in the reverse wave field simulation, effective compensation for the attenuation of electromagnetic waves propagating in conductive media is achieved, which can restore the deep signal energy lost due to medium absorption, improve the imaging intensity of reflected signals in high-attenuation strata, and enhance the imaging signal-to-noise ratio and clarity of deep targets.

[0015] As a preferred embodiment of the ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation described in this invention, the step of obtaining the imaging value of each imaging point includes: By multiplying and accumulating the sparse forward propagation wavefield and the compensated reverse propagation wavefield at each imaging point within the corresponding local time window, and performing cross-correlation imaging condition calculation, the imaging value of each imaging point is obtained.

[0016] The beneficial effects of this preferred technical solution are as follows: by performing cross-correlation imaging conditions between the sparse forward propagation wavefield and the compensated reverse propagation wavefield within a local time window, the computational efficiency and imaging quality are unified. It not only inherits the storage advantages brought by sparse sampling, but also enhances the imaging energy of weak signals in the deep area through attenuation compensation operation. Ultimately, it improves the recognition ability of deep targets and the image signal-to-noise ratio while ensuring imaging accuracy.

[0017] As a preferred embodiment of the ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation described in this invention, the step of generating the reverse-time migration imaging profile includes: The imaging values ​​of all imaging points within the imaging area are combined in a grid to generate a reverse-time offset imaging profile.

[0018] The beneficial effects of this preferred technical solution are as follows: by combining the imaging values ​​of all imaging points in the imaging area into a grid, complete and continuous imaging of the underground structure is achieved. This combines the technical advantages of local wavefield sparse sampling and attenuation compensation calculation in the previous steps, generating a high-resolution, deep signal-enhanced reverse time-shifted imaging profile, providing intuitive and reliable image evidence for geological interpretation and target identification.

[0019] Secondly, the present invention provides a ground-penetrating radar reverse time migration imaging system based on local Nyquist attenuation compensation, comprising: The data acquisition and modeling module is used to acquire ground-penetrating radar observation data and build velocity models; The local time window calculation module is used to determine the imaging region and calculate the local time window for each imaging point using the velocity model and observation geometry. The forward propagation simulation module is used to perform forward propagation wavefield simulation and generate a sparse forward propagation wavefield based on the local time window and Nyquist sampling law. The anti-propagation simulation module is used to perform anti-propagation wave field simulation and generate a compensated anti-propagation wave field by reversing the sign of the conductivity term in the wave equation. The imaging calculation module is used to perform cross-correlation calculation on the sparse forward propagation wavefield and the compensated reverse propagation wavefield within a local time window to obtain the imaging value of each imaging point. The profile generation module is used to integrate all imaging values ​​and generate a reverse time-shifted imaging profile.

[0020] Thirdly, the present invention provides an electronic device, comprising: Memory and processor; The memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions, which, when executed by the processor, implement the steps of the ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation.

[0021] Fourthly, the present invention provides a computer-readable storage medium storing computer-executable instructions that, when executed by a processor, implement the steps of the ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation.

[0022] Compared with existing technologies, the beneficial effects of this invention are as follows: through the synergy of three mechanisms—local time window, Nyquist sparse sampling, and attenuation compensation operation—it reduces wavefield storage requirements and computational burden while maintaining imaging accuracy; it compensates for the energy attenuation of electromagnetic waves in conductive media, enhances the imaging energy and signal-to-noise ratio of deep reflection signals, solves the technical problems of high storage cost and insufficient deep imaging capability of traditional methods, and provides an efficient and reliable solution for high-precision ground-penetrating radar imaging under complex geological conditions. Attached Figure Description

[0023] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. 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.

[0024] Figure 1 This is a schematic diagram of the overall process of a ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation according to an embodiment of the present invention.

[0025] Figure 2 The diagram shows the distribution of dielectric constant and conductivity in the four-layer model.

[0026] Figure 3 This is the forward simulation profile and preprocessing results of the four-layer model.

[0027] Figure 4 The results are reverse time migrations under different imaging conditions. Detailed Implementation

[0028] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the protection scope of the present invention.

[0029] Example 1, referring to Figures 1-4 As an embodiment of the present invention, a ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation is provided, comprising: S100: Collect ground-penetrating radar observation data and construct a velocity model.

[0030] S200. Based on the velocity model and the observation geometry of the ground-penetrating radar observation data, the imaging region is obtained, and a corresponding local time window is calculated for each imaging point within the imaging region.

[0031] S300. Perform forward propagation wave field simulation according to the wave equation to obtain the full forward propagation wave field. Through the full forward propagation wave field and the local time window, perform sparse sampling and storage according to the Nyquist sampling law to obtain the sparse forward propagation wave field.

[0032] S400. Perform anti-propagation wave field simulation according to the wave equation, and reverse the sign of the conductivity term in the wave equation to achieve attenuation compensation, thereby obtaining the compensated anti-propagation wave field.

[0033] S500. The sparse forward propagation wavefield and the compensated reverse propagation wavefield are cross-correlated within the local time window to obtain the imaging values ​​of each imaging point. The imaging values ​​of all imaging points are integrated to generate a reverse time-shifted imaging profile.

[0034] It should be noted that the traditional full-storage ground-penetrating radar reverse time migration method requires the preservation of full-wave field data under the traditional cross-correlation imaging conditions, resulting in huge storage overhead and limiting its application in large-scale two-dimensional and three-dimensional data. In addition, electromagnetic waves will experience energy attenuation when propagating in highly conductive media, resulting in weak deep reflection signals. Existing sparse sampling methods have insufficient imaging stability under attenuation environments and are difficult to effectively recover effective deep signals.

[0035] Therefore, to address the issues of high storage costs and insufficient deep imaging capabilities, the following steps (S100-S500) are employed: first, a velocity model incorporating dielectric constant and conductivity is constructed; then, the imaging region is determined based on the observation geometry of the ground-penetrating radar (GPR) data, and a local time window is calculated; a sparse forward propagation wavefield is obtained through forward propagation wavefield simulation combined with Nyquist sparse sampling; in the reverse propagation simulation, attenuation compensation is achieved by reversing the sign of the conductivity term; finally, the sparse forward propagation wavefield and the compensated reverse propagation wavefield are cross-correlated and imaged within the local time window, achieving an effective balance between storage efficiency and deep imaging quality, and providing a reliable technical solution for high-precision GPR detection under complex geological conditions.

[0036] Example 2, refer to Figures 1-4 As an embodiment of the present invention, based on the above embodiment, a ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation is provided.

[0037] In this embodiment of the application, taking urban underground space exploration as the application scenario, the steps of constructing a velocity model based on the ground-penetrating radar observation data collected in S100 (A1~A2) are specifically implemented as follows: A1. Collect ground-penetrating radar observation data.

[0038] Specifically, a radar antenna with a center frequency of 900MHz is used to collect ground-penetrating radar observation data along the survey line at fixed track spacing, recording the raw radar signal containing radar direct waves, radar reflected waves, and radar noise.

[0039] A2. Based on the ground-penetrating radar observation data, construct a velocity model of dielectric constant and conductivity.

[0040] Specifically, a velocity model is constructed using the dielectric constant inversion method. This is achieved through travel time analysis of radar direct waves or reflectors of known depth, using the formula: in, The speed of electromagnetic wave propagation. The speed of light in a vacuum Given the relative permittivity, the permittivity distribution is obtained by inversion. .

[0041] conductivity model The initial conductivity of each layer is assigned based on the conductivity parameter values ​​in a typical medium. Utilizing the amplitude attenuation characteristics of the original radar signal in the medium, the energy attenuation curves of radar reflected waves at different depths are fitted, and the attenuation coefficient is used as the basis for the calculation. With conductivity Relationship: in, The attenuation coefficient is... For electrical conductivity, Angular frequency, Permeability, The absolute permittivity of the medium. Given the vacuum permittivity, the conductivity model is constrained and optimized to obtain a value containing... and The velocity model of the parameters.

[0042] In an optional implementation, step S100 may further employ a metal reflection calibration method to construct a velocity model. The steps are as follows: place a metal reflector at a known depth, collect ground-penetrating radar observation data, measure the two-way travel time of the radar reflected wave, and then apply the formula... Calculate the dielectric constant at the known depth, where, For two-way travel, The depth is used as the control point for spatial interpolation or constraint inversion to construct a velocity model.

[0043] In another optional implementation, step S100 can also use common center point (CMP) velocity analysis to construct a velocity model. The steps are as follows: acquire CMP gathers, scan and fit the hyperbolic radar reflection waves at different offsets, convert the optimal root mean square velocity obtained by fitting into layer velocity through the Dix formula, and then obtain the dielectric constant distribution to construct a velocity model.

[0044] In this embodiment of the application, step S200, which involves calculating the corresponding local time window, includes steps B1 to B2: B1. Calculate the principal reflection time corresponding to each imaging point within the imaging area using the velocity model.

[0045] Specifically, the imaging region refers to the two-dimensional spatial range jointly determined by the velocity model and the observation geometry of ground-penetrating radar (GPR) data. The horizontal range is determined by the start and end points of the radar survey lines, and the depth range is determined by the set maximum imaging depth. The imaging region is discretized into regular grid points (i.e., imaging points), based on the imaging points in the velocity model. electromagnetic wave propagation speed The principal reflection time is calculated using the formula: in, Principal reflection time, This represents the distance from the imaging point to the transmitting antenna; calculations are performed for each imaging point. The wave propagation time to the transmitting antenna is taken as the principal reflection time at the imaging point. .

[0046] B2. Based on the principal reflection time, determine the range of the local time window corresponding to each imaging point to obtain the corresponding local time window.

[0047] Specifically, based on the principal reflection time and source wavelet width The local time window range for each imaging point is determined using the following formula: in, For local time windows, To account for the time tolerance of uncertainties in the velocity model, a value is taken as... ; The source wavelet width is 2ns, with the time width of the 900MHz Ricker wavelet as the reference value.

[0048] In an optional implementation, step S200 may further involve determining a local time window based on the maximum and minimum speeds, the step being: based on the maximum speed in the speed model... and minimum speed Through the formula: Determine the range of the local time window.

[0049] In another optional implementation, step S200 may also employ the method of determining a local time window based on ray tracing, which involves: calculating the main reflection path and two-way travel time of each imaging point using the ray tracing method, and then determining the local time window based on the two-way travel time and wavelet width.

[0050] In this embodiment of the application, step S300, the step of obtaining the sparse forward propagation wave field, includes C1~C2: C1. Obtain the full forward propagation wavefield by performing forward propagation wavefield simulation, and extract the wavefield data of each imaging point within the corresponding local time window using the full forward propagation wavefield.

[0051] Specifically, in the full forward propagation wavefield obtained from the forward propagation wavefield simulation, for each imaging point... Based on the local time window corresponding to that point Extract the wavefield time series within this local time window. ,in .

[0052] C2. The wavefield data within the local time window is sparsely sampled using the Nyquist sampling theorem to obtain a sparse forward propagation wavefield.

[0053] Specifically, the Nyquist sampling interval is determined based on the maximum and minimum dielectric constants of the medium, and the sparse sampling interval is calculated. The formula is as follows: in, For sparse sampling interval, To meet the time step requirements for numerical stability, =0.01ns, The compression ratio is given by the formula: in, For maximum speed, For minimum speed, It is the maximum relative permittivity of the medium. Let be the minimum relative permittivity of the medium, and k be a safety factor greater than 1, used to ensure that the sampling law is still satisfied under complex wave field conditions. Let k = 5. , Calculated .

[0054] Therefore, the sampling interval The extracted wavefield data within the local time window are processed according to... By resampling at intervals, the wavefield data is compressed to 1 / β of its original value, resulting in a sparse forward propagation wavefield.

[0055] In an optional implementation, step S300 may also employ a sparse sampling method with a fixed compression ratio, comprising the following steps: setting a fixed compression ratio β=2, and sampling according to the corresponding sampling interval. Sparse sampling of wavefield data within a local time window is suitable for scenarios with high computational efficiency requirements.

[0056] In another optional implementation, step S300 may also employ an adaptive sparse sampling method, the steps of which are: based on the conductivity distribution in the velocity model... Adjust compression ratio With conductivity As a boundary, when During the high attenuation region, a smaller compression ratio is used. When σ≤0.005S / m is the low attenuation region, a larger compression ratio is used. This further improves compression efficiency.

[0057] In this embodiment of the application, step S400, the step of obtaining the compensated reverse propagation wavefield, includes D1~D3: D1. Perform a reverse propagation wave field simulation based on the wave equation that controls the propagation of electromagnetic waves.

[0058] Specifically, the wave equation adopts the two-dimensional TM mode Maxwell equation. Based on the two-dimensional TM mode Maxwell equation, the backpropagating wave field is solved using the finite difference time-domain method. The computational domain is set to 2.0m×2.0m, the grid spacing satisfies the Nyquist sampling law, the time step satisfies the Courant stability condition, and the boundary condition adopts a perfectly matched layer absorbing boundary.

[0059] D2. During the reverse propagation wave field simulation, the conductivity term in the wave equation is inverted, and attenuation compensation is performed.

[0060] Specifically, the physical basis of attenuation compensation calculation lies in the propagation characteristics of electromagnetic waves in conductive media. According to electromagnetic wave field theory, complex velocity... It can be represented as: From this, the phase velocity can be derived. With attenuation coefficient : in, For complex velocity, For phase velocity, For attenuation coefficient and conductivity Positive correlation. The attenuation term causes the wave field energy to decrease exponentially with propagation distance. When solving the wave equation controlling electromagnetic wave propagation, the conductivity parameter σ is replaced with -σ. For Maxwell's equations in the two-dimensional TM mode: The conductivity term in the third equation Modified to This achieves attenuation compensation.

[0061] D3. The compensated reverse propagation wave field is obtained through the attenuation compensation calculation.

[0062] Specifically, through attenuation compensation calculations, the energy of the backpropagating wavefield is enhanced during propagation, especially the energy of the deep reflected signal passing through a high-conductivity dielectric layer is recovered, resulting in a compensated backpropagating wavefield. .

[0063] In an optional implementation, step S400 may further employ a frequency-band attenuation compensation method, comprising the steps of: decomposing the backpropagation wavefield into multiple frequency bands, and for each frequency band... Assign a compensation coefficient Wherein, the compensation coefficient The range of values ​​and frequency band Negative correlation, compensation coefficient corresponding to high frequency band satisfy The compensation coefficient corresponding to the low frequency band satisfy .

[0064] In another optional implementation, step S400 may also employ an adaptive compensation method, the steps of which are: based on the conductivity distribution in the velocity model... Adjust the compensation intensity: in the high conductivity region σ>0.008S / m, use full compensation, such as replacing σ with -σ; in the low conductivity region σ<0.008S / m, use partial compensation, such as replacing σ with -0.5σ, to balance the compensation effect and numerical stability.

[0065] In this embodiment of the application, step S500, the step of generating the reverse time-lapse imaging profile, includes E1~E2: E1. By multiplying and accumulating the sparse forward propagation wavefield and the compensated reverse propagation wavefield at each imaging point within the corresponding local time window, and performing cross-correlation imaging condition calculation, the imaging value of each imaging point is obtained.

[0066] Specifically, for each imaging point , its sparse forward propagation wave field With compensation reverse propagation wave field In the corresponding local time window The imaging values ​​are obtained by performing cross-correlation calculations within the internal system. The formula is as follows: in, For the image values, the summation operation is performed within a local time window. It will be carried out internally.

[0067] E2. Combine the imaging values ​​of all imaging points within the imaging area into a grid to generate a reverse-time offset imaging profile.

[0068] Specifically, the calculated imaging values ​​I(x,z) are arranged in a grid according to their coordinate positions in the calculation area to form a two-dimensional imaging data matrix. For example, in a calculation area of ​​2.0m×2.0m, a grid resolution of 0.01m×0.01m is used to interpolate the imaging values ​​onto regular grid points to generate a reverse time-shifted imaging profile.

[0069] In an optional implementation, step S500 may also employ a cross-correlation imaging method based on amplitude normalization. The steps are as follows: before cross-correlation calculation, amplitude normalization processing is performed on the sparse forward propagation wavefield and the compensated reverse propagation wavefield, using the following formula: It eliminates the influence of amplitude differences on imaging results and is suitable for imaging complex media with drastic amplitude changes.

[0070] In another optional implementation, step S500 can also employ a multi-scale fusion profile generation method, which involves generating reverse time migration imaging profiles with different frequency components, and then fusing these reverse time migration imaging profiles through wavelet transform or other multi-scale analysis methods. This preserves the detailed information of the high-frequency components while utilizing the stability of the low-frequency components to generate a final imaging profile of higher quality.

[0071] In summary, this invention provides a physical parameter basis for wavefield simulation and attenuation compensation by constructing a velocity model that includes dielectric constant and conductivity; it reduces the storage of invalid information by using an adaptive local time window to extract wavefield data; it reduces storage overhead and input / output burden by combining the Nyquist sampling law to perform structured sparse sampling of the local wavefield; it recovers the energy of deep reflection signals by reversing the sign of the conductivity term in the backpropagation wavefield simulation; and it achieves efficient and high-precision imaging through cross-correlation imaging conditions within the local time window. This invention provides an innovative solution to the technical problems of high storage cost and insufficient deep imaging capability of traditional reverse time migration methods, and provides reliable technical support for fields such as urban underground space exploration, geological exploration, and infrastructure safety assessment.

[0072] Example 3 illustrates a schematic scheme of a ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation. It should be noted that the technical solution of this ground-penetrating radar reverse-time migration imaging system based on local Nyquist attenuation compensation is based on the same concept as the aforementioned ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation. Details not described in detail in this embodiment can be found in the description of the aforementioned ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation.

[0073] This embodiment also provides a ground-penetrating radar reverse time migration imaging system based on local Nyquist attenuation compensation, including: The data acquisition and modeling module is used to acquire ground-penetrating radar observation data and build velocity models; The local time window calculation module is used to determine the imaging region and calculate the local time window for each imaging point using the velocity model and observation geometry. The forward propagation simulation module is used to perform forward propagation wavefield simulation and generate a sparse forward propagation wavefield based on the local time window and Nyquist sampling law. The anti-propagation simulation module is used to perform anti-propagation wave field simulation and generate a compensated anti-propagation wave field by reversing the sign of the conductivity term in the wave equation. The imaging calculation module is used to perform cross-correlation calculation on the sparse forward propagation wavefield and the compensated reverse propagation wavefield within a local time window to obtain the imaging value of each imaging point. The profile generation module is used to integrate all imaging values ​​and generate a reverse time-shifted imaging profile.

[0074] This embodiment also provides an electronic device applicable to ground-penetrating radar reverse time migration imaging based on local Nyquist attenuation compensation, comprising: a memory and a processor; the memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions to implement the ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation as proposed in the above embodiment.

[0075] This embodiment also provides a storage medium storing a computer program that, when executed by a processor, implements the ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation as proposed in the above embodiments.

[0076] The storage medium proposed in this embodiment and the ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation proposed in the above embodiments belong to the same inventive concept. Technical details not described in detail in this embodiment can be found in the above embodiments, and this embodiment has the same beneficial effects as the above embodiments.

[0077] Based on the above description of the implementation methods, those skilled in the art can clearly understand that the present invention can be implemented using software and necessary general-purpose hardware, and of course, it can also be implemented using hardware. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as a computer floppy disk, read-only memory (ROM), random access memory (RAM), flash memory, hard disk, or optical disk, etc., including several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods of the various embodiments of the present invention.

[0078] Example 4: Based on the previous three examples, this example provides a ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation. In order to verify the beneficial effects of the present invention, scientific demonstration is carried out through economic benefit calculation and simulation experiment.

[0079] To verify the effectiveness and stability of the local Nyquist attenuation compensation ground-penetrating radar reverse time migration imaging proposed in this paper, two sets of numerical experiments are designed in this section: (1) forward modeling and imaging experiments based on a four-layer medium model, used to analyze the energy recovery and structure recognition capabilities of various RTM imaging conditions under noise-free conditions; (2) noise resistance performance analysis under different noise levels, used to evaluate the imaging stability of the algorithm in complex noise environments.

[0080] All experiments were based on Maxwell's equations in two-dimensional TM mode, solved using the finite-difference time-domain (FDTD) method. The grid spacing and time step size of the computational domain satisfied the Nyquist sampling theorem and the Courant stability condition. A perfectly matched layer (PML) was used as the absorbing boundary. The computational domain was 2.0 m × 2.0 m in size, the source signal was a 900 MHz Ricker wavelet, and the time window length was set to 30 ns, corresponding to 1800 time steps.

[0081] Figure 2 (a)~(b) show the distribution of dielectric constant and conductivity of the four-layer model. There are differences in dielectric constant and conductivity between the layers to simulate typical non-uniform underground media. The parameters of each layer are shown in Table 1.

[0082] Table 1 Medium Parameter Table

[0083] In this four-layer model, the third layer is a high-conductivity medium, which leads to energy attenuation and waveform distortion, making it a critical region for verifying the effectiveness of attenuation compensation. To examine the algorithm's adaptability to velocity model errors, the velocity model used for migration was Gaussian smoothed (e.g., ...). Figure 2 (c)~(d) are shown to simulate the velocity uncertainty commonly found in actual GPR imaging.

[0084] Figure 3 This is the forward simulation profile and preprocessing results of the four-layer model. Figure 4 The reverse time migration results under different imaging conditions are shown. The overall imaging effects of CIC, NCIC, and LCIC ((a)~(c)) are basically consistent, indicating that moderately sparse sampling does not lose the main reflection information. EAIC imaging profile ( Figure 3 (d) The appearance of an arc-shaped artifact above the interface indicates that imaging conditions relying solely on single-moment amplitude are insufficient to reliably extract effective reflections. In contrast, ACCIC and NLACCIC demonstrate superior performance in deep energy recovery: through attenuation compensation, both enhance the energy of reflected waves passing through the high-conductivity layer, resulting in clearer imaging of deep interfaces. The imaging results of NLACCIC and ACCIC are almost identical, but NLACCIC only stores local wavefield data satisfying the Nyquist condition, reducing storage requirements by approximately 95%.

[0085] Table 2 summarizes the comparative analysis of storage cost, running time, NRMSE, PSNR and DRR of different algorithms, as shown in Table 2.

[0086] Storage cost and runtime were compared with conventional full-storage CIC. The results show that while EAIC has the lowest storage cost (0.31 MB), its runtime is slightly longer (4.22 min), mainly because it needs to search for the maximum energy density point and perform I / O operations at each time step, increasing computational latency. NLACCIC has the shortest runtime (3.72 min) and a storage cost of 23.50 MB, only slightly higher than EAIC. Compared to CIC, NLACCIC saves 95.72% of storage space and 11.70% of runtime. Furthermore, NLACCIC has a NRMSE of 0.0176, a PSNR of 12.6588, and a DRR of 29.9661, demonstrating superior overall performance. This indicates that the proposed local Nyquist sparse sampling strategy can reduce wavefield preservation without incurring additional computational burden.

[0087] Table 2 Comparison of storage cost, runtime, and NRMSE for different RTM algorithms

[0088] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation, characterized in that, include: Collect ground-penetrating radar observation data and construct a velocity model; The imaging region is obtained by combining the velocity model with the observation geometry of the ground-penetrating radar observation data, and a corresponding local time window is calculated for each imaging point within the imaging region. A forward propagation wavefield simulation is performed based on the wave equation to obtain a full forward propagation wavefield. The full forward propagation wavefield and the local time window are then sparsely sampled and stored according to the Nyquist sampling law to obtain a sparse forward propagation wavefield. Based on the wave equation, a reverse propagation wave field simulation is performed, and the conductivity term in the wave equation is reversed to achieve attenuation compensation, thus obtaining a compensated reverse propagation wave field. The sparse forward propagation wavefield and the compensated reverse propagation wavefield are cross-correlated within the local time window to obtain the imaging values ​​of each imaging point. The imaging values ​​of all imaging points are then integrated to generate a reverse time-shifted imaging profile.

2. The ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation as described in claim 1, characterized in that, The steps for collecting ground-penetrating radar observation data and constructing a velocity model include: Collect ground-penetrating radar observation data; Based on the ground-penetrating radar observation data, a velocity model of dielectric constant and conductivity is constructed.

3. The ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation as described in claim 2, characterized in that, The steps to calculate the corresponding local time window include: The principal reflection time corresponding to each imaging point within the imaging area is calculated using the velocity model. Based on the principal reflection time, the range of the local time window corresponding to each imaging point is determined, and the corresponding local time window is obtained.

4. The ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation as described in claim 3, characterized in that, The steps to obtain a sparse forward propagation wave field include: The full forward propagation wavefield is obtained by performing forward propagation wavefield simulation, and the wavefield data of each imaging point in the corresponding local time window is extracted through the full forward propagation wavefield. The wavefield data within the local time window is sparsely sampled using the Nyquist sampling theorem to obtain a sparse forward propagating wavefield.

5. The ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation as described in claim 4, characterized in that, The steps to obtain the compensated backpropagation wave field include: Based on the wave equation controlling the propagation of electromagnetic waves, perform a reverse propagation wave field simulation; In the reverse propagation wave field simulation process, the conductivity term in the wave equation is inverted, and attenuation compensation calculation is performed. The compensated reverse propagation wavefield is obtained through the attenuation compensation calculation.

6. The ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation as described in claim 5, characterized in that, The steps to obtain the imaging value for each imaging point include: By multiplying and accumulating the sparse forward propagation wavefield and the compensated reverse propagation wavefield at each imaging point within the corresponding local time window, and performing cross-correlation imaging condition calculation, the imaging value of each imaging point is obtained.

7. The ground-penetrating radar reverse-time migration imaging method based on local Nyquist attenuation compensation as described in claim 6, characterized in that, The steps for generating a reverse time-lapse imaging profile include: The imaging values ​​of all imaging points within the imaging area are combined in a grid to generate a reverse-time offset imaging profile.

8. A ground-penetrating radar reverse-time migration imaging system based on local Nyquist attenuation compensation, employing the method described in any one of claims 1-7, characterized in that, include: The data acquisition and modeling module is used to acquire ground-penetrating radar observation data and build velocity models; The local time window calculation module is used to determine the imaging region and calculate the local time window for each imaging point using the velocity model and observation geometry. The forward propagation simulation module is used to perform forward propagation wavefield simulation and generate a sparse forward propagation wavefield based on the local time window and Nyquist sampling law. The anti-propagation simulation module is used to perform anti-propagation wave field simulation and generate a compensated anti-propagation wave field by reversing the sign of the conductivity term in the wave equation. The imaging calculation module is used to perform cross-correlation calculation on the sparse forward propagation wavefield and the compensated reverse propagation wavefield within a local time window to obtain the imaging value of each imaging point. The profile generation module is used to integrate all imaging values ​​and generate a reverse time-shifted imaging profile.

9. An electronic device, comprising: Memory and processor; The memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions. When the computer-executable instructions are executed by the processor, they implement the steps of the ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation as described in any one of claims 1 to 7.

10. A computer-readable storage medium storing computer-executable instructions that, when executed by a processor, implement the steps of the ground-penetrating radar reverse time migration imaging method based on local Nyquist attenuation compensation as described in any one of claims 1 to 7.

Images