According to the seismic convolution model, the seismic signal is generated by the convolution of the seismic reflection coefficient with the seismic wavelet in the time domain. In the frequency domain, the seismic signal spectrum is calculated by multiplying the amplitude spectrum and the phase spectrum of the seismic signal as follows: s = w ⊗ r Time ⇔ S = A ⋅ e − i ⋅ A n g Frequency (1) where the left side of the equivalent is defined as Eqs 1-1, and w is seismic wavelet. In this paper, it is assumed that the seismic wavelets are zero-phase, and r is a seismic reflection coefficient, which is sparse. s is a seismic record. The right side of the equivalent is defined as Eqs 1-2, and S is the spectrum of the seismic record. A is the amplitude spectrum of the seismic record. e − i * A n g is the phase spectrum of the seismic record and, A n g is the phase.

According to Eq. 1, the deconvolution procedure entails the compression of the seismic wavelet, leading to an expansion of the seismic frequency band in the frequency domain. Therefore, determining the range of frequency broadening is very important. The technique of seismic record spectrum scanning can help determine the effective frequency bandwidth. The limited band information is 4–64 Hz, as shown in Figure 1. In practical applications, it is advisable to use a narrower scanning interval band. Then, the distribution of effective information within the band can be more accurately determined.

According to the results of the spectrum scanning, the spectrum fitting method is employed to expand the frequency range within the effective frequency band. Then, the octave can be increased to enhance seismic resolution. In this paper, the spectral fitting method adopts a frequency division weighted superposition approach. At first, Gaussian functions are constructed in different frequency bands. These functions are subsequently employed for frequency division processing, as depicted in Figure 2A. The Gaussian functions of different frequency bands are expressed as f i , where i is the frequency and its number is the total number in the effective frequency bandwidth. The size of f i is typically approximately 10 Hz in bandwidth and its number is the same as that of i . Figure 2B shows the amplitude spectrum of a real seismic record.

Gaussian functions with different frequencies are employed to limit the amplitude spectrum of the seismic record, leading to the acquisition of seismic frequency division amplitude spectra (illustrated by the colored line in Figure 3A). According to the results of spectrum scanning, the desired range of wavelet amplitude spectrum can be determined accordingly. The amplitude spectra of the frequency division are weighted and superimposed, and the expansion of the spectrum is performed while adhering to the constraints of the desired wavelet amplitude spectrum, as below: F r e q α = argmin α i ∑ i A d i a g f i α i − W 2 2 = argmin α A f α − W 2 2 , (2) where α is the weight coefficient of the frequency division amplitude spectrum, which is a vector. α i i = 1 , 2 , ⋯ , N are the elements of α . W is the expected wavelet spectrum, which is determined by the effective frequency bandwidth, as shown by the blue dashed line in Figure 3A. A f = A d i a g f i i = 1 , 2 , ⋯ , N is the matrix constructed by different frequency division amplitude spectra. Within the confines of the expected wavelet amplitude spectrum, the amplitude spectra of various frequency bands are amalgamated through the utilization of weight coefficients. This effectively expands the frequency range of the seismic record. Simultaneously, this method effectively maintains the attributes of the seismic spectrum curve, as shown in Figure 3B. In short, the spectral fitting method maintains the seismic spectrum pattern by superimposing and dividing frequencies, while also upholding the desired wavelet constraint.

Following the acquisition of high-resolution seismic data and wavelets through the spectral fitting method, it is possible to conduct sparse deconvolution. The conventional method for sparse deconvolution involves constructing a sparse objective function when the wavelet is known. This can be done as follows: L r = W r − s 2 2 + λ r 1 (3) where W is the seismic wavelet matrix. r is the seismic reflection coefficient, s is the seismic record, and λ is the scale coefficient. Sparse deconvolution requires an accurate seismic wavelet. Figures 4A,B show the results of sparse deconvolution. Evidently, the sparse deconvolution method yields superior spectral recovery within the effective band. However, there is a significant discrepancy in the spectrum between 100 and 200 Hz outside the frequency range (as depicted in Figure 4B), resulting in a substantial difference between the calculated reflection coefficient and the true seismic reflection coefficient (as illustrated in Figure 4A). The primary factor is that the seismic spectral energy in high frequencies is diminished, leading to reduced accuracy in the recovery of this portion through sparse deconvolution.

