[发明专利]一种基于X波段航海雷达的海面风场测量方法

权利要求书

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：在图像时间间隔内，假设图像灰度不变，那么灰度值满足下式：

∀X∈Ω:|I(X+w)-I(X)|2=0---(35)]]>

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

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

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

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

∀X∈Ω:|▿I(X+w)-▿I(X)|2=0---(37)]]>

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

ED2(w)=∫ΩΨ(|▿I(X+w)-▿I(X)|2)dX---(38)]]>

式中，为图像空间导数，为当前图像的空间导数，为下时刻图像的空间导数，图像空间导数由二维最优化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)

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

10≤β/α≤2010≤α≤20---(42)]]>

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

∂E(w)∂u=∂ED1(w)∂u+α∂ED2(w)∂u+β∂ES(w)∂u=0∂E(w)∂v=∂ED1(w)∂v+α∂ED2(w)∂v+β∂ES(w)∂v=0---(43)]]>

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

Ix=∂xI(X+w)]]>

Iy=∂yI(X+w)]]>

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

Ixx=∂xxI(X+w)]]>(44)

Ixy=∂xyI(X+w)]]>

Iyy=∂yyI(X+w)]]>

Ixz=∂xI(X+w)-∂xI(X)]]>

Iyz=∂yI(X+w)-∂yI(X)]]>

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

ΨD1′(Iz2)·IzIx+αΨD2′(Ixz2+Iyz2)·(IxzIxx+Iyz+Ixy)+βdiv(ΨS′(|▿3u|2+|▿3v|2)·▿3u)=0ΨD1′(Iz2)·IzIy+αΨD2′(Ixz2+Iyz2)·(IxzIxy+IyzIyy)+βdiv(ΨS′(|▿3u|2+|▿3v|2)·▿3v)=0---(45)]]>

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

I*zk+1=I*(X+wk+1)-I*(X)]]>

≈I*(X+wk)+I*xkduk+I*ykdvk-I*(X)---(46)]]>

=I*zk+I*xkduk+I*ykdvk]]>

为公式表达方便，令：

I▿=(Ix,Iy,Iz)T]]>

S=I▿I▿T---(47)]]>

T=I▿xI▿xT+I▿yI▿yT]]>

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

ΨD1′k·(S11kduk+S12kdvk+S13k)+αΨD2′k(T11kduk+T12kdvk+T13k)+βdiv(ΨS′k▿(uk+duk))=0]]>(48)

ΨD1′k·(S12kduk+S22kdvk+S23k)+αΨD2′k(T12kduk+T22kdvk+T23k)+βdiv(ΨS′k▿(uk+dvk))=0]]>

式中

ΨD1′k:=Ψ′((duk,dvk,1)TSk(duk,dvk,1))]]>

ΨD2′k:=Ψ′((duk,dvk,1)TTk(duk,dvk,1))---(49)]]>

ΨS′k:=Ψ′(|▿(uk+duk)|2+|▿(vk+dvk)|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)变为下式：

ΨD1i′k·(S11ikduik+S12ikdvik+S13ik)+αΨD2i′k(T11ikduik+T12ikdvik+T13ik)]]>

+βΣl=x,yΣj∈Nl(i)ΨSi′k+ΨSj′k2·ujk+dujk-uik-duik(hlk)2=0]]>

ΨD1′k·(S12ikduik+S22ikdvik+S23ik)+αΨD2i′k(T12ikduik+T22ikdvik+T23ik)---(50)]]>

+βΣl=x,yΣj∈Nl(i)ΨSi′k+ΨSj′k2·vjk+dvjk-vik-dvik(hlk)2=0]]>

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

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

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

duk,n+1dvk,n+1=M11k,nM12k,nM21k,nM22k,n-1ruk,nrvk,n---(51)]]>

式中，du^k，n，dv^k，n为第k次迭代1、第n次迭代2时的光流增量，ru^k，n，rv^k，n定义如下：

M11k,n=ΨD1i′k,nS11ik,n+αΨD2i′k,nT11ik,n+βΣl=x,yΣj∈Nl(i)ΨSi′k,n+ΨSj′k,n2(hlk)2]]>

M22k,n=ΨD1i′k,nS22ik,n+αΨD2i′k,nT22ik,n+βΣl=x,yΣj∈Nl(i)ΨSi′k,n+ΨSj′k,n2(hlk)2---(52)]]>

M12k,n=ΨD1i′k,nS12ik,n+αΨD2i′k,nT12ik,n]]>

-ΨD1i′k,nS13ik,n-αΨD2i′k,nT13ik,n+βΣl=x,yΣj∈Nl-(i)ΨSi′k,n+ΨSj′k,n2ujk+dujk+1-uik(hlk)2]]>

+βΣl=x,yΣj∈Nl+(i)ΨSi′k,n+ΨSj′k,n2vjk+dvjk,n-vik(hlk)2=ruk,n]]>(53)

-ΨD1i′k,nS23ik,n-αΨD2i′k,nT23ik,n+βΣl=x,yΣj∈Nl-(i)ΨSi′k,n+ΨSj′k,n2ujk+dujk+1-uik(hlk)2]]>

+βΣl=x,yΣj∈Nl+(i)ΨSi′k,n+ΨSj′k,n2vjk+dvjk,n-vik(hlk)2=rvk,n]]>

式中，Nl-(i)={j∈Nl(i)|j<i}]]>表示已经迭代更新的像素点，Nl+(i)={j∈Nl(i)|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神经元；

输出层：            风速。