CN102681033A - 一种基于x波段航海雷达的海面风场测量方法 - Google Patents

Abstract

本发明公开了一种基于X波段航海雷达的海面风场测量方法，属于海洋动力环境遥感技术领域。所述的测量方法包含雷达图像前处理、风向测量和风速测量三部分，在风向测量指标中，将图像梯度、灰度和平滑项有机的结合起来，通过比例系数调节三者的比重，建立适合海面风场特征的模型，与现有技术相比，风向测量精度提高了68.4％以上。在风速测量指标中，当雷达单独测量时，将NRCS、实测风向、SNR作为BP网络输入，较传统算法风速精度提高了84％以上。在风速测量指标中，将海气边界层参数作为BP网络的附加输入，可进一步提高航海雷达测量风速的精度，考虑海气温差、盐度、潮位、大气压时，测量精度提高了48％以上。

Description

一种基于X波段航海雷达的海面风场测量方法

技术领域

本发明属于海洋动力环境遥感技术领域，特别涉及应用X波段航海雷达测量海面风场参数的方法。

背景技术

海面风场是海洋表面活动的主要动力来源，调节着海洋上层和低层大气之间的能量、动量、物质和水气的交换，直接影响着人类海上活动(海洋渔业、海上航行、海洋工程等)，它是海/气边界层以及生物圈系统的关键因素。

风掠过海面时，引起海面的粗糙不平，即微尺度波，它们与航海雷达的电磁波波长类似，当电磁波遇到微尺度波时，发生Bragg共振，产生后向回波，形成海杂波图像。2004年德国GKSS研究中心Dankert等人首次应用基于图像灰度不变模型(BCM)的光流法从海杂波图像序列中提取了海面风场，此方法简单，不需要校正雷达相位。见参考文献[1](王剑，段华敏.X波段雷达图像提取海洋表面风场.海洋技术，2010，29(3)：5-8页)和参考文献[2](Dankert H，Horstmann J，Rosenthal W.Ocean surface winds retrieved from marine radar image sequences.International Geosciences and Remote Sensing Symposium，2004，3：1903-1906P)，现有的BCM模型海面风场反演技术包含以下几个步骤：

步骤1、建立光流方程。假设短时间内图像灰度值不变，即满足下式：

G(X+f，t+Δt)＝G(X，t)

式中，G(X，t)为图像灰度值；X＝[x，y]^T为图像空间坐标；f＝[u，v]^T为光流；u、v为光流分量；Δt为两幅图的时间间隔。

上述的方程为非线性方程，为简化运算，常常将上述的方程按照一阶泰勒公式展开，令Δt→0，取极限并略去高阶项可得到基本约束方程，如下式：

G_xu+G_yv+G_t＝0

式中，

为灰度空间变化率，

为灰度时间变化率。

上式可简写为：

Gf+g＝0

式中，G＝(G_x，G_y)为灰度值空间导数；f＝(u，v)为光流；g＝G_t为灰度时间变化率；该式即为光流方程。

光流方程表明：海杂波图像中某一点的灰度时间变化率的大小是灰度空间导数与该点空间运动速度的乘积。

步骤2、应用加权最小二乘法求解光流方程，求取准则是：

J(f)＝(Gf+g)^TW(Gf+g)＝min

式中，W是适当取值的正定加权矩阵。

要使上述的求取准则方程成立，f应满足：

$\frac{\∂J}{\∂f}=G^{T}W(\mathrm{Gf}+g)+{\left[{(\mathrm{Gf}+g)}^T\mathrm{WG}\right]}^T=0$

从中解得：

f＝[G^T(W+W^T)G]^-1G^T(W+W^T)g

一般情况下，正定加权矩阵W取对称阵，即W＝W^T，所以加权最小二乘估计为：

f＝-(G^TWG)^-1G^TWg

为提高计算精度，取正定加权矩阵W为Sobel算子B。

G^TBGf＝-G^TBg

将上式展开：

$B\left[\begin{array}{c}\frac{\∂G}{\∂x}\\ \frac{\∂G}{\∂y}\end{array}\right]\left[\begin{array}{cc}\frac{\∂G}{\∂x}& \frac{\∂G}{\∂y}\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=-B\left[\begin{array}{c}\frac{\∂G}{\∂x}\\ \frac{\∂G}{\∂y}\end{array}\right]\frac{\∂G}{\∂t}$

$\left[\begin{array}{c}u\\ v\end{array}\right]=-{\left[\begin{array}{cc}B\frac{\∂G}{\∂x}\·\frac{\∂G}{\∂x}& B\frac{\∂G}{\∂x}\·\frac{\∂G}{\∂y}\\ B\frac{\∂G}{\∂y}\·\frac{\∂G}{\∂x}& B\frac{\∂G}{\∂y}\·\frac{\∂G}{\∂y}\end{array}\right]}^{-1}\left[\begin{array}{c}B\frac{\∂G}{\∂x}\·\frac{\∂G}{\∂t}\\ B\frac{\∂G}{\∂y}\·\frac{\∂G}{\∂t}\end{array}\right]$

记为如下形式：

$\left[\begin{array}{c}u\\ v\end{array}\right]=-{\left[\begin{array}{cc}B{\left(G_x\right)}^2& B{G}_x\·G_y\\ B{G}_y\·G_x& B{\left(G_y\right)}^2\end{array}\right]}^{-1}\left[\begin{array}{c}B{G}_x\·G_t\\ B{G}_y\·G_t\end{array}\right]$

假设在小邻域U内光流是常量，U内的像素点构成集合为{G₁，G₂，...，G_n}，运用加权最小二乘估计，此时有：

$G^T=\left[\begin{array}{cccc}\frac{\∂G_1}{\∂x}& \frac{\∂G_2}{\∂x}& \·\·\·& \frac{\∂G_n}{\∂x}\\ \frac{\∂G_1}{\∂y}& \frac{\∂G_2}{\∂y}& \·\·\·& \frac{\∂G_n}{\∂x}\end{array}\right],$ $g=\left[\begin{array}{cccc}\frac{\∂G_1}{\∂t}& \frac{\∂G_2}{\∂t}& \·\·\·& \frac{\∂G_n}{\∂t}\end{array}\right]$

那么光流分量的计算公式为：

$\left[\begin{array}{c}u\\ v\end{array}\right]=-{\left[\begin{array}{cc}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}B{\left(G_{\mathrm{ix}}\right)}^2& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}B{G}_{\mathrm{ix}}\·G_{\mathrm{iy}}\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}B{G}_{\mathrm{ix}}\·G_{\mathrm{iy}}& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}B{\left(G_{\mathrm{iy}}\right)}^2\end{array}\right]}^{-1}\left[\begin{array}{c}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}B{G}_{\mathrm{ix}}\·G_{\mathrm{it}}\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}B{G}_{\mathrm{iy}}\·G_{\mathrm{it}}\end{array}\right]$

式中 $G_{i*}=\frac{\∂G_i(X,t)}{\∂*}.$

步骤3、应用直方图统计得到的光流场，最大值即为主风向或主风速。

上述的反演方法中，由于海杂波图像序列具有较强的非线性，应用BCM模型计算海面风场会存在较大误差。

发明内容

为提高应用X波段航海雷达测量海面风场的测量精度，本发明公布了基于图像梯度、灰度不变模型(BGCM)能量函数的海面风场测量方法。

所述的海面风场测量方法包含三部分：雷达图像前处理、风向测量和风速测量，所述的风向测量具体步骤如下：

令I(X)表示三维海杂波图像序列，X＝[x，y，t]^T，[x，y]是空间区域Ω内坐标，t≥0代表时间序列，w＝[u，v，1]表示t时刻图像与t+1时刻图像间光流，u、v分别为光流的x，y分量，首先给出两种假设：

假设1：在图像时间间隔内，假设图像灰度不变，那么灰度值满足下式：

$\∀X\∈\mathrm{\Ω}:{|I(X+w)-I\left(X\right)|}^2=0---\left(35\right)$

其图像能量函数如式(36)：

E_D1(w)＝∫_ΩΨ(|I(X+w)-I(X)|²)dX                 (36)

式中，

ε＝10^-3为鲁棒函数，具有稳健性并抑制异常点；

假设2：当风场均匀变化或恒定时，此时雷达图像梯度不变，即：

$\∀X\∈\mathrm{\Ω}:{|\▿I(X+w)-\▿I\left(X\right)|}^2=0---\left(37\right)$

其图像能量函数满足下式：

$E_{D2}\left(w\right)={\∫}_{\mathrm{\Ω}}\mathrm{\Ψ}\left({|\▿I(X+w)-\▿I\left(X\right)|}^2\right)\mathrm{dX}---\left(38\right)$

式中，∨I＝(I_x，I_y)^T为图像空间导数，

为当前图像的空间导数，

为下时刻图像的空间导数，图像空间导数由二维最优化Sobel算子求得；

将式(36)、(38)模型有机的结合起来，得到数据项能量函数，如式(39)：

E_D(w)＝E_D1(w)+αE_D2(w)    (39)

式(39)中，α调节E_D1与E_D2的比重系数；

当上述两种假设不成立时，根据式(36)、(38)及(39)构造的基于BGCM模型的能量函数如下式：

E(w)＝E_D1(w)+αE_D2(w)+βE_S(w)    (41)

式中，β为平滑项调节系数，α、β值满足下式：

