基于改进的积分展开相位谱识别侧叠型砂体叠置区

An Improved Method for Identifying Side-Stack Type Sandbody Overlay Regions by Integral Expansion Phase Spectrum

摘要: 当前在河流相砂体叠置区域的识别方面，通常使用地震属性切片和基于地震沉积学的地层切片方法，或是通过瞬时属性和沿层平均属性等多种属性结合来识别叠置区域。对现有积分展开相位谱识别河道叠置区不连续性界限的方法进行改进，优化因分辨率与拾取方法导致提取位置不统一的时窗提取方法，对短时窗地震信号时间域补零从而细化相位频谱，以及在主值相位谱零频率处充零确保相位展开的稳定性。并用基于改进后的积分相位谱方法识别砂体叠置区，通过建立多种不同叠置情况的理论模型进行实验与胜利探区ZH地区实际资料测试。结果表明改进后的积分展开相位识别方法，时窗选取能更精确提取对应点位，展开相位谱细化且计算结果更准确，更稳定，积分展开相位谱识别叠置区域效果更好，通过实际资料验证该改进方法的有效性。

Abstract: In the identification of fluvial sandbody overlay areas at present, seismic attribute slicing and stratigraphic slicing based on seismic sedimentology are often used to identify overlapping areas, or through a combination of transient attributes and averages along the bed attributes. To improve the existing integral expansion phase spectrum identification fluvial facies sandbody overlay regions discontinuity boundary method, this paper optimizes the time window extraction method which had the problem that different time window positions are selected because of different resolution and picking methods, increases the phase resolution by adding zero to the time domain of short-time window seismic signal, and resets the zero frequency of the conventional phase spectrum to zero to ensure the stability of phase unfolding. The integral expansion phase spectrum method based on the improved method is used to identify the sandbody overlay area, as theoretical models of different superposition conditions were established to test the actual data in the ZH area of Shengli exploration area. The results show that the improved integral unexpansion phase recognition method can extract the corresponding points more accurately with time window selection, the resolution of the unfolded phase spectrum is higher, the calculation results are more accurate and more stable, and the integral expansion phase spectrum recognition effect is better. The effectiveness of the improved method is verified by actual data.

1. 引言

河流相砂体是我国陆相盆地重要的油气储集层，其沉积特征主要有：纵向上厚薄不均，砂泥互层，横向上连续性强弱不定，岩性变化快 [1] [2] [3] [4] 。这是由于河道的频繁改道与迁移，加之成岩改造及构造运动等因素，导致了河流相砂体易发育不同类型的砂体侧叠和堆叠，形成泥岩–砂岩–泥岩的薄互层叠置区 [5] [6] [7] [8] 。这种薄互层砂体的空间展布，以及这些叠置区域的识别是油气田开发地震的主要研究目标之一 [9] ，对于一个老油田的剩余油分布研究，以及开发方案的部署与调整都具有十分重要的意义。

目前在河流相砂体叠置区域的识别方面，通常使用地震属性切片和基于地震沉积学的地层切片进行砂体叠置的识别，或是通过瞬时属性和沿层平均属性等多种属性结合来识别叠置区域 [10] [11] [12] [13] 。而对于河流相砂体储层的油气藏开采，大都面对的是地震波调谐厚度以下的薄储层，这些薄储层为不同级次和不同尺度薄互层砂体的复合体，其间砂体具有复杂的叠置关系，在地震记录上通常都表现为一个周期的波形。对于这样一个周期的波形，常常利用其振幅、频率的变化来反映薄层厚度的变化，但对于两个河道砂体的叠置引起的地震波形变化，很难直接识别。井涌泉 [14] [15] 发现河流相薄互层砂体在地震剖面上产生的这些复合波的波形之间也存在着一些有规律的差异性特征，可能反映了薄互层砂体的相对厚薄及叠置关系。

