海浪的观测对海洋工程和海洋科学研究有重要意义。传统的浮标等观测方式只能获得单点的波浪信息; 随着卫星遥感技术的发展, 雷达高度计和合成孔径雷达等也被应用于海浪的观测, 它们能够获得大面积的波浪信息, 但是其时间和空间分辨率较差, 不适用于近岸海区。近年来, X波段海洋雷达被应用于海表物理参数的反演。1965年, Wright利用X波段雷达图像直接判读海浪的传播方向和波长^[1]。1985年, Young等^[2]提出利用三维傅里叶变换处理X波段雷达图像序列的算法, 并利用波数谱估算了海浪参数。我国学者在谱分析算法方面也做了很多研究和应用^[3-7]; 为了提取非均匀波浪场中的海浪信息, 不同的数据处理方法也被用于分析导航X波段雷达图像, 例如连续小波变换^[8]、经验正交函数分解^[9]等。

由于导航X波段雷达没有经过定标, 利用谱分析等算法反演海浪的有效波高时需要利用现场测量的浮标数据做定标^[9-10]。另外, 由于X波段的电磁波在传播时受降雨的衰减较大, 降雨时X波段雷达的观测精度受到很大影响, 一般情况下只能作为异常值剔除^[11], 这给实际应用带来很多不便, 而目前关于降雨情况下X波段雷达观测海浪的研究较少。本文提出了一种降低降雨影响的算法来反演海浪参数, 并且不需要浮标数据作定标。本文第1部分介绍实验数据的获取方法, 第2部分介绍本文提出的消除降雨对波高反演影响的方法, 第3部分利用实验数据对方法进行验证, 第4部分是对本文工作的讨论和总结。

1 实验和数据 1.1 实验介绍

2016年10月17日~10月18日, 我们在广东省博贺海洋气象观测站附近海域做了观测实验(图 1), 该海区的西北部为陆地、东部和南部为海洋, 其平均水深大约20 m。实验仪器包括一套导航X波段雷达测波系统和一个波浪浮标, 实验中雷达和浮标进行了同步观测, 其安装位置如图 1。浮标的主要测量参数的指标如表 1, 实验中设定浮标每30 min输出一组海浪数据, 包括海浪的有效波高、主波波向和主波周期等。

本实验选择自主研发的导航X波段雷达测波系统, 雷达的主要性能指标如表 2所示。实验中雷达每0.5~1 h观测一组数据, 每组数据包含64幅雷达图像。根据表 2, 天线的旋转周期为2 s, 所以每组数据的采集时间为128 s。雷达的观测范围大约是3 km, 雷达接收海面的后向散射后, 经过专用的数据采集卡保存为灰度图像, 图像中每个像素的分辨率为2.4 m。

1.2 实验数据

实验中观测的两幅雷达图像如图 2所示。图 2a是没有降雨时的雷达图像, 从中可以看到清晰的海浪条纹。图 2b是有降雨时的雷达图像, 由于降雨对电磁波的衰减, 雷达图像中的海浪回波明显减弱, 海杂波的条纹变得模糊。根据气象局的记录, 实验期间的天气为小到中雨。由于实验中没有测量降雨量, 本文根据雷达图像的模糊程度判断每组数据是否受到降雨的影响。共选取一次降雨过程前后的52组数据, 其中受降雨影响的数据共21组。

另外, 实验期间的风向为东南风或者南风, 主波波向为东南方或南方(图 1), 因此观测的波浪主要为风浪。

2 反演方法介绍

本节首先分析了无降雨和降雨时的导航X波段雷达图像, 通过谱分析研究了降雨的影响, 然后提出了消除降雨影响的方法。

2.1 数据分析

以图 2中的雷达图像为例说明降雨对反演海浪参数的影响。图 3a是从图 2a中选取的一个径向区域, 可以看出, 在没有降雨时雷达回波有明显的波峰和波谷, 灰度值的变化范围是0~140, 这是由于海面回波经波浪的倾斜调制、阴影调制等影响而形成的^[10]。对该径向作一维傅里叶变换得到其一维波数谱(图 3b), 谱有明显的峰值, 对应的主波波长为122.88 m, 即无降雨时可以从导航X波段雷达图像中有效提取波浪信息。

图 4a是从图 2b中选取的一个径向区域, 有降雨影响时的雷达回波中包含几个峰值, 但是相邻的波峰之间有多个极小值, 不易准确分辨真实的波谷, 灰度值的变化范围主要集中在40~120。图 4b的一维波数谱表明, 在波数较小的区域有明显的波峰(波数为0.0051 rad/m附近), 而这不是我们要观测的重力波。因此, 从受降雨影响的雷达图像中不易准确提取波浪信息。

为了从降雨情况下的导航X波段雷达图像中提取海浪的波高, 首先用主成分分析提取波浪变化的主成分, 经滤波减小降雨和系统噪声的影响, 用谱分析获得波浪的一维波数谱; 然后选取成熟的经验海浪谱作为标准谱, 建立反演的海浪谱与标准谱之间的模型, 从而估算海浪的有效波高。

