引用本文
张艳辉, 李世文, 翁爱华, 张素琴, 杨悦, 李建平, 唐裕. 中国异构地磁数据 C-响应估计 [J]. 地震地质, 2019,41(4): 979-995.
ZHANG Yan-hui, LI Shi-wen, WENG Ai-hua, ZHANG Su-qin, YANG Yue, LI Jian-ping, TANG Yu. C-RESPONSES ESTIMATION OF ISOMERIC GEOMAGNETIC DATA IN CHINA [J]. SEISMOLOGY AND GEOLOGY, 2019,41(4): 979-995.  
DOI:10.3969/j.issn.0253-4967.2019.04.011
Permissions
Copyright©2019, 《地震地质》编辑部
《地震地质》编辑部 所有
中国异构地磁数据 C-响应估计
张艳辉1), 李世文1), 翁爱华1), 张素琴2), 杨悦1), 李建平1), 唐裕1)
1) 吉林大学, 地球探测科学与技术学院, 长春 130026
2) 中国地震局地球物理研究所, 北京 100081
*通讯作者: 翁爱华, 男, 教授, E-mail: wengah@jlu.edu.cn。

〔作者简介〕 张艳辉, 男, 1991年生, 2015年于吉林大学地球探测科学与技术学院获应用地球物理专业学士学位, 现为吉林大学地球探测科学与技术学院地球探测信息与技术专业在读博士研究生, 主要从事电磁法及电磁勘探正反演理论研究, E-mail: zhangyh15@mails.jlu.edu.cn。

收稿日期: 2018-06-28
基金项目: 国家重大科研仪器专项(2011YQ05006010)资助
摘要

中国分布着世界上较为密集的地磁台站, 为利用地磁测深获得转换带及下地幔上部高分辨率电性结构奠定了基础, 但相应的 C-响应计算还需要进行系统研究。 文中针对中国地磁数据可能存在的噪声干扰、 记录长度长短不一、 数据记录间断以及数据观测种类不同(绝对和相对)等异构特点, 采用如下技术来提高 C-响应曲线的稳定性: 1)基于BIRRP软件包获得地磁各分量时间序列的频率响应, 使用整体光滑技术来压制数据噪声; 2)对于相对观测台站的数据, 采用近似估计法对水平分量进行分解, 进而估计 C-响应; 3)讨论非连续数据和短记录数据对 C-响应估计的影响; 4)对于沿海台站, 可使用比值法消除海洋效应对 C-响应函数的影响。 使用以上技术手段对中国100余个地磁台站的数据进行处理后, 最终得到42个台站的稳定 C-响应函数, 其中周期为1.3~113.7d的响应个数为24个; 周期为1.3~42.6d的响应个数为18个。 采用文中的技术手段能够很好地处理异构数据, 获得更多、 更稳定的 C-响应函数, 为中国高分辨率地磁测深反演研究提供了更为丰富的基础数据。

关键词: 地磁测深; C-响应; 整体光滑; 异构数据; 海洋效应
中图分类号:P315.72+1 文献标志码:A 文章编号:0253-4967(2019)04-0979-17
C-RESPONSES ESTIMATION OF ISOMERIC GEOMAGNETIC DATA IN CHINA
ZHANG Yan-hui1), LI Shi-wen1), WENG Ai-hua1), ZHANG Su-qin2), YANG Yue1), LI Jian-ping1), TANG Yu1)
1) College of Geo-Exploration Sciences and Technology, Jilin University, Changchun 130026, China
2) Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
Abstract

Geomagnetic depth sounding is an effective method for exploring deep structure of the earth. There are dense geomagnetic observatories in China, which lays a foundation to obtain the electrical structure of the transition zone and the upper part of the lower mantle beneath China. However, the corresponding C-responses estimation methods which are applied now cannot get the stable C-responses for many observatories. Thus, a large amount of geomagnetic data is wasted. Therefore, in order to make full use of the geomagnetic data, the estimation of C-responses needs to be systematically studied. Because of the heterogeneous characteristics of the data quality of China’s geomagnetic observation data, such as the quality of the data, the length of the record, the types of data(absolute and relative observation)and data discontinuity condition, many geomagnetic data are abandoned, this limits the resolution of mantle electrical structure studies. In this paper, the following techniques are used to improve the stability of the data and increase the number of the available geomagnetic observatories, in the meantime, the stability of the C-responses curves can be effectively improved: 1)obtaining the stable spectrums of the different components for each frequency by the BIRRP(Bounded Influence, Remote Reference Processing)software, and using the global smoothing technique to suppress data noise on geomagnetic data; 2)As for the geomagnetic data which only records the relative variation of the D, H and Z components and doesn’t have the baseline value, the horizontal component is decomposed by the approximate estimation method to obtain the C-responses of the relative variation data, and then the relative variation data is used directly for the C-responses estimation; 3)the effects of discontinuous data and short-record data on C-responses estimation are discussed. Under normal conditions, the discontinuity of the data has little influence on C-responses, and when the data length is shorter than 5 years, we can hardly get the available C-responses whose periods are longer than 40 days. All these experiments can provide a basis for the data processing of these kinds of observation data; 4)for coastal observatories, the ratio method is used to eliminate the influence of ocean effect on the C-responses functions. After carefully processing the data of more than 100 geomagnetic observatories in China by the above techniques, the stable C-responses function of 42 observatories is finally obtained, among them, the number of the observatories with C-responses ranging from 1.3 to 113.7 days is 24, and the observatories with C-responses ranging from 1.3 to 42.6 days are 18. The techniques of this paper can process heterogeneous data well and obtain more stable C-responses, which provides more basic data for high-resolution geomagnetic depth sounding inversion researches in China.

Keyword: geomagnetic depth sounding; C-responses; global smoothing; heterogeneous data; ocean effect
0 引言

地磁台站所观测到的磁场变化数据包含了关于地幔过渡带和下地幔电性结构的信息(Armadillo et al., 2001)。 通过对地球内部电性结构的研究可以探索地球的演化历史和动力学特征等问题(Booker et al., 2004; Khan et al., 2006, 2011; Pü the et al., 2015)。

Banks(1969)提出可将地磁测深数据转换为C-响应函数, 并通过反演C-响应函数来研究地球内部电性结构的分布情况(Banks, 1969; Schultz et al., 1987)。 目前, C-响应函数的定义有2种方法, 一种是Z/H方法, 其中垂直分量Z和水平磁场H均由单一台站提供, 该方法基于球谐函数分析求解, 并被应用于研究美国西南地区的导电结构(Schmucker, 1964)。 另一种方法是Z/Y方法(Schmucker, 1979, 1985), 在该定义中, 垂直分量由单一台站提供, 而水平分量由全球台站数据的球谐系数表示(Olsen, 1999)。 随后, Schmucker(1990)和Olsen(1992)使用校正因子对Z/Y方法进行了改进, 使C-响应的求取过程更加合理, 并拓宽了C-响应函数的周期(Olsen, 1992, 1998)。

地磁测深三维反演的分辨率与地磁台站的分布有着密切的关系。 中国有着较为密集的地磁台站, 而在现有的全球三维反演研究中(Kelbert et al., 2009; Semenov et al., 2012; Sun et al., 2015), 绝大部分地磁台站的数据并没有被使用, 这限制了对中国深部精细电性结构的认知。

究其原因, 主要在于中国利用地磁数据进行地球深部电性结构的研究相对较少, 没有形成完整的C-响应处理体系, 从而导致大量的地磁数据尚未应用于地球深部电性结构研究中。 同时, 由于地磁信号相对较弱, 比较容易受到环境和人为噪声干扰, 导致部分台站观测数据的质量较差, 加之近年新建的台站观测时间较短, 部分台站由于种种原因记录不连续。 此外, 据统计, 除国家级台站外, 大部分区域台站测量的地磁场数据类型为不包含基线值的相对变化数据。 对于这样的数据, 常规的C-响应处理方法难以得到理想的响应曲线。 因此, 如何充分利用中国丰富的地磁数据得到更多台站的稳定的C-响应, 是一个值得深入研究的问题。

为此, 本文基于BIRRP(Bounded Influence, Remote Reference Processing)方法(Chave et al., 2004), 将远参考方式改为基于本地台站数据的自参考方式, 求取了中国异构地磁记录的C-响应函数。 文中首先介绍地磁测深C-响应估计的基本原理, 然后讨论了压制数据噪声的整体光滑技术, 之后针对异构的地磁数据讨论了C-响应的求取技术, 最后将上述技术应用于中国地磁观测数据处理中, 得到了42个较为可靠的地磁测深C-响应数据, 为后续的地磁反演等研究奠定了基础。

1 地磁测深C-响应估计的基本原理

在磁层源可以近似表示为 P10的条件下, C-响应的计算公式可以表示为(Banks, 1969)

C(ω)=-a0tanθ2Hr(ω)Hθ(ω)(1)

式中, a0为地球的平均半径; Hr、 Hθ 分别为地表处垂直指向地心和水平N向的磁场; θ 为地磁余纬度; ω 为角频率。

公式(1)表明, 频率域的C-响应由位于某地磁纬度的台站的频率域磁场水平和垂直分量计算得到, 而地磁观测得到的是时间序列, 因此, 若利用实际数据计算地磁测深的C-响应, 大致需要进行如下步骤(Olsen, 1992, 1998; Semenov et al., 2012):

(1)去除原始时间序列中的长期变化场;

(2)由于在地磁坐标下求取C-响应, 故须旋转地磁分量时间序列的坐标系;

(3)将时间序列转换到频率域进行分析, 一般使用离散傅里叶进行变换;

(4)按照式(1)对不同周期的C-响应进行求取;

(5)如开展一维反演研究, 还需要消除海洋效应的影响。

需要说明的是, 由于IAGA(国际地磁学和高层大气物理协会)标准的地磁数据格式包含地磁偏角(D)、 水平分量(H)和垂直分量(Z)的绝对记录或相对记录值, 而式(1)需要使用水平N向磁场和垂直地心的磁场分量, 因此在实际计算前, 需要求取所需的地磁分量的频率响应。 此外, 式(1)是 P10源响应, 为外源场, 而观测的地磁场的主要组成部分为起源于地球内部的长期变化场, 计算时, 首先需要去除内源变化场(Olsen, 1999)。 内源场采用IAGA提供的代码进行求取 。

进行第二步处理的基本原因是由于磁层源与地球的磁赤道同心, 地磁测深数值模拟和反演均在该地磁坐标系中进行, 而实际观测数据大多数基于地理坐标系, 对于这部分数据需要将观测坐标旋转为地磁坐标系, 涉及的主要运算是坐标的欧拉旋转。 下文主要针对异构地磁观测数据, 介绍获取高可靠地磁测深C-响应的估计方法和结果。

2 噪声压制技术
2.1 BI(Bounded Influence)方法

由于噪声的干扰, 直接利用式(1)计算C-响应将影响计算的稳定性和效果。 为此, 与估计大地电磁测深阻抗相同, 本文采用最小二乘方法计算C-响应, 基本的计算公式为

C=-atanθ2< ZH* > < HH* > (2)

式中, H和Z为地磁场水平N分量和垂直分量的频率域响应, 二者均为一维矩阵, H* 为H的复共轭, * 表示矩阵的复共轭, ＜…＞为叠加自(互)功率谱, 表达式为(Huber, 1981)

< AB> =∑lAlBl(3)

A和B代表不同的磁场分量或者参考分量。 Semenov等(2012)通过对比验证发现, 在地磁测深的数据处理过程中, 采用本台站的自参考方法与远参考方法估计的C-响应具有较好的一致性。 因此在后面的计算中, 我们采用自参考法来进行C-响应的求取。

为了获得稳定的磁场分量的频率响应, 采用软件包BIRRP(Chave et al., 2004)计算频谱。 该软件包的核心是BI(Bounded Influence)技术, 其已成功应用于大地电磁测深数据分析和地磁测深C-响应估计等工作中。 其将投影矩阵统计分析和标准M估计技术相结合, 可同时消除Z和H中的干扰及二者中存在的相关噪声, 并且具有较高的计算效率, 故使用其计算磁场的频率响应。

C-响应估计的质量可通过求取相关系数和数据误差进行评价。 在实际进行数据反演时, 将相关系数较低或数据误差较大的C-响应数据舍弃。

C-响应不同周期的平方相关系数的计算公式为(Semenov et al., 2012)

coh2=|< ZH* > |2< ZZ* > < HH* > (4)

误差计算参考文献(Schmucker, 1999), 具体为

δC(ω)=|C|1-coh2coh21β1L-1(5)

设1-β 是置信水平, 与Semenov等(2012)的研究相同, 本文设置信水平为0.9, 因此β 为0.1。 L是不同周期计算时对应的时间序列的分段数目, 周期越大, L值越小。

Semenov等(2012)对ASP(Alice Springs)和HON(Honolulu)台站1957— 2007年的地磁数据 求取了C-响应, 为了验证BI方法求取地磁测深C-响应的正确性, 我们对这2个地磁台站的数据进行了重新处理, 由于数据资源的限制, 只能收集到1992— 2017年的地磁数据, 比原始的时间序列更短, 得到的C-响应曲线如图 1所示。 由图可见, 即使在更短时间序列的前提下, 通过本文方法处理得到的C-响应曲线在短周期方面与前人的结果保持着较好的一致性, 且在长周期方面更加稳定, “ 飞点” 情况明显减少。 从平方相关系数和误差棒也能明显看出该方法求得的C-响应相关系数更高、 误差棒更短, 从另外一个侧面反映了该方法的正确性和稳定性。

图 1 利用1992— 2017年地磁观测时间序列通过本文方法求取的C-响应与Semenov等(2012)利用1957— 2007年地磁观测时间序列计算结果(SK2012)的对比
上部为平方一致性曲线, 误差棒表示响应函数的可靠程度Fig. 1 Comparison of the C-responses obtained by this method with the results of Semenov et al., 2012.

2.2 整体光滑约束技术

在准静态近似条件下, 地磁测深C-响应具有体积效应(李世文等, 2018), 因此频率响应曲线应该是稳定光滑的。 图 2给出了根据Pü the等(2015)推导出的全球电导率模型计算的理论C-响应曲线。 C-响应的计算采用李世文等(2017)提出的方法。 由图可见, 理论C-响应曲线为复数, 包括实部和虚部, 并且其频率响应曲线平稳而光滑。

图 2 Pü the等(2015)由卫星数据计算的全球平均一维导电模型(a)及其C-响应曲线(b)
实线为C-响应的实部, 虚线为C-响应的虚部Fig. 2 One￣dimensional global average conduction model from Pü the et al., 2015(a) and its C-responses curve(b).

受到环境噪声的干扰导致地磁台站记录的数据质量较差时, 求取的响应曲线将存在上下震荡和 “ 飞点” 现象(Pü the et al., 2014), 这显然不符合理论数据的变化规律。 因此, 我们提出将不同频率的C-响应作为一个整体进行估计, 并要求其满足光滑条件, 即在式(1)中令ω i=i, C(ω i)=Ci, Hr(ω i)=Vi , - 2α0tanθHθ (ω i)=Hi, 则式(1)可以被表示为

Vi=CiHi, i=1, 2, …, N(6)

其中, N为频率个数。

考虑到C-响应的光滑特性, 对其进行光滑约束后, 将式(6)关于C-响应的求取问题转化为约束优化问题, 即:

min(ФC)=‖WC‖2(7)

式(7)满足:

HC=V(8)

其中, W是模型约束矩阵, 可选用一阶或二阶光滑矩阵(Constable et al., 1987)。 式(7)和式(8)可转化为正则化反演问题, 其目标函数为

Ф=Фd+λФC(9)

其中, λ 是平衡数据误差和模型光滑度的正则化参数(Tikhonov, 1963; 王家映, 2002), Ф C为模型光滑项; 数据预测误差Ф d被定义为

Фd=‖V-HC‖2(10)

式(9)取极小值时, ∂Ф∂C=0, 有:

⌊H* TH+λWTW」C=H* TV(11)

在式(11)中, 通过L曲线法选取最优正则化因子λ (Hansen, 1992), 从而可以得到经光滑约束后的C-响应。 由于该方法基于正则化反演(Regularized Inversion), 因此下文简称其为RI方法。 该方法与以往方法的本质区别是将整个频率范围作为一个整体求解, 并加入了光滑约束, 张艳辉等(2019)关于该方法有更加详细的讨论。

3 异构数据的处理

地磁观测数据中广泛存在着仪器和环境噪声干扰较大、 记录长度较短及数据记录不同程度间断的现象, 此外中国部分地区地磁台站记录数据为磁场各分量的相对变化, 我们将这些数据称为异构数据, 异构数据的存在严重影响着C-响应的求取。 本文采用一系列技术手段对含噪声数据、 相对数据、 短记录数据以及非连续数据等异构数据进行了测试, 同时对沿海台站的理论数据进行了海洋校正, 验证了使用异构数据求取稳定C-响应的可行性。

3.1 数据干扰的压制

大量的地磁台站遭受着越来越严重的环境噪声的干扰, 利用整体光滑约束方法求取C-响应能够更好地抵制噪声的干扰, 得到更真实的结果。 图 3给出了利用长春(CNH)和满洲里(MZL)台的数据计算得到的C-响应结果。 通过C-响应和平方相关系数曲线可见, 通过常规方法(BI)处理得到的C-响应曲线在长周期时存在明显的跳点, 而采用整体光滑技术(RI)得到的C-响应曲线则更加平稳、 光滑, 减少了数据噪声对结果的影响, 该方法对提高C-响应的质量和台站数据的利用率有着很重要的作用。

图 3 整体光滑约束技术处理结果(RI)和BIRRP方法处理(BI)结果的对比
上部为平方一致性曲线, 误差棒可反映C-响应函数的可靠程度Fig. 3 Comparison between the results of the global smooth constraint technique and the conventional processing.

3.2 相对数据

地磁台站记录的时间序列分为各分量绝对观测值(D、 H、 Z)和相对变化值(Δ D、 Δ H、 Δ Z), 前者为相对变化值加上各分量的基线值。 常规地磁数据处理方法认为对于不含基线值(即仅有相对变化值)的台站数据无法进行处理(Semenov et al., 2012), 从而造成了大量相对观测台站的数据被浪费。 本文尝试对地磁相对数据进行处理, 下文将详细介绍其原理。

式(1)中, Hr(ω )=F(Z), H0(ω )=F(X), 假设地磁场X分量的基线值为X0, 则X0+Δ X=X; 同时, 由傅里叶变换的线性性质可得:

FX0+ΔX=X0F[1]+F[ΔX]=2πδ(ω)X0+F[ΔX](12)

其中, δ (ω )为狄拉克函数, 在除零频之外的所有频率上的取值均为零。 从式(12)可见, 时间域的基线偏移在频率域中表现为零频的脉冲函数, 而由于C-响应研究中不考虑零频分量, 因此在进行地磁数据的时频转换时, 可以直接使用Δ X、Δ Y和Δ Z来代替X、 Y和Z分量进行C-响应估计。

此外, 经过数据分析发现, 地磁场水平分量变化以水平N分量为主, 因此水平分量总变化量与水平N分量的变化幅度相当, 故可采用近似法对相对数据进行处理(王桥, 2015), 即:

ΔX≈ΔH(13)

这里假设基线偏移在整个时间序列中是相同的。 如果基线值出现差异, 在进行资料处理时, 可通过人为加以调平, 保证整个时间序列的基线相同。

为了验证这种近似方法得到的C-响应的正确性, 本文对长春台站(CNH)和满洲里台(MZL)的绝对数据和相对数据(由绝对数据减去基线值得到)分别进行了常规处理和近似处理, 在得到了各个分量的时间序列后, 通过上述过程得到了对应的C-响应曲线(图4)。 由图可见, 采用绝对数据和相对数据得到的响应曲线几乎一致, 只在长周期时的个别频点存在细微的差别, 说明近似处理方法具有可行性。

图 4 使用绝对数据和不包含基线值的相对数据的C-响应处理结果对比
上部为平方一致性曲线, 误差棒可反映C-响应函数的可靠程度Fig. 4 Comparison of C-responses of the different data type.

3.3 短记录数据

陆续增加的新地磁观测台站的记录时间长度长短不一, 最短的只有5a。 本文使用满洲里台(MZL)1995— 2015年(2003、 2004年数据缺失)的数据, 人为缩短观测记录的时间长度, 以测试记录时间长度对C-响应值的影响, 所得结果如图 5所示。 由图可见, 记录时间越长, 对于每个周期可叠加的时间序列越多, 得到的C-响应越稳定; 对于相同的时间序列, 短周期可获得的叠加次数更多, 因此响应曲线更加稳定。 通过C-响应和平方一致性曲线可知, 连续记录时间> 15a时, 由于各个周期可叠加的时间序列足够多(> 40次), 因此可以得到较稳定的C-响应; 当连续记录时间< 10a或更短时, 长周期C-响应将出现明显的跳点, 而周期< 40d 的C-响应相对稳定。 这是由于短周期所需的时间序列叠加次数更多, 得到的响应曲线更稳定。 此外, 当连续记录时间< 5a时, 由于周期为40d 的C-响应所需要的时间序列少于40次叠加, 故无法得到稳定的响应曲线。

图 5 满洲里台站数据记录长度影响分析
a 平方一致性曲线; b 不同记录长度C-响应结果; c 数据记录长度示意图Fig. 5 Impact of data record length on C-response estimation at Manzhouli observatory.

因此, 在实际数据处理时, 短记录数据的短周期C-响应依然可被用于浅部地幔结构探测, 并约束长周期数据的反演, 从而提高可利用台站的密度。

3.4 非连续数据

由于种种原因, 大部分地磁台站的观测记录将出现不同程度的中断或数据丢失。 为测试记录不连续对C-响应值的影响, 同样使用满洲里台的数据, 经人为剪切部分观测数据之后进行计算, 所得结果如图 6所示。 由图可见, 数据的连续性越好, 得到的C-响应越稳定, 从C-响应和平方一致性曲线的趋势来看, 数据连续性较差对长周期的C-响应有一定的影响, 但主要在周期> 60d之后。 这是由于在同样长度的时间序列中, 短周期响应所需要的时窗长度更短, 因此短周期函数的叠加次数更多, 故对于相同的时间序列, 短周期的C-响应更加稳定。

图 6 数据中断影响分析
a 平方一致性曲线; b 不同中断记录的C-响应结果; c 中断程度示意图Fig. 6 Impact of data interruption on the estimation of C-response at MZL observatory.

当然, 当时间序列被切割的非常分散, 无法保证连续记录的长度满足该周期C-响应所需的时窗长度时, 显然是无法得到响应曲线的。 经过分析可知, 在通常存在数据中断及常规噪声干扰的情况下, 获得稳定C-响应的前提是需要完整的时窗长度进行至少40次叠加。 由于地磁观测台站每日都有数据质量监督, 不会出现大规模的数据中断(图6), 因此, 在实际数据处理中, 部分数据中断并不会导致求取的C-响应曲线出现较大的偏差, 可以将包含中断记录的数据加以利用, 从而提高数据的利用率。

3.5 海洋效应校正

地球表面约70%被高导海水覆盖, 海陆间明显的电导率差异将产生海洋效应(geomagnetic coast effect, GCE)(Ribaudo, 2011)。 进行三维正、 反演时可直接加入地球表层模型以消除海洋效应, 而在一维正、 反演时, 海洋效应的存在将影响地磁资料的解释(Parkinson et al., 1979; Fainberg et al., 1980)。

海洋效应的影响可以通过比值法进行校正(Kuvshinov et al., 2002), 具体公式为

Cobs, corr=Cobs·k=Cobs·C1DC1D+3Dsurface(14)

其中, Cobs, corr是校正后的C-响应, Cobs是观测数据, k称为校正系数, C1D是一维模型响应, 一维模型为中国平均电导率模型(李建平等, 2018); C1D+3D surface是包含地表电导分布的三维模型响应, 采用Kelbert等(2014)的表层电导率地图(S-map)及交错网格有限差分进行数值计算(李建平等, 2018)。

图 7为杭州台(HAZ)和大连台(DAL)全球平均理论模型的海洋效应影响结果。 由图可见, 在短周期时, 海洋效应对C-响应的影响较强, 且对不同台站的影响程度明显不同, 因此, 在进行一维正、 反演研究时, 对于沿海台站的海洋效应校正是很有必要的。 李建平等(2018)针对中国所有台站的模拟结果表明, 中国西部地区以及东部沿海台站的长周期C-响应曲线几乎没有受到海洋效应的影响, 而如图 8中黑色实线所界定的东部地区沿海台站的短周期C-响应则会受到海洋效应的影响。

图 7 基于Pü the等(2015)提出的全球平均模型下, 杭州台(HAZ)和大连台(DAL)海洋效应影响分析
一维模型响应为全球一维平均模型计算的正演响应, 三维模型响应为一维平均模型加上海洋-陆地电阻率模型后计算的正演响应(即消除了海洋效应的响应曲线)Fig. 7 Analysis of the impact of ocean effect at Hangzhou(HAZ) and Dalian(DAL)observatories.

图 8 中国地磁台站分布图
圆形为绝对观测台站; 三角形为相对观测台站; 实心图形为通过本文处理可得到稳定C-响应的台站。
对于黑色实线右侧的台站, 海洋效应对其短周期C-响应(1.3d)的影响幅度> 20%Fig. 8 Distribution map of geomagnetic station in China.

4 处理结果

中国分布着110余个地磁台站(图8), 但中东部台站相对密集, 在东北和西部地区较为稀疏, 且其中绝对测量台站有30余个, 其余为相对测量台站。 显然, 如果可以充分利用如此密集的台站所收集的地磁数据, 能够明显提高地磁测深三维反演的分辨率。

对于绝对测量台站的数据, 利用前文的讨论内容和相关技术, 经过仔细处理并分析后, 认为其中24个台站的C-响应曲线形态良好, 可以用于地磁测深研究。 表1和图 9给出了具体的台站和其相应的响应曲线。

表1 本文所使用的绝对测量地磁台站的名称、 地理与地磁经纬度以及时间长度信息 Table1 List of observatories with absolute observation data used in this paper, including station name, longitude and latitude of geographic and geomagnetic coordinates, and record length information

图 9 未经过海洋效应校正的绝对观测数据C-响应处理结果Fig. 9 C-responses curves of abosolute observation geomagnetic data without ocean effect correction.

对仅测量相对数据的二级地磁台站进行处理后, 最终得到18个台站的C-响应曲线, 具体台站见表2, 响应曲线见图10。 由于绝大部分二级台站为2008年以后才开始运行的新建台站, 记录时间较短, 有的甚至不足5a, 故在求取C-响应时将有效周期的范围缩小为1.33~42.67d。

表2 中国18个相对测量地磁台站的名称、 地理与地磁经纬度以及时间长度信息 Table2 List of 18 observatories with relative observation data used in this paper, including station name, longitude and latitude of geographic and geomagnetic coordinates, and record length information

图 10 未经过海洋效应校正的相对数据C-响应处理结果Fig. 10 C-responses curves of relative observation geomagnetic data without ocean effect correction.

其中, 受海洋校正影响的台站范围见图8。 对于该分界线以东区域的地磁台站, 海洋效应将导致C-响应曲线的短周期(1.3d)出现明显偏差, 因此在进行一维反演时需通过比值法进行海洋效应校正, 而在三维反演时可以加入地壳表层电导率模型来消除海洋效应的影响。

5 结论

为了更好地求取异构数据的响应曲线, 提高C-响应估计质量, 本文提出了一系列技术手段以提高地磁观测数据利用率。

(1)提出基于BI(Bounded Influence)方法的全局光滑约束技术以求取C-响应, 由此利用数据质量较差的台站依然能够获得较稳定响应曲线。

(2)对于相对观测数据, 本文使用近似方法来求取相对观测数据的C-响应, 处理得到中国18个相对观测地磁台的C-响应曲线, 为今后的研究工作提供了更多的基础数据。

(3)非连续数据和短记录数据对求取C-响应影响的测试结果表明, 非连续和记录时间> 10a的数据可以得到周期达到113.7d的C-响应曲线, < 10a且> 5a的短记录数据可通过处理得到周期达到42.6d的稳定的C-响应结果。

利用本文总结的技术对中国地磁台站观测数据进行处理, 将可获得可靠的C-响应曲线的台站数增加到了42个, 为开展中国高分辨率的地磁测深反演研究提供了更丰富的基础数据。

致谢 Chave A D 提供了BIRRP软件包; 中国地震局国家地磁台网中心提供了地磁数据; 审稿人对本文提出了宝贵的修改意见, 使本文的研究更加充实、 完整。 在此一并表示感谢!

参考文献
[1] 李建平, 翁爱华, 李世文, 等. 2018. 海洋效应对中国沿海地磁观测影响: 以广州台站为例[J]. 地球物理学报, 61(2): 649—658.
LI Jian-ping, WENG Ai-hua, LI Shi-wen, et al. 2018. The influence of ocean effect on geomagnetic observations in coastal areas of China: A case study of the Guangzhou observatory[J]. Chinese Journal of Geophysics, 61(2): 649—658(in Chinese). [本文引用:3]
[2] 李世文, 翁爱华, 李建平, 等. 2017. 浅部约束的地磁测深 C-响应一维反演[J]. 地球物理学报, 60(3): 1201—1210.
LI Shi-wen, WENG Ai-hua, LI Jian-ping, et al. 2017. 1-D inversion of C-response data from geomagnetic depth sounding with shallow resistivity constraint[J]. Chinese Journal of Geophysics, 60(3): 1201—1210(in Chinese). [本文引用:1]
[3] 李世文, 翁爱华, 唐裕, 等. 2018. 一维导电薄球层状模型的地磁测深 C-响应计算[J]. 地震地质, 40(2): 337—348. doi: 10.3969/j. issn. 0253-4967. 2018. 02. 004.
LI Shi-wen, WENG Ai-hua, TANG Yu, et al. 2018. C-response of geomagnetic depth sounding on a 1-D thin shell model[J]. Seismology and Geology, 40(2): 337—348(in Chinese). [本文引用:1]
[4] 王家映. 2002. 地球物理反演理论 [M]. 北京: 高等教育出版社.
WANG Jia-ying. 2002. Geophysical Inversion Theory [M]. Higher Education Press, Beijing(in Chinese). [本文引用:1]
[5] 王桥. 2015. 中国中东部地磁感应矢量研究 [D]. 北京: 北京大学.
WANG Qiao. 2015. Studies on the geomagnetic induction vectors of Mid-east China [D]. Peking University, Beijing(in Chinese). [本文引用:1]
[6] 张艳辉, 翁爱华, 李世文, 等. 2019. 基于全局光滑约束的地磁测深 C-响应估计[J]. 地球物理学报, 62(5): 1898—1907.
ZHANG Yan-hui, WENG Ai-hua, LI Shi-wen, et al. 2019. Estimation of C-responses of geomagnetic depth sounding based on global smooth constraint[J]. Chinese Journal of Geophysics, 62(5): 1898—1907(in Chinese). [本文引用:1]
[7] Armadillo E, Bozzo E, Cerv V, et al. 2001. Geomagnetic depth sounding in the northern Apennines(Italy)[J]. Earth, Planets and Space, 53(5): 385—396. [本文引用:1]
[8] Banks R J. 1969. Geomagnetic variations and the electrical conductivity of the upper mantle[J]. Geophysical Journal International, 17(5): 457—487. [本文引用:2]
[9] Booker J R, Favetto A, Pomposiello M C. 2004. Low electrical resistivity associated with plunging of the Nazca flat slab beneath Argentina[J]. Nature, 429(6990): 399—403. [本文引用:1]
[10] Chave A D, Thomson D J. 2004. Bounded influence magnetotelluric response function estimation[J]. Geophysical Journal of the Royal Astronomical Society, 157(3): 988—1006. [本文引用:2]
[11] Constable S C, Parker R L, Constable C G. 1987. Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data[J]. Geophysics, 52(3): 289—300. [本文引用:1]
[12] Fainberg E B, Zinger B S. 1980. Electromagnetic induction in a non-uniform spherical model of the Earth[J]. Annales de Geophysique, 36(2): 127—134. [本文引用:1]
[13] Hansen P C. 1992. Analysis of discrete ill-posed problems by means of the L-curve[J]. SIAM Review, 34(4): 561—580. [本文引用:1]
[14] Huber P. 1981. Robust Statistics[M]. Wiley, New York, US. [本文引用:1]
[15] Kelbert A, Kuvshinov A, Velímsky J, et al. 2014. Global 3-D electromagnetic forward modelling: A benchmark study[J]. Geophysical Journal International, 197(2): 785—814. [本文引用:1]
[16] Kelbert A, Schultz A, Egbert G. 2009. Global electromagnetic induction constraints on transition-zone water content variations[J]. Nature, 460(7258): 1003—1006. [本文引用:1]
[17] Khan A, Connolly J A D, Olsen N. 2006. Constraining the composition and thermal state of the mantle beneath Europe from inversion of long-period electromagnetic sounding data[J]. Journal of Geophysical Research: Solid Earth, 111(B10): 207—208. [本文引用:1]
[18] Khan A, Kuvshinov A, Semenov A. 2011. On the heterogeneous electrical conductivity structure of the earth’s mantle with implications for transition zone water content[J]. Journal of Geophysical Research: Solid Earth, 116(B1): B01103. [本文引用:1]
[19] Kuvshinov A V, Avdeev D B, Pankratov O V, et al. 2002. Chapter 3: Modelling electromagnetic fields in a 3D spherical earth using a fast integral equation approach[J]. Methods in Geochemistry & Geophysics, 35: 43—54. [本文引用:1]
[20] Olsen N. 1992. Day-to-day C-response estimation for Sq from 1 cpd to 6 cpd using the Z: Y method[J]. Journal of Geomagnetism and Geoelectricity, 44(6): 433—447. [本文引用:2]
[21] Olsen N. 1998. The electrical conductivity of the mantle beneath Europe derived from the C-responses from 3 to 720 hr[J]. Geophysical Journal International, 133(2): 298—308. [本文引用:2]
[22] Olsen N. 1999. Long-period(30 days -1 year)electromagnetic sounding and the electrical conductivity of the lower mantle beneath Europe[J]. Geophysical Journal International, 138(1): 179—187. [本文引用:2]
[23] Parkinson W D, Jones F W. 1979. The geomagnetic coast effect[J]. Reviews of Geophysics, 17(8): 1999—2015. [本文引用:1]
[24] Püthe C, Kuvshinov A. 2014. Mapping 3-D mantle electrical conductivity from space: A new 3-D inversion scheme based on analysis of matrix q-responses[J]. Geophysical Journal International, 197(2): 768—784. [本文引用:1]
[25] Püthe C, Kuvshinov A, Khan A, et al. 2015. A new model of earth’s radial conductivity structure derived from over 10 yr of satellite and observatory magnetic data[J]. Geophysical Journal International, 203(3): 1864—1872. [本文引用:1]
[26] Ribaudo J T. 2011. Flexible finite-element modeling of global geomagnetic depth sounding [D]. San Diego: University of California. [本文引用:1]
[27] Schmucker U. 1964. Anomalies of geomagnetic variations in the southwestern United States[J]. Journal of Geomagnetism and Geoelectricity, 15(4): 193—221. [本文引用:1]
[28] Schmucker U. 1979. Erdmagnetische Variationen und die elektrische Leitfähigkeit in tieferen Schichten der Erde[J]. Sitzungsbericht und Mitteilungen Braunschweigische Wiss. Gesellschaft, Sonderheft, 4: 45—102. [本文引用:1]
[29] Schmucker U. 1985. Magnetic and electric fields due to electromagnetic induction by external sources, electrical properties of the earth’s interior[J]. Land olt-Boernstein, New-Series, 5/2b. Berlin: Springer-Verlag. [本文引用:1]
[30] Schmucker U. 1990. Die eindringtiefen tagesperiodischer variationen [A]∥Haak V, Homilius, J eds. Protokoll Koll. Elektromagnetische Tiefenforschung, Hornburg. Niedersächsisches Land esamt, Hannover: 31—66. [本文引用:1]
[31] Schmucker U. 1999. A spherical harmonic analysis of solar daily variations in the years 1964—1965: Response estimates and source fields for global induction I. Methods[J]. Geophysical Journal International, 136: 439—454. [本文引用:1]
[32] Schultz A, Larsen J C. 1987. On the electrical conductivity of the mid-mantle: I. Calculation of equivalent scalar MT response functions[J]. Geophysical Journal International, 88(3): 733—761. [本文引用:1]
[33] Semenov A, Kuvshinov A. 2012. Global 3-D imaging of mantle conductivity based on inversion of observatory C-responses Ⅱ. Data analysis and results[J]. Geophysical Journal International, 191(3): 965—992. [本文引用:7]
[34] Sun J, Kelbert A, Egbert G D. 2015. Ionospheric current source modeling and global geomagnetic induction using ground geomagnetic observatory data[J]. Journal of Geophysical Research: Solid Earth, 120(10): 6771—6796. [本文引用:1]
[35] Tikhonov A N. 1963. Regularization of incorrectly posed problems[J]. Soviet Mathematics Doklady, 4(1): 1624—1627. [本文引用:1]