Skip to content

Program structure and algorithm

Jump to bottom Edit New page
wangyulu1999 edited this page Feb 12, 2025 · 3 revisions

Program structure

PRIDE PPP-AR software runs according the structure shown in Figure 1, the process procedures are divided into three modules, least-squares estimator and integer ambiguity resolution in addition to a data preparation and pre-processing module. The first part, data preparation and pre-processing, prepare table file and precise products for following data process. The spp (standard point positioning) module will be used in this part to calculate the prior positions of station. The function of sp3orb (SP3 orbit) is to transform SP3 orbit into a self-defined binary format. Then, the software can efficiently access the precise orbit products. If multipath correction is performed, mhm will generate the corresponding "mhm-file” to provide multipath correction values in the least-squares section. In least-squares estimator part, tedit is used to make data tentative pre-processing and generate “log-file” to record the RINEX (The Receiver Independent Exchange Format) health diagnosis information. Once got the “log-file”, lsq (least-square adjustment) module can realize parameter estimation and output results. Then used redig (a posteriori residual diagnosis) module, the residuals can be processed and new “log-file” can be generated. By the iteration of lsq and redig, data cleaning is completed. If ambiguity is not fixed, the ambiguity-float solution can be obtained. Otherwise, the module named arsig (ambiguity resolution at a single receiver) will be used to realize wide-lane and narrow-lane ambiguities resolution. In the next round lsq processing, these integer ambiguities will be introduced as hard constraints to achieve ambiguity-fixed solutions.

image9 Figure 1 Program structure of PRIDE PPP-AR

Modules of PRIDE PPP-AR

The functions and usages of each module of PRIDE PPP-AR are shown below.

  • spp is used to calculate initial positions of station, if the positioning mode is “S”, the initial position will output to the “sit.xyz” file. Otherwise, if the positioning mode is “K/P/L”, spp will also generate “kin_file” to record coordinates time series. Note: If the positioning mode is “F”, the station coordinates will be fixed to the IGS daily solution.

    ① S/F Mode:

    spp -elev 10 -trop saas -ts ${ts} -te ${te} -ti ${interval} ${rinexobs} ${rinexnav}

    ②K/P Mode:

    spp -elev 10 -trop saas -ts ${ts} -te ${te} -ti ${interval} -o kin_${ydoy_s[0]}${ydoy_s[1]}_${site} ${rinexobs} ${rinexnav}

    ③ L Mode:

    spp -elev 0 -trop non -ts ${ts} -te ${te} -ti ${interval} -o kin_${ydoy_s[0]}${ydoy_s[1]}_${site} ${rinexobs} ${rinexnav}

    Where, (Note: Symbol ${} denotes taking a value for the variable).

    Command Description
    [-elev] is optional, up to the elevation angle (unit: °)
    [-trop] is optional, Tropsphere correction model. NON: Not correct; saas: saastamoinen model
    [-ts] is optional, start time(format: year/month/day hour: minute:second)
    [-te] is optional, end time(format: year/month/day hour: minute:second)
    [-ti] is optional, sampleing rate (Unit: s)
    [-o] is optional, output file. Not specified for static calculation; specified for dynamic calculation, and the results are output to the specified file.
    ${rinexobs} is RINEX observation file
    ${rinexnav} is RINEX Broadcast emphemeris file
    [-?/-h] is optional, output information for help

    If you want to call the SPP module separately, you can input the command based on the following example.****

    Example: (Note: The folder must contain the corresponding O and N files, and the ts and te parameters should be within the time range of the observation files.)

    ① S/F Mode:

    Input (command line):

    spp -elev 10 -trop saas -ts 2023/06/30 00:00:00 -te 2023/06/30 23:59:59 -ti 30 abmf1810.23o brdm1810.23p

    Output (command line):

    Position :     2919786.1342   -5383745.6171    1774604.6673