2.2 提取波浪变化的主成分

在遥感应用领域, 主成分分析常用于把不同波段中的有用信息集中到数目尽可能少的新的主成分图像中, 达到信息综合与增强的目的。主成分分析得到的主成分之间互不相关, 同时使得原始图像的信息量损失最小^[12]。在降雨情况下, 导航X波段雷达图像中的海浪信息受到雨的影响而变模糊, 本研究用主成分分析提取雷达图像中海浪的主要形态。

选取雷达图像中主波波向附近的几个径向, 将其表示为矩阵:

$X = \left[{\begin{array}{*{20}{c}} {{x_{11}}}&{{x_{12}}}&{...}&{{x_{1n}}}\\ {{x_{12}}}&{{x_{22}}}&{...}&{{x_{2n}}}\\ {...}&{...}&{...}&{...}\\ {{x_{m1}}}&{{x_{m2}}}&{...}&{{x_{mn}}} \end{array}} \right]$

(1)

式中, X为原始数据矩阵, x_ij为雷达图像中第j(j=1, 2, …, m)个径向的第i(i=1, 2, …, n)个像元的灰度值, m为相邻径向数, n为每个径向的像元数。本文选取: m=11, n=512。

经过矩阵运算, 得到原始数据矩阵的主成分为^[13]:

$ \mathit{\boldsymbol{Y}} = {\mathit{\boldsymbol{U}}^\mathit{\boldsymbol{T}}}\mathit{\boldsymbol{X}} $

(2)

式中, Y为变换后的主成分, U为相关系数矩阵的单位特征向量矩阵。原始数据经主成分变换后, 得到一组m个新的变量(Y的各个列向量), 依次称为第一主成分, 第二主成分, …, 第m个主成分。由于第一主成分包含了最多的波浪信息, 本文选取第一主成分来表示波浪的变化。

2.3 滤波

海浪是由很多不同尺度的浪叠加而成的, 并且导航X波段雷达的图像中还有降雨和系统噪声的影响(例如图 3b中的低频噪声)。由于雷达的空间分辨率是2.4 m, 观测范围是3 km, 我们主要观测的波浪是波长在几十米至几百米的重力波。本研究中只保留波长为30~200 m的波浪, 而将雷达信号中的其他成分滤除。

为此, 选用Butterworth带通滤波器^[14]。通带的波数为0.03~0.21 rad/m, 阻带的波数为0.02~0.31 rad/m。

2.4 波数谱

对处理后的雷达信号做一维傅里叶变换, 得到其一维波数谱:

$S(k) = \int_{\; - \infty }^{\;\infty } {I(x) \cdot {{\rm{e}}^{ - ikx}}} {\rm{d}}x$

(3)

式中, I(x)为雷达回波强度的第一主成分经过滤波后的值, x为像素点的坐标, k为波数。

由于雷达图像没有经过定标, 方程(3)的谱能量是相对值, 并不是真实的能量。但是, 根据该波数谱(3)的峰值位置, 可以确定海浪的峰值波数k_peak。

根据海浪的频散关系:

${\omega ^{\rm{2}}} = gk\tanh \left( {kd} \right)$

(4)

式中, d是海浪传播区域的水深, $\omega = \frac{{2{\rm{ \mathsf{ π} }}}}{T}$是角频率, T是海浪的周期, g是重力加速度, k是波数。在深水情况下, 方程(4)可以简化为

${\omega ^{\rm{2}}} = gk$

(5)

此时波长和波速与波周期有关, 而与水深无关。根据方程(5)和峰值波数, 可以得到海浪的峰值周期T_peak。

2.5 海浪谱的数值仿真

JONSWAP谱是一种成熟的半经验海浪谱^[15], 适用于风浪和涌浪等多种海况。因此, 本文选用JONSWAP谱作为标准的海浪谱, 假定其适用于观测海区。JONSWAP谱的形式如公式(6)所示。

$S\left( \omega \right) = \alpha {g^2}{\omega ^{ - 5}}\exp [-1.25{({\omega _p}/\omega )^4}]{\gamma ^{\exp [-{{(\omega-{\omega _p})}^2}/(2{\sigma ^2}{\omega _p}^2)]}}$

(6)

式中, g为重力加速度; ω_p为谱峰频率; α为能量尺度参数; γ为谱升因子(γ的观测值介于1.5~6.0之间, 本文取其平均值3.3); σ为峰形参数, 其值等于

$\left. \begin{array}{l} \sigma = 0.07, \omega \le {\omega _p}\\ \sigma = 0.09, \omega > {\omega _p} \end{array} \right\}$

(7)

无因次常数α=0.081 7·F^–0.286, 无因次风区$F = (g \cdot f)/u_{10}^2$(f为风区, u₁₀为10 m高度处的风速), 其值范围在10^–1~10⁵, 无因次峰频率ω_p= 22(g/u₁₀)F^–0.33。

