A high-resolution processing technique for improving the energy of weak signal based on matching pursuit

This paper proposes a new method to improve the resolution of the seismic signal and to compensate the energy of weak seismic signal based on matching pursuit. With a dictionary of Morlet wavelets, matching pursuit algorithm can decompose a seismic trace into a series of wavelets. We abstract complex-trace attributes from analytical expressions to shrink the search range of amplitude, frequency and phase. In addition, considering the level of correlation between constituent wavelets and average wavelet abstracted from well-seismic calibration, we can obtain the search range of scale which is an important adaptive parameter to control the width of wavelet in time and the bandwidth of frequency. Hence, the efficiency of selection of proper wavelets is improved by making first a preliminary estimate and refining a local selecting range. After removal of noise wavelets, we integrate useful wavelets which should be firstly executed by adaptive spectral whitening technique. This approach can improve the resolutions of seismic signal and enhance the energy of weak wavelets simultaneously. The application results of real seismic data show this method has a good perspective of application.


Introduction
We usually assume, in seismic exploration, that the reflectivity series is sparse in time domain and a seismic trace can be modeled by convolving the reflectivity series with a wavelet.However, as wavelet spreads through different stratums, the amplitude, frequency etc. of wavelet will change; this lead to the assumption still has some deficiency and the conventional methods of improving resolution is not perfect.Actually, seismic data is nonstationary with varying frequency content in time domain, and the representation of it in frequency domain can illustrate many features which are difficult to visualize in time domain.Hence, in this paper, we try to develop a new method to increase the resolution of seismic data and enhance the energy of weak wavelets via decomposing the signal into a set of wavelets with different amplitude and frequency.
The technology of decomposition and reconstruction of seismic data could be used to analyse lithology and physical properties by taking advantage of the different characters of seismic data.Matching pursuit [Mallat andZhang 1993, Qian andChen 1993] is a kind of significant method which has been used to decompose a signal into a series of wavelets belonging to a comprehensive dictionary of functions like Gabor functions.The traditional atoms dictionary of matching pursuit algorithm was made up by scale, time delay, frequency, phase and amplitude.It was so big and redundant that can reduce the efficiency of decomposition of seismic trace.The search range of appropriate wavelets could be shrink by selecting a fixed scale parameter together with using complex-trace analysis to calculate attributes [Liu et al. 2004].Considering the Morlet wavelet is more suitable for energy and frequency quantification of seismic data and more appropriate for attenuation and resolution studies [Morlet et al. 1982a[Morlet et al. , 1982b]], they were employed as wavelets in the matching pursuit decomposition [Liu and Marfurt 2005].In additional, a useful way to improve the accuracy of matching pursuit is taking scale as a free parameter for making up more wavelets which were not considered by some of the matching pursuit implementation [Wang 2007].Even though taking scale as a free parameter can reduce residual, it also means that many more wavelets are included in the dictionary.
Matching pursuit was first used to analyze the seismic trace and it shows obviously advantages in timefrequency analysis [Chakraborty and Okaya 1995].And Marenco and Madisetti [1996] compared the difference between different time-frequency methods and used it to do deconvolution of seismic trace.Li et al. [1998] have integrated the matching pursuit with Kichhoff migration.In addition, it has been already applied to iden-tify gas reservoirs by low-frequency shadows in a timefrequency spectrum [Castagna et al. 2003, Sinha et al. 2005] and has been used in channel detection [Mohebian et al. 2013].In terms of improving resolution and enhancing the energy of weak signal, the differential resolution which uses a set of cascaded dipole filters to improve seismic resolution through spectral whitening [Claerbout 2004], a new algorithm [Sajid and Ghosh 2014] utilized a set of three cascaded difference operators to boost high frequencies together with a simple smoothing operator to boost low frequencies.In addition, an improved spectral whitening method [Lu et al. 2005] could be used to enhance the resolution of subdata obtained by SVD.
In this paper, we improve the matching pursuit algorithm by refining a local selecting range and increasing the accuracy of wavelets.In addition, we improve spectral whitening based on Lu [Lu et al. 2005] to enhance the resolution and the energy of weak signal.Results clearly illustrate the effectiveness of proposed algorithm.

