优化滤波方法及其在中国大陆布格重力异常数据处理中的应用

郭良辉，孟小红，石 磊，陈召曦

1 地下信息探测技术与仪器教育部重点实验室（中国地质大学，北京），北京 100083

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

3 中国地震局地球物理研究所，北京 100081

1 引 言

根据重力观测值整理得到的重力异常，包含了从地表到深部所有密度不均匀体引起的重力效应，信息非常丰富.实测重力异常是所有这些重力效应的叠加.要根据实测重力异常反演解释某个地质体（比如岩矿体、地质构造或物性分界面），必须首先从叠加异常中分离出单纯由这个地质体引起的异常，然后用这个异常进行反演解释.然而，从叠加异常中分离出某个地质体引起的异常，或者把叠加异常分解为几个地质体引起的单一异常，是比较困难的.重力异常分离是目前重力资料处理解释中没有很好解决的难题之一.

当前，重力异常分离方法众多，如最小二乘平滑法、带通滤波法、向上延拓法、匹配滤波法［1］、补偿圆滑滤波法［2］、正则化滤波法［3］、维纳滤波法［4-5］、小波变换法［6-8］、熵滤波法［9-10］、非线性滤波法［11］等.这些方法均有各自的特点和优势，但同时又有各自的局限性.

传统向上延拓法存在的一个主要问题是在进行向上延拓时，浅源短波长信息和深源长波长信息同时延拓，在浅源信息被压制的同时，深源信息也受到压制，将实测异常减去上延后的区域异常所得到的剩余异常中，仍包含有区域成分，异常分离不彻底［12-13］.针对该问题，Pawlowski［12］根据维纳滤波和格林等效层原理提出了优选延拓法，应用该方法进行向上延拓时，可以在压制浅源短波长信息的同时，保持深源长波长信息不衰减.此后，许德树等［14］、孟小红等［15］在优选延拓法基础上给出了基于优选向上延拓的延拓差值法，实现对某一指定频段或尺度的局部异常有选择地提取或分离，即类似于带通滤波.传统向上延拓法存在的另一个重要问题是延拓高度必须已知，然而实际应用中延拓高度往往是未知的［13］，优选延拓法同样存在此问题.针对这个问题，Zeng等［13］提出了向上延拓法的最佳延拓高度估计方法，该方法分别以一系列递增的延拓高度对实测异常做向上延拓，得到相应的延拓场，计算各前后不同延拓高度的延拓场的互相关，这些互相关值的最大偏差点所对应的高度即为最佳延拓高度.郭良辉等［16］和 Meng等［15］给出了优选向上延拓的最佳延拓高度估计办法.

由于延拓高度往往是未知的，这使得传统向上延拓法和优选延拓法的实际应用受到一定限制.本文在优选延拓法的理论基础上，进一步研究提出基于维纳滤波和格林等效层原理的优化滤波法，该方法分离重力异常不需要已知延拓高度，可对重力异常进行指定频段的低通、带通和高通滤波.

本文首先给出优化滤波法的方法原理和计算步骤，然后利用该方法对理论模型的含噪重力异常进行去噪和分离异常试验，并与传统向上延拓法和带通滤波法的处理结果进行对比，最后利用优化滤波法对中国大陆重力异常数据进行去噪和异常分离，并对分离出的区域重力异常做密度界面约束反演，得到中国大陆莫霍面深度分布特征.

2 优化滤波方法

2.1 优化滤波算子

假设实测重力异常gobs（x，y）是深源长波长异常gd（x，y）与浅源短波长异常gs（x，y）的叠加：

那么实测重力异常的功率谱密度函数P可看成是深源异常功率谱密度函数Pd与浅源异常功率谱密度函数Ps之和［5］，

Pawlowski［5］指出实测重力异常的傅里叶功率谱可通过引入格林等效层概念来模拟，即利用处于不同深度的m＋n层格林等效源薄层来建立实测重力异常的径向平均对数功率谱模型（假定各等效层互不相关），

