CN104316961A - 获取风化层的地质参数的方法 - Google Patents

Abstract

本发明提供一种获取风化层的地质参数的方法，所述方法包括：(a)获取目标区域的标准的初至数据；(b)确定最小偏移距；(c)确定每个点的风化层速度值；(d)计算选择的m个点中的任意两点之间的风化层速度值的变差函数值；(e)计算选择的n个点中的任意两点之间的微测井资料的速度值的变差函数值；(f)计算m个点中的任意一点与n个点中的任意一点之间的关于速度的交叉编译函数；(g)获得与m个点中的每个点的风化层速度值对应的m个加权系数和与n个点中的每个点的微测井数据的速度值对应的n个加权系数；(h)计算风化层速度的估计值。采用本发明获取风化层的地质参数的方法可得到更符合地质分布规律的浅层速度结构。

Description

获取风化层的地质参数的方法

技术领域

本发明总体来说涉及石油地震资料分析处理领域，更具体地讲，涉及一种基于协同克里金技术融合微测井资料的获取风化层的地质参数的方法。

背景技术

在地震勘探中，近地表的低降速带(即，风化层)的速度和厚度不均匀变化会引起地震波传播旅行时的不均匀延迟，致使实际观测得到的反射波到达时间发生畸变，极大地影响着地震剖面的成像和构造形态。静校正就是消除地表低降速带对勘探资料的影响，对于复杂地表条件下地震资料的处理，求取正确的静校正量是非常关键和必要的。如果不能很好的解决严重的静校正问题，就会在地震剖面上产生假象或不成像，就会给地质解释带来严重困难，甚至错误的地质构造认识。

目前，在解决静校正问题过程中，对于存在稳定折射层的区域，一般采用折射静校正方法(例如，基于迭代法的折射静校正方法和广义互换静校正方法)进行计算。但这两类方法在使用时，对于风化层速度，一般参照工作经验，给定常速度值进行计算。然而，很多地区的风化层速度横向变化剧烈，实际风化层速度与前述参照工作经验给定的常速度值存在较大差距。

此外，有的方法还会对实际风化层速度进行拟合计算或者结合微测井资料进行速度反演，其效果优于给定常速度值计算的结果。但采用速度拟合的方法仍存在地震剖面的成像横向不连续现象，需要进一步处理；而进行速度反演，则计算时间长，微测井资料的好坏将直接影响速度反演的结果，不确定程度高，从而影响处理后叠加剖面成像效果，使得构造形态不准确，信噪比大大降低，不能达到资料精细处理的目的。

发明内容

本发明的目的是提出一种基于协同克里金技术融合微测井资料的获取风化层的地质参数的方法，通过变异函数这一数学工具，分别刻画出微测井资料的速度值、厚度值和由初至数据拟合出的风化层速度、厚度的变化趋势，再运用交叉编译函数得到两类数据的整体变化趋势，并加以它们的相关系数约束，最终得到符合地质分布规律的浅层速度结构。

本发明的一方面提供一种获取风化层的地质参数的方法，所述方法包括：(a)获取目标区域的标准的初至数据；(b)确定最小偏移距；(c)在所述目标区域内选择m个点，从获取的初至数据中选择m个点中的每个点在最小偏移距范围内的初至数据，并对选择的m个点中的每个点在最小偏移距范围内对应的初至数据进行拟合，确定每个点的风化层速度值，m为大于等于1的自然数；(d)根据确定的每个点的风化层速度值，计算选择的m个点中的任意两点之间的风化层速度值的变差函数值；(e)在所述目标区域内选择n个点，从所述目标区域的微测井资料中提取n个点中的每个点的微测井资料的速度值，并根据提取的每个点的微测井资料的速度值，计算选择的n个点中的任意两点之间的微测井资料的速度值的变差函数值，n为大于等于1的自然数；(f)根据确定的m个点中的每个点的风化层速度值和n个点中的每个点的微测井资料的速度值，计算m个点中的任意一点与n个点中的任意一点之间的关于速度的交叉编译函数；(g)根据步骤(d)～步骤(f)的结果获得与m个点中的每个点的风化层速度值对应的m个加权系数和与n个点中的每个点的微测井数据的速度值对应的n个加权系数；(h)根据获得的与m个点中的每个点的风化层速度值对应的m个加权系数、与n个点中的每个点的微测井数据的速度值对应的n个加权系数、确定的m个点中的每个点的风化层速度值和n个点中的每个点的微测井资料的速度值，计算风化层速度的估计值。

可选地，步骤(g)可包括：根据步骤(d)～步骤(f)的结果利用协同克里金方程组获得与m个点中的每个点的风化层速度值对应的m个加权系数和与n个点中的每个点的微测井数据的速度值对应的n个加权系数。

可选地，步骤(g)可包括：利用下面的协同克里金方程组计算与m个点中的每个点的风化层速度值对应的m个加权系数和与n个点中的每个点的微测井数据的速度值对应的n个加权系数，

$\left\{\begin{array}{c}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}{C}_x(x_i,x_p)+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}{C}_c(x_i,x_j)+{\mathrm{\μ}}_x=C_x(0,x_i),p\∈[1,n]\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}{C}_c(y_j,x_i)+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}{C}_y(y_q,y_j)+{\mathrm{\μ}}_y=C_y(0,y_j),q\∈[1,m]\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}=1\\ \underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}=0\end{array}\right.$

其中，C_x(x_i,x_p)为n个点中的第i点的坐标位置x_i与n个点中的第p点的坐标位置x_p对应的微测井资料的速度值的变差函数值，C_x(0,x_i)为待求位置0与n个点中的第i点的坐标位置x_i对应的微测井资料的速度值的变差函数值，C_y(y_q,y_j)为m个点中的第q点的坐标位置y_q与m个点中的第j点的坐标位置y_j对应的风化层速度值的变差函数值，C_y(0,y_j)为待求位置0与m个点中的第j点的坐标位置y_j对应的风化层速度值的变差函数值，C_c(x_i,y_j)和C_c(y_j,x_i)为关于n个点中的第i点的坐标位置x_i对应的微测井资料的速度值与m个点中的第j点的坐标位置y_j对应的风化层速度值的交叉编译函数，为协同克里金方程组中与n个点中的第i点的微测井数据的速度值对应的加权系数，为协同克里金方程组中与m个点中的第j点的风化层速度值对应的加权系数，μ_x、μ_y为计算方程组的解与原始值的误差的参数。

可选地，步骤(h)可包括：利用下面的公式计算风化层速度的估计值Z^*(0)，

$Z^{*}\left(0\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}{Z}_x\left(x_i\right)+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}{Z}_y\left(y_j\right)$

其中，Z^*(0)为风化层速度在待求位置0处的估计值，Z_x(x_i)为n个点中的第i点的坐标位置x_i对应的微测井资料的速度值，为与n个点中的第i点的微测井数据的速度值对应的加权系数，Z_y(y_j)为m个点中的第j点的坐标位置y_j对应的风化层速度值，为与m个点中的第j点的风化层速度值对应的加权系数。

可选地，步骤(c)可还包括：对选择的m个点中的每个点在最小偏移距范围内对应的初至数据进行拟合，确定每个点的炮点延迟时的初始值，其中，所述方法可还包括：(o)确定最大偏移距；(p)从获取的初至数据中选择m个点中的每个点在最小偏移距至最大偏移距范围内的初至数据，对选择的m个点中的每个点在最小偏移距至最大偏移距范围内对应的初至数据进行拟合，确定每个点的高速层速度的初始值；(q)根据确定的m个点中的每个点的高速层速度的初始值和m个点中的每个炮点延迟时的初始值，利用折射初至时间公式依次迭代计算m个点中的每个点的检波点延迟时、炮点延迟时和高速层速度。

可选地，步骤(q)可包括：(q1)根据m个点中的每个点的高速层速度、炮点延迟时、初至数据中每个点对应的炮检对从激发到接收产生的时间、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算m个点中的每个点的检波点延迟时，其中，将m个点中的每个点的高速层速度的初始值和m个点中的每个点的炮点延迟时的初始值作为初值代入；(q2)根据m个点中的每个点的高速层速度、步骤(q1)得到的m个点中的每个点的检波点延迟时、初至数据中每个点对应的炮检对从激发到接收产生的时间、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算m个点中的每个点的炮点延迟时；(q3)根据步骤(q1)得到的m个点中的每个点的检波点延迟时、步骤(q2)得到的m个点中的每个点的炮点延迟时、初至数据中每个点对应的炮检对从激发到接收产生的时间、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算m个点中的每个点的高速层速度；(q4)根据步骤(q1)得到的m个点中的每个点的检波点延迟时、步骤(q2)得到的m个点中的每个点的炮点延迟时、步骤(q3)得到的m个点中的每个点的高速层速度、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算对应的初至数据中每个点对应的炮检对从激发到接收产生的时间；(q5)当所述时间小于等于设定的时间范围时，将步骤(q1)～步骤(q3)的结果作为m个点中的每个点的炮点延迟时、检波点延迟时和高速层速度，当所述时间大于设定的时间范围时，则返回执行步骤(q1)。

可选地，折射初至时间公式可为，

$t_{{\mathrm{AB}}_j}=T_A\left(y_j\right)+\frac{{\stackrel{\‾}{\mathrm{AB}}}_j}{S\left(y_j\right)}+T_B\left(y_j\right)$

其中，t_ABj为初至数据中m个点中的第j点对应的炮检对从激发到接收产生的时间，T_A(y_j)为m个点中的第j点的坐标位置y_j对应的炮点延迟时，T_B(y_j)为m个点中的第j点的坐标位置y_j对应的检波点延迟时，为m个点中的第j点对应的炮点与检波点之间的距离，S(y_j)为m个点中的第j点的坐标位置y_j对应的高速层速度。

可选地，所述方法可还包括：(k)根据m个点中的每个点的风化层速度值、m个点中的每个点的高速层速度、m个点中的每个点的炮点延迟时或者检波点延迟时，计算m个点中的每个点的风化层厚度值；其中，步骤(d)可还包括：根据确定的每个点的风化层厚度值，计算选择的m个点中的任意两点之间的风化层厚度值的变差函数值；步骤(e)可还包括：从所述目标区域的微测井资料中提取n个点中的每个点的微测井资料的厚度值，并根据提取的每个点的微测井资料的厚度值，计算选择的n个点中的任意两点之间的微测井资料的厚度值的变差函数值；步骤(f)可还包括：根据确定的m个点中的每个点的风化层厚度值和n个点中的每个点的微测井资料的厚度值，计算m个点中的任意一点与n个点中的任意一点之间的关于厚度的交叉编译函数；步骤(g)可还包括：根据步骤(d)～步骤(f)的结果获得与m个点中的每个点的风化层厚度值对应的m个加权系数和与n个点中的每个点的微测井数据的厚度值对应的n个加权系数；步骤(h)可还包括：根据获得的与m个点中的每个点的风化层厚度值对应的m个加权系数、与n个点中的每个点的微测井数据的厚度值对应的n个加权系数、确定的m个点中的每个点的风化层厚度值和n个点中的每个点的微测井资料的厚度值，计算风化层厚度的估计值。

