## 1 Introduction  
The Earth's gravitational field, a fundamental physical field, plays a vital role in numerous applications across various fields, rendering its determination imperative within the geodetic community. The gravitational potential (GP) field, a consequence of the Earth's density distribution, can be determined inside and outside the Earth, essentially across the entire spatial domain, by applying the Newtonian integral formula. As a result, the gravitational field throughout the space is established by applying the gradient operator to the GP field. However, the Earth's density distribution, as per the preliminary reference Earth model (PREM) [1], has been inadequately determined, making the gravitational field based on this distribution insufficient for general application requirements. Fortunately, the Earth's external gravitational field can be established by considering certain gravitational-related distributions, such as gravitation or potential distributions, over the Earth's surface (boundary) [2]. Determining a gravitational field (or GP field) based on the given distribution values on the boundary is termed the geodetic boundary value problem (GBVP). The conventional approach to addressing this boundary value problem involves using spherical harmonic analysis, which successfully represents the external gravitational field through comprehensive coverage of gravitation (or GP) measurements over the Earth's surface or a surface defined by the orbits of a satellite. Over the past few decades, satellite-based gravity measurement techniques have emerged to surmount challenges related to practical measurements on the Earth's surface, particularly in mountainous and oceanic regions. These techniques have proven instrumental in gravitational field determination. [PERSON] [3] introduced a method for formulating a gravity model by calculating the spherical harmonic coefficients based on the observations of orbit perturbations from artificial satellites. This seminal work sparked substantial efforts and led to significant advancements in delineating the Earth's gravitational field using diverse datasets, including satellite radar altimetry data, satellite tracking data, terrestrial gravity data, and combined data sources. The launch of dedicated satellite gravity mapping missions, such as the CHAMP mission in 2000 [4], the GRACE twin satellite mission in 2002 [5], and the GOCE mission in 2009 [6], heightened interest in the recovery of satellite gravitational fields. Furthermore, various techniques have been proposed, including satellite accelerations [7], orbital perturbations [8], harmonic analysis [9], and energy-integral methods [10; 11; 12]. These satellite-based gravity measurement approaches are particularly beneficial because they provide comprehensive coverage of the Earth. Recent advances in time and frequency science have facilitated the development of optical-atomic clocks (OACs), achieving stability and accuracy superior to the \(1\times 10^{-18}\) level over several hours in laboratory settings [13; 14; 15]. Notably, ultra-high-stability onboard and portable satellite clocks are anticipated to become available shortly [16; 17]. This development offers the prospect of determining the GP differences between a satellite and a ground station. Precise atomic-clock-related frequency signal links will be utilized, rooted in the principles of general relativity [18]. Consider a scenario in which a satellite transmits frequency signals received by two different ground stations. The GP difference between these stations can be ascertained by observing the frequency shift [19]. However, the precise extraction of the frequency shift attributed to the GP difference between a ground station and a satellite presents challenges. The primary reason is the contamination of observations by the Doppler effect and influences from the ionosphere and troposphere.