To validate the accuracy of the proposed methodology, we acquired 1,100 km² of seismic data and complete well data (logs, drilling, and interpretation results) from nine producing wells in the Yingmaili area in collaboration with oilfield operators. A feasibility analysis of these foundational datasets is critical for determining the viability of subsequent experiments. Table 1 summarizes the feasibility analysis conducted on the acquired data.

The acquired raw seismic data are in the time domain, with a trace spacing of 20 m × 20 m and a fold of 80. No time–depth conversion was applied. However, horizons interpreted in the time domain were subsequently converted to the depth domain for subsequent evaluation of sand body thickness. The figures in this study were generated using the iLOOP+ software, a reservoir geophysical platform developed by the research team.

Following the analysis of the foundational data, we conducted a 2-week investigation into the genetic mechanisms of lacustrine beach-bar thin sand bodies in the Cretaceous Shushanhe Formation within the Yingmaili area (Xia et al., 2019; Zhao et al., 2022; Kuang et al., 2017; Hu et al., 2025; Chen et al., 2024; Chen and Zong, 2022). As shown in Figure 1, an experimental workflow integrating the proposed methodology was subsequently developed.

A step-by-step explanation of the experimental process in Figure 1 is provided below.

Thin sand body pattern analysis and seismic response characterization. Sedimentary patterns and typical facies classifications from nine wells in the target interval were analyzed to identify the characteristic seismic responses of different sandstone–mudstone configurations. A sensitivity analysis of multi-parameter combinations was conducted to identify the key parameters exerting the greatest influence on the seismic response, guiding the subsequent workflow stages.

The concept of the image pyramid, initially proposed by Crowley (1981), establishes an intuitive multi-scale representation system by decomposing data into a hierarchical sequence of images preserving high-resolution details at the base while capturing low-resolution macroscopic trends at the apex. A decade later, Freeman and Adelson defined the theory of steerable filters, enabling linear combinations of base filters in arbitrary directions with full rotational capabilities (Chen, 2021). The steerable pyramid method employs advanced transformation techniques to decompose images into multiple scales and applies directional filters at each level (Morozov and Ma, 2009; Zhang et al., 2024). Figure 2 illustrates the decomposition principle of this approach. From a geophysical perspective, the method meticulously dissects input seismic data into layered images (termed “levels” in this study), where each level encapsulates geological information across different scales and orientations: finer scales reveal high-frequency details such as faults, fractures, and geological boundaries, while broader scales delineate sedimentary trends and structural frameworks. The steerable pyramid transform decomposes the seismic data into a series of multi-scale, multidirectional sub-bands. At each decomposition level, it generates (1) a low-resolution parent band (P) representing the large-scale background; (2) a set of oriented bandpass bands (B) capturing directional features such as sand body edges at that scale; and (3) a high-pass residual band (Q) containing the finest details. This decomposition enables targeted analysis and enhancement of specific directional components relevant to thin sand body geometry.

The steerable pyramid method significantly increases the accuracy of fault identification and geological boundary characterization, providing critical technical support for facies modeling, reservoir modeling, and seismic inversion.

In geophysical inversion, it is standard practice to incorporate prior information to constrain the inversion process, ensuring that the results conform to established physical and mathematical principles while adhering to trends predicted by geological models (Bai et al., 2025). This study employs integrated multisource facies models as key constraints; however, the inherent uncertainties in sedimentary facies characterization may introduce errors when applied as rigid constraints. Consequently, we implemented a soft-constrained approach for the facies-constrained seismic inversion. This methodology maintains alignment with seismic responses while ensuring geologically reasonable inversion outcomes, even when facies characterization is imperfect.

Conventional geostatistical inversion employs variogram functions to characterize spatial correlations within reservoir sand bodies. However, when these sand bodies exhibit rapid lateral variations and the available well control is sparse, variograms often fail to accurately capture the required spatial relationships, which can introduce significant errors into the inversion results. To address this methodological limitation, this study introduces an innovative approach that substitutes the variogram with a spatially more stable probability density function (PDF). The probability density function, derived from well-log statistics, is used as a constraint within the geostatistical stochastic optimization seismic inversion framework. This function is relatively straightforward to compute and demonstrates inherent stability across varying well densities. Consequently, this methodology establishes a novel inversion algorithm grounded in the depositional configurations of thin-bedded sandstone and mudstone sequences as its fundamental geological template.

