华北地区中东部位于燕山山脉以南，太行山山脉以东的广袤平原地区.该区域内包含多座分布于华北克拉通东块体与太行山区交界处的大型城市，受华北克拉通破坏活动的影响，区内拥有丰富的油气储备和地热资源(Yang et al., 2003; 高山等, 2009).富含自然资源对于城市建设和长期发展而言是积极的，但也表明该区的地质构造较为复杂.该区地震活动性较强，历史上发生过多次7级以上地震(图 1)，包括1626年7.0级灵丘地震，1679年8.0级三河—平谷地震，1966年7.2级邢台地震以及1976年7.8级唐山地震(赵博等，2013).参考历史地震分布，即将展开建设的雄安新区亦处于7.0级地震的空区内，1679年雄安新区内发生5.8级雄县地震，2006年邻近地区发生过文安5.1级地震(王武星等，2009).最新的地震危险性评估研究表明雄安新区抗震设防烈度应从7度调整为8度(杨百存等，2017).如上所述，华北地区中东部具有较强的地震活动性，同时又是我国的政治与经济发展重地，故该区域的地震活动性特征值得深入探讨.

图 1

Fig. 1

图 1 华北地区中东部的历史地震和构造环境 Fig. 1 The historical earthquakes and tectonic environment of Xiong′an and its adjacent area

地震活动的分布与断层的分布相一致，但其分布还应与岩石圈强度、均衡调整以及构造应力状态息息相关.岩石圈挠曲分析是指利用重力异常与地形数据推算岩石圈强度和以往均衡调整初始加载位置等信息的研究.岩石圈有效弹性厚度、均衡调整初始加载比分别代表着岩石圈强度与均衡调整初始加载的分布位置，是现今地形与重力异常数据所包含的岩石圈以往的均衡调整信息(Watts, 2001).垂向构造应力是现今岩石圈未达均衡状态而形成的地壳受力状态和其未来运动趋势(高尚华等，2016；佘雅文等，2017).研究岩石圈的力学属性，以及其过去的加载状态与未来的运动趋势，对于了解区域垂向构造活动具有重要意义，对城市的安全规划建设具有参考价值.

本文利用重力模型与地形数据，基于Fan小波的布格重力异常一致性方法(Kirby and Swain, 2011)与重力均衡调整方法(Fu and She, 2017; 佘雅文等, 2017)，研究华北地区中东部岩石圈有效弹性厚度、均衡调整初始加载比和垂向构造应力的分布特征，获取该区岩石圈均衡特性，并结合历史地震数据分析以上岩石圈参数与地震活动性的关系，为该区内的城市安全规划建设提供参考.

EIGEN-6C4是最新的基于重力卫星GRACE与GOCE以及地面和航空等重力数据融合而成的全球重力模型数据，模型精度较EGM2008重力模型有较大提高(Förste et al., 2014).ETOPO1为1′×1′分辨率的全球数字高程数据(Amante and Eakins, 2009)，被广泛的用于地学研究中.

首先，从EIGEN-6C4和ETOPO1中提取研究区域自由空气重力异常数据(图 2c)和地形数据(图 1).为了避免后续频域计算的边缘效应和保证小波方法足够的空间分辨率(Kirby and Swain, 2011)，计算过程中使用的数据范围相对于图 1向外延伸4°.同时，为了避免地球曲率的影响将数据转换到Lambert Conic投影(Snyder, 1987; Kirby and Swain, 2009)，投影中心点为东经116°北纬0°，数据采样率为20 km.

图 2

Fig. 2

图 2 华北地区重力异常数据及其地形与沉积层改正 (a)沉积层改正; (b)地形改正; (c)自由空气重力异常; (d)布格重力异常.黑色虚线为块体分界线. Fig. 2 The distribution of the gravity anomalies, sediment correction and terrain correction in North China (a) Sediment correction; (b) Terrain correction; (c) Free-air gravity anomalies; (d) Bouguer gravity anomalies. Theblack line denotes the boundaries of the blocks.

其次，对自由空气重力异常数据进行层间改正和地形改正，以获取研究区域的布格重力异常数据.地形改正基于分辨率为15″×15″的数字高程数据SRTM15_PLUS(Smith and Sandwell, 1997; Becker et al., 2009)，使用二维快速傅里叶变换计算方法(Li and Sideris, 1994; Kirby and Featherstone, 2002)，在保证精度小于1mGal的同时兼顾了计算效率，改正结果见图 2b.较大的地形改正值分布于地形剧烈变化的地区，最大改正值为15.1 mGal.

