基于參數分析法的GNSS測波浮標波浪方向譜估計

谷漢斌，王成成，單 瑞

(1.寧波大學 土木工程與地理環境學院/海洋工程研究院，寧波 315211；2.同濟大學 測繪與地理信息學院，上海 200092；3.中國地質調查局青島海洋地質研究所，青島 266237)

隨著全球導航衛星系統(Global Navigation Satellite System, GNSS)技術的不斷發展和進步，尤其是實時精密單點定位技術(Precise Point Positioning，PPP)的不斷完善，給海洋波浪觀測提供了便捷高效的手段[1-2]。GNSS浮標測波避免了原有的以慣性測波浮標為主的測波技術對浮標的高度的技術要求，如“波浪騎士”，如此小的浮標，直徑0.4 m、0.7 m和0.9 m，內部設液體穩定平臺，重力擺保持垂向加速度計始終處于垂立姿態；電子陀螺儀保持北向，其他兩個水平加速度計、兩個縱搖和橫搖傳感器、電池等密閉封存在這樣的小空間內，給高精度浮標測波帶來很高的技術要求。而高精度GNSS測波浮標可以通過PPP技術確定浮標的升沉和水平方向位移，動態測量精度達到厘米級，給海洋波浪測量提供了一種新的手段。其比GPS測波浮標采用的定位衛星多，抗干擾更強；當然，其算法和裝置也有優化的空間。

海洋波浪分析包括統計學分析和譜分析方法[3]。最傳統的波浪數據分析是統計分析法，針對水面升沉運動給出波浪統計參數，以便于工程應用。隨著對海洋波浪觀測和分析研究的深入以及數學上隨機過程的理論發展，譜分析方法越來越受到重視[3-5]。波浪譜分析方法包括頻譜分析和方向譜分析，頻譜分析只在頻域上進行，不考慮空間的分布特性，頻譜分析包括相關函數法、傅立葉變換法(FFT)和最大熵法(MEP)等，其中前兩種方法簡單易行；波浪的方向譜分析法不僅包括了波浪在頻域上的變化，還考慮方向上的分布，其中有直接傅立葉變換法(DFTM)、參數法、最大可能伸展法(EMLM)、最大熵法、貝葉斯法(BDM)、最大熵延申法(EMEP)等方法[5]，其中參數法是比較容易實現的方法，只需要升沉和水平方向三組位移過程就可以分析波浪的方向譜。當然，近些年波浪譜的研究也有很大發展，如采用小波變換分析海洋波浪方向譜[6]、采用聲學多普勒流速剖面儀(ADCP)估計波向[7]和海浪方向譜的數據測量和估計[8]及波浪特征分析[9]等。近年來，波浪方向譜的估計和觀測在國外也一直受到重視[10-11]。

波浪譜分析方法很多，它們各具優缺點。這里只討論基于參數法的波浪方向譜生成，以及支持該方法的傅立葉變換在頻域里的譜驗證。該方法便于應用于高精度GNSS浮標測量波浪的數據分析。

1 GNSS測波浮標簡介與處理程序包

GNSS浮標由1個主浮筒和3個輔浮筒組成，均為圓柱狀，3個輔浮筒彼此等間隔120°，主浮筒直徑約0.5 m，輔助浮筒直徑約0.35 m，其頂部為GNSS接收天線，主浮筒內安裝1臺TRIMBLE NETR9接收機,見圖1。GNSS浮標水上部分高約0.50 m，水下部分高約0.4 m，實驗時采樣頻率為5 Hz。GNSS單元可測得該浮標高精度位置及速度變化。本文擬為該類浮標測波提供數據分析方法。

圖1 GNSS測波浮標

基于Matlab程序設計了GNSS浮標的數據處理程序包，包括高頻和低頻數字濾波、波浪統計分析、波浪頻譜分析和波浪方向譜分析，以及相關列表和圖形顯示功能。本文主要介紹波浪統計分析和波浪頻譜與方向譜分析方法和結果。

2 波浪統計分析

