CN116774303B - 一种场源边界定位方法、装置及计算机可读存储介质 - Google Patents
一种场源边界定位方法、装置及计算机可读存储介质 Download PDFInfo
- Publication number
- CN116774303B CN116774303B CN202310760291.3A CN202310760291A CN116774303B CN 116774303 B CN116774303 B CN 116774303B CN 202310760291 A CN202310760291 A CN 202310760291A CN 116774303 B CN116774303 B CN 116774303B
- Authority
- CN
- China
- Prior art keywords
- boundary
- gradient tensor
- gravity gradient
- enhancement parameter
- derivative
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 79
- 238000003860 storage Methods 0.000 title claims abstract description 11
- 230000005484 gravity Effects 0.000 claims abstract description 105
- 238000004458 analytical method Methods 0.000 claims abstract description 49
- 230000002159 abnormal effect Effects 0.000 claims description 38
- 238000009933 burial Methods 0.000 claims description 7
- 230000004807 localization Effects 0.000 claims description 2
- 230000000694 effects Effects 0.000 abstract description 16
- 238000010586 diagram Methods 0.000 description 19
- 238000004364 calculation method Methods 0.000 description 16
- 230000005856 abnormality Effects 0.000 description 8
- 238000005516 engineering process Methods 0.000 description 8
- 238000004590 computer program Methods 0.000 description 7
- 230000006870 function Effects 0.000 description 7
- 230000008569 process Effects 0.000 description 7
- 238000012545 processing Methods 0.000 description 6
- 239000011435 rock Substances 0.000 description 5
- ATJFFYVFTNAWJD-UHFFFAOYSA-N Tin Chemical compound [Sn] ATJFFYVFTNAWJD-UHFFFAOYSA-N 0.000 description 4
- 230000008859 change Effects 0.000 description 4
- 238000005259 measurement Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 239000011159 matrix material Substances 0.000 description 3
- 238000010606 normalization Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 230000001133 acceleration Effects 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000009826 distribution Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000000704 physical effect Effects 0.000 description 2
- 238000012827 research and development Methods 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- 230000008901 benefit Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 230000009545 invasion Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开一种场源边界定位方法、装置及计算机可读存储介质,所述方法包括:在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位;计算重力场和/或重力梯度张量;所述重力场是所述引力位在笛卡尔坐标系下的一阶导数;所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数;所述重力梯度张量具有多个独立分量;获取所述重力场的解析信号模,并基于所述解析信号模对所述三维地质体进行边界定位,和/或获取所述重力场的重力梯度张量的方向解析信号的振幅,并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。通过本发明的方法对目标体的边界识别效果增强,且对三维地质体的识别和定位效果较好。
Description
技术领域
本发明涉及地球物理领域,尤其涉及一种场源边界定位方法及装置及计算机可读存储介质。
背景技术
重力梯度张量测量在国外已开展了十余年,并取得了一定的成果,然而这项测量技术及相关数据解释技术在我国的研究仍然较少。国内在全张量重力梯度仪的研发过程中还处于初级阶段。目前,有多家单位进行相关研究,主要是进行相关重力梯度仪的理论设计、误差分析和模型仿真,都处于实验室研发阶段,尚不能达到勘探要求。重力梯度数据的解释工作,可以进行单独分量的解释,也可以联合多个分量进行解释。在解释中,选择哪些分量或者如何结合不同的分量进行解释,便成为张量梯度数据解释工作所需要解决的重要问题之一。基于系统的发展,梯度测量数据的评估、分析也逐渐形成一套完善的理论体系。梯度数据和重力异常数据对比评估,梯度数据噪音的去除,以及梯度数据的处理方法等问题均成为重要的研究课题。与此同时,解释方法技术也被大量建立和发展,且和重力异常数据的解释方法相依托,形成了目前的重力异常及其梯度数据的解释现状。重力数据的解释工作,通常需要计算出重力场源体的三种基本信息:水平位置,深度范围,物性差异。水平位置的确定,常用到的解释方法即边界识别和增强方法,这种方法可以有效地突显出地质体边缘及构造的水平位置。深度计算通常需要根据异常的特征点及地质体形体来进行估算。物性参数的计算,有常用的物性反演计算,或者可采用快速成像算法快速给出与物性相关的量。重力异常边界识别技术利用了地质体在断裂位置、接触带附近、围岩周边位置的密度变化率较大,从而引起重力异常有明显变化的特点,通过异常的二次特征对变化率较大的位置进行判读,从而达到地质体边界识别(水平位置或在平面上的投影位置)的目的。
但目前基于张量数据进行边界识别的方法存在一些问题:目标体的边界识别效果不佳,对三维地质体的识别和定位效果并不好。因此,发明一种新的场源边界定位方法技术势在必行。
发明内容
本发明提供了一种场源边界定位方法及装置及计算机可读存储介质,用于解决现有技术中基于张量数据进行边界识别对目标体的边界识别效果不佳,对三维地质体的识别和定位效果并不好的问题。
本发明中的一种场源边界定位方法,所述方法包括:
在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位;
计算重力场和/或重力梯度张量;所述重力场是所述引力位在笛卡尔坐标系下的一阶导数;所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数;所述重力梯度张量具有多个独立分量;
获取所述重力场的解析信号模,并基于所述解析信号模对所述三维地质体进行边界定位,和/或获取所述重力场的重力梯度张量的方向解析信号的振幅,并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。
可选的,所述三维地质体为直立六面体模型。
可选的,所述方法还包括:计算所述重力梯度张量的一阶混合导数以对所述三维地质体进行边界定位。
可选的,所述方法还包括:根据下列公式获取第一边界增强参数|ED|:
其中,|Axz|为所述重力梯度张量在x方向解析信号Ax的垂向一阶导数模,|Ayz|为所述重力梯度张量在y方向解析信号Ay的一阶导数模;
根据所述第一边界增强参数|ED|对异常体进行边界识别。
可选的,所述方法还包括:根据下列公式获取第二边界增强参数|ED|_z:
所述第二边界增强参数|ED|_z为采用重力梯度张量在z方向解析信号的总水平导数;
根据所述第二边界增强参数|ED|_z对异常体进行边界识别。
可选的,所述方法还包括:根据下列公式获取第三边界增强参数|ED|_t:
根据所述第三边界增强参数|ED|_t对异常体进行边界识别。
可选的,所述第一边界增强参数、第二边界增强参数以及第三边界增强参数的极大值位置均对应所述异常体的边界。
可选的,所述方法还包括:根据下列公式获取归一化边界增强参数N|ED|:
其中,|ED|Wmax为包含计算点的移动窗口内各点第一边界增强参数中的最大值。
本发明中的一种场源边界定位装置,所述装置包括:
第一获取单元,用于在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位;
第二获取单元,用于计算重力场和/或重力梯度张量;所述重力场是所述引力位在笛卡尔坐标系下的一阶导数;所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数;所述重力梯度张量具有多个独立分量;
第三获取单元,用于获取所述重力场的解析信号模,并基于所述解析信号模对所述三维地质体进行边界定位,和/或获取所述重力场的重力梯度张量的方向解析信号的振幅,并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。
本发明中的一种计算机可读存储介质,所述计算机可读存储介质存储有一个或者多个程序,所述一个或者多个程序可被一个或者多个处理器执行,以实现如上任意一项所述的场源边界定位方法的步骤。
附图说明
图1是本发明实施例中一种场源边界定位方法的流程示意图;
图2是本发明实施例中笛卡尔坐标系下的长方体模型及地下空间剖分;
图3是本发明实施例中直立六面体示意图;
图4是本发明实施例中模型一的引力位一阶导数示意图;
图5是本发明实施例中模型一的重力梯度张量异常示意图;
图6是本发明实施例中模型一的部分重力梯度张量混合一阶导数示意图;
图7是本发明实施例中模型一的重力梯度张量解析信号模示意图;
图8是本发明实施例中模型一的重力梯度张量解析信号模的一阶导数示意图;
图9是本发明实施例中模型一|ED|、|ED|z、|ED|t与THD、ASMz、Tilt_ASM对比示意图;
图10是本发明实施例中模型二|ED|、|ED|z、|ED|t与THD、ASMz对比示意图;
图11是本发明实施例中模型三|ED|、|ED|z、|ED|t与THD、ASMz对比示意图;
图12是本发明实施例中模型四|ED|、|ED|-z、|ED|-t边界增强结果示意图;
图13是本发明实施例中模型四N|ED|、N|ED|-z、N|ED|-t边界增强结果对比示意图;
图14是本发明实施例中模型五Tilt、Tilt_Az和Tilt_Az′反边界增强结果示意图;
图15是本发明实施例中模型五Tilt_Azz与Tilt_Azz‘边界增强结果对比示意图;
图16是本发明实施例中模型六Tilt_Azz与Tilt_Azz‘边界增强结果对比示意图;
图17是本发明实施例中一种场源边界定位的结构示意图。
具体实施方式
下面结合附图和实施例对本发明作进一步的详细说明。可以理解的是,此处所描述的具体实施例仅仅用于解释本发明,而非对本发明的限定。另外还需要说明的是,为了便于描述,附图中仅示出了与本发明相关的部分而非全部结构。
应理解,在本文的各种实施例中,上述各过程的序号的大小并不意味着执行顺序的先后,各过程的执行顺序应以其功能和内在逻辑确定,而不应对本文实施例的实施过程构成任何限定。
利用重磁位场资料进行地质体的识别推断,如岩体范围、断裂的位置等场源边界的识别,是进一步深化研究的基础。边界识别或边界定位常指利用一系列的点或线段明确标定如断裂构造线或地质体的水平投影线位置的表达方式。但由于受到场源埋深、异常叠加以及噪声等多种因素的影响,直接进行场源的边界识别往往存在很大的误差,需要开展边界特征增强,以便于进行场源边界识别工作的开展。
重力梯度张量数据相对于传统的重力测量数据有着更高的频率信息,能更加准确、细致地研究地球内部构造、矿产资源分布等信息。本发明具体实施例是基于重力场和/或重力梯度张量重解析信号及其导数,并在此基础上,经过再求导、归一化等方法实现场源体边界增强,以提高边界定位的准确度。
本发明具体实施例提供一种场源边界定位方法,如图1所示,所述方法包括:
步骤100,在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位。具体的,在笛卡尔坐标系下,设地面某一点O为坐标原点,根据万有引力定律,一个剩余质量为m,剩余密度为ρ的地质体在其外部空间任意一点产生的引力位为:
其中,G=6.67×10-11m3/(kg·s)2为万有引力常量,为任意一点P(x,y,z)到场源点Q(ξ,η,ζ)的矢量半径,r=[(ξ-x)2+(η-y)2+(ζ-z)2]1/2,如图2所示。
步骤200,计算重力场和/或重力梯度张量;所述重力场是所述引力位在笛卡尔坐标系下的一阶导数;所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数;所述重力梯度张量具有多个独立分量。具体的实施例中,步骤200可以包括不同的实施方式。在第一种实施方式中,只计算重力场。其中,重力场是引力位在笛卡尔坐标系x,y,z三个方向上的一阶导数,即重力加速度矢量/>可以得到:
在另一具体的实施方式中,还需计算重力梯度张量,其中:重力梯度张量为地质体引力位的二阶导数,即为重力加速度矢量在笛卡尔坐标系x,y,z三个方向的一阶导数,也即重力梯度:
写成矩阵形式,即
其中:V为地质体的引力位,为引力位在x方向的二阶导数,即gxx;
为引力位在x、y方向的二阶导数,即gxy,其他参数的定义可参照如上。
步骤300,获取所述重力场的解析信号模,并基于所述解析信号模对所述三维地质体进行边界定位,和/或获取所述重力场的重力梯度张量的方向解析信号的振幅,并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。
在一具体的实施方式中,可以直接通过重力场进行边界定位,但定位边界较为模糊。在另一具体的实施方式中,还可以通过解析信号模对所述三维地质体进行边界定位。具体的,由于重力场是无源场,因此其旋度和散度均为零,即:
所以重力梯度张量为对称矩阵,并且迹为零:
gxx+gxy+gxz=0
gxy=gyx,gxz=gzx,gyz=gzy
因此,重力梯度张量的9个分量中,只有5个独立分量。
二维地质体位场φ(x)的解析信号A(x)可以写成复函数:
其中和/>是希尔伯特变换对。二维解析信号的振幅表示为:
将二维解析信号归纳到三维,并且任何位场的希尔伯特变换满足柯西—黎曼关系。将水平面观测的位场φ(x,y)的解析信号的定义扩展到三维,如下:
则A(x,y)的振幅表示为:
在三维情况下,重力梯度张量(GGT)第三列元素是第一和第二列元素的希尔伯特变换对。然后,可以定义在x、y和z方向的方向解析信号,其矩阵形式表示为:
因此,GGT的方向解析信号的振幅表示为:
对于多个相互影响的异常体,利用解析信号振幅的导数可以实现异常体分离。x、y、z方向上的方向性解析信号的导数可以表示为:
其中,α是x、y、z。
本发明具体实施例所述的方法从重力场和重力梯度理论出发,对重力场和重力梯度张量数据的解释技术进行探索性研究,基于重力场和重力梯度张量解析信号及其方向导数,提出了新的边界定位方法在目标体的边界识别上效果更好。对今后张量技术实际应用过程中的数据处理与解释提供一定的理论基础。
本发明具体实施例所述的一种场源边界定位方法,较佳的,所述三维地质体为直立六面体模型。具体的,直立六面体模型是一种非常典型和常见的三度体地质模型,在地球物理的实际工作中很常见的侵入岩、岩床、水平板、垂直板、岩脉、岩墙以及基底变质岩系都可以用直立六面体进行理论模拟。另外,通过将地下介质剖分成足够小的直立六面体,通过多个小直立六面体的组合,可以模拟三维复杂模型,实现三维复杂模型的正反演。因此对直立六面体模型进行模拟试算具有很强的实际意义。
本发明具体实施例所述的一种场源边界定位方法,较佳的,所述方法还包括:计算所述重力梯度张量的一阶混合导数以对所述三维地质体进行边界定位。
具体的实施例中,为了识别三维地质体中存在的异常体,建立如图3所示的笛卡尔坐标系,设直立六面体两角点坐标分别为(ξ1,η1,ζ1)和(ξ2,η2,ζ2),剩余密度为σ,在地面上任意观测点P(x,y,z)的重力异常、重力异常一阶导数及二阶导数的计算如下。
直立六面体引力位的一阶导数计算公式,一共有3个分量,如下:
通过分析可以看出,当上面三个式子中反正切函数的分子分母都为0时,这个分式无意义,此时会出现“奇点”。因此重新推导gx和gy,得到了不含解析奇点的计算公式:
建立单个异常体的模型一,异常体长宽高为500×500×500m3,剩余密度100kg/m3,顶面埋深20m,中心点坐标(0,0,270)。模型一的引力位一阶导数如图4所示。
直立六面体引力位的二阶导数计算公式,即直立六面体重力梯度张量,一共有9个分量,分别如下:
gzx=gxz,gzy=gyz,gyx=gxy
式中
同样该公式中gxx,gxy,gxz也存在解析奇点。因此推导出无解析奇点离散计算公式:
图5是模型一的重力梯度张量异常,可以看出gxx能有效地增强地质体的x方向边界,gyy能增强异常的y方向边界,gxy显示了立方体四个角点异常的峰值,能描述地质体四个方位角点信息,gzx的极大值可以用来识别x方向边界,gzy的极大值能识别地质体的y方向边界,不同张量异常可以反映地质体不同边界信息。
通过对引力位的二阶导数求导,即对重力梯度张量的混合一阶导数,可以得到引力位三阶混合导数,一共17个分量,分别为:
gxyz=gyxz=gxzy=gzxy
图6所示为模型一的部分重力梯度张量混合一阶导数,可以看出在各个方向上其具有更高的分辨率。
本发明一具体的实施例中,直立六面体的重力梯度张量在x、y、z方向上的方向解析信号可以直接利用重力梯度张量的6个分量计算获得。
图7为模型一的x、y、z方向解析信号模。可以看出,解析信号模的优点是不会出现假值,只有在边界位置出现高值;但是边界增强的横向分辨率不足,边界位置的过渡带较宽。只识别单个方向的边界,比如只识别南北边界,东西边界幅值非常弱或几乎没有增强的效果。因此本发明具体实施例所述的一种场源边界定位方法中,可以利用引力位的二阶导数、三阶导数、GGT解析信号计算获得,一共9个分量,分别为:
图8为模型一的解析信号一阶导数。相比解析信号模,垂向导数提高了边界极值的横向分辨率,边界过渡带缩小了很多,非边界区域也没有假值,但与解析信号模类似,只识别单个方向的边界,另一个方向的边界没有识别或者识别结果较为模糊。
基于图8分析重力梯度张量解析信号一阶导数的特征,可以看出Axz反应的x方向边界更清晰,Ayz反应的y方向边界更清晰,因此本发明具体实施例所述的一种场源边界定位方法,较佳的,所述方法还包括:根据下列公式获取第一边界增强参数|ED|:
其中,|Axz|为所述重力梯度张量在x方向解析信号Ax的垂向一阶导数模,|Ayz|为所述重力梯度张量在y方向解析信号Ay的一阶导数模;
根据所述第一边界增强参数|ED|对异常体进行边界识别。
具体的,第一边界增强参数|ED|的定义类似于传统解析信号的垂向导数ASz,其实质是重力梯度张量x方向解析信号Ax的垂向一阶导数模|Axz|与y方向解析信号Ay的一阶导数模|Ayz|的组合,通过第一边界增强参数|ED|可以增加对异常体的边界识别。
较佳的实施例中,除了Axz、Ayz具有较好的边界增强效果外,Azx、Azy也具有较好的边界增强效果,据此又定义了第二边界增强参数|ED|_z。较佳的,所述方法还包括:根据下列公式获取第二边界增强参数|ED|_z:
所述第二边界增强参数|ED|_z为采用重力梯度张量在z方向解析信号的总水平导数;
根据所述第二边界增强参数|ED|_z对异常体进行边界识别。
具体的,|ED|_z的定义则类似于传统的总水平导数,即采用重力梯度张量z方向的解析信号计算总水平导数。
本发明具体实施例所述的一种场源边界定位方法,联合|ED|、|ED|-z两个参量又定义了第三边界增强参数|ED|-t。较佳的,所述方法还包括:根据下列公式获取第三边界增强参数|ED|_t:
根据所述第三边界增强参数|ED|_t对异常体进行边界识别。
本发明具体实施例所述的一种场源边界定位方法,较佳的,所述第一边界增强参数、第二边界增强参数以及第三边界增强参数的极大值位置均对应所述异常体的边界。
上述|ED|、|ED|-z、|ED|-t均为边界增强参数,其极大值位置对应异常体的边界。采用四方位识别极大值的方法可以自动圈定异常体的边界,实现边界识别。
对上述三个边界识别参数与传统效果较好的THD、ASMz、Tilt_ASM方法进行对比,如图9所示。需要注意的是,理论上ASMz与Azz相同。另外,Tilt_ASM采用gz参数进行多次频域求导获得,即传统Tilt_ASM方法。|ED|、|ED|-z、|ED|-t、ASMz异常增强的极值范围较窄,THD、Tilt_ASM极值范围较宽;|ED|、|ED|-z、|ED|-t、ASMz、THD边界识别结果都能与真实模型吻合,|ED|-g在异常内部出现了虚假异常、Tilt_ASM由于经过多次求导计算,开展边界识别的假异常更多。因此,从边界识别的角度,传统的Tilt_ASM效果相对其它五种方法较差。
采用模型二再次进行比较,如图10所示。模型二与模型一相似,唯一的变化为深度变为350~850m。从图10可以看出,|ED|、ASMz已经不能识别边界,|ED|-t边界范围变小,|ED|-z、|ED|边界结果相似,范围更宽。
模型三为组合模型,该模型包含两个直立六面体,浅部异常体200×200×200m3,顶部埋深20m,剩余密度差100kg/m3;深部异常体750×750×750m3,顶部埋深220m,剩余密度差1300kg/m3。从图11可以看出:
(1)THD已经不能识别浅部的异常体,|ED|-t、|ED|-z、|ED|、ASMz均可以反映浅部异常体,其中|ED|-z的异常衬度最大;
(2)|ED|-z的异常形态反映深部异常体最明显,其次为|ED|-t,再次为|ED|,ASMz最差,THD虽然也能反映深部异常体,但是异常范围过宽;
(3)|ED|、ASM识别出的边界比真实模型范围稍小,|ED|-z范围稍大,|ED|-t识别的边界范围与真实模型更加接近。
通过上述三个模型的对比,得到如下结论:
(1)采用GGT解析信号一阶导数组合定义的|ED|、|ED|-z、|ED|-t总体效果好于THD、ASMz、Tilt_ASM;
(2)|ED|-z的定义和边界增强效果类似于HDT,|ED|-t的定义和边界增强效果类似于ASMz;
(3)|ED|浅部异常体的边界反映最好,但是深部异常体的边界反映较差;
(4)|ED|t、|ED|z浅部异常体的边界反映很好,也可以反映深部异常体的边界;
|ED|、|ED|-z、|ED|-t存在的主要问题是对弱异常不敏感,对深部地质体的异常信息缺少足够的反映。
本发明具体实施例所述的一种场源边界定位方法,较佳的,所述方法还包括:根据下列公式获取归一化边界增强参数N|ED|:
其中,|ED|Wmax为包含计算点的移动窗口内各点第一边界增强参数中的最大值。
具体的,在模型四中包含三个直立六面体,尺寸均为40×160m,剩余密度差均为300kg/m3,由西向东顶部埋深依次为10m、25m、40m。从图12可以看出|ED|、|ED|-z、|ED|-t与THD对浅部异常体都可以识别,但是深部两个异常体则很难识别。
因此,本发明具体实施例借鉴归一化总水平导数法,针对|ED|、|ED|-z、|ED|-t对大埋深场源的边界增强不足的特点,利用移动窗口内的最大值归一化计算点的值,并利用极大值标定场源边界,以期增强弱异常,来突出深部异常。归一化后的归一化边界增强参数NED表达式为:
其中:|ED|表示计算点的第一边界增强参数,|ED|Wmax为包含计算点的移动窗口内各点边界增强参数中的最大值,较佳的实施例中,|ED|也可以用|ED|-z、|ED|-t进行替换。
图13为模型四采用归一化的方法计算得到的N|ED|、N|ED|-z、N|ED|-t边界增强结果对比,计算窗口为80×80m。从图13中可以看出:
(1)采用归一化的方法确实能增强深部异常,但是深部异常边界的极值范围较宽;
(2)N|ED|的极值并不反映边界而是反映中心;
(3)N|ED|-z、N|ED|-t的极值反映边界,N|ED|_z的异常衬度更大。
较佳的实施例中,归一化类方法的边界增强效果与窗口的选取有关。
本发明具体实施例中,还提供一种解析信号倾斜角Tilt_ASM的方法采用gz的解析信号ASM替换gz进行边界增强,具体的计算方法为:
传统Tilt_ASM是针对gz参数开展多次频域求导。分析上式,AS就是重力梯度张量数据的解析信号Az。因此,可写为:
上式需要计算Az的x和y方向一阶导数,采用相邻点的有限差分代替求导,差分间距为两倍点距。但差分间距小于1/10~1/100倍点距时,有限差分的误差才可以忽略。本发明实施例中通过正演可以直接计算|Azx|2与|Azy|2,其精度比两倍点距有限差分计算结果高。
由于重力梯度张量数据提供了更加丰富的信息(Ax与Ay),理论上Axx与Ayy分别对x与y方向的边界更加敏感。|Axx|2与|Ayy|2的0值位置对应异常体边界。采用倾斜角方法,本发明实施例提出了重力梯度张量解析信号的Tilt方法,公式为:
模型五包含三个埋深不同、大小完全相同的直立六面体,埋藏深度分别为4km(棱柱体1)、1km(棱柱体2)和2.5km(棱柱体3),所有直立六面体的厚度都为1km,剩余密度为200kg/m3。
图14中的(a)和(b)为模型五Tilt、Tilt_Az边界增强结果与(c)的文献结果对比,从图中可以看出,传统Tilt方法异常体内为正值,异常体外为负值,边界为0值;Tilt_Az则可以根据极大值反映边界位置,且边界清晰;Tilt_Az′反映的边界更加规整,与模型的真实边界吻合更好。但是Tilt_Az和Tilt_Az′在异常体内部和异常体中间位置存在假异常。
解析信号的导数相对于解析信号具有更高的分辨率。利用解析信号的导数,对上两式做进一步的改进为:
另外,根据获得的异常形态及特征值分布情况,可以对Tilt_Azz、Tilt_Azz′进行大于0阈值的峰值提取。以Tilt_Azz为例,有:
图15的(a)和(b)为模型五Tilt_Azz与Tilt_Azz′边界增强结果对比,从图中可以看出,Tilt_Azz与Tilt_Azz′边界增强后的异常更加集中,与模型的真实边界吻合更好。
为了验证该方法对包含正、负异常数据的适用性,将模型五中东北部直立六面体的密度差修改为-200kg/m3,其它参数不变,构建模型六。图16的(a)和(b)为模型六Tilt_Azz与Tilt_Azz′边界增强结果对比,从结果可以看出,Tilt_Azz与Tilt_Azz′对负异常依旧适用。
本发明具体实施例还提供一种场源边界定位装置,如图17所示,所述装置包括:
第一获取单元1701,用于在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位;
第二获取单元1702,用于计算重力场和/或重力梯度张量;所述重力场是所述引力位在笛卡尔坐标系下的一阶导数;所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数;所述重力梯度张量具有多个独立分量;
第三获取单元1703,用于获取所述重力场的解析信号模,并基于所述解析信号模对所述三维地质体进行边界定位,和/或获取所述重力场的重力梯度张量的方向解析信号的振幅,并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。
本发明实施例还提供一种计算机可读存储介质,所述计算机可读存储介质存储有一个或者多个程序,所述一个或者多个程序可被一个或者多个处理器执行,以实现如上具体实施例中的任意一项所述的场源边界定位方法的步骤。
本发明实施例所述的方法通过对重力梯度张量解析信号与欧拉反褶积相结合,进行了地质体定位方法的研究,通过数值模拟进行试算,发现对于单个直立六面体与规则组合模型的边界识别和定位效果较好,对三维地质体具有较高的可行性,以期对今后张量技术实际应用过程中的数据处理与解释提供一定的理论基础。
本领域内的技术人员应明白,本申请的实施例可提供为方法、系统、或计算机程序产品。因此,本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器,使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
前述对本发明的具体示例性实施方案的描述是为了说明和例证的目的。这些描述并非想将本发明限定为所公开的精确形式,并且很显然,根据上述教导,可以进行很多改变和变化。对示例性实施例进行选择和描述的目的在于解释本发明的特定原理及其实际应用,从而使得本领域的技术人员能够实现并利用本发明的各种不同的示例性实施方案以及各种不同的选择和改变。本发明的范围意在由权利要求书及其等同形式所限定。
Claims (7)
1.一种场源边界定位方法,其特征在于,使边界定位结果更加准确,并且同时识别不同埋深的异常体,所述方法包括:
基于笛卡尔坐标系下的重力场和/或重力梯度张量数据,计算对异常体进行边界定位的参数;
其中,所述方法还包括:基于重力梯度张量解析信号模对场源体边界进行识别;
所述方法还包括:根据下列公式获取第一边界增强参数|ED|:
其中,|Axz|为所述重力梯度张量在x方向解析信号Ax的垂向一阶导数模,|Ayz|为所述重力梯度张量在y方向解析信号Ay的垂向一阶导数模;
根据所述第一边界增强参数|ED|对异常体进行边界识别;
其中,所述方法还包括:根据下列公式获取第二边界增强参数|ED|-z:
所述第二边界增强参数|ED|-z为采用重力梯度张量在z方向解析信号的总水平导数;
根据所述第二边界增强参数|ED|-z对异常体进行边界识别;
其中,所述方法还包括:根据下列公式获取第三边界增强参数|ED|_t:
根据所述第三边界增强参数|ED|_t对异常体进行边界识别;
或,
其中,所述方法还包括:基于重力梯度张量解析信号的倾斜角方法以对三维地质体进行边界定位;
其中,所述方法还包括:根据下列公式获取Tilt_Az’:
其中,所述方法还包括:根据下列公式获取Tilt_Azz:
其中,所述方法还包括:根据下列公式获取Tilt_Azz’:
2.根据权利要求1所述的一种场源边界定位方法,其特征在于,所述三维地质体为直立六面体模型。
3.根据权利要求1所述的一种场源边界定位方法,其特征在于,所述方法还包括:计算所述重力梯度张量的混合一阶导数以对所述三维地质体进行边界定位。
4.根据权利要求1所述的一种场源边界定位方法,其特征在于,所述第一边界增强参数、第二边界增强参数以及第三边界增强参数的极大值位置均对应所述异常体的边界。
5.根据权利要求1所述的一种场源边界定位方法,其特征在于,所述方法还包括:根据下列公式获取归一化边界增强参数N|ED|:
其中,|ED|Wmax为包含计算点的移动窗口内各点第一边界增强参数中的最大值。
6.一种场源边界定位装置,其特征在于,所述装置包括:
第一获取单元,用于在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位;
第二获取单元,用于计算重力场和/或重力梯度张量;所述重力场是所述引力位在笛卡尔坐标系下的一阶导数;所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数;所述重力梯度张量具有多个独立分量;
第三获取单元,用于获取所述重力场的解析信号模,并基于所述解析信号模对所述三维地质体进行边界定位,和/或获取所述重力场的重力梯度张量的方向解析信号的振幅,并基于所述方向解析信号的振幅对所述三维地质体进行边界定位;
其中,还根据重力梯度张量解析信号模对场源体边界进行识别;
其中,还根据下列公式获取第一边界增强参数|ED|:
其中,|Axz|为所述重力梯度张量在x方向解析信号Ax的垂向一阶导数模,|Ayz|为所述重力梯度张量在y方向解析信号Ay的垂向一阶导数模;
其中,还根据所述第一边界增强参数|ED|对异常体进行边界识别;
其中,还根据下列公式获取第二边界增强参数|ED|-z:
所述第二边界增强参数|ED|-z为采用重力梯度张量在z方向解析信号的总水平导数;
根据所述第二边界增强参数|ED|-z对异常体进行边界识别;
其中,还根据下列公式获取第三边界增强参数|ED|_t:
根据所述第三边界增强参数|ED|_t对异常体进行边界识别;
或,
其中,还根据重力梯度张量解析信号的倾斜角方法对三维地质体进行边界定位;
其中,还根据下列公式获取Tilt_Az’:
其中,还根据下列公式获取Tilt_Azz:
其中,还根据下列公式获取Tilt_Azz’:
7.一种计算机可读存储介质,其特征在于,所述计算机可读存储介质存储有一个或者多个程序,所述一个或者多个程序可被一个或者多个处理器执行,以实现权利要求1至5任意一项所述的场源边界定位方法的步骤。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310760291.3A CN116774303B (zh) | 2023-06-26 | 2023-06-26 | 一种场源边界定位方法、装置及计算机可读存储介质 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310760291.3A CN116774303B (zh) | 2023-06-26 | 2023-06-26 | 一种场源边界定位方法、装置及计算机可读存储介质 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN116774303A CN116774303A (zh) | 2023-09-19 |
CN116774303B true CN116774303B (zh) | 2024-05-07 |
Family
ID=88007735
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202310760291.3A Active CN116774303B (zh) | 2023-06-26 | 2023-06-26 | 一种场源边界定位方法、装置及计算机可读存储介质 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN116774303B (zh) |
Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
SU1084822A1 (ru) * | 1982-09-01 | 1984-04-07 | Куйбышевский ордена Трудового Красного Знамени политехнический институт им.В.В.Куйбышева | Устройство дл определени границ аналитического пика |
WO1995005614A1 (en) * | 1993-08-16 | 1995-02-23 | Noranda Inc. | Laplace gravity gradiometer |
CN1561455A (zh) * | 2001-09-27 | 2005-01-05 | 格雷维泰克仪器有限公司 | 用于重力梯度测量的装置 |
CN102466817A (zh) * | 2010-11-11 | 2012-05-23 | 中国石油天然气集团公司 | 归一化导数模法重力异常边界拾取方法 |
CN108508490A (zh) * | 2018-03-07 | 2018-09-07 | 吉林大学 | 一种基于解析信号的磁张量梯度数据均衡边界识别方法 |
CN110414060A (zh) * | 2019-06-28 | 2019-11-05 | 中国地质大学(武汉) | 一种基于四阶谱矩的位场边界识别方法 |
CN113886753A (zh) * | 2021-10-09 | 2022-01-04 | 中国自然资源航空物探遥感中心 | 基于张量特征值的Tilt法航磁边界检测方法、装置 |
CN114325870A (zh) * | 2021-12-14 | 2022-04-12 | 江苏省地质勘查技术院 | 基于三次样条函数计算位场梯度张量的方法及系统 |
CN116299740A (zh) * | 2022-12-08 | 2023-06-23 | 吉林大学 | 一种旋转矩形棱柱体的空间域重力多参量解析正演方法 |
-
2023
- 2023-06-26 CN CN202310760291.3A patent/CN116774303B/zh active Active
Patent Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
SU1084822A1 (ru) * | 1982-09-01 | 1984-04-07 | Куйбышевский ордена Трудового Красного Знамени политехнический институт им.В.В.Куйбышева | Устройство дл определени границ аналитического пика |
WO1995005614A1 (en) * | 1993-08-16 | 1995-02-23 | Noranda Inc. | Laplace gravity gradiometer |
CN1561455A (zh) * | 2001-09-27 | 2005-01-05 | 格雷维泰克仪器有限公司 | 用于重力梯度测量的装置 |
CN102466817A (zh) * | 2010-11-11 | 2012-05-23 | 中国石油天然气集团公司 | 归一化导数模法重力异常边界拾取方法 |
CN108508490A (zh) * | 2018-03-07 | 2018-09-07 | 吉林大学 | 一种基于解析信号的磁张量梯度数据均衡边界识别方法 |
CN110414060A (zh) * | 2019-06-28 | 2019-11-05 | 中国地质大学(武汉) | 一种基于四阶谱矩的位场边界识别方法 |
CN113886753A (zh) * | 2021-10-09 | 2022-01-04 | 中国自然资源航空物探遥感中心 | 基于张量特征值的Tilt法航磁边界检测方法、装置 |
CN114325870A (zh) * | 2021-12-14 | 2022-04-12 | 江苏省地质勘查技术院 | 基于三次样条函数计算位场梯度张量的方法及系统 |
CN116299740A (zh) * | 2022-12-08 | 2023-06-23 | 吉林大学 | 一种旋转矩形棱柱体的空间域重力多参量解析正演方法 |
Non-Patent Citations (15)
Title |
---|
《Three-directional analytic signal analysis and interpretation of magnetic gradient tensor》;Guo Can-wen et al;APPLIED GEOPHYSICS;第17卷(第2期);第285-296页 * |
《不同数据用于欧拉方程的模型计算》;范美宁等;地球物理学进展;第23卷(第4期);第1250-1253页 * |
《中国重力勘探的新进展》;王懋基;地球物理学报;第37卷(第S1期);第353-360页 * |
《全张量重力梯度数据的综合分析与处理解释》;袁园;中国博士学位论文全文数据库 基础科学辑(2015年第08期);第7-9、51-56、62页 * |
《地下目标体重力张量数据的边界识别与增强算法的研究》;张超等;测绘工程;第26卷(第4期);第16-21页 * |
《基于位场梯度张量数据的边界识别方法对比分析研究》;黄松等;东华理工大学学报(自然科学版);20200430;第43卷(第2期);第152-158页 * |
《增强解析信号技术的应用条件》;张季生;物探与化探;第23卷(第4期);第296-300页 * |
《川西藏东地区重磁处理解释与地壳结构研究》;代宏睿;中国优秀硕士学位论文全文数据库 基础科学辑(2022年第04期);第35-37页 * |
《油气模型的重力梯度张量研究》;蒋甫玉等;吉林大学学报(地球科学版);20110331;第41卷(第2期);第545-551页 * |
《重力场Tilt-Depth方法及其高阶推广》;刘鹏飞等;武汉大学学报(信息科学版);20170930;第42卷(第9期);第1236-1242、1277页 * |
《重力梯度张量解析信号的欧拉反褶积》;朱自强等;中南大学学报(自然科学版);第46卷(第1期);第217-222页 * |
《重力梯度张量解析信号的欧拉反褶积》;王灿;中国优秀硕士学位论文全文数据库 基础科学辑(2015年第02期);第1-47页 * |
《重力梯度张量解析信号的正演模拟》;王灿等;中国有色金属学报;第23卷(第09期);第2479-2497页 * |
《高精度重磁资料在松辽盆地古龙断陷火山岩气藏勘探中的应用》;王玉华等;石油地球物理勘探;第43卷(第1期);第107-112页 * |
袁园.《全张量重力梯度数据的综合分析与处理解释》.中国博士学位论文全文数据库 基础科学辑.2015,(2015年第08期),第7-9、51-56、62页. * |
Also Published As
Publication number | Publication date |
---|---|
CN116774303A (zh) | 2023-09-19 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Fedi et al. | Understanding imaging methods for potential field data | |
Ren et al. | Gravity gradient tensor of arbitrary 3D polyhedral bodies with up to third-order polynomial horizontal and vertical mass contrasts | |
Pilkington | Evaluating the utility of gravity gradient tensor components | |
Oruç | Depth estimation of simple causative sources from gravity gradient tensor invariants and vertical component | |
Ren et al. | Closed-form formula of magnetic gradient tensor for a homogeneous polyhedral magnetic target: A tetrahedral grid example | |
Martinez et al. | Denoising of gravity gradient data using an equivalent source technique | |
CN110333536A (zh) | 一种测距线性定位算法 | |
Li et al. | Using an equivalent source with positivity for low-latitude reduction to the pole without striation | |
Jiang et al. | A versatile solution for the gravity anomaly of 3D prism-meshed bodies with depth-dependent density contrast | |
Zhou | 3D vector gravity potential and line integrals for the gravity anomaly of a rectangular prism with 3D variable density contrast | |
Rim et al. | Gravity gradient tensor due to a cylinder | |
Ren et al. | Recursive analytical formulae of gravitational fields and gradient tensors for polyhedral bodies with polynomial density contrasts of arbitrary non-negative integer orders | |
Lindsay et al. | Inversion and geodiversity: searching model space for the answers | |
Ouyang et al. | Iterative magnetic forward modeling for high susceptibility based on integral equation and Gauss-fast Fourier transform | |
Dai et al. | The forward modeling of 3D gravity and magnetic potential fields in space-wavenumber domains based on an integral method | |
Abdelrahman et al. | Magnetic interpretation using a least-squares, depth-shape curves method | |
Fregoso et al. | Initializing cross-gradients joint inversion of gravity and magnetic data with a Bayesian surrogate gravity model | |
Dai et al. | Three-dimensional numerical modeling of gravity anomalies based on Poisson equation in space-wavenumber mixed domain | |
Huang et al. | Three-dimensional arbitrarily anisotropic modeling for time-domain airborne electromagnetic surveys | |
CN111123380B (zh) | 基于重磁梯度数据张量不变量的目标深度估计方法及系统 | |
CN112596113A (zh) | 一种基于重力不同阶梯度特征值交点的场源位置识别方法 | |
CN116774303B (zh) | 一种场源边界定位方法、装置及计算机可读存储介质 | |
Godio et al. | Integrated data processing for archeological magnetic surveys | |
Zhou et al. | Depth estimation of potential field by using a new downward continuation based on the continued fraction in space domain | |
Rim et al. | Borehole vector gravity: A numerical study |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |