东海陆架盆地西湖凹陷渐新统花港组年代标尺及层序界面定量识别

传统层序地层学研究主要是依据岩心、露头和生物化石，结合地震资料及测井资料进行不同级别的层序地层划分［1］。在缺乏详细岩心及生物化石或地球化学资料的情况下，肉眼识别测井曲线旋回尚存在一些多解性。通过井-震结合，充分运用数学手段对测井曲线进行分析，突显其旋回性和界面信息，可提高层序地层划分和识别的准确度及可信度［2］。傅里叶变换、沃尔什变换、小波变换、Hilbert‑Huang变换和最大熵谱分析技术等数学信号处理技术在层序地层划分、储层流体识别及储层岩性表征等方面进行了推广应用引发了广泛的关注［2-7］。傅里叶变换和沃尔什变换具有恒定窗口的平稳频率的特质；小波变换从不同尺度分析稳定和非稳定连续信号；Hilbert‑Huang变换是一种基于经验的数学分析方法，用于处理非线性和非稳定信号；最大熵谱分析是通过积分处理突出信号变化趋势［7-9］。这些数学分析方法成功地运用在四川盆地和塔里木盆地等层序研究中，显著地提高了层序地层划分的准确性。

西湖凹陷是东海陆架盆地最大的含油气性凹陷，约占东海陆架盆地油气资源的60 %［10］，具有巨大的勘探潜力，探明资源量为3 000 × 108～5 000 × 108 m3［10-11］，是近海天然气勘探重点区域之一［12-13］。花港组是西湖凹陷当前勘探开发的重点层位，近年来不断发现新的油气层［11，14］，其储量呈逐年快速增长的趋势。孙思敏等［15］通过岩心观察、测井与录井资料，将花港组划分为2个中期旋回；魏恒飞等［16］基于地震、钻测井及岩心资料综合分析的基础上，将花港组划分为2个长期基准面旋回和7个中期基准面旋回；张建培等［17］依据钻井地质、地球物理和分析化验等资料将花港组分为2套三级层序；高雁飞等［18］以层序地层学为指导依托测井数据和地震资料将花港组划分为3套三级层序；于兴河等［19］依据测井、岩心以及地震反射特征等资料将花港组划分为5套三级层序及12套四级层序；周瑞琦等［20］依托钻录井和地震等资料将花港组划分为3套三级层序。综上所述，花港组内部层序划分方案仍未统一，且层序叠置样式及时间尺度仍不清楚，严重制约有利相带预测、砂体等时对比及油气勘探进程。

本文以西湖凹陷中南部花港组为研究对象，基于层序地层学基本理论和方法，运用小波变换和最大熵谱分析技术，从测井曲线上获取反映沉积环境变迁的关键层序界面信息；结合地震、岩心等资料，逐级完成花港组层序界面定量识别并建立其“浮动天文标尺”，以期为花港组后续油气勘探提供地质理论基础。

西湖凹陷位于中国东海陆架盆地东北部，呈北东走向，东邻钓鱼岛隆起带，西接长江坳陷和海礁-渔山东隆起带，北靠富江凹陷，南连钓北凹陷；凹陷自西向东依次为西部斜坡带、西次凹、中央反转构造带、东次凹和东部断阶带，其面积约为5.18 × 104 km2；其中中央反转构造带自北向南分为嘉兴构造带、宁波构造带、玉泉构造带、黄岩构造带和天台构造带（图1a）。

研究区横跨西次凹和黄岩构造带，自下而上发育上始新统平湖组、渐新统花港组和中新统龙井组等。始新世晚期—渐新世早期发生玉泉运动，导致西湖凹陷西部和钓鱼岛褶皱带隆起及西湖凹陷中部快速沉降［21-22］；与此同时，全球发生变冷事件（Oi-1骤冷事件）［23］，海/湖平面快速下降使得部分平湖组露出水面遭受广泛的剥蚀［22］；渐新世晚期—早中新世早期，花港运动及全球变冷事件（Mi-1骤冷事件）导致湖平面下降［23］，造成花港组暴露遭受剥蚀，与上覆龙井组呈区域性角度不整合接触。28.1 Ma左右，花港组内部发生一定规模海退并造成高部位沉积间断；因此，花港组可分为花上段和花下段［19］。始新世西湖凹陷沉积物主要源于北部虎皮礁隆起，该物源大约贡献了总沉积物量的3/5；东部钓鱼岛隆起和西部海礁隆起为次要物源，大约贡献了总沉积物量的2/5［24］，该时期西湖凹陷主要发育河流-三角洲和湖泊沉积体系［25］（图1b）。花港组下伏平湖组发育优质烃源岩，自身砂体分布广泛且储层物性好，上覆区域性盖层，圈闭类型以滚动背斜和挤压背斜圈闭为主，具备良好的油气地质条件［13］。