Duration : 2023 06 30 00 00  0.00    86370.00

    ② K/P Mode:

    Input (command line):

    spp -elev 10 -trop saas -ts 2023/06/30 00:00:00 -te 2023/06/30 23:59:59 -ti 30 abmf1810.23o brdm1810.23p -o kin_2023181_abmf

    Output:

    The file kin_2023181_abmf contains the calculated results for all epochs.

    ③ L Mode: The same as in ② K/P mode.

  • otl is used to calculate the tide correction component of the station. The calculation code used in this module is the method from the website https://www.zcyphygeodesy.com/ provided by researcher ChuanYin Zhang, which is modularized into a separate executable program here. The usage is as follows. (Note: The results generated by this module are in binary files.)

    otl -b ${lat} -l ${lon} -h ${hgt} -s ${ts} -e ${te} -i ${ti} -o ${outfile}

    Where, (Note: Symbol ${} denotes taking a value for the variable).

    Command Description
    [-b] station latitude (unit: °)
    [-l] station longitude (unit: °)
    [-h] station ellipsoid height (unit: m)
    [-ts] start time(format:year/month/day hour:minute:second)
    [-te] end time(format:year/month/day hour:minute:second)
    [-ti] interval (Unit: s)
    ${outfile} output file name

  • sp3orb transforms SP3 orbit files into a self-defined binary format. Then, the software can efficiently access the precise orbit products. In addition, the reference frame is changed from an earth-fixed system into an inertial system through the ERP file.

    sp3orb ${sp3} -cfg ${ctrl_file}

    Where, (Note: Symbol ${} denotes taking a value for the variable).

    Command Description
    ${sp3} is SP3 file
    -cfg ${ctrl_file} is configuration file

  • tedit module is used to detect cycle jumps and receiver clock jumps, remove short piece of data and identify data gaps. This module can release the “log-file”, which can reflect the health condition of RINEX file. It should be noted that we do not repair cycle slips but just mark and record them in the ‘log’ file, for further processing.

    ① S/F Mode:

    tedit ${rinexobs} -time ${ymd[*]} ${hms[*]} -len ${session} -int ${interval} -xyz ${xyz[*]} -short 1200 -lc_check only -rhd ${rhd_file} -pc_check 300 -elev ${cutoff_elev} -rnxn ${rinexnav} -freq ${freq_cmb} -trunc_dbd ${tct_opt}

    ② K/P/L Mode:

    tedit ${rinexobs} -time ${ymd[*]} ${hms[*]} -len ${session} -int ${interval} -xyz kin_${year}${doy}_${site} -short 120 -lc_check no -elev ${cutoff_elev} -rhd ${rhd_file} -rnxn ${rinexnav} -freq ${freq_cmb} -trunc_dbd ${tct_opt}

    Where, (Note: Symbol ${} denotes taking a value for the variable).

    Command Description
    ${rinexobs} is RINEX observation file
    ${rinexnav} is RINEX broadcast emphemeris file
    [-time] is optional, strat time(format:year/month/day/hour/minute/second)
    [-len] is optional, session time of data processing (Unit:s)
    [-int] is optional, data processing interval (Unit:s)
    -xyz S/F Mode:Initial coordinates; P/K/L Mode:Initial ‘kin_*’ file
    [-short] is optional, data piece shorter than this value will be cut
    -lc_check is combination of cycle slips detection methods
    -pc_check checking the consistency of the receiver clock difference
    [-elev] is optional, cut-off elevation (Unit: °). The default value is 0
    [-rhd] is optional, the output RHD file
    [-freq] is optional, choose frequency combination. Default is “G12 R12 E15 C26 J12”
    [-trunc_dbd] is optional, whether to truncate ambiguities at the day-boundary. Options include "yes" or "no," with the default value being "no

  • mhm is used to construct the MHM model. Firstly, determine whether the “res-file” in the mhm folder was generated through F/S mode and whether it belongs to the same station. Secondly, divide the residuals of each satellite in each “res-file” into a 1 °× 1 ° grid according to their corresponding elevation and azimuth angles, and construct a 90 × 360 grid. Finally, perform quality control on the residuals in each 1 °× 1 ° grid, remove coarse errors, and output the mean value after removal as the multipath correction value for that grid to the “mhm-file”.

    mhm ${resdir}

    Where, (Note: Symbol ${} denotes taking a value for the variable).

    ${resdir} is residual folder for modeling

  • lsq is the estimator based on least-squares principle. In this module, we use the ionosphere-free combination to eliminate the first-order ionospheric delay. The module of lsq is used to process raw measurements and estimate unknown parameters such as positions, receiver clocks, tropospheric delays and ambiguities. It is the core part of PRIDE PPP-AR, after least-squares adjustment, we can obtain “pos-file” (or “kin-file”), “rck-file”, “ztd-file”, “htg-file”, “res-file” and “amb-file”.

    lsq ${config} ${rinexobs}

    Where, (Note: Symbol ${} denotes taking a value for the variable).

    Command Description
    ${config} is configuration file for current project
    ${rinexobs} is RINEX observation file

  • redig module can make new ambiguities according to the jumps within epochs. Besides, it can remove the huge residuals and remove the short periods. redig module updates “log-file” by reading “res-file” from lsq module. To clean the observations, the lsq and redig need to be iterated continuously.

    redig res_${year}${doy} -jmp $jump -sht $short -pce

    Where, (Note: Symbol ${} denotes taking a value for the variable).

    Command Description
    res_${year}${doy} is residual file generated by lsq
    [-jmp $jump] is optional, if the difference of residuals between adjacent epochs is greater than $jump, it is marked as a new ambiguity
    [-sht $short] is optional, if the effective time of ambiguity is less than $short, the corresponding observation data will be deleted
    [-pce] is optional, quality control for satellite’s data

  • arsig module is used to realize wide-lane and narrow-lane ambiguity fixed. Firstly, the float wide-lane ambiguities generated from lsq module are fixed to the nearest integer using a bootstrapping method. Secondly, used ionosphere-free combination ambiguities and fixed wide-land ambiguities, float narrow-lane ambiguities can be fixed by LAMBDA (Least-squares Ambiguity Decorrelation Adjustment) method in short observation sessions or be fixed directly to their nearest integers like wide-lane. Thirdly, with the fixed wide-lane and narrow-lane ambiguities (these constraint information is written in “cst-file”), users use lsq module to realize ambiguities fixed PPP solutions.

    arsig ${config}

    Where, (Note: Symbol ${} denotes taking a value for the variable).

    Command Description
    ${conifg} is configuration file for current project

