一种场源边界定位方法、装置及计算机可读存储介质与流程

1.本发明涉及地球物理领域，尤其涉及一种场源边界定位方法及装置及计算机可读存储介质。


背景技术：

2.重力梯度张量测量在国外已开展了十余年,并取得了一定的成果,然而这项测量技术及相关数据解释技术在我国的研究仍然较少。国内在全张量重力梯度仪的研发过程中还处于初级阶段。目前，有多家单位进行相关研究，主要是进行相关重力梯度仪的理论设计、误差分析和模型仿真，都处于实验室研发阶段，尚不能达到勘探要求。重力梯度数据的解释工作，可以进行单独分量的解释，也可以联合多个分量进行解释。在解释中，选择哪些分量或者如何结合不同的分量进行解释，便成为张量梯度数据解释工作所需要解决的重要问题之一。基于系统的发展，梯度测量数据的评估、分析也逐渐形成一套完善的理论体系。梯度数据和重力异常数据对比评估，梯度数据噪音的去除，以及梯度数据的处理方法等问题均成为重要的研究课题。与此同时，解释方法技术也被大量建立和发展，且和重力异常数据的解释方法相依托，形成了目前的重力异常及其梯度数据的解释现状。重力数据的解释工作，通常需要计算出重力场源体的三种基本信息：水平位置，深度范围，物性差异。水平位置的确定，常用到的解释方法即边界识别和增强方法，这种方法可以有效地突显出地质体边缘及构造的水平位置。深度计算通常需要根据异常的特征点及地质体形体来进行估算。物性参数的计算，有常用的物性反演计算，或者可采用快速成像算法快速给出与物性相关的量。重力异常边界识别技术利用了地质体在断裂位置、接触带附近、围岩周边位置的密度变化率较大，从而引起重力异常有明显变化的特点，通过异常的二次特征对变化率较大的位置进行判读，从而达到地质体边界识别(水平位置或在平面上的投影位置)的目的。
3.但目前基于张量数据进行边界识别的方法存在一些问题：目标体的边界识别效果不佳，对三维地质体的识别和定位效果并不好。因此，发明一种新的场源边界定位方法技术势在必行。


技术实现要素：

