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

Abstract

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

Description

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

技术领域

本发明涉及重力勘探技术领域，特别是涉及一种基于多变量核密度欧拉解概率密度成像方法。

背景技术

欧拉反褶积以Euler齐次方程为理论基础，利用总场及其梯度等位场数据，只需事先确定与场源性质有关的构造指数或枚举构造指数，便可快速有效地圈出异常源的基本轮廓，尤其适合于大面积位场数据的分析和解释。然而，使用欧拉反褶积对异常进行推断时，传统欧拉解图示方法一般使用色阶或不同大小圆来标定深度或构造指数，即色谱。但其难于标示异常源间及各单一欧拉解间的差异及相关性，尤其当欧拉解发散时，如对于在一处大量聚集(如欧拉解聚集程度大的地方)和在它处成零星分布(如多异常源间的伪迹)的两类欧拉解集，其色谱的空间分布上几乎没有差别。因而，如何在众多的欧拉解中提取优解、剔除谬解、分离欧拉解簇继而标定异常源一直是困扰人们的一个难题。

发明内容

鉴于以上所述现有技术的缺点，本发明的目的在于提供一种基于多变量核密度欧拉解概率密度成像方法，用于解决现有技术中存在的问题。

为实现上述目的及其他相关目的，本发明提供一种基于多变量核密度欧拉解概率密度成像方法，包括以下步骤：

确定待测区域范围；并测量重力矢量数据及重力梯度张量数据，或测量重力数据，通过离散余弦变换或傅里叶变换将重力数据换算到重力矢量数据及重力梯度张量数据；

根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程；

确定滑动窗口大小，并利用滑动窗口内的数据构建线性方程组；

利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解，并输出对应的欧拉解集；其中，所述欧拉解包括异常源空间位置、异常源类型和奇异值分解误差；

将欧拉解解集拆分成不同维度的组合，并根据所述欧拉解解集组合估算估计网格最大或最小带宽，及确定所述估计网格的大小；

将欧拉解数据作为独立同分布采样投影至所述估计网格，根据所述估计网络计算网格计数，并基于所述网格计数与核函数的卷积确定不同维度数据的异常源空间位置及类型。

可选地，所述根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程，包括：

式中，B_α为异常背景场，g_α为重力梯度，g_αβ为重力张量分量，{α,β∈x,y,z}； (x₀,y₀,z₀)为测点位置；(x,y,z)为待求解位场异常源位置；构造指数N是位场异常随场源深度衰减变化陡缓的量度，即异常源类型。

可选地，所述确定滑动窗口大小，并利用滑动窗口内的数据构建线性方程组，以及利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解，并输出对应的欧拉解集，包括：

确定滑动窗口大小为Wx×Wy；

将滑窗中的Wx×Wy个测点数据代入(1)～(3)构建线性方程组，如下：

Am＝b (4)；

式中，算子A为

b为

m为欧拉反褶积解集；其中，第i个欧拉反褶积解可以表示为m_i＝[x_i,y_i,z_i,N_i]；

利用奇异值分解求取异常源空间位置[x，y，z]及构造指数N。

逐一移动滑动窗口，构建线性方程组，并获得相应滑动窗口的欧拉解m_i，直至滑动窗口覆盖整个测区，构建欧拉解集{m_i}。

可选地，将欧拉解解集拆分成不同维度的组合，并根据所述不同欧拉解解集组合估算估计网格最大或最小带宽，及确定所述估计网格的大小，包括：

引入多维核密度估计，设独立同分布采样X＝(X₁,…,X_i,…,X_n)^T为取值于d维实数空间的独立同分布随机变量；

对于欧拉解集的概率密度估计，当1≤d≤4时，其密度分布函数F的核密度梯度估计

定义为：

式中，K为一维核密度基函数；

K_d,H为d个一维核密度基函数的内集；

χ为估计网格，可以写为(χ₁,…,χ_d)^T；

为梯度算子；

H为对角带宽矩阵；

d为数据样本维度数；

对于重力或重力梯度张量得欧拉解而言，第i个独立同分布随机变量X_i可以写为xⁱ,yⁱ,zⁱ,Nⁱ的不同维度组合；

采用高斯核作为一维核密度基函数，则有：

且

式(7)中第j维度估计网格χ_j＝{χ_jk|k＝1,…,M_j}，j＝1,…,d，及其相应带宽为

χ_jk为估计网格χ中第j维度方向χ_j的第k点；

对角带宽矩阵H，有：

式中，

std为标准差函数，iqr为四分位距函数，S_j为

边际样本方差；

当r＝2时，所述估计网格χ的相应带宽达到最大；当r＝0时，所述估计网格χ的相应带宽达到最小。

