*2.2. Stained VDWE*

To extract the wavefield associated with the weakly illuminated structure, we establish a governing equation. According to the principle of staining algorithm [9], the variables *V*, σ and *A* in VDWE are extended to complex domain, and possess the following forms:

$$\mathbf{V} = \overline{\mathbf{V}} + i\overline{\mathbf{V}} = \begin{bmatrix} \overline{\boldsymbol{\upsilon}\_{\mathbf{x}\prime}} \ \overline{\boldsymbol{\upsilon}\_{z}} \end{bmatrix}^{T} + i\begin{bmatrix} \overline{\boldsymbol{\upsilon}\_{\mathbf{x}\prime}} \ \overline{\boldsymbol{\upsilon}\_{z}} \end{bmatrix}^{T},\tag{5}$$

$$\mathfrak{o} = \overline{\mathfrak{o}} + i\overline{\mathfrak{o}} = \begin{bmatrix} \overline{\mathfrak{r}\_{xx}} \ \overline{\mathfrak{r}\_{zz}} \ \overline{\mathfrak{r}\_{xz}} \end{bmatrix}^T + i \begin{bmatrix} \overline{\mathfrak{r}\_{xx}} \ \overline{\mathfrak{r}\_{zz}} \ \overline{\mathfrak{r}\_{xz}} \end{bmatrix}^T,\tag{6}$$

$$A = \overline{A} + i\tilde{A} = \begin{bmatrix} \left(\overline{c\_P} + i\tilde{c\_P}\right)^2 & \left(\overline{c\_P} + i\tilde{c\_P}\right)^2 - 2\left(\overline{c\_S} + i\tilde{c\_S}\right)^2 & 0\\ \left(\overline{c\_P} + i\tilde{c\_P}\right)^2 - 2\left(\overline{c\_S} + i\tilde{c\_S}\right)^2 & \left(\overline{c\_P} + i\tilde{c\_P}\right)^2 & 0\\ 0 & 0 & \left(\overline{c\_S} + i\tilde{c\_S}\right)^2 \end{bmatrix} \tag{7}$$

where *<sup>i</sup>* <sup>=</sup> √−1, variables with <sup>−</sup> and <sup>∼</sup> symbols, respectively, represent the real and imaginary parts. Note that setting *cP* and *cS* is the core step in extracting the wavefield. *cP* and *cS* are consistent with the P- and S-wave velocity. *<sup>c</sup><sup>P</sup>* and *<sup>c</sup><sup>S</sup>* are obtained by a dimensionless multiplied *cP* and *cS* near the exploration target, while *<sup>c</sup><sup>P</sup>* and *<sup>c</sup><sup>S</sup>* are equal to 0 in the unconcerned area. The location where *<sup>c</sup><sup>P</sup>* and *<sup>c</sup><sup>S</sup>* are non-zero is referred to as the stained target.

By substituting Equations (5)–(7) into the VDWE, we obtain,

$$\begin{cases} \displaystyle \rho \frac{\partial (\nabla + i \tilde{\mathcal{V}})}{\partial t} = \mathcal{L} (\overline{\sigma} + i \overline{\sigma})\\ \displaystyle \frac{\partial (\overline{\sigma} + i \tilde{\mathcal{W}})}{\partial t} = \rho \left( \overline{\mathcal{A}} + i \overline{\mathcal{A}} \right) \mathcal{L}^T \left( \overline{\mathcal{V}} + i \overline{\mathcal{V}} \right) \end{cases} \tag{8}$$