可选地，步骤(k)可包括：利用下面的公式计算m个点中的每个点的风化层厚度值，

$D_y\left(y_j\right)=T\left(y_j\right)\·\frac{Z_y\left(y_j\right)\·S\left(y_j\right)}{\sqrt{S^2\left(y_j\right)-Z_y^2\left(y_j\right)}}$

其中，D_y(y_j)为m个点中的第j点的坐标位置y_j对应的风化层厚度值，T(y_j)为m个点中的第j点的坐标位置y_j对应的炮点延迟时或检波点延迟时，Z_y(y_j)为m个点中的第j点的坐标位置y_j对应的风化层速度值，S(y_j)为m个点中的第j点的坐标位置y_j对应的高速层速度。

可选地，所述方法可还包括：(m)根据风化层速度的估计值和风化层厚度的估计值，计算折射静校正量。

采用本发明所述的获取风化层的地质参数的方法是在拟合出风化层速度的基础上，进一步融合高精度的微测井资料的解释成果，得到优化后的风化层速度场，由于考虑了微测井资料的解释成果，则此速度场的速度分布更准确，更接近真实的速度分布。此外，在上述基础上还进行了迭代计算，得到最终的风化层速度场与高速层速度，从而求得更加准确的折射静校正量。

附图说明

通过下面结合附图进行的详细描述，本发明的上述和其它目的、特点和优点将会变得更加清楚，其中：

图1是示出根据本发明的示例性实施例的获取风化层的地质参数的方法的流程图；

图2是示出根据本发明的示例性实施例的计算任意距离对应的变差函数值的示意图。

具体实施方式

下面，将参照附图详细描述本发明的示例性实施例。

提供参照附图的以下描述以帮助对由权利要求及其等同物限定的本发明的示例性实施例的全面理解。包括各种特定细节以帮助理解，但这些细节仅被视为是示例性的。因此，本领域的普通技术人员将认识到在不脱离本发明的范围和精神的情况下，可对描述于此的示例性实施例进行各种改变和修改。此外，为了清楚和简洁，省略对公知的功能和结构的描述。

图1是示出根据本发明的示例性实施例的获取风化层的地质参数的方法的流程图。

在步骤S10中，获取目标区域的标准的初至数据。这里，标准的初至数据可由第三方初至拾取软件进行拾取获得。

在步骤S20中，确定最小偏移距和最大偏移距。这里，在目标区域的标准的初至数据中具有稳定的折射段的最小距离就是最小偏移距，在目标区域的标准的初至数据中具有稳定的折射段的最大距离就是最大偏移距。稳定的折射段一般出现在表层速度较低，在低速层下方具有稳定高速层的地区，对在这种地区获得的单炮资料进行拾取获得的初至数据，可以看出具有折射段现象。

这里，应该理解，每个工区(即，目标区域)的最小偏移距、最大偏移距是不同的。作为示例，对在不同的目标区域获得的初至数据，可通过观察以及在初至统计显示图上进行测量来确定获得的不同目标区域的初至数据的最小偏移距和最大偏移距。

在步骤S30中，在所述目标区域内选择m个点，从获取的初至数据中选择m个点中的每个点在最小偏移距范围内的初至数据，并对选择的m个点中的每个点在最小偏移距范围内对应的初至数据进行拟合，确定每个点的风化层速度值和每个点的炮点延迟时的初始值。这里，m可为大于等于1的自然数。

具体说来，在目标区域内选择m个炮点，每炮对应多道的初至数据，从每炮的所有初至数据中选择在最小偏移距范围内的道对应的初至数据，从选择的初至数据中获取每道对应的单炮偏移距H_b和每道拾取的初至时间G_b。这里，b表示道号，1≤b≤B，B为大于等于1的自然数。然后在横坐标为时间，纵坐标为偏移距的二维坐标系中，根据每道对应的单炮偏移距H_b和每道拾取的初至时间G_b利用线性拟合的方法得到一条最优的直线。此时，该最优直线的斜率则对应m个点中的每个点的风化层速度值的倒数，该最优直线的截距则对应m个点中的每个点的炮点延迟时的初始值。这里，利用线性拟合的方法得到最优直线的步骤为本领域的公知常识，本发明对此部分的内容不再赘述。

在步骤S40中，从获取的初至数据中选择m个点中的每个点在最小偏移距至最大偏移距范围内的初至数据，对选择的m个点中的每个点在最小偏移距至最大偏移距范围内对应的初至数据进行拟合，确定每个点的高速层速度的初始值。

具体说来，类似于步骤S30，从每炮的所有初至数据中选择在最小偏移距至最大偏移距范围内的道对应的初至数据，从选择的初至数据中获取每道对应的单炮偏移距和每道拾取的初至时间；然后在横坐标为时间，纵坐标为偏移距的二维坐标系中，根据每道对应的单炮偏移距和每道拾取的初至时间利用线性拟合的方法得到一条最优的直线。此时，该最优直线的斜率则对应于m个点中的每个点的高速层速度的初始值的倒数。

