New Mexico Tech
Earth and Environmental Science

ERTH 455 / GEOP 555 - Geodetic Methods

Lab 11: Modeling - Timeseries Analysis w/ CATS

"I like my crust deformed."
UNAVCO bumper sticker

Note that you DO have to work on redoubt for some things today!

Introduction

Last week, you wrote some code to estimate time series properties. This week, you'll use some other software to do that. Now that you know how to process data into position estimates and can convert them into, e.g., a local coordinate system, we make the assumption that you can repeat that for every day that you have data. If you concatenate the daily position solutions you get a time series and from that, we can finally start learning something about the long and short term geologic processes that affect the site.

The goal for this week is familiarize yourself a bit with a the program CATS (Create and Analyze Time Series) that analyzes time series by fitting a multi-parameter model to the data.

1. Plot the data

Change to your working directory and copy the test data:

    $> cd $GEOP555
    $> mkdir lab11
    $> cd lab11
    $> cp -r /usr/local/CATS/Bezy_2 ./
    

A directory listing should look similar to this:

    redoubt:lab11> ls Bezy_2/
    BEZD.cats  BEZR.cats  BZ01.cats  BZ03.cats  BZ05.cats  BZ07.cats  BZ09.cats  KAMD.cats  KLU.cats   MAYS.cats
    BEZH.cats  BZ00.cats  BZ02.cats  BZ04.cats  BZ06.cats  BZ08.cats  ES1.cats   KBG.cats   KLUC.cats  PETS.cats
    

For now, we will limit ourselves to working only on 1 file - BZ09.cats, which contains the timeseries in local NEU format for a station near Bezymianny Volcano, Kamchatka (beautiful spot!). The contents of the files look like this:

    2006.7644 0.00000 0.00000 0.00000 0.00110 0.00140 0.00320 
    2006.7671 -0.00579 0.00599 -0.00170 0.00100 0.00130 0.00290 
    2006.7699 -0.00412 0.00231 -0.00530 0.00100 0.00130 0.00290 
    2006.7726 -0.00468 0.00318 0.00380 0.00100 0.00130 0.00290 
    2006.7753 -0.00936 0.00212 -0.00100 0.00100 0.00130 0.00290 
    2006.7781 -0.00590 0.00050 -0.00310 0.00100 0.00130 0.00290 
    2006.7808 -0.00590 0.00412 0.00060 0.00100 0.00130 0.00290 
    2006.7836 -0.00479 0.00456 -0.00260 0.00100 0.00130 0.00290 
    2006.7863 -0.00590 0.00300 0.00160 0.00100 0.00130 0.00290 
    2006.7890 -0.00991 0.00225 -0.00180 0.00100 0.00130 0.00290 
    2006.7918 -0.00991 -0.00162 0.00620 0.00100 0.00130 0.00280
    ...
	

Where columns are:

	decimal year | North | East | Up | North error | East Error | Up Error (all meters)
	

Plot the north, east and up components of station BZ09 over time in 3 panels in a single figure. Note that this should be a discrete plot, i.e. disconnected dots. With an error bar function, you add the respective errors.

2. Analyze the Data

The benefit of CATS is that it allows you to choose between a range of stochastic models for noise approximation. Traditionally, we've been using white noise approximations to estimate our uncertainties in the site velocities. By now we know that the noise characteristics are more complicated than this.

You can find a lot of the detail on the various stochastic models in the manual. In addition to estimating noise parameters, the program also allows you to estimate slopes, and offsets at given times, as well as annual and semiannual parameters. Let's do that!

In your lab directory run:

    $> cats Bezy_2/BZ09.cats --model wh: --verbose
	

This estimates, for this station, the intercept, slope and white noise parameters according to the models in the manual, for all three components:

		...
+NORT  MLE    :    -2859.806080
+NORT  INTER  :        2.4353 +-      0.2396
+NORT  SLOPE  :      -13.6079 +-      0.0818
+NORT  WHITE NOISE
+NORT  WH     :        2.9409 +-      0.0615		
		...
	

(units are mm) Let's see what happens to these parameters with a more realistic noise model, when we add flicker noise.

    $> cats Bezy_2/BZ09.cats --model pl:k-1 --model wh: --verbose --output BZ09.out
	

The output is being redirected into BZ09.out. You can get the results to the screen with:

	$> grep "^+" BZ09.out
    

Having the slopes and intercepts for both models, add them to your plot using the line equation:

		y = a+b*t
	

Where t is time (set up a vector over the times given in the file or draw a line for 2 points in time), a is intercept and b is slope. Note that the intercept is defined at t0, which is the start of the first year of data (here: 2006.0), so you will have to subtract that value from your times in the above equation.

Some other parameters that we are often interested in are those for annual and semi-annual sinusoids, e.g., seasonal loading:

    $> cats Bezy_2/BZ09.cats --model pl:k-1 --model wh: --verbose --sinusoid 1y1 --output BZ09_seasonal.out		
	

The sinusoid parameter allows to add the period, time flag and harmonics. Here, the period is 1 year (1y) and we're interested in the first harmonic (1y1), i.e. the semiannual solution. Note that CATS appends to output files! So either remove the old one or create the new output file as I named it here: BZ09_seasonal.out.

The result file will now contain information on COS/SIN components. You could add these to your plots.

3. Add an "earthquake" and reanalyze

What happens to your velocity estimate if you have an actual offset in the data? Let's add an artificial offset:

    $> awk '{if ($1 > 2008.0) {print $1, $2+1, $3-0.75, $4+0.5, $5, $6, $7} else {print $0;}}' Bezy_2/BZ09.cats > Bezy_2/BZ09_steps.cats
	

Plot the new file BZ09_steps.cats

Now re-estimate the velocity, e.g,:

    $> cats Bezy_2/BZ09_steps.cats --model pl:k-1 --model wh: --verbose --output BZ09_steps.out
	

How has the velocity changed? Is this a realistic estimate for the long-term velocity of the site?

For known steps in time, CATS has a feature to help estimate the offsets and still get a long-term velocity. Create a new file Bezy_2/BZ09_step_fixed.cats and add information on the offset:

    $> cp Bezy_2/BZ09_steps.cats Bezy_2/BZ09_steps_fixed.cats
	
To the very top of Bezy_2/BZ09_steps_fixed.cats add the line:
    # offset 2008.0 7
	

In the manual, find what the 7 means

Now rerun the solution and compare the output. You should find information on the offset in all three components. How do the velocities compare to the "no-earthquake" data and the estimates without offset correction? What is the estimated offset for the components? Why isn't it exactly what we put in?

    $> cats Bezy_2/BZ09_steps_fixed.cats --model pl:k-1 --model wh: --verbose --sinusoid 1y1 --output BZ09_fixed.out
	

Deliverables: (submit via canvas!)

rg <at> nmt <dot> edu | Last modified: December 18 2017 17:46.