## 2 Methodology
### Principle of Reflectivity Sparse Inversion  
Reflectivity is a parameter that represents the difference in the physical properties of two rock layers. It is presented in the form of a pulse sequence and has a wide-band frequency spectrum. The seismic data collected are influenced by Earth filtering, and the seismic wavelet is band-limited. Reflectivity inversion eliminates the effects of stratigraphic filtering and estimates the true spatial structure of the subsurface layers from seismic data, aiming to enhance the resolution of seismic data. Assuming that the seismic wavelet is stationary within a certain spatiotemporal range, a convolution model can be used in the time domain to represent it:  
\[s(t)=w(t)*r(t)+n(t), \tag{1}\]  
where \(*\) denotes convolution operation, \(s(t)\) denotes the seismic record, \(w(t)\) represents the seismic wavelet, \(r(t)\) represents the impulse response, commonly referred to as the reflectivity, and \(n(t)\) denotes the noise. The equation can be rewritten in matrix form, involving matrix multiplication:  
\[\mathbf{S}=\mathbf{W}\mathbf{R}+\mathbf{N}, \tag{2}\]  
where \(\mathbf{S}\) represents the seismic record, \(\mathbf{W}\) represents the wavelet matrix, \(\mathbf{R}\) represents the reflectivity, and \(\mathbf{N}\) represents the noise. The wavelet matrix is formed by arranging seismic wavelets and reflects the inversion and time-shift processes in the convolution model Equation (1). In seismic exploration, inversion is generally expressed as follows:  
\[\mathbf{R}=\mathbf{G}\mathbf{S}, \tag{3}\]  
where \(\mathbf{G}\) is the inversion operator. For the noisy convolution model, there is an objective function:  
\[\mathbf{J}=\|\mathbf{S}-\mathbf{W}\mathbf{R}\|_{2}^{2}+\lambda\|\mathbf{R}\|_{p}^ {p}\, \tag{4}\]  
where \(\lambda\) is the regularization parameter, which balances the weight between the data fitting term and the sparsity constraint term. \(p\) is the norm order of the sparsity constraint, typically set as \(0<p\leq 1\), and controls the sparsity of the reflectivity. The symbol \(\|\cdot\|_{p}^{p}\) denotes the Lp-norm, which measures the sparsity of the reflectivity vector. By adjusting these parameters, the inversion results can be optimized to better align with the prior knowledge of geological structures. In the case of noisy data, based on the least squares method, the inversion reflectivity \(\mathbf{R}\) can be expressed as follows:  
\[\mathbf{R}=(\mathbf{W}\mathbf{W}^{\intercal}+\lambda\mathbf{Q})^{-1}\mathbf{ W}^{\intercal}\mathbf{S}, \tag{5}\]where \(\mathbf{Q}\) is a sparse matrix. According to Equation (3), the least squares inversion operator \(\mathbf{G}\) for data can be expressed as follows:  
\[\mathbf{G}=(\mathbf{W}\mathbf{W}^{\intercal}+\lambda\mathbf{Q})^{-1}\mathbf{W}^{ \intercal}, \tag{6}\]  
Based on different solution algorithms, both the inversion operator \(\mathbf{G}\) and the inversion reflectivity \(\mathbf{R}\) will be iteratively updated continuously. In this paper, the alternating direction method of multipliers (ADMM) solution algorithm [32] is taken as an example. We introduce an auxiliary variable \(\mathbf{Z}\) and rewrite the objective function as follows:  
\[\mathbf{min}\parallel\mathbf{S}-\mathbf{W}\mathbf{R}\parallel_{2}^{2}+ \boldsymbol{\lambda}\parallel\mathbf{Z}\parallel_{\mathbf{p}}\mathbf{s.t.R}= \mathbf{Z}, \tag{7}\]  
In this context, we formulate the Lagrangian equations to solve the aforementioned system of equations. The augmented Lagrangian function is defined as follows:  
\[\mathbf{L}(\mathbf{R},\mathbf{Z},\mathbf{U})=\parallel\mathbf{S}-\mathbf{W} \mathbf{R}\parallel_{2}^{2}+\boldsymbol{\lambda}\parallel\mathbf{Z} \parallel_{\mathbf{p}}+2\mathbf{p}\parallel\mathbf{R}-\mathbf{Z}+\mathbf{U} \parallel_{2}^{2}, \tag{8}\]  
where \(\mathbf{U}\) is the Lagrange multiplier, and \(\mathbf{\rho}>0\) is the penalty parameter.