用于区域重力场定量解释的多尺度刻痕分析方法

2015-12-12杨文采孙艳云侯遵泽于常青

地球物理学报 2015年2期

杨文采，孙艳云，侯遵泽，于常青

1大地构造与动力学国家重点实验室，中国地质科学院地质研究所，北京 100037

2中国地质大学地球物理与信息技术学院，北京 100083

1 引言

区域地球物理学旨在应用先进的物理场观测技术，按规定网站密度观测区域地球物理场，取得关于区域地壳上地幔组成结构的信息，最终揭示区域地壳三维结构与构造形态，为资源勘探、地质灾害预防及生态均衡提供科学依据.由于区域地球物理场观测网度均匀，不受岩芯分布及岩石出露的限制，而且信息来自各种深度的场源，可以根据物理学定律推算地壳三维结构与构造形态，因此区域地球物理学研究日益受到重视.

地球物理场包括重力场、磁场、反射地震波场、宽频地震波场、地温场等；观测的方式有地面、航空和人造卫星三种.航空磁场是观测成本最低的地球物理场，但是它随离开场源的距离的平方较快衰减，因此仅反映上地壳物质磁化强度信息，信息量少.反射地震波场是观测成本最高的地球物理场，虽然可以提供的信息量最多，但因要求太高的投入对大区域研究不实用.重力场是观测成本较低的地球物理场，它随离开场源的距离衰减，因此可反映全地壳物质的密度信息，信息量较多.尤其是我国大部分地区都完成了1∶20万区域地面重力调查，周边地区也有卫星重力场作参照，利用区域重力场的信息，揭示区域地壳三维密度结构与构造分区，就成为现今区域地球物理学研究的主攻方向.不过，要实现揭示区域地壳三维密度结构与构造的目标，还必须研究系统的数据处理、信息提取和反演解释方法.

经过多年研究，我们把小波多尺度分析、表面刻痕分析以及位场频率域解释理论和密度扰动反演方法有机地结合起来，形成了应用区域重力场刻画地壳三维密度结构构造的数据处理、反演解释和信息提取的方法系统，分为区域重力场按场源分层、小波变换多尺度分析、场源分层深度及密度扰动反演、分层刻痕分析和构造边界定位四个子系统，详见图1.这一方法系统简称为区域重力场多尺度刻痕分析方法，本文将扼要介绍其方法原理及应用效果.

图1 区域重力场多尺度刻痕分析方法框图Fig.1 The flow－chart of the multi－scale scratch analysis for quantitative interpretation of regional gravity field

2 布格重力场小波变换的尺度场源深度转换律

1990年以后，我们首先进行了重力异常小波变换多尺度分析研究，取得了较好的效果.我们应用Mallat塔式算法（Daubechies，1990；Mallat，1989），构造出专门用于重力异常分解的小波基（侯遵泽和杨文采，1997；侯遵泽等1998）.二维离散小波变换的原理和区域重力场分解应用已在前文作过介绍（侯遵泽和杨文采，1997，2011，2012；侯遵泽等，1998；杨文采等，2001），本文不再复述.区域重力场可表为二维单值实函数Δg（x，y），地壳三维密度结构模型可表为三维单值实函数Δρ（x，y）.如果不附加约束性规律，仅靠地面二维数据是不可能取得三维密度结构模型的.首先要附加的约束性规律称为布格重力场小波变换的尺度－场源深度转换律，推导如下.

区域重力场小波多尺度分析方法的理论依据之一，乃是布格重力异常场的尺度与场源埋藏深度呈现同步增大关联.换言之，场源埋藏深度越大，布格重力异常场的水平尺度就越大.例如，对于球形密度扰动物体，重力异常表为（曾华霖，2005）

式中G为万有引力常数，M为场源质量；h为场源中心深度，x为观测平面上测点离场源中心的水平距离.定义布格重力异常场的水平宽度w＝2x′，x′为异常从极大值衰减为1／2极大值的水平距离，W即为孤立异常的特征尺度.把x′代入（1）式可以导出由此可见，场源埋藏深度h越大，观测平面上布格重力异常场的水平尺度W就越大，二者之间近似呈正比例.当然，由于形状不同，比例系数会有所变化，即α＝0.2～0.9.例如上地壳脉形侵入体可取α＝0.8.

