EGU2020-10385_presentation LSTM algorithm

LSTM algorithm

1. Using AutoRegressive Integrated Moving Average
and Gaussian Processes with LSTM neural
networks to predict discrete geomagnetic signals
Authors: Laurentiu Asimopolos 1, Alexandru Stanciu 2,
Natalia-Silvia Asimopolos 1, Bogdan Balea 1,
Andreea Dinu 2, Adrian-Aristide Asimopolos 3
1 Geological Institute of Romania, Surlari Geomagnetic Observatory, BUCHAREST,
Romania ([email protected])
2 National Institute for Research and Development in Informatics BUCHAREST, Romania
([email protected])
3 Polytechnic University of Bucharest, Faculty of Transportation, BUCHAREST, Romania
([email protected])
Session EMRP2.5 –Electromagnetic induction in geophysics

2. Inthispaper,wepresenttheresultsobtainedforthegeomagneticdataacquired
attheSurlariObservatory,locatedabout30KmNorthofBucharest-Romania.
Theobservatorydatabasecontainsrecordsfromthelastsevensolarcycles,with
differentsamplingrates.
WeusedAR,MA,ARMAandARIMA(AutoRegressiveIntegratedMoving
Average)typemodelsfortimeseriesforecastingandphenomenological
extrapolation.ARIMAmodelisageneralizationofanautoregressivemoving
average(ARMA)model,fittedtotimeseriesdatatopredictfuturepointsinthe
series.
WemadespectralanalysisusingFourierTransform,thatgivesusarelevant
pictureofthefrequencyspectrumofthesignalcomponent,butwithoutlocatingit
intime,whilethewaveletanalysisprovidesuswithinformationregardingthetime
ofoccurrenceofthesefrequencies.
Waveletallowslocalanalysisofmagneticfieldcomponentsthroughvariable
frequencywindows.Windowswithlongertimeintervalsallowustoextractlow-
frequencyinformation,medium-sizedintervalsofdifferentsizesleadtomedium-
frequencyinformationextraction,andverynarrowwindowshighlightthehigh-
frequenciesordetailsoftheanalysedsignals.

3. Weextendthestudyofgeomagneticdataanalysisandpredictivemodellingby
implementingaLongShort-TermMemory(LSTM)recurrentneuralnetworkthatiscapable
ofmodellinglong-termdependenciesandissuitablefortimeseriesforecasting.Thismethod
includesaGaussianprocess(GP)modelinordertoobtainprobabilisticforecastsbasedonthe
LSTMoutputs.
TheevaluationoftheproposedhybridmodelisconductedusingtheReceiverOperating
Characteristic(ROC)Curvethatprovidesaprobabilisticforecastofgeomagneticstorm
events.
Inaddition,reliabilitydiagramsareprovidedinordertosupporttheanalysisofthe
probabilisticforecastingmodels.
Theimplementationofthesolutionforpredictingcertaingeomagneticparametersis
implementedintheMATLABlanguage,usingtheToolboxDeepLearningToolbox,which
providesaframeworkforthedesignandimplementationofdeeplearningmodels.
Also,inadditiontousingtheMATLABenvironment,thesolutioncanbeaccessed,
modified,orimprovedintheJupyterNotebookcomputingenvironment.

4. REMARKSABOUTSURLARIGEOMAGNETIC OBSERVATORY
GeomagneticfieldstudyinRomanianstationshasstartedwithirregularmeasurements,lateXIX
-th
century.In1943,thefoundationofSNGOmarksthebeginningofanewerainthesystematicstudyof
geomagneticfieldbyacontinuousregistrationofitsvariationsandbycarryingoutstandardabsolute
measurementsinafundamentalstation.Observatorylocationwasthoroughlyestablished,sothatmeetsthe
geomorphologicalandtechnologicalcriteria.SNGOislocatedinCăldăruşani-Surlariastro-geodeticpolygon,
inanareawithoutmagneticfieldanomaliesorsignificantlocalheterogeneityofelectricalconductivityinthe
basementandasufficientdistancefrommajorindustrialsourcesofdisturbance.
TheObservatorycoversanareaof3.56hectaresandcomprisessevenbuildingsandunderground
laboratoryforgeomagneticsensors.Insideareinstalledspecializedequipment’stomultiparametricmonitoring
offieldsofearth.Thesebuildingsweremadeduring1943-1969.Theundergroundlaboratoryandthemain
buildingswererenovatedandmodernizedin2006-2008.Thedesignofthespecialgeomagneticrecording
laboratorieswasmadeaftersomewellverified.Insidelaboratorieswerebuilt18speciallydesignedpillars
embeddeddeepintheground,whicharemountedhigh-resolutionsensors.

