应用随机过程课设论文

Harbin Institute of Technology

课程设计（论文）

课程名称：设计题目：院班级：设计者：学号：指导教师：设计时间：

哈尔滨工业大学

哈尔滨工业大学课程设计任务书

随机线性模型建模

--地震震级预测模型

本次建模，本人采用的数据是1970年1月1日至1982年12月31日期间的实测的四

川地区的地震震级（见表1，平均每15天进行1次测量，共有323个数据）。本文将利用表中的数据建立该地地震等级的随机线性模型，并对下几个时段的震级进行预测。

表1 1970年1月1日至1982年12月31日期间四川地区的地震震级实测值

为了清楚起见，可将上述数据画于图1中。易见该时间序列基本符合平稳随机序列的特征。

876

地震强

度

54320

100

150200测量次数

250

300

350

图1 上述时间范围内实测的地震等级示意图

1. 线性模型的建立：

确定平稳时间序列线性模型的步骤可归纳为五个步骤，以下结合本文的数据具体介绍建模的过程：

(1) 数据的观测：

对一个时间序列作n次测量得到一个样本Z1,Z2,,Zn，一般取n50。在本文中对于1970年1月1日至1982年12月31日期间的的四川地区的地震震级进行测量，得到了n=323个样本值Z1,Z2Z323(列于表1中)。

(2) 数据预先处理：

利用已知的样本值Z1,Z2Z323来计算样本的均值Z，并作变换tZtZ得n个数据

1,2,,n，n=323。

其中 Z

n

Zj=

j1

Z323

j1

323

=4.3152

同时得到一个n=323的新序列1,2,,n。

(3) 样本自协方差函数、自相关函数、偏相关函数的计算：

，，时，在计算样本自协方差函数样本自相关函数偏相关函数k0,1,2,,k，kkkk

一般取k

，常用k

n10

。在本文中取k



n10



32310

30。

 样本自协方差：

11k22knkn1k

nk



j1

jk,k0,1,2,k(kn)

其中 n=323, k=30。从而得到ˆk的相应数值，如下表所示：

表2 ˆk的数值

 样本自相关函数：

rkr0,k0,1,2,,k(kn)中，得到与ˆ相将上面求得的自协方差函数带入kk

ˆk的数值如下表所示，并绘出相应的样本自相关函数的曲线以方便分析。对应的

ˆk的数值

表3 

样本自相关函数

51015ｋ

202530

图2 样本自相关函数

 样本偏相关函数：

为求出样本的偏相关函数，需要求解如下的Yule-Walker方程：

112k1

11

21



k1

21

k1

2

1

kk

1

1



2

k

对上述方程两边去估计值仍然成立。当k分别取1,2,3„k时，分别求出ˆ11,ˆ22ˆkk的

值如下表所示(其中为方便起见，设ˆ00=1)。

表3 样本偏相关函数的值

为了便于分析与比较，同样画出ˆkk的曲线，见图3.

样本偏相关函数

图3 样本偏相关函数图

(4) 模型识别：

由图2、3或表2、3可以看出样本自相关函数与偏相关函数的取值具有如下特点：

ˆk拖尾：  样本自相关函数

ˆk的点图(图2)可见，当k逐渐增大时，ˆk越变越小，故可判断根据样本自相关函数ˆk具有拖尾性。 

 样本偏相关函数ˆkk在kp1处截尾：

中至多有一个使易见而当k5时，平均20个kkkk



0.1113，则认为

kk截尾在kp5处。

综上，可以判断本模型为AR(5)模型。

(5) 参数估计：

对于AR(5)模型计算参数估计值ˆ1,ˆ2,ˆ3,ˆ4,ˆ5，采用矩估计的方法，具体公式如下：

1ˆ1

1

ˆ22

ˆp

p1



1



p1

21

1

1

1

2

1

2p

1



2p

1

其中p=5。

将已知数据值带入上式可得：

ˆ1ˆ2ˆ3ˆ4ˆ5



1.0000

0.3411

0.33230.3181

0.3495

0.34111.00000.34110.33230.3181

0.33230.34111.00000.34110.3323

0.31810.33230.34111.0000 0.3411

0.3495



0.3181



0.3323



0.34111.0000

1

0.34110.1495

0.33230.13860.31810.1134 0.34950.16410.34570.1557

从而可得AR(5)模型为：

(Zt)=1(Zt1)2(Zt2)3(Zt3)4(Zt4)5(Zt5)at

tkl，

取

并在上式两边分别取估计值可得：

