应用频谱分解进行薄储层预测，源自于1973年Widess^[1]提出的利用地震振幅信息，来定量计算储层厚度的方法。实际上地震波是由一定带宽内一系列不同频率的地震子波组成，根据地震波调谐原理，能与有效频带内某一频率产生调谐效应厚度的地层，在调谐体上均可得到一定的地震响应，较高频率的调谐效应与薄层相对应，反之则相反。频谱分解方法从原来的短时窗傅利叶变换(STFT)到小波变换(CWT)再发展到S变换^[2]，后来为了获得比小波变换和S变换更高的时频分辨率，发展了匹配追踪算法(Matching pursuit，简称MP)等非线性时频分析方法。MP算法于1993年首先由Mallat和Zhang^[3]提出，之后很快被应用于地震信号处理领域；2003年Castagna等^[4]将MP算法用于地震谱分解；2007年Wang^[5]将MP算法应用于地震资料的时频谱分解，随后对MP算法进行了改进，提出了多道MP算法，取得了很好的效果；2013年张繁昌^{[6, 7]}提出了基于Morlet小波字典的复数域双参数快速MP算法，进而减小了动态匹配搜索参数的个数，并提高了计算效率；2016年李坤^[8]基于Morlet小波快速动态匹配追踪高分辨率谱分解算法与体现流体因子的频变AVO反射特征方程相结合，发展高分辨率时频域FAVO直接反演方法，能够精细刻画油气藏位置。

自20世纪70年代，朝鲜在西朝鲜湾发现了资源量达30亿t的油气储量，该油气区也包括北黄海盆地东部坳陷^{[9, 10]}，随后在西朝鲜湾钻探了24口井，发现了侏罗系烃源岩，并在中-上侏罗统和下白垩统中获得了工业性油流^[11]。自1966年以来中国地质调查局在北黄海盆地进行大量的地质、地球物理和地球化学勘探，目前仍未获得油气的重大突破^[12]，油气勘探的主要难点是薄砂体储层预测。北黄海盆地的储层段为中-上侏罗统和下白垩统，主要储集相带为扇三角洲平原、扇三角洲前缘和近岸水下扇^[13]，但由于碎屑岩储层厚度薄，横向变化大，储层物性较差，非均质性较强，纵向上为多套薄砂层相互叠置，给本区的油气勘探带来很大的技术难度^{[14, 15]}。

前人对北黄海盆地的储层预测做过一些研究，王强^[14]、梁杰^[16]和王改云^[17]通过对单井粒度分析、薄片观察、扫描电镜分析和压汞测试等技术，认为北黄海储层总体分选为中-差，为低孔、低渗储层，储层空间有粒间溶孔、粒内溶孔、微裂缝等次生孔隙，沉积相带、成岩作用和埋藏史是控制本区砂岩储集层的主要因素；袁书坤^[18]结合地震、测井和录井资料的综合分析和研究，揭示下白垩统底部砂体是中生界最好的储层，砂岩粒度较粗，孔、渗性好，其底部直接与上侏罗统优质烃源岩接触；刘金萍^{[19, 20]}和万晓明^[21]运用岩石物理分析、岩性曲线重构声波时差和约束稀疏脉冲反演手段，对中生界储层进行预测，刻画了目的层中砂体的空间分布，另外，针对本区纵波阻抗叠置特点，采用叠前同步反演方法得到了纵横波速度比和密度参数，圈闭了孔隙度大于5%的砂岩为有利储层。刘俊^[22]通过连续小波变换，对北黄海东部坳陷下白垩统进行频谱分解，结合地震属性和水平切片技术，实现了在频率域进行砂体描述和储层横向变化规律的研究。吴淑玉^[15]通过对比几种时频分析方法，优选S变换对北黄海东部坳陷下白垩统地震资料进行频谱分解，有效地对储集层的沉积环境变化描述及沉积相带进行划分，说明时频分析方法在该区碎屑岩薄储层预测和识别的适用性。

针对频谱分解能够有效预测薄储层特点^[23]，本文根据研究区块内最新钻探的3口井，对下白垩统自然伽玛(GR)和自然电位(SP)测井曲线进行了沉积相划分，然后通过分析北黄海下白垩统薄层砂体的预测难点，对原始地震资料进行提频处理，将频带拓展到适合薄储层预测的范围，采用时频分辨率最高的MP算法进行频谱分解和地震相的沉积旋回分析，在高频层序地层上对分频属性切片进行沉积相分析，结合砂岩百分含量和储层物性资料，预测和圈定有利储层。

1.   区域地质背景

