An analysis of the first-arrival times picked on the DSS and wide-angle seismic section recorded in Italy since 1968

We performed an analysis of refraction data recorded in Italy since 1968 in the frame of the numerous deep seismic sounding and wide-angle reflection/refraction projects. The aims of this study are to construct a parametric database including the recording geometric information relative to each profile, the phase pickings and the results of some kinematic analyses performed on the data, and to define a reference 1D velocity model for the Italian territory from all the available refraction data. As concerns the first goal, for each seismic section we picked the P-wave first-arrival-times, evaluated the uncertainties of the arrival-times pickings and determined from each travel time-offset curve the 1D velocity model. The study was performed on 419 seismic sections. Picking was carried out manually by an algorithm which includes the computation of three picking functions and the picking-error estimation. For each of the travel time-offset curves a 1D velocity model has been calculated. Actually, the 1D velocity-depth functions were estimated in three different ways which assume: a constant velocitygradient model, a varying velocity-gradient model and a layered model. As regards the second objective of this work, a mean 1D velocity model for the Italian crust was defined and compared with those used for earthquake hypocentre locations and seismic tomographic studies by different institutions operating in the Italian area, to assess the significance of the model obtained. This model can be used in future works as input for a next joint tomographic inversion of active and passive seismic data. Mailing address: Dr. Luciana De Luca, Dipartimento di Chimica e Fisica della Terra (CFTA), Università degli Studi di Palermo, Via Archirafi 26, 90123 Palermo, Italy; email: luciana.deluca@idpa.cnr.it


Introduction
In 2000 the three-year research project «Probable earthquakes in Italy between 2000 and 2030: determining priorities in seismic risk mitigation» was funded by Gruppo Nazionale Difesa dai Terremoti (GNDT).In the frame of such a project, we analysed the Deep Seismic Sounding (DSS) and Wide-Angle Reflection/Refraction (WARR) seismic sections recorded from 1968 to date in the Italian territory with the purpose of optimising local 2D velocity models and an average 3D crustal model for the whole area after integrating active seismic data with earthquake data.
An accurate determination of hypocenters and velocity model parameters requires a realistic reference velocity model.When one has little or no knowledge of the structure a common procedure tested by Ellsworth et al. (1991) and Kissling et al. (1994) is to begin the inversion using a best-fitting one-dimensional model and move progressively to more finely gridded threedimensional models.On the other hand, in wide regions with strong lateral variations of velocities large errors on earthquakes location can be introduced by the use of simple 1D layered mod-els.From active seismic experiments refined 1D and 2D velocity models can be inferred, which, being constructed independently of earthquake locations and earthquakes travel times, can provide an important constraint, in particular for the shallower portion of the crust, to obtain more accurate hypocenter locations than those produced using models based solely on earthquake data.In this regard, the first part of the analysis carried out on the set of refraction data available consisted in picking the P-wave first arrival times, calculating the 1D velocity model relative to each profile and defining an average 1D velocity model which could represent a reference model for Italian earthquake hypocenters location.
All the parameters derived from such an analysis, together with all the information connected with each profile, have been inserted into a parametric database, which is the first one dedicated to Italian active seismic data.

