Detection of geological anomalies in coal mining working faces using a scattered-wave imaging method

A concealed geological structure encountered during the excavation of a coal working face could connect the working face to high-pressure water in limestone strata, which can result in a serious or catastrophic water inrush accident. However, existing geophysical detection methods used to ensure the geological safety of working faces cannot detect small geological anomalies reliably. Based on the generalized theory of scattered waves, we have developed a novel and superior scattered wave imaging method for the detection at the roadway lateral wall, capable of wave vector extraction and multiwave imaging. In this method, the waves scattered from a geological anomaly can be dynamically and accurately extracted by the polarized filter function during the mapping processes of common scattering point (CSP) gathers. A numerical simulation was performed to study the seismic wave response characteristics of a small collapse column in a coal working face. The P and channel waves of the model were extracted and imaged using the novel imaging method. A field study of three-component seismic detection was performed in the Xuzhuang Coal Mine, demonstrating that the joint imaging of body and channel waves can detect small drop faults invisible to channel wave imaging alone. These results indicate that the proposed method can effectively image anomalous bodies on working faces in complex and noisy mine wavefields using multiwave information, providing a new approach for the reliable and timely detection of hazardous geological features hidden in working faces.


Introduction
According to China's recent "Dual Carbon Goals" policy, carbon dioxide emissions should peak before 2030, and the nation should achieve carbon neutrality by 2060. Because of urgent development needs and the time required to scale up replacement energy sources, coal will remain an essential part of China's energy mix for decades (Xie et al. 2021). Therefore, novel coal mining techniques still need to be developed, and any new safety concerns that 1 3 may arise with this shift in energy sources will need to be addressed.
This study addresses an emerging issue: water inrushes from limestones are increasingly threatening the lower measures of deep coal mines in China. At depths where groundwater pressures are high, the breaching of an unsuspected collapse column can often cause a catastrophic water inrush accident (Wu 2014;Yin et al. 2021). The ability to detect geological limestone structures hidden in the coal face benefits water control, safety, and production continuity in deep mining (Yuan 2020). Two types of methods are available for detecting geological anomalies hidden in working coal faces: exploration drilling and geophysical prospecting. Geophysical prospecting is becoming increasingly popular because of its low cost, quick results, and wide detection range (Liu et al. 2014).
Currently, the most common geophysical prospecting methods for working mine faces are ground-penetrating radar (GPR), DC resistivity, and seismic channel wave (Chen et al. 2019). Although convenient, the application of GPR in mines (Lavu et al. 2011;Emslie and Lagace 2016;Shahrabi and Hashemi 2021) faces limitations such as limited penetration depth, low resolution for strike faults, and serious interference from metallic objects . The DC resistivity technique reveals electrical differences between coal seams and other rock strata. However, the resolution and sensitivity of this method in the early detection of limestone inclusions are often unsatisfactory (Zhang et al. 2015;Hu et al. 2019). Seismic wave detection has proven effective in the high-precision detection of underground unknown geological bodies (Su et al. 2021;Ismail et al. 2021;Gammaldi et al. 2022). When coal seams are sandwiched between other rock strata in coal mines, seismic waves are often trapped in the coal seam, forming easily detectable channel waves (Krey 1963). Channel waves are shaped by reflection at strata boundaries, seam thickness and consequent interference, and other seam parameters; thus, they carry substantial information about the coal seams. The thickness and other parameters of the coal seam can be calculated using the dispersion curves of channel waves (Wang et al. 2012;Dresen and Rüter 2013;Ji et al. 2020). When channel waves encounter geological anomalies in coal seams, they reflect or scatter, and this property can be used to detect such structures (Buchanan et al. 1981;Pant et al. 1992;Teng et al. 2019;Guo et al. 2020). Therefore, channel wave exploration has become widely used for detecting geological anomalies in the working faces of coal seams Hu et al. 2022). However, channel waves cannot detect small geological structures such as faults with less than half the coal thickness and small diameter collapse columns; these small anomalies can render deeper anomalies undetectable by causing interference. Therefore, channel wave analysis cannot adequately detect geological anomalies in working faces .
In recent years, the analysis of diffracted and scattered waves has become an emerging trend in geophysics owing to its potential for imaging at higher resolutions than channel waves (Zhao et al. 2020;Xi and Huang 2020). The analysis of diffracted and scattered waves can effectively identify discontinuous geological bodies (Wu and Aki 2012;Schmandt et al. 2019;Mangalampalli et al. 2022) to a degree suitable for the detection of small structures in coal mines (Hu et al. 2010;Liu et al. 2022). However, seismic waves of different types and directions are mixed at a coal face. A key challenge in developing this method is effectively identifying the signatures of geological anomalies amid complex and mixed signals .
Based on the generalized theory of scattered waves and the practice of detecting small geological anomalies behind the roadway lateral wall of working faces, we developed a new imaging method that can achieve vector extraction and multiwave imaging of scattered waves to detect geological anomalies behind roadway lateral walls in the working faces of coal mines. The results of seismic numerical simulations and field studies conducted to validate the new method are also presented. Our results show that the novel imaging method effectively detects hidden geological anomalies in working faces.

