【基于时间序列的变形观测的数据处理方法的研究】

小谷加油

已于 2022-11-18 23:00:41 修改

阅读量546

点赞数 1

文章标签： python 时序数据库

于 2022-11-18 23:00:09 首次发布

本文链接：https://blog.csdn.net/weixin_48101151/article/details/127931100

版权

基于时间序列的变形观测的数据处理方法的研究

引言

通过采用基于时间序列的变形观测的数据处理方法的研究来研究沉降观测记录，是一个典型的时间序列类问题和线性回归预测模型，对此分三种方法进行研究，寻找最为合理的数据处理方法。

方法一：多项式拟合法

1.什么是多项式

$\begin{aligned} y &=a_{0}+a_{1} x+a_{2} x^{2}+\cdots+a_{n} x^{n} \\ &=\sum_{i=0}^{n} a_{i} x^{i} \end{aligned}$

这是多项式拟合的输出表达式

2 怎样拟合？

考虑输入样本数据如下：
$\begin{aligned} x &=\left[x_{1}, x_{2}, \ldots, x_{m}\right]^{T} \\ y &=\left[y_{1}, y_{2}, \ldots, y_{m}\right]^{T} \end{aligned}$
构建方程组：
$\begin{aligned} y_{1} &=a_{0}+a_{1} x_{1}+a_{2} x_{1}^{2}+\cdots+a_{n} x_{1}^{n} \\ y_{2} &=a_{0}+a_{1} x_{2}+a_{2} x_{2}^{2}+\cdots+a_{n} x_{2}^{n} \\ & \vdots \\ y_{m} &=a_{0}+a_{1} x_{m}+a_{2} x_{m}^{2}+\cdots+a_{n} x_{m}^{n} \end{aligned}$
令
$\begin{gathered} A=\left(\begin{array}{ccccc} 1 & x_{1} & x_{1}^{2} & \cdots & x_{1}^{n} \\ 1 & x_{2} & x_{2}^{2} & \cdots & x_{2}^{n} \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ 1 & x_{m} & x_{m}^{2} & \cdots & x_{m}^{n} \end{array}\right) \\ b=\left(\begin{array}{c} a_{0} \\ a_{1} \\ \vdots \\ a_{n} \end{array}\right) \\ f=\left(\begin{array}{c} y 1 \\ y_{2} \\ \vdots \\ y_{m} \end{array}\right) \end{gathered}$
则上面的方程组可以写为：
$A b = f$

3.拟合过程就是求b的过程

现在一个很直观的想法，就是：
$\begin{gathered} A b=f \\ b=A^{-1} f \end{gathered}$
但是，没这么简单。考虑一下矩阵维度
$A_{[m, n+1]} b_{[n+1,1]}=f_{[m, 1]}$

4.施密特正交化

$\begin{aligned} A b &=f \\ Q R b &=f \\ R b &=Q^{T} f \\ b &=R^{-1} Q^{t} f \end{aligned}$

拟合的原理就讲到这里

5.程序实现

采用matlab中自带的polyfit函数和ployval函数进行拟合程序如下：

clc,clear;
x0=[0 0 0 0 0 0 0 0 0 0 0 0 0 0];%2009年10月25日
x1=[10 12 9 13 13 12 14 13 14 12 14 13 10 12];%2010年3月30日数据
x2=[16 17 17 19 19 18 19 19 20 21 21 21 17 18];%2010年5月30日数据
x3=[23 24 23 24 26 27 24 24 26 25 26 26 25 25];%2010年7月25日
x4=[29 30 28 31 34 35 32 31 34 32 34 33 33 34];%2010年10月25日
x5=[32 35 30 40 42 40 39 35 39 40 45 41 39 40];%2011年1月5日
x6=[33 36 30 42 45 43 40 39 41 42 46 44 42 42];%2011年3月30日
x=[x0;x1;x2;x3;x4;x5;x6];
t=[0 156 217 273 365 437 521]; %根据时间段求得的
avar=mean(x,2)';%计算出平均累计沉降
y=-avar;

%%画图
xx=0:1:600;
[p,S ]=polyfit(t,y,2); 
yy=polyval(p,xx);
plot(xx,yy);
hold on;
for i=1:1:14
    plot(t,-x(:,i),'r*');
end
title('沉降观察记录函数预测')
xlabel('距离基础浇筑前的时间间隔 （天）');
ylabel('累计沉降深度 mm');

%%数据输入
year=input('请输入预测日期年份：');
mouth=input('请输入预测日期月份:');
data=input('请输入预测日期日:');
if year<2009&&mouth<10&&data<25
    disp('输入时间太早错误');
end
if mouth>12 or data>31
     disp('输入时间错误');
else 
   time=data-25+(mouth-10)*30+(year-2009)*365;%计算时间间隔
   h=polyval(p,time);
   disp(['预测的下降深度为:',num2str(h)]);
end

运行结果如下：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-KKPAGBiA-1668783365530)(C:\Users\小谷加油\AppData\Roaming\Typora\typora-user-images\1647531887228.png)]

可以输入需要预测的年份，月份，日期就可以计算距离第一次监测的日之差，通过时间点来预测接下来的沉降深度。

通过计算以及函数拟合，发现二次多项式拟合最好。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ti4x8bZh-1668783365530)(C:\Users\小谷加油\AppData\Roaming\Typora\typora-user-images\1647532637320.png)]

相关数据参数和范数

R: [3×3 double]
df: 4
normr: 3.5185

方法二:简单的移动平均法

1.相关推导过程

设时间序列为简单移动平均公式为
$Y_{1}, Y_{2}, \cdots, Y_{t}, \cdots$
简单移动平均公式为
$M_{t}=\left(Y_{t}+Y_{t-1}+\cdots+Y_{t-N+1}\right) / N \quad t \geqslant N$
式中，通常将 Mt 表示为 t 期移动平均数; 用 N 表示为移动平均数的项数。该公式表明: 当 t 向前移动一个时
期，就会增加一个新观测的数值，并且会去掉一个远期的数据，最后可以得到一个新的平均值。正是这些“吐故纳新”，逐期向前移动，所以把其称为移动平均法。
$M_{t-1}=\left(Y_{t-1}+Y_{t-2}+\cdots+Y_{t-N}\right) / N$
两式比较可知，这是移动平均法的递推公式。当Y 比较大的时候，利用移动平均法的递推公式可以在预测
计算中减少计算量，在编程中用这个递推公式可以简单很多。由于移动平均法可以平滑观测的数据，消除一些
不确定因素所产生的影响，使长期趋势显示出来，因而可以用于预测。预测公式为
$Y_{t+1}=M_{t}$

2.求解思路

（1）不考虑观测时间之间的差异，将时间看成连续累加的时间序列 1,2，3,4,5,6,7…,

（2）将相关累积沉积量取平均值

（3）动态求取2 3 4 5次项的MSE数值,取最小值下的次项作为移动次项

（4）根据公式计算预测值

3.相关程序

clc,clear; %%简单的移动平均法
x0=[0 0 0 0 0 0 0 0 0 0 0 0 0 0];%2009年10月25日
x1=[10 12 9 13 13 12 14 13 14 12 14 13 10 12];%2010年3月30日数据
x2=[16 17 17 19 19 18 19 19 20 21 21 21 17 18];%2010年5月30日数据
x3=[23 24 23 24 26 27 24 24 26 25 26 26 25 25];%2010年7月25日
x4=[29 30 28 31 34 35 32 31 34 32 34 33 33 34];%2010年10月25日
x5=[32 35 30 40 42 40 39 35 39 40 45 41 39 40];%2011年1月5日
x6=[33 36 30 42 45 43 40 39 41 42 46 44 42 42];%2011年3月30日
x=[x0;x1;x2;x3;x4;x5;x6];
t=1:1:7; %根据时间段求得的
y=mean(x,2)';

