2019, Vol. 62
地磁台站观测的磁场变化数据包含许多上地幔深部和下地幔上部电性结构的信息(Armadillo et al., 2001).利用这些地球内部的电性结构信息可以探索地球的演化历史、动力学特征等问题(Fukao et al., 2004；Booker et al., 2004；Khan et al., 2006, 2011；Püthe et al., 2015).
Banks(1969)首先提出将地磁测深数据转换为C-响应函数，并通过反演C-响应函数来研究地球内部的电性结构(Banks, 1969；Schultz and Larsen, 1987).目前被广泛应用来求取C-响应函数的方法有两种，一种是Z/H方法，该方法在假设磁层源满足P10的条件下，可以由单个台站地磁数据求取C-响应函数(Banks，1969；Schultz and Larsen, 1987)，并被应用于研究美国西南地区的导电结构(Schmucker，1970).另一种方法是Z/Y方法，是由Schmucker(1979, 1985)进一步提出的.在这种方法中，垂直分量Z由本台站提供，而水平分量Y由全球台站数据的球谐系数表示(Olsen，1999).随后，Schmucker(1990)和Olsen(1992)使用校正因子对Z/Y方法进行了改进，这使得C-响应的求取更加合理，并拓宽了C-响应的周期(Olsen, 1992, 1998).
C-响应求取通常采用最小二乘法.这里，我们借助Chave和Thomson(2004)的BIRRP(Bounded Influence，Remote Reference Processing)数据处理软件，结合Semenov和Kuvshinov(2012)提出的台站数据自参考方法，对地磁台站观测数据进行C-响应估计.由于BIRRP技术的核心是最小二乘法，因此下文中将其称为LS方法.
在使用LS方法估计C-响应时，当台站数据质量较好，观测时间较长，比如FUR和HER台站，求取的C-响应很稳定且与之前的研究结果类似(Olsen，1998；Semenov and Kuvshinov, 2012).然而，当台站数据质量较差时，得到的响应曲线存在上下震荡和“飞点”情况.Semenov和Kuvshinov(2012)认为这与仪器和环境噪声以及真实信号的相对微弱有关.
分析C-响应的求取过程，我们发现除了上述直接原因外，求取C-响应的方法本身也存在缺陷：目前的方法都是逐个频率求取的C-响应，没有考虑不同频率之间C-响应的连续性.事实上，和大地电磁测深(Magnetotelluric)相似，地磁测深的频率范围更低，因此可以忽略位移电流，电场和磁场在地球内部是以扩散的方式传播的(Chave and Jones, 2012).C-响应是地下导电信息的综合反映，应该表现为一条光滑的曲线.很明显，上述的C-响应求取方法没有将该因素考虑进去，导致求取的响应曲线会震荡、“飞点”和不连续，而这被认为是不可靠的(Püthe et al., 2015).
Püthe和Kuvshinov(2014)也注意到该问题.为得到更好的Q矩阵响应曲线采用了多频估计的方法.该方法在求解过程中强加了约束条件，它能有效减少噪声影响，得到更优的传递函数.
为了克服C-响应估计中的噪声干扰问题，并且由于C-响应与Q-响应实质上的一致性，本文将约束技术引进到C-响应估计中.首先，我们获取了纵向和横向磁场分量的频谱，然后基于正则化反演方法(RI)并通过加入光滑约束条件构建了目标函数，极小化目标函数可以求解约束后的C-响应，但其随着正则化因子的变化而变化.最优的C-响应曲线通过L-曲线和V-曲线来确定.最后，为了测试该方法的可靠性，分别使用合成数据和实际数据来进行测试，反演采用有限内存拟牛顿法(L-BFGS方法)，其反演精度和可靠性已被李世文等(2017)进行了分析验证.结果表明该方法能够得到更理想的响应曲线，增加反演结果的可靠性.
1 地磁测深C-响应在磁层源可以近似表示为P10的条件下，C-响应的计算公式可以表示为(Banks，1969)
|
(1) |
式中a0为地球的平均半径；Hr、Hθ分别为地表处垂直指向地心和水平北向的磁场；θ为地磁余纬度；ω为角频率.关于Hr和Hθ的求取参考Semenov和Kuvshinov(2012).
为了更好显示C-响应的特点，本文计算了Püthe等(2015)全球平均一维导电模型的C-响应曲线(见图 1).正演理论基于Banks(1969)和李世文等(2017)的文章.由图 1b可见，C-响应曲线的实部和虚部都相当平滑，这也反映了地磁测深的体积效应.因此，在实际地磁测深数据处理中，C-响应曲线在排除噪声干扰之后应该是光滑的，这也是本文计算C-响应的出发点.
|
图 1 (a) 全球平均一维导电模型(Püthe et al., 2015)及其C-响应曲线(b) Fig. 1 (a) Global 1D average conductivity model (Püthe et al., 2015); (b) Corresponding C-response curve |
在公式(1)中，将ωi记为i，Hr(ωi)写为Vi，
|
(2) |
其中，N为频率个数.通常C-响应估计时是由(2)式逐个频率计算.具体为，对于一个给定的频率ωi，考虑地磁信号的平稳性，取分布在其附近时窗内的K个频率作为频率ωi的地磁场分量.因此，对于单个频点，C-响应估计能够表示为如下的方程(Banks，1969)：
|
(3) |
逐频率通过(3)式可以计算得到所有频率的响应，从而获得整个C-响应曲线.虽然LS方法的应用是为了减少噪声的干扰来获得更稳定的C-响应.然而，当数据存在较大干扰时，得到的曲线受干扰影响仍会出现震荡或“飞点”.将北京台站(40.04°N, 116.175°E)2001年到2015年的地磁场观测数据进行处理，得到其C-响应曲线，见图 2.由图可见，C-响应曲线存在上下震荡和“飞点”，尤其是长周期时.除了难以解决的噪声和时间序列较短的问题，产生上述问题的另一个可能原因是，LS方法逐个频率求取响应值，割裂了相邻频率C-响应间的联系.
|
图 2 北京台站使用LS方法计算的C-响应 Fig. 2 C-responses of Beijing (BJI) observatory calculated by LS method |
因此，我们提出了新的求解C-响应的方法.该方法不仅仅利用(3)式来抑制噪声的影响，同时考虑了C-响应的光滑特性，要对其进行光滑约束.基于上述条件，将C-响应的求取转化为一个包含光滑约束的最优化问题，其中光滑项为：
|
(4) |
满足
|
(5) |
其中W是模型约束矩阵，也称为光滑矩阵.(4)式和(5)式可以转化为正则化反演问题，其目标函数为
|
(6) |
其中，λ是平衡数据误差和模型光滑度的正则化参数(Tikhonov，1963；王家映，2002)；数据预测误差Φd被定义为
|
(7) |
当(6)式取极小时
|
(8) |
从(8)式可以得到进行光滑约束后的C-响应，我们称之为RI方法.该方法与以往方法的本质区别是将整个频率范围作为一个整体求解，并加入了光滑约束.该方法将能够更好的压制噪声的干扰，得到更真实的C-响应.
3 影响因素从(8)式中可以看出，有两个因素会对C-响应计算产生影响：一是光滑矩阵，另一个是正则化参数.
3.1 光滑矩阵光滑矩阵相当于地球物理反演过程中的稳定因子，它是对模型解的空间进行限制，以减少多解性，求得稳定解(Constable et al., 1987；Zhdanov，1993).在本文中，我们使用光滑矩阵来对C-响应进行约束，使其接近于真实值而不至于出现解的偏离.对于公式中出现的W，本文选取了最大平滑稳定因子，又分为一阶导数(W1)和二阶导数(W2)(Constable et al., 1987).
一阶导数：
|
二阶导数：
|
其中，N为频率个数.
3.2 正则化参数正则化因子作为数据目标函数和模型约束目标函数的平衡系数，决定着数据处理结果的优劣(Tikhonov and Arsenin, 1977).正则化因子过小，则得到的结果主要拟合观测数据；若正则化因子过大，则会导致结果过于光滑从而造成数据的真实性降低.因此通常使用L-曲线法准则来选取最优正则化因子(Hansen，1992；Pardo and Torres-Verdín，2015).
Frasso和Eilers(2015)通过一系列的理论及实际数据的测试发现，在寻找最优光滑效果时L-曲线的形状与数据的噪声有着相当大的关系，当数据噪声水平较大时，L-曲线的拐点随着噪声水平的增加越发不明显(如图 3所示).由于L-曲线的拐点对应数据点最密的位置(Frasso and Eilers, 2015)，因此通过求取L-曲线上相邻两点之间距离得到的V-曲线被引入进来(如图 2b)，它可以很容易地定位数据分布最密的位置，从而找到最优的正则化因子.
|
图 3 不同噪声水平的L-曲线和V-曲线(Frasso and Eilers, 2015) (a) L-曲线；(b)对应的V-曲线.图中对应的最优正则化因子用黑点表示. Fig. 3 L-curves and V-curves at different noise levels (modified from Frasso and Eilers, 2015) (a) L-curves; (b) Corresponding V-curves. The optimal regularization parameters are denoted by black dots in the curves. |
为了验证L-曲线及V-曲线选取最优正则化参数的准确性和RI方法求取C-响应的正确性，本文以图 1的响应曲线为例进行讨论.理论模拟时假设H为1，依据公式(2)V=CH，那么V=C，其中C由理论模拟得到.假设噪声只干扰V，使得V′=V+e，因此e可以直接被加在理论C-响应上.将理论模型的响应曲线加入不同水平的高斯噪声形成合成数据，然后使用RI方法进行处理，将其结果与原始理论响应进行对比.
为了得到最优的正则化参数，其初始值设置为106，且每次计算变成0.8λ，设置计算次数为200.图 4a为合成数据一阶导数的L-曲线，σ为噪声幅度，与图 3相似的是，当噪声水平增加时，L-曲线的拐角变得不明显(如图 4a所示)，很难用其来确定最优正则化参数.因此，本文通过计算L-曲线中相邻数据点之间的距离从而得到了对应的V-曲线，如图 4b所示.由图可见，即使在噪声水平较高时，V-曲线也能轻易定位数据点最密的位置，从而找到最优正则化参数.图 4d和图 4e分别是不同噪声水平二阶导数的L-曲线和V-曲线，可以得到相似的结论.
|
图 4 合成数据处理结果 (a)一阶导数光滑矩阵不同噪声水平的L-曲线和(b)V-曲线；(d)二阶导数光滑矩阵不同噪声水平的L-曲线和(e)V-曲线；(c)合成数据一阶导数光滑矩阵得到的C-响应曲线；(f)合成数据二阶导数光滑矩阵得到的C-响应曲线. Fig. 4 Results of the synthetic C-responses (a) L-curves and (b) V-curves with different noise levels for the first-order derivative smoothness matrix; (d) L-curves and (e) V-curves with different noise levels for the second-order derivative smoothness matrix; (c) C-responses curves of the synthetic data for the first-order derivative smoothness matrix and (f) that of the second-order derivative smoothness matrix. |
采用上述方法估计出最优正则化参数后，将8%高斯噪声的合成数据分别选用一阶导数和二阶导数光滑矩阵进行C-响应重新估计，结果见图 4c和图 4f.由图可见，不同光滑矩阵处理得到的C-响应都能使受到噪声污染的响应曲线得到很好的恢复，这也说明了L-曲线和V-曲线选取最优正则化参数的正确性.
为了进一步测试影响C-响应估计的因素，本文对正演模型进行了重新设计，如表 1所示.将其正演响应加入不同水平的高斯噪声后进行了C-响应的重新求取，光滑矩阵选用二阶导数，求取结果如图 5所示.由图 5可见，当噪声水平相对较小时，RI方法得到的响应曲线与真实曲线更加接近；当噪声水平较大时，求取的响应曲线与理论曲线大致趋势相同，但存在一定程度的偏差(见图 5c实部)；同时当合成数据连续偏离理论曲线的一侧时，会造成求取结果向合成数据弯曲(如图 5a实部中部及图 5b虚部的末端).总体上，该方法在数据存在10%以内的高斯噪声时，能够取得较为理想的效果.
|
图 5 不同噪声水平的合成数据处理结果 (a) 5%高斯噪声合成数据；(b) 10%高斯噪声合成数据；(c) 15%高斯噪声合成数据. Fig. 5 Results of synthetic data processing at different noise levels (a) Synthetic data of 5% Gauss noise; (b) Synthetic data of 10% Gauss noise; (c) Synthetic data of 15% Gauss noise. |
|
表 1 用于模型测试的地电模型 Table 1 Geoelectric model for model test |
为了验证RI方法在实际应用中求取C-响应的可靠性，本文以兰州台站(LZH，36.087°N，103.845°E)为例进行分析.除部分数据缺失外，兰州台站已有超过20年的数据记录(1980—2011)，数据质量较好(Khan et al., 2011).为了获得不好的响应估计，本文截取2001—2009年的数据片段进行处理，图 6为时间序列，可见数据存在一些突变和“飞点”.
|
图 6 兰州台站磁场分量的时均值序列 (a)地理北分量；(b)东分量；(c)垂直分量. Fig. 6 Hourly mean time-series of the magnetic field components at LZH observatory (a) North component; (b) East component; (c) Vertical component. |
为得到最终的C-响应，首先需要估计出RI方法的最优正则化参数.图 7a和图 7b分别为一阶导数光滑矩阵处理得到的L-曲线和V-曲线.图 7c和图 7d分别为二阶导数光滑矩阵处理得到的L-曲线和V-曲线.由于数据受到噪声污染，L-曲线的拐点不太明显，因此，本文建议通过V-曲线的最小值来确定最优的正则化因子.
|
图 7 兰州台站处理结果的L-曲线和V-曲线 (a)一阶导数光滑矩阵的L-曲线和(b)V-曲线；(c)二阶导数光滑矩阵的L-曲线和(d)V-曲线. Fig. 7 L-curve and V-curve calculated using data from Lanzhou geomagnetic observatory (a) and (b) L-curve and V-curve of the first-order derivative smoothing factor; (c) and (d) L-curve and V-curve of the second-order derivative smoothing factor. |
图 8为对应图 7中最优正则化参数时不同光滑矩阵计算得到的C-响应曲线.由图可见，小周期时数据质量较高，因此通过LS和RI方法计算得到的响应曲线较为相似；然而在大周期时，两种方法计算得到的曲线大致趋势相同，但LS存在着明显的突变和“飞点”；从LS方法计算的数据误差棒也可以看出，响应曲线的可信度并不高(尤其在大周期时)；而通过RI计算得到的响应曲线能产生一个光滑的C-响应曲线.
|
图 8 兰州台站C-响应曲线的处理结果 (a) RI方法得到的平方一致性; (b) RI和LS方法处理得到的C-响应. Fig. 8 C-response curves based on data from Lanzhou observatory (a) Square consistence of RI method; (b) C-responses processed by RI and LS methods. |
为了更好的评估数据质量，本文使用平方一致性coh2进行质量衡量，其表达式为：
|
(9) |
其中，H*是H矩阵的共轭.平方一致性可以用来作数据质量评价指标，低的相干系数表示Z或H分量上的相对误差较大.因此，它能反映计算C-响应方法的抗噪能力(Semenov and Kuvshinov, 2012).
图 8上方为该台站数据的平方一致性.由图可见，一致性曲线比较震荡，且大周期时更小，说明在本文所用时间序列下，该台站的相对噪声较大.
另外，由图 8可见，RI方法计算的结果更加平滑，没有出现明显的拟合突变点的情况.同时，还可以看到，一阶光滑矩阵(黑色实线)和二阶光滑矩阵(黑色虚线)都能得到较光滑的C-响应，且结果较为相似；但前者在曲线的端点处未能拟合LS方法的趋势，曲线趋于平缓，而后者能更好的与LS方法保持一致.此外，由于二阶光滑矩阵的L-曲线的拐点更加明显，因此在下文中，统一采用二阶光滑矩阵来求取C-响应.
本文还对满洲里(MZL)台站和北京(BJI)台站进行了重新处理.MZL台站受噪声干扰较小，数据质量相对较好；而BJI台站受轨道交通、人文环境等噪声干扰严重，数据质量较差.计算结果如图 9所示.由C-响应可见，在数据质量较好的MZL台站，LS方法和RI方法计算的结果保持着很高的一致性，但后者更加平滑，这说明了RI方法的准确性.在数据质量较差的BJI台站，LS方法处理得到的响应曲线在大周期时存在明显的跳点，而采用RI方法得到的响应曲线更加平稳、光滑，极大减少了数据噪声对C-响应曲线的影响，这证明了该方法在求取C-响应时的优越性.
|
图 9 LS和RI方法实际台站处理结果 (a)满洲里台站；(b)北京台站. Fig. 9 Processing results of real data from two observatories using LS and RI methods (a) Manzhouli observatory; (b) Beijing observatory. |
为了进一步验证RI方法求取C-响应的合理性，本文从一维反演的角度，根据反演拟合情况和反演结果可靠性两个方面进行说明.反演采用一维L-BFGS方法(李世文等，2017).
实际数据仍以兰州台为例，采用图 8中的C-响应进行反演，结果见图 10.同时为了评价反演结果的可靠性，我们计算了反演模型的置信区间(Nabighian，1987)，具体的计算步骤见李世文等(2017).
|
图 10 兰州地磁台站C-响应反演结果 (a) LS和RI方法计算的C-响应及其反演结果；(b)反演的电导率-深度剖面和其对应的95%置信区间；(c) LS和RI方法得到C-响应的反演的拟合差对比. Fig. 10 Inversion results of data from C-responses at LZH observatory (a) C-responses estimated by LS and RI methods and their inversion results; (b) Conductivity-depth profiles from inversion and the corresponding 95% confidence intervals; (c) Inversion misfits comparison of the C-responses by LS and RI methods. |
由于噪声干扰，在图 10a中LS方法计算的C-响应曲线存在一些震荡和“飞点”，特别在大周期时更明显，因此反演结果只能拟合其变化趋势，尤其是虚部；而RI方法计算的响应曲线，其反演结果能够很好的拟合.图 10b为不同方法的反演结果及其对应的置信区间，其中红线为在更长时间序列的前提下，Semenov和Kuvshinov(2012)对兰州台站的反演结果.
由图可见，经过RI方法处理的反演结果置信区间更小，说明其反演可靠性更高.RI方法处理得到的结果与Semenov和Kuvshinov(2012)的结果更接近，这说明即使在数据长度较短或者质量较差的情况下，RI方法依然能够较好恢复实际响应曲线，从而得到更加可靠的反演结果.
虽然在深度范围400~900 km时，三者的反演结果很接近.然而在深度更大时，LS方法与Semenov和Kuvshinov(2012)的结果产生较大的偏差，这可能是由于LS计算得到的响应曲线不够光滑，大周期时，存在上下震荡和“飞点”，这影响了反演的拟合，导致反演结果较差.
此外，图 10c比较了反演的拟合差变化曲线，该拟合差为反演模型的C-响应与实测C-响应的差.由于RI方法计算得到的C-响应更加光滑，因此该反演拟合差下降更快，且最终拟合差更低.
6 结论在本文中，基于解的光滑约束和正则化反演理论，我们提出了计算地磁测深C-响应的正则化优化方法.通过合成数据和实际地磁台站数据的测试，说明了该方法在计算C-响应曲线上的优越性；反演结果的对比，则进一步体现了该方法的优势.
从理论响应和实测数据的测试结果中，可以得到以下结论：
(1) 使用正则化优化方法计算地磁测深C-响应曲线是可行的.该方法不再是单个频点计算，而是将所有频率响应作为一个整体且加上了光滑约束.通过该方法能够得到更真实的C-响应曲线，有效抑制由噪声干扰或其他因素引起的曲线偏差.
(2) 正则化参数在正则化优化方法中起着重要的作用.最优正则化因子可以通过V-曲线方法来进行确定.
(3) 该方法可以将数据质量不好的台站加以处理从而得到相对较好的响应数据，使得由于数据质量不好而无法参与反演的数据重新参与反演并得到稳定的反演效果，有效提高稀疏的地磁台站数据的利用率，从而有效提高反演结果的分辨率.
致谢 感谢Chave A D提供的BIRRP软件包，感谢国家台网中心提供的地磁数据，同时感谢审稿人对本文提出的宝贵修改意见，使本文的研究更加充实、完整.
Armadillo E, Bozzo E, Cerv V, et al. 2001. Geomagnetic depth sounding in the northern Apennines (Italy). Earth, Planets & Space, 53(5): 385-396. |
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. |
Chave A D, Jones A G. 2012. The Magnetotelluric Method: Theory and Practice. Cambridge University Press.
|
Constable S C, Parker R L, Constable C G. 1987. Occam′s inversion: a practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics, 52(3): 289-300. |
Frasso G, Eilers P. 2015. L- and V-curves for optimal smoothing. Statistical Modelling, 15(1): 91-111. DOI:10.1177/1471082X14549288 |
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. |
Hansen P C. 1992. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Review, 34(4): 561-580. DOI:10.1137/1034115 |
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 |
Li S W, Weng A H, Li J P, et al. 2017. 1-D inversion of C-response data from geomagnetic depth sounding with shallow resistivity constraint. Chinese Journal of Geophysics (in Chinese) (in Chinese), 60(3): 1201-1210. DOI:10.6038/cjg20170330 |
Nabighian M N. 1987. Electromagnetic Methods in Applied Geophysics: Volume 1, Theory. Tulsa: Society of Exploration Geophysics, doi: 10.1190/1.9781560802631.
|
Olsen N. 1992. Day-to-Day C-Response estimation for Sq from 1 cpd to 6 cpd using the Z:Y-method. Journal of Geomagnetism & Geoelectricity, 44(6): 433-447. |
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. |
Olsen N. 1999. Long-period (30 days-1 year) electromagnetic sounding and the electrical conductivity of the lower mantle beneath Europe. Geophysical Journal International, 138(1): 179-187. |
Pardo D, Torres-Verdín C. 2015. Fast 1D inversion of logging-while-drilling resistivity measurements for improved estimation of formation resistivity in high-angle and horizontal wells. Geophysics, 80(2): E111-E124. |
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. Geophysical Journal International, 197(2): 768-784. |
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. 1970. Anomalies of geomagnetic variations in the southwest United States. Bull. Scripps Inst. Ocean., Univ. Calif., 13: 1-165. |
Schmucker U. 1979. Erdmagnetische Variationen und die elektrische Leitfähigkeit in tieferen Schichten der Erde.// Sitzungsbericht und Mitteilungen Braunschweigische Wiss. Gesellschaft, Sonderheft, 4: 45-102.
|
Schmucker U. 1985. Magnetic and electric fields due to electromagnetic induction by external sources, electrical properties of the earth′s interior.// Landolt-Börnstein, New-Series, 5/2b. Berlin: Springer-Verlag.
|
Schmucker U. 1990. Die eindringtiefen tagesperiodischer variationen. // Haak V, Homilius J eds. Protokoll Koll. Elektromagnetische Tiefenforschung, Hornburg. Niedersächsisches Landesamt, Hannover, 31-66.
|
Schultz A, Larsen J C. 1987. On the electrical conductivity of the mid-mantle: ⅠCalculation of equivalent scalar MT response functions. Geophysical Journal International, 88(3): 733-761. |
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. |
Tikhonov A N. 1963. Regularization of incorrectly posed problems. Soviet Mathematics Doklady, 4(1): 1624-1627. |
Tikhonov A N, Arsenin V Y. 1977. Solutions of Ill-Posed Problems. New York: John Wiley & Sons.
|
Wang J Y. 2002. Inverse Theory in Geophysics (in Chinese). Beijing: Higher Education Press.
|
Zhdanov M S. 1993. Tutorial: Regularization in Inversion Theory. Colorado: Colorado School of Mines.
|
李世文, 翁爱华, 李建平, 等. 2017. 浅部约束的地磁测深C-响应一维反演. 地球物理学报, 60(3): 1201-1210. DOI:10.6038/cjg20170330 |
王家映. 2002. 地球物理反演理论. 北京: 高等教育出版社.
|