当多个源体产生的重力异常叠加时，叠加异常便不再具有类似于（2）式定义的特征尺度.例如，上地壳中心埋深4km的球体源的特征尺度为W＝6.13，下地壳中心埋深34km的球体源的特征尺度W＝52.1；它们叠加之后尺度可在6～52km的范围变化，没有特征尺度.小波多尺度分析方法利用小波基的特征尺度，把地面重力叠加异常场按尺度分解，使产生的异常集重新具有特征尺度.对上述例子，用尺度为6km，12km，24km和48km的小波基分解上述叠加异常，尺度为6km的一阶小波细节D1提取出上地壳中心埋深h＝4km的球体源的异常场，尺度为48km的四阶小波细节提取出下地壳中心埋深h＝34km的球体源的异常场，小波细节便重新恢复了特征尺度.因此，小波多尺度分析方法可用于刻画地壳三维密度结构.

设均匀方格网重力异常场数据的网眼间距为Δ，进行常规小波多尺度分析时，小波基取4～5个样点，它有峰有谷，特征尺度约为峰宽之半，取为Δ／2，小波细节的尺度呈2的整数幂递增.记小波细节阶为n，n阶小波细节特征尺度为

重力异常场小波多尺度分析要求小波细节特征尺度Ln和异常的特征尺度W匹配，即Ln＝W.代入（2）式有

代入（3）式有

（5）式表示重力异常源体埋深h与小波细节阶次n的关系式，称为重力场小波变换的尺度－源深度转换律.给出取样点间距Δ，不同比例尺的地面重力异常场数据在小波变换后由（3）式估算的小波细节特征尺度见表1.由此可见，取样点间距大于20km的重力异常场数据因小波细节特征尺度过大，可能不适合用于刻画精细地壳三维密度结构（表1中下划线数字）.取样点间距小于1km的重力异常场数据因小波细节特征尺度过小，只适合用于刻画上地壳三维密度结构.

表1 小波细节特征尺度Table 1 Characterized scale of wavelet details

综上所述，对于适当比例尺的区域重力数据而言，密度扰动源埋藏深度和地面重力异常场的水平尺度呈近似正比关系，这就是重力场小波变换的尺度－源深度转换律.因此，进行地面重力异常场多尺度分解以后，小尺度的小波细节反映埋藏深度小的场源分布，而大尺度的小波细节反映埋藏深度大的场源分布.自20世纪90年代侯遵泽和杨文采提出并采用小波多尺度分析方法进行重力场多尺度分解以来，重力异常的多尺度分析方法技术取得了良好的效果（侯遵泽和杨文采，1997，2011；侯遵泽等，1998；杨文采等，2001，2005，2014）.究其原因，正是由于小波细节反映不同埋藏深度的场源分布，可用于揭示地壳三维密度扰动的缘故.

3 重力场小波细节的等效层深度公式

区域重力场可表为二维单值实函数Δg（x，y），多尺度分解以后，小波细节为三维单值实函数ΔD（x，y，n）的取样数据集，而地壳三维密度结构模型表为密度扰动源埋藏深度的函数.必需附加约束性规律，把小波细节ΔD（x，y，n）转换成ΔD（x，y，h）的取样数据集，取得对应不同埋藏深度场源等效层的重力异常子集.在最理想的情况下，我们希望分解为对应上、中、下地壳和上地幔顶部起伏的四个重力异常子集.这时要附加布格重力场小波细节的尺度－深度转换约束，位场频率域解释理论为这个问题的解答提供了基础.对位场而言，“频率域”与“波数域”是同义的，指的都是异常随距离变化的频度.