The data set
The seismic-refraction data analysed in this study were recorded along the lines shown in fig. 1 during the DSS and WARR projects carried out in the Italian area and surrounding seas (table I).Table I reports the number of seismic sections recorded in each of the projects and used in the revision carried out, together with the indication of the geometric configuration used for data acquisition.The last column of table I reports the main references connected with the projects.
The data set includes both off-shore and onshore profiles in which the offset ranges from few kilometres to about four hundreds kilometres.In particular 40 of the 419 analysed lines were recorded as fans.
The data set also contains the high-density wide-angle reflection/refraction data of the Lisa and CropMare II experiments (de Franco et al., Several seismic sections relative to the other projects were recorded by analogical devices, then digitised and controlled by the former IR-RS Institute of Milan CNR (Biella et al., 1994;Osculati et al., 1995;Ferrari et al., 1998) and recently processed by the IDPA Institute (Corsi, 1999a,b).The procedure adopted for digitis- Fan Cassinis et al. (1990) ing the data was optimised so as to obtain an accurate temporal and spatial position of the seismograms and it led to significant improvements of data quality.The data of the Sicilia68 (Cassinis et al., 1969) project were available only as record sections manually reconstructed on paper by the analogical seismograms, therefore a digitising procedure including the scanning of the waveforms, the vectorialisation and the sampling of the traces was used.In this case the objective of the digitalisation was not to increase data quality, but only to preserve data from deterioration.

Database
A database containing the first-arrival times of the P-phases detected on each of the seismic sections, their uncertainties, the 1D velocitydepth functions and important information on the acquisition layout of each profile was constructed by Microsoft Access.It could prove a useful tool to obtain information on all the active seismic data recorded in Italy in the last forty years.Such a database will be available on the web site of the GNDT project (http://gndt.ingv.it/Att_scient/Attivita/progetti_GNDT.htm).
The data storing system consists in a software package implemented in Matlab and interfaced with the database, which reads the header of each seismic profile, the seismic data matrices, the files containing the picked times and the errors connected with them, the 1D velocity models calculated for each profile by three different methods and the mean-square-root of residuals relative to such modelling, and inserts them into the database.
The database consists of four tables named «profiles», «sections», «shots» and «models».In all the tables some general indications on the profiles have been inserted such as the name of the project, the name of the profile and its geometric configuration.This latter indicates if the acquisition was carried out in profile or fan layout, and in a common shot or common receiver configuration.
The geographic coordinates and the elevation of the starting and ending points of each profile are stored in the first table.
The table «sections» contains the information relative to each trace of the seismic section: offset, azimuth, geographical coordinates and elevation of the recording station or shot (for common shot or common receiver configuration respectively), first-arrival pickings, the associated uncertainty and a code which identifies the seismic phase.The last three fields are repeated for the travel-times relative to P-reflected waves picked as delayed phases.Similarly the last fields of this table should contain the S-waves travel times.However, as said before only the Pwave first arrivals has been picked so far.
The table «shots» contains all the parameters connected with the shot: latitude, longitude, elevation and timing; moreover, for off-shore experiments the sea-floor level sampled at each shot location is reported.
The table «models» lists the velocities, depths and gradients relative to each profile calculated by the three inversion methods of travel time-offset curves described in Section 4.1, together with the corresponding RMS values.

P-phase first arrivals picking algorithm
The arrival of a seismic phase produces a change in the frequency content and/or amplitude of the time series.The phase onset can be more reliably detected if a proper transformation of the seismic trace can be found, discolsing the signal features and their variability more clearly than the signal itself.The accuracy of pickings largely depends on the transformation carried out, which can be identified as a signal characteristic function.
The P-phase first-arrivals picking was carried out manually by an algorithm implemented in Matlab environment.
The procedure adopted for first-arrivals recognition is based on the analysis of the following characteristic functions: the Hilbert Transform modulus of each seismogram s i, which represents the envelope of the wave connected to the seismic phase, the characteristic function CF defined by Allen (1978) as and the signal itself (fig.2).The onset time in each seismic trace is picked independently on the three functions and the first-arrival time is calculated as their mean value.Such a procedure is not applicable to identify secondaryphases or S-wave arrivals.
In order to avoid over-or under-fitting of the data and to allow and appropriate weighting of the contribution to the final solution from relatively certain and uncertain picks when the inversion of these data are carried out, an uncertainty value was assigned to the first-arrival times picked.
Uncertainties are often assigned by inspection, taking into account the signal-to-noise-ratio and the frequency content of the signals, and likely using a constant value for each phase.Since signal-to-noise-ratio evaluation by automatic procedures is often a difficult task in case of very complex signals like those present in refraction seismic sections, it is also difficult to obtain reliable estimates of the arrival-times reading errors.For this reason, the picking algorithm applied to the data set performs the error estimate by two different methods.The first attempts to evaluate the noise level in the portions of the Hilbert Transform of the traces preceding the first-arrival time; the uncertainty is given by the time interval between the onset and the first subsequent sample where the Hilbert Transform exceeds an estimated level of noise.The second takes as error the largest difference between the onset time, which is the average of the arrival-times picked on the three functions, and each of these three arrival-times.

First-arrival times and errors
P-wave first arrivals and their uncertainties were identified on 419 DSS or wide-angle seismic sections.About 1900 of the 11 820 picked arrivals are relative to fans.The picking was not performed on extremely poor quality seismic sections or in case of a small number of seismic traces in the section.
The first-arrivals identification was carried out on the sections relative to the vertical component of motion after the application of bandpass frequency filters.
Figure 3 shows the travel times versus offset for all seismic profiles.The arrival times recog-   nised are mostly Pg, Pi and some Pn phases relative to diving waves in the upper crust, lower crust and in the upper mantle respectively.Figures 4 and 5 show the travel-times reduced using velocities 6 km/s and 8 km/s respectively: it can be observed that the branches of traveltime curves characterised by an apparent velocity typical of the Pn phase (fig.5) are those at offset greater than about 170 km.
As described in the previous paragraph an estimate of the errors connected with the picked travel-times was carried out using two methods.The mean value of the errors calculated by the method which evaluates the noise level on the Hilbert Transform is 0.025 s and was higher for the arrival times of the Pg phases (0.026 s) than for those of the Pn phases (0.014 s).The mean value of the uncertainties obtained by the statistical approach is 0.061 s and as in the previous case it is larger for the Pg arrivals than for Pn arrivals (respectively 0.103 s and 0.059 s).
Figure 6 shows the frequency distributions of the errors calculated by the two methods.

1D velocity models
For each of the travel time-offset curves correlated on the seismic sections acquired in a profile configuration, a 1D velocity model was calculated.The velocity-depth functions were estimated according to three different modelling methods which assume a constant velocity-gradient, a variable velocity-gradient with depth and a constant-velocity layered model, respectively.
The algorithms used for the 1D models calculation were implemented on Matlab platform and applied to the travel-time curves after smoothing them by cubic splines in order to make the timeoffset functions, T(x i), as much as possible matching monotonously increasing functions.

1D inversion
The first method of modelling used is a Monte Carlo-like technique.The estimated parameters were in this case the surface velocity v0 and the velocity gradient g, chosen among a set S{v0,g} of starting solutions to minimise the RMS of the residual times.According to the calculated parameters the maximum penetration depth h i is obtained as and then the velocity-depth function v(hi), where i is the index relative to the n experimental points on the travel-times curves.
The second method adopts at first an inversion of the smoothed travel time-offset functions by a modified Herglotz-Wiechert formula to derive a velocity-depth model.According to this formula, the velocity vi represents the instantaneous interval velocity at the depth hi corresponding to the turning depth of a diving wave or the refractor depth for a head wave; the velocity is obtained by These velocities are used to calculate the penetration depth hi as The velocity gradients corresponding to this function are next optimised by a non-linear least squares inversion together with the surface velocity v0 calculated by the previous method.
A third 1D modelling was carried out assuming constant-velocity layers.In particular the method is based on a dynamical inversion kernel which is automatically redesigned to simulate different correspondences between the data and the model space.Certainly, even though the number of correspondence was limited, the class of possible models is large and contains models with variable-thickness layers: the numbers of these latter can be from two to five.To obtain the optimal 1D model the algorithm searches for the inversion kernel to which corresponds the minimum of the Root-Mean-Squares (RMS) of the temporal residuals and according to a criterion of model simplicity.
Each of the travel-time curves was modelled by the three methods described both to evaluate the reliability of the solutions by comparing the models obtained and to provide different para-   meterisations of the 1D structure which can be chosen as starting models for future tomographic inversions of local and/or regional earthquakes.
The average RMS values calculated over the models obtained by each of the three methods are 0.1 s, 0.2 s and 0.1 s, respectively.Figure 7 shows some examples of the velocity-depth functions calculated for four seismic profiles.
In order to carry out regional or local tomographic analyses, average 1D velocity models relative to particular areas were determined.In these cases only constant gradient and variable gradient models were calculated.These were obtained both by inverting the average travel time-offset function relative to the particular area and by calculating the mean velocity-depth function over those relative to each profile of the same area.
Figure 8 shows an example of regional 1D velocity models (constant and variable gradient) calculated for the Southern Sicilian area and for the part of the Central Apennines crossed by profile 1 of Toscana78 project (Wigger, 1984).

Contour of 1D models
Figure 9 shows the velocity distribution in four depth ranges obtained by contouring over the Italian Peninsula the velocity functions calculated for each profile.In the contour map each 1D velocity model was assigned to the position corresponding to the centre of the profile.The contouring was carried out using a mesh-grid whose dimensions along the x-and y-directions are 8.4°a nd 8.9°, respectively and the distance among the nodes is about 0.17°in both directions.
The velocity-contour maps are relative to the constant gradient models and refer to the depth classes: 0-8 km, 8-16 km, 16-24 km and 24-32 km.They were obtained by extracting those functions for which the four depth intervals selected contain at least one observation each.The average velocities for each class were then calculated for all the extracted profiles and contoured.The contours relative to the 0-8 km and 8-16 km depth classes show similar velocity patterns clearly indicating two different domains, the Peri-Tyrrhenian domain and the Apenninic chain/foreland domain.The first domain is char-acterized by higher velocities than the second one with a percentage of velocity variation of ± 6% with respect to the mean value of the two classes.The classes 16-24 km and 24-32 km exhibit contours characterized by velocity values typical of the intermediate-lower crust and of the upper mantle, with an high velocity anomaly in the Central Italy.It must be emphasized that such contours represent an attempt at aerial extension of the 1D models, that can be affected by the artefacts caused by the sparse spatial coverage of the 1D velocity functions.

Mean 1D velocity model
The analysis of the velocity models relative to each profile allowed to define an average 1D velocity model (model A) for the Italian area based on the DSS/WARR data.In addition such a model was compared with the 1D models adopted for hypocentres location by the institutions operating in the Italian territory.
Each of the velocity-depth functions obtained characterises the crustal structures for depth intervals and with a resolution which depends on the maximum investigated offset and on the inter-distance between adjacent shots or receivers along the profiles, respectively.The following procedure was used to define the average 1D model: i) extraction of the velocitydepth functions for each profile (the modelling with a variable velocity gradient has been considered); ii) linear interpolation of the velocity-depth functions using a constant-depth step; iii) selection of the velocity-depth functions for which the maximum depth is larger than a given value (40 km has been fixed in this analysis); iv) calculation of the mean velocitydepth function and of its standard deviation by averaging the interpolated velocity-depth functions.
It has been verified that the effect on the mean velocity-depth function of the depth step used for the interpolation is negligible; a step equal to 1 km was therefore chosen.
We decided to obtain the mean 1D model performing a statistic analysis of all the calculated 1D models instead of performing an inversion of the whole travel-times, since the high dispersion connected with the data set because of the heterogeneity of the area would require to define an a priori data covariance which takes into account, besides of the picking error, also the difference in times due to local velocity variations.A future work will be dedicated to this kind of inversion.
Using the procedure described above the average function and its variability interval within ± σ were obtained (fig.10).The percentage of velocity variation connected with the model ranges between 6.5% and 13%.In particular it can be observed that the minimum variation value corresponds to a depth of about 15 km.In the depth-interval 15-40 km the variation values do not exceed 8.5%, while variations between 6.5% and 13% characterise the first 15 km of the function indicating a large variability of the velocity in the shallower part of the crust.
It is intriguing to observe that the percentage of velocity variation is comparable with the crustal velocity anomaly ranges detected by different authors who carried out tomographic studies in the Italian territory (Di Stefano et al.,1999, 2003).The velocity model was subsequently compared with some 1D models adopted to localise the Italian earthquakes.These models are: the constant-velocity layered model used for the compilation of the Italian earthquakes GNDT catalogue [model 1] (ING-GNDT, 2001); the constant-velocity layered model used for the hypocentres location in the North-Eastern Italy (Slejko et al., 1989) [model 2]; the 1D model obtained from the tomographic analysis in the North-Western Italy (Kissling et al., 1995) [model 3]; the 1D model used for the tomographic inversion of Italian earthquakes data (Di Stefano et al., 2003) [model 4]; and a constant-velocity layered model optimised for the Southern Tyrrhenian area (Capizzi et al., 2001) [model 5].
In order to quantify the differences among the models, we considered a time integral function (T j) defined as where N is the number of depth points characterising the velocity-depth function, V i and ∆zi are the velocity and thickness values of the models respectively.This quantity gives a measure of the delay feature of a model, since it represents the normal incidence travel-time between the surface and a generic depth point.
The values of the integrals Tj were calculated for each model in the depth interval 0-39 km and are reported in fig.10.Such values allow us to state that all the models considered in this analysis are comparable within the limits of variability of the average model obtained from DSS/WARR data.
Figure 11 also shows the integral functions Tj (zi) to compare the models in different range of depth.
It is important to refer in particular to the global models [model 1 and model 4].Model 1 is on average characterised by velocity values lower than those of model A (fig. 11) and in particular its integral function is the closest to that   II.In order to quantitatively compare the responses calculated for the two models, the mean-square-root of their residuals were determined with respect to the optimal smoothed travel-time curve (solid line with triangles) representing the experimental data.Such a smoothed curve was obtained by performing a cubic smoothing spline for the observed data set and optimising the smoothing parameter so as to achieve the minimum misfit with the observed data.The misfit proved to be equal to 0.69 s.The mean-square-root of residuals between the smoothed experimental data and the response respectively to model A and model 1 are 0.18 s and 0.42 s.These values indicate that model A has a significantly better agreement with the experimental data than model 1.

Conclusions
The analysis of 419 seismic refraction profiles recorded in the Italian territory since 1968 in the frame of numerous national and international projects allowed us to construct a parametric database containing very useful information relative to such profiles and to define a reference 1D velocity model for the Italian area.
The P-wave first-arrival times detected on each of the seismic sections, their uncertainties, the 1D velocity-depth functions calculated by three different modelling methods and other important information of each profile have been inserted into a parametric database.
The analysis of the velocity-depth functions relative to each profile allowed us to define an average 1D velocity model for the Italian area, which was compared with other 1D models generally adopted for hypocentres locations or tomographic studies in the Italian territory.
The model optimised from the analysis of refraction data can be considered a good estimate of a 1D model for the Italian area and can be used as input for a next joint inversion of active and passive seismic data.Such an inversion could provide a reference model for next earthquake locations on the Italian territory.

Fig. 1 .
Fig. 1.Location map of the DSS and wide-angle seismic profiles.

Fig. 2 .
Fig. 2.An example of the picking procedure applied to the seismic trace indicated by the arrow: the CF function and the modulus of the Hilbert Transform (see text for details) calculated for the selected seismic trace are shown on the left and on the right of the seismic section respectively; the dot indicates the travel time picked on the signal.Independently the picking is also carried out on the CF function and Hilbert Transform and the firstarrival time is calculated as the mean value of the three pickings.

