1 Introduction

The systematic improvement in geophysical data acquisition over the year’s demand for more innovative data processing methods. Traditional survey designs employ equidistant measurement on a regular grid (Dobróka et al. 2018). Unfortunately, measurements are sometimes taken out of the grid due to several obstacles encountered in the field of survey. Inaccessible sample locations are caused by natural (such as caves) or man-made (buildings) reasons which distorts already planned survey designs. Geophysical field surveys are categorized into regional or local surveys, depending on the size of the area to be covered. Regional surveys are grid-based since they cover a larger area while the local survey is target based involving ground crew. Several factors including the target type and its detectability, the geology of the area, and the equipment available are taken into consideration in the designing process. An acceptable sampling procedure should ensure the acquisition of those data that best resolve specific subsurface features or parameters of interest (Maurer and Boerner 1998; Curtis and Maurer 2000). This is important because no amount of subsequent data processing or analysis can ever compensate for inadequate or missing data that could have contributed significantly to resolving geological targets (Maurer et al. 2010). Appropriate survey design is therefore critical to justify an experiment in terms of the robustness, accuracy, and precision of recovered geological information. Although regular sampling along a survey line has always been the norm for geophysical field measurements, uncertainty sampling (sampling out of sample location) creates some level of statistical errors in processed data. Hence, the development of robust procedures with adequate tools for equidistant and non-equidistant sampling is of utmost importance.

Also, the development of advanced survey equipment which incorporates a global positioning system (GPS) enable easy navigation and even in random-walk data acquisition in recent times. This has necessitated the development of methods for the effective processing of equal in random-walk geophysical measurements. Data acquisition can be simplified and fastened by sampling at even in random intervals along a survey line but not necessarily at regular intervals. The possibility of sampling at irregular intervals with successful processing tools will eliminate the seeming distortions sampling out of grid creates in the data processing. To achieve this, a more robust inversion algorithm with higher noise reduction capabilities is required. Dobróka et al. (2015) presented an inversion based 1D Fourier transformation method known as the Iteratively Reweighted Least Squares Fourier Transform (IRLS-FT) which proved to be an effective tool for noise reduction. It was shown that the noise sensitivity of the continuous Fourier transform (and its discrete variants DFT and FFT) was sufficiently reduced by using robust inversion. The 1D Fourier transform was handled as a robust inverse problem using the IRLS algorithm with Cauchy–Steiner weights. The Fourier spectrum was further discretized using series expansion as a discretization tool. Series expansion based inversion methods were successfully used in the interpretation of borehole geophysical data (Szabó 2011, 2015) and also in processing Induced polarization data (Turai 2011). The method was generalized to 2D, and an application presented in solving reduction to the pole of a magnetic data set (Dobróka et al. 2017). In this paper, it is shown that the newly developed inversion-based Fourier transformation algorithm can also be used in processing non-equidistant (even in random walk) measurement geometry dataset.

2 Theory and methods

The robust Fourier Transformation algorithm investigated in this paper was developed in the framework of the series expansion based inversion methodology of the Department of Geophysics (University of Miskolc). In order to ensure an easily followable discussion, we give a short overview of the method.

2.1 Theoretical background of the IRLS-FT method

