2020, Vol. 63
2. 中国科学院大学, 计算地球动力学重点实验室, 北京 100049;
3. 中国地质科学院地质研究所, 北京 100037
2. Key Laboratory of Computational Geodynamics of Chinese Academy of Science, University of Chinese Academy of Sciences, Beijing 100049, China;
3. Lithosphere Research Center, Institute of Geology, Chinese Academy of Geological Sciences, Beijing 100037, China
随着大地测量资料日趋丰富，火山地区地表形变分析计算已成为目前研究地下岩浆囊参数的重要方法.最早的火山区地表形变计算模型是Mogi球状点源模型(Mogi, 1958)，给出了与岩浆房压力变化、深度和半径等参数相关的地表位移解析表达式.为更精确地描述岩浆房及其位移场，不同学者提出椭球状压力变形源计算方法：Davis(1986)提出半空间三轴椭球点状压力源解析表达式，Yang等(1988)提出下倾的有限长旋转椭球压力源计算方法.另外，Okada(1985, 1992)位错公式的扩张断层也可以用来模拟椭球压力源，例如朱桂芝等(2008)利用Okada三正交扩张点源-遗传算法反演长白山火山区变形源.此外，还有其他学者提出例如水平圆形破裂岩浆房压力源(Fialko et al., 2010)和管状压力源(Bonaccorso and Davis, 1999)的地表位移计算方法.但传统解析方法存在很大的局限性--难以考虑地形，尽管有学者提出一些解决办法，例如“变深度模型”(Williams and Wadge, 1998)和“地形更正模型”(Williams and Wadge, 2000).但“变深度模型”局限于对地表垂向位移的修正，无法精确考虑地形对地表水平位移的影响，“地形更正模型”需要满足地表坡度较小的假设，这些限制使得越来越多学者开始使用有限元方法在火山区地表形变中考虑地形和岩浆房的复杂几何(Trasatti et al., 2010; Charco and Del Sastre, 2014).但是传统有限元方法反演岩浆房压力源参数没有很好地解决网格问题，反演过程中每一次正演由于岩浆房位置和大小变化都需要重新生成一次网格，工作量巨大且耗时.据此，我们提出采用有限元“等效体力”方法来数值分析火山岩浆压力变形源，一方面可以克服传统火山区有限元反演每一次正演需要重新生成网格的难题，另一方面可考虑地形影响来分析地下压力变形源引起的地表形变和应力变化，并进一步计算岩浆应力扰动对周边断层稳定性的影响.
长白山天池火山是中国大陆仍在活动的新生代多成因复合火山(魏海泉等，1997；赵大鹏等，2004)，并具有特殊的动力学环境，不同走向断层交汇于此(刘若新等，1998).目前，研究学者已分别从深部结构(汤吉等，2001；张成科等，2002；赵大鹏等，2004；马学英等，2016)、地震活动性(吴建平等, 2005, 2007)、火山气体(Xu et al., 2012)和地表形变等手段(崔笃信等，2007)对长白山火山的活动性进行了研究.例如，在深部结构研究方面：大地电磁测深(汤吉等, 2001)揭示该区地壳存在低电阻层；深地震测深剖面(张成科等, 2002)表明该区壳内存在低速度、密度和电阻异常体；层析成像(赵大鹏等, 2004)显示该区存在显著低速异常.地震活动性研究方面：自2002年7月长白山天池火山地震显著增多，且震中环天池分布(图 1b-d)，深度在3.5~5 km (吴建平等, 2005, 2007)；Hong等(2004)认为这些小震丛集与美国圣海仑斯火山喷发前相似.地球化学研究方面，Xu等(2012)分析了2002-2006年长白山地区的二氧化碳、氦气和氢气等火山气体，获得了火山岩浆活动的证据.地表形变研究方面，长白山天池在2002年汪清地震之后建立了形变监测站(胡亚轩等，2007)，GPS、水准和地倾斜结果(崔笃信等, 2007)显示2002-2003年长白山火山区地表形变表现典型岩浆上涌特征，即垂向位移在靠近火山口处快速隆升，水平位移以火山口为中心向外呈现辐射状衰减，见图 1a，但在2002-2005年该区地震活动逐年减弱(胡亚轩等，2007).PS-InSAR资料显示(Ji et al., 2013)2006-2008年地表形变表现为岩浆侵入的膨胀特征，而2008-2010年表现为反向的沉降特征.
|
图 1 长白山地表形变测量和地震活动性 (a)长白山火山区形变监测; (b)月最大地震数; (c)月最大震级及累计地震矩; (d) 2002-2003年的地震震群. Fig. 1 Surface deformation and seismicity in the Changbaishan volcanic area (a) Monitoring network in Changbaishan volicanic area; (b) Maximum monthly number of earthquakes; (c) Time series of the monthly maximum magnitude and cumulative seismic moment; (d) Earthquake swarms detected by mobile seismic campaign during the summers of 2002 and 2003. |
前人利用不同的大地测量资料和理论模型进行了长白山火山区岩浆房源参数的反演(胡亚轩等，2007；朱桂芝等，2008；陈国浒等，2008；崔笃信等，2007; Ji et al., 2013；李宏卿等，2017)，但这些工作尚未考虑地表地形，并且长白山地下岩浆压力变形源应力扰动对周边断层稳定性影响的研究仍存在空白.为此，本文以长白山火山区为例，采用有限元“等效体力”方法，搭建高三维分辨率数值模型，分析不同地形坡度对地表位移的影响，并采用有限单元-遗传算法反演地下岩浆房压力变形源参数，最后根据岩浆房压力变形源参数计算岩浆应力扰动对周边断层稳定性的影响.
1 地下压力变形源的等效体力正反演方法 1.1 地下压力变形源的有限元等效体力方法模拟地下压力源的有限元方法通常需要在网格划分阶段刻画出岩浆房形状，并在腔壁上施加压力边界条件(Trasatti et al., 2010).该类方法由于腔壁并不是严格的球面，需要尽量地细分单元以逼近球面，并且这样的网格划分容易出现单元形状畸变，对建模要求较高.当用此类方法处理反演问题，由于每一次正演岩浆房的形状、大小发生改变，每一次正演需要重新生成网格或者进行网格重划分.为了克服上述局限，我们采用“等效体力”方法，结合网格自适应加密技术，省去大量网格建模时间，同时能够规避单元畸形问题.
对于球腔深度为d, 岩浆压力增量ΔP的岩浆房，经典Mogi解(Mogi, 1958)可表示为
|
(1) |
|
(2) |
其中，r为观测点相对岩浆源的径向距离，a为岩浆房半径，ΔV为岩浆房等效体积变化，μ为剪切模量.
Mogi模型代表的球状腔体压力源问题可以用3个正交扩张点源或3个正交张拉位错(Mindlin, 1936)来等效(见图 2)，这些概念上不同的源可以产生源外部相同的位移场.考虑如图 1所示的3对正交力偶Fh产生的位移为
|
(3) |
|
图 2 各向同性点源模型(Mindlin, 1936) (a)球腔压力源; (b)三正交力偶; (c)三正交张拉位错. Fig. 2 Isotropic point source models (a) Pressurized spherical cavity; (b) Three orthogonal dipoles; (c) Three orthogonal dislocations. |
其中，Gik(x)是格林函数，代表源处k方向单位力产生的接收点x处的i分量位移，并且矩张量为Mjk=Fhδjk.代入弹性半空间内部集中力的格林函数表达式(Mindlin, 1936)，当岩浆房半径a和岩浆房压力变化ΔP满足
|
(4) |
其中，δx=x′是狄拉克函数，代表x′处(源处)的点源集中力，Δδx=x′代表一对符号相反的脉冲(Burridge and Knopoff, 1964; Aki and Richards, 2002)，每一对力偶的强度可表示为(Bonafede and Ferrari, 2009)：
|
(5) |
由于狄拉克函数可用高斯分布函数的极限来近似，即
|
(6) |
其中，x1、x2、x3为网格积分点坐标，x′1、x′2、x′3为岩浆房点源中心坐标.αxi是高斯函数在xi方向的变量，且αxi取值主要取决于网格大小，αx1αx2αx3π3/2代表了单位体积内高斯函数的正则化系数，这种处理保证了等效体力函数足够光滑.并且，体力函数可以用来模拟不同的岩浆房形状，例如，当αx1=a，αx2=b，αx3=c，且α≠b≠c时，式(6)可以用来模拟椭球状岩浆房压力源，且椭球半轴分别为a、b和c.其他形状的压力源，可以根据弹性力学叠加原理，用式(6)的权重组合形式，用地震矩张量来表示：
|
(7) |
当式(7)应用于火山区地下岩浆房，岩浆房椭球的3主轴可能并不平行于“北-东-下”坐标轴，这时可以通过坐标变换(朱桂芝等，2008)由岩浆房主轴坐标系下的矩张量M变换得到笛卡尔坐标系下的矩张量Mij.
1.2 方法验证为了验证方法有效性，利用开源C++有限元库Deal II(Bangerth et al., 2007)编写并行有限元程序，并采用网格自适应加密技术(Kelly et al., 1983)来使网格疏密程度达到最佳，既优化求解精度又能够平衡计算代价.这里我们提供了一个球形压力源算例并和解析的Mogi模型对比，具体参数如下：岩浆房深度为5 km，半径为1 km，压力差为10 MPa，剪切模量μ=3×104MPa，泊松比为0.25.图 3给出了地表的三分量位移云图，图 3d给出了沿A-A′剖面的有限元解和解析解的对比，可以看出数值计算结果与解析解完全吻合，进而也验证了等效体力方法的正确性和程序编写的可靠性.
|
图 3 有限元与解析解的地表位移对比 (a)东西向位移；(b)南北向位移；(c)垂向位移；(d) A-A′测线有限元解与解析解对比. Fig. 3 Comparison of surface displacements between Mogi′s analytical solution and FEM solution (a) E-W displacement; (b) N-S displacement; (c) Vertical displacement; (d) Comparison of surface displacements between analytical solution and FEM solution in line A-A′. |
地球物理优化问题往往是非线性的，导致反演问题存在许多局部最优，因此一些局部优化算法例如线性化矩阵反演、最速下降法、共轭梯度法等容易过早地收敛于局部最优点，并且最优解极度依赖于初值的选择(Tiampo et al., 2004).遗传算法是具有稳健全局搜索性能的并行算法，具有不需要求目标函数倒数、不依赖初值的优点，已在地震定位(Billings et al., 1994)、速度结构(Bhattacharyya et al., 1999)等地球物理反演中得到广泛应用，越来越多学者采用遗传算法进行火山区岩浆房源参数反演(Fernández et al., 2001; Tiampo et al., 2004; 朱桂芝等，2008).
为改善遗传算法达到最优解区域后收敛速度变慢的问题，我们采用“移民机制”和自适应交叉、变异的改进措施(朱守彪, 2005).在反演中完全随机地产生初始种群，种群规模始终保持32，初始变异概率和交叉概率分别为0.2和0.9，迭代停止步数为100步，采用的目标函数为
|
(8) |
其中，Um是位移模拟值，Ui是形变观测值(GPS和水准测量).经过倒数变换，即构成遗传算法反演模型的适应度函数：Fit=1/F.待反演参数为椭球岩浆房的深度、经度、纬度、主轴长度、主轴走向、主轴倾角.
本文采用有限单元-遗传算法，即正演采用2.1节所述的有限元方法计算.每次正演过程中，在有限元的载荷向量组装阶段，每个单元积分点的体力项采用式(7)计算，然后形成总体刚度矩阵，求解线性方程组，最后得到位移场，GPS和水准测点上的位移和反演计算中的U是相对应的.相比解析方法，有限元计算代价较大，为了平衡计算代价，反演计算中采取逐步逼近方法(石耀霖和Assumpeao, 2000)，即每次迭代步数为100，连续计算10次，再根据这10个最佳模型，缩小搜索区间，给出新一轮反演的搜索范围，并在最终的一轮反演中适当加密网格.本文采用并行有限元技术，我们的所有模拟计算例子均在DELL T7920双CPU/28核的工作站上运行.一次粗网格反演耗时约9 h，细网格反演耗时约54 h.
2 地形影响为了定量化研究地形对膨胀压力源地表形变的影响，选取200 km×200 km×200 km的立方体为计算区域，计算如图 4所示不同深度和不同坡度条件下地形对地表位移的影响，材料参数和岩浆房压力同2.2节.
|
图 4 地形对地表位移的影响 (a)-(c):地表径向位移，坡度分别为10°, 20°和30°；(d)-(f):地表垂向位移，坡度分别为10°, 20°和30°. Fig. 4 Effects of topography on surface displacement (a)-(c) Surface radial displacements with slopes of 10°, 20° and 30°; (d)-(f) Surface vertical displacements with slopes of 10°, 20° and 30°. |
由于火山喷发的物理化学条件不同，因此火山在形态上有多种不同的类型，但火山总体上可以用圆锥来简化描述，我们分析了世界上一些知名火山圆锥底部半径在5~20 km，因此在有限元模拟中固定圆锥底部半径R=10 km，通过调整圆锥高度来改变火山坡度，同时为了减少变量个数，固定岩浆房半径a=1 km.
地形对地表形变的影响主要受埋深和地形坡度控制：当坡度固定，埋深d越浅，a/d越大，浅地表径向位移ur和垂向位移uz受地形影响越大.当a/d固定，坡度越大，地形对地表径向位移ur和垂向位移uz的影响越大.当埋深很浅且坡度很大，例如当a/d=0.5，坡度为30°时，垂向位移uz误差超过200%.
同时地形还影响地表位移最大值的位置，不论是否有地形，地表最大径向位移ur始终不在岩浆囊正上方.由图 4d-f可以看出埋深越浅，地表最大径向位移ur所在的位置越靠近岩浆囊中心，并且坡度越大，地表最大径向位移ur所在的位置越偏离岩浆囊中心.无地形以及当坡度较小时(图 4a-b)地表最大垂向位移uz所在的位置位于岩浆囊正上方，但当坡度达到30°时(图 4c)，最大垂向位移uz所在位置不再位于岩浆囊正上方.因此，当坡度超过30°时利用Mogi模型可能会给出错误的岩浆房位置.
3 长白山火山区岩浆房计算 3.1 岩浆房源参数反演长白山天池火山区已建成13个水准点和8个GPS站点，且均开展多期测量，地表形变(崔笃信等，2007；Ji et al., 2013)、地震活动性(吴建平等, 2005, 2007)和火山气体(Xu et al., 2012)等方面资料显示2002-2003年该区具有典型岩浆房活动.我们选用这一时间段该区域扣除东北地块整体运动背景场后的岩浆房引起的GPS水平位移和水准垂向数据(陈国浒等，2008)作为目标约束.值得注意的是，采用水准测量作为反演约束时，由于水准测量给出的是其他测点相对S0基准点的相对垂向距离，距离长白山火山口最远的S0测点的绝对垂向位移是一个较小的未知值，但并不为0.前人工作(陈国浒等，2008)往往假设S0测点垂向位移绝对值为0会引起微小误差，为此本文参考朱桂芝等(2008)方法，在反演过程中不做S0测点垂向位移为0的简化.
计算网格经过3次迭代自适应加密，最终的六面体单元总数为381529，网格包含地表起伏，见图 5.边界条件施加如下，地表自由，侧面和底边界无法向位移，岩浆房处施加等效体力.根据PREM模型，模型的介质参数为νp=6.7 km·s-1，νs=3.87 km·s-1, ρ=2900 kg·m-3.按照2.3节模拟反演方案对岩浆房压力变形源参数进行反演.
|
图 5 长白山火山区岩浆压力变形源有限元网格 Fig. 5 FEM grid of magma pressure sources beneath Changbaishan volcano area |
算例1，只考虑一个单独的椭球状岩浆房，反演参数见表 1，平均拟合残差为0.87 mm，但P3点观测值明显偏离火山口.朱桂芝等(2008)在反演中认为P3点的水平运动趋势相对其他测点明显偏离火山口中心，并把该点数据剔除.前人采用单独岩浆房的反演均难以拟合该点的水平运动(胡亚轩等, 2007; 朱桂芝等, 2008)，算例1只采用单独岩浆房也难以拟合P3点水平运动.对于水准测量垂直位移的拟合，整体趋势吻合，但水准测量的S2点比S3点距离火山口更远，垂直位移却更大，不排除S2点受局部张拉断层/岩浆侵入的影响.从地表形变的拟合情况看，椭球状压力膨胀源基本上能一阶近似描述长白山火山区地下的岩浆房.在反演中，压力差ΔP和岩浆房大小有关，根据式(1)和式(2)，对于Mogi模型，岩浆房等效体积变化ΔV=ΔPa3，唯一的地表形变对应唯一的岩浆房等效体积变化ΔV，但压力差ΔP和岩浆房半径a并不唯一.对于椭球岩浆房，压力差ΔP和岩浆房三主轴的确定以能够最好拟合地表形变为原则.前人的模拟认为压力差ΔP一般为10~30 MPa量级(Lisowski, 2007; Trasatti et al., 2010)，经计算，当压力差ΔP和岩浆房的三主轴如表 1所示时地表形变拟合最佳，同时岩浆房位置与电性结构等结果反映的岩浆房位置大体一致.
|
表 1 反演结果 Table 1 Inversion results |
对于P3点的位移异常，陈国浒等(2008)结合反演和间白山地震活动性，认为在间白山西侧即P3点附近存在另一个岩浆房压力变形源，并且推测存在隐伏的北西西向断裂连接长白山天池岩浆房和间白山岩浆房.据此，在算例2模拟中，在P3的附近东南方增加了一个球状岩浆房压力变形源(球岩浆房1)，并且认为该岩浆房与主椭球岩浆房是连通的，因此二者的压力取相同值.椭球岩浆房加上球岩浆房1得到的拟合结果见图 6，此时P3处的拟合明显改善，算例2的平均拟合残差为0.68 mm.
|
图 6 模拟结果和地表形变观测的比较 (a)-(c)分别为算例1~3模拟结果和GPS以及水准测量的比较；P0~P7为GPS测点，S0~S12为水准测点. Fig. 6 Comparison of displacements from modeling and observations Deformation modeling results from (a) case 1, (b) case 2 and (c) case 3. GPS sites P0~P7 and leveling sites S0~S12 are shown in Fig. 6. |
算例1、2均未能对S2点的垂直位移突变作出解释，因为垂向位移突变分布的范围较小，因此猜测S2点地下的干扰源应该较浅且较小.干扰源并不一定是岩浆房，也有可能是岩浆房压力作用下产生的微小张拉断层，但有限元处理张拉断层比较困难，此处算例3暂用一个半径较小的球源2来模拟S2点地下的干扰源，球源2的详细参数见表 1，算例3的拟合残差为0.65 mm.
最终以算例3确定的岩浆房作为最优模型，见图 7，椭球岩浆房与2002-2003年环长白山天池小震群空间分布(吴建平等, 2005, 2007)较为吻合，岩浆房深度范围为4~10 km之间，这和长白山火山区上地壳三维速度层析成像、电性结构等结果反映的岩浆房深度大体一致.
|
图 7 最优模型 (a)顶视图；(b)南视图；(c)东视图. Fig. 7 Optimal model of magama chamber (a) Top view; (b) Southward view; (c) Eastward view. |
西太平洋板块和欧亚板块的碰撞作用控制下，长白山火山区现今以东西向挤压为主，该区最大主压应力为NEE-SWW方向(李春锋等，2006).张春山等(2016)在长白山天池周边开展了水压致裂地应力测量，发现现今最大水平主应力方向在天池西部大青川测点(DQC)为近EW向为主, 与区域构造应力方向接近.但天池北部东清测点(DQ)和冰胡屯测点(BHT)以NW-NNW向为主与区域构造方向近乎正交.为了解释这种现象，张春山等(2016)认为长白山天池附近的应力环境复杂，是区域构造和岩浆活动共同作用的结果.于是提出了问题，即岩浆房扰动应力场Δσ和西太平洋板块俯冲主导的近NEE-SWW向区域构造应力场σ0的叠加能否解释长白山天池周边的应力测量结果？为了回答这个问题，我们给出了5 km深度的Δσxx、Δσyy、Δσxy云图(图 8a-c)，三分量均表现出明显花瓣特征，图 8d是根据有限元计算的扰动应力的最大主应力方向通过平滑方法(Hansen and Mount, 1990)得到的岩浆房引起的均匀网格点应力方向，整体趋势是天池西侧的Δσ主方向以近EW向为主，天池北部Δσ主方向以NNE-近NS向为主，因此岩浆房扰动应力场Δσ和区域构造应力场σ0的叠加有可能造成天池西部近EW向，天池北部以NW-NNW向为主的现今应力方向.由于很难获取准确的构造应力场量值，本文未定量计算构造应力场和扰动应力场的叠加，除了构造应力和岩浆房扰动，地形地貌、断裂等因素也会影响局部应力，张春山等(2016)的地应力测点均选在地形较平缓丘陵区，因此地貌对地应力方向的影响较小，天池西侧和北侧地应力方向差异更细致的原因还需要对区域构造应力、断裂作用等因素有更深入了解后才能得到真正的定量解释.可以确定的是距离岩浆房较近部位更容易受应力扰动的影响，很可能影响近处的断层稳定性或者地震活动.
|
图 8 岩浆房引起的扰动应力场 (a) Δσxx；(b) Δσyy；(c) Δσxy；(d)水平最大主应力方向. DQC:大青川测点，DQ:东清测点，BHT:冰胡屯测点. Fig. 8 Magma-induced local stress field (a) E-W normal horizontal stress; (b) N-S normal horizontal stress; (c) Horizontal shear stress; (d) Orientation of maximum horizontal stress. |
根据流动地震观测网记录(吴建平等，2007)，长白山天池火山地震自2002年7月显著增多, 震中环天池火山口分布，距天池水面深度小于5 km，且震中多发生在天池东北和西南部.震群高精度定位揭示，震群断层面走向NW，倾向西南，倾角约80°.岩浆房压力扰动是否有利于该节面的地震发生？我们通过库仑应力变化的计算来回答上述问题，根据张春山等(2016)本区现今NW向断裂具有左旋走滑的性质，因此库仑应力变化的发震节面选择为走向300°，倾角80°，滑动角0°.
我们采用库仑应力变化(ΔCFS)研究岩浆房扰动对断层稳定性影响，通常ΔCFS(Harris, 1998; 石耀霖和曹建玲, 2010)定义为
|
(9) |
式中，Δτ是断层面剪应力变化，Δσn是断层面上正应力变化，ΔP是孔隙压变化，μ是断层摩擦系数.当库仑应力为正，扰动应力有利于断层滑动, 反之则不利于断层滑动.
当岩石应力变化远快于岩石中流体压力扩散时，孔隙压变化ΔP可用Skemptons系数B表示，实际应用中常将视摩擦系数简化为断层摩擦特性和孔隙流体的综合参数，并取视摩擦系数μ′=μ(1-B).本文采用这种简化，并讨论不同视摩擦系数下的库仑应力变化：
|
(10) |
图 9a-b分别给出了岩浆房压力膨胀导致的断层接收面上的剪应力变化和正应力变化，图 9b-d，分别对应了视摩擦系数μ′为0、0.4、0.8的情形，随着视摩擦系数的增大，ΔCFS正值区有顺时针旋转的趋势.研究表明视摩擦系数的选取与断层类型(Parsons et al., 1999; Ali et al., 2008)以及断层滑动速率相关(Parsons et al., 1999).对于视摩擦系数的选择，通常滑动速率较大的走滑断层、正断层取0.2~0.4，速率较小的逆断层取0.6~0.8(Parsons et al., 1999).由于区域内NW向断层主要表现出现左旋走滑和张扭性特征，因此视摩擦系数取为0.4较为合理，μ′=0.4取值下的库仑应力变化能够解释天池火山口NW向震群在空间上主要分布于火山口的西南和东北部.类似地火山活跃期伴随地震事件的现象，在Vargas-Bracamontes和Neuberg(2012)发现火山喷发前往往伴随着VT(Volcano-Tectonic)地震，并且地震节面相对于构造应力场方向旋转近90°.
|
图 9 岩浆房引起的库仑应力变化 (a) Δσ；(b) μ′为0时的ΔCFS；(c) μ′为0.4时的ΔCFS；(d) μ′为0.8时的ΔCFS. Fig. 9 Magma-induced Coulomb failure stress changes (a) Δσ; (b)-(d) ΔCFS with friction coefficient 0, 0.4 and 0.8. |
Cayol和Cornet(1998)利用含地形有限元计算地表位移场并以此为目标约束利用Mogi模型反演地下岩浆房，以此研究地形对反演的影响，他们的结果显示Mogi模型反演的海平面以下岩浆房深度约等于有限元模型海平面以下岩浆房深度加上地表高程，即当把Mogi模型反演深度视为地表峰顶至岩浆房的垂向距离，Mogi模型可以给出近乎准确的岩浆房深度，但Mogi模型高估了岩浆房体积变化量, 例如地形坡度为20°时，岩浆房体积变化量可能被高估15%~20%，当地形坡度为30°时，岩浆房体积变化量可能被高估50%.朱桂芝等(2008)未考虑地表地形的反演结果显示，岩浆房深度为9.2 km，我们考虑地形反演的海平面以下岩浆房深度为6.9 km，考虑到长白山主峰约为2.4 km，根据Cayol和Cornet(1998)结论，二者反演的岩浆房深度并不矛盾，即岩浆房距离长白山主峰峰顶的深度约为9.3 km.陈国浒等(2008)利用单源和双源Mogi模型反演岩浆房位置，单源结果显示岩浆房体积变化为5.2×106m3，朱桂芝反演的等量三正交断层扩张量为2.9×106m3.二者的扩张量不一样，这是由于二者采用模型不一样，根据Feigl等(2000)，Okada位错理论的3正交等量扩张位错模型和Mogi模型引起相同地表位移时，二者的扩张量具有(1+ν)/(1-ν)的比例关系, 其中ν为泊松比.值得注意的是我们采用跟朱桂芝等(2008)相同的等量三正交断层扩张量约为2.5×106m3，这也反映了解析方法虽然可以大致限定深度，却会高估断层扩展量.由于反演的多解性，不同作者采用的方法不一样，有些反演方法依赖初值，因此造成了不同作者反演结果的差别，但不同作者大体上的岩浆房参数基本上与长白山火山区上地壳三维速度层析成像、电性结构等观测所刻画的岩浆房特征相吻合.同时，值得注意的是，陈国浒等(2008)的反演的双岩浆房分别位于长白山天池火山和间白山火山下方，本文的计算也支持在间白山底下存在和主岩浆房连通的小岩浆房时反演效果更佳，胡亚轩等(2004)的计算认为当考虑天池南侧的北西向隐伏断裂时反演效果更佳，三者的结论较为接近，即长白山天池岩浆房和间白山岩浆房之间可能存在隐伏断裂连接，这也证明了不同作者反演模型的一致性.
4.2 构造应力场的影响Vargas-Bracamontes和Neuberg(2012)在研究VT(Volcano-Tectonic)地震的时候，发现由于优势区域应力场的存在，火山区的断层滑动往往具有滑动方向和构造应力场方向一致或者接近的特点.并且他们的研究还考虑了岩浆房压力和构造应力场的不同相对大小时，构造应力场对库仑应力变化的影响.根据张春山等(2016)的结果，浅表水平主应力为2.31~12.39 MPa, 深部应力难以获知，我们参考Vargas-Bracamontes和Neuberg(2012)的方法，考虑构造应力偏应力σ1=0.5Δp和σ1=0.1Δp两种情况，分别计算库仑应力变化，见图 10.随着构造应力量值的增加，库仑应力变化的值增大，这是因为现今的构造应力场有利于长白山地区先存的NW向断裂活动.尽管深部构造应力场具体量值很难获得，很难获得定量结论，但我们仍然可以获得扰动应力场、构造应力场均有利于火山口附近NW向断裂的发生的定性结论.
|
图 10 构造应力场对库仑应力变化的影响 Fig. 10 Coulomb failure stress changes due to the interaction of regional maximum compression σ1 parallel to x-axis and internal pressurization |
值得注意的是，反映火山内部活动的还应有火山岩浆房的地震、重力、地磁探测等多种方法.本文的方法仅仅通过形变的方法来探索火山内部岩浆活动的趋势，对于缺乏形变数据的火山地区该方法有较大局限性.由于有限元方法的通用性，未来还可以求解岩浆房内部活动引起的重力变化(Cai and Wang, 2005; Battaglia and Hill, 2009)，也应该进一步考虑岩浆房由于温度差异引起的热应力(Cattin et al, 2005)和孔隙流体作用对岩浆房稳定性的影响(Gerbault et al., 2012).
5 结论通过“有限元等效体力”方法考虑地表地形影响，开展长白山火山区的地下岩浆房源参数反演，并计算岩浆应力扰动对周边断层稳定性的影响，得到以下结论：
(1) 在火山区地下压力变形源的地表形变计算中，地表地形影响难以忽略.岩浆房埋深越浅，地表位移受地形影响越大.埋深越浅，地表最大径向位移ur所在的位置越靠近岩浆房.当坡度达30°，最大垂向位移uz所在位置不再位于岩浆房正上方.
(2) 椭球状岩浆房压力源可较好模拟长白山火山地区2002-2003年间的地表形变测量.浅部可能存在局部张拉断层或岩墙侵入，可用较小的球状压力源近似拟合.
(3) 岩浆房扰动应力场和区域构造应力场的叠加有可能造成天池西部近EW向，天池北部以NW-NNW向为主的现今应力方向.
(4) 岩浆房压力引起的库仑应力变化能够解释天池火山口NW向震群在空间上主要分布于火山口的西南和东北部.
致谢 感谢绘图工具GMT(Wessel and Smith, 1995).
Aki K, Richards P G. 2002. Quantitative Seismology. 2nd ed. California: University Science Books.
|
Ali S T, Freed A M, Calais E, et al. 2008. Coulomb stress evolution in Northeastern Caribbean over the past 250 years due to coseismic, postseismic and interseismic deformation. Geophysical Journal International, 174(3): 904-918. DOI:10.1111/j.1365-246X.2008.03634.x |
Bangerth W, Hartmann R, Kanschat G. 2007. deal. Ⅱ-a general-purpose object-oriented finite element library. ACM Transactions on Mathematical Software, 33(4): 24. DOI:10.1145/1268776.1268779 |
Battaglia M, Hill D P. 2009. Analytical modeling of gravity changes and crustal deformation at volcanoes:The Long Valley caldera, California, case study. Tectonophysics, 471(1-2): 45-57. DOI:10.1016/j.tecto.2008.09.040 |
Bhattacharyya J, Sheehan A F, Tiampo K F, et al. 1999. Using genetic algorithm to model Broadband regional waveforms for crustal structure in the western United States. Bulletin of the Seismological Society of America, 89(1): 202-214. |
Billings S D, Kennett B L N, Sambridge M S. 1994. Hypocentre location:genetic algorithms incorporating problem-specific information. Geophysical Journal International, 118(3): 693-706. DOI:10.1111/j.1365-246X.1994.tb03994.x |
Bonaccorso A, Davis P M. 1999. Models of ground deformation from vertical volcanic conduits with application to eruptions of Mount St. Helens and Mount Etn. Journal of Geophysical Research:Solid Earth, 104(B5): 10531-10542. DOI:10.1029/1999JB900054 |
Bonafede M, Ferrari C. 2009. Analytical models of deformation and residual gravity changes due to a Mogi source in a viscoelastic medium. Tectonophysics, 471(1-2): 4-13. DOI:10.1016/j.tecto.2008.10.006 |
Burridge R, Knopoff L. 1964. Body force equivalents for seismic dislocations. Bulletin of the Seismological Society of America, 54(6A): 1875-1888. |
Cai Y G, Wang C Y. 2005. Fast finite-element calculation of gravity anomaly in complex geological regions. Geophysical Journal International, 162(3): 696-708. DOI:10.1111/j.1365-246X.2005.02711.x |
Cattin R, Doubre C, De Chabalier J B, et al. 2005. Numerical modelling of quaternary deformation and post-rifting displacement in the Asal-Ghoubbet rift (Djibouti, Africa). Earth and Planetary Science Letters, 239(3-4): 352-367. DOI:10.1016/j.epsl.2005.07.028 |
Cayol V, Cornet F H. 1998. Effects of topography on the interpretation of the deformation field of prominent volcanoes-application to Etna. Geophysical Research Letters, 25(11): 1979-1982. DOI:10.1029/98GL51512 |
Charco M, Del Sastre P G. 2014. Efficient inversion of three-dimensional finite element models of volcano deformation. Geophysical Journal International, 196(3): 1441-1454. DOI:10.1093/gji/ggt490 |
Chen G H, Shan X J, Moon W M, et al. 2008. A modeling of the magma chamber beneath the Changbai Mountains lcanic area constrained by InSAR and GPS derived deformation. Chinese Journal of Geophysics (in Chinese), 51(4): 1085-1092. |
Cui D X, Wang Q L, Li K, et al. 2007. Analysis of recent deformation of Changbaishan Tianchi volcano. Chinese Journal of Geophysics (in Chinese), 50(6): 1731-1739. |
Davis P M. 1986. Surface deformation due to inflation of an arbitrarily oriented triaxial ellipsoidal cavity in an elastic half-space, with reference to Kilauea volcano, Hawaii. Journal of Geophysical Research:Solid Earth, 91(B7): 7429-7438. DOI:10.1029/JB091iB07p07429 |
Feigl K L, Gasperi J, Sigmundsson F, et al. 2000. Crustal deformation near Hengill volcano, Iceland 1993-1998:coupling between magmatic activity and faulting inferred from elastic modeling of satellite radar interferograms. Journal of Geophysical Research:Solid Earth, 105(B11): 25655-25670. DOI:10.1029/2000JB900209 |
Fernández J, Tiampo K F, Jentzsch G, et al. 2001. Inflation or deflation? New results for Mayon Volcano applying elastic-gravitational modeling. Geophysical Research Letters, 28(12): 2349-2352. DOI:10.1029/2000GL012656 |
Fialko Y, Khazan Y, Simons M. 2010. Deformation due to a pressurized horizontal circular crack in an elastic half-space, with applications to volcano geodesy. Geophysical Journal of the Royal Astronomical Society, 146(1): 181-190. |
Gerbault M, Cappa F, Hassani R. 2012. Elasto-plastic and hydromechanical models of failure around an infinitely long magma chamber. Geochemistry, Geophysics, Geosystems, 13(3). DOI:10.1029/2011GC003917 |
Hansen K M, Mount V S. 1990. Smoothing and extrapolation of crustal stress orientation measurements. Journal of Geophysical Research:Solid Earth, 95(B2): 1155-1165. DOI:10.1029/JB095iB02p01155 |
Harris R A. 1998. Introduction to special section:Stress triggers, stress shadows, and implications for seismic hazard. Journal of Geophysical Research:Solid Earth, 103(B10): 24347-24358. DOI:10.1029/98JB01576 |
Hong H, Kadlec B J, Yuen D A, et al. 2004. Fast timescale phenomena at Changbaishan volcano as inferred from recent seismic activities.//American Geophysical Union, Fall Meeting 2004. AGU.
|
Hu Y X, Wang Q L, Cui D X, et al. 2004. Joint inversion of geometric deformation in Changbaishan volcanic area. Journal of Geodesy and Geodynamics (in Chinese), 24(4): 90-94. |
Hu Y X, Wang Q L, Cui D X, et al. 2007. Application of Mogi model at Changbaishan Tianchi volcano. Seismology and Geology (in Chinese), 29(1): 144-151. |
Ji L Y, Xu J D, Wang Q L, et al. 2013. Episodic deformation at Changbaishan Tianchi volcano, northeast China during 2004 to 2010, observed by persistent scatterer interferometric synthetic aperture radar. Journal of Applied Remote Sensing, 7(1): 073499. DOI:10.1117/1.JRS.7.073499 |
Kelly D W, Gago J P D S R, Zienkiewicz O C, et al. 1983. A posteriori error analysis and adaptive processes in the finite element method:Part I-error analysis. International Journal for Numerical Methods in Engineering, 19(11): 1593-1619. DOI:10.1002/nme.1620191103 |
Li C F, Zhang X K, Zhang Y, et al. 2006. Analysis of tectonic setting of Changbaishan Tianchi Volcano. Seismological and Geomagnetic Observation and Research (in Chinese), 27(5): 43-49. |
Li H Q, Zeng Z F, Liu S X, et al. 2017. Study on coupling response of multiphysical field around Tianchi lake of Changbaishan volcano based on COMSOL Multiphysics. Progress in Geophysics (in Chinese), 32(4): 1779-1783. DOI:10.6038/pg20170449 |
Lisowski M. 2007. Analytical volcano deformation source models.//Dzurisin D ed. Volcano Deformation. Berlin Heidelberg: Springer, 279-304.
|
Ma X Y, Teng J W, Liu Y S, et al. 2016. Physical property structure of the crust-mantle and deep geophysical feature in Changbaishan volcanic area. Progress in Geophysics (in Chinese), 31(5): 1973-1985. DOI:10.6038/pg20160512 |
Mindlin R D. 1936. Force at a point in the interior of a semi-infinite solid. Physics, 7(5): 195-202. DOI:10.1063/1.1745385 |
Mogi K. 1958. Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. Bulletin of the Earthquake Research Institute University of Tokyo, 36(2): 99-134. |
Okada Y. 1985. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America, 75(5): 1135-1154. |
Okada Y. 1992. Internal deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America, 82(2): 1018-1040. |
Parsons T, Stein R S, Simpson R W, et al. 1999. Stress sensitivity of fault seismicity:a comparison between limited-offset oblique and major strike-slip faults. Journal of Geophysical Research:Solid Earth, 104(B9): 20183-20202. DOI:10.1029/1999JB900056 |
Shi Y L, Assumpcao M. 2000. Genetic algorithm-finite element inversion of stress field of Brazil. Chinese Journal of Geophysics (in Chinese), 43(2): 166-174. |
Shi Y L, Cao J L. 2010. Some aspects in static stress change calculation-Case study on Wenchuan earthquake. Chinese Journal of Geophysics (in Chinese), 53(1): 102-110. DOI:10.3969/j.issn.0001-5733.2010.01.011 |
Tang J, Deng Q H, Zhao G Z, et al. 2001. Electric conductivity and magma chamber at the Tianchi volcano area in Changbaishan Mountain. Seismology and Geology (in Chinese), 23(2): 191-200. |
Tiampo K F, Fernández J, Jentzsch G, et al. 2004. Volcanic source inversion using a genetic algorithm and an elastic-gravitational layered earth model for magmatic intrusions. Computers & Geosciences, 30(9-10): 985-1001. |
Trasatti E, Giunchi C, Agostinetti N P. 2010. Numerical inversion of deformation caused by pressure sources:Application to Mount Etna (Italy). Geophysical Journal International, 172(2): 873-884. |
Vargas-Bracamontes D M, Neuberg J W. 2012. Interaction between regional and magma-induced stresses and their impact on volcano-tectonic seismicity. Journal of Volcanology and Geothermal Research, 243-244: 91-96. DOI:10.1016/j.jvolgeores.2012.06.025 |
Wei H Q, Liu R X, Li X D. 1997. Ignimbrite-forming eruptions from Tianchi volcano and their climate effect. Earth Science Frontiers (in Chinese), 4(1-2): 263-266. |
Wessel P, Smith W H F. 1995. New version of the generic mapping tools. Eos, Transactions American Geophysical Union, 76(33): 329. |
Williams C A, Wadge G. 1998. The effects of topography on magma chamber deformation models:Application to Mt. Etna and radar interferometry. Geophysical Research Letters, 25(10): 1549-1552. DOI:10.1029/98GL01136 |
Williams C A, Wadge G. 2000. An accurate and efficient method for including the effects of topography in three-dimensional elastic models of ground deformation with applications to radar interferometry. Journal of Geophysical Research:Solid Earth, 105(B4): 8103-8120. DOI:10.1029/1999JB900307 |
Wu J P, Ming Y H, Zhang H R. 2005. Seismic activity at the Changbaishan Tianchi volcano in the summer of 2002. Chinese Journal of Geophysics (in Chinese), 48(3): 621-628. |
Wu J P, Ming Y H, Zhang H R, et al. 2007. Earthquake swarm activity in Changbaishan Tianchi volcano. Chinese Journal of Geophysics (in Chinese), 50(4): 1089-1096. DOI:10.1002/cjg2.1126 |
Xu J D, Liu G M, Wu J P, et al. 2012. Recent unrest of Changbaishan volcano, northeast China:A precursor of a future eruption?. Geophysical Research Letters, 39(16): L16305. DOI:10.1029/2012GL052600 |
Yang X M, Davis P M, Dieterich J H. 1988. Deformation from inflation of a dipping finite prolate spheroid in an elastic half-space as a model for volcanic stressing. Journal of Geophysical Research:Solid Earth, 93(B5): 4249-4257. DOI:10.1029/JB093iB05p04249 |
Zhang C K, Zhang X K, Zhao J R, et al. 2002. Study on the crustal and upper mantle structure in the Tianchi volcanic region and its adjacent area of Changbaishan. Chinese Journal of Geophysics (in Chinese), 45(6): 812-820. |
Zhang C S, Wu M L, Zhang C Y, et al. 2016. Measurement of present-day stress and analysis of stress state in the Changbaishan mountains of Jilin province. Chinese Journal of Geophysics (in Chinese), 59(3): 922-930. DOI:10.6038/cjg20160314 |
Zhao D P, Lei J S, Tang R Y. 2004. Origin of the Changbai intraplate volcanism in Northeast China:evidence from seismic tomography. Chinese Science Bulletin, 49(13): 1401-1408. DOI:10.1360/04wd0125 |
Zhu G Z, Wang Q L, Shi Y L, et al. 2008. Modelling pressurized deformation source for Changbaishan volcano with homogenous expansion point source. Chinese Journal of Geophysics (in Chinese), 51(1): 108-115. |
Zhu S B. 2005. Improvement of genetic algorithm and its application on stress field inversion. Northwestern Seismological Journal (in Chinese), 27(2): 97-100. |
陈国浒, 单新建, Moon W M, 等. 2008. 基于InSAR、GPS形变场的长白山地区火山岩浆囊参数模拟研究. 地球物理学报, 51(4): 1085-1092. |
崔笃信, 王庆良, 李克, 等. 2007. 长白山天池火山近期形变场演化过程分析. 地球物理学报, 50(6): 1731-1739. |
胡亚轩, 王庆良, 崔笃信, 等. 2004. 长白山火山区几何形变的联合反演. 大地测量与地球动力学, 24(4): 90-94. |
胡亚轩, 王庆良, 崔笃信, 等. 2007. Mogi模型在长白山天池火山区的应用. 地震地质, 29(1): 144-151. |
李春锋, 张兴科, 张旸, 等. 2006. 长白山天池火山的地质构造背景. 地震地磁观测与研究, 27(5): 43-49. |
李宏卿, 曾昭发, 刘四新, 等. 2017. 基于COMSOL Multiphysics的长白山天池多物理场耦合响应研究. 地球物理学进展, 32(4): 1779-1783. DOI:10.6038/pg20170449 |
马学英, 滕吉文, 刘有山, 等. 2016. 长白山火山区壳幔物性结构与深部物理场特征. 地球物理学进展, 31(5): 1973-1985. DOI:10.6038/pg20160512 |
石耀霖, Assumpcao M. 2000. 巴西构造应力场的遗传有限单元法反演. 地球物理学报, 43(2): 166-174. |
石耀霖, 曹建玲. 2010. 库仑应力计算及应用过程中若干问题的讨论-以汶川地震为例. 地球物理学报, 53(1): 102-110. DOI:10.3969/j.issn.0001-5733.2010.01.011 |
汤吉, 邓前辉, 赵国泽, 等. 2001. 长白山天池火山区电性结构和岩浆系统. 地震地质, 23(2): 191-200. |
魏海泉, 刘若新, 李晓东. 1997. 长白山天池火山造伊格尼姆岩喷发及气候效应. 地学前缘, 4(1-2): 263-266. |
吴建平, 明跃红, 张恒荣, 等. 2005. 2002年夏季长白山天池火山区的地震活动研究. 地球物理学报, 48(3): 621-628. |
吴建平, 明跃红, 张恒荣, 等. 2007. 长白山天池火山区的震群活动研究. 地球物理学报, 50(4): 1089-1096. |
张成科, 张先康, 赵金仁, 等. 2002. 长白山天池火山区及邻近地区壳幔结构探测研究. 地球物理学报, 45(6): 812-820. |
张春山, 吴满路, 张重远, 等. 2016. 长白山地区现今地应力测量结果与应力状态分析. 地球物理学报, 59(3): 922-930. DOI:10.6038/cjg20160314 |
赵大鹏, 雷建设, 唐荣余. 2004. 中国东北长白山火山的起源:地震层析成像证据. 科学通报, 49(14): 1439-1446. |
朱桂芝, 王庆良, 石耀霖, 等. 2008. 各向同性膨胀点源模拟长白山火山区岩浆囊压力变形源. 地球物理学报, 51(1): 108-115. |
朱守彪. 2005. 遗传算法的改进及其在应力场反演中的应用. 西北地震学报, 27(2): 97-100. |