Point-spread function convolution to simulate prestack depth migrated images: A validation study

Seismic migration commonly yields an incomplete reconstruction of the Earth model due to restricted survey aperture, band-limited frequency content and propagation effects. This affects both illumination and resolution of the structures of interest. Through the application of spatial convolution operators commonly referred to as point-spread functions, simulated prestack depth-migrated images incorporating these effects may be obtained. Such simulated images are tailored for analysing dis-tortion effects and enhance our understanding of seismic imaging and subsequent interpretation. Target-oriented point-spread functions may be obtained through a va-riety of waveform and ray-based approaches. Waveform approaches are generally more robust, but the computational cost involved may be prohibitive. Ray-based approaches, on the other hand, allow for efficient and flexible sensitivity studies at a low computational cost, but inherent limitations may lead to less accuracy. To yield more insight into the similarities and differences between point-spread functions obtained via these two approaches, we first derive analytical expressions of both wave- and ray-based point-spread functions in homogeneous media. By considering single-point scatterers embedded in a uniform velocity field, we demonstrate the conditions under which the derived equations diverge. The accuracy of wave-based and ray-based point-spread functions is further assessed and validated at selected targets in a subsection of the complex BP Statics Benchmark model. We also compare our simulated prestack depth migrated images with the output obtained from an actual prestack depth migration (reverse time migration). Our results reveal that both the wave- and ray-based approaches accurately model illumination, resolution and amplitude effects observed in the reverse time-migrated image. Furthermore, although some minor deviations between the wave-based and ray-based approaches are observed, the overall results indicate that both approaches can be used also for complex models.

lost. An understanding of how these blurring and illumination effects occur during seismic imaging is important for proper interpretation of migrated seismic images.
For a given reflectivity model, m, we may compute synthetic seismic traces, d, via the equation d = Lm, where L is an applied forward modelling operator. A migrated image may then be obtained from the equation m mig = L T Lm (Schuster and Hu, 2000), with L T representing an approximation of the inverse L-operator. The blurring effects observed in migrated seismic data are now mathematically expressed through the L T L-operator, commonly referred to as the local point-spread functions (PSFs). Information about the local illumination, as well as the overall resolution, may be retrieved by analysing the PSFs, i.e. the point-scatterer response of the seismic imaging system (Beylkin et al., 1985;Lecomte and Gelius, 1998;Gelius et al., 2002a). This type of information can, in turn, be analysed for assessing how the choice of, e.g., survey geometry, source wavelet, migration aperture, imaging condition, amplitude corrections and other processing parameters affects the seismic image. When PSFs are used as 2(3)D convolution operators applied to an incident angle-dependent input reflectivity grid, simulated prestack depth-migrated (PSDM) images incorporating these resolution and illumination effects are obtained Toxopeus et al., 2003;Lecomte, 2008;Toxopeus et al., 2008). Thus, realistic seismic images can be easily and rapidly simulated based on possibly complex geomodels .
In addition to allowing for simulation of migration effects, accurate target-oriented PSFs may be used for other applications as well. Several studies demonstrate how the local PSF, or an estimate of its inverse, can be used for targetoriented migration and inversion (e.g. Schuster and Hu, 2000;Sjoeberg et al., 2003;Guitton, 2004;Xie et al., 2005;Valenciano et al., 2006;Yu et al., 2006;Aoki and Schuster, 2009;Tang, 2009;Ayeni and Biondi, 2010;Zhao and Sen, 2018;Jiang and Zhang, 2019). Fehler et al. (2005) propose to use PSFs to analyse how the choice of the migration operator affects image resolution. Thomson et al. (2016) demonstrate how the blurring function may be expressed as extended image gathers. This approach may be used to assess how blurring varies with the incidence slowness vector and reflectivity angle. Lecerf and Besselievre (2018) show how PSFs may be successfully applied in 4D time-lapse reservoir monitoring. Fast and efficient methods for the estimation of accurate targetoriented PSFs are therefore of ongoing interest.
Target-oriented PSFs are frequently estimated by various wave-equation-based approaches. Xie et al. (2005) illustrate how target-oriented PSFs may be obtained using a one-way wave-equation-based propagator to downward extrapolate the source functions and the receiver wavefields. Toxopeus et al. (2008) estimate PSFs by applying a one-way wave equation as previously outlined in Thorbecke et al. (2004). They demonstrate how these PSFs may be useful for both simulating defocusing in migration images and as inputs for inversion algorithms. Tang (2009) illustrates how PSFs may be retrieved using a one-way wave-equation-based Fourier finite-difference migration (Ristow and Rühl, 1994) and then used to approximate the inverse Hessian operator for targetoriented inversion.
As an alternative to wave-equation-based approaches, Gelius et al. (1991), Hamran and Lecomte (1993) and Lecomte and Gelius (1998) describe how target-oriented PSFs may be estimated at a much lower computational cost with the use of ray theory. Both ray-based and wave-equationbased approaches for estimating PSFs have their advantages and drawbacks. As for seismic modelling in general (see e.g. Carcione et al. 2002 for an in-depth review of advantages and disadvantages of various seismic modelling approaches), the choice of method is often made as a trade-off between accuracy and computing cost, although, in some cases of highly complex geology, ray-based approaches may not work properly. Ray-based methods generally suffer from drawbacks inherent in the high-frequency approximation used in ray tracing, such as incomplete wavefield, the requirement of a smooth velocity field, multipathing, the risk of breakdown in singular, caustic regions, etc. (Červený et al., 1977). However, ray approaches are often more robust than assumed (Gjøystdal et al., 2007), especially when applied to the rather smooth velocity models used for migration. Their low computational cost allows for fast and flexible modelling of how various migration parameters cause blurring and illumination effects in PSDM images (Lecomte et al., 2016;Jensen et al., 2018).
Although ray-based PSF convolution modelling has been applied in several geological case studies (e.g. Botter et al., 2014;Lecomte et al., 2015;Botter et al., 2017;Kjoberg et al., 2017;Eide et al., 2018;Rabbel et al., 2018;Grippa et al., 2019;Lubrano-Lavadera et al., 2019;Jensen et al., 2021), comparative studies between waveform-based modelled/migrated data and simulated seismic images obtained via ray-based and wave-based PSF convolution modelling are limited. Lecomte et al. (2003) illustrate how PSF-convolved images obtained from ray tracing match well with Kirchhoff PSDM images obtained from a model of the Gullfaks field located on the Norwegian Continental Shelf. Toxopeus et al. (2008) compare results obtained via wave-equationgenerated PSFs with actual time-migrated seismic images of sulphate dissolution and karst collapse-related deformation. The authors demonstrate an improvement in the modelled results compared with results obtained via repeated 1D convolution. Amini et al. (2020) compare prestack Kirchhoff time-migrated results with modelled results obtained after applying a lateral smoothing function aiming at reproducing lateral smoothing effects (Chen and Schuster, 1999), to repeated 1D convolved results. This approach, as illustrated by the authors, yields more similar migrated and modelled images. Yet, a more thorough analysis of how PSF-convolved images obtained via different approaches compare to actual PSDM images is needed. This applies particularly for ray-based approaches where few comparative studies exist.
In the present work, we start by deriving analytical expressions for wave-based and ray-based PSFs valid for homogeneous media. By considering the case of single-point scatterers embedded in a uniform background, analytical modelling based on either wave-or ray-based implementations can be compared with output from an actual PSDM. Possible limitations inherent in the ray-based approach are identified and discussed. We then employ a subsection of the complex BP Statics Benchmark model (Ellison and Innanen, 2016) and assess how wave-based and ray-based PSF-convolved seismic images at selected targets compare to the output from a full PSDM of the subsection obtained via reverse time migration (RTM). Similarities and differences between the results obtained from the two modelling approaches are also assessed and discussed. The results may add to our understanding of the potential applications, as well as possible limitations, of using wave-based and ray-based PSF convolution as an approach for simulating PSDM images.

