CN102466816B - 一种叠前地震数据地层弹性常数参数反演的方法 - Google Patents

Abstract

本发明涉及石油物探叠前地震数据地层弹性常数参数反演的方法，采集地震数据，对叠前地震数据进行处理，取得测井数据和提取角道集地震数据，形成全部叠加数据，将地震数据和测井数据层位标定，将不同入射角叠加的地震数据，反演得到波阻抗参数，计算对应的反射系数，形成测井与地震数据对，计算函数映射网络模型的权函数和模型参数、阻抗参数、弹性常数，绘制弹性常数剖面，用于储层岩性识别、油气预测、油水界面确定和油气藏的描述。本发明可利用常规地震数据和测井数据反演计算出弹性常数，对断层、尖灭带的反演有一定的适应能力，适应范围大、分辨率高、计算速度快、稳定性好、计算精度高、具有一定的抗噪能力的特点。

Description

一种叠前地震数据地层弹性常数参数反演的方法

技术领域

本发明涉及物探技术，属于油田的勘探、开发、开采过程中为储层预测、油水界面识别和油藏描述提供的一种叠前地震数据地层弹性常数参数反演的方法。

背景技术

地震勘探的过程，就是在地面上的一系列点上，利用人工激发地震波，地震波向地下传播，当遇到波阻抗(地震波在地层介质中向地下传播的速度与介质密度的乘积)界面(即上下地层波阻抗不相等面)时，在波阻抗界面上地震波产生反射，地震波传播方向发生改变，地震波开始向上传播，在地面上的一系列接收点上安置着接收器，接收向上传播的地震波数据，这是地震勘探的正过程(野外勘探过程)。而实际地面接收器接收到地震波数据中不但包含着地下地层波阻抗界面的信息，而且还包含着激发点和接收点空间位置和排列位置的信息，以及各种噪声干扰等。地震数据处理就是将野外勘探过程中接收的向上传播的地震波数据，经过处理，仅仅保留下反映地下地层波阻抗界面的信息，而消除其它的信息，这样得出的地震数据，就是叠后地震数据。地震波传播速度的大小反映了地下地层的结构和构造，而弹性常数反映了储层中油气分布情况以及油和水的分界面，是油气勘探、开发、开采过程中最重要的参数。弹性常数反演就是根据地面接收到的反映地下地层界面的反射信号，求取地下地层弹性常数的过程。

现行的地震技术可以求出相应的速度参数，但是由于地震记录中存在噪声，加上地震记录的分辨率较低，难于实现地质任务对地层纵横向分辨率的要求；测井技术虽然具有很高的纵向分辨率，但不具备横向分辨率，难于实现地质任务对井间地层参数变化的要求。

地震波速度反演方法要求地震波具有振幅真值，且激发地震波的震源为已知，保存有全套多次波信息，没有噪声干扰特别是没有决定性过程的规则干扰。在这样严格的条件下才能有效反演可靠的速度参数，但是这在实际地震数据采集过程中是无法保证的。因此实际地震波速度参数反演的根本问题是用于参数反演的已知信息严重不足，从而导致反演过程的失败。实际地震可以提供用于反演的数据，仅仅是叠后地震数据。叠后地震数据中，并不包含噪声信息，而仅仅包含着反映地下地层波阻抗界面(即速度)的信息和地震波传播的地震子波信息。仅仅已知叠后地震数据，要求取地下地层波阻抗界面的信息和地震波传播的地震子波信息，即用一个数据求两个未知数，因此求解需要的已知信息严重不足。通过对目标地区进行地质勘探，已知目标地区的地质层位信息，通过对叠后地震数据进行解释，将地质层位信息和地震层位信息对应起来，这些先验地质层位信息可以作为附加信息参加反演过程，这才是成功实现地震波速度参数反演的保证。

测井就是在已经钻探好的井中，将各种激发设备和接收设备放置在一根钢管内，激发设备和接收设备保持几米的距离远，用电缆连接这些激发设备和接收设备到地面。在地面上，通过匀速下降电缆和匀速提升电缆，同时利用激发设备进行信息激发，这些激发的信息穿过井壁地层，由接收设备接收，通过电缆传输到地面存储设备存储起来。一次可以同时得到许多测井参数。地震反演需要的测井参数是声波时差参数和密度参数。声波时差的倒数就是速度，因此利用这两个参数可以得出波阻抗和速度参数。声波时差参数是由声波测井测得的，密度参数是由密度测井测得的。声波测井设备包括一个声波脉冲发射器和一个声波脉冲接收器。由声波脉冲发射器发出的声波射向井壁，声波在地层中传播。声波脉冲接收器接收声波在地层中的传播。声波测井参数就是声波通过1英尺地层所需要的传播时间随深度变化的记录。密度测井设备包括一个屏蔽体内的放射性源和一个放射性探测器。由放射性源发出的伽马射线射向井壁地层，这些伽马射线可以看作为高速粒子，在地层中与电子碰撞。每次碰撞伽马射线传递能量给电子而失去一部分能量，能量减少后的伽马射线继续前进。放射性探测器接收到能量减少后的伽马射线。伽马射线能量的减少作为地层密度的指示记录下来，即密度参数。