The method uses series expansion base discretization of the Fourier spectrum using Hermite functions as basis functions. The entire process is robustified using the IRLS method by the application of Steiner weight instead of Cauchy, thereby enabling an internal iterative recalculation of the weights. Data conversion from the time domain to frequency domain is a common practice in geophysical data processing which improves interpretation, especially in signal processing. This transfiguration can be realized through the application of Fourier transformation. For discretely sampled time domain data sets, the Discrete Fourier Transformation (DFT) algorithm is often applied to determine its Discrete Frequency spectrum (Dobroka et al. 2018). It is well established in inverse problem theory that the simple least squares give optimal solution only when data noise follows Gaussian distribution. For sparsely distributed large errors such as outliers, the estimated model parameters may be highly inconclusive which constitute a limiting factor to the application of the least squares method since geophysical measurements normally contain outliers. To achieve statistical robustness, various methods have been developed over the years to deal with data outliers. A frequently used method in robust optimization is the Least Absolute Deviation (LAD) method which minimizes the L1 norm characterizing the misfit function between the observed and predicted data. Practical experience shows that inversion with minimization of the L1 norm gives more reliable estimates when a smaller number of large errors contaminate the data (Dobroka et al. 2017). Another viable solution involves the use of the Cauchy criterion, which assumes a Cauchy-distributed data noise. The use of data weights in the inversion is very important to ensure each data contribute to the solution based on its error margin. Cauchy inversion is also frequently used in geophysical inversion as a robust optimization method (Amundsen 1991). The integration of the IRLS algorithm with Cauchy weights is a very useful procedure but problematic since the scale parameter of the weights has to be prior known. This challenge has been adequately solved by Steiner (1988, 1997) who derived the scale parameters from the real statistics of the data set in the framework of the Most Frequent Value method (MFV). Szegedi and Dobróka (2014) proposed the use of Steiner weights in inversion based Fourier Transform method whilst (Szűcs et al. 2006) applied the method in groundwater modeling.

2.2 The 1D IRLS-FT algorithm

A change in data from the time domain to the frequency domain can be established using a Fourier transform. The connection enhances data interpretation since certain features are improved in one data format than the other. For the one-dimensional case, the Fourier transform is defined as

$$U(\omega ) = \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - \infty }^{\infty } {u(t)\,e^{ - j\omega t} dt} ,$$
(1)

where \(t\) denotes the time, \(\omega\) is the angular frequency and \(j\) is the imaginary unit. The frequency spectrum \(U(\omega )\) is the Fourier transform of a real-valued time function \(u(t)\) and it is generally a complex valued continuous function. Thus, the Fourier transform provides the frequency domain representation of a phenomenon investigated by the measurement of some quantity in the time domain. The inverse Fourier transform ensures a return from the frequency domain to the time domain.

$$u(t) = \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - \infty }^{\infty } {U(\omega )\,e^{j\omega t} d\omega }$$
(2)

In defining the Fourier transform as an inverse problem, the frequency spectrum \(U(\omega )\) should be described by a discrete parametric model. In order to satisfy this requirement, let us assume that \(U(\omega )\) is approximated with sufficient accuracy by using a finite series expansion

$$U(\omega ) = \sum\limits_{i = 1}^{M} {B_{i} \varPsi_{i} (\omega )} ,$$
(3)

where the parameter \(B_{i}\) is a complex-valued expansion coefficient and \(\varPsi_{i}\) is a member of an accordingly chosen set of real-valued basis functions. Using the terminology of (discrete) inverse problem theory, the theoretical values of time domain data (forward problem) can be given by the inverse Fourier transform

$$u^{theor} (t_{k} ) = u_{k}^{theor} = \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - \infty }^{\infty } {U(\omega )e^{{j\omega \,t_{k} }} d\omega } ,$$
(4)

where \(t_{k}\) is the k-th sampling time. Inserting the expression given in Eq. (1) one finds that

$$u_{k}^{theor} \cong \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - \infty }^{\infty } {\left( {\sum\limits_{i = 1}^{M} {B_{i} \varPsi_{i} (\omega )} } \right)} e^{{j\omega t_{k} }} d\omega = \sum\limits_{i = 1}^{M} {B_{i} \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - \infty }^{\infty } {\varPsi_{i} (\omega )} } e^{{j\omega t_{k} }} d\omega .$$
(5)

Let us introduce the notation

$$G_{ki} = \frac{1}{{\sqrt {2\pi } }}\int\limits_{ - \infty }^{\infty } {\varPsi_{i} (\omega )e^{{j\omega t_{k} }} d\omega } ,$$
(6)

where \(G_{ki}\) is an element of the Jacobian matrix of the size N-by-M. The Jacobian matrix is the inverse Fourier transform of the \(\varPsi_{i}\) basis function. Parameterization of the model is achieved by introducing a special feature of the Hermite functions, by making them the eigenfunctions of the Fourier transform as