T H E P O I N T -S P R E A D F U N C T I O N I N A H O M O G E N E O U S M E D I U M
A migrated image, m mig , is related to a model quantity, m, via the equation: where L is a forward modelling operator, and L −1 is a stable approximation of its inverse. Assuming the Born approximation holds, the modelling operator, L, may, for a single angular frequency value, be defined as where ω is angular frequency, S(ω) is the wavelet spectrum, G(ω, r |r s ) represents the Green's function from the source at Figure 1 Survey geometry and relevant parameters for a single shotreceiver pair in a homogeneous model r s to the image point at r , and G(ω, r g |r ) the Green's function from the image point at r to the receiver at r g . The Green's functions are, following the Born approximation, computed in the background media. A complete derivation of the Born approximation may be found in e.g. Bleistein et al. (2001) and Schuster (2017). We now consider a homogeneous medium with an image point of interest at r and a neighbouring point at r . Further, we assume a source located at position r s , and a receiver at position r g . The medium has a constant background velocity of c 0 . We denote the distance from r s to r as R s (r) = |r − r s |, and the distance from r g to r as R g (r) = |r − r g |. Figure 1 illustrates the survey geometry and the involved parameters.
Using these initial assumptions, and the survey geometry in Figure 1, we may now derive analytical expressions for the PSF generated at the point of interest in the homogeneous medium.

Wave-equation approach
Within a wave-equation approach, the inverse operator L −1 is usually approximated through the adjoint operator, L T . The PSF is then expressed via the operator L T L, which accounts for how limited illumination, overburden propagation and frequency bandwidth blur the migrated image. The migration equation may now, in the model space, be expressed as Here m(r ) represents the perturbed slowness at the point r : with c 0 representing the reference, or background, velocity, and c representing the actual medium velocity at r . PSF wave (r|r ) is the PSF which, for a single shot and using the definition for L in (2), is defined as (Schuster, 2017, p. 211) Here, a star represents the complex conjugate operation, and this equation can be interpreted as the migration response at point r due to a point scatterer at r'. The multiplication of the downgoing source field with the back-propagated reflections in (5) represents a standard cross-correlation imaging condition (Claerbout, 1971).
For a homogeneous 3D medium, we may define the farfield Green's functions analytically by where k is the wavenumber and R is the distance between two points of interest. Inserted into (5), and using the parameter definitions from Figure 1, as well as the relation k = ω/c 0 , we obtain the analytical expression for the PSF for a single shot and a point scatterer in a homogeneous 3D medium: For a homogeneous 2D medium, the Green's function is defined analytically by where k and R represent the same parameters as in (6), and H 1 0 (kR) is the zeroth order Hankel function of the first kind. Assuming a far-field approximation (for kR 1), this function is approximated as Through (8) and (9), we then obtain the analytical expression for the PSF for a single shot and a point scatterer in a homogeneous 2D medium: We now proceed to derive the corresponding ray-based approach.

