5.5. Generating Time Series and Velocity Solutions
To generate time series or estimate velocities from several years of data, you need to run glred or globk directly rather than using sh_glred.
As with creating irregular combined h-files, you can create the .gdl file using the ls command from the directory you intend to use for the multi-year analysis.
For example, for a multiyear analysis in pnw/vsoln/ using combined h-files from individual years,
$ ls ../pnw??/gsoln/*GLX > pnw_all.gdl
or, if you have not aggregated the data,
$ cat ../pnw??/gsoln/*.gdl > pnw_all.gdl
using cat with the .gdl files produced by sh_glred, rather than ls with the .glx and .GLX files in glbf/ avoids having to add the + signs for multiple files within a day when running glred for the long-term time series.
(When running globk to get velocities all files are combined regardless of whether the + signs are present or not.)
Your first step should always be to generate time series so that you can evaluate the statistics of your data and determine if there are any outliers that will distort the velocity solution. To run glred directly to produce a time series from this file, set (uncomment e.g.
$ glred 6 glred.prt glred.log pnw_all.gdl globk.cmd >& glred_rep.out
The solutions for all h-files, obtained independently, are written into a single print file (specified by prt_opt for finite constraints using globk, org_opt for generalized constraints using glorg) so you can plot the results with
$ sh_plot_pos -f glred*.org -r -t NONE
(Alternatively you can create/save residual .pos files for inspection and repeated plotting by running tssum first and then sh_plot_pos with *.pos instead of *.org for the input files.)
The first check of the solution for each epoch is the quality of the stabilization.
You can get a summary of the statistics by grep’ing on POS STAT in the print file.
The NRMS values for the stabilization should be near unity in all three components, and the WRMS 1–2 mm in horizontal and 3–10 mm in vertical (potentially higher if vertical was down-weighted more than the default factor of 10).
Check to make sure that no epoch has too few sites in the stabilization.
If you find any problems, look at the solution itself for that epoch, noting which sites were removed from the stabilization and/or have large adjustments in the position summary.
Next look at the histograms of nrms and wrms generated by sh_plot_pos to see if the distribution of scatters among your sites is approximately Gaussian with the median nrms ~ 1.0.
Finally, unless you have a network of hundreds of sites (in which case you are probably too experienced to need this chapter), you should view the time series for every site to be included in your solution, whether or not it is of direct scientific interest.
Look for outliers at individual epochs and sites with unusually high nrms (> 1.5–2.0, depending on how you have weighted the data).
If a large number of sites have outliers at a single epoch, review the GAMIT processing and the stabilization.
Once you have identified any problematic position estimates and created rename or sig_neu commands to remove them, you should repeat the time series to verify that your entries accomplished what you intended.
Once you are satisfied with your time series, you can run globk with the same .gdl file to obtain a velocity estimate:
$ globk 6 globk.prt globk.log. pnw_all.gdl globk..cmd VEL >&! glred_vel.out
where the VEL option invoked with the template globk.cmd file will turn on velocity estimates and printout:
apr_neu all 10 10 10 1 1 1
org_opt PSUM VSUM GDLF CMDS
and glorg.cmd, invoked from within globk.cmd includes stabilization for velocity as well as position:
* Set parameters to estimate in stabilization
pos_org xtran ytran ztran xrot yrot zrot
rate_org xtran ytran ztran xrot yrot zrot
The .org file output will contain stabilization information and a summary table for velocities as well as positions.
Also, the experiment list will now contain all of the files used and the computed chi-square increments as each successive h-file is added to the solution:
EXPERIMENT LIST from vel.srt
# Name SCALE Diag PPM Forw Chi2 Back Chi2 Status
1 H950820_PNW.GLX 1.000 0.000 0.115 -1.000 USED
2 H970607_PNW.GLX 1.000 0.000 0.220 -1.000 USED
3 H980813_PNW.GLX 1.000 0.000 0.234 -1.000 USED
5 H990516_PNW.GLX 1.000 0.000 0.184 -1.000 USED
6 H010719_PNW.GLX 1.000 0.000 0.295 -1.000 USED
With a weighting of the coordinates estimates such that the uncertainties from each (combined) h-file are ~ 1 mm horizontal and ~ 5 mm vertical, the chi-square increments are typically less than 0.5 for regional data alone and 0.4 to 1.0 for global data. If an h-file is based on a poor GAMIT solution or there is a model incompatibility with previous h-files, the chi-square increment will be anomalous. When this occurs, you will need to return to the time series to see which sites show outliers that would create an incompatibility. A common example is a wrong antenna height for a station, which will show up as a height anomaly in the time series. How much chi-square is increased by a problem at one site will depend both on the size of the error and how many sites are in the solution. Very high chi-square values can result if you have two sites in the solution with the same name but at different locations. To guard against this problem, you should run glist before you start the processing (see Section 5.1 of the GLOBK Reference Manual).
To plot the velocities, use
$ sh_plotvel -f globk_vel.org -R244.5/251/40.5/46.5
where the -f argument is the .prt, .org, or .vel file from which the SUMMARY VELOCITY is to be plotted, and the second command-line entry gives the longitude and latitude range for the plot using standard GMT commands.
See the documentation at the top of the script (or type sh_plotvel -help) for additional options for displaying site names, scaling the uncertainties, setting the confidence level for error ellipses, and adding topography, fault maps, and other geographic features.
The scripts sh_map_calif, sh_map_china, and sh_map_elements in ~/gg/com/ may be helpful as templates.
For large networks, particularly those including continuous stations with discontinuities and non-linear motion due to earthquakes, the interactive Matlab-based programs tsview and velview can be powerful tools.
Their use and installation ar documented on the GAMIT/GLOBK webpage (http://geoweb.mit.edu/~tah/GGMatlab/).
As a further aid to handling large data sets, we have recently developed programs (tssum, tsfit) that allow “prototyping” of GLOBK solutions by manipulating time series to identify discontinuities and non-linear effects quickly (see Chapter 4 of the GLOBK Reference Manual).
It will often be the case that the set of stations you use for stabilization of the time series for editing purposes define a reference frame that is inherently weaker in your region of interest than the final velocity solution—either because the editing frame encompasses a larger region and hence is more affected by spatially correlated errors, or because the number of well-known stations available for stabilization of individual h-files is much less than the total number of stations in the solution.
The time series that best represents your final velocity solution is one that is stabilized with all of the stations in the solution.
Thus, after you obtain this solution, you should rerun glred with an enlarged stab_site list and an .apr file generated from the velocity solution.
To get the .apr file, use sh_exglk, and to generate a stab_site list to include all sites with position and velocity uncertainties below a specified level program org2stab.