其中，深源场成分用m层格林等效层来模拟，而浅源场分量用n层格林等效层来模拟，误差E来源于对实测数据的功率谱的不完全拟合.显然，对于不同地区，充分拟合时等效层的数目将有所不同.

按照 Naidu［17］和 Dampney［18］的定义，任意深度h处格林等效层的功率谱密度函数可写成

其中，s是与期望的等效层强度成正比的常数，k是角频率.

设目标层为深部i→j（n＋m≥j≥i≥1）等效层，则优化滤波法的期望是保持目标层场源信息不变，而对其它层场源信息做压制，即，

设维纳滤波器的输入是实测异常gobs（x，y），期望输出是gpf（x，y），按照维纳滤波器设计的要求，滤波器的频率响应［19-20］为

其中，G代表频率域的频谱，〈〉代表数学期望，＊代表复共轭.式中的分子是滤波器期望输出和输入信号（实测异常）的互功率频谱密度函数，分母是滤波器输入的功率频谱密度函数.

结合公式（3），式（6）的维纳滤波转换函数可改写为

其中，P（i→j）为目标层场源信息的功率谱密度函数，P′（i→j）是目标层场源信息与其它层场源信息的互功率谱密度函数.

假设目标层与其它深度层场源信息互不相关，则P′（i→j）＝0，那么公式（7）的维纳滤波转换函数可进一步改写为

称公式（8）的维纳滤波转换函数为优化滤波算子.当i＝1时，i→j即为1→j，此时的优化滤波是分离1→j层场源信息的低通滤波；当j＝m＋n时，i→j即为i→m＋n，此时的优化滤波是分离i→m＋n层场源信息的高通滤波；当n＋m＞j≥i＞1时，此时的优化滤波是分离i→j层场源信息的带通滤波.因此，利用公式（8）可对实测重力异常进行不同形式的优化滤波，分离出目标层场源在原始观测面上的异常信息.

2.2 优化滤波计算步骤

优化滤波法的计算步骤如下：

（1）对实测异常做傅里叶变换，计算径向平均对数功率谱；

（2）用分段线性拟合法拟合实测径向平均对数功率谱；

首先分析实测径向平均对数功率谱形状特征，确定格林等效源层数和各等效层分布的频率范围.一般地，径向频率越低，频段分割越细窄，段数越多，而径向频率越高，频段分割越粗糙，段数越少.每频段对应着一个格林等效层，低频段主要对应区域深源等效层，中频段主要对应浅源等效层，而高频段主要对应噪声［1，4-5］.然后利用直线拟合实测径向平均对数功率谱的每个分段功率谱.最后利用（9）式估计各等效层的近似深度［1，21-22］：

其中，r1和r2分别是拟合直线的径向频率起点和终点，P（r1）和P（r2）分别是r1和r2对应的实测径向平均对数功率谱.

（3）利用公式（4）表示的功率谱密度函数模型和所估计出的各等效层深度值来拟合实测径向平均对数功率谱；

设a＝e－2kh，公式（4）对于第i个径向频率点可写成

其中，Pi是第i个径向频率点的径向平均对数功率谱，ei是功率谱拟合误差.

对于l个功率谱样点，有l线性方程构成的方程组，写成矩阵如下：

通过共轭梯度法等迭代优化算法对公式（12）求解可得到向量S，把它代入公式（4）可计算任意第i个径向频率点或第i个格林等效层的理论径向平均对数功率谱.

（4）利用公式（8）构建优化滤波算子；

（5）对实测异常进行优化滤波处理；

（6）进行反傅里叶变换，得到空间域的目标异常.

3 理论模型数据试验

图1 （a）理论模型各长方体在水平面位置分布，（b）理论重力异常等值线图，（c）含噪重力异常等值线图.重力异常单位：mGal.Fig.1 The outline of the model on the horizontal plane（a），the theoretical gravity anomaly of the model（b），and the noised gravity anomaly（c）.The unit of gravity anomaly is mGal.

理论模型由三个不同深度层、不同大小和剩余密度值的11个直立长方体组合而成，图1a显示了各直立长方体在水平面上的投影分布.假设深层长方体A1、A2和A3产生区域重力异常，中层长方体B1、B2、B3、B4、B5和浅层长方体C1、C2和C3共同产生局部重力异常（或剩余重力异常），理论模型的理论重力异常等值线图见图1b，对理论重力异常加入基准值4%的高斯随机噪声见图1c.

对含噪重力异常计算径向对数功率谱（图2），通过功率谱形状分析，用0～0.2745cycles·km-1（频段1）、0.2745～1.2157cycles·km-1（频段2）和1.2157～5cycles·km-1（频段3）三个等效层分段拟合（图2）.这里，频段1的功率谱主要对应A层长方体的低频异常信息，频段2的功率谱主要对应B和C层长方体的中高频信息，而频段3的功率谱主要对应高斯噪声的高频干扰信息.

下面对含噪重力异常进行优化滤波去噪和分离异常试验.首先进行频段1和2的优化滤波（低通滤波），目的是压制高斯噪声，优化滤波算子的频率响应见图3a，滤波结果见图3d.从图可见，优化滤波法有效压制了高斯噪声干扰，除C层有效信号保留的不是很好之外（因为其频率与噪声的较接近），A和B层的有效信号都较好地保留了.然后进行频段1的优化滤波（低通滤波），目的是分离出A层各长方体引起的区域异常，优化滤波算子的频率响应见图3b，滤波结果见图3e.从图可见，优化滤波法有效压制了高斯噪声干扰及B和C层的局部异常，A层的有效信号被较好地保留.最后进行频段2的优化滤波（带通滤波），目的是分离出B和C层各长方体共同引起的局部异常，优化滤波算子的频率响应见图3c，滤波结果见图3f.从图可见，优化滤波法有效压制了高斯噪声干扰及A层的区域异常，B和C层的有效信号得到较好的保留.

图2 含噪重力异常径向对数功率谱及其分段拟合黑色：含噪重力异常功率谱；蓝色：分频段拟合；红色：模型拟合的功率谱.Fig.2 The radially averaged logarithm power spectrum of the noised gravity anomaly（black line）and its fitting by using piece－wise linearizationThe blue lines show the fitted straight lines in each radial frequency section，and the red line indicates the fitted radial logarithm power spectrum by using power spectrum density function models.

图3 压制高斯噪声的优化滤波算子频率响应（a）和滤波结果（d），分离A层长方体异常的优化滤波算子频率响应（b）和滤波结果（e），分离B和C层长方体异常的优化滤波算子频率响应（c）和滤波结果（f）.重力异常单位：mGal.Fig.3 The amplitude response of the preferential filtering operator for suppressing Gaussian noise （a），separating regional anomalies of layer A （b），and separating local anomalies of layers B and C（c）.（d），（e）and（f）show the corresponding anomalies of（a），（b）and（c）.The unit of gravity anomaly is mGal.

图4 带通滤波法的压制高斯噪声（a）、A层的区域异常（b）、B和C层的局部异常（c），传统向上延拓法的压制高斯噪声（d）、区域异常（e）和局部异常（f）.重力异常单位：mGal.Fig.4 The denoised anomalies（a），the separated regional anomalies（b）and the separated local anomalies（c）by using the low－pass filtering，and those（d，e，f）of the upward continuation.The unit of gravity anomaly is mGal.

下面，分别利用带通滤波法和传统向上延拓法对理论含噪重力异常进行同样的去噪和异常分离试验.压制高斯噪声时，带通滤波法的滤波波长宽度为800m，向上延拓高度为300m，结果分别见图4（a，d）.从图可见，带通滤波法有效压制了高斯噪声干扰，保留有效信号，而传统向上延拓法在压制高斯噪声干扰的同时，也部分压制了有效信号.分离A层的区域异常时，带通滤波法的滤波波长宽度为4000m，向上延拓高度为1000m，结果分别见图4b和（e）.从图可见，带通滤波法和传统向上延拓法在有效压制高频噪声干扰及B和C层局部异常的同时，也部分压制了A层的区域异常，尤其是传统向上延拓法的区域异常被压制较多.分别用两种方法压制高斯噪声后的异常减去区域异常，即可得到分离B和C层的局部异常，结果分别见图4（c，f）.从图可见，带通滤波法、传统向上延拓法都能分离出B和C层的局部异常，但也遗留了A层的部分区域异常，尤其是传统向上延拓法的区域异常遗留较多，也就是说异常分离不彻底.

图5显示了这三种不同方法分离出的异常与理论重力异常（图1b，没有噪声干扰）沿X＝5680m剖面的对比图.从图4和图5可见，在压制高斯噪声方面，三种方法分离效果基本一致，但传统向上延拓法部分压制了有效信号；在分离A层区域异常方面，优化滤波法较好保留了区域异常，带通滤波法部分压制了区域异常，而传统向上延拓法较多压制了区域异常；在分离B和C层局部异常方面，优化滤波法结果最接近理论值，而带通滤波法部分遗留了区域异常，传统向上延拓法较多遗留了区域异常.

图5 不同方法在X＝5680m剖面的异常分离结果与理论重力异常对比（a）压制高斯噪声；（b）A层的区域异常；（c）B和C层的局部异常.蓝色：理论值；红色：优化滤波法；绿色：带通滤波法；黄色：传统向上延拓法.Fig.5 Comparisons between the theoretical gravity anomaly（blue lines）and the separated anomalies by using the three methods along profile X＝5680m（a）Denoised anomalies；（b）Separated regional anomalies；（c）Separated local anomalies.Red lines correspond to the preferential filtering，green lines to the low－pass filtering，and yellow lines to the upward continuation.

4 中国大陆重力异常数据试验

依托“深部探测技术与试验研究”项目—第五子课题（Sinoprobe－01－05），开展了中国大陆重力异常数据的应用实验研究.利用最新的地球重力场模型EGM2008（Earth Gravitational Model 2008）［23］计算了中国大陆自由空气重力异常.EGM2008是由美国国家地理空间情报局（NGA）推出的官方版高阶地球重力场模型，该模型的阶次完全至2159（另外球谐系数的阶扩展至2190，次为2159），空间分辨率约为5′.由于采用了GRACE卫星跟踪数据、卫星测高数据及地面重力数据等，使得该模型无论在精度还是在分辨率方面均取得了巨大的进步.该模型在海域的精度可达到1∶100万比例尺重力勘探的要求［23］，在中国大陆大部分地区精度可达10mGal［24］.因此，EGM2008可用于中等比例尺的海洋地质与资源调查和小比例尺的大陆重力编图与构造研究.

利用EGM2008，计算了中国大陆自由空气重力异常（网度为10km），然后按照区域重力调查规范做地形校正和布格校正（中间层校正），最终得到中国大陆布格重力异常数据（网度为10km），见图6所示.校正所用的地形高程数据（图10）来自于美国斯克里普斯海洋研究所网站（http：／／topex.ucsd.edu／cgi－bin／get＿data.cgi），空间分辨率为1′.从图6可见，校正后得到的布格重力异常在中国大陆范围内普遍存在高频噪声干扰，经分析认为，这些噪声干扰一部分是由EGM2008计算的自由空气重力异常数据自身携带而遗留下来的，另一部分来源于地形校正和布格校正，而后者主要是由于自由空气重力异常数据和地形数据的分辨率与精度不一致造成的.显然，高频噪声干扰对后续数据处理解释及应用等将产生严重影响，应首先对重力异常数据进行去噪处理.

对图6所示的中国大陆布格重力异常数据尝试采用了传统向上延拓法、带通滤波法、巴特沃斯滤波法等多种常规方法进行去噪处理，但结果都不理想，要么噪声没有压制干净，要么有效信号部分损失，利用本文的优化滤波法对中国大陆布格重力异常数据进行去噪试验，取得较满意结果.

图6 EGM2008计算并经过校正得到的中国大陆布格重力异常图（网度为10km）Fig.6 The Bouguer gravity anomaly of Chinese continent calculated from EGM 2008 with terrain correction and Bouguer correction.The grid size is 10km.

