Computer Modeling of the Outgoing GPR Signal

Iskakov, Kazizat; Tokseit, Dinara; Boranbaev, Samat; Akhmetov, Iskander; Gelbukh, Irina; Iskakov, Kazizat; Tokseit, Dinara; Boranbaev, Samat; Akhmetov, Iskander; Gelbukh, Irina

doi:10.13053/cys-27-1-4543

Services on Demand

Journal

Article

Indicators

Cited by SciELO
Access statistics

Computación y Sistemas

On-line version ISSN 2007-9737Print version ISSN 1405-5546

Comp. y Sist. vol.27 n.1 Ciudad de México Jan./Mar. 2023 Epub June 16, 2023

https://doi.org/10.13053/cys-27-1-4543

Articles

Computer Modeling of the Outgoing GPR Signal

Kazizat Iskakov¹

Dinara Tokseit¹

Samat Boranbaev¹

Iskander Akhmetov²

Irina Gelbukh³^*

¹1 Lev Nikolayevich Gumilev Eurasian National University, Astana, Kazakhstan. kazizat@mail.ru, tokseit1990@gmail.com, boranbaevsa@mail.ru.

²2 Institute of Information and Computational Technologies, Almaty, Kazakhstan. iskander.akhmetov@gmail.com.

³3 Instituto Politécnico Nacional, Centro de Investigación en Computación, Mexico.

Abstract:

This paper considers a mathematical model to reconstruct the shape and tabular value of the source based on real signal data from a Loza-V series GPR receiver. A geoelectric equation in a cylindrical coordinate system is chosen as the mathematical model. Experiments, using GPR, were conducted on a homogeneous area of clean river sand, with known geoelectric properties. In this case, the equation in question is reduced to the Riccati differential equation using a special function substitution. This allowed us to obtain an explicit expression linking the spectrum function describing the response of the medium (the real radar data) and the spectrum function describing the source behavior. From the found source spectra, using inverse Fourier transforms, the emitted source itself is reconstructed in tabular form. The methodology of source reconstruction was carried out at different locations of the receiver antenna from the source antenna. In practice, geophysicists are interested in the physical characteristics of heterogeneity depending on spatial coordinates. For numerical solution of inverse coefficient problem it is necessary to have tabular value of disturbance source and tabular values of reflected signals (GPR data) at measurement points. To solve these problems we have developed an algorithm of source reconstruction and, as a consequence, determination of media response corresponding to real GPR data at the points of observation. A series of numerical calculations demonstrating the effectiveness of the considered computer model for source recovery have been carried out.

Keywords: GPR; mathematical model; Ricatti equations; inverse Fourier transform; experimental studies; radar trace spectrum; source spectrum

1 Introduction

Electromagnetic investigation methods are used: for non-destructive examination of minerals in geology; diagnostics of objects in construction; condition of highways; in archaeological and other natural science tasks. Geophysical equipment is used for experimental research: ground penetrating radar (GPR).

Theoretical provisions and practical description for solution of such a class of tasks – georadar are proposed in [¹]. GPRs come with software, the output of this software product is a radarogram.

To interpret radarograms, the essence of which is to determine the geophysical section, physics-based formulas or a fitting method are used.

This method consists in comparing the obtained radarograms with the standard views available in the database. There is a different direction of radarogram interpretation, namely computer modelling of propagation and reflection of electromagnetic waves in a medium.

The radarogram provides information about the travel time to the heterogeneity, and the task of interpreting radarograms is to determine the physical characteristics of the heterogeneity.

In GPR studies, the measurement data obtained by the GPR antenna is a response of the medium and is a function of travel time. This data is then used as additional information to solve inverse coefficient problems.

To numerically solve this kind of inverse problem, we apply an optimization method, the essence of which is minimization of the quadratic functional of inconsistency of calculated and observed fields (data from the receiving antenna of GPR).

In computer simulations, to solve the inverse coefficient problem, a tabular value of the perturbation source as well as tabular values of the reflected signals at the measurement points are needed.

