Earth gravitational field interacts closely with a variety of natural phenomena. Therefore, studying and modeling the gravitational field and keeping tracks of its changes over time can make a huge contribution to the related geophysical and environmental researches. A well-proven method which is used in dynamic estimation of time varying state parameters is "Kalman filter". Kalman filter exploits dynamical properties of the states and noisy measurements to estimate the optimal state for every time epoch, in its filtering scheme. Smoothed properties of the involved dynamics, the accuracy of the noisy measurements, pay for noises in measurements and incompleteness of the interchangeable dynamics provide optimal estimates of state for the dynamic system. "Standard Kalman Filter" is performed to acquire the optimal solution when dealing with linear dynamics with white noises in both dynamics and measurements while a non-linear system cannot be dealt with in the same way. Various versions of Kalman filter have been introduced so far, including Extended Kalman Filter (EKF) and Unscented Kalman Filter (UKF). The former uses linearization to derive locally-linear dynamics while the latter is the approximation of the probability distribution of the states using sampling points of the state vector which are commonly referred to as "sigma points". Due to the linearization involved in EKF, the method does not give precise solutions when there are non-linearities of high orders while UKF solution is significantly better in the sense of precision, accuracy and optimality. In this paper, the 9-points Newton NUMERICAL DIFFERENTIATION is applied on the range values calculated from the relative state vector to reconstruct the so-called KBR range-rate observations and EKF is performed along with positions of observations to give an optimal estimate of the state vector at each time epoch. The calculated range-rates from both methods are then compared statistically. Then, UKF is used to fuse position and range-rate observations and to estimate state vector. Finally, the residuals of the reconstructed range-rate observations of three methods are compared statistically. The results show a relatively significant improvement in KBR observations reconstructions for UKF estimated states. To re-evaluate the superiority of UKF from the practicality point of view, the reconstructed range-rates are used to estimate the harmonics coefficients of geopotential model in a field recovery problem using "acceleration method" in a closed-loop simulation for EGM96. The degree-variances and differences in geoid heights are calculated and plotted. The results show a near 5 to 10-times fold improvement in gravitational field recovery using UKF results compared to the other two mentioned methods.