基于VP/VS波速比模型约束的张渤地震带深部电性结构研究

吴萍萍， 丁志峰* ， 谭捍东， 杨歧焱， 王兴臣

1 中国地震局地球物理研究所， 北京 100081 2 中国地质大学(北京)地球物理与信息技术学院， 北京 100083 3 河北地质大学地球科学学院， 石家庄 050031

0 引言

华北克拉通中新生代以来经历了独特的岩石圈改造和破坏过程，并伴随着广泛的岩浆活动和大规模的构造伸展盆地(翟明国等, 2005; 段永红等, 2016)，是中国大陆内部地震活动强烈地区之一(王椿镛等, 2017).地震活动分布图(图1)反映华北地区的地震活动具有明显的分带性，主要有四条重要的地震带：张家口—渤海地震带(简称张渤地震带)、河北平原地震带、汾渭地震带和鄂尔多斯西缘地震带.其中张渤地震带位于华北平原的北部，整体呈现北西向的地震活动带特征，该区域分布着一系列雁行排列的北西至北西西向断裂(图2)，该断裂带与多条NNE-NE向断裂带交汇，M6.0以上的强震主要群集于NE向断裂交汇地段.王椿镛等(2017)统计该区域的发震深度主要集中在8～25 km范围，在深度大于25km的地震数量急剧减少.深部介质的电阻率对地球内部温度、流体、熔融/半熔融等异常反映灵敏，已有研究表明地震孕育和发生与地下高导体(例如火山高导岩浆囊、壳内滑脱层、含水高导层等区域)和高阻体的接触边界息息相关(赵国泽等, 2001; 魏文博等, 2010; Zhang et al., 2016; Ye et al., 2018; 叶涛等, 2018).因此通过获取张渤地震带壳内高导和高阻异常体的分布特征，可为认识研究区的物质状态、地壳运动过程、深部孕震环境等科学问题提供深部电性结构依据.

大地电磁法以天然电磁场为场源，具有探测深度大、对低阻体反映灵敏、不受高阻屏蔽、设备轻便、成本低等优势，成为研究地球深部电性结构的主要地球物理方法之一.华北克拉通也开展了大量的大地电磁探测研究工作(魏文博等, 2006, 2008, 2010; Ouro-Djobo等, 2018; 詹艳等, 2011; 徐光晶等, 2015)，但因在张渤地震带大地电磁测点数量少，深部电阻率结构成像分辨率有限.同时大地电磁法因受限于方法本身的局限性，采集易受电磁干扰、体积效应、纵向分辨随深度增加而降低等因素影响，反演获取的电阻率模型在深部电阻率经常会出现明显的拖尾现象，对于深部异常体的底界面形态的分辨率不足(金胜等, 2010).

已有研究表明基于交叉梯度结构耦合约束(Gallardo and Meju, 2003, 2004)的大地电磁和地震资料的联合反演，可以实现两种方法的优势互补，克服单方法的局限性，提高反演精度(彭淼等, 2013; Bennington et al., 2015; Ogunbo et al., 2018; Peng et al., 2019; Wu et al., 2020).但地震数据和大地电磁数据的联合反演需要考虑两种不同数据的耦合性，在联合反演过程中需要权衡多种权重参数(例如，交叉梯度项权重因子、正则化因子、数据权重因子等)，在实际操作中存在困难.

目前华北地区已有很多地震资料的研究成果，包括人工地震深反射地壳结构(段永红等，2016)、天然地震体波层析成像(黄金莉和赵大鹏, 2005; 杨婷等, 2012; 杨歧焱等, 2018)、地震面波和噪声面波成像(鲁来玉等, 2009; 何正勤等, 2009; Fang et al., 2010; 唐有彩等, 2011; 潘佳铁等, 2011; 房立华等, 2013)、地壳厚度和泊松比(嵇少丞等, 2009; Lei, 2012)等.这些研究成果为认识该区域深部速度结构提供很好的素材.地下介质速度和电阻率模型在一些特殊构造环境(例如局部熔融半熔融异常体、含水的裂隙或断裂、壳内滑脱层等)在三维结构分布上具有相似性(Newman et al., 2008; Berdichevsky et al., 2008)，根据这一特征可以相互验证各自反演结果的可靠性，这也为本文基于地震学模型结构耦合约束的大地电磁二维反演算法开发提供思路.

弹性参数VP、VS、VP/VS、泊松比等属性的空间分布特征对揭示地球内部介质的非均匀性有重要的意义，影响弹性参数的因素主要有矿物组成、温度、压力和空隙流体等.根据已有的实验资料(陈颙等, 2009)可知，不同的地震波属性对介质温度和成分变化的反应灵敏程度不同，例如当温度发生变化时，S波速度对温度异常的响应能力比P波速度强.VP/VS波速比模型(泊松比可以通过波速比换算得到)同时兼顾了介质VP和VS属性特点，可以提供比单纯P波和S波速度参数更为丰富的信息(李志伟等, 2009).在熔融/半熔融体中，波速比值随着熔融程度的增加而增大(刘琼林等, 2011)，在这个地质环境下电导率通常也表现出类似的特征.虽然很难定量分析波速比与电导率的物性关系式，但是可以根据这两种物性在这些特殊环境中表现出来的结构相似性进行耦合约束反演，减小反演的多解性.

