CN115236664A - 一种航海雷达图像反演有效波高的方法 - Google Patents

Abstract

本发明公开了一种航海雷达图像反演有效波高的方法，首先利用离线的浮标数据和雷达，通过最小二乘法拟合得到海浪的波峰周期和平均过零周期之间的关系，再通过3D‑FFT谱分析法得到二维海浪谱，从二维海浪谱中获取海浪的波峰周期，然后计算平均过零周期，通过阴影法计算海浪的均方根波陡，最后通过平均过零周期和均方根波陡参数计算海浪的有效波高。本发明能独立地通过雷达数据反演海浪的有效波高参数，使得利用雷达测量海浪有效波高参数的过程完全脱离了浮标的支持，不需要借助于其他仪器，增加了算法的适用性和独立性。

Description

一种航海雷达图像反演有效波高的方法

技术领域

本发明属于利用遥感手段进行海浪参数反演的领域，涉及一种航海雷达图像反演有效波高的方法，特别是一种基于组合策略的航海雷达图像反演有效波高的方法。

背景技术

近年来利用X波段航海雷达对海浪进行实时遥测是海洋监测领域的一个研究热点。目前，利用X波段航海雷达图像反演有效波高的两种主要方法分别是谱分析法和阴影法。

谱分析法是综合利用了海杂波图像的时间和空间信息的算法，该方法通过对连续的海杂波图像进行FFT变换、滤波、积分等操作，最终得到海浪的频率方向谱和信噪比，从海浪的频率方向谱中可进一步得到海浪的波向和频率信息，通过信噪比和有效波高之间的线性关系可计算得到有效波高信息。

谱分析法的优点在于可以同时得到海浪的多种参数——波向、频率、有效波高，遮挡对谱分析法的影响可以忽略不计，这是许多方法不可比拟的优点。但是谱分析法中信噪比和有效波高之间的系数需要通过浮标实验来确定，这增加了谱分析法的工程应用成本，另外谱方法需要连续的多副图像，算法计算时间较长，这也限制了该方法的应用。

2014年，挪威海洋物理学家Rune[13]综合前人的研究，明确指出了海杂波主要是由阴影调制产生，并且提出了不需要标定参数的阴影法来计算有效波高。阴影法相比于传统的谱分析法，其最大的优点在于阴影法的计算完全建立在理论推导的基础上，不需要浮标实验来标定相关的系数。另外，阴影法中只需要一幅雷达图像即可进行计算，算法时间成本低，增强了算法的适用性。但是在目前的研究过程中阴影法中海浪的平均过零周期均是由浮标辅助提供的，这一条件限制了阴影法的应用和发展。

发明内容

针对上述现有技术，本发明要解决的技术问题是提供一种基于3D-FFT谱分析法和阴影法组合策略的航海雷达图像反演有效波高的方法，解决3D-FFT谱分析法算法计算时间长、工程应用成本高和阴影统计法中波高参数计算依赖浮标数据的问题。

为解决上述技术问题，本发明的一种航海雷达图像反演有效波高的方法，包括以下步骤：

步骤1：利用离线的浮标数据和雷达，通过最小二乘法拟合得到海浪的波峰周期和平均过零周期之间的关系式；

步骤2：选定研究区域，通过3D-FFT谱分析法得到选定研究区域的二维海浪谱，从二维海浪谱中获取海浪的波峰周期；

步骤3：将步骤2得到的海浪的波峰周期代入步骤1得到的海浪波峰周期和平均过零周期之间的关系式，计算得到平均过零周期T_m02；

步骤4：在原始雷达图像中选取研究区域，在选定的研究区域中将原始雷达图像径向维等分为m个子区域；计算区分阴影区与非阴影区的阈值τ_S；根据τ_S计算各子区域的阴影比例函数；然后根据阴影比例函数计算得到均方根波陡σ_RMS；

步骤5：通过平均过零周期T_m02和均方根波陡σ_RMS计算海浪的有效波高H_S，H_S满足：

式中，g为重力加速度。

进一步的，海浪的波峰周期和平均过零周期之间的关系具体为：

T_m02＝aT_M+b

式中，T_M为波峰周期，T_m02为平均过零周期；

记离线实验中浮标的平均过零周期序列为T_m02i，i＝1,2,3...n，对应的通过3D-FFT谱分析法利用雷达数据得到的波峰周期序列记为T_Mi，i＝1,2,3...n，根据最小二乘法得到：

式中，

和

分别是平均过零周期序列和波峰周期序列的平均值。

进一步的，步骤2包括：

步骤2.1：通过最近点插值，将原始雷达图像中极坐标下的回波强度值转换到笛卡尔坐标系下对应的点上，记原始雷达图像中某一点的坐标为(r₀,θ₀)，其转换后对应的笛卡尔坐标为(x₀,y₀)，插值后得到的雷达图像序列记为η(x,y,t)，转换公式为：

式中，r₀为原始雷达图像中某一点的径向距离，θ₀为原始雷达图像中某一点的方位角，x₀,y₀分别为笛卡尔坐标系的横坐标值和纵坐标值，sqrt( )、round( )和rem( )分别是开平方函数、取整函数、和求余函数；

步骤2.2：选取的研究区域大小为L_x*L_y的笛卡尔框，L_x,L_y为笛卡尔框沿x和y方向的范围，雷达图像序列的总时间为T，对雷达图像序列η(x,y,t)进行三维傅里叶变换得到三维图像谱，然后利用变换公式的离散形式将图像序列从时间空间域转化到波数频率域：

式中，k_x和k_y分别为x和y方向的波数，ω为角频率，则分辨率的计算公式为：

