205
1 INTRODUCTION
Theionosphereisthemajorerrorsourcethataffects
positioning accuracy of GPS users. GPS satellites
broadcast the Klobuchar parameters to correct the
ionospheric delay. But its accuracy is notenough to
meet the required navigation performance in many
cases. Many countries have been developed
augmentationsystemssuch
asGBAS,SBAStoachieve
higher navigation performance. In SBAS(satellite‐
based augmentation system), like WAAS, they
estimate local ionospheric delays more accurately
usingsparenetworkofreferencestationsandprovide
estimated result to users. To estimate the local
ionospheric delays, grid‐based models are generally
used. Grid‐based model estimates
the ionospheric
delaysatthepre‐definedgridpointsandprovidethe
valuestouser.Thenusersgetionosphericdelaysby
interpolating received values close to their
ionospheric pierce point. This grid model usually
works well. But,in ionosphericstorm condition,the
ionosphericspatialgradientbecomeshigherandlocal
ionospheric delays
have irregularity structures.
Because of its interpolating method, grid‐based
Modeling of Ionospheric Delay for SBAS Using
Spherical Harmonics Functions
D.Han,H.Yun&C.Kee
M
echanicalandAerospaceEngineeringandtheInstituteofAdvancedAerospaceTechnology,SeoulNationalUniversity,
Korea
ABSTRACT: In SBAS (satellite‐based augmentation system), it is important to estimate ionospheric delay
accurately to guarantee userʹs accuracy and integrity. Grid based ionospheric models are generally used to
estimateionosphericdelayforSBAS.In
gridbasedmodel,SBASbroadcastsverticalionosphericdelaysatthe
gridpoint,andusersgettheir ionospheric delay byinterpolatingthosevalues.Ionosphericmodel based on
spherical harmonics function is another method to estimate ionospheric delay. This is a function based
approachandsphericalharmonicsfunctionisa2 ‐Dfourier
series,containingtheproductoflatitudedependent
associated Legendre functions and the sum of the longitude dependent sine and cosine terms. Using
ionospheric delay measurements, coefficients for each spherical harmonics functions are estimated. If these
coefficients are known, user can reconstruct ionospheric delay. In this paper, we consider the spherical
harmonics based model and propose a ionospheric delay estimation strategy for SBAS that can be used to
mitigateionosphericdelayestimationerror,especiallyinstormcondition.First,coefficientsareestimatedunder
initialorderanddegree.Thenresidualerrorsforeachmeasurementaremodeledbyhigherorderanddegree
terms, then coefficients
for these terms are estimated. BecauseSBAS messagecapacity is limited, innormal
condition,initialordertermsareonlyusedtoestimateionosphericdelay.Ifionosphericstormisdetectedand
thereisneedtomitigatetheerror,higherordertermsarealsousedanderrorcanbedecreased.Tocompare
the
accuracy of spherical harmonics based model with grid based model, some post‐processing test results are
presented. Raw observationdata isobtained fromRINEX format and the root mean square(RMS) and max
valueofresidualerrorsarepresented.
http://www.transnav.eu
the International Journal
on Marine Navigation
and Safety of Sea Transportation
Volume 7
Number 2
June 2013
DOI:10.12716/1001.07.02.07
206
models have problem estimating those irregularity
structuresandmodelingerrorsareincreased.
Spherical harmonics model is another option for
modelingionosphericdelay.Incontrasttogrid‐based
model, Spherical harmonics model is analytic
model[1]. Spherical harmonics are used in many
theoretical and practical applications such as
gravitational fields modeling and
geoids modeling.
There are some researches that use spherical
harmonicsforionospheremodeling[1,2,5,6].Inthis
paper,wecarriedoutpost‐processingtesttocompare
the performance of spherical harmonics model and
grid model. In this test, we only conducted the
accuracyperformancetest.
2 MEASUREMENTOUTLIER
DETECTION
SBAScollectsdual‐frequencyGPSobservationsfrom
reference stations. Sometimes, observed data have
outliers and these outliers affect the modeling
performance.Toimprovetheperformance,itmustbe
detectedandremoved.Inthispaper,polynominalfit
method is used to detect the outliers [3]. For each
satellite,previous30
epochdata are fittedin second
order.Residualsarecalculatedfromfittedvalueand
current epoch data. If calculated residuals exceed
threshold,thesedataareclassifiedasoutlier.Figure1
shows raw vertical ionospheric delay of PRN 16
satelliteandfigure2showsverticalionosphericdelay
afterremovingoutliers.
Figure1. vertical ionospheric delay of PRN 16 satellite
(beforeremovingoutlier)
Figure2.verticalionosphericdelayofPRN16satellite(after
removingoutlier)
3 IONOSPHERICDELAYMODELINGWITH
SPHERICALHARMONICSFUNCTION
3.1 Sphericalharmonicsfunction
SphericalharmonicsfunctionisbasedontheFourier
seriesexpansioninlongitudeandLegendrefunction
in latitude as the orthogonal basis functions[2]. The
mathematicalexpressionofverticalionosphericdelay
usingsphericalharmonicsisasfollows.
00
( , ) (sin )( cos( ) sin( ))
nm
vnmnmnm
nm
If P fC m S m




