本发明属于勘探地球物理领域，特别涉及一种地层几何曲率估计方法。
背景技术：
石油及天然气是一种重要的战略资源，在我国国民经济及国家安全中占有举足轻重的地位。地震属性是对地震勘探数据的一种度量，能够在视觉上增强具有解释意义的特征，对地层特征以及储层特性有直接反映。因此，地震属性可在油气藏勘探开发中发挥重要的作用。几何属性是地震资料解释中最为重要的地震属性之一，主要用于增强和显示地震层位的几何形态，主要包括倾角\方位角、连续性和地层曲率等。倾角\方位角属性体包含有重要的地震地层学和地理学信息，不但可直接用于勘探工区的解释，还可为后续处理及解释提供基础数据(例如结构曲率计算、结构导向滤波。而地层结构曲率可通过对视倾角求取一阶导数而获取,可用于描述地质体的几何形态变化，不但对地层的弯曲、褶皱及断层等结构反映敏感，而且还可用于监控地震资料处理的质量。因此，准确稳定地从三维叠后数据中估计地层曲率具有重要的理论价值及应用价值。
目前已有多项技术可用于地层曲率估计，主要技术如下所示：
现有技术1：基于层位拾取的地层曲率估计方法
该类方法首先利用软件人工或者自动拾取目的层位，然后计算拾取到层位的二阶空间导数。基于拾取层位的二阶空间导数可以估计地层的曲率。
现有技术1的特点：
优点：简单易行。
缺点：1、拾取层位需要软件的协助。
2、只能针对目的层位计算地层曲率，不能针对这个三维地震数据体计算。
3、受人为干扰因素影响较大。
现有技术2：基于地震道互相关的地层曲率估计方法
该类方法首先计算三维地震数据各地震道与相邻道之间的互相关，进而基于互相关值估计相邻道之间的时差以获得视倾角，最后对视倾角计算空间偏导数以估计地层曲率。
现有技术2的特点：
优点：实现简单、计算量小且不受人为因素干扰。
缺点：1、由于仅采用相邻道计算互相关，受噪声影响较大。
2、振幅横向变化及跨越断层时误差较大。
现有技术3：基于倾角扫描的地层曲率估计方法
该类方法首先按照一定的采样间隔以及范围设定一系列视倾角，然后沿每个视倾角计算相似度或者相干值，选取具有最大的相似度(或相干值)的视倾角做为分析点的倾角，最后对倾角计算空间偏导数以估计地层曲率。
现有技术3的特点：
优点：较为稳定、抗噪性能好且不受人为因素干扰。
缺点：1、由于每个设定的视倾角下都要计算相似度或，导致计算量较大。
2、倾角估计的精度受限于视倾角的采样间隔。
3、分析窗跨越断层时误差较大。
现有技术4：基于梯度结构张量的地层曲率估计方法
该类方法首先在分析窗内计算所有地震数据的梯度向量，然后由梯度向量构建梯度协方差矩阵，接着对协方差矩阵进行特征分解以得到主特征向量，利用主特征向量估计地层的倾角，最后对倾角计算空间偏导数以估计地层倾角。
现有技术4的特点：
优点：估计准确、抗噪性能好且不受人为因素干扰。
缺点：振幅横向变化及跨越断层时曲率估计误差较大。
技术实现要素：
本发明的目的在于提供一种联合梯度结构张量和多窗分析的地层结构曲率估计方法，以解决上述技术问题。
为了实现上述目的，本发明采用如下技术方案：
联合梯度结构张量和多窗分析的地层结构曲率估计方法，包括以下步骤：
步骤01：采集三维叠后地震数据并利用Hilbert变换获得每道地震信号的复地震道；
步骤02：在复地震道的基础之上计算得到三维瞬时相位数据体，同时计算瞬时相位的梯度向量；
步骤03：针对三维数据体的每个样点附近选取众多相邻分析窗；
步骤04：在每个相邻分析窗内构建相位梯度的协方差矩阵并对其进行特征分解以得到三个特征值及特征向量，利用主特征向量计算视倾角，同时利用三个特征值计算相似度；
步骤05：在所有相邻分析窗内选择相似度最大的分析窗，将其对应的视倾角指定为分析样点的视倾角；
步骤06：判断三维数据中的所有样点是否都计算得到了视倾角：若没有返回步骤03；若完成进入步骤07；
步骤07：计算三维视倾角数据体的偏导数；
步骤08：利用视倾角的偏导数计算地层的结构曲率。
进一步的，步骤01具体包括：
若用f(x,y,t)表示三维地震数据，其中x和y表示三维地震数据在两个空间方向的指标，t表示三维地震数据在时间方向的指标；固定x和y，分别对每道信号f(x,y,t)进行Hilbert变换得到h(x,y,t)为：
其中τ为积分临时变量；则对应的复地震道c(x,y,t)为
c(x,y,t)＝f(x,y,t)+jh(x,y,t)＝A(x,y,t)ejφ(x,y,t)，
其中j为虚数单位，A(x,y,t)为瞬时振幅，φ(x,y,t)为瞬时相位。
进一步的，步骤02具体包括：
利用复地震道c(x,y,t)的实部f(x,y,t)和虚部h(x,y,t)计算得到三维瞬时相位体φ(x,y,t)如下：
同时可由三维瞬时相位体φ(x,y,t)计算得到瞬时相位的梯度向量如下：
进一步的，步骤03具体包括：
假设当前的分析点为(x,y,t)，在其附近选取众多相邻分析窗；假设分析窗大小为Mx×My×Mt；其中Mx，My及Mt均为奇数，且Mx＝2Wx+1，My＝2Wy+1，Mt＝2Wt+1；得到MxMyMt个相邻分析窗，且这些相邻分析窗均必须包含当前分析点(x,y,t)。
进一步的，步骤04具体包括：
在步骤03得到的每个相邻分析窗内利用步骤02得到的相位梯度构建协方差矩阵，构建当前分析点为(x,y,t)的相位梯度协方差矩阵C(x,y,t)：
其中mx，my及mt均为求和临时指标；
对相位梯度协方差矩阵C(x,y,t)进行特征值分解得到三个特征值μ1，μ2，μ3及对应的三个特征向量v1，v2，v3，并有μ1≥μ2≥μ3：
其中μ1及v1为协方差矩阵C(x,y,t)的主特征值及主特征向量；利用主特征向量v1的三个元素v1x，v1y及v1t估计在此相邻分析窗内地层沿x方向的视倾角px(x,y,t)及沿y方向的视倾角py(x,y,t)：
利用三个特征值μ1，μ2及μ3计算相似度similarity(x,y,t)如下：
在每个相邻分析窗内均计算视倾角和相似度；若分析点(x,y,t)共有M个相邻分析窗，记第m个相邻分析窗内估计的沿x及y方向的视倾角分别为px(m,x,y,t)及py(m,x,y,t)，记第m个相邻分析窗内估计的相似度为similarity(m,x,y,t)。
进一步的，步骤05具体包括：
分析点(x,y,t)周围共有M个相邻分析窗，选出相似度最大的相邻分析窗序号为win_opt(x,y,t)：
则分析点(x,y,t)处的沿x方向及y方向的视倾角分别为Px(x,y,t)及Py(x,y,t)：
Px(x,y,t)＝px(win_opt(x,y,t),x,y,t)，
Py(x,y,t)＝py(win_opt(x,y,t),x,y,t)。
进一步的，步骤07具体包括：
针对沿x方向的视倾角数据体Px(x,y,t)分别计算其沿x方向偏导数Pxx(x,y,t)及y方向的偏导数Pxy(x,y,t)：
针对沿y方向的视倾角数据体Py(x,y,t)分别计算其沿x方向偏导数Pyx(x,y,t)及y方向的偏导数Pyy(x,y,t)：
Δx和Δy为三维地震数据体在x和y方向上的采样间隔。
进一步的，步骤08具体包括：
利用步骤07中计算得到的三维视倾角数据体的偏导数Pxx(x,y,t)、Pxy(x,y,t)、Pyx(x,y,t)及Pyy(x,y,t)计算各点的最大正曲率MP(x,y,t)和最大负曲率MN(x,y,t)：
其中
a(x,y,t)＝0.5×Pxx(x,y,t)，
b(x,y,t)＝0.5×Pyy(x,y,t)，
c(x,y,t)＝0.5×[Pxy(x,y,t)+Pyx(x,y,t)]。
相对于现有技术，本发明具有以下有益效果：
本发明可更加稳定地估计三维地震数据的结构曲率属性，且减少断层等不连续结构对曲率属性估计的干扰。
附图表说明
图1为本发明流程图；
图2为分析点附近的相邻分析窗在空间方向及空间-时间方向示意图；其中，图2(a)为沿二维空间方向观察相邻分析窗示意图；图2(b)为沿空间-时间方向观察相邻分析窗示意图；
图3为分析点附件的所有相邻分析窗示意图；
图4为利用商业软件计算的最大正曲率沿层切片及利用本发明方法计算的的最大正曲率沿层切片；其中，图4(a)为利用商业软件计算的最大正曲率沿层切片；图4(b)为利用本发明方法计算的的最大正曲率沿层切片；
图5为利用商业软件计算的最大负曲率沿层切片及利用本发明方法计算的的最大负曲率沿层切片；其中，图5(a)为利用商业软件计算的最大负曲率沿层切片；图5(b)为利用本发明方法计算的的最大负曲率沿层切片。
具体实施方式
下面结合附图及表格对本发明做进一步详细的说明。
本发明为针对三维叠后地震数据利用梯度结构张量和多窗分析估计地层结构曲率的方法。本发明首先对地震数据进行复地震道分析以获得瞬时相位，然后在瞬时相位数据体上构建梯度结构张量并结合多窗分析技术进行地层倾角估计，进而对倾角计算空间导数，最终利用倾角的空间导数估计三维叠后地震数据的结构曲率，为后续地震资料解释人员提供更为可靠的结构曲率属性。
请参阅图1所示，本发明一种联合梯度结构张量和多窗分析的地层结构曲率估计方法，包括以下步骤：
步骤01：采集三维叠后地震数据并利用Hilbert变换获得每道地震信号的复地震道；
步骤02：在复地震道的基础之上计算得到三维瞬时相位数据体，同时计算瞬时相位的梯度向量；
步骤03：针对三维数据体的每个样点附近选取众多相邻分析窗；
步骤04：在每个相邻分析窗内构建相位梯度的协方差矩阵并对其进行特征分解以得到三个特征值及特征向量，利用主特征向量计算视倾角，同时利用三个特征值计算相似度；
步骤05：在所有相邻分析窗内选择相似度最大的分析窗，将其对应的视倾角指定为分析样点的视倾角；
步骤06：判断三维数据中的所有样点是否都计算得到了视倾角：若没有返回步骤03；若完成进入步骤07；
步骤07：计算三维视倾角数据体的偏导数；
步骤08：利用视倾角的偏导数计算地层的结构曲率。
步骤01中采集三维叠后地震数据并利用Hilbert变换获得每道地震信号的复地震道，具体包括：
若用f(x,y,t)表示三维地震数据，其中x和y表示三维地震数据在两个空间方向的指标，t表示三维地震数据在时间方向的指标。固定x和y，分别对每道信号f(x,y,t)进行Hilbert变换得到h(x,y,t)为：
其中τ为积分临时变量；则对应的复地震道c(x,y,t)为
c(x,y,t)＝f(x,y,t)+jh(x,y,t)＝A(x,y,t)ejφ(x,y,t)，
其中j为虚数单位，A(x,y,t)为瞬时振幅，φ(x,y,t)为瞬时相位。
步骤02中在复地震道的基础之上计算得到三维瞬时相位数据体并计算瞬时相位的梯度向量，具体包括：
由步骤01得到的复地震道c(x,y,t)的实部f(x,y,t)和虚部h(x,y,t)计算得到三维瞬时相位体φ(x,y,t)如下：
其中arctan表示反正切函数。由三维瞬时相位体φ(x,y,t)计算得到瞬时相位的梯度向量如下：
而φ(x,y,t)的偏导数由下列方法计算：
其中max表示取最大值，ε表示阻尼因子，取值较小(一般取0.001-0.01之间)。另外，计算f(x,y,t)以及h(x,y,t)的偏导数通过中心差分实现：
Δx，Δy及Δt为三维地震数据体在x，y及t方向上的采样间隔。
步骤03中针对三维数据体的每个样点附件选取众多相邻分析窗，具体包括：
假设当前的分析点为(x,y,t)，在其附近选取众多相邻分析窗。假设分析窗大小为Mx×My×Mt(其中Mx，My及Mt均为奇数，且Mx＝2Wx+1，My＝2Wy+1，Mt＝2Wt+1)，可以得到MxMyMt个相邻分析窗，且这些相邻分析窗均必须包含当前分析点(x,y,t)。假设分析窗大小为3*3*3，若用黑点表示当前分析点，则沿二维空间方向观察相邻分析窗如图2(a)所示，共有1个中心点分析窗和8个非中心点分析窗；沿空间-时间方向观察相邻分析窗如图2(b)所示，共有1个中心点分析窗和2个非中心点分析窗。综上，分析窗大小为3*3*3时，共有27个相邻分析窗(如图3所示1个中心点分析窗和26个非中心点分析窗，其中第1个为中心点分析窗)。
步骤04中在每个相邻分析窗内构建相位梯度的协方差矩阵并对其进行特征分解以得到三个特征值及主特征向量，利用主特征向量计算视倾角，同时利用三个特征值计算相似度，具体包括：
在步骤03得到的每个相邻分析窗内利用步骤02得到的相位梯度构建协方差矩阵，以中心点分析窗为例(即图3中的第一个窗)，构建当前分析点为(x,y,t)的相位梯度协方差矩阵C(x,y,t)：
其中mx，my及mt均为求和临时指标；
对相位梯度协方差矩阵C(x,y,t)进行特征值分解得到三个特征值(μ1，μ2及μ3，并有μ1≥μ2≥μ3)及对应的三个特征向量(v1，v2及v3)：
其中μ1及v1为协方差矩阵C(x,y,t)的主特征值及主特征向量。利用主特征向量v1的三个元素v1x，v1y及v1t估计在此相邻分析窗内地层沿x方向的视倾角px(x,y,t)及沿y方向的视倾角py(x,y,t)：
利用三个特征值μ1，μ2及μ3计算相似度similarity(x,y,t)如下：
在每个相邻分析窗内均计算视倾角和相似度。若分析点(x,y,t)共有M个相邻分析窗，记第m个相邻分析窗内估计的沿x及y方向的视倾角分别为px(m,x,y,t)及py(m,x,y,t)，记第m个相邻分析窗内估计的相似度为similarity(m,x,y,t)。
步骤05：在所有相邻分析窗内选择相似度最大的分析窗，将其对应的视倾角指定为分析样点的视倾角，具体包括：
分析点(x,y,t)周围共有M个相邻分析窗，选出相似度最大的相邻分析窗序号为win_opt(x,y,t)：
则分析点(x,y,t)处的沿x方向及y方向的视倾角分别为Px(x,y,t)及Py(x,y,t)：
Px(x,y,t)＝px(win_opt(x,y,t),x,y,t)，
Py(x,y,t)＝py(win_opt(x,y,t),x,y,t)。
步骤06：判断三维数据中的所有样点是否都计算得到了视倾角：若没有返回步骤03；若完成进入步骤07，具体包括：
对三维地震数据体中每一点均进行步骤3至步骤6的处理，可得到沿x方向的视倾角数据体Px(x,y,t)及沿y方向的视倾角数据体Py(x,y,t)。
步骤07中计算三维视倾角数据体的偏导数，具体包括：
针对沿x方向的视倾角数据体Px(x,y,t)分别计算其沿x方向偏导数Pxx(x,y,t)及y方向的偏导数Pxy(x,y,t)：
同样地，针对沿y方向的视倾角数据体Py(x,y,t)分别计算其沿x方向偏导数Pyx(x,y,t)及y方向的偏导数Pyy(x,y,t)：
上述求偏导数过程中，Δx和Δy为三维地震数据体在x和y方向上的采样间隔。
步骤08中利用视倾角的偏导数计算地层的结构曲率，具体包括：
为简化计算，利用步骤07中计算得到的三维视倾角数据体的偏导数Pxx(x,y,t)、Pxy(x,y,t)、Pyx(x,y,t)及Pyy(x,y,t)计算各点的各种曲率(最大正曲率MP(x,y,t)，最大负曲率MN(x,y,t))：
其中：
a(x,y,t)＝0.5×Pxx(x,y,t)，
b(x,y,t)＝0.5×Pyy(x,y,t)，
c(x,y,t)＝0.5×[Pxy(x,y,t)+Pyx(x,y,t)]。
以某油田的三维叠后地震数据为例，分别采用商业软件及本发明方法计算最大正曲率和最大负曲率。商业软件计算得到的最大正曲率和最大负曲率的时间切片如图4(a)及图5(a)所示，可以看出受噪声影响较大。利用本发明方法计算得到的最大正曲率及最大负曲率的时间切片如图4(b)及图5(b)所示，可以看出相对于商业软件的结果，本发明的结果反映地质结构更清晰且抗噪性能较好。
最后需要说明的是，以上算例对本发明的目的，技术方案以及有益效果提供了进一步的验证，这仅属于本发明的具体实施算例，并不用于限定本发明的保护范围，在本发明的精神和原则之内，所做的任何修改，改进或等同替换等，均应在本发明的保护范围内。