引言

采用随机有限断层法进行地震动模拟时，地震动传播过程中的衰减包括两部分，其一是路径衰减，即随距离$R$的增加而使振幅减小的几何衰减，以及由于地球介质对波能的吸收和耗散产生的非弹性衰减；其二是场地效应的影响，即需要考虑近地表地震动幅值影响和与路径无关的地震动高频衰减。Anderson和Hough （1984）引入kappa滤波器，即exp （−πfκ）（f为频率，κ为高频衰减系数），来模拟地震动高频部分的衰减，局部场地的地震动高频衰减则可采用0 km处的κ值κ₀表示。 κ表征在某一频率之后，地震动傅里叶幅值谱会随着频率的增大而减小。 如何估计特定场址的κ值，基于不同的研究目标已有很多计算方法，例如：Anderson和Humphrey （1991）提出的宽带反演方法，其关键在于假定震源谱形状，此种方法不受震级大小的约束；Kilb等（2012）提出固体应力降方法，适用于各种震级，主要用于拐角频率不确定、不易测定或探究拐角频率影响。工程地震领域场址地震动模拟中，普遍采用由Anderson和Hough （1984）提出的加速度反应谱法估计κ值，直接从加速度记录的傅里叶幅值谱的高频部分测量κ，在半对数坐标系下，高频衰减系数κ可以用最简单的线性近似得到。但这些估计κ值的方法普遍依赖于实际地震记录，在预测未来某场址的地震动时，很难直接从实际地震记录来获取研究区域的κ值，因此需要构建一个用于估计κ值的实用模型，使其能够在对不同类型场址的地震动预测中合理地模拟地震动高频衰减。 然而，κ值的影响因素很多，例如：Hanks （1982）认为高频衰减现象与场地性质有关；Boore （1983）认为κ值与震源或场地条件有关，而与传播路径无关；Anderson和Hough （1984）提出κ值与震中距呈正相关，且主要受近地表地质构造影响；Atkinson （1984）则认为κ值与场地和震源都存在一定的关系。

虽然不同学者对影响κ值的因素尚未达成一致的结论，但对局部场地地震动高频衰减系数κ₀影响因素的研究表明，κ₀与平均剪切波速 v_S30之间存在比较显著的相关性（Silva et al，1998；Chandler et al，2006；Drouet et al，2010；van Houtte et al，2011；Laurendeau et al，2013；Edwards et al，2015）。 比如，Anderson等（ 1996 ）模拟了S波在土层中的传播，结果表明κ₀对地震动峰值振幅的影响与平均剪切波速 v_S30有相同的量级，并建议完备的场地分类应该同时考虑κ₀和 v_S30的影响。为了便于工程应用中能快速进行κ₀的估计，本研究试图在κ₀和 v_S30之间寻找合适的统计关系，用于近似模拟工程场址的地震动高频衰减。本文拟基于国内外学者研究得到的546个κ₀值，分析κ₀值随v_S30增大的变化趋势，并选取合适的函数形式拟合得到适合于工程应用的κ₀-v_S30模型，进而得到κ₀估计模型，以期为工程场址地震动模拟给出合理的高频衰减系数估计。

1.   地震动模拟中局部场地的高频衰减

随机有限断层法是基于随机理论发展起来的地震动模拟方法，其主要思想是将一个发震规模较大的断层划分为一系列子断层（子源），然后将每一个子断层视为点源，得到每个子源在场址产生的地震动，最终累加所有子源对场址的贡献即可得到整个断层破裂在场址产生的地震动。如果在频域综合考虑震源、传播路径及场地效应的影响，则每个子源产生的傅里叶幅值谱F_A（M₀，f，R）可以写为：

式中： M₀为子源的地震矩； S（M₀，f）表示震源谱； P（R，f）描述传播路径衰减影响；f为频率；R表示震源距；$I(f)= （ 2\text{π}f{\rm{i}} ） ^{n}$为地震动类型因子，当输出为位移、速度和加速度时， n的取值分别为0，1，2；$G （ f ） $描述场地效应影响，可以表示为

式中：A（f）为近地表幅值放大因子，表示研究区域场地对地震谱的放大系数随频率的变化关系；K（f）为与路径无关的高频衰减项，可用kappa滤波器来表示，即

此时$K （ f ） $种对距离的依赖被认为是线性的，表示为

