regularized least-squares inversion

To force the energy in the SODCIGs to concentrate at the zero-offset location, we can pose the problem as a regularized inversion process, and the objective function is defined as follows:  

$$J({\bf m}) = \Vert {\bf W_d}({\bf Lm} - {\bf d}) \Vert _2 + f\left( D({\bf m}) \right),$$ (27)

where $\bf d$ is the recorded data. $\bf L$ is a 2-D/3-D wave-equation modeling operator that transforms the model to prestack data; here I use the adjoint of the Double Square Root (DSR) migration operator. $\bf W_d$ is a mask weight, which enables us to minimize the data residuals only at known locations, and $\bf m$ is the model space in terms of the SODCIGs, a 3-D image cube in the 2-D case and a 5-D image cube in the 3-D case. Operator $D({\bf \cdot})$ is defined as follows:

$$\begin{eqnarray} D({\bf m}) &=& {\bf diag}(\vert{\bf h}\vert) \bf m \ {\bf diag... ...=& {\bf diag}(\vert h_1\vert,\vert h_2\vert,\cdots,\vert h_M\vert)\end{eqnarray}$$ (28) (29)

which is the DSO operator acting along the offset dimension to penalize energy far from zero-offset locations. Near-offset energy, especially that around zero-offset locations where $\vert{\bf h}\vert\approx 0$, will not be affected. After applying DSO, the model-dependent sparseness transform operator $f(\cdot)$, which minimizes model residuals in the L₁ norm or Cauchy norm, is performed. The sparseness constraint is applied depth-by-depth and CMP-by-CMP along the offset dimension. The purpose of adding such a sparseness constraint is to penalize noise which is typically incoherent and weak, and consequently enhance the resolution of the final inverted result.

In fact, adding the DSO regularization term in the SODCIGs is similar to adding a roughener to smooth along the offset-ray parameters in the ADCIGs. As offset-ray parameters are connected to the offset wavenumbers via the following equation:

$${\bf p_h} = \frac{{\bf k_h}}{\omega},$$ (30)

for a single frequency,

$$\frac{\partial}{\partial {\bf p_h}} = \omega \frac{\partial}{\partial {\bf k_h}},$$ (31)

therefore, for a specific CMP location $\bf m_0$, roughening along ray-parameters can be expressed as follows:

$$\begin{eqnarray} 0 & \approx & \frac{\partial}{\partial {\bf p_h}}P(\tau,{\bf m_... ...{\partial}{\partial {\bf k_h}} P(\omega, {\bf m_0}, {\bf k_h}, z),\end{eqnarray}$$ (32) (33) (34)

where S is the slant-stack operator and ${\mathcal F}^{-1}_{\omega}$ is the inverse Fourier transform over frequencies. Considering the Fourier duality, convolving the wavefield $P(\omega, {\bf m_0}, {\bf k_h}, z)$ with the differential operator in the offset-wavenumber domain has the same effect as multiplying the wavefield $P(\omega, {\bf m_0}, {\bf h}, z)$ with ${\bf diag}({\bf h})$ in the offset-space domain:

$$\frac{\partial}{\partial {\bf k_h}}P(\omega, {\bf m_0}, {\bf... ...rightarrow {\bf diag}({\bf h})P(\omega, {\bf m_0}, {\bf h}, z).$$ (35)

So smoothing along the offset-ray parameters acts like the DSO regularization term.

Following the discussion in the previous section, the objective function (27) can be rewritten as follows:

$$J({\bf m}) = \Vert {\bf W_d(Lm - d)} \Vert _2 + \Vert {\bf W_s}D({\bf m}) \Vert _2,$$ (36)

or more concisely in terms of fitting goals,

$$\begin{eqnarray} 0 & \approx & {\bf W_d(Lm-d)} \ 0 & \approx & \epsilon {\bf W_s} D({\bf m}) ,\end{eqnarray}$$ (37) (38)

where $\bf W_s$ is the diagonal weighting matrix that forces the sparseness constraints. Here I use the Cauchy norm, thus

$$W_{si} = \frac{1}{\sqrt{1+ ( D(m_i)/\sigma)^2}}.$$ (39)