首先计算中国大陆布格重力异常的径向平均对数功率谱（图7），根据功率谱形状特征，用0～0.0014cycles·km-1（频段 1）、0.0014～0.0039 cycles·km-1（频段2）、0.0039～0.0099cycles·km-1（频段3）、0.0099～0.0174cycles·km-1（频段4）、0.0174～0.0249cycles·km-1（频段5）、0.0249～0.05cycles·km-1（频段6）6个等效层分段拟合.分析认为，频段5和频段6主要对应高频噪声干扰，频段1～4主要对应有效信号.然后进行频段1～4的优化滤波（低通滤波），目的是压制高频噪声干扰，得到的去噪布格重力异常见图8所示.从图可见，优化滤波法有效压制了噪声干扰，又较好地保留了深部有效信息.通过对比，去噪后的布格重力异常与前人成果［25］基本一致.之后进行频段1～2的优化滤波（低通滤波），目的是分离出中国大陆区域重力异常，结果见图9所示.

从图8和9可见，中国境内布格重力异常值自东向西逐渐减小.布格重力异常值在海岸线为0mGal左右；由海岸线向西，异常值缓慢递减，并进入负值区；沿大兴安岭—太行山—武陵山一带显示为过渡带，这是我国一条重要的重力梯级带，走向为NE向；再往西，至青藏高原周边地区，异常值迅速减小，这是另一条重要的重力梯级带，围绕青藏高原呈弧形展布；青藏高原大部分地区的异常值小于-400mGal.可见，由优化滤波法去噪和分离出的中国大陆布格重力异常，能较好反映出全国重力场变化特征，提供关于地壳和上地幔中物质质量分布的信息.其中，区域布格重力异常（图9）主要由莫霍面引起的，反映了莫霍面起伏的基本特征.

图7 中国大陆布格重力异常径向对数功率谱及其分段拟合黑色：中国大陆布格重力异常功率谱；蓝色：分频段拟合；红色：拟合的功率谱.Fig.7 The radially averaged logarithm power spectrum of the Bouguer gravity anomaly of Chinese continent（black line）and its fitting by using piece－wise linearization The blue lines show the fitted straight lines in each radial frequency section，and the red line indicates the fitted radial logarithm power spectrum by using the power spectrum density function models.

图8 优化滤波去噪后的中国大陆布格重力异常（网度为10km）Fig.8 The denoised Bouguer gravity anomaly of Chinese continent by the preferential filtering method.The grid size is 10km

图9 优化滤波分离出的中国大陆区域重力异常（网度为10km）Fig.9 The separated regional gravity anomaly of Chinese continent by the preferential filtering method.The grid size is 10km

为认识中国大陆莫霍面分布特征，利用优化滤波法分离出的区域重力异常（图9）反演了中国大陆莫霍面深度，其中反演方法采用带先验信息约束的非线性回归法［26］，以减少反演的多解性.自20世纪60年代以来我国开展了大量的深地震测深、深地震反射和宽频地震观测等地震探测，其具有较大的探测深度和较高的垂向分辨率，为研究中国大陆深部结构和构造格架提供了及其重要的信息.对这些人工地震探测推断的莫霍面深度图进行了收集、整理和矢量化，最终得到中国大陆莫霍面深度控制点数据库［27］，数据点分布见图10所示，它们将作为本文莫霍面深度约束反演的已知控制点信息.

以图10所示的已知控制点信息作为约束，采用非线性回归法对优化滤波法分离出的区域重力异常（图9）进行莫霍面深度约束反演，得到中国大陆莫霍面深度分布，见图11所示.反演结果的均方差为4.58km.通过对比分析，本文反演得到的莫霍面深度分布与前人成果［28］基本一致，但细节稍多些.从图11可见，以东经105°～110°的南北地震带为界，东西两侧的中国大陆莫霍面深度明显不同，以东地区相对较浅，介于30～45km，莫霍面等值线主要呈北东向展布，以西地区相对较深，在45～74km之间，莫霍面等值线主要呈东西向展布，体现出东西不同的深部结构和构造动力学背景.由于篇幅限制，重力约束反演得到的中国大陆莫霍面深度特征（图11）将另文详细描述.

