张佳丽，张安民，孙朝辉，张学峰，张 亮*

(1.天津大学 海洋科学与技术学院，天津 300072；2.自然资源部 第二海洋研究所, 浙江 杭州 310012)

0 引言

内波指发生在密度稳定层结流体中、最大振幅出现在流体内部的波动，是一种普遍存在的自然现象，这种自然现象存在于海洋内部而不易被观测到[1-2]。内孤立波是一种特殊的内波，具有振幅大、传播速度快和非线性的特点[3]，能够进行单向长距离传播。内孤立波经过时，水体的垂向流速、水平流速和垂向剪切均很强，对海上石油钻机和水下航行器安全运行造成了威胁。因此，开展内孤立波的研究对海洋科学研究、海洋资源开发和国防安全保障均具有重要意义[4]。GUO et al[5]综述了近10 a来国内外对内孤立波及其相关动力学过程的研究状况，对遥感图像、现场测量、数值模拟 3种研究方法进行了总结归纳并分析了中国南海的内孤立波，其中现场测量的方法大多使用了ADCP( Acoustic Doppler Current Profiler，声学多普勒流速剖面仪)。ADCP是对内波进行声学监测的有效手段，利用它测得流速、水回波信号强度等数据，可以反演出内波信息[6]。HUANG et al[7]使用ADCP和CTD(Conductivity-Temperature-Depth, 温盐深仪)，通过现场测量捕捉到了南海的一个极端内孤立波，并进行了成因分析。XIE et al[8]使用ADCP对河口重力流产生的内孤立波破碎过程进行了研究。然而，在实际使用中，ADCP的测量精度在很大程度上取决于海洋环境条件及船载辅助设备记录精度，特别是走航式ADCP受到噪声的影响，获得的数据十分复杂[9]，测得流速数据中包含有各种测量噪声和一些流速异常值。移动中值、三阶平滑低通滤波、窗函数及几种滤波方法组合的综合滤波方法是目前已有的滤波方法[10]，但是对于混有流速异常值及各种噪声误差的ADCP数据，这几种方法均不能很好地实现自适应滤波处理，存在截止频率设置不合理导致滤波失真、流速异常值剔除效果不理想等问题。在一些ADCP处理软件中，流速异常值需要人肉眼观察判断，根据测区的实际情况，利用计算机设置相关检测阈值来识别和剔除[11]。

1 抗差Vondrak滤波方法研究

1.1 传统Vondrak滤波方法

Vondrak滤波方法需要综合考虑拟合度与平滑度两方面，其中拟合度反映观测信息，而平滑度反映非观测信息。拟合度与平滑度公式，用最小二乘法求解。对于一组有N个观测值的数值序列x(ti)，ti为测量时刻(i=1,2,...,N)，Vondrak滤波需要满足准则Q：

(1)

其中：Q→min表示Vondrak滤波准则达到最小；F表示滤波结果与原始观测数据的拟合度；S表示滤波结果的平滑度；ε∈(0,+∞)，是一个无量纲的正数，定义为平滑因子。

Vondrak滤波的本质是通过确定合理的平滑因子，达到最佳滤波效果。当ε→∞时，必须有F→0，即此时滤波结果与原始观测数据几乎重合，而当F=0时，滤波结果与原始观测数据完全重合，失去滤波功能; 当ε→0时，必须有S→0,即此时平滑度趋向于0，滤波结果近似为抛物线，导致滤波过度。HVF法是通过反复迭代计算拟合函数与平滑函数的方差，用两者比值反复修复权系数来确定平滑因子,其基本原理如下：

(1)建立两类误差观测方程，其中V1为拟合度误差，V2为平滑度误差：

(2)

(2)令两类观测误差的权值分别为P1和P2，根据最小二乘原理得到：

(3)

将式(2)代入式(3)得到误差方程：

(4)

式(3)的法方程为：

(5)

式中：

(6)

(7)

其中：B1为n阶单位阵,f1=y,f2=0,B2为一带型矩阵[13]。

