CN104156629A - 一种基于相对辐射校正的导航雷达图像反演海面风向方法 - Google Patents

Abstract

本发明的目的在于提供一种基于相对辐射校正的导航雷达图像反演海面风向方法，包含导航雷达图像序列采集、导航雷达图像预处理、导航雷达图像相对辐射较正和海面风向反演四个部分。本发明在应用风条纹反演风向中增添导航雷达图像相对辐射校正环节，有效消除了导航雷达回波强度径向衰减对风向反演造成的影响。导航雷达图像相对辐射校正采用一种自适应拉格朗日最小二乘分段拟合校正法，既保证校正后不会破坏海杂波图像特征，又提高了工程的适用性。

Description

一种基于相对辐射校正的导航雷达图像反演海面风向方法

技术领域

本发明涉及的是一种反演海面风向方法。

背景技术

海面风场信息是海洋动力学重要参数，是海洋与大气能量和气体交换的主要驱动力。因此，了解和掌握海面风场信息对渔业、海运及气象监测都有深远的意义。现有海面风向信息获取方法主要分成两类：站点式现场测量和遥感测量。导航雷达是遥感测量手段的一种，因具有不受光线影响、不受天气影响、实时连续反馈、高分辨率和使用便捷等优点，成为现阶段提取海面风向信息的热门课题。

现阶段国内外应用导航雷达图像测量海面风向，主要方法为基于海面风速分布不均匀导致的导航雷达图像中呈现出的风条纹来计算的。风条纹存在尺度200～500m，频率接近静态或准静态，风条纹方向与风向平行等特征。目前，应用风条纹特征反演海面风向的方法主要有以下几种：2003年Dankert等人提出局部梯度法，2004年由Dankert等人提出的光流法。2010年中国海洋大学段华敏使用和Dankert2004年相同的方法提取出了海表面风场信息。2012年哈尔滨工程大学的贾瑞才博士和2013年中国海洋大学的硕士生李金凤都对Dankert2004年的方法进行了改进，提高了风向反演精度。以上方法虽都已提取出海面风向信息，但应用以上方法都没有考虑导航雷达回波强度径向衰减对风向提取造成的影响。

导航雷达径向回波强度存在非线性衰减的特征会影响风条纹的分布特征，实测导航雷达一条线上径向回波强度分布如图12。Maurice W.Long指出在均匀粗糙海面的回波强度在近距离和远距离分别正比于R-3和R-7。在近区时，充满天线波束照射海面区域内的回波强度服从R-3的关系；在远区时，由于入射角减少，出现海浪反射的双路径效应，造成回波强度与探测距离成R-7的关系；因此，理论分析和实验都证明了导航雷达径向回波强度存在非线性衰减的特征。通过实验还发现径向衰减特征是由100～500m不同波长的波组成，而风条纹的尺度为200～500m，并且两者都属于静态信号，两种特征相近的信号是无法分离的，因此，导航雷达图像径向衰减特征会影响风条纹的分布特征，从而降低风向信息提取的精度。

现有导航雷达图像相对辐射校正方法分为两大类：非线性校正法和线性回归法。非线性校正法主要有直方图匹配法，此方法是对需要校正图像和参考图像作直方图均衡化，然后用参考图像的直方图均衡化去匹配需要校正图形的。此方法需要存在参考图像，匹配过程中图像的特征信号会被淹没。

线性回归校正是通过对样本点建立线性模型来进行图像校正的，具有精度高、速度快和易于工程实现等优点。

传统线性回归法主要分为三步实现，第一步选取适当样本点，第二步应用样本点确定线性方程回归参数建立线性拟合方程，第三步根据线性拟合方程，对图像进行相对辐射校正。第一步样本点选取，现在主要有图像回归法(IR)、伪不变特征法(PIF)、暗级-亮级法和未变化级辐射归一法(NC)。第二步确定线性拟合方程现阶段有两种方法，第一种是建立二阶最小二乘多项式，第二种为求取多幅导航雷达图像均值。结合导航雷达径向回波具有非线性衰减特征，和风条纹静态特征，发现第二步现有的两种方法都不适用。

发明内容

本发明的目的在于提供提高数据拟合速度，简化拟合方程的复杂度，且不会破坏海杂波图像的一种基于相对辐射校正的导航雷达图像反演海面风向方法。

本发明的目的是这样实现的：

本发明一种基于相对辐射校正的导航雷达图像反演海面风向方法，其特征是：

(1)导航雷达图像序列采集：采集连续N幅导航雷达图像，将其定义为一组导航雷达图像序列，并收集对应时间和位置的实际风向与风速值；

(2)导航雷达图像预处理：对导航雷达图像序列极坐标图像进行归一化，使导航雷达图像序列中的每幅导航雷达图像的线数固定，再对单幅导航雷达图像通过观测去除遮挡区域的影响，只保留海杂波图像区域，应用图像中值滤波抑制同频信号对导航雷达图像的干扰；

(3)导航雷达图像相对辐射校正：对步骤(2)处理得到的导航雷达图像的回波强度值按方位向进行平均，得到关于距离的回波强度均值，对导航雷达图像径向回波强度进行相对辐射校正，得到相对辐射校正图像序列；

(4)海面风向反演：对相对辐射校正图像序列首先利用低通滤波器得到海面静态特征图像，然后基于波数能量谱尺度分离方法得到风条纹能量谱，最终获得海面风向。

本发明还可以包括：

1、导航雷达图像相对辐射校正的具体步骤为：

(1)对导航雷达图像沿方位向进行平均，令导航雷达图像回波强度值为σ(θ_p,r_q)，其中θ_p(p＝1,2,…P)为方位向角度，共P个，r_q(q＝1,2,…Q)为径向距离，共Q个，生成一维径向平均回波强度值y：

$y\left(r_q\right)=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}\mathrm{\σ}({\mathrm{\θ}}_p,r_q)---\left(1\right)$

径向距离r_q定义为x，则x和y形成的数据集为数据点总数为n；

(2)选取x中间位置x₀，将n分成两段数据，个数分别为n₁和n₂，若K为间断点数，则K＝1，生成两个数据集进行分段拟合；

(3)假设数据共有K个间断点，则形成K+1个数据集，即 ${\mathrm{\Ω}}_k={\left\{(x_i^k,y_i^k)\right\}}_{i=1}^{n_k},k=\mathrm{1,2},\·\·\·K+1,$ n_k为第k个数据集的点数， $x_{k-1}{\≤x}_i^k\≤x_k,$ x_k是数据集Ω_k和Ω_k+1的分段点；

(4)拟合函数f(x)的形式如下：

$f\left(x\right)=\left\{\begin{array}{c}{f}_1\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^1{h}_j\left(x\right),1\≤x\≤x_1\\ f_2\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^2{h}_j\left(x\right),x_1\≤x\≤x_2\\ \·\\ \·\\ \·\\ f_k\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^2{h}_j\left(x\right),x_{k-1}\≤x\≤x_k\\ \·\\ \·\\ \·\end{array}\right.$

式中f_k-1(x_k)＝f_k(x_k)，是数据集Ω_k的基函数，基函数选为1，x，x2；m为基函数的个数，即m＝3；为待确定回归系数；为了使拟合误差达到最小且在x_k处连续，建立最小二乘回归模型如下：

$\begin{array}{c}\underset{{\mathrm{\α}}_j^k}{\min }\underset{i=1}{\overset{n_1}{\mathrm{\Σ}}}{(f_1\left(x_i^1\right)-y_i^1)}^2+\underset{i=1}{\overset{n_2}{\mathrm{\Σ}}}{(f_2\left(x_i^2\right)-y_i^2)}^2+...+\underset{i=1}{\overset{n_k}{\mathrm{\Σ}}}{(f_k\left(x_i^k\right)-y_i^k)}^2+...\underset{i=1}{\overset{n_k}{\mathrm{\Σ}}}{(f_{K+1}\left(x_i^{K+1}\right)-y_i^{K+1})}^2,\\ s.t.f_{k-1}\left(x_k\right)=f_k\left(x_k\right)\end{array}$

令基函数h(x)＝[1 x x²]，拟合参数行向量分别为： $X_k=\left[\begin{array}{c}h\left(x_1^k\right)\\ h\left(x_2^k\right)\\ \·\\ \·\\ \·\\ \text{h}\left(x_{n_k}^k\right)\end{array}\right],$ $y_k=\left[\begin{array}{c}{y}_1^k\\ y_2^k\\ \·\\ \·\\ \·\\ y_{n_k}^k\end{array}\right],$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_1^k\\ {\mathrm{\α}}_2^k\\ {\mathrm{\α}}_3^k\end{array}\right],$ k＝1,2,…K+1；则多段数据集拟合参数可以表示为：

