基于三维大地电磁AR-QN反演的长白山天池火山区电性结构

引用本文

阮帅, 汤吉, 董泽义, 王立凤, 邓琰, 韩冰. 基于三维大地电磁AR-QN反演的长白山天池火山区电性结构[J]. 地震地质, 2020,42(6): 1282-1300.
RUAN Shuai, TANG Ji, DONG Ze-yi, WANG Li-feng, DENG Yan, HAN Bing. ELECTRIC STRUCTRUE MODEL OF TIANCHI VOLCANO IN CHANGBAI MOUNTAINS BASED ON THREE-DIMENSIONAL AR-QN MAGNETOTELLURIC INVERSION[J]. SEISMOLOGY AND GEOLOGY, 2020,42(6): 1282-1300. 复制到剪切板

DOI:10.3969/j.issn.0253-4967.2020.06.002
Permissions

《地震地质》编辑部所有

基于三维大地电磁AR-QN反演的长白山天池火山区电性结构

阮帅¹⁾, 汤吉¹⁾, 董泽义¹⁾, 王立凤¹⁾, 邓琰²⁾, 韩冰¹⁾

1)中国地震局地质研究所, 地震动力学国家重点实验室, 北京 100029

2)中国科学院青藏高原研究所, 大陆碰撞与高原隆升重点实验室, 北京 100101

*通讯作者: 汤吉, 男, 1956年生, 研究员, 博士生导师, 主要从事大地电磁理论、大地电磁网(Network-MT)的理论和应用研究, E-mail: tangji@ies.ac.cn。

〔作者简介〕阮帅, 男, 1983年生, 2015年于成都理工大学获固体地球物理学专业博士学位, 现为中国地震局地质研究所博士后, 主要从事地球物理电磁法数值模拟和三维反演技术研究, E-mail: marshall.ruanel@qq.com。

收稿日期: 2020-01-20

基金项目: 国家自然科学基金(41704078, 41674081)资助

摘要

长白山天池火山是具有潜在喷发危险的巨型活火山, 对其深部低阻岩浆囊的精确探测和准确定位是人们关注的一个重要问题。应用大地电磁方法研究深部电阻率结构是确定岩浆囊的最有效途径之一。文中采用自主研发的三维大地电磁自适应正则化拟牛顿反演程序, 对长白山天池火山区大地电磁数据进行了带地形三维反演, 获取了较合理的深部三维电性结构模型。反演除拟合大地电磁主阻抗外, 同时增加了对横向不均匀性更灵敏的倾子数据, 并通过定性分析对比、模型验证的方法确认了三维反演的电性结构可靠性。三维电性结构模型比以前的二维解释结果更加准确合理, 在有限稀疏的数据集下尽量精细地勾画了壳内低阻岩浆囊及通道的三维空间分布。三维电性结构显示壳内岩浆囊位于天池东北部深10~30km处, 整体上岩浆囊和岩浆通道系统呈 “V”字形分布, 天池火山口的岩浆通道在天池以北10km的EW向剖面图上清晰可见。

关键词: 长白山天池火山; 电性结构; 大地电磁; 拟牛顿反演; 自适应正则化

中图分类号:P317.2 文献标志码:A 文章编号:0253-4967(2020)06-1282-19

ELECTRIC STRUCTRUE MODEL OF TIANCHI VOLCANO IN CHANGBAI MOUNTAINS BASED ON THREE-DIMENSIONAL AR-QN MAGNETOTELLURIC INVERSION

RUAN Shuai¹⁾, TANG Ji¹⁾, DONG Ze-yi¹⁾, WANG Li-feng¹⁾, DENG Yan²⁾, HAN Bing¹⁾

1)State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China

2)Key Laboratory of Continental Collision and Plateau Uplift Institute of Tibetan plateau Research,Chinese Academy of Sciences, Beijing 100101, China

Abstract

The Tianchi volcano in Changbai Mountains is a giant active volcano with potential risk of eruption, so more detailed researches of deep magma chamber as well as its closely related three-dimensional electric structure model are being concerned by many scholars. Based on adaptive regularized three-dimensional quasi-Newton inversion program developed by the author, this paper carries out three-dimensional magnetotelluric inversion including topography using the existing data observed decades ago. Our three-dimensional magnetotelluric adaptive regularization quasi-Newton inversion method improves the conventional quasi-Newton method by approximating Hessian matrix of the data misfit using LBFGS formula instead of total Hessian matrix which is approximated in conventional quasi-Newton inversion. This not only guarantees the precise regularization-term Hessian matrix, but also allows the regularization parameter to be adaptively updated on every inversion iteration without destroying the descent of total objective function. Thus, the regularization parameter's adaptive updating strategy on every iteration can be established based on L2-norm ratio of the data misfit function and regularization function, and our AR-QN inversion algorithm was implemented. The model synthetic data inversion results indicate that AR-QN inversion has very stable iteration flow, and is weakly dependent on initial model selection, which tends to be more suitable for the three-dimensional inversion task when MT survey sites are sparsely distributed on the surface.
Besides the impedances data, the tipper data which is more sensitive to horizontal conductive discontinuity was fitted in this three-dimensional inversion as well. Through comparison of qualitative analysis to measured data and anomalous-body-erasing forward modeling test, our new three-dimensional electric structure was proved to be a reliable model under the constraints of current magnetotelluric data. The three-dimensional electric structure is more accurate and rational compared with the electric structure obtained by previous researches based on two-dimensional inversion, and clearly characterizes the distribution of crustal magma chamber and magma channels with high resolution, even inverted by the limited spare-sites-distributed magnetotelluric dataset. Our result locates the crustal magma chamber on the northeast of Tianchi within the depths of 10~30km, and the general magma chamber along with magma channels system is represented as “V” shape. The magma channel of Tianchi volcano can be clearly shown on the EW profile 10km away from the north of Tianchi volcano. The magma chamber and channels system can be divided into three layers from shallow to deep: Magma channel with a depth of less than 5km and connecting to the Tianchi crater; alkaline fluid magma layer with a depth of 5~10km; and trachytic magma layer with a depth of 10~30km. As the coverage of the existing MT survey sites is still incomplete on the whole Tianchi volcano area, adding more MT survey sites in the northwest, southwest directions and especially in North Korea will help us get more reliable, higher-resolution three-dimensional electrical structure model for Changbaishan Tianchi volcano area.

Keyword: Tianchi volcano in Changbai Mountains; electric structure; magnetotelluric; quasi-Newton inversion; adaptive regularization

文章图片

0 引言

长白山天池火山区位于中国吉林省东部的中朝边界地带, 是长白山山脉的最高峰, 鸭绿江、图们江和松花江的发源地。根据岩石年代学、岩石化学的研究结果, 天池火山的演化可分为 “ 玄武岩造盾” 、 “ 粗面岩造锥” 和 “ 碱流岩爆发” 3个发展阶段(樊祺诚等, 2006)。近代最强烈的一次喷发发生在距今约1 000a之前(刘若新等, 1998; 魏海泉, 2014), 在日本北海道南部及九州岛北部发现的该次喷发的火山灰厚达2~5km(Machida et al., 1983, 1990), 这次喷发也被称为 “ 千年大喷发” 。 2002—2005年天池火山的频繁活动引起了国内外学者的广泛关注(刘国明等, 2006)。自2010年起, 温泉水温、火山区形变监测结果出现异常(冯恒栋等, 2014); 2014年后, 火山气体He、 H₂和CO₂的监测结果均出现异常(刘国明等, 2018), 表明天池火山活动可能更加剧烈。 20世纪末, 中国地震局 “ 八五” 、 “ 九五” 重点项目在天池火山区进行了61个大地电磁测深(MT)点的观测, 并对数据进行了倾子感应矢量的定性解释(汤吉等, 1999)和二维快速松弛(RRI)反演解释(汤吉等, 2001), 发现了位于天池北、东地区深10~30km处存在低阻异常体, 认为其与壳内岩浆囊对应。 2014年, 国土资源部公益性行业科研专项与国家 “ 深部探测技术与实验研究” 专项再次对研究区的南北线进行了MT的数据观测和反演, 也证实了天池北部存在低阻异常体, 但低阻体上方(约深15km)的局部形态与早期研究不同(仇根根等, 2014)。由于受观测点分布稀疏和当时MT反演方法技术的限制, 前期2次研究均使用了当时较先进的二维反演算法。然而, 二维反演假定电阻率沿走向不变, 参与拟合的不同极化模式数据选择、主轴旋转角的估计等都存在一定的人为性, 加上偏离MT剖面的测点投影方式、参与拟合的MT响应类型、所用反演方法、是否考虑地形等各不相同, 使得2次研究得到的电性结构存在差异。但2次研究的结果都支持天池东北深部的低阻体为岩浆囊的观点。除MT资料的证据外, 主动源深地震探测(张先康等, 2002)发现天池火山口9km以下P波出现低速异常, 背景噪声成像(王武等, 2017)显示天池火山口下方深9~30km的范围内存在2%~7%的S波低速异常, 这种壳内P波、 S波低速异常应该是熔体存在的表现(Watanabe, 1993; 李天觉等, 2019)。截至目前, 天池火山口下方深10~30km的范围内存在低阻、高波速比、低波速的壳内岩浆囊是该区地球物理研究的共同看法, 但在岩浆囊的精确定位以及是否存在岩浆通道及其分布形状等方面仍存在一些争议。

