0 引言

华北地区是中国重要的地震活动区，一直被作为地震重点监测区域。1975年以来，华北地区各省地震台网陆续完善，产出大量地震波形和到时数据。这些数据主要以2种形式记录：①首都圈数字记录台站投入使用之前模拟形式记录；②首都圈数字记录台站投入使用后的数字化记录。如何更好利用天然地震资源研究地球内部速度结构成为一项具有重要意义的工作。河北作为地震多发省份，曾发生多次破坏严重的大地震，如1831年磁县7.5级、1966年邢台7.2级破坏性大地震造成重大的人员伤亡和财产损失。目前，在邢台震区和磁县震区小地震频发。因此，研究河北南部及邻区的地壳速度结构，有助于揭示该区孕震机理。

地球物理学者对华北地区地壳结构的研究从未停止（王志砾，2005；Huang et al，2009），但因受到不同因素的影响，如地震定位精度、地震台站分布、数据质量等，研究结果一致性稍差，尤其对于河北南部及邻区的大震区域，地壳速度结构仍需进一步研究。不同的地球物理方法有其自身不足，比如：人工爆破是探测地壳速度结构的有效方法，精度高，但研究过程中需使用大量钻井工程和炸药物品，成本大幅度提高，应用受限制。天然地震具有能量大、深度大等特点，受地震位置精度及分布影响，短期内发生的地震很难满足高分辨率地壳结构研究的需求，可利用历史地震记录弥补震相数据不足、位置分布不均的不足。鉴于此，整理1980—2012年河北省南部及邻区的地震到时数据，采用震源位置和速度结构联合反演方法，获得该区地壳P波速度结构，并对中小地震震源位置重新定位，分析地壳三维速度结构特征，为研究地震孕育和地震动力学过程提供参考。

1 反演方法及解的评价

当地震发生时，利用Aki等学者的研究成果（Aki and Lee, 1976；刘福田，1984），给出地震波走时残差δt用线性化方程，即

$ \delta t=\Delta t+\frac{\partial t}{\partial x}\Delta x+\frac{\partial t}{\partial y}\Delta y+\frac{\partial t}{\partial z}\Delta z+\sum\limits_{n=1}^{N}{\frac{\partial t}{\partial {{v}_{n}}}\Delta {{v}_{n}}} $

(1)

式中：Δt为发震时刻，Δx为经度，Δy为纬度，Δz为深度，Δv_n为速度扰动量，N为速度参数总个数。若有i个地震及j个地震台站，则式（1）可表示为

$ \delta t=A\delta v+B\delta x $

(2)

式中：δt为地震射线走时残差，δv为节点处速度扰动量，δx是震源参数扰动量。

反演采用正交投影算子方法，把式（2）中的速度参数和震源参数解耦，给出方程组后求解，即

$ \left(I-{{P}_{B}} \right)A\delta v=\left(I-{{P}_{B}} \right)\delta t $

(3)

$ B\delta x={{P}_{B}}\left(\delta t-A\delta v \right) $

(4)

式中，P_B为与震源参数有关的从Rm到B的像空间R(B)上的正交投影算子，正交算子只与地震震源参数有关。分析速度参数和地震震源位置解耦后的方程组，可知速度的扰动量仅与其初值相关，与震源位置扰动量无关，而震源位置扰动量则与速度扰动量相关。

采用LSQR算法进行，反演无法直接给出解的分辨率矩阵。因此，对解的评价使用检测板方法（Humphreys et al，1988；Zhao et al，1992, 1994；周龙泉等，2007）。检测板测试主要过程为：给定初始速度模型参数，把各个网格节点值通过给定范围的扰动值进行正负相间扰动（本文计算扰动值取实际值的±3%），由实际地震射线数据正演计算获得相应的理论走时，加上给定的随机误差，作为新的射线走时数据进行反演（要求检测板测试方法反演过程与实际成像过程一致），比较反演结果与分辨率检测板的相似程度，作为解的评价。

2 资料选取

选取河北南部及邻区（34.0°—38.0°N，112.0°—118.0°E）作为研究区域，整合1980—2012年晋冀鲁豫交界区域地震数据，选取研究区域内62个地震台4 540个地震事件，获得27 709条质量较高的P波到时数据, 采用速度结构与震源位置联合反演方法，反演该区P波三维速度结构。

结合人工勘探结果获得研究区地壳速度结构初始模型，具体数值见表 1（祝治平等1995；房立华，2009；武岩，2011）。反演采用网格节点法对初始速度模型进行参数化，将研究区域以1.0°×1.0°均匀划分网格。采用连续函数表示模型速度结构，网格内波速度利用内插计算（Zhao et al，1992）。

3 P波速度结构特征

图 1展示了反演使用的地震事件和台站分布，研究区内台站分布比较理想，尤其在地震多发区域，台站密集；地震主要发生在邢台震源区、邯郸磁县震源区、山西地震带的太原附近及山东菏泽震源区，地震射线基本覆盖研究区域全部，有足够地震射线穿透震源区，为反演精度提供保障。图 1中绿色五角星分别为1966年3月22日邢台7.2级地震、1831年6月12日磁县7.5级地震、1983年11月7日菏泽5.9级地震震中位置。

