*3.1. Attribute-Based Predicted Multiparameter Analysis for Goaf Prediction*

The geophysical characteristics of the goaf and its upper strata have changed compared with those before mining, including differences in density, velocity and other elastic parameters, which provide the physical basis for detecting goafs by the seismic method [24]. A large amount of data proves that the existence of goafs changes the characteristics of coal-seam reflection waves, which is manifested in the distortion, dislocation, and even disappearance of coal-seam reflection-wave groups, the reduction in seismic amplitude, frequency, phase reversal, and waveform changes. These characteristics are the main signs for identifying goafs. Currently, seismic attributes are among the most effective ways to identify goafs.

The extraction of different seismic attributes is essentially the selection of different attribute parameters. Different geological characteristics correspond to different physical parameters. Therefore, selecting reasonable seismic attributes in a wide range of seismic attributes depends on the characteristics of geological bodies or geological phenomena. To improve the accuracy of identifying goafs, we analyzed the geological characteristics of goaf development and its seismic response and selected several layer attributes that are extremely sensitive to goafs and small faults from more than 20 attributes, including average amplitude energy, instantaneous frequency, dominant frequency, instantaneous phase, and waveform-clustering attributes.

Vakman compared various methods for extracting instantaneous parameters of signals and concluded that the analytical signal method is the most reasonable [25]. The analytical signal method uses the Hilbert transform to obtain the corresponding analytical signal to form a complex signal and then extracts the instantaneous parameters of the signal.

Assuming that the real signal is *x*(*t*), the corresponding imaginary part can be obtained by the Hilbert transform [26]:

$$y(t) = HT[\mathbf{x}(t)] = \frac{1}{\pi} \int\_{-\infty}^{\infty} \frac{\mathbf{x}(\tau)}{t - \tau} d\tau \tag{1}$$

*HT*[•] represents the Hilbert transform. Then, we can construct an analytic signal *s*(*t*) by *x*(*t*) and *y*(*t*):

$$s(t) = x(t) + iy(t) = A(t)e^{i\theta(t)}\tag{2}$$

According to the analytic signal, the instantaneous parameters of the original signal *x*(*t*) can be extracted, including

$$A(t) = \sqrt{\mathbf{x}^2(t) + y^2(t)}\tag{3}$$

$$\text{plhi}(t) = \tan^{-1}[y(t)/x(t)]\tag{4}$$

$$
\omega(t) = \frac{1}{2\pi} \frac{dphi(t)}{dt} = \frac{1}{2\pi} \frac{\mathbf{x}(t)y'(t) - y(t)\mathbf{x}'(t)}{\mathbf{x}^2(t) + y^2(t)}\tag{5}
$$

where *A*(*t*), *φ*(*t*), and *ω*(*t*) represent the instantaneous amplitude, instantaneous phase, and instantaneous frequency, respectively.

The instantaneous amplitude is a measure of the reflective intensity, which is the square root of the total energy of the real and imaginary parts of the analytical signal. It mainly reflects the change in energy and can highlight the change in special rock strata. The instantaneous phase is a measure of the continuity of seismic events in seismic profiles, and its phase can be displayed regardless of the intensity of energy, even if the weak-amplitude effective wave can be well displayed on the instantaneous phase diagram. When the wave propagates in an anisotropic homogeneous medium, its phase is continuous when the wave is in an abnormal medium. Its phase will change significantly in the abnormal position, and it is obviously discontinuous in the profile. The instantaneous frequency is the rate of time variation of the phase, which can reflect the lithological changes in the stratum and help to identify the stratum.

For the same detection object, the obvious change in three instantaneous pieces of information in the same location may reflect a change in the physical properties of the detection object. Among the three parameters, the instantaneous phase spectrum has the highest resolution. The change in the instantaneous frequency and amplitude spectrum is also more intuitive and can be used to determine the general location of underground anomalies and the instantaneous phase to determine the stratified boundary.

The C3 algorithm is a coherence algorithm based on feature structure analysis [27]. The calculation of coherence values in the algorithm can use multichannel seismic records. When calculating the coherence value at the sampling point (*x*, *y*, *t*), first define a data body containing the point (controlled by a two-dimensional spatial window and a onedimensional time window), assuming that the time window contains J channels. For the seismic record, the number of sampling points in the time window is 2 N + 1; let *dT <sup>n</sup>* (*x*, *y*, *t*)=(*dn*1, *dn*2, ··· , *dn J*), *n* = −*N*, −*N* + 1, ··· , *N* − 1, *N*. The seismic trace in the data body can be represented as matrix.

