一种基于新型海浪色散关系带通滤波器进行X波段导航雷达反演海浪参数方法

摘要

本发明属于海浪遥感技术领域，具体涉及一种利用导航雷达获取的海杂波图像进行海浪参数反演的基于新型海浪色散关系带通滤波器进行X波段导航雷达反演海浪参数方法。本发明包括：雷达图像数据采集；雷达图像预处理；对笛卡尔坐标系下的图像序列应用傅里叶变换，得到雷达图像的三维波数频率图像谱；海浪谱信息提取；海浪信息反演。本发明的滤波器保留了船速对色散关系的影响，有效地解决了传统宽带通滤波器的带宽会随运动速度的增大而增大这一问题，使得可以在雷达平台随舰船运动情况下进行滤波；本发明的新型滤波器的带通边界推导中没有对任何量取近似值，减小了计算误差，不会对带通边界产生影响，带宽计算更加精确，提升了海浪反演精度。

说明书

技术领域

本发明属于海浪遥感技术领域，具体涉及一种利用导航雷达获取的海杂波图像进行海浪参数反演的基于新型海浪色散关系带通滤波器进行X波段导航雷达反演海浪参数方法。

背景技术

对X波段雷达的回波图像进行分析，可以获得海浪的参数。但是该回波图像中除了含有海浪信号外还含有大量的噪声，如同频干扰、背景噪声、雨雪干扰及目标物回波等。这些噪声将严重影响海浪的参数反演，使结果失真。所以需要从雷达的回波图像中去除这些噪声信号并提取出有用的海浪信号。由于该信号符合重力波的色散关系，所以可以利用色散关系来构造滤波器，将信号提取出来。具体的滤波流程如下：

海杂波图像谱I(k_x,k_y,ω)→色散关系滤波→信号图像谱F(k_x,k_y,ω)

海浪符合重力波的色散关系这一理论是由LeBlond和Mysak于1978年提出的。见参考文献[1](Paul H.LeBlond,Lawrence A.Mysak.Waves in the ocean[M].1978,Elsevier.)

假定在选定的分析区域内海面浪场和流场是空间均匀和时间稳定的，则一阶重力波的色散关系描述如下

$>{\omega }_0^2=g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)---\left(1\right)$>

其中，d为水深，g为当地重力加速度，为波数，λ为海浪波长。若存在一个相对于雷达天线平台速度的表面流其中那么引入一个多普勒频移的色散关系变为：

$>\omega ={\omega }_0+\overrightarrow{k}·\overrightarrow{u}---\left(2\right)$>

则

$>\omega =\sqrt{g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)}+\overrightarrow{k}·\overrightarrow{u}---\left(3\right)$>

上述为雷达与海表面不存在相对流速与存在相对流速下的色散关系曲线。曲线a为流速为0时的色散关系曲线，它关于纵轴对称。曲线b为存在多普勒频移的色散关系曲线。等高线代表能量的分布。

根据文献显示，目前已有针对海杂波图像海浪信号提取的基于色散关系的带通滤波器主要分为两大类，即窄带和宽带。在带通滤波器的发展史上，最早出现的是窄带形式的滤波器，随后为了解决窄带滤波器的不足而出现了宽带滤波器。

窄带滤波器的原理是寻找合适的方法计算出色散关系公式(3)中的流速矢量这样色散关系即为已知。由于海浪能量绝大部分分布在色散关系曲线附近，因此通过在色散关系曲线附近设置窄带通，便可以获取真实的海浪信息，滤除噪声干扰。

最早的窄带通色散关系滤波器是1999年Reichert、Nieto等人提出的。Reichert等人分别在文献[2-4]中给出了两种思想相似的带通滤波器定义，如公式(4)(5)所示。