for k=1:5
m=length(y);
n=[2 3 4 5]; %n为移动平均项数
for i=1:length(n) 
    for j=1:m-n(i)+1
        yhat{i}(j)=sum(y(j:j+n(i)-1))/n(i); %求y的预测值
    end
    MSE(i)=1/(m-n(i))*sum((y(n(i)+1:m)-yhat{i}(1:m-n(i))).^2); %求MSE 
end
[ans,p]=min(MSE);
p=p+1;
data=6;
data=data+k;
y=cat(2,y,[yhat{p}(m-n(p)+1)]);
fprintf('第%d次记录',data);
fprintf('下降的深度：%f\n',yhat{p}(m-n(p)+1));
fprintf('最优的MSE为：%f\n',ans);
fprintf('移动项数为：%d\n',p);
end
times=1:1:length(y);
for i=1:1:14
    plot(t,-x(:,i),'r*');
    hold on;
end
plot(times,-y);
title('沉降观察记录函数预测')
xlabel('测量次数');
ylabel('累计沉降深度 mm');

4.相关运行结果

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-O6uSLAaK-1668783365531)(file:///C:\Users\小谷加油\AppData\Roaming\Tencent\Users\3043037377\QQ\WinTemp\RichOle\GO9QTJJ8FIQU{XP9FTPNUGA.png)]

第7次记录下降的深度：36.952381
最优的MSE为：95.536480
移动项数为：2
第8次记录下降的深度：38.555556
最优的MSE为：80.577546
移动项数为：2
第9次记录下降的深度：38.621693
最优的MSE为：69.067874
移动项数为：2
第10次记录下降的深度：38.043210
最优的MSE为：60.528508
移动项数为：2
第11次记录下降的深度：38.406820
最优的MSE为：53.836171
移动项数为：2

方法三:指数平滑法

1.基本原理

基于一次 / 二次指数平滑法的预测模型
1) 基于一次指数平滑法的模型。对风功率观测数据建立一次指数平滑模型，其预测公式为：
$p_{t+1}^{(1)}=\alpha * p_{t}+(1-\alpha) * p_{t}^{(1)}$

式中：为 t+1 时间点的一次功率预测值，即当前时间点的一次平滑值； Pt 为 t 时间点的功率实际值，为 t 时间点的一次功率预测值，即上个时间点的一次平滑值； α 为平滑系数。一般选择。

