CN113138426B - 一种基于多变量核密度欧拉解概率密度成像方法 - Google Patents

一种基于多变量核密度欧拉解概率密度成像方法 Download PDF

Info

Publication number
CN113138426B
CN113138426B CN202110390937.4A CN202110390937A CN113138426B CN 113138426 B CN113138426 B CN 113138426B CN 202110390937 A CN202110390937 A CN 202110390937A CN 113138426 B CN113138426 B CN 113138426B
Authority
CN
China
Prior art keywords
euler
data
grid
dimensional
probability density
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
Application number
CN202110390937.4A
Other languages
English (en)
Other versions
CN113138426A (zh
Inventor
朱自强
杨博
鲁光银
曹书锦
毛雅静
张亮
张东兴
郭晓旺
陈新跃
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202110390937.4A priority Critical patent/CN113138426B/zh
Publication of CN113138426A publication Critical patent/CN113138426A/zh
Application granted granted Critical
Publication of CN113138426B publication Critical patent/CN113138426B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00Measuring gravitational fields or waves; Gravimetric prospecting or detecting
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/15Correlation function computation including computation of convolution operations
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/18Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Operations Research (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computing Systems (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Probability & Statistics with Applications (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Complex Calculations (AREA)

Abstract

本发明提供一种基于多变量核密度欧拉解概率密度成像方法,首先通过测量或位场变换获得重力矢量及重力梯度张量数据,构建三维重力/重力梯度张量欧拉反褶积方程,逐一移动滑动窗口,获得一系列含异常源空间位置[x,y,z,N]欧拉解集,针对常规欧拉反褶积仅通过欧拉反褶积解的构造指数的“色谱”难于评价、分离异常,而不综合考虑欧拉解的空间分布特性及构造指数的相关性。采用三维核密度估计对欧拉反褶积解集综合分析,从而实现了欧拉反褶积对相邻异常有效分离、快速定位。

Description

一种基于多变量核密度欧拉解概率密度成像方法
技术领域
本发明涉及重力勘探技术领域,特别是涉及一种基于多变量核密度欧拉解概率密度成像方法。
背景技术
欧拉反褶积以Euler齐次方程为理论基础,利用总场及其梯度等位场数据,只需事先确定与场源性质有关的构造指数或枚举构造指数,便可快速有效地圈出异常源的基本轮廓,尤其适合于大面积位场数据的分析和解释。然而,使用欧拉反褶积对异常进行推断时,传统欧拉解图示方法一般使用色阶或不同大小圆来标定深度或构造指数,即色谱。但其难于标示异常源间及各单一欧拉解间的差异及相关性,尤其当欧拉解发散时,如对于在一处大量聚集(如欧拉解聚集程度大的地方)和在它处成零星分布(如多异常源间的伪迹)的两类欧拉解集,其色谱的空间分布上几乎没有差别。因而,如何在众多的欧拉解中提取优解、剔除谬解、分离欧拉解簇继而标定异常源一直是困扰人们的一个难题。
发明内容
鉴于以上所述现有技术的缺点,本发明的目的在于提供一种基于多变量核密度欧拉解概率密度成像方法,用于解决现有技术中存在的问题。
为实现上述目的及其他相关目的,本发明提供一种基于多变量核密度欧拉解概率密度成像方法,包括以下步骤:
确定待测区域范围;并测量重力矢量数据及重力梯度张量数据,或测量重力数据,通过离散余弦变换或傅里叶变换将重力数据换算到重力矢量数据及重力梯度张量数据;
根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程;
确定滑动窗口大小,并利用滑动窗口内的数据构建线性方程组;
利用奇异值分解所述线性方程组,获取所述线性方程组的欧拉解,并输出对应的欧拉解集;其中,所述欧拉解包括异常源空间位置、异常源类型和奇异值分解误差;
将欧拉解解集拆分成不同维度的组合,并根据所述欧拉解解集组合估算估计网格最大或最小带宽,及确定所述估计网格的大小;
将欧拉解数据作为独立同分布采样投影至所述估计网格,根据所述估计网络计算网格计数,并基于所述网格计数与核函数的卷积确定不同维度数据的异常源空间位置及类型。
可选地,所述根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程,包括:
Figure BDA0003016673840000021
Figure BDA0003016673840000022
Figure BDA0003016673840000023
式中,Bα为异常背景场,gα为重力梯度,gαβ为重力张量分量,{α,β∈x,y,z}; (x0,y0,z0)为测点位置;(x,y,z)为待求解位场异常源位置;构造指数N是位场异常随场源深度衰减变化陡缓的量度,即异常源类型。
可选地,所述确定滑动窗口大小,并利用滑动窗口内的数据构建线性方程组,以及利用奇异值分解所述线性方程组,获取所述线性方程组的欧拉解,并输出对应的欧拉解集,包括:
确定滑动窗口大小为Wx×Wy;
将滑窗中的Wx×Wy个测点数据代入(1)~(3)构建线性方程组,如下:
Am=b (4);
式中,算子A为
Figure BDA0003016673840000024
b为
Figure BDA0003016673840000025
m为欧拉反褶积解集;其中,第i个欧拉反褶积解可以表示为mi=[xi,yi,zi,Ni];
利用奇异值分解求取异常源空间位置[x,y,z]及构造指数N。
逐一移动滑动窗口,构建线性方程组,并获得相应滑动窗口的欧拉解mi,直至滑动窗口覆盖整个测区,构建欧拉解集{mi}。
可选地,将欧拉解解集拆分成不同维度的组合,并根据所述不同欧拉解解集组合估算估计网格最大或最小带宽,及确定所述估计网格的大小,包括:
引入多维核密度估计,设独立同分布采样X=(X1,…,Xi,…,Xn)T为取值于d维实数空间的独立同分布随机变量;
对于欧拉解集的概率密度估计,当1≤d≤4时,其密度分布函数F的核密度梯度估计
Figure BDA0003016673840000031
定义为:
Figure BDA0003016673840000032
Figure BDA0003016673840000033
式中,K为一维核密度基函数;
Kd,H为d个一维核密度基函数的内集;
χ为估计网格,可以写为(χ1,…,χd)T
Figure BDA0003016673840000034
为梯度算子;
H为对角带宽矩阵;
d为数据样本维度数;
对于重力或重力梯度张量得欧拉解而言,第i个独立同分布随机变量Xi可以写为xi,yi,zi,Ni的不同维度组合;
采用高斯核作为一维核密度基函数,则有:
Figure BDA0003016673840000035
Figure BDA0003016673840000036
式(7)中第j维度估计网格χj={χjk|k=1,…,Mj},j=1,…,d,及其相应带宽为
Figure BDA0003016673840000041
χjk为估计网格χ中第j维度方向χj的第k点;
对角带宽矩阵H,有:
Figure BDA0003016673840000042
式中,
Figure BDA0003016673840000043
std为标准差函数,iqr为四分位距函数,Sj
Figure BDA0003016673840000044
边际样本方差;
当r=2时,所述估计网格χ的相应带宽达到最大;当r=0时,所述估计网格χ的相应带宽达到最小。
可选地,所述将欧拉解数据作为独立同分布采样投影至所述估计网格,根据所述估计网络计算网格计数,包括:
将观测数据样本X投影至所述估计网格χ,并以观测数据样本Xi在所述估计网格χ各点上的权重u之和网格计数C,作为所述观测数据样本的近似;
设所述估计网格χ上存在一条边χjk,j(k+1),并将边χjk,j(k+1)简化为(χjkj(k+1));
若Xi正好位于由边
Figure BDA0003016673840000045
围成的,且含2d面和2d节点的超矩形中,则设定边χjk,j(k+1)的两个端点的序号分别为j-1和2d-1+j-1;其中,该超矩形的角标为l=(l1,…ld);
将Xi投影至边χjk,j(k+1)计算该两个端点的权重
Figure BDA0003016673840000046
Figure BDA0003016673840000047
则权重
Figure BDA0003016673840000048
可以定义为:
Figure BDA0003016673840000049
式中,Xij为Xi第j个变量,
Figure BDA00030166738400000410
向下取整运算符号;
则观测数据样本Xi与超矩形l的第q个节点的网格计数可以为:
Figure BDA00030166738400000411
可选地,所述基于所述网格计数与核函数的卷积确定不同维度数据的异常源空间位置及类型,包括:
基于所述网格计数确定多变量核密度估计值,则相应概率密度结果可以由下式计算得到:
Figure BDA0003016673840000051
将上式改写为卷积形式,则有:
Figure BDA0003016673840000052
Figure BDA0003016673840000053
式中,
Figure BDA0003016673840000054
采用傅立叶变换计算Cj-l
Figure BDA0003016673840000055
的卷积计算量
Figure BDA0003016673840000056
获得观测数据样本X于估计网格χ的概率密度分布
Figure BDA0003016673840000057
继而,通过所述概率密度分布确定不同维度数据的异常源空间位置及类型。
可选地,利用式(12)或式(13)计算得到观测数据样本X于估计网格χ上概率密度值
Figure BDA0003016673840000058
若所述观测数据样本X为一维数据,即d=1,继而通过概率密度曲线峰值与横轴的坐标确定一维数据的异常源空间位置及类型;
若所述观测数据样本X为二维数据,即d=2,通过概率密度图二维峰值与两横轴的坐标确定二维数据的异常源空间位置及类型;
若所述观测数据样本X为三维数据,即d=3,通过概率密度等值面确定三维数据的异常源空间位置及类型;
若所述观测数据样本X为高于三维的数据,即d>3,则将对应的概率密度结果降维至三维,并通过概率密度等值面确定三维数据的异常源空间位置及类型。
如上所述,本发明提供一种基于多变量核密度欧拉解概率密度成像方法,具有以下有益效果:
本发明充分利用欧拉反褶积优解能有效确定异常源中心,具有聚焦的特性,计算各欧拉解的概率密度,通过概率密度值的大小,从而达到快速定位异常源目的。相比于原有欧拉解剔除策略,该方法更易区分异常形态、判断异常位置。而且本发明不用剔除欧拉解缪解,减少了数据处理步骤。另外,本发明对于如{x}, {y},{z},{N}等一维欧拉解子集而言,可以极小的计算量,通过概率密度曲线的峰值于坐标横轴的位置,快速确定地下地质构造相应的空间位置和构造指数。对于如{x,y},{y,z},{x,z},{x,N}等二维欧拉解组合而言,通过概率密度图的二维峰值于坐标轴的位置,快速确定地下地质构造相应的空间位置和构造指数。对于如{x,y,z},{y,z,N},{x,y,N},{x,z,N}等三维欧拉解组合而言,通过三维概率密度等值面的数量级及等值面的多寡,快速确定地下地质构造相应的空间位置和构造指数。本发明在某一维度方向的估计网格χj采用等间距的估计点,可以将观测数据样本{Xj}整体投影至χj,从而易于构建矢量化(于MATLAB而言)/ 高效并行的算法(于C/C++等编程语言而言),进而快速获得网格计数C。同时,本发明与传统重力/重力梯度张量欧拉解解释技术难于剔除欧拉解缪解、确定欧拉解优解、分离欧拉解簇和标识欧拉反褶积解相比,本发明以各欧拉解的概率密度值为基础,可通过各概率密度的峰值点,进而快速有效地确定各异常源。
附图说明
图1为一实施例提供的基于多变量核密度欧拉解概率密度成像方法的流程示意图;
图2为一实施例提供的圆柱体正演含噪声结果示意图;
图3为一实施例提供的一维概率密度曲线{x}的结果示意图;
图4为一实施例提供的一维概率密度曲线{y}的结果示意图;
图5为一实施例提供的一维概率密度曲线{z}的结果示意图;
图6为一实施例提供的一维概率密度曲线{N}的结果示意图;
图7为一实施例提供的二维概率密度曲线{x,y}的结果示意图;
图8为一实施例提供的二维概率密度曲线{x,z}的结果示意图;
图9为一实施例提供的二维概率密度曲线{x,N}的结果示意图;
图10为一实施例提供的二维概率密度曲线{y,z}的结果示意图;
图11为一实施例提供的二维概率密度曲线{y,N}的结果示意图;
图12为一实施例提供的二维概率密度曲线{z,N}的结果示意图;
图13为一实施例提供的三维概率密度曲线{x,y,z}的结果示意图;
图14为一实施例提供的三维概率密度曲线{x,y,N}的结果示意图;
图15为一实施例提供的三维概率密度曲线{x,z,N}的结果示意图;
图16为一实施例提供的三维概率密度曲线{y,z,N}的结果示意图;
图17为一实施例提供的实际观测布格重力异常数据图;
图18为一实施例提供的第一待测区域中欧拉解及概率密度结果{x,y}的示意图;
图19为另一实施例提供的第一待测区域中欧拉解及概率密度结果{x,y}的示意图;
图20为一实施例提供的第一待测区域中欧拉解及概率密度结果{x,y,z}的对比示意图;
图21为另一实施例提供的第一待测区域中欧拉解及概率密度结果{x,y,z}对比示意图。
具体实施方式
以下通过特定的具体实例说明本发明的实施方式,本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用,本说明书中的各项细节也可以基于不同观点与应用,在没有背离本发明的精神下进行各种修饰或改变。需说明的是,在不冲突的情况下,以下实施例及实施例中的特征可以相互组合。
需要说明的是,以下实施例中所提供的图示仅以示意方式说明本发明的基本构想,遂图式中仅显示与本发明中有关的组件而非按照实际实施时的组件数目、形状及尺寸绘制,其实际实施时各组件的型态、数量及比例可为一种随意的改变,且其组件布局型态也可能更为复杂。
请参阅图1至图21所示,本发明提供一种基于多变量核密度欧拉解概率密度成像方法,包括:
根据勘探或研究工作要求,确定测区范围,根据测量设备条件,测量重力矢量数据及重力梯度张量数据,或通过离散余弦变换或傅里叶变换将重力数据换算到重力矢量数据及重力梯度张量数据。
根据所述重力矢量数据或所述重力梯度张量数据,利用公式(1-3)构建三维重力梯度张量欧拉反褶积方程;
确定滑动窗口大小,如wx=wy=3,并利用滑动窗口内的数据,根据公式(4) 构建线性方程组;
利用奇异值分解所述线性方程组,获取所述线性方程组的欧拉解 mi=[xi,yi,zi,Ni],并输出对应的欧拉解集{mi};其中,所述欧拉解包括异常源空间位置xi,yi,zi、异常源类型和奇异值分解误差Ni
在公式(7)计算前,将欧拉解解集拆分成不同维度的组合,如{xi},{xi,yi}, {xi,yi,zi}等;并根据所述欧拉解解集组合,利用公式(8—9)估算估计网格最大或最小带宽,及确定所述估计网格的大小;
根据公式(10)将欧拉解数据作为独立同分布采样投影至所述估计网格;根据所述估计网络,利用式(11)计算网格计数,基于公式(12)所述网格计数与核函数的卷积,计算不同维度欧拉解数据的概论密度值,进而确定异常源空间位置及类型。
可选地,所述根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程,包括:
Figure BDA0003016673840000081
Figure BDA0003016673840000091
Figure BDA0003016673840000092
式中,Bα为异常背景场,gα为重力梯度,gαβ为重力张量分量,{α,β∈x,y,z}; (x0,y0,z0)为测点位置;(x,y,z)为待求解位场异常源位置;构造指数N是位场异常随场源深度衰减变化陡缓的量度,即异常源类型。
可选地,所述确定滑动窗口大小,并利用滑动窗口内的数据构建线性方程组,以及利用奇异值分解所述线性方程组,获取所述线性方程组的欧拉解,并输出对应的欧拉解集,包括:
确定滑动窗口大小为Wx×Wy;
将滑窗中的Wx×Wy个测点数据代入(1)~(3)构建线性方程组,如下:
Am=b(4);
式中,算子A为
Figure BDA0003016673840000093
b为
Figure BDA0003016673840000094
m为欧拉反褶积解集;其中,第i个欧拉反褶积解可以表示为mi=[xi,yi,zi,Ni];
利用奇异值分解求取异常源空间位置[x,y,z]及构造指数N。
可选地,将欧拉解解集拆分成不同维度的组合,并根据所述欧拉解解集组合估算估计网格最大或最小带宽,及确定所述估计网格的大小,包括:
引入多维核密度估计,设独立同分布采样X=(X1,…,Xi,…,Xn)T为取值于d维实数空间的独立同分布随机变量;
对于欧拉解集的概率密度估计,当1≤d≤4时,其密度分布函数F的核密度梯度估计定义为:
Figure BDA0003016673840000095
Figure BDA0003016673840000096
式中,K为一维核密度基函数或多维核密度函数;
Kd,H为一维核密度基函数或多维核密度函数的内集;
χ为估计网格,可以写为(χ1,…,χd)T
Figure BDA0003016673840000101
为梯度算子;
H为对角带宽矩阵;
d为观测数据样本维度数;
对于重力或重力梯度张量得欧拉解而言,第i个独立同分布随机变量Xi可以写为xi,yi,zi,Ni的不同维度组合。
选取高斯核作为一维核密度基函数,则有:
Figure BDA0003016673840000102
Figure BDA0003016673840000103
式(7)中第j维度估计网格χj={χjk|k=1,…,Mj},j=1,…,d,及其相应带宽为
Figure BDA0003016673840000104
χjk为估计网格χ中第j维度方向χj的第k点;
对角带宽矩阵H,有:
Figure BDA0003016673840000105
式中,
Figure BDA0003016673840000106
std为标准差函数,iqr为四分位距函数,Sj
Figure BDA0003016673840000107
边际样本方差;
当r=2时,所述估计网格χ的相应带宽达到最大;当r=0时,所述估计网格χ的相应带宽达达到最小。
可选地,所述将不同组合欧拉解数据作为独立同分布采样X投影至所述估计网格χ,根据所述估计网络计算网格计数C,包括:
将观测数据样本X投影至所述估计网格χ,并以观测数据样本Xi在所述估计网格χ各点上的权重u计算网格计数C,并将此网格计数C作为所述观测数据样本的近似;
设所述估计网格χ上存在一条边χjk,j(k+1),并将边χjk,j(k+1)简化为(χjkj(k+1));
若Xi正好位于由边
Figure BDA0003016673840000111
围成的,且含2d面和2d节点的超矩形中,则设定边χjk,j(k+1)的两个端点的序号分别为j-1和2d-1+j-1;其中,该超矩形的角标为l=(l1,…ld);
将Xi投影至边χjk,j(k+1),以计算Xi在该边χjk,j(k+1)两个端点的权重
Figure BDA0003016673840000112
Figure BDA0003016673840000113
则权重
Figure BDA0003016673840000114
可以定义为:
Figure BDA0003016673840000115
式中,Xij为Xi第j个变量,
Figure BDA0003016673840000116
向下取整运算符号;
则观测数据样本Xi与超矩形l的第q个节点的网格计数可以为:
Figure BDA0003016673840000117
可选地,所述基于所述网格计数与核函数的卷积确定不同维度数据的异常源空间位置及类型,包括:
基于所述网格计数确定多变量核密度估计值,则相应概率密度结果可以由下式计算得到:
Figure BDA0003016673840000118
将上式变为卷积形式,则有:
Figure BDA0003016673840000119
Figure BDA00030166738400001110
式中,
Figure BDA00030166738400001111
采用傅立叶变换计算Cj-l
Figure BDA00030166738400001112
的卷积计算量
Figure BDA00030166738400001113
获得观测数据样本X于估计网格χ的概率密度分布
Figure BDA00030166738400001114
继而,通过所述概率密度分布确定不同维度数据的异常源空间位置及类型。
可选地,若所述观测数据样本X为一维数据,即d=1,通过概率密度曲线峰值与横轴的坐标确定一维数据的异常源空间位置及类型;
若所述观测数据样本X为二维数据,即d=2,通过概率密度图二维峰值与两横轴的坐标确定二维数据的异常源空间位置及类型;
若所述观测数据样本X为三维数据,即d=3,通过概率密度等值面确定三维数据的异常源空间位置及类型;
若所述观测数据样本X为高于三维的数据,即d>3,则将对应的概率密度结果降维至三维,并通过概率密度等值面确定三维数据的异常源空间位置及类型。
上述实施例仅例示性说明本发明的原理及其功效,而非用于限制本发明。任何熟悉此技术的人士皆可在不违背本发明的精神及范畴下,对上述实施例进行修饰或改变。因此,举凡所属技术领域中具有通常知识者在未脱离本发明所揭示的精神与技术思想下所完成的一切等效修饰或改变,仍应由本发明的权利要求所涵盖。

Claims (7)

1.一种基于多变量核密度欧拉解概率密度成像方法,其特征在于,包括以下步骤:
确定待测区域范围,并测量重力矢量数据及重力梯度张量数据;或测量重力数据,通过离散余弦变换或傅里叶变换将重力数据换算到重力矢量数据及重力梯度张量数据;
根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程;
确定滑动窗口大小,并利用滑动窗口内的数据构建线性方程组;
利用奇异值分解所述线性方程组,获取所述线性方程组的欧拉解,并输出对应的欧拉解集;其中,所述欧拉解包括异常源空间位置、异常源类型和奇异值分解误差;
将欧拉解解集拆分成不同维度的组合,并根据所述欧拉解解集组合估算估计网格最大或最小带宽,及确定所述估计网格的大小;
将欧拉解数据作为独立同分布采样投影至所述估计网格,根据所述估计网格计算网格计数,并基于所述网格计数与核函数的卷积确定不同维度数据的异常源空间位置及类型。
2.根据权利要求1所述的基于多变量核密度欧拉解概率密度成像方法,其特征在于,所述根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程,包括:
Figure FDA0003506768210000011
Figure FDA0003506768210000012
Figure FDA0003506768210000013
式中,Bα为异常背景场,gα为重力梯度,{α∈x,y,z};(x0,y0,z0)为观测点位置坐标;(x,y,z)为待求解位场异常源位置;构造指数N是位场异常随场源深度衰减变化陡缓的量度,即异常源类型。
3.根据权利要求2所述的基于多变量核密度欧拉解概率密度成像方法,其特征在于,所述确定滑动窗口大小,并利用滑动窗口内的数据构建线性方程组,以及利用奇异值分解所述线性方程组,获取所述线性方程组的欧拉解,并输出对应的欧拉解集,包括:
确定滑动窗口大小为Wx×Wy;
将滑窗中的Wx×Wy个测点数据代入(1)~(3)构建线性方程组,如下:
Am=b (4);
式中,算子A为
Figure FDA0003506768210000021
b为
Figure FDA0003506768210000022
m为欧拉反褶积解集;其中,第i个欧拉反褶积解表示为mi=[xi,yi,zi,Ni];
利用奇异值分解求取异常源空间位置[x,y,z]及构造指数N;
逐一移动滑动窗口,构建线性方程组,并获得相应滑动窗口的欧拉解mi,直至滑动窗口覆盖整个测区,构建欧拉解集{mi}。
4.根据权利要求3所述的基于多变量核密度欧拉解概率密度成像方法,其特征在于,将欧拉解解集拆分成不同维度的组合,并根据不同欧拉解解集组合估算估计网格最大或最小带宽,及确定所述估计网格的大小,包括:
引入多维核密度估计,设独立同分布采样X=(X1,…,Xi,…,Xn)T为取值于d维实数空间的独立同分布随机变量;
对于欧拉解集的概率密度估计,当1≤d≤4时,其密度分布函数F的核密度梯度估计
Figure FDA0003506768210000023
定义为:
Figure FDA0003506768210000024
Figure FDA0003506768210000025
式中,K为一维核密度基函数;
Kd,H为d个一维核密度基函数的内积;
χ为估计网格,写为(χ1,…,χd)T
Figure FDA0003506768210000031
为梯度算子;
H为对角带宽矩阵;
d为数据样本维度数;对于重力或重力梯度张量得欧拉解而言,第i个独立同分布随机变量Xi写为xi,yi,zi,Ni的不同维度组合;
采用高斯核作为一维核密度基函数,则有:
Figure FDA0003506768210000032
Figure FDA0003506768210000033
式(7)中第j维度估计网格χj={χjk|k=1,…,Mj},j=1,…,d,及其相应带宽为
Figure FDA0003506768210000034
χjk为估计网格χ中第j维度方向χj的第k点;
对角带宽矩阵H,有:
Figure FDA0003506768210000035
式中,
Figure FDA0003506768210000036
std为标准差函数,iqr为四分位距函数,Sj
Figure FDA0003506768210000037
边际样本方差;Xij为Xi第j个变量;
当r=2时,所述估计网格χ的相应带宽达到最大;当r=0时,所述估计网格χ的相应带宽达到最小。
5.根据权利要求4所述的基于多变量核密度欧拉解概率密度成像方法,其特征在于,所述将欧拉解数据作为独立同分布采样投影至所述估计网格,根据所述估计网格计算网格计数,包括:
将观测数据样本X投影至所述估计网格χ,并以观测数据样本Xi在所述估计网格χ各点上的权重u之和网格计数C,作为所述观测数据样本的近似;
设所述估计网格χ上存在一条边χjk,j(k+1),并将边χjk,j(k+1)简化为(χjkj(k+1));
若Xi正好位于由边
Figure FDA0003506768210000041
围成的,且含2d面和2d节点的超矩形中,则设定边χjk,j(k+1)的两个端点的序号分别为j-1和2d-1+j-1;其中,该超矩形的角标为l=(l1,…ld);
将Xi投影至边χjk,j(k+1)计算该两个端点的权重
Figure FDA0003506768210000042
Figure FDA0003506768210000043
则权重
Figure FDA0003506768210000044
定义为:
Figure FDA0003506768210000045
式中,Xij为Xi第j个变量,
Figure FDA0003506768210000046
向下取整运算符号;
则观测数据样本Xi与超矩形l的第q个节点的网格计数为:
Figure FDA0003506768210000047
6.根据权利要求5所述的基于多变量核密度欧拉解概率密度成像方法,其特征在于,所述基于所述网格计数与核函数的卷积确定不同维度数据的异常源空间位置及类型,包括:
基于所述网格计数确定多变量核密度估计值,则相应概率密度结果由下式计算得到:
Figure FDA0003506768210000048
将上式改写为卷积形式,则有:
Figure FDA0003506768210000049
Figure FDA00035067682100000410
式中,
Figure FDA00035067682100000411
采用傅立叶变换计算Cj-l
Figure FDA00035067682100000412
的卷积计算量
Figure FDA00035067682100000413
获得观测数据样本X于估计网格χ的概率密度分布
Figure FDA00035067682100000414
继而,通过所述概率密度分布确定不同维度数据的异常源空间位置及类型。
7.根据权利要求6所述的基于多变量核密度欧拉解概率密度成像方法,其特征在于,利用式(12)或式(13)计算得到观测数据样本X于估计网格χ上概率密度值
Figure FDA0003506768210000051
若所述观测数据样本X为一维数据,即d=1,继而通过概率密度曲线峰值与横轴的坐标确定一维数据的异常源空间位置及类型;
若所述观测数据样本X为二维数据,即d=2,通过概率密度图二维峰值与两横轴的坐标确定二维数据的异常源空间位置及类型;
若所述观测数据样本X为三维数据,即d=3,通过概率密度等值面确定三维数据的异常源空间位置及类型;
若所述观测数据样本X为高于三维的数据,即d>3,则将对应的概率密度结果降维至三维,并通过概率密度等值面确定三维数据的异常源空间位置及类型。
CN202110390937.4A 2021-04-12 2021-04-12 一种基于多变量核密度欧拉解概率密度成像方法 Active CN113138426B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110390937.4A CN113138426B (zh) 2021-04-12 2021-04-12 一种基于多变量核密度欧拉解概率密度成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110390937.4A CN113138426B (zh) 2021-04-12 2021-04-12 一种基于多变量核密度欧拉解概率密度成像方法

Publications (2)

Publication Number Publication Date
CN113138426A CN113138426A (zh) 2021-07-20
CN113138426B true CN113138426B (zh) 2022-04-01

Family

ID=76811179

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110390937.4A Active CN113138426B (zh) 2021-04-12 2021-04-12 一种基于多变量核密度欧拉解概率密度成像方法

Country Status (1)

Country Link
CN (1) CN113138426B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113806686B (zh) * 2021-11-19 2022-02-08 中南大学 大规模复杂地质体重力梯度快速计算方法、装置和设备

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20080023946A (ko) * 2006-09-12 2008-03-17 한국지질자원연구원 오일러 해를 이용한 지하공동의 3차원 중력역산 방법 및이를 이용한 3차원 영상화 방법
CN108732622A (zh) * 2018-05-18 2018-11-02 吉林大学 一种不同高度数据融合联合反演地质体几何形态的方法
CN112462442A (zh) * 2020-11-30 2021-03-09 山东大学 重磁位场场源位置估计方法、系统、介质及电子设备

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130018588A1 (en) * 2011-07-11 2013-01-17 Technolmaging, Llc. Method of real time subsurface imaging using gravity and/or magnetic data measured from a moving platform

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20080023946A (ko) * 2006-09-12 2008-03-17 한국지질자원연구원 오일러 해를 이용한 지하공동의 3차원 중력역산 방법 및이를 이용한 3차원 영상화 방법
CN108732622A (zh) * 2018-05-18 2018-11-02 吉林大学 一种不同高度数据融合联合反演地质体几何形态的方法
CN112462442A (zh) * 2020-11-30 2021-03-09 山东大学 重磁位场场源位置估计方法、系统、介质及电子设备

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
基于自适应带宽的快速动态高斯核均值漂移算法;周芳芳 等;《控制理论与应用》;20080831;第25卷(第4期);第608-612页 *
基于自适应模糊聚类分析的重力张量欧拉反褶积解;曹书锦 等;《中南大学学报(自然科学版)》;20120331;第43卷(第3期);第1033-1039页 *
重力场欧拉反褶积最优解提取;周勇 等;《物探化探计算技术》;20170930;第39卷(第5期);第598-604页 *

Also Published As

Publication number Publication date
CN113138426A (zh) 2021-07-20

Similar Documents

Publication Publication Date Title
CN105427300B (zh) 一种基于低秩表示和学习字典的高光谱图像异常探测方法
Sauber et al. Multifield-graphs: An approach to visualizing correlations in multifield scalar data
Nascimento et al. Vertex component analysis: A fast algorithm to unmix hyperspectral data
Altmann et al. Unsupervised post-nonlinear unmixing of hyperspectral images using a Hamiltonian Monte Carlo algorithm
Bardsley et al. Randomize-then-optimize: A method for sampling from posterior distributions in nonlinear inverse problems
Plaza et al. Spatial/spectral endmember extraction by multidimensional morphological operations
AU2011248992B2 (en) Windowed statistical analysis for anomaly detection in geophysical datasets
Cardiel et al. Using spectroscopic data to disentangle stellar population properties
Nascimento et al. Comparing edge detection methods based on stochastic entropies and distances for PolSAR imagery
CN103955926B (zh) 基于Semi-NMF的遥感图像变化检测方法
Bremer et al. Derivation of tree skeletons and error assessment using LiDAR point cloud data of varying quality
Bourennane et al. Improvement of target-detection algorithms based on adaptive three-dimensional filtering
CN113138426B (zh) 一种基于多变量核密度欧拉解概率密度成像方法
Jerez et al. Fast target detection via template matching in compressive phase retrieval
Morháč et al. Complete positive deconvolution of spectrometric data
Kappel MSR, a multi-spectrum retrieval technique for spatially-temporally correlated or common Venus surface and atmosphere parameters
Tashlinskii et al. Analysis of methods of estimating objective function gradient during recurrent measurements of image parameters
Uribe et al. A hybrid Gibbs sampler for edge-preserving tomographic reconstruction with uncertain view angles
CN113011086B (zh) 一种基于ga-svr算法森林生物量的估测方法
Tate Estimating the fractal dimension of synthetic topographic surfaces
US11634175B2 (en) Dip angle-steering median filtering method based on a niche differential evolution algorithm
US6961677B1 (en) Method and apparatus for categorizing unexplained residuals
CN105571716B (zh) 一种基于差分与卷积核的线采样高光谱数据目标探测方法
Clifford et al. A process to colorize and assess visualizations of noisy x-ray computed tomography hyperspectral data of materials with similar spectral signatures
CN113177306B (zh) 一种基于b样条函数的欧拉解集概率密度成像方法

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