在地震波相位的应用方面，主要是基于傅氏变换后的主值相位(2π周期内)和基于希尔伯特变换的瞬时相位及其展开相位，目前主要应用于地层厚度估算 [16] [17] 、地震层位追踪 [18] [19] ，以及断层、地层不整合面等识别 [20] [21] ，少有将相位应用到河道砂体间叠置区域的识别。刘阳 [22] ，王彬霓 [23] [24] 在总结前人对河流相不连续性界限预测方法与相位谱应用研究的基础上，考虑到地质结构变化往往更易于影响地震反射波的相位特征 [25] [26] ，提出了使用积分展开相位谱的方法来识别河流相砂体的不连续性界限，展现了积分展开相位谱在河流相叠置砂体识别方面的一定潜力，但是该方法受信号的时窗选取与噪音大小的影响较大。

本文对上述积分展开相位谱进行了目标层位时窗提取、短时窗地震信号时间域补零，以及相位谱频率域充零等多方面的改进实验，并通过理论模型与实际数据共同验证，改进后的方法识别河流相砂体的叠置区更准确有效。

2. 方法原理

2.1. 积分展开相位谱的计算方法

相位谱常见的两种表示方法分别是主值相位谱和展开相位谱。主值相位谱使用范围广，主要通过傅里叶变换或希尔伯特变换等直接求取。而展开相位谱是在主值相位谱的基础上进行求取，表现为对主值相位解包裹后的形态。

一般来说，地震记录 $x\left(t\right)$ 可以表示为地震子波 $b\left(t\right)$ 与反射系数 $r\left(t\right)$ 的褶积，即： $x\left(t\right)=b\left(t\right)\ast r\left(t\right)$ 。利用傅里叶变换可以得到其频谱，从而求取得到该信号不同频率的主值相位谱，即：

${\theta }_p\left({\omega }_i\right)=\arcsin \left(\frac{\mathrm{Im}\left\{X\left({\omega }_i\right)\right\}}{\sqrt{{\mathrm{Re}}^2\left\{X\left({\omega }_i\right)\right\}+{\mathrm{Im}}^2\left\{X\left({\omega }_i\right)\right\}}}\right)$ (1)

其中， $X\left(\omega \right)$ 为地震记录 $x\left(t\right)$ 的频谱， ${\omega }_i$ 为第i点的频率， ${\theta }_p\left({\omega }_i\right)$ 为主值相位谱，其范围为：

$-\pi \le {\theta }_p\left({\omega }_i\right)\le \pi$ (2)

通常在该范围内 ${\theta }_p\left({\omega }_i\right)$ 的值是不连续的，即相位在−π到+π之间折叠而呈周期性变化。因此，求取展开相位谱，就是将不连续的相位谱展开成一个连续的相位谱。首先确定展开相位的展开点位置，通过判断相邻频率采样点的相位差 $\Delta \theta \left({\omega }_i,\Delta \omega \right)$ 来确定展开点的位置：

$\Delta \theta \left({\omega }_i,\Delta \omega \right)={\theta }_p\left({\omega }_i+\Delta \omega \right)-{\theta }_p\left({\omega }_i\right)$ (3)

其中， $\Delta \omega$ 为相邻两个离散点的频率差，即 $\Delta \omega ={\omega }_{i+1}-{\omega }_i$ 。根据伊藤条件可知，相邻采样点的相位差应满足：

$-\pi <\Delta \theta \left({\omega }_i,\Delta \omega \right)\le \pi$ (4)

若相位差不满足伊藤条件时，认为此时为不连续点，从而对相位谱进行展开处理。为了消除不连续现象，通常在常规相位谱上加上适当的 $2k\pi$ 。即

${\text{e}}^{j{\tilde{\theta }}_p\left({\omega }_i\right)}={\text{e}}^{j\left[{\theta }_p\left({\omega }_i\right)+2\pi k\right]}$ (5)

其中， ${\tilde{\theta }}_p\left({\omega }_i\right)$ 为 ${\omega }_i$ 频率时的展开相位谱，k为展开相位的次数，依次从1开始增加。

将同一地震道不同频率的展开相位谱进行一个积分运算，来进一步突出不同地震道之间的展开相位谱的差异，称为积分展开相位谱：

${\theta }_H={\displaystyle \int \left|{\tilde{\theta }}_p\left(\omega \right)\right|\text{d}\omega }$ (6)

2.2. 时窗选取

信号时窗的选择决定是否能够最大程度提取到有用的地震信息，特别是对地震信号的相位计算有很大的影响。选取时窗过大，波形中含有其他层位信息，会干扰目的层位的属性解释结果；选取时窗过小，提取出的波形含有的层位信息小于计算需要的目的层位信息，导致信息不全，也会影响属性解释效果。

实验提取多种不同的时窗，包括顶底界面自定义的固定时窗、峰谷极值时窗、波谷时窗、波峰时窗、极大值时窗和周期时窗。发现周期时窗在地震相位属性计算时能够最大程度保留信号波形的信息，但是由于地震信号采样精度的影响，周期时窗提取的地震波形的零点实际上不存在于地震资料中，各种方法提取短时窗信号时，提取的起始点位置分布在“理论零点”附近，如图1所示，分别是“理论起点”前后的A、B点，以及“理论终点”前后的C、D点。

Figure 1. Schematic diagram of waveform extracted by time window

图1. 时窗提取的波形示意图

为了明确时窗提取时，“理论零点”周围不同的两点对时窗提取后，处理得到的相位谱的影响，将周期时窗起点设定为A、B两点，终点设定为C、D两点，由此组合成AD、AC、BD、BC四种不同的周期时窗提取模式。

Figure 2. Four types of wedge-shaped sand body models (yellow-sandstone, brown-mudstone)

图2. 四种楔形砂体模型(黄色砂岩，褐色泥岩)

将如图2所示的四种楔形砂体模型，模型设定各参数如表1所示，其中单砂体的最大厚度为调谐厚度(λ/4 = 11.7 m)。利用雷克子波进行褶积正演得到的合成地震记录后，根据上文描述的四种周期提取模式，分别提取该模型地震记录的周期时窗，处理得到各自的积分展开相位谱如图3所示，其中BD时窗影响最大，所有楔形砂体出现的区域，均表现为异常高值，难以区分叠置和非叠置区域；AC时窗的模型(4)部分在异常高值与背景低值之间不断变换，难以准确识别；对比分析发现，AD时窗与BC时窗这种表现出一定对称性的时窗提取模式下，AD时窗效果好于BC时窗，所以最终确定AD时窗作为一个周期时窗。

Table 1. Model parameter table

表1. 模型参数表

(a) AD时窗处理得到的积分谱 (b) AC时窗处理得到的积分谱(c) BD时窗处理得到的积分谱 (d) BC时窗处理得到的积分谱

Figure 3. Integral unwrapped phase spectrum obtained from post-processing of four different time windows extraction

图3. 四种不同时窗提取后处理得到的积分展开相位谱

2.3. 周期时窗提取方法改进

考虑到低信噪比下，噪音不仅仅通过影响波形变化来影响相位谱的计算结果，也会影响到传统周期时窗的提取，所以有必要对时窗提取方法进行改进，由传统的直接在一定范围内寻找波形正负变换的点，作为周期时窗的起始点，改进为如图4所示，先寻找较易追踪的峰谷极值层位，再以此为基点，在限制范围内寻找一个周期的起始点或结束点，如图5所示。

Figure 4. Peak valley extremum time window

图4. 峰谷极值时窗

Figure 5. Improved cycle time window

图5. 改进后的周期时窗

(a) 5%噪音的时窗提取图(改进前)(b) 5%噪音的时窗提取图(改进后)

Figure 6. Schematic diagram of cycle time window extraction before and after improvement under 5% random noise

图6. 5%随机噪音下改进前后的周期时窗提取示意图

(a) 10%噪音的时窗提取图(改进前)(b) 10%噪音的时窗提取图(改进后)

Figure 7. Schematic diagram of cycle time window extraction before and after improvement under 10% random noise

图7. 10%随机噪音下改进前后的周期时窗提取示意图

验证该时窗提取的改进方法是否更优，实验有效信号含5%、10%随机噪音下，分别利用改进前后的周期时窗提取方法，对含噪音的合成地震记录提取周期时窗，如图6(a)，图7(a)所示，为改进前的方法，无法有效提取出噪音影响下的部分道的一个周期波形，而图6(b)，图7(b)所示为改进后的方法，可以更加准确、有效地提取出对应的周期波形。

2.4. 时间域信号补零

