SPECTRAL METHOD FOR NUMERICAL CALCULATION OF DERIVATIVES IN DIGITAL PROCESSING OF SUBSURFACE RADAR SOUNDING SIGNALS

The practical realization of the morphologic analysis of ra dar subsurface as the processing of input digital data is presented here together with real data processing results. The analysis depends on a new numerical differentiation alg orithm in frequency domain. It is compared with a central finite differences scheme and resu lts of computations demonstrate approximately 80-times decrease of the maximum errors.


Introduction
The task of high accuracy numerical differentiation of discrete data is important not only in such known areas of the applied sciences as solution of partial differential equations [4], but also in solution of some specific tasks, namely, for the analysis of signals during subsurface radio location of ground surface segments. The main goal of this processing is to define from radar sounding signals the information about targets, e.g. ground inner objects or a localization of geological structures layers and its properties. In this paper we propose the so-called morphologic approach [3].
The main idea of morphologic analysis of subsurface radar signals [1] is to use the methods of calculation geometry instead of formal equations describing signal deformation. The axiom of space homogeneity and isotropy asserts that solutions of equation defining the propagation of disturbances in environment have the same form for all coordinate systems that are obtained by translation and rotation of axes. On the other hand, according to the Curie principle, the electrical field in environment has a characteristic symmetry, which is a subgroup of cylindrical transformation group [1].
Violation of the homogeneity and isotropy leads to a non-symmetry, i.e. to typical deformation of signal form. When a signal is presented by a smooth curve, then the problem is reduced to the classification of characteristic points on the curve, the number of such points is always finite. By a smooth curve we mean that at each point there exists a derivative and, therefore, it is always possible to construct a tangent and normal vectors, i.e. to introduce a local coordinate system.
There are three different cases of mutual location of points relatively to the local coordinate system: in Fig.1 (I), (II) and (III) correspond to the linear, quadric and cubic forms, respectively. The other types of points cannot exist for a signal set in the finite interval. A relative distribution of points of each type is changing depending on a sounding environment, and points of cubical type are connected with reflections from radar targets. Indeed, according to the Newton principle, the deformation of the curve is proportional to the second derivative, and the cubic form satisfies this condition (see Fig.1 (III)).
Let us consider reflection of the separate realization of sounding pulse S(t) (Fig.2) onto the space curve given in the form of complex analytical signal where the radar signal is a real part u re (t) = S(t) = Re − → U (t)], and its Hilbert transformation u im (t) = H{S(t)} = Im − → U (t) is an imaginary part of the complex signal. Here U 0 (t) is the envelope, and Φ(t) is the instantaneous phase on the interval (0 < t < T ) of one development of duration T , and H{. . .} is the fast Hilbert transformation [5]. In [3] it is defined that the analysed process is described by three parameters: two external parameters, i.e. the radial υ R and tangential υ T components of the instantaneous velocity  is connected to the fact that a reflection group always generates convolution-type (A 3 -type) catastrophe [1] (Fig.3). Then, parameters of the form a(t), b(t) describe a function with twice degenerated special points [2], this allows us to localize the position of targets during subsurface sounding The points belonging to the domain outside of the coloured region defined by curve (1.2) indicate the target. It is the basis for understanding the morphologic method.

Central Finite Differences Scheme and Spectral Method
As we see from (1.1), the quality of the analysis depends directly on the accuracy of numerical approximation of derivatives. Let us consider this question in detail. The Central Finite Difference (CFD) approximation is defined for a discrete data: where △t is a time step, and t = k △ t, (k = 1, 2, . . . , N t − 1) is the discrete time. Such approximation has the second order of accuracy: Spectral Method for Numerical Calculation of Derivatives in Radar Signals 35 The CFD scheme has its analog in the frequency domain. Let u(t) be a discrete signal, its spectrum is obtained by the Fast Fourier Transformation (FFT): where ω = ν △ f, ν = 0, 1, . . . , N f − 1 is a discrete variable in frequency domain with step △f and F {. . .}is the FFT transformation. Every sample of the signal can be considered as a shift by one step forward or backward. According to the Fourier shift theorem [5] it is equivalent to the multiplication of the spectrum by a complex exponent: thus the analog of CFD formula in the frequency domain can be written as By using the Euler formula for the complex exponent we obtain: where j = √ −1 is the imaginary unit. Thus in frequency domain the CFD scheme means a multiplication of the spectrum by imaginary Sine-function (similar to digital filtering methods [5]): By repeating the same rule after triple multiplication we get the approximation of the third order derivative: Using these formulas for a digital signal we can eliminate the second order error term in (2.2) and reduce high order terms by using a compensation in the frequency domain: 3) The implementation of the method consists of two parallel processes, which use two separate inputs of the signal for FFT in order to reduce a computation time. One of them is filtered to approximate the first derivative, and the second approximates the third derivative. Then, before application of the inverse FFT, they are added to each other.

Numerical Experiment
In order to compare the accuracy of the proposed method and CFD scheme we applied both methods for numerical computation of derivatives of the artificial signal, generated by the modulated Gauss pulse (see Fig.5), which is similar to radar main pulse given in Fig.4: Discrete data is given at t = 1, . . . , 512 and N f = 1024. The first derivative of S is defined as In order to reduce possible boundary effects we process the signal on a double time interval with its mirror imaging. Fig.6 presents results obtained by the spectral method. The errors of numerical approximations of the first derivative by the CFD scheme and the spectral method are presented in Fig. 7.  The results show essential difference between CFD and spectral methods. Comparing both methods we see a substantial difference in location of the peak value of the test signal. The error was computed by using the analytical solution for derivative. We conclude that the maximum error of spectral approximation is equal to 3.7 × 10 −6 and it is 80-times smaller than the error for CFD method (which is equal to 2.7 × 10 −4 ). This result depends on a shape of the signal and the accuracy of the spectral approximation can be improved by further compensation of high order terms in series (2.2) in the same manner as it was demonstrated above.

Radar Sounding Data Processing
We implemented the proposed method to process data sets obtained by subsurface radar sounding of various ground places. Each separate signal from data set given in Fig.4 is shown in Figures 8-11 as a vertical line, coloured according to its value. The current number of each separate signal (or trace) depends on the spatial position of the radar antenna. Thus, all separate signals produce a two-dimensional picture or image called a radar "profile". In Fig.8 we present the profile of the test polygon with three rows of small-sized targets, an individual target is considered in Fig.9, a system of underground pipelines is shown in Fig.10, and, finally, the geological section of some ground surface strip is considered in Fig.11.
The morphologic method of signal processing leads to elimination of cross-reverberation and interference of reflection images from various targets and geological layers. This property of the spectral method and computation of derivatives (1.1) with high accuracy gives a possibility to increase a visual clarity of radar sounding profiles and location of subsurface targets or ground geological layers. Figure 10. Subsurface radar sounding profile: before (left) and after (right) processing. Underground pipelines system. Figure 11. Subsurface radar sounding profile: before (left) and after (right) processing. Geological section of some ground strip.

Conclusion
The spectral method requires much larger computation time in comparison with finite central differences method because application of direct and inverse fast Furrier transformation procedures are time consuming. Thus this method is not useful for real-time processing applications. On the other hand, we can eliminate errors produced by finite differences schemes. If high accuracy of derivatives is necessary, one can implement the algorithm in such a way that most of problems, connected with processing quality, stability of the method and convergence, may be successfully solved.