Equivalent offset migration in the working face
Based on the Huygens-Fresnel principle, the equivalent offset migration (EOM) method was first proposed in 1994 (Bancroft and Geiger 1994) to use the CSP gather formed by the equivalent offset for prestack migration to realize the imaging of scattered waves (Bancroft et al. 1998). The method was later improved by measuring the EOM of converted waves (Wang et al. 1999), rugged topography (Geiger and Bancroft, 1996), and vertical receiver arrays (Bancroft and Xu 1999). Furthermore, the application of this method in metal mines (Gou et al. 2007) and urban metro projects ) has been successful with good detection. A normal two-dimensional vertical detecting direction can be changed into a bedding horizontal detecting direction for in-plane abnormal geological structure detection performed at the lateral wall of the roadway in a mining face. Based on the aforementioned concept, the signal of a scatter point can be obtained at each receiving point when a seismic wave is excited on the lateral wall. In the EOM method (Bancroft and Geiger 1994;Bancroft et al. 1998), the double-squareroot equation of the time t of the seismic wave, propagating from a source point S at the roadway lateral wall to a scattering point SP at depth Z 0 and then to the receiver point R at the roadway lateral wall, can be converted into a single square root equation by assuming the existence of a selfexcited self-receiving point E on the roadway lateral wall (Fig. 1). For convenient calculation, we define the midpoint of S and R as M. At this time, the equivalent offset h e , which is the distance from E to the lateral wall projection point CSP of the scattering point SP in the working face, is given as where h is the distance between S and M; x is the distance from M to CSP; and v mig represents the root-mean-square velocity at the scattering point.
The scattered waves are mapped from the input traces to the CSP gathers using Eq. 1 to extract the scattered waves to the maximum theoretical extent without any loss of reflected waves. Finally, after obtaining the velocity analysis results for the CSP gather along the hyperbolic trajectory, we use Kirchhoff migration on the gathers to complete the scattered wave imaging. Compared with the conventional reflected wave method, the EOM method, with its strong extraction capability and high folds of scattered waves, is superior in the detection of small-sized structures in coal mines.