测井参数数据采集很密集，采样率很小，因此测井参数的视分辨率很高，所得到的波阻抗和速度曲线中具有丰富的高频成分。但是这些高频成分既不对应反射界面，也不对应地层界面。因此利用测井参数得到的波阻抗和速度曲线，必须根据地质解释人员所掌握测量地区丰富的已知地质层位信息和储层油汽水信息，在满足所研究的目的层的对层和标定的条件下，对测井参数进行交互解释后，才能作为约束用的附加信息参加反演过程。对一层砂岩一层泥岩的砂泥岩薄互层而言，必须在测井参数曲线上进行详细的解释，识别出砂泥岩薄互层。如果没有可靠的井中薄互层解释，就谈不上薄互层反演；同样没有井参数的有效约束，要想从地面地震记录中直接反演薄互层的可能性是很小的。

地震数据纵向分辨率虽然很低，不能满足储层预测和油藏描述的要求，但是地震数据横向具有一定的连续性和分辨率。而测井资料纵向分辨率虽然很高，能满足储层预测和油藏描述的要求，但是测井资料横向不具有分辨能力。目前不能有效的利用地震数据横向连续性和测井资料纵向高分辨率的互补特性，并使井中数据的高分辨率特征拓展到井间地震波反演的数据中。

发明内容

本发明目的在于提供一种直接利用常规叠前(纵波)地震数据和测井数据，反演反映地下地层岩性、油气水界面变化的叠前地震数据地层弹性常数参数反演的方法。

本发明通过以下技术方案实现，包括以下步骤：

1)激发及记录地震波，采集地震数据，对叠前地震数据进行处理；

步骤1)所述的处理包括地表一致性振幅处理和地表一致性反褶积处理，速度分析、动校正和剩余静校正、剩余振幅补偿，叠加和叠前偏移处理，形成偏移归位的叠前道集数据。

2)在工区进行测井，取得测井数据和地质、岩心数据，确定入射角范围，进行角道集分析处理，提取角道集地震数据，最后进行角道集叠加，形成角道集叠加数据和全部入射角的叠加数据；

步骤2)所述的测井就是使用测井技术和设备，测量得到声波时差数据和密度数据，得到测井数据。

3)将地震数据和测井数据层位标定与对层，同时生成井中数据，采用常规的测井方法得到测井数据，得到声波时差曲线和密度曲线，并根据叠后地震数据和已采集的先验地质地层层位信息和钻井地层信息，把深度域测井的声波时差曲线和密度曲线标定为时间域，与叠后地震数据和解释的层位对应一致，同时生成井中波阻抗、井中纵波速度、井中横波速度、井中地层密度数据；

步骤3)所述的先验地质地层层位信息是实际地质探查得到的已知工区地层构造。

步骤3)所述的钻井地层信息是通过钻井的方法取出大小地层岩心样本反映的地质层位构造。

4)将步骤2)的不同入射角叠加的地震数据，反演得到不同入射角叠加地震数据的波阻抗参数；

步骤4)所述的反演包括以下步骤：

(1)在地震数据上，根据已知的先验地质地层层位信息和钻井地层信息，对地震数据体进行层位解释，拾取地质地层层位对应的地震地层层位，并对拾取的层位进行检验和校正处理以及平滑和内插处理，使得层位闭合和平滑；

(2)对测井的声波时差曲线和密度曲线进行分析，编辑和剔除其中异常值，并根据叠后地震数据和解释的层位，已知的先验地质地层层位信息和钻井地层信息，对测井的声波时差曲线和密度曲线进行标定，把深度域测井的声波时差曲线和密度曲线标定为时间域，与叠后地震数据和解释的层位一致起来，并生成井中波阻抗；

(3)对输入的叠后地震数据进行小波奇性分析与小波特征提取。选择小波尺度函数，构造滤波器，采用逐次分离的方法进行分解，得到了地震道的奇性特征；

(4)利用解释的地震层位信息和井中的波阻抗，生成初始波阻抗模型；

(5)把要反演的区间分成若干个子区间，并且使得相邻区间之间有一定的重叠；

(6)在第一个区间上，依据地震数据和初始波阻抗模型，利用一维波动方程，采用非线性最小二乘拟合方法，求解第一个区间上的波阻抗；

(7)利用一维波动方程和第一个区间上计算的波阻抗，将波场向下延拓到整个第一个区间，舍去区间的重叠部分，从第二个区间的起始部分开始，将第一个区间向下延拓到第二个区间的波场作为第二个区间的初始条件；依次类推，求得整个区间上的波阻抗；

(8)重复(6)-(7)直到所有区间求解完毕，得到一个地震道整个区间的波阻抗参数；

(9)对于所有的地震道，重复(3)-(8)过程，得到所有地震道的最终波阻抗参数。

5)对不同入射角叠加的地震数据的波阻抗参数，计算其对应的反射系数；

6)将井位置处测井得到的纵波速度、横波速度、地层密度和井位置处地震道的波阻抗参数、反射系数、地震数据、层位信息，构成速度参数、波阻抗参数、反射系数、地震数据、层位数据对；

7)根据数据对，计算函数映射网络模型的权函数和模型参数；

步骤7)所述的函数映射网络模型是：

$y=g\left(x\right)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{W}_k\left(x\right)y_k$

$r_k\left(x\right)={\left[{(x-x_k)}^T{A}^{T}A(x-x_k)\right]}^{\frac{1}{2}}---\left(1\right)$

