ADCP方向譜估計(jì)的傾斜修正算法研究

鄭威，楊立，惠力，劉敏，王志

(山東省科學(xué)院海洋儀器儀表研究所，山東青島266001)

聲學(xué)多普勒流速剖面儀(ADCP)是用來(lái)測(cè)量不同深度上海水流速的儀器。美國(guó)RD公司“駿馬”系列ADCP上安放有姿態(tài)傳感器，可以測(cè)量得到俯仰角(pitch)和橫滾角(roll)，當(dāng)ADCP安放在船上的時(shí)候，姿態(tài)角時(shí)刻不斷變化，需要對(duì)ADCP的每組測(cè)量值進(jìn)行傾斜修正。而當(dāng)坐底式ADCP安放在水底的固定的支架上時(shí)，姿態(tài)角變化較小，不用對(duì)每組數(shù)據(jù)進(jìn)行傾斜修正，只有當(dāng)姿態(tài)角較大時(shí)才需要進(jìn)行傾斜修正。

RD公司技術(shù)文獻(xiàn)沒(méi)有給出傾斜狀態(tài)下方向譜的計(jì)算方法，本文通過(guò)旋轉(zhuǎn)矩陣的引入給出了傾斜狀態(tài)下方向譜估計(jì)方法。RD公司商用軟件WAVESMON選取當(dāng)傾斜姿態(tài)角大于10°時(shí)進(jìn)行修正［1－4］，而當(dāng)使用坐底式ADCP時(shí)，實(shí)測(cè)數(shù)據(jù)姿態(tài)角大于10°的情況不多，小于10°的姿態(tài)角同樣會(huì)帶來(lái)較大誤差，所以為了得到盡可能準(zhǔn)確的估計(jì)值，通常是對(duì)每組數(shù)據(jù)進(jìn)行傾斜修正［5］。我們研究發(fā)現(xiàn)當(dāng)ADCP處于傾斜條件下，方向譜估計(jì)誤差不僅需要考慮姿態(tài)角，同時(shí)還必須要考慮到ADCP的安放深度。本文重點(diǎn)研究不同傾斜狀態(tài)下，ADCP傾斜姿態(tài)角和安放深度對(duì)于方向譜估計(jì)的影響，并以RD公司JANUS配置為例進(jìn)行說(shuō)明。

1 方向譜估計(jì)基本原理

Isobe在1984年利用方向譜和互譜間的關(guān)系來(lái)進(jìn)行方向譜估計(jì)，這也是目前為止最常用的方向譜估計(jì)基本理論［1］。假定波浪振幅很小且由多數(shù)組成波迭加而成，可推導(dǎo)出任意兩個(gè)波浪特性間的互譜等于相應(yīng)波浪特性與波面間的傳遞函數(shù)的乘積的傅里葉變換，即

f代表海浪頻率，k代表波數(shù)，φmn(f)代表第m個(gè)測(cè)點(diǎn)上和第n個(gè)測(cè)點(diǎn)波浪特性之間的互譜，Hm是波浪特性m的傳遞函數(shù)，S(f，θ)代表海浪方向譜，xmn，ymn分別代表m和n測(cè)點(diǎn)之間的距離在x和y軸向上的投影距離。

在實(shí)際海洋方向譜估計(jì)過(guò)程中，記錄不同測(cè)點(diǎn)位置處的時(shí)間序列，然后求得不同測(cè)點(diǎn)之間的互譜。由方程(1)可知，方程左側(cè)的互譜是已知量，當(dāng)測(cè)量物理量選定后，不同物理量對(duì)應(yīng)的傳遞函數(shù)H通過(guò)線性波浪理論是已知的，方程右側(cè)只有方向譜是未知的，方程(1)求逆運(yùn)算即可進(jìn)行方向譜估計(jì)，求逆運(yùn)算過(guò)程中使用不同理論可以得到不同的方向譜估計(jì)值，所以方向譜估計(jì)算法也不是唯一的［1］。

2 傾斜條件下ADCP方向譜估計(jì)方法

2.1 ADCP旋轉(zhuǎn)矩陣

坐底式ADCP俯視圖如圖1所示:

(1)儀器坐標(biāo)系(X，Y，Z):當(dāng)從ADCP的表面上看的時(shí)候，逆時(shí)針?lè)较蛴?-1-4-2。X軸的方向從波束1指向波束2，Y軸是波束4指向波束3，Z軸在4個(gè)波束的中心，從水面表層指向水底。