(a) 不补零的AD时窗积分谱 (b) 补零101的AD时窗积分谱

(c) 补零401的AD时窗积分谱 (d) 补零1001的AD时窗积分谱(e) 补零1201的AD时窗积分谱 (f) 补零2001的AD时窗积分谱

Figure 8. Integral unwrapped phase spectrum obtained by processing seismic signals with different lengths of zeros extracted from time windows

图8. 时窗提取的地震信号补零不同长度后处理得到的积分展开相位谱

地震信号受单砂体厚度、叠置砂体层数以及其他影响，不同道之间通过时窗提取得到的地震信号长短不一，会对后续不同道相位对比产生影响，考虑使用时间域补零法消除该影响 [27] [28] 。

主值相位求取时，对拾取出的有限长度的地震数据，在数据尾端补零至合适的点数，首先使得不同周期长度的所有地震道的数据长度得到了统一。另外如图8所示，对时窗提取出的地震信号补零不同点数的对比实验发现，在补零不超过1001点时，随着补零点数的增加，积分谱的形状变化稍微变好，且突变值与背景值差距越来越大，更有利于该方法的识别；补零点数超过1001点后，补零点数增加反而会使积分谱的形状变差。

所有道补零至1001点时，利用积分展开相位谱识别叠置区域效果最好，此时依据采样率计算发现，地震记录对应的奈奎斯特频率为1000 Hz，即补零1001点对应频率域中的相位谱数据在每1 Hz频率上都有相应的相位变化数据。

2.5. 频率域零频率充零

因为时窗提取出的波形并非真正从零点开始，导致相位处理时，计算出的主值相位不同道的零频率处的相位往往会杂乱无章，部分道的零频率位置以0开始，部分道则跳变到π开始，且使得整道结果偏离一定位置，从而对后续展开谱和积分谱的识别造成影响。为了消除这种影响，选择频率域零频率处充零的方法实验。

如图9理论模型所示，设置多种不同叠置情况的模型，模型参数根据表1基本不变，仅修改单砂体最大厚度为λ/8，即使砂体叠置的最大累计厚度也不超过调谐厚度。利用雷克子波褶积后得到如图10所示的地震记录，按照改进后的时窗提取方法，提取图中红绿线标出的地震信号，计算得到主值相位谱后，分别按照零频率处充零与不充零两种方法处理得到其展开相位谱与积分展开相位谱。

Figure 9. A sand mudstone model with multiple overlapping relationships

图9. 多种叠置关系的砂泥岩模型

Figure 10. Schematic diagram of synthetic seismic records and time window extraction corresponding to multiple stacked sand and mudstone models

图10. 多种叠置关系的砂泥岩模型对应的合成地震记录与时窗提取示意图

根据理论模型可知，7~10道全部为砂体叠置区域，如图11所示，(a)、(b)分别是提取出不充零与充零下处理得到的展开相位谱中的7~10道。可以看到黑色的两道充零前后没有变化；绿色的两道充零前以π开始，变化趋势与黑色两道一致，但在开始位置出现一截异常增大区域，充零后异常增大区域消失，整体趋势充零前一致，且正好相差2π位置。所以充零并不会改变相位展开的整体趋势，而会改进一些杂乱无章的变化，增加相位计算结果的稳定性。

(a) 充零前7~10道展开相位谱(b) 充零后7~10道展开相位谱

Figure 11. Extract the comparison image of the unfolded phase spectrum of channels 7~10 before and after zero charging

图11. 提取出7~10道的展开相位谱充零前后对比图

如图12所示，对比展开相位谱，可以看到所有地震道的展开相位按照同等间隔分开后，从左往右黑、绿、粉、紫、绿、黑、红7个颜色的数据依次为图9的模型(1)~模型(7)对应的展开相位，充零前因为部分道从π开始，整道的相位数据位移2π，导致对比显得杂乱无章，充零后相位展开更稳定，所有道相位谱变得更整齐，数据整体趋势不变而所有道对比更加统一；对比积分展开相位谱，可以看到充零后背景道数据更加整齐，异变道数据更加突出，更有利于利用积分展开相位谱识别砂体叠置区域。