4.本发明提供了一种场源边界定位方法及装置及计算机可读存储介质，用于解决现有技术中基于张量数据进行边界识别对目标体的边界识别效果不佳，对三维地质体的识别和定位效果并不好的问题。
5.本发明中的一种场源边界定位方法，所述方法包括：
6.在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位；
7.计算重力场和/或重力梯度张量；所述重力场是所述引力位在笛卡尔坐标系下的一阶导数；所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数；所述重力梯度张量具有多个独立分量；
8.获取所述重力场的解析信号模，并基于所述解析信号模对所述三维地质体进行边
界定位，和/或获取所述重力场的重力梯度张量的方向解析信号的振幅，并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。
9.可选的，所述三维地质体为直立六面体模型。
10.可选的，所述方法还包括：计算所述重力梯度张量的一阶混合导数以对所述三维地质体进行边界定位。
11.可选的，所述方法还包括：根据下列公式获取第一边界增强参数|ed|：
[0012][0013]
其中，|a
xz
|为所述重力梯度张量在x方向解析信号a
x
的垂向一阶导数模，|a
yz
|为所述重力梯度张量在y方向解析信号ay的一阶导数模；
[0014]
根据所述第一边界增强参数|ed|对异常体进行边界识别。
[0015]
可选的，所述方法还包括：根据下列公式获取第二边界增强参数|ed|
_z
：
[0016][0017]
所述第二边界增强参数|ed|
_z
为采用重力梯度张量在z方向解析信号的总水平导数；
[0018]
根据所述第二边界增强参数|ed|
_z
对异常体进行边界识别。
[0019]
可选的，所述方法还包括：根据下列公式获取第三边界增强参数|ed|_
t
：
[0020][0021]
根据所述第三边界增强参数|ed|
_t
对异常体进行边界识别。
[0022]
可选的，所述第一边界增强参数、第二边界增强参数以及第三边界增强参数的极大值位置均对应所述异常体的边界。
[0023]
可选的，所述方法还包括：根据下列公式获取归一化边界增强参数n|ed|：
[0024][0025]
其中，|ed|
wmax
为包含计算点的移动窗口内各点第一边界增强参数中的最大值。
[0026]
本发明中的一种场源边界定位装置，所述装置包括：
[0027]
第一获取单元，用于在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位；
[0028]
第二获取单元，用于计算重力场和/或重力梯度张量；所述重力场是所述引力位在笛卡尔坐标系下的一阶导数；所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数；所述重力梯度张量具有多个独立分量；
[0029]
第三获取单元，用于获取所述重力场的解析信号模，并基于所述解析信号模对所述三维地质体进行边界定位，和/或获取所述重力场的重力梯度张量的方向解析信号的振幅，并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。
[0030]
本发明中的一种计算机可读存储介质，所述计算机可读存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现如上任意一项所述的场源边界定位方法的步骤。
附图说明
[0031]
图1是本发明实施例中一种场源边界定位方法的流程示意图；
[0032]
图2是本发明实施例中笛卡尔坐标系下的长方体模型及地下空间剖分；
[0033]
图3是本发明实施例中直立六面体示意图；
[0034]
图4是本发明实施例中模型一的引力位一阶导数示意图；
[0035]
图5是本发明实施例中模型一的重力梯度张量异常示意图；
[0036]
图6是本发明实施例中模型一的部分重力梯度张量混合一阶导数示意图；
[0037]
图7是本发明实施例中模型一的重力梯度张量解析信号模示意图；
[0038]
图8是本发明实施例中模型一的重力梯度张量解析信号模的一阶导数示意图；
[0039]
图9是本发明实施例中模型一|ed|、|ed|z、|ed|
t
与thd、asmz、tilt_asm对比示意图；
[0040]
图10是本发明实施例中模型二|ed|、|ed|z、|ed|
t
与thd、asmz对比示意图；
[0041]
图11是本发明实施例中模型三|ed|、|ed|z、|ed|
t
与thd、asmz对比示意图；
[0042]
图12是本发明实施例中模型四|ed|、|ed|-z
、|ed|-t
边界增强结果示意图；
[0043]
图13是本发明实施例中模型四n|ed|、n|ed|-z
、n|ed|-t
边界增强结果对比示意图；
[0044]
图14是本发明实施例中模型五tilt、tilt_az和tilt_az′
反边界增强结果示意图；
[0045]
图15是本发明实施例中模型五tilt_a
zz
与tilt_a
zz
‘
边界增强结果对比示意图；
[0046]
图16是本发明实施例中模型六tilt_a
zz
与tilt_a
zz
‘
边界增强结果对比示意图；
[0047]
图17是本发明实施例中一种场源边界定位的结构示意图。
具体实施方式
[0048]
下面结合附图和实施例对本发明作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释本发明，而非对本发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本发明相关的部分而非全部结构。
[0049]
应理解，在本文的各种实施例中，上述各过程的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本文实施例的实施过程构成任何限定。
[0050]
利用重磁位场资料进行地质体的识别推断，如岩体范围、断裂的位置等场源边界的识别，是进一步深化研究的基础。边界识别或边界定位常指利用一系列的点或线段明确标定如断裂构造线或地质体的水平投影线位置的表达方式。但由于受到场源埋深、异常叠加以及噪声等多种因素的影响，直接进行场源的边界识别往往存在很大的误差，需要开展边界特征增强，以便于进行场源边界识别工作的开展。
[0051]
重力梯度张量数据相对于传统的重力测量数据有着更高的频率信息，能更加准确、细致地研究地球内部构造、矿产资源分布等信息。本发明具体实施例是基于重力场和/或重力梯度张量重解析信号及其导数，并在此基础上，经过再求导、归一化等方法实现场源体边界增强，以提高边界定位的准确度。
[0052]
本发明具体实施例提供一种场源边界定位方法，如图1所示，所述方法包括：
[0053]
步骤100，在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位。具体的，在笛卡尔坐标系下，设地面某一点o为坐标原点，根据万有引力定律，一个剩余
质量为m，剩余密度为ρ的地质体在其外部空间任意一点产生的引力位为：
[0054][0055]
其中，g＝6.67
×
10-11
m3/(kg
·
s)2为万有引力常量，为任意一点p(x,y,z)到场源点q(ξ,η,ζ)的矢量半径，r＝[(ξ-x)2+(η-y)2+(ζ-z)2]
1/2
，如图2所示。
[0056]
步骤200，计算重力场和/或重力梯度张量；所述重力场是所述引力位在笛卡尔坐标系下的一阶导数；所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数；所述重力梯度张量具有多个独立分量。具体的实施例中，步骤200可以包括不同的实施方式。在第一种实施方式中，只计算重力场。其中，重力场是引力位在笛卡尔坐标系x，y，z三个方向上的一阶导数，即重力加速度矢量可以得到：
[0057][0058][0059][0060][0061]
在另一具体的实施方式中，还需计算重力梯度张量，其中：重力梯度张量为地质体引力位的二阶导数，即为重力加速度矢量在笛卡尔坐标系x，y，z三个方向的一阶导数，也即重力梯度：
[0062][0063]
写成矩阵形式，即
[0064][0065]
其中：v为地质体的引力位，为引力位在x方向的二阶导数，即g
xx
；
[0066]
为引力位在x、y方向的二阶导数，即g
xy
，其他参数的定义可参照如上。
[0067]
步骤300，获取所述重力场的解析信号模，并基于所述解析信号模对所述三维地质体进行边界定位，和/或获取所述重力场的重力梯度张量的方向解析信号的振幅，并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。
[0068]
在一具体的实施方式中，可以直接通过重力场进行边界定位，但定位边界较为模糊。在另一具体的实施方式中，还可以通过解析信号模对所述三维地质体进行边界定位。具体的，由于重力场是无源场，因此其旋度和散度均为零，即：
[0069][0070][0071]
所以重力梯度张量为对称矩阵，并且迹为零：
[0072]gxx
+g
xy
+g
xz
＝0
[0073]gxy
＝g
yx
,g
xz
＝g
zx
,g
yz
＝g
zy
[0074]
因此，重力梯度张量的9个分量中，只有5个独立分量。
[0075]
二维地质体位场φ(x)的解析信号a(x)可以写成复函数：
[0076][0077]
其中和是希尔伯特变换对。二维解析信号的振幅表示为：
[0078][0079]
将二维解析信号归纳到三维，并且任何位场的希尔伯特变换满足柯西—黎曼关系。将水平面观测的位场φ(x,y)的解析信号的定义扩展到三维，如下：
[0080][0081]
则a(x,y)的振幅表示为：
[0082][0083]
在三维情况下，重力梯度张量(ggt)第三列元素是第一和第二列元素的希尔伯特变换对。然后，可以定义在x、y和z方向的方向解析信号，其矩阵形式表示为：
[0084][0085]
因此，ggt的方向解析信号的振幅表示为：
[0086][0087][0088][0089]
对于多个相互影响的异常体，利用解析信号振幅的导数可以实现异常体分离。x、y、z方向上的方向性解析信号的导数可以表示为：
[0090][0091][0092][0093]
其中，α是x、y、z。
[0094]
本发明具体实施例所述的方法从重力场和重力梯度理论出发，对重力场和重力梯度张量数据的解释技术进行探索性研究，基于重力场和重力梯度张量解析信号及其方向导数，提出了新的边界定位方法在目标体的边界识别上效果更好。对今后张量技术实际应用过程中的数据处理与解释提供一定的理论基础。
[0095]
本发明具体实施例所述的一种场源边界定位方法，较佳的，所述三维地质体为直立六面体模型。具体的，直立六面体模型是一种非常典型和常见的三度体地质模型，在地球物理的实际工作中很常见的侵入岩、岩床、水平板、垂直板、岩脉、岩墙以及基底变质岩系都可以用直立六面体进行理论模拟。另外，通过将地下介质剖分成足够小的直立六面体，通过多个小直立六面体的组合，可以模拟三维复杂模型，实现三维复杂模型的正反演。因此对直立六面体模型进行模拟试算具有很强的实际意义。
[0096]
本发明具体实施例所述的一种场源边界定位方法，较佳的，所述方法还包括：计算所述重力梯度张量的一阶混合导数以对所述三维地质体进行边界定位。
[0097]
具体的实施例中，为了识别三维地质体中存在的异常体，建立如图3所示的笛卡尔坐标系，设直立六面体两角点坐标分别为(ξ1,η1,ζ1)和(ξ2,η2,ζ2)，剩余密度为σ，在地面上任意观测点p(x,y,z)的重力异常、重力异常一阶导数及二阶导数的计算如下。
[0098]
直立六面体引力位的一阶导数计算公式，一共有3个分量，如下：
[0099][0100][0101][0102]
通过分析可以看出，当上面三个式子中反正切函数的分子分母都为0时，这个分式无意义，此时会出现“奇点”。因此重新推导g
x
和gy，得到了不含解析奇点的计算公式：
[0103][0104][0105]
建立单个异常体的模型一，异常体长宽高为500
×
500
×
500m3，剩余密度100kg/m3，顶面埋深20m，中心点坐标(0，0，270)。模型一的引力位一阶导数如图4所示。
[0106]
直立六面体引力位的二阶导数计算公式，即直立六面体重力梯度张量，一共有9个分量，分别如下：
[0107][0108][0109][0110][0111][0112][0113]gzx
＝g
xz
,g
zy
＝g
yz
,g
yx
＝g
xy
[0114]
式中
[0115]
同样该公式中g
xx
,g
xy
,g
xz
也存在解析奇点。因此推导出无解析奇点离散计算公式：
[0116][0117][0118][0119]
[0120][0121][0122]
图5是模型一的重力梯度张量异常，可以看出g
xx
能有效地增强地质体的x方向边界，g
yy
能增强异常的y方向边界，g
xy
显示了立方体四个角点异常的峰值，能描述地质体四个方位角点信息，g
zx
的极大值可以用来识别x方向边界，g
zy
的极大值能识别地质体的y方向边界，不同张量异常可以反映地质体不同边界信息。
[0123]
通过对引力位的二阶导数求导，即对重力梯度张量的混合一阶导数，可以得到引力位三阶混合导数，一共17个分量，分别为：
[0124][0125][0126][0127][0128][0129][0130][0131][0132][0133]
[0134][0135][0136][0137][0138][0139][0140]gxyz
＝g
yxz
＝g
xzy
＝g
zxy
[0141]
图6所示为模型一的部分重力梯度张量混合一阶导数，可以看出在各个方向上其具有更高的分辨率。
[0142]
本发明一具体的实施例中，直立六面体的重力梯度张量在x、y、z方向上的方向解析信号可以直接利用重力梯度张量的6个分量计算获得。
[0143]
图7为模型一的x、y、z方向解析信号模。可以看出，解析信号模的优点是不会出现假值，只有在边界位置出现高值；但是边界增强的横向分辨率不足，边界位置的过渡带较宽。只识别单个方向的边界，比如只识别南北边界，东西边界幅值非常弱或几乎没有增强的效果。因此本发明具体实施例所述的一种场源边界定位方法中，可以利用引力位的二阶导数、三阶导数、ggt解析信号计算获得，一共9个分量，分别为：
[0144][0145][0146][0147][0148][0149][0150]
[0151][0152][0153]
图8为模型一的解析信号一阶导数。相比解析信号模，垂向导数提高了边界极值的横向分辨率，边界过渡带缩小了很多，非边界区域也没有假值，但与解析信号模类似，只识别单个方向的边界，另一个方向的边界没有识别或者识别结果较为模糊。
[0154]
基于图8分析重力梯度张量解析信号一阶导数的特征，可以看出axz反应的x方向边界更清晰，ayz反应的y方向边界更清晰，因此本发明具体实施例所述的一种场源边界定位方法，较佳的，所述方法还包括：根据下列公式获取第一边界增强参数|ed|：
[0155][0156]
其中，|a
xz
|为所述重力梯度张量在x方向解析信号a
x
的垂向一阶导数模，|a
yz
|为所述重力梯度张量在y方向解析信号ay的一阶导数模；
[0157]
根据所述第一边界增强参数|ed|对异常体进行边界识别。
[0158]
具体的，第一边界增强参数|ed|的定义类似于传统解析信号的垂向导数asz，其实质是重力梯度张量x方向解析信号a
x
的垂向一阶导数模|a
xz
|与y方向解析信号ay的一阶导数模|a
yz
|的组合，通过第一边界增强参数|ed|可以增加对异常体的边界识别。
[0159]
较佳的实施例中，除了a
xz
、a
yz
具有较好的边界增强效果外，a
zx
、a
zy
也具有较好的边界增强效果，据此又定义了第二边界增强参数|ed|
_z
。较佳的，所述方法还包括：根据下列公式获取第二边界增强参数|ed|
_z
：
[0160][0161]
所述第二边界增强参数|ed|
_z
为采用重力梯度张量在z方向解析信号的总水平导数；
[0162]
根据所述第二边界增强参数|ed|
_z
对异常体进行边界识别。
[0163]
具体的，|ed|
_z
的定义则类似于传统的总水平导数，即采用重力梯度张量z方向的解析信号计算总水平导数。
[0164]
本发明具体实施例所述的一种场源边界定位方法，联合|ed|、|ed|-z
两个参量又定义了第三边界增强参数|ed|-t
。较佳的，所述方法还包括：根据下列公式获取第三边界增强参数|ed|
_t
：
[0165][0166]
根据所述第三边界增强参数|ed|
_t
对异常体进行边界识别。
[0167]
本发明具体实施例所述的一种场源边界定位方法，较佳的，所述第一边界增强参数、第二边界增强参数以及第三边界增强参数的极大值位置均对应所述异常体的边界。
[0168]
上述|ed|、|ed|-z
、|ed|-t
均为边界增强参数，其极大值位置对应异常体的边界。采用四方位识别极大值的方法可以自动圈定异常体的边界，实现边界识别。
[0169]
对上述三个边界识别参数与传统效果较好的thd、asmz、tilt_asm方法进行对比，
如图9所示。需要注意的是，理论上asmz与a
zz
相同。另外，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效果相对其它五种方法较差。
[0170]
采用模型二再次进行比较，如图10所示。模型二与模型一相似，唯一的变化为深度变为350～850m。从图10可以看出，|ed|、asmz已经不能识别边界，|ed|-t
边界范围变小，|ed|-z
、|ed|边界结果相似，范围更宽。
[0171]
模型三为组合模型，该模型包含两个直立六面体，浅部异常体200
×
200
×
200m3，顶部埋深20m，剩余密度差100kg/m3；深部异常体750
×
750
×
750m3，顶部埋深220m，剩余密度差1300kg/m3。从图11可以看出：
[0172]
(1)thd已经不能识别浅部的异常体，|ed|-t
、|ed|-z
、|ed|、asmz均可以反映浅部异常体，其中|ed|-z
的异常衬度最大；
[0173]
(2)|ed|-z
的异常形态反映深部异常体最明显，其次为|ed|-t
，再次为|ed|，asmz最差，thd虽然也能反映深部异常体，但是异常范围过宽；
[0174]
(3)|ed|、asm识别出的边界比真实模型范围稍小，|ed|-z
范围稍大，|ed|-t
识别的边界范围与真实模型更加接近。
[0175]
通过上述三个模型的对比，得到如下结论：
[0176]
(1)采用ggt解析信号一阶导数组合定义的|ed|、|ed|-z
、|ed|-t
总体效果好于thd、asmz、tilt_asm；
[0177]
(2)|ed|-z
的定义和边界增强效果类似于hdt，|ed|-t
的定义和边界增强效果类似于asmz；
[0178]
(3)|ed|浅部异常体的边界反映最好，但是深部异常体的边界反映较差；
[0179]
(4)|ed|
t
、|ed|z浅部异常体的边界反映很好，也可以反映深部异常体的边界；
[0180]
|ed|、|ed|-z
、|ed|-t
存在的主要问题是对弱异常不敏感，对深部地质体的异常信息缺少足够的反映。
[0181]
本发明具体实施例所述的一种场源边界定位方法，较佳的，所述方法还包括：根据下列公式获取归一化边界增强参数n|ed|：
[0182][0183]
其中，|ed|
wmax
为包含计算点的移动窗口内各点第一边界增强参数中的最大值。
[0184]
具体的，在模型四中包含三个直立六面体，尺寸均为40
×
160m，剩余密度差均为300kg/m3，由西向东顶部埋深依次为10m、25m、40m。从图12可以看出|ed|、|ed|-z
、|ed|-t
与thd对浅部异常体都可以识别，但是深部两个异常体则很难识别。
[0185]
因此，本发明具体实施例借鉴归一化总水平导数法，针对|ed|、|ed|-z
、|ed|-t
对大埋深场源的边界增强不足的特点，利用移动窗口内的最大值归一化计算点的值，并利用极大值标定场源边界，以期增强弱异常，来突出深部异常。归一化后的归一化边界增强参数ned表达式为：
[0186][0187]
其中：|ed|表示计算点的第一边界增强参数，|ed|
wmax
为包含计算点的移动窗口内各点边界增强参数中的最大值，较佳的实施例中，|ed|也可以用|ed|-z
、|ed|-t
进行替换。
[0188]
图13为模型四采用归一化的方法计算得到的n|ed|、n|ed|-z
、n|ed|-t
边界增强结果对比，计算窗口为80
×
80m。从图13中可以看出：
[0189]
(1)采用归一化的方法确实能增强深部异常，但是深部异常边界的极值范围较宽；
[0190]
(2)n|ed|的极值并不反映边界而是反映中心；
[0191]
(3)n|ed|-z
、n|ed|-t
的极值反映边界，n|ed|
_z
的异常衬度更大。
[0192]
较佳的实施例中，归一化类方法的边界增强效果与窗口的选取有关。
[0193]
本发明具体实施例中，还提供一种解析信号倾斜角tilt_asm的方法采用gz的解析信号asm替换gz进行边界增强，具体的计算方法为：
[0194][0195]
传统tilt_asm是针对gz参数开展多次频域求导。分析上式，as就是重力梯度张量数据的解析信号az。因此，可写为：
[0196][0197]
上式需要计算az的x和y方向一阶导数，采用相邻点的有限差分代替求导，差分间距为两倍点距。但差分间距小于1/10～1/100倍点距时，有限差分的误差才可以忽略。本发明实施例中通过正演可以直接计算|a
zx
|2与|a
zy
|2，其精度比两倍点距有限差分计算结果高。
[0198]
由于重力梯度张量数据提供了更加丰富的信息(a
x
与ay)，理论上a
xx
与a
yy
分别对x与y方向的边界更加敏感。|a
xx
|2与|a
yy
|2的0值位置对应异常体边界。采用倾斜角方法，本发明实施例提出了重力梯度张量解析信号的tilt方法，公式为：
[0199][0200]
模型五包含三个埋深不同、大小完全相同的直立六面体，埋藏深度分别为4km(棱柱体1)、1km(棱柱体2)和2.5km(棱柱体3)，所有直立六面体的厚度都为1km，剩余密度为200kg/m3。
[0201]
图14中的(a)和(b)为模型五tilt、tilt_az边界增强结果与(c)的文献结果对比，从图中可以看出，传统tilt方法异常体内为正值，异常体外为负值，边界为0值；tilt_az则可以根据极大值反映边界位置，且边界清晰；tilt_az′
反映的边界更加规整，与模型的真实边界吻合更好。但是tilt_az和tilt_az′
在异常体内部和异常体中间位置存在假异常。
[0202]
解析信号的导数相对于解析信号具有更高的分辨率。利用解析信号的导数，对上
两式做进一步的改进为：
[0203][0204][0205]
另外，根据获得的异常形态及特征值分布情况，可以对tilt_a
zz
、tilt_a
zz
′
进行大于0阈值的峰值提取。以tilt_a
zz
为例，有：
[0206][0207]
图15的(a)和(b)为模型五tilt_a
zz
与tilt_a
zz
′
边界增强结果对比，从图中可以看出，tilt_a
zz
与tilt_a
zz
′
边界增强后的异常更加集中，与模型的真实边界吻合更好。
[0208]
为了验证该方法对包含正、负异常数据的适用性，将模型五中东北部直立六面体的密度差修改为-200kg/m3，其它参数不变，构建模型六。图16的(a)和(b)为模型六tilt_a
zz
与tilt_a
zz
′
边界增强结果对比，从结果可以看出，tilt_a
zz
与tilt_a
zz
′
对负异常依旧适用。
[0209]
本发明具体实施例还提供一种场源边界定位装置，如图17所示，所述装置包括：
[0210]
第一获取单元1701，用于在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位；
[0211]
第二获取单元1702，用于计算重力场和/或重力梯度张量；所述重力场是所述引力位在笛卡尔坐标系下的一阶导数；所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数；所述重力梯度张量具有多个独立分量；
[0212]
第三获取单元1703，用于获取所述重力场的解析信号模，并基于所述解析信号模对所述三维地质体进行边界定位，和/或获取所述重力场的重力梯度张量的方向解析信号的振幅，并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。
[0213]
本发明实施例还提供一种计算机可读存储介质，所述计算机可读存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现如上具体实施例中的任意一项所述的场源边界定位方法的步骤。
[0214]
本发明实施例所述的方法通过对重力梯度张量解析信号与欧拉反褶积相结合，进行了地质体定位方法的研究，通过数值模拟进行试算，发现对于单个直立六面体与规则组合模型的边界识别和定位效果较好，对三维地质体具有较高的可行性，以期对今后张量技术实际应用过程中的数据处理与解释提供一定的理论基础。
[0215]
本领域内的技术人员应明白，本技术的实施例可提供为方法、系统、或计算机程序产品。因此，本技术可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本技术可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
[0216]
本技术是参照根据本技术实施例的方法、设备(系统)、和计算机程序产品的流程
图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0217]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0218]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0219]
前述对本发明的具体示例性实施方案的描述是为了说明和例证的目的。这些描述并非想将本发明限定为所公开的精确形式，并且很显然，根据上述教导，可以进行很多改变和变化。对示例性实施例进行选择和描述的目的在于解释本发明的特定原理及其实际应用，从而使得本领域的技术人员能够实现并利用本发明的各种不同的示例性实施方案以及各种不同的选择和改变。本发明的范围意在由权利要求书及其等同形式所限定。

