CN104297753B - 一种基于自适应缩减算子的导航雷达图像反演海面风向方法 - Google Patents

Abstract

本发明公开了一种基于自适应缩减算子的导航雷达图像反演海面风向方法。包括以下几个步骤：采集N幅导航雷达图像形成一组导航雷达图像序列；对导航雷达图像进行中值滤波处理；对滤波后的导航雷达图像进行方位向归一化，固定每幅滤波后的导航雷达图像的方位向线数；对方位向归一化后的导航雷达图像进行全局低通滤波处理，得到包含风条纹的海面静态特征图像；根据海面静态特征图像，进行海面风向反演得到海面风向。本发明设计了一种自适应算法，可以根据缩减后图像计算得到的判定参数判断缩减分辨率与风条纹尺度是否适宜，提高了海面风向反演精度。

Description

一种基于自适应缩减算子的导航雷达图像反演海面风向方法

技术领域

本发明属于海面风向遥感技术领域，尤其涉及一种基于自适应缩减算子的导航雷达图像反演海面风向方法。

背景技术

海面风场信息是海洋动力学重要参数，主要包含风向和风速信息两个方面，因此，了解和掌握海面风向信息具有深远的意义。现有测量海面风向信息方法有站点式现场测量和遥感测量两类。导航雷达是遥感测量手段的一种，因其具有不受光线影响、不受天气影响、实时连续反馈、高分辨率和使用便捷等优点，成为现阶段获取海面风向信息的重要手段。

现阶段应用导航雷达图像获取海面风向的方法主要有两类：一类是基于导航雷达回波强度和风向分布关系，另一类是基于风条纹反演海面风向。基于导航雷达回波强度和风向分布关系获取海面风向需要360°全幅无遮挡的导航雷达图像，但无论岸基、塔基还是航海导航雷达都是无法实现全幅无遮挡探测，因此，此方法在工程和实验中都是无法实现的。风条纹是由海面风速分布不均匀导致的导航雷达回波图像中呈现出来的条纹特征，本发明属于基于风条纹反演海面风向这类。见参考文献Moeng C H,Sullivan P.AComparison of Shear-and Bouyancy-Driven Planetary Boundary Layer Flower[M].American Meterorological Society,1994.Hatten H,Seemann J,Horstmann J,etal.Azimuthal dependence of the radar cross section and the spectralbackground noise of a nautical radar at grazing incidence[C]//Geoscience andRemote Sensing Symposium Proceedings,1998.IGARSS'98.1998IEEEInternational.IEEE,1998,5:2490-2492。

风条纹具有尺度200～500m，频率接近静态或准静态，风条纹方向与风向平行等特征。目前，国内外基于风条纹特征反演海面风向的算法主要有以下几种：2004年由Dankert等人提出的光流法，2003年Dankert等人提出局部梯度法。2010年中国海洋大学段华敏应用光流法从小麦岛导航雷达数据中提取出了海面风场信息。2012年哈尔滨工程大学的贾瑞才博士和2013年中国海洋大学的硕士生李金凤都对光流法进行了改进，用以提高了风向反演精度。光流法要求在连续的图像上存在灰度一致的风条纹，但根据风条纹的成因在同一序列图像中无法获得风条纹的移动，因此光流法在理论和实验中都证明是无法反演海面风向的。见参考文献Dankert H，Horstmann J，etal.Ocean wind fields retrieved fromradar-image sequences.International Geoscience and Remote Sensing Symposium(IGARSS),v4,p2150-2152,2002.Dankert H,Horstmann J,Rosenthal W.Ocean surfacewinds retrieved from marine radar image sequences.International Geosciencesand Remote Sensing Symposium,2004,3:1903-1906P.Dankert H,Horstmann J,Rosenthal W.Ocean wind fields retrieved from radar image sequences.Journal ofGeophysical Research,2003,108(C11):16-1-16-11P。

2003年Dankert应用局部梯度法反演海面风向，首先要对海面静态特征图像进行平滑和缩减三次得到风条纹图像，然后根据风条纹与风向平行这个特征，通过计算像元主梯度方向，其垂直方向则是海面风向，具体步骤流程如图2所示；2003年Dankert应用局部梯度法利用Ekofish2/4k平台数据，得到海面风向与参考风向标准差为14.24°。2005年Dankert利用FINO-I平台数据，得到海面风向与参考风向标准差为12.77°。2007年Dankert利用FINO-I平台更多数据，得到海面风向与参考风向标准差为12.77°，通过实验和理论分析可知应用局部梯度法是可以从导航雷达图像中反演出海面风向的。见参考文献DankertH,Horstmann J.Wind measurements at FINO-I using marine radar-image sequences[C]//Geoscience and Remote Sensing Symposium,2005.IGARSS'05.Proceedings.2005IEEE International.IEEE,2005,7:4777-4780.Dankert H.A marine radar windsensor.Journal of Atmospheric and Oceanic Technology,2007,24:1629-1642P。