(2)船舶坐標(biāo)系(S，F(xiàn)，M)，選用右舷，前進(jìn)，桅桿方向分別代表三個(gè)坐標(biāo)系軸。

(3)大地坐標(biāo)系(E，N，U)，分別為正東，正北，垂直向上，坐標(biāo)系的方向都是固定的，所以這種坐標(biāo)系在ADCP中是最常用的參考坐標(biāo)系。

姿態(tài)傳感器測(cè)量的是艦船坐標(biāo)系相對(duì)于大地坐標(biāo)系的角度(見(jiàn)圖2)。在右手原則下，分別相對(duì)于S軸和F軸逆時(shí)針轉(zhuǎn)動(dòng)，Tilt1(pitch)和Tilt2(roll)值為正值。

圖1 坐底式ADCP俯視圖Fig.1 Top view of upward-looking ADCP

圖2 ADCP姿態(tài)測(cè)量Fig.2 Posture measurement of ADCP

大地坐標(biāo)系的三個(gè)軸向是固定的，一般選為參考的坐標(biāo)系。當(dāng)無(wú)傾斜的時(shí)候，儀器坐標(biāo)系和大地坐標(biāo)系是在一個(gè)水平面內(nèi)，但ADCP存在傾斜時(shí)，就不在同一個(gè)水平面內(nèi)，所以大地坐標(biāo)系和儀器坐標(biāo)系存在一個(gè)旋轉(zhuǎn)矩陣。因?yàn)閞oll測(cè)量的是船舶坐標(biāo)系下相對(duì)于大地坐標(biāo)系的夾角而不是儀器坐標(biāo)系和大地坐標(biāo)系的夾角，所以當(dāng)使用坐底式 ADCP 時(shí)，roll首先要加上 180°［6－8］。

將姿態(tài)角組成旋轉(zhuǎn)矩陣

其中 CH=cos(H)，SH=sin(H)，CP=cos(P)，SP=sin(P)，CR=cos(R)，SR=sin(R)，其中 H，P，R 分別代表姿態(tài)角heading，pitch，roll。由旋轉(zhuǎn)矩陣的定義可知，應(yīng)該滿足［1］

2.2 ADCP測(cè)點(diǎn)傾斜修正方法

ADCP可以利用不同深度上速度測(cè)點(diǎn)組成測(cè)點(diǎn)陣列，然后利用速度測(cè)點(diǎn)陣列通過(guò)復(fù)雜的反演算法來(lái)進(jìn)行方向譜估計(jì)。ADCP一般選取接近表層的3層12個(gè)測(cè)點(diǎn)來(lái)進(jìn)行方向譜的估計(jì)，12個(gè)測(cè)點(diǎn)組成了一個(gè)空間測(cè)點(diǎn)陣列，可以利用傳統(tǒng)傳感器陣列方法來(lái)進(jìn)行方向譜的估計(jì)。如圖3，當(dāng)ADCP存在傾斜時(shí)，ADCP的4個(gè)波束 1，2，3，4 的位置變化到了 1'，2'，3'，4'，對(duì)應(yīng)的 ADCP 測(cè)點(diǎn)陣列變成了一個(gè)傾斜的陣列。

圖3 傾斜的ADCP測(cè)點(diǎn)陣列Fig.3 Measurement points array of declined ADCP

一般情況下方向譜估計(jì)值都是在大地坐標(biāo)系下，所以首先要求傾斜后陣列中各點(diǎn)在大地坐標(biāo)下的坐標(biāo)值。假設(shè)無(wú)傾斜的ADCP的4個(gè)波束和垂直方向的夾角為α，選取的3層測(cè)點(diǎn)深度 depth分別為 bin1，bin2，bin3。從圖3可知，對(duì)應(yīng)depth有l(wèi)ength=depth*tan α，則該層波束 1，2，3，4 在大地坐標(biāo)系下的測(cè)點(diǎn)坐標(biāo)［x0，y0，z0］分別為(－ length，0，depth)，(length，0，depth)，(0，length，depth)，(0，－length，depth)。當(dāng) ADCP 傾斜時(shí)，對(duì)應(yīng)的測(cè)點(diǎn)位置變化到 1'，2'，3'，4'，當(dāng)傾斜較大時(shí)，ADCP 的測(cè)點(diǎn)位置偏移會(huì)導(dǎo)致方向譜估計(jì)出現(xiàn)較大誤差，所以測(cè)量點(diǎn)的位置坐標(biāo)要進(jìn)行修正［8－9］，傾斜后位置坐標(biāo)［x0'，y0'，z0'］可以通過(guò)下列表達(dá)式獲取

