Wavefield extraction using multi-channel chirplet decomposition

In most acoustical and seismic applications, the wavefield is a combination of waves (e.g., direct arrival and reflections, normal modes, and compressional and shear waves) As knowledge of wave characteristics (i.e., time arrival, waveform, amplitude, and phase) is useful for solving the inverse problem (e.g., source localization and characterization, and estimation of environment parameters), wavefield extraction has always been a crucial issue. Wavefield decomposition can be achieved in different ways. A common way consists in projecting the signal on a set of functions. This set can be a basis of functions (e.g., time, frequency, time-frequency, Tau-P, frequency-wavenumber, and basis of matrix methods) where wave components overlap is small and upon which the filter is applied. For instance, the frequency domain achieved by the Fourier representation corresponds to the projection of the signal based on sinusoids. The set can also be a frame (overcomplete). In this case, the functions often correspond to a well-studied type (e.g., Gabor functions and wavelets) forming a dictionary of atoms and the projection is called atomic decomposition. The decomposition can be either adaptive, by searching a collection of different atoms in order to best match the inner structure of the signal, or not adaptive (i.e., fixed) as the well known continuous wavelet transform. After the decomposition, the signal can be written as a linear combination of atoms selected in the dictionary. Depending on the chosen method, this processing has limitations. First, changing the basis often requires a receiver configuration which consists of many sensors generally aligned and regularly spaced along an array. It is the case for frequency-wavenumber filtering,¹ Tau-P decomposition,² and matrix method of decomposition.³ Moreover, the new basis is not always successful to avoid overlap between wave components. This is particularly true for frequency and time-frequency domains, even if some methods attempt to prevent these intrinsic drawbacks.⁴ Those limitations prevent decomposition for numerous applications. For instance, in the case of reservoir monitoring using information derived from induced microseismicity, despite the short distance between the sources and the sensors positioned in boreholes, because the array is limited to a small number of sensors not necessarily regularly spaced, it is difficult to go beyond the classical location map. Moreover, the surrounding zone of microseismic sources is the place where body waves and dispersive waves are greatly interfering.

We propose in this paper a method to overcome these limitations. The objective is that after the decomposition each atom corresponds unambiguously to one wave component. For this reason, an adaptive decomposition on a large overcomplete dictionary is chosen. Two issues are linked to this type of methods: the choice of elementary functions (called atoms) that build up the dictionary and the algorithm of decomposition that selects from the dictionary the atoms that best suit and that calculates the projection coefficient. We use a dictionary of seven-parameter chirplets initially developed by Bardainne et al.,⁵ which allows a great adaptability. This dictionary is described in Sec. II. The algorithm of decomposition that we have developed includes an array processing based on the matching pursuit (MP) algorithm followed by an optimization step. It is described in Sec. III. Results on synthetic data are presented in Sec. IV.

In order to associate unambiguously one atom to one wave component, we sought a dictionary with two constraints. First, a wave component can be perfectly described by one atom of the dictionary. Then, the combination of several wave components (two or more) cannot be described by only one atom. The first assumption requires a dictionary extremely flexible and thus very large. We chose a seven-parameter chirplet dictionary, a chirplet being a time-limited function which corresponds to a piece of chirp. Like wavelets, chirplets are based on a mother function generally Gaussian. The two-parameter Gabor dictionary (i.e., time translated and frequency shift version of the mother function) or the wavelet dictionary (i.e., time translated and dilated/contracted) are the most conventional chirplet dictionaries.

A third parameter (i.e., the linear coefficient of modulation) has been introduced.⁶ Acoustic or seismic wave components can vary not only in arrival time (depending on travel path for a given wave source), in frequency (depending on source bandwidth and filtering attenuation by the medium), in duration (depending on source and dispersion), in linear or nonlinear frequency modulation (depending on source and dispersion) but also in overall waveform which characterizes the type of waves (e.g., body waves, guided waves, or interface waves). The symmetric Gaussian envelope of the mother wavelet is unsuitable for most of these waves. As a consequence, Gabor, wavelet, and even four-parameter dictionaries are not flexible enough to fully describe physical events and to fulfill assumptions formulated above. This is the reason why we used a normalized chirplet dictionary introduced by Bardaine et al.⁵ in a seismic context where atoms are described by seven appropriate parameters to suit the physical configuration, i.e., the time delay $τi$ and six morphological parameters. The first four parameters $fi$⁠, $oi$⁠, $q1i$⁠, and $q2i$ represent, respectively, the frequency shift, the chirplet duration (i.e., number of oscillations), the linear frequency modulation coefficient, and the nonlinear modulation coefficient. The last two parameters $e1i$ and $e2i$ control the envelope shape by acting, respectively, on the symmetry and on the duration of the constant flat shape between the left and right Gaussian halves. Once established the range of definition for each parameter, we build up the dictionary by accounting for all the possible combinations of the parameters. For details on the formulation of the chirplet dictionary, see Ref. 5. Finally, note that the larger the number of parameters and the range of solutions, the larger the computational complexity and time, but the more the assumptions can be fulfilled.

A decomposition associating unambiguously one atom to one wave component corresponds to a “sparse” problem because the number of atoms involved in the decomposition is “small” in comparison to the number of sample points needed to store the signals. Over the last years the signal processing community has been interested in this problem and algorithms have been developed to find sparse solutions. Unfortunately, the problem of constructing the best approximation of a signal as a linear superposition of $K$ atoms drawn from a coherent dictionary becomes an intractable combinatorial problem, even if the dictionary is finite. As a general rule, we face a trade-off between signal approximation and constraint relaxation. Therefore, it is addressed by heuristic stepwise algorithms which do not seek the optimal solution. The selected atoms and their coefficient of projection depend on the choice of the decomposition method. The well known MP algorithm introduced by Mallat et al.⁷ determines iteratively a set of atoms. Numerous methods have been developed to improve or to complete the MP method (e.g., basis pursuit⁸ and method of frames⁹). All these methods fail when one signal (i.e., one single sensor) is only considered. The originality of our decomposition method consists in taking advantage of the common information provided by a same wave (corresponding to the same physical phenomenon) at each sensor, leading us to extend the MP to the multi-component array configuration. Indeed, the basis pursuit is not extensible to array processing and the method of frames is not efficient in the case of sparse problem.

We first detail the decomposition of observed signals into a linear combination of $K$ chirplets by the MP. The starting point of the first iteration $(k=0)$ is the signal itself: $x0=x$⁠. For each iteration, the MP finds the best match between the signal and the atom dictionary $Dchirp={hi}$⁠, $i=1,…,M$⁠, by computing the maximum modulus inner product (i.e., the projection) of the signal among all atoms:

with the coefficient $Ak$ given by

and where $xk$ is the residual after the $kth$ iteration $xk=xk−1−Ak−1hk−1$⁠. After $K$ iterations, the signal is described by the linear combination:

where $RKx$ is the residual. This iterative procedure is performed until some predefined requirements are satisfied. The stopping criterion is generally the percentage of energy given by the residual between the decomposition result and the original signal. The MP is easy to implement besides being the fastest algorithm of decomposition. Note that if the seven-parameter chirplet dictionary is used, each selected atom $k$ is fully described by nine parameters, i.e., the seven of the associated chirplet $hk$ plus two parameters of the decomposition coefficient $Ak$ (i.e., the amplitude and the phase). To take into account the mutual information shared by sensor and component signals while preserving a degree of freedom to reflect wave changes, we propose the following two-step algorithm: a multi-component array MP process followed by an optimization process.

We first assume that a same wave is present on the whole array. The fulfillment of this assumption depends on the configuration of propagation and reception. This array can be divided into several parts to make sure this assumption stays valid. Depending on the configuration, we also assume that some of the seven parameters of chirplets describing the same wave across the array and some of the two projection parameters (i.e., amplitude and phase) are fixed, whereas the others vary independently. Still considering the wave, the same assumption is made for the outputs of a three-component sensor. The repartition between fixed and varying parameters is not necessarily identical but the set of fixed-sensor parameters must be included into the set of fixed-component parameters. The choice of fixed and varying parameters depends on the way of propagation of the waves through the complex media. Note that the amplitude and the phase are always supposed to be varying parameters across the array (for sensors and components). Parameters can thus be divided into three classes: the sensor and the component fixed parameters $F$⁠, the sensor and component varying parameters $V$⁠, and the mixed parameters $M$⁠. By assumption all the $M$ parameters are component fixed and sensor varying. The chirplet dictionary indexed by these classes of parameters is written $Dchirp={hF,M,V}$⁠. For instance, in a seismic context the chirplet is assumed to possess the same six morphological parameters across the array, whereas time delay and both the projection parameters are allowed to vary. For the outputs of a three-component geophone, only the projection parameters are allowed to vary. Finally, we can write $F={f,o,q1,q2,e1,e2}$⁠, $M=τ$⁠, and $V={|Ap|,arg(Ap)}$⁠.

