采用自适应泰森多边形参数化的瞬变电磁二维贝叶斯反演方法

1.本发明涉及地球探测信息技术领域，尤其涉及一种采用自适应泰森多边形参数化的瞬变电磁二维贝叶斯反演方法。


背景技术：

2.瞬变电磁法(transient electromagnetic method,tem)又可称时间域电磁法(time-domain electromagnetic method,tdem)是地球物理电磁法勘探的重要分支，一般用接地导线或回线向地下空间发射一次电磁场，从而激发地下介质而产生二次电磁场，关断一次场后观测纯二次场来推断地下电性结构分布。瞬变电磁收发装置灵活多变，如有地面回线源装置、接地导线源装置、半航空装置以及航空装置等。因此其应用领域广泛，深部的地球结构研究、中部的矿产资源勘探以及浅部的环境工程勘察均可覆盖。
3.瞬变电磁三维反演是目前主要的研究热点，近年来国内外开发了多个瞬变电磁三维正反演算法，然而瞬变电磁三维反演的实际应用受到计算资源与计算时间的双重影响而并没有较好推广。拟二维或拟三维反演是目前应用最广也是最适合瞬变电磁快速解释的方案，此时采用一维模型的模拟结果去拟合数据，但模型通过二维或三维形式表示。这类方法在许多二、三维结构影响不强的实际勘探中获得了有效应用，尤其在沉积环境下更加适合。中心回线作为瞬变电磁主流的勘探装置，由于发射源的尺寸一般较小，信号受二、三维地电结构的影响较小。因此中心回线瞬变电磁的拟二维或三维反演结果多数情况下能较好地对地下电性结构进行有效反映。
4.目前拟二维反演最主流的为采用横向约束的反演(lci)，通常沿测线对相邻模型的电阻率或层厚进行约束。当采用最光滑或者最小结构模型约束时，如果出现模型过度参数化的情况，就需要采用模型正则化来稳定反演结果。这种反演结果将会受到初始模型和所加的模型约束的影响，例如为了获得一个更加可靠的反演结果，需要采用多次试错的方式去确定横向约束的程度。lci反演为线性化的最优化问题，通过在一定范围内对定义的目标函数进行求极小值而获得一个或一套反演结果模型。尽管可以通过对反演结果模型进行线性化的灵敏度分析来估计模型的不确定度，但是线性化过程会丢失一些有价值的信息，导致反演的不确定度估计出现一定偏差。由于缺少足够的信息去评价反演结果模型的质量，将影响瞬变电磁数据的解释。
5.为严格估计反演结果的不确定度，可以采用贝叶斯方法去求解非线性反演问题。在贝叶斯方法框架下，模型参数被作为随机变量处理，反演结果采用模型参数的后验概率密度(pdf)来表示。通过探索后验的pdf，可以对模型参数的不确定度进行有效量化。基于这一优势，近十年来，贝叶斯方法在电磁数据的反演问题中获得越来越多的应用，部分研究人员也将贝叶斯单独反演或联合反演引入到了tem数据的反演中。
6.贝叶斯方法需要产生大量的模型去估计反演模型参数的pdf，会涉及大量的正演运算。因此瞬变电磁的贝叶斯反演绝大多数是基于一维假设，二、三维正演的计算资源和时
间消耗都远大于一维情况，难以获得实际应用。


技术实现要素：

