CN109283523B - Geological radar B-scan data processing method - Google Patents

Geological radar B-scan data processing method Download PDF

Info

Publication number
CN109283523B
CN109283523B CN201810867342.1A CN201810867342A CN109283523B CN 109283523 B CN109283523 B CN 109283523B CN 201810867342 A CN201810867342 A CN 201810867342A CN 109283523 B CN109283523 B CN 109283523B
Authority
CN
China
Prior art keywords
scan
data
frequency
section
imaging
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810867342.1A
Other languages
Chinese (zh)
Other versions
CN109283523A (en
Inventor
周艳辉
师振盛
王晓凯
陈文超
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201810867342.1A priority Critical patent/CN109283523B/en
Publication of CN109283523A publication Critical patent/CN109283523A/en
Application granted granted Critical
Publication of CN109283523B publication Critical patent/CN109283523B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/885—Radar or analogous systems specially adapted for specific applications for ground probing
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

The invention discloses a geological radar B-scan data processing method, which utilizes a signal processing method to process B-scan data so as to separate and obtain a data profile for detecting tunnel lining abnormity and railway roadbed diseases. The method comprises the steps of firstly improving the time resolution of signals by adopting a continuous spectrum whitening spectrum expansion technology, then carrying out offset imaging on high-resolution data by adopting efficient phase shift offset, improving the transverse resolution of a section, focusing reinforcing steel bars inside a lining or a roadbed on the imaged section into a point-shaped structure, and finally separating different morphological components such as a reinforcing steel bar structure and a reflecting interface by adopting a sparse separation method based on morphological component analysis based on the difference of layer characteristics such as a reinforcing steel bar point target and the reflecting interface (a concrete and clay layer or surrounding rock interface). The method is beneficial to improving the accuracy of tunnel lining and railway roadbed detection, and provides a powerful basis for tunnel and railway maintenance.

Description

Geological radar B-scan data processing method
Technical Field
The invention belongs to the field of signal processing of geological radar signals, and particularly relates to a geological radar B-scan data processing method.
Background
Tunnel construction is an important problem in road and railway construction, wherein tunnel lining construction is an important link in tunnel construction at present. Later-stage tunnel lining detection is an important judgment on construction quality. On the other hand, the railway subgrade diseases induce subgrade instability, which is an important factor influencing the safe operation of railways, and the development of an effective railway subgrade detection method has important significance. A geological Radar (GPR) disposed on the surface of a detection target emits a broadband short-pulse high-frequency electromagnetic wave through a transmitting antenna, and a reflected signal generated when the electromagnetic wave passes through a dielectric body or an interface having an electrical difference in the propagation process of a lining, a roadbed, and a deeper medium is received by a receiving antenna disposed on the surface layer. The geological radar has the unique advantage of ultra-shallow exploration, is a rapid high-resolution nondestructive testing method, and has wide application in tunnel and roadbed detection.
At present, the signal processing method based on geological radar B-scan data mainly comprises
Prior art 1: mean value subtraction method
(1) Selecting a plurality of continuous channels on the B-scan section, accumulating the values of the channels corresponding to the time points and dividing the values by the number of the channels to obtain an average channel;
(2) subtracting the average track in the step (1) from each track on the B-scan to obtain a final data processing result;
disadvantages of prior art 1:
because actual geological radar data are complex, the number of tracks for calculating an average track is difficult to select, and when the selection is not correct, a useful signal is easy to be damaged, so that the subsequent lining abnormity detection is influenced; the method assumes that the clutter and the interface reflection signal are in a strict horizontal structure on the B-scan section, and the actual received data generally does not meet the assumption, so the processing effect of the method is limited;
prior art 2: filtering method
On a B-scan section, the main components of clutter and interface reflected signals are generally represented as horizontal banded structures, the main components of clutter are generally at the top of the section, the steel bar reflected signals are in a parabolic state, and a time-space domain filter or a frequency-wave number domain filter can be designed by utilizing the characteristic difference of the clutter and the interface reflected signals in a time-space domain or different distribution areas in a frequency-wave number domain, so that the two components are separated;
disadvantages of prior art 2:
these filters usually set filter parameters based on numerical simulation, and there is a certain deviation between the actual B-scan profile of the geological radar and the simulated profile, because the set filter parameters are used for the actual B-scan profile processing, which causes a problem that the filter parameters are difficult to set properly, thereby affecting the subsequent separation effect and detection effect.
Disclosure of Invention
Based on this, the present disclosure discloses a geological radar B-scan data processing method, which is characterized by comprising the following steps:
s100, processing each channel of data (A-scan) on the B-scan section by adopting a spectral whitening frequency broadening method to obtain a B-scan section with high time resolution;
s200, performing two-dimensional offset imaging on the high-time-resolution B-scan data based on a phase shift offset imaging method to obtain a B-scan imaging section;
s300, carrying out iterative sparse separation on the B-scan imaging section obtained in the step S200 based on a morphological component analysis method;
s400, analyzing the point-like characteristics and the layered structure of the results obtained through iterative sparse separation.
The invention has the following beneficial effects:
1) the spectral whitening spectrum expanding method can effectively improve the time resolution of the B-scan section and enhance the intensity of deep weak reflection;
2) the method can be used for efficiently imaging the B-scan section by adopting phase shift offset, so that on one hand, hyperbolic type scattering signals generated by the steel bars are focused to the underground position of the steel bars, the steel bars are imaged to be in a point-shaped characteristic, and reflected waves generated by the interface of a concrete layer and a clay layer or surrounding rock are still in a layered structure after being imaged, so that the two characteristics are effectively distinguished, which shows that the phase shift offset can improve the transverse resolution of the B-scan section and is beneficial to transverse fine detection of a target;
3) the separation method disclosed by the invention utilizes the difference of sparse representation characteristics of different structures to carry out sparse separation on different morphological components, specifically adopts two-dimensional discrete wavelet transform without downsampling to carry out sparse representation on a point-shaped structure, and carries out sparse representation on a layered structure by two-dimensional curvelet transform. This achieves consistency of the transformation basis with the analyzed features, and thus the point-like structures and the layered structures can be effectively separated.
Drawings
FIG. 1 is a geological radar B-scan data processing flow in one embodiment of the present disclosure;
FIG. 2 is a flow chart of a Fourier domain B-scan data spectrum expansion based on a spectral whitening method in an embodiment of the present disclosure;
FIG. 3(A) is a sparse separation process of morphological component analysis in one embodiment of the present disclosure;
FIG. 3(B) is an atomic illustration of a downsamplionless wavelet in one embodiment of the present disclosure;
FIG. 3(C) is a curvelet atomic schematic in one embodiment of the disclosure;
FIG. 4(A) is a schematic representation of a geometric model according to an embodiment of the present disclosure;
FIG. 4(B) is a B-scan profile generated using a finite difference method in one embodiment of the present disclosure;
FIG. 4(C) is a B-scan cross section after deflection according to one embodiment of the present disclosure;
FIG. 4(D) is a diagram of a separated dot-like structure component in one embodiment of the present disclosure;
FIG. 4(E) is a separate lamination component in one embodiment of the present disclosure;
FIG. 5(A) is a measured B-scan profile according to an embodiment of the present disclosure;
FIG. 5(B) is a high resolution B-scan profile after spectral whitening frequency broadening processing in one embodiment of the present disclosure;
FIG. 5(C) is a B-scan cross section after deflection according to one embodiment of the present disclosure;
FIG. 5(D) is a diagram of a separated dot-like structure component in one embodiment of the present disclosure;
FIG. 5(E) is a separated layered structure component in one embodiment of the present disclosure;
FIG. 6(A) is a measured B-scan profile according to an embodiment of the present disclosure;
FIG. 6(B) is a diagram of a separated dot-like structural component in one embodiment of the present disclosure;
fig. 6(C) is a separate lamination component in one embodiment of the present disclosure.
Detailed Description
The present invention will be described in further detail with reference to the accompanying fig. 1 to 6(C) and the embodiments.
In one embodiment, the present disclosure discloses a method for geo-radar B-scan data processing, the method comprising the steps of:
s100, processing each channel of data (A-scan) on the B-scan section by adopting a spectral whitening frequency broadening method to obtain a B-scan section with high time resolution;
s200, performing two-dimensional offset imaging on the high-time-resolution B-scan data based on a phase shift offset imaging method to obtain a B-scan imaging section;
s300, carrying out iterative sparse separation on the B-scan imaging section obtained in the step S200 based on a morphological component analysis method;
s400, analyzing the two components subjected to iterative sparse separation to obtain the dot characteristics and the layered structure.
Further, if the time resolution of a certain piece of data of the B-scan profile is high, step S100 is omitted.
In order to solve the problems in the prior art, the embodiment provides a B-scan profile signal processing method for a geological radar to serve for later lining or roadbed detection according to target signal characteristics required to be detected. Firstly, a signal frequency spectrum expanding (frequency broadening) technology based on spectral whitening is adopted to improve the resolution of a two-dimensional geological radar B-scan data signal in the time direction, and if the time resolution of the received two-dimensional geological radar data is higher, namely the bandwidth of a certain data is more than 5 octaves, the step can be omitted. And finally, respectively obtaining targets such as a reinforcing mesh and a reflection interface by adopting a sparse separation method based on morphological component analysis based on the difference of layer characteristics such as a reinforcing steel bar isopoint target and a reflection interface (a concrete and clay layer or surrounding rock interface).
In an embodiment, the frequency spreading method in step S100 specifically includes the following steps:
s101, calculating one-dimensional Fourier transform of each data on a B-scan section to obtain an amplitude spectrum and a phase spectrum of each data;
s102, setting a plurality of band-pass filters in the effective bandwidth range of each channel of data, and dividing an amplitude spectrum into a plurality of sub-bands;
s103, calculating a data amplitude spectrum weighting coefficient in each sub-band;
s104, multiplying the amplitude spectrum in each sub-band by a corresponding weighting coefficient, keeping the phase spectrum unchanged, and obtaining a frequency spectrum after frequency extension;
and S105, calculating one-dimensional Fourier inverse transformation of each data on the B-scan section based on the frequency extended frequency spectrum to obtain each data on the frequency extended B-scan section.
In this embodiment, the number of bandpass filters in step S102 is 5-8.
In one embodiment, the phase shift offset imaging method in step S200 includes the steps of:
s201, performing two-dimensional Fourier transform on the B-scan data after the high-time resolution processing, wherein the corresponding relation is as follows:
Figure BDA0001750045610000061
wherein f (x, z, t) & gtnon phosphorz＝0Representing the B-scan data received at the surface layer of the detection target, x representing the spatial coordinates, z representing the depth, t representing the time,
Figure BDA0001750045610000062
is an imaginary unit; f (k)x，z，ω)|z＝0Is a two-dimensional Fourier transform of the B-scan signal, kxRepresents a spatial frequency (spatial wave number), ω is an angular frequency;
s202, applying a phase shift operator to a wave field F (k) with the depth zxZ, ω) to obtain a wavefield at a depth of z + Δ z;
and S203, obtaining a phase shift offset imaging result of the B-scan based on two-dimensional inverse Fourier transform and making the time t equal to 0.
In one embodiment, the iterative sparse separation in step S300 is determined by the following relation:
Figure BDA0001750045610000063
wherein f is a B-scan section after offset imaging, f1Representing structural components of the bars and other points therein, f2Representing the horizontal or near horizontal lamellar component thereof,
Figure BDA0001750045610000071
is a two-dimensional non-down-sampling discrete wavelet transform base,
Figure BDA0001750045610000072
is f1The L1 norm of the wavelet transform coefficients,
Figure BDA0001750045610000073
is a two-dimensional curvelet transform basis,
Figure BDA0001750045610000074
is f2The L1 norm of the curvelet transform coefficients,
Figure BDA0001750045610000075
and lambda is a regularization parameter which is an error term corresponding to the error and satisfies lambda > 0.
In one embodiment, the sparsely separated result in step S400 includes two components: one is a profile that includes rebar or other point-like features; the other is a flat layer structure, wherein the shallowest layer is a layer obtained by imaging direct waves from a transmitting antenna to a receiving antenna and reflected waves of the transmitting antenna on the surface layer of the detection target, and the deep layer is an interface between a lining and surrounding rocks or clay and other medium layers.
In one embodiment, the step S103 specifically includes: and calculating the amplitude root mean square value of each sub-band, and dividing a fixed constant by the root mean square value to obtain a weighting coefficient of the amplitude spectrum in the corresponding sub-band.
In this embodiment, the fixed constant is an overall scale factor, and may be selected from 2000 and 3000, without any fixed requirement.
In one embodiment, the calculation formula in step S202 is:
Figure BDA0001750045610000076
in the formula (I), the compound is shown in the specification,
Figure BDA0001750045610000077
for the phase shift operator, Δ z is the continuation step length, and θ is half of the propagation speed of the electromagnetic wave in the tunnel lining.
In one embodiment, the imaging result calculation formula of B-scan in step S202 is as follows
Figure BDA0001750045610000078
Then, according to the conversion relation between the depth z and upsilon, the imaging section is divided
Figure BDA0001750045610000081
Converted into time-space profiles
Figure BDA0001750045610000082
In one embodiment, the invention provides a geological radar data signal processing method for tunnel lining abnormity and roadbed disease detection, which comprises the steps of high-resolution frequency broadening processing, offset imaging and sparse separation based on morphological component analysis on geological radar data, so as to obtain a point structure component and a layer structure component for abnormity detection.
As shown in fig. 1, the geological radar data processing method provided by the present invention is specifically implemented by the following steps:
preparing a geological radar B-scan section, and processing the section as follows:
step 1: for each track on the B-scan profile, i.e. each A-scan, f0(t) performing Fourier domain frequency spectrum broadening processing to improve the resolution ratio;
referring to fig. 2, the fourier domain spectrum expansion process includes the steps of:
(1) for A-scan f0(t) performing a one-dimensional Fourier transform, the formula being as follows:
Figure BDA0001750045610000083
wherein t represents time, i is an imaginary unit, F0And (omega) is Fourier transform of the A-scan signal, and omega is angular frequency. Because the acquired data is discrete data, the implementation of the one-dimensional Fourier transform can be realized by adopting a fast algorithm of the discrete Fourier transform so as to improve the calculation efficiency.
(2) Setting K narrow band-pass filters H in effective frequency band range of signalkEach band pass filter having a frequency range of (ω)k-1，ωk) K is 1, 2, L, K, and the signal is frequency filtered to obtain K sub-band spectrums Hk{F0(ω)}，k＝1，2，L，K；
(3) Setting a constant factor S at each band pass filter HkFrequency range (ω) of (K ═ 0, 1, L, K)k-1，ωk) First, the average energy density P of the signal spectrum is calculatedkIs composed of
Figure BDA0001750045610000091
Wherein, | Hk{F0(ω) } | represents Hk{F0(ω) } modulo value;
then, a weighting factor W of the spectrum in the frequency range is calculatedkIs composed of
Figure BDA0001750045610000092
(4) For each band pass filter HkThe signal spectrum in the (K-0, 1, L, K) frequency range is multiplied by a corresponding weighting factor WkAnd obtaining an expanded frequency spectrum, wherein the formula is as follows:
Figure BDA0001750045610000093
(5) merging the frequency spectrum after frequency extension into
Figure BDA0001750045610000094
Performing one-dimensional Fourier inverse transformation to obtain a frequency extended signal, wherein the formula is as follows:
Figure BDA0001750045610000095
step 2: performing two-dimensional offset imaging on the B-scan data subjected to the high-resolution frequency broadening processing based on a phase shift offset imaging method;
first, pre-counting B-scan signal f (x, z, t)z＝0As an input signal, where z denotes depth, where z ═ 0 denotes B-scan data received at the surface layer, and x denotes a spatial coordinate, the phase shift processing mainly includes the steps of:
(1) to f (x, z, t) & gtnon-z＝0Performing two-dimensional Fourier transform, wherein the corresponding relation is as follows:
Figure BDA0001750045610000096
in the formula, F (k)x，z，ω)|z＝0Is a two-dimensional Fourier transform of the B-scan signal, kxRepresenting spatial frequencies (spatial wavenumbers). Because the acquired data is discrete data, the implementation of the two-dimensional Fourier transform is realized by adopting a fast algorithm of the discrete Fourier transform, and the data needs to be zero-filled in the space and time directions under the normal condition, so that the number N of space sampling pointsxNumber of sampling points in time NtAre all powers of 2, thereby improving computational efficiency. dt and dx are the sampling intervals in the temporal and spatial directions, corresponding to kxAnd ω are discrete sequence values with a sampling interval dkx＝1/(dxgNx)，dω＝2π/(dtgNt)，kxHas a value of 0, dkx，2dkx，L，(Nx-1)dkxThe value of omega is 0, d omega, 2d omega, L, (N)t-1)dω。
(2) Applying a phase-shift operator to the wavefield F (k) at depth zxK Δ z, ω), the wavefield at depth (k +1) Δ z can be extended to:
Figure BDA0001750045610000101
in the formula, Δ z is an extension step, and usually takes a value of 0.1 cm or 0.2 cm, etc., and the depth is discretized into z of 0, Δ z, 2 Δ z, L, K Δ z, L, (K-1) Δ z, and the number of extensions is K, and the value thereof depends on the detection depth. Upsilon is half of the propagation speed of the electromagnetic wave in the tunnel lining, is related to the relative dielectric constant of the reinforced concrete, and is set as
Figure BDA0001750045610000102
Wherein c is0299792458 m/s is the propagation velocity of light in vacuum, εrIs the relative dielectric constant of the medium. EpsilonrUsually taken as the existing empirical parameters, but subject to the constraints of specific construction conditions and the variation of concrete medium, in specific calculations, the epsilon is usually required to be finely adjustedrThe criterion is to make the reflected hyperbolic signal from the steel bars in the B-scan focus the most intense energy.
Figure BDA0001750045610000103
For the phase shift operator, in the specific calculation, when
Figure BDA0001750045610000104
When it is used, order
Figure BDA0001750045610000105
(C) Using an inverse two-dimensional fourier transform and letting time t equal to 0, we can reconstruct the B-scan signal as:
Figure BDA0001750045610000111
spatial frequency k in formulaxThe inverse Fourier transform is realized by a fast Fourier transform algorithm, and the superposition operation is directly carried out in the frequency omega direction. This is done for the wavefield at each depthAnd (5) operating.
And step 3: performing iterative sparse separation on the B-scan based on morphological component analysis, wherein the iterative sparse separation based on the morphological component analysis is determined by the following relational expression:
Figure BDA0001750045610000112
wherein f and
Figure BDA0001750045610000113
known as f1，f2Is the component to be solved. f is the B-scan section after offset, f1For the component 1 to be separated, the component of the steel bar structure at equal points is shown, f2The component 2 represents the horizontal layer component, and these are vectorized representations.
Figure BDA0001750045610000114
The wavelet function is a two-dimensional non-downsampling discrete wavelet transform base, is selected as a symmlet wavelet of 6 orders, can also be selected as a symmlet wavelet of other orders or a Daubechies wavelet and the like,
Figure BDA0001750045610000115
the number of the scales is 8 for the two-dimensional curvelet transform basis, and a basic atom schematic and a two-dimensional curvelet atom schematic of the two-dimensional non-downsampling discrete wavelet transform are respectively given in fig. 3(B) and fig. 3 (C).
Figure BDA0001750045610000116
Is f1The L1 norm of the wavelet transform coefficient vector,
Figure BDA0001750045610000117
is f2The L1 norm of the coefficient vector of the curvelet transform, defined as the sum of the absolute values of all the elements of the coefficient vector,
Figure BDA0001750045610000118
is defined as f-f1-f2The sum of the squares of all elements. And lambda is a regularization parameter and satisfies lambda > 0. As shown in FIG. 3(A), f1Is defined as f, f2Initialized to 0. In each iteration, the following steps are performed:
(1) fixed f2Solving an inverse problem defined by the following relation by adopting an Iterative soft-threshold (Iterative soft-threshold) method:
Figure BDA0001750045610000119
to obtain f1；
(2) Fixed f1Solving an inverse problem defined by the following relation by adopting an Iterative soft-threshold (Iterative soft-threshold) method:
Figure BDA0001750045610000121
to obtain f2。
As shown in fig. 3(a), the iteration termination condition is defined by the maximum iteration number, and is selected to be 300 times, so that the requirement of the calculation accuracy can be met, and the calculation efficiency can be remarkably improved.
And 4, step 4: analyzing the two separated components;
two components f separated in step 31And f2Respectively representing a point-shaped structure component and a layered structure component in the B-scan section after the offset, and according to the structural arrangement of the tunnel lining, the point-shaped component section mainly comprises point-shaped characteristics such as reinforcing steel bars, the layered structure component mainly comprises shallow layered clutter and mainly comprises a direct wave layer from a transmitting antenna to a receiving antenna and a reflected wave on the surface of the lining, and the deep layered structure is an interface of the lining and a backfill layer or a clay layer and surrounding rock.
Processing (1) simulated tunnel lining geological radar B-scan data and (2) actually measured B-scan data by using a data processing method provided by the disclosure; for the finite difference synthetic model data, because the medium is uniform, the attenuation of the reflected wave is reduced, and the time resolution is high, the high-resolution frequency extension processing of the step 1 is omitted, and only the processing of the step 2 to the step 4 is carried out. The measured B-scan data was analyzed for the treatments from step 1 to step 4.
Fig. 4(a) -4 (E) show model data and processing results thereof:
as shown in fig. 4(a), which is a geometric schematic of the model, the thickness of the concrete layer is 0.45 m, and a clay layer is arranged below the concrete. The distance between the steel bars and the surface layer is 0.05 meter, the radius of each steel bar is 0.01 meter, and the distance between the transverse centers of the adjacent steel bars is 0.1 meter. The signal of the emission source is a pulse source, the center frequency of the transmitting antenna is 1GHz, the transverse distance between the transmitting antenna and the receiving antenna is 0.01 m, and the transverse movement interval of the antenna is 0.01 m. Fig. 4(B) is a generated B-scan section with horizontal axis and vertical axis representing the receiving time, and the early received strong horizontal clutter and hyperbolic waveform of the reinforcement bar reflection can be seen. The hyperbolic waveforms reflected by the steel bars interfere with each other, so that the transverse resolution of the section is reduced, and the detection of the steel bars is not facilitated. The horizontal reflection layer energy of the interface of the concrete and the clay layer near 7ns is weaker; and part of the reflected signals are covered by the hyperbolic wave reflected by the steel bars, as shown by the black frame in the figure, so that the effective identification cannot be carried out. Fig. 4(C) is a B-scan section after phase shift imaging, with the hyperbolic waveform of the rebar reflection focused into an enhanced point-like structure, unlike the horizontal structure of clutter and deep reflections. Based on these two different structural components, the morphological component analysis in step 5 is used to obtain two components in fig. 4(D) and fig. 4(E) by sparse separation. Fig. 4(D) shows a dot structure component, and fig. 4(E) shows a horizontal layer structure component. It can be seen that the point-like structure of the reinforcing bars is mainly distributed in fig. 4(D), the horizontal clutter and the deep reflecting interface are mainly distributed in fig. 4(E), and the horizontal reflecting layer covered in fig. 4(B) is clearly shown in fig. 4 (E).
Fig. 5(a) -5 (E) show the data processing results of the measured tunnel lining geological radar:
fig. 5(a) is actually measured B-scan data for detecting lining anomaly, the main frequency of the transmitting antenna is 900MHz, the thickness of the tunnel segment is 35 cm, two layers of reinforcing steel bars are distributed therein, the radius of the main reinforcing steel bar is 0.01 m, and the radius of the reinforcing steel bar is 0.005 m. The shallow, intense energy laminar noise and the hyperbolic wave form of the bar reflections in the duct piece are more clearly seen in fig. 5 (a). Due to the dispersion effect of the actual medium, the interface reflection of the duct piece and the surrounding rock behind is weak on the graph 5(A), and the method is difficult to be used for effectively detecting the tunnel lining abnormity. Fig. 5(B) shows the result of the high-resolution frequency extension processing, and it can be seen that the time-direction resolution of the profile is effectively improved, and the energy of the deep reflected wave is effectively enhanced. Fig. 5(C) shows the offset imaging result, the hyperbolic wave of the steel bar reflection is effectively focused, and the deep reflection interface is more smoothly represented, so that the two represent different morphological structures. Fig. 5(D) and 5(E) show the results of separation using step 3. Fig. 5(D) shows a dot-shaped structural component, and the distribution of the upper and lower sets of reinforcing steel bars can be clearly seen. Fig. 5(E) is a horizontal layered structure component, and the reflection interface of the deep duct piece and the surrounding rock is clearly depicted.
Fig. 6(a) -6 (C) show the data processing results of the measured railway roadbed geological radar:
fig. 6(a) is actually measured B-scan data for railway roadbed detection, wherein the main frequency of the transmitting antenna is 900MHz, and a mesh-shaped steel bar structure is distributed underground. The shallow, highly energetic clutter and the reflective hyperbolic structure of the mesh reinforcement can be seen. Furthermore, the layered structure reflected by the underground interface cannot be clearly identified due to the mutual interference with the hyperbolic structure reflected by the steel bars. By using the method of the invention, the B-scan data of FIG. 6(A) is firstly subjected to high-resolution processing based on spectral whitening, and then sparse separation is carried out by using step 3. Fig. 6(B) and 6(C) are the corresponding separation results. Fig. 6(B) is a separated dot structure, which can clearly distinguish the distribution of the reinforcing bars. Fig. 6(C) is a separate layered structure, and the subsurface reflection interface is significantly enhanced in addition to the shallow strong clutter. This provides an effective treatment result for the following roadbed disease detection.
In practical application, with the increase of depth, the propagation energy of electromagnetic waves is attenuated quickly, and the frequency is lowered, so that the signal energy reflected from a lining interface is weaker, the main frequency is lower, and the difficulty in identification directly on a B-scan section by using the prior art is high. The invention carries out efficient phase shift offset imaging on the B-scan after frequency extension, so that parabolic scattering signals reflected by the steel bars are focused back to the underground position of the steel bars and are expressed as point-shaped characteristics, the spatial resolution is effectively improved, and reflected waves generated by a lining interface are still in a linear structure after being imaged. This illustrates that phase shift offset algorithms can make them appear as more distinct structural features than prior art, providing good premises for subsequent separation. The invention utilizes a morphological component analysis sparse separation method to separate the point structures such as reinforcing steel bars and the like from the layered structure of the interface, wherein the point structures are sparsely represented by two-dimensional discrete wavelet transform without downsampling, and the layered structure is sparsely represented by two-dimensional curvelet transform. The separation method of the invention utilizes the sparse representation characteristics of different structures, thereby effectively separating the point-like structure from the layered structure. Compared with the prior art 1, the data processing method can adapt to the condition of medium and non-uniform media; compared with the fixed parameter filtering method in the prior art 2, the separation method can be self-adaptive to separation data, automatically adjusts parameters, and does not need to manually set the parameters.
Finally, it should be noted that the above models and practical data examples provide further verification for the purpose, technical solution and advantages of the present invention, which only belong to the specific embodiments of the present invention, and are not used to limit the scope of the present invention, and any modification, improvement or equivalent replacement made within the spirit and principle of the present invention should be within the scope of the present invention.