In order to achieve a precise sparse reflection coefficient, it is essential to capture accurate spectral properties of high-frequency signals. According to Eq. 2 and Eq. 1, it is easy to establish a relationship between the seismic spectrum and the seismic reflection coefficient. s = W ⋅ r ⇔ F − 1 A f α e − i ⋅ A n g (4)

Where F − 1 is the Fourier inverse matrix. According to Eq. 2, W is the wavelet matrix, and the wavelet spectrum can be accurately set. So, W is deterministic. Thus, Eq. 4 can be derived as a formula for the reflection coefficient. r = W − 1 F − 1 A f α e − i ⋅ A n g (5)

Where W − 1 is the inverse matrix of W . According to Eqs 3, 5, a new form of sparse deconvolution objective function can be constructed. F r e q _ L 1 α = A f α − W 2 2 + λ W − 1 F − 1 A f e − i ⋅ A n g α 1 (6)

After obtaining the weight coefficient α , the seismic reflection coefficient can be obtained according to Eq. 5.

The recently developed objective function offers two benefits in comparison to conventional sparse deconvolution techniques. The initial point to consider is the elimination of the necessity to extract the wavelet, thereby mitigating potential errors in the extraction of the seismic wavelet. The second point is that the expected wavelet encompasses a wider spectrum of frequencies, resulting in a more accurate restoration of spectral characteristics. Compared with Figures 4A,C, the two sparse deconvolution methods show significant differences in the reflection coefficient. Figure 4C depicts the outcome of the application of the proposed method, demonstrating a significantly improved accuracy in the estimation of the inverted reflection coefficient. The primary factor is that the conventional sparse deconvolution method cannot effectively handle the spectral characteristics of high-frequency seismic records (green curve in Figure 4B). This results in a notable disparity between the calculated reflection coefficient and the true value. The sparse deconvolution method proposed in this paper is closer to the true seismic reflection coefficient because it can more effectively restore the spectral characteristics of high-frequency seismic signals (as depicted by the green curve in Figure 4D). This demonstrates the effectiveness of the objective function in this paper.

There are differences in the reflection coefficients obtained by the sparse deconvolution method in different seismic channels, encompassing discrepancies in time bias and amplitude fluctuations. These differences contradict the expected gradual lateral transition of geological features. So, the spatial continuity constraints are essential for achieving sparse solving while maintaining constructive control.

The first-order differential matrix, commonly referred to as total variation (TV) regularization, lead to a step-like effect that is not appropriate for spatial control constraints. The second-order differential matrix possesses smooth characteristics, enabling it to preserve the surface of three-dimensional data and ensure the spatial continuity of seismic data. In this paper, the method of regularizing the Hessian matrix is used to impose spatial continuity constraints. The Hessian matrix is composed of multiple second-order differential matrices. Assuming the three-dimensional seismic data is s , then its Hessian matrix is: H = L x , x s L x , y s L x , t s L y , x s L y , y s L y , t s L t , x s L t , y s L t , t s (7) where H is the Hessian matrix. L i , j s : = ∂ s ∂ i ∂ j = D i j ⊛ s , D i j is the differential filter. x , y , t are the three different directions of the 3D seismic data—the direction of the survey line, the survey trace, and the time.

The regularization of the Hessian matrix can be expressed as: U s = ∑ i , j ∈ x , y , t L i , j s 2 2 = ∑ i , j ∈ x , y , t l i , j ⊛ s 2 2 (8) where l i , j are second-order differential operators in different directions (Gholami and Sacchi, 2013). The regularization of the Hessian matrix is computed for the 3D geological model data, as depicted in Figure 5. The figure demonstrates that the Hessian matrix exhibits surface smoothing characteristics. Therefore, it can be employed to limit the spatial coherence of the seismic reflection coefficient.

In order to be able to optimize the solution for seismic data = F − 1 A f α e − i ⋅ A n g , the convolution operation of Hessian in Eq. 8 should be changed to a matrix operation (Guo et al., 2022). U s = ∑ i , j ∈ x , y , t L i , j s 2 2 = ∑ i , j ∈ x , y , t P l i , j ⋅ s 2 2 (9) where P l i , j is the matrix form of the differential operators.

