土石坝由于其投资少、施工简单、易就地取材等优点，已成为中国坝工建设中的主要坝型之一，对其进行位移成因分析及预测可更好地掌握坝体变形规律，了解坝体变化趋势，从而实现对坝体的安全监控，对大坝的安全运行与预警具有重要意义^[1–3]。目前，大坝变形预测最广泛的方法是传统数学模型（统计模型、确定性模型和混合模型）。该模型简单易行，但各因子间存在多重相关性和非线性映射不强等问题，从而使得模型拟合精度难以保证，预测结果可能产生偏差^[4]。时间序列分析方法可在缺少环境量监测数据的情况下，仅通过分析效应量监测数据的历史变化过程，寻找效应量的变化特征、趋势和发展规律，其中自回归移动平均模型（ARMA）和差分自回归移动平均模型（ARIMA）在大坝位移预测研究中取得了一定成果^[5–7]，但两模型并不适用于具有周期性的监测数据，有一定的局限性。

众所周知，土石坝主要由松散土体构成，具有一定的弹塑性，其位移主要受到水压、温度和时效等因子的共同作用。通常运行期大坝的水压因子与温度因子具有一定周期性变化，而时效因子变化为趋势性的，所以大坝的位移变化通常呈一定的周期性和趋势性。季节差分自回归移动平均模型（SARIMA）可适用于此类数据的预测^[8–10]，但该模型通过普通差分消除原序列的趋势性可能会出现过差分的问题，且普通差分形式更适合描述短期状态或非均衡状态，不适合描述长期状态，且会导致原始信息丢失过多，预测结果出现偏差^[11]。因此，首先采用HP滤波将实测数据序列分解为趋势项和周期项，这可消除SARIMA模型中普通差分导致的信息丢失问题；然后，对趋势项建立MLR预测模型，对周期项建立SARIMA预测模型；最后，将以上两模型的预测结果组合分析，即MLR–SARIMA预测模型。MLR–SARIMA模型突出了MLR模型在趋势性数据上的预测优势和SARIMA模型在周期性数据上的预测优势，且该模型仅从实测位移数据分析预测，可适用于缺少环境量数据的情况。

1 HP滤波

土石坝位移变化通常呈一定的周期性和趋势性，HP滤波可有效地将其分解为周期项和趋势项两部分，进而分别进行建模分析。HP滤波是由Hodrick和Prescott于1980年提出的，在宏观经济学中应用HP滤波研究序列中的趋势性和周期性已十分普遍^[12–13]，之后被学者广泛应用于各行业数据序列的趋势性研究^[14–16]，其主要原理为：

设 ${Y_t}$ 为包含趋势成分 $Y_t^{\rm{T}}$ 和周期成分 $Y_t^{\rm{C}}$ 的时间序列，可以表示为：

$ {Y_t} = Y_t^{\rm{T}} + Y_t^{\rm{C}}\;\;\left( {t = 1,2,\;\cdots\;,T} \right) $

(1)

HP滤波计算即从 ${Y_t}$ 中将 $Y_t^{\rm{T}}$ 分离出来。一般地， $Y_t^{\rm{T}}$ 常被定义为式（2）最小化问题的解：

$\min \sum\limits_{t = 1}^T{[{{({Y_t} - Y_t^{\rm{T}})}^2} + \lambda {{[(Y_{t + 1}^{\rm{T}} - Y_t^{\rm{T}}) - (Y_t^{\rm{T}} - Y_{t - 1}^{\rm{T}})]}^2}]} $

(2)

式中， $\lambda$ 为平滑参数，一般地，年度、季度和月度数据的 $\lambda$ 取值分别取100、1 600和14 400。

通过对式（2）中不同的 $Y_t^{\rm{T}}$ 求偏导，并令导数为0，可得1阶条件：

${Y_t} = Y_t^{\rm{T}} + \lambda (Y_{t - 2}^{\rm{T}} - 4Y_{t - 1}^{\rm{T}} + 6Y_t^{\rm{T}} - 4Y_{t + 1}^{\rm{T}} + Y_t^{\rm{T}})$

(3)

