基于參數(shù)分析法的GNSS測波浮標(biāo)波浪方向譜估計(jì)

谷漢斌，王成成，單 瑞

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

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

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

波浪譜分析方法很多，它們各具優(yōu)缺點(diǎn)。這里只討論基于參數(shù)法的波浪方向譜生成，以及支持該方法的傅立葉變換在頻域里的譜驗(yàn)證。該方法便于應(yīng)用于高精度GNSS浮標(biāo)測量波浪的數(shù)據(jù)分析。

1 GNSS測波浮標(biāo)簡介與處理程序包

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

圖1 GNSS測波浮標(biāo)

基于Matlab程序設(shè)計(jì)了GNSS浮標(biāo)的數(shù)據(jù)處理程序包，包括高頻和低頻數(shù)字濾波、波浪統(tǒng)計(jì)分析、波浪頻譜分析和波浪方向譜分析，以及相關(guān)列表和圖形顯示功能。本文主要介紹波浪統(tǒng)計(jì)分析和波浪頻譜與方向譜分析方法和結(jié)果。

2 波浪統(tǒng)計(jì)分析

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

圖2 波浪出現(xiàn)方向和頻次

(1)

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

表1 數(shù)據(jù)前處理結(jié)果對比

表2 統(tǒng)計(jì)波高和周期對比

3 頻譜生成方法與驗(yàn)證

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

(2)

其逆變換為

(3)

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

圖3 JONSWAP計(jì)算譜與估計(jì)譜對比

圖4 波面垂向位移時(shí)間過程

4 方向譜估計(jì)方法

前文提到波浪方向譜有多種估計(jì)方法，其中參數(shù)估計(jì)法是一種簡單易行的方法，該方法由LONGUET-HIGGINS等人提出[12]，BORGMAN[13]、PANICKER和BORGMAN[14]及HASSELMANN等人[15]對該方法進(jìn)行了不斷完善，分別提出截?cái)嗟母盗⑷~系列、余弦函數(shù)冪次系列函數(shù)、圓形正態(tài)分布和正態(tài)分布等截?cái)嗪瘮?shù)表示方向分布。其中截?cái)嗟母盗⑷~系列形式最簡單

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

(4)

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

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

圖5 某次浮標(biāo)升沉和水平方向時(shí)間系列

圖6 波浪頻譜

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

(5)

或表示為

(6)

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

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

(7)

(8)

于是，可計(jì)算

(9)

(10)

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

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

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

表3 統(tǒng)計(jì)波高和周期對比

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

5 在GNSS測波浮標(biāo)中的使用

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

圖9 高精度定位浮標(biāo)位移時(shí)間過程

圖10 數(shù)據(jù)分析結(jié)果

6 結(jié)論

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