由于海杂波图像谱由三维离散傅里叶变换得到，因此波数及频率ω均是离散化后的点。式(4)中，I(k_x,k_y,ω)为海杂波图像谱；F(k_x,k_y,ω)为带通滤波后的信号图像谱，该滤波器对所有的波数进行遍历得到对应的频率ω，并只保留该处的能量作为海浪能量，其他部分作为噪声滤除[2][3]。式(5)中d为水深，由估流得到，ω_p代表某一频率片的频率。该滤波器同样对所有的波数进行遍历获取对应的频率ω，若ω∈[ω_p-1,ω_p+1]，则保留频率片中第p、p-1或p、p+1片上该波数下的能量[4]。因此，该滤波器带通比(4)式所述带通略宽。见参考文献[2-4](Jose>

2012年沈继红等对上述滤波器进行了改进，考虑了高阶色散关系后认为海浪能量不仅仅分布在一阶色散关系曲线上。高阶色散关系公式为

$>{\omega }_m=m\sqrt{\frac{g\left|\overrightarrow{k}\right|}{m}·\tanh \left(\frac{d\left|\overrightarrow{k}\right|}{m}\right)}+\overrightarrow{k}·\overrightarrow{u}m=\mathrm{1,2},...\left(6\right)$>

这里m为色散关系的阶次，ω_m为m阶次波的频率。可见当m＝1时公式(6)即为一阶色散关系方程(3)，且高阶色散关系求取的ω_m会增大。利用高阶色散关系边界生成新的滤波器，该滤波器定义如下

其中，为滤波器上界，由高阶色散关系曲线求极限得到，比公式(4)(5)中的结果增大；为滤波器下界，与公式(4)(5)中的结果相同。该滤波器仍属于窄带滤波器，需要在带通滤波之前估计流速流向，且同样是对每个离散化的波数进行遍历。见参考文献[5](沈继红，李英，王淑娟，戴运桃，王晓莺.X波段雷达中基于高阶色散关系的滤波器设计[J].中国航海,2012，4(35):4-7.)

以上为目前已有文献中窄带滤波器的发展现状，它的滤波效果完全依赖于估流精度。由于图像分辨率有限、谱能量分布不均匀、噪声特性复杂等因素影响，目前国内外尚缺乏精确估流的技术手段，估流效果都不够理想，流速流向估计误差严重影响了窄带滤波器的性能，基于窄带通的不够真实的色散关系曲线会导致滤波器保留了大量不属于海浪的噪声而滤掉了真实海浪信息。

针对这一问题，Rune于2002年首先设计了基于最大流速U_max的宽带通滤波器[6]。该方法不需要精确估流，解决了窄带滤波器最大的不足。目前该方法已经成为被广泛使用的主流滤波器。考虑在深水的情况下，该滤波法被定义如下：

其中

$>\left|B_n\right|\frac{{(\omega -\frac{\mathrm{\Delta \omega }}{2}-U_{\max }\frac{{\omega }^2}{g})}^2}{g}-\frac{\mathrm{\Delta k}}{\sqrt{2}}$>$>\left|B_P\right|\frac{{(\omega -\frac{\mathrm{\Delta \omega }}{2}-U_{\max }\frac{{\omega }^2}{g})}^2}{g}+\frac{\mathrm{\Delta k}}{\sqrt{2}}$>

Δω为频率分辨率，Δk为波数分辨率，U_max为正常情况下海表面可能达到的最大流速。该滤波器边界公式B_n、B_P由带流速的基本色散关系公式(3)推导得到。在推导过程中对部分项取了近似值，首先当水深大于等于二分之一波长时，可以设其次由于的值小于1，的值大于值，且海浪波长越大，二者差距越大。因此为简化计算的同时尽可能地减少误差，可将带入公式(3)，那么得到

$>\omega =\sqrt{g\left|\overrightarrow{k}\right|}+\frac{{\omega }^2}{g}·\left|\overrightarrow{u}\right|·\cos \theta ---\left(9\right)$>

因此当cosθ＝1时，

$>{\left|\overrightarrow{k}\right|}_{\min }=\frac{{(\omega -\frac{{\omega }^2}{g}·|\overrightarrow{u}\left|\right)}^2}{g}---\left(10\right)$>

当cosθ＝-1时，

$>{\left|\overrightarrow{k}\right|}_{\max }=\frac{{(\omega +\frac{{\omega }^2}{g}·|\overrightarrow{u}\left|\right)}^2}{g}---\left(11\right)$>

它们分别对应了滤波器的下边界和上边界。同时考虑到频率分辨率和波数分辨率的影响便可以得到B_n、B_p滤波器的边界。见参考文献[6](Rune>

2012年，贾瑞才等人又对上述基于色散关系的宽带通滤波方法进行了改进[7]。该方法以假设的最大波数k_m区别于上述滤波方式。该滤波器定义为

其中，

$>\left|k_P\right|\frac{{(\omega -\frac{\mathrm{\Delta \omega }}{2}-U_{\max }|k_m\left|\right)}^2}{g}+\frac{\mathrm{\Delta k}}{\sqrt{2}}$>$>\left|k_n\right|\frac{{(\omega -\frac{\mathrm{\Delta \omega }}{2}-U_{\max }|k_m\left|\right)}^2}{g}-\frac{\mathrm{\Delta k}}{\sqrt{2}}$>

k_m的计算方法是首先计算出谱I(k_x,k_y,ω)中的最大值C，然后获取谱值大于0.9C的所有点的波数模值，再将波数模值算术平均，即得到|km|。该滤波器同B_n、B_p滤波器相似，用主导波数k_m代替了ω²/g，使得滤波器带宽放大，但并不如ω²/g精确。见参考文献[7](袁赣南，贾瑞才，王淑娟，戴运桃.色散关系带通滤波器设计[J].华中科技大学学报(自然科学版)，2012，3(40):36-30.)

现有的这两种宽带滤波器均通过对频率ω进行遍历，求取该ω对应的波数模值的上界和下届后，对每个ω下的波数进行筛选。这种对ω进行遍历的滤波器与对波数k进行遍历的滤波器具有相同的时间复杂度，均为0(n²)。但是，从上面的分析可见，这两种滤波器中带入的近似值将导致该滤波器存在两项缺陷：一是仅适用于水深大于等于二分之一波长情况；二是仅在大波长下该滤波器滤波效果较为精确。

除此之外，两种宽带滤波方法均无法直接在平台运动状态下使用，从算法定义中可见当雷达平台处于运动状态时，滤波器的带宽将被大幅度的放大，起不到滤波效果，使得反演结果中包含了大量的噪声干扰，影响精度。以B_n、B_p滤波器为例，当雷达平台静止时，U_max＝3m/s，随着船速的增大，U_max也应越大，带通随之变宽，且上边界的变化速度大于下边界。另外，该滤波器中的U_max也将导致运动状态下带宽增加，不能起到滤波效果。

由此可见，虽然宽带滤波器不需要估流，使用最大期望流速U_max解决了窄带通下估流不准确带来的问题，但是在动态下滤波器的带通会被无限放大，无法起到滤波效果。

针对这一问题，本发明设计了一种新型的色散关系带通滤波器。该滤波器对波数进行遍历，并且属于宽带的范畴，但是又有效地解决了传统宽带通滤波器的带宽会随运动速度的增大而增大这一问题。将雷达平台的运动速度考虑到带通的计算方程中，使运动状态下的带宽与静态相同。运动状态下船速的增大会导致海浪能量跟随色散关系产生偏移，而新带通中添加了这种偏移量，使得新带通具有跟随能量分布的特点。并且该滤波器继承了宽带滤波器无需事先估流的特点。通过实验数据证明该新型色散关系带通滤波器适用于动态环境下海杂波图像的滤波。

发明内容

本发明的目的在于提供一种提升海浪反演精度的基于新型海浪色散关系带通滤波器进行X波段导航雷达反演海浪参数方法。

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

(1)雷达图像数据采集：采集N幅空间海域杂波连续图像，并同步记录当时船舶的航向和航速信息；

(2)雷达图像预处理：选取距离船艏向角度为，像素为N1*N1的矩形区域作为分析的子图像序列，对矩形区域进行中值滤波，使用最近点插值的方法获取图像笛卡尔坐标；

(3)对笛卡尔坐标系下的图像序列应用傅里叶变换，得到雷达图像的三维波数频率图像谱；

(4)海浪谱信息提取：使用色散关系带通滤波器对三维波数频率图像谱进行滤波处理，得到三维海浪图像谱；对该图像谱中的频率域进行积分，得到二维海浪波数图像谱；使用调制传递函数MTF获取二维海浪谱。

(5)海浪信息反演：利用海浪谱反演计算得到海浪参数，包括有效波高、波峰周期、波峰峰向。

所述步骤(2)包括：

(2.1)选取距离船艏向角度为像素为N1*N1的矩形区域作为分析的子图像序列；

(2.2)对N幅图像的分析区域按照如下方式进行旋转：

其中，为第一幅图像的框心角；为N幅图像对应的航向与第一幅航向的夹角；为第i幅图像旋转后的框心角；

(2.3)对步骤(2.1)选取的雷达图像分析序列应用3×3模板的2-D非线性平滑中值滤波；

(2.4)对于雷达原始坐标系中的每一点，利用最近点插值方法获取图像的笛卡尔坐标。

所述步骤(4)包括：

(4.1)将步骤(1)中记录的航速数据转化到分析区域所在的笛卡尔坐标系下：

其中，为选取的框心角，u_船为船速；

(4.2)对步骤(3)中求取的三维波数频率能量谱I(k_x,k_y,ω)进行滤波操作

其中，

$>{\omega }_n=\sqrt{g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)}+k_x·u_x+k_y·u_y-\left|\overrightarrow{k}\right|·U_{\max };$>

$>{\omega }_p=\sqrt{g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)}+k_x·u_x+k_y·u_y-\left|\overrightarrow{k}\right|·U_{\max };$>

