范颖霏 姜云木 刘章军

摘要：工程上常将规范反应谱拟合的确定性取值作为地震动演变功率谱参数，这虽然便于工程应用，但同时存在明显弊端：一方面，参数并非从实测强震记录中获得，具有较强的经验性;另一方面，确定性参数生成的地震动样本过于规则，且工程特性单一，无法全面地反映地震动的随机性。为了克服上述难题，从PEER的NGAGWest2地震动数据库中筛选1 766条主轴方向的实测强震记录，并根据《中国地震动参数区划图》建议的场地类别和聚类方法划分为15组;随后，识别实测强震记录的演变功率谱参数，并结合K-S检验和BIC信息准则，确定每个参数最优概率模型;最后，根据演变功率谱参数统计建模结果，结合降维方法，生成了Ⅱ类场地与Ⅲ类场地的代表性时程。与《建筑抗震设计规范》中建议取值不同，研究参数来源于实测地震动，且具有随机性，避免地震动样本工程特性的单一化和规则化。

关键词：参数识别; 全随机性地震动; 实测记录; 统计建模; 降维模拟

中图分类号： P315.9;TU311.3      文献标志码：A   文章编号： 1000-0844（2024）01-0126-11

DOI：10.20000/j.1000-0844.20220130001

Statistical modeling of evolutionary power spectral density parameters of ground motions

Abstract： In engineering practices， the deterministic values of the response spectrum in accordance with seismic design codes are generally applied as evolutionary power spectral density （EPSD） parameters of ground motions. Although this treatment is convenient for engineering applications， it comes with two evident disadvantages. On the one hand， parameters are not obtained from measured strong motion records， which makes them highly empirical. On the other hand， the ground motion samples generated by deterministic parameters are extremely regular and exhibit unitary engineering characteristics， which cannot fully reflect the randomness of ground motions. To overcome the abovementioned challenges， we selected 1 766 strong motion records in the main axis direction from PEER and divided them into 15 groups based on the site categories suggested by the Seismic Ground Motion Parameters Zonation Map of China. Subsequently， the EPSD parameters of strong motion records were identified. Combined with the K-S test and BIC information criterion， the optimal probability models for each parameter were determined. Finally， based on statistical modeling results of EPSD parameters， representative time histories of site categories II and III were generated through the dimension reduction method. In contrast to the recommended values from the Code for Seismic Design of Buildings， the parameters suggested in this study were derived from measured records and possessed significant randomness， thus avoiding the simplification and regularization of engineering characteristics of ground motion samples.

Keywords：parameter identification; fully stochastic ground motion; measured records; statistical modeling; dimension reduction simulation

0 引言

相對于震源机制等地震学因素和“震源-传播途径-局部场地”的物理模型，工程上更加关注地震动本身的工程特性。演变功率谱模型从地震动的幅值、持时、频谱等工程特性出发，是地震工程中应用最广泛的随机地震动模型。一般来说，演变功率谱模型由强度调制函数和平稳功率谱两部分组成，强度调制函数对平稳信号进行调幅;平稳功率谱描述局部场地作用在白噪声过程下的能量分布规律。在确定了演变功率谱模型之后，便可用谐波合成法［1］模拟地震动加速度过程。

目前演变功率谱的常用参数是根据规范反应谱拟合得到的确定性取值，这虽然便于工程应用，但是忽略了地震震级、传播途径、场地条件等因素的不确定性，从而忽略了演变功率谱参数显著的随机性。这将导致生成的地震动样本过于规则，且工程特性单一，不能反映地震动丰富的概率信息和细部概率结构。

对于参数取值的确定性研究，已有较多学者以抗震规范为基础，取得了较为完善的成果［2-5］。对于参数的随机性建模，丁艳琼等［6］基于工程随机地震动物理模型，对7 000条地震动进行了精细的分组与参数识别，并对其进行了统计建模;李杰等［7］同样基于工程随机地震动物理模型，对5 000条地震动进行了参数识别与统计建模。虽然随机性参数与功率谱的概念似有矛盾，但却更符合工程上对地震动认知与应用。

