学術雑誌掲載論文等

gravimetry

A uthor(s )

X u, Peiliang

C itation

C ommunications in Nonlinear S cience and Numerical

S imulation (2018), 59: 515-543

Is s ue D ate

2018-06-01

UR L

http://hdl.handle.net/2433/230356

R ig ht

©

2017 T he A uthor(s). Published by E lsevier B .V . T his is an

open access article under the C C B Y -NC -ND

license.(http://creativecommons.org/licenses/by-nc-nd/4.0/)

T ype

J ournal A rticle

T extvers ion

publisher

ContentslistsavailableatScienceDirect

Commun

Nonlinear

Sci

Numer

Simulat

journalhomepage:www.elsevier.com/locate/cnsns

Research

paper

Measurement-based

perturbation

theory

and

differential

equation

parameter

estimation

with

applications

to

satellite

gravimetry

Peiliang

Xu

DisasterPreventionResearchInstitute,KyotoUniversity,Uji,Kyoto611-0011,Japan

a

r

t

i

c

l

e

i

n

f

o

Articlehistory:

Received3August2016 Revised18August2017 Accepted23November2017 Availableonline2December2017

Keywords:

Differentialequationparameterestimation Earth’sgravityfield

Satellitegravimetry

Measurement-basedperturbation Conditionadjustmentwithparameters Nonlineardifferentialequations NonlinearVolterra’sintegralequations

a

b

s

t

r

a

c

t

The numerical integration method has been routinely used by major institutions world- wide, for example, NASA Goddard Space Flight Center and German Research Center for Geosciences (GFZ), to produce global gravitational models from satellite tracking mea- surements of CHAMP and/or GRACE types. Such Earth’s gravitational products have found widest possible multidisciplinary applications in Earth Sciences. The method is essentially implemented by solving the differential equations of the partial derivatives of the orbit of a satellite with respect to the unknown harmonic coefficients under the conditions of zero initial values. From the mathematical and statistical point of view, satellite gravimetry from satellite tracking is essentially the problem of estimating unknown parameters in the New- ton’s nonlinear differential equations from satellite tracking measurements. We prove that zero initial values for the partial derivatives are incorrect mathematically and not permit- ted physically. The numerical integration method, as currently implemented and used in mathematics and statistics, chemistry and physics, and satellite gravimetry, is groundless, mathematically and physically. Given the Newton’s nonlinear governing differential equa- tions of satellite motion with unknown equation parameters and unknown initial condi- tions, we develop three methods to derive new local solutions around a nominal reference orbit, which are linked to measurements to estimate the unknown corrections to approx- imate values of the unknown parameters and the unknown initial conditions. Bearing in mind that satellite orbits can now be tracked almost continuously at unprecedented accu- racy, we propose the measurement-based perturbation theory and derive global uniformly convergent solutions to the Newton’s nonlinear governing differential equations of satel- lite motion for the next generation of global gravitational models. Since the solutions are global uniformly convergent, theoretically speaking, they are able to extract smallest possi- ble gravitational signals from modern and future satellite tracking measurements, leading to the production of global high-precision, high-resolution gravitational models. By directly turning the nonlinear differential equations of satellite motion into the nonlinear integral equations, and recognizing the fact that satellite orbits are measured with random errors, we further reformulate the links between satellite tracking measurements and the global uniformly convergent solutions to the Newton’s governing differential equations as a con- dition adjustment model with unknown parameters, or equivalently, the weighted least squares estimation of unknown differential equation parameters with equality constraints,

E-mailaddress:[email protected]

https://doi.org/10.1016/j.cnsns.2017.11.021

for the reconstruction of global high-precision, high-resolution gravitational models from modern (and future) satellite tracking measurements.

© 2017 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license. ( http://creativecommons.org/licenses/by-nc-nd/4.0/)

1. Introduction

The historyofgeodesyhaschangedwiththelaunch ofthefirstartificial satelliteSputnik-1on 4October1957 bythe formerSovietUnion.Satelliteshavebroughtarevolutionarychangeofthewaywe measuretheEarth, bothgeometrically andphysically.Sixtyyearsofsatellitegravimetryhassincewitnessedaprofoundadvanceinbothsatellitegravitytheoryand practicalproductionofEarth’sgravitationalmodels,inparticular,inthepastdecadeorso,thankstothelaunchesofthree dedicatedsatellite gravity missions:the ChallengingMini-satellitePayload(CHAMP)launched in2000,theGravity Recov-eryandClimateExperiment(GRACE)in2002,andtheGravity fieldandsteady-stateOcean CirculationExplorer(GOCE)in 2009.Thededicatedsatellitegravitymissionshaveprovidedfrontierchallengesandsciencetoolstoexploreandunderstand solidandfluidgeophysicalprocessesanddynamicsoftheEarth(seee.g.,[19,78,87,94,104,109]).From thegeodeticpointof view,themostexpectedcelebrationandcontinuation ofthesesixtyyearsofgreatachievementofsatellitegravimetry may culminateinthelaunchoftheGRACEFollow-onmissionscheduledin2017.FormoreinformationonGRACEFollow-on,the readerisreferredto[15,22]andthewebsitehttps://gracefo.jpl.nasa.gov/.

TocomputetheEarth’sgravitationalfieldfromsatellitetrackingmeasurements,anumberofmathematicalmethodshave beenproposedtoestablishthelinksbetweenmeasurementsandtheforceparametersoftheEarth’sgravitationalfield.The majorclassesofmethodsinclude:(i)linearperturbationmethods;(ii)thedynamicalnumericalintegrationmethod;(iii)the orbit-energy-basedmethod;(iv)thetwo-point(orbital)boundaryvalueproblemtheory,whichwasfirstsolvedbySchneider

[120,121]andfurther developed byIlk et al.[40,41](seealso[74]);(v)theorbit-invertedaccelerationapproach.Theideaof computingaccelerationsfromGPS-derivedcoordinateswasfirstproposedbyJekeliandGarcia[48]forairbornegravimetry andthenappliedtoreconstructtheEarth’sgravitationalfieldfromCHAMPmissionbyReubeltetal.[85](seealso[10,21]); and(vi)satellitegradiometry(seee.g.,[93,95]).Recently,Xu[111,114]proposed a measurement-based perturbation method, whichisgloballyconvergentuniformly.

Inthispaper,wewillbemainly concernedwiththefirsttwotypesofmethodstoreconstructtheEarth’sgravitational field fromsatellite trackingmeasurements,because Earth’sgravitationalmodelswere firstderived by usinglinear pertur-bation methods, andbecause the numericalintegration method hasnow been routinely used to produce globalEarth’s gravitationalmodels.Althoughtheorbit-energy-basedmethodismathematicallyrigorous,asfirstproposedbyBjerhammer

[7,8]andlikelyalsoindependentlyby Wolff[110]andfurthermodifiedbyJekeli[47]inordertoaccountfortechnological advanceofspaceobservation,itisstill notabletofullyutilizetheunprecedentedaccuracyofall modernspace measure-ments (see also [36]). Neither the two-point boundary value problemnor the orbit-inverted accelerationapproach have beenusedbyanymajorinstitutionssuchasNASA GoddardSpaceFlightCenterandGFZ toproduceglobalEarth’s gravita-tionalmodels,though theyareindeedusedtocompute gravitationalmodels fromGRACE andCHAMPmeasurements (see e.g.,[10,21,74,85]).The two-pointboundaryvalueproblemalways hasbeenbasedonshortarcsand, asaresult,cannever fullyutilizeunprecedentedaccuracyandcontinuityofmissionsofGRACE/CHAMPtypestotheirpotentiallimitof measure-menttechnology.Nevertheless,themethodisonlyofquasi-linearaccuracybutiswellknowntobedifficulttoimplement, inparticular, forsparsetracking data.Orbit-inverted accelerationscan beunreliableandinaccurate, sincethe operatorof differentiationis ill-posed in nature.As a direct consequence,position-invertedacceleration signalscan be distorted and their noise canbe substantially amplified, theextent ofwhichdependson thenoise level ofsatellite positions,the time intervalusedtoderive accelerationsandregularization.Wheneverpossible,oneshouldavoidsolvingthisintermediate in-verseill-posedproblemanduseitsdistortedsignalsofaccelerationtofurtherinvertforglobalgravitationalmodels.Inthe caseofsatellitegradiometry,sincewedirectlymeasurethesecondgradientsofthegravitationalpotential,theobservational equationsofthegravitationaltensorsaremathematicallystraightforward([93,95]).

Estimatingunknown differentialequation parameters hasbeen essential in manyareas ofscience and engineering. A keymathematical componentofsuch an estimationprocedurehasbeen tosolvethe deriveddifferential equationsofthe partialderivativeswithrespecttotheunknowndifferentialequationparametersundertheassumptionofzeroinitialvalues forthepartialderivatives,asoriginallypublishedbyGronwall[28]almost100yearsago(seealso[26,37]).Thisestimation techniquehasfound widest possible applications,forexample, inmathematics, statistics, chemistry, physics andsatellite gravimetry,which,however,hasnowbeenbestknownasthedynamicalnumericalintegrationmethodingeodesy.Actually, inthecommunityofgeoscience, thismethodseemedtobefirsthintedatby Anderle[2],likelyindependentofwhathad beenpublished by then, since no referenceswere madeto the mathematical literature of [26,28] and[37]. The method wasthen(re-)publishedinamathematicalpaperbyRileyetal.[90](seealso[3,5,76]),again,withoutanyreferencetothe above-mentionedmathematicalliterature.

Thedynamicalnumericalintegrationmethodhasgainedwidespreadacceptancewithoutchallengewiththepublications ofAnderle[2]andRileyetal.[90],andhassincebeenroutinelyusedbyalmostallmajorinstitutionsworldwidetoproduce globalgravitational models fromsatellite tracking measurements, likely partly attributedto the fact that NASA Goddard SpaceFlightCenterusedandimplementedthisnumericalintegrationideabyAnderle[2]asthemathematicalfoundationto computeEarth’sgravitationalmodels(seee.g.,[67,70]).Amongthemostimportantgravitationalmodelsbeforethededicated satellitegravitymissionsCHAMP,GRACEandGOCEaretheGEMseriesofgravitationalmodelsfromtheGoddardSpaceFlight Center(seee.g.,[72,73])andthosefromthejointGerman–Frenchteam(seee.g.,[6,97]).Forreviewsonprogressinsatellite gravimetrybeforethelaunchesofthededicatedsatellitegravitymissions,thereaderisreferredto[64](seealso[65])for abrief progressreport from1958 to1982 andto[80] for asummarized report,retrospective andprospectiveon gravity observation.

Althoughthenumericalintegration methodhasbeenwidely acceptedbyalmost all majorinstitutions worldwideasa standardmethodformakingglobalgravitationalmodels fromsatellite trackingmeasurements,thisisan unbelievable sci-entificfallacy, since themethod isbased on groundlessclaim andproved to be incorrect, mathematically andphysically

[112,115,116].Unfortunately,themethodhasbeenusedforalmost60yearsbyallimportantinstitutionsworldwideto com-puteglobalgravitationalmodelsfromsatellitetrackingmeasurementsforuseingeodesy,solidgeophysics,oceandynamics, hydrology,interactionof ocean andsurfacewaterwithatmosphere, andfarmore beyond, bearinginmind that no tech-nical documentationis available to providea complete andmathematically rigoroussupport forthe method,to my best knowledge,afteranextensivesearchofliterature.

Inthe1965U.S.NavalWeaponsLaboratorytechnicalreport,Anderle[2]wrotetwosentencesinthesectionofProcedure tomeanthismethodby saying“Numericalintegration...wasusedtocomputetheorbitofthesatellite.Thepartialderivatives ofsatellitepositionwithrespecttoorbitandgravityparameters werealsoobtainedbynumericalintegrationoftheperturbation equations”.Nomathematicalformulationwasgiventoprovidefurthertechnicalsupportandexplanationoftheidea.In con-nectionwiththeideabyAnderle[2],acompletepublicationaboutthemathematicsofthemethodwasgivenbyRileyetal.

[90]fromHughesAircraftCompanyandAerospaceCorporation,whichwassoonwellknownandacceptedamonggeodesists worldwide.AlthoughRileyetal.[90]correctlyderivedthedifferentialequationsoftheorbitandvelocityofasatellitewith respecttodifferentialequationparametersintheirmathematicspaper,theysimplyclaimedatthebeginningofthesecond pageofthepaper“Theinitialvalues...ingeneralwillbezeroif

β

kisadifferentialequationparameter”.Asamatteroffact,the

differentialequationsforthe partialderivativeswithrespectto theequationunknownsandthezeroinitialvaluesforthe partialderivativesdatedbackmuchearliertoGronwall[28]andRitt[91](see also[26,37]).Mathematically,this zero-initial-valuestatementdoesnotderivefromtheoriginaldifferentialequationsbutisnothingmorethanaclaim.Obviously,akey elementtodecidetheparticularsolutionofdifferentialequationsisclaimedwithoutanymathematical/physicaljustification. InthetechnicalreportonGoddardEarthmodels(5and6),Lerchetal.[67]reportedthattheyusedthenumerical integra-tionmethodfromtheideaofAnderle[2].Moreprecisely,intheappendixoforbittheoryforthesoftwaresystemGEODYN, Lerchetal.[67] wrote“Thepartialderivatives...areobtained bydirect numericalintegrationofthevariationalequations” on pageA1–10and“Initially,...therestofthematrix(correspondingtothepartialderivatives– notesaddedbytheauthor)iszero” onpage A1–22. To support theproduction ofglobal gravitationalmodels fromsatellite tracking measurements, Goddard SpaceFlight Center,together withComputerSciences Corporation,preparedalengthy technicalreportofabout700pages

[70].InSection 6.1.4 onpage 6–11,although Long etal.[70]correctlyrealized thattheir differential Eqs.(6–49)required initialconditions,andeventhoughtheinitialvaluesare thekey tosolvethedifferentialequations,theydidnottouchthe issueofhow to determine theinitial values butchose toshow howto numericallysolve the differential equationsasif theinitialvalueshadbeengiven.ThenumericalintegrationmethodwasfollowedbythejointGerman/Frenchteamaswell (seee.g., [86]). Recently,Xu [112,115,116]mathematically provedrigorouslythat assigningzerovaluesto theinitial partial derivativesviolatesthephysicsofmotionofcelestialbodies.MorewillcomeinSection2.2.

Profoundtechnologicaladvancesinspaceobservationhavebeenachievedincomparisonwiththosein1950sand1960s whenlinearperturbationsweresignificantlydevelopedandthenumericalintegrationmethodwashintedatandpublished. Twomostimportantfeaturesoftheseadvancesare: (i)lowEarthorbiting(LEO)satellitescannowbetrackedanddirectly measured,preciselyandalmostcontinuously,byusingGlobalNavigationSatelliteSystems(GNSS).Inotherwords,wehave preciseandcontinuousorbitsofLEO satelliteswithan arcofarbitrarylength;and(ii) trackingmeasurements areof un-precedentedaccuracy. With GNSS,the orbital precision ofLEO satellites can now routinely reach the levelof 1 cm(see e.g.,[101]) andeventhelevel ofmillimetresover ashortperiodoftime, asdemonstratedby experimentsontheground

e.g.,[53]).Muchhigheraccuracycanbeexpectedwhenthenewgenerationsoflaserinstrumentsareoperational(seee.g.,

[4,81,96,98,107]).

Bothlinearperturbationmethodsandthenumericalintegrationmethod(evenifitwere correct)are notabletoutilize alongorbital arcofcontinuoustrackingandunprecedentedaccuracyofmodernspaceobservationtechnology.Linear per-turbationmethodsareonlyvalidinasmallneighbourhoodofmeanorbitalelementsandalmostbreakdowninthecaseof smalldivisorssuchasresonances,criticalinclinationsandcircularorbits.Inthecaseofthenumericalintegrationmethod, letustreatitasifitwere correctfornow,commonpracticeisalwaystodividealong arcintomanysmallpieces,sayin hoursoroneday,becausethemodellingerrorwillincreasewithtime.Uptoacertainepoch,themodellingerrorwill dom-inatesuchthatagravitationalsolutionwouldnotbephysicallymeaningfulanymore.Tocontroltheincreaseofmodelling errors,onewouldhavetodividealongarcintomanyshortarcs.Adirectconsequenceofthiscommonpracticeisthatwe willnotbeabletoextractsmallgravitationalsignalsfromsatellitetrackingmeasurements,sincesmallgravitationalsignals wouldtaketimetoshowuptheireffectsontheorbit.Thus,tosummarize,wewouldconcludethatbothlinearperturbation methodsandthenumericalintegrationmethod(asifitwerecorrect)aretooapproximatetobenefitfromanddonotmatch profoundtechnologicaladvancesinmodernandfuturespaceobservation.Actually,LambeckandColeman[65]pointedout thatfurtherimprovementsbothinsatellite gravitytheoryanddataevaluationmethodswererequiredbeforethenext gen-erationofsatellitegravitymissionswaslaunched.

Thepurposesofthispaperarethreefold:(i)tobriefly reviewthemethodsoflinearperturbationandtoprovethatthe numericalintegrationmethodisgroundless,mathematicallyandphysically,thoughthe numericalintegrationmethodwas essentially first published inthe mathematical literature by Gronwall [28] (see also [26,37,90]) andthen widely used in chemistryandphysics[20,39,68],statistics[84,108],andlikelyindependentlydevelopedandappliedinsatellitegravimetry

[2,3,5,67,70,76,86,90]; (ii) given differential equationswithunknown parameters andunknown initial conditions, we will construct the linearizedsolutions to the differential equations in termsof the unknown corrections to approximate val-uesoftheunknowndifferentialequation parametersandtheunknowninitialconditionsbyextendingEulerandmodified Eulernumericalintegrationmethods.Theselinearizedsolutions areoflocalnature,since theyarederived withanominal referencetrajectory.Theyaremathematicallyrigorousandrequirenoassumptionofzeroinitialvaluesforthepartial deriva-tiveswithrespecttotheunknown differentialequationparameters, asotherwise incorrectlydocumented intheliterature ofmathematics, statistics,chemistry, physicsandsatellite gravimetry.Theselocalsolutions shouldhelp betterunderstand theadvantages,disadvantagesandlimitations/problemsoflinearperturbationmethodsandthefoundationalerroneousness ofthenumericalintegration method;and(iii)to constructmathematicallyimprovedandglobaluniformlyconvergent so-lutions to thegoverning nonlineardifferential equationsof LEO satellite motionsuch that they can take fullcurrent and futuretechnological advancesof spaceobservationto extractsmallestpossible gravitationalsignalsfromsatellite tracking measurements.Asaresult,weexpecttoproduceglobalhigh-precisionhigh-resolutiongravitationalmodels,whichcanalso be calledthenext generationof globalgravitationalmodels. From themathematicalpoint ofview,the accuracyand res-olution ofthe next generation ofglobal gravitationalmodels can be sufficiently highup to the limit that modern space observationcanprovide.

Thepaperisorganizedasfollows.Section 2willbriefly reviewlinearperturbationmethods,withemphasisforthe de-terminationofgravitationalmodels.Sincethenumericalintegrationmethod,thoughfirstpublishedbyGronwall[28]almost 100yearsago(seealso[26,37,90]),hasbeenwidelyusedbyalmostallmajorinstitutionsworldwidetocomputeglobal grav-itationalmodelsfromsatellitetrackingmeasurements(ofCHAMPandGRACEmissions),wewillfirstoutlinethemethodand thenfollow[112]toprovethatassigningzeroinitialvaluestothepartialderivativesofsatellitepositionandvelocitywith respecttothegravitationalunknownparameters,namely,theharmoniccoefficients,ismathematicallyerroneousand physi-callynotpermitted.WewilldevelopEulerandmodifiedEulernumericalintegrationmethodstosolvenonlineardifferential equations withunknown parameters andunknown initial conditions in Section 3.As a result, we can representthe lin-earizedlocalsolutionstotheoriginalnonlineardifferentialequationsintermsoftheunknowncorrectionstoapproximate valuesoftheunknownparameters andtheunknown initialconditions,whichcanthen beusedtoestimate theunknown differentialequationparametersfrommeasurements.Thecontentofthissectionshouldhelpunderstandcorrect implemen-tationofnumericalintegrationtechniquesforgravitationalmodelling.InSection4,byassumingthatLEOsatelliteorbitsare precisely measured withGNSS,we willpresenta measurement-basedperturbation theory,asoriginally developedby Xu

[111],whichguarantees mathematicallyglobaluniformconvergenceofthesolutionstotheNewton’sdifferentialequations ofsatellitemotionforsatelliteorbitsofarbitrarylength.Finally,inSection5,wewillproposethemethodof measurement-based condition adjustment with unknown parameters to reconstruct global gravitationalmodels fromsatellite tracking measurements.

2. Linearperturbationandthestandard-implementednumericalintegrationmethod

2.1. Linearperturbationmethods

ThemotionofartificialsatellitesisgovernedbyNewton’slawofgravitation.Almostallearlierworksonsatellite gravime-tryandcelestialmechanicsarebasedonLagrange’splanetaryequations:

da dt =

2

na dT

de dt =

1−e2 na2 e

∂

T

∂

M−

(

1−e2

)

1 /2

na2 e

∂

T

∂ω

(1b)

d

ω

dt =−

cosi na2

(

1−e2

)

1 /2 _sini

∂

T

∂

i +

(

1−e2

)

1 /2

na2 e

∂

T

∂

e (1c)

di dt =

cosi na2

(

1−e2

)

1 /2 _sini

∂

T

∂ω

−

1

na2

(

1−e2

)

1 /2 _sini

∂

T

∂

(1d)

d

dt =

1

na2

(

1−e2

)

1 /2 _sini

∂

T

∂

i (1e)

dM dt =n−

1−e2 na2 e

∂

T

∂

e −

2

na

∂

T

∂

a (1f)

(seee.g.,[12,34,52,102]),whereK=[a, e,

ω

, i,

, M] arethesixKeplerianorbitalelements,whichstandforthesemi-major axisoftheorbitalellipse,the eccentricity,theargumentoftheperigee,theinclination oftheorbitalplane, thelongitude ofthe ascending node andthe mean anomaly,respectively; n is themean motion, andT is the disturbingpotential (of unknownforce orequationparameters p), whichisusually asmallquantity.Incelestialmechanics ofthesolarsystem, T

cancome fromdisturbingplanets ofextremely smallmasseswhen comparedwiththe solarmass(see e.g., [34]); inthe caseof artificial Earth’s satellites,T can be mainly due to the disturbing potential of the Earth, celestial bodies of very largedistancessuch astheSunandtheMoon,and/orother disturbingforcessuch asthesolidearthandoceantides,the radiationpressureoftheSunandtheairdragoftheatmosphere(see e.g.,[47,52]).Sincethegeneralrelativisticdragging, namely, Lense–Thirringeffect, has been shown to be of significant impact on satellite orbits (see e.g., [42–44,88,89]), it shouldbe fullytakenintoaccount inthefuturecomputationofhighprecision, highresolution gravitationalmodelsfrom trackingmeasurements.One mayequivalentlyrewrite Lagrange’splanetaryEqs.(1)insixother orbitalelements(seee.g.,

[12,96,102]).

Withoutlossofgenerality,wedenotethegeneralsolutiontoLagrange’splanetaryEqs.(1)byK(t,p,ck),whereckstands

forsixarbitraryintegrationconstants.Differentintegrationconstantsspecifythemotionsofdifferentsatellites.Inprinciple, ifthegeneralsolutionK(t,p,c_k)wouldbeanalyticallyavailable, giventheparameters pandthesixintegrationconstants

c_k0 (oralternatively an initial pointKt0 or anysixindependentvalues onK(t, p,ck)),one canthen obtain theparticular

solutionK(t, p, ck0 ) and useit tocompute andpredict theorbitof the celestialbody atany time t.On the other hand,

giventhegeneralsolution K(t,p,ck) anda sufficientnumber ofmeasurements onK(t, p,ck),one can thenestimate the

(unknown)forceparameterspfromthemeasurements.Unfortunately,thegeneralanalyticalsolutionK(t,p,c_k)canonlybe obtainedfortheidealizedtwo-bodyprobleminwhichaparticleofnegligiblemassisattractedbyanotherpointmass(see e.g.,[12,52,83,102]).Ingeneral,noanalyticalsolutionto(1)canbepossible,fortworeasons:(i)Lagrange’splanetaryEqs.(1) arenonlinearandextremelydifficulttosolveanalytically.Thus,perturbationtheoryhasbeenplayingafundamentalrolein celestialmechanicsandsatellite gravimetry(seee.g.,[11,14,34,52,60,102,123]),whichattemptstoconstructanapproximate solutiontoLagrange’splanetaryEqs.(1)throughthe procedureofsuccessiveapproximation.Actually,therearetwo types ofperturbationmethods.Oneistoconstructanapproximatesolutionthroughthemathematicalstandardapproachofsmall parameterperturbation(seee.g.,[11,14,34,79,102,123]).Theothermethodistotreatthevariablesontherighthandsideof Lagrange’splanetaryEqs. (1)asconstants,exceptforthe meananomaly M,andthen integrate thedifferentialEqs. (1)to constructanapproximatesolution.Thislatterapproachhasbeenwidelyapplied insatellitegeodesy(see e.g.,[49,52,60]); and(ii)thedisturbingpotentialfunctionTitselfmaynotbeexactand/orsufficientlyprecise.Actually,perturbationtheory, togetherwithastronomicalmeasurements,historicallyplayedadecisiveroleincorrectlyidentifyinganunknowndisturbing celestialbodyto explainthe deviationsofthetheoreticalpredictions frommeasurements by Adams(https://en.wikipedia. org/wiki/John_Couch_Adams), Le Verrier (https://en.wikipedia.org/wiki/Urbain_Le_Verrier) andLowell [71] (see also https: //en.wikipedia.org/wiki/Planets_beyond_Neptune),successfullyleadingtothegreatdiscovery ofbothNeptuneandPlutoin 1846and1930,respectively(see e.g.,[29,34,66]), thoughthedata ofNeptunegivenbyLe VerrierandAdamsare inlarge errors(seee.g.,[38]).

Asan inverseproblemofcelestialmechanics,satellite gravimetryis toreconstructtheunknownforce parametersp in thenonlinearLagrange’splanetaryEqs.(1)withunknowninitialconditionsfromasufficientnumberofmeasurements on

K(t,p,c_k).Inthiscase,weassumethatthedisturbingpotentialfunctionTitselfispreciselygivenbutcancontainanumber

ofunknownparameters p,thoughpartofTmaybedirectlymeasuredandcorrected.Wealsoimplicitlyassume thatthere existnootherunknownsourcesthatcancontributetoTinanon-negligibleway.Insatellitegeodesy,thedisturbingpotential

TismainlyattributedtotherotatingEarth,theattractionbytheSun,theMoonandlargeplanets,thesolidEarthandocean

thedisturbingpotentialToftheEarthinthenon-inertialearth-fixedreferenceframeasfollows:

T=GM r

Nmax

l=2

l

m=0

_R

r

l

[Clmcos

(

m

λ

)

+Slmsin

(

m

λ

)

]Plm

(

cos

θ

)

, (2)

(seee.g.,[35,52,124]),whereNmax isamaximumnumberofdegreesandorders,RisthemeanradiusoftheEarth,Clm and

Slm are theunknown normalized,dimensionlessharmoniccoefficientswhich willbecollectedintothe unknownvector p

andtobeestimatedfromsatellitetrackingmeasurements,

λ

and

θ

arethelongitudeandcolatitude,respectively,andPlm(t)

isthenormalizedLegendrefunction.However,insatellite gravimetryandcelestialmechanics, thedifferential equationsof motionofasatellitearealmostalwaysgivenintheinertialreferenceframe(seee.g.,[52,96,102]).Inthiscase,wewillneed thetransformationofcoordinates fromtheearth-fixed referenceframeintotheinertialreferenceframethrough rotations (seee.g.,[96]).Morespecifically,ifweusesphericalcoordinatesystems,thenweneedthefollowingtransformation:

λ

=

α

+

δα

−

ω

et, (3a)

θ

=

ζ

+

δζ

, (3b)

(seee.g.,[47]),where

δα

and

δζ

arethecorrectionsto

α

and

ζ

,whichdependonprecession,nutation,theEarth’srotation andpolarmotionandcanbecomputedfromtheoryandmeasurements(seee.g.,[96]),

α

and

ζ

aretherightascensionand co-declinationintheinertialreferenceframeofepochJ2000.0,

ω

eistherateoftheEarth’srotation.Thedisturbingpotential

Tintheinertialreferenceframeisclearlyafunctionoftime.

TorigorouslydeterminemathematicallytheunknownforceparametersClmandSlmfromsatellitetrackingmeasurements,

wehavetofirstexactlysolveLagrange’splanetaryEqs.(1),withthedisturbingpotential Tgivenby (2),linktheexact so-lutiontothesatellitetrackingmeasurementsandfinallyestimateClm andSlm.Unfortunately,thenonlineardifferentialEqs.

(1)aretoocomplicatedtoexactlysolveanalytically. Thus,perturbationtheoryisalwaysusedtoconstructanapproximate solutionto(1)(seee.g.,[11,34,49,52,60,102,124]).Themostcompleteperturbationtheorywasfullydevelopedforthe deter-mination ofClm andSlm fromsatellite trackingmeasurements by Kaula[49,52],following theapproachof Kozai[60] and

giventherepresentationofTexpressedintermsofthesixKeplerianorbital elementsbyGroves[124].Itisnowadayswell

knownasKaulalinear perturbationtheory.Morespecifically,sincetheexact analyticalsolutionisgenerallyhardoralmost impossibletoobtain,Kaula[49,52]derived thelinearperturbationsolutiontoLagrange’splanetaryEqs.(1)bytreating all the orbitalelementsasconstants and/orby replacing themwiththe meanorbital elements, exceptforthe rapidly time-varyingelementofthemeananomalyM,andthenintegratingallthetermsontherighthandsideof(1).

The major advantage ofKaula’s linear perturbation solutionis its suitability to analyze the physicalproperties ofthe solution.Since T isexpressed interms ofsixKeplerian orbital elements,the physicalfeatures of Tcan be further classi-fiedinto secular,long-periodic andshort-periodic(see e.g.,[49,52,60,124]). Moreprecisely, thetermsasa functionofthe meananomalyMwillchangeperiodically(andrapidly)andareshort-periodic;thetermsasafunction of

ω

butnot Mare

long-periodic; andthetermsirrelevant of

ω

norMare secular.Inother words,secular termschanges slowlybut approxi-matelylinearlywithtime. Inaddition,one canalsoidentify,in thelinearperturbationsolution,the physicallyinteresting phenomenon of orbital mean-motion resonancewhen the rotation rateof the Earthand the meanmotion of a satellite are commensurable(see e.g.,[9,16,34,52,57,102,119]). Nevertheless,Kaula’slinearperturbationsolutionisonlyvalidlocally around the neighbourhood ofthe meanorbital elements andwilldiverge withthe increase oftime. Thus, Kaula’slinear perturbationsolutionwillnotbe abletofullyutilizealongorbital arcandhighprecision ofmodern spaceobservationto estimatep.

2.2. Thestandard-implementednumericalintegrationmethod

IntheCartesiancoordinatesystem,themotionofanartificialsatellite canalsobemathematicallywrittenalternatively bythefollowingnonlinearvectordifferentialequation:

¨

x=aE

(

t,x,x˙,p

)

+aM

(

t,x,x˙,pM

)

, (4)

(see e.g.,[12,52,96,102]),wherex istheposition vectorof thesatellite inthe inertialreferenceframe,aE

(

t, x, x˙, p

)

is the

Earth’sgravitationalattraction exertedon thesatellite, p isthe vector ofunknown parameters,andaM

(

t, x, x˙, pM

)

stands

forallother forceswhichmayincludesolidearthandoceantides,atmosphericaldrag, solarradiationpressure and third-body effects (see e.g., [47,86,96,102]), with p_M standingfor the unknown parameters (if any) of theseforce models. For the next generationofhighprecision, highresolution globalgravitationalmodels fromsatellite tracking,theforce model

aM

(

t, x, x˙, pM

)

shouldincludethegeneralrelativisticdragging,namely,theLense–Thirringeffect,whichhasbeenshownto

haveasignificanteffectonsatelliteorbits,satellite-to-satelliterangesandrange-ratesforasufficientlylengthyarc(seee.g.,

[42–44,88,89])buthasnotbeenconsidereduptothepresentintheproductionofglobalgravitationalmodels.

Intheremainder ofthispaper,wewillusex¨ andx˙ tostandforthefirstandsecondderivativesofxwithtime, respec-tively.(Weusethenotationdx/dttostandfortimederivativein(1),sincethedotofitheredoesnotlookgood)Sincethe accelerationaE

(

t, x, x˙, p

)

isindependentofx˙, itcanberewrittenasfollows:

aE

(

t,x,p

)

=− GM

r3 x+

∂

T

whereGMistheproductoftheEarth’smassMandthegravitationalconstantG,r=

x

, andTisthedisturbingpotentialof theEarth’sgravitationalfieldwiththeparametersofharmoniccoefficients.Because(4)isformulatedinaninertialreference system(seee.g.,[47,102]),theearth-fixed

λ

and

θ

inthedisturbingpotentialTof(2)mustfirstbetransformedthrough(3) into theinertial reference frame (see e.g., [47]). Since the second termaM

(

t, x, x˙, pM

)

of (4) addsno new mathematical

difficultytosolvethenonlineardifferentialEq.(4),wewilllimitourselvestotheEarth’sgravitationalfieldintheremainder ofthispaper.Thus,thenonlineardifferentialEq.(4)canbesimplifiedasfollows:

¨

x=aE

(

t,x,p

)

, (6)

wherepisanunknownvectorofequationparameters.Initialconditionsto(6)areunknownaswell.

We should note, however,that ifaM

(

t, x, x˙, pM

)

exists but isneither estimated together withp norcorrected with a

sufficientlyprecisemodel,theeffectofaM

(

t, x, x˙, pM

)

willbeabsorbedintotheestimateofp.Suchaneffectistheoretically

systematicand of signal nature. It cannot be filtered out, since, in principle, a filter can only reduce the level of noise butwould mess up orsmear the truesignal ofconcern, ifthe signals are not constant inside the window ofthe filter. Inotherwords,filtering,ascurrentlyusedinsatellitegravimetry,willdefinitelydistortthetruesignalofinterest,because gravity signalsare clearly not constant inside the window offiltering. Totest the Lense–Thirring dragging fromsatellite measurements,forexample,oneshouldprobablyhavetoeitherestimatetheeffecttogetherwiththegravitationalmodelor todevelopstatisticalhypothesistestingmethodsforweakcontinuousfunctions;thelatterisoftheoreticalinterestbyitself anddeservesaseparatefullresearch.

Givenanumberof(geometrical)trackingmeasurementstothesatellite suchaspositionsofthesatellite,ranges,range ratesanddirectionsto thesatellite, denoted byy₁ , y₂ , . . . , yn or inthevector formy,theproblemof satellite gravimetry

istousethetrackingmeasurements ytodeterminethegravitationalparametersp.Mathematically,thisisessentially the problemofestimatingtheunknownparameterspofthenonlineardifferentialEq.(6)withunknowninitialconditionsfrom satellitetrackingmeasurements. IfthereareanumberofLEO satellites,sayssatellites,themotionofeach satellitebeing governedbythesamedifferentialequationsoftype (6)withthesamegravitationalparameterspbutwithdifferentinitial conditionsordifferentintegrationconstants.Ifwecollectsatellitetrackingmeasurementsyiontheithsatellite,thenwewill

havetocombineallthesemeasurementsy₁ , y₂ , . . . , ystogethertosolvefortheparametersp.Inthefollowingdevelopment

ofthemethod,withoutlossofgenerality,wewillconfineourselvestoonesatellite.

Todeterminepfromy,oneofthemostimportantstepsistorepresenteachyiintermsofp.Sinceageometricalsatellite

trackingmeasurementyiisgenerallya functionoftheposition andvelocity ofthesatellite attheithepoch,inprinciple,

wehavetofirstsolvethenonlineardifferentialEq.(6).Letx(t,p,c)denotethegeneralsolutiontothenonlineardifferential

Eq.(6), withc standingforsixarbitraryintegrationconstants. Thesesixintegrationconstantsc aremathematically inde-pendentoftheequationparametersp.Inother words,thegeneralsolutionx(t,p,c)mathematicallyrepresentsaninfinite numberofsolutions tothe differentialEq.(6),which canphysicallydescribe themotions ofdifferentsatellites.As faras initialvaluesareproperlygiven,one canthenusethemtofixsixarbitraryintegrationconstantscandobtain thespecific solutiontouniquelydescribethemotionofthesatellite.Mathematically,thevectorcforthisparticularsolutioncannowbe expressedasthefunctionsofp andtheinitialconditions.Insatellitegeodesy,initialvaluesareoftenthreeinitialposition

coordinatesx₀ andthreeinitialvelocitycomponentsv₀ ofthesatelliteattheinitialepocht₀ .Asaresult,theorbitalposition solutionofmotionofthesatellite canbe implicitlywrittenasx(t,p,c(p,x₀ ,v₀ )).Weemphasizethatsixarbitrary integra-tionconstantsccanalsobealternativelydeterminedfromanysixindependentvaluesonthesolution,insteadoftheinitial positionandvelocityconditions,sinceanyspecificvalueontheparticularsolutioncontainstheinformationonc.Exceptfor theidealizedtwo-bodyproblemwiththepointmassmodel,itisalmostimpossibletoobtainananalyticalsolutiontothe nonlineardifferentialEq.(6)withtheinitial conditionsx₀ andv₀ .Thus,wecanonlyusetheimplicitorbitx(ti,p,c(p,x0 , v₀ ))todevelopobservationalequationsforgeometricaltrackingmeasurementsofanykind.Takinga(velocity-independent) geometricaltrackingmeasurementyiasanexample,wecansymbolicallywriteitsobservationalequationasfollows:

yi=f

(

x

(

ti,p,c

(

p,x0 ,v0

)))

+

ǫ

i, (7)

wheref(·)standsforanonlinearfunctionaland

ǫ

iistherandomerrorofthemeasurementyiwithzeromean.

Undertheframeworkofthenumericalintegrationmethod,asdescribedinLerchetal.[67],Longetal.[70]andReigber

[86]forsatellitegravimetry,wewillhavetofirstlinearizetheobservational Eq.(7).Givena setofapproximatevaluesp0 ,

x0 ₀ andv0 ₀ forp,x₀ andv₀ ,respectively,onecanthennumericallyintegratethenonlineardifferentialEq.(6)andobtainthe approximatepositionofthesatelliteattimeepochti,whichisdenotedbyx0 i.Thus,theobservationalEq.(7)canbeformally

linearizedasfollows:

δ

yi=aix

x0 +aiv

v0 +aip

p+

ǫ

i (8a)

attheapproximatevaluesofp0 ,x0 ₀ andv0 ₀ , where

δ

yi=yi−f

(

x0 i

)

,

x0 =x0 −x0 ₀ ,

and

p=p−p0 .

Thethreerowvectorsa_ix,a_ivandaipare allcomputedattheapproximatevaluesp0 ,x0 0 andv0 0 , anddefined,respectively,

asfollows:

aix=

∂

f

(

x

(

ti,p,c

(

p,x0 ,v₀

)))

∂

xT

0

=

∂

f

(

x

(

ti,p,c

(

p,x0 ,v0

)))

∂

xT

∂

x

∂

xT

0

, (8b)

aiv=

∂

f

(

x

(

ti,p,c

(

p,x0 ,v0

)))

∂

vT

0

=

∂

f

(

x

(

ti,p,c

(

p,x0 ,v0

)))

∂

xT

∂

x

∂

vT

0

, (8c)

aip=

∂

f

(

x

(

ti,p,c

(

p,x0 ,v0

)))

∂

pT

=

∂

f

(

x

(

ti,p,c

(

p,x0 ,v0

)))

∂

xT

∂

x

∂

pT. (8d)

Asakey steptoestimatepfromy,we havetocompute thevectorsaix,aivandaip.Thecommonmatrixofthepartial

derivatives

∂

f(·)/

∂x

T_{in(8b)–(8d)canbereadilyobtained,aswidelyavailable(seee.g.,[49,52,67,70,105]).Ifwewouldtreat}

p as ifit were givenfornow, theproblem isturned into a standard problemof statisticalorbitdetermination withthe nonlineardifferentialEq.(6),theinitialconditionsx₀ andv₀ ,andthemeasurementsy.Inthiscase,computingthematrices ofthepartialderivatives

∂

x/

∂

xT

0 and

∂

x/

∂

vT0 istheoretically equivalenttofinding thestatetransitionmatrixforthestate

ofpositionandvelocityfromtheinitialepocht₀ tothecurrentepochti.Thisproblemhasbeencompletelysolvedandwell

documentedin,forexample,[32,70,103,105].

Now,tocompletethefinalconstructionoftheobservationalEq.(8a),thekeyissueistocomputethepartialderivativesof

x(ti,p,c(p,x0 ,v0 ))withrespecttopin(8d),namely,

∂

x/

∂

pT₌

_∂

_x

₍

_t

i, p, c

(

p, x0 , v0

))

/

∂p

T _{forconcisenessofnotations.Since}

wedonothaveananalyticalsolutionx(ti,p,c(p,x0 ,v0 )),itisnotpossibletodirectlycomputeits partialderivativeswith

respecttotheparametersp.Instead,onehasattemptedtoobtainthesepartialderivativesthroughsolvingtheirdifferential equations.

Tostartwith,letuscollectthesatellitepositionxandvelocityvattimetinthevectorzanddenotethepartial deriva-tivesofzwithrespecttopbyS(t,p),namely,

z=

(

xT_,vT

₎

T_{, (9a)}

S

(

t,p

)

=

∂

z

∂

pT. (9b)

Thepartialderivatives

∂

x/

∂

pT_{isobviouslypartofamoregeneralmatrixS(t,p)ofpartialderivatives.Ithasbeenrigorously}

shownmathematicallythatS(t,p)canbedirectlyderivedfromtheoriginalnonlineardifferentialEq.(6)(see,e.g.,[3,76,90]) andisgovernedbythefollowingsystemofdifferentialequations:

∂

S

(

t,p

)

∂

t =

0 I

∂

aE

(

t,x,p

)

/∂xT 0

S

(

t

)

+

0

∂

aE

(

t,x,p

)

/∂pT

. (10)

Obviously,the derived Eq.(10)ofthe partial derivativesdonot mathematicallyadd anynew informationonthe original problemofsatellitegravimetry.Actually,forageneraldifferentialequationwithunknownparameters,thedifferential equa-tionsoftype(10)werealreadygivenbyGronwall[28]andRitt[91]almost100yearsago(seealso[26]).

Although we have the differential Eq.(10) for the partial derivatives S(t, p), they can still be useless, unlessthe ini-tial conditionsofS(t,p) at thetime epoch t₀ are available. Weuse theitalic font for“the” before “initial conditions” to emphasizethattheinitialconditionsofS(t,p) cannotbearbitrarilygivenbutmustcomplywiththeoriginalproblem. Un-fortunately,theoriginalproblemofsatellitegravimetry,namely,thegoverningdifferentialEq.(6)andthesatellitetracking measurementsy,donotprovideanydirecthint/clueonwhatvaluesS(t0 ,p)cantakeon.

mathematicallyincorrectandphysicallynotpermitted.Itisunbelievabletoseethatthenumericalintegrationmethodends uponmathematicallyvaingroundwithoutrhymeorreason,thoughithasbeenwidelyusedbyalmostallmajorinstitutions worldwideto produceglobalgravitationalmodels fromsatellite trackingmeasurements, withwidestpossible applications inEarthSciences.

2.3.NozeroinitialvaluesforS(t0,p)permittedmathematicallyandphysically

Inthispartofthepaper,wewillprovethatnozeroinitialvaluesforS(t₀ ,p)canbepermittedfromboththe mathemat-icalandphysicalpointsofview.WewilluseacounterexamplereportedinXu[112,115]for this purpose,though one can readilyconstructmanyothercounterexamples.Thenwe willusestrictlylogicalreasoningstoexplain whyS(t₀ ,p)cannot bezero.Forsomeofthearguments,thereaderisreferredto[112]fordetails.

Asageneralruleinmathematics,weneednothingmorethanacounterexampletodisprovethatsomethingisincorrect. Letusstartwiththesecondcounterexamplein[112],whichisrewrittenasfollows:

¨

y+p2 ₁ y−p2 cos

(

p1 t

)

=0, (11)

wherep₁ andp₂ aretwoequation(unknown)parameters.BydirectlysolvingthedifferentialEq.(11),weobtainthegeneral

solution:

y

(

t

)

= p2

2p2 ₁

{

cos

(

p1 t

)

+p1 tsin

(

p1 t

)

}

+c1 sin

(

p1 t

)

+ c2

p1 cos

(

p1 t

)

, (12)

wherec₁ andc₂ are two arbitraryintegrationconstants. Mathematically,integration constantsc₁ andc₂ are independent ofp₁ andp₂ .As farasc₁ andc₂ aregivenspecificvalues,whichcanbeimplicitlydefined, forexample,throughassuming twovaluesofy(t) attwo differenttimeepochs, wewillthen obtaintheparticularsolutionof(12).Withthetruegeneral solution(12),wecaneasily computeandobtainthetruevaluesofthederivativesofy(t)withrespecttop₁ andp₂ ,which aresimplygivenasfollows:

∂

y

(

t

)

∂

p1 =−

p2

p3 ₁

{

cos

(

p1 t

)

+p1 tsin

(

p1 t

)

}

+ p2

2p1 t 2 _cos

₍

_p

1 t

)

+c₁ tcos

(

p₁ t

)

− c2

p2 ₁ cos

(

p1 t

)

− c2 t

p₁ sin

(

p1 t

)

, (13a)

∂

y

(

t

)

∂

p2 =

1

2p2 ₁

{

cos

(

p1 t

)

+p1 tsin

(

p1 t

)

}

. (13b)

Foranarbitraryt₀ ,thederivatives

∂

y(t)/

∂

p₁ of(13a)and

∂

y(t)/

∂

p₂ of(13b)clearlycannotbezero.Actually,ifthe deriva-tives(13a)and(13b)wouldbeequaltozeroatthetimeepocht₀ ,wewouldreadilyhavetwoequationsfortwounknowns

p1 andp2 andwouldbeabletosolvethemwithoutanyvalueofy(t).Forexample,iftheinitialvaluesforthesederivatives

wouldbeallowedtobeequaltozero,thenthesecondderivative(13b)wouldbeturnedintothefollowingequation:

cos

(

p1 t0

)

+p1 t0 sin

(

p1 t0

)

=0, (14) foranynon-zerop₁ .Bysolvingthisequation,wecouldobtainthevalue(s)ofp₁ (ifsolutionsexist).Obviously,thisislogically ridiculous,since thiswould indicate that we would be able to determine the unknown parameter p₁ in the differential

Eq.(11) withoutanyinformationony(t). Inthecaseoft0 =0, (14)becomes1=0— anevenmoreridiculousexpression.

Thesourceoferrorsisclearlywiththeassumptionofsettingthevaluesofpartialderivatives(13a)and(13b)tozeroatthe initialepocht₀ .

Forthesatellite gravimetry problem(6)withthe geometrical satellite trackingmeasurements y, we willsimply make some logicalreasonings andexplanations. The mathematicalproof andphysical explanations ofno zeroinitial valuesfor

S(t₀ ,p)canbefoundin[112].

Remark1:Forthegeneralorbitsolutionx(t,p,c)of(6)withoutgiveninitialconditions,theproblemofsatellite

gravime-tryistodetermineboth pandcfromthegeometricalsatellite trackingmeasurements y.The orbitalpositionatthetime epochticanbemathematicallywrittensymbolicallyasx(ti,p,c).x(ti,p,c)isonlyapointofthegeneralsolutionx(t,p,c)

ofthesatelliteatthetimeepochti.Logically,noorbitalpositionatonetimeepochissuperiortoanyotherorbitalpositions

ofthesameorbitbutatdifferenttimeepochs.Sincet₀ isarbitrary,anyorbitalpositioncanequallyserveasaninitial con-dition.IfS(t₀ ,p)couldbe settozero,thenalltheother S(ti,p)(ti=t0 ) couldbetreatedinthesamemanneraszerofrom

themathematicalpointofview,implyingphysicallythatx(t,p,c)wouldnotbeafunctionofp.Obviously,thisconclusion violatesourstartingdifferentialEq.(6)withp.Thus,S(t0 ,p)cannotbeequaltozero.

Remark2: According to[115,116],forthe particularorbitsolution x(t, p,cx0 ) of(6),withtheintegration constantscx0

fixed/given,forexample,fromtwoorbitalpositionsattwodifferentepochs,thepartialderivativesS(t,p)atanytimeepoch

arbitrarilydifferentfromeithert₀₁ ort₀ .BytheclaimofRileyetal.[90](seealso[67,70]),wehaveS

(

t_{0 i}, p

)

=0, withwhich

wecanfurtherobtainS(t,p,t₀ i)bysolvingthedifferentialEq.(10).Forthreearbitrarilydifferenttimeepochst0 ,t01 andt0 i,

their correspondingpartial derivativesatthesametime epocht, namely,S(t,p,t₀ ), S(t,p,t₀₁ ) andS(t,p,t₀ i),willnot be

equalto eachother.Thisobviouslycontradictsthefact thatS(t,p)isunique forthisparticularorbit. Thesource oferrors againcertainlycomesfromtheincorrectassumptionofzeroinitialvaluesforS(t,p).

Remark3:IfS(t,p)couldbesettozeroattheinitialepocht₀ ,byfollowingthesamelogicalreasoningsasin(13)and/or

(14),wewouldbeabletosolveforpfromthesystemofequationsS

(

t0 , p

)

=0, sincethenumberofequationsS

(

t0 , p

)

=0is exactlyequaltothatoftheunknownparametersp,implyingthatwecoulddeterminetheunknownharmoniccoefficientsp

withoutanysatellitetrackingmeasurements;thisisagainanunacceptableresult.Actually,ontheotherhand,ifS

(

t₀ , p

)

=0,

thenwe couldsolvethedifferentialequations(10)andobtaintheorbitalsolution,denotedbyxt0 (t,p,c(p,x0 ,v0 )),which

impliesthat we do not needany initial conditionsx₀ and v₀ tofind the particularorbital solution.Since t₀ is arbitrary, wecouldobtainan infinitenumberofdifferentsolutionsforthesamesatellitegravimetryproblem. Allthesearecertainly incorrectmathematically,againwiththesourceoferrorsintheassumptionofS

(

t₀ , p

)

=0.

Theabovecounterexample,togetherwithallthemathematical,physicalandlogicalreasonings,hasall clearlynullified theclaimbyRileyetal.[90]thattheinitialvaluesofthepartialderivativeswithrespecttoequationparametersare gener-allyzero.Actually,thisclaimisalsousedasastartingpointfortheGoddardSpaceFlightCentersoftwaresystemGEODYN (seee.g.,[67,70])andinEurope(seee.g.,[86])fortheproductionofglobalgravitationalmodelsfromsatellitetracking mea-surements. Bearingin mind that a variety ofglobal gravitationalmodelproducts hasbeen widely used inand farmore beyondgeodesysuch assolidgeophysics,hydrology,continentalwatervariation andso on,we believethat software sys-temsmustbefirstupdatedontoasolidmathematicalfoundationrightnowbeforecontinuingtoproduceandcirculatesuch gravitationalproductsfromsatellitetrackingmeasurements.Finally,westateatheoremin[112]toconcludethissubsection asfollows:

Theorem1. Giventhegoverningvectordifferentialequations(6)withtheunknownharmoniccoefficients pin(2),thensetting theinitialvaluesofthepartialderivativesoftheorbitandvelocitywithrespecttotheunknownharmonic coefficientsptozero atanyspecifiedinitialtimeepocht0 isnotpermitted,mathematicallyandphysically.

Beforeclosingthissection,Ishouldpointoutthatzeroinitialpartialderivativeswithrespecttotheunknownparameters ofdifferentialequationshasbeenroutinelyused beyondgeodesy.Completely independentofthedevelopmentinsatellite geodesy,forsimplicitybutwithoutlossofgenerality,givenanordinarydifferentialequationy˙=f

(

t, y, p

)

, Gronwall[28] cor-rectlyderived thedifferentialequation ofywithrespecttotheunknown equationparameterp,asinthecaseof(10)(see also[91]),butincorrectly claimedthat its initial value isequalto zerowithout providinganyreasons orarguments.The work of Gronwall [28] wasthen further spread through the book on ordinary differential equationsby Goddington and Levinson[26].Actually,thesolutiontothegivenordinarydifferentialequationcanbesymbolicallywrittenasfollows:

y

(

t,p

)

=y

(

t0 ,p

)

+ t

t0

f

(

t,y,p

)

dt, (15)

fromwhichwecanonlyobtainthefollowingidentity:

dy

(

t,p

)

dp

_t

=t0

= dy

(

t0 ,p

)

dp , (16)

butcertainly not thezero initial derivative,aswe have provedin thispaper.Since all thegeodetic literature on satellite gravimetryhasnotcitedormentionedanyofthesemathematicalpublications,itseemsthattheclaimofzeroinitialpartial derivativeswithrespectto theequation parameters hasbeentakenforgranted everywhereforalmost 100 years,though incorrectly,aswehaveprovedinthispaperandXu[112].Theincorrectclaimofzeroinitialderivativesnowstillcontinues tospread,ascanbeseen,forexample,in[20,37,39,68,84,108].

3. Linearizationandnumericalintegrationtechniquesforestimationofunknowndifferentialequationparameters frommeasurements

IfS(t0 ,p)=0andisunknown,thenthedifferentialEq.(10)arenotuseful.Asaconsequence,wearenotabletocompute

aipof(8d)tocompletetheconstructionoftheobservationalEq.(8a).Thequestionnowishowwecanproperlyimplement

numericalintegrationtechniquestodeterminetheEarth’sgravitationalfieldfromsatellitetrackingmeasurements.In prin-ciple,giventhesatellitetrackingmeasurementsywithacorrespondingweightingmatrixW,wecanwritetheleastsquares objectivefunctionasfollows:

min:[y−f

(

x

(

ty,p,c

))

]TW[y−f

(

x

(

ty,p,c

))

] (17)

subjectto theequalityconstraintdefinedbythedifferentialEq.(6),wheref(·)are thetheoreticalvaluesofthe

measure-mentsy,andeachx(tyi,p,c) satisfies(6)andstandsforthetheoretical orbitalposition ofthe satelliteatthetime epoch

t_yi whenthe tracking measurement y_i iscollected. If initial conditionsare available, the constantsc can be alternatively expressedintermsoftheseinitial conditions.Ifthetrackingmeasurement yiisinvolvedwithtwosatellites,we needthe

theoreticalorbitalpositions(andlikelyalsothevelocities)ofthesetwosatellitestocomputethetheoreticalvalueofyi.For

moredetailsonthistypeofobservationalequations,thereaderisreferredtoXu[111].

Sincetheequalityconstraintsare givenintheformofdifferentialequations,we cannotuse conventionaloptimization methods to solve the minimization problem(17), subjectto (6). We have to use numerical techniques to discretize the differentialEq.(6)such thatwecan representx(tyi,p,c))intermsoftheunknown differentialequationparameters p.As