Claims (9)

1. A geological radar B-scan data processing method is characterized by comprising the following steps:
s100, processing each channel of data on the B-scan section by adopting a spectral whitening frequency broadening method to obtain the B-scan section with high time resolution;
s200, performing two-dimensional offset imaging on the high-time-resolution B-scan data based on a phase shift offset imaging method to obtain a B-scan imaging section;
s300, carrying out iterative sparse separation on the B-scan imaging section obtained in the step S200 based on a morphological component analysis method;
s400, analyzing the point-like characteristics and the layered structure of the results obtained through iterative sparse separation;
the phase shift offset imaging method in step S200 includes:
s201, performing two-dimensional Fourier transform on the B-scan data after the high-time resolution processing, wherein the corresponding relation is as follows:
Figure FDA0002636698300000011
wherein f (x, z, t) & gtnon phosphorz＝0Representing the B-scan data received at the surface layer of the detection target, x representing the spatial coordinates, z representing the depth, t representing the time,
Figure FDA0002636698300000012
is an imaginary unit; f (k)x,z,ω)|z＝0Is a two-dimensional Fourier transform of the B-scan signal, kxRepresenting the spatial wave number, ω being the angular frequency;
s202, applying a phase shift operator to a wave field F (k) with the depth zxZ, ω) to obtain a fourier transform of the wavefield at a depth z + Δ z;
and S203, obtaining a phase shift imaging result of the B-scan based on two-dimensional inverse Fourier transform and making the time t equal to 0.
2. The method of claim 1, wherein step S100 is omitted if the temporal resolution of a certain piece of data of the B-scan profile is high.
3. The method according to claim 1, wherein the frequency spreading method in step S100 specifically includes the following steps:
s101, calculating one-dimensional Fourier transform of each data on a B-scan section to obtain an amplitude spectrum and a phase spectrum of each data;
s102, setting a plurality of band-pass filters in the effective bandwidth range of each channel of data, and dividing an amplitude spectrum into a plurality of sub-bands;
s103, calculating a data amplitude spectrum weighting coefficient in each sub-band;
s104, multiplying the amplitude spectrum in each sub-band by a corresponding weighting coefficient, keeping the phase spectrum unchanged, and obtaining a frequency spectrum after frequency extension;
and S105, calculating one-dimensional Fourier inverse transformation of each data on the B-scan section based on the frequency extended frequency spectrum to obtain each data on the frequency extended B-scan section.
4. The method of claim 1, wherein the iterative sparse separation in step S300 is calculated by the following relation:
Figure FDA0002636698300000021
wherein f is a B-scan section after offset imaging, f1Representing structural components of the bars and other points therein, f2Representing the horizontal or near horizontal lamellar component thereof,
Figure FDA0002636698300000022
is a two-dimensional non-down-sampling discrete wavelet transform base,
Figure FDA0002636698300000023
is f1The L1 norm of the wavelet transform coefficients,
Figure FDA0002636698300000024
is a two-dimensional curvelet transform basis,
Figure FDA0002636698300000025
is f2The L1 norm of the curvelet transform coefficients,
Figure FDA0002636698300000026
and lambda is a regularization parameter which is an error term corresponding to the error and satisfies lambda > 0.
5. The method of claim 1, wherein the sparsely separated result of step S400 includes two components: one is a section of rebar or other point-like feature;
the other is a flat layer structure, wherein the shallowest layer is a layer obtained by imaging direct waves from a transmitting antenna to a receiving antenna and reflected waves of the transmitting antenna on the surface layer of the detection target, and the deep layer is an interface between a lining and surrounding rocks or clay and other medium layers.
6. The method according to claim 3, wherein the step S103 is specifically: and calculating the amplitude root mean square value of each sub-band, and dividing a fixed constant by the root mean square value to obtain a weighting coefficient of the amplitude spectrum in the corresponding sub-band.
7. The method of claim 1, wherein: the calculation formula in step S202 is:
Figure FDA0002636698300000031
in the formula (I), the compound is shown in the specification,
Figure FDA0002636698300000032
for the phase shift operator, Δ z is the continuation step size, and θ is half of the propagation velocity of the electromagnetic wave in the medium.
8. The method of claim 1, wherein: the formula for calculating the imaging result of B-scan in step S202 is as follows
Figure FDA0002636698300000033
Then according to the conversion relation between the depth z and the propagation time t
Figure FDA0002636698300000034
Sectioning the image
Figure FDA0002636698300000035
Converted into time-space profiles
Figure FDA0002636698300000036
9. The method of claim 3, wherein: the number of the band pass filters in the step S102 is 5 to 8.
CN201810867342.1A 2018-08-01 2018-08-01 Geological radar B-scan data processing method Active CN109283523B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810867342.1A CN109283523B (en) 2018-08-01 2018-08-01 Geological radar B-scan data processing method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810867342.1A CN109283523B (en) 2018-08-01 2018-08-01 Geological radar B-scan data processing method

