专利名称：一种基于x波段航海雷达的海面风场测量方法
技术领域：
本发明属于海洋动力环境遥感技术领域，特别涉及应用X波段航海雷达测量海面风场参数的方法。
背景技术：
海面风场是海洋表面活动的主要动力来源，调节着海洋上层和低层大气之间的能量、动量、物质和水气的交换，直接影响着人类海上活动(海洋渔业、海上航行、海洋工程等)，它是海/气边界层以及生物圈系统的关键因素。风掠过海面时，引起海面的粗糙不平，即微尺度波，它们与航海雷达的电磁波波长类似，当电磁波遇到微尺度波时，发生Bragg共振，产生后向回波，形成海杂波图像。2004年德国GKSS研究中心Dankert等人首次应用基于图像灰度不变模型(BCM)的光流法从海 杂波图像序列中提取了海面风场，此方法简单，不需要校正雷达相位。见参考文献[I](王剑，段华敏.X波段雷达图像提取海洋表面风场.海洋技术，2010，29 (3) :5-8页)和参考文献[2] (Dankert H, Horstmann J, Rosenthal ff. Ocean surface winds retrieved frommarine radar image sequences. International Geosciences and Remote SensingSymposium, 2004，3 1903-1906P)，现有的BCM模型海面风场反演技术包含以下几个步骤步骤I、建立光流方程。假设短时间内图像灰度值不变，即满足下式G (X+f, t+ A t) = G(X, t)式中，G(X,t)为图像灰度值；X= [x,y]T为图像空间坐标；f = [u,v]T为光流；U、v为光流分量；At为两幅图的时间间隔。上述的方程为非线性方程，为简化运算，常常将上述的方程按照一阶泰勒公式展开，令At — 0，取极限并略去高阶项可得到基本约束方程，如下式Gxu+Gyv+Gt = 0式中，G1 = dG/dx ,Gy = SG/+为灰度空间变化率，G1 = dG/dt为灰度时间变化率。上式可简写为Gf+g = 0式中，G = (Gx,Gy)为灰度值空间导数；f = (u,v)为光流；g = Gt为灰度时间变化率；该式即为光流方程。光流方程表明海杂波图像中某一点的灰度时间变化率的大小是灰度空间导数与该点空间运动速度的乘积。步骤2、应用加权最小二乘法求解光流方程，求取准则是J (f) = (Gf+g) tW (Gf+g) = min式中，W是适当取值的正定加权矩阵。要使上述的求取准则方程成立，f应满足^ = GTW(Gf + g) + [(Gf + gfWGf = 0从中解得
f = [Gt (ff+ffT) G] ^1Gt (ff+ffT) g一般情况下，正定加权矩阵W取对称阵，即W = Wt,所以加权最小二乘估计为f = - (GtWG) ^1G1Wg为提高计算精度，取正定加权矩阵W为Sobel算子B。G1BGf = -G1Bg将上式展开
'SG][SG'
D dx「aG ^GITm] dx SGB— — =~B Pir ~T
dG dx Dy vdG dt
DyDy
\dG dG dG DGT1 [dG dGl
D---D---B---
「 I Udx dx dx dydx dt
L 」v dG dG dG dG ^dG dGB---B----D----—
dy dx 办办」L办沒_ 记为如下形式
r n 「m"|BG-GvV [BGx-GiI
「00291=- 、xJ x y x 1
L J [vj [BGy-Gx B(Gyf \ [BGy-Gt_ 假设在小邻域U内光流是常量，U内的像素点构成集合为(G1, G2, ...，GJ，运用加权最小—乘估计，此时有
'DG1 BG2 dGn~
厂r 办办dx =么么，Sgl
_ Sy Sy &」~dT ... ~dT那么光流分量的计算公式为
r n [ YjBiGic f YjBGic -G1J [YjBGic. GitM = - :=1':1
卜」YjBG6l-Giy YjB(Giy )2 ^BGiy-Gtt _ /=1/=1」U=I_式中
O*步骤3、应用直方图统计得到的光流场，最大值即为主风向或主风速。上述的反演方法中，由于海杂波图像序列具有较强的非线性，应用BCM模型计算海面风场会存在较大误差。

发明内容
为提高应用X波段航海雷达测量海面风场的测量精度，本发明公布了基于图像梯度、灰度不变模型(BGCM)能量函数的海面风场测量方法。所述的海面风场测量方法包含三部分雷达图像前处理、风向测量和风速测量，所述的风向测量具体步骤如下令I (X)表示三维海杂波图像序列，X = [X, Y，t]T，[x, y]是空间区域Q内坐标，t彡0代表时间序列，w = [u, V, I]表示t时刻图像与t+1时刻图像间光流，U、V分别为光流的x，y分量，首先给出两种假设假设I :在图像时间间隔内，假设图像灰度不变，那么灰度值满足下式
VXe Q: \I{X + w)-I{Xf = 0(35)其图像能量函数如式(36)Em (w) = f 叫(11 (X+w) -I ⑴ 12) dX(36)式中，￥(s2)= ^s2 +s2，e = 10 3为鲁棒函数,具有稳健性并抑制异常点；假设2 :当风场均匀变化或恒定时，此时雷达图像梯度不变，即VXeQ-.\VI(X + w)- V/(X)|2 = 0(37)其图像能量函数满足下式ED2(w) = [/(|W(X + ^)- VI(Xf)dX(3 8)
式中，V I= (Ix, Iy)T为图像空间导数，W(JQ为当前图像的空间导数，W(X + *v)为下时刻图像的空间导数，图像空间导数由二维最优化Sobel算子求得；将式(36)、(38)模型有机的结合起来，得到数据项能量函数，如式(39)Ed(w) = Edi (w) + a Ed2 (w) (39)式(39)中，a调节Em与Ed2的比重系数；当上述两种假设不成立时，根据式(36)、(38)及(39)构造的基于BGCM模型的能量函数如下式E(w) = Edi (w) + a Ed2 (W) + 3 Es (w) (41)式中，0为平滑项调节系数，a、^值满足下式
_] j 10, ,20 (42) 使用正则化方法求解光流W，光流w为下式的解
'SE(W) _ dEm(w) | JEm(W) | ^dEs(W) _Q I dudududu
dE(w) dEm(w) | dEm(w) | dEs(w) Q(43)为公式表达方便，给出如下形式的简写Ix=dxI{X + w)Iy=dyI{X+w)Iz = I (X+w) -I (X)Ixx=DJiX+ w)ixy=8xyI(X+ w)(44)Iyy=dyyI{X+w)Ixz=dxI{X + w)-dxI{X)Iyz=SyIiX+w)-dyI{X)求解式(43)得到欧拉-拉格朗日方程组如下式1K此 + O. Q1 丄 +IyJxy) +f)div(W^\W,u\2 +|V3v|2).V3M) = 0(的)
I1K1CO'"y+0.(4人 +Iyj ) + f)div(^(\v4 +|v3v|2).v3v) = o应用迭代法求解欧拉-拉格朗日非线性方程组，计算光流W ;假设第k步迭代时，光流Wk = [uk,vk, 1]T，令迭代初始条件《° =
T;将式(44)表示为U ;为消除f的非线性，根据泰勒公式，用一阶近似得下式Il*1 =I,{X +wk+l)~ I,{X)
I,{X +wk) + Ilxduk + Ik,ydvk -L(X) (46)= Itz + Iiduk + llydvh为公式表达方便，令7V =(IxJyJz)TS=IvI/(47)T = IvJX经局部线性化后式(45)变为■ (S^duk +S^dvk +S^ + aW'^iT^du +T^dvk +T11I) + ^div{W'sV{u + du)) = 0^ . (S^duh + Sk22dvk + Skn) + a W* (T1lIduk + T^dvk + T2\) + pdivQFpiv11 +dvk)) = 0(48)式中d := lFXidulc, dv\l)T SlcIdulcJvlc,!))%k2 := lF'{{duk,dvk,\)TTk{duk,dvk,\))(49)lFf :=f (|v<y +duk)\2 +|V(v4 +dvk)[)·
对式(48)非线性系统进行离散化；第k次迭代时，离散化时的一些参数k为迭代次数;图像总数目为m ;像素大小 < ;图像I (X+wk)中，像素集合{illNk为总像素点数；光流增量duk，dvk;在I方向上邻域队(1)，1 G {x，y};经离散化后式(48)变为下式Ki .Of +S^ + aniiTuA- +T1IM +T1D
r n^ ^ lPsi +lPsi Uki +du) - U^ -du^ ,、+Al 2 , j -~TTTTTL = 0 (50)
I=X^jeNl(I)LKnI )K-iS^M +Sk11M +S^ + ani^du- + W +TL)
^ ^ lFlk +1F'11 vk +dvk -vk -dvk+PY. 2 7 -~77T-2L=0
I=X^yjeNlH) 1V1I )为求解式(50)，需要经过两次迭代过程，迭代2过程嵌入在迭代I过程中，具体如下(I)迭代初始化，令光流分量和光流增量均为零，即(u°，V0) = (0,0), (du0, dv°)=(0,0)；(2)根据下式求解光流增量Knt11I = K'" Mkf] 1KI(51)
{dvk-"+1) IvM2*;" Mk22") {rvk-")v ,式中，duk’ n, dvk’n为第k次迭代I、第n次迭代2时的光流增量，Mku"MgMKruk’n，rvk’nS义如下MklC =+aW^f +(SJ^ X
I=X^jeNl(X) l^hI )Mk22" =+ <2力;+d H(52)
I=X^jeNl(X) l^hI )Mf2" = K;"S=
,…，… ^ ^ W'f-n+W'f-n Uki+du11；1-^-^S^-aW^+fS^ X 9 1] (hL
l=x,y jeN~(02^nl )
权利要求
1. 一种基于X波段航海雷达的海面风场测量方法，其特征在于所述的测量方法包括雷达图像的前处理、风向測量和风速测量三部分内容，所述的风向測量具体步骤如下 令I (X)表示三维海杂波图像序列，X = [x，y，t]T, [x,y]是空间区域Ω内坐标，t彡O代表时间序列，w = [u, V, I]表示t时刻图像与t+Ι时刻图像间光流,u、v分别为光流的X,y分量，首先给出两种假设 假设I :在图像时间间隔内，假设图像灰度不变，那么灰度值满足下式 /ΧεΩ·. \I{X + w) _ 7(X)|2 = 0 (3 5) 其图像能量函数如式(36)Em (w) = / Ω Ψ (11 (X+w) -I ⑴ 12) dX (36) 式中，^^ = #77，ε = 10_3为鲁棒函数，具有稳健性并抑制异常点； 假设2 :当风场均匀变化或恒定时，此时雷达图像梯度不变，即 VX e Ω: \VI(X + w) - W(X)|2 =O (37) 其图像能量函数满足下式 Ed2(W) = \ψφ {Χ+w)- VI{X)\2W(3 8 式中，▽/ = (/ノブ为图像空间导数，为当前图像的空间导数，W(X +め为下时刻图像的空间导数，图像空间导数由ニ维最优化Sobel算子求得； 将式(36)、(38)模型有机的结合起来，得到数据项能量函数，如式(39)Ed (w) = Edi (w) + a Ed2 (w) (39) 式(39)中，α调节Em与Ed2的比重系数； 当上述两种假设不成立时，根据式(36)、(38)及(39)构造的基于BGCM模型的能量函数如下式E (w) = Edi (w) + a Ed2 (W) + β Es (w) (41) 式中，β为平滑项调节系数，α β值满足下式 |10<^/ <20 (42) I 10<α<20 、 ' 使用正则化方法求解光流W，光流W为下式的解
2.根据权利要求I所述的一种基于X波段航海雷达的海面风场测量方法，其特征在于所述的雷达图像的前处理通过如下步骤实现 第一步，雷达图像方位和角度校正，以及应用中值滤波器抑制噪声； 所述的雷达图像方位和角度校正就是要消除方位和角度对回波强度的影响；所述的中值滤波器抑制噪声就是应用相邻像素的灰度中值来替代该像素的灰度值
3.根据权利要求2所述的一种基于X波段航海雷达的海面风场测量方法，其特征在于所述的主导波数根据下面方法确定 将滤波后海浪图像谱F(3)(k，co)按波数角度进行积分，得到二维波数模频率域海浪图像谱I(2)(|k|，《)，如下式
4.根据权利要求2所述的一种基于X波段航海雷达的海面风场测量方法，其特征在于所述的步骤4. 2具体通过如下步骤实现 首先，根据滤波后海浪图像谱F(3)(k，co)和信噪比snr计算海浪谱总能量P，如下式
5.根据权利要求2所述的一种基于X波段航海雷达的海面风场测量方法，其特征在于所述的步骤4. 3构造目标函数的具体步骤为 (1)、计算海浪谱if#, )； 应用调制传递函数MTF对三维波数频率图像谱进行非线性校正，得到海浪谱, MTF定义如下式 MTF = |k|e (17) MTF校正方法如下式
6.根据权利要求I所述的一种基于X波段航海雷达的海面风场测量方法，其特征在于BP网络结构设计如下 输入层包含三个神经元雷达散射截面NRCS、实测风向0和海浪的信噪比snr，其中海浪信噪比snr计算公式如下
全文摘要
本发明公开了一种基于X波段航海雷达的海面风场测量方法，属于海洋动力环境遥感技术领域。所述的测量方法包含雷达图像前处理、风向测量和风速测量三部分，在风向测量指标中，将图像梯度、灰度和平滑项有机的结合起来，通过比例系数调节三者的比重，建立适合海面风场特征的模型，与现有技术相比，风向测量精度提高了68.4％以上。在风速测量指标中，当雷达单独测量时，将NRCS、实测风向、SNR作为BP网络输入，较传统算法风速精度提高了84％以上。在风速测量指标中，将海气边界层参数作为BP网络的附加输入，可进一步提高航海雷达测量风速的精度，考虑海气温差、盐度、潮位、大气压时，测量精度提高了48％以上。
文档编号G01W1/02GK102681033SQ20121012850
公开日2012年9月19日 申请日期2012年4月27日 优先权日2012年4月27日
发明者刘利强, 卢志忠, 戴运桃, 贾瑞才 申请人:哈尔滨工程大学