Multi-constraint collaborative seismic imaging and interpretation method for fine characterization of small faults in complex coal mining areas

1 Introduction

The fine characterization of small faults in complex coal mining areas is directly related to roadway layout, working face design, and production safety, and constitutes a key issue in the geological assurance of coal mining (Li D. et al., 2024). Faults that remain undetected or are mislocated at the exploration stage carry consequences well beyond a loss of geometric accuracy: mining across or adjacent to a fault alters the local stress state and can trigger fault reactivation, which is a recognised source of coal bursts under the superposition of mining-induced static and dynamic stresses (Cai et al., 2021; Chen et al., 2018), and can open pathways between the working face and confined aquifers in the floor, leading to water inrush (Song and Liang, 2021; Meng et al., 2012). The reliability with which faults of small throw can be located is therefore a safety constraint on mining deployment rather than a refinement of interpretation alone. For mining areas with relatively deep burial, small fault throws, and complex surface conditions, the seismic responses of small faults and local no-coal zones are usually weak and easily affected by shallow low-velocity structures, the distribution of surface obstacles, and the irregularity of the observation system (Peng et al., 2026; Dong et al., 2023). As a result, conventional seismic interpretation has limited capability for identifying micro-tectonic features and seam morphology variations (Cheng et al., 2023; Seelam et al., 2024). The mining-revealed data of the Daqiang Coal Mine indicate that the actual structural complexity is significantly higher than that previously interpreted from seismic data, suggesting that existing results can no longer satisfy the requirements for fine characterization of small faults, no-coal zones, and seam occurrence morphology in subsequent mining deployment. Therefore, establishing a high-confidence seismic imaging and interpretation method for complex coal mining areas has clear theoretical significance and engineering value (Li D. et al., 2024; Peng et al., 2026).

Progress relevant to the fine identification of small faults can be organised into four technical branches—acquisition, regularization, imaging and interpretation—which correspond to the successive links of the seismic workflow and are compared in Table 1. In the acquisition branch, high-density and wide-azimuth observation improves the completeness of wavefield recording and the preservation of structural responses, and provides the sampling basis for transparent exploration of complex coal mining areas (Peng et al., 2026; Dong et al., 2023; Cheng et al., 2023); work here is concerned mainly with bin size, fold and azimuth distribution, and is normally evaluated by geometry statistics rather than by the imaging result. In the regularization branch, offset-vector-tile (OVT) organisation and five-dimensional interpolation mitigate the non-uniform offset–azimuth sampling produced by variable observation under complex surface conditions, and supply a stable data geometry for azimuth-preserving imaging (Cary, 1999; Zhang et al., 2021; Vermeer, 2012; Zeng et al., 2023); this step is usually treated as self-contained, its output not being explicitly constrained by the subsequent migration operator. In the imaging branch, high-precision prestack migration improves the focusing and repositioning of reflection and diffraction energy adjacent to complex faults, and thereby enhances the clarity and continuity of small-scale fracture images (Han et al., 2022), although such studies generally assume that the input gathers are already geometry-consistent. In the interpretation branch, structure-oriented attributes, multi-scale curvature, coherence analysis and multi-attribute fusion have been applied to fault-boundary enhancement, characterization of mined-out areas and prediction of coal-body structure (Seelam et al., 2024; Li et al., 2023; Cui et al., 2023; Zhou et al., 2024; Ercoli et al., 2023; Li H. et al., 2024); the stability of this branch depends directly on the quality of the migrated volume, which is controlled by the three upstream branches.

TABLE 1. Comparison of representative approaches in the four technical branches relevant to fine identification of small faults.

Table 1 compares representative studies of the four branches in terms of the technical link that is optimised, the scale of the targets resolved, and whether the workflow is closed by engineering revelation. The comparison shows that most published work improves one or two links, whereas the coordinated optimisation of all four links, verified against mining-revealed data, is seldom reported for coal mining areas that combine deep burial, small fault throw and dense surface obstacles. This indicates that the identification capability for small structures in structurally complex areas depends on the coordinated optimisation of acquisition, processing, imaging and interpretation, rather than on the local enhancement of any single technical link.

Despite these advances, existing studies tend to focus on observation-system design, improvements in individual processing techniques, or attribute-enhanced interpretation as isolated optimizations (Li D. et al., 2024; Peng et al., 2026; Dong et al., 2023). For complex coal mining areas characterized by deep burial, small fault throws, dense surface obstacles, and shallow low-velocity structures, a closed-loop technical workflow that integrates high-density sampling, geometry-consistent regularization, velocity-field-constrained imaging, robust interpretation, and engineering-revealed verification is still lacking. Increased sampling density does not automatically translate to higher imaging resolution: if the regularization and migration processes cannot effectively preserve offset–azimuth consistency, the sampling gain cannot be converted into an imaging gain. Meanwhile, under conditions of weak structural response and overlapping anomalies, the stability of single-attribute interpretation remains limited, and the reliability of interpretation results in supporting mining deployment and production safety needs further improvement (Li D. et al., 2024).

Based on the above considerations, this study takes the Daqiang Coal Mine as the research subject and, with the goal of fine characterization of small faults, no-coal zones, and seam morphology, constructs an integrated technical workflow that combines high-density 3D acquisition, OVT-domain regularization, velocity-field-constrained azimuth-preserving prestack time migration, and multi-attribute fusion interpretation. The objectives are to improve the consistency of wavefield sampling and coverage under complex surface conditions, enhance the repositioning accuracy of reflections and diffractions associated with small faults, and improve the continuity and credibility of fault-boundary interpretation. The contribution of this study lies not in proposing an isolated new formulation detached from engineering context, but in systematically reorganizing high-density sampling, OVT-domain processing, high-precision imaging, and multi-attribute interpretation under the multiple engineering constraints of complex coal mining areas, thereby forming a verifiable, reusable, and transferable fine-characterization method for small faults that provides high-confidence seismic-geological evidence for mining deployment and safety assessment.