Ray-based approach
In our ray-based approach, we take inspiration from the local imaging method derived from local plane-wavenumber diffraction tomography (Hamran and Lecomte, 1993). We consider a small volume V 0 around the image point of interest. We assume that the Born approximation holds in this region, and that the background slowness is locally homogeneous. Moreover, we assume that the incident and scattered wavefields are plane within the volume (Hamran and Lecomte, 1993;Gelius et al., 2002a). We let r denote the location of an image point of interest, and r a reference point within the same volume.
Using a ray-based approach, we may define the Green's functions via the high-frequency approximation: G ω, r j |r = A ω, r j |r exp iωτ r j |r , j = s, g, where A(ω, r j |r) represents the complex amplitude, and τ is the traveltime between a source (or receiver) and the image point. Assuming constant amplitude and a linear phase within the local volume V 0 , we can Taylor expand the Green's function around the reference point, r : G ω, r j |r ∼ = G ω, r j |r exp iω∇τ r j |r (r − r) .
We now introduce the scattering wavenumber vector K defined by K = −ω ∇τ s (r|r s ) + ∇τ g r g |r = −k s + k g , where τ s is the traveltime from source point to r, and τ g is the traveltime from r to receiver point. It can be shown that K represents the Fourier vector of the model space (Gelius, 1995), meaning the Fourier transform of the model space coordinates. This means that for each shot-receiver combination over the model space, a scattering wavenumber vector can be defined in the wavenumber domain. The entire collection of scattering wavenumber vectors will then yield the local wavenumber domain representation of the PSF. Thus, the expression in (13) links the model and acquisition domains.
We assume now that the phase of the wavefield varies significantly faster than the amplitude changes within the local volume of interest. At the same time, unlike in (10), we ensure that the PSF is normalized with respect to the reflectivity. The following ray-based approximation of the inverse of the modelling operator L may then be derived (Gelius et al., 1991;Gelius, 1995): Here, J[K|ω, r g ] represents a Jacobian, derived in full in the Appendix, that maps between the Fourier-transformed model (scattering wavenumber) domain and the acquisition domain illustrated in Figure 1. In the general case, all quantities needed to compute it can be obtained from dynamic ray tracing. The inverse operator expressed in (14) is obtained from the use of a deconvolution imaging condition (Claerbout, 1971). Through (11)-(14) and the definition of the forward modelling operator L from (2), we obtain an expression for the ray-based PSF, under a Born-scattering assumption, valid for both 3D and 2D media (Gelius et al., 2002a): where dK = J(K|ω, r g )dωdr g . Here, the PSF is represented as a sum over all the scattering wavenumber vectors at the image point of interest, and the source spectrum |S(ω)| is included as an obliquity factor for the Born case. The equation represents an inverse Fourier transform, yielding a real-valued PSF in the spatial domain.
To allow for direct comparisons between wave-generated and ray-generated PSFs, we consider the analytical expression for the ray-based PSFs valid for a single shot and a point scatterer in a homogeneous 2D medium. By inserting the parameters defined in Figure 1 into Equation (15), the following analytical solution, derived in full in the Appendix, is obtained: where

