东沙海域内孤立波浅化过程及含圈闭涡核内孤立波的反射地震初探

宋海斌， 邝芸艳， 杨胜雄， 关永贤， 龚屹， 范文豪， 张锟

1 海洋地质国家重点实验室， 同济大学， 上海 200092 2 海洋与地球科学学院， 同济大学， 上海 200092 3 南方海洋科学与工程广东省实验室(广州)， 广州 511458 4 中国地质调查局广州海洋地质调查局， 广州 510760

0 引言

内孤立波向陆坡陆架传播是一个复杂的动力学过程，该过程称为内孤立波的浅化过程(Cai et al., 2012; Guo and Chen, 2014; Lamb, 2014; Alford et al., 2015).当内孤立波传播到浅水时，由于非线性和频散效应之间平衡被打破，波形会发生变化.在浅化过程中，内孤立波可能会从下凹型转变为上抬型(Liu et al., 1998)，并且可能由于不稳定而发生破碎(Aghsaee et al., 2010)，或在内部形成涡核(Lamb, 2002).这些动力过程已经通过数值模拟(Holloway et al., 1997; Liu et al., 1998; Zhao et al., 2003; Grimshaw et al., 2010)和物理模拟(Boegman et al., 2005; Cheng and Hsu, 2010)进行研究，并通过锚系、高频声学和遥感等手段进行了现场观测(Orr and Mignerey, 2003; Zhao et al., 2003; Lynch et al., 2004; Bourgault et al., 2007; Shroyer et al., 2009; Fu et al., 2012；Chang et al., 2021; Chang, 2021).

地震海洋学(Holbrook et al., 2003；宋海斌，2012)近期一个主要进展在于内孤立波的研究(Song et al., 2021a).如拜阳等(2015a)利用地震海洋学方法研究了南海东北部东沙海域内孤立波的结构特征，并比较了内孤立波水平波数谱与GM模型谱，发现其能量在低波数段比GM谱大两个数量级左右.Tang等(2016)使用反射地震数据与局部水文数据通过马尔科夫链—蒙特卡洛反演方法获得了内孤立波的温度和盐度分布.Bai等(2017)利用反射地震数据系统分析了南海东北部浅化内孤立波的极性转换.Tang等(2018)用高分辨率地震成像和海水结构特征联合反演研究了南海北部陆架次级内孤立波.Geng等(2019)用南海北部众多反射地震数据识别出来的内孤立波，研究了内孤立波振幅与水深的关系，结果表明内孤立波最大振幅随水深呈对数减小.范文豪等(2020)在南海东北部海盆首次利用地震海洋学方法发现了海盆中的内孤立波.Gong等(2021a)研究了浅化内孤立波的振幅垂向结构.Gong等(2021b)研究了极性转换内孤立波对湍流混合增强的影响.Song等(2021b)提出了可以研究细结构变化的共炮检距剖面叠前偏移方法，并给出了一个浅化内孤立波波形变化的实例.Fan等(2021,2022)研究了中美洲太平洋海域的第二模态内孤立波特征.

海洋中的湍流混合过程是间歇性的和瞬时的.因此，需要大量观测才能准确评估湍流混合的空间分布(Whalen et al., 2012; Waterhouse et al., 2014; Kunze, 2017).地震海洋学具有观测范围宽、空间分辨率高(Ruddick et al., 2009；Song et al., 2021a)的优点，适合观测湍流混合的空间分布.Sheen等(2009)使用反射地震数据给出了南大西洋海洋锋面的混合参数剖面.Holbrook等(2013)详细介绍了利用反射地震数据估算湍流混合的理论基础.随后，越来越多的学者利用反射地震方法研究了不同海洋区域湍流混合，以及不同海洋现象引起的湍流混合的空间分布特征(如拜阳等，2015b；Fortin et al., 2016; Dickinson et al., 2017; Mojica et al., 2018).南海东沙海域上陆坡区域内孤立波浅化过程复杂，湍流混合时空分布很不均匀.利用高分辨率的地震海洋学方法可以给出内孤立波波形结构和一系列定量参数(如最大振幅、传播速度、半高宽、振幅垂向结构等)，估算混合参数(耗散率、扩散率)的空间分布.通过叠合波形结构与混合参数空间分布，可以较好地分析内孤立波浅化过程与能量耗散、混合增强的关系(Gong et al., 2021b)，弥补常规物理海洋观测的不足.

