抗差Vondrak濾波方法在內(nèi)孤立波提取中的應(yīng)用研究

張佳麗，張安民，孫朝輝，張學峰，張 亮*

(1.天津大學 海洋科學與技術(shù)學院，天津 300072；2.自然資源部 第二海洋研究所, 浙江 杭州 310012)

0 引言

內(nèi)波指發(fā)生在密度穩(wěn)定層結(jié)流體中、最大振幅出現(xiàn)在流體內(nèi)部的波動，是一種普遍存在的自然現(xiàn)象，這種自然現(xiàn)象存在于海洋內(nèi)部而不易被觀測到[1-2]。內(nèi)孤立波是一種特殊的內(nèi)波，具有振幅大、傳播速度快和非線性的特點[3]，能夠進行單向長距離傳播。內(nèi)孤立波經(jīng)過時，水體的垂向流速、水平流速和垂向剪切均很強，對海上石油鉆機和水下航行器安全運行造成了威脅。因此，開展內(nèi)孤立波的研究對海洋科學研究、海洋資源開發(fā)和國防安全保障均具有重要意義[4]。GUO et al[5]綜述了近10 a來國內(nèi)外對內(nèi)孤立波及其相關(guān)動力學過程的研究狀況，對遙感圖像、現(xiàn)場測量、數(shù)值模擬 3種研究方法進行了總結(jié)歸納并分析了中國南海的內(nèi)孤立波，其中現(xiàn)場測量的方法大多使用了ADCP( Acoustic Doppler Current Profiler，聲學多普勒流速剖面儀)。ADCP是對內(nèi)波進行聲學監(jiān)測的有效手段，利用它測得流速、水回波信號強度等數(shù)據(jù)，可以反演出內(nèi)波信息[6]。HUANG et al[7]使用ADCP和CTD(Conductivity-Temperature-Depth, 溫鹽深儀)，通過現(xiàn)場測量捕捉到了南海的一個極端內(nèi)孤立波，并進行了成因分析。XIE et al[8]使用ADCP對河口重力流產(chǎn)生的內(nèi)孤立波破碎過程進行了研究。然而，在實際使用中，ADCP的測量精度在很大程度上取決于海洋環(huán)境條件及船載輔助設(shè)備記錄精度，特別是走航式ADCP受到噪聲的影響，獲得的數(shù)據(jù)十分復雜[9]，測得流速數(shù)據(jù)中包含有各種測量噪聲和一些流速異常值。移動中值、三階平滑低通濾波、窗函數(shù)及幾種濾波方法組合的綜合濾波方法是目前已有的濾波方法[10]，但是對于混有流速異常值及各種噪聲誤差的ADCP數(shù)據(jù)，這幾種方法均不能很好地實現(xiàn)自適應(yīng)濾波處理，存在截止頻率設(shè)置不合理導致濾波失真、流速異常值剔除效果不理想等問題。在一些ADCP處理軟件中，流速異常值需要人肉眼觀察判斷，根據(jù)測區(qū)的實際情況，利用計算機設(shè)置相關(guān)檢測閾值來識別和剔除[11]。

1 抗差Vondrak濾波方法研究

1.1 傳統(tǒng)Vondrak濾波方法

Vondrak濾波方法需要綜合考慮擬合度與平滑度兩方面，其中擬合度反映觀測信息，而平滑度反映非觀測信息。擬合度與平滑度公式，用最小二乘法求解。對于一組有N個觀測值的數(shù)值序列x(ti)，ti為測量時刻(i=1,2,...,N)，Vondrak濾波需要滿足準則Q：

(1)

其中：Q→min表示Vondrak濾波準則達到最小；F表示濾波結(jié)果與原始觀測數(shù)據(jù)的擬合度；S表示濾波結(jié)果的平滑度；ε∈(0,+∞)，是一個無量綱的正數(shù)，定義為平滑因子。

