Compressive Sensing-Based Spectral Inversion Algorithm and Its Application in Detailed Fault Characterization

1. Introduction

With the continuous advancement of mechanization and automation in coal mining, the industry has enjoyed remarkable improvements in productivity and efficiency. However, as mining operations reach greater depths, the geological and coal seam conditions become increasingly complex and challenging to manage [1,2,3,4,5]. This added complexity introduces a host of hidden geological factors that, if not properly addressed, can rapidly escalate into severe disasters [6,7,8]. Faults are among the most prevalent and impactful geological factors in coal mining, significantly affecting production efficiency and posing hidden disaster risks. Fault zones exhibit pronounced tectonic stress, extensive coal seam fractures, poor rock integrity, and low strength. These conditions severely impact the deployment of mining faces and mechanized extraction. Additionally, faults often serve as direct channels for water inrush disasters, posing substantial threats to mining safety [9]. Effective management strategies, including advanced geophysical detection and risk mitigation techniques, are essential to address these challenges and ensure safe, efficient mining operations.

Currently, 3D seismic exploration technology has established a comprehensive technical framework for studying disaster risks in coal mining. It has become the optimal method for identifying significant geological features such as major faults, collapse columns, mined-out areas, and other problematic geological structures [10,11,12,13,14]. This advanced technology provides critical insights essential for ensuring the safety and efficiency of mining operations. However, it is important to note that while 3D seismic exploration offers superior lateral resolution compared to conventional drilling methods, its vertical resolution is considerably lower. This limitation poses a challenge, as the technology's ability to accurately detect and resolve smaller faults and intricate geological structures remains insufficient. Therefore, there is an urgent need for high-resolution processing methods to improve the vertical resolution of seismic data, thereby enhancing the geological support capabilities of 3D seismic exploration for safe coal mining operations. In recent years, scholars have recognized this need and proposed various high-resolution processing techniques for seismic data. These methodologies include deconvolution [15,16,17], time-frequency analysis [18,19,20,21,22,23,24,25,26,27], inverse Q-filtering [28,29,30,31,32], spectral whitening [33,34,35], broadband constrained inversion [36,37,38,39,40,41], among others. Each of these techniques has demonstrated varying degrees of success in practical applications, offering promising prospects for the advancement of seismic exploration technology.

Yilmaz [15] introduced the class of deconvolution methods, which is designed to enhance the resolution of seismic data by compressing the seismic wavelet. Since then, a variety of practical algorithms have been developed to achieve this goal. Among these algorithms are pulse deconvolution, predictive deconvolution, homomorphic deconvolution, and surface-consistent deconvolution. Although these deconvolution methods have proven to be effective in enhancing the resolution of seismic data, they come with certain limitations and challenges. One of the primary issues is that these methods are based on several assumptions that may not always hold true in real-world scenarios. For instance, they often assume that the seismic wavelet is minimum phase or that the subsurface is horizontally layered, which may not be accurate. Moreover, the application of deconvolution can lead to a reduction in the signal-to-noise ratio (SNR).

Time-frequency analysis methods have garnered significant attention in geophysics, primarily for enhancing the resolution of seismic data. Huang [24] introduced a high-resolution processing technique for seismic data utilizing the generalized S-transform. Mao [22] proposed an adaptive continuous wavelet domain method tailored specifically for high-resolution seismic data processing. They provided detailed guidelines for selecting reference frequencies and calculating the frequency extension range, which enable effective enhancement of seismic data resolution while maintaining a high signal-to-noise ratio.

The application of inverse Q filtering techniques in seismic data processing is aimed at mitigating the attenuation effects caused by subsurface materials. By compensating for amplitude decay and frequency loss, this method effectively improves the phase characteristics and continuity of seismic records. Currently, three main types of inverse Q filtering methods are employed both internationally and domestically: series expansion-based approximate high-frequency compensation methods, wavefield extrapolation-based inverse Q filtering methods, and other inverse Q filtering techniques [32]. Accurate estimation of the subsurface quality factor Q is crucial for the successful application of inverse Q filtering methods. However, further research is needed to overcome the challenges of accurately and efficiently applying inverse Q filtering in complex underground structures.