$y=\left[\begin{array}{c}{y}_1^T\\ y_2^T\\ \·\\ \·\\ \·\\ y_{K+1}^T\end{array}\right]$ $\mathrm{\α}=\left[\begin{array}{c}{\mathrm{\α}}_1^T\\ {\mathrm{\α}}_2^T\\ \·\\ \·\\ \·\\ {\mathrm{\α}}_{K+1}^T\end{array}\right]$ $X=\left[\begin{array}{cccc}{X}_1& 0& ...& 0\\ 0& X_2& ...& 0\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ 0& 0& ...& X_{K+1}\end{array}\right]$

$Z=\left[\begin{array}{cccccc}h\left(x_1\right)& -h\left(x_1\right)& 0& ...& 0& 0\\ 0& h\left(x_2\right)& -h\left(x_2\right)& ...& 0& 0\\ \·& \·& \·& & \·& \·\\ \·& \·& \·& & \·& \·\\ \·& \·& \·& & \·& \·\\ 0& 0& 0& ...& -h\left(x_{K-1}\right)& 0\\ 0& 0& 0& ...& h\left(x_K\right)& -h\left(x_K\right)\end{array}\right],$

则最小二乘回归模型的二次范数形式为：

$\underset{\mathrm{\α}}{\min }{\left|\right|\mathrm{X\α}-y\left|\right|}^2,s.t.\mathrm{Z\α}=0,$

对上式建立拉格朗日函数为：

$L(\∂,\mathrm{\λ})={\left|\right|\mathrm{X\α}-Y\left|\right|}^2+2{\mathrm{\λ}}^T\mathrm{Z\α},$

其中λ为列向量约束参数，即λ＝[λ₁, λ₂, … ,λ_K]^T，由多元函数极值必要条件，得到:

$\left\{\begin{array}{c}\frac{\∂L}{{\∂\mathrm{\α}}_1}=2X^T(\mathrm{X\α}-y)+2Z^T\mathrm{\λ}=0\\ \frac{\∂L}{\∂\mathrm{\λ}}=\mathrm{Z\α}=0\end{array}\right.,$

由上式可得：

$\mathrm{\α}=\widehat{\mathrm{\α}}-{\left(X^{T}X\right)}^{-1}{Z}^T\mathrm{\λ},$

其中 $\widehat{\mathrm{\α}}={\left(X^{T}X\right)}^{-1}{X}^{T}y,$ 将 $\mathrm{\α}=\widehat{\mathrm{\α}}-{\left(X^{T}X\right)}^{-1}{Z}^T\mathrm{\λ},$ 带入 $\left\{\begin{array}{c}\frac{\∂L}{{\∂\mathrm{\α}}_1}=2X^T(\mathrm{X\α}-y)+2Z^T\mathrm{\λ}=0\\ \frac{\∂L}{\∂\mathrm{\λ}}=\mathrm{Z\α}=0\end{array}\right.,$ 可解得：

$\mathrm{\λ}={\left[Z{\left(X^{T}X\right)}^{-1}{Z}^T\right]}^{-1}Z\widehat{\mathrm{\α}},$

将式带入获得最小二乘回归系数α的解，根据导航雷达回波机制数据集Ω_k可以保证式中的是非退化的，Z是非零的，因此，α的解是唯一的，从而获得有多段拟合函数f(x)；

(5)计算两个数据集绝对值误差均值S_k,k＝1,2…K+1，每点的绝对值误差e_i,i＝1,2,…n，若e_i>S_k&e_i+1>S_k&e_i+2>S_k，则在e_i的处位置再次分段，执行步骤(6)；其他情况或分段内点数少于十个，则不分段，执行步骤(7)；

(6)在分段点处对数据再次进行分段，返回从步骤(3)开始重新拟合，K为分段点的数量，再执行步骤(5)的判断；得到数据集最优分段函数拟合函数值；

(7)对导航雷达图像中每个θ_p的回波强度值σ(θ_p,r_q)按r_q分段的减去拟合函数值，实现导航雷达图像相对辐射校正；校正遍历导航雷达图像序列的每幅图像，最终得到相对辐射校正后图像序列f′(θ,r,t)，θ为方位向角度，r为径向距离，t为图像序列时间轴。

2、海面风向反演的具体步骤为：

(1)构建全局低通滤波器以获得包含风条纹特征的海面静态特征图像：

$f(\mathrm{\θ},r)=\frac{\underset{t=1}{\overset{N_t}{\mathrm{\Σ}}}{f}^{\′}(\mathrm{\θ},r,t)}{N_t}$ f(θ,r)为海面静态特征图像，在笛卡尔坐标下可以表示为 和分别为笛卡尔坐标下图像沿x和y轴的坐标，N_x和N_y为坐标位置，N_t为时间序列；

(2)应用二维离散快速傅里叶变换将海面静态特征图像变换到能量谱域，即：

$F(k_x,k_y)=\underset{N_x}{\mathrm{\Σ}}\underset{N_y}{\mathrm{\Σ}}f(x_{N_x},y_{N_y})\exp [-i*2\mathrm{\π}(\frac{z_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}{*k}_y}{\max \left(y_{N_y}\right)})]$

F(k_x,k_y)为f(x,y)的傅里叶系数，复指数项为：

$\begin{array}{c}\exp [-i*2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})]=\cos 2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})\\ +i\sin 2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})\end{array}$

其中，

$k_x=\frac{2\mathrm{\π}}{\max \left(x_{N_x}\right)*d_x}$

$k_y=\frac{2\mathrm{\π}}{\max \left(y_{N_y}\right)*d_y}$

(k_x,k_y)为在能量谱Τ下的坐标，d_x、d_y为导航雷达图像分辨率；根据二维谱性质得到的海面静态特征图像能量谱为A(k_x,k_y)：

$A(k_x,k_y)=\sqrt{\left[\mathrm{Re}\left(F(k_x,k_y)\right)\right]+\left[\mathrm{Im}\left(F(k_x,k_y)\right)\right]}$

(3)基于尺度分离获得风条纹能量谱：导航雷达图像呈现的海面条纹波长L与能量谱波数k的关系式为：

$k=\frac{2\mathrm{\π}}{L}$

其中，(k_x,k_y)为能量谱域Τ的坐标，即k_x和k_y为k在x和y轴的波数分量；

设L_d为风条纹波长尺度下限，L_t为风条纹波长尺度上限，根据得到风条纹能量谱波数下限为k_d：

$k_d=\frac{2\mathrm{\π}}{L_d}$

风条纹能量谱波数上限为k_t：

$k_t=\frac{2\mathrm{\π}}{L_t}$