Relative to existing integrated approaches, the specific increments of this work are threefold. First, the acquisition parameters are not selected empirically but are derived from the joint quantitative constraints of spatial anti-aliasing, coverage gain and coverage uniformity, and the anti-aliasing condition is evaluated over the full velocity and dip range of the study area instead of at a single parameter set. Second, OVT-domain regularization, velocity-field-constrained analysis and azimuth-preserving migration are coupled through one common data organisation, so that the offset–azimuth consistency recovered upstream is carried into the migration operator rather than being re-mixed at the imaging stage. Third, the chain is closed by engineering revelation, the interpreted fault system being checked against 41 boreholes and mining exposure, which allows the sampling gain to be traced quantitatively through the imaging gain to the interpretation gain. These three increments jointly distinguish the present workflow from integrations that optimise acquisition, processing or interpretation as separate links.

2 Materials and methods

The method developed in this study is organised as a single chain in which each module is constrained by the output of the preceding one, so that the sampling gain obtained at acquisition is carried forward rather than being lost at a later stage. The chain comprises five modules, and the quantities transmitted between them are as follows.

The first module is high-density 3D acquisition (Section 2.2), which uses a 24L6S150T orthogonal geometry with vibroseis excitation over 3–96 Hz and a variable-observation design in obstacle areas. It delivers to the next module a bin size of 5 m × 5 m, a nominal fold of 300, the coverage-uniformity coefficient, and a highest effective frequency of about 75 Hz. The second module is OVT/COV data organisation and five-dimensional regularization (Section 2.3.1), which parameterizes the traces by midpoint and half-offset following Equation 9 and reconstructs the missing geometry through Equations 10, 11. It transmits the midpoint, offset magnitude and azimuth, the OVT index, and the regularization operators. The third module is velocity-field-constrained analysis (Section 2.3.1), which builds the semblance velocity spectrum of Equation 12 and converts to interval velocity through Equation 13, passing the root-mean-square velocity field, the interval velocity, and the azimuth-dependent velocity used by the migration operator. The fourth module is azimuth-preserving Kirchhoff prestack time migration (Section 2.3.2), which computes traveltimes under the root-mean-square velocity approximation of Equations 14, 15 and performs the weighted diffraction stack of Equation 16, delivering the traveltime operator, the migration aperture and weighting function, and the image volume. The fifth module is interpretation, in which structure-oriented attributes along the T1 horizon are fused by principal component analysis and thresholded to delineate faults and no-coal zones; its outputs are the interpreted structural elements, which are then checked against 41 boreholes and mining exposure in Sections 3, 4.

A quality-control loop is applied across the chain rather than within any single module. Fold uniformity, gather flatness and residual moveout are examined after regularization and after velocity analysis, and any departure from the tolerances established in Section 2.2.4 is fed back to the corresponding upstream module before migration proceeds. This feedback is what keeps the offset–azimuth consistency recovered during regularization from being re-mixed at the imaging stage, and it distinguishes the present workflow from the sequential application of the same individual techniques.

2.1 Study area and seismic-geological conditions

2.1.1 Survey area location and geological tasks

The survey area is located at the boundary between Kangping County, Liaoning Province, and Tongliao City, Inner Mongolia. It has a polygonal planimetric extent with an area of approximately 4.8 km² (Figure 1; Table 2). The main mineable target seam in the area is the No. One coal seam, with a burial depth of approximately 740–1,189 m. The deep burial of the target seam, the high structural complexity, and the development of local small faults, seam thinning, and no-coal zones exert direct controls on subsequent mining deployment.

TABLE 2. Survey area boundary coordinates (Beijing 1954 coordinate system, 1985 national elevation datum, Gauss 3° projection, central meridian 123°).

The seismic exploration tasks in this work have a clear dual nature of structural interpretation and engineering service, including the following three aspects.

• Identify all faults in the No. 1 coal seam with throw greater than 10 m, with a verification rate of no less than 90%; identify the nature, attitude, and extent of faults with throw ≥ 5 m, with a verification rate greater than 80%; and provide seismic-response interpretations for break points and suspected points with throw less than 5 m.
• Determine the floor undulation morphology of the No. 1 coal seam and produce a floor contour map with a contour interval of 10 m, with depth errors controlled within ±15 m.
• Identify folds with amplitude ≥20 m, as well as seam-thinning zones and no-coal zones with width ≥20 m within the area.

The core of seismic exploration in this area is not the conventional clarification of regional structural patterns, but rather the high-precision, small-scale, and strongly constrained fine interpretation centered on the target coal seam. In particular, the requirement to identify faults with throws on the order of 5 m, as well as smaller break points and suspected points, dictates that this study must simultaneously focus on vertical resolution, lateral continuity, structural-imaging sharpness, and interpretation reliability.

2.1.2 Analysis of seismic-geological conditions

The surface of the survey area is a sandy plain with villages, forest belts, and rivers (Figure 2). The shallow low-velocity layer has a thickness of approximately 3–5 m, the groundwater table is approximately 7–8 m deep, and the near-surface seismic-geological conditions are complex. At depth, there is a significant acoustic-impedance contrast between the No. 1 coal seam and the surrounding rocks, which can form a stable target reflector T1 (Figure 3).