本文利用南海东沙海域的一条地震剖面，通过地震数据处理、利用叠前数据估算内孤立波传播速度、利用地震数据估算混合参数等方法，详细分析浅化内孤立波/波包的特征，初步探究浅化过程中形成的含圈闭涡核内孤立波与第二模态内孤立波的相关参数以及对湍流混合的影响，以深化对东沙海域复杂的内孤立波浅化过程的认识.

1 数据与方法

1.1 地震数据采集参数及处理方法

本文主要使用二维多道反射地震数据观测并分析内孤立波浅化的特征，该数据是由广州海洋地质调查局的“探宝号”地球物理综合调查船于2009年夏季在南海东沙海域采集.本文使用了二维测网中的12号多道地震测线(图1).采集过程中使用的拖缆全长6 km，包含480道水听器组合，道间距为12.5 m，采样间隔为2 ms.气枪震源容量是5080 in3(1 in=2.54 cm)，震源主频35 Hz.炮间距为25 m，最小偏移距为250 m，放炮间隔时间约为10 s.

图1 测线位置图

常规的地震数据处理流程主要包括定义观测系统、噪声和直达波压制、速度分析、动校正、水平叠加和偏移.首先对于一个原始的地震数据我们需要根据采集信息定义一个观测系统.第二步需要压制地震数据中的噪声和直达波.地震数据中通常会掺杂大量噪声，如随机噪声，海洋生物或他船干扰等等.与反射地震波相比，这些噪声的频率通常分布在低频区域或者高频区域.因此，我们采用一个带通滤波器去除8 Hz以下的低频噪声和100 Hz以上的高频噪声.此外，海水的反射界面通常较浅，这使得反射波容易受到直达波的干扰.在常规地震数据处理中通常直接切除直达波部分，但是在海水地震数据处理中这样操作会失去相当大一部分有效数据，尤其是浅部的海水层反射.为了解决这一问题我们采用压制直达波的方式来减弱直达波的影响.具体做法是先用一个中值滤波器对数据进行滤波，由于直达波和反射波振幅差异较大，中值滤波器得到的主要是直达波成分.然后在用原始数据减去中值滤波的结果便可以压制大部分直达波成分，这样就达到了恢复浅部反射波数据的目的.第三个步骤是速度分析，这一步的目的是获得海水中的声速分布，为下一步动校正做准备.速度分析完成后，利用得到的速度模型将共中心点(Common Mid-Point, CMP)道集数据中的反射同相轴拉平，这一步就是动校正(Normal move-out, NMO).动校正会使大偏移距的地震波形被拉长产生畸变，需要切除这些畸变的数据.传统的线性动校拉伸切除了浅部大量反射数据，这不利于浅部的观测.我们采用手动切除动校拉伸的方式，保留了部分浅部数据.

随后，将每个共中心点道集中的各个道叠加，再按照共中心点的位置排列就可以得到叠后地震剖面.这个叠后地震数据容易受到垂向噪声和残留直达波的影响.垂向噪声主要是震源鸣震引起的，带通滤波器不易去除.为了提高数据的信噪比，我们使用倾角滤波的方法去除数据中的垂向噪声和残留直达波.此外，当反射界面倾斜较大(不是水平界面)时叠加剖面上的同相轴与界面真实的位置有偏移，我们使用叠后时间偏移使反射同相轴归位到正确的成像位置.完成上述步骤后我们就可以得到一个高信噪比的叠后地震剖面.有关地震数据处理的原理和细节可以参考Ruddick等(2009).

1.2 利用叠前地震数据与气候态温盐数据估算内孤立波传播速度

