## 1 Introduction  
The modeling of Earth's gravity field is an essential task in the physical geodesy, and the constructed gravity information can be used for the fields of solid geophysics, oceanography, geodesy, and glaciology [1]. Beginning in the year 2000, several dedicated satellite gravimetry missions, such as the CHAllenging Minisatellite Payload (CHAMP), Gravity Recovery And Climate Experiment (GRACE), and Gravity field and steady-state Ocean Circulation Explorer (GOCE), have significantly improved the accuracy of the static gravity field and its temporal variations by several orders of magnitude [2; 3; 4]. Although these missions were carried out with different measurement principles, they had a common feature: all were equipped with global positioning system (GPS) receivers to realize the concept known as satellite-to-satellite tracking in high-low mode (SST-hl). SST-hl is a highly sophisticated technique to map the long-wavelength spectrum of the gravity field [5; 6; 7]. There are four methods for gravity field recovery that use SST-hl observations: celestial mechanics, energy balance, short-arc, and acceleration methods [8]. According to the types of differentiated acceleration, two acceleration methods exist: the average acceleration method and point-wise acceleration method [9]. [PERSON] et al. [8] and [PERSON] et al. [10] compared the different methods for modeling gravity field models with the GOCE SST-hl data and pointed out that, with the exception of the energy balance method, the other three methods had a comparable performance. The celestial mechanics approach is confronted with a large computational effort due to the integration of the variation equations and a nonlinear system of equations; in the case of the energy balance method, space gravity spectroscopy information is exclusively scalar and not available in all three directions [11; 12; 13]. The acceleration method is based on [PERSON]'s equation of motion, which balances the gravitational vector with satellite accelerations in the inertial reference frame [14; 15]. The principle of acceleration approach is simple and has been successfully applied to CHAMP, GRACE, and GOCE SST-hl data analysis [14; 15; 16; 17; 18; 19]. The key to recovering the gravity model by the acceleration method is to estimate the satellite accelerations with a high degree of precision. Several approaches have been used to estimate the GPS-based acceleration including orbit differentiation, Doppler differentiation, and carrier phase differentiation approaches [20]. The most commonly used method is performed by double time-differentiating successive trajectories of the moving satellite, also known as the orbit differentiation algorithm. A major drawback to this method is that the accuracy of the differentiated accelerations is strongly dependent on the position precision; an increase or decrease in the number of visible satellites can lead to discontinuities [21]. Another method is the Doppler measurement, which can be used afterward to obtain the satellite's velocity before the satellite accelerations are computed by the first-order derivative of the satellite's velocity. One shortcoming is that the raw Doppler observables may be heavily contaminated by measurement noises [22]. An alternative method is to numerically differentiate the GPS carrier phase observations then obtain both the range rate and range acceleration, which is called the carrier phase differentiation method. This method has been successfully applied in airborne gravimetry [23; 24] and satellite gravimetry [25]. The most prominent advantage of this approach is the avoidance of the cycle ambiguity problem and the elimination of systematic errors (e.g., hardware delays from the receiver and GPS satellites) that are common to adjacent epochs [23; 25]. [PERSON] et al. [25] adopted the average acceleration approach to recover the gravity field model from satellite accelerations derived from onboard GPS phase measurements.