由此可以得到 ADCP 傾斜后的測(cè)點(diǎn)位置為［x0'，y0'，z0'］［9－11］

當(dāng)知道傾斜后的坐標(biāo)［x0'，y0'，z0'］，就可以得到對(duì)應(yīng)［x0'，y0'，z0'］的傳遞函數(shù) H'，ADCP 方向譜的傾斜修正本質(zhì)是計(jì)算修正后的測(cè)點(diǎn)位置［x0'，y0'，z0'］和相應(yīng)傳遞函數(shù)H'，然后進(jìn)行方向譜計(jì)算。ADCP利用不同深度上測(cè)點(diǎn)組成測(cè)點(diǎn)陣列進(jìn)行方向譜估計(jì)是ADCP目前反演海洋表面波最常用的方法，這種方法具備高精度，同時(shí)，ADCP波束方向的速度值幾乎在任何情況下都能得到，這種方法可以在大部分海洋條件下進(jìn)行工作［12－15］。

2.3 算法實(shí)現(xiàn)流程

當(dāng)ADCP傾斜時(shí)，修正的目的就是為了降低傾斜帶來(lái)的方向譜估計(jì)誤差。本文仿真和實(shí)測(cè)數(shù)據(jù)表明，ADCP姿態(tài)角pitch，roll和ADCP安放深度都能影響到方向譜估計(jì)值。所以在方向譜的傾斜修正中，RD公司只考慮固定姿態(tài)角10°門限(圖4a)，顯然是不能滿足實(shí)際需要的，所以為了得到準(zhǔn)確的估計(jì)，考慮姿態(tài)角pitch，roll的同時(shí)不應(yīng)該考慮ADCP的安放深度。圖4b是本文給出的一種傾斜修正算法，計(jì)算如下:

(1)首先設(shè)立傾斜修正原則，本文設(shè)立的原則是當(dāng)ADCP傾斜帶來(lái)有義波高H1/3變化大于指定常數(shù)時(shí)，需要傾斜修正。其中指定常數(shù)要大于ADCP儀器的有義波高測(cè)量精度，本文設(shè)定該指定常數(shù)不小于1%。

(2)在對(duì)一組ADCP數(shù)據(jù)進(jìn)行方向譜估計(jì)之前，首先要考慮到ADCP安放深度，以及姿態(tài)角pitch，roll值(其中姿態(tài)角pitch，roll要求大于儀器的測(cè)量精度，一般為2°)。

(3)通過(guò)仿真得到給定的ADCP安放深度和姿態(tài)角pitch，roll條件下的方向譜估計(jì)誤差RMSE，然后通過(guò)方程(8)得到有義波高變化ΔH1/3/H1/3。

(4)將仿真得到的有義波高變化根據(jù)傾斜修正原則判斷是否需要進(jìn)行傾斜修正。當(dāng)需要修正時(shí)，首先對(duì)各測(cè)點(diǎn)位置進(jìn)行傾斜修正，然后計(jì)算方向譜;而不需要修正時(shí)，不需要進(jìn)行測(cè)點(diǎn)位置修正，可直接計(jì)算方向譜。

圖4 ADCP傾斜條件下方向譜計(jì)算Fig.4 Wave directional spectrum estimation for tilting ADCP

3 仿真分析

仿真條件:(1)在水深15 m和25 m的情況下，選取陣列形式為JANUS配置的3層12點(diǎn)陣列，每層測(cè)點(diǎn)之間的為0.5 m;(2)假設(shè)姿態(tài)角 pitch，roll從0°～10°變化，計(jì)算不同水深下方向分布估計(jì)誤差RMSE隨著傾斜姿態(tài)角pitch，roll的變化，S(f，θ)和^S(f，θ)分別代表了方向譜真實(shí)值和方向譜估計(jì)值，H1/3代表有義波高。