其中：

m表示样本个数，

表示地震道波阻抗参数，

表示反射系数，

表示地震数据，

表示层位信息。y_k分别表示纵波速度参数、横波速度参数和地层密度参数；

向量x＝(x¹，x²，…，xⁿ)^T，这里n＝4表示向量x有四个参数即：

波阻抗参数，

反射系数，

地震数据，

层位信息；

x∈Rⁿ为输入，y∈R为输出，{(x_k,y_k)，k＝1，2，...，m}为经验样本集，权函数W_k和模型参数a，A＝{a_ij}。

步骤7)所述的模型参数计算算法为：

$E=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{(y_k-{\hat{y}}_k)}^2---\left(2\right)$

其中 ${\hat{y}}_k=g\left(x_k\right),$ k＝1，2，Λ，m；

$a_{\mathrm{ij}}^t=a_{\mathrm{ij}}^{t-1}-\mathrm{\β}\frac{\∂E}{\∂a_{\mathrm{ij}}}---\left(3\right)$

其中β＞0为步长，t为迭代步数，且

$\frac{\∂E}{\∂a_{\mathrm{ij}}}=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}({\hat{y}}_k-y_k)\frac{\∂{\hat{y}}_k}{\∂a_{\mathrm{ij}}},$ $\frac{\∂\hat{y}}{\∂a_{\mathrm{ij}}}=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{y}_k\frac{\∂W_k}{\∂a_{\mathrm{ij}}},$

$\frac{\∂r_l}{\∂a_{\mathrm{ij}}}=\left[(x^j-x_l^j)\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{a}_{\mathrm{ik}}(x^k-x_l^k)\right]/r_l---\left(4\right)$

式中：l＝1，2，Λ，m，i，j＝1，2，Λ，n。

8)计算纵波速度参数、横波速度参数和地层密度参数。

步骤8)所述的计算是根据步骤7)过程确定的地震道波阻抗参数、反射系数、地震数据、层位等与速度参数的函数映射网络模型，使用函数映射网络模型计算算法，得到所有地震道的地层纵波速度参数V_p、地层横波速度参数V_s和地层密度参数ρ。

9)计算弹性波阻抗参数；

步骤9)弹性波阻抗计算采用以下公式：

$\mathrm{EEI}=A{I}_0{\left(\frac{\mathrm{AI}}{{\mathrm{AI}}_0}\right)}^{\cos \mathrm{\χ}}{\left(\frac{\mathrm{GI}}{{\mathrm{GI}}_0}\right)}^{\sin \mathrm{\χ}}---\left(6\right)$

其中，

tanχ＝sin ²θ， ${\mathrm{GI}=V}_p{V}_s^{-8K}{\mathrm{\ρ}}^{-4K}$

${G{I}_0=V}_{p0}{V}_{s0}^{-8K}{\mathrm{\ρ}}_0^{-4K},$ $K=\frac{1}{2}(\frac{V_s^2}{V_p^2}+\frac{V_{s0}^2}{V_{p0}^2})---\left(7\right)$

其中，θ表示地层入射角，AI表示地层波阻抗参数，AI₀表示地层参考波阻抗参数，V_p表示地层纵波速度参数，V_p0表示地层参考纵波速度参数，V_s表示地层横波速度参数，V_s0表示地层参考横波速度参数，ρ表示地层密度参数，ρ₀表示地层参考密度参数。

10)使用反射系数与入射角对应方程拟合系数参数；

步骤10)所述的反射系数与入射角的对应计算是：

R_i＝A+Bsin²θ_i+Csin²θ_itg²θ_i           (8)

其中：R_i是地震反射系数，θ_i是入射角。已知一组(R_i，θ_i)，i＝1，2，K，N，拟合出系数A、B和C。

步骤10)所述的拟合系数A、B和C采用下列矩阵方程计算：

$\left(\begin{array}{ccc}{u}^{T}u& u^{T}v& u^{T}w\\ v^{T}u& v^{T}v& v^{T}w\\ w^{T}u& w^{T}v& w^{T}w\end{array}\right)\left(\begin{array}{c}A\\ B\\ C\end{array}\right)=\left(\begin{array}{c}{r}^{T}u\\ r^{T}v\\ r^{T}w\end{array}\right)---\left(9\right)$

对于不同入射角叠加地震数据的反射系数，求出对应的系数A、B和C；

其中：

u^T＝(1，1，K，1)，v^T＝(sin²θ₁，sin²θ₂，K，sin²θ_N)

w^T＝(sin²θ₁tg²θ₁，sin²θ₂tg²θ₂，K，sin²θ_Ntg²θ_N)，r^T＝(R₁，R₂，K，R_N)是已知数据向量，“T”表示向量转置。

11)计算弹性常数；

步骤11)所述的弹性常数包括：体积模量κ、拉梅常数λ、剪切模量μ、杨氏模量E和泊松比σ。

步骤11)所述的计算公式是：

对体积模量κ、拉梅常数λ和剪切模量μ：

$\mathrm{\Δ}{\mathrm{\κ}}_i=\frac{3A+B+2C}{1.5}{Z}_i{V}_i,$ Δλ_i＝(2A+B+C)Z_iV_i， $\mathrm{\Δ}{\mathrm{\μ}}_i=\frac{C-B}{2}{Z}_i{V}_i---\left(10\right)$

