Author: ggan076 (page 1 of 2)

Checker-board inversion data fitting

Our results in regards to the checker-board tests show convincing observations. However there are still work to be done in terms of the data fitting. The joint inversion misfits do not reflect the good misfit of local and teleseismic event inversion individually.

One sample

In order to see the behaviour of the joint inversion, we start with the simplest case of jointly invert one local and one teleseismic event.


Individually inverting each event shows that the misfit can be minimized. For a single source this is theoretically expected. For the local source inversion, after six iterations, the rms misfit reaches 0.46 ms. Whereas for teleseismic source, after six iterations, the rms misfit reaches 1.84 ms.


The corresponding models are as follow:

Where the first three are the initial model, the middle three are for the local event inversion, and the bottom three are for the teleseismic event inversion.


Next is to see the trend in misfit when the two events are inverted jointly. Since the events are coming from different directions, there should not be any conflicting information with regards to raypaths, so we expect the misfit to be similarly small.


Results shows that this is indeed the case, where the rms misfit after six iterations reaches 1.37 ms. This observation does not indicate a deterioration in data fit while performing a joint inversion of a single local and teleseismic source.

joint inversion – rms residual (millisecond) vs iteration


The recovered model is also improved accordingly.


Adding events

Adding one more local source to be inverted does not reduce the misfit significantly. After six iterations, rms misfit reaches 1.52 ms. On the other hand, adding another teleseismic event on top of that doubles the misfit to 3.22 ms, albeit this value is still very small.

4 events inverted


We keep adding more events. Additional local source (5 in total) produce a misfit of 3.02 ms, slightly less than 4 sources inversion. If on top of that we invert an extra teleseismic source (6 in total), the rms value goes to 4.25 ms. So far, it seems that the highest influence in terms of increasing the misfit has been in relation to the addition of teleseismic source.

6 events

Misfit trend

We would like to show the progression of the misfit of the data with the addition of events for local sources only. We show here the rms residual in millisecond after six iterations using varying number of local sources starting from 1 event increasing to 2, 5, 10, 50, and 100 events. Consistently in all cases we turn the regularization off by making eta and epsilon equal to zero.



(table) 6 iterations rms values for varying number of sources

1 2 5 10 50 100
0 302.25 232.17 245.38 231.66 250.63 261.3
1 61.99 42.75 45.35 58.04 64.75 75.62
2 13.85 10.25 12.42 18.14 25.81 32.31
3 3.76 3.38 5.5 8.79 14.77 19.17
4 1.49 1.8 3.28 5.96 10.98 13.91
5 0.8 0.77 2.25 4.39 9.61 11.42
6 0.46 0.61 1.58 3.29 8.46 9.78


Visually, we see that apart from the one source inversion, from 2 sources to 100 sources the lines corresponded to the rms residual improvement are shifted upward with the addition of more sources. This seems to imply that progressively the misfit shy away from zero with the inclusion of extra sources in the inversion.

Evaluating synthetic test misfit

One of the mystery in our synthetic tomography test is the data misfit shown by the joint inversion of local and teleseismic sources. In summary, inversion of local and teleseismic earthquakes individually shows a good data misfit for 600+ sources each, which is not reflected in the joint inversion using the same sets of sources.

In this post we perform the simplest case of joint inversion using one source each and with no regularization to examine the aforementioned observation. For the local inversion we use one event from central North Island, while the teleseismic source is an event in the Sumatra subduction zone.




Of the local source, eight iterations of inversion step shows that the root-mean-squared residual misfit decrease significantly from 212 to 0.31 ms.


rms residual (ms) vs iteration


The associated model recover a pair of checker-board cube around 37 S/175 E.




Similarly, eight iterations of the teleseismic source inversion shows that the rms residual decrease from 191 to 0.71 ms.


rms residual (ms) vs iteration


The recovered model is highly smeared especially in the EW direction. However it can still resolve patterns at shallow depth <30 km.




The ray-paths projected shows that the paths of the two sources are hardly overlapped.



Eight iterations of the joint inversion shows that the misfit is reduced as well as the two individual cases (rms residual from 202 to 0.43 ms).


rms residual (ms) vs iteration


The resulting model recover three boxes of checker-board pattern at shallow depth beneath the AVF, improving from both local and teleseismic inversion by increasing the depth penetration and lateral recovery.



This outcome tells that fundamentally the approach works as we expected. Next we test whether the addition of numerous sources and/or the usage of regularization contribute to the high misfit of the data.

AVF P-wave travel time tomography: 100+ teleseismic events

We implement P-wave travel time tomography using FMTOMO, to image the subsurface structure beneath the Auckland Volcanic Field. The input information, which is the residual travel time is obtained using an adaptive stacking code, tcas by Nick Rawlinson. Inversion is then implemented using more than one hundred teleseismic sources.

