High-precision GPS positioning using MADOCALIB

Introduction

By receiving radio waves from positioning satellites such as the US GPS, you can know your current coordinates. In addition to GPS, other global navigation satellite systems (GNSS) include Russia’s GLONASS, Europe’s Galileo, Japan’s quasi-zenith satellite system (QZSS), China’s BeiDou, and India’s Navigation with Indian Constellation (NavIC).

This positioning signal contains a ranging signal to measure the pseudo-distance from the satellite to the receiver, and a navigation message signal to transmit the satellite position, time, and radio path delay. The receiver estimates its own coordinates based on these two signals from the satellite.

However, the satellite coordinates, satellite time, and ionospheric delay on the radio wave path in the navigation message may differ from reality due to such as the modeling process. Such differences from reality lead to positioning errors. The process of reducing such errors and improving positioning accuracy by estimating these correction amounts using a satellite radio wave observation network on the ground and transmitting them to the receiver is called augmentation, and the high-precision positioning achieved in this way is called PPP (precise point positioning).

Japan’s Quasi-Zenith Satellite, Michibiki, broadcasts three types of augmentation signals. One of these is MADOCA-PPP (multi-GNSS advanced orbit and clock augmentation - precise point positioning). And the reference program for using MADOCA-PPP is MADOCALIB (MADOCA library).

An upgraded version of MADOCALIB was released on July 1, 2024, enabling it to process augmentation signals broadcast by future Michibiki satellites. Using this new MADOCALIB, we will try to experience future high-precision positioning.

Augmentation information broadcast from positioning satellites includes Michibiki’s CLAS (centimeter level augmentation service) and MADOCA-PPP, Galileo’s HAS (high accuracy service), and BeiDou’s PPP-B2b. Since CLAS and MADOCA-PPP include a carrier phase bias in the augmentation information, it is expected that the AR (amibiguity resolution) method will be able to obtain a fix solution with higher accuracy than the PPP solution.

Downloading and preparing MADOCALIB

Here, we will use a Windows PC to process sample data in MADOCALIB.

The Cabinet Office of Government of Japan has released MADOCALIB as open source software on GitHub since June 5, 2024. On July 1, 2024, MADOCALIB was updated to version 1.2. You can download MADOCALIB by clicking the link below and clicking Code → Download ZIP (or by using git client software).

For example, create a folder called madocalib on your desktop and collect the necessary files there. From the downloaded MADOCALIB, copy the following to this folder:

• rnx2rtkp.exe in the bin folder,
• sample.conf in the app\consapp\rnx2rtkp\gcc_mingw folder (you will edit the contents later), and
• 2024162all.200.l6 in the sample_data folder, etc. (files with bat in them are not necessary).

MADOCALIB_Detailed manual_ver001.pdf in the doc folder is the MADOCALIB manual, and MADOCALIB_manual_ver002.pdf is the manual for the ionospheric delay estimation function newly added in v.1.2.

Therefore, we will modify this sample.conf following the procedure written in the latter manual.

In addition, for data conversion and display, we will also use rtkplot.exe, rtkconv.exe, and convbin.exe from the Windows binary file of RTKLIB v.2.4.3b34. Copy these to this folder as well.

Running the included sample

Then, open a command prompt, navigate to this folder, and perform positioning. Using rnx2rtkp.exe, confirm that positioning can be performed using only the ranging data TSK200JPN_S_20241620000_01D_30S_MO.rnx and the navigation data TSK200JPN_S_20241620000_01D_MN.rnx. This method is called single positioning.

rnx2rtkp.exe -p 0 -sys G,R,E,J -f 2 TSK200JPN_S_20241620000_01D_30S_MO.rnx TSK200JPN_S_20241620000_01D_MN.rnx

This filename is quite long, but the command prompt will complete the rest of the filename if you type part of it and press Tab. If a different filename is completed, you can press Tab again until the correct filename is found.

The options specified here are as follows. These are options to perform single positioning without using sample.conf.

• -p 0 specifies independent positioning.
• -sys G,R,E,J specifies using GPS (G for GPS), GLONASS (R for Russia), Galileo (E for Europe), and QZSS (J for Japan) for positioning. We specify this option to use the same satellite system as the file specified in sample.conf.
• -f 2 specifies using L1 and L2 frequency bands. This is the default behavior, so there is no need to specify it.

This ranging data was collected every 30 seconds on June 10, 2024 (162 days from January 1), and covers 24 hours of data. From the file name, it can be assumed that the data is from the Geospatial Information Authority of Japan’s GNSS observation point in Tsukuba City, Ibaraki Prefecture (coordinates 36.1056N, 140.0871E). If you create an account on Mirai, you can also download real-time data and archive data for other locations and dates and times. You can create a Mirai account for free.

These files are written in RINEX (receiver independent exchange format), a plain text format, so you can view the contents with a text editor.

If you see a bunch of text and what appear to be coordinates like 36.105576920 or 140.087107164 then everything is working properly. Add the -o single.pos option to save the results to the file single.pos.

