一种基于三维地震资料表征的地下火山通道识别方法
技术领域
本发明属于地震资料综合解释及地质分析，具体涉及一种基于三维地震资料表征的地下火山通道识别方法。
背景技术
盆地的形成与演化过程往往伴随不同规模的火山活动，火山岩是盆地充填系列的重要组成部分，在多种类型盆地中均有发育。油气勘探实践表明，火山岩储层物性与常规储层差异明显，其受深度影响较小，因此在盆地深部也具备良好的油气储集条件，可形成油气藏，盆地深部火山岩具有良好的勘探前景[毛治国，朱如凯，王京红等.中国沉积盆地火山岩储层特征与油气聚集.特种油气藏，2015，22(5)]。
针对深部火山岩油气勘探，首先需要落实火山岩空间分布特征。有效识别火山机构，是落实火山岩空间分布的基础，而火山通道是火山机构的关键组成部分，因此，火山通道的有效识别意义重大，正确有效识别火山通道是落实火山机构和火山岩油气藏的关键。由于火山通道位于地下深处，且近似垂直于水平面，常规的几何属性如相干、均方根属性等均只能反映小范围时窗信息，往往存在较多的干扰项，多解性强，无法有效反映火山通道分布情况。在以往的工作中，火山通道的识别仅通过二维地震剖面人工观察识别，在一块三维地震资料中设置一定的道间距逐条剖面观察，找出火山通道分布，这种方法一方面需要耗费大量的人力，且有一定的多解性，另一方面可能造成火山通道识别遗漏且无法精细刻画通道结构特征。
因此，特别需要一种方法能快速准确地识别火山通道。
发明内容
本发明的目的是能够快速准确地识别火山通道。
为了实现上述目的，本发明提供一种基于三维地震资料表征的地下火山通道识别方法，包括：步骤1：对叠后三维地震数据体进行各项异性扩散滤波，获取扩散滤波后的地震数据体；步骤2：基于所述扩散滤波后的地震数据体，获取相似性数据体；步骤3：对所述相似性数据体进行二值化处理，获得二值属性数据体；步骤4：基于所述二值属性数据体，计算目的层以下的二值属性累计值；步骤5：对所述二值属性累计值进行增益计算，获得增益后属性；步骤6：搜索所述增益后属性，确定局部最值点集，通过所述局部最值点集验证地下火山通道预测结果，获得地下火山通道坐标集；步骤7：基于所述地下火山通道坐标集和所述二值属性数据体，获得地下火山通道点集空间分布。
优选的，所述步骤2包括：采用纵横比大于第一阈值的时窗计算所述扩散滤波后的地震数据体的相似性，获得相似性数据体。
优选的，所述步骤3包括：采用二值截断法对所述相似性数据体进行二值化处理，获得所述二值属性数据体。
优选的，所述步骤4包括：基于所述二值属性数据体，采用大于或等于第一预设时间的时窗计算目的层以下的二值属性累计值，获得二值体累计属性平面分布。
优选的，所述步骤5包括：采用滑动增强算法对所述二值属性累计值进行增益计算，如以下公式(1)所示：
其中，F(x,y)为增益后属性，L为滑动窗口边长的一半，V(x,y)为t0至t1时窗内二值属性累计值，V(x,y,t)为t0至t1时间内任一时刻t的二值属性平面切片，(x，y)为叠后三维地震数据体的道坐标。
优选的，所述搜索所述增益后属性，确定局部最值点集包括：步骤61：将所述增益后属性按搜索窗大小划分为m×n个网格，其中，m为横向格数，n为纵向格数；步骤62：分别搜索每个网格的局部最大值点，获得局部最值点列表，所述局部最值点列表包括每个网格的局部最大值点的线道号和增益后属性；步骤63：分别比较所述局部最值点列表中的每个局部最大值点的增益后属性与预设截断阈值，当所述局部最大值点的增益后属性小于所述预设截断阈值时，则删除所述局部最大值点，否则保留所述局部最大值点，形成局部最值点集。
优选的，在所述步骤62中，根据以下步骤搜索一个所述网格的局部最大值点：步骤621：针对所述网格，以当前搜索窗在所述网格内进行最大值搜索，获得第一最大值点；步骤622：以所述第一最大值点为中心构建新的搜索窗，以所述新的搜索窗在所述网格内进行最大值搜索，获得第二最大值点，若所述第一最大值点与所述第二最大值点重合，将所述第二最大值点作为所述网格的局部最大值点，否则以所述第二最大值点作为所述第一最大值点，重复所述步骤622。
优选的，所述通过所述局部最值点集验证地下火山通道预测结果，获得地下火山通道坐标集具体包括：基于所述局部最值点集的地震剖面，删除非火山通道坐标，获得地下火山通道坐标集。
优选的，所述步骤7包括：根据以下步骤获得所述地下火山通道坐标集中的每个坐标点的空间分布：以所述二值属性数据体的任意平面时间切片为基础，采用圆形表征任意平面时间切片上地下火山通道的等效聚合分布形态，以所述地下火山通道坐标点为中心构建搜索窗，在所述搜索窗内计算所述地下火山通道坐标点的数据汇聚中心及所述圆形的半径，获得所述地下火山通道点的空间分布。
优选的，通过以下公式(2)计算每个所述地下火山通道坐标点在所述搜索窗内所述任意平面时间切片中的等效聚合分布形态：
其中，(xt，yt)为在所述搜索窗内t时间切片时的数据汇聚中心，V(x,y,t)为(x，y)道在t时间切片时的二值属性值，s为搜索窗数值；
所述地下火山通道坐标点的数据汇聚中心为所述圆形的圆心，通过以下公式(3)计算所述圆形的半径：
其中，Rt为所述圆形的半径，即t时间切片的火山通道坐标点的等效聚合分布形态的半径。
本发明的有益效果在于：本发明的基于三维地震资料表征的地下火山通道识别方法针对深部三维地震数据处理可有效获得火山通道在目的地层以下的空间分布状况，准确地预测火山通道发育位置，有效刻画火山通道空间分布形态，快速准确地识别火山通道，为深层火山岩分布预测和火山岩油气藏勘探提供依据。
本发明具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施例中将是显而易见的，或者将在并入本文中的附图和随后的具体实施例中进行详细陈述，这些附图和具体实施例共同用于解释本发明的特定原理。
附图说明
通过结合附图对本发明示例性实施方式进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施方式中，相同的参考标号通常代表相同部件。
图1示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的流程图。
图2示出了原始地震资料。
图3示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的各项异性滤波扩算后地震资料。
图4示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的二值化处理的结果剖面示意图。
图5示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的计算目的层以下的二值属性累计值示意图。
图6示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的增益后属性示意图。
图7示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的火山通道平面分布示意图。
图8示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的火山通道局部平面分布示意图。
图9示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的火山通道空间分布剖面示意图。
具体实施方式
下面将更详细地描述本发明的优选实施方式。虽然以下描述了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。
根据本发明的一种基于三维地震资料表征的地下火山通道识别方法，包括：步骤1：对叠后三维地震数据体进行各项异性扩散滤波，获取扩散滤波后的地震数据体；步骤2：基于扩散滤波后的地震数据体，获取相似性数据体；步骤3：对相似性数据体进行二值化处理，获得二值属性数据体；步骤4：基于二值属性数据体，计算目的层以下的二值属性累计值；步骤5：对二值属性累计值进行增益计算，获得增益后属性；步骤6：搜索增益后属性，确定局部最值点集，通过局部最值点集验证地下火山通道预测结果，获得地下火山通道坐标集；步骤7：基于地下火山通道坐标集和二值属性数据体，获得地下火山通道点集空间分布。
具体地，叠后三维地震数据体是一般的地球物理处理流程后形成的成果或纯波数据体，其资料信噪比一般不高，反映不连续性信号反映较差。以叠后三维地震数据体为输入数据，利用专业地震解释软件对地震资料进行各项异性扩算滤波，能够提高目的层之下地震资料信噪比，获得扩散滤波后的地震数据体，扩散滤波后的地震数据体比原始资料清晰且能够反映横向大尺度不连续性，对扩散滤波后的地震数据进行处理获得相似性数据体，相似性体能够反映地震资料空间不连续性特征的地震数据体，其信号强弱代表地震资料的不连续新程度，对相似性数据体进行二值化处理，获得二值属性数据体，计算二值属性数据体目的层以下的二值属性累计值，对二值属性累计值进行增益计算，获得增益后属性，对增益后属性进行搜索，确定局部最值点集，通过局部最值点集验证地下火山通道预测结果，进而获得地下火山通道坐标集，基于地下火山通道坐标集和二值属性数据体，获得地下火山通道点集空间分布。
根据示例性的基于三维地震资料表征的地下火山通道识别方法针对深部三维地震数据处理可有效获得火山通道在目的地层以下的空间分布状况，有效刻画火山通道空间分布形态，快速准确地识别火山通道，为火山机构预测提供直接依据。
作为优选方案，步骤2包括：采用纵横比大于第一阈值的时窗计算扩散滤波后的地震数据体的相似性，获得相似性数据体。
具体的，通过相似性属性计算获得适合大尺度裂缝研究的相似性数据体，，采用纵横比大于第一阈值的时窗计算扩散滤波后的地震数据体的相似性，优选的，第一阈值为6。
作为优选方案，步骤3包括：采用二值截断法对相似性数据体进行二值化处理，获得二值属性数据体。
具体的，采用二值截断法对相似性数据体进行二值化处理获得反映(断裂、火山通道等)等不连续信息的二值属性体，采用95％-99％之间的某一截断阈值对相似性数据体进行二值化处理，获得二值属性数据体。例如，采用95％截断阈值表示信号能量自强至弱5％信号时为1，95％信号时为0。
作为优选方案，步骤4包括：基于二值属性数据体，采用大于或等于第一预设时间的时窗计算目的层以下的二值属性累计值，获得二值体累计属性平面分布。
具体的，对二值属性数据体，采用大于或等于第一预设时间的时窗计算目的层(火山岩发育层)以下的二值属性体平面累积不连续属性，获得二值累计属性平面分布，该属性反映目的层面以下地震资料的整体不连续性。优选的，第一预设时间为1.5秒，若目的层之下地震资料的时窗不足1.5秒将影响预测效果，不能采用本方法。
作为优选方案，步骤5包括：采用滑动增强算法对二值属性累计值进行增益计算，如以下公式(1)所示：
其中，F(x,y)为增益后属性，L为滑动窗口边长的一半，V(x,y)为t0至t1时窗内二值属性累计值，V(x,y,t)为t0至t1时间内任一时刻t的二值属性平面切片，(x，y)为叠后三维地震数据体的道坐标。
具体的，采用滑动平均增强算法增益二值属性体平面累积不连续属性，采用30-50道间距正(奇数道)方形窗口增强平面上不连续属性，具体算法如公式(1)所示，获得不连续性增益平面属性，即增以后属性。
作为优选方案，搜索增益后属性，确定局部最值点集包括：步骤61：将增益后属性按搜索窗大小划分为m×n个网格，其中，m为横向格数，n为纵向格数；步骤62：分别搜索每个网格的局部最大值点，获得局部最值点列表，局部最值点列表包括每个网格的局部最大值点的线道号和增益后属性；步骤63：分别比较局部最值点列表中的每个局部最大值点的增益后属性与预设截断阈值，当局部最大值点的增益后属性小于预设截断阈值时，则删除局部最大值点，否则保留局部最大值点，形成局部最值点集。
作为优选方案，在步骤62中，根据以下步骤搜索一个网格的局部最大值点：步骤621：针对网格，以当前搜索窗在网格内进行最大值搜索，获得第一最大值点；步骤622：以第一最大值点为中心构建新的搜索窗，以新的搜索窗在网格内进行最大值搜索，获得第二最大值点，若第一最大值点与第二最大值点重合，将第二最大值点作为网格的局部最大值点，否则以第二最大值点作为第一最大值点，重复步骤622。
作为优选方案，通过局部最值点集验证地下火山通道预测结果，获得地下火山通道坐标集具体包括：基于局部最值点集的地震剖面，删除非火山通道坐标，获得地下火山通道坐标集。
具体的，根据火山机构规模确定平面搜索窗尺度，搜索增益后属性，确定局部最值点集，该点集即为识别的火山通道坐标集。采用适当的搜索窗搜索增益后属性的局部最值点，搜索窗一般采用正方形窗口，窗口尺度根据单个火山机构规模确定，通常选取0.5-1Km，(约25-50道间距)，具体搜索流程及算法如下：
(1)将增益后平面属性按搜索窗大小划分为m×n网格，m＝Crossline/g，n＝inline/g，其中m为横向格数、n为纵向格数，g为搜索窗边长，平面网格数据超出范围处按0填充。
(2)逐个网格按以下算法搜索局部最大值，并获得局部最值点列表：
a.针对一个网格，搜索当前网格内的第一最大值点，获得第一最大值点线道号(x，y)；
b.以当前获得的第一最大值点为中心构建新的搜索窗，在新的搜索窗内进行最大值点搜索，将搜索的最大值点作为第二最大值点，若第二最大值点与第一最大值点重合则停止搜索当前网格，若不重合则将第二最大值点作为第一最大值点，重复执行步骤b；
c.将当前网格搜索到的第二最大值点线道号(x，y，z)添加至搜索列表中，z为改点属性值，若列表中已存在该点则不添加；
d.跳转至下一网格，进行下一个网格的最大值搜索，若当前网格为最后一个网格则结束局部最大值搜索流程；
(3)最值点截断筛选：选取适当的z值预设截断阈值，分别比较局部最值点列表中的每个局部最大值点的增益后属性与预设截断阈值，当最值点的增益后属性小于预设截断阈值，则删除该最值点，否则保留该最值点。预设截断阈值选取H/h×10％至H/h×40％，其中H为第一预设时间的时窗，h为地震数据纵向采样时间间隔。
基于局部最值点集的地震剖面验证火山通道预测结果，删除非火山通道坐标，获得地下火山通道坐标集。
作为优选方案，步骤7包括：根据以下步骤获得地下火山通道坐标集中的每个坐标点的空间分布：以二值属性数据体的任意平面时间切片为基础，采用圆形表征任意平面时间切片上地下火山通道的等效聚合分布形态，以地下火山通道坐标点为中心构建搜索窗，在搜索窗内计算地下火山通道坐标点的数据汇聚中心及圆形的半径，获得地下火山通道点的空间分布。
作为优选方案，通过以下公式(2)计算每个地下火山通道坐标点在搜索窗内任意平面时间切片中的等效聚合分布形态：
其中，(xt，yt)为在搜索窗内t时间切片时的数据汇聚中心，V(x,y,t)为(x，y)道在t时间切片时的二值属性值，s为搜索窗数值；
地下火山通道坐标点的数据汇聚中心为圆形的圆心，通过以下公式(3)计算圆形的半径：
其中，Rt为圆形的半径，即t时间切片的火山通道坐标点的等效聚合分布形态的半径。
具体的，围绕火山通道坐标集，以任意火山通道刻画为例，以属性数据体的任意平面时间切片为基础，采用圆形表征任意平面时间切片上地下火山通道的等效聚合分布形态，通过各平面时间切片上火山通道的等效聚合分布形态表征火山通道空间展布。等效聚合分布形态包括数据汇聚中心(圆心)和圆形的半径，以一个地下火山通道坐标点为中心构建搜索窗，在搜索窗内计算地下火山通道坐标点的数据汇聚中心及圆形的半径，采用公式(2)计算地下火山通道坐标点在搜索窗内任意平面时间切片中的数据汇聚中心，采用公式(3)计算圆形的半径，进而获得一个地下火山通道点的空间分布。通过上述方法分别获取地下火山通道坐标集中的每个坐标点的空间分布，对比实际地震剖面，识别的火山通道与地震剖面和地质认识吻合性好，识别效果明显，可以为火山机构分布预测提供依据，为火山岩油气藏勘探提供了重要支撑。
实施例
图1示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的流程图。图2示出了原始地震资料。图3示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的各项异性滤波扩算后地震资料。图4示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的二值化处理的结果剖面示意图。图5示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的计算目的层以下的二值属性累计值示意图。图6示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的增益后属性示意图。图7示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的火山通道平面分布示意图。图8示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的火山通道局部平面分布示意图。图9示出了根据本发明的一个实施例的基于三维地震资料表征的地下火山通道识别方法的火山通道空间分布剖面示意图。
如图1所示，根据本发明的一种基于三维地震资料表征的地下火山通道识别方法，包括：
步骤1：对叠后三维地震数据体进行各项异性扩散滤波，获取扩散滤波后的地震数据体；
例如，针对某探区的三维地震满覆盖约5400km2，纵向采样率2毫秒，道间距25米，工区发育二叠系火山岩，埋深大，岩性多解性强，火山机构识别较困难，如图2所示为上述探区的原始三维地震资料，资料信噪比低。对叠后三维地震数据体开展各项异性扩算滤波，各项异性扩算滤波后的地震资料信噪比明显提高，大尺度断点更加清晰，如图3所示。
步骤2：基于扩散滤波后的地震数据体，获取相似性数据体；
其中，采用纵横比大于第一阈值的时窗计算扩散滤波后的地震数据体的相似性，获得相似性数据体；
例如，针对上述三维地震资料获得的扩散滤波后的地震数据体，采用纵横比>6时窗计算扩散滤波后的地震数据体的相似性，获得相似性数据体。
步骤3：对相似性数据体进行二值化处理，获得二值属性数据体；
其中，采用二值截断法对相似性数据体进行二值化处理，获得二值属性数据体；
例如，针对上述三维地震资料获得的相似性数据体，采用97％-99％之间某一截断阀值获得二值属性数据体，二值化处理的结果剖面如图4所示。
步骤4：基于二值属性数据体，计算目的层以下的二值属性累计值；
其中，基于二值属性数据体，采用大于或等于第一预设时间的时窗计算目的层以下的二值属性累计值，获得二值体累计属性平面分布；
例如，针对上述三维地震资料获得的二值属性数据体，采用2秒时窗计算二值属性累积平面属性，二值属性累计值如图5所示。
步骤5：对二值属性累计值进行增益计算，获得增益后属性；
其中，采用滑动增强算法对二值属性累计值进行增益计算，如以下公式(1)所示：
其中，F(x,y)为增益后属性，L为滑动窗口边长的一半，V(x,y)为t0至t1时窗内二值属性累计值，V(x,y,t)为t0至t1时间内任一时刻t的二值属性平面切片，(x，y)为叠后三维地震数据体的道坐标；
例如，针对上述三维地震资料获得的二值属性累计值采用35-45道间距(奇数道)正方形窗口增强地震资料的整体不连续性，增益后属性如图6所示。
步骤6：搜索增益后属性，确定局部最值点集，通过局部最值点集验证地下火山通道预测结果，获得地下火山通道坐标集；
其中，搜索增益后属性，确定局部最值点集包括：
步骤61：将增益后属性按搜索窗大小划分为m×n个网格，其中，m为横向格数，n为纵向格数；
步骤62：分别搜索每个网格的局部最大值点，获得局部最值点列表，局部最值点列表包括每个网格的局部最大值点的线道号和增益后属性；
其中，根据以下步骤搜索一个网格的局部最大值点：步骤621：针对网格，以当前搜索窗在网格内进行最大值搜索，获得第一最大值点；步骤622：以第一最大值点为中心构建新的搜索窗，以新的搜索窗在网格内进行最大值搜索，获得第二最大值点，若第一最大值点与第二最大值点重合，将第二最大值点作为网格的局部最大值点，否则以第二最大值点作为第一最大值点，重复步骤622；
步骤63：分别比较局部最值点列表中的每个局部最大值点的增益后属性与预设截断阈值，当局部最大值点的增益后属性小于预设截断阈值时，则删除局部最大值点，否则保留局部最大值点，形成局部最值点集；
其中，通过局部最值点集验证地下火山通道预测结果，获得地下火山通道坐标集具体包括：基于局部最值点集的地震剖面，删除非火山通道坐标，获得地下火山通道坐标集；
例如，针对上述三维地震资料获得的增益后属性，根据资料中地区火山作用特征选取0.8-1km内某一数值设置搜索窗，搜索搜索窗内的局部最大值点，针对局部最大值点列表，采用H/h×25％至H/h×35％内某一截断阈值，筛选点集，获得局部最值点集，对局部最值点集的地震剖面对比最终获得24个火山通道分布位置，即地下火山通道坐标集，火山通道平面分布如图7所示，局部平面放大图如图8所示。
步骤7：基于地下火山通道坐标集和二值属性数据体，获得地下火山通道点集空间分布；
其中，根据以下步骤获得地下火山通道坐标集中的每个坐标点的空间分布：以二值属性数据体的任意平面时间切片为基础，采用圆形表征任意平面时间切片上地下火山通道的等效聚合分布形态，以地下火山通道坐标点为中心构建搜索窗，在搜索窗内计算地下火山通道坐标点的数据汇聚中心及圆形的半径，获得地下火山通道点的空间分布；
其中，通过以下公式(2)计算每个地下火山通道坐标点在搜索窗内任意平面时间切片中的等效聚合分布形态：
其中，(xt，yt)为在搜索窗内t时间切片时的数据汇聚中心，V(x,y,t)为(x，y)道在t时间切片时的二值属性值，s为搜索窗数值；
地下火山通道坐标点的数据汇聚中心为圆形的圆心，通过以下公式(3)计算圆形的半径：
其中，Rt为圆形的半径，即t时间切片的火山通道坐标点的等效聚合分布形态的半径；
例如，上述某探区的地下火山通道坐标集分别获取每个坐标点各平面时间切片上火山通道的等效聚合分布形态，即地下火山通道点的空间分布，火山通道空间分布剖面如图9所示。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。
