Seismic data analysis using local time-frequency decomposition

Local time-frequency (LTF) decomposition

The Fourier series is by definition an expansion of a function in terms of a sum of sines and cosines. Letting a causal signal, $f(x)$ , be in range of $[0,L]$ , the Fourier series of the signal is given by

$$f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left[a_n\cos\left(\frac{2\pi nx}{L}\right)+ b_n\sin \left(\frac{2\pi nx}{L}\right)\right]\;.$$ (1)

The notion of a Fourier series can also be extended to complex coefficients as follows:

$$f(x) = \sum_{n=-\infty}^{\infty}A_n \Psi_n(x)\;,$$ (2)

where $A_n$ are the Fourier coefficients and $\Psi_n (x) = e^{i(2\pi nx/L)}$ .

Nonstationary regression allows the coefficients $A_n$ to change with $x$ . In the linear notation, $A_n(x)$ can be obtained by solving the least-squares minimization problem

$$\min_{A_n} \Vert f(x)-\sum_n A_n(x) \Psi_n(x) \Vert _2^2\;.$$ (3)

The minimization problem is ill posed because there are a lot more unknown variables than constraints. Our solution is to include additional constraints in the form of regularization, which limits the allowed variability of the estimated coefficients (Fomel, 2009). Tikhonov's regularization (Tikhonov, 1963) can modify the objective function to

$$\widetilde{A_n}(x) = \arg\min_{A_n}\Vert f(x)-\sum_n A_n(x)\Psi_n(x)\Vert _2^2 + \epsilon^2 \sum_n \Vert\mathbf{D}[A_n(x)]\Vert _2^2\;,$$ (4)

where $\mathbf{D}$ is the regularization operator and $\epsilon$ is a scaling parameter. One can define $\mathbf{D}$ , for example, as a gradient operator that penalizes roughness of $A_n(x)$ .

We use shaping regularization (Fomel, 2007b) instead of Tikhonov's regularization to constrain the least-squares inversion. Shaping is a general method for imposing constraints by explicit mapping the estimated model to the desired model, eg., smooth model. Instead of trying to find and specify an appropriate regularization operator, the user of the shaping-regularization algorithm specifies a shaping operator, which is often easier to design.

The absolute value of time-varying coefficients $\vert A_n(x)\vert$ provides a time-frequency representation, and equation 2 provides the inverse calculation. In the discrete form, a range of frequencies can be decided by the Nyquist frequency (Cohen, 1995) or by the user's assignment. In a somewhat different approach, Liu et al. (2009) minimized the error between the input signal and each frequency component independently. Their algorithm and the proposed algorithm are equivalent when the decomposition is stationary (or using a very large shaping radius), because they both reduce to the regular Fourier transform. In the case of nonstationarity, their approach does not guarantee invertability, because it processes each frequency independently.

Seismic data analysis using local time-frequency decomposition