最后，基于Crust1.0(Laske et al., 2013)全球地壳结构数据，使用快速计算重力异常方法(Parker, 1973)，获取研究区域沉积层引起的重力异常(图 2a).从图 2a中可见，沉积层造成的重力异常最大值为-130 mGal左右，分布于渤海地区.为了消除沉积层对后续岩石圈挠曲计算的影响，故将其从布格重力异常中剔除(Chen et al., 2014)，得到沉积层改正后的布格重力异常如图 2d所示.图 2d中布格重力异常梯度带沿太行山脉东侧分布，其两侧的重力异常值相差近300 mGal，与该地区的地形(图 1)形成明显的负相关分布关系，这为后续的岩石圈挠曲分析研究提供了有力支持.

本研究基于CRUST1.0模型展开，该模型分辨率为1°×1°，对本研究计算结果的分辨率有一定程度的制约.部分学者依据地震波资料推算了该区域的沉积层厚度和地壳密度结构分布(罗艳等, 2008; Chen et al., 2009; 胡刚等, 2017；Zheng et al., 2017)，但相关研究对于浅层构造的约束不足，又受限于波速结构和密度结构之间经验关系的精度，给出的结果各不相同，不能有效提升该区域沉积层覆盖模型的分辨率.未来我们计划在该地区展开高密度重力和GPS联合观测，实质提高相关计算结果的分辨率和精度.

依据重力与地形数据之间存在转换函数关系，可计算岩石圈有效弹性厚度T_e.目前比较主流的方法是由Forsyth(1985)提出的利用布格重力异常数据与地形数据的一致性随波数的变化关系反演岩石圈有效弹性厚度.该方法首先将重力与地形数据转换到频域，并利用公式(1)计算观测数据的一致性；其次将两组数据和岩石圈有效弹性厚度作为输入，利用Forsyth(1985)提出的岩石圈弹性板均衡模型计算预测一致性；最后将观测和预测一致性进行拟合，利用全局搜索或者最小二乘的方法获取研究区域的岩石圈有效弹性厚度.

Pérez-Gussinyé等(2009a, b)使用滑动窗口方法计算岩石圈有效弹性厚度的分布，由于采用不同尺度数据窗口的反演结果不尽相同，最终的结果是利用不同尺度的窗口结果进行加权平均得到的.Kirby和Swain(2008, 2009, 2011)提出了基于Fan小波反演岩石圈的有效弹性厚度与加载比分布，该方法克服了采用固定尺度窗口而造成的数据频域信息不完整的缺陷，并且可以计算研究区域任意网格数据节点上的有效弹性厚度与加载比.

(1)

本研究采用公式(1)计算一致性Γ(k_f, X)，〈〉_θ代表对以θ为变量的数据求平均，θ为组成Fan小波的莫莱特小波的参数(Kirby, 2005).k_f为波数，其与波长λ的关系为λ=2π/k_f，X为空间位置. 和分别为布格重力异常与地形数据的小波系数，*代表取共轭.反演计算的目标函数为

(2)

其中，ReΓ²为观测一致性实部平方. ε_ReΓ²为ReΓ²的误差，以θ为变量采用剪刀误差估计方法(Pérez-Gussinyé et al., 2009a, b; Kirby and Swain, 2013)计算误差. ReΓ_c²为预测一致性实部的平方，在此本文使用Kirby and Swain(2011)总结的方法计算预测一致性Γ_c.在计算得到观测一致性后，设定岩石圈有效弹性厚度T_e的初始值，结合公式(2)使用阻尼最小二乘拟合方法反演T_e值.

加载比f为均衡调整的莫霍面与地表初始加载的比值，Kirby和Swain(2011)给出的表达式为

(3)

其中和分别为均衡调整初始地表与地幔加载的小波系数，ρ_c和ρ_m分别为地壳与地幔密度，取值为2800和3300 kg·m^-3.在此平均转换波长附近的三个f值作为该空间节点的加载比(Kirby and Swain, 2008)，转换波长是预测一致性为0.5时所对应的波长(Simons et al., 2003; Kirby and Swain, 2008).需要注意的是，当加载完全来自莫霍面时f应为无限大，为了表述方便本文的加载比为F=f/(1+f)；当F为1时代表加载完全来自莫霍面，为0时代表加载来自地表，为0.5时代表地表与莫霍面加载相等(Mckenzie, 2003, Kirby and Swain, 2009, She et al., 2016).