Validation in the homogeneous case
To validate the derived equations, we consider two targets (targets 1 and 2) in the homogeneous 2D model illustrated in Figure 2. Two versions of the model were used, one with a constant velocity of 2 km/s (low-velocity model) and one with a constant velocity of 4 km/s (high-velocity model). We use two versions to demonstrate how the background velocity affects the PSF estimation. Using a 2D acoustic finite-difference (FD) approach on a regular grid, we performed forward modelling was performed for a single shot and a fixed receiver array located along the top row (Fig. 2). To implement the FD code, the domain was discretized with a second-to fourth-order scheme (second-order in time and fourth-order in space) following the numerical implementation presented in Youzwishen and Margrave (1999). A standard zero-phase Ricker wavelet with a peak frequency of 10 Hz was used as a source pulse, and the temporal sample interval was 1 ms. Although a combination of high velocity and low-frequency wavelet is geologically unrealistic for shallow targets, which we get for target 2 in the v = 4 km/s case, we deliberately chose these parameters to illustrate where the divergence between the wave-based and ray-based approaches occurs from a theoretical point. Half-spaces were added to all edges of the models to avoid reflections from the edges. Finally, a reverse time migration (RTM) was performed on the seismograms to obtain the responses from the point scatterers.
The obtained PSFs were extracted at both targets from the migrated images using a window of size 0.4 km times 0.4 km (41 × 41 grid points), with the PSF centred at the middle of the window. Correspondingly, PSFs for both targets were calculated from Equations (10) and (16) using the same  (10) and analytic implementation of equation (16) for target 2. (j-l) Corresponding wavenumber spectra for the PSFs survey parameters. Following the cross-correlation imaging condition inherent in RTM, an extra |S(ω)| factor was included in the numerator of Equation (16) to obtain the |S(ω)| 2 term. All PSFs were further normalized to have a maximum amplitude of one to allow for relative comparison of the seismic responses. The PSFs, and their corresponding wavenumber spectra, are plotted in Figure 3 (low-velocity model) and Figure 4 (high-velocity model). Figure 5 illustrates a comparison of the centre traces obtained for all PSFs, both vertically and horizontally.
For target 1, all PSFs (Fig. 3a-c and Fig. 4a-c) have the same resolution patterns. The same observation applies in the corresponding wavenumber domains (Fig. 3d-f and Fig. 4d-f), where the spectra match in both overall coverage and in the amplitude range. We do, however, notice that the ray-based PSFs diverge slightly at the shallower target 2 (Figs 3i, 3l, 4i and 4l). This is further illustrated in Figure 5. Here, we notice how for target 1, the centre traces from the analytically computed PSFs align almost perfectly with the traces obtained from the PSF extracted from a complete migration. For target 2, we notice slightly more deviations between the ray-based PSF traces and the others, both in the vertical and horizon-tal directions. This is particularly seen for the shallow PSFs computed in the high-velocity model ( Fig. 5d and h).
To understand the observed differences, we consider Equations (10) and (16) again. First, we keep in mind that both equations are derived under a far-field approximation. Some deviances from the PSFs extracted from RTM are therefore to be expected for both approaches, particularly at the shallower target 2. To quantify this issue, we computed the PSF at 49 equidistant grid points in the low-velocity model using the Green's function defined in (8) with full Hankel functions rather than their asymptotic approximation. At all 49 targets, the absolute error spectra obtained from subtracting the PSFs obtained with the Green's function defined in (8) from the PSFs obtained under the far-field assumption in (10) were calculated. The standard deviation was then computed for each error spectrum, and 2D linear interpolation was applied to obtain expected standard deviation values for the entire model. Figure 6 illustrates the obtained result. We notice that the deviation is largest close to the shot point, as expected considering that the criteria kR 1 might not be fulfilled here. We do, however, also notice that the standard deviation values are close to zero for most parts of the model, meaning that the  (10) and analytic implementation of equation (16) for target 1. (d-f) Corresponding wavenumber spectra for the PSFs. (g-i) PSFs obtained from RTM, analytic implementation of equation (10) and analytic implementation of equation (16) for target 2. (j-l) Corresponding wavenumber spectra for the PSFs far-field approximation should be valid throughout most of the model.
Regarding the two approaches themselves, we note that the phase terms are identical in both equations. Differences between the two approaches thus depend on the amplitude factors. We notice that in the ray approach, as expressed by (16), we only consider the distance to the image point at r, in question, and not the surrounding scattering points located at r . In essence, the amplitude term in (16) may be considered as a simple scaling factor valid throughout the entire target of interest. This follows from (14) where we deconvolve the amplitude factors in accordance with the plane-wave assumption. The amplitude factors in (16) are therefore not attributable to the Green's functions, but merely result from the Jacobian mapping between the scattering wavenumber and acquisition domains.
We illustrate this effect by again considering the two targets in the low-velocity model. Using the analytical expression in (10), we computed the PSF amplitude responses (assuming a constant phase of value one) obtained at the point scatterers. For simplicity, we let ω = 1 and set |S(ω)| = 1 during the computations to extract amplitude variations only. By selecting the same PSF window of size 0.4 km × 0.4 km and normalizing the centre point of the PSFs to a value of one in both cases, we obtained the amplitude responses presented in Figure 7. The corresponding amplitude responses for the ray-based approach obtained via (16) would simply, following the planewave assumption, yield a constant value throughout the entire PSF window.
As amplitudes vary within the target area for the wavebased approach, slight asymmetry around the centre point along the high-resolution axis (in this case the vertical axis) is expected for traces obtained from the wave-based PSF. We also notice that the amplitude range is larger for the shallow target (Fig. 7b). This also follows from (10) as a smaller distance between shot point and target point yields smaller values for the r and r parameters in the denominator in (10). The asymmetry is usually negligible if the PSF operator only spans a small amount of grid points around the centre point. However, for low-resolution PSFs resulting from low-frequency bandwidth, and/or high velocity at the target area, the effect would be more noticeable. This is for instance observed in Figures 4(h) and 5(d), where the combination of small distance between shot point and target point, high velocity, and low wavelet frequency, yields a noticeable asymmetry around the PSF centre point for the wave-based PSF. Furthermore, the plane-wave Standard deviation of absolute error spectra obtained between normalized PSFs estimated with and without the far-field approximation expressed in equation (9) at 49 equidistant points followed by 2D interpolation. The computations were performed in the homogeneous model with constant velocity of v = 2 km/s assumption inherent in the far-field ray-based approach may produce differences in the PSF side lobes compared with wavebased PSFs where curvature effects are preserved. These effects should also be more prominent for shallow targets.
To illustrate in greater detail how the wave-generated PSFs may diverge from ray-generated PSFs in the homogeneous case, we estimated amplitude ranges as illustrated in Figure 7 for normalized PSFs obtained via (10) at 81 equidistant grid points in the low-velocity model. The stan-  Figure 8. The results confirm that the divergence between the two approaches will decrease with increasing distance from shot point in a homogeneous model.
In summary, divergence between wave-based and raybased approaches is attributable to the ray-based local planewave assumption in amplitude, normalization with respect to reflectivity, the far-field approximation and the Jacobian mapping taking into account irregular sampling of the model space. The divergence is mostly observed for shallow targets characterized by high-velocity and/or low-frequency bandwidth. However, as such a combination is, as pointed out, geologically unrealistic, we should, within the validity of ray theory and even for complex models, not expect much divergence between PSFs regardless of whether they are estimated with a wave-or ray-based approach.

