Angle-domain common-image gathers from Fresnel volume migration

. In complex geological settings, such as in hard-rock environments, Fresnel volume migration (FVM) has been suc-cessfully applied and found to deliver superior image quality compared to conventional imaging techniques. However, previous studies on FVM have mainly focused on obtaining kinematic seismic images, and the analysis of the migrated amplitudes has not received major attention. Therefore this study presents a method for constructing angle-domain common-image gathers (ADCIGs) and common-angle stacks from FVM, which can facilitate prestack amplitude analysis from the migrated seis- 5 mic data in the angle- domain. These ADCIGs were constructed inside the migration loops using phase slowness vectors derived from traveltime gradient ﬁelds. We then tested this method on synthetic and ﬁeld seismic data and investigated the reliability of the output for amplitude versus angle (AVA) analysis. The test results obtained showed that the AVA responses from the common-angle stacks resemble that of the input synthetic shot gather of migration relatively well, indicating the promising feasibility of AVA analysis from common-angle stacks. When implemented on ﬁeld data acquired from a hard-rock 10 environment, the proposed method can provide common-angle stacks with a higher signal-to-noise ratio and better reﬂection coherency compared to the common-angle stacks from the standard Kirchhoff prestack depth migration. This study extends the implementation of FVM toward amplitude analysis, which can help improve the feasibility of hard-rock characterization. of test the in Hence, results this study help


Fresnel Volume Migration
Using the exact Kirchhoff weighting function, FVM obtains the image value U of a subsurface image point at the lateral distance x and depth z by weighted summation of the recorded wavefield over the corresponding diffraction surface A for this 60 image point (Buske et al., 2009): where: 2.2 Angle-domain common-image gather and common-angle stack 80 In prestack depth migration, a subsurface image point can be understood as the coincidence of the incident wavefront coming from the source with the emergent wavefront reflected at the image point and propagating toward the receiver (Audebert et al., 2002). These wavefronts can be parameterized by the source and receiver isochrons, wherein each isochron has its associated attributes, such as a traveltime gradient and a phase slowness vector. The resultant of the source and receiver phase slowness vectors produces an illumination slowness vector. Furthermore, four types of diffraction angles can be defined with the source 85 and receiver phase slowness vectors (Audebert et al., 2002;Liu et al., 2018), as shown in Fig. 2. The four types of diffraction angles are: 1. scattering opening angle (α), that is, the angle between the source (P s ) and receiver (P r ) phase slowness vectors; 2. scattering azimuth (φ), that is, the azimuth of the plane containing both the source and receiver phase slowness vectors; 3. illumination dip angle (Θ), that is, the angle between the normal and the illumination slowness vector (P m ); and 90 4. illumination azimuth (Φ), that is, the azimuth of the illumination slowness vector. Fig. 2 shows that at the subsurface image point with a pair of source and receiver isochrons, the scattering opening angle α can be calculated from the phase slowness vectors P s and P r (Liu et al., 2018), and the incident angle θ can be calculated as half the scattering opening angle α as follows: . Geometric relationship between the phase slowness vectors and the diffraction angles (modified from Liu et al., 2018). In this figure, S and R represent the source and receiver on the surface, respectively; Is and Ir represent the source and receiver isochrons, respectively; and N, E, W, and Z denote the spatial axes. The black dot at the center of the spatial axes marks the subsurface image point. In addition, Ps, Pr, and Pm represent the source phase, receiver phase, and illumination slowness vectors, respectively; whereas α, φ, Θ, and Φ denote the scattering opening angle, scattering azimuth, illumination dip angle, and illumination azimuth, respectively.
Here, P s and P r can be determined from the source (T s ) and receiver (T r ) traveltime gradients: An ADCIG at the subsurface image point can be constructed by binning the prestack images into the corresponding incidentangle bins θ. Therefore, Eq. 1 can be reformulated as follows: Note that the implementation of Eq. 10 requires the binning interval and size to be properly chosen. The smaller the binning interval and size are, the more precise the final prestack images are in the incident-angle bins in representing the reflection amplitude from a particular incident angle and the more memory resource (or longer runtime) needed for the computation.

