Satellite Navigation-derived Instantaneous Velocities
Brendan Crowell Pacific Northwest Seismic Network University of Washington
If you use this code for your research, please cite either manuscript
Crowell, B. W. (2021), Near-field strong ground motions from GPS-derived velocities for 2020 Intermountain Western United States Earthquakes, Seismological Research Letters, doi: 10.1785/0220200325.
Crowell, B., DeGrande, J., Dittmann, T., & Ghent, J. (2023). Validation of Peak Ground Velocities Recorded on Very-high rate GNSS Against NGA-West2 Ground Motion Models. Seismica, 2(1). https://doi.org/10.26443/seismica.v2i1.239
This package computes GNSS velocities using the broadcast ephemeris and high-rate RINEX files. This loosely follows the method of Colosimo et al. (2011). It uses both the L1 and L2 phase observables and inverts for velocity on an epoch by epoch basis with the narrow lane combination. There is no single frequency mode currently (i.e. Grapenthin et al. (2018)), however I have tested it using just L1 and it produces almost identical results.
Dependencies
You will need to download crx2rnx and teqc and place both in the main working directory where SNIVEL is located. In some distributions of crx2rnx, the executable is capitalized. Simply rename this file by doing:
mv CRX2RNX > crx2rnx
The following python3 packages are required:
georinex, numpy, scipy, obspy, and urllib
You should be able to get all of these through pip, brew, or conda depending on your Python install.
One complexity with georinex is that it is not available through conda directly. You need to install pip on conda and then install georinex by calling the anaconda version of pip3 such that
/path/to/anaconda/bin/pip3 install georinex
Running the code
The logic flow is rather simple:
(1) in the file sites_process.txt, you simply need to add 4-character station ids in lowercase. If the station is not located in UNAVCO or CWU archives, you will need to put your files into a directory called 'rinex_hr'. This folder is autogenerated, so run SNIVEL.py once and you will get the directory. When you put your RINEX file in this folder, ensure that it follows standard naming convention, "SITEDOY0.YRo". If you use uppercase for SITE, enter it that way in the sites_process.txt file. If your RINEX is Hatanaka compressed, use the crx2rnx command to uncompress it (i.e. go from YRd to YRo). For files downloaded, the uncompression is performed automatically.
(2) in the file dates_process.txt, you enter in the dates and times you want to process. There are 4 columns in this file:
YEAR DOY START_TIME MINUTES
YEAR should always be 4 digits.
DOY is the day of year, three characters between 001 and 365/366 for leap years. Check the NGS GPS Calendar to determine your DOY.
START_TIME is the time in HR:MN:SC format that you want to start processing. Note that RINEX files use GPS time, so there is currently an 18 second difference between UTC time (GPS-UTC=18) as of Feb 2020.
MINUTES is simply the number of minutes you wish to process. Decimals are allowed here.
(3) edit the top line of SNIVEL.py with your environment path to python. Or you can just run "python3 SNIVEL.py"
(4) The code will do the following things once it is started:
First, it will download the broadcast orbit file for the first date given. This file will be put into a directory called nav/.
Next it will try to download the RINEX file for the first site at the first date given. If the RINEX file is already in the rinex_hr/ directory, it will skip this step.
Georinex will then read the RINEX header to determine the a priori position of the station. If this is terribly off, you should manually change this in the RINEX file and rerun the code. I haven't tested how bad it needs to be before this is a problem, but many 10s of meters would not be great.
After this, the RINEX file is trimmed to remove non-GPS satellites and to the START_TIME and MINUTES specified in the dates_process.txt file. In subsequent versions, I will add capabilities for GLONASS, Galileo and Beidou, but for now, GPS is sufficient.
Next, from the broadcast ephemeris file, it reads the Klobuchar ionospheric corrections.
At this point, a file is created in the output/ folder that will host all of the observables to perform the velocity calculation. The file naming is output/observables_SITE_DOY_YEAR.txt. The columns of this file are as follows:
index, gps time, gps week, gps second of week, GPS satellite number, x direction cosine, y direction cosine, z direction cosine, L1 (m) with sat clock correction, L2 (m) with sat clock correction, azimuth, elevation angle, L1 ionosphere correction, L2 ionosphere correction, Niell slant dry troposphere delay, dx between receiver and satellite, dy between receiver and satellite, dz between receiver and satellite, Niell wet tropospheric delay mapping function
While this file is being created, it is finding the proper locations of the satellites at the broadcast time, computing the Niell hydrostatic tropospheric delay, computing the relativistic clock correction, and computing the Klobuchar ionospheric corrections.
Once the observables file is created, the velocity calculations begin. Starting from the second epoch, it computes the time difference in the satellite locations (from i to i-1) and computes the time difference in the narrow-lane combination. L1 and L2 observations are pre-corrected with Klobuchar and slant tropospheric delay values. Every epoch, a new Green's function matrix is created that contains the direction cosines and a single value. A single row looks like
[rx, ry, rz, 1]
where the 1 allows for drifts in the clock bias. It then simply performs a bounded least squares problem where the dLC values are the data vector and we are solving for a vector of dx/dt, dy/dt, dz/dt, dclock/dt. The dxyz values are then converted to dn/dt, de/dt, du/dt and output to the file output/velocities_SITE_DOY_YEAR.txt.
Example
After you have installed dependencies, if you run this script directly out of the box by typing:
python3 SNIVEL.py
it will compute the velocities at the nearest sites during the Mw 7.1 Ridgecrest, CA earthquake on July 6, 2019. If you plot the values for station P595, you should get something that looks like:
Thats it, there is nothing else to really do, simply change the site list and the dates file and start computing some velocities!