尽管前期基于MT研究的二维反演剖面上显示有对应岩浆囊的低阻异常体, 但实际的地下电性结构十分复杂, 特别是长白山天池火山区, 不仅在地形上表现出强三维性, 而且地下介质也为明显的三维电性结构。人们对MT的三维模型用二维近似反演解释进行的研究(Berdichevsky et al., 1998)表明, 二维反演不一定能很好地反映地下三维体真实的电阻率空间分布, 早期长白山天池火山二维MT反演解释的准确性必然受到了这一效应的影响。获取更可靠解释结果的途径是对研究区MT观测数据进行三维反演, 求取能够全局拟合研究区所有MT观测数据的整体三维电性结构模型, 尽量压制因物理假设不适用火山构造、参数设置具有人为性、对各测线进行局部拟合的二维反演所带来的误差。

近年来, 随着计算机硬件和MT理论的发展, 三维反演已逐渐走向成熟。数据空间OCCAM反演(Siripunvaraporn et al., 2005)迭代稳定, 不太依赖初始模型, 但需要显式计算并存储雅可比矩阵, 内存消耗巨大且计算速度慢, 不太适应大型网格的全频带实测数据反演。非线性共轭梯度(NLCG)方法对内存的需求小、计算速度快(Newman et al., 2000), 但实际应用表明其更依赖于初始模型, 在正则化因子和初始模型选择不合理时往往会得到很多假异常体。为了结合两者优点, 我们自主开发了自适应正则因子的拟牛顿(QN)三维MT反演算法(AR-QN)对长白山天池火山区的实测数据进行反演。 AR-QN基于有限内存拟牛顿(LBFGS)更新公式(Byrd et al., 1994)单独逼近数据拟合项的海塞阵, 避免了常规NLCG反演正则罚项海塞矩阵的精度损失, 按拟合项与正则项梯度二范数平权的原则调整正则因子。理论模型合成数据的反演测试表明该方法迭代稳定, 对初始模型的依赖性不强, 反演结果受观测数据噪声的影响低于常规QN和NLCG反演方法(邓琰等, 2019)。

本文利用AR-QN三维反演方法对长白山天池火山区的61个MT测点资料重新进行了反演解释, 并与前期的二维反演剖面进行对比, 对长白山天池火山区岩浆囊系统做出更精细可信的解释。

1 数据定性分析及反演数据选择

长白山天池火山区是重要的风景名胜区, 随着景区的开发, 研究区内各种电磁干扰越来越强。 1997年的MT观测(点位如图 1 所示)电场信号受环境干扰严重, 主要表现在全区0.1~1Hz频段大多出现视电阻率曲线呈45°上升、阻抗相位曲线趋向于0的近源干扰现象(图 2 给出了几个典型测点的近源干扰情况)。干扰频段的倾子数据质量大多优于阻抗数据, 倾子的实部、虚部曲线比较光滑, 数值大多在-0.5~0.5的可信范围, 据此特征前期研究主要使用倾子感应矢量对研究区的电性结构进行定性解释(汤吉等, 1999)。随着周边工业和旅游设施的持续开发, 2014年的MT观测干扰更加严重。从已发表的数据曲线看, 不仅存在视电阻率曲线超45°上升的现象, 且100~1 000s的低频带也受到了强干扰(如C03和C14测点)(仇根根等, 2014)。在反演之前如不剔除这些干扰数据, 有可能会因拟合近源效应的视电阻率而出现大量假高阻体。

	Figure Option View Download New Window
	图 1 长白山天池火山区大地电磁测点分布图图中典型MT测点的数据曲线见图2Fig. 1 Topography of the study area and MT sites around Tianchi volcano in Changbai Mountains.

	Figure Option View Download New Window
	图 2 长白山天池火山区近源干扰的典型MT主阻抗和倾子数据曲线Fig. 2 Typical near-field noise influenced MT primary impedance and tipper curves of Tianchi volcano in Changbai Mountains.

前期研究表明, 该区0.01~1Hz(按100Ω ·m背景电阻率估计探深大致为5km)的二维偏离度偏高, 周期为5.7s的单频倾子实感应矢量指示天池火山口东北部存在大型低阻体(汤吉等, 1999)。为减小电磁干扰的影响, 我们对全区数据重新进行了分析, 剔除了干扰较大的倾子频点, 计算了高、中、低3个频带的平均倾子实感应矢量(图 3), 中频段(按100Ω ·m计算的趋肤深度为5~30km)的计算结果再次证实了天池东北存在深部低阻体(图3b)。

	Figure Option View Download New Window
	图 3 长白山天池火山区高、中、低频段平均倾子的实感应矢量(箭头指向低阻)Fig. 3 Tipper real inductive arrows of Tianchi volcano in Changbai Mountains in three high, middle and low frequency bands(arrows point to conductor).

前期二维偏离度分析和平均倾子实感应矢量都表明研究区电性结构的三维性很强。从图 2 所示的MT视电阻率曲线也可知天池火山区的地下电阻率变化范围非常宽, 给反演初始模型的选择带来困难, 此时常规的NLCG和QN反演方法很难得到可靠的反演结果。我们分析了AR-QN方法和常规方法在原理上的异同, 用模型合成数据测试测点分布稀疏、观测频点较少时的模型恢复效果, 以评价AR-QN是否更适合长白山天池火山区的MT数据反演。

2 AR-QN三维MT反演方法的原理及模型数据反演测试

一般来说, 三维MT反演问题是求解式(1)的最优化问题:

Φ = $\frac{1}{2}$ ||W_d[F(m)-b]|| $_{2}^{2}$ + $\frac{λ}{2}$ ||W_m(m-m_ref)|| $_{2}^{2}$ (1)

其中, 等式右边第1项为数据拟合项, F(m)为模型向量m的三维MT正演响应, b为观测数据向量, W_d一般为对角阵, 各元素为相应观测数据的预期观测误差的倒数; 第2项为正则罚项, m_ref为先验参考模型, 正则罚项对m-m_ref起全局约束作用, W_m为正则算子, 本文取单位梯度算子(最平坦模型期望); 超参数λ 被称为正则因子(或正则参数), 控制反演解偏向于拟合数据还是尽量相似于正则期望。由于Φ 为经典加权非线性最小二乘形式, 一般用迭代法求其极小值。对于第k次迭代, 构造使目标函数下降的方向矢量d_k(被称为下降方向), 则第k+1次的迭代点可用式(2)计算:

m_k₊₁=m_k+α d_k (2)

其中, α 为迭代步长, 一般通过满足Wolfe条件的不精确线搜索确定(Moré et al., 1994)。下降方向d_k有不同的构造方式, 分别对应不同的反演方法, 如非线性共轭梯度(NLCG)、高斯牛顿(GN)、和常规型拟牛顿(QN)等。

NLCG利用共轭方向原理建立式(3):

d_k=-g_k+β _kd_k_-₁(3)

其中, g_k为总目标函数梯度, β _k为根据共轭方向构建的标量(如Ploak-Ribiere-Polyar公式), 若F(m_k)对m_k的复雅可比矩阵为J_k, g_k可按式(4)求取(Newman et al., 2000):

$\{\begin{array}{l} {g^{d}}_{k} = Re {{J^{T}}_{k} {W^{T}}_{d} W_{d} [F (m_{k} {) - b]}^{*}} \\ {g^{m}}_{k} = {W^{T}}_{m} W_{m} (m_{k} - m_{ref}) \\ g_{k} = {g^{d}}_{k} + λ {g^{m}}_{k} \end{array}$ (4)

其中, Re{…}表示取实部, 上角标的星号为取共轭, g_k可利用 “ 拟正演” 求取以避免矩阵J_k的显式计算(Rodi et al., 2001)。 NLCG算法因此计算速度极快, 内存需求很小, 其每次迭代需要一次附加的 “ 拟正演” 计算, 用式(5)确定步长α 的线搜索初值:

α ₀= $\frac{{d^{T}}_{k} g_{k}}{{(W_{d} J_{k} d_{k})}_{2}^{2} + λ d_{k}^{T} {W^{T}}_{d} W_{d} d_{k}}$ (5)

牛顿型方法(如GN和常规QN等)在构建d_k时基于Φ 函数在点m_k处的二阶泰勒展开:

Φ (m_k₊₁)≈ Φ (m_k)+(m_k₊₁-m_k)^Tg_k+ $\frac{1}{2}$ (m_k₊₁-m_k)^TB_k(m_k₊₁-m_k)(6)

其极小值要求式(6)对m_k₊₁的偏导数为0, 从而得到牛顿法的下降方向方程:

B_kd_k=( ${B^{d}}_{k}$ +λ ${W^{T}}_{m}$ W_m)d_k=- ${g^{d}}_{k}$ -λ ${g^{m}}_{k}$ (7)

其中, ${B^{d}}_{k}$ 为数据拟合项的海塞矩阵, 在三维MT反演中难以求取。通过一阶近似, 牛顿法可以退化为GN, 其下降方向方程为

[Re( ${J^{T}}_{k} {W^{T}}_{d}$ W_d $J_{k}^{*}$ )+λ ${W^{T}}_{m}$ W_m]d_k=- ${g^{d}}_{k}$ -λ ${g^{m}}_{k}$ (8)

常规QN的下降方向为

${\hat{B}}_{k}$ d_k=- ${g^{d}}_{k}$ -λ ${g^{m}}_{k}$ (9)

其中, ${\hat{B}}_{k}$ 通常为总海塞的有限内存QN逼近(Byrd et al., 1994), 它可基于前n次(3≤ n≤ 20)迭代的向量对s和y进行高效计算。 s和y的表达式为

$\{\begin{array}{l} s_{k} = m_{k + 1} - m_{k} \\ y_{k} = g_{k + 1} - g_{k} \end{array}$ (10)

从迭代次数看, GN反演收敛最快, 常规NLCG反演最慢。但GN(或OCCAM)每次迭代需显式计算J_k, 这需要大量的 “ 拟正演” 操作(Siripunvaraporn et al., 2005)。因此, 从计算时间上比较, 三维MT反演的常规QN和NLCG方法远快于OCCAM反演。

但如果从实际应用中的迭代稳定性和初始模型依赖度上进行比较, 常规NLCG和QN方法要比OCCAM方法逊色许多。原因在于若迭代中保持λ 不变, 则正则项在早期几乎对目标函数无贡献, 当初始模型选取不合理时, 很可能一开始就因强行拟合数据而进入了奇异点。当拟合项逐渐减小时, 会更倾向于降低正则项而导致拟合项发散, 造成迭代不稳定, 甚至陷入局部奇异点。

OCCAM方法更加稳定的主要原因在于在式(8)保留了正则项海塞矩阵的精确形式, 每次迭代的λ 可以不同而并不破坏目标函数下降性, 且算法在每次迭代中自动搜索拟合项下降到期望水平的最大λ , 以保证每次迭代的线性子问题都 “ 充分正则” , 从而提升了收敛稳定性并降低了对初始模型依赖度。常规QN和NLCG方法的下降方向方程则不具备这一能力。

因此, AR-QN方法改进了常规QN的式(9), 利用向量 ${y^{d}}_{k}$ = ${g^{d}}_{k + 1}$ - ${g^{d}}_{k}$ 而非y_k单独逼近数据拟合项的海塞矩阵, 形成的改进型下降方向方程为

[ ${\hat{B}}_{d}^{k}$ +λ ${W^{T}}_{m}$ W_m]d_k=- ${g^{d}}_{k}$ -λ ${g^{m}}_{k}$ (11)

这种做法的好处有: 1)保留了正则项海塞矩阵的精确形式, 总海塞矩阵的逼近精度更高; 2)若 ${\hat{B}}_{d}^{k}$ 是QN的成功逼近, 正则因子可如OCCAM那样灵活更新。 AR-QN方法借鉴了 “ 比值法” (陈小斌等, 2005), 随迭代根据 ${g^{d}}_{k}$ 和 ${g^{m}}_{k}$ 的二范数之比确定正则因子(邓琰等, 2019)。当总目标函数收敛到极小时, ${g^{d}}_{k}$ +λ ${g^{m}}_{k}$ → 0, 且 $({g^{d}}_{k})$ =λ $({g^{m}}_{k})$ , 可认为这样的模型是平衡数据拟合差极小与正则罚函数极小2个条件的期望解。

为验证AR-QN反演算法的迭代稳定性和收敛速度, 建立一个简单模型(名为CUBES)检验算法在测点稀疏、频点较少时的反演效果。 CUBES模型在100Ω ·m均匀半空间中含有1个1 000Ω ·m高阻体和1个10Ω ·m低阻体, 其三维网格剖分和“ 观测频率” 如表1所示, 网格电阻率分布及测点位置如图 4 所示。对其正演MT响应的全阻抗和倾子数据进行了4个反演测试, 反演均沿用正演网格, 按视电阻率相对误差5%、倾子绝对误差0.01设置各参与拟合数据的误差门槛值(目标函数的数据误差矩阵W_d), 其他反演参数见表2。表2中, CB-10、 CB-100和CB-1000用于测试AR-QN算法对初始模型的依赖性, CB-FIX用于测试初始模型偏离真实模型背景时固定正则因子的反演效果, CB-1000和CB-FIX的对比结果可以分析自适应正则因子策略在初始模型偏离真实背景时是否能提升模型的恢复效果。

表1 CUBES模型三维网格剖分参数及合成数据频率表 Table1 3D mesh parameters and frequency list for CUBES model

	Figure Option View Download New Window
	图 4 CUBES模型的网格电阻率及测点分布示意图黑点表示MT测点, 红色虚线为水平、垂直切片的位置Fig. 4 Mesh resistivity and sites' distribution of CUBES model.

表2 CUBES模型的AR-QN反演测试参数对照表 Table2 Inversion parameters for CUBES model synthetic data AR-QN test

4个反演测试的收敛曲线如图 5 所示。 CB-10、 CB-100和CB-1000的RMS、平均粗糙度曲线很光滑(见图5a—c), 均在约第17次迭代时快速稳定地收敛到RMS=1, 且三者的最终平均模型粗糙度和正则因子基本相同, 表明AR-QN反演方法有很强的迭代稳定性。而CB-FIX在前5次迭代时RMS快速减小, 在随后的第10~30次迭代中陷入停滞, 一直到第48次迭代才使RMS收敛到1(图5d)。在正则因子固定时, 反演模型的平均粗糙度从一开始就急剧升高, 随着迭代的进行出现剧烈变化, 直到最终也没有降低到前3个反演测试的水平。

	Figure Option View Download New Window
	图 5 CUBES模型合成数据的AR-QN反演收敛曲线Fig. 5 Convergence curves of CUBES model synthetic data AR-QN inversions.

4个反演测试的迭代初期模型和最终结果如图 6 和图 7 所示。 CB-10、 CB-100和CB-1000模型在早期比较光滑(图6a—c), 每次迭代都是在近 “ 充分正则” 的模型空间里降低拟合差。而CB-FIX反演一开始就因初始模型偏离背景、测点分布稀疏而陷入奇异解, 虽然经过3倍于CB-1000的迭代次数后勉强能恢复2个异常体的位置, 但是反演结果的冗余构造(图7d)明显多于前3个反演。

	Figure Option View Download New Window
	图 6 CUBES模型合成数据的AR-QN反演初期迭代模型各子图的电阻率切片绘制方式同图4Fig. 6 Early iteration model of CUBES model synthetic data AR-QN inversions.

	Figure Option View Download New Window
	图 7 CUBES模型合成数据的AR-QN反演最终结果各子图的电阻率切片绘制方式同图4Fig. 7 Final results of CUBES model synthetic data AR-QN inversions.