$$\begin{cases} \begin{array}{l} \frac{\partial \left(\overline{\tau^{p}} + i\overline{\tau^{p}}\right)}{\partial t} = \rho \left(\overline{c\_{P}} + i\overline{c\_{P}}\right)^{2} \left(\frac{\partial \left(\overline{v\_{\varepsilon}^{\prime}} + i\overline{v\_{\varepsilon}^{\prime}}\right)}{\partial x} + \frac{\partial \left(\overline{v\_{\varepsilon}^{\prime}} + i\overline{v\_{\varepsilon}^{\prime}}\right)}{\partial z}\right) \\ \rho \frac{\frac{\partial \left(\overline{v\_{\varepsilon}^{\prime\prime}} + i\overline{v\_{\varepsilon}^{\prime\prime}}\right)}{\partial t} = \frac{\partial \left(\overline{\tau^{p}} + i\overline{\tau^{p}}\right)}{\partial z} \end{array}}{\frac{\partial \left(\overline{v\_{\varepsilon}^{\prime\prime}} + i\overline{v\_{\varepsilon}^{\prime\prime}}\right)}{\partial z}} \end{cases} \tag{9}$$

Taking the real parts of Equations (8) and (9), we acquire an equation expressed as,

$$\begin{cases} \begin{array}{c} \rho \frac{\partial \overline{V}}{\partial t} = \mathcal{L} \overline{\sigma} \\ \frac{\partial \overline{\sigma}}{\partial t} = \rho \left( \overline{A} \mathcal{L}^T \overline{V} - \overline{A} \mathcal{L}^T \overline{V} \right) \\ \frac{\partial \overline{\tau'}}{\partial t} = \rho u \left( \frac{\partial \overline{\tau'}}{\partial x} + \frac{\partial \overline{\sigma}\_i}{\partial z} \right) - 2\rho \beta \left( \frac{\partial \overline{\upsilon'}\_i}{\partial x} + \frac{\partial \overline{\upsilon}\_i}{\partial z} \right) \\ \rho \frac{\partial \overline{\upsilon''\_i}}{\partial t} = \frac{\partial \overline{\tau''}}{\partial x} \\ \rho \frac{\partial \overline{\upsilon''\_i}}{\partial t} = \frac{\partial \overline{\tau''}}{\partial z} \end{array} \tag{10}$$

where *α* = *cP* <sup>2</sup> <sup>−</sup> *<sup>c</sup><sup>P</sup>* <sup>2</sup> and *<sup>β</sup>* <sup>=</sup> *cPcP*. The real parts of S-waves are defined as *<sup>v</sup><sup>S</sup> <sup>x</sup>* = *vx* − *v<sup>P</sup> <sup>x</sup>* and *vS <sup>z</sup>* = *vz* − *v<sup>P</sup> <sup>z</sup>* . The imaginary parts of Equations (8) and (9) can be written as,

$$\begin{cases} \rho \frac{\partial \bar{V}}{\partial t} = \mathcal{L}\bar{\sigma} \\ \frac{\partial \bar{\sigma}}{\partial t} = \rho \left( \overline{A} \mathcal{L}^T \bar{V} - \overline{A} \mathcal{L}^T \overline{V} \right) \\ \frac{\partial \bar{\tau}^{\overline{p}}}{\partial t} = \rho u \left( \frac{\partial \bar{v}\_t^{\overline{\tau}}}{\partial x} + \frac{\partial \bar{v}\_t^{\overline{\tau}}}{\partial z} \right) + 2\rho \beta \left( \frac{\partial \overline{\tau}\_t}{\partial x} + \frac{\partial \overline{\tau}\_t}{\partial z} \right) \\ \rho \frac{\partial \overline{v}\_t^{\overline{p}}}{\partial t} = \frac{\partial \overline{\tau}^{\overline{p}}}{\partial x} \\ \rho \frac{\partial \overline{v}\_t^{\overline{p}}}{\partial t} = \frac{\partial \overline{\tau}^{\overline{p}}}{\partial z} \end{cases} \tag{11}$$

The imaginary parts of S-waves are defined as *v* -*S <sup>x</sup>* <sup>=</sup> *<sup>v</sup><sup>x</sup>* <sup>−</sup> *<sup>v</sup>* -*P <sup>x</sup>* and *v* -*S <sup>z</sup>* <sup>=</sup> *<sup>v</sup><sup>z</sup>* <sup>−</sup> *<sup>v</sup>* -*P <sup>z</sup>* . Equations (10) and (11) jointly constitute the stained VDWE (SVDWE). According to Equation (11), we find that *<sup>V</sup>* is only excited when *<sup>V</sup>* arrives at the stained target. The imaginary part *<sup>V</sup>* is named as stained wavefield. The term *<sup>A</sup>*L*T<sup>V</sup>* can be considered as the source signature of *<sup>V</sup>*. *<sup>A</sup>* and *<sup>β</sup>* are dimensionless due to the principle of *<sup>c</sup><sup>P</sup>* and *<sup>c</sup><sup>S</sup>* setting. It is evident that Equation (10) is approximately equivalent to the VDWE.

### *2.3. Inner Product Imaging Condition with Filters*

Imaging condition can convert seismic data into structure images through the focus of the source and receiver wavefield. For the ERTM, the inner product imaging condition has an advantage in terms of definite physical meanings. However, applying such imaging condition inevitably results in migration artifacts because of backscatter. To suppress the unwanted seismic events in the migration profile, we add filters to the imaging conditions. The final images can be expressed as:

$$
\widehat{I^{PP'}} = \sum\_{t=1}^{T} \left( F\_{So}^{P} \overline{\widehat{V\_{So}^{P}}} \right) \cdot \left( F\_{Rx}^{P} \overline{\widehat{V\_{Rx}^{P}}} \right) \tag{12}
$$

$$
\widehat{I^{\widetilde{PS}}} = \sum\_{t=1}^{T} \left( F\_{\text{So}}^{\text{P}} \widehat{\mathbf{V}\_{\text{So}}^{\widetilde{\mathbf{P}}}} \right) \cdot \left( F\_{\text{Rx}}^{\text{S}} \widehat{\mathbf{V}\_{\text{Rx}}^{\widetilde{\mathbf{S}}}} \right) \tag{13}
$$

where .*IPP* and *I* -*PS* are PP and PS images, the symbol · denotes the inner product operation. *<sup>V</sup>* with the subscripts *So* and *Re* denotes the stained wavefields of seismic sources and receivers. The superscripts *P* and *S* represent P- and S-wave, respectively. *T* is the maximum time of seismic records. *FP So* is a filter with the mathematical form shown below:

$$F\_{So}^{P} = \begin{cases} \begin{array}{c} 1, \ \widehat{D\_{So}^{P}} \cdot \mathfrak{n}\_{So}^{P} > 0 \\ 0, \ \mathcal{D}\_{So}^{P} \cdot \mathfrak{n}\_{So}^{P} < 0 \end{array} \tag{14}$$

where *D* .*<sup>P</sup> So* = *d* -*P <sup>x</sup>* ,*d* -*P z* is the Poynting vector, *n<sup>P</sup> So* is the direction vector. If the result of the inner product between *D* .*<sup>P</sup> So* and *<sup>n</sup><sup>P</sup> So* is positive (negative), then the wavefield is retained (muted). Other filters have a similar expression. The Poynting vector can be calculated by:

$$\begin{cases} \begin{array}{c} \widehat{d}^{\mathfrak{p}}\_{a} = -\widehat{\tau}^{\widetilde{P}} \widehat{\upsilon}^{\widetilde{p}}\_{a} \\ \widehat{d}^{\widetilde{S}}\_{a} = -\left[ \widehat{\tau}^{\widetilde{P}}\_{ab} - \widehat{\tau}^{\widetilde{P}} \delta(a-b) \right] \widehat{\upsilon}^{\widetilde{S}}\_{b} \end{array} \end{cases} \tag{15}$$

where *a* and *b* are the *x*- or *z*-component of the Poynting vectors, and *δ* is the Dirac function. From formula (15), the Poynting vector is 0 if either particle-velocity and stress is 0. There is a position between the crest and the trough where the wavefield is 0, but the Poynting vector at that position should not be 0. To improve the stability of equation (15), the local least squares strategy [25] is applied.

### **3. Numerical Examples**

In the first example, we employ a layered velocity model to probe the performance of the above equations. The P-wave migration velocity model is sketched in Figure 1a. The S-wave migration velocity is obtained using an empirical formula *cS* = *cP*/ <sup>√</sup>3. The density is set as 2000 kg/m3. The model is dispersed by 650 × 400 grid numbers, with a grid size of 10 × 10 m2. A Ricker wavelet with the dominant frequency of 25 Hz is set as the source signature at the location of (1 km, 0). The wavefields are simulated by the finite difference method with the second-order in time and tenth-order in space. In addition, we add a perfectly matched layer around the simulated area to absorb spurious reflections from the truncated boundary.

Based on the VCWE, the *x*-component of the particle-velocity is shown in Figure 2a, in which the P- and S-waves are coupled together. When the coupled wavefield is extrapolated by solving Equation (4), it can be completely decoupled into pure P- and S-waves (Figure 2b,c). Notably, the seismic energy of the deep structure information is relatively weak in the decoupled wavefield.

**Figure 1.** Migration model: (**a**) P-wave velocity; (**b**) Stained target.

**Figure 2.** *x*-component of particle-velocity wavefield: (**a**) Coupled wave; (**b**) Pure P-wave; (**c**) Pure S-wave.

To reduce the shielding effect of the overburden rock on seismic waves, we simulate the wavefield based on SVDWE. Here, the second interface is set as the stained target. *<sup>c</sup><sup>P</sup>* is obtained by *cP* multiplied 10−<sup>6</sup> in the target area while is zero in other locations. The setting of *<sup>c</sup><sup>S</sup>* is similar to that of *<sup>c</sup>P*. Considering the similarity between *<sup>c</sup><sup>P</sup>* and *<sup>c</sup><sup>S</sup>* models, only one of them is shown in Figure 1b. By solving SVDWE, the stained wavefield *<sup>V</sup>* is generated after the *<sup>V</sup>* hits the stained target. The stained wavefield *<sup>V</sup>* only contains the reflection and transmission about the second interface as exhibited in Figure 3. Consequently, we infer that imaging with the stained wave can reduce interference from non-target structures.

**Figure 3.** *x*-component of the stained wavefield: (**a**) Pure P-wave; (**b**) Pure S-wave.

Furthermore, we utilize a contrastive approach to illustrate the adaptability of the ERTM for different equations. We evenly deploy 10 shots and 650 receivers on the surface of the layered velocity model. The total length of seismic records is 3 s, with a step time of 1ms. Figure 4 shows the multi-shot stacked images based on VDWE. Moreover, we implement the ERTM using the SVDWE, the PP and PS images are sketched in Figure 5. The settings of the stained target are the same as those

of previous numerical experiments. Thanks to the stained wavefields, the primary information in the images is the second interface. However, some unwanted migration artifacts may interfere with subsequent interpretations.

**Figure 4.** ERTM image obtained based on the VDWE: (**a**) PP imaging; (**b**) PS imaging.

**Figure 5.** ERTM image based on the SVDWE: (**a**) PP imaging; (**b**) PS imaging.

Next, we demonstrate that the proposed imaging condition plays an important role in eliminating migration artifacts. The filters in the imaging condition need to be prepared in advance, and their properties are determined by the direction vector and Poynting vector. For the snapshot shown in Figure 3a, the Poynting vector (Figure 6) can be obtained using Equation (15). Here, we exhibit four filtered wavefields (Figure 7) that, respectively, correspond to different direction vectors: *n<sup>P</sup> So* = (1, 0), *nP So* <sup>=</sup> (−1, 0), *<sup>n</sup><sup>P</sup> So* <sup>=</sup> (0, <sup>−</sup>1), and *<sup>n</sup><sup>P</sup> So* = (0, 1).

**Figure 6.** Poynting vector: (**a**) *x*-component; (**b**) *z*-component.

**Figure 7.** Filtered wavefield: (**a**) downgoing; (**b**) upgoing; (**c**) Leftgoing; (**d**) rightgoing.

For flat strata, images with fewer artifacts can be obtained in the inner product of the upgoing stained source wavefield and the down-going stained receiver wavefield (Figure 8). Figure 9 quantitatively compares PP and PS images, and the two sets of single-channel data are extracted from Figure 8 (from 2 km to 3 km in depth), and then mapped to the wavenumber domain through the Fourier transform. It is prominent that PS images are characterized by abundant high wavenumber information, robustly confirming that PS images have a higher spatial resolution than PP images. The basic reason for such a phenomenon is that S-wave travels slower than P-wave in the subsurface. Namely, the converted S-waves possess a shorter wavelength.

**Figure 8.** ERTM profile with fewer artifacts compared with Figure 5: (**a**) PP imaging; (**b**) PS imaging.

**Figure 9.** Comparison between PP and PS images: (**a**) Single-channel data; (**b**) Wavenumber spectrum.

In this example, a salt model is applied to evaluate the efficiency of the proposed approach in complex geologic regions. The true P-wave velocity model is exhibited in Figure 10a, and S-wave velocity is also calculated using the previous empirical formula. The model size is 400 × 400 grids, with a grid size of 10 × 10 m<sup>2</sup> in both horizontal and vertical directions. A total of 40 sources are evenly excited on the surface of the model. Each source has 400 receivers with a spacing of 10 m. The maximum recording time is 6 s. Figure 11 shows the PP and PS images obtained by the conventional ERTM, and it is obvious that the shallow structure and the salt boundary can be effectively identified. However, the high-velocity overburden rock leads to insufficient illumination for the deep subsurface, which ultimately blurs the subsalt image. We set the stained target in Figure 10b. Figure 12 displays

the PP and PS images with the proposed imaging condition. We can see that the two complementary images are favorable for the identification of subsalt faults.

**Figure 10.** Migration model: (**a**) True P-wave velocity; (**b**) Stained target.

**Figure 11.** ERTM images based on the VDWE: (**a**) PP imaging; (**b**) PS imaging.

**Figure 12.** ERTM images based on the SVDWE: (**a**) PP imaging; (**b**) PS imaging.

In the last experiment, a field data is utilized to confirm the effectiveness of our method. Figure 13 displays the migration velocity model that has 940 × 300 grid points with a 20 × 20 m<sup>2</sup> element. There are 240 sources unevenly distributed at the free surface. Figure 14 show some shot gathers whose maximum recording time is 5 s. Since we only have the vertical component of the

seismic record where P-wave information dominates, the valuable PS images are not obtained. The PP images based on VDWE are sketched in Figure 15, from which we can see the main structures. To enhance the quality of imaging in weakly illuminated areas, we set the box in Figure 13 as the stained targets. By using the proposed method, the migration events (Figure 16) relevant to the stained targets are highlighted, and the reflections of nontarget structures are muted.

**Figure 14.** Two shot gathers.

**Figure 15.** PP imaging based on VDWE.

**Figure 16.** PP imaging based on SVDWE.

### **4. Discussion**

Our approach shows the satisfying performance for imaging weakly illuminated structures. When an overall migrated profile is acquired, it is proper to achieve locally enhanced imaging by artificially setting exploration targets in advance. Since the SVDWE contains two systems of equations, the algorithm of such an approach is twice as computationally intensive as the conventional one. The staining strategy is intrinsically related to the imaging results and is worthy of further investigation. The filters in the inner product imaging condition screen different wavefields according to the direction of wave propagation. Theoretically, it is possible to achieve wavefield decomposition in an arbitrary direction. However, in this paper, we decompose the wavefield into up- and downgoing categories under the assumption that the target structure is relatively flat. To enhance the imaging adaptability for complex geological situations, the wavefield screening method according to stratigraphic dip angle needs to be persistently explored.