pdp3 batch script

pdp3 is the batch script for PPP processing by PRIDE PPP-AR, which automatically processes GNSS data according to command line parameters. The user needs to ensure that the command line parameters are entered correctly and the configuration file is modified. pdp3 contains the information of processing procedures, you can read the script for more details.

The three basic configuration parameters in the "Basic Setting" section of the script can be modified by users as needed.

  1. Debugging Settings

    Setting Setting
    DEBUG=NO Should results files be kept in case of program execution failure.
    YES Results files are kept in case of program execution failure.
    NO(Defult) Results files are deleted in case of program execution failure.

  2. Offline Settings

    Setting Description
    OFFLINE=NO Should requests for updating/downloading files be skipped to save processing time in offline mode.
    YES Skip requests, automatic downloading of broadcast ephemeris and precise products, etc., is not possible.
    NO(Defult) Do not skip requests, automatic downloading of broadcast ephemeris and precise products, etc., is possible.

  3. Product Settings

    Setting Description
    USECACHE=YES Whether to use precision products and some table files under the local path
    YES(Defult) If there are corresponding files in the product directory/table directory, copy them to the working directory, otherwise, the corresponding files will be downloaded.
    NO Download the corresponding file directly and do not match under the local path.

  • The main() function is the entry of the script, and its flow is as follows:

    1. Analyze command line parameters, check executable programs and required system tools;
    2. Define and initialize variables;
    3. Output configuration information to the screen;
    4. Determine the project directory structure based on the input time span;
    5. Determine whether to call the PrepareMhmModel() function to generate a “mhm-file” based on the settings;
    6. Call the ProcessSingleSession() function based on input parameters to start the computation.

  • The PrepareMhmModel() function is used to generate a “mhm-file” for the station, and its flow is as follows:

    1. Obtain the satellite constellation and determine the modeling period based on the configuration file;
    2. Check if the “res-file” for the modeling period exists; if it does not exist, use the Rinex OBS file and configuration file of the modeling period to calculate the “res-file” in S mode;
    3. Call the mhm module to calculate the “mhm-file”.

  • The ProcessSingleSession() function is used to handle the observation data for a single session, and its flow is as follows:

    1. Initialization, including variable definition and assignment, copying configuration files to the current directory, etc;
    2. Call PrepareTables() function to prepare the required table files; call PrepareRinexNav() function to prepare the broadcast ephemeris; call PrePareProducts() to prepare the required precision products;
    3. If the positioning mode is for low Earth orbit satellite tracking, then prepare the necessary low Earth orbit satellite antenna data;
    4. Call sp3orb to convert the orbit product to binary;
    5. Call ProcessSingleSite() for single station data processing.

  • PrepareTables() for preparing table files, and its flow is as follows:

    1. Link the local files in the tale directory to the current directory;

    2. Prepare the leap seconds file “leap.sec”;

      (1) Checking whether the leap seconds file in the current directory accompanies the software (marked with * in the first line);

      (2) Download the leap seconds file if it does not match or if the leap seconds file does not exist in the current directory;

      (3) If the download is failed, copy the leap seconds file from the table directory to the current directory;

      (4) If the download is successful and does not match the leap seconds file in the table directory, replace the leap seconds file in the table directory with the downloaded leap seconds file.

    3. Prepare satellite parameters file “sat_parameters”.

      (1) Check the time lag of the satellite parameter files in the current directory;

      (2) Download satellite parameter files if it's lagging or does not exist in the current directory;

      (3) If the download is failed, copy the satellite parameter file from the table directory to the current directory;

      (4) If the download is successful and does not match the satellite parameter file in the table directory, replace the satellite parameter file in the table directory with the downloaded satellite parameter file.

  • PrepareRinexNav() for preparing broadcast ephemeris:

    1. If processing today's data, download and merge the hourly GPS broadcast ephemeris and hourly GLONASS broadcast ephemeris;

    2. If there is no broadcast ephemeris in short naming format in the data directory, match the broadcast ephemeris in long naming format starting with “BRDC00IGS_R_”, “BRDC00IGN_R_” and “BRDM00DLR_S_”;

    3. If neither short-name format nor long-name format is available, then download the broadcast ephemeris;

      (1) Downloading the multi-system broadcast ephemeris if it is after 2016;

      (2) Failure of multi-system broadcast ephemeris downloads or processing of pre-2016 data, downloading of GPS and GLONASS broadcast ephemeris and merging them.

    4. Checking that the required satellite systems are available in the broadcast ephemeris.

  • PrepareProducts() for preparing products, and its flow is as follows:

    1. Determine the directory where the precise products are located according to the “Product directory” in the configuration file, if it is not modified, i.e. Default, create a “product” directory under the year directory, where the precise products are located in the “common” subdirectory, the VMF1/VMF3 required grid files are located in the “vmf” subdirectory, the ionosphere grid files are located in the “ion” subdirectory, the SINEX files are located in the “ssc” subdirectory, and the low earth orbit satellite products are located in the “leo” subdirectory;

    2. Prepare precise orbit, precise clock error and ERP products, and merge them into one file if there are multiple files;

      (1) Copy the corresponding product from the products directory to the current directory if it is not the default product;

      (2) If it is the default product and there are corresponding files in the products directory, copy them to the current directory;

      (3) If it is the default product and there are no corresponding files in the products directory, download the corresponding files; the default download of Wuhan University Rapid Product (WUM0MGXRAP) product after 001 days in 2020; For 2019 and earlier, the default download is the IGS Third Reprocessing Combined Product (IGS2R03FIN);

      (4) If the Wuhan University Rapid Product (WUM0MGXRAP) is not available, then download the Wuhan University Real-time Archive Product (WUM0MGXRTS). The update frequency for the Wuhan University Real-time Archive Product is 3 hours, which is shorter than the 24 hours for the Rapid Product, but it has slightly lower accuracy and does not currently provide quaternion products and phase bias products.

    3. Prepare quaternions products and phase bias products, and merge them into one file if there are multiple files;

      (1) Copy the corresponding product from the products directory to the current directory if it is not the default product;

      (2) If it is the default product and there are corresponding files in the products directory, copy them to the current directory;

      (3) If it is the default product and there are no corresponding files in the products directory, download the corresponding files; the default download of Wuhan University Rapid Product (WUM0MGXRAP) product after 001 days in 2020; For 2019 and earlier, the default download is the IGS Third Reprocessing Combined Product (IGS2R03FIN).

    4. Check Low Earth Orbit Attitude Files (only available in positioning mode L);

      (1) If you are not using default products, copy the corresponding products from the product directory to the current directory;

      (2) If you are using default products and there are corresponding files in the product directory, copy them to the current directory;

      (3) If there are no corresponding files when processing GRACE/GRACE-FO, you need to manually invoke prepare_leodata.sh to download the low Earth orbit products.

    5. Prepare antenna correction files. (the Antenna Exchange Format, ANTEX);

      (1) The antenna correction file is related to the precise product being used;

      (2) If using standard precise products, use the ANTEX file defined in the "SYS / PCVS APPLIED" section of the clock bias header file;

      (3) If using CODE products (COD0MGXRAP/COM) in IGS14 framework, use M14.ATX or M20.ATX;

      (4) If none of them of (2) or (3), the latest IGS ANTEX file will be used by default.

    6. Copy or download the solution file if the positioning mode is F; copy or download the IONEX file if higher-order ionosphere correction is performed; download the atmospheric corresponding grid file if the VMF1/VMF3 mapping function is used.

  • The process of ProcessSingleSite() function is used to process data for a single station, and its flow is as follows:

    1. Initialization, including variable definition and assignment, getting configuration options, etc.;
    2. Call the ComputeInitialPos() function to calculate an initial coordinate and output the obtained coordinates to sit.xyz.
    3. If the positioning mode is not L, call the xyz2blh() function to check if the h-coordinate of the sit.xyz result for the station is within the range of -4km to 20km from the Earth's surface, which is the normal ground station coordinate range.
    4. Invoke tedit for data pre-processing based on the positioning mode and related configuration information
    5. Call lsq and redig iterations for residual editing to identify residual cycle slips until no new ambiguities and observations are censored
    6. If the command-line parameters are not specified for float solutions and signal bias products exist, call arsig to fix ambiguity, and then call lsq to adjust again; otherwise, end the calculation to get only the ambiguity-float solution.