ω_n、ω_p为滤波器频带的上、下边界；U_max是最大水流速度；是已知量；

ω_n、ω_p由以下计算过程得出：

因此色散关系

海表面流速最大可以取值为U_max，滤波器的带通边界为

由于则

$>{\omega }_n=\sqrt{g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)}+k_x·u_x+k_y·u_y-\left|\overrightarrow{k}\right|·U_{\max };$>

$>{\omega }_p=\sqrt{g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)}+k_x·u_x+k_y·u_y-\left|\overrightarrow{k}\right|·U_{\max }.$>

本发明的有益效果在于：由于目前国内外尚缺乏精确估流的技术手段，估流效果都不够理想，流速流向估计误差严重影响了窄带滤波器的性能，工程上无法使用。本发明的滤波器继承了宽带滤波器无需估流的优点，使得滤波效果不受估流精度影响，具有很高的工程使用价值和应用前景；当雷达平台处于运动状态时，传统宽带滤波方法均无法直接在平台运动状态下使用，其滤波器的带宽随着船速增加被迅速放大，在高航速下算法完全失效，使得反演结果中包含了大量的随机背景噪声干扰。本发明的滤波器保留了船速对色散关系的影响，有效地解决了传统宽带通滤波器的带宽会随运动速度的增大而增大这一问题，使得可以在雷达平台随舰船运动情况下进行滤波；本发明的新型滤波器的带通边界推导中没有对任何量取近似值，减小了计算误差，不会对带通边界产生影响，带宽计算更加精确，提升了海浪反演精度。

