基于BVAR模型的施工期库岸边坡变形预测方法

基于bvar模型的施工期库岸边坡变形预测方法
技术领域
1.本发明属于变形监测技术领域，涉及基于bvar模型的施工期库岸边坡变形预测方法。


背景技术：

2.水利工程库岸边坡安全稳定关系到水利工程的整体安全，一旦出现滑坡，后果将十分严重。针对水利工程库岸边坡的各种安全监测项目中，变形监测最为直观和重要，是判断边坡安全稳定的直接依据。因此，基于库岸边坡变形监测数据，建立边坡变形监控模型，对判断边坡变形趋势和安全性态而言极为重要。施工期边坡变形监测数据较少，并且没有环境量(如水位、降雨、气温)监测数据，难以建立传统的以各种环境量为影响因子的统计模型或基于机器学习的非线性模型。相对于传统建模方法而言，时间序列模型通过数据自身所蕴含的结构与规律来预测监测数据的未来值，具有高效、稳定的特点。目前已有许多学者将时间序列模型应用于边坡工程领域，准确预报了边坡的变形趋势。传统的时间序列模型在施工期边坡变形监测数据较少的情况下，其建模效果往往较差。在众多时间序列模型中，贝叶斯向量自回归(bvar)模型在小样本数据建模方面具有明显的优势。通过众多学者的研究成果可知，bvar模型的精准度相对更高，采用该方法对施工期边坡变形监测数据进行建模分析能够解决传统方法中建模困难、建模效果差等关键问题。综上所述，施工期边坡变形监测数据的建模和预测分析是目前亟待解决的一个难题。


技术实现要素：