Combine Eqs 6, 9 to get the final objective function: F α = A f α − W 2 2 + λ r 1 + γ ∑ i , j ∈ x , y , t P l i , j ⋅ r 2 2 (10) where r = W − 1 F − 1 A f e − i ⋅ A n g α . λ and γ are the scale factors, which can be preferentially determined by model testing. The objective function (Eq. 10) is solved using the Split Bregman algorithm (Guo et al., 2022). The weight coefficient α can be obtained by solving the objective function. Then, the high-resolution seismic record and the seismic reflection coefficient can be obtained through Eq. 11, respectively. s g = F − 1 A f e − i ⋅ A n g α r g = W − 1 F − 1 A f e − i ⋅ A n g α (11)

Where, s g is the high-resolution seismic record. r g is the seismic reflection coefficient.

For testing, 3D seismic data from the GST area of the Sichuan Basin was utilized. The targeted stratum is the Qixia Formation, and its thickness remains relatively stable at approximately 110 m. The depth of the targeted stratum is 4,500–4,700 m. The reservoir type is a dolomite pore reservoir, characterized by low porosity and low permeability. The physical parameters of the reservoir are very similar to those of the surrounding rock. The reservoir thickness is thin, with each group of reservoirs ranging from 6 to 10 m in thickness. The dominant frequency of the seismic data is 25 Hz. The response of the reservoir within the Qixia Formation is disrupted by the strong reflection from the upper and lower boundaries of the targeted stratum. Low-resolution seismic data in the Sichuan Basin presents two primary challenges: firstly, it results in unclear seismic structural characteristics, and secondly, it obscures thin reservoir responses.

The reservoir of the Qixia Formation in the GST area of the Sichuan Basin is predominantly situated within the central portion of the targeted stratum. The forward geological model is designed to depend on reservoir characteristics, as shown in Figure 9A. The arrow points to the dolomite reservoir, which has a designed thickness of 10 m. The position of the reservoir gradually shifts from the left to the right, moving toward the center of the targeted stratum. The seismic forward data is shown in Figure 9B. The forward simulation process (from the geological model to the seismic record) uses the theory of the Zoeppritz equation. The dominant frequency of the seismic wavelet is 35 Hz. It can be observed that the reservoir exhibits a strong bright spot response when it is located near the center of the targeted stratum. The bright spot response of the reservoir gradually weakens as it approaches the upper boundary of the targeted stratum. The forward model test shows that the presence of a reservoir leads to a bright spot response within the targeted stratum. However, in practice, the bright spot response of thin reservoirs may not be readily apparent due to limitations in the resolution and signal-to-noise ratio of seismic data. High-resolution seismic data is essential for accurately predicting the presence of thin reservoirs.

Figure 10 illustrates a comparison of the efficacy of different methods on real seismic data. Figure 10A shows the original seismic section. There are two industrial gas wells named Well 1 and Well 2, but there is no bright spot response from the reservoir at Well 1. This is primarily due to the scarcity of seismic data, leading to an extended duration of the wavelet. The identification of the reservoir is challenging due to the obstruction caused by the side lobes of the seismic-reflected wavelets at the top and bottom. Simultaneously, the stratum contact (P1m) consists of low-speed argillaceous limestone above and high-speed limestone below. It is observed as a peak response on the seismic section. Nevertheless, the presence of seismic noise and constraints in frequency ranges, combined with the impact of intricate lithological formations in the upper section of the targeted stratum, result in the overlapping of seismic wavelets from different reflection interfaces. Consequently, it is not feasible to precisely track the position of the stratum. Hence, the interpretation of the P1m strata (as shown in Figure 10A, at position ③) presents a challenging task. Two methods are employed to improve the resolution of seismic data in order to address this issue. Figure 10B shows the effects of spectral modeling deconvolution with spatial continuity constraints. It can be seen that the resolution of seismic data is effectively improved, particularly in the vicinity of the reservoir (Figure 10B, at positions ① and ②). The bright spot response of the reservoir is very clear. However, there is a sudden increase or decrease in the energy level at the upper boundary of the targeted stratum (P1m), and the continuity of the stratum deteriorates, as shown at the circled position in Figure 10B (③). The primary factor is that the upper part of the targeted stratum contains a complex combination of lithologies. Furthermore, a flaw exists in the approach employed to retrieve the energy of the reflection interface, leading to an impact on the seismic reflection within the upper portion of the targeted stratum and causing a disruption in the energy convergence of P1m. Figure 10C shows the application effect of the proposed method. The reservoir is also clearly highlighted, and the bright spot response corresponds more accurately to the industrial gas wells. The convergence of the formation interface (P1m) energy and the interpretation of the structure becomes easier. The interpretation results are consistent with the principles of geology and logging.

