The fidelity of depth seismic imaging depends on the accuracy of the velocity models used for wave-7 field reconstruction. Models can be decomposed in two components corresponding to large scale 8 and small scale variations. In practice, the large scale velocity model component can be estimated 9 with high accuracy using repeated migration/tomography cycles, but the small scale component 10 cannot. Therefore, wavefield reconstruction does not completely describe the recorded data and 11

more »

... images are perturbed by artifacts. 12 There are two possible ways to address this problem: improve wavefield reconstruction by 13 estimating more accurate velocity models and image using conventional techniques (e.g. wavefield 14 cross-correlation), or reconstruct wavefields with conventional methods using the known smooth 15 velocity model, and improve the imaging condition to alleviate the artifacts caused by the imprecise 16 reconstruction, as suggested in this paper. 17 In this paper, the unknown component of the velocity model is described as a random function 18 with local spatial correlations. Imaging data perturbed by such random variations is characterized 19 by statistical instability, i.e. various wavefield components image at wrong locations that depend on 20 1 the actual realization of the random model. Statistical stability can be achieved by local wavefield 21 averaging either in spatial windows defined in the vicinity of the data acquisition locations, or in 22 local windows defined in the vicinity of image points. We use the latter approach and show that 23 the technique is effective in attenuating imaging artifacts without being hampered by some of the 24 limitations of data-space alternatives. 25 Seismic imaging in complex media requires accurate knowledge of the medium velocity. Assuming 26 single scattering, imaging requires propagation of the recorded wavefields from the acquisition sur-27 face, followed by the application of an imaging condition highlighting locations where scattering 28 occurs, i.e. where reflectors are present. 29 The main requirement for good-quality imaging is accurate knowledge of the velocity model. 30 Errors in the model used for imaging lead to inaccurate reconstruction of the seismic wavefields and 31 to distortions of the migrated images. In a realistic seismic experiment the velocity model is not 32 known exactly. Migration velocity analysis produces large scale approximations of the model, but 33 fine scale variations remain elusive. For example, when geology includes complicated stratigraphic 34 structures, the rapid velocity variations on the scale of the seismic wavelength and smaller cannot 35 be estimated correctly. Therefore, even if the broad kinematics of the seismic wavefields are recon-36 structed correctly, the extrapolated wavefields also contain distortions that lead to image artifacts 37 obstructing the image of the geologic structure under consideration. While it is certainly true that 38 even the recovery of a long-wave background may prove to be a challenge in some circumstances, 39 we do not attempt to address that issue in this paper. Instead, we concentrate solely on the problem 40 of dealing with the effect of a small scale random variations. 41 There are two ways in which we can approach this problem: 42 • The first option is to improve the velocity analysis methods to estimate the small-scale vari-43 ations in the model. Such techniques take advantage of all information contained in seismic 44 wavefields and are not limited to kinematic information of selected events picked from the 45 data. Examples of techniques in this category are waveform inversion (Tarantola, 1987; Pratt 46 and Worthington, 1990; Pratt, 1990), wave-equation tomography (Woodward, 1992) or wave-47 3 equation migration velocity analysis (Sava and Biondi, 2004a,b; Shen et al., 2005). A more 48 accurate velocity model allows for more accurate wavefield reconstruction. Then, wavefields 49 can be used for imaging using conventional procedures, e.g. cross-correlation. 50 • The second option is to concentrate on the imaging condition, rather than concentrate on 51 wavefield reconstruction. Assuming that the large-scale component of the velocity models 52 is known (e.g. by iterative migration/tomography cycles), we can design imaging conditions 53 that are not sensitive to small inaccuracies of the reconstructed wavefields. Imaging artifacts 54 can be reduced at the imaging condition step, despite the fact that the wavefields incorporate 55 small kinematic errors due to velocity fluctuations. 56 The two options are complementary to each other, and both can contribute to imaging accuracy. In 57 this paper, we concentrate on the second approach. For purposes of theoretical analysis, it is con-58 venient to model the small-scale variations velocity fluctuations as random but correlated variations 59 superimposed on a smooth known velocity. We assume that we know the smooth background, but 60 we do not know the random fluctuations. The goal is to design an imaging condition that alleviates 61 artifacts caused by those random fluctuations. 62 Conventional imaging consists of cross-correlations of extrapolated source and receiver wave-63 fields at image locations. Since wavefield extrapolation is performed using an approximation of the 64 true model, the wavefields contain random time delays, or equivalently random phases, which lead 65 to imaging artifacts. 66 One way of mitigating the effects of the random model on the quality of the resulting image 67 involves using techniques based on acoustic time reversal (Fink, 1999). A well-known result is 68 that, under certain assumptions, a signal sent through a random medium, recorded by a receiver 69 array, time reversed and sent back through the same medium, refocuses at the source location in 70 4 a statistically stable fashion. Statistical stability means that the refocusing properties (i.e. image 71 quality) are independent of the actual realization of the random medium (Papanicolaou et al., 2004; 72 Fouque et al., 2005). 73 An equivalent formulation of this result is that statistical stability can be achieved if random 74 phase shifts between signals recorded at nearby locations are removed by cross-correlation, prior 75 to back propagation into the medium. This observation lies at the heart of coherent interferometric 76 imaging (Borcea et al., 2005, 2006c,a,b) or imaging and velocity analysis in presence of uncer-77 tain models (Dussaud, 2005). With this approach, local spatial cross-correlations of data traces on 78 the acquisition surface are computed and extrapolated to the image location using standard tech-79 niques. Thus, random phase shifts that cause imaging artifacts are removed prior to extrapolation 80 and imaging quality in the smooth background medium increases. This method assumes the ex-81 istence of a continuous acquisition array with dense receiver layout to ensure coherency of local 82 cross-correlations. 83 We investigate an alternative way of using time reversal to increase imaging statistical stability. 84 Instead of coherent interferometry applied to data on the acquisition surface, we first extrapolate 85 wavefields to all locations in the imaging volume and then apply local spatial cross-correlations in 86 the vicinity of every image point. Correlations in the image-space damps small random fluctuations 87 in the extrapolated wavefields. The cross-correlations do not relocate energy in space, but simply 88 produce local averages of the extrapolated wavefields. 89 The procedure closely resemble conventional imaging procedures where wavefields are extrap-90 olated in the image volume and then cross-correlated in time at every image location. Our method 91 uses averaging in three-dimensional local windows around image locations. From implementa-92 tion and computational cost points of view, our technique does not differ much from conventional 93 5 imaging, although the imaging properties are improved. Therefore, we use the name interferomet-94 ric imaging condition for our technique to contrast it with conventional imaging condition which 95 represents a special case of this method for infinitely small local windows. 96 One advantage of using the interferometric imaging condition in wave-equation migration is 97 that it also makes efficient use of the data obtained by a sparse array, since the cross-correlation is 98 performed at an image point on wavefields extrapolated from all data traces simultaneously. Fur-99 thermore, local averaging around the image locations is inherently three-dimensional, in contrast 100 with the two-dimensional averaging typical for interferometric imaging parametrized on the sur-101 face. This increases signal-to-noise ratio and improves random phase cancellation, although it also 102 increases computational cost proportionally. Processing with image-space coordinates is simpler 103 than processing using data-space coordinates because, after extrapolation, wavefields are simpler 104 since wave propagation complications have been unraveled in the extrapolation process (Stolk and 105 Symes, 2004). 106 Finally, we note that the proposed method is not an extension, but an alternative to coherent in-107 terferometry. When assumptions permit, the two methods can be used simultaneously in an imaging 108 functional, thus taking advantage of the best qualities of both approaches. 109 SEISMIC IMAGING CONDITIONS Conventional seismic imaging methods share the fundamental assumption of single scattering in 110 the subsurface (Born approximation). Under this assumption, we can represent waves recorded at 111 the surface as a convolution of Green's functions (G) corresponding to sources on the surface and 112 scattering points in the subsurface. Assuming an impulsive source, we can write 113 P (x m , ω m ) = G (x s , y m , ω m ) G (y m , x m , ω m ) , where P denotes recorded acoustic data at coordinates x m , x s are coordinates of the source, y m are 114 coordinates of the scattering points, and ω m is the frequency of the propagating wave. Relation (1) 115 can be written in an equivalent form using time instead of frequency variables, but for simplicity 116 we use the frequency-domain notation throughout this paper. The Green's function G characterizes 117 data propagation in the real medium of velocity v. 118 Imaging with recorded data P (x m , ω m ) is a two-step procedure: 119 • The first step consists of extrapolation of source and receiver wavefields from the recording 120 surface to image locations. The source wavefield corresponds to simulated waves propagating 121 forward in time from the source location x s , and the receiver wavefield corresponds to waves 122 propagating backward in time from recording locations x m . 123 • The second step consists of an imaging condition evaluating whether the two extrapolated 124 wavefields match kinematically, which indicates whether a reflector is present in the medium. 125 Wavefield extrapolation and imaging can be implemented in different space and time domains, 126 for example downward continuation with the one-way wave-equation implemented in frequency-127 wavenumber, frequency-space or mixed domains, or wavefield extrapolation with the two-way 128 wave-equation in time-space, etc. Wavefield extrapolation can also be performed using Kirchhoff 129 integral methods followed by conventional imaging. However, the actual extrapolation method is 130 irrelevant for the discussion in this paper. For simplicity, we represent wavefield extrapolation with 131 Green's functions in the frequency-domain, although our examples use time-domain two-way finite-132 difference solutions to the acoustic wave-equation, process typically known as reverse-time migra-133 tion (Kosloff and Baysal, 1983). Any other implementation of Green's functions can be substituted 134 in our discussion with no change to the conclusions. 135 157 summation over temporal frequency, indicates the presence of reflectors. Here and for the rest of 158 the paper, summation over multiple seismic experiments (sources) is assumed. 159 The propagation geometry of the imaging condition (2) is summarized in figure 1(a): source 160 waves are propagated forward in time from the source at x s to the scatterer location y m , and receiver 161 waves are propagated backward in time from the receiver coordinates x m to the scatterer location 162 y m . 163 352 e.g. downward continuation, Kirchhoff integral methods, etc. 353 The parameters used in our examples, explained in Appendix A, are: 354 • seismic spatial wavelength λ = 20 m, 355 • wavelet central frequency ω = 150 Hz, 356 • background velocity: v 0 = 3000 + k z m/s (k = 1s −1 ), 357 • random fluctuations parameters: r a = 30, r c = 5, α = 2, and 358 • random velocity v constructed from the background velocity v 0 . 359 Consider the model depicted in figures 10(a)-10(f). As in the preceding example, the left panels 360 depict the known smooth velocity v 0 , and the right panels depict the model with random variations. 361 The imaging target is represented by the oblique lines, figure 9 (b), located around z = 700 m, which 362 simulate a cross-section of a complex stratigraphic model (e.g. meandering channels). 363 We model data with random velocity and image using the smooth model. Figures 10(a)-10(f) 364 show wavefield snapshots in the two models for different propagation times, one before the source 365 wavefields interact with the target reflectors and one after this interaction. Figures 11(a) and 11(b) 366 show the corresponding recorded data on the acquisition surface located at z = λ, where λ represents 367 the wavelength of the source pulse. 368 Migration with conventional imaging condition of the data simulated in the smooth model us-369 ing the smooth velocity produces the images in figures 12(a)-12(c). The targets are well imaged, 370 although the image also shows artifacts due to truncation of the data on the acquisition surface. In 371 contrast, migration with the conventional imaging condition of the data simulated in the random 372 19 model using the smooth velocity produces the images in figures 12(b)-12(d). Those images are 373 distorted by the random variations in the model that are not accounted for in the smooth migration 374 velocity. The targets are hard to discern since they overlap with many truncation and defocusing 375 artifacts caused by the inaccurate migration velocity. Finally, figures 13(a)-13(c) show the migrated 376 images for the same situation as the one depicted in figures 12(a)-12(c), except that migration uses 377 the interferometric imaging condition from equation (5). In this situation, since we are using the 378 same model for modeling and migration, the interferometric imaging condition is not expected to 379 change the image much. Interestingly, some of the truncation artifacts are attenuated, but otherwise 380 the images are similar and the targets are easy to identify. Similarly, figures 13(b)-13(d) show the 381 migrated images for the same situation as the one depicted in figures 12(b)-12(d), with migration 382 using the interferometric imaging condition in equation (5). Many of the artifacts caused by the 383 inaccurate velocity model are suppressed and the imaging targets are more clearly visible and easier 384 to interpret. Furthermore, the general patterns of amplitude variation along the imaged reflectors 385 are similar between figures 12(d) and 13(d). We note that the reflectors are not as well imaged as 386 the ones obtained when the velocity is perfectly known. This is because the interferometric imaging 387 condition described in this paper does not correct kinematic errors due to inaccurate velocity. It only 388 acts on the extrapolated wavefields to reduce wavefield incoherency and add statistical stability to 389 the imaging process. Further extensions to the interferometric imaging condition can improve fo-390 cusing and enhance the images by correcting wavefields prior to imaging. However, this topic falls 391 outside the scope of this paper and we do not elaborate on it further. 392 DISCUSSION The interferometric imaging conditions defined by equations (5) or (7) represent extensions of con-393 ventional cross-correlation imaging conditions. Statistical stability in presence of model random 394 20 variations is achieved by averaging in local windows in space and time. The computational cost of 395 interferometric imaging condition is proportional to the number of samples in the image multiplied 396 by the number of cross-correlation lags involved in the local space and time averaging, i.e. the 397 size of the image times the sizes of the decoherence length and time. When the size of the deco-398 herence windows is reduced to zero, the interferometric imaging condition becomes similar to the 399 conventional autocorrelation imaging condition with identical properties and implementation cost. 400 As indicated earlier, the actual method used for wavefield reconstruction is not relevant for our 401 discussion about this imaging condition. Any type of wavefield extrapolation can be substituted 402 when reconstructing Green's functions in the velocity model without random variation. Thus, this 403 form of imaging condition can be applied equally well to migration by downward continuation using 404 space, wavenumber or mixed space-wavenumber extrapolation, reverse-time migration or Kirchhoff 405 migration. 406 Areas of application for interferometric imaging include geologic areas with rapid velocity 407 variations caused by stratigraphic features unaccounted for in typical velocity analysis or imaging 408 through highly scattering media, e.g. basalt, anhydrite or inhomogeneous salt. 409 CONCLUSIONS Conventional seismic imaging conditions based on wavefield cross-correlations are extended to 410 achieve statistical stability for models with rapid, small-scale velocity variations. Random velocity 411 variations on a scale comparable with the seismic wavelength are modeled by correlated Gaussian 412 distributions. 413 Statistical stability is achieved by local averaging of cross-correlated wavefields at image lo-414 cations. The proposed interferometric imaging condition is a natural extension of and reduces to 415 21 the conventional cross-correlation imaging condition when the averaging window is made infinitely 416 small. 417 The main characteristic of the method is that it operates on extrapolated wavefields at image 418 positions (thus the name interferometric imaging condition), in contrast with alternative approaches 419 involving migration of interferograms obtained by local data cross-correlations on the acquisition 420 surface. 421