CN104915555A - 均质土坝溃坝分级预警指标的提取方法 - Google Patents

Abstract

本发明公开了一种均质土坝溃坝分级预警指标的提取方法，采用滑动窗口方法分割测值序列，基于主成分分析（PCA）或核主成分分析算法（KPCA）提取土坝监测物理量滑动窗口矩阵，在85%贡献率的基础上确定主成分个数/主成分矩阵范数随土坝性态的演化。实际溃坝试验表明上述两个指标的增大和突变分别对应了土坝性态恶化和最终溃决，充分说明主成分个数或主成分矩阵范数可以作为土坝安全状况的一、二级预警指标；两者之间可以相互校核，其中一级预警指标对应土坝性态恶化，二级预警指标对应土坝溃决。

Description

均质土坝溃坝分级预警指标的提取方法

技术领域

本发明涉及一种基于PCA和KPCA的均质土坝溃坝分级预警指标的提取方法。

背景技术

目前国内外土石坝溃坝试验主要包括室内小比尺模型坝溃坝试验与实体坝现场溃坝试验，室内小比尺模型坝溃坝试验由于模型与原型比尺差别大，不能反映原型坝实际溃决过程。

监测数据是土坝性态的真实反映，但也具有单个测点信息不完备、包含仪器测量误差、测点代表性不够等问题。而现有的预警指标提取方法都是针对单独测点，没有考虑单个测点在评价整个土坝性态方面的信息不完备性。同时通过单个测点信息处理容易受误差影响。

本发明就是通过综合考虑全部实测资料，这样既能充分利用各个测值序列包含的土坝信息，同时也能避免信息冗余，而且能考虑各个变量之间的相关性。为检验本发明方法的合理性。本发明通过对某均质土坝实体试验监测数据分析进行检验。

采用主成分提取土坝性态演化过程中的综合信息，分析比较综合信息随土坝性态演化的关系。首先采用主成分分析(PCA)方法进行主成分分析，考虑到PCA算法不能考虑变量之间的非线性关系。考虑到土坝性态演化期间，各因素之间存在复杂的非线性关系，因此，本发明在PCA方法的基础上应用了非线性降维-核主成分分析算法(KPCA)，有效的解决了监测变量多、各变量之间存在相关性的问题。

在研究过程中，针对KPCA核函数及核参数进行研究，选择一定累积贡献率基础上的主成分数量和主成分矩阵范数作为土坝安全预警指标，实际溃坝试验表明本发明的方法是合理的。

发明内容

目的：为了克服现有技术中存在的不足，本发明提供一种基于主成分分析(PCA)和核主成分分析(KPCA)算法的均质土坝溃坝分级预警指标的提取方法。

技术方案：为解决上述技术问题，本发明采用的技术方案为：

一种均质土坝溃坝分级预警指标的提取方法，包括以下步骤：

步骤A)监测数据初步分析和处理，包括时间序列长度、采样频次统一；

步骤B)拟定合理的窗口大小，根据采样频次和测点数量确定实测资料窗口分割；

步骤C)采用主成分分析算法或核主成分分析算法提取主成分，得到累积贡献率≥85％的主成分个数和主成分矩阵；当测值空间线性时，采用主成分分析算法；当测值空间非线性时，采用核主成分分析算法；

步骤D)采用滑动窗口法：主成分曲线图按照时间段连续列出，绘制由各个窗口提取的主成分个数随时间的演化过程线图；对比实际土坝性态演化情况说明主成分个数可以作为溃坝多级预警指标；

步骤E)计算各个窗口主成分矩阵范数，绘制主成分矩阵范数随时间演化过程线图；实际溃坝试验表明主成分范数变化与土坝性态演化相对应，因此，主成分范数可以作为溃坝分级预警指标；

步骤F)主成分个数和主成分矩阵范数两个指标可以相互校核并同时反映了土坝性态演化，可以提供土坝性态演化的预警指标，即主成分个数/主成分矩阵范数第一个范数峰值作为土坝第一级预警指标，第二突变峰值作为土坝二级监控指标，二级指标对应溃坝状态；

评价研究成果：主成分个数和主成分矩阵范数可以作为溃坝多级预警指标，也说明应用PCA和KPCA提取溃坝预警指标的合理性。

所述的均质土坝溃坝分级预警指标的提取方法，其特征在于：监测数据的参数包括环境量参数、渗流场参数、温度场参数和变形场参数；

监测数据初步分析的研究变量包括：上游水位数据、下游水位数据、渗压水位数据、渗流量、含水率、表面、内部变形数据、坝体(光纤)温度数据。

所述的均质土坝溃坝分级预警指标的提取方法，其特征在于：所述监测数据初步分析包括：

1a)监测数据预处理，统一各测点测值序列长度和采样频次；

2a)过程分析:绘制各监测物理量随时间变化的过程线，了解测值随时间变化情况，通过测值过程线初步分析变量之间的相关性；

3a)相关分析：通过线性相关分析确定变量两两之间的相关程度，为采用主成分分析或核主成分分析提供基础；如果各变量之间是独立的，采用主成分分析，如果变量之间是相关的，则必须采用核主成分分析；

衡量相关密切程度的指标称为相关系数，用R表示，设有a，b两种变量，同步资料共有k对，其相关系数R的计算公式为：

