一种基于频率域协模拟的叠前地震随机反演方法与流程

本发明属于地球物理油气勘探，具体地说涉及一种基于频率域协模拟的叠前地震随机反演方法。

背景技术：

1、地震反演可以利用地震数据来预测储层弹性参数、物性参数、岩性、孔隙度、含水饱和度等相关信息，为油气资源评价奠定基础。反演问题通常是一个非线性问题，由于地震反演方法本身的局限性，在许多情况下使用常规线性反演方法来对反演问题进行求解存在较大误差，且反演结果的分辨率难以满足薄储层预测的要求。与线性反演方法相比，随机反演不需要获得反演问题的解析解，且可以有效地利用测井资料中的高频信息，提高反演结果的垂直分辨率，为识别薄储层提供有力的技术支撑。

2、随机反演方法通常将反演问题视为一种迭代优化问题，其利用随机模拟从先验分布中提取样本，并采用优化算法（例如metropolis-hastings、模拟退火）对样本进行优选，最终获得收敛于后验概率分布的反演结果。随机模拟是随机反演的核心，其制约了反演的精度和效率。按照模拟方法来说，主要分为两类，其一即顺序模拟方法，主要包括序贯高斯模拟、直接序贯模拟等。顺序模拟方法在测井数据的基础上，通过分析参数在空间上展布趋势，利用加权平均法计算每一个网格点处的模拟均值与方差，从而建立高分辨率储层参数模型。然而，该类方法计算量大，在缺乏井数据的情况下难以建立可靠的储层参数模型。另一种随机模拟方法是频率域模拟方法，如快速傅里叶滑动平均（fft-ma）模拟、随机介质建模等。该类方法通过从地震数据中估计协方差和自相关矩阵，可以建立反映地震空间结构特征的储层参数模型。与传统的序列高斯模拟相比，该类方法具有更高计算效率。然而，频率域模拟方法难以有效实现多参数之间的协模拟，这制约了其在叠前ava反演中的应用。

3、优化算法可以对随机模拟所生成的样本进行优选，通过不断迭代，使其逐步收敛于后验概率分布。马尔科夫链蒙特卡洛（mcmc）算法及其改进算法、遗传算法和模拟退火等算法是在随机反演中常用的优化方法。然而，这些算法均为全局优化算法，计算效率较低。

技术实现思路

1、针对现有技术的种种不足，发明人在长期实践中研究设计出一种基于频率域协模拟的叠前地震随机反演方法，将fft-ma协模拟和迭代更新方法相结合，提高叠前随机反演的精度与计算效率。

2、本发明的随机反演方法包括以下步骤：

3、步骤一，利用测井资料，分别建立纵横波速比与纵波模量、纵横波速比与密度的二元联合概率分布。

4、步骤二，利用fft-ma模拟建立初始纵横波速比模型。

5、步骤三，在纵横波速度比模型的基础上，结合纵横波速比与纵波模量、纵横波速比与密度的联合概率分布，计算纵波模量与密度的取值范围。

6、步骤四，以纵波模量与密度的取值范围作为约束，利用fft-ma协模拟建立初始纵波模量与密度模型。

7、步骤五，根据贝叶斯定理，建立目标函数公式，并计算初始目标函数值。

8、步骤六，在现有纵横波速比模型的基础上，利用迭代更新算法产生新的纵横波速比模型。

9、步骤七，在新纵横波速比模型的基础上，计算每个网格点处纵波模量与密度的取值范围。

10、步骤八，在现有纵波模量与密度模型的基础上，将新的纵横波速比模型作为软数据，利用迭代更新算法产生新的纵波模量与密度模型。

11、步骤九，计算新模型所对应的目标函数值，并与原有模型对应的目标函数值进行对比，利用metropolis准则来判断是否保留新模型。

12、步骤十，重复步骤六~步骤九，直到达到最大迭代次数或目标函数值小于门槛值。

13、进一步地，步骤一中，分别建立纵横波速比与纵波模量、纵横波速比与密度的二元联合概率分布，具体方法为：

14、根据测井纵波模量、纵横波速度比与密度曲线，分别建立纵横波速比与纵波模量的二元联合概率分布、纵横波速比与密度的二元联合概率分布。

15、进一步地，步骤二中，利用fft-ma模拟建立初始纵横波速比模型，具体方法为：

