Summary
This post explores possibilities for estimating satellite orbits based on Doppler variations of their signals over time. Selected orbits parameters have been estimated using the Python's SciPy optimisation library. Initial experiments have been promising but more remains to be done.
Background
There are many useful applications of the Doppler effect. I found this recent example particularly interesting in the use of Doppler measurements to detect exoplanets.
Some years ago I was trying detect low earth orbit (LEO) satellites in the radio amateur bands by looking for their characteristic Doppler profile. As a satellite approaches a ground station (GS) its radio signal frequency appears higher, while overhead there's no frequency offset, and as it recedes the carrier frequency continues to drop. Several satellites can be seen in the "waterfall" plot below, with time increasing from top to bottom and frequency from left to right. The vertical scale is 60 minutes but the plot superimposes results from a few hours. Although this spectrum region is dedicated to radio amateur satellites, there are often sources of terrestrial interference which obscure the small satellite signals, so part of the task is to automatically remove the nonsatellite signals (a few remnants can still be seen below). In addition the spectral width of LEO signals varies, so a filtering and median tracking approach has been used to only estimate the carrier frequency.

Showing spectral signatures of several LEO satellites in the 70 cm band, over 60 minutes. 
An Orbit Estimation Experiment. While the "blind" detection of LEO satellites is interesting in its own right, this post asks whether some orbital parameters can be estimated just using their Doppler shapes. Doppler shift is proportional to the satellite's velocity in the direction of the GS, so overhead passes show greater Doppler variations than low elevation passes. Initially I assumed it would be possible to estimate the minimum distance from the GS to the satellite and perhaps it's orbital height. Many LEO satellites have almost circular orbits and I've made that assumption below. In addition it seemed likely that collecting the same signal from two different locations might provide information regarding orbit inclination. During the investigations I realised that the Earth's motion must be taken into account, which provides some additional possibilities, even for one GS.
Others have used Doppler measurements to refine existing orbit estimates, which are usually disseminated by the TLE standard. Mark alerted me to the strf software and I changed my data collection format to use their file standard. So far I've used signals collected from two satellites, by myself and Mark. On April 16th we each recorded a pass from LUSAT (LO19) and FalconSat3. Both of these provide relatively strong continuous downlink transmissions, and they use different orbits and signal characteristics. The figure shows a summary of the satellite passes used. These GSs are located near Adelaide and are almost 90 km apart.
The figure below used the strf code and shows one of the FalconSat3 signals. In this case the signal bandwidth is significant and I've written some Python code to implement a simple comb filter to track the carrier frequency.
Given a set of carrier frequency measurements over several minutes, plus an orbit model with unknown parameters, it's a nontrivial problem to estimate the orbital unknowns. Fortunately SciPi provides a versatile
library function to address this type of problem. My model uses the following unknowns: time offset, frequency offset, angle between GS and the satellite's crossing of the equator, orbit height and orbit offset from 90 degrees inclination.
The figure below shows an example of Doppler measurements and fitted Doppler trajectory. The blue points are automatically extracted from a selected portion of the raw strf file which contains spectral estimates from a simple RTL2832 software defined radio, connected to a low gain antenna. The orange curve represents the modelled trajectory after the MSE optimisation. In this case the estimated orbit height of 487 km was quite accurate, plus the time of minimum range when compared to TLE predictions using
gpredict.

FalconSat3 Path Estimate


LUSAT Path Estimate

Also shown is the model's estimated path in 3D of FalconSat3 (blue dots) relative to the GS (yellow X). In both cases a red marker shows the last time position. Note this path illustrates a descending node pass (i.e. moving away from the equator in green), like the figure above. When the satellite direction in the model was assumed to be ascending, the fitted solution was poor. The model didn't make a good prediction of inclination  ~45 degrees compared to 35 actual. In the LUSAT case, the minimum range estimate was accurate to within a few km, but the orbit height was quite poor (1109km compared to about 790 km). In this case the ascending model gave a slightly better fit than the descending node model.
Conclusions
At times the Doppler trajectory seems able to give surprisingly accurate orbit estimates. Note that the sample recordings used strong signals and captured more of the trajectory than a low SNR pass would provide. In the nearpolar LUSAT case, with similarly good data, the results are not so close and seem somewhat dependent on the initial estimates given to the estimator. Perhaps this is not surprising given the nonlinear nature of the model and the reduced affect of GS motion compared to a nonpolar orbit. A simpler model, without inclination, seems to give robust results but doesn't provide as much information. This MSE solver has many parameters options and I've barely scratched the surface in its use.
The original idea of using two GS recordings from the same pass hasn't been explored in any detail. I expect a larger distance separation will be needed to get useful results. A model which jointly fits (say) the frequency difference between the two signals is likely to work better than making separate estimates then combining the results. In addition, at present the signal strength measurements are not used at all! The use of (say) two antennas (whose gain patterns are known), pointed in different directions at one GS, might open many other options for orbit estimation.