位场频率域解释理论证明（Bhimasankaram et al.，1977；杨文采等，1978，1979；杨文采，1985，1986a），对单个重力场源，重力异常对数振幅谱的斜率与源体顶面或重心埋藏深度成反比.以水平圆柱体为例，对数功率谱为LnS（ω）＝－Ln（2π2Gr2ρ）ωh，其中G为万有引力常数，ρ为源体与围岩密度差，ω为圆频率（波数）；r为半径，h为中心埋藏深度，对数功率谱与频率（波数）成反比，比例系数β＝Ln（2π2Gr2ρ）.对区域众多重力场的场源，区域重力异常对数功率谱的斜率与源体重心平均埋藏深度成反比.因此，如果地壳中重力场的场源集中在几个不同的深度段，区域重力场对数功率谱对应不同频段将呈现有可分辨的几个不同斜率的线段，埋藏深度越大的源体对应斜率越陡的线段.换言之，如果区域重力场对数功率谱对应有可分辨的几个不同斜率的线段，则表明它可以按场源深度分解为几个重力异常子集，它们就是区域重力场小波细节或细节组合，对应不同埋藏深度的场源等效层.因此，布格重力场小波细节的等效层深度公式可近似表为

式中hn为n阶小波细节对应的场源等效层埋藏深度，为功率谱最大斜率；β为比例系数，随源体形状有一定变化.

小波细节合并之后得到了N个异常子集，要把它们分别作二维Fourier变换，求它们的平均功率谱，即各个方向功率谱曲线的平均.由公式（6）可以计算出异常子集对应的等效层平均埋藏深度hn，然后把小波细节集ΔD（x，y，n）转换成ΔD（x，y，h）的取样数据集.必需指出，区域重力场能否按场源深度分解，取决于重力场的频率域属性，并非所有地区都可以作有意义的分解.笔者多年来的研究表明，中国大陆的区域重力场可以按场源深度分解为多个异常子集，详见图2.图2a—2d分别给出了中国大陆、青藏区、扬子区和华北区的区域重力场对数功率谱，区域的范围和数据比例尺见表2.由图2可见，四个区域重力场对数功率谱不同频段呈现有可分辨的几个不同斜率的线段，说明四个地区重力场都可以作有意义的分解.

区域重力场的频率分析得出两个结果：一是重力场可否用作有意义的分解，二是如果可以重力场应分解为几个等效层的异常子集.图2中划出了中国大陆、青藏区、扬子区和华北区重力场对数功率谱不同斜率的线段，除青藏区地壳构造比较复杂可分为5～6层外，其他试验区确定都可分为4～5层.区域重力场范围和数据长度都有限，尤其在高频段，仅进行Fourier变换并不能取得准确的功率谱，有限离散的区域位场功率谱的计算需要特殊的技巧（杨文采，1985，1986a）.注意在图2b—2c中青藏区和扬子区重力场对数功率谱中选择不同斜率的线段时，躲开了最高频段，因为这里含有明显噪声，不符合频率增高重力场对数功率谱斜率明显减少的规律.

表2 中国试验区的范围和数据比例尺Table 2 Number of equivalent layers in different test area

图2 区域重力场对数功率谱（a）中国大陆；（b）青藏区；（c）扬子区；（d）华北区.Fig.2 Logarithmic power spectra in some test regions（a）Chinese continent；（b）Qinghai－Tibetan area；（c）Yangtze craton area；（d）Northern China.

进行重力异常多尺度分析要注意正确选择小波变换最高阶次m.区域重力场分解时选择小波变换最高阶次m的原则是：m必须大于根据功率谱特征决定的等效层数N，等效层数见表2最右列.以青藏区为例，根据布格重力异常的功率谱特征决定重力异常分解成5个层次，根据m＞N的原则，可将m选为8.然后，探讨将8阶小波变换多尺度分析应用于重力异常的分解，可计算出小波细节D1—D8和8阶小波逼近S8.接下来的问题是如何组合m个小波细节形成对应N个等效层场源的重力异常子集.根据小波变换多尺度分析理论可知，小波细节的尺度是随2m次方增加的.因此，小波细节D2的尺度为D1的2倍，小波细节D3的尺度为D1的4倍.由于小波细节D1和D2的尺度相差不大，通常可把（D1＋D2）作为地壳最浅等效层的重力异常对待；它可能代表盆地上部沉积层，也可能代表上地壳结晶基底，还要通过计算场源等效层的平均深度来决定.另一方面，如果某个小波细节的异常幅度很小，也要通过与邻近小波细节合并使等效层场源的重力异常子集相对均衡.

