CN103941257A - 一种基于波数能量谱的导航雷达图像反演海面风向的方法 - Google Patents

Abstract

本发明属于利用遥感手段反演海洋表面风向的领域，特别涉及一种基于波数能量谱的导航雷达图像反演海面风向的方法。本发明包括：雷达图像预处理；提取海面静态特征；反演海面风向信息。本发明提出的一种基于波数能量谱的导航雷达图像反演海面风向的方法，不需要空间域缩减，能自动适应不同尺度的风条纹，可以在不同海况条件下测量风向，误差结果稳定，工程可用度高。本发明提出的一种基于波数能量谱的导航雷达图像反演海面风向的方法，通过波数能量谱波数尺度分离只将风条纹能量谱提取出来，提高了风向反演精度，降低了其他尺度噪声的干扰。

Description

一种基于波数能量谱的导航雷达图像反演海面风向的方法

技术领域

本发明属于利用遥感手段反演海洋表面风向的领域，特别涉及一种基于波数能量谱的导航雷达图像反演海面风向的方法。

背景技术

海面风场是海-汽交界面的重要动力来源，对航海作业和海洋动力学研究有重要的意义。海面风场信息主要包括海面风向和海面风速两个方面，本发明是用于反演海面风向信息的。

传统风向信息获取方法主要分为两大类：站点式现场测量和遥感测量。

站点式现场测量主要有岸基风向标单点测量、浮标、观测船上测量等，但这种测量方式易受周围遮挡环境和天气影响，海上测量受船舶和浮标运动等影响，并且对于远洋海域很难放置设备获取测量数据[1-3]。见参考文献[1-3]（段华敏.X波段雷达提取海面风场.中国海洋大学硕士论文.2009，1-3页.杨劲松.合成孔径雷达海面风场、海浪和内波遥感技术.北京：海洋出版社,2005：1-3,24-44页.Wu L C,Doong D J,Wu C C,etal.Wave and current extractionfrom marine radar images.Ocean Waves and Analysis,Proceedings of the Fifth InternationalSymposium on Ocean Wave Measurement and Analysis,Madrid Spain,2005.）

遥感测量主要有机载或星载的微波散射计、SAR等，这种手段可以测量大面积海面风场数据。但仍存在一些无法克服的缺陷：对于机载遥感技术，飞机航行受天气影响较大，在恶劣海洋气象状态下是无法工作的。对于卫星遥感测量技术，其空间分辨率低，无法获取近岸数据；时间采样率低，无法获得实时海洋数据；通过研究发现卫星云图上获得的多为云顶风而不是海表面风场信息[4-8]。见参考文献[4-8](Horstmann J,Koch W,Lehner S,etal.Oceanwinds from RADARSAT-1ScanSAE.Can J Remote Sens.2002,28(3):524-533P.Horstmann J,Lehner S and Rosenthal W.Detection of wave groups in SAR images and radar image sequences.IEEE Trans Geosci Remote Sens,2003,41:1437-1446P.Horstmann J,Koch W.Measurement ofocean surface winds using synthetic aperture radars.IEEE Journal of Oceanic Engineering,2005,30(3):508-515P.Horstmann J,Schiller H,Stellenfleth J S,etal.Global wind speed retrieval fromSAR.IEEE Transactions on Geoscience and Remote Sensing,2003,41(10):2277-2286P.Plant W J,Keller W C and Hayes K.Simulaneous measurement of ocean winds and waves with an airbornecoherent real aperture radar.Journal of Atmospheric and Oceanic Technology–special section,2004,22:832-846P.)

导航雷达是一种具有高时空分辨率的脉冲式雷达系统，通过旋转发射X波段微波，获取与海面发生后向散射的海面信号。X波段导航雷达因具有不受光线影响、不受天气影响、实时连续反馈、高分辨率和使用便捷等优点成为海面气象数据获取研究的热点。

应用X波段导航雷达图像测量海面风向国内外主要分为两类：一类是基于雷达回波强度和风向分布关系提取海面风向信息；一类是基于海面风速分布不均匀导致的雷达回波图像中呈现出的风条纹来计算海面风向信息。

