基于航海微波雷达的波浪重建及预测方法

魏鋆宇; 高思哲; 谭啸; 郭鹏程; 尹子明; 王阔

doi:10.12404/j.issn.1671-1815.2506550

科学技术与工程 ›› 2026, Vol. 26 ›› Issue (19) :8233-8242 DOI: 10.12404/j.issn.1671-1815.2506550

电子技术、通信技术

基于航海微波雷达的波浪重建及预测方法

作者信息 +

Wave Reconstruction and Prediction Method Based on Marine Microwave Radar

Author information +

文章历史 +

摘要

波浪重建及预测旨在实现海洋表面瞬时波高的还原和预报,该方法能够为船舶航行提供前瞻性信息,降低海上作业的风险。航海微波雷达是目前用于精确波浪测量的新兴工具,可用于推进波浪重建及预测的发展。提出了一种基于正则化最小二乘的波浪重建及预测算法。首先根据航海微波雷达获取的水质点径向速度空时矩阵和传播算子矩阵,利用正则化最小二乘法估计波浪模型参数矩阵,正则化参数的引入保证了参数求解的稳定性;然后根据线性波浪传播模型,重建或预测不同空时位置的海面瞬时波高。在数值仿真实验中,分析了波浪可预测区的分布规律,并评估了该方法在不同海况下的重建及预测性能,结果显示在可预测区内参考波高与预测波高吻合程度很高,证明了所提方法具有良好的波浪重建及预测性能。

Abstract

Wave reconstruction and prediction aims to realize the reduction and forecasting of transient wave heights on the ocean surface, and the method can provide forward-looking information for ship navigation and reduce the risk of offshore operations. Marine microwave radar is an emerging tool currently used for accurate wave measurements that can be used to advance wave reconstruction and prediction. The wave reconstruction and prediction algorithm based on regularized least squares was proposed. Firstly, the wave model parameter matrix was estimated using the regularized least squares method based on the radial velocity space-time matrix and the propagation operator matrix of the water particle acquired by the marine microwave radar, and the introduction of the regularization parameter ensures the stability of the parameter solution. Then, based on the linear wave propagation model, the instantaneous wave heights at the sea surface were reconstructed or predicted for different airtime positions. In the numerical simulation experiments, the distribution law of the wave predictable zone was analyzed, and the reconstruction and prediction performance of the proposed method was evaluated under different sea conditions, and the results show that the reference wave height matches the predicted wave height very well within the predictable zone, which proves that the proposed method has a good wave reconstruction and prediction performance.

Graphical abstract

关键词

波浪重建及预测 / 航海微波雷达 / 正则化最小二乘 / 线性波浪传播模型 / 波浪可预测区

Key words

wave reconstruction and prediction / marine microwave radar / regularized least squares / linear wave propagation model / wave predictable zone

引用本文

引用格式 ▾

魏鋆宇, 高思哲, 谭啸, 郭鹏程, 尹子明, 王阔. 基于航海微波雷达的波浪重建及预测方法[J]. 科学技术与工程, 2026, 26(19): 8233-8242 DOI:10.12404/j.issn.1671-1815.2506550

WEI Yun-yu, GAO Si-zhe, TAN Xiao, GUO Peng-cheng, YIN Zi-ming, WANG Kuo. Wave Reconstruction and Prediction Method Based on Marine Microwave Radar[J]. Science Technology and Engineering, 2026, 26(19): 8233-8242 DOI:10.12404/j.issn.1671-1815.2506550

登录浏览全文

4963

注册一个新账户忘记密码

波浪重建及预测是一种根据海浪几何模型对不同空时位置的海浪瞬时波高进行估计和预测的方法^[1]。近年来,航海微波雷达依靠其较高的时空分辨率逐渐成为重要的海浪监测工具,其回波信号可作为波浪重建和预测的初始条件,实现对当前或未来海浪演变过程的准确监测和预测^[2-3]。因此,研究基于航海微波雷达的波浪重建及预测方法,对于海洋灾害预警、保障海上作业安全等方面具有重要意义。

根据雷达体制不同,可将航海微波雷达分为相参微波雷达和非相参微波雷达两大类。非相参微波雷达通常通过对雷达照射区域的反向散射图像应用三维快速傅里叶变换(fast Fourier transform,FFT)来获取海浪信息。Wright^[4]使用X波段航海雷达的回波图像对海浪的起伏特征进行了研究。之后,Keller等^[5]提出了复合表面散射理论,该理论为航海微波雷达的波浪探测研究奠定了理论基础。随后也出现了较为成熟的非相参微波雷达系统,例如WAVEX系统、WaMoS Ⅱ系统和OWME系统^[6-10]。不过,非相参微波雷达回波信号处理精度依赖于调制转换函数(modulation transfer function,MTF)的准确获取^[6-8], 因此非相参微波雷达的海浪探测精度、稳定性和置信度可能难以达到要求。与非相参微波雷达不同,相参微波雷达可以分辨回波信号的相位,从而能够获得海面的多普勒信息,因此相参微波雷达可以测得海面水质点的速度信息。Plant 等^[11]拉开了利用相参雷达探测海浪的序幕。随后不同国家的学者都针对相参微波雷达的海洋观测领域展开了研究。武汉大学无线电海洋遥感实验室于2011年独立研发了一款高精度的S波段海洋遥测系统MORSE(microwave ocean remote sensor),并进行了多次海上比测实验^[3,12]。相参微波雷达无需校正的优势,使得其在未来将被用于更多的海浪探测研究中,因此,针对相参微波雷达开展波浪重建及预测方法的研究具有重大意义。

