CN113096107A - 基于傅立叶变换的b样条函数欧拉解概率密度成像方法 - Google Patents

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

Info

Publication number
CN113096107A
CN113096107A CN202110424241.9A CN202110424241A CN113096107A CN 113096107 A CN113096107 A CN 113096107A CN 202110424241 A CN202110424241 A CN 202110424241A CN 113096107 A CN113096107 A CN 113096107A
Authority
CN
China
Prior art keywords
euler
probability density
estimation
solution
fourier transform
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.)
Granted
Application number
CN202110424241.9A
Other languages
English (en)
Other versions
CN113096107B (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.)
Hunan University of Science and Technology
Original Assignee
Hunan University of Science and Technology
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 Hunan University of Science and Technology filed Critical Hunan University of Science and Technology
Priority to CN202110424241.9A priority Critical patent/CN113096107B/zh
Publication of CN113096107A publication Critical patent/CN113096107A/zh
Application granted granted Critical
Publication of CN113096107B publication Critical patent/CN113096107B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20076—Probabilistic image processing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Complex Calculations (AREA)

Abstract

本发明提供一种基于傅立叶变换的B样条函数欧拉解概率密度成像方法，确定待测区域范围，测量重力矢量数据及重力梯度张量数据；构建三维重力梯度张量欧拉反褶积方程，确定滑动窗口大小，构建线性方程组，并获取欧拉解，输出欧拉解集；将欧拉解解集拆分成不同维度的组合，并估算相应组合的估计区间的上/下界及估计网格带宽，以及确定估计网格的大小；构建基于B样条函数的多变量概率密度估计，并将欧拉解解集组合作为独立同分布采样投影至估计网格，计算网格计数，获取欧拉解概率密度分布结果，确定不同维度数据的异常源空间位置及类型。本发明引入分箱近似方法，将样本数据投影至估计网格，基于快速傅里叶变换实现估计网格与密度函数的离散卷积。

Description

