Search

SUMMARY Severely crooked 2-D reflection profiles generate a swath distribution area of common midpoint (CMP) traces. These traces are scattered in all directions and no longer correspond to CMP gathers directly below the acquisition line, effectively posing a pseudo-3-D problem. Standard seismic processing techniques neglect the 3-D character of crooked seismic transects and often yield suboptimal images. The out-of-plane energy ignored in 2-D stacking techniques causes amplitude smearing effects due to anomalous reflection traveltimes. The imaging of dipping structures in geologically complex areas, such as in Larder Lake, are particularly affected because of complex energy scattering. These out-of-plane reflected wavefield components degrade energy focusing and signal alignment, worsening as the obliqueness of the profile with respect to local geology increases. To mitigate the contribution of scattered out-of-plane energy, we implement azimuthal binning. Pseudo-processing lines, parallel to the maximum geological dip, are generated by using locally optimum bin orientations. We use a priori geological information to devise optimum bin orientations in accordance to the local strike. We find that the optimized binning method significantly outperforms the standard approach in terms of in-phase amplitude stacking and reflection continuity. This is because the standard approach is based on a bin geometry always perpendicular to the seismic line traverse, ignoring 3-D effects, whereas in the optimum binning approach each individual bin may have a different orientation along the line. Moreover, we enhance the final reflectivity image by performing seismic migration with a velocity model that incorporates near-surface velocity estimations from full-waveform inversion. We demonstrate that the combination of these two innovations results in a significant impact on the imaging quality of structures at depth. In addition to imaging structures at depth, we also illustrate the significance of their continuation towards the near surface. This is especially important in the Larder Lake area where tectonic activity may have emplaced ore deposits in the first few kilometres of the subsurface.

We successfully apply the acoustic Laplace–Fourier waveform tomography method to delineate P-wave velocity structures of the mega-splay fault system in the central part of the seismogenic Nankai subduction zone offshore Japan, using densely sampled wide-angle ocean bottom seismograph (OBS) data originally acquired in 2004. Our success is due to new and carefully designed data preconditioning and inversion strategies to mitigate (i) the well-known non-linearity of waveform inversion, (ii) the challenges arising from crustal-scale survey designs (e.g. undersampling of the OBSs), and (iii) modelling errors due to the use of the acoustic assumption.We identify a sixfold set of key components that together lead to the success of the high-resolution waveform tomography image: (i) Availability of low-frequency components (starting at 2.25 Hz) reducing the non-linearity, and access to large offset data (up to 55 km) increasing the depth of illumination and the recovery of low wavenumber components. (ii) A highly accurate traveltime tomography result (with an rms error of approximately 60 ms) that further mitigates the non-linearity. (iii) A hierarchical inversion approach in which phase spectra are inverted first to reduce artefacts from the acoustic assumption, and amplitude information is only incorporated in the final stages. (iv) A Laplace–Fourier domain approach that facilitates a multiscale approach to mitigate non-linearity by restricting the inversion to the low frequency components and early arrivals first, and sequentially including higher frequencies and later arrivals. (v) A pre-conditioning strategy for eliminating undesirable high wavenumber components from the the gradient. (vi) A strategy for source estimation that reduce the influence of the instrumental design.In the OBS case study used for illustration purposes, Laplace–Fourier waveform tomography retrieves velocity anomalies as small as 700 m (horizontally) and 350 m (vertically) above the top of the Philippine Sea Plate. The resulting velocity structures include low-velocity zones and thrust structures which have not been previously identified clearly. The velocity models are validated by scrutiny of synthetic and observed waveforms, by evaluating the coherency of source estimates, and by comparison with 3-D pre-stack migrated (PreSDM) images. Chequerboard tests and point-scatter tests demonstrate both the reliability and the limitations of the acoustic implementation.

SUMMARY A comprehensive validation of 2-D, frequency-domain, acoustic wave-equation tomography was undertaken in a ‘blind test’, using third-party, realistic, elastic wave-equation data. The synthetic 2-D, wide-angle seismic data were provided prior to a recent workshop on the methods of controlled source seismology; the true model was not revealed to the authors until after the presentation of our waveform tomography results. The original model was specified on a detailed grid with variable P-wave velocity, S-wave velocity, density and viscoelastic Q-factor structure, designed to simulate a section of continental crust 250 km long and 40 km deep. Synthetic vertical and horizontal component data were available for 51 shot locations (spaced every 5 km), recorded at 2779 receivers (spaced every 90 m), evenly spread along the surface of the model. The data contained energy from 0.2 to 15 Hz. Waveform tomography, a combination of traveltime tomography and 2-D waveform inversion of the early arrivals of the seismic waveforms, was used to recover crustal P-velocity structure from the vertical component data, using data from 51 sources, 1390 receivers and frequencies between 0.8 and 7.0 Hz. The waveform tomography result contained apparent structure at wavelength-scale resolution that was not evident on the traveltime tomography result. The predicted (acoustic) waveforms in the final result matched the original elastic data to a high degree of accuracy. During the workshop, the exact model was revealed; over much of the model the waveform tomography results provided a good correspondence with the true model, from large- to intermediate-(wavelength) scales, with a resolution limit on the order of 1 km. A significant, near-surface low-velocity zone, invisible to traveltime methods, was correctly recovered; the results also provided a high-resolution image of the complex structure of the entire crust, and the depth and nature of the crust–mantle transition. Some inaccuracies were observed near the edges of the images due to the limited ray coverage inherent to the footprint of the survey geometry. Several aspects of the waveform tomography strategy were critical to the success of the acoustic method with realistic, synthetic, viscoelastic data: (i) the accuracy of the starting model from traveltime tomography, (ii) implementation in the frequency domain, (iii) the use of complex-valued frequencies to effect time damping of the data residuals, (iv) the selection of a suitable subset of data and data frequencies, (v) progressive inversion of low- to high-frequency components of the data, (vi) initial, pre-inversion matching of the amplitudes between observed and modelled data, and (vii) sufficient preconditioning of both the data and the update images. Combined, these strategies were effectively equivalent to a multiscale approach that mitigated the non-linearity of the seismic inverse problem. During the inversion we carried out repeated forward modelling to ensure our modelled waveforms matched the observed data as closely as possible in both frequency and time domains. The synthetic data set used in this paper provides a benchmark for future testing of modelling, inversion, and imaging algorithms for wide-angle lithospheric imaging.

SUMMARY We provide a series of numerical experiments designed to test waveform tomography under (i) a reduction in the number of input data frequency components (‘efficient’ waveform tomography), (ii) sparse spatial subsampling of the input data and (iii) an increase in the minimum data frequency used. These results extend the waveform tomography results of a companion paper, using the same third-party, 2-D, wide-angle, synthetic viscoelastic seismic data, computed in a crustal geology model 250 km long and 40 km deep, with heterogeneous P-velocity, S-velocity, density and Q-factor structure. Accurate velocity models were obtained using efficient waveform tomography and only four carefully selected frequency components of the input data: 0.8, 1.7, 3.6 and 7.0 Hz. This strategy avoids the spectral redundancy present in ‘full’ waveform tomography, and yields results that are comparable with those in the companion paper for an 88 per cent decrease in total computational cost. Because we use acoustic waveform tomography, the results further justify the use of the acoustic wave equation in calculating P-wave velocity models from viscoelastic data. The effect of using sparse survey geometries with efficient waveform tomography were investigated for both increased receiver spacing, and increased source spacing. Sampling theory formally requires spatial sampling at maximum interval of one half-wavelength (2.5 km at 0.8 Hz): For data with receivers every 0.9 km (conforming to this criterion), artefacts in the tomographic images were still minimal when the source spacing was as large as 7.6 km (three times the theoretical maximum). Larger source spacings led to an unacceptable degradation of the results. When increasing the starting frequency, image quality was progressively degraded. Acceptable image quality within the central portion of the model was nevertheless achieved using starting frequencies up to 3.0 Hz. At 3.0 Hz the maximum theoretical sample interval is reduced to 0.67 km due to the decreased wavelengths; the available sources were spaced every 5.0 km (more than seven times the theoretical maximum), and receivers were spaced every 0.9 km (1.3 times the theoretical maximum). Higher starting frequencies than 3.0 Hz again led to unacceptable degradation of the results.

SUMMARY Modern wide-angle surveys are often multi-fold and multi-channel, with densely sampled source and receiver spacings. Such closely spaced data are potentially amenable to multi-channel techniques involving wavefield propagation methods, such as those commonly used in reflection data processing. However, the wide-angle configuration requires techniques capable of handling very general wave types, including those not commonly used in reflection seismology. This is a situation analogous to that faced in cross-borehole seismics, where similar wave types are also recorded. In a real cross-borehole example, we compare pre-stack migration, traveltime tomography and wavefield inversion. We find that wavefield inversion produces images that are quantitative in velocity (as are the tomograms) but are of significantly higher resolution; the wavefield inversion results have a resolution comparable to that of the (qualitative) pre-stack migration images. We seek to extend this novel development to the larger-scale problem of crustal imaging. An essential element of the approach we adopt is its formulation entirely within the temporal frequency domain. This has three principal advantages: (1) we can choose to ‘decimate’ the data, by selecting only a limited number of frequency components to invert, thus making inversion of data from large numbers of source positions feasible; (2) we can mitigate the notorious non-linearity of the seismic inverse problem by progressing from low-frequency components in the data to high-frequency components; and (3) we can include in the model any arbitrary frequency dependence of inelastic attenuation factors, Q(ω), and indeed solve for the spatial distribution of Q. An initial synthetic test with an anomaly located within the middle crust yields a velocity image with the correct structural features of the anomaly and the correct magnitude of velocity anomaly. This is related to the fact that the reconstruction is obtained from forward-scattered waves. Under these conditions, the method thus behaves much like tomography. A second test with a deeper, more extensive anomaly yields an image with the correct velocity polarity and the correct location, but with a deficiency in low and high wavenumbers. In this case, this is because the reconstruction is obtained from backscattered waves; under these conditions the method behaves not like tomography, but like migration. A more extensive test, based on a large wide-angle survey in south-eastern California and western Arizona, demonstrates a real potential for high-resolution imaging of crustal structures. Although our results are limited by the acoustic approximation and by the relatively low frequencies that we can model today, the images are sufficiently encouraging to warrant future research. The problem of local minima in the objective function is the most significant practical problem with our method, but we propose that appropriate ‘layer’ stripping methods can handle this problem.