图3 （a）青藏区布格重力异常图；（b）最深等效层异常S8，等效层深度74.8kmFig.3 （a）The Bourger gravity anomalies in the Qinghai－Tibet area；（b）The wavelet approximation S8for the deepest equivalent layer

图4 （a）最浅等效层（D1＋D2）异常图，等效层深度3.03km；（b）最浅等效层脊形化系数图Fig.4 （a）The wavelet detail（D1＋D2）of the Bourger gravity anomalies；refers to the equivalent layer of average depth 3.03km；（b）The corresponding ridge coefficient image of the layer

可以青藏区为例，根据布格重力异常的功率谱特征决定重力异常分解成5个层次，而m＝8，计算出小波细节D1—D8和S8.小波细节合并为（D1＋D2）、（D3＋D4）、（D5＋D6）、（D7＋D8）和S8五个异常子集.青藏高原主要数据来自地面实测，高原外围的部分数据来自卫星观测，数据融合后的青藏区布格重力场见图3a.小波多尺度分解和细节合并之后得到各等效层的重力异常见图4a—7a.由平均功率谱可计算出各层的场源平均埋藏深度，最浅等效层异常（D1＋D2）为3.03km，浅等效层异常（D3＋D4）为12.83km，中等效层异常（D5＋D6）为19.52km；深等效层异常（D7＋D8）为37.93km，最深等效层异常S8为74.8km.最浅等效层、浅等效层和中等效层位于青藏区上地壳顶部、中部和底部；深等效层位于青藏区中地壳底部.最深等效层位于青藏区上地幔顶部，可能主要反映莫霍面起伏，由于本文只研究地壳内部密度成像不予讨论.如果要专门研究下地壳密度扰动，还可以单独研究异常子集D8，或者把D7和D8分别当成独立的等效层.这时平均埋藏深度计算结果为D7对应34km，D8对应52km.

图5 （a）浅等效层（D3＋D4）异常图，等效层深度12.83km；（b）浅等效层脊形化系数图Fig.5 （a）The 3rd and 4th order wavelet details（D3＋D4）refers to the equivalent layer of average depth 12.83km；（b）The corresponding ridge coefficient image of the layer

图6 （a）中等效层（D5＋D6）异常图，等效层深度19.52km；（b）中等效层脊形化系数图Fig.6 （a）The wavelet details（D5＋D6）refers to the equivalent layer of average depth 19.52km；（b）The corresponding ridge coefficient image of the layer

4 等效层的密度扰动反演

在完成重力场小波变换多尺度分解之后，要由小波细节转换为地壳三维密度结构模型，用广义线性反演方法求取各等效层的密度扰动（Yang and Du，1993；杨文采，1986b，1997，1993a，1993b，1995）.位场正演问题是求下述Poisson方程边值问题的解：

其中v（x，y，z＝0）为地表位场值，u（x，y，z）为下半空间任意点位场函数；ρ（x，y，z）为场源项.由数学物理方法知，在地表（7）式的解为（杨文采，1997）

在本文中v用小波细节ΔD（x，y，n＝1，2，3，…）代入；ρ用等效层体密度扰动Δρ＝ρ－ρ′代入；（8）式可用广义线性反演求解（杨文采，1986b）.

等效层的密度扰动指平面上变化的场源密度与等效层平均密度的差.除小波多尺度分析的小波细节数据外，广义线性反演还要求输入各等效层的平均密度ρ′和场源平均埋藏深度z′.各等效层的平均密度主要根据岩石物理性质测定结果给出.以青藏区为例，如最浅等效层平均密度取为2.60g·cm－3，上地壳中部浅等效层密度取为2.67g·cm－3，上地壳底部中等效层密度2.81g·cm－3；中地壳底部深等效层密度2.95g·cm－3；上地幔顶部最深等效层密度3.20g·cm－3.由于反演出的各等效层的横向相对密度扰动图外观与小波细节异常（图4a—7a）差别不大，本文不再罗列.