1998年，Hatten等人提出基于雷达回波强度和风向分布关系提取海面风向信息，依据雷达观测方向与海面风向夹角不同时，雷达回波强度呈现非线性分布，从而获取海面风向信息。文献指出HH极化X波段雷达观测方向是逆风向时雷达回波强度达到最大，从而确定逆风方向来获取海表面主风向[9]。但是此方法需要360°全幅无遮挡的X波段雷达图像，但无论岸基、塔基还是航海雷达都是无法实现全幅无遮挡探测的，因此，此方法在工程和实验中都是无法实现的。见参考文献[9]（Hatten H,Seemann J,Horstmann J,etal.Azimuthal dependence ofthe radar cross section and the spectral background noise of a nautical radar at grazing incidence.International Geoscience and Remote Sensing Symposium.Seattle,USA,1998,5:2490-2492P）

另一类是应用雷达回波强度与风速的关系，主要是基于海面风速分布不均匀导致的雷达回波图像中呈现出的风条纹来计算海面风向信息。雷达回波上存在的风条纹成因如下：

当风掠过海面时会产生风切应力，海-气界面上会受到不均匀的风切应力作用，由于风切应力的不同会产生分布不均匀的阵风风速场叠加在海面上，由于风速的作用在海面上会产生不同尺度的微湍流，这些微湍流改变了海表面的粗糙度。因此，根据粗糙度的分布在X波段导航雷达图像上会形成亮暗相间的条纹，即风条纹。风条纹的尺度为100～500m，频率接近静态或准静态，风条纹方向与风向平行[10]。见参考文献[10]（Moeng C H,Sullivan P.AComparison of Shear-and Bouyancy-Driven Planetary Boundary Layer Flower[M].）。

2004年，Dankert等人提出应用光流法从X波段导航雷达图像序列中提取海面风向信息，利用提取的风条纹在多幅图像上的移动来确定风场信息，这就要求在连续的图像上存在灰度一致的风条纹[11]。但风条纹信号是静态信号，根据风条纹的成因在同一序列图像中无法获得风条纹的移动，因此光流法在理论和实验中都证明是无法反演海面风向的。见参考文献[11]（Dankert H,Horstmann J,Rosenthal W.Ocean surface winds retrieved from marine radar imagesequences.International Geosciences and Remote Sensing Symposium,2004,3:1903-1906P）

2003年，Dankert等人提出应用局部梯度法从X波段导航雷达图像序列中测量海面风向。此算法首先需要从X波段导航雷达图像序列中提取出海面静态特征图像，再对海面静态特征图像进行平滑和缩减得到风条纹图像，然后应用Sobel算子来求取风条纹图像梯度方向，由于风向与风条纹方向平行，梯度的垂直方向就是海面风向。Dankert等人已经应用来自Ekofish2/4和FINO-I平台的数据验证了应用此算法计算海面风向信息的可行性[12]。见参考文献[12]（Dankert H,Horstmann J,Rosenthal W.Ocean wind fields retrieved from radar image sequences.Journal of Geophysical Research,2003,108(C11):16-1-16-11P）

局部梯度法是对风条纹图像进行空间域处理的，通过实验发现此方法计算海面风向的关键在于对海面静态特征图像缩减和平滑的次数，只有缩减后图像分辨率为风条纹尺度的1/16～1/4时，才能计算出正确的风条纹的梯度，从而获得准确的海面风向。但风条纹的尺度是在100～500之间不定的，在工程中无法自动寻找到合适的缩减图像分辨率，而固定的缩减图像分辨率会导致风向反演误差不稳定，进而导致此算法无法实际应用。

本发明提出了一种基于波数能量谱的导航雷达图像反演海面风向的方法，此方法不需要对海面静态特征图像进行空间域的缩减和平滑处理，只需根据风条纹的特性就可以将其能量谱特征提取出来，从而计算出海面风向，既可以满足工程对海面风向的精度又可以程序化应用到工程中去。

发明内容

本发明的目的在于提供一种保证海面风向反演精度，又能适用工程实时反馈海面风向信息的基于波数能量谱的导航雷达图像反演海面风向的方法。

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

（1）雷达图像预处理：

对实测X波段导航雷达图像应用2-D非线性平滑中值滤波，

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