基于叠前地震数据估算反射体移动速度的方法已成功应用于第一模态和第二模态内孤立波的传播速度估计上(Tang et al., 2015; Fan et al., 2021；邝芸艳等，2021).我们首先从一系列共炮检距道集COG(Common Offset Gather)剖面中追踪与同一个内孤立波波形相对应的反射同相轴(图2).然后，我们记录了这些同相轴上固定反射点的炮号(#shot)和CMP号(#CMP)，炮号和CMP号可以从叠前地震数据中读取.需要注意的是，选择的反射点必须是不受地震成像影响的特征，例如内孤立波的波谷或波峰(图2中的黑点).反射点炮号的变化代表时间的变化，CMP号的变化代表距离的变化.假设内孤立波保持匀速运动，则内孤立波的传播速度v可表示为

图2 利用叠前地震数据计算内孤立波传播速度示意图(Song et al., 2021b)

(1)

其中dCMP表示反射点CMP号的变化距离量，ds表示反射点炮号的变化量， dt表示放炮时间间隔.

1.3 地震数据估算混合参数

Klymak和Moum(2007)的研究表明，水平波数谱可以表示成内波部分和湍流部分的叠加.其中高波数成分由湍流主导，湍流部分的谱能量与波数的-5/3次方成正比关系.水平波数谱的湍流部分可以通过简化的Batchelor模型表示(式(2))，因此可以根据观测的水平波数谱估算湍流耗散率.另外，湍流扩散率可以根据耗散率由式(3)计算(Osborn, 1980).

(2)

Kρ=Γε/N2,

(3)

式中φT表示水平波数谱，Γ=0.2为混合系数，N为浮频率，CT=0.4为Kolmogorov常数，ε为湍流耗散率，kx表示水平波数，Kρ为湍流扩散率.

(4)

利用测线附近的温盐再分析数据，计算了浮频率随深度变化曲线，用于混合参数的估计.如前所述，再分析数据由欧盟哥白尼海洋环境监测中心提供，这里使用的是日平均数据.

2 结果

2.1 内孤立波、内孤立波包波形特征与相速度特征

根据垂向结构不同，海洋中观测到的内孤立波主要有第一模态内孤立波与第二模态内孤立波两种(Fan et al., 2021，2022).第一模态内孤立波有深水中常见的下凹型内孤立波与浅水中可能观测到的上抬型内孤立波，其等密面、等温面或反射同相轴具有单个下凹或上凸的形态.第二模态内孤立波分为凸型与凹型，实际海洋中常见的是凸型第二模态内孤立波，其等密面、等温面或反射同相轴具有上部上凸、下部下凹的形态.如图3所示，12号地震测线上分布着多个内孤立波，除ISW1为第二模态内孤立波外，ISW2-13均为下凹型第一模态内孤立波.其中，ISW2-10组成一个内孤立波包，我们称其为内孤立波包I.ISW11-13组成的内孤立波包称为内孤立波包II.12号地震测线采集方向由NW向SE，图3地震剖面上显示为自左至右.

内孤立波包II所处水深较大，深度为700～1000 m，三个内孤立波ISW11-13波形相对简单，受水深影响较小(图3，图4)，基于COG计算的传播速度(视速度)分别为2.57, 3.69, 1.94 m·s-1(方向向左-NW向，内孤立波传播方向与地震调查船行进方向相反，以下计算的传播速度方向均向左-NW向)，速度相对较大.ISW2-10内孤立波包(图3，图5)，所处水深310m左右，内孤立波形态较为复杂，基于COG计算的内孤立波传播速度分别为1.35，1.94, 1.68，1.91，1.24, 1.45，1.45, 2.36，2.02 m·s-1.与较深水的ISW11-13的3个内孤立波相比，传播速度较小，这基本符合传播速度与水深的反比关系规律，即水深越大，内孤立波传播速度越大.

图3 东沙海域12号地震测线(a)，(b)、(c)为该测线的两段剖面，分别含编号ISW1-10的内孤立波，编号ISW11-13的内孤立波

图4 含编号ISW11(a)、ISW12(b)、ISW13(c)内孤立波的测线段和含小振幅内波(d)的测线段地震剖面

如表1所示，ISW2-10孤立波包水深310～330 m，ISW11-13内孤立波水深730～900 m，两个波包水深差异较大，ISW2-10速度相对较小，ISW11-13速度相对较大.从计算的线性内波相速度可知，ISW11-13波包中的三个波因为水深的差异其相速度相差并不大，而振幅的影响有效果了，导致ISW11计算的传播速度在三个波中最大.虽然理论计算有这样的认识，但实际内孤立波瞬时速度的观测是非常困难的，也是受多种因素的影响.地震估算的相速度其实质是视相速度，真实速度Vtrue=Vseis|cosθ|，但波传播方向与测线方向的角度θ不是特别确定，如夹角小，视速度与真实速度差异不大，如角度为60°，则真实速度只是视速度的一半.此外，不同COG剖面上特征点的拾取也会有人为的误差.因此，需要有比较多的数据，进行系统统计分析，才会有观测方面较好、较全面的认识.ISW12的地震估算速度偏大，应该算是一个异常数据，假设它传播方向很不同则比较勉强，ISW13水深大速度小则主要是因为振幅小.ISW2-10的9个内孤立波，水深差异不大，振幅差异较大，但也没有速度随振幅增大而增大的清晰证据.两个波包中的内孤立波不稳定，可能是其地震估算速度变化复杂的一个原因.由于篇幅与本文数据有限，关于内孤立波速度的地震估算的进一步分析不再在本文赘述.

表1 内孤立波特征参数

ISW1为第二模态内孤立波(图5d)，水深340 m，基于COG计算的传播速度为0.81 m·s-1.根据邻近站位9月气候态温盐数据计算的第一模态线性相速度为1.41 m·s-1，第二模态线性相速度为0.65 m·s-1，第三模态线性相速度0.42 m·s-1.第二模态非线性系数为0.0039 s-1，内孤立波最大振幅约为20m，第二模态内孤立波理论传播速度为0.73 m·s-1.因此，从底部反射层的移动速度给出的0.81 m·s-1符合理论KdV方程计算的传播速度.

ISW2，ISW5和ISW9三个内孤立波，具有特殊的结构，180 m以上的内部结构较杂乱(图5)，推测可能存在圈闭涡核(trapped core).ISW2，ISW5顶部同相轴有上凸形状的迹象，似乎还可能是第二模态内孤立波，但从计算的传播速度来看，这几个内孤立波比较接近第一模态内孤立波的理论解而不是第二模态内孤立波传播速度的理论解.因此，上凸形状迹象可能是圈闭涡核存在的表现.

内孤立波包II与内孤立波包I之间，水深约500 m、水平距离65 km附近存在陡坎(图3)，较深水的内孤立波或波包通过该陡坎，发生了较大的浅化效应，因此，两个波包波形结构相差较大.根据Shroyer等(2010)的研究，内孤立波浅化可能形成第二模态内孤立波为首波，后续第一模态内孤立波波包为尾波的情形，但Shroyer等(2010)的第一模态内孤立波尾波振幅较小，与这里的观测结果不同，内孤立波包I振幅远远大于第二模态内孤立波ISW1的振幅.

2.2 水平斜率谱特征

对地震数据进行倾角滤波，拾取地震剖面中的同相轴，并根据数据和方法一节所述的方法可以获得同相轴起伏的水平斜率水平波数谱.图6展示了4段地震剖面拾取同相轴的平均水平斜率谱.与滤波前拾取同相轴计算波数谱的结果(这里没显示)不同，滤波后的波数谱计算结果表明，在高波数段即湍流子区间段，谱的斜率满足1/3定律(Kolmogorov幂指律).其中第1段、第2段剖面平均湍流扩散率约为1×10-4m2·s-1，比开阔海洋中的数值(10-5m2·s-1)大一个量级.第3段剖面、第4段剖面平均湍流扩散率约为1×10-3m2·s-1，比开阔海洋中的数值(10-5m2·s-1)大两个量级.且第4段剖面平均湍流扩散率与第3段剖面的扩散率相比稍大.这说明两个孤立波包发育区域湍流能量均非常强，因此湍流扩散率值很大.深水区(>700 m)的内孤立波包II与300～400 m水深处的内孤立波包I，虽然形态差异较大，但引发的湍流混合均很强.

图6 四段滤波后地震剖面(a-d)及其垂直位移水平斜率水平波数谱(e-h)，图(e-h)上三条1/3斜率直线段表示湍流扩散率分别为10-2m2·s-1, 10-3m2·s-1, 10-4m2·s-1的谱能量水平

2.3 湍流扩散率的空间分布特征

测线的扩散率剖面如图7所示，整个剖面的扩散率量级高于开阔大洋的扩散率量级(10-5m2·s-1).高扩散率值主要分布在三个区域，1)水深50～150 m之间特别是内孤立波包发育区；2)水深600～800 m之间，70～110 km处的陆坡海底上方；3)几个主要内孤立波发育区，自浅至深扩散率高值呈现条带状分布.内孤立波包、内孤立波所在区域，扩散率较大，这说明大振幅内孤立波对混合的贡献较大.对两个主要内孤立波包区域进行放大显示(图8，图9)，结果表明ISW2，ISW5，ISW9三个内孤立波所在位置的扩散率均较高，含圈闭涡核的内孤立波增强了混合.第二模态内孤立波ISW1发育区的100～160 m水深范围也有扩散率增大的现象.在海底附近有斑块状高扩散率值分布，可能是这些内孤立波引起的.内孤立波包II的前两个内孤立波ISW11，ISW12所在区域(图9)，垂向上均为扩散率高值，近海底存在横向展布的高扩散率区域，说明一方面内孤立波发育区混合增强，另一方面东沙陆坡海底上方为混合增强带，与内孤立波和海底相互作用密切有关.