In addition, for users who can't solve online, you can refer to the PrepareTables(), PrepareRinexNav() and PrepareProducts(), which has the download address of those external files. You can pick out these three functions and modify them as scripts for downloading files. Users can download and place them in the corresponding directory. See Appendix A for a brief description of these required external files.

It should be noted that the default precise products start with “WUMMGXRAP_”, in which the satellite clock products are located in the “${year}/clock” directory, the bias products are located in the “${year}/bias” directory, and the satellite ephemeris, attitude file and ERP file are located in the “${year}/orbit” directory; in case of high-order ionospheric correction, the IONEX file needs to be downloaded separately; if the mapping function is VM1/VM3, the tropospheric grid files of the current day and the hour before and after need to be downloaded for interpolation; the ANTEX file matched with the default precise product is recorded in the satellite clock product file header. If other products are used, the latest IGS14 ANTEX file at the time of software release is used by default. The SINEX file needs to be downloaded in F mode.

Algorithms for each module

spp

spp is used to compute the initial coordinates of a station and is the first module called by PRIDE PPP-AR software for data processing. It is an essential step for unfolding data preprocessing tasks. spp can calculate the initial coordinates of a station based on user-provided observation files and broadcast ephemeris files, and it can output the calculation results to the command line or in the form of files. When using this module, if there are multiple days of observation files for the same station and corresponding broadcast ephemeris files in the same directory, the software will automatically perform standard single-point positioning for multiple days. The output result files will include results for multiple days.