5 重力场多尺度刻痕分析

在反演计算了重力多尺度密度扰动图件之后，我们就获得了地壳三维密度扰动的映像.这些映像包含有大量地壳构造的信息，能否根据现代信息学理论用计算机把地壳构造的信息直接提取出来？要回答这个问题，首先要问你想提取什么地壳构造信息？我们首先想到的是现存的地壳变形带分布信息.

从物理学的观点看，地壳变形带是动力地质作用在地壳中的刻痕，它造成地壳岩石的密度在局部呈条带状变化.因此，不同强度的地壳变形带是岩石圈地质作用对应地壳中长条形的密度异常带，它们又在重力场中造成“刻痕”.这种刻痕的特征参数包括：重力异常梯度急剧变化，异常在平面上各向异性较强，各向异性长轴走向比较稳定，等等.把隐含在区域重力场中有关刻痕的所有特征参数都提取出来，可在综合分析的基础上识别不同时期形成的地壳变形带.计算了重力等效层多尺度密度扰动之后，从其中提取地壳变形带信息要有新的数学工具，这就是表面刻痕分析（孙艳云和杨文采，2014）.面对尺度变化巨大的自然界，应用数学中的局部分析和全局分析永远是准确认识自然的重要工具之一（Riley，1974）.表面刻痕分析就是局部分析的一个分支.

区域位场函数g＝g（x，y）是一个二维连续可微函数，可视为表面函数的一种特殊类型.把区域重力场网格化成一个个表面元，每个表面元的信息可用表面元的谱矩表征（Longust－Higgins，1962；Nayak，1973；Yanagi et al.，2001；黄逸云，1984，1985；杨叔子等，2007）.最简单的谱矩为二阶谱矩，它的三个元m20、m02和m11有各自的物理含义.m20、m02分别表示位场表面x和y方向的斜率的方差；m11是x和y方向的斜率的协方差.位场轮廓是指一个切平面与位场表面的交线.若已知表面的二阶谱矩为m20、m02和m11，位场任意表面元的一条轮廓在θ（与x轴方向的夹角）方向上的二阶谱矩m2（θ）由（9）式唯一确定

m2（θ）＝m20cos2θ＋2m11cosθsinθ＋m02sin2θ，（9）m2（θ）等于θ方向上轮廓的斜率的方差.若m2max和m2min是m2（θ）的最大值和最小值，则其所在的方向总是相互垂直的，这两个方向称为刻痕主方向.m2min的那条轮廓曲线与x轴的夹角称为主方向角θp，表示表面刻痕的走向

位场表面谱矩依赖于坐标系，如果坐标系旋转，尽管此时表面性质不变，表面的谱矩将会改变.因此，需要考虑一种与坐标系旋转无关的统计不变量.位场表面的二阶谱矩可以用统计不变量M2和Δ2来表达，它们是（Longust－Higgins，1962；Nayak，1973；Yanagi et al.，2001；黄逸云，1984，1985；杨叔子等，2007）

统计不变量M2表示位场表面元的二维斜率的方差，可以表征布格重力异常的刻痕强弱，定义为刻痕的强度系数；试验表明，M2着重表现幅度大的表面刻痕.统计不变量Δ2主要与位场表面元的各向异性有关，表示表面元刻痕脊形化的程度，且总有Δ2≥0.试验表明，脊形化反映的是场在局部区段各向异性的程度，Δ2可以较全面地反映表面刻痕随方向变化的形貌.

图7 （a）深等效层（D7＋D8）异常图，等效层深度37.93km；（b）深等效层脊形化系数图L—低密度扰动中心；H—高密度扰动中心.Fig.7 （a）The wavelet details（D7＋D8）refers to the equivalent layer of average depth 37.93km；（b）The corresponding ridge coefficient image of the layerLetter“L”denotes low－density units，and letter“H”denotes high－density units.