反演得到的华北地区中东部以及其周边区域的岩石圈有效弹性厚度和加载比分布如图 3(a, b)所示.图 3a中蓝线西侧为太行山地区，东侧为华北克拉通东块体.该区的岩石圈有效弹性厚度集中在10~65 km范围内，呈从东南向西北逐渐减小分布.华北克拉通东块体内部有效弹性厚度值较大，太行山地区的值较小，且其梯度带与块体边界有较高的一致性.郑勇等(2012)利用Fan小波方法计算的相同区域的岩石圈有效弹性厚度为5~40 km；杨亭等(2013)利用Moho地形导纳法反演的岩石圈有效弹性厚度结果为5~40 km；Chen等(2013)使用基于MTM谱分析的窗口方法计算了中国地区的岩石圈有效弹性厚度分布，对应本文研究区域的结果为10~60 km.以上结果表明，本研究反演的岩石圈有效弹性厚度分布结果与前人研究结果较为接近，其中相对于郑勇等(2012)和杨亭等(2013)的结果偏大.造成此差异的原因有两点：一是本研究使用了最新的EIGEN-6C4重力模型数据；二是本研究利用CRUST1.0数据去除沉积层的影响.

图 3

Fig. 3

图 3 华北地区中东部岩石圈有效弹性厚度、均衡调整初始加载比、垂向构造应力和大地热流分布 (a)岩石圈有效弹性厚度分布(Te); (b)均衡调整初始加载比分布(F); (c)垂向构造应力分布(σv)，箭头代表GPS垂向速度场; (d)大地热流分布.图中红色圆点代表公元前780年以来大于等于4级的历史地震，大小代表震级. Fig. 3 The distributions of the effective elastic thickness, the isostatic initial loading ratio, vertical tectonic stresses and the geothermal flow in middle east area of North China (a) Effective elastic thickness (Te); (b) Isostatic initial loading ratio (F); (c) Vertical tectonic stresses (σv), arrows denote the GPS vertical velocity field; (d) Geothermal flow. The red dots denote the epicenter of the earthquakes (Since 780 BC, M≥4), and the size denotes the magnitude.

图 3b为利用上述Fan小波方法计算的均衡调整初始加载比F分布，加载比高值主要分布于太行山以东的华北克拉通东块体内部，低值主要分布于西部太行山区和北部燕山山区.该结果表明太行山以东地区(华北克拉通东块体)均衡调整的初始加载主要来自于地幔，即现今的地形与重力异常的分布是由于之前的地幔物质上涌或地壳底部物质拆沉而形成的.同时，该加载比分布结果与朱日祥等(2012)厘定的由于太平洋板块俯冲为主因而造成华北克拉通破坏的位置(太行山以东地区)一致.

该区大地热流分布(图 3d)(汪集旸和黄少鹏, 1990; 胡圣标等, 2001; 姜光政等, 2016)表明，大地热流高值主要分布于太行山以东地区，低值主要分布于西北地区的太行山区，大地热流与均衡调整初始加载比分布较为一致.研究表明在地幔对流和俯冲板块后撤作用下，华北克拉通迅速减薄，地表热流升高(何丽娟和邱楠生, 2014; 张林友等, 2016)，因此大地热流和加载比的分布特征是互恰的，均是华北克拉通破坏的结果.

一般而言克拉通的大地热流相对于其他地区是较低的，大地热流较高的地区岩石圈强度应较低，而图 3a和3d表明，大地热流的分布随着有效弹性厚度的增大而增大.如渤海湾地区的高热流高有效弹性厚度的情况，该异常可能是华北克拉通破坏造成的.针对该区大地热流和有效弹性厚度的异常分布特征，值得在日后的工作中继续研究.