基于傅立叶变换的B样条函数欧拉解概率密度成像方法
技术领域
本发明涉及重力勘探技术领域，特别是涉及一种基于傅立叶变换的B样条函数欧拉解概率密度成像方法。
背景技术
欧拉反褶积是适合于大面积位场数据的一种半自动/自动估算场源位置的位场数据解释方法。源于欧拉超定方程组条件数大，致使欧拉解扰动大，且往往具有发散的趋势。因而，很多学者以欧拉超定方程组的标准差或截断误差作为标准，以过滤欧拉解集中的谬解；一些学者增加附加约束方程或约束条件，如利用重力梯度张量多分量含丰富异常信息(Euler deconvolution of gravity tensor gradient data)及利用“Worming”与欧拉方程之间的关系，构建混合欧拉反褶积方法以获取优解；或通过确定欧拉解集所构成的簇于欧拉解集关系，以剔除不隶属于各簇的欧拉解。然而对于相邻较近的异常源而言，传统簇分析方法难于区分各欧拉解簇。
发明内容
鉴于以上所述现有技术的缺点，本发明的目的在于提供一种基于快速傅立叶变换的B样条概率密度估计方法对欧拉解集进行处理，通过引入分箱近似方法，将样本数据快速投影至估计网格，以借助快速傅里叶变换实现估计网格与密度函数的离散卷积，在避免当数据样本过大或/和估计网格过大时，遍历估计网格上的每个节点是一个计算量巨大且伴随内存消耗过大的同时，获得较为聚焦的概率密度估计结果，以标示各异常源。
为实现上述目的及其他相关目的，本发明提供一种基于傅立叶变换的B样条函数欧拉解概率密度成像方法，包括有：
确定待测区域范围；测量重力矢量数据及重力梯度张量数据；或测量重力数据，通过离散余弦变换或傅里叶变换将重力数据换算到重力矢量数据及重力梯度张量数据；
根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程；
确定滑动窗口大小，并利用滑动窗口内的数据构建线性方程组；
利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解，并输出对应的欧拉解集；其中，所述欧拉解包括异常源空间位置、异常源类型和奇异值分解误差；
将欧拉解解集拆分成不同维度的组合，并根据所述欧拉解解集估算估计区间的上/下界及估计网格带宽，以及确定所述估计网格的大小；
构建基于B样条函数的多变量概率密度估计，并将欧拉解数据作为独立同分布采样投影至所述估计网格，计算网格计数；
对所述B样条函数和所述网格计数进行卷积，获取欧拉解概率密度分布结果，并根据所述欧拉解概率密度分布结果确定不同维度数据的异常源空间位置及类型。
可选地，所述根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程，包括以三维重力欧拉方程为例确定三维欧拉反褶积方程；
Figure BDA0003029211750000021
式中，α＝x,y,z，Bα为异常背景场，gα为重力梯度，gαβ为重力梯度张量分量；(x0,y0,z0)为观测点位置坐标；(x,y,z)为待求异常源的空间位置；构造指数N为异常幅值随距离变化的衰减率。
可选地，所述确定滑动窗口大小，并利用滑动窗口内的数据构建线性方程组，以及利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解，并输出对应的欧拉解集，包括：
确定滑动窗口大小为w＝wx＝wy；
利用公式(1)和滑动窗口构建矩形形式线性方程组：
Figure BDA0003029211750000031
式中，
Figure BDA0003029211750000032
nw＝wx×wy，α＝x,y,z；
欧拉反褶积解集为m＝{mi}；
第i个欧拉反褶积解mi＝[xi,yi,zi,Ni]T，且第i个欧拉反褶积解即为待求异常源空间位置及构造指数N。
可选地，所述构建基于B样条函数的多变量概率密度估计，并将欧拉解数据作为独立同分布采样投影至所述估计网格，计算网格计数的过程包括：
设定X1,...,Xi,...,Xn是d维概率密度f函数的独立同分布随机样本，且第i个样本可以表示为
Figure BDA0003029211750000033
其中，Xi在本文中可以是欧拉解xi,yi,zi,Ni的不同组合，包括：xi，[xi,yi]T，[xi,yi,zi]T，[xi,yi,zi,Ni]T等；
将样本Xi投影至大小为λ＝[λ1,...,λj,...,λd]T且带宽为h＝[h1,...,hj,...,hd]T的估计网格x＝[x1,...,xj,...,xd]T，则可以利用c次均匀B样条函数
Figure BDA0003029211750000034
表达d维B样条概率密度估计
Figure BDA0003029211750000035
即：
Figure BDA0003029211750000036
式中系数
Figure BDA0003029211750000037
用权重u表示为：
Figure BDA0003029211750000038
基于近似分箱技术的多变量概率密度估计方法重新构建公式(3)，有：
Figure BDA0003029211750000039
式中，
Figure BDA00030292117500000310
Cl为近似分箱方法计算得到的由估计网格xl上的网格计数。
可选地，将公式(5)中的Cl定义为在j维度方向存在边
Figure BDA0003029211750000041
并简写为
Figure BDA0003029211750000042
若Xi恰位于由边
Figure BDA0003029211750000043
构成的且含2d个节点和2d个面的超矩形中，则利用
Figure BDA0003029211750000044
在
Figure BDA0003029211750000045
上的投影表示边
Figure BDA0003029211750000046
两端的权重
Figure BDA0003029211750000047
和
Figure BDA0003029211750000048
有：
Figure BDA0003029211750000049
其中，
Figure BDA00030292117500000410
为向下取整运算符，两个端点的序号分别表示为j-1和2d-1+j-1。
可选地，设定公式(4)中的权重
Figure BDA00030292117500000411
为
Figure BDA00030292117500000412
于某超定矩阵在j维度方向的投影，则样本Xi于超矩形第q个节点的网格近似Cl可以写为：
Figure BDA00030292117500000413
其中，向量Q为q通过二进制运算获得并翻转的向量。
可选地，将公式(5)写为卷积的形式，则有：
Figure BDA00030292117500000414
式中，
Figure BDA00030292117500000415
对矩阵C和r进行补零及环绕排序，使C和r具有相同的大小，则有：
Figure BDA00030292117500000416
Figure BDA0003029211750000051
其中，C和r具有相同的维度λ1×,...,×λj×,...,×λd。
可选地，根据离散卷积定理，令：
Figure BDA0003029211750000052
其中，F代表傅里叶正变换，F-1为傅里叶逆变换；
待求的密度估计
Figure BDA0003029211750000053
可由方程(11)中S的子集除以λ1×...×λj×...×λd的乘积得到，有：
Figure BDA0003029211750000054
式中，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样条概率密度估计方法，本申请实施例的计算更为高效，且更易区分异常形态、判断异常位置。
根据上述记载，在一示例性实施例中，根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程，包括以三维重力欧拉方程为例确定三维欧拉反褶积方程，有：
Figure BDA0003029211750000091
式中，α＝x,y,z，Bα为异常背景场，gα为重力梯度，gαβ为重力梯度张量分量；(x0,y0,z0)为观测点位置坐标；(x,y,z)为待求异常源的空间位置；构造指数N为异常幅值随距离变化的衰减率。
根据上述记载，确定滑动窗口大小，并利用滑动窗口内的数据构建线性方程组，以及利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解，并输出对应的欧拉解集的过程包括：
确定滑动窗口大小为w＝wx＝wy；
利用公式(1)和滑动窗口构建矩形形式线性方程组：
Figure BDA0003029211750000092
式中，
Figure BDA0003029211750000093
nw＝wx×wy，α＝x,y,z；
欧拉反褶积解集为m＝{mi}；
第i个欧拉反褶积解mi＝[xi,yi,zi,Ni]T，且第i个欧拉反褶积解即为待求异常源空间位置及构造指数N。
根据上述记载，构建基于B样条函数的多变量概率密度估计，并将欧拉解数据作为独立同分布采样投影至所述估计网格，计算网格计数的过程包括：
设定X1,...,Xi,...,Xn是d维概率密度f函数的独立同分布随机样本，且第i个样本可以表示为
Figure BDA0003029211750000094
其中，Xi在本文中可以是欧拉解xi,yi,zi,Ni的不同组合，包括：xi，[xi,yi]T，[xi,yi,zi]T，[xi,yi,zi,Ni]T；
将样本Xi投影至大小为λ＝[λ1,...,λj,...,λd]T且带宽为h＝[h1,...,hj,...,hd]T的估计网格x＝[x1,...,xj,...,xd]T，则可以利用c次均匀B样条函数
Figure BDA0003029211750000101
表达d维B样条概率密度估计
Figure BDA0003029211750000102
即：
Figure BDA0003029211750000103
式中系数
Figure BDA0003029211750000104
用权重u表示为：
Figure BDA0003029211750000105
基于近似分箱技术的多变量概率密度估计方法重新构建公式(3)，有：
Figure BDA0003029211750000106
式中，
Figure BDA0003029211750000107
Cl为近似分箱方法计算得到的由估计网格xl上的网格计数。
具体地，将公式(5)中的Cl定义为在j维度方向存在边
Figure BDA0003029211750000108
并简写为
Figure BDA0003029211750000109
若Xi恰位于由边
Figure BDA00030292117500001010
构成的且含2d个节点和2d个面的超矩形中，则利用
Figure BDA00030292117500001011
在
Figure BDA00030292117500001012
上的投影表示边
Figure BDA00030292117500001013
两端的权重
Figure BDA00030292117500001014
和
Figure BDA00030292117500001015
有：
Figure BDA00030292117500001016
其中，
Figure BDA00030292117500001017
为向下取整运算符，两个端点的序号分别表示为j-1和2d-1+j-1。
设定公式(4)中的权重
Figure BDA00030292117500001018
为
Figure BDA00030292117500001019
于某超定矩阵在j维度方向的投影，则样本Xi于超矩形第q个节点的网格近似Cl可以写为：
Figure BDA00030292117500001020
其中，向量Q为q通过二进制运算获得并翻转的向量。
将公式(5)写为卷积的形式，则有：
Figure BDA00030292117500001021
式中，
Figure BDA0003029211750000111
对矩阵C和r进行补零及环绕排序，使C和r具有相同的大小，则有：
Figure BDA0003029211750000112
Figure BDA0003029211750000113
其中，C和r具有相同的维度λ1×,...,×λj×,...,×λd。
根据离散卷积定理，令：
Figure BDA0003029211750000114
其中，F代表傅里叶正变换，F-1为傅里叶逆变换；
待求的密度估计
Figure BDA0003029211750000115
可由方程(11)中S的子集除以λ1×...×λj×...×λd的乘积得到，有：
Figure BDA0003029211750000116
式中，S[a:b，c:d]表示矩阵S从a到b的行的子集和从c到d的列的子集。
根据上述记载，本发明提供一种基于快速傅立叶变换的B样条函数重力/重力张量欧拉解概率密度成像方法的实施例，如图1所示，包括以下步骤：
确定待测区域范围；测量重力矢量数据及重力梯度张量数据；或测量重力数据，通过离散余弦变换或傅里叶变换将重力数据换算到重力矢量数据及重力梯度张量数据；
根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程；
确定滑动窗口大小，并利用滑动窗口内的数据构建线性方程组；
利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解，并输出对应的欧拉解集；其中，所述欧拉解包括异常源空间位置、异常源类型和奇异值分解误差；
将欧拉解解集拆分成不同维度的组合，并根据所述欧拉解解集组合估算估计区间的上/下界及估计网格带宽，以及确定所述估计网格的大小；
构建基于B样条函数的多变量概率密度估计，并将欧拉解解集组合作为独立同分布采样投影至所述估计网格，计算网格计数；
对所述B样条函数和所述网格计数进行卷积，获取欧拉解概率密度分布结果，并根据所述欧拉解概率密度分布结果确定不同维度数据的异常源空间位置及类型。具体地，
1)确定三维欧拉反褶积方程，以三维重力欧拉方程为例：
Figure BDA0003029211750000121
式中，α＝x,y,z，Bα为异常背景场，gα为重力梯度，gαβ为重力梯度张量分量。
式中，T＝T(x,y,z)为位场函数；(x0,y0,z0)为观测点位置坐标；(x,y,z)为待求异常源的空间位置；N为异常幅值随距离变化的衰减率，特定的地质体对应特定的衰减率，即构造指数。
2-1)构建欧拉线性方程组：利用一个大小为w＝wx＝wy的正方形滑动窗口逐步覆盖整个测区，每个滑动窗口可以由方程(1-1)组成的如下矩阵形式方程组：
Figure BDA0003029211750000122
式中，
Figure BDA0003029211750000131
nw＝wx×wy，α＝x,y,z；欧拉反褶积解集为m＝{mi}，第i个欧拉反褶积解mi＝[xi,yi,zi,Ni]T，即为待求异常源空间位置及构造指数。
3)欧拉解集的B样条密度估计，设X1,...,Xi,...,Xn是d维概率密度f函数的独立同分布随机样本，第i个样本可以表示为
Figure BDA0003029211750000132
特别地，Xi在本文中可以是欧拉解xi,yi,zi,Ni的不同组合，如xi，[xi,yi]T，[xi,yi,zi]T，[xi,yi,zi,Ni]T等。将样本Xi投影至大小为λ＝[λ1,...,λj,...,λd]T且带宽为h＝[h1,...,hj,...,hd]T的估计网格x＝[x1,...,xj,...,xd]T，则可以利用c次均匀B样条函数
Figure BDA0003029211750000133
表达d维B样条概率密度估计
Figure BDA0003029211750000134
即：
Figure BDA0003029211750000135
式中系数
Figure BDA0003029211750000136
可以用权重u表示为
Figure BDA0003029211750000137
4)基于快速傅里叶变换的B样条密度估计，通过公式(3-1)可以计算各个欧拉解于估计网格x上的概率密度值，然而当数据样本过大或/和估计网格过大时，遍历估计网格上的每个节点是一个计算量且伴随内存消耗过大的之一问题，为此，借鉴基于近似分箱技术的多变量概率密度估计方法，重构式(3-1)：
Figure BDA0003029211750000138
式中，
Figure BDA0003029211750000139
Cl为近似分箱方法计算得到的由估计网格xl上的网格计数。
为计算式(5-1)中Cl，将其定义在j维度方向存在边
Figure BDA00030292117500001310
并简写为
Figure BDA00030292117500001311
假定Xi恰位于由边
Figure BDA00030292117500001312
构成的且含2d个节点和2d个面的超矩形中，则可以利用
Figure BDA00030292117500001313
在
Figure BDA00030292117500001314
上的投影表示边
Figure BDA00030292117500001315
两端的权中
Figure BDA00030292117500001316
和
Figure BDA00030292117500001317
即：
Figure BDA0003029211750000141
其中，
Figure BDA0003029211750000142
为向下取整运算符。该边两个端点的序号分别表示为j-1和2d-1+j-1。令式(4-1)中权重
Figure BDA0003029211750000143
为
Figure BDA0003029211750000144
于某超定矩阵在j维度方向的投影。
则样本Xi于超矩形第q个节点的网格近似Cl可以写为
Figure BDA0003029211750000145
这里，向量Q为q通过二进制运算获得并翻转的向量，如以d＝3且q＝4为例，Q＝[0,0,1]。
根据Wand算法，需将公式(5-1)写为卷积的形式：
Figure BDA0003029211750000146
式中
Figure BDA0003029211750000147
式中，C和r的卷积可使用离散卷积定理计算，针对矩阵C和r大小不同的情况，根据特殊结构矩阵特性，对其补零及环绕排序，使C和r具有相同的大小，如矩阵(9-1)和(10-1)所示，为了便于理解，这里只给出二维的情况。
Figure BDA0003029211750000148
Figure BDA0003029211750000149
完成上述步骤后，C和r也就具有相同的维度λ1×,...,×λj×,...,×λd，根据离散卷积定理，令
Figure BDA0003029211750000151
其中F代表傅里叶正变换，F-1为傅里叶逆变换。则待求的密度估计
Figure BDA0003029211750000152
可由方程(11-1)中S的子集除以λ1×...×λj×...×λd的乘积得到
Figure BDA0003029211750000153
式中，对于二维情况，S[a:b，c:d]表示矩阵S从a到b的行的子集和从c到d的列的子集。
在一个具体实施例中，利用三个正态分布
Figure BDA0003029211750000154
构建一组样本数为n＝3000的一维随机数据，并引入高斯核平滑估计算法作为对比，以验证B样条密度估计算法的可靠性，如图2所示。在图2中，true为真值，BSS为B样条概率密度估计方法，KS为高斯核平滑估计方法，BSSFFT为基于快速傅立叶变换的B样条概率密度估计方法。
建立理论模型，模型的立方体大小1000m×1000m×1000m，质心为(-1000，-2000，1500)，剩余密度均为0.36g/cm3。观测高度为50m，测区范围为x：-10000m～10000m，y：-10000m～10000m；测网大小为100m×100m。由理论模型产生的地面重力数据，如图3及含3％噪声数据的图4所示，共有200×200＝40000个采样点。本发明针对这个模型实例进行详细说明，具体地，
确定三维欧拉方程，通过测量或位场转换获得重力矢量/重力梯度张量值；如图3和图4所示。
确定滑动窗口大小w＝wx＝wy，利用该正方形滑动窗口逐步覆盖整个测区，每个滑动窗口组成如下矩阵形式方程组：
Figure BDA0003029211750000155
式中，
Figure BDA0003029211750000156
nw＝wx×wy，α＝x,y,z；欧拉反褶积解集为m＝{mi}，第i个欧拉反褶积解mi＝[xi,yi,zi,Ni]T，即为待求异常源空间位置及构造指数，如图5和图6所示。
将不同欧拉反褶积解集的数据组合，如{x}，{x,y}，{x,y,z}，{x,y,z,N}等，分别对不同数据组合进行估计，确定估计带宽h＝[h1,...,hj,...,hd]T，估计网格大小λ＝[λ1,...,λj,...,λd]T。
借鉴基于近似分箱技术的多变量概率密度估计方法，重构式
Figure BDA0003029211750000161
得到：
Figure BDA0003029211750000162
式中，
Figure BDA0003029211750000163
Cl为近似分箱方法计算得到的由估计网格xl上的网格计数。
利用分箱技术，按照式
Figure BDA0003029211750000164
将观测样本某一维度信息{Xj}直接投影至估计网格xj，计算各观测数据样本于各估计网格点上的权重，并计算各观测数据样本与估计网格上的网格计数Cl。
利用快速傅立叶变换计算Cl与B样条基函数间的卷积，获得欧拉解概率密度结果，以概率密度曲线图(即图7和图8)、概率密度分布图(及图9至图14)或概率密度等值面(即图19至图22)显示，并对这些结果进行解释，剔除欧拉解集中的缪解，分离欧拉解所构成的簇，以标识相邻异常源。
如图15至图22所示，相比于基于规范化B样条函数重力欧拉解概率密度成像结果，本发明获得概率密度成像结果更为聚焦。其中，图15至图18为基于规范化B样条函数重力欧拉解概率密度成像结果，图19至图22为本发明中基于傅立叶快速变换的B样条函数重力欧拉解概率密度成像结果。
本发明采用基于规范化B样条的非参数概率密度估计方法，以各个欧拉解间的相似性及其集聚程度，分离各欧拉解簇以标示各异常源。针对B样条概率密度估计方法存在不聚焦及遍历估计网格计算量巨大等问题，本发明充分利用估计网格同一维度方向具有相同的间隔之一特点，构建了基于B样条基函数的多变量概率估计方法，并借助快速傅里叶变换实现估计网格与B样条基函数的离散卷积，计算各欧拉解的概率密度，通过概率密度值的大小，从而达到快速定位异常源目的。相比于原有欧拉解剔除策略及B样条概率密度估计方法，本发明的计算更为高效，且更易区分异常形态、判断异常位置。且本发明还具有以下有益效果：
(1)本发明针对B样条概率密度估计方法，随着估计网格的增大存在不聚焦的问题，构建基于B样条概率密度的多变量密度估计方法；
(2)为避免当数据样本过大或估计网格过大时，遍历估计网格上的每个节点是一个计算量巨大且伴随内存消耗过大的之一问题，本发明引入分箱近似方法，将样本数据快速投影至估计网格，以借助快速傅里叶变换实现估计网格与密度函数的离散卷积；
(3)在无先验信息的情况下，本发明以各个欧拉解间的相似性及其集聚程度，作为计算欧拉解概率密度值的依据，进而区分由欧拉解组成的多个簇；
(4)与传统重力/重力梯度张量欧拉解解释技术难于剔除欧拉解缪解、确定欧拉解优解、分离欧拉解簇和图示欧拉反褶积解相比，本发明以各个欧拉解的概率密度值为基础，可通过各概率密度的峰值点，进而有效快速确定各异常源。
上述实施例仅例示性说明本发明的原理及其功效，而非用于限制本发明。任何熟悉此技术的人士皆可在不违背本发明的精神及范畴下，对上述实施例进行修饰或改变。因此，举凡所属技术领域中具有通常知识者在未脱离本发明所揭示的精神与技术思想下所完成的一切等效修饰或改变，仍应由本发明的权利要求所涵盖。

