基于界面反演增强的位场边缘识别方法

冯旭亮，魏泽坤

(1.西安石油大学 陕西省油气成藏地质学重点实验室，陕西 西安 710065； 2.西安石油大学 地球科学与工程学院，陕西 西安 710065)

0 引言

油气资源、深部矿体等的形成与分布往往与断裂、火成岩体等有着直接或间接的关系，因此利用地球物理方法探测断裂或岩体的位置在油气藏勘查和深部找矿中起着至关重要的作用[1-4]，这些断裂或岩体的位置可统称为地质边界。重、磁方法在研究地质体横向不均匀性、特别是地质边界方面有独特的优势。利用重、磁位场资料识别地质体边界的方法可分为数理统计、数值计算和其他方法[5]，其中数值计算类方法研究和使用最多，比较有代表性方法主要有垂向导数[6-7]、总水平导数[8-9]、解析信号振幅[10-12]、倾斜角(Tilt导数)[13-15]、θ图[16]5种基本方法，以及在这5种基本方法上发展起来的其他方法，如倾斜角总水平导数[17]、改进Tilt梯度水平总梯度[18]等。

王万银等[5]概括了各类位场边缘识别方法的研究现状，重点对数值计算类中垂向导数、总水平导数、解析信号振幅、倾斜角、θ图这5种基本方法及在其基础上发展的各方法原理和应用效果进行了研究。Ekinci等[19]对以上5种常用的边缘识别方法进行了系统研究，并处理了土耳其安纳托利亚半岛地区的实际资料，与已知构造的对比结果证实了以上边缘识别方法的有效性。为进一步增强位场边缘识别方法的处理效果，特别是对深部小构造的识别能力，一种措施为对边缘识别方法处理结果进行滤波增强[20-23]。不同基本方法的组合也是常用的提高边缘识别效果的措施，如倾斜角解析信号模法[24]、解析信号模倾斜角法[25]等。英高海等[26]对比了15种具有代表性的边界增强方法，认为解析信号模垂向导数、总水平导数、解析信号模倾斜角法是对于磁异常适用性较强、应用效果较好的3种方法。

为提高弱异常(深部异常体)的边缘识别效果，一些学者研究了均衡类方法。马国庆等[27]利用水平与垂直导数的相关系数进行位场边缘识别，并基于不同阶水平导数构造了多个增强型的位场边缘识别方法[28-29]。此外，王彦国等[30]利用不同阶次的垂向导数的比值进行位场边缘检测，Du等[31]则利用垂向导数的平均值与标准差的相关系数进行位场边缘检测，于平等[32]将水平方向解析信号及其垂向导数与传统的均衡边界识别方法结合形成了位场边界探测方法。郭灿文等[33]利用水平导数与垂向导数的标准差的相关系数识别了磁源边界。Pham等[34]利用总水平导数构建逻辑函数并据此识别位场边界，可明显增强边界识别的效果。

归一化类方法也是增强深部场源弱异常识别能力的主要措施，如归一化总水平导数法[35-37]、归一化θ图法[38]。Wang等[39]提出了归一化总水平导数垂向导数(NVDR_THDR)边缘识别技术，该技术不但具有边缘探测功能，并且具有边缘增强功能，在划分断裂[40-42]、解释岩体边界[43]等方面取得了较好的应用效果。我们曾在识别裂陷盆地基底不同区域界面形态时，尝试将直接迭代密度界面反演方法[44-45]加入到NVDR_THDR位场边缘识别技术之中，提高了基底深部小规模断裂的识别能力[46]，但并未对该方法的细节及应用效果进行详细研究。本文将详细阐述这一技术的原理，并通过不同类型的模型试算及实际资料处理，以分析该方法的正确性及应用效果。

1 方法原理

基于界面反演增强的位场边缘识别方法由界面反演技术和位场边缘识别技术组成。首先利用密度界面反演方法对重力异常进行处理(若为磁异常，可利用磁性界面反演方法进行，或转换为磁源重力异常)，之后利用位场边缘识别技术对密度界面反演结果进行处理。其中密度界面反演技术采用直接迭代法计算，计算公式为[45,47]：

i=1,2,…,N

(1)