A moving-window approach based on seismic horizon interpretation is applied to calculate segmental variance, which defines the probable distribution range of acoustic impedance. In scenarios with sufficient well data, the extreme values derived from stratified statistics can be directly adopted as the spatial constraint distribution. Subsequently, a stochastic algorithm is employed to derive an acoustic impedance sequence that satisfies both these statistical characteristics and the seismic response. Beyond enabling the inversion of acoustic impedance sequences, this method can also be extended to the inversion of other reservoir petrophysical parameters.

Facies-constrained seismic inversion technology essentially represents an advanced development of well-log constrained stochastic optimization inversion. This method uses statistical characteristics of well impedance and sedimentary facies parameters to establish constraints, eliminating the need for complex modeling processes. By addressing the inverse problem through a forward modeling framework, implementing constraints becomes more straightforward.

We implemented an iterative stochastic inversion workflow constrained by facies models using the following procedure. As a constraint, the initial inversion cycle utilized a preliminary facies model, established based on regional geological understanding and geophysical data. This facies model was then refined by integrating the first round inversion results with well-based facies interpretations. The optimized facies model subsequently served as the constraint for the next inversion cycle. This iterative process continued until the discrepancy between the facies model and the inversion outcomes was minimized (Morozov and Ma, 2009). Within this framework, the geological facies model, inversion data volume, and original seismic data form a ternary closed-loop system; the detailed implementation workflow is illustrated in Figure 3. Specifically, this system creates an optimized workflow in which the facies model constrains the inversion process, while the inversion results in turn drive updates to the facies model. Ultimately, this ensures that both the facies model and the inversion results conform to geological principles while remaining consistent with the seismic response characteristics. The aforementioned facies-model constraint serves as a parametric constraint within the inversion workflow, while the probability density function constitutes an algorithmic constraint in the stochastic optimization inversion method. These represent two distinct types of constraints. The integration of facies-model constraints with the novel stochastic optimization inversion represents a significant innovation of this study.

Situated in the western section of the North Tarim Uplift in West China, the study area, the Yingmaili region, is flanked by the Wensu uplift to the west and the Halahatang depression to the east (Figure 4) (Zhang et al., 2018; Sun et al., 2020). During the initial sedimentary phase of the Cretaceous Shushanhe Formation, a proximal fan delta lake system with small sand bodies emerged in the northern region, while the southern region saw the development of a distal braided river delta lake system with larger sand accumulations. The middle sedimentary period of the Shushanhe Formation marked the peak of lake flooding, characterized predominantly by shallow lake mudstones. In its late sedimentary phase, the formation’s dynamics were influenced by a blend of near-inundation paleo-lift and east–west proximate Yingmaili submarine low uplift, leading to the formation of extensive shallow lake beach-bar deposits in the Yingmaili–Donghetang area. The Shushanhe Formation represents a thin-bedded sand body reservoir within a lacustrine beach bar, situated against a backdrop of northward-dipping monoclinal structures. This formation is notably capable of facilitating long-distance communication with distant oil sources via unconformities and fractures. The upper part of the structure features an upward-dipping, sharply delineated extinction zone of the sand body, identified as a hydrocarbon-bearing area, with the sandstone thickness ranging between 10 m and 25 m. This lithological reservoir is estimated to contain a billion tons of reserves (Chen et al., 2019). The seismic survey spanning the study area encompasses 1,100 km², including nine wells; it exhibits a low signal-to-noise ratio (SNR) characterized by a predominant frequency of 22 Hz and a bandwidth ranging from 4 Hz to 60 Hz within the targeted interval. The presence of thin-bedded features alongside seismic data with low SNR and subdued dominant frequencies presents substantial challenges to accurately characterizing the clastic sand reservoir.

The seismic survey in the study area covers approximately 1,100 km² with a sampling interval of 2 ms. The dominant frequencies vary significantly, ranging from approximately 20 Hz to 43 Hz, while the central frequency is 35 Hz. Consequently, a 35 Hz zero-phase Ricker wavelet was employed for this project. The Ricker wavelet comprises one central peak and two side lobes, as illustrated in Figure 5.

