CN112444859A - 协同异变蚂蚁体的页岩储层断裂识别方法及系统 - Google Patents
协同异变蚂蚁体的页岩储层断裂识别方法及系统 Download PDFInfo
- Publication number
- CN112444859A CN112444859A CN201910794914.2A CN201910794914A CN112444859A CN 112444859 A CN112444859 A CN 112444859A CN 201910794914 A CN201910794914 A CN 201910794914A CN 112444859 A CN112444859 A CN 112444859A
- Authority
- CN
- China
- Prior art keywords
- attribute
- coherent
- ant
- data
- seismic
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 80
- 230000001427 coherent effect Effects 0.000 claims abstract description 133
- 238000005516 engineering process Methods 0.000 claims abstract description 61
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 22
- 239000003016 pheromone Substances 0.000 claims description 67
- 238000005070 sampling Methods 0.000 claims description 48
- 241000257303 Hymenoptera Species 0.000 claims description 44
- 230000003044 adaptive effect Effects 0.000 claims description 42
- 238000004364 calculation method Methods 0.000 claims description 27
- 230000000007 visual effect Effects 0.000 claims description 19
- 238000012545 processing Methods 0.000 claims description 14
- 230000009466 transformation Effects 0.000 claims description 13
- 238000003708 edge detection Methods 0.000 claims description 12
- 238000004590 computer program Methods 0.000 claims description 10
- 238000010606 normalization Methods 0.000 claims description 9
- 238000000638 solvent extraction Methods 0.000 claims description 7
- 230000033001 locomotion Effects 0.000 claims description 5
- 238000001228 spectrum Methods 0.000 claims description 5
- 238000003860 storage Methods 0.000 claims description 4
- 238000012804 iterative process Methods 0.000 claims description 3
- 238000010835 comparative analysis Methods 0.000 claims description 2
- 238000004422 calculation algorithm Methods 0.000 abstract description 36
- 238000001514 detection method Methods 0.000 abstract description 14
- 230000004927 fusion Effects 0.000 abstract description 13
- 238000011161 development Methods 0.000 abstract description 12
- 238000006243 chemical reaction Methods 0.000 abstract 1
- 206010017076 Fracture Diseases 0.000 description 115
- 208000010392 Bone Fractures Diseases 0.000 description 104
- 230000006870 function Effects 0.000 description 21
- 230000000694 effects Effects 0.000 description 16
- 230000008569 process Effects 0.000 description 16
- 239000011159 matrix material Substances 0.000 description 15
- 230000018109 developmental process Effects 0.000 description 11
- 238000010586 diagram Methods 0.000 description 8
- 238000012360 testing method Methods 0.000 description 8
- 230000015572 biosynthetic process Effects 0.000 description 7
- 238000005755 formation reaction Methods 0.000 description 7
- 239000003079 shale oil Substances 0.000 description 6
- 238000009826 distribution Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 239000011435 rock Substances 0.000 description 4
- 230000005012 migration Effects 0.000 description 3
- 238000013508 migration Methods 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000007480 spreading Effects 0.000 description 3
- 238000003892 spreading Methods 0.000 description 3
- YBJHBAHKTGYVGT-ZKWXMUAHSA-N (+)-Biotin Chemical compound N1C(=O)N[C@@H]2[C@H](CCCCC(=O)O)SC[C@@H]21 YBJHBAHKTGYVGT-ZKWXMUAHSA-N 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 238000005452 bending Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000012512 characterization method Methods 0.000 description 2
- 230000036039 immunity Effects 0.000 description 2
- 230000035772 mutation Effects 0.000 description 2
- 239000000243 solution Substances 0.000 description 2
- 239000004215 Carbon black (E152) Substances 0.000 description 1
- 230000035508 accumulation Effects 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000013142 basic testing Methods 0.000 description 1
- 229960002685 biotin Drugs 0.000 description 1
- 235000020958 biotin Nutrition 0.000 description 1
- 239000011616 biotin Substances 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 239000012141 concentrate Substances 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000008021 deposition Effects 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 239000005712 elicitor Substances 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 229930195733 hydrocarbon Natural products 0.000 description 1
- 150000002430 hydrocarbons Chemical class 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000002347 injection Methods 0.000 description 1
- 239000007924 injection Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000005192 partition Methods 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 230000004304 visual acuity Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
- G01V1/302—Analysis for determining seismic cross-sections or geostructures in 3D data cubes
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/004—Artificial life, i.e. computing arrangements simulating life
- G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/63—Seismic attributes, e.g. amplitude, polarity, instant phase
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/64—Geostructures, e.g. in 3D data cubes
- G01V2210/642—Faults
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Theoretical Computer Science (AREA)
- Biomedical Technology (AREA)
- Molecular Biology (AREA)
- Biophysics (AREA)
- Computational Linguistics (AREA)
- Data Mining & Analysis (AREA)
- Evolutionary Computation (AREA)
- General Health & Medical Sciences (AREA)
- Artificial Intelligence (AREA)
- Computing Systems (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Health & Medical Sciences (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提出了一种协同异变蚂蚁体的页岩储层断裂识别方法及系统，通过利用变模态分解、Hilbert变换、倾角扫描及相干体技术形成自适应变模态相干体技术，提高了相干检测的精度；还通过相干属性、边缘属性及曲率属性构成多属性融合启发因子，结合协同异变的迭代方式的智能蚂蚁算法形成协同异变蚂蚁体技术，实现页岩储层多尺度断裂系统智能一体化识别，在识别大尺度断裂的同时提高小尺度断裂测检测精度，进而指导页岩储层油气的勘探与开发。
Description
技术领域
本发明涉及地震几何属性解释技术领域，尤指一种协同异变蚂蚁体的页岩储层断裂识别方法及系统。
背景技术
目前国内外石油工作的重心之一页岩油气的勘探与开发。地下的复杂的断裂系统对页岩油气的形成、运移、成藏起着重要的控制作用，其分布和形态直接影响油气藏的识别和描述，因此，精确预测断裂展布形态对页岩油藏勘探非常重要。
在三维地震勘探中，对断裂的精细刻画主要借助于相干体技术和蚂蚁体技术。相干体技术是由Amoco石油公司于1995年提出，其是利用相邻道地震信号之间的相干性来描述地层、岩性等的横向非均匀性，特别对识别断裂以及了解与储集层特征密切相关的砂体展布等方面非常有效。当遇到地下存在断裂或某个局部区域地层不连续变化时，一些地震道的反射特征就会与其附近地震道的反射特征出现差异，导致地震道局部的不连续性。这样，通过检测各地震道之间的相干程度，即可检测出断裂或不连续变化的信息。
蚂蚁体技术是由Schlumberger公司推出的一种基于蚁群算法的断裂、裂缝自动分析和识别的智能系统。蚂蚁体技术大大提高了断裂解释的效率和精度，它是基于数据体本身进行分析，减少了人为主观因素对断裂解释结果的影响，克服了解释者对断裂解释的主观性。运用蚂蚁自动追踪技术进行构造精细解释，解释结果更加符合地质实际情况，能够对地下构造有客观的认识，大大减少了常规解释的时间。
地震相干体技术从20世纪90年代开始推广使用，是目前用于地震断裂和裂缝检测非常有用的属性之一，相干体技术原理是在三维地震数据体中，对每一道每一样点求得与周围数据相关性，相干值较低的点反映了反射波波形不连续性，与地质上的断裂或地层边界有关。
相干体技术至今主要有三代算法：
第一代(C1)算法是计算主测线、联络测线方向的相关系数，合成主联方向相关系数；
第二代(C2)算法是对任意多道地震数据计算相干，是一种基于水平切片或层位上一定时窗内计算方法；
第三代(C3)算法是本征值相干算法，它把多道地震数据组成协方差矩阵，应用多道特征分解技术求得多道数据之间的相关性。目前使用最多、效果较好的是C3相干算法，其计算流程如下：
首先，沿采样点所在位置抽取附近的道的数据组成一个以计算点所在道为中心的D矩阵，如式(1-1)所示：
其中，dij是地震道(i，j)的数据，i＝1,2,…,N为Line，j＝1,2,…,M为Crossline。针对实际需要，数据点长度N大于一个波长。
其次，计算矩阵D的协方差矩阵Cov，计算公式如式(1-2)所示：
最后，对协方差矩阵Cov求特征向量，找出其中最大的特征值λmax，并利用式(1-3)计算相干值：
其中，C为C3相干算法的相干值，λmax为协方差矩阵Cov最大的特征值，λ为协方差矩阵Cov所有的特征值。
前述提到的相干体技术的理论研究和实际应用都已十分成熟，虽然在一定的地质条件下，可以解决断裂刻画问题，但是它也存在一定的局限性：C3相干算法遍历整个三维数据体使用固定大小时窗，无法适应频率随着深度增加和方向改变会发生相应变化，且没有消除轻微扰动的影响，易会产生噪音和假象，进而导致相干体技术只适合反射特征明显、全区沉积连续性较好的地层，无法一体化识别不同尺度组合的断裂系统。
蚂蚁体技术是通过在地震数据体中播撒大量的“蚂蚁”，在地震属性体中发现满足预设断裂条件的断裂痕迹的“蚂蚁”将“释放”某种信号，召集其他区域的蚂蚁集中在该断裂处对其进行追踪，直到完成该断裂的追踪和识别。而其它不满足断裂条件的断裂痕迹将不进行标注。人工蚂蚁具有一定的记忆能力，能够记忆已经访问过的节点，而且人工蚂蚁的选择机制是按照一定的算法规律有意识的寻找最短路径，而非像真实蚂蚁一样盲目选择。
每只蚂蚁在寻路的过程中需要满足如下的特性：
1、每只蚂蚁在行走过了的路径上都会释放信息素(生物素)，而且信息素是不断挥发的；
2、蚂蚁访问下一点是一个伪随机状态，其选择的时候与蚂蚁当前行走路径的距离和当前蚂蚁视野范围内点有关；
3、蚂蚁的一次循环运动中不能去到已经访问过的点，此时就需要一个禁忌搜索索引表，将已经访问过的点索引加入到禁忌搜索表，蚂蚁不会去到禁忌搜索表中的点去。
蚂蚁体算法的迭代过程可以表述如下：初始化启发因子ηij、信息素τij；算法开始将一定数量的蚂蚁随机放入地震属性体中的某个数据点，将当前点加入禁忌搜索表；迭代过程中每只蚂蚁根据自身的位置按照概率选择公式随机选择下一个数据点，然后将走过的点加入禁忌搜索表，每次每只蚂蚁走过了所有的点为一次迭代过程，其中第k只蚂蚁的从i到j点的概率为：
其中，k为蚂蚁的序号；i为蚂蚁当前所在位置；j为蚂蚁下一路径所在位置；t为当前时刻；τij(t)为t时刻在节点i与j连线上残留的信息素；ηij(t)为由位置i到位置j的启发因子；α、β分别表征信息素和启发因子的相对重要程度；Ω表示k蚂蚁下一次能去的所有数据点的集合；s为集合Ω中的所有数据点。
当蚂蚁完成一次迭代后，信息素按照如式(1-5)进行更新：
其中，τij(t+n)为经过n次迁移后路径ij上的信息素；ρ为信息素的残留程度；1-ρ为信息素的挥发程度；Δτij为一次循环中留在路径ij上的信息素，具体的计算方式如式(1-6)。
其中，Q都是一个恒定值，通常表示一次完整迭代过程中蚂蚁信息素的总增加量；Lk则表示第k只蚂蚁本次寻路所走过的总路径长度，若该路径满足迭代终止要求则停止。否则，需用式(1-5)重新计算各路径的信息素浓度，进行第二轮的搜索，通过信息素的不断更新最终收敛于最优断裂路径上。
相比于相干体技术，蚂蚁体技术对断裂的识别精度有所提高，但是它也只能大尺度的断裂进行精准识别，对于小尺度的断裂则识别精度不高，无法识别页岩储层中地震同相轴只有轻微扭动的未断开的微小裂缝，不能满足勘探开发的地质需求，需要寻找多尺度组合的断裂系统一体化精准识别方法。
通常，可以根据断裂持续活动时间、规模等因素，将断裂划分出若干级别，不同尺度级的断裂控制着不同规模的油气分布特征。综上来看，现有的相干体技术和蚂蚁体技术得到的结果信噪比和分辨率的都比较低，对于大尺度断裂成像基本满足要求，对于小尺度断裂(微小裂缝)识别精度较低。
为此亟需一种可以实现多尺度断裂系统的一体化识别及精细刻画的页岩储层断裂识别方案。
发明内容
为解决上述问题，本发明提出了一种协同异变蚂蚁体的页岩储层断裂识别方法及系统，可以提高了相干检测的精度，协同异变蚂蚁体技术实现页岩储层多尺度断裂系统智能一体化识别，识别效果显著，分辨率高，抗噪性好，能精准刻画不同尺度的断裂展布特征，对于页岩气储层断裂组合、解释具有推广应用价值。
在本发明一实施例中，提出了一种协同异变蚂蚁体的页岩储层断裂识别方法，该方法包括：
获取页岩储层的叠后地震数据；
利用变模态分解处理所述叠后地震数据，进一步根据Hilbert变换得到每个地震采样点的视周期；
设定信噪比特征对应的权重系数，根据所述地震采样点的视周期及权重系数，得到每个地震采样点的自适应相干时窗；
按照所述自适应相干时窗对叠后地震数据进行自适应网格剖分，得到相干数据，进一步结合地层倾角扫描得到的倾角属性，构建自适应变模态相干体技术，得到三维地震相干属性；
根据所述相干数据，通过八方位边缘检测得到边缘属性；
根据所述叠后地震数据，通过二次曲面拟合得到最大正曲率属性；
将所述到三维地震相干属性、边缘属性及最大正曲率属性进行归一化处理，得到归一化后的相干数据、边缘属性及最大正曲率属性，并构建多属性启发因子；
设定蚂蚁种群数量，将归一化后的相干数据进行分块，在分块数据中散入指定个人工蚂蚁，通过伪随机搜索起点，开始追踪检测断裂信息；
根据所述多属性启发因子，蚂蚁体按照异步实现方式进行迭代，更新全局路径，根据所述全局路径更新信息素；
根据迭代终止条件终止所述蚂蚁体，根据所述信息素浓度进行断裂线刻画，得到页岩储层断裂数据。
在本发明一实施例中，还提出了一种协同异变蚂蚁体的页岩储层断裂识别系统，该系统包括：
数据获取模块，用于获取页岩储层的叠后地震数据；
视周期计算模块，用于利用变模态分解处理所述叠后地震数据，进一步根据Hilbert变换得到每个地震采样点的视周期；
自适应相干时窗计算模块，用于设定信噪比特征对应的权重系数，根据所述地震采样点的视周期及权重系数，得到每个地震采样点的自适应相干时窗；
相干属性计算模块，用于按照所述自适应相干时窗对叠后地震数据进行自适应网格剖分，得到相干数据，进一步结合地层倾角扫描得到的倾角属性，构建自适应变模态相干体技术，得到三维地震相干属性；
边缘属性计算模块，用于根据所述相干数据，通过八方位边缘检测得到边缘属性；
曲率属性计算模块，用于根据所述叠后地震数据，通过二次曲面拟合得到最大正曲率属性；
归一化模块，用于将所述到三维地震相干属性、边缘属性及最大正曲率属性进行归一化处理，得到归一化后的相干数据、边缘属性及最大正曲率属性，并构建多属性启发因子；
起点搜索模块，用于设定蚂蚁种群数量，将归一化后的相干数据进行分块，在分块数据中散入指定个人工蚂蚁，通过伪随机搜索起点，开始追踪检测断裂信息；
信息素更新模块，用于根据所述多属性启发因子，蚂蚁体按照异步实现方式进行迭代，更新全局路径，根据所述全局路径更新信息素；
断裂识别模块，用于根据迭代终止条件终止所述蚂蚁体，根据所述信息素浓度进行断裂线刻画，得到页岩储层断裂数据。
在本发明一实施例中，还提出了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现协同异变蚂蚁体的页岩储层断裂识别方法。
在本发明一实施例中，还提出了一种计算机可读存储介质，所述计算机可读存储介质存储有执行协同异变蚂蚁体的页岩储层断裂识别方法的计算机程序。
本发明提出的协同异变蚂蚁体的页岩储层断裂识别方法及系统通过利用变模态分解、Hilbert变换、倾角扫描及相干体技术形成自适应变模态相干体技术，提高了相干检测的精度；还通过相干属性、边缘属性及曲率属性构成多属性融合启发因子，结合协同异变的迭代方式的智能蚂蚁算法形成协同异变蚂蚁体技术，实现页岩储层多尺度断裂系统智能一体化识别，在识别大尺度断裂的同时提高小尺度断裂测检测精度，进而指导页岩储层油气的勘探与开发。
附图说明
图1是本发明一实施例的协同异变蚂蚁体的页岩储层断裂识别方法流程图。
图2是本发明一实施例的生成地震采样点的自适应相干时窗的步骤流程图。
图3是本发明一实施例的构建自适应变模态相干体技术的步骤流程图。
图4是本发明一实施例的计算边缘属性的流程图。
图5是本发明一实施例的蚂蚁追踪断裂线进行刻画的流程图。
图6是本发明一实施例的协同异变蚂蚁体的页岩储层断裂识别系统结构示意图。
图7A是本发明一具体实施例的实际地震剖面图。
图7B是利用现有相干体技术得到的剖面图。
图7C是利用本发明的自适应变模态相干体技术得到的剖面图。
图8A是利用现有相干体技术的沿层切片效果图。
图8B是利用本发明的自适应变模态相干体技术得到的岩层切片效果图。
图9A是本发明一具体实施例的F1函数的函数值与迭代次数的对应关系图。
图9B是本发明一具体实施例的F2函数的函数值与迭代次数的对应关系图。
图9C是本发明一具体实施例的F3函数的函数值与迭代次数的对应关系图。
图10是本发明一具体实施例的协同异变多属性融合蚂蚁体的剖面效果图。
图11是本发明一具体实施例的协同异变多属性融合蚂蚁体的岩层切片效果图。
具体实施方式
下面将参考若干示例性实施方式来描述本发明的原理和精神。应当理解，给出这些实施方式仅仅是为了使本领域技术人员能够更好地理解进而实现本发明，而并非以任何方式限制本发明的范围。相反，提供这些实施方式是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。
本领域技术人员知道，本发明的实施方式可以实现为一种系统、装置、设备、方法或计算机程序产品。因此，本公开可以具体实现为以下形式，即：完全的硬件、完全的软件(包括固件、驻留软件、微代码等)，或者硬件和软件结合的形式。
在页岩储层的地震解释中，断裂解释一直重点和难点，其直接控制油气藏圈闭，对注采开发有直接的影响作用；在断裂的附近往往伴生裂缝的发育，裂缝是页岩油气的形成、运移、成藏的主要通道；研究断裂的展布特征对页岩油气的压裂开发具有一定的知道作用。因此，开展断裂检测研究，提高断裂识别精度，对页岩储层的勘探与开发非常重要。目前三维地震相干体技术和蚂蚁体技术在一定程度上可以完成地下大尺度断裂展布及组合特征的刻画，但对于小尺度裂缝的分辨能力不足；此类小尺度的断距在地震剖面上表现，常为同相轴扭曲扰动，没有明显的错断，在钻井资料中也不易观察识别，常规方法识别难度较大。
本发明利用边缘检测属性、倾角属性、曲率属性等先验约束下，融合自适应变模态分相干体技术与协同异变蚂蚁体技术研究更高精度的多尺度断裂在不同尺度不同方向的一体化识别方法，识别大尺度断裂的同时，提高小尺度断裂测检测精度，进而指导页岩储层油气的勘探与开发。
根据本发明的实施方式，提出了一种协同异变蚂蚁体的页岩储层断裂识别方法及系统。
下面参考本发明的若干代表性实施方式，详细阐释本发明的原理和精神。
图1是本发明一实施例的协同异变蚂蚁体的页岩储层断裂识别方法流程图。如图所示，该方法包括：
步骤S1，获取页岩储层的叠后地震数据；
步骤S2，利用变模态分解处理所述叠后地震数据，进一步根据Hilbert变换得到每个地震采样点的视周期；
步骤S3，设定信噪比特征对应的权重系数，根据所述地震采样点的视周期及权重系数，得到每个地震采样点的自适应相干时窗；
步骤S4，按照所述自适应相干时窗对叠后地震数据进行自适应网格剖分，得到相干数据，进一步结合地层倾角扫描得到的倾角属性，构建自适应变模态相干体技术，得到三维地震相干属性；
步骤S5，根据所述相干数据，通过八方位边缘检测得到边缘属性；
步骤S6，根据所述叠后地震数据，通过二次曲面拟合得到最大正曲率属性；
步骤S7，将所述到三维地震相干属性、边缘属性及最大正曲率属性进行归一化处理，得到归一化后的相干数据、边缘属性及最大正曲率属性，并构建多属性启发因子；
步骤S8，设定蚂蚁种群数量，将归一化后的相干数据进行分块，在分块数据中散入指定个人工蚂蚁，通过伪随机搜索起点，开始追踪检测断裂信息；
步骤S9，根据所述多属性启发因子，蚂蚁体按照异步实现方式进行迭代，更新全局路径，根据所述全局路径更新信息素；
步骤S10，根据迭代终止条件终止所述蚂蚁体，根据所述信息素浓度进行断裂线刻画，得到页岩储层断裂数据。
以下将对前述每一步骤进行详细解释说明。
首先，结合步骤S1-步骤S4来看，先要构建自适应变模态相干体技术，针对常规相干体技术的局限性，本发明利用变模态分解、Hilbert变换，并结合地层倾角扫描改进得到自适应变模态相干体技术。
执行步骤S1，获取某一工区的叠后地震数据；
通过步骤S2及步骤S3，得到地震采样点的自适应相干时窗，具体过程如图2所示，包括：
步骤S201，将所述叠后地震数据中的原始地震信号经过变模态分解后，得到多个子模式信号。
步骤S202，将所述子模式信号进行Hilbert变换得到Hilbert谱，根据所述Hilbert谱得到每个地震采样点的视周期。
步骤S203，选取
至
作为自适应的相干时窗长度，其中T为地震采样点的视周期。
步骤S204，利用所述信噪比特征对应的权重系数，调节所述自适应的相干时窗长度，得到每个地震采样点的自适应相干时窗；其中，信噪比特征对应的权重系数与所述信噪比成正比。
通常，地震波随着深度加深，高频成分严重损失，频率会变得越来越低，因此，要计算采样点的时窗大小，不同的深度选择不同的时窗长度。变时窗长度利用Hilbert变换计算得到，但是Hilbert变换只有针对窄带信号，并且一个时间点只含有单个频率分量的情况下，通过Hilbert变换得到的瞬时频率才有意义，否则就失去了瞬时频率属性的物理意义。基于这个问题，可以通过经验模态分解算法，即将一个信号不同频率的波形或趋势逐级地分解得到一系列具有不同特征频率的窄带信号或固有模态函数，再利用Hilbert变换计算瞬时频率。由于这种Hilbert变换的每个固有模态函数不一定是正交的，会导致运算多解性出现，造成频率成分泄露，无法完成等价逆变换。为避免这种情况，本发明步骤S2采用变模态分解算法来改进计算。
变模态分解作为一种非递归的、完全自适应的方法，可以以数学上健壮的方式稀疏地分解信号。地震信号的变模态分解可以表示为一个有约束的变分问题，如式(2-1)所示：
其中，f(t)为原始地震信号；k＝1,2,…,K；K为分解的模式(信号)个数；uk(t)为分解的第k个离散子模式(信号)；
为uk(t)的解析形式；ωk为第k个模式(信号)对应的中心频率；i为虚数单位；αk为平衡参数。
利用拉格朗日乘子法可以将式(2-1)转化为无约束问题，如式(2-2)所示：
其中，λ是拉格朗日算子。
式(2-2)需要利用乘数的交替变方向法进行求解，可以得到三个未知参数的迭代更新，如式(2-3)至(2-5)所示。
其中，ω为角频率；n是迭代次数；
分别是uk(t)、f(t)、λ(t)的傅里叶变换；Ωk是模态分解频率的半平面。
经过变模态分解后得到一系列的子模式(信号)，再利用Hilbert(希尔伯特)变换计算Hilbert谱，进而求取当前采样点的视周期T，随后利用
至
作为自适应的相干时窗长度。当信噪比高时，时窗长度会靠近
计算效果更理想，而信噪比低时，时窗长度会靠近
计算效果更理想。为此，如式(2-6)，可以设置随信噪比变化的权重系数，在上述操作的基础上进一步调节时窗长度，得到每个地震采样点的自适应相干时窗，应用于后面步骤的相干计算：
其中，Tc(t)为地震采样点的自适应相干时窗长度；T(t)是根据变模态分解和Hilbert变换得到的采样点的视周期；κ是与信噪比成正比的权重系数。
通过步骤S4，构建自适应变模态相干体技术，得到三维地震相干属性，具体过程如图3所示，包括：
步骤S401，按照所述自适应相干时窗对叠后地震数据进行自适应网格剖分，剖分为M×N网格的相干数据；例如，5×5网格。
步骤S402，扫描某一采样点的倾角属性，使所述采样点所在道不动，周围的n道的叠后地震数据每道进行线性数据细分插值，经过上下移动与中心道对比分析，得到相关系数并选取相关系数最大的移动距离，所述采样点的倾角属性为所述相关系数最大的移动距离的方向，该方向指的是倾角方向。
步骤S403，将所述相干数据与倾角属性融入自适应相干体技术中，构建自适应变模态相干体技术，利用所述自适应变模态相干体技术得到三维地震相干属性。
从地质的角度上来看，可以将地层简化成微小平面的集合，通常用倾角和方位角来描述地层的产状。在地震处理过程中，同样可以用倾角和方位角来描述一个反射同相轴所体现的地质层面。因此，获取反射同相轴的倾角信息尤为重要。
反射界面倾斜或扭曲在地震记录采样点上体现为上下位置的相对移动，但是轻微的倾斜可能移动距离不够一个采样间隔的时间，这时候就需要对数据进行重采样，在拟合算法不够准确的情况下可能造成数据特征被破坏。本发明采用计算过程快速简单的基于信号相似度的计算方法，主要过程为：对于要计算的采样点所在道不动，其周围的n道地震数据每道先进行线性数据细分插值，然后上下移动与中心道对比分析，找出相关系数最大的一个移动距离，而这个点的倾角方向就是该相关系数最大的移动距离的方向，该倾角方向即为倾角属性。相关系数的具体表达为式(2-7)：
其中，X，Y代表不同倾角方向的时间的信号；ρxy为信号X与Y的相关系数；Cov(X,Y)为信号X与Y的协方差；Var(X)是信号X的方差；Var(Y)是信号Y的方差。
将自适应的相干时窗、倾角扫描属性融入C3相干体技术，构建得到高分辨率高精度的自适应变模态相干体技术。利用变模态分解算法得到的瞬时属性更准确，得到的自适应相干时窗可以提高相干算法的分辨率；倾角扫描在一定程度上可以降低构造起伏所引起的断裂相干假象，提高相干算法检测的精度和准确性。
进一步的，结合步骤S5、计算边缘属性，具体过程如图4所示，包括：
步骤S501，在所述相干数据的网格的8个方向上设置Sobel算子权值，分别得到8个方向的一阶导数；
步骤S502，分别将8个方向的一阶导数的最大值作为所述边缘属性。
边缘属性是指图像中局部强度变化最显著的部分，主要存在于目标与目标、目标与背景、区域与区域(RGB不同)之间，沿着边缘矢量方向图像变化剧烈，垂直矢量方向图像变化缓慢。
边缘检测不仅在图像分割中有重要作用，在地震解释分析上同样有很大作用。地震数据突变的地方边缘属性强，比如，断裂处、裂缝处、岩性突变处；虽然在断裂发育很强的地方边缘的方向是紊乱的，但是借鉴图像处理中边缘检测的思想来相同处理所得到地震属性一类的数据体，可以取得很好的应用效果。边缘属性主要是由数据的一阶导数或者二阶导数计算而来，且二阶导数对噪音敏感。针对地震数据有噪音干扰的特点，本发明选择一阶导数Sobel边缘检测算子计算地震属性的边缘值。
对于Sobel算子如果只计算X、Y两个方向，不能很好的包涵所有的边缘信息。为了提高检测精度，本发明在相干体数据的5×5网格上，采用八方向Sobel算子计算得到8个方向的一阶导数：
其中，Gs为第s个方向上的一阶导数；cij为相干数据网格中第i行第j列的相干属性值；
为第s个方向上第i行第j列的Sobel算子权值。
具体的，不同方向的Sobel算子权值矩阵可以表示如下：
其中，m1是0°方向的Sobel算子矩阵；m2是22.5°方向的Sobel算子矩阵；m3是45°方向的Sobel算子矩阵；m4是67.5°方向的Sobel算子矩阵；m5是90°方向的Sobel算子矩阵；m6是112.5°方向的Sobel算子矩阵；m7是135°方向的Sobel算子矩阵；m8是157.5°方向的Sobel算子矩阵。随后，将8个方向的一阶导数的最大值作为边缘属性输出，即：
E＝max{Gs},s＝1,2,...,8；(2-9)
其中，E为边缘属性。
结合步骤S6，计算最大正曲率属性，具体过程如下：
曲率属性描述的是地震数据体的几何变换，与地震反射体的弯曲程度相对应，其对褶皱和断裂的反映敏感，是用于地质体构造特征的有效手段。曲率在数学上用于度量曲线的弯曲程度，曲率属性是应用曲率方法来计算地质体在几何空间上的分布形态，定量的描述构造特征，从而实现对断裂等几何构造的有效识别。
按照不同的算法，曲率可以分为很多种。不同算法的曲率属性以及用不同的参数计算的某一种曲率属性，可以不同程度反映断裂、线性特征、局部形状等信息。曲率属性有最大曲率、最小曲率、平均曲率、高斯曲率、最大正曲率、最小负曲率、走向曲率、倾向曲率、等值曲率和弯曲度，其中最大正曲率属于法向曲率中最大的正曲率，该属性可以突出边界，方法层面中的断裂信息。本发明在计算曲率属性时，需要事先对曲面进行M×N(例如，5×5)网格剖分，使需要计算曲率的点位于网格中心，对该曲面的局部进行二次曲面拟合，拟合公式如下：
z(x,y)＝ax2+by2+cxy+dx+ey+f；(2-10)
其中，z(x,y)为曲面的深度；x、y分别为离散后的网格的横、纵坐标；a、b、c、d、e、f是二次曲面拟合系数。根据拟合结果，可以计算最大正曲率属性：
其中，K为最大曲率属性。
结合步骤S7，构成多属性启发因子，具体过程为：
为了使得不同尺度、量级的多属性公共构成蚂蚁体算法的启发因子，需要对三维属性相干属性、边缘属性和最大正曲率属性做归一化处理。利用归一化的多属性构成多属性启发因子，具体如式(2-12)所示：
η＝w1·(1-C)+w2·E+w3·K+；(2-12)
其中，η为多属性启发因子；C为归一化后的相干数据；E为归一化后的边缘属性；K+为归一化后的曲率属性；w1、w2、w3分别为三种属性对应的权系数。
经过步骤S7，启发因子不再是单一因素决定，而是由几项指示断裂的地震属性线性加权而来，能够在保留全局信息的同时，将蚂蚁体的最优解限定在各项合理范围内。
进一步的，结合步骤S8-步骤S10，通过协同异变多属性融合蚂蚁体技术，实现页岩储层的断裂识别。
在本发明中，蚂蚁体技术选择启发因子为主导因素(即α小于β，α、β分别表征信息素和启发因子的相对重要程度)，且利用的是经过前述步骤S7得到的相干属性、边缘检测属性与曲率属性共同构建多属性启发因子。
在相干体的基础上，蚂蚁体可以设置了多个蚁群，并按照异步实现、协同优化的方式进行迭代，每个蚁群采用逆序的方式从高维到低维分维计算，适应度最差的维度将被对应蚁群中的其他维度代替。同时，蚂蚁之间的运动相互不影响，只与路径上的信息素和启发因子有关，每迭代完成一次，其对应的信息素、个体最优、全局最优将会对应更新。此外，多属性启发因子是利用相干、边缘及曲率等多属性融合构建的。信息素在迭代过程中也可能发生突变、重置或者说蚂蚁个体的信息素认知发生重置，以此来增加算法处理高维问题的能力和全局搜索的能力。
具体的，结合步骤S8，通过分块伪随机选择起点，过程为：
常规蚂蚁体技术散布蚂蚁时往往随机选择起点，这样会引入很强的不确定性，如果蚂蚁在某一区域的数量很少或者没有，该区域的断裂将不会被检测到。因此，以二维情况为例，本发明将归一化后的相干数据按照一定的大小分块，在每个块的独立区域里边放入一只人工蚂蚁；
在块的独立区域中，利用伪随机方式选择起点，如式(2-13)所示：
其中，选择一随机数P，P∈(0,1)，当P＞Pi时，第i个点为这个区域的起点；Ci为第i个点的归一化后的相干数据。
另外，为了方便快速收敛，可以将该小块区域内的最佳相干值所在位置作为起点。
结合步骤S9，蚂蚁追踪断裂线进行刻画，详细过程如图5所示，包括：
步骤S901，在所述蚂蚁体在选定起点之后，且在达到终止条件结束之前，在所述块的独立区域内自动寻找下一点，下一点的概率公式如式(2-14)所示，这个概率只和信息素及启发因子有关：
其中，P(i,j)为下一点的概率；R为蚂蚁体在位置(i，j)能见度范围内所有点的集合；τa是位置(i，j)的信息素；ηβ是位置(i，j)的多属性启发因子；α、β分别表征信息素和多属性启发因子的相对重要程度；
步骤S902，选择一随机数k，k∈(0,1)，当P(i,j)＞k时，将对应的(i，j)点作为蚂蚁体搜索的下一点；
步骤S903，经过迭代过程，记录所述蚂蚁体的行走路径，所述路径是所述蚂蚁体刻画出来的页岩储层断裂线；
在迭代过程中，假如蚂蚁一直没有被强制终止，蚂蚁就会沿着它认为合理的地方走下去，在被迫终止以后如果行走的路径满足要求，记录下本次的行走路径，这条路径也就是蚂蚁刻画出来的断裂线。
步骤S904，根据所述蚂蚁体刻画出来的页岩储层断裂线，比对本次路径和历史最优或全局最优路径，更新信息素，其中，信息素的更新公式如下：
τn+1(i,j)＝Δτ+(1-ρ)·τn(i,j)；(2-15)
τn+1为n+1次迭代的信息素；τn为n次迭代的信息素；ρ为信息素的挥发比例；Δτ为信息素的增量。
通过异步实现的思路，信息素的更新也可以按单个蚂蚁个体进行。
在步骤S9中，信息素更新有一个默认的规则，增加的信息素要比挥发的信息素要少，整个数据中信息素的浓度会越来越低，但是断裂处的信息素浓度随着迭代次数的增加相对其他位置会越来越强。一次迭代完成之后，人工蚂蚁的路径可能不是最优，但是随着迭代的增加，蚂蚁的最优路径会无限接近真实断裂线。为了增加收敛速度和全局搜索能力，每一次迭代结束后，信息素的增加只对最优路径进行；对此，可以对最优路径的判断来确定信息素的增量，从而更新信息素，具体判断过程如式(2-16)所示：
其中，Lbest为蚂蚁体的最优路径长度。在实际操作时，为了避免算法陷入局部最优，需要限定信息素的范围[τmin，τmax]。
最后，结合步骤S10，根据迭代终止条件，蚂蚁体终止行走，进而根据信息素的浓度进行断裂线的刻画，从而得到页岩储层断裂数据，该数据可以清晰地展示断裂平面的连续组合关系，显示出大、小尺度断裂情况以及断裂剖面的情况。
在本发明提出的协同异变蚂蚁体的页岩储层断裂识别方法中，可以按照需要进行蚂蚁体参数的设置，包括：
每只蚂蚁能够行走的最大步数MaxStep；
蚂蚁的视野范围View；
原始活动边界Origin，表示最开始分块的大小；
非法步数Fstep，表示蚂蚁行走的当前步不满足断裂的要求，一般的判定准则为这一步能见度范围内的相干值都很大(此处指的不是断裂)；
合法步数Rstep，表示一次非法步之后最少的合法步数，行走次数小于最小次数则认为这一次寻路是不合理的；
最小步进数MinStep，表示经过一次迭代之后蚂蚁最少行走的次数；
行走次数大于MaxStep时，强制禁止蚂蚁继续行走；
当非法步数大于Fstep时，强制禁止蚂蚁继续行走；
当合法步数小于Rstep时，强制禁止蚂蚁继续行走；
当一次迭代完成后行走次数小于MinStep，判定本次路径不合法，不进行蚂蚁信息素的更新。
本发明提出的协同异变蚂蚁体的页岩储层断裂识别方法应用自适应变模态相干体技术和协同异变多属性融合蚂蚁体技术在三维地震数据中实现多尺度断裂系统的一体化识别，效果显著，分辨率高，抗噪性好，计算结果忠实于地震资料，能精准刻画不同尺度的断裂展布特征，对于页岩气储层断裂组合、解释具有推广应用价值。
需要说明的是，尽管在上述实施例及附图中以特定顺序描述了本发明方法的操作，但是，这并非要求或者暗示必须按照该特定顺序来执行这些操作，或是必须执行全部所示的操作才能实现期望的结果。附加地或备选地，可以省略某些步骤，将多个步骤合并为一个步骤执行，和/或将一个步骤分解为多个步骤执行。
基于同一发明构思，本发明还提出了一种协同异变蚂蚁体的页岩储层断裂识别系统，如图6所示，该系统包括：
数据获取模块601，用于获取页岩储层的叠后地震数据；
视周期计算模块602，用于利用变模态分解处理所述叠后地震数据，进一步根据Hilbert变换得到每个地震采样点的视周期；
自适应相干时窗计算模块630，用于设定信噪比特征对应的权重系数，根据所述地震采样点的视周期及权重系数，得到每个地震采样点的自适应相干时窗；
相干属性计算模块604，用于按照所述自适应相干时窗对叠后地震数据进行自适应网格剖分，得到相干数据，进一步结合地层倾角扫描得到的倾角属性，构建自适应变模态相干体技术，得到三维地震相干属性；
边缘属性计算模块605，用于根据所述相干数据，通过八方位边缘检测得到边缘属性；
曲率属性计算模块606，用于根据所述叠后地震数据，通过二次曲面拟合得到最大正曲率属性；
归一化模块607，用于将所述到三维地震相干属性、边缘属性及最大正曲率属性进行归一化处理，得到归一化后的相干数据、边缘属性及最大正曲率属性，并构建多属性启发因子；
起点搜索模块608，用于设定蚂蚁种群数量，将归一化后的相干数据进行分块，在分块数据中散入指定个人工蚂蚁，通过伪随机搜索起点，开始追踪检测断裂信息；
信息素更新模块609，用于根据所述多属性启发因子，蚂蚁体按照异步实现方式进行迭代，更新全局路径，根据所述全局路径更新信息素；
断裂识别模块610，用于根据迭代终止条件终止所述蚂蚁体，根据所述信息素浓度进行断裂线刻画，得到页岩储层断裂数据。
应当注意，尽管在上文详细描述中提及了协同异变蚂蚁体的页岩储层断裂识别系统的若干模块，但是这种划分仅仅是示例性的并非强制性的。实际上，根据本发明的实施方式，上文描述的两个或更多模块的特征和功能可以在一个模块中具体化。反之，上文描述的一个模块的特征和功能可以进一步划分为由多个模块来具体化。
基于同一发明构思，还提出了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现协同异变蚂蚁体的页岩储层断裂识别方法。
基于同一发明构思，还提出了一种计算机可读存储介质，所述计算机可读存储介质存储有执行协同异变蚂蚁体的页岩储层断裂识别方法的计算机程序。
为了对上述协同异变蚂蚁体的页岩储层断裂识别方法及系统进行更为清楚的解释，下面结合一个具体的实施例来进行说明，然而值得注意的是该实施例仅是为了更好地说明本发明，并不构成对本发明不当的限定。
以某一实际工区为例，利用现有相干体技术与本发明的自适应变模态相干体技术对实际工区数据进行处理，剖面结果显示如附图7A-图7C所示，沿层切片显示如附图8A及图8B所示。
图7A为实际地震剖面图，图7B为利用现有相干体技术得到的剖面图，图7C为利用本发明的自适应变模态相干体技术得到的剖面图；图中横坐标为CDP道集的道号，纵坐标为时间。
图8A为利用现有相干体技术的沿层切片效果图，图8B为利用本发明的自适应变模态相干体技术得到的岩层切片效果图；图中横坐标为测线的线号，纵坐标为CDP道集的道号。
对比后可以看出(图7A-图7C可以着重参考图中圈出的部分)，本发明自适应变模态相干体技术在每个采样点处自动计算视频率，获取视周期进行相干时窗的选择，相比于现有相干技术，本发明检测效果更精确、更贴近实际，且展现更多的细节，可以为后面的蚂蚁体追踪多尺度断裂系统识别做准备。
接下来，利用现有蚂蚁算法与本发明的协同异变蚂蚁算法对三个标准测试函数进行理论试验，如表1及表2所示，表1是三个标准测试函数的具体参数，表2是现有蚂蚁体技术和本发明的协同异变蚂蚁体技术的基本参数设置。
试验时，标准函数的维度都设置为300，迭代次数规定为1000，每个标准函数利用不同的蚂蚁算法各自行40次独立试验，试验统计的最优值、最差值、平均值及标准差如表3所示。
表1标准测试函数
| 基本测试函数 | 代号 | 取值范围 | 最优值 |
| Sphere | F1 | [-5.12,5.12] | 0 |
| Weighted sphere | F2 | [-5.12,5.12] | 0 |
| Rotated hyper-ellipsoid | F3 | [-65.536,65.536] | 0 |
表2蚂蚁算法基本参数
| 算法 | 参数这是 |
| 现有蚂蚁算法 | а＝1，β＝5，τ∈[0.001,0.999]，ρ＝0.02，C<sub>Ao</sub>＝20 |
| 本发明 | а＝1，β＝5，τ∈[0.001,0.999]，ρ＝0.02，C<sub>Ao</sub>＝20，C<sub>o</sub>＝10 |
表3理论试验结果
图9A为F1函数的函数值与迭代次数的对应关系图，图9B为F2函数的函数值与迭代次数的对应关系图，图9C为F3函数的函数值与迭代次数的对应关系图；其中，圈1(实线)及圈2(虚线)分别表示常规蚂蚁算法和本发明的算法。
结合表2、表3及图9A至图9C可以清晰看到，相比于现有蚂蚁算法，本发明的协同异变蚂蚁算法在求解精度、收敛速度方面具有很大的优势，可以很大程度提高蚂蚁体技术的性能。结合多属性融合启发因子，利用协同异变多属性融合蚂蚁体技术检测实际地震数据，效果如图10及图11所示。
图10显示了协同异变多属性融合蚂蚁体的剖面效果图；图11显示了协同异变多属性融合蚂蚁体的岩层切片效果图。
与图7B、图8A相比(可以着重参考图中圈出的部分)，现有的相干体技术检测结果粒度较粗、连续性差，可以基本大致确定大尺度断裂发育的位置，但无法确定断裂组合样式，对小尺度断裂无法识别。而本发明提出的协同异变多属性融合蚂蚁体技术相对于现有相干体技术检测精度大幅度提高，断裂剖面检测结果与地震剖面相互吻合，断裂平面展示连续、组合关系清晰，大尺度断裂和小尺度断裂一体化显示，对比明显，因此可以证明该方法的有效性和先进性。
本发明提出的协同异变蚂蚁体的页岩储层断裂识别方法及系统通过利用变模态分解、Hilbert变换、倾角扫描及相干体技术形成自适应变模态相干体技术，提高了相干检测的精度；还通过相干属性、边缘属性及曲率属性构成多属性融合启发因子，结合协同异变的迭代方式的智能蚂蚁算法形成协同异变蚂蚁体技术，实现页岩储层多尺度断裂系统智能一体化识别，在识别大尺度断裂的同时提高小尺度断裂测检测精度，进而指导页岩储层油气的勘探与开发。
虽然已经参考若干具体实施方式描述了本发明的精神和原理，但是应该理解，本发明并不限于所公开的具体实施方式，对各方面的划分也不意味着这些方面中的特征不能组合以进行受益，这种划分仅是为了表述的方便。本发明旨在涵盖所附权利要求的精神和范围内所包括的各种修改和等同布置。
Claims (14)
1.一种协同异变蚂蚁体的页岩储层断裂识别方法，其特征在于，该方法包括：
获取页岩储层的叠后地震数据；
利用变模态分解处理所述叠后地震数据，进一步根据Hilbert变换得到每个地震采样点的视周期；
设定信噪比特征对应的权重系数，根据所述地震采样点的视周期及权重系数，得到每个地震采样点的自适应相干时窗；
按照所述自适应相干时窗对叠后地震数据进行自适应网格剖分，得到相干数据，进一步结合地层倾角扫描得到的倾角属性，构建自适应变模态相干体技术，得到三维地震相干属性；
根据所述相干数据，通过八方位边缘检测得到边缘属性；
根据所述叠后地震数据，通过二次曲面拟合得到最大正曲率属性；
将所述到三维地震相干属性、边缘属性及最大正曲率属性进行归一化处理，得到归一化后的相干数据、边缘属性及最大正曲率属性，并构建多属性启发因子；
设定蚂蚁种群数量，将归一化后的相干数据进行分块，在分块数据中散入指定个人工蚂蚁，通过伪随机搜索起点，开始追踪检测断裂信息；
根据所述多属性启发因子，蚂蚁体按照异步实现方式进行迭代，更新全局路径，根据所述全局路径更新信息素；
根据迭代终止条件终止所述蚂蚁体，根据所述信息素浓度进行断裂线刻画，得到页岩储层断裂数据。
2.根据权利要求1所述的协同异变蚂蚁体的页岩储层断裂识别方法，其特征在于，利用变模态分解处理所述叠后地震数据，进一步根据Hilbert变换得到每个地震采样点的视周期，包括：
将所述叠后地震数据中的原始地震信号经过变模态分解后，得到多个子模式信号；
将所述子模式信号进行Hilbert变换得到Hilbert谱，根据所述Hilbert谱得到每个地震采样点的视周期。
3.根据权利要求1所述的协同异变蚂蚁体的页岩储层断裂识别方法，其特征在于，设定信噪比特征对应的权重系数，根据所述地震采样点的视周期及权重系数，得到每个地震采样点的自适应相干时窗，包括：
选取
至
作为自适应的相干时窗长度，其中T为地震采样点的视周期；
利用所述信噪比特征对应的权重系数，调节所述自适应的相干时窗长度，得到每个地震采样点的自适应相干时窗；其中，信噪比特征对应的权重系数与所述信噪比成正比。
4.根据权利要求1所述的协同异变蚂蚁体的页岩储层断裂识别方法，其特征在于，按照所述自适应相干时窗对叠后地震数据进行自适应网格剖分，得到相干数据，进一步结合地层倾角扫描得到的倾角属性，构建自适应变模态相干体技术，得到三维地震相干属性，包括：
按照所述自适应相干时窗对叠后地震数据进行自适应网格剖分，剖分为M×N网格的相干数据；
扫描某一采样点的倾角属性，使所述采样点所在道不动，周围的n道的叠后地震数据每道进行线性数据细分插值，经过上下移动与中心道对比分析，得到相关系数并选取相关系数最大的移动距离，所述采样点的倾角属性为所述相关系数最大的移动距离的方向；
将所述相干数据与倾角属性融入自适应相干体技术中，构建自适应变模态相干体技术，利用所述自适应变模态相干体技术得到三维地震相干属性。
5.根据权利要求4所述的协同异变蚂蚁体的页岩储层断裂识别方法，其特征在于，根据所述相干数据，通过八方位边缘检测得到边缘属性，包括：
在所述相干数据的网格的8个方向上设置Sobel算子权值，分别得到8个方向的一阶导数；
分别将8个方向的一阶导数的最大值作为所述边缘属性。
6.根据权利要求5所述的协同异变蚂蚁体的页岩储层断裂识别方法，其特征在于，根据所述叠后地震数据，通过二次曲面拟合得到最大正曲率属性，包括：
对所述叠后地震数据的曲面进行M×N网格剖分，使需要计算曲率的点位于网格中心，对所述曲面进行二次曲面拟合，拟合公式如下：
z(x,y)＝ax2+by2+cxy+dx+ey+f；
其中，z(x,y)为曲面的深度；x、y分别为离散后的网格的横、纵坐标；a、b、c、d、e、f是二次曲面拟合系数；
根据拟合结果，利用下述公式计算最大正曲率属性：
其中，K为最大曲率属性。
7.根据权利要求6所述的协同异变蚂蚁体的页岩储层断裂识别方法，其特征在于，将所述到三维地震相干属性、边缘属性及最大正曲率属性进行归一化处理，得到归一化后的相干数据、边缘属性及最大正曲率属性，并构建多属性启发因子，其中，构建多属性启发因子的公式为：
η＝w1·(1-C)+w2·E+w3·K+；
其中，η为多属性启发因子；C为归一化后的相干数据；E为归一化后的边缘属性；K+为归一化后的曲率属性；w1、w2、w3分别为三种属性对应的权系数。
8.根据权利要求7所述的协同异变蚂蚁体的页岩储层断裂识别方法，其特征在于，设定蚂蚁种群数量，将归一化后的相干数据进行分块，在分块数据中散入指定个人工蚂蚁，通过伪随机搜索起点，开始追踪检测断裂信息，包括：
将归一化后的相干数据按照一定的大小分块，在每个块的独立区域里边放入一只人工蚂蚁；
在块的独立区域中，利用伪随机方式选择起点，如下式所示：
其中，选择一随机数P，P∈(0,1)，当P＞Pi时，第i个点为这个区域的起点；Ci为第i个点的归一化后的相干数据。
9.根据权利要求8所述的协同异变蚂蚁体的页岩储层断裂识别方法，其特征在于，根据所述多属性启发因子，蚂蚁体按照异步实现方式进行迭代，更新全局路径，根据所述全局路径更新信息素，包括：
所述蚂蚁体在选定起点之后，且在达到终止条件结束之前，在所述块的独立区域内自动寻找下一点，下一点的概率公式为：
其中，P(i,j)为下一点的概率；R为蚂蚁体在位置(i，j)能见度范围内所有点的集合；τa是位置(i，j)的信息素；ηβ是位置(i，j)的多属性启发因子；α、β分别表征信息素和多属性启发因子的相对重要程度；
选择一随机数k，k∈(0,1)，当P(i,j)＞k时，将对应的(i，j)点作为蚂蚁体搜索的下一点；
经过迭代过程，记录所述蚂蚁体的行走路径，所述路径是所述蚂蚁体刻画出来的页岩储层断裂线；
根据所述蚂蚁体刻画出来的页岩储层断裂线，比对本次路径和历史最优或全局最优路径，更新信息素，其中，信息素的更新公式如下：
τn+1(i,j)＝Δτ+(1-ρ)·τn(i,j)；
τn+1为n+1次迭代的信息素；τn为n次迭代的信息素；ρ为信息素的挥发比例；Δτ为信息素的增量。
10.根据权利要求9所述的协同异变蚂蚁体的页岩储层断裂识别方法，其特征在于，根据所述蚂蚁体刻画出来的页岩储层断裂线，比对本次路径和历史最优或全局最优路径，更新信息素，还包括：
在每次迭代结束后，对最优路径进行的信息素的增加，信息素的增量的公式如下：
Lbest为蚂蚁体的最优路径长度。
11.根据权利要求1-10中任一项所述的协同异变蚂蚁体的页岩储层断裂识别方法，其特征在于，该方法还包括：设置蚂蚁体参数，所述蚂蚁体参数包括：
每只蚂蚁能够行走的最大步数MaxStep；
蚂蚁的视野范围；
原始活动边界，表示最开始分块的大小；
非法步数Fstep，表示蚂蚁行走的当前步不满足断裂的要求；
合法步数Rstep，表示一次非法步之后最少的合法步数；
最小步进数MinStep，表示经过一次迭代之后蚂蚁最少行走的次数；
行走次数大于MaxStep时，强制禁止蚂蚁继续行走；
当非法步数大于Fstep时，强制禁止蚂蚁继续行走；
当合法步数小于Rstep时，强制禁止蚂蚁继续行走；
当一次迭代完成后行走次数小于MinStep，判定本次路径不合法，不进行蚂蚁信息素的更新。
12.一种协同异变蚂蚁体的页岩储层断裂识别系统，其特征在于，该系统包括：
数据获取模块，用于获取页岩储层的叠后地震数据；
视周期计算模块，用于利用变模态分解处理所述叠后地震数据，进一步根据Hilbert变换得到每个地震采样点的视周期；
自适应相干时窗计算模块，用于设定信噪比特征对应的权重系数，根据所述地震采样点的视周期及权重系数，得到每个地震采样点的自适应相干时窗；
相干属性计算模块，用于按照所述自适应相干时窗对叠后地震数据进行自适应网格剖分，得到相干数据，进一步结合地层倾角扫描得到的倾角属性，构建自适应变模态相干体技术，得到三维地震相干属性；
边缘属性计算模块，用于根据所述相干数据，通过八方位边缘检测得到边缘属性；
曲率属性计算模块，用于根据所述叠后地震数据，通过二次曲面拟合得到最大正曲率属性；
归一化模块，用于将所述到三维地震相干属性、边缘属性及最大正曲率属性进行归一化处理，得到归一化后的相干数据、边缘属性及最大正曲率属性，并构建多属性启发因子；
起点搜索模块，用于设定蚂蚁种群数量，将归一化后的相干数据进行分块，在分块数据中散入指定个人工蚂蚁，通过伪随机搜索起点，开始追踪检测断裂信息；
信息素更新模块，用于根据所述多属性启发因子，蚂蚁体按照异步实现方式进行迭代，更新全局路径，根据所述全局路径更新信息素；
断裂识别模块，用于根据迭代终止条件终止所述蚂蚁体，根据所述信息素浓度进行断裂线刻画，得到页岩储层断裂数据。
13.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至10任一所述方法。
14.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有执行权利要求1至10任一所述方法的计算机程序。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910794914.2A CN112444859A (zh) | 2019-08-27 | 2019-08-27 | 协同异变蚂蚁体的页岩储层断裂识别方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910794914.2A CN112444859A (zh) | 2019-08-27 | 2019-08-27 | 协同异变蚂蚁体的页岩储层断裂识别方法及系统 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112444859A true CN112444859A (zh) | 2021-03-05 |
Family
ID=74741275
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910794914.2A Pending CN112444859A (zh) | 2019-08-27 | 2019-08-27 | 协同异变蚂蚁体的页岩储层断裂识别方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112444859A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110333285A (zh) * | 2019-07-04 | 2019-10-15 | 大连海洋大学 | 基于变分模态分解的超声兰姆波缺陷信号识别方法 |
| CN115932967A (zh) * | 2023-01-06 | 2023-04-07 | 东北石油大学三亚海洋油气研究院 | 一种基于蚂蚁体算法的裂缝识别方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103675913A (zh) * | 2012-09-04 | 2014-03-26 | 中国石油天然气集团公司 | 一种裂缝型储层检测方法 |
| WO2016041189A1 (zh) * | 2014-09-19 | 2016-03-24 | 杨顺伟 | 一种评价页岩气储层及寻找甜点区的方法 |
| CN108931814A (zh) * | 2017-05-26 | 2018-12-04 | 中国石油化工股份有限公司 | 一种基于多属性融合的基岩裂缝预测的方法 |
| CN109001801A (zh) * | 2018-07-30 | 2018-12-14 | 中国石油化工股份有限公司 | 基于多次迭代蚁群算法的断层变尺度识别方法 |
-
2019
- 2019-08-27 CN CN201910794914.2A patent/CN112444859A/zh active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103675913A (zh) * | 2012-09-04 | 2014-03-26 | 中国石油天然气集团公司 | 一种裂缝型储层检测方法 |
| WO2016041189A1 (zh) * | 2014-09-19 | 2016-03-24 | 杨顺伟 | 一种评价页岩气储层及寻找甜点区的方法 |
| CN108931814A (zh) * | 2017-05-26 | 2018-12-04 | 中国石油化工股份有限公司 | 一种基于多属性融合的基岩裂缝预测的方法 |
| CN109001801A (zh) * | 2018-07-30 | 2018-12-14 | 中国石油化工股份有限公司 | 基于多次迭代蚁群算法的断层变尺度识别方法 |
Non-Patent Citations (1)
| Title |
|---|
| 王鹏;龚伟;: "断裂裂缝识别技术在塔中卡塔克1区块的应用", 新疆地质, no. 1, 15 December 2013 (2013-12-15) * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110333285A (zh) * | 2019-07-04 | 2019-10-15 | 大连海洋大学 | 基于变分模态分解的超声兰姆波缺陷信号识别方法 |
| CN110333285B (zh) * | 2019-07-04 | 2021-07-27 | 大连海洋大学 | 基于变分模态分解的超声兰姆波缺陷信号识别方法 |
| CN115932967A (zh) * | 2023-01-06 | 2023-04-07 | 东北石油大学三亚海洋油气研究院 | 一种基于蚂蚁体算法的裂缝识别方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CA2279266C (en) | Method for determining barriers to reservoir flow | |
| Hale | Methods to compute fault images, extract fault surfaces, and estimate fault throws from 3D seismic images | |
| AU2008239658B2 (en) | Inverse-vector method for smoothing dips and azimuths | |
| US20100131205A1 (en) | Method for identifying and analyzing faults/fractures using reflected and diffracted waves | |
| US20070168133A1 (en) | Computer-based method for while-drilling modeling and visualization of layered subterranean earth formations | |
| US6487502B1 (en) | System for estimating the locations of shaley subsurface formations | |
| US10073190B2 (en) | Method and system for geophysical modeling of subsurface volumes based on computed vectors | |
| WO1997033184A1 (en) | Method for geophysical processing and interpretation using seismic trace difference for analysis and display | |
| CN111596978A (zh) | 用人工智能进行岩相分类的网页显示方法、模块和系统 | |
| CN106772592B (zh) | 绕射波聚焦能量的分析方法及装置 | |
| GB2401453A (en) | Method of determining the petrophysical properties of a borehole and the region surrounding it | |
| WO2012139082A1 (en) | Event selection in the image domain | |
| CN111506861B (zh) | 一种目的层有利区域裂缝强度计算方法 | |
| Lagos et al. | Microseismic event location using global optimization algorithms: An integrated and automated workflow | |
| CN111399048A (zh) | 一种对断溶体计算相关属性及数据加权重构的方法 | |
| US6850864B2 (en) | Method for analyzing dip in seismic data volumes | |
| EP4337993A1 (en) | Method and system for seismic imaging using s-wave velocity models and machine learning | |
| CN112444859A (zh) | 协同异变蚂蚁体的页岩储层断裂识别方法及系统 | |
| Vasco et al. | Nonuniqueness in traveltime tomography: Ensemble inference and cluster analysis | |
| Ding et al. | Reliability analysis of seismic attribute in the detection of fault-karst | |
| Phelps et al. | Exploring viable geologic interpretations of gravity models using distance-based global sensitivity analysis and kernel methods | |
| CN114861515A (zh) | 层速度数据体的计算方法、装置、设备及介质 | |
| US11022710B2 (en) | Identifying anomalies in a subterranean formation based on seismic attributes | |
| US11346968B2 (en) | Estimation of reservoir flow properties from seismic data | |
| CN112444871A (zh) | 基于散射波地震响应特征的裂缝间距定量分析方法和设备 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |