2019, Vol. 62
2. 中国地质调查局广州海洋地质调查局, 广州 510075;
3. 吉林省有色金属地质勘查局, 长春 130021
, WENG AiHua1
, ZHANG YanHui1, LI JianPing2, YANG Yue1, TANG Yu1, ZOU ZongLin1, LI ChunCheng3
2. Guangzhou Marine Geological Survey, China Geological Survey, Guangzhou 510075, China;
3. Jilin Nonferrous Metal Geological Exploration Bureau, Changchun 130021, China
地幔导电性与地幔物质成分、温度、水含量、熔融状态等密切相关(Karato，1990；Xu et al., 2000；Huang et al., 2005；Yoshino et al., 2008；Manthilake et al., 2016)，通过对地球内部电性结构的研究可以探索地球的演化历史、动力学特征等问题(Fukao et al., 2004；Booker et al., 2004；Khan et al., 2006, 2011；Püthe et al., 2015).
区域性和全球范围的一维地磁测深数据反演结果显示，地幔的导电性存在明显横向不均匀性(Peyronneau and Poirier, 1989；Olsen，1998；赵国泽等，2001；Semenov and Jozwiak, 2006；汤吉等，2006；Püthe et al., 2015；徐光晶等，2015；王桥和黄清华，2016；Munch et al., 2018)；最近的全球地磁台站数据三维反演结果进一步揭示在地幔中存在多处电性异常区域，比如中国东北地区存在明显的高导体，南欧和北非地区存在明显的高阻体(Kelbert et al., 2009；Semenov and Kuvshinov, 2012；Sun et al., 2015)等.这些结果对认识地幔转换带及下地幔上部的热结构、板块俯冲过程和地幔水循环等起到了重要作用.
地磁测深反演，一般是认为在磁层源可以由单一球谐函数P10近似表示的前提下，对地磁台站的时间序列数据进行处理并计算得到C-响应(Banks, 1969)；之后通过求解基于L2-范数的目标函数，得到一个既满足光滑渐变又较好拟合观测数据的模型(Egbert and Kelbert, 2012).反演中假设数据噪声和模型参数均符合高斯分布.然而，在对全球中低纬度地区地磁台站数据进行处理后，我们发现，大部分台站的C-响应曲线相对光滑；但也有一些台站，由于存在数据记录不连续、信号所受干扰大、连续观测时间短等问题，估计得到的C-响应曲线振荡，存在明显飞值，此时数据噪声不符合高斯分布(张艳辉等，2019).如果继续选择以L2-范数度量数据预测误差，反演将会受到这些异常数据的影响，导致反演结果的可靠性降低(Menke，1984；Oldenburg，1984；Farquharson，2008).
如果能在地磁测深三维反演中充分利用目前分布十分稀疏的地磁台站资料，保留尽可能多的观测信息，增加可利用的台站分布密度，则可以提高反演对地幔电性结构的分辨能力(Kelbert et al., 2008).因此，对于C-响应数据不稳定的台站，如何有效降低“飞点”数据对反演可靠性的影响，尽可能多的保留数据有效信息，是地磁测深反演需要研究的一个重要问题.
Menke(1984)及Oldenburg(1984)等指出，如果采用L1-范数对数据拟合差进行度量，可以明显压制“飞点”数据对反演结果的影响，得到相对可靠的反演模型.Farquharson和Oldenburg(1998)利用数值试验，对数据残差项分别按照L1-范数和L2-范数度量进行一维反演，验证了L1-范数在压制“飞点”数据影响方面有较明显的优势.
本文将L1-范数应用于地磁测深C-响应数据的三维反演.通常认为地幔的导电性横向上在一个数量级以内缓慢变化(Kelbert et al., 2009；Semenov and Kuvshinov, 2012)，符合L2-范数反演的基本假设，故本文在构建反演的目标函数时，对模型粗糙度仍采用L2-范数进行度量，而只对稳定性存在较大差别的C-响应数据采用L1-范数度量数据拟合(后文均称为L1-范数反演).为模拟具有“飞点”的稳定性较差C-响应数据，我们按照指数概率密度分布函数计算随机数据噪声.本文首先介绍地磁测深C-响应的数值模拟和数据处理方法，并介绍了基于有限内存拟牛顿法的L1-范数反演理论；然后通过反演叠加指数分布噪声的合成数据，研究L1-范数的反演效果；最后，将L1-范数反演应用于全球实际数据的三维反演，得到新的地球电性结构模型.
1 地磁测深C-响应 1.1 C-响应的定义为了进行地磁数据的反演，我们通常将地磁台站记录的时间序列信号转换为C-响应或其他响应函数.本文研究基于广泛采用的C-响应函数展开，在磁层源可以近似表示为P10的条件下，C-响应的计算公式可以表示为
|
(1) |
式中a0为地球的平均半径；Hr、Hθ分别为地表处垂直指向地心和水平北向的磁场；θ为地磁余纬度；ω为角频率.C-响应的单位为km，一般情况下实部为正值，可以作为电磁场穿透深度的近似值；虚部为负值，其幅值在一定程度上反映了电阻率的变化.C-响应的周期范围一般为n~102天，可用于研究地下400~1600 km深度的电性结构(Schmucker，1987；Kelbert et al., 2008；Püthe et al., 2015).
1.2 C-响应的数值模拟本文中任意模型的C-响应计算采用Uyeshima和Schultz(2000)提出的交错网格有限差分法对式(2)进行数值求解.取正时谐的条件下，关于磁场的亥姆霍兹方程为
|
(2) |
式中H为磁场，ρ为电阻率，ω为频率，μ0为真空磁导率，i为虚数单位.(2)式的计算域包括高阻空气层、地壳和地幔.计算中空气层的电阻率取1010Ωm，厚度选择为2倍地球半径，以保证上边界的感应磁场已衰减至足够小；下边界为核幔边界，并认为此处的磁场切向分量为0.与直角坐标系中的网格单元不同，地磁测深C-响应的数值模拟在球坐标系中进行，对应的网格单元为曲边六面体，并将磁场定义在网格的棱边上.求解(2)式得到磁场分量后，通过(1)式和位置插值函数即可得到地表任意位置测点的C-响应.
在交错网格有限差分中，采用迭代求解的方法(Kelbert et al., 2008)计算(2)式，并进行磁场的散度校正以加快收敛并保证计算精度.三维数值模拟的详细介绍参考Uyeshima和Schultz(2000)及李建平等(2018a).
1.3 观测数据的C-响应计算进行地磁数据处理时，我们首先剔除原始时间序列数据中存在记录错误的点，然后基于国际地磁参考场移除地球磁场的长期变化.由于地磁测深的正演计算是在地磁坐标系中进行的，因此需要将观测数据由地理坐标系旋转到地磁坐标系中，得到对应的磁场三分量.之后，本文基于Chave和Thomson(2004)的BIRRP数据处理软件，采用自参考方法(Semenov and Kuvshinov, 2012)替代远参考，计算各台站周期为3~113天的C-响应，并选择0.95的置信度求取平方相关系数和数据误差.
图 1给出了国际地磁台站FUR(Fürstenfeldbrück)和LNP(仑坪)的原始观测数据曲线及处理得到的C-响应和0.95置信度对应的平方相关系数(squared coherency/coh2)曲线.一般来说，平方相关系数越大，数据质量越好(Chave and Thomson, 2004).从图 1可见FUR的观测连续，记录时间较长，达到了65年，数据质量也较好，磁场的X和Z分量没有明显的异常值；数据处理得到的各周期响应曲线较为光滑，误差棒短，平方相关系数值大，噪声符合高斯统计规律.LNP的磁场各分量虽然也较稳定，但是观测存在多个间断，总记录时长不足20年；数据处理得到的响应值稳定性较差，误差棒较长，平方相关系数值小，尤其是大周期的响应，受到数据时长的制约，波动较大，其噪声为非高斯分布.
|
图 1 国际地磁台站FUR(a、c)和LNP(b、d)的原始观测数据及处理得到的C-响应 原始数据记录长度：FUR(1950—2015年)、LNP(1965—1969、1980—1984及1989—2000年). Fig. 1 Observed data in time-series and estimated C-responses of FUR (a, c) and LNP (b, d) Length of records: FUR (1950—2015), LNP (1965—1969, 1980—1984 and 1989—2000). |
表 1给出了参与实际数据反演的129个位于中低纬度的台站各周期平方相关系数的平均值，大部分台站的平均系数大于0.5；但其中64个台站的系数平均值小于0.5，表明C-响应曲线稳定性较差，尤其在大周期存在明显“飞点”，数据误差为非高斯噪声分布，不符合L2-范数反演的基本假设.地磁台站分布原本就十分稀疏，如果将相关系数相对较小但部分频率响应合理的台站直接全部丢弃，显然不利于提高对地幔电性结构的反演分辨能力和可靠性.因此，本文提出采用L1-范数度量数据拟合项的方法，有效降低“飞点”数据的影响，以在反演中保留尽可能多的测点数据.
|
表 1 本文选取地磁台站的位置信息及平方相关系数(coh2) Table 1 Location details and squared coherencies (coh2) of observatories used in this paper |
地磁测深三维反演问题通常表示为寻找目标函数
|
(3) |
的最优解，其中Φ(m, λ)一般定义为
|
(4) |
式中的Φd(m)和Φm(m)分别为数据残差项和模型粗糙度项；λ为正则化参数，用来平衡Φd(m)和Φm(m)的权重.这里，m不是Kelbert等(2008, 2009)描述地球内部导电性分布的球谐系数矢量，而是表示球坐标系中地球内部各剖分网格单元的平均电导率分布.采用这种网格化模型参数的目的是提高反演模型的自由度，避免模型复杂度受球谐函数系数最高阶数控制的问题.
如前文所述，在反演中，我们对数据拟合项采用L1-范数的形式度量，对模型粗糙度项依然采用L2-范数的形式度量，那么目标函数可以写为
|
(5) |
式中d为观测数据向量；m和m0分别为模型参数和初始模型参数向量；ψ(m)为模型的正演响应；Cd-1和Cm-1分别为数据和模型的协方差算子.根据Kelbert等(2008)及Egbert和Kelbert(2012)，引入新参数
|
(6) |
由L1-范数定义的数据拟合项
|
(7) |
式中，Cmodel(ω)和Cobs(ω)分别为模型的响应数据和合成/观测数据，δC(ω)为数据误差.为了避免(7)式及其导数在L1-范数反演求解中可能出现的零值和奇点，本文采用Ekblom(1988)提出的扰动Lp-范数定义式
|
(8) |
这里，φ(x)为通式，表示Φd或者Φm求和计算中的某一项，本文只对Φd进行L1-范数的计算；ε为一很小的值，以确保在p=1，x=0时L1-范数度量函数的偏导数存在.
2.2 有限内存拟牛顿法反演为求解目标函数(4)式的最优化解，我们选择有限内存拟牛顿法(L-BFGS)进行反演(Avdeev and Avdeeva, 2009；刘云鹤和殷长春，2013；翁爱华等，2015).
有限内存拟牛顿法由牛顿法演化而来，其迭代基本格式为
|
(9) |
|
(10) |
|
(11) |
式中，k为迭代次数；αk是搜索步长；T为矩阵的转置运算；pk是搜索方向；Hk-1是海森矩阵的逆.在根据(9)式求解搜索方向pk时，利用了目标函数的二阶导数信息，因此能快速达到收敛条件.由式(5)可知，目标函数对当前迭代模型的梯度∇Φk包括∇Φd和∇Φm两部分.在数值计算时，∇Φm项有解析表达式，容易计算；∇Φd的求解则采用伴随矩阵的方法(参考2.3节中介绍).
牛顿法在求解搜索方向时需要计算和存储海森矩阵并求其逆，对存储空间和计算时间有很高的要求.拟牛顿法改进了这一点，它通过构造一个对称正定矩阵Bk来近似表示海森矩阵(Nocedal and Wright, 1999)，从而(10)式可写成
|
(12) |
为了使Bk合理的逼近Hk，要求其满足
|
(13) |
式中sk=mk+1-mk，yk=∇Φk+1-∇Φk.通常给定一个初始单位矩阵B0，然后通过修改Bk得到Bk+1.显然Bk+1依赖于之前各次迭代时的Bi(i=1, …, k)，因此对存储空间要求较高，并且耗时较长.有限内存拟牛顿法在拟牛顿法的基础上，由一个对称正定初始矩阵B0∈Rn×n开始，Bk+1利用下面的校正公式
|
(14) |
计算.反演中，预先确定正整数M，从而只利用最近M次的迭代信息计算Bk+1，极大的减少了内存需求，节省了时间，提高了计算速度.
2.3 目标函数的梯度通过2.1节的分析，我们可以得到目标函数梯度计算公式为
|
(15) |
|
(16) |
式中
由(15)式可以看出，在有限内存拟牛顿法反演的梯度计算过程中，需要求解灵敏度矩阵；但根据Kelbert等(2008)、Egbert和Kelbert(2012)的思路，可以采用伴随正演的方法直接计算灵敏度矩阵转置和数据向量的乘积，避免了显式计算灵敏度时所占用的大量时间和计算机资源.
3 棋盘模型试验理论模型采用棋盘渐变模型，对其数值模拟结果添加随机指数噪声得到稳定性较差的合成数据.本文的地磁测深三维正演模拟基于Uyeshima和Schultz(2000)提出的球坐标系中交错网格有限差分法，并假设外部场源可以P10的形式表示(Banks，1969)；观测点设计为120个(纬向：56°N—56°S，间隔16°；经向：0°~360°，间隔24°)；选择13个周期，等对数间隔的分布在3~113天之间；空气层厚度为2倍地球半径.正演的横向网格为5°×5°，在靠近南北极和赤道时进行更细的剖分，以增加数值计算精度；反演采用10°×10°的横向网格，在保证计算精度的同时也保证了反演效率.
3.1 模型描述三维地球模型在Kelbert等(2008)的一维背景模型(表 2)基础上进行层内扰动得到.根据地磁测深可探测的深度范围，理论模拟时，我们假设670~900 km内存在导电性异常，并用Kelbert等(2008)的棋盘模型模拟该异常.在第k层，该模型的电阻率分布计算公式为
|
|
表 2 一维背景模型，根据Kelbert等(2008) Table 2 The 1-D background model based on Kelbert et al. (2008) |
|
(17) |
式中ρk(θ, ϕ)表示电阻率；θ和ϕ分别为地磁坐标系中的纬度和经度；Slm为连带勒让德函数，Pl为零阶连带勒让德函数；l和m分别为球谐函数的阶数和次数，L为阶数的最大值；alm为球谐系数.为保证电阻率在背景模型基础上存在一个数量级范围内的变化，取a53值为1.6，其他alm均为零.670~900 km的电阻率异常分布如图 2.
|
图 2 用于反演测试的理论模型 棋盘模型位于670~900 km范围内，模型电导率由球谐系数a53=1.6的球谐函数描述；lg(σ)扰动幅值的范围为[-1，1]. Fig. 2 Theoretical model used for inversion test The checkerboard model between depth 670 km and 900 km and the corresponding coefficient a53 equals 1.6; the perturbed magnitude of lg(σ) is in the range [-1, 1]. |
为了更好的模拟真实C-响应的观测误差水平，依据Kelbert等(2008)和Fujii和Schultz(2002)的分析，数据误差δC(ω)可以采用
|
(18) |
计算，式中θ为地磁余纬度；β0和β1为经验系数，本文选择β0和β1的值分别为10和80.得到δC(ω)之后，计算遵循指数分布的随机噪声，并叠加到图 2模型的理论响应Cmodel(ω)上，以合成具有“飞点”的稳定性较差数据Cobs(ω)，其中C-响应的实部和虚部分别作为独立数据信息添加噪声.
指数概率密度分布函数p(d)的计算公式为
|
(19) |
式中，d表示数据集合，σ为数据标准差，〈d〉为数据平均数.为了更直观地展示非高斯噪声数据的特点，我们将3.3节中第六层存在棋盘模型时合成C-响应的误差信息按照顺序绘制在图 3中.可以发现，指数噪声较为复杂，变化范围大，存在多个“飞点”，并产生较多的异常大值，因此合成数据满足我们对稳定性较差数据的要求.
|
图 3 理论上生成的指数随机噪声，这些噪声将被按照顺序叠加到下文的理论C-响应上产生具有“飞点”的合成数据 横坐标为数据序号；纵坐标为生成的随机噪声与其标准差的比值，该比值越大说明噪声幅值越大. Fig. 3 The theoretical random noise data generated by exponential probability density function, and they will be added to the forward C-response in sequence to obtain the synthetic data with "outliers" The x axis is the series of data, the y axis is the ratio of generated noise and standard error. The noise level improves with the ratio increasing. |
由于我们选择的周期最大为113天，根据Kelbert等(2008)的研究表明，此时C-响应难以反映1600 km之下的电性结构.因此，我们仅反演该深度之上的电导率.反演的初始模型借用表 2的圈层结构，反演允许所有层导电性都可以改变，并允许电导率在410 km、520 km和670 km的界面上发生突变.图 4给出了部分深度的反演结果.由图可以看出，L1-范数的反演结果虽然对异常体的具体电阻率值反演不够精确，为异常幅值的60%左右；但棋盘模型样式得到了较好恢复，异常没有明显变形，且异常主要集中在真实存在异常的670~900 km层内.这说明即使在数据稳定性较差，存在“飞点”数据时，L1-范数反演也没有受到明显影响，仍具有良好的反演效果.
|
图 4 叠加指数噪声的合成数据L1-范数三维反演结果 图中白色圆圈为测点；我们采用反演得到的电导率与初始电导率比值的形式以突出异常分布. Fig. 4 Inversion result of 3-D L1-norm of the synthetic data with exponential noise The white circles are locations of observatories. In order to stand out the anomalies, we draw the distribution maps with conductivity ratio of inversion result and starting model. |
图 5给出了反演结果的响应数据与合成数据的交汇图，交汇点基本都分布在45°线附近，说明模型响应较好的拟合了合成数据.反演的拟合差RMS(Root Mean Square)最终值略大于1.0(图 6a)，这是L1-范数反演中没有过度拟合不稳定数据的结果；模型粗糙度随着迭代进行逐渐稳定，这说明迭代结束时模型参数已基本不发生变化(图 6b)，拟合已经较为充分.反演迭代过程中，终止条件采用拟合差、正则化因子和反演迭代次数等多个参数进行评价；如果第k次迭代与k-1次迭代的拟合差变化不超过0.002，就按照正则化因子递减的方法(吴小平和徐果明，1998)更新正则化因子，进行新的最优化求解，直到满足某终止条件.此次反演的初始正则化因子λ0=100，正则化因子对应的终止条件为λk < 10-4.
|
图 5 叠加指数噪声的合成数据与L1-范数三维反演模型响应的交汇图 (a)、(b)分别为C-响应实部和虚部的交汇图，图中虚线为45°直线，在完全拟合的情况下，所有点均应落在虚线上. Fig. 5 Cross plots of inversion response obtained by L1-norm and synthetic data with exponential noise (a) and (b) are the real and imaginary component of C-response, respectively. The dash line is the 45° diagonal line.In the case of full fitting, all points should fall on the dash line. |
|
图 6 反演数据拟合差(a)、模型粗糙度(b)、目标函数(c)和正则化因子(d)的迭代变化曲线 Fig. 6 Iterative curves of RMS (a), roughness of model (b), objective function (c) and regularization parameter (d) |
我们在1950年之后记录的全球地磁台站数据中，剔除总观测时长小于5年及数据观测质量极差的台站，得到了全球范围内129个位于中低纬度的地磁台站(表 1)，按照前文所述方法处理各台站数据得到C-响应用于三维反演.
三维反演采用10°×10°的横向网格，在表层12.65 km内按照1°×1°的网格计算以更好的消除海洋效应影响(Kelbert et al., 2009；李建平等，2018b)，海陆分布造成的电导率横向变化根据Chave和Thomson(2004)中的电导数据计算得到；反演初始模型采用表 2的背景模型，电性在12.65~1600 km内变化，1600 km之下的模型保持不变.反演采用有限内存拟牛顿法，对数据拟合误差的度量采用L1-范数，模型粗糙度度量采用L2-范数.计算数据拟合时所用的方差为数据处理时按照Jackknife方法得到的误差.反演中，正则化因子的初始值和变化方式与理论模型测试时的相同.由图 7各台站平均拟合差分布及图 8的部分台站数据拟合曲线可以看出，反演结果较好的拟合了观测数据.
|
图 7 全球C-响应数据L1-范数反演完成后各台站的数据拟合差分布 Fig. 7 The distribution of RMS in every observatory of the L1-norm inversion result of global C-responses |
|
图 8 全球地磁台站L1-范数反演部分台站的数据拟合对比曲线 各图中的三字母为台站代码；黑色实心圆为观测C-响应的实部，黑色菱形为观测C-响应的虚部；实线为反演结果C-响应的实部，虚线为反演结果C-响应的虚部. Fig. 8 Data fitting curves of global data′s L1-norm inversion in some stations The three letters in figures are the code of station and the black solid circles and rhombus are the real and imaginary component of observed C-response, while the solid and dash line are the real and imaginary component of C-response of inversion result. |
从图 9的反演模型对比可见，采用L1-范数计算数据拟合项的方法进行反演得到的结果和Kelbert等(2009)采用传统L2-范数对59个地磁台站的C-响应数据进行反演得到的结果基本一致，但也存在部分差别.比如，在地幔转换带深部，中国东北地区为高导异常，南欧和北非均为高阻；夏威夷在900 km以下均为高导；菲律宾海及以东地区在转换带表现为明显的高阻.而Kelbert等(2009)的反演结果中非洲东南在转换带地区存在的大范围高导异常在本文模型中的异常范围显著减小.
|
图 9 全球地幔转换带及下地幔上部电导率分布图 (a)采用本文L1-范数反演得到的结果；(b) Kelbert等(2009)采用传统L2-范数反演得到的结果.图中白色圆圈为各反演所采用的台站位置. Fig. 9 Global electrical conductivity model in mantle transition zone and upper part of lower mantle (a) is the result of L1-norm inversion and (b) is the result of L2-norm inversion of Kelbert et al. (2009). The white circles are observatories used in this paper and Kelbert et al. (2009), respectively. |
我们的反演结果也显示出新的导电异常，在北美地区520~900 km深度范围内呈现大面积高导，且主要集中在北美大陆的东侧，我们认为这与东太平洋板块俯冲进入美洲板块下方的地幔相关.由于采用的台站数量更多，本文反演结果还发现在南美洲南端，转换带表现为电导率的增高；在澳大利亚东南部，地幔转换带深部，存在一个明显的高导异常，该异常已经被Koyama等(2014)的研究发现.同时，Zhao(2004)由P-波层析成像得到的全球地幔波速结构也显示，在转换带深度范围内，南美洲南端和澳大利亚东南部均呈现低波速异常，与本文电性的高导异常有较好的对应关系，预示着转换带中可能发生了部分熔融.
5 结论本文实现了数据拟合误差采用L1-范数度量的地磁测深三维反演技术，并利用合成数据和全球实际观测C-响应数据进行了数值试验.通过反演结果及对比分析，我们得到如下结论：
(1) 当数据误差稳定性较差，遵从非高斯分布，比如指数分布时，对数据残差项采用L1-范数度量能有效降低“飞点”的影响，反演结果基本不受“飞点”影响.
(2) 实际数据三维反演得到的全球电性结构与Kelbert等(2009)的结果基本一致，反演结果均显示：在地幔转换带深部，中国东北地区为高导异常，南欧和北非均为高阻；夏威夷在900 km以下均为高导；菲律宾海及以东地区在转换带表现为明显的高阻.
(3) 本文结果显示，在520~900 km深度内，北美地区东侧存在大范围的高导异常体，可能与东太平洋板块俯冲进入美洲板块有关.
(4) 由于采用的台站数量更多，我们的反演结果还显示：在南美洲南端，转换带表现为明显的高导；在澳大利亚东南部，地幔转换带深部，也存在一个明显的高导异常.这些新发现的高导异常特征，与已发表的全球地震P-波层析成像结果中的低速异常区对应，预示着这些地区可能发生了部分熔融.
(5) L1-范数反演能够充分利用全球C-响应数据信息，提高地磁测深对地球深部电性结构的分辨能力，更好的研究全球地幔导电结构.
(6) 本文提出的基于L1-范数度量数据拟合项的三维反演思路，并不局限于地磁测深，可用于其他类似问题的正则化反演.
致谢 感谢国际地磁数据中心(World Data Centre for Geomagnetism)提供的全球地磁台站观测数据下载服务(http://www.wdc.bgs.ac.uk/data.html)，感谢Chave A D提供的BIRRP数据处理程序，感谢GMT中文社区相关资源对本文绘图的重要帮助(https://gmt-china.org/).感谢两位审稿专家给出的建设性建议，使本文更加严谨、完善.
Avdeev D, Avdeeva A. 2009. 3D magnetotelluric inversion using a limited-memory quasi-newton optimization. Geophysics, 74(3): F45-F57. DOI:10.1190/1.3114023 |
Banks R J. 1969. Geomagnetic variations and the electrical conductivity of the upper mantle. Geophysical Journal International, 17(5): 457-487. DOI:10.1111/j.1365-246X.1969.tb00252.x |
Booker J R, Favetto A, Pomposiello M C. 2004. Low electrical resistivity associated with plunging of the Nazca flat slab beneath Argentina. Nature, 429(6990): 399-403. DOI:10.1038/nature02565 |
Chave A D, Thomson D J. 2004. Bounded influence magnetotelluric response function estimation. Geophysical Journal International, 157(3): 988-1006. DOI:10.1111/gji.2004.157.issue-3 |
Egbert G D, Kelbert A. 2012. Computational recipes for electromagnetic inverse problems. Geophysical Journal International, 189(1): 251-267. DOI:10.1111/gji.2012.189.issue-1 |
Ekblom H. 1988. A new algorithm for the Huber estimator in linear models. BIT Numerical Mathematics, 28(1): 123-132. DOI:10.1007/BF01934700 |
Farquharson C G, Oldenburg D W. 1998. Non-linear inversion using general measures of data misfit and model structure. Geophysical Journal International, 134(1): 213-227. DOI:10.1046/j.1365-246x.1998.00555.x |
Farquharson C G. 2008. Constructing piecewise-constant models in multidimensional minimum-structure inversions. Geophysics, 73(1): K1-K9. |
Fujii I, Schultz A. 2002. The 3D electromagnetic response of the earth to ring current and auroral oval excitation. Geophysical Journal International, 151(3): 689-709. DOI:10.1046/j.1365-246X.2002.01775.x |
Fukao Y, Koyama T, Obayashi M, et al. 2004. Trans-Pacific temperature field in the mantle transition region derived from seismic and electromagnetic tomography. Earth & Planetary Science Letters, 217(3-4): 425-434. |
Huang X G, Xu Y S, Karato S. 2005. Water content in the transition zone from electrical conductivity of wadsleyite and ringwoodite. Nature, 434(7034): 746-749. DOI:10.1038/nature03426 |
Karato S. 1990. The role of hydrogen in the electrical conductivity of the upper mantle. Nature, 347(6290): 272-273. DOI:10.1038/347272a0 |
Kelbert A, Egbert G D, Schultz A. 2008. Non-linear conjugate gradient inversion for global EM induction: Resolution studies. Geophysical Journal International, 173(2): 365-381. DOI:10.1111/j.1365-246X.2008.03717.x |
Kelbert A, Schultz A, Egbert G. 2009. Global electromagnetic induction constraints on transition-zone water content variations. Nature, 460(7258): 1003-1006. DOI:10.1038/nature08257 |
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. Journal of Geophysical Research: Solid Earth, 111(B10): 207-208. DOI:10.1029/2006JB004270 |
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. Journal of Geophysical Research, 116(B1): B01103. DOI:10.1029/2010JB007458 |
Koyama T, Khan A, Kuvshinov A. 2014. Three-dimensional electrical conductivity structure beneath Australia from inversion of geomagnetic observatory data: evidence for lateral variations in transition-zone temperature, water content and melt. Geophysical Journal International, 196(3): 1330-1350. DOI:10.1093/gji/ggt455 |
Li J P, Weng A H, Li S W, et al. 2018a. 3-D forward method for geomagnetic depth sounding based on finite difference method in spherical coordinate. Journal of Jilin University (Earth Science Edition) (in Chinese), 48(2): 411-419. |
Li J P, Weng A H, Li S W, et al. 2018b. The influence of ocean effect on geomagnetic observations in coastal areas of China: A case study of the Guangzhou observatory. Chinese Journal of Geophysics (in Chinese), 61(2): 649-658. DOI:10.6038/cjg2018L0433 |
Liu Y H, Yin C C. 2013. 3D inversion for frequency-domain HEM data. Chinese Journal of Geophysics (in Chinese), 56(12): 4278-4287. DOI:10.6038/cjg20131230 |
Manthilake G, Bolfan-Casanova N, Novella D, et al. 2016. Dehydration of chlorite explains anomalously high electrical conductivity in the mantle wedges. Science Advances, 2(5): e1501631. DOI:10.1126/sciadv.1501631 |
Menke W. 1984. Geophysical Data Analysis: Discrete Inverse Theory. New York: Academic Press.
|
Munch F D, Grayver A V, Kuvshinov A, et al. 2018. Stochastic inversion of geomagnetic observatory data including rigorous treatment of the ocean induction effect with implications for transition zone water content and thermal structure. Journal of Geophysical Research: Solid Earth, 123(1): 31-51. DOI:10.1002/2017JB014691 |
Nocedal J, Wright S J. 1999. Numerical Optimization. New York: Springer.
|
Oldenburg D W. 1984. An introduction to linear inverse theory. IEEE Transactions on Geoscience & Remote Sensing, GE-22(6): 665-674. |
Olsen N. 1998. The electrical conductivity of the mantle beneath Europe derived from C-responses from 3 to 720 hr. Geophysical Journal International, 133(2): 298-308. DOI:10.1046/j.1365-246X.1998.00503.x |
Peyronneau J, Poirier J P. 1989. Electrical conductivity of the Earth′s lower mantle. Nature, 342(6249): 537-539. DOI:10.1038/342537a0 |
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. Geophysical Journal International, 203(3): 1864-1872. DOI:10.1093/gji/ggv407 |
Schmucker U. 1987. Substitute conductors for electromagnetic response estimates. Pure & Applied Geophysics, 125(2-3): 341-367. |
Semenov A, Kuvshinov A. 2012. Global 3-D imaging of mantle conductivity based on inversion of observatory C-responses-Ⅱ: Data analysis and results. Geophysical Journal International, 191(3): 965-992. DOI:10.1111/j.1365-246X.2012.05665.x |
Semenov V Y, Jozwiak W. 2006. Lateral variations of the mid-mantle conductance beneath Europe. Tectonophysics, 416(1-4): 279-288. DOI:10.1016/j.tecto.2005.11.017 |
Sun J, Kelbert A, Egbert G D. 2015. Ionospheric current source modeling and global geomagnetic induction using ground geomagnetic observatory data. Journal of Geophysical Research: Solid Earth, 120(10): 6771-6796. DOI:10.1002/2015JB012063 |
Tang J, Zhao G Z, Wang J J, et al. 2006. Study of the formation mechanism for volcanism in Northeast China based on deep electric structure. Acta Petrologica Sinica (in Chinese), 22(6): 1503-1510. |
Uyeshima M, Schultz A. 2000. Geoelectromagnetic induction in a heterogeneous sphere: a new three-dimensional forward solver using a conservative staggered-grid finite difference method. Geophysical Journal International, 140(3): 636-650. DOI:10.1046/j.1365-246X.2000.00051.x |
Wang Q, Huang Q H. 2016. The spatio-temporal characteristics of geomagnetic induction vectors in North China. Chinese Journal of Geophysics (in Chinese), 59(1): 215-228. DOI:10.6038/cjg20160118 |
Weng A H, Li D J, Li Y B, et al. 2015. Selection of parameter types in controlled source electromagnetic method. Chinese Journal of Geophysics (in Chinese), 58(2): 697-708. DOI:10.6038/cjg20150230 |
Wu X P, Xu G M. 1998. Improvement of Occam′s inversion for MT data. Chinese Journal of Geophysics (Acta Geophysica Sinica) (in Chinese), 41(4): 547-554. |
Xu G J, Tang J, Huang Q H, et al. 2015. Study on the conductivity structure of the upper mantle and transition zone beneath North China. Chinese Journal of Geophysics (in Chinese), 58(2): 566-575. DOI:10.6038/cjg20150219 |
Xu Y S, Shankland T J, Poe B T. 2000. Laboratory-based electrical conductivity in the Earth′s mantle. Journal of Geophysical Research: Solid Earth, 105(B12): 27865-27875. DOI:10.1029/2000JB900299 |
Yoshino T, Manthilake G, Matsuzaki T, et al. 2008. Dry mantle transition zone inferred from the conductivity of wadsleyite and ringwoodite. Nature, 451(7176): 326-329. DOI:10.1038/nature06427 |
Zhang Y H, Weng A H, Li S W, et al. 2019. Estimation of C-responses of geomagnetic depth sounding based on global smooth constraint. Chinese Journal of Geophysics (in Chinese), 62(5): 1898-1907. DOI:10.6038/cjg2019M0199 |
Zhao D P. 2004. Global tomographic images of mantle plumes and subducting slabs: insight into deep Earth dynamics. Physics of the Earth & Planetary Interiors, 146(1-2): 3-34. |
Zhao G Z, Tang J, Liang J G, et al. 2001. Measurement of network-MT in two areas of NE China for study of upper mantle conductivity structure of the back-arc region. Seismology and Geology (in Chinese), 23(2): 143-152. |
李建平, 翁爱华, 李世文, 等. 2018a. 基于球坐标系下有限差分的地磁测深三维正演. 吉林大学学报(地球科学版), 48(2): 411-419. |
李建平, 翁爱华, 李世文, 等. 2018b. 海洋效应对中国沿海地磁观测影响——以广州台站为例. 地球物理学报, 61(2): 649-658. DOI:10.6038/cjg2018L0433 |
刘云鹤, 殷长春. 2013. 三维频率域航空电磁反演研究. 地球物理学报, 56(12): 4278-4287. DOI:10.6038/cjg20131230 |
汤吉, 赵国泽, 王继军, 等. 2006. 基于地下电性结构探讨中国东北活动火山形成机制. 岩石学报, 22(6): 1503-1510. |
王桥, 黄清华. 2016. 华北地磁感应矢量时空特征分析. 地球物理学报, 59(1): 215-228. DOI:10.6038/cjg20160118 |
翁爱华, 李大俊, 李亚彬, 等. 2015. 数据类型对三维地面可控源电磁勘探效果的影响. 地球物理学报, 58(2): 697-708. DOI:10.6038/cjg20150230 |
吴小平, 徐果明. 1998. 大地电磁数据的Occam反演改进. 地球物理学报, 41(4): 547-554. DOI:10.3321/j.issn:0001-5733.1998.04.013 |
徐光晶, 汤吉, 黄清华, 等. 2015. 华北地区上地幔及过渡带电性结构研究. 地球物理学报, 58(2): 566-575. DOI:10.6038/cjg20150219 |
张艳辉, 翁爱华, 李世文, 等. 2019. 基于全局光滑约束的地磁测深C-响应估计. 地球物理学报, 62(5): 1898-1907. DOI:10.6038/cjg2019M0199 |
赵国泽, 汤吉, 梁竞阁, 等. 2001. 用大地电磁网法在长春等地探测上地幔电导率结构. 地震地质, 23(2): 143-152. DOI:10.3969/j.issn.0253-4967.2001.02.003 |