$$\mathbf{D}(x, y, t) = \begin{bmatrix} d\_{-N\prime} d\_{-N+1\prime} \cdots d\_{\prime} d\_{N} \end{bmatrix}^T \tag{6}$$

The correlation matrix is calculated as follows:

$$\mathbf{C}(\mathbf{x},y,t) = \mathbf{D}^{T}\mathbf{D} = \sum\_{n=-N}^{n=N} d\_{n}d\_{n}^{T} = \begin{pmatrix} \sum\_{n=-N}^{n=N} d\_{n1}^{2} & \sum\_{n=-N}^{n=N} d\_{n1}d\_{n2} & \dots & \sum\_{n=-N}^{n=N} d\_{n1}d\_{nf} \\ \sum\_{n=-N}^{n=N} d\_{n2}d\_{n1} & \sum\_{n=-N}^{n=N} d\_{n2}^{2} & \dots & \sum\_{n=-N}^{n=N} d\_{n2}d\_{nf} \\ \vdots & \vdots & \ddots & \vdots \\ \sum\_{n=-N}^{n=N} d\_{nf}d\_{n1} & \sum\_{n=-N}^{n=N} d\_{nf}d\_{n2} & \dots & \sum\_{n=-N}^{n=N} d\_{nf}^{2} \end{pmatrix} \tag{7}$$

The coherence value *E*c (*x, y, t*) at point (*x, y, t*) is defined as

$$E\_{\mathsf{c}}(\mathbf{x}, y, t) = \frac{\lambda\_{\max}}{\operatorname{Tr}(\mathsf{C})} = \frac{\lambda\_{\max}}{\sum\_{j=1}^{J} \mathsf{C}\_{jj}} = \frac{\lambda\_{\max}}{\sum\_{j=1}^{J} \lambda\_{j}} \tag{8}$$

where *λmax* is the maximum value in *λj*.

It is easy to know from Equation (7) that the correlation matrix is positive definite or semipositive definite, so the value of *Ec* falls in the range of 1/J~J. When the waveforms of different traces are the same, *Ec* is equal to 1; when the waveforms are different, *Ec* is less than 1, and the greater the difference between waveforms, the smaller the value of Ec. The difference between waveforms can reflect the changes in geological bodies, so *Ec* can depict the discontinuities of geological bodies.

#### *3.2. TEM-Based Apparent Resistivity Prediction for Water Accumulation Identification*

The transient electromagnetic method (TEM) is a geophysical method that uses the principle of electromagnetic induction [28,29]. It emits pulsed electromagnetic waves as an excitation field source (called 'primary field') through an ungrounded return line or ground-line-source electric dipole. It also observes and studies the spatial distribution and time characteristics of the electromagnetic field (called 'secondary field') of the stratum or geological target body (underground goaf) generated by the excitation field (called 'primary field'), which can be used to infer and explain the spatial structure and physical characteristics of the stratum or goaf water accumulation.

Transient electromagnetic measurements are carried out by using the secondary eddycurrent field excited by a step-pulse current in a large fixed source loop [30,31]. A rectangular emitter loop is placed on the surface of a homogeneous isotropic ground with magnetic conductivity and conductivity. The area of the loop is S, and the inner loop provides a step-pulse current. Before disconnecting the current, the emission current produces a stable magnetic field in the space around the loop and in the ground, as shown in Figure 3a.

When the current is turned off, the magnetic field generated by the emission current also disappears immediately. This change is quickly transmitted to the ground around the return line in the primary field with the aid of air and underground conductive medium and induces a current in the earth to maintain the magnetic field before the emission current is disconnected, so that the magnetic field of the space does not immediately disappear. However, because of the ohmic loss of the medium, the magnetic field generated by it will decay rapidly with the decay of the induced current. The rapid decay of the magnetic field

will induce a new weaker eddy current in the surrounding medium. This process continues until the earth's ohmic losses have exhausted the magnetic field energy (see Figure 3b), forming a geodetic electromagnetic field as the target of detection.

**Figure 3.** Rectangular loop lines and equivalent eddy-current ring, (**a**) establishment of stable magnetic field; (**b**) attenuation of magnetic field (Tx: transmission line frame. t: decay time).

#### **4. Application and Discussion**

