内孤立波具有短周期，强振幅的特点，在世界海洋中分布广泛，特别是在边缘海区域和海岸带的起伏地形地区．一般认为，在有利的海洋条件下，正压潮与海底地形的相互作用是产生内孤立波最常见的机制之一．当正压潮从产生地向外辐射，它的波前面变陡，变得更加非线性，最终导致高频非线性内孤立波的形成．在世界海洋中，南海北部的内孤立波最为活跃．南海东北部深海盆地中内孤立波起源于吕宋海峡西部的M2内潮，而全日潮对内孤立波的生成起着次要作用等（2004）将产生于吕宋海峡的内孤立波分为两类，分别是较大、较窄，每天同一时间到达的type-a型内孤立波，及较弱、较宽，每天晚一小时到达的type-b型内孤立波．

对于南海海盆东北部内孤立波的研究，Huang等（2014）从2010年3月到8月在该区域现场观测到了137个内孤立波，其中type-a型内孤立波大约出现在每天晚上的同一时间，Batan岛南部是type-a型内孤立波的主要产生源，Itbayat岛南部和Batan岛南部是type-b型内孤立波的主要产生源．Xie等（2015）利用MITgcm模拟了孤立的气旋或反气旋对内孤立波的影响，发现涡旋在内孤立波两侧减速和加速的影响会使内孤立波波前面发生扭曲，反气旋还能形成次一级的尾部内孤立波包．Huang等（2016）在该区域捕捉到了一个最大振幅达240m，西向流速为2.55m·s-1的巨型内孤立波，分析其由吕宋海峡半日潮产生，并与Batan岛附近黑潮加强的层结有关．Huang等（2017）研究了中尺度涡旋对吕宋海峡产生、穿过南海北部深海盆的内孤立波的影响，通过观察发现反气旋（AE）和气旋（CE）对会较大改变内孤立波的振幅、传播方向和传播速度．