Vondrak濾波的本質(zhì)是通過確定合理的平滑因子，達到最佳濾波效果。當ε→∞時，必須有F→0，即此時濾波結(jié)果與原始觀測數(shù)據(jù)幾乎重合，而當F=0時，濾波結(jié)果與原始觀測數(shù)據(jù)完全重合，失去濾波功能; 當ε→0時，必須有S→0,即此時平滑度趨向于0，濾波結(jié)果近似為拋物線，導致濾波過度。HVF法是通過反復迭代計算擬合函數(shù)與平滑函數(shù)的方差，用兩者比值反復修復權(quán)系數(shù)來確定平滑因子,其基本原理如下：

(1)建立兩類誤差觀測方程，其中V1為擬合度誤差，V2為平滑度誤差：

(2)

(2)令兩類觀測誤差的權(quán)值分別為P1和P2，根據(jù)最小二乘原理得到：

(3)

將式(2)代入式(3)得到誤差方程：

(4)

式(3)的法方程為：

(5)

式中：

(6)

(7)

其中：B1為n階單位陣,f1=y,f2=0,B2為一帶型矩陣[13]。

(8)

(9)

(10)

(11)

其中：C為任意常數(shù)。

(12)

其中：P為初始權(quán)值，這里取P=1。

1.2 抗差Vondrak濾波方法

在實際ADCP觀測數(shù)據(jù)中，粗差的存在會嚴重影響數(shù)據(jù)處理結(jié)果的精度，在ADCP測量粗差點，即流速異常值附近的數(shù)據(jù)擬合結(jié)果是不可靠的，因此，需要首先對粗差進行處理。抗差Vondrak濾波方法分為兩步，首先基于IGG3方法去除原始觀測數(shù)據(jù)中的流速異常值，然后采用HVF方法進行時域濾波得到平滑的流速曲線。

IGG3算法是一種加權(quán)迭代的方法，根據(jù)殘差和有關(guān)參數(shù)，按照所選的權(quán)函數(shù)計算每個觀測值在下一次迭代中的權(quán)。經(jīng)過計算，含有粗差的觀測值的權(quán)會越來越小或降為0。IGG3的算法模型為：

(13)

(14)

其中：V為觀測值殘差;σ為方差因子，這里用均方根誤差表示為：

(15)

(16)

其中：P為初始權(quán)值，這里取P=1。

當海流觀測中不含流速異常值時，僅采用傳統(tǒng)的Vondrak濾波方法就可以得到很好的平滑曲線，船載ADCP在測量過程中出現(xiàn)的噪聲可以被有效濾除。當流速異常值存在時，僅通過平滑因子來控制濾波性能的傳統(tǒng)Vondrak濾波方法不能對流速異常值進行有效的剔除，在流速異常值附近，數(shù)據(jù)擬合結(jié)果會出現(xiàn)偏差，這對之后內(nèi)孤立波的提取與分析造成不良影響。為了改善濾波效果，引入IGG3方法對原始數(shù)據(jù)進行流速異常值自適應(yīng)剔除，基本操作步驟如下：

(2)將步驟(1)中計算得到的值帶入IGG3方法的權(quán)因子函數(shù)中，得到各觀測值的新權(quán)；

(3)通過新權(quán)與原始觀測數(shù)據(jù)相乘得到剔除流速異常值的觀測數(shù)據(jù)。

1.3 經(jīng)典濾波方法

FFT、小波分析和滑動平均是3種經(jīng)典的濾波降噪方法。FFT[20]是離散傅里葉變換(Discrete Fourier Transform, DFT)的一種快速算法，N點序列x(N)的N點DFT可以表示為：

(17)

FFT將原始信號變換到頻域中進行分析，根據(jù)原始信號的頻譜，在頻域空間對噪聲成分進行有效抑制，再將結(jié)果進行傅立葉逆變換，得到降噪后的時域信號。FFT方法對周期性明顯的信號降噪效果較好，而對于非平穩(wěn)、非周期信號，特別是有異常值信號的降噪效果較差。此外，F(xiàn)FT方法需要利用信號的全部時域信息進行整體變換，缺乏時域分析功能，且信號不能同時在時域和頻域準確定位。