(1)
is geomagnetic latitude of an IPP
is geomagnetic longitude of an IPP
, are int eger degree and order of Legendre functionnm
, are unknown spherical harmonics coefficients
nm nm
CS
are normalized associated Legendre functions
nm
P
Like the grid model, this approach also assumes
that electron density is concentrated in a very thin
shellat350kmaltitude.Thedegreeandorder of the
Legendre function defines the resolution of the
model. Using above equation, local ionospheric
delays can be represented by spherical harmonics
coefficients.
3.2 Coefficient
estimationusingleastsquaremethod
To model the local ionospheric delay with spherical
harmonics functions, vertical ionospheric delay
measurements are collected. Then the least square
methodisusedtoestimatetheunknowncoefficients.
Thesystemcanberepresentedasthesematrices.
zHx
(2)
is vertical ionospheric delay measurementz
H is geometry matrix
is unknown coefficientsx
   
 
00 1 1 00 1 1 1 1 1 1
00 2 2 00 2 2 2 2 2 2
00 00
sin cos 0 sin sin 0 sin cos sin sin
sin cos 0 sin sin 0 sin cos sin sin
sin cos 0 sin sin 0 sin cos sin sin
nm nm
nm nm
NN NNnmN NnmN
PP PmPm
PP PmPm
H
PP PmPm