Theoretically, the reflected energy of the coal seam is mainly controlled by the acoustic-impedance contrast between the two sides of the interface. Under the condition of normal incidence, the reflection coefficient is given by Equation 1:where Z is acoustic impedance, ρ is density, and v is the P-wave velocity. This expression shows that as long as a stable impedance contrast exists between the roof, floor, and surrounding rocks of the coal seam, a traceable reflection interface can be formed, which is the physical basis for using T1 as the main interpretation horizon.

R 0 = Z 2 − Z 1 Z 2 + Z 1 , Z = ρ v (1)

In practice, seismic recordings are not pure interface responses but the result of convolving the reflection-coefficient series of the strata with the source wavelet. The 1D convolution model can be written as Equation 2:where s t is the seismic trace, w t is the wavelet, r t is the reflection-coefficient series, and n t is the noise. This expression shows that if the shallow region is complex, the noise is strong, or sampling is inadequate, then the target reflection in the convolution result will be broadened, contaminated, or even aliased. Therefore, the essence of the complex seismic-geological conditions in this area is not limited to structural complexity but is manifested as a combined result of weakened effective signals, loss of high-frequency information, and superimposed spatial-sampling risk.

s t = w t * r t + n t (2)

2.2 High-density 3D seismic acquisition system

In response to the characteristics of the Daqiang Coal Mine—deeply buried target seam, small fault scale, complex shallow low-velocity structure, and densely distributed surface obstacles—this study designed a fully digital high-density 3D seismic observation system with the goals of improving wavefield sampling density, coverage uniformity, and data fidelity. The observation-system design, source-and-receiver parameter selection, system adjustment in obstacle areas, and quality control all serve subsequent OVT-domain regularization, prestack time migration, and fine interpretation.

2.2.1 Observation system design

Based on the seismic-geological conditions of the study area, the burial depth of the target seam, and the requirements for fine identification of small faults, a 24-line, 6-shot orthogonal-bundle observation system with mid-point excitation was used in this acquisition. Specific parameters are listed in Table 3. The system uses 24 lines × 150 channels per line for a total of 3,600 channels; common-depth-point (CDP) grid of 5 m × 5 m; fold of 300; trace spacing of 10 m; receiver-line spacing of 60 m; shot spacing of 10 m; and shot-line spacing of 30 m.

TABLE 3. Observation system parameters.

Compared with the conventional 3D seismic acquisition in 2006, the CDP bin in this acquisition was reduced from 10 m × 10 m–5 m × 5 m, increasing the planar sampling density by a factor of 4, while the fold was raised from 40 to 300. This observation system achieves small bins, high fold, and regular wide-azimuth sampling, the quantitative basis for which is analysed in Section 2.2.4.

2.2.2 Source and receiver parameters

Based on the shallow structure, target burial depth, and operational conditions of the study area, point and line tests were used to optimize the acquisition parameters. Excitation was performed using a BV620LF vibroseis truck (1 truck, 2 sweeps), with a sweep frequency of 3–96 Hz, sweep length of 16 s, and drive level of 70%. Reception was carried out using SD7X10Hz single-point digital geophones, and a Sercel 428XL digital seismic recorder was used with a 1 m sampling interval and 2 s record length.

This combination of parameters ensures the excitation and reception of effective reflection information from the deep target seam while balancing stable operation under complex surface conditions and the ability to preserve weak signals. Single-point digital reception enhances the recording capability of local wavefield differences, providing high-quality raw data for subsequent identification of small faults and fine imaging.

2.2.3 System adjustment in obstacle areas and quality control

Villages, forest belts, and other obstacles are widely distributed in the study area (Figure 4), which made it impossible to implement the regular grid in some local areas. To address this issue, an “in-area shifting and out-area densification” variable-observation design was adopted to jointly adjust the observation system in obstacle areas and their adjacent areas, so as to maintain coverage continuity and sampling uniformity. The adjusted observation system was analyzed and verified for coverage using KLseisII software, with results shown in Figures 5–10.

A total of 31,635 production physical points were completed in this acquisition, with a Class-A record rate of 85.88% and a finished-product rate of 100% (Table 4). The acquisition produced raw data with relatively uniform coverage, high signal-to-noise ratio (SNR), and good continuity. The results indicate that under obstacle-area conditions, local observation-system adjustments combined with whole-process quality control can effectively ensure the acquisition quality of high-density 3D data in complex coal mining areas.

TABLE 4. Quality indicators of the Daqiang 3D acquisition.

2.2.4 Quantitative design of high-density acquisition parameters

The requirements for fine characterization of small faults in complex coal mining areas are reflected in the integrated constraints on spatial-sampling capability, coverage gain, and sampling-geometry stability. For an exploration area such as the Daqiang Coal Mine—with deep target seams, small fault throws, complex shallow low-velocity structures, and dense obstacles—the selection of acquisition parameters cannot remain at the level of empirical densification but should be based on a quantitative analysis of spatial anti-aliasing conditions and coverage constraints.

From the perspective of spatial sampling, high-dip reflections and diffractions adjacent to faults in 3D prestack imaging are most sensitive to the sampling interval. To ensure alias-free imaging of the target wavefield during migration, the spatial-sampling interval should satisfy the anti-aliasing condition:where Δ x is the effective spatial-sampling interval, v is the effective velocity near the target seam, f max is the highest effective frequency to be preserved, and θ max is the maximum dip of the target reflection or diffraction event. Taking the effective velocity near the target seam as v = 2,400 m/s, f max = 75 Hz, and θ max = 70°, we obtain Equation 4:

Δ x ≤ Δ x c r i t = v 4 f max ⁡ sin ⁡ θ max (3)

Δ x c r i t = 2400 4 × 75 × sin ⁡ 70 ° ≈ 8.5 m (4)

Equation 3 is an analytical constraint, so its implication for bin-size selection can be examined over the whole parameter range of the study area rather than at a single design point. Within the survey area the effective velocity near the No. 1 coal seam varies between about 2,200 and 2,600 m/s over the 740–1,189 m burial-depth range of the target seam, and the maximum dip of the reflections and fault-plane diffractions to be preserved varies between about 30° and 70°. Table 5 lists the critical sampling interval computed from Equation 3 over this parameter space for three levels of the highest effective frequency to be preserved, and Figure 11 shows the corresponding curve family together with the distribution of the critical sampling interval in the velocity–dip plane.

TABLE 5. Sensitivity of the anti-aliasing critical sampling interval to effective velocity, maximum dip and highest effective frequency. Each cell gives the range in metres obtained for effective velocities of 2,200–2,600 m/s.

The sensitivity analysis shows that the critical sampling interval varies approximately linearly with velocity and approximately as the reciprocal of the sine of the maximum dip, so that the two factors together move the constraint over a range of about 6.5–21.7 m across the parameter space examined. A 10 m bin satisfies the anti-aliasing condition only in the low-frequency, low-dip, high-velocity corner of that range: it is already violated at 75 Hz for dips steeper than about 50°, and at 90 Hz for dips steeper than about 40°. The 5 m × 5 m bin adopted in this study, by contrast, satisfies the condition over the entire parameter space, retaining a margin of at least 1.5 m even for the most demanding combination of 2,200 m/s, 90 Hz and 70°, which gives a critical interval of 6.5 m. The design point used above therefore does not represent a favourable special case but lies in the middle of the admissible range, which indicates that the choice of a 5 m bin is robust to the velocity and dip uncertainty of the study area and is more conducive to preserving small-fault boundaries and diffraction responses.

From the perspective of coverage gain, under random-noise conditions, the SNR gain after stacking can be approximated as Equation 5:where F is the fold. Taking the conventional 2006 3D observation system, with a fold of 40, as the reference, and with a fold of 300 for the present design, the relative SNR gain is given by Equation 6:

G S N R ≈ F (5)

G S N R , 2 G S N R , 1 = F 2 F 1 = 300 40 ≈ 2.74 (6)

As shown in Figure 12, the 300-fold coverage in this study can improve random-noise suppression capability by approximately 2.74 times relative to the original 40-fold coverage, which is conducive to enhancing the identifiability of weak deep reflections, small-fault break points, and local anomaly boundaries.

To further characterize the integrated sampling capability of different observation systems, a per-bin coverage-density indicator is introduced (Equation 7):where Δ x and Δ y are the bin dimensions in the two directions. For the present observation system, D c o v,2 = 300/(5 × 5) = 12.0 traces/m², while for the conventional 2006 3D observation system, D c o v,1 = 40/(10 × 10) = 0.4 traces/m². The ratio D c o v,2/D c o v,1 = 30. This indicates that the effective coverage density per bin in the present high-density acquisition is about 30 times that of the conventional 3D system. The corresponding parameter constraints, values, and computed results are summarized in Table 6.

D c o v = F Δ x · Δ y (7)

TABLE 6. Quantitative constraints for high-density acquisition design.

To quantify the effect of observation-system adjustment in obstacle areas, the coverage-uniformity coefficient (CVF) can be used (Equation 8):where F i is the fold of the i-th bin, F ¯ is the mean fold, and M is the total number of bins. A smaller C V F indicates more uniform coverage. For complex obstacle areas, the “in-area shifting and out-area densification” variable-observation design essentially aims to reduce C V F, maintaining the continuity of local coverage and the stability of sampling geometry.

C V F = 1 M Σ i F i − F ¯ 2 F ¯ 2 (8)

In summary, the selection of high-density acquisition parameters in this area has a clear quantitative basis: the 5 m × 5 m small bin satisfies the requirements for high-frequency information preservation and spatial anti-aliasing in complex structures; the 300-fold coverage significantly improves random-noise suppression; per-bin coverage density is increased by approximately 30 times relative to conventional 3D; and the variable-observation design in obstacle areas improves coverage uniformity and enhances geometric compatibility for subsequent OVT-domain regularization and prestack migration. Therefore, the high-density 3D observation system in this study is not a simple empirical densification, but a quantitatively designed result under the joint constraints of spatial sampling, coverage gain, and sampling-geometry stability.

2.3 OVT-domain regularization and azimuth-preserving prestack time migration

The key to processing high-density 3D seismic data in complex coal mining areas lies not in the local enhancement of a single processing module, but in the stable recovery of reflection–diffraction information related to small faults under conditions of complex near-surface disturbances, obstacle-area variable observation, and weak deep structural responses, and in effectively converting the sampling gain into an imaging gain. Based on the high-density sampling foundation established above, the key processing flow in this study is organized into three steps—OVT-domain regularization and reconstruction, velocity-field-constrained analysis, and azimuth-preserving Kirchhoff prestack time migration—forming a joint processing framework oriented toward fine characterization of small faults in complex coal mining areas. The theoretical foundations of this framework draw on OVT/common-offset-vector (COV) data organization, OVT-domain 5D interpolation and Fourier sparse reconstruction, semblance velocity analysis, Dix velocity conversion, and Kirchhoff prestack time migration. The contribution of this study is the targeted coupling and parameterized implementation of these mature methods under the multiple engineering constraints of complex coal mining areas.