E S T I M AT I N G P O I N T -S P R E A D F U N C T I O N S I N I N H O M O G E N E O U S M O D E L S
In a homogeneous model, it is straightforward to analytically implement (10) and (16) and compute the seismic response due to a point scatterer. However, the Green's functions are not so easily computed for more complex models and must be estimated via other methods.
When applying a wave-based approach, one option is to estimate the PSFs by implementing the wave equation to retrieve the required Green's functions in (5). The true velocity model is then used for estimating the forward modelled Green's functions, and the-often smooth-migration velocity model is used for estimating the back-propagation . An alternative approach is to perturbate the velocity or density value at a single point in the smooth velocity model, followed by forward modelling and migration over the model. The extracted point-scatterer response would then yield the PSF (Cao, 2013).
If a ray-based approach is implemented, we may, following Lecomte and Gelius (1998) and Lecomte (2008), estimate the scattering wavenumber vectors in (13) from ray tracing or similar (e.g. eikonal solvers) in a smooth background velocity model (as done in actual migration). First, a target point is selected for PSF computation. With a plane-wave assump-tion, ray tracing is then used to calculate the incident and scattered slowness vectors (p S and p R ) at the target point for each shot-receiver combination in a seismic survey. The so-called illumination vectors, I SR = p R − p S , are then computed for all these combinations (Fig. 9a). It can be shown mathematically that any reflector perpendicular to an illumination vector will be well illuminated (Gelius et al., 2002a). As such, the fan of resulting illumination vectors obtained from all shotreceiver combinations in the survey (Fig. 9b) contains information about the local geological-dip range which may be imaged at the point of interest. Furthermore, the local acrossreflector resolution depends on the magnitude of each illumination vector, which is estimated as a function of opening angle and medium velocity via (Lecomte, 2008): where θ SR is the opening angle between the incident and scattered wavefield, c is the velocity of the incoming and scattered wavefields at the image point, and u SR is a unit vector pointing in the direction of I SR . By mapping the properly weighted wavelet spectrum (e.g. squared in the case of a cross-correlation imaging condition following equations (7) and (10)) along each illumination vector, scattering wavenumber vectors as defined in (13) are obtained for the target point (Fig. 9c). The entire collection of scattering wavenumber vectors now represents the PSF in the wavenumber domain. A Fourier transform of the spectrum yields the PSF in the space domain (Fig. 9d). An important point to consider is that the procedure illustrated in Figure 9 assumes that the wavenumber spectrum is defined directly on a regular grid in the wavenumber domain. Intrinsically, however, scattering wavenumber vectors are best mapped in the polar domain in 2D (spherical in 3D). But a direct implementation of a 2D FFT in such a domain, followed by a polar-to-Cartesian transformation, is cumbersome because seismic data are typically characterized in the polar domain by very irregular sampling of the angular coordinate (corresponding to irregular -and possibly lacking -illumination). As such, various non-uniform interpolation and resampling strategies are used when converting raw data from polar coordinates to Cartesian coordinates. This is seen in, e.g., synthetic aperture radar (SAR) imaging (e.g. Carrara et al., 1995;Jakowatz et al., 1996;Doerry, 2012). Note, however, that the mentioned studies contain far less irregular sampling than in seismic data as these studies deal with a homogeneous background velocity field (i.e. the air). In our approach, a nearest-point interpolation technique is simply used. This technique is easy to implement and highly efficient: each polar sample is assigned to the nearest sample of the selected regular grid. However, such a mapping of points along vectors all attached to the same origin or pole (here, the reference point for which the PSF is calculated) yields an irregular hit count, i.e. with a higher number of hits close to the origin (where K x = K z = 0) than further away. To compensate for this effect, the value at each grid cell is normalized by dividing it with the number of hits at that cell (Lecomte et al., 2005). This simple and efficient procedure both compensates for the polar-to-Cartesian mapping (Jacobian) and a possibly irregular illumination, including local redundancy leading to coherent noise in the migrated image (see example in Lecomte et al., 2005 ; Fig. 16). The simulated prestack depth-migrated (PSDM) image obtained this way is, in essence, representing a 'perfect' PSDM image where the illumination has been regularized, as it should. If the goal, however, is to simulate PSDM images obtained from a migration algorithm with a different imaging condition which does not properly compensate for irregular illumination, in particular in the case of redundant illumination, the mapping may first be done directly in the polar domain, where one can easily control range and sampling while keeping the actual illumination hit. The final mapping over to the Cartesian domain could then be done by interpola-tion after accounting for aliasing via, e.g., adequate smoothing (Lecomte et al., 2005).

