首都圈极低频电磁台网区地下电性结构探测

董泽义 汤 吉* 赵国泽 陈小斌 崔腾发 韩 冰 姜 峰 王立凤

1)中国地震局地质研究所，地震动力学国家重点实验室，北京 100029

2)应急管理部国家自然灾害防治研究院，北京 100085

3)中国地震局地震预测研究所，北京 100036

4)中国科学院南海海洋研究所，边缘海与大洋地质重点实验室，广州 510301

0 引言

地震是对人类生命与财产安全造成严重危害的突发性自然灾害之一。20世纪以来，中国共发生6级以上地震800多次，死于地震的人数达55万之多，占全球因地震死亡人数的50%以上。显然，开展地震预测预报研究既是科学发展的需要，也是人类社会可持续发展的需要(陈颙等，2003；陈运泰，2009)。

由于地震发生的机制十分复杂，地震预测预报问题迄今尚未得到解决，仍然是世界重要的前沿科学技术难题之一。尽管如此，世界各国特别是地震灾害严重的国家，如中国、美国、日本、意大利、希腊、俄罗斯等针对地震预测预报投入了大量人力物力，并开展了长期的探索研究，力求减轻地震灾害(赵国泽等，2012)。经过上百年的努力探索，全球建立了大量地震监测台站，积累了大量观测资料，也取得了一些重要进展。大量观测资料表明，强震前的电磁异常现象确实是存在的，在众多地震预报预测研究手段中，电磁观测被认为是最有可能首先取得地震预测突破性进展的方法之一(Bleieretal.，2005；Zhaoetal.，2009；Uyeda，2015)，其已成为中国地震短临预测的最有效手段之一。

中国从1966年邢台地震开始对地震电磁观测开展研究，至今已建设了数百个地电场、地磁场和地电阻率观测台站，得到了不少地震前后电、磁场和电阻率对异常变化信息(赵玉林等，1978；黄清华等，2006；高曙德等，2010)。但这些传统的电磁法观测手段主要关注恒定场成分，采样率较低(最高采样频率为1Hz)，未涉及对地震事件敏感的超低频、极低频(SLF、ELF，3～300Hz)及其附近频带的电磁场观测(赵国泽等，2012，2015)。相比传统的电法或磁法，交变电磁场法在地震观测中具有自身优势与潜力，也被应用于地震电磁异常的观测研究，在国内外取得了一些重要的研究成果(Reddyetal.，1976；张云琳等，1994；汤吉等，1998)。但大多数观测关注的是天然源电磁场，由于受到日益严重的人文干扰，其应用效果也受到了限制。为了提高对各种电磁异常现象的观测识别和捕捉能力，建立中国地震电磁观测新技术系统，人工源极低频电磁技术CSELF(Control Source Extremely Low Frequency)被提出和发展起来，经过约十年的研究和试验证明，CSELF法是当前最有发展前途的地震短临前兆监测方法之一，它不仅具有抗干扰能力强、信噪比高、覆盖范围大、探测深度大、有利于识别和捕捉地震异常信号等许多优点，同时其既可用于研究空间电磁场异常，也可用于研究地壳电性结构异常，使观测成本极大降低(赵国泽等，2003a，b，2007，2010a，2012)。

“十一五”国家重大科学技术设施工程项目“极低频探地(WEM)工程”地震预测分系统在地震活动较强烈的南北地震构造带南段的川滇地区和地震重点监视区首都圈地区分别建设了15个CSELF电磁台站，组成了中国乃至世界的第1个既观测天然源电磁场(频带为1000～0.001Hz)也观测人工源电磁场(频带为300～0.1Hz)的大陆CSELF电磁台网(图1a)。该台网既能监测地震引起的空间电磁场异常，也能监测震源区附近地下电性结构的变化，实现了真正的四维立体电磁异常信息的动态监测，提升了川滇与首都圈地区中强以上地震电磁异常的捕捉能力，同时也提高了地震预测预报水平，推进了地震预警立体监测系统的建设。