g(s,t)为雷达图像像元点在极坐标位置(s,t)的图像回波强度值；f'(r,θ)为滤波后图像在极坐标位置(r,θ)的灰度值；N(r,θ)为中心点在(r,θ)处的像元点，(s,t)取以(r,θ)为中心的8个点，将3×3模板中值滤波器的模板中心与极坐标图像的像元位置重合，与周围8个相邻像元点的回波强度值进行排列，取中间的回波强度值赋给中心位置的像元，模板遍历全幅雷达图像获得中值滤波后的图像序列；

（2）提取海面静态特征：

应用图像序列建立全局低通滤波器，滤除高频信号，保留静态或准静态频率信号：

（2.1）极坐标图像归一化，使图像序列中的每幅图像的线数固定，固定像元点数和像元位置：

（2.1.1）建立有固定N条线和M个同心圆的极坐标新图像；

（2.1.2）将原图像角度值与新角度值相等，或大于原图像角度值的第一条线上的雷达回波强度赋到新图像的线上；

（2.1.3）重复（2.1.2）直到新图像上的N条线上都具有原雷达回波强度值，获得归一化的新极坐标图像；

（2.2）应用X波段导航雷达图像序列建立全局低通滤波器，提取海面静态特征信号：

（2.2.1）根据海面风条纹属于静态频率特性，在单幅雷达图像上无法获得风条纹，应用归一化后的极坐标图像序列，对其在同一位置的像元点上按时间进行积分平均，构建图像空间全局低通滤波器：

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

f(θ,r)为极坐标海面静态特征图像，f'(θ,r,t)图像序列中时间为t的单幅图像，N_t为时间序列；

（2.2.2）构建N条线和M个同心圆的二维极坐标图像；

（2.2.3）将获得的M×N个低通滤波后的新像元灰度值按位置赋予新构建的二维极坐标图像上，得到二维极坐标海面静态特征图像；

（3）反演海面风向信息：

（3.1）极坐标图像区域选取和插值为笛卡尔坐标图像，在二维极坐标海面静态特征图像上选取适当大小的扇形极坐标图像，在选取区域内的图像灰度值为f(r_i,θ_j)(i,j=1,2...n)：

（3.1.1）根据选取区域的极坐标角度θ_j及半径r_i信息构建笛卡尔坐标系 $y_{N_y}=r_j\sin \left({\mathrm{\θ}}_i\right)(N_x{N}_y=\mathrm{1,2}...n);$

（3.1.2）选取区域所在的中心角度θ_c和半径r_c，计算选取极坐标图像中中心点在全幅极坐标图像中的笛卡尔坐标位置(x_c,y_c)，以中心点所在位置为中心计算选取区域内的像元点在全幅极坐标图像中的笛卡尔坐标的x-y轴分量(x_i,y_j)；

（3.1.3）应用最近点插值，当两个坐标(x_i,y_j)和位置最近时将f(r_i,θ_j)赋到上，将选取的极坐标区域图像上的灰度值都插值到新生成的笛卡尔坐标上，横纵坐标为像元点所在位置，得到笛卡尔坐标下的二维静态特征图像；

（3.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{x_{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)为的傅里叶系数，复数指数项为：

$\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)})\\ +\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)为f(x,y)在能量谱域Τ的坐标，d_x、d_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]}$

A(k_x,k_y)为二维静态特征能量谱值；

（3.3）应用能量谱域波数尺度分离带通滤波器提取出风条纹能量谱图：

根据海面条纹波长与能量谱波数存在谱域中存在的关系：

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

k为能量谱波数，(k_x,k_y)为能量谱域Τ的坐标，即k为在x和y轴的波数分量；L为雷达图像呈现的海面条纹波长，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)为风条纹能量谱；

（3.4）根据风条纹能量谱提取海面主风向：

得到的风条纹能量谱I(k_x,k_y)关于一三或二四象限镜像对称，两个能量集中区域连线方向为风条纹的垂直向，而连线的垂直方向则为风条纹平行方向，由于风条纹方向与风向平行，风条纹方向就为海面主风向，

N_w=|θ_c|+|α|+|β|；