$R=\frac{\underset{i=1}{\overset{k}{\Σ}}(a_i-\stackrel{\‾}{a})(b_i-\stackrel{\‾}{b})}{\sqrt{\underset{i=1}{\overset{k}{\Σ}}{(a_i-\stackrel{\‾}{a})}^2}\sqrt{\underset{i=1}{\overset{k}{\Σ}}{(b_i-\stackrel{\‾}{b})}^2}}---\left(1\right)$

式中：a_i、b_i为两种变量对应的相关点测值；为对应监测资料的均值。

所述的均质土坝溃坝分级预警指标的提取方法，其特征在于：主成分分析算法，包 括以下步骤：

每个滑动窗口由n个p维随机向量X_i的组成，其中p对应测点数，n对应测次数；由于各个监测项目单位不一，且各个变量测值变化范围相差比较大，因此还必须进行相应的标准化；如果样本矩阵中含有逆向指标，需要把它转化成正向指标，可以取其倒数来处理，此时协方差矩阵就变为相关矩阵；

1b)对原始数据进行标准化处理：将X中心化为(就是每列元素减去所在列元素的平均值)：

$\begin{array}{cc}{\stackrel{\‾}{x}}_j=\frac{1}{n}{X}_j^{T}I& (j=1,2,\mathrm{...},p)\end{array}$

标准化后的数据为：

$\begin{array}{cc}{x}_{ij}^{*}=\frac{x_{ij}-{\stackrel{\‾}{x}}_j}{\sqrt{Var\left(x_j\right)}}& (i=1,2,\mathrm{...},n;j=1,2,\mathrm{...},p)\end{array}$

其中， $\sqrt{Var\left(x_j\right)}=\sqrt{\frac{1}{n}{(x_{ij}-{\stackrel{\‾}{x}}_j)}^2}$ 是第j列变量的标准差，Ι是n×1维列向量；

得到中心化的数据(即)，求出标准化后的数据矩阵  $X^{*}={\left(x_{ij}^{*}\right)}_{n\×p},i=1,2,\mathrm{...},n;j=1,2,\mathrm{...},p;$ (注：中心化后协方差矩阵不变)；

2b)计算数据矩阵X^*的协方差矩阵，记为Σ，这里称为X^*的相关矩阵，Σ为p×p维矩阵，Σ的每一个元素Σ(s,t)为：

$\mathrm{\Σ}(s,t)=\frac{\underset{k=1}{\overset{n}{\Σ}}{x}_{ks}^{*}\·x_{kt}^{*}}{n-1},s,t=1,2,\mathrm{...},p$

$r_{st}=\frac{\underset{i=1}{\overset{n}{\Σ}}(x_{is}-{\stackrel{\‾}{x}}_s)(x_{it}-{\stackrel{\‾}{x}}_t)}{(n-1)\sqrt{\underset{i=1}{\overset{n}{\Σ}}{(x_{is}-{\stackrel{\‾}{x}}_s)}^2\underset{i=1}{\overset{n}{\Σ}}{(x_{it}-{\stackrel{\‾}{x}}_t)}^2}}$

此相关矩阵为对称矩阵，所以在以后的计算中，我们取上三角矩阵；

3b)求Σ^up的特征值及相应的特征向量：

在这里把Σ^up的p个特征值记为λ₁≥λ₂≥…≥λ_p≥0，相应的正交特征向量记为ω₁,ω₂,…,ω_p；对协方差矩阵Σ^up进行特征分解，使得Σ^upW＝WΛ；其中，Λ是矩阵Σ^up的特征值矩阵，W是矩阵Σ^up的特征向量：

4b)求主成分： 

用特征值累积贡献率来作为选择标准在已经确定的p个主成分中，合理选取前m个主成分；累积贡献率达到一定的数值(如80％-90％，本发明选用85％)，保证特征提取的主成分中包含原始变量的绝大多数信息。

所述的均质土坝溃坝分级预警指标的提取方法，其特征在于：核主成分分析算法，包括以下步骤：

1c)选择核函数及核参数

根据第一主成分贡献率的大小，从多项式核函数、高斯核函数、神经网络核函数中选择合适核函数及优化核参数；

①高斯径向基核函数：

$\begin{array}{cc}K(x_1,x_2)=\exp \left(-\frac{\left|\right|x_1-x_2|{|}^2}{2{\mathrm{\σ}}^2}\right)& \mathrm{\σ}\≥0\end{array}$

②多项式核函数：

K(x,y)＝(s<x₁,x₂>+b)^d,s,d∈R且c≥0

③Sigmoid核函数：

K(x,y)＝tanh(a(x·y)+b),a,b∈R

2c)输入含有p个监测变量，n个监测样本的数据集X＝[x₁,x₂,…,x_n]，标准化原始指标：将X中心化为(就是每列元素减去所在列元素的平均值)：

$\begin{array}{cc}{\stackrel{\&OverBar;}{x}}_j=\frac{1}{n}{X}_j^{T}I,& (j=1,2,\mathrm{...},p)\end{array}$

标准化后的数据为：

$\begin{array}{cc}{x}_{ij}^{*}=\frac{x_{ij}-{\stackrel{\&OverBar;}{x}}_j}{\sqrt{\mathrm{var}\left(x_j\right)}}& (i=1,2,\mathrm{...},n;j=1,2,\mathrm{...},p)\end{array}$

其中，是第j列变量的标准差，I是n×1维列向量；

3c)应用核函数计算核矩阵，得到n×n维核矩阵

K_ij＝(x_i,x_j),i,j＝1,2,…,n

4c)计算中心化后的核矩阵

$\tilde{K}=K-L_{n}K-{\mathrm{KL}}_n+L_n{\mathrm{KL}}_n$

其中是n×n维元素为的矩阵；

5c)求解矩阵的特征值λ和特征向量β：

$\mathrm{\&lambda;}\mathrm{\&beta;}=\frac{\tilde{K}}{n}\mathrm{\&beta;}$

6c)找出累积贡献率大于一定数值(本发明为85％)的特征值λ和对应的特征向量β，单位化后得到特征向量

7c)核主成分提取：将矩阵在上进行投影；

$Y={\stackrel{\&OverBar;}{\mathrm{\&beta;}}}^T\tilde{K}.$

所述的均质土坝溃坝分级预警指标的提取方法，其特征在于：计算主成分矩阵范数可采用通过行范数、列范数、F范数；

(1) $\left|\right|x|{|}_{\mathrm{\&infin;}}=\underset{i}{\max }\underset{j=1}{\overset{n}{\&Sigma;}}\left|x_{ij}\right|$ (行范数)；

(2) $\left|\right|x|{|}_1=\underset{j}{\max }\underset{i=1}{\overset{m}{\&Sigma;}}\left|x_{ij}\right|$ (列范数)；

(3) $\left|\right|x|{|}_F={(x_1^2+x_2^2+\mathrm{...}+x^2)}^{1/2}$ (F范数)；

计算各个窗口经过主成分获核主成分分析后，通过一定累积贡献率得到的矩阵范数值，绘制范数变化过程线图；主成分范数变化与土坝性态演化相对应，因此，主成分范数作为溃坝分级预警指标；即第一个范数峰值作为第一级预警指标，对应大坝性态恶化；第二个范数突变峰值作为二级监控指标，二级指标对应溃坝状态。

有益效果：本发明提供的均质土坝溃坝分级预警指标的提取方法，本发明就是通过综合考虑全部实测资料，这样既能充分利用各个测值序列包含的土坝信息，同时也能避免信息冗余，而且能考虑各个变量之间的相关性。为检验本发明方法的合理性，本发明通过对某均质土坝实体试验监测数据分析进行检验，采用主成分提取土坝性态演化过程中的综合信息，分析比较综合信息随土坝性态演化的关系。首先采用主成分分析(PCA)方法进行主成分分析，考虑到PCA算法不能考虑变量之间的非线性关系，考虑到土坝性态演化期间，各因素之间存在复杂的非线性关系，因此，本发明在PCA方法的基础上应用了非线性降维-核主成分分析算法(KPCA)，有效的解决了监测变量多、各变量之间存在相关性的问题。

附图说明

图1为本发明的主成分分析流程图；

图2为本发明中核主成分分析流程图；

图3为竖向位移变化过程线；

图4为测点M1渗压水位变化过程线；

图5为渗流量变化过程线；

图6为各时段PCA提取主成分累计贡献率(96×63/窗口)；

图7为σ参数对应累计贡献率；

图8为连续滑动窗口时间段对应表；

图9为各时段主成分累积贡献率≥85％对比图(96×63/窗口)；

图10为各时段累积贡献率≥85％对比图(96×28/窗口)；

图11为累积贡献率≥85％对比图(144×63/窗口)；

图12为各时段累积贡献率≥85％对比图(144×28/窗口)；

图13为核主成分矩阵随时间变化过程线。

具体实施方式

下面结合具体实施例对本发明作更进一步的说明。

如图1所示，为本发明中主成分分析流程图。

如图2所示，为本发明中核主成分分析流程图。

以均质土坝溃坝实体试验来验证：

本发明是对某均质土坝实体试验监测数据进行分析研究，以检验本发明方法的合理性。

1工程概况及试验过程

1.1工程概况

该试验水库最大库容10万m³，挡水建筑物为均质粘性土坝，最土坝高9.7m，坝顶长120m、宽3m，平均坝宽27.25m，上下游边坡为1:2.5，坝体土料为中粉质壤土，黏粒含量为11.5％，压实度97％，内聚力C值范围为7.5～39.5kPa。

1.2试验过程

试验过程包括土坝填筑(仪器埋设)、首次蓄水、水位稳定、第二次蓄水和溃决等过程，参见表1。

表1试验过程历时

2测试参数及测点布设

测试参数包括环境量参数、渗流场参数、温度场参数和变形场参数，共布设上游水位测点1个、下游水位测点1个、气温测点1个、库水温测点1个、大气压力测点1个、渗压水位测点13个、渗流量测点1个、坝体含水率测点2个，表面变形测点6个、内部变形测点6个。温度采用分布式光纤测量，空间分辨率1m，坝体表面温度测线20m，库水温测线20m，坝体内部温度测线343.1m。

3监测数据初步分析处理

本次试验研究变量初步选择上游水位数据(1)、下游水位数据(1)、渗压水位数据(13)、渗流量(1)、含水率(2)、表面、内部变形数据(12)、坝体(光纤)温度数据(33)，共计63个监测变量。

3.1监测数据预处理

保证各个测点序列长度和采样频次的统一。因此对于粗大误差数据必须进行数据剔除，对采样频次过低的时间序列采用插值方法补齐相应的插值频率。

3.2过程分析

