## Appendix A Post-Stack Seismic Inversion Technology  
In the petroleum industry, post-stack seismic inversion technology is a conventional and mature technology to characterize reservoirs. Ignoring multiples, transmission loss, spherical spreading, and attenuation, the objective function of post-stack seismic inversion can be expressed as  
\[\underset{\mathbf{z}}{\mathbf{min}}\|\mathbf{d}-\mathbf{W}\mathbf{r}(\mathbf{ z})\|_{2}^{2}, \tag{1}\]  
where \(\mathbf{z}=\left[z_{1},z_{2},\ldots,z_{g}\right]^{T}\) is the acoustic impedance vector, subscript \(g\) is the sampling number, superscript \(T\) is the transpose operator, \(\mathbf{d}=\left[d_{1},d_{2},\ldots,d_{g+q-2}\right]^{T}\) is the real seismic record, superscript \(q\) is the length of seismic wavelet vector \(\mathbf{w}=\left[w_{1},w_{2},\ldots,w_{q}\right]^{T}\), and \(\mathbf{W}\) is the seismic wavelet convolution matrix  
\[\mathbf{W}=\left[\begin{array}{cccc}w_{1}&&&&\\ w_{2}&w_{1}&&\\ \vdots&w_{2}&\ddots&\\ w_{q}&\vdots&\ddots&w_{1}\\ &w_{q}&\ddots&w_{2}\\ &&\ddots&\vdots\\ &&&w_{q}\end{array}\right]_{(q+g-2)\times(g-1)}, \tag{2}\]  
and \(\mathbf{r}(\mathbf{z})\) is  
the Longmaxi Formation, the number of time sampling points during the inversion process is about 50 (Figures 10a and 11a). Due to these properties, we approximate density of Longmaxi Formation to be a constant and then the Equation (A3) can be can be converted to be  
\[\mathbf{r}(\mathbf{v})=\left[r_{1},r_{2},\dots,r_{g-1}\right]^{T}=\left[\frac {v_{2}-v_{1}}{v_{2}+v_{1}},\frac{v_{3}-v_{2}}{v_{3}+v_{2}},\dots,\frac{v_{g}-v _{g-1}}{v_{g}+v_{g-1}}\right]^{T}.\] (A4)  
where \(\mathbf{v}=\left[v_{1},v_{2},\dots,v_{g}\right]^{T}\) is the velocity vector.  
Due to limitation in the frequency bandwidth of seismic data, directly solving for the objective function (13) is strongly unstable. Based on the Tikhonov regularization theory and the Bayes theory, a stable and reliable impedance inversion result can be obtained [30]. The impedance update during iterative process of seismic inversion is expressed as  
\[\mathbf{v}^{(y+1)}=\mathbf{v}^{(y)}+\left[\left(\mathbf{G}^{(y)}\right)^{T} \mathbf{G}^{(y)}+\lambda^{(y)}\right]^{-1}\left(\mathbf{G}^{(y)}\right)^{T} \left[\mathbf{d}-\mathbf{W}\mathbf{r}\left(\mathbf{v}^{(y)}\right)\right],\] (A5)  
where superscript \(y=1\), \(2\),..., \(Y\) is the iteration number, \(\mathbf{v}^{(y)}\) is the velocity at the \(y\)th iteration, \(\mathbf{v}^{(1)}\) is the initial velocity model from sonic logging data, \(\mathbf{G}^{(y)}\) is the Jacobi matrix at the \(y\)th iteration  
\[\mathbf{G}^{(y)}=\mathbf{W}\left[\begin{array}{cccc}\frac{-2v_{2}^{(y)}}{ \left(v_{2}^{(y)}+v_{1}^{(y)}\right)^{2}}&\frac{v_{1}^{(y)}}{\left(v_{2}^{(y)} +v_{1}^{(y)}\right)^{2}}&&&\\ &\frac{-2v_{2}^{(y)}}{\left(v_{3}^{(y)}+v_{2}^{(y)}\right)^{2}}&\frac{v_{2}^{( y)}}{\left(v_{3}^{(y)}+v_{2}^{(y)}\right)^{2}}&&&\\ &&&\ddots&\\ &&&\frac{-2v_{Z}^{(y)}}{\left(v_{g}^{(y)}+v_{g-1}^{(y)}\right)^{2}}&\frac{v_{g- 1}^{(y)}}{\left(v_{g}^{(y)}+v_{g-1}^{(y)}\right)^{2}}\end{array}\right]\] (A6)  
\(\lambda^{(y)}=\chi_{\Delta t}^{(y)}/\chi_{\Delta t}^{(y)}\) is the adaptive regularization parameter at the \(y\)th iteration, \(\chi_{\Delta t}^{(y)}\) is the variance of the difference between the real and synthetic seismograms at the \(y\)th iteration, \(\chi_{\Delta t}^{(y)}\) is the variance of the velocity residual before and after the \(y\)th iteration. The expression of regularization parameter implies that when the noise level in the seismic data is high, the value of regularization parameter is large, thus serving to suppress noise to obtain a relatively stable inversion result. On the contrary, when the noise approaches zero, the value of the regularization parameter is close to zero, and the inversion result is then determined mainly from the seismic data, which leads to a high-resolution solution.