getting the residual time: tcas

The code tcas works by taking the predicted time (eg. using ak135 1-D Earth model) to a reference station to be compared with other stations. The difference between the predicted time is the base for the approximate initial alignment. Cross correlation is then performed to create the final improved alignment from which we get the residual time from the 1D-Earth model that we used.



teleseismic inversion

We apply to the tcas code 131 earthquakes limited to epicentral distance greater than 27 degrees to avoid triplication of the P-wave arrival.


100+ events


The output residuals at every stations are plotted on a map to gauge the underlying trend in velocity field. Red color on the map shows positive residual which means larger than average value thus a slow field and vice versa for blue. The size of the circle shows how big is the residual. The indication is that there is a prevalent slow region to the east of the AVF.



Inversion improve the preliminary trend by showing with more precision the split from positive to negative anomaly. The split is striking Northwest-Southeast extending to 80 km deep. We interpret this to be related to the Junction Magnetic anomaly of the Dun Mountain Ophiolite Belt.


Teleseismic inversion


Obtaining seismic travel-time residuals with tcas


In this post we have shown that the adaptive stacking tcas can reflect pickings done manually. Here we would like to present how tcas may help us in our seismic tomograpy study.

We use a catalog of three events from similar hypocentres travelling to an array of receivers.



The tcas result shows the following improvement in trace alignment. Each row shows the progression from raw data to initial alignment to the final alignment for one event.



We then plot the residual patterns at every stations that came out from tcas calculation. We superimpose pattern from all events. Red color shows positive residual which means larger than average value thus a slow field and vice versa for blue. The preliminary indication is that we have a slow region to the east of the AVF and fast to the westward.



Inversion uses a model with velocity grid cells spread roughly every 0.1 degrees (both NS and EW) and 10 km in depth. The picking error is estimated to be 0.1 s. Regularization is used with damping and smoothing terms of 1 and 2. And we have epsilon equals 1 and eta equals 10. The result shows misfit and chi^2 that are decreasing with iterations. With chi^s converge to value around 1 after two iterations.



The resulting model presents a perturbation of velocity field around the average with similar pattern to the presumption.



more events

We move forward this time by adding more earthquakes as our sources. As a qualification, we consider teleseismic events further than 27 degrees epicentrally from Auckland. Also we only consider events with good quality data by plotting and judging the trace output. We end up with 13 events.



There are reversal in polarity at station HIZ for events prior to 2014, and at OUZ before 2004. This has been addressed prior to applying tcas.

Final tcas alignment for all events shows good results. The output of its residuals/final time-shifts too are well contained between +/- 2 s. So we think tcas has performed as we expect.



The map of the residuals across the array shows more consistency. We again see slow velocity to the east of AVF.




Using the same model, inversion is done with standard deviation set to 0.1 s. If we decide to make all regularization terms equal to 0, chi^2 seems to converge around 7 after 4 iterations. This is equivalent to a variance reduction of 65%.



With regularization parameters for damping and smoothing set to 1 and 1 with epsilon of 4 and eta of 8, velocity anomaly start to show a pattern resembling the dt residual map. Chi^2 converge to 8 after 2 iterations.



station TLZ and TOZ

It seems like the conflicting information at station TLZ and TOZ prevent the fit to go below chi^2 equals to 7. Once we removed these two stations from the inversion, chi^2 converge around 5 after 4 iterations with regularization terms set to zero. The variance is reduced from 0.20 to 0.05 s^2 or equal to 75% reduction.



Similarly with the same regularization applied, chi^2 converge to 7 after 2 iterations but with more coherent resulting model.



inversion with catalog of 100+ events

We keep adding events for consideration in the inversion to get a better image of the subsurface presented by the resulting model. We would like to limit our teleseismic event to be not closer than 27 degrees epicentrally from our most Northwesterly station to triplication of arrival. The resulting catalog consist of 166 events:



We would also like for the catalog to contain only data with reasonable signal quality after lowpass filtering of 2.0 Hz has been done. On top of data availability, this cut the usable source data further but still we have 131 earthquakes which we have better confidence to work on.

After performing tcas to these events, we plot the residual dt map using the output from tcas. The map shows that the trend of slow to the east of AVF is retained. Residuals at station TLZ and TOZ are more consistent too. There are some stations near the centre where the prevalent residuals can be negative and positive. But these are fairly close to the orientation of the presumed split from fast to slow, so we think this is an understandable outcome.



The output of the inversion with zero regularization shows the following.



While by using the same regularization parameters as the previously,


Manual picks and TCAS: Comparison of residuals of three events with the same epicentre

