Ground penetrating radar concurrent backprojection imaging method and device for multi-layer media

Abstract

The present invention provides a method and device for concurrent back-projection imaging of a multi-layered ground-penetrating radar. The method comprises: setting an imaging region according to the detection region and discretizing it to obtain a three-dimensional imaging matrix; dividing the three-dimensional imaging matrix into a plurality of sub-imaging matrices equal to the number of layers of the detection region medium according to the distribution of the detection region medium; creating a thread pool according to the number of layers of the detection region medium, wherein the number of core threads and the maximum number of threads of the thread pool are equal to the number of layers of the detection region medium; having each thread in the thread pool concurrently perform back-projection imaging operations on the corresponding sub-imaging matrices according to the echo data of the ground-penetrating radar, and obtaining the imaging value of each imaging unit in each sub-imaging matrix; normalizing and imaging the three-dimensional imaging matrix to obtain a concurrent back-projection imaging result. Through the above methods, the present invention can reduce the computational complexity of the electromagnetic wave propagation path in the multi-layered medium background in the traditional back-projection imaging algorithm, and significantly shorten the computational time of three-dimensional imaging.

Description

Technical Field

[0001] The present invention belongs to the technical field of ground penetrating radar imaging, and in particular relates to a ground penetrating radar concurrent backprojection imaging method and device for multi-layer media. Background Art

[0002] Ground Penetrating Radar (GPR) is a device that uses an antenna to transmit and receive high-frequency electromagnetic waves to detect the properties and distribution of materials in underground space. During propagation, electromagnetic waves scatter and refract when encountering objects or structures with electromagnetic properties different from those of the background medium. Based on this characteristic, GPR can effectively estimate key parameters such as the spatial position, geometric structure, and material properties of underground targets by analyzing echo data, enabling non-destructive detection of subsurface targets. At construction sites, GPR can be used to inspect wall structures, detect defects such as cracks and voids in concrete, and locate pipelines and cables, thereby preventing and minimizing damage to underground facilities. In archaeological exploration, GPR can be used to detect buried artifacts and historical sites. In geological exploration and mining, GPR can be used to detect ore bodies, strata, and rock structures, helping explorers better understand the subsurface and improving the accuracy of mining exploration. GPR also has broad application prospects in groundwater detection, highway and bridge inspection, and glacier detection.

[0003] In the field of ground-penetrating radar (GPR) data analysis, the "delay-sum" backprojection method is widely used in GPR signal processing due to its effectiveness in imaging point-scattering targets in complex, multi-layered media. Prior art backprojection imaging of GPR in multi-layered media requires traversing the imaging area and iteratively solving all imaging units to obtain imaging values. To solve the imaging value for each imaging unit, it is necessary to traverse all aperture points and calculate the propagation path of the electromagnetic wave between the aperture point and the imaging unit.

[0004] Now assume that the aperture point is 0m away from the ground, the line dimension coordinate is C0, the imaging unit is located in the nth layer of the medium, and the line dimension coordinate is C n , the relative dielectric constant of the i-th layer medium is ε i , thickness is The intersection point of the electromagnetic wave propagation path and the upper surface of the i-th layer of medium is C i-1 , and the intersection point with the lower surface of the i-th layer of medium is C i , assuming that the intersection point C1 of the electromagnetic wave propagation path and the lower surface of the first layer of medium is the unknown number x, then based on Snell's law, we can get:

[0005]

[0006] After solving for x, Snell's law can be used to sequentially determine the line-dimensional coordinates of the intersections of the electromagnetic wave propagation path and each medium interface, thereby determining the electromagnetic wave propagation path and time delay. Based on the calculated time delay, the signal amplitude corresponding to the A-scan echo data at each aperture point is extracted. All obtained signal amplitudes are combined into a one-dimensional imaging signal sequence for that imaging unit, and the sum of these sequences is the imaging value for that imaging unit. The imaging operation is performed on all imaging units in the imaging area, and then normalized to obtain the backprojection imaging result.

[0007] It can be seen that in the prior art, when performing backprojection imaging calculations on ground penetrating radar in multi-layer media, there are problems of poor imaging quality and long calculation time. Summary of the Invention

[0008] The present invention provides a concurrent back-projection imaging method and device for a ground-penetrating radar in a multi-layer medium, which solves the problems of poor imaging quality and long calculation time in the back-projection imaging operation of the ground-penetrating radar in a multi-layer medium.

[0009] Based on the above objectives, the present invention proposes a concurrent back-projection imaging method for a ground-penetrating radar of a multi-layer medium, comprising: setting an imaging area according to a detection area and discretizing it to obtain a three-dimensional imaging matrix; dividing the three-dimensional imaging matrix into a plurality of sub-imaging matrices equal to the number of medium layers in the detection area according to the medium distribution in the detection area, wherein one sub-imaging matrix corresponds to one layer of medium; creating a thread pool according to the number of medium layers in the detection area, wherein the number of core threads and the maximum number of threads in the thread pool are equal to the number of medium layers in the detection area; having each thread in the thread pool concurrently perform a back-projection imaging operation on the corresponding sub-imaging matrix according to the echo data of the ground-penetrating radar, and obtaining an imaging value of each imaging unit in each of the sub-imaging matrices; normalizing and imaging the three-dimensional imaging matrix to obtain a concurrent back-projection imaging result.

[0010] Optionally, the method of concurrently performing back-projection imaging operations on corresponding sub-imaging matrices according to the echo data of the ground penetrating radar by each thread in the thread pool to obtain imaging values ​​of each imaging unit in each sub-imaging matrix includes: for any thread, for any aperture point, creating an imaging unit queue, and placing any imaging unit in the sub-imaging matrix as an initial imaging unit into the imaging unit queue; sequentially taking out imaging units from the imaging unit queue to calculate an estimated electromagnetic wave propagation path between the aperture point and the imaging unit; calculating a two-way electromagnetic wave propagation delay between the aperture point and each imaging unit according to the estimated electromagnetic wave propagation path and an electromagnetic wave velocity formula; extracting a signal amplitude of the A-scan echo data corresponding to the aperture point according to the two-way electromagnetic wave propagation delay, and accumulating it to the imaging value of the corresponding imaging unit; adding imaging units adjacent to the imaging unit that have not undergone iterative operations to the imaging unit queue, and repeatedly iteratively calculating the imaging values ​​of each imaging unit until the imaging unit queue is empty; and traversing each aperture point to obtain a final imaging value of each imaging unit in the corresponding sub-imaging matrix.

