引言
地震危险性分析是工程地震工作的一个重要组成部分，主要目的是为工程的抗震设计提供依据．目前国际上地震危险性评价有两种方法：确定性法和概率法．其中概率法应用最为广泛．概率地震危险性分析方法最早由Ｃｏｒｎｅｌｌ（１９６８）提出，ＭｃＧｕｉｒｅ（１９７６）根据这一方法编制了Ｆｏｒｔｒａｎ程序，因此该法又称Ｃｏｒｎｅｌｌ－ＭｃＧｕｉｒｅ法．目前我国的概率地震危险性分析采用两级划分的原则，即先划分大的地震带以求取地震活动性参数，然后在地震带内划分潜在震源区，根据空间分布函数将地震带的地震年发生率分配到各个潜在震源区中去，以计算地震危险性．由于对地震构造认识的局限性，Ｆｒａｎｋｅｌ（１９９５）及Ｆｒａｎｋｅｌ等（２０００）在美国中东部地震区划中使用了空间光滑地震活动性的方法，其最显著的特点为不根据地震构造来划分潜在震源区，而直接使用根据地震目录进行空间光滑后的点源来进行地震危险性计算（Ｃａｏｅｔａｌ，１９９６）．由于该法简单易行，在世界各国的地震危险性分析中得到了广泛应用（Ｌａｐａｊｎｅｅｔａｌ，１９９７，２００３；ＰｅｌáｅｚＭｏｎｔｉｌｌａａ，Ｃａｓａｄｏ，２００２；ＰｅｌáｅｚＭｏｎｔｉｌｌａａｅｔａｌ，２００３；Ｈａｇｏｓｅｔａｌ，２００６；杨勇等，２００８）．地震学家们在使用这一方法时往往都是把整个研究区域作为统计区域来计算地震活动性参数，并且直接根据空间格点统计得到的地震频数来计算地震年发生率并进行地震危险性计算．显然，由于地震活动在空间上具有差异性，地震活动性参数（如地震发生率、ｂ值、震级上限等）在整个研究区域内并不是相同的，因此应当根据地震活动特征的不同，划分不同的地震构造单元（地震带、地震区），统计每个单元的地震活动性参数．本文将根据我国划分的地震带，以地震带为统计区域计算地震活动性参数，然后根据Ｆｒａｎｋｅｌ（１９９５）的空间光滑法将地震事件分配到空间格点中去，以每个格点中所占的地震频度作为权重来分配地震年发生率，最后计算地震危险性．
１研究区域及资料
研究区域为华北地震区（图１）．选择华北地震区的６个地震带（区），以《中国地震动参数区划图》宣贯教材（胡聿贤，高孟潭，２００１）中“综合潜在震源区方案”的地震活动性参数作为所选６个地震带的地震活动性参数（表１）．以华北地区的仪器地震记录和历史地震记录作为计算目录，使用基于震级的余震删除法（Ｏｇａｔａｅｔａｌ，１９９５）删除了地震目录中的余震和前震．根据Ｆｒａｎｋｅｌ等（１９９６）地震目录完整性分析方法，对华北地区仪器记录的地震目录进行完整性分析（图２ａ，ｂ），即从累积Ｇ－Ｒ关系曲线和［Ｍｃ到Ｍｃ＋ｄｍ］（Ｍｃ为起始震级，ｄｍ为震级区间，一般取０．５）的时间累积频度曲线判断目录完整性．从图２ａ，ｂ中可以看出，华北地区１９７０年以来震级在Ｍ３．０—３．５的地震记录中，Ｇ－Ｒ关系曲线和时间累积频度曲线线性趋势良好，因此可以认为该时间段华北地区Ｍ≥３．０的地震记录是基本完整的．黄玮琼等（１９９４）研究表明，我国华北地震区历史破坏性地震记录自公元１５００年以来是基本完整的．图２ｄ中历史地震时间累积曲线在公元１７００—１８００年间斜率变小，这并非是由地震缺失造成的，而是该时间段为地震活动平静期（马宗晋，１９７５）．华北地区历史地震Ｇ－Ｒ关系曲线线性趋势良好（图２ｃ），因此可以认为华北地区自１５００年以来Ｍ≥４．７地震是基本完整的．
２地震活动率
计算本文中建立地震活动性模型的方法为Ｆｒａｎｋｅｌ（１９９５）提出的．即首先将研究区域划分成一定大小的网格，统计并计算网格中的地震活动率，然后根据光滑函数将网格中的值光滑到空间其它格点中去，从而形成地震活动性模型．在以往的研究中，人们往往是根据实际地震频度统计地震活动率，这样就会导致Ｍ４．０地震与Ｍ５．０地震对地震活动率的贡献是一样的（Ｆｒａｎｋｅｌｅｔａｌ，１９９６）．为了避免这种情况，人们往往采用多个活动性模型加权的方法来计算地震活动性模型．Ｌａｐａｊｎｅ等（１９９７，２００３）采用基于地震释放能量的方法来计算地震活动率．然而这一方法，会过度高估高震级所在区域的地震活动率，低估低震级区域的地震活动率．为了达到一个平衡，在本文中，我们将使用基于地震应变的地震活动率计算方法．
３模型
本研究中，根据华北地区１９７０年以来Ｍ≥３．０和１５００年以来Ｍ≥４．７的地震记录，考虑不同震级地震对地震危险性影响的差异，我们建立了４个地震活动性模型作为空间分布函数．将研究区域划分为０．１°×０．１°的网格，统计每个网格中大于等于起算震级ｍ０的地震频数ｎｉ，然后通过光滑函数将网格中的地震光滑到其它空间格点中，从而得到光滑后每个网格中的地震活动率ｎ珘ｉ．
３．１模型１使用１９７０年以来Ｍ≥３．０的地震记录，采用Ｆｒａｎｋｅｌ（１９９５）高斯空间光滑法式中，ｎ珘ｉ为第ｉ个网格进行光滑后的地震事件数目；Δｉｊ为第ｉ个网格与第ｊ个网格的距离；ｃ为相关距离．一般取ｉ网格周围３ｃ距离内网格的数据来计算ｎ珘ｉ．在高斯光滑中，相关距离ｃ的选择非常重要．在先前的研究中往往是根据地震的定位误差来确定ｃ值（Ｌａｐａｊｎｅｅｔａｌ，１９９７，２００３；Ｆｒａｎｋｅｌｅｔａｌ，１９９６，２００２；Ｐｅｔｅｒｓｅｎｅｔａｌ，２００８），ｃ值的确定具有很大的主观性和不确定性（Ｆｒａｎｋｅｌｅｔａｌ，１９９６）．特别是在地震目录较多的地区，ｃ的选择具有空间统计意义，单纯的使用地震定位误差可能是不够充分的．文中我们将在考虑地震误差的前提下，使用空间核密度估计中的交叉验证法来求取相关距离ｃ（徐伟进，高孟潭，２０１２）．图３ａ为模型１的最优相关距离ｃ值的交叉验证曲线，纵坐标为交叉验证系数（为平均均方误差），横坐标为距离，使交叉验证系数取最小值的距离为最优带宽．相关距离ｃ在２１ｋｍ处时交叉验证系数达到最小值，因此本文中模型１的最优相关距离ｃ＝２１ｋｍ．
３．２模型２模型２为华北地区１５００年以来Ｍ≥４．７地震，同样采用高斯空间光滑．图３ｂ为交叉验证曲线，相关距离ｃ在４３ｋｍ处时交叉验证系数达到最小值，则模型２的最优相关距离为ｃ＝４３ｋｍ．
３．３模型３模型３使用１５００年来Ｍ≥６．０的地震，由于６级以上的地震能够产生地表破裂，故模型３使用Ｌａｐａｊｎｅ等（２００３）提出的基于断层破裂方向的点椭圆高斯光滑，平滑公式为，高孟潭，２００１）中的“综合潜在震源区方案”来确定．用震级－破裂长度经验公式ｌｇＬ＝ａ＋ｂｍｕ（Ｗｅｌｌｓ，Ｃｏｐｐｅｒｓｍｉｔｈ，１９９４）求取地震所产生的断层破裂长度Ｌ，则椭圆长轴为σ＝κＬ，κ≥１；短轴为τ＝ωＬ，ω＜１．这一方法明确了光滑中相关距离、平滑半径选取的物理意义，并将断层因素引进光滑函数中，使光滑后的地震空间分布更趋合理，物理意义也更加明确．
３．４模型４为了强调Ｍ７．０以上地震在地震危险性计算中的重要性，我们将模型４的地震选取为华北地区记录到的所有Ｍ７．０以上的地震．由于Ｍ７．０以上地震地表破裂广，用模型３中的点椭圆光滑方法已不能满足，我们将用断层模型来计算模型４的地震活动率．模型４的算法与模型３相同，只是在模型３中将地震进行点椭圆光滑．在模型４中我们将地震看成一条断层，根据《中国地震动参数区划图》宣贯教材（胡聿贤，高孟潭，２００１）中的“综合潜在震源区方案”来确定断层方向，使用Ｗｅｌｌｓ和Ｃｏｐｐｅｒｓｍｉｔｈ（１９９４）的经验公式计算断层长度．对断层上的每个点进行与模型３相同的光滑（图４）．图４中σ与τ分别为椭圆的长短轴．用震级－破裂宽度公式ｌｇＷ＝ａ＋ｂｍｕ（Ｗｅｌｌｓ，Ｃｏｐｐｅｒ－ｓｍｉｔｈ，１９９４），求取地震破裂宽度Ｗ值，则σ＝κＷ，τ＝ωＷ，κ＞ω＞０．事实上，对于不同破裂机制的地震，其破裂宽度是具有一定差异的．然而，文中所使用的空间光滑方法是具有空间统计意义的，因此不同类型的破裂宽度的差异对最终结果的影响并不明显．此外，具体确定每次地震产生的破裂长度和宽度是一项非常复杂的工作，超出了本文的研究范围．在实际地震区划或地震危险性分析应用中，应考虑不同类型地震破裂范围的差异．关于上述模型中地震活动率的计算，模型１和模型２是直接根据每个网格中统计到的地震频度来计算，而模型３和模型４的地震活动率则是根据文中第２节描述的方法来计算的．
４地震危险性计算
地震年发生率是计算地震危险性非常重要的参数．本文中我们将华北地区各地震带的地震年发生率（ν４）根据前面描述的地震活动性模型分配到空间各格点中去，用于地震危险性计算．对于某一场点，计算场点处地震动参数值ｕ超过给定地震动参数值ｕ０的年发生率为（Ｒｅｉｔｅｒ，１９９０）式中，ｍ０为起算震级；ｎｉ（ｍ０）为第ｉ个潜源区ｍ≥ｍ０的地震年发生率，ｐｉ（ｍ）为震级的概率密度函数，ｐｉ（ｒ）为场点到潜源之间距离的概率密度函数；Ｐ（ｕ＞ｕ０｜ｍ，ｒ）为在距离场点ｒ处，震级为ｍ的地震产生的地震动值ｕ超过给定值ｕ０的概率．Ｐ（ｕ＞ｕ０｜ｍ，ｒ）中的ｕ值是根据衰减关系计算得到的．本文中使用的衰减关系为《中国地震动参数区划图》宣贯教材（胡聿贤，高孟潭，２００１）中使用的中国东西部衰减关系模型（汪素云等，２０００）．由于本研究直接使用空间格点计算地震危险性，则上式可简写为
５结果分析
图５为采用４个模型计算得到的地震活动性模型（图５ａ－－ｄ）和根据相应的地震活动性模型计算得到的５０年超越概率１０％的峰值加速度（ＰＧＡ）分布（图５ａ′－－ｄ′）．地震活动性模型中值为每个格点（０．１°×０．１°）内Ｍ４．０—４．１的地震活动率．图５ａ及图５ａ′为使用模型１计算得到的地震活动性模型和ＰＧＡ分布．可以看出，地震发生率较高的地区为环鄂尔多斯地块区的汾渭地震带和河套—银川地震带．在有仪器记录以来发生过强破坏性地震的邢台、唐山、营口地区地震发生率也较高．显然根据模型１得到的地震活动性模型主要反映了各地区现今有仪器记录以来的地震活动性水平．图５ａ′为与图５ａ对应的５０年超越概率１０％水平下的ＰＧＡ分布．可以看到，在鄂尔多斯地块区东缘的汾渭地震带及邢台、唐山、营口等地区ＰＧＡ可达２００ｃｍ／ｓ２以上，与现今地震危险性分析结果基本一致．然而在历史上曾发生过特大破坏性地震的汾渭地震带南端及郯庐地震带中南部，由于现今地震活动较弱，计算得到的ＰＧＡ值偏小．图５ｂ及图５ｂ′为使用模型２计算得到的地震活动性模型和ＰＧＡ分布．由于模型２使用的是公元１５００年以来Ｍ≥４．７的地震目录，采用的是圆形高斯光滑．因此根据模型２得到的华北地区地震活动性模型和ＰＧＡ分布，主要反映的是华北地区中强地震的活动水平和地震危险性水平．为了突出６级和７级以上地震的构造特征，我们采用模型３和模型４计算了华北地区的地震活动性模型和ＰＧＡ分布（图５ｃ，ｄ，ｃ′，ｄ′）．地震活动性模型和ＰＧＡ分布能够较好地反映华北地区的地震构造特征，并适当突出了具有地震构造地区的地震危险性水平．然而在地震构造区，构造往往较为复杂，构造各段的活动水平也不尽相同．本文使用的模型较为简单，不能够反映复杂的地震构造特征及构造各段活动水平的差异．图６ａ为根据上述４个模型的５０年超越概率１０％的ＰＧＡ等权相加得到的华北地区ＰＧＡ分布．可以看出，该分布与使用《中国地震动参数区划图》宣贯教材（胡聿贤，高孟潭，２００１）中的综合潜源方案计算得到的华北地区５０年超越概率１０％的ＰＧＡ分布（图６ｂ）具有一定的相似性，它们都能够反映华北地区的地震危险性水平．然而，二者在某些方面也存在着较为明显的差异：前者计算得到的符合ＰＧＡ≥１００ｃｍ／ｓ２条件的区域面积明显要比后者的大，而符合ＰＧＡ≥２５０ｃｍ／ｓ２条件的区域面积则比后者的要小．
６讨论与结论
本文以地震带为统计单元计算地震活动性参数，以空间光滑的地震活动性模型为空间分布函数．采用４个模型计算了华北地区的地震活动性模型和地震危险性．得出如下结论：１）采用仪器记录地震计算得到的地震活动性模型和地震危险性结果能够反映华北地区现今的地震活动水平和地震危险性水平，符合人们对现今华北地区地震危险性的认识．２）采用历史强震（Ｍ≥４．７）计算的地震活动性模型和地震危险性结果，较好地反映了华北地区中强地震活动区的地震危险性水平．３）以地震应变计算地震活动率，并根据点椭圆模型和线椭圆模型计算得到的地震活动性模型，能够较好地反映空间大地震的活动水平和构造特征．根据采用４个地震活动性模型计算得到的地震危险性结果，加权平均得到的华北地区ＰＧＡ分布，并与采用四代图中综合潜源方案计算得到的华北地区ＰＧＡ分布作了比较，发现二者具有一定的相似性，即均能反映华北地区的地震危险性水平．然而，二者在某些方面也存在着较为明显的差异：首先二者几何形状的差异主要是由于潜在震源区的类型不同造成的；其次，前者计算得到的符合ＰＧＡ≥１００ｃｍ／ｓ２条件的区域面积明显要比后者的大，而符合ＰＧＡ≥２５０ｃｍ／ｓ２条件的区域面积则比后者的要小．这种差异可能是由以下几种原因造成的：①在我国现行地震危险性分析中，潜在震源区的划分主要有两条原则，即原地复发和构造类比．采用地震活动性模型计算得到的ＰＧＡ分布能够较好地反映根据原地复发原则划分潜源区的地震危险性水平；而对于根据构造类别原则划分的潜源区，由于现今地震活动水平较低，故使用地震活动性模型计算得到的ＰＧＡ不能反映这些地区的地震危险性水平．②Ｂｅａｕｖａｌ等（２００６）的研究发现，如果将潜在震源内部的地震活动水平看成是均匀的，则计算得到的地震危险性水平偏高的点要比偏低的点所占的比率大．③华北地震地区历史地震记录较短，而华北地震区强震复发间隔很长，对于Ｍ７．０以上大震，一般在千年以上，因此采用原地复发原则也会低估地震危险性水平．通过本文研究可以看出，根据地震目录，以空间光滑的地震活动性模型为空间分布函数，计算得到的地震活动性模型能够反映华北地区地震活动的空间不均匀性和地震构造特征，计算得到的地震危险性结果能够反映华北地震区的地震危险性水平．然而考虑到大地震的构造分布往往比较复杂，本文中采用的模型相对比较简单，不能够充分反映大地震活动的构造特征以及构造各段的地震活动水平．因此，在地震危险性分析中，对于高震级的潜在震源区只使用地震活动性资料来计算是不够的，应充分利用地质、地球物理、大地测量等方面的定量资料来综合描述．