N_w为相对正北向的风向；θ_c为选取区域的中心角；β为笛卡尔坐标下计算得到的风向；α为船艏向；将计算出来的方向与风向标测得的风向所在象限比较，保留与其象限一致的方向去除180°模糊的方向，获得准确的海面风向。

本发明的有益效果在于：

本发明提出的一种基于波数能量谱的导航雷达图像反演海面风向的方法，不需要空间域缩减，能自动适应不同尺度的风条纹，可以在不同海况条件下测量风向，误差结果稳定，工程可用度高。本发明提出的一种基于波数能量谱的导航雷达图像反演海面风向的方法，通过波数能量谱波数尺度分离只将风条纹能量谱提取出来，提高了风向反演精度，降低了其他尺度噪声的干扰。较现有的局部梯度法反演海面，风向反演精度提高了58.97%，完全达到了工程要求。本发明提出的一种基于波数能量谱的导航雷达图像反演海面风向的方法，本发明方法只针对的是风条纹能量谱，不需要适应条纹尺度变化，可以提高运算时间。本发明提出的一种基于波数能量谱的导航雷达图像反演海面风向的方法，通过全局低通滤波器将包含条纹信号的海面静态特征信号提取出来，采用极坐标图像归一化技术解决了雷达图像像元数目不定的工程问题。

附图说明

图1本发明提供的海面风向测量方法流程图；

图2为本发明中值滤波前后雷达图；

图3本发明提取出的静态特征雷达图像；

图4为极坐标图像区域划分和区域选取；

图5为本发明中极坐标插值笛卡尔坐标示意图；

图6为笛卡尔坐标下的风场信息图像；

图7为本发明中得到风条纹功率谱图；

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

图9为现有局部梯度法测得风向与参考风向误差分布结果。

具体实施方式

下面结合附图对本发明做进一步描述。

本发明公开了一种基于波数能量谱的导航雷达图像反演海面风向的方法，属于利用遥感手段反演海洋表面风向领域。所述的测量方法包含雷达图像预处理、海面静态特征提取和海面风向信息反演三个部分，通过建立全局低通滤波器将雷达图像中海面静态特征信号提取出来，为后续海面风场的研究奠定了基础。通过二维离散快速傅里叶变换得到二维静态特征的能量谱图，结合风条纹的尺度为100～500m，建立波数能量谱带通滤波器将风条纹能量谱分离出来，根据风条纹能量谱分布特性确定出测量海域的平均主风向。较现有的局部梯度算法比，精度提高了58.97%，抑制了反演方法的不稳定，反演海面风向适用性达到了99%，足以达到工程和环境监测要求。

本发明的海面风向信息提取方法主要分为三个部分：雷达图像预处理、海面静态特征提取和海面风向信息反演，三个部分中海面风向信息反演是重点发明内容。所述的海面风向信息提取方法具体步骤如下：

（一）雷达图像预处理：

在应用X波段导航雷达图像序列反演海面风向前需要对实测图像进行预处理，对实测雷达图像序列应用3×3模板的2-D非线性平滑中值滤波，从而抑制同频对海面风向提取的影响。

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

式(1)中g(s,t)为雷达图像像元点在极坐标位置(s,t)的图像回波强度值；f'(r,θ)为滤波后图像在极坐标位置(r,θ)的灰度值；N(r,θ)为中心点在(r,θ)处的像元点，(s,t)取以(r,θ)为中心的8个点。

（二）海面静态特征提取：

为了反演海面风向信息需要获得海面风切应力导致的风条纹图像，根据海面风条纹属于静态频率特性，在单幅雷达图像上无法获得风条纹，需要对雷达图像序列低通滤波处理。

获取海面静态特征图像需要对雷达图像序列进行积分平均，这就要求图像序列中的每幅图像的线数固定，因此要对极坐标图像归一化。首先建立固定线数的极坐标，再将与线所在角度最近的原图像线的回波强度值赋到新线上，从而达到像元点数和像元位置都是固定的雷达图像。

本发明构建全局低通滤波器可以将高频的海浪信号等滤除，只保留静态或半静态的海面信号，主要包含风条纹信号。为获得海面静态特征构建的全局低通滤波器的核心算法为：

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