北黄海盆地位于黄海北部，盆地的北部为辽东-海洋岛隆起区，南部为胶北-刘公岛隆起区，西隔渤海口与渤海湾盆地相望，东邻安州盆地，面积约2.4×10⁴ km²，该盆地是印支-燕山运动时期在中朝板块东部(华北地块东延部分)区域隆起背景之上发育的以中、新生代为主陆内断陷盆地。盆地自东向西、由北往南可划分为东部坳陷、东部隆起、中部坳陷、中西部隆起、西部坳陷和南部凹陷群6个彼此分割的二级构造单元^[24](图 1)，其中东部坳陷位于北黄海盆地的东北部，总体构造为NE走向，是北黄海盆地最有油气远景的区域^[25-28]。

图  1  北黄海盆地构造区划图

(据文献[30]，箭头所指的蓝色边框为本文研究区块)

Figure  1.  Geotectonic map of the north Yellow Sea basin

(modified from reference[30], arrow aim at the blue frame is this study area)

下载: 全尺寸图片 幻灯片

晚侏罗世时期，由于太平洋板块沿NW方向向欧亚板块的俯冲，加上西部郯庐断裂的左旋走滑双重作用，整个盆地形成张扭性断陷。早白垩世时期，盆地处于干旱的亚热带区，气候干燥，盆地水体逐步变浅，地层减薄，岩性逐步变粗，盆地内发育了陆相红层和火山碎屑岩，在局部还出现短期性浅盐湖^[29]。白垩纪晚期，盆地由早期的张扭逐步转向压扭抬升，因区域抬升而露出水面，下白垩统首先遭受到区域剥蚀，抬升较高部位可剥蚀到上侏罗统，由于抬升剥蚀时间长，整个盆地缺失上白垩统地层。

根据岩心资料分析，下白垩统主要发育扇三角洲-滨湖相的沉积环境^[13](图 2)。扇三角洲平原由分支河道和分支河道间构成。以分支河道沉积相为代表的NYS1井位于河道中心，岩性以薄互层的灰色砂砾岩或浅灰色细砂岩与褐灰色泥岩相间，砂体发育好，分选为差-中等，可见交错层理、块状层理及冲刷-充填构造，GR和SP测井曲线为箱型，具有下粗上细的正旋回韵律。以分支河道间沉积相为代表的NYS2井位于河道边部，岩性由泥岩和粉砂岩组成，砂体减少、泥岩增加，发育水平层理，GR表现为薄层砂岩的指状特征，GR和SP测井测线为箱状锯齿形，具有下细上粗的反旋回韵律。扇三角洲前缘由水下分流河道、水下河道间、河口坝和席状砂构成，以NYS3井为代表，水下分流河道为灰色砂砾岩及浅灰色细砂岩，成分以石英为主，分选差-中等，GR呈钟形上细下粗的正旋回韵律；水下河道间以厚层的泥岩为主，间夹薄层的粉砂岩，GR为平直状，未见明显旋回性；河口坝以粉砂岩或细砂岩组成，GR呈漏斗型，可见反旋回，形成下细上粗的垂向层序；席状砂为粉砂或灰质粉砂岩，砂岩和泥岩呈薄互层，GR呈指状或锯齿状，未见明显的旋回性。

图  2  北黄海东部坳陷下白垩统地层联井对比分析

Figure  2.  Well constrast of the lower Cretaceous of the Eastern Depression of the North Yellow Sea Basin

下载: 全尺寸图片 幻灯片

北黄海盆地东部坳陷面积小，沉积物源搬运距离有限，储层岩性较为混杂，发育有砂砾岩、粗砂岩、细砂岩和粉砂岩等多种岩性，砂体储层整体较薄。白垩纪期间还发育大量的火成岩，储层的横向和纵向非均质性强，由于受机械压实作用和后期成岩改造作用的影响，储集层较为致密，物性普遍较差，下白垩统目的层埋藏较深，孔隙度一般在0~20%，平均在10%之内，属于低孔特低渗储集层。

2.   匹配追踪算法原理

谱分解技术是近年来发展起来的一种基于频率谱分解的储层特色解释技术，利用时频变换，将三维地震数据生成纵向上频率连续变化的振幅体，通过某一频率对应的调谐振幅的平面特点，描述储层的横向变化，能够更好地确定储层几何形态，改善地震分辨率^[31]。目前具有分辨率较高的频谱分解技术的是小波变换和MP算法，小波变换虽然克服了短时窗傅利叶变换的单分辨率分析的不足，通过引入了尺度因子a进行时频分析，但是由于尺度因子a与频率f没有直接的联系，而且频率在小波变换中没有明显表现出来，因此小波变换的结果不是一种真正的时频谱，且小波变换的另外一个缺点是其不具有自适应的特点，一旦基本小波被选定，就必须用它来分析所用的待分析数据。而MP算法的基函数是用超完备的冗余函数系统取代，在每一步迭代中选择最适应信号的近似部分的波形，在波束所有基内选择最佳适应全局特性的基，匹配追踪分离了与已给出的字典相一致的信号结构，这些信号结构特征信息的和是最能紧凑地表示原信号。