[0011] Optionally, the method of obtaining the estimated electromagnetic wave propagation path between the aperture point and the imaging unit includes: setting an electromagnetic wave propagation profile based on the aperture point and the imaging unit and establishing a two-dimensional rectangular coordinate system; obtaining the initial values ​​of the intersection sequence of the electromagnetic wave propagation path and the interface of each layer of the medium through an initialization method; and performing iterative optimization and solution based on a greedy idea to obtain the estimated electromagnetic wave propagation path between the aperture point and the corresponding imaging units in the sub-imaging matrix.

[0012] Optionally, the initialization method is used to obtain the initial value of the sequence of intersections between the electromagnetic wave propagation estimation path and the interfaces of each layer of media, including: traversing the adjacent imaging units in the upper, lower, left, right, front and rear directions of the imaging unit; if any of the adjacent imaging units has completed the solution of the electromagnetic wave propagation estimation path with the aperture point, then the calculation result of any of the solved adjacent imaging units is taken as the initial value of the sequence of intersections between the electromagnetic wave propagation estimation path and the interfaces of each layer of media; otherwise, the intersection of the electromagnetic wave propagation estimation path and each media interface is initialized to the intersection of the line from the aperture point to the imaging unit and each media interface, to obtain the initial value of the intersection sequence.

[0013] Optionally, the iterative optimization solution is performed based on the greedy idea to obtain the estimated electromagnetic wave propagation path between the aperture point and the corresponding imaging units in the sub-imaging matrix, including: calculating the deviation of each intersection in the intersection sequence and the total deviation of the intersection sequence according to the electromagnetic wave propagation estimated intersection deviation evaluation formula; setting the estimation interval according to each intersection in the intersection sequence and the deviation of each intersection; iteratively updating the intersection with the largest absolute value of deviation in the intersection sequence to the midpoint of the estimation interval, updating the intersection sequence, and calculating the total deviation of the updated intersection sequence; if the updated total deviation is greater than the total deviation before the update, then modifying the intersection with the largest absolute value of deviation in the upper and lower limits of the estimation interval to update it to the midpoint, until the updated total deviation is less than the total deviation before the update; updating the intersection with the largest absolute value of deviation to the midpoint, repeating the iteration, and stopping the iteration when the total deviation of the intersection sequence is less than the preset threshold.

[0014] Optionally, calculating the deviation of each intersection in the intersection sequence and the total deviation of the intersection sequence according to the electromagnetic wave propagation estimated intersection deviation evaluation formula includes: calculating the deviation of each intersection in the intersection sequence according to the electromagnetic wave propagation estimated intersection deviation evaluation formula:

[0015]

[0016] Among them, Δ i is the deviation between the estimated path of electromagnetic wave propagation and the intersection point of the i-th layer medium interface, C i 、C i-1 、C i+1 are the intersection points of the estimated electromagnetic wave propagation path and the interfaces of the i-th, i-1-th, and i+1-th layers in the intersection sequence, respectively. i , ε i+1 are the dielectric constants of the i-th layer and the i+1-th layer respectively, are the thicknesses of the i-th layer and the i+1-th layer of medium respectively; and the sum of the absolute values ​​of the deviations of the intersections in the intersection sequence is calculated to obtain the total deviation of the intersection sequence.

[0017] Optionally, the step of extracting the signal amplitude of the A-scan echo data corresponding to the aperture point according to the electromagnetic wave round-trip propagation delay and accumulating it to the imaging value of the corresponding imaging unit includes: jf The jth aperture point A j A-scan echo data j The signal amplitude in (t) Among them, round(·) means rounding to the nearest integer, Δt is the time window sampling interval, T is the time window length, W is the number of sampling points in the time dimension of the A-scan echo data; the signal amplitude The imaging value is accumulated to the imaging unit.

[0018] Based on the same inventive concept, the present invention also proposes a concurrent back-projection imaging device for a ground-penetrating radar of a multi-layer medium, comprising: an imaging matrix division unit, for setting an imaging area according to a detection area and discretizing it to obtain a three-dimensional imaging matrix, and dividing the three-dimensional imaging matrix into a plurality of sub-imaging matrices equal to the number of medium layers in the detection area according to the medium distribution in the detection area, wherein one sub-imaging matrix corresponds to one layer of medium; a thread pool creation unit, for creating a thread pool according to the number of medium layers in the detection area, wherein the number of core threads and the maximum number of threads in the thread pool are equal to the number of medium layers in the detection area; a back-projection imaging unit, for concurrently performing back-projection imaging operations on corresponding sub-imaging matrices according to the echo data of the ground-penetrating radar through each thread in the thread pool, and obtaining imaging values ​​of each imaging unit in each sub-imaging matrix; a normalization unit, for normalizing and imaging the three-dimensional imaging matrix, and obtaining concurrent back-projection imaging results.

[0019] Based on the same inventive concept, the present invention also proposes an electronic device, including a memory, a processor, and a computer program stored in the memory and runnable on the processor. When the processor executes the program, the ground-penetrating radar concurrent backprojection imaging method for multi-layer media as described above is implemented.

[0020] Based on the same inventive concept, the present invention also proposes a computer storage medium, wherein the storage medium stores at least one executable instruction, and the executable instruction enables a processor to execute the above-mentioned multi-layer medium ground penetrating radar concurrent backprojection imaging method.

[0021] As can be seen from the above, the beneficial effects of the technical solution provided by the present invention are: the present invention provides a method and device for concurrent backprojection imaging of a multi-layered medium by a ground-penetrating radar, the method comprising: setting an imaging area according to the detection area and discretizing it to obtain a three-dimensional imaging matrix; dividing the three-dimensional imaging matrix into a plurality of sub-imaging matrices equal to the number of medium layers in the detection area according to the medium distribution in the detection area, and one sub-imaging matrix corresponding to one layer of medium; creating a thread pool according to the number of medium layers in the detection area, the number of core threads and the maximum number of threads of the thread pool being equal to the number of medium layers in the detection area; using each thread in the thread pool to concurrently perform backprojection imaging operations on the corresponding sub-imaging matrices according to the echo data of the ground-penetrating radar, and obtaining imaging values ​​of each imaging unit in each sub-imaging matrix; normalizing and imaging the three-dimensional imaging matrix to obtain concurrent backprojection imaging results, which can reduce the computational complexity of the electromagnetic wave propagation path in the multi-layered medium background in the traditional backprojection imaging algorithm, realize concurrent backprojection imaging in three-dimensional space through the thread pool technology, and greatly shorten the calculation time of three-dimensional imaging. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0023] Figure 1 Schematic diagram of the flow of a method for concurrent backprojection imaging of a multi-layered medium by ground penetrating radar according to an embodiment of the present invention;