$\left\{\begin{array}{c}10\≤\mathrm{\β}/\mathrm{\α}\≤20\\ 10\≤\mathrm{\α}\≤20\end{array}\right.---\left(42\right)$

使用正则化方法求解光流w，光流w为下式的解：

$\left\{\begin{array}{c}\frac{\∂E\left(w\right)}{\∂u}=\frac{\∂E_{D1}\left(w\right)}{\∂u}+\mathrm{\α}\frac{\∂E_{D2}\left(w\right)}{\∂u}+\mathrm{\β}\frac{\∂E_S\left(w\right)}{\∂u}=0\\ \frac{\∂E\left(w\right)}{\∂v}=\frac{\∂E_{D1}\left(w\right)}{\∂v}+\mathrm{\α}\frac{\∂E_{D2}\left(w\right)}{\∂v}+\mathrm{\β}\frac{\∂E_S\left(w\right)}{\∂v}=0\end{array}\right.---\left(43\right)$

为公式表达方便，给出如下形式的简写：

$I_x={\∂}_{x}I(X+w)$

$I_y={\∂}_{y}I(X+w)$

I_z＝I(X+w)-I(X)

$I_{\mathrm{xx}}={\∂}_{\mathrm{xx}}I(X+w)$

$I_{\mathrm{xy}}={\∂}_{\mathrm{xy}}I(X+w)---\left(44\right)$

$I_{\mathrm{yy}}={\∂}_{\mathrm{yy}}I(X+w)$

$I_{\mathrm{xz}}={\∂}_{x}I(X+w)-{\∂}_{x}I\left(X\right)$

$I_{\mathrm{yz}}={\∂}_{y}I(X+w)-{\∂}_{y}I\left(X\right)$

求解式(43)得到欧拉-拉格朗日方程组如下式：

$\left\{\begin{array}{c}{\mathrm{\Ψ}}_{D1}^{\′}\left(I_z^2\right)\·I_z{I}_x+\mathrm{\α}{\mathrm{\Ψ}}_{D2}^{\′}(I_{\mathrm{xz}}^2+I_{\mathrm{yz}}^2)\·(I_{\mathrm{xz}}{I}_{\mathrm{xx}}+I_{\mathrm{yz}}+I_{\mathrm{xy}})+\mathrm{\β}\mathrm{div}({\mathrm{\Ψ}}_S^{\′}({\left|{\▿}_{3}u\right|}^2+{\left|{\▿}_{3}v\right|}^2)\·{\▿}_{3}u)=0\\ {\mathrm{\Ψ}}_{D1}^{\′}\left(I_z^2\right)\·I_z{I}_y+\mathrm{\α}{\mathrm{\Ψ}}_{D2}^{\′}(I_{\mathrm{xz}}^2+I_{\mathrm{yz}}^2)\·(I_{\mathrm{xz}}{I}_{\mathrm{xy}}+I_{\mathrm{yz}}{I}_{\mathrm{yy}})+\mathrm{\β}\mathrm{div}({\mathrm{\Ψ}}_S^{\′}({\left|{\▿}_{3}u\right|}^2+{\left|{\▿}_{3}v\right|}^2)\·{\▿}_{3}v)=0\end{array}\right.---\left(45\right)$

应用迭代法求解欧拉-拉格朗日非线性方程组，计算光流w；假设第k步迭代时，光流w^k＝[u^k，v^k，1]^T，令迭代初始条件w⁰＝[0，0，1]^T；将式(44)表示为

为消除

的非线性，根据泰勒公式，用一阶近似得下式：

$I_{*z}^{k+1}=I_{*}(X+w^{k+1})-I_{*}\left(X\right)$

$\≈I_{*}(X+w^k)+I_{*x}^{k}d{u}^k+I_{*y}^{k}d{v}^k-I_{*}\left(X\right)---\left(46\right)$

$=I_{*z}^k+I_{*x}^{k}d{u}^k+I_{*y}^{k}d{v}^k$

为公式表达方便，令：

$I_{\▿}={(I_x,I_y,I_z)}^T$

$S=I_{\▿}{I_{\▿}}^T---\left(47\right)$

$T=I_{\▿x}{I}_{\▿x}^T+I_{\▿y}{I}_{\▿y}^T$

经局部线性化后式(45)变为：

${\mathrm{\Ψ}}_{D1}^{\′k}\·(S_{11}^{k}d{u}^k+S_{12}^{k}d{v}^k+S_{13}^k)+\mathrm{\α}{\mathrm{\Ψ}}_{D2}^{\′k}(T_{11}^{k}d{u}^k+T_{12}^{k}d{v}^k+T_{13}^k)+\mathrm{\β}\mathrm{div}({\mathrm{\Ψ}}_S^{\′k}\▿(u^k+d{u}^k))=0$

${\mathrm{\Ψ}}_{D1}^{\′k}\·(S_{12}^{k}d{u}^k+S_{22}^{k}d{v}^k+S_{23}^k)+\mathrm{\α}{\mathrm{\Ψ}}_{D2}^{\′k}(T_{12}^{k}d{u}^k+T_{22}^{k}d{v}^k+T_{23}^k)+\mathrm{\β}\mathrm{div}({\mathrm{\Ψ}}_S^{\′k}\▿(u^k+d{v}^k))=0---\left(48\right)$

式中

${\mathrm{\Ψ}}_{D1}^{\′k}:={\mathrm{\Ψ}}^{\′}\left({(d{u}^k,{\mathrm{dv}}^k,1)}^T{S}^k(d{u}^k,{\mathrm{dv}}^k,1)\right)$

${\mathrm{\Ψ}}_{D2}^{\′k}:={\mathrm{\Ψ}}^{\′}\left({(d{u}^k,d{v}^k,1)}^T{T}^k({\mathrm{du}}^k,{\mathrm{dv}}^k,1)\right)---\left(49\right)$

${\mathrm{\Ψ}}_S^{\′k}:={\mathrm{\Ψ}}^{\′}({|\▿(u^k+d{u}^k)|}^2+{|\▿(v^k+d{v}^k)|}^2)$

对式(48)非线性系统进行离散化；第k次迭代时，离散化时的一些参数：k为迭代次数；图像总数目为m；像素大小图像I(X+w^k)中，像素集合{i|1 ，...，N^k}，N^k为总像素点数；

光流增量du^k，dv^k；在l方向上邻域N_l(i)，l∈{x，y}；经离散化后式(48)变为下式：

${\mathrm{\Ψ}}_{D1i}^{\′k}\·(S_{11i}^{k}d{u}_i^k+S_{12i}^{k}d{v}_i^k+S_{13i}^k)+\mathrm{\α}{\mathrm{\Ψ}}_{D2i}^{\′k}(T_{11i}^{k}d{u}_i^k+T_{12i}^{k}d{v}_i^k+T_{13i}^k)$

$+\mathrm{\β}\underset{l=x,y}{\mathrm{\Σ}}\underset{j\∈N_l\left(i\right)}{\mathrm{\Σ}}\frac{{\mathrm{\Ψ}}_{\mathrm{Si}}^{\′k}+{\mathrm{\Ψ}}_{\mathrm{Sj}}^{\′k}}{2}\·\frac{u_j^k+d{u}_j^k-u_i^k-d{u}_i^k}{{\left(h_l^k\right)}^2}=0$ (50)

${\mathrm{\Ψ}}_{D1}^{\′k}\·(S_{12i}^{k}d{u}_i^k+S_{22i}^{k}d{v}_i^k+S_{23i}^k)+\mathrm{\α}{\mathrm{\Ψ}}_{D2i}^{\′k}(T_{12i}^{k}d{u}_i^k+T_{22i}^{k}d{v}_i^k+T_{23i}^k)$

$+\mathrm{\β}\underset{l=x,y}{\mathrm{\Σ}}\underset{j\∈N_l\left(i\right)}{\mathrm{\Σ}}\frac{{\mathrm{\Ψ}}_{\mathrm{Si}}^{\′k}+{\mathrm{\Ψ}}_{\mathrm{Sj}}^{\′k}}{2}\·\frac{v_j^k+d{v}_j^k-v_i^k-d{v}_i^k}{{\left(h_l^k\right)}^2}=0$

为求解式(50)，需要经过两次迭代过程，迭代2过程嵌入在迭代1过程中，具体如下：

(1)迭代初始化，令光流分量和光流增量均为零，即(u⁰，v⁰)＝(0，0)，(du⁰，dv⁰)＝(0，0)；

(2)根据下式求解光流增量：

$\left(\begin{array}{c}{\mathrm{du}}^{k,n+1}\\ {\mathrm{dv}}^{k,n+1}\end{array}\right)={\left(\begin{array}{cc}{M}_{11}^{k,n}& M_{12}^{k,n}\\ M_{21}^{k,n}& M_{22}^{k,n}\end{array}\right)}^{-1}\left(\begin{array}{c}{\mathrm{ru}}^{k,n}\\ {\mathrm{rv}}^{k,n}\end{array}\right)---\left(51\right)$

式中，du^k，n，dv^k，n为第k次迭代1、第n次迭代2时的光流增量，

ru^k，n，rv^k，n定义如下：

$M_{11}^{k,n}={\mathrm{\Ψ}}_{D1i}^{\′k,n}{S}_{11i}^{k,n}+\mathrm{\α}{\mathrm{\Ψ}}_{D2i}^{\′k,n}{T}_{11i}^{k,n}+\mathrm{\β}\underset{l=x,y}{\mathrm{\Σ}}\underset{j\∈N_l\left(i\right)}{\mathrm{\Σ}}\frac{{\mathrm{\Ψ}}_{\mathrm{Si}}^{\′k,n}+{\mathrm{\Ψ}}_{\mathrm{Sj}}^{\′k,n}}{2{\left(h_l^k\right)}^2}$

$M_{22}^{k,n}={\mathrm{\Ψ}}_{D1i}^{\′k,n}{S}_{22i}^{k,n}+\mathrm{\α}{\mathrm{\Ψ}}_{D2i}^{\′k,n}{T}_{22i}^{k,n}+\mathrm{\β}\underset{l=x,y}{\mathrm{\Σ}}\underset{j\∈N_l\left(i\right)}{\mathrm{\Σ}}\frac{{\mathrm{\Ψ}}_{\mathrm{Si}}^{\′k,n}+{\mathrm{\Ψ}}_{\mathrm{Sj}}^{\′k,n}}{2{\left(h_l^k\right)}^2}---\left(52\right)$

$M_{12}^{k,n}={\mathrm{\Ψ}}_{D1i}^{\′k,n}{S}_{12i}^{k,n}+\mathrm{\α}{\mathrm{\Ψ}}_{D2i}^{\′k,n}{T}_{12i}^{k,n}$

$-{\mathrm{\Ψ}}_{D1i}^{\′k,n}{S}_{13i}^{k,n}-\mathrm{\α}{\mathrm{\Ψ}}_{D2i}^{\′k,n}{T}_{13i}^{k,n}+\mathrm{\β}\underset{l=x,y}{\mathrm{\Σ}}\underset{j\∈N_l^{-}\left(i\right)}{\mathrm{\Σ}}\frac{{\mathrm{\Ψ}}_{\mathrm{Si}}^{\′k,n}+{\mathrm{\Ψ}}_{\mathrm{Sj}}^{\′k,n}}{2}\frac{u_j^k+d{u}_j^{k+1}-u_i^k}{{\left(h_l^k\right)}^2}$

$+\mathrm{\β}\underset{l=x,y}{\mathrm{\Σ}}\underset{j\∈N_l^{+}\left(i\right)}{\mathrm{\Σ}}\frac{{\mathrm{\Ψ}}_{\mathrm{Si}}^{\′k,n}+{\mathrm{\Ψ}}_{\mathrm{Sj}}^{\′k,n}}{2}\frac{v_j^k+d{v}_j^{k,n}-v_i^k}{{\left(h_l^k\right)}^2}=r{u}^{k,n}$ (53)

$-{\mathrm{\Ψ}}_{D1i}^{\′k,n}{S}_{23i}^{k,n}-\mathrm{\α}{\mathrm{\Ψ}}_{D2i}^{\′k,n}{T}_{23i}^{k,n}+\mathrm{\β}\underset{l=x,y}{\mathrm{\Σ}}\underset{j\∈N_l^{-}\left(i\right)}{\mathrm{\Σ}}\frac{{\mathrm{\Ψ}}_{\mathrm{Si}}^{\′k,n}+{\mathrm{\Ψ}}_{\mathrm{Sj}}^{\′k,n}}{2}\frac{u_j^k+d{u}_j^{k+1}-u_i^k}{{\left(h_l^k\right)}^2}$

$+\mathrm{\β}\underset{l=x,y}{\mathrm{\Σ}}\underset{j\∈N_l^{+}\left(i\right)}{\mathrm{\Σ}}\frac{{\mathrm{\Ψ}}_{\mathrm{Si}}^{\′k,n}+{\mathrm{\Ψ}}_{\mathrm{Sj}}^{\′k,n}}{2}\frac{v_j^k+d{v}_j^{k,n}-v_i^k}{{\left(h_l^k\right)}^2}=r{v}^{k,n}$

式中， $N_l^{-}\left(i\right)=\{j\&Element;N_l\left(i\right)|j<i\}$ 表示已经迭代更新的像素点， $N_l^{+}\left(i\right)=\{j\&Element;N_l\left(i\right)|j>i\}$ 表示未迭代更新的点；

(3)根据相对误差判断结果是否满足要求，若相对误差小于2％则满足要求，输出光流增量；否则，将最新光流增量作为已知，继续迭代2过程；

(4)根据下式计算该像素点的光流分量：

u^k+1＝u^k+du^k，n+1    (54)

v^k+1＝v^k+dv^k，n+1

(5)根据相对误差判断光流结果是否满足要求：若相对误差小于2％则满足要求，输出此像素点光流，否则判断迭代1的次数是否已达上限；若迭代1次数未达到上限则将最新光流作为已知，继续迭代1过程；若已达上限则采取高分辨率到低分辨率策略，应用高斯低通滤波器平滑，下采样，得到低分辨率图像，使用新图像进行迭代，完成存在大偏移量时风向的提取；

(6)对得到的光流场进行直方图统计，求取频率最大值C，将频率＞0.95C的角度值矢量平均得到主风向。

所述的风速测量包括两种方式的测量方法，第一种为雷达图像单独测量风速，第二种为采用辅助信息测量风速；第一种，雷达单独测量风速时，与风速有关的信息包括归一化的雷达散射截面、实测风向、海浪的信噪比，应用BP网络确定它们之间的关系，再根据训练好的BP网络计算风速。第二种为多信息辅助雷达测量风速。

本发明的优点在于：

(1)、在风向测量指标中，将图像梯度、灰度和平滑项有机的结合起来，通过比例系数调节三者的比重，建立适合海面风场特征的模型，通过实验结果可知，与现有技术相比，风向测量精度提高了68.4％以上。

(2)、在风速测量指标中，当雷达单独测量时，将NRCS、实测风向、SNR作为BP网络输入，通过实验结果可知，较传统算法，风速精度提高了84％以上。

(3)、在风速测量指标中，当存在海气边界层参数时，让它们作为BP网络的附加输入，可进一步提高航海雷达测量风速的精度。如考虑海气温差、盐度、潮位、大气压时，精度提高了48％以上。

附图说明

图1为本发明中的迭代法求解光流场的流程图；

图2为本发明是图1中的迭代2过程示意图；

图3a、3b为本发明中改进的BGCM与现有的BCM模型分别与参考风向比对散点图；

图4a、4b为本发明中BP网络不同阶段雷达单独测量风速散点图；

图5为本发明中多信息辅助雷达测量风速散点图；

图6为本发明提供的海面风场测量方法流程图；

图7为本发明中海流计算流程图；

图8为BP网络的输入输出示意图。

具体实施方式

下面结合附图对本发明提供的基于X波段航海雷达的海面风场测量方法详细介绍。

本发明提供的基于X波段航海雷达的海面风场测量方法，如图6所示流程图，本发明提供的海面风场测量方法包括雷达图像的前处理、风向测量和风速测量三部分内容，下述步骤中第一步～第八步为雷达图像的前处理，第九步为风向测量，第十步为风速测量，具体步骤如下：

第一步，雷达图像方位和角度校正，以及应用中值滤波器抑制噪声。

1.1、雷达图像方位和角度校正：雷达图像是极坐标图像，由于雷达信号强度在距离R上衰减满足式(1)，且回波强度随着方位角θ不同有所不同，满足式(2)，因此在对雷达图像进行风速和风向测量之前，首先需要对雷达图像进行方位和角度校正。所述的方位和角度校正就是要消除方位和角度对回波强度的影响(具体消除方法可参见参考文献MIROS A/S.WAVEX4.0Technical Handbook，2002)。

σ⁰∝R^-n，n∈[2，4]           (1)

σ⁰∝A₂+B₂*cosθ+C₂*cos2θ    (2)

式(2)中，A₂、B₂、C₂为常数，需根据实测数据标定，每部雷达均不同；R为雷达天线与目标之间的距离，σ⁰是回波强度。

1.2、中值滤波器抑制噪声：应用相邻像素的灰度中值来替代该像素的灰度值：

$I(x,y)=\underset{(s,t)\&Element;S_{\mathrm{xy}}}{\mathrm{median}}\left\{g(s,t)\right\}---\left(3\right)$

式(3)中，S_xy表示中心在(x，y)点，尺寸为m×n的矩形窗口的坐标组，g(s，t)为雷达图像，I(x，y)为中值滤波后图像的灰度值。

第二步，在第一步处理校正和滤波后的雷达图像中选取测量区域。由于岸边以及建筑物遮挡的影响，并不是整幅雷达图像均含有丰富的海杂波信号，因此需要选取感兴趣的区域，称为子图像。选取测量区域的原则是：海浪在时间和空间上分布均匀，并且海浪受地形影响较小。为便于处理，一般应用最近点插值方法将选取的测量区域图像由极坐标转换为直角坐标。多幅子图像组成子图像序列σ(x，y，t)，总时间长度在40～100s范围。

第三步，确定三维波数频率图像谱：对子图像序列σ(x，y，t)做三维离散傅里叶变换如下式所示，变换后得到三维波数频率图像谱：

$I^{\left(3\right)}(k_x,k_y,\mathrm{\&omega;})=\underset{y=0}{\overset{L_y}{\mathrm{\&Sigma;}}}\underset{x=0}{\overset{L_x}{\mathrm{\&Sigma;}}}\underset{t=0}{\overset{T}{\mathrm{\&Sigma;}}}\mathrm{\&sigma;}(x,y,t)\exp [-2\mathrm{\&pi;i}(k_{x}x/L_x+k_{y}y/L_y+\mathrm{\&omega;t}/T)]---\left(4\right)$

式(4)中，L_x、L_y为子图像空间大小，T为总时间长度；k为波数，k_x、k_y分别为波数k在x和y方向的波数分量。图像功率谱为频谱的平方，可通过下式求取：

P(k_x，k_y，ω)＝|I⁽³⁾(k_x，k_y，ω)|²＝R_e ²(k_x，k_y，ω)+I_m ²(k_x，k_y，ω)    (5)

式(5)中R_e、I_m分别代表I⁽³⁾(k_x，k_y，ω)的实部和虚部。

图像谱分辨率受图像时间和空间尺度限制，即：

$\mathrm{\&Delta;}{k}_x=\frac{2}{L_x},$ $\mathrm{\&Delta;}{k}_y=\frac{2\mathrm{\&pi;}}{L_y},$ $\mathrm{\&Delta;\&omega;}=\frac{2\mathrm{\&pi;}}{T}---\left(6\right)$

式(6)中，Δk_x与Δk_y为波数分量分辨率；Δω为频率分辨率。

第四步，计算海表层流(简称海流)：

结合图7海流计算流程图海流计算做详细说明：

4.1、对于第三步中的三维波数频率图像谱，应用如下形式的色散关系带通滤波器抑制噪声干扰：

$\left|k_p\right|=\frac{{(\mathrm{\&omega;}+\frac{\mathrm{\&Delta;\&omega;}}{2}+|k_m|\&CenterDot;|U_{\max }\left|\right)}^2}{g}+\frac{\mathrm{\&Delta;k}}{\sqrt{2}}---\left(8\right)$

$\left|k_n\right|=\frac{{(\mathrm{\&omega;}-\frac{\mathrm{\&Delta;\&omega;}}{2}-|k_m|\&CenterDot;|U_{\max }\left|\right)}^2}{g}-\frac{\mathrm{\&Delta;k}}{\sqrt{2}}---\left(9\right)$

上式中，|k_n|、|k_p|为带通滤波器频带边界；Δω为频率分辨率，ω为海浪频率，Δk为波数分辨率；U_max为粗略估计的雷达天线与海浪场间的最大相对速度，当雷达载体静止时，U_max为粗略估计的最大海表面流速，在文献[4](2012年华中科技大学学报(自然科学版)期刊，袁赣南等人的文章色散关系带通滤波器设计)中指出，雷达静止时，U_max取值为1.5m/s，考虑到每个海域海流速度不同，本发明中要求U_max取值在1.2～2m/s；F⁽³⁾(k，ω)为滤波后海浪图像谱。k_m为主导波数，根据下面方法确定：

将滤波后海浪图像谱F⁽³⁾(k，ω)按波数角度进行积分，得到二维波数模频率域海浪图像谱I⁽²⁾(|k|，ω)，如下式：

$I^{\left(2\right)}\left(\right|k|,\mathrm{\&omega;})={\&Integral;}_0^{2\mathrm{\&pi;}}{F}^{\left(3\right)}\left(\right|k|,{\mathrm{\&theta;}}_1,\mathrm{\&omega;})d{\mathrm{\&theta;}}_1---\left(10\right)$

式(10)中，F⁽³⁾(|k|，θ₁，ω)＝F⁽³⁾(k，ω)为三维波数频率图像谱；θ₁为海浪谱波数角度。在I⁽²⁾(|k|，ω)中，C为最大能量值，选取能量大于0.95*C的能量点，将它们的波数模进行算术平均，即得到主导波数|k_m|。

4.2、自适应选取阈值构造低通滤波器，抑制能量低于阈值的噪声。

首先，根据滤波后海浪图像谱F⁽³⁾(k，ω)和信噪比snr计算海浪谱总能量P，如下式：

$P=\frac{\mathrm{snr}+1}{\mathrm{snr}}\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}\underset{k}{\overset{N_{\mathrm{\&omega;}}}{\mathrm{\&Sigma;}}}{F}^{\left(3\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}},{\mathrm{\&omega;}}_k)\&CenterDot;\mathrm{\&Delta;}{k}_x\mathrm{\&Delta;}{k}_y\mathrm{\&Delta;\&omega;}---\left(11\right)$

式中：

$\mathrm{snr}=\frac{\mathrm{sig}}{\mathrm{bgn}}---\left(12\right)$

$\mathrm{sig}=\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}{F}^{\left(2\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}})\mathrm{\&Delta;}{k}_x\mathrm{\&Delta;}{k}_y---\left(13\right)$