Prior to establishing models based on depositional patterns, petrophysical parameter statistics were conducted for the target interval. Excluding two wells with anomalous data in the study area, sonic transit time (DT) analyses for sandstones and mudstones were performed on the remaining seven wells. As shown in Figure 6, the DT distributions of sandstones and mudstones exhibit substantial overlap, with minimal distinction observed in their frequency distribution plots. The average DT value for sandstones is approximately 260.26 μs/m, while that for mudstones is approximately 269.8 μs/m. Density logs were available from three wells, all of which were used to analyze sandstone and mudstone densities. Based on the statistical characteristics presented in Figure 7, the density values were taken as 2.402 g/cm³ for sandstone and 2.414 g/cm³ for mudstone.

Based on the compiled boundary conditions and acoustic characteristics of the target interval in the Yingmaili Block, this study first established a sandstone wedge model to analyze the seismic response characteristics of the individual sand bodies. As illustrated in Figure 8, which comprises four subplots corresponding to mudstone interlayer thicknesses of 2 m, 4 m, 6 m, and 8 m, the model demonstrates relatively low seismic resolution for sand body detection. Root mean square (RMS) attributes were extracted from each pattern model. Figure 9 shows that when the overall thickness of the wedge with mudstone interlayers is thin, the seismic reflection energy exhibits minimal variation between models, with no distinct waveform differentiation, rendering seismic identification unfeasible. However, as the individual sand body thickness increases, the energy contrast increases, and waveform variations become pronounced, with clear seismic responses emerging near the mudstone interlayers. Under these conditions, individual sand bodies become seismically identifiable. A specific data analysis revealed the following detection limits,

With a 2 m mudstone interlayer, individual sand bodies thinner than 19 m and stacked sand bodies thinner than 7 m are undetectable.

Given the multifactor influences on the seismically identifiable thicknesses of sand bodies, we developed advanced pattern models that incorporate multiple geological factors. An analysis revealed that the sandstone-mudstone configurations can be categorized into five fundamental types across all sedimentary facies.

Thick configuration (∼25 m): Box-shaped, with 1–3 mudstone interlayers exhibiting diverse spatial distributions.

These pattern models are detailed in Figure 10. Forward modeling was conducted for each configuration, followed by a sensitivity analysis of the extracted seismic attributes. Using the maximum amplitude as an example, attribute variations across the configuration types were analyzed. The results indicate the following sensitivity hierarchy for thin sand bodies in the Shushanhe Formation in the Yingmaili area: sand thickness > interlayer position > number of interlayers. Five attributes were selected for the sensitivity analysis: maximum amplitude, RMS, integral of the absolute amplitude, total waveform energy, and amplitude kurtosis. The analytical results are summarized in Table 2. The selection of the five sensitive seismic attributes was based on preliminary reservoir prediction research conducted in the study area. These attributes represent the most frequently employed predictive parameters in prior academic studies. Their sensitivity ranking was established through a comprehensive statistical analysis of each mechanism model, in which a higher quantitative sensitivity value indicates better attribute sensitivity, and a lower value indicates poorer sensitivity. Finally, a corresponding sensitivity radar chart integrating all datasets was constructed, as illustrated in Figure 11.

As individual sand body thickness increases, the energy contrast increases, the energy contrast is increased, and waveform variations become more pronounced, enabling the effective seismic identification of single sand bodies. By populating an equal-thickness sandstone model with acoustic parameters representing different sedimentary rhythms, the seismic response effects of various rhythmic patterns were investigated. When the sand body thickness was progressively reduced from 50 m to 1 m, a comparative analysis of the waveform and energy characteristics between the rhythmic models and homogeneous models led to the following conclusions: when individual rhythmic sand bodies are thinner than 25 m, the reflection energies remain essentially consistent with the minimal rhythmic influence on the reflection intensity; however, when individual rhythmic sand bodies exceed 25 m in thickness, the composite rhythms exhibit a slightly stronger seismic reflection intensity than other rhythm types.

According to the steerability theorem, given an N th order filter, the number of basis filters M must satisfy the following inequality (Equation 1): M ≥ N + 1 N + 2 2 . (1)

Because we employ second-order filters, N = 2. The theorem implies that six basis filters are required.