inthecaseof(8a),givensomeapproximatevaluesp0 ,x0 ₀ andv0 ₀ ofp,x₀ andv₀ ,respectively,wecanobtainthenumerical solutionx0

(

t, p0 , c

(

p0 , x0 ₀ , v0 ₀

))

(orsimplyx0 (t)forconcisenessofnotations)bynumericallysolvingthefollowingnonlinear

differentialequations:

¨

x0 =aE

(

t,x0 ,p0

)

(18)

underthegiveninitialvaluesofx₀ 0 andv0 ₀ .Accordingly,thesolutionofvisdenotedbyv0

(

t

)(

=x˙0

(

t

))

.

SincethedifferentialEq.(6)arenonlinear,we mayattempttofindtheirapproximatesolutionsintermsofp byeither directlylinearizing(6)orusing numericalintegrationmethods. Forconvenience,we rewrite thesecond order differential

Eq.(6)asanequivalentsystemoffirstorderdifferentialequations:

˙

z

(

t

)

= x˙_˙

(

t

)

v

(

t

)

= v

(

t

)

aE

(

t,x,p

)

. (19)

Asinthecaseof(6),boththeequationparameterspandinitialconditionsto(19)areunknown.

3.1. Thelinearizedlocalsolution

Subtractingz0 (t)from(19),wehave