波浪重建是对雷达探测的空时域的回波信号的瞬时波高进行估计。由线性波浪理论可知,海面水质点径向速度与海浪模型参数之间存在定量关系,根据水质点径向速度空时分布即可实现波浪重建。文献[13-14]提出势函数法进行波浪重建,并在相参接收的准相参微波雷达实验中得到初步验证,但该方法由于没有求解海浪模型参数,无法实现波浪预测。随后,Connell等^[15]提出最小二乘法,用于估计海浪模型参数,再根据海面的几何模型叠加实现波浪重建及预测。Støle-Hentschel等^[16]基于相参X波段雷达重建海面,并将重建海面的有效波高、峰值波周期和浮标数据进行对比,该方法间接证明了波浪重建的准确性。

波浪预测是在波浪重建的基础上对未来时间和雷达探测区域外瞬时波高的预测分析。由于受到波浪传播过程中色散特性的影响,随着波浪的传播和演化,观测区内的波浪成分会在越来越小的时空范围内重叠,通常将该区域称为可预测区。Morris等^[17]定点测量长峰波波浪后,提出在一定时空范围内存在一个区域可实现对波浪的准确重建,这是对理论可预测区域概念的首次阐述。Wu等^[18]提出了波浪预测的可预测区域大小取决于波浪群速度的观点,该观点由多次数值仿真实验得到,并被Naaijen等^[19]验证。上述研究均通过接触式点测量方法实现,针对遥感设备的理论可预测区研究也在随后得到发展。基于线性波浪理论,Kosleck^[20]分析和仿真了空间采样下波浪的理论可预测区域,并且在水池实验中得到了验证。在中国,马学文等^[21]基于线性波浪理论,通过哈尔滨工程大学的水池实验对波浪场的理论可预测区进行了验证。目前,针对相参微波雷达的波浪重建及预测方法研究受到国内外学者的广泛关注,且尚处于快速发展阶段。探索一种适用于航海微波雷达的重建及预测方法,对海面波浪信息的准确度量、海上作业安全、海洋能源探测、舰船耐波性设计等领域均具有重要意义。

现提出基于正则化最小二乘的波浪重建及预测方法,该方法首先分别沿时间维度和距离维度提取水质点径向速度,得到水质点径向速度空时矩阵,并基于线性波浪理论构建传播算子矩阵。然后根据L曲线法获取最优的正则化参数,利用正则化最小二乘法求解波浪模型参数矩阵。接着依据相参微波雷达的特点,将海面近似成多个带宽极窄的分量叠加,并计算所有分量的交集,即可划分波浪的可预测区。在可预测区内,将模型参数代入极坐标下线性波浪模型中,可实现不同空时位置的海面高度重建或预测。利用数值仿真实验验证所提方法在不同海况下的波浪重建及预测性能。

1 方法

1.1 线性波浪理论

海面的几何模型是根据线性波浪理论生成的,通常利用已知海浪方向谱推导产生的波幅,然后叠加多个振幅、频率、相位和传播方向不同的正弦波,从而得到海面上某一时空位置的水面高度,即Longuet-Higgins模型^[22-23]。

设t时刻在位置(x,y)处的波面方程为η(x,y,t),其波面模型可表示为

(1)η(x,y,t)=$\stackrel{\mathit{M}}{\sum _{\mathit{m}=1}}\stackrel{\mathit{N}}{\sum _{\mathit{n}=1}}$A_m_,_ncos(k_nxcosθ_m+k_nysinθ_m-ω_nt+ε_mn)

式(1)中:η(x,y,t)为t时刻平面(x,y)处的浪高;A_m_,_n为角频率为ω_n且传播方向与x轴夹角为θ_m的组成波的振幅;k_n为波数;ε_mn为初始相位,理论上符合[0,2π)区间的均匀分布;M为对波传播方向的离散取样数量;N为对波频率(或波数)的离散取样数量。根据深水条件下的频散关系^[24]可知,

${\mathit{\omega }}_{\mathit{n}}^{2}$

=k_ng(g为重力加速度)。

航海相参微波雷达以掠入射方式照射海面,假设其位置为极坐标原点,r为雷达照射点到海面观测点的径向距离,ϕ为照射方位角,则极坐标下波面模型可表示为

(2)η(r,ϕ,t)=$\stackrel{\mathit{M}}{\sum _{\mathit{m}=1}}\stackrel{\mathit{N}}{\sum _{\mathit{n}=1}}$A_m_,_ncos[k_nrcos(ϕ-θ_m)-ω_nt+ε_mn]

深海中海面上水质点运动的轨迹是以波高为直径的圆,当水质点运动到最高位置(即波峰位置)时,其运动方向与波向一致,当水质点运动到最低位置(即波谷位置)时,其运动方向与波向相反^[25]。

令模型参数对(a_mn,b_mn)=(A_mncosε_mn,A_mnsinε_mn),且相位项Ω_mn=k_nrcos(ϕ-θ_m)-ω_nt。海面水质点运动水平速度v_h和垂直速度v_v可表示为

(3)

