公元1679-09-02北京东北部三河-平谷一带发生M8.0地震，震中烈度高达Ⅺ度，对包括三河-平谷在内的北京地区造成十分严重的地震灾害^[1-2]。地震之后1个月内发生10余次强余震，3个月后才逐渐减弱，且有感余震延续时间长达1 a多，因三河-平谷地震区特殊的地理位置及复杂的地质构造背景，一直以来受到学者的广泛关注^[3-7]。与其他地球物理学方法相比，重力测量具有方法简单、成本低、覆盖面积广、快速高效等优点，是区域构造和动力学研究的重要方法之一。

目前，大多数反演方法都是建立在参数反演的基础上，该类反演方法一般选择较简单的地质模型来拟合观测数据，但必须收集反演区域大量的先验信息，且该方法对层状结构的地壳模型的反演效果一般。而基于块体生长模式的重力三维反演算法不需要大量先验信息，对层状结构的地壳模型也具有较好的反演效果。

为研究三河-平谷地震区浅层构造背景、地震孕育机理以及地震与构造的关系，本文利用高精度重力异常数据，采用基于块体生长模式的重力三维反演算法对1679年三河-平谷8.0级地震区进行三维密度反演，得到地下密度异常体的分布特征及断裂构造的空间展布特征，结合区内强震活动分析地震的孕育及发生与断裂构造的关系，为该地震区的地质构造研究、资源勘探、地震预报及防灾减灾等工作提供依据。

1 重力三维反演 1.1 反演原理

为确定填充的规定密度差并拟合观测数据异常体的几何结构，将Δρ^-和Δρ⁺定义为地下岩体适当的规定密度差，对于整个异常结构而言，Δρ^-和Δρ⁺可以为相同值，也可以为相反值。将靠近拟合区域的地下岩体分解为小平行六面体(棱柱)的单元三维网格，以规定密度差填充选定的单元格，然后构建由m个填充单元组成的三维聚合定义的异常体模型^[8]。

由于单元格和规定密度差的选择具有非常高的自由度，必须用更快的方法代替对模型空间的一般探索。本文采用基于块体生长模式的反演算法，逐步选择单元格，将其填充并聚集到模型中。对于全局模型的拟合，重力异常与模型参数和残差v_i的关联方程可表示为：

$\begin{array}{*{20}{c}} {\Delta g_{i}-\left(\Delta g_{i}^{c}+A_{i j} \Delta \rho_{j}\right) f-\left(g_{0}+g_{x}\left(x_{i}-x_{M}\right)+\right.}\\ {\left.g_{y}\left(y_{i}-y_{M}\right)\right)=v_{i}, i=1, \cdots, n} \end{array} $

(1)

式中，A_ij为第j个棱柱单元在第i个观测点处单位密度的垂直引力^[9]；g₀、g_x和g_y为主要未知数；x_M和y_M为测量点坐标(x_i，y_i)(i=1, …, n)的平均值；Δg_i^c为与先前填充的单元格相对应模型的重力值；Δρ_j为指定的密度差；f≥1为未知的比例因子，用于将模型异常(Δg_i^c+A_ijΔρ_j)与观测异常(Δg)相拟合。

基于模型“拟合度”(残差的最小二乘最小化)和模型“平滑度”(总异常质量的L₂范数最小化)，采用混合最小化条件：

$\boldsymbol{v}^{\mathrm{T}} \boldsymbol{Q}_{D}^{-1} \boldsymbol{v}+\lambda f^{2} \boldsymbol{m}_{k}^{\mathrm{T}} \boldsymbol{Q}_{M}^{-1} \boldsymbol{m}_{k}=\min $

(2)

式中，m_k=(Δρ₁^k, …, Δρ_m^k)^T为第k个模型密度差的列向量；v为第k个残差的列向量；Q_D为重力数据不确定性的先验协方差矩阵；Q_M为模型参数不确定性的协方差矩阵，Q_M=A^TQ_D^-1A，该协方差矩阵能够使反演的模型位于适当深度，为反演方法的关键；λ为平衡模型拟合度和模型平滑度的因子。

对于每个新单元参数f、g₀、g_x、g_y，可以在每个求解步骤中进行调整，再通过式(2)评估单元格(填充有规定密度差)的适用性，最后选择产生最小值的第j个棱柱进行最终填充并聚合到生长体中。对于每个连续步骤，调整后的比例值f减小，附加参数g₀、g_x、g_y几乎达到稳定值。当f接近1时该过程停止，同时生成完整的三维异常体^[10]。

1.2 理论模型计算

在利用实际观测数据进行反演之前，需先建立理论模型来进行正演模拟计算，并将正演得到的理论数据用于反演，以验证反演方法的有效性及稳定性。

图 1(a)为简单的理论模型，该模型由2部分组成：1)密度差为500 kg/m³的长方形岩体，呈EW向分布，长600 m，宽400 m，高200 m，中心坐标为(1 200，1 000，300)；2)密度差为-500 kg/m³的负形球体，半径为200 m，球心坐标为(500，1 000，300)。模拟测点采用平面有规则布设(点间距为50 m，线长2 000 m)，假设所有重力数据精度相同，且不添加任何区域趋势(g₀=g_x=g_y=0)，正演计算得到的重力异常见图 2，该重力异常左边为负异常区域，右边为正异常区域。将测区地下岩体划分为60 m×50 m×40 m的棱柱，考虑密度差为0的初始模型，即无任何特殊初始信息，直接搜索规定密度差为500 kg/m³和-500 kg/m³的棱柱异常岩体，且异常体不断扩展，当f接近1时扩展停止^[11]。图 1(b)为调整后异常岩体的体积。

对比图 1(a)和1(b)可以看出，反演的密度体中心坐标跟拟合模型较为一致，且反演的密度体大小与拟合模型相近，表明该反演方法具有可行性，能够较为准确地确定地质体的空间位置。

2 三河-平谷8.0级地震区重力异常及其反演分析 2.1 研究区重力异常特征

三河-平谷地区(116.6°~117.5°E，39.8°~40.2°N)位于河北省平原地震带和张家口-渤海地震带交汇区附近，图 3为研究区布格重力异常(根据北京市地质矿产局物化探队综合研究室1986-12编制完成的1:10万《北京平原区重力编图报告》改绘)、主要断裂构造分布及1960~2019年三河-平谷地区地震分布。研究区内构造环境较为复杂，自北向南主要分布有二十里长山断裂(F1)、通县断裂(F2)、夏垫断裂(F3)及香河断裂(F4)，其中F1位于华北平原北缘，为燕山山脉和华北平原的过渡带；F2为顺义盆地和通县隆起的分界线；F3为盆地边界断裂，控制着大厂凹陷的形成和演化，为大厂重力低异常带和通县重力高异常带的分界线；F4为大厂重力低异常带和新集-侯家营重力高异常带的分界线。

区域布格重力异常自西北向东南呈低-高-低-高相间展布，重力异常高值区主要分布在F2以南、F3以北(布格幅值为5~20 mGal)及F4以南(布格幅值为8~20 mGal)，沿F1存在局部高重力异常(布格幅值为3~10 mGal)；其他地区均呈低布格异常分布(布格幅值为-24~0 mGal)。重力高异常区主要分布在通县-二十里长山和上仓镇-侯家营-宝坻一带，前者为大兴重力高异常带，可分为燕郊高异常区和三河-马坊高异常区，其中燕郊高异常区被F2和F3围限；后者称为宝坻重力高异常带。两个重力高异常带之间为大厂重力低异常带，主要分布于西集-大厂-三河一带，为NE向展布，长约30 km，异常中心位于大厂附近，极小值为-15.4 mGal。研究区西北部顺义重力低异常区属于北京重力低异常带，该异常带南翼缓、北翼陡。三河-平谷M8.0地震震中位于燕郊重力高异常区、三河-马坊重力高异常区与大厂重力低异常带之间的交汇过渡低值区域。

2.2 三维反演结果

反演得到的模型由大约50 000个单元组成，每个单元的平均边长为1 000 m，平衡系数为18.3，残差为159 μGal。重力异常主要反映沉积盆地的形态，图 4为三河-平谷地区不同深度的三维密度结构水平切片，在成像结果中，反演区域的密度对比在-300~200 kg/m³之间，其中该地区的主要地质构造如大厂盆地、顺义盆地、通县隆起等在成像结果中都清晰可见。

图 4(a)为2 km深度的密度异常分布。从成像结果可以看出，负密度异常区主要分布在顺义、平谷、大厂及邦均周边地区，正密度异常区主要分布在齐心庄、潘各庄、侉子店、燕郊以北及二十里长山等地区，该深度地震活动较少。

图 4(b)为4 km深度的密度异常分布。从成像结果来看，负密度异常区域主要分布在顺义-杨镇一带、徐辛庄西北部、大厂、皇庄、平谷及邦均周边地区，正密度异常区域主要分布在通县-燕郊一带、侉子店、二十里长山、三河东北部以及侯家营地区。该深度地震活动主要分布在宝坻重力高异常带和大厂重力低异常带之间，地震活动较少。结合图 4(a)和4(b)可以看出，正密度异常区域随深度的增加而逐渐增大。

图 4(c)为6 km深度的密度异常分布。从成像结果可以看出，负密度异常区域主要分布在顺义以南、徐辛庄西北部和大厂、皇庄、邦均及其北部地区，平谷地区的负密度异常在此深度已经消失；正密度异常区域主要分布在燕郊西北部、二十里长山、侯家营及三河东北部地区。该深度处地震呈带状分布，基本都沿夏垫断裂密集成带^[12]，且地震活动主要集中于马坊-齐心庄附近，位于夏垫断裂、通县断裂、二十里长山断裂交汇处，表明三河-平谷M8.0地震区的地震构造复杂，有可能由多方向的断裂活动造成^[13]。

图 4(d)和4(e)分别为8 km和10 km深度的密度异常分布。从成像结果来看，正密度异常区域增大，负密度异常区域逐渐减小甚至消失，说明此深度已为某些主要构造的最深处。该深度处地震震中也主要分布在马坊-齐心庄附近，小震活动比较密集^[14]。

图 5(a)为NW-SE向垂向剖面，该剖面穿过通县断裂(F2)、夏垫断裂(F3)及香河断裂(F4)。其中，夏垫断裂为首都圈东部地区一条NNE向岩石圈尺度的区域性深断裂带，全长超过100 km，总体走向为30°^[1]。该断裂带延伸较长，切割较深，是控制大厂凹陷形成和演化的边界断层，构成大厂凹陷和大兴隆起的分界线^[15]。在剖面上，该断裂表现为上陡下缓的铲式结构，断裂上小震密集分布，中强震较多，震源深度大多集中在8~10 km，且地震大多沿断裂走向分布^[15]，两侧较大的密度差表明其地质结构具有较大的反差且构造复杂，与三河-平谷M8.0地震的震源深度约10 km及震源区位于低密度体与高密度体交汇部位^[16]的结论一致。

图 5(b)为116.96°E处SN向垂向剖面，该剖面穿过二十里长山断裂(F1)及夏垫断裂(F3)。二十里长山断裂在地貌上表现为一排NW向的串珠状孤山，为北部燕山山脉和北京平原的过渡带^[4]，走向NW，倾向NE，倾角约为70°~80°。该断裂在形成和发展过程中受NE向断裂的限制向东南延伸，与夏垫断裂北段的分支断裂交汇。

2.3 小结

由三维密度异常分布图可以看出，研究区内断裂均为高密度异常带与低密度异常带的分界线，夏垫断裂与香河断裂深度约为10 km，并未向下切穿上地壳而进入下地壳。研究区内地震的震源深度较浅，大多为5~10 km，其他深度发生地震较少^[14]，地震分布最密集的区域为齐心庄-马坊附近。若将这些地震看成是三河-平谷M8.0地震产生的主破裂及次级破裂面的分布，可以发现，三河-平谷M8.0地震产生的破裂面具有多向性，说明该区的地震极有可能由多方向的断裂活动造成。

3 结语

本文对基于块体生长模式的重力三维密度反演算法进行模拟，并采用重力异常数据对三河-平谷地震区进行三维密度反演，得到地震区的浅部密度结构，能较好地反映地下异常体的密度分布、断裂构造及地下分层结构，并得出以下结论：

1) 通过对理论模型进行模拟，反演得到的异常体位置与理论模型相近，异常体的中心坐标与理论模型较为一致，表明反演方法能够准确地确定地质体的空间位置；

2) 研究区内断裂构造复杂，并在马坊-齐心庄一带形成复杂的交汇区，该区地震活动较为密集，三河-平谷M8.0地震就发生在断裂构造带的复杂交汇区；

3) 研究区位于浅部地壳块体的隆、凹相间处，地壳结构具有较大反差，该震区地壳仍处于剪切-挤压变形的构造环境中，这种地壳的差异变化控制着地震的孕育和发生；

4) 研究区内齐心庄-马坊地区为小震活动密集区，震源深度大多集中在5~10 km，重力反演初步结果显示，夏垫断裂切割深度也约为10 km，说明三河-平谷M8.0地震的震源深度为10 km左右。