频谱分析和滤波技术作为物探资料数据处理中的一种重要手段，是与现代计算机的发展紧密相关的。早在20世纪50年代初期，滤波技术已在多种物探资料的数据处理中起着日益重要的作用。当时二维资料处理还受到一定条件的限制，许多方法只局限于剖面解释，或局限于简单模型、简单滤波的研究，其中局部场与区域场的分离，求导和解析延拓是最早发展起来的。从20世纪50年代末期到60年代，物探数据处理开始转向重视Fou-rier变换的滤波方法和波谱的研究。计算机等值线图绘制方法和显示方法的发展使二维资料的数字处理逐渐广泛地采用。1965年快速Fourier变换(FFT)的问世，使重磁资料数据处理中的波数域方法成为主要方法。
余弦变换和Fourier变换一样都属于正交变换中的正弦类变换，其存在的条件与Fourier积分收敛条件相同，并且在某些方面具有与Fourier变换相似的性质。但余弦变换有其自身独特的优越性，对于实连续信号，能避免复数运算，而且与K-L变换(Karhunen-Loèvetransform)具有相似的性能，能够去除原信号的相关性，从而保留原信号的最大能量。自Ahmed等(1974)提出了离散余弦变换(DCT)的定义后，DCT在语音、图像编码以及数据压缩等信号处理方面得到了广泛的应用(Rao等，1990;Dinstein等，1990)。然而，目前在物探数据处理中，DCT除用于地震数据和图像压缩(Wang等，2000;Averbuch等，2001)外，在国内外还没有发现用于重力异常数据处理的相关文献。
近30年来，Hilbert变换在重磁异常正反演解释中的应用获得了较大的发展。Nabig-hian(1972)最早借助于Hilbert变换由磁场的水平分量(垂直分量)求垂直分量(水平分量);Stanley(1976，1977)提出一种以磁场水平和垂直梯度为基础的解释方法;Mo-han和Sundararajan(1982，1983)把Hilbert变换用于位场定量解释中;Sundararajan等(1996)利用改进的Hilbert变换研究了关于自然电位解释理论中场源定位问题。Hilbert变换具有可利用位场资料的全部信息以及受背景场影响较小等特点，因此可以提高物探资料数据处理的精度。利用Hilbert变换计算重力归一化总梯度是一种新的尝试。
1.异常导数的计算
图7-1给出了两种方法计算的与理论垂向、水平一阶导数对比分析图，图中可以清晰地看到，利用余弦变换计算的异常导数(图7-1中c)与理论异常导数a拟合效果非常好，除边界几个数据因重力异常的有限截断产生的吉布斯效应残留使误差较大外，数据的计算精度均很高，误差为-0.09%～5%.而利用Fourier变换计算的异常导数b尽管与理论导数曲线走向相似，但其偏离程度非常大。而余弦变换计算的异常导数和理论异常导数拟合效果非常好。
图7-1 不同方法计算的无限长水平圆柱体一阶导数对比分析图
2.密度界面反演
图7-2为采用DCT法和Parker-Oldenberg法(Parker，1973;Oldenberg，1974)反演的二维常密度单界面模型深度对比分析图。图中采用目前公认的精度较高的Parker-Olden-berg法反演的界面深度相对于理论模型深度的计算点最大误差和均方差分别为0.148km、0.013km;而DCT法反演的为0.041km、0.003km，最大误差和均方差分别降低了0.107km、0.010km。这说明DCT法的反演精度明显高于Parker-Oldenberg法，其反演精度提高了3倍多。
图7-2 二维常密度单界面模型深度反演对比分析图
3.断层断点位置反演
在直角坐标系中，取重力异常的观测面为z=0，z轴向下为正，设重力异常g(x，y，z)为此坐标系中的函数，其z、x和y方向水平一阶导数分别表示为gz(x，y，z)、gx(x，y，z)和gy(x，y，z)。则根据Thompson1982年给出的Eouler齐次方程，重力异常g(x，z)的齐次关系可写为:
(x－x0)gx(x，y，z)+(y－y0)gy(x，y，z)+(z－z0)gz(x，y，z)=－Ng(x，y，z)(7-1)式中:(x0，y0，z0)为场源的位置;(x，y，z)为重力异常计算点位置;N为构造指数。
由上式可知，若确定场源的位置，须知重力异常的水平和垂向一阶导数以及构造指数N的数值解。N可以通过模拟地质体的简单模型的理论异常和理论异常导数的关系获得，如果反演的模式确定了，则该模型的构造指数一般来说是确定的;而实测重力异常导数的计算，不同的方法所获得的数值有一定的差异，因此考虑异常一阶导数的计算精度是十分必要的，只有获得了高精度的重力异常导数，利用超定方程组的最小二乘法近似数值解才具有更高的精度。
图7-3为分别用DFT和DCT法计算的垂直台阶重力异常垂向一阶导数(图7-3a)和水平一阶导数(图7-3b)，很显然基于DFT的异常一阶导数与理论导数相比，垂向和水平导数曲线显得舒缓，与理论导数的偏差很大;用DCT法计算的重力异常垂向一阶导数的最大误差为0.460×10－9/s2，均方差为0.189×10－9/s2，水平导数最大误差为0.182×10－9/s2，均方差为0.028×10－9/s2;可见用DCT计算的异常导数具有很高的精度。因此采用DCT法必然能够获得高精度的反演结果(张凤旭，2006，2007)。
图7-3 铅垂台阶模型重力异常一阶导数
图7-4的和中圆心位置为台阶模型断面的中心点位，“+”为反演结果，从图中可以看出，无论是铅垂台阶还是倾斜台阶，均出现了有规律的噪声，这是由于导数计算精度的影响所至，但干扰点是分散的，有70%的有效点恰好集中在圆心处，如果用点密度的概念来说明反演结果，则点密度最大处便为反演结果。该结果与台阶的中心位置(圆心)恰好重合，这说明，尽管在反演中出现了噪声，但用基于DCT的Euler法仍然可以获得高精度的反演结果。
图7-4 台阶模型反演特征
综上所述，采用DCT法能够获得高精度的重力异常数据处理结果。因此本次野外重力测量资料的处理，可以称为高精度的数据处理。
4.利用Hilbert变换计算重力归一化总梯度
重力归一化总梯度法(GH法)是由前苏联学者别列兹金(В.М.Березкин)于20世纪60年代末提出的，是一种利用在较高精度下测量的重力异常来确定场源、断裂位置及密度分界面的方法，可以用于寻找贮油气藏的构造。目前，计算GH场方法主要有Fourier级数法和Fourier变换法。在前人工作的基础上，提出用Hilbert变换计算重力归一化总梯度(称Hilbert变换法)，同时在同一计算环境下，研究三种方法计算的GH场识别异常的分辨率问题。
(1)Hilbert变换的特性
给定一实连续信号f(t)，其Hilbert变换定义为
东北地球物理场与地壳演化
式中:*为卷积符号;t为时间域(空间域)变量;τ的意义同t;
东北地球物理场与地壳演化
(t)可以看成是f(t)通过一个滤波器的输出，该滤波器的单位冲击响应h(t)=1/πt。由Fourier变换的理论可知，ih(t)=i/πt的Fourier变换是符号函数sgn(ω)(ω为角频率)，因此Hilbert变换的频率响应
如果记H(ω)=H(ω)exp［iφ(ω)］，那么
东北地球物理场与地壳演化
以上分析说明，Hilbert变换是幅频特性为1的全通滤波，信号通过变换后，频率成分做90°相移，而频谱的幅度不发生变化。
(2)实连续信号Hilbert变换的通式
利用前面Hilbert变换公式和性质以及Fourier变换公式和性质，Thomas推导出实连续信号第一类Hilbert变换通式
东北地球物理场与地壳演化
Mohan等为研究场源精确定位问题，定义了改进的Hilbert变换通式(本文称第二类Hil-bert变换)
东北地球物理场与地壳演化
式中:ReF(ω)和ImF(ω)表示f(x)Fourier变换的实部和虚部。
公式(7-5)服从Hilbert变换幅频特性，也就是通过变换后，信号频率成分做90°相移，而频谱的幅度不发生变化。公式(7-6)在公式(7-5)的基础上，相位继续相移，振幅依然保持不变。
根据以上公式，可以推导出利用Hilbert变换计算重力归一化总梯度的
东北地球物理场与地壳演化
(3)模型实验及分辨率对比分析
不含油气的背斜可视为均匀密度体，其GH场中只有一个极大值，也就是只存在一个奇点;顶部含油气的背斜是非均匀密度体，它的顶部有密度亏损，由它引起GH场中有两个极大值，即背斜双侧各存在一个奇点，两者间有一个中心在油气藏内的相对极小值，即“两高夹一低”。如果探测区是已知含油气区，此低密度体就可能是油气藏反映。因此可以把GH场中“两高夹一低”的特征作为探测油气藏的解释标志。
为了探讨利用Hilbert变换研究GH场的有效性及其对异常的分辨能力，利用Fourier级数法、Fourier变换法和Hilbert变换法三种方法计算GH场值进行对比分析。其中，Fou-rier变换法和Hilbert变换法使用的圆滑滤波方法均为前文所提的组合滤波法，而且用三种方法计算GH场值进行对比分析时，谐波数N的选取均为反映油气藏的解释标志的最佳效果取值。
图7-5 计算GH场的三度体球冠断面图
采用非均匀密度的三度体球冠模型近似表示三度背斜型油气藏。计算重力异常的三度体球冠模型(截取球冠的球体半径为2.5km)的不变参数见图7-5;h1和h2分别为模型的贮油气藏厚度及底层厚度，它们为可变参数;在计算中，通过逐渐减小模型的贮油气藏厚度h1(h1与h2的和不变)，研究三种方法计算的GH场等值线变化特征，从而实现三种方法的分辨率对比分析。
计算图7-6的三度体模型不变参数均在图7-5中给定，等值线距在图名中示出，其中两侧存在的非等距等值线值见等值线上标注。
图7-6中a示出了h1=0.2km和h2=0.3km时，三种方法计算的GH场等值线图。a1为谐波数N=44时，采用Fourier级数法计算的GH场等值线图;a2为N=65时，Fourier变换法计算的GH场等值线图;a3为N=56时，Hilbert变换法计算的GH场等值线图。
当所计算的三度体储油球冠厚度(h1=0.2km，其余地质参数不变)较大时，三种方法计算的GH场在A中均明显表现出“两高夹一低”的典型标志。图7-6a1中奇点特征值分别为:‘两高’极大值5.4，‘一低’相对极小值1.6，它们的差为3.8;图7-6a2中对应的特征值分别为6.7和1.4，差为5.3;图7-6a3中分别为8.6和1.3，差为7.3。三种方法所获得的GH场奇点中，Hilbert变换法计算的奇点极大值最大，Fourier变换法的其次，Fourier级数法的最小，Hilbert变换法计算的比Fourier级数法大3.2;在双侧极大值中心的相对极小值表现的特征恰好和极大值相反，Hilbert变换法的最小;而且极大值与相对极小值的差值仍然是Hilbert变换法的最大，Fourier级数法的最小。由此可见，与其他两种方法相比，Hilbert变换法计算的GH场中，双峰异常(“两高夹一低”)表现趋势最明显，这说明，Hilbert变换法具有更高的分辨率。
在计算中还发现，当N在图中给定值上下浮动较大幅度(甚至超过15)时，三种方法计算的GH场中“两高夹一低”的标志都很明显。
图7-6中b给出了h1=0.1km，h2=0.4km与谐波数为N=49(Fourier级数法)、N=67(Fourier变换法)，N=57(Hilbert变换法)时，三种方法计算的GH场等值线图。图中奇点极大值、相对极小值以及它们的差值分别为:图b1为5.8、1.3、4.5，图b2为6.8、1.2、5.6，图b3为8.7、1.0、7.7，其变化规律和图7-6a相似。对比分析图7-6a和图7-6b，图7-6b中“两高夹一低”标志没有图7-6a中的表现明显，这和低密度体厚度减小有关。但通过分析图7-6a、图7-6b的GH场奇点极大值、相对极小值和它们的差值，并充分考虑三种方法获得的GH场等值线特征，更明显地看出，利用Fourier变换计算GH的方法分辨率优于别列兹金法，Hilbert变换法的分辨率最高，它明显优于前两种方法。但和a相比能够得到明显的“两高夹一低”标志的谐波数N的取值变化范围缩小，三种方法的变化范围分别为:Fourier级数法44～52，Fourier变换法58～74，Hilbert变换法44～69。
图7-6 三种方法计算的GH场等值线特征图
图7-6中c给出了h1=0.05km、h2=0.45km与谐波数N=49、N=67、N=57时，三种方法计算的GH场等值线对比分析图。图中奇点极大值的变化规律和前面相关叙述相似。但和图7-6a、7-6b相比，图7-6c中c1和c2只存在一个极大值和一个相对极小值，这说明，当油气藏厚度减小到一定值时，Fourier级数法和Fourier变换法已经不能识别“两高夹一低”标志。大量模型实验证实，当三度体储油球冠油气藏部分(低密度体)厚度减小到低于球冠总厚度的十分之一时，若采用Fourier级数法和Fourier变换法计算，则无论N取何值都得不到能充分识别“两高夹一底”标志的GH场特征图(图7-6c1、7-6c2为最佳效果图)，而图7-6c3在N取40～59之间仍然可识别“两高夹一低”标志。这进一步证实用Hilbert变换法计算的GH场，其分辨油气藏的能力明显高于其他两种方法。
内容全部来源于网络收集，如有侵权，请联系网站删除：QQ：24596024