First arrival Q tomography based on an adjoint-state method

Under the assumption of invariant ray path in a weakly dissipative (high quality factor Q ) subsurface medium, a tomographic inversion approach composed of two cascading applications of first arrival traveltime and Q tomography is proposed for compensating amplitude loss caused by near-surface anomalies, such as unconsolidated soils or the overburden gas cloud. To improve the computational efficiency, these two related tomography methods were adopted with an adjoint-state technique. First, arrival traveltime tomography will be performed to provide an inverted velocity model as one of the inputs for the following first arrival Q tomography. Then, the synthetic first break generated by the inverted velocity model will be used as a stable guidance of accessing the scopes of first arrival waveforms in the time domain where the potential attenuated time information is contained. The attenuated time will be estimated through a logarithmic spectral ratio linear regression corresponding to frequency-dependent propagation responses of different wave types. All these estimated attenuated times will be applied with reference signals to generate synthetic attenuated seismic data in the time domain, and their discrepancies with real data will be evaluated using similarity coefficients. The ones with larger values will be selected as optimal attenuated time inputs for the following Q tomographic inversion. Examples of both synthetic and field data reveal the feasibility and potential of this method.


Introduction
During the propagation of seismic waves inside the Earth, they suffer amplitude attenuation and dispersion due to the inelasticity and heterogeneities of the subsurface medium (Bremaecker & Jean-Claude 1977;Müller et al. 2010). The overlying unconsolidated soils or the overburden gas clouds severely attenuate seismic energy and eventually lead to images with degraded illumination and resolution (Rice et al. 1991;Abercrombie 1997). This seismic absorption property of a medium is quantified and designated by a quality factor (Q). A proper estimation and compensation of Q for seismic data are essential operations to improve the seismic image resolution and obtain correct amplitude and phase information to make their identification and interpretation more feasible, which is a crucial ability to accurately predict reservoir properties (Best et al. 1994;Carcione et al. 2003).
The early attempts to remove these attenuation effects were performed in the data domain by using time-variant deconvolution or time-variant deconvolution (Hargreaves & Calvert 1991;Varela et al. 1993) to recover the frequency components of seismic waves. Other methodologies aiming to further stabilise such inverse Q filtering have also been developed (Wang 2002(Wang , 2006(Wang , 2008b. These methods are based on a 1D backward propagation algorithm and could be incorporated into process flow to restore the resolution in some degree; however, this has failed to handle with any real geological complexity. On the other hand, model-driven methods to mitigate attenuation effects during wave propagation adopting a wave-equation-based technique have proved to be more promising in the presence of strong Q heterogeneities and high structural complexity. A number of Q migration methods established on multi-dimensional seismic absorption compensation belong to this category (Wang & Guo 2004;Wang 2008a). Traynin et al. (2008) developed Kirchhoff Q migration by introducing amplitude and phase compensating terms. Zhang et al. (2010), Zhu et al. (2014) proposed a reverse time Q reverse time migration approach using a new formulated time-domain visco-acoustic wave equation with separated operators dominating amplitude loss and phase dispersion, respectively.
There also have been various approaches of Q tomography used to achieve promising results from either migrated or prestack seismic data. Xin et al. (2008) proposed 3D tomographic amplitude inversion based on ray theory to compensate for the amplitude loss of reflection data using prestack depth migrated common image gathers. Shen et al. (2018) proposed wave-equation-based Q tomography to handle complex wave propagations for obtaining more accurate Q factors along with both vertical and lateral direction variations. However, it is still challenging to estimate a Q model at the near surface by reflection-based Q tomographic inversion due to the fact that available offset information of shallow depth is limited (He et al. 2016). As a result, we propose a Q tomography method based on first arrival to build a near-surface Q Model for using information better from a far offset dataset.
A reliable Q factor model estimation is one essential step for successful application Q compensation imaging. The amplitude of seismic waveforms is one commonly used attribute for reconstructing a Q model. Earlier approaches, such as that of Brzostowski & McMechan (1992), used amplitude at the main frequency to invert a near-surface Q model. However, such a surface-consistent Q algorithm frequently failed in practice, since the remaining uncertainties from geometrical spreading and scattering factors normally lead to inaccurate Q model estimation.
Another prevalently used attribute is attenuated time (or dissipation time), and can be calculated according to the change of amplitude spectrum at different times or locations of seismic data. Under the assumption of a frequencyindependent Q model, the attenuated time is contained in a frequency-dependent term, which is feasible to separate from the effects of geometrical spreading and other frequencyindependent wave propagation responses in the medium.
The estimation of attenuated time is normally done by either exploiting amplitude information at different arrival times or the frequency shift caused by absorption. Ganley & Kanasewich (1980) proposed a spectrum ratio method to estimate attenuated time. Quan & Harris (1997) used the centroid-frequency shift information varied along the ray path to reconstruct the Q distribution tomographically. Hu et al. (2011) defined a frequency-weighted-exponential function to fit various asymmetric source amplitude spectrum more accurately without losing the simple closed-form relation between a CFS and Q profile. An adaptive correction corresponding to different absorption effects is implemented in the CFS method by Xin et al. (2014) to further improve the accuracy in both Q estimation and the compensation process. Similar to the CFS method, Zhang & Ulrych (2002) proposed a peak-frequency-shift (PFS) method by assuming a Ricker-like amplitude spectrum of the source signature to estimate Q factor directly from common midpoint gathers in a form of layer stripping. Combined with a PFS method, a redatuming operator was adapted to optimise the attenuated traveltime approximation by Oliveira et al. (2017), and the accuracy of the PFS method was further improved. All these attenuated traveltime estimation approaches are established on the same impulse response function describing the seismic attenuation effects. However, the influences of different responses of waveform types should also be considered. In particular, when it comes to manipulating Q tomography with first arrivals, it is crucial to identify different wave types from first arrivals blended with direct waves, turning waves and refraction for obtaining attenuated time more accurately. To address this issue, we proposed an adaptive correction method to improve the accuracy of attenuated time estimation, which we will discuss in details in this article.
The preliminary estimated attenuated traveltimes would be used for deriving both velocity and Q models by implementation into a general traveltime tomography scheme. Nowadays, for a high dense seismic acquisition survey, this typically involves over 10 000 source-receiver pair deployments, which leads to millions of traveltime picks and parameterisations from the velocity model (Vidale 1988). A conventional ray-based tomographic method has difficulty in handling with such large datasets due to extensive memory occupations and computational load resulting from the calculation of the misfit function gradient (Fréchet derivatives). Furthermore, the efficacy of a conventional ray tracing method would be compromised in a complex medium due to shadow zones and multipathing issues. Since then, numbers of traveltime tomographic inversion algorithms adopting highly stable grid-based Eikonal equation solvers and adjoint-state technique have been proposed to circumvent these issues. The early attempts (Sei & Symes 1994;Leung & Qian 2006) derived the adjoint equations of Eikonal solvers applied to transmission traveltime tomographic inversion. 578 Taillandier et al. (2009) proposed a first arrival traveltime tomography based on the adjoin-state method implemented with the fast sweeping method (FSM) proposed by Leung & Qian (2006) to solve for both forward modelling and an adjoint equation. Huang & Bellefleur (2012) further expanded the application to a joint transmission and reflection traveltime tomography combined with Huygens' Principle to bring benefits of revealing deeper subsurface structures. The original isotropic Eikonal equation solver has also been reformulated to an anisotropic one by Waheed et al. (2016) to extend the adjoint-state method solution to the anisotropic case.
Here, we describe implementation of a new nonlinear Q tomography method using attenuated traveltimes from first arrivals based on an adjoint-state technique, to maintain the efficiency and feasibility of practical application of a large seismic dataset. First, an attenuated time associated governing equation of the adjoint-state method is formulated, by which the gradient of misfit function could be calculated during Q tomographic inversion iterations. Then, to guarantee the accuracy of estimated attenuated traveltime, we investigate the characteristics of different wavefield frequency responses corresponding to refraction and other first arrival waveforms. Following that, an adaptive correction method is developed, using similarity coefficients to automatically select the optimal attenuated time. Finally, we test the algorithm on both synthetic data and field data, and obtain validated results to demonstrate the feasibility and potential of our method.

Theory
We assume the wave propagation in a dissipative homogeneous medium can be described by a linear system for a linear frequency-attenuation model (Aki & Richards 2002). If the quality factor Q is frequency-independent, then a constant Q model over the frequency band is observed in the seismic signal (Kjartansson 1979;Ganley & Kanasewich 1980). Ignoring velocity dispersion (i.e. c (f ) = c), for a particular frequency f , the corresponding wavefield solution after propagating for a certain time can be determined bŷ whereÛ ref is the wavefield frequency response at reference location, and the G factor is assumed to be frequencyindependent and responsible for the effects of geometrical spreading, reflection/transmission coefficients and so on (Červený 2001). Here, P(f ) is a frequency-dependent amplitude term of the wavefield response varied with different wave types propagating in the medium. The fourth and fifth terms in the fourth exponent are responsible for the phase shift of wavefield response, and the fifth exponent is responsible for attenuation effects on amplitude decay (Ward & Toksöz 1971). This is linearly proportional to frequency indicating that higher frequencies contents suffer more attenuation than lower ones, and this can be characterised in terms of attenuated time t * (Romero et al. 1997), defined by where the integral is evaluated along the ray path l and two points source and receiver situated on it. Q (l) and c(l) are the quality factor and seismic velocity, respectively, defined along each point of the ray path. The first arrival Q tomography requires estimation of attenuated traveltime based on equation (1). Then the attenuated traveltime will be used to invert the Q factor model according to equation (2) through tomographic inversion. Conventionally, the computation cost of ray-based tomography is almost linear with number of source-receiver pairs, whose ray path forward modelling is both time and memory consuming when coming for a large seismic dataset. To improve computation efficiency, we adopt our first arrival Q tomographic inversion with an adjoint-state technique based on the original work of Sei & Symes (1994) and Taillandier et al. (2009). The attenuated traveltime is obtained by solving for a new formulated Eikonal equation as discussed next.

Governing equation for attenuated time
Converting equation (2) from integral form to differential form, yields where d ∕ dl is the directional derivative along the direction of the ray path. Keers et al. (2012) has shown that attenuation only affects the waveform through the complex and frequency-dependent traveltime when attenuation is small (Q −1 ≪ 1), other than changing ray paths. In other words, the ray path will remain the same regardless of whether the media is viscose or not. Therefore, it is feasible for us to separate the velocity and Q tomographic inversion into two cascading applications in our proposed method. In isotropic media, the propagation direction of ray is consistent with the gradient of the first arrival traveltime field. Once we obtaine the inverted velocity model by traveltime tomography, we can calculate the first arrival by solving the Eikonal equation, which also provides the magnitude of the gradient of the first arrival traveltime field, For the convenience of derivation, denoting Q −1 (x) as m(x) to represent the parameterised attenuation model. The leftside differential term dt * ∕ dl in equation (3)  a Laplace operator applied, and multiplied it with equation (4), Equation (5) is the governing equation of attenuated time, and could be solved by the FSM algorithm, which is finite difference method performed in model parametric space. For more detailed numerical implementations of FSW, we refer to Leung & Qian (2006). In this case, the computational cost of equation (5) is linear with the number of shots rather than the number of shot-receiver pairs and there is no need to store the ray path in memory. According to equation (5), the attenuated time t * is linearised by the model parameter m(x), which implies the corresponding Q tomographic inversion is also a linear process.
To verify the accuracy of numerical solutions calculated by equation (5), we conducted a synthetic test to compare them with those of generated by conventional ray-based method based on equation (2). Both forward modelling algorithms are performed in a 2D constant background velocity model as shown in figure 1a. The model width is 2000 m and the depth 800 m. The value of rectangular attenuation anomaly (Q −1 ) is unitless with a magnitude of 20. The source located right beneath Q is an anomaly at depth of 700 m and lateral position of 1000 m. The receivers are all located along the top surface of the model. From figure 1b, we can see that two calculated attenuated times are almost identical to each other, except the small misfit above the boundary corners of the Q anomalies. These errors are caused by the limited precision of finite difference method solving for equation (5), and they are acceptable without any significant impact on the final tomographic inversion result.

Q Tomography based on adjoint-state method
Q tomography aims to obtain the estimated Q model that minimises the difference between the observed and modeled where s is shot number, r is receiver number, ns is the number of shots, nr is the number of receivers in each shot, Ω is the whole model space, w(s, r) and w(x) are weighting coefficients for each receiver, t * cal is the calculated attenuated time and t * obs is the observed attenuated time. x r is the set of all the receiver locations in one shot.
Here, we simply use the conjugate gradient method to invert for the attenuation model iteratively. The Q factor at the kth iteration will be updated by where To avoid the intensive computational burden of the Fréchet derivative matrix approximation, we adopted the adjoint-state technique to efficiently obtain the gradient g of misfit function (6) by solving for a new formulated adjoint Partial Differential Equation, and it turns the original inverse modelling problem into a forward modelling problem. Here we derive the mathematical representation of g based on an adjoint-state technique (Plessix 2006;Huang & Bellefleur 2012).
Suppose a perturbation m of the attenuation model m induces a perturbation t * of attenuated time t * , and also induces a perturbation O(m) of objective function O(m). By ignoring the higher order term of t * , O(m) can be expressed by Since equation (5) is a linear equation, m and t * are also controlled by equation (5). For expression convenience, we will only derive objective function perturbation for a single source as O s (m), and the total perturbation O(m) will be a simple summation from all sources. Introducing the relationship between m and t * into equation (9) The terms in equation (11) associated with m should be kept and the other terms can be eliminated by choosing adjoint-state variable that satisfies and the boundary condition then equation (11) can be further simplified to Now we build the relationship between m and O(m), and the gradient can be expressed as Adjoint-state variable over parametric space can be solved by equations (12) and (13)  can be explained as the back propagation of the attenuated time residual from attenuation model to the source. Once we obtained the misfit function gradient g(x), the Q factor model m(x) can be iteratively updated according to equation (7), until the misfit O(m) reduces to a given level.

Estimation of attenuated time and adaptive correction method
The attenuated traveltime is implicitly blended in the seismic waveforms, and it is difficult to pick manually from the seismogram as a conventional traveltime tomography process. As a result, we developed a feasible scheme to extract attenuated traveltime from first arrivals. In the first stage of our Q tomography method, an inverted velocity model c(x) is generated by any general traveltime tomography. It is used as one of inputs for the following Q tomographic inversion to calculate the misfit function gradient, as indicated by equation (15). Meanwhile, the synthetic first arrival traveltime for each receiver is calculated according to c(x). These are used as guidance for constraining the estimation of attenuated traveltime within certain time windows from attenuated seismic traces.
In practical terms, such implementation has the benefit of obtaining stable synthetic first arrival traveltimes for all receivers without involving manual picking errors. This is also one aspect of our proposed first arrival Q tomography that is superior to those using seismic reflection data or other types of data. According to the propagation properties described by equation (1), to avoid the complex and wrapping issues caused by the phase shifting term, we simply use its amplitude information to estimate attenuated time. From equation (1), we have the amplitude spectrum of attenuated wavefield response propagating for a certain time expressed as where A ref (f ) is the amplitude spectrum of wavefield at reference location, P g is frequency-independent geometric factor including the effects of geometrical spreading, P p (f ) is frequency-dependent propagation response for a certain wave type and e −f t * is the attenuation response. To estimate attenuated time t * , one needs to remove the effects from other frequency-dependent terms of both  A ref (f ) and P p (f ). Then the left attenuation response term e −f t * can be used for estimating t * . Now we will discuss the exact forms of P p (f ) corresponding to different wave types. Three main types of wave exist in the first arrival waveform: direct wave, turning wave and refraction wave. It is not easy to separate them directly in an apparent way. Fortunately, the direct and turning waves have the same amplitude spectra and waveforms as source wavelets, while the waveform of refraction is the integral of a source wavelet (Aki & Richards 2002). These frequency-dependent propagation responses can be expressed as for direct wave and turning wave, Figure 2 illustrates how frequency-dependent propagation responses varied with these three wave types. All these waveforms were modeled by a 3D acoustic equation finite difference solver and their waveforms of first arrival were windowed out manually. The source wavelet we used in this experiment is a Ricker wavelet at a main frequency of 15 Hz. In the left column of figure 2, we see that the direct and turning wave have an identical waveform to the original Ricker wavelet, while the refraction is different. One should note that the refraction preserves a low peak frequency even when the attenuation effect is not considered in this acoustic modelling case. Such a low peak frequency phenomenon is caused by the frequency-dependent propagation responses as indicated by equation (17).
Inspired by this observation, to identify and separate refraction from other wave types among first arrivals, we will calculate the attenuated time twice by equation (16) according to assumptions of direct/turning wave and refraction, respectively. Here we use logarithmic spectral ratio method to estimate attenuated time. The natural logarithm applied to the both sides of equation (16) gives a linear function Here, ln P g can be treated as a constant. Thus, within the selected frequency band, the attenuated time can be estimated by applying the least-square linear regression, which calculates the average slope of spectrum ratio ln These two estimated attenuated times will be applied with reference signal A ref for all selected frequency bands based on equation (16), to generate synthetic attenuated seismic data S syn in the time domain. Then we can evaluate the similarity between synthetic and real attenuated seismic data S real for each sample point (i = 0, 1, 2, ..., n.) in the time domain. The similarity coefficient is defined as Between the two estimated attenuation times, the one with the larger similarity coefficient will be selected as an optimal attenuated time for further first arrival Q tomography. Furthermore, a similarity coefficient can be served for quantifying the reliability of the attenuated time. In practice, many of the traces in the seismic data were severely contaminated with noise, making them unreliable for use in attenuated time extraction. Based on our experiment, we prefer to remove those noisy traces whose similarity coefficients are less than 585 Figure 10. Survey geometry map: source (red) and receiver (blue) positions. 0.5 and use the similarity coefficients as weighting factors in tomographic inversion for the remaining traces.

2D synthetic data test: simple model
A 2D numerical seismic survey is simulated over a 10 000 m long by 1000 m deep model comprising two layers (figure 3). The velocity of the top layer linearly increases with depth from 1500 to 3000 m s −1 . The bottom layer becomes homogeneous with a constant velocity of 4500 m s −1 . The corresponding attenuation model is shown in figure 4. The background attenuation values are zero, and two attenuation anomalies are located beneath the surface. The shapes of anomalies are both ellipses with semi-major lengths of 200 m and semi-minor axis lengths of 150 m. The attenuation value is 0.02 (Q = 50) at the center of the anomaly and reduces linearly to zero from the center. A total of 490 sources and 2000 receivers are deployed along the topographic surfaces with source and receiver intervals of 20 and 5 m, respectively. The synthetic data are generated by a visco-acoustic wave-equation finite difference solver, which is based on standard linear solid method (Carcione et al. 1988). The source wavelet is a Ricker wavelet with 15-Hz dominant frequency.
We compared the attenuated traveltime estimated by our proposed method with the theoretical attenuated traveltime generated by solving equation (5) as shown in figure 5. The almost identical result indicates the capability and validity of estimating attenuated traveltime by our proposed method in such noise-free data. In practice, we will not use all estimated attenuated traveltimes for the following Q tomography. The amplitude spectra of seismic traces will be affected by the superpositions of both turning waves and refraction. In such a case, the corresponding similarity coefficients are small enough to lead to zero weights in Q tomographic inversion.
We first use the true velocity to perform Q tomography. The inverted Q model using true velocity is shown in figure 6. Compared with the true attenuation model (figure 4), we find the inverted Q model is almost the same as the true one except for smearing 'tails' below the anomalous bodies. Such artifacts are mainly caused by the limited acquisition and non-uniform directions of the ray path near the highvelocity layer.
We then use the inverted velocity model generated by traveltime tomography (figure 7) to perform Q tomography. The inverted velocity model is structurally consistent with the true model but has a lower resolution. Velocities at the left and right boundaries are less reliable because of limited ray coverage. The attenuation model inverted by Q tomography using the inverted velocities is shown in figure 8. The shapes of the two attenuation anomalies are almost same as the one shown in figure 6 except for the appearance of insignificant 'tails' artifacts caused by the inaccurate velocity model. It is evident that the near-surface Q anomalies could be successfully reconstructed by our proposed two-step first arrival Q tomography.

2D synthetic data test: Marmousi model
The first arrival Q tomographic inversion based on an adjoint-state method is tested on the Marmousi model with strong heterogeneous presence ( acquisition geometry on the surface with 221 sources and 884 receivers at intervals 20 and 5 m, respectively. The source wavelet is same as in the previous simple 2D model experiment, whose dominant frequency is 20 Hz. The corresponding attenuation model is presented in figure 9b, and it shows relatively high attenuated trends in the left-side region compared to the right-side. As shown in figure 9d, the inverted attenuation model estimated by the true velocity model produces a good Q tomographic inversion result in the near surface, ranging from top surface down to a depth of around 300 m. It preserves a similar trend compared to the true answer from figure 9b, and even some detailed features of attenuated anomalies are slightly reconstructed, e.g. in between distances of 2400 and 3000 m. Relevantly, the result using an inverted velocity model (figure 9c) shows a similar near-surface attenuation distribution as presented in figure 9e. However, the result resolution is degraded, losing some edge and detailed feature depictions compared to figure 9d. This indicates that the accuracy of the inverted model we used for Q tomography has effects on the quality of the inverted attenuation model by our proposed method. A precise velocity model could bring benefits to both attenuated traveltime estimation and Q tomographic inversion.
With a limited offset, we could only retrieve a reasonable near-surface attenuation structure based on first arrival Q tomography. The ray coverage on the edges of model is poor and becomes worse with increasing depth, leading to unreliable inverted results. This could be one limitation of our first arrival Q tomography. Other than that, there is good agreement with the true and inverted attenuation model through the near-surface profile. This figure indicates that our tomographic method is capable of accommodating such strong heterogenous.

Field data test
The third example is an application to a 3D land field dataset acquired in a mountainous area in southern China, and the whole seismic survey area is around 350 km 2 . We selected one part of the data to test the feasibility and effectiveness of our proposed method. The acquisition geometry is shown in figure 10. The total source and trace number are 7407 and 32 558 862 respectively. The maximum offset is approximately 8500 m. The seismic recording length is 6 s with sample rate of 2 ms. The elevation map of our chosen survey area is shown in figure 11, in which the maximum elevation difference is 1040 m, including different complex landforms such as river channel deposit, piedmont alluvial fan and others. Since then, the velocity lateral variation and attenuation effect at the near surface is expected to be intensive. Additionally, the poor condition of seismic data acquisition leads to poor-quality seismic records with low SNR, which is very challenging for obtaining a satisfactory tomographic solution.
One shot gather is selected as shown in figure 12 to further illustrate the details of this dataset. The seismic traces are contaminated by noise and shallow reflections can be barely identified from the original shot gather. Division frequency scanning analysis is performed on the record to compare the low frequency (10-25Hz) and high frequency (25-60 Hz) features in figure 13 parts a and  high-frequency components in left panel is more apparent than the right one as shown in figures 14b and 14c. All these observations prove the existence of some strong attenuation anomalies in this area. More specifically, the attenuation effects in the left region are much stronger than the right region as depth increases.
A 3D near-surface velocity model is estimated by first arrival traveltime tomography as presented in figure 14. The most of recovered near-surface velocity structure of the left region shows a low velocity anomaly, while a high anomaly presents on right. According to the relationship of strata velocity and medium Q factor, a low velocity rock has less compaction and, consequently, the attenuation effect will be stronger. This observation is consistent with CMP stacking (figure 14a). The inverted velocity model is used for generating first arrivals, by which the attenuated first arrival waveforms will be estimated as described in section 3.3. The waveform length is determined by the dominant frequency from current seismic traces statistics. In this case of a poor-quality (low SNR) data condition, to avoid unstable attenuated time estimation by a single channel (trace), we employed a multi-channel weighted average function to apply to several high-quality reference traces. We choose 50 shot gathers with the largest dominant frequencies, and then apply a weighted averaging on those near offset channels within 100 m to maintain high quality of data. The averaged reference trace will be eventually used in attenuated traveltime estimation through a spectral ratio algorithm. However, the weighted averaging over multiple traces will diminish the local absorption anomaly features on a finer scale, which leads to a global 'smoothing effect' and degrades the resolution of the inverted Q model. The inverted near-surface Q model is shown in figure 16. The abrupt subsurface velocity variation at the near surface (as shown in figure 15) and limited acquisition aperture results in a shallow effective imaging depth of the inverted Q model. Below the effective imaging region, we simply fill in the same maximum Q value of each trace to enhance the visibility of the resulting image. In the inverted Q model, the low Q anomalies are mainly distributed in the left region with large thicknesses, and the high Q anomalies are in right region with small thicknesees. This is consistent with the velocity anomaly distribution (figure 15) as we discussed previously.
To further assess the validity of the inverted Q model, we perform the following procedures. First, the elevation of a high-velocity top interface is determined by the inverted velocity model. Above this interface, we assume that the outgoing seismic ray travels vertically, and hence the exact timespatial attenuation effect could be evaluated by equations (1) and (2) in CMP stacking. Eventually, these attenuation effects are removed from the original stacking, and the corrected CMP stacking is presented in figure 17a. The lateral continuity of seismic events is enhanced with respect to both amplitude and geological structural in the lateral direction. The resolution of the whole stacking image improved and the frequency bandwidths of the two panel regions broadened by recovering both low-and high-frequency contents, as shown in figures 17b and 17c.

Discussion
In this paper, we developed a new first arrival Q tomography based on an adjoint-state method. The calculation of misfit function gradient was implemented with an adjoint-state method, which turned the inverse modelling problem into a forward modelling problem. This Q tomography has high computational efficiency and its computational cost is linear with size of model parameterisation. It is feasible this can handle a large dataset problem . Furthermore, an adaptive correction method was designed to select the optimal estimated attenuated traveltime according to the frequency-dependent wave propagation responses for refraction, turning and direct wave types.
In a synthetic data test of the simple model, first arrival Q tomography recovered both shallow and deep attenuation anomalies and provided a good image resolution. Its success was partly attributed to an accurate attenuated traveltime estimation based on a good inverted velocity model. The results from the Marmousi model demonstrated some benefits and robustness in recovering near-surface attenuation anomalies in the presences of heterogeneousness. The ac-curacy of the inverted attenuation model was affected by the precision of the velocity model we used in Q tomography. Also, the capability to further improve the quality of the Q model relies significantly on how successfully we could calculate the attenuated first arrival, according to the waveforms of reference traces used in adaptive correction. Practically, apart from elaborately processing prestack data to improve the SNR of recorded traces, we implemented a statistical method for multi-channels with different weighting to produce qualified candidates as our pilot traces. The real data test results showed that this method worked effectively in low SNR circumstances and brought benefit to stabilising tomographic solutions. Even though the quality of the raw seismic record was poor and only partial offset data with limited acquisition were used, the inverted Q model at the near surface showed a reasonable consistency corresponding to the velocity model.
The resolution of the final inverted velocity and Q model are limited, since our tomographic inversion is ray-theory based and only uses the first arrival. To improve the deeper subsurface imaging, a joint inversion of both transmission and reflection could be an efficient way to enhance the final 590 inversion result. Also, tomographic inversion based on the wave equation would be another way to further improve resolution if one had enough computing resource to afford the heavy computational cost in practical applications.