引言

震前电磁异常现象已经被大量的破坏性地震震例所证实，且已成为地震日常监测预报特别是短临预报最有效的手段之一（丁鉴海等，2006）。依据岩石破裂实验可知，孕震区介质的电导率与应力应变关系密切（郝锦绮等，1989）。当有流体存在时，孕震区介质在应力作用下将经历微破裂阶段，介质中产生很多微破裂，其体积发生膨胀（扩容），此时将有大量流体流入介质中，从而改变介质的电性结构（Nur，1972；Scholz et al，1973；车用太，鱼金子，2006）。因此，基于地磁台站观测的时变记录，根据电磁感应理论，利用地磁测深法可以研究孕震过程中地下电性结构的变化，亦可为地震预测提供判据。

Parkinson （1959）于1959年提出了帕金森矢量的概念，指出其大小和方向可以有效地反映地下电性结构的特征。Schmucker （1970）引入了转换函数的概念，将频率域地磁短周期变化磁场分为正常场部分和异常场部分，认为在外源场均匀的情况下，正常场部分由外源场和水平横向均匀介质的感应场组成，异常场部分仅与介质的横向不均匀性相关。转换函数以及由转换函数组成的感应矢量大小反映了地下电性结构的横向不均匀程度，感应矢量方向指向高导体方向，因此可以通过研究转换函数的变化研究地下电性结构（陈伯舫，1974；侯作中，史铁生，1984）。自20世纪70年代以来，国内外研究人员通过转换函数来研究地震伴随的地下电性结构变化，已取得较多成果（Yanagihara，1972；徐文耀等，1978；Rikitake，1979；Beamish，1982；龚绍京，吴占峰，1986；曾小苹，林云芳，1995；邢西淳等，2008；龚绍京等，2015）。总结前人对震前转换函数随时间变化的研究认为，震前转换函数异常周期主要集中在10 —48 min周期，对于周期超过90 min的转换函数基本无明显异常变化。转换函数的频率特性曲线在地震发生前或临震前后出现上升或大的起伏而后有一定程度的恢复。由转换函数获得的感应矢量方位角存在震前逐渐指向震中，震后回到正常方位或继续变化至新的方位的现象。

前人对震前转换函数变化的研究主要基于多个固定地磁台站资料，台站间距离和震中距均较大，无法描绘转换函数异常变化的范围信息。本文拟利用滇西北地磁台阵8个磁通门磁力仪测点记录的短周期地磁信号资料，使用Chave和Thomson （2004）提出的有界影响估计方法估计该台阵区域内2013年洱源M_S5.5地震前后的地磁转换函数，探讨地震前后不同震中距下台站转换函数的变化。

1.   观测资料

滇西北地区位于川滇地块边缘的一个特殊位置，构造活动强烈。不同走向、规模和活动强度的断裂纵横交错，地震活动比较频繁。为了开展震磁关系研究，滇西北地区呈“十”字型布设了由8个GM4磁通门磁力仪测点组成的小台阵，地磁观测时间为2012年3月至2013年12月。台阵仪器采样间隔为1秒，峰-峰值噪声小于0.1 nT。台阵记录的原始数据为秒采样的地磁垂直分量（Z）、水平分量（H）和磁偏角（D）的变化值。将原始数据中明显的非天然磁场的尖峰和台阶剔除，生成预处理秒采样数据，然后对其进行高斯滤波处理，转换成为预处理分数据。本文的研究工作便基于预处理分数据开展。

在台阵观测期间，台阵区内发生过两次M_S≥5.0地震：2013年3月3日云南洱源M_S5.5地震，震中（25.9°N，99.7°E）位于维西—乔后断裂带附近，震源深度约为9 km；2013年4月17日云南洱源M_S5.0地震，震中（25.9°N，99.8°E），震源深度约为5 km。这两次地震为主-余型地震序列，震源机制解显示该震源断层为正断兼右旋走滑型破裂（常祖峰等，2014；黄小龙等，2015）。图1给出了滇西北地磁台阵测点分布以及两次地震的震中和余震的分布，表1列出了台阵各测点与洱源两次M_S≥5.0地震震中的距离。

图  1  滇西北地磁台阵测点分布及洱源两次M_S≥5.0地震的震中和余震分布图

F₁：怒江断裂；F₂：澜沧江断裂；F₃：维西—乔后断裂；F₄：龙蟠—乔后断裂；F₅：鹤庆—洱源断裂； F₆：小金河—丽江断裂；F₇：红河断裂；F₈：程海—宾川断裂

Figure  1.  Distribution of measuring points of geomagnetic array in northwest Yunnan Province and epicenters and aftershocks of two M_S≥5.0 earthquakes in Eryuan

F₁：Nujiang fault；F₂：Lancangjiang fault；F₃：Weixi-Qiaohou fault；F₄：Longpan-Qiaohou fault；F₅：Heqing-Eryuan fault； F₆：Xiaojinhe-Lijiang fault；F₇：Honghe fault；F₈：Chenghai-Binchuan fault

下载: 全尺寸图片 幻灯片

表  1  滇西北地磁台阵各测点与洱源两次M_S≥5.0地震震中的距离

Table  1.  The distances between the epicenters of two M_S≥5.0 earthquakes in Eryuan and the measuring points of the geomagnetic array in northwest Yunnan Province

测点名称	与洱源MS≥5.0地震震中的距离/km			测点名称	与洱源MS≥5.0地震震中的距离/km
	MS5.5	MS5.0			MS5.5	MS5.0
云龙	37.15	42.85		牛街	47.12	50.29
剑川	52.73	52.55		鹤庆	76.53	77.34
炼铁	8.34	11.52		下关	64.50	55.88
洱源	22.94	27.38		丽江	116.05	118.71

下载: 导出CSV 

| 显示表格

2.   计算方法及资料处理

2.1   转换函数

在中低纬度地区，单台垂直场转换函数可以近似表示为（Schmucker，1970）

式中，$\omega $为周期，Z（ω），H_x（ω）和H_y（ω）分别为地磁场垂直分量、南北向水平分量和东西向水平分量的频谱值，A和B为转换函数，ε为残差，式中各元素均为与周期有关的复数。当某一周期ω内有多个地磁短周期事件时，可将上式改写为矩阵形式，即

式中，Z为垂直分量矩阵，H为水平分量矩阵，T＝[A，B]^T为转换函数矩阵，r为残差。利用最小二乘估计（Sims et al，1971），可以得到转换函数的估计值为

式中，“*”表示共轭转置矩阵。

利用转换函数可以得到感应矢量，感应矢量用于描述横向地下电性结构，其中实感应矢量指向地下高导体方向，其大小可描述地下电性结构的横向不均匀性，计算公式为

式中：Ar，Br分别为转换函数A，B的实部；α为方位角，其反方向为实感应矢量方向；I为地磁变化优势面的倾角；L为实感应矢量的大小。

根据高斯-马尔可夫理论，当残差互不相关且方差相同时，最小二乘估计是最优无偏估计。但是实际的地磁数据中噪声并不满足高斯分布假设，同时最小二乘估计对少量的异常数据、磁场变量中不相关噪声以及模型本身的不足非常敏感，这会导致转换函数估计结果产生严重的偏移。为了消除不相关噪声对转换函数估计的影响，Gamble等（1979）提出了远参考大地电磁数据处理方法，该方法基于远参考道信号与测点信号相关而噪声不相关的原理，利用远参考道数据参与转换函数估计，其转换函数估计值为

式中，H_r代表远参考道磁场。

为了提高转换函数估计的稳定性，许多科研人员提出多种稳健估计方法（Robust方法）来改善转换函数估计的稳定性。其中，M估计（Egbert，Booker，1986）是一种加权最小二乘估计，其根据观测误差和剩余功率谱的大小对数据进行加权，降低离群点对转换函数估计的影响，其迭代形式为

式中，k为迭代次数，V为权重矩阵，一般选择胡贝尔（Huber）权重函数进行权重矩阵计算。

Chave和Thomson （2003，2004）将M估计与帽子矩阵对角线统计分析相结合，发展出有界影响估计法。该方法可以同时消除输出道和输入道中的干扰及二者中存在的相关噪声，并且能够最大程度地减小离群点和杠杆点（leverage）的影响，从而在强噪声背景下增强结果的可靠性和稳定性。其迭代形式为

式中，V为残差权重矩阵，W为杠杆权重矩阵。杠杆权重矩阵W的计算方法参考Chave和Thomson （2003）的文章。

有界影响估计法在处理大地电磁数据时效果较好，王桥和黄清华（2016）在处理华北地磁数据时验证了有界影响估计法同样适用于处理地磁测深数据。本文参考其验证思路设计了有界影响估计稳定性的验证方法。首先选择台阵内数据质量较高的云龙测点连续60天的地磁变化三分量数据，求取6天长度的平均时间序列，得到测试序列${\bar h_{{x}}}\!\!\!\!{\text{（}}\!{{t}} \!{\text{）}}\!\!\!\!$，${\bar h_y}\!\!\!\!{\text{（}}\!t \!{\text{）}}\!\!\!\!$ 和 $\bar z\!\!\!\!{\text{（}}\!t \!{\text{）}}\!\!\!\!$；然后分别在水平变化分量${\bar h_{{x}}}\!\!\!\!{\text{（}}\!{{t}} \!{\text{）}}\!\!\!\!$，${\bar h_y}\!\!\!\!{\text{（}}\!t \!{\text{）}}\!\!\!\!$ 和垂直变化分量 $\bar z\!\!\!\!{\text{（}}\!t \!{\text{）}}\!\!\!\!$ 中添加5%的高斯白噪声；再利用有界影响估计、最小二乘估计和加权最小二乘估计三种方法估计转换函数A和B。图2给出了所得转换函数A的实部Ar的对比结果。由图中可以看出：当输入道中存在噪声时，三种方法的结果较为一致且稳定性较好；当输出道中存在噪声时，三种方法所得转换函数估计的稳定性均有所下降，但有界影响估计的稳定性明显高于其它两种估计方法。

图  2  不同求取方法所得地磁转换函数的稳定性对比

（a） 在水平分量中加入5%高斯白噪声的转换函数实部Ar；（b） 在垂直分量中加入5%高斯白噪声的转换函数实部Ar

Figure  2.  Comparison of geomagnetic transfer function stability estimated by three different methods of BI estimator，least squares and weighted least squares

（a） Real component Ar of geomagnetic transfer function result with 5% Gaussian white noise in horizontal component；（b） Real component Ar of geomagnetic transfer function result with 5% Gaussian white noise in vertical component

下载: 全尺寸图片 幻灯片

2.2   资料处理

根据王辉等（2018）对长周期大地电磁远参考道处理的研究，在中低纬度利用相距较远的固定地磁台站数据，对地磁台阵测点进行远参考道处理是可行的，且选取的远参考点与测点相距数百千米。因此本文在转换函数估计过程中，选取云南省固定地磁台站通海台（24.1°N，102.75°E）作为远参考台，利用远参考道原理有效地消除各测点水平分量信号中不相关噪声的干扰。采用Chave和Thomson （2004）提出的有界影响估计，并使用其编写的BIRRP程序对各测点进行转换函数估计。该程序首先使用一个短的自回归滤波器进行滤波，并根据感兴趣的周期长度对观测数据进行分段，然后用斯莱伯恩（Slepian）窗（Thomson，1982）对观测数据加窗进行傅里叶变换。选择研究的频点，通过有界影响方法估计，得到不同周期的转换函数，最后采取刀切法（Jackknife）计算转换函数的标准误差（Eisel，Egbert，2001）。

由于转换函数在一个较短时间段内可认为是相对稳定的，因此本文选择7天为一个时间窗进行一组转换函数估计，将时间窗向后滑动四天，进行下一组转换函数估计。在时间窗内人工选取干扰小且尽可能连续的观测数据进行处理，对于选取时段内因地电观测干扰或其它原因造成的短暂缺数，本文采用线性插值的方法进行插值以保证观测记录的连续性。若连续缺数时间较长则将时间序列分段进行转换函数估计以保证转换函数估计的可靠，并选取每个时间窗中间的时刻作为转换函数记录时间点。

3.   震例分析

本文计算了滇西北地磁台阵8个测点自2012年6月1日至2013年6月30日的转换函数A和B，计算周期范围为14.2—51.2 min。根据各测点的转换函数估计结果，牛街、下关测点的转换函数实部Ar在0附近变化，其它测点的转换函数实部Ar的绝对值均较小，约为0.1。同时在地磁中纬度地区，转换函数A受季节变化影响比转换函数B更明显（Vargas，Ritter，2016），同时也会对地下电性结构的判断产生干扰，因此本文主要研究各测点转换函数实部Br的变化。

图3和图4分别为云龙测点转换函数A，B的实部Ar和Br随时间的变化曲线。根据异常成片原则（龚绍京等，2015）：① 要有多个周期同时出现异常；② 要有连续多个点同时出现异常。云龙测点在2012年6月至9月转换函数较为稳定，没有明显异常变化。2012年10月至2013年1月，由于仪器问题，该测点大量缺数，无法进行转换函数估计。在M_S5.5地震震前1—2个月，该测点多个周期转换函数实部Br的绝对值逐渐变大，在地震发生前转换函数实部Br发生转折，其绝对值逐渐减小，异常幅度超过1倍标准差，而转换函数实部Ar未见成片异常出现。云龙测点在研究周期范围内，随着转换函数周期的增大，其估计值的误差也整体随之变大，其它测点也有相同的误差变化趋势，表明该区长周期信号信噪比较低。

图  3  云龙台不同周期T转换函数实部Ar的时间变化曲线

Figure  3.  The time curves of the real component Ar of transfer function for different periods T at the Yunlong station

下载: 全尺寸图片 幻灯片

图  4  云龙台不同周期T转换函数实部Br的时间变化曲线

Figure  4.  The time curves of the real component Br of transfer function at the Yunlong station for different periods T

下载: 全尺寸图片 幻灯片

图5以T＝18.3 min为例给出5个测点转换函数实部Br随时间的变化曲线（剑川、炼铁、鹤庆测点缺数较多，故未给出）。可以看到云龙测点转换函数实部Br在M_S5.5地震前存在明显的变化过程，洱源测点在M_S5.5地震震前1个月左右，其转换函数实部Br的绝对值也存在逐渐变大的异常，异常幅度超过1倍标准差，至临震前转换函数实部Br的绝对值则减小；而牛街、下关、丽江测点的转换函数实部Br却并无明显的异常变化。

图  5  周期T为18.3 min时各测点转换函数实部Br的时间曲线

Figure  5.  The time curves of Br at the period of 18.3 minutes for five geomagnetic stations

下载: 全尺寸图片 幻灯片

研究区域在2013年3月3日洱源M_S5.5地震前8个月未发生M_S＞3.0地震，认为研究区域处于地震平静期（高琼等，2014）。并且根据各测点数据的实际记录情况，选择M_S5.5地震前六个月（2012年8月18日）和震前半个月（2013年2月16日）两个时间节点绘制各测点不同周期的实感应矢量并进行对比（图6），由于测点数据因仪器故障导致数据记录缺失，个别测点选用相邻时间节点的实感应矢量。结果显示，各测点不同周期实感应矢量无明显方向变化。云龙、剑川测点不同周期的实感应矢量大致指向东向，与其余测点的实感应矢量方向大致相反。

图  6  不同周期实感应矢量随时间的变化

Figure  6.  The real induction vector varies with time at different periods

（a） T＝18.3 min；（b） T＝25.6 min；（c） T＝36.6 min；（d） T ＝51.2 min

下载: 全尺寸图片 幻灯片

图7给出了不同测点在T＝18.3 min时的实感应矢量长度随时间的变化曲线，其它周期类似，故未给出。可以看到云龙测点在M_S5.5地震前1—2个月实感应矢量长度存在明显的上升-下降变化趋势，且在M_S5.0地震前也存在这一变化趋势。而其它测点在M_S5.5地震前实感应矢量长度在均值附近浮动，并未出现明显的异常变化。

