重力异常是地壳或上地幔物质密度变化及构造运动特征的反映，从地表到地球深处所有密度分布的不均匀性是引起重力异常的主要地质因素^[1]。申重阳等^[2]研究发现，断层位错会引起断层周围局部重力场变化，而地震的孕育和发生往往与活动断层的运动紧密相关。2008年，美国国家地理空间情报局(NGA)发布超高阶全球重力场模型EGM2008，该模型基于GRACE卫星重力资料、地面实测重力资料及DTM地形等数据进行调和分析，完全阶次至2 159，空间分辨率为5′(约9 km)^[3]，使得全球重力场模型的空间分辨率提升到空前高度。许多学者利用该模型在全国各地开展了研究^[4-5]，而河北省作为受地震灾害影响严重的地区之一，目前尚无利用该重力资料对区域深部断裂开展的全面、系统性研究。本文将在EGM2008全球重力场模型基础上，利用向上延拓与导数换算等数据处理方法，结合地质构造资料对河北省张家口地区深部构造特征进行解释，讨论基于EGM2008模型研究河北省深部断裂系统的可行性。

1 区域构造背景与布格重力异常提取

研究区位于中朝准地台，主要发育有尚义-平泉断裂、张家口断裂、怀安-万全断裂、洗马林断裂等18条活动断裂。其中，近EW向构造形成较早，多被其他方向断裂切割或改造；NE向断裂是区域主要的控制性构造；NW向断裂规模较小，属配套构造，但活动较新，以断块差异运动和断裂活动为主。新构造运动时期，华北平原强烈裂陷，形成大型坳陷，裂陷区基岩面与山区北台面分布高程相比，差异升降最大幅度达14 000 m。控制断陷盆地的主断裂具有水平剪切性质，现代水系沿断裂线水平扭动，表明华北地区第四纪活动伴有明显的水平运动分量。第四系沉积等厚线与新近系等厚线形态基本相同^[6]，新构造运动在华北平原北部发育形成NW向沙河-武清-渤海第四纪活动断陷盆地带。隆起区表现为间歇式抬升，在华北山地的河谷内普遍发育3~4级阶地，在山前形成级数不等的台地^[7]。

张家口地区自有记录以来发生破坏性地震(包括余震)33次，其中M4.7~4.9地震14次，M5.0~5.9地震14次，M6.0~6.9地震5次。影响较大的地震分别为1337-09河北怀来6级地震、1618-11河北蔚县6级地震、1628-10河北怀安西洋河堡6级地震、1720-07河北沙城6级地震及1998-01-10河北张北6.2级地震。

