Transcript Document

An Updated Statistical Prediction
Model for Solar Energetic Particles:
New Approaches
Outline
Christopher Balch
Space Environment Center
SHINE Workshop
Midway, UT
31 July – 4 August 2006
Contact Info:
email: [email protected]
phone: 303-497-5693
mailing address: 325 Broadway, Boulder, CO 80305
• Abstract/Introduction
• Event database
• Current model verification measures
• Density estimation in brief
• Results: one-variable density models
• Results: two-variable density models
• Summary
Abstract
Solar Energetic Particle (SEP) events are an important type of space
weather phenomena and have real, significant impacts on human
activities and systems in space as well as communication in the polar
caps. The Space Environment Center (in Boulder Colorado) has the
responsibility to be the official source for alerts, warnings, and forecasts
of SEP's in the U.S. Although a physics based model for these
predictions would be ideal, current forecasts are limited by making use
of what solar phenomena are actually observed and available in real
time. This has led to an emphasis on development of empirical
statistical prediction methods. This paper will discuss an updated
database of SEP events and a corresponding non-SEP event database
that have been constructed for the last two solar cycles, 1986-2004.
The paper will focus on new approaches to predicting SEP probability
based on this data, including multi-variate density estimation and
discriminant analysis. The performance of these methods will be
compared with the current operational model using standard forecast
verification measures.
Introduction: SEP Effects & Affected Activities
• Impacts of SEP’s
– Human health hazard
– Spacecraft Electronics
– The Polar Ionosphere
• Affected Industries
–
–
–
–
Open
Closed
Open
Manned Spaceflight
Spacecraft Operations
HF communication
Airline operations
Introduction: SEC nowcasts, forecasts, products
• GOES real-time particle flux data
– Range of energies 0.6 – 500 MeV
– Event defined as flux of 10 p cm-2 s-1 ster-1 (PFU) at > 10 MeV
• Daily forecast
– Three day probabilities for SEP events
• Short term warnings
– Based primarily on flare observations
– Thresholds: 10 PFU  10 MeV and 1 PFU  100 MeV
– Prediction for onset time, maximum flux, time of maximum flux, and expected
event duration
• Alerts
– Real-time reports of an observed event
– Issued shortly after threshold has been attained
• Additional notices at 100, 1000, and 10000 PFU.
• Event Summaries
– Important for users to have an ‘all-clear’ indicator
Introduction: Practical Considerations
for Operational Prediction Models
• Forecasters must know the necessary conditions
for an SEP to occur
• Input parameters for the model must be valid for a
specific parent solar event
• The observations that provide the parameters must
be routinely available in real-time
• Assessment of a potential SEP must be done in
real-time – prior to arrival of the first particles
(particles may arrive before solar event is over)
• Today’s prediction models are based on empirical
and statistical analyses
Event Database
• We consider two solar cycles of proton events observed during 1986-2004.
• A review of all proton events from 1986-2004 was carried out by an examination of corrected, archived 5minute GOES integrated particle flux values.1
• A total of 165 events were identified and parameterized
• Events originating from distinct solar sources were separated (even if a pre-existing event was still in
progress)
• Associated GOES XRS events were identified by comparison of 1-minute x-ray flux time series with the
GOES particle data.1
• Associated ground-based H-alpha flare reports and radio sweep events were also carefully reviewed and
included as appropriate2
• Associated CME events from the CDAW CME catalog were identified and associated as appropriate3
• Characteristics of the proton event were noted and recorded (e.g. late, shock-driven maximum (or
secondary maximum), multiple-injection events, well-connected events, events which began when a
previous event was still in progress, suspected backside source events, slow-riser events, events with
associated erupting prominence)
• Events whose maximum was driven by the energetic storm particle (ESP) effect were flagged as such
1 Source:
NOAA/SEC CDF archive files and NOAA/NGDC GOES archives
2 Source:
NOAA/SEC event archives (collected from USAF-SEON network). Additional contributions from the NOAA/NGDC
archives are gratefully acknowledged
3
Source: We gratefully acknowledge the CDAW Data Center (The CME catalog is generated and maintained by NASA and The Catholic
University of America in cooperation with the Naval Research Laboratory. SOHO is a project of international cooperation between ESA and NASA)
Control Event Database
• We find 27/165 of the SEP events are suspected to be from a backside source. We drop these from
analysis because the radiative signals are absent or significantly affected.
• We find that 4/165 of the SEP events resulted from the ESP effect exclusively. The occurrence of these
events is related to the physics of particle trapping in interplanetary shocks – which differs from most of the
other events which have contributions of particles well ahead of the shock. These events are also dropped
from the analysis
• The remaining 134 events are found to establish the following necessary conditions for an SEP event:
– Peak X-ray Flux ≥ C2.4 (2.4 x 10-6 W m-2)
– Integrated X-ray Flux ≥ 0.008 J m-2
– Background Subtracted Integrated X-ray Flux ≥ 0.00595 J m-2
• A careful review of GOES 1-minute XRS data and XRS event archives shows that there are 4038 events
that meet the necessary conditions for the same time period (1986-2004). 134 are associated with proton
events and 3904 are not.
• 3.3% of events meeting the necessary conditions result in a proton event
• H-alpha and Radio event data was reviewed and associated as appropriate
• CME event data was also reviewed and associated as appropriate
• Patrol information for Radio observations and CME observations were included for each event (to indicate
whether observations were available at the time of the event)
• Association (or not) with an SEP event is indicated for each control event
Event Parameters
Proton Events parameters
XRS/CME event parameters
Onset
Time when 10 MeV flux begins to rise
Onset/Max/End
XRS event times (End at ½ power point)
Threshold
Time when 10 MeV flux reaches 10 PFU
Peak Flux
Maximum 1-8 Å x-ray flux (from XRS)
Maxtime
Time of 10 MeV maximum flux
OptClass
H-alpha optical class of associated flare
Endtime
Time when 10 MeV flux drops below 10 PFU
OptLocation
H-alpha location of associated flare
EventType
Characteristics of the event
Type II
Identify association of type II radio sweep (Yes or No)
P10Max
Maximum flux at 10 MeV
Type IV
Identify association of type IV radio sweep (Yes of No)
P30Max
Maximum flux at 30 MeV
CME onset
Time of onset of associated CME
P60Max
Maximum flux at 60 MeV
CME speed
Leading edge CME speed (linear fit)
P100Max
Maxmum flux at 100 MeV
Integrated XRS flux
Integral of XRS flux from onset to end
P10Fluence
Integral of 10 MeV flux from threshold to end
P30Fluence
Integral of 30 MeV flux from threshold to end
Bkgd Subtracted
integrated XRS flux
Integral of background subtracted XRS flux from onset
to end (background is taken to be the pre-event level)
P60Fluence
Integral of 60 MeV flux from threshold to end
Temperature1
Derived Temperature using ratio of two XRS channels
P100Fluence
Integral of 100 MeV flux from threshold to end
Emission Measure1
Derived Emission Measure using ratio of two XRS
channels
Associated
XRS event
Identifies associated XRS event
Associated SEP
event
Identifies associated proton event (or set to none if there
is no association)
Associated
CME event
Identifies associated CME event
1 Temperature
and Emission measure were derived using the SolarSoft library routines
GOES_CHIANTI_TEM.PRO and GOES_MEWE_TEM.PRO
Assessing Probability Forecasts: Verification Measures
• Key Measures of forecast performance: Hsu and Murphy, 1986, describe the
following attributes in assessing probability forecasts: quality, skill, reliability, and resolution
• Quality: QR  (1 / N )N ( f  o ) 2
i
i
i 1
QR is the ‘quadratic score’, where fi are the forecast probabilities and oi are the observations
(1=event, 0=no event). It is also the mean square deviation of forecasts from observations.
• Skill:
•
SS  (QR*  QR) / QR*
Measures quality of QR relative to QR*, where QR* is the quality of forecasts that use the overall
event occurrence rate as the prediction for every event (e.g. 0.033 for proton events, given the
necessary conditions)
Partitioning of QR: The quality can be separated into three components as follows:
QR  QR *  (1 / N )i 1 N i ( f i    oi  ) 2  (1 / N )i 1 N i ( oi   o ) 2
T
•
•
T
Where the forecast probabilities are divided into T distinct intervals, Ni is the number of forecasts
(and corresponding observations) for each interval, <fi> is the average forecast probability for
each interval, <oi> is the average occurrence rate of events for each interval, and ō is the overall
average event occurrence rate.
Reliability: The second term of the equation measures the correspondence between the
forecast probability and the occurrence rate and is referred to as the reliability (REL) of the
predictions. Smaller values of REL indicate more reliable forecasts and improve the overall
quality.
Resolution: The third term measures the ability of the forecasts to separate the sample into
subsamples for which the respective relative frequencies differ from the average occurrence
frequency and is referred to as the resolution (RES) of the predictions. Larger values of RES
indicate forecasts with better resolution and improve the overall quality
• An equivalent form for the equation is QR = QR* + REL - RES
Current Model Performance Attributes
• Current Operational Model based on Balch
(1999) study
• Using the database, the probability
predictions from the model were tested
• Attributes diagram (Hsu and Murphy, 1986)
shows forecast performance along with
reference lines for perfect reliability, perfect
resolution, zero skill, zero resolution
• Quality: QR=0.0246
(rms error = 0.157)
• Skill: SS = (QR*-QR)/QR*=0.234
• Reliability: REL = 0.0006
• Resolution: RES = 0.0081
• Fairly good reliability with one exception
(0.25-0.35 range)
• Problem forecast bin 0.25-0.35: occurrence
rate for a type of event(≥X7, Type II/IV, Xint≥0.895)
from 1999 study was 3/10 but new data shows
10/12
Verification Measures: Categorical Forecasts
A = number of hits
B = number of false alarms
C = number of missed events
D = number of correct nulls
• The following statistical measures provide
information about the quality of these categorical
forecasts:
POD = A / (A+C) – probability of detection
FAR = B/(A+B) – false alarm rate
PC = (A+D)/N – percent correct
HSS = (A+D–E)/(N-E) – Heidke skill score. HSS measures the
fraction of correct forecasts adjusted by E, the number of
forecasts we would expect to be correct by chance.
• For probability forecasts, probability threshold
can be treated as an independent variable
ranging from 0.0 to 1.0, hence each of the
categorical quality measures (POD, FAR, PC,
HSS) can be considered to be a function of pt
Event Observed
Event Forecast
• Suppose we define a probability threshold pt
• Suppose that we will issue a warning whenever
the forecast fi ≥ pt. However, if fi < pt then no
warning is issued
• Then the forecasts and observations can be
analyzed in terms of a 2x2 contingency table
shown to the right, where we observe that:
Yes
No
Yes
A
B
nw
No
C
D
nn
np
nc
N
Derivation of E
Probability of (event = Y) = (A+C)/N
Probability of (forecast = Y) = (A+B)/N
Probability of a chance hit
P(event=Y and forecast=Y)
(A+B)*(A+C)/N2
Probability of (event = N) = (B+D)/N
Probability of (forecast = N) = (C+D)/N
Probability of a chance correct null
P(event=N and forecast=N)
(B+D)*(C+D)/N2
Probability of chance hit or chance correct null:
[(A+B)*(A+C) + (B+D)*(C+D)]/N2
Number of correct forecasts by chance:
E = [(A+B)*(A+C) + (B+D)*(C+D)]/N
Categorical quality measures for the current model
FAR falls as the
threshold level is
increased
However, POD also
decreases with
increasing threshold
Skill score (HSS)
optimization is achieved
for range of probabilities
20-30%
At optimal point:
POD = 56% (75/134)
FAR = 55% (91/166)
PC = 96% (3888/4038)
HSS = 0.48
White – POD
Red – FAR
Green – PC
Yellow – PCH
Blue – GSS
Cyan - HSS
Raw categorical quality measures for single parameters
• Proton prediction performance
based on a single parameter:
z = log (background subtracted
integrated x-ray flux)
• Prediction Rule: issue a
warning if z  zt, where zt is a
threshold level.
• Performance is calculated as a
function of various threshold
levels
• Skill score (HSS) optimization
is achieved for zt = -0.86
(background subtracted integrated
x-ray flux = 0.14 J m-2)
• At optimal point:
POD = 59%
FAR = 59%
PC = 96%
HSS = 0.46
White – POD
Red – FAR
Green – PC
Yellow – PCH
Blue – GSS
Cyan - HSS
Results for Single Parameter Categorical Forecasts
Parameter
Optimal
Threshold
POD
FAR
PC
HSS
Bg Subtr Int
XRS flux
0.14 J m-2
59%
59%
96%
0.46
Int XRS flux
0.17 J m-2
51%
58%
96%
0.44
Log (Emission
Measure cm-3)
50.3
36%
50%
97%
0.40
Peak XRS flux
X1.9
37%
58%
96%
0.37
CME speed
1532 km/s
40%
60%
92%
0.36
Type II & IV
Both occur
57%
77%
92%
0.30
Derived T
22 x 106 K
34%
82%
93%
0.20
Longitude
≥W79
9%
92%
92%
0.04
Single Variable Density Estimation
• Density Estimation (Silverman, 1998) is a method of
deducing continuous probability density functions from
observed data
• We consider the distribution of one of the parameters in
our data set, for example integrated x-ray flux
• We expect to have different probability distributions for
this parameter, depending on whether it is from the set of
proton associated events, or from the set of events not
associated with proton events
• Once smooth distributions are found for these two cases,
we can deduce a continuous function for the probability
for an event as a function of the parameter
Description of the Method
• Let {Xi} be the set of n observed values of the parameter, i=1, 2,
…, n
• The probability density function fˆ is defined on the continuous
domain of values x for this parameter, such that:
n
1
 x  Xi 
