时间序列与灰色系统组合模型

§7.2时间序列与灰色系统组合模型7.2.1灰色系统概述灰色系统理论是华中科技大学教授邓聚龙教授于20世纪70年代末至80年代初提出，已广泛应用于社会、经济、农业、生态、生物等各个领域。灰色系统是指信息部分明确、部分不明确的系统，已知的信息称为白色，未知的信息称为黑色。它通过对原始数据的重新生成，特别没有规律的原始数据序列通过累加或累减处理而成为具有较强规律性的新数列，再用微分方程来描述这一新的数列，解此微分方程即得到自变量与因变量的关系。7.2.2GM（1,1）模型建立GM（1,1）模型的实质是对原始序列做一次累加生成的序列呈现一定的规律，然后建立一阶线性微分方程模型，求得拟合曲线对系统进行预测。（1）GM（1,1）预测模型设有原始序列（7-40）将其累加生成新数列，i=1，2，…，n（7-41）(0)(),1,2,,Xiin(1)()Xi其中，(1)(0)1()()ijXiXj相应的微分方程为（7-42）(1)(1)dXaXudt累加矩阵为(1)(1)(1)(1)(1)(1)1[(1)(2)]121[(2)(3)]121[(1)()]12XXXXBXnXn（7-43）常数向量为(0)(0)(0)(2),(3),,()TnYXXXn（7-44）应用最小二乘法求得解的系数得（7-45）1()TTnaBBBYu并带入微分方程的解，得到时间函数（7-46）(1)(1)()((0))atuuXtXeaa再求导还原得（7-48）(0)(0)ˆ(1)((1))auuXtaXea这两个方程就是GM(1,1)模型的预测方程。此时其实际的预测值可由下式得（7-49）(0)(1)(1)ˆˆˆ(1)(1)()XtXtXt（2）GM（1,1）预测模型的检验方法根据GM（1,1）模型的预测方程可采用3种检验方法，残差的大小检验、关联度检验和后验方差检验。设t时刻的残差为（7-50）(0)(0)ˆ()()()etXtXt残差的均值为（7-51）11()nteetn残差的方差为（7-52）22111(())ntSeten原始序列的方差值为（7-53）2(0)(0)2211(()())ntSXtXtn后验差比值和小误差概率为后验方差检验的两个重要数据。显然，C越小，表示越大而越小；大说明原始数据的方差大，即原始数据的离散度大；小说明残差方差小，残差的离散程度小。因12SCS2()0.6745PPekeS2S1S2S1S此，C小表示尽管原始数据的离散程度高，但模型计算所得的值与实际值的差并不太离散。根据C和P的大小可以综合评定预测模型的精度，具体指标见下表精度等级PC好0.950.35较好0.800.50合格0.700.65不合格0.700.657.2.3组合预测模型在时间序列预测实践中，对于某一时间序列预测问题，可用各种预测模型进行预测。一般来说，采用预测模型不同，预测结果也不同；一种更为科学的做法是，讲不同的预测方法进行适当的组合，这就是组合预测方法。为了有效地利用各种模型所提供的信息提出了组合预测方法。组合预测的类型一般分为两种综合类型：一种是权重组合；另一种是区域综合。（1）权系数组合预测模型组合预测的关键是恰当地确定各个预测模型的权系数。权重模型可以用下式表示（7-54）1ˆˆ()()JjjjXtXt为了保证模型的无偏性，应满足如下约束条件j（7-55）11Jjj其中，为不同模型组合预测值，为不同模型的预测值，为不同模型的权系数。ˆ()Xtˆ()jXtj（2）区域综合组合预测模型区域综合组合模型解求的是多种预测值置信区间的交集。设J种预测值有置信区间（7-59）ˆ()()jjXtt则的置信区间是这J个区间的交集ˆNlx（7-60）1ˆˆ()()()JjjjjjXXtt若上式是空集，则依次排除最大与最小的预测值置信区间，若剩余的模型超过半数扔由上式进行区域预测，否则重新建模。7.2.4AR(p)模型与灰色系统线性组合预测模型下面我们对AR(p)模型与灰色系统的最优组合预测模型作出分析。利用（7-54）式，J个模型组合预测的形式为j=1,2,…，Jt=1,2,…，N（7-61）1ˆˆ()()JjjjXtXt而组合预测模型的预测误差可以表示为（7-62）11ˆˆ()()()(()())()JJjjjtjjetXtXtXtXtej在极小化准则minQ，下，可得21(())NtQetmin,..1TTTQQeeWEWstRW如果各个预测模型的预报误差是不相关的，则E矩阵是可逆矩阵，按最小二乘法求的最优权向量为（7-64）1110()TWRERER可以证明最优组合预测模型的精度高于任何一个单一预测模型的精度。对于只有两个预测模型的情况，矩阵E可表示为：（7-65）11122122EEEEE§7.3频域分析方法信号频域分析方法是傅里叶变换及离散傅里叶变换快速算法，它把一个信号分解为各个不同的频率分量，使信号的时域特征与频域特征联系起来，成为信号分析处理的有力工具。如对一个实测变形时间序列进行分析时，为了求得变形中的主频率与振幅，可以先用傅里叶变换初步确定时间序列中的主频率，根据主频率利用最小二乘法模拟时间序列求出系数，在对模拟值与实际观测的残差序列重复上述过程，直到残差序列中不在出现突出的主频率。7.3.1经典谱分析与现代谱分析经典的傅里叶分析有两种方法：一种是直接法，又称周期图法，它是直接对数据X（N）进行FFT，然后取其幅度平方得到信号的功率谱（7-66）21()()SXN这种方法计算方便，是目前常用的方法。另一种是间接法，又称自相关法，它是先对数据求自相关函数，然后对自相关函数作FFT得到功率谱，即（7-67）()xrm()()MimmMSrme1mN这两种方法本质是一样的，都把有限时段的数据看做是无限抽样序列的开窗截取后的结果。但这将出现两个问题：一是频率分辨率的极限取决于抽样数据的长度；二是发生信号频谱“泄漏”现象，即功率谱主瓣内的能量泄漏到旁瓣内。7.3.2短时傅里叶分析和小波分析经典的傅里叶分析和现代谱分析方法都是见建立在平稳信号的处理基础上，给出的结果显示了信号总体所包括的各种频率成分。短时傅里叶变换是对信号加上有限的时窗后作傅里叶分析。（7-69）(,)()()jtFxtWtedt对于非平稳信号的分析，短时傅里叶变换存在时频分辨率固定不变和产生干扰项的缺陷，因而应用受到局限。总之，频域分析处理技术的发展，使信号分析处理能力上升到一个更高的层次。§7.4变形动态响应分析7.4.1动态响应分析在分析监测桥梁在风荷载作用下的震动，高层建筑在风力、温度作用下的摆动等这类变形，一般采用连续的自动记录装置，得到一组时间序列数据。在进行时间序列数据处理时，往往要解释时序数据出现的原因，如在变形分析中探讨变形量与变形原因之间的关系，也就是寻求动态变形原因与变形之间的关系。因此，为了把分析变形描述成具有系统理论计算的动态模型系统，除了把原有变形体的变形描述成时间过程外，还把变形原因描述成时间过程；因此一个动态系统的性能可以表达成（7-70）()[(),(),]XtfutXtt其中，为变形体的变形，为激励。()Xt()ut7.4.2最小二乘动态响应分析动态响应分析除了时域分析外还可以用频谱分析。最小二乘响应分析是对输入信号中所包含的谐波与分量用频谱分析法分析，然后用这些频率采用三角函数多项式来模拟输入和输出信号。对所选频率在模拟输入、输出信号中是否重要，可以用显著性检验方法。如用选择的k个频率模拟的输入信号为，为了检验某频率jw是否重要，可以先除去jw，由k-1个信号对u(t)进行模拟，则可以求得模拟信号，若满足下式()utˆ()ut()ˆ()jut（7-73）()ˆˆˆ()()()()jutututut则说明频率jw对，模拟u(t)是很重要的。§7.5变形时序分析的应用实例为了说明时域分析方法，我们选取某一做坝高75m、坝顶长220m的面板堆石坝的安全监测观测数据，下表为数据序号位移残差序号位移残差序号位移残差序号位移残差10.0-0.611316.30.992529.2-0.803752.520.7-1.131417.20.662630.9-0.323854.332.2-0.861517.70.052733.51.043955.943.8-0.481619.90.912835.61.924056.755.4-0.111720.20.002938.03.094158.067.20.461821.20.233038.22.074259.379.91.931921.41.253138.74360.0811.01.812021.62.283241.24461.5912.21.792123.51.603347.84562.61012.71.062225.01.333447.54663.01112.8-0.052324.62.953552.24764.11213.9-0.182425.33.483652.84865.37.5.1利用AR（p）模型与GM（1,1）模型及其组合模型预测分析（1）AR(p)模型显然1号点的垂直位移是非平稳的时间序列，它具有趋势特征，因此要首先消除趋势项。我们采用拟合函数法消除趋势项，拟合函数采用如下多项式模型（7-74）2012()kkXtaatatat拟合多项式时，对k采用添项法建模，从k=1开始逐步添加，比较每次的残差平方和，作统计检验有()rsk（1，N-k-1）（7-75）(1)()~()/(1)rskrskFrskNk当新添不能使残差平方和显著下降时，则说明拟合的多项式较好地表达了与t间的关系，但一般不要采用高次多项式。经过计算采用如下多项式1kt（7-76）()0.61241.2247Xtt根据建立AR(p)模型的要求，先对残差进行平稳性检验与正态性检验以满足要求，然后对残差进行零化处理和标准化处理。在利用处理后的数据建立AR(p)模型，经过模型辨识，按FPE原则，AR(p)模型的阶数p=2,模型系数为{-0.6726,0.2988}，即（7-77）12()0.6726()0.2988()nnnXtXtXt其中，这里为零化处理和标准化处理后的数据。()nXt（2）GM（1,1）模型利用第1期到第30期的垂直位移数据，建立GM(1,1)模型，并对第31期到第48期的位移进行预测，预测利用等维新信息模型，增加相同个数的新息与去掉相同个数的旧息同时进行。（3）组合模型为了建立组合预测模型，我们将AR(p)模型与灰色系统的预测模型进行组合预测，为了求得组合预测的系数，利用式（7-64）对AR(p)模型与灰色系统的预测模型预测的第31期到第42期数据计算组合预测系数，并利用结果组合预测第43期到第48期。7.5.2利用傅里叶变换进行频域分析为了说明傅里叶变换对时间时间序列处理方法，我们按下式（7-78）模拟的时间序列，模拟偶然误差为2，采样间隔为0.001，根据申农采集定律，可以分辨最高频率不超过500的信号。得到原始时间序列。（7-78）sin(250)sin(2120)ytte显然，直接从原始时间序列图中识别出规律是困难的。利用傅里叶变换，将时域信号转换到频域进行分析，得到下图的频谱。从变换后的频谱图可以看到两个明显的波峰频率分别为：50和120左右。在确定了主要频率之后，可以依下式，按最小二乘法来估计振幅和相位。01122sin(250)sin(2120)yAAtAt011sin(250)cos(250)AAtAAt22sin(2120)cos(2120)At