可选地，所述将欧拉解数据作为独立同分布采样投影至所述估计网格，根据所述估计网络计算网格计数，包括：

将观测数据样本X投影至所述估计网格χ，并以观测数据样本X_i在所述估计网格χ各点上的权重u之和网格计数C，作为所述观测数据样本的近似；

设所述估计网格χ上存在一条边χ_jk,j(k+1)，并将边χ_jk,j(k+1)简化为(χ_jk,χ_j(k+1))；

若X_i正好位于由边

围成的，且含2^d面和2^d节点的超矩形中，则设定边χ_jk,j(k+1)的两个端点的序号分别为j-1和2^d-1+j-1；其中，该超矩形的角标为l＝(l₁,…l_d)；

将X_i投影至边χ_jk,j(k+1)计算该两个端点的权重

和

则权重

可以定义为：

式中，X_ij为X_i第j个变量，

向下取整运算符号；

则观测数据样本X_i与超矩形l的第q个节点的网格计数可以为：

可选地，所述基于所述网格计数与核函数的卷积确定不同维度数据的异常源空间位置及类型，包括：

基于所述网格计数确定多变量核密度估计值，则相应概率密度结果可以由下式计算得到：

将上式改写为卷积形式，则有：

式中，

采用傅立叶变换计算C_j-l与

的卷积计算量

获得观测数据样本X于估计网格χ的概率密度分布

继而，通过所述概率密度分布确定不同维度数据的异常源空间位置及类型。

可选地，利用式(12)或式(13)计算得到观测数据样本X于估计网格χ上概率密度值

若所述观测数据样本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采用等间距的估计点，可以将观测数据样本{X_j}整体投影至χ_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)构建三维重力梯度张量欧拉反褶积方程；

确定滑动窗口大小，如w_x＝w_y＝3，并利用滑动窗口内的数据，根据公式(4) 构建线性方程组；

利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解 mⁱ＝[xⁱ,yⁱ,zⁱ,Nⁱ]，并输出对应的欧拉解集{mⁱ}；其中，所述欧拉解包括异常源空间位置xⁱ,yⁱ,zⁱ、异常源类型和奇异值分解误差Nⁱ；

在公式(7)计算前，将欧拉解解集拆分成不同维度的组合，如{xⁱ}，{xⁱ,yⁱ}， {xⁱ,yⁱ,zⁱ}等；并根据所述欧拉解解集组合，利用公式(8—9)估算估计网格最大或最小带宽，及确定所述估计网格的大小；

根据公式(10)将欧拉解数据作为独立同分布采样投影至所述估计网格；根据所述估计网络，利用式(11)计算网格计数，基于公式(12)所述网格计数与核函数的卷积，计算不同维度欧拉解数据的概论密度值，进而确定异常源空间位置及类型。

可选地，所述根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程，包括：

式中，B_α为异常背景场，g_α为重力梯度，g_αβ为重力张量分量，{α,β∈x,y,z}； (x₀,y₀,z₀)为测点位置；(x,y,z)为待求解位场异常源位置；构造指数N是位场异常随场源深度衰减变化陡缓的量度，即异常源类型。

可选地，所述确定滑动窗口大小，并利用滑动窗口内的数据构建线性方程组，以及利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解，并输出对应的欧拉解集，包括：

确定滑动窗口大小为Wx×Wy；

将滑窗中的Wx×Wy个测点数据代入(1)～(3)构建线性方程组，如下：

Am＝b(4)；

式中，算子A为

b为

m为欧拉反褶积解集；其中，第i个欧拉反褶积解可以表示为mⁱ＝[xⁱ,yⁱ,zⁱ,Nⁱ]；

利用奇异值分解求取异常源空间位置[x，y，z]及构造指数N。

可选地，将欧拉解解集拆分成不同维度的组合，并根据所述欧拉解解集组合估算估计网格最大或最小带宽，及确定所述估计网格的大小，包括：

引入多维核密度估计，设独立同分布采样X＝(X₁,…,X_i,…,X_n)^T为取值于d维实数空间的独立同分布随机变量；

对于欧拉解集的概率密度估计，当1≤d≤4时，其密度分布函数F的核密度梯度估计定义为：

式中，K为一维核密度基函数或多维核密度函数；

K_d,H为一维核密度基函数或多维核密度函数的内集；

χ为估计网格，可以写为(χ₁,…,χ_d)^T；

为梯度算子；

H为对角带宽矩阵；

d为观测数据样本维度数；

对于重力或重力梯度张量得欧拉解而言，第i个独立同分布随机变量X_i可以写为xⁱ,yⁱ,zⁱ,Nⁱ的不同维度组合。

选取高斯核作为一维核密度基函数，则有：

式(7)中第j维度估计网格χ_j＝{χ_jk|k＝1,…,M_j}，j＝1,…,d，及其相应带宽为

χ_jk为估计网格χ中第j维度方向χ_j的第k点；

对角带宽矩阵H，有：

式中，

std为标准差函数，iqr为四分位距函数，S_j为

边际样本方差；

当r＝2时，所述估计网格χ的相应带宽达到最大；当r＝0时，所述估计网格χ的相应带宽达达到最小。

可选地，所述将不同组合欧拉解数据作为独立同分布采样X投影至所述估计网格χ，根据所述估计网络计算网格计数C，包括：

将观测数据样本X投影至所述估计网格χ，并以观测数据样本X_i在所述估计网格χ各点上的权重u计算网格计数C，并将此网格计数C作为所述观测数据样本的近似；

设所述估计网格χ上存在一条边χ_jk,j(k+1)，并将边χ_jk,j(k+1)简化为(χ_jk,χ_j(k+1))；

若X_i正好位于由边

围成的，且含2^d面和2^d节点的超矩形中，则设定边χ_jk,j(k+1)的两个端点的序号分别为j-1和2^d-1+j-1；其中，该超矩形的角标为l＝(l₁,…l_d)；

将X_i投影至边χ_jk,j(k+1)，以计算X_i在该边χ_jk,j(k+1)两个端点的权重

和

则权重

可以定义为：

式中，X_ij为X_i第j个变量，

向下取整运算符号；

则观测数据样本X_i与超矩形l的第q个节点的网格计数可以为：

可选地，所述基于所述网格计数与核函数的卷积确定不同维度数据的异常源空间位置及类型，包括：

基于所述网格计数确定多变量核密度估计值，则相应概率密度结果可以由下式计算得到：

将上式变为卷积形式，则有：

式中，

采用傅立叶变换计算C_j-l与

的卷积计算量

获得观测数据样本X于估计网格χ的概率密度分布

继而，通过所述概率密度分布确定不同维度数据的异常源空间位置及类型。

可选地，若所述观测数据样本X为一维数据，即d＝1，通过概率密度曲线峰值与横轴的坐标确定一维数据的异常源空间位置及类型；

若所述观测数据样本X为二维数据，即d＝2，通过概率密度图二维峰值与两横轴的坐标确定二维数据的异常源空间位置及类型；

若所述观测数据样本X为三维数据，即d＝3，通过概率密度等值面确定三维数据的异常源空间位置及类型；

若所述观测数据样本X为高于三维的数据，即d＞3，则将对应的概率密度结果降维至三维，并通过概率密度等值面确定三维数据的异常源空间位置及类型。

上述实施例仅例示性说明本发明的原理及其功效，而非用于限制本发明。任何熟悉此技术的人士皆可在不违背本发明的精神及范畴下，对上述实施例进行修饰或改变。因此，举凡所属技术领域中具有通常知识者在未脱离本发明所揭示的精神与技术思想下所完成的一切等效修饰或改变，仍应由本发明的权利要求所涵盖。

Claims (7)

1.一种基于多变量核密度欧拉解概率密度成像方法，其特征在于，包括以下步骤：

确定待测区域范围，并测量重力矢量数据及重力梯度张量数据；或测量重力数据，通过离散余弦变换或傅里叶变换将重力数据换算到重力矢量数据及重力梯度张量数据；

根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程；

确定滑动窗口大小，并利用滑动窗口内的数据构建线性方程组；

利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解，并输出对应的欧拉解集；其中，所述欧拉解包括异常源空间位置、异常源类型和奇异值分解误差；

将欧拉解解集拆分成不同维度的组合，并根据所述欧拉解解集组合估算估计网格最大或最小带宽，及确定所述估计网格的大小；

将欧拉解数据作为独立同分布采样投影至所述估计网格，根据所述估计网格计算网格计数，并基于所述网格计数与核函数的卷积确定不同维度数据的异常源空间位置及类型。

2.根据权利要求1所述的基于多变量核密度欧拉解概率密度成像方法，其特征在于，所述根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程，包括：

式中，B_α为异常背景场，g_α为重力梯度，{α∈x,y,z}；(x₀,y₀,z₀)为观测点位置坐标；(x,y,z)为待求解位场异常源位置；构造指数N是位场异常随场源深度衰减变化陡缓的量度，即异常源类型。

3.根据权利要求2所述的基于多变量核密度欧拉解概率密度成像方法，其特征在于，所述确定滑动窗口大小，并利用滑动窗口内的数据构建线性方程组，以及利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解，并输出对应的欧拉解集，包括：

