0.   引言

近年来，海洋开发和利用不断加深，特别是中国海上风电、跨海大桥、海上石油平台、海底隧道、海底管道等重大海洋工程投入使用。与此同时，中国海域及邻区所处地震地质环境复杂，地震活动性较强，海域地震频发^[1]，海洋工程的抗震设防能力及重要性凸显，海域场地工程特性得到充分关注^[2]。针对中国海域及邻近区域的场地地震反应研究工作仍不足^[1]，对海域土层条件，场地特性及分类等研究还处在起步阶段^[3]，所以研究海域场地效应和海域场地分类可以为海洋工程抗震设防，海洋结构安全性评价提供帮助。

Nakamura^[4]提出基于地脉动的水平和竖向谱比（horizontal-to-vertical spectral ratio, HVSR）作为场地特征评估方法，国内外涌现很多基于HVSR的场地地震效应和场地分类的研究。其中Boore等^[5]认为水平向地震动与竖向地震动之比（H/V）可以消除震源、断层类型、传播路径等因素，用以分析场地效应；同时Boore等^[5]认为H/V谱比法可同时应用于陆域和海域场地的研究。Field等^[6]利用HVSR法研究了1994年北岭地震场地非线性放大效应。Bonilla等^[7]利用南加州地区的井下阵列钻孔记录对比了HVSR法和传统的场地特征研究方法，发现两种方法得到的场地卓越周期基本一致。Wen等^[8]利用HVSR法分析了集集地震的场地非线性反应。任叶飞等^[9]将广义线性反演法和HVSR法进行对比，得出HVSR方法可以应用于评估场地卓越周期。荣棉水等^[10]利用美国Garner Valley Downhole Array（GVDA）竖向观测台阵的记录，对比了传递函数法和HVSR法的差异，并指出在可忽略场地竖向放大的情况下，可以利用HVSR研究场地特性。Ren等^[11]分析了2008年汶川地震的场地非线性效应，并提出了5个关键参数来表征场地非线性。Dhakal等^[12]利用K-NET六座海底台站记录到的315次地震事件分析了海底台站的场地效应，当PGA > 50gal时，海底部分台站所在场地出现非线性。Kubo等^[13]基于日本DONET1观测网记录到的2016年日本Off-Mie地震，利用HVSR法分析了台站所在场地的场地效应，认为部分台站场地的非线性反应抑制了地震动的放大效应。Ji等^[14]将HVSR计算出的非线性响应程度（DNL）和非线性响应百分比（PNL）联合PGA和V_S30（或V_S20）作为变量，采用k-means方法识别非线性场地响应。Yaghmaei-Sabegh等^[15]利用2017年伊朗Ezgeleh地震所记录到的地震数据，提出了基于广义回归神经网络（generalized regression neural networks，GRNN）的HVSR分类方法，该方法仅利用地震记录就可进行场地类别的划分，并将其与Vs30场地分类方法进行对比，展现了GRNN方法的可靠性。周旭彤等^[16]利用HVSR方法对日本DONET1台网20座海底台站场地效应进行分析，研究结果表明，部分台站具有相似的场地特征，海底台站的布设方式和海底地形会对谱比结果产生影响。

本文基于日本K-NET台网的海底台站收集的地震数据，利用HVSR联合动态聚类的方法对海底台站的场地特征进行识别，给出了台站所在场地的分类结果，为研究海域场地效应提供参考。

1.   海域强地震动观测数据与谱比方法

1.1   海域强地震动观测现状及数据处理

海域强地震动观测与陆域相比，存在设备安置难度大、回收效率低及观测稳定性较差等技术难题，发展远落后于现有陆域观测水平及海洋工程建设需求。美国和日本较早开展海域地震动观测，中国台湾地区也布设了海底地震观测设备^[17]，美国在20世纪70年代（1977年）针对海底地震记录观测需求，提出了SEMS计划（seafloor earthquake monitor system，以下简称SEMS）^[3]，日本于90年代（1996年）在东京圈近海开展ETMC计划（earthquake and tsunami monitoring cable，以下简称ETMC）^[18]。美国SEMS计划开发了4代海底地震仪，安装在南加州海域，水深区间为数10~300 m。日本ETMC系统布置在相模湾，观测系统由海底强震仪及海啸压力传感器组构，水深区间在900~2000 m不等。本文收集、整理了日本K-net台网6个海底台站2000—2018年的海域地震动数据。海底台站信息见表 1，分布见图 1。

