引言

龙门山断裂带位于青藏高原东缘、四川盆地西北缘，长约500 km，宽约30—50 km，沿NE−SW向展布，断层滑动以逆冲为主兼具右旋走滑分量。如图1所示，龙门山断裂带的主要断裂自西向东分布有后山断裂（汶川—茂县断裂）、中央断裂（映秀—北川断裂）、前山断裂（江油—灌县断裂）（唐荣昌等，1991；邓起东等，1994；贾东等，2003）。在该断裂带上相继发生了2008年汶川M_S8.0地震和2013年芦山M_S7.0地震，两次地震的地震成因和该断裂带的后续地震危险性研究备受关注。构建三维地质模型能够实现浅层地质调查和深部构造资料的综合，对岩石圈（地壳）结构进行直观描述，是进行孕震环境和地震机制分析研究的前提。

图  1  龙门山断裂带地质构造背景及剖面分布

Ⅰ ：竹巴龙—资中剖面；Ⅱ ：花石峡—简阳剖面；Ⅲ ：黑水—三台剖面；Ⅳ ：唐克—蒲江—阆中三角剖面. F_1-1：汶川—茂汶断裂；F_1-2：平武—青川断裂；F₂：映秀—北川断裂；F₃：灌县—安县断裂

Figure  1.  The geological structure of Longmenshan fault zone and profile position

Ⅰ ：Zhubalong-Zizhong profile；Ⅱ ：Huashixia-Jianyang profile；Ⅲ ：Heishui-Santai profile；Ⅳ ：Tangke-Pujiang-Langzhong profile. F_1-1：Wenchuan-Maowen fault；F_1-2：Pingwu-Qingchuan fault； F₂：Yingxiu-Beichuan fault；F₃：Guanxian-Anxian fault

下载: 全尺寸图片 幻灯片

近二三十年关于龙门山断裂带的研究已经累积了丰富的地球物理勘探资料，如深地震测深剖面（陈学波等，1988；韩渭宾，蒋国芳，2000；Wang et al，2007；滕吉文等，2008；徐朝繁等，2008；朱介寿，2008；Li et al，2009；Jia et al，2014；王夫运等，2015；张新彦等，2017）、基于接收函数反演或层析成像的速度结构模型构建（王椿镛等，2002，2003；刘启元等，2009；郭飚等，2009；雷建设等，2009；吴建平等，2009；Li et al，2009；邓文泽等，2014）、基于重磁数据的密度模型反演（邓起东，1994；滕吉文等，2008；张季生等，2009；唐新功等，2012；张恩会等，2015；杨光亮等，2015）以及大地电磁探测（李立，金国元，1987；孙洁等，2003；陈高等，2006；徐朝繁等，2008；王绪本等，2008，2009）。

龙门山区域的研究成果虽然丰富，但呈现单一分散的状态，针对这些已有成果进行客观整理的报道甚少。为此，本文拟将龙门山区域具有代表性的二维剖面作为数据源，通过专业的三维建模软件构建该区域深达70 km的三维地质模型，并通过实测的布格重力异常数据检验模型，继而通过重力资料约束反演对该模型进一步予以修正，为有效预测龙门山未知区域提供参考。

1.   构造背景

龙门山断裂带处于华北地块、华南地块与羌塘地块的结合部位，西北侧是松潘—甘孜古残留洋盆，东南侧是稳定的扬子克拉通（Jia et al，2014）。龙门山断裂带及邻近区域的速度结构反演结果表明：该区域的速度结构从青藏高原东南部的山地到扬子地台具有强烈的不均匀性；四川盆地所处沉积层较厚，且其沉积环境自晚古生代以来变化不大（王椿镛等，2002），地壳结构整体稳定，地壳及上地幔的平均速度均明显偏高（王椿镛等，2015），其中上地壳呈现低速异常，表明该地区处于坳陷状态（胥颐，2009）；松潘—甘孜地块地壳因介质岩性变化而呈增厚状态（王椿镛等，2015），其上地壳速度较高，中下地壳速度偏低，上地幔也表现为明显低速异常（胥颐，2009；王椿镛等，2015）。

2.   数据与方法

三维地质建模是一门通过空间信息管理、分析与地穴统计技术构造三维地质体，并进行地质解释的多学科教程技术（潘懋等，2007）。地质单元具有多尺度的非均匀性和连续性，根据地质单元的连续性可以从测得的部分信息推算其相邻区域的趋势。但是，这种推算方法的主观性较强，建模过程会掺杂建模人员的主观判断，因此建模的结果具有多解性，即模型的随机性。

建模是基于一定的算法将二维资料推演成三维空间信息的过程。首先，在模型表面上确定剖面与断裂带的位置，在模型表面约束剖面的地表迹线、断层的走向及影响范围；然后，由剖面迹线建立对应的剖面，对剖面信息进行数字化，剖面信息包含断层在剖面位置向下延展轨迹数据、层位在该剖面的具体埋深及形态信息；最后，依据适合的算法，由二维剖面生成三维模型。实际建模过程中需要人为干预，协调各剖面信息以确保数据的一致性。对模型赋予物性参数并进行正演，与实测的重力数据进行对比，从而验证模型的可信度。

本文使用GeoModeller软件（4.0.5版）进行三维地质建模，该软件支持输入不同类型的数据，如二维剖面、地质图、钻孔数据等，支持对数据的可视化和编辑，定义不同地质单元的接触关系后，使用势场插值的方式构建三维模型（Calcagno et al，2011）。本文取研究区范围为（102°E—107°E，29°N—34°N），建模内容包括断层模型和地层模型。用于比对的布格重力原数据是（98°E—108°E，26°N—34°N）范围内的采样间隔为5 km的川滇地区格网数据。

2.1   断层参数

用于三维建模的断层参数包括地表迹线和断裂倾向、倾角。断层地表迹线根据模型分辨率进行了适当的简化。

本文只考虑前山断裂带、中央断裂带和后山断裂带，忽略其它断层的影响。龙门山后山断裂包括北段的平武—青川断裂（F_1-2）和南段的汶川—茂汶断裂（F_1-1）。平武—青川断裂走向N60°—70°E （贾东等，2003）；汶川—茂汶断裂中段走向为N25°—45°E，大约10 km深度以上产状的倾向为NW，倾角约为55°，在5 km深度处开始变缓，倾角约为35°，在北端与平武—青川断裂斜交（唐荣昌等，1991；马保起等，2005；冯杨洋，2016），汶川—茂汶断裂南段的走向为N45°E （徐锡伟等，2008）。

龙门山中央断裂，即映秀—北川断裂，分为北段和南段两部分。北川—映秀断裂北段（F_2-2）的走向为N45°E，倾向NW，整体倾角约为50°—70°，在平均19 km深度以上倾向NW，倾角约为55°，在6 km深度处开始变缓；北川—映秀断裂南段（F_2-1）的走向为N35°—45°E，倾向为NW，平均17 km深度以上倾向NW，倾角约为55°，在9 km深度处开始变缓（邓起东等，1994；赵小麟等，1994；冯杨洋，2016）。

龙门山前山断裂（F₃）又称灌县—安县断裂，总体走向呈N35°—45°E，断面倾向NW，倾角为50°—70°，由南到北依次为F_3-1，F_3-2，F_3-3，F_3-4：南段的F_3-1和F_3-2走向为N43°E，倾向为NW，倾角45°—65°不等，在10 km深度以上倾向NW，倾角约为50°，在4 km处开始变缓，倾角为25°；中段F_3-3 （灌县—江油断裂）主要发育在中生代地层中，倾向NW，倾角较陡，产状为北东陡南西缓；北段F_3-4在平均14 km深度以上倾向NW，倾角约为45°，在8 km处开始变缓（邓起东等，1994；冯杨洋，2016）。

2.2   地层参数

地层参数包括地层界面形态和地层物性参数。地层密度根据密度与速度的转化公式（冯锐，1985）和实测剖面的地震波速度推算，结合重力反演最终确定，详见表1。地层界面形态根据收集到的剖面获取，经过筛选，最终使用以下四个剖面作为建模数据源：唐克—蒲江—阆中三角测线DSS剖面Ⅳ（陈学波等，1988；邓阳凡等，2011）、黑水—三台剖面Ⅲ（邓阳凡等，2011）、花石峡—简阳剖面Ⅱ（王有学等，2005）和竹巴龙—资中剖面Ⅰ（王椿镛，2002，2003；Jia et al，2014），各剖面的地表迹线见图1，图2为经过修改和数字化的剖面。

表  1  地层密度

Table  1.  Formation density

	地层	密度/（g·cm−3）
	沉积层	2.36
	上地壳	2.67
	中地壳	2.9
	下地壳	3.14
	地幔	3.35

下载: 导出CSV 

| 显示表格

图  2  本文采用的五个数字化剖面

（a） 竹巴龙—资中剖面（改自Wang et al，2007）；（b） 唐克—蒲江剖面（改自韩渭宾，蒋国芳，2000）；（c） 蒲江—阆中剖面（改自韩渭宾，蒋国芳，2000）；（d） 阆中—唐克剖面（改自韩渭宾，蒋国芳，2000）；（e） 花石峡—简阳剖面（改自王有学等，2005）；（f） 黑水—三台剖面（改自王成，2012）

Figure  2.  Five digital profiles used in this paper

（a） Zhubalong-Zizhong profile （revised from Wang et al，2007）；（b） Tangke-Pujian profile （revised from Han and Jiang，2000）；（c） Pujiang-Langzhong profile （revised from Han and Jiang，2000）；（d） Langzhong-Tangke profile （revised from Han and Jiang，2000）；（e） Huashixia-Jianyang profile （revised from Wang et al，2005）；（f） Heishui-Santai profile （revised from Wang，2012）

下载: 全尺寸图片 幻灯片

2.3   重力数据

本文重力数据购买于全国地质资料馆。中国地质调查局发展研究中心于2000年对原有的各行业间的数据进行整合，按照统一的重力基本网（1985）、坐标系（北京54）和国家高程标准（1985）、正常重力场公式、地改半径和中间层密度，建立了规范的、安全的全国范围重力数据库。该数据库收录原地矿部共计255个1:20万图幅的1:20万数据、84个1:20万图幅的1:50万数据、27个1:100万图幅的1:100万数据以及1999年以来 “国土资源大调查” 专项部署完成的青藏高原6个1:100万图幅和遍布全国的78个1:20万图幅的区域重力调查成果数据。

2.4   方法

本文利用Geomodeller进行建模。该软件支持多种地质数据的导入、可视化与编辑，在定义了描述不同地层单元之间关系的地层柱状图后，利用势场插值法计算三维模型。Geomodeller基于势场理论进行三维插值，使用协同克里金算法计算接触点（地质、断层）和方向数据，并将其视作三维标量场的等势面增量和导数（Lajaunie et al，1997）。根据不同接触关系定义地层层序，定义曲面为不同场的特殊等值面，根据地层层序的定义生成三维地质模型（Calcagno et al，2011）。

地球物理正演是指在地球物理资料解释理论中，由地质体的形状、产状、空间位置等赋存状态和密度、磁性、电性、弹性、速度等物性参数计算该地质体引起的场异常过程。正演模型是利用计算机模拟方法来获得已知地质体的地震响应特征，从而指导地震资料的构造和岩性解释。与反演相比，正演的优势在于改变模型参数方便，选用的计算方法灵活，模型制作简单。

正演算法采用多线程空间卷积算法：给定分辨率，将连续的地质模型离散化为多个正则或半正则三维体元；根据已知的观测值计算每个体元的权重，并将其作为卷积算子；对模型进行空间卷积，统计不同体元的权重总和，得到模型各点的重力异常；评估每一次正演结果的残差，迭代上述计算，当残差值达到最小即停止迭代（Gibson et al，2013）。

3.   讨论与分析

3.1   模型及重力异常残差

最终建成的地层模型和断层模型分别如图3a和图3b所示。由图3a可见，沉积层整体较薄，受分辨率的影响并未表现出较大的变化；地壳层以龙门山断裂带为界，两侧有明显的厚度差异，其中下地壳贡献的厚度变化值最大，上地壳次之，中地壳厚度变化最不明显；本模型的地幔层主要反映上地幔，在龙门山断裂带西北侧的松潘—甘孜地块埋深较大，经过龙门山断裂带的抬升，在东南侧的四川盆地上地幔约有10—20 km的抬升。

图  3  地层模型（a）和断层模型 （b）

Figure  3.  Fault model （a） and formation model （b）

下载: 全尺寸图片 幻灯片

图  4  模型正演结果及残差

（a） 实际观测值；（b） 模型正演计算值；（c） 模型残差

Figure  4.  Model forward result and misfit

（a） Observed gravity；（b） Computed gravity from model forward；（c） Model misfit

下载: 全尺寸图片 幻灯片

实测剖面（王成，2012）显示，位于四川盆地的龙泉驿、成都、彭州至都江堰段莫霍面埋深处于45—48 km之间，郸县、双流、资阳、简阳、成都至汉川段莫霍面埋深介于42—44 km之间；位于龙门山断裂带的耿达（埋深48 km）处莫霍面陡降至达维的68 km，经小金、丹巴缓慢上升至64 km，越过鲜水河断裂，道孚附近的莫霍面变浅至60 km，在茂县（45 km）陡降至黑水的58 km，莫霍面出现了13 km的陡降带，随后，经壤口、红原、唐克、麦溪缓慢上升至50 km。模型整体趋势与实测数据吻合。

对模型进行正演以验证其可靠性，参数见表1。布格重力异常原件的重力值范围为−76.17×10⁻⁵—−435.15×10⁻⁵ m/s² （图4a），计算的重力值范围为−172.88×10⁻⁵—−424.23×10⁻⁵ m/s² （图4b）。根据总体形态，计算场与观测场具有一致性，能反映断裂沿东北—西南对角线分布，但根据数值分析，残差较大（图4c）。如图所示，观测值西侧边缘表现为带状的重力低值不规则圈闭，计算值表现为存在于西北角的小面积重力低值圈闭，根据竹巴龙—资中测线（Wang et al，2007）的解释剖面，本文推断西南角的重力低值与埋深17—20 km处存在的延伸至龙门山断裂带的低速异常体有关；南侧中部偏北的重力低值与花石峡—简阳剖面（王有学等，2005）上埋深25—40 km处的低速异常体、阿坝—双流测线（朱介寿，2008）上埋深20—30 km处低速异常体有关。至于龙门山断裂带，观测值中显示出一条明显的重力梯度带，与计算值形态基本吻合，但断裂带西南端和东北端的重力梯度带表现为明显的向内侧收缩，出现这种差异我们推测原因主要有两点，其一是由于断裂带两端缺少实测剖面数据约束，使得曲面插值算法在模型边缘产生了误差，其二是模型正演在边缘处存在一定的误差。

根据实际建模中遇到的误差问题，总结如下：① 剖面数据由人工反演获取，反演的算法与人工解释的过程均存在一定的误差；② 在人工地震测深等试验下，剖面位置均不是规则的线，后期处理中经过一定程度的加工；③ 模型建造中需要将不同剖面进行匹配，为了保证整个模型的一致，需要人为修改部分剖面上的信息。

从整体来看，模型正演得到的重力异常展布与观测值基本一致，即整体重力值呈现为自西北至东南单调上升的重力异常梯度带，龙门山断裂带体现为高梯度变化带，这表明速度剖面构造的正演模型能很好地反映重力场的变化趋势。

3.2   各地层正演结果

通过调整物性参数的设置，排除其它各层的影响从而实现各地层的单独正演，结果如图5所示。

图  5  构建模型各地层的正演结果

（a） 沉积层；（b） 上地壳；（c） 中地壳；（d） 下地壳；（e） 地幔

Figure  5.  Forward results of each formations involved in the built model

（a） Sediment；（b） Upper crust；（c） Middle crust；（d） Lower crust；（e） Mantle

下载: 全尺寸图片 幻灯片

1） 沉积层。沉积层引起的布格重力异常整体表现为−45.53×10⁻⁵—−444.63×10⁻⁵ m/s²范围内的负异常，西北高、东南低，分区明显。沿龙门山断裂带形成了一条重力梯度带，该梯度带将整个区域划分为东西两侧，松潘—甘孜盆地沉积层薄、四川盆地沉积层厚。

2） 中上地壳。上地壳和中地壳皆显示为负异常值，前者变化范围较大，在−49.29×10⁻⁵—−352.03×10⁻⁵ m/s²之间，后者变化范围较小，在−89.14×10⁻⁵—−117.68×10⁻⁵ m/s²；两者形态差异较大，上地壳大致呈现为自西北至东南的梯度带，龙门山断裂带表现为明显且粗糙的高梯度变化带，而中地壳呈现为多个无规律的圈闭，或重力高圈闭，或重力低圈闭。