小波分析[21-22]能有效地從原始信號中提取有效信息，并通過伸縮和平移對信號在時域和頻域進行多尺度分析。降噪是小波分析在信號處理領(lǐng)域的主要應(yīng)用之一，該方法的降噪步驟如下：(1)對信號進行N層小波分解，根據(jù)小波分解的特點，將信號分為低頻部分和高頻部分；(2)選擇1個閾值，對高頻系數(shù)進行抑制；(3)將處理過的系數(shù)進行重構(gòu)，得到降噪后的信號。小波變換在時域和頻域均具有很好的局部化特征，但是在小波變換方法降噪的過程中，小波基函數(shù)、分解層數(shù)以及閾值的選擇都很關(guān)鍵，直接影響到對真實信號的降噪效果。

滑動平均是一種簡單的數(shù)據(jù)平滑方法，通過取N個相鄰數(shù)據(jù)的平均值來表示所取數(shù)據(jù)的端點或中點值。對于N個含有噪聲的實測數(shù)據(jù){Zi}，取n個相鄰數(shù)據(jù)的平均值，即滑動窗大小為n，滑動平均的數(shù)學表達式為：

(18)

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

滑動平均方法計算簡單，處理數(shù)據(jù)速度快，但是需要對滑動窗口進行選擇。

2 對ADCP數(shù)據(jù)的處理

2.1 傳統(tǒng)Vondrak濾波方法與經(jīng)典濾波方法對比

為了檢驗Vondrak濾波方法在海洋內(nèi)孤立波提取中的適用性，選取了某海域75 kHz走航式 ADCP 流速數(shù)據(jù)進行分析驗證。

圖1為包含流速異常值的走航式ADCP流速數(shù)據(jù)，共1 000個采樣點，流速值單位為mm/s。由圖1可見，該段流速數(shù)據(jù)的正常值應(yīng)集中在1 000 mm/s以下，而在0～100、340～350和750～950三個區(qū)間的采樣點都包含有明顯高于正常流速的流速異常值，最大異常值可達12 000 mm/s以上，需要剔除。

圖1 原始流速數(shù)據(jù)Fig.1 Raw velocity data

圖2為使用傳統(tǒng)Vondrak濾波方法和幾種經(jīng)典濾波方法得到的濾波結(jié)果。經(jīng)過多次實驗，對幾種經(jīng)典濾波方法的濾波參數(shù)進行設(shè)置。其中小波分析方法選取db4小波函數(shù)，對原始流速數(shù)據(jù)進行5層分解后再進行閾值重構(gòu)，從而完成濾波處理；FFT方法是根據(jù)數(shù)據(jù)的頻譜分析結(jié)果，選擇30 Hz作為截止頻率，得到了圖2c所示的濾波結(jié)果；滑動平均方法采用了N=50的滑動窗對原始流速數(shù)據(jù)進行了時域濾波。

圖2 4種濾波方法的濾波結(jié)果Fig.2 Filtering results of four filtering methods

從圖2可見，傳統(tǒng)Vondrak濾波方法和小波分析方法的濾波結(jié)果較為接近；FFT方法與前兩種方法相比，保留了更多的波動信息；而滑動平均方法得到的濾波曲線在[0, 100]區(qū)間效果不理想，且存在大量明顯的流速異常值，因此滑動平均方法并不適用于直接處理含有大量流速異常值的ADCP數(shù)據(jù)。對于僅存在測量噪聲的數(shù)據(jù)段，傳統(tǒng)Vondrak濾波、小波分析和FFT 3種濾波方法均具有較合理的濾波效果，濾波后的流速曲線光滑。但在具有流速異常值的數(shù)據(jù)段，濾波效果受到流速異常值的影響，濾波曲線朝著異常值方向“凸起”，特別是在[850, 900]區(qū)間采樣點的流速異常值較為連續(xù)，“凸起”現(xiàn)象尤為明顯。原始流速觀測數(shù)據(jù)在[340, 350]區(qū)間的流速異常值最大，但異常值的個數(shù)很少且不是成片連續(xù)存在，雖然經(jīng)過濾波后也存在“凸起”現(xiàn)象，但與[850, 900]區(qū)間的“凸起”相比較小，這說明傳統(tǒng)Vondrak濾波、小波分析和FFT 3種濾波方法對連續(xù)的粗差更為敏感。