式中f(θ,r)为海面静态特征图像，在笛卡尔坐标可以表示为f'(θ,r,t)为归一化图像序列中时间为t的单幅雷达图像，N_t为时间序列。

（三）海面风向信息反演：

本发明的风向信息反演是基于风条纹图像波数能量谱特征的一种方法。对二维图像进行二维离散快速傅里叶变换（2D FFT）及可得到图像的能量谱图，从中可以获得图像像元灰度值频率分布特性。所述的风向信息反演是基于图像能量谱特征的，首先要对获得包含风条纹的图像进行二维离散快速傅里叶变换，即：

$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{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})]---\left(3\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)})\\ +\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(4\right)$

其中，

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

$k_y=\frac{2\mathrm{\π}}{\max \left(y_{N_y}\right)*d_y}---\left(6\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(7\right)$

海面条纹波长与能量谱波数存在的关系式为：

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

式(8)中k为能量谱波数，(k_x,k_y)为能量谱域Τ的坐标，即k为在x和y轴的波数分量；L为雷达图像呈现的海面条纹波长。

由于风条纹的波长存在一定的范围，假设L_d为风条纹波长尺度下限，L_t为风条纹波长尺度上限，根据公式(8)可以获得风条纹能量谱波数下限k_d：

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

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

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

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

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

由于傅里叶变换的周期性，得到的风条纹能量谱关于一三或二四象限镜像对称，两个能量集中区域连线方向为风条纹的垂直向，而连线的垂直方向则为风条纹平行方向，由于风条纹方向与风向平行，所以可得到海面风向。

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

N_w=|θ_c|+|α|+|β|    (12)

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

由于垂线向有两个方向，因此，计算出的风向存在180°模糊问题。由于单点风向标测得的海面风向所属的象限是固定的，将计算出来的方向与风向标测得的风向所在象限比较，保留与其象限一致的方向去除180°模糊的方向，从而获得准确的海面风向。

本发明提供的基于波数能量谱的导航雷达图像反演海面风向的方法，流程图如1所示，本发明的海面风向提取方法包含雷达图像预处理，海面静态特征信号提取和海面风向信息反演三部分，下面第一步为雷达图像预处理，第二步到第三步为海面静态特征信号提取，第四步到第九步为海面风向信息反演，具体步骤如下：

第一步，对实测雷达图像序列应用3×3模板的2-D非线性平滑中值滤波，抑制同频对海面风场研究的影响。

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

式(1)中g(s,t)为雷达图像像元点在极坐标位置(s,t)的图像回波强度值；f'(r,θ)为滤波后图像在极坐标位置(r,θ)的灰度值；N(r,θ)为中心点在(r,θ)处的像元点，(s,t)取以(r,θ)为中心的8个点。

将3×3模板中值滤波器的模板中心与极坐标图像的某个像元位置重合，将其与周围8个相邻像元点的回波强度值进行排列，取中间的回波强度值赋给中心位置的像元，模板遍历全幅雷达图像获得中值滤波后的图像序列，图2为本发明中值滤波前后雷达图像。

第二步，极坐标图像归一化，使图像序列中的每幅图像的线数固定，从而达到像元点数和像元位置固定。

本发明应用的X波段导航雷达图像是由大约3600条射线状线（3600个角度值）和600个同心圆组成，间隔角度大约是0.1°，由于雷达工作时脉冲数量不定和天线旋转时受到各种外界环境的干扰，所以线数是不固定的。为了使图像的线数固定即间隔角度固定使用极坐标图像归一化算法，步骤如下：

①建立一个1800条线（1800个角度值）和600个同心圆组成的极坐标新图像；

②将原图像角度值与新角度值相等，或大于原图像角度值的第一条线上的雷达回波强度赋到新图像的线上；

③不断重复②知道新图像上的N条线上都具有原雷达回波强度值，从而获得归一化的新极坐标图像。

第三步，应用雷达图像序列建立全局低通滤波器，滤除海浪等高频信号，仅保留静态或准静态频率特征信号，主要为风场导致的风条纹信号。具体步骤如下：

①根据海面风条纹属于静态频率特性，在单幅雷达图像上无法获得风条纹，对归一化的雷达图像序列在同一位置的像元点上按时间(80s)进行积分平均，构建一个图像空间全局低通滤波器：

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

式中f(θ,r)为风条纹图像，f'(θ,r,t)雷达图像序列时间为t的单幅雷达图像，N_t为时间序列；

②构建1800条线和600个同心圆的二维极坐标图像；

③将获得的1800×600个低通滤波后的新像元灰度值按位置赋予新构建的二维极坐标图像上，得到二维极坐标海面静态特征图像，如图3所示。

第四步，极坐标图像区域划分，将二维极坐标图像的像元点按照128*128的数量将图像分成14个小区域，如图4所示。计算每个小区域的海面风向，最终海面平均风向为每个小区域风向的平均值。图4中红框区域的中心角度为203°，像元灰度值为f(r_i,θ_j)(i,j=1,2...128)。

第五步，图4中红框区域的极坐标应用最近点插值为笛卡尔坐标，主要是将原区域的雷达回波强度值应用位置最近的方式，将值赋予到已建立的笛卡尔坐标下，具体步骤如下：

①根据图4中扇形区域的极坐标角度θ_j及半径ri信息构建笛卡尔坐标系 $y_{N_y}=r_j\sin \left({\mathrm{\θ}}_i\right)(N_{x,}{N}_y=\mathrm{1,2}...128);$

②选取区域所在的中心角度θ_c和半径r_c，计算选取极坐标图像中中心点在全幅极坐标图像中的笛卡尔坐标位置(x_c,y_c)，以中心点所在位置为中心计算选取区域内的像元点在全幅极坐标图像中的笛卡尔坐标的x-y轴分量(x_i,y_j)；

③应用最近点插值，当两个坐标(x_i,y_j)和位置最近时将f(r_i,θ_j)赋到上，将选取的极坐标区域图像上的灰度值都插值到新生成的笛卡尔坐标上，横纵坐标为像元点所在位置，这样就得到笛卡尔坐标下的二维静态特征图像，插值方法如示意图5。图6为图4中红框区域插值结果图，横纵坐标为128*128像元点所在位置。

第六步，获取海面静态特征二维能量谱图，对二维图像进行二维离散快速傅里叶变换（2DFFT）及可得到图像的能量谱图。本发明对二维笛卡尔坐标下的海面静态特征图像应用二维离散快速傅里叶变换（2D FFT）获得海面静态特征图像能量谱，数学模型为：

$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{x_{N_x}*k_x}{\max \left(x_{N_x}\right)}+\frac{y_{N_y}*k_y}{\max \left(y_{N_y}\right)})]---\left(3\right)$