3.本发明的目的是提供基于bvar模型的施工期库岸边坡变形预测方法，解决了现有施工期边坡变形监测数据量少、无环境量监测数据的问题。
4.本发明所采用的技术方案是，基于bvar模型的施工期库岸边坡变形预测方法，具体按照以下步骤实施：
5.步骤1，获取边坡变形监测原始数据作为初始样本，对原始数据进行adf检验；
6.步骤2，构建bvar模型；
7.步骤3，通过吉布斯采样器获取bvar模型待估计参数和模型预测值；
8.步骤4，获得边坡变形的最终预测值及预测区间。
9.本发明的特点还在于：
10.步骤1采用matlab软件中自带的adftest函数对边坡原始变形监测数据进行adf检验，通过adf检验判断监测数据的平稳性。
11.步骤2具体按照以下步骤实施：
12.步骤2.1，检验bvar模型最优滞后阶数；
13.步骤2.1具体为：采用matlab库bayvar中的order_var_d.m程序计算监测数据的bic值，将bic最小值作为bvar模型的最优滞后阶数；为避免模型自由度过多，将order_var_d.m程序中最大滞后阶数取为5阶；
14.步骤2.2，设定bvar模型超参数组合；
15.首先，给定超参数的可能性取值，超参数中，衰减参数d＝{0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1}，总体紧度γ＝{0.01,0.05,0.1,0.2,0.3,0.4,0.5,1}，相对紧度w＝{0.1,0.2,0.3,0.4,0.5,1}；
16.步骤2.3，根据theil u统计量确定最优超参数组合；
17.采用bayvar中的err_d.m脚本程序对上述所有超参数组合进行试算，得到模型精度评价指标，即theil u统计量；theil u统计量越小证明预测误差越小；将试算得到的所有超参数组合的theil u统计量进行比较，将theil u统计量最小的超参数组合作为最优超参数组合。
18.步骤3具体为：采用bayvar中的bvar_gs_d.m脚本程序获取bvar模型待估计参数和模型预测值。
19.步骤3将吉布斯预迭代次数设定为5000，以确保参数的收敛性，然后舍弃原来的迭代，将吉布斯有效迭代次数设定为1000；bvar模型系数矩阵β中各元素及协方差矩阵σ中各元素均为待估计参数，吉布斯采样器估计参数。
20.步骤3吉布斯采样器估计参数的具体流程为：
21.步骤3.1，指定先验分布的均值作为系数矩阵β和协方差矩阵σ的初始值参数的初始值；
22.步骤3.2，在初始值和监测数据已知的条件下，结合系数的先验信息得到系数矩阵的后验分布，并在满足稳定性条件的基础上对系数进行抽样，以得到系数的估计值；
23.步骤3.3，在更新后的系数矩阵和监测数据已知的条件下，结合协方差矩阵中参数的先验信息计算其后验分布，并在满足稳定性条件的基础上对其进行抽样，以得到协方差矩阵的估计值；
24.步骤3.4，把迭代过程中所产生的bvar模型参数、模型参数的迭代序列存储起来；
25.步骤3.5，重复步骤3.2-步骤3.3，经过足够次数的迭代后，模型参数的迭代序列收敛且平稳；
26.步骤3.6，通过吉布斯迭代序列的概率分布得到模型预测值样本。
27.步骤4中由吉布斯采样得到1000个边坡变形预测值样本，通过求取边坡变形预测值样本的中值得到bvar模型的最终预测值，并将样本中的最大值、最小值作为bvar模型的预测区间。
28.步骤1中若监测数据是非平稳的，则adftest函数输出变量h＝0；若监测数据是平稳的，则adftest函数输出变量h＝1；若原始监测数据非平稳，则将原始监测数据进行n阶差分后(n＝1,2,3
…
)再次进行adf检验，直至转变为平稳的数据。
29.本发明的有益效果是：基于预测区间来考虑边坡变形监测中的不确定性因素，对边坡变形监测数据进行建模和预测分析，具有较好的效果。有效解决了传统建模方法中施工期边坡变形监测数据量少、无环境量监测数据而造成的建模困难这一关键问题，为边坡的建模分析提供新的可行方法。bvar模型解决了传统方法建模困难、建模效果差的难题，在边坡的建模和预测分析方面具有较好的预测效果，推进水利工程中的边坡变形预测方面的分析研究。
附图说明
30.图1是本发明基于bvar模型的施工期库岸边坡变形预测方法的流程图；
31.图2是本发明中多点位移计m07zpr原始监测数据过程线；
32.图3是本发明中测点m07zpr(5m)滞后阶数检验结果；
33.图4是本发明中测点m07zpr(10m)滞后阶数检验结果；
34.图5是本发明中测点m07zpr(20m)滞后阶数检验结果；
35.图6是本发明中测点m07zpr(30m)滞后阶数检验结果；
36.图7是本发明中采用bvar模型对测点m07zpr(5m)监测数据进行建模得到的十步预测值及预测区间；
37.图8是本发明中采用bvar模型对测点m07zpr(10m)监测数据进行建模得到的十步预测值及预测区间；
38.图9是本发明中采用bvar模型对测点m07zpr(20m)监测数据进行建模得到的十步预测值及预测区间；
39.图10是本发明中采用bvar模型对测点m07zpr(30m)监测数据进行建模得到的十步预测值及预测区间。
具体实施方式
40.下面结合附图和具体实施方式对本发明进行详细说明。
41.本发明基于bvar模型的施工期库岸边坡变形预测方法的流程，如图1所示，具体按照以下步骤实施：
42.步骤1，获取边坡变形监测原始数据作为初始样本，对原始数据进行adf检验；
43.采用matlab软件中自带的adftest函数对边坡原始变形监测数据进行adf检验，通过adf检验判断监测数据的平稳性。若监测数据是非平稳的，则adftest函数输出变量h＝0；若监测数据是平稳的，则adftest函数输出变量h＝1。若原始监测数据非平稳，则将原始监测数据进行n阶差分后(n＝1,2,3,...)再次进行adf检验，直至转变为平稳的数据。
44.步骤2，构建bvar模型；
45.步骤2.1，检验bvar模型最优滞后阶数；
46.步骤2.1具体为：采用matlab库bayvar中的order_var_d.m程序计算监测数据的bic值，将bic最小值作为bvar模型的最优滞后阶数；为避免模型自由度过多，将order_var_d.m程序中最大滞后阶数取为5阶。
47.步骤2.2，设定bvar模型超参数组合；
48.首先，给定超参数的可能性取值，超参数中，衰减参数d＝{0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1}，总体紧度γ＝{0.01,0.05,0.1,0.2,0.3,0.4,0.5,1}，相对紧度w＝{0.1,0.2,0.3,0.4,0.5,1}。
49.步骤2.3，根据theil u统计量确定最优超参数组合；
50.采用bayvar中的err_d.m脚本程序对上述所有超参数组合进行试算，得到模型精度评价指标，即theil u统计量。theil u统计量越小证明预测误差越小。将试算得到的所有超参数组合的theil u统计量进行比较，将theil u统计量最小的超参数组合作为最优超参数组合。
51.步骤3，通过吉布斯采样器获取bvar模型待估计参数和模型预测值。
52.步骤3具体为：采用bayvar中的bvar_gs_d.m脚本程序获取bvar模型待估计参数和模型预测值。将吉布斯预迭代次数设定为5000，以确保参数的收敛性，然后舍弃原来的迭代，将吉布斯有效迭代次数设定为1000。
53.bvar模型系数矩阵β中各元素及协方差矩阵σ中各元素均为待估计参数，吉布斯采样器估计上述参数的流程包括：
54.步骤3.1，指定先验分布的均值作为系数矩阵β和协方差矩阵σ的初始值参数的初始值；
55.步骤3.2，在初始值和监测数据已知的条件下，结合系数的先验信息得到系数矩阵的后验分布，并在满足稳定性条件的基础上对系数进行抽样，以得到系数的估计值；
56.步骤3.3，在更新后的系数矩阵和监测数据已知的条件下，结合协方差矩阵中参数的先验信息计算其后验分布，并在满足稳定性条件的基础上对其进行抽样，以得到协方差矩阵的估计值；
57.步骤3.4，把迭代过程中所产生的bvar模型参数、模型参数的迭代序列存储起来；
58.步骤3.5，重复步骤3.2-步骤3.3，经过足够次数的迭代后，模型参数的迭代序列收敛且平稳；
59.步骤3.6，通过吉布斯迭代序列的概率分布得到模型预测值样本。
60.步骤4，获得边坡变形监测数据的最终预测值及预测区间。
61.由吉布斯采样得到1000个边坡变形预测值样本，通过求取边坡变形预测值样本的中值得到bvar模型的最终预测值，并将样本中的最大值、最小值作为bvar模型的预测区间。
62.本发明中bvar模型原理如下：假设一个p阶滞后的var模型共包含m个变量：{y1,y2,
…
,ym}，并且在该模型中不存在常数项，则var模型系统表示为如下形式
[0063][0064]
式中，a
i,j,τ
为第i个方程中变量yj的τ阶滞后项y
j-τ
的系数，u
it
为第i个方程中的随机扰动向量。上述var模型中需要确定2m2p个超参数：包括m2p个先验均值δ
i,j,τ
和m2p个先验方差
[0065]
贝叶斯方法可以把先验信息引入到var模型中，大大减少了模型系统中需要估计的参数数量。本发明采用明尼苏达先验，包含以下多个基本假定：
[0066]
(1)正态性：随机扰动向量u
t
＝(u1,u2,
…
,um)
t
服从多元正态分布nm(0,σ)；
[0067]
(2)独立性：协方差矩阵σ和模型系数a
i,j,τ
相互独立；
[0068]
(3)协方差矩阵σ的先验分布取为扩散先验分布，即π(σ)
∝
|σ|-(m+1)/2
；
[0069]
(4)模型系数a
i,j,τ
相互独立，并且a
i,j,τ
服从正态分布
[0070]
(5)a
i,j,τ
的先验期望是δ
i,j,τ
，它的取值由公式(2)确定；
[0071][0072]
(6)标准差s
i,j,τ
由公式(3)确定，其中γ为总体紧度；
[0073][0074]
其中，g(τ)为τ阶滞后变量相对一阶变量的紧度，本发明中的g(τ)函数采用公式(4)所示的调和紧度函数；
[0075]
g(τ)＝τ-d
,d＞0
ꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0076]
f(i,j)为相对紧度函数，它的取值按照公式(5)确定，其中w
ij
是介于0到1之间的常数。因此，参数f(i,j)的选取转化为确定超参数w的大小；
[0077][0078]
si/sj为第i个方程与第j个方程的自回归残差标准差之比，si为第i个方程对常数项和其p阶滞后项的最小二乘回归的残差标准差得到。因此，通过公式(3)，所有参数的先验方差的大小转化为确定三个超参数(γ、d和w)的大小。
[0079]
超参数(γ、d和w)大小的选取标准由预测精度确定。在bvar模型系统中，预测精度是由theil u统计量决定的。theil u统计量是基于所研究的预测模型的均方根误差与基于随机游走预测模型的均方根误差之比，即：
[0080][0081]
式中，a
t
表示监测数据在t处的实际值，f
m,t
表示基于预测方法m的预测值，而f
rw,t
是基于随机游走模型的预测值。如果u＝1.0，则说明方法m的预测精度与随机游走模型的预测精度基本相同；如果u《1.0，说明方法m的预测精度比随机游走模型的预测精度更高；如果u》1.0，则说明方法m的预测精度比随机游走模型的预测精度低。u值越接近于0说明预测效果越好。故u值达到最小时，即可确定最优超参数组合。
[0082]
如果把式(1)中的模型转换为如下形式的不相关模型：
[0083]
yi＝zβi+εi,i＝1,2,
…
,m
ꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0084]
其中：
[0085]
[0086][0087]
那么，模型参数的明尼苏达先验分布简记为：
[0088]
(β|σ)～n(μ0,v0),π(σ)
∝
|σ|-(m+1)/2
ꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0089]
此处μ0的分量为0或1，而协方差矩阵v0是由n块对角矩阵构成的，其中的每一个对角矩阵的元素由组成，非对角块元素均为零。根据贝叶斯公式，参数(β,σ)的联合后验分布密度函数为：
[0090][0091]
其中：
[0092]
上述后验分布包含了对参数的信息，如均值、方差，但密度函数太复杂，所以对每个参数的后验估计通过吉布斯采样来实现，每次只更新一个分量，把高维问题转化为低维。
[0093]
吉布斯采样次数足够大时，迭代序列才能达到平稳分布，即去除初始不稳定状态的迭代序列样本点。因此，需设定预迭代次数和有效迭代次数，以确保迭代序列的收敛性。通过吉布斯迭代序列的概率分布可以得到bvar模型预测值样本，然后将预测值样本的中值作为bvar模型的最终预测值，将样本中的最大值、最小值作为bvar模型的预测区间。
[0094]
实施例
[0095]
本发明基于bvar模型的施工期库岸边坡变形预测方法的实施例，选取黄金峡水利枢纽左岸边坡多点位移计m07zpr监测数据建立bvar模型，其原始监测数据过程线，如图2所示。图2中多点位移计m07zpr包含m07zpr(5m)、m07zpr(10m)、m07zpr(20m)和m07zpr(30m)四个测点。
[0096]
通常情况下，将已有的样本数据t分成两个时期t1和t～t1。t1时期内的数据用于估计bvar模型参数和预测，t～t1时期内的数据用于计算比较预测误差并确定最终的超参数数值。为了有足够的数据来进行模型的估计与预测，每个测点选用前122个监测数据作为模型的起始估计区间段，向前预测十期。
[0097]
对原始数据进行(augmented dickey-fuller，adf)增广迪基-福勒检验，是为防止边坡变形监测数据在建模过程中出现“伪回归”现象，保证结果的有效性和稳健性。
[0098]
采用matlab软件中自带的adftest函数对m07zpr(5m)、m07zpr(10m)、m07zpr(20m)和m07zpr(30m)四个测点中的原始监测数据进行adf检验。由于原始数据无法通过adf检验，因此对所有数据进行一阶差分，其adf检验结果如表1所示，
[0099]
表1为本实施例中边坡变形监测数据的adf检验结果
[0100][0101]
注：adf检验中的c为常数项，t为趋势项，k为滞后阶数；滞后阶数k的选择标准是以贝叶斯信息准则(bic)最小的原则。
[0102]
一阶差分处理后的数据均在5％的显著水平下平稳，可用于建立bvar模型。
[0103]
采用bayvar中的order_var_d.m脚本程序计算监测数据的bic值，检验m07zpr(5m)、m07zpr(10m)、m07zpr(20m)和m07zpr(30m)四个测点的最优滞后阶数。为避免模型自由度过多，此处将最大滞后阶数设定为5。根据bic最小值确定上述四个测点的最优滞后阶数，因此m07zpr(5m)测点最优滞后阶数为2阶，如图3所示，m07zpr(10m)测点最优滞后阶数为3阶，如图4所示，m07zpr(20m)和m07zpr(30m)两个测点最优滞后阶数均为1阶，如图5～图6所示。
[0104]
设定bvar模型超参数(衰减参数d、总体紧度γ和相对紧度w)可能性取值的集合。其中，d＝{0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1}，γ＝{0.01,0.05,0.1,0.2,0.3,0.4,0.5,1}，w＝{0.1,0.2,0.3,0.4,0.5,1}。m07zpr(5m)、m07zpr(10m)、m07zpr(20m)和m07zpr(30m)四个测点的超参数组合将从上述集合中选取。
[0105]
采用bayvar中的err_d.m脚本程序对上述四个测点的所有超参数组合进行试算，检验模型精度，即theil u统计量。theil u统计量越小证明预测误差越小。将试算得到的所有超参数组合的theil u统计量进行比较，得到theil u统计量最小的超参数组合，如表2所示，并作为bvar模型的最优超参数组合。
[0106]
表2是本发明中bvar模型theilu统计量最小的超参数组合
[0107][0108]
采用bayvar中的bvar_gs_d.m脚本程序获取bvar模型待估计参数和模型预测值。本发明中迭代次数采用程序中的默认设置，吉布斯预迭代次数设定为5000，以确保参数的
收敛性，然后舍弃原来的迭代，将吉布斯有效迭代次数设定为1000。在吉布斯迭代序列收敛且平稳后，得到边坡变形监测数据的十步预测值样本(共10
×
1000个预测样本数据)，计算样本数据的中值得到边坡变形监测数据的十步预测值，计算样本数据的的最大、最小值得到边坡变形监测数据的预测区间。m07zpr(5m)、m07zpr(10m)、m07zpr(20m)和m07zpr(30m)四个测点的监测数据十步预测结果如图7～图10所示。由图7～图10可知，在上述四个测点的预测结果中，有90％的预测值位于bvar模型预测区间内，表明bvar模型在区间预测方面具有较好的短期预测效果。
[0109]
实施例总结
[0110]
为了综合评价所建立bvar模型的预测效果，将传统var模型和bvar模型十步预测的残差平方和(sse)、平均绝对误差(mae)和均方根误差(rmse)进行比较，如表3所示。
[0111]
表3是本发明中var模型和bvar模型十步预测误差的比较结果
[0112][0113]
由表3可知，m07zpr(5m)、m07zpr(10m)、m07zpr(20m)和m07zpr(30m)四个测点的bvar模型十步预测误差普遍小于传统var模型。bvar模型的三种误差指标(sse、mae和rmse)均在0.1以下，说明bvar模型的整体预测效果较好。通过上述实施例预测结果表明，采用bvar模型对施工期边坡变形监测数据进行短期预测具有较好的效果，能为实际工程提供依据。