航空重力梯度测量系统及数据预处理方法

郑 崴，陈学锋

（1.中国地质大学 （北京）地球物理与信息技术学院，北京 100083；2.河南高等机电专科学校，河南 新乡 453000；3.河南科技学院信息工程学院，河南 新乡 453003）

1 绪言

重力勘探是地球物理勘探的基本方法之一。它是在地球表面上利用重力仪，测量由于地下物质密度分布不均匀引起的重力变化（或重力异常），这种物质密度的分布与地质构造及矿产分布有密切的关系。因此，通过研究重力异常，可以了解和推断地球的结构、地壳的构造，以及勘探矿产资源等等。最早用于重力勘探的设备是测量重力梯度的扭称，但是由于设备笨重，效率低下，在20世纪50年代，被更轻便，效率较高的重力仪所取代［1］。可以说重力勘探是由重力梯度测量开始的。20世纪70年代以后，出于军事方面应用的需求，重力梯度仪获得了飞速的发展，美国、法国等国家先后研制出了不同原理、不同结构的重力梯度仪。而现代电子技术、计算机技术、超导量子干涉技术、低温微波空腔谐振技术、超导负弹簧技术的发展和应用，重力梯度仪的灵敏度和稳定性取得了突破性的进展。由于重力梯度相对于重力异常具有更高的分辨率，随着重力梯度仪的发展，它在地质资源勘查、地质调查等领域发挥出重要的作用。

相对于地面测量技术来讲，航空重力梯度测量具有快速、经济和高精度的特点，它几乎可以到达任何测量区域。除了应用于精细地质勘探外，还应用于国防、军事领域。因此，航空重力梯度测量技术是重力勘探中最新、最先进的技术。本文首先介绍航空重力梯度测量的基本概念与发展现状，然后讨论标准傅里叶变换在航空重力梯度数据预处理中的应用。

2 重力和重力梯度

在地球表面任一点的重力是由地球引力与地球自转产生的惯性离心力的合力，它与在地球表面的位置和物体的质量有关。三度体在u＝（x，y，z）处的重力场可以用重力位的标量函数

式中，G为万有引力常数，为u'＝（x'，y'，z'）物质的密度。

在笛卡尔坐标系中，可以从重力位推导出重力和重力梯度的表达式。重力是重力位的一阶导数，即

重力仪测量的是垂直方向分量gz，即

而重力梯度是重力位的二阶导数，即重力加速度矢量的三个分量分别沿直角坐标系三个方向的一阶导数。重力梯度是一个张量，可以表示为式（4）中的形式。图1位重力场矢量及其梯度示意图。

在无源的空间内，引力的散度和旋度均为0，所以重力梯度是一个对称的张量，有TXY＝TYX，TXZ＝TZX，TYZ＝TZY，且有TXX＋TYY＋TZZ＝0，因此，重力梯度张量的9个分量中，只有5个是独立的分量。

在球坐标系中，设有一质量为M的均质球体，则有

图1 重力场矢量及其梯度示意图

式中，r为测量点到目标体中心的距离，由式（5）可见，物体的密度与距离r成反比，重力与r2成反比而重力梯度与距离的3次方r3成反比，可见，重力梯度比重力衰减的快，利用重力梯度可以更好的估计物体边缘的位置，所以说重力梯度相比于重力来讲具有更高的分辨率。而重力梯度异常值反映的是由于地下密度分布不均匀的重力异常在空间的变化率，重力梯度张量的各分量，可以反映出目标地质体不同的地质特征，如图2（a）是针对已知地质目标完成的一次航空重力梯度测量所获取的结果，已知目标地质体为盐类的盖层，表现为TZZ图中圆圈部分所显示的正值异常。而在图2（b）中的各个张量分量则可以明显的分辨出目标体的地质特征，其中TXX和TXY可以显示出目标体的南北向的边缘位置，从图中能够明显的分辨出目标体的南北向的位置轮廓；同样的，由TYY和TYZ图可以分辨出目标体的东西向的边缘位置。而由TXY图则显示出盖层具有典型的2高点和2低点的特点4极异常。

图2 已知盐类盖层航空重力梯度测量张量分量特征

3 航空重力测量的基本原理和航空重力梯度测量系统

根据牛顿第二定律，作用于单位质点上的总加速度矢量fi（称为比力），与载体运动加速度矢量ri和引力加速度矢量Gi的关系可以表示为