$\mathrm{bgn}=\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}\underset{l}{\overset{N_{\mathrm{\&omega;}}}{\mathrm{\&Sigma;}}}{I}^{\left(3\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}},{\mathrm{\&omega;}}_l)\mathrm{\&Delta;}{k}_x\mathrm{\&Delta;}{k}_y\mathrm{\&Delta;\&omega;}-\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}{F}^{\left(2\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}})\mathrm{\&Delta;}{k}_x\mathrm{\&Delta;}{k}_y---\left(14\right)$

$F^{\left(2\right)}\left(k\right)=2{\&Integral;}_{\mathrm{\&omega;}>0}{I}^{\left(3\right)}(k_x,k_y,\mathrm{\&omega;})\mathrm{\&delta;}(\mathrm{\&omega;}-{\stackrel{\&OverBar;}{\mathrm{\&omega;}}}_0)\mathrm{d\&omega;}---\left(15\right)$

式中，sig为海浪图像谱能量，bgn为背景噪声能量，I⁽³⁾(k，ω)是三维图像谱，F⁽²⁾(k)是二维海浪图像谱，N_kx、N_ky、N_ω为三维波数频率图像谱的波数在x方向和y方向的范围，以及频率的范围，Δk_x、Δk_y为波数分辨率，Δω为频率分辨率，

是色散关系带通滤波器

然后，对三维波数频率图像谱中的所有的能量点按照从大到小的顺序进行排序，排序后第i个能量点的能量为P_i，由于海浪能量较背景噪声能量大，前m个能量点的总和为P，即：

$P=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{P}_i---\left(16\right)$

最后，由上式即可求出数目m，第m个能量点的能量值为P_m，那么设P_m为阈值C_it。

根据上述的阈值C_it构造低通滤波器，选取能量大于阈值C_it的能量点，得到低通滤波后的三维波数频率图像谱。

4.3、构造目标函数。

(1)、计算海浪谱

由于受雷达成像过程非线性和海浪本身弱非线性的影响，海杂波图像谱与真实海浪谱间存在较大差异，大量论文均指出可应用调制传递函数(MTF)对三维波数频率图像谱进行非线性校正，得到海浪谱

MTF定义如下式：

MTF＝|k|^β    (17)

MTF校正方法如下式：

$F_M^{\left(3\right)}(k,\mathrm{\&omega;})={\left|k\right|}^{\mathrm{\&beta;}}\&CenterDot;I^{\left(3\right)}(k,\mathrm{\&omega;})---\left(18\right)$

式中，F⁽³⁾(k，ω)为图像谱，图像谱经调制传递函数校正后变为海浪谱

MTF的参数β需根据长时间大量的实验数据确定，它取决于雷达成像的非线性程度，对于同一部雷达只需要标定一次即可，在文献[5](崔利民、何宜军，IEEE International Geoscience and RemoteSensing Symposium，2009年一文Measurement of ocean wave spectra with vertical polarizationx-band radar imgae sequences)中，作者标定出X波段雷达水平极化下β取值1.17；在文献[6](周蓓硕士论文，2008年一文X波段雷达海表面流信息提取技术研究)中，标定出X波段雷达水平极化下β取值1.21；在文献[7](Lzquierdo P等，Jouranl of Waterway，Port，Coastal andOcean Engineering期刊，2005年一文Comparison of wave spectra from nautical radar images andscalar buoy data)中使用β值1.2。本发明中，考虑到β值取决于雷达成像机制的非线性，要求β值范围大于等于117小于等于1.21。

(2)、计算隶属度。根据谱分量到中心的距离来度量隶属度函数大小，谱分量距离中心越近，隶属度越大，反之，则越小。对于固定频率ω_i上有N个谱分量组成集合{k₁ k₂...k_N}，那么中心k_i0，最大半径r_i定义如下：

$k_{i0}=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{k}_j---\left(19\right)$

$r_i=\underset{k_j}{\max }\left|\right|k_j-k_{i0}\left|\right|---\left(20\right)$

式中，k_j为谱分量；i＝1，2，...，N；i＝1，2...，N_ω。根据所述的谱分量到中心的距离确定隶属度时，集合中谱分量隶属度为模糊隶属度μ(k_ij)：

$\mathrm{\&mu;}\left(k_{\mathrm{ij}}\right)=1-\frac{\left|\right|k_j-k_{i0}\left|\right|}{r_i+\mathrm{\&delta;}}---\left(21\right)$

式中，δ为很小的正数，目的是避免出现μ(k_ij)＝0的情况。

(3)、构造目标函数。将海浪谱

和模糊隶属度μ(k_ij)的乘积作为最小二乘法的权值，那么海流反演计算的目标函数为：

式中，ω_i为实测数据海浪谱

中第i个频率分量，ω_p为理论上色散关系方程给出的频率分量；N_ω为三维波数频率图像谱的频率的范围，N为ω_i频率上具有的谱分量数目；

