The September 2000 Analysis Committee Report describes the SCIGN Analysis Committee's history and previous work.

At the September 2000 meeting, the Analysis Committee presented compared JPL and SIO results for six weeks before and after the Hector Mine earthquake. Although these results are already obsolete and will be replaced by the SCIGN combined solution, Nancy King continued to use this data set to develop scripts and compare positions and baselines. The comparison scripts will be useful for the combined solution.

In this experiment, there are 112 stations processed by both JPL and SIO. For each station, there are six time series: X_SIO, X_JPL, Y_SIO, Y_JPL, Z_SIO, and Z_JPL. Each time series is 12 weeks long. For each station, we obtained a SIO-JPL difference series by subtracting each day's JPL component from the SIO component. Previous work presented results for the mean and rms of each SIO-JPL difference series. These results highlighted the problem of reference frames. Each institution attempted to define a reference frame across the time of the Hector Mine earthquake. The first comparison of XYZ positions showed that the two realizations were not compatible. A second comparison, using the same stations with identical velocities and a priori positions, yielded much better results. Time series for 112 stations with SIO and JPL using different and "identical" reference frames show that the SIO-JPL difference changes at the time of the Hector Mine earthquake. This change is obvious when SIO and JPL use different reference frames. With the "same" reference frame, the change is much less marked and disappears altogether for many stations.

In this report, we present the SIO-JPL differences in a new way. For each day of the 12-week period, we computed the mean of the SIO-JPL differences for all 112 stations. Plots of mean SIO-JPL difference by day, show a clear change at the time of the earthquake. There is also an odd transient bias signal after the earthquake. However, the important feature is the change, of about 2 mm in X, 3 mm in Y, and 2 mm in Z, that occurs at the time of the earthquake. Clearly the reference frames are still not really "the same." This report's section on the SCIGN combined solution, describes how we decided to define a common reference frame.

It it useful to compare baselines because they are independent of reference frame. Since this data set contains 112 stations common to SIO and JPL, there are 112*111/2 = 6216 possible baselines. In practice there are fewer, because not all pairs of stations have enough common observation days and because several stations in this data set had antenna height errors (since corrected, but still a problem in this particular data set). For about 5800 baselines, we compared baseline length and transformed XYZ components to NEU. We used code provided by Yehuda Bock for this transformation. The scripts that automate this comparison will be used for the combined solution too.

For each baseline, there are eight time series: L_SIO, L_JPL, N_SIO, N_JPL, E_SIO, E_JPL, U_SIO, U_JPL. Each time series is 12 weeks long. For each baseline, we obtained a SIO-JPL difference series by subtracting each day's JPL component from the SIO component. We then computed the mean and rms of each difference time series.

The results are excellent. The mean baseline length difference is 0.2 mm, with rms scatter about the mean of 1.4 mm. The mean differences for the horizontal baseline components are on the order of 0.1 mm, with rms scatter about the mean of about 1.5 mm. For the vertical baseline component, the mean difference is about 0.5 mm and the rms scatter about the mean is about 5 mm.

SIO-JPL Baseline Differences

	Mean difference, mm	RMS, mm
Length	-0.20	1.4
N	-0.08	1.3
E	0.12	1.6
U	0.53	4.9

The table above shows average results for baselines of all lengths. Since geodetic baselines generally show proportional error, we plotted the baseline component differences as a function of baseline length. The table below contains links to these plots.

The table below shows the parameters of the straight line fits to the plots of SIO-JPL baseline component difference and RMS versus baseline length.

Best-fitting straight lines for SIO-JPL Baseline Difference Plots

	Mean difference	Mean difference	RMS	RMS
	Y-intercept, mm	Slope, mm/km	Y-intercept, mm	Slope, mm/km
L	0.02	0.0014	1.21	0.0012
N	-0.06	-0.0001	1.21	0.0009
E	0.08	0.0003	1.29	0.0025
U	0.52	0.0001	4.74	0.0010

The slopes and intercepts in the table are parameters of a simple model. The intercept is the mean difference or RMS for a zero-length baseline, and is, loosely speaking, the irreducible error in our system. In geodetic terminology, this is our "zero error." The slope is the proportional error. For short baselines, JPL and SIO results agree to better than 0.1 mm in the horizontal, and about 0.5 mm in the vertical. Corresponding RMS scatter is about 1.3 mm in the horizontal and 4.7 mm in the vertical. For baseline lengths of 700 km, the longest in this data set, this simple model predicts that the mean JPL-SIO difference is -0.15 mm, 0.31, and 0.58 mm in N, E, and U, respectively. The slope of the best-fitting straight line of RMS versus baseline length yields an estimate of the length-dependent standard deviation in the measurements: about 1 part per billion in length, north, and vertical, and 2.5 parts per billion in east.

To obtain the first official SCIGN product from SIO, JPL, and USGS solutions, the Hector Mine coseismic displacements, we simply computed the weighted mean and standard deviation of the coseismic displacement vector for each station. However, since each institution's results are derived from the same RINEX files and are therefore correlated, we were not sure how to obtain realistic uncertainties. The simple weighted mean will not work for positions, since it doesn't make sense to average positions (as opposed to position differences from one day to the next) calculated relative to different realizations of a local ITRF97 reference frame.

It turns out that even the agreement between SIO and JPL coseismic vectors was fortuitous, since these early results were obtained before JPL and SIO attempted to redefine the reference frame across the earthquake. As described above in the Continuation of September 2000 Comparison section, SIO and JPL eventually used the same reference stations with the same a priori positions and velocities. And yet, as the plots of mean SIO-JPL difference by day show, the mean SIO-JPL differences in position change at the time of the earthquake; the coseismic vectors will now not agree. This is probably because the details of SIO and JPL reference frame stabilization differ and the reference frames are not really the same.

This discouraging result shows us that, at some point, SIO h files and JPL stacov files must be adjusted together. MIT's GLOBK/GLORG software, designed to combine solutions in this way, is an obvious choice for forming a SCIGN combined solution. There may be other methods, such as QOCA, and we may investigate other software later.

To obtain the combined solution we have to do the following, bootstrapping as we iterate on these steps:

• (1) Resolve problems with antenna heights, etc.
  
• (2) Choose the stations that will define a preliminary SCIGN reference frame.
• (3) Rerun data at JPL and SIO, with common antenna heights and appropriate constraints.
  
• (4) Figure out how best to use of GLOBK/GLORG to combine SIO's h files and JPL's stacov files.
  
• (5a) Use GLOBK/GLORG to combine loosely constrained daily h and stacov files for the SCIGN reference stations, with tight constraints on the a priori IGS97 epoch positions and velocities. This step yields positions and velocities for a preliminary reference frame that we designate "SCIGN97 v.1."
• (5b) Understand and resolve any offsets in the reference stations. Repeat (5a) as necessary.
• (5c) Using a rigid-body Helmert transformation, transform all stations for that day into the SCIGN97-v.1 reference frame, based on the Helmert transformation defined using the preliminary reference stations.
• (5d) Combine all stations on all days. All "well-behaved" stations are then used to define the final SCIGN97 reference frame.
• (5e) Redo the rigid-body Helmert transformation on each day using the final SCIGN97 reference frame defined in (5d).
• (6) Set up an automated system to get h and stacov files in the same place and run GLOBK/GLORG.
  
• (7) Set up a Web page and ftp area where users can obtain the results in numerical and graphical form. Compare SIO, JPL, and the combined solution there.

Tasks (1) through (4) are done, although we have had to iterate some of these steps and may have to do so again. We thank Tom Herring for his help with task (4). Tasks (5a) and (5b) are in progress, and we are finding it frustrating and time-consuming to resolve the offsets. Task (6) is in progress also. Tasks 5(c-e) and (7) are not done yet.