$$\fancyscript{F}\left\{{\left. {H_{n}^{(0)} (t)} \right\} =} \right.(- j)^{n} H_{n}^{(0)} (\omega),$$
(7)

and respectively for the inverse Fourier transform

$$\fancyscript{F}^{{-{1}}} \left\{{\left. {H_{n}^{(0)} (\omega)} \right\} =} \right.(j)^{n} H_{n}^{(0)} (t).$$
(8)

They are modified by scaling because, in geophysical applications, the frequency covers wide ranges. The theoretical values can, therefore, be written in the linear form

$$u_{k}^{theor} = \sum\limits_{i = 1}^{M} {B_{i} G_{k,i} } .$$
(9)

2.3 The 2D IRLS-FTalgorithm

The 2D Fourier transform of a function u(x,y) can be calculated by the integral

$$U(\omega_{x} ,\omega_{y} ) = \frac{1}{2\pi }\int\limits_{ - \infty }^{\infty } {\int\limits_{ - \infty }^{\infty } {u\,(x,y)e^{{ - j\,(\omega_{x} x + \omega_{y} y)}} dxdy} } ,$$
(10)

its inverse is given by the formula

$$u(x,y) = \frac{1}{2\pi }\int\limits_{ - \infty }^{\infty } {\int\limits_{ - \infty }^{\infty } {U(\omega_{x} } ,\,\omega_{y} )} e^{{j\,(\omega_{x} x + \omega_{y} y)}} d\omega_{x} d\omega_{y} ,$$
(11)

where x, y are the spatial coordinates, U(ωxy) is the 2D spatial-frequency spectrum and ωx, ωy indicate the spatial-angular frequencies. The discretization of the continuous spectrum can be done through series expansion,

$$U(\omega_{x} ,\omega_{y} ) = \sum\limits_{n = 1}^{N} {\sum\limits_{m = 1}^{M} {B_{n,m} \varPsi_{n,m} (\omega_{x} ,\omega_{y} )} } ,$$
(12)

where Ψn,mxy) are frequency dependent basis functions, Bn,m are the expansion coefficients which represent the model parameters of the inverse problem. The basis function system should be square integrable in the interval (− ∞, ∞). The Hermite functions meet this criterion with an additional advantage as discussed in the 1D. Dobróka et al. (2015) showed that the elements of the Jacobian matrix can be considered as the inverse Fourier transform of the basis function system. Therefore, they can be calculated more easily if the basis functions are chosen from the eigenfunctions of the inverse Fourier transformation. It can be shown, that the normed and scaled Hermite functions given by

$$H_{n} (\omega_{x} ,\alpha ) = \frac{{e^{{ - \frac{{\alpha \,\omega_{x}^{2} }}{2}}} h_{n} (\omega_{x} ,\alpha )}}{{\sqrt {\sqrt {\frac{\pi }{\alpha }} \,n\,!\,(2\alpha )^{n} } }}\,,$$
(13)
$$H_{m} (\omega_{y} ,\beta ) = \frac{{e^{{ - \frac{{\beta \,\omega_{y}^{2} }}{2}}} h_{m} (\omega_{y} ,\beta )}}{{\sqrt {\sqrt {\frac{\pi }{\beta }} \,m\,!\,(2\beta )^{m} } }}\,,$$
(14)

are eigenfunctions of the inverse Fourier transformation and the Jacobian matrix of the inverse problem can be written as

$$G_{k,l}^{n,m} = \,\frac{{(j)^{n + m} }}{{\sqrt[4]{\alpha \beta }}}H_{n}^{(0)} \left( {\frac{{x_{k} }}{\sqrt \alpha }} \right)\,H_{m}^{(0)} \left( {\frac{{y_{l} }}{\sqrt \beta }} \right)\,\,.$$
(15)

here \(H_{n}^{(0)} ,H_{m}^{(0)}\) denote the non-scaled Hermite functions and provides a fast solution to the forward problem

$$u(x_{k} ,y_{l} ) = \sum\limits_{n = 1}^{N} {\,\sum\limits_{m = 1}^{M} {B_{n,\,m} } } G_{k,l}^{n,m} .$$
(16)