F(k_x,k_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)})\\ +\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(4\right)$

其中，

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

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

(k_x,k_y)为f(x,y)在能量谱Τ的坐标，d_x、d_y为雷达图像分辨率，本发明中公式d_x=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(7\right)$

第七步，对静态特征图像能量谱应用波数能量谱尺度分离，将风条纹信号的谱能量从静态特征图像能量谱中分离出来，主要是依据风条纹波长L和频域坐标尺度k的关系来提取。

由于风条纹的尺度为100～500m，计算得到风条纹能量谱波数下限k_d：

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

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

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

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

I(k_x,k_y)为风条纹能量谱，图7波数尺度分离后I(k_x,k_y)在(k_x,k_y)坐标下得到的风条纹等高线能量谱图，虚线为(k_x,k_y)坐标下的波长范围。

第八步，获得海面方向信息。由于傅里叶变换的周期性，得到的风条纹能量谱关于一三或二四象限镜像对称，两个能量集中区域连线方向为风条纹的垂直向，而连线的垂直方向则为风条纹平行方向，由于风条纹方向与风向平行，即可得到海面方向。

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

N_w=|θ_c|+|α|+|β|(11)

N_w为相对正北向的海面风向；θ_c为选取区域的中心角，选取为110^°；β为笛卡尔坐标下计算得到的风条纹方向；α为船艏向，为93°。图6中得到N_w为33°或213°。

由于单点风向标测得的海面风向为36°属于第一个象限，将计算出来两个方向与第一个象限比较，保留第一个象限角度，去除第三象限的角度，从而获得的海面风向为33°。