The high-density three-dimensional seismic and transient electromagnetic combined detection method is employed to identify the goafs in the detection area. Three-dimensional seismic design has 5843 survey points in the whole area, including 2399 shot points and 3444 detection points. The exploration area is 0.25 km2, and the coverage area of 24 times is 0.2 km2. The normal distance of the transient electromagnetic method is 26~35 m. The engineering survey line is arranged along the design roadway and its surroundings. The point distance of key areas above the roadway is 10 m, some areas are densified to 5 m, and the point distance of other areas is 20 m. The survey lines are arranged along the east–west axis, and a total of six survey lines are arranged. Two survey lines are arranged above the No. 4 coal roadway, numbered N0 and N35 from south to north. Four survey lines are arranged along the No. 9 coal roadway, numbered l0, L33, l59, and L85 from south to north, and the survey points are numbered 0, 20, and 40 from west to east according to the point distance, as shown in Figure 4.

**Figure 4.** Layout of 3D seismic and transient electrical detection.

The buried depth of Coal Seam 4 is 0–150 m, with an average of 70 m. The 9th coal seam has a buried depth of 40–210 m, with an average of about 110 m. The extremely shallow burial depth of the target layer and the small survey area led to the limited length of survey lines. Therefore, seismic exploration adopts the method of small array and high-density acquisition, with high shot-point density and high coverage times. After the field collection test, the shot interval has changed from 24 m to 6 m, and the coverage times have changed from 24 to 64. Meanwhile, due to the restriction of the field acquisition and construction conditions (vehicle walking platform), the distribution of shot and detection

points is extremely uneven (Figure 5a), resulting in extremely uneven coverage of the study area (Figure 5b).

**Figure 5.** The distribution of shot and receiver points (**a**), the actual coverage times (**b**), and coverage times after bin homogenization (**c**) in the study area.

Due to the limitation of terrain conditions, such as large drop and steep slopes between the platforms, the original coverage times are extremely uneven. Therefore, the surface homogenization method is used to eliminate the adverse effects of uneven coverage, including static bin homogenization and dynamic bin homogenization. Static bin homogenization is used to make only one track in each offset group in a bin. If there are N tracks in an offset group in this bin, the N tracks in the offset are weighted and superimposed to represent the offset group; if there are N tracks in an offset group in the bin, then the N tracks are weighted and superimposed, and the superimposed track represents the offset group. Dynamic surface equalization is mainly used to find the offset group in the adjacent element. Define a large rectangular bin centered on the current bin. Any missing offset set in the current bin searches for the corresponding offset set in the adjacent bin within a circle, where the radius of the circle is a linear function of the offset set and its range is limited to the large bin.

During the processing, the seismic trace is edited on the single-shot record to eliminate the strong interference channel and the invalid redundant channel. On this basis, dynamic panel homogenization is used to supplement the missing offset set, and static panel homogenization is used to make the coverage times and energy more balanced (Figure 5c).

Figure 6 shows the CDP gathers before and after bin homogenization, where Figure 6a,b show the processing effects of bin-homogenizing the number of insufficient superpositions traced to the average number of superpositions. Figure 6c, d show the processing effects of bin-homogenizing the number of excessive superpositions traced to the average number of superpositions. From the processing effects, we find that the basic shape of the seismic profile has not changed. The overall energy of the seismic section is more balanced.

**Figure 6.** Comparison of stacked gathers before and after bin homogenization, (**a**) before bin homogenization, (**b**) after bin homogenization, (**c**) before bin homogenization, (**d**) after bin homogenization.

Figure 7 shows the comparison of stacked gathers before and after bin homogenization. We can determine that the energy of the missing trace is recovered (at the red elliptical mark), and the amplitude of the homogenized bin can reflect the real situation of the underground media.

**Figure 7.** Comparison of stacked gathers before and after bin homogenization.

Due to the influence of engineering construction, noise has a great influence on shallow seismic exploration. Suppressing noise is very important for the accurate identification of goafs. Data analysis shows that the main noise sources include linear interference, surface waves, and abnormally strong amplitudes. In this study, a multidomain joint denoising method is employed through experiments to classify and eliminate different types of noise. This follows the law of progressive, step-by-step denoising and gradually improves the signal-to-noise ratio. For linear interference, the apparent speed and frequency range of linear interference is analyzed, and the noise attenuation method with less waveform modification is adopted to eliminate the noise. For the outliers and abnormal amplitudes, the regional anomalous amplitude attenuation method and the frequency-division abnormal amplitude attenuation method are used for suppression. The strategy of multiple-trace statistics, single-trace denoising, and frequency-division suppression is adopted to suppress random abnormal amplitude noise and surface waves. Strong energy interference in seismic records is identified automatically in different frequency bands, and the spatial position of noise is determined. According to the defined threshold value and attenuation coefficient, time-varying and space-varying methods are adopted to suppress random abnormal amplitude noise. Figure 8 shows the comparison of single-shot records before and after suppressing noise. Within the scope of the red rectangular box, we can determine that the surface waves and other noises are effectively suppressed, and the energy of the primary reflection wave is obviously enhanced. Figure 9 is the comparison of stacked gathers before and after suppressing noise. Within the scope of the red elliptical frame, it can be seen that the weak waveform after noise suppression is displayed, which can more effectively reflect the subtle changes in underground media.