图7 12号测线扩散率与地震剖面叠合图

图8 12号测线含内孤立波包I段扩散率与地震剖面叠合图

3 讨论

3.1 含圈闭涡核的内孤立波

含圈闭涡核内孤立波是近期内孤立波浅化研究中的一个热点问题，在数值模拟和观测研究方面均有了些初步成果(如Lamb，2002； Lien et al., 2012, 2014; Chang et al., 2021; Chang, 2021)，一般认为流速超过波形传播速度时，会形成涡核，这是非线性程度较高的表现，也为物质横向输送提供了重要途径.最近的研究认为流速与波形传播速度相当时，也能出现涡核(Chang et al., 2021; Chang, 2021).

Lien等(2012)首次在东沙海域观测到了含圈闭涡核的内孤立波，该内孤立波是在传播速度从2 m·s-1急剧下降到1.3 m·s-1的时候形成的，该文图7—9给出的结果表明，该内孤立波位于173.23°，21°N，水深400 m左右.涡核引起的密度翻转为10～50 m，基于细尺度参数化方法估算的湍流动能耗散率为O(10-4)W·kg-1，涡扩散率为O(10-1)m2·s-1.

Lien等(2014)进一步研究了含圈闭涡核的内孤立波，利用448 m和500 m深的两个锚系及走航资料，观测了五个第一模态下凹型大振幅ISW.研究表明沿波传播方向流速约为2 m·s-1，垂直速度约为0.7 m·s-1，最大垂直位移为167 m，占总深度的0.37.对于破碎波，由下层厚度归一化的最大垂直位移达到0.4，与实验室实验结果一致(Helfrich, 1992).至少有三个内孤立波达到破碎极限准则(Umax>C,Umax为流速，C为波传播速度，参见Zhang和Alford(2015))，表明存在对流不稳定.在波核内出现100 m的垂直翻转，持续达48 min，湍流动能耗散率估计为1.5×10-4W·kg-1.