为了体现震级、距离以及场地因素对地震动造成的影响，本文对实测强震记录进行了精细的分组。利用最小二乘法，识别演变功率谱的模型参数;并结合拟合优度检验和贝叶斯信息准则，确定每个参数的最优概率模型;最后引入随机函数约束［8-10］，仅用6个随机变量便可实现地震动过程在概率层面上的模拟。由于地震动代表性时程具有赋得概率，可进而与概率密度演化理论相结合［11-12］，对结构进行精细的动力响应和可靠度分析。

1 地震动工程随机模型及人工合成方法

根据文献［13］，非平稳地震动过程的演变功率谱密度函数可表示为：

SUg（ω，t;λUg）=f（t;λf）2S（ω;λS） （1）

式中：SUg（ω，t;λUg）为非平稳地震动加速度过程Ug（t）的单边演变功率谱;f（t;λf）为强度调制函数;S（ω;λS）为单边功率谱。

对于强度调制函数，采用工程上常用的Amin-Ang模型［14］。该模型优点是可以较为全面地刻画地震动的波形;缺点是采用了分段形式，形式复杂且无法准确地反映地震动峰值到达时刻［15］。

式中：t1、t2和c分别代表非平稳地震动过程的平稳段起始时间、平稳段结束时间和衰减段的衰减速度。因此，Amin-Ang模型的参数向量可表示为λf=（t1，t2，c）。

對于单边功率谱，采用工程上常用的Kanai-Tajimi模型［16］。该模型优点是从动力学角度推导，参数具有明确的物理意义;缺点是S（0;λS）≠0，不符合实际情况：

式中：ωg和ξg分别为场地土的卓越圆频率和阻尼比;S0为白噪声激励的强度，可表示为［17］：

式中：Amax表示地震动峰值加速度;r为峰值因子;ωe为S0=1时S（ω;λS）与坐标轴围成的面积。因此，Kanai-Tajimi模型的参数向量可以表示为λS=（ωg，ξg，r）。

这样，演变功率谱密度函数的参数向量λUg便可表示为：

λUg=（λf λS）=（t1，t2，c，ωg，ξg，r） （5）

确定演变功率谱密度函数SUg（ω，t;λUg）之后，便可采用谐波合成法生成地震动过程。根据文献［10］，非平稳地震动过程的源谱表达为：

式中：Δω为频率步长;Δω=ωu/N;N为截断项数;ωk=k×Δω（k=1，2，…，N）;ωu为截断频率。

式（6）中，Xk，Yk为一组零均值的标准正交随机向量，满足以下基本条件：

E［Xk］=E［Yk］=0，E［XjYk］=0

E［XjXk］=E［YjYk］=δjk （7）

式中：E［·］为数学期望;δjk为Kronecker-delta记号。

传统的谐波合成法在本质上属于Monte Carlo方法：一方面，生成的样本数量巨大，极大地增加了结构动力反应分析的计算量;另一方面，Monte Carlo方法在本质上属于随机抽样方法，所生成的样本概率信息不完备，无法进行结构精细化动力反应分析和动力可靠性评价。为克服传统Monte Carlo方法的上述挑战，本文引入了随机函数的降维思想［8-10］，将标准正交随机变量集{Xk，Yk}（k=1，2，…，N）定义为基本随机向量的正交函数形式：

2 实测强震记录

本文从太平洋地震动工程研究中心（PEER）NGA-West2地震动数据库中筛选了1 766条主轴方向的实测强震记录，且严格遵循以下原则挑选：

（1） 断层距离应大于10 km，以减少近断层地震动的数量。

（2） 实测记录的矩震级应大于5，以排除对结构影响较小的地震动。