(8)

(9)

(10)

(11)

其中：C为任意常数。

(12)

其中：P为初始权值，这里取P=1。

1.2 抗差Vondrak滤波方法

在实际ADCP观测数据中，粗差的存在会严重影响数据处理结果的精度，在ADCP测量粗差点，即流速异常值附近的数据拟合结果是不可靠的，因此，需要首先对粗差进行处理。抗差Vondrak滤波方法分为两步，首先基于IGG3方法去除原始观测数据中的流速异常值，然后采用HVF方法进行时域滤波得到平滑的流速曲线。

IGG3算法是一种加权迭代的方法，根据残差和有关参数，按照所选的权函数计算每个观测值在下一次迭代中的权。经过计算，含有粗差的观测值的权会越来越小或降为0。IGG3的算法模型为：

(13)

(14)

其中：V为观测值残差;σ为方差因子，这里用均方根误差表示为：

(15)

(16)

其中：P为初始权值，这里取P=1。

当海流观测中不含流速异常值时，仅采用传统的Vondrak滤波方法就可以得到很好的平滑曲线，船载ADCP在测量过程中出现的噪声可以被有效滤除。当流速异常值存在时，仅通过平滑因子来控制滤波性能的传统Vondrak滤波方法不能对流速异常值进行有效的剔除，在流速异常值附近，数据拟合结果会出现偏差，这对之后内孤立波的提取与分析造成不良影响。为了改善滤波效果，引入IGG3方法对原始数据进行流速异常值自适应剔除，基本操作步骤如下：

(2)将步骤(1)中计算得到的值带入IGG3方法的权因子函数中，得到各观测值的新权；

(3)通过新权与原始观测数据相乘得到剔除流速异常值的观测数据。

1.3 经典滤波方法

FFT、小波分析和滑动平均是3种经典的滤波降噪方法。FFT[20]是离散傅里叶变换(Discrete Fourier Transform, DFT)的一种快速算法，N点序列x(N)的N点DFT可以表示为：

(17)

FFT将原始信号变换到频域中进行分析，根据原始信号的频谱，在频域空间对噪声成分进行有效抑制，再将结果进行傅立叶逆变换，得到降噪后的时域信号。FFT方法对周期性明显的信号降噪效果较好，而对于非平稳、非周期信号，特别是有异常值信号的降噪效果较差。此外，FFT方法需要利用信号的全部时域信息进行整体变换，缺乏时域分析功能，且信号不能同时在时域和频域准确定位。

小波分析[21-22]能有效地从原始信号中提取有效信息，并通过伸缩和平移对信号在时域和频域进行多尺度分析。降噪是小波分析在信号处理领域的主要应用之一，该方法的降噪步骤如下：(1)对信号进行N层小波分解，根据小波分解的特点，将信号分为低频部分和高频部分；(2)选择1个阈值，对高频系数进行抑制；(3)将处理过的系数进行重构，得到降噪后的信号。小波变换在时域和频域均具有很好的局部化特征，但是在小波变换方法降噪的过程中，小波基函数、分解层数以及阈值的选择都很关键，直接影响到对真实信号的降噪效果。

滑动平均是一种简单的数据平滑方法，通过取N个相邻数据的平均值来表示所取数据的端点或中点值。对于N个含有噪声的实测数据{Zi}，取n个相邻数据的平均值，即滑动窗大小为n，滑动平均的数学表达式为：

(18)

其中：n=2m+1；k为滑动平均后的序列项，k=m+1,m+2,...,N-m。

滑动平均方法计算简单，处理数据速度快，但是需要对滑动窗口进行选择。

2 对ADCP数据的处理

2.1 传统Vondrak滤波方法与经典滤波方法对比

为了检验Vondrak滤波方法在海洋内孤立波提取中的适用性，选取了某海域75 kHz走航式 ADCP 流速数据进行分析验证。