To determine suitable orientations for the six basis filters, analysis can be conducted through an icosahedron. Prior to 3D filter design, a three-dimensional coordinate system must be defined, as shown in Figure 12 (Equations 2–5). Given the application of the basis directional filter bank in the frequency domain, coordinates are defined within the 3D frequency–wavenumber domain. Let vector V j a represent the symmetry axis of the basis filter B j , with its direction described by angles   θ j and φ j . Let vector V extend from the origin ( k x = k y = k z = 0 ) to point M, with its direction similarly described by angles θ and φ . Let Ω denote the angle between vectors V j and V .

We have also defined the direction cosines for each coordinate axis: k x , k y , k z , α = sin ⁡ θ ⁡ cos ⁡ φ , β = sin ⁡ θ ⁡ cos ⁡ φ , γ = cos ⁡ θ , (2) α j = sin ⁡ θ j ⁡ cos ⁡ φ j , β j = sin ⁡ θ j ⁡ cos ⁡ φ j , γ j = cos ⁡ θ j . (3)

In the three-dimensional case, we define six directional filters, β j θ , φ = cos 2 ⁡ Ω   θ , φ , θ j , φ j , j = 0 , 1 , … 5 . (4)

The expression for B j can be directly computed using the dot product V · V j , β j θ , φ = V · V j V · V j 2 . (5)

In the three-dimensional case, the impulse responses of the six filters are functions of θ and φ . Because the filters are symmetric, each possesses two maxima. For instance, the basis filter β 0 has maxima at θ = 90 ∘ , φ = 45 ∘ , with another maximum located at the antipodal point θ = 90 ∘ , φ = 225 ∘ . Given the symmetric nature of the basis filters employed, our analysis focuses on the angular range of θ = 0 ∘ − 90 ∘ and φ = 0 ∘ − 180 ∘ , as shown in Table 3.

After determining the decomposition levels and directional filter parameters, steerable pyramid processing is employed, guided by geophysical insights from forward modeling analysis to enhance directional feature identification and clarify structural information. The filter coefficient is set to 0–0.6. Planar enhancement processing is prioritized due to its ability to meet these requirements, ultimately yielding steerable pyramid results that align with the preset goals. The steerable pyramid method excels in multidirectional and multi-scale decomposition and processing, offering flexibility in selecting output sub-band datasets at various scales. Figure 13 displays the decomposed and optimized sub-band volumes from Level 0 to Level 5. Depending on specific geological needs, these sub-bands can be partially or fully reconstructed, thereby providing a rich set of data sub-bands to support the analysis of individual seismic datasets.

Spectral analysis was performed on the original seismic data and various reconstructed datasets from steerable pyramid processing. As shown in Figure 14, the spectral characteristics of the fully stacked seismic data are nearly identical to those of the original seismic data, indicating that the processed information is largely preserved. The stacked results of sub-bands 0, 1, and 2 exhibit higher dominant frequencies, which aid in identifying geological and structural information that is difficult to discern from the original seismic data. Consequently, the combined sub-band 0 + 1 + 2 data can be used as the foundational dataset for facies-constrained stochastic optimization inversion.

Through steerable pyramid processing, the dominant frequency of the target interval’s seismic data was enhanced from 22 Hz to 35 Hz. Based on the Rayleigh criterion and the average P-wave velocity (3,798 m/s) calibrated by well logs in this area, the vertical resolution of the seismic data improved from approximately 43 m before processing to approximately 27 m after processing. The resolution enhancement factor is approximately 1.59 times, which provides a critical data foundation for the subsequent identification and characterization of thin-bedded beach-bar sand bodies, whose thicknesses are generally less than 25 m in the study area.

Relative to the original seismic data, the reconstructed pyramid data maintain a consistent seismic frequency while enhancing the clarity and definition of the seismic homogeneous axis characteristics. Compared with the original seismic data, the medium–high-frequency, sub-band data derived from pyramid processing exhibit increased seismic frequencies, providing clearer structural insights and additional seismic details, as illustrated in Figure 15.