As reference stations (task 2), we wanted stations with
a long and uncomplicated time history,
high-quality data,
good geographic distribution.


Of course, there is no set of stations matching all these criteria. Station histories that are long also tend to be complicated, with many equipment changes and offsets. Stations from Mexico and the BARGEN network, with desirable positions relative to the rest of SCIGN, tend to have short time histories. PVEP/PVE3 and CICE/CIC1 are pairs of different stations in roughly the same location. Balancing all these considerations, we selected the stations in the table below as the best set. It is very likely that we will revise this list later.

SCIGN97 v.1 reference frame map


Ref Station	Start Date	Antenna Ht, m	Events	JPL Offsets	SIO Offsets
ALAM (for now)	05/08/1999	0.1283 BPA	list	none	none
BLYT	01/13/1994	list	list	list	none
CAT1	06/25/1995	0.1614 BPA	07/27/1999 Receiver swap	list	none
CICE	03/24/1995	0.0814 BPA	none	none	none
CIC1	02/09/1999	0.07927	08/18/1999 Receiver swap	none	none
COSO	08/16/1995	list	list	list	none
DYER	04/28/1999	0.1386 BPA	12/16/1999 Firmware upgrade	none	none
ECHO	05/07/1999	0.1288 BPA	list	none	none
FERN	05/02/1999	0.1245 BPA	12/23/1999 Firmware upgrade	none	none
FRED	03/25/1999	0.1295 BPA	list	none	none
GUAX	01/19/2001	0.0083	none	not in list	none
PVEP	05/17/1993	list	list	list	01/03/1998
PVE3	09/27/2000	0.0083 BPA	none	not in list	none
RAIL (replace ALAM in future?)	04/30/1999	0.1298 BPA	list	not in list	none
SIO3	05/06/1993	list	list	list	04/12/2000
SPMX	10/14/1998	0.0083 BPA	list	none	none
TONO (replace ALAM in future?)	01/22/1999	0.1376 BPA	list	not in list	none
VNDP	05/25/1992	list	list	list	none

Here are results for the reference stations as of September 21, 2001. Since we are still working on this, results here may change.

JPL/SIO	JPL/Combined	SIO/Combined
ALAM	ALAM	ALAM
BLYT	BLYT	BLYT
CAT1	CAT1	CAT1
CIC1	CIC1	CIC1
CICE	CICE	CICE
COSO	COSO	COSO
DYER	DYER	DYER
ECHO	ECHO	ECHO
FERN	FERN	FERN
FRED	FRED	FRED
PVE3	PVE3	PVE3
PVEP	PVEP	PVEP
SIO3	SIO3	SIO3
SNI1	SNI1	SNI1
SPMX	SPMX	SPMX
VNDP	VNDP	VNDP

To estimate offsets, we used code written by John Langbein. This program, gsp_rate_off_3.f, uses our best estimate of the noise model to estimate offsets correctly. We fit a linear rate of change, annual and semi-annual terms, white noise, and flicker noise. The noise parameters are:

	White, mm	Flicker, mm
N	1.5	0.6
E	2.5	1.0
U	5.0	2.0

These parameters come from John's study of the PANGA data. The rule of thumb is that the flicker noise is 40% of the white noise.

Below are offsets found in the results for the SCIGN reference stations.

SIO Offsets

Site	Date	Component	Estimated offset, m	Estimated offset, m	Cause
			Standard least squares	John's program gps_rate_off_3
alam	1999 12 19	N?	not estimated	-0.0032 +/- 0.0022 (Note 1)	Hector Mine?
alam	continued				Addition of dome on 1999 12 21?
pvep	1997 01 01	U	0.1063	0.1087 +/- 0.0021	???
pvep	1997 04 23	U	-0.1117	-0.1089 +/- 0.0021	???
pvep	1998 01 13	E	0.0113	0.0105 +/- 0.0011	Antenna swap 1998 01 03?
sio3	1998 01 13	N	-0.0067	-0.0002 +/- 0.0009 (Note 2)	Receiver swap
sio3	1998 01 13	E	0.0057	0.0010 +/- 0.0011	Receiver swap
sio3	2000 04 12	U	-0.0684	-0.0553 +/- 0.0019	Antenna swap
sio3	2001 03 01	U	0.0301	0.0281 +/- 0.0024	???
sni1	2000 12 21	U	not estimated yet	-0.0140 +/- 0.0021	???