本发明首先绘制了竖向位移监测数据、渗压水位监测数据、渗流量随时间变化的过程线，由图3-图5可知各测点测值之间存在明显的相关性。

通过单个变量过程线分析可以看出土坝运行期间的状态变化，但是这几个变量中，监测变量的时效变化各有不同，竖向位移监测数据对水位变化十分敏感，测点M1渗流压力的变化在首次蓄水和水位稳定期并不明显，并且对第二次蓄水期水位变化有着明显的滞后现象。渗流量的变化随着水位变化比较敏感，但是仍存在明显的滞后现象。值得注意的是，在土坝竖向位移的同时，坝体位移的变化同时会影响渗流压力的变化，而渗流压力的变化也影响着渗流量的变化，仅仅通过单个监测变量并不能很好的反映土坝的运行状态。下面对三者进行线性相关分析，进一步探讨两两之间的线性相关性，借此分析线性相关在度量两过程相关方面的特点。

3.3相关分析

相关分析法是一种比较简单的分析两个随机变量相关性的计算方法，在土石坝渗流研究方法多有应用，相关分析把变量都看作随机变量，其目的是确定变量之间的相关联系程度如何。衡量相关密切程度的指标称为相关系数，用R表示，设有a，b两种变量，同步资料共有k对，其相关系数R的计算公式为：

$R=\frac{\underset{i=1}{\overset{k}{\&Sigma;}}(a_i-\stackrel{\&OverBar;}{a})(b_i-\stackrel{\&OverBar;}{b})}{\sqrt{\underset{i=1}{\overset{k}{\&Sigma;}}{(a_i-\stackrel{\&OverBar;}{a})}^2}\sqrt{\underset{i=1}{\overset{k}{\&Sigma;}}{(b_i-\stackrel{\&OverBar;}{b})}^2}}---\left(1\right)$

式中：a_i、b_i为两种变量对应的相关点测值；为对应监测资料的均值。

在相关分析中，一般要求同步观测资料系列在10对以上；相关系数R的绝对值越大，两者关系越密切。当R>0时，为正相关，表示b随a的增大而增大；R<0时，为负相关，表示b随a的增大而减小。

接分析结果，选择渗流量和渗流压力作为随机变量，取30个样本，含水率为a，渗流压力为b，通过式(1)计算二者的相关系数：R＝0.357438613，所以渗流压力是随着渗流量的增大而增大的。

同样计算竖向位移和渗流压力的相关系数：R_D1＝0.003596547，R_D2＝0.0119854547，R_D3＝0.0281269917，R_D4＝0.066507065，R_D5＝0.077217307，R_D6＝0.055317066。

通过以上相关性分析，竖向位移、渗流量、渗流压力两两之间存在非常明显的相关性，在主成分分析的基础上必须考虑核主成分分析。

4监测数据主成分分析

监测数据初步分析：试验研究变量选择上游水位数据(1)、下游水位数据(1)、渗压水位数据(13)、渗流量(1)、含水率(2)、表面、内部变形数据(12)、坝体(光纤)温度数据(33)，共计63个监测变量。监测时间跨度为10月16日到11月17日，每30min选择一个监测样本，对监测数据应用主成分算法进行主成分提取。初步拟定窗口大小：(96×63/窗口)，窗口大小解释为63个测试变量，表示为列；96个样本，表示每连续 两天作为一个窗口，一共是96行，如图6所示。

应用PCA对监测数据进行主成分提取，10月16日～10月29日主成分个数一直增加，10月29日～11月10日逐渐下降，在11月17日突变。通过和试验过程对比，随着首次蓄水期，水位稳定期水位变化，主成分个数也在变化，但是在第二期蓄水期主成分个数从11月7日逐渐减少，在土坝溃决时，主成分个数突变。

综上所述，主成分分析是一种较好的特征提取技术，可以提取土坝溃坝预警指标，两个峰值分别对应坝体性态恶化和溃决。

5监测数据的核主成份分析

5.1核函数选择及核参数优化

不同的核函数会导致评价的效果大不相同，因此如何根据具体数据指标选择恰当的核函数是KPCA应用领域遇到的一个重大难题。至今为止，对核函数的研究大多是在核函数的性质和核函数的构造上面，但是还没有比较完善的理论可以针对具体数据选取恰当的核函数和核参数，使得数据分析可以更有效，这是我们还需要进一步研究的问题。

当选取核函数解决数据问题时，通常采用的方法有：一是利用Cross-Validation方法，选取不同的核函数，找出结果误差最小的核函数，这个称之为有效核函数；二是应用专家的通过多次实验而给定的核函数；三是利用由Smits等人提出的混合核函数方法。

为了选择合适的核函数及合适的核参数，根据本次监测变量包括6个测点(D1-D6)的表面水平位移、表面竖向位移；14个坝体渗流压力测点A1-A6，B7-B14的倾度数据(M)，渗流压力数据(B)、温度数据(TF)、含水率(msvl)；上游水位监测数据(UPWL)，渗流量数据(Q)、雨量数据(PR)，共计63个监测变量；每30分钟选择一个测试样本。这样，选择10月16日、17日两天作为测试样本，进行主成分提取。选择三种常用核函数对比分析，通过第一主成分累积贡献率的大小，选择最优核函数及核参数。

5.1.1多项式核函数 

Κ(x,y)＝(s<x₁,x₂>+b)^d，s，d∈R且c≥0

