基于经验模态分解的高能电子通量小时预报模型建立方法与流程

本发明涉及一种空间高能电子通量的预报方法，特别涉及一种基于经验模态分解的高能电子通量小时预报模型建立方法。

背景技术

在磁暴恢复相期间，导致卫星无法正常运行或完全损坏。地球同步轨道位于外辐射带区域，该区域分布了大量高能带电粒子(相对论电子)。同时，数百颗同步轨道卫星该区域内运行(参考文献1)。当发生大的磁暴后，高能电子通量会在短时间下降，随后会升高3-4个量级，这些高能粒子具有很高的能量，大量高能电子通量从磁层的外辐射带渗透到地球同步轨道(geo)，在这其中>2mev的高能电子能够穿透卫星表面并聚积在材料内部，引起材料的深层充放电，严重威胁卫星安全与稳定工作(参考文献2)。根据统计,1992年3月到1994年4月,因高能带电粒子的聚积引发的geo卫星故障共有50多次(参考文献3)。因此，对磁暴期间高能电子通量的预报,能够提前对卫星采取必要措施降低高能电子对卫星的危害，具有重要的科学价值与应用价值。

目前，高能电子通量的预报主要采用统计的和人工智能的方法等。在统计模型当中，以线性模型与非线性模型为主。早期的研究当中，paulikas等(参考文献4)首先发现高能电子通量与太阳风速度有着良好的相关性。在此基础上，baker等通过统计的方法和线性预测滤波技术(简称plf)发现太阳风速度与geo的高能电子通量具有最直接的联系(参考文献5)。该方法以历史探测到的太阳风速度作为输入，预报后一天的高能电子通量，结果表明在太阳风速度达到最大值时，电子通量输出的脉冲响应时间系数在1-2天左右达到最大值。美国nasa的空间天气预报中心在此基础上发展了相对论电子通量的预报模型(简称refm)。refm模式利用太阳风参数作为输入，输出基本预报值，然后对预报值进行修正，为预测-校正法，输出后面1-3天的高能电子通量值，第一天的预报效率约为0.71，后面两天预测值很差，原因是外辐射带在磁扰动期处于快速变化过程中(参考文献5)。由于磁层对太阳风的响应是非线性的，单用线性滤波技术预报有局限性(参考文献6)。rigler(参考文献7)在2004年利用卡尔曼滤波改进线性滤波模型，使滤波系数可以动态地改变，结果表明，预报效率得到了较大的提高，该方法被运用到了许多模型[sakguchietal.2012；hetian,2013；](参考文献8、9)。

除了线性滤波模型，非线性滤波模型也得到了广泛发展。非线性滤波模型主要基于高能电子突然增加的加速机制出发。目前，高能电子的加速机制有两种观点：第一，径向扩散机制；第二，波粒相互作用。基于经向扩散机制，li(参考文献10、11)[2001；2004]提出了径向扩散模型(简称rdf模型)，该模型用太阳风参数和行星际磁场分量作为输入，预报1-2天后geo上>2mev的高能电子通量。预报效率为0.64。对太阳活动高年的预报精度不理想。turner在2008年发展了low-e模型，该模型以前一日与当日的低能量电子通量和相对论电子通量作为输入，预报后一天的电子通量值，平均预报率在0.73(参考文献12)。基于波粒相互作用机制，alexanderet.al.[2017]集合了径向扩散机制与波粒相互作用机制(参考文献13、14)，运用多元线性回归的方法建立了同步轨道区>2mev电子的日积分预报模式，输入参数为太阳风速度、密度、动压，ulf波强度(pc4-5)，>600kev电子通量(种子电子)和行星际电场值，提前一天的预测标准差约为0.58，该模型的特点是电子通量极值预报值提前于测量值。

由于高能电子通量变化与各参数之间的关系是非线性的，高能电子通量值的变化也是非稳定的，所以输入参数与输出值之间难以用函数关系描述。神经网络方法有很好的学习能力，是解决这些非线性关系的较好方法。fukataetal.(参考文献15)[2002]和ling(参考文献16)[2010]建立神经网络模型预报高能电子通量，fukata的模型预报效率约为0.6，ling的多层反馈模型预报效率要好于前者，预报效率约为0.7左右，输入参数为地磁扰动指数，忽略了太阳风参数。此外，径向基函数、支持向量机等智能算法也被运用于高能电子通量的预报中[xue&ye,2013；guoetal,2013等](参考文献17、18)。

