IEICE TRANS. INF. & SYST., VOL.E92–D, NO.7 JULY 2009

1382

PAPER

Special Section on Large Scale Algorithms for Learning and Optimization

Estimation of a Long-Term Variation of a Magnetic-Storm Index
Using the Merging Particle Filter
Shin’ya NAKANO†∗a) , Nonmember and Tomoyuki HIGUCHI†∗ , Member

SUMMARY
The D st index is the most popular measure of a scale ofmagnetic storms, and it is widely used as a monitor of the conditions of the
Earth’s magnetosphere. Since the D st index contains contributions frommultiple magnetospheric phenomena, it is important to distinguish each of
the contributions in order to obtain meaningful information about the conditions of the magnetosphere. There have been several eﬀorts which modeled temporal evolution of the D st index empirically, and these empiricalmodels considers some contributions separately. However, they take only
short-term varations into accout, and contributions from phenomena whichshow long-term variations are neglected. In the present study, we have developed a technique for estimating the component of long-term variationsof the D st index using solar wind data and a nonlinear empirical model. The
newly-developed technique adopts an algorithm which is similar to the particle filter. This algorithm allows an on-line processing of a long sequenceof D st data, which would enable a real-time estimation of system variables
in a nonlinear system model. The estimates of the long-term variations
can be used for accurate estimation of other contributions to the D st index,
which would provide credible information about the conditions of the magnetosphere. The framework proposed in the present study could be applied
for the purpose of continuous real-time monitoring of the environment of
the magnetosphere.
key words: merging particle filter, time-series analysis, geomagnetic data

1.

Introduction

Magnetic storm is a global phenomenon which causes great
decrease in the magnetic field on the ground at low and midlatitudes. The degree of the decrease is commonly measured
by D st index [1], which is derived from magnetic-field data
at four diﬀerent ground observatories at mid latitudes. The
D st index represents a deviation of the magnetic field at low
latitudes from a normal state. The value of D st is thus nearzero in a normal state. During magnetic storm, it decreases
typically by 50–100 nT and by hundreds of nT for strong
magnetic storm. The D st index is calculated for every hour,and its hourly data have been accumulated and published forabout 50 years since 1957.
It is widely accepted that the decrease of D st index during magnetic storm is caused by an electric current flowing westward in the magnetosphere, which is called the ring
current (Fig. 1; see also, e.g., [2]). It is also well known
that an electric current flowing eastward on the outer boundary of the magnetosphere, which is called the magnetopause
Manuscript received December 5, 2008.
Manuscript revised March 7, 2009.
†
The authors are with the Institute of Statistucal Mathematics,
ROIS, Tokyo, 106–8569 Japan.
∗
The authors are also with Japan Science and Technology
Agency, Kawaguchi-shi, 332–0012 Japan.
a) E-mail: shiny@ism.ac.jp
DOI: 10.1587/transinf.E92.D.1382

Fig. 1 Schematic picture of the sun-earth system which is focused on in
the present paper. This picture is a projection onto the Earth’s equatorial
plane.

current, also has a significant eﬀect on the value of D st as
depicted in Fig. 1. The enhancement of the magnetopause
current causes the increases of D st , while the enhancementof the ring current causes the decrease of D st . Hence, thevariations of D st are mostly a mixture of the ring current effect (RC eﬀect) and the magnetopause current eﬀect (MPCeﬀect).
The variations of the ring current and the magnetopause current are closely associated with the conditionsof the solar wind, which streams outward from the sun, although these two currents depend on diﬀerent parameters of
the solar wind. Burton et al. [3] have empirically modeled
the temporal evolution of the D st index using some solar
wind parameters. They decomposed the D st variations into
the RC eﬀect and the MPC eﬀect, and then consider the
temporal variation of each of two eﬀects separately. They
assumed the MPC eﬀect to be given as a function of one of
the solar-wind parameters, the solar-wind dynamic pressure
Pd . The RC eﬀect was described as a nonlinear evolutive
system which requires the solar-wind electric field as an input.
Their model (hereafter we will refer to it as the
Burton’s model) was constructed under the assumption that
all the variations of D st are due to the RC eﬀect and/or the
MPC eﬀect, and they assumed other contributions to the D st
index to be constant, that is, independent from time. However, it is not trivial whether a variation of another contribution is really negligible or not. In particular, a long-term
magnetic field variation which can not be represented by the
Burton’s model might make some significant influence overthe D st index. As a matter of fact, we can estimate whether
another factor may significantly contribute to D st variations

