基于小波分析的层序地层划分及识别——以川东地区上二叠统龙潭组为例

余瑜，林良彪，蓝彬桓，洪薇，郭炎

(1.油气藏地质及开发工程国家重点实验室(成都理工大学)，四川 成都 610059；2.成都理工大学沉积 地质研究院，四川 成都 610059；3.国家电网遂宁供电公司，四川 遂宁 629000)

层序地层学最初由VAIL等人于20世纪70年代提出(VAIL P R, et al., 1977)，经过这数十年的发展，在得到广泛应用的同时也被地质学家赋予了新的地质意义。近年来，层序地层学在沉积演化、油气勘探领域被越来越多的使用(陈洪德等，2011；林良彪等，2007；于波，2016；尹青等，2015)，陈洪德等(2011)通过分析中上扬子克拉通华力西—印支期的层序剖面结构、体系域叠加样式等，揭示了中上扬子地区华力西—印支期的盆地层序充填特征；林良彪等(2007)通过野外剖面、钻井和地震资料，在识别层序界面的基础上进行层序划分，将四川盆地上三叠统划分为3个构造层序。技术方法的更新进步使得层序地层学的研究向着高频层序的方向发展，受冰川型海平面变化和气候影响的高频程序(四级—六级层序)引起越来越多的地质学家重视(尹青等，2015；张运波等，2011；赵宗举等，2010；徐伟等，2012)。高频层序的时间跨度不到1Myr，甚至可能仅数个10Ka。传统的层序地层学研究方法与技术手段在高频层序的研究中效果甚微，不少的学者利用频谱分析技术对地层的测井曲线进行高频旋回层序研究，更新和完善了高频层序识别与划分相关的方法与技术(尹青等，2015；张运波等，2011；赵宗举等，2010；徐伟等，2012；陈钢花等，2007；王起琮，2009)。

随着四川盆地下寒武统筇竹寺组、上奥陶统五峰组—下志留统龙马溪组的页岩气勘探相继获得重大突破(管全中等，2015；张海全等，2011)，同样作为四川盆地古生代三套泥质烃源岩之一的龙潭组(吴家坪组)相应的研究却相对缺乏。前人对于龙潭组/吴家坪组的研究主要集中在层序格架、沉积环境及烃源岩的评价上(徐强等，2004；林良彪等，2009；李凤杰等，2009；张廷山等，2017；腾格尔等，2010；赵宗举等，2012)，且对于地层层序的划分存在着一定的争议(徐强等，2004；林良彪等，2009；李凤杰等，2009)。基于此，结合前人研究成果，笔者以川东地区上二叠统龙潭组的钻测井资料为主要依据，利用频谱分析和Fischer图解对龙潭组进行高频层序地层识别、划分的研究，分析其可容纳空间变化的规律及对米兰科维奇沉积旋回的记录。

1 地质背景

研究区位于四川盆地东部，构造单元上属于上扬子地区(四川油气区石油地质编写组，1989)。根据现今地表地形，四川盆地划分为川西低陡构造区、川北低平构造区等6个构造单元，研究区位于川东高陡构造区(图1)。中晚二叠世之交的东吴运动致使扬子板块西南缘的康滇古陆隆升成陆，并成为四川盆地的主要物源区之一，四川盆地整体上呈现出西南地势相对较高、向北逐渐降低的特点，在盆地内西部和南部发育了一套海陆交互相地层，岩性以暗色泥岩、碳质泥岩为主，夹有煤线，即龙潭组(P3l)；向盆地内北部、东部相变为海相地层，即吴家坪组(P3w)(四川油气区石油地质编写组，1989)，笔者为简便，统称为龙潭组。

1.井位；2.研究区位置；3.盆地边界；4.构造分区图1 研究区地质构造位置图Fig.1 Regional tectonics location of study area

2 数据的选择与处理

对于高频层序的研究，不同学者采用的原始数据有不小的差别，夏国清等(2010)根据野外露头的地层颜色、岩性等特征进行了米级沉积旋回的划分；赵宗举等(2010)利用GR及LLD测井曲线，结合钻井岩心、地震反射特征及碳氧同位素曲线等资料，研究了塔中—巴楚地区上奥陶统碳酸盐岩发育的高频层序特征；尹青等(2015)、张运波等(2011)、王起琮(2009)、彭博等(2013)均利用GR测井资料进行了高频(米兰科维奇旋回)层序的研究，并得到了较为理想的结果。研究区内龙潭组野外露头出露情况较差，为了追求更加精确的结果，本次研究采用川东地区龙潭组钻井的测井资料进行分析。研究区内龙潭组整体上是一套以泥页岩为主的夹碎屑岩或碳酸盐岩的地层，相对于砂岩和碳酸盐岩来说，泥岩颗粒更为细小，比表面积相对更大，更易吸附放射性元素。如K、Th、U等元素，从而具有相对更高的伽马值，因此，在研究中选择相对于泥质含量较为敏感的自然伽马(GR)曲线进行分析研究。