In order to solve these issues, we have developed a source reconstruction algorithm in this paper. The physical characteristics of the investigated objects include: dielectric and magnetic permeability and conductivity of media.

The issues of numerical method for solving inverse coefficient problems for the geoelectric equation are discussed in the monograph by S. I. Kabanikhin [¹⁰].

Application of optimization methods for numerical solution of a class of coefficient inverse problems is presented in the monograph by K.T. Iskakov et al. [⁸].

The use of engineering methods for interpreting radarograms, the essence of which consists in determining the geoelectric section: dielectric and magnetic permeability; conductivity and depth of heterogeneity are considered in a series of works [²⁴, ²⁰, ²³, ¹²].

The method of layer-by-layer recalculation, for inverse problems in the case of layered media, is widely used.

The idea of layer-by-layer recalculation is that the differential equation describing this process is reduced to Riccati differential equation by special function substitution, for which the solution can be written out in analytical form [¹⁵].

The theoretical foundations and practical applications of subsurface georadolocation are outlined in [², ¹⁶, ⁶, ³, ⁹, ⁴].

The theory of uncorrected and inverse problems received a rapid development in the 20th century and goes back to the first works in this direction by academician A.N. Tikhonov [¹⁸].

The theoretical aspects and numerical methods for solving inverse coefficient problems for the geoelectric equation are devoted to a monograph by V.G. Romanov and S.I. Kabanikhin [¹⁷].

The application of optimization methods for the numerical solution of the coefficient inverse and its theoretical justifications are given in the monograph by S.I. Kabanikhin and K.T. Iskakov. [¹¹].

One of the first algorithms of layer-by-layer recalculation for solving a second-order differential equation for a horizontally layered homogeneous medium is described in the work by A.N. Tikhonov and D.N. Shakhsuvarov. [¹⁹].

For solving the practical task of electrical prospecting, this algorithm was used in the work by Dmitriev V.I. [⁵]. In works of Karchevsky A.L. [¹³, ¹⁴] described an algorithm of solution of seismic task for a horizontally layered medium, and in [14] this technique is applied to Maxwell equations in case of horizontally layered medium of any anisotropy. The above-mentioned studies of direct and inverse problems of electrical prospecting refer to a series of problems when observations are carried out on the daytime surface of plane-parallel media.

If the surface has a relief, for example, in archaeological problems of investigation of barrows, then considering the influence of relief of the surface on the electric field and on the results of 2D and 3D inversion of electric tomography (ET) survey data is an important task.

In [²⁵] the effectiveness of elimination of topographical effects in 2D inversion programs that include topography in the grid in the presence of daylight surface topography has been studied. Modeling of the apparent resistivity curves has been performed by the method of integral equations.

Experimental studies were carried out on the laboratory polygon using Loza-V GPR. Real data are used for numerical calculations to determine the shape and tabular value of the emitted signal.

The work also uses the methodology for signal cleaning from noise and interference and signal recovery given in papers [⁷, ²²].

2 Mathematical Model for Source Recovery

The propagation of electromagnetic waves in a medium is described by a system of Maxwell’s equations [¹⁷].

With a special choice of a perturbation source, namely if the source is a long cable located along the coordinate axis and assuming that the functions describing the geoelectric section depend on one variable, namely the depth, then the system of Maxwell equations is simplified and we have a formulation problem for the geoelectric equation [¹⁷].

This model is suitable for the study of direct and inverse problems in the case of slime media.

Let us consider the geoelectric equation in cylindrical coordinate system:

εwtt+σwt=1μ(wrr+1twr+wzz)+1rδ(r)δ(z−z∗)q(t). (1)

The initial conditions are of the form:

w(r,z,0)=0, wt(r,z,0)=0. (2)

Assume that ∂w∂r is bounded at r→0. Let the following relation take place:

u(ξ,z,ω)=∫−∞∞eiωt∫0∞w(r,z,t)rJ0(ξr)drdt. (3)

Then, task (1)-(2) will take the form:

uzz−b2(z) u=δ(z−z∗)q^(ω), (4)

[uz]z=0=0=0, [u]z=0=0. (5)

Let us use the definition of the generalized derivative: uz={uz}+[u]z∗δ(z−z∗), in order to transfer the right-hand side of equation (4) to the boundary condition, then finally write problem (4)-(5), differently:

uzz−b2(z)u=0, (6)

[uz]z∗=0=0, [u]z∗=0=0, (7)

u[z]z=z∗=0. (8)

Using standard techniques, by analogy as in [¹³, ¹⁴], let us go to the Riccati equation, in order to derive an analytical solution, we have:

uz=yu⇒y′+y2=b2, (9)

z∈(−∞,z∗), z∈[0,∞), (10)

y(z)=b0, (11)

y(z)=−b1=y0. (12)

Let us calculate:

y∗=b0(y0+b0)e2b0z∗+(y0−b0)(y0+b0)e2b0z∗−(y0−b0), (13)

u|z=z∗+0y∗−u|z=z∗−0b0=q^(ω), (14)

u|z=z∗+0−u|z=z∗−0=0, (15)

u∗=q^(ω)y∗−b0, (16)

u(0)=2b0eb0z∗(y0+b0)e2b0z∗−(y0−b0)∗q^(ω)y∗−b0. (17)

Let:

z∗→0⇒y∗=y0=−b1, (18)

u(ξ,0,ω)=q^(ω)b1+b0. (19)

And we finally get it:

u(r0,0,ω)=q^(ω) ∫0∞ξJ0(ξr0)ξ2−(ω2μ0ε1−iωμ0σ1)+ξ2−ω2μ0ε0dξ. (20)

Thus, gathering all the calculations, let’s describe the solution scheme for source recovery in steps:

1. Let GPR data be known - response of the medium in a homogeneous medium: w(r0,0,t) - real radar data at the measurement point, where r₀is the distance of the antenna location from the source. The source power is on the order of 200 kHz.
2. Let us determine the spectrum, function w(r0,0,t) based on the Fourier transform, we obtain S(r0,0,ω). Solving the Riccati equation analytically, we have:

S(r0,0,ω)=q^(ω)∫0∞ξJ0(ξr0)ξ2−(ω2μ0ε1−iωμ0σ1)+Rdξ, (21)

where: R=ξ2−ω2μ0ε0.

3. By numerically calculating the integral on the right-hand side of formula (2), we obtain that the value of the source in the frequency domain is finally calculated according to the formula:

q^(ω)=−S(r0,0,ω)I, (22)

where I is the approximate value of the integral.

4. The source value in the frequency domain q^(ω) (22) consists of technical and computational errors that result from the averaging of the signal during gating and the conversion of the signal from analog to digital.
To clean the signals from noise and interference, wavelet transform and Harr, Dobesh filters are applied by analogy as in [⁷]. Translation of signals from binary to another format.
5. We construct a calibration function F^(ω;α,ω0), using the theory of uniform function approximation.
6. Let the calibration function F^(ω;α,ω0) have the form:

F^(t)=e−αtcos⁡ (ω0t). (23)

7. Based on the Fourier transform, we get:

F^(ω;α,ω0)=α−i(ω−ω0)α2+(ω−ω0)2. (24)

8. To find the unknown parameters α, ω0 consider the quadratic deviation:

ψ=∑ω|q^(ω)−F^(ω;α,ω0)|2, ω¯≤ω≤ω¯, (25)

where the values of ω¯, ω¯ depend on the radar characteristics.

9. Calculating the derivatives needed to minimise the quadratic deviation (25), we obtain:

ψ=∑ω[q^r+iq^i−F^]2=∑ω{[q^r−αα2+(ω−ω0)2]2+[q^i+ω−ω0α2+(ω−ω0)2]2} (26)

Let us define the following variables in order to make the following expressions easier to understand:

A=q^r−αα2−(ω−ω0)2, (27)

B=(−1)(α2+(ω−ω0)2)−2α2(α2−(ω−ω0)2)2, (28)

C=−(ω−ω0)2α[α2+(ω−ω0)2]2, (29)