如图5(a)—(c)所示，ISW2，ISW5与ISW9三个内孤立波核心部位，没有清晰的反射同相轴，与Gong等(2021a)研究内孤立波振幅垂向结构的多个内孤立波形态明显不同，从波形图像来看，最突出的含圈闭涡核特征的内孤立波是ISW9内孤立波.地震数据估算的ISW2，ISW5, ISW9三个内孤立波传播速度分别为1.35，1.91，2.36 m·s-1.从传播速度来看，最有可能的含圈闭涡核的内孤立波是ISW2内孤立波，因为它最容易满足流速大于波形传播速度的条件.图10给出的耗散率分布结果表明，ISW2，ISW5, ISW9三个内孤立波核心部位均为高值区，约为O(10-5)W·kg-1，与Lien等(2012, 2014)的细结构参数方法估算值相比小一个数量级，这与ISW2，ISW5, ISW9三个内孤立波振幅相对较小有关.但这些内孤立波的耗散率、扩散率已经很大了，增强的湍流混合对区域水文、生态以及沉积物悬浮搬运的影响值得重点关注.

图10 12号测线含内孤立波包I段耗散率与地震剖面叠合图

3.2 第二模态内孤立波

值得指出的是，具有复杂波形的内孤立波包I的前方为第二模态内孤立波ISW1.与第一模态内孤立波相比，东沙海域第二模态内孤立波的研究很少.Yang等(2009)首次详细研究了东沙海域第二模态内孤立波，研究所用的主要锚系水深350 m.第二模态ISW通常显示上部(下部)水柱中的等温线向上(向下)位移，以及波致流速的三层结构.夏季第二模态 ISW引起的等温线位移在75 m处为20±14 m，在240 m处为22±15 m，特征时间尺度约为8.0±4.3 min(Yang et al., 2009).在夏季的第二模态 ISW中，90%出现在第一模态 ISW之后，表明第二模态 ISW可能与第一模态 ISW有关.Yang等(2010)在Korteweg-de Vries(KdV)方程的基础上，使用一个解析性三层模型来研究两种(凸型与凹型)第二模态ISW类型的特征，并认为中间层厚度是主要控制因素.因为具有较厚中间层的分层海洋非常罕见，因此凹型波很少出现，这一分析结果与东沙海域观测结果吻合得很好(Yang et al., 2010).