在其他更复杂的模型、不同噪音水平的数据下, AR-QN和常规NLCG三维MT反演算法的对比成果也已发表(邓琰, 2019), 因篇幅所限不再赘述。模型合成数据的反演测试表明, AR-QN相对于常规方法具有迭代高效稳定、不太依赖初始模型的优点。由于现阶段OCCAM反演仍然受到内存、计算时间的限制而无法在大、中型三维MT问题中使用, 故对于长白山天池火山区这种测点稀疏不均匀分布且电性结构复杂而难以确定初始模型的MT资料, AR-QN方法比常规QN(或NLCG)方法更适用。

3 研究区实测MT资料反演及可靠性验证

为生成研究区的三维反演网格, 在核心区域横向按1.5km进行等距剖分, 纵向网格剖分从最高地形点向下进行, 全部地形网格按100m进行等间距分层, 小于最高频趋肤深度的1/3(按256Hz, 100Ω ·m计算), 最低地形格以下以1.2倍逐渐扩大至海拔-100km。为满足三维正演的边界条件, 在侧向和垂向按等比扩展附加网格, 具体参数如表3所示。

表3 长白山天池火山区的三维反演网格扩展参数 Table3 3D mesh extending parameters of Tianchi volcano area in Changbai Mountains

使用Google Earth地形数据得到每个水平格的中心高程, 之上的网格以空气电阻率填充, 设定其不参与反演, 其余网格均填充为100Ω ·m(见第2节的定性分析), 最终得到如图 8 所示的三维带地形反演初始模型。

Figure Option
View Download New Window

图 8 长白山天池火山区三维MT带地形反演网格及初始模型示意图
白色为空气网格; 横纵坐标为网格编号; 色标为常用对数电阻率。a 海拔1.58km的水平切片; b 海拔1.28km的水平切片; c SN向穿过天池火山口的剖面; d EW向穿过天池火山口的剖面Fig. 8 3D topographic mesh and initial model of three-dimensional MT inversion for Tianchi volcano in Changbai Mountains.

三维MT反演的输入数据根据第1节的分析结果剔除了受近源干扰的主阻抗、绝对值> 1.0或曲线不光滑的倾子数据, 反演目标函数的误差加权矩阵W_d由原始EDI文件的方差确定, 并按视电阻率相对误差5%、倾子数据绝对误差0.01设定了误差门槛值。经AR-QN三维MT反演迭代100次后收敛。如图9a所示, 实测数据的反演平均粗糙度上升缓慢, RMS平滑稳定降低直至收敛, 每次迭代都是相对光滑的模型修正, 最终收敛解应和CUBES模型的前3个反演结果(图5c, 7c)具有相似的可靠性。

	Figure Option View Download New Window
	图 9 长白山天池火山区的AR-QN三维MT反演收敛曲线 a AR-QN反演, 初始模型为图8; b 固定小正则因子的反演, 初始模型为图a的结果Fig. 9 Convergence curves of AR-QN 3D MT inversion for Tianchi volcano in Changbai Mountains.

由于倾子数据对电阻率的横向突变更敏感, 与正则项要求模型处处光滑的期望相悖, 为了在AR-QN反演结果的邻域里找到分辨率更高的模型, 可以AR-QN的反演结果为初始模型进行固定小正则化因子(根据试验选择为0.001)的常规QN反演, 其收敛曲线如图9b所示。与图9a相比, 模型的平均粗糙度从0.5上升至约1.6, RMS大致从3.7降为3.6, 表明经过图9a的100次迭代, 模型已经收敛于观测数据噪音允许的拟合差函数极小值(EDI文件的方差和预设误差门槛和实际噪音不等价)。将正则因子降低至0.001之后, 模型的平均粗糙度增加了约3倍, 而RMS只降低了0.1(图9b)。在同样的拟合差水平下, 更简单的模型往往更可靠, 为避免出现过多的冗余构造降低模型的可靠性, 在67次迭代后中止反演。

最终得到的长白山天池火山区三维电性结构如图 10 所示。三维电性结构中的壳内低阻体分布与倾子感应矢量(图 3)的定性分析一致, 深部低阻体位于天池东北, 10Ω ·m低阻等值面更清晰地显示了 “ 囊状” 和 “ 通道状” 的低阻体, 它的形状及和围岩的接触关系都表明其应是壳内岩浆囊系统。 “ 通道” 状低阻体在原东测线以北约10km处(图10b), 这是前期二维反演无法给出的电性结构信息, 也是三维反演的优势所在。

	Figure Option View Download New Window
	图 10 长白山天池火山区的AR-QN三维反演地壳电性结构示意图 a 海拔-15km的切片组合10Ω ·m等值面; b 天池以北约10km的剖面组合10Ω ·m等值面Fig. 10 Crustal electrical structure from AR-QN 3D MT inversion for Tianchi volcano in Changbai Mountains.

为进一步验证三维反演结果的可靠性, 尤其是位于天池东北部可能与壳内岩浆囊对应的深度为10~30km的低阻体的可信度, 将图 10 中红色等值面内深度> 10km的区域电阻率重新填充为100Ω ·m, 得到验证模型, 进行三维正演计算, 以得出其理论MT响应。实测数据、三维电性结构的正演响应、验证模型的正演响应如图 11 所示。从反演的拟合情况看(实线), 近源干扰较大的测点N12在1~10Hz的yx向视电阻率拟合较差(图11d), 其他基本都在可接受范围内。验证模型的理论响应在N12测点(天池以东)、 NE04测点(天池东北)明显受到 “ 抹去” 低阻异常的影响(图11c, d), 视电阻率、倾子数据的拟合情况明显变差。

	Figure Option View Download New Window
	图 11 长白山天池火山典型测点的三维MT反演数据拟合及低阻体模型验证曲线Fig. 11 3D AR-QN MT inversion data fitting on typical MT sites of Tianchi volcano in Changbai Mountains and theoretical response of conductive-body-erased verification.

图10显示, 壳内低阻体总体呈 “ V” 字形结构, 在5~15km深度范围内形成2个分叉, 这一形状是否真实将影响对壳内岩浆囊的进一步解释。为此分别再对2个分叉进行了上述验证过程。计算了每个分叉上方测点在 “ 抹去” 低阻后的数据拟合情况: 1)“ 抹去” 天池火山口一侧的分叉后, 局部测点平均拟合差从2.5增加到3.3; 2)“ 抹去” 东北角的分叉后, 局部测点的平均拟合差从3.7增加到4.0。

通过模拟验证, 结合倾子实感应矢量(图3b)的定性分析结果可知, 三维电性结构揭示的低阻体位置和形状整体可靠。观测数据对天池火山口一侧的低阻分叉很灵敏, 且因其位于测区中部, 数据约束完备, 故其结构信息应是可信的。虽然观测数据对NE侧的低阻分叉灵敏性稍弱, 但模型验证结果也证明了其存在。考虑第2节的CUBES模型CB-FIX反演算例(图 7d), 当初始模型偏离背景、正则因子选择不合理时, 测区边缘的反演结果可能因数据约束不足而影响真实形态。故仅通过现有观测数据只能确认NE侧存在低阻分叉, 其更精细的电性结构需加密、外扩MT测点并重新进行三维反演后才能更可靠地测定。

4 与前期研究成果的对比及三维电性结构的解释

理想情况下, 二维结果和三维结果的同测线剖面往往整体相似, 仅在细节上存在差异。将第3节得到的三维电性结构沿原二维反演测线切片, 图 12 和图 13 对比了早期二维反演解释结果与三维电性结构的不同。

Figure Option
View Download New Window

图 12 长白山天池火山区三维电性结构和前期二维反演结果对比(EW向)
EW1剖面(红色)为2001年的二维RRI反演剖面线; EW2剖面(绿色)为过三维结构低阻体中心的剖面线。a 剖面位置示意图; b EW1剖面的前期二维RRI反演结果(汤吉等, 2001); c 三维反演结果的EW2剖面切片; d 三维反演结果的EW1剖面切片Fig. 12 Comparison of 3D electrical structure and previous 2D inversion results of the Tianchi volcano in Changbai Mountains(EW direction).

Figure Option
View Download New Window

图 13 长白山天池火山区三维电性结构和前期二维反演结果对比(SN向)
SN1剖面(红色)为2001年的二维RRI反演剖面线; SN2剖面(绿色)为过三维结构低阻体中心的剖面线。 a 剖面位置示意图; b SN1剖面的前期二维RRI反演结果(汤吉等, 2001); c 三维反演结果的SN2剖面切片; d 三维反演结果的SN1剖面切片Fig. 13 Comparison of 3D electrical structure and previous 2D inversion results of the Tianchi volcano in Changbai Mountains(SN direction).

EW1剖面二维结果(图12b)在天池以东5km、 10~30km深度范围内显示出一个孤立低阻体, 被解释为壳内岩浆囊(汤吉等, 2001)。三维反演结果(图12d)与二维结果形似, 但电阻率和背景差异极小。而在EW1剖面以北5km处的EW2剖面(图12c)上岩浆囊C(虚线框内)的显示更加明显, 其位置相对EW1的二维结果的位置更靠东。另外, EW2剖面在深度15km以上出现了通道状的低阻结构P1。通过对比分析, 原二维反演结果上的岩浆囊应为EW2剖面附近低阻的旁侧效应, 三维反演压制了这一效应, 在恢复岩浆囊C真实位置的同时, 还揭示了其与低阻体P1(可能为岩浆通道)的接触关系。

图 13a所示的SN1剖面是由2段测线组成的折线, 沿此折线的前期二维反演结果(图13b)和三维结果模型的垂向切片(图13d)均反映在天池火山口北侧存在壳内低阻体, 但二者的形状和范围有所区别。在二维剖面上, 该低阻体为一个封闭孤立的块体, 与底面的连接性不明显。 SN1剖面以东的SN2剖面(图13c)上低阻体的形态和二维RRI反演结果类似, 这同样表明二维反演结果受到了旁侧效应的影响。

为更好地分析C、 P1、 P2这3个低阻体的形态, 对三维电性结构进行不同深度的水平切片, 如图 14 所示。

	Figure Option View Download New Window
	图 14 长白山天池火山区三维电性结构不同深度的水平切片Fig. 14 Horizontal slices of 3D electrical structure of different depths in Tianchi volcano of Changbai Mountains.

图14清晰地展示了天池东北低阻体的横向分布及其与围岩的接触关系。在深度为30km(图14c)的水平切片上, 低阻体C的核心位置位于测线SN1的NE向。由于此区域MT观测点较多, 且在第3节进行了低阻体可靠性验证, 可以确定在现有MT数据点的约束下壳内低阻体底部水平位置应比较准确。随着深度减小, 在深15km(图14b)、 10km(图14a)处低阻体逐渐分开为P1、 P2 2个独立结构。

结合EW2剖面(图12c)、 SN1剖面(图13d)及低阻等值面的形状(图10a), 并综合该区已有的大地电磁、主动源深地震探测、背景噪音成像等研究成果进行分析可知, 低阻体C为长白山天池火山的壳内岩浆囊。在岩浆囊上方, 低阻体P1从深到浅逐渐延伸至天池火山口, 它应为长白山天池火山与壳内岩浆囊的岩浆通道。现有观测数据仍不足以反映低阻体P2的精细结构, 需综合其他学科的研究成果进行解释。深地震测深揭示的P波速度结构在与P1、 P2平面位置相同、深约8km处发现2个低速异常, 据此认为马鞍山-三道白河断裂带与岩浆的输运、上涌过程密切相关(张先康等, 2002)。天池火山岩的K-Ar同位素年龄表明, P2西南角老房子小山的粗面玄武岩表面年龄(约0.87Ma)比天池火山口附近的碱流岩表面年龄(约0.023Ma)约老0.85Ma(樊祺诚等, 2006)。也有学者根据火山岩的结晶分异研究提出了 “ 分层岩浆房” 的观点(陈晓雯等, 2017), 认为壳内深度5~8km为全新世喷发碱流质岩浆, 深度10~30km为造锥阶段粗面质岩浆。综合以上观点, 可认为在长白山天池火山区 “ 粗面岩造锥” 阶段之后, P2区域5km以上的古岩浆通道逐渐凝固, 受岩浆结晶分异及周边断裂带的互相作用, 在 “ 碱流岩爆发” 阶段P2的深部已成为上层碱流质岩浆的贮存区。

综上所述, 我们利用AR-QN三维MT反演方法得到的三维电性结构应比前期的二维反演结果更准确可靠, 长白山天池火山区的三维电性结构模型显示在天池以北约15km、以东10~15km、深度10~30km的范围内存在低阻岩浆囊。整体岩浆系统大致呈 “ V” 字形分布, 从浅到深可分为3层: 1)深度5km以上为与天池火山口连通的岩浆通道(仅P1可见); 2)深度5~10km为碱流质岩浆层(P1和P2); 3)深度10~30km为粗面质岩浆层(C)。必须说明的是, 目前长白山天池火山区现有的MT测点覆盖程度严重不足, 我们得到的三维电性结构是能总体拟合全部有效MT响应的 “ 最简单” 解, 这是在有限数据集下尽量压制反演多解性的最可靠手段, 但真实电性结构的准确性从根本上必须由MT有效测点的覆盖程度决定。因此, 扩大观测范围, 补充西北、西南、尤其东南朝鲜境内的MT观测, 用更完备的观测数据集进行三维反演, 才能得到更可靠、分辨率更高的长白山天池火山电性结构, 为中国东北火山群的监控、防灾工作提供更准确、全面的地球物理参考资料。

5 结论

本文使用AR-QN三维MT反演方法对长白山天池火山区已有MT观测资料进行带地形反演, 获得了研究区的三维电性结构。主要得到以下结论:

(1)火山区的电性结构往往具有强三维性, 带地形的三维反演能较好地克服这些影响, 更真实地恢复地下电阻率模型。反演前需要仔细挑选合理的MT响应, 剔除受干扰的频点。

(2)对AR-QN三维MT反演方法进行了模型合成数据反演试验, 结果表明该方法迭代稳定、不太依赖初始模型。对于测点分布稀疏、难以选择初始模型的情况, AR-QN的适应性更强。

(3)三维反演克服了测线外的异常体旁侧效应, 利用 “ 伞状” 稀疏分布的有限MT观测点获得了尽量可靠的长白山天池火山三维电性结构。与早期二维反演结果相比, 三维反演考虑了地形影响, 增加了倾子数据的拟合。收敛曲线、定性分析及模型验证表明其具有可靠性。

(4)三维电性结构结果与早期研究结果略有不同, 深10~30km的岩浆囊应位于天池火山口以北约15km、以东10~15km的范围内。岩浆囊和岩浆通道系统整体呈 “ V” 字形分布, 与长白山天池火山区演化进程紧密相关。三维电性结构也清晰地显示了物质向天池火山口运移的通道, 而前期的二维反演结果未见显示。

随着国家对长白山天池景区及周边的不断开发, 电磁干扰将日益严重, 建议应尽快对研究区进行 “ 抢救式” 的MT观测, 在空白区补充MT控制测点, 争取国际合作以获取朝鲜境内的MT资料, 以提升现有三维电性结构的准确性和完备性。

参考文献

文献列表

[1]	陈小斌, 赵国泽, 汤吉, 等. 2005. 大地电磁自适应正则化反演算法[J]. 地球物理学报, 48(4): 937—946. CHEN Xiao-bin, ZHAO Guo-ze, TANG Ji, et al. 2005. An adaptive regularized inversion algorithm for magnetotelluric data[J]. Chinese Journal of Geophysics, 48(4): 937—946(in Chinese). [本文引用:1]
[2]	陈晓雯, 魏海泉, 杨良锋, 等. 2017. 长白山天池火山的岩石矿物学特征: 对结晶分异过程和岩浆混合作用的启示[J]. 地球学报, 38(2): 177—192. CHEN Xiao-wen, WEI Hai-quan, YANG Liang-feng, et al. 2017. Petrological and mineralogical characteristics of Tianchi volcano, Changbai Mountains: Implications for crystallization differentiation and magma mixing[J]. Acta Geoscientia Sinica, 38(2): 177—192(in Chinese). [本文引用:1]
[3]	邓琰. 2019. 玉树地震区三维电性结构及孕震环境研究[D]. 北京: 中国地震局地质研究所. DENG Yan. 2019. The study of three-dimensional electrical structure and the seismogenesis environment in Yushu earthquake area[D]. Institute of Geology, China Earthquake of Academy, Beijing(in Chinese). [本文引用:2]
[4]	邓琰, 汤吉, 阮帅. 2019. 三维大地电磁自适应正则化有限内存拟牛顿反演[J]. 地球物理学报, 62(9): 3601—3614. DENG Yan, TANG Ji, RUAN Shuai. 2019. Adaptive regularized three-dimensional magnetotelluric inversion based on the LBFGS quasi-Newton method[J]. Chinese Journal of Geophysics, 62(9): 3601—3614(in Chinese). [本文引用:1]
[5]	樊祺诚, 隋建立, 王团华, 等. 2006. 长白山天池火山粗面玄武岩的喷发历史与演化[J]. 岩石学报, 22(6): 1449—1457. FAN Qi-cheng, SUI Jian-li, WANG Tuan-hua, et al. 2006. Eruption history and magma evolution of the trachybasalt in the Tianchi volcano, Changbaishan[J]. Acta Petrologica Sinica, 22(6): 1449—1457(in Chinese). [本文引用:2]
[6]	冯恒栋, 南颖, 杜会石. 2014. 长白山天池火山监测与火山灾害研究动态[J]. 延边大学学报(自然科学版), 40(3): 265—268. FENG Heng-dong, NAN Ying, DU Hui-shi. 2014. Current research trends of Changbaishan Tianchi volcano monitoring and volcanic disasters[J]. Journal of Yanbian University(Natural Science Edition), 40(3): 265—268(in Chinese). [本文引用:1]
[7]	李天觉, 陈棋福. 2019. 利用接收函数方法研究中国东北东南部地区不同构造体地壳特征[J]. 地球物理学报, 62(8): 2899—2917. LI Tian-jue, CHEN Qi-fu. 2019. Crustal structure of different tectonic units in southeastern part of Northeast China using receiver functions[J]. Chinese Journal of Geophysics, 62(8): 2899—2917(in Chinese). [本文引用:1]
[8]	刘国明, 孙玉涛, 李婷, 等. 2018. 长白山天池火山气体地球化学监测及数据处理方法探讨[J]. 矿物岩石地球化学通报, 37(4): 621—628. LIU Guo-ming, SUN Yu-tao, LI Ting, et al. 2018. The geochemical monitoring of volcanic gases in the Tianchi area of the Changbaishan Mountains and discussion of their data processing methods[J]. Bulletin of Mineralogy, Petrology and Geochemistry, 37(4): 621—628(in Chinese). [本文引用:1]
[9]	刘国明, 于洪茌, 谭雨文. 2006. 长白山天池火山监测工作进展[J]. 岩石学报, 22(6): 1491—1493. LIU Guo-ming, YU Hong-chi, TAN Yu-wen. 2006. Progress in monitoring the Tianchi volcano, Changbaishan[J]. Acta Petrologic Sinica, 22(6): 1491—1493(in Chinese). [本文引用:1]
[10]	刘若新, 魏海泉, 李继泰. 1998. 长白山天池火山近代喷发 [M]. 北京: 科学出版社. LIU Ruo-xin, WEI Hai-quan, LI Ji-tai. 1998. Neoteric Eruption of Changbai Tianchi Volcano [M]. Science Press, Beijing(in Chinese). [本文引用:1]
[11]	仇根根, 裴发根, 方慧, 等. 2014. 长白山天池火山岩浆系统分析[J]. 地球物理学报, 57(10): 3466—3477. QIU Gen-gen, PEI Fa-gen, FANG Hui, et al. 2014. Analysis of magma chamber at the Tianchi volcano area in Changbai Mountains[J]. Chinese Journal of Geophysics, 57(10): 3466—3477(in Chinese). [本文引用:2]
[12]	汤吉, 邓前辉, 赵国泽, 等. 2001. 长白山天池火山区电性结构和岩浆系统[J]. 地震地质, 23(2): 191—200. TANG Ji, DENG Qian-hui, ZHAO Guo-ze, et al. 2001. Electric conductivity and magma chamber at the Tianchi volcano area in Changbaishan Mountains[J]. Seismology and Geology, 23(2): 191—200(in Chinese). [本文引用:2]
[13]	汤吉, 晋光文, 赵国泽, 等. 1999. 感应矢量及其在长白山天池火山区的应用[J]. 地质论评, 45(S1): 294—303. TANG Ji, JIN Guang-wen, ZHAO Guo-ze, et al. 1999. Induction arrow and its application in Tianchi Volcano, Changbai Mountains[J]. Geological Review, 45(S1): 294—303(in Chinese). [本文引用:3]
[14]	王武, 陈棋福. 2017. 长白山火山区地壳S波速度结构的背景噪声成像[J]. 地球物理学报, 60(8): 3080—3095. WANG Wu, CHEN Qi-fu. 2017. The crust S-wave velocity structure under the Changbaishan volcano area in Northeast China inferred from ambient noise tomography[J]. Chinese Journal of Geophysics, 60(8): 3080—3095(in Chinese). [本文引用:1]
[15]	魏海泉. 2014. 长白山天池火山 [M]. 北京: 地震出版社. WEI Hai-quan. 2014. Changbaishan Tianchi Volcano [M]. Seismological Press, Beijing(in Chinese). [本文引用:1]
[16]	张先康, 张成科, 赵金仁, 等. 2002. 长白山天池火山区岩浆系统深部结构的深地震测深研究[J]. 地震学报, 24(2): 135—143. ZHANG Xian-kang, ZHANG Cheng-ke, ZHAO Jin-ren, et al. 2002. Deep seismic sounding investigation into the deep structure of the magma system in Changbaishan-Tianchi volcanic region[J]. Acta Seismologica Sinica, 24(2): 135—143(in Chinese). [本文引用:2]
[17]	Berdichevsky M N, Dmitriev V I, Pozdnjakova E E. 1998. On two-dimensional interpretation of magnetotelluric soundings[J]. Geophysical Journal International, 133(3): 585—606. [本文引用:1]
[18]	Byrd R H, Nocedal J, Schnabel R B. 1994. Representations of quasi-Newton matrices and their use in limited memory methods[J]. Mathematical Programming, 63(1-3): 129—156. [本文引用:2]
[19]	Machida H, Arai F. 1983. Extensive ash falls in and around the sea of Japan from large Quaternary eruptions[J]. Journal of Volcanology and Geothermal Research, 18(1-4): 151—164. [本文引用:1]
[20]	Machida H, Moriwaki H, Zhao D C. 1990. The recent major eruption of Changbai Volcano and its environmental effects[J]. Geographical reports of Tokyo Metropolitan University, (25): 1—20. [本文引用:1]
[21]	Moré J J, Thuente D J. 1994. Line search algorithms with guaranteed sufficient decrease[J]. ACM Transactions on Mathematical Software, 20(3): 286—307. [本文引用:1]
[22]	Newman G A, Alumbaugh D L. 2000. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients[J]. Geophysical Journal International, 140(2): 410—424. [本文引用:2]
[23]	Rodi W, Mackie R L. 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion[J]. Geophysics, 66(1): 174—187. [本文引用:1]
[24]	Siripunvaraporn W, Egbert G, Lenbury Y, et al. 2005. Three-dimensional magnetotelluric inversion: Data-space method[J]. Physics of the Earth and Planetary Interiors, 150(1-3): 3—14. [本文引用:2]
[25]	Watanabe T. 1993. Effects of water and melt on seismic velocities and their application to characterization of seismic reflectors[J]. Geophysical Research Letters, 20(24): 2933—2936. [本文引用:1]

2005

0.0

陈小斌, 赵国泽, 汤吉, 等. 2005. 大地电磁自适应正则化反演算法[J]. 地球物理学报, 48(4): 937—946.

CHEN

Xiao-bin

, ZHAO

Guo-ze

, TANG

, et al. 2005. An adaptive regularized inversion algorithm for magnetotelluric data[J]. Chinese Journal of Geophysics, 48(4): 937—946(in Chinese).

In this paper, a new inversion method, Adaptive Regularized Inversion Algorithm (ARIA), is presented to overcome the difficulty in determination of regularized factors for magnetotelluric (MT) inversion. In ARIA, first a new data variance disposing method, data variance normalization method, is put forward. This method uses a new way to calculate the influence matrix of data variance in inversion. Thus, the data variance only influences data fitting, and has no influence on the weight between the data object function and the model constraint object function. So the influence factors in determination of the regularized factor are reduced. Second, the definition of roughness kernel matrix is presented in the course of constructing the model constraint object function, and a concise equation of it is derived. Thus the construction of the model object function becomes very simple and direct. Third two adaptive methods of regularized factors are put forward based on the relations of data object function, model constraint object function, and regularized factor. ARIA is used to solve the one_dimensional inversion by MT data by the constraint of the flattest model. Several examples are given to illustrate the effectiveness of ARIA.

1 Department of Geophysics, Peking University, Beijing 100871, China 2 Institute of Geology, Earthquake Adminstration, Beijing 100029, China

针对大地电磁正则化反演中正则化因子的选取困难问题提出了自适应正则化反演算法（Adaptive Regularized Inversion Algorithm, ARIA）. 在该算法中， ①提出了一种新的数据方差处理方法：数据方差规范化，使得数据方差的大小只对数据的拟合发生影响，不对数据目标函数和模型约束目标函数的权重产生影响，从而减少了正则化因子取值的影响因素；②提出了粗糙度核矩阵的概念，并给出了由基本结构插值基函数计算粗糙度核矩阵的公式，使得模型目标函数的构建更为简便、直接；③根据数据目标函数、模型约束目标函数和正则化因子之间的关系，提出了两种正则化因子自适应调节方法. 本文详细阐述了最平缓模型约束下的大地电磁一维连续介质反演的ARIA实现，以几个算例的分析比较来说明ARIA的有效性.

... (陈小斌等, 2005), 随迭代根据 gdk和 gmk的二范数之比确定正则因子(邓琰等, 2019) ...

2017

0.0

... 的观点(陈晓雯等, 2017), 认为壳内深度5~8km为全新世喷发碱流质岩浆, 深度10~30km为造锥阶段粗面质岩浆 ...

2019

0.0

... 理论模型合成数据的反演测试表明该方法迭代稳定, 对初始模型的依赖性不强, 反演结果受观测数据噪声的影响低于常规QN和NLCG反演方法(邓琰等, 2019) ...

... 在其他更复杂的模型、不同噪音水平的数据下, AR-QN和常规NLCG三维MT反演算法的对比成果也已发表(邓琰, 2019), 因篇幅所限不再赘述 ...

2019

0.0

... (陈小斌等, 2005), 随迭代根据 gdk和 gmk的二范数之比确定正则因子(邓琰等, 2019) ...

2006

0.0

樊祺诚, 隋建立, 王团华, 等. 2006. 长白山天池火山粗面玄武岩的喷发历史与演化[J]. 岩石学报, 22(6): 1449—1457.

FAN

Qi-cheng

, SUI

Jian-li

, WANG

Tuan-hua

, et al. 2006. Eruption history and magma evolution of the trachybasalt in the Tianchi volcano, Changbaishan[J]. Acta Petrologica Sinica, 22(6): 1449—1457(in Chinese).

In this paper we report new chronological data and chemical analysis results on the Tianchi Volcano,and provide further constrains on the eruption history and the changes of magma types of the shield-forming stage and the cone-forming stage.At the beginning of the early Pleistocene (ca.2Ma),the Tianchi Volcano started the trachybasaltic activity of the shield-forming stage,and then the trachybasaltic magma changed to the trachytic and trachyandesitic magmas at the end of the early Pleistocene (ca.1Ma);the cone-forming stage of trachyte mainly occurred in the middle Pleistocene,and the trachyte magma changed to the pantelleritic magma at the late Pleistocene (ca.0.1Ma).In the middle-late Pleistocene,volcanic activity of the trachyte and pantellerite from the crustal chamber dominated the cone-forming stage,but the episodic volcanism of the trachybasah from the mantle chamber also occurred. Since the mantle chamber supplied trachybash influx from below to the crsutal chamber,the Tianchi Volcano is a long-life volcano. Fractional crystallization and magma mixing were two important processes in the Tianchi Volcano,the former determined the bimodal characteristics of the volcanic rocks,and the latter triggered the volcanic eruption.In the late Pleistocene and the Holocene,volcanic activity of pantellerite dominated the Tianchi Volcano,with episodic activity of basaltic trachyandesite,trachyandesite and trachyte throughout the cone-forming stage.This fact indicated that there is a multi-layered magma chamber system beneath the Tianchi Volcano.The trachybasalte influx from the mantle chamber disordered the magma layers and caused the magma mixing,then triggered the eruptions.

本文新提出的年代学和岩石化学结果，进一步从天池火山与区域火山活动的关系，论述了天池火山造盾、造锥历史和岩浆结晶分异转型的时间约束，早更新世早期（2Ma前）开始粗面玄武岩的造盾，早更新世晚期（约1Ma）粗面玄武岩向粗安岩、粗面岩演化，中更新世是粗面岩造锥的主阶段，到了晚更新世（约0．1Ma）粗面岩向碱流岩演化。在中-晚更新世来自地壳岩浆房的粗面岩、碱流岩造锥过程中，来自地幔的粗面玄武岩浆喷发活动始终没有间断过。由于来自地幔粗面玄武质岩浆持续向地壳岩浆房补给，所以天池火山是一座长寿命的火山。岩浆的结晶分异作用和混合作用是天池火山岩浆演化的两个最重要过程，前者形成天池火山双峰式火山岩分布特征，后者成为天池火山喷发的触发机制。天池火山在晚更新世-全新世碱流质岩浆主喷发期兼有少量玄武质粗安岩、粗安岩或粗面质岩浆的交替喷出，揭示了天池火山的地壳岩浆房熔体的分层结构特点，由于来自地幔粗面玄武质岩浆注入地壳岩浆房，导致不同层位岩浆的扰动和混合作用，触发天池火山的喷发。

... 3个发展阶段(樊祺诚等, 2006) ...

... 85Ma(樊祺诚等, 2006) ...

2014

0.0

... 自2010年起, 温泉水温、火山区形变监测结果出现异常(冯恒栋等, 2014) ...

2019

0.0

... 李天觉等, 2019) ...

2018

0.0

... 2014年后, 火山气体He、 H₂和CO₂的监测结果均出现异常(刘国明等, 2018), 表明天池火山活动可能更加剧烈 ...

2006

0.0

刘国明, 于洪茌, 谭雨文. 2006. 长白山天池火山监测工作进展[J]. 岩石学报, 22(6): 1491—1493.

LIU

Guo-ming

, YU

Hong-chi

, TAN

Yu-wen

. 2006. Progress in monitoring the Tianchi volcano, Changbaishan[J]. Acta Petrologic Sinica, 22(6): 1491—1493(in Chinese).

The Tianchi Volcano, which is located at the border of NE China and North Krea, is one of the most dangerous volcanoes in China. In this paper, we introduced the development in establishment of the Tianchi monitoring station, and reported the progress in recent years, including monitoring data on seismology, on terrain deformation and on geochemistry.

长白山天池火山是我国最具有潜在喷发危险的活火山。本文介绍了天池火山监测系统建设的基本情况，报道了近些年在火山地震监测、大地形变监测以及地球化学监测等方面的工作进展。

... 2002—2005年天池火山的频繁活动引起了国内外学者的广泛关注(刘国明等, 2006) ...

1998

0.0

... 近代最强烈的一次喷发发生在距今约1 000a之前(刘若新等, 1998 ...

2014

0.0

仇根根, 裴发根, 方慧, 等. 2014. 长白山天池火山岩浆系统分析[J]. 地球物理学报, 57(10): 3466—3477.

QIU

Gen-gen

, PEI

Fa-gen

, FANG

Hui

, et al. 2014. Analysis of magma chamber at the Tianchi volcano area in Changbai Mountains[J]. Chinese Journal of Geophysics, 57(10): 3466—3477(in Chinese).

According to South Korean media's comments that the Tianchi volcano in Changbai Mountain is likely to erupt in coming 2 years, we carried out an MT sounding observation profile about 103 km long at Tianchi volcano area for detecting the deep electrical structures. Since the electromagnetic interference from unknown sources is very strong in the research area, we used a number of advanced data processing techniques including remote reference, Robust processing, Rhoplus analysis, impedance tensor decomposition, and a mutual approximate calculation method between apparent resistivity and phase, finally obtained high-quality and more reliable electromagnetic data in the strong interference areas. Then we analysed the two dimensional structural trend and characteristics of induction vector at Tianchi volcano area in Changbai Mountain, using NLCG inversion technique to interpret the inversion data, then compared the inversion results with previous exploration results. The results show that there is an obvious magma channel at the bottom of Tianchi volcano, the magma channel is closed in depth about 5 8 km. At the bottom of the crater, near the north direction about 7 km deep exists an obvious low resistance body, the resistivity is less than 10 Omega m, connected with the magma channel. It is inferred that it may be the developing magma chamber in shallow surface. The abnormal region is between stations C07-009 near the entrance gate of Changbai mountain and between C04-C05 station, there are two obvious nearly vertical low resistivity zones existing in the depth of about 7 similar to 17 km, connecting with below low resistance body directly, it is inferred that the low resistivity zone contains the magma. As the depth increases, the crust low resistivity anomaly bodies are widely developed in the depth of 13km to 30km, north to the position of 20 km south from Tianchi volcano, which may be the active magma chamber. Comparing the inversion results with the previous exploration results, the electrical characteristics are overall similar, but locally different.

[Qiu Gen-Gen; Pei Fa-Gen; Fang Hui; Du Bing-Rui; Zhang Xiao-Bo; Zhang Peng-Hui; Yuan Yong-Zhen; He Mei-Xing; Bai Da-Wei] Inst Geophys & Geochem Explorat, Langfang 065000, Hebei, Peoples R China.;[Qiu Gen-Gen; Pei Fa-Gen; Fang Hui; Du Bing-Rui; Zhang Xiao-Bo; Zhang Peng-Hui; Yuan Yong-Zhen; He Mei-Xing; Bai Da-Wei] Minist Land & Resources, Electromagnet Detect Technol Key Lab, Langfang 065000, Hebei, Peoples R China.

针对外媒报道长白山天池火山在近2年内有可能喷发的言论,在长白山天池火山区布测了一条长度约为103 km的二维大地电磁测深观测剖面,对火山区深部电性结构进行探测研究.由于研究区内不明来源的电磁干扰非常强,对数据采用了远参考处理、Robust处理、Rhoplus分析、张量阻抗分解和基于一维层状介质电阻率与相位互算方法等先进处理技术,获取到一批在强干扰区质量较为可靠的电磁数据,利用数据计算分析了长白山天池火山区二维构造走向和感应矢量特征,采用NLCG二维反演技术对资料进行了二维反演解释,并将反演结果与前人探测结果进行了对比分析.探测结果表明:在天池火山口下方存在明显的直立型岩浆通道,岩浆通道在下方约5~8 km位置形成关闭;在火山口下方往北方向附近,在埋深位置约7 km深处存在一个明显的低阻异常体,电阻率小于10 Ωm,且与岩浆通道对接,推测其可能是地表浅部发育的岩浆囊;在长白山山门附近C07-C09号测点之间和C04-C05号测点之间,在埋深约7~17 km深处发现近直立型低阻带,低阻带与下方低阻体直接相连,推测低阻带内赋存有活动的岩浆;随着埋深的增加,从天池火山口南部约20 km位置往北方向,在埋深13~30 km之间壳内广泛发育明显的低阻异常体,推测其可能是活动的岩浆囊.反演结果与前人探测结果整体电性特征相似,但又局部不同.

... 专项再次对研究区的南北线进行了MT的数据观测和反演, 也证实了天池北部存在低阻异常体, 但低阻体上方(约深15km)的局部形态与早期研究不同(仇根根等, 2014) ...

... 从已发表的数据曲线看, 不仅存在视电阻率曲线超45°上升的现象, 且100~1 000s的低频带也受到了强干扰(如C03和C14测点)(仇根根等, 2014) ...

2001

0.0

汤吉, 邓前辉, 赵国泽, 等. 2001. 长白山天池火山区电性结构和岩浆系统[J]. 地震地质, 23(2): 191—200.

TANG

, DENG

Qian-hui

, ZHAO

Guo-ze

, et al. 2001. Electric conductivity and magma chamber at the Tianchi volcano area in Changbaishan Mountains[J]. Seismology and Geology, 23(2): 191—200(in Chinese).

Observations of magnetotellurics have been carried out at 46 sites along five profiles in shape of umbrella and along a northwest direction profile in the Tianchi volcano area and its vicinity in Changbaishan Mountain northeast China. There are 61 sites including 15 sites which were done in early time. The data processing has been made by using the Robust technique. High S/N ratio resistivity and phase have been obtained. The Bahr’s impedance tensor decomposition was applied for eliminating distortion. The induction arrow vector was calculated. The two-dimensional inversion method (RRI) is used for data interpretation. The results show that (1)There is low resistivity correspondence Tianchi water pool near the surface. There are four low resistivity zones corresponding to Jingjiang hot spring, Julongquan hot spring, gate of mountain and its vicinity and Shuangmu peak at the depth of above 5 km.(2)Low resistivity body at the depth about 12km is located the carter and its north and east area. The body is 20 to 30km thick and 15～25km wide. It is inferred that this anomaly body may be a magma chamber. (3)The crustal electric conductivity structure has obvious difference in direction of SN in this area. There is a high resistivity zone in the south part and low resistivity zone in the north part.

1. Institute of Geology, China Seismological Bureau, Beijing 100029; 2. Zhejiang Exploration Division of Petroleum, Hangzhou 310007; 3. Seismological Bureau of Jilin Province, Changchun 130022

在长白山火山区及其附近,沿5条伞状测线和1条北西向测线进行了61个大地电磁测深点的观测。用Robust等方法对观测资料进行了处理;应用阻抗张量分解方法消除了局部不均匀体的影响;计算了反映地下横向不均匀性的地磁感应矢量。利用RRI二维自动反演等技术对资料进行了解释。结果表明:在近地表存在低阻体,并分别与长白山天池水体、聚龙泉温泉、锦江温泉及长白山山门附近的地下水相对应。在长白山天池及其以北和以东地区,约12km深处存在电阻率很低的地质体,电阻率为10到几十Ω·m,可能是地壳岩浆囊。在测区范围内,地壳电性结构在南北方向存在着明显的横向变化,而东西方向相对均匀,变化比较平缓。

... 重点项目在天池火山区进行了61个大地电磁测深(MT)点的观测, 并对数据进行了倾子感应矢量的定性解释(汤吉等, 1999)和二维快速松弛(RRI)反演解释(汤吉等, 2001), 发现了位于天池北、东地区深10~30km处存在低阻异常体, 认为其与壳内岩浆囊对应 ...

... EW1剖面二维结果(图12b)在天池以东5km、 10~30km深度范围内显示出一个孤立低阻体, 被解释为壳内岩浆囊(汤吉等, 2001) ...

1999

0.0

... 5的可信范围, 据此特征前期研究主要使用倾子感应矢量对研究区的电性结构进行定性解释(汤吉等, 1999) ...

... 7s的单频倾子实感应矢量指示天池火山口东北部存在大型低阻体(汤吉等, 1999) ...

2017

0.0

... 除MT资料的证据外, 主动源深地震探测(张先康等, 2002)发现天池火山口9km以下P波出现低速异常, 背景噪声成像(王武等, 2017)显示天池火山口下方深9~30km的范围内存在2%~7%的S波低速异常, 这种壳内P波、 S波低速异常应该是熔体存在的表现(Watanabe, 1993 ...

2014

0.0

... 魏海泉, 2014), 在日本北海道南部及九州岛北部发现的该次喷发的火山灰厚达2~5km(Machida et al ...

2002

0.0

... 深地震测深揭示的P波速度结构在与P1、 P2平面位置相同、深约8km处发现2个低速异常, 据此认为马鞍山-三道白河断裂带与岩浆的输运、上涌过程密切相关(张先康等, 2002) ...

1998

0.0

... 人们对MT的三维模型用二维近似反演解释进行的研究(Berdichevsky et al ...

1994

0.0

... AR-QN基于有限内存拟牛顿(LBFGS)更新公式(Byrd et al ...

... 其中, B^k通常为总海塞的有限内存QN逼近(Byrd et al ...

1983

0.0

... 魏海泉, 2014), 在日本北海道南部及九州岛北部发现的该次喷发的火山灰厚达2~5km(Machida et al ...

1990

0.0

... , 1983, 1990), 这次喷发也被称为 #cod#x0201c ...

1994

0.0

... 为迭代步长, 一般通过满足Wolfe条件的不精确线搜索确定(Mor#cod#x000e9 ...

2000

0.0

... 非线性共轭梯度(NLCG)方法对内存的需求小、计算速度快(Newman et al ...

... _k为根据共轭方向构建的标量(如Ploak-Ribiere-Polyar公式), 若F(m_k)对m_k的复雅可比矩阵为J_k, g_k可按式(4)求取(Newman et al ...

2001

0.0

... 求取以避免矩阵J_k的显式计算(Rodi et al ...

2005

0.0

... 数据空间OCCAM反演(Siripunvaraporn et al ...

... 操作(Siripunvaraporn et al ...

1993

0.0