圖5選取了深度為10 m，25 m的2種情況?？梢钥闯霎?dāng)ADCP深度為10 m或25 m時(shí)，ADCP傾斜姿態(tài)角pitch，roll越大，方向譜估計(jì)誤差RMSE越大;在相同的姿態(tài)角pitch，roll下，ADCP安放深度25 m的方向譜估計(jì)誤差RMSE要大于ADCP安放深度為10 m的方向譜估計(jì)誤差RMSE，即隨著深度的增大，相同的姿態(tài)角pitch，roll帶來(lái)的方向譜估計(jì)誤差變大。這種現(xiàn)象在理論上是可以解釋的，假設(shè)無(wú)傾斜時(shí)的ADCP中心點(diǎn)處的水平面內(nèi)坐標(biāo)為(0，0)，ADCP放置深度為h，傾斜條件下姿態(tài)傳感器pitch，roll值使得中心點(diǎn)位置(0，0)在水平面內(nèi)分別移動(dòng)h×sin(pitch)，h×sin(roll)，此變化位移不只跟姿態(tài)角pitch，roll值有關(guān)，同樣跟ADCP放置深度h有關(guān)，當(dāng)ADCP深度h越大，相同姿態(tài)角pitch，roll帶來(lái)的測(cè)點(diǎn)位置偏移越大，導(dǎo)致方向譜估計(jì)誤差RMSE越大。

圖5 方向譜估計(jì)誤差RMSE隨姿態(tài)角變化Fig.5 Variation of wave directional spectrum estimation error RMSE with tilting angles

圖6給出了由公式(8)計(jì)算得到的有義波高變化ΔH1/3/H1/3隨pitch，roll的變化關(guān)系曲線，從中可以看出:(1)有義波高變化ΔH1/3/H1/3與方向譜估計(jì)誤差RMSE變化規(guī)律相似;(2)當(dāng)姿態(tài)角pitch，roll等于10°時(shí)，ADCP安放深度為25 m時(shí)，對(duì)應(yīng)有義波高變化ΔH1/3/H1/3分別為0.052和0.05，而當(dāng)ADCP安放深度為10 m時(shí)，對(duì)應(yīng)有義波高變化ΔH1/3/H1/3分別為0.18和0.18，即ADCP安放深度不同，相同的姿態(tài)角帶來(lái)的有義波高變化不同。RD公司給出傾斜修正門限為10°，沒(méi)有考慮到ADCP安放深度，這是不合理的。

圖6 有義波高百分比ΔH1/3/H1/3隨傾斜角度變化Fig.6 Variation of significant wave height percentage ΔH1/3/H1/3with tilting angles

考慮實(shí)測(cè)數(shù)據(jù)的有義波高的精度一般為1%，若假定實(shí)測(cè)數(shù)據(jù)因?yàn)閮A斜帶來(lái)的有義波高變化ΔH1/3/H1/3為1%作為傾斜修正門限時(shí)，仿真表明當(dāng)水深為25 m的情況下，姿態(tài)角pitch，roll大于4°即需要進(jìn)行傾斜修正，此時(shí)RD公司的姿態(tài)角10°門限顯然是不能滿足要求的。

4 實(shí)測(cè)數(shù)據(jù)分析

選取RD公司ADCP的實(shí)測(cè)數(shù)據(jù)進(jìn)行傾斜修正分析。數(shù)據(jù)1和2來(lái)源于北緯34.34°，西經(jīng)76.42°的美國(guó)卡羅萊納州lsrb試驗(yàn)點(diǎn)，選取的是600 kHz坐底式的ADCP，置于水下25.4 m，采樣頻率是2 Hz。數(shù)據(jù)3和4來(lái)源于美國(guó)加利福尼亞州，1 200 kHz坐底式的ADCP位于水下10.9 m，采樣頻率為2 Hz。

圖7中數(shù)據(jù)1和數(shù)據(jù)2的深度為25.4 m，同仿真中25 m進(jìn)行對(duì)比，可以看出實(shí)測(cè)數(shù)據(jù)1，2與仿真情況的有義波高百分比ΔH1/3/H1/3隨著姿態(tài)角pitch，roll變化的曲線是吻合的。同樣在圖8中數(shù)據(jù)3和數(shù)據(jù)4的深度為10.9 m，與仿真中10 m進(jìn)行對(duì)比分析，實(shí)測(cè)數(shù)據(jù)3和實(shí)測(cè)數(shù)據(jù)4隨著姿態(tài)角pitch，roll變化曲線與仿真10 m的曲線是吻合的。實(shí)測(cè)數(shù)據(jù)1、2、3和4與仿真結(jié)論吻合，證明了仿真分析的正確性。