以上的模型在预报高能电子通量方面取得了巨大的成功，但精度还有提升的空间。由于高能电子通量在磁暴期间变化非常剧烈，几十分钟之内就会快速地上升，一天之内可从10上升到10⁵(electron/cm²srs)(参考文献8)。因此，从高能电子通量随时间变化的趋势来看，其实质是一个非平稳非线性的时间序列，且非平稳性的特点尤为突出。以往的模型利用统计的方法解决了非线性给预报带来的影响(参考文献19、20)，但非平稳性的影响没有被重视，给预报带来困难。

涉及到的参考文献如下：

[1]刘帅，李智，基于系统辨识技术的地球辐射带高能电子通量预报模型研究综述，装备学院学报[j]，2015,26(4),82-88.

[2]wrenn,g.l.,rodgersd.j.,k.a.rydexl,asolarcycleofspacecraftanomaliesduetointernalcharging,ann.geophys.,2002,22,953-956

[3]何甜,刘四清,薛炳森等，利用地磁脉动预报地球同步轨道相对论电子通量方法的研究，地球物理学报，2009,52(10)，2419-2427.

[4]paulikasga,blakejb.effectsofsolarthewindonmagnetosphericdynamics:energeticelectronsatthesynchronousorbit[c]//quantitativemodelingofmagnetosphericprocesses,americangeophysical,1979；180-202

[5]bakerd.n.,mcpherronr.l.,etal.,linearpredictionfilteranalysisofrelativisticelectronpropertiesat6.6re,j.geophy.res.,1990,95,15133-15140.

[6]lisheng,huangwenggeng,dynamicpredictionmodelofrelativisticelectrondifferentialfluxesatthegeosynchronousorbit.chin.j.2017,37(3):298-311.

[7]rigler,ej,baker,d.n.,weigel,r.s.,vassiliadis,d.,&klimas,a.j.(2004).adaptivelinerpredictionofradiatuionbeltelectronsusingthekalmanfilterspaceweathertheinternationaljournalofresearch&applications2(3)2004

[8]sakguichk.,y.miyoshi,s.satio.,etal.,elativeelectronfluxforecastatgeostationaryorbitusingkalamanfilterbasedonultivariateautoregressivemodelspaceweather,2013,11,79-89,doi:10.1002/swe.20020.

[9]何甜,刘四清,沈华，龚建村，使用地磁脉动参数定量预报地球同步轨道相对论电子通量的建模研究，空间科学学报，2013，33(1)，20-27.

[10]lix.,m.temerin,etal.,quantitativepredietionofradiationbeltelectronsatgeostationaryorbitbasedonsolarwindmeasurements,geophys.res.lett.,2001,28,1887-1890.

[11]lix.,variationsof0.7-60.mevelectronsatgeosynchronousorbitarafunctionofsolarwind,spaceweather,2004,2,s03006,doi:10.1029/2003sw0000017.

[12]turnerd.l.,lix.l.,quantitativeforecastofrelativisticelectronfluxatgeosynchronousorbitbasedonlow-energyelectronflux,spaceweather,2008,6,s05005,doi:10.1029/2007sw000354.

[13]alexanderpotapov,l.ryzhakova,b.tsegmed,anewapproachtopredictandestimateenhancementsof“killer”electronfluxatgeosynchronousorbit.actaastronautica,2016,126,47～51.

[14]alexanderpotapov,b.tsegmed,l.v.ryzhakova,solarcyclevariationof“killer”electronsatgeosynchronousorbitandelectronfluxcorrelationwiththesolarwindparametersandulfwavesintensity,actaastronaut.2014,93,55–63.

[15]futakam.,taguchis.,etal.,neuralnetworkpredictionofenergeticelectronsatgeosynchronousorbitduringthestormrecoveryphase:effectsofrecurringsubstorms,ann.,geophys.2002,20(7),947-951.

[16]linga.g.,ginetg.p.,aneuralnetworkbasedgeosynchronousenergeticelectronfluxforecastingmodel.spaceweather.2010,82,s509003,doi:10.1029/2010sw000576.

[17]薛炳森，叶宗海,地球同步轨道高能电子增强事件预报方法[j].空间科学学报，2004,24(4):283～288.

[18]郭策，薛炳森，林兆祥，地球同步轨道高能电子通量预报方法研究，空间科学学报，2013，33(4)，418-426.

