Structure-adapted Multichannel Matching Pursuit for Seismic Trace Decomposition

While the single-channel matching pursuit decomposes a seismic trace into a series of wavelets, the multichannel matching pursuit examines the lateral coherence of the seismic traces as a constraint to improve the lateral continuity of the decomposition. However, the presence of structures in the subsurface negatively affects the performance of the multichannel matching pursuit. We proposed a structure-adapted matching pursuit method that combines with a dynamic time-warping (DTW) algorithm to estimate the similarity between adjacent seismic traces and extract an optimal wavelet along the dip plane. This structure-adapted implementation would significantly speed up the convergence of the decomposition process. We modified the DTW algorithm by combining the Euclidean distance of the seismic trace and the first-order temporal derivative of the seismic trace. We also updated the amplitudes of all extracted wavelets simultaneously based on the least-squares principle. This DTW-based, structure-adapted, least-squares, multichannel matching pursuit method would improve the robustness and accuracy of seismic trace decomposition.


Introduction
In matching pursuit, a single seismic trace is decomposed into a stack of wavelets, with each optimal wavelet selected to have the highest correlation coefficient with the seismic trace (Wang, 2007). These decomposed wavelets are represented in the time-frequency domain generated by the Wigner distribution, and the time-frequency spectrum has a relatively high resolution compared to conventional time-frequency analysis techniques. Therefore, the matching pursuit method has gained great popularity in seismic applications, including data compression (Hu et al., 2015), non-stretching trace correction (Zhang et al., 2013), wavefield reconstruction (Ozebk et al., 2010;Vassallo et al., 2010) and hydrocarbon reservoir detection (Xue et al., 2020;Zeng et al., 2017).
Wavelets are selected from a given dictionary, and hence matching pursuit is inherently a greedy solution (Mallet & Zhang, 1993). In practice, the size of the given dictionary can be minimized by assuming that the seismic trace can be represented by specified wavelets. Wang (2007) suggested using the analytical forms of the Morlet wavelet. Liu and Marfurt (2007) tried to speed up the process by extracting multiple basic wavelets in parallel computation. Liu et al. (2014) proposed a method to improve the efficiency of the matching pursuit by using a fast exponential matching pursuit based on the Ricker wavelet. Liu and Zhang (2017) extended the parallel computation by considering the number of envelope peaks versus the number of computer nodes. Semnani et al. (2019) introduced the quantum swarm evolutionary matching pursuit to increase the efficiency of the conventional matching pursuit.
In addition to the high computational cost, the matching pursuit suffers from the non-uniqueness of decomposition and is very sensitive to random noise in the seismic trace. To overcome these limitations, Wang (2010) proposed a multichannel matching pursuit (MC-MP) method that uses lateral coherence as a constraint in the decomposition. The MC-MP method is based on linearity theory, which states that the basic wavelet shared by a group of adjacent traces is the best match with the average of these local traces. The MC-MP method has been shown to improve lateral continuity and suppress random noise.
In MC-MP, the local seismic trace is based on an average of the local neighbouring traces (Wang, 2010). However, time delays of wavelets in neighbouring traces lead to a loss of accuracy of the average and make wavelet extraction more complicated (Durka et al., 2005). Moreover, suboptimal estimation of the wavelet leads to more iterations in the decomposition and inconsistent lateral time-frequency spectra. Therefore, in this paper, we propose a method to implement MC-MP in a structure-adapted way.
The dynamic time-warping (DTW) algorithm can be used to estimate the similarity between two time series and is widely used in geophysics (Chen et al., 2018;Hale, 2013;Venstad, 2014;Yao & Wang, 2022). However, conventional DTW only uses the Euclidean distance between two seismic traces and ignores the trend. Keogh and Pazzani (2001) used the derivative of the time series to accelerate the convergence of DTW. Choi and Kim (2018) modified the DTW based on directional similarity. In this paper, we modify the DTW algorithm by combining the Euclidean distance and trend information from two seismic traces.
As with the standard matching pursuit, the residual in each iteration is only orthogonal to the selected basic wavelet (Mallet & Zhang, 1993;Tropp and Gilbert, 2007). To speed up the convergence of the matching pursuit, we update the amplitudes of all wavelets simultaneously using the least-squares principle instead of calculating the amplitude of the extracted wavelet traces individually. Pati et al. (1993) suggested projecting the residual trace orthogonally onto the individual wavelets. Feng et al. (2017) and  applied this method to seismic trace decomposition and called it orthogonal matching pursuit. However, the term 'orthogonal' seems misleading here, as  points out that a currently selected wavelet is not necessarily orthogonal to the previously selected wavelets. We would like to emphasize that the simultaneous calculation of the wavelet amplitudes is based on the principle of least squares. Due to the simultaneous least-squares estimation, either in a single or multichannel implementation, the proposed method significantly improves the convergence of the decomposition.
In this paper, we first define the least-squares matching pursuit method. Then we present the structure-adapted multichannel implementation for matching pursuit. The best match of the basic wavelet in a local neighbourhood is obtained by applying the DTW algorithm to the original seismic traces. We demonstrate the validity and effectiveness of the proposed algorithm through studies of coal-seam reflection extraction and time-frequency spectrum computation in a complex subsurface structural environment.