在步骤S50中，根据确定的m个点中的每个点的高速层速度的初始值和m个点中的每个点的炮点延迟时的初始值，利用折射初至时间公式依次迭代计算m个点中的每个点的检波点延迟时、炮点延迟时和高速层速度。

可选地，根据确定的m个点中的每个点的高速层速度的初始值和m个点中的每个点的炮点延迟时的初始值，利用折射初至时间公式依次迭代计算m个点中的每个点的检波点延迟时、炮点延迟时和高速层速度的步骤可包括：

(a)根据m个点中的每个点的高速层速度、炮点延迟时、初至数据中每个点对应的炮检对从激发到接收产生的时间、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算m个点中的每个点的检波点延迟时。这里，可将m个点中的每个点的高速层速度的初始值和m个点中的每个点的炮点延迟时的初始值作为初值代入。

(b)根据m个点中的每个点的高速层速度、步骤(a)得到的m个点中的每个点的检波点延迟时、初至数据中每个点对应的炮检对从激发到接收产生的时间、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算m个点中的每个点的炮点延迟时。

(c)根据步骤(a)得到的m个点中的每个点的检波点延迟时、步骤(b)得到的m个点中的每个点的炮点延迟时、初至数据中每个点对应的炮检对从激发到接收产生的时间、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算m个点中的每个点的高速层速度。

(d)根据步骤(a)得到的m个点中的每个点的检波点延迟时、步骤(b)得到的m个点中的每个点的炮点延迟时、步骤(c)得到的m个点中的每个点的高速层速度、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算对应的初至数据中每个点对应的炮检对从激发到接收产生的时间。

(e)当所述时间小于等于设定的时间范围时，将步骤(a)～步骤(c)的结果作为m个点中的每个点的炮点延迟时、检波点延迟时和高速层速度，当所述时间大于设定的时间范围时，则返回执行步骤(a)。

这里，应该理解，可当所述时间小于等于设定的时间范围时，保留步骤(a)～步骤(c)此时计算得到的结果。然而本发明不限于此，还可以预先设定步骤(a)～步骤(c)的循环次数，当所述时间大于设定的时间范围时，将当前的循环次数加一，并判断当前的循环次数是否达到预先设定的循环次数，如果没达到，则返回执行步骤(a)，如果达到，则保留步骤(a)～步骤(c)此时计算得到的结果。

可选地，折射初至时间公式可为，

$t_{{\mathrm{AB}}_j}=T_A\left(y_j\right)+\frac{{\stackrel{\‾}{\mathrm{AB}}}_j}{S\left(y_j\right)}+T_B\left(y_j\right)---\left(1\right)$

公式(1)中，t_ABj为初至数据中m个点中的第j点对应的炮检对从激发到接收产生的时间，T_A(y_j)为m个点中的第j点的坐标位置y_j对应的炮点延迟时，T_B(y_j)为m个点中的第j点的坐标位置y_j对应的检波点延迟时，为m个点中的第j点对应的炮点与检波点之间的距离，S(y_j)为m个点中的第j点的坐标位置y_j对应的高速层速度。

具体说来，可将步骤S30中得到的m个点中的第j点的炮点延迟时的初始值和步骤S40得到的m个点中的第j点的高速层速度的初始值作为初值代入到公式(1)中，由于t_ABj和可从初至数据中获得，因此，可由公式(1)求得m个点中的第j点的检波点延迟时T_B(y_j)；再将高速层速度的初始值、T_B(y_j)、t_ABj和代入公式(1)中，可求得炮点延迟时T_A(y_j)；再将T_B(y_j)、T_A(y_j)、t_ABj和代入公式(1)中，可求得高速层速度S(y_j)；然后再将T_A(y_j)、T_B(y_j)、S(y_j)和代入公式(1)中，可求得时间t_ABj；最后，对时间t_ABj进行判断，当此时的时间t_ABj小于等于设定的时间范围时，则保留当前计算得到的T_A(y_j)、T_B(y_j)、S(y_j)，当此时的时间t_ABj大于设定的时间范围时，则利用上述的步骤再重新求取新的炮点延迟时、检波点延迟时和高速层速度，直至时间t_ABj小于等于设定的时间范围。

在步骤S60中，根据m个点中的每个点的风化层速度值、m个点中的每个点的高速层速度、m个点中的每个点的炮点延迟时或者检波点延迟时，计算m个点中的每个点的风化层厚度值。

可选地，可利用下面的公式计算m个点中的每个点的风化层厚度值，

$D_y\left(y_j\right)=T\left(y_j\right)\·\frac{Z_y\left(y_j\right)\·S\left(y_j\right)}{\sqrt{S^2\left(y_j\right)-Z_y^2\left(y_j\right)}}---\left(2\right)$

公式(2)中，D_y(y_j)为m个点中的第j点的坐标位置y_j对应的风化层厚度值，T(y_j)为m个点中的第j点的坐标位置y_j对应的炮点延迟时或检波点的延迟时，Z_y(y_j)为m个点中的第j点的坐标位置y_j对应的风化层速度值，S(y_j)为m个点中的第j点的坐标位置y_j对应的高速层速度。