其中：Z_i是地震波阻抗参数，V_i是地震纵波速度参数；杨氏模量E和泊松比σ：

$E=\frac{\mathrm{\λ}(3\mathrm{\λ}+2\mathrm{\μ})}{\mathrm{\λ}+\mathrm{\μ}},$ $\mathrm{\σ}=\frac{\mathrm{\λ}}{2(\mathrm{\λ}+\mathrm{\μ})}---\left(11\right)$

12)绘制弹性常数剖面，用于储层岩性识别、油气预测、油水界面确定和油气藏的描述。

本发明有效地利用了常规地震数据和测井数据反演计算弹性常数，可利用常规地震数据和测井数据反演计算出弹性常数，对断层、尖灭带的反演有一定的适应能力，本发明具有适应范围大、分辨率高、计算量小、计算速度快、稳定性好、计算精度高、具有一定的抗噪能力的特点。

附图说明

图1不同入射角叠加数据叠加剖面对比图

(a)入射角0度(全部入射角)

(b)入射角5度

(c)入射角12度

(d)入射角18度

(e)入射角24度

(f)入射角30度

图2不同入射角叠加数据反演的波阻抗参数对比图

(a)入射角0度(全部入射角)

(b)入射角5度

(c)入射角12度

(d)入射角18度

(e)入射角24度

(f)入射角30度

图3反演的弹性波阻抗参数对比图

(a)入射角5度

(b)入射角12度

(c)入射角18度

(d)入射角24度

(e)入射角30度

图4拟合系数参数对比图

(a)拟合系数A剖面

(b)拟合系数B剖面

(c)拟合系数C剖面

图5弹性常数对比图

(a)体积模量κ

(b)剪切模量μ

(c)拉梅常数λ

(d)杨氏模量E

(e)泊松比σ

具体实施方式

本发明首先利用叠前地震数据分析方法，对叠前地震数据按照入射角大小不同，进行不同入射角数据叠加，形成多个不同入射角的叠加地震数据；对这些不同入射角的叠加地震数据进行分析和对比，得出能够很好反映地下储层变化的不同入射角叠加地震数据。然后利用地震数据、测井资料的互补特性，研究地震、测井资料联合波动方程反演技术，以实现井间地层波阻抗参数的空间分布。再利用函数映射网络模型标定技术，在地震数据和测井资料的约束下，把井间地层波阻抗参数标定成为速度参数。最后利用不同炮检距叠加地震数据反演速度参数，计算弹性常数，在油田的勘探、开发、开采过程中为储层预测、油水界面识别和油藏描述提供高分辨率的弹性常数的技术的方法。本发明反演出的地层物性参数包括：地层弹性波阻抗和体积模量κ、拉梅常数λ、剪切模量μ、杨氏模量E和泊松比σ五个弹性常数参数。

本发明的一种叠前地震数据地层弹性常数参数反演的方法，包括以下步骤：

1)采用常规的地震勘探方法采集地震数据，对叠前地震数据进行处理；

步骤1)所述的处理包括地表一致性振幅处理和地表一致性反褶积处理，速度分析、动校正和剩余静校正、剩余振幅补偿，叠加和叠前偏移处理，形成偏移归位的叠前道集数据。

2)对叠前地震数据进行入射角分析，得出能够反映地下储层变化的不同方位角叠加地震数据；

步骤2)所述的入射角分析包括建立层速度场(速度分析)、角道分析、角道叠加。对于经过叠前去噪、叠前提高分辨率处理、叠前偏移归位处理的叠前地震数据，进行剩余振幅补偿处理，并且计算层速度，然后结合测井数据和地质、岩心数据分析，确定入射角范围，进行角道集分析处理，提取角道集地震数据，最后进行角道集叠加，形成角道集叠加数据。

3)地震数据和测井数据层位标定与对层，同时生成井中数据；

步骤2)所述的地震数据和测井数据层位标定与对层，同时生成井中数据，就是采用常规的测井方法得到测井数据，得到声波时差曲线和密度曲线，并根据叠后地震数据和已采集的先验地质地层层位信息和钻井地层信息，把深度域测井的声波时差曲线和密度曲线标定为时间域，与叠后地震数据和解释的层位对应一致，同时生成井中波阻抗和井中速度数据(纵波速度、横波速度、地层密度)；

4)将步骤2)的不同入射角叠加的地震数据，反演得到不同入射角叠加地震数据的波阻抗参数；