步骤2.3：假设波数处对应的理论频率为ω₀，真实频率ω通过色散关系表示，则两者之间插值的平方和为：

将三维图像谱中的能量值作为权重，得到目标函数为：

当SSE关于的偏导数均为零时，SSE有最小解，解算出流速为：

式中，h为水深，g为重力加速度，F为三维图像谱F(k_x,k_y,ω)，海表层流

波数

其中，k_x,k_y分别为x和y方向的波数；

步骤2.4：利用基于色散关系的带通滤波器对三维图像谱进行滤波，采用的滤波器是基于最大流速的带通滤波器：

式中，E(k_x,k_y,ω)为滤波后得到的三维图像谱，其中，

式中，k为波数，dω为频率分辨率，dk＝(dk_x＝dk_y)为波数分辨率，U_max为最大流速，当B_n取值为负时，则令B_n＝0；

再对滤波后得到的三维图像谱E(k_x,k_y,ω)进行频率ω维积分得到二维图像谱I(k_x,k_y)：

I(k_x,k_y)＝∫_ω＞0E(k_x,k_y,ω)dω

步骤2.5：通过经验调制传递函数将二维图像谱I(k_x,k_y)转换为二维海浪谱E(k_x,k_y)：

E(k_x,k_y)＝|M(k_x,k_y)|²·I(k_x,k_y)

经验调制传递函数为：

|M(k_x,k_y)|²≈k^-β

其中，β为经验系数；

步骤2.6：通过二维海浪谱E(k_x,k_y)利用色散关系计算海浪的波峰周期T_M，计算公式为：

式中，g为重力加速度，k_m是海浪谱E(k_x,k_y)中能量值最大处的波数大小，计算公式为：

其中，k_xm，k_ym分别是波数最大处的横坐标与纵坐标。

进一步的，步骤四中计算区分阴影区与非阴影区的阈值τ_S包括：

步骤4.2.1：图像边缘检测，通过对原始雷达图像在相邻t个不同方向上分别基于差分算子进行卷积运算，得到边缘梯度图像I_Ei(r,θ)：

式中，I(r,θ)为原始雷达图像，r为距离，θ为方位角，

表示卷积运算，I_Ei(r,θ)为经过不同的差分算子卷积运算后得到的边缘梯度图像，H_i(r,θ)中i＝1,2,3...t，表示t个不同方向上的像素点差分算子；

步骤4.2.2：将差分图像I_Ei(r,θ)中强度值大小在前N％像素点做为阴影区和非阴影区的交界处像素点，得到子边缘图像I_Ti(r,θ)，其中i＝1,2,3...t，具体公式为：

步骤4.2.3：将I_Ti(r,θ)进行叠加，得到原始边缘图像I_T(r,θ)，其中i＝1,2,3...t：

再利用边缘比较法将原始边缘图像I_T(r,θ)中的每个点与阈值τ_F比较来剔除孤立噪声点，得到边缘图像I_F(r,θ)：

步骤4.2.4：对于边缘图像中每一个孤立的1，找到在该位置处原始雷达图像的像素η，通过边缘图像与原始雷达图像之间的对应关系建立一个图像直方图函数F_H(η)，对F_H(η)取众数即可得到区分阴影区与非阴影区的阈值τ_S：

τ_S＝mode(F_H(η))。

进一步的，步骤4中根据τ_S计算各子区域的阴影比例函数包括：

步骤4.3.1：利用阈值τ_S将原始雷达图像I_T(r,θ)分割成阴影区和非阴影区，阴影图像I_s(r,θ)为：

其中，I(r,θ)为原始雷达图像，I_s(r,θ)为分割后的阴影图像；

步骤4.3.2：计算各子区域的掠射角γ_i，子区域中心点对应的掠射角即为该子区域的掠射角，假设雷达天线高度为H，雷达图像距离分辨率离为R，起始位置基准点数为h₀，第i个子区域起始相对点数为p_i，末端相对点数为q_i，其中i＝1,2,3...m，则第i个子区域掠射角的计算公式为：

步骤4.3.3：计算各子区域的阴影比例函数S(γ_i)，其中i＝1,2,3...m，则第i个子区域阴影比例函数的计算公式为：

式中，P_i表示第i个子区域对应的阴影图像I_s(r,θ)中像素值为1的像素点个数，Q_i表示第i个子区域对应的阴影图像I_s(r,θ)中所有的像素点个数。

进一步的，步骤4中根据阴影比例函数计算得到均方根波陡包括：

步骤4.4.1：同一区域的阴影比例函数S(γ_i)和照明比L(γ_i)满足：

L(γ_i)＝1-S(γ_i)

利用最小二乘法，通过Smith函数拟合各子区域中照明比L(γ_i)与掠射角γ_i之间的关系，其中i＝1,2,3...m，估计得到均方根波陡σ_RMS，Smith函数为：

式中，

μ_i＝tan(γ_i)

其中，erfc( )是高斯误差函数，erfc( )函数是由高斯曲面概率密度函数推导得到，高斯曲面概率密度函数的表达式为：

式中，ξ是高斯海平面距离水平面的高度，υ为高度偏差的均方根，Δξ为水平轴方向增加Δτ下高斯面的高度增量。

进一步的，航海雷达为X波段航海雷达。

本发明的有益效果：本发明设计了一种基于航海雷达图像，特别是X波段航海雷达图像的将3D-FFT谱分析法和阴影统计法相结合的有效波高组合反演算法，弥补了3D-FFT谱分析法算法计算时间长、工程应用成本高和阴影统计法中波高参数计算依赖浮标数据，即需要外部参考设备进行标定的不足，增加算法的适用性和独立性。