式中，震源距的影响可以通过设置R为0来消除，即κ₀表示地震波在场地下方地质结构中的衰减。Anderson和Hough （1984）认为对于不同观测点的近地表地质条件， κ₀可以表示地震波在观测台站下几百米甚至几千米的衰减，即κ₀是一个与局部场地条件相关的参数。

由于式（1）所示的G（f）主要描述场地效应的影响，因此在地震动的模拟中，式（3）所述的高频衰减主要考虑局部场地的影响，因此代表与路径无关的高频衰减项K（f）可表示为

如前所述，κ₀是一个与局部场地条件相关的参数，已有研究也表明，平均剪切波速 v_S30作为场地条件参数研究基础深厚，其与场地放大效应的相关性显著且得到了学界的普遍认可，作为一个独立参数应用便捷意义清晰，在地震工程领域应用广泛（Castellaro et al，2008；Lee， Trifunac， 2010；Boore et al， 2011）。 因此，本文在后续的研究中，将以 v_S30作为场地条件划分的依据，讨论对局部场地高频衰减系数κ₀的估计。

2.   κ₀数据来源和分析

本文用于统计和分析的κ₀如表1所示，包含了亚洲、欧洲、北美洲和大洋洲的546个数值。其中，平均剪切波速v_S30的变化范围在106.8—2 394.0 m/s之间，κ₀的取值范围在0.002 9—0.087 5 s之间，为了详细分析κ₀随v_S30的变化情况，本文在v_S30＜1500 m/s范围内，以Δv_S30＝100 m/s为间隔，当v_S30大于1 500 m/s时，以300 m/s为间隔，对每个分级档内κ₀的数量、最小值、最大值、标准差和均值进行了统计分析（图1，表2）。可以看出，当v_S30的变化范围在100—1 000 m/s时，虽然数据分布较多，但数据标准差介于0.011—0.018之间，分布相对较为离散；当 v_S30大于1 000 m/s时，数据分布相对较少，标准差也相对较为集中，在0.003 9—0.093 6之间。图2直观地给出了表2中的统计数据的变化情况，图中表示了各组数据的最小值、下四分位数、中位数、上四分位数以及最大值的分布情况。可以看出，虽然κ₀中位数在不同的v_S30分级档内有一定的波动，但总体上呈随着 v_S30的增大逐渐减小的趋势。

表  1  本文研究采用的κ₀数目、来源及相应的v_S分布范围

Table  1.  The number and source of κ₀ used in this study and the corresponding distribution range of v_S

序号	数据个数	vS30/（m·s−1）	地区	数据来源
1	60	106.8—904.2	日本	Cabas等（2017）
2	16	213.2—744.1	日本	Cabas等（2017）
3	50	507.7—1 433.4	日本	van Houtte等（2011）
4	27	1 106.8— 2 394.0	日本	van Houtte等（2011）
5	4	515.7—1 301.3	日本	Laurendeau等（2013）
6	14	170.6—1 428.1	法国	Drouet等（2010）
7	24	192.1—747.1	瑞士	Edwards等（2015）
8	8	1 174.0—1 810.5	瑞士	Edwards等（2011）
9	16	380.1—1 811.5	瑞士	Edwards等（2011）
10	54	160.1—942.8	中国台湾	Huang等（2017）
11	4	233.1—684.8	中国台湾	Lai等（2016）
12	10	167.5—496.4	中国台湾	Lai等（2016）
13	29	191.8—746.9	土耳其	Bora等（2017）
14	16	142.6—1 029.6	意大利	Bora等（2017）
15	5	1 054.7—1 392.5	克罗地亚	Stanko等（2017）
16	4	854.1—953.3	克罗地亚	Stanko等（2017）
17	11	516.6—715.5	中国台湾	van Houtte等（2011）
18	38	299.6—652.9	中国	傅磊和李小军（2017）
19	10	435.7—1 518.3	新西兰	van Houtte等（2018）
20	9	401.0—661.8	亚利桑那	Kishida等（2014）
21	6	550.9—1 000.3	加利福尼亚	van Houtte等（2011）
22	117	170.9—1 531.2	中国台湾	Chang等（2019）
23	7	263.3—660.3	中国	郑旭等（2019）
24	7	531.2—912.1	日本	朱百慧（2016）