附图说明

图1流速为零与流速非零情况下的色散关系曲线图以及其能量分布；

图2运动状态下B_n、B_p滤波效果不明显；

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

图4分析区域位置示意图；

图5最近点插值示意图；

图6船速坐标系转换示意图；

图7航行状态下本发明反演海表面有效波高与WAMOS有效波高对比；

图8航行状态下本发明反演海表面波峰峰向与WAVEX波峰峰向对比；

图9航行状态下本发明反演海表面波周期与WAVEX波周期对比；

图10频率为0.74rad/s时的海杂波图像谱；

图11经B_n、B_P滤波后的海浪图像谱(U_max＝10m/s)；

图12经ω_n、ω_P滤波后的海浪图像谱(U_max＝3m/s)。

具体实施方式

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

本发明公开了一种利用新型海浪色散关系带通滤波器进行X波段导航雷达反演海浪参数的方法，反演计算的海浪参数包括有效波高、波峰周期、波峰峰向等，所述方法包含雷达图像的采集，雷达图像预处理、海杂波图像谱获取、海浪信息提取和海面信息反演五个部分。通过傅里叶变换获取海杂波图像谱，并使用新型的带通滤波器将海浪能量从海杂波图像谱中提取出来，为后续反演奠定了基础。与传统方法相比较，本发明所公开的滤波方法可以准确的提取海浪信号。利用海上航行测量数据反演得到的波峰峰向值与WAVEX系统测量值的偏差仅为7.02o，波周期值的相对误差仅为2％，说明该方法不仅客服了传统方法估流不准确的缺陷而且弥补了传统滤波器无法适用于动态环境下海浪参数反演的不足。

本发明包括：

步骤1，数据采集。采集N幅空间海域杂波连续图像，并同步记录当时船向及船速。

步骤2，图像预处理。选取距离船艏向角度为像素为N1*N1的矩形区域作为分析的子图像序列。对该分析区域进行中值滤波，并使用最近点插值的方法获取图像笛卡尔坐标。

步骤3，对笛卡尔坐标系下的图像序列应用傅里叶变换，得到三维波数频率图像谱。

