Ground Penetration Radar is a controlled source geophysical method which uses high frequency electromagnetic waves to study shallow layers. The use of this method began in 1956 and has been developed since the 1970s. GPR system has been commercially available since the 1980s, and has been widely used since the mid-decade. In the GPR method, high frequency electromagnetic waves (from 12.5 to 2300 MHz) transmit to the earth and these waves are reflected by objects or clear boundaries between the underground layers. Electromagnetic reflections are created by the difference in electrical conductivity (dielectric constant) between materials that pass electromagnetic waves. GPR electromagnetic waves pass through materials that are electrically conductive, but they are heavily absorbed when passing through high-conductivity materials such as clays, organic acid soils and saturated saline water  .
Resolution of GPR varies from depth of centimeters to a several meters with maximum depth of about 100 meters. Resolution of this method depends on difference of electrical properties between target and surrounding electrical medium, target geometry and used bandwidth. Resolution of this method can be high which can recognized fine layers in near surface structures and objects in the ground can be well  .
The wavelet transform is used extensively in signal analysis and noise attenuation. In addition, wavelet domain allows local precise descriptions of signal behavior. The Fourier coefficient represents a component for all time and therefore local events must be described by the phase characteristic which can be abolished or strengthened over a large period of time. Wavelet expansion coefficients represents a component which is local itself and convenient for interpretation. The wavelet Transform may allow separation of signal components which overlaps in time and frequency. Wavelets can be designed to fit the unique applications because they are also customizable and tunable and there is not just one wavelet. They are ideal for adaptive devices that adjust themselves to the signal. Finally, production of wavelets and the calculation of discrete wavelet transforms (DWT) are well suited to digital computing  .
Finally basis of Auto Regression (AR) is the fitting of an appropriate model on data, which in practice results in more information from data process. Estimation of the parameters of the regression model (AR) is very important. Also use of a suitable model is very important in approximation of random signals spectrum, especially in parametric methods. Parametric methods of power spectrum estimation are based on selecting the appropriate model for the data  .
In order to obtain a higher-resolution spectral estimation than other models, recursive operator is a suitable tool. Generally, it is much easier to work with an Auto Regression model. For this reason, suitable techniques are presented for estimating of model parameters  .
2. Problem Definition and Research Method
GPR signal in receiver X(t) can be formulated as equation (usually infected with noises):
In Equation (1), the earth reflection series s(t) is matched in the wavelet due to the spring w(t) and the noise n(t) is added to the data, which must be corrected. Since it is impossible to eliminate all of the noise, then the denoising is done with the aim of obtaining X(t) is as close as possible to s(t)  .
In this research, noise attenuation method of auto regressive on wavelet and f-x domains was used to reduce the noise of the GPR signal. First AR operator and definition of f-x domain will illustrate and explained using some examples. Ultimately, all methods will be applied to real and synthetic data of the GPR signal and their results will be compared together.
3. Wavelet Transform with Tunable Quality Factor (TQWT)
In this section, a wavelet transform method is described where it easily uses the Q factor and wavelet transform process can be controlled using the periodic nature of the signal. The function of scale has real values. The Q factor is a ratio of the frequency center to signal bandwidth. Most wavelet transform methods have a poor ability to control behavior of signal frequencies. The TQWT method is completely discrete, it has a complete signal reconstruction ability and mean redundancy.
3.1. Low Pass Scale
The low-frequency region of the signal is analyzed using the low-pass function. If the scale parameter is chosen for the low-pass scale function part, the sampling changes to , where is sampling frequency of the input signal. If , then the function of the low-pass signal changes:
And if, low pass function changes a signal like:
Note that low-pass function keeps signal behavior around a zero frequency according to Figure 1.
3.2. High Pass Scale
The high-frequency signal region is analyzed by high-frequency scale function. If the scale parameter is selected for the high-pass section function, the sampling is changed to , which is sampling frequency of input signal. If , the of the high-pass scale function will change such a signal:
And if , high pass scale function changes like:
Figure 1. Low-pass function keeps signal behavior around the zero frequency  . (a) Low-pass scaling block diagram; (b) low-pass scaling with α < 1; (c) low-pass scaling with α > 1.
The high-pass scale function keeps signal around the Nyquist frequency according to Figure 2.
Using a frequent dual channel filter, the TQWT method is applied. This transform is represented by three steps in Figure 3.
Using a frequent dual channel filter, the TQWT method is applied. This transform is represented by three steps in Figure 3  .
denotes the subband in jth step and it’s generated for j = 1 high-pass subband in the first stage. The subband bandwidth of the jth frequency is Also is a sampling frequency of input data.
Figure 4 shows frequency response for for 4 values of . By changing these variables, spectral frequency decomposition can be controlled.
Figure 2. The high-pass scale function keeps signal around the nyquist frequency  . (a) High-pass scaling block diagram; (b) High-pass scaling with β < 1; (c) High-pass scaling with β > 1.
Figure 3. Using a frequent dual channel filter, the TQWT method is applied. This transform is represented by three steps in figure  .
Figure 4. Frequency response H_1^j (ω) for 1 ≤ j ≤ J for 4 values of (β, α)  .
The dual-channel filter in Figure 3 with the factor is redundant. If this filter is repeated in the low-pass, the Redundancy factor is:
This relation is deduced from the fact that a sampling frequency under the jth subband for is determined by , So sampling frequency is equivalent to (6).
3.4. Central Frequency
The frequency response in the jth step is nonzero at the distance .
The central frequency of the jth stage is approximately average .
This relationship is based on the sampling frequency:
The bandwidth of the frequency response under the jth subband is about a half of bandwidths of all frequencies that have a non-zero response.
If the Q factor is expressed on the basis of
Note that Q is not dependent on the stage. A Transform method is Q constant and depends only on filter bank parameters. To determine α and β, they are rewritten based on r and Q:
A Q factor must be greater or equal to 1, and if Q is equal to 1, wavelet transform builds a second derivative of a Gaussian function. If Q is greater than one that means signal has a higher periodic behavior and the Oversampling rate should be greater than one. If r is near 1, the transition region is very thin and the time response of the wavelet is not properly determined. If makes a sinc wavelet, then it’s therefore optimal to .
When the signal is staircase oscillatory, a number of vanishing moments of the wavelet transform is considered. This is not the case when we talk about a periodic signal like a GPR signal. The number of steps for this type of wavelet transform is limited by the length of the signal. After a certain stage, the signal will be very short to become two samples. If countless steps are taken, conversion of these steps is difficult to reconstruct. J must be determined so that the wavelet does not extend beyond the length of the category.
4. F-X Domain
Generally, a signal in time domain is seen as a mixed signal in frequency domain. From frequency domain, information such as phase spectrum, the amplitude spectrum and energy spectrum can be obtained. The energy distribution is also extracted in terms of frequency  .
Given that in time domain, input data or signal, for each frequency is shown as sinusoidal functions in the location, the separation of the signal from the noise is made easier. Therefore, we will investigate predictive filter in this domain. This method assumes that a class is composed of delayed impacts in accordance with Figure 5.
The f-x prediction filter method is very useful method to attenuate random noise. 1. Selesnick, (I.W. Sparse signal representations using the tunable Q-factor wavelet transform. in Wavelets and Sparsity XIV. 2011. International Society for Optics and Photonics.)
In this method, it is assumed that trends are linear. The seismic section can be divided into smaller windows in cases where trends are not linear, in order to satisfy the assumption given that trends in that window are linear. Assume that a seismic is a sequence of impacts with different domains defined as  :
where t is the time, x is the lateral position, is the jth impact, and is the delay function that represents the event’s form. Considering Fourier transform of (15) we have:
Since exponential functions can be written as a set of sinuses and cosines, in (16), the hypothesized seismic traces are converted into a set of sinuses and cosines as functions of ω and x. Since the f-x filter only predicts linear data, events should be linear to x. So the function must be linear  . Given the linearity of in the frequency domain, (16) will be:
Figure 5. A schematic diagram for delayed traces  .
is a conjugate constant, which depends on the source power and the reflection coefficient, and is the gradient of the linear event. (17) Shows that by linearization assuming, the function is completely sinusoidally in x, which means that this function is predictable and predictable. Therefore, the signal is a complex exponential function in terms of x in the time-frequency domain and it will be predictable  . The Wiener filter is an optimal stable linear filter for noise-damaged images. The Wiener filter needs to assume that the signal and noise are secondarily stable. For this purpose, noises are considered frames with a mean zero. The procedure used in this method is to minimize the sum of squares of the difference between output signal and the actual output signal from the filter. The Wiener filter is incapable of reconstructing the frequency components that are damaged by noise, and can only eliminate them. This filter is incapable of neutralizing the image of the image and has a low velocity. To speed up this filter, we can use this filter to obtain a shock response, reversed the fast Fourier transform  . To speed up this filter, we can use inverse fast Fourier transform and obtain the impulse response.
5. Autoregressive Vector
Classical models of time series are divided into stationary and non-stationary sections. Autoregressive vector is part of classic static time series model that we will discuss in this article. One of the limitations of our models is that they impose a one-way relationship whose predictor variables are influenced by predicted variables which this process must be reversed. However, in many cases, the reversal of this action must also be made where the variables affect each other. In this framework, all variables behave symmetrically  . In an Autoregressive model, the variables are predicted using a linear combination of the previous values of the variables. The term Autoregressive represents a regression of the variable against itself  . The regression model is defined as follows:
where c is constant and is a white noise with mean zero and the variables and are the parameters of the model. In this case, is called the p-th Autoregressive model with and is represented by AR (p). An Autoregressive structure is simple, useful, and easy to understand in a wide range of fields. The first-order Autoregressive is defined as follows  .
Firstly, the residue is expected to be zero
And secondly, errors are not self-correlated.
Automatic regression models are considerably flexible in controlling a wide range of different time series. Figure 6 shows main steps of Autoregressive vector.
Autoregressive models do not allow us to explain on causal relationships. This is especially applies when Autoregressive models are merely generalized for the processing of unknown time series. While a causal interpretation requires a basic cost model. However Autoregression provides a proper interpretation of the variables shown  .
In order to increase the signal-to-noise ratio of a multi-point signal (S/N), firstly the Autoregression vector of the noisy data is estimated and shown as a . Then, the forward estimated of denoised data is obtained by the following equation  .
Similarly, the estimated inverse equation of de-noised data are obtained:
represents the complex conjugate of Autoregressive vector operator. The final equation of the de-noised data are obtained by mean for forward and inverse estimates:
Figure 6. The AutoRegressive vector analysis algorithm  .
6. Applying Autoregressive Filter in f-x and Tunable Quality Wavelet Transform (TQWT) Domain on GPR Data
Given that the purpose here is to reduce the random noise, in the first step a random noise should be added to the section which obtained from synthetic models. For synthetic data generation, the ground model is firstly considered with arbitrary coefficients. Given the similarity in wave propagation of electromagnetic waves and seismic waves, synthetic modeling of GPR and seismic data are analogous. Random noise is added to synthetic model and various de-noising methods which presented in this article apply to it. After denoising of synthetic model, real GPR data will be denoised. Here, by applying linear Autoregression modeling on wavelet and f-x domain, random noise in actual data will be weakened. The procedure is applying these steps on real data in wavelet domain and then transforms the denoised data waveform from wavelet domain. After these steps and in next section, results will be compared and evaluated.
6.1. Wavelet Function of Synthetic Data Generator
The wavelet used to generate this synthetic data is the Ricker Wavelet. This wavelet is formulated as  :
is Ricker wavelet, t is time and f is the central frequency. Ricker wavelet is symmetric in time domain and has zero mean ( ).
The generated synthetic signal is formulated as:
In (26), is a received that is correlated with the Ricker wavelet (as a synthetic wavelet source). is a random noise that is added to the input signal.
6.2. Synthetic Model Denoising
6.2.1. TQWT Denoising in Soft and Hard Thresholding
To evaluate efficiency of mentioned denoising methods, an artificial section considered which is consisting of five layers and three pipes according to Figure 7. Based on that section, GPR data will be generated as Figure 8-1 and random noise will be added according to Figure 8-2. The sampling frequency is 40 ns and the noise is the 40 dB white Gaussian noise. As shown in Figure 8-2, the path of the layers are somewhat lost particularly thinner layers. Also, the prominence and boundary of layers and pipes throughout the noisy cross-section have been damaged.
Figure 8-1 shows the GPR cross section in the time domain and in Figure 8-2 represents the noisy section which is noise totally dominated and covers the tremors. Figure 8-3 shows the application of the TQWT filter with a soft thresholding method and Figure 8-4 after the hard thresholding. Figure 8-4 shows that this filter causes a slight softness to be signaled but has not been successful in noise attenuation and that noises are still present.
The only difference between noisy and the hard thresholding denoised signals are in a deformation of the noises. The noises are somewhat diminished, but it shows another form of noise as an elongation along the axis of time. The boundaries of layers and hyperbolas are still confusing and effect of noise is not diminished. Figure 8-3 shows application of the TQWT with a soft thresholding method. This method has been far more successful in noise attenuation. A significant part of the noise has been attenuated, and boundaries of layers and hyperbolas are significantly clear.
Figure 7. Synthetic model which consists of five layers and two shapes of a pipe and a duct. Their conductivity and dielectric constant are shown.
Figure 8. 1―Synthetic data, 2―Synthetic data with white noise, 3―Denoised data by soft thresholding TQWT, 4―Denoised data by Hard Thresholding TQWT, 5―Denoised by Autoregressive-FX.
High-frequency noises have been far more attenuated, but the elongation portion is also included in this filter. But this elongation is significantly less than the hard thresholding. Therefore, signal is significantly improved by soft thresholding method compared to the hard thresholding.
6.2.2. Investigation of TQWT Power Spectrum Results in Soft and Hard Thresholds
Figure 9 shows the results of the power spectrum. Soft thresholding power spectrum results (blue) at low to moderate frequencies are the most reduced in comparison to the power spectrum of input data (red). Its reduction compared to the hard thresholding (black) is remarkable. The gradient of soft thresholding power spectrum is also confirmed with the power spectrum of signal which indicates the sustaining of signal in this method. The highest level of noise attenuation in soft thresholding was also achieved with random noise at high frequencies. But the point is that after a certain frequency (2.10 MHz) and at high frequencies, noise gain is quite clear in both the soft and hard thresholding.
6.2.3. Autoregressive-FX Filter
To denoising in this section, Autoregressive filtering is applied in time domain. Figure 8-5 shows data after denoising by Autoregressive model in the time domain. After applying the Autoregressive filter is time domain, it is observed that the noise has been attenuated and the boundaries of the layers are specified much well. In Figure 8-5 after applying the Autoregressive filter, layer boundaries are well detectable, softened, and the noise is well reduced. The point is elongation in direction of axis of time is happening due to the application of the filter. But in general, it has successfully restored the boundary between layers and parabolas.
7. Real Data Denoising
In this section the filtering on real data will be discussed. Firstly the Autoregressive-FX filter and after that the TQWT in soft and hard thresholding ways.
7.1. Autoregressive-FX Filter
With respect to Figure 10, it can be concluded that the AR filter has been relatively successful in noise attenuation. Layer and hyperbolic boundaries have been improved and some of the high-frequency noise has been eliminated. It can be said that the signal is denoised, but some of the high-frequency noise is still dominant. As we said before, our noise is random does not have a definite origin. In Figure 10, we see that the right side contains a lot of noise, especially at high arrivals and in most parts of the data it is completely polluted by noise which that part of the traces are covered.
7.2. Real Data Noise Attenuation Using TQWT
Figure 11-1 shows the main signal and the Figure 11-2 and Figure 11-3 shows TQWT filters applied to signal by soft and hard thresholding. In Figure 11-2 and Figure 11-3, it is observed that in hard thresholding method, significant amounts of high-frequency noise are observed. Also, the high-gain times are totally infected with high-frequency noise and no noise reduction has been made. Also high arrivals in traces are completely infected with high-frequency noise and noise attenuation is not successful. But the soft thresholding completely attenuated the noise and like raw synthetic data, the Background is completely out of high-frequency noise. Of course, part of the elongations marks still shows up, but the boundary of the deep layers is well recovered. Reflection of the hyperbolas and their interference points are clear and none of them have been removed. So that it can be said that the noise in this section of the signal is completely eliminated and the noise attenuation is very well done. This results are confirmed by power spectrum of real data by TQWT hard and soft thresholding in Figure 12.
Figure 9. Power spectrum in terms of synthetic data frequency. 1―Red: raw input signal. 2―Black: Input noisy data. 3―Blue: Denoised output data using soft thresholding TQWT. 4―Green: Output data denoised by TQWT Hard Thresholding. 5―Dashed Brown: Denoised synthetic data using Autoregressive-FX filter.
Figure 10. Real data noise attenuation using Autoregressive-FX.
Figure 11. 1―Real data, 2―Denoised by soft thresholding TQWT, 3―Denoised by Hard thresholding TQWT.
Figure 12. Power spectrum of real data. 1―Red: Raw signal, 2―Black: Denoised by soft thresholding TQWT, 3―Denoised by Hard thresholding TQWT.
As stated several factors make the GPR data noisy such as mobile networks, power stations, etc. Considering that the increase of signal-to-noise ratio is the main objective of this research a method was chosen that to attenuating the noise and minimizes the damage to the signal. In this paper synthetic and real GPR data was denoised with the TQWT in two soft and hard thresholding and the Autoregressive-FX filter. According to the results, it can be concluded that the hard thresholding method is not suitable. Results do not show a good improvement, especially at low frequencies, most of the noise remains, and the elongations which caused by hard threshold has also been added to these noises.
However, in the soft thresholding problem of elongation noise is observed. But a significant part of the noise is attenuated and high-frequency noise has been eliminated. After the deactivation, the spikes from the real signal are well observed. Also boundaries of layers and hyperbolas are also returned. The Autoregressive-FX method has succeeded to denoising of high-frequency noise. The layer boundaries are also clear, but less denoising has been achieved compared to the soft thresholding TQWT. Also layer boundaries are less detected and part of high-frequency noise is still seen.
As a result, the wavelet domain provides a suitable frequency resolution for noise attenuation. But the type of filter applied to the wavelet transform also affects denoising. The TQWT method has a good frequency resolution because it operates in the wavelet domain, but the results of hard thresholding indicate that the selection of the appropriate filter for de-noising substantially changes results. Finally, by comparing all of results, it can be said that the TQWT soft thresholding method has been much more successful compared to the hard thresholding and Autoregressive-FX in denoising and converting the actual signal.