传统局部梯度法对海面静态特征图像应用2×2缩减算子固定的缩减三次，图像分辨率降低为原来的8倍，即得到固定的缩减分辨率。Dankert实验指出将图像分辨率缩减为风条纹尺度的1/16～1/4时，才能计算出正确的风条纹的梯度，从而获得准确的海面风向。但风条纹的尺度在200～500m之间，固定的缩减次数导致缩减分辨率不一定在与风条纹合适的比例范围内，导致大量实验数据不可应用。分辨率与风条纹的比例是在一个范围内，并且风条纹尺度会根据不同的海况而改变，很难寻找到与风条纹尺度合适的比例分辨率，这样也会导致反演得到的海面风向精度降低。因此，传统的局部梯度法很难满足工程的要求，需要设计一种可以根据不同风条纹尺度来调节缩减次数及分辨率大小的方法。

发明内容

本发明的目的是提供一种具有高精度的，基于自适应缩减算子的导航雷达图像反演海面风向方法。

一种基于自适应缩减算子的导航雷达图像反演海面风向方法，包括以下几个步骤：

步骤一：采集N幅导航雷达图像形成一组导航雷达图像序列；

步骤二：对导航雷达图像进行中值滤波处理；

步骤三：对滤波后的导航雷达图像进行方位向归一化，固定每幅滤波后的导航雷达图像的方位向线数；

步骤四：对方位向归一化后的导航雷达图像进行全局低通滤波处理，得到包含风条纹的海面静态特征图像；

步骤五：根据海面静态特征图像，进行海面风向反演得到海面风向。

本发明一种基于自适应缩减算子的导航雷达图像反演海面风向方法，还可以包括：

1、对滤波后的导航雷达图像进行方位向归一化的方法为：

(1)读取极坐标下的每个滤波后的导航雷达图像的方位向线数和径向点数，方位向线数为3600条，间隔角度为0.1°，有N＝3600个角度值Ω_i,i＝1,2,…N，径向点数为220个；

(2)建立方位向固定为1800条的新极坐标导航雷达图像，间隔角度为0.2°，有N_new＝1800个新的角度值，θ_j,j＝1,2,…N_new，径向点数为220个；

(3)为新极坐标导航雷达图像赋灰度值，若Ω_i＝θ_j，或第一个Ω_i>θ_j，则将角度值Ω_i对应的方位向线的灰度值赋给新极坐标导航雷达图像的新角度值θ_j对应的方位向线上；

重复步骤(3)直到所有新极坐标导航雷达图像的所有方位向线都具有灰度值，得到方位向归一化后的导航雷达图像。

2、根据海面静态特征图像，进行海面风向反演得到海面风向的方法为：

(1)将极坐标下的海面静态特征图像插值为笛卡尔坐标下海面静态特征图像；

(2)将笛卡尔坐标下海面静态特征图像进行平滑处理，得到一次平滑图像；

笛卡尔坐标下海面静态特征图像为：

f(x_i,y_j) i＝1,2,…N_x,j＝1,2,…N_y

其中(x_i,y_j)为沿x和y轴的笛卡尔坐标，N_x和N_y为笛卡尔坐标下沿x和y轴所取像元总数，

一次平滑图像为：

H^r(m,n)为二项式卷积核，r为二项式卷积核的阶数，(m,n)为二项式卷积核的坐标；

(3)应用自适应缩减算子对一次平滑图像进行缩减，得到缩减后图像；

缩减后图像为：

F^(K)＝C_(↓K)*F

自适应算子C_(↓K)为：

进一步得到：

其中，(x_α,y_β)为图像缩减后新生成的坐标，α＝1,2,…N_x-4/2，β＝1,2,…N_y-4/2，K为缩减率，缩减后图像分辨率变为K*γ，γ为图像缩减前图像分辨率；

(4)将缩减后图像进行平滑处理，得到二次平滑图像；

二次平滑图像为：

其中，(x_α′,y_β′)二次平滑图像坐标，H^R为二项式卷积核，R为二项式卷积核阶数，(M,N)为二项式卷积核坐标，

(5)对二次平滑图像应用优化Sobel梯度算子得到像元梯度方向直方图；

优化Sobel梯度算子为：

其中，D_x和D_y为优化Sobel算子分别沿x，y轴的梯度算子，

二次平滑图像每个像元点沿x和y轴的梯度值为：