˙

z

(

t

)

−z˙0

(

t

)

= x˙

(

t

)

−x˙ 0

₍

_t

₎

˙

v

(

t

)

−v˙0

(

t

)

= v

(

t

)

−v

0

₍

_t

₎

aE

(

t,x,p

)

−aE

(

t,x0 ,p0

)

. (20)

Denoting

z˙

(

t

)

=z˙

(

t

)

−z˙0

(

t

)

,

x

(

t

)

=x

(

t

)

−x0

(

t

)

,

v

(

t

)

=v

(

t

)

−v0

(

t

)

,

andthenlinearizingtherighthandsideof(20),wehave

z˙

(

t

)

=

v

(

t

)

Fax

(

t

)

x

(

t

)

+Fap

(

t

)

p

= 0 I

Fax

(

t

)

0

x

(

t

)

v

(

t

)

+ 0

Fap

(

t

)

p

= 0 I

Fax

(

t

)

0

z

(

t

)

+ 0

Fap

(

t

)

p, (21a)

Fax

(

t

)

=

∂

aE

(

t,x,p

)

∂

xT

_x

=x0_(t),p=p0

, (21b)

Fap

(

t

)

=

∂

aE

(

t,x,p

)

∂

pT

_x

=x0(t),p=p0

, (21c)