图8 多尺度密度扰动边界刻痕系数图像（a）为中等效层，等效层深度为19.52km；（b）为深等效层，等效层深度为37.63km.L—低密度扰动中心；H—高密度扰动中心.Fig.8 The ridge－boundary coefficient image of the studied area（a）The middle layer of depth 19.52km；（b）The deep layer of depth 37.63km.Letter“L”denotes low－density units，and letter“H”denotes high－density units.

把隐含在区域重力场中有关刻痕的特征参数M2和Δ2提取出来后，便可在综合它们特征的基础上识别局部结构为脊索形的地壳变形带.为了对位场表面元的刻痕准确定位，我们综合这两个统计不变量定义刻痕的脊形化系数Λ为

Λ值变化区间为［0，1］，体现了表面形态的各向异性特征；与统计不变量M2成反比，可增加对位场表面弱刻痕的识别能力.脊形化系数Λ反映的是面元内场源密度的各向异性度的变化，Λ＝0表明面元内场源密度各向同性，Λ＝1表明面元内场源密度呈线性高度各向异性.场源密度呈线性高度各向异性就是地壳变形带的典型特征.

提取区域位场函数刻痕信息就是要由地壳三维密度结构模型计算等效层密度扰动的脊形化系数Λ（x，y，n＝1，2，3，…）.孙艳云和杨文采（2014）详细阐述了提取脊形化系数的理论与方法，并将这种方法应用于中国大陆重力场刻痕信息识别，对比了不同类型的重力场刻痕与已知的地壳变形带的对应关系.与之前不同的是，本文不是用区域重力场，而是用区域重力场多尺度分解之后取得的异常子集数据，或者多个等效层密度扰动数据作为刻痕信息识别的输入数据，所以叫多尺度刻痕分析.以青藏区为例，把图4a—7a中布格重力等效层或等效层密度扰动作为输入，求得的脊形化系数Λ分别见图4b—7b.

6 边界刻痕系数和构造边界划分

迄今为止，大陆构造单元的划分都是大地构造学家各自进行的，不同大地构造学家划分的大陆构造单元各自不同，缺乏客观的标准.如果能够用数学方法把构造单元信息从地球物理场中识别和提取出来，用计算机自动识别构造单元就成为可能.在提取区域位场函数的刻痕信息，即计算等效层密度扰动的脊形化系数Λ之后，我们还可以进一步提取密度扰动的刻痕边界信息，因为它们也是大陆构造单元的边界.

在提取构造单元边界信息时，边界越窄单元划分的定位越准确.脊形化系数Λ＞0.5的条带常常较宽，为了准确地对边界定位，要再定义一个参数以取得窄的边界线，这个参数称为边界刻痕系数MΛ，它对应脊形化系数Λ较大区域的边界，可以从脊形化系数Λ作进一步增强处理计算出，详见参考文献（Longust－Higgins，1962；孙艳云和杨文采，2014；Riley，1974；黄逸云，1984，1985；杨叔子等，2007）.边界刻痕系数MΛ提取连续性好的脊形化系数分布区段，并将它的边界锐化定位.试验表明，在结晶岩分布区，脊形化系数Λ图像反映显生宙地壳变形带，边界刻痕系数图像反映构造单元边界位置（孙艳云和杨文采，2014）.以青藏区为例，把图6—7中等效层密度扰动和脊形化系数Λ作为输入，求得的中、深等效层边界刻痕系数分别见图8a，8b.

提取等效层密度扰动边界刻痕系数MΛ之后，不仅可以取得自动划分的等效层构造单元边界位置，还可以取得等效层构造单元的密度扰动属性.以青藏区为例，首先研究图7a中的等效层密度扰动极性正负，在图中用字母“L”标明等效层低密度扰动中心位置；用“H”标明等效层高密度扰动中心位置.对比图7a和8a可见，自动划分的等效层构造单元都与标明“L”或“H”的等效层密度扰动中心对应，因此可知等效层构造单元属于高密度地块还是低密度地块.

以上介绍了应用区域重力场刻划地壳三维密度结构和构造单元分区的数据处理、反演解释和信息提取的方法研究的主要结果，青藏区的实际应用仅作结果图示，未对其地质含义进行深入讨论.由于青藏区的重力场多尺度刻痕的地质含义仍然要继续研究，今后再另外拟文发表.