Improved matching pursuit
The conventional matching pursuit algorithm is an expensive iteration procedure.It can decompose any signal into a linear expansion of wavelets by selecting parameters appropriately from a large and redundant dictionary with low efficiency.Inspired by Wang [2007], Liu and Marfurt [2005], we use Morlet wavelets to compose a wavelets family.Exploiting complex-trace analysis can obtain three complex-trace attributes so that we could shrink the search range of the time delay (u), the mean frequency (~) and the phase (U).Since scale (v) can be recognized as a free parameter which can control the width of a wavelet in time and the frequency bandwidth [Morlet et al. 1982b, Wang 2001], we take advantage of the level of correlation between constituent wavelets with different scale and average wavelet obtained from well-seismic calibration to design the search range of scale.
A Morlet wavelet m(t) is a kind of complex function in time domain which can be define as [Morlet et al. 1982b] (1) where ~m is the mean angular frequency, v is the scale parameter, u is the time delay and U is the phase.According to the seismic data is real signal, the Morlet wavelet can be written as (2) Let f as signal and D as a dictionary composed of Morlet wavelets.Matching pursuit is implemented iteratively, with the purpose of best matching the inner structures of signal, we compute a linear expansion of f over a set of vectors chosen from D. Each iteration can be regarded as a process of orthogonal decomposition.If the first optimal wavelet is m r0 , we then have an expression of f as where R (1) f is the residual vector after decomposing f in the direction of g r0 .The g r0 is orthogonal to R (1) f, so (4) If the iterations can meet the stop condition, after M iterations, the signal f could be expressed as (6) where a m is the amplitude of m rm , the final residual R (M) f is regarded as noise.During each iteration, a m as well as m rm can be determined by the following stages.
(1) Calculate the complex seismic trace by using Hilbert transform.The use of the complex trace p(t) make it possible to define instantaneous amplitude, phase, and frequency [Taner et al. 1979, Liu et al. 2007].The complex trace is (7) Calculate the envelope, instantaneous phase and instantaneous frequency of the complex seismic trace [Liu et al. 2004].We set the time of the maximum envelope of complex trace to be the time delay u n (8) set the instantaneous phase of the time delay u n to be the phase and set the instantaneous frequency to be the center frequency ~n (2) An average wavelet (Figure 1a) can be abstracted after well-seismic calibration.At the same time we construct a series of Morlet wavelets with fix u n , ~n, U n and different scale which is uniformly distributed by small step size.We calculate every correlation coefficient between the average wavelet and every Morlet wavelet (Figure 1b), and then obtain the changing trend of correlation coefficients.Figure 1b shows that scale parameter can dramatically impact the level of correlation coefficients.Specifically, the tiny variation of scale can bring the significant change in sidelobe of wavelets.After analyzing the changing trend of correlation coefficients, the search range can be constructed with the scale responding to the maximum correlation coefficient as the center of it and small step size as searching step.
(3) Once we obtain preliminary estimates of these parameters, we can construct a Morlet dictionary in user-defined time delay, frequency, phase and scale range [Gribonval 2001], where Du is the time-sampling interval, Dv is the sampling interval used for v, D~is the frequency-sampling interval, and DU is 5 degrees [Wang 2007].An optimal Morlet wavelet can be got by following equation during every iteration: (11) (4) Since m rm is orthogonal to R m+1 , the amplitude of the optimal wavelet m rm (t) can be estimated by (12) where normalizes optimal wavelet m rm .At last the final optimal wavelet g rm (t) is calculated by a m multiplying m rm .