下载: 导出CSV 

| 显示表格

图  1  本文研究采用的κ₀值随平均剪切波速 v_S30的分布

Figure  1.  The distribution of κ₀ value used in this study with mean shear wave velocity v_S30

下载: 全尺寸图片 幻灯片

表  2  κ₀值在不同 v_S30范围的分组统计

Table  2.  Group statistics of κ₀ values in different v_S30 ranges

vS30/（m·s−1）	κ0值
	数量	最小值	最大值	标准差	均值
100—200	37	0.019 6	0.076 7	0.014 85	0.054 48
200—300	65	0.022 2	0.072 6	0.014 58	0.051 33
300—400	95	0.009 2	0.077 2	0.015 38	0.043 41
400—500	95	0.008 8	0.087 5	0.016 06	0.040 43
500—600	73	0.003 8	0.067 4	0.015 16	0.035 60
600—700	46	0.003 8	0.080 9	0.015 17	0.029 28
700—800	24	0.013 0	0.071 1	0.014 11	0.033 48
800—900	19	0.004 9	0.073 9	0.017 84	0.034 12
900—1000	16	0.014 3	0.052 8	0.010 99	0.027 61
1 000—1 100	7	0.015 3	0.027 0	0.003 93	0.022 73
1 100—1 200	7	0.006 0	0.025 8	0.007 52	0.015 47
1 200—1 300	10	0.013 6	0.026 0	0.004 61	0.019 67
1 300—1 400	6	0.010 1	0.031 5	0.009 06	0.020 13
1 400—1 500	9	0.009 6	0.026 3	0.006 15	0.015 72
1 500—1 800	16	0.006 9	0.027 1	0.006 75	0.016 06
1 800—2 100	11	0.002 9	0.029 3	0.009 36	0.013 71
2 100—2 400	10	0.002 9	0.025 0	0.009 30	0.011 93

下载: 导出CSV 

| 显示表格

图  2  κ₀值统计数据箱型图

Figure  2.  κ₀ value statistics box plot

下载: 全尺寸图片 幻灯片

3.   κ₀变化特点及模型选择

由图2可以看出，虽然整体数据分布较为离散，但κ₀的中位数值随平均剪切波速 v_S30的增大而减小的趋势还是很明显的。为了更加详细地分析κ₀随 v_S30的变化，本文选用滑动平均法对数据进行了处理。滑动窗长度为200 m/s，步长为20 m/s，按照式（6）所示的方法计算每个窗内全部数据的均方根值κ_rms。图3中三角形是按照式（6）计算得到的115个均方根值，可以很清晰地看出，κ_rms呈现出随着 v_S30的增加逐渐减小的趋势。

图  3  四种函数对κ_rms的拟合情况

Figure  3.  Fitting of four functions to κ_rms

下载: 全尺寸图片 幻灯片

式中：κ_i表示不同滑动窗内的κ₀值，N为κ₀的个数.

为了拟合κ₀随v_S30的变化，本文选用了四个常用的函数形式，即线性函数、多项式函数、对数线性函数和双对数线性函数，对上述115组数据进行拟合，拟合曲线如图3所示。函数表达形式和拟合结果列于表3。分析四个函数的拟合结果可知，相比于线性和双对数线性函数，多项式和对数线性函数拟合的和方差SSE最接近于0，拟合优度R²最接近于1，拟合效果相对更为理想。但从图3所示的曲线变化趋势来看，多项式拟合结果在 v_S30＞2 000 m/s的时候有逐渐增大的趋势，这和本文对κ₀变化的认识不符。考虑到对数线性函数拟合无论在拟合优度和曲线的变化趋势上，都能很好地模拟κ₀随 v_S30的变化，与整体数据变化趋势最为接近，因此本文选用对数线性函数作为模型形式对原始数据进行拟合。

表  3  κ_rms与 v_S30的经验关系

Table  3.  The empirical relationship between к_rms and v_S30

拟合函数	模型参数
	a	b	SSE	R2	c
线性函数：κrms＝avS30＋b	−1.557×10−5	4.478×10−2	3.61×10−3	7.734×10−1
多项式函数：κrms＝${av^2_{{\rm{S}}30}} $＋bvS30＋c	1.284×10−8	−4.741×10−5	6.559×10−4	9.587×10−1	5.886×10−5
对数线性函数：κrms ＝algvS30＋b	−3.439×10−2	1.286×10−1	1.466×10−3	9.07 8×10−1
双对数线性函数：lgκrms＝algvS30＋b	−4.488×10−1	−2.72×10−1	2.245×10−3	8.587×10−1
注：表中a，b和c为模型拟合参数，SSE 表示和方差，R2为拟合优度。

下载: 导出CSV 

| 显示表格

4.   基于v_S30的κ₀估计模型

4.1   模型参数拟合及实用性分析

由第3节的讨论可知，κ₀与平均剪切波速 v_S30之间呈现出很强的负相关性，本文采用式（7）所示的对数线性函数作为κ₀-v_S30的模型函数，即

为了能够合理地对模型参数a和b进行拟合，本文对原始数据的分布情况进行分析。数据整体分布情况如图4所示，可以看出数据整体呈近似对数正态分布形式，因此本文采用最小二乘法对式（7）所示的模型函数进行参数拟合。最小二乘法是一种数学优化技术，它通过最小化误差的平方和来寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小，从而保证拟合模型与实际数据之间的误差最小，使所得到的拟合模型符合要求。需要说明的是，由于部分原始数据偏离样本均方根值较远，为了保证拟合参数的合理性，本文对原始数据以偏离统计均方根值的程度进行了适当的筛选，最后选取477个原始数据用于模型拟合。

图  4  数据量分布情况图

Figure  4.  Data distribution diagram

下载: 全尺寸图片 幻灯片

利用最小二乘法，对选取的477个原始数据进行如式（7）所示的模型参数拟合，拟合得到的模型为

模型拟合得到的和方差SSE＝0.066 58，拟合优度R²＝0.489 9，标准差σ＝0.011 8。由于和方差接近于0，标准差也相对较小，说明对数线性函数模型对原始数据拟合的效果较为理想。式（8）所示模型对原始数据的拟合情况如图5所示，可以看出，拟合曲线能够较好地穿过原始数据，且能够较好地表示κ₀随 v_S30的变化趋势。此外，图5中列出了三个学者提出的κ₀模型：van Houtte （2011）模型是基于日本地区数据提出的，通过比较可以看出，虽然在 v_S30＞700 m/s部分与本文的模型较为贴近，但模型在v_S30较小的范围内对κ₀的估值偏大；模型Edwards （2011） 是基于瑞士地区数据提出的，v_S30在100—1 500 m/s范围内对κ₀的估计普遍偏低，在v_S30＞1 500 m/s时与本文得到的结果相近；朱百慧 （2016） 模型也是基于日本地区数据提出的，该模型在v_S30＜500 m/s部分与本文研究结果较为符合，但在v_S30大于500 m/s的范围内，对κ₀的估计明显偏低。

图  5  对数线性函数拟合以及与其它模型对比

Figure  5.  Log-linear function fitting and comparison with other models

下载: 全尺寸图片 幻灯片

为进一步探究对数模型的拟合效果，本文采用式（9）对拟合结果进行残差分析，从而确定拟合结果的优劣，即：

式中，ε表示由估计值与实测值确定的残差，κ_p表示κ₀的估计值，κ_t表示κ₀的实测值。

图6所示为计算得到的残差ε随着 v_S30的变化，可以看出，残差在零线附近 ［ 0，0.02 ］范围内呈正态分布，结果介于 ［ −0.025，0.025 ］ 之间，进一步说明拟合结果较为理想。

图  6  对数线性函数拟合残差图

Figure  6.  Log-linear Function fit residual plot

下载: 全尺寸图片 幻灯片

为了检验式（8）所示的模型对整体数据的包容情况，本文分析了原始数据在拟合均值±1倍标准差范围内的分布情况，计算结果表明，76.7%的原始数据点落在拟合均值±σ区间范围内（图7）。对结果进一步分析表明，本文建立的模型能够较好地表征κ₀与v_S30之间的关系，拟合模型的标准差取值也能较好地表示数据的离散程度，式（8）所示的模型可用来估计κ₀均值。

图  7  整体数据校验和MTN台站数据校验情况

Figure  7.  Overall data check and MTN station data check

下载: 全尺寸图片 幻灯片

