海洋内波是指发生在密度稳定层化的海水内部的一种重力波，是各种尺度海水运动过程中的一个重要环节.探讨内波的产生、分布和衰减，进一步认识其形成机制和能量耗散机制，是20世纪70年代以来国内外物理海洋学研究的热点.而内波研究的进展和水平很大程度上取决于所获取观测资料的数量和质量，但无论传统的接触式温-盐-深（CTD）测量还是卫星遥感等观测手段，都无法在短时间内对整个海水垂向剖面的温盐结构进行成像，给出内波波动的空间形态^[1].2005年，Holbrook等^[2]在利用反射地震方法研究挪威海海水温盐结构时，发现反射层存在着波长从几十米到几千米，振幅几十米的类似正弦波的起伏.经过分析认为，这种波状起伏反映了海洋内波的形态特征；计算其水平波数能量密度谱（以下简称水平波数谱或内波谱）后，发现其谱形与著名的Garett-Munk内波模型的GM76模型谱^[3~5]基本一致，从而验证了反射地震方法在海洋内波研究中的可行性.

南海东北部具有复杂的地形和潮流分布，内波非常发育，被国际公认为优良的天然内波试验场^[1].对南海东北部内波已有较多研究，但定量分析内波特征的研究工作则较少^{[1, 6~12]}.最近我们对国家基础研究发展规划（“973”）项目“中国边缘海形成演化及其重大资源的关键问题”在南海东北部采集的一条地震剖面（包括海洋水层部分）进行重新处理，获得水柱的反射结构剖面^[13]；通过分析，确认地震叠加剖面上同相轴呈现的起伏变化反映了内波的总体形态.该地震数据给出了低频反射地震可以对内波清晰成像的新的证据，从地震剖面上，能够清楚地观察到内波所造成的反射特征变化^[13].基于该海水层地震反射数据计算的水平波数能量密度谱与GM76模型谱基本一致，但在低波数段和高波数段中，两者的振幅及斜率存在着一定差异；经分析认为这种差异主要与内潮波和复杂海底地形的强烈非线性相互作用以及内波破碎等因素有关^[14].

目前，地震海洋学研究内波主要基于波数谱的特征研究，本文拟利用经验模态分解（Empirical Mode Decomposition，EMD）^[15]方法对我国南海东北部海域的内波特征进行研究.EMD是Huang等（1998）^[15]提出的一种新的数据分析方法，与HSA（Hilbert Spectrum Analysis）算法结合，被美国国家航空与航天局（NASA）称为HHT（Hilbert Huang Transform）方法.这种分析方法不仅能够像小波变换那样对信号进行时频分析，同时也具有类似于傅里叶变换那样的谱分析能力；更值得一提的是，HHT还是一种自适应数据分析方法，并且针对非线性、非平稳信号具有非常好的应用效果^{[16, 17]}.

Huang等在其1998年的论文^[15]中对“瞬时频率”进行了深入讨论.瞬时频率通常可通过Hilbert变换的方法得到，然而，并非所有的信号（或数据序列）通过这种方法都可以得到有物理意义的瞬时频率，只有满足“窄带”要求的信号才可以，这里所谓的窄带具体是指信号的极值个数要与其零值个数相等.但不幸的是，自然界中很少有信号是天然满足窄带要求的，于是就必须对信号做某种处理以使其满足这个“窄带”要求.Huang等^[15]同时也给出了解决办法，他们认为信号都是由一系列的固有模态函数（Intrinsic Mode Function，IMF）组成的.IMF的定义如下：（1）对于整个数据序列来说，其极值个数应与跨零点个数相等，或至多相差一个；（2）对于序列的每一点来说，由所有极大值组成的上包络线和由所有极小值组成的下包络线的平均值为0.而IMF的获得，就是由经验模态分解来完成的.具体实现方法如下：

使用三次样条插值的方法将所有的极大值连成与数据序列长度相等的曲线，此为数据的上包络线；对所有的极小值进行同样操作，得到下包络线.上、下包络线应该将所有的数据都包含在其范围内，其平均值记为m₁，并将数据与m₁的差值记为h₁：

理想情况下，h₁应该就是IMF，然而事实情况并非如此.通常来说，h₁并不满足IMF的定义，具体来说就是在h₁中仍然有很多比零小的极大值，和比零大的极小值，因此也就不能满足极值点数目与跨零点数目相等或至多相差一个.这时需要将h₁当作数据，重复以上步骤：

Huang等人给出了这个迭代过程的终止条件：定义标准偏差SD（StandardDeviation）为

当SD小于某一值时，h_1k就是我们所要的IMF了.此时，记

将c₁从原数据中分离出来，得到：

一般来说，r₁仍然可以再分解，重复以上分解过程，有：

有两个条件可以给出分解结束的判断：

（1）c_n或者r_n变得非常小，小于预定值时，分解结束；

（2）r_n变为单调函数，也就是说从r_n中无法再分离出IMF分量，分解结束.