105
However, a very small binning interval and size will increase the risk of failing to capture the prestack images into the corresponding incident-angle bins. Hence, a binning interval and size of 1 • and 0.5 • , respectively, are considered reasonable choices for most cases. These parameters will result in incident angles of 0 • ≤ θ ≤ 0.5 • in the first bin, 0.5 • < θ ≤ 1.5 • in the second bin, 1.5 • < θ ≤ 2.5 • in the third bin, and so forth. Depending on the binning size, a correction factor needs to be multiplied to the prestack images in the first and last incident-angle bins to compensate for the different sizes of these bins. For example, for 110 a binning size of 0.5 • , the prestack images in the first and last incident-angle bins are multiplied by a factor of 2.
In practice, the final images from KPSDM and FVM are obtained by stacking the migrated amplitudes of the same image points from all the migrated shot gathers. With the constructed ADCIGs, the stacking can now be performed for different incident angles at every image point. The stacking can be viewed as filling the available incident-angle bins at every image point with the corresponding migrated amplitudes from all the possible source-receiver offsets. The stacking results in a common-115 angle stack at each spatial location of the migrated seismic volume, which consists of optimal amplitude illumination and all the possible prestack images from all the possible incident angles at that location.

Amplitudes Versus Angle analysis
AVA analysis aims to examine the partitioning of incident wave energy at an interface formed by two media as a function of the incident angle and the contrasts of the P-wave and S-wave velocities and the density between the two media. The partitioning of 120 the wave energy in reflection seismic data can be described by angle-dependent reflection coefficients and formulated using the Zoeppritz equations (Zoeppritz, 1919). In this study, we performed an AVA analysis using the two-term Aki-Richard equations (Richards and Frasier, 1976;Aki and Richards, 1980) as an approximation to the exact solution of the Zoeppritz equations: where R denotes the reflection coefficient as a function of the incident angle θ. V P , V S , and ρ in the equations denote the averages of the P-wave velocity, S-wave velocity, and density of the two media, respectively. The operator ∆ denotes the difference of the corresponding parameters between the media above and below the interface. The A and B are the AVA intercept and gradient, respectively. While the AVA intercept corresponds to the normal incidence reflectivity, the AVA gradient 3 Synthetic data test

135
In general, the accuracy of estimating reflection coefficients from field data is restricted by the inaccurate knowledge of the source impulse and downward continuation process, resulting from unknown absorption effects and the true medium velocity (Temme, 1984). We can avoid these problems by adopting synthetic seismic data in which the medium velocity is known, and the source impulse can be accurately estimated.

Synthetic seismic data 140
We generated a synthetic shot gather for a model consisting of two isotropic and homogeneous layers, representing a shale overlaying a gas-sand lithology (Fig. 3). The reflector lies in the model at a depth of 2000 m, and the synthetic shot gather was generated using a finite-difference modeling code (Hellwig, 2017) with a zero-phase wavelet as the source signal. We performed the modeling for an off-end cable-spread acquisition geometry with 301 receivers, including a zero-offset receiver, a 5 s trace length, and 10 m source and receiver intervals. We then sorted these 301 synthetic traces to obtain a whole 2D 145 seismic data set containing 301 synthetic shot gathers, in which each shot has a 3000 m far offset (Fig, 5a).