We would like to compare the result we obtained in this post to the same set of events but by using the adaptive stacking TCAS approach.



We plot the raw data, initial alignment by removing the predicted arrival travel-time, and the final alignment after the cross-correlation. We allow the final time-shift to be within +/- 2 seconds.



The residuals contained in the output are consistent considering the same event path. Furthermore they are in relative agreement with the residuals obtained by manual picking with the mean removed from the observed minus predicted values.



Firstly the result shows that tcas is robust as it gives a consistent output for equal path. Secondly it shows that it reflects our manual work by showing similar perturbation pattern in the stations. Therefore there is the feasibility to use tcas to automate pickings using more sources.

Comparing residuals for three events with the same epicentre


We have three events in the Celebes Sea south of Philippine that have virtually the same epicenter. So presumably these would yield equal or at least similar paths to our stations.



In details, events locations are:




The ray-paths are projected to be indeed of uniform trajectory. Where we have basically a single path to each station from all three events.



arrivals and residuals

We would like to compare the residuals out of these three events i.e. the differential time of predicted ak135 model and observed arrival.

We will pick the arrival and find the residual by computing the time between that picked time (blue line) and predicted time (red line). We would expect them to be fairly close to one another for each respective station.



The results are combined in a spreadsheet:




Of course since the depth of the event is not exactly the same, this yields to a difference of two seconds for the amount of time it took to reach the stations. We think this is a reasonable difference in consideration.

The observed arrivals are all late arrivals in respect to the predicted arrivals. Value wise, the differences in residuals between events are fairly consistent.

inversion and misfit

We proceed with the inversion using this data. Regularization is turned on. It seems that the most satisfying result is obtained by making the damping in the regularization to equal zero. While we do so, smoothing need only to have a factor of 1 to have this progression in misfit.



The model correspond to this inversion steps is the following.



The negative anomaly is completely expected considering all of our arrival are late with respect to the 1D-Earth model.

large scale real data tomography: teleseismic events only

We perform the inversion using teleseismic events only for comparison. The stations and picks are identical to the data for teleseismic earthquake in this attempt.


Eighteen stations are used, divided into 11 in central Auckland, and 7 surrounding the former array. Low resolution channels are removed.


There are 121 events considered, though 13 are discarded due to bad data. We also have issue with the time window while we were picking.


The inversion shows negative velocity anomaly in most area in the Northland. The negative velocity in the North is a feature of our joint inversion.

There isn’t any change in travel-time residual distribution between initial and final model. The distribution is centred more in the positive direction which explains the mainly negative velocity anomaly.

AVF tomography: local events + northern illumination

Relative travel-time tomographic inversion is performed to the Auckland Volcanic Field (AVF), extending to include the northern parts of the North Island New Zealand. First P-wave arrivals are picked manually from earthquake traces that can be ranged in terms of quality. The travel-time residual between the picked arrival and the background ak135 1D-Earth model is the basis of this tomographic inversion.


Eighteen stations are used in the inversion. The stations can be differentiated into the 11 that are centred to the city of Auckland, and the remaining 7 that surround the former array. For all of the stations, we are not using the low resolution channels.


The events that we use here can be separated into local and teleseismic events. The local events make up the sources that are originated from the South of the stations array. Meanwhile the teleseismic events, mostly from the pacifics, cover the illumination from the North.

There are 116 events sourced locally that we looked at. Most of these events have good trace data, but there are 4 events in which picking is not applicable due to the noise. From the remaining 112, some stations don’t show earthquake occurrence in the traces. This happened to the stations on the outside periphery, and the most times with station OUZ. We think that this may be related to the location of the stations. Station OUZ is the furthest station away from the sources, and thus the time window in which the traces are shown could be out of range here. All data that exhibit these situations are not considered for the inversion. Other traces that are regarded too noisy are also not included.

local events

Similarly with the teleseismic events, there are 121 events taken into account, 13 of which are discarded due to bad data. The same situation with the local data is also observed here. Station HIZ and OUZ at the two ends of the array are the most often disregarded from the inversion because of this. Again on top of all of this, if a trace data is too noisy it would also be discarded.

teleseismic events


To assess the inversion performance, we plot the evolution of standard deviation (in seconds) and chi^2 throughout the iterative steps. We see that after 3 iterations, we won’t get any more significant improvement to the model.The graphs show that the overall improvement from the initial model is less than 10%. Relaxing the damping parameter does reduce the std deviation and chi^2 but only slightly. However doing so might introduce redundant features.

standard deviation and chi^2

inversion result

The following result is obtained after performing 3 iteration steps. We see that the Auckland Volcanic Field is at the edge of a predominantly negative velocity anomaly (in km/s). We are still unsure on what might have caused this observation. But the negative anomaly in Auckland is a characteristic for areas of active volcanic activity. For example in the southern extent of the model, the negative anomalies reflect the locations of Mt. Taranaki, and the Taupo Volcanic Zone.