This study identifies goaf and water accumulation characteristics by three steps: seismic forward modeling assistant interpretation, goaf seismic-wave-anomaly characteristics interpretation, and seismic- and electrical-method comprehensive-anomaly interpretation. First, a goaf geological model suitable for the exploration area is established, and wave forward modeling is carried out to study the seismic response characteristics of various goafs (Figure 10) and assist in anomaly interpretation.

**Figure 8.** Comparison of single-shot records before (**left**) and after (**right**) noise suppression.

**Figure 9.** Comparison of stacked gathers before (**upper**) and after (**lower**) suppressing noise.

**Figure 10.** Forward modeling and seismic profile of goaf of different types and in different positions.

Due to the air or water content in the goaf, the absorption coefficient of the goaf is significantly higher than that of the normal coal seam, resulting in no reflection or weak reflection at the bottom interface of the goaf and substratum [32]. There is a fault in the wave group; the reflected wave is characterized by energy enhancement, frequency reduction, and an obvious concave phase (Figure 11a), or the energy becomes weaker and the continuity becomes worse (Figure 11b).

**Figure 11.** Seismic response of goaf in actual seismic profiles. (**a**) Increased energy and low frequency, and (**b**) weak energy and worse continuity.

Comparing the reflected waves, seismic attributes, and drilling verification results of coal seams in the study area, a total of four coal-seam goafs are explained. The goaf labeled <sup>1</sup> is located at the northern boundary of the exploration area, and the goaf anomaly is irregular. Combined with the interpretation results of the northern transient electromagnetic method, the goaf is produced by mining from north to south in the Yuling coal mine, and the abnormal area in the seismic work area is approximately 1359 m2.

The goaf marked with <sup>2</sup> is located outside the southeast exploration boundary, in front of the Xiaoxiyao mine roadway, and the boundary of the abnormal area is irregular. The goaf is caused by mining from east to west in the Xiaoxiyao mine, and the abnormal area in the area is approximately 7352 m2.

The goaf marked with <sup>3</sup> and <sup>4</sup> is located southeast of the exploration boundary, and the abnormal boundary is irregular. The goaf is formed by the westward mining of the Xiaoxiyao mine. However, the northern part of the abnormal area is greatly disturbed by the raw coal transportation channel, and the southern part may be affected by the deterioration of data quality due to its close proximity to the boundary of the detection area, so the reliability of the results needs to be verified by more drilling results.

There is a seismic anomaly zone of coal-seam structure change in coal seam 9, which is located in the west of the exploration area. The plane shape is an irregular polygon, and the east side is connected with the seismic anomaly of the collapse column. The scope is delineated by a dotted line in Figure 12. The drilling data confirmed that the reason was that coal seams 9–1 and 9–2 merged into one layer.

A collapse column was identified via the drilling and log data from coal seam 4 (Figure 13). The lithology below K2 sandstone, as the marker bed of the No. 4 coal roof, is breccia, accompanied by a small amount of coal chips and blocks, with poor cementation. The shape of the collapse column on the plane is an approximate ellipse, especially in the instantaneous frequency attribute, and the characteristics of the collapse column are the most obvious. The collapsed column is a good water channel and water-accumulating area, so it is very important to identify the development of the coal-seam collapse column.

Based on the seismic interpretation of the development characteristics of special geological bodies such as goafs, faults, and collapse columns, the transient electromagnetic method is employed to predict the water accumulation in the key areas of the study area. The transient electromagnetic survey lines covering this area are N0 and N35. Due to

terrain limitations, we adopted small coils and multiple circles of central loop to detect water accumulation in the goaf using the sensitivity of electrical resistivity to water.

**Figure 12.** Distribution plan of seismic-interpretation abnormal areas of coal seams 4 (**left**) and 9 (**right**), -<sup>1</sup> –-4 are the number of the goafs.

**Figure 13.** Comparison of logging and drilling properties of collapsed columns.

The length of survey line N35 is 910 m, and it is located on the 1320 platform. There are two small-scale low-resistance anomalies with relatively small ranges and amplitudes, which are located near point 310 and point 570 (at the red circle of the Figure 14). Among them, the low-resistivity anomaly of measuring point 570 is also shown near the corresponding position of line N0, and the anomaly is caused by the enhanced water yield of the formation near the fault. The abnormal position of point 310 is consistent with the