步骤4)所述的反演包括以下步骤：(1)在地震数据上，根据已知的先验地质地层层位信息和钻井地层信息，对地震数据体进行层位解释，拾取地质地层层位对应的地震地层层位，并对拾取的层位进行检验和校正处理以及平滑和内插处理，使得层位闭合和平滑；(2)对测井的声波时差曲线和密度曲线进行分析，编辑和剔除其中异常值，并根据叠后地震数据和解释的层位，已知的先验地质地层层位信息和钻井地层信息，对测井的声波时差曲线和密度曲线进行标定，把深度域测井的声波时差曲线和密度曲线标定为时间域，与叠后地震数据和解释的层位一致起来，并生成井中波阻抗；(3)对输入的叠后地震数据进行小波奇性分析与小波特征提取。选择小波尺度函数，构造滤波器，采用逐次分离的方法进行分解，得到了地震道的奇性特征。(4)利用解释的地震层位信息和井中的波阻抗，生成初始波阻抗模型；(5)把要反演的区间分成若干个子区间，并且使得相邻区间之间有一定的重叠；(6)在第一个区间上，依据地震数据和初始波阻抗模型，利用一维波动方程，采用非线性最小二乘拟合方法，求解第一个区间上的波阻抗；(7)利用一维波动方程和第一个区间上计算的波阻抗，将波场向下延拓到整个第一个区间，舍去区间的重叠部分，从第二个区间的起始部分开始，将第一个区间向下延拓到第二个区间的波场作为第二个区间的初始条件；依次类推，求得整个区间上的波阻抗.(8)重复(6)-(7)直到所有区间求解完毕，得到一个地震道整个区间的波阻抗参数。(9)对于所有的地震道，重复(3)-(8)过程，得到所有地震道的最终波阻抗参数。

5)对不同入射角叠加地震数据的波阻抗参数，计算其对应的反射系数；

如果已知地震波阻抗参数，那么可以计算反射系数参数：

$R_i=\frac{Z_{i+1}-z_i}{Z_{i+1}+Z_i}---\left(1\right)$

这里，R_i是地震反射系数，Z_i是地震波阻抗参数。对于不同入射角叠加地震数据的波阻抗参数，使用方程(1)可以计算其对应的反射系数。

6)将井位置处测井得到的速度参数(纵波速度、横波速度、地层密度)和井位置处地震道的波阻抗参数、反射系数、地震数据、层位信息，构成速度参数(纵波速度、横波速度、地层密度)、波阻抗参数、反射系数、地震数据、层位数据对；

地层速度参数的变化，必然引起地震特征参数的变化，包括波阻抗参数的变化；也就是说地震特征与储层速度参数之间，必定存在某种映射关系。只要准确的建立起这种映射关系，就可以由地震特征来预测储层速度参数的空间分布。而函数映射网络模型具有非常强的非线性映射功能，可以自动归纳、总结出隐蔽的规律，这为储层速度参数横向预测提供了可能。函数映射网络模型过程包括网络模型参数计算和利用网络模型参数进行预测两个过程。我们通过函数映射网络模型，来建立地层波阻抗参数等地震特征与储层速度参数之间的这种映射关系。

函数映射网络模型包括经验样本集{(x_k，y_k)，k＝1，2，...，m}，四组功能函数：距离函数r_k(x)，活化函数，权函数w_k(x)，输出函数g(x)，以及模型参数a，A＝{a_ij}。

建立速度参数(纵波速度、横波速度、地层密度)、波阻抗参数、地震数据、层位数据对。把井位置处测井得到的速度参数(纵波速度、横波速度、地层密度)和井位置处地震道的波阻抗参数、反射系数、地震数据、层位等，构成速度参数、波阻抗参数、反射系数、地震数据、层位数据对，用于建立井位置处地层速度参数(纵波速度、横波速度、地层密度)与地震道的波阻抗参数、反射系数、地震数据、层位等地震特征参数之间的映射关系。

在函数映射网络模型的结构中，这一步就是经验样本数据集：{(x_k，y_k)，k＝1，2，...，m}，其中m表示井中速度参数的个数，y_k表示井中第k速度参数值(纵波速度、横波速度、地层密度)，x_k＝(地震道波阻抗参数，反射系数，地震数据，层位参数)，这样就构成了速度参数(纵波速度、横波速度、地层密度)、波阻抗参数、反射系数、地震数据、层位数据对，用于建立井位置处地层速度参数与地震道的波阻抗参数、反射系数、地震数据、层位等地震特征参数之间的映射关系。

7)根据速度参数(纵波速度、横波速度、地层密度)、波阻抗参数、反射系数、地震数据、层位数据对，使用函数映射网络模型计算算法，得到函数映射网络模型的权函数和模型参数；

本发明所述的一种叠前地震数据地层弹性常数参数反演方法，函数映射网络模型计算算法采用基于梯度的最速下降法。建立目标函数

$E=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{(y_k-{\hat{y}}_k)}^2---\left(2\right)$

其中 ${\hat{y}}_k=g\left(x_k\right),$ k＝1，2，Λ，m。那么

$a_{\mathrm{ij}}^t=a_{\mathrm{ij}}^{t-i}-\mathrm{\β}\frac{\∂E}{\∂a_{\mathrm{ij}}}---\left(3\right)$

其中β＞0为步长，t为迭代步数。且

$\frac{\∂E}{\∂a_{\mathrm{ij}}}=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}({\hat{y}}_k-y_k)\frac{\∂{\hat{y}}_k}{\∂a_{\mathrm{ij}}},$ $\frac{\∂\hat{y}}{\∂a_{\mathrm{ij}}}=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{y}_k\frac{\∂W_k}{\∂a_{\mathrm{ij}}},$

$\frac{\∂r_l}{\∂a_{\mathrm{ij}}}=\left[(x^j-x_l^j)\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{a}_{\mathrm{ik}}(x^k-x_l^k)\right]/r_l---\left(4\right)$

这里l＝1，2，Λ，m，i，j＝1，2，Λ，n。根据函数映射网络模型计算算法的这些公式，对函数映射网络模型进行计算，以确定地震道的波阻抗参数、反射系数、地震数据、层位等与速度参数(纵波速度、横波速度、地层密度)的映射关系，计算算法就是根据经验样本集和函数映射网络模型结构，确定函数映射网络模型的权函数W_k和模型参数a_ij。

8)对于所有的地震道，根据步骤7)过程确定的地震道波阻抗参数、反射系数、地震数据、层位等与速度参数(纵波速度、横波速度、地层密度)的函数映射网络模型，计算得到所有地震道的速度参数(纵波速度、横波速度、地层密度)；

本发明所述的一种叠前地震数据地层弹性常数参数反演方法，函数映射网络模型计算方法是：

令x＝(x¹，x²，…，xⁿ)^T，

函数映射网络模型计算如下：

$y=g\left(x\right)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{W}_k\left(x\right)y_k$

$r_k\left(x\right)={\left[{(x-x_k)}^T{A}^{T}A(x-x_k)\right]}^{\frac{1}{2}}---\left(5\right)$

其中x∈Rⁿ为输入，y∈R为输出，{(x_k，y_k)，k＝1，2，...，m}为经验样本集，模型参数由模型参数计算算法确定。

根据函数映射网络模型计算算法，对所有的地震数据道，由地震道的波阻抗参数、反射系数、地震数据、层位参数，可以计算出对应位置处的纵波速度参数、横波速度参数、密度参数。这样可以得到所有地震道的纵波速度参数、横波速度参数、密度参数。

9)计算弹性波阻抗参数；

根据扩展弹性波阻抗的公式：

$\mathrm{EEI}={\mathrm{AI}}_0{\left(\frac{\mathrm{AI}}{{\mathrm{AI}}_0}\right)}^{\cos \mathrm{\χ}}{\left(\frac{\mathrm{GI}}{{\mathrm{GI}}_0}\right)}^{\sin \mathrm{\χ}}---\left(6\right)$

其中

tanχ＝sin ²θ， $\mathrm{GI}=V_p{V}_s^{-8K}{\mathrm{\ρ}}^{-4K}$

${\mathrm{GI}}_0=V_{p0}{V}_{s0}^{-8K}{\mathrm{\ρ}}_0^{-4K},$ $K=\frac{1}{2}(\frac{V_s^2}{V_p^2}+\frac{V_{s0}^2}{V_{p0}^2})---\left(7\right)$

这里，θ表示地层入射角，AI表示地层波阻抗参数，AI₀表示地层参考波阻抗参数，V_p表示地层纵波速度参数，V_p0表示地层参考纵波速度参数，V_s表示地层横波速度参数，V_s0表示地层参考横波速度参数，ρ表示地层密度参数，ρ₀表示地层参考密度参数。

根据步骤4)反演确定地层波阻抗参数AI，根据步骤6)---步骤8)反演计算确定地层纵波参数V_p、地层横波速度参数V_s和地层密度参数ρ，给定地层参考纵波参数V_p0、地层参考横波速度参数V_s0和地层密度参数ρ₀，由方程(6)和方程(7)可以计算不同入射角的弹性波阻抗参数。

10)使用反射系数与入射角对应方程，拟合系数参数；

反射系数与入射角的对应方程为：

R_i＝A+Bsin²θ_i+Csin²θ_itg²θ_i           (8)

这里，R_i是地震反射系数，θ_i是入射角。已知一组(R_i，θ_i)，i＝1，2，K，N，使用方程(8)可以拟合出系数A、B和C。为了拟合系数A、B和C，首先建立目标函数：

$Q=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(R_i-A-B{\sin }^2{\mathrm{\θ}}_i-C{\sin }^2{\mathrm{\θ}}_{i}t{g}^2{\mathrm{\θ}}_i)}^2---\left(9\right)$

令

r^T＝(R₁，R₂，K，R_N)

u^T＝(1，1，K，1)

v^T＝(sin²θ₁，sin²θ₂，K，sin²θ_N)

w^T＝(sin²θ₁tg²θ₁，sin²θ₂tg²θ₂，K，sin²θ_Ntg²θ_N)

其中“T”表示向量转置，则方程(9)可以表示为

Q＝(r-Au-Bv-Cw)^T(r-Au-Bv-Cw)

 ＝r^Tr+A²u^Tu+B²v^Tv+C²w^Tw         (10)

 -2Ar^Tu-2Br^Tv-2Cr^Tw+2ABu^Tv+2ACu^Tw+2BCv^Tw

目标函数Q对系数A、B和C分别求导数，且等于零，有

$\frac{\∂Q}{\∂A}=2A{u}^{T}u-2r^{T}u+2B{u}^{T}v+2C{u}^{T}w=0$

$\frac{\∂Q}{\∂B}=2B{v}^{T}v-2r^{T}v+2A{u}^{T}v+2C{v}^{T}w=0$

$\frac{\∂Q}{\∂C}=2C{w}^{T}w-2r^{T}w+2A{u}^{T}w+2B{v}^{T}w=0---\left(11\right)$

整理，有矩阵方程

