一种基于高斯分布M-H采样的反演模型更新方法与流程

本发明涉及地球物理反演及油气储层预测领域，尤其是一种基于高斯分布m-h采样的反演模型更新方法。

背景技术：

地震反演是预测油气储层的重要步骤，它根据探测器已知的地震记录数据，通过地震记录数据与待求物理量的数学模型来建立最优化问题，并通过反演方法进行求解最优化问题，从而得到待求物理量最优估计的过程。基于蒙特卡洛-马尔科夫链的地震反演技术是地震反演的重要方法，它通过随机采样来更新反演模型，并通过马尔科夫链来完成整个反演的过程。

m-h(metropolis-hastings)采样，是随机采样的重要方法，该方法能保证采样结果具有收敛性；该方法是通过引入状态转移函数对模型进行采样，然后引入判定函数判定采样是否接受，从而通过不断的采样得到最终的采样结果的过程，该方法由metropolis提出，由hastings改进，两人从数学上论证了采样结果收敛的合理性，lavielle将基于均匀分布的m-h采样用于基于高斯分布的贝叶斯地震反演中，证明了m-h采样的在地震反演中的可行性，但现有技术中m-h采样采用均匀分布，均匀分布采用的是0-1之间随机取值，结果呈任何取值分布，导致值域范围过大，收敛速度慢，同时对地震实际数据缺乏分析，导致陷入局部最优；因此需要一种反演模型更新方法可以克服以上问题。

技术实现要素：

本发明的目的在于：本发明提供了一种基于高斯分布m-h采样反演模型更新方法，解决了现有更新反演模型采用基于均匀分布的m-h采样、缺乏对实际数据分析导致收敛值域过大、收敛速度慢、陷入局部最优的问题。

本发明采用的技术方案如下：

一种基于高斯分布m-h采样的反演模型更新方法，包括如下步骤：

步骤1：对地震信息进行预处理获取待反演参数的初始模型；

步骤2：选择第t道待反演参数的初始模型计算初始波阻抗对数l1，将初始波阻抗对数l1作为基于高斯分布m-h采样的初始模型；

步骤3：计算初始波阻抗对数l1的方差，根据方差建立服从高斯分布的状态转移函数，根据状态转移函数的高斯分布对初始波阻抗对数l1进行反复更新并获取更新后的第t道待反演参数；

步骤4：判断t是否大于地震信息中的地震道，若是，则结束更新；若否，令t＝t+1跳至步骤2继续反演。

优选地，所述步骤1包括如下步骤：

步骤1.1：获取地震信息，所述地震数据包括地震数据和测井数据；

步骤1.2：根据地震数据提取层位，对测井数据采用克里金插值建模得到待反演参数初始模型，根据地震数据和测井数据进行井震匹配提取子波。

优选地，所述步骤2包括如下步骤：

步骤2.1：选择第t道待反演参数的初始模型，根据待反演参数和波阻抗对数递推关系求初始波阻抗对数l1，待反演参数和波阻抗对数递推关系如公式1所示：

lk＝f(xk)公式(1)

其中，xk表示待反演参数，lk表示波阻抗对数；

步骤2.2：将初始波阻抗对数l1作为基于高斯分布m-h采样的初始模型。

优选地，所述步骤3包括如下步骤：

步骤3.1：计算初始波阻抗对数l1的均值μk，根据均值μk计算其方差δk，计算公式如公式2所示：

其中，n表示向量l1的元素个数，l1i表示l1中第i个元素；

步骤3.2：根据方差δk建立服从高斯分布的状态转移函数，计算如公式3所示：

f(l1,l2)∝n(0,δk)公式(3)

其中，l1表示初始波阻抗对数，l2表示更新后的波阻抗对数，n(0，δk)表示均值为0，方差为δk的高斯分布；

步骤3.3：根据状态转移函数的高斯分布更新波阻抗对数获得更新后的波阻抗对数l2，计算如公式4所示：

l2＝l1+delta·n公式(4)

其中，l1表示初始波阻抗对数，l2表示更新后的波阻抗对数，n表示与波阻抗对数具有相同元素的服从高斯分布的随机向量，delta表示步长；

步骤3.4：计算判别函数a(l1,l2)，计算如公式5所示：

a(l1,l2)＝min(1,p(l2)/p(l1))公式(5)

其中，p(l1)表示l1的后验概率密度函数，p(l2)表示l2的后验概率密度函数；

步骤3.5：在0-1之间按均匀分布随机选取u，并判断u是否小于a(l1,l2)，若否，则将更新后的波阻抗对数l2作为初始波阻抗对数l1后跳至步骤3.6；若是，则不进行操作；

步骤3.6：判断迭代次数是否大于迭代次数设定值，若否，则跳至步骤3.1进行重复更新；若是，则根据更新后的波阻抗对数计算更新后的第t道待反演参数。

优选地，所述步骤3.4中后验概率密度函数计算如公式6所示：

其中，g表示根据子波获取的常数矩阵，l表示根据层位获取的第t道波阻抗对数，s0表示第t道原始地震数据，c1为常数项。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1.本发明基于均匀分布m-h采样改为基于正态分布m-h采样，高斯分布通过统计井插值外扩后的波阻抗对数数据的分布得到，符合实际数据的分布，解决了现有更新反演模型采用基于均匀分布的m-h采样、缺乏对实际数据分析导致收敛值域过大、收敛速度慢、陷入局部最优的问题，达到了缩小收敛值域范围、快速寻优的效果；

2.本发明通过计算波阻抗对数方差，基于方差建立服从高斯分布的状态转移函数，通过基于高斯分布的m-h采样进行更新，高斯分布采样最后的结果仍然服从高斯分布，值域不变，避免了均匀分布采样结果为任何值域导致收敛值域过大的缺点，加快收敛。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1为本发明的方法流程框图；

图2为本发明的单道井旁道波阻抗的直方图分布与对应的高斯分布示意图；

图3为本发明的含噪声的地震剖面；

图4为本发明的多道波阻抗初始模型示意图；

图5为本发明的多道波阻抗更新结果示意图；

图6为本发明的单道井旁道波阻抗更新结果及对比示意图；

图7为本发明的状态转移函数采用均匀分布和高斯分布的误差对比图；

图8为本发明的方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明，即所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，术语“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

下面结合图1-8对本发明作详细说明。

本发明解决的技术问题：解决了现有更新反演模型采用基于均匀分布的m-h采样、缺乏对实际数据分析导致收敛值域过大、收敛速度慢、陷入局部最优的问题

技术手段：一种基于高斯分布m-h采样的反演模型更新方法，包括如下步骤：

步骤1：对地震信息进行预处理获取待反演参数的初始模型；

步骤2：选择第t道待反演参数的初始模型计算初始波阻抗对数l1，将初始波阻抗对数l1作为基于高斯分布m-h采样的初始模型；

步骤3：计算初始波阻抗对数l1的方差，根据方差建立服从高斯分布的状态转移函数，根据状态转移函数的高斯分布对初始波阻抗对数l1进行反复更新并获取更新后的第t道待反演参数；

步骤4：判断t是否大于地震信息中的地震道，若是，则结束更新；若否，令t＝t+1跳至步骤2继续反演。

步骤1包括如下步骤：

步骤1.1：获取地震信息，所述地震数据包括地震数据和测井数据；

步骤1.2：根据地震数据提取层位，对测井数据采用克里金插值建模得到待反演参数初始模型，根据地震数据和测井数据进行井震匹配提取子波。

步骤2包括如下步骤：

步骤2.1：选择第t道待反演参数的初始模型，根据待反演参数和波阻抗对数递推关系求初始波阻抗对数l1，待反演参数和波阻抗对数递推关系如公式1所示：

lk＝f(xk)公式(1)

其中，xk表示待反演参数，lk表示波阻抗对数；

步骤2.2：将初始波阻抗对数l1作为基于高斯分布m-h采样的初始模型。

步骤3包括如下步骤：

步骤3.1：计算初始波阻抗对数l1的均值μk，根据均值μk计算其方差δk，计算公式如公式2所示：

其中，n表示向量l1的元素个数，l1i表示l1中第i个元素；

步骤3.2：根据方差δk建立服从高斯分布的状态转移函数，计算如公式3所示：

f(l1,l2)∝n(0,δk)公式(3)

其中，l1表示初始波阻抗对数，l2表示更新后的波阻抗对数，n(0，δk)表示均值为0，方差为δk的高斯分布；

步骤3.3：根据状态转移函数的高斯分布更新波阻抗对数获得更新后的波阻抗对数l2，计算如公式4所示：

l2＝l1+delta·n公式(4)

其中，l1表示初始波阻抗对数，l2表示更新后的波阻抗对数，n表示与波阻抗对数具有相同元素的服从高斯分布的随机向量，delta表示步长；

步骤3.4：计算判别函数a(l1,l2)，计算如公式5所示：

a(l1,l2)＝min(1,p(l2)/p(l1))公式(5)

其中，p(l1)表示l1的后验概率密度函数，p(l2)表示l2的后验概率密度函数；

步骤3.5：在0-1之间按均匀分布随机选取u，并判断u是否小于a(l1,l2)，若否，则将更新后的波阻抗对数l2作为初始波阻抗对数l1后跳至步骤3.6；若是，则不进行操作；

步骤3.6：判断迭代次数是否大于迭代次数设定值，若否，则跳至步骤3.1进行重复更新；若是，则根据更新后的波阻抗对数计算更新后的第t道待反演参数。