CA S E S T U DY: B P S TAT I C S B E N C H M A R K M O D E L
To assess the validity of simulated prestack depth-migrated (PSDM) images obtained via point-spread function (PSF) convolution modelling, we consider a slightly modified subsection of the 1994 BP Statics Benchmark P-velocity model (Ellison and Innanen, 2016). The velocity model is illustrated in Figure 10(a) with three selected target areas highlighted. The subsection is dominated by a large high-velocity intrusion in a complex geological setting characterized by uneven layers and faults. A thin, homogeneous water layer with constant P-velocity of 1.5 km/s was added on top of the model for simulation of a marine-type survey. Model and survey parameters are defined in the Figure 10 caption.
Using a zero-phase Ricker wavelet with a peak frequency of 20 Hz, sampled at 1 ms, finite-difference forward modelling was performed using the same approach as in the homogeneous case. Similarly, half-spaces were added at all four boundaries to avoid unwanted boundary reflections. Next, RTM was applied on the forward modelled traces. For the migration itself, a velocity model smoothed over slowness was Wave-based PSFs and corresponding wavenumber spectra for PSFs generated at the targets highlighted with white boxes in Figure 10(a). Wavenumber spectra are plotted for K x between −32.9 and +32.9 km −1 , and for K z between −37.8 and −3.7 km −1 applied (Fig. 10b). The final image obtained after RTM is illustrated in Figure 10(c). A phase-shift operator was applied to the migrated traces to obtain approximate zero-phase traces.
Following the procedure in Cao (2013), wave-based PSFs were estimated at the centre of the three target areas by adding point scatterers to the smooth velocity model presented in Figure 10(b) and extracting their responses after forward modelling and RTM. All wave-based PSFs, with corresponding wavenumber spectra, are illustrated in Figure 11. For target 3, an additional PSF was estimated just below the boundary of the high-velocity intrusion. This was done as targets characterized by high-velocity contrasts and complex geology may need to account for space-variant PSFs for accurate modelling results, as both the illumination and across-reflector resolution may vary substantially throughout the target area. For target 3, we therefore illustrate the different results obtained Figure 12 Ray-based PSFs and corresponding wavenumber spectra for PSFs generated at the targets highlighted with white boxes in Figure 10(a). Wavenumber spectra are plotted for K x between −32.9 and +32.9 km −1 , and for K z between −37.8 and −3.7 km −1 when employing a single PSF versus employing two PSFs. In the combined image, the centre PSF was convolved with the reflectivity grid above the high-velocity intrusion, and the other PSF was convolved with the reflectivity grid at and below the intrusion boundary. For simplicity, the two modelled images were then simply added together to yield the final image.
Using a wavefront construction approach presented in Vinje et al. (1993), ray tracing was next performed in the same smooth velocity model in order to estimate ray-based PSFs at the centre of each target area via the procedure outlined in Figure 9. For target 3, as we did for the wave-based approach, an additional PSF was estimated below the intrusion boundary. P-wave reflections were estimated based on the average reflectivity and average angle range values obtained from ray tracing. This was done for simplicity, as such an approach only requires one modelling run instead of multiple modelling runs for combination of angle ranges. All ray-based PSFs, with corresponding wavenumber spectra, are illustrated in Figure 12.
For each target, the seismic amplitude values were normalized by extracting the amplitude value from the same point at a well-imaged reflector and setting this value equal in the migrated and modelled images. The results obtained at all targets are illustrated in Figures 13-16. We will now discuss the overall findings with focus on how well the modelled images capture illumination, resolution and amplitude effects observed in the PSDM image.

Illumination effects
Targets 1 and 2 are located in well-illuminated parts of the model. The migrated results (Figs 13a and 14a) capture all the geological features observed in the reflectivity grids (Figs 13b and Fig. 14b). We further observe that the modelled results ( Fig. 13c-d and Fig. 14c-d) match the migrated results well. For target 2, however, there is a slight dimming of the modelled seismic response at the steepest part of the upper reflector (yellow box in Fig. 14c and d). By comparing the traces obtained at x = 2.15 km (Fig. 14e), we do observe in more detail how the amplitudes of the modelled traces weaken at depths around z = 0.66 km -0.7 km compared with the migrated trace. This could result from minor changes in illumination affecting different parts of the target area.
Illumination effects are also captured well by the modelled images for target 3, whether through the use of one PSF (Fig. 15) or two PSFs (Fig. 16). Target 3 is characterized by two primary reflectors, with a fault crossing the layer between the reflectors. The fault is not imaged in the migrated image ( Fig. 15a) due to lack of illumination. This is also captured through PSF convolution modelling using both wave-based and ray-based PSFs (Fig. 15c-d and Fig. 16c-d). As such, this illustrates how PSF convolution modelling may properly account for limited illumination of dipping geological features such as faults, even when the latter have an elastic impedance contrast across.

