Simultaneously computing the 3D dip attributes for all the time samples of a seismic trace

The seismic dip attribute is regularly used to aid structural interpretation and is commonly adopted as a compulsory input for computing other seismic geometric attributes. One disadvantage of current dip computation algorithms is that interpreters compute the dip attribute time sample by time sample and do not consider the relationship between dip values of nearby samples. The classic convolution theory suggests one formation boundary should have the corresponding seismic event. However, the seismic wavelet always has a certain time duration. As a result, one formation boundary has a corresponding seismic event that consists of several time samples. Ideally, the time samples, which belong to the same boundary, should have approximately the same dip attributes. In this research, a sample by sample computation procedure is treated as an independent optimisation procedure. Then, simultaneously computing the seismic dip of time samples of one seismic trace can be regarded as a multi-objective optimisation procedure. The proposed method is based on analysing features of seismic waveform within user-defined windows. Considering that nearby time samples should have continuous dip values, we the dynamic time warping to simultaneously compute seismic reflectors’ dip values of a seismic trace. We applied our method to a field seismic data to demonstrate its effectiveness.


Introduction
Interpreters (e.g. Hoecker & Fehmer 2002;Qi et al. 2014;Alfarraj, et al. 2017;Batnagar et al. 2019) have used the seismic volumetric dip attribute in interpreting geological objects at the subsurface, such as chaotic, slumps and fan deposition, and discontinuities, such as faults and channels. Other seismic geometric attributes such as coherence and curvature need a seismic dip attribute as one of the compulsory inputs (Marfurt et al. 1998). The seismic dip is computed by comparing the seismic amplitude or derived seismic attributes between nearby seismic traces (e.g. instantaneous phase) within a user-defined window. There are two main categories for computing seismic dip attributes. The first category analyses the seismic amplitude within a user-defined window. Bahorich & Farmer (1995) computed seismic reflectors' dips by comparing cross-correlation coefficients between truncated seismic data by a set of preset time-lagging windows. Marfurt et al. (1998) proposed that cross-correlation coefficient-based analytical seismic traces are more robust to noise than those of seismic traces. Marfurt (2006) further improved the accuracy of dip value of seismic events nearby the discontinuities using Kuwahara's multiple windows searching. Based on the assumption that the value of seismic amplitude has minimal changes along the dip formation surface, Bakker et al. (1999) used the structure gradient tensor (GST), which is computed from seismic amplitude, to compute the seismic reflector's dip. To improve instability, Luo et al. (2006) computed seismic reflector dips and azimuths using a weighted structural tensor. Alfarraj et al. (2017) computed the seismic reflector dips by analysing the seismic waveform curvature/texture. However, Lou et al. (2019) found that it is very difficult for GST to compute an accurate seismic reflector's dip if the seismic survey has steep dipping seismic events. The second category analyses derived seismic attributes. Luo et al. (1996) computed seismic reflector dip by using GST analysis that was applied to the seismic instantaneous phase. Barnes (2000) pointed out that smoothing the instantaneous phase attribute before dip computation could further improve the stability of dip computation. Wang et al. (2019) further improved the accuracy of seismic dip by applying Kuwahara's multiple windows searching on the dip computed using an instantaneous phase attribute.
All the developed algorithms compute the seismic dip attribute time sample by time sample independently, and rarely consider the relationship between a seismic reflector's dip values that belong to the same seismic reflection event. As a result, the algorithms may produce artifacts that are associated with seismic stratigraphy boundaries if we use an improper temporal analysis window. In this research, we propose a novel algorithm to compute a seismic reflector's dip by considering the dip values that belong to the same seismic events. We treat the seismic reflector's dip computation of one analysis point as a single optimisation problem. Thus, simultaneously determining the dip of a seismic trace can be regarded as optimising multiple objectives problems. We begin with illustrating the two most popular dip computation methods: (i) dip computing by comparing cross-correlation coefficients between seismic traces truncated by a set of preset time-lagging windows, and (ii) dip computation using GST analysis. Then, we propose our objective problem, a constrained optimisation procedure. The optimising procedure chooses a set of dip values for all time samples along which we obtain the maximum value of summarised correlation cross-correlation coefficients. The constraint is that the dip variation of nearby time samples from a seismic trace should be smaller than a user-defined threshold. Next, we use dynamic time warping (DTW) to simultaneously compute a lagged time set for a seismic trace. Finally, the lagged time set is converted to the seismic reflector's dip in term of ms m −1 . To demonstrate the effectiveness of the proposed method, we applied it to the poststack seismic data from a F3 seismic survey acquired in New Zealand.