Next, we add the satellite orbit and satellite clock augmentation information 2024162all.204.l6 to determine coordinates. We also add the ionospheric delay augmentation information 2024162all.200.l6, which is newly available from v.1.2 of MADOCALIB.

rnx2rtkp.exe -k sample.conf TSK200JPN_S_20241620000_01D_30S_MO.rnx TSK200JPN_S_20241620000_01D_MN.rnx 2024162all.204.l6 -mdciono 2024162all.200.l6 -o mdcar200.pos

The options specified here have the following meanings:

• -k sample.conf reads options from sample.conf instead of defaults,
• -mdciono 2024162all.200.l6 uses 2024162all.200.l6 as augmentation for ionosphere delays, and
• -o mdcar200.pos outputs results to file mdcar200.pos instead of screen.

In this sample, 2024162all.200.l6 and 2024162all.201.l6 are prepared as augmentation information for ionospheric delay. It seems that these two pieces of information cannot be used at the same time, so we will use the former first.

Let’s plot these results and compare them. Start rtkplot.exe and select single.pos from File → Open solution-1... menu and select mdcar200.pos from Open solution-2... menu.

Here, the lines are spaced 50 centimeters apart, the red plots are standalone positioning results, and the light blue plots are positioning results using the new MADOCA-PPP. Here, to enable the display of statistical information in the upper right corner, we have turned on Show Statistics in the Edit → Options... menu.

While the single positioning plots vary widely over the course of 24 hours, the MADOCA-PPP results plots are more concentrated. According to the statistics, the coordinates are 36.105591478N, 140.087128046E. Although the exact coordinates of the observation points are unknown, the MADOCA-PPP augmentation information has reduced the variation.

Next, to see the changes over time, change the menu in the upper left corner of RTKPLOT from Gnd Trk (ground track) to Position.

This mode change will display three graphs against time. From the top to bottom, they show coordinate changes in easting, northing, and up-down, in meters. The horizontal axis shows two-hour intervals. With single positioning, the position changes greatly over time. But with MADOCA-PPP positioning, there is almost no change except in the initial state.

Now, click on button 1 to display only the MADOCA-PPP results. For about 15 minutes after the start of positioning, the result will be light blue (PPP solution), but after that it will turn dark blue (Fix solution). The Fix solution is a self-evaluation of the state in which the system continues to resolve discrepancies in the ranging augmentation using AR, and achieves a higher accuracy of positioning than PPP, which uses corrections based on augmentation information.

There are two files for ionospheric delay augmentation data, and so far we have been using the one with PRN (pseudo random noise) number 200. Similarly, perform positioning calculations with PRN201, and plot these two position changes on the same graph. Rotate the mouse wheel on the time axis to expand the time axis. The horizontal axis ruled interval is set to 5 minutes. In the PPP positioning state 10 minutes after the start of positioning, the purple plots are for PRN200, and the green plots are for PRN201.

The time course of the positioning coordinates differed between when PRN200 was used and when PRN201 was used. However, after about 10 minutes, both showed almost identical coordinates.

MADOCALIB high-precision positioning using our own observation data

Next, we will perform MADOCA-PPP positioning using my GNSS observation station data and the L6 signal archive published on the Cabinet Office’s Michibiki website.

Here, we use the raw data 20240703a.gps recorded with the NovAtel OEM729 receiver, one of the two receivers, for one hour from 2023-07-03 00:00:00 UTC (UTC is coordinated universal time). This raw data includes both navigation data and ranging data. In fact, I also tried using raw data from the u-blox ZED-F9P receiver, but only about eight satellites were available, so I gave up using this.

We will also need the satellite orbit and clock data and ionospheric delay data from MADOCA-PPP for the same time period. These can be downloaded from the archives on the Michibiki website. We must agree to the terms of use to download the data. Then, we search for the date and download the necessary files.

Of these archived data, those related to satellite orbits and clocks are currently being broadcast by Michibiki, but those related to ionospheric delays have not yet been broadcast by Michibiki. The following table summarizes the PRN number and content correspondence table on page 6 of the signal specification IS-QZSS-MDC-002, excluding CLAS information and adding the Michibiki satellite number.

Since Michibiki 5, 6, and 7 have not yet been launched, the possible satellite numbers are written in parentheses. The satellite orbit and satellite clock augmentation information content is currently the same for all satellites, so download one of the PRN204 to PRN211. Here, we download 2024185A.209.l6, which is for PRN209. Here, 2024 stands for 2024, 185 stands for 2023-07-03, which is the 185th day from January 1st, A stands for the start time of the one-hour data from 00:00:00 UTC, and 209 stands for PRN209.

On the other hand, it seems that different information about ionospheric delay augmentation is distributed depending on the PRN. Here, we will download both 2024185A.200.l6 and 2024185A.201.l6. By searching the ionospheric delay information by date, I found that the ionospheric delay data will be available from 00:00 UTC on June 28, 2024.

We put these downloaded files into the madocalib folder created earlier.

Next, we use rtkconv.exe to extract the navigation message file and the ranging file from 20240703a.gps. Click the Convert button to create 20240703a.nav, which is the navigation data file, and 20240703a.obs, which is the ranging data file.

We can also run this process on the command prompt using convbin.exe. Here, we exclude GLONASS, SBAS, BeiDou, and NavIC data with the -y option.

convbin.exe -y R,S,C,I 20240703a.gps

These two files are used for single positioning. The mask angle is the default value of 10 degrees, and the frequency bands are also the default values ​​of L1 and L2. The satellite systems are limited to three systems of GPS, Galileo, and Michibiki.

rnx2rtkp -p 0 -t -sys G,E,J 20240703a.nav 20240703a.obs -o 20240703a-single.pos

This output file 20240703a-single.pos was plotted using RTKPLOT.

There are multiple discontinuities in the coordinates over the course of an hour. This was not good observation data. Next, we will display the changes in the coordinates over time.

Although there are discontinuities in the coordinates, we use the MADOCA-PPP augmentation information for this data. First, we calculate and plot the position using only the augmentation information of the PRN209 satellite orbit and clock.

rnx2rtkp.exe -k sample.conf 20240703a.nav 20240703a.obs 2024185A.209.l6 -o 20240703a-mdc.pos

Here, the spacing between the lines is 20 centimeters. The position changes significantly for about 15 minutes after the start of positioning, but then converges to a constant coordinate. Next, the ionospheric delay information of PRN200 was added and the positioning calculation was performed. Here, -x 2 is the trace option for debugging, and a debug file called 20240703a-mdc200.pos.trace is created.

rnx2rtkp.exe -k sample.conf 20240703a.nav 20240703a.obs 2024185A.209.l6 -mdciono 2024185A.200.l6 -o 20240703a-mdc200.pos -x 2

From this result, the presence or absence of ionospheric delay correction information made little difference to the plot. The manual MADOCALIB_manual_ver002.pdf mentioned above describes how to read the debug file. When searching for a line containing miono_sel_area in this debug file, the only line that was found was 2 miono_sel_area: lat=34.44,lon=132.41. This indicates that ionospheric delay information was not used.

Let’s try using the ionospheric delay information from PRN201 in the same way. In the debug file for PRN201, in addition to 2 miono_sel_area: lat=34.44,lon=132.41, there was also a line that read 2 miono_sel_area: closest RegionID=5,AreaNo=5,dist=267.559. It seems that the ionospheric delay information was used.

So, let’s compare the Gnd Trk of PRN200 and PRN201. The positioning result of PRN201 had smaller coordinate change.

Previous graph was wrongly plotted and is corrected on 2024-07-09.

Next, we compare the time-dependent changes in coordinates for both. Both converge to almost the same coordinates in about 15 minutes, but in this case, the coordinate changes using PRN201 were smaller.

We found that the ionospheric delay information newly added to MADOCA-PPP has the effect of providing positioning results with little time change in coordinates from the start, and that the results may vary depending on the ionospheric augmentation information provided.

In fact, during the course of this experiment, we noticed that the time synchronization (NTP) settings on the data collection server at the GNSS observatory at our workplace were incomplete, resulting in a time delay of about one minute. I have already corrected this.

MADOCA-PPP augmentation information distribution via the Internet

The Cabinet Office of Government of Japan has started distributing MADOCA-PPP on the Internet. To receive the Internet distribution, application is required. Distribution is done in NTRIP format using the Mirai system. When usage is permitted, the NTRIP caster address, ID, password, and mount point name are notified. Nine mount points are prepared for PRNs, and the number of simultaneous usage permitted sessions is nine. In other words, we can use the augmentation information of all PRNs at the same time. The expiration date of usage is the last day of March each year, and application is required for continued use.

With internet distribution, you can get information faster because there are no transmission delays caused by such as a satellite repeater (transponder). The screen on the left below shows the augmentation message of PRN205 distributed via the internet, and the screen on the right shows the augmentation message of PRN205 receiving the Michibiki L6E signal with an Allystar receiver. The information was available about 9 seconds faster via internet distribution.

For example, QZNMA G06 (CNAV) on the 7th line of the left screen is on the 16th line of the right screen. Since MADOCA-PPP transmits one data part per second, this corresponds to a delay of 9 seconds for the L6E signal.

Conclusion

MADOCALIB has been updated to v.1.2, and in addition to the high-precision positioning using AR, it can now handle ionospheric delay augmentation information, which allows coordinates to be determined in a short time. This ionospheric delay information will be broadcast by the Michibiki 5, 6, and 7 satellites, which are scheduled to be launched, and can be obtained from the official archive.

The PPP (precise point positioning) requires about 15 minutes for initial convergence of the coordinates, during which the coordinate output is unstable. In MADOCA-PPP, the newly available ionospheric delay information can reduce this coordinate change, so it can be said that coordinates can be obtained more quickly. In the future, I will consider the differences in the content of the two ionospheric delay augmentation data.

With the successful launch of the H3 rocket No. 3 (H3F3) by JAXA, the early start of operation of Michibiki 5, 6, and 7 has become a reality. I am very excited about the future of Michibiki.