用于多普勒流速剖面仪测波浪的方向谱反演算法研究

李 晨 吴建波 高 超 王道龙 李 俊 邓 锴 王长红

①(中国科学院声学研究所 北京 100190)

②(国家海洋局第一海洋研究所 青岛 266003)

1 引言

海浪的方向谱给出海浪能量在频率和方向上的2维分布，细致地描述了海浪能量的传播信息。使用声学多普勒流速剖面仪(Acoustic Doppler Current Profiler, ADCP)测波浪是目前精度较高、应用较广的一种海浪测量手段，它利用各换能器空间位置不同造成的相位差来计算波向[1]。

ADCP测量波浪计算方向谱包括流速反演和波面高度反演两种方法，其中波面高度反演由于直接测量波面变化，不需任何中间转换，也无需考虑洋流影响[2]，而且水气界面声波反射强度明显比水质点大，因此测量数据精确度较高，与雷达测波仪[3]相比测量精度与测量范围均有一定优势。

目前常见的几种基于仪器阵列的海浪方向谱反演算法有最大似然法(MLM)、扩展本征矢(EEV)法和贝叶斯方向谱估计法(BDM)等。最大似然法是波浪方向谱分析中最常用的方法，具有计算速度快、计算结果稳定、对噪音不敏感等特点，同时具有一定的数值精度；但当仪器阵列距离反射面较近时，其方向分辨率会受到一定影响[4]。扩展本征矢法能够得到能量分布更加集中于主波向的方向谱[5,6]，但其适用范围有待进一步论证。贝叶斯法能够使能量在方向上的分布更平滑，获得更加可靠、精度更高的方向谱，前提是阵列元素个数不小于 4，交叉谱质量较好且仪器布放深度不宜过大[7]；其主要问题是计算时间长。美国RDI公司[1]和挪威Nortek公司[8]的ADCP波浪数据反演都使用了最大似然法。另外还有 SUV(Surface height, transverse velocity component(U)and longitudinal velocity component(V))法等，由于其只适用于含有中央垂向波束的情况[9]，因此本文暂未讨论。

本文通过计算机仿真和实测数据处理两方面工作，将目前常见的几种基于仪器阵列的海浪方向谱反演算法用于4波束Janus配置的ADCP波面高度测量数据的处理并进行分析比较，确定它们在ADCP波浪方向谱反演中能够取得的反演效果和适用范围。本文选择最具普遍性的 MLM、国内专家提出的EEV和国际公认精确度较高的BDM进行比较，具有典型意义；本文对3种算法进行了理论上的误差分析与对比；在计算机仿真中，引入不同的仪器布放深度(影响波束投影位置、范围与间距)作为变化因素进行理论分析和仿真结果的比较，使仿真更具实际意义；在实测数据处理中，一方面将RDI ADCP数据分别使用WavesMon软件和我们自主编写的算法进行反演并与比测仪器进行结果比较，另一方面将声学所ADCP数据的反演结果与比测仪器进行比较，证明了我们自主研发的仪器和配套算法具有良好的波浪测量和反演精度。

2 算法原理

使用有限数据估计海浪方向谱的基本思想是，由同一物理量(波面高度或波生流速)在不同位置的时间序列(波束阵列，如ADCP情形)或者浪场中某一点处不同物理量的时间序列(定点测量，如浮标情形)经傅氏变换构成频域交叉谱。波束阵列的交叉谱矩阵包含了同一频率的波浪对于不同位置波束的相对相位信息，这正是方向谱反演的关键。

波束阵列的交叉谱与方向谱的关系可以表示为

2.1 最大似然法

适用于波束阵列的最大似然法(Maximum Likelihood Method, MLM)基本思想是将方向谱估计值表示为交叉谱的线性叠加，并将此估计值最小化来接近真实值。方向谱估计值(k,ω)的表达式为

为了增加能量的集中程度，提高波向辨识度，在最大似然法中引入迭代。具体做法是用MLM算得的方向谱估计值0代入式(1)还原交叉谱，再用此交叉谱重新估计方向谱得到，最后通过

2.2 扩展本征矢法

2.3 贝叶斯法