2.4 The inversion algorithm

The Gaussian least squares method (LSQ) which minimizes the \(L_{2}\)-norm of the deviation vector between the observed and calculated data is normally applied when the data noise follows the regular distribution. Unfortunately, most geophysical data contains irregular noise with randomly occurring outliers making the least squares method (LSQ) less effective for processing. Dobroka et al. (2012) emphasized the possibilities of obtaining a good result in inverse problem solution when the data is weighted. To develop a robust algorithm, we minimized the weighted norm of the deviation vector using Cauchy weights which were further modified to Cauchy–Steiner weights. The minimized weighted norm is given as

$$E_{w} = \sum\limits_{k = 1}^{N} {w_{k} e_{k}^{2} }$$
(17)

where \(e_{k}^{{}}\) is the difference between the measured and predicted data, \(w_{k}\) is the Cauchy weight, given by

$$w_{k} = \frac{{\sigma^{2} }}{{\sigma^{2} + e_{k}^{2} }}$$
(18)

for the k-th datum (\(\sigma\) being the scale parameter). Inversion procedures involving Cauchy weights are problematic since the scale parameter has to be determined a priori but easily solvable by using Steiner-weights (Steiner 1988). In the context of Steiner’s most frequent value method (MFV), the scale parameter \(\sigma^{2}\) can be determined from the data set in an internal iteration loop. The stop criterion is defined by experience from a fixed number of iterations. After this, the Cauchy weights are calculated by using the Steiner’s scale parameter. The so-called Cauchy–Steiner weights at the last step of the internal iterations are given by

$$w_{k} = \frac{{\varepsilon^{2} }}{{\varepsilon^{2} + e_{k}^{2} }},$$
(19)

where \(\varepsilon\) is the Steiner’s scale factor called dihesion which is determined iteratively.

In practice, the misfit function is non-quadratic in the case of Cauchy–Steiner weights (because \(e_{k}\) contains the unknown expansion coefficients) and so the inverse problem is nonlinear which can be solved again by applying the method of the iteratively reweighted least squares (Scales et al. 1988). In the framework of this algorithm, a 0-th order solution \(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {B}^{(0)}\) is derived by using the non-weighted LSQ method and the weights are calculated as

$$w_{k}^{(0)} = \frac{{\varepsilon^{2} }}{{\varepsilon^{2} + (e_{k}^{(0)} )^{2} }}$$
(20)

with \(e_{k}^{(0)} = u_{k}^{measured} - u_{k}^{(0)}\), where \(u_{k}^{(0)} = \sum\limits_{i = 1}^{M} {B_{i}^{(0)} G_{ki} }\) and the expansion coefficients are given by the LSQ method. In the first iteration the misfit function

$$E_{w}^{(1)} = \sum\limits_{k = 1}^{N} {w_{k}^{(0)} e_{k}^{(1)2} }$$
(21)

is minimized resulting in the linear set of normal equations

$${\mathbf{G}}^{T} {\mathbf{W}}^{(0)} {\mathbf{G}}\vec{B}^{(1)} = {\mathbf{G}}^{T} {\mathbf{W}}^{(0)} \vec{u}^{measured}$$
(22)

The minimization of the new misfit function

$$E_{w}^{(2)} = \sum\limits_{k = 1}^{N} {w_{k}^{(1)} e_{k}^{(2)2} }$$
(23)

gives \(\vec{B}^{(2)}\) which serves again for the calculation of \(w_{k}^{(2)} .\) This procedure is repeated giving the typical j-th iteration step

$${\mathbf{G}}^{T} {\mathbf{W}}^{(j - 1)} {\mathbf{G}}\,\vec{B}^{(j)} = {\mathbf{G}}^{T} {\mathbf{W}}^{(j - 1)} \vec{u}^{measured}$$
(24)

with the \({\mathbf{W}}^{(j - 1)}\) weighting matrix

$$W_{kk}^{(j - 1)} = w_{k}^{(j - 1)}$$
(25)