D=α[2(ω−ω0)](α2+(ω−ω0)2)2, (30)

K=q^i+ω−ω0α2+(ω−ω0)2, (31)

N=(−1)[α2+(ω−ω0)2]+(ω−ω0)2[α2+(ω−ω0)2]2. (32)

With these definitions in place, we can proceed with the following:

∂ψ∂α=∑{2AB+2[+ω−ω0α2+(ω−ω0)2]C}, (33)

∂ψ∂ω0=∑ω{2[A]D+2[K]N}. (34)

Apply the method of conjugate gradients [¹⁷]:

[αω0]k+1=[αω0]k−γk Pk. (35)

Here Pk=∇ψ[αk,ω0k]−βkPk−1:

∇ψ=[∂ψ∂α∂ψ∂ω0], (36)

βk=−‖∇ψ[αk,ω0k]‖2‖∇ψ[αk−1,ω0k−1]‖2. (37)

10. The next approximations αk+1,ω0k+1 will be calculated using conjugate gradient formulas (35)-(37).
11. As a result of the calibration, the final original source is as follows:

F(t)=e−α∗tcos⁡(ω0∗t), (38)

where: α∗,ω0∗ are the found values for which the functional (25), reached a minimum.

The found value of the source, can be used to solve the forward problem, and the inverse problem of determining the dielectric permittivity and conductivity of the media.

3 Experimental Results

Experimental studies on reception of reflected signals from a homogeneous site (river sand) were carried out. The experiment was carried out with a Loza-V series GPR, with antenna sweep: 100 cm.

Research task: geophysical investigation of homogeneous clean sand underlayer structure: simulation of pulse source from “Loza-V” series device; determination of spectral characteristics of signals emitted by the antenna.

Calculation of the source signal in tabular form, are necessary in the future for solving inverse coefficient problems in determining the geoelectric section. A homogeneous section, river sand, measuring 8 metres in width and in length, was chosen for the experiments.

An experiment was carried out in which the GPR source was placed at point B0 and the antenna was placed at a distance of 1 metre and is shown in Fig. 1. Trace plots, i.e. the response of the medium is shown in Fig. 2.

Fig. 1 Measurement scheme

Fig. 2 Radarogram trace plot. (The antenna is positioned from the source at a distance of - 1 metre)

The actual radarogram data, which has its own format, has been converted into a text format of tabulated signal values. Signals received from Loza-V series GPR are output in ”geo” format. For signal analysis and visualization we convert from the binary format of the file ”geo” in the format ”txt”.

File ”txt” consists of three columns: x, t, alg. Where: x - number of data received from each observation point, by default x=0 ... 10; T -time, T0=0 ns, T512=255,5 ns, step 0.5 ns; Alg -amplitude values.

To encode the amplitude value of an analogue signal, an 8 or 16 bit representation of the amplitude values is used.

In the case of 8-bit encoding, the amplitude measurements of the analogue signal will be obtained with an accuracy of 1/256 of the dynamic range of the digital device, since 8 bits allow 28 numbers to be represented-256).

The digital signal in this encoding is then a set of numbers from 0 to 255 (or -128 to 127). This accuracy is not sufficient to reconstruct the original signal because there will be errors of non-linear distortion.

If you encode the amplitude of the analogue signal to 16 bits, you will have 265 times the accuracy. Digit capacity of 16 bits allows to encode 216=65536 values of amplitude.

A digital signal is then a collection of numbers from 0 to 65535 (or -32768 to 32767). This coding approach keeps nonlinear distortion to a minimum. A digital signal is a dimensionless set of numbers and has no common units of measurement.

The radarogram data, which has its own format, is converted into a text format of tabulated signal values. Fig. 2 shows a graphical representation of the response.

4 Related Work

The spectrum of the radarogram trace (medium response), is shown in Fig. 3. The frequency scale of the graph varies from 0 to 1000MHz.

Fig. 3 Spectrum of the radarogram trace. (The antenna is located at a distance of 1 metre)

The modulus of the frequency distribution of the signal u(r0,0,ω) can be seen in Fig. 3. Having the radarogram trace spectra, using equation (2.)