Regarding the algorithm implementation of the spp module, it is adapted from the open-source software RTKLIB. If users are interested in the algorithm implementation, they can visit the following link: https://github.com/tomojitakasu/RTKLIB to learn more details. Here, we extend our heartfelt thanks to the author, Tomoji Takasu, for their contributions.

otl

otl is used to calculate the tide correction component of the station. Tide correction is necessary to be considered for high-precision PPP positioning solution, and the tide correction has a great impact on the coastal station, and its component can be up to 5cm in U direction, so it needs to be eliminated in the solution.

There are two methods for tide correction calculation in the software. One is the conventional Scherneck method (the model is FES2004), but this method requires that the tide correction coefficients for the corresponding station coordinates are available in the oceanload file under the table directory. If users want to continue to use the conventional Scherneck's method to calculate the tide correction, they need to refer to section General operation steps to get the tide correction coefficients of the stations (they need to submit the coordinates on the web site, and then get the tide correction coefficients by email), and then add the obtained coefficients into the oceanload file under the table directory, which is a cumbersome operation; In another way we borrowed the tide component calculation method from the geophysical geodesy large-scale scientific computing platform (www.zcyphygeodesy.com) developed by researcher ChuanYin Zhang of the China Academy of Surveying and Mapping Sciences.

mhm

mhm is used to calculate the multipath correction values of stations. Multipath delay is one of the important error sources of high-precision PPP, which is difficult to eliminate or weaken through differential methods or observation value combinations, and restricts the precision of multi-system GNSS precision positioning.

The software currently uses the multipath hemispherical map to compensate multipath delay, which is mainly divided into two parts: first, check whether the header information of the “res-file” used for modeling meets the requirements, then obtain the signal frequency of each constellation and determine whether it belongs to overlap-frequency signal, and then divide the satellite residuals into a 360 ° × 90 ° grid according to azimuth and elevation angles. Next, calculate the average and standard deviation of each 1 °× 1 ° grid, and perform quality control to eliminate errors. Finally, calculate the RMS and mean STD of the entire model as the model indicator output.

image10

Figure 2 mhm algorithm flow.

tedit