因此，要获得引力加速度Gi，必须获得载体运动加速度ri，再从测得的比力fi中将其减去。在地面静止状态下，重力仪的运动加速度ri＝0，此时重力仪测得的就是重力场矢量。而在航空重力梯度测量中，由于运动加速度不断的变化且不为零，就必须在测得的比力fi中分离出运动加速度ri，而分离方法是通过对共用基线的两个加速度计的输出求差，以消除载体运动的影响。如果共用基线是旋转稳定的，由差值读数可以获得重力梯度分量，这也是重力梯度测量的基本原理。除了要消除载体运动加速度带来的影响，还要保证在测量过程当中，加速度计保持在平稳、稳定的工作状态，这一点主要是通过硬件来实现，如可以将加速度计安装在INS（Inertial Navigation System）惯导平台等。

这里主要介绍两种已经投入商用的航空重力梯度测量系统。

3.1 基于旋转加速度计的重力梯度系统

旋转加速度计重力梯度系统是市面上唯一投入商用的航空重力梯度系统，又分为全张量测量系统和部分张量测量系统。

全张量重力梯度测量系统是由美国Bell Aerospace（现已并入Lockheed Martin公司）自主研究制造的Air-FTG系统。它的核心器件都是图3所示的旋转加速度仪（GGI，Gravity Gradient Instrument），它是将两对相互正交的加速度计安装在一个圆底盘的边缘上。GGI旋转工作时，每对加速度计感应的圆盘的运动加速度大小相等方向相反，可以互相抵消。通过计算每对加速度计间的距离和底盘的旋转频率，可以测得一个平面上的两个梯度分量。一个全张量梯度测量系统包括三个GGI，这三个GGI的轴线互相垂直，每个GGI的轴线沿铅垂线以同样的角度相交，从上向下看，三条轴线的投影角度为120度，如图4所示。为保证测量精度，GGI必须以恒定的速度绕垂直轴方向旋转。同时，为了隔离载体运动过程对系统的影响，可以将旋转加速度仪安装在由陀螺仪稳定的惯性平台上。另外，系统配备独立的坐标系，通过将三个GGI在外部坐标上的梯度分量进行线性组合，可以计算出重力梯度张量的5个独立分量，进而得到重力梯度张量的全部分量值。

图3 加速度计安装及转盘运动示意图

图4 Air-FTG的三轴GGI示意图

部分重力梯度张量测量系统的代表是由澳大利亚的BHB Billiton公司出资研制的FALCON系统，此系统于1997年研制成功，1999年投入商业航空物理勘探。与全张量系统相比，FALCON只配备一个 GGI，它可以直接测量出TXY，（TXX－TYY）/2两个梯度分量，再通过面积性测量和数学变换得到重力水平分量gx和梯度分量TXX，根据实际需要，它还可以由TXY计算出其他的张量分量，也就是说FALCON系统可以通过数学计算的方法得到全部张量分量。

3.2 超导航空重力梯度系统

超导航空重力梯度系统（SGG）的研制源于20世纪70～80年代，近十年来，超导航空重力系统作为下一代的主力航空重力勘探系统，已经成为目前研究的重点和热点。超导航空重力梯度系统的设计指标是具有更高的灵敏度和更精确的测量结果，设计噪声指标优于1E。

英国的ARKeX公司正在改造用于卫星上超导重力梯度系统，由于卫星的飞行比较平稳，而飞机的加速度很大，ARKeX改造的重点就在克服飞机的飞行干扰上。ARKeX公司的最终目标是生产经营世界上首台超导勘探重力梯度仪EGG（Exploration Gravity Gradiometer）。EGG的设计可配备六个加速度计模块，用来测量全部的梯度张量分量，但在实验过程中制造出的第一台EGG只用了两个加速度计模块，只能测量垂直方向的张量元素TZZ。

4 航空重力梯度测量的数据预处理

全张量重力梯度测量系统是目前应用较广泛的航空重力梯度测量系统，以它的数据处理过程来介绍航空重力梯度的数据预处理方法。FTG的数据处理过程主要有两个阶段：数据预处理阶段和后任务处理阶段。第一个阶段是实现对加速度计的输出值因外界测量环境的变化而进行的补偿。加速度计是在飞机高速运动的状态下完成测量，因此，在这个阶段，需要对量测量进行一些必要的校正，这些校正主要针对以下几种影响：①每一个AGG的底部磁盘高速旋转产生的向心力的影响；②飞机飞行方向变化和飞行姿态变化的影响；③飞机飞行过程中燃油消耗量变化的影响④飞行过程中加速度不断变化的影响。

当这些校正完成后，还要对校正后的数据进行插值，以供第二个数据处理阶段——后任务处理阶段使用。在预处理的最后，需要将校正后的数据重新插值到0.5Hz，同时解算出每一个张量元素值，并将其换算到外部现实世界的直角坐标系统。而后任务处理阶段的主要任务是降噪，即在存在大量随机干扰的测量数据中获得有用数据的信息，使得输出的张量分量具有较为满意的信噪比。这需要多种滤波技术的支持来实现。