第九步，获取主海面风向信息。对第四步得到的每个小区域重复执行第五步到第八步，得到每个小区域的相对正北向的海表面风向最后将结果进行矢量平均得到X波段导航雷达图像序列测得的平均海面风向值：

$N_{\mathrm{wa}}=\underset{i}{\mathrm{\Σ}}{N}_w^i(i=\mathrm{1,2,3}...14)---\left(12\right)$

N_wa为最终计算得到的海面主风向，区域划分的越多则平均风向精度会有所增加，但要保证在进行二维离散快速傅里叶变换时像元的数量，不然在谱图像上无法获得风条纹尺度的波数k，这样就无法获得准确的海面风向信息。

本发明公开的是一种基于波数能量谱的导航雷达图像反演海面风向开展实验，实验数据来源于福建平潭县海坛岛，雷达平台安装在海坛岛岸基，此海域平均海深25m，受地形影响此海域经常出现高海况（>4级）。雷达安装高度为40m，旋转周期为2.4s，发射9.5GHz的X波段微波信号，采用HH极化，掠入射角为1°～5°，当遇到微波波长1/2的毛细波时产生Bragg散射，回波信号将返回到雷达接收机上形成海杂波图像。雷达图像覆盖海域为2km，每32幅图像为一个序列（76.8s），径向和幅向分辨率都为7.5m。实验数据从2010年10月到11月，共采集1634组数据。由于实验期间台风“鲶鱼”的影响，风向多为东北风，出现短暂的西南风。参考数据来源于同位置放置的风速计和风向标，以每分钟记录实时海面风向和风速值。

为验证实验结果将本发明测得结果以下称“ESM(energy spectrum method)风向”，现有方法局部梯度法测得结果称为“LGM(local gradient method)风向”，局部梯度法选用固定缩减图像分辨率为30m，风向标测得风向以下称为“参考风向”。两两海面风向结果误差分布相比较结果如图8、图9，两个图比较可以发现ESM风向与参考风向的偏差更小，说明ESM风向与参考风向更为吻合。LGM风向与参考风向有存在40°以上的误差，主要是由于空间域处理时缩减的次数固定，没有与风条纹的尺度进行比对，导致误差不稳定。两种方法海面风向误差统计如表1。

表1海面风向误差统计

本发明提供的海面风向测量方法，是基于风条纹图像波数域的海面风向测量方法。首先，设计了全局低通滤波器从雷达图像序列中提取出包含风条纹信息的海面静态特征图像；其次，应用二维离散快速傅里叶变换得到海面静态特征能量谱；最后，针对风条纹特征，设计了频域尺度分离带通滤波器，将从海面静态特征能量谱中提取出风条纹能量谱，根据谱峰分布位置确定海面风向。较现有的LGM测量海面风向精度有所提高，避免了对风条纹图像进行空间域缩减和平滑处理不当造成反演精度的降低。如表1所示，本发明测得的海面风向与参考风向的相关系数达到了0.99，偏差为-0.64°，本发明测得的海面风向完全达到工程应用需求。

Claims (4)

1.一种基于波数能量谱的导航雷达图像反演海面风向的方法，其特征在于： 

（1）雷达图像预处理： 

对实测X波段导航雷达图像应用2-D非线性平滑中值滤波， 

g(s,t)为雷达图像像元点在极坐标位置(s,t)的图像回波强度值；f'(r,θ)为滤波后图像在极坐标位置(r,θ)的灰度值；N(r,θ)为中心点在(r,θ)处的像元点，(s,t)取以(r,θ)为中心的8个点，将3×3模板中值滤波器的模板中心与极坐标图像的像元位置重合，与周围8个相邻像元点的回波强度值进行排列，取中间的回波强度值赋给中心位置的像元，模板遍历全幅雷达图像获得中值滤波后的图像序列； 

（2）提取海面静态特征： 

应用图像序列建立全局低通滤波器，滤除高频信号，保留静态或准静态频率信号： 

（2.1）极坐标图像归一化，使图像序列中的每幅图像的线数固定，固定像元点数和像元位置： 

（2.1.1）建立有固定N条线和M个同心圆的极坐标新图像； 

（2.1.2）将原图像角度值与新角度值相等，或大于原图像角度值的第一条线上的雷达回波强度赋到新图像的线上； 

（2.1.3）重复（2.1.2）直到新图像上的N条线上都具有原雷达回波强度值，获得归一化的新极坐标图像； 

（2.2）应用X波段导航雷达图像序列建立全局低通滤波器，提取海面静态特征信号： 

（2.2.1）根据海面风条纹属于静态频率特性，在单幅雷达图像上无法获得风条纹，应用归一化后的极坐标图像序列，对其在同一位置的像元点上按时间进行积分平均，构建图像空间全局低通滤波器： 

f(θ,r)为极坐标海面静态特征图像，f'(θ,r,t)图像序列中时间为t的单幅图像，N_t为时间序列； 

（2.2.2）构建N条线和M个同心圆的二维极坐标图像； 

（2.2.3）将获得的M×N个低通滤波后的新像元灰度值按位置赋予新构建的二维极坐标图像上，得到二维极坐标海面静态特征图像； 

（3）反演海面风向信息： 

（3.1）极坐标图像区域选取和插值为笛卡尔坐标图像，在二维极坐标海面静态特征图像上选取适当大小的扇形极坐标图像，在选取区域内的图像灰度值为f(r_i,θ_j)(i,j＝1,2...n)： 

（3.1.1）根据选取区域的极坐标角度θ_j及半径r_i信息构建笛卡尔坐标系

（3.1.2）选取区域所在的中心角度θ_c和半径r_c，计算选取极坐标图像中中心点在全幅极坐标图像中的笛卡尔坐标位置(x_c,y_c)，以中心点所在位置为中心计算选取区域内的像元点在全幅极坐标图像中的笛卡尔坐标的x-y轴分量(x_i,y_j)； 

（3.1.3）应用最近点插值，当两个坐标(x_i,y_j)和位置最近时将f(r_i,θ_j)赋到 上，将选取的极坐标区域图像上的灰度值都插值到新生成的笛卡尔坐标上，横纵坐标为像元点所在位置，得到笛卡尔坐标下的二维静态特征图像； 

（3.2）获取静态特征二维能量谱，将图像从空间域转化到谱域进行处理，通过二维离散快速傅里叶变换实现，对笛卡尔坐标下的二维静态特征图像进行二维离散快速傅里叶变换，即： 

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

其中， 

(k_x,k_y)为f(x,y)在能量谱域Τ的坐标，d_x、d_y为雷达图像分辨率，二维谱性质得到 的能量谱为： 

A(k_x,k_y)为二维静态特征能量谱值； 

（3.3）应用能量谱域波数尺度分离带通滤波器提取出风条纹能量谱图： 

根据海面条纹波长与能量谱波数存在谱域中存在的关系： 

k为能量谱波数，(k_x,k_y)为能量谱域Τ的坐标，即k为在x和y轴的波数分量；L为雷达图像呈现的海面条纹波长，L_d为风条纹波长尺度下限，L_t为风条纹波长尺度上限，获得风条纹能量谱波数下限k_d： 

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

将风条纹信号功率谱提取出来： 

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

（3.4）根据风条纹能量谱提取海面主风向： 

得到的风条纹能量谱I(k_x,k_y)关于一三或二四象限镜像对称，两个能量集中区域连线方向为风条纹的垂直向，而连线的垂直方向则为风条纹平行方向，由于风条纹方向与风向平行，风条纹方向就为海面主风向， 

N_w=|θ_c|+|α|+|β|； 

N_w为相对正北向的风向；θ_c为选取区域的中心角；β为笛卡尔坐标下计算得到的风向；α为船艏向；将计算出来的方向与风向标测得的风向所在象限比较，保留与其象限一致的方向去除180°模糊的方向，获得准确的海面风向。 

2.根据权利要求1所述的一种基于波数能量谱的导航雷达图像反演海面风向的方法，其特征在于：计算风条纹波长尺度下限L_d取500m。 

3.根据权利要求1所述的一种基于波数能量谱的导航雷达图像反演海面风向的方法，其特征在于：计算风条纹波长尺度上限L_t取100m。 

4.根据权利要求1所述的一种基于波数能量谱的导航雷达图像反演海面风向的方法，其特征在于：船艏向角α取93°。