To confirm the precision of the high-resolution data, synthetic records from the two wells depicted in Figure 10 were utilized to compare with the processed data near the same wells. Figure 11A shows the logging data from Well 1 and its synthetic seismic record with a wavelet of dominant frequency 35 Hz. The high-resolution data obtained through the conventional spectral modeling approach and the novel method introduced in this paper were extracted in close proximity to the well in order to be compared with the synthetic seismic record, as shown in Figure 11A (3) and (4), respectively. The waveforms of both methods closely resemble the synthetic record at the targeted stratum position, and both methods exhibit a strong response from the reservoir. However, in the upper portion of the targeted stratum, the spectral modeling method does not match well with the synthetic record in terms of the time-shift and amplitude bias. The correlation coefficient of the method proposed in this paper is 0.71. While the waveform recovery method in this paper is more accurate and better matched with the synthetic record, and its correlation coefficient reaches 0.85. Figure 11B shows the synthetic record of Well 2 and the comparison of high-resolution seismic data near the well. It can be observed that the waveform generated by the traditional method in the targeted stratum poorly matches the synthetic seismic record, especially in stratum contact (P1m) where a significant time shift is evident. While the method described in this paper can accurately capture the reservoir information. The seismic stratum and the interpreted logging stratum correspond better.

According to the correspondence of the bright spot response to the reservoir, seismic bright spot responses indicate the presence of reservoirs. We selected a time window of 8 milliseconds below the upper boundary and 8 milliseconds above the lower boundary of the targeted stratum in Figure 10. The maximum peak amplitude in the time window is then extracted from various seismic data, as shown in Figure 12. Figure 12A shows the amplitude property extracted from the original seismic data. The bright spot response is not visible at Well 1, and there is a weak bright spot response at Well 2. In Figure 12B, the amplitude property extracted using the spectral modeling method with a spatial continuity constraint is shown. The bright spot response at Well 1 is enhanced, and a new bright spot appears at Well 2. But in the southern area (indicated by the white dashed line), there is a structural interpretation error caused by the absence of energy convergence at the upper boundary of the targeted stratum. This error leads to the occurrence of false bright spots and an unclear pattern. Figure 12C shows the amplitude property extracted from the data processed using the method described in this paper. The bright spot response is observed in Well 1 and Well 2, which is consistent with the forward analysis and logging interpretation. The explanation of the lower right corner of the figure is accurate, and the regularity of the highlighted amplitude is stronger.

High-resolution seismic data also contributes positively to the detection of minor faults. To assess the effectiveness of this approach in identifying minor defects, a particular area within the SN district in the Sichuan Basin was chosen for evaluation. The work area is characterized by a syncline structure, with the east and west wings exhibiting upturned formations and faults. Reservoirs primarily consist of lithological formations and are located in the central part of the syncline. Minor faults in this area have significant implications for reservoir reconstruction. The characterization of minor faults provide a foundation for accurate reservoir prediction, but its imaging is blurred due to the low resolution of the seismic data. The dominant frequency of seismic data is approximately 26 Hz, which makes it difficult to detect minor faults. Figure 13A shows the original seismic section, which demonstrates poor seismic resolution. So, minor faults and micro-tectonic morphology are similar, making it difficult to accurately identify them, especially in the y1 area. Figure 13B shows the seismic section processed using the method described in this paper. The figure shows the internal section of the syncline. The minor fault shown in Figure 13B is precisely delineated, with clear indication of the fault’s orientation and angle of inclination. Currently, seismic data has identifiable minor fault breaks of approximately 8 milliseconds. Based on the velocity of around 4,000 m/s here, the identifiable minor fault break is approximately 16 m. This confirms the accurate identification ability of the proposed method.

In light of these insights, the coherence property along the stratum interface was extracted, as shown in Figure 14. The faults on both sides of the work area developed, but the reservoir was primarily located in the syncline. So, the characterization of faults within the syncline was more significant, while micro-faults were indistinct within the syncline. Figure 14A shows the coherence map extracted from the original seismic data. The wells y1 and y2 in the work area are industrial gas wells, and the imaging logging shows that both wells have faults. There are faults in well y2, but no faults in well y1 on the original seismic coherence map. Figure 14B shows the coherence property of high-resolution data processed using the proposed method. The two wing faults can be clearly described, and the smaller faults are more visible within the syncline. The development of micro-faults can be clearly identified at well y1, which confirms the effectiveness of the proposed method in characterizing micro-faults.