在风区、峰形参数、谱升因子等参数一定时, 输入风速就可以确定公式(6)。由于不同形状的谱的峰值周期不同, 可以建立以最小化经验谱与反演谱的峰值周期为目标的约束函数

$\Delta = {T_{\rm{sp}}} - {T_{\rm{peak}}}$

(8)

式中, T_sp为仿真海浪谱的峰值周期; T_peak为反演的雷达图像谱的峰值周期。

本文中的海浪谱仿真只与2个参数有关: (1)10 m风速u₁₀; (2)峰值周期T_peak。求解海浪谱的方法如下: (1)输入不同的风速值, 迭代求解方程(8), 找出使目标函数最小的风速值, 即:首先设置风速的范围是1~25m/s, 以0.2 m/s为步长选取不同的风速, 得到与之对应的ω_p, 从而得到不同风速时的峰值周期T_sp, 根据公式(8)确定最优风速。(2)将求得的风速值代入公式(6), 作为真实的海浪谱。

2.6 有效波高的反演

确定海浪谱(公式(6))后, 根据海浪理论可以得到有效波高。对海浪谱S(ω)关于ω积分, 可得海面高度的方差为:

$\left\langle {{\varsigma ^{\rm{2}}}} \right\rangle = \int_{\;{\rm{0}}}^{\;\infty } {S\left( \omega \right){\rm{d}}\omega } \;$

(9)

从而可以得到海浪的有效波高:

${H_{{\rm{1/3}}}} = 4{\left\langle {{\varsigma ^{\rm{2}}}} \right\rangle ^{{\rm{1/2}}}}$

(10)

综上所述, 本文提出了一种新的利用降雨情况下的导航X波段雷达图像反演海浪的有效波高的算法, 流程图如图 5所示。

3 结果分析

为了验证提出的算法, 下面首先以图 2为例进一步说明算法的步骤, 然后用浮标测量的有效波高对算法进行了验证。

3.1 数据处理

以图 2为例进一步说明算法的基本步骤。选取原始雷达图像中相邻的11个(图 2中的白色虚线附近)径向灰度值, 作主成分分析, 如图 6。根据图 6a和图 6b, 第一主成分包含了明显的海浪变化, 而其他主成分的噪声较大(振幅在0附近), 所以本文只考虑第一主成分是合理的。与图 3a和图 4a相比, 主成分的波峰和波谷更加清晰, 即能更好地反映波浪的变化。

然后对图 6中的主成分做带通滤波, 如图 7, 可以看出滤波后的数据有效消除了原数据的低频和高频分量, 而保留了需要的海浪信息。图 8是对图 7作一维傅里叶变换后的波数谱, 无降雨时谱的峰值对应的波长为112 m(图 8a)、降雨时谱的峰值波长为102.4 m(图 8b), 与真实的海浪符合较好。这说明本文提出的方法可以有效从降雨时的雷达图像中提取海浪的峰值信息。

3.2 有效波高的反演结果

利用第2部分的算法, 分别利用每组导航X波段雷达图像反演海浪的有效波高, 然后与同步观测的浮标测量的有效波高做比较, 结果如图 9。图 9表明反演的有效波高与浮标测量的有效波高有较好的一致性, 二者的均方根误差为0.23 m、相关系数为0.57。这说明本文提出的利用降雨情况下的雷达图像估计有效波高的算法是合理的。

4 讨论与展望

降雨等恶劣天气条件下的海浪观测有重要应用, 但是降雨对导航X波段雷达信号的衰减影响较大, 导致不能正常观测海浪参数, 而且实际应用时一般没有浮标等观测数据, 只能利用雷达图像。本文提出了一种新的利用降雨时的雷达图像反演海面的有效波高与海浪谱的方法, 与经典的波高反演算法^[10]相比, 本文算法的一个优点就是不需要利用浮标等辅助数据做参考, 独立应用雷达图像反演海浪的波高。通过对比实验数据和雷达反演数据, 对反演方法进行了验证, 有效波高反演的均方根误差为0.23 m、相关系数为0.57, 证明了该方法的可行性。

该方法也存在一些不足之处。第一, 当降雨量过大时, 雷达图像中几乎没有海浪条纹, 无法用主成分分析提取波峰和波谷信息, 而实验中没有实测降雨量, 所以本文没有定量考虑降雨的影响; 第二, 本实验中观测的波浪适用于JONSWAP谱, 但是其他海区的海浪谱可能与JONSWAP谱有偏差, 算法的适用性有待更多验证; 第三, 主成分分析法在对数据的处理过程中只考虑了图像数据的二阶统计信息, 未能利用数据中的高阶统计信息, 忽略了多个像素间的非线性相关性, 所以变换后的数据间仍有可能存在高阶冗余信息。

在未来的工作中, 将继续寻找有效方法提高有效波高反演的精度, 尤其是浅海情况下的反演。另外, 需要采用多地点长时间序列的观测数据对方法进行进一步验证, 希望能将其应用于业务化的波浪观测中。