自新生代以来，青藏高原受印度板块和欧亚板块碰撞挤压俯冲的影响出现了强烈的构造变形和物质运移.青藏高原及邻近地带地质构造、岩石层结构和地球物理效应等相关研究引起了国内外地学界的广泛关注(滕吉文等，1980; 曾融生，1992；傅容珊，1999；Kind et al., 2002; Clark, 2005).在国内外地学专家的努力下，对青藏高原岩石圈的三维结构、物质特性及深部动力过程的研究有了很大进展.

地震波层析成像显示壳内速度呈现明显的东西向条带分布，以几条缝合线为界，不同地块显示不同特征.接近高原内部地壳较厚，地壳中存在低速流变物质，是软地壳，高原周边区域地壳较薄，波速比小，可能是硬地壳(沈旭章等，2013).岩石圈地幔在印度大陆与青藏高原东、西部结点处较厚，中部减薄，平均速度较高，具有克拉通特性.大地电磁探测结果发现高原内部普遍存在低阻地壳及上地幔(赵国泽等，2004；金胜等，2010).地震反射资料及GPS观测数据显示印度板块向北挤压俯冲仍是青藏高原缩短隆升的主要动力(Zhang et al., 2004).高原东西部岩石圈变形特征存在明显差异，西部以北东向地壳缩短运动为主，东部以顺时针旋转为主，不同的物质结构导致了这种差异(嘉世旭，2008).从深反射地震剖面发现，塔里木块体的岩石圈与青藏高原岩石圈面对面双向俯冲碰撞(高锐等，2002).

前人对青藏高原岩石圈结构框架、物理性质及深部动力学过程的研究取得了很多成果，然而对壳幔结构、动力学机制及高原变形机制的认识依旧存在争议.地壳俯冲模式认为印度板块仅俯冲到青藏高原南部2/3距离，北部1/3可能发生了拆沉和热减薄作用；拆沉-底侵模式认为青藏高原3 Ma以来的快速隆升是拆沉作用导致软流圈物质上涌的结果(杨立强等，2001).青藏高原地壳广泛存在的低速、低阻、高泊松比物质，张乐天等(2012)认为是岩石圈拆沉和软流圈底侵造成的上地幔顶部低速物质侵入的结果；而Nelson等(1996)则认为是俯冲摩擦生热和含盐流体导致岩石部分熔融的结果；Royden等(1997)认为是源于下地壳流.上述争议和分歧，尚需要提供更多的深部地球物理证据.探测青藏高原地壳、上地幔结构, 揭露其内部构造是理解和建立青藏高原隆升机制与地球动力学过程的必然途径.密度构造是其中必不可缺少的一部分.

重力位场数据是地球内部密度结构最直接的反映.在青藏高原密度结构研究方面，方剑(1997)和柯小平等(2009)基于地面较为稀疏的重力数据，给出过横向分辨率分别为2°×2°和1°×1°密度结构，密度结构分辨率都比较低.王新胜等(2013)和Bai等(2013)综合地面重力和地震波速资料，得到过青藏高原东北缘地壳0.5°×0.5°，地幔1°×1°的岩石圈三维密度结构.但只是采用了重力，而没有加入重力梯度.重力梯度作为重力位的二阶导数能较好地反映高频的密度变化，提高密度结构的分辨率；重力梯度全张量中有五个独立分量可有效增加观测数据.

由于观测条件的限制，高原重力地面观测资料缺乏限制了该区域内部密度结构的研究.随着重力卫星的发展，具有灵敏度高、环境干扰小、高效、高精度、全球覆盖等优点的卫星重力梯度测量为地球内部密度的研究提供了一种有效手段(祝竺等，2015).通过频谱分析，重力梯度异常能反映相对波长短的高频信号，重力异常能反映相对波长长的低频信号.联合梯度张量数据的高频信号成分和重力异常数据的低频信号成分可以得到加强的重力异常数据，该数据能同时保留重力异常数据及其梯度张量数据的所有信息，拓宽重力异常数据的频宽(袁园，2015).因此，本文利用GOCE卫星Level 2沿轨道重力梯度测量数据产品，经过平均轨道、正常场改正，延拓内插后得到分辨率为15′×15′的高精度近地面重力梯度异常联合EGM2008地球重力场模型得到的同精度重力自由空气异常共同研究了青藏高原及邻区岩石圈密度结构.

