中国东北地区在地质构造上位于西伯利亚克拉通与华北克拉通之间，属古亚洲构造域、蒙古—鄂霍茨克构造域和滨太平洋构造域叠合作用的部位，是大陆块体作用较为活跃的区域，经历多期构造作用，固体矿产和油气资源丰富，长期以来是地学研究的热点区域(杨宝俊等，2006；Kelty et al., 2008；Zhou，2009).涂广红等(2006)利用重磁场资料对东北地区地质构造进行分析，总结出剩余重磁异常特征与该区地质构造有良好的对应关系.杜晓娟等(2009)根据该区域重力场的异常特征，将东北地区分为兴安、松嫩及张广才岭三大异常区，并划分出与重力异常分区对应的一、二级构造.众多的天然地震成像结果获得的速度结构表明东北地区的地幔转换带(410~660 km)内分布着大规模的高速异常体，并且进一步讨论了区域构造演化和东北地区火山起源问题，说明其与该地区大地幔楔内软流圈上涌和滞留的西太平洋板片的关系(赵大鹏等，2004；田有等，2011；徐峣等，2014；王开燕等，2015).索奎等(2015)通过重震反演的方法得到了东北地区地壳上地幔0~500 km的密度分布，研究了其与断裂的关系，以及火山的深部起源.

天然地震层析成像技术是一种研究地球内部速度结构的有效方法，其提供的三维P波速度结构能够较好地反映深大断裂位置、确定深部岩浆活动区域.地幔中不均匀的物质质量分布会形成导致地幔对流、地幔柱和板块俯冲等动力学过程的应力，因此了解上地幔的密度结构是理解地球内部构造运动和动力学机制的关键(Kaban et al., 2016).重力反演是根据观测异常对地下密度分布的定量计算，得到的三维密度结构能够反映地下密度分布状态.Braitenberg等(2006)年利用重力反演研究南海海岭的地壳结构和岩石圈的弹性厚度，研究中地幔结构的变化没有被考虑.Bai等(2013)利用喜马拉雅山和西藏的重力数据构建了一个详细的地壳模型，暗示了一个均匀的上地幔结构.有些研究假设地震速度和密度之间有直接的线性关系，尝试拟合重力和地震两个数据集(Maceira and Ammon, 2009；Tiberi et al., 2008)，但这些方法的一个明显的缺点在于没有考虑到速度和密度可能会以不同的方式取决于温度条件和物质组成，速度强烈的取决于温度条件(Goes and Lee, 2002)，而密度可以同时受两种因素共同影响(Kaban et al., 2015).目前对中国东北地区地下深部密度结构研究还依旧较少，对该地区的深部三维密度结构的研究有助于揭示其构造演化(如岩浆活动、火山活动、岩石圈演化等)与西太平洋板片俯冲的关系，以及深部物质与能量交换和动力学过程，对该地区矿产、油气资源的突破具有很大的指导意义.

本次研究区域(34—55°N，99—146°E)主要包括中国东北地区、华北部分地区、蒙古部分地域及日本海海域.中国东北地区处于中朝克拉通与西伯利亚克拉通过渡的复杂构造演化带，属古亚洲构造域和环太平洋构造域构造叠加作用最为显著的地区，地质构造错综复杂，主要地质构造单元走向呈北东向和北北东向展布.区域内分布有大兴安岭和张广才岭两条巨大岩浆岩带和镜泊湖、五大连池、长白山等新生代火山以及多个盆地(松辽盆地等)(潘佳铁等，2014)，油气和固体矿产资源比较丰富.此外，欧亚大陆下西太平洋板片的俯冲角度和速度随时间在变化(Koppers et al., 2001, 2003；Sharp and Clague, 2006)，其俯冲板片是否穿过地幔转换带(410~660 km)一直是研究的热点(Liu et al., 2015；Wei et al., 2012；Zhang et al., 2016)，获取该区的深部的密度结构，对于研究板块俯冲、火山活动等动力学过程具有重要的科学意义.

本文利用卫星重力异常，去除沉积层和莫霍面等产生的重力场得到剩余布格重力异常数据，反演得到整个区域的上地幔三维密度结构，刻画地下构造的几何形态，同时结合天然地震层析成像得到的三维速度模型，共同揭示西太平洋板片在研究区域的俯冲形态，为东北地区深部构造演化和动力学过程的研究提供依据.

根据1:200万卫星重力数据(图 1a)，利用数据径向对数功率谱斜率与深度的对应关系去除50 km以上沉积层和莫霍面起伏所引起的异常(Blakely，1995；穆石敏等，1990；贺日政等，2007)，参考深度依据前人研究(Chen，2010；Guo et al., 2015)，利用剩余的区域重力异常(图 1b)反演地下的密度不均匀分布.

图 1

Fig. 1

图 1 研究区卫星重力异常 (a)原始重力异常；(b)区域重力异常.黑色曲线代表大兴安岭断裂带，1到16号线为反演剖面位置. Fig. 1 Satellite gravity anomaly of study area (a) Original gravity anomaly; (b) Regional gravity anomaly; Black curve represents the Da Hinggan Mountains fault zone, and lines 1 to 16 are inversion profiles.

三维重力反演是将目标区域地下划分为有限多的密度一定的棱柱体单元.为了与地震解释相对应，反演的网格大小水平方向为1°×1°，垂向网格剖分间距为50 km，这样将研究区地下划分为48×22×20个棱柱体.地面上任意观测点的重力值是地下各棱柱体对此产生重力影响的总和，设观测点坐标为O=(x₀，y₀，z₀)，密度差为ρ的单一棱柱体对观测点产生的重力异常为(Li and Chouteau, 1998)：

(1)

式中：x_i=x₀－ξ_i，y_j=y₀－η_j，z_k=z₀－ζ_k，i，j，k=1，2；ξ，η，ζ定义长方体单元的八个角点；r_ijk= ，u_ijk=(－1)ⁱ(－1)^j(－1)^k；γ为万有引力常数，γ=6.67×10^－11m³/(kg·s²).因此，由密度分布产生的重力异常可以近似为所有长方体单元响应的叠加，如果有n个观测点和m个棱柱体单元，其线性关系表达式为

(2)

其中，g是地表观测重力异常值，ρ是地下空间密度分布，G是重力正演核函数.

由于重力反演问题是不适定问题，地下网格个数远大于观测数点个数，反演方程是病态的，另外，反演需要大量计算资源，计算时间也较长.据此，本文选取正则化方法约束Tikhnov参数方程，通过约束项平衡数据拟合差与模型目标函数，减小反演多解性，建立目标函数为

(3)

s(ρ)是稳定函数，φ(ρ)是观测数据与预测数据的拟合函数，α为正则化参数，其值对反演结果有很大影响，本文选取自适应方法(Zhdanov，2002)选择正则化参数，即随着反演结果的收敛逐渐减小稳定函数的贡献值，

(4)

(5)

深度加权函数的加入可有效解决反演密度集中于地表的“趋肤效应”，增强反演结果在深度方向的分辨率.我们取Li和Oldenburg(2000)采用的一种基于灵敏度矩阵的加权函数，其优点是仅仅用到核函数矩阵的值，不需要先验信息及避免了其他参数的影响，表达式为

(6)

式中，G_ij表示正演核函数中的元素，i为观测点编号，j为地下网格编号，β是常数，一般取1，调整其值，来控制加权函数的强度.将深度加权函数引入正则化方程后，式(3)可以转化为

(7)

此方法能大大增加反演结果的光滑性，并一定程度上克服重力反演深度分辨率不好的缺点.由于观测数据量大，离散的网格数比较多，反演需要耗费大量的计算资源和时间，本文选取预处理共轭梯度算法快速求解大型线性方程组(Pilkington，1997；Wang et al., 2017)，其只需要进行矩阵向量乘积，需要低存储和低计算成本，能够大大提高反演的效率.

本次重力反演得到了研究区域0到1000 km深度范围的三维密度结构，截取了如图 1所示剖面位置的密度反演结果(图 2和图 3).综合图 6的三维速度结构结果，发现，在浅部(100~300 km处)大兴安岭造山带区域主要表现低密度低波速异常，松辽盆地表现较弱的高密高速异常并伴随着部分低密度低波速异常区域；在深部(300~700 km处)松辽盆地和大兴安岭造山带附近仍以高密高速异常为主.欧亚大陆下的太平洋板块俯冲形成了日本海沟到东亚大陆广泛地区的沟槽弧后系统，以及强烈频繁的地震活动和大规模构造岩浆和火山活动.许多全球和区域成像研究揭示了太平洋板块从日本海沟向地幔转换带俯冲，停滞在地幔转换带中的太平洋板片与其上面上地幔区域形成一个巨大的地幔楔(Zhao，2004；Lei and Zhao, 2005；Huang and Zhao, 2006；Wei et al., 2012；Chen et al., 2017).

图 2

Fig. 2

图 2 反演的剩余密度异常剖面(37—44°N) 黑色箭头指示大兴安岭断裂带位置，虚线对应410 km和660 km深度面. Fig. 2 Inverted residual density profiles (Latitude 37—44°N) The black arrow indicated the location of the Da Higaan Mountains fault zone and dashed line corresponds to the depth of 410 km and 660 km.

图 3

Fig. 3

图 3 反演的剩余密度异常剖面(45—52°N) 黑色箭头指示大兴安岭断裂带位置，虚线对应410 km和660 km深度面. Fig. 3 Inverted residual density profiles (Latitude 45—52°N) The black arrow indicated the location of the Da Hinggan Mountains fault zone and dashed line corresponds to the depth of 410 km and 660 km.

