本发明涉及高拱坝变形时空序列预测技术领域，具体涉及一种结合BP网络、 GRU神经网络以及时空克里格方法对高拱坝变形时空序列进行预测。
背景技术：
随着现代传感器数据种类和数量的快速增长，如何合理有效的利用这些收集的海量信息预测大坝整体的变形值，进而把握高拱坝整体变形状态，已经成为水利工作者亟待解决的问题。只有了解大坝整体的状态，才能做到对发现的异常情况及时处理，防患于未然。才能保证大坝始终运行在安全高效的环境中。
大坝变形值通过埋设在坝体内部的位移传感器测量得到，位移传感器限于大坝实际结构，只能布设在有限的几处重点监测位置。由于测点布置相对离散，无法把握大坝空间整体状态。目前常见的做法是将大坝看作由一群三维空间点组成的有机整体，利用空间插值算法插补出未知点位移值。这一方法虽然达到了了解大坝整体状态的目的，但是它的弊端是只能了解当前时间点或历史时间上的变形状态，无法对整体变形状态做出预测。要预测大坝整体变形值，就是在空间插值的基础上加入时间维度，预测构成大坝的任意一点在任一时间的变形值。由于大坝变形的过程本质上是时空过程，所以设计出许多数据收集方法来记录数据中每个测点的空间和时间信息，这些信息被称为时空序列(Spatio-Temporal Series, ST)数据。本发明研究的问题可以理解为通过有限测点的大坝历史变形数据来预测未来任意时间任意位置的变形值，即时空序列预测。
时空数据由于在空间和时间两个维度上同时存在依赖关系，这就使得它与经典数据挖掘文献中研究的其他数据不同。许多广泛使用的数据挖掘方法都建立在数据实例独立且分布相同(Independent Identically Distributed,IID)的假设基础上。然而，该假设在处理时空数据时可能并不成立，这是因为数据实例在空间和时间结构上彼此相关，并且在不同的空间区域或时间段中会显示不同的属性。在数据分析过程中忽略这些依赖关系可能导致预测结果准确性不高，预测结论的可解释性较差。在时空数据研究领域，除限制传统预测算法的有效性之外，空间和时间信息的存在还使得在时空数据挖掘中应用新兴技术成为可能，如人工神经网络(Artificial Neural Network,ANN)和门限循环神经网络(Gated Recurrent Unit, GRU)等。因此应用神经网络技术，结合时间和空间维度，综合考虑大坝变形整体时空趋势和局部时空趋势，对未观测时空位置变形值进行建模和预测进而展示出大坝整体变形状态，是本发明所要解决的问题。
技术实现要素：
发明目的：基于时空一体化的高拱坝变形时空序列预测方法，解决现有技术中无法综合考虑时空序列预测中的时空相关性的问题。为综合考虑时间和空间因素，将时空演变过程分为整体和局部两部分，应用简单神经网络拟合整体时空趋势，采用传统时空预测模型拟合局部随机过程，构造一种新型时空序列预测算法。由于简单神经网络无法准确捕捉时间序列中的长期依赖性，为解决该问题引入门限循环网络，利用门限循环网络对时间序列数据高精度预测的优势优化简单神经网络。
技术方案：一种基于时空一体化的特高拱坝变形时空序列预测方法，包括以下三个方面：
(1)大坝传感器在监测数据时，由于各种现实因素，实测数据中时常带有随机误差项，假设时空过程分解为两部分：Zi(t)＝Mi(t)+ei(t)，式中Mi(t)代表确定性时空变化，ei(t)代表去除确定趋势后的小规模随机变化并且满足 E(ei(t))＝0，即期望为零。利用简单BP神经网络拟合整体时空趋势；
(2)去除整体时空趋势后得到残差项ei(t)，对残差项进行线性无偏估计，选用时空克里格方法拟合局部时空趋势；
(3)引入门限循环神经网络对大坝各测点时间序列进行预测，预测出相关测点变形值，将预测后的顺河向位移数据作为训练数据优化BP网络。
所述利用简单BP神经网络拟合整体时空趋势的具体步骤如下：
(1.1)确定BP网路拓扑结构：确定神经网络输入向量和目标输出向量，拟设输入层中神经元数目为4，输入向量分别对应各个移传感器的三维坐标和传感器测得位移时的时间。设输出层的节点数目为1，目标输出向量对应大坝测点顺河向位移值。网络结构中包含多个隐含层，关于隐含层个数、各隐含层中神经元个数设置的问题，目前没有合适的算法直接求出最佳个数，按神经元个数从小到大依次尝试后计算出最佳个数；
(1.2)确定网络结构权重系数：神经网络无需深入了解大坝具体受力情况，而是通过训练网络结构中的系数(权重)来达到时空预测的目的。BP网络在初始化时以[-1,1]之间的随机数对权重赋值；
(1.3)生成激活后的估计值：网络中的每个神经元可以在输出端使用非线性激活函数，产生到其他神经元的平滑信号。激活函数有很多种类，其中一种常用的非线性激活函数是由逻辑函数定义的S形传递函数；
(1.4)多次训练得到理想模型：将输入向量带入神经网络中，经过前向传播和反向传播算法的不断训练，将误差控制在理想范围内。即每次训练后计算期望输出与实际输出误差，误差通过反向传播算法不断调整神经元权重。所述内容 (2)去除整体时空趋势后得到残差项ei(t)，对残差项进行线性无偏估计，选用时空克里格方法拟合局部时空趋势。拟合局部时空趋势的具体步骤如下：
(2.1)残差数据应该尽可能充分，对于残差数据，两两计算时空距离；
(2.2)寻找一个拟合曲线拟合距离与时空半方差的关系，从而能根据任意时空距离计算出相应的半方差；半方差函数定义为rij＝0.5*E[(Zi-Zj)2]，根据地理学第一定律，即空间上相近的属性相近。rij表达了属性的相似度，空间上的相似度就用距离dij来表达。时空克里金插值假设rij和dij之间存在着函数关系，这种关系可以是线性、二次函数、指数、对数关系等等。为了确认这种关系，我们需要首先对观测数据集计算任意两个点的距离dij和半方差rij，这样会得到多个 (dij，rij)的数据对，得到函数关系r＝r(d)，将所有的d和r绘制成散点图，寻找一个最优的拟合曲线拟合d和r的关系，得到函数关系式r＝r(d)。
(2.3)计算出所有已知点之间的时空半方差rij；
(2.4)对于未知点Zo，计算它到所有已知点Zi的半方差rio；
(2.5)求解代价函数得到最优系数λi；
(2.6)使用最优系数对已知点的属性值(顺河向位移值)进行加权求和，得到未知点顺河向位移值的估计值，结束时空克里格方法。
所述内容(3)引入门限循环神经网络对大坝各测点时间序列进行预测，预测出相关测点变形值，将预测后的顺河向位移数据作为训练数据优化BP网络。 GRU网络进行大坝各测点时间序列预测的具体步骤如下：
(3.1)将大坝各测点时间序列预测问题转化为监督学习问题：使用Keras 深度学习框架实现GRU网络。Keras中的GRU模型将数据分为两类：输入向量和输出向量。应用到时间序列预测问题上，以上一个时间点(t-1)的观察值作为输入，当前时间点(t)的观察值作为输出来实现问题转换的目的。GRU网络在本次时间序列预测中应用多输入GRUs网络。多输入GRUs网络输入向量中不仅包括测点自身顺河向位移时间序列，还包括环境量数据。具体而言，输入变量中包括滞后时间点(t-1)的观察值，还包括(t-1)时刻环境变量值；
(3.2)变换时间序列数据，使其平稳。先判断时间序列是否平稳，如果不平稳需要去除时间趋势使其平稳。平稳化时间序列有消除时间序列趋势的方法是对数据进行差分运算，即(t)时刻位移值减去(t-1)时刻位移值，得到差分序列；
(3.3)将观测结果标准化；
(3.4)建立GRU神经网络：GRU循环网络包含输入层、隐含层和输出层。在添加GRU层的时候，首先需要确定输入层输入向量，输出层输出向量。初始时，随机生成权重矩阵，如步骤1.1、1.2；
(3.5)对于GRU来说，由前向传播计算出重置门r、更新门z、当前缓存记忆内容g、当前时间步的最终记忆h这四个向量的值。计算过程依次为：
rt＝σ(Wrxt+Urht-1)
zt＝σ(Wzxt+Uzht-1)
gt＝tanh(Wg(ht-1⊙rt)+xtUg)
ht＝yt＝ht-1⊙(1-zt)+zt⊙gt
以上式中，σ为Sigmoid激活函数(即S型函数)；xt表示时刻t的输入；h 代表GRU网络的输出，ht代表时间t时刻的最终输出，ht-1保存t-1时刻的记忆； W和U分别代表GRU网络中的权重系数，；
(3.6)通过BPTT计算误差项δj值：利用BPTT算法反向计算每个神经元的误差项δj本质上是误差函数E对神经元j输入量netj的偏导。BPTT算法将第l 层t时刻的误差项值沿两个方向传播：一是沿时间的反向传播，即从当前t时刻开始，计算每个时刻的误差项，与权重系数U相关；另一传播方向是是将误差项向上一层传播，与权重系数W相关；
(3.7)计算每个权重系数的梯度：GRU网络在训练之前对权重w赋予较小的随机值，根据BPTT算法不断训练权重。网络的权重通过纠正神经元的输出变量和理想输出之间的误差来改变并且及时调整。不断重复执行这一过程，使得每训练依次网络的权重就更新一次，直到最后所有训练数据的输出项达到理想的精度，结束训练。
本发明采用上述技术方案，具有以下有益效果：
1.选择BP神经网络拟合整体时空趋势而不是直接利用传统时空预测方法，目的是解决传统时空预测方法对非平稳时空序列预测精度不高的问题。
2.现有方法普遍将时间和空间数据性质分开来考虑，无法结合时空相关的特性对时空序列做出统一预测。本发明的大坝变形时空序列预测算法能结合时间和空间两方面因素考虑，提高了时空序列预测的精确度。
3.由于简单BP神经网络无法准确捕捉时间序列中的长期依赖性，为解决该问题引入门限循环网络GRU，利用GRU网络对时间序列数据高精度预测的优势优化简单神经网络。
附图说明
图1为部署在坝体内的传感器分布图；
图2为时空序列预测问题形式化示意图；
图3为本发明方法实施例的总体框架图；
图4 BP网络结构示意图；
图5为GRU误差反向传播过程。
具体实施方式
下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
如图1所示，已知有m个位移传感器测点埋设在坝体内部，每个传感器获取到的位移数据可视作一条时间序列，则搜集到的监测数据可以表示为m个时间序列。这些数据组成时间序列集合S＝{s1,s2,...,si}(i＝1...m)。其中
si＝{X1,X2,...,Xj}(j＝1...n)表示单个站点i的监测数据序列，序列中表示共有n 个时间点，且时间序列中的每个观测值Xj都是一个二维向量，包含横河向位移和顺河向位移(本文仅以顺河向位移为例说明)。将大坝视作由众多空间点组成，大坝连续坝段顺河向位移监测数据可视为由组成大坝的各点的时间序列组成，即全集A，
定义本发明中的问题为：在给定的时刻t(t＜n)，基于有限测点集合S，使用集合S中(1,2,...,t)有限时刻数据，预测全集A(S及其补集CAS)中任意一点在未来任意时刻(t+1,t+2,...,t+F)的对应值。具体而言，模型的输入数据是已知测点的坐标和历史位移数据，输出结果是任意测点未来任意时刻的径向位移预测值。问题示意图可参照表1所示，以表1为例具体说明。表1中横向为各个测点时间序列数据，纵向为各时间节点顺河向位移数据，X值为补集待预测值。以大坝实测数据为例，用随机方式在其中选取4个时间序列数据作为数据测试集，比如测试集随机选取为C4A25IP01，C4A25PL06，C4A29IP01和C4A29PL05，即CAS＝{C4A25IP01C4A25PL06C4A29IP01C4A29PL05}，S为剩余项点的集合，t为2015年11月1日，F为310。
表1 某高拱坝时空预测问题形式化示意表
求解问题为利用S中从2012年6月1到2015年10月31号日的数据预测补集CAS从2015年11月1号到2016年9月6号对应值，即表1中X处数值。图 2表示假设当前时间点为t3，要以测点s1,s2,s3,s4在时间t1,t2,t3上的数据预测t4时刻空间任意一点值。
本发明提供基于时空一体化的高拱坝变形时空序列预测算法的实施例，图3 是高拱坝变形时空序列预测的总体框架图。可以看出，本发明实施例的主要步骤如下：
(1)大坝传感器在监测数据时，由于各种现实因素，实测数据中时常带有随机误差项，假设时空过程可以分解为两部分：Zi(t)＝Mi(t)+ei(t)，式中Mi(t) 代表确定性时空变化，ei(t)代表去除确定趋势后的小规模随机变化并且满足 E(ei(t))＝0，即期望为零。利用简单BP神经网络拟合整体时空趋势；
(2)去除整体时空趋势后得到残差项ei(t)，对残差项进行线性无偏估计，选用时空克里格方法拟合局部时空趋势；
(3)引入门限循环神经网络对大坝各测点时间序列进行预测，预测出相关测点变形值，将预测后的顺河向位移数据作为训练数据优化BP网络。
其中步骤(1)的具体步骤为：
(1.1)确定BP网路拓扑结构：确定神经网络输入向量和目标输出向量，本实验中拟设输入层中神经元数目为4，分别对应三维坐标和时间，如图4所示。设输出层的节点数目为1，对应大坝测点顺河向位移值。网络结构中关于隐含层个数、各隐含层中神经元个数设置的问题，按神经元个数从小到大依次尝试后计算出最佳个数；由于大坝变形过程复杂，为较精确的诠释大坝变形的内在规律，设置BP网络隐含层个数为2，第一层中隐含单元个数分别为5、10、20、30、 40、50、60、70，第二层中隐含层单元个数分别为0、5、10、20，遍历第一层隐含单元个数的不同情况，在遍历的同时对第二层隐含单元个数做相同的事情，即两个循环嵌套。
(1.2)确定网络结构权重系数：神经网络无需深入了解大坝具体受力情况，通过训练网络结构中的系数来达到时空预测的目的。BP网络不同层神经元之间通过权重实现完全连接。BP网络在初始化时以[-1,1]之间的随机数对权重赋值；
(1.3)生成激活后的估计值。激活函数有很多种类，其中一种常用的非线性激活函数是由逻辑函数定义的S形传递函数；
(1.4)多次训练得到理想模型：将输入向量带入神经网络中，经过前向传播和反向传播算法的不断训练，将误差控制在理想范围内，使用Matlab软件中 BP网络预设误差0.001，训练次数为1000，当误差达到预设精度或学习次数大于设定的最大次数，则结束算法。
其中步骤(2)的具体步骤为：
(2.1)残差数据应该尽可能充分，残差项有34*1247，即42398个残差数据，可认为残差数据充分，对于残差数据，两两计算时空距离；
(2.2)寻找一个拟合曲线拟合距离与时空半方差的关系，从而能根据任意时空距离计算出相应的半方差；
(2.3)计算出所有已知点之间的时空半方差rij；
(2.4)对于未知点Zo，计算它到所有已知点Zi的半方差rio；
(2.5)求解代价函数得到最优系数λi，代表未知点的估计值；
(2.6)使用最优系数对已知点的属性值进行加权求和，得到未知点的估计值，结束时空克里格方法。
其中步骤(3)的具体步骤为：
(3.1)将时间序列预测问题转化为监督学习问题：使用Keras深度学习框架实现GRU网络。Keras中的GRU模型将数据分为两类：输入向量和输出向量。应用到时间序列预测问题上，以上一个时间点(t-1)的观察值作为输入，当前时间点(t)的观察值作为输出来实现问题转换的目的。GRU网络在本次时间序列预测中应用多输入GRUs网络。多输入GRUs网络输入向量中不仅包括测点自身顺河向位移时间序列，还包括环境量数据。具体而言，输入变量中包括滞后时间点 (t-1)的观察值，还包括(t-1)时刻环境变量值，如表2所示。表2中WatLU(t-1) 表示相对于时间点(t)滞后一天的上游水位值，单位为m；WatLD(t-1)表示相对于时间点(t)滞后一天的下游水位值，单位为m；Rain(t-1)表示相对于时间点(t)滞后一天的降雨量，单位为mm；Flow(t-1)表示相对于时间点(t)滞后一天的大坝流量值，单位为L/s；
表2 测点C4-A25-PL-02GRU多输入实力表
(3.2)变换时间序列数据，使其平稳。具体的操作流程是先从观测值中去除时间趋势，然后对时间序列建模得到预测值，最后将时间趋势再加回到预测值中。平稳化时间序列有消除时间序列趋势的一般方法是对数据进行差分运算，即 (t)时刻位移值减去(t-1)时刻位移值，得到差分序列；
(3.3)将观测结果标准化；
(3.4)建立GRU神经网络：GRU循环网络包含输入层、隐含层和输出层。在添加GRU层的时候，首先需要确定输入层输入向量，输出层输出向量。初始时，随机生成权重矩阵；
(3.5)对于GRU来说，由前向传播计算出重置门r、更新门z、当前缓存记忆内容g、当前时间步的最终记忆h这四个向量的值。计算过程依次为：
rt＝σ(Wrxt+Urht-1)
zt＝σ(Wzxt+Uzht-1)
gt＝tanh(Wg(ht-1⊙rt)+xtUg)
ht＝yt＝ht-1⊙(1-zt)+zt⊙gt
以上式中，σ为Sigmoid激活函数(即S型函数)；xt表示时刻t的输入；ht-1保存t-1时刻的记忆；W和U分别代表GRU网络中的权重系数；
(3.6)通过BPTT计算误差项δj值：利用BPTT算法反向计算每个神经元的误差项δj本质上是误差函数E对神经元j输入量netj的偏导。BPTT算法将第l 层t时刻的误差项值沿两个方向传播：一是沿时间的反向传播，即从当前t时刻开始，计算每个时刻的误差项，与权重系数U相关；另一传播方向是是将误差项向上一层传播，与权重系数W相关。BPTT具体应用到GRU网络上的流程示意如下图5所示，每一步残差计算的过程将在表3中得到具体阐述。表3中δ1、δ2为前向传播输出值与期望值之间的误差；
表3
(3.7)计算每个权重系数的梯度：GRU网络在训练之前对权重w赋予较小的随机值，根据BPTT算法不断训练权重。网络的权重通过纠正神经元的输出变量和理想输出之间的误差来改变并且及时调整。不断重复执行这一过程，使得每训练依次网络的权重就更新一次，直到最后所有训练数据的输出项达到理想的精度，结束训练。