Hashimoto等人[7]将贝叶斯模型引入海浪方向谱的估计，称其为贝叶斯方向谱估计法(Bayesian Directional spectrum estimation Method, BDM)，简称贝叶斯法。考虑误差ε，将交叉谱矩阵中的元素φi改写为

其中εi相互独立，且遵从均值为 0、方差为σ2的高斯分布，N=M×M,K为方向离散点的个数。则最优的超参数u和方差σ2通过将赤池贝叶斯信息量标 准(Akaike’s Bayesian Information Criterion,ABIC)最小化求得

其中L(x,σ2)为给定φi时x和σ2的似然函数，p(x|u2,σ2)为x的先验分布。

2.4 理论误差分析

(1)在方向谱反演算法中，将频率和方向分别进行离散化选点，所得方向谱实际上是由离散点数决定行(频率)数和列(方向)数的2维矩阵。本文3种算法均采用逐行计算的方法，即依次对每个频率离散点的方向分布进行解算；但在计算方向分布时，MLM和EEV对各个方向离散点的计算是孤立的，而BDM是通过解一个形如Ax=b的非齐次线性方程组来实现本频率点的方向分布解算。显然，在BDM的求解方法下，方向分布的相关性更强，随机误差的影响更小。

(2)为保证方向分布光滑连续，BDM 要求方向分布向量x(f)=[x1(f),x2(f), …,xK(f)]中的元素满足xk-2xk-1+xk-2≈0,2＜k≤K；这一限定条件即包含在方程组Ax=b中。而MLM和EEV并无类似的平滑操作。

(3)BDM 计算x(f)是采用迭代法逐步逼近真实值，连续两个x(f)的差向量的模值小于 10-3时迭代结束，以此计算ABIC，并通过改变超参数u的值重新计算x(f)，即在外层也建立迭代，选出最小的ABIC值和对应的x，从而降低误差水平；如2.1节所述，MLM 也使用迭代对方向谱结果进行修正；而EEV并未引入修正算法。

(4)BDM 中非齐次线性方程组Ax=b求解x时，方程组有解的条件是增广矩阵(A,b)的秩等于系数矩阵A的秩，而A是行数和列数均为K的方阵；显然当K值较大时，方程组无解的可能性(即无法得到当前频率点的方向分布)较大。这是BDM的一个显著缺陷。

(5)EEV把交叉谱矩阵分成信号和噪声两部分，认为本征值明显较大的部分为信号，其余为噪声；MLM 则认为整个交叉谱矩阵均含有噪声。显然MLM的划分方法更符合实际情况；而EEV将一部分噪声也当作信号，这增加了信号部分的能量，使得方向谱反演结果的谱峰更加尖锐，但相对实际情况而言引入了人为误差。

由上述分析可知，BDM采取的误差消除措施最多，在方向分布有解的情况下，反演结果应最逼近真实值；MLM 进行了迭代计算，没有平滑处理；EEV未引入修正和平滑算法，其结果误差应为最大。BDM的缺陷是方向分布会有无解的情况发生，且由于其需要进行内外两层迭代，因此运算耗时较长；而EEV具有谱峰尖锐的优势。

3 计算机仿真

本文仿真和试验采用的波浪观测方法是，ADCP布放在海底，向上发射声波波束，通过接收回波获得水面高度变化的时间序列。

计算机仿真的步骤是：构造海面3维高度模型；根据仪器布放深度确定波束阵列水面投影的相对位置和范围，采集投影内的水面高度数据；使用不同反演算法处理数据获得方向谱结果。

3.1 海面模型

海浪不仅波高、频率不尽相同，而且会沿各个方向传播，除主波向以外，在其两侧±π/2范围内都有谐波的扩散。根据线性海浪理论，可用双叠加法模型(来源于Pierson模型)模拟3维波面[10]：

其中S(ωi,θj)为 2维海浪波谱，m,n分别为频率和方向的离散化点数，ki为波数，ωi为角频率，θj为方向角，ξ和η分别为海平面上各个离散点的横坐标和纵坐标，εij为随机初相位，是在0～2π间均匀分布的随机变量。线性海浪理论应用于工程时，通常认为海浪能量的频率分布和方向分布相互独立，因此2维海浪谱可用1维海浪谱函数S(ω)和方向谱函数D(θ)的乘积来表示。

