受印度洋板块与欧亚板块碰撞的影响，青藏高原产生了极其复杂的地球物理场，出现了地壳强烈形变、上地幔物质上涌以及地壳物质重新分异、组合与运移^[1]，成为研究大陆动力学的理想试验场，多年来受到国内外地学家的广泛关注.青藏高原东北缘地处青藏高原块体、鄂尔多斯块体和阿拉善块体交汇地带，区内主要构造为柴达木地块和祁连地块.自20世纪70、80年代以来，地学家通过地球物理探测对青藏高原及边缘地区地壳上地幔地球物理场特征和动力学特征进行了研究，对该区的构造变形，物质运动等进行了讨论.研究表明，青藏高原东北缘下地壳为低阻层^[2]，地壳整体上和上地壳一样呈酸性，巴颜喀拉块体和秦岭块体南部曾发生过下地壳拆沉，并导致高原的加速隆升^[3].青藏高原北部地区广泛发育新生代的火山岩^[4-5]，并普遍具有低Pn波速、缺失Sn波的特点^[6]，说明该地区岩石圈地幔具有异常的热结构，部分熔融正在发生，Sn波缺失区对应着青藏高原现代岩石圈拆沉和软流圈物质底侵的场所^[7].大地电磁探测结果显示，青藏高原中、下地壳普遍存在大规模的高导层，其成因可能与印度板块俯冲过程中由于摩擦生热导致的岩石部分熔融和脱水作用有关，也可能与岩石圈拆沉造成幔源物上涌有关^[8].各向异性研究表明，青藏高原东北缘SKS快波偏振方向呈现出顺时针旋转，这预示该区上地幔物质在区域构造应力场的作用下发生了顺时针旋转的形变以至流动^[9].青藏高原东部岩石圈地幔各向异性强度可达10%，综合分析流变学实验和岩石物理学实验研究成果，这种异常大的SKS分裂延时是由于熔体定向分布强化的结果^[10].熊熊等^[11]指出青藏高原东部地壳增厚主要发生在下地壳，下地壳低强度分布可能是导致该地区地壳、地幔运动解耦的重要原因.

综上可以看出，前人在探讨青藏高原东北缘壳幔结构及动力学机制时，涉及到了岩石圈拆沉、软流圈物质作用和部分熔融等因素.拆沉作用是由于岩石圈内密度差异而导致的重力坍塌.在拆沉过程中，加厚的下地壳转变成密度较大的榴辉岩与下伏的岩石圈地幔一起拆沉进入软流圈地幔.上涌的软流圈物质和部分熔融物质一般表现为较低密度异常.由此可见，岩石圈三维密度结构对验证青藏高原东北缘已有研究成果及进一步深化认识该区壳幔动力学机制具有非常重要的作用.本文综合重力资料和地震波走时资料，对青藏高原东北缘岩石圈三维密度和速度结构进行了反演，对密度结构进行了重点讨论.

在青藏高原三维密度结构研究方面，方剑和柯小平等基于地震层析成像所得三维速度结构，通过速度和密度之间的转换关系建立三维初始密度模型，然后反演得到了岩石层^[12]和地壳^[13]三维密度结构.两者反演结果分辨率都依赖于已有层析成像结果和重力观测数据的分辨率，所得密度结构横向分辨率分别为2°×2°和1°×1°，都比较低，且对反演结果仅作了初步分析.随着重力测量技术的发展和地面重力资料的积累，重力资料分辨率越来越高，为反演高分辨率、高精度的地下密度结构提供了条件.本研究采用5′×5′高分辨率网格重力数据.重力反演具有很强的非唯一性，必须加入有效约束才能得到改善.将层析成像所得速度模型转换为密度模型作为重力反演的初始模型来对重力反演进行约束，在国内外得到广泛应用^[12-14]，取得了较好的效果.本研究沿用以往研究方法，基于地震层析成像所得速度模型建立研究区初始密度模型，对重力反演进行约束.最终反演所得密度模型与速度模型具有较强的相似性，在某些区域和深度有些差异，这表明地震资料的加入对重力反演进行了有效地模型约束.

本研究综合重力资料和地震波走时资料，获取可靠的岩石圈三维密度结构.选取33°N-40°N，99°E-106°E范围作为本文的研究区域.图 1给出了研究区地形分布和区域构造背景.研究区由南到北分别分布着巴颜喀拉块体、柴达木块体、祁连块体和阿拉善块体.研究区东部覆盖我国南北地震带部分地区，东临古老的鄂尔多斯块体.

图 1

Fig. 1

图 1 研究区地形及主要构造特征 Fig. 1 Topography and major tectonics of the study area

本研究在进行地震层析成像时采用了区域地震走时和远震相对走时残差两部分资料.由于反演地球深部速度结构需要利用较长路径的地震射线^[15]，加之研究区区域台站和区域地震较少，故将地震层析成像的范围扩展到33°N-43°N，99°E-122°E（有关结果的重点讨论还是针对33°N-40°N，99°E-106°E的范围）.收集了上述区域585个地震台站接收的1980-2005年之间的7940个近震和12657个远震事件，总计149054个P波近震到时和193085个远震相对走时残差数据用于最终反演.其中，走时资料读取误差小于0.3s，大多数初至P波的拾取误差在0.1~0.2s，远震震中距均在30°~90°之间.本研究所用数据与田有进行华北克拉通地壳上地幔成像所用数据一致，反演结果可达到0.5°×0.5°的分辨率，详细的震中、台站分布，不同深度射线覆盖密度分析和不同尺度分辨率检测结果见文献[16].

收集了研究区5′×5′格网自由空气重力异常资料和30″×30″GTOPO30格网地形数据，经过地形改正和布格改正得到5′×5′格网布格重力异常.其中，地形改正采用严格积分算法.地表实测重力异常是地下各层界面起伏与密度不均匀的综合反应.因此，在进行岩石圈密度结构反演之前需要先把界面起伏引起的重力效应扣除.由于地表重力场包含了地表到地核物质密度不均匀的信息，所以还需要将岩石圈以下物质密度不均匀产生的重力效应扣除.首先收集了研究区0.5°×0.5°格网莫霍面和沉积层起伏资料，利用频率域Parker正演公式^[17]计算出莫霍面、沉积层界面起伏产生的重力异常.然后采用Bowin^[18]给出的场源埋深与重力位球谐函数阶数之间的表达式，利用EEGM2008模型计算得到2~33阶卫星重力场，即岩石圈以下物质密度不均匀产生的重力异常.最后从布格重力异常中扣除莫霍面、沉积层界面起伏产生的重力异常和2~33阶卫星重力场，得到岩石圈内密度不均匀产生的重力异常（如图 2所示），以此来进行岩石圈密度反演.

图 2

Fig. 2

图 2 岩石圈密度不均匀产生的重力异常 Fig. 2 Bouguer gravity anomalies of lithosphere result from density anomaly

首先利用地震波走时资料进行层析成像，得到三维P波速度结构；然后利用速度-密度经验关系式，将P波速度扰动转化为密度扰动，以此作为重力反演的初始密度模型；最后利用布格重力异常进行重力反演求取三维密度结构.

本研究采用赵大鹏发展的地震层析成像算法进行近震、远震联合反演.该算法采用有效的三维射线追踪方法计算走时及射线路径，在反演过程中对震源进行重新定位，并允许模型中存在速度间断面，目前在国内外已得到广泛应用^[19-21].图 3为本研究所用一维P波初始速度模型，该模型是在IASP91模型基础上，参考文献[22-24]的研究结果修改后给出的.

图 3

Fig. 3

图 3 一维P波初始速度模型 Fig. 3 The initial 1D P-wave velocity model used in this study

层析成像时深度方向设置了13层节点，深度分别为10、25、40、60、80、100、140 km和180 km，水平方向上节点间距为0.5°.利用上述资料和方法反演得到了研究区8个具有代表性深度上的P波速度扰动，如图 4所示.

图 4

Fig. 4

图 4 青藏高原东北缘岩石圈P波速度结构 Fig. 4 P-wave velocity structures beneath northeastern margin of the Tibetan Plateau

首先设置与地震层析成像相同的地下格网模型.重力反演分辨率随着深度增加而降低，为了改善这一情况，我们将初始密度模型中深部4层设置为水平方向1°×1°的格网，格网点值是由0.5°×0.5°的格网值平滑得到，然后利用Birch定律给出的速度-密度之间经验关系式^[25]，结合冯锐等^[26]的研究成果，将层析成像所得三维P波速度扰动转化为三维密度扰动，建立研究区三维初始密度模型.然后运用代数重构算法对线性方程组进行反演解算.该算法采用行迭代技术，每次修正只针对一个观测方程进行，避开了大型稀疏矩阵求逆计算，与传统的最小二乘方法相比，大大提高了计算效率.这样，经过6次迭代最终反演得到了青藏高原东北缘岩石圈8个不同深度上的密度结构.图 5为最终密度结果正演所得重力异常理论值与实测值间的差值.