在步骤S70中，根据确定的每个点的风化层速度值或厚度值，计算选择的m个点中的任意两点之间的风化层速度值或厚度值的变差函数值。

可选地，可利用下面的公式计算风化层速度值的变差函数值，

${\mathrm{\γ}}_y^{*}\left(h\right)=\frac{1}{2N\left(h\right)}\underset{j=1}{\overset{N\left(h\right)}{\mathrm{\Σ}}}{[Z_y\left(y_j\right)-Z_y(y_j+h)]}^2---\left(3\right)$

公式(3)中，表示风化层速度值的变差函数值，h为两点之间的距离，N(h)为m个点中满足两点之间的距离为h的点对数目，Z_y(y_j)为m个点中的第j点的坐标位置y_j对应的风化层速度值，1≤j≤N(h)，Z_y(y_j+h)为与第j点的坐标位置y_j的距离为h的坐标位置y_j+h对应的风化层速度值。

同样地，可将确定的m个点中的每个点的风化层厚度值替代m个点中的每个点的风化层速度值代入到上述公式(3)中，则可计算得到选择的m个点中的任意两点之间的风化层厚度值的变差函数值。

在步骤S80中，在所述目标区域内选择n个点，从所述目标区域的微测井资料中提取n个点中的每个点的微测井资料的速度值或厚度值，并根据提取的每个点的微测井资料的速度值或厚度值，计算选择的n个点中的任意两点之间的微测井资料的速度值或厚度值的变差函数值，n为大于等于1的自然数。

同样地，可将确定的n个点中的每个点的微测井资料的速度值或厚度值替代m个点中的每个点的风化层速度值代入到上述公式(3)，则可计算得到选择的n个点中的任意两点之间的微测井资料的速度值或厚度值的变差函数值。

在步骤S90中，根据确定的m个点中的每个点的风化层速度值和n个点中的每个点的微测井资料的速度值，计算m个点中的任意一点与n个点中的任意一点之间的关于速度的交叉编译函数。

可选地，计算交叉编译函数的公式可为，

${\mathrm{\γ}}_{\mathrm{xy}}\left(h\right)=\frac{1}{2}E[Z_x(s+h)Z_y(s+h)-Z_x(s+h)Z_y\left(s\right)-Z_x\left(s\right)Z_y(s+h)+Z_x\left(s\right)Z_y\left(s\right)]---\left(4\right)$

公式(4)中，γ_xy(h)表示交叉编译函数，h为两点之间的距离，E为数学期望，Z_x(s)为任意坐标位置s对应的微测井资料的速度值，Z_x(s+h)为与任意坐标位置s的距离为h的坐标位置s+h对应的微测井资料的速度值，Z_y(s)为任意坐标位置s对应的风化层速度值，Z_y(s+h)为与任意坐标位置s的距离为h的坐标位置s+h对应的风化层速度值。

可选地，本发明所述的方法可还包括：分别设置n个点的微测井资料的速度值Z_x(x_i)和m个点的风化层速度值Z_y(y_j)的相关系数。这里，一般来讲，相关系数用于描述微测井资料的速度值、风化层速度值对协同克里金方程组的贡献程度，如果微测井资料的速度值的贡献程度高，则相关系数大，如果风化层速度值的贡献程度高，则相关系数小。相关系数的范围在0-1之间，极端的若相关系数为0，则协同克里金方程组对应的计算结果与仅利用风化层速度值进行计算得到的结果相等，若相关系数为1，则协同克里金方程组对应的计算结果仅利用微测井资料的速度值进行计算得到的结果相等。例如，相关系数的值一般可手动设置，大小可在0.6-0.8之间。

在步骤S100中，根据步骤S70～步骤S90的结果获得与m个点中的每个点的风化层速度值或厚度值对应的m个加权系数和与n个点中的每个点的微测井数据的速度值或厚度值对应的n个加权系数。

可选地，可根据步骤S70～步骤S90的结果利用协同克里金方程组获得与m个点中的每个点的风化层速度值对应的m个加权系数和与n个点中的每个点的微测井数据的速度值对应的n个加权系数。

可选地，可利用下面的协同克里金方程组计算与m个点中的每个点的风化层速度值对应的m个加权系数和与n个点中的每个点的微测井数据的速度值对应的n个加权系数，

$\left\{\begin{array}{c}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}{C}_x(x_i,x_p)+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}{C}_c(x_i,x_j)+{\mathrm{\μ}}_x=C_x(0,x_i),p\∈[1,n]\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}{C}_c(y_j,x_i)+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}{C}_y(y_q,y_j)+{\mathrm{\μ}}_y=C_y(0,y_j),q\∈[1,m]\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}=1\\ \underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}=0\end{array}\right.---\left(5\right)$

公式(5)中，C_x(x_i,x_p)为n个点中的第i点的坐标位置x_i与n个点中的第p点的坐标位置x_p对应的微测井资料的速度值的变差函数值，C_x(0,x_i)为待求位置0与n个点中的第i点的坐标位置x_i对应的微测井资料的速度值的变差函数值，C_y(y_q,y_j)为m个点中的第q点的坐标位置y_q与m个点中的第j点的坐标位置y_j对应的风化层速度值的变差函数值，C_y(0,y_j)为待求位置0与m个点中的第j点的坐标位置y_j对应的风化层速度值的变差函数值，C_c(x_i,y_j)和C_c(y_j,x_i)为关于n个点中的第i点的坐标位置x_i对应的微测井资料的速度值与m个点中的第j点的坐标位置y_j对应的风化层速度值的交叉编译函数，为协同克里金方程组中与n个点中的第i点的微测井数据的速度值对应的加权系数，为协同克里金方程组中与m个点中的第j点的风化层速度值对应的加权系数，μ_x、μ_y为计算方程组的解与原始值的误差的参数。这里，加权系数 需满足无偏估计前提条件，即，无偏估计就是系统误差为零的估计，是协同克里金方程组的前提假设。