我们采用的1维海浪谱为 PM(Pierson-Moskowitz)谱，由于其数据基础较为坚实，迄今仍是最常用的海浪频谱模型之一。采用ISSC(国际船舶结构会议)建议的方向谱函数：

其中θ为组成波传播方向相对风向的角度。

3.2 仪器布放

波束阵列测波浪根据各波束的相对位置和测得波浪的相对相位确定波长和波向等参数。当所测波浪的波长小于足印最小间距的2倍时，波束阵列将无法解出波向，这就是空间的奈奎斯特定理。在方向谱中，波长过小(或频率过高)的波浪在正确波向以外的区域(即谱峰波向区域以外)会出现能量分布，这是波向解算错误在方向谱中的表现形式。而当各组成波的波长一定时，随着足印间距的增加，能够解出波向的波浪最小波长λmin随之增加，方向谱反演结果中谱峰波向区域以外的能量分布范围由高频向低频扩展，此区域能量幅值在频率上的分布规律也与谱峰区域一致。这就是足印间距变化对方向谱能量分布造成的影响，在后面的仿真和实测结果图中都得到了体现。

凹陷路面是典型的问题路面，一旦路面出现了因压迫或是其他原因所产生的凹槽，便必然会出现承重问题。凹槽对于路基形成了一定程度上的破坏，路基的整体性能便会受到影响。所以不能只是简单的对裂缝进行修补，而是要针对基层被损坏的部分去进行二次修建。要真正做好养护，必须要对路基的凹陷段以及翻浆的路段采取相应的治理措施去进行优化处理，在底层部分，需要使用水泥去对碎石砂（6%）开展垫层操作，稳定结构，并且要进行反复的进行碾压，直到结构完全成型。在基层的修复结束后，可以利用再生沥青对表层部分进行处理，这样能够得到更加理想的修复效果，尽可能避免凹陷问题的产生。

根据各算法的原理可知，足印间距只是式中的一个常数项，其值的改变并不会对算法本身引起的误差造成影响；但由 2.4节可知，MLM, EEV和BDM对整体误差的消除程度不同，因此由足印间距变化引起能量分布的变化对各算法的影响不同，各算法反演方向谱的谱峰区域和谱峰以外区域的能量分布都会存在差异。

本文仿真采用的ADCP波束阵列分布方式是经典的4波束Janus配置，波束的倾角θBeam=20°，波束开角θSensor=2°。如果令D为每个波束在水面投影(即“足印”)中心到整个波束阵列中心点的距离，则仪器的布放深度z与D的关系为D=z·tanθBeam。在每个以“足印”中心为圆心，半径r=z·tanθSensor的水平圆面内采集波面高度数据，然后计算均值作为该测点当前时刻的水面高度测量值，因此仪器布放深度决定了各个波束投影的位置和间距。

由上面分析知，仪器布放深度与足印间距成正比，因此其与方向谱反演效果之间存在联系，可通过改变其值对各种反演算法进行较全面的性能比较。布放深度z与足印最小间距Bmin的关系为

计算机仿真中的参数设置如表1所示，其中λ为所测波浪的谱峰波长。

表1 计算机仿真各参数设置

3.3 仿真结果

波浪方向谱以极坐标的形式表示，用颜色变化区分能量值的大小，单位为 m2/Hz/rad；为了便于比较，规定谱峰(即整个方向谱的能量最大值点)处能量为 1，以此将各方向谱进行能量归一化，则所得谱图即表示能量的分布情况和谱峰的尖锐程度。矢径表示频率，单位为Hz；极角表示方向，取值范围为0～360°, 0°为正北，90°为正东，为符合日常观察习惯，将坐标作了适当翻转。

图1为仿真数据的方向谱比较，图2为仿真数据的能量在方向上1维分布比较。

图1 仿真数据的方向谱比较，风速10 m/s

图2 仿真数据的能量在方向上1维分布比较，风速10 m/s

3.4 结果分析

通过计算机仿真结果的比较可看出：

(2)EEV反演结果的谱峰最为尖锐，谱峰波向准确度最高；但峰值超出理论值较多，λ/2Bmin＜2时谱峰区域以外的能量分布较多。