Publications (2)

Publication Number Publication Date
CN109283523A CN109283523A (en) 2019-01-29
CN109283523B true CN109283523B (en) 2021-04-13

Family

ID=65183364

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810867342.1A Active CN109283523B (en) 2018-08-01 2018-08-01 Geological radar B-scan data processing method

Country Status (1)

Country Link
CN (1) CN109283523B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110397446B (en) * 2019-06-26 2020-12-22 北京交通大学 Method for determining optimal detection time of tunnel lining cavity diseases
CN112014836B (en) * 2020-09-21 2022-03-04 四川长虹电器股份有限公司 Short-range personnel target tracking method based on millimeter wave radar

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102495402A (en) * 2011-12-14 2012-06-13 中国民航大学 Method for detecting and identifying road surface disaster target under interference of reinforcement high-reflection echo

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7675454B2 (en) * 2007-09-07 2010-03-09 Niitek, Inc. System, method, and computer program product providing three-dimensional visualization of ground penetrating radar data
CN102112894B (en) * 2008-08-11 2015-03-25 埃克森美孚上游研究公司 Estimation of soil properties using waveforms of seismic surface waves
CN107356967B (en) * 2017-07-26 2019-04-12 西安交通大学 A kind of compacting seismic data shields by force the sparse optimization method of interference
CN107992901A (en) * 2017-12-18 2018-05-04 武汉大学 A kind of borehole radar image rock stratum sorting technique based on textural characteristics

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102495402A (en) * 2011-12-14 2012-06-13 中国民航大学 Method for detecting and identifying road surface disaster target under interference of reinforcement high-reflection echo