步骤4，海浪信息提取。使用本发明设计的新型带通滤波器对第三步中得到的三维波数频率图像谱滤波，得到三维海浪图像谱。

步骤5，海浪信息反演。对第四步的三维海浪图像谱中的频率域进行积分，得到二维海浪图像谱。使用调制传递函数(MTF)获取海浪谱。利用海浪谱反演得到波高、波峰峰向、波周期等海浪参数。

所述步骤2包括以下步骤：

步骤2.1，选取距离船艏向角度为像素为Nx*Ny的矩形区域作为分析的子图像序列；

步骤2.2，对N幅图像的分析区域按照如下方式进行旋转：

其中，为第一幅图像的框心角；为N幅图像对应的航向与第一幅航向的夹角；为第i幅图像旋转后的框心角。

步骤2.3，对步骤2.1选取的雷达图像分析序列应用3×3模板的2-D非线性平滑中值滤波；

步骤2.4，对于雷达原始坐标系中的每一点，利用最近点插值方法获取图像的笛卡尔坐标。

所述步骤4包括以下步骤：

步骤4.1，将步骤1中记录的航速数据转化到分析区域所在的笛卡尔坐标系下。转化方式如下

其中，为选取的框心角，u_船为船速。

步骤4.2，对步骤3中求取的三维波数频率能量谱I(k_x,k_y,ω)进行滤波操作

其中，

$>{\omega }_n=\sqrt{g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)}+k_x·u_x+k_y·u_y-\left|\overrightarrow{k}\right|·U_{\max };$>

$>{\omega }_p=\sqrt{g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)}+k_x·u_x+k_y·u_y-\left|\overrightarrow{k}\right|·U_{\max };$>

ω_n、ω_p为滤波器频带的上、下边界；U_max是假设的最大水流速度，可以对其进行设定；是已知量。

所述步骤4.2中ω_n、ω_p由以下计算过程得出：

基本色散关系方程如公式(4)所示，其中因此色散关系方程变为

假设正常情况下海表面流速最大可以取值为U_max，则本发明滤波器的带通边界为

由于则

$>{\omega }_n=\sqrt{g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)}+k_x·u_x+k_y·u_y-\left|\overrightarrow{k}\right|·U_{\max }---\left(18\right)$>

$>{\omega }_p=\sqrt{g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)}+k_x·u_x+k_y·u_y-\left|\overrightarrow{k}\right|·U_{\max }---\left(19\right)$>

由所述步骤三可知，本发明选择N幅图像进行反演，对应着N/2个频率ω_i(i＝1...N)。每个频率片ω_i中有N_x*N_y个每个都可以对应的分解为唯一的k_x、k_y，即(k_x,k_y,ω_i)。也就是说，每个都有N个频率ω_i与之相对应。频率的分辨率dω定义如下：

$>\mathrm{d\omega }=\frac{2\pi }{T}$>

其中，时间尺度T为N幅图像对应的总时间。N_x*N_y为所选分析区域的像素。

本发明在所选分析区域内以波数分辨率为步长对k_x、k_y进行遍历，波数的分辨率dk_x、dk_y定义如下：

$>{\mathrm{dk}}_x=\frac{2\pi }{L_x},{\mathrm{dk}}_y=\frac{2\pi }{L_y}$>

其中，L_x、L_y选择的分析区域矩形框的边长。按照公式(18)、(19)对每个求出相应的ω_n、ω_p，若ω_n≤ω_i≤ω_p则保留该点能量，其他部分的能量置为0。

下面将结合附图对本发明提出的适用于动态环境下的新型色散关系带通滤波器作进一步的详细说明。具体实施方式流程图见图3。其中，第一步为图像信息采集；第二步是对雷达图像进行预处理；第三步通过三维傅里叶变换获取波数频率谱；第四步为利用本文发明的带通滤波器对波数频率谱进行滤波处理；第五步到第七步为从滤波后的图像谱中计算得到海浪的波高、波峰峰向和波周期等参数。

第一步，采集32幅连续的海杂波图像，记录其时间总长度为T(约1.5分钟)，同步记录航向(船艏向)、航速(相对速度)值。

第二步，对所选分析区域进行图像预处理。

(1)选取距船艏75°方向、距离雷达平台的实际距离为600米、实际空间尺寸为960m*960m的矩形区域(象素分辨率为7.5m，即128*128个象素点)构成图像分析序列。分析区域位置如图4所示。由于航行过程中船艏会有一定的摆动，导致每幅图像的选框位置略有差异。为了使选择的分析区域不变，需要将32幅图像以其中的第一幅的航向为基准进行旋转。这使得选框的中心角会在75°附近摆动。公式如下

其中，为第一幅图像的框心角，本文中设定为75°；为32幅图像对应的航向与第一幅航向的夹角；为第i幅图像旋转后的框心角。

(2)对所选雷达图像分析序列应用3×3模板的2-D非线性平滑中值滤波

$>f^{\prime }(r,\theta )=\underset{(s,t)\in (r,\theta )}{\mathrm{median}}\left\{g(s,t)\right\}---\left(21\right)$>

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

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

(3)利用最近点插值方法获取图像笛卡尔坐标。将雷达原始坐标系中的每一点以极坐标的形式表示为(r,θ,z)，转换为笛卡儿坐标后表示为(x,y,z)。有

$>\left(\begin{array}{c}x=r^{*}\cos \theta \\ y=r*\sin \theta \end{array}\right)---\left(22\right)$>

矩形框点的坐标设为(x_i,y_i)。利用最近点插值就是让矩形框中的每一点都在扇形区域中找到离它最近的一点得到其对应的极坐标(r_i,θ_i)，并将扇形上的这一点的回波强度赋值给矩形框中的点(x_i,y_i)。任取一点记为(x₀,y₀)，得到其最近点的极坐标方式为：

$>\left(\begin{array}{c}{r}_0=\mathrm{round}\left(\mathrm{sqrt}({x_0}^2+{y_0}^2)\right)\\ {\theta }_0=\mathrm{round}\left(\mathrm{rem}(\arctan (y_0,x_0)+2\pi ,2\pi )\right)\end{array},-,-,-,\left(23\right)\right)$>

其中，rem()为求余函数，round()是向最近点取整函数。将极坐标下(r₀,θ₀)点对应的回波强度值赋值给矩形框里的点(x₀,y₀)便完成了雷达图像分析区域的坐标转换。图5为最近点插值示意图。

第三步，对图像序列用用傅里叶变换，将空间时间域的图像变为波数频率域，得到三维波数频率图像谱。

对海杂波序列η(x,y,t)进行三维离散傅里叶变换，如下:

$>F(k_x,k_y,\omega )=\underset{0}{\overset{L_y}{\Sigma }}\underset{0}{\overset{L_x}{\Sigma }}\underset{0}{\overset{T}{\Sigma }}\eta (x,y,t)\exp [-2\mathrm{\pi i}(k_{x}x/L_x+k_{y}y/L_y+\mathrm{\omega t}/T)---\left(24\right)$>

其中，F(k_x,k_y,ω)为变换后得到三维波数—频率图像谱。本文选择的分析区域是边长为L_x*L_y的矩形框，其中L_x＝L_y＝960m；时间尺度T为32幅图像对应的总时间；k_x和k_y分别为x和y方向的波数，ω为频率。

三维波数频率能量谱由下面公式定义：

$>I(k_x,k_y,\omega )=\frac{1}{L_x{L}_{y}T}{\left|F(k_x,k_y,\omega )\right|}^2---\left(25\right)$>

第四步，利用本文发明的基于色散关系的带通滤波器进行滤波。

(1)读取第一步信息采集中获取的船速数据，将其转化到分析区域所在的笛卡尔坐标系下。按照流速转化示意图，图6，转化方式如下

其中，为选取的框心角带入75°，u_船为船速。由于船艏向朝向x轴正方向，则相当于水向x轴负方向流。因此图6中船速朝向x轴负方向。

(2)对第三步中求取的三位波数频率能量谱I(k_x,k_y,ω)进行滤波操作。为了将相对于雷达平台的实时海表面流速引入到色散关系带通滤波器的边界中，使滤波器的边界能随流速的改变而改变，根据海浪的色散关系设计带通滤波器

ω_n是滤波器的频率下边界，ω_p是滤波器的频率上边界，与是已知项，U_max是正常状况下水流可以达到的最大速度本文设定为3m/s。且则因此ω_n、ω_p变为

$>{\omega }_n=\sqrt{g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)}+k_x·u_x+k_y·u_y-\left|\overrightarrow{k}\right|·U_{\max };$>

$>{\omega }_p=\sqrt{g\left|\overrightarrow{k}\right|\tanh \left(\right|\overrightarrow{k}\left|d\right)}+k_x·u_x+k_y·u_y-\left|\overrightarrow{k}\right|·U_{\max };$>

