Radiation pattern analyses for seismic multi-parameter inversion of HTI anisotropic media

In seismic waveform inversion, selecting an optimal multi-parameter group is a key step to derive an accurate subsurface model for characterising hydrocarbon reservoirs. There are three parameterizations for the horizontal transverse isotropic (HTI) media, and each parameterization consists of five parameters. The first parameterization (P-I) consists of two velocities and three anisotropy parameters, the second (P-II) consists of five elastic coefficients and the third (P-III) consists of five velocity parameters. The radiation patterns of these three parameterizations indicate a strong interference among five parameters. An effective inversion strategy is a two-stage scheme that first inverts for the velocities or velocity-related parameters and then inverts for all five parameters simultaneously. The inversion results clearly demonstrate that P-I is the best parameterization for seismic waveform inversion in HTI anisotropic media.


Introduction
For reservoir characterisation, seismic waveform inversion needs to extract multiple physical parameters. Because all the parameters involved in the inversion have a determinate influence on the seismic waveform data, there are potential cross-talks between the different parameters in the inversion. These unavoidable cross-talks will affect the resolution of the inverted geophysical parameters (Gholami et al. 2013a). Therefore, for multi-parameter inversion, it is necessary to select an optimal model parameterization and to design an efficient inversion strategy to solve such a sophisticated nonlinear inverse problem.
The radiation pattern analysis is an appropriate method for parameterization selection. It calculates the Fréchet kernels with respect to the parameters of the subsurface point scatterer. The amplitudes of the Fréchet kernels vary with the variation of the scattering angle. The shape of such amplitude variations is a radiation pattern, which reveals the sensitivity of a parameter in the inversion objective function (Eaton & Stewart 1994;Alkhalifah & Plessix 2014;Kamath & Tsvankin 2016;Yao et al. 2018). Ideally, the radiation patterns for different parameters should be as different as possible (Tarantola 1986). Fewer overlaps between the radiation patterns are required to suppress cross-talks in multi-parameter inversion.
Parameterization can be composed of parameters that directly describe the media, such as the P-wave velocity, or parameters that are generated by combining known parameters; for example, Kim & Min (2014) defined the generated parameters as the functions of Lamé parameters and anisotropy parameters. Tarantola (1986) compared the diffraction patterns of three parameterizations (density-impedances, density-velocities, density-Lamé parameters) for elastic isotropic media and concluded that density-impedance parameterization is an adequate choice for seismic reflection data, since the diffraction patterns of the three parameterizations are different from each other. Köhn et al. (2012) tested the aforementioned three parameterizations on the Cross-Triangle-Squares (CTS) model and the modified Marmousi-II model, and obtained wellreconstructed velocity and density models.
In this paper, we study the multi-parameter inversion for the elastic anisotropic media. The vertical transverse isotropic (VTI) media is one of the commonly studied cases. Lee et al. (2010) and Gholami et al. (2013aGholami et al. ( , 2013b investigated the elastic parameters in the VTI media, while Kim & Min (2014) introduced a parameterization with Lamé parameters and two combined parameters for the anisotropic property to replace the conventional elastic parameters. Alkhalifah & Plessix (2014) used the normal moveout (NMO) and horizontal velocities together with other anisotropy parameters, for the VTI media, to compose different parameterizations and to analyse their perturbation radiation patterns, obtaining a group of parameters that have a limited trade-off among them. Other researchers, such as Silva et al. (2016), also proposed a parameterization for VTI media, whereas Kamath & Tsvankin (2016) gave the general radiation equations for parameters in heterogeneous VTI media. Masmoudi & Alkhalifah (2016) studied the radiation patterns for orthorhombic anisotropic media, which involves more parameters than the VTI media, and they composed the known parameters to form a new parameterization and azimuthal-dependent radiation patterns.
In this paper, we study parameterization in seismic waveform inversion for horizontal transverse isotropic (HTI) media. For the HTI media, Pan et al. (2016) inverted for the elastic parameters when the symmetry axis and survey line are parallel, and included the Hessian matrix in the inversion to mitigate the cross-talks between different parameters. Gao et al. (2018) considered a general case where an angle exists between the symmetry axis and survey line and inverted for the velocities and the anisotropy parameters with given azimuth angles. In the present study, we conduct an analysis of radiation patterns and attempt to select an optimal parameterization for the multi-parameter inversion for the HTI anisotropic media. We present three parameterizations for the HTI media in section 2 and analyse the radiation patterns corresponding to these three parameterizations, respectively, in section 3. Finally, we present the inverse theory, inversion strategy and demonstrations in section 4.

Three parameterizations for the HTI media
We consider three groups of parameterizations, each of which has five parameters, to represent the HTI media. The first parameterization (P-I) includes five Thomsen parameters, P-velocity, S-velocity and three anisotropy parameters. The second parameterization (P-II) considers the elastic coefficients in the stiffness tensor that constitutes the stress-strain relation. The third parameterization (P-III) is composed of five velocity parameters. The parameters in P-II and P-III are all related to the parameters in P-I. We show their relationships in what follows.
Parameterization P-I refers to five Thomsen parameters (v P,v , v S,v ,̃( E) ,̃( E) ,̃( E) ). In this parameterization, v P,v and v S,v are Thomsen's P-velocity and the slow S-velocity, respectively, and the two subscript vs here have different meanings: 'v' in v P,v means P-wave polarisation normal to the symmetrical axis, while 'v' in v S,v means S-wave polarisation perpendicular to the plane of the cracks and parallel to the symmetrical axis in the HTI media. The rest of the parameters (̃( E) ,̃( E) ,̃( E) ) are related to the three Thomsen's anisotropy parameters ( (E) , (E) , (E) ) that describe the characteristics of the velocity anisotropy (Thomsen 1986), as̃( (1) Superscript '(E)' indicates that these parameters are the equivalence parameters. All five parameters in P-I are positively valued since, in most weak anisotropic cases, the values of ( (E) , (E) , (E) ) are in range (−0.2, 0] (Rüger 1997;Tsvankin 1997).
Parameterization P-II refers to five elastic coefficients (c 11 , c 13 , c 33 , c 44 , c 55 ). The relations between these elastic coefficients and the five parameters in P-I are where is the density. Parameterization P-III consists of (v P,v , v S,v , v P,vn , v S,vn , v S,vpa ), where subscript 'n' in v P,vn and v S,vn indicates that they are the NMO velocities and v S,vpa is the fast S-velocity. The last three parameters in P-III can be obtained through the parameters in P-I by Combining (2) and (3), we obtain the relations between the elastic coefficients in P-II and the velocity parameters in Radiation patterns corresponding to parameterization P-II (c 33 , c 55 , c 11 , c 13 , c 44 ). (c) Radiation patterns corresponding to parameterization P-III (v P,v , v S,v , v P,vn , v S,vn , v S,vpa ). P-P and P-SV radiation patterns are presented in blue and red curves, respectively. Reflection angle ranges corresponding to P-P and P-SV waves are [0°, 90°] and [0°, 33.3°], respectively.

P-III as
For convenience in the comparison stage, the inversion results of any parameterizations can be converted to each of the other parameterizations.

Radiation pattern analyses
In the radiation pattern analyses, we only compare twodimensional radiation patterns along a given azimuth for the selection of parameterizations.
The equations of radiation patterns for parameters in P-I and P-II can be found in Gao et al. (2018, equations (B11) and (B13)). The equation of radiation patterns for parameters in P-III is given here as where can be either the P-wave reflection or the SV-wave reflection and, for example, R P−P is the P-P radiation pattern. The radiation patterns R P− (c ij ) for the elastic coefficients in P-II can be found in Gao et al. (2018, equation (B11)).  Thomsen parameters (v P,v , v S,v ,̃,̃,̃). (b) Second parameterization (P-II), five elastic coefficients: (c 33 , c 55 , c 11 , c 13 , c 44 ). (c) Third parameterization (P-III), velocity parameters (v P,v , v S,v , v P,vn , v S,vn , v S,vpa  When calculating the radiation patterns of three parameterizations, we set the background model to be isotropic with a scatter in the middle. The ratio of S to P velocities is 0.55. In figure 1, refers to the reflection angle, indicating the angle between the norm line and the reflection wave. The reflection angle ranges for the reflected P and S waves are [0°, 90°] and [0°, 33.3°], respectively. For the radiation patterns of parameterization P-I (figure 1a), we can make the following observations: (1) In the radiation patterns for the P-P-wave (blue curves), there are overlaps between the P-velocity v P,v (the first panel) and the rest of the parameters (the remaining four panels).
(2) For the slow S-velocity v S,v (the second panel in figure 1a), its P-P radiation patterns are similar to the P-S radiation patterns. The two patterns are also overlapped at the middle offsets, for the angle within the range of 15-75°for P-P reflection and of 5-30°for the P-S-wave reflection.
(3) The two radiation patterns for̃( E) and the slow S-velocity v S,v seem similar, and they are overlapped at the middle offsets. (4) The radiation patterns for̃( E) and̃( E) are not obvious at the far offsets.
(5) The radiation patterns for̃( E) are overlapped with (E) at the far offsets and̃( E) at the middle offsets.
For the radiation patterns of parameterization P-II (figure 1b), we can make the following observations: (1) The P-P-wave radiation patterns for c 33 are different from the other four parameters in the near offsets (small reflection angles).
(2) The radiation patterns show overlaps between c 55 and c 11 in the far offset, between c 55 and c 13 in the middle and far offset, and between c 55 and c 44 in the middle offsets.
(3) The radiation patterns for the P-S-wave are similar for all parameters except c 11 .
For the radiation patterns of parameterization P-III (figure 1c), we can make the following observations: (1) The P-P radiation patterns for P-velocity v P,v and for the remaining four parameters are different in the near offsets.
(2) The P-P radiation pattern for the slow S-velocity v S,v is different from the P-P radiation patterns for the P-wave NMO velocity v P,vn , the S-wave NMO velocity v S,vn and the fast S-wave velocity v S,vpa in the near offsets.
(3) The P-S radiation patterns for the P-wave NMO velocity v P,vn and for the S-wave NMO velocity v S,vn are almost identical, whereas the corresponding P-P radiation patterns are different, influencing the reflections with reflection angle larger than 30 and 60°, respectively. (4) The P-P radiation pattern for the fast S-velocity v S,vpa mainly concentrates on the 30-60°of the reflection angle, and has overlaps with the P-wave NMO velocity v P,vn . Its P-S radiation pattern is similar to that for S-wave velocity v S,v , hence the overlaps are unavoidable.
These analyses indicate that there are overlaps between the radiation patterns, either for P-P or P-SV waves, for different parameters. These overlaps will lead to cross-talks between different parameters in the multi-parameter inversion. Meanwhile, the radiation patterns also show different shapes for some parameters. These shape differences are helpful for separating the parameters in the inversion. Therefore, we definitely need an appropriate parameter-inversion sequence and an effective inversion strategy for the multi-parameter inversion.

Inversion experiments
In the following inversion experiments, we consider a general case where the survey line has an azimuth angle from the symmetry axis of the HTI media. While defining the intrinsic coordinate system with the x-axis parallel to the symmetry axis, the survey coordinate system is azimuthally rotated by a constant angle. The elastic coefficients (c ij and c ′ ij ) that are defined in these two coordinate systems, respectively, can be connected through the Bond transform (Bond 1943).
In the inversion, the gradients of each parameter with respect to the objective function are calculated by following the chain rule: first obtaining the derivative for c ′ ij , and then calculating c ij and the parameters in P-I or P-III sequentially according to their relations. The detailed equations for the derivatives with respect to c ′ ij , c ij and parameters in P-I can be found in Gao et al. (2018).
The first order derivatives with respect to the P-III parameters are related to the objective function with respect to the P-II parameters by The gradients of the objective function with respect to the velocities in P-III are calculated following the chain rule.
Considering the parameter magnitude differences, each parameter is normalised in the following way: where i indicates the ith parameter in each parameterization, m i is the original parameter vector, m i min and m i max are the minimum and maximum values, respectively, of m i andm i is the normalised parameter vector. The magnitude range of each parameter to [0, 1] is converted by (7).
To balance the sensitivity difference between the parameters in each parameterization, a tuning factor is used to balance the magnitude difference among the gradients of different parameters (Wang & Pratt 1997;Gao & Wang 2016). The  tuned gradient vectors will be used for updating the normalised model. All the five parameters in each parameterization share the same step length, which is calculated by the parabolic method (Vigh et al. 2009). Finally, the parameters in the normalised model are converted to their real values by an inverse application of (7).
To demonstrate the efficiencies of inversion with different parameterizations, we use a modified SEG/EAGE overthrust model to test the feasibilities of the three parameterizations (Aminzadeh et al. 1997). In the test, the S-velocity model is calculated from the P-velocity model through an empirical formula (Castagna et al. 1985). That is, the geological structures of the P-and S-velocities are assumed to be similar. In the inversion, the density is assumed to be known and it is related to the P-velocity using Gardner's formula (Gardner et al. 1974). However, the structures of the anisotropy parameters are set to be different from the velocity models, to test the resolution of the multi-parameter inversion. The structures of different anisotropy parameters are the same though, since the anisotropy parameters are influenced by the same pattern of fractures. Figure 2a presents the true models for P-I, while figures 2b and 2c shows the true models for P-II and P-III, respectively, which are converted from the models in P-I according to their relations.
Each model is divided into 110 × 200 cells, and the cell size is 10 × 10 m 2 . There are 36 shots at a burial depth of 40 m and 96 receivers at the surface for each shot. The source is a Ricker wavelet with the main frequency at 15 Hz (Wang 2015). A wave equation that constitutes the two-dimensional three-component (2D3C) displacements defined in the survey coordinate system is solved numerically by a rotated staggered-grid finite-difference method (Saenger et al. 2000;Saenger & Bohlen 2004). The time step length and the grid size follow the numerical stability condition and the dispersion condition, respectively, for the anisotropic media .
During seismic waveform inversion, a shot-encoding technique is used to save the computation time (Castellanos et al. 2015;Krebs et al. 2009;Rao & Wang 2017;. During the shot-encoding waveform tomography, all shot gathers are encoded into one super-shot gather. Therefore, this technique needs less number of forward modelling runs and it obviously reduces the computational time of the inversion procedure. According to Wang & Rao (2006), seismic waveform inversion uses seismic data of one frequency band at a time, rather than a single frequency component, and is implemented in a multi-scale fashion. First, the inversion of lowfrequency seismic data recovers the large-scale component of the model parameters. Then, the sequence inversions of higher frequency seismic data refine the model parameters gradually. Five frequency bands are used in the test, and they are [0, 6], [5, 11], [10, 16], [15,21] and [20,26] Hz.
Meanwhile, the inversion adopts a two-stage framework. When using parameterizations P-I and P-III, the first stage inverts for only the P-and S-velocities and keeps the remaining parameters as the starting models. When using parameterization P-II, the first stage only inverts for c 33 and c 55 , which are directly related to the P-and S-wave velocities. The second stage inversion, using the inverted P-and S-velocities or c 33 and c 55 models as their starting models, updates all five parameters simultaneously. In this two-stage strategy, all model parameters in the inversion of each frequency band are the same except the maximum iterations which are 100 and 200 for the first and second stages, respectively.
For the waveform inversion, the starting models in P-I are the smoothed versions of the true models. The starting models for parameterizations P-II and P-III are converted from the starting models for parameterization P-I in the same way as their true models. Figure 3 displays the final inversion models. Figure 3a shows that the P-and S-velocity models are recovered well, as with the anisotropic̃( E) and̃( E) models. However, the anisotropic̃( E) model is not reconstructed, since it is highly sensitive to the errors in the data that are caused by the model error of the other four parameters. Figure 3b shows the reconstructed elastic coefficients in P-II. The P-wave velocity and S-wave velocity and the directly related models c 33 and c 55 are well inverted in the shallow parts, whereas the other three models are not well updated. Figure 3c presents the recovered velocity models in P-III. It shows that the P-wave velocity, S-wave velocity and fast S-wave velocity have been reconstructed close to the true velocities, while the structures of the P-and S-wave NMO velocities have been contaminated by the artefacts.
In figure 4, the inverted models that are shown in figures 3b and 3c are converted to those that are shown in figures 4b and 4c, respectively, and they are compared to the inverted models of parameterization P-I that are shown in figure 4a. The P-and S-wave velocities in figures 4b and 4c are well recovered, especially the shallow part, while the converted anisotropic models are badly inverted. There are obvious artefacts that are shown in the converted anisotropic models. However, the recovered model structures appear to be similar to the true structures. Compared to the badly inverted anisotropic̃( E) model in parameterizations P-I (figure 4a) and P-II (figure 4b), the model shapes of̃( E) in figure 4c are better.
We extracted a trace at distance 1000 m from the true models and the inverted models in parameterizations P-I, P-II and P-III (figure 5). It can be seen that the inverted P-wave and S-wave velocities in P-I (red lines) and P-III (blue lines) are closer to the true models than their counterparts in P-II (green lines). For̃( E) and̃( E) models, parameterization P-I gives better results compared to parameterizations P-II and P-III, where more artefacts can be seen. Parameterization P-I fails to recover thẽ( E) model, while P-II and P-III are able to invert the trends of thẽ( E) models, although the values are incorrect.
We also show the shot gather differences for the x-component between the observed seismic data and the synthetic data by using three parameterizations P-I, P-II and P-III as an example, as shown in figure 6. The three figures are plotted using the same scale. It can be seen that the shot gather in parameterization P-I is closer to the observed shot gather than the other two parameterizations.
To check the effectiveness of the inversions, we calculate the correlation coefficients between the true model and the inverted model using: where m true and m inv are the true and inverted models, the ranges of l and k are [i-nwind, i + nwind] and [j-nwind, j + nwind], respectively, and nwind is half of the window length. In the example described here, the half window size is set as 8.
The values of each parameter are positive, and their corresponding correlation coefficients are in the range of 0 to 1, in which 1 means that the true model and the inverted model are the same. For each parameterization, the average of the correlation coefficients of the five parameters is shown in figure 7. It can be seen that parameterization P-I can recover the shallow parts with depth 0-450 m very well, and parameterization P-II also can give acceptable results in the shallow layers except for the layers between offsets 800 and 1000 m. For deep layers, with depths greater than 600 m, near the left and right boundaries the inverted models in parameterization P-III seem to be better than the results in parameterizations P-I and P-II. However, for the complex structures with offsets between 600 and1600 m, parameterization P-I can give results that are superior to those for parameterizations P-II and P-III.

Conclusion
For multi-parameter inversion in HTI media, we have compared three parameterizations based on the 2D3C wave equations. The radiation patterns indicate that there are overlaps among the sensitivities of different physics parameters, but this is also true regarding the sensitivity differences between the model parameters. In all three parameterizations, the Pand S-velocity models, or those parameters related to P-and S-velocities, can be well reconstructed, while the anisotropy parameter models are reconstructed differently depending on the parameterizations.
The results of the correlation coefficients that are between the true models and the inverted models show that in the shallow layers, parameterization P-I gives the best results. For deep layers near the left and right boundaries, parameterization P-III can recover the models better. However, parameterization P-III is unable to recover the complex parts of the model compared to the other two parameterizations. Therefore, for the parameters in the HTI anisotropic media, parameterization P-I is the best choice for the multi-parameter waveform inversion.