The red boxes signify the extent of our model that we’re most confident on. However it doesn’t mean that features outside of this box are all artefacts, they are simply less accurate due to the configuration of our sources and stations.

We add here the travel-time residual distribution for the initial and final model after 3 iterations. There is little change as explained by the standard-deviation analysis. The histograms are centred near zero, although we do have travel-times that are very late or very early.

large scale real data tomography: 100+ local events

In this travel-time tomographic inversion we plan to use more events to get a significant result in terms of AVF subsurface velocity field.


Eighteen stations are used in the inversion. These are divided into 11 stations located around the central area, and 7 outer surrounding stations. We removed the low resolution channels associated with the station.


events and picks

We have 116 earthquakes around the region that we plan on using.

Almost all events have good trace data, which means that we have a good quality source input that we can work with. Our next step is to determine the first arrival pick with impeccable accuracy. The automatic trigger does a good job with picking the arrivals to a couple of seconds or so, but to consistently get an accuracy to the second or even less, picking by hand is still the go to method.

While we are doing the picking, some number of traces don’t show the event. We presume that this is mostly due to the location of the related station. The relative distances can make the time window in which the traces are shown becoming too short. For example, the northernmost station OUZ is the one that mostly show this condition. For the following results, the corresponding picks are excluded from the inversion, but we can revisit these picks and include them later on.

Even though the number of traces that don’t show an earthquake is small, it is still small things that we can improve on.


The split in the positive and negative anomaly (in km/s) which we also saw in the inversion with less events can be seen more distinctly this time around. Positive anomaly is prominent to the North and similarly, negative in the South. Another highlight is the negative anomaly right below Auckland (36.6-37.2 S, 174.6-175.2 E) at 50+ km.

We can even see a dip in the anomaly although we cannot be sure how much influence does the smearing effect has here. It is consistent with the primary paths direction from the South-East. However unlike any smearing that we have seen before in previous results, that are usually constricted to the outer and deep part of the model, here it is rather central. For example, the deep lineaments in the NS slice are likely to be smearing, and are in agreement with the three deepest events. Furthermore, our previous analysis from the checker-board test shows that at least for the first 70 km or so, the smearing is not too dominant. Therefore, we may be seeing a jump or a continuation of structure to depth commencing somewhere in the south Auckland.

Here are the plots showing how the root-mean-square travel-time residuals, standard deviation, and chi^2 change through each iteration. The first iteration is when the misfit drops drastically, afterward the plots quickly becoming stable. These plots justify the number of iterations we did as after three or four iterations we already have a convergence.

We have here the histogram of the travel-time residuals for initial and final model. There is no change between the two, but that is consistent with the standard-deviation plot. Most residuals are within a couple of seconds, but we do have some that can be more than 20 s late.

Resolution test on the AVF

We propose to perform P-wave tomography to infer the structure beneath the AVF, with Nick Rawlinson’s FMTOMO. Our first step is to establish the resolution that we can expect from this exercise. Below, we plot the locations of the earthquakes — a total of 1330 cataloged events, 681 local + 649 teleseismic — used in this test:



The local events made up the sources coming in from the southeast of the receivers array. Whereas teleseismic sources complement this with incoming rays from other directions but most notably from the Northwest.

We compute synthetic seismic P-wave arrival times for the AVF seismic stations, as if the earth structure under the AVF has a checker-board pattern. For this test, this will be our “true” model of the Earth. And so the objective of this test is to invert these synthetic travel times to recover the “true” model checker-board pattern. Regions where the inverted results show a checker-board pattern suggest a good resolution. And the size of each checker-board pattern imply the size of the smallest feature resolvable.

Coarse Checker-board Test

The first checker-board model shall have a wide grid spacing creating a fairly large checker-board patterns. In the following image, the ensuing checker-board model is represented by the three panels on the left while its recovery after the inversion is represented by the panels on the right.


Large checker-board test



  • Dimension of the region where resolution is most exemplary is signified by the red box, ranging from Northland to Waikato, and to 280 km depth.
  • Outside of the box, the smearing artefact is too prominent. Caused by ray paths having similar orientation.
  • Within the box, the size of the features that are resolvable, corresponds to the size of each checker-board, is in the order of 60 km.

Fine Checker-board Test

We reduce the grid spacing, creating a smaller sized checker-board model, inferring that smaller features now have the potential to be resolved.



  • The red box where resolution is commendable, is now span over the Auckland region, down to about 80 km deep.
  • Within this box, subsurface features on the scale of 30 km can be resolved.