andIisa(3×3)identitymatrix.

Giventhe initialconditions x₀ andv₀ fortheoriginal problemof satellitegravimetry, wecan havethe corresponding initial conditions

z₀ for

z(t). Thus, accordingto Stengel [99] andGrewal andAndrews[27],we can readilywrite the solutiontothelineardifferentialEq.(21a)asfollows:

z

(

t

)

=

(

t,t0

)

z0 + t

t0

(

t,

τ

)

0

Fap

(

τ

)

d

τ

p, (22)

orequivalently,

z

(

t

)

=z0

(

t

)

+

(

t,t0

)

z0 + t

t0

(

t,

τ

)

0

Fap

(

τ

)

d

τ

p, (23a)

where

(t,t₀ )isthestatetransitionmatrixandisequalto

(

t,t0

)

=

(

t

)

−1

₍

_t

0

)

, (23b)

(seee.g.,[27]),andthefundamentalmatrix

(t)isthesolutiontothefollowingmatrixdifferentialequations:

˙

(

t

)

= 0 I

Fax

(

t

)

0

(

t

)

, (23c)

undertheinitialmatrixconditions

(

t0

)

=I₆ , (23d)

withI₆ isa(6×6)identitymatrix.Formorepropertiesabout

(t),includingitsuniquenessandnon-singularity,thereader

isreferredto[27].

It is clearthat the solution (23a)of the satellite orbitandvelocity isa linear vector function ofthe corrections

z₀

tothe approximateinitialvalues[x0 ₀ , v0 ₀ ]andthe corrections

pto theapproximatevaluesp0 .Therefore,we can readily linearizetheoriginalsatellitetrackingmeasurement(7)withrespectto

z₀ and

p.

3.2. Numericalintegrationmethodstoconstructlocalsolutions

When numerical integrationmethods are required, one always assumes that the functionsto be integratedare given and/or known,andthe targetis tousesuch methods tonumericallycompute the integrationofthe functions,ascan be foundinanystandardtextbooksonnumericalanalysisandnumericalintegration(seee.g.,[100,106]).Thesetechniquescan be directly used to compute an approximatereference orbitof a satellite, given initial conditions [x0 ₀ , v0 ₀ ] and p0 . How-ever,insatellite gravimetry fromsatellite trackingmeasurements, sinceboth initialconditions[x₀ ,v₀ ]andthedifferential equationparameters pareunknown,itisimpossibletoexactlycomputesatellite orbitsby directlyimplementinganywell documentednumericalintegrationmethods.

Inthispartofthepaper,unlikestandardtextbooksonnumericalintegrationtocomputetheintegralofagivenfunction (without anyunknown parameters)(see e.g., [100,106]),our basicidea isto constructa solutionto nonlineardifferential equationswithunknownequationparametersandunknowninitialconditions,withtheaidofnumericalintegration meth-ods.Moreprecisely, inthe caseofsatellitegravimetry fromtrackingmeasurements, wewillrepresenttheorbitalsolution