Polarized filter function
Adaptive instantaneous polarization analysis (Diallo et al. 2006), a Hilbert transform-based method, accurately selects a dynamic covariance calculation time window using the instantaneous frequency characteristics of a signal and then obtains the instantaneous polarization parameters of each sampling point in three-component (X, Y, Z) seismic records. The procedure for solving the covariance matrix is described in detail in the cited article and will not be repeated here.
Below is a list of the relevant equations of the polarization parameters when using the adaptive instantaneous polarization analysis method.
( In this formula, the range of T e is (0,1). When T e is equal to 1, the signal shows characteristics of complete linear polarization; when T e is equal to 0, it shows characteristics of ellipsoidal polarization. (3) M a j o r p o l a r i z a t i o n a z i m u t h a n g l e = arctan ℜ y 1 (t) ℜ x 1 (t) . Here, ℜ represents the real part of the signal, and x 1 (t) and y 1 (t) represent the seismic records of two horizontal components in a single receiver.
(4) M a j o r p o l a r i z a t i o n d i p a n g l e Here, z 1 (t) represents the seismic records of a vertical component in a single receiver. (5) Major polarization direction �� ⃗ V . The eigenvector corresponding to eigenvalue λ 1 is the major polarization direction.
Working faces contain different types of seismic waves in different directions. To accurately image the scattered seismic waves in working faces, extracting specific types of seismic signals at specific locations is necessary. To this end, a polarization filtering function with a spatial direction filtering attribute is constructed using Eq. 2: For this function, G 1 (t) = T e (t) or 1−T e (t) is the polariza- is the directional filter factor. In the equation of G 2 (t), i is a placeholder for the subscript corresponding to each type of objective imaging wave, such as p for the P wave, sh for the SH wave, and sv for the SV wave, among others. The superscripts a and b represent the exponentially weighted values of the two factors, which adjust the system bias for different filtering purposes. The polarization strength factor G 1 is mainly used to suppress nonlinear waves such as noise signals. The parameter a is usually assigned the value of 1, and increasing the value of a can further suppress nonlinear wave signals.
(2)  Moreover, G 1 is represented by 1 − T e when extracting nonlinear waves, such as Rayleigh waves. The directional filter factor G 2 filters the waves based on direction, using the angle between the theoretical direction of vibration and the main polarization direction of the objective wave to complete the dynamic selection and filtering of seismic waves of different types and directions. The parameter b is typically set to 1, and increasing the value of b can further suppress objectivetype echoes in non-objective directions. For a known scattering point in space, the propagation direction vector of the seismic echo signal at a receiving point is � ⃗ L , and the major polarization direction of the three-component seismic record is �� ⃗ V . In the formula for the directional filter factor, �� ⃗ L i represents the theoretical vibration direction vector of the different types of seismic waves, where the P-wave vector ��� ⃗ L p is consistent with � ⃗ L . The SHwave vector ���� ⃗ L sh and the SV-shear-wave vector ���� ⃗ L sv are perpendicular to � ⃗ L in the shot-receiver plane and the normal line plane, respectively. The Love channel wave vector �� ⃗ L l is consistent with ���� ⃗ L sh . It is evident that the filter function can reduce the influence of nonlinear interference waves and improve the signal-to-noise ratio (SNR) of seismic records. In addition, the manner in which its filtering characteristics change with the spatial position of the imaging points and the type of imaging wave can improve the filtering accuracy and, to a certain extent, enable vector extraction of the wavefield.

Scattered-wave imaging method
Figure 2 depicts an observational arrangement for seismic scattered wave detection, consistent with the usual deployment of sensors for working faces. For a horizontal working face, a survey line with multiple source and receiver points is arranged along the lateral wall of the roadway, where X is the roadway strike direction of the working face, Y is the bedding direction of the working face, and the scattered wave imaging region with multiple CSP gathers (CSP1 to CSPn) is in the XY plane. A two-dimensional vector in the XY plane can now express the propagation direction and theoretical polarization direction of seismic waves from any scattering point (SP) to the receiver. Figure 2 shows the theoretical polarization directions as green lines.
Using the aforementioned observational setup, the equivalent offset formula and polarization filter function for detection at roadway lateral walls are dynamically combined to map the CSP gathers in the detection region, allowing for the dynamic extraction of scattered waves of different types and directions. Further, the velocity analysis and Kirchhoff migration of the CSP gathers of the different types of scattered waves are performed to enable detection at roadway lateral walls. Note that the velocity is usually unknown. We set the initial velocity v 0 (usually referring to the velocity of direct waves) to map the CSP gathers and perform velocity analysis. Then, the multiple iterations of mapping CSP gathers are performed using the values of velocity analysis. In the case of uncomplicated formation, a stable velocity profile can be obtained via two to three iterations (Bancroft et al. 1998), which not only provides accurate velocity information for migration imaging but also improves the accuracy of mapping the CSP gathers. Each gather outputs a stacked trace using Kirchhoff prestack migration based on the velocity profile. Finally, all the stacked traces are combined to form a seismic profile.
A basic flowchart of this method is shown in Fig. 3, wherein the amplitude of the main polarization direction A(t) represents the composite value of the three-component signals in the main polarization direction, and the output value A out (t) represents the amplitude of the objective imaging wave after polarization filtering. Defining the included angles between the main polarization direction and the three axes of the receiver, X, Y, and Z, as α, β, and γ, respectively, the amplitude of the main polarization direction A(t) can be calculated using Eq. 3: Here, A x (t), A y (t), and A z (t) represent the amplitude values of seismic signals in the X, Y, and Z components, respectively, at time t.

Model parameters
The geological structure of a coal mine working face was simplified to a three-layer structure of rock-coal-rock. A three-dimensional (3D) model of a working face with a (3)  Figure 4c shows a cross sectional view of the XY plane at Z = 50 m, demonstrating that the collapse column has a diameter of 10 m and that the midpoint is located at point (100, 155). To obtain realistic petrophysical parameters for the working face model, we used the lithology of common coal measure strata, as shown in Table 1.   The seismic survey line was set up near the floor on the lateral wall of the working face and extended along the wall in the X direction. A total of 21 three-component receivers at 10 m intervals and 21 sources at 10 m intervals were installed to achieve wave acquisition. A mesh size of 0.5 × 0.5 × 0.5 m, a Ricker wavelet of 200 Hz, point source excitation, and a sampling interval of 0.05 ms were used for numerical calculations, with a recording length of 350 ms. A high-order, staggered-grid, finite difference scheme was used for elastic wave forward modeling (Xia et al. 2004). Figure 5 shows the combined snapshots of the 11th shot signal of the 3D working face model, comprising one YZ section (X = 100 m), one XZ section (Y = 100 m), and two XY sections (Z = 50 m and 100 m). The YZ section of Fig. 5 clearly shows that the seam acts as a waveguide for seismic waves and that the signals entering and exiting the coal seam are different. Three types of waves, P, S, and channel waves, exist on the XY-section at Z = 50 m, while only body waves exist on the XY-section of non-coal seam rock at Z = 100 m. It is also apparent that the presence of the coal seam alters the body waves, with the P wave train becoming longer and the S and channel waves overlapping in the coal seam. Figure 6 shows XY-section snapshots of three components of the simulated wavefield at Z = 50 m at different times. In the 60 ms snapshots (Fig. 6a, b, and c), the incident P wave scatters at the position of the collapse column. The P wave is stronger in the Y component but much weaker than the channel wave in general. Under the influence of the coal seam, the P wave has a long wave train. In the 110 ms snapshots (Fig. 6d, e, and f), channel waves with different phase velocities are scattered at the collapse column in ways that show stronger energy in the Y and Z components. The reason for these phenomena is that the incident waves generated by the point source are mainly P and converted SV waves, which would form Rayleigh channel waves under the waveguide action of the coal seam (Yang et al. 2010). Rayleigh channel waves are elliptically polarized waves that vibrate primarily in the Y and Z directions. Figure 7 shows the synthetic three-component records of the model in the 11th shot. Figure 7a shows weak scattered wave signals in the X component, indicating that the main component of the signal consists of seismic waves that vibrate in non-X directions. As shown in Fig. 7b, the Y component of the P waves scattered from the collapse column is the most visible of the three components. The P wave train scattered from the distal face of the collapse column is longer than that scattered from the proximal face owing to multiple scattering effects within the collapse column. The Y component of the column channel wave scattered from the collapse column is also strong. As shown in Fig. 7c, the scattered P wave in the Z component is not visible, whereas the channel wave in the Z component is visible. If the S wave signal is not visible in all three axes, we may regard the P and channel waves as the target waves for vector separation of the wavefields, completing the imaging of small anomalies in the complex wavefields of coal seams.

Polarization analysis
For convenient calculation, the synthetic seismic records were resampled to reduce the sampling interval to 0.1 ms. As shown in Fig. 8, the adaptive instantaneous polarization analysis method was used to analyze the threecomponent seismic response records of 300 sampling points in the 94-124 and 170-200 ms intervals of the 11th receiver in the 11th shot. Based on the polarization coefficient results shown in Fig. 8a, the polarization coefficient of the 94-124 ms interval (the black line) is significantly greater than that of the 170-200 ms interval (the red line), which indicates that the signal in the earlier interval is more linearly polarized and conforms to the polarization characteristics of P waves. However, owing to the interference from the coal seams, the polarization linearity of the signals in the 94-124 ms interval has local jitter, implying that it is not as stable as the theoretical signal. The polarization coefficient of the 170-200 ms signal, concentrated near the value of 0.2, displays a degree of ellipsoidal polarization, conforming to the polarization characteristic of the Rayleigh channel waves. Clearly, the polarization strength factor G 1 of the polarization filter function developed can effectively separate the two types of waves. Figure 8b shows the polarization dip angle of the signals from the two intervals, and Fig. 9 presents the polarization azimuth of the signals from the two intervals. Combining the two intervals for comparative analysis, the following is observed: (1) the signal in the 94-124 ms interval is the P wave generated right above the collapse column, with polarization dip angles near 0° and polarization azimuth angles near 90°; These results show that the directional filter factor G 2 of the polarization filter function helps extract the objective imaging wave in the objective direction. Thus, the adaptive instantaneous polarization analysis method accurately obtains the dynamic polarization parameters of objective imaging waves in complex signals. The vector wavefield separation of the P and Rayleigh channel waves can be realized using the polarization filter function.

Imaging results
The CSP gathers of the model data were mapped using the conventional EOM method and the lateral wall scattered wave imaging method described in this paper. A total of 41 CSP gathers were arranged with a spacing of 5 m, a maximum equivalent offset distance of 100 m, an equivalent offset step of 5 m, and a root-mean-square velocity v mig equal to the P wave velocity of the model. Figure 10 displays the No. 21 CSP gather (at X = 100 m) mapped using the two methods. Figure 10a shows that there are many types of waves scattered from the collapse column in the CSP gather, mapped by the conventional EOM method, in which  Figure 10b, c, and d show the three single wave CSP gathers of the P, Rayleigh channel, and Love channel waves mapped by the processing steps shown in Fig. 3. When extracting linearly polarized P and Love channel waves, the polarization strength factor G 1 is derived using the value of T e . When extracting ellipsoidally polarized Rayleigh channel waves, G 1 is derived using the value of 1 − T e . The directional filter factor G 2 is selected according to the difference in the theoretical vibration directions of the three waves (refer to Sect. "Polarized filter function" for details). Compared with Fig. 10a, the P-wave CSP gather shown in Fig. 10b only reveals the P waves scattered at the collapse column, while the other waves are wholly suppressed. The Rayleigh channel wave CSP gather shown in Fig. 10c only reveals the Rayleigh channel waves scattered at the collapse column, while the other types of waves are entirely suppressed. The Love channel wave CSP gather shown in Fig. 10d reveals no significant signal because the Love channel wave is not developed in the simulation.
Based on the aforementioned observations, the conventional EOM method is unable to effectively map the CSP gathers of the objective single wave, producing an imaging error if that method is used. However, this issue does not apply to the novel imaging method, which is capable of accurate imaging.
In accordance with the time-distance relationship of the seismic waves in the CSP gathers shown in Fig. 10, the velocity analysis was performed using stacked energies with a velocity interval of 0.1 m/ms and a velocity scanning range of 1.0-4.5 m/ms. Figure 11 shows the velocity analysis spectrum of the P-and channel wave CSP gathers after separating the direct waves. The P-wave velocity is ~ 3.2 m/ms, and the channel wave velocity is ~ 1.7 m/ms. The P-wave velocity is consistent with the roof and floor of the model, and the channel wave velocity is lower than the S-wave velocity.
Following the velocity analysis results, migration imaging was performed on the CSP gathers of P and Rayleigh channel waves, and the results are shown in Fig. 12. Figure 12a  and b present the P and Rayleigh channel waves, respectively, and the green circle represents the actual position of the collapse column. The scattered wave imaging method for the working face can be used to image different types of waves from the roadway lateral wall. The imaging of the two wave types can suppress the circular arc phenomena that easily appear in conventional reflection wave imaging. The imaging of the proximal face of the collapse column is consistent with the model, and the reason for the deviation in the imaging of the distal face depth is that a low-velocity medium exists inside the collapse column. Although the deviation in the imaging of the distal face affects the judgment of the collapse column diameter, accurate imaging of the proximal face shows a great value for the correct prediction of the nearest position of abnormal bodies in practical detections.

Field overview
The primary mining coal seam in the Xuzhuang Coal Mine, Jiangsu Province, is the No. 7 coal seam. It has an average thickness of 5 m and is located in the middle and lower part of the Shanxi formation. As shown in Fig. 13, the unformed working face 7335 is situated in the deep part of the mining area. In this area, there is a suspected fault labeled F117, with a drop of 15 m (inferred by 3D seismic exploration in the ground). If the unknown strike of the fault crosses into the designed region, it would significantly impact the layout of the working face. In addition, the fault named f 38 72 , with a drop of 1.5 m exposed in another roadway, may also affect the layout. Therefore, determining the distribution of the faults in advance is necessary to guide the subsequent safe and efficient construction of working face 7335.

Data acquisition
As shown in Fig. 14, a seismic survey line was arranged along the lateral wall of part of roadway 7335 according to the detection requirements of the numerical model. The X direction represented the roadway strike direction, and the Y direction represented the detection direction. A total of 21 sources at 10 m intervals and 21 three-component down-hole geophones at 10 m intervals were installed along a 200 m survey line. The X and Y directions of the three-component geophones were parallel to and horizontally perpendicular to the survey line, respectively, whereas the Z direction was vertically perpendicular to the survey line. Explosive charges containing 100 g of explosive were placed into blast holes at depths of 2 m to serve as field sources. To ensure effective excitation, the blast holes were blocked with mud. Figure 15 shows the original three-component records of a single shot in the field under an automatic gain control (AGC) display. The abandoned traces were removed to facilitate analysis, and the record was arranged in the order of the X, Y, and Z components. The SNR of the physically observed wave field in the mine is low, and the component distribution of the different wave types is clearly different. The linear arrival direct waves (the magenta line in the figure) show three forms: P, S, and channel waves. Comparing travel distances and arrival times, we found that the apparent velocities of the direct P wave, direct S wave, and direct channel wave were ~ 3.7−4.0, 2−2.2, and 0.9−1.1 m/ms, respectively, consistent with the geological background. The roadway acoustic wave exists mainly in the X and Z components owing to its rolling propagation along the roadway strike. The hyperbolic weakscatterer wave in the Y component (indicated by the blue lines), which cannot be seen in the X and Z components, is inferred to be P-wave signals scattered from the detection area. The hyperbolic channel wave in the Y component (indicated by the red line) is the channel wave scattered from a geological anomaly in the detection region; however, it displays poor continuity and is only visible in a portion of the traces. Hence, effective detection requires fully extracting the relevant signals from the wavefields of working coal mines, which are more complex and display a lower SNR than our numerical simulation. As shown in Fig. 16, after filtering and denoising the acquired data, vector extraction and imaging of the P waves and Rayleigh channel waves were performed using the novel imaging method in the context of the numerical simulation. A total of 41 CSP gathers were arranged with an interval of 5 m, a maximum equivalent offset distance of 100 m, an equivalent offset step of 5 m, and rootmean-square velocities v mig of 3.8 m/ms for the P waves and 1.0 m/ms for the channel waves, based on empirically observed direct wave velocities. The parameters of the polarization filtering function were the same as in the numerical simulation. Figures 16a and b show the depth profiles of the P and Rayleigh channel waves, respectively. The X direction represents the roadway strike direction, and the Y direction represents the detection direction and depth of observation.

Results and analysis
The following observations were made: (1) in Fig. 16a, owing to the low SNR of the original records and the lack of clarity in the record of the scattered P waves, the phase axis in the imaging section is not sufficiently continuous. These results imply that there are two strong abnormal interfaces in the section. The one located at a detection depth of 220-230 m (the red line in the figure) has a relatively strong continuous phase axis and is inferred to be the fault F117; the other, located at a detection depth of 130-150 m (the magenta ellipse in the figure) has a relatively weak continuous phase axis and is inferred to be fault f 38 72 with a small drop.
(2) An envelope processing of the mapped CSP gather was performed in the channel wave imaging, as shown in Fig. 16b. This can reduce the problem of superposition difference caused by the phase inconsistency of the channel wave (Ge et al. 2008). The channel wave section also has an abnormal interface at 220-240 m, which is inferred to be the F117 fault. The absence of an interface in the middle of the channel wave section is because channel waves were poorly developed there. Comparing the two images revealed that the channel wave cannot respond to a fault with a drop of less than half the coal thickness and that the imaging of scattered body waves can effectively compensate for this defect. The novel imaging method combines the detection of body and channel waves and improves the accuracy of geological structure detection in working faces.

Discussion
According to the Huygens-Fresnel principle, any change of seismic waves caused by the inhomogeneity of a 3D geological space can be called a scattered seismic wave (Wu and Aki 1988), a term that broadly includes direct waves, reflected waves, refracted waves, diffracted waves, among others. The scattered wave has the advantages of high accuracy and a low requirement in SNR, rendering it highly suitable for detecting small-scale geological anomalies in the complex wave field of coal mines. Therefore, conducting corresponding research on the mine scattered wave imaging method is necessary to improve the accuracy of geological structure detection in working faces.
The EOM method is a complete and applied scatteredwave migration imaging method proposed by Bancroft (Bancroft and Geiger 1994;Bancroft et al. 1998). A detailed derivation and method can be found in the references given in this article. Under the full space effect of the coal mine, different types of seismic waves in different directions are mixed. The detections at the roadway lateral wall of working faces focus on geological anomalies in the same strata, in which the seismic wave propagation path is simple, and the polarization angle of the three-component signal strongly correlates with the positions of the anomalies. Therefore, we constructed the instantaneous polarization filter function of the three-component signal and integrated it into the CSP gather mapping process to accurately and vectorially extract scattered waves of objective types and avoid the false appearance caused by the coexistence of different types of waves.
Although the instantaneous polarization analysis method adopted herein can accurately obtain the polarization parameters at each point, there may be a sudden change or overflow of the instantaneous frequency value given complex signals, which may significantly affect the extraction effect of the target wave. Therefore, in practical applications, manmade constraints are required to reduce the effect of the abnormal values on the final result according to the possible signal polarization characteristics of the detection target.
The polarized filter function has the dual functions of linearity and directional filtering. The polarization strength factor G 1 , represented by T e or 1 − T e , can suppress nonlinear waves (such as acoustic waves) and extract elliptical polarization waves (such as Rayleigh waves). The directional filter factor G 2 uses the cosine value of the angle between the theoretical direction of vibration and the main polarization direction of the objective wave to ensure filtering stability. The parameters a and b are usually assigned the value 1 to prevent the filtering from changing the authenticity of the signal. When facing a signal with a low SNR, the values of a and b can be increased to enhance the filtering effect. Most existing seismic detection methods arranged along the roadway lateral wall of working coal faces capture data from channel waves in coal seams and have achieved remarkable results . However, in practice, channel waves may not be fully developed because of differences in the structures and petrophysical properties of coal seams. Channel waves do not respond to small geological structures and hinder the use of body waves to detect deep geological anomalies. Therefore, extracting body and channel waves separately for joint detection significantly improves detection accuracy in coal mines. This was the goal behind developing the multiwave imaging method.

Conclusions
Geological anomalies hidden in the working faces of coal mines constitute major sources of danger. A novel imaging method capable of vector extraction and multiwave imaging of scattered waves was developed for detecting geological anomalies behind roadway lateral walls in the working faces of coal mines, yielding the following results.
(1) This novel method can dynamically and accurately extract the scattered waves of objective types in the target area using the polarized filter function in the CSP gather mapping process. (2) In the numerical simulation, the waveguide effect of coal seams would generate channel waves and cause changes in the length of the wave train. The P and Rayleigh channel waves are the wave types with structural response abilities, and both can be effectively extracted and imaged for a small collapse column. (3) In the field study, the joint detection of body waves and channel waves can account for the defect caused by the inability of the channel wave to respond to a fault with a drop of less than half the coal thickness. This result shows that this novel imaging method can improve the accuracy of detecting geological anomalies in working faces.
This study is the first to systematically apply the scattered wave imaging method to detect hidden anomalies in coal mining working faces. This method can effectively realize the imaging of different types of scattered waves generated at anomalies in working faces under complex mine wavefield conditions and provide a new approach for the reliable and timely detection of hazardous geological features hidden in working faces.

Conflict of interest
On behalf of all authors, the corresponding author states that there are no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.