其中，G_x和G_y分别为所有像元点沿x轴和y轴的梯度方向，由G_x和G_y得到每个像元点的梯度方向G_θ为：

对得到的所有像元的梯度方向进行直方图统计，得到像元梯度方向直方图；

(6)通过像元梯度方向直方图得到稳定系数，基于自适应算法确定最优缩减率K；

从像元梯度方向直方图中选取从从的所有梯度方向，组成方向值数组G_P＝[G₁,G₂,…G_p]，

稳定系数为：

其中，σ为G_P的标准差，μ为G_P的均值，

选取阈值为10，缩减率K的取值范围为：K＝2,3,…8，

在缩减率K<8，并且稳定系数小于阈值时，选取当前缩减率K作为最优缩减率

在缩减率K<8，并且稳定系数大于阈值时，K＝K+1，重复步骤(2)～步骤(6)，

在缩减率K＝8，并且稳定系数大于阈值时，选取稳定系数最小时对应的缩减率K作为最优缩减率；

(7)求得对应最优缩减率K的主梯度方向；

将当前的对应最优缩减率K的像元梯度方向直方图中选取的所有梯度方向，组成的新的方向值数组G_M＝[G₁,G₂,…G_m]，包含的梯度方向数量为m，对新的方向值数组求平均值得到主梯度方向G_main：

(8)根据主梯度方向求得海面风向W_d：

W_d＝G_main±90^。

将计算出来的方向与风向标测得的风向所在象限比较，保留与其象限一致的方向去除180°模糊的方向，从而获得准确的海面风向。

3、将极坐标下的海面静态特征图像插值为笛卡尔坐标下海面静态特征图像方法为：

(1)在极坐标下的海面静态特征图像中选取扇形区域；

(2)得到扇形区域中像元点的极坐标(r,θ)对应的笛卡尔坐标(x,y)：

(3)建立与扇形区域相对应的笛卡尔坐标下的矩形区域，矩形区域的笛卡尔坐标为(x_i,y_i)，找到与矩形区域的笛卡尔坐标(x_i,y_i)距离最近的扇形区域的笛卡尔坐标(x,y)，将扇形区域的笛卡尔坐标(x,y)对应的极坐标(r,θ)的像元灰度值赋给矩形区域的笛卡尔坐标(x_i,y_i)，得到笛卡尔坐标下海面静态特征图像。

本发明的有益效果为：

1、本发明提出基于自适应缩减算子的导航雷达图像反演海面风向方法，设计了一种自适应缩减算子，将图像分辨率可以缩减为风条纹尺度的任意比例，扩大了缩减分辨率适应风条纹尺度的范围，提高了导航雷达图像反演海面风向的适用性；

2、本发明提出基于自适应缩减算子的导航雷达图像反演海面风向方法，设计了一种自适应算法，可以根据缩减后图像计算得到的判定参数判断缩减分辨率与风条纹尺度是否适宜，提高了海面风向反演精度；

3、本发明提出应用稳定系数作为自适应判断依据，稳定系数可以准确判断梯度方向的分布情况是否达到风向反演要求，增强了算法在工程上的适用性；

4、本发明提出基于自适应缩减算子的导航雷达图像反演海面风向方法，增添稳定系数阈值判断，提高了算法在工程中的运行速率；

5、本发明提出基于自适应缩减算子的导航雷达图像反演海面风向方法，稳定系数阈值是通过大量统计实际导航雷达数据得到的，增强了算法的可靠性。

附图说明

图1 具体实施方式流程图；

图2 局部梯度法反演海面风向的流程图；

图3a 导航雷达图像中值滤波前图像；图3b导航雷达图像中值滤波后图像；

图4 方位向归一化导航雷达图像序列；

图5 海面静态特征图像；

图6 最近点插值示意图；

图7 笛卡尔坐标海面静态特征图像；

图8 部分笛卡尔坐标下的海面静态特征图像平滑过程示意图；

图9 一次平滑图像；

图10 F(x_i,y_j)的部分一次平滑图像灰度值；

图11 缩减后图像；

图12 部分缩减后图像平滑过程；

图13 二次平滑图像；

图14 亚采样率K＝2梯度方向概率分布图；

图15 最优亚采样率梯度方向概率分布图；

图16 导航雷达设备及探测具体参数表；

图17 传统风向与参考风向误差分布结果；

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

图19 海面风向误差统计表。

具体实施方式

下面将结合附图对本发明做进一步详细说明。

本发明一种基于自适应缩减算子的导航雷达图像反演海面风向方法，包括以下几个步骤：

步骤1，导航雷达图像序列采集。采集一组导航雷达图像序列，导航雷达图像序列中包含N幅导航雷达图像，同时收集对应时间和位置的实际风向与风速值。