Resolution effects
Differences in resolution between the migrated and modelled images are primarily observed at target 3, which is characterized by high-velocity contrasts. No discernible differences in resolution are observed for targets 1 and 2, but for target 1 we notice that the wave-based traces plotted in Figure 13(e, f)  Fig. 13b). (f) Traces obtained at x = 2.73 km from migrated and modelled results (yellow line marked (f) in Fig. 13b). The horizontal dashed lines in (e) and (f) represent exact reflector positions align slightly closer to the traces obtained from the migrated image than the ray-based traces do. As illustrated in Figure 8 for the homogeneous case, this might be a result of the shallowness of the target. The lack of curvature effects captured by the ray-based PSFs may also yield slightly different responses. At target 2, the traces obtained from wave-based and raybased convolution (Fig. 14e, f) are almost identical and overall align better than for the shallower target 1, indicating that the greater target depth leads to less divergence between the wave-based and ray-based PSFs.
For target 3, however, we observe that the lower reflector, representing the intrusion boundary, suffers from slightly worse resolution in the migrated image (Fig. 15a) than in the modelled images (Fig. 15c-d). This is not surprising as the single PSFs used for convolution modelling in Figure 15(c, d) were estimated at the centre of the target area, located above the intrusion boundary. The lower velocity here yields PSFs with greater resolution than a PSF estimated at, or below, the boundary. By comparing traces extracted at x = 1.61 km (Fig. 15f), we do notice how the resolution of the migrated trace deviates from the modelled traces at the lower reflector. In order to mitigate this issue, multiple PSFs may be used to capture the spatial variability of the PSFs at the target. As illustrated in Figure 16(c, d), we indeed observe that when a sec-ond PSF is estimated below the intrusion boundary and used for convolution with the lower part of the reflectivity grid, the combined image obtained from the two PSF convolutions better capture the resolution observed in the migrated image. This is also evident if we compare the traces plotted in Figure 16(f). What is of note is that the traces obtained from wave-based and ray-based PSFs overlap almost perfectly (Fig. 16e, f), indicating that both methods yield approximately identical results at this target.

Amplitude effects
The use of a single PSF convolution operator may, in some cases, not fully capture amplitude effects caused by transmission loss. For target 1, we notice that the loss of energy of the transmitted wave across the main reflector results in slightly greater amplitude contrast between the reflectors in the migrated image (Fig. 13a) compared with the modelled images ( Fig. 13c-d). As the modelled images were obtained by convolving only one PSF with the input target reflectivity grid, this transmission effect is not accounted for. This is also observed in the traces obtained at x = 2.45 km (Fig. 13e). We observe here that the modelled traces, whether a wave-based or ray-based PSF is applied, indeed deviate slightly from the  Fig. 14b). (f) Traces obtained at x = 2.33 km from migrated and modelled results (yellow line marked (f) in Fig. 14b). The horizontal dashed lines in (e) and (f) represent exact reflector positions migrated trace in terms of amplitude. The same phenomenon is observed for target 3 and is particularly evident in the trace obtained at x = 1.61 km (Fig. 15f). When two PSFs are used, this problem is mitigated (Fig. 16f).
Amplitude differences between migrated and modelled images could also, in general, be caused by migration artefacts, interference issues and/or wave-energy originating from other parts of the model outside the considered target area. Such effects may not be fully captured by the (local) PSF convolution operators. The amplitude dimming observed in the part of target 2 highlighted by the yellow box in Figure 14(c, d) could, in addition to illumination sensitivity at the target, be explained by this phenomenon. Sharp boundaries in areas with high-velocity contrasts are prone to migration artefacts, and we do also notice for target 3 that the migrated image (Fig. 15a) appears to suffer from such artefacts at the diffraction point where the intrusion boundary begins to dip.
We further notice that, although most traces at all targets match well in terms of peak and trough locations, some deviances between the migrated and modelled traces are observed. Examples include the lower reflector in Fig. 14f (target 2) and the lower reflector in Figures 15(f) and 16(f) (target 3). This is most likely attributable to slight misalignments of reflector locations in the PSDM image, resulting from the smoothing of the velocity model prior to back-propagation of the wavefield in the RTM algorithm.