The spectral whitening method aims to flatten the high-frequency portion of seismic data by adjusting the amplitude spectral distribution to compensate for missing components. When employing this method, the phase spectrum of seismic data remains unaltered, with only the high-frequency amplitude spectrum being elevated accordingly. This approach is known for its simplicity and robustness. However, it is worth noting that it may have the drawback of relatively poor fidelity.

Spectral inversion technology is widely recognized as an effective method for enhancing resolution in wideband constrained inversion techniques. Puryear and Castagna [36] proposed the reflection coefficient parity decomposition theory, which decomposes the reflection coefficient into odd and even components. By leveraging the distinct resolution capabilities of these components for thin layers, spectral inversion target functions yield a high-resolution reflection coefficient. This theory has significantly advanced the application of spectral inversion in enhancing the resolution of seismic data processing. However, conventional spectral inversion methods usually solve the objective function under $L_2$ norm constraints, resulting in limited effectiveness for resolution enhancement.

Spectral inversion relies on the crucial assumption that reflection coefficients are sparse. To address this, we introduce compressed sensing theory [42,43,44,45] into the spectral inversion. Spectral inversion is typically an underdetermined optimization problem, where the number of unknowns exceeds the number of equations. In this study, we employ the $L_0$ norm for sparse constraints, given its effectiveness in promoting sparsity. We propose a novel sparsity-adaptive hard thresholding pursuit algorithm to solve the sparse-constrained optimization problem, aiming to accurately estimate the reflection coefficients. It does not require the use of the Matching Pursuit algorithm to obtain initial sparse reflection coefficients [46]. Furthermore, we present a comprehensive method to enhance the resolution of stabilized seismic data. First, the sparsity-adaptive hard thresholding pursuit algorithm is used to solve the $L_0$ norm-constrained optimization problem, yielding the reflection coefficients. Following this, we perform wavelet extraction based on shaping regularization to ensure that the extracted wavelets are stable and consistent, followed by wavelet decomposition and high-frequency wavelet replacement, ultimately achieving the objective of improving seismic data resolution. Applying this processing workflow to actual coal mine seismic data significantly enhances data resolution. This improvement aids in the identification of small-scale fault systems, which are often difficult to detect with conventional methods. Consequently, this enhanced capability provides better geological support, greatly contributing to the safe production and operation of coal mines.

2. Methods

2.1. Spectral Inversion Objective Function

Based on the convolution model, seismic records can be represented as the convolution of reflection coefficients and a wavelet:

$$d(t)=w(t)\ast r(t)+n(t)$$

(1)

where $n(t)$ is noise. By applying the Fourier transform to the above equation, we get:

$$D(f)=W(f)\ast R(f)+N(f)={\displaystyle \sum _{n=1}^N{A}_n\cdot W(f)e^{-i\cdot 2\pi f{\tau }_n}}+N(f)$$

(2)

Here, $A_n$ represents the amplitude of the reflection coefficients. For a specific frequency component $f_m(k=1,\dots ,M)$, it can be written as:

$$D(f_m)={\displaystyle \sum _{n=1}^N{A}_n\cdot W(f_m)e^{-i\cdot 2\pi f_m{\tau }_n}}+N(f)$$

(3)

Letting $Q=e^{-i\cdot 2\pi f_m{\tau }_n}$, the equation can be rewritten in matrix form as:

$$d=WQx+N=\Phi x+N$$

(4)

Where $\Phi =WQ$, $d$ is the observed seismic data, $W$ is the diagonal matrix of seismic wavelets, $Q$ is the Fourier basis matrix, and $x$ is the sequence of reflection coefficient amplitudes. Applying the $L_0$ norm for sparse constraint, the optimization problem can be formulated as:

$$\widehat{x}=\arg \min {\Vert x\Vert }_0,s.t.\Phi x=d$$