在筛选的实测强震记录中，Rjb（Joyner-Boore距离）范围为0～540 km;vS30（地下30 m平均剪切波速）范围为106.83～1 525.85 m/s;M（矩震级）范围为：5.0～7.9。

《中国地震动参数区划图》［18］中建议了Ⅰ0～Ⅳ五种场地类别。本文结合文献［9］，根据每条实测强震记录的vS30对实测强震记录进行分类。同时，为了进一步体现震级与距离因素对地震动工程特性造成的影响，本文以M和Rjb为指标，采用K均值聚类的方法，将每个场地的实测强震记录另分为3组，可概括为：远场大震、近场小震和近场大震［6］。由于小震无法传播到较远的距离，因此没有远场小震这一项。

Ⅱ类场地的聚类结果如图1所示。可以看出，聚类方法的分组结果与《建筑抗震设计规范（GB 50011—2010）》［19］中的地震分组有着相似性。

表1给出了不同聚类分组和场地对应的实测强震记录数量：

此外，实测强震记录还需要四阶Butterworth滤波处理以及1%～99%范围的能量截取，以保证研究结果的可靠性［20］。

图2和图3分别给出了不同场地的反应谱均值对比图和Ⅱ类场地不同分组的反应谱均值对比图。需要说明的是，为了对比清晰，将实测强震记录均调幅至0.2g·m/s2。由图2可知，随着场地土变软，反应谱峰值逐渐向长周期移动，且长周期成分也逐渐增多。由图3可知，不同分组的反应谱峰值和长周期成分有着显著差异，充分证明了震级和距离对地震动工程特性的显著影响。

3 演变功率谱参数识别

对于第i条实测记录ai（t），其随时间变化的归一化能量曲线Ii（t）可表示为［21］：

式中：Ti为第i条实测强震记录的持时。

对于非平稳地震动过程，其随时间变化的归一化模型能量曲线P（t）为：

以Ii（t）为目标值，采用最佳平方逼近原则，便可识别ai（t）的参数向量λf，i：

根据帕萨瓦尔定理，信号在时域上的能量与频域上的能量相等，可以得到：

令

于是：

同时结合式（3）和式（4）得：

由于ai（t）的峰值加速度Amax，i是已知的，将参数向量λf，i代入式（14）中便可得到Ei，再根据式（15b）便可得到与ai（t）对应的峰值因子ri。

对于平稳地震动功率谱的场地参数，采用拟合反应谱的方法识别。为简便起见，在对场地参数进行识别时，将地震动视为等效平稳过程。Vanmarcke将随机过程的反应谱定义为单质点体系反应峰值系数的平均值与反应方差的乘積［22］，因此，地震动反应谱与功率谱转换公式可以定义为：

式中：（ω0;λf，i）为等效平稳过程的峰值因子;σ（ω0，ξ;λS，i）为等效平稳过程的反应方差。ω0与ξ分别为结构的固有圆频率和阻尼比，在本文中ω0≥1.05 rad/s，ξ=0.05;Td，i为等效平稳过程的持续时间，即为强度超过峰值50%的持续时间，对于Amin-Ang模型，Td，i的表达式为：

以ai（t）前6 s反应谱αi（ω0）为目标值，根据式（16），采用最佳平方逼近原则，对参数向量λS，i进行识别：

这样，便完成了与ai（t）对应的演变功率谱参数向量λUg，i的识别。类似地，可依次识别1 766组实测强震记录的参数向量λUg，i。

以SFERN地震为典型实例，根据式（11）与式（16），反应谱和地震动归一化能量的结果见图4（a）

与图4（b），可以看出，模型与实测记录均拟合良好，验证了本文采用的识别方法的有效性。参数识别结果如表2所列。

4 演变功率谱参数统计建模

为了体现演变功率谱参数的随机性，本文选取了5种备选概率模型对每个参数的概率分布进行拟合。这5种备选概率模型分别为：对数正态分布（LOGN）、耿贝尔分布（GUM）、广义极值分布（GEV）、韦布尔分布（WEI）、伽马分布（GAM），它们的概率密度函数如表3所列。

以Ⅱ类场地的参数ωg为例，概率分布拟合结果如图5。可以看出，难以直接判断ωg的最优概率模型。

为此，本文引入了K-S检验。一般来说，当K-S检验为0时，表示接受该备选概率模型;为1时，表示拒绝该备选概率模型。进一步地，当有多个备选概率模型通过K-S检验时，则根据贝叶斯信息准则进一步筛选。贝叶斯信息准则可表示为［24］：

BIC=ln（n）u-2ln（L） （21）

式中：u为概率模型中参数的个数;L为模型最大似然函数;n为样本容量。一般地，BIC值越小，代表模型拟合度越好。

结合K-S检验与BIC信息准则，每个演变功率谱参数的最优概率模型以及对应的模型参数如表4所列。

在表4中，当概率模型为LOGN时，Par1为均值，Par2为标准差;概率模型为GUM时，Par1为位置参数，Par2为形状参数;概率模型为GEV时，Par1为形状参数，Par2为尺度参数，Par3为位置参数;概率模型为WEI时，Par1为比例参数，Par2为形状参数;概率模型为GAM时，Par1形状参数，Par2为尺度参数。

为了更清晰地表现参数随场地变化的趋势，图6展示了每个参数均值不同场地的变化趋势图。

图6（a）可知，随着土质变软，场地土卓越圆频率逐渐变小;而峰值因子和场地土阻尼比没有明显的变化。对于场地土阻尼比，文献［6］指出，阻尼比不仅与场地土的软硬有关，还与其厚薄、孔隙率、形状、体积等因素有关。因此不能通过简单的场地类别与分组反映其影响的强弱。对于峰值因子，文献［3］指出，峰值因子对场地变化不敏感。

图6（b）可知，随着土质的变弱，地震动的平稳段起始时间和地震动衰减因子逐渐缩小，而地震动的平稳段持续时间逐渐延长。这说明随着土质的变弱，地震动的总持时逐渐延长。

综上，建议将场地土阻尼比和峰值因子近似做均值处理，而其他参数做随机变量处理。

5 数值算例

5.1 降维模拟的实现

由前文可知，模拟地震动共需要6个随机变量{t1，t2，c，ωg，Xk，Yk}。且根据式（8），Xk，Yk均是同一个基本随机变量的函数形式。这样，模拟地震动过程共需要5个基本随机变量。

式中：l=1，…，nsel;F-1ωg、F-1t1、F-1t2及F-1c分别代表任意分组下ωg、t1、t2及c最优累计分布函数的反函数。这样，便得到了具有随机性的演变功率谱模型参数。

式中：l=1，…，nsel。

这样，便得到了地震动过程的基本随机变量集Θ={Θ1，Θ2，Θ3，Θ4，Θ5}，将基本随机变量Θ1～Θ4，场地土阻尼比均值和峰值因子均值代入式（1）中，便可得到非平稳地震动的演变功率谱模型SUg（ω，t;λUg）;将基本随机变量Θ5代入式（8）中，得到随机变量{Xk，Yk}。再根据式（6），便可模拟地震动样本。

5.2 数值模拟

在本文算例中，以Ⅱ类和Ⅲ类场地为例。根据表4可知，场地土阻尼比ξg和峰值因子r的取值如表5所列。

其他模拟参数为：样本数量144;频率截断项数N=1 600;截断频率ωu=240 rad/s;频率步长Δω=0.15 rad/s;模拟持时40 s ;时间步长Δt=0.01 s。

图7分别给出了Ⅱ类和Ⅲ类场地的代表性样本。可以看出，样本在频谱和持时上均具有明显的随机性，更符合实际地震动。

图8分别给出了Ⅱ类和Ⅲ类场地的模拟值均值和标准差与各自目标值的对比结果。可以看出，模拟值与目标值对应效果良好，证明了降维方法的有效性。为了进一步证明降维方法的精确性，表6中分别给出了Ⅱ类场地144、377和987条样本的模拟误差，可以看出随着样本数量的增加，均值及标准差误差逐渐变小，体现了降维方法的收敛性;987条样本的模拟误差小于5%，符合工程要求，体现了降维方法的工程适用性。

图9分别给出了Ⅱ类和Ⅲ类场地的模拟反应谱均值与实测反应谱均值的对比。可以看出，实测反应谱均值在模拟反应谱均值的±1倍标准差范围内，且与其拟合较为一致，证明了本文参数模型的正确性和有效性。

6 结论

本文基于实测强震记录，对演变功率谱模型进行了统计建模，得到以下结论：

（1） 为了反映震级以及传播距离对地震动参数的影响，对实测强震记录进行精细化分组，通过对比各组反应谱均值可知，震级以及传播距离等因素对地震动有直接且显著的影响。

（2） 对演变功率谱参数识别与统计建模可知，震级、传播距离等不确定因素均间接地体现在了演变功率谱参数上。因此，有必要考虑演变功率谱参数的随机性。

（3） 本文的参数模型来源于实测强震记录，且具有显著的随机性，克服了传统经验化和固定化的参数取值而导致地震动样本工程特性单一的缺陷;同时，采用降维方法生成的样本具有良好的效率和精度。因此，将本文的参数模型和降维方法相结合，可高效且精确地模拟出能够全面反映地震动工程特性的样本。

参考文献（References）

［1］ SHINOZUKA M，JAN C M.Digital simulation of random processes and its applications［J］.Journal of Sound and Vibration，1972，25（1）：111-128.

［2］ 任文杰，赵利军，穆蒙蒙.现行抗震规范下地震动随机模型参数研究［J］.世界地震工程，2017，33（2）：33-38.

REN Wenjie，ZHAO Lijun，MU Mengmeng.Study on parameters of random earthquake ground motion model based on current seismic code［J］.World Earthquake Engineering，2017，33（2）：33-38.

［3］ 张治勇，孙柏涛，宋天舒.新抗震规范地震动功率谱模型参数的研究［J］.世界地震工程，2000，16（3）：33-38.

ZHANG Zhiyong，SUN Baitao，SONG Tianshu.Study on the parameters of seismic power spectrum model based on the new seismic code［J］.World Information on Earthquake Engineering，2000，16（3）：33-38.

［4］ 白國良，朱丽华.基于现行抗震规范的Kanai-Tajimi模型参数研究［J］.世界地震工程，2004，20（3）：114-118.

BAI Guoliang，ZHU Lihua.Study on the parameters of Kanai-Tajimi model based on the code （GB 50011-2001）［J］.World Information on Earthquake Engineering，2004，20（3）：114-118.

［5］ 欧进萍，朱荻涛.设计用随机地震动的模型及其参数确定［J］.地震工程与工程振动，1991，11（3）：45-54.

OU Jinping，ZHU Ditao.Random earthquake ground motion model and its parameter determination used in aseismic design［J］.Earthquake Engineering and Engineering Vibration，1991，11（3）：45-54.

［6］ 丁艳琼，李杰.工程随机地震动物理模型的参数识别与统计建模［J］.中国科学：技术科学，2018，48（12）：1422-1432.

DING Yanqiong，LI Jie.Parameters identification and statistical modelling of physical stochastic model of seismic ground motion for engineering purposes［J］.Scientia Sinica （Technologica），2018，48（12）：1422-1432.

［7］ 李杰，王鼎.工程随机地震动物理模型的参数统计与检验［J］.地震工程与工程振动，2013，33（4）：81-88.

LI Jie，WANG Ding.Parametric statistic and certification of physical stochastic model of seismic ground motion for engineering purposes［J］.Journal of Earthquake Engineering and Engineering Vibration，2013，33（4）：81-88.

