一种基于物理空间树状结构的跨维贝叶斯采样器的设计方法

本发明涉及地球物理反演，尤其涉及一种基于物理空间树状结构的跨维贝叶斯采样器的设计方法。

背景技术：

1、地球物理大地电磁法勘探是一种利用天然电磁场源获取地下物理响应的地球物理方法，广泛应用于油气、矿产和地下深层地电结构等领域的研究。由于实际观测中存在干扰，反演结果也表现出非唯一性。因此，随着勘探的深入和计算机软硬件发展水平的快速提高，大地电磁反演结果的可靠性评价逐渐成为国内外学者研究的重点和热点。对反演结果进行不确定性分析，给出地下电性可能存在的分布情况，将极大的提高人们对反演结果量化评价，并在实际工程中给出指导意见。

2、现阶段，为了体现反演结果的不确定性并量化评价反演结果，都是进行基于贝叶斯框架下的概率反演。利用马尔科夫蒙特卡洛（markov chain montecarlo，mcmc）采样方法进行模型提取，并进行参数的统计推断。

3、在贝叶斯反演中，由于实际勘探对地下构造的了解是有限的，无法提前确定模型参数的维度。文献(green, p.j., 1995. reversible jump markov chain monte carlocomputation and bayesian model determination, biometrika,82(4), 711-732.) 将模型参数视为未知量，并将其纳入贝叶斯框架中，这种方法被称为跨维贝叶斯反演。跨维贝叶斯反演使得模型的复杂性由数据驱动，无需人为选择优化参数，从而避免了人为因素对反演结果的影响，该方法成为一种主流方法。文献(malinverno, a., 2002. parsimoniousbayesian markov chain monte carlo inversion in a nonlinear geophysicalproblem, geophys. j. int., 151(3),675–688.)首次将其直流电阻反演中，该方法给出了一种贝叶斯参数估计的自然扩展，即将地球模型参数化为层状介质，将层数作为参数进行反演，并从后验分布中获取分层介质的样本。借助这种自动参数化方法，可以使得在数据覆盖更好、噪声低的地区提供更高的模型分辨率。然而，该参数化方法无法处理高维的大地电磁反演问题。

4、目前，在高维参数空间中进行参数跨维采样的方法中，voronoi单元法是一种常用的方法。文献(bodin, t., sambridge, m., rawlinson, n. & arroucau, p., 2012.transdimensional tomography with unknown data noise, geophys. j. int., 189,1536–1556.)通过连接相邻原子核的垂直平分线构成一组连续多边形，从而将地下电导率剖面离散成一系列voronoi多边形。通过增加或减少原子核的数量来实现对模型参数数目的跨维采样。该方法的优点在于可以通过指定原子核的位置和每个细胞内部的属性值来表示地球物理模型空间的参数位置和物性参数，从而实现对模型的灵活划分。然而，voronoi单元法需要在一个单元中包含三个未知量。随着地球物理模型数据量的增加，其复杂度也会大大增加，其时间成本也会大大增加。此外，voronoi单元剖分方法的网格剖分是不均匀的任意剖分，这对于基于有限差分的大地电磁正演计算方法带来了挑战。文献(hawkins,r. & sambridge, m., 2015. geophysical imaging using transdimensional trees,geophys. j. int, 203, 972–1000.) 提出了一种小波树状结构的采样方法。该方法将地球物理模型在物理空间中的反演转化为在基于小波基函数的树状结构中进行研究和操作，可以用少量的参数表示复杂的地球物理模型。但该方法由于树状结构节点保存的是小波系数，而不代表任何物理意义，因此在选择先验分布时面临困难，无法与物理空间中的电导率范围相对应。同时，不同小波基函数的选择会直接影响最终的反演结果，因此该方法的应用条件也有一定的限制。因此，有必要提供一种新的树状结构来解决上述技术问题。

技术实现思路