表  1  海底台站信息

Table  1.  Information of seabed stations

台站编号	东经/(°)	北纬/(°)	海水深度/m
KNG201	139.9183	34.5956	2197
KNG202	139.8393	34.7396	2339
KNG203	139.6435	34.7983	902
KNG204	139.5711	34.8931	933
KNG205	139.4213	34.9413	1486
KNG206	139.3778	34.0966	1130

下载: 导出CSV 

| 显示表格

图  1  海底台站与水深分布图

Figure  1.  Distribution map of submarine stations and water depth

下载: 全尺寸图片 幻灯片

为研究6个台站所在海域的场地效应，考虑噪声对地震记录的影响，本文对上述数据进行筛选，筛选原则如下：①震级大于4级且PGA < 20 cm/s²；②震源位置均在海域。据此选取338组地震记录，每组记录均包含三分量，震级震中距分布见图 2。由图 2可知，震级主要分布在4.0~5.5级，震中距主要分布在0~300 km。

图  2  震级震中距分布图

Figure  2.  Distribution map of magnitude and distance

下载: 全尺寸图片 幻灯片

针对原始地震记录可能出现的基线漂移、仪器噪声及环境噪声等信号偏差，对筛选后的海域地震动记录采用基线校正和滤波的方法进行处理，其中采用Boore提出的分段式基线对地震动记录进行校正^[19]，构建4阶Butterworth带通滤波器，截止频率为0.1~20 Hz，处理前后的地震记录对比见图 3，由图可知，未处理地震动将出现基线漂移，中、长周期信号失真。处理后的地震记录更符合实际物理过程。

图  3  地震记录处理

Figure  3.  Processing of seismic records

下载: 全尺寸图片 幻灯片

1.2   海域强地震动记录HVSR计算

在利用强震记录研究场地效应时，需尽量选取S波时段^[20]。姚鑫鑫等^[21]关于S波的截取对HVSR谱比结果影响的研究表明，截取S波可以更好地体现高频段场地放大效应。因此本文对筛选后的海域地震动记录截取S波到时，将Arias强度能量^[22]达到90%所对应的持时作为S波的结束时间。

对截取后的S波数据进行快速傅里叶变换（FFT），选取0.2~20 Hz频带范围，并利用Konno平滑方法^[23]对傅里叶振幅谱进行平滑，平滑系数b取20。图 4（a）表明，利用该平滑方法既不低估场地放大效应的幅值，又不影响卓越频率的识别，可以较好地体现各频段的信息。将平滑后的地震记录S波的EW、NS向利用几何平均合成水平向，并除以相应的UD向平滑后的傅里叶振幅谱，得到的HVSR结果见图 4（b）。

图  4  HVSR计算过程

Figure  4.  Calculation process of HVSR

下载: 全尺寸图片 幻灯片

2.   动态聚类学习原理与评价

2.1   动态k-means聚类学习原理

Zhang等^[24]完善了动态$k$-means聚类方法，针对反应谱数据开展了聚类分析，本文通过添加谱比序列标签，采用动态聚类算法实现HVSR谱比分类。动态聚类学习原理如下：针对$m$个数据序列$\boldsymbol{Z} =$ $\left\{\boldsymbol{z}_1, \boldsymbol{z}_2, \cdots, \boldsymbol{z}_m\right\}^{\mathrm{T}} \in \mathbb{R}^{m \times n}$，各序列由$n$个样本点组成行向量，即${{\bf{z}}_i} = \{ {z_{i1}}, {z_{i2}}, \cdot \cdot \cdot , {z_{in}}\}$，其中$i = 1, 2, \cdot \cdot \cdot , m$。通过聚类、迭代方法，将数据序列实现分类，即将计算序列${\bf{Z}}$和各类数据中心$\boldsymbol{C}=\left\{\boldsymbol{c}_1, \boldsymbol{c}_2, \cdots, \boldsymbol{c}_k\right\}^{\mathrm{T}} \in \mathbb{R}^{k \times n}$的欧式距离聚类：