为了检验式（8）所示模型的适用性，本文选用了根据中国台湾地区TSMIP强震台站（包括48个MTN强震台站）1989—2020年的48条强震记录计算得到的κ₀值，v_S30为 391.0—954.1 m/s的范围之间，κ₀在0.026 1—0.055 6 s的范围之间变化。图7三角形数据点所示为48个κ₀数据落在本文建立的模型中的情况，可以看出，48个数据中有91.7%的数据点落在κ₀均值±σ范围之内，进一步说明本文给出的κ₀-v_S30模型可以用来合理地估计地震动的高频衰减。

4.2   对κ₀估计的讨论

虽然本文得到了一个可以近似估计局部场地地震动高频衰减系数的模型，但是在进行数据的统计分析中，还存在一些问题需要进一步讨论。

1） 数据区域性的影响. 在本文的统计分析中，没有区分数据来源区域，将亚洲、欧洲、北美洲、大洋洲等不同区域的数据放在一起进行统计分析。分析图1所示的分类数据和图7所示的κ₀均值的估计值，可以看出，κ₀的估计具有一定的区域差异性，例如Drouet 等（2010） 研究结果中的法国数据分布于拟合曲线下方，van Houtte等 （2011）研究采用的日本数据则相对均匀地分布在拟合曲线上下两侧，Bora等（2017）研究中给出的意大利数据集中分布于拟合曲线下方，而Chang等（2019）研究中给出的中国台湾数据则集中分布于拟合曲线上方，由此可见数据的区域性分布较为明显，对κ₀的估计也会产生一定的影响。

2） 场地条件的影响. 本文给出的κ₀的估计模型仅考虑其与平均剪切波速v_S30之间的关系，由于κ₀描述地震波在观测台站下几百米甚至几公里的衰减，因此，对κ₀取值的影响可能会与场址下方的介质相关。此外，对v_S30的计算方法可能也有影响，v_S30和场地效应之间的联系是统计性的，而不是物理性的直接预测，因而不同的计算方法会产生不同的结果，从而间接影响κ₀的估计。

3） κ₀计算方法的影响. κ₀计算方法有很多种，不同的计算方法会导致结果的偏差。由于本文的研究数据来源于不同学者的研究结果，在κ₀计算方法上有差异，从而导致κ₀计算结果离散性较大，最终影响κ₀的估计。

4） 模型选择的影响. 由于回归模型本身具有多样性，不同回归模型得到的κ₀的估计可能会有差异。

5.   讨论与结论

本文选取了国内外学者计算得到的对应不同场地条件的546个高频衰减系数κ₀，运用最小二乘法进行回归分析，研究了κ₀均值与场地平均剪切波速v_S30之间的关系，得到了如下主要结论：

1） 运用均值法和滑动平均法对κ₀值随平均剪切波速v_S30的变化进行了分析，结果表明，虽然κ₀具有明显的区域差异性，但随着v_S30的增大，呈现出逐渐减小的趋势，且这种变化趋势可以用合适的函数进行模拟；然后采用线性函数、多项式函数、对数线性函数和双对数线性函数对κ₀的均方根值进行拟合，通过分析和比较，最终选取对数线性函数作为κ₀的拟合模型。

2） 基于筛选得到的477个数据，采用最小二乘法对模型参数进行了拟合，得到了适合于工程应用的κ₀-v_S30实用模型。对模型的适用性分析表明，本文构建的模型在预测场址地震动参数时能够合理估计地震动的高频衰减。

3） 本文通过对拟合模型的和方差、拟合优度、标准差和残差分析，验证了该模型的可行性。并基于中国台湾地区MTN强震台站的κ₀落在模型中的情况检验了κ₀-v_S30模型的适用性。结果表明，本文构建的模型能够合理估计不同场地条件下的κ₀均值，因此，在预测场址地震动参数时，可用本文构建的κ₀-v_S30模型估计场址的κ₀系数，进而近似模拟地震动的高频衰减效应。

需要特别指出的是，虽然在模型的构建中，在数据区域性、场地条件、κ₀计算方法以及模型选择及回归模型等方面都存在不确定性，并且数据的离散性较大，但总体来看，κ₀随 v_S30增大逐渐减小的趋势符合地震动在局部场地的衰减特征，并且在±1倍标准差的范围内，可以大致描述参数的变化特征，从而能够实现κ₀的近似估计。