[19]xia0f.l.，zhangs.eta1,rapidaccelerationofradiationbeltenergeticelectronsbyz-modewaves.geophys.res.lett.，2010,39:l03103.

[20]zhangs,，xiaofl.chorus-drivenouterradiationbeltelectrondynamicsatdifferentl·shells.chinesephyslett,2010,27,12:129401.

[24]a.c.kellermanandy.y.shprits,ontheinfluenceofsolarwindconditionsontheouter-electronradiationbelt.journalofgeophysicalresearch,vol.117,a0521,doi:10.1029/2011ja017253,2012

[25]yousrfimr,kasmaeibs,vahabiea,etal.inputselectionbasedoninformationtheoryforconstructingpredictormodelsofsolarandgeomagneticactivityindices[j].solarphysics,2009,258(2):297–318.

[26]huangne,shenz,longsr.theempiricalmodedecompositionandthehilbertspectrumfornonlinearandnon-stationarytimeseriesanalysis.proceedingsoftheroyalsocietysoclond,1998,454(1971):903-995.

[27]vapnikvn.thenatureofstatisticallearningtheory.newyork:springer-verlag,2000:35-39.

[28]brown,r.g.andp.y.c.hwang.1992.introductiontorandomsignalandappliedkalmanfiltering,secondedition,johnwiley&sons,inc.

[29]yangpei-cai,houxiu-ji.onnon-sataionarybehaviorsandpredictiontheoryofclimatesystems.actameteorologicasinica,2005,(5)；556-570

[30]kim,k.c.,d.-y.lee,h,-j.kim,l.r.lyons,e.s.lee,m.k.ozturk,andc.r.choi(2008),numericalcalculationsofrelativisticelectrondriftlosseffect,j.geophys.res.,113,a09212,doi:10.1029/2007ja013011.

[31]satio,s.,y.miyoshi,andk.seki(2010),asplitintheouterradiationbymagnetopauseshadowing:testparticlesimulations,j.geopres.,115,a08210,doi:10.1029/2009ja014738

技术实现要素：

本发明的目的，在于提供一种基于经验模态分解的高能电子通量小时预报模型建立方法，克服了非平稳性给高能电子通量带来的预报困难，其分解后的数据序列与时间系数的原始数据序列相比具有更强的规律性，可以显著提高其预测的精度。

为了达成上述目的，本发明的解决方案是：

一种基于经验模态分解的高能电子通量小时预报模型建立方法，包括如下步骤：

步骤1，通过经验正交函数对高能电子通量的时间系数进行简化；

步骤2，采用经验正交函数的第一阶基函数对高能电子通量时间系数的经验正交函数重构展开；

步骤3，利用已知的电子通量以及基函数得到电子通量的时间系数；

步骤4，选择输入参数；

步骤5，采用经验模态分解emd算法对时间系数进行分解；

步骤6，对时间系数的各个分量进行拟合。

上述步骤1中，经验正交函数的表达式是：

其中，d表示天数，nt表示电子通量观测的世界时，1≤t≤24,m表示eof展开的阶数，ak表示第k阶时间系数，ek表示第k阶基函数。

上述步骤2中，经验正交函数重构展开的表达式是：

yreco(d,nt)＝a1(d)e1(nt)

其中，yreco为重构后的电子通量，a1(d)为随天数变化的第一阶时间系数，e1(nt)为电子通量随时间变化第一阶基函数。

上述步骤3中，时间系数的表达式是：

其中，d表示天数，nt表示电子通量观测的世界时，1≤t≤24,yreco为重构后的电子通量，a1(d)为随天数变化的第一阶时间系数，e1(nt)为电子通量随时间变化第一阶基函数。

上述步骤4中，选择输入参数采用正交最小二乘算法。

上述步骤5的具体过程是：

步骤51，识别出时间系数序列a1(t)中所有极大值点并拟合其包络线eup(t)，识别出序列a1(t)中所有极小值点并拟合其包络线elow(t)；

步骤52，计算上下包络线的平均值m1(t)：

步骤53，将a1(t)减去m1(t)得到h1(t)，再把h1(t)看做新的序列a1(t)；

步骤54，重复步骤51-53，经过n次的计算，直到h1(t)＝a1(t)-m1(t)满足imf条件为止，记作a1(t)＝h1(t)，则a1(t)为序列的第一个imf分量，并且是a1(t)中周期最短的分量；