Note 1: Because of the gap in SIO data, this offset is the same as the ALAM N 1999 10 25 offset in JPL data.
Note 2: After deleting glitches, offset is at 1998 01 21.


JPL Offsets

Site	Date	Component	Estimated offset, m	Cause
			John's program gps_rate_off_3
alam	1999 10 25	N	-0.0024 +/- 0.0010 (Note 1)	Hector Mine
blyt	Starts Aug 1999	N "slow offset"	Not estimated	???
blyt	1999 10 24	E	0.0066 +/- 0.0012	Hector Mine
blyt	1996 05 05	U	0.0084 +/- 0.0021	???
cice	1996 05 05	N	0.0026 +/- 0.0009	???
cice	1996 05 05	U	-0.0070 +/- 0.0023	???
coso	1996 05 05	N	-0.0028 +/- 0.0009	???
coso	1996 05 05	U	-0.0038 +/- 0.0021	???
pvep	1998 01 13	N	-0.0070 +/- 0.0009	Antenna swap 1998 01 03?
pvep	1998 01 13	E	0.0168 +/- 0.0013	Antenna swap 1998 01 03?
pvep	1996 05 10	U	0.0208 +/- 0.0020	???
pvep	1998 01 13	U	-0.0440 +/- 0.0024	Antenna swap 1998 01 03?
sio3	1999 10 16	N	0.0046 +/- 0.0008	Hector Mine
sio3	1998 03 05	E	0.0032 +/- 0.0011	???
sio3	2000 04 12	U	-0.0974 +/- 0.0019	Antenna swap
sni1	2000 12 20	U	-0.0178 +/- 0.0021	???
vndp	1996 05 07	U	0.0090 +/- 0.0022	Receiver swap

Note 1: Offset date doesn't match SIO date because SIO data has a gap.Offset date doesn't match Hector Mine date because JPL data has a short gap.

1) In analyzing noise for a subset of the Panga array, I found that there is a relation between the amount of flicker noise relative to white noise. This is important since once we know the ratio between these components, we can confidently estimate the size of the offset (and other parameters) but be incorrect with regards to the absolute size of the standard errors. The plot pl_fl_wh illustrates the relation between flicker and white noise. I used 7 sites with differing styles of monuments (2 are IGS? style monuments (drao and albh), 2 are braced monuments, one is a pin in bedrock, and the other 2 are CORS (coast guard?) sites; as expected the CORS sites have about twice the noise as the other sites; Also, the time series lengths are between 2+ years to 7 years).

2) Given that we know the relative weights between white and flicker noise (flicker is about 40% of the white), I then explored the weighting used to estimate each parameter, ie;

x_i= W_ij * d_j

where x_i are the model parameters, W_ij is the weighting matrix (The inverse of the design matrix), and d_j are the data (GPS time series)

The plot m1 shows the data weighting for our case for 4 parameters; the nominal value (or average), and 3 offsets, one around 90 days, the second at 180, and the 3rd at 270. For the offset around 90 days, the weights are approximately exponential with about a 20 day time constant. However, because the data have temporial covariance, the weights "see" the effect of the offsets at other times. For the case where the data are uncorrelated, see the results in m1a; here, the results are more intuitive; for the 1st offset, average the 90 days before and after the offset.

The next complication is adding the rate to the list of unknowns. This is shown in m2. Now, for the 2nd offset, you'll see a slight slope to the weights.