图10 中国大陆地形图及已知莫霍面深度控制点分布（红点）Fig.10 The topography of Chinese continent and positions of the Moho depth control points deduced by the deep seismic detection（red points）

图11 重力反演得到的中国大陆莫霍面深度分布（网度为10km）Fig.11 The Moho depth distribution of Chinese continent from gravity inversion.The grid size is 10km

5 结 论

本文在优选延拓法的理论基础上，研究提出了基于格林等效层概念和维纳滤波器的优化滤波法，用于对重力异常去噪和分离，给出了优化滤波算子的构建和优化滤波的计算步骤.理论模型数据试验表明本文方法有效，异常分离效果优于传统的带通滤波法和向上延拓法.通过中国大陆重力异常数据的优化滤波去噪和异常分离，得到有效的布格重力异常数据和区域重力异常数据，对后者进行密度界面约束反演得到中国大陆莫霍面深度分布.

与优选延拓法相比，优化滤波法与延拓高度无关，其异常分离不受延拓高度困扰，具有一定的优势.当然，优化滤波法也有其缺点，如利用格林等效层近似地下场源，且假设目标层与其他深度层场源信息互不相关，而实际地质中，不同深度层场源可能存在部分相关，比如浅层场源有可能向下延伸到深部；另外，不同深度层场源的频谱都是宽频带的，且互相重叠，这样，很难从叠加的频谱中彻底分离出目标层场源的单一频谱.

（References）

［1］Spector A，Grant F S.Statistical models for interpreting aeromagnetic data.Geophysics，1970，35（2）：293－302.

［2］侯重初.补偿圆滑滤波方法.石油物探，1981，（2）：22－29.Hou Z C.Filtering of smooth compensation.Geophysical Prospecting for Petroleum （in Chinese），1981，（2）：22－29.

［3］安玉林，管志宁.滤除高频干扰的正则化稳定因子.物探化探计算技术，1985，7（1）：13－23.An Y L，Guan Z N.The regularized stable factors of removing high frequency disturbances.Computing Techniques for Geophysical and Geochemical Exploration （in Chinese），1985，7（1）：13－23.

［4］Pwalowski R S，Hansen R O.Gravity anomaly separation by Wiener filtering.Geophysics，1990，55（5）：539－548.

［5］Pwalowski R S.Green′s equivalent－layer concept in gravity band－pass filter design.Geophysics，1994，55（5）：69－76.

［6］侯遵泽，杨文采.中国重力异常的小波变换与多尺度分析.地球物理学报，1997，40（1）：85－95.Hou Z Z，Yang W C.Wavelet transform and multi－scale analysis on gravity anomalies of China.Chinese J.Geophys.（in Chinese），1997，40（1）：85－95.

［7］Fedi M，Quarta T.Wavelet analysis for the regional－residual and local separation of potential field anomalies.Geophysical Prospecting，1998，46（5）：507－525.

［8］杨文采，施志群，侯遵泽等.离散小波变换与重力异常多重分解.地球物理学报，2001，44（4）：534－541.Yang W C，Shi Z Q，Hou Z Z，et al.Discrete wavelet transform for multiple decomposition of gravity anomalies.Chinese J.Geophys.（in Chinese），2001，44（4）：534－541.

［9］Nikitin A A，Vasov O K，Belov A P，et al.Vozmozhnosti kompleksnoy geofizicheskoy interpretatsii na baze entropiynogo fil′tra.Izvestiya Akademii Nauk Turkmenskoy SSR.Seriya Fiziko－Tekhnicheskikh，Khimicheskikhi Geologicheskikh Nauk（in Russia），1984，2：79－82.

［10］郭良辉，孟小红，石磊.磁异常ΔT三维相关成像.地球物理学报，2010，53（2）：435－441.Guo L H，Meng X H，Shi L.3Dcorrelation imaging for magnetic anomaly ΔT data.Chinese J.Geophys. （in Chinese），2010，53（2）：435－441.