Claims (8)

1.一种基于傅立叶变换的B样条函数欧拉解概率密度成像方法，其特征在于，包括以下步骤：
确定待测区域范围；测量重力矢量数据及重力梯度张量数据；或测量重力数据，通过离散余弦变换或傅里叶变换将重力数据换算到重力矢量数据及重力梯度张量数据；
根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程；
确定滑动窗口大小，并利用滑动窗口内的数据构建线性方程组；
利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解，并输出对应的欧拉解集；其中，所述欧拉解包括异常源空间位置、异常源类型和奇异值分解误差；
将欧拉解解集拆分成不同维度的组合，并根据所述欧拉解解集组合估算估计区间的上/下界及估计网格带宽，以及确定所述估计网格的大小；
构建基于B样条函数的多变量概率密度估计，并将欧拉解解集组合作为独立同分布采样投影至所述估计网格，计算网格计数；
对所述B样条函数和所述网格计数进行卷积，获取欧拉解概率密度分布结果，并根据所述欧拉解概率密度分布结果确定不同维度数据的异常源空间位置及类型。
2.根据权利要求1所述的基于傅立叶变换的B样条函数欧拉解概率密度成像方法，其特征在于，所述根据所述重力矢量数据或所述重力梯度张量数据构建三维重力梯度张量欧拉反褶积方程，包括以三维重力欧拉方程为例确定三维欧拉反褶积方程；
Figure FDA0003029211740000011
式中，α＝x,y,z，Bα为异常背景场，gα为重力梯度，gαβ为重力梯度张量分量；(x0,y0,z0)为观测点位置坐标；(x,y,z)为待求异常源的空间位置；构造指数N为异常幅值随距离变化的衰减率。
3.根据权利要求2所述的基于傅立叶变换的B样条函数欧拉解概率密度成像方法，其特征在于，所述确定滑动窗口大小，并利用滑动窗口内的数据构建线性方程组，以及利用奇异值分解所述线性方程组，获取所述线性方程组的欧拉解，并输出对应的欧拉解集，包括：
确定滑动窗口大小为w＝wx＝wy；
利用公式(1)和滑动窗口构建矩形形式线性方程组：
Figure FDA0003029211740000021
式中，
Figure FDA0003029211740000022
nw＝wx×wy，α＝x,y,z；
欧拉反褶积解集为m＝{mi}；
第i个欧拉反褶积解mi＝[xi,yi,zi,Ni]T，且第i个欧拉反褶积解即为待求异常源空间位置及构造指数N。
4.根据权利要求3所述的基于傅立叶变换的B样条函数欧拉解概率密度成像方法，其特征在于，所述构建基于B样条函数的多变量概率密度估计，并将欧拉解解集组合作为独立同分布采样投影至所述估计网格，计算网格计数的过程包括：
设定X1,...,Xi,...,Xn是d维概率密度f函数的独立同分布随机样本，且第i个样本可以表示为
Figure FDA0003029211740000023
其中，Xi在本文中可以是欧拉解xi,yi,zi,Ni的不同组合，包括：xi，[xi,yi]T，[xi,yi,zi]T，[xi,yi,zi,Ni]T等；
将样本Xi投影至大小为λ＝[λ1,...,λj,...,λd]T且带宽为h＝[h1,...,hj,...,hd]T的估计网格x＝[x1,...,xj,...,xd]T，则可以利用c次均匀B样条函数
Figure FDA0003029211740000024
表达d维B样条概率密度估计
Figure FDA0003029211740000031
即：
Figure FDA0003029211740000032
式中系数
Figure FDA0003029211740000033
用权重u表示为：
Figure FDA0003029211740000034
基于近似分箱技术的多变量概率密度估计方法重新构建公式(3)，有：
Figure FDA0003029211740000035
式中，
Figure FDA0003029211740000036
Cl为近似分箱方法计算得到的由估计网格xl上的网格计数。
5.根据权利要求4所述的基于傅立叶变换的B样条函数欧拉解概率密度成像方法，其特征在于，将公式(5)中的Cl定义为在j维度方向存在边
Figure FDA0003029211740000037
并简写为
Figure FDA0003029211740000038
若Xi恰位于由边
Figure FDA0003029211740000039
构成的且含2d个节点和2d个面的超矩形中，则利用
Figure FDA00030292117400000310
在
Figure FDA00030292117400000311
上的投影表示边
Figure FDA00030292117400000312
两端的权重
Figure FDA00030292117400000313
和
Figure FDA00030292117400000314
有：
Figure FDA00030292117400000315
其中，
Figure FDA00030292117400000316
为向下取整运算符，两个端点的序号分别表示为j-1和2d-1+j-1。
6.根据权利要求5所述的基于傅立叶变换的B样条函数欧拉解概率密度成像方法，其特征在于，设定公式(4)中的权重
Figure FDA00030292117400000317
为
Figure FDA00030292117400000318
于某超定矩阵在j维度方向的投影，则样本Xi于超矩形第q个节点的网格近似Cl可以写为：
Figure FDA00030292117400000319
其中，向量Q为q通过二进制运算获得并翻转的向量。
7.根据权利要求6所述的基于傅立叶变换的B样条函数欧拉解概率密度成像方法，其特征在于，将公式(5)写为卷积的形式，则有：
Figure FDA00030292117400000320
式中，
Figure FDA0003029211740000041
对矩阵C和r进行补零及环绕排序，使C和r具有相同的大小，则有：
Figure FDA0003029211740000042
Figure FDA0003029211740000043
其中，C和r具有相同的维度λ1×,...,×λj×,...,×λd。
8.根据权利要求7所述的基于傅立叶变换的B样条函数欧拉解概率密度成像方法，其特征在于，根据离散卷积定理，令：
Figure FDA0003029211740000044
其中，F代表傅里叶正变换，F-1为傅里叶逆变换；
待求的密度估计
Figure FDA0003029211740000045
可由方程(11)中S的子集除以λ1×...×λj×...×λd的乘积得到，有：
Figure FDA0003029211740000046
式中，S[a:b，c:d]表示矩阵S从a到b的行的子集和从c到d的列的子集。
CN202110424241.9A 2021-04-20 2021-04-20 基于傅立叶变换的b样条函数欧拉解概率密度成像方法 Active CN113096107B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110424241.9A CN113096107B (zh) 2021-04-20 2021-04-20 基于傅立叶变换的b样条函数欧拉解概率密度成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110424241.9A CN113096107B (zh) 2021-04-20 2021-04-20 基于傅立叶变换的b样条函数欧拉解概率密度成像方法