根据风条纹能量谱波数上下限设计能量谱带通滤波器，将风条纹信号能量谱提取出来，数学模型为：

I(k_x,k_y)为风条纹能量谱；

(4)海面风向提取：根据傅里叶变换具有周期性，得到的风条纹能量谱I(k_x,k_y)关于一三或二四象限镜像对称，取对称的两个能量集中区域，连线的垂直方向则为风条纹平行方向，得到海面风向；将计算出来的方向与风向标测得的风向所在象限比较，保留与其象限一致的方向去除180°模糊的方向；

对计算的风向进行校正后，得到相对正北方向的海面风向，校正公式为：

N_w＝|θ_c|+|α|+|β|

其中，N_w为相对正北向的风向；θ_c为选取区域的中心角；β为笛卡尔坐标下计算得到的风向；α为船艏向。

3、风条纹波长尺度下限L_d为200m，风条纹波长尺度上限L_t为500m。

本发明的优势在于：

1、本发明提出在风向反演过程中增添一个相对辐射校正环节，消除了导航雷达回波强度径向非线性衰减对风条纹提取造成的影响，大大提高了风向反演精度和适用性；

2、本发明校正方法采用分段拟合校正，拟合的分段函数更符合导航雷达回波强度随径向距离非线性衰减的特征，保证校正后不会破坏海杂波特征；

3、本发明校正方法采用绝对值误差自适应判断分段点，实现了校正算法的自动操作，更有利于将算法应用到实际工程中；

4、本发明校正方法通过建立拉格朗日函数，应用多元函数极值求解出最小二乘回归系数，保证了在分段处函数的连续性；

5、本发明校正方法采用最高为二次多项式最小二乘拟合，提高了数据拟合速度，并简化了拟合方程的复杂度，增强了算法的适用性；

6、本发明基于波数能量谱尺度分离，应用相对辐射校正导航雷达图像反演海面风向。校正降低了径向回波强度的非线性对反演的影响，波数能量谱尺度分离反演风向不需要空间域缩减，能自动适应不同尺度的风条纹，提高了风向反演精度和计算速度，完全达到工程要求。

附图说明

图1为本发明的具体实施方式流程图；

图2a为导航雷达图像中值滤波前图像，图2b为导航雷达图像中值滤波后图像；

图3为去除遮挡区域导航雷达图像序列；

图4为径向灰度值分段拟合结果；

图5为未校正极坐标海面静态特征图像；

图6为校正后极坐标海面静态特征图像；

图7为最近点插值示意图；

图8为笛卡尔坐标海面静态特征图像；

图9为风条纹能量谱和能量谱带通滤波器；

图10为本发明风向与参考风向误差分布结果；

图11为传统风向与参考风向误差分布结果；

图12为实测导航雷达一条线上径向回波强度分布。

具体实施方式

下面结合附图举例对本发明做更详细地描述：

结合图1～12，本发明流程分为导航雷达序列采集、导航雷达图像预处理、导航雷达图像相对辐射校正和海面风向反演这四大块。

具体包括步骤如下：

(1)导航雷达图像序列采集：采集连续N幅导航雷达图像，将其定义为一组导航雷达图像序列，并收集对应时间和位置的实际风向与风速值；

(2)导航雷达图像预处理：对导航雷达图像序列极坐标图像进行归一化，使导航雷达图像序列中的每幅导航雷达图像的线数固定，再对单幅导航雷达图像通过观测去除遮挡区域的影响，只保留海杂波图像区域，应用图像中值滤波抑制同频信号对导航雷达图像的干扰；

(3)导航雷达图像相对辐射校正：对步骤(2)处理得到的导航雷达图像的回波强度值按方位向进行平均，得到关于距离的回波强度均值，对导航雷达图像径向回波强度进行相对辐射校正，得到相对辐射校正图像序列；

(4)海面风向反演：对相对辐射校正图像序列首先利用低通滤波器得到海面静态特征图像，然后基于波数能量谱尺度分离方法得到风条纹能量谱，最终获得海面风向。

步骤导航雷达图像相对辐射校正具体为：

步骤3.1，对导航雷达图像沿方位向进行平均，令导航雷达图像回波强度值为σ(θ_p,r_q)，其中θ_p(p＝1,2,…P)为方位向角度，共P个，r_q(q＝1,2,…Q)为径向距离，共Q个，生成一维径向平均回波强度值y：

$y\left(r_q\right)=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}\mathrm{\σ}({\mathrm{\θ}}_p,r_q)---\left(1\right)$

径向距离r_q定义为x，则x和y形成的数据集为数据点总数为n。

步骤3.2，首先选取x中间位置x₀，将n分成两段数据，个数分别为n₁和n₂，若K为间断点数，则K＝1，生成两个数据集进行分段拟合。

步骤3.3，假设数据共有K个间断点，则形成K+1个数据集，即 ${\mathrm{\Ω}}_k={\left\{(x_i^k,y_i^k)\right\}}_{i=1}^{n_k},k=\mathrm{1,2},...K+1,$ n_k为第k个数据集的点数， $x_{k-1}{\≤x}_i^k\≤x_k,$ x_k是数据集Ω_k和Ω_k+1的分段点。

步骤3.4，拟合函数f(x)的形式如下：

$f\left(x\right)=\left\{\begin{array}{c}{f}_1\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^1{h}_j\left(x\right),1\≤x\≤x_1\\ f_2\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^2{h}_j\left(x\right),x_1\≤x\≤x_2\\ \·\\ \·\\ \·\\ f_k\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^2{h}_j\left(x\right),x_{k-1}\≤x\≤x_k\\ \·\\ \·\\ \·\end{array}\right.---\left(2\right)$

式中f_k-1(x_k)＝f_k(x_k)，是数据集Ω_k的基函数，基函数一般选用简单的函数形式方便计算，本文选为1，x，x2；m为基函数的个数，即m＝3；为待确定回归系数。为了使拟合误差达到最小且在x_k处连续，则建立最小二乘回归模型如下：

$\begin{array}{c}\underset{{\mathrm{\α}}_j^k}{\min }\underset{i=1}{\overset{n_1}{\mathrm{\Σ}}}{(f_1\left(x_i^1\right)-y_i^1)}^2+\underset{i=1}{\overset{n_2}{\mathrm{\Σ}}}{(f_2\left(x_i^2\right)-y_i^2)}^2+...+\underset{i=1}{\overset{n_k}{\mathrm{\Σ}}}{(f_k\left(x_i^k\right)-y_i^k)}^2+...\underset{i=1}{\overset{n_k}{\mathrm{\Σ}}}{(f_{K+1}\left(x_i^{K+1}\right)-y_i^{K+1})}^2,\\ s.t.f_{k-1}\left(x_k\right)=f_k\left(x_k\right)\end{array}---3$

令基函数h(x)＝[1 x x²]，拟合参数行向量分别为： $X_k=\left[\begin{array}{c}h\left(x_1^k\right)\\ h\left(x_2^k\right)\\ \·\\ \·\\ \·\\ \text{h}\left(x_{n_k}^k\right)\end{array}\right],$ $y_k=\left[\begin{array}{c}{y}_1^k\\ y_2^k\\ \·\\ \·\\ \·\\ y_{n_k}^k\end{array}\right],$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_1^k\\ {\mathrm{\α}}_2^k\\ {\mathrm{\α}}_3^k\end{array}\right],$ k＝1,2,…K+1；则多段数据集拟合参数可以表示为：

