聚类算法是发掘数据的类别标记、揭示数据内在结构及规律的无监督学习算法. 传统的聚类算法主要包括K-means法^[1]、层次聚类法^[2]，皆具备较好的数据处理效果. 近年来，在识别数据典型特征时，一些学者将聚类算法引入到了地震工程研究领域中. Stojadinović等^[3]基于K-means法改进结构震后损伤评估的计算效率. Mariniello等^[4]将遗传算法与K-means法相结合，证实在对结构开展抗震安全评估时，聚类算法可以优化解决方案的可行性. 在以往的抗震研究中，研究人员经常优先研究结构自身优化、震后评估问题，再对结构的抗震安全性能做出合理评价. 鉴于实际震害资料的匮乏，如果得到的震害数据有限，也将很难全面评估结构的抗震效果. 为此，地震易损性的研究受到了学者们的广泛关注.

地震易损性是指在不同地震动激励下，结构可能发生的破坏状态^[5]，研究方法通常分为专家法、经验法及解析型研究法等^[6].目前，在解析型地震易损性的研究中，通常需要增量动力分析法（incremental dynamic analysis，IDA）、多样条带分析法（multiple stripe analysis，MSA）和耐震时程分析法（endurance time analysis，ETA）等^[7]执行大规模的非线性计算. 虽然这些方法简便、稳定且发展成熟，能得到较精确的结构性能评估效果，但是计算成本较高.

为了解决这个问题，一些学者提出以规范标准反应谱法^[8]、应用技术委员会法（applied technology council，ATC）^[9]、改进程序法^[10-11]等筛选地震波子集以降低地震易损性计算成本. 这些方法能大幅度降低计算量，但研究人员仍然面临着所选地震波过多及无法合理选择具有代表性的地震波问题. 此外，传统的地震波筛选方式较难涉及复杂的地震动强度指标（intensity measures，IM），在地震波随机性的问题上也较难考虑全面，从而使得地震易损性的结果可能不够准确.

学者们采用聚类算法来解决上述问题. 例如，K-means、层次聚类法皆可以用于筛选地震波子集以降低易损性计算成本^[12]. 然而，K-means法仅适用于样本分布为近似“超球”状的数值型数据，层次聚类法不能回溯调整已构成的簇结构，具有不可逆性. 并且，两者在直接面对高维数据时也容易陷入“维数灾难”的问题.

谱聚类（spectral clustering，SPC）是基于谱图划分的聚类算法^[13]. 该算法适合对任意呈流行状分布的数据进行聚类，适用于高维数据且具备收敛到全局最优解的特点. 不过，要解决各种情形下的数据聚类问题，已有算法都是不完美的，SPC法也不例外. 传统的SPC算法仅适用于原始各簇数据分布较均衡的情形. 当面对各簇样本数相差悬殊的类不平衡问题时，SPC算法将不再适用.

基于上述分析，为了筛选出考虑复杂多指标且适用于评估所选结构抗震性能的典型地震波，本研究基于SPC算法划分了地震波的相似属性及特征. 同时，考虑到传统SPC算法的缺陷，提出将聚类隶属度矩阵的近似正交约束与SPC算法相结合的平衡化谱聚类算法（balanced spectral clustering，BSPC）.在该算法中，利用基于图Laplacian 的隶属度矩阵识别地震波的归属类别. 以高拱坝工程为代表，构建了体现拱坝动力破坏机制的聚类易损性曲线，探索结构破坏模式及失效性态，以在减少地震波输入的同时保证易损性结果的准确性.

SPC算法通过无向带权图 ${\boldsymbol{G}}{\text{ = }}\left( {V{\text{,}}E{\text{,}}{\boldsymbol{W}}} \right)$表示样本点集 $X = \left\{ {{x_i}} \right\}_{i = 1}^n$及其邻近矩阵，从而达到“簇内相近，簇间分离”的效果^[13]. 其中，V为顶点集， $V = \left\{ {{v_i}} \right\}_{i = 1}^n$；E为连接顶点 ${v_i}$与 ${v_j}$之间边的集合， $E = \left\{ {{e_{ij}}} \right\}$；W为图 ${\boldsymbol{G}}$的权重矩阵， ${\boldsymbol{W}} = {\left[ {{w_{ij}}} \right]_{n \times n}}$.该算法要求图的顶点集V与样本点集X一一对应. 同时，采用一种保留数据局部特征的流行降维方法，即转换图Laplacian矩阵，达到降维效果. 定义图Laplacian矩阵为 $ {\boldsymbol{L}} = {\boldsymbol{D}} - {\boldsymbol{W}} $. 其中，D为度矩阵， ${\boldsymbol{D}} = {\rm{diag}}\;\left( {{{\left[ {{d_{ii}}} \right]}_{n \times n}}} \right)$，d_ii为顶点 ${v_i}$的度， ${d_{ii}} = \displaystyle \sum\nolimits_j {{w_{ij}}} $.传统SPC算法具体可参考文献[13].SPC的优化问题表达式如下：