[0024] Figure 2 A schematic diagram of a method for performing a backprojection imaging operation on a corresponding sub-imaging matrix by any thread according to an embodiment of the present invention;

[0025] Figure 3 This is an example diagram of aperture point arrangement according to an embodiment of the present invention;

[0026] Figure 4 An example diagram of a three-dimensional model according to an embodiment of the present invention;

[0027] Figure 5 This is an example diagram of the electromagnetic wave propagation path during the back-projection imaging operation according to an embodiment of the present invention;

[0028] Figure 6 This is an example diagram of a three-dimensional back-projection imaging result obtained by back-projection imaging calculation according to an embodiment of the present invention;

[0029] Figure 7This is a schematic structural diagram of a multi-layer medium ground penetrating radar concurrent backprojection imaging device according to an embodiment of the present invention;

[0030] Figure 8 FIG. 4 is a schematic diagram of the hardware structure of an electronic device according to an embodiment of the present invention. DETAILED DESCRIPTION

[0031] In order to make the objectives, technical solutions and advantages of the present disclosure more clearly understood, the present disclosure is further described in detail below in conjunction with specific embodiments and with reference to the accompanying drawings.

[0032] It should be noted that, unless otherwise defined, the technical terms or scientific terms used in the embodiments of the present invention should have the usual meanings understood by people with ordinary skills in the field to which the present disclosure belongs. The "first", "second" and similar words used in the embodiments of the present invention do not indicate any order, quantity or importance, but are only used to distinguish different components. "Include" or "comprise" and similar words mean that the elements or objects appearing before the word include the elements or objects listed after the word and their equivalents, without excluding other elements or objects. "Connect" or "connected" and similar words are not limited to physical or mechanical connections, but may include electrical connections, whether direct or indirect. "Up", "down", "left", "right" and the like are only used to indicate relative positional relationships. When the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0033] The embodiment of the present invention implements a method for concurrent back-projection imaging of a multi-layered medium by ground penetrating radar. Figure 1 As shown, the concurrent backprojection ground penetrating radar imaging method for multi-layer media includes:

[0034] Step S11: setting an imaging area according to the detection area and discretizing it to obtain a three-dimensional imaging matrix. According to the medium distribution in the detection area, the three-dimensional imaging matrix is ​​divided into a plurality of sub-imaging matrices equal to the number of medium layers in the detection area, and one sub-imaging matrix corresponds to one layer of medium.

[0035] Ground-penetrating radar (GPR) has multiple data recording methods, the most basic of which is called A-scan echo data. A-scan echo data is single-channel echo data obtained by the fixed radar's transmitting and receiving antennas after completing a single detection. It records information such as the instantaneous amplitude, instantaneous phase, instantaneous frequency, and round-trip propagation time of the electromagnetic wave during its propagation through the detection area. B-scan echo data is a data format that combines multiple channels of A-scan echo data in a spatial aperture point sequence. Compared with A-scan echo data, B-scan echo data is richer in information and makes it easier to accurately and effectively analyze the detection area. C-scan echo data is a data format that combines the B-scan echo data corresponding to multiple equally spaced parallel survey lines on the survey plane in a spatial sequence.

[0036] In the embodiment of the present invention, C-scan echo data is acquired and pre-processed to obtain imaging parameters. The measurement line arrangement direction is denoted as x direction, there are K measurement lines in total, the measurement line direction is denoted as y direction, the measurement line range is [0, Y], there are N measurement points in total, and the coordinates of the i-th measurement point are y i , the A-scan echo data obtained at the i-th measuring point is recorded as The number of sampling points in the time dimension of the A-scan echo data is W, and the B-scan echo data are combined into C-scan echo data according to the spatial relationship. The C-scan echo data can be recorded as E0(x, y, t), which is a three-dimensional matrix with dimensions of K×N×W.

[0037] The operating parameters of the detection radar include the antenna step interval Δx, the time window length T, and the number of time dimension sampling points W. The time window sampling interval Δt is calculated using the following formula:

[0038]

[0039] The C-scan echo data is preprocessed to remove the direct coupled wave and obtain the preprocessed C-scan echo data. According to the detection area, the imaging area is set and discretized to obtain the three-dimensional imaging matrix O(x,y,z). According to the number of medium layers in the imaging area, O(x,y,z) is divided into n sub-imaging matrices O1(x,y,z),...,O n (x, y, z), where n is the number of medium layers in the imaging area. i (x, y, z) corresponds to the range of the i-th layer of medium in the imaging area.

[0040] Step S12: creating a thread pool according to the number of medium layers in the detection area, wherein the number of core threads and the maximum number of threads in the thread pool are equal to the number of medium layers in the detection area.

[0041] In the embodiment of the present invention, a thread pool is created by using the ThreadPoolExecutor open source tool under the concurrent open source framework in Python. The number of core threads and the maximum number of threads in the thread pool are both the number of media layers n in the imaging area, and the i-th thread t i Responsible for the i-th sub-imaging matrix O i Back projection imaging of (x,y,z).

[0042] Step S13: each thread in the thread pool concurrently performs a back-projection imaging operation on the corresponding sub-imaging matrix according to the echo data of the ground penetrating radar, and obtains the imaging value of each imaging unit in each sub-imaging matrix.

[0043] In the embodiment of the present invention, optionally, for any thread, the following method is used to perform back-projection imaging operation on the corresponding sub-imaging matrix to obtain the imaging value of each imaging unit in the corresponding sub-imaging matrix, such as Figure 2 Shown, including:

[0044] Step S131: For any aperture point, first create an imaging unit queue, and put any imaging unit in the sub-imaging matrix as an initial imaging unit into the imaging unit queue.

[0045] For example, when traversing the jth aperture point A j When creating an imaging unit queue Q j (C), and O i Any imaging unit C0=(x0,y0,z0) in (x,y,z) is put into the imaging unit queue Q j (C) in.

[0046] Step S132: sequentially taking out the imaging units in the imaging unit queue, and calculating the estimated electromagnetic wave propagation path between the aperture point and the imaging unit.

[0047] When electromagnetic waves propagate in multilayer media with different relative dielectric constants, they follow Snell's law, that is, electromagnetic waves propagate from a medium with a relative dielectric constant of ε to a medium with a relative dielectric constant of ε. i θ i The incident angle is the relative dielectric constant ε j The angle of refraction in the medium is θ j , and the relationship is satisfied:

[0048] In step S132, optionally, first, the electromagnetic wave propagation profile is set according to the aperture point and the imaging unit and a two-dimensional rectangular coordinate system is established. Then, the initialization method is used to obtain the initial values ​​of the intersection sequence of the electromagnetic wave propagation estimated path and each layer of the medium interface. For example, from the imaging unit queue Q j (C) Take out the imaging unit C at the head of the teamf , according to the aperture point A j With imaging unit C f The coordinates are set to define the electromagnetic wave propagation profile and establish a two-dimensional rectangular coordinate system. The initial values ​​of the intersection sequence I(C) between the estimated electromagnetic wave propagation path and each layer of the medium interface are obtained through active initialization or passive initialization methods. The length of the sequence is the number of medium interface layers m+1 between the aperture point and the imaging unit. Specifically, the thickness and relative dielectric constant of the aperture point, imaging unit and each layer of medium through which the electromagnetic wave propagates are traversed, and the adjacent imaging units in the directions above, below, left, right, front and back of the imaging unit are traversed; if any of the adjacent imaging units has completed the solution of the electromagnetic wave propagation estimation path with the aperture point, the calculation result of any of the adjacent imaging units that has been solved is taken as the initial value of the intersection sequence of the electromagnetic wave propagation estimation path and the interface of each layer of medium, that is, the initial value of the intersection sequence is obtained by the active initialization method; otherwise, the difference in the relative dielectric constant of each layer of medium is ignored, and the intersection of the electromagnetic wave propagation estimation path and each medium interface is initialized to the intersection of the line from the aperture point to the imaging unit and each medium interface, and the initial value of the intersection sequence is obtained, that is, the initial value of the intersection sequence is obtained by the passive initialization method. If the coordinates of the aperture point are (C0,0) and the coordinates of the imaging unit are (C m ,d m ), the depth of the interface of the i-th layer medium is d i , there are m layers of medium between the aperture point and the imaging unit, then through the passive initialization method, the intersection point C of the electromagnetic wave propagation estimation path and the interface of the i-th layer of medium i for:

[0049]

[0050] The sequence of intersection points between the estimated path of electromagnetic wave propagation and the interface of the medium is recorded as I(C)=[C0,...,C i ,...,C m ], the sequence length is m+1, where C0 is the aperture point coordinate, C m is the imaging unit coordinate.

[0051] After obtaining the initial value of the intersection point sequence I(C), an iterative optimization solution is performed based on a greedy approach to obtain the estimated electromagnetic wave propagation path between the aperture point and each imaging unit in the corresponding sub-imaging matrix. Specifically, a local optimal and global optimal solution method is designed based on the greedy approach, and I(C) is iteratively optimized to obtain the estimated electromagnetic wave propagation path.

[0052] In an embodiment of the present invention, the deviation of each intersection in the intersection sequence and the total deviation of the intersection sequence are first calculated according to the electromagnetic wave propagation estimated intersection deviation evaluation formula. Optionally, the deviation of each intersection in the intersection sequence is calculated using the following electromagnetic wave propagation estimated intersection deviation evaluation formula:

[0053]

[0054] Among them, Δ i is the deviation between the estimated path of electromagnetic wave propagation and the intersection point of the i-th layer medium interface, C i 、C i-1 、C i+1 are the intersection points of the estimated electromagnetic wave propagation path and the interfaces of the i-th, i-1-th, and i+1-th layers in the intersection sequence, respectively. i , ε i+1 are the dielectric constants of the i-th layer and the i+1-th layer respectively, are the thickness of the i-th layer and the i+1-th layer respectively. i If it is greater than 0, it means that C i Should move toward the aperture point on the medium interface, if Δ i If it is less than 0, it means C i Should move toward the imaging unit on the medium interface. |Δ i The larger the | is, the greater the C i The greater the deviation is. Since C0 is the coordinate of the aperture point, C m is the imaging unit coordinate, so Δ0=Δ m = 0. Calculate the sum of the absolute values ​​of the deviations of the intersections in the intersection sequence to obtain the total deviation Δ of the intersection sequence:

[0055] Then, the estimation interval is set according to each intersection in the intersection sequence and the deviation of each intersection. The intersection with the largest absolute value of the deviation of the intersection sequence I(C) is recorded as the kth element. According to the designed local optimal solution method, the local optimal solution is performed on the intersection I(k). The estimation interval is set as (ρ l ,ρ r ),in,

[0056]

[0057]

[0058] Iteratively update the intersection point with the largest absolute value of deviation in the intersection point sequence as the midpoint of the estimated interval, update the intersection point sequence, and calculate the total deviation of the updated intersection point sequence; if the updated total deviation is greater than the total deviation before the update, modify the intersection point in the upper and lower limits of the estimated interval that is not equal to the largest absolute value of deviation and update it as the midpoint until the updated total deviation is less than the total deviation before the update. That is, perform iterative operation, and take the estimated interval (ρ l ,ρ r ) midpoint ρ mid is the intersection point between the electromagnetic wave propagation path and the interface of the k-th layer medium, and the midpoint ρ mid Substitute the updated intersection sequence I(C) into the intersection sequence I(C) and calculate the total deviation Δ′ of the updated intersection sequence I′(C). If the updated total deviation Δ′ is greater than or equal to the total deviation Δ before the update, it means that the current selected intersection point is not the local optimal solution. Then modify the estimation interval (ρ l ,ρ r ) The upper and lower limits are not equal to C k The endpoint is the midpoint ρ mid , and enter the iteration again; if the updated total deviation Δ′ is less than the total deviation Δ before the update, it means that the selected intersection point is the local optimal solution, then the intersection point with the largest absolute value of the deviation is updated to the midpoint ρ mid , completing the local optimal solution. Repeat the iteration and stop when the total deviation of the intersection sequence is less than the preset threshold. The local optimal solution at this time is the global optimal solution.

[0059] Step S133: calculating the round-trip propagation delay of the electromagnetic wave between the aperture point and each imaging unit according to the estimated electromagnetic wave propagation path and the electromagnetic wave velocity formula.

[0060] The formula for the electromagnetic wave velocity is:

[0061]

[0062] Where c is the propagation speed of electromagnetic waves in a vacuum, and ε is the relative dielectric constant of the medium.

[0063] The aperture point A is obtained according to the electromagnetic wave propagation estimated path and the electromagnetic wave velocity formula. j With each imaging unit C f The two-way propagation delay τ of the electromagnetic wave between jf .