7 结论

（1）把小波多尺度分析、表面刻痕分析以及位场频率域解释理论和反演方法有机地结合起来，可形成一个数据处理、反演解释和信息提取的方法系统，这一方法系统简称为区域重力场多尺度刻痕分析方法，应用于刻画地壳分层的三维密度结构、地壳变形带和构造单元分区.

（2）多尺度刻痕分析包含频率域重力场按场源分层、小波变换多尺度分析、场源分层深度及密度扰动反演、分层刻痕分析和构造边界定位四个子系统.本文扼要地介绍其方法原理及应用效果.

（3）多尺度刻痕分析除应用小波变换和表面谱矩之外，还运用了地球重力场本身的规律，即以前导出的重力场小波细节的等效层深度公式和本文导出的重力场小波变换的深度－尺度转换律.

（4）从地球物理探测到大地构造学发现，是一个多学科综合研究的探索过程.要取得重大研究成果，必须研发和组合来自不同学科的多个新方法技术，使多学科综合研究有宽厚的理论支撑.本文介绍的四个子系统组合的理论支撑分别来自应用数学、地球物理学和信息科学.

Bhimasankaram V L S， Nagendra R，Rao S V S.1977.Interpretation of gravity anomalies due to finite inclined dikes using Fourier transformation.Geophysics，42（1）：51－59.

Daubechies I.1990.The wavelet transform，time frequency localization and signal analysis.IEEETransactiononInformationTheory，36（5）：961－1005.

Hou Z Z，Yang W C.1997.Wavelet transform and multiscale analysis of gravity field in China.ChineseJ.Geophys.（in Chinese），40（1）：85－95.

Hou Z Z，Yang W C，Liu J Q.1998.Multi－scale inversion of density contrast within the crust of China.ChineseJ.Geophys.（in Chinese），41（5）：642－651.

Hou Z Z，Yang W C.2011.Multi－scale inversion of density structure from gravity anomalies in Tarim Basin.Sci.inChina（Ser.D），54（3）：399－409.

Hou Z Z，Yang W C.2012.Applications of Wavelet Multi－Scale Analysis（in Chinese）.Beijing：Science Press.

Huang Y Y.1984.The characterization of three－dimensional Radom surface topography.JournalofZhejiangUniversity（in Chinese），18（2）：138－148.

Huang Y Y.1985.Geometrical interpretation and graphical solution of second order spectrum moments and statistical invariants for random surface characterization.JournalofZhejiangUniversity（in Chinese），19（6）：143－153.

Longust－Higgins M S.1962.The statistical geometry of random surfaces in hydrodynamic instability.／／13th Sympo.Applied Mathem.Amer.Math.Soc.，105－143.

Mallat S G.1989.Multifrequency channel decompositions of images and wavelet models.IEEETransactionsonAcoustics，Speech，andSignalProcessing，37（12）：2091－2110.

Nayak P R.1973.Some aspects of surface roughness measurement.Wear，26（2）：165－174.

Riley K F.1974.Mathematical Methods for the Physical Sciences：An Informal Treatment for Students of Physics and Engineering.New York：Cambridge Univ.Press.

Sun Y Y，Yang W C.2014.Recognizing and extracting the information of crustal deformation belts from the gravity field.ChineseJ.Geophys.（in Chinese），57（5）：1578－1587，doi：10.6038／cjg20140521.

Yanagi K，Hara S，Endoh T.2001.Summit identification of anisotropic surface texture and directionality assessment based on asperity tip geometry.InternationalJournalofMachine ToolsandManufacture，41（13）：1863－1871.

Yang S Z，Wu Y，Xuan J P，et al.2007.Identification on Surface Topography，in Time Series Analysis in Engineering Application （in Chinese）.Wuhan：Huazhong University of Science ＆ Technology Press.

Yang W C，Guo A Y，Xie Y Q.1978.Theory and methods for interpretation of gravity and magnetic anomalies in the frequency domain.BullInst.Geophys.Geochem.Prospect（in Chinese），（2）：134－179.

