Gravity field modeling using space frequency signal transfer technique between satellites
mmanuscript submitted to
JGR: Solid Earth
Gravity field modeling using space frequency signaltransfer technique between satellites
Ziyu Shen , Wen-Bin Shen , , Xinyu Xu , Shuangxi Zhang School of Resource, Environmental Science and Engineering, Hubei University of Science andTechnology, Xianning, Hubei, China Time and Frequency Geodesy Research Center, Department of Geophysics, School of Geodesy andGeomatics, Wuhan University Key Lab of Surveying Eng. and Remote Sensing, Wuhan University
Key Points: • Low- and high-orbit satellites are connected by frequency links controlled by pre-cise clocks. • Gravitational potential along the orbit of low-orbit satellite is determined. • Optic-clock-based Earth’s gravity field could be established.
Corresponding author: WenBin Shen, [email protected]
Corresponding author: Xinyu Xu, [email protected] –1– a r X i v : . [ phy s i c s . g e o - ph ] A ug anuscript submitted to JGR: Solid Earth
Abstract
Here we provide an alternative approach to determine the Earth’s external gravitationalpotential field based on low-orbit target satellite (TS), geostationary satellites (GS), andmicrowave signal links between them. By emitting and receiving frequency signals con-trolled by precise clocks between TS and GS, we can determine the gravitational poten-tial (GP) at the TS orbit. We set the TS with polar orbits, altitude of around 500 kmabove ground, and three evenly distributed GSs with equatorial orbits, altitudes of around35000 km from the Earth’s center. In this case, at any time the TS can be observed viafrequency signal links by at least one GS. In this way we may determine a potential dis-tribution over the TS-defined sphere (TDS), which is a sphere that best fits the TS’ or-bits. Then, based on the potential distribution over the TDS, an Earth’s external grav-itational field can be determined. Simulation results show that the accuracy of the po-tential filed established based on 30-days observations can achieve decimeter level if op-tical atomic clocks with instability of 1 × − τ − / are available. The formulation pro-posed in this study may enrich the approachs for determining the Earth’s external grav-ity field. The Earth’s gravity field is a fundamental physical field of the Earth. Since grav-ity field has various and significant applications in many fields and branches, its deter-mination is one of the main tasks in geodetic community. If the density distribution ofthe Earth is given, one may determine the gravitational potential (GP) field both insideand outside the Earth (namely in whole space domaion) by Newtonian integral formula,and consequently the gravity field in whole space is determined by applying the gradi-ent operator to the gravity potential (geopotential) field, where geopotential is the sumof GP and the centrifugal force potential generated by the Earth’s rotation. However,since the Earth’s density distribution (e.g. the preliminary reference Earth model, PREM)(Dziewonski & Anderson, 1981) was poorly determined, the gravity field determined basedon density distribution cannot satisfy the general application requirements. Fortunately,one can determine the external gravity field with successive accuracy requirement if somekind of distribution related to gravity (for instance the gravity distribution or gravitypotential distribution) over the Earth’s surface (boundary) is given (Hofmann-Wellenhof& Moritz, 2005). How to determine a gravity field (or equivalently gravity potential field)based on the given distribution values on the boundary (e.g. the Earth’s surface) is re-ferred to as the geodetic boundary value problem (GBVP). A generally accepted approachto solve this boundary problem is to use the spherical harmonic analysis, which can suc-cessfully express the external gravity field once a full coverage of the gravity (or geopo-tential) measurements over the Earth’s surface or a surface (e.g. a surface defined by theorbits of a flying satellite) enclosing the whole solid Earth is provided. To overcome dif-ficulties in practical measurements on the Earth’s surface, especially in mountain areasand ocean areas, there appeared different satellite-based gravity measurement techniques,which have their own advantages especially in the aspect of full coverage over the Earth.Kaula (1966) proposed the method of establishing gravity model by observing theorbit perturbation of artificial satellites, and solve the coefficient of GP. Since the newsatellite gravity missions appeared (e.g., the CHAMP mission (Reigber et al., 2002) launchedon 2000, the GRACE twin satellite mission (Tapley, Bettadpur, Watkins, & Reigber, 2004)launched on 2002, the GOCE mission launched on 2009), scholars have paid extensiveattention on the recovery of satellite gravity field and various methods have been pro-posed, such as orbital perturbation (Hwang, 2001), harmonic analysis (Reubelt, Austen,& Grafarend, 2003), satellite accelerations (Ditmar & Sluijs, 2004) and energy integral(Han, Jekeli, & Shum, 2002; Jekeli, 1999; Visser, Sneeuw, & Gerlach, 2003). These satellite-based gravity measurement techniques have their own advantages especially in the as-pect of full coverage over the Earth. –2–anuscript submitted to
JGR: Solid Earth
For example, the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE)mission, which aims to make detailed measurements of Earth’s gravity field, leading todiscoveries about gravity field determination and ocean circulation investigations (Drinkwa-ter, Floberghagen, Haagmans, Muzi, & Popescu, 2003; Hirt, Kuhn, Featherstone, & G¨ottl,2012). GOCE satellite system flies in a near-polar orbit with an altitude of about 250km above the ground, consisting of an on-board three-axis gravity gradiometer, GNSSreceiver, satellite-to-satellite tracking and relevant equipments (Bock et al., 2011; Hirtet al., 2012). Concerning the global gravity field determination aspect, the GOCE’s mis-sion may map gravity field features with 1 to 2 cm accuracy for geoid undulations andabout 1 mgal for gravity, down to scales of about 100 km, or spherical harmonic degreeabout 200 (Hirt et al., 2012; Pail et al., 2011). Since the GOCE satellite retired in 2013,the current satellite gravity mission on-going is the GRACE Follow-On (Kornfeld et al.,2019), which is a twin satellite system with the height of about 500 km. These satellite-based gravity measurements have greatly improved our understanding of the Earth’s grav-ity field (Flechtner et al., 2016; Kvas & Mayer-G¨urr, 2019; Pail et al., 2019).In recent years, thanks to the quick development of time and frequency science, theoptical-atomic clocks (OACs) with stability and accuracy better than 1 × − levelin several hours have been developed in laboratory environment (Huang, Guan, Zeng,Tang, & Gao, 2019; McGrew et al., 2018; Mehlst¨aubler, Grosche, Lisdat, Schmidt, & Denker,2018; Oelker et al., 2019). Especially, portable and on-board-satellite clocks with ultra-high stability will be available in the near future (Altschul et al., 2014; Hannig et al., 2019;Schiller et al., 2012). This provides potential realization in the near future to determinethe GP differences between a satellite and a ground station using precise atomic-clock-related frequency signal links based upon general relativity theory (GRT) (Einstein, 1915).Suppose a satellite sends frequency signals and two receivers on different ground stationsreceive the signals, then the geopotential difference between the two stations can be de-termined by observing the frequency shift (W. Shen, Ning, Liu, Li, & Chao, 2011). How-ever, how to practically and precisely extract the frequency shift signals caused by thegeopotential difference between a ground station and a satellite is a challenging prob-lem, due to the fact that Doppler effects, ionosphere and troposphere effects contami-nate the observations seriously. In order to overcome these difficulties, recently a moreprecise formulation of the satellite frequency signal transfer (SFST) approach based ontri-frequency combination technique was established (Z. Shen, Shen, & Zhang, 2016, 2017),which aims to determine the GP difference between a satellite and a ground station orbetween two satellites at an accuracy level of several centimeters if high-precise frequencysignal links are established. To precisely compare the frequency signals, the relative sta-bility of clocks should reach about 10 − in several hours (corresponding to about 1 cmin height) for the practical applications of SFST method in geodesy.Based on the SFST technique, in this study we formulate an alternative approachto determine the Earth’s external gravity field, which is completely different from theconventional ones. The basic idea, which was put forward several years ago by our group(Z. Shen & Shen, 2017) is that the GP along a low-orbit target satellite (TS) can be de-termined using frequency signal links between TS and geostationary satellites (GSs). Insection 2 we briefly introduce the relativistic geodesy and SFST technique, by which theGP difference between the TS and GS could be determined. Then we formulate an ap-proach to show how to determine a GP distribution over a TS-altitude defined sphere(TSS), which is bounded by the flights of the TS and defined as the sphere that can bestfit the TS’s orbits. Given the GP distribution over the TSS, a global gravity field (or Earth’sgravity model) could be determined. In sections 3 and 4, we conducted simulation ex-periments, and results show that the proposed approach in this study is prospective. Insection 5 we summarize the main results and discuss relevant potential issues. –3–anuscript submitted to JGR: Solid Earth
The GRT predicts that the frequency (or tick rate) of a clock is related to the geopo-tentials at the place where the clock is located. Specifically, suppose two clocks are lo-cated at different positions P and Q where the geopotential values are W P and W Q re-spectively, and accurate to c − , the frequencies f P and f Q of the two clocks satisfy thefollowing equation (Bjerhammar, 1985; Weinberg, 1972) W P − W Q = f P − f Q f · c + O ( c − ) , (1)where c is the speed of light in vacuum, f = ( f P + f Q ) / O ( c − ) are high order termswhich can be neglected in the case that the two stations are in the vicinity of Earth. Ifthe clock frequencies f P and f Q are precisely measured and compared, the geopotentialdifference W P − W Q between P and Q can be derived. The study of geodesy problems(such as geopotential determination) by the method of clock comparison is regarded asrelativistic geodesy (Flury, 2016; Puetzfeld & L¨ammerzahl, 2019).Currently, there are three kinds of method to compare clocks located at differentplaces: (1) clock transportation (Grotti et al., 2018; Kopeikin et al., 2016), (2) transferfrequency signals via optical fibre links (Z. Shen, Shen, Peng, et al., 2019; Takano et al.,2016; Wu, M¨uller, & L¨ammerzahl, 2019), and (3) transfer frequency signals via satelliteand free-space links (Deschˆenes et al., 2016; Z. Shen et al., 2017). The first two meth-ods are suitable for clocks comparison on ground, while the third method is designed forsatellite clocks comparison. But transferring frequency signals via satellite is much morecomplex than Eq. (1). For example, the satellite is in high-speed motion state which givesrise to Doppler effects; the mediums in space (such as ionosphere and troposphere) willalso cause frequency shifts during a microwave or optical signal’s propagation throughthem. In order to address these problem, Kleppner, Vessot, and Ramsey (1970) proposeda method to transfer microwave frequency signal between a satellite and a ground site;and it is successfully applied for verifying Einstein’s equivalence principle (Vessot & Levine,1979; Vessot et al., 1980). The main idea of the frequency transfer method is that a satel-lite and a ground site are connected by 3 microwave links simultaneously (as depictedin Fig. 1). In this case the first order Doppler effect and most of the medium influencewill be canceled out in the output beat frequency ∆ f ∆ ff e = f (cid:48) s − f s f e − ( f (cid:48)(cid:48) e − f (cid:48) e ) + ( f (cid:48) e − f e )2 f e . (2)where f e and f s are emitted frequency from ground site and satellite respectively; theyare received as f (cid:48) e and f (cid:48) s at satellite and ground site respectively. For each microwavelink, the emitted frequency values are different from the received values. When the fre-quency signal f (cid:48) e is received at satellite, it is transmitted immediately and received atground site as f (cid:48)(cid:48) e . Details can be referred to Vessot and Levine (1979).Kleppner’s method was later improved and introduced to relativistic geodesy forGP determination (Z. Shen et al., 2016, 2017), and was regarded as satellite frequencysignal transmission (SFST) method. According to SFST, the GP difference between asatellite and a ground site is given in the following form (Z. Shen et al., 2017)∆ φ es c ≡ φ s − φ e c = ∆ ff e − v s − v e c − (cid:88) i =1 q ( i ) + Λ f + O ( c − ) , (3)where ∆ φ es is the GP difference between the satellite and the ground station, v s and v e are velocities of satellite and ground site respectively, q ( i ) ( i = 1 , , ,
4) are quantitiesrelated to the positions and velocities of the satellite and ground site, second Newtonian –4–anuscript submitted to
JGR: Solid Earth
Figure 1.
Ground station E emits a frequency signal f e at time t , denoted by uplink (blueline). Satellite S transmits the received signal f (cid:48) e (the downlink denoted by blue line) and emits anew frequency signal f s at time t (the downlink denoted by dark-blue line). The ground stationreceives signals f (cid:48)(cid:48) e and f (cid:48) s at time t at position E (cid:48) . φ is gravitational potential (GP). potential, vector potential, and third- and forth-order terms, Λ f is the correction termsfor ionospheric, tropospheric and tidal effect, O ( c − ) denote high order terms that canbe neglected. Details can be referred to Z. Shen et al. (2017).The theoretical precision of Eq. (3) is at the 10 − level, much better than the orig-inal formula applied in GP-A experiment where its theoretical precision is limited to 10 − .The precision of SFST method for determining GP is about several centimeters, providedthat the stability of OACs can reach 10 − level (Z. Shen et al., 2017). Suppose the GP value of a ground station is given, then we can determine the GPvalues of a satellite by establishing SFST links between them. If the GP distribution (GPD)over a TSS (e.g. GOCE-type or GRACE-type satellite) is determined, the Earth’s ex-ternal gravity field can be determined correspondingly (see Sect. 4). However, since the –5–anuscript submitted to
JGR: Solid Earth orbit of a satellite for the purpose of determining the gravity field is relatively close toground (e.g., the height of GRACE satellite is about 500 km), only a short arc lengthof the orbit is visible to a certain ground station. If we want to determine the GPD overthe TDS, hundreds of ground datum stations with given GP values are needed in orderto guarantee that the satellite can connect to at least 1 ground station at any time (Z. Shen,Shen, & Zhang, 2018), which is impractical for the foreseeable future.Although the SFST method described in Sect. 2.1 was originally designed for de-termining the GP difference between a satellite and a ground site, it can also be usedfor determining the GP difference between two satellites after some modification (Z. Shen,Shen, & Zhang, 2019). Suppose a TS (with low earth orbit) is connected to a GS (withhigh earth orbit), then the setup of the SFST links between them are depicted as Fig.2.
Figure 2.
Geostationary satellite ( GS ) emits a frequency signal f G at time t , denoted byright arrow (blue line). Target satellite ( T S ) transmits the received signal f (cid:48) G (the left arrowdenoted by blue line) and emits a new frequency signal f T at time t (the left arrow denoted bypurple line). The geostationary satellite receives signals f (cid:48)(cid:48) G and f (cid:48) T at time t at position GS (cid:48) . φ denotes gravitational potential. An emitter of the geostationary satellite GS emits a frequency signal f G at time t . When the signal is received by the target satellite T S at time t , it immediately trans-mits the received signal f (cid:48) G and emits a frequency signal f T simultaneously. These twosignals transmitted and emitted from the satellite are received by a receiver at geosta-tionary satellite GS at time t , which are noted as f (cid:48)(cid:48) G and f (cid:48) T , respectively. During theperiod of the emitting and receiving, the position of the geostationary satellite in spacehas been changed from GS to GS (cid:48) . Since the target satellite transmits and emits sig-nals at the same instant it receives signal; its position in the signal links is supposed tobe the point T S at time t . If we set f G = f T , the gravitational difference between the GS and T S can be expressed as∆ φ GT c ≡ φ T − φ G c = ∆ ff G − v T − v G c − (cid:88) i =1 q ( i ) + Λ f + O ( c − ) , (4) –6–anuscript submitted to JGR: Solid Earth where the foot mark G and T denote the GS and T S respectively, and the beat frequency∆ f is given by. ∆ ff G = f (cid:48) T − f T f G − ( f (cid:48)(cid:48) G − f (cid:48) G ) + ( f (cid:48) G − f G )2 f G . (5)Note that there might be a small amount of latency during the transmitting, therebythe positions of the satellite is slightly different at the time it receives and emits signals.Suppose the delay of the signal transponder is about 800 ns (Pierno & Varasi, 2013), andthe orbit height of TS is about 500 km (the velocity is about 7.6 km/s); then the satel-lite moves only 0.62 mm between receiving and emitting signals, and this influence canbe neglected for the SFST links (Z. Shen et al., 2016).Similar to the satellite-to-ground link, if the GP difference between the GS and TSare measured by SFST method, and the absolute GP values of the GS is given, then theGP values of the TS can be derived. The height of a GS is about 35790 km above theequator. If the TS is a low orbit satellite satellite whose height is about 500 km abovethe geoid (e.g., the case of GRACE-FO satellite), then the GS can cover more than halfof the TS’s orbit sphere, as depicted in Fig. 3. Therefore only two evenly distributed GSis sufficient for incessantly SFST links between a GS and the TS, and the GP values ofthe TS’s orbit sphere can be determined correspondingly. However, in practice it is bet-ter to apply three evenly distributed GSs for more reliable and stable connections, as shownin Fig. 5 Figure 3.
Suppose the height of TS is 500 km above the geoid and the Earth’s radius R isabout 6370 km, then the distance between Earth’s center and the target satellite is about 6870km. The distance between Earth’s center and a geostationary satellite is about 42170 km. Whenthe TS is located at the poles (N or S), it can connect the GS without being blocked by theEarth, for the TS is visible until it reach the P point (if block threshold OH = 6400 km, theangle of PON θ = 12 . ◦ ). Therefore the two satellites are inter-visible for more than half of theorbit period of the target satellite, and consequently two evenly distributed GSs are sufficient forincessantly SFST links between a GS and the TS. If the GPD over the TSS is given, one can derive the gravity field outside the TSS,and according to the spherical harmonic expansion formula, the determined gravity fieldoutside the TSS can also be expanded to the Earth’s surface (Heiskanen & Moritz, 1967).
In this section we conducted several simulation experiments to verify the SFST methodfor satellite gravity model establishment. Currently, the most precise atomic clock on-board a satellite is only about 10 − τ − / ( τ in second) in stability (Laurent, Masson- –7–anuscript submitted to JGR: Solid Earth net, Cacciapuoti, & Salomon, 2015; Liu et al., 2018). While the best optical atomic clockon ground have reached the stability of 4 . × − τ − / (Oelker et al., 2019). In theprospect of much better clocks onboard satellites in the future, our experiments will adoptdifferent clock stability levels from 10 − τ − / to 10 − τ − / ; and the results can showus the minimum requirements of clock stability for a satellite gravity model in a certainprecision.The scheme of a simulation experiment is comparing a prior satellite gravity modelto a recovered one, as depicted in Fig. 4; details are explained in the following subsec-tions. Figure 4.
The scheme of the simulation experiment.
In Sect. 2.2 we have shown that two GSs are sufficient for incessant SFST links toa TS. But in practice it is more reliable to adopt three evenly distributed GSs above theequator. Therefore in our experiments we chose the meteorological satellite METEOSAT-9 of EU (at 9.2 ◦ E), the communication satellite CHINASAT-1A of China (at 130.0 ◦ E),and the communication satellite ECHOSTAR-10 of US (at 110.2 ◦ W) as the GSs; andthe GRACE-FO 1 satellite (orbit height is about 500 km) as the TS. The setup of theexperiment is depicted in Fig 5.The orbital period of TS (GRACE-FO 1) is about 1.6 h, and the inclination of TSis about 89 ◦ . For the purpose of obtaining the GP data over its orbit sphere with a res-olution of 30 (cid:48) × (cid:48) , the observations should continue for at least 576 hours (24.0 day),with observation interval being smaller than 8.0 second. Therefore, we set the observa-tional time span as 30.0 days, and the observation interval as 1 second, which fully sat-isfies the resolution requirements. For every observation the TS is connected to a near-est visible GS with SFST links. For the purpose of simulation, the orbits data of thesefour satellites (one TS and three GSs) can be generated from two-line element set (TLE)data by simplified general perturbations models 4 (SGP4) (Croitoru & Oancea, 2016).The GPs at orbits of these satellites can be calculated by EGM2008 model (Pavlis, Holmes,Kenyon, & others, 2012). Then the GP difference between the TS and one proper GScan be obtained. These data are all regarded as true values, hence the errors of orbit dataand gravitational potential model EGM2008 are not considered.The frequency of a microwave signal will be affected by ionosphere and troposphereenvironment. We adopt the International Reference Ionosphere Model (Bilitza et al., 2017;Rawer, Bilitza, & Ramakrishnan, 1978) to obtain the electron density values to estimatethe ionospheric influence (Namazov, Novikov, & Khmel’nitskii, 1975). Since the height –8–anuscript submitted to JGR: Solid Earth
Figure 5.
The TS (target satellite) is connected with 3 GSs (geostationary satellites). Whenthe TS flies around the Earth, it is continuously connected with one or two GSs via SFST links.Given the gravitational potentials at orbits of GSs, the gravitational potential at the orbit of TScan be obtained. of TS is about 500 km which is much higher than the troposphere layer (typically fromground to 60 km height), the influence of troposphere can be neglected. The GP at thesatellites’ orbits will also be influenced by periodical tidal effects, which are well mod-eled (Voigt et al., 2017) and can be removed by some mature softwares such as ETERNA(Wenzel, 1996) or Tsoft (Van Camp & Vauterin, 2005). In our experiment we use ETERNAto generate and analyze tide signals. These tidal signals also include the influences ofother planets (such as Venus, Jupiter etc.) besides the Sun and the Moon.
After setting the input data, the next step focuses on determining the GP valuesat TS’s orbit. There are 3 GSs, denoted as GS i ( i = 1 , ,
3) respectively. As we takea SFST measurement (every 1 second), we can obtain an observed GP difference value∆ ˆ φ GT i ( t ) according to Eq. (4). If the GP of GS i ( t ) is given, the observed GP ˆ φ T ( t ) canbe derived as ˆ φ T ( t ) = GS i ( t ) − ∆ ˆ φ GT i ( t ). –9–anuscript submitted to JGR: Solid Earth
Table 1.
The input datas used in simulation experiments.
Entities Values of ParametersGS Satellite METEOSAT-9CHINASAT-1AECHOSTAR-10TS Satellite GRACE FO 1Gravity field model EGM2008Ionospheric model International Reference IonosphereTide correction ETERNAObservation duration Jan 01 ∼ Jan 30, 2020Mearsurement interval 1 s
The observed values ˆ φ T ( t ) are different from true GP value φ T ( t ) because they areinfluenced by various error sources. In this simulation experiment we have consideredclock error e clk , ionosphere residual error e ion , satellite’s position and velocity errors e pos and e vel , GS’s potential errors e pot and tidal correction residual error e tide . The abovementioned various errors are considered as noises, which are added to the true values.The total errors e all are expressed in the following form e all = e clk + e ion + e pos + e vel + e pot + e tide , (6)and the observed values ˆ φ T ( t ) can be expressed asˆ φ T ( t ) c = φ G ( t ) c + ∆ f ( t ) f G − v T ( t ) − v G ( t ) c − (cid:88) i =1 q ( i ) + Λ f ( t ) + e all ( t ) , (7)The magnitude and behavior of each kind of error play important role in this experiment;thereby we need to investigate different error models based on different error sources tomake the simulation case more close to the real case.We first set the clock error magnitude of e clk as 1 . × τ − / , which is achiev-able currently. Considering the present best clocks with stability of 1 × − in sev-eral hours, we reduce the magnitude of e clk to 1 . × − τ − / and 1 . × − τ − / respectively to improve the observations. Although there are many kinds of random noisesthat affect Atomic clocks’ signals (Major, 2013), the most prominent components are whitefrequency modulation and random walk frequency modulation (Galleani, Sacerdote, Tavella,& Zucca, 2003). Correspondingly the behaviors of clock errors are modeled as follow-ing equation e clk ( t ) = a clk + b clk · t + c clk · φ ( t ) + d clk · (cid:90) t ξ ( t ) dt, (8)where a clk , b clk , c clk and d clk are constant coefficients, φ ( t ) and ξ ( t ) are both standardwhite Gaussian noises. Each term in the right side of Eq.(8) has clear physical mean-ing; specifically a clk denotes the initial frequency difference, b clk · t is the drift term, c clk · φ ( t ) is the white noise component, and d clk · (cid:82) t ξ ( t ) dt represents the random walk ef-fect. As we set proper values of constant coefficients in accordance with the performanceof OACs in Oelker et al. (2019), a series of frequency comparison data with errors em-bedded can be generated. The statistic property of three clock error series are shown inFig. 6.Other error sources are discussed in detail in Z. Shen et al. (2017). Although Z. Shenet al. (2017) focus on the satellite to ground case, we can use similar methods to ana- –10–anuscript submitted to JGR: Solid Earth
Figure 6.
The total Allan deviation for three different clocks. The instabilities of the clocksare about 10 − τ − / , 10 − τ − / and 10 − τ − / for case 1, case 2 and case 3 respectively. lyze the errors in the satellite to satellite case, and most of them demonstrate the samemagnitude. The magnitudes of these error sources are listed in Table 2.As for the mathematical model of these errors, we adopt a general error model whichcontains systematic (initial) offset, drift and white Gaussian noises for each of the er-ror source, expressed as the following equation e j ( t ) = a j + b j · t + c j · φ i ( t ) , ( j = ion, tro, pos, vel, tide, asy ) (9)where a j , b j and c j are constant coefficients, which are randomly set in accordance withthe error magnitudes listed in Table 2According to Eqs. (8) and (9), we can generate the noise signals e all ( t ) term in Eq.(7) based on the magnitudes and nature of the error sources at any time. Noted that thefirst 4 terms or the right side of Eq. (7) are true values, the 5th terms (Λ f ( t )) is the cor-rections of ionosphere and tide effect. The values of these 5 terms can be directly cal-culated. Therefore we can get a set of relevant ”Observed” values, which constitute timeseries of TS’s gravitational potential ˆ φ T ( t ) as the left side of Eq. (7). Since the TS fliesover the whole Earth in a period of about 30 d, these values are corresponding to thegravitational potential at different time points on the TS’s orbits and we obtain a setof values ˆ φ T ( x, y, z ) related to orbit data. After a continuous observation of 30.0 days (720.0 hours), there are 324,000 observedGP points ˆ φ T ( x, y, z ) distributed over the TSS enclosing the Earth (see Fig. 7). The cor-responding disturbing potentials of the three different experiment cases are depicted inFig. 8, where the first subfigure demonstrates the disturbing potential of EGM2008 as –11–anuscript submitted to JGR: Solid Earth
Table 2.
Error magnitudes of different error sources in determining the gravitational potentialdifference between a satellite and a ground station. They are transformed to relative frequency.Details can be referred to Z. Shen et al. (2017).
Influence factor (Residual) Error magnitude in ∆ f/f e ionospheric correction residual δf ion ∼ − tidal correction residual δf tide ∼ − position & velocity δf vepo ∼ − (10 mm and 0.1mm/s a )clock error δf osc ∼ − τ − / a Satellite’s position errorsare assumed as 10 mm (Kang et al., 2006), velocity errors are assumed as 0.1mm/s (Sharifi, Seif, & Hadi,2013);
Figure 7.
The trace of target satellite (TS) in Earth-Centered and Earth-Fixed (ECEF)coordinate for (a) (b) true values. we can see that the observe results of case 1 (clock instability of 10 − τ − / )seems to be useless, but the observe results of case 3 (clock instability of 10 − τ − / )are almost identical with the true values. In the next section, we will calculate the spher-ical harmonic expansion coefficients for 3 different recovered Earths gravity field mod-els (REGMs) based on the observed GP values in our simulation experiments. These co-efficients will be compared with the true value of EGM2008 to evaluate their accuracy. Based on the determined GPD over the TDS using frequency links, we can deter-mine the gravitational potential field outside the solid Earth by least-squares (LS) ap-proach. The Earths gravitational potential V at a point ( r, θ, λ ) outside the Earth canbe expanded into a series of spherical harmonics (Hofmann-Wellenhof & Moritz, 2005) V ( r, θ, λ ) = GMa N max (cid:88) n =0 n (cid:88) m =0 (cid:16) ar (cid:17) n +1 (cid:0) C nm cos mλ + S nm sin mλ (cid:1) P nm (cos θ ) , (10)where the spherical coordinates ( r, θ, λ ) represent a 3-D position in the Earth-Centered,Earth-Fixed (ECEF) reference frame, r is the geocentric radius, θ and λ are the spher- –12–anuscript submitted to JGR: Solid Earth
Figure 8.
The disturbing potentials (m s − ) of (a) EGM2008 as true values; (b) experimentcase 1, using clocks with instability of about 10 − τ − / ; (c) experiment case 2, using clocks withinstability of about 10 − τ − / ; (d) experiment case 3, using clocks with instability of about10 − τ − / . ical co-latitude and longitude respectively, GM is the geocentric gravitational constant, a is the semi-major axis of the reference ellipsoid, C nm and S nm the (fully-normalized)geopotential coefficients which describe the external gravitational field of the Earth, P nm are the (fully-normalized) associated Legendre functions of degree n and order m , and N max is the maximum degree of the harmonic expansion.For the linear observation equation Eq. (10), the functional and statistical mod-els of the gravitational field recovery from the GPD observations are defined by a stan-dard Gauss-Markov model as follows: y = Ax + (cid:15) , E { y } = Ax , D { y } = σ Q = σ P − , (11)where y is the vector of GP observations, A is the design matrix, x is the vector of (un-known) geopotential coefficients C nm , S nm to be estimated, (cid:15) is the vector of observa-tion errors, D { y } is the error variance-covariance matrix, P is the weight matrix, Q isthe inverse of the weight matrix, and σ is the variance component.Based on the data processing method described above, we estimated three REGMsup to degree and order 200 from GP values distributed over the TSS in three differentcases, corresponding respectively the clocks instabilities of 10 − τ − / , 10 − τ / and10 − τ − / . Here we set the weight matrix P as a unit matrix by considering that thenoise in GP observations is white noise. The absolute values of the coefficient differences(logarithm representation) between the recovered harmonic expansion coefficients of Earthsgravity field and that of EGM2008 are illustrated globally in Fig. 9 as (a), (b) and (c).The results show that if the clock instability is poorer than 10 − τ − / (as cases(a) and (b)),the precision of recovered Earth’s gravity fields is poor. When the clock in- –13–anuscript submitted to JGR: Solid Earth
Figure 9.
The offset between recovered harmonic expansion coefficients and that ofEGM2008. The clock instabilities of different experiments are about (a) − τ − / ; (b) − τ − / ; (c) − τ − / . (d) and (e) are another two demonstration experiments which weset the total error magnitude for a SFST link is about 10 − and 10 − respectively. (f ) denotesthe recovered coefficients based on true values of GPD over the TDS, determined by EGM2008–14–anuscript submitted to JGR: Solid Earth
Table 3.
The statistic information of the GP offset at target satellite’s orbit between valuescalculated by EGM2008 (true values) and by recovered Earths gravity field models (REGMs).
Case Clock precision Mean offset (m / s ) STD (m / s )(a) 10 − τ − / − τ − / − τ − / -0.0009 0.1142(d) 10 − τ − / -0.0001 0.0114(e) 10 − τ − / -2.5e-6 0.0013(f) NA 2.4e-8 0.0006 stability reach 10 − τ − / as case (c), we can obtain a fairly good coefficients for ordersand degrees lower than 50. In addition, the zonal and near-zonal coefficients of REGMsare worse than other kinds of coefficients when the clocks instability reaches 10 − τ − / .This is due to the fact that we did not use any regularized technique to deal with theill-posed problem caused by the polar gap of GOCE mission. However, even if we usethe regularized technique to deal with the ill-posed problem, the above mentioned prob-lem still exists (Baur et al., 2014).The performance of REGMs can also be evaluated by the calculated GPs at theTSS (TS-orbit defined spherical surface). The mean offset and standard deviation (STD)of the calculated GP distribution difference between EGM2008 and REGMs over TSSare illustrated in Table 3. The STDs of GPD over the TSS are respectively 1135.3232m / s , 11.3758 m / s and 0.114 m / s for cases (a), (b) and (c).If the clock instabilities can be improved to 10 − τ − / or even 10 − τ − / level,some other error sources, such as ionospheric residual and velocity error, will take dom-inant and become the bottle neck for the precision of recovered coefficients. In that case,further detailed analysis or corrections for various error sources is required for establish-ing better correction models. Since there might be an extended period for us to set on-board clocks of 10 − τ − / instability level, in this paper we will leave that researchesfor future works. However, in order to show the potential of this method, we conductedtwo simplified experiments that The total error (the sum of clock errors and various othererror sources) of SFST links are set to 10 − and 10 − respectively as illustrated in Fig.9 (d) and (e). We can see that if the total error magnitudes can reduced to 10 − , therecovered harmonic expansion coefficients show fairly good quality, close to that recov-ered from the true values as illustrated in Fig. 9 (f). The REGM’s precision of case (d),(e) and (f) are shown in table 3. In this paper we formulated an alternative method to determine satellite gravityfield based on precise clocks and frequency signals transfer. It is a new application of gen-eral relativistic theory in geodesy, and the gravity field can be determined at the pre-cision levels of about 10 m / s , 10 m / s and 10 − m / s , given the clock stabilitiesof 10 − τ − / , 10 − τ − / and 10 − τ − / respectively. Currently the stability of a satel-lite’s onboard clock is about 10 − τ − / , and it is the main error influence for GP de-termination and EGM establishment. However, precise optical atomic clocks have reached10 − τ − / level under laboratory environments (Oelker et al., 2019). It is foreseeablethat in the near future the stability of onboard atomic clocks can achieve a similar level,and the precision of the EGM established by intersatellite SFST method can reach 1 cm –15–anuscript submitted to JGR: Solid Earth level. Compared to the conventionally used methods of establishing satellite gravity modelsuch as using gravimeter and gravity gradiometer to measure the first-order and second-order derivative of potential, the SFST method may directly determine the GPs, sim-plifying in some sense the estimation of the harmonic coefficients of Earth’s gravity field.According to this study, once the onboard clocks’ stabilities reach the level of 10 − τ − / ,the relativistic method will be applicable for high precision satellite gravity field deter-mination, and can be used to provide a decimeter level Earth gravity model with a res-olution of around 1 ◦ × ◦ . If clock stabilities are better than 10 − τ − / (10 − τ − / or even 10 − τ − / for instance), various other error sources (ionosphere and tropospherecorrection residual errors, satellite position errors, et al.) will be the bottle-neck for de-termining a precise Earth’s gravity field, and more precise error correction models needto be established. Acknowledgments
This study is supported by National Natural Science Foundation of China (NSFC) (grantNos. 41721003, 41631072, 41874023, 41804012, 41429401, 41574007), and Natural Sci-ence Foundation of Hubei Province (grant No. 2019CFB611).
References
Altschul, B., Bailey, Q. G., Blanchet, L., Bongs, K., Bouyer, P., Cacciapuoti, L., . . .Wolf, P. (2014, January). Quantum tests of the einstein equivalence principlewith the STE-QUEST space mission.
Adv. Space Res. , (1), 501–524.Baur, O., Bock, H., H¨ock, E., J¨aggi, A., Krauss, S., Mayer-G¨urr, T., . . . Zehentner,N. (2014, October). Comparison of GOCE-GPS gravity fields derived bydifferent approaches. J. Geodesy , (10), 959–973.Bilitza, D., Altadill, D., Truhlik, V., Shubin, V., Galkin, I., Reinisch, B., & Huang,X. (2017, February). International reference ionosphere 2016: From ionosphericclimate to real-time weather predictions. Space Weather , (2), 418–429.Bjerhammar, A. (1985, January). On a relativistic geodesy. Bull. Am. Assoc. Hist.Nurs. , (3), 207–220.Bock, H., J¨aggi, A., Meyer, U., Visser, P., van den IJssel, J., van Helleputte, T., . . .Hugentobler, U. (2011, May). GPS-derived orbits for the GOCE satellite. J.Geodesy , (11), 807–818.Croitoru, E.-I., & Oancea, G. (2016, June). Satellite tracking using norad two-lineelement set format. AFASES 2016 , (1), 423–432.Deschˆenes, J.-D., Sinclair, L. C., Giorgetta, F. R., Swann, W. C., Baumann, E.,Bergeron, H., . . . Newbury, N. R. (2016, May). Synchronization of distantoptical clocks at the femtosecond level. Phys. Rev. X , (2), 021016.Ditmar, P., & Sluijs, A. A. v. E. v. d. (2004, September). A technique for model-ing the earth’s gravity field on the basis of satellite accelerations. J. Geodesy , (1), 12–33.Drinkwater, M. R., Floberghagen, R., Haagmans, R., Muzi, D., & Popescu, A.(2003). GOCE: ESA’s first earth explorer core mission. In G. Beutler,M. R. Drinkwater, R. Rummel, & R. Von Steiger (Eds.), Earth gravity fieldfrom space — from sensors to earth sciences: Proceedings of an ISSI work-shop 11–15 march 2002, bern, switzerland (pp. 419–432). Dordrecht: SpringerNetherlands.Dziewonski, A. M., & Anderson, D. L. (1981, June). Preliminary reference earthmodel.
Phys. Earth Planet. Inter. , (4), 297–356.Einstein, A. (1915, January). Die feldgleichungen der gravitation. Sitzungsberichteder K¨oniglich Preußischen Akademie der Wissenschaften (Berlin), Seite 844-847. , , 844–847.Flechtner, F., Neumayer, K.-H., Dahle, C., Dobslaw, H., Fagiolini, E., Raimondo, –16–anuscript submitted to JGR: Solid Earth
J.-C., & G¨untner, A. (2016, March). What can be expected from the GRACE-FO laser ranging interferometer for earth science applications?
Surv. Geophys. , (2), 453–470.Flury, J. (2016, July). Relativistic geodesy. J. Phys. Conf. Ser. , (1), 012051.Galleani, L., Sacerdote, L., Tavella, P., & Zucca, C. (2003, June). A mathematicalmodel for the atomic clock error. Metrologia , (3), S257.Grotti, J., Koller, S., Vogt, S., H¨afner, S., Sterr, U., Lisdat, C., . . . Calonico, D.(2018, May). Geodesy and metrology with a transportable optical clock. Nat.Phys. , (5), 437–441.Han, S.-C., Jekeli, C., & Shum, C. K. (2002, August). Efficient gravity field recov-ery using in situ disturbing potential observables from CHAMP: GRAVITYFIELD RECOVERY FROM CHAMP. Geophys. Res. Lett. , (16), 36–1–36–4.Hannig, S., Pelzer, L., Scharnhorst, N., Kramer, J., Stepanova, M., Xu, Z. T., . . .Schmidt, P. O. (2019, May). Towards a transportable aluminium ion quantumlogic optical clock. Rev. Sci. Instrum. , (5), 053204.Heiskanen, W. A., & Moritz, H. (1967, December). Physical geodesy. Bull. Geode-sique , (1), 491–492.Hirt, C., Kuhn, M., Featherstone, W. E., & G¨ottl, F. (2012, May). Topo-graphic/isostatic evaluation of new-generation GOCE gravity field models:TOPOISOSTATIC EVALUATION OF GOCE GRAVITY. J. Geophys. Res. , (B5).Hofmann-Wellenhof, B., & Moritz, H. (2005). Physical geodesy . Springer.Huang, Y., Guan, H., Zeng, M., Tang, L., & Gao, K. (2019, January). Ca + ion op-tical clock with micromotion-induced shifts below 10 − Phys. Rev. A , (1),011401.Hwang, C. (2001, May). Gravity recovery using COSMIC GPS data: application oforbital perturbation theory. J. Geodesy , (2), 117–136.Jekeli, C. (1999, October). The determination of gravitational potential differencesfrom satellite-to-satellite tracking. Celest. Mech. Dyn. Astron. , (2), 85–101.Kang, Z., Tapley, B., Bettadpur, S., Ries, J., Nagel, P., & Pastor, R. (2006, Jan-uary). Precise orbit determination for the GRACE mission using only GPSdata. J. Geodesy , (6), 322–331.Kaula, W. M. (1966). Theory of satellite geodesy: Applications of satellites togeodesy . Dover Publications.Kleppner, D., Vessot, R. F. C., & Ramsey, N. F. (1970, January). An orbiting clockexperiment to determine the gravitational red shift.
Astrophys. Space Sci. , (1), 13–32.Kopeikin, S. M., Kanushin, V. F., Karpik, A. P., Tolstikov, A. S., Gienko, E. G.,Goldobin, D. N., . . . Hanikova, E. A. (2016, July). Chronometric measure-ment of orthometric height differences by means of atomic clocks. GravitationCosmol. , (3), 234–244.Kornfeld, R. P., Arnold, B. W., Gross, M. A., Dahya, N. T., Klipstein, W. M., Gath,P. F., & Bettadpur, S. (2019, May). GRACE-FO: The gravity recovery andclimate experiment Follow-On mission. J. Spacecr. Rockets , (3), 931–951.Kvas, A., & Mayer-G¨urr, T. (2019, December). GRACE gravity field recovery withbackground model uncertainties. J. Geodesy , (12), 2543–2552.Laurent, P., Massonnet, D., Cacciapuoti, L., & Salomon, C. (2015). TheACES/PHARAO space mission. C. R. Phys. , (5), 540–552.Liu, L., L¨u, D.-S., Chen, W.-B., Li, T., Qu, Q.-Z., Wang, B., . . . Wang, Y.-Z. (2018,July). In-orbit operation of an atomic clock based on laser-cooled 87rb atoms. Nat. Commun. , (1), 2760.Major, F. G. (2013). The quantum beat: The physical principles of atomic clocks .Springer Science & Business Media.McGrew, W. F., Zhang, X., Fasano, R. J., Sch¨affer, S. A., Beloy, K., Nicolodi, D., –17–anuscript submitted to
JGR: Solid Earth . . . Ludlow, A. D. (2018, December). Atomic clock performance enablinggeodesy below the centimetre level.
Nature , (7734), 87–90.Mehlst¨aubler, T. E., Grosche, G., Lisdat, C., Schmidt, P. O., & Denker, H. (2018,June). Atomic clocks for geodesy. Rep. Prog. Phys. , (6), 064401.Namazov, S. A., Novikov, V. D., & Khmel’nitskii, I. A. (1975, April). Doppler fre-quency shift during ionospheric propagation of decameter radio waves (review). Radiophys. Quantum Electron. , (4), 345–364.Oelker, E., Hutson, R. B., Kennedy, C. J., Sonderhouse, L., Bothwell, T., Goban, A.,. . . Ye, J. (2019, October). Demonstration of 4 . × stability at 1 s for twoindependent optical clocks. Nat. Photonics , (10), 714–719.Pail, R., Bruinsma, S., Migliaccio, F., F¨orste, C., Goiginger, H., Schuh, W.-D., . . .Tscherning, C. C. (2011, November). First GOCE gravity field models derivedby three different approaches. J Geod , (11), 819.Pail, R., Yeh, H.-C., Feng, W., Hauk, M., Purkhauser, A., Wang, C., . . . Xu, H.(2019, November). Next-Generation gravity missions: Sino-European numeri-cal simulation comparison exercise. Remote Sensing , (22), 2654.Pavlis, N. K., Holmes, S. A., Kenyon, S. C., & others. (2012). The development andevaluation of the earth gravitational model 2008 (EGM2008). Journal of Geo-physical Research: Solid Earth , (B4), B04406.Pierno, L., & Varasi, M. (2013, June). Switchable delays optical fibre transponderwith optical generation of doppler shift (No. 8466831).Puetzfeld, D., & L¨ammerzahl, C. (Eds.). (2019).
Relativistic geodesy: Foundationsand applications . Springer, Cham.Rawer, K., Bilitza, D., & Ramakrishnan, S. (1978). Goals and status of the interna-tional reference ionosphere.
Rev. Geophys. , (2), 177.Reigber, C., Balmino, G., Schwintzer, P., Biancale, R., Bode, A., Lemoine, J.-M.,. . . Zhu, S. Y. (2002, July). A high-quality global gravity field model fromCHAMP GPS tracking data and accelerometry (EIGEN-1S): A GLOBALGRAVITY FIELD MODEL. Geophys. Res. Lett. , (14), 37–1–37–4.Reubelt, T., Austen, G., & Grafarend, E. W. (2003, August). Harmonic analysisof the earth’s gravitational field by means of semi-continuous ephemerides ofa low earth orbiting GPS-tracked satellite. case study: CHAMP. J. Geodesy , (5), 257–278.Schiller, S., Gorlitz, A., Nevsky, A., Alighanbari, S., Vasilyev, S., Abou-Jaoudeh, C.,. . . Levi, F. (2012, April). The space optical clocks project: Development ofhigh-performance transportable and breadboard optical clocks and advancedsubsystems. In (pp. 412–418). IEEE.Sharifi, M. A., Seif, M. R., & Hadi, M. A. (2013). A comparison between numericaldifferentiation and kalman filtering for a leo satellite velocity determination. Artificial Satellites , (3), 103–110.Shen, W., Ning, J., Liu, J., Li, J., & Chao, D. (2011, January). Determination of thegeopotential and orthometric height based on frequency shift equation. Nat.Sci. , (5), 388–396.Shen, Z., Shen, W., & Zhang, S. (2018, April). Determination of the gravitationalpotential at GOCE-type satellite orbit using frequency signal transmissionapproach. In (p. 19415).Shen, Z., Shen, W., & Zhang, S. (2019, April). Determination of gravitational po-tential distribution over a geocentric quasi-sphere based on satellite-to-satellitefrequency transmitting. In (p. 4591).Shen, Z., & Shen, W.-B. (2017). Determination of gravitational potential distri-bution over a geocentric quasi- sphere based on links between GRACE- andGNSS-type satellites. In (p. 11250).Shen, Z., Shen, W.-B., Peng, Z., Liu, T., Zhang, S., & Chao, D. (2019, April).Formulation of determining the gravity potential difference using Ultra-Highprecise clocks via optical fiber frequency transfer technique. J. Earth Sci. , –18–anuscript submitted to JGR: Solid Earth (2), 422–428.Shen, Z., Shen, W.-B., & Zhang, S. (2016, August). Formulation of geopotentialdifference determination using optical-atomic clocks onboard satellites andon ground based on doppler cancellation system. Geophys. J. Int. , (2),1162–1168.Shen, Z., Shen, W.-B., & Zhang, S. (2017, July). Determination of gravitationalpotential at ground using Optical-Atomic clocks on board satellites and onground stations and relevant simulation experiments. Surv. Geophys. , (4),757–780.Takano, T., Takamoto, M., Ushijima, I., Ohmae, N., Akatsuka, T., Yamaguchi, A.,. . . Katori, H. (2016, August). Geopotential measurements with synchronouslylinked optical lattice clocks. Nat. Photonics , (10), 662.Tapley, B. D., Bettadpur, S., Watkins, M., & Reigber, C. (2004, May). The gravityrecovery and climate experiment: Mission overview and early results: GRACEMISSION OVERVIEW AND EARLY RESULTS. Geophys. Res. Lett. , (9),L09607.Van Camp, M., & Vauterin, P. (2005, June). Tsoft: graphical and interactive soft-ware for the analysis of time series and earth tides. Comput. Geosci. , (5),631–640.Vessot, R. F. C., & Levine, M. W. (1979, February). A test of the equivalence prin-ciple using a space-borne clock. Gen. Relat. Grav. , (3), 181–204.Vessot, R. F. C., Levine, M. W., Mattison, E. M., Blomberg, E. L., Hoffman, T. E.,Nystrom, G. U., . . . Wills, F. D. (1980, December). Test of relativistic gravita-tion with a Space-Borne hydrogen maser. Phys. Rev. Lett. , (26), 2081–2084.Visser, P. N. A. M., Sneeuw, N., & Gerlach, C. (2003, June). Energy integralmethod for gravity field determination from satellite orbit coordinates. J.Geodesy , (3), 207–216.Voigt, C., F¨orste, C., Wziontek, H., Crossley, D., Meurers, B., P´alink´aˇs, V., . . . Sun,H. (2017). The data base of the international geodynamics and earth tideservice (IGETS). In (p. 4947).Weinberg, S. (1972). Gravitation and cosmology: Principles and applications of thegeneral theory of relativity . New York: Wiley.Wenzel, H.-G. (1996). The nanogal software: Earth tide data processing packageETERNA 3.30.
Bull. Inf. Mar´ees Terrestres , , 9425–9439.Wu, H., M¨uller, J., & L¨ammerzahl, C. (2019, March). Clock networks for height sys-tem unification: a simulation study. Geophys. J. Int. , (3), 1594–1607.(3), 1594–1607.