16、基于fft-ma模拟方法，模拟模型可以表征为以下形式：

17、                                  公式1

18、                                      公式2

19、式中，为模拟模型的均值，为测井数据的标准差，表示协方差函数的共轭根，表示高斯白噪声，其均值为0，方差为1。

20、进一步地，步骤三中，计算纵波模量与密度的取值范围，具体方法为：

21、将纵横波速比作为次级变量，将纵波模量与密度作为主变量。在纵横波速比模拟结果的基础上，定义一个区间，因此主变量的条件概率分布可以表征为：

22、              公式3

23、其中表示给定的次级变量区间长度。对条件概率分布进行高斯分布拟合，即

24、                                 公式4

25、式中，为高斯分布的均值，为方差。因此，当次级变量模拟结果为时，主变量的取值范围为。

26、进一步地，步骤四中，利用fft-ma协模拟建立初始纵波模量与密度模型，将纵横波速度比作为次级变量，分别将纵波模量和密度作为主变量，分别进行以下计算：

27、子步骤4.1，首先通过fft-ma模拟对次级变量进行模拟，获得次级变量的模拟结果，然后建立以下的协模拟公式来获得主变量的模拟结果，实现多个反演参数的协同模拟：

28、                                公式5

29、                          公式6

30、式中，为主变量模拟结果的均值，与分别为主变量与次级变量的方差，为主变量与次级变量的相关系数；

31、子步骤4.2，根据主变量的取值范围为，此时最终协模拟结果可写为：

32、              公式7

33、进一步地，步骤五中，计算初始目标函数值，包括以下子步骤：

34、子步骤5.1：根据zeoppritz近似式与地震褶积模型，建立反演参数与地震记录的正演关系。在aki-richards近似式基础上，推导基于纵波模量、纵横波速比与密度的地震反射系数方程：

35、 公式8

36、式中，为入射角；、、分别代表地下纵波模量、纵横波速比和密度的平均值；、、分别代表地层分界面两侧的纵波模量、纵横波速比和密度的变化量；

37、地震反射系数与地震数据的数学关系可以利用褶积公式来表征：

38、                         公式9

39、式中，表示子波矩阵；

40、将公式8与公式9相结合，即可建立纵波模量、纵横波速比和密度等弹性参数与地震记录的地震正演关系；

41、子步骤5.2：根据贝叶斯定理，后验概率分布可以写为先验分布与似然函数的乘积，即：

42、                               公式10

43、假设先验分布与似然函数均满足高斯分布，即：

44、           公式11

45、     公式12

46、因此，根据公式10、11与12，目标函数可以写为：

47、   公式13

48、式中，表示正演方程，即公式8与公式9；表示地震噪声的协方差矩阵，表示反演参数的协方差矩阵；表示先验平滑约束模型，通常通过测井数据插值与外推来获得；表示采样点数。

49、进一步地，步骤六中，利用迭代更新算法产生新的纵横波速比模型，具体方法为：

50、根据公式1与公式7，新模型可以表示为：

51、                              公式14

52、式中，为权重系数；为当前模型。对于次级变量纵横波速比，。对于主变量纵波模量与密度，。

53、进一步地，步骤九中，利用metropolis准则来判断是否保留新模型，包括以下子步骤：

54、子步骤9.1，结合反演似然函数与先验约束模型，建立如下目标函数：

55、  公式15

56、式中，表示正演方程，即公式8与公式9；表示地震噪声的协方差矩阵，表示反演参数的协方差矩阵；

57、子步骤9.2，为了避免反演目标值陷入局部极值，利用metropolis准则对产生的新模型进行优选，建立如下接受概率表达式：

58、               公式16

59、式中，表示反演目标函数，为可调参数。产生一个0到1之间的随机数，当接受概率大于随机数时，更新当前最优解，当接受概率小于随机数时，保留原来的最优解。

60、本发明的有益效果是：

61、将快速傅里叶滑动平均协模拟方法与迭代更新方法相结合，与传统的快速傅里叶滑动平均模拟相比，快速傅里叶滑动平均协模拟通过充分利用反演参数之间的相关性，提高了反演参数模拟的精度；而迭代更新方法可以提高反演迭代过程中的稳定性，有效提高反演的收敛性和精度。