［11］Keating P，Pinet N.Use of non－linear filtering for the regional－residual separation of potential field data.Journal of Applied Geophysics，2011，73（4）：315－322.

［12］Pwalowski R S.Preferential continuation for potential－field anomaly enhancement.Geophysics，1995，60（2）：390－398.

［13］Zeng H，Xu D，Tan H.A model study for estimating optimum upward－continuation height for gravity separation with application to a Bouguer gravity anomaly over a mineral deposit，Jilin province，northeast China.Geophysics，2008，72（4）：145－150.

［14］许德树，曾华霖.优选延拓技术及其在中国布格重力异常图处理上的应用.现代地质，2000，14（2）：215－222.Xu D S，Zeng H L.Preferential continuation and its application to Bouguer gravity anomaly in China.Geoscience（in Chinese），2000，14（2）：215－222.

［15］Meng X H，Guo L H，Chen Z X，et al.A method for gravity anomaly separation based on preferential continuation and its application.Applied Geophysics，2009，6（3）：217－225.

［16］郭良辉，孟小红，陈召曦.优选向上延拓及其延拓高度估计.CPS／SEG Beijing 2009International Geophysical Conference＆Exposition，2009，ID1164.Guo L H，Meng X H，Chen Z X.Preferential upward continuation and the estimation of its continuation height.CPS／SEG Beijing 2009International Geophysical Conference＆Exposition，2009，ID1164.

［17］Naidu P.Spectrum of the potential field due to randomly distributed sources.Geophysics，1968，33（2）：337－345.

［18］Dampney C N G. The equivalent source technique.Geophysics，1969，34（1）：39－53.

［19］Wiener N.Extrapolation，interpolation，and smoothing of stationary time series.New York：John Wiley ＆Sons Inc，1949.

［20］Clarke G K C.Optimum second－derivative and downward continuation filters.Geophysics，1969，34（3）：424－437.

［21］Boler F M.Aeromagnetic measurements，magnetic source depths，and the Curie point isotherm in the Vale－Owyhee，Oregon［Master thesis］.Corvallis：Oregon State University，1978.

［22］Connard G，Couch R，Gemperle M.Analysis of aeromagnetic measurements from the Cascade Range in central Oregon.Geophysics，1983，48（3）：376－390.

［23］张明华，张家强.现代卫星测高重力异常分辨能力分析及在海洋资源调查中应用.物探与化探，2005，29（4）：295－298，303.Zhang M H，Zhang J Q.Resolution of modern satellite altimetric gravity anomaly and its application to marine geological survey.Geophysical ＆ Geochemical Exploration（in Chinese），2005，29（4）：295－298，303.

［24］杨金玉，张训华，张菲菲等.EGM2008地球重力模型数据在中国大陆地区的精度分析.地球物理学进展，2012，27（4）：1298－1306.Yang J Y，Zhang X H，Zhang F F，et al.On the accuracy of EGM2008earth gravitational model in Chinese Mainland.Progress in Geophysics （in Chinese），2012，27（4）：1298－1306.

［25］袁学诚.中国地球物理图集.北京：地质出版社，1996.Yuan X C.Chinese Geophysical Atlas（in Chinese）.Beijing：Geological Publishing House，1996.

［26］杨永.利用重、磁异常研究东海南部中生界分布［硕士论文］.西安：长安大学，2010.Yang Y.Research on Mesozoic strata in the southern part of east China sea by gravity and magnetic anomalies［Master thesis］（in Chinese）.Xi′an：Chang′an University，2010.

［27］熊小松.中国大陆莫霍面深度与变化特征及其地球动力学意义［博士论文］.北京：中国地质科学院，2010.Xiong X S.Moho depth and variation of the continent in China and its geodynamic implications［Ph.D.thesis］（in Chinese）.Beijing：China Academy of Geological Science，2010.

［28］Li S，Mooney W D，Fan J.Crustal structure of mainland China from deep seismic sounding data.Tectonophysics，420（1－2）：239－252.