重力梯度全张量测量作为一种有效的地球物理勘探手段，能够更好地对地下密度异常体进行特征描述和位置确定^[1]。与航空重力测量相比，航空重力梯度测量最大的吸引力在于对飞机加速度的不敏感性^{[2, 3]}。此外，重力梯度测量也比重力测量具有更高的灵敏度和分辨率。随着航空/卫星重力测量系统的成熟商用，以及在全球各地持续不断地飞行实验，今后的最大挑战在于对海量重力梯度数据的分析、处理和解释^{[4, 5]}。最近二十几年来，重力梯度全张量数据的处理和解释技术也发展十分迅速。Pedersen等对重力和磁梯度张量数据获取和处理进行了研究，提出张量不变量可以视为非线性滤波器从而增强场源信号，并就重力测量实验中，为达到实验目标所需的飞行高度和飞行测线间距给出了建议^[6]。Bell介绍重力梯度测量在油气勘探中的应用，并从功率谱的角度与重力数据进行了比较^[7]。Zhang提出张量欧拉反褶积法来处理测线，格网和非格网全张量重力梯度数据，与传统欧拉反褶积法相比利用了重力梯度的全部6个分量，得到地质结构更加明确紧凑的解，但该方法需要预先知道场源体形状^[8]。Valentin根据重力梯度数据，利用张量反褶积法实现等效源体的位置确定,并与常规反褶积法进行了对比^[9]。Mataragio等介绍了航空重力梯度测量在加拿大纽芬兰地区硫矿物的开采前景，并利用重力数据的张量不变量进行场源的轮廓增强^[10]。Beiki利用全张量重力梯度的解析信号完成反演工作，在地质结构位置边缘增强方面取得较好的效果^[11]。Beiki等利用全张量重力、磁力梯度数据的特征向量来估计地质体的深度，但需要预先确定地质体梯度变化最大的平面位置^[12]。Qruc利用重力梯度数据的倾角图和张量不变量的最大值来进行地质体的边缘探测和深度估计^{[13, 14]}。蒋甫玉等对三度体球冠模型重力梯度张量进行了研究，为寻找储油背斜构造提供理论依据^[15]。马国庆利用张量局部波数法解释位场数据，并与常规局部波数法进行了比较^[16]。Fitz Gerald等利用重力梯度多尺度边缘探测技术在巴西某地进行盆地断层网络的三维反演工作，成功解释了超过50%的地形数据信息^[17]。 Wedge等利用全张量重力梯度数据进行模型质量异常体的深度估计^{[18, 19]}。本文利用重力梯度特征向量分析法来估算模型密度异常体的深度位置，通过对重力梯度信号中添加噪声，利用多尺度分析法将重力梯度信号分解在不同频带，准确探测到了不同埋深、不同波长反映的密度异常体位置，验证了本文提出的方法在探测密度异常深度信息时，在有噪声和干扰场源环境下的有效性和准确性。最后将该方法应用于澳大利亚西澳洲Kauring地区实测重力梯度场，给出了测区内密度异常体的深度信息。

为了论述清晰，我们首先给出重力梯度的数 学推导。根据文献[23],由空间位置 x ′处的质量体v(该质量体体积仍用v表示)在空间 x 处引起的重力位为：

式中，G是牛顿万有引力常数；| x-x ′|= 是定义在局部东北天坐标系中积分点(x′ ₁,x′ ₂,x′ ₃)与观测点(x₁,x₂,x₃)之间的空间距离；ρ( x ′)是积分点处的质量源密度；dv′=dx′ ₁dx′ ₂dx′ ₃是质量源的体积单元。重力梯度张量 Γ 是重力位函数的二阶导数，有：

式中,g 是重力加速度向量。在质量源外部空间，有下面两个条件成立：

(1) ∇× g =0 表示重力梯度张量对角线分量为零，即满足拉普拉斯方程；

(2) ∇× g =0 表示重力梯度张量是对称的。

所以在质量源外部空间，重力梯度张量只有5个独立的分量。