The steerable pyramid technique extracts planar attributes from the raw seismic data, significantly amplifying the directional qualities of the data. Subsequent steerable pyramid two-dimensional (2D) processing can further sharpen these directional qualities, aligning the analysis with geological exploration objectives, as depicted in Figure 16. As observed in Figure 15a, the original seismic data extraction of the planar attributes clearly outlines the contours of the geological formations. In Figure 15b, the locally reconstructed seismic data from sub-bands 0, 1, and 2 exhibit planar attributes with enhanced directional qualities, clearer details, and accentuated boundary delineations.

The advancing and receding patterns of the delta front are rendered significantly more discernible in the seismic profiles, with the transition zone between the delta and lacustrine facies exhibiting a consistent response in the original seismic attributes. In the profiles processed using the 3D steerable pyramid technique, the structural delineations are exceedingly clear and pronounced, especially at junctures where facies abruptly transition. Simultaneously, the 2D enhanced attribute map highlights sudden shifts in the directional attributes at the phase boundaries, which are attributed to the refined directional features. Utilizing the integrated profile features, the precise delineation of the facies boundary between the delta front and the lacustrine facies is markedly enhanced, as illustrated in Figure 17.

The research fully leverages the integration of well and seismic data. During attribute analysis, the spatial distribution of seismic attributes derived from steerable pyramid-processed data was examined. Guided by well-log characteristics at the well points and 3D steerable pyramid profiles, a combined well–seismic and map-section analysis was implemented. The planar distribution of sedimentary facies was determined based on the maximum amplitude attribute distribution of each sand group within the target interval, as illustrated in Figure 18. Five sedimentary microfacies were delineated within the target intervals. The resulting planar distribution of sedimentary facies for individual sand layer groups was then utilized as the initial model for the facies-constrained volume. Subsequently, by correlating the sedimentary facies maps with seismic attribute profiles, the initial model was transformed into a 3D volume, thereby integrating sedimentary facies and seismic facies and constructing the facies-constrained model for inversion. The methodology for building this facies-constraint model is illustrated in Figure 19. The sedimentary facies distribution map was compiled at a scale of 1:2500, with a theoretical cartographic accuracy of approximately 0.25 m. However, the actual delineation accuracy of facies boundaries is primarily governed by the lateral resolution of the seismic data, meaning the effective spatial resolution falls short of 0.25 m. As the vertical resolution of the processed seismic data has been enhanced, an improvement in spatial resolution compared to the original data can be expected. The accuracy of phase calibration is determined by the alignment between seismic horizons and well-log stratigraphic divisions, with geological constraints established based on well-log curve trends.

The target Cretaceous Shushanhe Formation in the study area comprises lacustrine beach-bar sandstone deposits characterized by sandstone–mudstone interbeds. Prior to delineating thin sand body reservoir boundaries, non-reservoir sections must be excluded, necessitating pre-inversion screening of sensitive parameters that effectively characterize the reservoir. Sensitive inversion parameters were optimized from well-log characteristic curves of the sand bodies. As shown in Figure 18, when impedance fails to adequately distinguish sandstone from mudstone, spontaneous potential (SP) and natural gamma ray (GR) logs provide better differentiation. Therefore, petrophysical parameter SP or GR curves must be inverted jointly to identify effective reservoirs. Figure 20 indicates that effective reservoirs are screened using thresholds of SP < 140 and GR < 85. However, due to inconsistencies in SP baseline and value ranges across wells, pre-inversion processing requires baseline correction, value-range standardization, and normalization of SP logs for all wells in the study area. Additionally, given the limited availability of density logs in the study area, velocity is directly inverted as a pseudo-impedance volume.

To address the challenge of limited well coverage in the target area, this study fully leverages the complementary advantages of well logs, seismic data, and geological interpretations to achieve high-precision inversion results. Consequently, implementing a joint well–seismic–facies integrated seismic inversion is crucial. This approach combines the superior vertical resolution of well-log data with the continuous lateral coverage provided by 3D seismic data to achieve the set high-resolution inversion objectives.

Given the limited number of wells and the presence of thin reservoirs within the study area, the reservoir prediction uses geostatistical stochastic inversion techniques, with a particular emphasis on facies-constrained inversion derived from sedimentary models. This approach entails applying sedimentary insights to the inversion process, guided by well-log data; thereby enhancing both the reliability and the precision of reservoir predictions.