上地壳西北角表现为低重力异常，与奔子栏—唐克埋深在10—25 km范围内速度为5.7—5.9 km/s的低速异常体有关（王椿镛等，2003）；东侧中部的阿坝至理县（朱介寿，2008）、花石峡附近（王有学等，2005；邓阳凡等，2011）及黑水附近（蔡学林，2007）的中地壳存在低速异常体，与正演结果所显示的重力低值圈闭相对应。

3） 下地壳。较中上地壳和沉积层，下地壳引起的布格重力异常变化趋于光滑，整体为正异常，西北高、东南低，范围为26.37×10⁻⁵—52.73×10⁻⁵ m/s²，变化较小，龙门山断裂带表现为光滑的异常梯度带。下地壳异常变化小，说明厚度变化小，整体厚度较大，平均厚度达30 km。

4） 莫霍面。莫霍面的重力异常逐渐趋于简单，局部重力异常现象消失，反映了上地幔速度结构的简单化。重力异常沿龙门山断裂带表现为一条重力梯度带，梯度值介于75.08×10⁻⁵—257.82×10⁻⁵ m/s²。模型内布格重力异常数据呈现为规则的自西北至东南单调上升的重力异常梯度带。地幔层的正演结果与利用重力数据反演获取的4阶小波逼近图（深度等同于地幔层）的趋势一致（姜文亮，张景发，2011）。

龙门山断裂带西北侧的松潘—甘孜地块夹持于昆仑断裂带与龙门山断裂带之间，莫霍面在整体上呈凹陷形态，且地壳自西北向东南方向逐渐变薄，表现为台阶状抬升，莫霍面深度为60—70 km，东南侧四川盆地的莫霍面深度经龙门山的抬升，其埋深陡增至40—50 km （王有学等，2005）。

3.3   模型反演

为降低残差，在重力数据约束下对模型进行反演，根据残差计算（图6a−c）可知反演结果良好，图6d−f为反演后的模型展示，沿断层走向莫霍面的抬升情况很明显，与上述分析一致。

图  6  反演分析及模型展示

（a） 实际观测值；（b） 模型反演计算值； （c） 反演残差；（d） 反演模型；（e） 南北向切片展示；（f） 东西向切片展示

Figure  6.  Inversion and model presentation

（a） Observed gravity；（b） Computed gravity from model inversion；（c） Inversion misfit； （d） Inversion model；（e） North-south slice；（f） West-east sections

下载: 全尺寸图片 幻灯片

4.   讨论与结论

本文挑选整合了龙门山区域的经典文献研究资料，从中挑选了部分二维剖面，经过数字化、检校编辑、地层定义等步骤，构建了龙门山区域的三维地质模型，利用该区域实测重力数据验证，结果证明了模型的可用性；之后通过重力数据反演，修正了模型的误差，形成了龙门山区域基本可用的三维地质模型。

模型显示，龙门山区域的莫霍面埋深约40—66 km，分区明显，西低东高、北低南高；下地壳厚度处于9—31 km的变化范围，西北端厚，向东南端逐渐变薄；中地壳厚度变化无明显规律，存在明显的隆起与凹陷，整体厚度值在7—19 km；上地壳厚度约为22 km，自南至北呈增厚趋势，其中部分区域存在隆起与凹陷；沉积层埋深呈现东北浅、西南深的分布状态，埋深范围为1—9 km。其中中上地壳存在的隆起、凹陷与前人研究中所揭示的低速异常体吻合，地幔的变化趋势与前人反演的自东北向西南单调上升重力梯度带吻合。

本文通过建模对龙门山区域的地球物理勘探成果进行了有选择的统一整理，试图综合管理现有的成果并以三维的形式展示，是对地球物理数据统一管理与共享方式的有益探索。模型的精度受原数据精度限制，未能充分还原地壳结构的细节，建模过程忽略了一些特殊地质体，若需对特殊地质体进行建模需要高精度数据支撑；建模过程中主要存在不同剖面数据间的不一致性，受建模人员的主观影响较大，这是三维地质建模下一步需要解决的问题。