μ(k_ij)分别为ω_i频率上第j个海浪谱和隶属度。

4.4、海流计算：

(1)、初始估流：在此仅考虑0阶次波，根据0阶色散关系方程计算海流，步骤与传统算法的初始估流类似，具体如下：

①取阈值C_fg，使得能量大于阈值C_fg的谱数量n₂，一般选择阈值C_fg的量值为0.2或n₂值在数十个。

②将式(22)的目标函数Q²对海流U的两个分量U_x、U_y求偏导，并使其为零，得到海流的计算公式如下：

U＝D·B    (23)

$D={\left[\begin{array}{cc}\mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x^2& \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x{k}_y\\ \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x{k}_y& \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_y^2\end{array}\right]}^{-1}---\left(24\right)$

$B=\left[\begin{array}{c}\mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x({\mathrm{\&omega;}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k|d})\\ \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_y({\mathrm{\&omega;}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k|d})\end{array}\right]---\left(25\right)$

式中，F_M为海浪谱；μ为隶属度；k_x、k_y为k_ij的分量；ω₀为0阶波频率，d为海深，其它变量的定义同前。经过初始估流阶段可得到粗略精度的海流。

(2)、迭代估流：根据0、1阶次波(大于1阶的海浪能量很少)确定海流。

①应用实时色散关系带通滤波器(如式(30))计算信噪比snr；然后应用自适应阈值选取技术得到阈值C_it，选取大于阈值C_it的图像谱，数目为n₃；再根据初始估流得到的较低精度的海流，应用色散关系方程式(26)，计算0、1阶次波频率：

${\mathrm{\&omega;}}_p=(p+1)\sqrt{\frac{g\left|k\right|}{p+1}\tanh \left(\frac{\left|k\right|d}{p+1}\right)}+\left|k\right|\left|u\right|\cos \left(\mathrm{\&theta;}\right)---\left(26\right)$

式中，p为海浪能量的阶次，p＝0时，该式就变成了基本色散关系方程，p＞1时的海浪被称作高阶次波，ω_p为p阶次波的频率。

②对这n₃个图像谱进行如下讨论：若图像谱均在[0，ω_N]范围内，再判断它们属于0阶次还是1阶次波，若|ω_i-ω₀(k_i)|＜|ω_i-ω₁(k_i)|则属于0阶次波，若|ω_i-ω₀(k_i)|＞|ω_i-ω₁(k_i)|则属于1阶次波。若有图像谱超出[0，ω_N]频率范围的，需要按照傅里叶变换的性质(原点对称性和偶数性)对信号进行重构，然后再判断它属于0阶还是1阶次波。其中ω_N为奈奎斯特频率。将判断好的数据根据下式计算海流：

U＝D·B    (27)

$D={\left[\begin{array}{cc}\mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x^2& \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x{k}_y\\ \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x{k}_y& \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_y^2\end{array}\right]}^{-1}---\left(28\right)$

$B=\left[\begin{array}{c}\underset{n_{30}}{\mathrm{\&Sigma;}}{F}_M\mathrm{\&mu;}{k}_x({\mathrm{\&omega;}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k\left|\right)d})+\underset{n_{31}}{\mathrm{\&Sigma;}}{F}_M\mathrm{\&mu;}{k}_x({\mathrm{\&omega;}}_1-2\sqrt{\frac{g\left|k\right|}{2}\tanh (\frac{\left|k\right|d}{2}})\\ \underset{n_{30}}{\mathrm{\&Sigma;}}{F}_M\mathrm{\&mu;}{k}_y({\mathrm{\&omega;}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k|d})+\underset{n_{31}}{\mathrm{\&Sigma;}}{F}_M\mathrm{\&mu;}{k}_y({\mathrm{\&omega;}}_1-2\sqrt{\frac{g\left|k\right|}{2}\tanh (\frac{\left|k\right|d}{2}})\end{array}\right]---\left(29\right)$

式中，n₃₀、n₃₁分别为n₃个谱分量中0、1阶次波数目，ω₀、ω₁分别为0、1阶次波频率，其它变量定义同前。

③一般迭代10至15次可满足工程上精度需要，判断是否迭代结束(迭代次数是否达到要求)，如果没有结束则将新表面流代入式(26)，得到新的0阶次和1阶次波频率，返回步骤①；若迭代结束则停止计算，输出海流y。

第五步，在确定实时海流后，针对应用实时色散关系带通滤波器抑制海浪噪声对海面风场的干扰，滤波器形式如下：

式中，k_mn为实测波数；F⁽³⁾(k_mn，ω_v)为图像谱，P_D ⁽³⁾(k_mn，ω_v)为滤波后图像谱；从上式可看出，频率为ω_v时，滤波器带宽为δk_v，带宽是由频率决定的，在频带外部的能量被保留下来，其它能量均被抑制掉。

δk_v＝|k(ω_v+1)|-|k(ω_v-1)|；v＝0，...，N_t/2    (31)

式中，N_t为图像数目；频率ω_v时理论上的波数为k(ω_v)。

第六步，对滤波后的海浪图像谱P_D ⁽³⁾(k_mn，ω_v)进行三维反傅里叶变换，得到滤波后图像序列σ′(x，y，t)，如下式：

${\mathrm{\&sigma;}}^{\&prime;}(x,y,t)=\frac{1}{L_x{L}_{y}T}\underset{y=0}{\overset{L_y}{\mathrm{\&Sigma;}}}\underset{x=0}{\overset{L_x}{\mathrm{\&Sigma;}}}\underset{t=0}{\overset{T}{\mathrm{\&Sigma;}}}{P}_D(k_x,k_y,\mathrm{\&omega;})\exp \left[2\mathrm{\&pi;i}(k_{x}x/L_x+k_{y}y/L_y+\mathrm{\&omega;t}/T)\right]---\left(32\right)$

式中，P_D(k_x，k_y，ω)为滤波后图像谱。

第七步，对滤波后图像序列σ′(x，y，t)进行滑动平均，抑制在时间上变化频率较高的噪声，保留静态信号和周期较大的信号，形式如下：

$G(x,y,,\mathrm{\&tau;})=\frac{1}{n}\underset{t=i}{\overset{i+n-1}{\mathrm{\&Sigma;}}}{\mathrm{\&sigma;}}^{\&prime;}(x,y,t)---\left(33\right)$

式中，G(x，y，τ)为滑动平均后新序列，新序列图像数量为N_t-n+1，B_t为图像数目，n为参与平均的图像数量，τ为滑动平均后新序列的图像的时间。

第八步，在时间上，第七步采取措施抑制了噪声干扰，本步骤从图像空间角度采取措施抑制噪声干扰。应用低通滤波器进行图像的缩减和平滑，抑制空间尺度小的噪声，低通滤波器设计步骤如下：

(1)用4阶二项式系数卷积核B⁴去平滑图像；

(2)对子图像取2×2的抽样平均；

(3)用2阶二项式系数卷积核B²去平滑图像。

上述操作可用算子R_|2＝B²S_|2B⁴来完成，S_|2表示图像2×2抽样平均，经图像缩减后尺寸为原图像1/4。

$B^2=\frac{1}{16}\left[\begin{array}{ccc}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right];$ $B^4(m,n)=\frac{1}{256}\left[\begin{array}{ccccc}1& 4& 6& 4& 1\\ 4& 16& 24& 16& 4\\ 6& 24& 36& 24& 6\\ 4& 16& 24& 16& 4\\ 1& 4& 6& 4& 1\end{array}\right]---\left(34\right)$

图像缩减次数可根据图像分辨率要求确定，图像空间分辨率7.5×7.5m，若要得到60×60m的分辨率，图像缩减次数为3次，即(R_|2)³。

第九步，对于经过上述的前处理过的雷达图像，进行海面风向测量。

令I(X)表示三维海杂波图像序列，X＝[x，y，t]^T，[x，y]是空间区域Ω内坐标，t≥0代表时间序列，w＝[u，v，1]表示t时刻图像与t+1时刻图像间光流，u、v分别为光流的x，y分量。首先给出两种假设：

假设1：在图像时间间隔内，假设图像灰度不变，那么灰度值满足下式：

$\&ForAll;X\&Element;\mathrm{\&Omega;}:{|I(X+w)-I\left(X\right)|}^2=0---\left(35\right)$

其图像能量函数如式(36)：

E_D1(w)＝∫_ΩΨ(|I(X+w)-I(X)|²)dX    (36)

式中，

ε＝10^-3为鲁棒函数，具有稳健性并抑制异常点。

假设2：当风场均匀变化或恒定时，此时雷达图像梯度不变，即：

$\&ForAll;X\&Element;\mathrm{\&Omega;}:{|\&dtri;I(X+w)-\&dtri;I\left(X\right)|}^2=0---\left(37\right)$

其图像能量函数满足下式：

$E_{D2}\left(w\right)={\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&Psi;}\left({|\&dtri;I(X+w)-\&dtri;I\left(X\right)|}^2\right)\mathrm{dX}---\left(38\right)$

式中，∨I＝(I_x，I_y)^T为图像空间导数，

为当前图像的空间导数，

为下时刻图像的空间导数，图像空间导数可由二维最优化Sobel算子求得。

假设2很好的克服了假设1对于雷达回波变化过于敏感问题，同时能够提取出静态或均匀变化目标。

将式(36)、(38)模型有机的结合起来，得到数据项能量函数，如式(39)：

E_D(w)＝E_D1(w)+αE_D2(w)    (39)

式(39)中，α可调节E_D1与E_D2比重。

当图像中噪声或其它异常情况引起上述两种假设不成立时，再使用式(5)能量函数计算会存在较大误差，因此需要能量平滑项函数，考虑到光流场的连续性，将局部非平稳运动包含进来，1992年Rudin等人给了如下平滑项：

$E_S\left(w\right)={\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&Psi;}({\left|{\&dtri;}_{3}u\right|}^2+{\left|{\&dtri;}_{3}v\right|}^2)\mathrm{dX}---\left(40\right)$

根据式(36)、(38)及(39)构造的基于BGCM模型的能量函数如下式：

E(w)＝E_D1(w)+αE_D2(w)+βE_S(w)    (41)

式中，β为平滑项调节系数。在文献[3](2005年ICCV论文集Bruhn A的文章Towardsultimate motion estimation combing highest accuracy with real time performance)中指出α值为16.5，β为160，本发明中考虑到梯度、灰度、平滑能量项与海面风场的关系，要求αβ值满足下式：

$\left\{\begin{array}{c}10\&le;\mathrm{\&beta;}/\mathrm{\&alpha;}\&le;20\\ 10\&le;\mathrm{\&alpha;}\&le;20\end{array}\right.---\left(42\right)$

在式(41)中有两个未知数(w＝[u，v])，然而只有一个方程，因此方程的解不是唯一的，这个问题被称作是光流计算的瓶颈问题。解决这个问题的方法主要有两类：正则化法(全局法)和局部邻域法。

本发明中使用正则化方法求解光流w，正则化法是将光流方程的求解过程转换为一全局能量泛函的最小化过程，光流矢量w应满足使E(w)能量最小，即w为下式的解：

$\left\{\begin{array}{c}\frac{\&PartialD;E\left(w\right)}{\&PartialD;u}=\frac{\&PartialD;E_{D1}\left(w\right)}{\&PartialD;u}+\mathrm{\&alpha;}\frac{\&PartialD;E_{D2}\left(w\right)}{\&PartialD;u}+\mathrm{\&beta;}\frac{\&PartialD;E_S\left(w\right)}{\&PartialD;u}=0\\ \frac{\&PartialD;E\left(w\right)}{\&PartialD;v}=\frac{\&PartialD;E_{D1}\left(w\right)}{\&PartialD;v}+\mathrm{\&alpha;}\frac{\&PartialD;E_{D2}\left(w\right)}{\&PartialD;v}+\mathrm{\&beta;}\frac{\&PartialD;E_S\left(w\right)}{\&PartialD;v}=0\end{array}\right.---\left(43\right)$

为公式表达方便，给出如下形式的简写：

$I_x={\&PartialD;}_{x}I(X+w)$

$I_y={\&PartialD;}_{y}I(X+w)$

I_z＝I(X+w)-I(X)

$I_{\mathrm{xx}}={\&PartialD;}_{\mathrm{xx}}I(X+w)$ (44)

$I_{\mathrm{xy}}={\&PartialD;}_{\mathrm{xy}}I(X+w)$

$I_{\mathrm{yy}}={\&PartialD;}_{\mathrm{yy}}I(X+w)$

$I_{\mathrm{xz}}={\&PartialD;}_{x}I(X+w)-{\&PartialD;}_{x}I\left(X\right)$

$I_{\mathrm{yz}}={\&PartialD;}_{y}I(X+w)-{\&PartialD;}_{y}I\left(X\right)$

求解式(43)可得到欧拉-拉格朗日方程组如下式：

$\left\{\begin{array}{c}{\mathrm{\&Psi;}}_{D1}^{\&prime;}\left(I_z^2\right)\&CenterDot;I_z{I}_x+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2}^{\&prime;}(I_{\mathrm{xz}}^2+I_{\mathrm{yz}}^2)\&CenterDot;(I_{\mathrm{xz}}{I}_{\mathrm{xx}}+I_{\mathrm{yz}}+I_{\mathrm{xy}})+\mathrm{\&beta;}\mathrm{div}({\mathrm{\&Psi;}}_S^{\&prime;}({\left|{\&dtri;}_{3}u\right|}^2+{\left|{\&dtri;}_{3}v\right|}^2)\&CenterDot;{\&dtri;}_{3}u)=0\\ {\mathrm{\&Psi;}}_{D1}^{\&prime;}\left(I_z^2\right)\&CenterDot;I_z{I}_y+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2}^{\&prime;}(I_{\mathrm{xz}}^2+I_{\mathrm{yz}}^2)\&CenterDot;(I_{\mathrm{xz}}{I}_{\mathrm{xy}}+I_{\mathrm{yz}}{I}_{\mathrm{yy}})+\mathrm{\&beta;}\mathrm{div}({\mathrm{\&Psi;}}_S^{\&prime;}({\left|{\&dtri;}_{3}u\right|}^2+{\left|{\&dtri;}_{3}v\right|}^2)\&CenterDot;{\&dtri;}_{3}v)=0\end{array}\right.---\left(45\right)$

应用迭代法求解欧拉-拉格朗日非线性方程组，计算光流w。假设第k步迭代时，光流w^k＝[u^k，v^k，1]^T，令迭代初始条件w⁰＝[0，0，1]^T；将式(44)表示为

为消除

的非线性，根据泰勒公式，用一阶近似得下式：

$I_{*z}^{k+1}=I_{*}(X+w^{k+1})-I_{*}\left(X\right)$

$\&ap;I_{*}(X+w^k)+I_{*x}^{k}d{u}^k+I_{*y}^{k}d{v}^k-I_{*}\left(X\right)---\left(46\right)$

$=I_{*z}^k+I_{*x}^{k}d{u}^k+I_{*y}^{k}d{v}^k$

为公式表达方便，令：

$I_{\&dtri;}={(I_x,I_y,I_z)}^T$

$S=I_{\&dtri;}{I_{\&dtri;}}^T---\left(47\right)$

$T=I_{\&dtri;x}{I}_{\&dtri;x}^T+I_{\&dtri;y}{I}_{\&dtri;y}^T$

经局部线性化后式(45)变为：

${\mathrm{\&Psi;}}_{D1}^{\&prime;k}\&CenterDot;(S_{11}^{k}d{u}^k+S_{12}^{k}d{v}^k+S_{13}^k)+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2}^{\&prime;k}(T_{11}^{k}d{u}^k+T_{12}^{k}d{v}^k+T_{13}^k)+\mathrm{\&beta;}\mathrm{div}({\mathrm{\&Psi;}}_S^{\&prime;k}\&dtri;(u^k+d{u}^k))=0$

${\mathrm{\&Psi;}}_{D1}^{\&prime;k}\&CenterDot;(S_{12}^{k}d{u}^k+S_{22}^{k}d{v}^k+S_{23}^k)+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2}^{\&prime;k}(T_{12}^{k}d{u}^k+T_{22}^{k}d{v}^k+T_{23}^k)+\mathrm{\&beta;}\mathrm{div}({\mathrm{\&Psi;}}_S^{\&prime;k}\&dtri;(u^k+d{v}^k))=0---\left(48\right)$

式中

${\mathrm{\&Psi;}}_{D1}^{\&prime;k}:={\mathrm{\&Psi;}}^{\&prime;}\left({(d{u}^k,{\mathrm{dv}}^k,1)}^T{S}^k(d{u}^k,{\mathrm{dv}}^k,1)\right)$

${\mathrm{\&Psi;}}_{D2}^{\&prime;k}:={\mathrm{\&Psi;}}^{\&prime;}\left({(d{u}^k,d{v}^k,1)}^T{T}^k({\mathrm{du}}^k,{\mathrm{dv}}^k,1)\right)---\left(49\right)$

${\mathrm{\&Psi;}}_S^{\&prime;k}:={\mathrm{\&Psi;}}^{\&prime;}({|\&dtri;(u^k+d{u}^k)|}^2+{|\&dtri;(v^k+d{v}^k)|}^2)$

对式(48)非线性系统进行离散化。第k次迭代时，离散化时的一些参数：k为迭代次数；图像总数目为m；像素大小

图像I(X+w^k)中，像素集合{i|1，...，N^k}，N^k为总像素点数；光流增量du^k，dv^k；在l方向上邻域N_l(i)，l∈{x，y}。经离散化后式(48)变为下式：

${\mathrm{\&Psi;}}_{D1i}^{\&prime;k}\&CenterDot;(S_{11i}^{k}d{u}_i^k+S_{12i}^{k}d{v}_i^k+S_{13i}^k)+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k}(T_{11i}^{k}d{u}_i^k+T_{12i}^{k}d{v}_i^k+T_{13i}^k)$

$+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k}}{2}\&CenterDot;\frac{u_j^k+d{u}_j^k-u_i^k-d{u}_i^k}{{\left(h_l^k\right)}^2}=0$

${\mathrm{\&Psi;}}_{D1}^{\&prime;k}\&CenterDot;(S_{12i}^{k}d{u}_i^k+S_{22i}^{k}d{v}_i^k+S_{23i}^k)+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k}(T_{12i}^{k}d{u}_i^k+T_{22i}^{k}d{v}_i^k+T_{23i}^k)---\left(50\right)$

$+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k}}{2}\&CenterDot;\frac{v_j^k+d{v}_j^k-v_i^k-d{v}_i^k}{{\left(h_l^k\right)}^2}=0$

为求解式(50)，需要经过两次迭代过程，迭代2过程嵌入在迭代1过程中。

下面结合图1，对求解过程进行详细叙述。

(1)迭代初始化，令光流分量和光流增量均为零，即(u⁰，v⁰)＝(0，0)，(du⁰，dv⁰)＝(0，0)。

(2)根据下式求解光流增量：

$\left(\begin{array}{c}{\mathrm{du}}^{k,n+1}\\ {\mathrm{dv}}^{k,n+1}\end{array}\right)={\left(\begin{array}{cc}{M}_{11}^{k,n}& M_{12}^{k,n}\\ M_{21}^{k,n}& M_{22}^{k,n}\end{array}\right)}^{-1}\left(\begin{array}{c}{\mathrm{ru}}^{k,n}\\ {\mathrm{rv}}^{k,n}\end{array}\right)---\left(51\right)$

式中，du^k，n，dv^k，n为第k次迭代1、第n次迭代2时的光流增量，

ru^k，n，rv^k，n定义如下：

$M_{11}^{k,n}={\mathrm{\&Psi;}}_{D1i}^{\&prime;k,n}{S}_{11i}^{k,n}+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k,n}{T}_{11i}^{k,n}+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2{\left(h_l^k\right)}^2}$

$M_{22}^{k,n}={\mathrm{\&Psi;}}_{D1i}^{\&prime;k,n}{S}_{22i}^{k,n}+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k,n}{T}_{22i}^{k,n}+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2{\left(h_l^k\right)}^2}---\left(52\right)$

$M_{12}^{k,n}={\mathrm{\&Psi;}}_{D1i}^{\&prime;k,n}{S}_{12i}^{k,n}+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k,n}{T}_{12i}^{k,n}$

$-{\mathrm{\&Psi;}}_{D1i}^{\&prime;k,n}{S}_{13i}^{k,n}-\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k,n}{T}_{13i}^{k,n}+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l^{-}\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2}\frac{u_j^k+d{u}_j^{k+1}-u_i^k}{{\left(h_l^k\right)}^2}$

$+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l^{+}\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2}\frac{v_j^k+d{v}_j^{k,n}-v_i^k}{{\left(h_l^k\right)}^2}=r{u}^{k,n}$ (53)

$-{\mathrm{\&Psi;}}_{D1i}^{\&prime;k,n}{S}_{23i}^{k,n}-\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k,n}{T}_{23i}^{k,n}+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l^{-}\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2}\frac{u_j^k+d{u}_j^{k+1}-u_i^k}{{\left(h_l^k\right)}^2}$

$+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l^{+}\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2}\frac{v_j^k+d{v}_j^{k,n}-v_i^k}{{\left(h_l^k\right)}^2}=r{v}^{k,n}$

式中， $N_l^{-}\left(i\right)=\{j\&Element;N_l\left(i\right)|j<i\}$ 表示已经迭代更新的像素点， $N_l^{+}\left(i\right)=\{j\&Element;N_l\left(i\right)|j>i\}$ 表示未迭代更新的点。

(3)根据相对误差判断结果是否满足要求，若相对误差小于2％则满足要求，输出光流增量；否则，将最新光流增量作为已知，继续迭代2过程。

(4)根据下式计算该像素点的光流分量：

u^k+1＝u^k+du^k，n+1    (54)

v^k+1＝v^k+dv^k，n+1

(5)根据相对误差判断光流结果是否满足要求：若相对误差小于2％则满足要求，输出此像素点光流，否则判断迭代1的次数是否已达上限(上限为m-1，即图像总数减1)；若迭代1次数未达到上限则将最新光流作为已知，继续迭代1过程；若已达上限则采取高分辨率到低分辨率策略，应用高斯低通滤波器平滑，下采样，得到低分辨率图像，使用新图像进行迭代，完成存在大偏移量时风向的提取，如图2。

(6)对得到的光流场进行直方图统计，求取频率最大值C，将频率＞0.95C的角度值矢量平均得到主风向。

第十步，风速测量，包括两种方式的测量方法，第一种为雷达图像单独测量风速，第二种为采用辅助信息测量风速。第一种，雷达单独测量风速时，仅应用航海雷达测量，我们可知的，与风速有关的信息包括：归一化的雷达散射截面(NRCS)、实测风向、、海浪的信噪比(snr)，可应用BP网络确定它们之间的关系，再根据训练好的BP网络计算风速，下面给出BP网络的设计。

BP网络结构设计：

输入层包含三个神经元：NRCS、实测风向θ和snr。一般情况下，选取航海雷达图像中的某一区域进行分析，我们认为区域内电磁波传播方向是相同的，因此在设计输入层时没有将电磁波传播方向作为输入；实测风向由基于BGCM模型海面风场反演算法计算得到；海浪信噪比snr计算公式如下：

$\mathrm{snr}=\frac{\&Integral;F^{\left(2\right)}\left(k\right)d^{2}k}{\&Integral;I_{\mathrm{bgn}}^{\left(3\right)}(k,\mathrm{\&omega;})d^2\mathrm{kd\&omega;}}---\left(55\right)$

式中，F⁽²⁾(k)为海浪谱，

为色散关系滤波器频带外部图像谱；

BP网络采用双隐层结构，隐层1、2分别有5、2个神经元，按照BP网络设计的一般原则，中间层神经元的传递函数为S型函数。每个神经元的输出可有下式得到：

$N_{\mathrm{out}}=S(-N_{\mathrm{bias}}+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_i{x}_i)---\left(56\right)$

式中，N_bias为具体到每个神经元的偏差值；n为神经元的输入总数；w_i输入权值；x_i为输入值；S为非线性转移函数，其特点是函数本身和导数都是连续的，单极型S函数定义如下：

$f\left(x\right)=\frac{1}{1+e^{-x}}---\left(57\right)$

网络的输出层只有一个神经元为风速，根据上面的分析可知，雷达单独测量风速时，BP网络结构为3×5×2×1，图中x₁、x₂、x₃分别为NRCS、归一化的实测风向、归一化的SNR，y为风速，如图8所示。

雷达单独测量时，如果仅将NRCS和实测风向作为输入，可减少输入层数量，隐层和输出层不变。

网络训练与测试：