The facies-constrained seismic inversion technique represents an advanced evolution of the logging-constrained stochastic optimization inversion approach. This innovative method harnesses the statistical properties of well impedance and sedimentary facies parameters to set constraints, streamlining the process without necessitating complex modeling. By transforming the inversion challenge into a solvable problem, this technique simplifies the application of constraints. The procedural flowchart is illustrated in Figure 21.

The target Cretaceous Shushanhe Formation in the Yingmaili area comprises lacustrine beach-bar sand bodies, characterized by interbedded sandstone and mudstone. To accurately delineate the effective reservoir sand body boundaries, non-reservoir zones must be excluded. Prior to seismic inversion, we used a screening process to identify sensitive parameters characterizing effective reservoirs. Based on a comprehensive analysis, the spontaneous potential log was optimized as a sensitive parameter for inverting thin sand body reservoirs.

Using stochastic optimized inversion techniques, this approach mitigates inversion errors stemming from limited well control by enhancing the geological model constraints, thereby significantly boosting the sand layer identification precision between wells. To facilitate the facies-constrained inversion, it is essential to construct spatial depositional facies belts via comprehensive multi-attribute and geological analyses to develop an accurate depositional facies model. Determining the constraint space for each facies type in this model is crucial to ensure the correlation of data within the same facies category across a region. This approach guarantees that the inversion outcomes align with the seismic attributes and the depositional patterns, yielding geologically meaningful insights, as demonstrated in Figure 22. Based on the analysis of continuity and discontinuity in thin sand body inversion results, the outcomes incorporating facies-constrained inversion demonstrate superior performance in this regard.

As can be seen from Figures 22, 23, the stochastic optimization inversion results incorporating facies-model constraints demonstrate superior accuracy when characterizing thin sand body reservoirs, especially with respect to their top and bottom boundaries. In terms of the planar attribute distributions, the facies-constrained stochastic optimization inversion provides a sharper delineation of the thin sand body boundaries in the map view. We also identified amplitude-dependent variations in the inversion results, a phenomenon consistent with theoretical expectations and geological principles. The synthetic forward model closely aligns with the original seismic data, with minimal local discrepancies. Furthermore, the consistency in the sand body distribution patterns and sedimentary characteristics across multidirectional inversion profiles provides an additional validation of the inversion reliability.

As the constraint boundaries within the inversion process become more relaxed, the model experiences fewer restrictions, leading to smoother generated feature curves at the expense of diminishing details with respect to the limited logging characteristics. Because the reservoir in the study area is thin and non-homogeneous, the accuracy of the thin interbedded identification by the inversion will be impacted if a set spatial constraint is consistently applied. In this article, the constraint-space-finding workflow of different facies models was constructed by subphase type to narrow the perturbation range of the inversion values, which can effectively solve the problem of loose inversion constraints.

Well-based validation of the facies-model-constrained stochastic optimization inversion for thin sand bodies was performed using two typical development wells (YM106H and YM103) in the Yingmaili area as blind verification wells. The vertical distribution patterns of the sand bodies were compared with the inversion profiles, and the results are shown in Figure 24. Based on the time–depth relationship established through well-to-seismic calibration, the time-domain SP inversion volume was converted to the depth domain. Subsequently, using the reservoir cutoff value of SP < 140 (determined from the sensitivity analysis of reservoir inversion parameters discussed earlier), sand body thickness was separately calculated for wells YM103 and YM106HD. The validation results demonstrate prediction accuracy rates of 63% for sand bodies thicker than 4 m and 84% for those thicker than 10 m, with supporting data provided in Tables 4 and 5 (Equation 6). The calculation method for the sand body thickness prediction accuracy rate is as follows: R = 100 % − b − a a × 100 % , (6) where R is the conformity rate, a is the logging-interpreted thickness of thin sand bodies, and b is the inversion-predicted thickness of thin sand bodies. The discrepancy between the inversion-predicted and logging-interpreted thicknesses of thin sand bodies is defined as the prediction error thickness. The ratio of the prediction error thickness to the logging-interpreted thickness is termed the inaccuracy prediction ratio. The prediction conformity rate is then calculated as 100% minus this inaccuracy prediction ratio.

This case study demonstrates excellent predictive performance and confirms the practical viability of the proposed methodology for thin sand body reservoir prediction in various global subtle hydrocarbon reservoir plays.