Yang W C.1979.Theorem of equivalence between the twodimensional and three－dimensional potential fields in the frequency domain.ActaGeophys.Sin.（in Chinese），22（1）：89－96.

Yang W C，Guo A Y，Xie Y Q，et al.1979.Interpretation of gravity anomalies in frequency domain （A）.ComputingTechniquesfor GeophysicalandGeochemicalExploration（in Chinese），1（1）：1－16.Yang W C.1985.The power spectrum analysis of the field data（A）.ComputingTechniquesforGeophysicalandGeochemical Exploration（in Chinese），7（3）：188－196.

Yang W C.1986a.The power spectrum analysis of the field data（B）.ComputingTechniquesforGeophysicalandGeochemical Exploration（in Chinese），8（1）：14－22.

Yang W C.1986b.A generalized inversion technique for potential field data processing.ChineseJ.Geophys.（ActaGeophys.Sin.）（in Chinese），29（3）：283－291.

Yang W C.1993a.Nonlinear chaotic inversion of seismic tracesⅠ：Theory and numerical experiments（Ⅰ）.ChineseJ.Geophys.（in Chinese），36（2）：222－232.

Yang W C.1993b.Nonlinear chaotic inversion of seismic tracesⅡ：Lyapunov exponents and attractors.ChineseJ.Geophys.（in Chinese），36（2）：376－387.

Yang W C，Du J Y.1993.Approaches to Solve Nonlinear Problems of the Seismic Tomography.／／Wei Y，Gu B L eds.Acoustic Imaging Volume 20.New York：Plenum Press，591－603.

Yang W C.1995.BG－Inverse scattering method for inversion of seismic wavefield.ChineseJ.Geophys.（in Chinese），38（3）：358－366.

Yang W C.1997.Theory and Methods in Geophysical Inversion（in Chinese）.Beijing：Geological Publishing House.

Yang W C，Shi Z Q，Hou Z Z，et al.2001.Discrete wavelet transform for multiple decomposition of gravity anomalies.ChineseJ.Geophys.（in Chinese），44（4）：534－541，doi：10.3321／j.issn：0001－5733.2001.04.012.

Yang W C.2009.Tectonophysics in Paleo－Tethys Domain （in Chinese）.Beijing：Petroleum Industry Press.

Yang W C，Wang J L，Zhong H Z，et al.2012.Analysis of regional magnetic field and source structure in Tarim Basin.ChineseJ.Geophys.（in Chinese），55（4）：1278－1287，doi：10.6038／j.issn.0001－5733.2012.04.023.

Zeng H L.2005.Gravity Field and Gravity Exploration （in Chinese）.Beijing：Geological Publishing House.

附中文参考文献

侯遵泽，杨文采.1997.中国重力异常的小波变换与多尺度分析.地球物理学报，40（1）：85－95.

侯遵泽，杨文采，刘家琦.1998.中国大陆地壳密度差异多尺度反演.地球物理学报，41（5）：642－651.

侯遵泽，杨文采.2011.塔里木盆地多尺度重力场反演与密度结构.中国科学（D辑），41（1）：29－39.

侯遵泽，杨文采.2012.小波多尺度分析应用.北京：科学出版社.

黄逸云.1984.三维随机表面形貌的识别.浙江大学学报，18（2）：138－148.

黄逸云.1985.随机表面二阶谱矩和统计不变量的集合解释及其图解法.浙江大学学报，19（6）：143－153.

孙艳云，杨文采.2014.从重力场识别与提取地壳变形带信息的方法.地球物理学报，57（5）：1578－1587，doi：10.6038／cjg20140521.

杨叔子，吴雅，轩建平等.2007.时间序列分析的工程应用（下）.武汉：华中科技大学出版社.

杨文采，郭爱缨，谢玉清.1978.重磁异常频率域解释的理论与方法.物化探研究报导，（2）：134－179.

杨文采.1979.二维与三维位场在频率域的等价定理.地球物理学报，22（1）：89－96.