$y=\left[\begin{array}{c}{y}_1^T\\ y_2^T\\ \·\\ \·\\ \·\\ y_{K+1}^T\end{array}\right]$ $\mathrm{\α}=\left[\begin{array}{c}{\mathrm{\α}}_1^T\\ {\mathrm{\α}}_2^T\\ \·\\ \·\\ \·\\ {\mathrm{\α}}_{K+1}^T\end{array}\right]$ $X=\left[\begin{array}{cccc}{X}_1& 0& ...& 0\\ 0& X_2& ...& 0\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ 0& 0& ...& X_{K+1}\end{array}\right]$

$Z=\left[\begin{array}{cccccc}h\left(x_1\right)& -h\left(x_1\right)& 0& ...& 0& 0\\ 0& h\left(x_2\right)& -h\left(x_2\right)& ...& 0& 0\\ \·& \·& \·& & \·& \·\\ \·& \·& \·& & \·& \·\\ \·& \·& \·& & \·& \·\\ 0& 0& 0& ...& -h\left(x_{K-1}\right)& 0\\ 0& 0& 0& ...& h\left(x_K\right)& -h\left(x_K\right)\end{array}\right]---\left(4\right)$

则公式(3)的二次范数形式为：

$\underset{\mathrm{\α}}{\min }{\left|\right|\mathrm{X\α}-y\left|\right|}^2,s.t.\mathrm{Z\α}=0---\left(5\right)$

对(5)式建立拉格朗日函数为：

$L(\∂,\mathrm{\λ})={\left|\right|\mathrm{X\α}-Y\left|\right|}^2+2{\mathrm{\λ}}^T\mathrm{Z\α}---\left(6\right)$

其中λ为列向量约束参数，即λ＝[λ₁, λ₂, … ,λ_K]^T。由多元函数极值必要条件，得到:

$\left\{\begin{array}{c}\frac{\∂L}{{\∂\mathrm{\α}}_1}=2X^T(\mathrm{X\α}-y)+2Z^T\mathrm{\λ}=0\\ \frac{\∂L}{\∂\mathrm{\λ}}=\mathrm{Z\α}=0\end{array}\right.---\left(7\right)$

由上式可得：

$\mathrm{\α}=\widehat{\mathrm{\α}}-{\left(X^{T}X\right)}^{-1}{Z}^T\mathrm{\λ}---\left(8\right)$

其中将公式(8)带入公式(7)的第二个等式可解得：

$\mathrm{\λ}={\left[Z{\left(X^{T}X\right)}^{-1}{Z}^T\right]}^{-1}Z\widehat{\mathrm{\α}}---\left(9\right)$

将(9)式带入(8)式可获得最小二乘回归系数α的解，根据导航雷达回波机制数据集Ω_k可以保证公式(8)中的是非退化的，Z是非零的，因此，α的解是唯一的，从而获得有多段拟合函数f(x)。

步骤3.5，计算两个数据集绝对值误差均值S_k,k＝1,2…K+1，每点的绝对值误差e_i,i＝1,2,…n，若e_i>S_k&e_i+1>S_k&e_i+2>S_k，则在e_i的处位置再次分段，执行步骤3.6；其他情况或分段内点数少于十个，则不分段，执行步骤3.7。

步骤3.6，在分段点处对数据再次进行分段，返回从步骤3.3开始重新拟合，K为分段点的数量，再执行步骤3.5的判断。最终根据以上步骤得到数据集最优分段函数拟合函数值。

步骤3.7，对导航雷达图像中每个θ_p的回波强度值σ(θ_p,r_q)按r_q分段的减去拟合函数值，实现导航雷达图像相对辐射校正。校正遍历导航雷达图像序列的每幅图像，最终得到相对辐射校正后图像序列f′(θ,r,t)，θ为方位向角度，r为径向距离，t为图像序列时间轴。

本发明应用波数能量谱尺度分离方法来提取海面风向信息，步骤海面风向反演具体为：

步骤4.1，为获得包含风条纹特征的海面静态特征图像，需要构建的全局低通滤波器，应用下式实现：

$f(\mathrm{\θ},r)=\frac{\underset{t=1}{\overset{N_t}{\mathrm{\Σ}}}{f}^{\′}(\mathrm{\θ},r,t)}{N_t}---\left(10\right)$

式中f(θ,r)为海面静态特征图像，在笛卡尔坐标下可以表示为 和分别为笛卡尔坐标下图像沿x和y轴的坐标，N_x和N_y为坐标位置，N_t为时间序列。

步骤4.2，本发明是基于图像能量谱特征的，将海面静态特征图像变换到能量谱域，这里应用二维离散快速傅里叶变换(2DFFT)，即：

$F(k_x,k_y)=\underset{N_x}{\mathrm{\Σ}}\underset{N_y}{\mathrm{\Σ}}f(x_{N_x},y_{N_y})\exp [-i*2\mathrm{\π}(\frac{z_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}{*k}_y}{\max \left(y_{N_y}\right)})]---\left(11\right)$

其中，F(k_x,k_y)为f(x,y)的傅里叶系数。复指数项可以写为：

$\begin{array}{c}\exp [-i*2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})]=\cos 2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})\\ +i\sin 2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})\end{array}---\left(12\right)$

其中，

$k_x=\frac{2\mathrm{\π}}{\max \left(x_{N_x}\right)*d_x}---\left(13\right)$

$k_{\text{y}}=\frac{2\mathrm{\π}}{\max \left(y_{N_y}\right)*d_y}---\left(14\right)$

(k_x,k_y)为在能量谱Τ下的坐标，d_x、d_y为导航雷达图像分辨率。根据二维谱性质得到的海面静态特征图像能量谱为A(k_x,k_y)：

$A(k_x,k_y)=\sqrt{\left[\mathrm{Re}\left(F(k_x,k_y)\right)\right]+\left[\mathrm{Im}\left(F(k_x,k_y)\right)\right]}---\left(15\right)$

步骤4.3，基于尺度分离获得风条纹能量谱。导航雷达图像呈现的海面条纹波长L与能量谱波数k的关系式为：

$k=\frac{2\mathrm{\π}}{L}---\left(16\right)$

其中，(k_x,k_y)为能量谱域Τ的坐标，即k_x和k_y为k在x和y轴的波数分量的。

设L_d为风条纹波长尺度下限，L_t为风条纹波长尺度上限，根据公式(16)得到风条纹能量谱波数下限为k_d：

$k_d=\frac{2\mathrm{\π}}{L_d}---\left(17\right)$

风条纹能量谱波数上限为k_t：

$k_t=\frac{2\mathrm{\π}}{L_t}---\left(18\right)$

根据风条纹能量谱波数上下限设计能量谱带通滤波器，将风条纹信号能量谱提取出来，数学模型为：

I(k_x,k_y)为风条纹能量谱。

步骤4.4，海面风向提取。根据傅里叶变换具有周期性，得到的风条纹能量谱I(k_x,k_y)关于一三或二四象限镜像对称。取对称的两个能量集中区域，连线的垂直方向则为风条纹平行方向，由于风条纹方向与风向平行，所以可得到海面风向。计算出的风向存在180°模糊问题，为了解决这个问题将计算出来的方向与风向标测得的风向所在象限比较，保留与其象限一致的方向去除180°模糊的方向，从而获得准确的海面风向。

由于选取区域和船艏向的影响，计算的风向需要校正后才能得到相对正北方向的海面风向，校正公式为：

N_w＝|θ_c|+|α|+|β|       (20)

其中，N_w为相对正北向的风向；θ_c为选取区域的中心角；β为笛卡尔坐标下计算得到的风向；α为船艏向。

下面给出具体实例：

具体实施步骤共分为十五步，第一步为导航图像序列采集；第二步到第四步是导航雷达图像预处理；第五步到底十步为导航雷达图像相对辐射校正；第十一步到第十五步为基于波数能量谱域尺度分离提取出海面风向信息。具体步骤如下：

第一步，采集32幅位置相同、时间连续的导航雷达图像，定义为一个导航图像序列，记录其时间总长度为T(约1.5分钟)，同步记录同时空由风向标获得的风向θ_w、风速计测得的风速U_w。选取2010年10月22日，10：35的导航雷达图像序列为例，此时风向为36°，风速为17.1M/s。

第二步，对导航雷达图像应用3×3模板的2-D非线性平滑中值滤波，抑制同频信号对海面风向的影响。

$g^{\′}(r,\mathrm{\θ})\underset{(s,t)\∈N(r,\mathrm{\θ})}{\mathrm{median}}\left\{g(s,t)\right\}---\left(21\right)$

式中g(s,t)为极坐标(s,t)处的回波强度值；g'(r,θ)为滤波后在极坐标(r,θ)处的回波强度值；N(r,θ)为(r,θ)为中心的像元点，(s,t)取以(r,θ)为中心相邻的8个像元点。

将3×3模板中值滤波器的模板中心与极坐标图像的N(r,θ)中心重合，将(r,θ)与周围8个相邻像元点(s,t)的回波强度值进行比较，取中间回波强度值更新N(r,θ)中心点值，模板遍历所选导航雷达图像序列，最终获得中值滤波后的导航雷达图像序列，图2为本发明中值滤波前后导航雷达图像。

第三步，由于导航雷达工作时脉冲数量不定和天线旋转时受到各种外界环境的干扰，导致导航雷达图像方位向线数不固定，本发明提出方位向归一化来解决这个问题，具体步骤如下：

①读取滤波后导航雷达图像方位向线数和径向点数，方位向大约3600条射线状线，即大约N＝3600个角度值，Ω_i,i＝1,2,…N，径向600个点，方位向间隔大约是0.1°。

②建立极坐标新图像，固定1800条线，即N_new＝1800个角度值，θ_j,j＝1,2,…N_new，径向600个点，方位向固定间隔0.2°；

③若Ω_i＝θ_j，或第一个Ω_i>θ_j，则将Ω_i所对应线的值赋到新图像θ_j所对应的线上；

④不断重复③直到新图像上的N_new条线上都具有读取滤波后导航雷达图像的回波强度值g'(θ,r)，新导航雷达图像灰度值为σ(θ,r)，从而获得方位向归一化的新极坐标图像序列。

第四步，对新极坐标下的导航雷达图像序列进行观测，去除岸基、遮挡对导航雷达图像的影响。图3为去除遮挡后的导航雷达图像序列，方位向选取106°～291°、径向选取600m～2250m,图像中白色区域为去除的区域。

第五步，对新极坐标下的导航雷达图像沿方位向进行平均，令导航雷达回波强度值为σ(θ_p,r_q)，其中θ_p(p＝1,2,…P)为方位向角度，P为选取区域线数926条，r_q(q＝1,2,…Q)为径向距离，Q为径向点数238个，生成一维径向平均回波强度值y：

$y\left(r_q\right)=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}\mathrm{\σ}({\mathrm{\θ}}_p,r_q)---\left(22\right)$

径向距离r_q定义为x，则x和y形成的数据集为数据点总数为n＝238。

第六步，首先选取x中间位置x₀，将n分成两段数据，个数分别为n₁＝119和n₂＝120，若K为间断点数，则K＝1，生成两个数据集进行分段拟合。

假设数据共有K个间断点，则形成K+1个数据集，即 ${\mathrm{\Ω}}_k={\left\{(x_i^k,y_i^k)\right\}}_{i=1}^{n_k},k=\mathrm{1,2},...K+1,$ n_k为第k个数据集的点数， $x_{k-1}{\≤x}_i^k\≤x_k,$ x_k是数据集Ω_k和Ω_k+1的分段点。

第七步，拟合函数f(x)的形式如下：

$f\left(x\right)=\left\{\begin{array}{c}{f}_1\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^1{h}_j\left(x\right),1\≤x\≤x_1\\ f_2\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^2{h}_j\left(x\right),x_1\≤x\≤x_2\\ \·\\ \·\\ \·\\ f_k\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^2{h}_j\left(x\right),x_{k-1}\≤x\≤x_k\\ \·\\ \·\\ \·\end{array}\right.---\left(23\right)$

式中f_k-1(x_k)＝f_k(x_k)，是数据集Ω_k的基函数，基函数一般选用简单的函数形式方便计算，本文选为1，x，x²；m为基函数的个数，即m＝3；为待确定回归系数。为了使拟合误差达到最小且在x_k处连续，则建立最小二乘回归模型如下：

$\begin{array}{c}\underset{{\mathrm{\α}}_j^k}{\min }\underset{i=1}{\overset{n_1}{\mathrm{\Σ}}}{(f_1\left(x_i^1\right)-y_i^1)}^2+\underset{i=1}{\overset{n_2}{\mathrm{\Σ}}}{(f_2\left(x_i^2\right)-y_i^2)}^2+...+\underset{i=1}{\overset{n_k}{\mathrm{\Σ}}}{(f_k\left(x_i^k\right)-y_i^k)}^2+...\underset{i=1}{\overset{n_k}{\mathrm{\Σ}}}{(f_{K+1}\left(x_i^{K+1}\right)-y_i^{K+1})}^2,\\ s.t.f_{k-1}\left(x_k\right)=f_k\left(x_k\right)\end{array}---\left(24\right)$

令基函数h(x)＝[1 x x²]，拟合参数行向量分别为： $X_k=\left[\begin{array}{c}h\left(x_1^k\right)\\ h\left(x_2^k\right)\\ \·\\ \·\\ \·\\ \text{h}\left(x_{n_k}^k\right)\end{array}\right],$ $y_k=\left[\begin{array}{c}{y}_1^k\\ y_2^k\\ \·\\ \·\\ \·\\ y_{n_k}^k\end{array}\right],$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_1^k\\ {\mathrm{\α}}_2^k\\ {\mathrm{\α}}_3^k\end{array}\right],$ k＝1,2,…K+1；则多段数据集拟合参数可以表示为：

$y=\left[\begin{array}{c}{y}_1^T\\ y_2^T\\ \·\\ \·\\ \·\\ y_{K+1}^T\end{array}\right]\mathrm{\α}=\left[\begin{array}{c}{\mathrm{\α}}_1^T\\ {\mathrm{\α}}_2^T\\ \·\\ \·\\ \·\\ {\mathrm{\α}}_{K+1}^T\end{array}\right]X=\left[\begin{array}{cccc}{X}_1& 0& ...& 0\\ 0& X_2& ...& 0\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ 0& 0& ...& X_{K+1}\end{array}\right]$

$Z=\left[\begin{array}{cccccc}h\left(x_1\right)& -h\left(x_1\right)& 0& ...& 0& 0\\ 0& h\left(x_2\right)& -h\left(x_2\right)& ...& 0& 0\\ \·& \·& \·& & \·& \·\\ \·& \·& \·& & \·& \·\\ \·& \·& \·& & \·& \·\\ 0& 0& 0& ...& -h\left(x_{K-1}\right)& 0\\ 0& 0& 0& ...& h\left(x_K\right)& -h\left(x_K\right)\end{array}\right]---\left(25\right)$

则公式(4)的二次范数形式为：

$\underset{\mathrm{\α}}{\min }{\left|\right|\mathrm{X\α}-y\left|\right|}^2,s.t.\mathrm{Z\α}=0---\left(26\right)$

对(6)式建立拉格朗日函数为：

$L(\∂,\mathrm{\λ})={\left|\right|\mathrm{X\α}-Y\left|\right|}^2+2{\mathrm{\λ}}^T\mathrm{Z\α}---\left(27\right)$

其中λ为列向量约束参数，即λ＝[λ₁, λ₂, … ,λ_K]^T。由多元函数极值必要条件，得到:

$\left\{\begin{array}{c}\frac{\∂L}{{\∂\mathrm{\α}}_1}=2X^T(\mathrm{X\α}-y)+2Z^T\mathrm{\λ}=0\\ \frac{\∂L}{\∂\mathrm{\λ}}=\mathrm{Z\α}=0\end{array}\right.---\left(28\right)$

由上式可得：

$\mathrm{\α}=\widehat{\mathrm{\α}}-{\left(X^{T}X\right)}^{-1}{Z}^T\mathrm{\λ}---\left(29\right)$

其中将公式(29)带入公式(28)的第二个等式可解得：

$\mathrm{\λ}={\left[Z{\left(X^{T}X\right)}^{-1}{Z}^T\right]}^{-1}Z\widehat{\mathrm{\α}}---\left(30\right)$

将(30)式带入(29)式可获得最小二乘回归系数α的解，根据导航雷达回波机制数据集Ω_k可以保证公式(8)中的是非退化的，Z是非零的，因此，α的解是唯一的，从而获得有多段拟合函数f(x)。

第八步，计算两个数据集的绝对值误差均值S_k,k＝1,2…K+1，及每点的绝对值误差e_i,i＝1,2,…n，若e_i>S_k&e_i+1>S_k&e_i+2>S_k，则在e_i的处位置再次分段，执行第九步；其他情况或分段内点数少于十个，则不分段，执行第十步。

第九步，在分段点处对数据再次进行分段，返回从第六步开始重新拟合，K为分段点的数量，再执行第八步的判断。最终根据以上步骤得到数据最优分段函数。本发明提出使用新型相对辐射校正方法来校正导航雷达径向回波强度，应用的是自适应拉格朗日最小二乘分段拟合，结果见图4。与平均径向灰度值相比偏差0.83，标准误差22.19，相关性0.9986，拟合接近实际分布。

第十步，对新极坐标下的导航雷达图像每个θ_i的回波强度值σ(θ_i,r_j)按r_j分段的减去拟合函数值，实现导航雷达图像相对辐射校正。校正遍历导航雷达图像序列的每幅导航图像，最终得到相对辐射校正后导航雷达图像序列f′(θ,r,t)，θ为方位向角度，r为径向距离，t为图像序列时间轴。

第十一步，为获得包含风条纹特征的海面静态特征图像，需要构建的全局低通滤波器，应用下式实现：

$f(\mathrm{\θ},r)=\frac{\underset{t=1}{\overset{N_t}{\mathrm{\Σ}}}{f}^{\′}(\mathrm{\θ},r,t)}{N_t}---\left(31\right)$

式中f'(θ,r,t)为相对辐射校正后图像序列，f(θ,r)为海面静态特征极坐标图像，N_t为图像序列中包含雷达图像的个数，N_t＝32。

图5为未校正极坐标海面静态特征图像，图6为校正后海面静态特征图像。从图中可以看出，相对辐射校正不但没有改变风条纹尺度分布特征，反而降低了导航雷达图像径向衰减对风条纹成像的影响，使得风条纹相对背景信号更突出的表现出来。

第十二步，为了进行图像频域转换，需要将极坐标下的海面静态特征图像插值成笛卡尔坐标下的图像，具体步骤如下：

①在极坐标海面静态特征图像中选取相对北向203°方向、距离导航雷达平台距离630m、空间尺寸为960m*960m的近似矩形区域，由于图像像元分辨率为7.5m，即选取128*128个像元点，选取区域见图6中黑色框区域。

②近似矩形区域海面静态特征图像像元点极坐标的形式表示为(r,θ)，转换为笛卡儿坐标后表示为(x,y)，两者存在下式关系：

$\left\{\begin{array}{c}x=r*\cos \mathrm{\θ}\\ y=r*\sin \mathrm{\θ}\end{array}\right.$

③设建立的矩形区域的笛卡尔坐标为(x_i,y_i)，利用最近点插值就是(x_i,y_i)找到距离最近的(x,y)，并将(x,y)对应的极坐标(r_i,θ_i)像元灰度值赋值给建立的矩形区域的(x_i,y_i)点。图7为最近点插值示意图。

④若任取建立的矩形笛卡尔坐标记为(x₀,y₀)，则离其最近点的图像极坐标(r₀,θ₀)为：

$\left\{\begin{array}{c}{r}_0=\mathrm{round}\left(\mathrm{sqrt}({x_0}^2+{y_0}^2)\right)\\ {\mathrm{\θ}}_0=\mathrm{round}\left(\mathrm{rem}(\arctan (y_0,x_0)+2\mathrm{\π},2\mathrm{\π})\right)\end{array}\right.$

其中，rem()为求余函数，round()为向最近点取整函数。图8为插值后得到的笛卡尔坐标下海面静态特征图像和分别为笛卡尔坐标下图像沿x和y轴的坐标，N_x和N_y为坐标位置，

第十三步，本发明是基于图像能量谱特征的，将笛卡尔坐标下的海面静态特征图像变换到能量谱域，这里应用二维离散快速傅里叶变换(2DFFT)，即：

$F(k_x,k_y)=\underset{N_x}{\mathrm{\Σ}}\underset{N_y}{\mathrm{\Σ}}f(x_{N_x},y_{N_y})\exp [-i*2\mathrm{\π}(\frac{z_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}{*k}_y}{\max \left(y_{N_y}\right)})]---\left(34\right)$

其中，F(k_x,k_y)为f(x,y)的傅里叶系数。复指数项可以写为：

$\begin{array}{c}\exp [-i*2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})]=\cos 2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})\\ +i\sin 2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})\end{array}---\left(35\right)$

其中，

$k_x=\frac{2\mathrm{\π}}{\max \left(x_{N_x}\right)*d_x}---\left(36\right)$

$k_y=\frac{2\mathrm{\π}}{\max \left(y_{N_y}\right)*d_y}---\left(37\right)$

其中，(k_x,k_y)为在能量谱Τ下的坐标，d_x＝7.5m、d_y＝7.5m为导航雷达图像分辨率。

根据二维谱性质得到的笛卡尔坐标下海面静态特征图像能量谱为A(k_x,k_y)：

$A(k_x,k_y)=\sqrt{\left[\mathrm{Re}\left(F(k_x,k_y)\right)\right]+\left[\mathrm{Im}\left(F(k_x,k_y)\right)\right]}---\left(38\right)$

第十四步，基于尺度分离获得风条纹能量谱。导航雷达图像呈现的海面条纹波长L与能量谱波数k的关系式为：

$k=\frac{2\mathrm{\π}}{L}---\left(39\right)$

其中，(k_x,k_y)为能量谱域Τ的坐标，即k_x和k_y为k在x和y轴的波数分量的，由第十三步可以获得。

根据风条纹在导航雷达图像上的成像特征，可知L_d＝200m为风条纹波长尺度下限，L_t＝500m为风条纹波长尺度上限，根据公式(39)得到风条纹能量谱波数下限为k_d：

$k_d=\frac{2\mathrm{\π}}{L_d}---\left(40\right)$

风条纹能量谱波数上限为k_t：

$k_t=\frac{2\mathrm{\π}}{L_t}---\left(41\right)$

根据风条纹能量谱波数上下限设计能量谱带通滤波器，将风条纹信号能量谱提取出来，数学模型为：

I(k_x,k_y)为风条纹能量谱，图9为滤波后的风条纹能量谱和能量谱带通滤波器的尺度带。

第十五步，海面风向提取。根据傅里叶变换具有周期性，得到的风条纹能量谱I(k_x,k_y)关于一三或二四象限镜像对称。取对称的两个能量集中区域，连线的垂直方向则为风条纹平行方向，由于风条纹方向与风向平行，所以可得到海面风向。

由于选取区域和船艏向的影响，计算的风向需要校正后才能得到相对正北方向的海面风向，校正公式为：

N_w＝|θ_c|+|α|+|β|       (43)

其中，N_w＝203°为相对正北向的风向；θ_c＝110°为选取区域的中心角；β为笛卡尔坐标下计算得到的风向，通过图9计算β＝0°，得到；α＝93°为船艏向；

计算出的风向存在180°模糊问题，为了解决这个问题将计算出来的方向与风向标测得的风向所在象限比较，保留与其象限一致的方向去除180°模糊的方向，从而获得准确的海面风向。N_w去除180°模糊得到海面风向为23°，而此时实际风向为36°，相差13°符合工程要求。下面通过大量数据序列来验证算法的有效性和适用性。

本发明公开的是一种基于导航雷达图像径向拟合校正的海面风向反演方法开展实验，导航雷达平台及风速计和风向标安装在海坛岛岸基，海坛岛位于福建平潭县，此海域平均海深25m，受地形影响此海域经常出现高海况(>4级)。导航雷达具体参数见表1。实验数据从2010年10月到11月，共采集1634组数据。由于台风“鲶鱼”的影响，此段时间风向主要是东北风，短暂出现西南风。参考数据为同位置的风速计和风向标采集的每分钟海面风向和风速值，为了与雷达图像时间段对比，对两分钟内的海面风向和风速值进行平均，得到‘参考风向’和‘参考风速’。

表1  导航雷达具体参数

为验证本发明的有效性，将由导航雷达图像径向拟合校正后反演的海面风向以下称“本发明风向”，由局部梯度法反演的海面风向以下称“传统风向”，分别统计本发明风向与传统风向和参考风向的误差分布，见图10、11。通过对两个图比较可以发现本发明风向相对参考风向的误差更小，说明本发明风向比传统风向更接近参考风向。传统风向与参考风向有存在20°左右的误差，主要是由于传统风向反演过程中在进行空间缩减时，将图像分辨率统一缩减到30m，没有针对不同风条纹的尺度进行判断。而本发明风向与参考风向的误差都在10°左右，本发明方法明显优于传统方法。具体两种方法海面风向误差统计如表2。

表2  海面风向误差统计

本发明提供了一种基于导航雷达图像径向校正的海面风向反演方法。首先，利用自适应拉格朗日最小二乘分段拟合，拟合出导航雷达图像径向平均值符合的分段函数；然后，将导航雷达图像依据分段函数进行相对辐射校正，从而去除导航雷达成像机制导致的径向非线性衰减；最后，利用波数能量谱尺度分离法将海面风向提取出来。较传统局部梯度法反演风向，不仅风向反演精度有所提高，而且避免了对海面静态特征图像进行空间缩减时带来的误差。本发明反演的海面风向与参考风向的相关系数达到了0.9924，标准差为8.05°，偏差为-0.85°，完全符合工程要求。

Claims (4)

1.一种基于相对辐射校正的导航雷达图像反演海面风向方法，其特征是：

(1)导航雷达图像序列采集：采集连续N幅导航雷达图像，将其定义为一组导航雷达图像序列，并收集对应时间和位置的实际风向与风速值；

(2)导航雷达图像预处理：对导航雷达图像序列极坐标图像进行归一化，使导航雷达图像序列中的每幅导航雷达图像的线数固定，再对单幅导航雷达图像通过观测去除遮挡区域的影响，只保留海杂波图像区域，应用图像中值滤波抑制同频信号对导航雷达图像的干扰；

(3)导航雷达图像相对辐射校正：对步骤(2)处理得到的导航雷达图像的回波强度值按方位向进行平均，得到关于距离的回波强度均值，对导航雷达图像径向回波强度进行相对辐射校正，得到相对辐射校正图像序列；

(4)海面风向反演：对相对辐射校正图像序列首先利用低通滤波器得到海面静态特征图像，然后基于波数能量谱尺度分离方法得到风条纹能量谱，最终获得海面风向。

2.根据权利要求1所述的一种基于相对辐射校正的导航雷达图像反演海面风向方法，其特征是：导航雷达图像相对辐射校正的具体步骤为：

(1)对导航雷达图像沿方位向进行平均，令导航雷达图像回波强度值为σ(θ_p,r_q)，其中θ_p(p＝1,2,…P)为方位向角度，共P个，r_q(q＝1,2,…Q)为径向距离，共Q个，生成一维径向平均回波强度值y：

$y\left(r_q\right)=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}\mathrm{\σ}({\mathrm{\θ}}_p,r_q)$

径向距离r_q定义为x，则x和y形成的数据集为数据点总数为n；

(2)选取x中间位置x₀，将n分成两段数据，个数分别为n₁和n₂，若K为间断点数，则K＝1，生成两个数据集进行分段拟合；

(3)假设数据共有K个间断点，则形成K+1个数据集，即 ${\mathrm{\Ω}}_k={\left\{(x_i^k,y_i^k)\right\}}_{i=1}^{n_k},k=\mathrm{1,2},...K+1,$ n_k为第k个数据集的点数， $x_{k-1}{\≤x}_i^k\≤x_k,$ x_k是数据集Ω_k和Ω_k+1的分段点；

(4)拟合函数f(x)的形式如下：

$f\left(x\right)=\left\{\begin{array}{c}{f}_1\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^1{h}_j\left(x\right),1\≤x\≤x_1\\ f_2\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^2{h}_j\left(x\right),x_1\≤x\≤x_2\\ \·\\ \·\\ \·\\ f_k\left(x\right)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j^2{h}_j\left(x\right),x_{k-1}\≤x\≤x_k\\ \·\\ \·\\ \·\end{array}\right.$

式中f_k-1(x_k)＝f_k(x_k)，是数据集Ω_k的基函数，基函数选为1，x，x2；m为基函数的个数，即m＝3；为待确定回归系数；为了使拟合误差达到最小且在x_k处连续，建立最小二乘回归模型如下： $\begin{array}{c}\underset{{\mathrm{\α}}_j^k}{\min }\underset{i=1}{\overset{n_1}{\mathrm{\Σ}}}{(f_1\left(x_i^1\right)-y_i^1)}^2+\underset{i=1}{\overset{n_2}{\mathrm{\Σ}}}{(f_2\left(x_i^2\right)-y_i^2)}^2+...+\underset{i=1}{\overset{n_k}{\mathrm{\Σ}}}{(f_k\left(x_i^k\right)-y_i^k)}^2+...\underset{i=1}{\overset{n_k}{\mathrm{\Σ}}}{(f_{K+1}\left(x_i^{K+1}\right)-y_i^{K+1})}^2,\\ s.t.f_{k-1}\left(x_k\right)=f_k\left(x_k\right)\end{array}$

令基函数h(x)＝[1 x x²]，拟合参数行向量分别为： $X_k=\left[\begin{array}{c}h\left(x_1^k\right)\\ h\left(x_2^k\right)\\ \·\\ \·\\ \·\\ \text{h}\left(x_{n_k}^k\right)\end{array}\right],$ $y_k=\left[\begin{array}{c}{y}_1^k\\ y_2^k\\ \·\\ \·\\ \·\\ y_{n_k}^k\end{array}\right],$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_1^k\\ {\mathrm{\α}}_2^k\\ {\mathrm{\α}}_3^k\end{array}\right],$ k＝1,2,…K+1；则多段数据集拟合参数可以表示为：

$y=\left[\begin{array}{c}{y}_1^T\\ y_2^T\\ \·\\ \·\\ \·\\ y_{K+1}^T\end{array}\right]$ $\mathrm{\α}=\left[\begin{array}{c}{\mathrm{\α}}_1^T\\ {\mathrm{\α}}_2^T\\ \·\\ \·\\ \·\\ {\mathrm{\α}}_{K+1}^T\end{array}\right]$ $X=\left[\begin{array}{cccc}{X}_1& 0& ...& 0\\ 0& X_2& ...& 0\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ 0& 0& ...& X_{K+1}\end{array}\right]$

$Z=\left[\begin{array}{cccccc}h\left(x_1\right)& -h\left(x_1\right)& 0& ...& 0& 0\\ 0& h\left(x_2\right)& -h\left(x_2\right)& ...& 0& 0\\ \·& \·& \·& & \·& \·\\ \·& \·& \·& & \·& \·\\ \·& \·& \·& & \·& \·\\ 0& 0& 0& ...& -h\left(x_{K-1}\right)& 0\\ 0& 0& 0& ...& h\left(x_K\right)& -h\left(x_K\right)\end{array}\right],$

则最小二乘回归模型的二次范数形式为：

$\underset{\mathrm{\α}}{\min }{\left|\right|\mathrm{X\α}-y\left|\right|}^2,s.t.\mathrm{Z\α}=0,$

对上式建立拉格朗日函数为：

$L(\∂,\mathrm{\λ})={\left|\right|\mathrm{X\α}-Y\left|\right|}^2+2{\mathrm{\λ}}^T\mathrm{Z\α},$

其中λ为列向量约束参数，即λ＝[λ₁， λ₂， … ,λ_K]^T，由多元函数极值必要条件，得到:

$\left\{\begin{array}{c}\frac{\∂L}{{\∂\mathrm{\α}}_1}=2X^T(\mathrm{X\α}-y)+2Z^T\mathrm{\λ}=0\\ \frac{\∂L}{\∂\mathrm{\λ}}=\mathrm{Z\α}=0\end{array}\right.,$

由上式可得：

$\mathrm{\α}=\widehat{\mathrm{\α}}-{\left(X^{T}X\right)}^{-1}{Z}^T\mathrm{\λ},$

其中 $\widehat{\mathrm{\α}}={\left(X^{T}X\right)}^{-1}{X}^{T}y,$ 将 $\mathrm{\α}=\widehat{\mathrm{\α}}-{\left(X^{T}X\right)}^{-1}{Z}^T\mathrm{\λ}$ 带入 $\left\{\begin{array}{c}\frac{\∂L}{{\∂\mathrm{\α}}_1}=2X^T(\mathrm{X\α}-y)+2Z^T\mathrm{\λ}=0\\ \frac{\∂L}{\∂\mathrm{\λ}}=\mathrm{Z\α}=0\end{array}\right.,$ 可解得：

$\mathrm{\λ}={\left[Z{\left(X^{T}X\right)}^{-1}{Z}^T\right]}^{-1}Z\widehat{\mathrm{\α}},$

将式带入获得最小二乘回归系数α的解，根据导航雷达回波机制数据集Ω_k可以保证式中的是非退化的，Z是非零的，因此，α的解是唯一的，从而获得有多段拟合函数f(x)；

(5)计算两个数据集绝对值误差均值S_k,k＝1,2…K+1，每点的绝对值误差e_i,i＝1,2,…n，若e_i>S_k&e_i+1>S_k&e_i+2>S_k，则在e_i的处位置再次分段，执行步骤(6)；其他情况或分段内点数少于十个，则不分段，执行步骤(7)；

(6)在分段点处对数据再次进行分段，返回从步骤(3)开始重新拟合，K为分段点的数量，再执行步骤(5)的判断；得到数据集最优分段函数拟合函数值；

(7)对导航雷达图像中每个θ_p的回波强度值σ(θ_p,r_q)按r_q分段的减去拟合函数值，实现导航雷达图像相对辐射校正；校正遍历导航雷达图像序列的每幅图像，最终得到相对辐射校正后图像序列f′(θ,r,t)，θ为方位向角度，r为径向距离，t为图像序列时间轴。

3.根据权利要求1或2所述的一种基于相对辐射校正的导航雷达图像反演海面风向方法，其特征是：海面风向反演的具体步骤为：

(1)构建全局低通滤波器以获得包含风条纹特征的海面静态特征图像：

$f(\mathrm{\θ},r)=\frac{\underset{t=1}{\overset{N_t}{\mathrm{\Σ}}}{f}^{\′}(\mathrm{\θ},r,t)}{N_t}$

f(θ,r)为海面静态特征图像，在笛卡尔坐标下可以表示为 和分别为笛卡尔坐标下图像沿x和y轴的坐标，N_x和N_y为坐标位置，N_t为时间序列；

(2)应用二维离散快速傅里叶变换将海面静态特征图像变换到能量谱域，即：

$F(k_x,k_y)=\underset{N_x}{\mathrm{\Σ}}\underset{N_y}{\mathrm{\Σ}}f(x_{N_x},y_{N_y})\exp [-i*2\mathrm{\π}(\frac{z_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}{*k}_y}{\max \left(y_{N_y}\right)})]$ F(k_x,k_y)为f(x,y)的傅里叶系数，复指数项为：

$\begin{array}{c}\exp [-i*2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})]=\cos 2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})\\ +i\sin 2\mathrm{\π}(\frac{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})\end{array}$

其中，

$k_x=\frac{2\mathrm{\π}}{\max \left(x_{N_x}\right)*d_x}$

$k_y=\frac{2\mathrm{\π}}{\max \left(y_{N_y}\right)*d_y}$

(k_x,k_y)为在能量谱Τ下的坐标，d_x、d_y为导航雷达图像分辨率；根据二维谱性质得到的海面静态特征图像能量谱为A(k_x,k_y)：

$A(k_x,k_y)=\sqrt{\left[\mathrm{Re}\left(F(k_x,k_y)\right)\right]+\left[\mathrm{Im}\left(F(k_x,k_y)\right)\right]}$

(3)基于尺度分离获得风条纹能量谱：导航雷达图像呈现的海面条纹波长L与能量谱波数k的关系式为：

$k=\frac{2\mathrm{\π}}{L}$

其中，(k_x,k_y)为能量谱域Τ的坐标，即k_x和k_y为k在x和y轴的波数分量；

设L_d为风条纹波长尺度下限，L_t为风条纹波长尺度上限，根据得到风条纹能量谱波数下限为k_d：

$k_d=\frac{2\mathrm{\π}}{L_d}$

风条纹能量谱波数上限为k_t：

$k_t=\frac{2\mathrm{\π}}{L_t}$

根据风条纹能量谱波数上下限设计能量谱带通滤波器，将风条纹信号能量谱提取出来，数学模型为：

I(k_x,k_y)为风条纹能量谱；

(4)海面风向提取：根据傅里叶变换具有周期性，得到的风条纹能量谱I(k_x,k_y)关于一三或二四象限镜像对称，取对称的两个能量集中区域，连线的垂直方向则为风条纹平行方向，得到海面风向；将计算出来的方向与风向标测得的风向所在象限比较，保留与其象限一致的方向去除180°模糊的方向；

对计算的风向进行校正后，得到相对正北方向的海面风向，校正公式为：

N_w＝|θ_c|+|α|+|β|

其中，N_w为相对正北向的风向；θ_c为选取区域的中心角；β为笛卡尔坐标下计算得到的风向；α为船艏向。

4.根据权利要求3所述的一种基于相对辐射校正的导航雷达图像反演海面风向方法，其特征是：风条纹波长尺度下限L_d为200m，风条纹波长尺度上限L_t为500m。