选取研究区内未发育断层及不发育火山岩的井位的GR曲线，应用Matlab软件的小波分析工具箱(Wavelet Toolbox)对钻井自然伽马数据进行预处理。首先，对井位的GR数据进行0.125 m的等间距数据采集，即在厚度为1m的地层采集8个数据；其次，因为测井曲线可能记录了与地质作用不相关的数据，为消除人为因素及环境因素对自然伽马测井曲线的影响，在小波工具箱中对原始的GR数据进行降噪(D-noise)处理；之后，根据前人的研究方法(王起琮，2009；彭博等，2013)，将降噪处理后的GR数据放入Matlab软件的小波工具箱进行dbN小波变换。dbN小波能够有效地分解含多种频率成分的复杂信号(王起琮，2009)，最终得到了多层的小波分解曲线；最后，再对降噪处理后的GR数据进行连续小波变换(Continuous Wavelet)，选取morl小波对数据进行分析处理，得到了相应的二维能量谱图(图2、图3)。

3 龙潭组旋回划分

3.1 地层旋回划分

以普光5井为例，得到小波分析的图解和层序地层划分(图2)。结合二维能量谱图和具较强的旋回性质的dbN曲线进行龙潭组内部的三级层序划分。

普光5井龙潭组厚度为197 m，岩性以灰岩和泥岩为主，morl小波分析得到的二维能量谱图显示，普光5井龙潭组底部具有相对较强的能量，这与底部泥岩相对发育致使GR测井曲线数据具正高值有关。据能量谱图和d9、d8小波分解曲线，普光5井龙潭组可划分为3个三级层序，分别是SQ1、SQ2、SQ3。其中，SQ1深度分布范围为5 696.5～5 747.5 m，SQ2为5 659.5～5 696.5 m，SQ3为5 550.5～5 659.5 m，泥岩、灰质泥岩主要分布于SQ1、SQ2。据前人研究表明，受东吴运动和海平面下降的影响，龙潭组底界面为造山侵蚀不整合界面，属于I型层序界面(林良彪等，2009)，龙潭期早期四川盆地东部广泛发育低位体系域(LST)。因此，将SQ1划分为低位体系域、海侵体系域(TST)及高位体系域(HST)，厚度上TST与HST大致相等，LST厚度较小。而SQ2、SQ3底界面均属于II型层序界面，均缺乏陆棚边缘体系域沉积，只识别出了海侵体系域与高位体系域，且TST与HST厚度近似，发育TST、HST近完全对称型的层序。前人研究表明，这类型层序是形成于较高的可容纳空间及较低的沉积物供给的沉积环境下，水动力条件相对较弱(郑荣才等，2010)。

据前人研究显示，川东北的吴家坪组与川南、川西的龙潭组呈渐变过渡关系，具低能环境沉积特征，发育碳酸盐缓坡沉积。吴家坪组后期，由于持续的海侵导致的快速沉积作用及碳酸盐岩生产率的增加，碳酸盐缓坡向台地转变，最终形成了晚二叠世长兴期的碳酸盐台地沉积(张帆等，1993)。沉积环境的变化控制着龙潭组层序厚度及岩性的变化，发育泥页岩的SQ1、SQ2厚度明显小于以灰岩沉积为主的SQ3，SQ1和SQ2主要发育深水碳酸盐缓坡沉积，主要为低能环境的泥页岩、薄层灰岩等沉积。SQ3发育浅水碳酸盐缓坡沉积，沉积了一套以灰岩为主夹少量泥质灰岩的岩石，泥岩几乎不发育，沉积环境整体上水动力条件依然较弱，灰岩以泥-微晶灰岩为主。

应用同样的方法对自深1井进行小波分析与层序划分(图3)，得到SQ1、SQ2、SQ3三个层序。自深1井龙潭组岩性以泥岩、页岩夹薄层状砂岩为主，在层序划分上，与普光5井不同的是，自深1井的SQ1层序未发育LST,这也是造成自深1井与普光5井发育米级旋回数量存在差异的原因。

图2 川东地区上二叠统龙潭组小波变换与层序划分图(普光5井)Fig.2 Sequence division and wavelet transform of the Longtan Formation in eastern Sichuan Basin (Puguang 5 Well)

3.2 Fischer图解

Fischer图解是一个以时间为横坐标，以米级旋回经过线性沉降校正的累计厚度偏移为纵坐标的一个变化图解(苏德辰等，1995)。曲线代表着沉积时可容纳空间的变化，因此Fischer图解被广泛应用于以碳酸盐为主的海相地层的旋回地层划分及海平面升降变化的研究中(彭博等，2013；苏德辰等，1995)。在这需要指出的是米级旋回是指地层记录的在野外露头上可以直接识别出的、几十厘米到几米厚的地层堆积岩石单元(ANDERSON E J.，et al.，1990；梅冥相等，2000)。米级旋回又称之为高频旋回(MITCHUM R M.et al.,1991)，在经典层序地层学的划分中通常指的是四级层序，即准层序(VAIL P R.et al.,1991)，是在米兰科维奇天文周期旋回的影响下形成的岩石地层。在构造相对稳定地区，米兰科维奇旋回常受控于气候的影响，而气候的变化与地球轨道参数有着紧密的联系。

图3 川东地区上二叠统龙潭组小波变换与层序划分图(自深1井)Fig.3 Sequence division and wavelet transform of the Longtan Formation in eastern Sichuan Basin (Zhishen 1 Well)

在野外剖面上，识别米级旋回的方法主要根据露头的岩性、岩相、沉积构造及垂向沉积序列进行识别。在本次研究中，笔者采用对自然伽马测井曲线进行小波分解，提取显示相对高频的d4曲线进行旋回划分(图2)，划分依据是将d4分解曲线上相邻的两个波峰之间的部分划分为一个米级旋回。应用这个方法，将普光5井龙潭组划分为51个米级旋回，普光5井龙潭组厚度为197 m，因此得到每个米级旋回平均厚度为3.86 m；自深1井龙潭组划分为46个米级旋回，自深1井厚度为121.5 m，每个米级旋回萍姐厚度为2.64 m。

普光5井龙潭组的Fischer图解如图4，可以看出与前文三级层序的划分具有良好的契合程度：曲线由左向右可以划分为3个旋回，其中旋回1至旋回14划分对应层序SQ1，曲线先后经历了降低、升高、再降低的趋势，累计方差大部分情况下小于0，表明可容纳空间经历了降低—升高—降低；旋回15至旋回24对应层序SQ2，累计方差均小于0，曲线在旋回20升至最高，之后持续降低，曲线的上升和下降对称性较明显，表明海平面的升降变化较稳定；旋回25至旋回51对应为SQ3，发育了大套的灰岩及泥灰岩，可容纳空间经过短暂的下降之后快速上升，在旋回34升至最高，之后进入可容纳空间下降阶段。对比自深1井与普光5井的Fischer图解，二者间存在着较大的差异，普光5井的SQ1旋回1-旋回3可容纳空间存在着降低的现象，即SQ1发育有LST，而自深1井未见LST的发育，反映可容纳空间的曲线快速上升，然后经过震荡波动后再次下降；普光5井SQ2呈现对称性的上升—下降形态，自深1井SQ2从旋回13至旋回33呈震荡下降，岩性剖面上表现为泥页岩夹薄层砂岩；普光5井SQ3岩性以碳酸盐岩为主，厚度显著大于以泥页岩为主的自深1井SQ3。

图4 普光5井、自深1井龙潭组Fischer图解与层序划分图Fig.4 Fischer plots and sequence division of Puguang 5 Well and Zhisheng 1 Well

普光5井与自深1井的Fischer图解与三级层序具有良好的对应关系，对于可容纳空间变化的记录有着较好的显示，表明利用小波变换对自然伽马测井曲线进行三级层序及米级层序的划分是一种有效地方法。

4 频谱分析

结合PAST2.0软件Spectral analysis程序对普光5井及自深1井龙潭组自然伽马曲线进行频谱分析(夏国清等，2010；伊海生，2011)，得到普光5井及自深1井的频谱分析如图5所示。频谱分析结果显示具有多个能量峰值点，普光5井的频谱分析在置信度90%以上的分别对应的频率数值为0.002 93、0.009 76、0.017 56、0.026 37、0.033 20、0.052 73、0.092 773；自深1井的频谱分析分别对应的频率数值为0.012 295、0.015 369、0.025 615、0.032 787、0.048 156、0.056 352、0.146 27。频谱分析中能量峰值对应的频率的倒数为旋回平均厚度(平均旋回厚度=1/频率/采样密度)，表1列出了频率及对应的平均旋回厚度。已知普光5井龙潭组厚度为197 m，自深1井龙潭组厚度为121.5 m，用龙潭组厚度除以相应的旋回平均厚度可得相应的旋回个数。

井位序号频率平均旋回厚度(m)旋回个数(个)旋回周期(ka)沉积速率(m/ka)普光5井10.002 9342.662 12520.009 7612.807 3815偏心率长周期0.032 030.017 567.118 452840.026 374.740 244250.033 203.765 0652偏心率短周期0.037 760.052 732.370 578370.092 771.347 38146斜率周期0.035 5自深1井10.012 2910.166 731220.015 368.133 2615偏心率长周期0.020 330.025 624.879 952540.032 793.812 493250.048 162.595 7347偏心率短周期0.026 060.056 352.218 205570.146 270.854 58142斜率周期0.022 5

BERGER A et al.(1994)认为地球在不同时期轨道周期存在着周期性的变化，轨道周期据时间跨度的不同包括了偏心率周期、斜率周期、岁差周期，其中偏心率周期分为长周期和短周期，长周期为413 ka，短周期为123 ka、95 ka(平均为100 ka)，斜率周期为50.6 ka、38.0 ka，岁差周期为21.3 ka、17.6 ka。普光5井中，如果①点对应的偏心率长周期(时间跨度413 ka)，即f=0.002 93时，则龙潭组时限为5(个)×413 ka=2.65 Myr。根据最新的国际年代地层表(ICS，2016)，晚二叠世包括了吴家坪阶[Wuchiapingian，(259.8±0.4)～(254.14±0.07 Ma)、长兴阶[Changhsingian，(254.14±0.07)～(252.17±0.06)Ma]。其中，吴家坪阶时限大致为5.19～6.13 Myr，显然f=0.002 93的峰值不能代表普光5井的偏心率长周期。若②点即f=0.009 76对应的是偏心率长周期，则龙潭组时间跨度为15×413 ka=6.195 Myr，接近于当今国际公认的吴家坪介时限。将f=0.033 20峰值点确认为偏心率短周期，即发育了52个偏心率短周期旋回。而⑦点f=0.092 773峰值点对应的是斜率周期。自深1井分别发育15个偏心率长周期、47个偏心率短周期及142个斜率周期。

通过频谱分析，识别出了普光5井及自深1井龙潭组的偏心率长、短周期及斜率周期，表明川东地区龙潭组对于米兰科维奇旋回具有较为完整的保存。值得注意的是，普光5井和自深1井龙潭组自然伽马的频谱分析均发育了15个偏心率长周期，表明即使整体上沉积环境相异，但是对于米兰科维奇旋回的偏心率长周期的记录依然较为完整，而短周期及斜率周期旋回的不同则是受沉积环境的控制，普光5井龙潭组沉积速率为0.032 0～0.037 7 m/ka，自深1井沉积速率为0.020 3～0.02 60 m/ka，普光5井SQ3沉积环境由碳酸盐缓坡向碳酸盐台地转变，碳酸盐生产率大大的提高。如果依据BERGER等的研究成果，偏心率长周期为413 ka，则频谱分析得到的数据计算吴家坪介时间跨度为6.195 Myr，虽然接近于6.13 Myr，但显然是不合理的，同时这是建立在地层完整、没有沉积间断出现的情形下。LASKAR J et al.,(2004)研究了北纬38°夏至日30～40 Ma的地球轨道周期的理论值，得到偏心率周期为405 ka、125 ka、99 ka、95 ka，表明在地质历史演化过程中，地球的轨道周期是有小幅度波动的，并非固定不变的，上扬子地区晚二叠世偏心率长周期应略小于413 ka。

遗憾的是，普光5井和自深1井的频谱分析中，在f>0.1之后(除自深1井的⑦点外)，频谱分析数据均在90%置信度曲线之下，无法进一步的识别出龙潭组的岁差周期。

5 结论

(1)应用Matlab小波分析工具对川东地区龙潭组钻井的自然伽马曲线进行分析，得到相应的小波分解曲线及二维能量谱图，可将龙潭组划分为3个三级层序，其中SQ1发育低位体系域、海侵体系域及高位体系域，而SQ2、SQ3可分别识别出海侵体系域与高位体系域。

(2)基于小波分解曲线得到的Fischer图解与三级层序具有良好的对应关系，对于可容纳空间的变化具有良好的显示效果。

(3)通过对龙潭组自然伽马曲线进行频谱分析，显示龙潭组明显受到了米兰科维奇天文周期的控制，对于米兰科维奇旋回具有较为完整的记录。川东龙潭组可识别出偏心率周期、斜率周期。其中，龙潭组发育15个偏心率长周期，而受沉积环境不同的影响，发育的偏心率短周期及斜率周期有所不同。