基于Tanimoto相似的重力异常相关成像方法研究

doi:10.11720/wtyht.2019.0163

基于Tanimoto相似的重力异常相关成像方法研究

王颖舜, 李军, 简兴祥

成都理工大学地球物理学院,四川成都 610059

Research on gravity anomaly correlation imaging method based on Tanimoto similarity

WANG Ying-Shun, LI Jun, JIAN Xing-Xiang

College of Geophysics, Chengdu University of Technology, Chengdu 610059, China

责任编辑: 王萌

收稿日期: 2019-03-25 修回日期: 2019-05-23 网络出版日期: 2019-12-20

基金资助:

国家重点基础研究发展计划“深部矿产地球物理三维建模与预测”项目. 2017YFC0601504

Received: 2019-03-25 Revised: 2019-05-23 Online: 2019-12-20

作者简介 About authors

王颖舜(1997-),男,本科,目前在读于成都理工大学勘查技术与工程专业。Email:408692063@qq.com 。

摘要

相关成像方法主要用于求解异常源空间赋存位置,但其成像结果不理想,在深度方向上存在发散现象。由于正演对成像的验证性作用,笔者设计了简单模型、多个异常源叠加模型,对传统基于Pearson相似算法、余弦相似算法的相关成像方法进行研究,指出传统相关成像中利用Pearson相似计算重力异常时出现假异常的现象,同时提出了基于Tanimoto相似的相关成像,并针对该算法改进了核函数,对多源异常采取窗口化的异常分离方式。研究表明笔者提出的方法具有异常源纵、横向分辨率较传统相关成像方法有所提高,异常源形态勾勒准确的特点,为重力资料的解释提供了有力信息。

关键词： 重力异常 ; Tanimoto相似 ; 核函数 ; 虚假异常源 ; 窗口化

Abstract

The correlation imaging method is mainly used to solve the location of anomalous source space, but its imaging results are not ideal, and there is divergence in the longitudinal direction. Due to the confirmatory effect of forward modeling on imaging, this paper designs a simple model and multiple anomaly source overlay models. The traditional imaging methods based on Pearson similarity algorithm and cosine similarity algorithm are studied. It is pointed out that Pearson similarity is used to calculate gravity anomalies in traditional correlation imaging. The phenomenon of false anomalies occurs, the related imaging based on Tanimoto similarity is proposed, and the kernel function is improved for the algorithm. The research shows that the method proposed in this paper has an improved vertical and horizontal resolution of the source compared with the traditional correlation imaging method, and the abnormal source form outlines the accurate features, which provides powerful information for the interpretation of gravity data.

Keywords： gravity anomaly ; Tanimoto similar ; kernel function ; false anomaly source ; windowing