INPEFA（integrated prediction error filter analysis）曲线可以突显出原始测井曲线中肉眼不易观察到的变化趋势［26］。通过最大熵频谱分析（MESA-GR）处理测井曲线得到估计值，然后用测井曲线实际值减去估计值得到预测误差滤波分析曲线（PEFA， prediction error filter analysis）［27］，再对PEFA曲线进行积分处理得到趋势更加明显的INPEFA曲线。趋势的变化由多个拐点隔开，INPEFA值向上增加表现为正趋势，向上减小则对应负趋势。INPEFA曲线的拐点和曲线变化趋势反映了不同级次的层序界面和沉积旋回变化；正趋势往往代表气候开始湿润，水位升高；负趋势则代表气候开始干旱，水位降低。

Daubechis（DB）小波变换是离散小波变换（DWT， discrete wavelet transform）的一种，本文采用Matlab中的db6小波变换分析该井伽马曲线，通过DB小波分析对GR进行重构，得到12条近似信号曲线（图2a1—a12）和12条细节信号曲线（图2d1—d12），其信号可定义为GR = a n + d1 + d2 + d3 + d4 + d5 + d6 + d7 + d8 + d9 + d10 + d11 + d12［28］。离散小波公式如下：

式中：Dn （t）为细节信号；An （t）为近似信号；n为分解程度。

近似信号曲线（图2a1—a12）是一种形态接近GR曲线的重构信号，其过滤掉了GR曲线上多余的细节，最大程度地保留了GR曲线的原始信息并将GR值变化的趋势完美地展现出来，选择合适阶次的近似信号曲线对识别GR曲线中所包含的最大湖泛面（MFS）信息有很好的辅助作用（近似信号曲线的最大值一般对应洪泛面）。

细节信号曲线（图2d1—d12）是一种突出GR周期变化细节的重构信号，把GR曲线中包含的界面信息以异常震动的形式凸显出来，洪泛面和层序界面都是对应于细节信号的异常震动。值得注意的是，单一阶次DB小波变换结果可能会造成界面位置发生轻微错位，从而不能精确地锁定界面位置［29］，综合运用多阶次的细节信号曲线可更加准确地确定层序界面位置。

Morlet 小波变换是连续小波变换（CWT， continuous wavelet transform）的一种，通过Morlet小波变换分析，将一维测井曲线（如GR）转换成二维频谱和小波系数曲线。

连续小波公式如下［30］：

式中：α为拉伸或压缩的尺度参数（α>0），β为小波位置偏移参数，f（t）是输入信号；Ψ为小波函数。α值越大，分辨率越低，用于划分长期旋回；反之，α值越小，分辨率越高，用于划分短期旋回或长期旋回。

离散小波（DWT）对层序界面和趋势变化响应灵敏，选择不同级次的小波重构曲线可以识别出潜在的不同级次的层序界面［31］；连续小波（CWT）将一维测井曲线转换为二维频谱，其对测井曲线的旋回性识别度较高且稳定性更强，通过改变尺度参数α可以得到能识别出不同级次旋回的二维频谱。

本文依据旋回地层学原理对花港组进行天文轨道周期的识别，对研究区的三级层序界面进行年龄厘定，将层序地层学与旋回地层学有机结合，避免了层序地层学的“等时而不定时”和旋回地层学的“定时而不等时”的问题。研究区构造活动稳定，海（湖）平面升降产生的沉积物可记录天文轨道周期［32］。GR曲线记录了K，Th和U放射性元素含量的变化，准确反映了沉积物中泥质含量的变化。泥质含量与湖/海平面变化和沉积物供给有关，因此GR曲线可作为古气候良好替代指标用于建立天文时间标尺［33］。