因为重力梯度张量是实对称的，所以总可以找到正交矩阵 P 将其对角化

式中，Λ =diag(λ₁,λ₂,λ₃)； P =(P₁,P₂,P₃)分别是重力梯度张量的特征值及相对应的特征向量，上标T表示矩阵转置。式(3)说明某一观测点的重力梯度张量在以特征向量 P 的3个分量为坐标轴的坐标系中，可以消除交叉重力梯度分量，化成对角线的形式，即只有主对角线重力梯度分量。

若 λ₁是某观测点重力梯度张量对应的主特征值^[6]，即|λ₁|=max(|λ₁|,|λ₂|,|λ₃|)，则λ₁对应的特征向量 P ₁的方向则平行于观测点x和质量源质心x^mc定义的直线的方向，即有：

式中,R为观测点x与质量源质心x^mc的空间距离。

重力梯度特征向量法的基本思想是，根据空间重力梯度观测值，在观测面以下，将地形区域离散成紧密排列的单位体元阵(体元可以是棱柱体或其他多面体)，空间每个观测点重力梯度张量主特征值对应的特征向量所定义的直线必将穿过体元阵中的单位体元。如果一个单位体元被直线穿过，则该单位体元自身的值累加一次，当所有的观测点重力梯度主特征值对应的特征向量全部计算完成时，最大幅值的单位体元则表征密度异常的质心位置，该体元所在的深度就是密度异常体所在的深度^[18]。考虑到每个观测点对单位体元的贡献大小不同，我们给体元单元不同观测点累加值赋予不同的权函数。

式中, M、N分别是单位体元和重力梯度观测点的总个数；C_k^voxel是第k个体元的幅值；c_n^k是第n个观测点在第k个体元的累加单值；w_n^k为权函数 (本文定义累加单值为1，权函数为相应观测点对应的重力梯度垂直分量)。

根据位场理论，当观测点与密度异常体之间距离增加时，重力梯度高频分量比低频分量衰减的更快，损失的能量更多，这使得区分小而浅的密度异常体和大而深的密度异常体变得困难，因为它们可能在空间产生相似的重力梯度分布特征^[19]。

为了解决这个问题，本文提出了多尺度分析法，根据重力梯度的频域特征将其分解成不同的频带进行分析处理，该方法的一个优点在于可以将高频噪声和低频信号拆分在不同的频带，并且避免浅层密度异常高频信息对深层密度异常的干扰，从而准确分离解析出不同埋深密度异常体的位置信息。

文献[20]根据两个Shannon函数的乘积构造了一种连续小波低通滤波器，在过渡带内衰减速度很快，性能优越。其时频响应公式如下：

式中，ω为截止频率；a为小于截止频率的一个常数频率^[20]。

将上述小波滤波器扩展到二维空间得到：

式中，L_ηlp、L_ηhp和L_ηbp分别是新型二维小波低通、高通和带通滤波器，η是归一化频率；η_lp、η_hp分别是新型小波低通和高通滤波器归一化截止频率；η_bps、η_bpl是新型二维小波带通滤波器上下限截止频率； 是归一化常数频率。

图 1给出了不同截止频率的二维新型小波低通、带通、高通滤波器。其中低通和高通二维新型小波滤波器通带截止波长是10 m,带通滤波器1的通带上下限截止波长是25 m和15 m,带通滤波器2的通带上下限截止波长是35 m和25 m。

图 1 二维新型小波滤波器归一化频率响应 Fig. 1 Normalized Frequency Response of Modified Two-Dimensional Wavelet Filter

图选项

构建的模型是1 000×1 000 m²区域范围内不同埋深的两个棱柱体，模型统计信息见表 1，计算了10 000个点的重力梯度值(行和列方向各100个点值，数据间隔为10 m)。

解算得到的空间重力梯度见图 2，设置解算面位于地表表面，即x₃=0。从图 2中可以看出，浅层密度异常体A在地表重力梯度反映上比深层密度异常体B更加明显，说明由于埋深的不同，密度异常体B比密度异常体A损失了更多的高频信息。