$\left(\begin{array}{ccc}{u}^{T}u& u^{T}v& u^{T}w\\ v^{T}u& v^{T}v& v^{T}w\\ w^{T}u& w^{T}v& w^{T}w\end{array}\right)\left(\begin{array}{c}A\\ B\\ C\end{array}\right)=\left(\begin{array}{c}{r}^{T}u\\ r^{T}v\\ r^{T}w\end{array}\right)---\left(12\right)$

这是一个主对角线占优的三元一次线性方程组，因此方程有唯一解。由方程(12)可以求出系数A、B和C。对于不同入射角叠加地震数据的反射系数，可以求出对应的系数A、B和C。

11)利用系数A、B和C，计算弹性常数；

系数A、B和C与体积模量κ、拉梅常数λ和剪切模量μ变化之间有如下的关系，即

$\mathrm{\Δ}{\mathrm{\κ}}_i=\frac{3A+B+2C}{1.5}{Z}_i{V}_i$

Δλ_i＝(2A+B+C)Z_iV_i

$\mathrm{\Δ}{\mathrm{\μ}}_i=\frac{C-B}{2}{Z}_i{V}_i---\left(13\right)$

这里Z_i是地震波阻抗参数，由步骤4)计算取出，V_i是地震纵波速度参数，由步骤8)计算取出。对方程(13)分别进行求和，可以得出体积模量κ、拉梅常数λ和剪切模量μ。另外根据剪切模量和拉梅常数，还可以导出杨氏模量E和泊松比σ剖面：

$E=\frac{\mathrm{\λ}(3\mathrm{\λ}+2\mathrm{\μ})}{\mathrm{\λ}+\mathrm{\μ}}$

$\mathrm{\σ}=\frac{\mathrm{\λ}}{2(\mathrm{\λ}+\mathrm{\μ})}---\left(14\right)$

这样就可以计算出五个弹性常数参数：体积模量κ、拉梅常数λ、剪切模量μ、杨氏模量E和泊松比σ。

12)绘制弹性波阻抗和弹性常数剖面，用于储层岩性识别、油气预测、油水界面确定和油气藏的描述。

本发明实施情况如下：

首先对叠前地震数据进行不同入射角道集分析，抽取不同入射角的数据，抽取入射角范围分别是全部角度、0-10度、8-16度、14-22度、20-28度、大于22度(30度)，构成全部入射角(0度)、5度、12度、18度、24度和30度角道集，将其叠加产生0度、5度、12度、18度、24度和30度角道集叠加数据。图1是不同入射角叠加数据叠加剖面对比图，(a)入射角0度，(b)入射角5度，(c)入射角12度，(d)入射角18度，(e)入射角24度，(f)入射角30度。对这些不同入射角叠加地震数据，使用地震数据波阻抗反演方法，得到不同入射角叠加地震数据的波阻抗参数。图2为不同入射角叠加数据反演的波阻抗参数对比图，(a)入射角0度，(b)入射角5度，(c)入射角12度，(d)入射角18度，(e)入射角24度，(f)入射角30度。然后对这些不同入射角叠加数据反演的波阻抗参数，分别反演得到纵波速度参数、横波速度参数和密度参数。使用纵波速度参数、横波速度参数和密度参数，反演计算弹性波阻抗参数。图3是反演的弹性波阻抗参数对比图，(a)入射角5度，(b)入射角12度，(c)入射角18度，(d)入射角22度，(e)入射角30度。对不同入射角叠加数据反演的波阻抗参数，计算其对应的反射系数序列，对不同入射角反射系数序列，拟合计算拟合系数A、B、C剖面。图4是拟合系数参数剖面对比，(a)拟合系数A剖面，(b)拟合系数B剖面，(c)拟合系数C剖面。使用拟合系数A、B、C参数和纵波速度参数、密度参数，计算地层弹性常数参数，图5是弹性常数对比剖面，(a)体积模量κ，(b)剪切模量μ，(c)拉梅常数λ，(d)杨氏模量E，(e)泊松比σ。

Claims (6)

1.一种叠前地震数据地层弹性常数参数反演的方法，特征是通过以下技术方案实现，包括以下步骤： 

1）激发及记录地震波，采集地震数据，对叠前地震数据进行处理； 

2）在工区进行测井，取得测井数据和地质、岩心数据，确定入射角范围，进行角道集分析处理，提取角道集地震数据，最后进行角道集叠加，形成角道集叠加数据和全部入射角的叠加数据； 

3）将地震数据和测井数据层位标定与对层，同时生成井中数据，采用常规的测井方法得到测井数据，得到声波时差曲线和密度曲线，并根据叠后地震数据和已采集的先验地质地层层位信息和钻井地层信息，把深度域测井的声波时差曲线和密度曲线标定为时间域，与叠后地震数据和解释的层位对应一致，同时生成井中波阻抗、井中纵波速度、井中横波速度、井中地层密度数据； 

4）将步骤2）的不同入射角叠加的地震数据，反演得到不同入射角叠加地震数据的波阻抗参数； 

5）对不同入射角叠加的地震数据的波阻抗参数，计算其对应的反射系数； 

6）将井位置处测井得到的纵波速度、横波速度、地层密度和井位置处地震道的波阻抗参数、反射系数、地震数据、层位信息，构成速度参数、波阻抗参数、反射系数、地震数据、层位数据对； 

7）根据数据对，计算函数映射网络模型的权函数和模型参数； 