A review of algorithms for the seismic dip computation
Geologists and geophysicists prefer different ways to describe the orientation of a planar plane in the subsur- face. Geologists choose dip and strike, while geophysicists use inline dip and crossline dip that are oriented along the inline and crossline directions of a seismic survey to define the orientation of a planar surface (Qi et al. 2014). Figure 1 shows a schematic diagram illustrating a seismic reflector's dip computation using a discrete scanning algorithm. The yellow and blue lines (figure 1) represent userdefined dips and the preset extrema dips, respectively. First, the algorithm computes coherence using algorithms such as semblance, variance, principal components or some other statistical measurements along candidate dips. In this example, the discrete scanning strategy obtains the maximum coherence value when the truncated window is parallel with the orientation shown with a green line (figure 1). The dip indicated by the green line is regarded as the seismic reflector's dip for the point that is centered at the analysis window. Figure 2 illustrates a schematic diagram showing the dip computation using GST-based algorithm. First, the algorithm truncates the seismic trace using a user-defined window. Then, the algorithm constructs a two by two structure tensor matrix using the extracted seismic traces. Next, the algorithm computes the eigenvalues and vectors of the matrix. The first eigenvector (the eigenvector corresponding to the largest eigenvalue) is converted to the reflector's dip (the first eigenvector is regarded perpendicular to local seismic orientation).

The seismic dip computation using dynamic time warping
The current seismic dip computation algorithms compute the seismic reflector's dip time sample by sample. However, seismic amplitudes not only vary along the horizontal for the same seismic event, but they also vary along the time axis. Some samples have a strong seismic amplitude while others may have a weak seismic amplitude. As a result, the current seismic reflector's dip algorithm would generate an artifact 692 Figure 2. The seismic reflector's dip computation using a GST-based algorithm. The element g i is the gradient of the extracted data. The parameters associated with weak seismic reflection between two strong seismic events. Note that each formation boundary in subsurface should correspond to a seismic reflection, and each seismic reflection event consists of several seismic time samples. Thus, the time samples that belong to the same seismic reflection events should have similar dip values unless the seismic events were terminated by discontinuities such as unconformities or faults. As a result, the seismic dip should also demonstrate a gradual variation along the time axis, except the time sample located at the discontinuity locations. In this case, we propose to simultaneously compute the seismic dip for all time samples of a seismic trace. The objective can be expressed as follows: where j is time index of the seismic trace; coh [j, (p, q)] is the coherence value of the jth time sample for truncated seismic data along a local plane that has the inline dip p and crossline dip q; n is the total sample number of seismic trace; p(j) is the inline dip; q(j) is crossline dip and p thre and q thre are the user-defined threshold values. In this paper, we choose dynamic time warping to solve the optimisation problem defined in equation (1). Several researchers have successfully use DTW to solve geophysical problems such as flattening the seismic events (Jin et al. 2017) and archiving non-stretching normal moveout correction (Chen, et al. 2018). DTW has been applied to compute a seismic reflector's dip attribute: Arias (2016), Zhang et al. (2019) and Lou et al. (2021) have applied DTW to compute the seismic reflector's dip attribute. Let f (i) and g(j) denote two time series. In a signal processing field, DTW is commonly used to align two signals f (i) and g(j) as follows (Sakoe & Chiba 1978): where is a user-defined threshold that controls the stretching or squeezing rate; d(1 : N) is a sequence of integer time lag l that minimises the Euclidean distance between the referred and target signals as follows: To obtain the optimal sequence u(1 : N), DTW first recursively computes the accumulation error accu_e(i, l) from the alignment error e(i, l) as follows (Sakoe & Chiba 1978): Then, DTW obtains the minimising path, the shift sequence u(1 : N) through a backtracking step as follows (Sakoe & Chiba 1978): The trial time lag l(i) in equation (2) can be treated as the preset dip values shown in figure 1 (Arias 2016). We should obtain the minimum Euclidean distance between the seismic amplitudes of two seismic traces if time lag u(i) equals the dip 693 of seismic reflection event. Thus, we can use DTW to simultaneously find the optimal dip set for all the time samples of a seismic trace. Note that the optimal lag may not locate at the integer time index of the seismic trace. To obtain a fractional time lag, we fine interpolate both the refereed and target seismic traces L times before alignment. Then, we obtain the following relationship between seismic reflector's dip and time lag: where dt is the time increment and dx and dy are the grid interval along inline and crossline directions, respectively. Note that the seismic reflection event usually varies across the whole seismic survey due the horizontal changes of properties of subsurface such as impedance. Thus, the algorithm should compare the similarity between seismic waveforms instead of the error of Euclidean distance. In this paper, we use semblance to represent the degree of similarity between seismic waveforms (Qi et al. 2014). We should obtain a maximum coherence value if the truncated window is parallel with the local seismic reflection window (Marfurt 2006). To covert the 'maximum coherence' to the 'error' matrix that is needed by DTW, we replace the error of Euclidean distance e(i, l) using the equation defined as follows: where j is the time index of the seismic data; M t is the half size of the predefined defined time window in terms of sample number; m t stands the two-way travel time; f H is the imaginary part of analytic seismic data f (Taner, et al. 1979); x is the inline sequence number; y is the crossline sequence number; M is the total trace number along inline and crossline directions within the predefined window; p is candidate inline dip and q is the candidate crossline. Note that the maximum semblance is obtained between extracted seismic data if the truncation window is parallel with the local seismic reflector's dip (Marfurt 2006). Then, the dip computation problem defined in equation (2) is a maximisation problem and we convert this to a minimisation problem through equation (7) Note that the candidate semblance coefficients defined in equation (7) is a 3D matrix. To construct a 3D accumulated matrix for the back tracking, we recursively compare the semblance values as follows: accu_e (1, p, q) = e(1, p, q) accu_e (j, p, q) = e(j, p, q) Compared to the 2D accumulated matrix in the traditional DTW (figure 3a), our new accumulated error is a 3D matrix in the dip computation procedure (figure 3b). Figure 3 demonstrates that the new error accumulation procedure needs to compare the error values of seven neighbors that are different to the traditional error accumulation procedure (the traditional accumulation procedure only needs to compare the error values of three neighbors). Finally, we computed a set of optimal inline lag d p [1 : N] and crossline lag d q [1 : N] in the backtracking step: where Ω p and Ω q are the searching scope defined in equation (9). Note that the searching scope in the dip computation is a 3D sub-volume. Figure 4 parts a and b show the backtracking searching scope for the traditional DTW and modified DTW, respectively. Instead of comparing the accumulated error values of four neighbors, our new backtracking step needs to compare the accumulated error values of seven neighbors. In this paper, the inline dip inline and crossline dip crossline are defined in the unit of ms m −1 according to the computed optimal inline time lag set u p (1 : N) and crossline time lag set u q (1 : N) as follows: where dx and dy are the common mid-point grid intervals along inline and crossline directions, respectively.