2.3.1 Joint processing framework and key constraints

The main difficulties of data processing in this area include: variable observation in obstacle zones causes uneven sampling in offset–azimuth space; shallow low-velocity structures and strong surface-wave interference weaken effective deep reflections; the target seam is deeply buried and the diffraction responses adjacent to small faults are weak; and if regularization, velocity analysis, and migration are implemented independently, geometric inconsistencies accumulate and are amplified during imaging. To address these issues, this study adopts a joint processing framework that couples OVT-domain regularization, velocity-field-constrained analysis and azimuth-preserving Kirchhoff prestack time migration through a common data organization, in which the output of each step is used as an explicit constraint on the next, so that the offset–azimuth consistency recovered during reconstruction is preserved through velocity analysis and carried into the migration operator rather than being re-mixed at the imaging stage. The theoretical basis of each step is given below.

For a 3D observation system, source point s and receiver point r can be parameterized as midpoint–half-offset vector form (Cary, 1999):where m is the midpoint coordinate, h is the half-offset vector, h is the offset, and φ is the azimuth. This parameterization reorganizes the original 3D acquisition data into data units with similar propagation geometry, providing a consistent input geometry for regularization reconstruction and azimuth-preserving migration.

m = s + r 2 , h = r − s 2 , h = h x 2 + h y 2 , φ = tan − 1 h y h x (9)

In the regularization step, the irregular-sampling model in the frequency domain can be written as (Zhang et al., 2021; Liu and Sacchi, 2004):where d ω is the observed data at frequency ω, P Ω is the irregular-sampling operator, F − 1 is the inverse Fourier transform operator, c ω are the regular-grid Fourier coefficients, and n ω is the noise term. Following OVT-domain 5D interpolation with minimum weighted norm interpolation (MWNI)/Fourier sparse reconstruction, the regularization objective can be written as (Zhang et al., 2021; Liu and Sacchi, 2004):where the first term is the data-consistency term, the second term is the sparsity-constraint term, W is the frequency–wavenumber domain weight matrix, and λ is the regularization parameter. The role of this process is to recover the consistency of offset–azimuth sampling and to mitigate spatial holes, footprints, and local sampling distortions caused by variable observation in obstacle areas.

d ω = P Ω F − 1 c ω + n ω (10)

c ^ ω = argmin d ω − P Ω F − 1 c ω 2 2 + λ W c ω 1 (11)

For velocity analysis, the classical semblance velocity spectrum is used as the basic criterion (Taner and Koehler, 1969):where u j is the corrected amplitude of the j-th trace, M is the number of traces involved in stacking, N t and Δ t j is the correction determined by velocity v and zero-offset time t₀. To further estimate interval velocities from root-mean-square (RMS) velocities, the Dix velocity conversion is used (Dix, 1955):

S v , t 0 = Σ k Σ j u j t k + Δ t j v , t 0 2 M · Σ k Σ j u j 2 t k + Δ t j v , t 0 (12)

v int , n = v r m s , n 2 t n − v r m s , n − 1 2 t n − 1 t n − t n − 1 (13)

Therefore, the key points of Section 2.3.1 can be summarized as: OVT-domain regularization is responsible for recovering geometry-consistent sampling, semblance–Dix analysis is responsible for constructing a stable velocity field, and the two together establish reliable input conditions for subsequent azimuth-preserving migration.

2.3.2 Velocity-field-constrained azimuth-preserving Kirchhoff PSTM

After completing OVT-domain regularization and velocity-field construction, the imaging step uses azimuth-preserving Kirchhoff prestack time migration. Its core lies in the fact that, for data within the same OVT unit, the offset and azimuth information are explicitly preserved, so traveltime computation and weighted stacking can be performed under similar propagation-geometry conditions, thereby reducing the geometric aliasing caused by complex variable observation (Cary, 1999; Han et al., 2022).

For an image point x, the two-way traveltime from the source side to the receiver side can be written as (Schneider, 1978):

T x ; s , r = T s x ; s + T r x ; r (14)

In the time-migration framework, the commonly used RMS-velocity approximate form is:where ρ s and ρ r are the planar distances between the surface projection of the image point and the source and the receiver, respectively, τ is the one-way vertical time of the image point, and v x , φ is the imaging velocity field obtained from Section 2.3.1. The Kirchhoff-type prestack time-migration imaging expression is (Han et al., 2022; Schneider, 1978):where I x is the imaging result at image point x, W is a composite weight function (geometric spreading, dip compensation, aperture limitation and amplitude weighting), and d is the input data. This integral form itself is the classical Kirchhoff migration and is not a newly proposed fundamental migration formula. The methodological feature of this study is that the migration input data come from OVT-domain 5D regularization, and the traveltime computation uses an imaging velocity field consistent with the OVT azimuth organization, so the migration operator and data geometry maintain higher consistency. Prestack denoising, surface-consistent predictive deconvolution, OVT-domain five-dimensional regularization, velocity analysis and azimuth-preserving Kirchhoff prestack time migration were all implemented on the BGP GeoEast seismic data processing system.

T x ; s , r ≈ ρ s 2 v 2 x , φ + τ 2 + ρ r 2 v 2 x , φ + τ 2 (15)