(3)MLM 反演结果与理论能量分布的拟合度略差于BDM，但也能够得到较为准确的谱峰区域和明确的谱峰波向，结果比较稳定；只是λ/2Bmin＜2时谱峰区域以外会出现较为明显的能量分布。

4 实测数据处理

由于计算机仿真与实际情况之间必然会存在差异，因此为了考察各反演算法对于ADCP实测数据的适用性，我们分别使用RDI公司生产的600 kHz ADCP和中国科学院声学研究所自行研制的 150 kHz ADCP工程样机在青岛近海进行波浪观测试验，通过分析其波面高度测量数据，对3种方向谱反演算法进行比较。仪器布放时考虑了空间的奈奎斯特定理，保证波束在水面投影的最小间距小于所测波浪波长的一半。两台仪器的换能器分布方式均为4波束Janus配置，波束的倾角为20°，波束开角为2°。

由于布放海域很难保证底部完全平坦，因此仪器实际测量时会产生布放偏差，即仪器布放倾斜导致各个波束在水面投影的相对位置出现偏移。为消除此误差，在处理实测数据时，我们考虑仪器的实际姿态，根据两个倾斜角纵摇(pitch)，横摇(roll)计算各波束在水面的实际投影位置；只要保证各“足印”位置计算准确，其在一定程度内的相对位移变化(如纵摇，横摇均在10°以内)就不会对反演结果的准确性产生影响。我们选择波浪大小适宜、倾斜角较小的两组实测数据进行分析比较，以保证其与计算机仿真结果具有对照意义。

试验的海浪波长、仪器布放深度等参数如表 2所示，可见相对2011年1月5日的试验条件而言，1月9日波浪的波长和λ/2Bmin的值都较大。

表2 试验相关参数

由于本章方向谱的对比中使用了其他波浪数据处理软件的截图以及比测仪器反演结果的转化图，因此本章中的方向谱图未进行统一尺度或能量归一化操作，能量幅值的单位均为m2/Hz/rad。

4.1 ADCP实测数据反演结果

首先利用RDI公司的600 kHz ADCP测得的波浪高度数据，对3种方向谱反演算法进行比较。

因为 RDI有自己的波浪数据处理软件 Waves Mon，所以我们可使用其反演结果作为比较标准；同时布放的还有 Datawell公司生产的波浪浮标(DWR MKIII波浪骑士)，由于其测量精度高，被称为“测波的标准仪器”，因此也可使用其结果作为方向谱反演算法的性能衡量指标。结果如图3所示。

其次利用中国科学院声学研究所自行研制的150 kHz ADCP工程样机测得的波浪高度数据，对3种方向谱反演算法进行比较。使用同时布放的波浪骑士的反演结果作为各算法的性能衡量指标。结果如图4所示。

图3 实测数据的方向谱和能量在方向上1维分布比较，2011年1月5日13:10

图4 实测数据的方向谱和能量在方向上1维分布比较，2011年1月9日12:00

4.2 结果分析

综合分析上述两组实测结果可看出：

(1)实际测量中 BDM 仍能取得最佳的反演结果，能量分布最集中，谱峰区域以外的能量分布较少。但出现无解频率点的现象仍会发生；其计算时间仍为MLM和EEV的数百倍。

(2)实际测量中 EEV 在λ/2Bmin较大的情况下体现出其谱峰较为尖锐的特性(但在能量在方向的1维分布中其谱峰尖锐程度与BDM相比略差)，且波向准确度较高。在λ/2Bmin较小的情况下其谱峰尖锐程度较差，能量分布最为分散，说明EEV的适用范围小于BDM和MLM。

(3)MLM 的反演结果在波浪波长较小的情况下比BDM略差而明显好于EEV，在波长较大的情况下与EEV相比各有优势，因此在实测数据的反演计算中表现出较高的计算效率和稳定性。

5 结束语

综上所述，通过MLM，EEV和BDM 3种波浪方向谱反演算法处理ADCP波面高度数据的结果对比，可以得出以下结论：

(1)根据计算机仿真和实测数据处理的结果得出的结论并不完全相同(如 EEV 方向谱的尖锐程度)，这是因为：一方面仿真模拟的是理想状况，与实际海况相比存在差异，也说明现实海浪的复杂性；另一方面并未考虑声波在传播过程中的噪声、散射等干扰因素对数据质量造成的影响。

