Time-frequency analysis of seismic data using local attributes

SHAPING REGULARIZATION FOR INVERSE PROBLEMS

In this appendix, we review the theory of shaping regularization in inversion problems. Fomel (2007b) introduces shaping regularization, a general method of imposing regularization constraints. A shaping operator provides an explicit mapping of the model to the space of acceptable models.

Consider a system of linear equation given as $\mathbf{Ax} = \mathbf{b}$, where $\mathbf{A}$ is a forward-modeling operator, $\mathbf{x}$ is the model, and $\mathbf{b}$ is the data. By equation 8, we can find that the proposed time-frequency decomposition can be written as the form $\mathbf{Ax} = \mathbf{b}$. The standard regularized least-squares approach to solving this equation seeks to minimize $\left\vert \left\vert \mathbf{Ax} = \mathbf{b} \right\vert \right\vert _{2}^{2}... ...epsilon ^{2} \left\vert \left\vert \mathbf{Dx} \right\vert \right\vert _{2}^{2}$, where $\mathbf{D}$ is the Tikhonov regularization operator (Tikhonov, 1963) and $\varepsilon$ is scaling.

The formal solution, denoted by $\hat{\mathbf{x}}$, is given by

$$\hat{x}=\left( \mathbf{A}^{T}\mathbf{A}+\varepsilon ^{2}\mathbf{D}^{T}\mathbf{D}\right)^{-1}\mathbf{A}^{T}\mathbf{b},$$ (14)

where $\mathbf{A}^{T}$ denotes the adjoint operator. Fomel (2007b) defined a relation between a shaping operator S and a regularization operator $\mathbf{D}$ as

$$\mathbf{S}=\left(\mathbf{I}+\varepsilon ^{2}\mathbf{D}^{T}\mathbf{D}\right)^{-1}.$$ (15)

Substituting equation A-2 into equation A-1 yields a formal solution of the estimation problem regularization by shaping

$$\mathbf{\hat{x}}=\left(\mathbf{A}^{T}\mathbf{A}+\mathbf{S}^... ...bf{I}\right)\right]^{-1}\mathbf{S}\mathbf{A}^{T}\mathbf{b}.$$ (16)

Introducing scaling of $\mathbf{A}$ by $1∕λ$ in equation A-3 , we obtain

$$\hat{\mathbf{x}}=\left[\lambda ^{2}\mathbf{I}+\mathbf{S}\le... ...bf{I}\right)\right]^{-1}\mathbf{S}\mathbf{A}^{T}\mathbf{b}.$$ (17)

The conjugate-gradient method can be used to compute the inversion in equation A-4 iteratively. As shown by Fomel (2009), the iterative convergence for inversion in equation A-4 can be dramatically faster that the one in equation A-1.

The main advantage of shaping regularization is the relative ease of controlling the selection of $\lambda$ and $S\mathbf{S}$ in comparison with $\varepsilon$ and $\mathbf{D}$ cite[]Fomel2009. In this paper, we choose $\lambda$ to be the median value of $\mathbf{\Psi_{k}(t)}$ and the shaping operator to be a Gaussian smoothing operator. The Gaussian smoothing operator is a convolution operator with the Gaussian function that is used to smooth images and remove detail and noise. In this sense it is similar to the mean filter, but it uses a different kernel to represent the shape of a Gaussian (bell-shaped) hump. The Gaussian function in 1D has the form,

$$G(x)=\frac{1}{\sqrt{2 \pi \sigma}} \exp \left( -\frac{x^{2}}{2\sigma^{2}}\right).$$ (18)

One can implement Gaussian smoothing by Gaussian filtering in either the frequency domain or time domain. Fomel (2007b) shows that repeated application of triangle smoothing can also be used to implement Gaussian smoothing efficiently, in which case the only additional parameter is the radius of the triangle smoothing operator. In this paper, we use the repeated triangle smoothing operator to implement the Gaussian smoothing operator.

The computational cost of generating a time-frequency representation with our method is $O\left(N_{t}N_{f}N_{iter}/N_{p}\right)$, where $N_{t}$ is the number of time samples, $N_{f}$ is the number of frequencies, $N_{iter}$ is the number of conjugate-gradient iterations, $N_{p}$ is the number of processors used to process different frequencies in parallel, and $N_{iter}$ decreases with the increase of smoothing and is typically around 10.

Time-frequency analysis of seismic data using local attributes