3.1 水平层速度结构

选取的地震经过5次迭代反演，获得该区地壳P波速度结构。在此仅给出3 km、8 km、12 km、20 km深度水平层切片P波速度结构及检测板测试结果，见图 2。反演主要是使用直达P波数据，而震源深度30 km以下，地震发生相对较少，对于不同地壳深度P波射线穿透稀疏，分辨率检测板结果较差，因此主要研究30 km以上深度的地壳速度结构。

（1）在深度为3 km时，P波速度分辨率检测板测试结果较好，台站和地震密集区，分辨率效果较好，而东北部区域分辨率检测结果相对较差。该深度P波速度基本在4.0—5.0 km/s，平均约4.7 km/s，呈现2个高速区域波速约5.0 km/s：①邢台地震到磁县地震的震中区域，高速区面积较大；②菏泽震源区以南较小区域。

（2）在深度为8 km时，P波速度分辨率检测板测试结果较好，尤其在地震多发区域。该深度P波速度为5.5—6.2 km/s，平均约5.9 km/s，变化趋势与3 km基本一致，2块高速区域仍存在，与周边区域速度差值更加明显，高速体P波速度值约6.2 km/s，分布范围与3 km深度基本一致。

（3）在深度为12 km时，P波速度分辨率检测板测试结果较好。该深度P波速度为5.8—6.5 km/s，平均约6.1 km/s，高速区域仍存在，异常块体增大，位置略微南移，说明在不同深度高速区域不同。

（4）在深度为20 km时，P波速度分辨率检测板测试结果，除周边区域外，在地震多发区域较好。该深度P波速度为5.9—6.5 km/s，水平变化趋势与3 km、8 km、12 km深度明显不同，主要表现在：邢台震源区到磁县震源区下方高速区域转为低速体，低速体与太行山走向基本一致，P波速度值约5.9 km/s，区域周围约6.3 km/s；菏泽震源区下方高速体逐渐消失，且区域变小，位置向北偏移；位于山西断裂带太原附近的地震多发区下方出现高速异常体，P波速度值约6.5 km/s，呈现明显的横向不均匀性。

从图 2中不同深度的水平层速度分布可见，研究区P波速度呈现横向、纵向不均匀性，表现为高、低速体间隔形态，尤其在太行山前断裂带及苍东断裂带上，P波速度结构显示更加不均匀；邢台地震震源区下方存在低速体且上地幔上隆，王椿镛等(1993, 1994)、嘉世旭等（2005）的研究成果也佐证该特点，磁县震源区下方存在岩浆侵入，与该区反演结果中出现的低速体结果一致。

3.2 垂向剖面P波速度结构

根据地震震中分布的条带形态和断层分布产状，沿图 1中所示的剖面方向，绘制研究区域垂直方向的P波速度结构剖面，即沿磁县—邢台NNE向及垂向，沿菏泽NNE向及垂向，沿苍东断裂NNE向5条垂直剖面，基本代表研究区主要地震分布带。总体来看，每个方向的速度结构剖面均存在明显的不规则速度异常区域，特别是中地壳速度分布更加不均匀，明了地壳中地质构造的复杂性，沿太行山断裂带、聊城—兰考断裂带、菏泽断裂带在速度结构的垂直剖面均能较明显体现。

图 3中AA′剖面沿邢台地震到磁县地震方向，是研究区地震带，呈NNE向，与太行山走向基本一致；BB′剖面在邢台震源区下方通过，可更好展示其P波速度结构变化；CC′穿过研究区密集的NNE断层分布带；DD′剖面呈NWW向，经太原地震带、磁县震源区、菏泽震源区3个主要地震多发区域；EE′剖面穿过菏泽震源区，断层分布较复杂。

（1）AA′剖面在邢台—邯郸震源区下方8—16 km深度范围内存在一个高速异常体，长度约200 km，且高速体产状不均匀，根据位置判断，在邢台震源区和磁县震源区下方的高速体显著，P波速度值达6.4 km/s，比周边区域高约0.3 km/s。徐锡伟等（2000）利用邢台震区地震资料，结合人工地震勘探资料（邵学钟等，1993；王椿镛, 1993, 1994；武岩，2011），研究表明，在邢台震源区和磁县震源区下方的地壳水平分层，P波存在高低速间隔形态；对于磁县震源区下方P波速度结构，江娃利等（1994）、刁桂苓等（1999）、张小涛等(2008, 2011)，研究结果均表明：磁县震中深度约15 km，且在震源区下方存在P波速度高速体，磁县地震发生在高速体下。据历史地震资料，该区中小地震空间分布与断层一致，震源深度分布与速度异常体深度相关。

（2）BB′剖面。邢台震源区附近下方8—23 km深度范围内，地壳P波速度结构呈明显的高倾角态势，据前人（林真明等，1990；顾梦林等，2000）对邢台震源区下方速度结构的研究，莫霍面呈局部平缓隆起状态，隆起上部岩石主要表现为脆性破裂性质。因此，邢台震源区的余震主要发生在深度23 km以上部位，与本文反演结果一致。

（3）CC′剖面。P波速度呈现明显横向不均匀性，与前人（徐锡伟等，2000；嵇少丞等，2009；杨文采等，2011；张路，2012）研究成果一致。

（4）DD′剖面。菏泽震源区下方深度15—20 km处存在低速异常体，P波速度值为6.2 km/s，在太行山山脉下方，地壳呈高速隆起状，表明太行山区地下速度结构与东部平原差异较大，P波速度结构横向差异明显。

（5）EE′垂直剖面。明显可见菏泽震源区下方存在2个P波速度异常结构体，表现为高低速相间隔分布特征。

反演结果表明，研究区P波速度呈现明显的横向差异，且不同深度差异不同，太行山区下方速度结构与东部平原差异明显，表明西部的太行山区下方P波速度结构变化更加复杂；在邢台地震、磁县地震、菏泽地震的震源区下方均存在速度异常体，且大地震基本发生在速度异常体内部或边缘区域。在约6—10 km深度上，存一个P波高速异常体，长度约200 km，在深度10 km处基本消失；磁县台震区西部和南部存在明显低速体；邢台震区下方7—13 km处存在P波速度异常体，震中位于异常体边界，磁县震区下方同样具有复杂的P波速度分布，1931年磁县7.5级地震震中位于高低速交界处。

4 地震定位

利用反演得到的地壳三维P波速结构度模型，对研究区中小地震进行重新定位，结果显示，地震P波走时残差和震源位置有较大改善。图 4给出重新定位前后震源深度分布，可见定位前震源深度在10 km、15 km的地震约1/4，且深度分布不合理，经重定位，震源深度基本成正态分布，主要分布在10 km左右，深度分布更加合理，集中在5—16 km；重新确定震源参数后，地震波走时均方根残差（RMS）从1.68 s降到0.82 s，震源位置平均偏差为：东西方向1.4 km，南北方向1.3 km，垂直方向2.0 km。

研究区内中小地震重新定位后平均深度为11.9 km，比定位前的9.5 km增加2.4 km，这是因为定位前较多地震的模拟记录对深度控制能力较弱，出现较多直接给定10 km作为震源深度的地震。在参与反演的4 500多个地震中，有700多个地震震源深度给定为10 km，使用反演后的P波速度模型重新定位，定位结果精度提高，震源深度接近正态分布。

由图 1可知所选地震集中分布在AA′DD′剖面两侧，图 5给出沿2个剖面的震源深度分布，绘图时将AA′、DD′剖面两侧35 km内震源位置分别投影到剖面上。图中不同的颜色表示不同震级，其中黑色表示1.0—1.9级、绿色表示2.0—2.9级、黄色表示3.0—3.9级、粉色表示4.0—4.9级、红色表示5.0—5.9级。沿AA′剖面，震源深度自西南向东北逐渐加大，磁县震源区中小地震震源深度小于邢台震源区，且地震数量更少，由地震地质构造可知，2个震源区的地质构造不同，王椿镛(1993, 1994)研究表明，邢台震区下方上地幔上隆，可能存在岩浆浸入，从而造成岩石脆弱，导致地震多发。沿DD′剖面，太原、磁县、菏泽3个地震多发区，震源深度明显不同，自西北向东南呈现逐渐减小态势，联合该区地震地质构造可知，3个区域地质构造分别形成于古生界、新近系、全新统，呈现出地质年代越早，震源深度分布越浅的负相关特征。

所选地震重新定位结果表明：①P波走时残差、震源位置有较大改善。地震定位前主要集中在0—10 km内，重定位后震源深度分布更加合理，基本成正态分布，主要分布在5—16 km，且地震波走时均方根残差（RMS）下降0.86 s；②研究区中小地震条带分布显著，其中：第1条沿邢台—磁县震区呈NEE向展布，震源深度主要分布在7—14 km；第2条呈NNW方向，垂直于邢台—磁县地震带，震源深度主要分布在8—18 km；第3条沿交城断裂—系舟山北麓断裂呈NNE方向条带分布，震源深度主要在7—15 km。

5 结论

本文在整合历史地震到时数据基础上，采用P波速度结构与地震定位联合反演方法，获得河北省南部及邻区地壳P波三维速度结构，并对中小地震进行精确定位。研究区反演结果表明，地壳P波速度结果可信，可较好分辨各速度结构体。

P波速度结构分析表明：①研究区地壳P波速度结构在不同深度上均呈现横向不均匀性，尤其在10—25 km深度范围内，横向不均匀性更加突出；②深度10 km处，在邢台地震—邯郸磁县震源区下方存在水平层状高速异常体，该异常体速度分布在太行山区与东部平原差异明显；③该区断裂带分布与P波速度异常体的产状明显相关，大地震基本发生在异常体或P波高低速交界区域。

中小地震重新精定位结果显示：①重新定位后地震P波走时残差和震源位置均有较大改善；②中小地震条带分布明显，震源深度分布与地质构造年代具有一定的负相关性。