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

摘要

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

主权项

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神经元；输出层：            风速。