CN114814774B - 一种高频雷达探测船尾波的信号建模方法 - Google Patents

Abstract

本发明公开了一种高频雷达探测船尾波的信号建模方法，应用Wigley船体模型，通过仿真舰船Kelvin尾波波高图，得到其波高谱。Kelvin尾波波高谱与海浪谱相加得到叠加谱，代入高频雷达海面后向散射模型进行海洋回波多普勒谱的仿真。对结果加窗平滑，并添加噪声后，最终仿真得到了高频雷达船尾波信号。首次完成了该信号的仿真工作，验证了在海况相对平静，船体尺寸较大，船速较快的前提下，高频雷达在海杂波叠加下探测船尾波是可行的。

Description

一种高频雷达探测船尾波的信号建模方法

技术领域

本发明属于海洋遥感领域，尤其涉及一种高频雷达探测船尾波的信号建模方法。

背景技术

船只在海洋中航行会产生船尾波，探测船尾波对于船只航行管理和海洋环境监测都具有重要的意义，最近受到了广泛关注。

船尾波的特殊性在于，无论舰船使用何种动力，只要在行进中，船体都会与水体发生相互作用，在船尾产生扇形区域的波动。同时，通过大量的观测和实验测量，人们发现船尾波具有范围大、持续时间长的特点。船尾波探测不仅可以弥补雷达传统探测方法的局限不足，而且可能获得有助于目标分类的信息以及给未来舰船的设计提供参考。

海洋雷达探测船尾波的研究一直不曾间断，不过绝大部分是基于微波合成孔径雷达，关于高频雷达船尾波探测的工作很少。综上所述，研究高频雷达船尾波探测是有意义的，而高频雷达船尾波信号的建模方法也缺乏系统的研究。

发明内容

本发明针对现有技术的不足，提供一种高频雷达探测船尾波的信号建模方法。该方法的提出填补了船尾波探测研究的一个空白。

本发明的技术方案为一种高频雷达探测船尾波的信号建模方法，包含以下步骤：

步骤1，Kelvin尾波波高图仿真。以船体吃水线所在平面中心点为原点，原点指向船尾为x轴正方向，垂直于吃水线平面向上为z轴正方向建立右手坐标系。兼顾计算复杂度和贴近实际的要求，选定仿真所用船体几何模型为Wigley椭圆船型。根据船尾波的特征，选定船尾一个边长为L的正方形区域作为仿真区域。仿真区域四个角分别为(O,±L/2)和(L,±L/2)。

步骤1.1，给定仿真参数。即船体几何参数和船速U。由船速可算得纯横断波波数k₀＝g/U²，g是重力加速度。

步骤1.2，在参数选定后，进行离散化操作。首先，给定离散化点数，对仿真区域x轴数值区间和y轴数值区间进行离散化得到两个行向量和/>其次，给定离散化点数，对船体中心平面长度方向x轴数值区间和吃水深度方向z轴数值区间进行离散化得到两个行向量和/>最后，给定离散化点数，离散波传播角数值区间得到行向量/>

步骤1.3，仿真计算。根据Kelvin尾波波高公式，先通过对船体中心平面的双重积分得到舰船自由波谱，再对波传播角进行积分得到仿真区域离散点的波高。最终计算得到二维矩阵ζ，即Kelvin尾迹波高图。

步骤2，Kelvin尾波波高谱仿真。对步骤1所得二维矩阵ζ做二维离散傅里叶变换得到二维矩阵ζ_k。ζ_k取模的平方再除以L²得到二维矩阵S_wake，即Kelvin尾迹波高谱。

步骤3，高频雷达海面回波谱仿真。选择散射模型为Barrick根据微扰法建立的单基地高频雷达海面后向散射模型，仿真输入的有向海浪谱由PM谱和Longuet-Higgins方向因子组成。

步骤3.1，参数赋值。包括海浪谱参数和雷达参数。

步骤3.2，参数归一化。给出归一化波矢量、归一化多普勒频率、归一化耦合系数、归一化海浪谱和归一化回波谱的定义。

步骤3.3，离散化。给定离散化点数，离散归一化频率区间得到行向量将海浪谱离散化得到S_wave矩阵，并与步骤2得到的S_wake矩阵相加，得到叠加谱矩阵S。

步骤3.4，归一化一阶谱计算。根据狄利克雷函数限制条件得出一阶谱值不为零的归一化频率点。通过查表S计算得到对应点的一阶谱数值。最终得到的一阶谱向量记为σ₁。

步骤3.5，归一化二阶谱计算。二阶谱狄利克雷函数限制较为复杂，需要进一步拆解计算。

步骤3.5.1，降重积分。通过定义新变量，将二阶谱计算的双重积分变为单重积分。

步骤3.5.2，对于每一个归一化频率值，得出积分限，离散积分变量，通过解方程和查表得到对应点二阶谱数值。最终得到二阶谱向量记为σ₂。

步骤3.6，加窗加噪。归一化一阶谱与二阶谱相加得到归一化回波谱，进一步得到高频雷达海面回波多普勒谱σ_f。而为了仿真更接近实测。要对σ_f加窗进行平滑以拓宽一阶峰，以及按照指定信噪比SNR添加噪声，最终得到仿真结果σ_E。

步骤1所述船体几何模型如下：

其中l是半船长，b是半船宽，D是吃水深度，(x,z)是船体中心平面坐标点，Y(x,z)是对应坐标点的船体偏移量；

步骤1.2，离散化得到两个行向量和/>如下：

其中，仿真区域x轴数值区间是[0,L]，y轴数值区间是[-L/2,L/2]，M和N是离散点数；

离散化得到的行向量和/>如下：

其中，船体中心平面长度方向x轴数值区间是[-l,l]，船体中心平面吃水深度方向z轴数值区间是[-D,0]，N_x和N_z是离散点数；

离散化得到的行向量如下：

其中，波传播角θ的数值范围是[-π/2,π/2]，N_θ是离散点数；

步骤1.3，对每一个点(x_i,y_i)(i＝0,1,...,M,j＝0,1,...,N)，波高为：

其中，对每一个θ_k(k＝0,1,...,Nθ)，有

A(θ)也称自由波谱，是船的特有属性。其中

而对每一个x′_p(p＝0,1,...,N_x)有

步骤3，根据单基地高频雷达海面后向散射模型，一阶谱和二阶谱如下：

其中，ω是多普勒角频率，k₀是雷达波数，是雷达波矢量，ω_B是布拉格角频率，Γ^L是耦合系数，m,m′是多普勒频移标志，/>和/>是海浪波矢量，/>和/>是有向海浪谱；p和q分别是平行和垂直于雷达波束方向的波矢量平面坐标。

当雷达工作频率为f₀，角频率ω₀＝2πf₀，波数c是光速；布拉格频率布拉格角频率ω_B＝2πf_B；耦合系数Γ^L＝Γ^H+Γ^EM，Γ^H和Γ^EM分别是流体力学耦合系数和电磁耦合系数；

有向海浪谱公式如下：

其中，f(k)是无向海浪谱，g(θ)是扩散函数；

仿真中使用的无向海浪谱模型PM谱公式如下：

其中，k是海浪波波数，是截止波数，g是重力加速度，U_19.5是距离海面19.5m高度处的风速；

仿真中使用的扩散函数Longuet-Higgins方向因子公式如下：

其中，θ是海浪波与雷达波束的夹角，s是扩展因子，是主浪向，；

步骤3.2，归一化波矢量定义为：

归一化多普勒频率：

η＝ω/ω_B

归一化耦合系数：

γ^L＝γ^H+γ^EM

归一化海浪谱：

F(K)＝(2k₀)⁴f(k)

归一化方向谱：

归一化一阶谱：

归一化二阶谱：

步骤3.3，离散化得到的行向量如下：

其中，归一化频率η的数值范围是[-η_lim,η_lim]，离散点数为N_η；

步骤3.4，当η_a＝1,对矩阵S查表得到又有

步骤3.5.1，定义如下变量：

则二阶谱可重新表示为：

其中，由/>给出；

步骤3.5.2，对每一个η_a(a＝0,1,...,N_η)，有积分限

离散化得到的行向量如下：

其中，角θ的数值范围是[-θ_L,θ_L]，离散点数为N_θ；

对每一个θ_b(b＝0,1,...,N_θ)，解如下方程：

可求出进而求出：

θ′_b＝sin^-1(K_b sinθ_b/K′_b)+π

又有海浪谱公式：

其中，和/>可查表S可得到；

归一化耦合系数计算公式如下：

至此，计算

步骤3.6，由归一化操作，有