与现有技术相比，其优点在于：

(1)本发明利用阴影法计算海浪的平均波陡，通过平均过零周期和平均波陡参数计算海浪的有效波高，不需要浮标实验来标定相关的系数，算法时间成本低，能够降低研发成本和缩短研发周期。

(2)本发明通过对谱分析法和阴影统计法两种算法的组合，能独立地通过雷达数据反演海浪的有效波高参数，使得利用雷达测量海浪有效波高参数的过程完全脱离了浮标的支持，不再需要借助于其他仪器，从而提高了算法的工程适用性和降低了额工程成本。

附图说明

图1是雷达波峰周期与浮标平均过零周期拟合图；

图2是3D-FFT谱分析法研究区域选择示意图；

图3是3D-FFT谱分析法笛卡尔框选择示意图；

图4是3D-FFT谱分析法获取海浪参数流程图；

图5是有效波高2.0m对应的二维海浪谱；

图6是平均过零周期折线图；

图7是阴影法区域分割示意图；

图8是有效波高折线图；

图9(a)和图9(b)是雷达反演有效波高与浮标有效波高比较图；

图10是本发明流程图。

具体实施方式

下面结合说明书附图和实施例对本发明做进一步说明。

结合图10，本发明包括以下步骤：

步骤1，离线数据处理，通过最小二乘法拟合得到海浪的波峰周期和平均过零周期之间的关系，包括以下步骤：

步骤1.1，利用离线的浮标数据和雷达，通过最小二乘法拟合得到海浪的波峰周期和平均过零周期之间的关系。假设两者之间线性关系，并记对应的线性关系为：

T_m02＝aT_M+b

式中，T_M为波峰周期，T_m02为平均过零周期。

记离线实验中浮标的平均过零周期序列为T_m02i，i＝1,2,3...n，对应的通过3D-FFT谱分析法利用雷达数据得到的波峰周期序列记为T_Mi，i＝1,2,3...n，由最小二乘法的相关知识有：

式中，

和

分别是平均过零周期序列和波峰周期序列的平均值。

评价线性关系的最常用参数为相关系数r，r的取值范围是[0,1]，其反映了两个序列之间线性关系的密切程度，r越接近1表示拟合效果越好，相反则表示拟合效果越差。r的定义为：

步骤2，通过3D-FFT谱分析法得到二维海浪谱，从二维海浪谱中获取海浪的波峰周期，包括以下步骤：

步骤2.1，选择的研究区域大小为L_x*L_y，雷达的空间分辨率为R，从而研究区域的像素大小为M_x*M_y。

步骤2.2，通过3D-FFT谱分析法得到二维海浪谱，从二维海浪谱中获取海浪的波峰周期。3D-FFT谱分析法分为图像预处理，3D-FFT，色散关系估流，带通滤波、积分，经验调制传递函数，计算波峰周期等步骤，所述步骤2.2包括以下步骤：

步骤2.2.1，图像预处理：通过最近点插值，将原始雷达图像中极坐标下的回波强度值转换到笛卡尔坐标系下对应的点上。记原始雷达图像中某一点的坐标为(r₀,θ₀)，其转换后对应的笛卡尔坐标为(x₀,y₀)，插值后得到的雷达图像序列记为η(x,y,t)，转换公式为：

式中，r₀为原始雷达图像中某一点的径向距离，θ₀为原始雷达图像中某一点的方位角，x₀,y₀分别为笛卡尔坐标系的横坐标值和纵坐标值，sqrt( )、round( )和rem( )分别是开平方函数、取整函数、和求余函数。

步骤2.2.2，对由步骤2.2.1插值得到的雷达图像序列η(x,y,t)进行三维傅里叶变换得到三维图像谱F(k_x,k_y,ω)，变换公式为：

选取的研究区域大小为L_x*L_y的笛卡尔框，L_x,L_y为笛卡尔框沿x和y方向的范围，雷达图像序列的总时间为T，则上述变换公式可以写成：

利用变换公式的离散形式将图像序列从时间空间域转化到波数频率域，其离散形式为：

式中，k_x和k_y分别为x和y方向的波数，ω为角频率，则分辨率的计算公式为：

式中，L_x,L_y是所选研究的笛卡尔框沿x和y方向的范围，T为雷达图像序列的总时间。

步骤2.2.3，色散关系估流：假设波数处对应的理论频率为ω₀，真实频率ω通过色散关系表示，则两者之间插值的平方和可表示为：

将三维图像谱中的能量值作为权重，从而得到目标函数为：

当SSE关于的偏导数均为零时，SSE有最小解：

从而可以解算出流速为：

式中，h为水深，g为重力加速度，F为由步骤2.2.2得到的三维图像谱F(k_x,k_y,ω)，海表层流

波数

其中，k_x,k_y分别为x和y方向的波数。

步骤2.2.4，带通滤波、积分：利用基于色散关系的带通滤波器对由步骤2.2.2得到的三维图像谱进行滤波，采用的滤波器是基于最大流速的带通滤波器：

式中，E(k_x,k_y,ω)为滤波后得到的三维图像谱，F(k_x,k_y,ω)为由步骤2.2.2得到的三维图像谱，其中，

式中，k为波数，dω为频率分辨率，dk＝(dk_x＝dk_y)为波数分辨率，U_max为最大流速，当B_n取值为负时，则令B_n＝0。

其中，

式中，L_x,L_y是所选研究的笛卡尔框沿x和y方向的范围，T为雷达图像序列的总时间。

再对步骤2.2.4滤波后得到的三维图像谱E(k_x,k_y,ω)进行频率ω维积分得到二维图像谱I(k_x,k_y)，积分公式为：