x(t,p,c(p,x₀ ,v₀ ))totheNewton’snonlineardifferentialEq.(6)intermsofitsapproximatevalue,theunknowncorrections

z₀ totheapproximateinitialvalues[x0 ₀ , v0 ₀ ]andtheunknowncorrections

poftheharmoniccoefficientsbynumerically solving thenonlinear differentialEq.(6) undertheinitial (unknown) conditionsx₀ andv₀ . Recall thatforeach measure-mentyi atthetimeepochtyi,weobtainthenominalapproximateorbitx0

(

tyi, p0 , c

(

p0 , x0 0 , v0 0

))

bynumericallysolving the

nonlineardifferentialEq.(18)undertheinitialconditionsx0 ₀ andv0 ₀ attheinitialtimeepocht₀ .Theprocedureofnumerical integrationhastopartitionthetimeinterval[t₀ ,t_yi]intoanumberofsub-intervals,usuallyequidistantsuchthat

tj=t0 +jh, j=1,2,...,myi

whereh=

(

tyi−t0

)

/myi.Onecanthenapplynumericalintegrationmethodstoprogressivelycomputeallthenominal

refer-encepositionsx0

(

tj, p0 , c

(

p0 , x0 ₀ , v0 ₀

))

.

withalinearized,one-iterationsolution.Thus,wewillfocusonexplicitnumericalintegrationmethodstoprogressivelysolve thenonlineardifferentialEq.(6)intheremainderofthissection.