步骤3.4中后验概率密度函数计算如公式6所示：

其中，g表示根据子波获取的常数矩阵，l表示根据层位获取的第t道波阻抗对数，s0表示第t道原始地震数据，c1为常数项；技术效果：本发明基于均匀分布m-h采样改为基于正态分布m-h采样，高斯分布通过统计井插值外扩后的波阻抗对数数据的分布得到，符合实际数据的分布，解决了现有更新反演模型采用基于均匀分布的m-h采样、缺乏对实际数据分析导致收敛值域过大、收敛速度慢、陷入局部最优的问题，达到了缩小收敛值域范围、快速寻优的效果；通过计算波阻抗对数方差，基于方差建立服从高斯分布的状态转移函数，通过基于高斯分布的m-h采样进行更新，高斯分布采样最后的结果仍然服从高斯分布，值域不变，避免了均匀分布采样结果为任何值域导致收敛值域过大的缺点，加快收敛。

实施例1

如图1-8所示，待反演参数为波阻抗时反演如下：

步骤1.1：输入地震信息，根据地震数据提取层位，根据测井数据采用克里金插值建模得到波阻抗初始模型，通过地震信息采用井震匹配提取子波；

步骤2.1：选择第t道模型的波阻抗aik，通过波阻抗与波阻抗对数的递推关系求初始波阻抗对数，波阻抗与波阻抗对数的递推关系如公式1所示：

lk＝ln(aik)公式(1)

其中，aik表示波阻抗，lk表示波阻抗对数；

步骤2.2：将初始波阻抗对数l1作为更新初始模型；

步骤3.1：计算波阻抗对数lk的均值μk，根据均值μk计算方差δk，计算公式如公式2所示：

其中，n表示向量l1的元素个数，l1i表示l1中第i个元素；

步骤3.2：建立服从高斯分布的状态转移函数，计算如公式3所示：

f(l1,l2)∝n(0,δk)公式(3)

其中，l2表示更新后的波阻抗对数，n(0，δk)表示均值为0，方差为δk的高斯分布；

步骤3.3：根据状态转移函数的高斯分布更新波阻抗对数获得更新后的波阻抗对数l2，计算如公式4所示：

l2＝l1+delta·n公式(4)

其中，l1表示初始波阻抗对数，l2表示更新后的波阻抗对数，n表示与波阻抗对数具有相同元素的服从高斯分布的随机向量，delta表示步长；

步骤3.4：计算判别函数a(l1,l2)，计算如公式5所示：

a(l1,l2)＝min(1,p(l2)/p(l1))公式(5)

其中，p(l1)表示l1的后验概率密度函数，p(l2)表示l2的后验概率密度函数；

后验概率密度函数计算如公式6所示：

其中，g表示根据子波获取的常数矩阵，l表示根据层位获取的第t道波阻抗对数，s0表示第t道原始地震数据，c1为常数项；

步骤3.5：在0-1之间按均匀分布随机选取u，并判断u是否小于a(l1,l2)，若否，则将更新后的波阻抗对数l2作为初始波阻抗对数l1后跳至步骤3.6；若是，则不进行操作；

步骤3.6：判断迭代次数是否大于迭代次数设定值，若否，则跳至步骤3.1进行重复更新；若是，则根据更新后的波阻抗对数计算更新后的第t道待反演参数。

步骤4包括如下步骤：

步骤4.1：判断t是否大于输入数据中的地震道数据，若是，则结束更新得到的更新后的待反演参数剖面；若否，令t＝t+1跳至步骤2继续反演。

效果分析：如图2所示，高斯分布能够很好地模拟插值后的波阻抗对数，采用高斯分布模拟m-h采样具有合理性；如图3-5所示，更新后的波阻抗剖面的结果相较于波阻抗剖面初始模型能够更好地服从地震记录的趋势，反应了本申请的正确性，如图6所示initial表示波阻抗初始模型，inversion表示更新后的波阻抗结果，well表示真实的波阻抗，更新结果比初始模型更接近真实数据，进一步反应了本申请的正确性；如图7所示u表示采用均匀分布m-h采样下的收敛次数与误差的关系，n表示高斯分布m-h采样下的收敛次数与误差的关系，图中所示：本申请的方法减小了收敛值域，能更快的寻优，本发明首先通过分析源于测井数据建模获得的模型数据的统计特性，结合高斯分布构建m-h采样中的状态转移函数，对初始模型采用上述中的状态转移函数形成的高斯分布进行采样并通过判别函数更新，通过反复迭代，输出最优的反演结果，本发明解决了现有地震反演采用均匀分布作为状态转移函数导致收敛值域过大、收敛速度慢、陷入局部最优的问题，达到了缩小收敛值域范围、快速寻优的效果。