一次指数平滑法适用于风功率短期预测。
2）基于二次指数平滑法的模型。二次指数平滑法是在一次指数平滑的基础上再平滑，同时也结合了历史平均和变化趋势，对风功率观测数据建立二次指数平滑模型，其预测公式为：
$\begin{gathered} p_{t+T}=a_{t}+b_{t} * T \\ p_{t+1}^{(2)}=\alpha * p_{t+1}^{(1)}+(1-\alpha) * p_{t}^{(2)} \\ \left\{\begin{array}{l} a_{t}=2 * p_{t+1}^{(1)}-p_{t+1}^{(2)} \\ b_{t}=\frac{\alpha}{1-\alpha}\left(p_{t+1}^{(1)}-p_{t+1}^{(2)}\right) \end{array}\right. \end{gathered}$
式中： P^t+T为 t+T 时间点的功率预测值， T 为由 t 时间点向后推移时间点数； P（2） t+1 为 t+1 时间点的二次功率预测值，即当前时间点的二次平滑值； P（2） t 为 t 时间点
的二次预测值，即上个时间点的二次平滑值。一般选择
$P_{0}^{(2)}=P_{0}^{(1)}=P_{10}$

2.代码实现

from turtle import onclick
import matplotlib.pyplot as plt
# 支持中文
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

# 一次指数平滑预测
def es1(list3, t, a):
    if t == 0:
        return list3[0]  # 初始的平滑值取实际值
    return a * list3[t - 1] + (1 - a) * es1(list3, t - 1, a)  # 递归调用 t-1 → 12
 
# 二次指数平滑预测
def es2(list3, t, a):
    if t == 0:
        return list3[0]
    return (a * es2(list3, t - 1, a) + (1 - a) * list3[t - 1])
 
def answer3(list2):
    # 指数平滑法
    a = 0.8  # 平滑常数
    listES = []  # 指数平滑值的列表
    em=[]
    for i in range(len(list2)):
        if i == 0:
            listES.append(list2[i])
            continue
        s = a * list2[i] + (1 - a) * listES[-1]
        listES.append(s)
    print("指数平滑值的列表：{}".format(listES))
 
    # 一次指数平滑预测
    t = len(list2)  # 预测的时期 13
    x = es1(list2, t, a)
    print("下一个数的一次指数平滑预测：{}".format(x))
 
    # 二次指数平滑预测
    m=1  # 预测的值为之后的第m个
    yt = es2(list2, t - 1, a)
    ytm = listES[t - 2]
    esm = ((2 * ytm - yt) + m * (ytm - yt) * a / (1 - a))
    print("之后的第{}个数的二次指数平滑预测：{}".format(m, esm))
    
    # 画图
    plt.scatter(list(range(len(listES))), listES)
    plt.title('指数平滑法的趋势调整')
    plt.xlabel('测量 次数')
    plt.ylabel('累计沉降量 mm')
    plt.show()
    
if __name__ == '__main__':
    list1 = [ 0,	-12.2	,-18.7	,-24.8,	-32.14,	-38.3	,-40.3]  # 
    n = 3  # 移动平均期数
    answer3(list1)  # 指数平滑法

3.相关运行结果

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-tCvPz7Q0-1668783365532)(D:\desktop\Figure_1.png)]

数平滑值的列表：[0, -9.76, -16.912, -23.222400000000004, -30.356480000000005, -36.711296, -39.5822592]
下一个数的一次指数平滑预测：-39.5822592
之后的第1个数的二次指数平滑预测：-125.81225600000002

数平滑值的列表：[0, -9.76, -16.912, -23.222400000000004, -30.356480000000005, -36.711296, -39.5822592]
下一个数的一次指数平滑预测：-39.5822592
之后的第1个数的二次指数平滑预测：-125.81225600000002

总结

本文通过对三种方法的对比，简单的研究了基于时间序列的变形规则的数据处理方法，有多项式拟合法，平移移动法和指数平滑法。

多项式拟合法：具有直观和精确的特点，对近期数据内的点有良好的拟合特性，也可以直观算出某一天的具体累计沉降量。但是缺点是对长时间跨度无法精准的求解，且误差较大。

平移移动法：忽略了测量时间，以次数作为时间序列，可以直观精确预测后几次的数据，其残差较小。但缺点是无法与测量日期相连，缺少数据可信度。

一次和二次指数平滑法，构建了建筑物沉降预测模型及方法，利用 python 语言编写了相应的预测程序。以对不同平滑系数下的风功率预测结果进行了对比分析，并计算了预测误差，验证了方法的可行性。

综合判断，其各有优缺点，因此需要综合考虑找到合适的沉降预测模型。

参考文献

[1]池其才,周世健,王奉伟.基于时间序列的变形监测数据处理方法比较研究[J].测绘与空间地理信息,2015,38(07):193-195.

[2]徐彦农,王一帆,王浩淳,田辰蔚,张兆欣,李昕潞.基于一次/二次指数平滑法的风功率预测方法[J].南方农机,2021,52(21):26-28.

[3]江东,付晶莹,黄耀欢,庄大方.地表环境参数时间序列重构的方法与应用分析[J].地球信息科学学报,2011,13(04):439-446.