以上两个条件只要满足其一，就可以判定分解结束.最终，我们可以得到：

本文实际计算中，SD以及c_n阈值均取0.05，最大迭代次数设定为100.

2001年8月4~6日，国家973“中国边缘海形成演化及其重大资源的关键问题”项目委托广州海洋地质调查局，利用“探宝号”调查船在南海东北部进行了地震数据采集，其中一段测线的位置如图 1所示，总长463 km.用道间距为12.5 m的240道水听器记录地震数据，采样间隔为2 ms，记录长度为10 s.震源采用容量为3000 in³套筒枪阵，炮间距为50 m，近炮检距为250 m.对该测线进行重新处理获得南海东北部地震水柱成像剖面，发现研究区有大量的热盐细结构与内波发育^[13].层状结构在南海海盆上方最显著，在上陆坡、增生楔（增生杂岩）、吕宋火山弧处内波非常发育，指示海底地形对内波形成与演化的影响^[13].本文利用南海海盆上方CMP37000~39000段（总长12.5 km）的水柱反射剖面（图 2）进行南海东北部内波特征的EMD方法应用研究.

图 1

Fig. 1

图 1 南海东北部位置图 粗黑线代表地震测线的位置，测线上黑点处的数字表示反射地震数据处理的共中心点（CMP）序号. Fig. 1 Bathymetry of the northeastern South China Sea Thick black line indicates the position of the seismic profile, and the CMP numbers are also shown along the line.

用董崇志等^[14]给出的方法拾取了图 2所示水层中的57条起伏反射的同相轴.反射同相轴拾取方法如下：采用Strata软件的同相轴自动追踪模块，以6.25 m的水平间隔（CMP间距）对一系列反射层进行数字化.在用户选定的时间窗（我们选用的是15 ms）内，从初始解释反射点向左右两边自动追踪波形，选取最近的波峰或波谷.在拾取层位过程中，采取人机结合的方式，针对自动拾取层位的不合理跳跃进行手工校正，使其符合反射层的真实位置.拾取的反射同相轴数据给出每一个反射层的双程走时（单位ms），假定地震波在海水中的传播速度为1.5k m/s，可得到相应的深度数据（单位m）.对每个反射层的深度数据取平均，认为是反射层的平均深度，然后利用各层的深度数据减去相应层的平均深度，并且认为通过这种距平计算所得到的就是海水各个质点的位移数据.因为平均深度的分布是不均匀的，为了对位移数据进行二维成像，需要对平均深度插值来使其均匀分布，这里采用的是线性插值方法.最终获得垂向100层的二维位移图如图 3所示.

图 2

Fig. 2

图 2 南海东北部地震反射剖面揭示的水柱内起伏反射---内波剖面 紫红色线为拾取的反射同相轴，D1~D4为拾取的上部4个反射层位；TWT-双程走时，CMP Number-共中心点号. Fig. 2 Water-column undulating reflections from one reflection seismic section for showing internal waves in the northeastern South China Sea Aubergine curves indicate picked reflection events, D1~D4 show upper four picked reflectors, TWT:Two-way Travel Time, CMP Number:Common Midpoint Number.

图 3

Fig. 3

图 3 根据地震剖面计算获得的内波垂直位移分布 Fig. 3 Distribution of vertical displacements of the internal wave field calculated from the seismic section

图 3给出根据地震剖面计算获得的内波垂直位移分布（位移向下为正）.这一区域垂直位移变化于-35~40 m之间.考虑到共中心点叠加技术的平均效应，真实的位移最大值可能还稍高一些.垂直位移沿这12.5 km剖面呈现较大变化.由于海水层深部资料有限，我们无法确定内波的确切模态，但可知剖面右侧比左侧存在内波高模态的效应.在CMP37000附近，位移以正的为主，在CMP37800附近，位移随着水深增大由正变负，在CMP38500左右，位移随着水深增大由负变正，而在CMP39000附近，位移随着水深增大由负变正再变负，显示高模态的变化特征.

图 3只能大致显示海水内波的一般特征，为了进一步提取数据中的信息，我们采用EMD方法来对各层位移数据（插值前57个反射同相轴、插值后的100层）进行分解，这样就可以看到海水位移在不同尺度的波形特征.利用EMD方法对数据进行分解，所得到的各个分量的波数分布范围大致是以二进的方式递减，并且先分离出来的分量主要包含原始信号中的高波数信息，后面分量的波数逐渐降低. 图 4为前四层数据的EMD分解图，IMF1到IMFn表示对原始数据分解所得到的各个分量，而最后一个曲线表示的就是原始数据减去所有分量所得到的残差（RES），一般来讲，这个残差所代表的是原始数据的整体变化趋势.另外需要说明的一点是，利用EMD方法对数据进行分解所得到的IMF个数完全是由数据本身的性质所决定的，并不是所有分解都会得到相同个数的IMF，图 4的四个反射层（D1~D4）数据分解所得到的IMF都是8个，这只是一个巧合而已.从图 4中可以看出，各层数据的IMF7和IMF8分量，以及第2层和第4层的IMF6分量，都呈现出比较好的类似于正弦振动的波形特征.