As the data pre-processing module of PRIDE PPP-AR, tedit is mainly used to check the original observation file and detect the cycle slips and outliers. The main process can be divided into two parts: firstly, check the original observation file and construct the test volume, then identify the outliers in the test volume to determine whether the cycle slips occurs. Finally, the check information is output to the log-file. The data processing flow is shown in the figure below

image11

Figure 3 tedit algorithm flow. (The light blue columns indicate the constructed test volume; the blue line indicates the input of each test volume; the black dashed line indicates the running process controlled by the pdp3 script and the process is performed only once; SDBS indicates the single difference between satellites)

  1. Check the observation file and construct the test volume

    Read the original dual-frequency observations on an epoch-by-epoch basis. Calculate satellite elevation, distance between station and satellite, and satellite clock error based on the broadcast ephemeris. Then construct the geometry-free combination $L_G=\frac{L_1-L_2}{\lambda_2-\lambda_1}$. In this process, the $L_G$ is initially checked for ionospheric anomalies and cycle slips are detected. Then remove data based on data availability and satellite elevation, etc.

    Due to the instability of the receiver clock, the receiver clock error needs to be checked in advance, including the clock jump detection and gross error rejection. Construct the receiver clock jump test volume as follows:

$$ \begin{cases} R_{i,k}=(P_{i,k}-P_{i,k-1})-(L_{i,k}-L_{i,k-1}) \tag{1} \end{cases} $$

where, $i=1, 2$ denotes the frequency number; $k$ and $k-1$ denotes the adjacent non-rejected current and previous epochs. The magnitude of the clock jump check for all satellites in the current epoch is used to determine whether a receiver clock jump occurs in that epoch, and if it exists, the receiver clock jump correction is calculated to the original observation and recorded for subsequent data processing.

The receiver clock error check is required in S/F mode. The receiver clock error test volume is constructed in the following form:

$$ \begin{cases} \bar P_{0,k}=P_{0,k}-(\rho-c\mathrm{\Delta}T^s) \tag{2} \end{cases} $$

where, $\rho$ and $\mathrm{\Delta}T^s$ are the distance between station and satellite and the satellite clock error calculated based on the broadcast ephemeris, respectively. This test volume eliminates the geometric distance, satellite clock error, and ionospheric delay, and leaves only the receiver clock error and multi-path effects for pseudo-range, and is therefore used to check for gross error in the receiver clock error. The weighted mean value of this test volume is obtained by median-based robust estimation and compared with a given threshold value to locate and reject the gross error in receiver clock error.

The M-W combination and the phase ionosphere-free combination are constructed, and the cycle slips are initially detected based on the change rate of the M-W combination of adjacent epoch. Based on the results of the above data rejection and cycle slips detection, short piece of data or big gap of adjacent non-rejected ephemerides are identified and marked for subsequent examination.

  1. Cycle slips detection and data rejection according to the test volume

    The pdp3 script passes different control parameters to tedit according to the positioning mode, and different combinations of methods are used in tedit for data preprocessing, with only SDBS polynomial fit in S/F mode, and only M-W combination tests in K mode.

    For the SDBS polynomial fit, the fitted values are chosen from the difference between the epoch results to facilitate better positioning of the cycle slips and gross error. The process is mainly divided into two parts: firstly, we fit each satellite LG combination, calculate its fitted residuals and RMS(Root Mean Square), mark the satellites that cannot be successfully fitted, count the number of available epoch for each satellite in the fitted arc, and take the satellite with the highest number of available epoch as the reference satellite. Next, the SDBS fitting term is constructed as follows:

$$ \begin{cases} L_{c,k}=\left(L_{0,k}^{si}-\left(\rho^{si}-c\mathrm{\Delta}T^{si}\right)\right)-(L_{0,k}^{sr}-(\rho^{sr}-c\mathrm{\Delta}T^{sr})) \tag{3} \end{cases} $$

where, $si$ and $sr$ indicates the current satellite and the reference satellite, respectively. The test volume eliminates the effects of receiver clock error, satellite clock error and geometric distance, leaving only tropospheric delay, ambiguity, multi-path effect and observation noise, etc. The difference between the epoch can eliminate ambiguity and weaken tropospheric delay. The residuals and RMS are calculated by fitting the difference between the epoch results, and the residuals are used to determine whether the cycle slip occurs and mark it. Finally, the fitting results are statistically calculated: (1) to determine whether the current satellite has a cycle slip based on the LC combination SDBS fitting results; (2) to determine whether the reference satellite has a cycle slip based on the LG combination fitting results and the LC combination fitting results.