步骤2，导航雷达图像预处理。对导航雷达图像序列中导航雷达图像进行中值滤波，抑制同频信号对海面风向反演的干扰。

步骤3，海面静态特征图像提取。对滤波后极坐标导航雷达图像进行方位向归一化，使导航雷达图像序列中的每幅导航雷达图像的方位向线数固定。再对归一化的导航图像序列通过建立全局低通滤波器，提取出包含风条纹的海面静态特征图像。

步骤4，海面风向反演方法。首先将极坐标下的海面静态特征图像插值到笛卡尔坐标下，依据稳定系数确定自适应缩减算子，从而得到最优梯度直方图统计图，应用主梯度方向和风条纹特征得到海面风向。

本发明具体实施方式流程图见图1，分为导航雷达序列采集、导航雷达图像预处理、海面静态特征图像提取和海面风向反演这四大块。具体实施步骤共分为十三步，第一步为导航图像序列采集；第二步是导航雷达图像预处理；第三步到第四步为海面静态特征图像提取；第五步到第十三步为海面风向信息反演。具体步骤如下：

第一步，以2010年10月22日10：35的导航雷达图像序列为例，采集此刻导航雷达图像序列，其中包含位置相同、时间连续的32幅导航雷达图像，经历时间总长度为T(约1.5分钟)，同步记录由风向标和风速计获得的风向θ_w、风速U_w，此刻风向θ_w＝36^°，风速为U_w＝17.1m/s。

第二步，通过导航雷达图像序列预处理抑制同频信号对海面风向反演的影响，对导航雷达图像序列中每幅导航雷达图像都应用3×3模板的2D非线性平滑中值滤波，滤波后图像灰度值g'(r,θ)为：

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

将中值滤波器3×3模板中心与极坐标图像的N(r,θ)中心重合，通过(r,θ)与周围8个相邻像元点(s,t)的回波强度值比较，选取回波强度中间值来更新N(r,θ)处的回波强度值，模板以步长单位1遍历极坐标导航雷达图像，最终获得中值滤波后的导航雷达图像，图3为本发明中值滤波前后导航雷达图像。

第三步，由于雷达本身和外界环境的干扰，导致导航雷达图像序列中每幅导航雷达图像方位向线数都不一致。为了使图像序列的方位向线保持一致，本发明提出应用方位向归一化来解决，具体步骤如下：

①读取导航雷达图像方位向线数和径向点数，方位向大约3600条射线状线，即大约N＝3600个角度值，线与线间隔角度大约为0.1°Ω_i,i＝1,2,…N，径向选取220个点；

②建立方位向固定1800条线新极坐标图像，即N_new＝1800个角度值，θ_j,j＝1,2,…N_new，径向220个点，方位向固定间隔0.2°；

③为新极坐标图像赋灰度值，若Ω_i＝θ_j，或第一个Ω_i>θ_j，则将Ω_i线所对应的灰度值赋到新图像θ_j所对应的线上；

④重复步骤③直到新图像上的N_new条线上都具有中值滤波后g'(θ,r)灰度值，得到新雷达图像灰度值为f'(θ,r)，从而获得方位向归一化的新极坐标导航图像序列f'(θ,r,t)。图4为方位向归一化导航雷达图像序列，由于岸基、遮挡对图像的影响，方位向选取106°～291°，径向选取220个点，即600m～2250m，导航雷达图像中白色区域为去除的区域。

第四步，为获得风条纹图像，先要从导航雷达图像序列中提取包含风条纹特征的海面静态特征图像。本发明通过构建的全局低通滤波器来实现，全局低通滤波器的构建如下：

其中，f'(θ,r,t)为方位向归一化后导航雷达序列，f(θ,r)为极坐标海面静态特征图像，N_t为导航图像序列中包含导航雷达图像的个数，N_t＝32。

图5为海面静态特征图像，从海面静态特征图像中可以观测到明暗相见的条纹特征，即为风条纹。图中观测可以得到条纹尺度在300m左右，风条纹条纹方向几乎与海面风向平行，这为海面风向反演奠定了基础。

第五步，为了使海面静态特征图像便于空间运算，在保证不破坏其空间特征的前提下，选取适当极坐标下的海面静态特征图像插值为笛卡尔坐标下的海面静态特征图像，具体步骤如下：

①在极坐标海面静态特征图像中选取相对正北向203°、距离雷达630m、198*198个像元点，由于雷达图像分辨率为7.5m，即选取区域为近似1485m*1485m的矩形区域，如图4中黑色框区域所示。

②若近似矩形区域的海面静态特征图像像元点极坐标表示为(r,θ)，插值为笛卡儿坐标为(x,y)，则根据三角定理两者存在下式关系：