I(k_x,k_y)＝∫_ω＞0E(k_x,k_y,ω)dω

步骤2.2.5，经验调制传递函数：通过经验调制传递函数将由步骤2.2.4得到的二维图像谱I(k_x,k_y)转换为二维海浪谱E(k_x,k_y)，转换公式为：

E(k_x,k_y)＝|M(k_x,k_y)|²·I(k_x,k_y)

经验调制传递函数是通过实验数据拟合得到，目前使用较为普遍的经验调制传递函数的形式为：

|M(k_x,k_y)|²≈k^-β

其中，β为经验系数。

步骤2.2.6，计算波峰周期：通过由步骤2.2.5得到的二维海浪谱E(k_x,k_y)利用色散关系计算海浪的波峰周期T_M，计算公式为：

式中，g为重力加速度，k_m是海浪谱E(k_x,k_y)中能量值最大处的波数大小，计算公式为：

其中，k_xm，k_ym分别是波数最大处的横坐标与纵坐标。

步骤3，计算平均过零周期，包括以下步骤：

步骤3.1，由步骤1可得到海浪波峰周期和平均过零周期之间的关系式，由步骤2通过3D-FFT谱分析法计算可得到海浪的波峰周期，将波峰周期带入关系式计算平均过零周期。

步骤4，通过阴影法计算海浪的均方根波陡，包括以下步骤：

步骤4.1，选择研究区域，步骤4.1包括以下步骤：

步骤4.1.1，在原始雷达图像中选取方位角度范围为θ₁-θ₂，径向距离范围为r₁-r₂的扇形区域作为研究区域，在选定的研究区域中将原始雷达图像径向维等分为m个子区域。

步骤4.2，计算区分阴影区与非阴影区的阈值，步骤4.2包括以下步骤：

步骤4.2.1，图像边缘检测，通过对原始雷达图像在相邻t个不同方向上分别基于差分算子进行卷积运算，得到边缘梯度图像I_Ei(r,θ)，运算公式为：

式中，I(r,θ)为原始雷达图像，r为距离，θ为方位角，

表示卷积运算，I_Ei(r,θ)为经过不同的差分算子卷积运算后得到的边缘梯度图像，H_i(r,θ)中i＝1,2,3...t，表示t个不同方向上的像素点差分算子。

步骤4.2.2，将由步骤4.2.1.得到的差分图像I_Ei(r,θ)中强度值大小在前N％像素点做为阴影区和非阴影区的交界处像素点，得到子边缘图像I_Ti(r,θ)，其中i＝1,2,3...t，具体公式为：

步骤4.2.3，将由步骤4.2.2得到的子边缘图像I_Ti(r,θ)进行叠加，得到原始边缘图像I_T(r,θ)，其中i＝1,2,3...t：

再利用边缘比较法将原始边缘图像I_T(r,θ)中的每个点与阈值τ_F比较来剔除孤立噪声点，得到边缘图像I_F(r,θ)，具体公式为：

步骤4.2.4，对于边缘图像中每一个孤立的1，找到在该位置处原始雷达图像的像素η，通过边缘图像与原始雷达图像之间的对应关系建立一个图像直方图函数F_H(η)，对F_H(η)取众数即可得到阴影阈值τ_S。

τ_S＝mode(F_H(η))

步骤4.3，计算各子区域的阴影比例函数，步骤4.3包括以下步骤：

步骤4.3.1，利用步骤4.2求得的阴影阈值τ_S将原始雷达图像I_T(r,θ)分割成阴影区和非阴影区，阴影图像I_s(r,θ)为：

其中，I(r,θ)为原始雷达图像，I_s(r,θ)为分割后的阴影图像。

步骤4.3.2，计算各子区域的掠射角γ_i，子区域中心点对应的掠射角即为该子区域的掠射角。假设雷达天线高度为H，雷达图像距离分辨率离为R，起始位置基准点数为h₀，第i个子区域起始相对点数为p_i，末端相对点数为q_i，其中i＝1,2,3...m，则第i个子区域掠射角的计算公式为：

步骤4.3.3，计算各子区域的阴影比例函数S(γ_i)，即计算各子区域对应步骤4.2.2得到的阴影图像中像素值为1的像素点个数占其所在子区域内所有像素点个数的比例，其中i＝1,2,3...m，则第i个子区域阴影比例函数的计算公式为：

式中，P_i表示第i个子区域对应的阴影图像I_s(r,θ)中像素值为1的像素点个数，Q_i表示第i个子区域对应的阴影图像I_s(r,θ)中所有的像素点个数。

步骤4.4，计算均方根波陡，步骤4.4包括以下步骤：

步骤4.4.1，同一区域的阴影比例函数S(γ_i)和照明比L(γ_i)存在以下关系：

L(γ_i)＝1-S(γ_i)

利用最小二乘法，通过Smith函数拟合各子区域中照明比L(γ_i)与掠射角γ_i之间的关系，其中i＝1,2,3...m，估计得到均方根波陡σ_RMS，Smith函数为：

式中，

μ_i＝tan(γ_i)

其中，erfc( )是高斯误差函数，erfc( )函数是由高斯曲面概率密度函数推导得到，高斯曲面概率密度函数的表达式为：

式中，ξ是高斯海平面距离水平面的高度，υ为高度偏差的均方根，Δξ为水平轴方向增加Δτ下高斯面的高度增量。

步骤5，通过平均过零周期和均方根波陡参数计算海浪的有效波高，步骤5包括以下步骤：

步骤5.1，利用有效波高与均方根波陡、平均过零周期之间的关系计算得到有效波高，具体公式为：