u(r0,0,ωi)=q^(ωi)∫0∞ξJ0(ξr0)ξ2−(ωiμ0ε1−iωiμ0σ1)+ξ2−ωi2μ0ε0dξ. (39)

For i=1,2,3,…,n/2. Let us represent the integral of the complex function in equation (39) as:

∫0∞ξJ0(ξr0)ξ2−(ω2μ0ε1−iωμ0σ1)+ξ2−ω2μ0ε0dξ=∫0∞(p(ξ)+iq(ξ))dξ∫0∞f(ξ)dξ. (40)

The numerical calculation of the integral (40) is realised by the iterative rectangular (or trapezoidal) method:

∫f(ξ)bξ≈h∑j=0∞f(ξj). (41)

The summation is carried out to a given accuracy ε:

|f(ξj)|<ε.

The following parameters were used in the calculation: ε0=1, ε1=5, μ0=1, r0=1, σ1=1/500. The graph of the source spectrum is shown in Fig 4.

Fig. 4 Source spectrum graph. (The antenna is located at a distance of - 1 metre)

By inverse Fourier transform of the spectrum q(ω) we find the absolute (complex) source signal. The plot of the real part of the source signal, in the case of our experiment is shown in Fig. 5.

Fig. 5 Graph of the recovered signal source, from experiment

Fig. 5 below shows a graph of the reconstructed source, when the antenna is located at a distance of 1 metre from the source.

A certificate of authorship No 31827 dated ”17” January 2023 [²¹] was received for GPR signal source detection program.

5 Conclusion

We considered a computer model to reconstruct the shape and tabular value of the source on the basis of real data of signals from the Loza-V series GPR receiver.

As a result of strictly justified mathematical calculations, we obtained an explicit expression linking the spectrum function describing the response of the medium (real radar data) and the spectrum function describing the source behavior.

Based on the found source spectrum on the basis of inverse Fourier transform, the emitted GPR source itself is reconstructed in tabular form and its graphical representation.

The results of researches of this article can be applied to numerical solution of inverse coefficient problems, as it is necessary to have a tabular value of the source of disturbance, as well as tabular values of reflected signals (GPR data), in the points of measurements.

A series of numerical calculations were carried out to show the effectiveness of the considered computer model on source recovery.

References

1. Aleksandrov, P. N. (2017). Theoretical foundations of the GPR method. Editorial de Literatura de Física y Matemáticas de Moscú. [ Links ]

2. Andriyanov, A. V. (2005). Issues of subsurface radiolocation. Collective monograph, Moscow: Radiotekhnika, pp. 416. [ Links ]

3. Conyers, L. B., Goodman, D. (1997). Ground-penetrating radar: An introduction for archaeologists. AltaMira press. [ Links ]

4. Daniels, D. J. (2004). Ground penetrating radar. The Institution of Electrical Engineers, London, UK. DOI: 10.1049/PBRA015E. [ Links ]

5. Dmitriev, V. (1968). General method of electromagnetic field calculation in layered medium. Computational Methods and Programming, Vol. 10, pp. 55–65. [ Links ]

6. Finkelstein, M. I., Karpukhin, V. I., Kutev, V. A., Metelkin, V. N. (2017). Subsurface radiolocation. Moscow, Radio and Communications, pp. 216. [ Links ]

7. Iskakov, K. T., Boranbaev, S. A., Uzakkyzy, N. (2017). Wavelet processing and filtering of the radargram trace. Eurasian Journal of Mathematical and Computer Applications, Vol. 5, No. 4, pp. 43–54. DOI: 10.32523/2306-3172-2017-5-4-43-54. [ Links ]

8. Iskakov, K. T., Romanov, V. G., Karchevsky, A. L., Oralbekova, J. O. (2014). Research of inverse problems for differential equations and numerical methods of their solution. L.N. Gumilev Eurasian National University Press. [ Links ]

9. Jol, H. M. (2009). Ground penetrating radar. Theory and applications, Elsevier sciences, pp. 524. [ Links ]