波浪統計分析是最常用的波浪觀測數據處理方法，通常按照工程應用需求，統計出平均波高、三分之一大波的平均波高即有效波高、十分之一大波的平均波高、最大波高等，以及對應的波浪周期。這里只討論針對一個波浪系列的統計分析，波浪統計分析還包括如波浪重現期等相關的分析，不是這里討論的重點。對于實測的波面位移時間變化 ，i=1,2……n，計算其平均值，該值為平均水平面，采用這個平均值對波面位移時間過程在位移時間坐標系內切割位移時間過程線，把位移零點平移到這個均值線，得到一系列跨零點，按上跨零點法分析波高和周期，以從大到小的順序排列，找出前述各特征波高。

圖2 波浪出現方向和頻次

(1)

將波浪騎士浮標原始數據采用式(1)計算各特征值，與波浪騎士處理軟件W@ves21的分析結果進行對比，只有標準差比較接近，其他數據存在一些偏差(表1)，可能波浪騎士處理時，采用了不同的濾波或分段方式。對波浪騎士測得的升沉和兩個水平方向位移進行相關分析，相關系數計算值與其處理結果一致，升沉與南北向位移相關系數0.13，與東西向位移為0.004，南北向位移與東西向位移相關系數0.76。采用升沉數據進行波高統計分析的對比結果見表2，波高和周期的統計值相差基本接近。

表1 數據前處理結果對比

表2 統計波高和周期對比

3 頻譜生成方法與驗證

假定波面垂直向上位移時間過程η(t)是一個平穩隨機過程，傅立葉變換將其變換到頻域X(ω)，即

(2)

其逆變換為

(3)

這樣頻域和時域是一個相互變換的過程。因此，波浪頻譜已知，可以通過譜分割的辦法構造時間域的隨機過程[3]。由于隨機相位的引入，一個波譜可以生成無數個隨機過程，但是這些隨機過程對應一個頻譜。目前常用的波譜有JONSWAP譜、我國港口工程規范譜、PM譜等。給出波譜參數，波譜就能夠表達出來；例如給定有效波高1.0 m、有效波浪周期8.5 s、水深50 m，采用JONSWAP譜，譜型見圖3星點所示。對該譜進行分割200份，最小截斷頻率fp/5，最大取5fp，采用Longuet-Higgins不規則波模型構建出一個隨機波面過程見圖4。對這個隨機波面過程進行逆變換分析，即為譜估計，采用快速傅立葉變換進行譜估計[3]，得到估計譜如圖3黑線所示?？焖俑盗⑷~變換過程中，要注意(1)變換前使用截尾窗將整個隨機波面過程看作一個長周期進行周期化，變換后去掉截尾窗影響；(2)快速傅立葉變換后的粗譜需要進行平滑處理；(3)傅立葉變換時，針對時間步長△t=1的時間系列，變換后的頻域幅值要換算成對應實際△t的值；(4)由于很多計算機語言提供了傅立葉變換函數庫，針對一般性問題，頻域函數為復數，做波浪頻譜分析時常采用實部，輸出頻譜時對頻率輸出點數進行了截半處理。

圖3 JONSWAP計算譜與估計譜對比

圖4 波面垂向位移時間過程

4 方向譜估計方法

前文提到波浪方向譜有多種估計方法，其中參數估計法是一種簡單易行的方法，該方法由LONGUET-HIGGINS等人提出[12]，BORGMAN[13]、PANICKER和BORGMAN[14]及HASSELMANN等人[15]對該方法進行了不斷完善，分別提出截斷的傅立葉系列、余弦函數冪次系列函數、圓形正態分布和正態分布等截斷函數表示方向分布。其中截斷的傅立葉系列形式最簡單

S(f，θ)=S(f)[a0(f)+Σn{an(f)cos(nθ)+bn(f)sin(nθ)}]

(4)

式中：S(f，θ)和S(f)分別為波浪方向譜和頻譜，后面方括號內的項為方向分布函數G(f，θ)，其中a0、an和bn為系數，f為頻率。該式適用于浮標觀測方向譜，并取得不少測波資料。理解其處理方法，更有助于理解觀測數據。

在波浪騎士用戶手冊中，頻譜采用升沉時間系列(圖5)的傅立葉變換得到。采用本文驗證的程序對該系列進行分析，結果見圖6。與波浪騎士的處理結果基本一致，波浪騎士數據振蕩更大一些，Smax=1.18，主頻略小，本文計算Smax=1.2；低頻部分也基本一致。這里有兩點需要注意：(1)波浪騎士手冊提到對原始數據采用兩種濾波，周期分別為43 s和256 s；本文對數據系列采用濾波周期47 s。(2)頻譜平滑處理時，平滑點數P[3]本文取N/160，N為傅立葉變換系列處理點數。