垂向构造应力是由于岩石圈未达到艾里均衡状态，地壳自重与地壳受到流体地幔浮力的应力差，其为地壳的均衡调整提供动力.垂向构造应力正值对应该区域具有隆升的运动趋势，负值代表该区域的运动趋势为下沉，为零代表该区域处于均衡状态.本文基于Fu和She(2017)以及佘雅文等(2017)的方法，利用CRUST1.0地壳密度分层数据计算了研究区域的垂向构造应力分布.该区垂向构造应力约为-20~20 MPa，中西部地区垂向构造应力向上，东北和西南地区向下. 图 3c中垂向构造应力表明，研究区域内中部和西部地区的运动趋势为隆升，其他地区以沉降为主.除中南部由于地下水开采造成的严重沉降地区外，已有的GPS垂向速度场(赵斌等, 2014)与垂向构造应力的分布较为一致.需要注意的是，GPS观测的沉降量比隆升量大一个量级(图 3c)，华北地区沉降主要由于地下水开采造成，构造活动的信号被淹没且不易提取，须在后续的工作中继续研究.

对华北地区中东部公元前780年以来4级以上的地震(数据来源为中国地震台网中心，其中4级以上地震为180个，5级以上地震为114个，6级以上地震为28个)及其对应的岩石圈有效弹性厚度、均衡调整加载比和垂向构造应力进行统计分析，结果如图 4所示.图 4(a, e, i)分别为研究区域岩石圈有效弹性厚度、加载比和垂向构造应力的统计分布.图 4(b，f，j)、图 4(c，g，k)和图 4(d，h，l)分别为该区公元前780年以来4≤M＜5、5 ≤M＜6和M≥6的地震震中所对应的岩石圈有效弹性厚度、加载比与垂向构造应力的统计分布，可见不同震级的分布特征较为一致.该区的岩石圈有效弹性厚度集中于10~15 km(图 4a)，地震震中所对应的岩石圈有效弹性厚度集中于10~25 km(图 4b，4c，4d).加载比分布集中于0.3~0.8(图 4e)表明该区均衡调整多是地表与地幔加载共同作用的结果，且地幔加载的情况占比大于地表加载；该区地震震中所对应的加载比分布(图 4f、4g、4h)表明，地震多发于加载比为0.4~0.6的区域.研究区域垂向构造应力分布(图 4i)集中于-15~20 MPa，地震集中分布于垂向构造应力向下的区域，且-10 MPa对应的地震最多，零值区域的地震最少(图 4j, 4k, 4l).

图 4

Fig. 4

图 4 华北地区中东部岩石圈有效弹性厚度、均衡调整初始加载比和垂向构造应力与地震活动性的统计对比 (a, e, i)为研究区域有效弹性厚度，均衡调整初始加载比和垂向构造应力的统计分布；(b, f, j)、(c, g, k)、(d, h, l)分别为4≤M＜5、5≤M＜6和M≥6的地震震中所对应的有效弹性厚度，均衡调整初始加载比和垂向构造应力的统计分布. Fig. 4 The statistical comparison of the effective elastic thickness, isostatic initial loading ratio, vertical tectonic stresses and seismicity in middle eastern area of North China Fig a, e, i show the statistical distributions of effective elastic thickness, isostatic initial loading ratio and vertical stresses. Fig b, f, j, Fig c, g, k and Fig d, h, l show the statistical distributions of the three parameters above corresponding to the epicenter of historical earthquakes with 4≤M < 5, 5≤M < 6 and M≥6, respectively.

从以上统计结果可知，研究区域的地震主要分布于岩石圈有效弹性厚度适中(10~25 km)的区域，而岩石圈有效弹性厚度低值或高值区地震较少.Watts(2001)和Chen等(2013)的研究结果亦得到类似的结论.该对应关系的可能原因是岩石圈有效弹性厚度低值代表岩石圈强度较低，表明该区域接近于塑性即在外力作用下易于变形，故不足以积攒产生强地震的能量；而高值代表岩石圈强度较大，表明该区域近于刚性，故不易发生地震.

加载比为0.7~0.8的区域总体占比较大(图 4e)，但地震的分布在该区域相对于加载比为0.4~0.6的区域较少(图 4f, 4g, 4h)，即表明均衡调整初始加载为莫霍面加载的区域中地震较少，而地表与莫霍面加载近似一致的地区地震较多.故据以上分析可推断，研究区域的地震活动性随着加载比的增加而降低.需要注意的是，研究区域内加载比为0.1~0.3的区域较小，该加载比区间对应的地震数量亦小，故该加载区间地震分布并未符合以上规律.