本文基于Acycle 2.1［34］软件，输入GR数据（软件自动对数据中空值和无效值进行插值），对其进行去趋势化处理；在此基础上，利用多窗口频谱分析方法对去趋势化GR数据序列进行频谱分析，其波峰代表着不同的天文轨道周期，重点选择95 %置信度的数据结果，对90 %置信度以下的数据结果则有选择地使用；对数据序列进行滑动窗口频谱分析及高斯带通滤波处理，利用滤波获取的405 kyr长偏心率周期调谐短偏心率周期和岁差周期［35-36］，并最终建立花港组“浮动天文时间标尺”。

本文以地震接触关系、INPEFA曲线及小波分析结果识别花港组顶底不整合界面（二级层序界面），并利用INPEFA曲线、小波重构曲线及频谱划分三级及四级层序界面，具体研究结果如下。

西湖凹陷渐新统顶底发育两个区域不整合界面，分别为平湖组与花港组分界面（底界面/SSB1）和花港组与龙井组分界面（顶界面/SSB2）；SSB1和SSB2均对应二级层序界面。

SSB1形成于始新世末期玉泉构造运动（33.9 Ma），该地质事件标志着盆地由裂陷阶段向拗陷阶段转变［22，37］，对应于地震界面T30（图4）。T30界面在地震上表现连续性强的反射轴，与上覆地层呈上超接触关系，与下伏地层呈削截接触关系（图3）。花港组以单束粉属-双束粉属-栎粉属-桤木粉属-高腾粉属组合为主，平湖组以桤木粉属-高腾粉属组合为主［38］。从平湖组到花港组，元素含量Ga/Rb比值和Sr/Cu比值均明显升高，表明花港组沉积时期较平湖组沉积时期寒冷［39］。

该界面处岩性由泥岩向厚层砂岩转变，GR曲线则在界面处发生从高幅线型向低幅钟形或锯齿状的变化。在INPEFA曲线上对应于次要负拐点，表明该界面处沉积环境发生转变，岩性由灰色泥岩向细砂岩及粉砂岩互层过渡。该界面在小波变换细节信号曲线d6（-5～5）和d7（3～16）上为低振幅振动，在d8细节信号曲线上为强振幅振动（-19～19），在小波变换近似信号曲线a7上为低值区（图4）。

SSB2形成于渐新世末期的花港运动（23.03 Ma）［22］，对应于地震界面T24，该界面在地震剖面上与上覆地层呈上超接触关系，与下伏地层呈削截接触关系（图3），该界面对应的地震反射轴连续性好且全区易于追踪（图3）。其在INPEFA曲线上位于次要正拐点（图4），对应由粉砂岩与灰色泥岩互层向灰色泥岩过渡，表明沉积环境发生转变。在细节信号曲线（d6，d7，d8）上呈高频弱振幅振动（-15～15），在近似信号曲线a7上对应相对低值区（82）（图4，图5）。

根据Catuneanu等［40］所倡导的层序地层学划分体系，三级层序通常包含一个完整的水进-水退旋回，三级层序界面通常对应不整合面或与之对应的整合界面；四级层序边界则主要为洪泛面或次级水进-水退旋回转换面。

以花港组顶底界面的深度为起止点确定窗口系数，对B3井GR曲线进行积分处理得到INPEFA曲线，通过识别其趋势和拐点分析花港组内部的不同级次层序界面和旋回结构，自下而上识别出了4个三级层序界面（SB1—SB5）。如图4所示，B3井INPEFA曲线上可识别4个明显的正拐点，对应深度分别为4 194.2，3 870.0， 3 662.3和3 468.8 m。这些正拐点对应182～400 m厚度的旋回变化转折，代表湖平面下降阶段结束转为上升阶段，岩性由泥岩过渡为砂岩；这些拐点深度可解释为三级层序界面。值得注意的是，SQ3的底界面（SB2）处发育一厚层砂岩（40 m），对应强烈河道侵蚀面，造成INPEFA曲线在该界面呈现假性负拐点的特征（图4）。根据INPEFA曲线的相对大尺度和明显拐点和趋势划分旋回期次，能提高层序地层研究的准确性和可行性。

