Localized rank-reduction method

Without the loss of generalization, we introduce the theory of the rank-reduction method in the case of a 3D problem. The 2D problem is just special case of the 3D problem with a single crossline trace. In the rank-reduction method, a pre-transformed matrix, e.g., the block Hankel matrix, is assumed to be of low rank. Then the goal becomes to extract the principal components of the pre-transformed matrix, i.e., the reflection signals, and separate the diffractions that usually corresponds to the inessential components due to lower spatial coherence.

Assuming that the tensor form of the input 3D seismic data is $D(x,y,t)$, the corresponding form in frequency-space domain is obtained after Fourier transform and can be expressed as $D(x,y,f)$. For each frequency slice, we map the frequency domain data into several block Hankel matrices. This process is referred as Hankelization (Siahsar et al., 2017; Oropeza and Sacchi, 2011; Huang et al., 2016). Given a frequency $f$, the Hankel matrix can be constructed as

$$\mathbf{H}_i(f)= \begin{bmatrix} D(i,1,f)&D(i,2,f)&\cdots&D(i,Q_y... ...vdots&\ddots&\vdots\\ D(i,P_y,f)&D(i,P_y+1,f)&\cdots&D(i,N_y,f) \end{bmatrix},$$ (1)

where the subscript $i$ is the row and $\mathbf{H}_i(f)$ is the Hankel matrix in frequency domain corresponding to the data matrix in frequency-space domain $D(x,y,f)$. $P_y=\lfloor N_y/2\rfloor+1$, $Q_y=N_y-P_y+1$ and $\lfloor \cdot \rfloor$ denotes the integer part of an input argument. $N_y$ denotes the number of traces in the $y$ direction.

Then, the corresponding block Hankel matrix can be inserted with equation 1

$$\mathbf{B}(f)= \begin{bmatrix} \mathbf{H}_1(f)&\mathbf{H}_2(f)&\c... ...hbf{H}_{P_x}(f)&\mathbf{H}_{P_x+1}(f)&\cdots&\mathbf{H}_{N_x}(f) \end{bmatrix},$$ (2)

where $P_x=\lfloor N_x/2\rfloor+1$, $Q_x=N_x-P_x+1$. $N_x$ denotes the number of traces in the $x$ direction.

The rank-reduction can be achieved by minimizing the following objective function

$$\begin{align*}\begin{split} \text{min} \quad &O=\Vert\mathbf{B}(f)-\mathbf{L}(f)... ...\ \text{subject to} \quad&\text{rank}(\mathbf{L}(f))=L, \end{split}\end{align*}$$ (3)

where $\mathbf{L}(f)$ denotes the low-rank section corresponding to $\mathbf{B}(f)$. $L$ denotes the optimal rank parameter. The above objective function can be conveniently minimized by the well-known singular value decomposition (SVD). Therefore, the SVD of block Hankel matrix $\mathbf{B}(f)$ can be expressed as

$$\mathbf{B}(f)=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T,$$ (4)

where $\mathbf{U}=[\mathbf{u}_1,\mathbf{u}_2,\cdots,\mathbf{u}_N]$ is a left singular matrix, $\mathbf{\Sigma}=[{\sigma}_1,{\sigma}_2,\cdots,{\sigma}_N]$ is a singular matrix and $\mathbf{V}=[\mathbf{v}_1,\mathbf{v}_2,\cdots,\mathbf{v}_N]$ is a right singular matrix. $\sigma_i$ represents the $i$th singular value, and $N$ represents the number of columns of $\mathbf{B}(f)$. Next, the rank-reduced matrix $\mathbf{L}(f)$ can be calculated by extracting the first $L$ principal components, i.e.,

$$\mathbf{L}(f)=\sum^L \limits_{i=1} {\sigma}_i\mathbf{u}_i\mathbf{v}_i^T.$$ (5)

Then, the solved low-rank Hankel matrix is rearranged to a vector form according to inverse Hankelization process. The rank-reduction method can be summarized as follows:

$$\begin{split} \mathbf{R} &=\mathcal{F}^{-1}\mathbf{T}\mathcal... ...\\ \mathbf{T} &=\mathcal{A}\mathcal{L}\mathcal{H}, \end{split}$$ (6)

where $\mathbf{R}$ denotes the separated reflections. The separated diffractions can directly obtained by $\mathbf{D}-\mathbf{R}$. $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote forward and inverse Fourier transforms.

In equation 6, filter $\mathbf{T}$ is referred to as the frequency-domain rank-reduction operator, containing the Hankelization $\mathcal{H}$, principal component extraction $\mathcal{L}$ and inverse Hankelization process $\mathcal{A}$. Because of the complexity of seismic data, the plane-wave assumption is only valid locally. Thus, it is more appropriate to apply the rank-reduction method locally, i.e., in local windows:

$$\mathbf{R} =\mathcal{W}^{-1}\mathcal{F}^{-1}\mathbf{T}\mathcal{F}\mathcal{W}\mathbf{D},$$ (7)

where $\mathcal{W}$ and $\mathcal{W}^{-1}$ denote a pair of windowing and reconstruction operators.