5. General view in SNGO
Absolute measurement laboratory in SNGO
The first magnetogram in October 16, 1943
in SNGO

7. ShortdescriptionofnetworkandinfrastructurefeaturesattheSNGO:
• NATlocalnetworkaccessimplementednativelybytheCisco-880router
• VPNRemoteaccessserverintheObserver'slocalnetworkforuserslocated
inotherlocationsimplementedthroughCiscoIOS12.4T
• DHCPautomaticconfigurationofnetworksettingsimplementedwithISC
dhcpd3.1.1
• DNSaddressresolutionbynameandviceversaimplementedwithISCBIND
9
• IntegratedLDAPauthenticationofusersinthelocalnetworkimplemented
withopenldap22.4.12
• NTPclocksynchronizationserviceforallcomputersconnectedtothelocal
networkimplementedwithntpd4.2.4
• DatabaseserverforstoringdatafromallSurlariGeomagneticObservatory
acquisitionsystemsimplementedwithMySql5.0.67
• Webserverhoststheobservatory'swebsitedeployedwithApache22.2.10
2)Automatically,periodically(e.g.,averagedover5secondsor60seconds)tothe
Webserver.
•FTPserverfiletransferservicesforIntuitiveusersusingIntranetvsftpd2.0.7
•FileserverserverforgeneralusefilestoragewithSamba3.2.4

8. Forimplementation,theseserviceshavebeengroupedinto4categories,takinginto
accounttheirspecificity,theneedforhardwareresourcestorun,andcriticaldependency
betweenthem.Forimplementation,thesolutionwaschosenaseachgrouptorunona
separateserver.Cisco-880router-runsNATandVPNserverservices(alreadyimplemented
andactive).
netserv*-serverforcommonnetworkservices:DHCP,DNS,LDAP
dbserv*-databaseserver
filer*-fileserver
webserv*-WebserverandFTP
Toreducecostsandreducethenumberofcomputersinstalled,allarevirtualmachines
runningunderVMwareServer2onasinglephysicalserverforahighreliability(RAID1for
dualharddrives,dualNICmountedinfail-overarchitecture,UPSwithmonitoringsoftware).
Modernizationofthenetworkinfrastructurehasbeencarriedoutsothatitcanprovideall
theservicesnecessaryfortheexploitation.
Description of each unit within the data server
Schematic overview of the database
acquisition and management program.

9. Thesoftwarepackageconsistsofthefollowingcomponents:
-SQLServer.VIonthedatabaseserverisdesignedtotransmitmeasuredvaluesin
twoways:
1)uponrequest,toapplicationsinstalledoncomputersofdifferentusers;users
canloginwithausernameandpasswordtoaccessthedatastoredonthedatabase
server;
2)Automatically,periodically(e.g.,averagedover5secondsor60seconds)tothe
Webserver.
ClientVI,withmultipleinstallationsondifferentusers'computers,hastheroleof
allowingthemtocreateontheirowncomputerstextfilescontainingmeasured
values​​retrievedfromthedatabaseserver(figure3and4);
Thisallowstheusertoselectthestartandendmomentsofthetimeperiodforwhich
theSQLServerVI,applicationwillberequested,copythemeasuredvalues;
Theusercanspecify,whereappropriate,thepathandfilenameinwhichthe
measuredmeasurementvalues​​transmittedbytheSQLServerVIapplicationare
saved;
-WebServerVIisinstalledontheWebserverandreceivesthemeasuredvalues,
periodicallytransmittedbySQLServerVIandwillsavethemintextfiles;Thefiles
aredaily,thefilenameisDD-MM-YY,andthefileswithinamonthwillbestoredin
thesamefolder,themonthlyfoldernamebeingMM-YY.