Each step of these iterations contains an internal loop for the determination of the Steiner’s scale parameter which is repeated until a proper stop criterion is met.

3 The applicability of IRLS-FT method in non-equidistant sampling

In applying surface geophysical methods the measurement points are often positioned along a regular grid. This requires expensive geodetic work. There is an old purpose to reduce this expense by developing data processing methods applicable also in case of non-regular (even randomly collected) data sets (Sauerländer et al. 1998). In the following, we investigate the applicability of the IRLS-FT method in processing non-equidistantly sampled data sets.

3.1 The 1D application of the IRLS-FT method

A Morlet waveform (Fig. 1) was created to test the applicability of the new method in one dimension. The noiseless time function of the test data can be described by the formula below

$$u(t\,) = \kappa {\kern 1pt} t^{\eta } e^{ - \lambda t} \sin (\omega t + \varphi ),$$
(26)

where the Greek letters represent the parameters of the signal. We specified fixed values for the signal parameters as follows: \(\eta = 2\), \(\lambda = 2 0\), \(\omega = 40\pi\), \(\varphi = \pi /4\). The time domain signal can be seen in Fig. 1 (\(\kappa\) is determined so that \(u_{\hbox{max} } = 1).\)

Fig. 1
figure 1

The calculated Morlet waveform

As a first step, the Morlet waveform was sampled equidistantly in 401 points ranging over the time interval of [− 1, 1] and processed using the traditional DFT method to give both the real and imaginary parts of Fourier Transform. The same equidistantly sampled waveform was also processed using the IRLS-FT method. The result is shown in Fig. 2. Both the traditional DFT and the IRLS-FT gave similar real and imaginary parts of the spectrum (all the applied softwares were written in Matlab).

Fig. 2
figure 2

The spectrum of the equidistantly sampled Morlet waveform processed with DFT and also with IRLS-FT method

To test the efficiency of the IRLS-FT method in random-walk measurements, the Morlet waveform was sampled randomly for processing. In this experiment, the same number of samples were used with randomly selected positions in the whole time interval. The randomness of the sampling is demonstrated in Fig. 3. Figure 4 shows the IRLS-FT spectrum of the signal. It can be seen, that both its real and imaginary parts are exactly the same as that, found in regularly sampled case (Fig. 2). This proves that the inversion based Fourier Transform method gives the same results in processing both regularly and non-regularly sampled data set.

Fig. 3
figure 3

The randomly selected sample points used for sampling the Morlet waveform

Fig. 4
figure 4

The IRLS-FT spectrum of the randomly sampled Morlet waveform

3.2 The 2D application of the IRLS-FT method

To prove the applicability of the 2D inversion based IRLS FT method, it was tested on synthetic magnetic data sets sampled at regular equidistant and non-equidistant intervals. In all, 41 × 41 “measurement points” were sampled along a 5 m × 5 m regular grid and further randomized to obtain non-equidistant measurements. Data were generated for a surface between ± 100 m both in the x and y directions above a ‘CL’ shaped magnetic body (inclination I = 63°, declination D = 3°, magnetization 200 nT). The surface magnetic data were calculated by the Kunaratnam (1981) method and was subsequently reduced to the pole (I = 90°) by applying the formula in the frequency domain

$$R(u,v) = T(u,v)S(u,v),$$
(27)

where T(u,v) is the 2D Fourier transform of the magnetic data set, S(u,v) is the frequency domain operator of pole reduction and R(u,v) is the reduced data set after the data reduction process. First, the reduction to the pole was performed by using the conventional DFT algorithm on equidistantly sampled magnetic data. The map of noiseless magnetic data on the equidistant grid and its version reduced to pole are given in Fig. 5 (left) and (right). The 2D spectrum of the regularly sampled data set was calculated by means of DFT (Fig. 6 (left)) and IRLS-FT (Fig. 6 (right)) (both 2D softwares were written in Matlab).

Fig. 5
figure 5

Noise-free magnetic map calculated on equidistant grid (left), noise-free equidistant magnetic map reduced to the pole using the conventional DFT (right)

Fig. 6
figure 6

DFT spectrum of the magnetic data set calculated on equidistant grid (left), IRLS-FT spectrum of the magnetic data set calculated on equidistant grid (right)

(For numeric reasons the calculations were made on the data set transformed to [− 1,1] in both x and y coordinates resulting in an appropriate scale in the wavenumber domain.) The IRLS-FT spectrum was calculated using Hermite functions of the (maximal) order of M = 28. The similarity of the two results are obvious and can be increased by using higher order basis functions. This has a consequence in increased computation time and in a rapid change of the condition number. The M = 28 was found as a good compromise between accuracy and stability.

The total magnetic intensity data were calculated on a randomized set of sample points. The randomly selected x and y coordinates are shown in Fig. 7, the sample points are defined by ordering all y coordinates to all of the x position resulting in 41 × 41 sample points. The magnetic field data calculated in these non-regular positions were processed using the 2D IRLS-FT algorithm, the resulting spectrum is given in Fig. 8 (for the sake of comparability the IRLS-FT spectrum calculated in the regular grid is also shown).

Fig. 7
figure 7

The randomly selected x and y coordinates used for generating 41 × 41 sample points

As it was observed in the above 1D example, the inversion based 1D IRLS-FT algorithm can successfully be applied in processing non-regularly sampled data set. This important observation is extended to the two-dimensional case by testing the 2D IRLS-FT method on the randomly selected set of sample points. The result is shown in Fig. 8, where for the sake of comparison the 2D IRLS spectrum of the regularly sampled data sets (same as in Fig. 5 on the right) is also presented.

Fig. 8
figure 8

IRLS-FT spectrum of the magnetic data set calculated on equidistant grid (left), IRLS-FT spectrum of the magnetic data set calculated on the non-equidistant grid (right)

As can be seen, the two spectra are the same (both calculations were performed with the same inversion parameters using series expansion by means of Hermite functions of M = 28 order). The produced spectra prove that the 2D IRLS-FT algorithm has also a very important feature that it works on both regularly or irregularly sampled data sets. From this fact it is straightforward to expect, that the method gives the same result in reducing to pole for both regularly sampled (in Fig. 9, left) and randomly sampled data sets (Fig. 9, right).

Fig. 9
figure 9

The pole reduced magnetic data sets using IRLS-FT. Equidistant sampling was used in the left, non-equidistant (random) sampling was used in the right

As it was observed in the above 1D example, the inversion based 1D IRLS-FT algorithm was successfully applied in processing non-regularly sampled data set. This important observation has been extended to the two-dimensional case by testing the 2D IRLS-FT method on the randomly selected set of sample points to give a similar result.

4 Conclusion

We present a further theoretical development of the robust inversion based IRLS-FT method for application in random walk geophysical measurement data processing. This new applicability enables field measurements to be taken at unequal intervals, in that, measurements are obtained without preliminary geodetic work. The introduced IRLS-FT method treats the Fourier transform as an inverse problem. The complex spectrum is discretized by series expansion and the inversion problem is solved for the series expansion coefficients. Taking advantage of the fact, that the Hermite functions are eigenfunctions of the Fourier Transform, they were chosen as basis functions making the algorithm quicker for computation of the Jacobian matrix even in 2D problems.

Both one and two-dimensional applications of the IRLS- FT showed favorable results for the new method. In one dimension, IRLS-FT processed waveforms were similar for both equidistant and non-equidistant (even random) sampling. An application on magnetic data showed similar results indicating that the new FT method is applicable irrespective of the sampling protocol applied in the field survey or data acquisition process. The data processing strength of the IRLS- FT method quickens field data acquisition as measurements are not necessarily taken on a grid.

Comparatively, the inversion-based Fourier transformation algorithm proved to have a higher noise rejection capability than the traditional DFT method as demonstrated in reduction to pole of magnetic data (Dobróka et al. 2017). In this paper, it was further shown that the inversion-based Fourier transformation algorithm can be effectively used in processing data set collected in non-equidistant (even in random walk) measurement geometry.