具体说来，在公式(5)的方程组中，仅加权系数为未知量，其他所有的参数均为已知量。这里，C_x(x_i,x_p)和C_y(y_q,y_j)可由公式(3)计算得到，C_c(x_i,y_j)和C_c(y_j,x_i)可由公式(4)计算得到，C_x(0,x_i)和C_y(0,y_j)则可根据m个点中的每个点的风化层速度值和n个点中的每个点的微测井资料的速度值，利用拟合算法得到。

具体说来，利用拟合算法计算待求位置0与m个点中的第j点的坐标位置y_j对应的风化层速度值的变差函数值C_y(0,y_j)的步骤可为：根据m个点的每个点的风化层速度值，利用公式(3)可求得m个点中的任意两点之间的风化层速度值的变差函数值；然后利用拟合算法(例如，多元线性回归法)求出关于风化层速度值的变差函数值的拟合曲线；最后基于所述拟合曲线就可得到任意距离对应的变差函数值。这里，利用拟合算法求出拟合曲线的步骤为本领域的公知常识，本发明对此部分的内容不再赘述。

图2是示出根据本发明的示例性实施例的计算任意距离对应的变差函数值的示意图。

参照图2，横坐标表示距离h，纵坐标表示变差函数值。如图2所示，利用公式(3)可求得已知的m个点中的任意两点之间的风化层速度值的变差函数值(例如，)，如图2中所示的散点，然后再利用多元线性回归方法得到图中所示的拟合曲线1，这样就可基于该拟合曲线1得到任意距离h对应的关于风化层速度值的变差函数值。

例如，协同克里金方程组(5)中的C_y(0,y_j)为待求位置0与m个点中的第j点的坐标位置y_j对应的风化层速度值的变差函数值，由于第j点的坐标位置y_j和待求位置0的坐标位置都是已知的，也就是说，待求位置0与坐标位置y_j之间的距离是已知，则基于变差函数值的拟合曲线就可得到C_y(0,y_j)的值。

同样地，可利用上述方法计算待求位置0与n个点中的第i点的坐标位置x_i对应的微测井资料的速度值的变差函数值C_x(0,x_i)。同样地，还可利用上述方法计算任意距离对应的风化层厚度值的变差函数值、微测井资料的厚度值的变差函数值。

在步骤S110中，根据获得的与m个点中的每个点的风化层速度或厚度值对应的m个加权系数、与n个点中的每个点的微测井数据的速度值或厚度值对应的n个加权系数、确定的m个点中的每个点的风化层速度值或厚度值和n个点中的每个点的微测井资料的速度值或厚度值，计算风化层速度或厚度的估计值。

可选地，可利用下面的公式计算风化层速度的估计值，

$Z^{*}\left(0\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}{Z}_x\left(x_i\right)+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}{Z}_y\left(y_j\right)---\left(6\right)$

公式(6)中，Z^*(0)为风化层速度在待求位置0处的估计值，Z_x(x_i)为n个点中的第i点的坐标位置x_i对应的微测井资料的速度值，为与n个点中的第i点的微测井数据的速度值对应的加权系数，Z_y(y_j)为m个点中的第j点的坐标位置y_j对应的风化层速度值，为与m个点中的第j点的风化层速度值对应的加权系数。

同样地，可将确定的m个点中的每个点的风化层厚度值、n个点中的每个点的微测井数据的厚度值、与n个点中的第i点的微测井数据的厚度值对应的加权系数、与m个点中的第j点的风化层厚度值对应的加权系数代入到上述公式中，则可计算得到风化层厚度的估计值。

在步骤S120中，根据风化层速度的估计值和风化层厚度的估计值，计算折射静校正量。

可选地，可利用下面的公式计算折射静校正量

公式(7)中，为在待求位置0处对应的炮点或检波点折射静校正量，τ为井深或检波器初至时间，k为低降速层的层数，1≤k≤f，f为大于等于1的自然数，D₀为给定的厚度，Z^*(0)为风化层速度在待求位置0处的估计值，为风化层厚度在待求位置0处的估计值，H_d为基准面高程，H_g为高速层顶界面高程，V_R为基准面校正速度。

具体说来，在利用公式(2)计算风化层厚度值时，如果T(y_j)代入的为m个点中的第j点的坐标位置y_j对应的炮点延迟时，则由公式(7)得到的就是在待求位置0处对应的炮点折射静校正量，同样地，如果T(y_j)代入的为m个点中的第j点的坐标位置y_j对应的检波点延迟时，则由公式(7)得到的就是在待求位置0处对应的检波点折射静校正量。

这里，在现有技术中很难对低降速层进行分层，因此，在大部分计算中，一般都默认低降速层的层数为1层。