(a) 充零前的展开相位谱 (b) 充零后的展开相位谱(c) 充零前的积分相位谱 (d) 充零后的积分相位谱

Figure 12. Comparison diagram of unfolded phase and integrated phase obtained before and after zero frequency charging of the principal phase (5 * (n − 1) is added to all the data for each expanded phase, and n is the number of channels, so as to interval the data of each channel evenly. The different color regions in the unwrapping phase correspond to the model (1)~model (7) in Figure 9 from left to right)

图12. 主值相位零频率出充零前后得到的展开相位、积分相位对比图(展开相位每一道的所有数据都在原始数据的基础上加5*(n − 1)，n为道号数，以此将各道数据均匀的间隔开，展开相位中不同颜色区域从左往右依次对应图9的模型(1)~模型(7))

3. 模型测试分析

3.1. 正演模型建立

Figure 13. A sand mudstone model with multiple overlapping relationships

图13. 多种叠置关系的砂泥岩模型

在图9模型的基础上，增加理论模型的砂体叠置类型如图13所示，模型各参数按照表1所示，楔形单砂体的长度为50 m，最大厚度为λ/8，即砂体叠置的累计厚度最大为调谐厚度。使用雷克子波褶积正演得到的合成地震记录如图14所示，所有地震道的波形均因干涉形成单一的复合波，无法直接识别砂体厚度和叠置关系。

模型(1)两个楔形砂体反向完全叠置，最左侧10 m贴合叠置，实际上没有泥岩间隔层，可以看作厚度为8.19 m的非叠置砂体；模型(2)~(5)依旧是两个反向楔形砂体叠置，两砂体的叠置区域依次减少10 m，高程依次减少1.17 m，右侧砂体尖端的10 m始终保持与左侧砂体如模型(1)的贴合叠置，至模型(5)，因为两砂体接触位置无泥岩间隔层，可以看作一个连通的砂体；模型(6) (7)两个楔形砂体之间的泥岩间隔层始终保持相同厚度，模型(6)泥岩间隔为1.17 m，模型(7)泥岩间隔为2.34 m是模型(6)的一倍；模型(8)两楔形单砂体方向完全一致叠置，左侧10 m贴合叠置视作连通，右侧40 m上下砂体厚度完全相同。

Figure 14. The seismic response synthesized by the sand-mudstone model in Figure 13

图14. 图13砂泥岩模型合成的地震响应

3.2. 改进后积分展开相位谱分析

Figure 15. Figure 14 shows the integral unfolded phase spectrum obtained from the seismic response calculation

图15. 图14地震响应计算得到的积分展开相位谱

模型按照改进后的方法处理得到的积分展开相位谱如图15所示，积分谱表现出平稳的背景值，以及部分异常高值，对应模型的叠置区域，模型(1)~(4)砂体叠置区域变少，积分谱上异常高值区域也减少。部分叠置区域在积分谱上也表现为背景值，未能有效识别出来。

归纳分析发现，一部分叠置区域的两层薄砂体厚度相同且均小于调谐厚度的一半，推测对正演时地震波形干涉叠加形成复合波的影响具有对称性，从而使波形含有的相位信息变为单一砂体；另一部分叠加区域的下层薄砂体厚度还要小于泥岩间隔层厚度，在小于调谐厚度时，下层薄砂体对波形的影响过小，导致无法识别。

4. 实际工区应用

在实际资料应用方面，选取胜利探区ZH地区馆上段5砂组的地震资料来测试方法的准确性、适用性。该工区沉积环境复杂多变，以曲流河沉积为主，储集体以细砂岩为主，且厚度较薄，横向变化快，易形成多个砂体叠合构造。本文方法在该工区具有较好的应用价值。

为了验证改进后的积分展开相位谱方法，依据图16所示为ZH地区馆上段5砂组均方根振幅平面属性图，提取图中红线标出的CrossLine1096、1100和InLine1400三条不同测线，以及各自测线的馆上段5砂组下部砂体层位，分别选取图中黄线所示的，在地质人员与地震解释人员为工区结合地质、地震和测井资料在实际开发中已经验证了的河道砂体叠置区带。处理得到最终的积分展开相位谱。