［8］ LIU Z J.Dimension reduction of Karhunen-Loeve expansion for simulation of stochastic processes［J］.Journal of Sound and Vibration，2017，408：168-189.

［9］ 刘章军，刘子心，阮鑫鑫，等.地震动随机场的POD降维表达［J］.中国科学：技术科学，2019，49（5）：589-601.

LIU Zhangjun，LIU Zixin，RUAN Xinxin，et al.POD-based dimension reduction representation of stochastic ground motion fields［J］.Scientia Sinica （Technologica），2019，49（5）：589-601.

［10］ LIU Z J.Random function based spectral representation of stationary and non-stationary stochastic processes［J］.Probabilistic Engineering Mechanics，2016，45：115-126.

［11］ LI J，CHEN J B.Stochastic dynamics of structures［M］.Chichester，UK：John Wiley & Sons，Ltd，2009.

［12］ LI J.The principle of preservation of probability and the generalized density evolution equation［J］.Structural Safety，2008，30（1）：65-77.

［13］ 刘章军，刘子心.基于规范反应谱的全非平稳地震动过程模拟［J］.振动工程学报，2017，30（3）：457-465.

LIU Zhangjun，LIU Zixin.Simulation of fully non-stationary ground motion based on seismic design response spectrum［J］.Journal of Vibration Engineering，2017，30（3）：457-465.

［14］ AMIN M，ANG A H S.Nonstationary stochastic models of earthquake motions［J］.Journal of the Engineering Mechanics Division，1968，94（2）：559-584.

［15］ BOGDANOFF J L，GOLDBERG J E，BERNARD M C.Response of a simple structure to a random earthquake-type disturbance［J］.Bulletin of the Seismological Society of America，1961，51（2）：293-310.

［16］ KANAI K，YOSHIZAWA S，SUZUKI T.An empirical formula for the spectrum of strong earthquake motions Ⅱ［J］.東京大学地震研究所彙報，1961，41：261-270.

［17］ CLOUGH R W，PENZIEN J.Dynamic of structures［M］.New York：McGraw-HILL Inc，1993.

［18］ 国家质量监督检验检疫总局，中国国家标准化管理委员会.中国地震动参数区划图：GB 18306—2015［S］.北京：中国标准出版社，2016.

General Administration of Quality Supervision，Inspection and Quarantine of the People's Republic of China，Standardization Administration of the People's Republic of China.Seismic ground motion parameters zonation map of China：GB 18306—2015［S］.Beijing：Standards Press of China，2016.

［19］ 中華人民共和国住房和城乡建设部，国家质量监督检验检疫总局.建筑抗震设计规范：GB 50011—2010［S］.北京：中国建筑工业出版社，2010.

Ministry of Housing and Urban-Rural Development of the People's Republic of China，General Administration of Quality Supervision，Inspection and Quarantine of the People's Republic of China.Code for seismic design of buildings：GB 50011—2010［S］.Beijing：China Architecture & Building Press，2010.

［20］ ASCE.Minimum design loads for buildings and other structures （ASCE/SEI 7-10）［S］.Reston：American Society of Civil Engineers，2010.

［21］ VANMARCKE E H，LAI S S P.Strong-motion duration and RMS amplitude of earthquake records［J］.Bulletin of the Seismological Society of America，1980，70（4）：1293-1307.

［22］ VANMARCKE E H.Structural response to earthquakes[J].Developments in Geotechnical Engineering，1976，15：287-337.

［23］ KOTZ S，NADARAJAH S.Extreme value distributions：theory and applications［M］.Singapore：Upper Saddle River：Prentice Hall，2000.

[24] SCHWARZ G.Estimating the dimension of a model[J].The Annals of Statistics，1978，6（2）：461-464.

［25］ SHIELDS M D，ZHANG J X.The generalization of Latin hypercube sampling［J］.Reliability Engineering & System Safety，2016，148：96-108.