Introduction¶
This exercise is intended to introduce the concept of GNSS-Acoustic travel time residuals, how they can introduce bias into a GNSS-A positioning inversion, and how we can use them to refine an array position.
A more formal definition of GNSS-A residuals may be found in the GNATSS documentation, but the short version is that these residuals are the difference between two two-way travel time (TWTT) estmates between the sea surface transducer and the seafloor transponders, that measured in the field and that modelled after the fact when preparing the positioning inversion.
The model used to estimate TWTTs is generally some flavor of an acoustic raytracing operation between the sea surface transducer and the seafloor transponder, and relies on information about the sound velocity structure of the water column as well as the position of the transducer and seafloor transponders. In principle, if we had perfect knowledge of these components we should be able to replicate the measured TWTT, but this is never achieved in practice because the ocean is dynamic, the position of the transducer cannot be perfectly estimated with GNSS, and the transponders move with the seafloor away from their initial positions (albeit slowly). From a different perspective, if we can minimize the signal in the residuals contributed by the oceanographic and transducer positioning uncertainties then the remaining varitations should be the result of seafloor transponder displacements.
The GNSS-A array positioning performed by GNATSS is based upon this observation, but in order to work it depends on the user minimizing the contributions of undesired signals to the residuals through a process called residual flagging.
Files Present¶
In the directory ../data/exercise-2 you will find the following data:
- A configuration file
NCB1_YYYY_config.yamland output folderoutput_YYYYfor each of the years 2018, 2021, 2022, 2023, and 2024 - A folder
ctdthat contains sound velocity information for these data
In this exercise we will consider data collected at a GNSS-Acoustic array NCB1 offshore Coos Bay, OR, from 2018-2024 for the Near-Trench Community Geodesy Experiment. The inversions will be managed by the configuration files and saved in output directories, one for each year the array was surveyed during. The configuration files have been pre-made with all of the required information and should not need to be modified.
To expedite the processing, we have already generated gps_solution.csv files for each survey using the posfilter module of GNATSS, which are present in the output folders. Because of this, users will only need to run the solver module when processing these data. In fact, posfilter parameters are missing from the configuration files, as are the required data files for the posfilter, so users will be required to skip the posfilter during this exercise. The simplest way of doing this is to invoke the --solver flag if running GNATSS from the command line, in python the user may set skip_posfilter=True when executing run_gnatss().
Starting a GNATSS Inversion¶
We will start by considering the data collected during the most recent survey in 2024. If you are running GNATSS on the command line, process the data using the command gnatss run --solver ../data/exercise-2/NCB1_2024_config.yaml. Otherwise, we may execute the equivalent in python with the following:
from gnatss.main import run_gnatss# Set input variables
config_yaml = "../data/exercise-2/NCB1_2024_config.yaml"
skip_posfilter = Trueoutput = run_gnatss(
config_yaml,
skip_posfilter=skip_posfilter
)Starting GNATSS ...
Gathering sound_speed at ../data/exercise-2/ctd/CTD_NCB1_Ch_Mi
Gathering gps_solution at ../data/exercise-2/output_2024/gps_solution.csv
Loading sound_speed from ../data/exercise-2/ctd/CTD_NCB1_Ch_Mi
Loading gps_solution from ../data/exercise-2/output_2024/gps_solution.csv
Loading deletions from ../data/exercise-2/output_2024/
Filtering out data outside of distance limit...
Pre-filtering data with fewer than 3 replies...
Computing harmonic mean...
lat=43.58017583 lon=-125.097547806 height=-1076.6562 alt=-1076.6562 internal_delay=0.08 sv_mean=1481.682 pxp_id='NCB1-3' azimuth=57.44 elevation=39.67
lat=43.566349513 lon=-125.097616709 height=-1116.0454 alt=-1116.0454 internal_delay=0.36 sv_mean=1481.695 pxp_id='NCB1-2' azimuth=-58.75 elevation=39.38
lat=43.57274798 lon=-125.115119432 height=-1140.9439 alt=-1140.9439 internal_delay=0.64 sv_mean=1481.704 pxp_id='NCB1-1' azimuth=-175.59 elevation=39.52
Finished computing harmonic mean
Preparing data inputs...
Perform solve...
--- 30521 epochs, 91563 measurements ---
After iteration: 1, rms residual = 152.28 cm, error factor = 80.191
NCB1-3
D_x = -8.135402e-02 m, Sigma(x) = 5.604466e-05 m
D_y = -6.186801e-01 m, Sigma(y) = 5.608682e-05 m
D_z = 4.513089e-01 m, Sigma(z) = 5.206306e-05 m
NCB1-2
D_x = -8.135402e-02 m, Sigma(x) = 5.604466e-05 m
D_y = -6.186801e-01 m, Sigma(y) = 5.608682e-05 m
D_z = 4.513089e-01 m, Sigma(z) = 5.206306e-05 m
NCB1-1
D_x = -8.135402e-02 m, Sigma(x) = 5.604466e-05 m
D_y = -6.186801e-01 m, Sigma(y) = 5.608682e-05 m
D_z = 4.513089e-01 m, Sigma(z) = 5.206306e-05 m
After iteration: 2, rms residual = 104.39 cm, error factor = 47.617
NCB1-3
D_x = -1.384231e-04 m, Sigma(x) = 5.602523e-05 m
D_y = 3.227520e-07 m, Sigma(y) = 5.607893e-05 m
D_z = 9.380519e-05 m, Sigma(z) = 5.205427e-05 m
NCB1-2
D_x = -1.384231e-04 m, Sigma(x) = 5.602523e-05 m
D_y = 3.227520e-07 m, Sigma(y) = 5.607893e-05 m
D_z = 9.380519e-05 m, Sigma(z) = 5.205427e-05 m
NCB1-1
D_x = -1.384231e-04 m, Sigma(x) = 5.602523e-05 m
D_y = 3.227520e-07 m, Sigma(y) = 5.607893e-05 m
D_z = 9.380519e-05 m, Sigma(z) = 5.205427e-05 m
After iteration: 3, rms residual = 104.39 cm, error factor = 47.617
NCB1-3
D_x = -1.509118e-10 m, Sigma(x) = 5.602523e-05 m
D_y = 4.309075e-10 m, Sigma(y) = 5.607893e-05 m
D_z = -7.087753e-10 m, Sigma(z) = 5.205426e-05 m
NCB1-2
D_x = -1.509122e-10 m, Sigma(x) = 5.602523e-05 m
D_y = 4.309086e-10 m, Sigma(y) = 5.607893e-05 m
D_z = -7.087751e-10 m, Sigma(z) = 5.205426e-05 m
NCB1-1
D_x = -1.509119e-10 m, Sigma(x) = 5.602523e-05 m
D_y = 4.30908e-10 m, Sigma(y) = 5.607893e-05 m
D_z = -7.087755e-10 m, Sigma(z) = 5.205426e-05 m
---- FINAL SOLUTION ----
NCB1-3
x = -2660361.453 +/- 5.602523e-05 m del_e = 0.289 +/- 5.971492e-05 m
y = -3785656.3584 +/- 5.607893e-05 m del_n = -0.0543 +/- 6.148844e-05 m
z = 4373657.0134 +/- 5.205426e-05 m del_u = 0.7118 +/- 4.057857e-05 m
Lat. = 43.580175341632085 deg, Long. = -125.09754422671408, Hgt.msl = -1075.9443863983117 m
NCB1-2
x = -2660958.2954 +/- 5.602523e-05 m del_e = 0.289 +/- 5.971492e-05 m
y = -3786495.9768 +/- 5.607893e-05 m del_n = -0.0541 +/- 6.148821e-05 m
z = 4372517.1234 +/- 5.205426e-05 m del_u = 0.7118 +/- 4.057891e-05 m
Lat. = 43.56634902617636 deg, Long. = -125.09761313050257, Hgt.msl = -1115.3335735795908 m
NCB1-1
x = -2661822.702 +/- 5.602523e-05 m del_e = 0.2892 +/- 5.971572e-05 m
y = -3785267.4497 +/- 5.607893e-05 m del_n = -0.0541 +/- 6.148772e-05 m
z = 4373014.9389 +/- 5.205426e-05 m del_u = 0.7118 +/- 4.057848e-05 m
Lat. = 43.57274749300744 deg, Long. = -125.11511585101815, Hgt.msl = -1140.2321436132925 m
------------------------
There are 0 outliers found during this run.
Exporting residuals to residuals.csv ...
Successfully exported data to ../data/exercise-2/output_2024/residuals.csv
Exporting outliers to outliers.csv ...
No data found for outliers, skipping export!
Exporting distance_from_center to dist_center.csv ...
Successfully exported data to ../data/exercise-2/output_2024/dist_center.csv
Exporting process_dataset to process_dataset.nc ...
Successfully exported data to ../data/exercise-2/output_2024/process_dataset.nc
Successfully exported data to ../data/exercise-2/output_2024/residuals.png
Successfully exported data to ../data/exercise-2/output_2024/residuals_enu_components.png
Finished GNATSS.
Let’s talk briefly about what we are looking at in these plots.
The upper figure, saved as residuals.png, is a plot of the TWTT residuals. These residuals have been converted from units of seconds to units of centimeters by mulitplying them by the harmonic mean sound velocity. There is one time series of residuals for each transponder (so 3 in total for this array) and these time series are all plotted on the same axis. Note that “0 cm” is relative here and defined by the average of the three time series.
The lower figure, saved as residuals_enu_components.png is a rough estimate of the offset of the array center over the span of the survey. It is calculated by assuming the TWTT residuals are vectors oriented along the ray path from the sea surface above the array center to the indivdual transponders, decomposing these vectors into E, N, and U components, and averaging them.
We can see in these plots that most of the signal in the residuals is focussed around 0 cm, so we are doing a decent job of predicting the TWTTs. However, there are a number of “outliers” focussed along vertical lines that correspond to large offsets in the ENU plot. These are data artifacts, usually the result of GNSS cycle slips but in rare and extreme cases can be improperly recorded travel times (such as if a digit is dropped from the recording). Because these points are clearly non-physical measurements, we should remove them from the inversion in order to prevent them from biasing our result.
Flagging Residuals with --residual-limit¶
GNATSS has the funtionality to identify and remove these outlier residuals using a few different strategies. Two means that a user can set are a distance threshold from the center of the array with --distance-limit and by setting a gps_sigma_limit in the configuration file. However, the primary means for flagging residuals is to set an outlier threshold, which directs GNATSS to flag any TWTT residual for removal if it has a magnitude greater than the threshold. Using the example above, we want to remove the outliers that are reflected in the vertical residual spikes but not the residuals close to the main time series near 0. We should be able to do this by setting an outlier threshold of a few hundred cm, which can be done by re-running the previous GNATSS inversion using the --residual-limit XXX flag with XXX being the desired threshold in cm.
In python:
distance_limit = 150
residual_limit = 250
output = run_gnatss(
config_yaml,
distance_limit=distance_limit,
residual_limit=residual_limit,
skip_posfilter=skip_posfilter
)Starting GNATSS ...
Gathering sound_speed at ../data/exercise-2/ctd/CTD_NCB1_Ch_Mi
Gathering gps_solution at ../data/exercise-2/output_2024/gps_solution.csv
Loading sound_speed from ../data/exercise-2/ctd/CTD_NCB1_Ch_Mi
Loading gps_solution from ../data/exercise-2/output_2024/gps_solution.csv
Loading deletions from ../data/exercise-2/output_2024/
Filtering out data outside of distance limit...
Pre-filtering data with fewer than 3 replies...
Computing harmonic mean...
lat=43.58017583 lon=-125.097547806 height=-1076.6562 alt=-1076.6562 internal_delay=0.08 sv_mean=1481.682 pxp_id='NCB1-3' azimuth=57.44 elevation=39.67
lat=43.566349513 lon=-125.097616709 height=-1116.0454 alt=-1116.0454 internal_delay=0.36 sv_mean=1481.695 pxp_id='NCB1-2' azimuth=-58.75 elevation=39.38
lat=43.57274798 lon=-125.115119432 height=-1140.9439 alt=-1140.9439 internal_delay=0.64 sv_mean=1481.704 pxp_id='NCB1-1' azimuth=-175.59 elevation=39.52
Finished computing harmonic mean
Preparing data inputs...
Perform solve...
--- 30521 epochs, 91563 measurements ---
After iteration: 1, rms residual = 152.28 cm, error factor = 80.191
NCB1-3
D_x = -8.135402e-02 m, Sigma(x) = 5.604466e-05 m
D_y = -6.186801e-01 m, Sigma(y) = 5.608682e-05 m
D_z = 4.513089e-01 m, Sigma(z) = 5.206306e-05 m
NCB1-2
D_x = -8.135402e-02 m, Sigma(x) = 5.604466e-05 m
D_y = -6.186801e-01 m, Sigma(y) = 5.608682e-05 m
D_z = 4.513089e-01 m, Sigma(z) = 5.206306e-05 m
NCB1-1
D_x = -8.135402e-02 m, Sigma(x) = 5.604466e-05 m
D_y = -6.186801e-01 m, Sigma(y) = 5.608682e-05 m
D_z = 4.513089e-01 m, Sigma(z) = 5.206306e-05 m
After iteration: 2, rms residual = 104.39 cm, error factor = 47.617
NCB1-3
D_x = -1.384231e-04 m, Sigma(x) = 5.602523e-05 m
D_y = 3.227520e-07 m, Sigma(y) = 5.607893e-05 m
D_z = 9.380519e-05 m, Sigma(z) = 5.205427e-05 m
NCB1-2
D_x = -1.384231e-04 m, Sigma(x) = 5.602523e-05 m
D_y = 3.227520e-07 m, Sigma(y) = 5.607893e-05 m
D_z = 9.380519e-05 m, Sigma(z) = 5.205427e-05 m
NCB1-1
D_x = -1.384231e-04 m, Sigma(x) = 5.602523e-05 m
D_y = 3.227520e-07 m, Sigma(y) = 5.607893e-05 m
D_z = 9.380519e-05 m, Sigma(z) = 5.205427e-05 m
After iteration: 3, rms residual = 104.39 cm, error factor = 47.617
NCB1-3
D_x = -1.509118e-10 m, Sigma(x) = 5.602523e-05 m
D_y = 4.309075e-10 m, Sigma(y) = 5.607893e-05 m
D_z = -7.087753e-10 m, Sigma(z) = 5.205426e-05 m
NCB1-2
D_x = -1.509122e-10 m, Sigma(x) = 5.602523e-05 m
D_y = 4.309086e-10 m, Sigma(y) = 5.607893e-05 m
D_z = -7.087751e-10 m, Sigma(z) = 5.205426e-05 m
NCB1-1
D_x = -1.509119e-10 m, Sigma(x) = 5.602523e-05 m
D_y = 4.30908e-10 m, Sigma(y) = 5.607893e-05 m
D_z = -7.087755e-10 m, Sigma(z) = 5.205426e-05 m
---- FINAL SOLUTION ----
NCB1-3
x = -2660361.453 +/- 5.602523e-05 m del_e = 0.289 +/- 5.971492e-05 m
y = -3785656.3584 +/- 5.607893e-05 m del_n = -0.0543 +/- 6.148844e-05 m
z = 4373657.0134 +/- 5.205426e-05 m del_u = 0.7118 +/- 4.057857e-05 m
Lat. = 43.580175341632085 deg, Long. = -125.09754422671408, Hgt.msl = -1075.9443863983117 m
NCB1-2
x = -2660958.2954 +/- 5.602523e-05 m del_e = 0.289 +/- 5.971492e-05 m
y = -3786495.9768 +/- 5.607893e-05 m del_n = -0.0541 +/- 6.148821e-05 m
z = 4372517.1234 +/- 5.205426e-05 m del_u = 0.7118 +/- 4.057891e-05 m
Lat. = 43.56634902617636 deg, Long. = -125.09761313050257, Hgt.msl = -1115.3335735795908 m
NCB1-1
x = -2661822.702 +/- 5.602523e-05 m del_e = 0.2892 +/- 5.971572e-05 m
y = -3785267.4497 +/- 5.607893e-05 m del_n = -0.0541 +/- 6.148772e-05 m
z = 4373014.9389 +/- 5.205426e-05 m del_u = 0.7118 +/- 4.057848e-05 m
Lat. = 43.57274749300744 deg, Long. = -125.11511585101815, Hgt.msl = -1140.2321436132925 m
------------------------
There are 213 outliers found during this run.
This is 0.7% of the total number of data points.
Please re-run the program again with `--remove-outliers` flag to remove these outliers.
Exporting residuals to residuals.csv ...
Successfully exported data to ../data/exercise-2/output_2024/residuals.csv
Exporting outliers to outliers.csv ...
Successfully exported data to ../data/exercise-2/output_2024/outliers.csv
Exporting distance_from_center to dist_center.csv ...
Successfully exported data to ../data/exercise-2/output_2024/dist_center.csv
Exporting process_dataset to process_dataset.nc ...
Successfully exported data to ../data/exercise-2/output_2024/process_dataset.nc
Successfully exported data to ../data/exercise-2/output_2024/residuals.png
Successfully exported data to ../data/exercise-2/output_2024/residuals_enu_components.png
Finished GNATSS.
If you look closely at this new solution, you will notice that it hasn’t actually changed yet but in the residual plot a number of residuals have been highlighted as “To be removed”. This allows the user to preview flagged residuals in order to verify that nothing they want to keep will be accidentally removed.
The flagged residuals have been saved in a new output file called outliers.csv, which contains a record of time stamps and magnitudes for each flagged residual. This file is updated each time GNATSS is run, so if there are residuals flagged that you do not want to remove simply running GNATSS again with a different outlier threshold will update this file.
If the flagged residuals are satisfactory you can remove them from the inversion by running GNATSS with the --remove-outliers flag. This concatenates any residuals in the outliers.csv file to a new file, deletions.csv. GNATSS automatically reads the deletions file when run and will ignore any residuals with time stamps logged in that file. This is a non-destructive operation, so if for some reason you wish to include data previously removed back into the inversion you do not need to re-do any GNSS-A pre-processing; simply alter or remove the deletions file and you can re-interpret the residuals as needed.
Let’s go ahead and remove the residuals that have been flagged.
output = run_gnatss(
config_yaml,
distance_limit=distance_limit,
remove_outliers=True,
skip_posfilter=skip_posfilter
)Starting GNATSS ...
Gathering sound_speed at ../data/exercise-2/ctd/CTD_NCB1_Ch_Mi
Gathering gps_solution at ../data/exercise-2/output_2024/gps_solution.csv
Loading sound_speed from ../data/exercise-2/ctd/CTD_NCB1_Ch_Mi
Loading gps_solution from ../data/exercise-2/output_2024/gps_solution.csv
Loading deletions from ../data/exercise-2/output_2024/
Found /Users/johndesanto/Documents/gnatss-workshop/docs/../data/exercise-2/output_2024/outliers.csv file. Including into cuts...
Filtering out data outside of distance limit...
Pre-filtering data with fewer than 3 replies...
Computing harmonic mean...
lat=43.58017583 lon=-125.097547806 height=-1076.6562 alt=-1076.6562 internal_delay=0.08 sv_mean=1481.682 pxp_id='NCB1-3' azimuth=57.44 elevation=39.67
lat=43.566349513 lon=-125.097616709 height=-1116.0454 alt=-1116.0454 internal_delay=0.36 sv_mean=1481.695 pxp_id='NCB1-2' azimuth=-58.75 elevation=39.38
lat=43.57274798 lon=-125.115119432 height=-1140.9439 alt=-1140.9439 internal_delay=0.64 sv_mean=1481.704 pxp_id='NCB1-1' azimuth=-175.59 elevation=39.52
Finished computing harmonic mean
Preparing data inputs...
Perform solve...
--- 30308 epochs, 90924 measurements ---
After iteration: 1, rms residual = 121.11 cm, error factor = 69.832
NCB1-3
D_x = -7.482188e-02 m, Sigma(x) = 5.615829e-05 m
D_y = -6.268656e-01 m, Sigma(y) = 5.620100e-05 m
D_z = 4.544635e-01 m, Sigma(z) = 5.216588e-05 m
NCB1-2
D_x = -7.482188e-02 m, Sigma(x) = 5.615829e-05 m
D_y = -6.268656e-01 m, Sigma(y) = 5.620100e-05 m
D_z = 4.544635e-01 m, Sigma(z) = 5.216588e-05 m
NCB1-1
D_x = -7.482188e-02 m, Sigma(x) = 5.615829e-05 m
D_y = -6.268656e-01 m, Sigma(y) = 5.620100e-05 m
D_z = 4.544635e-01 m, Sigma(z) = 5.216588e-05 m
After iteration: 2, rms residual = 43.65 cm, error factor = 25.163
NCB1-3
D_x = -1.458872e-04 m, Sigma(x) = 5.613856e-05 m
D_y = 2.119286e-06 m, Sigma(y) = 5.619312e-05 m
D_z = 9.633313e-05 m, Sigma(z) = 5.215706e-05 m
NCB1-2
D_x = -1.458872e-04 m, Sigma(x) = 5.613856e-05 m
D_y = 2.119286e-06 m, Sigma(y) = 5.619312e-05 m
D_z = 9.633313e-05 m, Sigma(z) = 5.215706e-05 m
NCB1-1
D_x = -1.458872e-04 m, Sigma(x) = 5.613856e-05 m
D_y = 2.119286e-06 m, Sigma(y) = 5.619312e-05 m
D_z = 9.633313e-05 m, Sigma(z) = 5.215706e-05 m
After iteration: 3, rms residual = 43.65 cm, error factor = 25.163
NCB1-3
D_x = -3.193048e-10 m, Sigma(x) = 5.613856e-05 m
D_y = 6.11597e-10 m, Sigma(y) = 5.619311e-05 m
D_z = -2.902998e-10 m, Sigma(z) = 5.215706e-05 m
NCB1-2
D_x = -3.193037e-10 m, Sigma(x) = 5.613856e-05 m
D_y = 6.115983e-10 m, Sigma(y) = 5.619311e-05 m
D_z = -2.902997e-10 m, Sigma(z) = 5.215706e-05 m
NCB1-1
D_x = -3.193041e-10 m, Sigma(x) = 5.613856e-05 m
D_y = 6.115972e-10 m, Sigma(y) = 5.619311e-05 m
D_z = -2.903000e-10 m, Sigma(z) = 5.215706e-05 m
---- FINAL SOLUTION ----
NCB1-3
x = -2660361.4465 +/- 5.613856e-05 m del_e = 0.2991 +/- 5.983838e-05 m
y = -3785656.3666 +/- 5.619311e-05 m del_n = -0.054 +/- 6.161674e-05 m
z = 4373657.0166 +/- 5.215706e-05 m del_u = 0.7161 +/- 4.064859e-05 m
Lat. = 43.58017534394847 deg, Long. = -125.09754410234194, Hgt.msl = -1075.9400771835633 m
NCB1-2
x = -2660958.2889 +/- 5.613856e-05 m del_e = 0.2991 +/- 5.983839e-05 m
y = -3786495.985 +/- 5.619311e-05 m del_n = -0.0538 +/- 6.161651e-05 m
z = 4372517.1265 +/- 5.215706e-05 m del_u = 0.7161 +/- 4.064893e-05 m
Lat. = 43.56634902850221 deg, Long. = -125.09761300615807, Hgt.msl = -1115.3292644370006 m
NCB1-1
x = -2661822.6954 +/- 5.613856e-05 m del_e = 0.2993 +/- 5.983919e-05 m
y = -3785267.4579 +/- 5.619311e-05 m del_n = -0.0538 +/- 6.161602e-05 m
z = 4373014.9421 +/- 5.215706e-05 m del_u = 0.7161 +/- 4.064851e-05 m
Lat. = 43.57274749534801 deg, Long. = -125.11511572664887, Hgt.msl = -1140.227836665118 m
------------------------
There are 0 outliers found during this run.
Exporting residuals to residuals.csv ...
Successfully exported data to ../data/exercise-2/output_2024/residuals.csv
Exporting outliers to outliers.csv ...
No data found for outliers, skipping export!
Exporting distance_from_center to dist_center.csv ...
Successfully exported data to ../data/exercise-2/output_2024/dist_center.csv
Exporting process_dataset to process_dataset.nc ...
Successfully exported data to ../data/exercise-2/output_2024/process_dataset.nc
Successfully exported data to ../data/exercise-2/output_2024/residuals.png
Successfully exported data to ../data/exercise-2/output_2024/residuals_enu_components.png
Finished GNATSS.
Flagging Residuals with --residual-range-limit¶
We can see quite a bit more detail in the residuals now that we have removed some outliers. There is a pretty clear oceanograhpic signal, but believe it or not this is not necessarily a problem. A fundamental assumption of the GNATSS inversion is that the ocean behaves like a horizontally layered medium with any variation occuring in the upper water column above the thermocline. If this is the case, all of the oceanographic variation will occur in a section of the water column where the acoustic rays to each transponder are geometrically close to each other as long as we survey at the array center. The result is that each ping sees a similar oceanographic delay to each transponder and the residuals plot nearly on top of each other, and as long as this holds the oceanographic signal will be interpreted as a vertical offset rather than a horizontal offset. (This is why we can achieve cm-level horizontal positioning accuracy but worse or even unconstrained vertical positioning accuracy.)
However, there are also still some of the outliers from cycle slips and these outliers overlap in magnitude with the “main” residual time series at some points. Because of this, an outlier threshold is not appropriate to remove all of these points. An alternative tool in GNATSS for outlier flagging is to set a --residual-range-limit flag. This measures the difference between the minimum and maximum residual for each epoch, if the range is greater than the threshold that epoch is flagged. This routine takes advantage of the key point described above that the signals we want to remove are those where the residual time series for each transponder diverge, since that is evidence of our fundamental assumptions of the ocean structure starting to break down. Similar to the --residual-limit flag, the --residual-range-limit XXX flag takes a value in cm as input.
Let’s demonstrate this:
residual_range_limit = 40
output = run_gnatss(
config_yaml,
distance_limit=distance_limit,
residual_range_limit=residual_range_limit,
skip_posfilter=skip_posfilter
)Starting GNATSS ...
Gathering sound_speed at ../data/exercise-2/ctd/CTD_NCB1_Ch_Mi
Gathering gps_solution at ../data/exercise-2/output_2024/gps_solution.csv
Loading sound_speed from ../data/exercise-2/ctd/CTD_NCB1_Ch_Mi
Loading gps_solution from ../data/exercise-2/output_2024/gps_solution.csv
Loading deletions from ../data/exercise-2/output_2024/
Filtering out data outside of distance limit...
Pre-filtering data with fewer than 3 replies...
Computing harmonic mean...
lat=43.58017583 lon=-125.097547806 height=-1076.6562 alt=-1076.6562 internal_delay=0.08 sv_mean=1481.682 pxp_id='NCB1-3' azimuth=57.44 elevation=39.67
lat=43.566349513 lon=-125.097616709 height=-1116.0454 alt=-1116.0454 internal_delay=0.36 sv_mean=1481.695 pxp_id='NCB1-2' azimuth=-58.75 elevation=39.38
lat=43.57274798 lon=-125.115119432 height=-1140.9439 alt=-1140.9439 internal_delay=0.64 sv_mean=1481.704 pxp_id='NCB1-1' azimuth=-175.59 elevation=39.52
Finished computing harmonic mean
Preparing data inputs...
Perform solve...
--- 30308 epochs, 90924 measurements ---
After iteration: 1, rms residual = 121.11 cm, error factor = 69.832
NCB1-3
D_x = -7.482188e-02 m, Sigma(x) = 5.615829e-05 m
D_y = -6.268656e-01 m, Sigma(y) = 5.620100e-05 m
D_z = 4.544635e-01 m, Sigma(z) = 5.216588e-05 m
NCB1-2
D_x = -7.482188e-02 m, Sigma(x) = 5.615829e-05 m
D_y = -6.268656e-01 m, Sigma(y) = 5.620100e-05 m
D_z = 4.544635e-01 m, Sigma(z) = 5.216588e-05 m
NCB1-1
D_x = -7.482188e-02 m, Sigma(x) = 5.615829e-05 m
D_y = -6.268656e-01 m, Sigma(y) = 5.620100e-05 m
D_z = 4.544635e-01 m, Sigma(z) = 5.216588e-05 m
After iteration: 2, rms residual = 43.65 cm, error factor = 25.163
NCB1-3
D_x = -1.458872e-04 m, Sigma(x) = 5.613856e-05 m
D_y = 2.119286e-06 m, Sigma(y) = 5.619312e-05 m
D_z = 9.633313e-05 m, Sigma(z) = 5.215706e-05 m
NCB1-2
D_x = -1.458872e-04 m, Sigma(x) = 5.613856e-05 m
D_y = 2.119286e-06 m, Sigma(y) = 5.619312e-05 m
D_z = 9.633313e-05 m, Sigma(z) = 5.215706e-05 m
NCB1-1
D_x = -1.458872e-04 m, Sigma(x) = 5.613856e-05 m
D_y = 2.119286e-06 m, Sigma(y) = 5.619312e-05 m
D_z = 9.633313e-05 m, Sigma(z) = 5.215706e-05 m
After iteration: 3, rms residual = 43.65 cm, error factor = 25.163
NCB1-3
D_x = -3.193048e-10 m, Sigma(x) = 5.613856e-05 m
D_y = 6.11597e-10 m, Sigma(y) = 5.619311e-05 m
D_z = -2.902998e-10 m, Sigma(z) = 5.215706e-05 m
NCB1-2
D_x = -3.193037e-10 m, Sigma(x) = 5.613856e-05 m
D_y = 6.115983e-10 m, Sigma(y) = 5.619311e-05 m
D_z = -2.902997e-10 m, Sigma(z) = 5.215706e-05 m
NCB1-1
D_x = -3.193041e-10 m, Sigma(x) = 5.613856e-05 m
D_y = 6.115972e-10 m, Sigma(y) = 5.619311e-05 m
D_z = -2.903000e-10 m, Sigma(z) = 5.215706e-05 m
---- FINAL SOLUTION ----
NCB1-3
x = -2660361.4465 +/- 5.613856e-05 m del_e = 0.2991 +/- 5.983838e-05 m
y = -3785656.3666 +/- 5.619311e-05 m del_n = -0.054 +/- 6.161674e-05 m
z = 4373657.0166 +/- 5.215706e-05 m del_u = 0.7161 +/- 4.064859e-05 m
Lat. = 43.58017534394847 deg, Long. = -125.09754410234194, Hgt.msl = -1075.9400771835633 m
NCB1-2
x = -2660958.2889 +/- 5.613856e-05 m del_e = 0.2991 +/- 5.983839e-05 m
y = -3786495.985 +/- 5.619311e-05 m del_n = -0.0538 +/- 6.161651e-05 m
z = 4372517.1265 +/- 5.215706e-05 m del_u = 0.7161 +/- 4.064893e-05 m
Lat. = 43.56634902850221 deg, Long. = -125.09761300615807, Hgt.msl = -1115.3292644370006 m
NCB1-1
x = -2661822.6954 +/- 5.613856e-05 m del_e = 0.2993 +/- 5.983919e-05 m
y = -3785267.4579 +/- 5.619311e-05 m del_n = -0.0538 +/- 6.161602e-05 m
z = 4373014.9421 +/- 5.215706e-05 m del_u = 0.7161 +/- 4.064851e-05 m
Lat. = 43.57274749534801 deg, Long. = -125.11511572664887, Hgt.msl = -1140.227836665118 m
------------------------
There are 670 outliers found during this run.
This is 2.21% of the total number of data points.
Please re-run the program again with `--remove-outliers` flag to remove these outliers.
Exporting residuals to residuals.csv ...
Successfully exported data to ../data/exercise-2/output_2024/residuals.csv
Exporting outliers to outliers.csv ...
Successfully exported data to ../data/exercise-2/output_2024/outliers.csv
Exporting distance_from_center to dist_center.csv ...
Successfully exported data to ../data/exercise-2/output_2024/dist_center.csv
Exporting process_dataset to process_dataset.nc ...
Successfully exported data to ../data/exercise-2/output_2024/process_dataset.nc
Successfully exported data to ../data/exercise-2/output_2024/residuals.png
Successfully exported data to ../data/exercise-2/output_2024/residuals_enu_components.png
Finished GNATSS.
We can see that much of the cycle slips have now been flagged as well as some places where the time series diverged. Once again, we can remove these outliers with the --remove-outliers flag.
output = run_gnatss(
config_yaml,
distance_limit=distance_limit,
remove_outliers=True,
skip_posfilter=skip_posfilter
)Starting GNATSS ...
Gathering sound_speed at ../data/exercise-2/ctd/CTD_NCB1_Ch_Mi
Gathering gps_solution at ../data/exercise-2/output_2024/gps_solution.csv
Loading sound_speed from ../data/exercise-2/ctd/CTD_NCB1_Ch_Mi
Loading gps_solution from ../data/exercise-2/output_2024/gps_solution.csv
Loading deletions from ../data/exercise-2/output_2024/
Found /Users/johndesanto/Documents/gnatss-workshop/docs/../data/exercise-2/output_2024/outliers.csv file. Including into cuts...
Filtering out data outside of distance limit...
Pre-filtering data with fewer than 3 replies...
Computing harmonic mean...
lat=43.58017583 lon=-125.097547806 height=-1076.6562 alt=-1076.6562 internal_delay=0.08 sv_mean=1481.682 pxp_id='NCB1-3' azimuth=57.44 elevation=39.67
lat=43.566349513 lon=-125.097616709 height=-1116.0454 alt=-1116.0454 internal_delay=0.36 sv_mean=1481.695 pxp_id='NCB1-2' azimuth=-58.75 elevation=39.38
lat=43.57274798 lon=-125.115119432 height=-1140.9439 alt=-1140.9439 internal_delay=0.64 sv_mean=1481.704 pxp_id='NCB1-1' azimuth=-175.59 elevation=39.52
Finished computing harmonic mean
Preparing data inputs...
Perform solve...
--- 29638 epochs, 88914 measurements ---
After iteration: 1, rms residual = 120.62 cm, error factor = 69.542
NCB1-3
D_x = -7.165169e-02 m, Sigma(x) = 5.677350e-05 m
D_y = -6.269408e-01 m, Sigma(y) = 5.681506e-05 m
D_z = 4.537265e-01 m, Sigma(z) = 5.273603e-05 m
NCB1-2
D_x = -7.165169e-02 m, Sigma(x) = 5.677350e-05 m
D_y = -6.269408e-01 m, Sigma(y) = 5.681506e-05 m
D_z = 4.537265e-01 m, Sigma(z) = 5.273603e-05 m
NCB1-1
D_x = -7.165169e-02 m, Sigma(x) = 5.677350e-05 m
D_y = -6.269408e-01 m, Sigma(y) = 5.681506e-05 m
D_z = 4.537265e-01 m, Sigma(z) = 5.273603e-05 m
After iteration: 2, rms residual = 42.71 cm, error factor = 24.593
NCB1-3
D_x = -1.468406e-04 m, Sigma(x) = 5.675354e-05 m
D_y = 2.336929e-06 m, Sigma(y) = 5.680714e-05 m
D_z = 9.625314e-05 m, Sigma(z) = 5.272716e-05 m
NCB1-2
D_x = -1.468406e-04 m, Sigma(x) = 5.675354e-05 m
D_y = 2.336929e-06 m, Sigma(y) = 5.680714e-05 m
D_z = 9.625314e-05 m, Sigma(z) = 5.272716e-05 m
NCB1-1
D_x = -1.468406e-04 m, Sigma(x) = 5.675354e-05 m
D_y = 2.336929e-06 m, Sigma(y) = 5.680714e-05 m
D_z = 9.625314e-05 m, Sigma(z) = 5.272716e-05 m
After iteration: 3, rms residual = 42.71 cm, error factor = 24.593
NCB1-3
D_x = -1.195506e-10 m, Sigma(x) = 5.675354e-05 m
D_y = 3.175827e-10 m, Sigma(y) = 5.680713e-05 m
D_z = -1.803113e-10 m, Sigma(z) = 5.272716e-05 m
NCB1-2
D_x = -1.195515e-10 m, Sigma(x) = 5.675354e-05 m
D_y = 3.175828e-10 m, Sigma(y) = 5.680713e-05 m
D_z = -1.803108e-10 m, Sigma(z) = 5.272716e-05 m
NCB1-1
D_x = -1.195506e-10 m, Sigma(x) = 5.675354e-05 m
D_y = 3.175824e-10 m, Sigma(y) = 5.680713e-05 m
D_z = -1.803115e-10 m, Sigma(z) = 5.272716e-05 m
---- FINAL SOLUTION ----
NCB1-3
x = -2660361.4433 +/- 5.675354e-05 m del_e = 0.3017 +/- 6.049460e-05 m
y = -3785656.3667 +/- 5.680713e-05 m del_n = -0.0533 +/- 6.229238e-05 m
z = 4373657.0158 +/- 5.272716e-05 m del_u = 0.7143 +/- 4.108808e-05 m
Lat. = 43.58017535007013 deg, Long. = -125.09754406969948, Hgt.msl = -1075.9418608896467 m
NCB1-2
x = -2660958.2857 +/- 5.675354e-05 m del_e = 0.3017 +/- 6.049461e-05 m
y = -3786495.9851 +/- 5.680713e-05 m del_n = -0.0531 +/- 6.229215e-05 m
z = 4372517.1258 +/- 5.272716e-05 m del_u = 0.7144 +/- 4.108843e-05 m
Lat. = 43.56634903462006 deg, Long. = -125.0976129735229, Hgt.msl = -1115.3310483098217 m
NCB1-1
x = -2661822.6923 +/- 5.675354e-05 m del_e = 0.3019 +/- 6.049541e-05 m
y = -3785267.458 +/- 5.680713e-05 m del_n = -0.0532 +/- 6.229166e-05 m
z = 4373014.9413 +/- 5.272716e-05 m del_u = 0.7143 +/- 4.108799e-05 m
Lat. = 43.572747501472655 deg, Long. = -125.11511569401678, Hgt.msl = -1140.229621044797 m
------------------------
There are 0 outliers found during this run.
Exporting residuals to residuals.csv ...
Successfully exported data to ../data/exercise-2/output_2024/residuals.csv
Exporting outliers to outliers.csv ...
No data found for outliers, skipping export!
Exporting distance_from_center to dist_center.csv ...
Successfully exported data to ../data/exercise-2/output_2024/dist_center.csv
Exporting process_dataset to process_dataset.nc ...
Successfully exported data to ../data/exercise-2/output_2024/process_dataset.nc
Successfully exported data to ../data/exercise-2/output_2024/residuals.png
Successfully exported data to ../data/exercise-2/output_2024/residuals_enu_components.png
Finished GNATSS.
We’ve now made significant progress with this particular data set. In general, I prefer to flag to a standard where the outlier threshold is such that it borders the main residual time series (this varies according to the individual data set) and the residual range threshold is 20 cm.
However, it is worth noting that the mean of the residuals changes as you remove outliers, particularly large outliers. This shifts the relative position of “0 cm”, so I prefer to flag conservatively over many iterations as I approach these standards in order to minimize the risk of removing good data from the inversion.
Further Exercises¶
- Continue flagging the data in the 2024 NCB1 survey until you are satisfied with the results.
output = run_gnatss(
config_yaml,
distance_limit=distance_limit,
residual_limit=,
residual_range_limit=,
skip_posfilter=skip_posfilter
)- There are four other data sets collected at this array between the years 2018-2023. Try to process these data and flag them on your own. Like with the 2024 survey, each of these data sets has already been preprocessed so you only have to run the GNATSS solver module.
config_yaml = "../data/exercise-2/NCB1_2023_config.yaml"
output = run_gnatss(
config_yaml,
distance_limit=distance_limit,
residual_limit=,
residual_range_limit=,
skip_posfilter=skip_posfilter
)