MP算法的基本思想是基于信号的可分解和重构，是在一个确定的函数集合中自适应地选择一些函数来表示一个信号的计算过程，函数集合中的每个函数都称为原子。MP算法是通过函数 f 垂直投影到函数集合D的元素上来进行重复的迭代估算，算法的基本表公式为：

式中[f, g_r0为第1次迭代时所选取的子波与原始地震信号的内积，R¹f为第1次迭代所产生的残差，其中R¹f和g_r0是正交的，满足以下关系：

为了使[f, g_r0] g_r0尽可能地接近原始信号f，必须使迭代最终产生的残差项足够小，使得内积项[ f, g_r0]足够大，即匹配追踪信号的实质是要找到一个与f最接近的原子。采用归纳法可解释匹配追踪如何得到进一步的应用。首先令R⁰f=f，假定已知第n次的余项Rⁿf，n≥0。由选择公式C，选择一个能够很好地匹配余项Rⁿf的元素g_γn∈D，

余项Rⁿf在下式中是一个子分解：

该式定义了第n+1次的余项，由于Rⁿ⁺¹f与g_γn正交，满足：

将这个分解一直进行到第m次，分解f为一个级联求和的形式：

由公式(4)得：

同样地，${\left\| f \right\|^2} $被分解为一个级联求和的形式：

由公式(5)得能量守恒公式：

原始矢量f被分解成被选取为最匹配其余项的字典元素的和。尽管此分解为非线性分解，但是保证能量守恒，假定其为正交线性分解。

而Morlet小波与地震子波具有相似性，适合对地震信号的能量和频谱进行分析，Morlet小波作为母小波创建子波字典的表达式为：

式中t为时间，τ为中心位置时间，φ为小波相位，f为主频，k为控制Morlet小波带宽的常数。原子的构建通过对不同的参数调控，包括对窗函数进行伸缩、平移和调制得到，尺度因子可以控制窗函数的伸缩(尺度因子越大，原子时宽越大；频宽越窄，频率分辨率越高)，时频和频率控制原子的平均时间和平均频率，通过三者相互作用使得信号在参数空间的能量分布更加集中。

基于Morlet小波的MP算法的步骤包括：(1)读取地震信号作为初始的剩余信号，将信号作为Hilbert变换，得到复地震道；(2)求取复地震道瞬时包络，确定包络中满足阈值条件的峰值；(3)将峰值处的时间，瞬时频率作为Morlet子波原子的初始值，进一步求取振幅和相位等信息；(4)建立改进的原子库，选取最佳原子；(5)计算最佳原子的时频分布；(6)从当前剩余信号中减去所选Morlet小波原子的实部；(7)重复上述步骤，直到满足预定的收敛条件。

3.   实际应用效果分析

3.1   地震提频处理

北黄海下白垩统薄层单砂体厚度大多为5~20 m^[32]，研究目的层深度段为1 600~2 800 m，地层速度为4 000~5 000 m/s，地震数据频带范围为8~60 Hz，主频为18 Hz，根据ΔZ=ν/(4×f)计算，可识别的最薄砂体厚度为25 m，通过对原始地震进行合成地震记录发现，薄层段砂体在地震上无明显响应，一个波峰或者波谷里面包含有几套薄层砂体。如NYS3井2 950~2 962 m深度段的灰色砂砾岩，与下段的砂泥岩互层，在地震反射上对应一个波谷(图 3)。一般情况下，利用叠后地震阻抗反演仅可识别砂岩厚度达几十米，因此，对于薄层的砂体，利用传统的地震阻抗反演方法都难以分辨单个砂体厚度，由图 4可以看出NYS3井3 041~3 062 m深度段21 m厚的河口坝薄层砂砾岩为荧光显示，而波阻抗反演显示与结果刚好相反，通过分析发现其上2 950~2 977 m段27 m厚的大套砂砾岩层段，阻抗值大，而3 041~3 062 m深度段薄层砂岩段含有一定的油气孔隙，阻抗值减小，因此强阻抗厚层砂砾岩对下面这套含油气薄层砂岩段有明显的岩性遮挡作用。

图  3  原始地震对薄层砂体合成地震记录的响应

Figure  3.  Thin sand body of the original seismic synthetic seismogram response

下载: 全尺寸图片 幻灯片

图  4  常规波阻抗反演体反映的碎屑岩砂体

Figure  4.  Conventional seismic inversion reflect the clastic sand body

下载: 全尺寸图片 幻灯片

频谱分解之前首先对地震资料进行脉冲反褶积提频处理，将地震的频宽拓展到8~80 Hz，提频后的地震带宽比常规地震数据带宽在高频端高出20 Hz，低频端保持不变，主频从原先的18 Hz拓展到25 Hz，拓频后的地震数据比常规地震细节丰富了很多。

对原始地震数据提频后，再进行分频处理，然后对高频段进行地震合成记录分析，可以看出地震响应能够反映薄层砂岩段的特征(图 5)，说明通过提频后的谱分解手段能够解决比原始地震主频更高的薄地层。

图  5  地震数据提频后高频段对薄层砂体的合成地震记录响应

Figure  5.  Thin sand body of the high frequency seismic synthetic seismogram after frequency-broadening processing

下载: 全尺寸图片 幻灯片

3.2   薄层MP算法频谱分解

将提频后的地震数据按照5~80 Hz的频率范围，频谱间隔为5 Hz进行MP算法进行频谱分解的剖面图(见图 6)，可以看出，20 Hz的频谱剖面在下白垩统目的层反应较为微弱，主要在下白垩统下部的地层有反应，30 Hz的频谱剖面在下白垩统反应最强烈，40 Hz也有所反应，但是相对30 Hz的频谱剖面较弱，代表着下白垩统的砂体以20~40 m厚度为主，50 Hz的频谱剖面显示较弱，但是在局部也能够显示，说明存在岩性较细的砂岩。通过对NYS2井进行MP算法频谱分解的地震相沉积旋回分析(图 7)，可以看出下白垩统存在2个上升半旋回和3个下降半旋回，粗砂体主要在下白垩统的下半部分，为扇三角洲砂体快速堆积的结果。根据等时沉积序列模式，将下白垩统地层按照旋回特征划分为5个沉积序列，从下到上分别为SSⅠ、SSⅡ、SSⅢ、SSⅣ和SSⅤ，针对每个沉积旋回进行频谱调谐切片，从平面上分析地层岩性厚度对调谐频率的振幅响应。

图  6  地震剖面分频显示

Figure  6.  Difference frequency profile

下载: 全尺寸图片 幻灯片

图  7  NYS2井点处的MP时频分析图

Figure  7.  Matching pursuit time- frequency analysis diagram in well NYS2 location

下载: 全尺寸图片 幻灯片

通过调谐数据体的振幅切片能够获得内部调谐细节，比常规全频带的振幅切片更能有效地展示不同的沉积特征^[15]。对提频地震数据进行频段扫描，得到一系列的频率域的调谐振幅数据体。图 8为下白垩统SSII层序内部的20、30、40和50 Hz的频率体均方根属性体切片(属性已归一化)，其中红黄色代表砂体，蓝紫色代表泥质。20 Hz的均方根属性平面图显示NYS1井砂体显示最为强烈，根据岩性资料，NYS1井的岩性最粗，最厚的细砂岩单层有62 m厚，位于河道边缘的NYS2和NYS4为弱反射，而位于扇三角洲前缘的NYS3井和NYS5井由北向南均方根能量逐渐降低，砂体逐渐由粗变细；30 Hz的均方根能量显示最强，沿着NYS4井、NYS2井和NYS3井分布呈强能量指状分布，由NYS3井和NYS5井往西南的扇三角洲前缘也呈现出砂体环指状减薄；40 Hz的砂体呈指状分布，NYS3井指示能量最强，单层砂体最厚为15 m，而40 Hz能够分辨18.75 m厚左右的砂体，NYS1井、NYS4井和NYS2井为弱能量；50 Hz整体能量较弱，NYS2井和NYS3井有所显示，指示着两口井存在薄层砂体。

图  8  SSⅡ层序分频地震数据的均方根属性切片

Figure  8.  RMS attributes slice map with different frequency seismic of SSⅡ sequence

下载: 全尺寸图片 幻灯片

根据表 1，下白垩统的扇三角洲砂岩百分含量最大值在NYS1井(65.9%)处，此处为靠近物源的扇三角洲平原沉积相，其中分支河道的以巨厚的细砂岩、粉砂岩夹薄层的泥岩为主，物性分选较差，河道底部为几米到十几米厚的砂砾岩滞留沉积；NYS2井砂岩百分含量为21.7%，距离物源稍远的扇三角洲平原沉积，由于河道的多次冲刷和充填，分支河道和河道间不断变迁，表现在巨厚泥岩间夹薄层砂岩；NYS3井为扇三角洲前缘沉积，砂岩百分含量为25.4%岩性较细，由于坡度减缓加上水体的承托力，流速减慢，下切作用减弱，河道变宽而浅，厚度减小，同时由于可容空间增大，侧向侵蚀作用相对变强。

表  1  下白垩统扇三角洲沉积的地层厚度和砂岩百分含量

Table  1.  Thickness and sandstone percentage of the lower Cretaceous fan delta deposits

井名	NYS1	NYS2	NYS3	NYS6	NYS7
厚度/m	195	235	337	103	151
砂岩百分比/%	65.9	21.7	25.4	10.22	11.3

下载: 导出CSV 

| 显示表格

结合时频属性切片可以看出砂体的分布由NW向SE逐渐发散，呈现扇状，推断物源方向在研究区的西北边，由西北到东南，水体逐渐加深，岩性由粗变细，能量从高能向低能转换，说明沉积环境从扇三角洲平原过渡到扇三角洲前缘，往南为前三角洲沉积。根据岩性分析，从古隆起上NYS1井到斜坡上的NYS2井为扇三角洲平原沉积，骨架以分支河道为主，而沉积中心的NYS3井为扇三角洲前缘沉积，骨架以水下分流河道为主(图 2)，由此说明频谱分解结果与录井结果是相吻合的。

3.3   时频分析的油气检测

根据取心井段下白垩统层段282块岩心样品分析，得出岩石物性数据表(表 2)，孔隙度最大值在NYS1井处。总体来看下白垩统平均孔隙度6.7%~14.5%，渗透率0.73~21.5 mD，从储层物性资料来看，研究区的下白垩统以低孔特低渗储层为主，而且储层物性随深度降低趋势，中生代地层(2 500 m以下)孔、渗条件普遍变差，储层类型为孔隙型和裂缝型^[33]。

表  2  下白垩统砂岩测井解释物性数据统计表

Table  2.  Petrophysical date of lower Cretaceous sandstone

井号	深度/m	孔隙度/%	渗透率/mD
NYS1	2 292.98~2 302.17	0.2~21	0.000 1~30
NYS2	2 333.15~2 624.97	0.2~15	0.01~7.8
NYS3	2 960.7~3 084.44	0.2~12	0.000 1~46.2
NYS4	2 229.6~2 301.8	5.1~17.3	0.3~31.5

下载: 导出CSV 

| 显示表格

地震波穿透地下介质时会发生衰减，经过含油气地层时，主频向低频方向移动，高频能量迅速衰减(图 9a)，由此可以采用时频分析进行油气检测^[34]。储层相对不含油气的围岩对地震波有更强的吸收作用，而且频率越高吸收越强，因此通过高频衰减能够指示地层是否含油气。图 9b为对地震数据体经过S变换得到的30 Hz的离散能量体(左)和10~80 Hz的频率道集(右)。可以看出，NYS2井2 295~2 385 m为9套含油气细砂岩，其中2 323~2 351 m日产油9.6 m³，频率道集在垂向上频率迅速衰减，低频能量相对较强；2 515~2 520 m和2 589~2 592 m段的含油气细砂岩，日产油4.14 m³，同样表现为高频迅速衰减，但是衰减的幅度没有顶部的层段快，说明含气层段比较薄，谱分解都能够与实钻结果吻合。

图  9  NYS2井谱分解进行烃类检测

(a)含油气频谱示意图；(b)谱分解成果图

Figure  9.  Spectrum analysis of NYS2 well for hydrocarbon detection

(a) Sketch of oil and gas spectrum; (b) Spectrum analysis results

下载: 全尺寸图片 幻灯片

4.   结论

(1) 北黄海盆地东部坳陷下白垩统为薄砂体、横向变化大和非均质性的储层特征，利用传统波阻抗反演方法难以区分有利储集层，而频谱分解技术是寻找薄储层的有利手段。

(2) 通过对地震数据进行拓频处理，采用匹配追踪算法进行频谱分解和地震相沉积旋回分析，在高频层序上对薄砂体储层进行分频解释，有效地进行薄砂体识别和烃类检测。

(3) 根据频谱分析结果，结合砂岩百分含量和孔渗等资料分析，下白垩统虽为低孔特低渗储层，但以指状分布的扇三角洲平原砂岩百分含量高，且烃类检测效果比较好，为下白垩统砂岩储集层圈闭的有利区域。