图 2 模型重力梯度全张量 Fig. 2 Model Gravity Gradient Tensor

图选项

利用§1.1介绍的重力梯度特征向量法，在重力梯度观测面以下建立了平面范围为1 000×1 000 m²,深度为300 m的单位体元阵，每个体元用10×10×10 m³的立方体来表示。这样通过计算每个观测点重力梯度全张量的主特征值对应的特征向量，然后根据式(5)，将特征向量所在直线穿过的单位体元的幅值累加，得到了幅值不等的不同深度单位体元阵。将该体元阵取不同深度平面横截面，取深度范围-275~0 m、深度间距为25 m的12个平面横截面，结果见图 3。沿纵轴500 m竖切面，不同埋深模型累加幅值随深度变化关系见图 4。

图 3 不同埋深模型探测结果 Fig. 3 Detection Results of Buried Models of Different Depth

图选项

图 4 不同埋深模型累加幅值随深度变化关系 Fig. 4 Accumulative Amplitudes Variation with Depth of Different Buried Models

图选项

通过图 3和图 4可以看到，重力梯度特征向量法能够准确而直观地确定密度异常体的深度坐标，根据不同深度体元层的幅值大小，可得到不同埋深地质异常体的质心位置，如图 3(e)和图 3(k)与相邻深度体元层幅值相比，幅值更大，锁定了模型中两个密度异常体A和B的正确质心位置[200 m,500 m,-100 m]和[800 m,500 m,-250 m]。另外，对比图 3(e)和图 3(k)可发现，浅层密 度异常体A所反映的幅值要远大于深层密度异常体B的幅值，这是因为深层密度异常比浅层密度异常高频信息衰减的更多造成的。在实际重力梯度测量中，不同埋深密度异常体由于能量随深度衰减程度的不同，在重力梯度的反映上可能互相干扰，从而给数据分析和解释带来难度，此外由于背景场和测量仪器的噪声误差影响，也会给密度异常体的深度探测带来挑战。

通过向模型各重力梯度分量中加入不同大小水平的高斯噪声来模拟实际重力梯度测量中背景场或系统噪声对密度异常信号探测的影响。高斯噪声为Gn(μ,σ)，其中μ为噪声期望水平，σ为噪声标准差，设置μ=0，σ=3E。利用添加噪声后的模型重力梯度全张量进行特征向量法分析，结果见图 5。

图 5 含噪声的不同埋深模型探测结果 Fig. 5 Detection Results of Buried Models Within Noise of Different Depth

图选项

从图 5中可以看出，由于噪声水平相对于密度异常体A的信号较小，图 5(e)仍然能够分辨出密度异常体A的位置信息，然而密度异常体B由于噪声的影响，已经不能准确显示其深度位置,见图 5(k)。造成这种结果的原因是在重力梯度中加入噪声，当噪声水平与重力梯度信号在观测面幅值相当时，将使得重力梯度主特征值对应的特征向量的指向偏离密度异常体质心位置，从而引起位置探测结果的不准确。

利用第§1.2介绍的二维小波滤波器对含有噪声信号的重力梯度各分量进行不同频带滤波，然后根据重力梯度特征向量法对不同频带重力梯度信号进行分析处理，得到不同埋深对应不同长度波段信号的探测结果。图 6给出了分频滤波后埋深250 m处的密度异常平面横截面图，从图 6中可以看出，密度异常体B的重力梯度分量通过分频滤波后，位置信息能够明确显示在低频带上，这与我们前面的推断是符合的，也说明了利用重力梯度特征向量与多尺度分析相结合的方法能够有效对不同密度异常体进行深度位置探测。

图 6 分频滤波后的模型B埋深探测 Fig. 6 Depth Detection of Model B After Different Frequency Band Filtering

图选项