(5)

The $L_0$ norm is defined as the number of non-zero elements in a vector. By incorporating the $L_0$ norm sparse constraint, the spectral inversion method can obtain reflection coefficients that better approximate actual conditions.

2.2. High-Resolution Seismic Data Processing Workflow

This paper introduces a robust approach for improving seismic data resolution using the spectral inversion algorithm. The proposed method consists of two crucial steps:

(1) Employ the sparsity adaptive hard thresholding pursuit algorithm to solve the norm-constrained optimization problem and obtain accurate reflection coefficients.

(2) Apply shaping regularization for wavelet extraction, followed by wavelet decomposition and high-frequency wavelet replacement, thereby effectively enhancing the overall resolution.

2.2.1. Spectral Inversion Algorithm Based on Compressed Sensing Theory

As shown in the equation (5), we solved the optimization problem to obtain the reflection coefficients.

While threshold algorithms [47,48] provide fast computational results for solving sparse constrained problems, they often yield subpar outcomes compared to matching pursuit class algorithms [49,50]. The latter has gained significant attention due to its superior practical application results. The orthogonal matching pursuit algorithm, for instance, selects only one atom in each iteration that best matches the residuals. However, this selection process becomes inefficient when dealing with high-dimensional signals. Building upon the IHT algorithm, Foucart [51] introduced the Hard Thresholding Pursuit (HTP) algorithm. In each iteration, HTP selects the atoms corresponding to the largest K elements in the sparsity coefficients of the signal. The sparsity coefficients of the signal are then obtained by solving a minimization problem. The iterative process of the algorithm can be described as follows:

$$s^{n+1}=\{\mathrm{indices}\mathrm{of}\mathrm{K}\mathrm{largest}\mathrm{entries}\mathrm{of}{x}^n+{\Phi }^T(d-\Phi x^n)\}$$

(6)

$$x^{n+1}=argmin\{{\Vert d-\Phi z\Vert }_2,supp(z)\subseteq S^{n+1}\}$$

(7)

where $supp(z)$ express the support set of $z$. During the second step, the minimization problem can be solved using either the least squares algorithm or the gradient descent algorithm. However, when the signal's dimensionality is large or when it contains a significant number of non-zero elements, the least squares algorithm becomes computationally challenging due to the calculation of a high-dimensional inverse matrix. To address this issue, we employ the Nesterov accelerated gradient (NAG) algorithm, known for its fast convergence and high efficiency, to solve the minimization problem in Equation (5). The process of reconstructing sparse signals using the HTP algorithm can be summarized as follows:

Input: sensing matrix $\Phi \in R^{M\times N}(M\le N)$ , measurement vector $d\in R^{M\times 1}$ , sparse level $K$;

Initialization: sparse coefficient approximation $x\in 0$ , index set ${\Lambda }_0=\varnothing$ , atomic matrix ${\Phi }_0=\varnothing$;

Loop:

(1) calculate ${\Upsilon }^n=x^n+{\Phi }^T(d-\Phi x^n)$(2) remain the index ${\Lambda }_n$ of the maximum K elements in ${\Upsilon }^n$, and the atomic matrix in this index is ${\Phi }_n$;

(3) solve $\min {\Vert d-\Phi x\Vert }_2^2,s.t.{\Vert x\Vert }_0\le K$ by the Nesterov accelerate algorithm:

Input: learning rate $\eta$, attenuation rate $\epsilon$Initialization: initial velocity $v_0=0$, objective function $J(x)=\min {\Vert d-\Phi x\Vert }_2^2,s.t.{\Vert x\Vert }_0\le K$

Loop:

$$① \mathrm{Calculate}{v}_t=\epsilon v_{t-1}+\eta {\nabla }_{x}J(x-\epsilon v_{t-1})$$

$$②\mathrm{Update}{x}_t=x_{t-1}-v_t$$

If the condition is met, stop iteration; else , return step① continue iteration.

(4) if the condition is met, stop iteration; else n = n +1, return step①continue iteration.

In addition, the sparsity of reflection coefficients in the transform domain is often unknown, and the accuracy of capturing this sparsity directly impacts the reconstruction accuracy. Therefore, it is crucial to develop a reconstruction algorithm that can automatically adapt to the signal's sparsity. In this paper, we propose a sparsity-adaptive hard threshold pursuit algorithm. The algorithm consists of the following steps:

Input：sensing matrix $\Phi \in R^{M\times N}(M\le N)$, measurement vector $d\in R^{M\times 1}$, noise level $\epsilon$ ;

Initialization: $K^1=0$， $K^2=\max (0.001N,1)$，$x^1=0$, $F(x^1)=\Vert \Phi x^1-d\Vert /\Vert d\Vert =1$, solve $x^2$ in the condition of $\Vert x_0\Vert \le K^2$, $F(x^2)=\Vert \Phi x^2-d\Vert /\Vert d\Vert$, $t=1$;

Loop:

(1) $K^{t+2}=K^{t+1}+(\epsilon -F(x^{t+1}))(K^{t+1}-K^t)/(F(x^{t+1})-F(x^t))$(2) solve $x^{t+2}$ by the Hard Thresholding pursuit algorithm

(3) calculate $F(x^{t+2})=\Vert \Phi x^{t+2}-d\Vert /\Vert d\Vert$

(4) if the condition is met, stop iteration; else, $t=t+1$, return step (1) continue iteration;

Output: sparse level $K$ , sparse coefficient $x$ .

To compare the effects of solving the objective function under different norm constraints, we replace the $L_0$ norm constraint in equation (5) with an $L_2$ norm constraint of the optimization problem. We present the inversion results under the two norm constraints using a seismic data example, as shown in Figure 1. The comparison results reveals that the $L_0$ norm results exhibit a higher level of sparsity and signal-to-noise ratio compared to the $L_2$ norm constraint. This enables a broader bandwidth seismic reflection coefficient model.

2.2.2. Continuous Wavelet Transform

The Morlet wavelet $\psi (t)$ can be viewed as a complex sinusoid modulated by a Gaussian envelope (Morlet J, 1982):

$$\psi (t)=e^{-\frac{t^2}{2}}{e}^{j\sigma t+j\phi }$$

(8)

Here, $\sigma$represents the modulation frequency of the Morlet wavelet, which typically needs to be greater than 5.33. The Fourier transform of the Morlet wavelet is:

$${\widehat{\psi }}_A(\omega )=\sqrt{2\pi }{e}^{j\phi }{e}^{-\frac{{(\omega -\sigma )}^2}{2}}$$

(9)

Subsequently, the optimal scale of the Morlet wavelet is obtained by optimizing the Fourier transform of the seismic trace $D(\omega )$ and the Fourier transform of the reflection coefficient $R(\omega )$ . The optimal scale $a_{opt}$ can be determined by minimizing the objective function in equation (10):

$$\arg \min {\left|\left|D(\omega )\right|-\left|R(\omega ){\widehat{\psi }}_A(a\omega )\right|\right|}^2$$

(10)