当BP网络结构设计好后，为使网络满足实际应用要求，需要大量的输入和输出数据进行训练。在训练过程中不断调整权值，通过调整，使得BP网络的输出误差达到最小。学习速率是训练过程中重要因子，它决定每次循环权值的变化率，在一般情况下，倾向于选择较小的速率保证学习的质量。采用不同的训练函数也会影响网络性能，如收敛速度、网络推广能力等等。在考虑以上问题的基础上，BP网络的训练过程安排如下：

基于仿真图像的海面风场实验中，共模拟雷达图像350组，将它们按3∶1的比例分为训练样本和测试样本，训练次数为1000次，训练函数为trainlm，学习函数为learnlm。模拟图像过程中的输入的风速为训练阶段的参考风速。

基于实测图像的海面风场实验中，共获取雷达数据和风速计数据1577组，将它们按2∶1比例分为训练样本和测试样本，训练次数同样为1000次，训练函数为trainlm，学习函数为learnlm。风速计测量的风速为训练阶段的参考风速。

接下来对训练好的网络进行测试，在模拟图像海面风场实验中，将模拟图像过程中输入风速作为测试阶段的参考风速；在实测图像海面风场试验中，将风速计测量的风速作为测试阶段参考风速。

第二种为多信息辅助雷达测量风速，此时基于BP网络的风速测量方法如下：

由海洋物理模型函数可知，在理论上，海面风场、电磁波方向会影响NRCS的强度，除此之外海气边界层稳定与否也会影响NRCS的强度，当考虑海气边界层信息时可提高风速反演精度。表征海气边界层(MABL)的参数包括：海气温差C_(a，s)，海水盐度、大气压强P，潮位等。多信息辅助雷达测量风速时，设计的BP网络结构如下表1所示：

表1多信息辅助雷达测量时的BP网络结构

应用普通温度计测量大气温度，CDH45型海水温盐度计测量海水温度和盐度，水银气压表测量大气压，验潮站获取潮位信息。在实测实验中共获取各类数据716组，由于数据量较小，将它们全部用于网络训练。训练和学习函数、训练次数同单独雷达图像测量情况，训练阶段的参考风速由风速计测量得到。

在工程应用中可根据附加输入的种类，改变输入层的数目，隐层和输出层数目不变。

应用本发明公开的一种基于X波段航海雷达的海面风场测量技术开展实验，实验地点在国家海洋局某水文观测站附近海域。航海雷达天线架设高度为40米，平均旋转周期为2.36秒，工作在短脉冲模式，雷达探测海域半径约为2千米。2010年下半年，共获取不连续雷达图像数据1577组，获取其它多种信息716组(包括：参考风速、风向，海水盐度，大气温度，海水温度，大气压强，潮汐数据)。航海雷达测量区域中心位于船艏75度方向，近点离岸600米，面积960*960m²，此海域平均水深25米；风速计、风向标高度距离海面55m，它们的测量结果作为参考。

应用本发明方法测量的方向称为“BGCM风向”，现有技术测量的风向称为“BCM风向”，将它们分别与风向标测量风向(“参考风向”)比对，散点图和误差统计结果如图3、表2。从图中可直观看出：较传统BCM风向，BGCM风向具有更好的一致性和稳定性。表2定量的给出了风向误差统计情况，BGCM风向误差为8°远优于传统BCM风向的25.3°。

表2风向误差统计

将本发明方法测量的风速称为“实测风速”，BP网络训练阶段称为“训练风速”，测试阶段称为“测试风速”。当仅使用雷达测量时，不同阶段的实测风速与参考风速散点图如图4a和图4b，误差统计结果如表3。由结果可知，当加入海浪信噪比snr后，提高了风速测量精度，由BP网络计算的风速精度远优于光流法。

表3雷达单独测量风速误差统计

当有多种信息辅助雷达测量风速时，实测风速与参考风速散点图如图5，误差统计结果如表4。由表4可看出，已知的海气边界层参数越多，测量的风速越准确，风速误差仅为0.26m/s。

表4多信息辅助下雷达测量风速误差统计

本发明提供的海面风场测量方法，基于图像梯度和灰度不变模型进行海面风向测量方法，首先根据图像灰度和梯度不变的假设，建立图像的能量函数，并应用正则化方法求解，得到欧拉-拉格朗日非线性方程组；然后应用迭代方法求解方程组，得到光流场；最后应用直方图统计光流场的方向，获取主风向。其中基于BP神经网络的海面风速测量方法中，设计了计算风速的两种BP网络：当航海雷达单独测量时，将归一化的雷达散射截面、实测风向和海浪信噪比作为BP网络输入，计算风速；当已知海气边界层参数时，让它们作为BP网络的附加输入，计算风速。本发明提供的方法与现有技术相比，风向和风速测量精度都得到提高。

Claims (6)

1.一种基于X波段航海雷达的海面风场测量方法，其特征在于：所述的测量方法包括雷达图像的前处理、风向测量和风速测量三部分内容，所述的风向测量具体步骤如下：

令I(X)表示三维海杂波图像序列，X＝[x，y，t]^T，[x，y]是空间区域Ω内坐标，t≥0代表时间序列，w＝[u，v，1]表示t时刻图像与t+1时刻图像间光流，u、v分别为光流的x，y分量，首先给出两种假设：

假设1：在图像时间间隔内，假设图像灰度不变，那么灰度值满足下式：

$\&ForAll;X\&Element;\mathrm{\&Omega;}:{|I(X+w)-I\left(X\right)|}^2=0---\left(35\right)$

其图像能量函数如式(36)：

E_D1(w)＝∫_ΩΨ(|I(X+w)-I(X)|²)dX    (36)

式中，

ε＝10^-3为鲁棒函数，具有稳健性并抑制异常点；

假设2：当风场均匀变化或恒定时，此时雷达图像梯度不变，即：

$\&ForAll;X\&Element;\mathrm{\&Omega;}:{|\&dtri;I(X+w)-\&dtri;I\left(X\right)|}^2=0---\left(37\right)$

其图像能量函数满足下式：

$E_{D2}\left(w\right)={\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&Psi;}\left({|\&dtri;I(X+w)-\&dtri;I\left(X\right)|}^2\right)\mathrm{dX}---\left(38\right)$

式中，

为图像空间导数，

为当前图像的空间导数，

为下时刻图像的空间导数，图像空间导数由二维最优化Sobel算子求得；

将式(36)、(38)模型有机的结合起来，得到数据项能量函数，如式(39)：

E_D(w)＝E_D1(w)+αE_D2(w)    (39)

式(39)中，α调节E_D1与E_D2的比重系数；

当上述两种假设不成立时，根据式(36)、(38)及(39)构造的基于BGCM模型的能量函数如下式：

E(w)＝E_D1(w)+αE_D2(w)+βE_S(w)    (41)

式中，β为平滑项调节系数，αβ值满足下式：

$\left\{\begin{array}{c}10\&le;\mathrm{\&beta;}/\mathrm{\&alpha;}\&le;20\\ 10\&le;\mathrm{\&alpha;}\&le;20\end{array}\right.---\left(42\right)$

使用正则化方法求解光流w，光流w为下式的解：

$\left\{\begin{array}{c}\frac{\&PartialD;E\left(w\right)}{\&PartialD;u}=\frac{\&PartialD;E_{D1}\left(w\right)}{\&PartialD;u}+\mathrm{\&alpha;}\frac{\&PartialD;E_{D2}\left(w\right)}{\&PartialD;u}+\mathrm{\&beta;}\frac{\&PartialD;E_S\left(w\right)}{\&PartialD;u}=0\\ \frac{\&PartialD;E\left(w\right)}{\&PartialD;v}=\frac{\&PartialD;E_{D1}\left(w\right)}{\&PartialD;v}+\mathrm{\&alpha;}\frac{\&PartialD;E_{D2}\left(w\right)}{\&PartialD;v}+\mathrm{\&beta;}\frac{\&PartialD;E_S\left(w\right)}{\&PartialD;v}=0\end{array}\right.---\left(43\right)$

为公式表达方便，给出如下形式的简写：

$I_x={\&PartialD;}_{x}I(X+w)$

$I_y={\&PartialD;}_{y}I(X+w)$

I_z＝I(X+w)-I(X)

$I_{\mathrm{xx}}={\&PartialD;}_{\mathrm{xx}}I(X+w)$ (44)

$I_{\mathrm{xy}}={\&PartialD;}_{\mathrm{xy}}I(X+w)$

$I_{\mathrm{yy}}={\&PartialD;}_{\mathrm{yy}}I(X+w)$

$I_{\mathrm{xz}}={\&PartialD;}_{x}I(X+w)-{\&PartialD;}_{x}I\left(X\right)$

$I_{\mathrm{yz}}={\&PartialD;}_{y}I(X+w)-{\&PartialD;}_{y}I\left(X\right)$

求解式(43)得到欧拉-拉格朗日方程组如下式：

$\left\{\begin{array}{c}{\mathrm{\&Psi;}}_{D1}^{\&prime;}\left(I_z^2\right)\&CenterDot;I_z{I}_x+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2}^{\&prime;}(I_{\mathrm{xz}}^2+I_{\mathrm{yz}}^2)\&CenterDot;(I_{\mathrm{xz}}{I}_{\mathrm{xx}}+I_{\mathrm{yz}}+I_{\mathrm{xy}})+\mathrm{\&beta;}\mathrm{div}({\mathrm{\&Psi;}}_S^{\&prime;}({\left|{\&dtri;}_{3}u\right|}^2+{\left|{\&dtri;}_{3}v\right|}^2)\&CenterDot;{\&dtri;}_{3}u)=0\\ {\mathrm{\&Psi;}}_{D1}^{\&prime;}\left(I_z^2\right)\&CenterDot;I_z{I}_y+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2}^{\&prime;}(I_{\mathrm{xz}}^2+I_{\mathrm{yz}}^2)\&CenterDot;(I_{\mathrm{xz}}{I}_{\mathrm{xy}}+I_{\mathrm{yz}}{I}_{\mathrm{yy}})+\mathrm{\&beta;}\mathrm{div}({\mathrm{\&Psi;}}_S^{\&prime;}({\left|{\&dtri;}_{3}u\right|}^2+{\left|{\&dtri;}_{3}v\right|}^2)\&CenterDot;{\&dtri;}_{3}v)=0\end{array}\right.---\left(45\right)$

应用迭代法求解欧拉-拉格朗日非线性方程组，计算光流w；假设第k步迭代时，光流w^k＝[u^k，v^k，1]^T，令迭代初始条件w⁰＝[0，0，1]^T；将式(44)表示为为消除

的非线性，根据泰勒公式，用一阶近似得下式：

$I_{*z}^{k+1}=I_{*}(X+w^{k+1})-I_{*}\left(X\right)$

$\&ap;I_{*}(X+w^k)+I_{*x}^{k}d{u}^k+I_{*y}^{k}d{v}^k-I_{*}\left(X\right)---\left(46\right)$

$=I_{*z}^k+I_{*x}^{k}d{u}^k+I_{*y}^{k}d{v}^k$

为公式表达方便，令：

$I_{\&dtri;}={(I_x,I_y,I_z)}^T$

$S=I_{\&dtri;}{I_{\&dtri;}}^T---\left(47\right)$

$T=I_{\&dtri;x}{I}_{\&dtri;x}^T+I_{\&dtri;y}{I}_{\&dtri;y}^T$

经局部线性化后式(45)变为：

${\mathrm{\&Psi;}}_{D1}^{\&prime;k}\&CenterDot;(S_{11}^{k}d{u}^k+S_{12}^{k}d{v}^k+S_{13}^k)+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2}^{\&prime;k}(T_{11}^{k}d{u}^k+T_{12}^{k}d{v}^k+T_{13}^k)+\mathrm{\&beta;}\mathrm{div}({\mathrm{\&Psi;}}_S^{\&prime;k}\&dtri;(u^k+d{u}^k))=0$ (48)

${\mathrm{\&Psi;}}_{D1}^{\&prime;k}\&CenterDot;(S_{12}^{k}d{u}^k+S_{22}^{k}d{v}^k+S_{23}^k)+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2}^{\&prime;k}(T_{12}^{k}d{u}^k+T_{22}^{k}d{v}^k+T_{23}^k)+\mathrm{\&beta;}\mathrm{div}({\mathrm{\&Psi;}}_S^{\&prime;k}\&dtri;(u^k+d{v}^k))=0$

式中