式中pk(xi,yi)和pk-1(xi,yi)分别为第i个计算点处第k次和第k-1次的反演结果；G为万有引力常量；Δρ(xi,yi)为第i个计算点处的密度差；go(xi,yi)和gcal(xi,yi)分别为第i个点处实测的重力异常和第k-1次迭代时该点的正演拟合重力异常。计算拟合重力异常时将密度界面与参考界面之间的物质剖分为垂直并置的棱柱体，重力计算点位于棱柱体正上方，其个数与计算点数相同，水平尺寸与计算点间距一致，计算公式为：

(2)

式中非线性函数fi(pj)用来计算第j个棱柱体在第i点处引起的重力异常，其表达式为：

(3)

式中：(dξj,dηj,dζj)为剖分棱柱体内任一微元的坐标，ξj1～ξj2、ηj1～ηj2、ζj1～ζj2分别为第j个棱柱体在3个方向的坐标范围。

利用式(1)反演密度界面深度时，通常令p0为零开始迭代，并利用go(xi,yi)与gcal(xi,yi)的均方差小于重力数据的噪声水平而终止迭代计算。然而，本文重点将其应用于深部场源体异常的增强，并不要求反演结果有较高的准确性，因此一般迭代3～5次即可。需要说明的是，对于多个孤立地质体引起的异常，从几何形态而言，并非密度界面，但其引起的异常也可利用该方法进行增强。

本文选用归一化总水平导数垂向导数(NVDR_THDR)技术作为位场边缘识别方法，用该技术直接处理密度界面的反演结果。结合本文处理流程，其方法原理[39]如下：

计算密度界面反演结果总水平导数(THDR)：

(4)

式中：f(x,y,z)为密度界面反演结果。

计算总水平导数的n阶垂向导数：

(5)

式中：n为垂向导数的阶次，阶次越高，计算结果的水平分辨率越高，但计算结果不稳定，阶次一般选1或2即可。

计算总水平导数的峰值：

PTHDR(x,y,z)=

(6)

计算总水平导数峰值与总水平导数的比值：

VDR_THDR(x,y,z)=

(7)

计算总水平导数垂向导数的极大值并利用其对总水平导数垂向导数归一化：

(8)

2 模型试算

2.1 裂陷盆地基底模型

裂陷盆地模型如图1a所示，该盆地由7个大小不一的次盆组成(编号为A～G)，基底最大深度为6 km。假设沉积层与基底的密度差为-0.4×103kg/m3，将盆地的沉积层剖分为水平尺寸为1 km×1 km垂直并置的棱柱体计算其重力异常，并加入均值为0、标准差为0.1×10-5m/s2的高斯白噪声，作为盆地模型的实测重力异常，如图1b所示。理论重力异常基本能反映各次盆的位置及形态，但距离较近的几个次盆的边界(实质为断裂)无法区分。

利用本文方法对图1b所示的重力异常进行处理，其中在利用密度界面反演技术时，给定密度差为-0.4×103kg/m3，迭代5次的结果如图2a所示。与原始重力异常(图1b)相比，图2a的反演结果中可以明显区分7个次盆的位置及形态，边界较为清晰，利用NVDR_THDR对密度界面反演结果做进一步处理，结果如图2b所示，其极大值位置与理论模型中次盆的边界较为吻合，并且也能在一定程度上区分距离较近的几个次盆的范围，如次盆B、C、D以及次盆E、F。利用NVDR_THDR方法直接对重力异常进行处理，结果如图2c所示。与图2b相比，图2c的结果中次盆B的右边界不完整，无法呈现E和F的南北边界及分界线。为提高小规模构造的识别能力，在NVDR_THDR技术中垂向导数的阶次选为2重新处理，结果如图2d所示，与垂向一阶导数结果(图2c)相比，垂向导数阶次为2时计算结果幅值更窄，但其仍然无法有效识别上述次盆的边界。

本文方法的实质是将重力异常转换为一种与真实构造更为接近的“异常”从而提高边缘识别能力，因此若能找到一种与理论构造更为接近的处理结果或对重力异常局部增强结果并用NVDR_THDR进行边缘识别，是否也能得到理想的效果？图2e～h分别为采用位场分离、稳定向下延拓、倾斜角和解析信号振幅对重力异常处理并利用NVDR_THDR对异常进行提取的结果。其中，剩余重力异常的NVDR_THDR结果(图2e)极大值幅值较窄，也能清晰地呈现次盆B的东部边界，但无法识别次盆D、E和F的南北边界。采用稳定向下延拓技术向下延拓3 km重力异常的NVDR_THDR结果(图2f)与剩余异常结果较为相似，但其仍然无法反映E和F的南北边界。倾斜角的NVDR_THDR结果(图2g)与垂向导数阶次为2的NVDR_THDR结果(图2d)较为类似，但其幅值更宽一些。解析信号振幅的NVDR_THDR结果(图2h)对于两个规模较大的次盆A和G识别较好，但对于规模较小且距离较近的次盆无法有效识别其边界。图2的结果表明，本文提出的方法更有利于盆地基底断裂的识别，尤其对于小规模断裂构造更为有效。

