Research on application of focal mechanism and site conditions in judgment of Hebei earthquake influence field
-
摘要: 震后应急工作中地震影响场的判定和快速给出较为合理的地震烈度分布图，是震后应急救援的重要依据，对于政府了解灾情、部署工作以及估算灾害损失都尤为重要，所以本文以此为研究目的，力求震后快速给出准确的地震烈度分布图。本文收集整理了河北地区中强地震的实际等震线图，将其与加入了震源机制解影响参数的烈度衰减关系计算得到的理论等震线图进行对比。结果显示：随着震级的增大，由衰减关系计算得到的等震线图与实际地震等震线图在高烈度区（≥Ⅶ度）相似度更高。另外，根据震后24小时内余震频度的空间变化，对极震区理论等震线修正后，其与实际等震线更加贴合，即理论计算烈度与实际调查烈度值更加接近。最后，对河北地区划分网格，根据地震动衰减关系计算震例对各个网格中心点产生的影响—基岩PGA。提取场地类别属性，考虑场地放大因子，完成基岩PGA到地表PGA的转换。将地表PGA换算成烈度，并与实际地震等震线图进行对比分析。结果表明，考虑了场地放大效应的地震影响场在高烈度区与实际等震线相似度很高，且相似度超过基于震源机制解的烈度衰减关系方法。Abstract: The determination of seismic influence field in post-earthquake work, a quickly gived reasonable map of seismic intensity distribution, was an important basis for emergency rescue after the earthquake and was important for the government to understand the disaster situation, deploy work and estimate the disaster loss. In this paper, the isoseismal maps of moderate strong earthquakes in Hebei Province were collected and sorted out, it was compared with the theoretical isotherm map generated by the regional earthquake intensity attenuation relationship based on the focal mechanism solution. The results show that, with the increase of magnitude, the similarity was higher between the theoretical isoseismal map and the actual isoseismal map in the high intensity area （≥Ⅶ）, the theoretical isoseismal map was calculated by the attenuation relationship of seismic intensity with focal mechanism solution. In addition, according to the spatial variation of aftershock frequency within 24 hours after the earthquake, the theoretical isoseismal line in the polar region is modified, which is more consistent with the actual isoseismal line. That is to say, the theoretical calculated intensity is closer to the actual investigation intensity value. Finally, the grid of Hebei area was divided, and the bedrock PGA of the earthquake case on each grid center point was calculated according to the attenuation relationship of ground motion. Then, the site category attributes were extracted, considered site amplification factor, and the conversion of bedrock PGA to surface PGA was completed. The surface PGA was converted into intensity and compared with the actual earthquake isoseismal map. The results show that the similarity between the seismic influence field calculated by considering the site amplification effect and the actual isoseismal line was very high in the high intensity area, and the similarity is higher than the intensity attenuation relation method based on the focal mechanism solution.
-
引言
极端事件在随机现象的研究过程中出现频率很低，比如大地震、洪水、飓风等，但这些事件一旦发生，将会给生产或生活带来巨大的影响，因此对数据样本中的极端变异性数据信息进行深入研究成为众多学者的共识，在此基础上形成的极值统计分析方法成为灾害损失评估与风险评价的重要工具。
von Bortkiewicz （1922）是近代第一个明确提出极值问题的学者，其在研究正态分布样本的极差时，发现来自正态分布的样本最大值具有新的分布；von Mises （1923）首次对正态样本极值的渐近分布进行了研究；Fréchet （1927）的研究表明，来自不同分布但有某种共同性质的最大值可以有相同的渐近分布；Fisher和Tippett （1928）将次序统计量规范化，按中心极限定理的思想，提出了极值分析渐近原理的基础性定理，完成了极值统计分析基础理论的搭建；Gnedenko （1943）对极值分布的渐近原理给出了严格的证明。在此基础上，众多学者从理论分析（Jenkinson，1955；De Haan，1970，1971）与实际应用（史道济，2006；Bhunya et al，2012）两个方面开展了对极值理论的大量研究。
最具破坏力的自然灾害之一的地震，尤其是超强震，对社会经济造成的损失是巨大的，因此对特定区域作地震风险估计对于宏观决策和微观管理都是必要的。Nordquist （1945）将极值理论引入到震级预测过程；多名研究人员（Epstein，Lomnitz，1966；陈培善，林邦慧，1973；Yegulalp，Kuo，1974）对最大震级预测模型作了系统的研究；高盂潭和贾素娟（1988）对极值理论在工程抗震中的应用作了详细论述；陈虹和黄忠贤（1995）及钱小仕等（2012，2013）给出了基于统计分析的地震危险性评价指标及其计算方法。
对于广义极值模型，对模型分类和统计规律的表达起主要作用的是形状参数，其数值变化是研究者比较关心的。极值分布主要用极大似然估计和矩估计法确定参数取值，对于参数的区间估计，特别是地震重现水平的区间估计，传统的Delta法（Rao，1965）得到的对称区间，随着震级的提高，与预测结果在大于重现水平取值时有更大的不确定性的实际情况不符。为此，一些科研专家尝试用轮廓似然估计代替极大似然估计，用以确定极值分布的参数值。Murphy和Van der Vaart （2000）证明了在通常情况下轮廓似然估计与极大似然估计是等价的。Tajvidi （2003）、Gilli和Këllezi （2006）及鲁帆等（2013）对各类实际问题的风险评估，证明了轮廓似然区间估计的不对称性恰是对极值变量的不确定性随取值增加而变大的较好表达。
综上所述，在基于广义极值分布对地震危险性进行评价的过程中，本文拟用轮廓似然函数进行参数估计，以便更合理地表达强震预测的不确定性。
1. 极值分布相关理论简述
作为风险管理和安全评价的重要分析工具，极值分析主要针对随机变量中的极端变异性数据进行系统建模，进而对极端事件在未来发生的可能性给出预测评价。为了更清晰地表述基于轮廓似然估计的广义极值分布模型的构建过程，对重要理论和概念概述如下：
1.1 广义极值分布的基本理论
极值分析的基础定理（Fisher，Tippett，1928）：若X1，X2，···，Xn，···是分布函数为F（x）的独立同分布的随机变量序列，对于任意n∈Z＋，令Yn＝max{X1，X2 ，···，Xn}，如果存在常数列{αn＞0}和{βn}，以及非退化函数Λ（x），使得
$$ \underset{n\to \mathrm{\infty }}{\mathrm{lim}}\left(\frac{{Y}_{n}－{\beta }_{n}}{{\alpha }_{n}}{\text{≤}} x\right)＝{\varLambda } （ x ） ，x\in {{\bf{R}}} $$ (1) 成立，则称Λ（x）为极值分布，并称F（x）属于极值分布Λ（x）的最大值吸引场，记为F（x）∈MDA（Λ）。
吸引场的概念说明相互独立且服从相同分布函数F（x）的随机变量序列，在满足一定的条件下，序列的区组最大值近似服从广义极值（generalized extreme value，缩写为GEV）分布。对于验证条件αn和βn难于确定的问题，Fréchet （1927）证明了区组最大值分布规律是稳定的：即使随机变量分布不同，若区组最大值的近似分布存在，则其服从GEV分布，这是极值理论应用的基础。
Jenkinson （1955）将Λ（x）的三种极值分布形式统一为一个表达式，称之为广义极值分布。具有位置参数和尺度参数的广义极值分布的分布函数记为
$$ {{Gev}} （ x；\xi ，\mu ，\sigma ） ＝\exp\left[-{\left(1＋\xi \frac{x－\mu }{\sigma }\right)}^{-\frac{1}{\xi} }\right]，1＋\frac{\xi （ x－\mu ） }{\sigma } ＞ 0 \text{，} $$ (2) 式中，
$\xi $ 为形状参数，μ为位置参数，σ为尺度参数。GEV分布的类型由形状参数ξ的符号决定：当ξ＞0时，Gev（x，ξ ，μ，σ）表示Fréchet分布；当
$\xi＝0 $ 时，Gev（x，$\xi $ ，μ，σ）表示耿贝尔（Gumbel）分布；当$\xi ＜0$ 时，Gev（x，$\xi $ ，μ，σ）表示具有有限上端点的威布尔（Weibull）分布。推导过程史道济（2006）一文中有详细描述。应用极值理论进行安全评价时，常用到分位数概念，定义如下：
设随机变量X∽F（x），若xp满足条件
$$ P\left\{X{\text{≤}} {x}_{p}\right\}＝p，0 ＜ p ＜ 1 \text{，} $$ (3) 则称xp为F（x）的p分位数。
从分位数的定义可知：xp＝F−1（p），因此，由GEV分布函数（式2）可得
$$ {x}_{p}＝\mu －\frac{\sigma }{\xi } ［ 1－{ （ -\ln p ） }^{-\xi } ］ ，\xi {\text{≠}} 0 \text{，} $$ (4) 若
$\xi ＜0 $ ，令p→1可得GEV分布理论的上端点$$ {x}^{*}＝\mu －\frac{\sigma }{\xi } {\text{.}} $$ (5) 为了尽可能地满足样本之间相互独立的条件，对时间序列进行区组划分时，应设置足够长的时间间隔，以每个时间段中最大值序列作为GEV建模样本，最大限度地满足样本的渐近独立性。在此基础上，进行模型参数估计和适用性检验，最后利用模型进行安全性评价.
1.2 GEV分布参数的轮廓似然估计
GEV分布的形状参数
$ \xi $ 可以确定模型的分类及密度曲线形状，而分位数是安全评价的重要指标，本节将利用轮廓似然估计给出这两个重要参数的估计.设随机变量X∽f（x，θ），θ∈Θ，其中概率密度f（x，θ）的形式已知，θ＝（θ1，θ2，···，θk）包含k个未知参数，Θ为θ的可能取值范围。X1，X2，···，Xk是样本，x1，x2，···，xk是样本值，则（X1，X2，···，Xk）的联合概率密度在（x1，x2，···，xk）的函数值L（θ）称为样本的似然函数。
取
$\Theta_0 \subset \Theta $ ，令$ {\widehat{\theta }}_{0} $ 和$ \widehat{\theta } $ 分别为θ在Θ0和Θ中的极大似然估计值，称λ为对数似然比值，$\lambda （ x ） ＝-2\ln ［ {L （ x；{\widehat{\theta }}_{0} ） }/ {L （ x；\widehat{\theta } ） } ］ $ ，相应的λ（X）称为对数似然比统计量。史宁中（2008）的研究表明，在一定条件下，当Θ0＝{θ0}$ \mathrm{时} $ ，有$ {\widehat{\theta }}_{0} $ ＝θ0，此时有$ \lambda （ X ） \stackrel{L}{\to }{\chi }^{2} （ k ） $ ，即对数似然比近似服从自由度为k的χ2分布。这是轮廓似然估计法进行参数区间估计的理论基础。若将θ中的未知参数分成两类
$ {\theta }_{i} 和 {\overline{\theta }}_{i} $ ，则似然函数写作L（$ {\theta }_{i}，{\overline{\theta }}_{i} $ ），其中θi是研究者重点关注的分量，相应的$ {\overline{\theta }}_{i} $ 是θ的其它未知分量。则θi的轮廓似然函数定义为$L_{\rm{p}} （ {\theta }_{i} ） ＝\underset{{\overline{\theta }}_{i}}{\mathrm{max}}L （ {\theta }_{i}，{\overline{\theta }}_{i} ） $ ，即取定θi时，函数值Lp（θi）是L（$ {\theta }_{i}，{\overline{\theta }}_{i} $ ）的最大值。基于上述理论，在利用GEV分布进行安全评价过程中主要处理ξ≠0的情况，所以下文仅对ξ≠0的情况进行表述，ξ＝0的推导过程与之类似。
由GEV的概率密度公式
$$ {{gev}} （ x；\xi ，\mu ，\sigma ） ＝\frac{1}{\sigma }{{gev}} （ x；\xi ，\mu ，\sigma ） { （ 1＋\xi \frac{x－\mu }{\sigma } ） }^{-\left(1＋\frac{1}{\xi }\right)} ，1＋\frac{\xi （ x－\mu ） }{\sigma } ＞ 0 $$ (6) 可得，GEV的对数似然函数为
$$ \ln L （ \xi ，\mu ，\sigma ） ＝n\ln\frac{1}{\sigma }－\left(1＋\frac{1}{\xi }\right)\sum\limits _{i＝1}^{n}\ln\left[1＋\xi \left(\frac{{x}_{i}－\mu }{\sigma }\right)\right]－\sum\limits _{i＝1}^{n}{\left[1＋\xi \left(\frac{{x}_{i}－\mu }{\sigma }\right)\right]}^{-\frac{1}{\xi }} \text{，} $$ (7) 式中，
$ 1＋\xi \left(\dfrac{{x}_{i}-\mu }{\sigma }\right) ＞ 0，i＝\mathrm{1，2}，\cdots ，n$ 。因此GEV关于形状参数
$ \xi $ 的轮廓似然函数可表示为$L_p （ \xi ） ＝\underset{\mu ，\sigma }{\mathrm{max}}\ln L （ {\xi }_{i}，\mu ，\sigma ） $ ，即对于在可能取值范围内ξ的每个值，其函数值Lp（ξ）取在μ和σ定义范围内，使得lnL（ξ ，μ，σ）取得最大值，即GEV轮廓似然函数对应的点集为$$ \left\{ （ {\xi }_{i}，L_{\rm{p}} （ {\xi }_{i} ） ） |L_{\rm{p}} （ {\xi }_{i} ） ＝\ln L （ {\xi }_{i}，{\mu }_{i}，{\sigma }_{i} ） ＝\underset{\mu ，\sigma }{\mathrm{max}}\ln L （ {\xi }_{i}，\mu ，\sigma ） \right\} \text{，} $$ (8) 由此可得
$\xi $ 的轮廓似然估计值为$$ \hat{\xi }{＝\xi }_{i0}，\mathrm{其}\mathrm{中}{L_{\rm{p}}} （ {\xi }_{i0} ） ＝\underset{i}{{\max}}L_{\rm{p}} （ {\xi }_{i} ） {\text{.}} $$ (9) 如果
$ \hat{\xi } $ ＜0，则可得GEV分布的上端点为$$ {x}^{*}＝{\mu }_{i0}－\frac{{\sigma }_{i0}}{{\xi }_{i0}} {\text{.}} $$ (10) 因为对数似然比近似服从χ2分布，所以当
$\xi ＝\xi_i$ 成立时，λ（$\xi $ ）＝2{Lp（$ \hat{\xi } $ ）－Lp（$\xi $ ）}∽χ2（1），进而可得$\xi $ 的置信水平为1－α的置信区间为$$ I{E}_{1－\alpha } （ \xi ） ＝\left\{{\xi }_{i}|\lambda （ {\xi }_{i} ） {\text{≤}} {{\chi }}^{2}_{1-\mathrm{\alpha }} （ 1 ） \right\} {\text{.}} $$ (11) 为了求分位数的置信区间，必须在似然函数中引入参数xp 。由分位数的计算（式4）可得
$\; \mu ＝{x}_{p}＋{\sigma }/{\xi } ［ 1－{ （ -\ln p ） }^{-\xi } ］ $ ，将μ代入到对数似然表达式（式7）得lnL（ξ， xp， σ）。因此，GEV关于xp的轮廓似然函数$L_{\rm{p}} （ {x}_{p} ） ＝\underset{\xi ，\sigma }{\mathrm{max}}\ln L （ \xi ， {x}_{p}， \sigma ） $ 。即对于在可能取值范围内xp的每个值，其函数值Lp（xp）取在ξ，σ定义范围内，使得lnL（ξ， xp， σ）取得最大值，即轮廓似然函数对应的点集为$$ \left\{ （ {x}_{{p}_{i}}， L_{\rm{p}} （ {x}_{{p}_{i}} ） ） |L_{\rm{p}} （ {x}_{{p}_{i}} ） ＝\ln L （ {\xi }_{i}， {x}_{{p}_{i}}， {\sigma }_{i} ） ＝\underset{\xi ， \sigma }{\mathrm{max}}\ln L （ \xi ， {x}_{{p}_{i}}， \sigma ） \right\} ，$$ (12) 由此可得xp的轮廓似然估计值为
$$ {\hat{x}}_{p}＝{x}_{{p}_{i0}}，{L_{\rm{p}}} （ {x}_{{p}_{i0}} ） ＝\underset{i}{{\max}}L_{\rm{p}} （ {x}_{{p}_{i}} ） ，$$ (13) 类似于
$\xi $ ，xp的置信水平为1－α的置信区间为$$ I{E}_{\alpha } （ {x}_{p} ） ＝\left\{{x}_{{p}_{i}}|\lambda （ {x}_{{p}_{i}} ） {\text{≤}} {{\chi }}^{2}_{1－{\alpha }} （ 1 ） \right\} {\text{.}} $$ (14) 1.3 地震发生的重现期与重现水平
假设X1，X2 ，···，Xn为某一特定区域以半年为间隔的地震震级最大值样本，在进行余震删除工作后，鉴于时间间隔比较长，可认为样本满足渐近独立性，且服从GEV分布。设第一次出现超阈值u的时间为τ1＝min{k：Xk＞u}，令P{X＞u}＝1－Gevξ（u）＝q，则有
$$ P\left\{{\tau }_{1}＝k\right\}＝P\left\{{X}_{1}{\text{≤}} u， {X}_{2}{\text{≤}} u， \cdots ， {X}_{k－1}{\text{≤}} u， {X}_{k} ＞ u\right\}＝q{ （ 1-q ） }^{k－1}， k＝\mathrm{1， 2}， \cdots {，} $$ (15) 即第一次出现超阈值的时间服从几何分布。由分布的性质可得
$E （ {\tau }_{1} ） ＝\sum\limits_{k＝1}^{n}kq{ （ 1－q ） }^{k－1}＝ {1}/{q}$ ，也就是说第一次出现超阈值u的平均时间$T （ u ） ＝{1}/{q}＝{1}/ ［ {1－{{{Gev}}}_{\xi } （ u ） } ］ $ 。另由几何分布的性质可得：相邻两次超阈值的时间与第一次出现超阈值的时间理论上是相等的，这是危险事件重现期应具备的基本属性。依上述讨论，重现期定义为
给定时间序列X1，X2，···，Xn及阈值u，若序列中变量第一次出现超阈值u的平均时间为T（u），则称u为重现水平，相应的T（u）称为重现期。由上述分析可得重现水平为u的重现期为
$$ T （ u ） ＝\frac{1}{1-{{{Gev}}}_{\xi } （ u ） } \text{，} $$ (16) 相应地，给定重现期
$ T $ ，求重现期的反函数，可得相应的重现最大震级（重现水平）为$$ u （ T ） ＝{{{{Gev}}}_{\xi }}^{-1}\left(1－\frac{1}{T}\right) {\text{.}} $$ (17) 实际上重现期为T的重现水平u（T）就是GEV分布的
$ 1－{1}/{T} $ 分位数，即$u （ T ） ＝{x}_{1－{1}/{T}}$ 。重现期与重现水平是进行地震预报分析的两个重要指标，也是进行地震应急预案制定的重要参考因素。下文将以轮廓似然估计法为工具，对东昆仑地震带的地震危险性进行综合分析。2. 基于轮廓似然估计的东昆仑地震带地震危险性分析
2.1 地震数据样本完整性分析
依据对大陆活动地块划定的东昆仑地震带区域边界（张国民等，2005），以从国家地震科学数据共享中心获取的正式地震目录为源数据，提取（82.8°E—104.2°E，33.5°N—37.1°N）自公元前193年至2019年12月的4 385条记录作为研究样本。为了最大程度地满足地震记录样本的相互独立性，采用震级相关时空窗法（陈凌等，1998）对研究样本中的余震进行剔除，得到东昆仑地震带的震例数据3 616条。本文对地震大小的描述采用面波震级MS，如果面波震级缺失，通过公式MS＝1.13ML－1.08 （汪素云等，2010）将近震震级ML转换为MS，进行数据填充，进而得到用于地震危险性分析的完整数据样本。
东昆仑地震带的地震空间分布如图1a所示，整个地震带内断裂带密集排布，在高海拔地区大体沿西北向东南展布。地震带内发震点的分布呈现东北部频率高、震级低，中部和西南部频率低、震级高的特点。MS7.0以上强震基本上都发生在断裂带上，而在地壳隆起边缘的断裂带上，地震发生频率较高。
![]()
图 1 东昆仑地震带的地震分布规律（a） 地震空间分布；（b） M-t图Figure 1. Distribution law of earthquakes in East Kunlun seismic zone（a） Spatial distribution of earthquakes；（b） M-t diagram东昆仑地震带MS5.0以上地震从1937年开始记录（黄玮琼等，1994），MS6.0上地震从1926年开始记录，1970年中国地震台网建立后，发震情况记录则比较完整。1930—2019年的MS3.0以上震级与时间关系如图1b所示。M-t图表明东昆仑地震带MS4.0和MS3.0以上地震分别从1950年和1965开始才有相对完整的记录。为了尽可能地保留数据信息的同时满足模型分析对数据连续性的要求，本文以半年为区组间隔，提取了1950年后的震级最大值为地震危险性分析样本。
2.2 基于轮廓似然估计的东昆仑地震带GEV分布模型构建
按1.2节GEV分布参数的轮廓似然估计理论，以Matlab软件为计算平台，以遗传算法作为数值计算的主要工具，按如下步骤估计GEV分布的形状参数
$\xi $ ：首先，设定$\xi $ 的可能取值范围为 ［$\xi_l $ ，$\xi_u $ ］ ，对于任意$\xi_i $ ∈ ［$\xi $ l ，$\xi_u $ ］ ，依据式（8），搜索μi，σi，使得$L_{\rm{p}} （ {\xi }_{i} ） ＝\mathop{\mathrm{max}}\limits_{{\mu ，\sigma }}\ln L （ {\xi }_{i}，\mu ，\sigma ） $ ，并将满足条件的点（$\xi_i $ ，μi，σi，Lp（$\xi_i $ ））加入轮廓似然估计点集${{PL}}_{\xi} $ ；其次，依据式（9）确定i0，由$\xi $ 的轮廓似然估计值$\xi_{i0} $ 和相应的μi0，σi0，确定GEV分布函数GEV（x；$\xi_{i0} $ ，μi0，σi0）；最后，依据式（11）确定$ \xi $ 的置信水平为1－α的置信区间。根据上述步骤得到的轮廓似然估计点集PLξ如图2所示。图2表明轮廓似然函数与形状参数
$\xi $ 存在类似抛物线的关系，在ξ＝−0.204 0时轮廓似然函数取得最大值，此时μ＝0.847 5，σ＝4.834 5。进而可得，GEV的概率密度函数中参数（$\xi $ ，μ，σ）的估计值为（−0.204 0，0.847 5，4.834 5）。![]()
图 2 形状参数与轮廓对数似然函数之间的关系Figure 2. Relationship between shape parameters and profile log likelihood function基于轮廓似然估计构建的GEV分布模型对样本数据分析的适用性检验如图3所示。直方图轮廓趋势线与概率密度曲线基本吻合；P-P图的散点在56°线附近小幅波动，表明样本与理论分布的分位数特征基本是一致的；且
$ \hat{\xi } $ ＜0表示震级的重现水平有理论上限。上述检验表明，本文构建的GEV分布模型适用于对昆仑山地震带作地震危险性分析。![]()
图 3 GEV分布模型适应性检验图（a） 密度曲线与直方图；（b） P-P 检验Figure 3. Adaptability test of GEV distribution model（a） Density curves and histograms；（b） P-P test由于
$ \xi $ 的轮廓置信上限小于0，则依据式（10）可求震级的理论上限.同时，在区组震级X∽GEV （x；ξi0，μi0，σi0）的条件下，利用$E （ X ） ＝{{\mu }_{i0}－{{\sigma }_{i0}}{ ［ 1－ \mathrm{\Gamma } （ 1－{\xi }_{i0} ） ］ }/{{\xi }_{i0}}}，{\xi }_{i0} ＜ 1$ ，可求得最大震级理论平均值，其中Γ（x）是伽马函数。为了比较参数估计方法的差异，表1分别列出了轮廓似然法和极大似然法（钱小仕等，2012）的估计结果，在进行模型主要参数估计时，两种估计方法的效果基本相同。
表 1 轮廓似然估计与极大似然估计GEV的分布结果对比Table 1. Comparation of the profile likelihood estimation and maximum likelihood estimationof GEV distribution模型分析项目 轮廓似然估计 极大似然估计 参数估计 （−0.204 0，0.847 5，4.834 5） （−0.204 4，0.847 8，4.834 8） 形状参数置信区间 ［ −0.259 0，−0.134 0 ］ ［ −0.268 8，−0.140 1 ］ 震级理论上限 MS8.989 1 MS8.981 9 地震带最大震级均值 MS5.178 9 MS5.179 0 通过上述分析可知，东昆仑地震带每半年的最大震级平均约为MS5.2，理论上的最大震级约为MS9.0，属于巨震级别，说明这一区域的地质条件很不稳定，属于大地震发生危险性极高的地区。
2.3 东昆仑地震带重现水平及置信区间的轮廓似然估计
上文估计的半年最大震级理论均值和上限，是对地区发震情况的简要概述，对应急预案制定参考意义不是很大。重现期和重现水平是评估地震危险性的两个核心指标，在估计了GEV分布重要参数后，给定重现期T，按式（17）可以确定理论重现水平u（T）。为了满足决策需要，有时需要求出重现水平的置信区间。利用信息矩阵获得的置信区间（Rao，1965）是关于置信水平对称的，但实际上随着预测震级的提高，在高于置信水平的时候，震级选择会更加离散，也就是说置信区间关于置信水平对称是不合理的。为此，下文将利用轮廓似然法来估计重现水平的置信区间。
在给定重现期的情况下，确定重现水平及置信区间可按如下步骤进行：首先，对于给定重现期T，按式（17）可以确定理论重现水平u（T），根据u（T）值设定一个相对保守的重现水平取值范围 ［ XPl，XPu ］ 。对于任意
$x_{p_i}$ ∈ ［ XPl，XPu ］ ，依据式（12），搜索ξi ，σi，使得$L_{\rm{p}} （ {x}_{{p}_{i}} ） ＝ \mathop{\mathrm{max}}\limits_{{\xi ，\sigma }}\ln L （ \xi ，{x}_{{p}_{i}}，\sigma ） $ ，并将满足条件的点（$\xi_i $ ，$x_{p_i} $ ，σi，Lp（$x_{p_i} $ ））加入轮廓估计点集$PL_{x_p} $ ；其次，依据式（13）确定i0，得到xp的轮廓似然估计值$x_{p_{i0}} $ ；最后，依据式（14）确定xp的置信水平为1－α的轮廓置信区间。根据上述步骤得到的重现期分别为20年、50年、100年、500年的重现水平和置信区间如图4所示。
![]()
图 4 重现水平及置信区间的轮廓似然估计（a） 20年重现期；（b） 50年重现期；（c） 100年重现期；（d） 500年重现期Figure 4. The reappearance level and confidence interval of the profile likelihood estimation（a） 20-year return period；（b） 50-year return period；（c）100-year return period；（d） 500-year return period用Delta法和轮廓似然估计得到的估计结果对比分析详见表2。估计结果从数值计算的角度佐证了对于重现水平的估计，轮廓似然法与极大似然法是等效的（Murphy，Vander Vaart，2000）.
表 2 极大似然估计与轮廓似然估计的重现水平对比Table 2. Comparation of the recurrence level between profile likelihood estimation and maximum likelihood estimation重现期
/年极大似然估计
重现水平极大似然估计95%
置信区间轮廓似然估计
重现水平轮廓似然估计95%
置信区间轮廓估计重现水平
两侧区间长度比1 MS5.13 ［ 4.97，5.29 ］ MS5.13 ［ 4.98，5.29 ］ 1.07 5 MS6.36 ［ 6.15，6.57 ］ MS6.36 ［ 6.17，6.59 ］ 1.21 10 MS6.72 ［ 6.48，6.96 ］ MS6.72 ［ 6.51，7.00 ］ 1.33 20 MS7.03 ［ 6.75，7.30 ］ MS7.03 ［ 6.79，7.37 ］ 1.42 50 MS7.36 ［ 7.03，7.69 ］ MS7.36 ［ 7.10，7.80 ］ 1.69 100 MS7.58 ［ 7.20，7.95 ］ MS7.58 ［ 7.29，8.10 ］ 1.79 500 MS7.97 ［ 7.48，8.46 ］ MS7.97 ［ 7.63，8.68 ］ 2.09 表2的对比结果说明，对于重现期10年以下的重现水平置信区间的估计，轮廓似然法与极大似然法的估计结果基本相同，即针对短期地震危险性评估，两种方法是等效的。但在进行中长期地震危险性分析时，轮廓似然估计得到的估计区间整体向右偏移，即置信下限和上限比相应的极大似然估计要高，而且重现水平右侧的宽度随着重现期的提高，与左侧区间宽度之比越来越大，说明随着重现期的变长，重现水平随之提高，而且发生超重现水平的震级取值更分散，是对震级区间更为保守的预测，也是与实际情况相吻合的。两种方法估计结果的直观差异如图5所示。
![]()
图 5 重现水平的轮廓似然估计与极大似然估计对比Figure 5. Comparation of the reproduction level between profile likelihood estimation and maximum likelihood estimation本文所用源数据是东昆仑地震带截止到2019年12月的数据，为了对比两种方法的分析效果，首先截取1920年后100年的震例数据，其中大于极大似然估计下限（MS7.20）和轮廓估计下限（MS7.29）的样本数都是5个，分别为1924年和1973年的MS7.3、1937年和1997年的MS7.5以及2001年的MS8.1，其中MS8.1地震已经超过极大似然估计上限MS7.95，但在轮廓似然估计上限之内。如果估计精度为MS0.1，则大于轮廓似然估计下限的样本为3个，从数量来说，优于极大似然估计.如果截取1420年后的500年的地震数据，大于极大似然估计下限（MS7.48）的样本数为3个，分别是1937年和1997年的MS7.5，以及2001年的MS8.1，而超过轮廓似然估计下限（MS7.63）的样本只有1个MS8.1震例，重现数量上优于极大似然估计。以上实例分析表明，轮廓似然估计作为震级重现水平区间估计的理论工具，比以信息矩阵为基础的Delta法更加合理和适用。
3. 讨论与结论
本文对轮廓似然估计法在广义极值（GEV）分布模型的参数估计中的应用从理论分析到数值计算进行了详细地阐述，并将构建的极值分布模型应用于东昆仑地震带的地震预报。在对参数的点估计过程中，轮廓似然估计与极大似然估计效果相当；在进行重现水平置信区间估计过程中，如果重现期比较短，两种方法估计效果几乎无差异，但在进行中长期地震预报时，轮廓似然估计法得到的重现水平置信区间对于不确定性的表达更加充分，较基于信息矩阵的Delta法得到的对称置信区间更为合理。
基于轮廓似然估计的GEV分布模型能够对地震危险性作出相对比较客观的评价，但该模型所用数据为区组最大值信息，在模型构建过程中对观测信息利用不够充分，使得由模型得到的预测结果与实际情况存在偏差，作者后续将主要从历史信息挖掘和模型调整与优化两个方向入手，以构建更加有效的地震预报预测模型。
-
![]()
图 1 M≥6.0地震的理论计算等震线与实际等震线的对比图
Figure 1. Comparison between theoretical calculation isoseismal lines and actual isoseismal lines for M≥6.0 earthquakes
（a） M6.0—6.5；（b） M6.6—7.0；（c） M7.1—7.5；（d） M7.8；（e） M8.0
![]()
图 1 M≥6.0地震的理论计算等震线与实际等震线的对比图
Figure 1. Comparison between theoretical calculation isoseismal lines and actual isoseismal lines for M≥6.0 earthquakes
（a） M6.0—6.5；（b） M6.6—7.0；（c） M7.1—7.5；（d） M7.8；（e） M8.0
![]()
图 2 1976年唐山M7.8地震主震后2小时（a）、6小时（b）、12小时（c）和24小时（d）内余震频度分布图
Figure 2. Distribution of aftershock frequency within 2 hours （a），6 hours （b），12 hours （c） and 24 hours （d） after Tangshan M7.8 earthquake in 1976
![]()
图 3 理论等震线极震区修正前（a）后（b）对比图
Figure 3. Comparison of theoretical isoseismal areas before （a） and after （b） correction
![]()
图 5 1966年宁晋东南M7.6地震（a）和1976年唐山M7.8地震（b）的地震影响场对比图
Figure 5. Comparison of influence fields between the 1966 Ningjin southeast M7.2 earthquake （a） and the 1976 Tangshan M7.8 earthquake （b）
表 1 研究区内选取的震例
Table 1 The selected earthquake cases in the studied area
编号 发震时间
年-月-日震中位置 精度类别 震中烈度 震源深度/km M 震源机制类型 震中参考位置 北纬/° 东经/° 1 1618-11-16 39.80 114.50 2 Ⅷ $6\tfrac{1}{2}$ 正断 河北蔚县附近 2 1624-04-17 39.50 118.80 3 Ⅷ $6\tfrac{1}{2}$ 走滑 河北滦县 3 1628-10-07 40.70 114.20 2 Ⅷ $6\tfrac{1}{2}$ 走滑兼正断 河北怀安西洋河堡 4 1665-04-16 39.90 116.60 2 Ⅷ $6\tfrac{1}{2}$ 逆冲-走滑 北京通县西 5 1679-09-02 40.00 117.00 2 Ⅺ 8.0 走滑 河北三河平谷 6 1720-07-12 40.40 115.50 2 Ⅸ $6\tfrac{3}{4}$ 走滑 河北沙城 7 1730-09-30 40.00 116.20 2 Ⅷ $6\tfrac{1}{2}$ 走滑 北京西北郊 8 1830-06-12 36.40 114.30 2 Ⅹ $7\tfrac{1}{2}$ 走滑 河北磁县 9 1882-12-02 38.10 115.50 2 Ⅷ 6.0 正断 河北深县 10 1945-09-23 39.70 118.70 Ⅷ $6\tfrac{1}{4}$ 走滑 河北滦县 11 1966-03-08 37.35 114.92 Ⅸ+ 6.8 走滑 河北隆尧东 12 1966-03-22 37.50 115.10 2 Ⅹ 9 7.2 走滑 河北宁晋东南 13 1966-03-26 37.68 115.27 1 Ⅶ+ 15 6.2 走滑 河北束鹿南 14 1967-03-27 38.50 116.50 Ⅶ 6.3 走滑 河北河间、大城 15 1976-07-28 39.60 118.20 Ⅺ 22 7.8 走滑 河北唐山 16 1976-07-28 39.90 118.70 Ⅸ 22 7.1 走滑 河北滦县 17 1976-11-15 39.33 117.50 Ⅷ 17 6.9 逆断 天津宁河西 18 1977-05-12 39.20 117.70 1 Ⅶ 19 6.2 走滑 天津汉沽附近 19 1998-01-10 41.12 114.43 Ⅷ 10 6.2 逆断 河北张北 
下载: 导出CSV
表 2 研究区选用的地震烈度衰减系数
Table 2 Seismic intensity attenuation coefficients used in the studied area
发震方式 C1 C2 C3 a0 σ 走滑型 长轴 5.291 0 1.438 0 −4.305 4 25 0.622 4 短轴 3.148 8 1.338 7 −3.272 4 14 0.649 2 全震例 长轴 5.861 9 1.390 2 −4.451 5 25 0.586 2 短轴 2.954 9 1.349 4 −3.106 4 10 0.615 3 注：C1，C2，C3为回归系数，a0表示距离饱和因子，σ表示衰减关系的误差。
下载: 导出CSV
表 3 华北地区地震动衰减关系系数
Table 3 Attenuation coefficients of ground motion in North China
系数与
方差椭圆长轴 椭圆短轴 M＜6.5 M≥6.5 M＜6.5 M≥6.5 A 2.024 3.565 1.204 2.789 B 0.673 0.435 0.664 0.420 C 2.329 2.329 2.016 2.016 D 2.088 2.088 0.944 0.944 E 0.399 0.399 0.447 0.447 σ 0.245 0.245 0.245 0.245
下载: 导出CSV
表 4 适用于中国场地分类的场地放大系数
Table 4 Site magnification factors for site classification in China
场地
类型不同基岩地震动加速度PGA （cm·s−2）下的放大系数 PGA≤100 PGA＝200 PGA＝300 PGA＝400 PGA≥500 Ⅰ 1.0 1.0 1.0 1.0 1.0 Ⅱ 1.4 1.3 1.2 1.1 1.0 Ⅲ 2.1 1.6 1.2 1.0 1.0 Ⅳ 2.5 1.7 1.2 0.9 0.9
下载: 导出CSV
表 5 烈度与水平向峰值加速度的对应关系
Table 5 Corresponding relationship between intensity and horizontal peak acceleration
烈度 水平向地面峰值加速度/（cm·s−2） Ⅵ 45—89 Ⅶ 90—177 Ⅷ 178—353 Ⅸ 354—707
下载: 导出CSV
-
白仙富,戴雨芡,赵恒. 2014. 地震影响场应急评估方法研究[J]. 自然灾害学报,23(4):91–102. Bai X F,Dai Y Q,Zhao H. 2014. Study on the emergency evaluation of earthquake influence field[J]. Journal of Natural Disasters,23(4):91–102 (in Chinese).
陈鲲,俞言祥,高孟潭. 2010. 考虑场地效应的ShakeMap系统研究[J]. 中国地震,26(1):92–102. doi: 10.3969/j.issn.1001-4683.2010.01.009 Chen K,Yu Y X,Gao M T. 2010. Research on ShakeMap system in terms of the site effect[J]. Earthquake Research in China,26(1):92–102 (in Chinese).
国家地震局震害防御司. 1995. 中国历史强震目录[M]. 北京: 地震出版社: 158–414. Department of Earthquake Disaster Prevention, State Seismological Administration. 1995. Catalogue of Historical Strong Earthquakes in China[M]. Beijing: Seismological Press: 158–414 （in Chinese）.
郝敏,谢礼立. 2006. 集集地震等震线和PGA、PGV等值线关系的研究[J]. 地震工程与工程振动,26(1):18–21. doi: 10.3969/j.issn.1000-1301.2006.01.003 Hao M,Xie L L. 2006. Study on the relationship between isoseismal and isolines of PGA and PGV for Chi-Chi earthquake[J]. Earthquake Engineering and Engineering Vibration,26(1):18–21 (in Chinese).
胡聿贤,张敏政. 1984. 缺乏强震观测资料地区地震动参数的估算方法[J]. 地震工程与工程振动,4(1):1–11. Hu Y X,Zhang M Z. 1984. A method of predicting,ground motion parameters for regions with poor ground motion data[J]. Earthquake Engineering and Engineering Vibration,4(1):1–11 (in Chinese).
李志强,袁一凡,李晓丽,何萍. 2008. 对汶川地震宏观震中和极震区的认识[J]. 地震地质,30(3):768–777. doi: 10.3969/j.issn.0253-4967.2008.03.015 Li Z Q,Yuan Y F,Li X L,He P. 2008. Some insights into the macro-epicenter and meizo seismal region of Wenchuan earthquake[J]. Seismology and Geology,30(3):768–777 (in Chinese).
龙锋,闻学泽,徐锡伟. 2006. 华北地区地震活断层的震级-破裂长度、破裂面积的经验关系[J]. 地震地质,28(4):511–535. doi: 10.3969/j.issn.0253-4967.2006.04.001 Long F,Wen X Z,Xu X W. 2006. Empirical relationships between magnitude and rupture length,and rupture area,for seismogenic active faults in North China[J]. Seismology and Geology,28(4):511–535 (in Chinese).
吕红山,赵凤新. 2007. 适用于中国场地分类的地震动反应谱放大系数[J]. 地震学报,29(1):67–76. doi: 10.3321/j.issn:0253-3782.2007.01.008 Lü H S, Zhao F X. 2007. Site coefficients suitable to China site category. Acta Seismologica Sinica, 29（1）: 67–76 （in Chinese）. 沈正康,万永革,甘卫军,李铁明,曾跃华. 2004. 华北地区700年来地壳应力场演化与地震的关系研究[J]. 中国地震,20(3):211–228. doi: 10.3969/j.issn.1001-4683.2004.03.001 Shen Z K,Wan Y G,Gan W J,Li T M,Zeng Y H. 2004. Crustal stress evolution of the last 700 years in North China and earthquake occurrence[J]. Earthquake Research in China,20(3):211–228 (in Chinese).
石建梁,闫庆民,葛秋莹. 2011. 用椭圆衰减关系模型计算任意场点烈度及地震动参数的数值方法[J]. 内陆地震,25(1):21–28. doi: 10.3969/j.issn.1001-8956.2011.01.004 Shi J L,Yan Q M,Ge Q Y. 2011. An algorithm for arbitrary engineering site earthquake intensity or motion parameter using ellipsoid attenuation model[J]. Inland Earthquake,25(1):21–28 (in Chinese).
孙艳萍,陈文凯,周中红,张苏平,何少林,李少华. 2014. 甘肃地区地震烈度影响场评价及其与地形因子的关系[J]. 地震工程学报,36(3):697–704. doi: 10.3969/j.issn.1000-0844.2014.03.0697 Sun Y P,Chen W K,Zhuo Z H,Zhang S P,He S L,Li S H. 2014. Evaluation of seismic intensity influence field in Gansu Province and its relationship with terrain factors[J]. China Earthquake Engineering Journal,36(3):697–704 (in Chinese).
田家勇,兰晓雯,谢周敏,陆鸣,时振梁. 2010. 地形对地震烈度衰减的影响[J]. 震灾防御技术,5(3):281–287. doi: 10.3969/j.issn.1673-5722.2010.03.002 Tian J Y,Lan X W,Xie Z M,Lu M,Shi Z L. 2010. Influence of topography on seismic intensity attenuation[J]. Technology for Earthquake Disaster Prevention,5(3):281–287 (in Chinese).
王想,王亚茹,宫猛,郭蕾,王晓山. 2017. 华北第5活动幕平静时段背景下平静异常综合判定[J]. 地震地磁观测与研究,38(1):21–27. doi: 10.3969/j.issn.1003-3246.2017.01.004 Wang X,Wang Y R,Gong M,Guo L,Wang X S. 2017. Comprehensive judgment of calm anomaly under the background of quiet time of the fifth activity in North China[J]. Seismological and Geomagnetic Observation and Research,38(1):21–27 (in Chinese).
巫建伟,陈崇成,吴小竹,林剑峰,黄昭,张锦福,郑师春,张颖. 2013. 基于GeoKSCloud的地震影响场分析云服务研究—以福建省为例[J]. 地球信息科学学报,15(5):695–704. Wu J W,Chen C C,Wu X Z,Lin J F,Huang Z,Zhang J F,Zheng S C,Zhang Y. 2013. Cloud service for seismic influence field analysis based on GeoKSCloud:A case study in Fujian Province[J]. Journal of Geo-Information Science,15(5):695–704 (in Chinese). doi: 10.3724/SP.J.1047.2013.00695
许卫晓. 2011. 烈度分布快速评估方法研究[D]. 哈尔滨: 中国地震局工程力学研究所: 21–22. Xu W X. 2011. The Study of Methods for Rapid Assessment of Seismic Intensity Distribution[D]. Harbin: Institute of Engineering Mechanics, China Earthquake Administration: 21–22 （in Chinese）.
张彦琪,范柱国,陈坤华,崔建文,李世成,冉华. 2013. Shakemap场地校正方法及其在云南地震动强度和烈度速报中的应用[J]. 地震研究,36(1):108–115. doi: 10.3969/j.issn.1000-0666.2013.01.017 Zhang Y Q,Fan Z G,Chen K H,Cui J W,Li S C,Ran H. 2013. Shakemap site correction method and its application in the rapid prediction of ground motion intensity and seismic intensity in Yunnan[J]. Journal of Seismological Research,36(1):108–115 (in Chinese).
张彦琪. 2014. 基于强震动台网的云南烈度速报场地影响校正研究[D]. 昆明: 昆明理工大学: 12–13. Zhang Y Q. 2014. The Study on Site Conditions and Amplification of Seismic Intensity Based on Ground Motion Intensity Network in Yunnan[D]. Kunming: Kunming University of Science and Technology: 12–13 （in Chinese）.
郑韵. 2015. 震源机制和余震序列在地震应急烈度快速判定中的应用研究[D]. 北京: 中国地震局地震预测研究所: 1–2. Zheng Y. 2015. The Application of Focal Mechanism and Aftershock Sequence in Fast Judgement of Earthquake Emergency Intensity[D]. Beijing: Institute of Earthquake Science, China Earthquake Administration: 1–2 （in Chinese）.
郑韵,姜立新,杨天青,刘杰. 2016. 基于震源机制解的分区地震烈度衰减关系研究[J]. 震灾防御技术,11(2):349–359. doi: 10.11899/zzfy20160218 Zheng Y,Jiang L X,Yang T Q,Liu J. 2016. Study on seismic intensity attenuation relationship with regions via focal mechanism solutions[J]. Technology for Earthquake Disaster Prevention,11(2):349–359 (in Chinese).
中国国家标准化管理委员会, 中华人民共和国国家质量监督检验检疫总局. 2008. 中国地震烈度表（GB/T17742—2008）[S]. 北京: 中国标准出版社: 8. China standardization administration, General administration of quality supervision, inspection and quarantine of the People’s Republic of China. 2008. The Chinese Seismic Intensity Scale（GB/T 17742—2008）[S]. Beijing: China Standard Press: 8.
中华人民共和国住房和城乡建设部, 中华人民共和国国家质量监督检疫总局. 2010. 建筑抗震设计规范（GB50011—2010）[S]. 北京: 中国建筑工业出版社: 5 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. 2010. Code for Seismic Design of Buildings （GB50011—2010）[S]. Beijing: China Building Industry Press: 5.
-
期刊类型引用(1)
1. 王继鑫，荣棉水，符力耘，傅磊. 用微动台阵记录联合反演场地浅层速度结构——以唐山响嘡台3~#场地为例. 地震地质. 2020(06): 1335-1353 . 
百度学术
其他类型引用(0)
计量
- 文章访问数: 368
- HTML全文浏览量: 230
- PDF下载量: 66
- 被引次数: 1
