*2.2. A Simple Example*

We took a circular cave model (Figure 2a) as an example to illustrate the workflow of MVSS beamforming. The seismic data of the 10th shot obtained by finite-different time-domain (FDTD) numerical simulations (Liu et al., 2018 [29]) are shown in Figure 2b. The receivers were arranged from 0 m to 30 m with an interval of 0.2 m. There were 151 receivers in total.

For example, we took point A (16 m, 7 m) as the target point (Figure 2). The delay time was calculated as Equation (A1), and the results are shown in Figure 3. The 151 black dots represent the delay time of each receiver, and the 20 red circles represent the delay time of each shot. The first source shared the same spatial position and delay time with the 36th receiver. We moved the target point and repeated the delay time calculation using Equation (A1).

As per the definition of a subarray in Equation (A11), we divided the array of 151 receivers into 76 subarrays, and each subarray was composed of 75 receivers. Taking the first subarray as an example, we calculated the weight of each receiver in the subarray (Figure 4b) according to the DL and MV method using Equation (A13). We multiplied the 75 delayed signals (Figure 4a) with the weight of each receiver (Figure 4b) and then added them together to obtain the amplitudes on the axis line *x* = *xA* (Figure 4c). Two reflection interfaces were observed from the curve. The delay time was reduced by *t* = *tm* + *tn* = 0.0186 s to exclude the delay time of point A itself. Before the 0.0186 s was excluded, the delay time is not absolutely correct. If we see point A as the target point, the delay time of the receiver which is above point A should be 0. In delay time calculations according to Equation (A1), the reference point is not point A (16,7) but it is the point on x axis (16,0). This results in the 0.0186 s of delay time of point A itself. So, the 0.0186 s should be excluded from the delay time and the correct time-axis position of point A can be obtained. This can also be found in Equation (A3) as (*t* − *τm*). Then, we repeated the process above with each subarray to obtain the amplitude curve on the axis line *x* = *xA* (Figure 4d). In this certain case, the amplitude of a reflected wave was 0.005 times that of a surface wave in seismic signal. Eventually, the amplitude was 0.05 times that of a surface wave on the superposed curve for all subarrays. Then, we picked the peak amplitude and the corresponding estimated time (0.0093 s) as the amplitude result of point A.

**Figure 2.** (**a**) Circular cave model. (**b**) Signal of the 10th shot with 20 stacked traces.

**Figure 3.** Time delay of each receiver and source.

We repeated the delaying, weighting, and superposition processes with the subarrays at every point in the imaging area to obtain the preliminary image for this shot (Figure 5a). Then, the CF (Figure 5b) was obtained according to Equation (A15). The CF matrix was used as another weight to enhance the SNR by multiplying it with the preliminary image. In this step, surface wave artifacts were suppressed because the adopted methods were developed for reflected waves to only amplify the reflected signals with focus. Other interferences were suppressed for the same reason. Then, we obtained the final image of the 10th shot.

**Figure 4.** Delayed signal and weighted superposition of point A. (**a**) Delayed signal of the 1st subarray. (**b**) Weight of each receiver in the 1st subarray according to the MV method. (**c**) Superposed curve of weighted signals. (**d**) Amplitude curve on axis line *x* = *xA*.

**Figure 5.** Preliminary image multiplied with the CF matrix to obtain the time-domain image of one shot. (**a**) Preliminary image. (**b**) CF matrix. (**c**) Image result of the 10th shot.

The time-domain image was converted into a depth-domain image according to the background velocity. The same processing scheme, including beamforming imaging and time-depth conversion, was applied to the remaining shots. The final image was a stacked result of all shots (Figure 6).

### *2.3. Related Works*

### 2.3.1. Kirchhoff Migration

Kirchhoff migration was adopted as a related work in this paper. The methodology is mainly from the authors of [30]. The involved work migrates a single shot record using prestack time migration. The algorithm is a simple travel time path summation with few options. Travel time from source to scatter point (also known as a target point in beamforming methods) is approximated by a Dix equation using the RMS velocity from the model at the lateral position halfway between the source and receiver and at the vertical travel time of the scatter point. Similarly, from the scatter point to a receiver, a Dix equation using the RMS velocity at halfway between the scatter point and the receiver is used. There

is no topographic compensation. The source and receivers are assumed to be on the same horizontal plane.

**Figure 6.** Stacked MVSS beamforming image of 20 shots.