本文首先对观测异常经过低阶项、地形效应和沉积层界面起伏效应的改正后，得到由壳、幔密度不均一所引起的剩余异常.然后，以改正剩余量作为观测值，以改正剩余量归一化权重组成的对角矩阵作为观测权阵，基于Tikhonov正则化理论建立反应目标函数.反演过程中，地震层析速度转化密度提供初始约束，通过非等权最小二乘迭代法计算得到青藏高原及邻区下方0~120 km深度范围内，0.5°×0.5°密度分布.最后，结合已有研究成果，分析所得密度结果反映的构造和运动机制，为该区岩石圈壳、幔运动特征及动力学特征的认识提供参考.

在区域三维密度结构反演中，通常将地下结构划分为直立长方体.在空间直角坐标系中，长方体在其外任意点(x₀, y₀, z₀)处产生的重力异常为

(1)

积分可得：

(2)

同理可得：

(3)

上面观测方程可简化并向量化为：

(4)

其中，T_p=[T_xx; T_yy; T_xy; T_xz; T_yz; G₀]为参与反演的联合观测向量；T_xx为重力异常X分量对X轴的变化率，即水平重力梯度异常T_xx观测量，同理，T_yy为重力异常Y分量对Y轴的变化率，即水平重力梯度异常T_yy观测量，T_xy为重力异常X分量对Y轴的变化率，即重力梯度异常T_xy观测量，T_xz为重力异常X分量对Z轴的变化率，即重力梯度异常T_xz观测量，T_yz为重力异常Y分量对Z轴的变化率，即重力梯度异常T_yz观测量.G₀为重力自由空气异常观测量；G_p为模型空间密度分布重力及重力梯度正演核函数矩阵；ρ为岩石圈剩余密度列向量；G为万有引力常数.

首先利用2009年11月1日至2011年6月9日时间段采集的GOCE卫星GO_CONS_EGG_TRF_2沿轨重力梯度观测数据(图 1a)，经过5点模型梯度法平均轨道归算、球谐系数为2~20阶的GRS80参考椭球正常重力梯度改正(图 1b)，迭代法重力梯度向下延拓(Sebera et al., 2014)和最临近点数据内插计算，得到局部指北坐标系(LNOF)(Z指向地心，X指北向，Y指西向)下，22°N—47°N, 75°E—114°E范围内，观测高度20 km，观测分辨率15′×15′的重力梯度异常(图 1c).

图 1

Fig. 1

图 1 (a) GOCE L2 Tzz; (b)经过平均轨道归算和正常场改正后的Tzz; (c) 20 km高度处的Tzz Fig. 1 (a) GOCE L2 gravity gradient Tzz; (b) Tzz after average orbital reduction and normal field correction; (c) Tzz after downward continuation to altitude 20 km

经过对比计算，GOCE L2延拓近地面重力梯度异常分辨率大约对应EGM2008模型720阶的重力异常，采用EGM2008地球重力场模型计算相同观测范围和观测高度，球谐系数为2~720阶的自由空气重力异常G₀.

以往重力及重力梯度联合反演研究中很少顾及重力和重力梯度两种不同类观测数据权重的问题(秦朋波和黄大年，2016).观测数据权重的合理确定不仅可以有效避免观测数据数量级相差过大，计算中大数变化会掩盖掉小数变化的问题，而且可以有效提高目标函数的收敛速度.根据赫尔默特方差分量估计思想，依据验前估计确定观测值权重的方法，选取任意观测量标准差作为单位标准差，并以此单位标准差与其余观测量标准差的比值作为各个观测量的权重可以有效解决上述问题.本文选择t_xx分量标准差std(t_xx)作为标准差，以此获得改正剩余量观测权阵

其中std(t_xx), std(t_yy)，std(t_xy)，std(t_xz)，std(t_yz)和std(g₀)分别为改正剩余重力梯度异常t_xx, t_yy，t_xy，t_xz，t_yz和g₀分量标准差.根据Tikhonov正则化理论和正演矩阵方程(4)，选取数据拟合残差最小和模型长度最小作为稳定泛函，最终建立顾及相对观测权重的联合反演模型Φ(ρ)(Tikhonov and Arsenin, 1977)：

(5)

其中，μ为正则化参数，由L曲线得到.当时，模型求解达到最优，求得.

由于反演问题的不适定性，运用非等权的迭代最小二乘法进行迭代求解：

(6)

其中ρ为迭代计算密度异常矢量，首次迭代密度由地震层析S波速度V_S，根据S波速度-密度转换关系式和标准的岩石圈密度结构得到.初始密度提供了有效的形状约束，减少了反演的多解性，提高了反演结果的可靠性.

迭代的过程就是求解最佳拟合观测数据异常模型参数的过程.以数据拟合差t_p－G_pρ的均方差小于或等于一个小的常数ε₁，或是密度改正比例Δρ/ρ不大于ε₂，其中Δρ为每次迭代密度修正量，ε₂为某一小的常数，且ε₂＞0，或是迭代次数不超过最大迭代次数k_max作为迭代终止条件.

最后，运用模型分辨率矩阵R_ρ对最终解算结果进行评价.根据反演基本理论，当R_ρ=I模型参数会被唯一确定.R_ρ矩阵对角线元素越接近于1，理论解算能力越好，结果越优.

(7)

研究区域范围为75°E—114°E, 22°N—47°N(图 2).几乎覆盖了青藏高原的全部，以及高原东部的四川盆地、东北部的鄂尔多斯地块、北部的昆仑—柴达木地体和东南部扬子地块的一部分.所选区域由不同地质时期的构造单元拼合而成(滕吉文等，1980).

图 2

Fig. 2

图 2 研究区域地形及主要构造特征 红色框内是主要密度研究区域; 运动速率参考http://sps.unavco.org/crustal_motion). HM—喜马拉雅块体；LSB—拉萨块体；QTB—羌塘块体；QB—柴达木盆地；SG—松潘—甘孜块体；SCB—四川盆地(任纪舜等, 1997). Fig. 2 Topography and major tectonics of the study area Red box is the main area for density study. HM—Himalaya block; LSB—Lhasa block; QTB—Qiangtang block; QB—Qaidam block; SG—Songpan-Garzê block; SCB—Sichuan basin (Ren et al., 1997).

第3节数据处理所得观测重力及重力梯度异常值是地球内部密度异常和密度界面起伏的综合影响的反映.为了得到0~120 km岩石圈密度不均匀的重力异常和重力梯度异常，需要分别扣除地表地形、沉积层密度界面起伏和120 km岩石圈以下物质不均匀的影响.

(8)

其中：T_p，T_p¹，T_p²，T_p³和t_p分别为重力和重力梯度观测值、低阶项影响、地形效应、沉积层效应和改正剩余量组成的列向量.

利用场源深度与重力位球谐系数阶数关系表达式(Bowin, 1983)：其中R为地球平均半径；n为球谐系数阶数；d_n为场源最大埋深.扣除2~44阶球谐系数计算的岩石圈以下物质不均匀引起的重力及重力梯度异常.收集研究区30″×30″GTOPO30格网地形数据和0.5°×0.5°格网沉积层起伏资料，利用直立长方体正演公式(2) (3) 计算表面地形质量和沉积层界面起伏产生的异常.然后利用改正得到的剩余量参与密度反演计算.改正剩余量如图 3a—3f.

图 3

Fig. 3

图 3 (a)重力梯度异常Txx; (b)剩余重力梯度异常Tyy; (c)重力梯度异常Txy; (d)剩余重力梯度异常Txz; (e)重力梯度异常Tyz; (f)剩余重力梯度异常G0 Fig. 3 (a) Gravity gradient anomalies Txx; (b) Residual gravity gradient anomalies Tyy; (c) Gravity gradient anomalies Txy; (d) Residual gravity gradient anomalies Txz; (e) Gravity gradient anomalies Tyz; (f) Residual gravity gradient anomalies G0

模型在深度方向上设置10, 25, 40, 60, 80, 100 km 6层节点，水平方向上节点间距设置0.5°×0.5°，建立格网模型.基于直立长方体正演公式(2) (3)，根据公式(5) 建立反演目标函数.反演流程如图 4所示.

图 4

Fig. 4

图 4 反演流程图 Fig. 4 Flowchart of inversion

由S波速度V_S和S波速度-密度转换关系式ρ=0.049V_S^0.5(Miller et al., 1991)以及标准的岩石圈密度结构

ρ_ref=2.7 g·cm^-3，0＜depth≤15 km；

ρ_ref=2.7 g·cm^-3，15＜depth≤35 km；

ρ_ref=2.8 g·cm^-3，35＜depth≤50 km；

ρ_ref=3.3 g·cm^-350＜depth≤120 km，

(Artemieva, 2006)得到初始密度异常.其中V_S的单位是m·s^-1，转换密度ρ的单位是g·cm^-3.

然后，应用非等权最小二乘方法对改正后的剩余异常进行三维密度迭代反演.第一步，给定估计密度异常的初值ρ，源于地震波速转换；第二步，将ρ代入公式(6)，计算剩余异常和计算异常的差值；第三步，选取数据拟合中误差ε₁≤10^-3，密度改正比例ε₂≤10^－4，最大迭代次数k_max≤100，判断计算结果的收敛性.如果是，即停止计算；否则计算修正密度值Δρ和新的密度估计值ρ，转到第二步，以本次迭代的值作为初始值进行下一步迭代.最终迭代收敛拟合残差分布如图 5a—5f.

图 5

Fig. 5

图 5 最终密度结果正演所得异常值与实测值差值分布直方图及标准差 Fig. 5 Histogram and standard deviation (STD) of the differences between observed gravity and gravity gradient anomalies and those from forward modeling

从图中可以看出，最终密度异常正演所得异常理论值与实测值之差均呈现正态分布.剩余残差主要是高频误差.对角线分量T_xx、T_yy的残差值(±0.2E，±0.3E)大于非对角线分量T_xy、T_xz、T_yz的残差值(±0.06E，±0.04E，±0.1E).重力Gz残差值为最大，约为±2 mGal.T_xx, T_xy, T_xz, T_yy, T_yz, Gz各分量均方根误差分别为0.07E，0.016E，0.015E，0.11E，0.029E和0.61 mGal.

由图 6a和图 6b可见，10 km深度，高、低密度异常相间分布.塔里木盆地、柴达木盆地和印度板块俯冲前陆盆地均表现为强烈的低密度(~2.55 g·cm^-3)特征，表明该区域沉积层较厚，这个低密度异常是沉积盆地低密度与周围结晶基底高密度之间差异的体现.25 km深度，羌塘块体北部存在中地壳低密度层(~2.75 g·cm^-3).接收函数反演(李永华，2006)和体波层析成像(丁志峰等，1999)均揭示了该区低速层的存在，这些低速层与其下部的高温上地幔有关.

图 6

Fig. 6

图 6 (a) 10 km深度密度；(b) 25 km深度密度；(c) 40 km深度密度；(d) 60 km深度密度；e) 80 km深度密度，蓝色虚线表示密度起伏界面，红色实线是95°E经线；(f) 100 km深度密度 Fig. 6 (a) Density section of 10 km depth; (b) Density section of 25 km depth; (c) Density section of 40 km depth; (d) Density section of 60 km depth; (e) Density section of 80 km depth. The blue dash line is density undulation interface. The red solid line is longitude 95°E; (f) Density section of 100 km depth

图 6c中，40 km深度，青藏高原内部为中地壳，表现为低密度(~3.05 g·cm^-3)，邻区为中下地壳，表现为高密度(~3.15 g·cm^-3).高原内中下地壳强烈的低密度(~3.0 g·cm^-3)围绕高原边界分布.鄂尔多斯块体和四川盆地块体密度异常较地震波速异常明显偏低.龙门山断裂带两侧高低密度差异更加明显并一直向下延伸到60 km.图 6d所示，60 km深度，青藏高原密度约为3.15 g·cm^-3，邻区密度较高约为3.3 g·cm^-3，原因是在该深度上，青藏高原为下地壳，周围邻区已进入地幔.

如图 6e所示，80 km深度，高原内部东、中、西密度特征差异明显.低密度以95°E为中心线呈东西对称分布.以班公—怒江缝合带为中心，在拉萨块体和羌塘块体内从北向南出现了“低-高-低”的密度分布特征.青藏高原上地幔的顶部显示出低密度的特征.四川盆地和鄂尔多斯盆地内，地壳高密度异常较地震波速异常明显偏低.由图 6f可见，100 km深度，高原内部密度与邻区逐渐接近.塔里木盆地内密度随深度均匀增加，从地壳到上地幔没有发现低密度层.

体波层析成像显示存在羌塘中地壳低速层，并指出这些低速层与其下部的高温上地幔造成的部分熔融有关(丁志峰等，1999).25 km深度密度结果表明该低密度层主要分布在该块体北部，暗示了北部比南部存在更广泛的壳内部分熔融.

壳内低密度层代表了壳内较为软弱的区域，我们密度结果显示壳内强低密度主要分布在高原边界，反映了高原边界软弱层的存在.张希(2004)对1993—2001GPS水平运动资料分析表明，印度板块每年以约5 cm的速度北向俯冲，强大作用力持续传到高原内部，受北侧阿拉善、东侧鄂尔多斯和南侧四川盆地三个坚硬块体阻挡，使得高原边缘形成了强大的区域应力场，高原边界为高应变积累闭锁区.这为深层热物质上涌、地壳物质部分熔融等提供了条件，从而导致高原边界与周围区域相比形成强低密度和密度梯度带.大地震多发生在构造边界和密度正负异常梯级带(陈石，2011).

80km深度，拉萨、羌塘块体内密度起伏和莫霍界面起伏特征一致，根据大地构造研究结果，羌塘块体向北倾斜，柴达木块体和巴颜喀拉块体向南倾斜.印度和羌塘块体不同时期分别向喜马拉雅和拉萨块体数次挤入并形成多个双向俯冲断层(曾融生等，1998).这种起伏特征也验证了印度、羌塘块体向喜马拉雅、拉萨地块挤入的双向俯冲模式.

塔里木盆地内从地壳到上地幔没有发现低密度层，对应地震波柴达木盆地接触部位基底及盖层介质均表现为高速异常的研究结果(赵俊猛等，2008)，表明此处应该是深部致密地幔物质向上迁移的结果.

本文利用地震波速转换的密度作为初始值，利用卫星重力和重力梯度联合反演了青藏高原及邻区岩石圈三维密度结构并对其动力学意义进行了讨论.所得主要结论如下：

(1) 青藏高原内部为低密度的中地壳，邻区为高密度的中下地壳.中地壳强低密度层主要分布在高原边界.印度板块俯冲和周围坚硬块体阻挡作用导致在高原边界形成了高应变积累闭锁区，为壳内低密度软物质的形成提供了条件.

(2) 80 km深度，青藏高原上地幔顶部显示出低密度的特征.高原内部东、中、西密度特征差异明显，低密度以95°E为中心线呈东西对称分布.以班公—怒江缝合带为中心，在拉萨块体和羌塘块体内从北向南出现了“低-高-低”的密度分布起伏特征.该特征与GRACE得到的莫霍面起伏特征一致，结合大地构造结果，这种起伏特征验证了印度、羌塘块体从南北两侧分别向喜马拉雅、拉萨地块挤入的双向俯冲模式.

(3) 塔里木盆地和柴达木盆地内壳幔高密度的结构特征，对应地幔物质上涌.四川盆地和鄂尔多斯盆地内，地壳高密度异常较地震波速异常明显偏低，表明古老的四川盆地和鄂尔多斯盆地比想象中更冷、更坚硬.

感谢S波速度的分享.感谢European Space Agency(ESA)提供GOCE数据产品.本文大部分图件使用GMT(Wessel et al., 2013) 绘制.