σ(η)＝ω_Bσ(f)

即

一种高频雷达探测船尾波的信号建模方法，其意义在于为高频雷达探测船尾波的后续研究提供支持。通过仿真证实了，在海况温和、船体尺寸和船速比较大的前提下，叠加于海杂波的船尾波强度属于高频雷达可探测的级别。后续利用此方法，还可以验证高频雷达船尾波信号目标分类的能力。此外，此方法还可延拓至潜艇船尾波探测的研究，有很大的探索空间。

附图说明

图1是本发明提供的技术方案流程图。

图2是本发明技术方案Kelvin尾波波高图仿真流程图。

图3是本发明技术方案高频雷达船尾波信号仿真流程图。

图4是本发明技术方案二阶谱仿真流程图。

图5是本发明实施例船体模型图。

图6是本发明实施例Kelvin尾波波高图仿真结果。

图7是本发明实施例Kelvin尾波波高谱仿真结果。

图8是本发明实施例高频雷达船尾波信号仿真结果。

具体实施方式

本发明提供的方法能够用计算机软件技术实现流程。下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

步骤1，Kelvin尾波波高图仿真。以船体吃水线所在平面中心点为原点，原点指向船尾为x轴正方向，垂直于吃水线平面向上为z轴正方向建立右手坐标系。仿真所用船体几何模型为鉴于船尾波的特征，选定船尾一个正方形区域作为仿真区域，边长为L＝800m，四个角分别为(O,±400)和(800,±400)。

步骤1.1，给定仿真参数。即半船长l＝31、半船宽b＝4.8m、吃水深度D＝2m和船速U＝15knot。纯横断波波数k₀＝g/U²＝9.8/(15*0.514)²＝0.165m^-1，g＝9.8m/s²是重力加速度。附图5是此参数下船体模型图。

步骤1.2，在参数选定后，进行离散化操作。仿真区域x轴数值区间是[0,800]，y轴数值区间是[-400,400]。离散化点数分别为M＝800和N＝800，分别离散化得到行向量(x₀,x₁,...,x_i,...,x_M-1,x_M)和(y₀,y₁,...,y_j,...,y_N-1,y_N)。船体中心平面长度方向x轴数值区间是[-31,31]，离散点数为N_x＝1860，离散化得行向量船体中心平面吃水深度方向z轴数值区间是[-2,0]，离散点数为N_z＝60，离散化得到行向量波传播角θ的数值范围是/>离散点数为N_θ＝500，离散化得到行向量/>

步骤1.3，仿真计算。对每一个点(x_i,y_i)(i＝0,1,...,800,j＝0,1,...,800)，有波高其中，对每一个θ_k(k＝0,1,...,500)，有/>A(θ)也称自由波谱，是船的特有属性。其中，/> 而对每一个x′_p(p＝0,1,...,1860)有/>最终计算得到大小为800×800的二维矩阵ζ，即Kelvin尾迹波高图。附图6即所绘制Kelvin尾迹波高图。

步骤2，Kelvin尾波波高谱仿真。对步骤1所得二维矩阵ζ做二维离散傅里叶变换得到二维矩阵ζ_k。ζ_k取模的平方再除以L²得到二维矩阵S_wake，即Kelvin尾迹波高谱。附图7即所绘制Kelvin尾迹波高谱。

步骤3，高频雷达海面回波谱仿真。根据单基地高频雷达海面后向散射模型，一阶谱表示为：二阶谱表示为：。仿真使用的有向海浪谱由PM谱/>和Longuet-Higgins方向因子组成。有向海浪谱/> 是波矢量。耦合系数Γ^L＝Γ^H+Γ^EM。m,m′是多普勒频移标志。

步骤3.1，参数赋值。g＝9.8m/s²是重力加速度，U_19.5＝8m/s是距离海面19.5m高度处的风速；s＝2是扩展因子，/>是主浪向。雷达工作频率为f₀＝15MHz，角频率ω₀＝2πf₀，波数/>c＝3×10⁸m/s是光速，/>是雷达波矢量。布拉格频率/>布拉格角频率ω_B＝2πf_B。

步骤3.2，参数归一化。波矢量归一化波矢量：

归一化多普勒频率：η＝ω/ω_B；

归一化耦合系数

归一化海浪谱：F(K)＝(2k₀)⁴f(k),K_c＝k_∈/2k₀；

归一化方向谱：