4.1 傅里叶变换解算重力梯度张量分量

不难看出，在数据预处理阶段，解算出重力梯度张量的各个分量并实现坐标轴的变换，在数据预处理阶段是非常重要的，为了获得符合解释要求的数据，还要保证解算出的张量具有可观的信噪比（S/N）。这种数据预处理的要求，可以由标准傅里叶变换的方法来实现。因为重力梯度张量实际是重力位的在各个坐标轴方向的二阶导数，利用重力梯度的傅里叶变换就可以导出相应的重力位。重力位是重力测量中很重要的物理量，一旦确定了重力位，所有的重力分量都可以由它得到。并且，在进行傅里叶变换的同时，可以设计合适的加权函数来抑制噪声的影响，达到提高信噪比的目的。

设在测量过程中飞机的飞行高度是不变的，在高度为Z的X-Y水平面内，重力位的傅里叶变换对为

对式（9）的傅里叶变换部分求积分，得

不难看出，重力的傅里叶变换和重力位的傅里叶变换存在着一定的关系，再经过反变换，可得其重力位

所以，在高度为Z水平面内，只要测出gx，gy，TXX，TYY，TXY五个量中的任意一个分量，就可以通过相应的变换关系得出重力位，此时，所有X，Y方向的重力的导数就都可以计算出来［11］。而在实际的数据处理过程中，设需要处理的FTG测量数据已经实现网格化，按照网格的水平方向和垂直方向对数据进行傅里叶变换，而垂直方向的傅里叶变换相对于水平方向的傅里叶变换更容易实现，此时式（7）的关系可表示为式（12）。

在进行傅里叶变换时，可以将整个待处理数据划分为多个小的数据块，每个小的数据块又有多个网格组成。对这些由多个网格组成的数据进行积分运算，可推出相应的重力位，再求重力位的均值，可以提高输出的信噪比，再进行微分变换，则可以推导出其他的张量分量，如TZZ。在这个计算过程中，为导出具有较高信噪比的重力分量的均值，可以通过将重力分量转换到频率域内，根据量测数据的频谱结构，设计相应的加权函数滤除域内的高频噪声或者放大有用信号的功率，以此来达到提高信噪比的目的。

所以，在理论上来讲，由重力位或其任意导数的傅里叶变换都可以导出梯度张量的各个分量，反之亦然。但在实际应用中，有些变换是无法直接实现的，比如说利用TXX求垂直方向的重力分量gz，理论上可以由式（14）直接计算得到。

但是，在式（14）中，由于TXX是沿X轴方向分量，它不包含垂直方向Z的信息，这就给数据的处理带来两个问题。第一个问题，式（14）中的分母kx≈0，即无法利用TXX的傅里叶变换求得重力分量gz。第二个问题，接近ky轴的地方，TXX的幅值衰减的很快，而噪声信号的幅值却很大，若在k域内直接进行傅里叶变换，相当于对噪声信号进行放大，这必然导致信噪比的降低。

对于噪声的影响，可以有多种方法来解决。可以将滤波器的截止频率设置在远离ky轴的位置上，或者设计维纳滤波器来减弱噪声的影响。而对于第一个问题，常用的方法是，利用已测得的包含垂直方向的张量分量（如TXZ，TYZ）来导出gz，即进行傅里叶变换前对这些测得的张量分量进行线性组合。经过线性组合后，有用信号的频谱会较为完整的保留下来，再去利用傅里叶变换去导出重力位，则可以计算出具有较高信噪比重力分量gz。但是这种方法的有两方面的问题，首先，进行数据的线性组合需要完成大量的运算。其次，线性组合后会给降噪带来一定的难度，如果有用的数据与噪声数据组合在一起，数据变换后很难用滤波器去削弱这部分噪声的影响。所以，利用线性组合来解决这个问题并不是最好的选择。

而FTG系统则是用硬件设计的方式来解决这一问题，一个FTG系统包含三个GGI，每个GGI的轴线与垂直方向有相同的夹角120°。也就是说，每个GGI的输出的张量分量当中都包含有关垂直方向的梯度元素TZZ的分量，因此，GGI输出结果中包含垂直方向的信息。设三个GGI输出的横差曲率分别表示为C1，C2和C3，则梯度分量TZZ可以直接有公式（15）计算得出

4.2 坐标系的变换

实测数据是基于测量坐标系（X轴指向地理东，Y轴指向地理北，Z轴指向垂直方向），由于在傅里叶变换中要用到旋转矩阵，这增加了数学计算的难度。为此，可以利用式（13）构建一个差分矩阵因子，实现由重力位到直角坐标系的梯度张量TXYZ的转换