[0064] Step S134: extracting the signal amplitude of the A-scan echo data corresponding to the aperture point according to the round-trip propagation delay of the electromagnetic wave, and accumulating it to the imaging value of the corresponding imaging unit.

[0065] Optionally, according to the electromagnetic wave round-trip propagation delay τjf The jth aperture point A j A-scan echo data j The signal amplitude in (t) Among them, round(·) means rounding to the nearest integer, Δt is the time window sampling interval, T is the time window length, W is the number of sampling points in the time dimension of the A-scan echo data; the signal amplitude The imaging value is accumulated to the imaging unit.

[0066] Step S135: adding the imaging units adjacent to the imaging unit that have not been iteratively calculated to the imaging unit queue, and repeatedly iteratively calculating the imaging value of each imaging unit until the imaging unit queue is empty.

[0067] The imaging units adjacent to the imaging unit in the up, down, left, right, front and back directions that have not been iterated are added to the imaging unit queue, and then the iteration is repeated in sequence according to the order of the imaging units in the imaging unit queue to calculate the imaging value of each imaging unit in the sub-imaging matrix corresponding to the thread until the imaging unit queue is empty. That is, the imaging value of each imaging unit in the sub-imaging matrix corresponding to the thread is completely calculated. In the embodiment of the present invention, a breadth-first search algorithm is used to search each imaging unit in the sub-imaging matrix, and the back-projection imaging calculation is performed in this order. Of course, in other embodiments of the present invention, a depth-first search algorithm can also be used to search each imaging unit in the sub-imaging matrix, and the back-projection imaging calculation is performed in sequence. At this point, the calculation of the imaging value of each imaging unit in the corresponding sub-imaging matrix from any aperture point is completed.

[0068] Step S136: traverse each aperture point to obtain the final imaging value of each imaging unit in the corresponding sub-imaging matrix.

[0069] Traverse each aperture point and complete the final imaging value of each imaging unit in the sub-imaging matrix corresponding to the thread.

[0070] Each thread concurrently calculates the final imaging value of each imaging unit in the corresponding sub-imaging matrix. The final imaging value of each imaging unit in each sub-imaging matrix is ​​obtained by calculation by each thread, and each sub-imaging matrix is ​​combined to obtain a three-dimensional imaging matrix.

[0071] Step S14: normalizing and imaging the three-dimensional imaging matrix to obtain a concurrent back-projection imaging result.

[0072] The three-dimensional imaging matrix obtained above is normalized and imaged to obtain a concurrent back-projection imaging result.

[0073] The following example illustrates the use of gprMax ground penetrating radar forward simulation software to obtain C-scan echo data. The data processing platform is equipped with an Intel i7-7700 central processing unit and an NVIDIA GeForce GTX 1060 graphics processor, with a main memory capacity of 16GB. Set the three-dimensional coordinate system Oxyz, and the plane where the aperture point is located is the Oxy plane. The survey line arrangement direction is parallel to the x-axis and has a length of 1 meter. The survey line direction is parallel to the y-axis and has a length of 1 meter. A total of 8 survey lines are set, and the interval between each survey line is 0.1 meters. The starting survey line is the straight line x = 0.15m. The starting point of each survey line is the intersection of the straight line y = 0.15m and the survey line, with a step of 0.1 meter. A total of 8 survey points are set. A total of 64 aperture points are set, and the aperture points are arranged as follows Figure 3 shown.

[0074] Forward 3D model Figure 4 As shown, the z-axis represents the depth dimension of the detection area, and the plane of the aperture point is z = 0 m. Three different layers of media are set within the detection area. The first layer has an upper surface at z = 0 m and a lower surface at z = 0.5 m, with a relative dielectric constant of 2.5. The second layer has an upper surface at z = 0.5 m and a lower surface at z = 0.8 m, with a relative dielectric constant of 4.5. The third layer has an upper surface at z = 0.8 m and a lower surface at z = 1.5 m, with a relative dielectric constant of 6. A metal spherical target with a radius of 0.02 m is placed at (0.4 m, 0.6 m, 1.2 m). A time window of 20 ns is set to acquire C-scan echo data and preprocess it.

[0075] Based on the detection area, the 3D imaging area was set to 1m×1m×1.5m, and the length, width, and height of the imaging unit were all set to 1cm, resulting in a 3D imaging matrix O(x,y,z) with dimensions of 100×100×150. Based on the distribution of the medium in the detection area, the 3D imaging matrix O(x,y,z) was divided into three sub-matrices: O1(x,y,z), O2(x,y,z), and O3(x,y,z). O1(x,y,z) corresponds to the first layer of medium, with dimensions of 100×100×50; O2(x,y,z) corresponds to the second layer of medium, with dimensions of 100×100×30; and O3(x,y,z) corresponds to the third layer of medium, with dimensions of 100×100×70.

[0076] The thread pool is created by the ThreadPoolExecutor open source tool under the concurrent open source framework in Python. The number of core threads and the maximum number of threads in the thread pool are both 3. The i-th thread t i Responsible for the i-th sub-imaging matrix O i Back projection imaging of (x, y, z). The following takes the third thread t3 as an example.

[0077] The third thread t in the thread pool i Responsible for back-projection imaging of O3(x,y,z). When performing back-projection imaging, first traverse each aperture point, and then traverse the jth aperture point A. j At (0.15m, 0.15m, 0m), create imaging unit queue Q j (C), and O i Any imaging unit C0 in (x, y, z) is placed in Q j In (C), the imaging unit stored in the imaging unit queue must not have performed imaging operations. Assume that the coordinates of C0 are (0.3m, 0.5m, 1.2m).

[0078] Then, from Q j (C) Take out the first element C0, determine the electromagnetic wave propagation profile based on the aperture point, the imaging unit and the projection point of the imaging unit to the plane where the aperture point is located, and create a two-dimensional rectangular coordinate system xOz with the aperture point as the origin, the straight line determined by the aperture point and the projection point of the imaging unit to the plane where the aperture point is located as the x-axis, and the direction perpendicular to the aperture point plane as the z-axis. Therefore, see Figure 5 In each figure, in this two-dimensional coordinate system, the coordinates of the aperture point are (0, 0) and the coordinates of the imaging unit are (0.38m, 1.2m).