$\left\{\begin{array}{l}{\mathit{v}}_{\mathrm{h}}(\mathit{r},\mathit{\varphi },\mathit{t})=\stackrel{\mathit{M}}{\sum _{\mathit{m}=1}}\stackrel{\mathit{N}}{\sum _{\mathit{n}=1}}{\mathit{\omega }}_{\mathit{n}}{\mathit{A}}_{\mathit{m},\mathit{n}}\mathrm{c}\mathrm{o}\mathrm{s}({\mathit{\Omega }}_{\mathit{m}\mathit{n}}+{\mathit{\epsilon }}_{\mathit{m}\mathit{n}})\\ {\mathit{v}}_{\mathrm{v}}(\mathit{r},\mathit{\varphi },\mathit{t})=\stackrel{\mathit{M}}{\sum _{\mathit{m}=1}}\stackrel{\mathit{N}}{\sum _{\mathit{n}=1}}{\mathit{\omega }}_{\mathit{n}}{\mathit{A}}_{\mathit{m},\mathit{n}}\mathrm{s}\mathrm{i}\mathrm{n}({\mathit{\Omega }}_{\mathit{m}\mathit{n}}+{\mathit{\epsilon }}_{\mathit{m}\mathit{n}})\end{array}\right.$

航海相参微波雷达的掠入射角为γ,则水质点径向速度为

(4)v(r,ϕ,t)=v_h(r,ϕ,t)cos(ϕ-θ_m)cosγ+v_v(r,ϕ,t)sinγ

令

(5)

$\left\{\begin{array}{l}{\mathit{\alpha }}_{\mathit{m}\mathit{n}}(\mathit{r},\mathit{\varphi },\mathit{t})={\mathit{\omega }}_{\mathit{n}}\mathrm{c}\mathrm{o}\mathrm{s}(\mathit{\varphi }-{\mathit{\theta }}_{\mathit{m}})\mathrm{c}\mathrm{o}\mathrm{s}\mathit{\gamma }\mathrm{c}\mathrm{o}\mathrm{s}{\mathit{\Omega }}_{\mathit{m}\mathit{n}}+\\ {\mathit{\omega }}_{\mathit{n}}\mathrm{s}\mathrm{i}\mathrm{n}\mathit{\gamma }\mathrm{s}\mathrm{i}\mathrm{n}{\mathit{\Omega }}_{\mathit{m}\mathit{n}}\\ {\mathit{\beta }}_{\mathit{m}\mathit{n}}(\mathit{r},\mathit{\varphi },\mathit{t})={\mathit{\omega }}_{\mathit{n}}\mathrm{c}\mathrm{o}\mathrm{s}(\mathit{\varphi }-{\mathit{\theta }}_{\mathit{m}})\mathrm{c}\mathrm{o}\mathrm{s}\mathit{\gamma }\mathrm{s}\mathrm{i}\mathrm{n}{\mathit{\Omega }}_{\mathit{m}\mathit{n}}-\\ {\mathit{\omega }}_{\mathit{n}}\mathrm{s}\mathrm{i}\mathrm{n}\mathit{\gamma }\mathrm{c}\mathrm{o}\mathrm{s}{\mathit{\Omega }}_{\mathit{m}\mathit{n}}\end{array}\right.$

则

(6)v(r,ϕ,t)=$\stackrel{\mathit{M}}{\sum _{\mathit{m}=1}}\stackrel{\mathit{N}}{\sum _{\mathit{n}=1}}$[a_mnα_mn(r,ϕ,t)+b_mnβ_mn(r,ϕ,t)]

1.2 波浪重建

航海相参微波雷达可探测一定范围内的水质点径向速度,若浪向已知,由水质点径向速度的空时序列,结合线性海面模型,求解模型参数a_mn和b_mn即可预测海面任意时空位置的波高。

航海相参微波雷达获得的径向速度空时矩阵V可以表示为矩阵乘法,即

(7)PX=V

(8)

$\left\{\begin{array}{l}\mathit{P}=\left[\begin{array}{ll}{\mathrm{\alpha }}_{\mathrm{m}\mathrm{n}}& {\mathrm{\beta }}_{\mathrm{n}\mathrm{m}}\end{array}\right]\\ \mathit{X}=\left[\begin{array}{l}{\mathrm{a}}_{\mathrm{m}\mathrm{n}}\\ {\mathrm{b}}_{\mathrm{m}\mathrm{n}}\end{array}\right]\end{array}\right.$

式(7)中:P为传播算子矩阵;X为待求解的参数矩阵;V为雷达在所有方位单元和距离单元观测到的径向速度空时矩阵。

由于雷达的探测距离元和照射方位角为离散值。若将雷达测得的某一径向距离r_i,照射方位角ϕ_j处的水质点径向速度记为v_e(r_i,ϕ_j,t_j),i∈[1,L],j∈[1,S],L和S分别为距离元和照射方位角的个数。LS≥2NM,N和M为估计时模型中考虑的波数分量个数和方向个数。则上述矩阵乘法可写为

(9)[$\begin{array}{ll}{\mathit{\alpha }}_{\mathit{n}\mathit{m}}({\mathit{r}}_{\mathit{i}},{\mathit{\varphi }}_{\mathit{j}},{\mathit{t}}_{\mathit{j}})& {\mathit{\beta }}_{\mathit{n}\mathit{m}}({\mathit{r}}_{\mathit{i}},{\mathit{\varphi }}_{\mathit{j}},{\mathit{t}}_{\mathit{j}})\end{array}$]$\left[\begin{array}{l}{\mathit{a}}_{\mathit{n}\mathit{m}}\\ {\mathit{b}}_{\mathit{n}\mathit{m}}\end{array}\right]$= $\left[\begin{array}{c}\begin{array}{l}\begin{array}{l}{\mathit{v}}_{\mathrm{e}}({\mathit{r}}_{1},{\mathit{\varphi }}_{1},{\mathit{t}}_{1})\\ {\mathit{v}}_{\mathrm{e}}({\mathit{r}}_{2},{\mathrm{\varphi }}_{2}^{},{\mathrm{t}}_{2}^{})\end{array}\\ \mathit{　}\mathit{　}\mathit{︙}\\ {\mathit{v}}_{\mathrm{e}}({\mathit{r}}_{\mathit{i}},{\mathit{\varphi }}_{\mathit{j}},{\mathit{t}}_{\mathit{j}})\end{array}\end{array}\right]$

该矩阵方程通常使用最小二乘法求解目标函数,即

(10)$X_{\mathrm{LS}}=\min _{x}\left(\|\boldsymbol{P} \boldsymbol{X}-\boldsymbol{V}\|_{2}^{2}\right)=\left(\boldsymbol{P}^{\mathrm{T}} \boldsymbol{P}\right)^{-1} \boldsymbol{P}^{\mathrm{T}} \boldsymbol{V}$

式(10)中:X_LS为常规的最小二乘解;$ \|\cdot\|^{2}$为欧几里得范数。虽然最小二乘法具有简单易行、计算量小、适用面广等优点,但是在相参微波雷达海浪探测中,受到高海况或硬目标杂波的影响,速度测量值中往往会存在许多异常值。最小二乘法对异常值较为敏感,因此利用最小二乘法求解会导致波浪预测的性能严重下降^[26-27]。

基于上述分析,用Tikhonov正则化求解PX=V中的目标函数,即

(11)$\boldsymbol{X}_{\lambda}=\min _{x}\left(\|\boldsymbol{P} \boldsymbol{X}-\boldsymbol{V}\|_{2}^{2}+\lambda^{2}\|x\|_{2}^{2}\right)$

式(11)中:X_λ为正则化最小二乘解;λ为正则化参数,目标函数也可表示为

(12)X_λ=(P^TP+λI)^-1P^TV

式(12)中:I为单位矩阵,相当于在矩阵P^TP的每个对角元素上加入一个扰动,从而使求逆的矩阵改变,降低条件数,使得最小二乘解更稳定。当λ=0时,正则化的方法退化为常规最小二乘。

正则化参数λ的最优值,可通过L曲线法求解。L曲线法实现简单,结果稳健,其通过结合图形特征来选择最优正则化参数。根据L曲线法原理,对于正则化参数的选择,需要在使残差和解的范数同时达到最小化之间进行权衡^[28-29]。

残差$\|P X-V\|_{2}$和解的范数$\left\|X_{\lambda}\right\|_{2}$的定义为

(13)$\zeta(\lambda)=\|P X-V\|_{2}$

(14)$\rho(\lambda)=\left\|X_{\lambda}\right\|_{2}$

Hansen等^[30]指出,由残差和解的范数构成的曲线γ=[lgζ(λ),lgρ(λ)]将会近似呈“L”形,最优的正则化参数λ则对应L曲线的拐角,如图1所示。

令

$\stackrel{\wedge }{\mathit{\zeta }}$

(λ)=lgζ(λ),

$\stackrel{\wedge }{\mathit{\rho }}$

(λ)=lgρ(λ),通过计算L曲线最大曲率可以确定L曲线的拐角位置。曲率κ(λ)可表示为

(15)κ(λ)=$\frac{\hat{\mathit{\rho }}\mathit{\text{'}}\hat{\mathit{\zeta }}″-\hat{\mathit{\rho }}″\hat{\mathit{\zeta }}\mathit{\text{'}}}{(\hat{\mathit{\rho }}{\mathit{\text{'}}}^{2}+\hat{\mathit{\zeta }}{\mathit{\text{'}}}^{2}{)}^{3/2}}$

1.3 可预测区划定

鉴于波浪传播具备色散特性,从航海相参微波雷达观测区传出的波浪成分会随着时间的流逝,在愈发狭小的时空范围内相互重叠。所以,借助观测区的回波信号来估计海浪模型参数并进行预测,仅能在有限的时空范围内实现对波浪的有效预测,这个区域即为可预测区。可预测区边界的确定是波浪精准预测的先决条件。

使用波数域上的积分表示二维波浪场任意时空位置(r,t)的连续海面高度,即:

(16)η(r,ϕ₀,t)=Re[${\int }_{{\mathit{k}}_{\mathrm{l}}}^{{\mathit{k}}_{\mathrm{h}}}$A_c(k)e^j(^kr^-^ωt⁾dk],r,t∈(-∞,+∞)

式(16)中:k_l和k_h分别为最小和最大的波数分量;ω为角频率;波数k对应的复振幅为A_c(k)。

若将海面近似成N_k个带宽极窄的分量叠加,对应的窄带分量为k_n。对于方向ϕ=ϕ₀,离散状态下海面高度表示为

(17)η(r,ϕ₀,t)=Re$\left[\stackrel{{\mathit{N}}_{\mathit{k}}}{\sum _{\mathit{n}=1}}{\int }_{{\mathit{k}}_{\mathrm{n}}-\frac{\mathrm{\Delta }\mathit{k}}{2}}^{{\mathit{k}}_{\mathrm{n}}+\frac{\mathrm{\Delta }\mathit{k}}{2}}{\mathit{A}}_{\mathrm{c}}\left(\mathit{k}\right){\mathrm{e}}^{\mathrm{j}(\mathit{k}\mathit{r}-\mathit{\omega }\mathit{t})}\mathrm{d}\mathit{k}\right]$

式(17)中:Δk=

$\frac{{\mathit{k}}_{\mathrm{h}}-{\mathit{k}}_{\mathrm{l}}}{{\mathit{N}}_{\mathit{k}}}$

,窄带分量k_n=k_l+(n-1/2)Δk。

假设在t₀时刻,观测空间范围r∈[R₀,R₁]的海面,C_gn代表波数分量k_n对应的群速,由深水条件下的色散关系可知:C_gn=

$\frac{\partial {\mathit{\omega }}_{\mathit{n}}}{\partial {\mathit{k}}_{\mathrm{n}}}$

$\frac{\mathit{\omega }}{2\mathit{k}}$

$\frac{\mathit{g}}{2\mathit{\omega }}$

(g为重力加速度)。

k_n对应的可预测区在空时域中呈现为一个长条带,如图2所示,群速度C_gn即为该区域斜率的倒数。

二维波浪场的可预测区如图3所示。图中蓝色区域表示可预测区,波浪重建也就是t₀时刻的预测, 对于时间小于t₀时刻的预测也被称为“后报”。C_g_l、C_g_h为观测区波浪的最大、最小群速,对应的波数为k_l、k_h。雷达有效探测距离为R=R₁-R₀。

2 重建及预测流程

根据上述理论分析,基于航海微波雷达的波浪重建及预测方法流程如图4所示,步骤如下。

(1)根据线性波浪模型,模拟不同海况下波浪的水质点径向速度及瞬时波高(作为参考瞬时波高)空时分布。

(2)分别沿时间维度和距离维度提取水质点径向速度,并构建传播算子矩阵。

水质点径向速度空时矩阵由雷达探测的部分方位单元和距离单元观测到的径向速度构成。假设雷达探测的距离元和照射方位角的个数分别为L和S,则某一径向距离r_i,照射方位角ϕ_j处的水质点径向速度记为

(18)v_e(r_i,ϕ_j,t_j), i∈[1,L],j∈[1,S]

因此,水质点径向速度空时矩阵可表示为

(19)V=$\left[\begin{array}{l}\begin{array}{l}{\mathit{v}}_{\mathrm{e}}({\mathit{r}}_{1},{\mathit{\varphi }}_{1},{\mathit{t}}_{1})\\ {\mathit{v}}_{\mathrm{e}}({\mathit{r}}_{2},{\mathit{\varphi }}_{2},{\mathit{t}}_{2})\end{array}\\ \mathit{　}\mathit{　}\mathit{︙}\\ {\mathit{v}}_{\mathrm{e}}({\mathit{r}}_{\mathit{i}},{\mathit{\varphi }}_{\mathit{j}},{\mathit{t}}_{\mathit{j}})\end{array}\right]$

通过线性波浪理论,可以构建传播算子矩阵P=[α_nm β_nm]。某一径向距离r_i,照射方位角ϕ_j处的传播算子矩阵可以表示为

(20)P=$\left[\begin{array}{llll}{\mathit{\alpha }}_{\mathit{n}\mathit{m}}({\mathit{r}}_{1},{\mathit{\varphi }}_{1},{\mathit{t}}_{1})& & {\mathit{\beta }}_{\mathit{n}\mathit{m}}({\mathit{r}}_{1},{\mathit{\varphi }}_{1},{\mathit{t}}_{1})& \\ {\mathit{\alpha }}_{\mathit{n}\mathit{m}}({\mathit{r}}_{2},{\mathrm{\varphi }}_{2}^{},{\mathit{t}}_{2})& & {\mathit{\beta }}_{\mathit{n}\mathit{m}}({\mathit{r}}_{2},{\mathit{\varphi }}_{2},{\mathit{t}}_{2})& \\ \dots & & \dots & \\ {\mathit{\alpha }}_{\mathit{n}\mathit{m}}({\mathit{r}}_{\mathit{i}},{\mathit{\varphi }}_{\mathit{j}},{\mathit{t}}_{\mathit{j}})& & {\mathit{\beta }}_{\mathit{n}\mathit{m}}({\mathit{r}}_{\mathit{i}},{\mathit{\varphi }}_{\mathit{j}},{\mathit{t}}_{\mathit{j}})& \end{array}\right]$

(3)通过L曲线法求解最优的正则化参数λ,对于二维长波峰波浪,令λ=0。

L曲线法通过选取最优正则化参数,可以平衡残差和解的范数的最小值。

残差$\|P X-V\|_{2}$和解的范数$\left\|X_{\lambda}\right\|_{2}$的定义为

(21)$\zeta(\lambda)=\|P X-V\|_{2}$

(22)$\rho(\lambda)=\left\|X_{\lambda}\right\|_{2}$

令

$\hat{\mathrm{\zeta }}$

(λ)=lgζ(λ),

$\hat{\mathrm{\rho }}$

(λ)=lgρ(λ),L曲线最大曲率位置即为L曲线的拐角位置。曲率κ(λ)可表示为

(23)κ(λ)=$\frac{\hat{\mathrm{\rho }}\mathit{\text{'}}\hat{\mathrm{\zeta }}″-\hat{\mathrm{\rho }}″\hat{\mathrm{\zeta }}\mathit{\text{'}}}{(\hat{\mathrm{\rho }}{\mathit{\text{'}}}^{2}+\hat{\mathrm{\zeta }}{\mathit{\text{'}}}^{2}{)}^{3/2}}$

(4)利用正则化最小二乘法,求解PX=V矩阵方程中的模型参数矩阵X=

$\left[\begin{array}{l}{\mathit{a}}_{\mathit{n}\mathit{m}}\\ {\mathit{b}}_{\mathit{n}\mathit{m}}\end{array}\right]$

。

基于正则化最小二乘法求解PX=V中的目标函数如式(10)所示,通过步骤(3)中L曲线法可获取正则化参数,即化简后的目标函数如式(11)所示。

(5)根据航海相参微波雷达的特点,划分波浪的可预测区。

当时间为t₀时,雷达观测海面空间范围为r∈[R₀,R₁]。若将海面近似成N_k个带宽极窄的分量叠加,对应的窄带分量为k_n。

对于窄带分量k_n,探测位置为

(24)X=r-C_gn(t-t₀), R₀≤X≤R₁

即

(25)R₀+C_gn(t-t₀)≤r≤R₁+C_gn(t-t₀)

波数k∈[k_l,k_h]的波浪对应的可预测区为所有分量可预测区的交集,即

(26)∏(t)=${\cap }_{\mathit{k}\in [{\mathit{k}}_{\mathrm{l}},{\mathit{k}}_{\mathrm{h}}]}$∏_k(t)

式(26)中:∏_k(t)为R₀+C_gk(t-t₀)≤r≤R₁+C_gk(t-t₀)。

(6)在可预测区内,可实现不同空时位置的海面瞬时波高的重建或预测。

模型参数矩阵为X=

$\left[\begin{array}{l}{\mathit{a}}_{\mathit{n}\mathit{m}}\\ {\mathit{b}}_{\mathit{n}\mathit{m}}\end{array}\right]$

,其中

$\left\{\begin{array}{l}{\mathit{a}}_{\mathit{m}\mathit{n}}={\mathit{A}}_{\mathit{m}\mathit{n}}\mathrm{c}\mathrm{o}\mathrm{s}{\mathit{\epsilon }}_{\mathit{m}\mathit{n}}\\ {\mathit{b}}_{\mathit{m}\mathit{n}}=-{\mathit{A}}_{\mathit{m}\mathit{n}}\mathrm{s}\mathrm{i}\mathrm{n}{\mathit{\epsilon }}_{\mathit{m}\mathit{n}}\end{array}\right.$

。将式(10)代入极坐标下线性波浪模型[式(2)]中,得

(27)η(r,ϕ,t)=$\stackrel{\mathit{M}}{\sum _{\mathit{m}=1}}\stackrel{\mathit{N}}{\sum _{\mathit{n}=1}}$A_m_,_ncos[k_nrcos(ϕ-θ_m)-ω_nt+ε_mn]

可实现在可预测区内,对海面上不同空时位置的瞬时波高的重建和预测。将该结果与步骤(1)中对应空时位置的参考瞬时波高进行对比,可分析波浪重建或预测的准确性。

3 结果分析

数值仿真实验中,将参考波高与重建及预测波高进行比较,用于验证波浪重建及预测方法的性能。对于短峰波波浪,波浪的振幅由海浪谱S(k)和方向分布函数G(ϕ)共同确定。模拟中,主波方向为90°,海面上方10 m处风速分别设置为8 m/s和15 m/s。假设雷达位于坐标原点,掠入射角度为5°,距离分辨率为7.5 m,雷达照射海面的径向距离覆盖范围为750 m,方位角角度分辨率为1°,雷达探测覆盖360°的范围。

为了便于确定最大、最小群速以及计算可预测区,海浪谱S(k)的最大、最小波数分量确定依据为

(28)S(k_min)=S(k_max)=δS(k_p)

式(28)中:k_p为谱峰波数值;δ为远小于1的值;k_min和k_max分别为最小、最大波数。

根据线性波浪模型和上述参数,模拟结果如图5所示。可以看出,不同风速下的海浪纹理方向存在明显差异,且海浪纹理反映的浪向与仿真设置数值大致相同。对比图5可以看出,当风速变大时,海况升高,波高和波周期都明显变大,这与实际海洋环境相符。

3.1 重建

在波浪重建研究中,常见的方法有势函数法和最小二乘法,其中势函数只能用来重建雷达探测区域的波浪场,而最小二乘法估计了海浪模型各分量的参数,结合线性波浪模型可实现对波浪的重建及预测。

势函数方法假设航海微波雷达在掠入射的条件下,其所测的水质点径向速度近似为水平速度,因此通过对水平速度在距离上积分可以得到速度势函数,由势函数在时间上求导实现波高重建。具体计算过程如下。

速度势ϕ_v(r,ϕ,t)可以通过对径向速度场v_r(r,ϕ,t)在径向距离r上积分得到,即

(29)ϕ_v(r,ϕ,t)=${\int }_{{\mathit{r}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{\mathit{r}}$v_r(r',ϕ,t)dr'

式(29)中:r_min为最小径向距离。

由小振幅波理论可知,动力学边界条件存在以下关系:

(30)gη+$\frac{\partial {\mathit{\varphi }}_{\mathit{v}}}{\partial \mathit{t}}$=0

因此,对势函数ϕ_v在时间上求导,实现波高η的重建,即

(31)η(r,ϕ,t)=-$\frac{1}{\mathit{g}}\frac{\partial {\mathit{\varphi }}_{\mathit{v}}}{\partial \mathit{t}}$

这种方法需要借助雷达同一方位角在短时间内的两次回波信号,以差分的形式实现前述求导,且差分时间越短越理想。虽然该方法可以提取出瞬时海面波高,但是无法计算海浪各分量的解相位模型参数,因此,势函数法无法实现波浪预测这一目标。

数值仿真实验中,水质点径向速度空时分布序列已知,可用于实现波浪重建方法验证。假设雷达位于坐标原点,雷达探测覆盖范围设置为0~1 500 m。如图6和图7所示为不同风速下t=0时刻(当前时刻)重建波高与参考波高的比较结果。

由图6和图7中可以看出,不同海况下重建的波浪与参考波浪吻合效果好,并且三种方法中正则化最小二乘法的重建效果显然最好。如表1所示为风速为8 m/s和15 m/s时重建波高与参考波高的相关系数以及均方根误差。相较于势函数法,正则化最小二乘法在不同风速下参考波高与重建波高的相关系数(correlation coefficient,COR)更接近1,且均方根误差(root mean squared error,RMSE)也更小,该结果进一步验证了正则化最小二乘法可提高波浪重建结果的准确性。

3.2 预测

对于二维波浪场,可预测区可通过对波浪进行时空间采样进行划定。根据1.3节中对于可预测区模型的分析可知,最小、最大波数k_g_l(r,t)、k_g_h(r,t)分别对应理论上可预测空时位置的最小、最大群速度C_g_l(r,t)、C_g_h(r,t),从雷达观测区域传播到任意空时位置的最大、最小群速度分别为

(32)C_g_l(r,t)=$\frac{\mathit{r}-{\mathit{R}}_{0}}{\mathit{t}-{\mathit{t}}_{0}}$

(33)C_g_h(r,t)=$\frac{\mathit{r}-({\mathit{R}}_{0}+\mathit{R})}{\mathit{t}-{\mathit{t}}_{0}}$

在仿真实验中,采用蒙特卡洛模拟的方法分析波高预测性能,并利用归一化失配误差(normalized mean error,NME)来评估预测波高与参考波高之间的误差,表达式为

(34)NME=$\frac{\stackrel{\mathit{M}}{\sum _{\mathit{s}=1}}\left|{\mathit{\eta }}_{\mathrm{P},\mathit{s}}(\mathit{r},\mathit{\varphi },\mathit{t})-\right.\left.{\mathit{\eta }}_{\mathrm{R},\mathit{s}}(\mathit{r},\mathit{\varphi },\mathit{t})\right|/{\mathit{H}}_{\mathit{s}}}{{\mathit{M}}_{\mathit{s}}}$

式(34)中:η_P_s(r,ϕ,t)为预测瞬时波高;η_R_s(r,ϕ,t)为参考瞬时波高;M_s为蒙特卡罗模拟实验次数;H_s为有效波高。

仿真实验中,设置雷达的探测距离约为750 m,距离分辨率为7.5 m,时间分辨率为0.5 s。如图8所示为不同风速下二维波浪瞬时波高预测误差的空时分布,其中白线划分的区域为理论可预测区域。结果显示使用正则化最小二乘方法计算的波高分布,分别在距离向和时间向都可以实现一定范围内的波高预测,在空时谱中可预测区域呈现菱形,其大小由群速度决定。在可预测区内,参考波高与预测波高非常接近,但随着空时位置远离可预测区域,预测误差逐渐增大,直至参考波高与预测波高完全不相关。

如图9和图10所示为海面上方10 m处风速分别为8 m/s和15 m/s时,位置r = 0 m处二维波浪场预测结果和预测误差的仿真结果。可以看到,在可预测区内,瞬时波高的预测误差均为0 m,因此使用本文所提方法可在一定空时范围内实现瞬时波高的准确预测,并且随着空时位置远离可预测区,误差逐渐增大,直至为1。不同风速对应不同的海况,图9中海面风速为8 m/s,对应低海况,图10中海面风速为15 m/s,即为高海况。图9中低海况下,可预测时长约为155 s,而图10中高海况下,可预测时长约为125 s。随着海况升高,海面复杂程度增大,可预测时长逐渐降低。

4 结论

基于线性波浪理论,分析了海面瞬时波高的几何模型,并使用正则化最小二乘法实现波浪重建及预测。正则化最小二乘法不仅克服了势函数法无法进行波浪预测的缺点,也解决了最小二乘法对异常值较为敏感的问题。对于波浪重建,分别介绍了基于势函数法和正则化最小二乘法的波浪重建方法,具体包括算法原理和仿真分析,以此来验证算法的准确性。结果显示在不同风速下,与势函数法的结果相比,正则化最小二乘法重建的瞬时波高与参考波高的相关性更高,且均方根误差更小。针对波浪预测,结合线性海面模型,求解模型参数,即可实现对可预测区内任意时空位置的波浪预测。因此,构建了二维波浪场的可预测区边界表达式,在可预测区域内使用正则化最小二乘法实现不同空时位置的海面瞬时波高预测。不同海况下的仿真结果表明,在可预测区内,参考波高与预测波高非常接近,但随着空时位置远离可预测区域,预测误差逐渐增大,直至参考波高与预测波高完全不相关。

参考文献

原文顺序 | 出版日期 | 本文引用

[1]	郑耿忠, 宋亚男, 徐荣华, 等. 基于图像序列的波浪三维重建算法优化[J]. 无线互联科技, 2022, 19(18): 135-138.

[2]	Zheng Gengzhong, Song Yanan, Xu Ronghua, et al. Optimization of wave 3D reconstruction algorithms based on image sequence[J]. Wireless Internet Technology, 2022, 19(18): 135-138.

[3]	Chen Z, Chen X, Zhao C, et al. Observation and intercomparison of wave motion and wave measurement using shore-based coherent microwave radar and HF radar[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 57(10): 7594-7605.

[4]	Chen Z, Chen X, Zhao C, et al. Wave height and wave period derived from a shipboard coherent S-band wave radar in the South China Sea[J]. Remote Sensing, 2019, 11(23): 1-22.

[5]	Wright J. Backscattering from capillary waves with application to sea clutter[J]. IEEE Transactions on Geoscience and Remote Sensing, 1966, 14(6): 749-754.

[6]	Keller W, Wright J. Microwave scattering and the straining of wind-generated waves[J]. Radio Science, 2012, 10(2): 139-147.

[7]	Hessner K, Nieto-Borge J, Bell P. Nautical radar measurements in Europe: applications of WaMos II as a sensor for sea state, current and bathymetry[J]. Current and Bathymetry, 2008, 90: 435-446.

[8]	Dankert H, Horstmann J, Lehner S, et al. Detection of wave groups in sar images and radar image sequences[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(6): 1437-1446.

[9]	Dankert H, Rosenthal W. Ocean surface determination from X-band radar-image sequences[J]. Journal of Geophysical Research: Oceans, 2004, 109(C4): 1-11.

[10]	Dankert H, Horstmann J, Rosenthal W. Wind and wave-field measurements using marine X-band radar-image sequences[J]. IEEE Journal of Oceanic Engineering, 2005, 30(3): 534-542.

[11]	Reichert K, Dannenberg J, Boom H. X-band radar derived sea surface elevation maps as input to ship motion forecasting[C]// OCEANS Sydney 2010 IEEE. Sydney: IEEE, 2010: 5603968.

[12]	Plant W, Keller W, Wright J. Modulation of coherent microwave backscatter by shoaling waves[J]. Journal of Geophysical Research: Oceans, 2012, 83(C3): 1347-1352.

[13]	魏鋆宇, 陈泽宗, 赵晨, 等. S波段船载测波雷达接收机设计及实现[J]. 科学技术与工程, 2021, 21(17): 7174-7179.

[14]	Wei Yunyu, Chen Zezong, Zhao Chen, et al. Design and implementation of S-band shipboard wave radar receiver[J]. Science Technology and Engineering, 2021, 21(17): 7174-7179.

[15]	Johnson J, Burkholder R, Toporkov J, et al. A numerical study of the retrieval of sea surface height profiles from low grazing angle radar data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(6): 1641-1650.

[16]	Nwogu O, Lyzenga D. Surface-wavefield estimation from coherent marine radars[J]. IEEE Geoscience and Remote Sensing Letters, 2010, 7(4): 631-635.

[17]	Connell B, Rudzinsky J, Brundick C, et al. Development of an environmental and ship motion forecasting system[C]// Proceedings of the ASME 2015 34th International Conference on Ocean(OMAE2015). St. John’s: Offshore and Arctic Engineering, Canada, 2015: 1-11.

[18]	Støle-Hentschel S, Seemann J, Nieto-Borge J, et al. Consistency between sea surface reconstructions from nautical X-band radar Doppler and amplitude measurements[J]. Journal of Atmospheric and Oceanic Technology, 2018, 35(6): 1201-1220.

[19]	Morris E, Zienkiewicz H, Belmont M. Short term forecasting of the sea surface shape[J]. International Shipbuilding Progress, 1998, 45(444): 383-400.

[20]	Wu G. Direct simulation and deterministic prediction of large-scale nonlinear ocean wave-field[D]. Cambridge: Massachusetts Institute of Technology, 2004.

[21]	Naaijen P, Trulsen K, Blondel-Couprie E. Limits to the extent of the spatio-temporal domain for deterministic wave prediction[J]. 2014, 61(3/4): 203-223.

[22]	Kosleck S. Prediction of wave-structure interaction by advanced wave field forecast[C]// Annual International Conference of the German Society for Maritime Technology(STG). Berlin: Society for Maritime Technology, 2013: 1-20.

[23]	马学文. 基于深度学习的相位解析波浪传播与船舶运动安稳期预报研究[D]. 哈尔滨: 哈尔滨工程大学, 2022.

[24]	Ma Xuewen. Phase-resolved wave propagation and ship motion quiet period prediction based on deep learning[D]. Harbin: Harbin Engineering University, 2022.

[25]	Kim I, Ducrozet G, Leroy V, et al. A real-time wave prediction in directional wave fields: strategies for accurate continuous prediction in time[J]. Ocean Engineering, 2024, 291: 116445.

[26]	李汉东, 陈泽宗, 赵晨, 等. 基于船舶运动的微波多普勒海面回波模拟与分析[J]. 科学技术与工程, 2021, 21(8): 3147-3152.

[27]	Li Handong, Chen Zezong, Zhao Chen, et al. Simulation and analysis of microwave doppler sea surface echo based on ship motion[J]. Science Technology and Engineering, 2021, 21(8): 3147-3152.

[28]	Ogunbiyi O, Massenat N, Quinete N. Dispersion and stratification of Per-and polyfluoroalkyl substances (PFAS) in surface and deep-water profiles: a case study of the Biscayne Bay area[J]. Science of the Total Environment, 2024, 909: 168413.

[29]	Schroeder K, Jacopo C. Oceanography of the mediterranean sea: an introductory guide[M]. New York: Elsevier, 2022.

[30]	国强, 李文韬. 基于正则化约束总体最小二乘的TDOA/FDOA无源定位方法[J]. 哈尔滨工业大学学报, 2022, 54(5): 81-87.

[31]	Guo Qiang, Li Wentao. Passive TDOA/FDOA location based on regularized constrained total least squares[J]. Journal of Harbin Institute of Technology, 2022, 54(5): 81-87.

[32]	于雨. 基于最大熵法的波浪浮标方向谱反演算法[J]. 气象水文海洋仪器, 2022, 39(2): 1-3.

[33]	Yu Yu. Wave buoy direction spectral inversion algorithm based on maximum entropy method[J]. Meteorology, Hydrology, and Marine Instruments, 2022, 39(2): 1-3.

[34]	Ben-Mhenni R, Bourguignon S, Ninin J. Global optimization for sparse solution of least squares problems[J]. Optimization Methods and Software, 2021, 37(5): 1740-1769.

[35]	Matysik O, Van Hulle M. Alternating step size method for solving ill-posed linear operator equations in energetic space[J]. Journal of Computational and Applied Mathematics, 2022, 416: 1-10.

[36]	Hansen P, O’Leary D. The use of the L-curve in the regularization of discrete ill-posed problems[J]. SIAM Journal on Scientific Computing, 1993, 14(6): 1487-1503.