Using the wavelet $\psi \left(\frac{t}{a_{opt}}\right)$ with the optimal scale $a_{opt}$ and the positions of the reflection coefficients, decompose the seismic data to obtain the coefficients of each reflected wavelet. Assuming the positions of the reflection coefficients are ${\tau }_1,{\tau }_2,\cdots ,{\tau }_q$, the coefficients ${\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_q$ of each reflected wavelet can be obtained by minimizing the objective function in equation (11):

$${\Vert d-{\displaystyle \sum _{m=1}^q{\alpha }_m\psi \left(\frac{t-{\tau }_m}{a_{opt}}\right)}\Vert }^2$$

(11)

Finally, replace and reconstruct $\psi (t)$ using high-frequency Morlet wavelets or other high-frequency wavelets, which can naturally extending low- and high-frequency energy of seismic data.

$$d_{new}={\displaystyle \sum _{m=1}^q{\alpha }_m{\psi }_{high}\left(t-{\tau }_m\right)}$$

(12)

Applying the aforementioned steps to seismic data processing can significantly improve the resolution of the seismic profile.

Building on the outlined key steps, this paper introduces a robust resolution enhancement technique for spectral inversion by leveraging compressive sensing theory (shown in Figure 2). To ensure stability during the wavelet calculation process, a shaping regularization wavelet extraction approach is employed. This technique enhances the estimation stability of wavelets, particularly in regions exhibiting low reflection coefficients or a low signal-to-noise ratio in seismic data.

3. Results

The effectiveness of the aforementioned method was validated by applying it to the processing of three-dimensional seismic data in a coal mine. The exploration area is currently in the design stage, and in order to ensure the smooth development of the exploration area and normal production after the mining area goes into operation, maximize the advantages of mechanized mining and ensure production safety, therefore, three-dimensional seismic exploration was implemented in the exploration area to gain early insights into underground structures and geological anomalies and to take necessary precautionary measures in advance. It was confirmed through drilling and actual mining operations that the 12th coal seam developed in the area is approximately 2 meters thick. Reflection character analysis revealed that this coal seam (T12) exhibits a relatively simple phase and strong energy, allowing for continuous tracking throughout the entire region. However, the coal seam is affected by small faults, which pose risks to mining safety. Therefore, it is necessary to characterize the distribution of fault planes in the T12 coal seam in this area. High-resolution processing of the seismic data was performed, and the obtained results are shown in Figure 3. It can be observed that this technique can recover thinner layer information (circled in yellow and blue) and provide clearer fault identification(circled in red) while preserving low frequencies effectively. It also improves the signal-to-noise ratio and enhances the temporal and spatial resolution of the data. The processing yields high-fidelity, high-resolution, and high signal-to-noise ratio results.

Further validation of the reliability and effectiveness of high-resolution seismic data processing was conducted by extracting attributes along layers and comparing them with those from the original data. From the amplitude attribute map along the T12 horizon of the survey area, it can be seen that the high-resolution processing results provide a more detailed depiction of formations and faults compared to the original seismic data (annotated with yellow dashed boxes), while also revealing more geological details, as shown in Figure 4. This indicates that, after enhancing resolution, broadband seismic data can better delineate stratigraphic information and small-scale structures.

Due to the extensive development and regional distribution of small faults within the coal seam, areas where the integrity of the coal seam is compromised exhibit weakened reflected energy and reduced relative impedance. Therefore, impedance attribute maps (shown in Figure 5) along T12 (annotated with red line) have also been created. The impedance attributes along the original seismic data can reflect the main area of small fault distribution, while the impedance data processed with improved resolution can display the impedance reduction phenomenon caused by the development of small faults with greater precision. The process of enhancing seismic data resolution allows for a more precise visualization of the subtle variations in impedance caused by the presence of small faults (annotated with black dashed boxes). Furthermore, this expanded resolution enables us to capture a wider range of small fault distribution, resulting in a more comprehensive and accurate depiction of their planar layout.

4. Discussion and Conclusions

The mechanization and intelligence of coal mining are continuously improving, along with increasing mining depths. Complex geological conditions and coal seam characteristics pose higher demands for safety in coal mine production. While 3D seismic exploration has evolved into a comprehensive technical system capable of detecting geological structures, spatial positioning, and lithology interpretation, its ability to identify small-scale structures remains limited. To address this limitation, this paper proposes a technical workflow based on compressed sensing spectral inversion, aimed at enhancing the resolution of seismic data. By effectively improving the resolution while preserving amplitude and fidelity, this method allows for a more precise characterization of small-scale geological information embedded within the seismic data. In order to validate the efficacy of this approach, the method was applied to real 3D seismic data obtained from a coal mining, resulting in enhanced data resolution and a more detailed depiction of stratigraphic information and faults. This enhancement provides robust support for subsequent safe production measures.