Publications (2)

Publication Number Publication Date
CN113096107A true CN113096107A (zh) 2021-07-09
CN113096107B CN113096107B (zh) 2024-08-02

Family

ID=76678900

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110424241.9A Active CN113096107B (zh) 2021-04-20 2021-04-20 基于傅立叶变换的b样条函数欧拉解概率密度成像方法

Country Status (1)

Country Link
CN (1) CN113096107B (zh)

Citations (6)

* 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차원 영상화 방법
WO2015154601A1 (zh) * 2014-04-08 2015-10-15 中山大学 一种基于无特征提取的紧致sfm三维重建方法
CN105203104A (zh) * 2015-09-16 2015-12-30 北京航空航天大学 一种适用于高精度惯导系统的重力场建模方法
CN108732622A (zh) * 2018-05-18 2018-11-02 吉林大学 一种不同高度数据融合联合反演地质体几何形态的方法
CN112462442A (zh) * 2020-11-30 2021-03-09 山东大学 重磁位场场源位置估计方法、系统、介质及电子设备
CN112668146A (zh) * 2020-12-03 2021-04-16 重庆科技学院 一种基于欧拉反褶积法实用性改进的场源位置估算方法

Patent Citations (6)

* 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차원 영상화 방법
WO2015154601A1 (zh) * 2014-04-08 2015-10-15 中山大学 一种基于无特征提取的紧致sfm三维重建方法
CN105203104A (zh) * 2015-09-16 2015-12-30 北京航空航天大学 一种适用于高精度惯导系统的重力场建模方法
CN108732622A (zh) * 2018-05-18 2018-11-02 吉林大学 一种不同高度数据融合联合反演地质体几何形态的方法
CN112462442A (zh) * 2020-11-30 2021-03-09 山东大学 重磁位场场源位置估计方法、系统、介质及电子设备
CN112668146A (zh) * 2020-12-03 2021-04-16 重庆科技学院 一种基于欧拉反褶积法实用性改进的场源位置估算方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
曹书锦;朱自强;鲁光银: "基于自适应模糊聚类分析的重力张量欧拉反褶积解", 中南大学学报. 自然科学版, vol. 43, no. 3, 31 December 2012 (2012-12-31) *
朱自强;王灿;鲁光银;曹书锦;: "重力梯度张量解析信号的欧拉反褶积", 中南大学学报(自然科学版), no. 01, 26 January 2015 (2015-01-26) *