步骤55，从a1(t)中分离出imf分量a1(t)，得到剩余分量：

r1(t)＝a1(t)-a1(t)

步骤56，将剩余分量r1(t)作为新的原始数据，重复步骤51-55得到其余的imf分量和一个余量，a1(t)被分解为：

其中，n表示时间系数经emd分解后分量的个数。

上述步骤6中，采用卡尔曼滤波算法对时间系数的各个分量进行拟合。

上述步骤6的具体过程是：

步骤61，假设线性模型如下：

ai＝φ^tθit

其中，ai为模型的输出，即后一天时间系数第i个分量的拟合值，φ^t为模型的输入，θit为第i个分量的线性滤波系数；r9为余量的拟合值，是余量的滤波系数，将各分量与余量的拟合值叠加得到后一天时间系数a1t的预报值；

步骤62，设定状态方程为：

θit+1＝θit+vt

θit为线性模型滤波系数，vt为随机的过程噪声；

步骤63，设定测量方程为：

zit＝φ^tθit+et

φ^t为线性模型中的输入，et为测量误差，zit为模型的输出；

假设随机噪声vt与测量误差et是相互独立且为正态分布的白噪声；

步骤64，采用卡尔曼滤波器更新算法，首先计算卡尔曼增益

由观测变量更新估计

θit+1＝θit+kit(zit-φ^tθit)

更新误差协方差

pit＝(i-kitφ^t)pit+q

将卡尔曼滤波器的目标转换成求误差的协方差p最小，得到p最小化时kit的取值k；

步骤65，模型的初始值θi0取由训练数据前30天做最小二乘回归得到的系数，初始协方差取010×10。

采用上述方案后，本发明采用经验正交函数的第一阶基函数对高能电子通量时间系数的经验正交函数(eof)重构展开；通过已知的电子通量以及基函数得到电子通量的时间系数；选择输入参数；采用经验模态分解emd算法对时间系数进行分解；对时间系数的各个分量进行拟合。本发明克服了非平稳性给高能电子通量带来的预报困难，其分解后的数据序列与时间系数的原始数据序列相比具有更强的规律性，可以显著提高其预测的精度。

附图说明

图1是使用eof分解后对数据进行简化，2001-2002年(训练集)高能电子通量前三阶基函数与所对应的时间系数；

图2是将数据进行emd分解，eofcoefficient为原始2001-2006年高能电子通量的时间系数，imf1-imf9以及r9时经过emd分解后时间系数的各分量；

图3是经卡尔曼滤波更新后2001-2006年高能电子通量时间系数各分量的预报值以及将各分量叠加之后eof系数的拟合值；

图4(a)、(b)、(c)是afteremd模型与beforeemd模型的时间系数预报值与真实值之间的对比，afteremd的预报结果更符合实际值；

图5(a)、(b)、(c)是为afteremd模型与beforeemd模型在2003和2006年误差绝对值对比；afteremd模型的误差绝对值低于emd模型；

图6(a)、(b)、(c)是分别为2005年(1月-6月)、2005年(7月-12月)、2006年(1月-5月)afteremd模型解与beforeemd模型重构后的正午12点电子通量对比；afteremd模型预报的电子通量更接近观测值；

图7是afteremd模型比beforeemd模型更能及时的反应高能电子通量在磁暴间的变化；

图8是本发明的流程图。

具体实施方式

以下将结合附图，对本发明的技术方案及有益效果进行详细说明。

本发明提供一种基于经验模态分解的高能电子通量小时预报模型建立方法，包括如下步骤：

步骤1，通过经验正交函数对高能电子通量的时间系数进行简化；

步骤2，采用经验正交函数的第一阶基函数对高能电子通量时间系数的经验正交函数重构展开；

步骤3，通过已知的电子通量以及基函数得到电子通量的时间系数；

本发明以2001-2006年的高能电子通量、太阳风参数以及地磁指数进行分析。高能电子通量来源于美国noaa网站上goes10卫星上>2mev电子通量5min数据，(http://goes.ngdc.noaa.gov/),为了本发明的研究把数据处理成每1h时间分辨率的高能电子通量数据。太阳风参数来源于nasaomni数据库，地磁指数来源于世界地磁数据中心日本memambetsu台站。

由于一天有24世界时，每小时均值电子通量数据的维数给预测带来困难，本发明采用eof简化数据，eof的基本思想是在尽量不损失原数据信息的条件下简化数据。将2001年-2002年的电子通量数据作为训练集，2003-2006年数据作为测试集(高能电子通量来源于美国noaa网站上goes10卫星上>2mev的真实数据)。电子通量由eof展开得：

其中，d表示天数，nt表示电子通量观测的世界时(1≤t≤24),m表示eof展开的阶数，ak表示第k阶时间系数，ek表示第k阶基函数。

图1为2001-2002年电子通量eof的展开结果，通过计算发现第一阶基函数的贡献率高达99.46％，且通过计算发现不同年份间的基函数相关性高达99％，结果表明，即使每年的电子通量随时间变化都不一样，但基函数基本保持不变。因此可以用第一阶基函数反应所有电子通量的变化。eof重构展开式为：

yreco(d,nt)＝a1(d)e1(nt)(2)

其中，yreco为重构后的电子通量，a1(d)为随天数变化的第一阶时间系数，e1(nt)为电子通量随时间变化第一阶基函数。

同时，通过已知的电子通量以及基函数，可以得到时间系数的表达式：

从时间序列的角度来看，通过对eof的简化，原时间序列中较平稳的一部分被剔除，剩余的时间系数更能体现高能电子通量随时间变化的非平稳、非线性的特点。因此，对高能电子通量的预报等价为对时间系数的预报。

步骤4，选择输入参数；

高能电子通量主要受太阳活动的影响，有准27天增强事件，其与磁暴有关。但是高能电子通量事件与磁暴并非简单的映射关系，统计发现只有一半磁暴会引起高能电子通量增强，这个特性增加了预报的难度(参考文献23)。由于目前的高能电子通量预报模型输入参数主要为磁暴指数ap和太阳风参数等，所以对于这种非线性和非稳定性的事件进行预报，目前的输入参数需要优化和筛选。

传统筛选参数的方法为观察参数与电子通量的相关系数，但是若输入参数过多，相关系数不能刻画众多参数与电子通量的相关性。为了解决线性相关系数的不足，本发明采用正交最小二乘算法(参考文献6)。其基本思想是在输入数据中按照回归的大小筛选出最重要的几个参数，使得回归结果的误差最小。

为了拟合后一天的时间系数，以太阳风参数(密度(ns)、速度(vs)、动压(p))(参考文献4、24)[e.g,paulikasandbalke,1979；a.c.kellermanandy.y.shprits,2012]地磁指数(ap、dst、ae)(参考文献25、8)[e.g,yousefimretal,2009；kaorisakaguichetal,2012]以及时间系数a1(参考文献6)[lishengetal；2017]的1-3天的历史数据作为输入参数的考虑范围总共21个输入参数，对训练数据进行正交最小二乘算法回归，得到下表前10个贡献度最高的参数为本模型输入参数。

表1由正交最小二乘算法筛析后模型的10个输入参数

步骤5，采用经验模态分解emd算法对时间系数进行分解；

emd是一种基于信号自身特点的信号分解算法。该算法不但完美地克服了小波变换中难于选取小波基与确定分解尺度的不足，还汲取了小波变换多分辨的优点，因此其更适合用于分析非线性非平稳信号。与此同时，其又是一种自适应的信号分解方法，可以根据时间系数变化的特点进行分解，保留了原数据当中较为重要的信息，所以高能电子通量的时间系数可以用该方法进行研究和分析。emd方法的基本思想：所有的复杂信号都是由简单的imf组成，并且这里每个imf都是彼此相互独立的(参考文献26、27)。也就是说，该方法可以把时间系数数据序列中真实存在的不同尺度的分量逐一的分解出来，产生一系列具有不相同特征尺度的数据序列，其分解后的数据序列与时间系数的原始数据序列相比具有更强的规律性，这可以显著提高其预测的精度。

时间系数的emd分解步骤如下：

1)首先，识别出时间系数序列a1(t)中所有极大值点(极小值点)并拟合其包络线eup(t)(elow(t))；

2)然后，计算上下包络线的平均值m1(t)，

3)最后，将a1(t)减去m1(t)得到h1(t)，再把h1(t)看做新的序列a1(t)，重复上述步骤，经过n次的计算，直到h1(t)＝a1(t)-m1(t)满足imf条件为止，记作a1(t)＝h1(t)，则a1(t)为序列的第一个imf分量，并且是原始序列中周期最短的分量。