Also Published As

Publication number Publication date
CN109283523A (en) 2019-01-29

Similar Documents

Publication Publication Date Title
Li et al. Denoising method of ground-penetrating radar signal based on independent component analysis with multifractal spectrum
Ni et al. Buried pipe detection by ground penetrating radar using the discrete wavelet transform
CN108037526B (en) Reverse-time migration method based on all-wave wave field VSP/RVSP seismic data
Luo et al. GPR imaging criteria
Wang et al. Noise suppressing and direct wave arrivals removal in GPR data based on Shearlet transform
CN109283523B (en) Geological radar B-scan data processing method
CN102565856B (en) Near-surface noise suppression method on basis of wave equation forward modeling
Zhang et al. Intensive interferences processing for GPR signal based on the wavelet transform and FK filtering
CN107390213B (en) A kind of time lag curve extracting method of the Ground Penetrating Radar record section based on sliding window
Aboudourib et al. A processing framework for tree-root reconstruction using ground-penetrating radar under heterogeneous soil conditions
Boiero et al. Estimating surface-wave dispersion curves from 3D seismic acquisition schemes: Part 1—1D models
JP4020365B2 (en) Frequency variable underground radar exploration method, apparatus and program
CN108254731A (en) The multiple dimensioned staged layer stripping full waveform inversion method of Coherent Noise in GPR Record
Xie et al. Identifying airport runway pavement diseases using complex signal analysis in GPR post-processing
RU2412454C2 (en) Method to process seismic data using discrete wavelet transform
Cui et al. The accurate estimation of GPR migration velocity and comparison of imaging methods
Lei et al. GPR detection localization of underground structures based on deep learning and reverse time migration
Wang et al. Automatic asphalt layer interface detection and thickness determination from ground-penetrating radar data
CN113433547A (en) Ground penetrating radar hidden crack offset imaging method, system, terminal and medium
CN108761449A (en) A kind of disaster target imaging method under reinforcement echo interference
Guan et al. Love wave full-waveform inversion for archaeogeophysics: From synthesis tests to a field case
CN109782346B (en) Acquisition footprint pressing method based on morphological component analysis
CN116540196A (en) Reinforced clutter suppression method based on distance compensation and low-rank sparse decomposition
CN109283524B (en) Method for improving geological radar signal resolution
CN110749923A (en) Deconvolution method for improving resolution based on norm equation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant