小波多尺度分析在青藏高原东缘重力场的应用

刘 威, 王绪本, 王向鹏, 张 翔

(成都理工大学 地球勘探与信息技术教育部重点实验室，成都 610059)

0 引言

青藏高原东缘是青藏高原与扬子板块接壤的过渡带,受印度板块北向地应力,青藏高原内物质在应力作用下向川滇区块运移[1]。该区域发生强烈变形,发育多条大型南北向深大断裂带,是我国大陆构造运动最为强烈的区域之一。近代以来，该区发生多起特大地震事件，是国内、外学者研究的热点区域。学者们先后利用不同方法对该区域构造背景成像及物质运移动力学进行相应的研究，丁志峰等[1]、王椿镛等[2]、白志明等[3]利用地震反射剖面资料对青藏高原东缘及邻近区域上地幔地壳反射结构及P波波速结构进行分析研究，探讨了地壳厚度变化规律及深部构造运移背景；孙洁等[4]、王绪本等[5]、白登海等[6]利用大地电磁剖面针对青藏高原东缘各构造地块电性结构特征进行研究，揭示了深部电性结构及提出了下地壳流模型，为动力学研究提供了科学依据；蒋福珍等[7]、孟小红等[8]、杨文采等[9]利用重磁资料对东缘及邻近区域重磁异常进行分析，基于重磁场研究了横向及纵向密度变化特征，并对邻区内主要构造背景进行综合分析。前人在此也取得了丰富的研究成果以及科学认知，但由于方法及技术的应用不同，产生了不同的看法，对此笔者利用小波多尺度分析方法对卫星重力数据进行位场分离处理，结合径向平均功率谱对不同界面的重力异常以及莫霍面起伏进行分析。重力异常是地球内部不均匀密度物质引起的不同深度重力场的叠加场效应，通过实现有效的位场分离，提取不同深度、不同介质因密度不均而引起的重力场，杨文采等[10-14]利用小波多尺度分析进行了有效的位场分离分析研究，并取得了较好地应用。

本文数据来源于WGM2012地球重力场模型(地形和重力资料来自美国加州大学圣迭戈分校斯克里普斯海洋研究所http://topex.ucsd.edu/)，利用多尺度小波分析提取不同尺度的重力异常，采用径向功率谱分析提取异常的位置，进行莫霍面界面起伏反演，从而可以分析深部构造[6-9]。

1 地质构造及重力场数据

在四川盆地的阻挡作用下,青藏高原东向运移的过程中,在东缘形成高原的强边界,在川滇地区形成弱边界[15],使得青藏高原构造背景复杂,地震活动频繁，研究区范围为98°E～104°E，26°N～34°N，区域地形起伏较大，向东地势急剧下降，构造复杂，自由空气重力异常与地形起伏相关。研究区内分布众多断裂带，分别为龙门山断裂，鲜水河断裂，哀牢山-红河断裂，澜沧江断裂，以及对青藏高原主要隆升起控制作用的金沙江缝合带，班公湖-怒江缝合带，雅鲁藏布江缝合带，该研究区地壳厚度显著增大，具有较明显的重力梯度带。

这里使用的WGM2012地球重力场模型，是基于EGM2008重力场模型和ETOPO1模型基础的，网格数据为2’×2’，球谐系数扩展至2159阶。根据位场理论，利用引力位来计算重力场，地面重力异常(自由空气重力异常)与扰动位为一阶偏导的关系，自由空气异常表示为式(1)。

(1)

图1 自由空气重力异常图

Δg中间层=-2πfMρΔh=-0.1118Δh

(2)

(3)

Δg布格=Δg自由+Δg中间层+Δg地形

(4)

式中：ρ为地表岩石密度，2.67g/cm3; Δx、Δy为方格网中节点网格距;rij为测点与节点的间距，hij为测点与节点的高程差，cij为数值积分系数，改正区内点系数cij为1，改正区内顶点系数cij为0.75，改正区边缘点系数cij为0.5，改正区外顶点系数cij为0.25。进行正常场改正和地形改正后的重力布格重力异常，可以较为精细地揭示研究区的深部结构，重力布格异常图如图2所示。

图2 布格重力异常图

经自由层校正后的重力布格异常，去除了地表起伏带来的异常影响，能较好地反映地下深部构造，研究区布格重力异常的异常值范围为-100 mGal～-680 mGal,重力异常值呈现自西北向西南逐渐增大的趋势，且异常呈串珠状分布，与研究区内断裂及地体分布有一定关系。其负低异常分布在扬子地块，川滇菱形块体，腾冲地块；高负异常分布在青藏高原的巴颜喀拉块体，松潘甘孜块体以及川滇北部部分区块，且存在条带状负异常。

2 方法原理

2.1 二维小波多尺度分析

地球重力场受地下各种密度不均匀介质体的影响，布格重力异常中包含了不同深度，不同尺度的各种重力场源。在研究不同尺度的重力场时，要对这种复杂的场进行场分离，分别得到不同研究尺度的场。小波多尺度分析可以将信号按各种不同频率进行分解，可以在任意细节上体现不同细节的地球物理意义[10,16]。为了研究不同地质体引起的对应重力场特征，可以利用小波多尺度分析方法分离出不同尺度的地质体在横向和纵向分布的重力异常场[8,17]。

根据小波分析理论，假设二维重力异常场为：

Δg(x,y)=f(x,y)

(5)

重力分解表达式可简化为：

Δg(x,y)=ANG+D1G+D2G+…DNG

例如：在《黄鹤楼送孟浩然之广陵》相关内容教学中，小学生由于自身经历的局限性，难以真正理解文章中描述的内容，从而难以融入自身的情感，朗读质量也不高。针对这种情况，教师可以利用多媒体工具，为学生尽可能还原当时的情景，通过视频、音频、图片等直观的刺激，让学生有更加直观的认识，更好地融入到教学情境中进行朗读，自身的情感也能更好地融入，朗读兴趣与积极性也会随之提升。

(6)

式中:D1～DNG为异常的1阶到N阶小波细节；ANG为异常N阶小波逼近。图3为小波分解示意图。

图3 小波分解示意图

2.2 径向平均功率谱

功率谱分析是由 Spector等[18]最早提出来的一种重磁场解析处理方法，用来估测地下异常体平均深度。经傅里叶变换后的重力异常对其频谱进行对数功率谱分析，在极坐标系下重磁异常的径向平均对数功率谱曲线图可以反映异常场源的平均深度，定量估计重磁异常的场源深度，径向功率谱及场源埋深深度表达式为：

E(r1)=A2(r1)

(7)

(8)

式中：E(r1)、lnE(r1)为异常波谱的功率谱与径向对数功率谱；A(r1)为频谱幅值；r为圆波数；h为场源埋深。径向对数功率谱值与径向频率线性相关，线性直线斜率常用来推断估计不同尺度重磁异常的场源平均深度。

2.3 Parker迭代反演

Parker-Oldenburg位场迭代公式常用于重磁反演计算中，由于其采用了快速傅里叶变化和反演迭代，能计算物性横向变化的连续界面，计算的速度较快。冯锐等[25]对Parker-Oldenburg迭代反演方法进行详细的叙述，通过反演迭代式(10)计算迭代修正量Δhij，利用Parker正演式(9)计算下界面引起的重力异常，当前后两次迭代结果差值小于阈值时，截止迭代，输出界面深度。

(9)

(10)

(11)

式中:F[h(r)]、F[Δg(r)]分别是介质界面深度h(r)和重力异常Δg的傅里叶变换;z0为界面平均参考深度;ρ为界面密度差;r为位置矢量;k为波数;G为万有引力常数;MN为总计算点;l为迭代次数;∈为设置阈值。

3 重力模型试算

3.1 建立正演模型

笔者建立了由不同尺寸长方体组合的三层地质异常模型，用5个同深度不同大小的小长方体模拟浅层地质体；用2个中等尺寸的长方体模拟中部穿插构造单元，用两个相接不同埋深的大长方体模拟深部区域背景。模拟背景区域为32 km×32 km，模拟示意图及模型参数分别如图4和表1所示。

图4 模型示意图

表1 正演模型参数

3.2 模型多尺度分析

利用不同尺度长方体组合模拟地下结构分布，组合的模型正演重力异常如图5所示。在进行多尺度小波分解的过程中，小波基函数的选择对小波分解的效果有一定的影响，通常选取的小波基函数不仅要满足一维分解重构的准确性，而且还要能体现二维数据成像的优良性[19],正交小波可以保证低阶细节不变形，对称性影响小波分解后的重构。笔者选取了具有正交性和对称性的coif3函数，即能够突出一维分解重构的准确性，还能体现二维数据成像的优良性[20-21]。图6为coif3的尺度函数和小波函数图像，对模型组合重力异常进行六阶多尺度分析，多尺度分析结果如图7所示。

图5 模型正演重力异常

图6 coif3小波基函数图

通过比较正演模型与小波多尺度分析结果，小波多尺度分析的各阶细节能清晰反应不同深度尺寸的重力异常体分布。如图7所示，1阶～5阶小波细节组合可以清晰看到5个不同位置的异常体，形态与正演形态有差别，但位置较一致，较好地体现了局部异常；小波6阶细节(图7(c))体现出了中部穿插异常体的位置，而小波6阶逼近较好地分离出了深部背景场，将1阶～6阶小波细节进行结合，可以分离出中浅部异常体的分布。因此小波多尺度分析既可以有效对区域场进行异常分离，还能针对不同深度异常体分布进行有效分离，有助于分析地下介质的重力场源情况。

图7 模型小波分析

3.3 径向功率谱分析

经小波多尺度分析后的各不同尺度场源异常，利用径向平均对数功率谱对其进行深度估计。对计算后的径向对数功率谱数据的最低波数段进行线性拟合，由式(8)可知，拟合直线二分之一斜率为场源的平均深度值。1阶～5阶小波细节异常图的波谱图8(a)拟合直线斜率为4.86，反应了深度约为2.43 km的浅层长方体；6阶小波细节异常图的波谱图8(b)拟合直线斜率为13.88，反应了深度约为6.94 km的中部长方体；6阶小波逼近异常图的波谱图8(c)拟合直线斜率为36.34，反应了深度约为18.17 km的深部背景长方体。与模型数据对比，径向对数功率谱估计的平均深度较准确。

图8 模型径向功率谱

4 多尺度分析结果

基于式(6)对青藏高原东缘的重力布格异常进行多尺度小波分解，小波细节(D1G～D5G)和小波逼近(A5G)结果如图9所示。对图9(a)～图9(e)小波细节分别进行径向平均功率谱计算分析，分别如图10(a)～图10(e)所示。径向平均功率谱计算结果表明，1阶和2阶细节的场源深度分别为4 km和11 km，反映了上地壳沉积物质密度不均匀体，3阶和4阶的场源深度分别为26 km和37 km，反映了中地壳和下地壳物质密度不均匀性，5阶小波的场源深度为54 km，反映了莫霍面的起伏。

图10 研究区多尺度小波分析重力异常径向平均功率谱图

小波多尺度分解的1阶小波细节(图9(a))和2阶小波细节(图9(b))在松潘甘孜地块、川滇地块、腾冲地块异常比较复杂且剧烈，表明浅部存在高密度体,横向上受构造影响具有不均匀性。在三江区域内显示多条南北向条带状负异常，条带状负异常的极大值与该处的哀牢山-红河断裂，澜沧江断裂等南北向断裂吻合性较好，鲜水河断裂在1阶和2阶小波细节中的负异常值吻合性也比较好，揭示了该区域沉积物质高频异常的特征，表明了该区复杂的构造背景。

随着小波分解阶数的增大，分解的重力异常的场源深度也随之变大，重力异常逐渐变得缓和与平滑。地壳内部不同物质密度分布不均会导致重力场出现大小不一的正负异常，通常情况下，正异常与正密度异常相关，负异常与负密度异常相关。3阶小波细节(图9(c))和4阶小波细节(图9(d))主要反映青藏高原东缘区域中地壳和下地壳物质密度的不均一性，如深大断裂引起的密度差异导致的异常，龙门山断裂附近存在正负异常错综复杂的条带，对研究区重力场的分区特征有所反映，可揭示青藏高原东缘深部存在不同性质的块体之间的碰撞与融合。

5阶小波细节(图9(e))的场源深度与研究区内莫霍面的平均深度接近，异常值从西至东呈现“高低高”现象，在局部区域出现团状负异常和团状正异常，正异常主要分布在松潘甘孜地块，四川盆地所在的扬子地块，巴彦克拉西侧，表现了刚性块体特征。负异常主要分布于川滇地块。正异常形成的原因有可能是上地幔高密度熔融物质沿断裂通道上涌，负异常可能是深部软、热物质向东运移所产生的。通过分析不同异常产生的原因，可以帮助我们更好地去理解深部块体间的物质运移以及深部大尺度的异常源。

图9 研究区多尺度小波分析图

小波5阶小波逼近(图5(g))可反映深部莫霍面起伏引起的重力异常，为研究区内区域场，1阶～5阶小波细节之和(图9(f))可作为场分离后的局部场。相比布格异常图，消除了区域背景场后的局部异常，异常值分布有所变化，变化呈现中间低两边高的异常趋势，相比重力布格异常的梯度变化而言，局部异常更趋向于片区状以及串珠状分布，在横向上变化明显，体现了莫霍面以上物质分布不均而引起的重力场变化，在川滇地块与扬子地块，甘孜地块与扬子地块的接连处重力异常皆有明显变化，且存在与川西地区地震事件具有紧密联系的鲜水河断裂和龙门山断裂，有可能是因为青藏高原深部热物质向东逃逸过程中延断裂通道向上运移且在不同深度环境下发生聚集，形成局部高密度异常体，在研究区域西南边出现一系列重力异常的过渡带、梯度带、突变，呈现局部狭窄的北西向异常带，极有可能是此处俯冲断裂所引起的成岩物质所导致的。

5 莫霍面结构分析

莫霍面为地球内部与地幔的分界面,是引起区域重力异常的深部主要分界面,通过对莫霍面的确定,可以为认识地壳与地幔的构造发育以及构造运动提供参考。5阶小波逼近见图5(g)，利用径向功率谱计算的平均场源深度为52 km，与测深地震识别的莫霍面深度相似[22]，异常的趋势与莫霍面的趋势具有相似性，可以反映莫霍面起伏形态。地球内部物质密度分布影响着地球外部重力场，地壳与地幔间的密度差为重力资料确定莫霍面的位置提供了良好的基础。综合前人在该区域及附近区域的研究成果及地质背景资料[9,23]，利用5阶小波逼近(图9(g))采用Parker-Oldenburg界面反演方法，对研究区莫霍面深度进行迭代反演，选取重力基本场为5阶小波逼近，参考深度z0=50 km，壳幔密度差的基本参数，反演结果与CRUST1.0地壳模型一致性较好[24]，莫霍面界面埋深范围约为37 km～59 km，莫霍面深度从西北至东南逐渐抬升，地壳厚度逐渐减薄(图11)。

图11 研究区莫霍面深度图

基于前人研究的深部地震以及大地电磁探测资料背景，将研究区域内的地壳结构分为沉积层、上地壳、中地壳、下地壳4 层，地壳各层及上地幔的初始密度分别为：沉积层为2.3 g/cm3、上地壳为2.60 g/cm3、中地壳为2.70 g/cm3、下地壳为2.80 g/cm3，上地幔为3.30 g/cm3[25]。使用Oasis Montaj平台中的GM-SYS模块对东西走向，横跨川滇地块和杨子地块的剖面AB进行人机交互重力密度拟合反演，调整界面深度以及模型密度参数，拟合效果较好(图12)。松潘甘孜，巴彦克拉地块的莫霍面深度在52 km以上，川滇地块莫霍面深度主要分布在47 km左右，与地震测深资料获得的地壳厚度图形态且厚度相一致，展现为西北至东南逐渐变浅的趋势，且青藏块体与川滇地块和扬子地块深度变化较大，形成了环青藏高原莫霍深浅变化的梯级带状。龙门山断裂附近莫霍面深度差较大，在该区域发生多起特大地震灾害，青藏高原深部物质向东逃逸过程中在刚性四川盆地的阻挡作用下，物质运移能量发生聚集且延龙门山构造进行释放，从而消耗掉了板块之间巨大的碰撞应力。川滇块体南部莫霍深度带较为轻缓地向东南渐变，与下地壳流方向基本一致，软弱的下地壳物质延南向、东南向流动，导致地壳厚度大于刚性体的四川盆地地壳厚度。青藏高原东缘地壳的东西部厚度差异，可以合理地解释青藏高原物质东流的channel flow模式，四川盆地刚性块体阻挡了高原物质的东流，高原物质在此堆积挤压，形成了龙门山局部挤压推覆构造带，形成的推覆构造带吸收了来自西部青藏高原的能量，这也是青藏高原东缘地震能量的主要来源。

图12 AB剖面密度模型拟合

6 结论

1)使用WGM2012地球重力场模型，解译计算后的重力数据基于多尺度小波分析应用于青藏高原东缘，结合径向平均功率谱计算不同尺度小波细节所揭示的地下场源平均深度，参考前人研究成果以及各种地质资料，分析不同尺度的异常与地质背景间的对应。

2)分解后的1阶小波细节反映了上部沉积层物质密度横向不均匀变化，3阶～4阶小波细节揭示了地壳内部物质在深部动力学影响下，物质纵向变化及横向扩展而导致的不均匀性；5阶小波细节及5阶逼近反映了莫霍面起伏引起的区域场异常变化，能够揭示地壳厚度总体性变化。

3)采用Parker-Oldenburg界面反演进行重力莫霍面反演，并结合深部地震研究进行重力剖面拟合验证，反演的莫霍面深度与测深地震所探测的地壳厚度趋势一致，呈梯度向东南减薄，间接说明在推覆挤压的过程中深部物质向东逃逸，与刚性块体接触为地应力集中及释放提供场所，导致龙门山等地活动构造激活，发生多起重大自然灾害。