从时间系数序列中分离出imf分量a1(t)，得到剩余分量：

r1(t)＝a1(t)-a1(t)(5)

将剩余分量r1(t)作为新的原始数据，重复步骤1-5可得到其余的imf分量和一个余量，原始序列a1(t)就可以被分解为：

其中，n表示时间系数经emd分解后分量的个数，例如若有10个分量，那么n＝10。

这里采用rilling等对huang等人提出的限定标准差(standarddeviation，sd)准则的改进的终止条件。

将时间系数通过emd方法分解如下图2，共得9个imf分量和1个余量，与原先变化剧烈的时间系数相比，分解后的各分量平稳程度越来越高，实现了时间系数序列的平稳化。emd方法是本发明的创新之处，将此方法运用于电子通量预报模型上对预测效率的提高会有较大的帮助。根据简化后的电子通量数据以及经emd分解后的时间系数，提出有关高能电子通量的预报模型：

其中，e1(nt)为第一阶基函数，a1-a9以及r9为后一天时间系数各分量的拟合值，y(d,nt)是后一天微分通量的1小时预报值。

步骤6，对时间系数的各个分量进行拟合。

利用线性模型以及筛选出的10个参数对时间系数的各分量进行拟合，得到十个简单的静态线性模型。但是在以往的研究中发现，例如太阳风与高能电子通量不是简单的映射关系，与其他参数之间也存在着非线性的关系，所以用静态线性模型拟合不能够刻画高能电子通量与各参数之间物理联系。由于卡尔曼滤波能动态的调整模型的参数，刻画动力学方程，因此使用卡尔曼滤波改进线性模型。

假设线性模型如下：

ai＝φ^tθit(8)

其中，ai(i＝1.....9)为模型的输出，即后一天时间系数第i个分量的拟合值，φ^t为模型的输入，这里是一个包含10个参数的向量，θit为第i个分量的线性滤波系数。r9为余量的拟合值，是余量的滤波系数，将各分量与余量的拟合值叠加得到后一天时间系数a1t的预报值。

卡尔曼滤波是一种递归处理数据的算法，通过估计过程的状态，使估计均方误差最小，以此得到状态的最优估计(参考文献28)。它的特别之处在于引入系统状态的概念并及时修正滤波系数，把卡尔曼滤波技术引用到时间系数分量的拟合，其对应的系统方程以及算法如下：

(1)状态方程：

θit+1＝θit+vt(11)

θit为线性模型滤波系数，vt为随机的过程噪声。

(2)测量方程：

zit＝φ^tθit+et(12)

φ^t为线性模型中的输入，即10个参数组成的向量，et为测量误差，zit为模型的输出。

假设随机噪声vt与测量误差et是相互独立且为正态分布的白噪声：

vt～n(0，q)(13)

et～n(0，r)(14)

实际操作中，过程噪声的协方差q和测量误差的协方差r会因计算迭代而变化，这里假设为常数。

(3)卡尔曼滤波器更新算法

计算卡尔曼增益

由观测变量更新估计

θit+1＝θit+kit(zit-φ^tθit)(16)

更新误差协方差

pit＝(i-kitφ^t)pit+q(17)

卡尔曼滤波器的目标是得到状态的最优估计，即误差最小化，因此可以转换成求误差的协方差p最小，式(13)就是p最小化时k的取值。对下一时刻的最佳估计是当前时刻的状态变量与测量变量的校正的结合，由式(14)体现。式(15)是有关误差协方差的更新。

(4)初始值与噪声的选择

初始值与噪声大小的选取是卡尔曼滤波预报的关键。针对本发明的模型，经多次尝试，模型的初始值θi0取由训练数据前30天做最小二乘回归得到的系数，初始协方差取010×10，测量误差的协方差由训练集的回归系数可得，过程噪声q一般是没有办法得到的，经过多次尝试，取10^-5i10×10时预报效果最好。