10. -Java Client VI, is installed on the Web server and is called by the actions of the
various web page visitors and allows viewing of measured values ​​in graphical form.
These software applications respect the SNGO internal network and public network
connections, thus centralizing data acquisition and online transmission on Web
Server.
Standard processing of geomagnetic data according to INTERMAGNET
requirements includes calculating the average at one minute of Geomagnetic
components.
Values and obtaining provisional data files, establishing the base level of records as
well as adopting the baseline level and making definitive data.
This infrastructure and organizing of database from SNGO allow us to monitoring
and make all processing procedures of data according to geomagnetic observatories
reequipments.

11. THE GEOMAGNETIC EQUIPMENT USED IN SURLARI OBSERVATORY
Animportantoperationinthefunctioningofmagnetometersisthecalibrationforestablishingthe
sensitivityofeachcomponentrecords(nT/mm).ThisisdonebymeansofDCpoweredcoilsplacedinthe
directionsperpendiculartotheaxissensors.TheDCpowerhaveaknownintensity.Thetwosystemshave
operatedcontinuouslyuntil2003,producinganaloguerecordsonphotographicpaper.
Itisalsoveryimportanttoobtainbaselevelofrecordswithabsolutemeasurementsmade​​withtheodolite
MattingWeissenberg,groundinductorandoscillationsboxinafirstphaseandthentheodoliteMatting
WeissenbergandquartzhorizontalmagnetometerforHcomponentQHM.Inadditiontothesedevicesin1968
wasbroughtVarianprotonprecessionmagnetometermeasuringthescalarvalueofthetotalfield.
TriaxialmagnetometerMAG-03MCallowssimultaneousrecordingofthecomponentsHx(North
direction),Hy(ontheEast)andHz(verticaldirectiondown)ofthegeomagneticfield.Thesensorsofthis
magnetometeraremagneticinductivetypeandaremadeofcoilswithalargenumberofturnsandamagnetic
corewithhighpermeability(Permalloy).
Characteristicresponse(floortype)ofthismagneticsensor,showsaverygoodfunctionalityforabandof
frequenciesbetween1Hzand2000Hz

12. MAG03DAMloggerhastwoconnectorsforanaloginputsignalfromtwofluxgatemagnetometersandaRS-
232connectorwith25pinsforoutputsignal.LoggeriscontrolledbysoftwaredevelopedinFORTRANthatallows
selectingthenumberofchannelstoberecorded(from1-6),thechoiceofsamplingrate(betweenonesecondand10
seconds)andthemeasuringrange.Thesamplingratereferstothefrequencyofpurchaseandnotthestorage,itis10
timeslower.
MAG-01hDIFluxgateMagnetometerproducedbyBartingtonInstrumentsLtd.,England,withTHEO010b
nonmagnetictheodolite,measuringdeclinationandinclinationofthegeomagneticfieldinabsoluteterms.These
measurementsareusedtoestablishthebaselevelofpermanentrecords.
ThistypeofmagnetometercanbeusedmountedonaWildT1theodolitenonmagnetic.Itcanmeasureboth
declinationandinclinationofthegeomagneticfieldwithanaccuracyof±onesecondsexagesimaldegree.
GeometricsG-856protonprecessionmagnetometerareusedforabsolutemeasurementsoftotalmagneticfield.
FGEvectormagnetometerwasbuiltbyDanishMeteorologicalInstituteusingthreecommercialfluxgate
sensors,mountedinablockofmarble12x12x12cm
3
throughthequartztubeisoffsetcoilsthatensuremaximum
stabilityordriftto3nT/year.Variationwithtemperatureofrecordedvaluesduesensors​​isbelow0.2nT/°Candofthe
electronics,as0.1nT/°C.
GSM90Overhauserprotonmagnetometerisascalarmagnetometerdesignedformagneticobservatories
andotherapplications(Volcanology),wherestabilityandaccuracyareabsolutelynecessary.
TorsionphotoelectricmagnetometerPSM

