1 Introduction

Body wave velocity measurements (P and S) are typically carried out via the seismic downhole (SDH) method [1]. The precision of the determination of the velocities of body waves affects the application of downhole data to other research related to, for example, liquefaction [2-4], elastic modulus determination [5,6], earthquake engineering (i.e., seismic zonation maps) [7,8], land subsidence

Received July 26 th, 2024, Revised September 13th, 2024, Accepted for publication October 2 nd, 2024 Copyright © 2024 Published by ITB Institut for Research and Community Service, ISSN: 2337-5760, DOI: 10.5614/j.math.fund.sci.2024.56.2.3

assessment [9], mining exploration [10], and data integration with other geotechnical engineering methods [11].

The measurement reads the travel times of body waves – primary and shear waves – from the seismic energy source to the geophone and continues to record the velocity of the P or S waves. The source is placed at some horizontal distance from the borehole and the depth of the geophone is changed every meter to collect measurement data until the required total depth (Figure 1). The real distance between the source and the receiver depends on the offset distance between the source and the axes of the borehole; then, the velocities are calculated by using the real distance and not by measurement depth.

Figure 1 Seismic downhole (SDH) measurement.

The ideal measured response is the arrival of the P waves followed by the arrival of the S waves propagating along the shortest path between the source and the geophone sensor. However, waves may also travel indirectly along the borehole via the borehole case and its filling material, which could be faster than the velocity of the soil layer, such as tube waves [12], or interact with their reflections [13], creating dynamically changed waveforms and attenuating as the depth of measurement increases [14]. This decreases the accuracy of seismic body wave picking [15].

The expected S wave can be separated from other waves using an F-K filter [16]. The F-K image is commonly obtained by performing the fast Fourier transform (FFT) on the time and distance axes simultaneously to get an image response for frequency versus the wave number = −1 dimension, where is the wavelength [17]. The slope in the F-K plane indicates the wave velocity where = . = . −1 . F-K filtering is done by zeroing the response value in the undesired velocity range. This also means that correct filtering results are only possible if the F-K response of the filtered wave group is separated from the K response of another wave group. However, the ratio between P and S wave

velocity \((v_p/v_s)\) for unsaturated rock/soil is always less than 4 [18]; this is further complicated by high P and S velocities, which widen the wavelength compared to the distance between measurements.

Standard F-K transform needs the data to be equally spaced, both in time and distance, and produces a biased response in SDH data transformation since the distance for each depth of measurement is not equally spaced (Figure 1). This can be corrected by implementing the non-uniform F-K transform (nuFK) in the distance dimension. In addition to its complexity, the source energy generation is not equal for each measurement, as the FFT process will produce wider spectral responses from non-sinusoidal signals or discontinuity [19]. The F-K response can also experience various responses in the frequency axis and wave number as shadow response and direct F-K filtering generally does not accommodate these wider responses to reconstruct the correct waveform. A common way to reduce this discontinuity is using normalization [20], i.e., reducing each trace to its mean value and dividing the results by the maximum value or standard deviation of each trace. This standard normalization method will be acceptable only if the highest amplitude of each trace is from the same wave tracing. However, the energy of the traveling wave decays over distance and the traveling wave may reflect or impose on another wave or high amplitude noise, which can lead to incorrect normalization, thus creating a shadow in its F-K transform response. We propose a two-step normalization method to reduce this shadow response.

2 Methods

2.1 Non-uniform F-K (nuFK) Transform

The F-K transform is performed by implementing the discrete Fourier transform (DFT) to each dimension in time (t)-distance (x) data to produce the frequency (f)-wave number (k) response. The standard implementation of DFT requires the data to have equal spacing, resulting in an equally spaced transform response. In time series data with length N, the sampling rate defines the spacing R and the transform presents the frequency response with the maximum range of half its sampling rate and frequency spacing equal to the sampling rate divided by N. The DFT is the discrete implementation of the Fourier transform as expressed in Eq. (1) [21]:

\[Y(\omega) = \int_{-\infty}^{\infty} y(t)e^{-j\omega t}dt;\] \[Y(f) = \sum_{n=0}^{N-1} Y_{fn} = \sum_{n=0}^{N-1} e^{-j2\pi \left(\frac{f}{N}\right)t_n} y_n\] (1)

where \(t_n = \frac{n-1}{R}\), \(f \in f_n = \frac{(n-1)R}{N}\), \(n = \{1, 2, ..., N\}\), and the full-size of Y is equal to N. A complete transformation Y is achieved by calculating the value \(Y_{fn}\) for each \(f_n\) in Eq. (2). Vice versa, the inverse transform will be completed by calculating each response \(y_n\) for all \(t_n\) responses.

\[Y = \begin{bmatrix} Y_{f_{1}} \\ \dots \\ Y_{f_{N}} \end{bmatrix} = \begin{bmatrix} e^{-j(2\pi f_{1}t_{1})} & \dots & e^{-j(2\pi f_{1}t_{N})} \\ \vdots & \ddots & \vdots \\ e^{-j(2\pi f_{N}t_{1})} & \dots & e^{-j(2\pi f_{N}t_{N})} \end{bmatrix} \begin{bmatrix} y_{1} \\ \dots \\ y_{N} \end{bmatrix}\] \[Y = \begin{bmatrix} e^{-j\left(2\pi \left(\frac{R.0}{N}\right)\left(\frac{0}{R}\right)\right)} & \dots & e^{-j\left(2\pi \left(\frac{R.0}{N}\right)\left(\frac{N-1}{R}\right)\right)} \\ \vdots & \ddots & \vdots \\ e^{-j\left(2\pi \left(\frac{R.(N-1)}{N}\right)\left(\frac{0}{R}\right)\right)} & \dots & e^{-j\left(2\pi \left(\frac{R.(N-1)}{N}\right)\left(\frac{(N-1)}{R}\right)\right)} \end{bmatrix} \begin{bmatrix} y_{1} \\ \dots \\ y_{N} \end{bmatrix}\] \[(2)\]

Manipulation operations such as frequency filtering can be done by leaving some specific frequency range \(f_x\) uncalculated. \(Y_{f_x} = 0\), i.e., a low-pass filter with frequency limit \(f_A\), will have \(Y_{f_x} = 0\) in the range \(f_A < f_x < f_N\).

DFT of non-equally spaced data can be performed to create an equally spaced frequency response of the data to a new R in Eq. (3), also known as adjoint non-uniform discrete Fourier transform (NDFT) [22], which nowadays is used to optimize the non-uniform fast Fourier transform (NFFT/NFFTN) by using a faster-adjusted algorithm:

\[Y(f) = \sum_{n=0}^{N-1} Y_{n,f} = \sum_{n=0}^{N-1} e^{-j2\pi \left(\frac{f}{N}\right)t_n} y_n w_n; \ w_n = \Delta t_n R\] \[Y = \begin{bmatrix} Y_{f_{1'}} \\ \cdots \\ Y_{f_{N'}} \end{bmatrix} = \begin{bmatrix} e^{-j(2\pi f_{1'}t_1)} & \cdots & e^{-j(2\pi f_{1'}t_N)} \\ \vdots & \ddots & \vdots \\ e^{-j(2\pi f_{N'}t_1)} & \cdots & e^{-j(2\pi f_{N'}t_N)} \end{bmatrix} \begin{bmatrix} y_1 w_1 \\ \cdots \\ y_N w_N \end{bmatrix}\] (3)