图1为包含流速异常值的走航式ADCP流速数据，共1 000个采样点，流速值单位为mm/s。由图1可见，该段流速数据的正常值应集中在1 000 mm/s以下，而在0～100、340～350和750～950三个区间的采样点都包含有明显高于正常流速的流速异常值，最大异常值可达12 000 mm/s以上，需要剔除。

图1 原始流速数据Fig.1 Raw velocity data

图2为使用传统Vondrak滤波方法和几种经典滤波方法得到的滤波结果。经过多次实验，对几种经典滤波方法的滤波参数进行设置。其中小波分析方法选取db4小波函数，对原始流速数据进行5层分解后再进行阈值重构，从而完成滤波处理；FFT方法是根据数据的频谱分析结果，选择30 Hz作为截止频率，得到了图2c所示的滤波结果；滑动平均方法采用了N=50的滑动窗对原始流速数据进行了时域滤波。

图2 4种滤波方法的滤波结果Fig.2 Filtering results of four filtering methods

从图2可见，传统Vondrak滤波方法和小波分析方法的滤波结果较为接近；FFT方法与前两种方法相比，保留了更多的波动信息；而滑动平均方法得到的滤波曲线在[0, 100]区间效果不理想，且存在大量明显的流速异常值，因此滑动平均方法并不适用于直接处理含有大量流速异常值的ADCP数据。对于仅存在测量噪声的数据段，传统Vondrak滤波、小波分析和FFT 3种滤波方法均具有较合理的滤波效果，滤波后的流速曲线光滑。但在具有流速异常值的数据段，滤波效果受到流速异常值的影响，滤波曲线朝着异常值方向“凸起”，特别是在[850, 900]区间采样点的流速异常值较为连续，“凸起”现象尤为明显。原始流速观测数据在[340, 350]区间的流速异常值最大，但异常值的个数很少且不是成片连续存在，虽然经过滤波后也存在“凸起”现象，但与[850, 900]区间的“凸起”相比较小，这说明传统Vondrak滤波、小波分析和FFT 3种滤波方法对连续的粗差更为敏感。

图3为传统Vondrak滤波、小波分析和FFT 3种滤波方法滤波结果对比图，3种滤波方法均可用于处理ADCP数据。从图中可以看出小波分析方法的滤波结果曲线在流速异常值附近“凸起”最大，FFT方法的滤波结果曲线虽然在流速异常值附近“凸起”最小，但是在含有其它测量噪声的数据段如[600, 800]采样点区间，滤波效果不是很理想，结果中含有更多的高频信息。表1为在[600, 800]采样点区间，3种方法滤波结果的相关系数与均方根误差对比，从表1中可知，传统Vondrak滤波方法在有效剔除测量噪声的同时，与原始信号具有最大相关性。

图3 3种滤波方法滤波结果对比图Fig.3 Comparison of filtering results of three filtering methods

表1 3种滤波方法滤波结果的相关系数与均方根误差对比Tab.1 Comparison of correlation coefficients and root mean square errors of filtering results obtained by three filtering methods

另外，小波分析方法的滤波效果与小波基函数、分解层数和截止频率等参数设置有关；FFT方法在使用之前也需要先对数据进行频谱分析，设置截止频率。而传统Vondrak滤波可以在对数据没有先验知识的情况下进行平滑滤波，较另外两种滤波方法具有一定优势。但对于具有流速异常值的数据区间，传统Vondrak滤波方法并不能进行很好的处理，滤波结果与实际流速值仍有较大差异，因此引入抗差Vondrak滤波方法来对ADCP数据进行更加合理的平滑滤波。

2.2 抗差Vondrak滤波方法的处理结果