${\mathrm{\&Psi;}}_{D1}^{\&prime;k}:={\mathrm{\&Psi;}}^{\&prime;}\left({(d{u}^k,{\mathrm{dv}}^k,1)}^T{S}^k(d{u}^k,{\mathrm{dv}}^k,1)\right)$

${\mathrm{\&Psi;}}_{D2}^{\&prime;k}:={\mathrm{\&Psi;}}^{\&prime;}\left({(d{u}^k,d{v}^k,1)}^T{T}^k({\mathrm{du}}^k,{\mathrm{dv}}^k,1)\right)---\left(49\right)$

${\mathrm{\&Psi;}}_S^{\&prime;k}:={\mathrm{\&Psi;}}^{\&prime;}({|\&dtri;(u^k+d{u}^k)|}^2+{|\&dtri;(v^k+d{v}^k)|}^2)$

对式(48)非线性系统进行离散化；第k次迭代时，离散化时的一些参数：k为迭代次数；图像总数目为m；像素大小

图像I(X+w^k)中，像素集合{i|1，...，N^k}，N^k为总像素点数；光流增量du^k，dv^k；在l方向上邻域N_l(i)，l∈{x，y}；经离散化后式(48)变为下式：

${\mathrm{\&Psi;}}_{D1i}^{\&prime;k}\&CenterDot;(S_{11i}^{k}d{u}_i^k+S_{12i}^{k}d{v}_i^k+S_{13i}^k)+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k}(T_{11i}^{k}d{u}_i^k+T_{12i}^{k}d{v}_i^k+T_{13i}^k)$

$+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k}}{2}\&CenterDot;\frac{u_j^k+d{u}_j^k-u_i^k-d{u}_i^k}{{\left(h_l^k\right)}^2}=0$

${\mathrm{\&Psi;}}_{D1}^{\&prime;k}\&CenterDot;(S_{12i}^{k}d{u}_i^k+S_{22i}^{k}d{v}_i^k+S_{23i}^k)+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k}(T_{12i}^{k}d{u}_i^k+T_{22i}^{k}d{v}_i^k+T_{23i}^k)---\left(50\right)$

$+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k}}{2}\&CenterDot;\frac{v_j^k+d{v}_j^k-v_i^k-d{v}_i^k}{{\left(h_l^k\right)}^2}=0$

为求解式(50)，需要经过两次迭代过程，迭代2过程嵌入在迭代1过程中，具体如下：

(1)迭代初始化，令光流分量和光流增量均为零，即(u⁰，v⁰)＝(0，0)，(du⁰，dv⁰)＝(0，0)；

(2)根据下式求解光流增量：

$\left(\begin{array}{c}{\mathrm{du}}^{k,n+1}\\ {\mathrm{dv}}^{k,n+1}\end{array}\right)={\left(\begin{array}{cc}{M}_{11}^{k,n}& M_{12}^{k,n}\\ M_{21}^{k,n}& M_{22}^{k,n}\end{array}\right)}^{-1}\left(\begin{array}{c}{\mathrm{ru}}^{k,n}\\ {\mathrm{rv}}^{k,n}\end{array}\right)---\left(51\right)$

式中，du^k，n，dv^k，n为第k次迭代1、第n次迭代2时的光流增量，

ru^k，n，rv^k，n定义如下：

$M_{11}^{k,n}={\mathrm{\&Psi;}}_{D1i}^{\&prime;k,n}{S}_{11i}^{k,n}+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k,n}{T}_{11i}^{k,n}+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2{\left(h_l^k\right)}^2}$

$M_{22}^{k,n}={\mathrm{\&Psi;}}_{D1i}^{\&prime;k,n}{S}_{22i}^{k,n}+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k,n}{T}_{22i}^{k,n}+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2{\left(h_l^k\right)}^2}---\left(52\right)$

$M_{12}^{k,n}={\mathrm{\&Psi;}}_{D1i}^{\&prime;k,n}{S}_{12i}^{k,n}+\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k,n}{T}_{12i}^{k,n}$

$-{\mathrm{\&Psi;}}_{D1i}^{\&prime;k,n}{S}_{13i}^{k,n}-\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k,n}{T}_{13i}^{k,n}+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l^{-}\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2}\frac{u_j^k+d{u}_j^{k+1}-u_i^k}{{\left(h_l^k\right)}^2}$

$+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l^{+}\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2}\frac{v_j^k+d{v}_j^{k,n}-v_i^k}{{\left(h_l^k\right)}^2}=r{u}^{k,n}$ (53)

$-{\mathrm{\&Psi;}}_{D1i}^{\&prime;k,n}{S}_{23i}^{k,n}-\mathrm{\&alpha;}{\mathrm{\&Psi;}}_{D2i}^{\&prime;k,n}{T}_{23i}^{k,n}+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l^{-}\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2}\frac{u_j^k+d{u}_j^{k+1}-u_i^k}{{\left(h_l^k\right)}^2}$

$+\mathrm{\&beta;}\underset{l=x,y}{\mathrm{\&Sigma;}}\underset{j\&Element;N_l^{+}\left(i\right)}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&Psi;}}_{\mathrm{Si}}^{\&prime;k,n}+{\mathrm{\&Psi;}}_{\mathrm{Sj}}^{\&prime;k,n}}{2}\frac{v_j^k+d{v}_j^{k,n}-v_i^k}{{\left(h_l^k\right)}^2}=r{v}^{k,n}$

式中， $N_l^{-}\left(i\right)=\{j\&Element;N_l\left(i\right)|j<i\}$ 表示已经迭代更新的像素点， $N_l^{+}\left(i\right)=\{j\&Element;N_l\left(i\right)|j>i\}$ 表示未迭代更新的点；

(3)根据相对误差判断结果是否满足要求，若相对误差小于2％则满足要求，输出光流增量；否则，将最新光流增量作为已知，继续迭代2过程；

(4)根据下式计算该像素点的光流分量：

u^k+1＝u^k+du^k，n+1    (54)

v^k+1＝v^k+dv^k，n+1

(5)根据相对误差判断光流结果是否满足要求：若相对误差小于2％则满足要求，输出此像素点光流，否则判断迭代1的次数是否已达上限；若迭代1次数未达到上限则将最新光流作为已知，继续迭代1过程；若已达上限则采取高分辨率到低分辨率策略，应用高斯低通滤波器平滑，下采样，得到低分辨率图像，使用新图像进行迭代，完成存在大偏移量时风向的提取；

(6)对得到的光流场进行直方图统计，求取频率最大值C，将频率＞0.95C的角度值矢量平均得到主风向。

所述的风速测量包括两种方式的测量方法，第一种为雷达图像单独测量风速，第二种为采用辅助信息测量风速；第一种，雷达单独测量风速时，与风速有关的信息包括归一化的雷达散射截面、实测风向、海浪的信噪比，应用BP网络确定它们之间的关系，再根据训练好的BP网络计算风速，

第二种为多信息辅助雷达测量风速，此时基于BP网络的风速测量方法如下：

表征海气边界层的参数包括：海气温差C_(a，s)，海水盐度、大气压强P，潮位，多信息辅助雷达测量风速时，设计的BP网络结构如下所示：

网络结构：          参数：

输入层：            雷达单独测量输入：NRCS，SNR，θ；

                    附加输入：C_(a，s)，盐度，P，潮位；

隐层：              隐层1：8神经元；

                    隐层2：5神经元；

                    隐层3：3神经元；

输出层：            风速。

2.根据权利要求1所述的一种基于X波段航海雷达的海面风场测量方法，其特征在于：所述的雷达图像的前处理通过如下步骤实现：

第一步，雷达图像方位和角度校正，以及应用中值滤波器抑制噪声；

所述的雷达图像方位和角度校正就是要消除方位和角度对回波强度的影响；所述的中值滤波器抑制噪声就是应用相邻像素的灰度中值来替代该像素的灰度值：

$I(x,y)=\underset{(s,t)\&Element;S_{\mathrm{xy}}}{\mathrm{median}}\left\{g(s,t)\right\}---\left(3\right)$

式(3)中，s_xy表示中心在(x，y)点，尺寸为m×n的矩形窗口的坐标组，g(s，t)为雷达图像，I(x，y)为中值滤波后图像的灰度值；

第二步，在第一步处理校正和滤波后的雷达图像中选取测量区域；

选取测量区域的原则是：应用最近点插值方法将选取的测量区域图像由极坐标转换为直角坐标，多幅子图像组成子图像序列σ(x，y，t)，总时间长度在40～100s范围；

第三步，确定三维波数频率图像谱：对子图像序列σ(x，y，t)做三维离散傅里叶变换如下式所示，变换后得到三维波数频率图像谱：

$I^{\left(3\right)}(k_x,k_y,\mathrm{\&omega;})=\underset{y=0}{\overset{L_y}{\mathrm{\&Sigma;}}}\underset{x=0}{\overset{L_x}{\mathrm{\&Sigma;}}}\underset{t=0}{\overset{T}{\mathrm{\&Sigma;}}}\mathrm{\&sigma;}(x,y,t)\exp [-2\mathrm{\&pi;i}(k_{x}x/L_x+k_{y}y/L_y+\mathrm{\&omega;t}/T)]---\left(4\right)$

式(4)中，L_x、L_y为子图像空间大小，T为总时间长度；k为波数，k_x、k_y分别为波数k在x和y方向的波数分量；

图像功率谱为频谱的平方，通过下式求取：

P(k_x，k_y，ω)＝|I⁽³⁾(k_x，k_y，ω)|²＝R_e ²(k_x，k_y，ω)+I_m ²(k_x，k_y，ω)    (5)

式(5)中R_e、I_m分别代表I⁽³⁾(k_x，k_y，ω)的实部和虚部；

图像谱分辨率受图像时间和空间尺度限制，即：

$\mathrm{\&Delta;}{k}_x=\frac{2}{L_x},$ $\mathrm{\&Delta;}{k}_y=\frac{2\mathrm{\&pi;}}{L_y},$ $\mathrm{\&Delta;\&omega;}=\frac{2\mathrm{\&pi;}}{T}---\left(6\right)$

式(6)中，Δk_x与Δk_y为波数分量分辨率；Δω为频率分辨率；

第四步，计算海流：

4.1、对于第三步中的三维波数频率图像谱，应用如下形式的色散关系带通滤波器抑制噪声干扰：

$\left|k_p\right|=\frac{{(\mathrm{\&omega;}+\frac{\mathrm{\&Delta;\&omega;}}{2}+|k_m|\&CenterDot;|U_{\max }\left|\right)}^2}{g}+\frac{\mathrm{\&Delta;k}}{\sqrt{2}}---\left(8\right)$

$\left|k_n\right|=\frac{{(\mathrm{\&omega;}-\frac{\mathrm{\&Delta;\&omega;}}{2}-|k_m|\&CenterDot;|U_{\max }\left|\right)}^2}{g}-\frac{\mathrm{\&Delta;k}}{\sqrt{2}}---\left(9\right)$

上式中，|k_n|、|k_p|为带通滤波器频带边界；Δω为频率分辨率，ω为海浪频率，Δk为波数分辨率；U_max为雷达天线与海浪场间的最大相对速度，当雷达载体静止时，U_max为最大海表面流速，U_max取值在1.2～2m/s；F⁽³⁾(k，ω)为滤波后海浪图像谱，k_m为主导波数；

4.2、自适应选取阈值构造低通滤波器，抑制能量低于阈值的噪声；

4.3、构造目标函数；

4.4、海流计算：

(1)、初始估流：在此仅考虑0阶次波，根据0阶色散关系方程计算海流，步骤如下：

①取阈值C_fg，使得能量大于阈值C_fg的谱数量n₂，选择阈值C_fg的量值为0.2或n₂值在数十个；

②将目标函数Q²对海流U的两个分量U_x、U_y求偏导，并使其为零，得到海流的计算公式如下：

U＝D·B    (23)

$D={\left[\begin{array}{cc}\mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x^2& \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x{k}_y\\ \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x{k}_y& \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_y^2\end{array}\right]}^{-1}---\left(24\right)$

$B=\left[\begin{array}{c}\mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x({\mathrm{\&omega;}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k|d})\\ \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_y({\mathrm{\&omega;}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k|d})\end{array}\right]---\left(25\right)$

式中，F_M为海浪谱；μ为隶属度；k_x、k_y为k_ij的分量；ω₀为0阶波频率，d为海深，其它变量的定义同前。经过初始估流阶段可得到粗略精度的海流；