where \(f \in f_{n'} = \frac{(n'-1)R}{N'}\), \(n' = \{1,2,...,N'\}\), \(t_n\) is an unequally spaced time series and the complete size of Y is equal to N'. Adjoint NDFT/NFFT will have transformed and inverse new single sample rate properties. Tighter data spacing will produce a better response and for a wider data spacing, of course, this will create an imperfect transformation, since the value for each frequency response is the summation of less information, which in our method is adjusted by using weight contribution \(w_n\), calculated from the neighboring spacing (nuFT). Based on Eq.(3), the new element number N' must be equal to or larger than \(t_N \times R\) to prevent circular stacking between the start time and the end time.

The inverse transform (inuFT) can now be performed to produce an equally spaced response y' with length N' and sample rate R as in Eq. (4), with an example shown in Figure 2. Imperfect transformation may also produce imaginary parts of the inverse transform, which can be neglected since the original data is completely real:

\[y_{t} = \Re\left\{\frac{1}{N'} \sum_{n'=0}^{N'-1} Y_{n',t}\right\} = \frac{1}{N'} \sum_{n'=0}^{N'-1} e^{-j2\pi \left(\frac{t}{N'}\right) f_{n'}} Y_{n'}\] \[\begin{bmatrix} y'_{n,t_{1}} \\ \cdots \\ y'_{n,t_{N'}} \end{bmatrix} = \Re\left\{\frac{1}{N'} \begin{bmatrix} e^{j(2\pi t_{1}f_{1})} & \cdots & e^{j(2\pi t_{1}f_{N'})} \\ \vdots & \ddots & \vdots \\ e^{j(2\pi t_{N'}f_{1})} & \cdots & e^{j(2\pi t_{N'}f_{N'})} \end{bmatrix} \begin{bmatrix} Y_{1} \\ \cdots \\ Y_{N'} \end{bmatrix} \right\}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\downarrow \text{time series}\] \[\text{[rumus tidak dapat ditampilkan dengan baik — lihat PDF asli]}\]

Figure 2 Left: inuFT result from nuFT (blue) and iNFFT (green) of non-uniform spaced data (black dots) compared to the uniform spacing of the reference signal (red dots). Right: uniform frequency response of nuFT (blue) and NFFT (green) from non-uniform time series data compared to the FFT response of the uniform reference signal (red).

In a two-dimensional case with time as the first dimension and distance as the second dimension, the non-uniform DFT (nuFT) is performed on each dimension sequentially to produce an \(M' \times N'\) size response, where M' and N' are the new sizes for each dimension, as expressed in Eq. (5):

\[Y_{m',n'} = \sum_{m'=1}^{M'} \sum_{n'=1}^{N'} \exp\left(\frac{-j2\pi f_{m'} t_m}{M'}\right) \exp\left(\frac{-j2\pi k_{n'} x_n}{N'}\right) y_{m,n} w_{m,n}\] (5)

with \(f_{m'}=(0:M'-1)R_1/M'\) is the equally spaced frequency response with length M' at sample rate \(R_1\); \(k_{n'}=(0:N'-1)R_2/N'\) is the equally spaced wave number response with length N' at sample rate \(R_1\); and \(Y_{m',n'}=Y(f_{m'},k_{n'})\). The \(Y(f_{m'},k_{n'})\) plane describes the velocity response or nuF-K response image calculated from non-uniformly spaced data. An F-K filter can be applied to a specific velocity range as low pass, high pass, band pass, or notch filter. The inverse F-K filter will have uniform spacing both in the time and the space dimension [23]. NuF-K in this paper was programmed using MATLAB in DFT basis and not using the faster algorithm (NFFT, [21,22,24]), since the latter may change the signal shape (Figure 2).

2.2 Two-Step Normalization

The records may have many different wave paths with different velocity and frequency characteristics. In the SDH method, the main wave is the wave group that travels directly from the source point to the receiver. Better normalization can be achieved by stricter normalization for this group.

Figure 3 Overall process of nuF-K transform with two-step normalization.

Assuming the main wave has the same shape and velocity, thus producing a highenergy F-K response, the normalization coefficient is calculated from this main wave's velocity band-passed F-K response. We used a Gaussian tapered window for the bandpass () function as the function of velocity limit and Gaussian shape factor as in formulated Eq. (6). The inverse Fourier transform of filtered F-K in the dimension results in a spectral vs distance response that is used to normalize the vs response before final filtering is applied. The complete procedure is shown in Figure 3.

\[FK_{BP} = \exp(-a(v - v_{lim})^2) \tag{6}\]

3 Results

3.1 Uniform vs Non-Uniform FK

As an example, we used synthetic data from the case for measurements on homogenous and three-layer models (Figure 4). The top surface layer is exposed to air. The depth spacing for the sensor ∆h is 1 m, and the total depth is 30 m.

Figure 4 Synthetic model configuration and acquired data for homogenous and three-layer system. Distant source to borehole/sensor axis (5 m) creating hyperbolic time arrival in one-layer homogenous data. The three-layer system with closer source to sensor axis distance is similar to our real measurement condition in the next sub-chapter.

By performing the nuF-K transform on the hyperbolic distance axis and setting the target spacing increment (x) to 0.76 m, the velocity response in the F-K plane will be correct when compared with the results of the uniform F-K (uFK) transform on the data (Figure 5). V^s in uF-K is slightly higher, at around 278 m (yellow), calculated from the slope between the two maxima of the uniform F-K image.

Figure 5 The nuF-K transform correcting the velocity response of the standard F-K transform. Each line is a constant velocity line indicating the layer velocity of the model. The yellow and magenta line in the left picture shows a slightly higher velocity response than the model calculated from the peak position of the uniform F-K () image.

3.2 Two-Step Normalization vs Per-Trace Normalization

In each source generation, the source waveform can change due to many factors, such as inconsistent power, strike slippage, compaction of layers due to the strike, or destruction of the strike base. To show this random condition, we used the synthetic data above, multiplied by the random amplitude with a maximum random value equal to half of the maximum value of each trace and the addition of low-frequency random noise (Figure 6).

Figure 6 uF-K in time-depth SDH noisy data ( > ) produced an incorrect velocity response with a shadow trail (left). uF-K corrects the velocity response (middle) and applies the second normalization to reduce the shadow response and noise (right).

The presence of noise in the F-K image depends on the velocity path of the noise. Thus, the random noise with no specific path will contribute less to the F-K amplitude. The shadow trail in the same frequency range as the main wave indicates that there is a discontinuity or non-sinusoidal response, where the Fourier transform response will be in the form of a point if the original signal is a pure sinusoid and widens very far in the Dirac delta function transformation [19]. This extreme change in distance response may occur in parts of the recording not only by random parameters but can also be formed by two or more

groups of waves with different speeds, such as in shallow recordings where the responses of P waves, S waves, and surface waves meet each other. Two-stage normalization aims to reduce this effect by adjusting the trace amplitude to follow the amplitude of groups of waves in the same F-K group, including traces that experience noise. When the shadowing effect is diminished or reduced, strict band-pass () filtering can be applied with less missing information left in the shadow (Figure 7).

Figure 7 By performing nuFKnorm2 to correct the trace position and reduce the shadowing effect (left), strict band-pass filtering will produce a clean S wave response (middle) and the P wave response will be left in its residue response as a notch filtered signal (right).

Figure 8 Left: The F-K responses for each velocity in the three-layer model system are too close to each other with a very large energy difference between P waves (red) and S waves (green). Such condition will create a good, filtered response in the higher energy velocity group (middle) while creating a poor response in the lower energy group (right).

The quality of the F-K image depends on the consistency of the waveform, both in time and distance, thus forming energy spectra that are similar in frequency and wave number. In real measurements, this shape consistency is quite challenging to obtain because the generation of the wave source can be different at every measurement depth. Strict band filtering on recordings of wave groups with inconsistent polarities can cause undesirable responses to satisfy the + = condition. Thus, the inconsistent polarity response of such an S wave is taken from the residue of one or several band-pass filters of the wave group with higher consistency in the data.

3.3 Real Measurement Data Filtering

The data comes from SDH measurements that we carried out in Yogyakarta. The surface of Yogyakarta's mountainous region consists of andesitic rock with high seismic velocities, covered by an alluvial layer from volcanic ash with low seismic velocity [25], which causes high-energy surface reflection. Complete measurement was performed from three source types at a 3-meter horizontal distance from the borehole axis: a vertical blow to generate a P wave, a first horizontal blow (S1), and a reverse horizontal blow (S2) to generate a reversal S wave source. We focused on the reversal S wave data, since the P wave always comes earlier, biasing the S wave arrival.

Figure 9 Comparison between uF-K with depth spacing ∆ = and nuF-K respacing to ∆ = . response from SDH measurements in Yogyakarta.

The dominant P wave velocity from the Yogyakarta measurements was 1,800 m/s, shown by the red line in the nuF-K response (Figure 9). This dominant velocity appeared from a distance of 16 to 31 m and could be predicted by its shape similarity. This behavior describes the same low frequency (around 30 Hz) source generation that interacts with the lower layer before reaching the sensor, which could be reduced by using nuF-K filtering. It also produced higher apparent velocity at around 1,800 m/s. Since the F-K energy was calculated from the signal's amplitude, the first arrival waves had a smaller amplitude than the interaction of waves that came after.

In contrast, the different shape of the arrival waves in shallow records shows partial surface layer interaction, where each sensor position receiving upward wave reflection at a different delay (negative velocity) that is not affected by F-K filtering in the positive range (Figures 10 and 11). Without filtering, the S wave arrival in each distance trace is hardly visible. Applying a notch velocity filter to remove the apparent velocity at 1000 − 2200m/s in the uniform F-K response (Figure 10Error! Reference source not found.), some polarity reversal of S

wave is visible and more can be acquired by applying the filter to the nuF-K response (Figure 11).

Figure 10 The S wave arrival time is calculated by finding the first arrival with reversing the polarity response between the S1 and the S2 data. F-K filtering is applied to reduce the P wave response, which may also give better visuals of P wave arrival, marked by a red cross. The visible S wave arrival in the notch filtered () response is marked by a blue cross. The green line and red line are from the three-layered model as a reference.

Figure 11 Notch filtered response to nuF-K transformed data produced better visual detection of S-wave arrival.

The depth versus velocity curve was acquired by recalculating the depth (ℎ) from source to sensor distance . By using a trigonometrical equation, the real-time delay between depths of measurement was calculated. The velocity response was validated by using bore-log data acquired from the same borehole, consisting of the layer description and the N-value from the standard penetration test (SPT) correlated to the S wave velocities [26] in Figure 12.

Figure 12 Correlation between acquired N-SPT value with calculated velocities for different types of soil layer.

4 Discussion

Non-uniform F-K transform for depth-to-distance conversion in SDH processing results in the correct velocity response of the F-K image. When the shape of the arrival waves is consistent along the measurement depth, the mean velocity response can be single-picked directly from the F-K image without detecting the arrival time for each depth. However, the main wave components (i.e., the P and S components) may not be easily distinguishable if the source has a very lowfrequency signal penetrating a high-velocity layer. Trace normalization is important to improve continuity in depth or distance to wave number transform to reduce the shadowing effect. The shadow response will create imperfect filtering since the shadow may overlap between filtered and unfiltered regions.

The second normalization step proposed in this paper uses the standard deviation of the tight band passed F-K response that will likely have the same wave shape in all distances. This normalization is best applied to the consistent shape of the wave group, i.e., the P wave component, which is not affected by sensor rotation inside the borehole. The first wave arrival shape may change drastically due to source and layer interaction; a comparison between the original data, the bandpassed and the notch-filtered signal is needed to select a successful main wave reduction filtering response.

5 Conclusion

The purpose of F-K filtering is to extract or reduce wave components in a specific range of velocities. The F-K filter modification described in this paper corrects the range and velocity response for seismic downhole processing. The hyperbolic distance between source and receiver producing an incorrect velocity response is corrected by implementing a non-uniform transform for depth-to-distance conversion. We also described a two-step normalization method for correcting

each trace's inconsistent energy and wave shape to reduce the F-K image's shadow response, thus correcting the filter response. However, the F-K filter will only be optimally applied if each wave component is distinguishable in the F-K image, far away from the F-K axes, i.e., when in the seismic downhole method each of the wave components has a huge difference in velocity and the source frequency is not too low.