Fig. 5 .
Fig. 5. Travel time versus offset of the P-wave first arrivals with reduction velocity 8 km/s.

Fig. 6 .
Fig. 6.Frequency distributions of the picked-travel times errors: Hilbert Transform method (upper) and statistical method (lower).

Fig. 8 .
Fig. 8. Average velocity models calculated for the Southern Sicilian area (on the top) and for the part of the Central-Apennines crossed by profile 1 of the Toscana78 project (on the bottom).The graphs on the left side are relative to the constant gradient model, those on the right side are relative to the variable gradient model.The points represent the velocity functions relative to each of the profiles; the solid line with diamonds represents the 1D model obtained inverting the average travel time-offset curve of the area; the solid line is the average of the velocity-depth functions relative to each profile.

Fig. 7 .
Fig. 7. Examples of velocity-depth functions calculated for four of the seismic sections of the database.The dotted line and the solid line with circles represent respectively the constant gradient and the variable gradient models; the solid and the dashed lines represent the constant-velocity layered model and its limits of variability.

Fig. 9 .
Fig. 9. Contour maps of mean velocities in four classes of depth, relative to constant gradient models.The stars represent the centres of the profiles which the 1D model is relative to.

Fig. 11 .
Fig. 11.Integral functions versus depth for the models shown in fig.10.In grey the range of variability of the integral function calculated for the 1D DSS/WARR model.

Fig. 10 .
Fig. 10.The average 1D velocity model obtained from DSS/WARR data is compared to 1D models used for earthquake location studies.

Fig. 12 .
Fig. 12. Experimental data set (points) and the responses calculated both for the model A (solid line) and for the model 1 (dash-dot line).The responses calculated for model A ±σ are shown by the dashed lines.The solid line with triangles is the optimal smoothed travel-time curve representing the experimental data.

Table I .
Number of the seismic sections recorded in each of the project and used in revision carried out and geometric configuration used for data acquisition.The last column lists the main references connected with the projects.

Table II .
Mean 1D model values.