圖3為傳統(tǒng)Vondrak濾波、小波分析和FFT 3種濾波方法濾波結(jié)果對比圖，3種濾波方法均可用于處理ADCP數(shù)據(jù)。從圖中可以看出小波分析方法的濾波結(jié)果曲線在流速異常值附近“凸起”最大，F(xiàn)FT方法的濾波結(jié)果曲線雖然在流速異常值附近“凸起”最小，但是在含有其它測量噪聲的數(shù)據(jù)段如[600, 800]采樣點區(qū)間，濾波效果不是很理想，結(jié)果中含有更多的高頻信息。表1為在[600, 800]采樣點區(qū)間，3種方法濾波結(jié)果的相關(guān)系數(shù)與均方根誤差對比，從表1中可知，傳統(tǒng)Vondrak濾波方法在有效剔除測量噪聲的同時，與原始信號具有最大相關(guān)性。

圖3 3種濾波方法濾波結(jié)果對比圖Fig.3 Comparison of filtering results of three filtering methods

表1 3種濾波方法濾波結(jié)果的相關(guān)系數(shù)與均方根誤差對比Tab.1 Comparison of correlation coefficients and root mean square errors of filtering results obtained by three filtering methods

另外，小波分析方法的濾波效果與小波基函數(shù)、分解層數(shù)和截止頻率等參數(shù)設(shè)置有關(guān)；FFT方法在使用之前也需要先對數(shù)據(jù)進行頻譜分析，設(shè)置截止頻率。而傳統(tǒng)Vondrak濾波可以在對數(shù)據(jù)沒有先驗知識的情況下進行平滑濾波，較另外兩種濾波方法具有一定優(yōu)勢。但對于具有流速異常值的數(shù)據(jù)區(qū)間，傳統(tǒng)Vondrak濾波方法并不能進行很好的處理，濾波結(jié)果與實際流速值仍有較大差異，因此引入抗差Vondrak濾波方法來對ADCP數(shù)據(jù)進行更加合理的平滑濾波。

2.2 抗差Vondrak濾波方法的處理結(jié)果

抗差Vondrak濾波方法首先采用IGG3方法剔除流速異常值。圖4為IGG3方法剔除流速異常值后的數(shù)據(jù)與原始數(shù)據(jù)對比圖。從圖中可以看出IGG3方法可以有效剔除流速誤差，包含在[0, 100]、[340, 350]和[750, 950] 3個采樣點區(qū)間的流速誤差基本被剔除。[0, 100]和[340, 350] 2個采樣點區(qū)間的流速異常值剔除效果較好，已十分接近正常流速值。[750, 950]采樣點區(qū)間的流速異常值經(jīng)過計算后減小很多，但還高于正常流速值。圖5為傳統(tǒng)Vondrak濾波結(jié)果與抗差Vondrak濾波結(jié)果對比圖，從圖中可以看出抗差Vondrak濾波結(jié)果要明顯優(yōu)于傳統(tǒng)Vondrak濾波結(jié)果?？共頥ondrak濾波方法對 “凸起”現(xiàn)象有明顯的抑制效果，在[340, 350]采樣點區(qū)間，“凸起”現(xiàn)象已經(jīng)消失，在[750, 950]采樣點區(qū)間的“凸起”現(xiàn)象也被很好地抑制，流速已在正常值附近。

圖4 IGG3方法剔除流速異常值后數(shù)據(jù)與原始數(shù)據(jù)對比圖Fig.4 Comparison of the data after eliminating abnormal velocity values by using IGG3 method with original data

圖5 傳統(tǒng)Vondrak濾波結(jié)果與抗差Vondrak濾波結(jié)果對比圖Fig.5 Comparison of traditional Vondrak filtering results with Robust Vondrak filtering results

3 內(nèi)孤立波提取結(jié)果與分析