Field application
The F3 seismic survey was acquired in North Sea and contains common typical geological structures such as salt domes, faults, folds and unconformities. The F3 seismic survey is used to demonstrate the effectiveness of our proposed method. The analysis window size along the horizontal 694  direction for GST is discrete scanning based, and our method is three inlines by three crosslines that is centered at the analysed seismic trace. Both GST and discrete scanning-based methods use a temporal window size of 20 ms. Figures 5  and 6 shows a representatively cropped inline seismic section and computed corresponding crossline dip using the proposed method, respectively. The size of temporal window in figures 6a and 4b are 4 and 20 ms, respectively. Note that our proposed algorithm produces the same results except the subtle differences pointed to by the white arrows in figure 6a and b. The phenomenon in figure 6 indicates that our proposed method is not sensitive to window size along the time axis. Thus, we set the window size along the time axis to 4 ms (one sample) in this paper. Figure 7 shows another representative inline seismic section. Figure 8a-d illustrates the computed crossline dip using the GST, discrete scanning, traditional DTW and proposed methods, respectively. The purple arrow in figure 8 indicates a channel. Note that the proposed method suc-cessfully produces inline dip anomalies on either side of the channel. However, both GST and discrete scanning methods only produce one anomaly for the channel. The red arrows in figure 8a-d indicate the locations where steeply dipping seismic reflection events exist. Note that GST fails to generate an accurate dip estimation for the steeply dipping seismic reflection events when compared to those of discrete scanning and proposed methods. The white arrows in figure 8 indicate the inline dip anomalies associated with faults. Note that the results of GST and discrete scanning are 'smeared' to a certain extent. However, our proposed method produces a 'sharp' anomaly for the inline dip attribute. Figure 8c and d indicates that the traditional DTW method (figure 8c) generates a lot of 'isolated' dip anomalies (representative 'isolated' anomalies are indicated by the black arrows) when compared to that of our proposed method (figure 8d). Figure 9 shows a representative time slice that is crossing a salt dome (the two-way travel time of the time slice is indicted by the yellow line in figure 7).   shows the time slice of computed inline dip using the GST, discrete scanning, traditional DTW and proposed methods, respectively. The red arrows (figure 10) indicate locations of salt features that can be better highlighted by dip attribute computation using our proposed method when compared to those of the GST and discrete scanning methods. The white arrows in figure 10 indicate the fault locations that are better highlighted by the computation using our method.  black arrows) when compared with that of proposed method. To further demonstrate the superiority of our method, we compare the curvature attribute that is computed from the seismic dip attribute (Chopra & Marfurt 2011). Figure 11 parts a-c show the time slices of most negative principal curvature computed from the dip attributes shown Figure 10a, b and d, respectively. The black arrows in Figure 11 indicate that the salt dome edges, which are heighted by our method, are missing in the results from the GST and discrete scanning methods.

Conclusions
In this paper, we treated seismic dip computation as a constrained optimisation problem and a novel method was proposed to simultaneously compute the seismic reflector's dip  of a seismic trace. To overcome the disadvantage of Euclidean distance in the computation of an error matrix, we proposed to use semblance to compute the error matrix for DTW. Semblance compared the similarity between seismic wave-forms. Thus, our error matrix was not sensitive to the seismic amplitude changes along horizontal direction. As a result, the proposed method was robust to noise and weak seismic reflection. Both the GST and discrete scanning methods need 698 interpreters to define a proper temporal window truncating the seismic traces. However, the real data example in this paper demonstrated that our method did not require a test in determining the temporal window size. The comparison of dip attribute on inline section and time slice illustrated that the computed dip attribute using the proposed method in this paper was superior in illuminating the fault locations and salt dome features. The comparison between derived curvature attribute further demonstrated that our method could produce higher-quality dip attributes than the GST and discrete scanning methods.