(2)3种算法中，BDM能够最为准确地再现波浪能量在方向上的分布，但由于其计算复杂度较高，因此会导致无解情况的发生(即对数据质量的要求较高)，而且计算时间远远大于另外两种算法，即运算效率较低，因此适用于精度要求较高的计算。

(3)EEV能够获得明确尖锐的谱峰，波向准确度较高，但波浪波长偏小时反演结果相对较差，因此适用范围比较有限，可在数据处理时作为参考。

(4)MLM 对计算机仿真数据和实测数据的反演都能够取得令人较为满意的结果，计算稳定度和效率最高，较适合应用于ADCP波面高度反演；其主要问题是谱峰区域以外会出现一定程度的能量分布，且要注意仪器布放时不能过分接近水面。

算法仍有改进的空间，如 BDM如何减少频率点无解情况的发生；EEV如何扩大适用范围；MLM如何减少谱峰区域以外的能量分布等。

[1]Brumley B H, Terray E A, and Strong B. System and method for measuring wave directional spectrum and wave height[P].US, Patent, No. 6052334, 2000.

[2]董兆乾, 蒋松年, 贺志刚. 南大洋船载走航式ADCP资料的技术处理和技术措施以及多学科应用[J]. 极地研究, 2010, 22(3):211-230.

Dong Zhao-qian, Jiang Song-nian, and He Zhi-gang. ADCP data technical processing of the CHINAREs and its multidisciplinary application[J].Chinese Journal of Polar Research, 2010, 22(3): 211-230.

[3]楚晓亮, 徐铭, 王峰, 等. 利用X波段雷达提取海浪信息的分析[J]. 中国海洋大学学报, 2011, 41(5): 110-113.

Chu Xiao-liang, Xu Ming, Wang Feng,et al.. Analysis of the wave information extracted by X band radar[J].Periodical of Ocean University of China, 2011, 41(5): 110-113.

[4]Isobe M and Kondo K. Method for estimating directional wave spectrum in incident and reflected wave field[C]. 19th International Conference on Coastal Engineering, Houston,TX, Sept. 1984: 467-483.

[5]管长龙, 文圣常, 张大错. 分析海浪方向谱的扩展本征矢方法I: 方法的导出[J]. 海洋与湖沼, 1995, 26(1): 58-62.

Guan Chang-long, Wen Sheng-chang, and Zhang Da-cuo. An extended eigenvector method for estimation of directional spectra of sea waves I: derivation of the method[J].Oceanologia et Limnologia Sinica, 1995, 26(1): 58-62.

[6]管长龙, 文圣常, 张大错. 分析海浪方向谱的扩展本征矢方法II: 方法的验证、比较和应用[J]. 海洋与湖沼, 1995, 26(3):241-246.

Guan Chang-long, Wen Sheng-chang, and Zhang Da-cuo. An extended eigenvector method for estimation of directional spectra of sea waves II: verification, comparison and application of the method[J].Oceanologia et Limnologia Sinica, 1995, 26(3): 241-246.

[7]Hashimoto N and Kobune K. Directional spectrum estimation from a Bayesian approach[C]. 21st International Conference on Coastal Engineering, Coata del Sol-Malaga,Spain, June 1988: 62-76.

[8]Siegel E, Pedersen T, and Maatje J. Real-time directional wave measurements: innovative engineering for subsurface acoustic Doppler current profilers[EB/OL]. http://www.nortekusa.com/lib/technical-notes/awac-introduction-in-sea-technology-february-2006, Feb. 2006.

[9]李亚光, 吴建波, 高超, 等. 基于 AWAC型 ADCP的波浪反演算法研究[J]. 海洋技术, 2010, 29(3): 55-58.

Li Ya-guang, Wu Jian-bo, Gao Chao,et al.. Research on wave direction measurement using AWAC ADCP[J].Ocean Technology, 2010, 29(3): 55-58.

[10]文圣常, 余宙文. 海浪理论与计算原理[M]. 北京: 科学出版社, 1984: 127-130, 276-278.

Wen Sheng-chang and Yu Zhou-wen. Ocean Theory and Computing Principles[M]. Beijing: Science Press, 1984:127-130, 276-278.