图1 首都圈极低频电磁台网区的电性结构测量Fig.1 Investigation of electrical resistivity structures beneath the CSELF network in the capital circle region.a 中国大陆极低频电磁台网；b 首都圈CSELF台网区的地质构造简图；c 大同台结构测量测点分布

了解CSELF台网区和各台站地下的电性结构背景，是研究相应地区地震预测监测问题的基础，这不仅对于认识川滇与首都圈2个台网区的深部结构特点具有重要意义，而且为识别地震活动前、后地下电性结构可能发生的动态变化、理解地震可能发生破裂的程度和动态过程、探索地震及地震异常现象可能的成因机制等提供了重要数据。该项目首次在川滇和首都圈2个台网区开展了深部电性结构探测，通过数据采集、处理与分析、反演获得了2个台网区各台站及附近区域的地下电性结构。本文将主要介绍如何开展首都圈CSELF台网区地下电性结构探测，包括野外数据采集、数据的维性分析、数据反演以及台网区的深部电性结构特征分析等方面。文中涉及的图件较多，由于篇幅有限，只能以部分台站为例给出示范图件，特此说明。

1 首都圈CSELF台网的基本情况与电性结构测量

1.1 台网的组成与地质构造背景

中国大陆CSELF台网依托中国地震局已有30个地震台站进行建设。从地理上看，首都圈CSELF台网基本覆盖了华北中、东部地区，由天津市宝坻台，河北省丰宁台、怀来台、文安台、涉县台，山东省无棣台、莱阳台、安丘台、莒县台，山西省大同台、代县台、夏县台，辽宁省大连台、营口台，内蒙古呼和浩特台共15个台站组成(图1a，b)。

中国华北是地球上最古老、形成于太古宙时期的克拉通块体之一(Liuetal.，1992)，它可分为东部华北盆地和西部鄂尔多斯地块2部分，中间由山西断陷带隔开，其南、北部分别为秦岭-大别造山带和阴山-燕山造山带，2条巨型线性构造(太行山重力梯度带与郯庐断裂带)穿其而过(图1b)。已有研究表明，华北东部地区自显生宙以来经历了独特的岩石圈改造和破坏过程，并伴随着大规模的构造伸展、断陷盆地的形成和广泛的岩浆活动(Griffinetal.，1998；Fanetal.，2000；Xu，2001；Gaoetal.，2004)。华北地区拥有独特的板内地震活动带，历史上曾发生多次7级以上破坏性大地震，包括震级最大的1679年三河-平谷8级地震和震惊世界的1976年7.8级唐山大地震，该区一直是研究大陆构造形变和动力学的天然实验室。首都圈CSELF电磁台网的大部分台站都建立在华北地块的构造边界和重要构造带上，对于提升整个华北区域强震的监测预报研究水平具有重要作用。

1.2 野外测量与数据

CSELF台网的观测装置与大地电磁法(Magnetotelluric sounding，简称MT)(陈乐寿等，1990)基本相同，观测的物理量也与MT相同，不同的是其既可观测天然源电磁信号，也可观测人工源电磁信号。为获得首都圈CSELF台网区15个台站及其附近的深部电性结构特征，我们以每个台站为中心，根据台站的区域构造特征，在其两侧共布设4个宽频带MT测点，形成1条与区域构造基本垂直的短MT剖面，以获取穿过台站的电性结构背景。图1c 以大同台为例给出了台站结构探测中台站与宽频带MT测点的空间位置分布关系。按照上述方案，在首都圈的15个台站共布设60个宽频带MT测点，野外数据采集工作于2016年完成。事实上，由于首都圈经济与社会发展太快，台站附近的电磁干扰较之前也有所增强，为获得高质量的宽频带MT数据，对许多测点进行了多次重测，因此工作量远大于60个测点。

大部分台站及宽频带MT测点的数据质量较好，有效频点的周期范围为0.001～2000s，但个别台站离城市中心太近，受电磁干扰影响较严重。作为例子，图2 给出了大同台、夏县台、大连台和丰宁台和附近MT测点的xy与yx向视电阻率与阻抗相位曲线，测点编号均按照从北至南(或从西至东)由小变大。分析曲线发现，大部分台站两侧的曲线形态发生了较大变化，可能代表着台站位于电性边界带上，这与前文所述的很多台站处于区域构造边界带上的事实相符。为了更直观地了解台网区不同频率的视电阻率分布特征，我们对台网各台站和所有宽频带MT测点数据进行插值，给出了首都圈CSELF台网区不同频率的xy、yx向视电阻率及两者几何平均的视电阻率平面分布云图，图3 给出了频率为100Hz、0.01Hz的xy、yx向以及两者几何平均的视电阻率平面图，按照趋肤深度的定义，这2个频率分别代表了浅部与深部的视电阻率分布情况。可以发现，在郯庐断裂带、太行山重力梯度带两侧，无论高频还是低频的视电阻率变化均较为明显，暗示其不仅是构造边界，同时也是电性边界；华北北部的阴山-燕山造山带、东部的华北平原的视电阻率从高频到低频均呈现出较好的继承性。此外，通过长期动态监测台网区的视电阻率分布云图也可观测台网区视电阻率是否发生异常变化，这可为区域的地震监测预报提供重要参考依据。

2 区域构造的维性分析

图3 首都圈CSELF电磁台网区不同频率的视电阻率分布图Fig.3 Apparent resistivity maps at 100Hz and 0.01Hz of the CSELF network in the capital circle region.ρxy xy向视电阻率；ρyx yx向视电阻率；ρave xy与yx向视电阻率的几何平均

图4 不同台站每个测点的二维偏离度Fig.4 The 2D skewness at each site for Datong，Xiaxian，Dalian and Fengning staions.

图5 首都圈CSELF台网各台站的电性主轴方位玫瑰图Fig.5 Rose diagrams of geoeletrical strike at each station of the CSELF network in capital circle region.

此外，根据阻抗张量分解结果，我们也给出了首都圈CSELF台站的电性主轴方位统计玫瑰图，如图5 所示，该图较直观地反映了各台站地下可能的电性主轴方位。由于玫瑰图的电性主轴方位角存在90°的模糊性，还需结合台站区域的地质构造走向确定具体的电性结构走向。总体上，各台站的电性主轴方位与所处地质构造单元构造的走向较一致，且同一构造单元的台站电性主轴方位具有一定的相似性。例如，同处在山西断陷带的大同台、代县台和夏县台，电性主轴方位基本一致。同样地，郯庐断裂带附近的安丘台与莒县台的电性主轴方位也基本一致。

磁感应矢量是MT数据中的一个重要参数，其定义为垂直磁场与水平磁场的复数比值，包括实感应矢量与虚感应矢量。实感应矢量主要反映地下介质导电性分布的横向不均匀性，矢量大小反映横向导电性差异的大小，在Parkinson规范下(Parkinson，1962)，矢量方向指向电流会聚的区域，即高导体所在的位置。此外，在理想的二维条件下，实感应矢量所指方向垂直于构造走向，可用于判断区域电性主轴方位。我们也分析了各台站MT剖面数据的实感应矢量，定性地了解了各台站的地下电性结构特征。图6 以大同台和丰宁台为例给出了沿剖面MT测点的实感应矢量随频率的分布。结果显示，如果用趋肤深度公式进行简单估算，2个台站在中上地壳表现出较强的横向不均匀性，而中下地壳的介质横向相对均匀，且2个台站的南段在中上地壳呈现相对高导，这与后文反演获得的电性结构揭示的特征具有一致性。

图6 大同台和丰宁台沿剖面的实感应矢量拟断面图Fig.6 Pseudo-sections of the real induced vector along the MT profile at Datong and Fengning stations.

3 台站及其MT数据反演与电性结构

3.1 一维与二维反演

为获得每个台站及其附近的深部结构，我们既对台站的每个测点进行了一维反演，也对台站及其MT测点组成的短剖面数据进行了二维反演。

一维反演研究每个台站正下方介质的电性结构随深度的变化情况。本文采用自适应正则化反演方法(陈小斌等，2005)对各台站观测的xy和yx向视电阻率与相位分别进行了反演，并对2个方向数据的几何平均也进行了反演，最后在每个台站获取3条电阻率随深度的变化曲线。图7 以大同台为例，给出了大同台及其MT测点的数据一维反演结果。我们发现，在1km以浅深度，3条反演曲线基本重合，表明该深度范围内基本以一维构造为主，可能反映大同盆地的沉积层厚度；大同台与其北侧2个MT测点10km深度以浅的一维结构比较一致，而与其南侧的MT测点的反演结构有明显差异，可能反映大同台两侧的深部电性结构存在较为明显的横向不均匀性，这一点在后面的二维反演结果中得到了证实。

图7 大同台及其MT测点的一维反演结果Fig.7 The 1-D resistiviy structures for Datong station and MT sites.

图8 大同台MT剖面的二维反演Fig.8 Two-diminsional inversion of MT profile at Datong station.a 不同光滑因子的L曲线；b 大同台沿MT剖面的地壳电性结构；c 二维反演数据与模型响应拟断面图。ρObs观测视电阻率；ρCalc 模型响应电阻率；φObs 观测相位；φCalc 模型响应相位

采用目前广泛使用的非线性共轭梯度(NLCG)法(Rodietal.，2001)进行二维反演，整个反演过程均在MT-Pionner大地电磁可视化集成软件系统下完成(陈小斌等，2004)。进行二维反演前，根据前文阻抗张量的分析结果，将每个台站的剖面MT数据全部旋转至各自的电性主轴方位，旋转后平行于电性主轴方向的数据为TE模式，与之垂直方向上的数据为TM模式。在对实测数据进行二维反演时，不同极化模式数据的反演结果差别较大。针对如何选择反演数据的问题，国内外学者已进行过大量研究，一般认为，TE极化模式对深部结构变化更灵敏，TM极化模式对表层结构变化更敏感，采用TE+TM联合反演能利用各自的优势，增加数据约束量，提高模型的可靠度(Berdichevskyetal.，1998)。但是，当数据受三维结构影响明显时，单独进行TM模式反演比利用TE模式或TE+TM模式联合反演更合理，反演结果中的虚假构造会明显减少，所得结果能够比较准确地反映地下电性结构特征。对于实测数据的二维反演，应优先考虑采用TM模式数据进行二维反演，其次为TM+TE模式，一般不宜采用TE模式(Ledo，2006；蔡军涛等，2010)。根据数据的定性分析，由于部分台站数据的低频段受三维效应的影响，我们仅利用每个台站剖面数据的TM单极化模式数据进行二维反演，以获得沿剖面的地下电性结构。

为获得可靠的二维反演模型，我们通过改变不同的反演参数进行了大量反演计算，如改变正则化因子，在0.1～1000间取十几个值进行多次反演，通过L曲线分析确定最优的正则化因子(图8a)；调整视电阻率与相位的误差门槛值；将数据旋转不同的角度等。所有反演采用的初始模型为100Ω·m的均匀半空间，利用测点中心网格的自动生成技术构建带地形的初始模型网格(陈小斌等，2009)。按照上述反演方案，对首都圈15个台站的剖面数据进行了大量反演计算，反演工作量较大。最终获得了每个台站及其附近地下稳定可靠的二维电性结构模型。

3.2 台站地下二维电性结构

反演得到的15个台站的地下电性结构模型较稳健，能较好地描述各台站附近地壳的结构特征，并与区域地质构造具有一定的对应关系。张继红等(2019)对山东省无棣台、安丘台和莒县台附近区域的地下二维电性结构模型进行了解译。无棣台的区域地壳电性结构相对简单，具有成层性，符合稳定地块区的电性结构特征(赵国泽等，2010b)；而安丘台和莒县台区域的地壳电性结构相对复杂，呈现出活动地块边界带的结构特点，也反映了郯庐断裂中段的活动性和电性特点。本文将以大同台为例，简要分析解译该台站附近的地下二维电阻率模型。

图9 首都圈CSELF电磁台网区二维反演模型不同深度的电阻率切片Fig.9 Horizontal slices of 2-D electrical resistivity structure of the CSELF network in capital circle region at different depths.

为了更好地了解台网区电阻率的空间分布特征，根据各台站剖面数据的二维反演结果，通过插值获得了首都圈CSELF台网区不同深部的电阻率分布特征。图9 分别给出了深度为1km、5km、10km、15km、20km和30km的电性结构平面图。将其与台网区视电阻率频率响应结构(图3)进行对比发现，反演获得的电性结构与视电阻率响应结构的基本特征大致相符。华北平原从浅至深，地壳均表现为低电阻率特征，而北部的阴山-燕山造山带、东部的辽东-胶东地块的地壳表现为相对高阻。值得注意的是，西部山西断陷带的地壳从北部的大同盆地到南部的夏县盆地均表现为相对低阻的特征。电性特征的差异与不同构造区的地形、地貌特征及表层出露岩石的电性具有较好的对应关系。华北平原地势较平坦，沉积层比较发育，而阴山-燕山、辽东-胶东地区多为山区或丘陵地貌，主要出露变质结晶岩和火成岩，按照岩石导电性的差异，沉积层的电阻率往往比变质岩和火成岩低得多。横向上，1～15km深度的电阻率横向差异较大；15km深度以深处，电阻率的横向差异小于浅部。太行山重力梯度带作为华北东部地块与西部地块的分隔带，在电性结构上表现为电性边界带，两侧的电阻率结构差异明显，西侧的电阻率小于其东侧的电阻率。同样地，郯庐断裂带作为中国东部的一条巨型走滑断裂，其中段穿过了首都圈CSELF台网区，两侧的地壳电性结构也差异明显，西侧低阻、东侧高阻，且切割深度较深，表现为超壳大断裂。前人在首都圈CSELF台网区也开展了不少大地电磁研究，如魏文博等(2006，2008)完成了应县—商河(HB-MT01)大地电磁探测；詹艳等(2011)在河北石家庄地区开展了大地电磁探测。这些工作揭示太行山前断裂表现为较大的电性边界带，西侧太行山隆起地壳结构为高阻区，而东侧的华北平原为低阻区。穿过郯庐断裂中段的大地电磁剖面(叶高峰等，2009；张继红等，2010)也反映郯庐断裂带东、西两侧地壳的电性结构特征差异明显，东低西高。以往的认识与本研究的结果基本一致。此外，华北大地电磁测深阵列观测实验研究所揭示的华北地区的地壳电性结构特征与本研究的结果也有较好的一致性(Ouro-Djobo等，2018)。由于台站比较稀疏，台站之间的电性结构主要通过插值获得，这里不做过多讨论。

4 台网区的三维电性结构特征

近年来，随着计算机硬件和MT理论的不断发展，大地电磁三维正、反演技术已逐渐走向成熟，基于数据空间OCCAM(Siripunvarapornetal.，2005)、非线性共轭梯度(NLCG)(Newmanetal.，2000)、自适应正则因子的拟牛顿算法(AR-QN)(阮帅，2015)等三维MT正、反演程序包均已实现，并被广泛应用于实测资料的反演解释中(Xiaoetal.，2015，2017；Caietal.，2017；阮帅等，2020)。