本发明所述的方法是在拟合出风化层速度的基础上，进一步融合高精度的微测井资料的解释成果，得到优化后的风化层速度场，由于考虑了微测井资料的解释成果，则此速度场的速度分布更准确，更接近真实的速度分布。

此外，在上述基础上还进行了迭代计算，得到最终的风化层速度场与高速层速度，从而求得更加准确的折射静校正量。本发明所述的方法能大大改善剖面的成像质量，本发明所述的方法技术操作方便、运行速度快，能够满足实际生产的需求。

本发明提供的基于协同克里金技术融合微测井资料的获取风化层的地质参数的方法主要通过变异函数这一数学工具，分别刻画出微测井资料的速度值、厚度值和由初至数据拟合出的风化层速度、厚度的变化趋势，再运用交叉编译函数得到两类数据的整体变化趋势，并加以它们的相关系数约束，最终得到符合地质分布规律的浅层速度结构，以提高折射静校正的计算精度，改善剖面的成像质量。

尽管已经参照其示例性实施例具体显示和描述了本发明，但是本领域的技术人员应该理解，在不脱离权利要求所限定的本发明的精神和范围的情况下，可以对其进行形式和细节上的各种改变。

Claims (10)

1.一种获取风化层的地质参数的方法，所述方法包括：

(a)获取目标区域的标准的初至数据；

(b)确定最小偏移距；

(c)在所述目标区域内选择m个点，从获取的初至数据中选择m个点中的每个点在最小偏移距范围内的初至数据，并对选择的m个点中的每个点在最小偏移距范围内对应的初至数据进行拟合，确定每个点的风化层速度值，m为大于等于1的自然数；

(d)根据确定的每个点的风化层速度值，计算选择的m个点中的任意两点之间的风化层速度值的变差函数值；

(e)在所述目标区域内选择n个点，从所述目标区域的微测井资料中提取n个点中的每个点的微测井资料的速度值，并根据提取的每个点的微测井资料的速度值，计算选择的n个点中的任意两点之间的微测井资料的速度值的变差函数值，n为大于等于1的自然数；

(f)根据确定的m个点中的每个点的风化层速度值和n个点中的每个点的微测井资料的速度值，计算m个点中的任意一点与n个点中的任意一点之间的关于速度的交叉编译函数；

(g)根据步骤(d)～步骤(f)的结果获得与m个点中的每个点的风化层速度值对应的m个加权系数和与n个点中的每个点的微测井数据的速度值对应的n个加权系数；

(h)根据获得的与m个点中的每个点的风化层速度值对应的m个加权系数、与n个点中的每个点的微测井数据的速度值对应的n个加权系数、确定的m个点中的每个点的风化层速度值和n个点中的每个点的微测井资料的速度值，计算风化层速度的估计值。

2.根据权利要求1所述的方法，其中，步骤(g)包括：根据步骤(d)～步骤(f)的结果利用协同克里金方程组获得与m个点中的每个点的风化层速度值对应的m个加权系数和与n个点中的每个点的微测井数据的速度值对应的n个加权系数。

3.根据权利要求2所述的方法，其中，步骤(g)包括：

利用下面的协同克里金方程组计算与m个点中的每个点的风化层速度值对应的m个加权系数和与n个点中的每个点的微测井数据的速度值对应的n个加权系数，

$\left\{\begin{array}{c}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}{C}_x(x_i,x_p)+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}{C}_c(x_i,y_j)+{\mathrm{\μ}}_x=C_x(0,x_i),p\∈[1,n]\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}{C}_c(y_j,x_i)+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}{C}_y(y_q,y_j)+{\mathrm{\μ}}_y=C_y(0,y_j),q\∈[1,m]\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}=1\\ \underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}=0\end{array}\right.$

其中，C_x(x_i,x_p)为n个点中的第i点的坐标位置x_i与n个点中的第p点的坐标位置x_p对应的微测井资料的速度值的变差函数值，C_x(0,x_i)为待求位置0与n个点中的第i点的坐标位置x_i对应的微测井资料的速度值的变差函数值，C_y(y_q,y_j)为m个点中的第q点的坐标位置y_q与m个点中的第j点的坐标位置y_j对应的风化层速度值的变差函数值，C_y(0,y_j)为待求位置0与m个点中的第j点的坐标位置y_j对应的风化层速度值的变差函数值，C_c(x_i,y_j)和C_c(y_j,x_i)为关于n个点中的第i点的坐标位置x_i对应的微测井资料的速度值与m个点中的第j点的坐标位置y_j对应的风化层速度值的交叉编译函数，为协同克里金方程组中与n个点中的第i点的微测井数据的速度值对应的加权系数，为协同克里金方程组中与m个点中的第j点的风化层速度值对应的加权系数，μ_x、μ_y为计算方程组的解与原始值的误差的参数。

4.根据权利要求1所述的方法，其中，步骤(h)包括：

利用下面的公式计算风化层速度的估计值Z^*(0)，

$Z^{*}\left(0\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_{x_i}{Z}_x\left(x_i\right)+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_{y_j}{Z}_y\left(y_j\right)$

其中，Z^*(0)为风化层速度在待求位置0处的估计值，Z_x(x_i)为n个点中的第i点的坐标位置x_i对应的微测井资料的速度值，为与n个点中的第i点的微测井数据的速度值对应的加权系数，Z_y(y_j)为m个点中的第j点的坐标位置y_j对应的风化层速度值，为与m个点中的第j点的风化层速度值对应的加权系数。

5.根据权利要求1所述的方法，其中，步骤(c)还包括：对选择的m个点中的每个点在最小偏移距范围内对应的初至数据进行拟合，确定每个点的炮点延迟时的初始值，

其中，所述方法还包括：

(o)确定最大偏移距；

(p)从获取的初至数据中选择m个点中的每个点在最小偏移距至最大偏移距范围内的初至数据，对选择的m个点中的每个点在最小偏移距至最大偏移距范围内对应的初至数据进行拟合，确定每个点的高速层速度的初始值；

(q)根据确定的m个点中的每个点的高速层速度的初始值和m个点中的每个点的炮点延迟时的初始值，利用折射初至时间公式计算m个点中的每个点的检波点延迟时、炮点延迟时和高速层速度。

6.根据权利要求5所述的方法，其中，步骤(q)包括：

(q1)根据m个点中的每个点的高速层速度、炮点延迟时、初至数据中每个点对应的炮检对从激发到接收产生的时间、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算m个点中的每个点的检波点延迟时，其中，将m个点中的每个点的高速层速度的初始值和m个点中的每个点的炮点延迟时的初始值作为初值代入；

(q2)根据m个点中的每个点的高速层速度、步骤(q1)得到的m个点中的每个点的检波点延迟时、初至数据中每个点对应的炮检对从激发到接收产生的时间、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算m个点中的每个点的炮点延迟时；

(q3)根据步骤(q1)得到的m个点中的每个点的检波点延迟时、步骤(q2)得到的m个点中的每个点的炮点延迟时、初至数据中每个点对应的炮检对从激发到接收产生的时间、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算m个点中的每个点的高速层速度；

(q4)根据步骤(q1)得到的m个点中的每个点的检波点延迟时、步骤(q2)得到的m个点中的每个点的炮点延迟时、步骤(q3)得到的m个点中的每个点的高速层速度、每个点对应的炮点与检波点之间的距离，利用折射初至时间公式计算对应的初至数据中每个点对应的炮检对从激发到接收产生的时间；

(q5)当所述时间小于等于设定的时间范围时，将步骤(q1)～步骤(q3)的结果作为m个点中的每个点的炮点延迟时、检波点延迟时和高速层速度，当所述时间大于设定的时间范围时，则返回执行步骤(q1)。

7.根据权利要求6所述的方法，其中，折射初至时间公式为，

$t_{{\mathrm{AB}}_j}=T_A\left(y_j\right)+\frac{{\stackrel{\‾}{\mathrm{AB}}}_j}{S\left(y_j\right)}+T_B\left(y_j\right)$

其中，t_ABj为初至数据中m个点中的第j点对应的炮检对从激发到接收产生的时间，T_A(y_j)为m个点中的第j点的坐标位置y_j对应的炮点延迟时，T_B(y_j)为m个点中的第j点的坐标位置y_j对应的检波点延迟时，为m个点中的第j点对应的炮点与检波点之间的距离，S(y_j)为m个点中的第j点的坐标位置y_j对应的高速层速度。

8.根据权利要求5所述的方法，所述方法还包括：

(k)根据m个点中的每个点的风化层速度值、m个点中的每个点的高速层速度、m个点中的每个点的炮点延迟时或者检波点延迟时，计算m个点中的每个点的风化层厚度值；

其中，步骤(d)还包括：根据确定的每个点的风化层厚度值，计算选择的m个点中的任意两点之间的风化层厚度值的变差函数值；

步骤(e)还包括：从所述目标区域的微测井资料中提取n个点中的每个点的微测井资料的厚度值，并根据提取的每个点的微测井资料的厚度值，计算选择的n个点中的任意两点之间的微测井资料的厚度值的变差函数值；

步骤(f)还包括：根据确定的m个点中的每个点的风化层厚度值和n个点中的每个点的微测井资料的厚度值，计算m个点中的任意一点与n个点中的任意一点之间的关于厚度的交叉编译函数；

步骤(g)还包括：根据步骤(d)～步骤(f)的结果获得与m个点中的每个点的风化层厚度值对应的m个加权系数和与n个点中的每个点的微测井数据的厚度值对应的n个加权系数；

步骤(h)还包括：根据获得的与m个点中的每个点的风化层厚度值对应的m个加权系数、与n个点中的每个点的微测井数据的厚度值对应的n个加权系数、确定的m个点中的每个点的风化层厚度值和n个点中的每个点的微测井资料的厚度值，计算风化层厚度的估计值。

9.根据权利要求8所述的方法，其中，步骤(k)包括：

利用下面的公式计算m个点中的每个点的风化层厚度值，

$D_y\left(y_j\right)=T\left(y_j\right)\·\frac{Z_y\left(y_j\right)\·S\left(y_j\right)}{\sqrt{S^2\left(y_j\right)-Z_y^2\left(y_j\right)}}$

其中，D_y(y_j)为m个点中的第j点的坐标位置y_j对应的风化层厚度值，T(y_j)为m个点中的第j点的坐标位置y_j对应的炮点延迟时或检波点延迟时，Z_y(y_j)为m个点中的第j点的坐标位置y_j对应的风化层速度值，S(y_j)为m个点中的第j点的坐标位置y_j对应的高速层速度。

10.根据权利要求8所述的方法，所述方法还包括：(m)根据风化层速度的估计值和风化层厚度的估计值，计算折射静校正量。