(1)当s＝0.05，c＝0，d＝3时，第一主成分的贡献率为75.71％。

(2)当s、c保持不变，d＝2时，第一主成分的贡献率为76.80％。

(3)当s、c保持不变，d＝3时，第一主成分的贡献率为78.68％。

(4)当s、c保持不变，当d＝4时，第一主成分的贡献率为79.60％。

MATLAB参数循环试值，求得多项式核函数第一主成分贡献率：

当s＝0.06、c＝0、d＝5时，第一主成分的贡献率最大为85.71％

5.1.2神经网络核函数

K(x,y)＝tanh(a(x·y)+b)，a，b∈R

经过MATLAB参数循环试值：

当a＝0.01，b＝1时，第一主成分贡献率最大为84.40％。

5.1.3高斯径向基核函数

$K(x_1,x_2)=\exp (-\frac{\left|\right|x_1-x_2|{|}^2}{2{\mathrm{\&sigma;}}^2}),\mathrm{\&sigma;}\&GreaterEqual;0$

参数序号对应的σ值

如图7，因为由Matlab7.0编程计算的高斯核函数中σ的不同取值对应的第一主成分贡献率可知，当σ≥75时，已经趋近于稳定值87.5％，故σ的取值定为75。

通过比较3中核函数，第一主成分贡献率高斯函数最大，且高斯函数仅一个参数，计算时稳定性更好。因此，本发明选择高斯径向基核函数进行核主成分提取，高斯核函数公式为：

$K(x_1,x_2)=\exp \left(-\frac{\left|\right|x_1-x_2|{|}^2}{2\&times;{75}^2}\right)---\left(2\right)$

5.2提取主成分 

窗口选择类别1：(96×63/窗口)，窗口大小解释为63个测试变量，表示为列；96个样本，表示每连续两天作为一个窗口，一共是96行；连续的32个窗口，相邻窗口包含相同一天的监测信息，连续滑动窗口如图8。

基于监测数据的每个窗口核主成分提取的具体步骤如下：

(1)输入含有63个监测变量，96个监测样本的原始矩阵x＝[x₁,x₂,…,x₆₃]_96×63，其中，x_i是列向量，x_(i)是行向量。将X中心化为(就是每列元素减去所在列元素的平均值)：

$\begin{array}{cc}{\stackrel{\&OverBar;}{x}}_j=\frac{1}{96}{X}_j^{T}I,& (j=1,2,\mathrm{...},63)\end{array}$

其中，I为96×1的单位向量。标准化后的数据为：

$\begin{array}{cc}{x}_{ij}^{*}=\frac{x_{ij}-{\stackrel{\&OverBar;}{x}}_j}{\sqrt{Var\left(x_j\right)}}& (i=1,2,\mathrm{...},96;j=1,2,\mathrm{...},63)\end{array}$

其中，是第j列变量的标准差，I是n×1维列向量。

(2)应用高斯径向基核函数(4-1)进行内积计算，得到96×96维核矩阵Κ＝(k_st)_96×96。

(3)矩阵Κ＝(κ_st)_96×96进行中心化处理：

$\tilde{K}=K-L_{96}K-{\mathrm{KL}}_{96}+L_{96}{\mathrm{KL}}_{96}---\left(4\right)$

其中：矩阵L₉₆为方矩阵

(4)求解矩阵的特征值λ和特征向量β：

$\mathrm{\&lambda;}\mathrm{\&beta;}=\frac{\tilde{K}}{96}\mathrm{\&beta;}---\left(5\right)$

其中：λ＝[λ₁，λ₂，…，λ₉₆]是核矩阵的特征值矩阵，β是核矩阵的特征值矩阵λ对应的特征向量矩阵。

(5)找出累积贡献率大于85％的特征值λ和对应的特征向量β，单位化特征向量

(6)核主成分提取：将矩阵在上进行投影。

$Y={\widetilde{\mathrm{\&beta;}}}^T\tilde{K}---(4-5)$

5.3窗口敏感性分析

5.3.1窗口(96×63)

绘制各时段累计贡献率对比图(96×63/窗口)和各时段主成分曲线图。主成分个数变化如图9所示，主成分曲线图按照时间段连续列出，称之为滑动窗口。其中，KPCA1-KPCA7分别代表提取的第一主成分、第二主成分、……、第七主成分的单个贡献率。

通过图9，可以看出在10月17日到11月11日，主成分个数先增加后减小，在10月30日主成分个数达到最大，在11月12日开始主成分个数明显减少，在11月17日突增。

通过主成分曲线图：10月17日到11月11日，主成分个数先增加后减小，在10 月30日主成分个数达到最大，在11月12日开始主成分个数明显减少，在11月17日突增。

表2试验时段水位变化

结合表2对主成分个数变化进行说明：首次蓄水期10月14日～10月25日，水位增加3.88m，主成分个数也逐渐增加；水位稳定期10月25日～11月7日，在渗流作用下水位下降0.2m，此阶段主成分个数先增加后减少，在10月30日达到最大值；11月7日第二期蓄水，水位继续增加，11月17日溃坝，溃坝时段为1h50min，但是此阶段主成分个数明显减少，在11月17日溃坝时发生突变。

通过主成分曲线图：从滑动窗口的变化，可以很直观的看出，在10月23号主成分曲线图波动性增大，初步认为土坝的运行状态开始不稳定；在10月30号主成分曲线波动性最大，对应的主成分个数也最大，可以理解为此时的土坝状态变化剧烈。在11月11号主成分曲线图趋于稳定，土坝运行稳定；11月17日主成分变化发生突变。

5.3.2窗口(96×28)

基于以上过程，再次选择窗口类别，在监测变量上进行修正，选择窗口类别2进行计算。窗口选择类别2：(96×28/窗口)，窗口大小解释为28个测试变量，96个样本，32个窗口，同样得到各时段主成分个数变化图以及可视化窗口，如图10所示。

通过主成分个数变化图：10月18日主成分开始变化，10月27日和30日主成分个数最多，11月10日主成分个数减少，11月17日主成分个数突变。

通过主成分曲线图：10月17日主成分开始变化，10月27日和30日主成分个数最多，11月12日主成分个数减少，11月17日主成分个数突变。

对比窗口(96×63/窗口)和(96×28/窗口)主成分个数变化的结果：

监测变量是63的窗口：10月17日到11月11日，主成分个数先增加后减小，在10月30日主成分个数达到最大，在11月12日开始主成分个数明显减少，在11月17日突增。

监测变量是28的窗口：10月18日主成分开始变化，10月27日和30日主成分个数最多，11月10日主成分个数减少，11月17日主成分个数突变。

5.3.3窗口(144×63)

对比主成分窗口(96×63/窗口)和(96×28/窗口)，两个窗口类别的主成分个数变化在反映土坝溃坝时段中，有同样的趋势，但是在土坝状态变化时段判别中多有差异。

为进一步验证改方法的合理性，在监测样本上进行修订，选择每相邻三天作为一个窗口进行计算、分析。窗口选择类别3：监测变量63个；样本数144个；连续窗口16个，记为(144×63/窗口)。

应用MATLAB对连续16个窗口(144×63/窗口)进行主成分提取，得到累积贡献率≥85％的主成分个数和主成分矩阵，各主成分累积贡献率见附录C(表4.5各时段累 积贡献率(144×63/窗口))。绘制主成分个数变化图(144×63/窗口)和各时段主成分曲线图，如图11所示。

通过主成分个数变化图：10月18日主成分开始变化，10月30日主成分个数最多，11月5日主成分个数减少，11月17日主成分个数突变。

通过主成分曲线图：10月26日主成分开始变化，10月30日主成分波动最剧烈，11月5日主成分逐渐平缓，11月17日主成分波动剧烈。

5.3.4窗口(144×28)

窗口类别4：监测变量28个；样本数144个；连续窗口15个，记为(144×28/窗口)。

应用MATLAB对连续15个窗口(144×28/窗口)进行主成分提取，得到累积贡献率≥85％的主成分个数和主成分矩阵，各主成分累积贡献率见附录C(表4.6各时段累积贡献率(144×28/窗口))。绘制主成分个数变化图(144×28/窗口)和各时段核主成分曲线图，如图12所示。

通过主成分个数变化图：10月22日开始主成分有明显变化，10月30日主成分个数最多，11月11日主成分个数减少，11月17日主成分个数突变。

通过主成分曲线图：10月24日主成分曲线图波动性增大，土坝的运行状态开始不稳定；10月27、28、29日和11月2、3、4日主成分波动均剧烈，土坝运行状态极不稳定；11月11日主成分曲线图趋于稳定，11月17日主成分曲线图突然剧烈变化。

5.3.5基于主成分个数的预警指标评价

综合上述四个类别窗口进行分析，对四个窗口类别反馈的溃坝时段信息和溃坝预警指标进行汇总，综合分析每个窗口反映的土坝运行状态。

表3四个窗口类别反馈的溃坝信息

可以总结出，在四个窗口都能反映土坝性态的演化，两级预警指标的变化与实际试验结果吻合，对窗口大小并不敏感，说明本发明所提出方法是合理的。四个窗口相比较而言，窗口(96×63)整体性和局部性都很好，反映的土坝溃坝时段信息最明确。

6基于核主成分矩阵范数的预警指标提取

6.1范数计算

本发明选择几种比较常见的矩阵范数进行主成分矩阵的计算。

(1) $\left|\right|x|{|}_{\mathrm{\&infin;}}=\underset{i}{\max }\underset{j=1}{\overset{n}{\&Sigma;}}\left|x_{ij}\right|$ (行范数)；

(2) $\left|\right|x|{|}_1=\underset{j}{\max }\underset{i=1}{\overset{m}{\&Sigma;}}\left|x_{ij}\right|$ (列范数)；

(3) $\left|\right|x|{|}_F={(x_1^2+x_2^2+\mathrm{...}+x^2)}^{1/2}$ (F范数)

本发明运用三种矩阵范数对4.4节核主成分算法计算出的主成分矩阵进行计算，在溃坝预警指标信息方面有所新发现。

6.2基于核主成分矩阵范数的预警指标评价

计算四个窗口类别每个主成分矩阵的范数值，绘制范数曲线图，如图13。通过范数曲线图可以很明显看出主成分范数变化与土坝性态演化相对应，和主成分个数变化及滑动窗口反映的溃坝时段信息相吻合，并且和溃坝监测过程十分接近，而且范数的变化对具体范数计算方法鲁棒，因此可以断定主成分范数可以作为溃坝分级预警指标，即第一个范数峰值作为第一级预警指标，第二个突变峰值作为二级监控指标，二级指标对应溃坝状态。

可见核主成分范数变化与主成分个数变化趋势一致，且与土坝性态演化过程相吻合。说明基于核主成分矩阵范数的预警指标对范数具体算法不敏感，可以作为土坝溃坝预警指标。

7总结

通过应用KPCA对监测数据进行主成分提取，绘制主成分个数对比图和主成分曲线图，对比二者反映的土坝性态演化与溃坝信息，说明主成分个数可以作为溃坝多级预警指标。通过行范数、列范数、F范数计算可知，主成分范数对范数计算方法不敏感，充分说明主成分范数可以作为土坝预警指标，即10月30日主成分范数的第一峰值对应溃坝一级预警指标，11月17日主成分峰值对应溃坝二级预警指标，达到二级指标时土坝对应溃坝状态。这一结论和溃坝试验过程吻合，进一步说明了主成分个数和主成分矩阵范数可以作为溃坝多级预警指标，也说明应用PCA和KPCA提取溃坝预警指标是合理性。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (6)

1.一种均质土坝溃坝分级预警指标的提取方法，包括以下步骤：

步骤A)监测数据初步分析和处理，包括时间序列长度、采样频次统一；

步骤B)拟定合理的窗口大小，根据采样频次和测点数量确定实测资料窗口分割；

步骤C)采用主成分分析算法或核主成分分析算法提取主成分，得到累积贡献率≥85％的主成分个数和主成分矩阵；当测值空间线性时，采用主成分分析算法；当测值空间非线性时，采用核主成分分析算法；

步骤D)采用滑动窗口法：主成分曲线图按照时间段连续列出，绘制由各个窗口提取的主成分个数随时间的演化过程线图；对比实际土坝性态演化情况说明主成分个数可以作为溃坝多级预警指标；

步骤E)计算各个窗口主成分矩阵范数，绘制主成分矩阵范数随时间演化过程线图；主成分范数变化与土坝性态演化相对应，因此，主成分范数可以作为溃坝分级预警指标；

步骤F)主成分个数和主成分矩阵范数两个指标可以相互校核并同时反映了土坝性态演化，可以提供土坝性态演化的预警指标，即主成分个数/主成分矩阵范数第一个范数峰值作为土坝第一级预警指标，第二突变峰值作为土坝二级监控指标，二级指标对应溃坝状态；

评价研究成果：主成分个数和主成分矩阵范数可以作为溃坝多级预警指标，也说明应用PCA和KPCA提取溃坝预警指标的合理性。

2.根据权利要求1所述的均质土坝溃坝分级预警指标的提取方法，其特征在于：监测数据的参数包括环境量参数、渗流场参数、温度场参数和变形场参数；

监测数据初步分析的研究变量包括：上游水位数据、下游水位数据、渗压水位数据、渗流量、含水率、表面、内部变形数据、坝体(光纤)温度数据。

3.根据权利要求2所述的均质土坝溃坝分级预警指标的提取方法，其特征在于：所述监测数据初步分析包括：

1a)监测数据预处理，统一各测点测值序列长度和采样频次；

2a)过程分析:绘制各监测物理量随时间变化的过程线，了解测值随时间变化情况，通过测值过程线初步分析变量之间的相关性；

3a)相关分析：通过线性相关分析确定变量两两之间的相关程度，为采用主成分分析或核主成分分析提供基础；如果各变量之间是独立的，采用主成分分析，如果变量之间是相关的，则必须采用核主成分分析；

衡量相关密切程度的指标称为相关系数，用R表示，设有a，b两种变量，同步资料共有k对，其相关系数R的计算公式为：

$R=\frac{\underset{i=1}{\overset{k}{\&Sigma;}}(a_i-\stackrel{\&OverBar;}{a})(b_i-\stackrel{\&OverBar;}{b})}{\sqrt{\underset{i=1}{\overset{k}{\&Sigma;}}{(a_i-\stackrel{\&OverBar;}{a})}^2}\sqrt{\underset{i=1}{\overset{k}{\&Sigma;}}{(b_i-\stackrel{\&OverBar;}{b})}^2}}---\left(1\right)$

式中：a_i、b_i为两种变量对应的相关点测值；为对应监测资料的均值。

4.根据权利要求1所述的均质土坝溃坝分级预警指标的提取方法，其特征在于：主成分分析算法，包括以下步骤：

每个滑动窗口由n个p维随机向量X_i的组成，其中p对应测点数，n对应测次数；由于各个监测项目单位不一，且各个变量测值变化范围相差比较大，因此还必须进行相应的标准化；如果样本矩阵中含有逆向指标，需要把它转化成正向指标，可以取其倒数来处理，此时协方差矩阵就变为相关矩阵；

1b)对原始数据进行标准化处理：将X中心化为就是每列元素减去所在列元素的平均值：

${\stackrel{\&OverBar;}{x}}_j=\frac{1}{n}{X}_j^{T}I,(j=1,2,\mathrm{...},p)$

标准化后的数据为：

$x_{ij}^{*}=\frac{x_{ij}-{\stackrel{\&OverBar;}{x}}_j}{\sqrt{Var\left(x_j\right)}},(i=1,2,\mathrm{...},n;j=1,2,\mathrm{...},p)$

其中，是第j列变量的标准差，Ι是n×1维列向量；