③设建立的矩形区域的笛卡尔坐标为(x_i,y_i)，找到与(x_i,y_i)距离最近的(x,y)点，利用最近点插值将(x,y)对应的极坐标(r_i,θ_i)像元灰度值，赋给新建的矩形区域的(x_i,y_i)点上。图6为最近点插值示意图。

④若取新建的矩形笛卡尔坐标为(x₀,y₀)，则求得离其最近点的图像极坐标(r₀,θ₀)为：

其中，rem()为求余函数，round()为向最近点取整函数，这样可将插值结果转为选取中心角度为正北向的笛卡尔坐标。图7为插值后得到的笛卡尔坐标下海面静态特征图像。

第六步，为了得到纯净的风条纹信息，本发明先要对海面静态特征图像进行平滑，去除高频信号对图像的干扰。

设笛卡尔坐标下的二维海面静态特征图像为f(x_i,y_j)(i＝1,2,…N_x,j＝1,2,…N_y)，其中(x_i,y_j)为沿x和y轴的笛卡尔坐标，N_x和N_y为笛卡尔坐标下沿x和y轴所取像元总数。对f(x_i,y_j)应用4阶杨辉三角滤波二项式卷积核，滤除高频细节信号。图8为部分笛卡尔坐标下的海面静态特征图像，则平滑过程可以简述为：

F(x_i,y_j)为平滑后像元灰度值，(x_i,y_j)为上图虚线相交中心坐标位置；4阶二项式卷积核为：

将上述步骤以步长单位1遍历整幅笛卡尔坐标下的海面静态特征图像，得到一次平滑图像F(x_i,y_j)(i＝1,2,…N_x-4,j＝1,2,…N_y-4)，如图9。图9中图像像元点数为194*194，相对图7高频信号被滤除了，从而平滑了图像中部分细节图像，但静态的风条纹信号并没有破坏。

第七步，对一次平滑图像F(x_i,y_j)(i＝1,2,…N_x-2,j＝1,2,…N_y-2)采用金字塔式降采样方法来增大采样间隔，降低图像分辨率，使风条纹信息更明显。本发明金字塔式降采样方法应用自适应缩减算子实现，自适应缩减算子如下：

C_(↓K)表示缩减率为K的缩减算子,K＝2,3,…8，K的最大值选取为8，是根据导航雷达分辨率而设定的，本发明使用雷达分辨率为7.5m，风条纹的Nyquist频率对应的波长为100m，K＝8缩减后图像分辨率变为60m，满足缩减后的图像分辨率不大于100m，这样保证缩减后在图像中可以保留完整的风条纹信息；C_(↓K)矩阵为由单位1组成的K×K矩阵。

缩减一次平滑图像F(x_i,y_j)过程中，经过C_(↓K)缩减后的图像F^(K)可以表示为：

F^(K)＝C_(↓K)*F

则缩减过程可以表达为：

式中F^(K)(x_α,y_β)为缩减后的像元灰度值，(x_α,y_β)为缩减后图像新生成的坐标，α＝1,2,…N_x-4/2，β＝1,2,…N_y-4/2，缩减后图像分辨率变为K*γ，γ为图像缩减前图像分辨率。第一次缩减运算时，令缩减率K选取为2，图10为F(x_i,y_j)的部分一次平滑图像灰度值，缩减算子为则此图像缩减后像元灰度值为：

此时由于C_(↓2)为二阶算子，(x_α,y_β)与(x_i,y_j)相同如图中虚线相交位置所示。将算子以步长K沿x和y轴遍历整个一次平滑图像F(x_i,y_j)，最终得到缩减后像元灰度值F⁽²⁾(x_α,y_β)(α＝1,2,…N_x-4/2,β＝1,2,…N_y-4/2)，缩减后图像分辨率变为15m，如图11。图11中图像像元数为97*97，增大了采样间隔，图像中高频信号被进一步滤除，相对于其他信号信息风条纹信号更加明显，并且整个图像中大部分都是风条纹图像，为利用风条纹反演海面风向提供了更好的风条纹信息。

第八步，为了得到纯净的风条纹信息，本发明对缩减后图像进行平滑，去除缩减后图像高频信号对风条纹信号的影响。

应用与一次平滑相同的4阶杨辉三角滤波二项式卷积核，对缩减后图像F⁽²⁾(x_α,y_β)(α＝1,2,…N_x-4/2,β＝1,2,…N_y-4/2)进行平滑，去除高频细节信号。图12为部分F⁽²⁾(x_α,y_β)图像，则平滑过程可以简述为：

G(x_α,y_β)为二次平滑像元灰度值，(x_α,y_β)为上图虚线相交中心坐标位置，H⁴如公式(6)所示。