[0079] The adjacent imaging units of the imaging unit are traversed, and the adjacent imaging units have not undergone backprojection imaging operations. Therefore, the initial value of the intersection sequence I(C) between the estimated electromagnetic wave propagation path and the interface of each layer of the medium is obtained by the passive initialization method: I(C) = [0, 0.16, 0.25, 0.38].

[0080] Based on the greedy approach, the intersection sequence I(C) is iteratively solved for local optimality to obtain a global optimal solution. At each iteration, the deviations of each element in I(C) and the total deviation Δ of I(C) are calculated using the electromagnetic wave propagation estimated intersection deviation evaluation formula and the electromagnetic wave propagation estimated path and medium interface intersection sequence total deviation evaluation formula. The electromagnetic wave propagation estimated intersection deviation evaluation formula and the intersection sequence total deviation evaluation formula are described above and will not be repeated here.

[0081] The local optimal solution is performed for the element with the largest absolute value of deviation in I(C). The maximum value of deviation in I(C) is the kth element. According to the relative dielectric constant ρ of the upper layer of the kth dielectric interface k , relative dielectric constant ε of the lower layer medium k+1 , the intersection point C of the estimated electromagnetic wave propagation path and the interface of the kth layer medium k and its deviation Δ k , set the estimation interval to (ρ l ,ρr ), where the upper limit of the estimation interval ρ r and the lower limit ρ l The calculation method of is described above.

[0082] Perform iterative operations, and in each iteration, take the midpoint of the estimated interval ρ mid is the intersection point to be selected between the electromagnetic wave propagation path and the interface of the k-th layer of medium, and the intersection sequence is calculated and the k-th element is updated to ρ mid The total deviation of the electromagnetic wave propagation path after Δ′ is calculated. If Δ′ is greater than or equal to Δ, it means that the current selected intersection point is not a local optimal solution. Then the upper and lower limits of the estimated interval are modified to be not equal to C. k The endpoint is ρ mid , and enter the iteration again; if Δ′ is less than Δ, it means that the selected intersection point is the local optimal solution, then modify C k The value of the current ρ mid The local optimal solution is completed.

[0083] Until the total deviation of I(C) is less than the threshold value 0.1, I(C) is the global optimal solution.

[0084] In this embodiment, a total of 6 iterations are performed, and the iterative solution process is shown in Table 1 below. During the iterative process, the electromagnetic wave propagation path is calculated as follows: Figure 5 shown.

[0085] Table 1 Iterative solution process

[0086] Number of iterations I(C) Deviation Total deviation Medium interface with maximum deviation value legend 1 [0,0.16,0.25,0.38] [0,0.35,0.30,0] 0.65 1 a 2 [0,0.18,0.25,0.38] [0,0.14,0.44,0] 0.58 2 b 3 [0,0.18,0.29,0.38] [0,0.28,0.25,0] 0.43 1 c 4 [0,0.21,0.29,0.38] [0,0.20,0.08,0] 0.28 1 d 5 [0,0.19,0.29,0.38] [0,0.11,0.11,0] 0.22 2 e 6 [0,0.19,0.28,0.38] [0,0.04,0.03,0] 0.07 1 f

[0087] Therefore, the intersection points of the estimated electromagnetic wave propagation path and the interfaces of each layer of media are [0, 0.19, 0.28, 0.38]. According to the electromagnetic wave velocity formula, the round-trip delay τ of electromagnetic wave propagation can be calculated. Based on this delay value, the signal amplitude of the aperture point A-scan echo data is added to the imaging value of the imaging unit.

[0088] If from Q j (C) Take out the first element of the queue as C1 (0.30m, 0.50m, 1.19m). Since the imaging unit C0 adjacent to C1 has completed the backward projection imaging operation, the aperture point A is initialized through the active initialization method. j The intersection sequence of the electromagnetic wave propagation estimation path between the imaging unit C1 and the medium interface is the aperture point A j The estimated electromagnetic wave propagation path between imaging unit C0 and the dielectric interface intersection sequence is calculated and globally optimized. Global optimization is achieved after only one iteration, with the estimated electromagnetic wave propagation path and the dielectric interface intersection sequence being [0, 0.17, 0.28, 0.38].

[0089] The final three-dimensional back projection image is as follows Figure 6 As shown, from Figure 6 It can be seen that the coordinates of the imaging focus point are (0.4m, 0.6m, 0.98m), the deviation distance from the preset target is 0.02 meters, and the total imaging calculation time is about 100 seconds, which is about 250 times faster than the traditional method.

[0090] The embodiment of the present invention uses thread pool technology to create multiple threads for concurrent backprojection imaging calculations based on concurrent programming and spatial locality, thereby achieving concurrent backprojection imaging. Based on spatial locality, the intersection of the electromagnetic wave propagation path and the medium interface is initialized. The intersection of the electromagnetic wave propagation path and each layer of interface is traversed to calculate the deviation of the intersection of the electromagnetic wave propagation path and the medium interface. The intersection of the medium interface with the maximum deviation is iteratively calculated using a bisection method until the maximum deviation is less than the imaging unit size, thereby obtaining a sequence of intersections between the electromagnetic wave propagation path and each layer of the medium interface. Therefore, compared with the prior art, the embodiment of the present invention does not need to calculate the electromagnetic wave propagation path under multi-layer media conditions through an equation solving method with high computational complexity. Instead, the electromagnetic wave propagation path can be obtained through iterative calculation. The spatial locality initialization method greatly reduces the number of iterative calculations, and combined with concurrent programming, the calculation time of three-dimensional backprojection imaging is greatly shortened.

[0091] In summary, the concurrent back-projection imaging method for multi-layer media by a ground-penetrating radar of an embodiment of the present invention obtains a three-dimensional imaging matrix by setting an imaging area according to the detection area and discretizing it, and divides the three-dimensional imaging matrix into a plurality of sub-imaging matrices equal to the number of medium layers in the detection area according to the medium distribution in the detection area, and one sub-imaging matrix corresponds to one layer of medium; a thread pool is created according to the number of medium layers in the detection area, and the number of core threads and the maximum number of threads of the thread pool are equal to the number of medium layers in the detection area; each thread in the thread pool concurrently performs back-projection imaging operations on the corresponding sub-imaging matrix according to the echo data of the ground-penetrating radar to obtain the imaging value of each imaging unit in each sub-imaging matrix; the three-dimensional imaging matrix is ​​normalized and imaged to obtain a concurrent back-projection imaging result, which can reduce the computational complexity of the electromagnetic wave propagation path in the multi-layer medium background in the traditional back-projection imaging algorithm, and realizes concurrent back-projection imaging in three-dimensional space through the thread pool technology, greatly shortening the calculation time of three-dimensional imaging.