地震海洋学自Holbrook等（2003）首次提出以来，已经广泛应用于各种海洋现象的调查研究，包括涡旋，洋流，内波，温盐阶梯不稳定，湍流和海底界面过程（陈江欣等，2017；耿明会等，2018）等．Buffett等（2012）提出SOG(地震偏移距组）的方法，将整个检波器排列分成若干组进行叠加，从而获得海水结构在采集时间内变化的图像．Sheen等（2012）将地震反射振幅表示为不同深度切片炮点（代表时间）和CMP(共中心点，代表距离）的函数，从而计算了平面中内波的水平速度．地震海洋学方法在研究水体结构方面具有独特优势（垂向和横向分辨率可达10m），也开始用于内孤立波的研究．如拜阳等（2015）利用地震海洋学方法研究了南海东北部东沙海域内孤立波的结构特征，并比较了内孤立波水平波数谱与GM模型谱，发现其能量在低波数段比GM谱大两个数量级左右．Tang等（2016）使用反射地震数据与局部水文数据通过马尔科夫链-蒙特卡洛反演方法获得了内孤立波的温度和盐度．Bai等（2017）利用反射地震数据系统分析了南海东北部浅化内孤立波的极性转换．黄晞桐等（2018）结合流体动力学数值模拟与反射地震正演分析了内孤立波浅化不同阶段海水层反射几何形态、反射结构等特征．Tang等（2018）用高分辨率地震成像和海水结构特征联合反演研究了南海北部陆架次一级内孤立波，该内孤立波包含两个第二模态内孤立波和一个上凸型内孤立波，并在地震剖面上直接观测到了剪切引起的Kelvin-Helmholtz不稳定．Geng等（2019）用南海北部众多反射地震数据识别出的内孤立波，研究了内孤立波振幅与水深的关系，内孤立波最大振幅随水深呈对数减小．孙绍箐等（2019）利用地震海洋学方法研究了地中海直布罗陀海峡附近内孤立波的结构特征．对于利用地震海洋学计算内孤立波视相速度的研究，Tang等（2014）使用叠前地震数据中内孤立波波谷对应的炮点-检波点对计算了南海东北部陆坡内孤立波的视相速度．为了减小拾取单个内孤立波波谷点的不确定性，Tang等（2015）使用网格搜索方法拾取内孤立波前翼具有最佳斜率的炮点-检波点对，计算了东沙岛东北部内孤立波包的视相速度．

到目前为止，对于南海海盆东北部内孤立波的研究大多是基于物理海洋观测、遥感观测或者数值模拟开展的．物理海洋学的观测方法只能定点对流经海水的物性结构进行刻画．在流向的下游，海水结构可能受其他海洋现象的影响而发生改变，而定点观测无法对此进行记录，这是物理海洋学定点观测的局限．这一局限可以被空间上连续观测的地震海洋学克服．如利用地震海洋学可以观察海水细结构的时间、空间变化．物理海洋学和地震海洋学各有优势，地震海洋学观测到的内孤立波，锚系可能没观测到，两者应该互为补充．以往利用地震海洋学方法发现的内孤立波大多在东沙岛附近．最近，我们利用地震海洋学方法对已有南海海盆东北部的地震数据重新进行了处理，在测线上发现了内孤立波的存在．该内孤立波是目前为止用地震海洋学方法在海盆中首次发现的内孤立波．本文将基于目前所获得的成果和前人的工作，主要研究南海盆地东北部附近内孤立波的细结构及其变化，内孤立波传播特征等．

1、研究区域概况

本文的研究区域位于南海东沙群岛以东，吕宋海峡以西，南海海盆东北部海域，水深1000～3000m（图1）．吕宋海峡是连接南海和西太平洋最深最宽的通道，且有着广泛的海水交换．经过吕宋海峡的水团具有“三明治”的垂向流动模式（Sheuetal.,2009;Tianetal.,2006），中层水向东流，而上、下层水向西流．Tang和Zheng(2011）将南海盆地东北部海域海水分为三层，上层水在550m之上，中间层水在反射地震剖面上具有中等至较弱反射的特点，下层水在大约1000m以下．冬季时黑潮在东北季风的作用下会从吕宋海峡入侵南海，而夏季时南海海水则会流出吕宋海峡（ShawandChao,1994），黑潮部分支流会入侵南海北部形成反气旋涡旋，且有继续由东向西移动的趋势．由于内潮与吕宋海峡的相互作用，东沙群岛附近海域发育有大量内孤立波（拜阳等，2015），而在南海盆地东北部海域则有大量中尺度涡旋发育（李燕初等，2003).

2、数据和方法

本文主要利用地震海洋学方法研究南海海盆东北部发育的内孤立波．地震数据是由MGDS海洋地球科学数据系统提供，依托于TAIGER项目由“R/VLangseth”于2009年4月14—16日内采集取得的二维多道反射地震数据．地震采集参数如下：采样率是2ms，炮间距为50m，道间距为12.5m,CMP间距为6.25m，最小偏移距为179m，每炮468道（大道号对应小偏移距）．海水层地震数据处理流程与常规反射地震数据处理流程基本相同，需要经过解编、观测系统定义、噪声衰减（包含直达波压制）、共中心点（CMP）选排、速度分析、动校正、叠加和叠后去噪等步骤得到叠加剖面，更详细的处理过程可以参考地震海洋学数据处理流程.

图1研究区测线位置分布

本文使用叠前偏移计算内孤立波细结构的变化．该方法的流程如下：首先对切除直达波的叠前数据进行速度分析，从而获得不同共中心点（CMP）位置的叠加速度．然后，对切除直达波的叠前数据抽取共偏移距道集，并利用速度分析获得的叠加速度在COG上做动校正，从而去除不同偏移距剖面上剩余时差的影响．之后，对动校正后的COG进行叠前Stolt偏移，使绕射波收敛，反射波归位，从而更准确的获得内孤立波细结构．最后利用倾角滤波对叠前偏移剖面进行滤波，提高叠前偏移剖面的信噪比．

本次计算内孤立波视相速度，借鉴了Tang等（2014）的方法，但有所改进（图2).Tang等（2014）在COG剖面上拾取内孤立波波谷对应的炮点和检波点对，通过拟合炮点-检波点对曲线，计算得到内孤立波视相速度（图2虚线框中）．需要指出的是原始的COG剖面信噪比较低，给拾取内孤立波波谷对应的炮点-检波点对带来困难．为了更好的计算内孤立波视相速度，本文利用速度分析获得的叠加速度场，对COG剖面进行叠前偏移，然后从信噪比较高的叠前偏移剖面上拾取内孤立波波谷对应的共中心点和炮点对，通过拟合共中心点-炮点对曲线，计算得到内孤立波视相速度以及视传播方向（图2点划线框中）．改进的内孤立波视相速度计算方法假设船速是固定的，如图3a中的坐标系统，内孤立波波谷（波峰）的水平速度v=(CMP2-CMP1）/T=(CMP2-CMP1）/[（s2-s1)*dt]，其中CMP1和CMP2是内孤立波在不同时刻的波谷（波峰）位置，s1和s2是CMP1和CMP2对应的采集炮号，dt是放炮的时间间隔．

改进的内孤立波视相速度计算方法在确定内孤立波视传播方向时，分以下两种情况．当船作业方向与内孤立波视传播方向相反时，无论船速Vship和内孤立波视相速度Vwater相对大小如何变化，随着偏移距的增大，运动过程中同一个内孤立波对应的CMP号减小，炮号增大（图3a).

当船作业方向与内孤立波视传播方向相同时，若Vship<Vwater，随着偏移距的增大，运动过程中同一个内孤立波对应的CMP号减小，炮号减小（图3b）；若Vship>Vwater，随着偏移距的增大，运动过程中同一个内孤立波对应的CMP号增大，炮号增大（图3c).

图2改进的内孤立波视相速度计算流程（点划线框中的，虚线框中是Tang等（2014）的流程）

3、结果与讨论

3.1内孤立波细结构及其变化特征

对地震数据进行解编、观测系统定义、噪声衰减（包含直达波压制（直达波压制后结果见图4b））、共中心点（CMP）选排、速度分析、动校正（动校正后结果见图4c）、叠加和叠后去噪等处理后，本研究在研究区20号地震测线上捕获到了一个波形较完整的内孤立波（图4a），该内孤立波为第一模态下沉型内孤立波，测线位置如图1中红线所示．20号测线上内孤立波前翼的反射同相轴略高（内孤立波视传播方向为从NW向SE，见3.2小节）．该内孤立波处的海底深度为3030m（图4a中地震剖面由于1125m以下基本没有信号，截取到1125m），最大振幅在344.5m处，最大振幅为86.1m（图5），属于大振幅内孤立波，内孤立波振幅随着深度增加呈先增大后减小的趋势．

在COG上做叠前偏移需要速度分析获得的叠加速度，为了指导速度谱能量团的拾取（叠加速度拾取的趋势），我们对20号测线附近的三个XBT计算了均方根速度（三个XBT站点位置见图1．在假设地层水平的条件下，均方根速度等于叠加速度），发现随着深度的增加，均方根速度呈现减小的趋势（图6),这与常规海底地层拾取叠加速度时，叠加速度随地层深度增加而增加的趋势不同．

图3计算内孤立波相速度示意图