图5d所示120 m为中间层中心深度，上部同相轴向上位移，下部同相轴向下位移，上、下层最大振幅约为20 m.该内孤立波所处水深为340 m.第二模态内孤立波研究(如Fan et al., 2021，2022)认为密跃层深度与水深一半深度偏离程度是衡量该内孤立波对称性、不稳定的一个指标.ISW1偏离程度较高，是一个不对称的凸型内孤立波.图10给出的耗散率分布结果表明，ISW1内孤立波核心部位也为高值区，约为O(10-6)W·kg-1.该第二模态内孤立波与后续的第一模态内孤立波包的先后位置关系，与Yang等(2009)总结的夏季第二模态内孤立波在第一模态内孤立波后面完全不同，进一步说明了东沙海域浅化过程是非常复杂的.

4 结论

南海东沙海域上陆坡区域内孤立波浅化过程复杂，湍流混合时空分布很不均匀.利用高分辨率的地震海洋学方法可以给出内孤立波波形结构和一系列定量参数，估算混合参数的空间分布.通过叠合波形结构与混合参数空间分布，可以较好地分析内孤立波浅化过程与能量耗散、混合增强的关系.

在东沙海域的一条地震剖面上，我们观测到上陆坡深水区(500～1000 m)、浅水区(300～400 m)两个内孤立波波包，有着明显不同的波形特征.深水区的内孤立波包的3个内孤立波波形结构相对简单.浅水区的内孤立波包的3个内孤立波波形复杂，在波核内部呈现杂乱反射特征，可能是含圈闭涡核的内孤立波，估算的耗散率、扩散率分别达O(10-5)W·kg-1、O(10-2.5)m2·s-1.该内孤立波包前方，则存在一个第二模态内孤立波.浅水区的第二模态内孤立波与复杂波形的内孤立波包是深水区第一模态内孤立波包通过500m左右深度处的陡坎演化形成的，东沙海域的陡坎对内孤立波浅化影响较大.

深水区内孤立波包虽然波形相对简单，但诱导的湍流混合强度也较强，形成了内孤立波发育区、陆坡上方厚度达200 m的高耗散率、扩散率区带，估算的耗散率、扩散率分别达O(10-5)W·kg-1、O(10-2.5)m2·s-1.因此，本文研究表明，在东沙海域上陆坡的深水区、浅水区，内孤立波持续耗散能量，增强了湍流混合.