[0092] The foregoing description is of specific embodiments of the present invention. In some cases, the actions or steps described in the embodiments of the present invention may be performed in an order different from that shown in the embodiments and still achieve the desired results. In addition, the processes depicted in the accompanying drawings do not necessarily require the specific order or sequential order shown to achieve the desired results. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0093] The embodiment of the present invention also provides a multi-layer medium ground penetrating radar concurrent back projection imaging device, such as Figure 7 As shown, the concurrent back-projection imaging device for ground penetrating radar of multi-layer media includes: an imaging matrix division unit, a thread pool creation unit, a back-projection imaging unit and a normalization unit.

[0094] An imaging matrix division unit is configured to set an imaging region according to the detection region and discretize the region to obtain a three-dimensional imaging matrix. Based on the distribution of media in the detection region, the three-dimensional imaging matrix is ​​divided into a plurality of sub-imaging matrices equal to the number of media layers in the detection region, with one sub-imaging matrix corresponding to one layer of media.

[0095] A thread pool creation unit, configured to create a thread pool according to the number of media layers in the detection area, wherein the number of core threads and the maximum number of threads in the thread pool are equal to the number of media layers in the detection area;

[0096] A backprojection imaging unit, configured to concurrently perform a backprojection imaging operation on a corresponding sub-imaging matrix according to the echo data of the ground penetrating radar through each thread in the thread pool, to obtain an imaging value of each imaging unit in each sub-imaging matrix;

[0097] The normalization unit is used to normalize and image the three-dimensional imaging matrix to obtain a concurrent back-projection imaging result.

[0098] For the convenience of description, the above device is described as being divided into various units according to their functions. Of course, when implementing the embodiments of the present invention, the functions of each unit can be implemented in the same or multiple software and / or hardware.

[0099] The apparatus of the above embodiment is applied to the corresponding method of the above embodiment and has the beneficial effects of the corresponding method embodiment, which will not be described in detail here.

[0100] Based on the same inventive concept, an embodiment of the present invention further provides an electronic device, which includes a memory, a processor, and a computer program stored in the memory and runnable on the processor, wherein when the processor executes the program, the method described in any one of the above embodiments is implemented.

[0101] An embodiment of the present invention provides a non-volatile computer storage medium, wherein the computer storage medium stores at least one executable instruction, and the computer executable instruction can execute the method described in any one of the above embodiments.

[0102] Figure 8805. A more specific schematic diagram of the hardware structure of an electronic device provided in this embodiment is shown. The device may include: a processor 801, a memory 802, an input / output interface 803, a communication interface 804, and a bus 805. The processor 801, the memory 802, the input / output interface 803, and the communication interface 804 are communicatively connected to each other within the device via the bus 805.

[0103] The processor 801 can be implemented as a general-purpose CPU (Central Processing Unit), a microprocessor, an application-specific integrated circuit (ASIC), or one or more integrated circuits, and is used to execute relevant programs to implement the technical solutions provided by the method embodiments of the present invention.

[0104] The memory 802 can be implemented in the form of ROM (Read Only Memory), RAM (Random Access Memory), static storage devices, dynamic storage devices, etc. The memory 802 can store an operating system and other application programs. When the technical solutions provided by the method embodiments of the present invention are implemented through software or firmware, the relevant program codes are stored in the memory 802 and called and executed by the processor 801.

[0105] The input / output interface 803 is used to connect to input / output modules to implement information input and output. The input / output modules can be configured as components in the device (not shown in the figure) or can be externally connected to the device to provide corresponding functions. Input devices may include a keyboard, mouse, touch screen, microphone, various sensors, etc., and output devices may include a display, speaker, vibrator, indicator light, etc.

[0106] The communication interface 804 is used to connect to a communication module (not shown) to enable communication between the device and other devices. The communication module can communicate via a wired method (such as USB, network cable, etc.) or a wireless method (such as mobile network, WIFI, Bluetooth, etc.).

[0107] The bus 805 comprises a pathway for transmitting information between the various components of the device (eg, the processor 801 , the memory 802 , the input / output interface 803 , and the communication interface 804 ).

[0108] It should be noted that although the above device only shows the processor 801, memory 802, input / output interface 803, communication interface 804, and bus 805, in a specific implementation, the device may also include other components necessary for normal operation. In addition, those skilled in the art will understand that the above device may only include the components necessary to implement the embodiments of the present invention, and does not necessarily include all the components shown in the figure.

[0109] Those skilled in the art should understand that the discussion of any of the above embodiments is merely illustrative and is not intended to imply that the scope of the present disclosure (including the claims) is limited to these examples. Within the scope of the present disclosure, the technical features in the above embodiments or different embodiments may be combined, the steps may be implemented in any order, and there are many other variations of the different aspects of the embodiments of the present invention as described above, which are not provided in detail for the sake of simplicity.

[0110] The embodiments of the present invention are intended to cover all such substitutions, modifications, and variations that fall within the broad scope of the appended claims. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the embodiments of the present invention should be included in the scope of protection of this disclosure.

Claims

1. A method for concurrent backprojection imaging of multi-layered ground penetrating radar, characterized in that: The method for concurrent backprojection imaging of a multi-layered medium by ground penetrating radar includes: An imaging region is set according to the detection region and discretized to obtain a three-dimensional imaging matrix. According to the distribution of the medium in the detection region, the three-dimensional imaging matrix is ​​divided into a plurality of sub-imaging matrices equal to the number of medium layers in the detection region, with one sub-imaging matrix corresponding to one layer of medium; Create a thread pool according to the number of media layers in the detection area, wherein the number of core threads and the maximum number of threads in the thread pool are equal to the number of media layers in the detection area; Each thread in the thread pool concurrently performs a backprojection imaging operation on the corresponding sub-imaging matrix according to the echo data of the ground penetrating radar, so as to obtain an imaging value of each imaging unit in each sub-imaging matrix; The three-dimensional imaging matrix is ​​normalized and imaged to obtain a concurrent back-projection imaging result.