Inwhatfollows,wewillusetheEulermethodandthemodifiedEulermethodtodemonstratetheconstructionofx(tj,p,

c(p,x₀ ,v₀ ))intermsof

z₀ and

p.OthernumericalintegrationmethodssuchasHeun’smethod,Runge–Kuttamethods ofanyorderand/ortheNewton–Cotesmethodcanbetreatedinthesamemannerandwillbeomittedhere.Theinterested readercanworkthemoutbyhimselforherself.Forconcisenessofnotations,wewilldenotetherighthandsideof(19)by

g(t,z(t),p)andrewrite(19)asfollows:

˙

z

(

t

)

=g

(

t,z

(

t

)

,p

)

(24)

underthe(unknown)initialconditionsz₀ (namely,x₀ andv₀ ). TostarttheEulermethod,wehave

z

(

t1

)

=z

(

t0

)

+hg

(

t0 ,z

(

t0

)

,p

)

, (25)

(seee.g.,[100,106]).Linearizingthevectorfunctionsg(·)at

(

z0 ₀ , p0

)

andbearinginmindtheapproximateorbitz0 (t,z0 (t),

p0 ),wecanrewrite(25)into:

z0

(

t1

)

+

z

(

t1

)

=z0 ₀ +

z0 +hg

(

t0 ,z0 ₀ ,p0

)