对于未知数 $Y_{ - 1}^{\rm{T}}$ 、 $Y_0^{\rm{T}}$ 、 $Y_{t + 1}^{\rm{T}}$ 、 $Y_{t + 2}^{\rm{T}}$ 无法通过1阶条件求解，此时假设趋势成分 $Y_t^{\rm{T}}$ 在观测值的最初和最终变化较慢，即设 $Y_t^{\rm{T}} - Y_{t - 1}^{\rm{T}} = Y_{t + 1}^{\rm{T}} - Y_t^{\rm{T}}$ （t=–1，0，t+1，t+2），从而求得其解。

将求得的 $Y_{ - 1}^{\rm{T}}$ 、 $Y_0^{\rm{T}}$ 、 $Y_{t + 1}^{\rm{T}}$ 、 $Y_{t + 2}^{\rm{T}}$ 代入式（3），并用矩阵表示方程组为：

$ \left[ {\begin{array}{*{20}{c}} {{Y_1}} \\ {{Y_2}} \\ {{Y_3}} \\ \vdots \\ {{Y_{t - 2}}} \\ {{Y_{t - 1}}} \\ {{Y_t}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {Y_1^{\rm{T}}} \\ {Y_2^{\rm{T}}} \\ {Y_3^{\rm{T}}} \\ \vdots \\ {Y_{t - 2}^{\rm{T}}} \\ {Y_{t - 1}^{\rm{T}}} \\ {Y_t^{\rm{T}}} \end{array}} \right] + \lambda \left[ {\begin{array}{*{20}{c}} 1&{ - 2}&1&0& \cdots &\cdots&\cdots&0 \\ { - 2}&5&{ - 4}&1&0& \cdots &\cdots&0 \\ 1&{ - 4}&6&{ - 4}&1&0& \cdots &0 \\ \vdots &\vdots&\vdots& \vdots &\vdots&\vdots&\vdots&\vdots\\ 0& \cdots &0&1&{ - 4}&6&{ - 4}&1 \\ 0&\cdots& \cdots &0&1&{ - 4}&5&{ - 2} \\ 0&\cdots & \cdots & \cdots &0&1&{ - 2}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {Y_1^{\rm{T}}} \\ {Y_2^{\rm{T}}} \\ {Y_3^{\rm{T}}} \\ \vdots \\ {Y_{t - 2}^{\rm{T}}} \\ {Y_{t - 1}^{\rm{T}}} \\ {Y_t^{\rm{T}}} \end{array}} \right] $

(4)

整理式（4）可得 ${Y_t} = Y_t^{\rm{T}} + \lambda {{G}}Y_t^{\rm{T}}$ ，G为系数矩阵，则可求得趋势成分：

$ Y_t^{\rm{T}} = {(\lambda {{G}} + 1)^{ - 1}}{Y_t} $

(5)

2 MLR–SARIMA模型 2.1 建模步骤

土石坝位移变化中的周期成分和趋势成分有各自的变化规律，MLR–SARIMA模型将MLR模型和SARIMA模型组合，可发挥MLR模型在趋势性预测和SARIMA模型在周期性预测方面的优势，建模步骤如下：

1）实测数据周期性与趋势性检验；

2）采用HP滤波计算，将实测数据序列分解为周期项和趋势项；

3）采用SARIMA模型对周期项进行建模；

4）采用MLR模型对趋势项进行建模；

5）将以上两模型的预测结果组合；

6）预测结果分析。

2.2 SARIMA模型 2.2.1 模型定义

SARIMA模型是一种适用于具有非平稳性和季节性时间序列的预测方法，通过普通差分和季节差分消除时间序列的非平稳性和季节性，将原时间序列转化为平稳且非季节性的序列，再进行ARMA建模分析，至今已广泛应用于各领域，模型的定义为^[17]：

$\varPhi (B)U({B^S}){\nabla ^d}\nabla _S^D{X_t} = \varTheta (B)V({B^S}){\alpha _t}$

(6)

其中：

$\varPhi (B) = 1 - {\phi _1}B - {\phi _2}{B^2} - \cdots - {\phi _p}{B^p}$

(7)

$U({B^S}) = 1 - {\varGamma _1}{B^S} - {\varGamma _2}{B^{2S}} - \cdots - {\varGamma _k}{B^{kS}}$

(8)

${{{\nabla}} ^d} = {(1 - B)^d}$

(9)

$\nabla _S^D = {(1 - {B^S})^D}$

(10)

$\varTheta (B) = 1 - {\theta _1}B - {\theta _2}{B^2} - \cdots - {\theta _q}{B^q}$

(11)

$V({B^S}) = 1 - {H_1}{B^S} - {H_2}{B^{2S}} - \cdots - {H_m}{B^{mS}}$

(12)

式中： $\varPhi $ 和U分别为非季节和季节AR模型的自回归算子； ${\phi _1},{\phi _2}, \cdots ,{\phi _p}$ 和 ${\varGamma _1},{\varGamma _2}, \cdots ,{\varGamma _k}$ 分别为 $\varPhi$ 和U回归算子的权重系数，p和k分别为 $\varPhi$ 和U回归阶数； ${\nabla ^d}$ 为常规差分算子； $\nabla _S^D$ 为季节差分算子；d为常规差分阶数；D为季节差分阶数；B为滞后算子；S为季节模型；B^S为季节滞后算子； $\varTheta$ 和V分别为非季节和季节MA模型的移动平均算子； ${\theta _1},{\theta _2}, \cdots ,{\theta _q}$ 和 ${H_1},{H_2}, \cdots ,{H_m}$ 分别为 $\varTheta$ 和V移动平均算子的权重系数，q和m分别为MA模型的移动平均阶数； ${X_t}$ 为实测序列； ${\alpha _t}$ 为均值为0、方差为 ${\sigma ^{\rm{2}}}$ 的白噪声序列。

2.2.2 模型识别

SARIMA模型识别是根据序列的自相关函数和偏自相关函数确定模型阶数p与q。通常模型阶数不唯一，需要根据赤池信息准则（Akaike’s information criterion，AIC）和施瓦兹信息准则（Schwartz Bayesian information criterion，BIC）进行模型识别，AIC和BIC的值越小表示模型越简洁；同时，根据德宾沃森统计量（Durbin Watson statistics，DW）范围准则，DW值一般应介于1.8～2.1之间，说明模型残差为正态分布且不相关^[18]。

$AIC = 2k + \ln\; \frac{{RS\!S}}{{{n_{\rm{obs}}}}}\quad\quad$

(13)

$BIC = {n_{\rm{obs}}}\ln\;{\sigma ^{\rm{2}}} + k\ln\;{n_{\rm{obs}}}$

(14)

式中，k为自由参数个数，n_obs为实测序列的总个数，RSS（residual sum of squared）为残差平方和， $\sigma^2$ 为残差的方差。

2.2.3 模型参数估计与检验

SARIMA模型采用最小二乘法进行参数估计，即真实值与估计值的残差平方和达到最小。同时，需要对方程作F检验，对变量作t检验，判断各自的显著性：若显著则符合要求；反之，则需重新进行参数估计。参数估计完成后需对模型残差进行检验，若残差为平稳的白噪声，说明残差中已无可用信息，检验通过；反之，需重新建立模型。

2.2.4 模型预测

通常将原始数据分为训练数据和测试数据，通过对训练数据进行分析建模，再用不同模型进行预测，将预测值与测试数据进行对比，以检验不同模型的优劣性。均方根误差（root mean square error，RMSE）和平均绝对误差百分比（mean absolute percentage error，MAPE）是最常用于评价预测模型优劣的指标。然而，MAPE指标有一点不足，当实测值很小时，即使实测值与预测值差异很小，也可能导致计算的MAPE系数很大。因此，可使用调整的平均绝对误差百分比（adapted mean absolute percentage error，AMAPE）^[19]。

$RMS\!E = \sqrt {\frac{{\displaystyle\sum\limits_{t = 1}^n {{{{\rm{(}}{{\hat {\textit{z}}}_t} - {{\textit{z}}_t}{\rm{)}}}^2}} }}{n}} $

(15)

$MAPE =\Bigg(\frac{1}{n}\sum\limits_{t = 1}^n {\left(\frac{{|{{\hat {\textit{z}}}_t} - {{\textit{z}}_t}|}}{{{{\textit{z}}_t}}}\right)}\Bigg)\times100\% $

(16)

$AMAPE =\left(\frac{1}{n}\sum\limits_{t = 1}^n {\left(\frac{{|{{\hat {\textit{z}}}_t} - {{\textit{z}}_t}|}}{{\dfrac{1}{n}\displaystyle\sum\limits_{t = 1}^n {{{\textit{z}}_t}} }}\right)}\right)\times100\%$

(17)

式中，t为时间，n为预测序列样本个数， ${{\textit{z}}_t}$ 为t时刻实测值， ${\hat {\textit{z}}_t}$ 为t时刻预测值。

2.3 MLR模型

多元线性回归分析模型是利用一个线性函数刻画多个自变量与某因变量之间的关系。由于时效因子的影响，土石坝位移通常会产生一定的趋势性，顾冲时等^[20]给出了时效因子的多种数学形式。可根据多元线性回归分析建立土石坝位移趋势项与时效因子的关系。设某因变量Y与k个自变量X_j（j=0,1, 2， $\cdots, \;k$ ）有响应关系，其多元线性回归模型为^[21]：

${Y_i} = {\beta _0} + {\beta _1}{X_1} + {\beta _2}{X_2} + \cdots + {\beta _k}{X_k} + \varepsilon $

(18)

式中： ${\beta _j}$ 为回归系数，表示在其他回归变量给定的条件下，第j个回归变量增加一个单位时，对从属变量平均响应的影响效果；ε为误差项，ε–N（0， $\sigma^2 $ ）。

MLR模型的参数估计与检验、预测均与SARIMA模型方法相同，这里不再赘述。

3 应用实例 3.1 工程概况

某枢纽工程由土石坝、溢洪道、左右岸泄洪洞、引水洞等组成。该土石坝最大坝高101.80 m，坝顶长度297.36 m，防浪墙顶高程715.30 m，心墙顶部高程710.00 m。水库校核洪水位708.80 m，正常高水位704.00 m，总库容5.21×10⁸ m³。选取水平位移测点D9从2000年1月到2007年12月的实测序列进行分析，将2000年到2006年的数据作为训练数据进行建模分析，将2007年的数据作为测试数据进行预测分析。D9测点的纵向桩号为坝下0+007.80 m，横向桩号为0+184.00 m，高程为708.96 m，约位于下游坝坡中上部位，该测点数据完整可靠，可反映土石坝位移过程的一般变化规律，具有一定代表性，其过程线见图1。

由图1可知，D9测点的水平位移具有一定的趋势性和周期性。趋势性位移主要来自时效因子的影响，由于水库中水位始终不能低于死水位，坝体会持续承受来自于水的推力，使坝体后移。但随着时间的变化，坝体的固结过程会使其变得愈加坚硬密实，水对坝体的作用力相对减小导致坝体后移的速率也逐渐减小。周期性位移主要来自于库水位的年周期性变化，土石坝由松散的土体构成，具有一定弹塑性，所以当上游水位升高时，水对坝体的推力较大，坝体后移，水位降低时，推力相对减小，坝体则恢复一部分位移。同时，由于热胀冷缩的原因，温度的年周期性变化也影响坝体水平位移的周期性。

3.2 MLR–SARIMA建模 3.2.1 HP滤波分解

针对土石坝水平位移的变化特征，采用HP滤波将实测数据序列分解为周期项和趋势项两部分，其结果见图2和3。

由图2、3可知：经HP滤波计算后，实测数据序列被分解为具有一定周期性和趋势性的两组序列。周期项均值为0，呈稳定的年周期性波动，最大值通常出现在每年11月至次年1月之间，最小值通常出现在每年5月至7月之间，年最大变幅约为7.00 mm；趋势项变化呈平缓上升的趋势，其位移从52.63逐渐增至59.71 mm，年最大变幅从1.42逐渐降低至0.51 mm。由以上分析可知，D9测点的水平位移主要是由周期因子引起的周期性位移，而由时效因子引起的趋势性位移较小，且呈越来越小的趋势，这符合土石坝水平位移的一般变化规律。针对以上两组数据的特征，分别采用SARIMA模型和MLR模型对其进行建模预测。

3.2.2 SARIMA建模

D9测点的周期项均值为0，无趋势性。仅对D9测点的周期项进行季节差分，季节差分后序列的自相关系数与偏自相关系数分布见图4。

由图4可知，经季节差分后周期项的自相关系数和偏自相关系数分布均在滞后12阶处显著，因此可判断其具有滞后12阶的相关性。初选模型有SARIMA(0，0，0)(1，1，1)₁₂、SARIMA(1，0，1)(1，1，1)₁₂、SARIMA(1，0，0)(1，1，1)₁₂和SARIMA(0，0，1)(1，1，1)₁₂。各模型的AIC、BIC和DW值见表1。

由表1可知，以上4个模型的DW统计量均在1.8～2.1之间，说明残差为平稳的白噪声。根据AIC和BIC越小越好、修正R²越大越好的准则，最终所选的模型为SARIMA(0，0，0)(1，1，1)₁₂，该模型的拟合结果见图5。

由图5可知，SARIMA的拟合数据保持了训练数据的周期性变化特征，该模型的拟合准确性高，误差较小，且其误差序列通过白噪声检验，无可用信息，说明该模型拟合能力较好。

3.2.3 MLR建模

由统计模型可知，趋势项的变化与时间的1次项及其对数项有关，其回归模型为：

$\delta = {c_1}\vartheta + {c_2}\ln\; \vartheta $

式中，δ为趋势项估计值，ϑ为时间参数，c₁、c₂为回归参数。

以时间的1次项及其对数项为自变量，对趋势项进行回归拟合，其结果见图6。

由图6可知，MLR模型对趋势项的拟合效果非常好，原始值与估计值几乎重叠，模型相关系数为0.998 6。表2中：常数项与时间1次项均为负数，表明其与拟合值呈负相关；时间对数项为正值，表明其与拟合值呈正相关，3个参数均通过了t检验；模型也通过了F检验。

3.2.4 MLR–SARIMA模型预测

由第3.2.2和3.2.3节可知，SARIMA模型和MLR模型对周期项和趋势项均具有较高的拟合能力。结合两种方法对2007年的12个月实测数据进行预测，结果见表3，变化规律见图7。

由表3可知：D9测点趋势项预测值从59.78到60.20缓慢增长，年最大变幅较小，为0.42，其预测结果符合趋势项的变化规律；周期项的预测值呈两端高、中间低的变化趋势，年最大变幅为7.05，与训练数据的变化趋势和幅度（图5）相似。图7中，预测数据与拟合数据的规律相似，与实测数据的拟合程度较好，能反映出D9测点水平位移的变化规律。表4为模型预测指标。由表4可知：MLR–SARIMA模型的RMSE、MAPE和AMAPE值分别为1.14、1.54和1.56，均小于SARIMA模型，说明该模型预测能力强于单一的SARIMA模型；这是由于该模型分别利用了MLR模型对趋势变化较强的预测能力和SARIMA模型对周期变化较强的预测能力。综上分析，MLR–SARIMA模型的预测值准确性高、误差较小，可进行土石坝位移预测。

4 结　论

针对土石坝水平位移具有趋势性和周期性的变化规律，采用HP滤波将实测序列分解为周期项和趋势项两组时间序列。通过分析两组序列6年的变化特征，研究该土石坝水平位移的成因；同时，对周期项采用SARIMA模型建模预测，对趋势项采用MLR模型建模预测，最后将两组序列的预测值结合并与实测数据进行比较。结果表明，该土石坝的水平位移主要是由周期因子影响的周期性位移，由时效因子影响的趋势性位移已趋于稳定。MLR–SARIMA模型的拟合程度较好，预测精确度较高，误差较小，该模型的RMSE、MAPE和AMAPE值均优于SARIMA模型，说明该模型汲取了SARIMA模型在周期性预测与MLR模型在趋势性预测的优势，可应用于具有周期性和趋势性土石坝位移的分析中。