Adaptive weak signal enhancement
According to the level of signal noise ratio of subspace to adjust the components of frequencies of signal adaptively, an adaptive spectral whitening method [Lu et al. 2005] can be used to enhance resolution of subdata obtained by SVD.We improved the method by analysing the advantages of it, and finally, formed an al-gorithm to strengthen the energy of weak wavelet got from matching pursuit.The spectral whitening [Gibson andLarner 1982, Coppens 1984], a basic method for spectrum balance, has been developed to enhance the bandwidth.It is an amplitude filter that can be realized in both time and frequency domains.The way realized in frequency domain, firstly, need to transform data in time domain to frequency domain, and then widen the bandwidth of the data in frequency domain, so the shortage of the whole process is lack consideration of the information in time domain.On the contrary, if the method will be executed in time domain, a band-pass filter is applied to filter the signal, and we sum the signal after balancing the amplitude of the signal.Even though this method takes the time information into account, it neglects frequency information [Wang 2012].
We develop an algorithm to adaptively enhance the energy of weak signal formed by a set of wavelets obtained from matching pursuit.The key step of it is to construct a spectral whitening filter [Li 2006].We need to calculate the complex seismic trace and get the envelope e(t), then the spectral whitening filter of the i wavelet is where o is the peak value of e(t), f(i) is the "whiten noise factor" of the i wavelet which is used to balance the energy of each wavelet and can be decided by the ( where f base and f wave are constants.The E i is the energy ration between the energy amplitude of i wavelet and the total energy amplitude of all wavelets in this trace: (15) If the seismic data s(t) can be decomposed into n wavelets, the instantaneous amplitude of a single wavelet g i (t), i =1,2,3…n, can be expressed as a i (t). a ~i(t) is the instantaneous amplitude estimated by a i (t), and then spectral whitening filter can be constructed point by point multiplication: (16) By multiplying a ~i(t) with Morlet wavelet got from matching pursuit, the reconstructed signal can be obtained by our method: (17) The flow chart of this paper is shown in Figure 2. The entire process can be summarized as three steps.
(1) Decompose the seismic data into a set of Morlet wavelets by matching pursuit.
(2) Calculate the instantaneous amplitude a ~i(t) to enhance the energy of weak signal component.
(3) Multiply the instantaneous amplitude a ~i(t) with Morlet wavelet to reconstruct new wavelet, and then sum these new wavelets as the final signal S ~(t).
The advantages of this method can be listed as following: (1) Using matching pursuit, it can not only decompose seismic data into a series of wavelets, but also can suppress spikes or random noise effectively.
(2) We use adaptive spectral whitening filter to every wavelet.It can consider the information both in time domain and frequency domain.

Application on a synthetic trace
Figure 3 illustrates the application of the algorithm to a synthetic trace.Figure 3a displays the synthetic trace comprised of additive noise and a set of wavelets with different frequency, amplitude, scale and time delay.We observe that, from the synthetic trace, the noise is obvious and the weak energy wavelets appear near the 400 sample point and below the 700 sample point.The synthetic trace is decomposed into a set of wavelets by using the matching pursuit in this paper and the result of decomposition is displayed in Figure 3b.
After removing the residual signal, we reconstruct the synthetic seismic trace with the wavelets (Figure 3c).By comparing Figure 3a and Figure 3c, we can see, using this method, the noise can be clean up effectively.If we use adaptive spectral whitening filter to every wavelet before reconstructing the synthetic trace, the result of reconstruction can be shown in Figure 3d.We can find the weak signal under the 700 sample point can be compensated obviously.
For a comparison, Figure 4 shows the time-frequency spectra of the signal in Figures 3a, 3c and 3d, respectively, which are generated by short-time Fourier transform with the same window [Gabor 1946, Qian andChen 1993].We can see the noise of original signal in time-frequency spectra.In addition, the energy near the 400 sample point and below the 700 sample point is too weak to be identified (Figure 4a). Figure 4b displays the time-frequency spectra of the reconstructed signal.We can see there is almost no noise in it.Therefore, it can be proved that our matching pursuit scheme may remove the noise effectively.If we apply adaptive spec-TECH TO IMPROVE THE WEAK SIGNAL VIA MP   tral whitening filter to every wavelet before reconstructing the signal, then we analyze the time-frequency spectra of the reconstructed signal.We can find that the resolution has been improved and the weak signal energy has been enhanced near the 400 sample point and below the 700 sample point (Figure 4c).Table 1 lists the energy ration, before and after applying adaptive spectral whitening to wavelets, between the energy amplitude of every wavelet and the total energy amplitude of all wavelets composing the trace.The number of wavelets is consistent with Figure 3b, and it can also prove that the adaptive spectral whitening could enhance the energy of weak wavelet effectively.
Figure 5 shows the application of synthetic seismic data.The synthetic seismic data is consisted of 60 traces and 1000ms per trace and the time interval of seismic data is 1ms. Figure 5a shows the seismic profile with additive noise.The seismic data includes two horizons, one is horizontal and the other is slanted.We can find the horizons are affected by the noise severely so that the relevant weak slanted horizon could not be recognized easily.The reconstructed profile using matching pursuit is displayed in Figure 5b and the residual error section is shown in Figure 5c.Obviously, the noise has been cleaned up effectively by the matching pursuit algorithm.After applying the adaptive spectral whitening filter to the reconstructed profile, the result can be showed in Figure 5d.It is easy to find that the energy of weak slanted horizon can be compensated effectively.

Application on real seismic data
The work area is in Junggar Basin, and the target formation belongs to Jurassic Sangonghe Formation J1s 2 1 Sand.The environment of deposition has been interpreted as Meandering river delta.Delta-front facies is the main sedimentologic facies in the reservoir [You Energy ration between the energy amplitude of every wavelet and the total energy amplitude of all wavelets (%)  Table 1.Energy ration between the energy amplitude of every wavelet and the total energy amplitude of all wavelets composing the trace.2006, Tang et al. 2009].The reservoir thickness is thin where sandstones and mudstones appears alternately in the vertical direction, and the area of sand changes rapidly in the horizontal direction, so it is difficult to determine the lapout in seismic section.Considering local accumulation of oil and gas can be formed here, some high-resolution work can be seen as the key step in this work area and it has practical meaning to develop oilfields.We apply our method to this work area.Figure 6a is a zoomed inline profile intersecting an oil well Zh1-2 through the original seismic data volume.The target layer is bounded between two horizons J1s1 and J1s 2 1 .We can see there is some noise in the seismic section and the resolution is so low that we cannot describe how the sands contact each other.Figure 6b shows the result obtained by our method.Here the number of wavelets to reconstruct the seismic trace could be determined by the experiment of matching pursuit.The parameters of our adaptive spectral whitening method are: f base = 0.01 and f wave = 0.001.In comparing with Figure 6a, we could observe that the weak energy of seismic data has been compensated effectively and the stratigraphic contact relation have become more clearly in the area cross the well.With the accuracy of it can be proved through the horizon calibrating of well Zh1-2, the stratigraphic trend and lapout could be determined relatively easily as well as precisely.Figure 7 illustrates the difference of the amplitude spectrum corresponding to the data shown in Figure 6.From Figure 6 and Figure 7, it is seen that our method not only can provide higher resolution, higher signal to noise ratio, more geology information, but also can compensate the weak energy of data.
Figures 8a, 8b and 8c show the profiles through the 25Hz, 35Hz and 45Hz isofrequency volumes corresponding to the seismic section shown in Figure 6a, respectively.The position of reservoir is indicated by Gamma etc. well-log curves.Even though the reservoir can be seen in the direction of the arrow in the three profiles, both the resolution and continuity of the data is poor.The 45Hz isofrequency section has higher resolution than others, but the existing of weak energy of signal and noise make the isofrequency section cannot display so well that we could conveniently identify the character of the reservoir.From Figure 8 we can see that the noise is abundant in original seismic data, especially in high frequency region.
The isofrequency profiles of 25Hz, 35Hz and 45Hz of the seismic section shown in Figure 6b

Conclusions
Matching pursuit is a powerful method to decompose a seismic trace into a series of wavelets.Considering the efficiency and accuracy of the iterative algorithm,   we have improved the algorithm by associating with well data from which we can confirm the central value of the search range of scale parameters and reduce the searching range.In addition, we have improved the adaptive spectral whitening through a special approach to get one of the "whiten noise factor" f(i).The signal could be reconstructed after the above process which is not only can remove the noise from the signal, but also can compensate the energy of weak signal and improve the resolution of stratum.The preliminary application results on synthetic and real data show the proposed method has a good perspective of applications.

Figure 1 .
Figure 1.(a) Average wavelet abstracted by well-seismic calibration.(b) Changing trend of correlation coefficients between average wavelet and Morlet wavelets with different scale, and the value of scale responding to the peak correlation is approximate 0.4.

Figure 3 .
Figure 3. Application of synthetic seismic trace: (a) original synthetic seismic trace with noise; (b) wavelets without noise decomposed by matching pursuit in this paper; (c) synthetic seismic trace reconstructed from (b); (d) synthetic seismic trace reconstructed from (b) after applying adaptive spectral whitening filter.

Figure 4 .
Figure 4. Short-time Fourier transform spectra of different signal with the same window: (a) time-frequency spectra of original synthetic seismic trace; (b) time-frequency spectra of the synthetic seismic trace reconstructed after using our matching pursuit scheme; (c) time-frequency spectra of the synthetic seismic trace reconstructed after matching pursuit and adaptively spectral whitening.

Figure 5 .
Figure 5. Application of synthetic seismic data: (a) noise seismic profile; (b) profile reconstructed by matching pursuit; (c) residual profile after matching pursuit; (d) result obtained by our method.
are exhibited in Figures 9a, 9b and 9c, respectively.Comparing Figure 8 and Figure 9, we found that the resolution of Figure 9 had been improved largely and the reservoir can be distinguished easily.The spatial continuity of 45Hz isofrequency profile became better than Figure 8c and the weak energy of signal compensated result in the reservoir more clearly.In addition, the noise can be removed well.
Figure 6.(a) Original seismic profile.(b) Seismic profile processed with our method.

Figure 7 .
Figure 7. (a) Amplitude spectrum of the original seismic profile in Figure 6a.(b) Amplitude spectrum of the seismic profile in Figure 6b (after processing with our method).