Zero-offset reflection coefficients
We began by analyzing the reliability of FVM and KPSDM using the exact Kirchhoff weighting function for calculating the zero-offset reflection coefficient. The zero-offset reflection coefficient R at an image point (x 0 , z 0 ) from the FVM output can be obtained by dividing the amplitude of the migrated wavefield U by that of the zero-offset incident wavefield A at the image 150 point (Temme, 1984): We obtained the migrated wavefield U (x 0 , z 0 , t s (x 0 , z 0 )) by implementing FVM (Eq. 10) to the input synthetic shot gathers.
The migration was conducted down to a depth of 4 km, with lateral migration grid intervals of 20 and 200 m, parallel and perpendicular to the 2D synthetic seismic line, respectively, and a vertical grid interval of 20 m. The true model velocity 155 (Fig. 3) was used to compute the migration traveltimes. Figure 4a shows the output migrated shot gather at the middle of the seismic line.
We calculated the amplitude of the zero-offset incident wavefield A(x 0 , z 0 , θ = 0 • ) at the image point by estimating beforehand the wavefield amplitude A 0 (θ = 0 • ) at the source normal to the image point from the input synthetic shot gather data. For a 2D isotropic and homogeneous medium, the amplitude decay as the wavefield travels from the source to a point at a distance 160 of x holds the following relationship (Aki and Richards, 1980): where A(x) is the amplitude of the direct wave at a distance of x from the source and A 0 is the wavefield amplitude at the source. We estimated A 0 by implementing Eq. 15 and utilizing the direct-wave amplitudes at the receivers. In particular, we used the receiver at a relatively large distance from the source, i.e., at 1000 m, to obtain the direct-wave amplitude in the 165 far-field of the source where the Eq. 15 holds. Finally, we calculated the amplitude of the zero-offset incident wavefield at the image point as follows: where A 0 (θ = 0 • ) denotes the wavefield amplitude at the source normal to the image point, and r is the distance from the source to the image point. Figure 4b shows the calculated reflection coefficients from the migrated shot gather at the middle 170 of the seismic line shown in Fig. 4a. Figure 4b also shows that the zero-offset reflection coefficient of the image point at the reflector below the source is −0.068 (rounded to three decimal places). For comparison, we also obtained the zero-offset reflection coefficient from the input synthetic shot gather (Fig. 5a) and the theoretical zero-offset reflection coefficient directly from the geologic model in Fig. 3. Let u x be the reflected amplitude recorded at the receiver with an offset of x meters from the source and A 0 be the wavefield amplitude at the source. The 175 reflection coefficient on the reflector at the midpoint between the source and the receiver R x is: where 2r represents the distance from the source to the image point and then to the receiver. The reflection coefficients calculated from the input synthetic shot gather are shown in Fig. 5b, where the reflection coefficient at the zero offset is −0.068 (rounded to three decimal places). The small spike observed on the reflection coefficient curve in Fig. 5b is due to a slightly 180 larger difference between the reflected amplitude at the zero and 10 m (the next receiver) offsets than the average differences between the reflected amplitudes at the other consecutive offsets. In particular, the difference between the reflected amplitudes at the zero and 10 m offsets is around 1.6e−09. On the other hand, the consecutive differences of the reflected amplitudes at the further offsets (the differences between the 10 and 20 m offsets, the 20 and 30 m offsets, etc.) are around 1.8e−12, 1.6e−12, and so forth. This peculiarity is mainly an artifact resulting from the forward modeling of the synthetic shot gather and is 185 unsubstantial. In general, we can determine the theoretical zero-offset reflection coefficient using the impedance contrast between the two layers, wherein the impedance is the product of the velocity and density of the medium. For a plane wave reflected at normal incidence from a boundary surface of medium 1 lying on top of medium 2, we can calculate the reflection coefficient R as follows (Chopra and Castagna, 2014): where V and ρ are the velocity and density, respectively, with indices corresponding to the two media. With Eq. 18, the theoretical zero-offset reflection coefficient was found to be −0.071 (rounded to three decimal places).
Overall, the results show that the zero-offset reflection coefficient calculated from the FVM (−0.068) is consistent with that calculated from the input synthetic shot gather and different from the theoretical one by only around 0.003, which is 195 reasonably negligible. This difference is mainly due to the numerical imprecision in the computation, especially given the fact that we used a Real precision data format instead of Double Precision during our calculations. The results suggest that FVM and KPSDM, along with the exact Kirchhoff weighting function, can handle the amplitudes of the input data very well. They can also provide an accurate zero-offset reflection coefficient at an image point, provided that the amplitude of the incident wavefield at the image point is accurately estimated, as in the case of synthetic seismic data. To gain a better understanding of the ADCIG profiles, we consider the ADCIGs in Figs. 6c and 6f as an example and discuss their construction through amplitude summation in the Kirchhoff integral of migrations. Figure 7a illustrates the amplitude summations within KPSDM and FVM. The amplitude summation starts from the first receiver at a distance of 300 m and 215 ends at the last receiver at a distance of 3300 m along the Y-axis (see the basemap in Fig. 7a). These first and last receivers correspond to the maximum incident-angle bin of 19 • at both sides from the source position in the middle of the seismic line.
Notably, the amplitude summation of KPSDM starts with an extremely small magnitude, which is very close to zero at the first receiver, and the summed amplitudes then start to grow toward the last receiver. Therefore, KPSDM produces amplitudes in the entire incident-angle bins up to 19 • , as shown in Fig. 6c. On the other hand, the amplitude summation of FVM accumulates 220 significant amplitudes only up to incident-angle bins of 7 • , as can be observed in Fig. 6f. Furthermore, the curvature observed from incident-angle bins of 0 • up to nearly 4.5 • in Fig. 7a corresponds to the part of the TWT isochron that coincides with the diffraction surface, in which the migrated amplitudes are supposed to be weighted and summed up to produce an image.
Beyond this angle range, the curves are relatively flat, indicating very little or no contribution from the amplitudes beyond the incident-angle bin of nearly 4.5 • in the summation.