图 2剖面显示俯冲的太平洋板片在日本海沟处呈现高密度异常，向西横向持续扩张，深度方向上有逐渐增加趋势.高密度体俯冲至地幔转换带(410~660 km)后在水平方向上向西俯冲，部分滞留板片进入下地幔.俯冲板片在横向上密度分布不均匀，与天然地震结果对应并且和前人研究一致(田有等，2011；Tian et al., 2016).

为了验证反演密度模型的可靠性，通过建立埋藏较深的板状体模型来模拟滞留板片进行正反演测试，板片模型大小为1000 km×600 km，深度范围为500~650 km，剩余密度设置为0.05 g·cm^-3.反演的垂向网格间距设为50 km，图 4a为板片模型的空间分布，图 4b为正演测试数据的反演结果垂向剖面图，图中白色方框为板片模型的真实位置，可以发现反演模型的空间形态恢复的较为准确.图 5给出模型和实际数据迭代反演的拟合均方根误差曲线，图 5a和图 5b分别是板片模型测试数据和东北实际观测数据反演的迭代残差曲线，能够反映反演过程的收敛情况.图 5c给出(如图 2和3)反演密度模型的计算预测数据与该地区实际观测数据的散点对比图，图中的直线为斜率为45°的辅助线，说明三维密度模型计算出的重力数据与观测重力异常拟合情况较好，进一步验证重力反演是有效及适用的.

图 4

Fig. 4

图 4 板片模型测试图 (a)模型空间分布图；(b)板片模型反演垂向剖面. Fig. 4 Test result of slab model (a) Spatial distribution of model; (b) Vertical profile of inverted slab density model.

图 5

Fig. 5

图 5 反演拟合误差分析图 (a)迭代均方根误差(模型); (b)迭代均方根误差(东北实际数据); (c)观测数据与反演预测数据对比(东北实际数据). Fig. 5 Fit error analysis of inversion (a) RMS error of inversion iteration (model); (b) RMS error of inversion iteration (actual data of NE China); (c) Comparison of estimated and observed data (actual data of NE China).

通过搜集到的1985个地震台站记录的14163次区域和局部地震的244180个P波走时数据和由694个台站记录的9988个远震事件的319857个相对走时残差数据，利用近震远震联合层析成像方法(Zhao et al., 1994；Ma et al., 2018)，反演过程中近震采用走时数据，远震采用相对走时残差数据，反演网格水平方向为1°×1°，清晰地得到了研究区域地壳与上地幔三维P波速度结构(图 6).

图 6

Fig. 6

图 6 P波速度结构剖面图(37—52°N)(Ma et al., 2018) 黑色箭头指示大兴安岭断裂带位置，虚线对应410 km和660 km深度面. Fig. 6 P-wave tomography profiles (Latitude 37—52°N)(Ma et al., 2018) The black arrow indicated the location of the Da Hinggan Mountains fault zone and dashed line corresponds to the depth of 410 km and 660 km.

天然地震层析成像结果与三维密度结构结果在形态上是对应的.从图 2的41—44°N反演结果剖面和图 3的45—46°N反演结果剖面可以看出近似水平运动的高密度体在松辽盆地(120—126°E)下方的200~600 km的区域，部分位置板片明显显示是不连续的，这一点与图 6对应位置的速度结构相吻合.从图 2和图 3的所有剖面都可以看出高密度异常在大兴安岭断裂带及以西部分地区的地下500 km附近的深处仍有高密度体存在，在天然地震的速度结构上体现不明显，我们推测这可能是大兴安岭断裂带下方的地幔岩石圈拆沉的结果.最终，我们推测得到板片俯冲动力学模型(图 7)，主要是来说明呈高密度高速异常的西太平洋板片深俯冲至地幔转换带后，继续西向俯冲，板片前端已到达大兴安岭造山带附近，部分滞留板片可能进入下地幔；板片深俯冲过程中可能会引发局部对流，导致幔源物质上涌.以上模型的推测与前人的研究结果相吻合(Chen，2010；Guo et al., 2016；Chen et al., 2017).

图 7

Fig. 7

图 7 板块俯冲动力学模型图 Fig. 7 Plate subduction dynamics model

本文利用东北亚大范围的重力数据反演获得了研究区域上地幔三维密度结构，与天然地震获得三维速度结构联合解释，共同揭示欧亚大陆下西太平洋板片的俯冲形态和动力学机制，能够提高单一方法的可靠性和解释准确性.反演获得上地幔密度结构表明在浅部(100~300 km处)大兴安岭造山带区域主要表现低密度低波速异常，松辽盆地表现较弱的高密高速异常并伴随着部分低密度低波速异常区域；在深部(300~700 km处)松辽盆地和大兴安岭造山带附近仍以高密高速异常为主，推测是下方地幔岩石圈拆沉的结果.板块俯冲的动力学模型展示呈高密度异常的不连续太平洋俯冲板片俯冲至地幔转换带后在水平方向上向西俯冲，前端已到达大兴安岭造山带附近，部分滞留板片可能进入下地幔.