式中，g为重力加速度，H_S为有效波高，σ_RMS为均方根波陡，T_m02为平均过零周期。

下面结合具体参数给出实施例。

本发明具体可以分为以下几步：第一步为离线数据处理，通过最小二乘法拟合得到海浪波峰周期和平均过零周期之间的关系，第二步为通过3D-FFT谱分析法得到二维海浪谱，从二维海浪谱中获取海浪的波峰周期，第三步为计算平均过零周期，第四步为通过阴影法计算海浪的均方根波陡，第五步为通过平均过零周期和均方根波陡参数计算海浪的有效波高。

本发明实施所用的航海雷达为X波段导航雷达，工作于短脉冲模式，回波数据数字化后以极坐标形式按线存储，两条相邻存储线间的时间间隔小于1ms，雷达天线扫描一周的时间约2.5s，一幅雷达图像的总线数为2048条，每根线上有2048个像元点，其方位向分辨率约为0.1°，径向分辨率约为2.5m。

本发明实例所用航海雷达的主要技术参数如表一所示：

表一 航海雷达的技术参数

结合图10，本发明具体实验步骤为：

第一步，离线数据处理，通过最小二乘法拟合得到海浪的波峰周期和平均过零周期之间的关系，包括以下步骤：

步骤1.1，利用离线的浮标数据和雷达，通过最小二乘法拟合海浪的波峰周期和平均过零周期之间的关系。图1是本例雷达波峰周期与浮标平均过零周期拟合图，图中横坐标为雷达计算得到的海浪波峰周期，纵坐标为浮标测量的海浪平均过零周期，拟合得到的关系式为：

T_m02＝0.6136T_M+2.538

式中，T_M为波峰周期，T_m02为平均过零周期，本实例中拟合结果的相关系数为0.75，这表明线性模型能较好地反映两者之间的关系。

第二步，通过3D-FFT谱分析法得到二维海浪谱，从二维海浪谱中获取海浪的波峰周期，包括以下步骤：

步骤2.1，选择的研究区域大小为L_x*L_y，雷达的空间分辨率为R，从而研究区域的像素大小为M_x*M_y。本实例中，L_x＝L_y＝960m，R＝7.5m，M_x＝M_y＝128，图2为研究区域选择示意图，图3为笛卡尔框选择示意图，图3中黑色区域为选取的笛卡尔框。

步骤2.2，通过3D-FFT谱分析法得到二维海浪谱，从二维海浪谱中获取海浪的波峰周期，3D-FFT谱分析法获取海浪参数流程图见图4，步骤2.2包括以下步骤：

步骤2.2.1，通过最近点插值，将原始雷达图像中极坐标下的回波强度值转换到笛卡尔坐标系下对应的点上。原始雷达图像中某一点的坐标记为(r₀,θ₀)，其转换后对应的笛卡尔坐标记为(x₀,y₀)，插值后得到的雷达图像序列记为η(x,y,t)，转换公式为：

式中，r₀为原始雷达图像中某一点的径向距离，θ₀为原始雷达图像中某一点的方位角，x₀,y₀分别为笛卡尔坐标系的横坐标值和纵坐标值，sqrt( )、round( )和rem( )分别是开平方函数、取整函数、和求余函数。

步骤2.2.2，对由步骤2.2.1最近点插值得到的雷达图像序列η(x,y,t)进行三维傅里叶变换得到三维图像谱F(k_x,k_y,ω)，变换公式为：

选取的研究区域大小为L_x*L_y的笛卡尔框，L_x,L_y是所选研究的笛卡尔框沿x和y方向的范围，雷达图像序列的总时间为T，则上述变换公式可以写成：

利用变换公式的离散形式将图像序列从时间空间域转化到波数频率域，其离散形式为：

式中，k_x和k_y分别为x和y方向的波数，ω为角频率，则分辨率的计算公式为：

其中，L_x,L_y是所选研究的笛卡尔框沿x和y方向的范围，T为雷达图像序列的总时间。本实例中，L_x＝L_y＝960m，T＝90s。

步骤2.2.3，基于色散原理和最小二乘原理估算海流信息，海流的计算公式为：

式中，ω₀为波数处对应的理论频率，h为水深，g为重力加速度，F为由步骤2.2.2得到的三维图像谱F(k_x,k_y,ω)，海表层流

波数

其中，k_x,k_y分别为x和y方向的波数。

步骤2.2.4，利用基于色散关系的带通滤波器对由步骤2.2.2得到的三维图像谱进行滤波，采用的滤波器是基于最大流速的带通滤波器：

式中，E(k_x,k_y,ω)为滤波后得到的三维图像谱，F(k_x,k_y,ω)为由步骤2.2.2得到的三维图像谱，其中，

式中，k为波数，dω为频率分辨率，dk＝(dk_x＝dk_y)为波数分辨率，U_max为最大流速，当B_n取值为负时，则令B_n＝0。

其中，

式中，L_x,L_y是所选研究的笛卡尔框沿x和y方向的范围，T为雷达图像序列的总时间。本实例中，L_x＝L_y＝960m，T＝90s。

再对步骤2.2.4滤波后得到的三维谱E(k_x,k_y,ω)进行频率ω维积分得到二维图像谱I(k_x,k_y)，积分公式为：

I(k_x,k_y)＝∫_ω＞0E(k_x,k_y,ω)dω

步骤2.2.5，通过经验调制传递函数将将由步骤2.2.4得到的二维图像谱I(k_x,k_y)转换为二维海浪谱E(k_x,k_y)，转换公式为：

E(k_x,k_y)＝|M(k_x,k_y)|²·I(k_x,k_y)