圖5 某次浮標升沉和水平方向時間系列

圖6 波浪頻譜

波浪騎士用戶手冊中方向分布函數表示為

(5)

或表示為

(6)

其中系數從浮標內部加速度傳感器二重積分得到的升沉、水平方向位移時間過程序列(圖5)的傅立葉變換的頻域系數計算得到。對圖5中3個位移時間系列進行傅立葉變換，從升沉、南北向和西東向位移時間系列分別得到頻域傅立葉系數的實部與虛部αvf、αnf、αwf和βvf、βnf、βwf，傅立葉系數表示為

Avf=αvf+iβvf，Anf=αnf+iβnf，Awf=αwf+iβwf

(7)

(8)

于是，可計算

(9)

(10)

本文采用上述方法對圖5的時間系列進行方向函數估計，結果見圖7-a，采用式(5)或式(6)計算有一定差別，式(5)形成的方向函數在圖中左側，式(6)形成的在右側。它們的共同點是方向函數的每個頻率在各方向積分為1，但是都存在方向函數為負的方向，與物理意義不符。盡管伯格曼(Borgman)提出的平滑函數可以修正方向函數，使之非負，但降低了峰頻，擴展了方向分布范圍[3]。本文建議，對方向函數修正的具體做法可采用如下步驟：(1)求方向函數為負的積分面積，方向函數為負值時調整為零；(2)把負數的積分面積按面積權數分配到方向函數為正的部分，修正正的方向函數；(3)檢驗修正后的方向函數的各方向積分為1。按這一修正方法，修正后方向函數分布見圖7-b。

7-a 由式(5)和式(6)分別估計的方向函數 7-b 式(5)和式(6)分別估計方向函數的修正結果

基于波浪的頻譜和方向函數分布，給出方向譜S(f,θ)，如圖8所示，其中圖8-a方向函數分布按式(5)計算，圖8-b按式(6)計算；基于兩個公式的計算相差不多。需要注意，按式(6)計算前需要計算波浪方向角θ0，該角度為方位角，而在浮標上角度是按北向a1和西向b1為正計算的角度，該角度要換算為方位角，否則方向譜譜峰位置會出錯。計算結果的主要參數與波浪騎士結果對比見表3，本文計算值都偏大，應該和濾波有關。

表3 統計波高和周期對比

8-a 按式(5)計算的方向譜 8-b 按式(6)計算的方向譜

5 在GNSS測波浮標中的使用

將上述波浪方向譜估計方法應用于GNSS測波浮標，用于波浪觀測，實測數據包括浮標升沉和兩個水平方向N和E的水平位移，在浮標受海流和風的影響比較小的情況下，該浮標可用于波浪數據觀測。圖9給出浮標受波浪作用產生的升沉及兩個水平位移時間過程。根據位移時間過程采用分析程序進行分析計算，結果見圖10。從實測波浪方向頻率圖可以看出，波浪主要出現在東北和西南方向；峰頻接近0.13 Hz，最大譜密度接近0.43；方向分布函數與實際波浪出現方向一致，同時給出了方向譜。統計H1/3=0.81 m，Hm0=0.98 m比統計值略大，統計平均波浪周期3.69 s，譜計算平均波浪周期4.00 s，也比統計值略大。總體看方向譜計算方法得到的結果基本可信。

圖9 高精度定位浮標位移時間過程

圖10 數據分析結果

6 結論

本文針對GNSS浮標應用于測波的波浪數據處理方法，參考波浪騎士手冊對實測結果進行了處理方法研究。給出了波浪統計處理方法、頻譜處理方法和方向譜處理方法，對數據處理方法進行了編程，形成一套適用于高精度GNSS測波浮標的程序包。在波浪方向譜分析中，方向函數易出現負值，與物理定義不符，提出了修正方法。數據處理結果基本與波浪騎士處理結果一致，方向譜計算結果也與統計分析結果基本一致。說明該程序包可應用于高精度GNSS浮標對波浪進行海上觀測。