EGM2008全球重力场模型自发布以来一直持续改进，本文通过世界重力图网站(http://bgi.omp.obs-mip.fr/)下载到最新的河北北部地区布格重力异常数据，结合相关活动断裂资料，绘制研究区布格重力异常分布图(图 1)。该异常分布图主要反映深部构造的演变趋势，从东南往西北异常幅值逐渐增大，说明莫霍面缓慢抬升^[8]，与该地区整体地势西北高、东南低的特征一致。

2 重力数据处理及断裂参数确定

布格重力异常转换处理和反演方法可分为褶积滤波(空间域)和FFT滤波(频率域)，在对深部断裂构造进行解释时通常使用FFT滤波，本文将采取向上延拓和导数换算两种方法对研究区深部断裂系统的构造特征进行解释。

2.1 重力场向上延拓

将观测面上的实测重力异常值换算到上半空间某一高度的方法称为向上延拓。该方法可压制浅表干扰异常或范围较小的局部异常，突出埋藏深度较大或规模较大的探测目标异常^[4]。

如果一个地质体在观测平面产生的重力异常Δg(ζ, η, 0)为已知，则可将这个观测平面看成一个无穷的单层物质面，其密度为：

$ \sigma(\zeta, \eta, 0)=\frac{1}{2 \pi G} \Delta g(\zeta, \eta, 0) $

(1)

该物质面在其上方空间任意点产生的重力异常与物质面下方地质体产生的重力异常等效，因此可利用该等效的单层物质面代替地质体，用于计算观测面上方的重力异常。一个密度不均匀的无限大物质面在其上方任意点(x, y, -z)产生的重力场为：

$ \begin{gathered} \Delta g(x, y, -z)= \\ G \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \frac{\sigma(\xi, \eta, 0) \mathrm{d} \xi \mathrm{d} \eta}{\left\{(\xi-x)^2+(\eta-y)^2+z^2\right\}^{\frac{3}{2}}}=\\ \frac{z}{2 \pi} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \frac{\Delta g(\xi, \eta, 0 z) \mathrm{d} \xi \mathrm{d} \eta}{\left\{(\xi-x)^2+(\eta-y)^2+z^2\right\}^{\frac{3}{2}}} \end{gathered} $

(2)

上式为实现空间域各项处理的基本公式^[4]。采用式(2)实现三度异常向上延拓的算法为：以平面上网格化的观测数据点为中心，点距Δx=2a为宽度，线距Δy=2a为长度，将计算区划分成M个矩形小区，假定每个水平矩形小区的密度均为常数，第k个小区的密度为σ_k，则式(2)可表示为：

$ \begin{array}{c} & \Delta g(x, y, -z)=\sum\limits_{k-1}^M \Delta g_k(x, y, -z)= \\ & \sum\limits_{k-1}^M\left\{\sigma_k \|-\arctan \right. \\ & \left.\frac{\left(-z \sqrt{(\xi-x)^3+(\eta-y)^3+z^3}\right)}{(\xi-x)(\eta-y)} \| \begin{array}{c} \xi t 2 \\ \xi t 1 \end{array}\right\} \end{array} $

(3)

上式就是实现三度体重力异常向上延拓的计算公式。

经查询相关资料，重力场上延高度和浅源顶深之间存在一定的经验关系：

1) 上延高度与浅源顶深总体上呈线性正比关系：当上延高度小于15 km时，浅源顶深约为上延高度的1/2；当上延高度为15~30 km时，浅源顶深变化较小，趋于10 km；当上延高度为30~40 km时，浅源顶深约为上延高度的2/5^[5]。

2) 上延高度越大，反映浅源顶深越大。

3) 上延到一定高度后，浅源顶深线性梯度变缓，该拐点上延高度应为区域深源场高度，该高度重力场为区域场。

根据反演结果结合本区实际情况，选取5 km、10 km、20 km、30 km这4个节点作为上延深度进行研究，利用Encom Modelvision Pro软件进行计算，结果如图 2所示。

2.2 重力场导数换算

导数换算是指将原布格重力异常换算为水平一次方向导数、垂向二阶导数或异常的总梯度模型、重力异常的解析信号等能更形象反映重力异常线性构造的一种数据处理方法。求取水平方向导数可突出线性构造在重力场中的反映，通常分析0°、45°、90°和135°等4个方向的水平导数，目的是突出与之垂直方向的线性构造，即0°反映近NW向，45°主要反映NW向，90°主要反映NS向，135°主要反映NE向线性(断裂等)构造。

对于一个密度不均匀且无限大的物质面，在其上方任意点产生的重力场为^[4]：

$ \begin{gathered} \Delta g(x, y, -z)= \\ G \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \frac{\sigma(\xi, \eta, 0) \mathrm{d} \xi \mathrm{d} \eta}{\left\{(\xi-x)^2+(\eta-y)^2+z^2\right\}^{\frac{3}{2}}}=\\ \frac{z}{2 \pi} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \frac{\Delta g(\xi, \eta, 0 z) \mathrm{d} \xi \mathrm{d} \eta}{\left\{(\xi-x)^2+(\eta-y)^2+z^2\right\}^{\frac{3}{2}}} \end{gathered} $

(4)

求其水平方向的导数，实质是求其在x方向的变化率，该变化率可近似表示为：

$ \left(V_{x x}\right)_0=\frac{\Delta g(\Delta x)-\Delta g(-x)}{2 \Delta x} $

(5)

令$\Delta g(x)=\sum\limits_{k-0}^m k a_x x^{k-1} $，并对x求导，可得重力异常的一阶导数表达式为：

$ V_{x z}(x)=\Delta g^t(x)=\sum\limits_{k=0}^m k(k-1) a_x x^{k-1} $

(6)

通过计算得到结果如图 3所示。

3 结果分析

根据上述计算结果，对研究区内5条主要断裂的规模、走向、切割深度、顶部埋深等进行分析，得出以下结论：

1) 尚义-平泉断裂：该断裂自尚义向东经崇礼、赤城、汤河口、下板城至平泉，活动时期从新近纪到第四纪，其中又以中晚更新世最为突出，全新世至今仍有明显活动。此外，该断裂南北两侧附近地震活动差异明显，断裂以北地震活动较弱，仅1998-01发生过张北6.2级地震，以南地震活动显著增强。在布格重力异常0°一阶方向导数灰度图(图 3(a))中表现为一条不连续的梯度带，可能受阴山断隆区与燕山断隆区两板块挤压干扰与互相切割所致；在红土梁乡与北新屯乡之间吻合度较高，断裂整体沿41°N线分布，与前人研究结果^[7]基本一致；在向上延拓20~30 km时，仍能看到较好反映，说明该断裂为一条深大断裂，切割深度为15 km以上。综上所述，尚义-平泉断裂为EW向的深部断裂带。

2) 怀安-万全北缘断裂：该断裂也称为怀安盆地北缘断裂，在135°一阶导数灰度图(图 3(d))中具有较好的反映，说明断裂总体走向NE。该断裂分为东西两段，其中西段断裂从沿渡口堡乡延伸至北沙城乡附近，走向近EW，向上延拓至5 km时已几乎无反映，说明该处断裂切割深度较浅；东部断裂由北沙城乡延伸至洗马林镇附近，走向NNE，向上延拓至10 km时仍能看到较明显反映，说明该部分比西段断裂切割深度大，属于较为活跃的部分，该特征与前人研究结果^[7]较为符合。

3) 怀安盆地南缘断裂：该断裂在布格重力异常0°一阶导数灰度图(图 3(a))中反映最为明显，在45°一阶导数灰度图(图 3(b))中也能看到微弱反映。经分析可知，该断裂大体走向EW，与相关资料描述一致。本文认为断裂西段部分呈NW走向，目前尚无资料对此进行描述，需进一步调查核实。在进行向上延拓时，通过观察可将该断裂分为东西两部分，向上延拓5 km时，该断裂西部等值线带几乎消失，反映已不明显，然而向上延拓10 km时断裂东部等值线带仍然存在，说明该断裂东部的切割深度远大于西部，属于相对活跃的部分，与前人研究^[7]中所描述的活动性有西弱东强的现象相呼应。

4) 阳原盆地北缘断裂：在布格重力异常图中能清晰找到该断裂所对应的等高线密集部分，分为东部、中部和西部3个部分，走向分别为近EW向、NE向、NEE向，与该断裂导数灰度图相对应。从布格重力等值线密集区与断裂位置对应的关系可判断，该断裂东部埋深最深，中部次之，西部埋深最浅。向上延拓5 km时，东部断裂已无反映；向上延拓10 km时，中部、西部断裂已几乎无反映，说明该断裂东部顶部埋深在2 km以内，中部、西部断裂顶部埋深在3~5 km之间。

5) 蔚广盆地南缘断裂：该断裂为蔚广盆地规模最大的断裂，对蔚广盆地的形成、发展起着决定性作用。布格重力异常图等值线密集区能明显地反映其NE向的整体走势，该特征在135°水平一阶导数灰度图(图 3(d))中也有明显反映。向上延拓5 km时，其西南段仍有明显反映，东北段反映则相对较弱；向上延拓10 km时，西南段反映已不明显，东北段反映消失，据此判断该断裂西南段顶部埋深约为5 km左右，西南段顶部埋深为2 km左右。

4 结语

本文通过获取的卫星重力资料计算得到布格重力异常，并在此基础上经过不同深度向上延拓计算与不同方向的一阶导数计算，得到研究区内5条断裂的规模、走向、切割深度等具体参数。结果表明，尚义-平泉断裂为EW走向的深部断裂带，阳原盆地北缘断裂东部顶部埋深在2 km以内，中部、西部断裂顶部埋深在3~5 km之间，怀安-万全北缘断裂与阳原盆地北缘断裂东部切割深度均比西部深，蔚广盆地南缘断裂西南段顶部埋深约为5 km左右，西南段顶部埋深为2 km左右。本文研究结果更加精细地刻画了张家口地区的深部构造环境，对以后河北省内重力推断断层工作具有参考意义。由于本文数据是基于卫星重力资料，虽测量范围大、资料齐全、成本低，但在测量精度方面无法与地面实测相媲美。若要对断裂性质进行更细致的探究，建议开展地面重力勘探，通过有针对性地布设更加密集的观测点来获取更高精度的地表实测重力资料，进而开展更加细致的分析研究。