Figure 16. RMS amplitude plane attribute of sand group 5 in Guanshang section of ZH area

图16. ZH地区馆上段5砂组均方根振幅平面属性图

如图17为CrossLine 1096测线的部分剖面，根据相关层位数据(图中红线所示)，提取出该剖面的一个周期时窗(图中绿线所示)进行处理。图18为解释人员根据工区地质资料、地震资料和测井资料确定的该位置两套砂体，绿色砂体与红色砂体在该区域出现侧叠型叠置，组成砂–泥–砂的薄互层叠置，与均方根振幅属性图中识别出的黄色叠置位置重合。处理得到的积分展开相位谱(图19)在该位置正好表现为突变，且周围区域没有过多干扰，有效识别出叠置区域。

Figure 17. CrossLine1096 line profile and window extraction map

图17. CrossLine 1096测线剖面及时窗提取图

Figure 18. CrossLine 1096 profile and sand body overlay diagram

图18. CrossLine 1096测线剖面及砂体叠置关系示意图

Figure 19. CrossLine 1096 line profile corresponding integrated phase spectra

图19. CrossLine 1096测线剖面对应积分相位谱

如图20为CrossLine 1100测线的部分剖面，根据相关层位数据(图中红线所示)，提取出该剖面的一个周期时窗(图中绿线所示)进行处理，该位置两层砂体的地震波形并未因为叠加而形成单一复合波，但是两砂体的波形还是因为叠加而产生变形，实际提取时窗为上层砂体受下层砂体产生形变的波形。图21为解释人员根据工区地质资料、地震资料和测井资料确定的该位置两套砂体侧叠区域，与均方根振幅属性图中识别出的黄色叠置位置重合。处理得到的积分展开相位谱(图22)在该区域表现为异常高值，其余区域为背景低值，有效识别出叠置区域。

Figure 20. CrossLine 1100 profile and time window extraction map

图20. CrossLine 1100测线剖面及时窗提取图

Figure 21. CrossLine 1100 profile and sand body overlay diagram

图21. CrossLine 1100测线剖面及砂体叠置关系示意图

Figure 22. CrossLine 1100 line profile corresponding integrated phase spectra

图22. CrossLine 1100测线剖面对应积分相位谱

图23为InLine 1400测线的部分剖面，蓝色为目标层位的层位数据，中间红圈标出位置，区域出现断层。图24为根据层位数据(图中红线所示)，提取出该剖面的一个周期时窗(图中绿线所示)进行处理。图25为解释人员根据工区地质资料、地震资料和测井资料确定的该断层右侧区域出现的多期叠置区域，与均方根振幅属性图中识别出的黄色叠置位置重合。处理得到的积分展开相位谱(图26)在该区域表现为连续的异常高值，其余区域为平稳的背景值，有效识别出叠置区域。

利用该工区三条横向、侧向测线对应含叠置区域的剖面，验证了不同叠置情况下，改进后的积分展开相位谱识别砂体叠置区域方法的有效性和适用性。

Figure 23. InLine 1400 survey line profile and fault area schematic diagram

图23. InLine 1400测线剖面及断层区域示意图

Figure 24. InLine 1400 line profile and window extraction map

图24. InLine 1400测线剖面及时窗提取图

Figure 25. Schematic diagram of InLine 1400 survey line profile and sand body stacking relationship

图25. InLine 1400测线剖面及砂体叠置关系示意图

Figure 26. Integrated phase spectrum corresponding to InLine 1400 line profile

图26. InLine 1400测线剖面对应积分相位谱

5. 结论

本文通过对积分展开相位识别砂体叠置区域的改进，取得一下效果：

1) 有效提高了提取地震波形一个周期时窗的抗噪性；

2) 提高了积分展开相位谱的计算结果的稳定性；

3) 提高叠置砂体识别的有效性和适用性。

如果进一步探究累计厚度小于调谐厚度时，部分叠置情况未能有效识别的原因，可以根据总结出的叠置砂体厚度与泥岩间隔层厚度关系，设置更多模型实验，从而找出不能识别的深层原因。

基金项目

四川省创新创业训练项目(S202210615160)资助。

参考文献

NOTES

^*第一作者。