基于傅立叶变换的B样条函数欧拉解概率密度成像方法

技术领域

本发明涉及重力勘探技术领域，特别是涉及一种基于傅立叶变换的B样条函数欧拉解概率密度成像方法。

背景技术

欧拉反褶积是适合于大面积位场数据的一种半自动/自动估算场源位置的位场数据解释方法。源于欧拉超定方程组条件数大，致使欧拉解扰动大，且往往具有发散的趋势。因而，很多学者以欧拉超定方程组的标准差或截断误差作为标准，以过滤欧拉解集中的谬解；一些学者增加附加约束方程或约束条件，如利用重力梯度张量多分量含丰富异常信息(Euler deconvolution of gravity tensor gradient data)及利用“Worming”与欧拉方程之间的关系，构建混合欧拉反褶积方法以获取优解；或通过确定欧拉解集所构成的簇于欧拉解集关系，以剔除不隶属于各簇的欧拉解。然而对于相邻较近的异常源而言，传统簇分析方法难于区分各欧拉解簇。

发明内容

鉴于以上所述现有技术的缺点，本发明的目的在于提供一种基于快速傅立叶变换的B样条概率密度估计方法对欧拉解集进行处理，通过引入分箱近似方法，将样本数据快速投影至估计网格，以借助快速傅里叶变换实现估计网格与密度函数的离散卷积，在避免当数据样本过大或/和估计网格过大时，遍历估计网格上的每个节点是一个计算量巨大且伴随内存消耗过大的同时，获得较为聚焦的概率密度估计结果，以标示各异常源。

为实现上述目的及其他相关目的，本发明提供一种基于傅立叶变换的B样条函数欧拉解概率密度成像方法，包括有：

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

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

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

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

将欧拉解解集拆分成不同维度的组合，并根据所述欧拉解解集估算估计区间的上/下界及估计网格带宽，以及确定所述估计网格的大小；

构建基于B样条函数的多变量概率密度估计，并将欧拉解数据作为独立同分布采样投影至所述估计网格，计算网格计数；

对所述B样条函数和所述网格计数进行卷积，获取欧拉解概率密度分布结果，并根据所述欧拉解概率密度分布结果确定不同维度数据的异常源空间位置及类型。

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

式中，α＝x,y,z，B

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

确定滑动窗口大小为w＝w

利用公式(1)和滑动窗口构建矩形形式线性方程组：

式中，

n

