One of the questions raised in The Doppler Shift of Satellite Radio Beacons was whether we could plot the actual orbit of a satellite using accurately timed Doppler shift data from multiple ground stations. I have outlined a method on this page which might give some direction to thinking about this idea.
I do not have actual Doppler shifted data from multiple stations so I used the PREDICT software to generate pass data at three arbitrary locations in the eastern part of the United States. The data returned from PREDICT gave me not only range data that I would otherwise determine using the Doppler shift analysis but the satellite location information to check my position estimates against as well. Using only the station locations and the ranges over a series of discrete times I calculated the satellite position at each time step using a method of trilateration. These positions are directly compared with those returned from PREDICT.
For this exercise I have used one of the precursor flights of a new generation of cubesats being developed in Europe and elsewhere, QB50P1.
The method I develop relies on the method of trilateration as described on Wikipedia. In this method we solve the system of equations describing three spheres. While there are several cases which do not yield real solutions, the solution of interest will result in two common points. This is explained in detail in the Wikipedia article.
To set up our problem I have chosen three ground stations at arbitrary locations in the eastern United States. The northwesternmost station is chosen as P0, the master station. The satellite position will be transformed back to P0 coordinates for comparison to the PREDICT data.
The trilateration method I used is simplified when the centers of the three spheres all lie in the same plane and that two of them lie on the x-axis. This will be accomplished with a series of coordinate system rotations. These are straightforward so without going into the details, I will just describe the flow of the calculation here.
- At ground station P0 Horizon Coordinate system, calculate the distance, bearing, and chordal distance to each of the other two ground stations. In the P0 Horizon Coordinate system the y-axis is North, the x-axis is East, and the z-axis points to the zenith.
- From those values, calculate the x, y, and z coordinates of each station
- Now rotate the P0 Horizon Coordinate system around the z-axis to align the x-axis with the bearing to the P1 ground station.
- Rotate the new coordinate system centered at P0 around the y-axis to place P1 on the x-axis.
- Finally rotate this coordinate system around the x-axis to place ground station P2 in the z=0 plane. This is now the preferred coordinate system for the trilateration calculation.
- For each time/satellite range data point determine the trilateration solution in this coordinate system.
- Now transform the satellite coordinates back to the original P0 Horizon Coordinates by reversing the previous coordinate system transformations.
- We now have the satellite positions [x,y,z] in the Horizon System at P0. This allows us to calculate the azimuth and elevation of the satellite to compare directly to the values presented by PREDICT. I also make estimates of the latitude, longitude of the ground track and altitude.
I coded these steps using Python and made the calculations for four different passes of QB50P1 within the range of the three arbitrary ground stations. Two passes were on the ascending portion of the orbit and two on the descending part. I show some of the representative results below.
Run 1 was a descending pass. The PREDICT which formed the artificial data set for ground station P0 is shown here:
PREDICT Generated Data Set for Ground Station P0, Run 1
Similar data for the same pass was obtained for the other two ground stations. The Time/Range data for all three ground stations were then processed through my trilateration script. The output is shown below:
Results of Trilateration of Ranges at Ground Station P0, Run 1
Ground Station Latitude and Longitude Station Latitude Longitude ------------------------------------- Station P0 40.8771 -81.3994 Station P1 39.1802 -76.6631 Station P2 37.3905 -79.1516 Input Data - Time and Range Data Time P0 Range P1 Range P2 Range --------------------------------------------- 04:29:01 2201.0 2321.0 2543.0 04:30:01 1791.0 1917.0 2131.0 04:31:00 1393.0 1526.0 1725.0 04:32:00 1023.0 1165.0 1335.0 04:33:01 728.0 875.0 982.0 04:34:01 628.0 745.0 726.0 04:35:00 800.0 853.0 687.0 04:36:00 1125.0 1133.0 894.0 04:37:01 1507.0 1489.0 1227.0 04:38:01 1910.0 1878.0 1610.0 04:39:00 2322.0 2281.0 2012.0 04:40:00 2738.0 2692.0 2424.0 Estimated Satellite Location parameters at Ground Station P0 Time range elevation azimuth latitude longitude altitude ----------------------------------------------------------------------------------------- 04:29:01 2201.0 9.2 14 58 -74 694 04:30:01 1791.0 13.8 14 55 -75 645 04:31:00 1393.0 21.7 15 51 -77 637 04:32:00 1023.0 34.3 16 47 -79 628 04:33:01 728.0 57.4 21 44 -80 624 04:34:01 628.0 82.2 155 40 -81 623 04:35:00 800.0 49.0 -173 37 -82 624 04:36:00 1125.0 29.6 -170 33 -83 625 04:37:01 1507.0 18.4 -169 29 -84 623 04:38:01 1910.0 11.4 -168 26 -85 633 04:39:00 2322.0 6.2 -168 22 -86 642 04:40:00 2738.0 3.1 -167 19 -87 697
Now the satellite positions determined from the range/time pairs can be plotted on a map along with the locations determined from the original PREDICT output for ground station P0. Runs 1 and 2 are the descending portion of the orbit and Runs 3 and 4 are the ascending portion. Here are the ground tracks for each of the four runs:
A few observations
- I have yet to determine an actual orbit for the satellite. As I develop an estimate for that I will present that on a new page.
- The crookedness of some of the ground tracks shown above, notably in Run 4, is most likely caused by PREDICT rounding the latitude and longitudes to the nearest degree. For comparison, I did so as well.
I ran Run #4 over again as Run 4hp with a single decimal point precision on my Range Only estimate. The ground track shown below shows a physically more realistic and presumably more accurate track:
- Interestingly, one degree along the earth’s surface is about 110 km. The range discrepancies that I calculated comparing PREDICT range data with Doppler shift range data on three different passes of the now decommissioned HAMSAT satellite were well under that value.
- The target orbital altitude for QB50P1 was to be 624 km. My altitude estimates were very close to this value at the times of closest approach but showed significant deviation at the beginning and end of the pass. I need to look at that a little closer.
- Remember that the information on this page is based on studies done on artificially created and internally consistent data sets using the PREDICT software. Having said that, I think the results are quite remarkable and perhaps even useful under certain circumstances. I am sure there are many opportunities to improve the process, methods, etc.
This looks like a great introduction into Doppler-based passive satellite ranging, thanks! Do you have any thoughts on factors limiting total precision? The ones that come to my mind include atmospheric refraction, signal sampling frequency, sampling oscillator jitter, antenna position precision, time precision at the ground station… Just as a thought experiment, how precise could this position determination be and which factors would be the ones limiting it?
Thanks for taking a look at the page. I have not really looked into the limits of precision for this method and in fact have never actually implemented the physical network to take the data. I got sidetracked on the hardware side trying to build a software defined radio and then sidetracked again by other things so this project has languished. I am sure that you could get estimates of the effects by passing known errors through the calculations but I have not done any of that.
OK, I’ll see what I can dig out. Thank you for a prompt answer!
Hi there! I know it’s been a while since you created this study, but I am doing an orbit determination project and thought this might be a great study to use. I was trying to plot the orbit of EO-79 to match the orbit in your plot for Run 1, but not getting much luck there (I found a TLE online from 2014 and propagated from there to see if I could match it to the starting point of your data). Are the times in UTC or just UT1? It would also be wonderful if you could potentially provide the TLE or elements of the orbit for the time of the orbit above. Thank you very much!
Hello. I am not going to be of much help to you on this as it has been so long since I worked on it. I found the TLE for qb50p1 downloaded Feb 22, 2015, a little later than the ones I used:
1 40025U 14033R 15053.67956488 .00002497 00000-0 28765-3 0 9992
2 40025 97.9649 313.7541 0012330 194.4611 165.6270 14.86584863 36780
The times are UTC. Sorry to say that’s about all I can say…Thanks for your interest.