I x = ∬ W x , s , r · d s , r , T x ; s , r d s d r (16)

2.4 Structure-oriented attributes and PCA multi-attribute fusion

The migrated volume obtained in Section 2.3 provides the input for interpretation. Because the seismic response of a small fault is weak and a single attribute is easily contaminated by noise and by overlapping anomalies, this study extracts a set of edge-sensitive attributes along the target reflector and fuses them by principal component analysis (PCA) before any fault is picked.

Attribute set. Structure-oriented (dip-steered, edge-preserving) filtering is first applied to the azimuth-preserving PSTM volume along the T1 reflector. Five edge-sensitive attributes are then extracted on the same 5 m × 5 m grid within a ±20 m window centred on T1: eigenstructure-based coherence (C3); multi-scale most-positive and most-negative curvature; dip deviation from the structure-oriented dip field; the amplitude edge (variance) attribute; and spectral-decomposition energy in the seam-controlling band of 30–50 Hz.

PCA fusion. Let p = 5 attributes be extracted at each of N grid nodes, forming an attribute matrix A of size N × p. To remove scale differences between attributes, each attribute is first standardized (Equation 17):where aᵢⱼ is the j-th attribute at the i-th node, and μⱼ and σⱼ are its mean and standard deviation over the survey area. The covariance matrix of the standardized data is then eigen-decomposed (Equations 18, 19): where C is the covariance matrix and λ_k and e_k are the k-th eigenvalue and eigenvector. The k-th principal-component score and its variance ratio are given by Equation 20:

a ∼ i j = a i j − μ j σ j (17)

C = 1 N − 1 A ∼ T A ∼ (18)

C e k = λ k e k , λ 1 ≥ λ 2 ≥ ⋯ ≥ λ p (19)

P C k = A ∼ e k , η k = λ k Σ i = 1 p λ i (20)

In this area the first three components cumulatively explain 88% of the variance (η₁ ≈ 62%, η₂ ≈ 17%, η₃ ≈ 9%) and are therefore retained (K = 3). The first component loads mainly on the coherence, curvature and edge attributes and dominates the fault response; a variance-weighted sum of the first K components gives the fused fault-enhancement attribute F (Equation 21):

F = Σ k = 1 K w k P C k , w k = λ k Σ i = 1 K λ i (21)

Threshold determination. F is normalized to the interval [0, 1] and the fault-anomaly threshold is set as Equation 22:where μ and σ in T is defined in Equation 22 are the area-wide mean and standard deviation of F, and k is a threshold coefficient. The coefficient is calibrated so that the fused-attribute anomalies coincide with the faults and break points already confirmed on the migrated profiles and at the control boreholes; k = 1.8 gives the best trade-off between detecting the known faults with throw ≥5 m and suppressing false anomalies. Fault traces are then picked where F exceeds T and remains laterally continuous, cross-checked against the coherence and curvature attributes to reject isolated spikes, and finally validated against the borehole-revealed data (Section 4.3).

T = μ F + k σ F (22)

Horizon interpretation and structure-oriented attribute computation were performed in GeoFrame (IESX); three-dimensional visualization was carried out in GeoViz; the multi-attribute fusion interpretation described above was additionally performed in Petrel and VVA; and structural mapping was completed in the CPS-3 geological drafting system.

3 Results

3.1 Processing results

The processing results show that multi-domain noise suppression effectively attenuates strong surface waves and linear interference, enhancing the visibility of effective deep reflections (Figures 13, 14).

Surface-consistent predictive deconvolution improves the lateral consistency of the wavelet, raises the dominant frequency, and stabilizes the phase (Figure 15).

OVT-domain 5D regularization restores the consistency of offset–azimuth sampling (Figure 16).

On this basis, OVT-domain prestack time migration achieves the stable repositioning of reflection–diffraction energy adjacent to complex faults (Figures 17, 18).

Therefore, the theoretical focus of this section is not on changing the integral form of Kirchhoff migration, but on establishing a consistent imaging framework of “OVT-domain regularization → velocity-field-constrained analysis → azimuth-preserving Kirchhoff prestack time migration.” For complex coal mining areas, this geometry–velocity–imaging unified processing approach is more suitable for stable imaging of small faults than locally enhancing any single processing step.

3.2 Quantitative evaluation of imaging quality and method advantages

To verify the effectiveness of the joint processing flow proposed in this study for fine imaging of small faults in complex coal mining areas, the evaluation is conducted at two levels: methodological applicability and imaging results. The former compares the differences between the conventional processing flow and the joint processing flow proposed in this study in terms of data organization, regularization, velocity analysis, and migration imaging; the latter uses indicators with clear statistical bases from the project results to provide a quantitative description.

Table 7 shows that the advantage of the proposed flow does not lie in the local enhancement of any single processing module, but rather in using OVT-domain data organization as a link to couple regularization, velocity control, and azimuth-preserving migration into a continuous processing flow.

TABLE 7. Comparison between the conventional processing flow and the joint processing flow of this study.

According to the project processing results, the final OVT prestack time-migration data volume (Figure 19) is superior to the conventional processing results (Figures 20–23) in terms of SNR, event continuity, and clarity of small-fault break points.

The total length of profiles in the entire area is 181.328 km, of which Class I + II profiles account for 96.258%, and the dominant frequency of migration results is stably above 40 Hz, indicating that the joint processing flow proposed in this study achieves stable imaging in high-density data processing of complex coal mining areas. Table 8 consolidates the core quantitative indicators of the whole chain, from the acquisition geometry through the processing quality to the interpretation output and the borehole verification, together with the corresponding values of the 2006 conventional 3D survey over the same area.