7.本发明要解决的技术问题在于针对现有技术中的缺陷，提供一种采用自适应泰森多边形参数化的瞬变电磁二维贝叶斯反演方法，将基于自适应泰森多边形网格的贝叶斯反演引入到瞬变电磁的数据处理中，旨在提供一种能更好估计反演模型参数不确定度的反演方法，为瞬变电磁数据的合理解释提供一种新方法。
8.本发明解决其技术问题所采用的技术方案是：
9.本发明提供一种采用自适应泰森多边形参数化的瞬变电磁二维贝叶斯反演方法，该方法包括以下步骤：
10.步骤1、根据已有信息及条件设定相关反演参数，包括：(1)电阻率范围[ρ
min
,ρ
max
]；(2)模型voronoi网格单元数:[k
min
,k
max
]；(3)模型水平位置范围：[x
min
,x
max
]；(4)模型深度范围：[z
min
,z
max
]；
[0011]
步骤2、根据步骤1中的反演参数作为先验条件随机生成初始模型m0；
[0012]
步骤3、设定模型的修改方式：(1)在网格中产生一个新的voronoi网格单元；(2)在网格中减少一个已有voronoi网格单元；(3)随机改变一个voronoi网格单元电阻率；(4)随机改变一个voronoi网格单元尺寸；
[0013]
步骤4、随机挑选步骤3中的一种修改方式，生成新的建议模型m1；
[0014]
步骤5、判断建议模型m1是否满足步骤1所定义的先验准则；
[0015]
步骤6、如果模型m1不满足先验准则，则拒绝该建议模型，设置m0＝m0，并跳回到步骤3重新生成新的建议模型；
[0016]
步骤7、如果模型m1满足先验准则，则进行一维正演模拟计算所有测点处模型m1的瞬变电磁响应f(m)；
[0017]
步骤8、采用步骤7中计算获得的瞬变电磁响应f(m)，计算数据的似然度p(d|m1)；
[0018]
步骤9、判断步骤8所计算的似然度p(d|m1)是否满足mh接受准则；若不满足，则跳到步骤6；
[0019]
步骤10、在步骤9中若似然度p(d|m1)满足mh接受准则，则接受该新建议模型，并设置m0＝m1。
[0020]
进一步地，本发明的所述步骤1中的方法具体为：
[0021]
根据采集数据、地质情况以及已有信息来设置相关先验信息，其中包括电阻率的分布范围[ρ
min
,ρ
max
]，设置模型的voronoi网格单元总数的范围[k
min
,k
max
]，设置模型的水平位置范围[x
min
,x
max
]和深度范围[z
min
,z
max
]；设置瞬变电磁正演的相关参数信息，包括回线尺寸大小、时间道范围。
[0022]
进一步地，本发明的所述步骤2中的方法具体为：
[0023]
根据步骤1中所设置的相关先验信息随机确定一个二维模型m0作为贝叶斯反演的起始模型，该二维模型采用voronoi多边形网格单元进行定义；对于给定的数据，模型网格的相关参数是未知的，因此voronoi单元的尺寸和数量在反演过程中是允许动态变化的；该二维模型的相关参数通过以下方式进行定义：
[0024]
m＝(x,z,ρ)
[0025]
其中x＝[x1,x2,
…
,xn]和z＝[z1,z2,
…
,zn]分别为每个voronoi网格单元核心的xz坐标位置；ρ＝[log
10
ρ1,log
10
ρ2,
…
,log
10
ρn]为每个网格单元的电阻率参数，为保证电阻率的非负性，采用了对数电阻率的形式；n为voronoi网格单元的数量；通过这种动态模型定义形式，voronoi网格单元在反演过程中表示任意形状和尺寸的二维模型；此时，模型的电阻率分布能自动适用于数据分辨率。
[0026]
所述步骤7中的具体方法为：
[0027]
对新的建议模型m1进行正演模拟，为保证正演计算效率，尽管建议模型m1为二维模型，但采用一维正演算法进行各测点的瞬变电磁响应计算；各测点处的一维模型采用最邻近插值方法从二维模型获取；
[0028]
一维瞬变电磁模拟，采用频率域转时间域算法进行计算；其中频率域响应采用偶极子响应结合高斯积分进行求解，然后采用正余弦变换转换到时间域获得瞬变电磁响应。
[0029]
进一步地，本发明的所述步骤8中的方法具体为：
[0030]
贝叶斯反演方法的结果为模型参数的后验概率分布；通过贝叶斯理论，后验分布要结合模型参数的先验信息以及观测数据的信息，定义为如下形式：
[0031][0032]
其中，p(m)和p(d|m)分别为先验分布与似然函数；p(d)为边际似然函数，作为后验分布的归一化常数；后验分布p(m|d)依赖于两种概率，即p(m)和p(d|m)；因此，在贝叶斯反演过程中要进行两种工作：(1)对先验信息进行概率性量化；(2)估计给定模型对观测数据的拟合程度概率；
[0033]
先验分布p(m)是与观测数据无关的模型参数的概率测度，用于控制模型空间范围；当没有明确的先验信息时，常假设对模型参数的先验信息在一定范围内服从于均匀分布或高斯分布；
[0034]
似然函数p(d|m)用于量化由于数据误差引起的模型模拟数据和观测数据之间的差异概率，其形式取决于观测误差的分布；假设数据误差服从高斯分布并且是不相关的，因此似然函数定义为如下形式：
[0035][0036]
其中n为观测数据的数量，cd为用于量化数据的不确定性和数据相关性的数据协方差矩阵，影响估计参数的不确定水平；对于不相关的数据误差，cd为数据方差的对角矩阵。
[0037]
进一步地，本发明的所述步骤9中判断似然度是否满足mh接受准则的方法具体为：
[0038]
贝叶斯方法的目标是通过合理的采样方式获得一系列拟合观测数据的模型，从而估计模型参数的目标后验概率分布；采用贝叶斯方法中的马尔科夫链蒙特卡洛mcmc方法；mcmc方法的实现算法采用mh采样算法，通过迭代的方式采用两步过程去构建马尔科夫链；在每次的迭代中，首先根据建议概率分布q(m1|m)对当前模型m进行随机扰动获得一个候选模型m1；然后，通过mh准则去判断是否接受该候选模型；采用由mh算法扩展而来的可逆跳动mcmc跨维抽样算法进行抽样，允许模型参数的维度变化；
[0039]
在每次迭代过程中由跨维mcmc算法产生的候选模型的结构概率可通过以下metropolis-hasting-green准则，即mhg准则来表示：
[0040][0041]
其中以及分别为模型m1和m的先验概率、似然概率以及建议概率的比值；|j|为从m到m1的模型转化雅克比矩阵，用于量化模型维度变化的影响，当两个模型的维度相同时，该雅克比矩阵为1；采用birth-death方案来控制模型的维度变化，在马尔科夫链的每一步中对模型维度的改变为1，此时雅克比矩阵为单位阵即|j|＝1，并不会对mhg准则方程产生影响；
[0042]
判断是否满足mh接受准则具体方式：利用mh公式计算建议模型m1的接受率α，然后利用随机函数生成一个在[0,1]间的均匀随机数，当接受率α大于该随机数时，接受该建议模型m1，否则拒绝该模型。
[0043]
本发明产生的有益效果是：本发明的采用自适应泰森多边形参数化的瞬变电磁二维贝叶斯反演方法，尤其适用于地下介质分层性较好的瞬变电磁探测数据解释，如沉积地质条件。本发明将自适应泰森多边形网格与贝叶斯方法相结合引入到瞬变电磁的二维反演中，旨在为瞬变电磁剖面数据提供一个可以有效计算模型参数不确定度的反演方法。较传统的横向约束反演方法本发明能在反演过程中自动调整模型参数的个数及尺寸去更好地适应观测数据的空间分辨率，从而避免了对模型参数与正则化参数的先验选择，同时贝叶斯方法能提供反演模型的不确定信息并更准确地确定实测数据的有效探测深度。
附图说明
[0044]
下面将结合附图及实施例对本发明作进一步说明，附图中：
[0045]
图1是本发明实施例的voronoi网格二维模型示意图；(a)voronoi网格单元的几何结构由黑色圆圈表示的核心定义；(b)由二维模型的voronoi网格单元定义的一维模型示意图；
[0046]
图2是本发明实施例的合成二维测试模型示意图；倒三角形为相距20m的tem测点；
[0047]
图3是本发明实施例的二维测试模型的tem响应；(a)部分时间道的tem响应沿测点的分布；(b)部分测点的tem响应衰减曲线；
[0048]
图4是本发明实施例的二维贝叶斯反演生成的后验模型合集中随机挑选的四个voronoi网格模型；各模型的voronoi网格数和数据拟合差(rms)均标注于图中；
[0049]
图5是本发明实施例的反演结果图，(a)二维贝叶斯反演生成的后验模型集合的rms统计柱状图；(b)模型集合的voronoi网格数的统计柱状图；
[0050]
图6是本发明实施例的反演结果图，(a-d)为二维贝叶斯反演的模型参数统计分布；(e-h)为一维贝叶斯反演的模型参数统计分布；
[0051]
图7是本发明实施例的不同深度处电阻率的边际概率分布，图中黑色虚线为真实模型电阻率；
[0052]
图8是本发明实施例的方法流程图。
具体实施方式
[0053]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
[0054]
本发明实施例中涉及以下名词解释：
[0055]
瞬变电磁法：瞬变电磁法又可称为时间域电磁法，其英文为transient electromagnetic method简称为tem，是地球物理电磁法的重要分支，一般有回线源发射和接地源发射两种形式。该方法具有勘探效率高、成本低以及精度较高等特点；
[0056]
贝叶斯方法：英文为bayesian approach，该方法是基于贝叶斯概率统计思想，将参数估计问题转化为概率推断问题，获得模型参数所服从的后验概率分布，从而对反演结果的不确定度进行定量的分析与评价。
[0057]
泰森多边形：英文为voronoi，是一组由连接两邻点线段的垂直平分线组成的连续多边形。一个泰森多边形内的任一点到构成该多边形的控制点的距离小于到其他多边形控制点的距离。在本发明中用于模型的参数化。
[0058]
正演：从模型参数到响应数据的正向计算过程；
[0059]
反演：从数据到模型参数估计的反向计算过程；
[0060]
其中下述大写加粗的字母表示向量或者矩阵。
[0061]
一、瞬变电磁一维正演模拟
[0062]
本发明需要涉及瞬变电磁法的一维正演模型，用于反演过程中计算模型的响应。瞬变电磁仅有少量特殊模型具有纯解析解，而普通的一维层状模型瞬变电磁响应需要采用频谱法进行计算。频谱法即频率域转时间域算法，本发明采用正余弦变换算法具体实现。对于频率域电磁场，考虑正谐时条件e
iωt
，下阶跃电流激发情况下的频率域响应转换到时间域响应的正余弦变换公式为
[0063][0064][0065]
其中ω为角频率，f(ω)和f(t)分别为频率域和时间域电磁场响应，re[]和im[]分别为求取复数的实部和虚部的算符。
[0066]
考虑正谐时条件e
iωt
，对于沿x方向的水平电偶极子，其频率域响应的垂直磁场分量可采用下式进行计算
[0067][0068]
其中p为电偶极子强度，p＝idx，dx为电偶极矩长度。其中p为电偶极子强度，p＝idx，dx为电偶极矩长度。k2＝ω2με-iωμσ；k
x
和ky分别为波数域波数λ在x和y轴方向上的分量；r
te
为一维层状模型介质表面的反射系数，与各层电导率和厚度相关，需要从底层逐渐往地面迭
代计算；为场点到偶极子源的水平距离，j1为一阶贝塞尔(bessel)函数。
[0069]
将发射回线进行剖分，采用(3)式计算各偶极子的频率域电磁场响应，再进行求和，即可得到发射回线的频率域电磁场响应，然后利用式(1)或(2)即可计算瞬变电磁响应。
[0070]
二、采用voronoi单元的二维模型参数
[0071]
本发明采用voronoi剖分对二维电阻率模型进行参数化。如附图1所示，一个二维模型的voronoi剖分是通过一组节点将二维电阻率模型划分为一系列不重叠的区域即voronoi单元。voronoi单元是一种非结构化剖分，能灵活地对复杂不规则的几何模型进行剖分。基于这一优势，voronoi布局在跨维贝叶斯反演中经常用于模型剖分和参数化。
[0072]
由于从观测数据难以先验地确定合适的模型剖分方式，因此voronoi单元的数目以及位置在反演过程中均作为可以变化的未知参数。因此，二维电阻率模型的模型参数定义为如下形式
[0073]
m＝(x,z,ρ)∈r
3n
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0074]
其中n为voronoi单元的数目，是一个需要确定的未知量。x＝[x1,x2,
…
,xn]和z＝[z1,z2,
…
,zn]分别表示voronoi单元核的x与z坐标。ρ＝[log
10
ρ1,log
10
ρ2,
…
,log
10
ρn]为每个网格单元的电阻率参数，为保证电阻率的非负性，采用了对数电阻率的形式。通过这种动态模型剖分方式，voronoi网格单元可以在反演过程中表示任意形状和尺寸的二维模型。此时，模型的电阻率分布能自动适用于数据分辨率。
[0075]
三、贝叶斯推断
[0076]
地球物理电磁法反演问题的目的是寻找一个或一组能够拟合观测数据且与先验信息吻合的电阻率模型。与基于最优化过程的传统反演方法不同，贝叶斯反演方法的结果为模型参数的后验概率分布。通过贝叶斯理论，后验分布需要结合模型参数的先验信息以及观测数据的信息，可以定义为如下形式
[0077][0078]
其中，p(m)和p(d|m)分别为先验分布与似然函数。p(d)为边际似然函数，作为后验分布的归一化常数。从上式可以看出，后验分布p(m|d)依赖于两种概率，即p(m)和p(d|m)。因此，在贝叶斯反演过程中需要进行两种工作：(1)对先验信息进行概率性量化；(2)估计给定模型的观测数据概率。
[0079]
由于贝叶斯反演过程需要利用观测数据的信息来更新先验概率，因此先验分布p(m)是不可或缺的，相当于梯度法反演中的正则化项，用于控制模型空间范围，同时表示为独立于观测数据的模型参数的概率测度。当没有明确的先验信息时，常假设对模型参数的先验信息在较大范围内服从于均匀分布或高斯分布。
[0080]
似然函数p(d|m)用于量化由于数据误差引起的模型模拟数据和观测数据之间的差异概率，其形式取决于观测误差的分布。例如，当观测数据误差服从于期望值为零的正态分布时，通常将似然函数表示为高斯分布，当不服从于高斯分布时，假设其服从于拉普拉斯分布更为合理。本专利中，假设数据误差服从与高斯分布并且是不相关的，因此似然函数可定义为如下形式
[0081][0082]
其中n为观测数据的数量，cd为用于量化数据的不确定性和数据相关性的数据协方差矩阵，影响估计参数的不确定水平。对于不相关的数据误差，cd为数据方差的对角矩阵。
[0083]
四、跨维mcmc采样
[0084]
贝叶斯方法的目标是通过合理的采样方式获得一系列能够拟合观测数据的模型从而去评价模型参数的目标后验概率分布。在各种不同的贝叶斯方法中，最常采用的为mcmc方法，由于这种方法的相对简单高效，被广泛应用于从模型空间中获得候选样本。mcmc方法将蒙特卡洛模拟与马尔科夫链技术联合起来，其中前者通过随机扰动模型参数获得模型样本，而后者通过一定的判断根据来指导这些随机模型的生成。
[0085]
最有名的mcmc实现算法为metropolis-hastings(mh)采样算法，这种算法通过迭代的方式采用两步过程去构建马尔科夫链。在每次的迭代中，首先根据建议概率分布q(m1|m)对当前模型m进行随机扰动获得一个候选模型m1。然后，通过mh准则去判断是否接受该候选模型。为了从跨维后验概率中进行抽样，需要mcmc算法能够在保持详细平衡条件的情况下能够不同维度的多模型子空间，而mh算法通常用于反演过程中模型参数不变的情况。因此，本专利采用由传统mh算法扩展而来的可逆跳动mcmc跨维抽样算法进行抽样，可以允许模型参数的维度变化。
[0086]
相较于传统的mh算法，在每次迭代过程中由跨维mcmc算法产生的候选模型的结构概率可通过以下metropolis-hasting-green准则(mhg)来表示
[0087][0088]
其中以及分别为模型m1和m的先验概率、似然概率以及建议概率的比值。|j|为从m到m1的模型转化雅克比矩阵，用于量化模型维度变化的影响，当两个模型的维度相同时，该雅克比矩阵为1。本专利中，我们采用较为常用的birth-death方案来控制模型的维度变化，在马尔科夫链的每一步中可以对模型维度的改变为1，此时雅克比矩阵为单位阵即|j|＝1，并不会对mhg准则方程产生影响。
[0089]
五、发明方法的模拟数据测试
[0090]
为验证本发明方法的应用效果，采用模拟数据测试，测试模型如图2所示为一深谷模型，相关的模型参数均已标注于图中。为获得合成模型数据，布置了一条点距为20m的测线，共61个测点。发射回线为边长为40m的矩形线框，观测时间道从2μs到1ms。采用瞬变电磁三维有限体积算法计算该模型各测点处的瞬变电磁响应，在响应中添加3％的高斯白噪声，并设置环境的绝对误差噪声水平为1nv/am2，从而构建了观测数据(如图3所示)。
[0091]
将本发明方法应用于该观测数据，初始参数设置voronoi单元数目：[30,150]；电阻率范围：[1,10000]ωm；voronoi单元核心的x坐标范围：[-650,650]m；voronoi核的z坐标范围[0,500]m。本次测试共设置了16条独立的马尔科夫链同时并行探索模型空间，每一条
链的初始模型均根据先验条件随机产生，并最终共生成100万个随机样本。每条链的前一半模型作为初始预热采样而被舍弃，之后每隔100个模型挑选一个模型形成聚合链，然后将所有聚合链集合起来估计模型参数的后验概率密度函数。
[0092]
贝叶斯反演结果的表示方式采用的是根据模型参数后验分布而构建的模型集合。由于本发明采用的是动态的模型参数，因此模型集合中将包括多种不同的voronoi数目模型，且模型复杂程度存在差异。附图4展示了4个随机抽选的后验模型集合中的模型，尽管4个模型具有相似的数据拟合差，但4个模型的voronoi单元数目，以及模型集合形态均有较大差别。这说明反演结果中存在多个能拟合数据的等效模型。图5展示了所有反演模型集合中数据拟合差位于1.3与1.4之间的模型分布直方图，以及这些模型的voronoi单元数目分布直方图。
[0093]
一旦获得反演的后验概率密度函数后，便可推断地下电阻率结构的不同统计特征。附图6所示，左侧为二维贝叶斯反演结果的均值模型(mean)、中值模型(median)、最大后验概率模型(mode)以及模型的标准差。尽管反演模型集合中的每个模型有其特定的模型特征，也许并不具有地质意义，但所有模型的统计量能对地下电阻率结构进行较好推断。这是因为不同参数化的结果模型整合了后验模型集合的信息，能够去除单个模型中那些不符合实际情况的特征，而且模型结果的统计信息能够提供更高的空间分辨率。作为对比，各测点的一维贝叶斯反演结果也在附图6中进行了展示。
[0094]
除了从反演结果模型中提取的统计信息外，反演模型的完整集合也能够给出地下结构界面信息。附图7展示了完整的后验概率密度分布在不同深度处的水平切片。图中的这些边际概率分布清晰显示了真实模型界面处的电阻率变化。
[0095]
应当理解的是，上面结合附图对本发明的实施进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。