1、本发明的主要目的是提供一种基于物理空间树状结构的跨维贝叶斯采样器的设计方法，旨在解决小波域树状结构的节点对先验分布选择困难，采样效率低的问题。

2、为实现上述目的，本发明提出的一种基于物理空间树状结构的跨维贝叶斯采样器的设计方法，所述基于物理空间树状结构的跨维贝叶斯采样器的设计方法包括：

3、s1：采用电磁测量仪器测量获得某区域的地下电磁场响应，并确定测量的地下深度；

4、s2：构建树状结构初始模型，构建马尔科夫链，具体是：

5、计算树状结构的最大深度，公式为：，其中，表示树状结构最大深度下的地下深度的剖分数量，表示树状结构的最大深度，并构建最大深度下的二叉树结构；

6、构建采样次数时的树状结构初始模型 ，其中，为初始模型中的节点数，为树状结构的深度等级，， 表示为深度下的个节点数所构成全部不同形状的树状结构的排列数，表示第个节点所构成的电导率参数；

7、从时开始构建马尔科夫链，对第次采样时的现有树状结构而言，可以将树状结构的节点分类为现有节点、可生成节点以及可消亡节点；

8、s3：任选树状结构的一类节点进行树状结构变化，形成新候选模型；若选择可生成节点，则表示新候选树状结构节点，新候选模型为；若选择可消亡节点，则表示新候选树状结构节点，新候选模型为；若选择现有节点中的一个时表示仅仅改变选中节点的电导率参数，新候选模型为；

9、s4：将s3中得到的树状结构下的新候选模型进行逆变换得到电导率结构的建议模型，并通过正演算子将建议模型映射到数据空间得到正演响应；

10、s5：基于s1中的电磁场响应和s4中得到的正演响应计算建议模型的似然函数；

11、s6：基于s5中的似然函数计算得到建议模型的接受率，取0-1之间的随机数，判断接受率是否小于此随机数，若是则接受新候选模型，即得到第次采样时的现有树状结构为：，否则拒绝新候选模型，即得到第次采样时的现有树状结构为：；

12、s7：令，重复s3至s6，并每进行次采样保存一次得到的采样模型，直至采样次数次结束；

13、s8：将s7中保存得到的所有采样模型进行统计计算，得到此区域的地下电导率的一维边缘概率密度分布。

14、可选地，所述s4中的逆变换具体包括：

15、构建尺度函数，以控制当前存在的节点，具体公式如下：

16、；

17、而对于树状结构深度为下的第个节点定义为：

18、；

19、其中，表示为在区间上的单位函数，表示为0-1的任意数无具体意义，表示为节点所在的等级，表示为节点所在该等级的位置；

20、构建解耦函数，以解耦上一等级节点与其下一等级节点的关系，即用其子节点的值覆盖子节点所在范围的父节点的值，具体公式如下：

21、；

22、同理，对于树状结构深度为下的第个节点定义：

23、；

24、其中，表示为在区间上的单位函数。

25、可选地，所述s5包括：

26、s5.1：根据树状结构候选模型逆变换得到的建议模型的正演计算响应和测量得到的电磁场响应进行最小二乘法计算，得到其拟合差，具体公式如下：

27、；

28、其中，代表的是模型向量映射到数据空间中的正演算子，是数据协方差矩阵，为实际测量到的电磁场响应数据；

29、s5.2：假设似然函数服从标准正态分布，似然函数的具体公式如下：

30、；

31、其中，为观测数据的数量。

32、可选地，所述s6包括：

33、s6.1：计算给定二叉树结构下的当前现有树状结构的先验分布 以及新候选模型的先验分布；

34、s6.2：计算第次采样时的现有树状结构模型逆变换得到的电导率模型的似然函数以及第次采样时的新候选树状结构模型经过逆变换得到的建议模型的似然函数 ；

35、s6.3：计算从现有树状结构模型到新候选树状结构模型的建议分布概率以及从新候选树状结构模型到现有树状结构模型的建议分布；