归一化一阶谱：

归一化二阶谱：

步骤3.3，离散化。取归一化多普勒频率η范围为[-2,2]，离散化点数为N_η＝801，离散化得到行向量将海浪谱离散化得到S_wave矩阵，并与步骤2得到的S_wake矩阵相加，得到叠加谱矩阵S。

步骤3.4，归一化一阶谱计算。

归一化一阶谱

对每一个η_a(a＝0,1,...,N_η)，计算σ₁(η_a)。

当η_a≠±1时，σ₁(η_a)＝0；

当η_a＝1,对矩阵S查表得到，又有/>

进而

如此得到的一阶谱向量记为σ₁。

步骤3.5，归一化二阶谱计算。归一化二阶谱

步骤3.5.1，降重积分。定义 二阶谱/>由/>给出。

步骤3.5.2，对每一个η_a(a＝0,1,...,N_η)，计算σ₂(η_a)。当η_a ²≤2，θ_L＝π；当η_a ²＞2，θ_L＝π-cos^-1(2/η_a ²)。给出离散点数为N_θ＝50000，离散化得到行向量则/>对每一个θ_b(b＝0,1,...,N_θ)，解方程求出/>进而求出/>θ′_b＝sin^-1(K_bsinθ_b/K′_b)+π。又海浪谱/>和可查表S得到。耦合系数

计算和

求和操作即可求得σ₂(η_a)。N_η个值组成的二阶谱向量记为σ₂。

步骤3.6，加窗加噪。根据步骤2归一化操作，有σ(η)＝f_Bσ(f)，

即

σ_f就是仿真所得高频雷达海面回波多普勒谱。而为了仿真更接近实测。要对σ_f加窗进行平滑以拓宽一阶峰，以及按照指定信噪比SNR添加噪声，最终得到仿真结果σ_E。附图8即为仿真所得高频雷达船尾波信号。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (1)

1.一种高频雷达探测船尾波的信号建模方法，其特征在于，包括以下步骤：

步骤1，Kelvin尾波波高图仿真；以船体吃水线所在平面中心点为原点，原点指向船尾为x轴正方向，垂直于吃水线平面向上为z轴正方向建立右手坐标系；兼顾计算复杂度和贴近实际的要求，选定仿真所用船体几何模型为Wigley椭圆船型；根据船尾波的特征，选定船尾一个边长为L的正方形区域作为仿真区域；仿真区域四个角分别为(O,±L/2)和(L,±L/2)；

步骤1.1，给定仿真参数；即船体几何参数和船速U；由船速可算得纯横断波波数k₀＝g/U²，g是重力加速度；

步骤1.2，在参数选定后，进行离散化操作；首先，给定离散化点数，对仿真区域x轴数值区间和y轴数值区间进行离散化得到两个行向量和/>其次，给定离散化点数，对船体中心平面长度方向x轴数值区间和吃水深度方向z轴数值区间进行离散化得到两个行向量/>和最后，给定离散化点数，离散波传播角数值区间得到行向量/>

步骤1.3，仿真计算；根据Kelvin尾波波高公式，先通过对船体中心平面的双重积分得到舰船自由波谱，再对波传播角进行积分得到仿真区域离散点的波高；最终计算得到二维矩阵ζ，即Kelvin尾迹波高图；

步骤2，Kelvin尾波波高谱仿真；对步骤1所得二维矩阵ζ做二维离散傅里叶变换得到二维矩阵ζ_k；ζ_k取模的平方再除以L²得到二维矩阵S_wake，即Kelvin尾迹波高谱；

步骤3，高频雷达海面回波谱仿真；选择散射模型为Barrick根据微扰法建立的单基地高频雷达海面后向散射模型，仿真输入的有向海浪谱由PM谱和Longuet-Higgins方向因子组成；

步骤3.1，参数赋值；包括海浪谱参数和雷达参数；

步骤3.2，参数归一化；给出归一化波矢量、归一化多普勒频率、归一化耦合系数、归一化海浪谱和归一化回波谱的定义；

步骤3.3，离散化；给定离散化点数，离散归一化频率区间得到行向量将海浪谱离散化得到S_wave矩阵，并与步骤2得到的S_wake矩阵相加，得到叠加谱矩阵S；

步骤3.4，归一化一阶谱计算；根据狄利克雷函数限制条件得出一阶谱值不为零的归一化频率点；通过查表S计算得到对应点的一阶谱数值；最终得到的一阶谱向量记为σ₁；

