地震的孕育及发生过程中,地下能量的积累与释放过程必然会引起岩石圈介质性质,如温度、应力等状态的改变。这一性质的改变将通过各种震磁效应引起岩石磁化率的改变(黄平章等,1990; 郝锦绮等,1992,1993; 侯登录等,2003; 李海燕,张世红,2005),从而在地面产生能观测到的岩石圈磁场异常变化。因此通过监测地面的岩石圈磁场异常变化研究地震活动性是地磁学科主要研究方向之一。我国的地震学者通过监测地磁台站的连续观测数据,发现了某些可能的地震地磁异常现象,取得了如“低点位移”“转换函数”等理论(冯志生等,2000,2004; 刘哲,2002; 史勇军等,2002; 王亶文,2003; 陈伯舫,2003; 丁鉴海等,2004; 王世昌等,2004; 张淑兰,江在森,2004; 吴小平等,2004; 李树华,2004; 李琪等,2005; 张继红,乔惠珍,2005; 蒋延林等,2005; 杨涛等,2005),但由于地磁台站的空间间隔较大,这一类研究对空间的辨识率不高。流动地磁研究者也使用空间高分辨率的流动地磁数据进行了地震地磁相关性的分析(解用明等,2001,2002,2005; 詹景立,黄健,2004),但仅重点分析以震中为中心的某一段时间-空间区域的变化,缺乏全时间-空间尺度的辨析,且不同地区、不同时间的地震活动对应的地磁变化形态并不一致,这可归因于其中的某些变化是磁场的正常变化。自2010年初以来,中国地震局流动地磁技术团队在大华北、南北地震带和南北天山3个区域持续开展了年尺度的矢量地磁复测,发现了一些与地震相关地磁异常现象,如围绕地震震中出现高梯度分布或者磁偏角零变线等现象(顾春雷等,2010,2012; 倪喆等,2014a,b,c; 苏树鹏等,2016)。本文主要阐述了流动地磁技术团队采用的的数据处理过程,并分析数据处理过程中误差的可能数值范围。

按照来源的不同,地球磁场可分为内源磁场和外源磁场。内源磁场可细分为源于地球液体外核磁流体动力学过程的主磁场和源于地壳及上地幔物质磁性的岩石圈磁场。外源场可细分为源于固体地球之外的各种电流体系产生的磁场和固体地球内部产生的感应磁场。

由于岩石圈磁场包含了地壳和上地幔物质状态和结构的丰富信息,因此通过对岩石圈磁场的研究可以推进对地球演化、地质构造演变及成矿研究等相关领域的基础研究。从地震监测预报目的来看,发生在地表60 km范围内的浅源地震往往比同级别的深源地震会带来更大破坏(傅淑芳等,1980; 许忠淮,2014),所以地磁场的地震监测目标重点应为对应浅源地震深度的岩石圈磁场。岩石磁学的实验室研究表明:岩石的状态,如温度、应力状态、岩石颗粒度及铁磁性物质含量等发生改变时,其磁性会发生相应的改变,即产生如压磁效应、膨胀磁效应、感应磁效应、热磁效应等震磁效应。地震的孕育及发生过程中,地下能量的积累与释放过程必然会引起岩石圈介质性质状态的改变,并产生能观测到的磁场异常变化。因此通过监测岩石圈磁场的空间分布及时间变化的异常特征来研究地震活动性是我们的目标。

合理、准确地描述岩石圈磁场的空间分布及时间变化需要尽可能精准地从观测数据中消除各种外源磁场以及地球主磁场的影响。因此,中国地震局制定了观测技术规范 中国地震局监测预报司.2015.关于印发《流动地磁测量基本技术要求(试行)》的通知.以及日变通化、长期变化改正的数据处理流程,以便得到尽可能精确的岩石圈磁场数据。

日变通化的主要目的是消除观测数据中的规则日变场和其他外源场成分,得到观测地磁场数据中的内源场部分。其计算公式为:

F_{site_int_Ti}=F_{site_mea_Ti}-F_{site_ext_Ti}

=F_{site_mea_Ti}-F_{obs_ext_Ti}

=F_{site_mea_Ti}-(F_{obs_mea_Ti}-F_{obs_int_Ti})

=F_{site_mea_Ti}-(F_{obs_mea_Ti}-F_{obs_int_T0})

=F_{site_mea_Ti}-(F_{obs_mea_Ti}-F_{obs_mea_T0})(1)

式中:F_{site_int_Ti}为T_i时刻观测的某一地磁要素的日变通化结果; F_{site_mea_Ti}为该地磁要素T_i时刻的测量值; F_{site_ext_Ti}为该地磁要素T_i时刻测量值中的外源场部分; F_{obs_ext_Ti}为距测量位置最近地磁台站的该地磁要素T_i时刻测量值中的外源场部分; F_{obs_mea_Ti}为距测量位置最近地磁台站的该地磁要素T_i时刻测量值; F_{obs_int_Ti}为距测量位置最近地磁台站的该地磁要素T_i时刻观测值中的内源场部分; F_{obs_int_T0}为距测量位置最近地磁台站的该地磁要素日变通化日地方时午夜00:00—03:00的测量值中的内源场部分; F_{obs_mea_T0}为距测量位置最近地磁台站的该地磁要素日变通化日地方时午夜00:00—03:00的测量值。