+hGgz0

z0 +hGgp0

p,

orequivalently,

z

(

t1

)

=

δ

z₀₁ 0 +

I6 +hGgz0

z0 +hGgp0

p, (26)

where

δ

z0 ₀₁ =z0 ₀ +hg

(

t0 ,z0 ₀ ,p0

)

−z0

(

t1

)

, (27a)

Ggz0 =

∂

g

(

t,z

(

t

)

,p

)

∂

zT

_z

(t)=z0 0,p=p0

, (27b)

and

Ggp0 =

∂

g

(

t,z

(

t

)

,p

)

∂

pT

_z(t)

=z0 0,p=p0

. (27c)

Toprogressfromt₁ tot₂ ,theEulermethodtakestheformofformula:

z

(

t2

)

=z

(

t1

)

+hg

(

t1 ,z

(

t1

)

,p

)

. (28)

Following the same procedure as described above, and neglecting all the terms of order h2 , we can rewrite the above formulaasfollows:

z

(

t2

)

=

δ

z0 ₁₂ +

I6 +hGgz1

z

(

t1

)

+hGgp1

p, (29)

where

δ

z0 ₁₂ =z0

(

t1

)

+hg

(

t1 ,z0

(

t1

)

,p0

)

−z0

(

t2

)

.

Inserting(26)into(29)andaftersomerearrangement,wehave

z

(

t2

)

=

δ

z0 ₁₂ +

δ

z0 ₀₁ +hGgz1

δ

z0 01 +

I6 +hGgz0 +hGgz1

z0 +[Ggp0 +Ggp1 ]

p

=

δ

z0 ₀₂ +hGgz1

δ

z0 01 +

I6 +h

1

j=0

Ggz j

z0 +h

1

j=0

Ggp j

p, (30)

where

δ

z0 ₀₂ =z0 ₀ +h

1

j=0

g

(

tj,z0

(

tj

)

,p0

)

−z0

(

t2

)

.

ThematricesG_gz1 andG_gp1 arecomputedinthesamemannerasin(27b)and(27c)butatthepointofz0 (t₁ ).

Repeatingthe same procedure as describedin the above, we can finally obtain the representationof the corrections

z(tyi)asfollows:

z

(

tyi

)

=

δ

z0 0 tyi+h

myi−1

j=1

Ggz j

δ

z0 ₀ j+

I6 +h myi−1

j=0

Ggz j

z0 +h myi−1

j=0