The multi-component array matching pursuit (MAMP) algorithm implements the MP across the array, through the assumption on fixed and varying parameters. At the $kth$ iteration, the algorithm selects the chirplet which maximizes the double sum of projections:

where $L$ is the number of sensors and $J$ is the number of components at each sensor. The projection coefficient is given by

MAMP proceeds by iterating the same way as the conventional MP. The extracted wave components are characterized by their nine parameters and the relevant part of the signal is fully characterized in a sparse way by 9 K information.

Overlap between two independent wave components, which may skew the decomposition, is different at each sensor, making the array processing particularly relevant for this situation. The summation gets partially rid of this problem and makes the result more reliable. To converge toward a more reliable result, we propose a recursive process. We first explain the subroutine cycle 2 described in Fig. 1. The starting point is the initial wavefield (IW) on which the MAMP process is applied with two iterations to determine a first approximation of the two first wave components (noted $W1$ and $W2$⁠). This corresponds to the conventional MAMP process. To improve the MAMP efficiency, $W2$ is subtracted from the initial wavefield IW and an additional MAMP iteration is applied on the residual $(R2=IW−W2)$⁠, giving a more reliable approximation of $W1$⁠. Then, this new estimate $W1$ is subtracted from the IW and again, another MAMP iteration is applied on the residual $(R1=IW−W1)$⁠, giving in turn a more reliable approximation of $W2$⁠. This process is iterated until the energy of the residual, $R12=WI−W1−W2$⁠, no longer decreases. This ends cycle 2 which corresponds to the decomposition of the IW into two atoms. Similarly, the process can be extended to a finite number of atoms. For $N$ atoms, the recursive process is a generalization of cycle 2 by replacing $W1$ by cycle $N–1$ and $W2$ by $WN$ (Fig. 1). After this optimized step, all significant wave components are supposed to be extracted.

The original method of wavefield decomposition described in this paper was applied on a synthetic seismic data set. Observations were simulated by the OASES software¹⁰ for a simple configuration where an omnidirectional source emitted a pulse in a three homogeneous layer medium. The advantage of using the OASES code instead of real data in a first test is that the code provides exact solutions of wave components. The signal to be filtered is recorded in the vicinity of the interface between layers 1 and 2. It is composed of two compressional waves, up-P and down-P. The two-component (vertical and radial) sensors were vertically aligned but not regularly spaced. Some parameters of the wavefield model were variables, so that the interference between the two wave components might sufficiently change across the array (Fig. 2). The decomposition performances have been evaluated by measuring the percentage of energy to rebuild the true wavefield by the extracted one. Table I shows the variations of the average percentage on both components for the five sensors for each method used [here, MP, MAMP, and optimized MAMP (OMAMP)]. Thanks to the array processing and the optimization step, the optimized MAMP significantly improves the efficiency of separation process. The optimization step is particularly efficient in the case of a large overlap between wave components (traces 1 and 2). When the two wave components overlap so much that they visually seem to be single, the performances of the method are expected to be lower, as it is the case for the trace 1 (down trace of the Fig. 2). Nevertheless, the OMAMP algorithm still gives satisfactory results.

In this article, we have presented an original method of wavefield decomposition based on a chirplet decomposition within a two-step algorithm. This algorithm accounts for the array configuration to better constrain the uniqueness of the solution while optimizing the conventional decomposition algorithms. Results show a large improvement of the decomposition efficiency, even in the case of a large overlap between signals. Beyond the wavefield separation aspect, this signal decomposition is a key to an effective compact representation for signal compression when signals may be repeated over time as in acoustics or seismicity.

This work was supported by the Agence Nationale de la Recherche (ANR) under the EMSAPCO2 project.