position of the collapse column interpreted by the seismic attribute, which is speculated to be caused by the increased water abundances of the collapse column.

**Figure 14.** One-dimensional inversion profile of the apparent resistivity of the N35 line.

The length of survey line N0 is 790 m, which is located on the 1350 platform. The apparent resistivity profile east of survey point 530 (in the direction of the large point) shows a large low-resistivity anomaly (at the red circle of Figure 15). Combined with the comprehensive interpretation of adjacent survey line N35 and three-dimensional seismic data, it is speculated that the anomaly here is caused by the fragmentation of the formation near the fault and the enhancement in water abundances.

**Figure 15.** One-dimensional inversion profile of the apparent resistivity of the N35 line.

Integrating the results of the seismic interpretation of the goaf and the transient electromagnetic interpretation of the low-resistance abnormal area of each survey line, the scope of the goaf and water accumulation, collapse column, and fault water-rich anomaly displayed in the detection area is delineated on the plane. The subsequent drilling and testing results are basically consistent with the prediction results; in particular, the accurate prediction of the scope of the goaf plays a vital role in safe mining and personnel safety. Figure 16 shows the verification hole DLT-18-6 passing through goaf No. 4. The verification hole DLT-18-6 in the goaf with number 4, compared with the information of the adjacent well, revealed that the section where no coal core was found happened to be the depth section of coal 4 and 9. Due to the presence of goaf, the diameter value of the coal-seam section shows significant fluctuations, and the gamma–gamma value is significantly lower than that of the No.3 coal seam. Combined with the prediction of water accumulation, relevant coal-mining measures can be effectively formulated to avoid production risks (Figure 17).

The local part of the study area is particularly affected by the on-site construction conditions, such as the mining party's heavy vehicles, gas stations, drilling areas, and blasting areas, which are close to the geophysical exploration work area or even overlap and have a large interference with the seismic signal. In addition, the local coverage times are low, and the signal-to-noise ratio of the data is poor, which directly leads to the poor reliability of the interpretation results. For example, the southwest corner of the study area is a typical disturbed area, with a large area of abnormal areas. However, we predicted

the possible results in advance, and the results were confirmed by a very small amount of drilling. Therefore, the judgment of the reliability of data is a very important link.

**Figure 16.** Borehole Histogram of the DLT-18-6 (The red line is GGFR, the blue line is CAL, and the green line is CAL).

**Figure 17.** Plan distribution map for predicting risk areas of goaf and water accumulation.

To ensure safe and efficient mining, it is recommended to deploy ground drilling or advance drilling and geophysical exploration in the abnormal area and poorly controlled area in advance during the mining process to further determine the nature and location of the abnormality and ensure safety [33]. Combined with relevant geological data such as tunneling and drilling in different stages, comprehensive geological research should be strengthened to further improve the reliability and accuracy of interpretation results.

#### **5. Conclusions**

For the shallow-buried, deep, and large-gap open-pit mine goaf and water accumulation prediction problem, we adopt an integrated geophysical prospecting strategy. To solve the problems of uneven coverage and noise caused by ultrashallow layers and high drops, we propose a direct bin homogenization and multidomain joint denoising method. The results show that the shallow three-dimensional seismic method can clearly reflect the reflected-wave variation features of the seam. The features of a junction between a normal seam and goaf vary obviously, the effect of the set goaf boundary and scope is good, and the resolution is high. With a detailed analysis of the survey-line repeated area of the shallow three-dimensional seismic and transient electromagnetic methods, along with the apparent resistivity threshold scope of the transient electromagnetic method in the normal seam, goaf and goaf water-accumulation area are determined and can be applied to guide the explanation work. The transient electromagnetic method could especially reflect sensitivity to the low resistance, judge the water accumulation in the goaf, and effectively avoid the limit of the shallow three-dimensional seismic method to judge the water accumulation in the goaf. These results provided useful references for predicting goaf and water accumulation of open-pit mines with ultrashallow and high-drop conditions.

**Author Contributions:** Conceptualization, S.Z.; methodology, S.Z. and W.G.; validation, S.C. and Q.C.; investigation, Q.C. and Q.M.; writing—original draft preparation, S.Z. and Q.M.; writing—review and editing, S.Z. and S.C.; visualization, W.G. and Y.D.; supervision, S.C.; project administration, S.C. and Q.M.; funding acquisition, S.Z. and S.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by the National Natural Science Foundation of China (42102212), the unveiling projects of the Department of Science and Technology of Shanxi Province, grant number (20191101018); and the Research and Development Project of Yangquan Coal Industry (Group) Co., Ltd. (GY18027), the National Natural Science Foundation of China (U21A20107).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.