c 2009 The Institute of Electronics, Information and Communication Engineers
Copyright NAKANO and HIGUCHI: ESTIMATION OF A LONG-TERM VARIATION

1383

on the basis of the Burton’s model itself using long-term
data of the D st index and the solar wind.
The purpose of the present paper is to show how we
can evaluate another contribution to D st variations and to
estimate how large such a contribution is. The estimation ofanother contribution is conducted using the merging particle
filter (MPF) [4]. This algorithm is applicable to nonlinear
system models like the Burton’s model and it is eﬀective
for applying to a very long-term sequence of data because
it is an on-line algorithm which can process long-term data
sequentially.
2.

Base Model of the Temporal Variation of the D st Index

The D st index is influenced by both the RC eﬀect and the
MPC eﬀect. Burton et al. took an approach to model the
MPC eﬀect as the first step and then eliminate it from D st to
obtain a pure RC eﬀect:
DRC (nT) = D st − b Pd + c(1)
where DRC represents the pure RC eﬀect, Pd denotes the
dynamic pressure of the solar wind which is observed by
spacecraft out of the magnetosphere, and b and c are theparameters which should be given a priori.√The second termin the right-hand side of this equation, b Pd , corresponds
to the MPC eﬀect, and the third term c corresponds to othercontributions than the RC eﬀect and the MPC eﬀect. Thismeans that the√Burton’s model assumes that the MPC eﬀect
is given by b Pd and that other contributions is given by
the constant c. The temporal evolution of the pure RC eﬀect
DRC can then be modeled as follows:
ΔDRC
DRC
= Q−
Δt
τ

(2)

where the parameters Q and τ should be given a priori.
The original paper by Burton et al. provided the optimal value of each constant coeﬃcient in Eqs. (1) and (2).
However, the optimal values have been revised by O’Brien
and McPherron (2000) [5] using a much longer term of thedata. According to them, the parameters b and c in Eq. (1)can be given as constants:
b = 7.26 nT/nPa−1/2
c = 11 nT.

(3)
(4)

The parameter Q, which represents the evolution of DRC ,can be given as follows:
Q (nT/hour) = −4.4 H(E − 0.49)
where H denotes the Heaviside function as⎧⎪⎪⎨1 if x ≥ 0
H(x) = ⎪⎪⎩0 otherwise.
The variable E is given as(5)

E = vBs(7)

where v is the solar-wind velocity and Bs is the southward
component of the solar-wind magnetic field. Both v and Bs
can be observed by spacecraft out of the magnetosphere aswell as Pd . Finally, τ is given as follows:
9.74
τ (hours) = 2.4 exp.
(8)
4.69 + H(E)
Once an initial value of the D st index is given, we canpredict the evolution of D st sequentially using the model
given by Eqs. (1) and (2) as long as the solar-wind data are
available. Hence, by comparing the prediction with the actual D st values, we can validate the Burton’s model and thevalues of coeﬃcients obtained by O’Brien and McPherron
and evaluate contributions from other factors than the RC
eﬀect and the MPC eﬀect. In the following sections, we
perform the prediction by putting the solar-wind data into
the Burton’s model.
3.

Modeling of the System

3.1

System Model

In order to perform the prediction of the D st index, we construct a state space model on the basis of the Burton’s model.
We decompose a state of D st at time t into three components
asDtst = DtRC + DtMPC + Dtres(9)