图  7  测点周期为18.3 min时实感应矢量长度L的时间曲线

Figure  7.  The time curves of L at the period of 18.3 minutes for five geomagnetic stations

下载: 全尺寸图片 幻灯片

4.   讨论与结论

本文利用转换函数方法分析了滇西北地磁台阵区内发生的两次M_S≥5.0地震，在所求取的周期范围内，云龙测点和洱源测点在M_S5.5地震前转换函数实部Br有明显的异常变化，震前转换函数Br的绝对值存在逐渐增大—转折—恢复的过程，异常幅度超过1倍标准差。同一测点不同周期的转换函数实部Br变化幅度、形态与异常发生时间并不完全相同，表明孕震过程中不同深度处的电性结构变化十分复杂。两测点转换函数实部Br符号相反，表明两测点间存在电导率异常。根据钱复业等（1982）对震磁异常范围的统计研究，对于M_S5.5地震，震磁异常范围大约为50—60 km。同时龚绍京和吴占峰（1986）通过研究唐山地震周边台站的转换函数认为，转换函数异常发生在余震区附近。除了距离震中最近的炼铁测点因缺数过多无法给出计算结果外，云龙、洱源以外的其它测点的震中距均接近或大于50 km。由洱源地震的余震分布可以看出（图1），余震区呈北东向的条带状，在余震区附近的云龙测点和洱源测点，其转换函数实部Br的绝对值出现异常变化，其它距离余震区较远的测点的转换函数实部Br无明显异常变化，表明转换函数异常区域可能在余震区周围，与龚绍京和吴占峰（1986）的结果相符合。

根据实感应矢量的定义，实感应矢量方向指向高导体，大小（即矢量长度）表示地下电性结构的横向差异程度。云龙测点和洱源测点在震前或临震时实感应矢量变大，在震后逐渐恢复，表明孕震过程中震源区地下电性横向结构差异增大，震后地下电性结构横向差异逐步恢复。对比平静期与震前的实感应矢量方向，并未发现实感应矢量在震前向震中偏转的变化。在计算周期范围内，云龙、剑川、炼铁、洱源、牛街、鹤庆、丽江测点的实感应矢量均向剑川、鹤庆方向汇聚；下关测点实感应矢量指向西，且感应矢量长度较大，由此推断在研究区域内沿乔后—维西断裂存在明显的高导带或高导体异常，其走向大致为NNW向。根据各测点不同周期实感应矢量的变化，高导体在不同深度的延伸是复杂的。该异常体的空间展布与大地电磁测深研究和短周期地磁变化（孙洁等，1989；袁伊人等，2015）结果较为吻合，认为该异常体以剑川鹤庆为中心，沿滇西北褶皱带呈NNW向展布，可能是由于上地幔物质隆起所致。洱源M_S5.5地震实感应矢量方向在震前没有明显变化，可能与研究区域的地质背景有关。

一般认为地磁短周期变化转换函数异常是由于膨胀区里的水侵入新产生的裂隙引起震源区及周围电导率发生变化从而导致转换函数的变化（Rikitake，1979；龚绍京，吴占峰，1986）。洱源M_S5.5地震震源区为富含热泉的高热流区域（吴乾蕃等，1988），在应力作用下，流体侵入微破裂中，可能导致震源区地下电性结构发生变化。虽然该震例转换函数异常的范围较小，但通过转换函数的变化特征，可以观察到震源区在孕震过程中地下电阻率先下降，临震或震后逐减增大的过程，与孕震理论研究的扩容模式（Nur，1972；Scholz et al，1973）给出的电阻率变化类似。因此，地磁转换函数是了解孕震区地下电性结构变化的有效途径之一。

本文地磁数据来自中国地震局地球物理研究所国家地磁台网中心，江苏省地震局张秀霞高级工程师、陈健工程师、郭宇鑫助理工程师、李玉助理工程师在地磁数据处理方面给予了热情帮助，作者在此一并表示诚挚的感谢。