PDF (3952KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

王颖舜, 李军, 简兴祥. 基于Tanimoto相似的重力异常相关成像方法研究. 物探与化探[J], 2019, 43(6): 1350-1357 doi:10.11720/wtyht.2019.0163

WANG Ying-Shun, LI Jun, JIAN Xing-Xiang. Research on gravity anomaly correlation imaging method based on Tanimoto similarity. Geophysical and Geochemical Exploration[J], 2019, 43(6): 1350-1357 doi:10.11720/wtyht.2019.0163

0 引言

相关成像作为一种重力异常几何形态反演方法,由国外学者Patella D^[1]于1997年首次提出,最初用于自然电场异常的解释,随后被Mauriello P和Patella D推广到了重磁领域^[2],目前被广泛应用到重力资料解释中。该方法计算速度快,能描述异常源空间位置,一度成为研究热点。

传统的相关成像方法存在深度方向分辨力低的问题^[3,4],并且在多个叠加异常源存在时基本失效,其成像结果不能刻画异常源赋存状态。学者们对此展开研究,传统的相关成像方法利用的是余弦相似^[5,6],这种相似算法在相关成像的应用范围最广,该方法对于单一异常源横向分辨较好,但深度方向发散,对多异常源成像效果差。在此基础上,发展了Pearson相似^[7,8],该算法对异常源深度方向分辨率有所提高,但存在虚假异常源这一现象,为解释工作带来困难。也有学者从相关成像中的核函数入手,对相关成像的核函数进行改造,这种方法能够在一定程度上提高深度方向、横向分辨率,但会改变对源位置的聚焦,改造的方式主要有对核函数求垂向导数^[9,10]、乘以加权因子^[11,12]。由于相关成像本质上可以看成一种滤波运算,所以有学者采取滤波^[13,14]、延拓^[15]、截取窗口等方法,滤波能够有效提高对源位置的分辨率,将浅源与深源分开,但滤波窗口的选择较为重要;延拓同样能够对异常源的位置求解较好,但是计算量有所增大,需将上下延拓数据一起进行处理;截取窗口对于深度向无叠加异常源时效果较好,但窗口大小以及窗口位置依赖于主观判断,不利于自动化反演。还有学者利用相关成像结果进行物性反演^[16,17],以重力异常几何形态的信息减少反演的多解性,但目前的相关算法对源的描述都有所偏差,深度方向基本都存在发散的现象,所以一直以来没有得到大范围的应用。

针对目前方法普遍存在对异常源深度方向发散的情况,笔者提出一种基于Tanimoto相似的成像算法,并在该相似下改进了核函数。通过对简单模型的正演结果进行传统的相关成像和新的相关成像,发现Pearson相似下的相关成像会出现虚假异常源,而新的相关成像方法的确能提高其对异常源位置的精确度,克服了相同模型下传统相关成像方法深度方向发散的问题,且对异常源形态勾勒更准确,具有很大的潜力。

1 相关成像方法

1.1 传统相关成像方法以及虚假异常的提出

在笛卡儿坐标系下的平面在基准面上,z轴取垂直向下为正,假设测区地下任意一点的坐标为(ξ,η,ζ),剩余密度为σ,则它在测区上任意点(x,y,z)处的重力位为V(x,y,z),该重力位^[18,19]可表示为:

(1)

V (x, y, z) = G \underset{v}{\int \int \int} \frac{σ d ξ d η d ζ}{r} d v,

其中:r= $\sqrt[]{{(ξ - x)}^{2} + {(η - y)}^{2} + {(ζ - z)}^{2}}$ 。该测点的重力异常表示为:

(2)

\begin{array}{r} Δ g & = \frac{\partial V (x, y, z)}{\partial z} \\ = G \underset{v}{\int \int \int} \frac{σ (ξ - z) d ξ d η d ζ}{r^{3}} d v 。 \end{array}

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 直角坐标系下异常源计算

Fig.1 Abnormal source calculation in Cartesian coordinate system,

将地下空间按图2网格化后取球形近似得其中任一网格P(x₀,y₀,z₀)产生的重力异常:

(3)

Δ g_{ij} (x_{0}, y_{0}, z_{0}) = \frac{Gσv z_{0}}{[(x_{0} {- x)}^{2} + (y_{0} {- y)}^{2} + z_{0}^{2}]^{\frac{3}{2}}},

其中:i代表剖分网格中的第i个网格,j代表第j个重力异常数据,σ代表该网格的剩余密度,v代表网格的体积,在某一网格中我们认为是均匀的。其中网格剖分数为m,重力异常数据量为n。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 相关网格化计算

Fig.2 related grid computing

则求解得正演网格数据可表示为:

(4)

Δ g = [\begin{array}{l} Δ g_{11} & \dots & Δ g_{1 n} \\ ︙ & ⋱ & ︙ \\ Δ g_{m 1} & \dots & Δ g_{mn} \end{array}],

同时重力异常数据可表示为向量d:

(5)

d = (d_{1}, d_{2}, \dots, d_{n - 1}, d_{n})^{T},

将计算出的异常值与重力异常分别归一化后求两者的相似系数,其中Pearson相似计算式为:

(6)

C_{i} = \frac{Cov (Δ g_{i}, d)}{\sqrt[]{Var [Δ g_{i}] Var [d]}},

余弦相似计算式为:

(7)

C_{i} = \frac{Δ g_{i} \cdot d}{\sqrt[]{‖ Δ g_{i} ‖_{2} \cdot ‖ d ‖_{2}}} 。

这样,计算所有网格对重力异常数据的相关系数C,即实现对重力异常的相关成像,成像结果表征地下密度异常地质体的空间赋存状态,勾勒地下密度异常体大致形态。

对于式(7),由许瓦兹不等式易知,(Δg_i·d)²≤‖Δg_i‖₂·‖d‖₂,即对于余弦相似C_i的取值范围为-1~1,同理对于式(2),不过是对去中心化的Δg_i、d求余弦相似,所以其取值范围也是-1~1。

相关系数的物理意义是:相关系数|C_i|大者则代表异常源在该网格处存在的可能性大,相关系数|C_i|小者则代表异常源在该网格处存在的可能性小。

将式(3)与式(6)、(7)联立可以发现Gσν这一项在计算过程中被消去了,于是为了简化计算,可以将式(3)的计算改为

(8)

ψ_{ij} (x_{0}, y_{0}, z_{0}) = \frac{z_{0}}{[(x_{0} {- x)}^{2} + (y_{0} {- y)}^{2} + z_{0}^{2}]^{3 / 2}},

式(8)的意义为,重力异常相关成像仅由所选网格在地表产生重力异常的曲线形态决定,与所选网格的密度无关,该函数在相关成像中被称为核函数,记为ψ(x₀,y₀,z₀)。导数相关成像方法^[20,21]的思路是通过对网格的重力异常核函数ψ(x₀,y₀,z₀)求各方向导数以此对核函数进行改造,通过不同方向上的信息求取地下空间异常源分布。

其中应用最广范的是对z方向求导,求导后的核函数形式为:

(9)

\begin{array}{r} ψ (x_{0}, y_{0}, z_{0}) & = \frac{\partial Δ g_{ij} (x_{0}, y_{0}, z_{0})}{\partial z} \\ = \frac{2 z_{0}^{2} - (x_{0} {- x)}^{2}}{[(x_{0} {- x)}^{2} + (y_{0} {- y)}^{2} + z_{0}^{2}]^{5 / 2}} 。 \end{array}

通过一个简单模型来说明传统相关算法成像方法的缺点。计算一个y方向无限延伸的矩形二度

体,其剩余密度为1 g/cm³,横向范围为-5~5 m,深度向范围为10~20 m,并附加2%的高斯噪声以模拟真实情况。正演结果如图3,模型空间及成像结果见图4。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 模型正演结果

Fig.3 Model forward results

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 二维模型及其相关成像结果

Fig.4 2D model and its imaging results

从图中4b可以看出由Patella D先生提出的相关成像对异常源横向范围划分较准确,但是深度向上有下移的趋势,且深度向衰减较慢,对异常源位置刻画不准确。图4c利用的是沿深度方向求导的核函数进行相关计算,这一方法由王绪本^[9]提出,原本用于大地电磁中。可以看到在此模型中这一方法效果并不好,深度向发散情况更严重,且完全偏离了异常源位置。

图4d中利用了Pearson相似,可以看到基于Pearson相似的相关成像方法对异常源位置刻画比传统余弦相似更为准确,关键在于其深度向发散大为改善。但笔者认为该方法存在缺陷。图5c中pearson相似的剖面图中在两侧出现了两个负概率,由于概率值在无穷远处是趋于0的,所以这里容易被错误地认为是一个具有负密度的异常体,所以笔者认为基于Pearson相似的相关成像会给解释工作带来困难。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 成像结果及其剖面

Fig.5 Imaging results and profiles

1.2 基于Tanimoto相似的相关成像

相关成像都是一种相似算法,受该思路的启发,笔者通过查阅文献以及模型验证,其深度向分辨率受相似计算算法影响较大,所以改进了原有的相似计算方法,提出了基于Tanimoto相似的相关成像方法。

Tanimoto相似计算式为:

(10)

C_{i} = \frac{Δ g_{i} \cdot d}{‖ Δ g_{i} ‖_{2} + ‖ d ‖_{2} - Δ g_{i} \cdot d},

由均值不等式可知‖Δg_i‖₂+‖d‖₂≥2Δg_i·d,带入式(10)有C∈[-1,1]。

同样对上述简单模型利用Tanimoto相似进行相关成像,其结果如图5。

从图5a成像结果中大概率圈闭中心与模型相对位置可以看出,这样的改进虽然能够提高深度向分辨率,但是其对异常源的位置刻画能力又有所下降。受对核函数的导数操作是改变位置分辨能力的有效方法的启发,笔者通过将核函数改造为:

(11)

ψ (x_{0}, y_{0}, z_{0}) = \frac{z_{0}}{[(x_{0} {- x)}^{2} + (y_{0} {- y)}^{2} + z_{0}^{2}]},

同样对简单模型进行成像,其成像结果如上图5b,这时的相关成像所得异常源位置与形态与相关概率圈闭与异常源的形态、位置基本一致,说明在该模型下有较高的成像精度。

为了更深入讨论上述5种方法的优劣,提取了上述5种方法位于模型中心的横向深度向两条剖面线分别是x=0 m;Depth=15 m。其结果如图5c、d。从Depth=15 m剖面可以看出,pearson相似在两端出现了-0.3左右的负概率,在无穷远处会趋于零,这表明了假异常的产生。对于横向,各方法分辨率相差不大。对于深度向,可以明显看到余弦相似以及余弦相似导数的发散性。对于改造的Tanimoto相似,其概率峰值位于异常范围10~20 m之间,且衰减速率最快,验证了该方法对深度向分辨率高。

2 复杂模型相关成像

2.1 二维复杂模型相关成像

为了验证该方法的实用性,设计了多个异常叠加的复杂模型并在正演中同样加入高斯噪声,模型参数及模型示意图见表1和图6。

表1 多异常源参数

Table 1 Multiple abnormal source parameters

参数	模型A	模型B	模型C
左顶点坐标(x,z)	(-17,2)	(15,2)	(-5,20)
二度体宽度/m	2	2	10
二度体厚度/m	2	2	10
剩余密度/(g·cm^-3)	0.21	0.21	0.08

新窗口打开| 下载CSV

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 复杂模型成像结果

Fig.6 Complex model imaging results

对于复杂模型而言,其重力异常叠加在一起不易分离,而相关成像又是对整条测线进行计算,导致相关成像方法对多源异常效果较差。这里我们采取马国庆^[7]提出的基于窗口数据的相关成像,并在该方法上做出改进,对每个窗口的大概率圈闭进行提取、扩边、联合成像,以达到排除叠加异常源影响,更准确的进行成像。其计算流程框图见图7。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 提取子域的Tanimoto相关成像算法流程

(k为阈值概率,在本算例中取0.8)

Fig.7 Extracting sub-domain Tanimoto related imaging algorithm flow

(Where k is the threshold probability, which is 0.8 in this example)

对表1中所述模型进行正演,并加入高斯噪声后,利用图7所示方法计算了改进核函数后的基于Tanimoto相似的相关系数,并且对传统的相关算法也进行了计算,其结果如图6所示。

从图6d中可以直观的看到,利用Tanimoto相似算法进行相关成像与该模型十分契合,同正演模型形态、位置基本一致,而传统相关成像(如图6c)深度方向更为发散,其深度分辨能力与源位置分辨能力均弱于本文所提出的方法。

2.2 三维多源模型相关成像

对于三维模型,将二维模型推广到三维中。模型参数见表2,模型空间展布如图8a。对表2中模型进行正演计算,利用本文方法对该模型进行相关成像,如图8和图9。

表2 三维模型参数

Table 2 3D model parameters

异常体	x方向位置/m	y方向位置/m	z方向位置/m	剩余密度/(g·cm^-3)
异常体1	15~17	-1~1	2~4	0.21
异常体2	-17~-15	-1~1	2~4	0.21
异常体3	-5~5	-5~5	20~30	0.08

新窗口打开| 下载CSV

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 三维模型及正演剖面

Fig.8 3D model and forward profile

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 三维模型正演及成像结果

a—模型正演结果;b—Depth=2 m剖面成像结果;c—Depth=20 m剖面成像结果;d—x=-16 m成像结果;e—x=0 m剖面成像结果;f—x=16 m剖面成像结果;g—y=-5 m剖面成像结果;h—y=0 m剖面成像结果;i—y=5 m剖面成像结果

Fig.9 3D model forward and imaging results

a—model forward performance; b—Depth=2 m profile imaging result; c—Depth=20 m profile imaging result; d—x=-16 m imaging result; e—x=0 m profile imaging result; f—x=16 m profile imaging results; g—y=-5 m profile imaging results; h—y=0 m profile imaging results; i—y=5 m profile imaging results

在相关成像过程中,依据剖面数据的特点,采取了图8a中的黑色窗口作为数据截取窗口,横向范围为13~19 m以及-13~-19 m,纵向范围为-3~3 m,选择高斯函数为拟合函数对数据进行延拓。对于切割出的大异常体阈值概率k取0.9,小异常体阈值概率取0.6,分别成像,再将数据进行融合。利用Tanimoto相似的相关成像方法,对图8a所示各剖面进行了成像(图9),可以发现其成像精度还是比较高的,由于高斯函数与异常源在地表造成异常并不一致,其反应深度往往更浅。图9e、g、h、i中这些剖面对大异常体3的纵向分辨不够准确。

由于给出的三维模型较为简单易于用窗口截取实现,所以得到了较好的成像效果。但对于难以分离的复杂叠加异常,人工截取窗口不易实现,这也是窗口截取这一方法的局限。

3 结论

1) 笔者从重力异常几何形态的成像理论出发,针对重力异常空间赋存状态成像方法存在的问题,提出了基于Tanimoto相似的重力异常相关成像方法,并结合窗口化的方法排除多源异常的互相干扰,设计优化了相应的算法。

2) 通过设计模型进行实验,应用本文所提出的方法,和传统的方法进行比较,可以看出本文方法的优越性,证明本方法的可行性。

3) 笔者所提出的多异常源分离方法采用了窗口选取这一思想,如何将窗口选取这一思想对复杂模型三维相关成像进行应用还有困难,但对于异常叠加效应不强,单一源特点仍保留的异常能取得较好的成像效果。

4) 本文所用模型横纵向比例都较为接近,而对于横纵向相差较大时,由于场源的等效性,本文所提出方法对这类模型分辨率将会下降。如何对横纵比相差大的模型也进行准确度较高的相关成像还需要努力。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Patella

Introduction to ground surface self-potential tomography

[J] Geophys. Prospect, 1997,45(4):653-681.