以花港组的顶底深度作为处理窗口的起止深度，对B3井GR进行尺度因子α = 512的连续小波变换（图4）和离散小波变换并选取细节信号曲线d6，d7，d8和近似信号曲线a7展开层序地层界面信号特征研究。其中，细节信号曲线d6整体呈现波峰波谷频繁交替但异常震荡显著，而d7和d8波峰波谷交替频率减慢并且曲线整体平缓，异常震荡较为明显；近似信号曲线a7最低值处发育大套的砂岩，最大值处发育砂泥互层，泥质含量高（图4）。细节信号曲线异常震荡处对应近似信号的最大值或最小值。频谱分析发现能量团尺度自下而上逐渐减小再增大，其颜色由暖色调逐渐变为冷色调再变为暖色调，代表一个大尺度的沉积旋回变化（图4）。这些能量团小波信号曲线的震荡位置（图4）与INPEFA曲线的拐点有良好的对应关系。因此，小波变换曲线周期性的震荡特征可为三级层序界面的定量识别提供佐证。

以SQ2为例，SQ2底界面（SB1）处细节信号曲线d6（-16～25）和d7（-43～48）均有明显的异常高频强振幅，d8振幅较弱，曲线表现为弱振幅（-11～16）向平滑过渡（图4）。SQ2顶界面（SB2）处d6有明显的异常振幅（-31～28），表现为高频强振幅（-31～28）向高频弱振幅（-8～13）转变，d7为强振幅振动（-43～48），d8振幅较弱，曲线从平滑过渡为低强度振幅（-11~18）（图4）。此外，SQ2底部发育一套薄的砂岩及冲刷面（图6），自下而上在频谱图上表现出大尺度能量团—小尺度能量团—大尺度能量团的变化趋势（图4），代表了水体由浅变深再变浅及水动力减弱再增强的基准面变化；对应于沉积物粒度由粗—细—粗的变化趋势，与INPEFA曲线的正趋势及负趋势变化具有高度的一致性。

SB3是花港组内部一个重要的三级层序界面，对应为T25地震界面。T25界面之上发育一系列地震上超现象（图3），界面对应于INPEFA曲线上正拐点位置（图4），岩性由灰色泥岩向灰色细砂岩转变（图6a）。该界面在小波细节信号d6曲线上表现为从弱振幅振动（-15～22）过渡为强振幅振动（-23～31），在d7上为强振幅振动（-46～34），在d8上强振幅振动（-2～3）（图4）。类似地，综合利用INPEFA曲线和小波变换进行层序地层研究，可将花港组自下而上划分为5套三级层序（SQ1—SQ5；图4，图5）。

在三级层序格架内，利用次一级地层叠置旋回信息可进一步识别四级层序（准层序组）。其中，利用INPEFA曲线负拐点及变化趋势，可很好地识别三级旋回内部的四级层序界面。

小波变换通过测井曲线的伸缩、平移把信号分解为不同尺度的成分，进而将测井曲线分解为不同周期的旋回［28］。近似信号曲线（图2a1—a12）是通过过滤测井曲线的细节特征而保留其变化趋势，最高值处一般对应洪泛面，最低值处对应层序界面［29］。细节信号曲线（图2d1—d12）则是屏蔽测井曲线的趋势信号，放大测井曲线的细节信息，异常震荡处对应层序界面。在单一层系内沉积岩岩性和物性变化具有方向性和连续性［41］，对应的频谱是通过能量团尺度的大小和颜色的强弱变换趋势显示测井曲线中包含的旋回信息以及反映地层的叠加样式。为了便于运用频谱进行体系域和地层叠加样式的研究，首先要构建体系域和地层叠加样式的频谱地质模型。

湖侵体系域在频谱上表现为自下而上由大尺度能量团向小尺度能量团变化且能量团凹面向上，其能量团的颜色逐渐由暖色调变为冷色调（图7a）；湖退体系域在频谱上表现为自下而上由小尺度能量团向大尺度能量团变化且能量团凹面向下，其能量团的颜色逐渐由冷色调变为暖色调（图7b）。