本发明选择32幅图像进行反演，对应着32个频率ω_i(i＝1...32)。本文中ω_i以频率分辨率dω为步长递增，频率的分辨率dω定义如下：

$>\mathrm{d\omega }=\frac{2\pi }{T}$>

其中，时间尺度T为32幅图像对应的总时间。

由于分析区域矩形框的实际空间尺寸为960m*960m，距离分辨率本文定为7.5m，因此每个频率ω_i中有128*128个

128*128即为所选分析区域的像素。

本发明在所选分析区域内以波数分辨率为步长对k_x,k_y进行遍历，k_x,k_y的取值范围为(-0.4189，0.4123)。波数的分辨率dk_x、dk_y定义如下：

$>{\mathrm{dk}}_x=\frac{2\pi }{L_x},{\mathrm{dk}}_y=\frac{2\pi }{L_y}$>

其中，L_x＝960m、L_y＝960m为选择的分析区域矩形框的边长。按照公式(18)、(19)对每个求出相应的ω_n、ω_p，若ω_n≤ω_i≤ω_p则保留该点能量，其他部分的能量置为0。

第五步，三维波数频率能量谱I(k_x,k_y,ω)经滤波后得到海浪三维图像谱F(k_x,k_y,ω)，将其对频率ω进行积分可以得到二维图像谱F(k_x,k_y)。由于傅里叶变换关于原点的对称性，若在整个ω域积分将失去相位速度方向信息，故只对ω＞0(或ω＜0)区域进行积分，定义如下：

$>I(k_x,k_y)={\int }_{\omega >0}E(k_x,k_y,\omega )\mathrm{d\omega }---\left(28\right)$>

第六步，使用调制传递函数MTF将图像谱转换为海浪谱

由于X波段航海雷达成像机制的非线性性，二维图像谱与实际的海浪谱之间会存在一定的差别，呈现非线性关系，在波数较大的区域表现得更为明显。因此，需要通过调制传递函数MTF对图像谱进行相应的变换得到海浪谱。计算公式如下：

$>E(k_x,k_y)={\left|M(k_x,k_y)\right|}^2·F(k_x,k_y)---\left(29\right)$>

式中E(k_x,k_y)为二维海浪，|M(k_x,k_y)|²为调制传递函数。

调制传递函数目前并没有对应的理论推导形式，而是利用大量试验数据拟合得到的。下面公式给出了经验的调制传递函数定义式：

|M(k_x,k_y)|²≈k^-β(30)

其中β为经验系数，一般在1.2左右取值，已被广泛使用。

第七步，求取波高、波峰峰向和波周期。

(1)波高的反演方法

本发明使用的X波段雷达计算有效波高的方法是基于有效波高与信噪比的平方根呈线性关系的原理来实现的，即

$>H_s=A+B\sqrt{\mathrm{SNR}}---\left(31\right)$>

其中，A，B是拟合直线的系数，SNR为信噪比。待定系数A和B是获取海浪有效波高的关键经验参数，需要通过长期大量的试验数据进行标定。

信噪比的定义和计算方式有很多种，在本文中，定义如下：

$>\mathrm{SNR}=\frac{\mathrm{SIG}}{\mathrm{BGN}}---\left(32\right)$>

其中，SIG为海浪谱能量，BGN为背景噪声能量。在深水情况SIG和BGN可以由下面公式定义：

$>\mathrm{SIG}=\underset{1}{\overset{N_x}{\Sigma }}\underset{1}{\overset{N_y}{\Sigma }}E(k_x,k_y){\mathrm{dk}}_x{\mathrm{dk}}_y---\left(33\right)$>

$>\mathrm{BGN}=\underset{1}{\overset{N_x}{\Sigma }}\underset{1}{\overset{N_y}{\Sigma }}\underset{i=1}{\overset{N}{\Sigma }}I(k_x,k_y,{\omega }_i)·\mathrm{d\omega }{\mathrm{dk}}_y{\mathrm{dk}}_x-\underset{1}{\overset{N_x}{\Sigma }}\underset{1}{\overset{N_y}{\Sigma }}I(k_x,k_y){\mathrm{dk}}_y{\mathrm{dk}}_x---\left(34\right)$>

(2)波峰峰向和波周期的反演方法