10. Kabanikhin, S. I. (2018). Inverse and uncorrected problems. Siberian Scientific Publishing House, pp. 511. [ Links ]

11. Kabanikhin, S. I., Iskakov, K. T. (2001). Optimization method of the solution of coefficient inverse problems. Novosibirsk State University. [ Links ]

12. Kabanikhin, S. I., Iskakov, K. T., Sholpanbaev, B., Shishlenin, M. A., Tokseit, D. K. (2018). Development of a mathematical model for signal processing using laboratory data. Bulletin of the Karaganda university-mathematics, Vol. 92, No. 4, pp. 148–157. [ Links ]

13. Karchevsky, A. L. (2005). The direct dynamic seismic problem for horizontally layered media. Siberian Electronic Mathematical Proceedings, Vol. 2, pp. 23–61. [ Links ]

14. Karchevsky, A. L. (2007). Analytical solution of maxwell equations in the frequency domain for horizontally layered anisotropic media. Geology and Geophysics, Vol. 48, No. 8, pp. 689–695. DOI: 10.1016/j.rgg.2006.08.005. [ Links ]

15. Karchevsky, L., Andrey (2020). Analytical solutions to the differential equation of transverse vibrations of a piecewise homogeneous beam in the frequency domain for the boundary conditions of various types. Journal of Applied and Industrial Mathematics, Vol. 14, pp. 648–665. DOI: 10.1134/S1990478920040043. [ Links ]

16. Kopeikin, V. V. (2007). Wave refraction in linear media with frequency dispersion. M:Nauka. [ Links ]

17. Romanov, V. G., Kabanikhin, S. I. (1991). Inverse Problems of Geoelectrics. Nauka. [ Links ]

18. Tikhonov, A. N., Arsenin, V. I. (1977). Solutions of Ill-posed Problems. Wiley. [ Links ]

19. Tikhonov, A. N., Shakhsuvarov, D. N. (1956). Method for calculation of electromagnetic fields excited by alternating current in layered media. Izv. AN SSSR, Geofizika, Vol. 3, pp. 251–254. [ Links ]

20. Tokseit, D. K., Boranbaev, S. A., Oralbekova, J., Nurzhanova, A. B. (2020). Determination of the geoelectric section by georadar data. Bulletin of D.Serikbayev East Kazakhstan State Technical University, Vol. 89, No. 3, pp. 154–160. [ Links ]

21. Tokseit, D. K., Iskakov, K., Boranbayev, S. A. (2023). GPR signal source identification software. Certificate of Record in the State Register of Rights to objects protected by copyright, No 9319, pp. 1. [ Links ]

22. Tokseit, D. K., Iskakov, K. T., Boranbaev, S. A. (2022). Techniques for processing source signals emitted by GPR. Universitet Enbekteri – University Proceedings, Vol. 86, No. 1, pp. 323–333. [ Links ]

23. Tokseit, D. K., Iskakov, K. T., Boranbayev, S. A., Shishlenin, M. A. (2020). Interpretation of radarograms of geological section based on experimental calculation formulas. Certificate of Record in the State Register of Rights to objects protected by copyright, No 9319, pp. 1. [ Links ]

24. Tokseit, D. K., Iskakov, K. T., Kabanikhin, S. I., Shishlenin, M. A. (2020). Program for calculating a mathematical model for determining the response of the medium and detecting heterogeneity. Certificate of Record in the State Register of Rights to objects protected by copyright, No 12894, pp. 1. [ Links ]

25. Turarova, M. K., Mirgalikyzy, T., Mukanova, B., Modin, I. N. (2022). Elimination of the ground surface topographic effect in the 2D inversion results of electrical resistivity tomography data. Eurasian Journal of Mathematical and Computer Applications, Vol. 10, No. 3, pp. 84–104. [ Links ]

Received: November 08, 2022; Accepted: January 12, 2023

^* Corresponding author: Irina Gelbukh, e-mail: ir.gelbukh@gmail.com

This is an open-access article distributed under the terms of the Creative Commons Attribution License