图1 裂陷盆地模型及其重力异常Fig.1 A rifted basin model and its gravity anomaly

利用密度界面反演方法增强位场边缘识别效果时，界面上下的密度差是一个关键参数。图2a为密度差是真实值的结果，若密度差取值不准确，是否会对位场边缘识别结果有影响？为此，分别给定密度差为-0.3×103kg/m3和-0.5×103kg/m3进行反演计算，结果分别如图3a和3c所示，相应地NVDR_THDR处理结果如图3b和3d所示。当密度差不准确时，直接迭代法均能得到密度界面的起伏变化，尽管其与理论模型的真实深度差别较大，但相对起伏变化是一致的。在此情形下，密度界面反演结果的NVDR_THDR计算结果与准确密度差时的处理结果几乎是一致的，均能较清晰地反映理论模型中次盆的边界位置。可见，密度差的选择几乎不会影响最终的位场边缘识别结果。

图2 裂陷盆地模型位场边缘识别结果Fig.2 The edge recognition results of the rifted basin model

图3 密度差不同时裂陷盆地模型位场边缘识别结果Fig.3 The edge recognition results of the rifted basin model with different density contrasts

2.2 孤立形体组合模型

孤立形体组合模型如图4所示，该组合模型由7个大小不一、埋深不同的立方体组成，组合模型的观测面为z=0 km的平面，x和y方向的范围均为0～20 km，观测面网格间距为0.25 km×0.25 km，共81×81个观测点。正演计算该组合模型引起的重力异常，并加入均值为0、标准差为0.05×10-5m/s2的高斯白噪声，最终得到的理论重力异常如图4c所示。

利用本文方法进行位场边缘识别时，需要将理论重力异常视为密度界面引起的重力异常并反演该“界面”的起伏变化，反演结果取决于界面上下密度差及界面平均深度。为保证反演的密度界面起伏变化与重力异常大小变化呈正相关以利于直观分析，密度差应取为负值。尽管盆地模型试验结果表明了密度差几乎不影响最终的位场边缘识别结果，但应尽可能地保证界面反演结果较浅以使其横向分辨力较高。实际资料处理时可根据岩石物性特征选择一个符合实际情况的密度差即可，这里选为-0.25×103kg/m3，界面平均深度取为2 km，反演结果如图5a所示，相应的NVDR_THDR处理结果见图5b。为便于对比，本文亦分别计算了模型理论重力异常的NVDR_THDR垂向导数阶次为1和2的结果，见图5c和5d。此外，裂陷盆地基底模型试验中，剩余重力异常和向下延拓重力异常的NVDR_THDR结果较佳，因此，本文也对孤立模型进行相同的处理，结果分别见图5e和5f。

与理论重力异常(图4c)相比，密度界面反演结果(图5a)的形态更接近理论模型的位置和形态，利用NVDR_THDR技术处理的结果(图5b)中，除模型E的左边界之外，浅部的5个模型C～G的边界较为清晰，深部模型B的边界位置也基本能识别出来。直接利用NVDR_THDR技术对理论重力异常的处理结果(图5c)中，深部模型B的位置有所显示，但其右侧边界偏差较大；浅部模型中，仅能完整的识别F和G的边界。垂向导数阶次为2时NVDR_THDR处理结果(图5d)中，异常的幅值较窄，更容易确定模型边界，且深部模型B的位置也更准确一些，但仍然无法完整识别浅部模型C、D、E的边界。剩余重力异常NVDR_THDR结果(图5e)与本文方法计算结果(图5b)较为相似，但异常极大值稍宽一些，模型C的边界也不太清晰。向下延拓重力异常NVDR_THDR结果(图5f)可以清晰地反映浅部模型C～G的边界， 甚至模型E的左边界也能显示出来，但对于深部模型B则无法有效反映，NVDR_THDR的幅值与噪声相当。不难理解，重力异常向下延拓之后，可以突出浅部异常、压制深部异常，因此，处理结果中只能突出浅部异常的边界。可见，对于孤立形体的组合模型，本文方法的优势不如在盆地基底断裂识别中明显，但其在一定程度上可平衡深部和浅部异常的识别。