emd将时间系数分解成10个分量，经卡尔曼滤波预报更新后结果如图3所示，从图中可以看出10个分量具有不同程度的非平稳性，且低频分量的平稳性比高频分量的平稳性有较大的提高，因为通过emd分解后，原始时间系数中不同特征尺度的分量被逐一分解出来，分解后的数据相比原数据序列平稳性得到了较大的提高，虽然这些分量还是存在着非平稳性，但是由于它们具有不同的特征尺度，它们之间的相互影响被消除(参考文献29)[yang；houetal.,2005]。尤其是在imf5之后，预测效果相当理想，误差接近0，预测精度接近1。将10个分量的拟合值进行叠加，得到后一天的时间系数的拟合值。相比于仅用卡尔曼滤波预报更新得到的预报值，基于emd方法得到的预报值，精度有较大的提高。为了验证emd方法的有效性，将仅用卡尔曼滤波拟合的预报值和基于emd方法拟合的预报值与原时间系数相对比。为了方便讨论，将仅用卡尔曼滤波拟合电子通量方法的模型称为beforeemd模型,基于emd方法拟合电子通量的方法为afteremd模型。图4是两种模型2003年与2006年时间系数的预报值与真实值之间的对比。从图4中可以看出经过卡尔曼滤波长时间的调整后，到2006年拟合效果已经较为理想。并且基于emd方法拟合的预报值相比于仅用卡尔曼滤波预报的拟合值更接近真实的时间系数，尤其是在高能电子通量变化剧烈的情况下，未经emd分解的时间系数由于数据序列剧烈的非平稳性，预测值与真实值之间误差较大，最严重的情况比如2003年5月的预测值比真实值相差几个量级，这给预报精度带来很大影响，通过emd方法，尽可能的解决了非平稳给预报带来的困难，在电子通量突然的情况时极大的降低了预测带来的误差，与实测值更加贴近。图5展示了它们所得预报值与真实值之间误差的绝对值。

通过计算发现，2003年beforeemd模型的预报值的误差绝对值的总和达到731，afteremd模型预报值的误差绝对值总和为470。2006年未使用emd方法的误差绝对值总和为302，使用emd方法后误差总和降低到239。可以看出非平稳性对时间系数的预测有较大的影响。

为了进一步说明emd方法的有效性，分别计算了、2003-2006年的均方根误差和线性相关系数。从表2可发现，经过emd处理后的时间系数均方根误差有较大的减小，线性相关系数也有较大的提升，进一步验证了emd方法的有效性。如表2所示，可以发现：

1)对于2003年的预测值，afteremd预测方法的均方根误差和相关性都要优于beforeemd预测方法，均方根误差由2.84减小到1.78，相关性由0.81增大到0.92；

(2)对于2004年的预测值，afteremd预测方法的均方根误差和相关性也都要优于beforeemd预测方法，均方根误差由2.61减小到2.16，相关性由0.84增大到0.92；

(3)对于2005年的预测值，afteremd预测方法的均方根误差和相关性同样都要优于beforeemd预测方法，均方根误差由3.52减小到2.45，相关性由0.79增大到0.89。

由2003年到2005年时间系数均方根误差与线性相关系数的对比，可以发现emd方法在解决非平稳时间序列上的有效性。

均方根误差的表达式为：

线性相关系数的表达式：

预报率的表达式：

其中,fi为预报值，fi为观测值，为预报值的均值，为测量值的均值，n为样本个数。

表2afteremd模型的均方根误差、线性相关性都高于beforeemd模型

高能电子通量预报结果及分析

将预报所得时间系数结合基函数得出后一天每1h>2mev高能电子通量的预报值。具体表达式如下:

其中，e1(nt)为第一阶基函数，a1-a9以及r9为后一天时间系数各分量的拟合值，y(d,nt)是后一天微分通量的1小时预报值。为了方便观测与讨论，且正午时间段的观测值更有实际应用价值，本发明选取了通过beforeemd模型与afteremd模型在2005年与2006年正午12时的预报值与观测值间的对比。图7中黑色实线表示观测到的高能电子通量的对数值，红色实线表示afteremd模型预报的高能电子通量对数值，蓝色虚点线表示beforeemd模型的预报值，无论是emd模型，还是afteremd模型预报的电子通量值与观测值在趋势上基本一致，这说明了经过正交最小二乘法筛选后的参数能够较好的拟合高能电子通量的预报值，验证了正交最小二乘法在筛选参数上的可行性。不过，高能电子的快速损失使得电子通量变化不规律[30][kimetal.,2008；saitoetal.,2010]，这也是造成高能电子通量非平稳性的主要原因，因此模型的预报由于实际电子通量的剧烈变化，很难跟的上电子通量观测值的变化，从图7中可以明显地观察到beforeemd模型的预测值的变化要比实测值滞后，这是因为高能电子通量数据序列具有极强的非平稳性。通过图7中的after模型可以看出，经过emd方法的处理，极大的减小了非平稳性给数据预测带来的影响。在高能电子突然变化的情况下，由于减小了非平稳性的影响，模型的预测值能尽可能的贴近了实际观测值的变化。在处理高能电子通量“拐点”的问题上有较好的效果。此外以往模型所做的日积分电子通量模型的时间分辨率过低，不能预测由于高压冲击造成电子损失的快速变化(参考文献31)[k.sakaguchetal.,2012]。因此本发明所建立的每小时均值预报模型，在时间分辨率上高于日积分电子通量的预报，可以通过每小时电子通量的变化，即时地观察到电子损失的变化。