因此，为了避免地震资料和大地电磁资料联合反演过程中多参数的权衡，同时充分利用已有的地震学模型，本文提出将VP/VS波速比模型作为固定模型约束，通过交叉梯度结构耦合函数带入大地电磁二维反演算法中，实现基于地震学模型约束的大地电磁二维反演法，该方法即克服了联合反演的局限，同时又可以有效克服大地电磁法深部电阻率分辨率不足的缺点.文中首先介绍大地电磁二维OCCAM(Constable et al.,1987)反演算法中加入VP/VS波速比结构耦合约束策略，通过理论模型合成数据检验算法的可靠性.然后将该方法应用于张渤地震带大地电磁测深剖面，对比分析大地电磁法无约束和有约束反演结果，检验约束反演算法在张渤地震带应用的优越性和实用性.根据张渤地震带反演结果设计一套理论模型，通过理论模型合成数据试算，检验算法在该研究区应用的有效性.最后结合该区域已有的地质-地球物理资料，探讨张渤地震带深部电性结构、孕震环境及其动力学机制等问题.

1 地质构造背景

图1为华北地区大地构造背景及地震分布图，华北克拉通主要由东、中和西三个主要构造区组成.西部为稳定的鄂尔多斯高原，地震活动性主要集中在地块的边缘地区，地块内部表现为稳定性特征，在块体西缘(鄂尔多斯地震带，图1中标注为D的地震带)曾在1739年发生宁夏平罗—银川M8.0地震.中部地区为太行山造山带，位于鄂尔多斯高原和华北平原之间，该区域构造复杂，地震频发，曾发生M8.0以上地震3次(1303年山西洪洞地震、1556年陕西华县地震、1695年山西临汾地震).东部地区华北平原是研究华北克拉通岩石圈活化的热点地区(Zhu et al., 2012; 翟明国, 2019)，该区域内地震分布呈现两条共轭地震条带(图1)，为北东向活动断裂的河北平原地震带(图1中标注为B的地震带)和北西向活动断裂的张渤地震带(图1中标注为A的地震带).本文研究区为张渤地震带.

图1 华北地区大地构造背景及地震分布图A 张渤地震带；B 河北平原地震带；C 汾渭地震带；D 鄂尔多斯地震带. 红色虚线表示主要断裂分布.Fig.1 Geological background of the North China and the earthquakes distributionA Zhangbo seismic belt; B Hebei plain seismic belt; C Fenwei seismic belt; D Ordos seismic belt. Red lines donate the distribution of the main faults.

张渤地震带西起太行山，横贯燕山和华北盆地，东入渤海，是中国东部地区一条重要的北西向地震活动带，主要发育NNE和NE向断裂，部分断层走向为NWW向.通过GPS资料分析，张渤地震带以左旋走滑为主，兼有挤压运动.华北平原块体和燕山块体相对运动时张渤带左旋走滑的直接动力源(陈长云, 2016).

本文大地电磁测线沿着该地震带布设(图2)，自东向西经过唐山、三河、北京、延庆到张家口地区，主要穿过9条断裂带，依次为宁河—昌黎断裂(F1)、唐山断裂(F2)、丰田—野鸡坨断裂(F3)、夏垫断裂(F4)、黄庄—高丽营断裂(F5)、南口断裂(F6)、延庆—矾山盆地北缘断裂(F7)、怀(来)—涿(鹿)盆地北缘断裂(F8)和张家口断裂(F9).其中在黄庄—高丽营断裂(F5)是一条京西山前大断裂，是华北平原和京西隆起的分界线，走向NE20°～50°，总长约110 km，断裂面倾向南东，倾角55°～75°.

图2 张渤地震带构造背景和MT测点分布图F1 宁河—昌黎断裂; F2 唐山断裂; F3 丰台—野鸡坨断裂; F4 夏垫断裂; F5 黄庄—高丽营断裂; F6 南口断裂; F7 延庆—矾山盆地北缘断裂; F8 怀(来)—涿(鹿)盆地北缘断裂; F9 张家口断裂. 红色五角星: 三河平谷M 8.0地震; 蓝色圆M 7.0～M 8.0之间地震; 黄色圆M 6.0～M 7.0之间地震; 绿色圆表示M 5.0～M6.0之间地震. 01～31数字为大地电磁测点编号.Fig.2 Geological settings of Zhangbo seismic belt and the locations of the MT sitesF1 Ninghe-Changli fault; F2 Tangshan fault; F3 Fengtai-Yejituo fault; F4 Xiadian fault; F5 Huangzhuang-Gaoliying fault; F6 Nankou fault; F7 North margin fault of Yanqing-Fanshan basin; F8 North margin fault of Huailai-Zhoulu basin; F9 Zhangjiakou fault. Red pentagram: Shanhe-Pinggu M8.0 earthquake; blue circles: M 7.0～M 8.0 earthquakes; yellow circles: M6.0～M7.0 earthquakes; green circles: M5.0 ～M6.0 eaqthquakes. The numbers from 01 to 31 are the number of the MT sites.