TABLE 8. Key quantitative indicators of the acquisition, processing, interpretation and verification stages, compared with the 2006 conventional 3D survey; the interpretation and verification entries are detailed in Sections 4.3,4.3.1.

The above results indicate that the proposed method achieves a relatively stable level in bandwidth preservation and overall imaging quality. Combining Figures 19, 21–23 shows that target-seam event continuity is enhanced, break-point positions are clearer, geometric distortion of profiles adjacent to obstacles is reduced, and the traceability of small-fault boundaries and local anomaly zones is markedly improved.

To facilitate a comprehensive comparison among different processing flows, the integrated imaging quality metric of this study is defined as Equation 23:where f d o m, S N R, C c o n t, and D f a u l t denote the normalized dominant frequency, normalized SNR, normalized target-seam continuity index, and normalized break-point clarity index, respectively, and α_i are the corresponding weights. This expression is an evaluation function defined in this study for unifying multiple imaging-quality indicators within the same comparison framework, rather than a fundamental geophysical equation.

Q i m g = α 1 f d o m + α 2 S N R + α 3 C c o n t + α 4 D f a u l t , Σ i α i = 1 (23)

Figure 24 summarizes the comparison between the conventional processing flow and the proposed method in terms of bandwidth preservation, SNR, break-point clarity, and processing adaptability. The results show that the proposed method outperforms the conventional flow on multiple indicators, with more stable bandwidth preservation, clearer effective deep reflections, sharper break-point boundaries, and more continuous target-horizon tracking. For fine characterization of small faults in complex coal mining areas, this consistent multi-indicator improvement is more indicative of method stability and engineering applicability than the visual enhancement of any single profile.

In summary, the advantages of the joint processing flow proposed in this study are mainly reflected in three aspects: first, OVT-domain data organization and 5D regularization improve sampling-geometry consistency under complex variable-observation conditions; second, velocity-field-constrained analysis improves the lateral stability of the velocity model in structurally complex regions; third, azimuth-preserving Kirchhoff prestack time migration enhances the focusing of reflection–diffraction information and the sharpening of break-point responses near small faults. Statistical results show that within 181.328 km of profiles, Class I + II profiles account for 96.258%, and the dominant frequency of migration results is stably above 40 Hz, indicating that the proposed method has good imaging stability and can support engineering interpretation in high-density 3D seismic data processing of complex coal mining areas.

4 Discussion

4.1 Mechanism of method effectiveness

Aiming at the requirements for identification of deeply buried, small-throw, small-scale faults in complex coal mining areas, this study establishes an integrated workflow combining high-density sampling, OVT-domain regularization, velocity-field-constrained azimuth-preserving prestack time migration, and structure-oriented interpretation. This workflow corresponds to a complete geophysical pipeline. High-density sampling guarantees the recording of small-scale structural information; OVT-domain regularization restores offset–azimuth sampling consistency; velocity-field constraints improve the stability of the traveltime field; azimuth-preserving migration enhances the repositioning of reflections and diffractions adjacent to faults; and structure-oriented interpretation improves the traceability of fault boundaries and anomaly zones. The results presented above show that this workflow can convert the sampling gain into an imaging gain, and further into an interpretation gain.

4.2 Imaging gain mechanism

The imaging gain in this study originates from the continuous coupling of sampling, regularization, and azimuth-preserving migration. The 5 m × 5 m CDP bin and 300-fold coverage improve the preservation of high-frequency reflections and diffractions adjacent to complex faults. OVT-domain regularization improves sampling-geometry consistency under variable-observation conditions in obstacle areas. Velocity-field-constrained azimuth-preserving prestack time migration improves the repositioning accuracy of break-point locations and boundaries. The results in Section 3 show that the final OVT prestack time-migration data volume outperforms the conventional processing results in SNR, event continuity, and break-point clarity. Within 181.328 km of profiles in the entire area, Class I + II profiles account for 96.258%, with the dominant frequency stably above 40 Hz. These results indicate that high-precision imaging of small faults in complex coal mining areas relies on the unified organization of sampling, regularization, and azimuth-preserving migration.

4.3 Engineering response of interpretation results

The improvement in interpretation accuracy results from the consistency between high-quality migration results and structure-oriented attribute responses. A total of 164 faults were interpreted across the entire area, all of which are normal faults, including 127 faults with throw ≥5 m. Compared with the conventional 3D seismic results from 2006, 89 newly added faults with throw ≥5 m were identified, and 11 previously interpreted faults were rejected. The study area as a whole exhibits a south-dipping monocline, with a NE-trending secondary syncline S1 having an amplitude of about 20 m. One no-coal zone with an area of 0.071 km² and two anomalous zones were also identified within the area. Borehole-revealed verification shows that the agreement rate between interpretation results and borehole revelations reaches 92.7%. These results demonstrate that the proposed method workflow has high external verifiability and can directly support working-face deployment, roadway design, and identification of local risk zones. The survey also satisfied all contractual geological tasks set out in Section 2.1.1: the Class I + II profile share reached 96.258% against a requirement of no less than 95%; 127 faults with throw ≥5 m were interpreted, 89 more than in the 2006 conventional 3D results; and the overall borehole agreement rate of 92.7% exceeds the contractual requirement of more than 80% for faults with throw ≥5 m. Per-throw-class verification rates and the achieved seam-floor depth error are given in the project completion report.

4.3.1 Typical single-well verification