得到中心化的数据(即)，求出标准化后的数据矩阵 $X^{*}={\left(X_{ij}^{*}\right)}_{n\&times;p},i=1,2,\mathrm{...},n;j=1,2,\mathrm{...},p;$ 中心化后协方差矩阵不变；

2b)计算数据矩阵X^*的协方差矩阵，记为Σ，这里称为X^*的相关矩阵，Σ为p×p维矩阵，Σ的每一个元素Σ(s,t)为：

$\mathrm{\&Sigma;}(s,t)=\frac{\underset{k=1}{\overset{n}{\&Sigma;}}{x}_{ks}^{*}\&CenterDot;x_{kt}^{*}}{n-1},s,t=1,2,\mathrm{...},p$

$r_{st}=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}(x_{is}-{\stackrel{\&OverBar;}{x}}_s)(x_{it}-{\stackrel{\&OverBar;}{x}}_t)}{(n-1)\sqrt{\underset{i=1}{\overset{n}{\&Sigma;}}{(x_{is}-{\stackrel{\&OverBar;}{x}}_s)}^2\underset{i=1}{\overset{n}{\&Sigma;}}{(x_{it}-{\stackrel{\&OverBar;}{x}}_t)}^2}}$

此相关矩阵为对称矩阵，所以在以后的计算中，我们取上三角矩阵；

3b)求Σ^up的特征值及相应的特征向量：

在这里把Σ^up的p个特征值记为λ₁≥λ₂≥…≥λ_p≥0，相应的正交特征向量记为ω₁,ω₂,…,ω_p；对协方差矩阵Σ^up进行特征分解，使得Σ^upW＝WΛ；其中，Λ是矩阵Σ^up的特征值矩阵，W是矩阵Σ^up的特征向量：

4b)求主成分：

用特征值累积贡献率来作为选择标准在已经确定的p个主成分中，合理选取前m个主成分；累积贡献率达到一定的数值，保证特征提取的主成分中包含原始变量的绝大多数信息。

5.根据权利要求1所述的均质土坝溃坝分级预警指标的提取方法，其特征在于：核主成分分析算法，包括以下步骤：

1c)选择核函数及核参数

根据第一主成分贡献率的大小，从多项式核函数、高斯核函数、神经网络核函数中选择合适核函数及优化核参数；

①高斯径向基核函数：

$K(x_1,x_2)=\exp (-\frac{\left|\right|x_1-x_2|{|}^2}{2{\mathrm{\&sigma;}}^2}),\mathrm{\&sigma;}\&GreaterEqual;0$

②多项式核函数：

K(x,y)＝(s<x₁,x₂>+b)^d,s,d∈R且c≥0

③Sigmoid核函数：

K(x,y)＝tanh(a(x·y)+b),a,b∈R

2c)输入含有p个监测变量，n个监测样本的数据集X＝[x₁,x₂,…,x_n]，标准化原始指标：将X中心化为(就是每列元素减去所在列元素的平均值)：

${\stackrel{\&OverBar;}{x}}_j=\frac{1}{n}{X}_j^{T}I,(j=1,2,\mathrm{...},p)$

标准化后的数据为：

$x_{ij}^{*}=\frac{x_{ij}-{\stackrel{\&OverBar;}{x}}_j}{\sqrt{var\left(x_j\right)}},(i=1,2,\mathrm{...},n;j=1,2,\mathrm{...},p)$

其中，是第j列变量的标准差，I是n×1维列向量；

3c)应用核函数计算核矩阵，得到n×n维核矩阵

K_ij＝(x_i,x_j),i,j＝1,2,…,n

4c)计算中心化后的核矩阵

$\tilde{K}=K-L_{n}K-{\mathrm{KL}}_n+L_n{\mathrm{KL}}_n$

其中是n×n维元素为的矩阵；

5c)求解矩阵的特征值λ和特征向量β：

$\mathrm{\&lambda;}\mathrm{\&beta;}=\frac{\tilde{K}}{n}\mathrm{\&beta;}$

6c)找出累积贡献率大于一定数值，的特征值λ和对应的特征向量β，单位化后得

到特征向量

7c)核主成分提取：将矩阵在上进行投影；

$Y={\stackrel{\&OverBar;}{\mathrm{\&beta;}}}^T\tilde{K}.$

6.根据权利要求1所述的均质土坝溃坝分级预警指标的提取方法，其特征在于：计算主成分矩阵范数可采用通过行范数、列范数、F范数；

$\left(1\right)---\left|\right|x|{|}_{\mathrm{\&infin;}}=\underset{i}{\max }\underset{j=1}{\overset{n}{\&Sigma;}}|x_{ij}|$ (行范数)；

$\left(2\right)---\left|\right|x|{|}_1=\underset{j}{\max }\underset{i=1}{\overset{m}{\&Sigma;}}|x_{ij}|$ (列范数)；

$\left(3\right)---\left|\right|x|{|}_F={(x_1^2+x_2^2+\mathrm{...}+x^2)}^{1/2}$ (F范数)；

计算各个窗口经过主成分获核主成分分析后，通过一定累积贡献率得到的矩阵范数值，绘制范数变化过程线图；主成分范数变化与土坝性态演化相对应，因此，主成分范数作为溃坝分级预警指标；即第一个范数峰值作为第一级预警指标，对应大坝性态恶化；第二个范数突变峰值作为二级监控指标，二级指标对应溃坝状态。