将上述步骤以步长单位1遍历整幅缩减图像F⁽²⁾(x_α,y_β)(α＝1,2,…N_x-4/2,β＝1,2,…N_y-4/2)，得到二次平滑图像如图13。图13图像像元数为93*93，对图11进一步平滑去掉了更多细节信号，使得风条纹特征更加明显。

第九步，经过以上缩减平滑后得到的二次平滑图像几乎都是风条纹信号，由于风条纹方向与风向平行，则所有像元梯度的主方向与风向垂直。为了获得像元的梯度方向，这里应用优化的Sobel梯度算子。

对二次平滑图像G(x_α,y_β)沿x和y轴应用优化Sobel算子计算图像的梯度方向，优化Sobel算子为：

从而得到二次平滑图像每个像元点沿x和y轴的梯度值：

G_x和G_y分别为所有像元点沿x轴和y轴的梯度方向，由G_x和G_y得到二次平滑图像像元每点的梯度方向G_θ为：

第十步，为了判定应用此缩减算子缩减后图像是否为纯净的风条纹图像，本发明应用梯度方向计算得到的稳定性系数来判定。

为了计算稳定系数，先要对计算得到的二次平滑像元梯度方向进行直方图统计，得到梯度方向概率分布图，如图14。图14为K＝2时梯度方向概率分布图。直接获得频率最大的梯度方向选取从的所有梯度方向，其中选取λ＝70，组成的方向值数组G_P＝[G₁,G₂,…G_p]，即包含的梯度方向数量为p＝41。稳定系数η计算公式为：

其中，σ为G_P的标准差，μ为G_P的均值。σ和μ计算公式为：

最终得到缩减率K＝2的稳定系数η＝13.4。

第十一步，为了提高计算效率，通过选取阈值提前判断计算得到的稳定系数是否达到稳定标准，即判定自适应缩减算子是否与风条纹尺度匹配。为了得到适当的阈值，通过大量的统计稳定系数与不同缩减率K缩减后得到海面风向精度的关系，最终，选取阈值为10，自适应循环条件为η<10。

在K<8时，当η<10时，执行第十二步；当η>10时，则K＝K+1，从第六步开始重新执行。运行完K＝8，并没有发现存在满足条件η<10时，则选取η最小时对应缩减率K，即：K＝K(min(η))，得到的梯度方向，执行第十二步。

第十二步，通过判断得到稳定系数满足条件的缩减率K，对第九步得到的满足条件K对应的梯度方向G_θ进行直方图统计，得到梯度方向概率分布图，如图15。图15为最优K得到的梯度方向概率分布图，相比图14像元总数比例提到0.025，说明梯度方向在此缩减率K下分布更集中，说明获得了纯净的风条纹，有利于计算出精确的海面风向。

从图15中获得的概率最大的梯度方向为选取从的所有梯度方向，组成的方向值数组G_M＝[G₁,G₂,…G_m]，即包含的梯度方向数量为m。对G_M求平均得到主梯度方向G_main：

得到主梯度方向G_main＝98.4°。

第十三步，由于海面风向与风条纹方向平行，而与风条纹的梯度方向垂直，则可以根据主梯度方向G_main得到的海面风向W_d：

W_d＝G_main±90°

计算出的海面风向W_d存在180°模糊问题，为了解决这个问题将计算出来的方向与风向标测得的风向所在象限比较，保留与其象限一致的方向去除180°模糊的方向，从而获得准确的海面风向，W_d＝8.4°。

本发明是在极坐标下选取的笛卡尔坐标，由于选取区域位置和船艏向的影响，计算的海面风向需要校正后才能得到相对正北方向的海面风向，校正公式为：

N_w＝|θ_c|+|α|+|W_d|

其中，N_w为相对正北向的海面风向；θ_c为选取区域的中心角115°；W_d为笛卡尔坐标下计算得到的风向；α为船艏向93°。将N_w换算为180°象限以内，最终得到海面主风向为36.4°，与风向标测得的36°只相差0.4°，算法完全符合工程要求，传统梯度法反演得到海面主方向为32.5°，相对传统局部梯度法精度提高了88.5％。下面通过大量导航雷达序列来验证算法的有效性和适用性。

本发明的基于自适应缩减算子的导航雷达图像反演海面风向方法是在福建平潭县海坛岛开展实验的，雷达实验平台及风向标和风速计安装在海坛岛岸基，探测海域平均海深25m，受地形影响此海域经常出现大于4级的高海况。导航雷达设备及探测具体参数见图16。实验数据选取2010年10月22日、24日、30日共180组导航雷达序列，由于台风“鲶鱼”的影响，22日海面风向主要是东北风，24日短暂出现西南风，台风过后30日海面风向转为东北风。参考数据为风速计和风向标每分钟采集的海面风向和风速值，为了对应每组导航雷达序列所经历时间，将由风速计和风向标采集的值进行两分钟平均，得到‘参考风向’和‘参考风速’。