N
00 00
T
nm nm
xC S C S
In general, this equation is over determined and
can be solved. The unknown coefficients are
computedasfollows
1
()
TT
x
HH Hz
(3)
Estimatedcoefficientswillbebroadcastedtouser,
thenusercanreconstructthelocalionosphericdelay
usingthesecoefficients.
207
4 TESTRESULTSANDDISCUSSION
4.1 Data
GPS RINEX format data taken from 32 stations are
consideredforanalysis.Thesedatawereobtainedvia
the ftp servers of CORS (ftp://
www.ngs.noaa.gov/cors/rinex). To investigate the
performanceinionosphericstormcondition,weselect
thedataobservedinsuchcondition.OnOctober29
th
,
2003, there is a severe ionospheric storm[4] and we
used these data. For example, figure 3 shows the
verticalionosphericdelaysonthatday.Inthisfigure,
2‐Dcubicpolynominalfunctioninmatlabisusedto
represent the raw measurements. The maximum
verticaldelayreach30metersand
spatialgradientis
veryhighinsomeregion.
Figure3. Vertical Ionospheric delay on October 29
th
, 2003.
(GPStime334320,2‐Dcubicpolynominal)
4.2 Methodology
At each epoch, vertical ionospheric delays are
estimatedusinggridmodel,andsphericalharmonics
model. Three methods are implemented, one is grid
modelandothersaresphericalmodel.
Method1 justusedGridmodeland64gridpoints
areusedinthismethod.
Method2usedsphericalharmonicsmodel.
When
the spherical harmonics model is used, the total
number of coefficients will be fixed depending on
pre‐determined order and degree. Because message
capacityislimited,itisimportanttodetermineproper
numberofparametersconsideringbothaccuracyand
message capacity. In this method , both order and
degree
aresetto5andtotalcoefficientsare36.
Method3usedsphericalharmonicsmodelwith81
coefficients(both order and degree are set to 8) to
estimate the ionospheric delay more accurately in
severecondition.IncontrasttoMethod2,estimation
processisdividedintotwostepsinthismethod.First
step is repetition of Method 2. In second step, the
residuals are calculatedfrom the measurements and
estimated value. Then higher order and degree
coefficientsareestimatedfromtheseresiduals.Atthis
time,thegeometrymatrixissameasinfirststepand
equationsareasfollows.
L
rzHx
(4)
H
H
x
1
()
TT
H
x
HH Hr
(5)
z is vertical ionospheric delay measurement
r is residual
is lower degree and order coefficients
L
x
is higher degree and order coefficients
H
x
Ratherthanestimatingallcoefficientsatonce,we
proposed this method to provide correction
informationtousermorefastermaintainingaccuracy
performance. When this method is used to estimate
spherical harmonics coefficients, users can roughly
compute ionospheric delay after they only receive
first 36 coefficients. After receiving all coefficients,
user
can improve the accuracy. In this method, our
proposed broadcasting strategy is as follows. Check
the residuals and if the accuracy of the first step is
enough, broadcast the only first 36 coefficients to
user.Ifresidualerrors aregrowandthereisaneedto
estimate the ionospheric delay more
accurately,
higher order and degree terms will also be
broadcastedtouser.
Figure4.Diagramformethod3
_______________________________________________
MethodModelNumberofParameters
_______________________________________________
1Grid64
2SHF36
3SHF81
_______________________________________________
Figure 5, 6, 7 respectively shows the estimated
verticalionosphericdelaysforeachmethod.
208
Figure5.EstimatedVerticalIonosphericdelay”
(GPStime334320,Method1)
Figure6.EstimatedVerticalIonosphericdelay
(GPStime334320,Method2)
Figure7.EstimatedVerticalIonosphericdelay
(GPStime334320,Method3)
4.3 Accuracyperformancetest
To compare the accuracy performance, residualsare
calculated for every measurement. Then maximum
value and root mean square(rms) value of residuals
arecomputedforeach epoch. This testisconducted
for five hours. In this interval, ionosphere condition
getsworseastimeflows.
Figure 8 and figure
9 shows the time history of
maximum and rms value of residual. When the
spherical harmonics model is used rather than grid
model,accuracyperformanceisimproved.Especially,
method 2 shows slightly better performance
compared with method 1, even though the total
number of parameters used in method 2 is
about a
half of method 1. The average improvement is 1.42
meters in maximum value and 0.24 meters in rms
value for method 2, 1.60 meters in maximum value
and 0.31 meters in rms value for method 3. The
improvementis more noticeableinlatter part of the
figurecorrespondingto
severeionospherecondition.
Figure8.Timehistoryofmaximumvalue
Figure9.Timehistoryofrmsvalue
5 CONCLUSIONS
Inthispaper,wemodeledthelocalionosphericdelay
using spherical harmonics functions. In ionospheric
storm condition, the vertical ionospheric delays are
collectedfrom32stations.Thenoutlierdetectionand
removal is conducted on these measurements.
Unknown spherical harmonics coefficients are
estimated by least square method. Using these
coefficients,
residuals for every measurement are
estimated to test the accuracy performance. These
209
resultsarecomparedwiththeestimatedresidualsof
gridmodel.Sphericalharmonicsmodelshowsbetter
accuracy performance with lower number of
parameters especially in more severe condition. We
also proposed conceptual estimation and
broadcasting strategy that can be used to provide
correction information more faster when we use
sphericalharmonics
modelforSBAS.
ACKNOWLEDGEMENT
This research was a part of the project titled “WA‐
DGNSS Development” funded by the Ministry of
Land,TransportandMaritimeAffairs,Korea.
REFERENCES
[1]D.VenkataRatnam,CH.Sufjatha,A.D.Sarma,Shudha
Ravindran, “Modelling of GAGAN TEC data using
SphericalHarmonicFunctions”, InternationalConverence
onComputersandDevicesforCommunication,2009
[2]Yongjin Moon, “Evaluation of 2‐Dimensional Ionosphere
Models for National And Regional GPS Networks in
Canada”,M.Sc.Thesis,UniversityofCalgary,2004
[3]Jiyun Lee, Sungwook Jung, Eugene Bang, Sam Pullen,
PerEnge,“Longtermmonitoringofionosphericanomaliesto
support the local area augmentation system”, ION GNSS
2010,Potland,September2010,pp.2651‐2660
[4]AlexandruEne,DiQiu,MingLuo,SamPullen,PerEnge,
“AComprehensive IonosphereStormData AnalysisMethod
to Support LAAS Threat Model Development”, ION NTM
2005,SanDiego,CA,24‐26January2005, pp.110‐130
[5]V.B.S.SRILATHAINDIRADUTT,G.SASIBHUSHANA
RAO,“2‐DElectronDensityProfileoftheIonosphereUsing
the Modified Spherical Harmonics Analysis”, The NEHU
Journal,VolIX,No.2,July2011,pp.87‐96
[6]Byung‐KyuChoi,WooKyungLee, Sung‐Ki Cho, Jong‐
Uk Park, Pil‐Ho Park, “Global GPS Ionospheric modeling
using spherical harmonics expansion approach”, Jounal of
AstronomyandSpaceSciences,2010,pp.359‐366