在体系域划分的基础上，在其格架内开展准层序组研究，准层序组对应四级层序。退积准层序组在频谱上表现为自下而上由大尺度能量团向小尺度能量团变化，能量团的颜色逐渐由暖色调变为冷色调（图8a）；进积准层序组在频谱上表现为之下而上由小尺度能量团向大尺度能量团变化，能量团的颜色逐渐由冷色调变为暖色调（图8b）；加积准层序组在频谱上表现为自下而上能量团尺度大小变化极小，能量团颜色基本没有变化（图8c）。

以SQ2为例，其岩性垂向序列向上砂岩含量先减少再增多，内部可识别明显的正旋回和负旋回转换界面，即对应INPEFA曲线的负拐点，因此SQ2可分为两套四级层序（H10和H9）。频谱（α = 512）图像分析显示，自下而上SQ2能量团尺度呈现出由大逐渐变小再变大的趋势，利用所建立的频率模型可分别解释为退积准层序组（H10）和进积准层序组（H9）。在小波变换曲线上，H9与H10界面（mfs2）处细节信号曲线d6显示高幅度振幅（-18～19）过渡为低幅度振幅（4～6）（图4）；近似信号曲线a7峰值指示洪泛面。通过以上分析，认为频谱成像、小波变换与INPEFA曲线的变化趋势在识别四级层序方面具有良好的一致性。

如图4所示，在α = 512的频谱中存在5个明显的旋回（图4），并与INPEFA曲线的变化趋势吻合。异常震荡处多为砂岩相变泥岩相或泥质含量快速变化或存在河道冲刷面等地层突变，反映了沉积环境的变化；频谱通过能量团变化的趋势反映沉积环境变化的过程。类似地，在SQ1，SQ2和SQ4中各识别了两个四级旋回及3个洪泛面（分别是mfs1，mfs2，mfs4），在SQ3和SQ5中各识别了3个四级旋回以及两个洪泛面（分别是mfs3，mfs5）。据此把研究区花港组进一步划分为12套四级层序（准层序组），自上而下分别命名为H1—H12（图4）。此外，研究区A1，A2，B1及B2井的层序划分结果与B3井的划分结果也基本一致（图5）。

本次将B3井的GR数据序列分为花上段和花下段，分别进行多窗口频谱分析，选取顶底界面为时间锚点做多窗口频谱分析，其频率的倒数即为旋回厚度。通过多窗口频谱分析，花港组下段GR数据序列识别了63.93，15.11，13.51，13.20，4.75，4.84，4.90，4.93和5.41 m旋回厚度，选定63.93 m旋回厚度为405 kyr长偏心率周期，则15.11，13.51和13.20 m分别代表96，86和84 kyr短偏心率周期；4.84，4.93和5.41 m旋回厚度分别代表31，32和35 kyr斜率周期（图9a）。

花港组上段GR数据序列共识别了85.84，20.04，19.24，17.81，8.56，6.94，6.87和6.76 m旋回厚度。选定85.84 m旋回厚度为405 kyr长偏心率周期，则20.04，19.24和17.81 m分别代表95，90和84 kyr短偏心率周期；8.56，6.94，6.87和6.76 m旋回厚度分别代表40，33，32和30 kyr斜率周期（图9b）。

基于GR频谱分析，识别了长偏心率、短偏心率、斜率周期及其旋回厚度。在此基础上，通过演化谱分析及滤波处理，识别出405 kyr长偏心率周期，再以405 kyr长偏心率周期调谐短偏心率及斜率周期。花港组下段GR数据序列通过滤波处理输出了10个405 kyr长偏心率周期（图10a），其高斯滤波频率（0.020 455 ± 0.001 504） cycle/m，31个94 kyr短偏心率周期，其高斯滤波频率为（0.073 999 ± 0.001 805）cycle/m，演化谱滑动窗口为132.95 m。花港组上段GR数据序列通过滤波处理输出了11个405 kyr长偏心率周期（图10b），其高斯滤波频率为（0.014 937 ± 0.032 87） cycle/m，35个100 kyr短偏心率周期，其高斯滤波频率为（0.049 888 ± 0.006 572） cycle/m，演化谱滑动窗口为133.875 m。