经验调制传递函数是通过实验数据拟合得到，目前使用较为普遍的经验调制传递函数的形式为：

|M(k_x,k_y)|²≈k^-β

其中，β为经验系数，本实例中β取值为1.2。

步骤2.2.6，计算波峰周期，通过由步骤2.2.5得到的二维海浪谱E(k_x,k_y)利用色散关系计算海浪的波峰周期T_M，计算公式为：

式中，g为重力加速度，k_m是海浪谱E(k_x,k_y)中能量值最大处的波数大小，计算公式为：

其中，k_xm，k_ym分别是波数最大处的横坐标与纵坐标。

本实例选取了有效波高为2.0m的一组数据，然后用3D-FFT谱分析法处理得到相应的二维海浪谱，结果见图5，经计算得到其对应的波峰周期为8.21s，从图中可以看出，二维海浪谱能较好地表征出海浪的波峰周期。

第三步，计算平均过零周期，包括以下步骤：

步骤3.1，由第一步可得到海浪波峰周期和平均过零周期之间的关系，由第二步通过3D-FFT谱分析法计算可得到海浪的波峰周期，将波峰周期带入关系式计算平均过零周期，计算结果详见图6，图中实线为浮标输出的平均过零周期，虚线是利用雷达数据通过3D-FFT谱分析算法得到的海浪平均过零周期。本实例选取了有效波高为2.0m的一组数据，经第二步计算得到其对应的波峰周期为8.21s，将波峰周期带入关系式计算得到其对应的平均过零周期为7.57s。

第四步，通过阴影法计算海浪的均方根波陡，包括以下步骤：

步骤4.1，选择研究区域，步骤4.1包括以下步骤：

步骤4.1.1，在原始雷达图像中选取方位角度范围为θ₁-θ₂，径向距离范围为r₁-r₂的扇形区域作为研究区域，在选定的研究区域中将原始雷达图像径向维等分为m个子区域。本实例中，m＝50，θ₁＝120°,θ₂＝230°,r₁＝200m,r₂＝2500m,阴影法区域分割示意图见图7。

步骤4.2，计算区分阴影区与非阴影区的阈值，步骤4.2包括以下步骤：

步骤4.2.1，图像边缘检测，通过对原始雷达图像在相邻t个不同方向上分别基于差分算子进行卷积运算，得到边缘梯度图像I_Ei(r,θ)，运算公式为：

式中，I(r,θ)为原始雷达图像，r为距离，＝为方位角，

表示卷积运算，I_Ei(r,θ)为经过不同的差分算子卷积运算后得到的边缘梯度图像，H_i(r,θ)表示t个不同方向上的像素点差分算子，其中，i＝1,2,3...t，本实例中t＝8，具体为：

步骤4.2.2，将由步骤4.2.1.得到的差分图像I_Ei(r,θ)中强度值大小在前N％像素点做为阴影区和非阴影区的交界处像素点，得到子边缘图像I_Ti(r,θ)，其中i＝1,2,3...t，具体公式为：

本实例中，t＝8，N取值为10。

步骤4.2.3，将由步骤4.2.2得到的子边缘图像I_Ti(r,θ)进行叠加，其中i＝1,2,3...t，得到原始边缘图像I_T(r,θ)：

本实例中，t＝8。

再利用边缘比较法将原始边缘图像I_T(r,θ)中的每个点与阈值τ_F比较来剔除孤立噪声点，得到边缘图像I_F(r,θ)，具体公式为：

本实例中，τ_F取值为5。

步骤4.2.4，对于边缘图像中每一个孤立的1，找到在该位置处原始雷达图像的像素η，通过边缘图像与原始雷达图像之间的对应关系建立一个图像直方图函数F_H(η)，对F_H(η)取众数即可得到阴影阈值τ_S。

τ_S＝mode(F_H(η))

步骤4.3，计算各子区域的阴影比例函数，步骤4.2包括以下步骤：

步骤4.3.2，利用步骤4.1求得的阴影阈值τ_S将原始雷达图像I_T(r,θ)分割成阴影区和非阴影区，阴影图像I_s(r,θ)为：

其中，I(r,θ)为原始雷达图像，I_s(r,θ)为分割后的阴影图像。

步骤4.3.3，计算子区域的掠射角γ_i，子区域中心点对应的掠射角即为该子区域的掠射角。假设雷达天线高度为H，雷达图像距离分辨率离为R，起始位置基准点数为h₀，第i个子区域起始相对点数为p_i，末端相对点数为q_i，i＝1,2,3...m，则第i个子区域掠射角的计算公式为：

本实例中，雷达天线高度H为45m，雷达图像距离分辨率R为7.5m，起始位置基准点数h₀为25，m＝50。

步骤4.3.4，计算各子区域的阴影比例函数S(γ_i)，即计算各子区域对应步骤4.2.2得到的阴影图像中像素值为1的像素点个数占其所在子区域内所有像素点个数的比例，i＝1,2,3...m，则第i个子区域阴影比例函数的计算公式为：

式中，P_i表示第i个子区域对应的阴影图像I_s(r,θ)中像素值为1的像素点个数，Q_i表示第i个子区域对应的阴影图像I_s(r,θ)中所有的像素点个数，本实例中m＝50。

步骤4.4，计算均方根波陡，步骤4.4包括以下步骤：

步骤4.4.1，同一区域的阴影比例函数S(γ_i)和照明比L(γ_i)存在以下关系：

L(γ_i)＝1-S(γ_i)

利用最小二乘法，通过Smith函数拟合各子区域中照明比L(γ_i)与掠射角γ_i之间的关系，i＝1,2,3...m，本实例中m＝50，估计得到均方根波陡σ_RMS，Smith函数为：

式中，

μ_i＝tan(γ_i)

其中，erfc( )是高斯误差函数，erfc( )函数是由高斯曲面概率密度函数推导得到，高斯曲面概率密度函数的表达式为：

式中，ξ是高斯海平面距离水平面的高度，υ为高度偏差的均方根，Δξ为水平轴方向增加Δτ下高斯面的高度增量。

第五步，通过平均过零周期和均方根波陡参数计算海浪的有效波高，包括以下步骤：

步骤5.1，利用有效波高与均方根波陡、平均过零周期之间的关系计算得到有效波高，具体公式为：

式中，g为重力加速度，H_S为有效波高，σ_RMS为均方根波陡，T_m02为平均过零周期。

本实例选取了有效波高为2.0m的一组数据，经步骤2计算得到其对应的波峰周期为8.21s，经步骤3将波峰周期带入步骤1拟合得到的关系式计算其对应的平均过零周期为7.57s，经步骤4计算得到其均方根波陡为2.38，本步骤5计算得到其有效波高为2.80m。

应用本发明提出的一种基于X波段航海雷达图像的有效波高反演方法，对2014年东海海域获取的大量雷达数据及相关时间段的海况信息进行实验分析。本实验验选取2014年11月16日18时-2014年11月19日0时在福建平潭采集到的岸基数据进行实验。分别利用阴影法、谱分析法和阴影统计法组合算法进行对比实验，阴影法、谱分析法和阴影统计法组合算法有效波高序列计算结果见图8。图中分别是浮标、阴影法和组合算法得到的有效波高。从图8中可以看出，以浮标输出的有效波高数据作为真值，阴影法和组合算法的计算出的有效波高数据相对整体偏高。图9(a)和图9(b)是通过阴影法和组合算法得到的有效波高与浮标有效波高之间比较的散点图，其中图9(a)是阴影法有效波高与浮标有效波高的散点图，图9(b)是组合算法有效波高与浮标有效波高的散点图。

在本发明实例中，首先分别将由雷达3D-FFT分析方法和浮标各自输出的平均过零周期T_m02进行比较分析，以浮标的波周期作为参考真值，并采用平均偏差BIAS和相关系数CORR作为评价算法的参数。假设X_i为该发明组合算法计算出的有效波高序列，Y_i是浮标输出的有效波高真值，其中i＝1,2,...N，则有：

通过以上公式计算可得表二：

表二 算法性能参数表

实验结果表明，组合算法的性能虽然相对于阴影法有所下降，但是组合算法能独立地利用雷达数据得到海浪的有效波高参数，使得利用雷达测量海浪有效波高参数的过程完全脱离了浮标的支持，不再需要借助于其他仪器，增加了算法的适用性和独立性。

Claims (7)

1.一种航海雷达图像反演有效波高的方法，其特征在于，包括以下步骤：

步骤1：利用离线的浮标数据和雷达，通过最小二乘法拟合得到海浪的波峰周期和平均过零周期之间的关系式；

步骤2：选定研究区域，通过3D-FFT谱分析法得到选定研究区域的二维海浪谱，从二维海浪谱中获取海浪的波峰周期；

步骤3：将步骤2得到的海浪的波峰周期代入步骤1得到的海浪波峰周期和平均过零周期之间的关系式，计算得到平均过零周期T_m02；

步骤4：在原始雷达图像中选取研究区域，在选定的研究区域中将原始雷达图像径向维等分为m个子区域；计算区分阴影区与非阴影区的阈值τ_S；根据τ_S计算各子区域的阴影比例函数；然后根据阴影比例函数计算得到均方根波陡σ_RMS；

步骤5：通过平均过零周期T_m02和均方根波陡σ_RMS计算海浪的有效波高H_S，H_S满足：

式中，g为重力加速度。

2.根据权利要求1所述的一种航海雷达图像反演有效波高的方法，其特征在于：所述海浪的波峰周期和平均过零周期之间的关系具体为：

T_m02＝aT_M+b

式中，T_M为波峰周期，T_m02为平均过零周期；

记离线实验中浮标的平均过零周期序列为T_m02i，i＝1,2,3...n，对应的通过3D-FFT谱分析法利用雷达数据得到的波峰周期序列记为T_Mi，i＝1,2,3...n，根据最小二乘法得到：

式中，

和

分别是平均过零周期序列和波峰周期序列的平均值。

3.根据权利要求1所述的一种航海雷达图像反演有效波高的方法，其特征在于：所述步骤2包括：

步骤2.1：通过最近点插值，将原始雷达图像中极坐标下的回波强度值转换到笛卡尔坐标系下对应的点上，记原始雷达图像中某一点的坐标为(r₀,θ₀)，其转换后对应的笛卡尔坐标为(x₀,y₀)，插值后得到的雷达图像序列记为η(x,y,t)，转换公式为：

式中，r₀为原始雷达图像中某一点的径向距离，θ₀为原始雷达图像中某一点的方位角，x₀,y₀分别为笛卡尔坐标系的横坐标值和纵坐标值，sqrt()、round()和rem()分别是开平方函数、取整函数、和求余函数；

步骤2.2：选取的研究区域大小为L_x*L_y的笛卡尔框，L_x,L_y为笛卡尔框沿x和y方向的范围，雷达图像序列的总时间为T，对雷达图像序列η(x,y,t)进行三维傅里叶变换得到三维图像谱，然后利用变换公式的离散形式将图像序列从时间空间域转化到波数频率域：

式中，k_x和k_y分别为x和y方向的波数，ω为角频率，则分辨率的计算公式为：

步骤2.3：假设波数处对应的理论频率为ω₀，真实频率ω通过色散关系表示，则两者之间插值的平方和为：

将三维图像谱中的能量值作为权重，得到目标函数为：

当SSE关于的偏导数均为零时，SSE有最小解，解算出流速为：

式中，h为水深，g为重力加速度，F为三维图像谱F(k_x,k_y,ω)，海表层流

波数

其中，k_x,k_y分别为x和y方向的波数；

步骤2.4：利用基于色散关系的带通滤波器对三维图像谱进行滤波，采用的滤波器是基于最大流速的带通滤波器：

式中，E(k_x,k_y,ω)为滤波后得到的三维图像谱，其中，

式中，k为波数，dω为频率分辨率，dk＝(dk_x＝dk_y)为波数分辨率，U_max为最大流速，当B_n取值为负时，则令B_n＝0；

再对滤波后得到的三维图像谱E(k_x,k_y,ω)进行频率ω维积分得到二维图像谱I(k_x,k_y)：

I(k_x,k_y)＝∫_ω＞0E(k_x,k_y,ω)dω

步骤2.5：通过经验调制传递函数将二维图像谱I(k_x,k_y)转换为二维海浪谱E(k_x,k_y)：

E(k_x,k_y)＝|M(k_x,k_y)|²·I(k_x,k_y)

经验调制传递函数为：

|M(k_x,k_y)|²≈k^-β

其中，β为经验系数；

步骤2.6：通过二维海浪谱E(k_x,k_y)利用色散关系计算海浪的波峰周期T_M，计算公式为：

式中，g为重力加速度，k_m是海浪谱E(k_x,k_y)中能量值最大处的波数大小，计算公式为：

其中，k_xm，k_ym分别是波数最大处的横坐标与纵坐标。

4.根据权利要求1所述的一种航海雷达图像反演有效波高的方法，其特征在于：步骤四所述计算区分阴影区与非阴影区的阈值τ_S包括

步骤4.2.1：图像边缘检测，通过对原始雷达图像在相邻t个不同方向上分别基于差分算子进行卷积运算，得到边缘梯度图像I_Ei(r，θ)：

式中，I(r，θ)为原始雷达图像，r为距离，θ为方位角，

表示卷积运算，I_Ei(r，θ)为经过不同的差分算子卷积运算后得到的边缘梯度图像，H_i(r，θ)中i＝1，2，3...t，表示t个不同方向上的像素点差分算子；

步骤4.2.2：将差分图像I_Ei(r，θ)中强度值大小在前N％像素点做为阴影区和非阴影区的交界处像素点，得到子边缘图像I_Ti(r，θ)，其中i＝1，2，3...t，具体公式为：

步骤4.2.3：将I_Ti(r，θ)进行叠加，得到原始边缘图像I_T(r，θ)，其中i＝1，2，3...t：

再利用边缘比较法将原始边缘图像I_T(r，θ)中的每个点与阈值τ_F比较来剔除孤立噪声点，得到边缘图像I_F(r，θ)：

步骤4.2.4：对于边缘图像中每一个孤立的1，找到在该位置处原始雷达图像的像素η，通过边缘图像与原始雷达图像之间的对应关系建立一个图像直方图函数F_H(η)，对F_H(η)取众数即可得到区分阴影区与非阴影区的阈值τ_S：

τ_S＝mode(F_H(η))。

5.根据权利要求1所述的一种航海雷达图像反演有效波高的方法，其特征在于：步骤4所述根据τ_S计算各子区域的阴影比例函数包括：

步骤4.3.1：利用阈值τ_S将原始雷达图像I_T(r，θ)分割成阴影区和非阴影区，阴影图像I_s(r，θ)为：

其中，I(r，θ)为原始雷达图像，I_s(r，θ)为分割后的阴影图像；

步骤4.3.2：计算各子区域的掠射角γ_i，子区域中心点对应的掠射角即为该子区域的掠射角，假设雷达天线高度为H，雷达图像距离分辨率离为R，起始位置基准点数为h₀，第i个子区域起始相对点数为p_i，末端相对点数为q_i，其中i＝1，2，3...m，则第i个子区域掠射角的计算公式为：

步骤4.3.3：计算各子区域的阴影比例函数S(γ_i)，其中i＝1，2，3...m，则第i个子区域阴影比例函数的计算公式为：

式中，P_i表示第i个子区域对应的阴影图像I_s(r，θ)中像素值为1的像素点个数，Q_i表示第i个子区域对应的阴影图像I_s(r，θ)中所有的像素点个数。

6.根据权利要求1所述的一种航海雷达图像反演有效波高的方法，其特征在于：步骤4所述根据阴影比例函数计算得到均方根波陡包括：

步骤4.4.1：同一区域的阴影比例函数S(γ_i)和照明比L(γ_i)满足：

L(γ_i)＝1-S(γ_i)

利用最小二乘法，通过Smith函数拟合各子区域中照明比L(γ_i)与掠射角γ_i之间的关系，其中i＝1，2，3...m，估计得到均方根波陡σ_RMS，Smith函数为：

式中，

μ_i＝tan(γ_i)

其中，erfc()是高斯误差函数，erfc()函数是由高斯曲面概率密度函数推导得到，高斯曲面概率密度函数的表达式为：

式中，ξ是高斯海平面距离水平面的高度，υ为高度偏差的均方根，Δξ为水平轴方向增加Δτ下高斯面的高度增量。

7.根据权利要求1所述的一种航海雷达图像反演有效波高的方法，其特征在于：所述航海雷达为X波段航海雷达。

Images