将由自适应缩减算子反演的海面风向以下称“本发明风向”，由传统局部梯度法反演的海面风向以下称“传统风向”。为验证本发明的精度，分别统计本发明风向与传统风向和参考风向的误差分布，见图17、18。通过对两个误差图对比发现本发明风向与参考风向的误差波动在20°以内，传统风向与参考风向的误差波动在40°以内，说明本发明风向比传统风向更接近参考风向。两种方法海面风向误差统计如图19。通过图19可以得出本发明风向与参考风向在各项指标上都由于传统风向，通过误差统计计算反演精度提高了58.3％。

根据风条纹尺度特征为200～500m，本发明提供了一种基于自适应缩减算子的导航雷达图像反演海面风向方法。较传统方法固定缩减次数，利用自适应缩减算子通过不同缩减率K的缩减海面静态特征图像，可以根据稳定系数适应不同的条纹尺度，从海面静态特征图像中提取出纯净的风条纹图像。而且本发明的自适应缩减算子可以将图像缩减为任意分辨率，从而提高了风条纹反演风向的利用率，提高了风向反演精度。本发明反演的海面风向与参考风向的相关系数达到了0.9956，标准差为7.62°，偏差为-1.04°，完全符合工程要求。

Claims (3)

1.一种基于自适应缩减算子的导航雷达图像反演海面风向方法，其特征在于，包括以下几个步骤：

步骤一：采集N幅导航雷达图像形成一组导航雷达图像序列；

步骤二：对导航雷达图像进行中值滤波处理；

步骤三：对滤波后的导航雷达图像进行方位向归一化，固定每幅滤波后的导航雷达图像的方位向线数；

步骤四：对方位向归一化后的导航雷达图像进行全局低通滤波处理，得到包含风条纹的海面静态特征图像；

步骤五：根据海面静态特征图像，进行海面风向反演得到海面风向；

所述根据海面静态特征图像，进行海面风向反演得到海面风向的方法为：

(1)将极坐标下的海面静态特征图像插值为笛卡尔坐标下海面静态特征图像；

(2)将笛卡尔坐标下海面静态特征图像进行平滑处理，得到一次平滑图像；

笛卡尔坐标下海面静态特征图像为：

f(x_i,y_j) i＝1,2,…N_x,j＝1,2,…N_y

其中(x_i,y_j)为沿x和y轴的笛卡尔坐标，N_x和N_y为笛卡尔坐标下沿x和y轴所取像元总数，一次平滑图像为：

$F(x_i,y_j)=\underset{m=-2}{\overset{2}{\&Sigma;}}\underset{n=-2}{\overset{2}{\&Sigma;}}{H}^r(m,n)f(x_{i-m},y_{j-n}),i=1,2,...N_x-4,j=1,2,...N_y-4$

H^r(m,n)为二项式卷积核，r为二项式卷积核的阶数，(m,n)为二项式卷积核的坐标；

(3)应用自适应缩减算子对一次平滑图像进行缩减，得到缩减后图像；

缩减后图像为：

F^(K)＝C_(↓K)*F

自适应算子C_(↓K)为：

进一步得到：

$F^{\left(K\right)}(x_{\mathrm{\&alpha;}},y_{\mathrm{\&beta;}})=\frac{1}{K^2}\underset{i=1}{\overset{K}{\&Sigma;}}\underset{j=i}{\overset{K}{\&Sigma;}}F(x_i,y_j)$

其中，(x_α,y_β)为图像缩减后新生成的坐标，α＝1,2,…(N_x-4)/2，β＝1,2,…(N_y-4)/2，K为缩减率，缩减后图像分辨率变为K*γ，γ为图像缩减前图像分辨率；

(4)将缩减后图像进行平滑处理，得到二次平滑图像；

二次平滑图像为：

$G(x_{{\mathrm{\&alpha;}}^{\&prime;}},y_{{\mathrm{\&beta;}}^{\&prime;}})=\underset{M=-2}{\overset{2}{\&Sigma;}}\underset{N=-2}{\overset{2}{\&Sigma;}}{H}^R(M,N)F^{\left(K\right)}(x_{{\mathrm{\&alpha;}}^{\&prime;}-M},y_{{\mathrm{\&beta;}}^{\&prime;}-N})$

其中，(x_α′,y_β′)二次平滑图像坐标，H^R为二项式卷积核，R为二项式卷积核阶数，(M,N)为二项式卷积核坐标，

