# Waveform inversion for orthorhombic anisotropy with P waves: feasibility and resolution

Waveform inversion for orthorhombic anisotropy with P waves: feasibility and resolution Summary Various parametrizations have been suggested to simplify inversions of first arrivals, or P waves, in orthorhombic anisotropic media, but the number and type of retrievable parameters have not been decisively determined. We show that only six parameters can be retrieved from the dynamic linearized inversion of P waves. These parameters are different from the six parameters needed to describe the kinematics of P waves. Reflection-based radiation patterns from the P–P scattered waves are remapped into the spectral domain to allow for our resolution analysis based on the effective angle of illumination concept. Singular value decomposition of the spectral sensitivities from various azimuths, offset coverage scenarios and data bandwidths allows us to quantify the resolution of different parametrizations, taking into account the signal-to-noise ratio in a given experiment. According to our singular value analysis, when the primary goal of inversion is determining the velocity of the P waves, gradually adding anisotropy of lower orders (isotropic, vertically transversally isotropic and orthorhombic) in hierarchical parametrization is the best choice. Hierarchical parametrization reduces the trade-off between the parameters and makes gradual introduction of lower anisotropy orders straightforward. When all the anisotropic parameters affecting P-wave propagation need to be retrieved simultaneously, the classic parametrization of orthorhombic medium with elastic stiffness matrix coefficients and density is a better choice for inversion. We provide estimates of the number and set of parameters that can be retrieved from surface seismic data in different acquisition scenarios. To set up an inversion process, the singular values determine the number of parameters that can be inverted and the resolution matrices from the parametrizations can be used to ascertain the set of parameters that can be resolved. Fourier analysis, Inverse theory, Waveform inversion, Body waves, Seismic anisotropy, Wave scattering and diffraction 1 INTRODUCTION Reservoir delineation and monitoring are among the major objectives of modern seismic imaging and inversion. The elastic properties of most hydrocarbon reservoirs can be well approximated using orthorhombic anisotropy (Schoenberg & Helbig 1997; Tsvankin 1997; Ivanov & Stovas 2016) on the scale of common exploration seismic wavelengths. Systems of orthogonal fractures or vertical layering within a reservoir can be described by macroscale orthorhombic media; see Fig. 1(a). Seismic imaging of a reservoir area aims to obtain the effective orthorhombic parameters, which are then converted into the microscale properties of the imaged reservoir using rock-physics modeling. The accuracy of the orthorhombic parameter estimates is crucial to a reliable conversion. Figure 1. View largeDownload slide Orthorhombic reservoir. (a) Anisotropic reservoir (Tsvankin 1997) in isotropic background. (b) Elements of stiffness matrix in Voigt notation. If the media is orthorhombic then only the red cells are non-zero in its main axes, monoclinic with horizontal symmetry plane by red and blue cells. Figure 1. View largeDownload slide Orthorhombic reservoir. (a) Anisotropic reservoir (Tsvankin 1997) in isotropic background. (b) Elements of stiffness matrix in Voigt notation. If the media is orthorhombic then only the red cells are non-zero in its main axes, monoclinic with horizontal symmetry plane by red and blue cells. Conventionally, a moveout velocity analysis of the first arrivals, that is, P waves, is used to provide an estimate of the orthorhombic medium parameters (Grechka & Tsvankin 1999). P-wave velocity analysis (Tsvankin 1997; Grechka & Tsvankin 1999; Stovas 2017) defines the limited resolution of classic ray-based methods in orthorhombic media. Recently, interest has grown in the inversion of dynamic wavefield components due to their higher resolution potential, despite the additional acquisition and computational requirements (e.g. Virieux & Operto 2009; Anikiev et al.2014; Kurzmann et al.2016). Iterative imaging techniques, such as full-waveform inversion (FWI) and least-squares reverse-time migration, provide high-resolution estimates of the subsurface parameters, but require more computational resources and above-average signal-to-noise ratios in the data. Marine data sets are usually able to meet the signal-to-noise ratio requirement with air guns, which excite P waves to illuminate the subsurface. To reduce computational costs, most industrial-scale 3-D FWI implementations of exploration inversion rely on pseudo-acoustic (qP wave) physics (e.g. Alkhalifah 2000; Gholami et al.2013a,b; Kamei et al.2013). However, the reflection coefficients cannot be well approximated by pseudo-acoustic modeling, especially in acquisitions limited by short offsets (e.g. Raknes & Arntsen 2014; Ovcharenko et al.2018). Since elastic anisotropic waveform inversions have higher computational costs to their isotropic counterparts, they only recently became the focus of attention (Köhn et al.2015; Podgornova et al.2015; Owusu et al.2016; Bergslid et al.2015; Oh & Alkhalifah 2016; Kamath & Tsvankin 2016). Feasibility and resolution studies for 3-D inversions are still computationally expensive, therefore they are often performed either in 2-D (Podgornova et al.2015) or with very few simultaneous (encoded) shots (Köhn et al.2015). Inversion for the full set of orthorhombic parameters in the whole subsurface model increases the likelihood of overfitting the data. If the orthorhombic anisotropy is restricted to the reservoir area of the subsurface model and the background is assumed to be isotropic, then the number of inverted parameters decreases and, thus, the inversion becomes more robust. When the anisotropy is restricted to the area of the reservoir, the scattering is well approximated by Born scattering in the isotropic background (Hudson & Heritage 1981). In Born approximation, scattering of plane waves can be fully characterized by a scattering function (Eaton & Stewart 1994), which is a function of the incident plane wave parameters and the perturbation of the elastic parameters. If the parameter perturbation and the types of incident and scattered waves are fixed, then the scattering function is reduced to a diffraction-based radiation pattern (Gholami et al.2013b; Alkhalifah & Plessix 2014). Finally, if the incident and scattered angles in the diffraction-based radiation pattern obey Snell’s law for horizontal reflectors, then the diffraction-based radiation pattern becomes a reflection-based radiation pattern (Alkhalifah & Plessix 2014). With the assumption of an isotropic background, the resolution can be investigated visually using the reflection-based radiation patterns (Alkhalifah & Plessix 2014; Oh & Alkhalifah 2016) and quantitatively using the singular value decomposition (SVD) analysis of those patterns (Eaton & Stewart 1994; de Hoop et al.1999; Podgornova et al.2015; Kazei & Alkhalifah 2017). Eaton & Stewart (1994) were the first to discover the importance of the scattering function in anisotropic media and established a resolution analysis for vertically incident qP waves in vertically transversally isotropic (VTI) media. de Hoop et al. (1999) derived a similar analysis for orthorhombic parameter perturbations that focused on reflections in an isotropic background. Reflection-based radiation patterns were used to analyse the ability of anisotropic elastic FWI to estimate VTI parameters using only P waves in various parametrizations (Gholami et al.2013b; Alkhalifah & Plessix 2014). Podgornova et al. (2015) proved that only three VTI parameters are invertible from P–P waves under the assumption of weak anisotropy, by using SVD analysis to investigate the FWI resolution possibilities with P–P and P–SV waves in Thomsen’s well-known parametrization for VTI media (Thomsen 1986). Thomsen’s parametrization was generalized to orthorhombic media by Tsvankin (1997) and rewritten with deviation parameters by Masmoudi & Alkhalifah (2016) and Oh & Alkhalifah (2016) to obtain a hierarchical anisotropy parametrization. Oh & Alkhalifah (2016) came to the conclusion that their parametrization minimizes trade-offs between parameters and, for P–P waves in most acquisition scenarios, allows the inversion of four out of the 10 independent elastic orthorhombic parameters. The traveltimes of P waves contain information on only six linear combinations of these parameters (Tsvankin 1997). A numerical resolution study (Köhn et al.2015) has shown that, with P-wave dominant sources, three-component measurements and perfect illumination, inversion for all 10 orthorhombic parameters may be feasible, but the general acquisition scenarios were not addressed. Here, we investigate the relation between the acquisition and the number of invertible parameters assuming an explosive source and an early arrival amplitude, namely, P waves recorded by hydrophones. In Section 2, we provide and discuss the derivation of an anisotropic scattering function and clarify its relation to radiation patterns. In Section 3, we map 3-D radiation patterns into the wavenumber domain and discuss the resolution qualitatively. In Section 4, we introduce the effective angle of illumination concept and use it to compare a recent hierarchical parametrization suggested by Masmoudi & Alkhalifah (2016) for acoustic case and by Oh & Alkhalifah (2016) for elastic with a classic stiffness matrix parametrization of orthorhombic media. In Section 5, we apply an SVD analysis to the sensitivities of P–P scattered waves to derive the number of parameters that can be inverted in different acquisition scenarios. 2 LINEARIZED SCATTERING, RADIATION PATTERNS The linearization of the wave equation with respect to parameter perturbation provides the crucial first-order sensitivity of the data to the model. The analysis of such sensitivity using plane wave and SVD provides us with the scattering radiation patterns, and insights into the model invertibility considering the limited illumination provided by conventional recording of seismic data. 2.1 The Born approximation Most methods for dynamic inversion of seismic wavefields, including linearized reverse-time migration and every iteration of FWI, rely on the Born approximation. Thus, the sensitivity of the wavefields, which defines the resolution of these methods, can be evaluated within the Born approximation. We analyse the resolution of dynamic inversions for orthorhombic parameters relying on the linearized scattering of an incident wavefield U0 on a perturbation of the density (δρ(x)) and stiffness tensor (δckjpq(x)) parameters. The frequency-domain Born approximation for the scattered wavefield δU in an anisotropic elastic medium is (e.g. Hudson & Heritage 1981; Beylkin & Burridge 1990; Shaw & Sen 2004)   \begin{eqnarray} \delta U_i (\mathbf {x}_s,\mathbf {x}_\mathbf {g},\omega ) &=& \int \limits _{V} \left(\delta \rho (\mathbf {x}) u^0_k \omega ^2\right.\nonumber\\ \left. - \delta c_{kjpq}(\mathbf {x}) \frac{\partial U^0_p}{\partial x_q} \frac{\partial }{\partial x_j} \right)G_{ik}(\mathbf {x}_g, \mathbf {x}) {\rm d}\mathbf {x}, \end{eqnarray} (1)where ω is the angular frequency, xs and xg are the coordinates of the source and receiver, respectively and Gik(xg, x) is the Green’s tensor for the background medium. Integration is performed over the volume of the domain where the perturbation is located. Born approximation (eq. 1) is valid for relatively small   \begin{eqnarray} \frac{|\delta \mathbf {m}|}{|\mathbf {m}|} << 1 \end{eqnarray} (2)and local   \begin{eqnarray} \frac{\omega d}{V_0}\frac{|\delta \mathbf {m}|}{|\mathbf {m}|} << 1 \end{eqnarray} (3)perturbations, where m is the vector of unperturbed elastic parameters, V0 is the smallest background velocity (Vs in the isotropic case) and d is the diameter of the perturbation. The above assumptions are best met in the case of low-contrast point scatterers, yet the Born approximation does not imply specific shape restrictions on the scattering perturbation. For instance, it is also valid for the perturbations within thin low-contrast layers and near-normal wave incidence angles (Shaw & Sen 2004). To obtain a simplified expression for the scattering of plane waves on an arbitrary perturbation, we first manipulate the right-hand side of eq. (1). Next, we discuss the cases of point scatterers and horizontal reflectors. Finally, we examine the general case through spectral sensitivities. For simplicity, we consider the far-field approximation from both the source and receiver sides, with an isotropic background. 2.2 Far-field plane-wave scattering Most FWI case studies focus on marine data sets due to their high signal-to-noise ratios. In marine acquisition, only P waves propagate through water and, thus, they dominate the incident wavefields in the target area of exploration. Therefore, we consider only P waves using the following far-field plane-wave approximation for the incident wavefield:   \begin{eqnarray} \mathbf {U}^0(\mathbf {x}_s,\mathbf {x},\omega ) &=& \alpha _s e^{i\frac{\omega }{v_p}||\mathbf {x}- \mathbf {x}_s||}\mathbf {s}\simeq \alpha _s e^{i\frac{\omega }{v_p}\mathbf {s}\cdot (\mathbf {x}- \mathbf {x}_g)}\mathbf {s},\, \alpha _s\nonumber\\ &=& \frac{1}{4\pi \rho v_p^2||\mathbf {x}- \mathbf {x}_s||} \end{eqnarray} (4)where ρ and vp are the background density and longitudinal wave velocity, respectively, and s is a unit vector pointing towards the source from the target area. The far-field part of the Green’s tensor in eq. (1) for isotropic background consists of P, SV and SH waves:   \begin{eqnarray} \mathbf {G}= \alpha _g (\mathbf {G}_P + \mathbf {G}_{SV} + \mathbf {G}_{SH}),\, \alpha _g = \frac{1}{4\pi \rho v_p^2||\mathbf {x}- \mathbf {x}_g||}. \end{eqnarray} (5)The coefficient αg is dropped in our resolution analysis, since it is related to geometrical spreading and can be effectively compensated using pseudo-Hessian techniques. Different components of the Green’s function are approximated locally by plane waves. For the first component GP, the approximation leads to:   \begin{eqnarray} \mathbf {G}_P = e^{i\frac{\omega }{v_p}||\mathbf {x}- \mathbf {x}_g||} \mathbf {g}\mathbf {g}\simeq e^{i\frac{\omega }{v_p}\mathbf {g}\cdot (\mathbf {x}- \mathbf {x}_g)}\mathbf {g}\mathbf {g} \end{eqnarray} (6)which is a plane P wave. (gg)ik = gigk is the outer product of vector g, pointing towards receiver from the perturbation, with itself. For the second and third components, GSV and GSH, similar expressions can be provided (Snieder 2002), yet we keep them out of the scope of this paper in order to focus on the early P-wave arrivals. By substituting eqs (4) and (5) into eq. (1) and projecting on g, we obtain the approximate amplitude of the scattered P–P waves:   \begin{eqnarray} \delta U_{PP} \equiv \frac{\delta \mathbf {U} \cdot \mathbf {g}}{\alpha _{PP}} \simeq \int \limits _{V}e^{i\mathbf {K}_{PP}\cdot \mathbf {x}} (\mathbf {s}\cdot \mathbf {g}\delta \rho - \mathbf {s}\mathbf {s}: \delta \mathbf {c}:\mathbf {g}\mathbf {g}) {\rm d}\mathbf {x} \end{eqnarray} (7)  \begin{eqnarray} =\mathbf {s}\cdot \mathbf {g}\int \limits _{V}e^{i\mathbf {K}_{PP}\cdot \mathbf {x}} \delta \rho {\rm d}\mathbf {x}- \mathbf {s}\mathbf {s}: \int \limits _{V}e^{i\mathbf {K}_{PP}\cdot \mathbf {x}} \delta \mathbf {c}{\rm d}\mathbf {x}: \mathbf {g}\mathbf {g}, \end{eqnarray} (8)where   \begin{eqnarray} \alpha _{PP} = \frac{\omega ^2}{v^2_p} \alpha _s \alpha _g e^{i\frac{\omega }{v_p}(-\mathbf {x}_s\mathbf {s}-\mathbf {x}_g\mathbf {g})}, \end{eqnarray} (9)  \begin{eqnarray} \mathbf {K}_{PP} = \frac{\omega }{v_p}(\mathbf {s}+\mathbf {g}). \end{eqnarray} (10)Our further resolution analysis is based on eq. (8), which serves as the simplified monochromatic scattering approximation. 2.3 Scattering function A monochromatic plane wave is determined by its phase, polarization and propagation direction. Since the phase is dropped, scattering depends only on the incident plane wave direction s, scattered plane wave direction g, their respective polarizations ς and ξ, the frequency ω, and on the perturbation of the parameters,   \begin{eqnarray} \delta \mathbf {m}= (\delta \rho (\mathbf {x}), \delta \mathbf {c}(\mathbf {x}))^T. \end{eqnarray} (11)Namely, the scattered wavefield from a point perturbation δm = δ(x)(δρ, δC)T in the first-order Born approximation is proportional to the scattering function (Beylkin & Burridge 1990; Eaton & Stewart 1994; Shaw & Sen 2004):   \begin{eqnarray} \delta U (\mathbf {s},\mathbf {g}, \boldsymbol {\varsigma }, \boldsymbol {\xi }, \delta \mathbf {m}) \propto R\equiv \varsigma _k \cdot \xi _k \delta \rho - s_i \varsigma _j \delta c_{ijkl} g_k \xi _l \end{eqnarray} (12)  \begin{eqnarray} = \boldsymbol {\varsigma }\cdot \boldsymbol {\xi }\delta \rho - \mathbf {s}\boldsymbol {\varsigma }: \delta \mathbf {c}: \mathbf {g}\boldsymbol {\xi }. \end{eqnarray} (13)The scattering function R is independent of the background medium, which can have up to triclinic anisotropy complexity. As there are too many possible cases of polarization and propagation relations in general anisotropic background, we consider scattering in an isotropic background medium. 2.4 Classic diffraction-based scattering radiation pattern—$$\mathcal {D}$$ For a given conversion type (P–P, P–SH, etc.), the scattered wavefield becomes a function of two vectors, s and g, called the scattering radiation pattern:   \begin{eqnarray} \mathcal {D}_{WI-WS, \delta \mathbf {m}}(\mathbf {s},\mathbf {g}) \equiv \frac{ R(\mathbf {s}, \mathbf {g}, \boldsymbol {\varsigma }(\mathbf {s}, WI), \boldsymbol {\xi }(\mathbf {g}, WS), \delta \mathbf {m})}{||\delta \mathbf {m}||}. \end{eqnarray} (14)Here, WI denotes the type of incident, and WS denotes the type of scattered waves (P − P, P − Si, Si − P or Si − Sj), for a given parameter perturbation δm. The radiation pattern characterizes the scattering of certain types of plane waves on a point-like perturbation of the parameters. For vertically propagating incident P waves, the scattering radiation patterns become functions of two real arguments, that is, g vector components, and can be presented as a 2-D plot (Eaton & Stewart 1994). We show the diffraction radiation pattern for P–P scattering on a density perturbation in Fig. 2(a) for the vertical incidence direction (e.g. Tarantola 1986). Figure 2. View largeDownload slide (a) The diffraction pattern for density perturbation. (b)  Narrow azimuth—ϕ acquisition and reflection from the reservoir with scattering angle θ. (c)  The reflection based radiation pattern. and (d)  its remapping into the spectral domain. The normalized wavenumber ||K||∝cos θmax. In all the patterns, the amplitude sign is shown through colour, ‘+’ corresponds to positive reflection coefficient—no phase shift after scattering, ‘−’—shows the flip of the phase—negative reflection/scattering coefficient. Figure 2. View largeDownload slide (a) The diffraction pattern for density perturbation. (b)  Narrow azimuth—ϕ acquisition and reflection from the reservoir with scattering angle θ. (c)  The reflection based radiation pattern. and (d)  its remapping into the spectral domain. The normalized wavenumber ||K||∝cos θmax. In all the patterns, the amplitude sign is shown through colour, ‘+’ corresponds to positive reflection coefficient—no phase shift after scattering, ‘−’—shows the flip of the phase—negative reflection/scattering coefficient. 2.5 Reflection-based radiation patterns—$$\mathcal {R}$$ Classic scattering radiation patterns have focused on diffractions from point scatterers. However, seismic surveys are commonly dominated by reflections (Fig. 2b), which prompted the development of linearized reflection coefficients (e.g. Rüger 1997) and, more recently, reflection-based radiation patterns (Shaw & Sen 2004; Gholami et al.2013b; Alkhalifah & Plessix 2014). The relation between the reflection patterns and linearized coefficients is well explained by Shaw & Sen (2004). In 3-D scattering experiments, scattering functions depend on four scalar parameters defined by the unit vectors, s and g. For a plane reflector, the scattered wave direction is determined by the incident direction and vice versa through Snell’s law. The dip and azimuth of the reflector determine how the s is tied to the g in particular, reducing the number of real variables that determine scattering from four to two:   \begin{eqnarray} \mathcal {R}_{WI-WS, \delta \mathbf {m}}(\mathbf {g}) \equiv \mathcal {D}_{WI-WS, \delta \mathbf {m}}(\mathbf {s}(\mathbf {g}),\mathbf {g}). \end{eqnarray} (15)The amplitudes of the reflections from perturbations of isotropic parameters do not depend on the reflector orientation; therefore, the reflection-based radiation patterns are the same for all reflector orientations. In anisotropic cases, the reflection-based radiation patterns depend on the dip and azimuth angles. In most seismic cases, the media is dominated by horizontal layering; for this reason, Gholami et al. (2013b) and Alkhalifah & Plessix (2014) focused on the horizontal reflectors. If we assume P–P scattering on a horizontal reflector, then Snell’s law can be expressed as (si and gi are projections of vectors s and g, respectively)   \begin{eqnarray} s_1&=& -g_1 = \sin \theta \cos \phi ,\, s_2=-g_2=\sin \theta \sin \phi ,\nonumber\\ s_3&=&g_3=\cos \theta /2. \end{eqnarray} (16)The radiation patterns become functions of two real variables, ϕ and θ, and can be plotted. The reflection-based radiation pattern for density and P–P scattering is represented in Fig. 2(c). We show the reflection-based radiation patterns for other cases mapped from the ϕ, θ domain into wavenumber domain (ϕ, kz = 2cos θ/2) as spectral sensitivities, to directly access the resolution in the normalized wavenumber domain. 3 SPECTRAL SENSITIVITIES Spectral sensitivities allow us to consider a broad class of anomalies that are small enough that we can neglect the second-order scattering. The relation between the wavenumbers illuminated in a perturbation and the directions of incident and scattered waves in a monochromatic plane-wave experiment has been recognized and successfully applied in the X-ray tomographic method of Ewald (1969) spheres. Devaney (1984) applied this method to geophysical diffraction tomography resolution quantification in the wavenumber domain. The technique was later applied in acoustic FWI scenarios for surface seismic acquisition (Mora 1989) and vertical seismic profiling (Wu & Toksöz 1987). Kazei et al. (2013a) introduced amplitude of the spectral sensitivity of data to the perturbations of various background media, and applied it to models with diving waves and multiples (Kazei et al.2015, 2013b). The elastic isotropic background case with VTI perturbations was considered by Podgornova et al. (2015). Here, we show the applicability of this technique to the general case of anisotropic perturbations of isotropic media. The possibility of recovering a set of parameters is defined by the sensitivity of the scattered wavefield to these parameters. Under the plane-wave assumption, the resolution is straightforward for the wavenumber-domain representation of perturbation. The concept of spectral sensitivities plays an important role in our analysis. 3.1 Definition The right-hand side of eq. (8) is, essentially, a Fourier transform of the perturbation and can be rewritten to reveal the nature of scattering in a far-field approximation as follows:   \begin{eqnarray} \delta U_{PP}(\mathbf {s},\mathbf {g},\omega ) = \mathbf {s}\cdot \mathbf {g}\delta \hat{\rho }(\mathbf {K}_{PP}) - \mathbf {s}\mathbf {s}: \delta \hat{\mathbf {c}}(\mathbf {K}_{PP}) :\mathbf {g}\mathbf {g}. \end{eqnarray} (17)If we consider the spectrum of perturbation parameters (δm(K)), then we see that:   \begin{eqnarray} \delta U_{PP} = \mathcal {R}_{PP} (\mathbf {s}, \mathbf {g}) ||\delta \mathbf {m}(\mathbf {K}_{PP})||, \end{eqnarray} (18)which allows us to represent radiation patterns in the normalized wavenumber domain:   \begin{eqnarray} \mathcal {S}_{\delta \mathbf {m}} (\mathbf {k}_{PP}, \varphi ) \equiv \mathcal {R}(\mathbf {s}(\mathbf {k}_{PP},\varphi ), \mathbf {g}(\mathbf {k}_{PP},\varphi )), \end{eqnarray} (19)where   \begin{eqnarray} \mathbf {k}_{PP} = \mathbf {K}_{PP}/k_0 = \mathbf {s}+\mathbf {g},\, k_{0} = \frac{\omega }{v_p}, \end{eqnarray} (20)Eq. (19) defines the spectral sensitivity for parameter δm. It is more convenient than the reflection-based radiation patterns as it quantifies the resolution for different normalized wavenumbers kPP directly. We normalize all of the wavenumbers, dividing them by ω/Vp, acknowledging that the sensitivities scale linearly with frequency in homogeneous backgrounds (Kazei et al.2013a). We focus on the inversion for vertical wavenumbers kPP = (0, 0, kz), or horizontally layered structures, which are dominant in most geological scenarios, returning us to the remapped reflection-based radiation patterns in a more convenient representation for resolution analysis domain. These remapped patterns are functions of two variables (kz, ϕ) and, therefore, can be easily plotted and examined. Snell’s law gives the following expressions in the case in which the incident and scattered directions are functions of ϕ and kz:   \begin{eqnarray} s_1&=& -g_1 = \sqrt{1-\frac{k_z^2}{4}}\cos \phi ,\, s_2=-g_2=\sqrt{1-\frac{k_z^2}{4}}\sin \phi ,\nonumber\\s_3&=&g_3=k_z/2. \end{eqnarray} (21)In Fig. 2(d), we show the spectral sensitivity for a perturbation in density when all of the elastic constants are kept fixed. There is no scattering at intermediate reflection angles (Fig. 2a, so the intermediate spatial wavenumbers for density do not contribute to the scattered wavefield and cannot be reconstructed. However, low and high wavenumbers can be reconstructed. This is different from the standard inversion for density where Vp is fixed and only the high wavenumbers show significant scattering. In the next section, we analyse the scattering of all the elastic constants and try to estimate the null space of the inversion. 3.2 Spectral sensitivity patterns Cij, ρ parametrization We examine the remapped radiation patterns with the newly introduced representation in Fig. 3. Tsvankin (1997) established that the traveltimes of P waves in orthorhombic media depend on six parameters. This is the number of parameters that kinematic methods can resolve in perfect illumination conditions. One might expect that the inversion of the amplitudes could recover more parameters, since more attributes of the wavefields are taken into account. However, the visual analysis of the P–P radiation patterns does not provide detailed information on how the Cij parameters are coupled. Figure 3. View largeDownload slide The P–P spectral sensitivity patterns for various Cij perturbations. The azimuth and Kz values are mapped the same way as for density which is shown in Fig. 2(d). It is clear that the sensitivity to all the orthorhombic parameters is non-zero if we start from isotropic background. Yet we can note some coupling effects already. C13 is very similar to C55, C12—to C66. Figure 3. View largeDownload slide The P–P spectral sensitivity patterns for various Cij perturbations. The azimuth and Kz values are mapped the same way as for density which is shown in Fig. 2(d). It is clear that the sensitivity to all the orthorhombic parameters is non-zero if we start from isotropic background. Yet we can note some coupling effects already. C13 is very similar to C55, C12—to C66. To analyse the scattering, we generate the spectral sensitivities for single frequency, infinite offset data. According to eq. (10), for single frequency infinite offset data, P waves cover vertical wavenumbers from 0 to $$\frac{2\omega }{V_p}$$. In Fig. 3, we see four non-orthorhombic elements with non-zero sensitivities. These non-zero sensitivities can be used to invert for monoclinic medium parameters with a horizontal symmetry plane. At the same time, we admit that the sensitivity to tilt-related parameters (rotations of perturbations around a horizontal axis) is zero in an isotropic background; therefore, the inversion for tilt from P waves does not work for pure vertical wavenumbers in isotropic backgrounds. Sensitivities, that is, remapped reflection-based radiation patterns, for different Cij parameters provide some information about which parameters cannot be inverted for; yet they provide almost no information about the coupling of parameters, as there are usually more than two parameters involved in the coupling, as we show in the next section. 3.3 Sensitivities in a hierarchical parametrization While necessary for the characterization of carbonate reservoirs, orthorhombic anisotropy is not observed in vast areas of the subsurface. An easy way to set up an inversion with different anisotropy restrictions in different domains of the subsurface is to split the parameters into three groups using a hierarchical parametrization based on the type of anisotropy that they describe. Such a parametrization, introduced by Oh & Alkhalifah (2016), uses Vp, Vs, ρ —for the isotropic part of the medium, where Vp and Vs are the velocities of P and S waves in isotropic media, respectively, and ρ is the density. In general orthorhombic media, the Vp and Vs considered here are the velocities with wave polarizations along the x1 horizontal axis. ε1, η1, γ1 —for the VTI part of the anisotropy, where ε1 and γ1 is the relative difference in velocities between vertically and horizontally propagating P and SH waves, respectively, and η1 is the Alkhalifah–Tsvankin anellipticity parameter in [x1, x3] plane. εd, ηd, γd, δ3 —are parameters that define deviation of the medium from VTI anisotropy. They are non-zero only if the medium has fully orthorhombic complexity. ηd = η2 − η1, where η2 is the Alkhalifah–Tsvankin (Alkhalifah 2000) anellipticity parameter in [x1, x3] plane. Similarly, γd = γ2 − γ1 and εd = ε2 − ε1 define differences in VTI anisotropy in [x2, x3] and [x1, x3] planes. δ3 is the well-known Tsvankin–Thomsen parameter which is zero in case the medium is VTI. Radiation patterns for any parametrization can be obtained through the chain rule   \begin{eqnarray} \boldsymbol {\mathcal {R}}_{\mathbf {m}, WT} = \frac{\partial \mathbf {C}}{\partial \mathbf {m}}:\boldsymbol {\mathcal {R}}_{\mathbf {C}, WT}. \end{eqnarray} (22)The Jacobian $$\frac{\partial \mathbf {c}}{\partial \mathbf {m}}$$ (eq. B15) reveals the actual perturbations of the Cij parameters corresponding to individual perturbations in the new parametrization. Therefore, we also represent the partial derivatives of Cij parameters graphically in Fig. 4(a) to understand the analysis better. Figure 4. View largeDownload slide (a) The partial derivatives of the Cij parameters (same scheme as in Fig. 1 b is used) in the Cij, ρ parametrization with respect to the Oh & Alkhalifah (2016) parameters. (b)  The spectral sensitivities of the new parameters. Figure 4. View largeDownload slide (a) The partial derivatives of the Cij parameters (same scheme as in Fig. 1 b is used) in the Cij, ρ parametrization with respect to the Oh & Alkhalifah (2016) parameters. (b)  The spectral sensitivities of the new parameters. Spectral sensitivities in the new parametrization elucidate some features of scattering that were not visible in the Cij, ρ parametrization. For example, the resemblance of radiation patterns for density and ε1 gives a hint that there is a trade-off between these parameters, which was hidden in several radiation patterns before. At the same time our knowledge is still limited by the choice of parametrization as we cannot analyse the coupling between three or more parameters. 4 EFFECTIVE ANGLE OF ILLUMINATION To perform a successful multiparameter inversion, we will need to retrieve all parameters in the whole domain of interest to avoid the dreaded null space. The capability of inverting of decoupled parameters can be investigated for a given vertical wavenumber, separately. The scattering of monochromatic plane waves for different wavenumbers in a perturbed medium is determined by the scattering angle–wavenumber relation (Mora 1989; Alkhalifah 2016; Kazei et al.2016) using a locally homogeneous medium assumption. Utilizing this property, we are going to examine the ability of the inversion in decoupling parameters for particular vertical wavenumber Kz. For single frequency and single azimuth data, there is a bijective relation between the wavenumber illuminated (two real parameters) and the incidence and scattered directions (two scalar angles), which means that the scattered wavefield maps directly to the parameter at this vertical wavenumber. Therefore, in principle, single frequency, single component data can resolve only one parameter. Unless we impose some relations between parameters at different wavenumbers with, for example, some type of regularization. In order to decouple VTI media parameters for a given wavenumber, we need a set of frequencies (Podgornova et al.2015) and azimuths contributing to this wavenumber illumination. While for all vertical wavenumbers, we assume the azimuth coverage is the same there might be additional angle-offset limitations due to frequency content of the data. These limitations are summarized into the effective angle of wavenumber illumination concept in this section. The scattering of plane waves can be expressed in the following system of equations:   \begin{eqnarray} \delta \mathbf {U}(\omega _i,\theta _j,\varphi _k) \propto \mathcal {R}(\theta _j,\varphi _k) \delta \mathbf {m}(\mathbf {K}(\theta _j,\omega _i)), \end{eqnarray} (23)which can be regrouped for a given wavenumber K, using the angle–wavenumber relation:   \begin{eqnarray} \mathbf {u}(\omega _i(\mathbf {K}),\theta _i(\mathbf {K}),\varphi _k) \propto \mathcal {R}(\theta _i(\mathbf {K},\omega ),\varphi _k) \delta \mathbf {m}(\mathbf {K}). \end{eqnarray} (24)Thus, data for each vertical wavenumber of the perturbation are available from different frequencies ωi, which gives us different opening angles θi = (K, ωi) for illumination of the wavenumber K. Having different angles of illumination for the same wavenumber K allows us to distinguish parameters at the scale given by this wavenumber through variation in the scattering patterns of parameters. In Fig. 5, we show that the set of angles ‘active’ in a radiation pattern for a given wavenumber are defined by the offset as well as by the frequency content of the signal. Figure 5. View largeDownload slide The same wavenumber can be illuminated from different opening angles and by different wave modes. Each angle availability is determined by the extent of the angle in the acquisition and existence of the appropriate frequencies to scale the source and geophone vectors. Figure 5. View largeDownload slide The same wavenumber can be illuminated from different opening angles and by different wave modes. Each angle availability is determined by the extent of the angle in the acquisition and existence of the appropriate frequencies to scale the source and geophone vectors. We examine the effective angle of illumination concept on P–P waves as the relations are simple in this case. According to eq. (20) and Fig. 5, the angles that contribute to the illumination of a given vertical wavenumber can be restricted by both the frequency and the aperture. If the recorded band of frequencies is narrow then the angles of illumination are determined by the bandwidth of the data. In this case, there is no benefit received from an aperture increase in decoupling the parameters; that is, there is no reason to increase the aperture while the frequency content of the data remains the same. For example, in $$K_z\in [\frac{2\omega }{v_p}\cos \frac{\theta _{{\rm max}}}{2};\frac{2\omega }{v_p}]$$ (the red region in Fig. 10a), the angles of illumination are not affected by the aperture at all. If a wide band of frequencies is recorded the effective angle of illumination is generally determined by the aperture. Namely, for the range of wavenumbers Kz between $$\frac{2\omega }{v_{{\rm min}}}$$ and $$\frac{2\omega _{{\rm max}}}{V_p}\cos \theta _{{\rm max}}$$, the full offset length can be exploited for inversion and the decoupling of parameters, that is, in the green area in Fig. 6(b), the resolution is defined by the aperture only. The effective angle of illumination concept suggests that the data frequency bandwidth and the aperture should be increased simultaneously to achieve better decoupling, needed for the inversion. In the next section, we assume that the bandwidth is sufficient and consider the resolution for the green area in Fig. 6(b). Figure 6. View largeDownload slide For a given vertical wavenumber, angles that contribute to its illumination can be restricted by frequency as well as by aperture. (a) A narrow band of frequencies. In the red area, the resolution of wavenumbers is defined by the frequency content. (b) A wide band of frequencies. In the green area, the resolution is defined by the illumination angles. Wavenumbers between $$\frac{2\omega _{{\rm min}}}{v}$$ and $$\frac{2\omega _{{\rm max}}}{v}\cos \theta _{{\rm max}}$$ can be inverted using the whole offset length. Figure 6. View largeDownload slide For a given vertical wavenumber, angles that contribute to its illumination can be restricted by frequency as well as by aperture. (a) A narrow band of frequencies. In the red area, the resolution of wavenumbers is defined by the frequency content. (b) A wide band of frequencies. In the green area, the resolution is defined by the illumination angles. Wavenumbers between $$\frac{2\omega _{{\rm min}}}{v}$$ and $$\frac{2\omega _{{\rm max}}}{v}\cos \theta _{{\rm max}}$$ can be inverted using the whole offset length. 4.1 Sensitivity matrices $$\mathbb {A}_{P-P}$$ – assembly The main outcome of the previous subsection is that particular wavenumber illumination can be limited not only by the lack of aperture in the data, but also by the band of the signal of the source and geophones. In this subsection, we discuss the matrix form of the simplified inversion problem for the vertical wavenumber K from the P–P scattering mode to prepare ourselves for quantitative resolution analysis. For an arbitrary perturbation of parameters, the scattered wavefield can be expressed in the wavenumber domain as   \begin{eqnarray} \delta U_{PP} (k_z, \phi , \omega (k_z, K_z)) \propto \mathbf {A}_{P-P}^T (k_z, \phi ) (\delta \mathbf {m}(\mathbf {K})), \end{eqnarray} (25)or in the vectorized form as   \begin{eqnarray} \delta \mathbf {U}_{PP}(K_z) \propto \mathbb {A}_{P-P} \delta \mathbf {m}(K_z), \end{eqnarray} (26)where kz ∈ [2cos (θmax/2), 2]. Vector A(kz, ϕ) represents the set of reflection-based radiation patterns remapped into the normalized wavenumber domain. For example, element A1(kz, ϕ) (Fig. 2d) is the radiation pattern of the density, while the radiation patterns for the 21 elastic constants (the total number of Voigt elastic parameters is 22 – 21 elastic constants and density), are shown in Fig. 3, which is ordered to correspond to the Voigt notation (Fig. 1b). We can clearly see couplings in pairs (resemblance in patterns) : C12 and C66, C13 and C55 and C23 and C44. We start by analysing the linearized problem that needs to be solved in order to reconstruct a low vertical wavenumber of perturbation, since low wavenumbers are considered the hardest and most essential wavenumbers for FWI to converge. First: let us assume that the frequency content is sufficient, and the illumination of the wavenumber is limited by the aperture (Fig. 6b). Then, there will be angles of illumination [0; θmax] available for a given wavenumber K, where θmax is determined by the aperture. First, sticking with the elastic coefficients and density parametrization, we seek to understand the trade-off and extract the linear combinations that expose such trade-off. Similarities of certain radiation patterns in Fig. 3 (e.g. C12 and C66) suggest the presence of a null space. The rank of the matrix $$A_{ik} = \mathcal {R}(\theta _i(\mathbf {K},\omega ),\varphi _k)$$ is equal to the number of parameters resolvable for a given wavenumber K value. The dimensionality of the null space is equal to the total number of parameters minus the rank of the matrix. Note that the part of the radiation patterns that plays role in this analysis is commonly attributed to high resolution in other studies (e.g. Alkhalifah & Plessix 2014), because the wavenumbers obtained from small offsets are higher than those from far offsets for a given frequency. The frequency is not an issue since we consider a wide band. The low frequencies in the data are useful for monoparameter FWI ten Kroode et al. (2013), yet for multiparameter they are necessary to resolve trade-off between parameters. When there is a lack of high frequencies, the effective angle range for inversion for large wavenumbers activates the high wavenumber part of the spectral sensitivity. Thus, inversion for low wavenumbers with limited offset data is similar to inversion for high wavenumbers with a limit in the highest frequency available in the data. Second: we assume that lack of low frequencies limits the resolution of low wavenumbers in the model. The offset is assumed to be effectively infinite in this case, for example, diving waves illuminate the whole region of interest (Kazei et al.2013b). Diving waves, propagating horizontally at the depth of perturbation, activate the ‘transmission’ part of the radiation patterns with wide opening angles (Alkhalifah & Plessix 2014), equivalent to low vertical wavenumbers in spectral sensitivities. The second case is common when the offsets are very far and most low wavenumbers are considered. For P–P scattering the normalized wavenumber kz ∈ [0, 2ϰ cos (θmin − eff − P)], where $$\cos (\theta _{{\rm min}-{\rm eff}-P}) = \frac{K_z v_p}{2\omega _{{\rm min}}}<1$$. Third: the most practical and the worst in terms of the data coverage case is when the illumination of a wavenumber is limited on both ends. This inevitably happens for very low wavenumbers when the aperture is limited, that is, there are no diving waves traveling at the depth of perturbation. For this case, there are two limiting variables and for analysis of the plots we assume that the azimuth is either full or single. 4.2 Matrix A—SVD, Hessian To perform a quantitative analysis of coupling between orthorhombic parameters at a certain resolution determined by vertical wavenumber, we compute the SVD of matrix $$\mathbb {A}$$:   \begin{equation*} \mathbb {A}= \mathbb {U}\mathbb {S}\mathbb {V}^T \end{equation*} Matrix $$\mathbb {V}$$ is a 10 × 10 matrix which consists of normalized eigenvectors of $$\mathbb {A}^T\mathbb {A}$$. If the standard quadratic misfit J of the observed and modeled data (e.g. Tarantola 1984) is optimized by inversion,   \begin{eqnarray} J = ||\delta \mathbf {U}- \mathbb {A}\delta \mathbf {m}||^2 \end{eqnarray} (27)then the Hessian of the fitting problem is $$\mathbb {A}^T \mathbb {A}$$. The columns of $$\mathbb {V}^T$$ are the eigenvectors of the Hessian. Diagonal elements of the matrix $$\mathbb {S}$$ are the eigenvalues of the Hessian. Thus, $$\mathbb {V}$$ tells us which linear combinations of parameters are invertible from the data, and $$\mathbb {S}$$ tells the difference between the hardships to invert for these linear combinations. One of the most important features of the eigenvectors of the Hessian is that they are orthogonal with respect to the misfit functional norm. Thus, if we choose the parametrization in which each singular vector has only one non-zero component, we could invert all of the parameters one by one without any cross-talk issues. The problem of such a parametrization is that it would be spatially dependent as the set of singular vectors depends on the illumination of a particular domain. Therefore, there is no perfect parametrization for all of the cases and it can only be tuned according to the target depth and the acquisition parameters. Matrix $$\mathbb {S}$$ can have less than 10 non-zero values on the diagonal and, in that case, some trade-offs between parameters are irresolvable. Ties between the parameters need to be introduced (i.e. anisotropy complexity reduction or some substantial parameters, such as Gardner’s law) or additional data must be taken into inversion, for example, other wave types, etc. Even when $$\mathbb {S}$$ has 10 non-zero values, the number of invertible parameters might be reduced by the presence of noise in the data. The noise is accounted for in our inversion resolution study by using the model resolution matrices (Menke 2012; Podgornova et al.2015) with different thresholds for the singular values. The resolution matrix is calculated using the formula   $\mathbb {R}= \mathbb {V}^T (\mathbb {S}\mathbb {S}^T)^{-1} \mathbb {V},$ where $$(\mathbb {S}\mathbb {S}^T)^{-1}$$ is the generalized inverse of the matrix with threshold for singular values, which is determined by the amplitude of noise in the data. In the next section, we apply the apparatus described above to the sensitivity patterns. 5 RESULTS In this section, we quantitatively explore the following questions: What is the null space of linearized dynamic inversions? How many orthorhombic parameters can be resolved in realistic cases? What are these parameters? Of course, answers to these questions depend on the acquisition and data signal-to-noise ratio. Generally, the resolution quantification can be applied this way: Choose Kz of interest at a target depth. From bandwidth and aperture, estimate the range of effective angle of illumination (Fig. 6b). Define the set of effective incidence vectors s (θ, ϕ → s) using the azimuth coverage and range effective angle of illumination. Assemble matrix $$\mathbb {A}$$, which maps parameters into radiation patterns. Compute the SVD of matrix $$\mathbb {A}$$. Interpret the SVD and compute the resolution matrices. We focus on low wavenumbers recovery and consider three options for effective angle limitations: Perfect illumination—there are diving waves up to the level of perturbation and there are frequencies low enough to illuminate the wavenumber by backscattering. The null space in case of perfect illumination defines principally irresolvable trade-offs in between parameters. Lack of long offsets—there are low frequencies, but only reflections illuminate the target zone. We first consider the lack of large offsets in the data, presuming that frequencies, low enough to cover the chosen wavenumber by backscattering, are present in the data. This means Kz > ωmin/vp in the case of P–P scattering. The maximum offset defines the maximum opening angle available for the wavenumber illumination, and we investigate the decoupling at small scattering angles part of the reflection-based patterns. This is also the case for high wavenumber recovery in general. Lack of low frequencies—there are diving waves, but there are not enough low frequencies. Next, we assume that there are not enough low frequencies in the data and, therefore, $$\cos \frac{\theta _{{\rm max}}}{2} \omega _{{\rm min}}/v_p < K_z < \omega _{{\rm min}}/v_p$$. In this case, we also assume that there are diving waves that illuminate the perturbation, that is, θmax = π. For each of the three cases specified above, we consider arbitrary azimuth coverage. Instead of comparing the spectral sensitivity patterns visually, we analyse them quantitatively using the SVD analysis. Workflow for the perfect illumination case: Do the SVD of the $$\mathbb {A}$$ matrix (each row represents a particular scattering experiment). Data at all angles in radiation patterns are inverted together. The number of parameters that are invertible is equal to the number of non-zero singular values in SVD. The principal null space is the span of singular vectors related to zero singular values. The graph of the singular values tells us how many singular values are above a certain threshold, which defines the number of invertible parameters in the presence of noise. Noise defines the threshold for realistic cases. The threshold depends on the signal-to-noise ratio in the gradient of a misfit functional (e.g. Artman & Witten 2010), not in the data. We set the threshold for the resolution matrix at 0.1 of the highest singular value (this corresponds to the eigenvalues of the Hessian above 0.01) and create the resolution matrix. The elements with high values at the diagonal of the resolution matrix are candidates for the possibly invertible parameters. Finally, we show the resolution matrix for inversion in the pseudo-acoustic approximation. 5.1 Perfect illumination 5.1.1 Cij, ρ parametrization We analyse SVD of the full matrix $$\mathbb {A}_{P-P}$$ in the Cij, ρ parametrization. Fig. 7(a) shows that the matrix has only six non-zero singular values and thus only six independent parameters are invertible in this case. Fig. 7(b) shows the set of singular vectors in the Cij, ρ parametrization. The numeric null space is the linear span of four singular vectors, corresponding to zero singular values (columns 7–10 in Fig. 7b). These vectors are:   \begin{eqnarray} 2\mathbf {C}_{12}-\mathbf {C}_{66}, \end{eqnarray} (28)  \begin{eqnarray} 2\mathbf {C}_{13}+\mathbf {C}_{55}, \end{eqnarray} (29)  \begin{eqnarray} 2\mathbf {C}_{23}+\mathbf {C}_{44}, \end{eqnarray} (30)  \begin{eqnarray} \mathbf {C}_{11}+2\mathbf {C}_{12}+\mathbf {C}_{22}-\mathbf {C}_{33}-\pmb {\rho }. \end{eqnarray} (31)Each Cij represents a unit perturbation vector of the respective stiffness parameters, similarly $$\pmb {\rho }$$ is unit density perturbation. In Appendix A, these combinations, which we found empirically, are substituted into the reflection-based radiation pattern expressions showing that inversion is absolutely not sensitive to them. We also find a particular combination of the stiffness matrix elements and density that does not induce reflections. Thus, density can be effectively represented by anisotropic parameters within the first-order scattering approximation. Identifiable combinations of parameters can be divided into two groups—the first three singular vectors show slightly higher sensitivity than the other three, which is in agreement with de Hoop et al. (1999). The most interesting fact is that, despite previously published works (de Hoop et al.1999; Oh & Alkhalifah 2016), in the case of perfect illumination, the gap in this standard parametrization is within one order of magnitude and all six parameters can principally be determined, granted good quality data. Finally, we set the threshold for the resolution matrix to 0.1 of the highest singular value (this corresponds to the eigenvalues of the Hessian above 0.01) and create the resolution matrix for this case. The elements with the highest values on the diagonal are candidates for a possible inversion. Off-diagonal elements indicate cross-talk between the parameters. Fig. 7(c) shows the highest values for the Cii elements and the density, which are interesting to compare with pseudo-acoustic approximation, which suggests that C44, C55 and C66 have little influence on P-wave scattering. Figure 7. View largeDownload slide (a) Singular values and (b) vectors of the matrix $$\mathbb {A}$$ for P–P scattering and Cij, ρ parametrization. Perfect illumination—maximum opening angle is equal to π, all azimuths are available. There are six non-zero singular values. Later on for resolution matrices two thresholds are used (0,3—black and 0,1—green), corresponding to noisy and clean data respectively. (c) Resolution matrix in Cij, ρ parametrization. Figure 7. View largeDownload slide (a) Singular values and (b) vectors of the matrix $$\mathbb {A}$$ for P–P scattering and Cij, ρ parametrization. Perfect illumination—maximum opening angle is equal to π, all azimuths are available. There are six non-zero singular values. Later on for resolution matrices two thresholds are used (0,3—black and 0,1—green), corresponding to noisy and clean data respectively. (c) Resolution matrix in Cij, ρ parametrization. 5.1.2 Hierarchical parametrization The radiation patterns are significantly simplified in the hierarchical parametrization and more features are distinguishable. The hierarchical parametrization is a modification to that given by Tsvankin (1997), which was shown to primary control the traveltime, so it would be natural to expect inversion to be focused on them. However, according to the SVD (Fig. 8), there are other parameters involved and the null space is different from that of the traveltime inversion, even though it has the same dimensionality. If the density is dropped out, we see (Fig. 8c) that there are four evidently distinguishable elements of the matrix; those are Vp, Vs, ε1 and εd. These four parameters can be inverted. With the simplification of inversion comes the price of losing two vectors that cannot be inverted with the assumed signal-to-noise ratio. Overall, we conclude that the inversion with perfect illumination is possible for six parameters, unless the parametrization is chosen so that the inversion is more sensitive to a subset of anisotropic parameters (e.g. Oh & Alkhalifah 2016) for realistic acquisition scenarios. Figure 8. View largeDownload slide Same as Fig. 7, but for the hierarchical parametrization. Figure 8. View largeDownload slide Same as Fig. 7, but for the hierarchical parametrization. In the next subsection, we try to eliminate the problem of null space by restricting updates to the pseudo-acoustic parameters set. Pseudo-acoustic approximation:   The inversion is sensitive to density perturbations (Fig. 7b). So, let us fix density as well as C44, C55 and C66, removing them from the inversion in order to analyse whether we can still invert for six parameters with a parameter set reduced in this way. To establish some connection with the acoustic case, we remove the horizontal shear velocity and its deviation parameters from our analysis by setting their updates to zero. This approach to pseudo-acoustic modeling is well known to produce shear wave artefacts during propagation (Duveneck et al.2008), but good enough to compare P–P scattering, since only P-wave plane waves are used for the calculations. The P–P scattering happens for perturbations of several parameters in the new parametrization, including Vs, but we deliberately keep the shear wave-related parameters (Vs, γ1, γd) and density fixed in the inversion. In the pseudo-acoustic approximation, there are still six independently determinable parameters. In the standard parametrization, we see two levels of singular values very clearly (de Hoop et al.1999), yet the levels are within one order of magnitude. All six independently determinable parameters have approximately the same impact on the data and, thus, all of them can be determined in well-illuminated areas (Fig. 9c). If the image is noisy, then only parameters corresponding to the first three singular vectors can be retrieved (Fig. 9). Figure 9. View largeDownload slide (a) Singular values and (b) vectors of the matrix $$\mathbb {A}$$ for P–P scattering. Perfect illumination—maximum opening angle is equal to π, all azimuths are available. There are six non-zero singular values. Simplified pseudo-acoustic approximation (C44, C55, C66 and ρ are not updated). This approximation creates S-wave artefacts (e.g. Duveneck et al.2008), yet allows to compare the results in our case, since we keep the plane wave in the background fixed. Figure 9. View largeDownload slide (a) Singular values and (b) vectors of the matrix $$\mathbb {A}$$ for P–P scattering. Perfect illumination—maximum opening angle is equal to π, all azimuths are available. There are six non-zero singular values. Simplified pseudo-acoustic approximation (C44, C55, C66 and ρ are not updated). This approximation creates S-wave artefacts (e.g. Duveneck et al.2008), yet allows to compare the results in our case, since we keep the plane wave in the background fixed. In the case of hierarchical parametrization, the three parameters that define an ellipsoidal medium (vp, ε1, εd) should be inverted first (Fig. 10c). The other parameters are not determinable when there is noise in the data, even if the noise level is low. When limited illumination is considered, these resolution estimates become even worse and only vp is invertible. Figure 10. View largeDownload slide Same as Fig. 9, but for the hierarchical parametrization. Figure 10. View largeDownload slide Same as Fig. 9, but for the hierarchical parametrization. Figure 11. View largeDownload slide (a) and (c) Number of resolvable parameters for different maximum opening angles and azimuth coverage in Cij, ρ parametrization, and (b) and (d)  in the new parametrization. We assume that the image-domain signal-to-noise ratio is high—threshold for singular values is set to 0.1 of the maximum singular value which is somewhat equivalent to the image-domain signal-to-noise ratio of about 100 for (a) and (b). The threshold for singular values is set to 0.3 of the maximum singular value, corresponding to image-domain signal-to-noise ratio of about 10 (c) and (d). Figure 11. View largeDownload slide (a) and (c) Number of resolvable parameters for different maximum opening angles and azimuth coverage in Cij, ρ parametrization, and (b) and (d)  in the new parametrization. We assume that the image-domain signal-to-noise ratio is high—threshold for singular values is set to 0.1 of the maximum singular value which is somewhat equivalent to the image-domain signal-to-noise ratio of about 100 for (a) and (b). The threshold for singular values is set to 0.3 of the maximum singular value, corresponding to image-domain signal-to-noise ratio of about 10 (c) and (d). 5.2 Lack of long offsets We first consider the lack of large offsets in the data, presuming that frequencies low enough to cover the wavenumber by backscattering are present. This means Kz > ωmin/vp in the case of P–P scattering. Maximum offset defines the maximum opening angle available for the wavenumber illumination in this case. Therefore, we investigate decoupling in the little scattering angles part of the reflection-based patterns. Figs 11(a) and (b) show the number of parameters that could be resolved in the Cij, ρ and the hierarchical parametrization, respectively when the signal-to-noise ratio in the image is high. Figs 11(c) and (d) show the same but when the signal-to-noise ratio is low. The azimuth is calculated from the first horizontal axis, the direction does not matter as the orthorhombic media are assumed to be symmetric with respect to x1, x3 plane. Fig. 12 shows the same as Fig. 11, but with pseudo-acoustic approximation enforced. Figure 12. View largeDownload slide Same as Fig. 11, but with pseudo-acoustic approximation enforced. Figure 12. View largeDownload slide Same as Fig. 11, but with pseudo-acoustic approximation enforced. We deduce the following observations from Figs 11 and 12: When azimuth is 0–10 deg only VTI parameters are resolved. The maximum number of those resolved is equal to three and is in agreement with (Podgornova et al.2015). At large depths when maximum opening angle is below approximately 40 deg only one parameter can be resolved. In order to resolve four parameters in the new parametrization, one needs 45 deg azimuth coverage and about 140 deg maximum opening angle and relatively clean data. More parameters can be resolved using the standard parametrization in general. When azimuth coverage is increased beyond 90 deg it does not provide new information because of the symmetry of the orthorhombic media with respect to x2, x3 plane. Basically, one can forget about multiparameter inversion when the signal-to-noise ratio is very low. When the illumination is uneven, azimuth wise, the condition number sometimes drops even though the illumination range increases. 5.3 Lack of low frequencies Here, we assume that the offset is sufficient to illuminate the target depth with diving waves, but low frequencies are either not recorded at all or contaminated by severe noise and therefore not useful for FWI. In this case, low spatial wavenumbers in the model can only be reconstructed from wide opening angles (Sirgue & Pratt 2004; Kazei et al.2013a; Alkhalifah & Plessix 2014). The minimum effective opening angle determines the active part of reflection-based radiation patterns and therefore the capability of inversion to decouple parameters for a given spatial wavenumber. Just like for the case with the lack of long offsets, we cut the singular values at two thresholds 0.1 and 0.3 (Figs 13 and 14) to represent relatively clean and noisy data acquired respectively. We observe the following effects from Figs 13 and 14: Starting from 60 deg, the number of invertible parameters (Figs 13a and b) stales with angle range increase if the data is relatively clean. From transmissions alone (180 deg scattering data), we can determine up to three parameters in the classic parametrization, and two parameters (Vp, εd) in the new hierarchical parametrization. This is, likely, due to the fact that ηd parameter, which principally is also determinable, has very low sensitivity. In case of low signal-to-noise ratio in the data, we observe some counterintuitive effects. With increase in azimuth in observations the number of determinable parameters can decrease from three to two and from four to three (Fig. 13c), and even from two to one (Fig. 13d). This happens if the problem is solved in the most straightforward way. For example, if one applies the singular values thresholding based on condition number limitations for reduced problems. The maximum number of parameters is invertible when azimuths are limited to 90 and 180 deg, similarly to the previous section. Essentially, it means that including more data into a data misfit under optimization in inverse problem can lead to less parameters inverted instead of more if the illumination is uneven. Generally, the number of invertible parameters in our simplified setup is increased when pseudo-acoustic approximation is used for Cij, ρ parametrization, and decreased for hierarchical parametrization. The latter probably happens due to deliberate exclusion of Vs from inverted parameters in the hierarchical parametrization, which balances the parametrization towards Vp sensitivity. For the Cij, ρ parametrization, the first singular vector has significant part outside of the pseudo-acoustic space (dominated by density) and therefore singular values become more balanced with density exclusion from inversion. Figure 13. View largeDownload slide The number of invertible parameters when low frequencies lack from the data, yet the whole target area is illuminated by diving waves. (a) and (c) In the standard parametrization, and (b) and (d) in the new parametrization. High S/N ratio is assumed in (a) and (b) pictures, that is, singular values below 0.1 of the highest are low-cut. Low signal-to-noise ratio is assumed in (c) and (d). Figure 13. View largeDownload slide The number of invertible parameters when low frequencies lack from the data, yet the whole target area is illuminated by diving waves. (a) and (c) In the standard parametrization, and (b) and (d) in the new parametrization. High S/N ratio is assumed in (a) and (b) pictures, that is, singular values below 0.1 of the highest are low-cut. Low signal-to-noise ratio is assumed in (c) and (d). Figure 14. View largeDownload slide Same as Fig. 13, but with pseudo-acoustic approximation enforced. Figure 14. View largeDownload slide Same as Fig. 13, but with pseudo-acoustic approximation enforced. 6 DISCUSSION According to the effective angle of illumination concept, several frequencies are needed in the misfit functional in order to resolve the trade-offs between parameters. In monoparameter inversion, the lack of low frequencies can to some extent be compensated by long offsets (Sirgue & Pratt 2004; Kazei et al.2013b), and vice versa. In multiparameter inversion, the increase of the offset will lead to the blending of parameters. Thus, there is no point to increasing the offset without also increasing the bandwidth. In 3-D linearized inversion, several orthorhombic parameters can be recovered from a single frequency, which can be opposed to the bands (sets) of frequencies needed to perform multiparameter FWI in 2-D. However, the inversion of a single frequency can only be monoparameter, unless: There are several components recorded; then the number of parameters for inversion can be up to the number of components recorded. The observations are 3-D; then the information about azimuthal variations of scattering and corresponding parameters can be recovered. There are constraints on the parameter relations (e.g. constant Poisson ratio); then the multiparameter inversion can become monoparametric. Assumptions about the spectrum (structure) of the model lead to coverage of a single wavenumber from several angles at a given frequency, for example, different sparsity regularizations are implemented. Our study shows that the optimal parametrization is very much dependent on the case; moreover, it is dependent on the background model. The newly introduced parametrization shows high sensitivity to the longitudinal waves velocity, comparatively to the anisotropic coefficients. Thus, its primary use could be in building a velocity model when the anisotropy parameters are of secondary importance or in helping to invert for velocity. At the same time, if anisotropic parameters are the target of inversion, then it would be possible to fix the velocity and focus on the anisotropic parameters. For simultaneous inversion, it might be more reasonable to use a straightforward stiffness, density parametrization. In the standard parametrization, it makes no sense to invert for a subset of parameters because there are no significant gaps in between and the inversion would be unstable (Appendix C). For simplicity and generality, we consider Born-linearized inversion in an isotropic background and plane incident/scattered waves. Ichinose et al. (2000) showed that the far-field term is good enough for scatterers at distances about half of the wavelength from the source and receiver locations, which validates our analysis for marine streamer surveys of the reservoir areas. He & Plessix (2017) demonstrated that, in weak VTI backgrounds, the radiation patterns remain very similar to those in isotropic backgrounds. Therefore, we expect that the results of our analysis would not change dramatically if a background with weak anisotropy was considered. Here, we consider weak anisotropy approximation and, thus, our results compliment the numerical tests for strong background anisotropy by Köhn et al. (2015). Estimation of image-domain signal-to-noise ratio is an important tool to apply these general results, as data-domain signal-to-noise ratio cannot be used to quantify the resolution without particular acquisition set up. In order to compare parametrizations we use rescaling, which is equivalent setting vp and ρ to one as well as to that used in Podgornova et al. (2015) and Oh & Alkhalifah (2016). Our approach also showed that the condition number of the sensitivity matrix $$\mathbb {A}^T\mathbb {A}$$ heavily depends on the scaling of the parameters. In particular, normalizing partial derivatives admits SVD of $$\mathbb {A}$$ properties similar to those of Cij, ρ parametrization. It shows that even scalar scaling of parameters can change multiparameter inversion dramatically. 7 CONCLUSIONS Linearized dynamic inversion of P waves for weakly anisotropic orthorhombic medium is fundamentally sensitive to six parameters only. Surprisingly, the number of parameters is the same as for the traveltime inversion. In particular, we found that vertical variations of density show the same scattering as a certain combination of elastic constants. Pseudo-acoustic inversion (i.e. inversion for pseudo-acoustic parameters only, keeping other elastic parameters fixed) constraint allows to completely remove the fundamental ambiguity in linearized inversion of P waves for orthorhombic media in the case of perfect illumination. The coverage of different wavenumbers in the velocity is investigated through spectral sensitivities. We introduced the effective angle of illumination concept, which allowed us to investigate illumination for particular wavenumbers in the velocity model. We estimated the number of invertible anisotropic orthorhombic parameters, given the acquisition setup, by using SVD of the spectral sensitivities. It turns out that the azimuth increase beyond 60 deg does not bring significant additional information on the resolvable parameters set when inverting for orthorhombic media. For monoparameter FWI, the lack of low frequencies can be compensated by long offsets in the data (Sirgue & Pratt 2004; Kazei et al.2013a), for multiparameter inversion both need to be present to allow for parameter decoupling at low wavenumbers. In order to decouple VTI parameters reliably, the maximum frequency included into the inversion needs to be at least twice higher than the minimum one. The best resolution is provided when the target area is covered by diving waves, yet even if the illumination is limited to 120 deg opening angles full set of parameters can be recovered in the classic parametrization. The azimuthally varying parameters are more or less well resolved starting from π/3 azimuth coverage. While the hierarchical parametrization reduces the trade-offs between parameters (Oh & Alkhalifah 2016), it also reduces the number of parameters that can be inverted focusing inversion on a subset of four parameters (Vp, Vs, ε1, εd) in case elastic inversion is performed and on a subset of three parameters (Vp, ε1, εd), which is covering the case of ellipsoidal anisotropy or prior knowledge anellipticity parameters (η1, ηd). If the goal of the inversion is to obtain a velocity model that explains most of the observed data then the hierarchical parametrization is optimal. If we need to extract more information about the anisotropy of a medium, considering optimal conditions and utilization of a Hessian approximation capable of treating the trade-off, and with a good initial model (such as that constructed using the hierarchical parametrization), then the classic Cij, ρ parametrization can help recover more anisotropic features of a reservoir at later stages. Acknowledgements We thank Alexey Stovas and Daniel Köhn for their constructive and helpful reviews of the manuscript. We are grateful to Ju-Won Oh, Nabil Masmoudi and Ramzi Djebbi of SWAG, KAUST for helpful discussions as well. VK is also grateful to Oleg Ovcharenko and Veronica Tremblay of KAUST for their suggestions on improving the manuscript readability. Finally, we thank KAUST for support and its HPC resources. GLOSSARY ω — angular frequency. x — imaging point vector. xs, xg — source and geophone radius vectors, respectively. s — unit vector pointing towards the source along the source wavepath direction. g — same as s but towards the geophone. Gik(x, x΄) ≡ Gik(x, x΄, ω) — Green’s tensor in the background medium for frequency ω. K — wavenumber domain vector. u(xs, xg) — monochromatic displacement wavefield at frequency ω. cijkl, c — stiffness tensor. Cij, C — stiffness matrix in Voigt notation. ρ, λ — density and first lame parameter. δf — variational perturbation of a function f. $$\hat{f}$$ — Fourier transform of a function f. WI − WS, δm — incident and scattered wave types and the perturbation parameter type, for example WI − WS, δm = P-P, λ or WI − WS, δm = P-SH, C13. $$\mathcal {R}_{WI-WS, \delta \mathbf {m}}(\mathbf {s},\mathbf {g})$$ — radiation pattern for plane waves and parameter perturbation of type WI − WS, δm. Einstein summation is assumed throughout the paper. (sg)ij ≡ sigj – the outer (dyadic) product of tensors (vectors). : — convolution of tensors on two indexes, for example c: sg ≡ cijklskgl. We choose Fourier convention following (Hudson & Heritage 1981):   \begin{eqnarray} u(t) = \int \limits _{-\infty }^{+\infty }\exp (-i\omega t) U(\omega ) {\rm d}\omega ,\, \hat{U}(\mathbf {K}) = \int \limits _{V}\exp (-i\mathbf {K}\mathbf {x}) U (\mathbf {x}) {\rm d}\mathbf {x}.\nonumber\\ \end{eqnarray} (32) REFERENCES Alkhalifah T., 2000. An acoustic wave equation for anisotropic media, Geophysics , 65( 4), 1239– 1250. https://doi.org/10.1190/1.1444815 Google Scholar CrossRef Search ADS   Alkhalifah T., Plessix R.-E., 2014. A recipe for practical full-waveform inversion in anisotropic media: an analytical parameter resolution study, Geophysics , 79( 3), R91– R101. https://doi.org/10.1190/geo2013-0366.1 Google Scholar CrossRef Search ADS   Alkhalifah T., 2016. Full-model wavenumber inversion: An emphasis on the appropriate wavenumber continuation, Geophysics , 81( 3), R89– R98. https://doi.org/10.1190/geo2015-0537.1 Google Scholar CrossRef Search ADS   Anikiev D., Kazei V., Kashtan B., Ponomarenko A., Troyan V., Shigapov R., 2014. Methods of seismic waveform inversion, Seism. Technol. , 11( 1), 1– 32. Artman B., Witten B., 2010. Image domain signal to noise estimate , US Patent App. 13/145,328. Bergslid T.S., Raknes E.B., Arntsen B., 2015. The influence of anisotropy on elastic full-waveform inversion, in SEG Technical Program Expanded Abstracts 2015 , pp. 1425– 1429, Society of Exploration Geophysicists. Beylkin G., Burridge R., 1990. Linearized inverse scattering problems in acoustics and elasticity, Wave Motion , 12( 1), 15– 52. https://doi.org/10.1016/0165-2125(90)90017-X Google Scholar CrossRef Search ADS   Cheverda V., Kostin V., 1995. R-pseudoinverses for compact operators in hilbert spaces: existence and stability, J. Inverse Ill-Posed Probl. , 3( 2), 131– 148. https://doi.org/10.1515/jiip.1995.3.2.131 Google Scholar CrossRef Search ADS   de Hoop M.V., Spencer C., Burridge R., 1999. The resolving power of seismic amplitude data: an anisotropic inversion/migration approach, Geophysics , 64( 3), 852– 873. https://doi.org/10.1190/1.1444595 Google Scholar CrossRef Search ADS   Devaney A.J., 1984. Geophysical diffraction tomography, IEEE Trans. Geosci. Remote Sens. , GE-22( 1), 3– 13. https://doi.org/10.1109/TGRS.1984.350573 Google Scholar CrossRef Search ADS   Duveneck E., Milcik P., Bakker P.M., Perkins C., 2008. Acoustic VTI wave equations and their application for anisotropic reverse-time migration, in SEG Technical Program Expanded Abstracts 2008 , pp. 2186– 2190, Society of Exploration Geophysicists. Eaton D.W.S., Stewart R.R., 1994. Migration/inversion for transversely isotropic elastic media, Geophys. J. Int. , 119( 2), 667– 683. https://doi.org/10.1111/j.1365-246X.1994.tb00148.x Google Scholar CrossRef Search ADS   Ewald P., 1969. Introduction to the dynamical theory of x-ray diffraction, Acta Crystallogr. A: Cryst. Phys. Diffract. Theor. Gen. Crystallogr. , 25( 1), 103– 108. https://doi.org/10.1107/S0567739469000155 Google Scholar CrossRef Search ADS   Gholami Y., Brossier R., Operto S., Prieux V., Ribodetti A., Virieux J., 2013a. Which parameterization is suitable for acoustic vertical transverse isotropic full waveform inversion? Part 2: synthetic and real data case studies from valhall, Geophysics , 78( 2), R107– R124. https://doi.org/10.1190/geo2012-0203.1 Google Scholar CrossRef Search ADS   Gholami Y., Brossier R., Operto S., Ribodetti A., Virieux J., 2013b. Which parameterization is suitable for acoustic vertical transverse isotropic full waveform inversion? Part 1: sensitivity and trade-off analysis, Geophysics , 78( 2), R81– R105. https://doi.org/10.1190/geo2012-0204.1 Google Scholar CrossRef Search ADS   Grechka V., Tsvankin I., 1999. 3-d moveout velocity analysis and parameter estimation for orthorhombic media, Geophysics , 64( 3), 820– 837. https://doi.org/10.1190/1.1444593 Google Scholar CrossRef Search ADS   He W., Plessix R.-É., 2017. Analysis of different parameterisations of waveform inversion of compressional body waves in an elastic transverse isotropic earth with a vertical axis of symmetry, Geophys. Prospect. , 65( 4), 1004– 1024. https://doi.org/10.1111/1365-2478.12452 Google Scholar CrossRef Search ADS   Hudson J., Heritage J., 1981. The use of the born approximation in seismic scattering problems, Geophys. J. Int. , 66( 1), 221– 240. https://doi.org/10.1111/j.1365-246X.1981.tb05954.x Google Scholar CrossRef Search ADS   Ichinose G.A., Goldstein P., Rodgers A.J., 2000. Relative importance of near-, intermediate-and far-field displacement terms in layered earth synthetic seismograms, Bull. seism. Soc. Am. , 90( 2), 531– 536. https://doi.org/10.1785/0119990134 Google Scholar CrossRef Search ADS   Ivanov Y., Stovas A., 2016. Upscaling in orthorhombic media: Behavior of elastic parameters in heterogeneous fractured earth, Geophysics , 81( 3), C113– C126. https://doi.org/10.1190/geo2015-0392.1 Google Scholar CrossRef Search ADS   Kamath N., Tsvankin I., 2016. Elastic full-waveform inversion for vti media: Methodology and sensitivity analysis, Geophysics , 81( 2), C53– C68. https://doi.org/10.1190/geo2014-0586.1 Google Scholar CrossRef Search ADS   Kamei R., Pratt R.G., Tsuji T., 2013. On acoustic waveform tomography of wide-angle obs data?strategies for pre-conditioning and inversion, Geophys. J. Int. , 194( 2), 1250– 1280. https://doi.org/10.1093/gji/ggt165 Google Scholar CrossRef Search ADS   Kazei V.V., Alkhalifah T., 2017. On the Resolution of Inversion for Orthorhombic Anisotropy, in 79th EAGE Conference and Exhibition 2017 . Kazei V.V., Troyan V.N., Kashtan B.M., Mulder W.A., 2013a. On the role of reflections, refractions and diving waves in full-waveform inversion, Geophys. Prospect. , 61( 6), 1252– 1263. https://doi.org/10.1111/1365-2478.12064 Google Scholar CrossRef Search ADS   Kazei V.V., Kashtan B.M., Troyan V.N., Mulder W.A., 2013b. Spectral sensitivity analysis of FWI in a constant-gradient background velocity model, in 75th EAGE Conference and Exhibition Incorporating SPE EUROPEC 2013 . Kazei V.V., Kashtan B.M., Troyan V.N., Mulder W.A., 2015. FWI spectral sensitivity analysis in the presence of a free surface, in SEG Technical Program Expanded Abstracts 2015 , pp. 1415– 1419, Society of Exploration Geophysicists. Kazei V., Tessmer E., Alkhalifah T., 2016. Scattering angle-based filtering via extension in velocity, in SEG Technical Program Expanded Abstracts 2016 , pp. 1157– 1162, Society of Exploration Geophysicists. Köhn D., Hellwig O., De Nil D., Rabbel W., 2015. Waveform inversion in triclinic anisotropic media-a resolution study, Geophys. J. Int. , 201( 3), 1642– 1656. https://doi.org/10.1093/gji/ggv097 Google Scholar CrossRef Search ADS   Kurzmann A. et al.  , 2016. Seismic Applications of Full Waveform Inversion, in High Performance Computing in Science and Engineering'16 , pp. 647– 665, eds Nagel W., Kröner D., Resch M., Springer, Cham. Google Scholar CrossRef Search ADS   Masmoudi N., Alkhalifah T., 2016. A new parameterization for waveform inversion in acoustic orthorhombic media, Geophysics , 81( 4), R157– R171. https://doi.org/10.1190/geo2015-0635.1 Google Scholar CrossRef Search ADS   Menke W., 2012. Geophysical Data Analysis: Discrete Inverse Theory , Vol. 45, Academic Press. Mora P., 1989. Inversion = migration + tomography, Geophysics , 54( 12), 1575– 1586. https://doi.org/10.1190/1.1442625 Google Scholar CrossRef Search ADS   Oh J.-W., Alkhalifah T., 2016. Elastic orthorhombic anisotropic parameter inversion: an analysis of parameterization, Geophysics , 81( 6), C279– C293. https://doi.org/10.1190/geo2015-0656.1 Google Scholar CrossRef Search ADS   Ovcharenko O., Kazei V., Peter D., Alkhalifah T., 2018. Variance-based salt body reconstruction for improved full-waveform inversion, Geophysics , submitted. Owusu J.C. et al.  , 2016. Anisotropic elastic full-waveform inversion of walkaway vertical seismic profiling data from the arabian gulf, Geophys. Prospect. , 64( 1), 38– 53. https://doi.org/10.1111/1365-2478.12227 Google Scholar CrossRef Search ADS   Podgornova O., Leaney S., Liang L., 2015. Analysis of resolution limits of vti anisotropy with full waveform inversion, in SEG Technical Program Expanded Abstracts 2015 , pp. 1188– 1192, Society of Exploration Geophysicists. Raknes E.B., Arntsen B., 2014. Time-lapse full-waveform inversion of limited-offset seismic data using a local migration regularization, Geophysics , 79( 3), WA117– WA128. https://doi.org/10.1190/geo2013-0369.1 Google Scholar CrossRef Search ADS   Rüger A., 1997. P-wave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry, Geophysics , 62( 3), 713– 722. https://doi.org/10.1190/1.1444181 Google Scholar CrossRef Search ADS   Sacks P.E., 1988. The inverse problem for a weakly inhomogeneous two-dimensional acoustic medium, SIAM J. Appl. Math. , 48( 5), 1167– 1193. https://doi.org/10.1137/0148070 Google Scholar CrossRef Search ADS   Santosa F., Symes W.W., 1989. Analysis of Least-Squares Velocity Inversion , Society of Exploration Geophysicists. Google Scholar CrossRef Search ADS   Schoenberg M., Helbig K., 1997. Orthorhombic media: modeling elastic wave behavior in a vertically fractured earth, Geophysics , 62( 6), 1954– 1974. https://doi.org/10.1190/1.1444297 Google Scholar CrossRef Search ADS   Shaw R.K., Sen M.K., 2004. Born integral, stationary phase and linearized reflection coefficients in weak anisotropic media, Geophys. J. Int. , 158( 1), 225– 238. https://doi.org/10.1111/j.1365-246X.2004.02283.x Google Scholar CrossRef Search ADS   Sirgue L., Pratt R., 2004. Efficient waveform inversion and imaging: a strategy for selecting temporal frequencies, Geophysics , 69( 1), 231– 248. https://doi.org/10.1190/1.1649391 Google Scholar CrossRef Search ADS   Snieder R., 2002. Chapter 1.7.1—general theory of elastic wave scattering, in Scattering , pp. 528– 542, eds Pike R., Sabatier P., Academic Press, London. Stovas A., 2017. Kinematic parameters of pure- and converted-mode waves in elastic orthorhombic media, Geophys. Prospect. , 65( 2), 426– 452. https://doi.org/10.1111/1365-2478.12420 Google Scholar CrossRef Search ADS   Tarantola A., 1984. Inversion of seismic reflection data in the acoustic approximation, Geophysics , 49( 8), 1259– 1266. https://doi.org/10.1190/1.1441754 Google Scholar CrossRef Search ADS   Tarantola A., 1986. A strategy for nonlinear elastic inversion of seismic reflection data, Geophysics , 51( 10), 1893– 1903. https://doi.org/10.1190/1.1442046 Google Scholar CrossRef Search ADS   ten Kroode F., Bergler S., Corsten C., de Maag J.W., Strijbos F., Tijhof H., 2013. Broadband seismic data ? The importance of low frequencies, Geophysics , 78( 2), WA3– WA14. https://doi.org/10.1190/geo2012-0294.1 Google Scholar CrossRef Search ADS   Thomsen L., 1986. Weak elastic anisotropy, Geophysics , 51( 10), 1954– 1966. https://doi.org/10.1190/1.1442051 Google Scholar CrossRef Search ADS   Tikhonov A.N., Arsenin V.I., John F., 1977. Solutions of Ill-posed Problems , Vol. 14, Winston, Washington, DC. Tsvankin I., 1997. Anisotropic parameters and p-wave velocity for orthorhombic media, Geophysics , 62( 4), 1292– 1309. https://doi.org/10.1190/1.1444231 Google Scholar CrossRef Search ADS   Virieux J., Operto S., 2009. An overview of full-waveform inversion in exploration geophysics, Geophysics , 74( 6), WCC1– WCC26. https://doi.org/10.1190/1.3238367 Google Scholar CrossRef Search ADS   Wu R., Toksöz M., 1987. Diffraction tomography and multisource holography applied to seismic imaging, Geophysics , 52( 1), 11– 25. https://doi.org/10.1190/1.1442237 Google Scholar CrossRef Search ADS   APPENDIX A: NULL SPACE OF P–P INVERSION Here, we show that some radiation patterns coincide analytically in isotropic background when inversion for vertical wavenumbers is considered. Let us start with eq. (28).   \begin{eqnarray} \mathcal {R}_{2\mathbf {C}_{12}-\mathbf {C}_{66},P-P} \nonumber \\ &=& \mathcal {R}_{2\mathbf {c}_{1122}-\mathbf {c}_{1212}-\mathbf {c}_{2121},P-P}(\mathbf {s},\mathbf {g}) \nonumber \\ &=& 2 s_1s_1g_2g_2 - 2s_1s_2g_1g_2 = 0. \end{eqnarray} (A1)This happens due to the fact that for horizontal reflectors or pure vertical wavenumbers:   \begin{eqnarray} (s_1=-g_1,s_2=-g_2). \end{eqnarray} (A2)Analogously the linear combinations:   \begin{eqnarray} 2\mathbf {C}_{13}+\mathbf {C}_{55}, \end{eqnarray} (A3)  \begin{eqnarray} 2\mathbf {C}_{23}+\mathbf {C}_{44}, \end{eqnarray} (A4)do not scatter P–P waves. The other Cij combinations fall into the null space in very similar way using the connection s3 = g3. For density, the relation is a little bit trickier:   \begin{eqnarray*} \mathcal {R}_{\vec{\rho },P-P} &=& \mathbf {s}\cdot \mathbf {g}\nonumber \\ &=& s_1g_1 + s_2g_2 + s_3g_3 \nonumber \\ &=& -s_1s_1-s_2s_2+s_3s_3 \nonumber \\ &=& -s_1s_1(s_1s_1+s_2s_2+s_3s_3) \nonumber \\ &&-\, s_2s_2(s_1s_1+s_2s_2+s_3s_3) \nonumber \\ &&+\, s_3s_3(s_1s_1+s_2s_2+s_3s_3) \nonumber \\ &=& \mathcal {R}_{-\mathbf {C}_{11}-\mathbf {C}_{12}-\mathbf {C}_{13},P-P} \nonumber \\ &&+\, \mathcal {R}_{-\mathbf {C}_{21}-\mathbf {C}_{22}-\mathbf {C}_{23},P-P} \nonumber \\ &&+\, \mathcal {R}_{\mathbf {C}_{31}+\mathbf {C}_{32}+\mathbf {C}_{33},P-P} \nonumber \\ &=& \mathcal {R}_{-\mathbf {C}_{11}-2\mathbf {C}_{12}-\mathbf {C}_{22}+\mathbf {C}_{33},P-P}. \end{eqnarray*} This finishes the proof of the possibility to represent density through anisotropic parameters perturbations. The following linear combination does not scatter:   \begin{eqnarray} \mathbf {C}_{11}+2\mathbf {C}_{12}+\mathbf {C}_{22}-\mathbf {C}_{33}-\rho . \end{eqnarray} (A5) APPENDIX B: HIERARCHICAL PARAMETRIZATION Stiffness matrix coefficients can be expressed in the hierarchical parametrization as follows (Oh & Alkhalifah 2016):   \begin{eqnarray} C_{11} =\rho V^2_{p}, \end{eqnarray} (B1)  \begin{eqnarray} C_{22} =(1+2\varepsilon _d)\rho V^2_{p}, \end{eqnarray} (B2)  \begin{eqnarray} C_{33} = \frac{1}{1+2\varepsilon _1} \rho V^2_{p}, \end{eqnarray} (B3)  \begin{eqnarray} C_{12} = \rho \sqrt{ \left( V^2_{p}-(1+2\gamma _1) V^2_{s} \right)} \end{eqnarray} (B4)  \begin{eqnarray} \times \sqrt{\left((1+2\delta _3)V^2_{p}-(1+2\gamma _1) V^2_{s} \right)} \end{eqnarray} (B5)  \begin{eqnarray} - (1+2\gamma _1 )\rho V^2_{s}, \end{eqnarray} (B6)  \begin{eqnarray} C_{13} = \rho \sqrt{\left(\frac{V^2_{p}}{1+2\varepsilon _1}-V^2_{s} \right) \left(\frac{V^2_{p}}{1+2\eta _1}-V^2_{s} \right)} \end{eqnarray} (B7)  \begin{eqnarray} -\rho V^2_{s}, \end{eqnarray} (B8)  \begin{eqnarray} C_{23} = \rho \left[\left(\frac{V^2_{p}}{1+2\varepsilon _1}(1+2\gamma _d)v{s}^2 \right)\right. \end{eqnarray} (B9)  \begin{eqnarray} \times \left.\left(\frac{V^2_{p}(1+\varepsilon _d)}{(1+2\eta _1)(1+2\eta _d)}-(1+2\gamma _d)V^2_{s} \right)\right]^{1/2} \end{eqnarray} (B10)  \begin{eqnarray} -&(1+2\gamma _d)\rho V^2_{s}, \end{eqnarray} (B11)  \begin{eqnarray} C_{44} =(1+2\gamma _d)\rho V^2_{s}, \end{eqnarray} (B12)  \begin{eqnarray} C_{55} =\rho V^2_{s}, \end{eqnarray} (B13)  \begin{eqnarray} C_{66} =(1+2\gamma _1)\rho V^2_{s}. \end{eqnarray} (B14)Which results in the following expression for the Jacobian $$\frac{\partial \mathbf {C}}{\partial \mathbf {m}}$$, where ϰ = Vs/Vp:   \begin{eqnarray} \delta \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}\rho \\ C_{11} \\ C_{22} \\ C_{33} \\ C_{12} \\ C_{13} \\ C_{23} \\ C_{44} \\ C_{55} \\ C_{66} \end{array}\right) = \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}1&0&0&0&0&0&0&0&0&0\\ 1&2&0&0&0&0&0&0&0&0\\ 1&2&0&0&2&0&0&0&0&0\\ 1&2&0&-2&0&0&0&0&0&0\\ 1-2\,{\it \varkappa }&2&-4\,{\it \varkappa }&0&0&0&0&1&-4\,{\it \varkappa }&0\\ 1-2\,{\it \varkappa }&2&-4\,{\it \varkappa }&-1&0&-1&0&0&0&0\\ 1-2\,{\it \varkappa }&2&-4\,{\it \varkappa }&-1&1&-1&-1&0&0&-4\,{\it \varkappa }\\ {\it \varkappa }&0&2\,{\it \varkappa }&0&0&0&0&0&0&2\,{\it \varkappa }\\ {\it \varkappa }&0&2\,{\it \varkappa }&0&0&0&0&0&0&0\\ {\it \varkappa }&0&2\,{\it \varkappa }&0&0&0&0&0&2\,{\it \varkappa }&0\end{array}\right) \delta \left(\begin{array}{c@{\quad}c@{\quad}c}\rho \\ V_{p}\\ V_{s}\\ \varepsilon _1\\ \varepsilon _d\\ \eta _1\\ \eta _d\\ \delta _3\\ \gamma _1\\ \gamma _d \end{array}\right) \end{eqnarray} (B15) APPENDIX C: ON STABILITY OF INVERSION BY SVD TRUNCATION In exploration geophysics and in seismic inversions in particular, we often face the problem of data not completely covering the model parameters therefore prior information through regularization is required (Tikhonov et al.1977). One of the possibilities to regularize inversion is through chopping off the singular values for generalized inverse (Santosa & Symes 1989; Sacks 1988). The question that naturally arises is ‘how many singular values should be kept?’. If one keeps r singular values the inverse is called r-pseudo-inverse and the solution obtained this way is called the r-solution (Cheverda & Kostin 1995). Here, we will refer to particular results from Cheverda & Kostin (1995) on the stability of r-pseudo-inverses to justify the choice of the thresholds. This is especially true in case of multiparameter inversions. The linearized system of equation:   \begin{eqnarray} \mathbb {A}\mathbf {m}= \mathbf {d}, \end{eqnarray} (C1)where m is the parameter perturbation, d is the observed data, does not have unique solution as there are parameter perturbations that do not scatter. Regularizing the system—introducing the prior information can, however provide us with some solutions. In this case, we look at the singular value decomposition (SVD) based pseudo-inverse of the matrix $$\mathbb {A}^{-1}$$, that provides us with a solution of minimum length to eq. (C1). With singular value decomposition of matrix $$\mathbb {A}$$:   \begin{eqnarray} \mathbb {A}= \mathbb {U}\mathbb {S}\mathbb {V}^*, \end{eqnarray} (C2)its pseudo-inverse is   \begin{eqnarray} \mathbb {A}^{-1} = \mathbb {V}\mathbb {S}^{-1}\mathbb {U}^*. \end{eqnarray} (C3)The solution is defined as follows:   \begin{eqnarray} \mathbf {m}_{{\rm est}} = \mathbb {A}^{-1}\mathbf {d}. \end{eqnarray} (C4)The matrix of singular values cannot be inverted for underdetermined problems, therefore in $$\mathbb {S}^{-1}$$ there are several singular values that can be inverted. In numerical applications, some singular values would not be zero exactly and therefore it is a key point to understand which singular values need to be treated as zero. There are several strategies to decide which singular values to drop: Singular values that are lower than a certain fraction of the first singular value. This allows to have certain condition number for the pseudo-inverse, which is important for stability of the inversion. Singular values that come after a large enough gap. This ensures that the singular values/vectors would not mix in between themselves and therefore contributes to the stability estimates. Keep singular values that combine into certain part of the Hilbert–Schmidt norm of $$\mathbb {S}$$. The latter makes sense for images storages—certain part of the energy contained in the image itself is than captured and the rest is dropped, which does not work for the operators. If the data d and the operator $$\mathbb {A}$$ are perturbed by δd and $$\delta \mathbb {A}$$, respectively, and, relative perturbations are bounded by ε and δ, respectively:   \begin{eqnarray} \frac{||\delta \mathbb {A}||}{||\mathbb {A}||} < \delta , \, \frac{||\delta d||}{||d||} < \varepsilon . \end{eqnarray} (C5)If also the relative gap in the singular values si of operator $$\mathbb {A}$$ is sufficiently large   \begin{eqnarray} \frac{|s_1|}{|s_r| - |s_{r+1}|} = d_r < 1/2\delta , \end{eqnarray} (C6)then the r-pseudo-inverse $$(\mathbb {A}+\delta \mathbb {A})^{-1}$$ exists (Cheverda & Kostin 1995) for compact operators (not only matrices). Moreover Cheverda & Kostin (1995) show that the relative perturbation δmest of the r-solution is bounded in the following way:   \begin{eqnarray} \frac{||\delta \mathbf {m}_{{\rm est}}||}{||\mathbf {m}_{{\rm est}}||} \le \mu _r(\mathbb {A}) \frac{(\rho +\tau )(1+\beta \delta )}{1-\mu _r(\mathbb {A})\rho }, \,{\rm if}\, \mu _r \rho < 1. \end{eqnarray} (C7)The following quantities are used in eq. (C7): the inconsistency parameter θ   \begin{eqnarray} \theta = \frac{||\mathbb {A}\mathbf {m}_{{\rm est}} - \mathbf {d}||}{s_r ||\mathbf {m}_{{\rm est}}||}, \end{eqnarray} (C8)the condition number of the system:   \begin{eqnarray} \mu _r = \frac{s_1}{s_r}, \end{eqnarray} (C9)the other quantities in eq. (C7)   \begin{eqnarray} \beta = \frac{\sqrt{2}d_r}{\sqrt{1-d_r \delta }\sqrt{1-d_r\delta \sqrt{1-2d_r\delta }}}, \end{eqnarray} (C10)  \begin{eqnarray} \rho = \delta (1+2\beta ), \end{eqnarray} (C11)  \begin{eqnarray} \tau = \sqrt{1+\theta ^2} \left(\varepsilon + \frac{d_r\delta }{1-d_r\delta } + \beta \delta \right) \end{eqnarray} (C12)The Cheverda–Kostin’s theorem, that we formulated above, essentially proposes conditions that are sufficient to have stable inversion for r parameters. Here, we just try to understand what is the limitation on singular values caused by the requirement:   \begin{eqnarray} \mu _r \rho < 1. \end{eqnarray} (C13)The latter is equivalent to   \begin{eqnarray} \frac{s_1}{s_r} \delta (1+2\beta ) < 1, \end{eqnarray} (C14)which poses the allowable mistakes in the estimates of the operator $$\mathbb {A}$$. We will strengthen the requirements by substituting β with   \begin{eqnarray} \beta _{up} \equiv 2 \sqrt{2} d_r > \beta , \end{eqnarray} (C15)which leads to bound,   \begin{eqnarray} \delta < \frac{s_r}{s_1(1+ 4\sqrt{2} d_r) }, \end{eqnarray} (C16)therefore it is enough to have   \begin{eqnarray} \delta < \frac{1}{7\mu _r d_r}. \end{eqnarray} (C17)Eq. (C17) suggests that it is sufficient to start very close to the real model in case of multiparameter inversions, to be sure that there is no ambiguity in results of multiparameter inversion. Unfortunately, condition (C17) is hard to verify in reality. If we consider linearized inversion for a single parameter though μr = 1, the stability of inversion is governed by the gap in singular values only. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Waveform inversion for orthorhombic anisotropy with P waves: feasibility and resolution

, Volume 213 (2) – May 1, 2018
20 pages

/lp/ou_press/waveform-inversion-for-orthorhombic-anisotropy-with-p-waves-RkuymC3ubB
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy034
Publisher site
See Article on Publisher Site

### Abstract

Summary Various parametrizations have been suggested to simplify inversions of first arrivals, or P waves, in orthorhombic anisotropic media, but the number and type of retrievable parameters have not been decisively determined. We show that only six parameters can be retrieved from the dynamic linearized inversion of P waves. These parameters are different from the six parameters needed to describe the kinematics of P waves. Reflection-based radiation patterns from the P–P scattered waves are remapped into the spectral domain to allow for our resolution analysis based on the effective angle of illumination concept. Singular value decomposition of the spectral sensitivities from various azimuths, offset coverage scenarios and data bandwidths allows us to quantify the resolution of different parametrizations, taking into account the signal-to-noise ratio in a given experiment. According to our singular value analysis, when the primary goal of inversion is determining the velocity of the P waves, gradually adding anisotropy of lower orders (isotropic, vertically transversally isotropic and orthorhombic) in hierarchical parametrization is the best choice. Hierarchical parametrization reduces the trade-off between the parameters and makes gradual introduction of lower anisotropy orders straightforward. When all the anisotropic parameters affecting P-wave propagation need to be retrieved simultaneously, the classic parametrization of orthorhombic medium with elastic stiffness matrix coefficients and density is a better choice for inversion. We provide estimates of the number and set of parameters that can be retrieved from surface seismic data in different acquisition scenarios. To set up an inversion process, the singular values determine the number of parameters that can be inverted and the resolution matrices from the parametrizations can be used to ascertain the set of parameters that can be resolved. Fourier analysis, Inverse theory, Waveform inversion, Body waves, Seismic anisotropy, Wave scattering and diffraction 1 INTRODUCTION Reservoir delineation and monitoring are among the major objectives of modern seismic imaging and inversion. The elastic properties of most hydrocarbon reservoirs can be well approximated using orthorhombic anisotropy (Schoenberg & Helbig 1997; Tsvankin 1997; Ivanov & Stovas 2016) on the scale of common exploration seismic wavelengths. Systems of orthogonal fractures or vertical layering within a reservoir can be described by macroscale orthorhombic media; see Fig. 1(a). Seismic imaging of a reservoir area aims to obtain the effective orthorhombic parameters, which are then converted into the microscale properties of the imaged reservoir using rock-physics modeling. The accuracy of the orthorhombic parameter estimates is crucial to a reliable conversion. Figure 1. View largeDownload slide Orthorhombic reservoir. (a) Anisotropic reservoir (Tsvankin 1997) in isotropic background. (b) Elements of stiffness matrix in Voigt notation. If the media is orthorhombic then only the red cells are non-zero in its main axes, monoclinic with horizontal symmetry plane by red and blue cells. Figure 1. View largeDownload slide Orthorhombic reservoir. (a) Anisotropic reservoir (Tsvankin 1997) in isotropic background. (b) Elements of stiffness matrix in Voigt notation. If the media is orthorhombic then only the red cells are non-zero in its main axes, monoclinic with horizontal symmetry plane by red and blue cells. Conventionally, a moveout velocity analysis of the first arrivals, that is, P waves, is used to provide an estimate of the orthorhombic medium parameters (Grechka & Tsvankin 1999). P-wave velocity analysis (Tsvankin 1997; Grechka & Tsvankin 1999; Stovas 2017) defines the limited resolution of classic ray-based methods in orthorhombic media. Recently, interest has grown in the inversion of dynamic wavefield components due to their higher resolution potential, despite the additional acquisition and computational requirements (e.g. Virieux & Operto 2009; Anikiev et al.2014; Kurzmann et al.2016). Iterative imaging techniques, such as full-waveform inversion (FWI) and least-squares reverse-time migration, provide high-resolution estimates of the subsurface parameters, but require more computational resources and above-average signal-to-noise ratios in the data. Marine data sets are usually able to meet the signal-to-noise ratio requirement with air guns, which excite P waves to illuminate the subsurface. To reduce computational costs, most industrial-scale 3-D FWI implementations of exploration inversion rely on pseudo-acoustic (qP wave) physics (e.g. Alkhalifah 2000; Gholami et al.2013a,b; Kamei et al.2013). However, the reflection coefficients cannot be well approximated by pseudo-acoustic modeling, especially in acquisitions limited by short offsets (e.g. Raknes & Arntsen 2014; Ovcharenko et al.2018). Since elastic anisotropic waveform inversions have higher computational costs to their isotropic counterparts, they only recently became the focus of attention (Köhn et al.2015; Podgornova et al.2015; Owusu et al.2016; Bergslid et al.2015; Oh & Alkhalifah 2016; Kamath & Tsvankin 2016). Feasibility and resolution studies for 3-D inversions are still computationally expensive, therefore they are often performed either in 2-D (Podgornova et al.2015) or with very few simultaneous (encoded) shots (Köhn et al.2015). Inversion for the full set of orthorhombic parameters in the whole subsurface model increases the likelihood of overfitting the data. If the orthorhombic anisotropy is restricted to the reservoir area of the subsurface model and the background is assumed to be isotropic, then the number of inverted parameters decreases and, thus, the inversion becomes more robust. When the anisotropy is restricted to the area of the reservoir, the scattering is well approximated by Born scattering in the isotropic background (Hudson & Heritage 1981). In Born approximation, scattering of plane waves can be fully characterized by a scattering function (Eaton & Stewart 1994), which is a function of the incident plane wave parameters and the perturbation of the elastic parameters. If the parameter perturbation and the types of incident and scattered waves are fixed, then the scattering function is reduced to a diffraction-based radiation pattern (Gholami et al.2013b; Alkhalifah & Plessix 2014). Finally, if the incident and scattered angles in the diffraction-based radiation pattern obey Snell’s law for horizontal reflectors, then the diffraction-based radiation pattern becomes a reflection-based radiation pattern (Alkhalifah & Plessix 2014). With the assumption of an isotropic background, the resolution can be investigated visually using the reflection-based radiation patterns (Alkhalifah & Plessix 2014; Oh & Alkhalifah 2016) and quantitatively using the singular value decomposition (SVD) analysis of those patterns (Eaton & Stewart 1994; de Hoop et al.1999; Podgornova et al.2015; Kazei & Alkhalifah 2017). Eaton & Stewart (1994) were the first to discover the importance of the scattering function in anisotropic media and established a resolution analysis for vertically incident qP waves in vertically transversally isotropic (VTI) media. de Hoop et al. (1999) derived a similar analysis for orthorhombic parameter perturbations that focused on reflections in an isotropic background. Reflection-based radiation patterns were used to analyse the ability of anisotropic elastic FWI to estimate VTI parameters using only P waves in various parametrizations (Gholami et al.2013b; Alkhalifah & Plessix 2014). Podgornova et al. (2015) proved that only three VTI parameters are invertible from P–P waves under the assumption of weak anisotropy, by using SVD analysis to investigate the FWI resolution possibilities with P–P and P–SV waves in Thomsen’s well-known parametrization for VTI media (Thomsen 1986). Thomsen’s parametrization was generalized to orthorhombic media by Tsvankin (1997) and rewritten with deviation parameters by Masmoudi & Alkhalifah (2016) and Oh & Alkhalifah (2016) to obtain a hierarchical anisotropy parametrization. Oh & Alkhalifah (2016) came to the conclusion that their parametrization minimizes trade-offs between parameters and, for P–P waves in most acquisition scenarios, allows the inversion of four out of the 10 independent elastic orthorhombic parameters. The traveltimes of P waves contain information on only six linear combinations of these parameters (Tsvankin 1997). A numerical resolution study (Köhn et al.2015) has shown that, with P-wave dominant sources, three-component measurements and perfect illumination, inversion for all 10 orthorhombic parameters may be feasible, but the general acquisition scenarios were not addressed. Here, we investigate the relation between the acquisition and the number of invertible parameters assuming an explosive source and an early arrival amplitude, namely, P waves recorded by hydrophones. In Section 2, we provide and discuss the derivation of an anisotropic scattering function and clarify its relation to radiation patterns. In Section 3, we map 3-D radiation patterns into the wavenumber domain and discuss the resolution qualitatively. In Section 4, we introduce the effective angle of illumination concept and use it to compare a recent hierarchical parametrization suggested by Masmoudi & Alkhalifah (2016) for acoustic case and by Oh & Alkhalifah (2016) for elastic with a classic stiffness matrix parametrization of orthorhombic media. In Section 5, we apply an SVD analysis to the sensitivities of P–P scattered waves to derive the number of parameters that can be inverted in different acquisition scenarios. 2 LINEARIZED SCATTERING, RADIATION PATTERNS The linearization of the wave equation with respect to parameter perturbation provides the crucial first-order sensitivity of the data to the model. The analysis of such sensitivity using plane wave and SVD provides us with the scattering radiation patterns, and insights into the model invertibility considering the limited illumination provided by conventional recording of seismic data. 2.1 The Born approximation Most methods for dynamic inversion of seismic wavefields, including linearized reverse-time migration and every iteration of FWI, rely on the Born approximation. Thus, the sensitivity of the wavefields, which defines the resolution of these methods, can be evaluated within the Born approximation. We analyse the resolution of dynamic inversions for orthorhombic parameters relying on the linearized scattering of an incident wavefield U0 on a perturbation of the density (δρ(x)) and stiffness tensor (δckjpq(x)) parameters. The frequency-domain Born approximation for the scattered wavefield δU in an anisotropic elastic medium is (e.g. Hudson & Heritage 1981; Beylkin & Burridge 1990; Shaw & Sen 2004)   \begin{eqnarray} \delta U_i (\mathbf {x}_s,\mathbf {x}_\mathbf {g},\omega ) &=& \int \limits _{V} \left(\delta \rho (\mathbf {x}) u^0_k \omega ^2\right.\nonumber\\ \left. - \delta c_{kjpq}(\mathbf {x}) \frac{\partial U^0_p}{\partial x_q} \frac{\partial }{\partial x_j} \right)G_{ik}(\mathbf {x}_g, \mathbf {x}) {\rm d}\mathbf {x}, \end{eqnarray} (1)where ω is the angular frequency, xs and xg are the coordinates of the source and receiver, respectively and Gik(xg, x) is the Green’s tensor for the background medium. Integration is performed over the volume of the domain where the perturbation is located. Born approximation (eq. 1) is valid for relatively small   \begin{eqnarray} \frac{|\delta \mathbf {m}|}{|\mathbf {m}|} << 1 \end{eqnarray} (2)and local   \begin{eqnarray} \frac{\omega d}{V_0}\frac{|\delta \mathbf {m}|}{|\mathbf {m}|} << 1 \end{eqnarray} (3)perturbations, where m is the vector of unperturbed elastic parameters, V0 is the smallest background velocity (Vs in the isotropic case) and d is the diameter of the perturbation. The above assumptions are best met in the case of low-contrast point scatterers, yet the Born approximation does not imply specific shape restrictions on the scattering perturbation. For instance, it is also valid for the perturbations within thin low-contrast layers and near-normal wave incidence angles (Shaw & Sen 2004). To obtain a simplified expression for the scattering of plane waves on an arbitrary perturbation, we first manipulate the right-hand side of eq. (1). Next, we discuss the cases of point scatterers and horizontal reflectors. Finally, we examine the general case through spectral sensitivities. For simplicity, we consider the far-field approximation from both the source and receiver sides, with an isotropic background. 2.2 Far-field plane-wave scattering Most FWI case studies focus on marine data sets due to their high signal-to-noise ratios. In marine acquisition, only P waves propagate through water and, thus, they dominate the incident wavefields in the target area of exploration. Therefore, we consider only P waves using the following far-field plane-wave approximation for the incident wavefield:   \begin{eqnarray} \mathbf {U}^0(\mathbf {x}_s,\mathbf {x},\omega ) &=& \alpha _s e^{i\frac{\omega }{v_p}||\mathbf {x}- \mathbf {x}_s||}\mathbf {s}\simeq \alpha _s e^{i\frac{\omega }{v_p}\mathbf {s}\cdot (\mathbf {x}- \mathbf {x}_g)}\mathbf {s},\, \alpha _s\nonumber\\ &=& \frac{1}{4\pi \rho v_p^2||\mathbf {x}- \mathbf {x}_s||} \end{eqnarray} (4)where ρ and vp are the background density and longitudinal wave velocity, respectively, and s is a unit vector pointing towards the source from the target area. The far-field part of the Green’s tensor in eq. (1) for isotropic background consists of P, SV and SH waves:   \begin{eqnarray} \mathbf {G}= \alpha _g (\mathbf {G}_P + \mathbf {G}_{SV} + \mathbf {G}_{SH}),\, \alpha _g = \frac{1}{4\pi \rho v_p^2||\mathbf {x}- \mathbf {x}_g||}. \end{eqnarray} (5)The coefficient αg is dropped in our resolution analysis, since it is related to geometrical spreading and can be effectively compensated using pseudo-Hessian techniques. Different components of the Green’s function are approximated locally by plane waves. For the first component GP, the approximation leads to:   \begin{eqnarray} \mathbf {G}_P = e^{i\frac{\omega }{v_p}||\mathbf {x}- \mathbf {x}_g||} \mathbf {g}\mathbf {g}\simeq e^{i\frac{\omega }{v_p}\mathbf {g}\cdot (\mathbf {x}- \mathbf {x}_g)}\mathbf {g}\mathbf {g} \end{eqnarray} (6)which is a plane P wave. (gg)ik = gigk is the outer product of vector g, pointing towards receiver from the perturbation, with itself. For the second and third components, GSV and GSH, similar expressions can be provided (Snieder 2002), yet we keep them out of the scope of this paper in order to focus on the early P-wave arrivals. By substituting eqs (4) and (5) into eq. (1) and projecting on g, we obtain the approximate amplitude of the scattered P–P waves:   \begin{eqnarray} \delta U_{PP} \equiv \frac{\delta \mathbf {U} \cdot \mathbf {g}}{\alpha _{PP}} \simeq \int \limits _{V}e^{i\mathbf {K}_{PP}\cdot \mathbf {x}} (\mathbf {s}\cdot \mathbf {g}\delta \rho - \mathbf {s}\mathbf {s}: \delta \mathbf {c}:\mathbf {g}\mathbf {g}) {\rm d}\mathbf {x} \end{eqnarray} (7)  \begin{eqnarray} =\mathbf {s}\cdot \mathbf {g}\int \limits _{V}e^{i\mathbf {K}_{PP}\cdot \mathbf {x}} \delta \rho {\rm d}\mathbf {x}- \mathbf {s}\mathbf {s}: \int \limits _{V}e^{i\mathbf {K}_{PP}\cdot \mathbf {x}} \delta \mathbf {c}{\rm d}\mathbf {x}: \mathbf {g}\mathbf {g}, \end{eqnarray} (8)where   \begin{eqnarray} \alpha _{PP} = \frac{\omega ^2}{v^2_p} \alpha _s \alpha _g e^{i\frac{\omega }{v_p}(-\mathbf {x}_s\mathbf {s}-\mathbf {x}_g\mathbf {g})}, \end{eqnarray} (9)  \begin{eqnarray} \mathbf {K}_{PP} = \frac{\omega }{v_p}(\mathbf {s}+\mathbf {g}). \end{eqnarray} (10)Our further resolution analysis is based on eq. (8), which serves as the simplified monochromatic scattering approximation. 2.3 Scattering function A monochromatic plane wave is determined by its phase, polarization and propagation direction. Since the phase is dropped, scattering depends only on the incident plane wave direction s, scattered plane wave direction g, their respective polarizations ς and ξ, the frequency ω, and on the perturbation of the parameters,   \begin{eqnarray} \delta \mathbf {m}= (\delta \rho (\mathbf {x}), \delta \mathbf {c}(\mathbf {x}))^T. \end{eqnarray} (11)Namely, the scattered wavefield from a point perturbation δm = δ(x)(δρ, δC)T in the first-order Born approximation is proportional to the scattering function (Beylkin & Burridge 1990; Eaton & Stewart 1994; Shaw & Sen 2004):   \begin{eqnarray} \delta U (\mathbf {s},\mathbf {g}, \boldsymbol {\varsigma }, \boldsymbol {\xi }, \delta \mathbf {m}) \propto R\equiv \varsigma _k \cdot \xi _k \delta \rho - s_i \varsigma _j \delta c_{ijkl} g_k \xi _l \end{eqnarray} (12)  \begin{eqnarray} = \boldsymbol {\varsigma }\cdot \boldsymbol {\xi }\delta \rho - \mathbf {s}\boldsymbol {\varsigma }: \delta \mathbf {c}: \mathbf {g}\boldsymbol {\xi }. \end{eqnarray} (13)The scattering function R is independent of the background medium, which can have up to triclinic anisotropy complexity. As there are too many possible cases of polarization and propagation relations in general anisotropic background, we consider scattering in an isotropic background medium. 2.4 Classic diffraction-based scattering radiation pattern—$$\mathcal {D}$$ For a given conversion type (P–P, P–SH, etc.), the scattered wavefield becomes a function of two vectors, s and g, called the scattering radiation pattern:   \begin{eqnarray} \mathcal {D}_{WI-WS, \delta \mathbf {m}}(\mathbf {s},\mathbf {g}) \equiv \frac{ R(\mathbf {s}, \mathbf {g}, \boldsymbol {\varsigma }(\mathbf {s}, WI), \boldsymbol {\xi }(\mathbf {g}, WS), \delta \mathbf {m})}{||\delta \mathbf {m}||}. \end{eqnarray} (14)Here, WI denotes the type of incident, and WS denotes the type of scattered waves (P − P, P − Si, Si − P or Si − Sj), for a given parameter perturbation δm. The radiation pattern characterizes the scattering of certain types of plane waves on a point-like perturbation of the parameters. For vertically propagating incident P waves, the scattering radiation patterns become functions of two real arguments, that is, g vector components, and can be presented as a 2-D plot (Eaton & Stewart 1994). We show the diffraction radiation pattern for P–P scattering on a density perturbation in Fig. 2(a) for the vertical incidence direction (e.g. Tarantola 1986). Figure 2. View largeDownload slide (a) The diffraction pattern for density perturbation. (b)  Narrow azimuth—ϕ acquisition and reflection from the reservoir with scattering angle θ. (c)  The reflection based radiation pattern. and (d)  its remapping into the spectral domain. The normalized wavenumber ||K||∝cos θmax. In all the patterns, the amplitude sign is shown through colour, ‘+’ corresponds to positive reflection coefficient—no phase shift after scattering, ‘−’—shows the flip of the phase—negative reflection/scattering coefficient. Figure 2. View largeDownload slide (a) The diffraction pattern for density perturbation. (b)  Narrow azimuth—ϕ acquisition and reflection from the reservoir with scattering angle θ. (c)  The reflection based radiation pattern. and (d)  its remapping into the spectral domain. The normalized wavenumber ||K||∝cos θmax. In all the patterns, the amplitude sign is shown through colour, ‘+’ corresponds to positive reflection coefficient—no phase shift after scattering, ‘−’—shows the flip of the phase—negative reflection/scattering coefficient. 2.5 Reflection-based radiation patterns—$$\mathcal {R}$$ Classic scattering radiation patterns have focused on diffractions from point scatterers. However, seismic surveys are commonly dominated by reflections (Fig. 2b), which prompted the development of linearized reflection coefficients (e.g. Rüger 1997) and, more recently, reflection-based radiation patterns (Shaw & Sen 2004; Gholami et al.2013b; Alkhalifah & Plessix 2014). The relation between the reflection patterns and linearized coefficients is well explained by Shaw & Sen (2004). In 3-D scattering experiments, scattering functions depend on four scalar parameters defined by the unit vectors, s and g. For a plane reflector, the scattered wave direction is determined by the incident direction and vice versa through Snell’s law. The dip and azimuth of the reflector determine how the s is tied to the g in particular, reducing the number of real variables that determine scattering from four to two:   \begin{eqnarray} \mathcal {R}_{WI-WS, \delta \mathbf {m}}(\mathbf {g}) \equiv \mathcal {D}_{WI-WS, \delta \mathbf {m}}(\mathbf {s}(\mathbf {g}),\mathbf {g}). \end{eqnarray} (15)The amplitudes of the reflections from perturbations of isotropic parameters do not depend on the reflector orientation; therefore, the reflection-based radiation patterns are the same for all reflector orientations. In anisotropic cases, the reflection-based radiation patterns depend on the dip and azimuth angles. In most seismic cases, the media is dominated by horizontal layering; for this reason, Gholami et al. (2013b) and Alkhalifah & Plessix (2014) focused on the horizontal reflectors. If we assume P–P scattering on a horizontal reflector, then Snell’s law can be expressed as (si and gi are projections of vectors s and g, respectively)   \begin{eqnarray} s_1&=& -g_1 = \sin \theta \cos \phi ,\, s_2=-g_2=\sin \theta \sin \phi ,\nonumber\\ s_3&=&g_3=\cos \theta /2. \end{eqnarray} (16)The radiation patterns become functions of two real variables, ϕ and θ, and can be plotted. The reflection-based radiation pattern for density and P–P scattering is represented in Fig. 2(c). We show the reflection-based radiation patterns for other cases mapped from the ϕ, θ domain into wavenumber domain (ϕ, kz = 2cos θ/2) as spectral sensitivities, to directly access the resolution in the normalized wavenumber domain. 3 SPECTRAL SENSITIVITIES Spectral sensitivities allow us to consider a broad class of anomalies that are small enough that we can neglect the second-order scattering. The relation between the wavenumbers illuminated in a perturbation and the directions of incident and scattered waves in a monochromatic plane-wave experiment has been recognized and successfully applied in the X-ray tomographic method of Ewald (1969) spheres. Devaney (1984) applied this method to geophysical diffraction tomography resolution quantification in the wavenumber domain. The technique was later applied in acoustic FWI scenarios for surface seismic acquisition (Mora 1989) and vertical seismic profiling (Wu & Toksöz 1987). Kazei et al. (2013a) introduced amplitude of the spectral sensitivity of data to the perturbations of various background media, and applied it to models with diving waves and multiples (Kazei et al.2015, 2013b). The elastic isotropic background case with VTI perturbations was considered by Podgornova et al. (2015). Here, we show the applicability of this technique to the general case of anisotropic perturbations of isotropic media. The possibility of recovering a set of parameters is defined by the sensitivity of the scattered wavefield to these parameters. Under the plane-wave assumption, the resolution is straightforward for the wavenumber-domain representation of perturbation. The concept of spectral sensitivities plays an important role in our analysis. 3.1 Definition The right-hand side of eq. (8) is, essentially, a Fourier transform of the perturbation and can be rewritten to reveal the nature of scattering in a far-field approximation as follows:   \begin{eqnarray} \delta U_{PP}(\mathbf {s},\mathbf {g},\omega ) = \mathbf {s}\cdot \mathbf {g}\delta \hat{\rho }(\mathbf {K}_{PP}) - \mathbf {s}\mathbf {s}: \delta \hat{\mathbf {c}}(\mathbf {K}_{PP}) :\mathbf {g}\mathbf {g}. \end{eqnarray} (17)If we consider the spectrum of perturbation parameters (δm(K)), then we see that:   \begin{eqnarray} \delta U_{PP} = \mathcal {R}_{PP} (\mathbf {s}, \mathbf {g}) ||\delta \mathbf {m}(\mathbf {K}_{PP})||, \end{eqnarray} (18)which allows us to represent radiation patterns in the normalized wavenumber domain:   \begin{eqnarray} \mathcal {S}_{\delta \mathbf {m}} (\mathbf {k}_{PP}, \varphi ) \equiv \mathcal {R}(\mathbf {s}(\mathbf {k}_{PP},\varphi ), \mathbf {g}(\mathbf {k}_{PP},\varphi )), \end{eqnarray} (19)where   \begin{eqnarray} \mathbf {k}_{PP} = \mathbf {K}_{PP}/k_0 = \mathbf {s}+\mathbf {g},\, k_{0} = \frac{\omega }{v_p}, \end{eqnarray} (20)Eq. (19) defines the spectral sensitivity for parameter δm. It is more convenient than the reflection-based radiation patterns as it quantifies the resolution for different normalized wavenumbers kPP directly. We normalize all of the wavenumbers, dividing them by ω/Vp, acknowledging that the sensitivities scale linearly with frequency in homogeneous backgrounds (Kazei et al.2013a). We focus on the inversion for vertical wavenumbers kPP = (0, 0, kz), or horizontally layered structures, which are dominant in most geological scenarios, returning us to the remapped reflection-based radiation patterns in a more convenient representation for resolution analysis domain. These remapped patterns are functions of two variables (kz, ϕ) and, therefore, can be easily plotted and examined. Snell’s law gives the following expressions in the case in which the incident and scattered directions are functions of ϕ and kz:   \begin{eqnarray} s_1&=& -g_1 = \sqrt{1-\frac{k_z^2}{4}}\cos \phi ,\, s_2=-g_2=\sqrt{1-\frac{k_z^2}{4}}\sin \phi ,\nonumber\\s_3&=&g_3=k_z/2. \end{eqnarray} (21)In Fig. 2(d), we show the spectral sensitivity for a perturbation in density when all of the elastic constants are kept fixed. There is no scattering at intermediate reflection angles (Fig. 2a, so the intermediate spatial wavenumbers for density do not contribute to the scattered wavefield and cannot be reconstructed. However, low and high wavenumbers can be reconstructed. This is different from the standard inversion for density where Vp is fixed and only the high wavenumbers show significant scattering. In the next section, we analyse the scattering of all the elastic constants and try to estimate the null space of the inversion. 3.2 Spectral sensitivity patterns Cij, ρ parametrization We examine the remapped radiation patterns with the newly introduced representation in Fig. 3. Tsvankin (1997) established that the traveltimes of P waves in orthorhombic media depend on six parameters. This is the number of parameters that kinematic methods can resolve in perfect illumination conditions. One might expect that the inversion of the amplitudes could recover more parameters, since more attributes of the wavefields are taken into account. However, the visual analysis of the P–P radiation patterns does not provide detailed information on how the Cij parameters are coupled. Figure 3. View largeDownload slide The P–P spectral sensitivity patterns for various Cij perturbations. The azimuth and Kz values are mapped the same way as for density which is shown in Fig. 2(d). It is clear that the sensitivity to all the orthorhombic parameters is non-zero if we start from isotropic background. Yet we can note some coupling effects already. C13 is very similar to C55, C12—to C66. Figure 3. View largeDownload slide The P–P spectral sensitivity patterns for various Cij perturbations. The azimuth and Kz values are mapped the same way as for density which is shown in Fig. 2(d). It is clear that the sensitivity to all the orthorhombic parameters is non-zero if we start from isotropic background. Yet we can note some coupling effects already. C13 is very similar to C55, C12—to C66. To analyse the scattering, we generate the spectral sensitivities for single frequency, infinite offset data. According to eq. (10), for single frequency infinite offset data, P waves cover vertical wavenumbers from 0 to $$\frac{2\omega }{V_p}$$. In Fig. 3, we see four non-orthorhombic elements with non-zero sensitivities. These non-zero sensitivities can be used to invert for monoclinic medium parameters with a horizontal symmetry plane. At the same time, we admit that the sensitivity to tilt-related parameters (rotations of perturbations around a horizontal axis) is zero in an isotropic background; therefore, the inversion for tilt from P waves does not work for pure vertical wavenumbers in isotropic backgrounds. Sensitivities, that is, remapped reflection-based radiation patterns, for different Cij parameters provide some information about which parameters cannot be inverted for; yet they provide almost no information about the coupling of parameters, as there are usually more than two parameters involved in the coupling, as we show in the next section. 3.3 Sensitivities in a hierarchical parametrization While necessary for the characterization of carbonate reservoirs, orthorhombic anisotropy is not observed in vast areas of the subsurface. An easy way to set up an inversion with different anisotropy restrictions in different domains of the subsurface is to split the parameters into three groups using a hierarchical parametrization based on the type of anisotropy that they describe. Such a parametrization, introduced by Oh & Alkhalifah (2016), uses Vp, Vs, ρ —for the isotropic part of the medium, where Vp and Vs are the velocities of P and S waves in isotropic media, respectively, and ρ is the density. In general orthorhombic media, the Vp and Vs considered here are the velocities with wave polarizations along the x1 horizontal axis. ε1, η1, γ1 —for the VTI part of the anisotropy, where ε1 and γ1 is the relative difference in velocities between vertically and horizontally propagating P and SH waves, respectively, and η1 is the Alkhalifah–Tsvankin anellipticity parameter in [x1, x3] plane. εd, ηd, γd, δ3 —are parameters that define deviation of the medium from VTI anisotropy. They are non-zero only if the medium has fully orthorhombic complexity. ηd = η2 − η1, where η2 is the Alkhalifah–Tsvankin (Alkhalifah 2000) anellipticity parameter in [x1, x3] plane. Similarly, γd = γ2 − γ1 and εd = ε2 − ε1 define differences in VTI anisotropy in [x2, x3] and [x1, x3] planes. δ3 is the well-known Tsvankin–Thomsen parameter which is zero in case the medium is VTI. Radiation patterns for any parametrization can be obtained through the chain rule   \begin{eqnarray} \boldsymbol {\mathcal {R}}_{\mathbf {m}, WT} = \frac{\partial \mathbf {C}}{\partial \mathbf {m}}:\boldsymbol {\mathcal {R}}_{\mathbf {C}, WT}. \end{eqnarray} (22)The Jacobian $$\frac{\partial \mathbf {c}}{\partial \mathbf {m}}$$ (eq. B15) reveals the actual perturbations of the Cij parameters corresponding to individual perturbations in the new parametrization. Therefore, we also represent the partial derivatives of Cij parameters graphically in Fig. 4(a) to understand the analysis better. Figure 4. View largeDownload slide (a) The partial derivatives of the Cij parameters (same scheme as in Fig. 1 b is used) in the Cij, ρ parametrization with respect to the Oh & Alkhalifah (2016) parameters. (b)  The spectral sensitivities of the new parameters. Figure 4. View largeDownload slide (a) The partial derivatives of the Cij parameters (same scheme as in Fig. 1 b is used) in the Cij, ρ parametrization with respect to the Oh & Alkhalifah (2016) parameters. (b)  The spectral sensitivities of the new parameters. Spectral sensitivities in the new parametrization elucidate some features of scattering that were not visible in the Cij, ρ parametrization. For example, the resemblance of radiation patterns for density and ε1 gives a hint that there is a trade-off between these parameters, which was hidden in several radiation patterns before. At the same time our knowledge is still limited by the choice of parametrization as we cannot analyse the coupling between three or more parameters. 4 EFFECTIVE ANGLE OF ILLUMINATION To perform a successful multiparameter inversion, we will need to retrieve all parameters in the whole domain of interest to avoid the dreaded null space. The capability of inverting of decoupled parameters can be investigated for a given vertical wavenumber, separately. The scattering of monochromatic plane waves for different wavenumbers in a perturbed medium is determined by the scattering angle–wavenumber relation (Mora 1989; Alkhalifah 2016; Kazei et al.2016) using a locally homogeneous medium assumption. Utilizing this property, we are going to examine the ability of the inversion in decoupling parameters for particular vertical wavenumber Kz. For single frequency and single azimuth data, there is a bijective relation between the wavenumber illuminated (two real parameters) and the incidence and scattered directions (two scalar angles), which means that the scattered wavefield maps directly to the parameter at this vertical wavenumber. Therefore, in principle, single frequency, single component data can resolve only one parameter. Unless we impose some relations between parameters at different wavenumbers with, for example, some type of regularization. In order to decouple VTI media parameters for a given wavenumber, we need a set of frequencies (Podgornova et al.2015) and azimuths contributing to this wavenumber illumination. While for all vertical wavenumbers, we assume the azimuth coverage is the same there might be additional angle-offset limitations due to frequency content of the data. These limitations are summarized into the effective angle of wavenumber illumination concept in this section. The scattering of plane waves can be expressed in the following system of equations:   \begin{eqnarray} \delta \mathbf {U}(\omega _i,\theta _j,\varphi _k) \propto \mathcal {R}(\theta _j,\varphi _k) \delta \mathbf {m}(\mathbf {K}(\theta _j,\omega _i)), \end{eqnarray} (23)which can be regrouped for a given wavenumber K, using the angle–wavenumber relation:   \begin{eqnarray} \mathbf {u}(\omega _i(\mathbf {K}),\theta _i(\mathbf {K}),\varphi _k) \propto \mathcal {R}(\theta _i(\mathbf {K},\omega ),\varphi _k) \delta \mathbf {m}(\mathbf {K}). \end{eqnarray} (24)Thus, data for each vertical wavenumber of the perturbation are available from different frequencies ωi, which gives us different opening angles θi = (K, ωi) for illumination of the wavenumber K. Having different angles of illumination for the same wavenumber K allows us to distinguish parameters at the scale given by this wavenumber through variation in the scattering patterns of parameters. In Fig. 5, we show that the set of angles ‘active’ in a radiation pattern for a given wavenumber are defined by the offset as well as by the frequency content of the signal. Figure 5. View largeDownload slide The same wavenumber can be illuminated from different opening angles and by different wave modes. Each angle availability is determined by the extent of the angle in the acquisition and existence of the appropriate frequencies to scale the source and geophone vectors. Figure 5. View largeDownload slide The same wavenumber can be illuminated from different opening angles and by different wave modes. Each angle availability is determined by the extent of the angle in the acquisition and existence of the appropriate frequencies to scale the source and geophone vectors. We examine the effective angle of illumination concept on P–P waves as the relations are simple in this case. According to eq. (20) and Fig. 5, the angles that contribute to the illumination of a given vertical wavenumber can be restricted by both the frequency and the aperture. If the recorded band of frequencies is narrow then the angles of illumination are determined by the bandwidth of the data. In this case, there is no benefit received from an aperture increase in decoupling the parameters; that is, there is no reason to increase the aperture while the frequency content of the data remains the same. For example, in $$K_z\in [\frac{2\omega }{v_p}\cos \frac{\theta _{{\rm max}}}{2};\frac{2\omega }{v_p}]$$ (the red region in Fig. 10a), the angles of illumination are not affected by the aperture at all. If a wide band of frequencies is recorded the effective angle of illumination is generally determined by the aperture. Namely, for the range of wavenumbers Kz between $$\frac{2\omega }{v_{{\rm min}}}$$ and $$\frac{2\omega _{{\rm max}}}{V_p}\cos \theta _{{\rm max}}$$, the full offset length can be exploited for inversion and the decoupling of parameters, that is, in the green area in Fig. 6(b), the resolution is defined by the aperture only. The effective angle of illumination concept suggests that the data frequency bandwidth and the aperture should be increased simultaneously to achieve better decoupling, needed for the inversion. In the next section, we assume that the bandwidth is sufficient and consider the resolution for the green area in Fig. 6(b). Figure 6. View largeDownload slide For a given vertical wavenumber, angles that contribute to its illumination can be restricted by frequency as well as by aperture. (a) A narrow band of frequencies. In the red area, the resolution of wavenumbers is defined by the frequency content. (b) A wide band of frequencies. In the green area, the resolution is defined by the illumination angles. Wavenumbers between $$\frac{2\omega _{{\rm min}}}{v}$$ and $$\frac{2\omega _{{\rm max}}}{v}\cos \theta _{{\rm max}}$$ can be inverted using the whole offset length. Figure 6. View largeDownload slide For a given vertical wavenumber, angles that contribute to its illumination can be restricted by frequency as well as by aperture. (a) A narrow band of frequencies. In the red area, the resolution of wavenumbers is defined by the frequency content. (b) A wide band of frequencies. In the green area, the resolution is defined by the illumination angles. Wavenumbers between $$\frac{2\omega _{{\rm min}}}{v}$$ and $$\frac{2\omega _{{\rm max}}}{v}\cos \theta _{{\rm max}}$$ can be inverted using the whole offset length. 4.1 Sensitivity matrices $$\mathbb {A}_{P-P}$$ – assembly The main outcome of the previous subsection is that particular wavenumber illumination can be limited not only by the lack of aperture in the data, but also by the band of the signal of the source and geophones. In this subsection, we discuss the matrix form of the simplified inversion problem for the vertical wavenumber K from the P–P scattering mode to prepare ourselves for quantitative resolution analysis. For an arbitrary perturbation of parameters, the scattered wavefield can be expressed in the wavenumber domain as   \begin{eqnarray} \delta U_{PP} (k_z, \phi , \omega (k_z, K_z)) \propto \mathbf {A}_{P-P}^T (k_z, \phi ) (\delta \mathbf {m}(\mathbf {K})), \end{eqnarray} (25)or in the vectorized form as   \begin{eqnarray} \delta \mathbf {U}_{PP}(K_z) \propto \mathbb {A}_{P-P} \delta \mathbf {m}(K_z), \end{eqnarray} (26)where kz ∈ [2cos (θmax/2), 2]. Vector A(kz, ϕ) represents the set of reflection-based radiation patterns remapped into the normalized wavenumber domain. For example, element A1(kz, ϕ) (Fig. 2d) is the radiation pattern of the density, while the radiation patterns for the 21 elastic constants (the total number of Voigt elastic parameters is 22 – 21 elastic constants and density), are shown in Fig. 3, which is ordered to correspond to the Voigt notation (Fig. 1b). We can clearly see couplings in pairs (resemblance in patterns) : C12 and C66, C13 and C55 and C23 and C44. We start by analysing the linearized problem that needs to be solved in order to reconstruct a low vertical wavenumber of perturbation, since low wavenumbers are considered the hardest and most essential wavenumbers for FWI to converge. First: let us assume that the frequency content is sufficient, and the illumination of the wavenumber is limited by the aperture (Fig. 6b). Then, there will be angles of illumination [0; θmax] available for a given wavenumber K, where θmax is determined by the aperture. First, sticking with the elastic coefficients and density parametrization, we seek to understand the trade-off and extract the linear combinations that expose such trade-off. Similarities of certain radiation patterns in Fig. 3 (e.g. C12 and C66) suggest the presence of a null space. The rank of the matrix $$A_{ik} = \mathcal {R}(\theta _i(\mathbf {K},\omega ),\varphi _k)$$ is equal to the number of parameters resolvable for a given wavenumber K value. The dimensionality of the null space is equal to the total number of parameters minus the rank of the matrix. Note that the part of the radiation patterns that plays role in this analysis is commonly attributed to high resolution in other studies (e.g. Alkhalifah & Plessix 2014), because the wavenumbers obtained from small offsets are higher than those from far offsets for a given frequency. The frequency is not an issue since we consider a wide band. The low frequencies in the data are useful for monoparameter FWI ten Kroode et al. (2013), yet for multiparameter they are necessary to resolve trade-off between parameters. When there is a lack of high frequencies, the effective angle range for inversion for large wavenumbers activates the high wavenumber part of the spectral sensitivity. Thus, inversion for low wavenumbers with limited offset data is similar to inversion for high wavenumbers with a limit in the highest frequency available in the data. Second: we assume that lack of low frequencies limits the resolution of low wavenumbers in the model. The offset is assumed to be effectively infinite in this case, for example, diving waves illuminate the whole region of interest (Kazei et al.2013b). Diving waves, propagating horizontally at the depth of perturbation, activate the ‘transmission’ part of the radiation patterns with wide opening angles (Alkhalifah & Plessix 2014), equivalent to low vertical wavenumbers in spectral sensitivities. The second case is common when the offsets are very far and most low wavenumbers are considered. For P–P scattering the normalized wavenumber kz ∈ [0, 2ϰ cos (θmin − eff − P)], where $$\cos (\theta _{{\rm min}-{\rm eff}-P}) = \frac{K_z v_p}{2\omega _{{\rm min}}}<1$$. Third: the most practical and the worst in terms of the data coverage case is when the illumination of a wavenumber is limited on both ends. This inevitably happens for very low wavenumbers when the aperture is limited, that is, there are no diving waves traveling at the depth of perturbation. For this case, there are two limiting variables and for analysis of the plots we assume that the azimuth is either full or single. 4.2 Matrix A—SVD, Hessian To perform a quantitative analysis of coupling between orthorhombic parameters at a certain resolution determined by vertical wavenumber, we compute the SVD of matrix $$\mathbb {A}$$:   \begin{equation*} \mathbb {A}= \mathbb {U}\mathbb {S}\mathbb {V}^T \end{equation*} Matrix $$\mathbb {V}$$ is a 10 × 10 matrix which consists of normalized eigenvectors of $$\mathbb {A}^T\mathbb {A}$$. If the standard quadratic misfit J of the observed and modeled data (e.g. Tarantola 1984) is optimized by inversion,   \begin{eqnarray} J = ||\delta \mathbf {U}- \mathbb {A}\delta \mathbf {m}||^2 \end{eqnarray} (27)then the Hessian of the fitting problem is $$\mathbb {A}^T \mathbb {A}$$. The columns of $$\mathbb {V}^T$$ are the eigenvectors of the Hessian. Diagonal elements of the matrix $$\mathbb {S}$$ are the eigenvalues of the Hessian. Thus, $$\mathbb {V}$$ tells us which linear combinations of parameters are invertible from the data, and $$\mathbb {S}$$ tells the difference between the hardships to invert for these linear combinations. One of the most important features of the eigenvectors of the Hessian is that they are orthogonal with respect to the misfit functional norm. Thus, if we choose the parametrization in which each singular vector has only one non-zero component, we could invert all of the parameters one by one without any cross-talk issues. The problem of such a parametrization is that it would be spatially dependent as the set of singular vectors depends on the illumination of a particular domain. Therefore, there is no perfect parametrization for all of the cases and it can only be tuned according to the target depth and the acquisition parameters. Matrix $$\mathbb {S}$$ can have less than 10 non-zero values on the diagonal and, in that case, some trade-offs between parameters are irresolvable. Ties between the parameters need to be introduced (i.e. anisotropy complexity reduction or some substantial parameters, such as Gardner’s law) or additional data must be taken into inversion, for example, other wave types, etc. Even when $$\mathbb {S}$$ has 10 non-zero values, the number of invertible parameters might be reduced by the presence of noise in the data. The noise is accounted for in our inversion resolution study by using the model resolution matrices (Menke 2012; Podgornova et al.2015) with different thresholds for the singular values. The resolution matrix is calculated using the formula   $\mathbb {R}= \mathbb {V}^T (\mathbb {S}\mathbb {S}^T)^{-1} \mathbb {V},$ where $$(\mathbb {S}\mathbb {S}^T)^{-1}$$ is the generalized inverse of the matrix with threshold for singular values, which is determined by the amplitude of noise in the data. In the next section, we apply the apparatus described above to the sensitivity patterns. 5 RESULTS In this section, we quantitatively explore the following questions: What is the null space of linearized dynamic inversions? How many orthorhombic parameters can be resolved in realistic cases? What are these parameters? Of course, answers to these questions depend on the acquisition and data signal-to-noise ratio. Generally, the resolution quantification can be applied this way: Choose Kz of interest at a target depth. From bandwidth and aperture, estimate the range of effective angle of illumination (Fig. 6b). Define the set of effective incidence vectors s (θ, ϕ → s) using the azimuth coverage and range effective angle of illumination. Assemble matrix $$\mathbb {A}$$, which maps parameters into radiation patterns. Compute the SVD of matrix $$\mathbb {A}$$. Interpret the SVD and compute the resolution matrices. We focus on low wavenumbers recovery and consider three options for effective angle limitations: Perfect illumination—there are diving waves up to the level of perturbation and there are frequencies low enough to illuminate the wavenumber by backscattering. The null space in case of perfect illumination defines principally irresolvable trade-offs in between parameters. Lack of long offsets—there are low frequencies, but only reflections illuminate the target zone. We first consider the lack of large offsets in the data, presuming that frequencies, low enough to cover the chosen wavenumber by backscattering, are present in the data. This means Kz > ωmin/vp in the case of P–P scattering. The maximum offset defines the maximum opening angle available for the wavenumber illumination, and we investigate the decoupling at small scattering angles part of the reflection-based patterns. This is also the case for high wavenumber recovery in general. Lack of low frequencies—there are diving waves, but there are not enough low frequencies. Next, we assume that there are not enough low frequencies in the data and, therefore, $$\cos \frac{\theta _{{\rm max}}}{2} \omega _{{\rm min}}/v_p < K_z < \omega _{{\rm min}}/v_p$$. In this case, we also assume that there are diving waves that illuminate the perturbation, that is, θmax = π. For each of the three cases specified above, we consider arbitrary azimuth coverage. Instead of comparing the spectral sensitivity patterns visually, we analyse them quantitatively using the SVD analysis. Workflow for the perfect illumination case: Do the SVD of the $$\mathbb {A}$$ matrix (each row represents a particular scattering experiment). Data at all angles in radiation patterns are inverted together. The number of parameters that are invertible is equal to the number of non-zero singular values in SVD. The principal null space is the span of singular vectors related to zero singular values. The graph of the singular values tells us how many singular values are above a certain threshold, which defines the number of invertible parameters in the presence of noise. Noise defines the threshold for realistic cases. The threshold depends on the signal-to-noise ratio in the gradient of a misfit functional (e.g. Artman & Witten 2010), not in the data. We set the threshold for the resolution matrix at 0.1 of the highest singular value (this corresponds to the eigenvalues of the Hessian above 0.01) and create the resolution matrix. The elements with high values at the diagonal of the resolution matrix are candidates for the possibly invertible parameters. Finally, we show the resolution matrix for inversion in the pseudo-acoustic approximation. 5.1 Perfect illumination 5.1.1 Cij, ρ parametrization We analyse SVD of the full matrix $$\mathbb {A}_{P-P}$$ in the Cij, ρ parametrization. Fig. 7(a) shows that the matrix has only six non-zero singular values and thus only six independent parameters are invertible in this case. Fig. 7(b) shows the set of singular vectors in the Cij, ρ parametrization. The numeric null space is the linear span of four singular vectors, corresponding to zero singular values (columns 7–10 in Fig. 7b). These vectors are:   \begin{eqnarray} 2\mathbf {C}_{12}-\mathbf {C}_{66}, \end{eqnarray} (28)  \begin{eqnarray} 2\mathbf {C}_{13}+\mathbf {C}_{55}, \end{eqnarray} (29)  \begin{eqnarray} 2\mathbf {C}_{23}+\mathbf {C}_{44}, \end{eqnarray} (30)  \begin{eqnarray} \mathbf {C}_{11}+2\mathbf {C}_{12}+\mathbf {C}_{22}-\mathbf {C}_{33}-\pmb {\rho }. \end{eqnarray} (31)Each Cij represents a unit perturbation vector of the respective stiffness parameters, similarly $$\pmb {\rho }$$ is unit density perturbation. In Appendix A, these combinations, which we found empirically, are substituted into the reflection-based radiation pattern expressions showing that inversion is absolutely not sensitive to them. We also find a particular combination of the stiffness matrix elements and density that does not induce reflections. Thus, density can be effectively represented by anisotropic parameters within the first-order scattering approximation. Identifiable combinations of parameters can be divided into two groups—the first three singular vectors show slightly higher sensitivity than the other three, which is in agreement with de Hoop et al. (1999). The most interesting fact is that, despite previously published works (de Hoop et al.1999; Oh & Alkhalifah 2016), in the case of perfect illumination, the gap in this standard parametrization is within one order of magnitude and all six parameters can principally be determined, granted good quality data. Finally, we set the threshold for the resolution matrix to 0.1 of the highest singular value (this corresponds to the eigenvalues of the Hessian above 0.01) and create the resolution matrix for this case. The elements with the highest values on the diagonal are candidates for a possible inversion. Off-diagonal elements indicate cross-talk between the parameters. Fig. 7(c) shows the highest values for the Cii elements and the density, which are interesting to compare with pseudo-acoustic approximation, which suggests that C44, C55 and C66 have little influence on P-wave scattering. Figure 7. View largeDownload slide (a) Singular values and (b) vectors of the matrix $$\mathbb {A}$$ for P–P scattering and Cij, ρ parametrization. Perfect illumination—maximum opening angle is equal to π, all azimuths are available. There are six non-zero singular values. Later on for resolution matrices two thresholds are used (0,3—black and 0,1—green), corresponding to noisy and clean data respectively. (c) Resolution matrix in Cij, ρ parametrization. Figure 7. View largeDownload slide (a) Singular values and (b) vectors of the matrix $$\mathbb {A}$$ for P–P scattering and Cij, ρ parametrization. Perfect illumination—maximum opening angle is equal to π, all azimuths are available. There are six non-zero singular values. Later on for resolution matrices two thresholds are used (0,3—black and 0,1—green), corresponding to noisy and clean data respectively. (c) Resolution matrix in Cij, ρ parametrization. 5.1.2 Hierarchical parametrization The radiation patterns are significantly simplified in the hierarchical parametrization and more features are distinguishable. The hierarchical parametrization is a modification to that given by Tsvankin (1997), which was shown to primary control the traveltime, so it would be natural to expect inversion to be focused on them. However, according to the SVD (Fig. 8), there are other parameters involved and the null space is different from that of the traveltime inversion, even though it has the same dimensionality. If the density is dropped out, we see (Fig. 8c) that there are four evidently distinguishable elements of the matrix; those are Vp, Vs, ε1 and εd. These four parameters can be inverted. With the simplification of inversion comes the price of losing two vectors that cannot be inverted with the assumed signal-to-noise ratio. Overall, we conclude that the inversion with perfect illumination is possible for six parameters, unless the parametrization is chosen so that the inversion is more sensitive to a subset of anisotropic parameters (e.g. Oh & Alkhalifah 2016) for realistic acquisition scenarios. Figure 8. View largeDownload slide Same as Fig. 7, but for the hierarchical parametrization. Figure 8. View largeDownload slide Same as Fig. 7, but for the hierarchical parametrization. In the next subsection, we try to eliminate the problem of null space by restricting updates to the pseudo-acoustic parameters set. Pseudo-acoustic approximation:   The inversion is sensitive to density perturbations (Fig. 7b). So, let us fix density as well as C44, C55 and C66, removing them from the inversion in order to analyse whether we can still invert for six parameters with a parameter set reduced in this way. To establish some connection with the acoustic case, we remove the horizontal shear velocity and its deviation parameters from our analysis by setting their updates to zero. This approach to pseudo-acoustic modeling is well known to produce shear wave artefacts during propagation (Duveneck et al.2008), but good enough to compare P–P scattering, since only P-wave plane waves are used for the calculations. The P–P scattering happens for perturbations of several parameters in the new parametrization, including Vs, but we deliberately keep the shear wave-related parameters (Vs, γ1, γd) and density fixed in the inversion. In the pseudo-acoustic approximation, there are still six independently determinable parameters. In the standard parametrization, we see two levels of singular values very clearly (de Hoop et al.1999), yet the levels are within one order of magnitude. All six independently determinable parameters have approximately the same impact on the data and, thus, all of them can be determined in well-illuminated areas (Fig. 9c). If the image is noisy, then only parameters corresponding to the first three singular vectors can be retrieved (Fig. 9). Figure 9. View largeDownload slide (a) Singular values and (b) vectors of the matrix $$\mathbb {A}$$ for P–P scattering. Perfect illumination—maximum opening angle is equal to π, all azimuths are available. There are six non-zero singular values. Simplified pseudo-acoustic approximation (C44, C55, C66 and ρ are not updated). This approximation creates S-wave artefacts (e.g. Duveneck et al.2008), yet allows to compare the results in our case, since we keep the plane wave in the background fixed. Figure 9. View largeDownload slide (a) Singular values and (b) vectors of the matrix $$\mathbb {A}$$ for P–P scattering. Perfect illumination—maximum opening angle is equal to π, all azimuths are available. There are six non-zero singular values. Simplified pseudo-acoustic approximation (C44, C55, C66 and ρ are not updated). This approximation creates S-wave artefacts (e.g. Duveneck et al.2008), yet allows to compare the results in our case, since we keep the plane wave in the background fixed. In the case of hierarchical parametrization, the three parameters that define an ellipsoidal medium (vp, ε1, εd) should be inverted first (Fig. 10c). The other parameters are not determinable when there is noise in the data, even if the noise level is low. When limited illumination is considered, these resolution estimates become even worse and only vp is invertible. Figure 10. View largeDownload slide Same as Fig. 9, but for the hierarchical parametrization. Figure 10. View largeDownload slide Same as Fig. 9, but for the hierarchical parametrization. Figure 11. View largeDownload slide (a) and (c) Number of resolvable parameters for different maximum opening angles and azimuth coverage in Cij, ρ parametrization, and (b) and (d)  in the new parametrization. We assume that the image-domain signal-to-noise ratio is high—threshold for singular values is set to 0.1 of the maximum singular value which is somewhat equivalent to the image-domain signal-to-noise ratio of about 100 for (a) and (b). The threshold for singular values is set to 0.3 of the maximum singular value, corresponding to image-domain signal-to-noise ratio of about 10 (c) and (d). Figure 11. View largeDownload slide (a) and (c) Number of resolvable parameters for different maximum opening angles and azimuth coverage in Cij, ρ parametrization, and (b) and (d)  in the new parametrization. We assume that the image-domain signal-to-noise ratio is high—threshold for singular values is set to 0.1 of the maximum singular value which is somewhat equivalent to the image-domain signal-to-noise ratio of about 100 for (a) and (b). The threshold for singular values is set to 0.3 of the maximum singular value, corresponding to image-domain signal-to-noise ratio of about 10 (c) and (d). 5.2 Lack of long offsets We first consider the lack of large offsets in the data, presuming that frequencies low enough to cover the wavenumber by backscattering are present. This means Kz > ωmin/vp in the case of P–P scattering. Maximum offset defines the maximum opening angle available for the wavenumber illumination in this case. Therefore, we investigate decoupling in the little scattering angles part of the reflection-based patterns. Figs 11(a) and (b) show the number of parameters that could be resolved in the Cij, ρ and the hierarchical parametrization, respectively when the signal-to-noise ratio in the image is high. Figs 11(c) and (d) show the same but when the signal-to-noise ratio is low. The azimuth is calculated from the first horizontal axis, the direction does not matter as the orthorhombic media are assumed to be symmetric with respect to x1, x3 plane. Fig. 12 shows the same as Fig. 11, but with pseudo-acoustic approximation enforced. Figure 12. View largeDownload slide Same as Fig. 11, but with pseudo-acoustic approximation enforced. Figure 12. View largeDownload slide Same as Fig. 11, but with pseudo-acoustic approximation enforced. We deduce the following observations from Figs 11 and 12: When azimuth is 0–10 deg only VTI parameters are resolved. The maximum number of those resolved is equal to three and is in agreement with (Podgornova et al.2015). At large depths when maximum opening angle is below approximately 40 deg only one parameter can be resolved. In order to resolve four parameters in the new parametrization, one needs 45 deg azimuth coverage and about 140 deg maximum opening angle and relatively clean data. More parameters can be resolved using the standard parametrization in general. When azimuth coverage is increased beyond 90 deg it does not provide new information because of the symmetry of the orthorhombic media with respect to x2, x3 plane. Basically, one can forget about multiparameter inversion when the signal-to-noise ratio is very low. When the illumination is uneven, azimuth wise, the condition number sometimes drops even though the illumination range increases. 5.3 Lack of low frequencies Here, we assume that the offset is sufficient to illuminate the target depth with diving waves, but low frequencies are either not recorded at all or contaminated by severe noise and therefore not useful for FWI. In this case, low spatial wavenumbers in the model can only be reconstructed from wide opening angles (Sirgue & Pratt 2004; Kazei et al.2013a; Alkhalifah & Plessix 2014). The minimum effective opening angle determines the active part of reflection-based radiation patterns and therefore the capability of inversion to decouple parameters for a given spatial wavenumber. Just like for the case with the lack of long offsets, we cut the singular values at two thresholds 0.1 and 0.3 (Figs 13 and 14) to represent relatively clean and noisy data acquired respectively. We observe the following effects from Figs 13 and 14: Starting from 60 deg, the number of invertible parameters (Figs 13a and b) stales with angle range increase if the data is relatively clean. From transmissions alone (180 deg scattering data), we can determine up to three parameters in the classic parametrization, and two parameters (Vp, εd) in the new hierarchical parametrization. This is, likely, due to the fact that ηd parameter, which principally is also determinable, has very low sensitivity. In case of low signal-to-noise ratio in the data, we observe some counterintuitive effects. With increase in azimuth in observations the number of determinable parameters can decrease from three to two and from four to three (Fig. 13c), and even from two to one (Fig. 13d). This happens if the problem is solved in the most straightforward way. For example, if one applies the singular values thresholding based on condition number limitations for reduced problems. The maximum number of parameters is invertible when azimuths are limited to 90 and 180 deg, similarly to the previous section. Essentially, it means that including more data into a data misfit under optimization in inverse problem can lead to less parameters inverted instead of more if the illumination is uneven. Generally, the number of invertible parameters in our simplified setup is increased when pseudo-acoustic approximation is used for Cij, ρ parametrization, and decreased for hierarchical parametrization. The latter probably happens due to deliberate exclusion of Vs from inverted parameters in the hierarchical parametrization, which balances the parametrization towards Vp sensitivity. For the Cij, ρ parametrization, the first singular vector has significant part outside of the pseudo-acoustic space (dominated by density) and therefore singular values become more balanced with density exclusion from inversion. Figure 13. View largeDownload slide The number of invertible parameters when low frequencies lack from the data, yet the whole target area is illuminated by diving waves. (a) and (c) In the standard parametrization, and (b) and (d) in the new parametrization. High S/N ratio is assumed in (a) and (b) pictures, that is, singular values below 0.1 of the highest are low-cut. Low signal-to-noise ratio is assumed in (c) and (d). Figure 13. View largeDownload slide The number of invertible parameters when low frequencies lack from the data, yet the whole target area is illuminated by diving waves. (a) and (c) In the standard parametrization, and (b) and (d) in the new parametrization. High S/N ratio is assumed in (a) and (b) pictures, that is, singular values below 0.1 of the highest are low-cut. Low signal-to-noise ratio is assumed in (c) and (d). Figure 14. View largeDownload slide Same as Fig. 13, but with pseudo-acoustic approximation enforced. Figure 14. View largeDownload slide Same as Fig. 13, but with pseudo-acoustic approximation enforced. 6 DISCUSSION According to the effective angle of illumination concept, several frequencies are needed in the misfit functional in order to resolve the trade-offs between parameters. In monoparameter inversion, the lack of low frequencies can to some extent be compensated by long offsets (Sirgue & Pratt 2004; Kazei et al.2013b), and vice versa. In multiparameter inversion, the increase of the offset will lead to the blending of parameters. Thus, there is no point to increasing the offset without also increasing the bandwidth. In 3-D linearized inversion, several orthorhombic parameters can be recovered from a single frequency, which can be opposed to the bands (sets) of frequencies needed to perform multiparameter FWI in 2-D. However, the inversion of a single frequency can only be monoparameter, unless: There are several components recorded; then the number of parameters for inversion can be up to the number of components recorded. The observations are 3-D; then the information about azimuthal variations of scattering and corresponding parameters can be recovered. There are constraints on the parameter relations (e.g. constant Poisson ratio); then the multiparameter inversion can become monoparametric. Assumptions about the spectrum (structure) of the model lead to coverage of a single wavenumber from several angles at a given frequency, for example, different sparsity regularizations are implemented. Our study shows that the optimal parametrization is very much dependent on the case; moreover, it is dependent on the background model. The newly introduced parametrization shows high sensitivity to the longitudinal waves velocity, comparatively to the anisotropic coefficients. Thus, its primary use could be in building a velocity model when the anisotropy parameters are of secondary importance or in helping to invert for velocity. At the same time, if anisotropic parameters are the target of inversion, then it would be possible to fix the velocity and focus on the anisotropic parameters. For simultaneous inversion, it might be more reasonable to use a straightforward stiffness, density parametrization. In the standard parametrization, it makes no sense to invert for a subset of parameters because there are no significant gaps in between and the inversion would be unstable (Appendix C). For simplicity and generality, we consider Born-linearized inversion in an isotropic background and plane incident/scattered waves. Ichinose et al. (2000) showed that the far-field term is good enough for scatterers at distances about half of the wavelength from the source and receiver locations, which validates our analysis for marine streamer surveys of the reservoir areas. He & Plessix (2017) demonstrated that, in weak VTI backgrounds, the radiation patterns remain very similar to those in isotropic backgrounds. Therefore, we expect that the results of our analysis would not change dramatically if a background with weak anisotropy was considered. Here, we consider weak anisotropy approximation and, thus, our results compliment the numerical tests for strong background anisotropy by Köhn et al. (2015). Estimation of image-domain signal-to-noise ratio is an important tool to apply these general results, as data-domain signal-to-noise ratio cannot be used to quantify the resolution without particular acquisition set up. In order to compare parametrizations we use rescaling, which is equivalent setting vp and ρ to one as well as to that used in Podgornova et al. (2015) and Oh & Alkhalifah (2016). Our approach also showed that the condition number of the sensitivity matrix $$\mathbb {A}^T\mathbb {A}$$ heavily depends on the scaling of the parameters. In particular, normalizing partial derivatives admits SVD of $$\mathbb {A}$$ properties similar to those of Cij, ρ parametrization. It shows that even scalar scaling of parameters can change multiparameter inversion dramatically. 7 CONCLUSIONS Linearized dynamic inversion of P waves for weakly anisotropic orthorhombic medium is fundamentally sensitive to six parameters only. Surprisingly, the number of parameters is the same as for the traveltime inversion. In particular, we found that vertical variations of density show the same scattering as a certain combination of elastic constants. Pseudo-acoustic inversion (i.e. inversion for pseudo-acoustic parameters only, keeping other elastic parameters fixed) constraint allows to completely remove the fundamental ambiguity in linearized inversion of P waves for orthorhombic media in the case of perfect illumination. The coverage of different wavenumbers in the velocity is investigated through spectral sensitivities. We introduced the effective angle of illumination concept, which allowed us to investigate illumination for particular wavenumbers in the velocity model. We estimated the number of invertible anisotropic orthorhombic parameters, given the acquisition setup, by using SVD of the spectral sensitivities. It turns out that the azimuth increase beyond 60 deg does not bring significant additional information on the resolvable parameters set when inverting for orthorhombic media. For monoparameter FWI, the lack of low frequencies can be compensated by long offsets in the data (Sirgue & Pratt 2004; Kazei et al.2013a), for multiparameter inversion both need to be present to allow for parameter decoupling at low wavenumbers. In order to decouple VTI parameters reliably, the maximum frequency included into the inversion needs to be at least twice higher than the minimum one. The best resolution is provided when the target area is covered by diving waves, yet even if the illumination is limited to 120 deg opening angles full set of parameters can be recovered in the classic parametrization. The azimuthally varying parameters are more or less well resolved starting from π/3 azimuth coverage. While the hierarchical parametrization reduces the trade-offs between parameters (Oh & Alkhalifah 2016), it also reduces the number of parameters that can be inverted focusing inversion on a subset of four parameters (Vp, Vs, ε1, εd) in case elastic inversion is performed and on a subset of three parameters (Vp, ε1, εd), which is covering the case of ellipsoidal anisotropy or prior knowledge anellipticity parameters (η1, ηd). If the goal of the inversion is to obtain a velocity model that explains most of the observed data then the hierarchical parametrization is optimal. If we need to extract more information about the anisotropy of a medium, considering optimal conditions and utilization of a Hessian approximation capable of treating the trade-off, and with a good initial model (such as that constructed using the hierarchical parametrization), then the classic Cij, ρ parametrization can help recover more anisotropic features of a reservoir at later stages. Acknowledgements We thank Alexey Stovas and Daniel Köhn for their constructive and helpful reviews of the manuscript. We are grateful to Ju-Won Oh, Nabil Masmoudi and Ramzi Djebbi of SWAG, KAUST for helpful discussions as well. VK is also grateful to Oleg Ovcharenko and Veronica Tremblay of KAUST for their suggestions on improving the manuscript readability. Finally, we thank KAUST for support and its HPC resources. GLOSSARY ω — angular frequency. x — imaging point vector. xs, xg — source and geophone radius vectors, respectively. s — unit vector pointing towards the source along the source wavepath direction. g — same as s but towards the geophone. Gik(x, x΄) ≡ Gik(x, x΄, ω) — Green’s tensor in the background medium for frequency ω. K — wavenumber domain vector. u(xs, xg) — monochromatic displacement wavefield at frequency ω. cijkl, c — stiffness tensor. Cij, C — stiffness matrix in Voigt notation. ρ, λ — density and first lame parameter. δf — variational perturbation of a function f. $$\hat{f}$$ — Fourier transform of a function f. WI − WS, δm — incident and scattered wave types and the perturbation parameter type, for example WI − WS, δm = P-P, λ or WI − WS, δm = P-SH, C13. $$\mathcal {R}_{WI-WS, \delta \mathbf {m}}(\mathbf {s},\mathbf {g})$$ — radiation pattern for plane waves and parameter perturbation of type WI − WS, δm. Einstein summation is assumed throughout the paper. (sg)ij ≡ sigj – the outer (dyadic) product of tensors (vectors). : — convolution of tensors on two indexes, for example c: sg ≡ cijklskgl. We choose Fourier convention following (Hudson & Heritage 1981):   \begin{eqnarray} u(t) = \int \limits _{-\infty }^{+\infty }\exp (-i\omega t) U(\omega ) {\rm d}\omega ,\, \hat{U}(\mathbf {K}) = \int \limits _{V}\exp (-i\mathbf {K}\mathbf {x}) U (\mathbf {x}) {\rm d}\mathbf {x}.\nonumber\\ \end{eqnarray} (32) REFERENCES Alkhalifah T., 2000. An acoustic wave equation for anisotropic media, Geophysics , 65( 4), 1239– 1250. https://doi.org/10.1190/1.1444815 Google Scholar CrossRef Search ADS   Alkhalifah T., Plessix R.-E., 2014. A recipe for practical full-waveform inversion in anisotropic media: an analytical parameter resolution study, Geophysics , 79( 3), R91– R101. https://doi.org/10.1190/geo2013-0366.1 Google Scholar CrossRef Search ADS   Alkhalifah T., 2016. Full-model wavenumber inversion: An emphasis on the appropriate wavenumber continuation, Geophysics , 81( 3), R89– R98. https://doi.org/10.1190/geo2015-0537.1 Google Scholar CrossRef Search ADS   Anikiev D., Kazei V., Kashtan B., Ponomarenko A., Troyan V., Shigapov R., 2014. Methods of seismic waveform inversion, Seism. Technol. , 11( 1), 1– 32. Artman B., Witten B., 2010. Image domain signal to noise estimate , US Patent App. 13/145,328. Bergslid T.S., Raknes E.B., Arntsen B., 2015. The influence of anisotropy on elastic full-waveform inversion, in SEG Technical Program Expanded Abstracts 2015 , pp. 1425– 1429, Society of Exploration Geophysicists. Beylkin G., Burridge R., 1990. Linearized inverse scattering problems in acoustics and elasticity, Wave Motion , 12( 1), 15– 52. https://doi.org/10.1016/0165-2125(90)90017-X Google Scholar CrossRef Search ADS   Cheverda V., Kostin V., 1995. R-pseudoinverses for compact operators in hilbert spaces: existence and stability, J. Inverse Ill-Posed Probl. , 3( 2), 131– 148. https://doi.org/10.1515/jiip.1995.3.2.131 Google Scholar CrossRef Search ADS   de Hoop M.V., Spencer C., Burridge R., 1999. The resolving power of seismic amplitude data: an anisotropic inversion/migration approach, Geophysics , 64( 3), 852– 873. https://doi.org/10.1190/1.1444595 Google Scholar CrossRef Search ADS   Devaney A.J., 1984. Geophysical diffraction tomography, IEEE Trans. Geosci. Remote Sens. , GE-22( 1), 3– 13. https://doi.org/10.1109/TGRS.1984.350573 Google Scholar CrossRef Search ADS   Duveneck E., Milcik P., Bakker P.M., Perkins C., 2008. Acoustic VTI wave equations and their application for anisotropic reverse-time migration, in SEG Technical Program Expanded Abstracts 2008 , pp. 2186– 2190, Society of Exploration Geophysicists. Eaton D.W.S., Stewart R.R., 1994. Migration/inversion for transversely isotropic elastic media, Geophys. J. Int. , 119( 2), 667– 683. https://doi.org/10.1111/j.1365-246X.1994.tb00148.x Google Scholar CrossRef Search ADS   Ewald P., 1969. Introduction to the dynamical theory of x-ray diffraction, Acta Crystallogr. A: Cryst. Phys. Diffract. Theor. Gen. Crystallogr. , 25( 1), 103– 108. https://doi.org/10.1107/S0567739469000155 Google Scholar CrossRef Search ADS   Gholami Y., Brossier R., Operto S., Prieux V., Ribodetti A., Virieux J., 2013a. Which parameterization is suitable for acoustic vertical transverse isotropic full waveform inversion? Part 2: synthetic and real data case studies from valhall, Geophysics , 78( 2), R107– R124. https://doi.org/10.1190/geo2012-0203.1 Google Scholar CrossRef Search ADS   Gholami Y., Brossier R., Operto S., Ribodetti A., Virieux J., 2013b. Which parameterization is suitable for acoustic vertical transverse isotropic full waveform inversion? Part 1: sensitivity and trade-off analysis, Geophysics , 78( 2), R81– R105. https://doi.org/10.1190/geo2012-0204.1 Google Scholar CrossRef Search ADS   Grechka V., Tsvankin I., 1999. 3-d moveout velocity analysis and parameter estimation for orthorhombic media, Geophysics , 64( 3), 820– 837. https://doi.org/10.1190/1.1444593 Google Scholar CrossRef Search ADS   He W., Plessix R.-É., 2017. Analysis of different parameterisations of waveform inversion of compressional body waves in an elastic transverse isotropic earth with a vertical axis of symmetry, Geophys. Prospect. , 65( 4), 1004– 1024. https://doi.org/10.1111/1365-2478.12452 Google Scholar CrossRef Search ADS   Hudson J., Heritage J., 1981. The use of the born approximation in seismic scattering problems, Geophys. J. Int. , 66( 1), 221– 240. https://doi.org/10.1111/j.1365-246X.1981.tb05954.x Google Scholar CrossRef Search ADS   Ichinose G.A., Goldstein P., Rodgers A.J., 2000. Relative importance of near-, intermediate-and far-field displacement terms in layered earth synthetic seismograms, Bull. seism. Soc. Am. , 90( 2), 531– 536. https://doi.org/10.1785/0119990134 Google Scholar CrossRef Search ADS   Ivanov Y., Stovas A., 2016. Upscaling in orthorhombic media: Behavior of elastic parameters in heterogeneous fractured earth, Geophysics , 81( 3), C113– C126. https://doi.org/10.1190/geo2015-0392.1 Google Scholar CrossRef Search ADS   Kamath N., Tsvankin I., 2016. Elastic full-waveform inversion for vti media: Methodology and sensitivity analysis, Geophysics , 81( 2), C53– C68. https://doi.org/10.1190/geo2014-0586.1 Google Scholar CrossRef Search ADS   Kamei R., Pratt R.G., Tsuji T., 2013. On acoustic waveform tomography of wide-angle obs data?strategies for pre-conditioning and inversion, Geophys. J. Int. , 194( 2), 1250– 1280. https://doi.org/10.1093/gji/ggt165 Google Scholar CrossRef Search ADS   Kazei V.V., Alkhalifah T., 2017. On the Resolution of Inversion for Orthorhombic Anisotropy, in 79th EAGE Conference and Exhibition 2017 . Kazei V.V., Troyan V.N., Kashtan B.M., Mulder W.A., 2013a. On the role of reflections, refractions and diving waves in full-waveform inversion, Geophys. Prospect. , 61( 6), 1252– 1263. https://doi.org/10.1111/1365-2478.12064 Google Scholar CrossRef Search ADS   Kazei V.V., Kashtan B.M., Troyan V.N., Mulder W.A., 2013b. Spectral sensitivity analysis of FWI in a constant-gradient background velocity model, in 75th EAGE Conference and Exhibition Incorporating SPE EUROPEC 2013 . Kazei V.V., Kashtan B.M., Troyan V.N., Mulder W.A., 2015. FWI spectral sensitivity analysis in the presence of a free surface, in SEG Technical Program Expanded Abstracts 2015 , pp. 1415– 1419, Society of Exploration Geophysicists. Kazei V., Tessmer E., Alkhalifah T., 2016. Scattering angle-based filtering via extension in velocity, in SEG Technical Program Expanded Abstracts 2016 , pp. 1157– 1162, Society of Exploration Geophysicists. Köhn D., Hellwig O., De Nil D., Rabbel W., 2015. Waveform inversion in triclinic anisotropic media-a resolution study, Geophys. J. Int. , 201( 3), 1642– 1656. https://doi.org/10.1093/gji/ggv097 Google Scholar CrossRef Search ADS   Kurzmann A. et al.  , 2016. Seismic Applications of Full Waveform Inversion, in High Performance Computing in Science and Engineering'16 , pp. 647– 665, eds Nagel W., Kröner D., Resch M., Springer, Cham. Google Scholar CrossRef Search ADS   Masmoudi N., Alkhalifah T., 2016. A new parameterization for waveform inversion in acoustic orthorhombic media, Geophysics , 81( 4), R157– R171. https://doi.org/10.1190/geo2015-0635.1 Google Scholar CrossRef Search ADS   Menke W., 2012. Geophysical Data Analysis: Discrete Inverse Theory , Vol. 45, Academic Press. Mora P., 1989. Inversion = migration + tomography, Geophysics , 54( 12), 1575– 1586. https://doi.org/10.1190/1.1442625 Google Scholar CrossRef Search ADS   Oh J.-W., Alkhalifah T., 2016. Elastic orthorhombic anisotropic parameter inversion: an analysis of parameterization, Geophysics , 81( 6), C279– C293. https://doi.org/10.1190/geo2015-0656.1 Google Scholar CrossRef Search ADS   Ovcharenko O., Kazei V., Peter D., Alkhalifah T., 2018. Variance-based salt body reconstruction for improved full-waveform inversion, Geophysics , submitted. Owusu J.C. et al.  , 2016. Anisotropic elastic full-waveform inversion of walkaway vertical seismic profiling data from the arabian gulf, Geophys. Prospect. , 64( 1), 38– 53. https://doi.org/10.1111/1365-2478.12227 Google Scholar CrossRef Search ADS   Podgornova O., Leaney S., Liang L., 2015. Analysis of resolution limits of vti anisotropy with full waveform inversion, in SEG Technical Program Expanded Abstracts 2015 , pp. 1188– 1192, Society of Exploration Geophysicists. Raknes E.B., Arntsen B., 2014. Time-lapse full-waveform inversion of limited-offset seismic data using a local migration regularization, Geophysics , 79( 3), WA117– WA128. https://doi.org/10.1190/geo2013-0369.1 Google Scholar CrossRef Search ADS   Rüger A., 1997. P-wave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry, Geophysics , 62( 3), 713– 722. https://doi.org/10.1190/1.1444181 Google Scholar CrossRef Search ADS   Sacks P.E., 1988. The inverse problem for a weakly inhomogeneous two-dimensional acoustic medium, SIAM J. Appl. Math. , 48( 5), 1167– 1193. https://doi.org/10.1137/0148070 Google Scholar CrossRef Search ADS   Santosa F., Symes W.W., 1989. Analysis of Least-Squares Velocity Inversion , Society of Exploration Geophysicists. Google Scholar CrossRef Search ADS   Schoenberg M., Helbig K., 1997. Orthorhombic media: modeling elastic wave behavior in a vertically fractured earth, Geophysics , 62( 6), 1954– 1974. https://doi.org/10.1190/1.1444297 Google Scholar CrossRef Search ADS   Shaw R.K., Sen M.K., 2004. Born integral, stationary phase and linearized reflection coefficients in weak anisotropic media, Geophys. J. Int. , 158( 1), 225– 238. https://doi.org/10.1111/j.1365-246X.2004.02283.x Google Scholar CrossRef Search ADS   Sirgue L., Pratt R., 2004. Efficient waveform inversion and imaging: a strategy for selecting temporal frequencies, Geophysics , 69( 1), 231– 248. https://doi.org/10.1190/1.1649391 Google Scholar CrossRef Search ADS   Snieder R., 2002. Chapter 1.7.1—general theory of elastic wave scattering, in Scattering , pp. 528– 542, eds Pike R., Sabatier P., Academic Press, London. Stovas A., 2017. Kinematic parameters of pure- and converted-mode waves in elastic orthorhombic media, Geophys. Prospect. , 65( 2), 426– 452. https://doi.org/10.1111/1365-2478.12420 Google Scholar CrossRef Search ADS   Tarantola A., 1984. Inversion of seismic reflection data in the acoustic approximation, Geophysics , 49( 8), 1259– 1266. https://doi.org/10.1190/1.1441754 Google Scholar CrossRef Search ADS   Tarantola A., 1986. A strategy for nonlinear elastic inversion of seismic reflection data, Geophysics , 51( 10), 1893– 1903. https://doi.org/10.1190/1.1442046 Google Scholar CrossRef Search ADS   ten Kroode F., Bergler S., Corsten C., de Maag J.W., Strijbos F., Tijhof H., 2013. Broadband seismic data ? The importance of low frequencies, Geophysics , 78( 2), WA3– WA14. https://doi.org/10.1190/geo2012-0294.1 Google Scholar CrossRef Search ADS   Thomsen L., 1986. Weak elastic anisotropy, Geophysics , 51( 10), 1954– 1966. https://doi.org/10.1190/1.1442051 Google Scholar CrossRef Search ADS   Tikhonov A.N., Arsenin V.I., John F., 1977. Solutions of Ill-posed Problems , Vol. 14, Winston, Washington, DC. Tsvankin I., 1997. Anisotropic parameters and p-wave velocity for orthorhombic media, Geophysics , 62( 4), 1292– 1309. https://doi.org/10.1190/1.1444231 Google Scholar CrossRef Search ADS   Virieux J., Operto S., 2009. An overview of full-waveform inversion in exploration geophysics, Geophysics , 74( 6), WCC1– WCC26. https://doi.org/10.1190/1.3238367 Google Scholar CrossRef Search ADS   Wu R., Toksöz M., 1987. Diffraction tomography and multisource holography applied to seismic imaging, Geophysics , 52( 1), 11– 25. https://doi.org/10.1190/1.1442237 Google Scholar CrossRef Search ADS   APPENDIX A: NULL SPACE OF P–P INVERSION Here, we show that some radiation patterns coincide analytically in isotropic background when inversion for vertical wavenumbers is considered. Let us start with eq. (28).   \begin{eqnarray} \mathcal {R}_{2\mathbf {C}_{12}-\mathbf {C}_{66},P-P} \nonumber \\ &=& \mathcal {R}_{2\mathbf {c}_{1122}-\mathbf {c}_{1212}-\mathbf {c}_{2121},P-P}(\mathbf {s},\mathbf {g}) \nonumber \\ &=& 2 s_1s_1g_2g_2 - 2s_1s_2g_1g_2 = 0. \end{eqnarray} (A1)This happens due to the fact that for horizontal reflectors or pure vertical wavenumbers:   \begin{eqnarray} (s_1=-g_1,s_2=-g_2). \end{eqnarray} (A2)Analogously the linear combinations:   \begin{eqnarray} 2\mathbf {C}_{13}+\mathbf {C}_{55}, \end{eqnarray} (A3)  \begin{eqnarray} 2\mathbf {C}_{23}+\mathbf {C}_{44}, \end{eqnarray} (A4)do not scatter P–P waves. The other Cij combinations fall into the null space in very similar way using the connection s3 = g3. For density, the relation is a little bit trickier:   \begin{eqnarray*} \mathcal {R}_{\vec{\rho },P-P} &=& \mathbf {s}\cdot \mathbf {g}\nonumber \\ &=& s_1g_1 + s_2g_2 + s_3g_3 \nonumber \\ &=& -s_1s_1-s_2s_2+s_3s_3 \nonumber \\ &=& -s_1s_1(s_1s_1+s_2s_2+s_3s_3) \nonumber \\ &&-\, s_2s_2(s_1s_1+s_2s_2+s_3s_3) \nonumber \\ &&+\, s_3s_3(s_1s_1+s_2s_2+s_3s_3) \nonumber \\ &=& \mathcal {R}_{-\mathbf {C}_{11}-\mathbf {C}_{12}-\mathbf {C}_{13},P-P} \nonumber \\ &&+\, \mathcal {R}_{-\mathbf {C}_{21}-\mathbf {C}_{22}-\mathbf {C}_{23},P-P} \nonumber \\ &&+\, \mathcal {R}_{\mathbf {C}_{31}+\mathbf {C}_{32}+\mathbf {C}_{33},P-P} \nonumber \\ &=& \mathcal {R}_{-\mathbf {C}_{11}-2\mathbf {C}_{12}-\mathbf {C}_{22}+\mathbf {C}_{33},P-P}. \end{eqnarray*} This finishes the proof of the possibility to represent density through anisotropic parameters perturbations. The following linear combination does not scatter:   \begin{eqnarray} \mathbf {C}_{11}+2\mathbf {C}_{12}+\mathbf {C}_{22}-\mathbf {C}_{33}-\rho . \end{eqnarray} (A5) APPENDIX B: HIERARCHICAL PARAMETRIZATION Stiffness matrix coefficients can be expressed in the hierarchical parametrization as follows (Oh & Alkhalifah 2016):   \begin{eqnarray} C_{11} =\rho V^2_{p}, \end{eqnarray} (B1)  \begin{eqnarray} C_{22} =(1+2\varepsilon _d)\rho V^2_{p}, \end{eqnarray} (B2)  \begin{eqnarray} C_{33} = \frac{1}{1+2\varepsilon _1} \rho V^2_{p}, \end{eqnarray} (B3)  \begin{eqnarray} C_{12} = \rho \sqrt{ \left( V^2_{p}-(1+2\gamma _1) V^2_{s} \right)} \end{eqnarray} (B4)  \begin{eqnarray} \times \sqrt{\left((1+2\delta _3)V^2_{p}-(1+2\gamma _1) V^2_{s} \right)} \end{eqnarray} (B5)  \begin{eqnarray} - (1+2\gamma _1 )\rho V^2_{s}, \end{eqnarray} (B6)  \begin{eqnarray} C_{13} = \rho \sqrt{\left(\frac{V^2_{p}}{1+2\varepsilon _1}-V^2_{s} \right) \left(\frac{V^2_{p}}{1+2\eta _1}-V^2_{s} \right)} \end{eqnarray} (B7)  \begin{eqnarray} -\rho V^2_{s}, \end{eqnarray} (B8)  \begin{eqnarray} C_{23} = \rho \left[\left(\frac{V^2_{p}}{1+2\varepsilon _1}(1+2\gamma _d)v{s}^2 \right)\right. \end{eqnarray} (B9)  \begin{eqnarray} \times \left.\left(\frac{V^2_{p}(1+\varepsilon _d)}{(1+2\eta _1)(1+2\eta _d)}-(1+2\gamma _d)V^2_{s} \right)\right]^{1/2} \end{eqnarray} (B10)  \begin{eqnarray} -&(1+2\gamma _d)\rho V^2_{s}, \end{eqnarray} (B11)  \begin{eqnarray} C_{44} =(1+2\gamma _d)\rho V^2_{s}, \end{eqnarray} (B12)  \begin{eqnarray} C_{55} =\rho V^2_{s}, \end{eqnarray} (B13)  \begin{eqnarray} C_{66} =(1+2\gamma _1)\rho V^2_{s}. \end{eqnarray} (B14)Which results in the following expression for the Jacobian $$\frac{\partial \mathbf {C}}{\partial \mathbf {m}}$$, where ϰ = Vs/Vp:   \begin{eqnarray} \delta \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}\rho \\ C_{11} \\ C_{22} \\ C_{33} \\ C_{12} \\ C_{13} \\ C_{23} \\ C_{44} \\ C_{55} \\ C_{66} \end{array}\right) = \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}1&0&0&0&0&0&0&0&0&0\\ 1&2&0&0&0&0&0&0&0&0\\ 1&2&0&0&2&0&0&0&0&0\\ 1&2&0&-2&0&0&0&0&0&0\\ 1-2\,{\it \varkappa }&2&-4\,{\it \varkappa }&0&0&0&0&1&-4\,{\it \varkappa }&0\\ 1-2\,{\it \varkappa }&2&-4\,{\it \varkappa }&-1&0&-1&0&0&0&0\\ 1-2\,{\it \varkappa }&2&-4\,{\it \varkappa }&-1&1&-1&-1&0&0&-4\,{\it \varkappa }\\ {\it \varkappa }&0&2\,{\it \varkappa }&0&0&0&0&0&0&2\,{\it \varkappa }\\ {\it \varkappa }&0&2\,{\it \varkappa }&0&0&0&0&0&0&0\\ {\it \varkappa }&0&2\,{\it \varkappa }&0&0&0&0&0&2\,{\it \varkappa }&0\end{array}\right) \delta \left(\begin{array}{c@{\quad}c@{\quad}c}\rho \\ V_{p}\\ V_{s}\\ \varepsilon _1\\ \varepsilon _d\\ \eta _1\\ \eta _d\\ \delta _3\\ \gamma _1\\ \gamma _d \end{array}\right) \end{eqnarray} (B15) APPENDIX C: ON STABILITY OF INVERSION BY SVD TRUNCATION In exploration geophysics and in seismic inversions in particular, we often face the problem of data not completely covering the model parameters therefore prior information through regularization is required (Tikhonov et al.1977). One of the possibilities to regularize inversion is through chopping off the singular values for generalized inverse (Santosa & Symes 1989; Sacks 1988). The question that naturally arises is ‘how many singular values should be kept?’. If one keeps r singular values the inverse is called r-pseudo-inverse and the solution obtained this way is called the r-solution (Cheverda & Kostin 1995). Here, we will refer to particular results from Cheverda & Kostin (1995) on the stability of r-pseudo-inverses to justify the choice of the thresholds. This is especially true in case of multiparameter inversions. The linearized system of equation:   \begin{eqnarray} \mathbb {A}\mathbf {m}= \mathbf {d}, \end{eqnarray} (C1)where m is the parameter perturbation, d is the observed data, does not have unique solution as there are parameter perturbations that do not scatter. Regularizing the system—introducing the prior information can, however provide us with some solutions. In this case, we look at the singular value decomposition (SVD) based pseudo-inverse of the matrix $$\mathbb {A}^{-1}$$, that provides us with a solution of minimum length to eq. (C1). With singular value decomposition of matrix $$\mathbb {A}$$:   \begin{eqnarray} \mathbb {A}= \mathbb {U}\mathbb {S}\mathbb {V}^*, \end{eqnarray} (C2)its pseudo-inverse is   \begin{eqnarray} \mathbb {A}^{-1} = \mathbb {V}\mathbb {S}^{-1}\mathbb {U}^*. \end{eqnarray} (C3)The solution is defined as follows:   \begin{eqnarray} \mathbf {m}_{{\rm est}} = \mathbb {A}^{-1}\mathbf {d}. \end{eqnarray} (C4)The matrix of singular values cannot be inverted for underdetermined problems, therefore in $$\mathbb {S}^{-1}$$ there are several singular values that can be inverted. In numerical applications, some singular values would not be zero exactly and therefore it is a key point to understand which singular values need to be treated as zero. There are several strategies to decide which singular values to drop: Singular values that are lower than a certain fraction of the first singular value. This allows to have certain condition number for the pseudo-inverse, which is important for stability of the inversion. Singular values that come after a large enough gap. This ensures that the singular values/vectors would not mix in between themselves and therefore contributes to the stability estimates. Keep singular values that combine into certain part of the Hilbert–Schmidt norm of $$\mathbb {S}$$. The latter makes sense for images storages—certain part of the energy contained in the image itself is than captured and the rest is dropped, which does not work for the operators. If the data d and the operator $$\mathbb {A}$$ are perturbed by δd and $$\delta \mathbb {A}$$, respectively, and, relative perturbations are bounded by ε and δ, respectively:   \begin{eqnarray} \frac{||\delta \mathbb {A}||}{||\mathbb {A}||} < \delta , \, \frac{||\delta d||}{||d||} < \varepsilon . \end{eqnarray} (C5)If also the relative gap in the singular values si of operator $$\mathbb {A}$$ is sufficiently large   \begin{eqnarray} \frac{|s_1|}{|s_r| - |s_{r+1}|} = d_r < 1/2\delta , \end{eqnarray} (C6)then the r-pseudo-inverse $$(\mathbb {A}+\delta \mathbb {A})^{-1}$$ exists (Cheverda & Kostin 1995) for compact operators (not only matrices). Moreover Cheverda & Kostin (1995) show that the relative perturbation δmest of the r-solution is bounded in the following way:   \begin{eqnarray} \frac{||\delta \mathbf {m}_{{\rm est}}||}{||\mathbf {m}_{{\rm est}}||} \le \mu _r(\mathbb {A}) \frac{(\rho +\tau )(1+\beta \delta )}{1-\mu _r(\mathbb {A})\rho }, \,{\rm if}\, \mu _r \rho < 1. \end{eqnarray} (C7)The following quantities are used in eq. (C7): the inconsistency parameter θ   \begin{eqnarray} \theta = \frac{||\mathbb {A}\mathbf {m}_{{\rm est}} - \mathbf {d}||}{s_r ||\mathbf {m}_{{\rm est}}||}, \end{eqnarray} (C8)the condition number of the system:   \begin{eqnarray} \mu _r = \frac{s_1}{s_r}, \end{eqnarray} (C9)the other quantities in eq. (C7)   \begin{eqnarray} \beta = \frac{\sqrt{2}d_r}{\sqrt{1-d_r \delta }\sqrt{1-d_r\delta \sqrt{1-2d_r\delta }}}, \end{eqnarray} (C10)  \begin{eqnarray} \rho = \delta (1+2\beta ), \end{eqnarray} (C11)  \begin{eqnarray} \tau = \sqrt{1+\theta ^2} \left(\varepsilon + \frac{d_r\delta }{1-d_r\delta } + \beta \delta \right) \end{eqnarray} (C12)The Cheverda–Kostin’s theorem, that we formulated above, essentially proposes conditions that are sufficient to have stable inversion for r parameters. Here, we just try to understand what is the limitation on singular values caused by the requirement:   \begin{eqnarray} \mu _r \rho < 1. \end{eqnarray} (C13)The latter is equivalent to   \begin{eqnarray} \frac{s_1}{s_r} \delta (1+2\beta ) < 1, \end{eqnarray} (C14)which poses the allowable mistakes in the estimates of the operator $$\mathbb {A}$$. We will strengthen the requirements by substituting β with   \begin{eqnarray} \beta _{up} \equiv 2 \sqrt{2} d_r > \beta , \end{eqnarray} (C15)which leads to bound,   \begin{eqnarray} \delta < \frac{s_r}{s_1(1+ 4\sqrt{2} d_r) }, \end{eqnarray} (C16)therefore it is enough to have   \begin{eqnarray} \delta < \frac{1}{7\mu _r d_r}. \end{eqnarray} (C17)Eq. (C17) suggests that it is sufficient to start very close to the real model in case of multiparameter inversions, to be sure that there is no ambiguity in results of multiparameter inversion. Unfortunately, condition (C17) is hard to verify in reality. If we consider linearized inversion for a single parameter though μr = 1, the stability of inversion is governed by the gap in singular values only. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial