*3.1. Numerical Experiments*

Six synthetic models were used to validate the proposed imaging method. For this purpose, we constructed an area with a height of 20 m and a width of 30 m (Figure 7) to simulate near-surface seismic exploration. The seismic data came from numerical simulations based on staggered-grid finite differences in the time domain [29,31]. A perfectly matched layer (PML) was used to absorb reflections on the left, right, and bottom boundaries of the model. The upper boundary of the model was a free boundary.

**Figure 7.** Sketch map of the geological model.

According to the near-surface geological conditions, the background P-wave velocity was set as 1500 m/s, the S-wave velocity was set as 700 m/s, and the density was set as 1.8 g/cm3. A 600 Hz Ricker wavelet was used as the source at the ground surface. Twenty shots were deployed with an interval of 1 m, and 151 receivers were evenly arranged at a 0.2 m intervals across the top of the model.

The signal time length was 0.05 s. According to the size of the area and background velocity, the estimated P-wave travel time of a vertically down and up path was 0.026 s and the travel time associated with the far-right source, reflection at bottom left, and farright receiver pair was 0.045 s (Figure 8). Therefore, the target P-wave signal integrity was guaranteed, which meant that a 0.05 s signal would be long enough to cover all the reflected objectives and the reflection of most part of the area could be received by multiple receivers.

**Figure 8.** Illustration of travel time.

All parameters in the numerical experiments are shown in Table 1. Note that the sampling frequency was 2 × 10<sup>5</sup> which is consistent with the time step of the numerical simulation. In the numerical experiments, adapting a shorter recorded signal can help save computational resources. However, in real cases, the recorded signal can be longer in order to obtain more information.

**Table 1.** Numerical experiment parameters.


The receivers are velocity sensors that collect vibrations perpendicular to the ground surface. The geological models are produced by a random near-surface model generation algorithm (Figure 9). For example, horizontal parallel layers (1), folds (2), random fluctuations (3), faults (4) and caves (5) were incorporated into the models. More details on the model generation methods can be found in Appendix B.

**Figure 9.** Geological model generation workflow.

The structural details of each model are shown in Table 2. The models in group A contained cave-type geological anomalies. The models in group B were composed of multiple strata. The models in group C were combinations of cavities and multiple strata. The medium parameters are shown in Table 3.
