In a homogeneous and isotropic medium, the elastic wave extrapolation (Aki and Richards, 1980) can be expressed as follows: ρ u ¨ = ( λ + 2 μ ) ∇ ( ∇ ⋅ u ) − μ ∇ × ∇ × u ,   (1) where u and u ¨ are the displacement vector wavefield and its second-order time derivative; λ , μ and ρ are the Lame’s moduli and density, respectively. Based on the Helmholtz theorem (Dellinger and Etgen, 1990), the elastic wavefield in an isotropic case can be separated into a curl-free P wavefield ( ∇ × u P = 0 ) and a divergence-free S wavefield ( ∇ ⋅ u S = 0 ) . u P and u S are the P-wave and S-wave displacement vector wavefields. Analogous to the separation of displacement wavefield, the second-order time derivative of displacement wavefield can be decomposed as u ¨ = u ¨ P + u ¨ S , where { u ¨ P = ( λ + 2 μ ) ∇ ( ∇ ⋅ u ) u ¨ S = − μ ∇ × ∇ × u u = u P + u S . (2) Decoupled Equation 2 is embedded in the update of the displacement wavefield. The P and S wavefields are constructed by the first two equations, respectively, and their summation can obtain the total elastic wavefield in the third equation (Ma and Zhu, 2003; Li, 2007; Zhu, 2017). In contrast to the summing of decoupled P wavefield and S wavefield, the decoupled S wavefield can be constructed by subtracting the P wavefield from the total elastic wavefield. Decoupled Equation 2 produces displacement vector wavefields of pure P- and S-waves. Correspondingly, the first-order stress-particle velocity wave equation has been proposed (Li, 2007; Du et al., 2017; Zhou et al., 2018): { τ ˙ = ( λ + 2 μ ) ∇ ⋅ v − μ ( ∇ v + ∇ v T ) ρ v ˙ = ∇ ⋅ τ τ ˙ p = ( λ + 2 μ ) ∇ ⋅ v ρ v ˙ P = ∇ τ p v S = v − v P , (3) where v and τ are particle velocity and stress of elastic wave, ∇ and ∇ ⋅ represent the operators of gradient and divergence, superscripted T represents the transpose, and subscripted P and S represent the P wave and S wave, respectively. Firstly, the particle velocity and stress tensor of elastic wave and the synthetic seismic records are computed by the first two equations, the conventional stress-particle velocity wave equation. Then, the auxiliary wavefield τ p can be constructed by the third equation τ ˙ p = ( λ + 2 μ ) ∇ ⋅ v and is used to compute the P-wave particle velocity v P . Finally, the S-wave particle velocity can be constructed by subtracting P wavefield particle velocity v P from total elastic wavefield particle velocity v . The source and receiver wavefield can be generated by the forward and backward extension, respectively, based on the decoupled wave equation. The decoupled wavefields are all vector, and their amplitude and phase are consistent with the original elastic wavefield.

For the decoupled vector wavefields, we obtain scalar imaging results by imaging conditions, including cross-correlation imaging condition of Cartesian components generating too many results to interpret, and dot product–based scalar imaging conditions. Regardless of source normalization, the dot product–based scalar imaging condition (Du et al., 2017; Zhu, 2017; Yang et al., 2018) for the PP wave can be written as follows: I p p ( x ) = ∫ s p ( x , t ) ⋅ r ¯ p ( x , t ) d t (4) in terms of source particle velocity vector s p and receiver particle velocity vector r p . Here, I p p is the migrated PP image by integrating the dot product over time t , symbol “ ⋅ ” denotes the dot product of two vectors, and tilde above wavefield variable denotes its conjugation.

Algebraically, the dot product is the sum of some related Cartesian components products, which means s p ⋅ r p = s x p r ¯ x p + s y p r ¯ y p + s z p r ¯ z p . Since the Cartesian components are independent over time t , the migrated PP image I p p can be disintegrated into three parts I p p x x , I p p y y , and I p p z z , where I p p x x = ∫ s x p r ¯ x p d t , I p p y y = ∫ s y p r ¯ y p d t , and I p p z z = ∫ s z p r ¯ z p d t represent cross-correlation imaging results of the x-axis, y-axis, and z-axis Cartesian components, respectively.

By introducing the opening angle θ shown in Figure 1, the dot product scalar imaging condition can be equivalently expressed as follows: I p p = ∫ | s p | | r ¯ p | cos ⁡ θ d t . (5) Here, | ⋅ | is the modulus of a vector. The amplitude of the PP image depends on the modulus of the incident wave, modulus of the reflected wave, and the extra weighting factor cos ⁡ θ . Depending on the seismic source and Green function between the source and scattering point, the modulus of the incident wave is desired in the PP image. The modulus of a reflected wave depends on the Green function between the receiver and scattering point and the reflection coefficient R p p . The reflection coefficient R p p quantitatively describes the amplitude and phase of the reflected wave while P wave is incidence on the interface. The modulus of the incident and reflected waves is desired information in an image to provide a reliable basis for seismic interpretation inversion. Regardless of wavefields modulus, the additional factor cos ⁡ θ will cause destructive interference in the final PP image. On the one hand, this extra factor cos ⁡ θ changes its sign when the opening angle θ > 90 ° or the incident angle α > 90 ° , which will cause the polarity reversal problem (Du et al., 2017; Zhou et al., 2018) at a large incident angle. On the other hand, I p p is also contaminated by backscattering artifacts with the incident angles close to 9 0 ∘ or opening angles near 18 0 ∘ .

In acoustic RTM or scalar-based elastic RTM, the Laplace filter (Youn and Zhou, 2001) has been used to suppress the backscattering artifacts in the PP image. As for the PP image in vector-based elastic RTM, the scalar imaging condition with a Laplace filter can be expressed as follows: I p p l a p = ∇ 2 I P P        = ∇ 2 ∫ ( s p ( x , t ) ⋅ r ¯ p ( x , t ) ) d t . (6) Here, I p p l a p is the migrated PP image with a Laplace filter, ∇ 2 = ∂ x 2 + ∂ y 2 + ∂ z 2 is the Laplace filter operator, ∂ x 2 , ∂ y 2 , and ∂ z 2 are the second-order partial derivatives along the x-axis, y-axis, and z-axis direction, respectively. Since partial derivatives are independent of time integration, the PP image I p p l a p can also be disintegrated by Cartesian components’ cross-correlation imaging results. As for 2D vector-based wavefields, the PP image I p p l a p can be separated into four items: I p p l a p = ∂ x 2 ∫ s x p r ¯ x p d t + ∂ x 2 ∫ s z p r ¯ z p d t + ∂ z 2 ∫ s x p r ¯ x p d t + ∂ z 2 ∫ s z p r ¯ z p d t         = ∂ x 2 I p p x x + ∂ x 2 I p p z z + ∂ z 2 I p p x x + ∂ z 2 I p p z z . (7) Here, ∂ x 2 I p p x x and ∂ z 2 I p p x x are second-order derivatives of x-axis Cartesian component cross-correlation imaging result along the x-axis direction and z-axis direction, respectively; ∂ x 2 I p p z z and ∂ z 2 I p p z z are second-order partial derivatives of z-axis Cartesian component cross-correlation imaging result along the x-axis direction and z-axis direction, respectively. Thereinto, ∂ x 2 I p p x x and ∂ z 2 I p p z z are parallel-oriented partial derivatives of Cartesian component cross-correlation imaging result. Meanwhile, ∂ z 2 I p p x x and ∂ x 2 I p p z z are normal-oriented partial derivatives of Cartesian component cross-correlation imaging result. The fault model (shown in Figure 2) has been introduced for reverse time migration to highlight the interaction characteristics of these decoupled migrated results on the flat and inclined interface. Figures 3A–D are the decoupled migrated results of ∂ x 2 I p p x x , ∂ x 2 I p p z z , ∂ z 2 I p p x x , and ∂ z 2 I p p z z , respectively.

Backscattering noise has been suppressed in all decoupled images. Moreover, these decoupled images have different migration capabilities. On the one hand, the decoupled items ∂ x 2 I p p x x (shown in Figure 3A) and ∂ x 2 I p p z z (shown in Figure 3B) related to second-order partial derivative associated with x-axis direction show similar migrated images sensitive to inclined structures. Compared with the ideal parallel-oriented result ∂ x 2 I p p x x , the decoupled normal-oriented item ∂ x 2 I p p z z would encounter severe crosstalks, causing destructive interference in the final stacked image. On the other hand, the decoupled items ∂ z 2 I p p x x (Figure 3C) and ∂ z 2 I p p z z (Figure 3D) related to second -order partial derivatives along the z-axis direction migrate good images in flat-layer structures. Compared with the ideal parallel-oriented migrated result ∂ z 2 I p p z z , the decoupled normal-oriented item ∂ z 2 I p p x x fails to image near-zero offset, contrary to the final PP wave image. Furthermore, the phase of ∂ z 2 I p p x x is opposite to ∂ z 2 I p p z z . The opposite phase between ∂ z 2 I p p z z and ∂ z 2 I p p x x will result in the Laplace filter’s disability in correcting the polarity reversal problem in the PP image. The scalar imaging condition with the Laplace filter successfully suppresses backscattering but fails to correct the polarity reversal.

Since normal-oriented partial derivative related items of the Laplace filter are greatly affected by crosstalk noise, the horizontal derivative related decoupled items have been selected to migrate the PP image. Analogous to the scalar imaging condition with the Laplace filter, we propose the scalar imaging condition with the pseudo-Laplace filter. To characterize the scalar imaging condition with a pseudo-Laplace filter mathematically, we first define the pseudo-Laplace operator ∇ ˜ 2 as follows: ∇ ˜ 2 = ( ∂ x 2 , ∂ y 2 , ∂ z 2 ) (8) and the PP wave’s Hadamard product image I P P as follows: I P P = ∫ s p ∘ r ¯ p d t = ( I p p x x , I p p y y , I p p z z ) , (9) where ∘ is the Hadamard operator (see also Appendix A). The pseudo-Laplace operator ∇ ˜ 2 comprises three array components, which are second-order partial derivatives along the x-, y- and z-axis, respectively. Their summation is just the Laplace operator. Meanwhile, the PP wave Hadamard product images I P P are composed of three array components: cross-correlation imaging results of the x-axis, y-axis, and z-axis Cartesian components. The summation of array components is just a PP wave scalar product image I p p . There is some specific connection between the pseudo-Laplace operator ∇ ˜ 2 and the Laplace operator ∇ ˜ 2 and between the PP wave Hadamard product image I P P and the PP wave scalar product image I p p . Unlike the scalar Laplace operator ∇ 2 and PP wave scalar product image I p p , the pseudo-Laplace operator ∇ ˜ 2 and PP wave Hadamard product image I P P are both vectors.

Combining the pseudo-Laplace operator with the PP wave Hadamard product image vector, we propose a scalar imaging condition with a pseudo-Laplace filter for vector-based wavefields as follows: I p p p s e − l a p = ∇ ˜ 2 ⋅ I P P             = ∇ ˜ 2 ⋅ ( ∫ s p ∘ r ¯ p d t ) . (10) Here, I p p p s e − l a p is the PP image by the scalar imaging condition with a pseudo-Laplace filter. As for 2D vector-based wavefield, the scalar imaging condition with a pseudo-Laplace filter for the PP wave should be simplified, and components related to y-axis should be neglected: I p p p s e − l a p = ∂ x 2 ∫ s x p r ¯ x p d t + ∂ z 2 ∫ s z p r ¯ z p d t             = ∂ x 2 I p p x x + ∂ z 2 I p p z z . (11) The scalar imaging conditions with a Laplace filter (Eq. 4) and a pseudo-Laplace filter (Eq. 8) depend on second-order partial derivatives of Cartesian components cross-correlation results. Different from the Laplace filter composed of parallel-oriented and normal-oriented items, only the parallel-oriented items are selected to form the pseudo-Laplace filter. As shown in Figure 4, we migrate the PP images of the fault model by the scalar imaging condition with the Laplace filter and with the pseudo-Laplace filter. Figure 4A is the migrated PP scalar image with a Laplace filter, the sum of Figures 3A–D. Similarly, the sum of Figures 3A,D is just one scalar migrated PP image with a pseudo-Laplace filter, as shown in Figure 4B. Compared with the PP image with a Laplace filter (Figure 4A), backscattering noise suppression (marked by red arrows) and polarity reversal correction (marked by red circles) have been shown in the PP image with a pseudo-Laplace filter (Figure 4B).

For the application in elastic RTM, we should migrate three Cartesian component’s images (or two images in a 2D case) by time integration and then sum second-order derivatives of three images together. Compared with the scalar imaging condition, additional storage of three Cartesian component’s images and computation of second-order derivative operation are introduced in scalar imaging conditions with a pseudo-Laplace filter. Compared with the storage and computation costs consumed by the wavefield extrapolation of elastic RTM, the additional calculation introduced by the proposed filter can be ignored to a certain degree. Thus, the scalar imaging condition with a pseudo-Laplace filter is easy to perform and does not require additional processing or storage, which is essential for elastic RTM.

Based on the above-tested fault model, the suppression of backscattering and correction of polarity reversal have been shown in the PP scalar image with the proposed pseudo-Laplace filter. The section will theoretically illustrate how the pseudo-Laplace filter suppresses backscattering noise and correct polarity reversal in PP images with the assumption of a plane wave.

For the vector-based elastic wavefields, the source-side particle velocity vector s p and receiver-side particle velocity vector r p are related to the polarization and propagation of the P wave. Decomposing the particle velocity wavefield of pure P wave in plane waves, we obtain the following: s p = | s p | p s e i k ( n s ⋅ x − v p t ) , (12A) and r ¯ p = | r p | R p p p r e i k ( v p t − n r ⋅ x ) . (12B) Here, | s p | and | r p | = | s p | R p p are the modulus of the incident and reflected wave, respectively. p and n are polarization unit vector and propagation unit vector, respectively. Superscripts s and r represent incident wave and reflected wave. v p , k = ω / v p , and ω are P wave’s propagated velocity, wavenumber and angular frequency, respectively. Assuming that vectors n and p vary slowly in the space-time domain, their temporal and spatial derivatives are small enough to be ignored.

Substituting the plane wave definitions (Eq. 12) into Eqs 4, 7, 11, we obtain an expression with the assumption of the plane wave as follows: I p p = ∫ | s p | | r p | R p p ( p s ⋅ p r ) e i k ( n s − n r ) ⋅ x d t , (13A) I p p l a p = − ∫ | s p | | r p | R p p ω 2 { [ ( n x s − n x r ) 2 + ( n z s − n z r ) 2 ] ( p x s p x r + p z s p z r ) } e i k ( n s − n r ) ⋅ x d t , (13B)

and I p p p s e − l a p = − ∫ | s p | | r p | R p p ω 2 [ ( n x s − n x r ) 2 p x s p x r + ( n z s − n z r ) 2 p z s p z r ] e i k ( n s − n r ) ⋅ x d t . (13C)

In the 2D case of incident pure P wave, the reflected P wave without wave conversion would be obtained in an observation coordinate system along the horizontal surface and vertical depth. As shown in Figure 5A, the geological structure information of the reflector has been introduced in the descriptions of particle velocity vectors in the observation coordinate system. A local coordinate system (as shown in Figure 5B) is constructed along with tangential and vertical directions of the reflector to simplify the representation of vectors in the source and receiver wavefield. As for pure P wave, the polarization vector p s is parallel to the propagation vector n s with the same positive direction. In the local coordinate system, the polarization vector p s and propagation vector n s of incident vector s p in source wavefield should be described by the incident angle α as follows: p s = sin ⁡ α i → + cos ⁡ α k → , (14A) and n s = sin ⁡ α i → + cos ⁡ α k → . (14B) Unlike source wavefield, the pure P wave in receiver wavefield should be described by conjugation of a reflected vector r ¯ p . Its propagation direction is the opposite to that of the reflected wave, and polarization direction is the same as the reflected wave. The polarization vector p s is parallel to the propagation vector n s with the same positive direction for reflected pure P wave r p . When it comes to reflected pure P wave, the polarization vector p r and propagation vector n r should be described in the local coordinate system by the reflected angle α as follows: p r = − sin ⁡ α i → + cos ⁡ α k → , (15A) and n r = sin ⁡ α i → − cos ⁡ α k → . (15B) According to the descriptions of incident vector (Eq. 14) and conjugation of reflected vector (Eq. 15), we can rewrite the scalar imaging condition, the scalar imaging condition with a Laplace filter, and the scalar imaging condition with a pseudo-Laplace filter in the local Cartesian coordinate system as follows: I p p = ∫ | s p | | r p | R p p ⁡ cos ⁡ θ e i k ( n s − n r ) ⋅ x d t , (16A) I p p l a p = − ∫ | s p | | r p | R p p ω 2 2 ⁡ cos ⁡ θ ( cos ⁡ θ + 1 ) e i k ( n s − n r ) ⋅ x d t , (16B)

and I p p p s e − l a p = − ∫ | s p | | r p | R p p ω 2 ( cos ⁡ θ + 1 ) 2 e i k ( n s − n r ) ⋅ x d t . (16C)

Here, θ = 2 α is the opening angle. The above-mentioned imaging algorithms need to be separated into terms related to amplitude and phase. As for the phase-related item, they agree with each other by the form of e i k ( n s − n r ) ⋅ x . The phase-related item is dependent on the illumination vector i s r = n s − n r , defined by Lecomte (2008). The illumination vector i s r satisfies the following relationship i s r = n s − n r = 2 ⁡ cos ( θ / 2 ) n ^ , where n ^ is a unit normal vector at each reflector. In a local coordinate system, a unit normal vector can be described as n ^ = ( 0,1 ) regardless of the inclination of a reflector. Phase-related items are dependent on the opening angle θ . In particular, i s r = 0 while θ = 180 ° . It indicates that i s r will be zero at the reflectors when the incident wave and reflected wave are with the same propagating path.

Otherwise, the amplitude-related items in I p p , I p p l a p , and I p p p s e − l a p imaging algorithms are different from each other. Factor | s p | | r p | R p p , coexisting in The above-mentioned imaging conditions, is useful for seismic inversion interpretation. Once the Laplace filter and pseudo-Laplace filter are introduced, the scalar imaging algorithms are influenced by angular frequency ω 2 . The introduction of ω 2 weakens the amplitude of low-frequency data and enhances the amplitude of high-frequency data, changing the spectrum of images and damaging effective low-frequency information. To maintain the spectrum of images and recover their effective low-frequency information, reasonable time integration is needed (Liu et al., 2010). Furthermore, the above-mentioned imaging algorithms are dependent on different weighting factors: w p p = cos ⁡ θ , (17A) w p p l a p = − 2 ⁡ cos ⁡ θ ( cos ⁡ θ + 1 ) , (17B)

and w p p p s e − l a p = − ( cos ⁡ θ + 1 ) 2 . (17C)

Here, w p p , w p p l a p , and w p p p s e − l a p are the introduced weighting factor of the scalar imaging condition, the scalar imaging condition with a Laplace filter, and the scalar imaging condition with a pseudo-Laplace filter, respectively. From Eq. 17, we can see that weighting factors vary with the opening angle θ .

The variation of the weighting factor for the opening angle θ is shown in Figure 6. For visual display, the amplitude is normalized by the corresponding max value. When the opening angle θ increases from 0 ∘ to 18 0 ∘ , the weighting factor w p p (represented by the blue curve) in the scalar imaging condition ranges from 1 through 0 to −1. The polarity of the image is reversed while the sign of weight factor changes from positive to negative near 9 0 ∘ , and the backscattering noise is generated by dot product cross-correlation of two wavefields with an opening angle of 18 0 ∘ or close to 18 0 ∘ . By introducing the Laplace filter, the weighting factor w p p l a p (represented by the red curve) will be 0 near 18 0 ∘ , which indicates that backscattering noise has been suppressed. However, the sign of the weighting factor w p p l a p still changes from positive to negative around 9 0 ∘ . In other words, the Laplace filter fails to correct the polarity reversal in PP images by scalar imaging conditions. Furthermore, the pseudo-Laplace filter has been introduced in scalar imaging conditions, and its weighting factor w p p p s e − l a p (represented by the yellow curve) is in the range of 1–0. Similar to the Laplace filter, for backscattering waves with 18 0 ∘ or near 18 0 ∘ opening angles, the weighting factor w p p p s e − l a p is zero. Unlike the Laplace filter, the weighting factor w p p p s e − l a p only ranges from 1 to 0, and its sign is always positive. Therefore, the weighting factor can suppress backscattering noise and correct the reversed polarity. As for the PP image by the scalar imaging condition with a pseudo-Laplace filter, the backscattering noise has been suppressed and polarity reversal has been corrected.

The two-layer flat model shown in Figure 7 is 10 × 2 k m . At a depth of 1 km, there is one flat interface. The first and second layer’s P-wave velocities would be 2400 m/s and 2700 m/s, respectively. The S-wave velocity is consistent with the relationship v s = v p / 1.73 , and density is set to 1.0 g/cm³. It contains 1000 points in the horizontal direction and 200 points in the vertical direction, with a space interval of 10 m. Figure 8 shows a synthetic seismic record generated using double receiving observation geometry, with a shot located at a depth of 10 m. At a depth of 10 m, there are 500 receivers with a 20 m receiver interval. As a result, the maximum offset is up to 5 km. The synthetic seismic data are generated using an explosive source of Ricker wavelet with a peak frequency of 20 Hz. The time interval is 0.8 ms, and the total record time is 2.4 s.

The migrated PP images by scalar imaging conditions, scalar imaging conditions with a Laplace filter, and scalar imaging conditions with a pseudo-Laplace filter are shown in Figure 9. It is evident that backscattering artifacts (marked by the red arrow) influence the PP image by the scalar imaging condition (shown in Figure 9A) and have been attenuated effectively in the PP image with the application of a Laplace filter (shown in Figure 9B) and with a pseudo-Laplace filter (shown in Figure 9C). Apart from backscattering noise, the other noise, such as polarity reversal, also occurs at the interface. As for the interface of 1 km depth, the maximum incident angle reaches 68 . 2 ∘ , over critical angle arcsin ( v p 1 / v p 2 ) = 62.73 ° and polarity reversed angle 4 5 ∘ . Thus, all three migrated PP images are similar and encounter phase-distorted inhomogeneous waves, such as refracted waves (marked by blue arrows) around 1.9 km distance. Besides, only PP images without or with the Laplace filter suffer from the polarity reversal around 1 km distance (marked by red circles), which would have been corrected in the PP image with a pseudo-Laplace filter. Therefore, the phase axis of the PP image with the pseudo-Laplace filter is more continuous than the PP image with the Laplace filter.

Additionally, we measure the amplitudes of these images at the interface at a depth of 1 km and convert the offset variable to the opening angle variable using geological and elastic parameters. Then, we compare these amplitudes (blue curves) from Figure 9 with the theoretical reflection coefficient R p p (yellow curves) and the corresponding weighting factor (red curves), respectively. The analytical solution R p p (yellow curves) is calculated by solving the Zoeppritz equation with the elastic parameters of the two-layer layer model. The opening angle ranges from 0 to 120° to avoid phase distortion when the incidence angle is greater than the critical angle of 62.73°. The extracted amplitudes (blue curves) match well with the weighting theoretical reflection coefficient (red curves) up to approximately 80°, verifying the correctness of the theoretical analysis.

What is more, the extracted amplitudes have higher values than the corresponding weighting theoretical reflection coefficient at large angles of incidence and then decline to zero due to the limited acquisition space. Both Laplace and pseudo-Laplace filters would fail to maintain the amplitude of images at large incidence. Figure 10A shows the extracted amplitudes from Figure 9A and weighting theoretical reflection coefficient R p p ⁡ cos ⁡ θ ; Figure 10B shows the extracted amplitudes from Figure 9B and weighting theoretical reflection coefficient R p p ⁡ cos ⁡ θ ( cos ⁡ θ + 1 ) . They both change their signs at approximately 90° angle of incidence. However, Figure 10C shows the extracted amplitudes from Figure 9C and weighting theoretical reflection coefficient R p p ( cos ⁡ θ + 1 ) 2 , and its sign remains unchanged at any opening angle. Therefore, the consistency of the pseudo-Laplace filter has been demonstrated by the analysis of amplitude variation versus the opening angle, which indicated that polarity reversal had been corrected.

The four-layer inclined model, as shown in Figure 11, is 10 × 4 k m . There are three inclined interfaces with a 1 0 ∘ dip angle at 1.5, 2.5, and 3.5 km depth, respectively. The P-wave velocity of the first layer, second layer, third layer, and fourth layer would be 2500 m/s, 2600 m/s, 2700 m/s, and 2800 m/s, respectively. The S-wave velocity satisfies the relationship v s = v p / 1.73 , and density is set to constant 1 g/cm³. It contains 1000 points and 400 points in the horizontal and vertical directions with a space interval of 10 m. Synthetic data are generated with double receiving observation geometry, where the shot is located at a depth of 10 m and a distance of 7 km. There are 500 receivers with a receiver interval of 20 m at a depth of 10 m. Therefore, the maximum offset is 7 km. The explosive source of the Ricker wavelet with 30 Hz peak frequency is set to generate the synthetic seismic data. The time interval is 1.0 ms, and the total record time is 3 s.

The migrated PP images by scalar imaging conditions with a Laplace filter and scalar imaging conditions with a pseudo-Laplace filter are shown in Figure 12. It is evident that backscattering artifacts (marked by the red arrow) influence the PP image by scalar imaging conditions (shown in Figure 12A) and have been attenuated effectively in the PP image with the application of a Laplace filter (shown in Figure 12B) and with a pseudo-Laplace filter (shown in Figure 12C). Apart from backscattering noise, the other noise, such as polarity reversal, also occurs at images along with the interfaces. As for the first interface, the maximum incident angle reaches 7 5 ∘ , which is over critical angle arcsin ( v p 1 / v p 2 ) = 74.05 ° , and polarity reversed angle 4 5 ∘ at the maximum offset of 7 km. Therefore, the polarity reversal around 5 km (marked by a red circle) and phase-distorted homogeneous wave such as refracted wave (marked by the blue arrow) would be introduced in PP images without or with a Laplace filter. As for the second interface, the maximum incident angle of 5 8 ∘ , equal to 11 6 ∘ opening angle, is less than the critical angle arcsin ( v p 2 / v p 3 ) = 74.35 ° and bigger than the polarity reversed angle of 4 5 ∘ . PP images without or with the Laplace filter of the second interface only encounter a polarity reversal problem around 4 km without phase aberration. As for the third interface, the maximum incident angle of 4 3 ∘ is near the polarity reversed angle of 4 5 ∘ and less than the critical angle arcsin ( v p 2 / v p 3 ) = 74.64 ° . The amplitude of the phase-reversed image is too little to influence the final stacked result. Furthermore, the polarity reversals at three interfaces have been corrected in the PP image with a pseudo-Laplace filter. Therefore, the phase axis of the PP image with a pseudo-Laplace filter is more continuous than the PP image with a Laplace filter.

To analyze the amplitude variation versus opening angle, we pick up the max amplitude of these images along the second interface and convert the offset variable to the opening angle variable with a geological structure. Since the first interface is affected by the direct wave and heterogeneous wave such as refracted wave, the second interface has been utilized. Once the amplitude variation has been picked up, the smoothing and normalization are required to avoid interference of other factors such as phase. As shown in Figure 13, the variations of normalized amplitude in three images match the variations of normalized weighting with an opening angle ranging from 0 to near 11 6 ∘ , the maximum opening angle of the geometry. Once the opening angle exceeds the maximum, normalized amplitudes would be zero. Around 9 0 ∘ opening angle, change the numerical symbols of amplitude that occurs in the PP image (represented by the blue curve) and the PP image with a Laplace filter (represented by the red curve). However, the yellow curve, representing the PP image with a pseudo-Laplace filter, has a consistent numerical symbol in amplitude. Therefore, the consistency of the pseudo-Laplace filter has been demonstrated by the analysis of amplitude variation versus the opening angle, which indicated that polarity reversal had been corrected.

The Marmousi 2 model is used in this example to show how the pseudo-Laplace filter can effectively suppress backscattering noise, resolve polarity reversal, and generate a high-quality PP image in complex geological structures.

As shown in Figure 14, the modified Marmousi 2 model (Martin et al., 2006) contains 1325 points in the horizontal direction with 10 m sample interval and 934 points in the vertical direction with 5 m sample interval. Thus, the size of the Marmousi 2 model is 13.25 × 4.67 k m . The P-wave velocity ranges from 1800 m/s to 4600 m/s, the S-wave velocity ranges from 1000 m/s to 2000 m/s, and density is the constant of 1 g/cm³. The full-receiving geometry, where 265 shots are excited with a 50 m shot interval at a depth of 10 m, and 1325 receivers are located with a 10 m receiver interval at the surface, is constructed to generate the seismic record. The source function is a Ricker wavelet with a peak frequency of 30 Hz. The recording time interval is 1.0 ms, and the recording length is 4.3 s.

Figure 15A is the migrated PP image by the scalar imaging condition with the Laplace filter, and Figure 15B is the migrated PP image by the scalar imaging condition with a pseudo-Laplace filter. The backscattering noise contamination has been suppressed successfully in these two images, while some backscattering noise still resides in the PP image with Laplace (marked by the red arrow). Furthermore, both of them have embodied the structural character of the model. Compared with the PP image with Laplace of Figure 15A, the overall appearance of Figure 15B is more apparent, including the shallow fault in detail. That is related to the weighting factor of the pseudo-Laplace filter is stronger than that of the Laplace filter. Furthermore, the events in the shallow layer of Figure 15B are more continuous than those of Figure 15A.

We further extracted the traces from PP images with a Laplace filter (Figure 15A) and a pseudo-Laplace filter (Figure 15B) at a depth of 1.32 km. As shown in Figure 16, we compare the amplitude of trace extracted from Figure 15A (blue curve) and Figure 15B (red curve) with the theoretical reflection coefficient R p p (yellow curves). The Zoeppritz equation calculates the theoretical reflection coefficient R p p at normal incidence. The variation trends of traces are consistent with the theoretical reflection coefficient. Moreover, the amplitude of trace from I p p p s e _ l a p (red curve) is generally slightly larger than that from I p p l a p (blue curve) because the weighting factor 2 ⁡ cos ⁡ θ ( cos ⁡ θ + 1 ) of PP images with a Laplace filter is stronger than the weighting factor ( cos ⁡ θ + 1 ) 2 of PP images with the pseudo-Laplace filter at small incident angles. However, there is a big difference between amplitudes and theoretical value in the fault area at 5–6 km and 7–8 km. The inaccuracy is caused by a false reflected wave where diffraction wave and multiple waves exist. Moreover, at the cover of anticline where sub-cover oil and gas reservoirs develop with 10–11 km distance, the phase of stacked I p p l a p (blue curve) is opposite to that of stacked I p p p s e _ l a p (red curve) and theoretical solution (yellow curve), as marked by a black rectangular box. The phase reversal originates from polarity reversal in the PP image with the Laplace filter and polarity correction in the PP image with the pseudo-Laplace filter at large incidence.

To observe clearly, we intercept the part of the Marmousi 2 model with a distance from 7–13 km and a depth from 0.5 to 4 km. Figures 17A,B are the partial enlargements of the P-wave velocity of the elastic Marmousi 2 model and S-wave velocity of the elastic Marmousi 2 model amplified, respectively. Correspondingly, Figure 18A is a partial enlargement of the PP image of the Marmousi 2 model by the scalar imaging condition with a Laplace filter (Figure 15A), and Figure 18B is a partial enlargement of the PP image of the Marmousi 2 model by the scalar imaging condition with the pseudo-Laplace filter (Figure 15B). Even as the polarity reversal described by Figure 16, the phase of the stacked PP image with a Laplace filter is opposite to that of the stacked PP image with a pseudo-Laplace filter and is no longer continuous at the cover of the anticline. As marked by the blue curve, the phase of the event in Figure 18B is more persistent than in Figure 18A, especially at the cover of an anticline. Furthermore, there is an oil and gas reservoir at sub-cover with the variation of S-wave velocity in Figure 17B. The disturbance at sub-cover (marked by the red arrow) succeeds to be imaged in Figure 18A but fails to be imaged in Figure 18A. Overall, the scalar imaging condition with a pseudo-Laplace filter generates a high-quality PP image in complicated geological structures.