步骤3.5，归一化二阶谱计算；二阶谱狄利克雷函数限制较为复杂，需要进一步拆解计算；

步骤3.5.1，降重积分；通过定义新变量，将二阶谱计算的双重积分变为单重积分；

步骤3.5.2，对于每一个归一化频率值，得出积分限，离散积分变量，通过解方程和查表得到对应点二阶谱数值；最终得到二阶谱向量记为σ₂；

步骤3.6，加窗加噪；归一化一阶谱与二阶谱相加得到归一化回波谱，进一步得到高频雷达海面回波多普勒谱σ_f；而为了仿真更接近实测；要对σ_f加窗进行平滑以拓宽一阶峰，以及按照指定信噪比SNR添加噪声，最终得到仿真结果σ_E；

步骤1所述船体几何模型如下：

其中l是半船长，b是半船宽，D是吃水深度，(x,z)是船体中心平面坐标点，Y(x,z)是对应坐标点的船体偏移量；

步骤1.2，离散化得到两个行向量和/>如下：

其中，仿真区域x轴数值区间是[0,L]，y轴数值区间是[-L/2,L/2]，M和N是离散点数；

离散化得到的行向量和/>如下：

其中，船体中心平面长度方向x轴数值区间是[-l,l]，船体中心平面吃水深度方向z轴数值区间是[-D,0]，N_x和N_z是离散点数；

离散化得到的行向量如下：

其中，波传播角θ的数值范围是[-π/2,π/2]，N_θ是离散点数；

步骤1.3，对每一个点(x_i,y_i)(i＝0,1,...,M,j＝0,1,...,N)，波高为：

其中，对每一个θ_k(k＝0,1,...,N_θ)，有

A(θ)也称自由波谱，是船的特有属性；其中

而对每一个x′_p(p＝0,1,...,N_x)有

步骤3，根据单基地高频雷达海面后向散射模型，一阶谱和二阶谱如下：

其中，ω是多普勒角频率，k₀是雷达波数，是雷达波矢量，ω_B是布拉格角频率，Γ^L是耦合系数，m,m′是多普勒频移标志，/>和/>是海浪波矢量，/>和/>是有向海浪谱；p和q分别是平行和垂直于雷达波束方向的波矢量平面坐标；

当雷达工作频率为f₀，角频率ω₀＝2πf₀，波数c是光速；布拉格频率/>布拉格角频率ω_B＝2πf_B；耦合系数Γ^L＝Γ^H+Γ^EM，Γ^H和Γ^EM分别是流体力学耦合系数和电磁耦合系数；

有向海浪谱公式如下：

其中，f(k)是无向海浪谱，g(θ)是扩散函数；

仿真中使用的无向海浪谱模型PM谱公式如下：

其中，k是海浪波波数，是截止波数，g是重力加速度，U_19.5是距离海面19.5m高度处的风速；

仿真中使用的扩散函数Longuet-Higgins方向因子公式如下：

其中，θ是海浪波与雷达波束的夹角，s是扩展因子，是主浪向；

步骤3.2，归一化波矢量定义为：

归一化多普勒频率：

η＝ω/ω_B

归一化耦合系数：

γ^L＝γ^H+γ^EM

归一化海浪谱：

F(K)＝(2k₀)⁴f(k)

归一化方向谱：

归一化一阶谱：

归一化二阶谱：

步骤3.3，离散化得到的行向量如下：

其中，归一化频率η的数值范围是[-η_lim,η_lim]，离散点数为N_η；

步骤3.4，当η_a＝1,对矩阵S查表得到又有

步骤3.5.1，定义如下变量：

则二阶谱可重新表示为：

其中，由/>给出；

步骤3.5.2，对每一个η_a(a＝0,1,...,N_η)，有积分限

离散化得到的行向量如下：

其中，角θ的数值范围是[-θ_L,θ_L]，离散点数为N_θ；

对每一个θ_b(b＝0,1,...,N_θ)，解如下方程：

可求出进而求出：

θ′_b＝sin^-1(K_bsinθ_b/K′_b)+π

又有海浪谱公式：

其中，和/>可查表S可得到；

归一化耦合系数计算公式如下：

至此，计算

步骤3.6，由归一化操作，有

σ(η)＝ω_Bσ(f)

即