fˆ x  
K

 
nh i 1  h 
Where we use the normal probability density function for the
z2 2
kernel function K:
e
K ( z) 
2
and h is a smoothing parameter
• Geometrically, each observed value contributes a ‘bump’ to the
overall density estimate, so that the resulting function is a
smooth, continuous probability density for the parameter Xi
Example
Density estimate for log of background
subtracted integrated x-ray flux for proton
events using the normal probability
density with a window width of 0.15
Density estimate for log of background
subtracted integrated x-ray flux for events
not associated with proton events (control
events)
Probability model for proton events, using the density functions
Probability model for proton events,
using the density functions
Probability = np*f1/(np*f1 +nc*f2),
np - number of proton events
f1 - density estimate - proton associated events
nc - number of control events
f2 - density estimate for the control events
The probability is set to constant
once maximum probability is reached
Density Model metrics:
• Accuracy = 0.0235
• Rms error = 0.153
• Skill = 0.269
• Reliability=2.59 x 10-4
• Resolution=8.94 x 10-3
For comparison - metrics
for operational model:
• Accuracy = 0.0246
• Rms error = 0.157
• Skill = 0.234
Density analysis for one parameter (background
subtracted integrated x-ray flux) gives better
accuracy, skill, reliability, and resolution !
However, the result is preliminary since the training
set is being used for the verification calculation
• Reliability=6.0 x 10-4
• Resolution=8.1 x 10-3
White – POD
Red – FAR
Green – PC
Yellow – PCH
Blue – GSS
Cyan - HSS
Performance of
categorical forecasts
using probability
thresholds
one parameter density
model
Parameter is background
subtracted integrated xray flux
Skill score optimization
is at probability
threshold of 0.18
At optimal point:
POD = 59% (79/134)
FAR = 59% (113/192)
PC = 96% (3773/4020)
HSS = 0.46
Two Parameter Density Estimation
• The method can be applied using more than one variable as a prediction
inputs
• In this case we consider a set of n observation vectors Xi which consist two
components
• Each component represents one of the observed variables (e.g. 1st component
could be integrated x-ray flux, 2nd component could be CME speed)
• The density estimate in this two-dimensional parameter space is defined to be:
1
fˆ x   2
nh
 x  Xi 