(a）船与内孤立波视传播方向相反；（b）船与内孤立波视传播方向相同，且船速小于内孤立波传播速度；（c）船与内孤立波视传播方向相同，且船速大于内孤立波传播速度．五角星和圆圈分别是震源和检波器，它们的共中心点（CMP1和CMP2）位于内波的波峰，虚线箭头表示以s2为参考点的坐标系．

图7a和图7c分别是CMP12060和CMP12500的速度谱，在速度谱上根据速度谱能量团和所得的该区域均方根速度变化趋势分别拾取相应的叠加速度（图7b和图7d），发现在内孤立波出现的位置（CMP12500）叠加速度有明显的增加．进一步将获得的叠加速度转成层速度．图8a是20号测线CMP9000-16000范围进行速度分析得到的层速度（经过平滑处理），图8b是对内孤立波部分进行速度分析得到的层速度（经过平滑处理）．可以看到内孤立波两翼的速度结构并不对称，前翼部位的速度场被抬升（内孤立波视传播方向为从NW向SE，见3.2小节）．

利用拾取的叠加速度场对COG剖面做叠前偏移，其结果如图9(a,b）所示．图9a是道号等于441对应COG的叠前偏移剖面，属于近偏移距的COG．图9b是道号等于409对应COG的叠前偏移剖面，属于远偏移距的COG．图9c是道号等于409对应的COG剖面（为便于与图9b比较，对原始COG剖面做了动校正，但未做叠前偏移）．比较图9b和图9c，可以看到偏移后剖面的信噪比得到较大提高，更易于拾取内孤立波波谷对应的共中心点和炮点对．由于本次研究水层中反射同相轴较平（图9c），偏移后的反射同相轴位置基本没有发生变化（图9b）．比较图9a和图9b发现在该两个偏移距的时间尺度下，随着时间推移，浅层（225m以上）的海水层混合得更加均匀，因而图9b浅层部分的反射较弱；图9a和9b中深层海水反射的特征变化不大．浅层海水相对深层海水反射变化较大，这与海水中温跃层变化剧烈的认识较一致．

3.2内孤立波的相速度特征

在计算内孤立波视相速度时使用如下公式．其中va是内孤立波的视相速度，dcdp是CMP间距，dt是放炮的时间间隔，a是拟合的CMP-炮点对曲线斜率．20号测线的道间距是12.5m，炮间距是50m．图10a和图10b分别展示了利用共偏移距叠前偏移去噪后的剖面拾取内孤立波波谷对应的CMP号和炮号．图10c是利用共偏移距叠前偏移去噪后的剖面拟合的CMP-炮点对曲线，计算的内孤立波视相速度va=0.82m·s-1，其误差为±0.28m·s-1．因为放炮的时间间隔等于20.33s，船速vship=2.46m·s-1,vship>va，随着偏移距增大，CMP号减小，炮号增大．所以船从SE向NW行驶（352°N方向，自北顺时针为正），内孤立波视传播方向与船运动方向相反（内孤立波视传播方向为172°N方向，自北顺时针为正）．

图10d是直接利用COG拟合的CMP-炮点对曲线，对比共偏移距道集叠前偏移剖面拟合的CMP-炮点对曲线（图10c），其线性较差．分析是由于原始COG信噪比较低，不能准确拾取CMP-炮点对造成的．为了验证共偏移距道集叠前偏移剖面计算的内孤立波视相速度准确性，我们利用KdV方程计算了该内孤立波理论的传播速度．具体过程如下：选取GDEM-Version3.0数据站位（119.5°E,20.5°N）的温盐数据计算浮频率，求解KdV本征方程，得到内孤立波各个模态的线性相速度C．站位（119.5°E,20.5°N）的前三个模态的线性相速度如图11a所示，第一模态线性相速度最大，为1.85m·s-1，后面模态的线性相速度减小．得到浮频率和线性相速度，就可以计算各个模态的本征函数．前三个模态的本征函数如图11b所示，我们使用第一模态的本征函数和线性相速度计算非线性系数α=-0.01s-1和非静力频散系数β=7.41×104m3·s-1，该内孤立波的最大振幅为86.1m，则内孤立波传播速度V为2.11m·s-1．该结果比利用共偏移距道集叠前偏移剖面计算的视相速度0.82m·s-1大很多，分析是由于20号测线上内孤立波西南和东南部发育有气旋（见下文），当内孤立波与气旋相互作用会使得内孤立波传播速度降低（Huangetal.,2017），而KdV理论没有考虑背景流（如涡旋）的影响，这有可能使其计算的内孤立波传播速度偏高．

图4（a)20号测线内孤立波部分地震叠加剖面（经过叠后FK倾角滤波，只显示了54～98km的水平范围以及0～1125m的垂向范围）．测线采集时间为2009年4月14日—2009年4月16日；（b）对炮集（2460炮）切除直达波；（c）对CMP道集（CMP12260）进行动校正，红线标记出了动校拉伸切除位置

图4（a)20号测线内孤立波部分地震叠加剖面（经过叠后FK倾角滤波，只显示了54～98km的水平范围以及0～1125m的垂向范围）．测线采集时间为2009年4月14日—2009年4月16日；（b）对炮集（2460炮）切除直达波；（c）对CMP道集（CMP12260）进行动校正，红线标记出了动校拉伸切除位置下载原图

图520号测线内孤立波振幅随深度变化