确定滑动窗口大小为Wx×Wy；

将滑窗中的Wx×Wy个测点数据代入(1)～(3)构建线性方程组，如下：

Am＝b (4)；

式中，算子A为

b为

m为欧拉反褶积解集；其中，第i个欧拉反褶积解表示为m_i＝[x_i,y_i,z_i,N_i]；

利用奇异值分解求取异常源空间位置[x，y，z]及构造指数N；

逐一移动滑动窗口，构建线性方程组，并获得相应滑动窗口的欧拉解m_i，直至滑动窗口覆盖整个测区，构建欧拉解集{m_i}。

4.根据权利要求3所述的基于多变量核密度欧拉解概率密度成像方法，其特征在于，将欧拉解解集拆分成不同维度的组合，并根据不同欧拉解解集组合估算估计网格最大或最小带宽，及确定所述估计网格的大小，包括：

引入多维核密度估计，设独立同分布采样X＝(X₁,…,X_i,…,X_n)^T为取值于d维实数空间的独立同分布随机变量；

对于欧拉解集的概率密度估计，当1≤d≤4时，其密度分布函数F的核密度梯度估计

定义为：

式中，K为一维核密度基函数；

K_d,H为d个一维核密度基函数的内积；

χ为估计网格，写为(χ₁,…,χ_d)^T；

为梯度算子；

H为对角带宽矩阵；

d为数据样本维度数；对于重力或重力梯度张量得欧拉解而言，第i个独立同分布随机变量X_i写为xⁱ,yⁱ,zⁱ,Nⁱ的不同维度组合；

采用高斯核作为一维核密度基函数，则有：

且

式(7)中第j维度估计网格χ_j＝{χ_jk|k＝1,…,M_j}，j＝1,…,d，及其相应带宽为

χ_jk为估计网格χ中第j维度方向χ_j的第k点；

对角带宽矩阵H，有：

式中，

std为标准差函数，iqr为四分位距函数，S_j为

边际样本方差；X_ij为X_i第j个变量；

当r＝2时，所述估计网格χ的相应带宽达到最大；当r＝0时，所述估计网格χ的相应带宽达到最小。

5.根据权利要求4所述的基于多变量核密度欧拉解概率密度成像方法，其特征在于，所述将欧拉解数据作为独立同分布采样投影至所述估计网格，根据所述估计网格计算网格计数，包括：

将观测数据样本X投影至所述估计网格χ，并以观测数据样本X_i在所述估计网格χ各点上的权重u之和网格计数C，作为所述观测数据样本的近似；

设所述估计网格χ上存在一条边χ_jk,j(k+1)，并将边χ_jk,j(k+1)简化为(χ_jk,χ_j(k+1))；

若X_i正好位于由边

围成的，且含2^d面和2^d节点的超矩形中，则设定边χ_jk,j(k+1)的两个端点的序号分别为j-1和2^d-1+j-1；其中，该超矩形的角标为l＝(l₁,…l_d)；

将X_i投影至边χ_jk,j(k+1)计算该两个端点的权重

和

则权重

定义为：

式中，X_ij为X_i第j个变量，

向下取整运算符号；

则观测数据样本X_i与超矩形l的第q个节点的网格计数为：

6.根据权利要求5所述的基于多变量核密度欧拉解概率密度成像方法，其特征在于，所述基于所述网格计数与核函数的卷积确定不同维度数据的异常源空间位置及类型，包括：

基于所述网格计数确定多变量核密度估计值，则相应概率密度结果由下式计算得到：

将上式改写为卷积形式，则有：

式中，

采用傅立叶变换计算C_j-l与

的卷积计算量

获得观测数据样本X于估计网格χ的概率密度分布

继而，通过所述概率密度分布确定不同维度数据的异常源空间位置及类型。

7.根据权利要求6所述的基于多变量核密度欧拉解概率密度成像方法，其特征在于，利用式(12)或式(13)计算得到观测数据样本X于估计网格χ上概率密度值

若所述观测数据样本X为一维数据，即d＝1，继而通过概率密度曲线峰值与横轴的坐标确定一维数据的异常源空间位置及类型；

若所述观测数据样本X为二维数据，即d＝2，通过概率密度图二维峰值与两横轴的坐标确定二维数据的异常源空间位置及类型；

若所述观测数据样本X为三维数据，即d＝3，通过概率密度等值面确定三维数据的异常源空间位置及类型；

若所述观测数据样本X为高于三维的数据，即d＞3，则将对应的概率密度结果降维至三维，并通过概率密度等值面确定三维数据的异常源空间位置及类型。

Images