为了验证以上加载比与地震分布的对应关系，使用分形模拟方法(Peitgen and Saupe, 1988; Pérez-Gussinyé et al., 2009a; Kirby and Swain, 2009;)模拟地形与莫霍界面起伏(图 5a, b)，并结合多层弹性地壳均衡调整模型(Simons et al., 2003; 佘雅文等，2017)模拟计算了不同加载比所对应的均衡调整前后的地壳界面差异(图 5c, F=0.5, T_e=20 km时地壳的均衡调整前后差异).均衡调整模型假设地壳不可压缩，故地表与莫霍面的均衡调整前后的差异一致.对应不同的加载比(0~1，间隔0.1)各进行100组模拟计算，并求平均值与标准差，得到均衡调整前后差异随加载比的变化趋势(图 5d，黑色圆点为均衡调整量级，误差棒为100组模拟结果的标准差)，有效弹性厚度值取20 km.从图 5d可见均衡调整前后的差异随着加载比的增大而减小，即表明莫霍面加载为主的均衡调整的量级小于以地表加载为主的情况.均衡调整的量级即为以往垂向构造活动中地壳的运动量级，该量级越大表明该地区地壳活动性越强，故对应的地震活动性亦越强.结合以上模拟结果与研究区域内加载比与地震分布的关系，可推断均衡调整初始加载比与地震活动性成反比例分布，即地震活动性随加载比的增加而减弱.

图 5

Fig. 5

图 5 均衡调整模拟结果 (a)模拟地形; (b)模拟莫霍面; (c)均衡调整前后差异; (d)均衡调整量级随加载比变化趋势. Fig. 5 The results of the simulation of isostatic adjustment (a) Simulated topography; (b) SimulatedMoho interface; (c) The difference between before and after isostatic adjustment; (d) The trend of magnitude of isostatic adjustment changes with loading ratio.

垂向构造应力反映了现今岩石圈内部加载情况，并表明了地壳的运动趋势.图 4j、4k、4l表明，垂向构造应力零值区附近地震活动最弱，4~6级地震多分布于垂向构造应力向下的区域，大于等于6级的地震多分布于垂向构造应力向上的区域.据以上分析可推断，地震活动性随着垂向构造应力绝对值的增大而增强.需要注意的是，垂向构造应力小于-20 MPa和大于20 MPa的区域占比较小，故该区域并未符合以上规律.

具体到研究区域内的城市，德州、沧州和天津地区的地震活动性明显弱于唐山、北京、衡水和邢台地区(图 3).即将开展大规模城市建设的雄安新区的岩石圈有效弹性厚度大约为15 km，均衡调整初始加载比为0.5~0.6，垂向构造应力为15~20 MPa，并且该区的大地热流亦较高.结合上文对历史地震(1679年5.8级雄县地震，2006年5.1级文安地震)和岩石圈参数的统计分析结果(图 4)表明，雄安地区的岩石圈参数对应的地震活动性较强，该区具有发生5级以上地震的可能性，故在城市建设初期应做好防震减灾工作.

首先，基于EIGEN-6C4重力模型和ETOPO1全球数字高程数据，利用基于Fan小波的布格重力异常一致性方法研究了华北地区中东部岩石圈有效弹性厚度与均衡调整初始加载比分布.其次，使用重力均衡调整方法计算了该区的垂向构造应力分布.最后，将以上结果与该区地震活动性进行了统计分析，得到以下主要认识：

(1) 华北地区中东部的岩石圈有效弹性厚度为10~65 km左右，自东南向西北逐渐减小分布；均衡调整初始加载比为0.5~0.8，表明现今的岩石圈挠曲状态主要由莫霍面加载形成；垂向构造应力为-20~20 MPa左右，中西部地区垂向构造应力向上，东北和西南地区向下.

(2) 华北地区中东部的地震活动性随着岩石圈有效弹性厚度和均衡调整初始加载比的增加而减弱.垂向构造应力零值区域地震活动性较弱，地震活动性随垂向构造应力绝对值的增大而增强.

(3) 雄安新区的岩石圈有效弹性厚度大约为15 km，均衡调整初始加载比为0.5~0.6，垂向构造应力为15~20 MPa，岩石圈参数对应的地震活动性较强，故在城市建设初期应做好防震减灾工作.