(1) $ \left.\begin{gathered} \mathop {\min }\limits_{{\boldsymbol{F}} \in {{\bf{R}}^{n \times k}}}{{\rm{tr}}} ({{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{LF}}); \\ {\text{ }}{\rm{s}}.{\rm{t}}.\;{\text{ }}{{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{F}} = {\boldsymbol{I}}. \\ \end{gathered} \right\}$

式中： ${\boldsymbol{I}} \in {{{\bf{R}}}^{n \times n}}$为单位矩阵； $ {\boldsymbol{F}} $为聚类隶属度矩阵， ${{{F}}_{ik}}$表示样本 ${x_i}$属于第 $ k $类的隶属度. 根据隶属度的含义， $ {\boldsymbol{F}} $第 $ i $行的 $ {{\boldsymbol{F}}_{i.}} $应属 $ k - 1 $维单纯形，即 ${{\boldsymbol{F}}_{i.}} \in {{{\varDelta }}_k} = \{ {\boldsymbol{z}}|{z_s} \geqslant 0,\displaystyle \sum\limits_{s = 1}^k {{z_s}} = 1.0,{\boldsymbol{z}} \in {{\bf{R}}^{1 \times k}}\}$.因此，本研究在式（1）的基础上提出的BSPC算法模型表达式为

(2) $ \left. \begin{aligned} \mathop {\min }\limits_{{\boldsymbol{F}} \in {{\bf{R}}^{n \times k}}} & {{\rm{tr}}} \left( {{{\boldsymbol{F}}^{\text{T}}}{{\boldsymbol{L}}_{{\text{sym}}}}{\boldsymbol{F}}} \right). \\ {\rm{s}}.{\rm{t}}. & \;\; {{\boldsymbol{F}}_{i.}} \in {{{\varDelta }}_k}, i = 1, \cdots ,n; {\text{ }}{{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{F}} = {\boldsymbol{I}}. \end{aligned}\right\} $

式中： $ {{\boldsymbol{L}}_{{\text{sym}}}} = {{\boldsymbol{D}}^{ - 1/2}}{\boldsymbol{L}}{{\boldsymbol{D}}^{ - 1/2}} = {\boldsymbol{I}} - {{\boldsymbol{D}}^{ - 1/2}}{\boldsymbol{W}}{{\boldsymbol{D}}^{ - 1/2}} $为正则化Laplacian矩阵， ${\boldsymbol{W}}$矩阵由K近邻法构建.

由式（2）可知：

1）约束“ ${{\boldsymbol{F}}_{i.}} \in {{{\varDelta }}_k}, {\kern 1pt} {\kern 1pt} i = 1, \cdots ,n$”仅使得数据点对各个聚类中心的隶属度之和为1.0，而对各类样本点的隶属度总和没有限制. 在处理非平衡数据时，容易出现“大类”样本的隶属度总和过高，“小类”样本的隶属度总和过低的问题，可能会导致“均匀效应”问题的产生. 因此，须对大类及小类的隶属度总和施加“平衡”约束.

2）约束“ ${\kern 1pt} {{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{F}} = {\boldsymbol{I}}$”对隶属度矩阵施加正交约束“平衡”大类和小类的隶属度总和^[14].若约束隶属度矩阵 $ {\boldsymbol{F}} $满足近似正交（ $ {{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{F}} \approx {\boldsymbol{I}} $），设 $ {{\boldsymbol{F}}_{ \cdot p}} $表示 $ {\boldsymbol{F}} $的第 $ p $列，则存在

(3) $ {\Vert {{\boldsymbol{F}}}_{.p}\Vert }_{\text{F}}^{2}={\displaystyle \sum _{i=1}^{n}{{\boldsymbol{F}}}_{ip}^{2}}\approx 1; p=1,\cdots ,n. $

(4) $ {{\boldsymbol{F}}_{ \cdot p}} \cdot {{\boldsymbol{F}}_{ \cdot q}} = \sum\limits_{i = 1}^n {{F_{ip}}} {F_{iq}} \approx 0;p \ne q. $

由式（3）可知，每个类（簇）所含样本的隶属度平方总和近似相等，保证了各类间隶属度总平方和的平衡. 例如，设 $ {\boldsymbol{U}} = \left[ {{u_{ij}}} \right] $为原始隶属度矩阵，考虑第 $ p $和第 $ q $这2类（簇），设第 $ p $类（簇）所含样本远大于第 $ q $类（簇）所含样本. 若将 $ {\boldsymbol{U}} $近似正交化为 $ {\boldsymbol{F}} $，即 $ {\boldsymbol{U}} \to {\boldsymbol{F}} $且 $ {{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{F}} = {\boldsymbol{I}} $，则第 $ i $个样本点属于 $ p $类（簇）的隶属度为 $ {F_{ip}} = {{{u_{ip}}} \mathord{\left/ {\vphantom {{{u_{ip}}} {\sqrt {\sum\nolimits_j {{u_{jp}}} } }}} \right. } {\sqrt {\displaystyle \sum\nolimits_j {{u_{jp}}} } }} $，属于 $ q $类（簇）的隶属度变为 $ {F_{iq}} = {{{u_{iq}}} \mathord{\left/ {\vphantom {{{u_{iq}}} {\sqrt {\sum\nolimits_j {{u_{jq}}} } }}} \right. } {\sqrt {\displaystyle \sum\nolimits_j {{u_{jq}}} } }} $.由于 $ \sqrt {\displaystyle \sum\nolimits_j {{u_{jp}}} } > \sqrt {\displaystyle \sum\nolimits_j {{u_{jq}}} } $，故 $ {F_{ip}} < {F_{iq}} $，即小 $ q $类（簇）被赋予了更大的隶属度，从而平衡了大类、小类的隶属度. 另外，式（4）则约束不同簇的样本隶属度直接呈近似正交关系，使得不同簇的样本点隶属度之间关联性变小，排它性增强，提高其类间分离性.

利用文献[15]中的迹惩罚函数求解式（2），则问题（2）转化为如下罚函数模型：

(5) $ \mathop {\min }\limits_{{{\boldsymbol{F}}_{i.}} \in {{{\varDelta }}_k},i = 1, \cdots, n} {\kern 1pt} {\text{ }}{\kern 1pt}{\kern 1pt}g({\boldsymbol{F}}) = \frac{1}{2}{\text{tr}}\;\left( {{{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{LF}}} \right)+\frac{\mu }{4}\left\| {{{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{F}} - {\boldsymbol{I}}} \right\|_{\text{F}}^2. $

式中： $ \mu $为罚参数. 经计算可知，式（5）等价于

(6) $ \mathop {{\text{min}}}\limits_{{{\boldsymbol{F}}_{i.}} \in {{{\varDelta }}_k},i = 1, \cdots ,n} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\text{ }} {g} ({\boldsymbol{F}}) = \frac{\mu }{4}\left\| {{\boldsymbol{F}}{{\boldsymbol{F}}^{\text{T}}} - \left( {{\boldsymbol{I}} - \frac{{\boldsymbol{L}}}{\mu }} \right)} \right\|_{\text{F}}^2+\varepsilon {\text{.}} $

式中： $ \varepsilon $为常数项.

令 $ {{\boldsymbol{A}}} = {\boldsymbol{I}} - \dfrac{{\boldsymbol{L}}}{\mu } $， $ R\left( {\boldsymbol{F}} \right) = {\boldsymbol{F}}{{\boldsymbol{F}}^{\text{T}}} - {\boldsymbol{A}} $，则式（6）可以改写为

(7) $ \underset{{{\boldsymbol{F}}}_{i.}\in {\varDelta }_{k},i=1,\cdots ,n}{\mathrm{min}}\text{ }\text{ }g({\boldsymbol{F}})=\frac{1}{2}{\Vert R\left({\boldsymbol{F}}\right)\Vert }_{\text{{{F}}}}^{2}. $

由于式（7）的目标函数是二次函数，考虑应用Gauss-Newton（GN）法^[16]对其求解. 记 $ S({{\boldsymbol{F}}^t}) $为 $ {g} ({\boldsymbol{F}}) $在点 $ {{\boldsymbol{F}}^t} $处的GN方向，则GN的迭代格式为

(8) $ {{\boldsymbol{F}}^{t+1}} = {{\boldsymbol{F}}^t}+{\alpha ^t}S({{\boldsymbol{F}}^t}). $

式中： ${\alpha ^t} $为步长.

根据文献[15]，应用GN方向的解析表达式求解式（7），其中，GN方向的解析表达式为

(9) $ S({\boldsymbol{F}}) = \left[{\boldsymbol{I}} - \frac{1}{2}{\boldsymbol{F}}{({{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{F}})^{ - 1}}{{\boldsymbol{F}}^{\text{T}}}\right][{\boldsymbol{AF}}{({{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{F}})^{ - 1}} - {\boldsymbol{F}}]. $

平衡化谱聚类算法的求解步骤如算法1所示.

算法1 　BSPC算法

输入：正则化Laplacian矩阵 $ {{\boldsymbol{L}}_{{\text{sym}}}} $、期望的聚类数 $ k $、罚参数 $ \mu $、最大迭代次数 $ T_{\rm{iter}} $、终止精度 $ \varepsilon $；

输出：聚类隶属度矩阵 $ {{\boldsymbol{F}}^*} $.

1) 设置 $ t = 0 $，初始点 $ {{\boldsymbol{F}}^0} $；

2) 根据式（9）计算GN方向 $ \begin{gathered} {{\boldsymbol{S}}^t} = S({{\boldsymbol{F}}^t}) \\ \quad \\ \end{gathered} $；

3) 利用线搜索求 $ {{\boldsymbol{S}}^t} $方向的歩长 $ {\alpha ^t} $，根据式（8）更新 $ {\boldsymbol{F}} $为 $ {{\boldsymbol{F}}^{t+1}} = {{\boldsymbol{F}}^t}+{\alpha ^t}S({{\boldsymbol{F}}^t}) $；

4) 将 $ {{\boldsymbol{F}}^{t+1}} $投影到单纯形，即 ${\boldsymbol{F}}_{i.}^{t+1} \in {{{\varDelta }}_k},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i = 1, \cdots ,n{\kern 1pt} $；

5) 若 $\frac{{\left| {{g} \left( {{{\boldsymbol{F}}^t}} \right) - {g} \left( {{{\boldsymbol{F}}^{t+1}}} \right)} \right|}}{{\max \;\{ {g} \left( {{{\boldsymbol{F}}^{t+1}}} \right),1\} }} > \varepsilon$或 $ t+1 < T_{\rm{iter}} $，则令 $ t = t+1 $重复步骤 2）；

输出： ${{\boldsymbol{F}}^*} = {{\boldsymbol{F}}^{t+1}}$.

为了准确反映拱坝结构真实性态，考虑震级、频谱特性及持时的影响，从美国太平洋地震研究中心地震动数据库中，根据水工抗震规范反应谱共选择了109条地震波（每条地震波有3个方向），如图1所示. 图中，β为谱加速度放大系数，T_time为时间. 具体筛选准则如下：1）参照文献[17]的选波标准，以震级5~8级，震中距0~100 km，剪切波速 $ {V_{\rm{s}}} $>500 m/s作为基本筛选条件；2）每条地震波以0.1g为间隔，调幅10次至1.0g；3）地震波归一化处理，在非线性分析时考虑3个方向的地震动作用.

目前，研究人员已提出了大量的IM用于实现结构地震易损性研究^[18].本研究参考上述文献共选择了50个具有代表性的IM，详细定义如表1所示.

考虑《建筑设计抗震规范》^[8]以峰值加速度（peak ground acceleration，PGA）为设防标准，本研究以拱坝100 a基准期超越概率2%的工程场地基岩地震峰值加速度0.325g作为标准化指标. 在聚类时设置最大迭代次数 $ T{\text{ = }}1 \times {10^5} $，精度 $ \varepsilon {\text{ = }}1 \times {10^{{{ - }}6}} $. 采用肘部判别法^[19]确定109条地震波的聚类类别数，如图2(a)所示. 图中，k为聚类数；E_s为反映聚类算法畸变程度的簇内误差平方和，该值越小，聚类效果越显著. 可知，随着k逐渐增大，E_s逐渐降低，当k=7时，E_s接近最低水平，斜率近似为0，表明选择聚类数k=7，聚类效果较好. 如图2(b)所示为聚类算法收敛性效果图. 图中，I_t为迭代次数，g(F)为目标函数值. 可以看出，该算法收敛速度较快，在迭代范围内能够快速达到稳定效果.

理论上，计算的地震波数量越多，性能评估结果越可靠. Luco等^[20]研究发现，选取10~20条地震波可以精确评估结构的抗震能力. 因此，本研究将109条地震波调幅10次后的有限元结果（1090次）作为参照基准（整体样本）. 根据BSPC算法的思想，将109条地震波聚类后（共7簇）按照各簇隶属度值降序排列，在每簇样本中均匀选取15条地震波作为聚类样本. 其中，均匀选取规则由各簇样本数与整体样本数（109条）的比例决定，如

(10) $ {N_{\rm{C}}} = \frac{{15}}{{109}} \times N. $

式中： ${N_{\rm{C}}}$为每簇所选聚类样本点数，N为每簇地震波样本点数.

另一个对比案例仍以1.3节筛选标准作为基本条件，参考《建筑设计抗震规范》^[8]中多组地震动时程分析的标准规定，并结合学者在大坝^[21]结构地震易损性分析时应用的规范反应谱法，同样在109条地震波中选择15条作为反应谱样本. 考虑到本研究所选的109条地震波反应谱的离散性较大，很可能导致反应谱样本计算的反应谱均值与标准规范反应谱出现较大误差，因此，参考文献[22]的选波方式，以反应谱在[0.1,T_g]平台段以及结构前5阶周期附近段为筛选标准，要求所选的15条地震波的反应谱均值在上述2段范围内与标准规范反应谱的误差不超过10%.

反应谱样本的15条地震波如图1所示，该15条地震波的反应谱均值，尤其是对拱坝抗震性能起关键作用的谱周期1.0 s以下对应的均值，与本研究标准规范反应谱最接近（15条地震波的均值反应谱与标准反应谱在上述两频段内的数值误差超过10%的少于5个点）.

如图3所示为各类样本的空间分布情况. 图中，F₁、F₂、F₃为隶属度值. 在图3的整体样本图中，聚类后的7类样本数分别为10、12、13、13、15、21、25，可以发现样本数量在簇间均匀分布（每类样本数皆集中在约15条），表明该算法避免了传统谱聚类在筛选样本时产生的均匀效应问题，提升了在类不平衡地震波数据集上的聚类效果，更符合工程实际.

与以往聚类算法不同的是，平衡化谱聚类算法通过转换正则化Laplacian矩阵达到降维的效果. 在将原始数据转换为Laplacian矩阵时，对 $ n $维Laplacian矩阵计算 $ k $个最小特征值，此时聚类个数 $ k $远远小于维数 $ n $，即降维到 $ k $，从而解决了地震动强度指标过多导致的“维数灾难”问题.

如图3所示的分布图表明，相较于聚类样本，反应谱样本存在过于集中的现象；聚类样本能够体现整体样本的各簇特征，覆盖整体样本的空间更广. 主要原因如下：规范反应谱法对频率特性的限制，导致筛选出的反应谱样本的类型过于单一. 考虑到地震易损性分析时所需地震波的随机性，如果所选地震波未能完整覆盖整体样本的关键特征，易损性结果很可能无法全面体现结构的震害规律.

多样条带分析法是根据最大似然估计原理得到结构的地震易损性概率. 靳聪聪等^[23]的研究表明，MSA法在一定结构分析数量下比IDA法能够更有效地估计易损性方程参数. 易损性函数描述结构在对应地震动强度下达到或超过规定损伤极限状态的概率，表达式为

(11) $ {p} = \varPhi \left(\frac{{\ln\; (X/\theta )}}{{{\sigma _{\rm{M}}}}}\right). $

式中： $ {p} $为结构在 $ {\text{IM}} = {X} $时发生破坏的概率， $\varPhi $为标准正态分布函数， $\theta$和 ${\sigma _{\rm{M}}}$分别为达到性能水准所需地震强度因子的平均值和对数标准差.

假定结构在不同地震动激励下，当 $ {\text{IM}} = {x_i} $时，存在 $ n_{i} $条地震波中有 $ z_{i} $条地震波导致结构破坏的概率由二项分布给出，表达式为

(12) $ P\left(z_{i}\right)=\left(\begin{array}{l} n_{i} \\ z_{i} \end{array}\right) p_{i}^{2}\left(1-p_{i}\right)^{n_{i}-z_{i}} . $

假设每条地震动激励下结构发生破坏的概率是独立的，则似然函数如下：

(13) $ L\left(\begin{array}{c} \theta \\ \sigma_{{\rm{M}}} \end{array}\right)=\prod_{i=1}^{m} P\left(z_{i}\right)=\prod_{i=1}^{m}\left(\begin{array}{l} n_{i} \\ z_{i} \end{array}\right) p_{i}^{z_{i}}\left(1-p_{i}\right)^{n_{i}-{{z}}_{i}}. $

将式（11）代入式（13），得到

(14) $\begin{split} L\left(\begin{array}{c}\theta \\ \sigma_{\rm{M}}\end{array}\right)=&\prod_{i=1}^m\left(\begin{array}{c}n_i \\ z_i\end{array}\right) \varPhi\left(\frac{\ln\; \left(x_i / \theta\right)}{\sigma_{\rm{M}}}\right)^{z_i}\times \\ &\left(1-\varPhi\left(\frac{\ln \;\left(x_i / \theta\right)}{\sigma_{\rm{M}}}\right)\right)^{n_i-z_i}. \end{split}$

式中：m为IM水平数量. 易损性方程参数 $ \theta $和 $ {\sigma _{\rm{M}}} $的最大似然估计如下：

(15) $ \begin{split} \left\{\hat{\theta }\text{，}{\hat{\sigma }}_{{\rm{M}}}\right\}=&\mathrm{arg}\mathrm{max}\;\;{\displaystyle \sum _{i=1}^{m}\left\{{\rm{ln}}\;\left(\begin{array}{c}{n}_{i}\\ {z}_{i}\end{array}\right)+{z}_{i}{\rm{ln}}\;\varPhi \left(\frac{{\rm{ln}}\;({x}_{i}/\theta )}{{\sigma }_{{\rm{M}}}}\right)+\right.}\\ &\left.({n}_{i}-{z}_{i}){\rm{ln}}\;\left[1-\varPhi \;\left(\frac{{\rm{ln}}\;({x}_{i}/\theta )}{{\sigma }_{{\rm{M}}}}\right)\right]\right\}. \\[-18pt] \end{split}$

以双曲型白鹤滩拱坝为例构建模型. 该拱坝工程场地类别为Ⅰ类，场地基本地震烈度为Ⅷ度，100 a基准期超越概率2%的工程场地基岩地震PGA为 0.325g^[17].拱坝坝顶高程834 m，坝高289 m，坝体上游正常蓄水位825 m，下游水位604 m，更多模型参数可参考文献[24]、[25]，如表2所示. 图中，E为弹性模量，ρ为密度，ν为泊松比. 考虑到工程实际，对现有数值模型进行适当简化，地基部分延伸为约2倍的拱坝坝体结构高度及宽度. 该有限元模型共36178个单元、40286个节点，坐标y轴负向为顺水流向，z轴正向向上. 坝体及地基均为Solid 65单元，少部分坝肩衔接处采用六节点棱柱体单元.

在进行非线性分析时在地基边界截断处设置弹簧-阻尼器实现黏弹性边界的施加^[26].坝体共计30条横缝，考虑到计算时拱坝收缩缝产生的开合-接触-滑移现象，选用商业软件ABAQUS提供的动态接触模型. 在接触设置时仅考虑横缝的张开与闭合，接触面之间的相互作用遵循接触-嵌入原则^[27].同时，引入Lee等^[28]提出的混凝土塑性损伤本构研究坝体损伤特性，本构设置参数可以参考文献[29]. 此外，由于本研究拱坝模型仅在坝面上下游铺设钢筋，坝面位置仍采用混凝土的弹性模量代替钢筋混凝土的弹性模量，并忽略钢筋的强化阶段，选择理想本构模型代替^[30]. 针对白鹤滩拱坝进行模态计算，选择Chen^[31]修正的动水压力模型实现正常蓄水位的模拟，获取的前5阶模态振型图如图4所示. 该模态结果与文献[32]中白鹤滩拱坝提取正常蓄水位下的振动数据响应模态较吻合，最大模态差异值仅6%.参考关于拱坝抗震的类似文献成果^[26-27]，在地震分析时选用瑞利阻尼，阻尼比为5%.

目前，学者们将峰值位移、应力分布、损伤等性能指标作为拱坝性能分析时的常用指标. 姚霄雯^[33]根据拱坝概率需求曲线的斜率突变特性，分别选择拱冠梁顶部位移、横缝开度及切向滑移指标评估拱坝抗震性能. Wang等^[34]结合大岗山拱坝的损伤分布及横缝开度，量化了性能指标的极限状态. 李静等^[35]认为坝面损伤面积比及坝体损伤体积比可以清晰地反映拱坝结构损伤规律. 范书立等^[17]建立以拱冠梁顶部位移、横缝开度、损伤体积比3项性能指标为基准的白鹤滩拱坝易损性曲面，Xu等^[36]证实上述3项指标具备评估拱坝抗震性能的能力. Song等^[37]针对白鹤滩拱坝开展了稳定性实验研究，表明拱坝坝顶位移值得关注.

在现行的各类结构抗震规范中，位移仍是常用的衡量标准. 从宏观角度来说，依照位移、损伤定义的性能指标也更容易直观地描述结构极限承载能力. 基于此，参考上述研究成果，分别以坝顶上游方向位移（DT）、横缝最大开度（C）、损伤体积比（DVR）作为拱坝结构性能指标，其中损伤体积比定义为

(16) $ {\text{DVR}} = \frac{{\displaystyle \sum {_e\int {_{{v_e}}{{{d}}_{\text{t}}}{\text{d}}{v_e}} } }}{{\displaystyle \sum {_e\int {_{{v_e}}{\text{d}}{v_e}} } }}. $

式中：d_t为宏观开裂损伤指数， $ {v_e} $为每个单元的体积，规定 $ {{{d}}_{\text{t}}}{\text{ > }}0.8 $.

结合文献[17]对白鹤滩拱坝性能指标的定义标准，将拱坝损伤发展历程的4个阶段作为极限状态划分的典型标志，分为基本完好、轻微破坏、中等破坏、严重破坏4类，并以每条地震波达到上述4种极限状态时的均值量化结构性能指标限值. 拱坝各类样本的性能水平划分取值如表3所示.

根据表3的极限状态取值，将1090次计算数据（整体样本）、15条聚类及反应谱样本的计算数据划分为4种极限区间，绘制如图5所示的统计图. 图中，L_s表示各性能指标的极限状态值.可以看出，在各极限状态下，整体样本与反应谱样本、聚类样本的各L_s(DVR)、L_s(DT)均存在相应误差；聚类样本的L_s( $ C $)与整体样本的基本吻合，表明在此类性能指标下，聚类样本的计算数据对拱坝地震响应L_s的影响与整体样本的具备一致性.

对比整体样本，反应谱样本的计算数据分布更紧凑，部分区域已出现空白. 其中，在SD极限状态下，L_s(DVR)的计算数据仅覆盖到0.4g，L_s(DT)与L_s( $ C $)指标值也仅覆盖到0.5g，未能体现整体样本在其他地震动强度下的计算数据点. 观察MD、ED区间，也同时存在着计算数据点分布不足的现象. 出现这种现象的主要原因很可能是应用反应谱法筛选的地震波未能完全体现整体样本的关键特征，从而会丢失一些能够反映整体样本的关键计算点. 相对来说，聚类样本的计算响应值在不同地震动强度下分布更稀疏，覆盖整体样本的效果更明显. 除了L_s(DVR)指标的少数计算数据存在缺失，L_s(DT)、L_s( $ C $)指标基本占据了整体样本计算数据的各分布区域，对各极限状态的区间皆有良好表现，这与聚类算法在筛选地震波时体现整体样本各簇特征的基本思想保持一致，表明聚类算法效果显著.

由于3类样本的有限元结果在不同地震动激励下的标准差及均值皆不相同，因此，为了评估3类样本在不同地震动强度下有限元结果的接近程度，选用变异系数CV进行分析：

(17) $ {\text{CV}} = \sigma /v .$

式中：σ、v分别表示不同地震波在相同地震动强度激励下响应值的标准差与均值.

各类样本的CV分布情况如图6所示. 可以看出，3类样本的CV(DVR)均分布在0.20~1.80，CV(DT)、CV( $ C $)主要集中在约0.75.反应谱样本、聚类样本与整体样本在CV(DVR)上的最大差异分别为0.148、0.142，在CV(DT)、CV( $ C $)上的最大差异分别为0.085、0.128与0.063、0.095，表明聚类样本与整体样本的误差更小. 此外，在所有地震动强度下，3项指标的聚类样本与整体样本的所有CV误差大多小于0.050，表明聚类算法筛选的地震波与整体样本的计算结果存在较高的一致离散性.

如图7所示为各类样本的地震易损性曲线. 图中，P_f为易损性概率. 可以看出，所有易损性曲线在地震动强度变化范围内皆呈现出斜率先增大后平稳的趋势；聚类及反应谱样本与整体样本的易损性结果基本保持一致，所示曲线的最终P_f均大于80%.

在SD极限状态下，当PGA=0.2g时，整体样本的P_f(DT)与P_f( $ C $)分别为56.46%、58.29%，同等条件下，反应谱样本的P_f(DT)与P_f( $ C $)分别为66.18%、60.90%，聚类样本的P_f(DT)与P_f(C)分别为53.58%、51.65%.当PGA=0.4g时，所有类别的P_f皆大于80%.可以看出，本研究拱坝结构在发生超过0.4g的地震时，坝踵弱部位产生小规模损伤及裂缝的可能性较大.

在MD极限状态下，当PGA=0.3g时，整体样本的P_f(DT)与P_f( $ C $)分别为55.70%、56.22%，此时，所有类别样本的易损性曲线上升趋势皆逐渐变缓，斜率降低. PGA＞0.6g时，所有易损性曲线的P_f皆逐渐超过80%，表明当PGA＞0.6g时，拱坝薄弱部位结构很可能已经发生宏观开裂现象. 例如，拱坝下游面中部及上部结构，均会产生不同程度的损伤，须进行少量修复.

在ED极限状态下，当PGA=0.1g时，所有类别样本的各P_f均保持在1%~7%，当PGA=0.5g时，3类样本的P_f(DT)分别为59.57%、63.22%、57.02%.当PGA=0.8g时，3类样本的P_f逐渐超过80%，表明发生大于0.8g的地震时，个别部位如坝踵及坝体下游中上部很可能已经出现严重的开裂现象，影响拱坝的安全稳定.

综上所述，3类样本在PGA＜0.4g时，所有SD极限状态的易损性曲线皆呈现出斜率较大的陡峭趋势，说明当发生类似强度的地震时，本研究拱坝将很难维持完好无损状态，容易超过性能指标的第一极限状态. 并且，所有易损性曲线的最终P_f皆能达到80%以上，表明本研究拱坝仍受到强震威胁.

从概率的角度审视聚类及反应谱样本的易损性曲线差异. 以整体样本的易损性曲线作为参照基准，分别将PGA=0.1g~1.0g的反应谱样本、聚类样本的P_f进行误差统计，分析不同地震动强度激励下两者的误差分布情况. 如图8所示为各类样本地震易损性曲线误差分析图. 图中，E_r为易损性概率误差所占百分比，P_D为误差概率密度值.

图8中左图表明，不同地震动强度下反应谱样本与整体样本的易损性概率差异点分布稀疏，E_r主要分布在0~15%，绝大多数E_r为5%~10%.相比而言，聚类样本与整体样本的易损性曲线更吻合，除极个别结果外，两者E_r主要集中在1%~5%.并且，在所有地震动强度下，聚类样本与整体样本的最大E_r仅为4.39%、3.84%、6.64%，明显小于反应谱样本的E_r（13.37%、11.39%、7.83%）.

图8右图为误差正态分布概率密度函数曲线. 可以看出，聚类样本与整体样本的E_r(DVR)，E_r(DT)小于5%的最小概率分别为92.24%、99.19%. 同等条件下，反应谱样本的计算结果仅为30.09%、82.14%. 在SD、MD极限状态下，聚类样本的E_r(DVR)、E_r(DT)不超过5%的概率皆接近100%，而此时反应谱样本的E_r(DVR)不超过5%的概率仅为49.72%、74.85%，E_r(DT)不超过5%的概率仅为38.88%、53.94%. 只有在SD极限状态下的E_r( $ C $)指标上，聚类样本该值不超过5%的概率小于90%，但仍远远大于该条件下的反应谱样本概率（53.63%）.

此外，聚类样本在SD极限状态下，与整体样本的最大E_r( $ C $)仅为6.64%，仍然小于同等条件下的反应谱样本E_r( $ C $)（7.83%）.并且，在另外2种极限状态MD、ED下，聚类样本误差E_r( $ C $)<5%的概率皆接近100%，也远远大于该状态下反应谱样本的概率（53.63%、74.16%）.

另外，较聚类样本而言，除了E_r( $ C $)的SD极限状态，反应谱样本的所有误差正态分布的均值及标准差值皆更大，证实聚类样本与整体样本的地震易损性曲线更接近，同时也表明采用聚类样本代替整体样本进行拱坝地震易损性评估可靠性更高. 仅通过少量的样本子集可以达到100多条整体样本的预测效果，验证了本研究聚类算法具备描述整体样本典型特征的能力，在保证地震易损性概率准确性的同时可以显著降低计算成本. 这一结论也能够根据图3、5得出，图5中就各类样本的地震响应值的分布情况预测了地震易损性结果，表明在地震易损性分析时，充分考虑地震动的随机性不容忽视.

值得说明的是，应用规范反应谱法可以获得多数满足要求的地震波，但该方法也容易忽略与结构或场地相关的一些具有代表性的地震波（如图3的边缘部分），无法充分考虑地震动的随机性，使得易损性的结果很可能不够准确. 在结构抗震性能评估时，聚类算法是对传统地震波筛选方式的一种补充. 该算法根据地震动强度指标的基本特性将原始地震波划分为不同簇，再从每簇样本中筛选出具备代表性的地震波，弥补了传统选波方式的不足.

（1）所提BSPC算法筛选出的聚类样本具备描述整体样本典型特征的能力，在损伤体比、坝顶位移、横缝最大开度这3个性能指标中，两者最大误差仅为4.39%、3.84%、6.64%，两者误差不超过5%的最小概率分别为92.24%、99.19%、81.75%，表明聚类算法筛选的样本与整体样本的易损性结果具备一致性.

（2）在SD极限状态时，3类样本的地震易损性曲线在PGA＜0.4g的条件下，皆呈现出斜率较大的陡峭趋势，说明当发生类似强度的地震时，本研究拱坝结构将很难维持完好无损状态，容易超过结构性能指标的第一极限状态.

（3）在所有极限状态下，聚类样本、反应谱样本与整体样本的地震易损性结果基本保持一致，最终的易损性概率均大于80%，表明本研究拱坝仍然容易受到强震威胁.

（4）本研究提出的聚类算法BSPC是对传统筛选地震波方式的一种补充. 该算法根据地震动强度指标的基本特性从原始地震波中筛选出了具备代表性的地震波子集，弥补了传统选波方式的不足.

（5）本研究侧重于谱聚类算法的改进及在地震工程领域中降低易损性计算量的应用，但本研究仅在MSA易损性方法下开展. 因此，后续应当研究聚类算法在不同易损性方法上的应用，充分考虑聚类算法对各易损性方法的潜在影响.