磁暴期间高能电子通量结果分析

高能电子尤其在磁暴期间变化十分剧烈，对磁暴间的高能电子通量的预报尤为关键，可以及时的对卫星采取保护措施，从而降低卫星损害的风险。因此本发明选取了2005年7月9日-7月19日磁暴期间的电子通量的变化。预报结果如图7，其中实线为观测值，虚线为afteremd模型预报值，虚点线为beforeemd模型预报值，从图7的比较中可以发现，beforeemd模型的预报值比实际观测值有至少1天的滞后，这是由于高能带电子的剧烈变化，非平稳性因素增强，预报结果受非平稳性的干扰，不能及时的反应实际高能电子的变化。通过emd方法的处理，极大的降低了非平稳性的影响，尽管是在高能电子剧烈变化的磁暴期间，afteremd模型能够及时的反应实际电子通量的变化，这从高能电子通量的预报意义上来说尤为重要。

下表3分别展示了afteremd模型与beforeemd模型在2003-2006年的预报率、线性相关系数以及均方根误差。从表中可以看出基于emd方法的模型在各方面预报结果上优于beforeemd，从而进一步验证了emd方法在预报高能电子通量上的有效性。

表3>2mev每小时均值高能电子通量(+1day)的预报效率、线性相关系数、均方根误差及均值；afteremd模型预报结果优于beforeemd模型

上述为本发明提出的每小时均值高能电子通量的预报模型，由于现有的模型主要以高能电子通量的日积分值进行预报，为了验证本发明所提出的模型的可行性，将日均值电子通量作为观测值，带入本发明提出的emd模型(日积分电子通量来源于美国nasaomni数据库)，表4为本发明提出的emd模型从2001年-2006年的日积分电子通量预报值的预报结果，经过长时间的调试，预报结果趋于稳定，到2006年，预报率达到最高，为0.86。其他年份预报结果也较为理想。为了验证预报结果的好坏，选取了以下四种模型进行比较：径向扩散模型[lietal；2004]、low-e模型[turnerandli；2008]、combo模型[turnerandli；2008]、nict模型(参考文献8)[k.sakaguchietal；2012],这四种模型较为经典或已在网上运行。结果如表5所示，从表五中的比较可以发现，在2005-2006年emd模型的预报率略高于其他模型，但2001-2004年提前一天的预报率有较大提高。这主要是由于每年高能电子通量的非平稳程度不一样，2001年-2004年电子通量的变化较为剧烈，磁暴次数较多，因此电子通量的变化的非平稳性很强，给电子通量的预报带来较为严重的干扰，因此以往的模型忽视了非平稳性对数据序列的影响，预测精度一直没有较大的突破，从本发明的预报结果可以进一步说明非平稳性对高能电子通量预报的影响，也验证了本发明所提出的emd模型的可行性。表4emd模型>2mev日积分电子通量(+1day)预报效率与均方根误差

表5geo高能电子通量(+1day)预报率比较；emd模型预报结果最为理想

本发明提出的emd模型通过emd方法极大的降低了高能电子通量中非平稳性给预报带来的影响，通过与仅用卡尔曼滤波所得预报结果的比较，发现每小时均值电子通联不管是预报效率、相关性都得到较大提高，其中2006年的最高预报效率可达0.84。为了与日积分电子通量预报模型做比较，同时利用emd模型预报了日积分电子通量值，结果表明在电子通量变化较为剧烈的年份，预报结果与以往经典模型相比有较大的提高。由于在以往的预报模型中几乎都忽视了非平稳性的影响，希望通过本发明的验证可以使更多的研究者关注高能电子通量变化的非平稳性给预报带来的影响。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。