RTM-based Teleseismic Reflection Tomography with Free Surface Multiples

Receiver function analysis of teleseismic converted and free surface reflected phases has long been a cornerstone of lithospheric studies. Discontinuities in elastic properties are revealed by deconvolving the incident wavefield from scattered phases and projecting the time differences to depth to form an image. The accuracy of the image is determined to a large extent by the accuracy of the method and background velocity model used, but popular approaches for projecting receiver functions to depth commonly rely on simplifying assumptions of a 1D velocity and planar discontinuities. In tectonically complex regions like subduction zones and rift systems, strong heterogeneity can create an ambiguous tradeoff between the background velocity and the depth of the discontinuities. Furthermore, such structures are apt to create caustics at high frequencies, rendering ray-based methods inadequate. In order to better constrain the background velocity and correctly place the discontinuities at depth, we employ a novel reverse-time migration (RTM) based reflection tomography method. We adapt our reflection tomography from exploration seismology for use with teleseismic phases. Active source methods for exploration have focused on the annihilation of extended images - image gathers formed with different subsurface angle or offset information - as a means of judging the accuracy of the model. Applying these approaches to teleseismic data is untenable because 1) the sparse and uneven distribution of earthquake sources leads to the incomplete construction of extended image, 2) the imperfect separation and source deconvolution of the scattered wavefield render previous error measurements unreliable, and 3) the planar geometry of incoming arrivals makes measures of subsurface offset insensitive to perturbations in the model. To overcome these obstacles, we have developed a flexible approach based on pairwise single-source image correlations. We determine the success of the RTM and thus the accuracy of the background velocity model by cross-correlation of the images produced using different teleseismic sources. Single-source images are created by propagating the incident and scattered wavefields to depth using a Helmholtz operator and combining the by applying an inverse scattering operator. The error function is then comprised of the weighted image correlation power at depth windows. The optimized velocity model is the one that minimizes power in the correlations away from zero depth shift. We develop our inversion scheme using the Augmented Lagrangian method and solve by conjugate gradient on a spline basis. We present details of the method and a 2D application to data from LA RISTRA in the western United States. In order to be effective in 2D, we require teleseismic phases arriving at the array at a broad sweep of incidence angles. With the Andean and Aleutian subduction zones along the strike of the array between 35° and 85° epicentral distance, LA RISTRA provides the ideal illumination for a tomographic inversion.