图4 孤立形体组合及其重力异常Fig.4 The combination model and its gravity anomaly

3 实际资料处理

为测试本文方法在实际资料处理中的效果，选择鄂尔多斯盆地北部重力数据(图6a)进行计算，该数据来源于区域1∶50万实测布格重力异常。图6b为利用NVDR_THDR技术直接对布格重力异常的处理结果，研究区北部NVDR_THDR极大值以近EW向为主，其次是NE向和NEE向，极大值分布较为密集。研究区中部和南部NVDR_THDR极大值以NW向、EW向和NE向为主，由于该区域重力异常变化相对较为平缓，因此NVDR_THDR的极大值分布较为稀疏。研究区西部NVDR_THDR极大值密集分布，走向以近SN向为主。

利用本文方法对研究区重力异常进行处理，反演时为得到深度较合理的密度界面，首先对布格重力异常添加常数的背景重力值，将其调整为(-120～0)×10-5m/s2之间。之后，根据整个鄂尔多斯盆地的地层密度资料，并为保证反演界面尽量较浅以利于突出局部异常，给定密度差为-0.5×103kg/m3，得到界面如图6c所示。利用NVDR_THDR技术对该界面处理结果如图6d所示，图中NVDR_THDR的极大值分布更密集一些，尤其在研究区中部和南部更为明显。研究区北部和西部的极大值分布与图6b中的较为相似，但异常幅值更窄，局部异常得到增强。

为进一步对比位场边缘识别结果，沿研究区的一条地震剖面[48](图6d中黑色实线)提取了图6b和6d的处理结果，如图7所示。位场边缘识别结果仅能反映断裂或岩性边界的平面位置，因此NVDR_THDR的极大值对应于断裂带中点偏上的某一个位置。图7中，断距较大的断裂附近重力异常NVDR_THDR和密度界面NVDR_THDR均为极大值，但重力异常NVDR_THDR无法有效识别断距较小的断裂。密度界面NVDR_THDR结果中，对于断距较小的断裂也能识别，例如剖面长度110、140、160 km处，这一特征尤为明显。可见，本文提出的基于界面反演增强的位场边缘识别方法在实际资料处理也较为有效。需要说明的是，理论上，当断裂埋深较浅、断距较大时，其引起的重力异常梯度也越大，相应的NVDR_THDR的值也越大。但实际地质情况较为复杂，NVDR_THDR的大小与断裂的规模、埋深以及断裂两侧的地层密度差等均有关，因此NVDR_THDR值的大小并不能完全反映断裂的埋深、规模等特征。

图5 孤立形体组合模型位场边缘识别结果Fig.5 The edge recognition results of the combination model

图6 鄂尔多斯盆地北部重力异常及位场边缘识别结果Fig.6 Gravity anomaly in northern Ordos Basin and the edge recognition results

图7 鄂尔多斯盆地北部位场边缘识别结果与地震剖面对比(地震剖面据文献[48])Fig.7 The comparison between the edge recognition results and the seismic profile in northern Ordos Basin (The seismic profile is from reference [48])

4 结论

为增强位场边缘识别方法对于深部小规模地质体的识别能力，本文提出了一种基于界面反演增强技术，首先利用密度界面反演方法对重力异常进行处理，之后将NVDR_THDR技术作为一种边缘提取和增强技术对密度界面反演结果进行提取，其比直接对重力异常进行NVDR_THDR计算的效果更佳。裂陷盆地模型试验证明了该方法能较好地用于盆地基底断裂的识别，密度差的大小几乎不影响最终的边缘识别结果。将该方法应用于孤立形体组合模型试验，处理结果也与理论模型较为吻合，与其他增强方法相比，这一方法具有一定的优势。鄂尔多斯盆地实际资料处理结果也表明了本文方法可以更好地识别盆地基底小规模断裂，实际应用效果较好。

基于密度界面增强的位场边缘识别方法可以直接用于重力异常， 但不能直接用于磁力异常处理。对于磁异常，可将其转化为磁源重力异常进行边缘识别，或将这一思路推广至磁力异常处理之中，将密度界面反演改为磁性界面反演，之后利用NVDR_THDR进行边缘提取和增强。