圖7 數(shù)據(jù)1和數(shù)據(jù)2的仿真和實(shí)測(cè)對(duì)比Fig.7 Simulation and experiment comparisons between data 1 and data 2

圖8 數(shù)據(jù)3和數(shù)據(jù)4的仿真和實(shí)測(cè)對(duì)比Fig.8 Simulation and experiment comparisons between data 3 and data 4

5 結(jié)語(yǔ)

當(dāng)ADCP存在傾斜時(shí)，ADCP的測(cè)點(diǎn)位置發(fā)生偏移，傾斜帶來(lái)的方向譜估計(jì)誤差影響ADCP波浪估計(jì)的準(zhǔn)確性。本文給出了一種傾斜修正算法，通過(guò)傾斜帶來(lái)的有義波高變化來(lái)判斷是否進(jìn)行測(cè)點(diǎn)傾斜修正，其中有義波高隨著姿態(tài)角和ADCP安放深度的定量變化是可以預(yù)先通過(guò)仿真得到的，新算法實(shí)現(xiàn)簡(jiǎn)單，且提高了當(dāng)ADCP傾斜時(shí)的方向譜估計(jì)精度。

［1］BRUMLEY B H，TERRAY E A，STRONG B S.System and method for measuring wave directional spectrum and wave height:US，US 6052334［P］.2001 －08 －28.

［2］LONHRMANN A，PEDERSEN T K.System and method for determining directional and non-directional fluid wave and current measurements:US，US 7352651［P］.2008 －10 －02.

［3］KROGSTAD H E，GORDON R L，MILLER M C.High-Resolution directional wave spectra from horizontal mounted acoustic doppler current profiler meters［J］.Journal of Atmospheric and Oceanic Technology.1988，5:340 － 352.

［4］ZEDEL L.Deep ocean wave measurement using a vertically oriented sonar［J］.Journal of atmpheric and oceanic technology.1994，11(1):182 －191.

［5］Teledyne RD Instruments.WaveMon v3.06 User’s Guide［Z］.US:RD Instruments，2009.

［6］SMITH J A，BULLARD G T.Directional surface wave estimates from Doppler sonar data［J］.Journal of atmospheric and oceanic technology，1995，12(3):617－632.

［7］VISBECK M，F(xiàn)ISCHER J.Sea surface conditions remotely sensed by upward-looking ADCPs［J］.Journal of atmpheric and oceanic technology.1995，12(1):141－149.

［8］TERRAY E A，KROGSTAD H E，CABRERA R，et al.Measuring wave direction using upward-looking Doppler sonar［C］//Proceedings of the IEEE 4th working conference on current measurement.US:IEEE，1990:252 －257.

［9］RD instruments.ADCP coordinate transformation:Formulas and calculations［Z］.US:RD Instruments，2010.

［10］TREVORROW M V，BOOTH L J.Extraction of Ocean wave directional wave spectra using steerable Doppler side-scan sonar［J］.Journal of atmospheric and oceanic technology，1995，10:1087 －1100.

［11］TREVORROW M V.Measurement of ocean wave directional spectra using Doppler side-scan sonar arrays［J］.Journal of atmospheric and oceanic technology.1995，16:603 －616.

［12］TRUMP C L，ALLAN N.，MARMORINO G S.Side-looking ADCP and Doppler radar measurements across a coastal front［J］.IEEE Journal of Oceanic Engineering，2000，25(4):423 －429.

［13］TERRAY E A，BRUMLEY B H，STRONG B.Measuring waves and currents with an upward-looking ADCP［C］//Proceedings of the IEEE Sixth Working Conference on Current Measurement.San Diego，CA，USA:IEEE ，1999:66 －71.

［14］STRONG B S.System and method for determining wave characteristics from a moving platform:US，US20100302908A1［P］.2010－12－02.

［15］TRENAMAN N，DEVINE P，STRONG B.ADCP-based multi-directional wave gauge and current profilinig［C］//OCEANS ＇02 MTS.US:IEEE ，2002:1763－1766.