Also Published As

Publication number Publication date
CN113096107B (zh) 2024-08-02

Similar Documents

Publication Publication Date Title
Park et al. Discrete sibson interpolation
Lawson et al. Spatial cluster modelling
Sarkar et al. An efficient approach to estimate fractal dimension of textural images
Bachthaler et al. Continuous scatterplots
Du et al. Meshfree, probabilistic determination of point sets and support regions for meshless computing
Smolik et al. Large scattered data interpolation with radial basis functions and space subdivision
Davis et al. Efficient 3D inversion of magnetic data via octree-mesh discretization, space-filling curves, and wavelets
Gesemann From particle tracks to velocity and acceleration fields using B-splines and penalties
CN101017476A (zh) 基于数据动态存取模型的产品点云型面特征分析方法
CN104616349A (zh) 基于局部曲面变化因子的散乱点云数据精简处理方法
Wang et al. A new point cloud simplification method with feature and integrity preservation by partition strategy
CN111738278B (zh) 一种水下多源声学图像特征提取方法及系统
Peters et al. Robust approximation of the Medial Axis Transform of LiDAR point clouds as a tool for visualisation
Mosegaard et al. Inverse methods: Problem formulation and probabilistic solutions
Karbauskaitė et al. Fractal-based methods as a technique for estimating the intrinsic dimensionality of high-dimensional data: a survey
CN117932974B (zh) 一种水库水下数字高程模型的构建方法
Kirveslahti et al. Representing fields without correspondences: the lifted euler characteristic transform
CN105719349B (zh) 基于最大化泊松圆盘采样的四面体网格化方法和系统
CN107748834B (zh) 一种计算起伏观测面磁场的快速、高精度数值模拟方法
Jjumba et al. Spatial indices for measuring three-dimensional patterns in a voxel-based space
CN113096107B (zh) 基于傅立叶变换的b样条函数欧拉解概率密度成像方法
Chang et al. Reverse engineering of a symmetric object
Kamatsuchi Turbulent flow simulation around complex geometries with cartesian grid method
Coles et al. Quantifying the topology of large-scale structure
CN113138426B (zh) 一种基于多变量核密度欧拉解概率密度成像方法

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