13. ARIMAmodelsarethemostgeneralclassofmodelsforforecastingatimeserieswhich
canbemadetobe“stationary”bydifferencing,perhapsinconjunctionwithnonlinear
transformations.Arandomvariablethatisatimeseriesisstationaryifitsstatisticalproperties
areallconstantovertime.Astationaryserieshasnotrend.Thisconditionmeansthatits
autocorrelations(correlationswithitsownpriordeviationsfromthemean)remainconstant
overtime,orequivalently,thatitspowerspectrumremainsconstantovertime.Arandom
variableofthisformcanbeviewedasacombinationofsignalandnoise,andthesignalcould
beapatternoffastorslowmeanreversion,orsinusoidaloscillation,orrapidalternationin
sign,anditcouldalsohaveaperiodicalcomponent.AnARIMAmodelcanbeviewedasa
“filter”thattriestoseparatethesignalfromthenoise,andthesignalisthenextrapolatedinto
thefuturetoobtainforecasts.
ARIMA(p,d,q)forecastingequationforastationarytimeseriesisalinear(i.e.,regression-
type)equationinwhichthepredictorsconsistoflagsofthedependentvariableand/orlagsof
theforecasterrors.
TheacronymARIMAstandsforAuto-RegressiveIntegratedMovingAverage.Lags
ofthestationarizedseriesintheforecastingequationarecalled"autoregressive"terms,lagsof
theforecasterrorsarecalled"movingaverage"terms,andatimeserieswhichneedstobe
differencedtobemadestationaryissaidtobean"integrated"versionofastationaryseries.
Random-walkandrandom-trendmodels,autoregressivemodels,andexponentialsmoothing
modelsareallspecialcasesofARIMAmodels.
Anon-periodicalARIMAmodelisclassifiedasan"ARIMA(p,d,q)"model,where:
•pisthenumberofautoregressiveterms,
•disthenumberofnon-periodicaldifferencesneededforstationarity,and
•qisthenumberoflaggedforecasterrorsinthepredictionequation.

14. The forecasting equation is constructed as follows. First, let ydenote the d
th
difference of
Y, which means:
If d=0: y
t= Y
t
If d=1: y
t= Y
t-Y
t-1
If d=2: y
t= (Y
t-Y
t-1) -(Y
t-1-Y
t-2) = Y
t-2Y
t-1+ Y
t-2
The second difference of Y(the d=2 case) is not the difference from 2 periods ago. Rather, it
is the first-difference-of-the-first difference, which is the discrete analogue of a second
derivative, i.e., the local acceleration of the series rather than its local trend.
In terms of y, the general forecasting equation is:
ŷ
t= μ + ϕ
1y
t-1+…+ ϕ
py
t-p-θ
1e
t-1-…-θ
qe
t-q
Here the moving average parameters (θ’s) are defined so that their signs are negative in the
equation, following the convention introduced by Box and Jenkins.

15. Matlab code:
STSFMs_ARIMA.AR = [];
STSFMs_ARIMA.MA = [];
Lags = num2cell(nan(1,STSFMs_ARIMA.Seasonality));
STSFMs_ARIMA.SAR = Lags;
STSFMs_ARIMA.SMA = Lags;
else
STSFMs_ARIMA = arima('SMALags',12,'D',1,'SARLags',12,...
'Seasonality',12);
D1 = LagOp({1 -1},'Lags',[0,1]);
D12 = LagOp({1 -1},'Lags',[0,12]);
D = D1*D12;
dFinal_test = filter(D,Final_test);
STSFMs_ARIMA_1 = arima('SMALags', 12, 'D',1,'SARLags',12,'MALags',1,'ARLags',1,...
'Seasonality',12);
STSFMs_ARIMA_Minitab = arima('SMALags', 12,'SMA',0.9852,
'D',1,'SARLags',12,'SAR',0.1264,'MALags',1,'MA',0.0500,'ARLags',1,'AR',0.9688,...
'Seasonality',12,'Constant',-0.0002256,'Variance',2);
max_order = 3;
[aic_matrix,bic_matrix] = ARMA_model_order_tuning(Final_test,max_order);
end

16. LongShortTermMemory(LSTM)networksareatypeofrecurrentneuralnetwork
thatiscapableoflong-termdependencies.
Conceptually,arecurringLSTMunittriesto"remember"allthepastknowledgeabout
whichthenetworkisseensofarandto"forget"theirrelevantdata.Thisisdoneby
introducingdifferentlayerswithactivationfunctionscalled"gates"fordifferentpurposes.
EachrecurrentLSTMunitalsomaintainsavectorcalledan"internalcellstate"that
conceptuallydescribestheinformationthatwaschosentoberetainedbytheprevious
recurrentLSTMunit.AnLSTMnetworkcomprisesfourdifferentgatesfordifferent
purposes,asdescribedbelow:
•Forgetgate:determinestheextenttowhichpreviousdatacanbeforgotten.
•Inputgate:determinestheinformationtobewrittenintheinternalcellstate.
•Inputmodulationgate:itisconsideredasapartoftheinputgateandisusedto
modulatetheinformationthattheinputgatewillwriteontheinternalcellstateinternally
addingnonlinearitytotheinformationandnormalizingtheinformation.Thisisdoneto
reducelearningtime,ensuringfasterconvergence.Althoughtheactionsofthisgateareless
importantthantheothersandareoftentreatedasaconceptthatoffersfinesse,itisgood
practicetoincludethisinthestructureoftheLSTMunit.
•Outputgate:determineswhichoutput(nexthiddenlayer)isgeneratedfromthe
internalstateoftheLSTMunit.

17. The architecture of an
LSTM cell
The implementation of the solution for predicting certain geomagnetic parameters is
implemented in the MATLAB language, using the Deep Learning Toolbox. It provides a
framework for the design and implementation of Deep Learning networks, with the
development of convolutional neural networks (ConvNets, CNN) and recurrent LSTM
networks for the classification or regression of image, time series and text data.
Also, in addition to using the MATLAB environment, the solution can be accessed,
modified or improved in the Jupyter Notebook computing environment.

18. Stages of implementing the solution
1. The data can be read as follows:
format short
omni_data= readtable('fisier.txt');
omni_data.Properties.VariableNames= {'Var1' 'Var2'
'Var3' 'Var4' 'Var5' 'Var6' 'Var7' 'Var8'};
data = table2cell(omni_data(:,8));
data = [data{:}];
2. Training data partitioning -70% and testing -30%:
numTimeStepsTrain= floor(0.7*numel(data));
dataTrain= data(1:numTimeStepsTrain+1);
dataTest= data(numTimeStepsTrain+1:end);
3. Data normalization
mu=mean(dataTrain);
sig = std(dataTrain);
dataTrainStandardized= (dataTrain-mu) / sig;

19. 4. Preparation of training data
XTrain= dataTrainStandardized(1:end-1);
YTrain= dataTrainStandardized(2:end);
5. Defining the Long Short Term Memory network
numFeatures= 1;
numResponses= 1;
numHiddenUnits= 150;
layers = [ ...
sequenceInputLayer(numFeatures)
lstmLayer(numHiddenUnits)
fullyConnectedLayer(numResponses)
regressionLayer];
options = trainingOptions('adam', ...
'MaxEpochs',150, ...
'GradientThreshold',1, ...
'InitialLearnRate',0.005, ...
'LearnRateSchedule','piecewise', ...
'LearnRateDropPeriod',125, ...
'LearnRateDropFactor',0.2, ...
'Verbose',0, ...
'Plots','training-progress');