图 4

Fig. 4

图 4 4个反射层（D1~D4）的经验模态分解结果 Fig. 4 Empirical mode decomposition results of 4 reflectors (D1~D4)

为了观察相近尺度分量的二维成像特征，我们需要将相近尺度的IMF分量组合在一起进行成像，但是由于各层数据分解出的IMF分量个数并不完全相同，所以这里实际采取的做法是从高波数至低波数取四个分量（IMF1~IMF4）分别进行成图（图 5），然后从较低波数至最低波数取四个分量（即IMFn-2，IMFn-1，IMFn，RES）分别进行成图（图 6），但这里需注意一点，从较低波数至最低波数时，将趋势值作为第一个分量.高波数分量部分（图 5）位移值较小，图像凌乱，没有给出有意义的信息.

图 5

Fig. 5

图 5 基于地震剖面得到的垂直位移的经验模态分解结果（高波数部分） （a）IMF1；（b）IMF2；（c）IMF3；（d）IMF4. Fig. 5 Empirical mode decomposition results of vertical displacements from seismic section (High wave-number parts)

图 6

Fig. 6

图 6 基于地震剖面得到的垂直位移的经验模态分解结果（低波数部分） （a）IMFn-2；（b）IMFn-1；（c）IMFn；（d）RES. Fig. 6 Empirical mode decomposition results of vertical displacements from seismic section (Low wave-number parts)

低波数部分（图 6）给出了南海东北部南海海盆上方区域有关内波较多的信息，图 6a呈现11个正的和其间10个负的垂直位移条带，图 6b给出5~6个正的和其间4~5个负的垂直位移条带，分别反映1.2 km与2.5 km波长的内波组成特征.图 7a~7c的波数谱给出了各分量的主波长特征，证实了这三个分量主要由1.2 km与2.5 km波长的内波组成.在900 m的深度范围内，这些波长的内波上、下波动基本耦合，表现出上、下基本一致的形态.图 6c和图 6d上部和下部有较大变化，分界面分别为600 m和700 m.图 6c下部由三正三负的垂直位移带组成，图 6d上部由一正一负的垂直位移带组成，分别反映波长约4 km与12.5 km的内波组成特征.图 7d的波数谱给出了这一分量主要由主波长约4 km的波组成（特别是600 m下方）的依据.而12.5 km的波长则是根据图 6d的波形估算的.图 6c和图 6d上、下部分的内波相位差90°左右，指示低波数内波能量的斜向传播特征.

图 7

Fig. 7

图 7 经验模态分解获得的垂直位移低波数部分的波数谱 （a）IMFn-3；（b）IMFn-2；（c）IMFn-1；（d）IMFn. Fig. 7 Wave-number spectra of low wave-number parts of vertical displacements by empirical mode decomposition

利用地震海洋学方法^{[13, 18~20]}研究海洋内波成为海洋地球物理学家与物理海洋学家共同关注的前沿问题.近年的研究主要在从起伏反射同相轴的拾取基础上，研究内波的波数谱，分析内波能量量级，对比与模型谱的相同与差异^{[2, 14, 21]}.内波是最重要的物理海洋现象之一，它对水下作业、海上设施等都有着重要的影响，因此了解内波的形态与分布特征具有重要的意义.然而传统的海洋地震剖面中，不同尺度的波叠合在一起，不利于对内波结构的精细认识.本文尝试利用当今时频分析的新手段---希尔伯特-黄变换中的经验模态分解方法对南海东北部地震数据处理获得的垂直位移数据进行分解，获得了颇有新意的结果.应用经验模态分解方法对南海东北部内波地震海洋学数据研究表明：南海海盆上方区域的内波包含波长约1.2、2.5、4、12.5 km的组成成分；其中波长约1.2、2.5 km的内波在200~1050 m的深度范围内上、下各层波动基本耦合，而波长约4 km与12.5 km的内波则以600~700 m的深度为界，上、下部分的内波相位差90°.这些研究表明，EMD方法可以清楚地揭示不同尺度的内波结构特征，这将有助于我们进一步认识内波的运动学特征；EMD方法有良好的应用前景，今后的工作是进一步?用二维EMD分解方法来处理地震海洋学数据.

感谢加拿大Dalhousie大学Ruddick教授、葡萄牙Aveiro大学Pinheiro教授、Lisbon大学的Ambar教授、Matias教授在地震海洋学合作研究中的帮助；感谢美国WoodsHole海洋研究所的黄瑞新教授和中国科学院南海海洋研究所甘子钧研究员、王东晓研究员、尚晓东研究员、周伟东研究员等在物理海洋学方面的指点.感谢973项目首席科学家李家彪研究员的支持与帮助.