SUMMARY In principle, crosshole, traveltime tomography is ideal for directly detecting and measuring seismic anisotropy. The traveltimes of multiple rays with wide angular coverage will be sensitive both to inhomogeneities and to anisotropy. In practice, the traveltimes will depend only on a limited number of the anisotropic velocity parameters, and the data may not be adequate even to determine these parameters uniquely. In addition, trade-offs may exist between anisotropy and inhomogeneities. In this paper, we use the linear perturbation theory for traveltimes in general, weakly anisotropic media to discuss the dependence of traveltimes in 2-D crosshole tomographic experiments on the anisotropic parameters. In a companion paper, we apply the results to synthetic and real data examples. We show that when measurements are restricted to a 2-D plane, the qP and qS traveltimes depend on subsets of the complete set of 21 anisotropic velocity parameters. Formulae are developed for the differential coefficients of the traveltimes with respect to these parameters in piecewise homogeneous and in linearly interpolated models. It is shown how in a generally oriented model element, the local parameters are related to the same parameters in the global model. The parameters that can be determined from 2-D tomographic data do not in general determine the full nature of the anisotropy. Rather, these parameters serve only to describe the intersection of the slowness sheet with the 2-D plane. Since many models may fit this description, additional information on symmetry properties and orientations is required. For example, if a priori information suggests that the anisotropy is transversely isotropic (TI), then we can determine some of the TI parameters and some information on the orientation of the axis of symmetry. Formulae are given relating the general parameters to those of a TI system with general orientation of the symmetry axis. The general formulae for qS traveltimes are intrinsically more complicated than those for qP. In the qS case, the traveltime perturbation depends on the polarization, which in turn depends on the perturbation. This makes the general problem non-linear even for small perturbations. However, the mean qS traveltime and the traveltime dependence on various subsets of parameters are linear. Although linear perturbation theory is invalid for qS rays, degenerate perturbation theory is valid for the calculation of the traveltimes and could be used in a non-linear inversion scheme.

SUMMARY Cross-borehole seismic data have traditionally been analysed by inverting the arrival times for velocity structure (traveltime tomography). The presence of anisotropy requires that tomographic methods be generalized to account for anisotropy. This generalization allows geological structure to be correctly imaged and allows the anisotropy to be evaluated. In a companion paper we developed linear systems for 2-D traveltime tomography in anisotropic media. In this paper we analyse the properties of the linear system for quasi-compressional waves and invert both synthetic and real data. Solutions to the linear systems consist of estimates of the spatial distributions of five parameters, each corresponding to a linear combination of a small subset of the 21 elastic, anisotropic velocity parameters. The parameters describe the arrival times in the presence of weak anisotropy with arbitrary symmetries. However, these parameters do not, in general, describe the full nature of the anisotropy. The parameters must be further interpreted using additional information on the symmetry system. In the examples in this paper we assume transverse isotropy (TI) in order to interpret our inversions, but it should be noted that this final interpretation could be reformulated in more general terms. The singular value decomposition of the linear system for traveltime tomography in anisotropic media reveals the (expected) ill-conditioning of these systems. As in isotropic tomography, ill-conditioning arises due to the limited directional coverage that can be achieved when sources and receivers are located in vertical boreholes. In contrast to isotropic tomography, the scalelength of the parametrization controls the nature of the parameter space eigenvectors: with a coarse grid all five parameters are required to model the data; with a fine grid some of the parameters appear only in the null space. The linear systems must be regularized using external, a priori information. An important regularization is the expectation that the elastic properties vary smoothly (an ad hoc recognition of the insensitivity of the arrival times to the fine-grained properties of the medium). The expectation of smoothness is incorporated by using a regularization matrix that penalizes rough solutions using finite difference penalty terms. The roughness penalty sufficiently constrains the solutions to allow the smooth eigenvectors in the null space of the unconstrained problem to contribute to the solutions. Hence, the spatial distribution of all five parameters is recovered. The level of regularization required is difficult to estimate; we advocate the analysis of a suite of solutions. Plots of the solution roughness against the data residuals can be used to find ‘knee points’, but for the fine tuning of the regularization one has little recourse but to examine a suite of images and use geological plausibility as an additional criterion. The application of the regularized numerical scheme to the synthetic data reveals that the roughness penalty should include terms that penalize high gradients addition to penalizing high second derivatives. Only when this constraint was included were the features of the original model recovered. The inversions of the field data yield good images of the expected stratigraphy and confirm previous estimates of the magnitude of the anisotropy and the orientation of the symmetry axis. The solutions further indicate an increase in anisotropy from the top to the bottom of the survey region that was not previously detected.