When testing the M-W combination, the mean and variance of the M-W combination are calculated recursively, and the M-W combination is marked according to the magnitude of the difference between the current epoch and the mean. tedit also has a polynomial fitting process for the LG combination to determine whether a cycle slip has occurred according to the fitting residuals and RMS, which is not used at present.

After running the above process, the short piece of data are checked again with the adjacent available epoch with larger intervals. Finally, the detection results and data rejection results are written to the log file, including the information related to the deleted observations and the newly added ambiguity.

lsq

lsq is based on the generalized least squares principle of the parameter elimination-recovery method for parameter estimation. As shown in the figure below, the data processing process can be divided into three parts: (1) initialization, obtaining configuration information, variable assignment, and the number of statistical parameters; (2) construction of the function model, constructing the parity mathematical model and filling the matrix by epoch, eliminating process parameters (such as receiver clock error) and state parameters (ambiguity parameters) in this process and compressing the normal equation matrix as needed; (3) adjustment, adjustment for the non-eliminated parameters, recovering the eliminated parameters and calculating the residuals, and inputting the results to different result files

image13

Figure 4 lsq algorithm flow. (The dashed box part is executed only when the ambiguity parameters is eliminated)

  1. Initialization

    First, get the configuration information required for lsq from the configuration file, such as satellite list, a priori constraints and process noise. S/F mode reads the initial coordinates of the station from the pos_ file or “sit.xyz” file; reads the data in the relevant files, such as the PCO/PCV of the receiver and satellite, and the coordinates of the antenna reference point for the station.

    Determine the information about the parameters and normal equation based on the positioning mode and relevant configuration information. The parameters are divided into three categories: constant parameters, process parameters and state parameters; the station coordinates are estimated as constant parameters in S/F mode, the process parameters include station position parameters in K mode, receiver clock error, zenith tropospheric delay parameters and horizontal tropospheric gradient parameters, and the state parameters are ambiguity parameters. If the integer ambiguity resolution method is rounding method, the ambiguity parameters need to be eliminated subsequently; otherwise, if the integer ambiguity resolution method is a search algorithm such as LAMBDA, the ambiguity parameters need to be retained before adjustment to obtain the variance-covariance matrix for integer ambiguity resolution. Thus, different matrix dimensions are assigned according to different integer ambiguity resolution strategies.

    Initialize the parameter vector with the normal equation matrix. Based on the stochastic process, the state transfer matrix of STO or PWC (Piece-Wise Constant) is assigned to the unit array and the white noise state transfer matrix to the zero matrix. Based on the generalized least squares principle, the corresponding diagonal elements of the parameters in the normal equation are assigned a priori weights, i.e.

$$ \begin{cases} N_{b,b}=diag([P_x\ P_yP_z···0]) \tag{4} \end{cases} $$

where the ambiguity parameter part will be filled later.

  1. Constructing the adjustment functional model

    Based on the initialized least squares estimator, construct the function model epoch by epoch and fill the corresponding matrix, and some parameters are eliminated in the estimator in this process. First, read the observation data of the current epoch and the corresponding OSB; read the log file with the deleted satellite and new ambiguity information of the corresponding epoch; if the positioning mode is K mode, read the initial coordinates of the current epoch in the kin file; if there is a receiver clock jump in the current epoch, read the receiver clock jump file generated in tedit and correct it in the observation value of the current epoch. Same as constant parameters and process parameters, update the information related to the ambiguity parameters of the current epoch to the estimator.

    The functional model corresponding to the original observation equation is established based on the satellite constellation selected in the configuration file and existing in the observation file, including the design matrix, OMC (Observed minus calculated) and the weights corresponding to the observations. In this process, the correction of each systematic error is carried out, and the calculated receiver clock error is used as the initial value of the receiver clock error parameter. The functional model of the ionosphere-free observation equation is constructed based on the functional model of the original observation equation, and the M-W combination is composed to calculate the initial value of the wide lane ambiguity, and then calculate the initial value of the ionosphere-free combination ambiguity. Calculate the elements of the normal equation matrix and fill them into the upper triangular matrix, and calculate the OMC weighted sum of squares in the functional model of the ionosphere-free combination for calculating the residual sum of squares.

    The state equation in the PPP is extended with the virtual observation equation as in equation $(5)$, where the state vector is the process parameter.

$$ \begin{cases} V=X_k-\mathrm{\Phi}X_{k-1}+\omega_k,P_w) \tag{5} \end{cases} $$