2. The method for concurrent back-projection ground penetrating radar imaging of a multi-layer medium according to claim 1, wherein: The method includes: for any thread, performing back-projection imaging operations on the corresponding sub-imaging matrices concurrently according to the echo data of the ground penetrating radar by each thread in the thread pool to obtain the imaging value of each imaging unit in each sub-imaging matrix. For any aperture point, create an imaging unit queue, and put any imaging unit in the sub-imaging matrix as an initial imaging unit into the imaging unit queue; Taking out imaging units in the imaging unit queue in turn, and calculating an estimated electromagnetic wave propagation path between the aperture point and the imaging unit; Calculating the electromagnetic wave round-trip propagation delay between the aperture point and each imaging unit according to the electromagnetic wave propagation estimated path and the electromagnetic wave velocity formula; Extracting the signal amplitude of the A-scan echo data corresponding to the aperture point according to the round-trip propagation delay of the electromagnetic wave, and accumulating it to the imaging value of the corresponding imaging unit; Adding the imaging units adjacent to the imaging unit that have not been iteratively calculated to the imaging unit queue, and repeatedly iteratively calculating the imaging value of each imaging unit until the imaging unit queue is empty; Traverse each aperture point to obtain the final imaging value of each imaging unit in the corresponding sub-imaging matrix.

3. The method for concurrent back-projection ground-penetrating radar imaging of a multi-layer medium according to claim 2, wherein: Determining the estimated electromagnetic wave propagation path between the aperture point and the imaging unit includes: Setting an electromagnetic wave propagation profile according to the aperture point and the imaging unit and establishing a two-dimensional rectangular coordinate system; The initial values ​​of the intersection sequence of the electromagnetic wave propagation estimation path and each layer of medium interface are obtained by initialization method; An iterative optimization solution is performed according to a greedy idea to obtain an estimated electromagnetic wave propagation path between the aperture point and each imaging unit in the corresponding sub-imaging matrix.

4. The method for concurrent back-projection ground-penetrating radar imaging of a multi-layer medium according to claim 3, wherein: The initialization method for obtaining the initial values ​​of the sequence of intersection points between the estimated electromagnetic wave propagation path and the interfaces of each layer of media includes: Traversing adjacent imaging units in directions above, below, left, right, front, and back of the imaging unit; If any of the adjacent imaging units has completed the solution of the electromagnetic wave propagation estimation path with the aperture point, then the calculation result of any of the adjacent imaging units that has been solved is taken as the initial value of the intersection sequence of the electromagnetic wave propagation estimation path and the interface between each layer of the medium; Otherwise, the intersection points of the estimated electromagnetic wave propagation path and each medium interface are initialized to the intersection points of the line from the aperture point to the imaging unit and each medium interface, to obtain the initial value of the intersection point sequence.

5. The method for concurrent back-projection ground-penetrating radar imaging of a multi-layer medium according to claim 3, wherein: The iterative optimization solution based on the greedy idea to obtain the electromagnetic wave propagation estimation path between the aperture point and each imaging unit in the corresponding sub-imaging matrix includes: Calculate the deviation of each intersection in the intersection sequence and the total deviation of the intersection sequence according to the intersection deviation evaluation formula estimated by electromagnetic wave propagation; Setting an estimation interval according to each intersection point in the intersection point sequence and the deviation of each intersection point; Iteratively updating the intersection point with the largest absolute value of deviation in the intersection point sequence as the midpoint of the estimation interval, updating the intersection point sequence, and calculating the total deviation of the updated intersection point sequence; If the updated total deviation is greater than the total deviation before the update, the intersection point between the upper and lower limits of the estimated interval that is not equal to the maximum absolute value of the deviation is updated to the midpoint until the updated total deviation is less than the total deviation before the update; The intersection point with the largest absolute value of the deviation is updated as the midpoint, and the iteration is repeated. The iteration is stopped when the total deviation of the intersection point sequence is less than a preset threshold.

6. The method for concurrent back-projection ground-penetrating radar imaging of a multi-layer medium according to claim 5, wherein: The calculation of the deviation of each intersection in the intersection sequence and the total deviation of the intersection sequence according to the intersection deviation evaluation formula estimated by electromagnetic wave propagation includes: The deviation of each intersection point in the intersection point sequence is calculated using the following electromagnetic wave propagation estimation intersection deviation evaluation formula: Among them, Δ i is the deviation between the estimated path of electromagnetic wave propagation and the intersection point of the i-th layer medium interface, C i 、C i-1 、C i+1 are the intersection points of the estimated electromagnetic wave propagation path and the interfaces of the i-th, i-1-th, and i+1-th layers in the intersection sequence, respectively. i , ε i+1 are the dielectric constants of the i-th layer and the i+1-th layer respectively, are the thicknesses of the i-th and i+1-th layers of dielectric respectively; The sum of the absolute values ​​of the deviations of the intersection points in the intersection point sequence is calculated to obtain the total deviation of the intersection point sequence.

7. The method for concurrent back-projection ground-penetrating radar imaging of a multi-layer medium according to claim 2, wherein: The extracting the signal amplitude of the A-scan echo data corresponding to the aperture point according to the round-trip propagation delay of the electromagnetic wave and accumulating the signal amplitude to the imaging value of the corresponding imaging unit includes: According to the electromagnetic wave round-trip propagation delay τ jf The jth aperture point A j A-scan echo data j The signal amplitude in (t) Among them, round(·) means rounding to the nearest integer, Δt is the time window sampling interval, T is the time window length, W is the number of sampling points in the time dimension of the A-scan echo data; The signal amplitude The imaging value is accumulated to the imaging unit.

8. A multi-layered ground penetrating radar concurrent back-projection imaging device, characterized in that: The multi-layer medium ground penetrating radar concurrent back projection imaging device comprises: An imaging matrix division unit is configured to set an imaging region according to the detection region and discretize the region to obtain a three-dimensional imaging matrix. Based on the distribution of media in the detection region, the three-dimensional imaging matrix is ​​divided into a plurality of sub-imaging matrices equal to the number of media layers in the detection region, with one sub-imaging matrix corresponding to one layer of media. A thread pool creation unit, configured to create a thread pool according to the number of media layers in the detection area, wherein the number of core threads and the maximum number of threads in the thread pool are equal to the number of media layers in the detection area; A backprojection imaging unit, configured to concurrently perform a backprojection imaging operation on a corresponding sub-imaging matrix according to the echo data of the ground penetrating radar through each thread in the thread pool, to obtain an imaging value of each imaging unit in each sub-imaging matrix; The normalization unit is used to normalize and image the three-dimensional imaging matrix to obtain a concurrent back-projection imaging result.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the computer program, the method for concurrent backprojection ground penetrating radar imaging of a multi-layer medium according to any one of claims 1 to 7 is implemented.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method for concurrent backprojection ground penetrating radar imaging of a multi-layer medium according to any one of claims 1 to 7 is implemented.

Images