20. 6. Training of the LSTM network
net = trainNetwork(XTrain,YTrain,layers,options);
7. Normalization of test data
dataTestStandardized= (dataTest-mu) / sig;
XTest= dataTestStandardized(1:end-1);
8. Effective prediction of the variable selected for analysis, based on the observed
data
net = resetState(net);
net = predictAndUpdateState(net,XTrain);
YPred= [];
numTimeStepsTest= numel(XTest);
for i= 1:numTimeStepsTest
[net,YPred(:,i)] =
predictAndUpdateState(net,XTest(:,i),'ExecutionEnvironment','gpu');
end

21. 9. Calculation of the square error
%unstandardizedata
YPred= sig*YPred+ mu;
YTest= dataTest(2:end);
rmse= sqrt(mean((YPred-YTest).^2))
10. Display of observed data (YTest) and forecasted data (YPred)
figure
plot(dataTrain(1:end-1))
hold on
idx=
numTimeStepsTrain:(numTimeStepsTrain+numTimeStepsTest);
plot(idx,[data(numTimeStepsTrain) YPred],'.-')
hold off
ylabel("Values")
title("Forecast")
legend(["Observed" "Forecast"])

22. Fortheexperimentalpart,weusedanLSTMmodeltopredicttheDstvalues1h,2h,
and3hahead,andaGaussianProcessclassifiertopredictthelevelofthegeomagnetic
storms.
ThetrainingdatasethasfourmonthsofhourlyDstvalues,startingfromMarch1st,
2001untilJune28th,2001.Thisperiodincludesthegeomagneticstormeventthat
occurredbetweenMarch31st–1stApril,asshowninnextfigure:
Training data,
1
st
March 2001 –28
th
June 2001

23. The Gaussian Process multiclass classifier is used for geomagnetic storm prediction
depending on the level of activity. For instance, for Dst < -250 nT there is super
geomagnetic storm, for -50 nT < Dst < -250 nT there is an intense storm, and for Dst > -
50 nT a moderate storm. For a 3h ahead prediction, we obtained a 0.98 ROC AUC, as
presented in next figure:
ROC curve of the
Gaussian Process classifier

24. WehavealsoanalyzedthegeomagneticstormthatoccurredbetweenDecember
14thandDecember16th,2006,aspresentedinnextfigure:
Geomagnetic storm
between 14
th
December 2006 –16
th
December 2006

25. Geomagnetic activity
predicted by the LSTM
model (3h ahead)
The3haheadestimationoftheDstispresentedinnextfigure:

26. ROC curve for the Gaussian
Process classifier (3h ahead)
Reliability diagram of the GP
classifier (3h ahead)
TheevaluationoftheGaussianProcessclassifierispresentedinthenexttwofigures:
TheLSTMmodelhadanRMSEof12.12andacorrelationcoefficientof0.83.

27. CONCLUSIONS
Forecastingtimeseriesdataisanimportantsubjectinmanydomains,inclusivein
geomagnetism.Traditionally,thereareseveraltechniquestoeffectivelyforecastthenext
lagoftimeseriesdatasuchasunivariateAutoRegressive(AR),univariateMovingAverage
(MA),SimpleExponentialSmoothing(SES),andmorenotablyAutoRegressiveIntegrated
MovingAverage(ARIMA)withitsmanyvariations.Inparticular,ARIMAmodelhas
demonstrateditsoutperformanceinprecisionandaccuracyofpredictingthenextlagsof
timeseries.
LongShortTermMemory(LSTM)isarecurrentneuralnetwork(RNN)thatenables
supportfortimeseriesandsequencedatainanetwork.LSTMperformsadditive
interactions,whichcanhelpimprovegradientflowoverlongsequencesduringtraining.
LSTMarebestsuitedforlearninglong-termdependencies.

28. ACKNOWLEDGEMENT
WeaddressourthanksbytheRomanianMinistryofEducationandResearchfor
financingoftheprojects:
•Therealizationof3Dgeological/geophysicalmodelsforthecharacterization
ofsomeareasofeconomicandscientificinterestinRomania,withContractno.
49N/2019.
•Institutionalcapacitiesandservicesforresearch,monitoringandforecasting
ofrisksinextra-atmosphericspace,acronymSAFESPACE,project
Nr.16PCCDI/2018,withinPNCDIII.

29. Thankyoufor your
attention!