欧拉反褶积解集为m＝{m

第i个欧拉反褶积解m

可选地，所述构建基于B样条函数的多变量概率密度估计，并将欧拉解数据作为独立同分布采样投影至所述估计网格，计算网格计数的过程包括：

设定X

将样本X

式中系数

基于近似分箱技术的多变量概率密度估计方法重新构建公式(3)，有：

式中，

可选地，将公式(5)中的C

若X

其中，

可选地，设定公式(4)中的权重

可选地，将公式(5)写为卷积的形式，则有：

式中，

对矩阵C和r进行补零及环绕排序，使C和r具有相同的大小，则有：

其中，C和r具有相同的维度λ

可选地，根据离散卷积定理，令：

其中，F代表傅里叶正变换，F

待求的密度估计

式中，S[a:b，c:d]表示矩阵S从a到b的行的子集和从c到d的列的子集。

如上所述，本发明提供一种基于傅立叶变换的B样条函数欧拉解概率密度成像方法，具有以下有益效果：

本发明采用基于规范化B样条的非参数概率密度估计方法，以各个欧拉解间的相似性及其集聚程度，分离各欧拉解簇以标示各异常源。针对B样条概率密度估计方法存在不聚焦及遍历估计网格计算量巨大等问题，本发明充分利用估计网格同一维度方向具有相同的间隔之一特点，构建了基于B样条基函数的多变量概率估计方法，并借助快速傅里叶变换实现估计网格与B样条基函数的离散卷积，计算各欧拉解的概率密度，通过概率密度值的大小，从而达到快速定位异常源目的。相比于原有欧拉解剔除策略及B样条概率密度估计方法，本发明的计算更为高效，且更易区分异常形态、判断异常位置。

附图说明

图1为一实施例提供的基于快速傅立叶变换的B样条函数重力欧拉解概率密度成像方法的流程示意图；

图2为一实施例提供的一维验证结果示意图；

图3为一实施例提供的无噪声的立方体重力正演图；

图4为一实施例提供的含3％噪声的立体重力正演图；

图5为一实施例提供的欧拉反褶积散点示意图；

图6为另一实施例提供的欧拉反褶积散点示意图；

图7为一实施例提供的基于快速傅立叶变换的一维欧拉解集{x}、{y}、{z}的一维B样条密度估计的概率密度曲线图；

图8为一实施例提供的基于快速傅立叶变换的一维欧拉解集{N}的一维B样条密度估计的概率密度曲线图；

图9为一实施例提供的单一立方体模型的二维欧拉解集{x,y}的二维B样条密度估计的概率密度分布图；

图10为一实施例提供的单一立方体模型的二维欧拉解集{x,z}的二维B样条密度估计的概率密度分布图；

图11为一实施例提供的单一立方体模型的二维欧拉解集{x,N}的二维B样条密度估计的概率密度分布图；

图12为一实施例提供的单一立方体模型的二维欧拉解集{y,z}的二维B样条密度估计的概率密度分布图；

图13为一实施例提供的单一立方体模型的二维欧拉解集{y,N}的二维B样条密度估计的概率密度分布图；

图14为一实施例提供的单一立方体模型的二维欧拉解集{z,N}的二维B样条密度估计的概率密度分布图；

图15为一实施例提供的立方体三维欧拉解集{x,y,z}的三维B样条密度估计的概率密度等值面图；

图16为一实施例提供的立方体三维欧拉解集{x,y,N}的三维B样条密度估计的概率密度等值面图；

图17为一实施例提供的立方体三维欧拉解集{x,z,N}的三维B样条密度估计的概率密度等值面图；

图18为一实施例提供的立方体三维欧拉解集{y,z,N}的三维B样条密度估计的概率密度等值面图；

图19为一实施例提供的立方体三维欧拉解集{x,y,z}基于快速傅立叶变换的三维B样条密度估计的概率密度等值面图；

图20为一实施例提供的立方体三维欧拉解集{x,y,N}基于快速傅立叶变换的三维B样条密度估计的概率密度等值面图；

图21为一实施例提供的立方体三维欧拉解集{x,z,N}基于快速傅立叶变换的三维B样条密度估计的概率密度等值面图；

图22为一实施例提供的立方体三维欧拉解集{y,z,N}基于快速傅立叶变换的三维B样条密度估计的概率密度等值面图。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，遂图式中仅显示与本发明中有关的组件而非按照实际实施时的组件数目、形状及尺寸绘制，其实际实施时各组件的型态、数量及比例可为一种随意的改变，且其组件布局型态也可能更为复杂。

请参阅图1至图22所示，本发明提供一种基于傅立叶变换的B样条函数欧拉解概率密度成像方法，包括以下步骤：

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

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

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

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

将欧拉解解集拆分成不同维度的组合，并根据所述欧拉解解集估算估计区间的上/下界及估计网格带宽，以及确定所述估计网格的大小；

构建基于B样条函数的多变量概率密度估计，并将欧拉解数据作为独立同分布采样投影至所述估计网格，计算网格计数；

对所述B样条函数和所述网格计数进行卷积，获取欧拉解概率密度分布结果，并根据所述欧拉解概率密度分布结果确定不同维度数据的异常源空间位置及类型。

本申请实施例采用基于规范化B样条的非参数概率密度估计方法，以各个欧拉解间的相似性及其集聚程度，分离各欧拉解簇以标示各异常源。针对B样条概率密度估计方法存在不聚焦及遍历估计网格计算量巨大等问题，本申请实施例充分利用估计网格同一维度方向具有相同的间隔之一特点，构建了基于B样条基函数的多变量概率估计方法，并借助快速傅里叶变换实现估计网格与B样条基函数的离散卷积，计算各欧拉解的概率密度，通过概率密度值的大小，从而达到快速定位异常源目的。相比于原有欧拉解剔除策略及B样条概率密度估计方法，本申请实施例的计算更为高效，且更易区分异常形态、判断异常位置。

根据上述记载，在一示例性实施例中，根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程，包括以三维重力欧拉方程为例确定三维欧拉反褶积方程，有：

式中，α＝x,y,z，B

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

确定滑动窗口大小为w＝w

利用公式(1)和滑动窗口构建矩形形式线性方程组：

式中，

n

欧拉反褶积解集为m＝{m

第i个欧拉反褶积解m

根据上述记载，构建基于B样条函数的多变量概率密度估计，并将欧拉解数据作为独立同分布采样投影至所述估计网格，计算网格计数的过程包括：

设定X

将样本X

式中系数

基于近似分箱技术的多变量概率密度估计方法重新构建公式(3)，有：

式中，

具体地，将公式(5)中的C

若X

其中，

设定公式(4)中的权重

将公式(5)写为卷积的形式，则有：

式中，

对矩阵C和r进行补零及环绕排序，使C和r具有相同的大小，则有：

其中，C和r具有相同的维度λ

根据离散卷积定理，令：

其中，F代表傅里叶正变换，F

待求的密度估计

式中，S[a:b，c:d]表示矩阵S从a到b的行的子集和从c到d的列的子集。

根据上述记载，本发明提供一种基于快速傅立叶变换的B样条函数重力/重力张量欧拉解概率密度成像方法的实施例，如图1所示，包括以下步骤：

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

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

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

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

将欧拉解解集拆分成不同维度的组合，并根据所述欧拉解解集组合估算估计区间的上/下界及估计网格带宽，以及确定所述估计网格的大小；

构建基于B样条函数的多变量概率密度估计，并将欧拉解解集组合作为独立同分布采样投影至所述估计网格，计算网格计数；

对所述B样条函数和所述网格计数进行卷积，获取欧拉解概率密度分布结果，并根据所述欧拉解概率密度分布结果确定不同维度数据的异常源空间位置及类型。具体地，

1)确定三维欧拉反褶积方程，以三维重力欧拉方程为例：

式中，α＝x,y,z，B

式中，T＝T(x,y,z)为位场函数；(x

2-1)构建欧拉线性方程组：利用一个大小为w＝w

式中，

3)欧拉解集的B样条密度估计，设X