抗差Vondrak滤波方法首先采用IGG3方法剔除流速异常值。图4为IGG3方法剔除流速异常值后的数据与原始数据对比图。从图中可以看出IGG3方法可以有效剔除流速误差，包含在[0, 100]、[340, 350]和[750, 950] 3个采样点区间的流速误差基本被剔除。[0, 100]和[340, 350] 2个采样点区间的流速异常值剔除效果较好，已十分接近正常流速值。[750, 950]采样点区间的流速异常值经过计算后减小很多，但还高于正常流速值。图5为传统Vondrak滤波结果与抗差Vondrak滤波结果对比图，从图中可以看出抗差Vondrak滤波结果要明显优于传统Vondrak滤波结果。抗差Vondrak滤波方法对 “凸起”现象有明显的抑制效果，在[340, 350]采样点区间，“凸起”现象已经消失，在[750, 950]采样点区间的“凸起”现象也被很好地抑制，流速已在正常值附近。

图4 IGG3方法剔除流速异常值后数据与原始数据对比图Fig.4 Comparison of the data after eliminating abnormal velocity values by using IGG3 method with original data

图5 传统Vondrak滤波结果与抗差Vondrak滤波结果对比图Fig.5 Comparison of traditional Vondrak filtering results with Robust Vondrak filtering results

3 内孤立波提取结果与分析

经过对ADCP数据的处理与实验分析，验证了抗差Vondrak滤波方法的有效性，现选一段某海域包含有内孤立波成分的走航式ADCP数据进行处理与分析，验证抗差Vondrak滤波方法在提取内孤立波中的适用性。该段数据为船载75 kHz ADCP采集获得，盲区深度42 m，垂直分辨率为16 m，时间分辨率为2 s，共15层流速数据。图6为1～5层流速数据，可以看出各层流速数据中都包含有流速异常值和大量噪声。图7为该段原始数据的流速效果图，左侧部分有一个形状类似于内孤立波的成分，但各层的水平流速分布、垂直剪切及其变化均无法获得，因此也无法对该类似内孤立波成分进行分析。图中右下部分有大量斑点，可能由流速异常值导致，需要进一步分析。

图6 1～5层原始流速数据Fig.6 Primary velocity data of first to 5th layers

图7 ADCP原始数据效果图Fig.7 Effect diagram of ADCP raw data

图8为抗差Vondrak滤波结果的水平流速等值线图，从中可以观察到一个清晰的内孤立波成分，位于[100, 300]采样点区间，图中右下部分没有出现大量流速异常值，取而代之的是均匀分布的水平背景流场。从图8中可以看出，波高由海面下58 m水深至波峰154 m水深附近，幅度可达96 m以上。经计算可知该内孤立波最大水平流速是1 890 mm/s,其核心位于水深100 m附近。 除此之外，图8中还可以清晰看到相邻层水平速度的垂向变化，水平流速从内孤立波的核心由内到外依次下降。抗差Vondrak滤波方法很好地从含有大量噪声和流速异常值的ADCP原始数据中提取出了内孤立波成分，通过观测点的时间周期及水平速度垂向分层变化清晰，为后续内孤立波的形成机制与传播等内容的研究奠定了基础。

图8 抗差Vondrak滤波结果等值线图Fig.8 Contour map of Robust Vondrak filtering result

4 结论

在移动观测平台上的矢量观测一直是海洋观测领域的研究前沿和难点。走航式ADCP观测的海流数据是研究海洋内波过程的重要而珍贵的数据来源。然而，该观测数据中普遍含有大量的噪声和流速异常值，严重影响了对海洋内波过程的认知和预测。本文通过理论分析与海洋观测数据处理，证明了传统Vondark滤波方法对处理走航式ADCP数据的有效性。为了剔除流速异常值，在传统Vondark滤波方法的基础上引入IGG3权函数因子，设计了一种抗差Vondark滤波方法，并采用该方法成功提取出了内孤立波。研究结果表明，流速异常值会对传统Vondark滤波方法造成影响，滤波曲线会向流速异常值方向“凸起”，其中连续分布的流速异常值影响最大，即“凸起”现象最为明显。抗差Vondark滤波方法很好地抑制了这种“凸起”现象，在获得了平滑的滤波曲线基础上实现了高精度的数据拟合，提取出的内孤立波成分较准确且各层水平流速清晰，是一种有效的自适应内孤立波提取方法。