The process parameters of the previous epoch are eliminated, and the information related to the process parameters of the current epoch is added to the estimator according to the state equation. If the integer ambiguity resolution method is the rounding method, the ambiguity parameter needs to be eliminated. If the lsq is performed after arsig, the SDBS ambiguity “cst_” file generated in arsig is read and appended to the estimator as a strong constraint.

  1. Adjustment

    At the end of the epoch cycle, the remaining integer ambiguity constraint is attached, the process parameters are eliminated, and the remaining ambiguity parameters are also eliminated if the integer ambiguity resolution method is the rounding method, and the normal equation matrix is compressed. The non-eliminated parameters are solved, and the variance of unit weight and the posterior variance of the non-eliminated parameters are calculated.

    Recover the eliminated parameters and calculate the residuals, and output the different results to the corresponding result files. S mode writes the position estimation and other information to the “pos_” file. Count the solvable ambiguity and write it to “amb_” file. if the integer ambiguity resolution method is LAMBDA method, write the information such as covariance matrix, residual sum of squares and degrees of freedom to “neq_” file for LAMBDA method to fix the narrow lane ambiguity.

redig

redig performs residual editing based on the “log_” file generated in tedit and the “res_” file generated in lsq, deletes the gross error and detects the residual small cycle slips, and its data processing flow is shown in Figure 5.

image14

Figure 5 redig algorithm flow. The black dashed line indicates that the process is executed only once

The initialization section first obtains the configuration information required by redig, such as the number of ephemerides, the satellites included and the command line input parameters, etc. Next, the phase residuals in the “res_” file are read, and the status of existing observations is identified (whether the ambiguity has been established, the observations have been deleted, etc.). The RMS of each satellite phase residual is calculated separately and written to the header of the “stt_” file, as well as output to the screen, and the time series of each satellite phase residual is written to the “stt_” file.

After the initialization is finished, the residuals are edited satellite by satellite, firstly checking whether there are short pieces of data and removing them. Secondly, check whether there is residual cycle slip or gross error in the residual time series: (1) calculate the residual difference between adjacent available epochs and its mean and standard deviation, and also calculate the mean and standard deviation after removing the absolute value of the maximum residual difference to detect the jump based on the above calculated value of chi-square test. Repeat the above process until it passes the chi-square test and the maximum residual difference value does not exceed the threshold; (2) if the adjacent available epochs both detect the cycle slips, distinguish the cycle slips from the gross error according to the residual difference value; (3) deal with special cases, such as the last available epoch detects the cycle slip and delete it. Finally, remove the short piece of data after residual editing and update the “log_” file according to the residual editing result.

arsig

arsig uses SDBS to eliminate the hardware delay at the receiver and resolution integer ambiguity for the wide lane and narrow lane, where the wide lane ambiguity is fixed by rounding, and the narrow lane ambiguity can be fixed by the same rounding method when the data processing period is long, otherwise it should be fixed by the LAMBDA method. Its data processing flow is shown in Figure 6.

image15

Figure 6. arsig algorithm flow. ZD indicates the zero difference

As can be seen from the figure, the key in arsig is the narrow-lane integer ambiguity resolution method, and all subsequent integer ambiguity resolution methods described refer to the narrow-aisle fuzziness fixation method. Firstly, initialization is performed and the parameters related to integer ambiguity resolution in the configuration file are read. If the integer ambiguity resolution method is rounding, the “amb_” file is read; if it is the LAMBDA method, the “neq_” file generated in lsq is read. Based on the available information, all satellite pairs are defined for SDBS, and the real values of SDBS ambiguity and its standard deviation are calculated. And calculate the standard deviation of iSDBS narrow lane ambiguity: LAMBDA method is calculated from the covariance matrix and unit weight variance, and the rounding method is empirical.

If the integer ambiguity resolution method is LAMBDA method, firstly the ZD covariance matrix is mapped to the corresponding SDBS covariance matrix. Secondly, narrow lane integer ambiguity resolution is performed according to LAMBDA algorithm, and acceptance test and ratio test are performed at the end of LAMBDA algorithm. If the test is failed, the number of candidate ambiguity is reduced and the ambiguity is fixed again until it can pass the test as well as meet a certain number of ambiguity. After the successful resolution, the information related to the narrow lane ambiguity is updated and the independent SDBS wide lane and narrow lane integer ambiguity are output to the “cst_” file.

Installation and validation - General operation steps - Contact Us - Website

Getting Started

Appendix

Information

Communication

Clone this wiki locally