3) So, I conclude that it is somewhat difficult to prescribe a weight function to estimate the size of the offset; rather, I'd suggest doing 'nearly' the full problem; assume that you know the amounts of white and flicker noise from part 1, compute the covariance function, take its inverse (the most cpu intensive part!), compute the design matrix; and with the inverse covariance matrix, estimate the size of the offsets and rates (or even rate changes). I can modify my existing code to make it easier to use. The question is the amount of data and cpu time.

For a prototype code, I got the following run-times using the code compiled with g77 (and misc options) under linux on a 500Mhz cpu 1 year 1.2 seconds
2 years 12 seconds
3 years 44 seconds
4 years 106 seconds


or cpu_time= 1.5 * (time_span_years)^3.08
This corresponds to about half an hour for 10 years of data, which is certainly not bad. It is about a factor of 3-4 times slower on my 200MHz ultra-sparc 2.

So, with lots of data, it bogs down quickly (One might chop the time-series in half with minimal lose of precision in offsets) Also, when estimating the size of offsets, one can decimate the time-series in sections far removed from the time of offset and achieve faster calculation of the offset (and rates and what-not).

4) The other issue is the size of offset that is detectable. Here are a few estimates (1 sigma).

For north where white noise is between 1 to 1.5 mm (and flicker is 40% of white) the standard error in estimating offsets is 0.44 to 0.65 mm

For east where white noise is between 2 and 2.5 mm (...), the standard error in offsets is 0.9 and 1.1 mm

and, for Up, with white noise is between 3 and 5 mm, the standard error in offsets is 1.3 and 2.2 mm.

Of course, the magnitude of these errors assume that each time-series has some sort of regional "filtering". Thus, one might do a visual scan of the time-series to look for offsets of these sizes that aren't due to EQs or some other human interference (antenna swaps).

5) Above, I made a comment about doing 'nearly the full problem' of estimating offsets. The most complete way is the Max. Likelihood estimation that both Hadly and I have done. It is very cpu intensive. (It takes about 1-day to grind through a 7 year long time-series!) But, once you think you know when offsets occurred and when the rate change might have occurred, you'll have both of these parameters plus having the amplitudes of the noise components. I'd recommend do this only on a few stations; perhaps the reference stations.....

The alternative is to identify the offsets and remove them with the proceedure in #3; then fill-in missing data with white noise (assuming there aren't too many gaps or too long of a gap....) and compute the power spectral density using no smoothing. After ignoring the power at the lowest frequencies, fit a curve representing flicker plus white noise; my experience with the Panga data indicates that this method does gives reliable results. If the ratio between the estimate of white and flicker noise deviates grossly from the linear trend in pl_fl_wh.ps, then one might re-estimate the size of the offsets using #3 again, but with a better noise model.

There are a few ways of normalizing the amplitude of flicker noise. Hadley did it one way, and I do it another way. The difference is a factor of 2 (Hadley would get numbers a factor of 2 bigger than I for flicker noise).

Each offset has a 5-character event label. The beginning of each time series is an event called START.

Master list of offsets. This list contains many non-SCIGN stations.
SCIGN offsets. This is a subset of the master list.

SIO offset list, put together by Rosanne Nikolaidis.

Photos of BGIS vandalism, by Aris Aspiotes
Monument and dome
Dome
Closeup of brick marks on dome
Closeup of nicks on antenna

Our first reaction was that the antenna must have been damaged. People at SOPAC (Yehuda Bock, Matt van Domselaar, and Rosanne Nikolaidis) used the baseline from DYHS to BGIS to establish the date, time, and size of the offset. Further analysis at SOPAC showed that the offset disappeared once a new radome was installed. Click here to see the single-epoch time series for the day on which the vandalism occurred.

We conclude that

• We need domes, so that they get smashed instead of antennas.
• Vandals have to work hard to break a dome, but they can do it.
• Damage to a dome can offset the time series. We need to figure out what's happening.
• Replacement with an intact dome reverses the offset. The Bell Gardens incident suggests that SCIGN domes are very uniform, and hence almost interchangable.