D I S C U S S I O N
The presented study illustrates the potential of target-oriented point-spread function (PSF) convolution modelling as a tool for simulating prestack depth-migrated (PSDM) images. For targets characterized by low-velocity contrasts, the use of one PSF estimated at the centre of the target area may be sufficient, thereby reducing computational cost. The cost may, however, still be significant if a wave-based approach is applied for PSF estimation. An alternative option is therefore to employ ray-based PSFs provided that the ray approach itself does not break down due to caustics, multipathing, etc. While the estimation of the wave-based PSFs in the BP Statics Benchmark Case required approximately ten hours of computation time on a standard workstation (3.40 GHz Intel core), the raybased PSFs were estimated in approximately 20 seconds on the same workstation.
Our findings illustrate that images obtained through PSF convolution modelling, whether we apply a wave-based or ray-based PSF, accurately capture most illumination, resolution and amplitude effects observed in PSDM images. Amplitude effects such as migration artefacts, interference issues and wave-energy originating from outside a selected target area may, however, be challenging to capture. Future studies could assess possible approaches for capturing such effects in the PSF operators.
Regarding the PSF operators themselves, we observed that simulated PSDM images obtained via wave-based or raybased PSFs were virtually indistinguishable for the applied synthetic models. The present study thus suggests that when PSF convolution modelling is applied for target-oriented seismic modelling, especially addressing the needs in, e.g. seismic interpretation, divergence effects between a wave-based and ray-based approach may be negligible. As such, ray-based PSFs may indeed offer an efficient and flexible alternative to wave-based PSFs given the very low computational cost involved. This is also confirmed from the analysis of the derived governing equations. However, as the derived equations only consider a 2D homogeneous isotropic medium, we acknowledge that further studies are required for more thorough assessments of how valid these findings are for complex cases. Furthermore, although the obtained results from the selected target areas of the BP Statics Benchmark model used in this study were found to be similar, other complex models with different geological challenges may yield greater divergence.
For more accurate seismic modelling, further refinements may be considered. One possibility is to include effects caused by anisotropy in the PSFs (Lecomte and Kaschwich, 2008). Another possibility is to include a more precise description on how the range of incidence angles impacts reflectivity, such as observed in amplitude versus offset (AVO)/amplitude versus angle (AVA) studies. Resolution of a target area of interest is determined by the combined effects of angle-dependent reflectivity and angle-dependent illumination/resolution (Lecomte, 2008). Future studies may therefore provide more accurate assessments of how the modelled amplitudes match migrated amplitudes without resorting to normalization, and which calibration strategies are most effective. Furthermore, the plane-wave assumption inherent in the ray-based PSF may be extended to a parabolic assumption to better capture curvature effects present in migrated images and wave-based PSFs (Gelius et al., 2002b). For the near-field, higher-order paraxial ray theory could also be applied for more precise wavefront approximations. Yet another possibility is to design PSFs for elastic modelling with PS-converted waves. This extension is straightforward as it merely involves applying a scaling factor to the individual scattering wavenumber vectors in the wavenumber domain (Gelius et al., 2002a). Such an approach could be useful to assess improvements in the imaging  Fig. 16b). (f) Traces obtained at x = 1.61 km from migrated and modelled results (yellow line marked (f) in Fig. 16b). The horizontal dashed lines in (e) and (f) represent exact reflector positions of small-scale features, as the lower S-velocities would yield increased resolution, though the latter is often counterbalanced by higher attenuation effects.
Furthermore, the 2D PSF convolution approach analysed in this study is easily extendable to 3D models as well, with only a marginal increase in the computational cost. Several studies already document the potential of applying ray-based PSFs for simulating PSDM images of complex 3D geology (e.g. Lecomte, 2008;Lecomte et al., 2015;Jensen et al., 2021). Due to the reduced cost, simulation of 3D PSDM images may be obtained quickly on detailed models with significantly denser sampling than what FD approaches may allow. As such, flexible sensitivity studies of how various model parameters (geology, petrophysical properties) and seismic parameters (illumination, survey geometry etc.) affect seismic images may be performed. For this study, our limitation to simpler 2D models was motivated from our need to compare the ray-based results to results obtained via costly wave-based approaches. However, further validation studies aimed towards 3D migrated data should be a potential future area of study.
In addition to improving the validity of PSF convolution modelling as a forward modelling approach, the abovementioned suggestions for future studies can also improve the application of PSFs as deconvolution operators. Even just a partial deconvolution of a migrated image may yield a substantial saving of computation time prior to, for instance, least-squares RTM (Aoki and Schuster, 2009). The speed in which ray-based PSFs can be designed may allow for efficient and flexible adjustments of the needed parameters, and further refinements are therefore desired for improved accuracy and validity. Overall, the choice of method is especially related to the application and data/method available. If one has full access to, and control of, the migration parameters, a wavebased PSF might be a better choice due to the guaranteed accuracy. However, when exact migration parameters are not available, ray-based approaches provide an easy, cheap and flexible alternative, providing one has a minimum of information about the velocity model and the migration needs.

C O N C L U S I O N S
Target-oriented point-spread functions (PSFs) used as spatial convolution operators allow for accurate seismic simulation of prestack depth-migrated (PSDM) images. Analytical modelling in a homogeneous velocity model based on governing equations derived for wave-and ray-based PSFs reveals negligible differences between PSFs obtained via the two approaches for the homogeneous case. Further comparisons between PSF-convolved images and actual PSDM images obtained at selected targets in the BP Statics Benchmark model suggest that PSF convolution modelling accurately captures resolution, illumination and amplitude effects observed in migrated images. The results indicate that regardless of whether the PSFs are estimated through a wave-based or ray-based approach, the modelled results are usually virtually indistinguishable. Some minor divergence may, however, occur at targets characterized by close proximity between shot point and target point, high-velocity and/or low-frequency bandwidth, and at targets where the ray tracing algorithm may not fully account for amplitude and curvature effects. These effects may be primarily attributable to the far-field approximation inherent in ray theory. Further quantitative analyses of differences in PSDM images simulated via wave-based and ray-based PSFs should assess whether the former approach yields sufficient improvements to justify the increased cost.

AC K N OW L E D G E M E N T S
The first author wishes to thank the University of Bergen for PhD project funding. Furthermore, the authors gratefully acknowledge NORSAR Innovation AS for providing an academic license of their modelling software. We acknowledge Amoco Tulsa Research Lab for developing and providing the BP Statics Benchmark model. We also thank Dr. Paul Lubrano-Lavadera for assistance in developing parts of our implemented computer codes. Finally, the authors would like to thank the Research Council of Norway for financial support through project #26763 (FOPAK).

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author (first author) upon reasonable request.