經(jīng)過對ADCP數(shù)據(jù)的處理與實驗分析，驗證了抗差Vondrak濾波方法的有效性，現(xiàn)選一段某海域包含有內(nèi)孤立波成分的走航式ADCP數(shù)據(jù)進行處理與分析，驗證抗差Vondrak濾波方法在提取內(nèi)孤立波中的適用性。該段數(shù)據(jù)為船載75 kHz ADCP采集獲得，盲區(qū)深度42 m，垂直分辨率為16 m，時間分辨率為2 s，共15層流速數(shù)據(jù)。圖6為1～5層流速數(shù)據(jù)，可以看出各層流速數(shù)據(jù)中都包含有流速異常值和大量噪聲。圖7為該段原始數(shù)據(jù)的流速效果圖，左側(cè)部分有一個形狀類似于內(nèi)孤立波的成分，但各層的水平流速分布、垂直剪切及其變化均無法獲得，因此也無法對該類似內(nèi)孤立波成分進行分析。圖中右下部分有大量斑點，可能由流速異常值導致，需要進一步分析。

圖6 1～5層原始流速數(shù)據(jù)Fig.6 Primary velocity data of first to 5th layers

圖7 ADCP原始數(shù)據(jù)效果圖Fig.7 Effect diagram of ADCP raw data

圖8為抗差Vondrak濾波結(jié)果的水平流速等值線圖，從中可以觀察到一個清晰的內(nèi)孤立波成分，位于[100, 300]采樣點區(qū)間，圖中右下部分沒有出現(xiàn)大量流速異常值，取而代之的是均勻分布的水平背景流場。從圖8中可以看出，波高由海面下58 m水深至波峰154 m水深附近，幅度可達96 m以上。經(jīng)計算可知該內(nèi)孤立波最大水平流速是1 890 mm/s,其核心位于水深100 m附近。 除此之外，圖8中還可以清晰看到相鄰層水平速度的垂向變化，水平流速從內(nèi)孤立波的核心由內(nèi)到外依次下降?？共頥ondrak濾波方法很好地從含有大量噪聲和流速異常值的ADCP原始數(shù)據(jù)中提取出了內(nèi)孤立波成分，通過觀測點的時間周期及水平速度垂向分層變化清晰，為后續(xù)內(nèi)孤立波的形成機制與傳播等內(nèi)容的研究奠定了基礎(chǔ)。

圖8 抗差Vondrak濾波結(jié)果等值線圖Fig.8 Contour map of Robust Vondrak filtering result

4 結(jié)論

在移動觀測平臺上的矢量觀測一直是海洋觀測領(lǐng)域的研究前沿和難點。走航式ADCP觀測的海流數(shù)據(jù)是研究海洋內(nèi)波過程的重要而珍貴的數(shù)據(jù)來源。然而，該觀測數(shù)據(jù)中普遍含有大量的噪聲和流速異常值，嚴重影響了對海洋內(nèi)波過程的認知和預測。本文通過理論分析與海洋觀測數(shù)據(jù)處理，證明了傳統(tǒng)Vondark濾波方法對處理走航式ADCP數(shù)據(jù)的有效性。為了剔除流速異常值，在傳統(tǒng)Vondark濾波方法的基礎(chǔ)上引入IGG3權(quán)函數(shù)因子，設(shè)計了一種抗差Vondark濾波方法，并采用該方法成功提取出了內(nèi)孤立波。研究結(jié)果表明，流速異常值會對傳統(tǒng)Vondark濾波方法造成影響，濾波曲線會向流速異常值方向“凸起”，其中連續(xù)分布的流速異常值影響最大，即“凸起”現(xiàn)象最為明顯?？共頥ondark濾波方法很好地抑制了這種“凸起”現(xiàn)象，在獲得了平滑的濾波曲線基礎(chǔ)上實現(xiàn)了高精度的數(shù)據(jù)擬合，提取出的內(nèi)孤立波成分較準確且各層水平流速清晰，是一種有效的自適應(yīng)內(nèi)孤立波提取方法。