基于X波段航海雷达的海面风场反演技术研究

王海霞，徐 进，王庆名，陈 丽，庞玺斌

(1.大连海事大学 航海学院，辽宁 大连 116026；2.大连海事大学 校友事务与合作处，辽宁 大连 116026)

0 引言

海洋风场是海洋表层运动的主要推力，它管控能源和动力的转移，能够将空气带到海洋，在气候变化中占有非常重要的作用。海面风场常被应用于各类海洋科学研究，如海洋环流、热带风暴、大气与海洋的耦合效应等[1]。因此，监测海面风场变化有着极其重要的意义。

目前，国外在航海雷达海上风场的反演方面已经存在一些研究成果。Ziemer等人[2]首次应用二维傅里叶变换，分析2张带有海杂波的雷达图像，计算的海浪谱与传统浮标获取的海浪谱值较为接近。Alpers等[3]和Vesecky等[4]对雷达海浪成像作了详细描述，但他们均未考虑流的影响。Lee等人[5]基于电磁波海面散射机理，指出海面的粗糙度与风速之间，呈现正相关特性。Hatten等人[6]发现航海雷达散射截面风向及电磁波照射方向存在一定关系，运用这种关系可以判断出风向，但此方法存在一定的局限性，往往不能满足所需的全方位航海雷达散射截面。Dankert等人[7]提出了一种基于光流运动估计的技术，从X波段航海雷达图像中获得风速和风向，该方法相对简单，不需对雷达系统的相位进行校正，也不存在角度模糊问题，所获得的海面风场具有较高的时空分辨率。Liu等人[8]针对未被遮挡的低散射强度雷达图像，提出了一种散射强度值选取法，并用对偶曲线拟合法减小了风向反演误差。

国内应用雷达监测海面风场研究的起步较晚，杨劲松[9]用合成孔径雷达(SAR)，采用二维傅里叶变换波数谱的方法，反演出了风向，并用COMD4模式函数算出了风速。金丽洁[10]利用浪高反演算法[11]，采用多重信号分类算法完成定向，实现了浪高场的反演。张璇等人[12]分析目前交叉极化对风场反演的优势，总结了多极化SAR图像反演海面风场的最新研究成果，并对多极化SAR图像反演海面风场的业务应用发展进行了初步展望。

本文基于X波段航海雷达反演海面风场的基本原理，以采集到的雷达原始图像为数据基础，基于梯度算法求得主风向，并与大连气象台的实测值进行对比，以期为海面风场的实时监测提供有效的技术手段。

1 海面风场反演方法

1.1 数据来源及设备参数

X波段船载雷达信号采集到的原始灰度图像(设备参数详见表1)监测范围为0.75海里，图像大小为1 024×1 024(相同地点连续获取的3幅图像如图1所示)。监测数据采集位置：大连港附近海域(121.830°E，38.913°N)。船载平台：大连海事大学“育鲲”轮号，采集时间为2015年8月14日20时4分。数据中含有大量的含杂波信息，用于反演海面海浪与风向等信息的研究。

表1 育鲲轮航海雷达设备参数

参数项参数值工作极化方式水平工作频段X波段雷达天线长度/ft8天线旋转速度/RPM28脉冲工作频率/Hz3 000/1 800 /785雷达控制台23Inch液晶工控机工作量程/海里1～12

图1 雷达图像原图

1.2 实验流程

海面风场反演过程包括雷达信号的实时采集和海面风场信息提取2部分。航海雷达风场反演流程如图2所示。

图2 航海雷达风场反演流程

1.3 图像预处理

1.3.1 图像求和

图像求和即将图像对应像素点的值进行求和，从而增强海浪信息所在位置的灰度强度。n幅图像中，每个像素点f(x，y)求和公式为

(1)

1.3.2 中值滤波法

中值滤波基于排序统计理论，将每一像素的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值，公式为

g(x，y)=median({f(x-i，y-j)，(i，j)∈W)}，

(2)