225
The following part focuses on AVA responses from ADCIGs in Figs. 6c and 6f. The AVA responses of the ADCIGs are shown in Fig. 7b and plotted together with the theoretical AVA response and AVA response from the synthetic shot gather for comparison. The intercept of each AVA response is normalized to the AVA intercept of the theoretical one. Differences between theoretical (exact) and seismic-based AVA responses for the same geologic model are common in practice, and Fig. 7b shows that our case is no exception. The difference observed is mainly due to the deviation of the underlying conditions for the 230 AVA response of the synthetic shot gather from the underlying conditions in the Zoeppritz equations, which are used for the theoretical AVA response. For example, the Zoeppritz equations imply that the amplitudes are equivalent to the reflection coefficients in the absence of modeling-and processing-related effects, such as transmission losses, attenuation, and divergence.
On the other hand, seismic amplitudes are not directly proportional to reflection coefficients because of the effects mentioned earlier. Certain distinguishing features of the Zoeppritz equations are further discussed at length by Allen and Peddy (1993), 235 Castagna (1993) and Chopra and Castagna (2014). Nevertheless, for qualitative AVA analysis, Fig. 7b shows  show that with the common-angle stack becoming closer to the middle of the seismic line, more amplitudes contribute to higher incident-angle bins, and generally fewer migration artifacts are found in the profiles. In Fig. 8a, KPSDM produces only a small number of significant amplitudes for incident-angle bins of not more than 3 • at the end of the seismic line. Figure  TWT isochrons. These artifacts are suppressed in the profile shown in Fig. 8b, mainly because of the destructive amplitude interference introduced by the stacking. Figure 8b also shows a higher contribution of amplitudes for incident-angle bins up to approximately 17 • . Beyond this angle, similar migration artifacts as in Fig. 8a  where the illumination of the amplitudes within the incident-angle bins in the common-angle stacks is at its maximum, the AVA responses from the common-angle stacks and synthetic shot gather are relatively well correlated.