直角坐标下的波数能谱与极坐标下的波数能谱的转换关系如下所示

E(k,θ)＝E(k_x,k_y)k>

其中

$>k=\sqrt{k_x^2+k_y^2},$>θ＝arctan2(k_y,k_x)+2π。

得到的角度范围在[-π,π]之间。

将海浪谱E(k,θ)转换为方向波谱E(f,θ)：

$>E(f,\theta )=E(k,\theta )\frac{\mathrm{dk}}{\mathrm{df}}=\frac{k}{2\pi }\frac{\partial k}{\partial \omega }E(k,\theta )---\left(36\right)$>

进一步针对方向频谱E(f,θ)中的方向θ积分，就可以获得一维频率谱，即

$>S\left(f\right)={\int }_0^{2\pi }E(f,\theta )\mathrm{d\theta }---\left(37\right)$>

进一步针对方向频谱E(f,θ)中的频率f积分，就可以获得一维方向谱，即

$>S\left(\theta \right)={\int }_0^{f_{\max }}E(f,\theta )\mathrm{df}---\left(38\right)$>

在得到S(θ)和S(f)的以后，可计算出海浪的周期，波向等参数。对于一维方向谱S(θ)，假设波谱峰值所在位置的坐标是θ，则有波向为

Pdir＝θ

从一维谱S(f)中得到可以得到波峰周期T_p和主波波长λ_P。

T_p＝1/f_p

其中，f_p为S(f)的最大值对应的频率，可以由下面公式算出

$>f_p=\frac{{\int }_{f_1}^{f_2}S\left(f\right)\mathrm{fdf}}{{\int }_{f_1}^{f_2}S\left(f\right)\mathrm{df}}---\left(39\right)$>

这里f₁和f₂由能量谱中的最大能量的80％来确定。再由

ω_p＝2π/T_P，k_p＝2π/λ_P

代入色散关系，可求得

$>{\lambda }_P=\frac{g}{2\pi ·{f_p}^2}---\left(40\right)$>

应用本发明提出的适用于动态环境下的基于色散关系的带通滤波器对2011年11月22-29日在渤海海域航行中的动态X波段雷达回波图像进行滤波处理。X波段雷达的架设高度为40米，平均旋转周期为2.39妙，工作在短脉冲模式，探测半径约为2千米海域。每3分钟采集一组数据，每组包含32幅图像。使用现场WAMOS系统提供的有效波高作为波高的参考值，WAMOS每2分钟输出以此结果；以WAVEX雷达提供的波周期和波峰峰向作为周期和波峰峰向的参考值，WAVEX每4分钟输出以此结果。实验图像区域中心位于船艏75度方向，面积960*960平方米，此海域平均水深20米。为尽量减小误差，该实验选取20分钟内得到的反演数据对波高，波周期与波峰峰向进行算术平均。并与WAMOS，WAVEX20分钟平均后的数值作对比。图7为应用本发明在动态环境下反演得到的有效波高与WAMOS测得有效波高的对比图，可见其走势基本相同，表二所示标准均方根误差为0.18低于设定的标准0.5m；图8为应用本发明反演得到的波峰风向与WAVEX所测对比，表二所示偏差为7.02°低于设定的标准20°；图9为应用本发明反演得到的海浪周期与WAVEX所测对比，表二所示相对标准误差为2％低于设定的标准10％。因此该试验有效，并说明达到了理想的滤波效果。

选取32个频率片中的某一片观察滤波效果，该片频率为0.74rad/s。图10为该频率下的海杂波能量谱，经过本发明设计的带通滤波器后得到的海浪信号图像谱如图11。比对图12可见目前的主流滤波器B_n、B_P由于运动状态下设置了较大的U_max已经无法起到滤波效果，而本发明所设计的滤波器滤波效果良好，保留了大部分真实的数与海浪的能量。

通过上述分析与实例计算，本发明设计的带通滤波器是简单有效的。

表一：实验参数设置

参数数值频率分辨率Δω0.082s^-1波数分辨率Δk0.0065m^-1空间x向分辨率Δx7.5m空间y向分辨率Δy7.5m非线性能量校正指数β1.2

表二：航行实验误差统计

误差统计波高波峰峰向波周期hrb均值1.3466.596.18WAMOS/WAVEX均值1.3469.26.23平均差02.050.05均方差0.0349.280.03标准均方根误差-偏差0.187.020.18相对标准误差2.00％