where Dtst denotes a model D st value at time t, DtRC denotes
the RC eﬀect on D st , DtMPC denotes the MPC eﬀect, and Dtres
denotes a residual eﬀect other than the RC and MPC eﬀects
on D st , which corresponds to −c in Eq. (1) and is of our maininterest in the present paper.
According to Eqs. (2) and (5), we describe a transitionof a state of DRC for an hour as
t−1
DtRC = Dt−1
− 0.49) −RC − 4.4 H(E

Dt−1
RC
.
τt

(10)

The MPC eﬀect DtMPC is assumed to be given as a functionof Ptd as
DtMPC = b

Ptd .

(11)

As for a transition of a state of Dres , two models are considered in order to evaluate whether a variation of Dres isnegligible or not. One model assumes Dres to be constant as
Dtres = −11.

(12)

according to Eq. (4). The other model is given as(6)

Dtres ∼ N(Dt−1
res , 0.01).

(13)

We also consider state transitions of the two solar-wind

IEICE TRANS. INF. & SYST., VOL.E92–D, NO.7 JULY 20091384

parameters Pd and E. As described above, these two parameters can be observed by spacecraft. However, spacecraft
basically observe local structures of the solar wind and the
local structures observed by spacecraft do not necessarily
agree with large-scale solar-wind structures which controls
the conditions of the Earth’s magnetosphere. Then, we assume that ‘eﬀective’ Pd and E are uncertain and we include
them in variables to be estimated. Thus, it is necessary tomodel the state transitions of them. The state of Pd at time t
is assumed to obey a log-normal distribution aslog Ptd ∼ N(log Pt−1
d , 0.02)(14)

because the solar-wind dynamic pressure Pd can not beless than zero. The variance of Pd was determined by amaximum-likelihood method using a set of OMNI2 hourly
solar-wind data P̄d . The OMNI2 solar-wind data were provided on the OMNIWeb database of National Space Science Data Center, NASA (http://omniweb.gsfc.nasa.gov/)
and they are also referred to Sect. 3.2. A transition of E
is represented using a Cauchy distribution in order to allow
large jumps which are sometimes observed in E. We assume
that E t obey the Cauchy distribution with a location parameter of E t−1 and a scale parameter of 1. The scale parameterfor the transition of E are given subjectively. However, since
it is associated with short-term variations, it does not makeany significant eﬀects on the estimates of Dres . In addition,
we assume that τ is time-dependent. We give τt as


9.74
t
(15)
τ = 2.4 exp
4.69 + H(Et−1 )where
Et−1 ∼ N(E t−1 , 0.25).

(16)

3.2 Observation Model
From the system described above, we can obtain observations of three variables. Space craft observations provide
the data of Pd and E at each hour, although these data aresometimes lost. We define the observation model of Pd andE at each hour asP̄td ∼ N(Ptd , 1)Ē t ∼ N(E t , 4)D̄tst = Dtst + wtDst= DtRC + b

Ptd + Dtres + wtDst(19)

where D̄tst denotes the observed D st index and wtDst is theobservation error contained in the D st data which obeys thenormal distribution N(0, 25).
3.3

State Space Model

For the convenience in the following section, we define astate vector as follows:⎛ t ⎞⎜⎜⎜ DRC ⎟⎟⎟⎜⎜⎜⎜ Dt ⎟⎟⎟⎟⎜⎜⎜ res ⎟⎟⎟xt = ⎜⎜⎜⎜⎜ Ptd ⎟⎟⎟⎟⎟ .
(20)
⎟⎟
⎜⎜⎜⎜
t ⎟
⎟
⎜⎜⎝ E ⎟⎟⎠
τt
Since the MPC eﬀect DtMPC is assumed to be given as a function of Ptd , DtMPC is not included in a state vector. Using thisstate vector, we can represent a transition of a state xt justby a conditional distribution as
xt ∼ p(xt |xt−1 ).

(21)

where p(xt |xt−1 ) can be defined by combining Eqs. (10)–(16). We also define an observation vector yt as⎛ t ⎞⎜⎜⎜ D̄ st ⎟⎟⎟⎟⎟⎜⎜yt = ⎜⎜⎜⎜ P̄td ⎟⎟⎟⎟ .
(22)⎝⎜ t ⎠⎟
Ē
However, the data of P̄td and Ē t are sometimes missing as
mentioned above. In such cases, we redefine an observation
vector yt as(23)
yt = D̄tst .
In either case, the observation model described in Sect. 3.2
can be written as
yt ∼ p(yt |xt ).

(24)

(17)
(18)

where P̄td and Ē t denote the observations of Pd and E ateach hour, respectively. In this study, we refer to the OMNI2
solar-wind hourly data as the data of P̄td and Ē t . Although
the original OMNI2 data do not contain Ē t data, we generate Ē t data from data of the solar-wind velocity v̄t and thesouthward component of the solar-wind magnetic field B̄tz
as Ē t = v̄t B̄tz . We can also use the values of the D st index ateach hour. The data of the D st index are provided by Data
Analysis Center for Geomagnetism and Space Magnetism,
Kyoto University. The observation model of the D st index isdefined as follows:

Fig. 2 Dependencies among the variables considered in the present
study. The variables which can be observed are indicated by black boxes
with white letters.

NAKANO and HIGUCHI: ESTIMATION OF A LONG-TERM VARIATION

1385

In ordinary cases, p(yt |xt ) can be defined by combining
Eqs. (17)–(19). If the data of either P̄td or Ē t are missing,
only Eq. (19) is evaluated; that is, p(yt |xt ) is determinedfrom Eq. (19).
Figure 2 illustrates the dependencies among the variables considered in the present study. The variables includedin the state vector xt are enclosed within a shaded box. The
observable variables are indicated by black boxes with whiteletters.
4.

An approximation of the filtered PDF p(xt |y1:t ) can beobtained by incorporating an observation yt into the particle
approximation of a predictive PDF {x(i)}N as follows:t|t−1 i=1p(xt |y1:t )p(xt |y1:t−1 ) p(yt |xt )=p(xt |y1:t−1 ) p(yt |xt )dxt
≈Algorithm

We conduct the estimation of the value of the state vector xt
at each hour from the observations yt on the basis of the system and observation models described in the previous section. As the system model described above is nonlinear, theestimation of a state is done using the merging particle filter (MPF) [4], which is applicable even to nonlinear systemmodels. The MPF is an algorithm based on the particle filter (PF) [6]–[8]. The MPF provides an approximation of theposterior probability density function given a sequence ofobservations {y1 , . . . , yT }, like the PF. However, it is much
more eﬃcient than the normal PF. In the following, the algorithm of the MPF is briefly reviewed.
In the MPF, a probability density function (PDF) p(x)is approximated by a set of N samples {x(1) , . . . , x(N) } as alsodone in the particle filter asp(x) ≈N
1  
δ x − x(i) .
N i=1

(25)

Each sample is called a ‘particle’ and an approximation ofa PDF using a set of particles is called a ‘particle approximation’. In the normal PF, a set of particles for representing
a posterior PDF contains many duplicates of the same particle, which causes highly redundant computational cost. On
the other hand, in the MPF, each particle is generated bymerging multiple samples in order to maintain the diversityof particles. This allows us to reduce the redundant computational cost.
We sequentially obtain particle approximations of apredictive PDF p(xt |y1:t−1 ) and a filtered PDF p(xt |y1:t ) at
each time step as follows. Suppose that a PDF p(xt−1 |y1:t−1 ),which is a posterior PDF given the observations until time}N ast − 1, is approximated by a set of particles {x(i)t−1|t−1 i=1
N

1  
p(x) ≈
δ xt−1 − x(i).
t−1|t−1
N i=1

(26)

sample taken from a conditional distribution p(xt |x(i)),
t−1|t−1
this set of samples {x(i)}N oﬀers a particle approximationt|t−1 i=1
of a predictive PDF at the next time, p(xt |yt−1 ), as
N

1  
δ xt − x(i)
.
t|t−1
N i=1

(27)

1


( j)

yt |xt|t−1

N


 

p yt |x(i)δ xt − x(i)
t|t−1
t|t−1
i=1



wi δ xt − x(i)
t|t−1

(28)

i=1

) is the likelihood of x(i)given the data yt ,where p(yt |x(i)
t|t−1
t|t−1
which can be calculated according to Eq. (24). The weight
wi is defined as)p(yt |x(i)t|t−1.
wi = ( j)j p(yt |xt|t−1 )(29)

Although Eq. (28) oﬀers an approximation of the filtered
PDF p(xt |y1:t ), we generate a new set of particles {x(i)}Nt|t i=1
which represents the filtered PDF in the following form:p(xt |y1:t ) ≈

N

1  
δ xt − x(i)
.
t|t
N i=1

(30)

In the MPF, each particle for representing the filtered PDF
is obtained by combining multiple particles taken from a
set of particles for representing the predictive PDF; thatis obtained from multiple samples fromis, a particle x(i)
t|t
{x(i)
}N . Here the number of particles to be combined ist|t−1 i=1
taken as 3. Then, in order to obtain a set of N particles}N , it is necessary to draw 3N samples with weights of{x(i)t|t i=1
wi { x̂(1,1), · · · , x̂(3,1), · · · , x̂(1,N), · · · , x̂(3,N)}. Each particle in
t|t
t|t
t|t
t|t
the new set, {x(i)}, is generated as a weighted sum of 3 samt|t
ples contained in this 3N samples as:
=
x(i)
t|t

n


α j x̂(t|tj,i) .

(31)

j=1

The set of weights {α j }nj=1 in Eq. (31) is set to satisfynαj = 1(32a)

α2j = 1

(32b)

j=1
n


If we obtain a set of particle where each particle x(i)
is a
t|t−1

p(xt |y1:t−1 ) ≈

=

j p

N




j=1

where α j ∈ R for all j. Eqs. (32a) and (32b) ensures thatthe new set of particles {x(i)}N has asymptotically the samet|t i=1
average and covariance of the filtered PDF p(xt |y1:t ) for
N → ∞.

IEICE TRANS. INF. & SYST., VOL.E92–D, NO.7 JULY 2009

1386

Fig. 3 Diﬀerence between the model D st and the actual D st for the casethat Dres is assumed to be constant.

5.

Fig. 4 Comparison between the model D st (black solid line) and the actual D st (gray solid line) from January 1999 to December 1999 for the casethat Dres is assumed to be constant.

ResultWe estimated the temporal evolution of the state vector xt
for the period from 1995 to 2005. The estimation was
performed for the period from 1995 because many of theOMNI2 solar-wind data are missed before 1995. An estimate of xt at each time step was provided from the ensem}N , which approximates a filtered
ble mean of particles {x(i)t|t i=1
PDF p(xt |y1:t ). The number of particles N used for the estimation was 1600. In order to evaluate whether the temporal
evolution was successfully estimated or not, we compare theD st evolution estimated using the MPF with a sequence ofthe real D st data.
Figure 3 shows diﬀerence between the real D st data
and the model D st for the case that Dres is assumed to beconstant where positive means that the model D st is larger
than the actual D st . Although the diﬀerence sometimes becomes larger than 20 nT, it mostly varies with the amplitude
of about 10 nT during the period from 1995 to 2005. However, some long-term variation is also seen in this figure.
In particular, the model D st is prominently smaller than theactual D st around the middle of 1999. Figure 4 compares
the model D st for constant Dres with the actual D st during
one year from January 1999 to December 1999. The grayline indicates the actual D st and the black line indicates themodel D st . From the end of May to July, the model D st underestimated the actual D st by about 10 nT.
Figure 5 compares the model D st for the case that Dres
is variable according to Eq. (13) with the actual D st for theperiod from January 1999 to December 1999. The grayline indicates the actual D st and the black line indicates themodel D st as in Fig. 4. The model D st well agree with theactual D st . Indeed, while the log-likelihood of the model
based on Eq. (12) was about −1.8 × 10−8 , the log-likelihoodof the model based on Eq. (13) was about −1.5 × 10−5 . The
dashed line shows the estimated Dres . The residual eﬀect
Dres exceeds zero from the middle of the May to the middleof the September, while the empirical model assumes that
Dres is negative (Dres = −11 nT). This result means that theresidual eﬀect is more variable than assumed by the previous studies. The Burton’s model sometimes overestimate
Dres and sometimes underestimate Dres by tens of nT, which

Fig. 5 Comparison between the model D st (black solid line) and the actual D st (gray solid line) from January 1999 to December 1999 for the casethat Dres is modeled according to Eq. (13). The estimated Dres (see text) is
also shown in a dashed line.

inevitably causes misestimation of DRC . As described in
Introduction, D st varies typically by 50–100 nT during magnetic storms. Thus, the improvement of the estimation ofDres is not negligible in order to evaluate the RC eﬀect using
the D st index.
6.

Concluding Remarks

We have presented a new technique to estimate the contribution from an unconsidered factor in the Burton’s model
which describes temporal evolution of a measure of magnetic storms, the D st index. The estimation is performed
using the merging particle filter which is applicable to nonlinear system models like the Burton’s model. By considering a long-term variation in addition to the RC eﬀect andthe MPC eﬀect which are already considered in the Burton’smodel, the agreement between the model D st value and theactually observed D st value is greatly improved. This factmeans that a long-term residual eﬀect Dres significantly contributes to the D st index. In order to accurately estimate thetemporal variation of the RC eﬀect DRC , we need to accurately estimate the variation of Dres . Thus, the techniquewhich we presented above would serve as an important tool
for analyzing the temporal evolution of DRC index.
Recently, continuous monitoring of magnetic storm activity becomes important because it is associated with theenvironment of the magnetosphere where many spacecraft

NAKANO and HIGUCHI: ESTIMATION OF A LONG-TERM VARIATION

1387

operate and with the conditions of the ionosphere. Thus,a real-time monitor of magnetic storm activity and the ring
current is important. In the present study, we adopted themerging particle filter (MPF). The MPF is an on-line algorithm applicable for processing even real-time data and it enables an on-line estimation of the long-term residual eﬀect
Dres . This on-line estimation provides valuable information
for real-time estimation of the ring current variation DRC .
Thus, the framework adopted in the present paper could be
useful for a real-time monitoring of the magnetospheric environment using the real-time edition of the D st index.
Acknowledgment
This study was partially supported by the Japan Science andTechnology Agency (JST) under the Core Research for Evolutional Science and Technology (CREST) program. We
are grateful to World Data Center for Geomagnetism, Kyoto
University for providing the D st index. We also thank National Space Science Data Center, NASA for providing theOMNI2 solar-wind data.
References[1] M. Sugiura and T. Kamei, Equatorial D st index 1957–1986, ISGI Publisher Oﬀ., Saint-Maur-des-Fosses, France, 1960.
[2] W. Baumjohann and R.A. Treumann, Basic space plasma physics, Imperial College Press, London, UK, 1987.
[3] R.K. Burton, R.L. McPherron, and C.T. Russell, “An empirical relationship between interplanetary conditions and Dst,” J. Geophys.
Res., vol.80, pp.4204–4214, 1975.
[4] S. Nakano, G. Ueno, and T. Higuchi, “Merging particle filter forsequential data assimilation,” Nonlin. Process. Geophys., vol.14,
pp.395–408, 2007.
[5] T.P. O’Brien and R.L. McPherron, “An empirical phase space analysis
of ring current dynamics: Solar wind control of injection and decay,”
J. Geophys. Res., vol.105, pp.7707–7719, 2000.
[6] N.J. Gordon, D.J. Salmond, and A.F.M. Smith, “Novel approach tononlinear/non-Gaussian Bayesian state estimation,” IEE Proc. F,
vol.140, pp.107–113, 1993.
[7] G. Kitagawa, “Monte Carlo filtering and smoothing method for
non-Gaussian nonlinear state space model,” Inst. Statist. Math. Res.
Memo., no.462, 1993.
[8] G. Kitagawa, “Monte Carlo filter and smoother for non-Gaussian nonlinear state space models,” J. Comp. Graph. Statist., vol.5, pp.1–25,
1996.

Shin’ya Nakano
received the Ph.D. degree in science with a specialty in space physics
and geophysics from Kyoto University in 2004.
Currently, he is a researcher at the Instituteof Statistical Mathematics, Research Organization of Information and Systems, funded by
Japan Science and Technology Agency. He isa member of the Society of Geomagnetism andEarth, Planetary and Space Sciences, the Japan
Geoscience Union, the American GeophysicalUnion.

Tomoyuki Higuchi
completed his doctoral course in the Department of Geophysics,
in 1989 at the University of Tokyo. He holds a
Ph.D. degree in Science. Currently, he is a professor in the Department of Statistical Modelingand Vice Director General at the Institute of Statistical Mathematics, Research Organization ofInformation and Systems. His specialty is statistical analysis, in particular, Bayesian modeling for knowledge discovery. He is a member of
the Japan Statistical Society, the Society of Geomagnetism and Earth, Planetary and Space Sciences, the American Geophysical Union, and the American Statistical Association.