To provide an intuitive basis for the 92.7% overall agreement of the 41 boreholes, Figure 25 shows the well-log curves, the synthetic seismogram and the well-tie migration profile of the typical borehole 1,001. The synthetic seismogram, obtained by convolving the wavelet with the reflection coefficients computed from the sonic (DT) and density (RHOB) logs (Seelam et al., 2024; Lu et al., 2026), matches the actual reflection package on the migration profile, providing an accurate calibration of the T1 reflection at the No. One seam floor and confirming the reliability of the horizon calibration and of the imaging. On this basis all 41 boreholes were compared with the interpretation, and representative wells are listed in Table 9. At boreholes of normal seam thickness the T1 reflection is continuous and strong, in full agreement with the drilling result; at boreholes 846, 914 and 906 the seam is markedly thinned by fault cutting, with drilled thicknesses of 1.1 m, 1.3 m and 2.15 m, and the corresponding T1 events are clearly offset and distorted; and no-seam boreholes such as 953 and 907 lie within fault zones where the events are clearly offset. Except for the three boreholes 838, 874 and 952, which are affected by positional offsets, the remaining 38 boreholes agree fully with the interpretation, giving an agreement rate of 92.7% (Zhou et al., 2024; Tong et al., 2025). These results show that the interpretation reproduces both the seam morphology and the effect of small-throw faults on the seam at the single-well scale.

TABLE 9. Representative single-well seismic verification, selected from the 41 boreholes. For boreholes 3–5 the seam thickness is the fault-reduced value.

4.4 Applicability conditions and research boundaries

The proposed method is applicable to coal mining areas with relatively deep target seams, small-scale faults, dense distributions of obstacles, complex shallow velocity structures, and high requirements on structural precision for mining deployment. The core conditions for method applicability include a high-density observation system, a relatively stable horizon-calibration foundation, and necessary borehole constraints. In imaging-quality evaluation, dominant frequency, profile quality grade, and borehole agreement rate have direct statistical bases, while SNR, continuity index, and break-point clarity index still require a unified area-wide standardized statistical system. Attribute-fusion interpretation has high dependence on migration-result quality and velocity-field stability. For deeper targets, more strongly heterogeneous backgrounds, and more complex wavefield conditions, parameter adjustment and workflow verification still need to be carried out in combination with specific seismic-geological conditions.

This study indicates that the key to fine identification of small faults in complex coal mining areas is to establish a closed technical workflow integrating sampling, regularization, velocity control, imaging, and interpretation. This technical workflow forms a high-precision geophysical interpretation framework oriented to mining engineering decision-making.

5 Conclusion

Aiming at the requirements for fine identification of deeply buried, small-throw, and small-scale faults in the No. One coal seam of the structurally complex Daqiang Coal Mine, this study establishes an integrated technical route comprising high-density sampling, OVT-domain regularization, velocity-field-constrained azimuth-preserving prestack time migration, and structure-oriented multi-attribute interpretation. The main conclusions are as follows.

• This study constructs a complete method workflow for fine identification of small faults in complex coal mining areas. High-density sampling guarantees the recording of small-scale structural information; OVT-domain regularization restores offset–azimuth sampling consistency; azimuth-preserving prestack time migration improves the repositioning quality of reflections and diffractions adjacent to faults; and structure-oriented multi-attribute interpretation enhances the traceability of fault boundaries and anomaly zones. This workflow forms a closed system among sampling, imaging, and interpretation.
• The coupling of high-density sampling and OVT-domain imaging is the main source of the imaging gain in this study. The 5 m × 5 m CDP bin and 300-fold coverage improve the preservation of high-frequency reflections and diffractions adjacent to complex faults. The processing results show that the final OVT prestack time-migration data volume outperforms conventional processing in SNR, event continuity, and break-point clarity. Within 181.328 km of profiles in the entire area, Class I + II profiles account for 96.258%, and the dominant frequency of migration results is stably above 40 Hz.
• Structure-oriented multi-attribute interpretation improves the identification capability of fault systems and anomaly zones in complex coal mining areas. A total of 164 faults were interpreted across the entire area, all normal faults, of which 127 have throw ≥5 m. Compared with conventional 3D seismic results from 2006, 89 new faults with throw ≥5 m were identified, and 11 previously interpreted faults were rejected. The No. One coal seam in the study area generally exhibits a south-dipping monocline, with a NE-trending secondary syncline S1 having an amplitude of about 20 m. One no-coal zone with an area of 0.071 km² and two anomalous zones were also identified within the area.
• Borehole-revealed verification shows that the interpretation results in this study have high reliability. Verification results from 41 boreholes within and around the area show that the agreement rate between interpretation results and borehole revelations reaches 92.7%. This indicates that the workflow established in this study can directly support working-face deployment, roadway design, and identification of local risk zones.
• The proposed method is applicable to coal mining areas with relatively deep target seams, small-scale faults, dense obstacle distributions, complex shallow velocity structures, and high requirements on structural precision for mining deployment. Method applicability for deeper targets, more strongly heterogeneous backgrounds, and more complex wavefield conditions still requires targeted verification combined with specific seismic-geological conditions.

In summary, the integrated technical route established in this study can produce stable, reliable, and engineering-applicable high-precision geophysical interpretation results in complex coal mining areas, providing a verifiable methodological basis for the fine identification of small faults and the deployment of mine production.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Zhixue Sun, China University of Petroleum, China

Reviewed by: Weichao Yan, Ocean University of China, China

Yi Liu, Jilin Jianzhu University, China