(5)对二次平滑图像应用优化Sobel梯度算子得到像元梯度方向直方图；

优化Sobel梯度算子为：

$\left\{\begin{array}{c}{D}_x=\frac{1}{32}\left[\begin{array}{ccc}3& 0& -3\\ 10& 0& -10\\ 3& 0& -3\end{array}\right]\\ D_y={D_x}^T\end{array}\right.$

其中，D_x和D_y为优化Sobel算子分别沿x，y轴的梯度算子，

二次平滑图像每个像元点沿x和y轴的梯度值为：

$\left\{\begin{array}{c}{G}_x=G*D_x\\ G_y=G*D_y\end{array}\right.$

其中，G_x和G_y分别为所有像元点沿x轴和y轴的梯度方向，由G_x和G_y得到每个像元点的梯度方向G_θ为：

$G_{\mathrm{\&theta;}}=arctan\frac{G_y}{G_x}$

对得到的所有像元的梯度方向进行直方图统计，得到像元梯度方向直方图；

(6)通过像元梯度方向直方图得到稳定系数，基于自适应算法确定最优缩减率K；

从像元梯度方向直方图中选取从的所有梯度方向，λ＝70，组成方向值数组G_P＝[G₁,G₂,…G_p]，

稳定系数为：

$\mathrm{\&eta;}=\frac{\mathrm{\&sigma;}}{\mathrm{\&mu;}}\&times;100\%$

其中，σ为G_P的标准差，μ为G_P的均值，

选取阈值为10，缩减率K的取值范围为：K＝2,3,…8，

在缩减率K<8，并且稳定系数小于阈值时，选取当前缩减率K作为最优缩减率

在缩减率K<8，并且稳定系数大于阈值时，K＝K+1，重复步骤(2)～步骤(6)，

在缩减率K＝8，并且稳定系数大于阈值时，选取稳定系数最小时对应的缩减率K作为最优缩减率；

(7)求得对应最优缩减率K的主梯度方向；

将当前的对应最优缩减率K的像元梯度方向直方图中选取的所有梯度方向，组成的新的方向值数组G_M＝[G₁,G₂,…G_m]，包含的梯度方向数量为m，对新的方向值数组求和得到主梯度方向G_main：

$G_{main}=\underset{i=1}{\overset{m}{\&Sigma;}}{G}_i;$

(8)根据主梯度方向求得海面风向W_d：

W_d＝G_main±90°

将计算出来的海面风向与风向标测得的风向所在象限比较，保留与其象限一致的海面风向去除180°模糊的海面风向，从而获得准确的海面风向。

2.根据权利要求1所述的一种基于自适应缩减算子的导航雷达图像反演海面风向方法，其特征在于：所述的对滤波后的导航雷达图像进行方位向归一化的方法为：

(1)读取极坐标下的每个滤波后的导航雷达图像的方位向线数和径向点数，方位向线数为3600条，间隔角度为0.1°，有N＝3600个角度值Ω_i,i＝1,2,…N，径向点数为220个；

(2)建立方位向固定为1800条的新极坐标导航雷达图像，间隔角度为0.2°，有N_new＝1800个新的角度值，θ_j,j＝1,2,…N_new，径向点数为220个；

(3)为新极坐标导航雷达图像赋灰度值，若Ω_i＝θ_j，或第一个Ω_i>θ_j，则将角度值Ω_i对应的方位向线的灰度值赋给新极坐标导航雷达图像的新角度值θ_j对应的方位向线上；

(4)重复步骤(3)直到所有新极坐标导航雷达图像的所有方位向线都具有灰度值，得到方位向归一化后的导航雷达图像。

3.根据权利要求1所述的一种基于自适应缩减算子的导航雷达图像反演海面风向方法，其特征在于：所述的将极坐标下的海面静态特征图像插值为笛卡尔坐标下海面静态特征图像方法为：

(1)在极坐标下的海面静态特征图像中选取扇形区域；

(2)得到扇形区域中像元点的极坐标(r,θ)对应的笛卡尔坐标(x,y)：

$\left\{\begin{array}{c}x=r*cos\mathrm{\&theta;}\\ y=r*sin\mathrm{\&theta;}\end{array}\right.;$

(3)建立与扇形区域相对应的笛卡尔坐标下的矩形区域，矩形区域的笛卡尔坐标为(x_i,y_i)，找到与矩形区域的笛卡尔坐标(x_i,y_i)距离最近的扇形区域的笛卡尔坐标(x,y)，将扇形区域的笛卡尔坐标(x,y)对应的极坐标(r,θ)的像元灰度值赋给矩形区域的笛卡尔坐标(x_i,y_i)，得到笛卡尔坐标下海面静态特征图像。