式中，g(x，y)为输出像素值；f(x-i，y-i)表示输入像素值；W是含有基数点的模板窗口，一般为3×3窗口或5×5窗口。

1.3.3 梯度法反演风向

风向反演主要是基于海面风引起的纹理可在雷达回波图像中成像，Dankertd等人[13]研究发现雷达图像中存在与风向平行的纹理信息，分析雷达图像中的海浪纹理可以获取风向[14]。梯度是用于描述图像灰度变化情况的二维列向量，能反映变化的剧烈程度[15]。通过将图像进行分割并统计梯度方向与大小得到该区域的梯度特征分布的统计结果[16]，梯度的角度为图像像素变化的方向，梯度的模为梯度方向的变化率。若求出海浪纹理变化最大的梯度方向，即可得到与之垂直的主风向。算法实施步骤如下。

(1)图像重采样

对雷达原始图像采用含有一定滤波作用的算子进行计算，平滑图像的同时，压缩图像像素数量。可用算子公式如下：

R=B2S2B4，

(3)

(4)

(5)

式中，S2表示2×2线性插值平均，先用四阶算子平滑，然后缩减图像，最后用二阶算子再次平滑图像。

(2)局部梯度的计算

用最优化的Sobel算子进行局部梯度计算：

(6)

Dy=DxT，

(7)

式中，T表示转置矩阵。此时图像的梯度G为

G=(Dx+iDy)*(A)，

(8)

式中，A表示图像；*表示卷积，G′与G″的计算方式如下所示：

G′=R*(G2)，

(9)

G″=R*(|G2|)。

(10)

由文献[17]可知，一致性参数C为：

(11)

(3)确定风向

应用一致性参数和梯度的模，做梯度方向角度的加权直方图。因为风向与梯度变化最大的方向垂直，在直方图中最大尖峰值加上90°，即为主风向。

2 实验结果与分析

2.1 雷达原始图像求和处理结果

在获取的雷达图像中包含许多信息，但在获取过程中由于各种原因会存在不同形式的干扰。因此在反演前应当先进行图像处理。首先将连续获取的雷达原始图像进行图像求和，增强海浪信息所在位置的灰度强度，如图3所示。

图3 雷达原始图像求和

2.2 中值滤波处理结果

采用5×5窗口，对雷达原始图像求和处理后的图像进行中值滤波处理，降低同频干扰，如图4所示。在进行了预处理后，为消除图像存在的尾迹效应，截取了没有尾迹效应的部分截取了没有尾迹效应的部分，即(275：402，445：572)中的图像窗口进行风向反演。

图4 中值滤波处理

2.3 局部梯度法反演风向

采用Matlab编写程序对图像进行梯度计算，并做出相应的梯度加权直方图，如图5所示。由图5可以看出，在梯度角度为88°时出现的频率最大。因为主风向与梯度角度直方图的峰值方向垂直，应当在获得的梯度角度加上90°，即为所求的主风向178°。

图5 梯度加权直方图

3 验证

根据大连气象台提供的2015年8月14日平均风向为183°，与实测结果相比风向偏差5°，反演的风向与预报结果是比较吻合的。因本文所用的雷达数据极化方式为水平极化，可以很好地消除180°模糊问题，因此，反演结果不存在180°模糊。

4 结束语

以大连海事大学教学实习船“育鲲”轮采集的航海雷达原始数据为基础，基于梯度算法进行风场反演。在实验过程中，存在多种因素可能会影响计算的结果。首先，在雷达成像过程中，获取到的风场信息包含了雷达重力波和毛细波的影响。此外，在处理同频干扰和尾迹效应问题时，虽然采用了中值滤波，但干扰并未完全消除。本文提出的方法还未考虑其他因素，如没有考虑到船舶运动对于雷达回波的影响；遇到遮挡物时，在遮挡物后方向形成盲区，不能准确显示出该区域的真实海面状况。这些对于风场反演的结果都将存在影响，使结果产生误差。而且，计算结果为风场的主风向，并不是该区域内的所有风向，各区域点之间也会存在差异。上述难点都需要在未来的工作中，继续研究。