Application to field data
We tested the proposed method on 2D seismic data acquired from an area in northern Finland that is dominated by ultramafic rocks and known to host mineral resources. The seismic data acquisition was carried out under the project Experiment of 280 Sodankylä Deep Exploration (XSODEX). Figure 9 shows the geometry of the seismic line, which has a total length of 22.26 km. The data acquisition was performed with a mixture of off-end and split-spread (roll-along) schemes, in which near and far offsets were set to 2 and 3300 m, respectively. A Vibroseis source was used, and the seismic signals were recorded with vertical-component geophones. The source and receiver intervals were set to 20 and 10 m, respectively. The data set consists of 452 shots with a maximum of 312 channels per shot, and 2001 samples per trace with a 2 ms sampling interval. Furthermore, we 285 processed the seismic data with the primary aim of suppressing the noise and retaining the energy lost as a result of geometrical spreading. Table 1 summarizes the key processing components and parameters. We performed KPSDM and FVM on the processed seismic data in a 3D fashion within a rectangular area of 8 × 22.5 km ( Fig. 9) and down to a depth of 8 km. Through the migrations, we constructed ADCIGs up to an incident angle of 35 • , and with a binning interval and binning radius of 1 • and 0.5 • from the corresponding bin center, respectively. A constant velocity 290 of 5600 m/s was used for the migration. Figure 10 shows the output migrated stack profiles. As a result of the focusing effect of the FVM, Fig. 10 shows that the migrated stack profile from FVM reveals noticeably more focused and coherent reflections a depth ranging from 10 to 18 km compared to KPSDM. On top of that, the FVM also reveals two distinct groups of coherent reflections at a depth ranging from 0 to 3 km between the distance ranges from 10 to12 km and from 13 to 18 km, which cannot be clearly distinguished in the profile from the KPSDM.

295
The red-dashed lines in Fig. 10 mark the position of the corresponding common-angle stacks shown in Fig. 11. Figure 11 also shows that the common-angle stack from FVM provides considerably less background noise and more coherent reflections than that from KPSDM. Furthermore, because of the lack of near-angle reflections, the nature of the field data does allow a proper AVA analysis. Nevertheless, the proposed method combined with the strength of FVM in image focusing has allowed us to obtain the common-angle stack from the area, with considerably less background noise than the common-angle stack from  that FVM would provide better common-angle stacks with a higher signal-to-noise ratio and reflection coherency for the AVA analysis as compared to those from KPSDM. Further investigations with more suitable field data are necessary to confirm this premise. geological settings, such as in hard-rock environments. 315 We also observed that the reliability of the proposed method for AVA analysis is strongly affected by the amplitude summation in the Kirchhoff integral, the incident-angle binning, and the amplitude summation in the common-angle stacking. These aspects suggest that the proposed method neither preserves the true amplitudes nor provides the true reflection coefficients.
Instead, it preserves the relative amplitudes in which an unknown multiplier may be present, but all the events are relatively correct (Hamspon, 2005). This finding agrees with the results of previous studies, which suggested that Kirchhoff migration 320 properly treats the relative amplitudes of the input data (e.g., Resnick et al., 1987). Furthermore, in contrast to other studies, which claim to produce true-amplitude prestacks (e.g., Mosher et al., 1996;Baina et al., 2002), the proposed method provides common-angle stacks that are potentially compatible with qualitative AVA analysis only. However, in practice, qualitative AVA analysis is a reliable tool for detecting AVA anomalies (see, e.g., Ostrander, 1984;Smith and Gidlow, 1987;Rutherford and Williams, 1989) and has been used in the industry as a valuable tool for directly detecting hydrocarbons (Chopra and Castagna,325 2014). On the other hand, quantitative AVA analysis relies primarily on deterministic corrections for true amplitudes, which can be rarely achieved in practice, and therefore, pose a great risk of seismic data misinterpretation (Chopra and Castagna, 2014). The misinterpretation risk is even more considerable when dealing with field data acquired from complex geological settings, such as hard-rock environments, which usually yield poor seismic data quality because of severe wavefield scattering.
Further studies are necessary to confirm the findings of this study. We suggest that further studies use synthetic data to 330 probe the implementation of this method on 3D geometries and with extra-wide offsets. The latter will be especially useful in examining the feasibility of AVA analysis directly on ADCIGs instead of common-angle stacks. When field data is used, we recommend the inclusion of fully preserved amplitude processing prior to the migration, and we also recommend comparing the AVA response of the output with that from the same data after FVM. Moreover, the field data used should be suitable for standard AVA analysis, such as field data obtained from sedimentary environments. Once the implementation of the method on 335 such field data is understood, we can extend the investigation to more challenging field data, such as field data acquired from hard-rock environments.