我们利用极低频地震预测分系统项目研发的“大尺度三维电性结构快速反演软件”(toPeak)，采用基于NLCG算法的ModEM并行程序包(Kelbertetal.，2014)对首都圈CSELF台网区的结构探测数据进行了三维反演。由于台站间距较大，考虑到网格剖分的问题，对每个台站仅选择一个测点参与反演。反演输入数据为15个测点在1000～0.0001Hz频段范围内所有49个频点的xy和yx向视电阻率和相位，视电阻率和相位数据的门限误差分别为5%与1.43°，自动构建的反演网格有55(NS向)×70(EW向)×85(垂直方向，包括7个空气层)个，反演初始模型采用100Ω·m 的均匀半空间。经过100次迭代，最终模型的数据拟合均方差RMS为2.36。图10 给出了首都圈CSELF台网区三维电性结构不同深度的平面图。与台网区二维反演电阻率平面分布相比(图9)，三维反演呈现的台网区不同深度的电阻率空间分布特征与之较为相似，差异主要存在于台站间的结构特征，这主要是由于台站间距较大，台站间没有数据约束所致，因此三维反演结果中台站间的电阻率都趋于初始模型的值。而二维反演结果中台站间的电阻率通过插值运算获得，虽然信息更为丰富，但可信度不高。

由于参与三维反演的数据有限，本文给出的首都圈CSELF电磁台网区三维电性结构背景是比较粗略的，对其结构特征不做细致讨论，但该模型的获取为首都圈CSELF电磁台网在地震监测预报中发挥作用提供了重要的基础资料。总而言之，通过对首都圈台网区观测数据的一维、二维和三维反演研究，本文揭示了各台站及其附近的地壳电性结构和整个台网区的电阻率空间分布特征，通过与前人的研究结果进行对比与分析，本研究的结果具有较高的可靠度。

图10 首都圈CSELF电磁台网三维反演模型不同深度的切片Fig.10 Horizontal slices of 3-D electrical resistivity structure of the CSELF network in capital circle region at different depths.

5 结论

本研究通过野外数据采集、处理与分析，一维、二维和三维反演，获得了可靠的首都圈CSELF电磁台网区各台站及其附近的地壳电性结构以及整个台网区不同深度的电性结构。结果表明，电性结构与台网区的构造具有很好的对应性，如北部阴山-燕山造山带、东部胶-辽地块以及中部太行山地区的台站地壳结构主要表现为高阻特征，而华北平原、山西断陷盆地区的台站地壳电性结构以高导特征为主。显然，各台站的电性结构特征与所处地区的地形、地貌以及区域地质构造特征具有较好的关联性。整个首都圈台网区的电阻率空间分布特征展示了台网区的中上地壳具有较强的横向不均匀性，特别是在太行山重力梯度带和郯庐断裂带两侧，这意味着它们不仅是重要的构造边界带，也是电性边界带。

首都圈CSELF电磁台网的大部分台站不仅位于重要的地质构造带上，同时也位于地震活动带上，如郯庐断裂带、张家口—渤海地震带以及山西裂谷地震带等。这些地区的台站区域地壳电性结构相对复杂，但基本具有上地壳多表现为高阻而中下地壳发育高导层的特征，并且区域的历史地震震中多位于中上地壳低阻与高阻交界处的高阻一侧，地震的空间分布与地下结构具有较大关系。由于文章篇幅有限，本文仅对大同台站及附近的地壳电性结构进行了解译，并讨论了大同-阳原地震的深部孕震结构。

首都圈CSELF电磁台网区的地下电性结构探测提供了整个区域的电性结构背景，对地震活动区的孕震环境、地震电磁异常现象产生的机理以及地震预测预报研究具有重要作用。未来，我们应长期不断地监测台网区的地下电性结构变化特征，真正地让它在防震减灾工作中发挥作用。

致谢野外工作人员在数据采集过程中的辛勤付出为本研究提供了基础数据；中国地震局地质研究所信息中心为大地电磁三维反演工作提供了超算平台。在此一并表示感谢！