K



h


i 1
n
where we use the standard, multivariate normal density function:

exp  12 z T z
K (z ) 
2

• In order to use a scalar smoothing parameter, h, it is necessary to
normalize the observation vectors – we rescale all of the parameters so
that their range is restricted to the interval [0,1]
2D Density Maps for Background
Subtracted Integrated X-ray flux
combined with H-alpha flare longitude
Red – proton event associated
Green – control event associated
Corresponding Probability
Map for Proton Events
Performance Statistics for this 2D probability model
(Log Background Subtracted Integrated X-ray Flux & Longitude)
QR = 0.0263
RMS = 0.162
Skill = 0.308
REL = 7.60 x 10-4
RES = 1.34 x 10-2
Performance for Categorical Forecasts – 2D density model
White – POD
Red – FAR
Green – PC
Yellow – PCH
Blue – GSS
Cyan - HSS
At optimal point, threshold probability = 21%,
POD = 53.7%, FAR = 41.5%, PC = 96.6%, HSS = 0.542
2D Density Maps for CME speed
combined with Emission Measure
Red – proton event associated
Green – control event associated
Corresponding Probability
Map for Proton Events
Performance Statistics for this 2D probability model
(CME speed and Emission Measure)
QR = 0.0423
RMS = 0.206
Skill = 0.316
REL = 1.86 x 10-3
RES = 2.14 x 10-2
Performance for Categorical Forecasts – 2D density model
White – POD
Red – FAR
Green – PC
Yellow – PCH
Blue – GSS
Cyan - HSS
At optimal point, threshold probability = 24%,
POD = 51.6%, FAR = 41.8%, PC = 94.3%, HSS = 0.517
Top Performers with respect to HSS
2D probability density models
Parameter1
Parameter2
BgSub Int Xray Flux
longitude
Int Xray Flux
QR
RMS
SS
REL
RES
0.0263
0.162
0.325
7.60E-04
1.34E-02
longitude
0.0269
0.164
0.308
1.19E-03
EM (Mewe)
CME speed
0.0423
0.206
0.316
CME speed
longitude
0.0544
0.233
EM (Chianti)
CME speed
0.0442
Log Max Xray Flux
CME speed
0.0443
HSS
Thresh
POD
0.543
0.21
53.7
41.5
96.6
1.32E-02
0.524
0.18
53.7
45.0
96.3
1.86E-03
2.14E-02
0.517
0.24
51.6
41.8
94.3
0.313
3.87E-03
2.86E-02
0.495
0.25
43.6
32.5
93.3
0.210
0.286
1.37E-03
1.91E-02
0.490
0.28
45.2
39.1
94.4
0.210
0.285
1.45E-03
1.91E-02
0.484
0.27
45.2
40.4
94.3
All of these pairs of parameters result in a prediction model
the HSS better than the existing operational model
FAR
PC
Summary
• Current practical SEP prediction models are limited by available observations and typically rely on
empirical/statistical relationships
• A comprehensive event database for 1986-2004 has been compiled for SEC proton events and
corresponding control events, where the control events are selected based on the necessary
conditions identified from the proton events
• The current operational model performance was evaluated using the data – standard verification
measures were derived for accuracy, reliability, resolution, and relative skill
• Categorical performance measures were also calculated as a function of probability thresholds
• Categorical performance based on optimal thresholds for input parameters were derived, but none
performed better than the operational model
• Density Estimation techniques were applied to selected parameters: the performance of the
probability forecasts were improved for one of the parameters relative to the operational model. The
categorical verification measures did not change significantly relative to the simple parameter
threshold method
• 2D densities and corresponding probability maps were derived. The following pairs of parameters
were found to have improved performance measures relative to the operational model:
Integrated X-ray Flux and Longitude, Emission Measure and CME speed, CME speed and
longitude, Peak X-ray flux and CME speed.
• The method shows promise – but a more rigorous approach needs to used to derive the verification
statistics
• We expect that adding more dimensions/parameters to the analysis will also improve the quality of
these important space weather predictions