式中系数

4)基于快速傅里叶变换的B样条密度估计，通过公式(3-1)可以计算各个欧拉解于估计网格x上的概率密度值，然而当数据样本过大或/和估计网格过大时，遍历估计网格上的每个节点是一个计算量且伴随内存消耗过大的之一问题，为此，借鉴基于近似分箱技术的多变量概率密度估计方法，重构式(3-1)：

式中，

为计算式(5-1)中C

其中，

则样本X

这里，向量Q为q通过二进制运算获得并翻转的向量，如以d＝3且q＝4为例，Q＝[0,0,1]。

根据Wand算法，需将公式(5-1)写为卷积的形式：

式中

式中，C和r的卷积可使用离散卷积定理计算，针对矩阵C和r大小不同的情况，根据特殊结构矩阵特性，对其补零及环绕排序，使C和r具有相同的大小，如矩阵(9-1)和(10-1)所示，为了便于理解，这里只给出二维的情况。

完成上述步骤后，C和r也就具有相同的维度λ

其中F代表傅里叶正变换，F

式中，对于二维情况，S[a:b，c:d]表示矩阵S从a到b的行的子集和从c到d的列的子集。

在一个具体实施例中，利用三个正态分布

建立理论模型，模型的立方体大小1000m×1000m×1000m，质心为(-1000，-2000，1500)，剩余密度均为0.36g/cm

确定三维欧拉方程，通过测量或位场转换获得重力矢量/重力梯度张量值；如图3和图4所示。

确定滑动窗口大小w＝w

式中，

将不同欧拉反褶积解集的数据组合，如{x}，{x,y}，{x,y,z}，{x,y,z,N}等，分别对不同数据组合进行估计，确定估计带宽h＝[h

借鉴基于近似分箱技术的多变量概率密度估计方法，重构式

利用分箱技术，按照式

利用快速傅立叶变换计算C

如图15至图22所示，相比于基于规范化B样条函数重力欧拉解概率密度成像结果，本发明获得概率密度成像结果更为聚焦。其中，图15至图18为基于规范化B样条函数重力欧拉解概率密度成像结果，图19至图22为本发明中基于傅立叶快速变换的B样条函数重力欧拉解概率密度成像结果。

本发明采用基于规范化B样条的非参数概率密度估计方法，以各个欧拉解间的相似性及其集聚程度，分离各欧拉解簇以标示各异常源。针对B样条概率密度估计方法存在不聚焦及遍历估计网格计算量巨大等问题，本发明充分利用估计网格同一维度方向具有相同的间隔之一特点，构建了基于B样条基函数的多变量概率估计方法，并借助快速傅里叶变换实现估计网格与B样条基函数的离散卷积，计算各欧拉解的概率密度，通过概率密度值的大小，从而达到快速定位异常源目的。相比于原有欧拉解剔除策略及B样条概率密度估计方法，本发明的计算更为高效，且更易区分异常形态、判断异常位置。且本发明还具有以下有益效果：

(1)本发明针对B样条概率密度估计方法，随着估计网格的增大存在不聚焦的问题，构建基于B样条概率密度的多变量密度估计方法；

(2)为避免当数据样本过大或估计网格过大时，遍历估计网格上的每个节点是一个计算量巨大且伴随内存消耗过大的之一问题，本发明引入分箱近似方法，将样本数据快速投影至估计网格，以借助快速傅里叶变换实现估计网格与密度函数的离散卷积；

(3)在无先验信息的情况下，本发明以各个欧拉解间的相似性及其集聚程度，作为计算欧拉解概率密度值的依据，进而区分由欧拉解组成的多个簇；

(4)与传统重力/重力梯度张量欧拉解解释技术难于剔除欧拉解缪解、确定欧拉解优解、分离欧拉解簇和图示欧拉反褶积解相比，本发明以各个欧拉解的概率密度值为基础，可通过各概率密度的峰值点，进而有效快速确定各异常源。

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