测线最东边为由宁河—昌黎断裂(F1)、唐山断裂(F2)、丰台—野鸡坨断裂(F3)组成的唐山强震区断裂构造.唐山断裂带是1976年唐山7.8级地震的发震构造，人工地震资料展示唐山断裂带走向NE30°，在地壳伴生复杂断裂构造带，表现为花状构造断裂特征，主干断裂沿伸到Moho面(刘保金等, 2011a)，属走滑断层性质.夏垫断裂(F4)是一条NNE向岩石圈尺度的区域性深断裂带，全长28 km，总体走向45°，倾向SE，倾角70～80°，是北京平原历史上唯一发生过8级地震的活动断裂.

测线最西侧为张家口断裂带(F9)，属于首都圈地震重点监视防御区，为张家口-渤海地震构造带与山西地震构造带的交汇区(赵博等, 2011)，由北西西向和近东西向两组、 多段断裂组成， 北西西向断层为断裂主体，倾向南东，倾角60°左右，全长约30 km，属正断层性质.怀涿盆地北缘断裂(F8)控制着怀来—涿鹿盆地的发育，总体走向北东，倾向南东，倾角50～70°，是重要的控盆活动正断裂之一.延庆—矾山盆地北缘断裂(F7)是一条控制延庆—矾山盆地北缘的主边界正断裂带，总体走向NE，倾向SE，倾角50～80°，全长约102 km.

2 基于VP/VS波速比模型约束的二维大地电磁反演方法

2.1 交叉梯度项

Gallardo和Meju(2003, 2004)通过定义两种不同物性参数的梯度的叉乘来定量描述两种物性的结构相似性，这是本文实现基于VP/VS波速比模型约束的大地电磁二维反演算法的关键.根据定义可将电阻率和波速比模型的交叉梯度函数表达式写成：

(1)

mr表示电阻率模型，ms表示VP/VS模型.根据求解旋度公式可将交叉梯度项分解为

(2)

(3)

(4)

ty=Wcg,s·mr.

(5)

Wcg,s为与网格和ms相关的常数矩阵. 当电阻率反演初始模型mr0为均匀半空间模型时，则Wcg,s·mr0=0，那么公式(5)可以改写成：

ty=Wcg,s·(mr-mr0).

(6)

取公式(6)的平方和构建VP/VS波速比模型的结构耦合约束项，则公式可写成：

(7)

2.2 基于VP/VS波速比模型约束的二维大地电磁反演目标函数

在大地电磁无约束反演目标函数中，通过权重因子η将VP/VS波速比模型的结构耦合约束项加入到无约束反演目标函数中，形成基于VP/VS波速比结构耦合约束的二维大地电磁反演目标函数，具体表示如下：

(8)

f(·)是二维大地电磁正演响应函数，λ为拉格朗日因子，Cd为数据协方差矩阵，Wm为模型光滑矩阵，dobs为观测数据矩阵.对目标函数φ求梯度，则其梯度表达式可以写成：

(9)

构建

(10)

根据OCCAM反演算法(Constable et al., 1987)的迭代公式，即可获得基于VP/VS波速比模型约束反演第k次模型更新公式：

(11)

公式中λ因子的选取采用自动搜索算法实现，每一次迭代都要求根据RMS值选取最优λ值，具体算法可参考Constable等(1987).图3为本文开发的基于VP/VS波速比模型约束的大地电磁二维反演算法流程图.

图3 基于VP/VS波速比模型约束的大地电磁二维反演流程图Fig.3 The flow chart of 2-D MT inversion algorithm with the VP/VS wave velocity ratio model constraints

3 理论模型合成数据反演试算

3.1 理论模型设计

为了检验基于VP/VS波速比模型约束的大地电磁反演算法的可靠性，本文设计一套复杂电阻率模型(图4a)，其中包含有A、B、C1、C2、F电阻率块体、近垂直构造体D和上涌构造体E.图4a中黑色三角形为大地电磁测点，基于有限单元法(Constable et al.,1987)正演计算10-3Hz～1103Hz 区间42个频点的TE和TM模式的视电阻率和相位，分别加入5%的高斯随机误差，作为观测数据.

图4b为采用的VP/VS波速比固定模型，波速比值范围为1.5～1.9.对比电阻率模型，波速比模型中设计了与电阻率模型结构相似(A、B、C1、F、D、E)和结构不相似(C2)的异常块体或构造体.将VP/VS波速比模型(图4b)作为固定模型约束，分别采用无约束和有约束大地电磁反演算法获得电阻率结构，对比分析这两种方法的反演结果，检验约束反演算法的可靠性和优越性.

图4 理论电阻率模型(a)和VP/VS波速比模型(b)(a) 理论电阻率模型; (b) 固定约束的VP/VS波速比地震学模型. A、B、C1、C2、F表示不同的物性块体；D为近垂直构造体；E为上涌构造体.Fig.4 The resistivity model (a) and the VP/VS wave velocity ratio model (b)(a) Synthetic resistivity model; (b) The fixed VP/VS wave ratio model. A, B, C1, C2 and F represent the different physical blocks; D represents the near vertical structure; E represents the upwelling structure.

3.2 理论模型反演结果

图5a、b分别为无约束和有约束反演结果.从图5a中可以发现大地电磁法在无约束反演情况下，浅部电阻率结构具有较高的分辨率，例如块体A的异常值和异常边界都恢复得与真实模型相似、C1和C2块体的上边界接近于真实模型边界、块体B的电阻率值在浅部恢复得与真实值相近.但无约束反演结果在深部分辨率明显不足，具体表现在：异常体C1和C2底边界缺失，形成了与F贯通的两个上涌构造体；近垂直构造体D在浅部有明显的边界信息，随着深度增加，构造体的边界信息减弱；上涌构造体E几乎没有体现出来.

图5 理论模型无约束和有约束的大地电磁二维反演结果(a) 无约束大地电磁反演结果; (b) 有约束大地电磁反演结果; (c) 无约束反演电阻率模型与VP/VS波速比模型的交叉梯度值空间分布; (d) 有约束反演电阻率模型与VP/VS波速比模型的交叉梯度值空间分布. RMS为数据拟合差.Fig.5 The results derived from the unconstrained and constrained inversion method(a) The resistivity model derived from the unconstrained inversion method; (b) The resistivity model derived from the constrained inversion method; (c) The spatial distribution of the cross-gradient value between the VP/VS ratio model and the unconstrained inverted resistivity model; (d) The spatial distribution of the cross-gradient value between the VP/VS ratio model and the constrained inverted resistivity model. RMS represents the data residual.

图5b为大地电磁法约束反演结果，从图中可以看出电阻率反演结果在深度分辨率有明显提高.在电阻率模型和VP/VS波速比模型结构相似处电阻率边界信息有明显改善，具体表现在：C1和C2块体由原来的连通状态恢复成两个独立圈闭的块体，其中C1的边界信息与VP/VS波速比模型相似，在C1的异常值和边界都恢复得与真实模型相似，而在C2位置处，因在VP/VS波速比模型中没有存在C2异常，虽然C2异常值有改善，但边界信息与真实情况仍有差距；上涌构造E和近垂直构造D的边界信息恢复得与真实模型相近.图5c和d分别为无约束反演结果、有约束反演结果与VP/VS波速比模型的交叉梯度值空间分布图，从图中可以发现在块体C1周边、近垂直构造D、上涌构造E的交叉梯度值有明显降低，这说明经过约束反演后的模型，在这几个区域的结构更相似.从反演结果中不难发现对于浅部电阻率值和形态(例如块体A)在无约束和有约束反演都可以很好地恢复，这说明基于地震学模型约束的大地电磁反演算法可以在满足数据拟合差的前提下，获取电阻率模型和地震学模型结构上更为相似的模型，有效提高大地电磁法深部电阻率分辨率.

4 张渤地震带大地电磁测深剖面

4.1 大地电磁数据采集和处理

图2为张渤带大地电磁测点布设图，测线跨越张家口、北京、唐山等城市，全长400 km，测线方向为北偏西54°.测点布设采用“+”字形正南北东西向布极方式，采集加拿大凤凰公司的MTU-5A宽频带大地电磁观测设备.从该测线中挑选出数据质量较高的31个MT测点，分别对每个测点数据进行滤波、重采样、频谱分析等预处理，基于WinGLink平台获取0.001 Hz到320.0 Hz的TE模式和TM模式的视电阻率和相位曲线.同时采用Berdichevskiy和Rhoplus方法检验视电阻率和相位曲线的一致性，将处理后的数据作为反演数据.

图6中黑色曲线为TE模式的视电阻率和相位，红色曲线为TM模式视电阻率和相位.根据地质构造特征，将测线上MT测点数据分成四段，分别为唐山断裂带区域(图6a，红色测点)、三河—平谷断裂带区域(图6a，黄色测点)、华北平原—京西丘陵过渡区(图6a，蓝色测点)和京西丘陵区域(图6a，黑色测点).从图中可以看出除了华北平原—京西丘陵过渡带的视电阻率和相位曲线形态差异比较大之外，其他三个区域的数据形态较为一致.其中唐山断裂带区域视电阻率曲线形态和京西丘陵区域曲线形态相似，但京西丘陵段的最高值出现在0.15 Hz左右，而唐山断裂带区域最高值出现在0.06 Hz左右，这两段的相位曲线差异比较大，说明地下电性结构山区和平原存在差别.三河—平谷断裂带视电阻率值比较低，最高值约300 Ωm，该区域的相位曲线形态与唐山断裂带相似.而华北平原—京西丘陵过渡带的视电阻率和相位曲线形态复杂，没有明显的规律，这也间接说明该区域地下电性结构复杂的特征.

图6 实测大地电磁资料视电阻率和相位曲线(a) 红色三角形为Segment 1测点位置; 黄色三角形为Segment 2测点位置; 蓝色三角形为Segment 3测点位置; 黑色三角形为Segment 4测点位置. 黑色曲线为TE模式的视电阻率和相位; 红色曲线为TM模式的视电阻率和相位.Fig.6 Measured MT apparent resistivity and phase data(a) Red triangles represent the MT sites corresponding to the data of segment 1; Yellow triangles represent the MT sites corresponding to the data of segment 2; Blue triangles represent the MT sites corresponding to the data of segment 3; Black triangles represent the MT sites corresponding to the data of segment 4. Black curves represent the apparent resistivity and phase data of TE mode; red curves represent the apparent resistivity and phase data of TM mode.

4.2 VP/VS波速比模型

杨歧焱等(2018)基于华北地区数字地震台网176个固定台站，采用近震体波层析成像法获取了张渤地震带VP/VS波速比地震学模型.本文提取了张渤地震带对应大地电磁测线剖面的VP/VS波速比模型(图7)，对该剖面进行平滑，去除模型的突变点，然后将该模型作为固定约束代入公式(5—11)中的Wcg,s，实现基于VP/VS波速比模型约束的大地电磁法二维反演.

图7 张渤地震带VP/VS波速比剖面图TS 唐山； SH 三河； BJ 北京； ZJK 张家口.Fig.7 The VP/VS wave velocity ratio profile along Zhangbo seismic beltTS Tangshan; SH Shanhe; BJ Beijing; ZJK Zhangjiakou.

4.3 反演网格

二维大地电磁反演网格剖分中，用均匀和非均匀剖分法分别对目标研究区域和边界区域进行网格化处理(图8).在目标研究区采用均匀网格剖分，X轴网格间距为3 km，Z轴网格间距为2 km.在边界网格采用的是指数递增网格间距.图中8红色框为VP/VS波速比模型约束区，在反演过程中只在红色框区域内进行结构耦合约束.

图8 大地电磁正、反演网格剖分黑色线为网格；黑色方框为约束反演区域.Fig.8 Mesh discretization for MT forward and inversionBlack lines represent the mesh; black rectangle shows the constrained inversion region.

4.4 张渤地震带反演结果

图9a为无约束和有约束的大地电磁法反演迭代RMS曲线图，从图中可知这两种方法都具有较好的收敛性.从RMS迭代变化趋势可以看出，当加入VP/VS波速比模型约束项时，RMS值迭代下降的速度变慢，但经过多次迭代后，有约束和无约束反演结果的RMS都收敛到相同水平上.

图9 RMS迭代曲线和约束项权重因子L-Curve曲线(a) 无约束和有约束反演RMS迭代曲线; 空心圈曲线为无约束反演(SP)迭代RMS曲线; 实心点曲线为有约束反演(CG)迭代RMS曲线. (b) 约束项权重因子L-Curve曲线; 实心圈曲线为不同权重因子的RMS值; 黑色箭头标注为本文选取的权重因子.Fig.9 Plots of RMS iteration curves and L-curve analysis for selecting constrained weight factors(a) RMS iteration curves for no-constrained and constrained inversion. The hollow circle curve represents the RMS iteration curve for unconstrained inversion. The solid dot curve represents the RMS iteration curve for constrained inversion. (b) L-curve anasys for selcting constrained weight factors; The solid dot curve represents the RMS value for different weighting factors; the black arrow represnets the optimal weighting value.

图9b为公式中VP/VS波速比模型约束项权重因子(公式8中η)的选择.本文采用“L-Curve”方法，选取了多个权重因子进行二维大地电磁法约束反演(图9b)，图中黑色箭头为本文最终确定的权重因子.图10展示了6个测点无约束和有约束反演模型响应曲线，从图中可以发现无约束和有约束反演模型响应曲线基本上重合，与观测数据有较好的一致性，说明这两个模型都可以很好地拟合实测数据.

图10 无约束和有约束反演模型在测点01、09、16、19、26和30处的正演响应红色点: 实测数据(real); 黑色圈: 无约束反演结果正演响应(SP); 蓝色叉: 有约束反演结果正演响应(CG).Fig.10 The synthetic data sets of the unconstrained and constrained inverted model at the sites of 01, 09, 16, 26, and 30.Red dots represent the measured data (real); black circles represent the synthetic data sets of the unconstrained inverted model (SP); blue crosses represent the synthetic data sets of the constrained inverted model (CG).

图11展示了大地电磁法二维无约束和有约束反演结果，图11a为测线对应的地形曲线，可以发现测线的西部测点主要位于京西丘陵区，测线的东部主要位于华北平原地区，测线上由东向西穿越了唐山断裂带、三河—平谷断裂带、华北平原—京西丘陵过渡带、延怀盆地和张家口断裂带.

图11 大地电磁法无约束和有约束反演结果(a) 地形剖面图; (b) 无约束反演电阻率模型; (c) 有约束反演电阻率模型. H1—H5为高阻异常体编号; L1—L3为低阻异常体编号.Fig.11 The resistivity models derived from the unconstrained and constrained inversion method(a) Topographic profile; (b) The resistivity model derived from the unconstrained inversion method; (c) The resistivity model derived from the constrained inversion method. H1—H5 represent the numbered high resistivity anomalies; L1—L3 represent the numbered low resistivity anomalies.

为了方便对比这两种方法反演结果，将图中5个高阻体标注为H1、H2、H3、H4和H5，3个低阻区域标注为L1、L2和L3.从图中可以看出，基于VP/VS波速比模型约束的反演结果(图11c)的横向分辨率和纵向分辨率都有明显的提高，主要表现有：在唐山强震区的高阻异常体H5沿伸到下地壳，经过约束反演后，高阻异常的拖尾现象有明显的改善，底部边界约在20 km左右；从夏垫断裂带到京西山前断裂带的高阻体异常体H4近垂直地沿伸到地幔，很难识别出该异常体在地壳内的异常形态，经约束反演后，深部的分辨率有明显的提升，高阻异常体H4在地壳中表现为圈闭异常形态；在京西丘陵区域的地壳高阻异常体H1、H2和H3在无约束反演中表现为整体高阻异常，深度延伸到下地壳，经约束反演后，可以明显看出H1、H2和H3这三个高阻异常体是独立分开的三个高阻异常，同时这三个异常体的底界面大约在20 km；经过约束反演后，发现深部高导区域L1和L2从30～40 km上涌，这个现象在无约束反演结果中没有发现；在夏垫断裂带区域无约束和有约束反演都可以发现明显的低阻异常体L3，但低阻异常的形态、规模有很大的差异.对比VP/VS波速比模型(图7)和大地电磁法约束反演结果(图11c)，在20 km范围内，电阻率模型与VP/VS波速比模型的异常体表现出明显的相似性.而在大于20 km深度，这两个模型存在比较大的差异，表现在约30 km处波速比模型出现低波速比异常层，而电阻率模型上未发现该特征.

图12为大地电磁法无约束和有约束反演的电阻率模型与VP/VS波速比模型的交叉梯度值空间分布图，mean|ty|为约束反演区域(图8黑色框)中每个网格的交叉梯度值的绝对值的平均值.从图中可以发现有约束反演结果与VP/VS波速比模型的交叉梯度值比无约束反演结果的交叉梯度值有明显下降(mean|ty|由2.94降到1.5)，这说明基于VP/VS波速比模型约束的大地电磁反演可以在满足数据拟合差的前提下，可以获取与地震学模型结构上更为耦合的电阻率模型.

图12 无约束、有约束反演电阻率模型与VP/VS波速比模型的交叉梯度值分布图(a) 无约束反演电阻率模型与VP/VS波速比模型的交叉梯度值空间分布; (b) 有约束反演电阻率模型与VP/VS波速比模型的交叉梯度值空间分布. mean|ty|为约束反演区域(图8红色框)中每个网格的交叉梯度值的绝对值的平均值.Fig.12 The spatial distribution of the cross-gradient value between the VP/VS ratio model and the non-constrained or the constrained inverted resistivity model.(a) The spatial distribution of the cross-gradient value between the VP/VS ratio model and the unconstrained inverted resistivity model; (b) The spatial distribution of the cross-gradient value between the VP/VS ratio model and the constrained inverted resistivity model. mean |ty| represents the average of the absolute value of the cross-gradient values in the common region (red box in Fig.8).

根据上述现象，不难看出要使得基于地震学模型约束反演算法能有效改善大地电磁反演结果的前提条件是地震学模型和电阻率模型存在结构相似性.当两种物性结构相似时，基于交叉梯度的约束反演可以有效地加强两种物性结构耦合关系，提高结构相似处的分辨率.在结构不相似的区域，交叉梯度项几乎不起作用，保留单方法反演的优势.Wu 等(2020)、Peng 等(2019)通过大量的理论模型算例充分地论证这一特征.

4.5 张渤地震带约束反演结果可靠性和有效性检验

通过对比无约束、有约束的大地电磁反演结果，可知在地震学模型的约束下，张渤地震带深部电性结构在横向和纵向的分辨率都有明显的提升.为了检验基于VP/VS波速比模型约束的大地电磁反演算法在张渤地震带应用的可靠性和有效性，本文考虑设计形态与反演结果相似的5个独立的高阻体(图11中H1—H5)和3个上涌的低阻构造体(图11中L1—L3)，通过理论模型合成数据试算，采用与实测数据反演一样的VP/VS波速比模型(图7)进行约束反演，对比无约束和有约束反演结果，从理论模型上检验该算法是否能提高张渤地震带深部电阻率横向和纵向分辨率.

为了能满足上述5个独立的高阻体和3个上涌的高导构造体的理论模型需求，本文将约束反演结果(图11c)中电阻率值大于103.5Ωm的电阻率值设置为恒定值103Ωm，将小于101.8Ωm的电阻率值设置为恒值10 Ωm，其他区域的值设置为102Ωm.图13a为本文设计的理论电阻率模型，从图中可以看出理论模型在20km以内设计有5个圈闭的高阻异常体(图13a 中H1—H5)，中下地壳3个低阻隆起(图13a 中L1—L3).

测点位置(图13a黑色三角形)和每个测点的频段选择与实际观测一样，正反演网格也跟上述一样.正演每个测点的理论视电阻率和相位曲线，同时加入与实测数据相同的误差，形成观测数据.

图13 理论电阻率模型(a)、无约束(b)和有约束(c)反演结果(a) 理论电阻率模型; (b) 无约束反演电阻率模型; (c) 有约束反演电阻率模型. H1—H5为高阻异常体编号; L1—L3为低阻异常体编号.Fig.13 Synthetic model (a), unconstrained (b) and constrained (c) inverted models(a) Synthetic model; (b) Unconstrained inverted model; (c) Constrained inverted model. H1—H5 represent the numbered high resistivity anomalies; L1—L3 represent the numbered low resistivity anomalies.

图13b和c分别为大地电磁法二维无约束和有约束的反演结果，对比这两种方法的反演结果，可以发现无约束反演的结果与上述实际数据无约束反演结果出现的问题相似，主要表现在：独立高阻异常体H1—H4在无约束反演(图13b)中表现为连接在一起、规模大的高阻异常；低阻体L1和L2几乎没有体现出来，L3的低阻异常体形态与真实异常体形态差别很大；同时异常体的底界面相比真实模型有明显的往下拖尾的现象.从无约束反演结果可以看出因大地电磁法受限于方法本身体积效应、深度分辨率随深度增加而降低等局限，无约束反演获得的电阻率模型在横向和垂向分辨率有限.

图13c为基于VP/VS波速比模型约束反演的结果，从图中可以明显的看出浅部高阻异常体(H1—H5)的横向分辨率和纵向分辨率都有明显的提高，5个高阻异常体可以清晰的展示出独立异常体的特征，异常体的拖尾现象得到明显的改善；对于低阻体L1、L2和L3异常体的反演精度有明显的提高.因此，通过这个理论模型合成数据试算，验证了基于VP/VS波速比模型约束的二维大地电磁反演法在张渤地震带应用可以提高深部电阻率的横向和纵向分辨率，从而验证了该算法在张渤地震带应用的可靠性和有效性.

5 讨论

根据反演得到的张渤地震带地壳电性结构，结合地质资料绘制电性构造解释图(断裂信息参考徐锡伟等(2002)).图14中断裂F1—F9的深部形态分别用黑色虚线展示，断裂地表位置与图2中断裂(F1—F9)位置对应.从整体上来看，电阻率结构表现出明显的横向不均匀性，与地表的各个构造相呼应.对比图11a地形曲线和构造分区，可以发现测线上断裂带在电性结构上表现为明显的电性分界面或电性梯度带，在张家口断裂带、延庆—怀柔盆地、夏垫断裂带和唐山断裂带在中上地壳普遍存在高阻异常体，而在中下地壳表现为低阻异常特征.从M≥3.0地震分布特征来看，发现该区域的地震活动与断裂带的分布息息相关，具有明显的分段性，以下将分段讨论剖面电性结构特征.

唐山震群主要有三条断裂带控制(宁河—昌黎断裂F1、唐山断裂F2和丰台—野鸡坨断裂F3).唐山7.8 级地震发生在主断裂带(F2)上，该区域为高阻块体与低阻块体的相交地带、偏向高阻体一侧.根据人工地震反射资料可知唐山断裂为主断裂，延伸至Moho面，地壳伴生的断裂构造表现为走滑断裂(花状构造)特征.图中可以发现高阻体的底界面大约在20 km左右，地震的发生主要集中在20 km以内，位于高阻区域或高低阻交界区.人工深反射剖面(刘保金等, 2011a)和天然地震观测(刘启元等, 2007; 杨歧焱等, 2018)发现唐山震区在深度<20 km内呈现高速异常，在约深20 km处存在明显的低速层(体)，该低速层(体)一直延伸到40 km以下的深度.刘国栋等(1983)发现唐山地震源区下约20 km深处存在高导异常.嘉世旭和张先康(2005)通过深地震测深资料也发现唐山强震区中下地壳存在低速异常.岩浆底侵和置换作用是下地壳和地幔物质交换、相互作用的重要形式，推测唐山震区中下地壳的低速低阻异常体与流体和热物质作用有关，大规模的幔源物质侵入，形成了中上地壳的孕震环境.

在三河—平谷断裂带区域，存在明显的壳浅部到深部低阻特征，中下地壳存在上涌低阻体(L3)，浅部和深部的高导体不连续，地震分布主要位于该不连续层中，说明该区域深浅高导异常的成因可能不同.邓前辉等(2001)通过大地电磁法也揭示夏垫断裂带中下地壳高导体和地幔上涌现象，浅部和深部的构造存在差异，推测震源体的形成和发震与深部高角度隐伏断裂有关.对应地质资料可知，三河—平谷断裂带为第四纪活动隐伏高倾角断裂带，地表覆盖有上新统-第四系沉积层，浅部低阻异常与沉积层对应.从人工地震深反射剖面(刘保金等, 2011b; 赵金仁等, 2004)可知该断裂带是一条深浅共存的断裂构造带，在Moho存在复杂的楔形反射带，说明该区域地壳发生过强烈的挤压、变形，同时也反映出岩浆活动对地下地壳结构经过强烈改造，形成了复杂的地壳深部结构.徐锡伟等(2002)也从人工地震反射剖面发现Moho面存在反射过渡带，推测是早期地壳伸展过程中形成断裂带，填充了来自上地幔的基性岩墙.从图14中深部上涌低阻异常(L3)可能反映的是上地幔热物质侵入，Zhu等(2012, 2015)指出西太平洋板块的俯冲时，俯冲板块的脱水交代上覆地幔楔使其发生熔融，俯冲带后撤后导致岩石圈强烈伸展，导致地幔热物质上涌，形成了中上地壳的孕震环境.