36、s6.4：根据s6.1至s6.3中计算得到的先验分布、似然函数和建议分布计算建议模型的接受率，具体公式如下：

37、；

38、其中，为雅克比矩阵，；

39、s6.5：取0-1之间的随机数，判断接受率是否小于此随机数，若是则接受建议模型，即，否则拒绝建议模型，即。

40、可选地，所述s6.1中包括：

41、s6.1.1：计算给定节点数范围下的均匀分布，具体公式如下：

42、；

43、其中，为树状结构的最大节点数，；为树状结构的最小节点数，一般设置为1，表示至少拥有一个节点；

44、s6.1.2：计算给定节点数下的排列数先验，即为在给定节点的情况下，树状结构的排列组合数的分布情况，具体公式如下：

45、；

46、其中，为个节点数所构成全部不同形状的树状结构的排列数，，为标准二项式系数；

47、s6.1.3：计算给定每个节点中参数的先验分布，具体公式如下：

48、；

49、s6.1.4：采用同一公式分别计算二叉树结构下的先验分布和，其中先验分布的具体公式如下：

50、。

51、可选地，所述s6.3包括：

52、s6.3.1：将当前树状结构上的现有节点构成可扰动节点组，可生成节点构成可生成节点组，可消亡节点构成可消亡节点组；

53、s6.3.2：若s3中选择的节点为现有节点，则计算扰动建议分布，若s3中选择的节点为可生成节点，则计算生成建议分布，若s3中选择的节点为可消亡节点，则计算消亡建议分布。

54、可选地，所述s6.3.2中扰动建议分布的计算包括：

55、计算给定选择第个节点的概率，具体公式如下：

56、；

57、其中是的元素个数；

58、计算给定在第个树节点上选择电导率值概率，具体公式如下：

59、；

60、其中，为参数扰动的正态分布的标准差，为第个树节点上选择电导率变化值；

61、所述s6.3.2中扰动建议分布的具体公式如下：

62、；

63、所述s6.3.2中生成建议分布的具体公式如下：

64、；

65、其中，,且表示选择放置新节点的位置的概率，具体公式如下：

66、；

67、其中，表示为当前模型的可生成节点集合数量；

68、所述s6.3.2中消亡建议分布的计算公式为：

69、；

70、其中，是在提议的出生后可能被删除的节点集，表示为候选模型的可消亡节点集合数量。

71、可选地，所述s6中：

72、选择现有节点的建议模型接受率的表达式为：

73、；

74、选择可生成节点的建议模型接受率的表达式为：

75、；

76、选择可消亡节点的建议模型接受率的表达式为：

77、；

78、其中，表示为当前模型的可消亡节点集合数量，表示为候选模型的可生成节点集合数量。

79、可选地，所述s7中，，不小于1000000，、取自然数。

80、本发明为实现物理空间树状结构，将小波域的树状结构进行改进，在原本的树状结构中其节点系数值表示为小波系数而不具备任何物理意义从而导致我们无法精确的给出有关系数的先验分布情况，这将导致反演得到的电导率值可能超出实际情况，因此本发明将小波域树状结构改为物理空间树状结构，使得树状结构的节点参数是满足具有电性分布特征的先验分布。本发明根据实际测量的地下电磁场响应和建议模型的正演响应计算得到建议模型的似然函数，计算出建议模型的接受率，从而对建议模型进行筛选得到与实际地下电导率模型更加接近的采样模型，且本发明能够任选一类节点改变其电导率参数，从而灵活调整反演的参数，以便于进行多次采样；本发明由于采用了物理空间树状结构，解决了小波域树状结构的节点对先验分布选择困难的问题，达到了自动化参数以及对反演区域多尺度剖分的效果，节省了马尔科夫链对低概率区的搜索，有效地提高了采样效率，也可以延展到对高维贝叶斯反演，为高维贝叶斯反演方法提供了一种自动参数化的思路。