式(1)给出了日变通化的数据推导过程,从中可看出在日变通化过程包含了3次近似处理,其理由为:

(1)用台站T_i时刻的外源磁场部分替代了同时刻测点的外源磁场部分。采用这种近似处理是由于包含规则日变场的外源磁场空间分布较为均匀,空间信号尺度较大(徐文耀,李卫东,1994)。

(2)用台站T₀时刻的内源场部分替代了台站T_i时刻的内源场部分,即认为在T_i至T₀的时间内,台站的内源场没有发生变化。这种近似要求T₀距离T_i的时间长度不宜较长,否则会引起较大的误差。

(3)用台站T₀时刻的测量值替代了台站T₀时刻的内源场值,这一近似要求对T₀的选择较为严苛,我们选择磁场日变形态极为平静日的当地午夜00:00—03:00作为T₀时刻,以求达到该时刻磁场测量值中包含的外源场部分最少。

对应以上3点近似处理的要求,日变通化过程中需要注意以下2点:

(1)当台站距离测点较远时,应考虑在测点附近架设磁场连续观测仪器,以减小日变场的空间不均匀引起的误差。

(2)日变通化日的选择应该非常仔细,需要选取距测量日期较近时间范围内的磁场平静日作为日变通化日。根据经验,通化日距测量日不宜超过30 d。

流动地磁的监测目标为岩石圈磁场的空间分布及时间变化,目前主要分析数据为岩石圈磁场的年度变化,即岩石圈磁场与去年的差异变化。其计算公式为:

ΔF_{Lith_T2-T1}=F_{Lith_T2}-F_{Lith_T1}

=(F_{Int_T2}-F_{Main_T2})-(F_{Int_T1}-F_{Main_T1})

=(F_{Int_T2}-F_{Int_T1})-(F_{Main_T2}-F_{Main_T1})(2)

式中:ΔF_{Lith_T2-T1}为某一地磁要素T₂时刻至T₁时刻的岩石圈磁场之差; F_{Lith_T2}、F_{Lith_T1}分别为该地磁要素在T₂和T₁时刻的岩石圈磁场; F_{Int_T2}、F_{Int_T1}分别为该地磁要素在T₂和T₁时刻的内源磁场; F_{Main_T2}、F_{Main_T1}分别为该地磁要素在T₂和T₁时刻的主磁场; T₂和T₁时刻的内源磁场F_{Int_T2}和F_{Int_T1}即为T₂和T₁时刻的日变通化结果; 而F_{Main_T2}-F_{Main_T1}为T₂和T₁时刻的主磁场之差,即为T₂和T₁时刻的地磁场长期规律变化。

流动地磁技术团队建立了1995年1月1日至今的中国大陆地磁场长期变化NOC模型,该模型基于1995年1月1日至今的中国地区36个台站的地磁场要素H、D、Z绝对连续观测时均值,并不断追加后续的台站观测数据。经过对台站原始数据的预处理,舍去了台站数据完整率与正确率较低的5个地磁台的记录数据,同时在计算区域的边界位置均匀地布设了15个边界控制点,用国际地磁参考场(IGRF)计算这15个边界控制点上的地磁场长期变化数据,最终将总计46个空间点位上的地磁数据作为NOC分析和建立长期变化模型的基础数据。在经过必要的处理后,将非磁扰日地方时 00:00—03:00每10 d的观测时均值的平均值作为具体研究对象。这批数据在空间上覆盖了中国大陆地区,在时间上跨越了接近2个太阳活动周期。

以各台站地磁要素X、Y、Z非磁扰日地方时00:00—03:00的每10 d的计算平均值分别构建地磁要素观测数据矩阵A^X、A^Y和A^Z,并进行了NOC计算:

A^X_46×758=U ^'X_46×46·Λ^'X_46×46·V ^'X_46×758^T

A^Y_46×758=U ^'Y_46×46·Λ^'Y_46×46·V ^'Y_46×758^T(3)

A^Z_46×758=U ^'Z_46×46·Λ^'Z_46×46·V ^'Z_46×758^T

式中:数据矩阵A^X、A^Y和A^Z的每一行为某一台站、相应地磁要素的时间变化; Λ^'X、Λ^'Y和Λ^'Z为观测数据矩阵A^X、A^Y和A^Z的本征值矩阵,是能量分配的描述; 由于观测数据矩阵A^X、A^Y和A^Z的构建所决定,U ^'X、U ^'Y和U ^'Z为矩阵A^X、A^Y和A^Z的空间域本征向量矩阵,是各地磁要素空间分布特征的描述; V ^'X、V ^'Y和V ^'Z为矩阵A^X、A^Y和A^Z的时间域本征向量矩阵,是各地磁要素时间变化特征的描述。

选用三次样条函数对各阶本征值能量的时间变化进行拟合插值,计算公式如下:

S_i(x)=(M_i-1)/(6h_i)(x_i-x)³+(M_i)/(6h_i)(x-x_i-1)³+((y_i-1)/(h_i)-(M_i-1)/6h_i)(x_i-x)+((y_i)/(h_i)-(M_i)/6h_i)(x-x_i-1)(4)

式中:M_i=S^″(x_i), i=0,1,2,…,n; h_i=x_i-x_i-1, i=1,2,…,n; x_i-1≤x≤x_i; M_i是待定参数,其中M₀=M_n=0,其余的M_i由下式得到:

{γ₁M₀+2M₁+α₁M₂=β₁

γ₂M₁+2M₂+α₂M₃=β₂

……

γ_n-1M_n-2+2M_n-1+α_n-1M_n=β_n-1(5)

式中: γ_i=(h_i)/(h_i+h_i+1), α_i=1-γ_i, i=1,2,…,n; β_i=6/(h_i+h_i+1)((y_i+1-y_i)/(h_i+1)-(y_i-y_i-1)/(h_i)), i=1,2,…,n。

选用曲面样条函数对各阶本征值能量的空间分布进行拟合插值,计算公式如下:

W(x,y)=a₀+a₁x+a₂y+∑^N_i=1F_ir²_iln(r²_i+ε)

∑^N_i=1F_i=∑^N_i=1x_iF_i=∑^N_i=1y_iF_i=0(6)

式中: r²_i=(x-x_i)²+(y-y_i)², ε=10^-7。

受限于原始数据的时间范围和空间范围,其有效时间范围为1995年1月1日至今,有效空间范围为(18°～54°N,73°～136°E)。将(2)式中的T₂和T₁代入(4)式,将测量数据的空间位置代入(6)式,并将结果带入(3)式,并选取适当的截断阶数,即可计算得到(2)式中F_{Main_T2}-F_{Main_T1}的数值,并完成长期变化改正,得到岩石圈磁场的差ΔF_{Lith_T2-T1}。

图1为各年度经过日变通化和长期变化改正后年度变化的中值(用黑色柱状图表示)和偏离中值的标准偏差(用黑色垂直线段表示),如2011年位置显示的为经过日变通化和长期变化改正后,2011年度数据减去2010年度数据的中值和标准偏差计算结果。按照传统的地磁学理论,除非发生地震、火山活动等剧烈的地壳物质运动,否则岩石圈磁场的基本不发生改变(徐文耀,2003)。因此,我们将岩石圈磁场的趋势性偏差,即中值作为判断数据处理方法流程的整体误差,将与中值的标准偏差作为数据偏离程度的衡量,即岩石圈年度改变量值的大致范围。从计算结果可得到如下几点:

(1)2011—2016各年度X、Y、Z要素数据处理误差均小于10 nT。

(2)除2012年度的X要素、2013年度的Y要素和2016年度的Z要素,其余大部分年度的各要素数据处理误差小于5 nT。

(3)各年度、各要素的数据处理标准偏差为10 nT左右。各年度磁场年度变化均正负值共存,在偏差较大的年份会出现正值较多(2016Z)或负值较多的情况(2012X、2013Y)。要素,其余大部分年度的各要素数据处理误差小于5 nT。

(3)各年度、各要素的数据处理标准偏差为10 nT左右。各年度磁场年度变化均正负值共存,在偏差较大的年份会出现正值较多

图1 年度误差柱状图及标准偏差
Fig.1 The annual error histogram and standard deviation

本文阐述了中国地震局流动地磁团队现有的数据处理流程,该流程依据地磁学的相关知识,利用台站连续观测数据,建立了独立的地磁学日通化和长期变化模型,形成了稳定可靠的数据处理流程,得到了高精度的岩石圈磁场年度变化数据。

选用测点最近地磁台日变通化和自然正交模型长期变化改正的数据处理方法过程中需要有诸多特别需要注意的地方。该处理流程带来的整体误差小于10 nT,大部分数据误差保持在5 nT之内。但该流程仍然有诸多可以改进的地方,如日变通化处理过程中运用的三点近似处理方法和长期变化改正过程中对数据的强烈依赖性都是今后值得仔细思考和改进之处。

本文是中国地震局流动地磁团队集体工作成果的阐述,所用流动地磁测量数据为全团队的工作结果,在此向全体团队成员表示感谢。日变通化和长期变化改正所使用的地磁台站数据来自国家地磁台网中心,感谢台网中心和全国地磁台站工作人员对流动地磁团队的长期大力支持。中国地震局流动地磁技术团队希望各位审者和读者给我们提出宝贵的意见和建议,以促进我们团队的健康进步,促进流动地磁监测乃至地磁学科的科学研究水平的提高,促进我国的防震减灾事业长远稳固发展,在此再次感谢。