Huang等（2014）利用2010年3月21日到2010年8月12日在吕宋海峡西侧深海盆地布置的锚系A2(21.1°N,119.1°E）观察到了137个内孤立波．在假设内孤立波运动年际变化不是很大的前提下（20号测线采集时间为2009年4月14—16日），其数据采集时间和地点与本研究基本相近．Huang等（2014）将137个内孤立波分成81个type-a类型内孤立波和56个type-b型内孤立波．Huang等（2014）发现所有的type-a型内孤立波在晚上（大约19∶00）被检测到，而type-b型内孤立波在每天早上到中午（大约5∶00—14∶00）被检测到．20号测线上的内孤立波采集时间约为2009年4月15日上午时间段（大约6∶35—8∶00），此时该区域type-b类型内孤立波盛行．初步推测20号测线捕捉到的内孤立波为type-b类型的内孤立波．

图11c中红色曲线是使用内孤立波延伸深度（H=1000m）通过KdV方程计算的内孤立波振幅的垂向变化，垂向模态最大值的深度为344m，与内孤立波最大振幅深度一致，振幅变化趋势与理论也较接近．说明在20号测线采集时间段内type-a型内孤立波和type-b型内孤立波振幅的垂向变化基本一致，因而未出现由于“平均效应”（KdV方程计算时使用的温盐数据是全天的，其计算的内孤立波振幅垂向变化是全天出现的type-a型内孤立波和type-b型内孤立波振幅垂向变化的平均，只有当采集时间段内type-a型内孤立波和type-b型内孤立波振幅的垂向变化基本一致，实际观测的type-b型内孤立波振幅垂向变化才和KdV方程计算的平均内孤立波振幅垂向变化基本一致）产生理论与观测结果存在较大误差的结果．

Huang等（2014）在没考虑涡旋影响的条件下，给出研究区type-b型内孤立波的平均传播方向为285°N（自北顺时针方向为正）．用共偏移距道集叠前偏移剖面计算的内孤立波视传播方向为从NW向SE(172°N方向，自北顺时针为正）．而内孤立波视传播方向是真实传播方向在地震测线上的投影，所以该内孤立波真实传播方向介于172°N—262°N之间（自北顺时针为正），其与Huang等（2014）给出的平均传播方向（285°N方向，自北顺时针为正）有差异．Huang等（2017）在研究反气旋和气旋对内孤立波传播路径影响时，将其分成四个阶段．其中第一个阶段内孤立波波前面中部位于反气旋的南部，此时内孤立波波前面的北部和南部分别折向南海东北部的北面和南面．在第四个阶段，内孤立波波前面中部位于气旋的北部，此时内孤立波波前面的北部和南部也分别折向南海东北部的北面和南面．我们利用CMEMS的物理海洋学模型，获得了该地区在2009年4月15日的地转流速（图12）．观察地转流速的分布，发现沿着20号测线，在内孤立波西北部发育有较强的反气旋（涡旋），西南和东南部发育弱的气旋．对应Huang等（2017）研究中反气旋和气旋对内孤立波传播路径影响的第一阶段和第四阶段的情况，即内孤立波波前面中部位于反气旋的南部，气旋的北部．此时，内孤立波波前面的南部折向南海东北部的南面．因而判断此处内孤立波位于内孤立波波前面的南部，使得其真实传播方向（介于172°N—262°N之间的方向，自北顺时针为正）有偏差．图13是CMEMS物理海洋模型给出的对应20号测线位置及采集时间的温度剖面，发现温跃层在西北侧变深，说明此处涡旋对包括内孤立波在内的海水温度场产生了影响，间接证明了本研究区内孤立波传播特征在很大程度上受涡旋的影响．

图6（a）T7＿39计算的均方根速度；（b）T7＿40计算的均方根速度；（c）T7＿41计算的均方根速度

图7（a）、（b）分别是CMP12060的速度谱和拾取的叠加速度；（c）、（d）同（a）、（b）但对应CMP12500

图8（a）CMP9000-16000范围速度分析得到的层速度；（b）CMP11060＿13694范围速度分析得到的层速度，背景是对应的地震叠加剖面

4、结论

以往利用地震海洋学方法发现的内孤立波大多是在东沙岛附近．本文首次利用地震海洋学方法在南海海盆东北部发现了海盆中的内孤立波，并对其进行了研究．该内孤立波在地震叠加剖面上表现为第一模态下沉型内孤立波，其最大振幅为86.1m，属于大振幅内孤立波．通过叠前偏移观察该内孤立波细结构的变化，发现内孤立波波形在采集时间段内整体较稳定，内孤立波浅层反射相对深层变化较大．

通过改进前人的方法，利用共偏移距道集叠前偏移剖面计算内孤立波视相速度．该方法比直接使用共偏移距道集拟合的共中心点-炮点对曲线线性更好．利用共偏移距道集叠前偏移剖面计算的该内孤立波视相速度为0.82m·s-1，内孤立波视传播方向为从NW向SE(172°N方向，0°指向北），其真实传播方向（介于172°N—262°N之间的方向，0°指向北）与Huang等（2014）给出的该地区内孤立波平均传播方向（285°N方向，0°指向北）有差异．通过KdV方程计算的内孤立波传播速度为2.11m·s-1，明显大于反射地震得到的视相速度（0.82m·s-1），而内孤立波垂向结构变化与理论值基本一致．20号测线上的内孤立波西北部有反气旋（涡旋），西南和东南部有气旋，反气旋和气旋对内孤立波传播路径产生影响，且此处内孤立波位于内孤立波波前面的南部，使得其真实传播方向（介于172°N—262°N之间的方向，0°指向北）有偏差．

由于判断该内孤立波传播路径受涡旋影响等是基于前人研究的成果，还需要结合遥感、物理海洋资料进一步研究，或开展相关的数值模拟工作．

图920号测线叠前偏移观察内孤立波细结构变化

图10（a）利用共偏移距叠前偏移去噪后的剖面拾取内孤立波波谷对应的CMP号；（b）利用共偏移距叠前偏移去噪后的剖面拾取内孤立波波谷对应的炮号；（c）是对COG进行叠前偏移拟合的CMP-炮点对曲线；（d）是直接利用COG拟合的CMP-炮点对曲线

图11（a)KdV计算的线性相速度；（b）不同模态对应的特征函数，其中黑线对应的是第一模态内孤立波，红线对应的是第二模态，绿线是第三模态；（c)KdV拟合内孤立波振幅，黑色线段是内孤立波振幅的垂向变化，红色曲线是使用内孤立波延伸深度（H=1000）获得的理论振幅分布，蓝线曲线是使用海底深度（H=3030m）获得的理论振幅分布

图1220号测线采集时对应的地转流速及海平面高度

图13CMEMS物理海洋模型给出的20号测线温度剖面

参考文献:

[1]拜阳,宋海斌,关永贤等.2015.利用地震海洋学方法研究南海东北部东沙海域内孤立波的结构特征.科学通报,60(10):944-951.

[2]陈江欣,宋海斌,关永贤等.2017.海底冷泉的地震海洋学初探.地球物理学报,60(2):604-616.

[3]耿明会,宋海斌,关永贤等.2018.南海北部雾状层分布和特征的地震海洋学研究.地球物理学报,61(2):636-648.

[4]黄晞桐,宋海斌,关永贤等.2018.基于流体动力学数值模拟的海水层反射地震研究.地球物理学报,61(7):2892-2904,doi:10.6038/cjg2018L0382.

[5]李燕初,蔡文理,李立等.2003.南海东北部海域中尺度涡的季节和年际变化.热带海洋学报,22(3):61-70.

[6]孙绍箐,张锟,宋海斌.2019.地中海直布罗陀海峡附近内孤立波的地球物理特征.地球物理学报,62(7):2622-2632.

[7]王浩森,徐行,刘青松,仲义,陈艇,强小科,周祐民,杨小强.四万年以来南海西北次海盆地磁场古强度演化及其地球动力学意义[J].地球物理学报,2020(07):2671-2681.

范文豪,宋海斌,张锟,许鹤翰,龚屹,孙绍箐.南海海盆东北部内孤立波的地震海洋学研究[J].地球物理学报,2020,63(07):2644-2657.

基金:国家重点研发计划(2018YFC0310000);全球变化专项国际合作项目(GASI-GEOGE-05);国家自然科学基金(41976048,41576047)资助.