频谱分析和演化谱分析表明，西湖凹陷花港组B3井沉积地层记录了长期稳定的405 kyr长偏心率周期和相对稳定的100 kyr短偏心率周期。以稳定的405 kyr长偏心率周期建立西湖凹陷B3井花港组上段和花港组下段“浮动天文标尺”，花港组三级层序持续时间跨度介于1～3 Myr。其中，SQ1沉积时间为1.2 Myr，H12持续时间约为0.405 Myr，H11持续时间约0.81 Myr。SQ2沉积时间跨度约2.0 Myr，H10和H9发育持续时间均为1.0 Myr左右；SQ3时间跨度为2.8 Myr，其中H8和H7共持续时长为2.0 Myr，H6持续时长约0.8 Myr。SQ4沉积时间为1.62 Myr，H5和H4持续时长分别为0.4 Myr和1.2 Myr；SQ5层序持续时间跨度长达 2.43 Myr，其中以H2和H1占据主导地位（约1.62 Myr）。

一级层序、二级层序及三级层序以不同规模不整合面或沉积间断面为界面，四级层序与体系域层序单元则以洪泛面为界面［42-43］。一级层序代表盆地从形成到衰亡的整个沉积序列，通常以板块尺度的构造不整合面为界面；二级层序一般包括完整的区域性水进水退沉积旋回，以区域性不整合及冲刷面为边界；始新世时期东海陆架盆地由拉张应力场转变为挤压应力场［22］，东部钓鱼岛褶皱带垂直隆起［44-45］。花港组顶底二级层序界面SSB1和SSB2分别对应玉泉运动（33.9 Ma）和花港运动（23.03 Ma）形成的不整合面，与此对应的是Oi-1变冷事件和Mi-1变冷事件及其海平面大幅度下降［24，46-47］。

三级层序界面以盆地或凹陷尺度的微角度不整合及与之对应整合面为界面，界面形成主要与海（湖）平面升降、气候变化和冰川消融等因素有关［48］。三级层序内部还可包含多个次一级水进水退的沉积旋回转换界面，可作为四级层序界面［42］。利用数学手段的测井资料分析或元素地球化学及同位素测试方法可划分更高级次的层序单元［7，49-50］。GR曲线记录了沉积物中U，Th和K等放射性元素含量，GR高值区泥质含量高，为湖平面上升至最高处；反之，则是湖平面下降到最低处。因此，基于GR数据获取的INPEFA曲线趋势变化可间接反推湖/海平面的变化［51］。西湖凹陷中南部花港组基于GR的INPEFA曲线拐点和趋势变化反馈了5次明显的水进-水退变化，代表了5个湖平面升降旋回（图4），与渐新世海侵-海退旋回具有高度的一致性［52］；这与本论文所建立的由体系域频谱模型所反馈的地层叠置样式亦具有良好的耦合性。综合研究表明，三级层序SQ1，SQ2和SQ4主要发育对称半旋回层序结构，而SQ3和SQ5发育非对称型半旋回层序结构（图4）。西湖凹陷花港组时期构造沉降缓慢且物源供给稳定［53］，因此湖平面升降周期及对称性决定了三级层序内部地层结构，并主导形成了多类型的四级层序旋回叠置样式（图4，图5）。

1） 通过对西湖凹陷花港组GR曲线进行小波变换、最大熵谱以及地震接触关系分析，并结合体系域与地层叠加样式地质模型，可将花港组划分为1套（准）二级层序、5套三级层序和12套四级层序。

2） 利用最大熵谱分析对GR曲线进行数学分析识别了花港组沉积时期5次明显的湖平面升降旋回，与5套三级层序有良好的耦合关系；花港组内部的三级层序和四级层序（准层序组）叠置样式主要受湖平面升降周期及对称性调配，并且西湖凹陷花港组记录了稳定的405 kyr长偏心率周期和100 kyr短偏心率周期，且三级层序持续时间介于1.01～3.04 Myr，四级层序持续时长集中在0.4～1.0 Myr。

3） 频谱对单井体系域和地层叠加样式的识别效果较好，DB小波中d5，d6和d7细节信号曲线组合对层序界面响应明显，INPEFA曲线跟更能体现出曲线的变化趋势；两种方法在层序划分中可相互验证，提高层序界面和旋回信息识别准确度。