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.
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.