Improved thermospheric neutral density models are required to better predict LEO satellite orbits. This research describes a new method to estimate the density using a tomography-based approach, where the orbit states of satellites serve as the measurements. The variational equation for the semimajor axis due to perturbing drag acceleration is used to relate the change in osculating specific mechanical energy of the orbit to the integrated density over the orbit. Using several such measurements from a number of satellites, one can estimate the density scale factor (i.e. a correction to an assumed density model). The problem considered here uses measurements from 100 satellites and solves for the spatially resolved global density scale factor discretized over 324 grid elements spanning 300 to 500 km altitude. This ill-posed problem is solved using Tikhonov regularization, with the 3D gradient chosen as the regularization operator, resulting in a penalty on the spatial smoothness of the estimated density. Preliminary simulation results show that the true time-averaged density can be reconstructed to within approximately 10%, using only simulated ground-based tracking measurements separated over 5 orbital revolutions.