技术特征：
1.一种场源边界定位方法，其特征在于，所述方法包括：在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位；计算重力场和/或重力梯度张量；所述重力场是所述引力位在笛卡尔坐标系下的一阶导数；所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数；所述重力梯度张量具有多个独立分量；获取所述重力场的解析信号模，并基于所述解析信号模对所述三维地质体进行边界定位，和/或获取所述重力场的重力梯度张量的方向解析信号的振幅，并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。2.根据权利要求1所述的一种场源边界定位方法，其特征在于，所述三维地质体为直立六面体模型。3.根据权利要求1所述的一种场源边界定位方法，其特征在于，所述方法还包括：计算所述重力梯度张量的一阶混合导数以对所述三维地质体进行边界定位。4.根据权利要求1所述的一种场源边界定位方法，其特征在于，所述方法还包括：根据下列公式获取第一边界增强参数|ed|：其中，|a
xz
|为所述重力梯度张量在x方向解析信号a
x
的垂向一阶导数模，|a
yz
|为所述重力梯度张量在y方向解析信号a
y
的一阶导数模；根据所述第一边界增强参数|ed|对异常体进行边界识别。5.根据权利要求4所述的一种场源边界定位方法，其特征在于，所述方法还包括：根据下列公式获取第二边界增强参数|ed|
_z
：所述第二边界增强参数|ed|
_z
为采用重力梯度张量在z方向解析信号的总水平导数；根据所述第二边界增强参数|ed|
_z
对异常体进行边界识别。6.根据权利要求5所述的一种场源边界定位方法，其特征在于，所述方法还包括：根据下列公式获取第三边界增强参数|ed|
_t
：根据所述第三边界增强参数|ed|
_t
对异常体进行边界识别。7.根据权利要求6所述的一种场源边界定位方法，其特征在于，所述第一边界增强参数、第二边界增强参数以及第三边界增强参数的极大值位置均对应所述异常体的边界。8.根据权利要求6所述的一种场源边界定位方法，其特征在于，所述方法还包括：根据下列公式获取归一化边界增强参数n|ed|：其中，|ed|
wmax
为包含计算点的移动窗口内各点第一边界增强参数中的最大值。9.一种场源边界定位装置，其特征在于，所述装置包括：第一获取单元，用于在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的
引力位；第二获取单元，用于计算重力场和/或重力梯度张量；所述重力场是所述引力位在笛卡尔坐标系下的一阶导数；所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数；所述重力梯度张量具有多个独立分量；第三获取单元，用于获取所述重力场的解析信号模，并基于所述解析信号模对所述三维地质体进行边界定位，和/或获取所述重力场的重力梯度张量的方向解析信号的振幅，并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现权利要求1至8任意一项所述的场源边界定位方法的步骤。

技术总结
本发明公开一种场源边界定位方法、装置及计算机可读存储介质，所述方法包括：在笛卡尔坐标系下获取三维地质体在其外部空间任意一点产生的引力位；计算重力场和/或重力梯度张量；所述重力场是所述引力位在笛卡尔坐标系下的一阶导数；所述重力梯度张量是所述引力位在笛卡尔坐标系下的二阶导数；所述重力梯度张量具有多个独立分量；获取所述重力场的解析信号模，并基于所述解析信号模对所述三维地质体进行边界定位，和/或获取所述重力场的重力梯度张量的方向解析信号的振幅，并基于所述方向解析信号的振幅对所述三维地质体进行边界定位。通过本发明的方法对目标体的边界识别效果增强，且对三维地质体的识别和定位效果较好。且对三维地质体的识别和定位效果较好。且对三维地质体的识别和定位效果较好。


技术研发人员：王萌 王珺璐 张广兵 张加洪 李瑜
受保护的技术使用者：中国自然资源航空物探遥感中心
技术研发日：2023.06.26
技术公布日：2023/9/20