(2)、迭代估流：根据0、1阶次波确定海流：

①应用实时色散关系带通滤波器计算信噪比snr；然后选取大于阈值C_it的图像谱，数目为n₃；再根据初始估流得到的较低精度的海流，应用色散关系方程，计算0、1阶次波频率：

${\mathrm{\&omega;}}_p=(p+1)\sqrt{\frac{g\left|k\right|}{p+1}\tanh \left(\frac{\left|k\right|d}{p+1}\right)}+\left|k\right|\left|u\right|\cos \left(\mathrm{\&theta;}\right)---\left(26\right)$

式中，p为海浪能量的阶次，p＝0时，该式就变成了基本色散关系方程，p＞1时的海浪被称作高阶次波，ω_p为p阶次波的频率；②对这n₃个图像谱进行如下讨论：若图像谱均在[0，ω_N]范围内，再判断它们属于0阶次还是1阶次波，若|ω_i-ω₀(k_i)|＜|ω_i-ω₁(k_i)|则属于0阶次波，若|ω_i-ω₀(k_i)|＞|ω_i-ω₁(k_i)|则属于1阶次波；若有图像谱超出[0，ω_N]频率范围的，需要按照傅里叶变换的性质对信号进行重构，然后再判断它属于0阶还是1阶次波；其中ω_N为奈奎斯特频率；

将判断好的数据根据下式计算海流：

U＝D·B    (27)

$D={\left[\begin{array}{cc}\mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x^2& \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x{k}_y\\ \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_x{k}_y& \mathrm{\&Sigma;}{F}_M\mathrm{\&mu;}{k}_y^2\end{array}\right]}^{-1}---\left(28\right)$

$B=\left[\begin{array}{c}\underset{n_{30}}{\mathrm{\&Sigma;}}{F}_M\mathrm{\&mu;}{k}_x({\mathrm{\&omega;}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k\left|\right)d})+\underset{n_{31}}{\mathrm{\&Sigma;}}{F}_M\mathrm{\&mu;}{k}_x({\mathrm{\&omega;}}_1-2\sqrt{\frac{g\left|k\right|}{2}\tanh (\frac{\left|k\right|d}{2}})\\ \underset{n_{30}}{\mathrm{\&Sigma;}}{F}_M\mathrm{\&mu;}{k}_y({\mathrm{\&omega;}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k|d})+\underset{n_{31}}{\mathrm{\&Sigma;}}{F}_M\mathrm{\&mu;}{k}_y({\mathrm{\&omega;}}_1-2\sqrt{\frac{g\left|k\right|}{2}\tanh (\frac{\left|k\right|d}{2}})\end{array}\right]---\left(29\right)$

式中，n₃₀、n₃₁分别为n₃个谱分量中0、1阶次波数目，ω₀、ω₁分别为0、1阶次波频率，其它变量定义同前；

③迭代10至15次满足需要，判断是否迭代结束，如果没有结束则将新表面流代入式(26)，得到新的0阶次和1阶次波频率，返回步骤①；若迭代结束则停止计算，输出海流y；

第五步，在确定实时海流后，针对应用实时色散关系带通滤波器抑制海浪噪声对海面风场的干扰，滤波器形式如下：

式中，k_mn为实测波数；F⁽³⁾(k_mn，ω_v)为图像谱，P_D ⁽³⁾(k_mn，ω_v)为滤波后图像谱；

第六步，对滤波后的海浪图像谱P_D ⁽³⁾(k_mn，ω_v)进行三维反傅里叶变换，得到滤波后图像序列σ′(x，y，t)，如下式：

${\mathrm{\&sigma;}}^{\&prime;}(x,y,t)=\frac{1}{L_x{L}_{y}T}\underset{y=0}{\overset{L_y}{\mathrm{\&Sigma;}}}\underset{x=0}{\overset{L_x}{\mathrm{\&Sigma;}}}\underset{t=0}{\overset{T}{\mathrm{\&Sigma;}}}{P}_D(k_x,k_y,\mathrm{\&omega;})\exp \left[2\mathrm{\&pi;i}(k_{x}x/L_x+k_{y}y/L_y+\mathrm{\&omega;t}/T)\right]---\left(32\right)$

式中，P_D(k_x，k_y，ω)为滤波后图像谱；

第七步，对滤波后图像序列σ′(x，y，t)进行滑动平均，抑制在时间上变化频率较高的噪声，保留静态信号和周期较大的信号，形式如下：

$G(x,y,,\mathrm{\&tau;})=\frac{1}{n}\underset{t=i}{\overset{i+n-1}{\mathrm{\&Sigma;}}}{\mathrm{\&sigma;}}^{\&prime;}(x,y,t)---\left(33\right)$

式中，G(x，y，τ)为滑动平均后新序列，新序列图像数量为N_t-n+1，N_t为图像数目，n为参与平均的图像数量，τ为滑动平均后新序列的图像的时间。

第八步，从图像空间角度采取措施抑制噪声干扰：应用低通滤波器进行图像的缩减和平滑，抑制空间尺度小的噪声，低通滤波器设计步骤如下：

(1)用4阶二项式系数卷积核B⁴去平滑图像；

(2)对子图像取2×2的抽样平均；

(3)用2阶二项式系数卷积核B²去平滑图像。

3.根据权利要求2所述的一种基于X波段航海雷达的海面风场测量方法，其特征在于：所述的主导波数根据下面方法确定：

将滤波后海浪图像谱F⁽³⁾(k，ω)按波数角度进行积分，得到二维波数模频率域海浪图像谱I⁽²⁾(|k|，ω)，如下式：

$I^{\left(2\right)}\left(\right|k|,\mathrm{\&omega;})={\&Integral;}_0^{2\mathrm{\&pi;}}{F}^{\left(3\right)}\left(\right|k|,{\mathrm{\&theta;}}_1,\mathrm{\&omega;})d{\mathrm{\&theta;}}_1---\left(10\right)$

式(10)中，F⁽³⁾(|k|，θ₁，ω)＝F⁽³⁾(k，ω)为三维波数频率图像谱；θ₁为海浪谱波数角度。在I⁽²⁾(|k|，ω)中，C为最大能量值，选取能量大于0.95*C的能量点，将它们的波数模进行算术平均，即得到主导波数|k_m|。

4.根据权利要求2所述的一种基于X波段航海雷达的海面风场测量方法，其特征在于：所述的步骤4.2具体通过如下步骤实现：

首先，根据滤波后海浪图像谱F⁽³⁾(k，ω)和信噪比snr计算海浪谱总能量P，如下式：

$P=\frac{\mathrm{snr}+1}{\mathrm{snr}}\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}\underset{k}{\overset{N_{\mathrm{\&omega;}}}{\mathrm{\&Sigma;}}}{F}^{\left(3\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}},{\mathrm{\&omega;}}_k)\&CenterDot;\mathrm{\&Delta;}{k}_x\mathrm{\&Delta;}{k}_y\mathrm{\&Delta;\&omega;}---\left(11\right)$

式中：

$\mathrm{snr}=\frac{\mathrm{sig}}{\mathrm{bgn}}---\left(12\right)$

$\mathrm{sig}=\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}{F}^{\left(2\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}})\mathrm{\&Delta;}{k}_x\mathrm{\&Delta;}{k}_y---\left(13\right)$

$\mathrm{bgn}=\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}\underset{l}{\overset{N_{\mathrm{\&omega;}}}{\mathrm{\&Sigma;}}}{I}^{\left(3\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}},{\mathrm{\&omega;}}_l)\mathrm{\&Delta;}{k}_x\mathrm{\&Delta;}{k}_y\mathrm{\&Delta;\&omega;}-\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}{F}^{\left(2\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}})\mathrm{\&Delta;}{k}_x\mathrm{\&Delta;}{k}_y---\left(14\right)$

$F^{\left(2\right)}\left(k\right)=2{\&Integral;}_{\mathrm{\&omega;}>0}{I}^{\left(3\right)}(k_x,k_y,\mathrm{\&omega;})\mathrm{\&delta;}(\mathrm{\&omega;}-{\stackrel{\&OverBar;}{\mathrm{\&omega;}}}_0)\mathrm{d\&omega;}---\left(15\right)$

式中，sig为海浪图像谱能量，bgn为背景噪声能量，I⁽³⁾(k，ω)是三维图像谱，F⁽²⁾(k)是二维海浪图像谱，N_kx、N_ky、N_ω为三维波数频率图像谱的波数在x方向和y方向的范围，以及频率的范围，Δk_x、Δk_y为波数分辨率，Δω为频率分辨率，

是色散关系带通滤波器；

然后，对三维波数频率图像谱中的所有的能量点按照从大到小的顺序进行排序，排序后第i个能量点的能量为P_i，前m个能量点的总和为P，即：

$P=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{P}_i---\left(16\right)$

最后，由上式求出数目m，第m个能量点的能量值为P_m，那么设P_m为阈值C_it；

根据上述的阈值C_it构造低通滤波器，选取能量大于阈值C_it的能量点，得到低通滤波后的三维波数频率图像谱。

5.根据权利要求2所述的一种基于X波段航海雷达的海面风场测量方法，其特征在于：所述的步骤4.3构造目标函数的具体步骤为：

(1)、计算海浪谱

应用调制传递函数MTF对三维波数频率图像谱进行非线性校正，得到海浪谱

MTF定义如下式：

MTF＝|k|^β    (17)

MTF校正方法如下式：

$F_M^{\left(3\right)}(k,\mathrm{\&omega;})={\left|k\right|}^{\mathrm{\&beta;}}\&CenterDot;I^{\left(3\right)}(k,\mathrm{\&omega;})---\left(18\right)$

式中，I⁽³⁾(k，ω)为图像谱，图像谱经调制传递函数校正后变为海浪谱

β值范围大于等于1.17小于等于1.21；

(2)、计算隶属度：

对于固定频率ω_i上有N个谱分量组成集合{k₁k₂...k_N}，那么中心k_i0，最大半径r_i定义如下：

$k_{i0}=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{k}_j---\left(19\right)$

$r_i=\underset{k_j}{\max }\left|\right|k_j-k_{i0}\left|\right|---\left(20\right)$

式中，k_j为谱分量；j＝1，2，…，N；i＝1，2…，N_ω，根据所述的谱分量到中心的距离确定隶属度时，集合中谱分量隶属度为模糊隶属度μ(k_ij)：

$\mathrm{\&mu;}\left(k_{\mathrm{ij}}\right)=1-\frac{\left|\right|k_j-k_{i0}\left|\right|}{r_i+\mathrm{\&delta;}}---\left(21\right)$

式中，δ为很小的正数，目的是避免出现μ(k_ij)＝0的情况；

(3)、构造目标函数：

将海浪谱

和模糊隶属度μ(k_ij)的乘积作为最小二乘法的权值，那么海流反演计算的目标函数为：

式中，ω_i为实测数据海浪谱

中第i个频率分量，ω_p为理论上色散关系方程给出的频率分量；N_ω为三维波数频率图像谱的频率的范围，N为ω_i频率上具有的谱分量数目；

μ(k_ij)分别为ω_i频率上第j个海浪谱和隶属度。

6.根据权利要求1所述的一种基于X波段航海雷达的海面风场测量方法，其特征在于：BP网络结构设计如下：

输入层包含三个神经元：雷达散射截面NRCS、实测风向θ和海浪的信噪比snr，其中海浪信噪比snr计算公式如下：

$\mathrm{snr}=\frac{\&Integral;F^{\left(2\right)}\left(k\right)d^{2}k}{\&Integral;I_{\mathrm{bgn}}^{\left(3\right)}(k,\mathrm{\&omega;})d^2\mathrm{kd\&omega;}}---\left(55\right)$

式中，F⁽²⁾(k)为海浪谱，

为色散关系滤波器频带外部图像谱；

BP网络采用双隐层结构，隐层1、2分别有5、2个神经元，按照BP网络设计的一般原则，中间层神经元的传递函数为S型函数，每个神经元的输出由下式得到：

$N_{\mathrm{out}}=S(-N_{\mathrm{bias}}+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_i{x}_i)---\left(56\right)$

式中，N_bias为具体到每个神经元的偏差值；n为神经元的输入总数；w_i输入权值；x_i为输入值；S为非线性转移函数，单极型S函数定义如下：

$f\left(x\right)=\frac{1}{1+e^{-x}}---\left(57\right)$

网络的输出层只有一个神经元为风速，雷达单独测量风速时，BP网络结构为3×5×2×1。

Images