ˆZ)=ˆ(Zˆˆˆ(ZZ)ˆ2(ZZ)ˆ5(ZZ) kl1kl1kl2kl5

将已知数据带入上式并整理得此AR(5)模型的预报公式为：

ˆ0.1495ZˆˆˆˆˆZ+ 0.1386Z+0.1134Z+0.1641Z+0.1557Z+1.2026klkl1kl2kl3kl4kl5

k

li,(n

式中当

ˆZZkli。 1i,时，有2kli

(6) 模型预测

由前面得到的AR(5)模型的预报公式可对以后几个时段的该地地震等级的值进行预

测。

利用Z319=3.7,Z320=4.6,Z321=3.1,Z322=4.4,Z323=4.3可得下几个时段（每15天为一个时段）的地震等级的预测值：

ˆZ0.14954.3+ 0.13864.4+0.11343.1+0.16414.6+0.15573.7323 +1.2026=4.1378

ˆZ0.14954.1378+ 0.13864.3+0.11344.4+0.16413.1+0.15574.6324 +1.2026=4.1411

ˆZ0.14954.1411+ 0.13864.1378+0.11344.3+0.16414.4+0.15573.1325 +1.2026=4.0875

ˆZ0.14954.0875+ 0.13864.1411+0.11344.1378+0.16414.3+0.15574.4326 +1.2026=4.2476



2. 非线性模型分析

在现实生活中更多的模型并不能很好的符合线性模型，这就需要引入非线性模型的概念。非线性模型的本质处理手段是将非线性过程线性化处理。结合我自身的课题,在非线性信号参量的估计当中采用卡尔曼滤波体系。卡尔曼滤波来源于winer过程，所以在利用

卡尔曼滤波时要求噪声参量是高斯白噪声。最早的Kalman滤波是卡尔曼（R.E.Kalman）于1960年提出的从与被提取信号有关的观测量中通过算法估计出所需要信号的一种滤波算法。他把状态空间的概念引入到随机估计理论中，利用系统噪声和观测噪声的统计特性，以系统的观测量作为滤波器的输入，以所要估计的值（系统状态或参数）作为滤波器的输出，滤波器的输入与输出之间是由时间更新算法（状态方程）和测量更新算法（测量方程）联系在一起的。时间更新由上一步的测量数据更新结果和设计卡尔曼滤波器时的先验信息确定，测量更新则在时间更新的基础上根据实时获得的测量数据确定。由于所用的信息都是时域内的量，所以不但可以对平稳的一维随机过程进行估计，也可以对非平稳的多维随机过程进行估计。卡尔曼滤波是最小方差估计。

Kalman滤波算法的计算流程如下：已知线性空间：

x(n)F(n,n1)x(n1)Q(n1)





y(n)H(n)x(n)R(n)

观测值n)

状态方程测量方程

图2-1 系统方程描述

出于对非线性滤波的需要，出现了EKF和UKF算法，均采用Kalman滤波为框架，区别在于对非线性函数的处理。

广义卡尔曼滤波（Extended Kalman Filter，EKF）是将非线性模型在最优状态估计处进行泰勒展开，相对于其它算法的综合性能较好，算法简单，易于实现，是比较常用的非线性滤波方法。如果可以精确得到状态方程和观测方程的高阶导数项，则可以得到真实的统计量。但在实际系统中，这一点很难满足，一般仅可获得一阶和二阶导数。因为EKF算法采用泰勒展开的线性化处理方式，所以只有当系统的状态方程和观测方程都接近线性且连续时，滤波结果才有可能接近真实值，对于强非线性系统，EKF滤波性能极不稳定，甚至发散。

自从1997年Julier和Uhlman提出无迹卡尔曼滤波（Unscented Kalman Filter，UKF）以来，已逐渐应用于导航、跟踪领域。其采用非线性变换（UT变换）代替传统的线性变换，体现了一种先进的思想。这种先进思想就是“非线性估计算法应更接近系统的非线性本质”。UKF算法以非线性变换UT变换为核心，通过某种采样策略在原先状态分布中选取一组Sigma点，使这些点的均值和方差等于原状态分布的均值和方差，对状态向量的后验概率密度函数（PDF）进行近似化。然后将这些点带入到非线性函数中，得到相应的非线性函数值点集，然后在测量的基础上调节样本点的位置，使得样本均值和样本方差分别以二次精度逼近实际分布的均值和方差（EKF只能达到一阶精度）。由此获得更多的观测假设，使得对系统状态的均值和协方差的估计更为准确。