黄庄—高丽营断裂带(F5)和南口断裂(F6)为京西丘陵的山前断裂带，是大型拆离断裂带，在电阻率结构上表现为明显的电阻率变化梯度带.在山前断裂带到夏垫断裂带之间浅部有几公里的反映沉积层特征的低阻层，在低阻层下方、延伸到下地壳为规模较大的高阻异常，未见低速体上涌现象.詹艳等(2011)从石家庄地区大地电磁剖面也发现太行山山前断裂从浅地表到下地壳表现为高电阻体，无壳内高导层发育.京西丘陵的山前断裂带是华北克拉通东部与中部的边界，在重力梯度带、地形、岩石圈厚度等资料显示在该区域附近变化显著(Zhu et al., 2012).规模大、延伸深的高阻异常(H4)可能是华北克拉通东部和中部地区构造运动形成的高密度的辉长岩、铁镁质岩石或者下地壳麻粒岩(翟明国等, 2005; 翟明国, 2019).华北平原地区强震主要发生在NE和NW两组构造带，黄庄—高丽营断裂带(F5)和南口断裂(F6)山前断裂带自新生代以来对北京凹陷起着重要的控制作用(焦青等, 2005;焦青和邱泽华, 2006)，是北京地区值得关注的强震发生重点监视防御区.

在南口断裂(F6)以西区域为京西丘陵地区，主要有延庆—怀来盆地和张家口断裂带.延怀盆地是由山西裂谷走滑运动形成的NE向次级张性构造区，地表延庆—矾山次级盆地北缘断裂(F7)和涿鹿—怀来盆地北缘断裂(F8)控制盆地构造特征.在电阻率结构成像中(图14)，在浅部表现明显的低阻异常，在中下地壳有明显的低阻体沿着盆地下方通道侵入的现象.从地震反射剖面(张先康等, 1996)可知在延庆—矾山盆地北缘断裂正下方存在低速体，其成因可能与该区域地壳伸展运动、幔源物质上涌有关，在地表形成盆岭构造.

图14 张渤地震带电性结构解释图白色圈为2007年至今M≥3.0地震(由中国地震台网中心提供). F1 宁河—昌黎断裂; F2 唐山断裂; F3 丰台—野鸡坨断裂; F4 夏垫断裂; F5 黄庄—高丽营断裂; F6 南口断裂; F7 延庆—矾山盆地北缘断裂; F8 怀(来)—涿(鹿)盆地北缘断裂; F9 张家口断裂. F1—F8走向为NE或NEE；F9走向为NW.Fig.14 The geology interpretation maps of the electrical structure beneath Zhangbo seismic beltWhite circles: M≥3.0 earthquakes since 2007 (provided by China earthquake Networks Center). F1 Ninghe-Changli fault; F2 Tangshan fault; F3 Fengtai-Yejituo fault; F4 Xiadian fault; F5 Huangzhuang-Gaoliying fault; F6 Nankou fault; F7 North margin fault of Yanqing-Fanshan basin; F8 North margin fault of Huailai-Zhoulu basin; F9 Zhangjiakou fault. The strike of F1—F8 is NE or NEE; the strike of F9 is NW.

张家口断裂带位于张渤地震带和山西地震带的交汇区(赵博等, 2011)，区域内断裂相交和切割，形成复杂的孕震构造格局，主要由北西西向和近东西向两组、 多段断裂组成， 北西西向断裂为断裂带主体(图14中F9)，与研究区数条北北东至北东向的断裂显著不同(高战武等, 2001)，是首都圈地震重点监视防御区.张学民等(2012)通过地震体波层析成像发现太行山北段多构造交界的张家口地区存在局部低速异常，推断该区域受到西太平洋板块俯冲的影响，表现为拉张应力场，深部地幔热物质活动及浅部应力累积引发活跃的构造活动.地质资料显示张家口地区出露中生代形成的源于古老下地壳或地幔的碱性花岗岩(杨进辉等, 2006; 韦忠良等, 2008; 汪洋和程素华, 2010)，从岩石学、地球化学角度论证了该区域经历了大规模的壳幔相互作用.从本文获取的电性结构图上可以发现在该断裂带上深部有低阻体上涌，在张家口中下地壳高导层的规模比怀来盆地大，这可能暗示张家口断裂带的深部物质作用更为强烈.

6 结论

本文实现了基于VP/VS波速比模型约束的大地电磁二维反演算法，通过理论模型合成数据试算，检验了算法的可靠性.将约束反演算法应用于张渤地震带31个大地电磁实测资料，研究发现唐山断裂带中上地壳表现为高阻异常，地震主要分布于高阻异常体周边，在下地壳底部有上涌的低阻异常体，推测该区域深部低阻异常体与幔源物质侵入有关.三河—平谷断裂带浅部表现为沉积层的低阻异常特征，深部高导异常与夏垫断裂带展布、壳幔物质相互作用有关，地震事件主要发生在浅部和深部高导异常不连续区域.太行山山前断裂电阻率结构上表现为明显的电阻率变化梯度带，在与华北平原过渡带中下地壳存在规模大、延伸深的高阻异常，推测与该区域因构造运动形成高密度的辉长岩、铁镁质岩石或者下地壳麻粒岩有关.延怀盆地在浅部表现明显的低阻异常，在中下地壳有明显的低阻体沿着盆地下方通道侵入的现象，其成因可能是地壳伸展运动，幔源物质上涌，形成地表的盆岭构造.张家口断裂带中下地壳高导异常层比怀来盆地深部高导异常规模大，这可能暗示张家口断裂带的深部物质作用更为强烈.

致谢感谢中国地质大学(北京)王瑜教授和中国地震局地球物理研究所尤惠川研究员对本文的地质构造解释提供的帮助.感谢审稿专家和编辑部对本文提出的重要修改意见.