Next: 38.3 NOR2DD Up: 38. Inversion of travel time Previous: 38.1 VELEST Contents Index
A simple way to estimate a crustal structure from travel time data is to use grid search. A grid search is performed by locating a set of events in a given area using many different crustal models which are varied in a systematic way and finding the model which gives the lowest average RMS for all events. HYP has an option to locate a data set for a large number of different models and then determine which model gives the lowest average RMS for the data set. This might be a useful option, particularly when a sparse data set is available. In order to use this option, an additional input parameter file h_models.par is used (see below for an example). When this file is in the working directory, HYP will switch to multiple model mode SO ONLY HAVE THIS FILE IN WORKING DIRECTORY IF MULTIPLE MODEL MODE IS INTENDED. When using this option, all events must use the same STATIONx.HYP file, otherwise the program fails. The input MUST be from a single file, NOT from the data base. THE PROGRAM MUST RUN IN NON INTERACTIVE MODE.
Difference with VELEST
VELEST makes an inversion for velocity in a given fixed layer thickness model. This is fast and give a good average velocity since all data is used together. Many different layer thicknesses must be searched (kind of grid search) and it is up to the user how to do this and often a script is made by the user. A disadvantage with VELEST is that only first arrivals can be used and other arrivals must be filtered out. First arrivals are not always easy to identify, however, it is up to the user to do this independently of VELEST. With HYP grid search, all phases are used.
Example of doing grid search
An example is shown using data from Western Norway. A data set of 127 events has been selected (available in test data set under CATALOGS/wn-grid.cat), see Figure 38.1. In order to sample the area uniformly, and not use too much data, only the 18 red events in the center are used for the grid search.
The data is in Nordic2 format so parameter NORDIC FORMAT in SEISAN.DEF must be set to 2.0.
![\begin{figure}\centerline{\includegraphics[width=0.9\linewidth]{fig/grid-search}}\end{figure}](img206.svg) |
The model used for this data is the standard model used for continental Norway and has been developed with data from S. Norway.
Vp [km/s], Depth [km] 6.2 0.0 6.6 12.0 7.1 23.0 8.05 31.0 8.25 50.0 8.5 80.0
This model has been used for more than 30 years and there is no strong indication that it is wrong. So it should be interesting to see if it can be improved or if the grid search finds the same model.
The first step is to locate the test data set with the current model. The distance weighting was set to 200-250 km to only sample the area near the events. Then a check is made of which layers are sampled. This is done with program HYP_COUNT_PHASES. The program counts the phases used for location as they appear in the print.out file. The phases counted are local phases as identified by the program and not the phases given as input. Zero output weight phases are not counted. In addition the sum of the weights used for the phases are counted. Only phases for up to 9 layers are counted. The output will give an idea of how many rays sample different layers. The program is started and no questions are asked. It is assumed that the print.out is in current directory. For the above data set of 18 events, the run is
hyp_count_phases Number of phases pg pn2 pn3 pn4 pn5 pn6 pn7 pn8 pn9 171 10 0 77 0 0 0 0 0 sg sn2 sn3 sn4 sn5 sn6 sn7 sn8 sn9 173 10 0 67 0 0 0 0 0 Sum of phase weights pg pn2 pn3 pn4 pn5 pn6 pn7 pn8 pn9 168.9 10.0 0.0 62.5 0.0 0.0 0.0 0.0 0.0 sg sn2 sn3 sn4 sn5 sn6 sn7 sn8 sn9 166.4 10.0 0.0 54.7 0.0 0.0 0.0 0.0 0.0 Output file is hyp_count_phases.out
pnx and snx are refracted phases for layer number x where the surface is layer 1 so the first refracted phase is pn2.
It is seen that there are mostly direct Pg and Sg and only a few refracted phases from layer 4 (see HYP program for explanation of phase names) which is Moho. There are no phases from deeper layers since the distance is limited to 250km so there is no reason to grid search for these layers.
Below is the h_model.par used for this test:
layer # start vp delta vp # delta start h delta h # delta start vs delta vs # delta st. vp/vs del vp/vs # delta
1 5.6 0.2 4 0.0 1.0 1 1.70 0.02 6
2 6.4 0.2 3 6.0 2.0 5
3 6.9 0.2 5 16.0 2.0 6
4 7.80 0.1 5 28.0 2.0 4
4 8.25 0.05 1 50.0 1.0 1
The first line is info only. Layer # is also only for information. For each layer, there is a start P-velocity, (start vp), increment in velocity (delta vp) and number of increments (# delta). The following inputs are then the same for layer depths. Optionally also the S-velocities can be varied independently of the P-velocity, here the option is left blank to indicate that Vp/Vs is used to get the test S-velocity. The Vp/Vs ratio can also be searched for, if left blank, the ratio is taken from the station file. Vp/Vs and S-velocity cannot be searched for at the same time. There must be an entry for each layer even if no variation is used. An example input file is given in DAT. The parameters for location not set in h_model.par like Lg velocity remain unchanged.
Ideally, the whole possible model space should be searched in order to avoid possible local minima. In practice this will take too long so the model space can first be search with a larger spacing and then refined in a smaller range. In our case a rather fine grid was used and also the Vp/Vs was varied giving a total of 216000 models. There is a check in HYP that models do not overlap. With 18 events this takes about 24h on a standard PC. The model space searched is summarized below:
Vp range Depth range Vp/Vs range 5.6-6.2 1.70-1.80 6.4-6.8 6-14 6.9-7.7 16-24 7.8-8.2 28-34
It is seen that he model space is searched quite extensively and hopefully including the range of possible models. The average RMS variation for these models (see below) is 0.3 to 1.2s so that many unlikely models are checked.
When HYP starts up, it will print out how many permutations are required so if more than a few thousand, reduce the number of models, of course depending on how many events are used. Using all parameters at the same time will often give too many combinations like searching for
Vp and Vs independently. In any case, it is an advantage to first try with just a few models to get a feeling for how sensitive the data is for model changes.
An output file h models.out is generated, see example below. For each model tested, one output line is given with the RMS and the model. In the example below only the last 5 models are shown. Since many models can have very similar average RMS, the best 30 models are printed at the end. These models are also printed out in the format used in STATION0.HYP.
Example of h_models.out with the 18 events
First is shown some lines from each test, in this case there were 216000. Then follows the best 30 models, both as one line and in format ready for the station file. The average hypocentral depth for the particular model is shown to get an idea of how much the depth depends on the model. Also shown is the range of RMS for all the test models. This gives an idea of the real variation between the models.
0.510 1.700 5.60 3.29 0.00 6.40 3.76 6.00 6.90 4.06 16.00 7.80 4.59 28.00 8.25 4.85 50.00 6.93 0.475 1.700 5.60 3.29 0.00 6.40 3.76 6.00 6.90 4.06 16.00 7.80 4.59 30.00 8.25 4.85 50.00 7.50 0.469 1.700 5.60 3.29 0.00 6.40 3.76 6.00 6.90 4.06 16.00 7.80 4.59 32.00 8.25 4.85 50.00 8.33 0.486 1.700 5.60 3.29 0.00 6.40 3.76 6.00 6.90 4.06 16.00 7.80 4.59 34.00 8.25 4.85 50.00 9.14 0.542 1.700 5.60 3.29 0.00 6.40 3.76 6.00 6.90 4.06 16.00 7.90 4.65 28.00 8.25 ... 0.584 1.800 6.20 3.44 0.00 6.80 3.78 14.00 7.70 4.28 26.00 8.10 4.50 34.00 8.25 4.58 50.00 12.56 0.639 1.800 6.20 3.44 0.00 6.80 3.78 14.00 7.70 4.28 26.00 8.20 4.56 28.00 8.25 4.58 50.00 10.70 0.619 1.800 6.20 3.44 0.00 6.80 3.78 14.00 7.70 4.28 26.00 8.20 4.56 30.00 8.25 4.58 50.00 10.99 0.605 1.800 6.20 3.44 0.00 6.80 3.78 14.00 7.70 4.28 26.00 8.20 4.56 32.00 8.25 4.58 50.00 11.49 0.590 1.800 6.20 3.44 0.00 6.80 3.78 14.00 7.70 4.28 26.00 8.20 4.56 34.00 8.25 4.58 50.00 12.13 Minimum and maximum rms 0.337170064 1.21897244 The best models rms vp/vs vp vs d vp vs d vp vs d vp vs d vp vs d av-h 0.337 1.720 5.60 3.26 0.00 6.40 3.72 6.00 6.90 4.01 26.00 7.90 4.59 34.00 8.25 4.80 50.00 14.21 0.337 1.720 5.80 3.37 0.00 6.40 3.72 6.00 6.90 4.01 26.00 7.90 4.59 34.00 8.25 4.80 50.00 14.24 0.338 1.740 6.20 3.56 0.00 6.60 3.79 10.00 6.90 3.97 26.00 8.00 4.60 32.00 8.25 4.74 50.00 11.72 0.338 1.740 6.20 3.56 0.00 6.60 3.79 10.00 6.90 3.97 24.00 8.00 4.60 32.00 8.25 4.74 50.00 11.51 0.338 1.720 5.60 3.26 0.00 6.40 3.72 6.00 6.90 4.01 26.00 8.00 4.65 34.00 8.25 4.80 50.00 13.13 0.338 1.740 6.20 3.56 0.00 6.60 3.79 10.00 7.10 4.08 26.00 8.10 4.66 34.00 8.25 4.74 50.00 11.88 0.338 1.740 6.20 3.56 0.00 6.60 3.79 10.00 6.90 3.97 26.00 8.20 4.71 34.00 8.25 4.74 50.00 11.82 0.338 1.720 5.80 3.37 0.00 6.40 3.72 6.00 6.90 4.01 26.00 8.00 4.65 34.00 8.25 4.80 50.00 13.32 0.338 1.740 6.20 3.56 0.00 6.60 3.79 10.00 7.10 4.08 26.00 8.00 4.60 32.00 8.25 4.74 50.00 11.14 .... 0.340 1.720 5.60 3.26 0.00 6.40 3.72 6.00 7.10 4.13 26.00 7.90 4.59 34.00 8.25 4.80 50.00 13.21 0.340 1.740 6.20 3.56 0.00 6.60 3.79 10.00 7.10 4.08 26.00 8.20 4.71 34.00 8.25 4.74 50.00 10.88 0.340 1.720 5.80 3.37 0.00 6.40 3.72 6.00 6.90 4.01 26.00 7.80 4.53 34.00 8.25 4.80 50.00 15.18 0.340 1.720 5.80 3.37 0.00 6.40 3.72 6.00 6.90 4.01 22.00 7.90 4.59 34.00 8.25 4.80 50.00 13.19 0.340 1.720 5.60 3.26 0.00 6.40 3.72 6.00 6.90 4.01 26.00 7.80 4.53 34.00 8.25 4.80 50.00 15.11 vp/vs= 1.720 rms= 0.337 5.60 0.00 3.26 6.40 6.00 3.72 6.90 22.00 4.01 7.90 34.00 4.59 8.25 50.00 4.80 vp/vs= 1.720 rms= 0.337 5.80 0.00 3.37 6.40 10.00 3.72 6.90 24.00 4.01 7.90 34.00 4.59 8.25 50.00 4.80
Evaluating the results
It is seen that for similar or same average RMS, there can be quite a bit of fluctuation in parameters where one parameter plays off another. In this case this is clearly seen in the velocity in first layer and the depth to the first layer so the data cannot distinguish between these models. Even if the models are similar, it is seen that the average hypocentral depths are affected by these small changes since the depth varies between 11 and 15km. So which model to choose? Program GRIDMIN calculates the average model for an RMS-range, below is an example using the above test.
gridmin
rms range to average(r) or
range determined by % increase in minimum rms(p=default)
percentage
2
Input file, def is h_models.out
test5
Minimum rms is 0.337170064
Number of iterations 216000
Number of values in rms range 171
RMS range 0.000 0.344
Vp/Vs 1.74 0.01
Vp SD Vs SD H SD
6.08 0.20 3.50 0.10 0.00 0.00
6.54 0.09 3.76 0.03 7.95 2.04
7.02 0.16 4.04 0.10 24.48 1.93
7.99 0.12 4.60 0.06 33.12 1.30
8.25 0.00 4.75 0.04 50.00 0.00
Averages for 30 best models
Vp/Vs 1.73
Vp SD Vs SD H SD
5.99 0.25 3.45 0.12 0.00 0.00
6.51 0.10 3.76 0.03 8.20 2.09
6.96 0.11 4.02 0.06 25.07 1.44
7.99 0.12 4.61 0.05 33.47 1.02
8.25 0.00 4.76 0.03 50.00 0.00
Output file is gridmin.out
So here the 171 best models (the 2and H (km) is depth to interface. The standard deviation gives an idea of the variation and uncertainty of the parameters and it is seen that the velocities are quite stable while layer thickness is less so.
In this case 18 events were used. Below is shown a comparison of results if fewer events are used:
Abbreviations: N: number of events, RMS: average rms for best model or maximum
RMS for NS number of models averaged, Vp/Vs: best of average Vp/Vs, Sd: standard deviation
Vp-x: velocity in layer x, dx: depth to layer x, P: test with P-readings only.
N RMS NS Vp/Vs Sd Vp-1 Sd Vp-2 sd d1 sd Vp-2 sd d3 sd Vp4 sd d3 sd 1 0.244 1.76 6.0 6.6 6.0 7.3 16.0 8.0 30.0 0.260 1.76 0.00 6.0 0.2 6.6 0.0 7.2 1.6 7.3 0.1 17.6 1.5 8.0 0.1 32.0 1.0 2 0.186 1.74 6.2 6.6 14.0 7.3 24.0 7.8 32.0 0.200 586 1.74 6.0 0.2 6.6 0.0 9.9 2.4 7.2 0.2 22.8 2.8 7.9 0.1 31.4 1.9 8 0.310 1.72 6.0 6.4 8.0 6.9 20.0 8.0 34.0 0.320 280 1.72 0.01 5.9 0.2 6.4 0.1 7.8 2.3 7.1 0.2 23.4 2.4 7.9 0.1 33.1 1.2 18 0.337 1.72 5.6 6.4 6.0 6.9 26.0 7.9 34.0 0.345 285 1.74 0.01 6.1 0.2 6.5 0.1 8.1 2.1 7.0 0.1 24.1 2.2 8.0 0.1 33.0 1.3 P 0.231 6.2 6.6 10.0 7.3 24.0 7.9 34.00 0.240 210 6.2 0.1 6.6 0.0 10.4 2.0 7.1 0.2 24.0 2.4 7.9 0.1 33.0 1.2 Old model 1.74 6.2 6.6 12.0 7.1 23.0 8.1 31.0 New model 1.74 6.1 6.5 8.1 7.0 24.1 8.0 33.0
It is seen that all tests find similar models but 8 and 18 events give the most similar result. The last lines show the comparison with new and old model and the largest difference is the depth to the first layer. For the 18 events, using the new model, the average RMS has decreased from 0.374 to 0.335. Using all the 127 events give a similar RMS reduction (0.435 to 0.394s) so the new model also improved the RMS for events not used for the grid search (except the 18 events), an indication that the model has been improved. The new model reduces the RMS of about 10 Since the readings have been done using the current model, there is a possibility that the readings are biased towards the current model. This is particularly the case for S-readings. A new grid search was done using only P-readings (see Table above) and the results are almost the same so it is conclude that there is not much bias. Note, there is no Vp/Vs search with P only.
Test with synthetic data
In order to test the sensitivity of the data set for different models, the grid search can be done with a test model using the theoretical arrival times for the phases picked but with arrival times
calculated with the test model. Different levels of noise can be added to the theoretical arrival times. The program used for this is TTLOCAl, also described elsewhere in the manual. The model is specified in STATIONx.HYP where x is taken for the first event in event file and station file can be in DAT or local directory. To get theoretical arrival times for data in file data.inp write
ttlocal data.inp 0.5
The 0.5 indicates the level of random noise added to the data. In this case +/-0.25s is added to the P-times and +/-0.25*Vp/Vs to the S-times. Vp/Vs is taken for the station file. The random noise is independent for P and S. All phase types and weights remain the same. Only phases P and S are modelled (P, Pg, Pb and Pn, same for S) so Lg and Rg are not changed. Testing with synthetic data gives an idea of how well the data can resolve a particular model with the given test data.