图 5

Fig. 5

图 5 最终密度结果正演所得重力异常理论值与实测值间的差值（单位：mGal） Fig. 5 The differences between observed gravity anomalies and that forwarded from final density results (unit:mGal)

图 6给出了青藏高原东北缘岩石圈三维密度结构.整体来看，研究区密度结构具有较强的横向非均匀性，并随深度变化表现出差异性；密度异常在深度上表现出“低-高-低”的分层特性.由图 4和图 6可以看出，速度和密度分布形态与地表构造有很好的对应，且向下延伸到180 km深度上，这表明地表构造格局的形成受深部动力学过程的影响.中地壳以下，除80 km和100 km深度之外，速度和密度具有相似的分布特征，即同一区域同一深度上高速异常对应高密度异常，低速度异常对应低密度异常，两者等值线走向和形态也基本一致.

图 6

Fig. 6

图 6 青藏高原东北缘岩石圈密度结构 Fig. 6 Lithosphere density structures beneath northeastern margin of the Tibetan Plateau

10 km深度上，研究区零星分布着高密度异常，舟曲和武都表现为低密度异常，异常值较小.25 km深度上，研究区中部转变为低密度异常，异常范围和强度有所增大.祁连块体主要表现为低密度异常，研究区东部高低密度异常呈相间分布.

40 km深度上，祁连块体和柴达木块体均表现为强烈的低密度异常，低密度异常区域上下边界分别对应于海原断裂和日月山断裂，异常等值线走向与地表断裂带走向基本一致.体波层析成像^[27]和接收函数^[28]研究均揭示了该区地壳低速层的存在.随着深度的增加，60 km深度上祁连块体由低密度异常转变为强烈的高密度异常.地震测深研究表明，祁连块体莫霍面深度小于60 km^[29].层析成像研究显示了该区由地壳底部强烈低速异常到上地幔顶部强烈高速异常的转变，表明此区莫霍面上、下速度差异较大^[30].高密度异常表明该深度已进入岩石圈地幔.Sn波走时成像显示该区上地幔顶部S波速度明显偏高^[31].

80 km深度上，除石岗附近为低密度异常之外，祁连块体、柴达木块体和阿拉善块体均表现为较强的高密度异常.研究区东部和巴颜喀拉块体北部则为低密度异常.在此深度上，密度异常等值线走向不再与地表断裂带走向一致（东南-西北走向），而是发生了顺时针旋转，近似于南北走向.值得注意的是，密度异常与速度异常形态仍然有很好的对应，但密度异常比速度异常整体偏高，出现了强高密度异常对应弱高速度异常和弱低密度异常对应强低速度异常现象.例如，石岗-门源及研究区东部为较弱的低密度异常，速度图像则显示为强烈的低速异常；岷县以西分布着强烈的低速异常，却表现为微弱的低密度异常.

100 km深度上，祁连块体、柴达木块体和阿拉善块体分布着强烈的高密度异常，巴颜喀拉块体则表现为低密度异常.研究区东部分布着两个低密度异常区域和两个高密度异常区域.在此深度上也出现了同一地区密度异常比速度异常偏高的现象，一些地区甚至表现为高密度、低速特征.例如，大同西部分布着高密度异常，但速度结构则显示为强烈的低速异常.

140 km深度上，秦岭北缘断裂和海原断裂之间存在两个强烈的低密度异常区域，阿拉善块体表现为高密度异常.到180 km深度上，低密度异常范围有所减小，强度有所减弱.与之相反，阿拉善块体的高密度异常范围和强度则都有所增强.

整体来看研究区中上地壳高、低密度异常相间分布，下地壳表现为强烈的低密度异常.大地震多发生在高、低密度（速度）异常边界区域或者下地壳低密度（速度）层之上，并偏向高密度（速度）异常一侧.青藏高原东北缘地壳密度结构特征为该区地震孕育、发生和下地壳物质侧向流动提供了条件.该区下地壳低速异常被认为与深层热物质向上迁移，地壳底部发生部分熔融有关^[27].而李永华等^[32]通过对S波速度结构和泊松比进行研究，推断该区地壳部分熔融并不存在.Rudnick^[33]指出，当地壳内热流值超过40mW/m²时，岩石P波速度随深度增加而减小，热流值越高，速度减小的幅度也就越大.大地热流研究表明^[34-35]，祁连块体平均热流值为68.340mW/m²，说明该区具有高热流背景.因此，下地壳低密度异常可能与该区高热流和强大的构造挤压背景有关.

GPS和各项异性研究均表明，青藏高原东北缘上地幔物质在区域构造应力场作用下发生了顺时针旋转形变以至运动^[36-37].大陆岩石圈密度非均匀性及其构造继承性对大陆内部构造变形起着重要作用，密度结构记录了青藏高原东北缘地壳上地幔物质运动特征.由图 6可以看出青藏高原东北缘地壳（60 km深度之上）密度异常等值线走向与地表断裂走向基本一致，呈西北-东南走向.而进入地幔后，密度异常等值线走向发生了顺时针偏转，在研究区东部接近于南北走向.密度异常等值线的顺时针旋转现象，反映了在印度板块碰撞、俯冲和阿拉善块体、鄂尔多斯块体阻挡联合作用下，青藏高原东北缘深部形成了强大的区域构造应力背景，使得物质运动发生了顺时针旋转.利用重力异常计算的高原及邻区地幔对流应力场与地壳运动格局的明显差异表征了青藏高原东部地壳与地幔物质运动的解耦^[11].地壳与地幔密度异常等值线顺时针旋转现象表明，青藏高原东北缘地壳和地幔具有不同的构造运动模式，暗示该区可能发生了壳幔解耦.壳幔解耦与该区强大的构造挤压应力背景有关.岩石圈S波速度结构和各向异性研究^[38]表明青藏高原岩石圈地幔的构造运动具有与地壳不同的模式，软弱的下地壳提供了壳幔运动解耦的条件，与本文结论互相印证.

80 km和100 km深度上，密度异常相对于P波速度异常整体显著偏高，甚至出现了高密度异常对应低速度异常的现象.当温度升高或发生部分熔融，亦或岩石含水量增加时，都会使得地震波传播速度降低.地质调查表明青藏高原东北部广泛发育新生代火山岩，并且现在岩浆作用仍在进行.加上该区低Pn波和Sn波缺失特征，说明青藏高原北部莫霍面以下岩石圈地幔具有异常高的热结构，部分熔融作用正在进行.岩石圈地幔P波速度整体偏低可能与该区高热流背景或者岩石含水量增加有关.高热流背景能够促使部分熔融的发生，从而降低P波的传播速度.速度异常相对密度异常的显著降低可能暗示该区部分熔融的发生.140 km和180 km深度上，研究区均分布着强烈的低密度异常，上地幔顶部的台站间Pn波到时差成像结果^[39]显示青藏高原东缘的上地幔为显著的低速异常.低密度异常暗示了深部软流圈热物质的上涌.鄂尔多斯块体至今仍保持着克拉通属性，具有巨厚的岩石圈根（约200 km），阻碍了青藏高原深部软流圈物质向东流动，从而导致软流圈物质的上涌.

印度板块与欧亚板块的持续碰撞和俯冲，为青藏高原地壳上地幔水平运动不断地提供动力.受阿拉善块体和具有巨厚岩石圈根的鄂尔多斯块体的阻挡，在青藏高原东北缘形成了强大的区域构造应力背景，导致了软流圈热物质的上涌，从而形成了高热流背景.强大的区域构造应力和高热流背景使得青藏高原东北缘可能发生了壳幔解耦和P波传播速度的显著偏低.

本文利用重力观测资料和地震波走时资料对青藏高原东北缘岩石圈三维密度结构进行了反演和讨论，所得主要结论如下：

（1）青藏高原东北缘中上地壳高低密度异常相间分布，下地壳表现为强烈的低密度异常.这种密度结构特征为该区地震孕育、发生和地壳物质侧向流动提供了条件，下地壳的低密度异常可能与该区高热流背景有关.

（2）青藏高原东北缘地幔密度等值线走向相对地壳发生了顺时针旋转，表明地壳和地幔具有不同的构造运动模式，暗示地壳和地幔物质运动发生了解耦.

（3）80~100 km深度上青藏高原东北缘岩石圈地幔P波速度相对于密度显著偏低，表明该区可能发生了部分熔融或者岩石含水量的增加.

（4）印度板块向欧亚板块的持续俯冲、碰撞和来自阿拉善块体、鄂尔多斯块体的阻挡，使得青藏高原东北缘形成了强大的区域构造应力场，导致了软流圈热物质的上涌，并最终导致了该区壳幔物质运动解耦和岩石圈地幔P波速度显著降低.

三维密度结构为进一步研究青藏高原壳幔变形，壳内物质侧向流动和高原隆升机制提供了依据.

田有博士提供了本研究所用地震波走时数据，在此表示感谢.同时感谢三位审稿专家对本文所提的建设性意见.