实测数据来源于澳大利亚西澳州地区2009年建立的航空重力/重力梯度测量Kauring试验场。该试验场距离珀斯(Perth)简达科特机场(Jandakot Airport)约115 km。见图 7。Fugro Airborne Surveys公司于2011年7月至2012年2月间，利用FALCON航空重力梯度仪(AGG)在 Kauring试验场飞行了总测线长度为1 265.3 km 的重力梯度数据，其中东北-西南方向共115条每50 m间隔的方向线(traverse line)，西北-东南方向共10条每1 000 m间隔的结线(tie line)。本文利用6个重力梯度分量，每个梯度分量164 624个测量数据进行测区地下密度异常的深度探测。注意到，测线数据是非格网分布的，这也是本文提出的算法能够处理离散非格网非同一高度水平面数据的优势之一，避免了格网化平滑过程中有用高频信号的丢失和损耗。试验场6个梯度分量(已经进行过地形改正，去除地形质量对地下密度异常探测的干扰)见图 8^{[21, 22]}。

图 7 Kauring重力梯度试验场位置 Fig. 7 Location of Kauring Gravity Gradient Test Site

图选项

图 8 Kauring试验场重力梯度全张量 Fig. 8 Gravity Gradient Tensor of Kauring Test Site

图选项

从图 8可以看出，在Kauring试验场，重力梯度测区与坐标系(东向)呈约30°分布，从图 8(f)重力梯度垂直分量Gdd中，我们发现比较明显的梯度反映位于测区中段东北和西南方向，即图中红色特征。通过利用§1提出的方法对Kauring实验场实测重力梯度数据进行处理，得到实验场内密度异常体的位置随深度变化信息，见图 9，我们分别在埋深为-10 m、-20 m、-50 m、-100 m、-250 m、-500 m、-1 000 m和-2 000 m深度截面对密度异常深度进行探测，从图 9中可以看出，试验场内密度异常主要分布在测区东北部和西南部。根据累加幅值大小，可以判定测区西南部密度异常深度较浅，主要分布在-50~-100 m深度范围；而东北部密度异常则在各层深度均有分布，在-500 m处达到幅值最大值。这与图 8(f)中重力垂直分布Gdd显示的重力梯度是一致的，西南部密度异常埋深较浅反映出较大的梯度值，东北部密度异常埋深较深，在观测面反映出中长波段信号较多。利用本文介绍的方法，对密度异常埋深可能性较大的深度范围还可以进行更为详细的划分，从而更加直观快速的判定密度异常体的深度信息。

图 9 Kauring试验场不同深度密度异常探测结果 Fig. 9 Detection Results of Different Depth Density Anomaly of Kauring Test Site

图选项

本文根据重力梯度全张量特征向量法对地下不同埋深密度异常进行了深度位置探测。通过构建模型重力梯度场验证了方法的有效性和准确性。由于不同埋深密度异常体随观测点距离不同，重力梯度信号能量的衰减程度也不同，埋深较深的密度异常体高频信号比低频信号衰减更快，所以地下埋深较浅的密度异常体会对埋深较深密度异常体产生干扰。此外，在实际重力梯度测量实验中，由于仪器和测区背景噪声等影响，也会对最后重力梯度结果的解释产生偏差。本文通过向模型重力梯度信号中添加不同水平的高斯噪声，利用多尺度分析方法对含噪声重力梯度信号进行处理，结果显示，将信号分解在不同频带后能够有效去除高频噪声或干扰场源对低频信号的影响，还可以区别出不同埋深密度异常的相互干扰，从而达到准确确定噪声背景下不同埋深密度异常深度位置的目的。最后，利用重力梯度特征向量和多尺度分析法对澳大利亚西澳洲Kauring地区实测重力梯度场进行了分析处理，给出了测区内不同埋深密度异常体的位置信息。该方法在重力梯度油气资源勘探，矿产资源调查中具有广泛的应用前景。得到的结果也可作为重力/重力梯度地质结构密度反演等的初始约束条件，从而进行更加深入的地质结构分析和解释工作。