自20世纪90年代早期华北岩石圈减薄概念提出以来,有关克拉通破坏的研究越来越热.尤其是近几年以来,国内外地球物理､地球化学以及地质学家针对华北克拉通破坏的相关问题开展了大量的研究工作^[1-6],使得我们对华北克拉通的深部结构､形成及演化有了更清晰的认识.然而,由于学科背景的差异和研究思路的不同,目前对华北克拉通的时空分布,破坏机制和地球动力学原因等问题科学家们尚存在争议^[7].华北克拉通地质构造较为复杂,周边分布着阴山—燕山造山带,祁连山造山带和秦岭—大别造超高压变质带.华北克拉通西部为鄂尔多斯块体及周围盆地,中部有太行山造山带和汾渭地堑,东部以华北盆地和郯庐断裂带为主.本文研究区域为郯庐断裂带及其以西克拉通地区,区内地形及主要构造如图 1所示.

图 1

Fig. 1

图 1 研究区地形及主要构造特征(据文献[23]和[32]修改) Fig. 1 Topography and major geological features of the study region (Modified from Ref. [23]and[32])

目前,就华北克拉通破坏机制而言有两种主流观点:拆沉作用和热侵蚀作用.无论哪种破坏机制,都涉及到了华北克拉通岩石圈物质的运移和转换,由此,获得高分辨率的岩石圈结构尤为重要.近年来,随着地震层析成像技术的发展,人们对华北克拉通及其周边地区地壳､上地幔速度结构进行了大量研究,并取得了一些显著的研究成果^[8-16].与之相比,地下深部密度结构研究成果相对较少^[17].地球内部密度构造继承性对研究岩石圈内部构造变形和地壳物质组分^[18]起着重要作用.拆沉作用泛指由于重力的不稳定性导致岩石圈地幔､大陆下地壳或大洋地壳沉入下伏软流圈或地幔的过程^[19].岩石圈加厚迫使冷的､密度较大的岩石圈地幔沉入热的､浮力更大的软流圈中^[20].从这个角度来讲,研究华北克拉通岩石圈密度结构有助于进一步探讨华北克拉通的破坏机制,尤其是拆沉机制.

布格重力异常很好地反映了地下介质密度的不均匀性.由于重力测量具有海陆空测量的优势,随着重力卫星技术^[21-22]的发展和地面观测资料的积累,重力资料具有分辨率高,数据结构好的特点.虽然重力资料具有较高的分辨率,但重力方法反演解的非唯一性很强.为了对重力反演进行约束,我们将地震资料引入到重力反演,以改善反演解的非唯一性.本文在地震层析成像所得高分辨P 波速度结构基础上,运用华北地区速度与密度之间的经验转换关系,建立了华北克拉通岩石圈三维初始密度模型.在此基础上利用布格重力异常采用代数重构算法快速反演得到了华北克拉通岩石圈三维密度分布.

选取33°N—43°N,100°E—120°E 范围作为本文的研究区域.收集了研究区0.5°×0.5°P 波三维层析成像结果^[23].冯锐等在对华北地区大量人工地震资料研究的基础上提出了适合该区的波速-密度关系^[24]为

(1)

利用(1)式将收集到的P 波速度扰动转换为密度扰动,并以此作为密度反演的初始模型.同时,收集了研究区5′×5′空气重力异常资料(图 2)和30″×30″GTOPO30格网地形数据(图 1),经过地形改正和布格改正得到布格重力异常.

图 2

Fig. 2

图 2 空气重力异常(mGal) Fig. 2 Air gravity anomalies (mGal)

地表实测重力异常是地下各层界面起伏和物质密度不均匀的综合反应.因此,研究岩石圈密度不均匀性之前,需要先将其它因素引起的重力效应扣除掉.首先,我们收集了研究区0.5°×0.5°莫霍面和沉积层界面起伏资料,利用频率域界面Parker正演公式^[25]正演计算出莫霍面､沉积层界面起伏产生的重力异常.Bowin给出了场源埋深深度与重力位球谐函数阶数之间的表达式^[26]

(2)

式中,Z为场源埋藏深度,R为地球半径,n为球谐函数阶数.利用EEGM96模型计算2~30阶球谐函数对应的长波长部分的重力异常,认为该异常是有岩石圈以下物质密度不均匀引起的.最后,从布格重力异常中扣除莫霍面､沉积层界面起伏产生的重力异常和2~30阶卫星重力异常,得到岩石圈剩余重力异常.该异常可以认为是由岩石圈密度不均匀所产生的,本文将依据它进行密度反演.

模型设置:深度方向上划分为6 层,分别为10km,25km,42km,60km,80km 和120km;水平方向上,上面4层划分为0.25°×0.25°等间距格网,下面2层设置为0.5°×0.5°等间距格网.这样的模型设置方式即增加了深层块体在地表处的异常效应,又有效地减少了反演模型参数的个数.

对于密度反演中线性方程组的求解,较常用的方法是最小二乘法^[27-29].该方法虽然能够同时计算出分辨率矩阵和协方差矩阵,但需要对大型稀疏系数矩阵进行求逆运算和大量的矩阵相乘运算,这不仅需要占用大量的计算机内存空间,还大大降低了计算速度.随着观测技术的提高,观测规模的扩大和观测资料的积累,观测方程系数矩阵会越来越大,奇异性也会随之增强.为了克服以上缺点,适应超大型稀疏系数矩阵,我们将代数重构算法用于观测方程的求解.该算法采用行迭代技术,对初始估计的解进行反复修正,每次修正只针对一个观测方程进行,避开了大型稀疏矩阵求逆计算,有效提高了计算效率.对于密度反演观测方程

(3)

代数重构算法迭代反演的计算公式为

(4)

其中,ρ 和g分别为密度和重力异常,G为系数矩阵,K和λ 分别为迭代次数和松弛因子,λ 一般介于0~2之间.代数重构法在迭代过程中逐步改善待重建图像ρ 的初始值ρ0,每一步修正一个观测方程,每m步迭代成为一轮迭代,修正的依据是实际观测量值gi与对第K步迭代产生的图像的假想测量值$\sum _{j=1}^{J}{{G}_{ij\rho }}_{j}^{k-1}$之间的差值$\left( {{g}_{i}}-\sum _{j=1}^{J}{{G}_{ij\rho }}_{j}^{k-1} \right)$,这个差值按一定的方式分配到图像向量中,使图像趋近于其最终解^[30].经过8次迭代最终反演得到了华北地区10~120km 深度范围内6个不同深度上的密度分布,如图 3所示.图 4为迭代过程中反演密度结果正演所得重力异常理论值与观测值之间的均方差.图 5为最终密度结果正演的重力异常值与实测值之间的差值.

图 3

Fig. 3

图 3 岩石圈密度异常分布(g/cm3) A—10 km 深度；B —25 km 深度；C一 42 km 深度；D —60 km 深度；E —80 km 深度；F —120 km 深度. Fig. 3 Density distributions of lithosphere beneath North China (g/cm3) A—10 km depth； B—25 km depth； C一42 km depth； D—60 km depth；E一80 km depth； F—120 km depth.

图 4

Fig. 4

图 4 迭代过程中,反演密度结果正演所得重力异常理论值与观测值间的均方差(mGal) Fig. 4 Standard deviations between gravity anomaly theoretical values forwarded from density results and gravity data during the iterations (mGal)

图 5

Fig. 5

图 5 最终密度结果正演的重力异常值与实测值之间的差值(mGal) Fig. 5 The discrepancy of gravity anomaly between forward model from final density results and observed values (mGal)

太行山构造带总体为NNE走向,北部与燕山构造带交汇,南部与秦岭—大别构造带相连,是华北及中国东部地区一条重要的构造带.反演结果显示:60km 以上,太行山造山带整体以低密度异常为主,异常随深度增加不断减弱.同一深度上,异常大小和分布在造山带南段､中段和北段有所差异,表明太行山造山带岩石圈密度分布具有分段特性,这种分段特性可能说明了太行山造山带南､北､中三段经历了不同的地质过程.造山带中部转变为高密度异常,且与其东边克拉通东部高密度异常区相连,这暗示了华北克拉通岩石圈破坏可能已经由华北盆地延伸到该地区.42km 上深度上,华北盆地存在较弱的高密度异常,随着深度增加,高密度异常范围扩大,异常迅速增强.80~120km 深度上,华北克拉通东部岩石圈地幔表现为南北非均匀的强烈的高密度异常,这表明该地区遭到破坏后地幔物质分布具有强烈的南北向非均匀性.

郯庐断裂带是我国东部一条巨型断裂带,整体走向为NNE,其两侧地球物理场特征有明显的差异.我们的反演结果很好的体现了这一点,以断裂带为界,东西两侧密度异常存在明显差异.10~42km深度上,断裂带为低密度异常,异常大小随深度增加而减小.该深度上,断裂带西侧低密度异常显著,东侧虽也显示为低密度异常,但低密度特征明显弱于西侧.60km 深度上,断裂带密度异分布比较平缓,异常值小于0.005g/cm3.80~120km 深度范围内,断裂带密度异常特征由低异常转变为高异常,其两侧异常特征仍然存在明显差异,西侧异常大于东侧,这表明该断裂可能穿透了岩石圈.

鄂尔多斯块体为一相对稳定块体,由于该区自中､新生代以来构造变形和岩浆强烈活动,形成了鄂尔多斯盆地与周边特殊的构造格局^[31].地壳范围内,鄂尔多斯块体整体表现为低密度异常,不同深度上夹杂着小范围的高密度异常.10km 深度上,块体中部有一显著低异常区域;25km 深度上,块体北边的呼包盆地为高密度异常;42km 深度上,鄂尔多斯块体中南部显示为高密度异常,北部为较弱的低密度异常.接收函数研究^[32]表明,鄂尔多斯块体中部地壳厚度为41km,向北则加深至45km.由此可推断42km 深度上鄂尔多斯块体中南部已进入岩石圈地幔部分,从而表现出了高密度异常,而北部则仍然为地壳部分,因此表现出了低密度异常.块体南部的渭河盆地有一小范围的较强高密度区域.从60km 深度开始,鄂尔多斯块体密度异常由低异常转变为高异常,且一直延伸到岩石圈底部,高密度特征随深度增加明显加强.这与层析成像研究结果^[33]相一致.60km 深度上,高密度异常最大值出现在块体中部,80~120km范围内,鄂尔多斯块体南北两端各出现一个高密度集中区域,分别与其南部的秦岭造山带和北部的阴山造山带的高密度异常分布相连.鄂尔多斯块体的这种南北特异密度异常分布可能与其受到南北两边秦岭造山带和阴山造山带的深部动力学过程有关.祁连块体地壳部分为低密度异常,地壳以下岩石圈部分表现为高密度异常.低､高密度异常分布分别与秦岭造山带的低､高密度异常分布相连,而在与鄂尔多斯块体连接处,密度异常分布则存在明显的突变,这表明了岩石圈特性在鄂尔多斯块体和祁连块体及秦岭造山带之间存在着很大的差异.

4.1 华北克拉通岩石圈密度分布具有很强的横向非均匀性,东､西部密度特征有所差异.密度分布形态与地表构造形态具有很好的相关性.地壳范围内,华北克拉通以低密度异常为主,80~120km 深度范围内以高密度异常为主.

4.2 鄂尔多斯块体地壳密度异常分布具有整体性,表现为较弱的低密度异常.岩石圈地幔则呈现出南北两端各存在一个高密度异常的态势,这表明鄂尔多斯块体在深部可能受到了来自南北两端深部动力学过程的影响.

4.3 华北克拉通岩石圈减薄可能已扩展到太行山造山带中段,太行山造山带岩石圈地幔密度分布具有分段特征.现今华北克拉通东部上地幔物质分布具有强烈的南北横向非均匀性.

4.4 随着观测资料分辨率的提高,重力反演观测方程系数矩阵会越来越大.发展和完善更成熟的行迭代算法用于解算线性反演方程组,从而避开大型稀疏矩阵求逆运算,是提高反演计算效率的有效途径.

本研究所得到的华北克拉通岩石圈三维密度结构,为进一步探讨华北克拉通破坏作用的空间分布､机制和驱动力提供了基础资料.

特别感谢吉林大学地球探测科学与技术学院田有博士提供了本研究所用的三维P 波速度结构资料.