式中：$\boldsymbol{c}_p=\left\{c_{p 1}, c_{p 2}, \cdots, c_{p n}\right\}$是$p$个类中心，由n个样本点组成的行向量。$\boldsymbol{U} \in {\mathbb{Z}^{z \times k}}$为二元矩阵，当${\lambda _{ip}} = 1$表示第$i$个序列分配给$p$类，并且每个序列仅可分配给同类，即$\sum\nolimits_{p = 1}^k {{\lambda _{ip}} = 1}$。通过设置权重矩阵$\boldsymbol{W}$，可更好识别局部数据跨度内时间序列的相关性，提高聚类性能，每个序列的权重之和为1，即$\sum\nolimits_{j = 1}^n {{w_{pj}} = 1}$。为了平滑在相邻序列间的权重，防止陷入局部最小，设置约束正则项，定义如下：

式中，$\alpha$为平滑参数，$\alpha = \sum\nolimits_{i = 1}^m {\sum\nolimits_{j = 1}^n {({z_{ij}} - \sum\nolimits_{i = 1}^m {{z_{ij}}/} } }$ $m{)^2}$。

综上，动态聚类可以被定义为具有正则化目标函数的优化问题：

式中，${m_0}$为每类所存在的最小序列数。$\sum\nolimits_{i = 1}^m {{\lambda _{ip}} \geqslant {m_0}}$代表$p$类中的每一类的序列数至少是m₀。式（3）的第一项使得所有类内的加权总距离最小，第二项使相邻类间隔的权重更平滑，即类间的区别更明显。通过设置计算序列到各中心的欧式距离，当此距离越近表示序列越相似。

2.2   动态聚类学习流程

该方法旨在最小化目标函数式（3），式（3）的变量由二元矩阵$\boldsymbol{U}$，权重矩阵$\boldsymbol{W}$和聚类中心矩阵$\boldsymbol{C}$构成。采用优化方法，即将其中两个变量作为已知量，采用信息更新方法优化第三变量，当迭代次数达到40次，目标函数（3）收敛条件为10^-6时停止更新。整个优化过程可以分为如下3步。

（1）给定假设的权重矩阵$\boldsymbol{W}$和聚类中心$\boldsymbol{C}$，当式（1）满足式（4），序列被分配到假设聚类中心的类中，各类中心与序列之间的权重距离被最小化。

（2）当确定权重矩阵$\boldsymbol{W}$和二元矩阵$\boldsymbol{U}$后，通过令式（1）的导数为0求解信息更新后的新聚类中心，具体过程如下：

变换得

通过式（6）可得到更新后的聚类中心。

（3）通过第一、二步确定二元矩阵$\boldsymbol{U}$和聚类中心$\boldsymbol{C}$，此时当式（3）满足

式（3）转化为计算权重矩阵$\boldsymbol{W}$的二次规划问题，具体求解方法可以参考文献[25]。

2.3   动态k-means聚类方法验证

为验证该方法在HVSR分类中的性能，选取5%阻尼比的反应谱来进行方法验证。依据《建筑抗震设计规范：GB50011—2010》^[26]，在场地分类为III类条件下，生成两组设计反应谱，两类反应谱见图 5，图 5中第一类反应谱的特征周期为0.9 s，地震影响系数最大值为0.9；第二类反应谱特征周期为0.3 s，地震影响系数最大值为0.79。

图  5  设计反应谱

Figure  5.  Design response spectra

下载: 全尺寸图片 幻灯片

依据上述不同特性反应谱，从PEER强震动数据库中^[27]，分别针对各分类选取50条地震记录。图 6给出了数据的选取情况。

图  6  不同记录反应谱特性

Figure  6.  Two clusters of response spectra

下载: 全尺寸图片 幻灯片

将两组数据利用动态聚类学习方法进行分类测试，结果见图 7。由图 7可知，采用动态聚类学习方法划分的中心和反应谱均值无偏差，展现出良好的聚类性能及鲁棒性。

图  7  反应谱均值与中心匹配

Figure  7.  Mean values of response spectra and cluster centroid

下载: 全尺寸图片 幻灯片

3.   海域地震动谱比法分类

3.1   海域地震动谱比法分类依据

将HVSR结果利用上述动态聚类方法进行聚类。利用误差平方和准则评价聚类性能，表示聚类的畸变程度。给定数据序列$\boldsymbol{Z}$，假设$\boldsymbol{Z}$中的序列可分为$k$类，即${\boldsymbol{Z}_1}, {\boldsymbol{Z}_2}, \cdot \cdot \cdot , {\boldsymbol{Z}_k}$；各类序列数分别为${m_1}, {m_2}, \cdot \cdot \cdot , {m_k}$；各类中心分别为${c_1}, {c_2}, \cdot \cdot \cdot , {c_k}$，则误差平方和准则函数为