如果在球坐标系中，有一重力梯度张量TUVW，则定义一个矩阵R，将其由直角坐标系换算到球坐标系中，此时的变换公式表示为式（17）

举个例子来讲，设测量所得的梯度分量为Tuw，已知它是一个包含有横差曲率的分量（类似前边的C3），则它的傅里叶变换表示为式（18）

只有在坐标原点上，式（17）的分母为0。而丢失的数据部分为一平均值，在整个数据处理过程它是一个常量，因此不会影响后续的处理结果。

利用式（18）对TUVW的各个张量元素进行推导，可以导出相应的重力位，同时，傅里叶变换后，平均运算将有用信号平均分配到每个张量分量中，以此来实现达到提高有用信号功率，抑制噪声的目的。

5 总结

地球重力梯度的研究一直是地球物理勘探领域最前沿、最热点的研究课题，重力梯度实际上是重力异常在空间内各个方向上的变化率，可以更好的反映出探测目标的地质信息，相对于重力异常而言重力梯度异常具有更高的分辨率。同时，依托于航空平台的重力梯度测量系统相对于传统的地面测量仪器来说，又具备快速，高效，经济的特点，具有广阔的研究前景和价值。

由于其是在高速动态的状态下完成实地的测量，航空重力梯度测量的数据处理步骤不同于传统重力测量的处理方法，在数据预处理的过程中，利用傅里叶变换可以快速的对重力梯度数据进行预处理，以获得相应的重力梯度张量分量和所需要的重力分量，同时也可以快速实现梯度张量在球坐标与直角坐标间的变换。

［1］曾华霖.重力场与重力勘探［M］.北京：地质出版社，2005.

［2］张金焕，李晓斌.航空重力梯度测量系统发展及现状［J］.中国矿业，2011，20（zk）：188-201.

［3］张永明，张贵宾，盛君.航空重力梯度测量技术及应用［J］.工程地球物理学报，2006，3（5）：375-380.

［4］张永明，张贵宾，盛君.航空重力梯度测量的基本理论与应用［J］.地质装备，2006，7（6）：23-27.

［5］张昌达.航空重力测量和航空重力梯度测量问题［J］.工程地球物理学报，2005，2（4）282-290.

［6］纪兵，刘敏，吕良，等.重力梯度仪在水下航行中的应用［J］.海洋测绘，2010，30（7）：23-25.

［7］孙中苗，翟振和.航空重力测量数据的小波阈值滤波［J］.武汉大学学报：信息科学版，2009，34（10）：1222-1225.

［8］Richard Lane.Airborne Gravity 2004：Abstract from the ACEG-PESA Airborne Gravity 2004workshop ［A］.Geoscience Australia，2004.

［9］Scott Hammond，Colm Murphy.Air-FTGTM：Bell Geospace's Airborne Gravity Gradiometer-A Description and Case Study.Preview，August 2003：24-26.

［10］王静波，熊盛青，周锡华，等.航空重力测量系统研究进展［J］.物探与化探，2009，33（4）：368-373.

［11］张永明，张贵宾，盛君.航空重力梯度测量的基本理论及应用［J］.地质装备，2006，7（6）：23-27.

［12］舒晴，周坚鑫，尹航，等.应用和研制中的航空重力梯度测量系统［C］.中国地球物理年会论文集，2011.

［13］Zhangxin Liu，Hai Lu.Research on Correlation Imaging of Gravity Gradient Based on FFT［C］.2011Xian International Conference on Fine Geological Exploration and Groundwater＆ Gas Hazards Control in Coal Mines，2011：203-209.

［14］Majid Beiki.Analytic Signals of Gravity Gradient Tensor and Their Application to Estimate Source Location ［J］.Geophysics，2010，75（6）：159-173.

［15］Gary Barnes，John Lumley.Processing Gravity Gradient Data［J］.Geophysics，2011，76（2）：133-147.

［16］Glennie，C.，and K.P.Schwarz.Airborne gravity by strapdown INS＝DPGS in a 100km by 100km area of the Rocky Mountains：Proceedings of the International Symposium on Kinematic Systems in Geodesy［J］.Geomatics and Navigation（KIS97），1997：619-624.

［17］Vassiliou，A.A.Stage II processing of airborne gravity gradiometer data using frequency domain techniques［C］.Proceedings of the 14th Annual Gravity Gradiometry Conference，1986：435-440.

［18］Dransfield，M.H.，Conforming Falcon gravity and the global gravity anomaly.Geophysical Prospecting，2010，58，469-483.

［19］Li Y.，and D.W.Oldenburg，Rapid construction of equivalent sources using wavelets［J］.Geophysics，2010，75（3），51-59.

［20］Saad，A.H.，Understanding gravity gradients—A tutorial［J］.The Leading Edge，2006（25），942-949.