The Least-Squares Multichannel Matching Pursuit
In the single-channel matching pursuit, a seismic trace sðtÞ is composed of N wavelets: where g c i is a basic wavelet (or atom) defined by a set of parameters c i , a i is the amplitude of the ith wavelet, and R ðNÞ fsðtÞg is the residual after extracting N wavelets. Matching pursuit decomposes a seismic trace into individual wavelets, starting at i ¼ 1, for which the original seismic trace is denoted as R ð0Þ fsðtÞg ¼ sðtÞ.
The basic wavelet g c i can be defined by the Morlet wavelet (Wang, 2007): where x m is the mean angular frequency, r is a parameter controlling the width of the wavelet, / is the phase, and t 0 is the time delay of the Morlet wavelet. Thus, the set of parameters defining a basic wavelet in Eq. (1) is where ; h i is the inner product operator, and g c i t ð Þ is the L 2 norm of g c i ðtÞ. Each wavelet is thus defined by a i ; c i f g. The decomposition is performed iteratively, and each iteration extracts a single wavelet.
The least-squares matching pursuit differs from the basic matching pursuit in that the amplitude a i is estimated. In the least-squares matching pursuit, not only the amplitude of the current wavelet is estimated, but also the amplitudes of the previously estimated wavelets are refined. The estimation and the refinement are performed simultaneously in the least-squares sense.
If we list N basic wavelets in matrix U N t ÂN , where each column is a basic wavelet, and list the corresponding amplitudes in vector a NÂ1 , we can represent N wavelets in a discrete form as U N t ÂN a NÂ1 . The implementation of the least-squares matching pursuit in a single iteration can be summarized as follows: 1. Extract the basic wavelet g c i that has an optimal correlation with the residual trace R iÀ1 where [ is the argument operator which combines matrix U N t ÂðiÀ1Þ and vector g c i to form matrix U N t Âi .
3. Estimate the amplitudes of all selected wavelets by minimization: Minimization not only estimates the amplitude of the current ith wavelet but also refines the amplitudes of all the previous (i -1) wavelets. The solution estimate in a least-squares sense may be expressed as where U N t Âi ð Þ y is the pseudo-inverse of U N t Âi , calculated by the SVD method. 4. Update the trace residual as The iteration stops once the preset criteria are satisfied.
In a least-squares solution, the optimization based on the minimum variational method may be expressed as This minimum variation indicates that the residual is orthogonal to all the selected wavelets U N t Âi . Since the wavelets are not necessarily orthogonal to each other, we refer to this method as least-squares matching pursuit, and not misleadingly as orthogonal matching pursuit.
We implement least-squares amplitude estimation in the MC-MP which was originally proposed by Wang (2010). The idea is to use local coherence in seismic data to improve lateral continuity in decomposition. The matched basic wavelet is the one shared by locally adjacent traces in a small group.
The initial estimate of a basic wavelet is extracted from an average trace of the local neighbouring traces. Based on the linear theory (Wang, 2010), the residual of an average trace y t ð Þ is the average of all residual traces: where L is the number of traces in a local group. The group of seismic traces can be expressed in matrix form as where S Nt ÂM trace is the matrix of M trace seismic traces, in which each column is a seismic trace vector, U N ð Þ is the dictionary matrix containing the atoms extracted in N iterations, A N ð Þ is the matrix of amplitudes in which each column is the amplitude of a wavelet corresponding to a trace, and R N ð Þ is the matrix of residual traces. The amplitude is calculated in the least-squares sense as follows: Vol. 180, (2023) Structure-Adapted Multichannel Matching Pursuit 853 where Á k k is the L2 norm of a matrix. The residual of the ith iteration is updated as follows: Figure 1a shows a group of synthetic seismic traces. Figure 1b, c compare the frequency spectra at 15 Hz estimated by the signal-channel least-squares matching pursuit method and the proposed MC-MP method. The decomposition ends when no new atom can be extracted from the seismic traces. The details of time-frequency spectra generation can be found in Wang (2007Wang ( , 2010. Figure 1b shows an obvious discontinuity between 2.0 and 2.5 s around trace 4, while Fig. 1c shows no discontinuity. Figure 1 clearly shows the advantage of spectral continuity in the time-frequency spectrum produced by leastsquares MC-MP.

Structure-Adapted Multichannel Implementation
The MC-MP method from the previous section uses lateral coherence but does not consider spatial structure. Here, we propose a structure-adapted multichannel implementation to improve the efficiency of the MC-MP method.
The DTW algorithm allows us to find the best match between two seismic series without explicitly estimating the local dip information (Chen et al., 2018;Hale, 2013;Venstad, 2014). The DTW between two seismic series can be expressed as follows: where q denotes the power, d denotes the cost function between two time series, the alignment path p is a sequence of K index pairs and Aðs 1 ; s 2 Þ is a set of all admissible paths. The admissible paths include paths that satisfy the following conditions: The problem with the conventional DTW algorithm is that it only considers the difference between the values of two seismic traces, while ignoring the trend of the seismic traces. Choi and Kim (2018) proposed the derivative DTW, which uses the trend of the data. However, the derivative DTW ignores the similarity between successive traces. Moreover, the derivative DTW is too sensitive to the noise in the seismic trace. Figure 2c shows the result of matching with the derivative DTW. It can be seen from the figure that the two traces are not correctly matched with the derivative DTW. Here we propose a method to modify DTW by combining the Euclidean distance value and the trend of the time series: where a is a balancing factor between the DTW of the original seismic traces and the DTW of the derivative of seismic traces. Figure 2d shows the matching result estimated with the proposed algorithm, which shows that the traces are correctly aligned. Figure 3 shows an example of how the mDTW is performed between two signals. These two seismic series are noisy traces with random noise with a level of 10 dB. Figure 3a shows that there is a singularity problem between the two traces, where a single point matches multiple points in another trace (in the red ellipse). Figure 3b shows the result obtained with the modified DTW algorithm, which shows that the singularity problem has been largely improved.
In summary, the implementation of the structureadapted, multichannel, least-squares matching pursuit method consists of the following steps: 1. Apply the modified DTW algorithm to the seismic traces of a group. 2. Determine the optimal basic wavelet using the weighted average along the local dip plane. 3. Update the dictionary of extracted atoms U i . 4. Based on the extracted dictionary of atoms, calculate the amplitudes of all basic wavelets in the group of unwarped seismic traces. 5. Calculate the updated residual energy. 6. Iterate until a predetermined criterion is met.

Applications
We apply the proposed structure-adapted leastsquares MC-MP method to both synthetic and field seismic data to test the validity and efficiency of the proposed algorithm. Figure 4a shows a series of synthetic seismic traces with constant time delays between adjacent traces. Figure 4b shows the wavelets extracted from the average trace without structure-adapted implementation. This simple average implementation causes not only an amplitude change, but also a phase distortion, and therefore makes the wavelet extraction inaccurate and time-consuming. An obvious proof of this is that after 100 iterations there are residuals between 500 and 600 ms (in red ellipses). The possible reason for the slow convergence of the original matching pursuit after 100 iterations is that the wavelets have cancelled out with the simple average and it is difficult to extract them from the residual traces. The extracted wavelets can be seen in the middle of Fig. 4b, which are not satisfactory. The right side of Fig. 4b shows the residuals of the original traces estimated without structure-adapted implementation within 100 iterations. This clearly shows that a large residual remains after 100 iterations.

Synthetic Study
In Fig. 4c, the wavelets extracted with the structure-adapted implementation match the original data very well. The wavelets at about 500 ms of the first trace are clearly separated. The residual traces obtained with the structure-adapted implementation show much better completeness within a smaller The structure-adapted implementation also accelerates the convergence of the decomposition. Figure 5 shows the normalized residual energy as a function of the number of iterations. Within 20 iterations, the residual energy in MC-MP with the structure-adapted implementation falls below -40 dB (the triangular dotted line). Without the structure adjustment, however, the number of iterations required for the normalized residual energy to fall below -20 dB is about 100. This makes sense because the wavelets extracted with the structureadapted implementation match the original trace better.

Coal-Seam Reflection Removal
The coal bed is always a good indicator for structural interpretation. However, the strong amplitude of the coal seam interferes with the weak reflection, making seismic interpretation difficult. We test the effectiveness of the proposed algorithm in decomposing seismic data with strong coal-seam reflections. Figure 6a shows the original data set consisting of 200 traces, recorded between 1200 and 1800 ms. This data set has very strong coal-seam reflections, close to the nearby reflections at about 1500 ms. This presents a challenge for detecting weak reflections below the coal-seam reflections. In the original data set, it is difficult to distinguish them from the neighbouring reflections. We apply the proposed algorithm to this data set and attempt to extract the strong coal-seam reflections from the original data set. Figure 6b shows the coal-seam reflections removed by the proposed algorithm. Figure 6c shows the area after we have removed the coal-seam reflections. It can be concluded that the proposed algorithm successfully removes the coalseam reflections while maintaining the lateral continuity of the seismic data set. To illustrate the advantage of the proposed algorithm in terms of effectiveness and efficiency, we show the extraction result with the conventional MC-MP in Fig. 6d. The comparison between Fig. 6c and d shows that the proposed algorithm performs better in extracting strong coal-seam reflections.

time-frequency Spectrum Estimation in Complex Structures
time-frequency spectra are commonly used in seismic data processing and interpretation to detect anomalies not found in the time domain. Now, we demonstrate the ability of the proposed algorithm to generate the time-frequency spectrum in an area with complex structures (Fig. 7). The data set has a time sampling interval of 4 ms. For simplicity, we show the range that includes the time between 1400 and 1848 ms for line numbers between 2.0 and 3.0 km. From Fig. 7, we can see that this range has curved structures under the surface, which pose a challenge to the conventional MC-MP algorithm. Figure 8 shows four frequency spectrum profiles estimated with the proposed algorithm. From the comparison of the spectrum profiles with the instantaneous envelope, it can be easily seen that a highresolution time-frequency spectrum with high lateral continuity can be obtained with the proposed algorithm in areas with complex structures.
Although the algorithm is more efficient in areas with complex structures, it involves an additional computational cost. One should limit the application of the structure-adapted implementation to domains with complex structures.

Conclusion
We proposed the structure-adapted, least-squares, multichannel matching pursuit (MC-MP) for seismic trace decomposition for cases with complex structures, which was also combined with a modified version of dynamic time warping. The modified DTW combines the derivative DTW and the standard DTW to achieve better matching of successive traces. Although the wavelet extraction is sequential, when a basic wavelet (or atom) is extracted, the amplitudes of all already extracted basic wavelets are updated simultaneously using the least-squares principle. This implementation has improved the efficiency and robustness of decomposition. While the MC-MP improves the uniqueness of wavelet extraction and the continuity of time-frequency spectra by exploiting the lateral coherence of local neighbouring traces, the structure-adapted implementation further improves the efficiency and convergence rate of