式中，$D$为数据序列分为$k$类时的每类数据序列与其中心序列的距离和。$D$越小表示每类数据序列到对应中心的距离越小，越紧凑，分类性能越好。分别令$k = 1, 2, \cdot \cdot \cdot , 8$进行聚类计算并计算$D$，计算结果见图 8。从图 8中可以看出，随着分类数增大，误差平方和趋于稳定，结合肘部法则^[28]确定最优分类数为$k = 3$。此外为保证每类结果具有统计意义，设置每类最少存在6组HVSR数据序列，即${m_0} = 6$。

图  8  误差平方和与分类数关系

Figure  8.  Identification of optimal number of cluster at the elbow point

下载: 全尺寸图片 幻灯片

3.2   海域场地特性分析

HVSR动态聚类结果见图 9，对应台站统计如图 10所示。由图 9，10可知：海域场地HVSR曲线无明显峰值或多峰值现象，结合肘部算法，相模湾海区最优化场地可分为三类，结合相模湾海区地貌^[29-30]，综合分为海底峡谷区（KNG203，KNG204）及海洋盆地区（KNG201，KNG202，KNG205，KNG206），两类地貌分类的HVSR谱比呈显著差异，海底峡谷区高频衰减慢且优势频段主要集中在中高频，台站KNG204代表的第一类主要集中在12 Hz，体现了海底峡谷中少量沉积的薄覆盖层环境；台站KNG203代表的第三类优势频带跨度较大（5~10 Hz），台站KNG201，KNG202，KNG205，KNG206则表现为中长周期发育，高频衰减快，呈现出典型的深厚覆盖层的海洋盆地特性。

图  9  HVSR聚类结果

Figure  9.  Clustering results of HVSR

下载: 全尺寸图片 幻灯片

图  10  各海域台站HVSR与类别对应

Figure  10.  HVSR of each station corresponding to category

下载: 全尺寸图片 幻灯片

其次，局部场地地形效应也显著影响海域地震动传播，以台站作为中心点，以类似地貌作为依据，通过采样点与标准点之间的距离（500 m）与高程计算局部场地地形坡度：

式中，h为采样点与标准点的平均高程差，l为采样点与标准点距离。

各台站所在局部场地地形坡度见表 2。由表 2可知：海底峡谷区台站KNG203较KNG204场地局部地形平缓度高，存在局部沉积形成的浅覆盖层；海洋盆地区台站KNG201，KNG205，KNG206所在场地局部地形坡度基本一致，KNG202所在场地处于盆地边缘，存在基岩面拉伸导致的局部场地不均匀分布特性。

表  2  局部场地地形坡度

Table  2.  Topographic slopes of local site

类别	台站	地貌	坡度	地形	场地特征
一	KNG204	海洋峡谷	4.62	陡峭	薄覆盖层
二	KNG201	海洋盆地	1.39	平缓	厚覆盖层
	KNG205		1.85
	KNG206		1.02
	KNG202		3.81	陡峭	薄覆盖层
三	KNG203	海洋峡谷	0.38	平缓	浅覆盖层

下载: 导出CSV 

| 显示表格

4.   结论

本文采用动态聚类学习方法结合海域场地HVSR谱比，分析了动态聚类方法在地震动分类的适用性；给出了相模湾海区场地工程特性，结合海洋工程地质给出了海域地形坡度与场地分类关系。希望通过对有限的海域场地地震记录研究，为近海以及海洋工程结构和设施的设计提供合理、科学的工程场地参考。得到以下2点结论。

（1）动态聚类学习方法针对时序数据的聚类性能及鲁棒性强，可满足局部及全局分辨率，应用于场地工程特性分析中，可结合HVSR谱比，最优化、自动化的给出场地分类特性；

（2）HVSR动态聚类的海洋场地特性表明：台站KNG203和KNG204位于海底峡谷区，台站KNG201、KNG202、KNG205、KNG206位于海洋盆地区，两类地貌分类对海域地震动传播影响较大，海底峡谷区高频衰减慢且优势频段主要集中在中高频，通过局部场地地形坡度分析给出了各台站所在场地的工程特性。