所述的函数映射网络模型是： 

其中： 

m表示样本个数，y_k分别表示纵波速度参数、横波速度参数和地层密度参数； 

向量x=(x¹,x²,…,xⁿ)^T，

这里n=4表示向量x有四个参数即：

波阻抗参数，反射系数，地震数据，

层位信息； 

x∈Rⁿ为输入,y∈R为输出,{(x_k,y_k),k=1,2,...,m}为经验样本集，权函数W_k和模型参数a,A={a_ij}； 

所述的模型参数计算算法为： 

其中

其中β>0为步长，t为迭代步数，且 

式中：l=1,2,…,m，i,j=1,2,…,n； 

8）计算纵波速度参数、横波速度参数和地层密度参数； 

9）计算弹性波阻抗参数； 

所述的弹性波阻抗参数计算采用以下公式： 

其中， 

其中，θ表示地层入射角，AI表示地层波阻抗参数，AI₀表示地层参考波阻抗参数，V_p表示地层纵波速度参数，V_p0表示地层参考纵波速度参数，V_s表示地层横波速度参数，V_s0表示地层参考横波速度参数，ρ表示地层密度参数，ρ₀表示地层参考密度参数； 

10）使用反射系数与入射角对应方程拟合系数参数； 

步骤10）所述的反射系数与入射角的对应计算是： 

R_i=A+Bsin²θ_i+Csin²θ_itg²θ_i     （8） 

其中：R_i是地震反射系数，θ_i是入射角；已知一组(R_i,θ_i)，i=1,2,...,N，拟合出系数A、B和C； 

所述的拟合系数A、B和C采用下列矩阵方程计算： 

对于不同入射角叠加地震数据的反射系数，求出对应的系数A、B和C； 

其中： 

u^T=(1,1,...,1)，v^T=(sin²θ₁,sin²θ₂,...,sin²θ_N) 

w^T=(sin²θ₁tg²θ₁,sin²θ₂tg²θ₂,...,sin²θ_Ntg²θ_N)，r^T=(R₁,R₂,...,R_N) 

是已知数据向量，“T”表示向量转置； 

11）计算弹性常数； 

所述的计算弹性常数公式是： 

对体积模量κ、拉梅常数λ和剪切模量μ： 

其中：Z_i是地震波阻抗参数，V_i是地震纵波速度参数； 

杨氏模量E和泊松比σ： 

12）绘制弹性常数剖面，用于储层岩性识别、油气预测、油水界面确定和油气藏的描述。 

2.根据权利要求1的方法，特征是步骤1）所述的处理包括地表一致性振幅处理和地表一致性反褶积处理，速度分析、动校正和剩余静校正、剩余振幅补偿，叠加和叠前偏移处理，形成偏移归位的叠前道集数据。 

3.根据权利要求1的方法，特征是步骤3）所述的先验地质地层层位信息是实际地质探查得到的已知工区地层构造。 

4.根据权利要求1的方法，特征是步骤3）所述的钻井地层信息是通过钻井的方法取出大小地层岩心样本反映的地质层位构造。 

5.根据权利要求1的方法，特征是步骤4）所述的反演包括以下步骤： 

（1）在地震数据上，根据已知的先验地质地层层位信息和钻井地层信息，对地震数据体进行层位解释，拾取地质地层层位对应的地震地层层位，并对拾取的层位进行检验和校正处理以及平滑和内插处理，使得层位闭合和平滑； 

（2）对测井的声波时差曲线和密度曲线进行分析，编辑和剔除其中异常值，并根据叠后地震数据和解释的层位，已知的先验地质地层层位信息和钻井地层信息，对测井的声波时差曲线和密度曲线进行标定，把深度域测井的声波时差曲线和密度曲线标定为时间域，与叠后地震数据和解释的层位一致起来，并生成井中波阻抗； 

（3）对输入的叠后地震数据进行小波奇性分析与小波特征提取，选择小波尺度函数，构造滤波器，采用逐次分离的方法进行分解，得到了地震道的奇性特征； 

（4）利用解释的地震层位信息和井中的波阻抗，生成初始波阻抗模型； 

（5）把要反演的区间分成若干个子区间，并且使得相邻区间之间有一定的重叠； 

（6）在第一个区间上，依据地震数据和初始波阻抗模型，利用一维波动方程，采用非线性最小二乘拟合方法，求解第一个区间上的波阻抗； 

（7）利用一维波动方程和第一个区间上计算的波阻抗，将波场向下延拓到整个第一个区间，舍去区间的重叠部分，从第二个区间的起始部分开始，将第一 个区间向下延拓到第二个区间的波场作为第二个区间的初始条件；依次类推,求得整个区间上的波阻抗； 

（8）重复（6）-（7）直到所有区间求解完毕，得到一个地震道整个区间的波阻抗参数； 

（9）对于所有的地震道，重复（3）-（8）过程，得到所有地震道的最终波阻抗参数。 

6.根据权利要求1的方法，特征是步骤8）所述的计算是根据步骤7）过程确定的地震道波阻抗参数、反射系数、地震数据、层位等与速度参数的函数映射网络模型，使用函数映射网络模型计算算法，得到所有地震道的地层纵波速度参数V_p、地层横波速度参数V_s和地层密度参数ρ。 

Images