(4) 由单波空间场及其梯度图像得到波数场;

　　(5) 由单波空间场以及其对应的波数场得到5?D时空频率场[Ix，y|k，ω];

　　(6) 由局部的3?D图像谱反演局部的水深及流速。

　　该算法是针对由岸基X波段雷达获取的时间序列的雷达图像，最终得到水深场。

3 数值模拟及分析

　　3.1 模拟非均匀波场及雷达图像

　　基于线性波理论，海浪可看成是各种不同的余弦波的线性叠加，该过程可利用频谱来模拟，本文选用与波浪相近的P?M谱。只有频谱还不足以描述海浪的特性，需要加入方向分布函数组成方向谱，才能符合实际的海面波场状况，本文的方向分布函数采用改进的光易型方向分布函数。同时考虑到波场的非均匀性，加入非等水深值及表面流速值，利用色散关系式(1)，可确定不同区域的波数与频率的关系，利用不同频率和传播方向的余弦波的叠加，可模拟出浅水区的非均匀波场的时间序列。图2所示是模拟的64幅时间序列的非均匀波场的前两幅(图像中像素点的个数为128×128个，每个像素点的分辨率为7.5 m×7.5 m)。

　　图2 模拟的64幅时间序列的非均匀波场的前两幅

　　根据雷达成像机理，利用起主要作用的阴影调制及倾斜调制模拟出时间序列的雷达图像。图3所示是模拟的64幅时间序列的雷达图像的前两幅。

　　图3 模拟的64幅时间序列的雷达图像的前两幅

　　3.3 对模拟数据进行处理

　　(1) 对64幅时间序列的雷达图像[G(Θ)]进行三维傅里叶变化得到复数值的三维波数?频率谱：

　　(2)

　　其中三维谱的谱分辨率为：

　　[Δkx=2πX， Δky=2πY， Δω=2πT] (3)

　　(2) 对得到的三维谱进行阈值滤波，滤除信号中包含的噪声，然后利用色散关系进行带通滤波，得到海浪信号。接下来将对滤波后的三维谱进行分解，得到单波成分的波谱，既进行频率分解和方向分解。在时间轴上进行的傅里叶变化使得频率分解被执行，既一系列不同频率所对应的二维波数谱，接着进行方向分解。本文采用了一组楔形滤波器，首先产生一个原型楔形滤波器，然后再通过旋转，双线性插值，得到一组滤波器，原型滤波器如图4(a)所示，旋转得到的部分滤波器如图4(b)～(d)所示。运用这一组方向滤波器对二维谱进行分解，最终得到一系列不同频率和传播方向所对应的单一成分的波谱[Gk|ω，?]。

　　(3) 对单一成分的波谱[ Gkω，?]进行二维反傅里叶变化得到复数值的单波空间场[ gx，y|ω，?]：

　　[gx，y|ω，?=2D IFFT(Gk|ω，?)] (4)

　　图4 一组方向滤波器中的前四个

　　单波空间场包含了幅值及相位模式信息：

　　[gx，y|ω，?=Ax，y|ω，?expi?x，y|ω，?=Regx，y|ω，?+iImgx，yω，?] (5)

　　与单波空间场对应的梯度图像：

　　[??x，??ygx，yω，?=2D IFFTi?kx，ky?Gkω，?] (6)

　　其中[kx，ky]代表复数值的波数向量，其实部代表局部的波数值。局部区域的大小选为8×8个像素点，因此要得到局部区域的波数，需要分析局部点所包含的所有像素点。

　　位于色散关系滤波器带宽内的背景噪声重新分布在了单波波数场中，因此为了消除噪声的影响，运用方差最小拟合法得到复数值的波数向量。

　　[kx=-i?v+?vxv2ky=-i?v+?vyv2] (7)

　　其中向量[v，][vx，][vy]通过行扫描局部区域内的单波空间场及其梯度图像获得，向量[v+]是向量[v]的共轭向量。

　　(4) 由一系列的单波空间场以及单波波数场可得到五维的时空频率谱 [Ix，yk，ω。]表面波信号的能量谱应分布在色散关系曲面上，将由色散关系式(1)得到的谱分量[ω]与图像谱[Ix，y|k，ω]中的分量[ωi]取加权方差，得到一个函数。本文利用该函数寻找最小值的方法求得局部的流速[ux，uy]及水深[d。]该加权方差函数表示为：

　　[fux，uy，d=i=0Ngkitanhkid+kx，iux+ky，iuy-ωi2?Ix，y|ki，ωi] (8)

　　式中：[N]表示谱坐标集[{kx，i，ky，i，ωi}]中元素的个数，通过设置阈值从局部能量谱中选取出谱坐标集：

　　[M0=(kx，i，ky，i，ωi)Ix，y|ki，ωiMAXIx，y|ki，ωi>ε] (9)

　　式中[ε]表示能量阈值。

　　加权方差函数是一个非线性函数，含有三个变量，求该函数最小值属于优化问题，本文采用拟牛顿法搜索最小值，并得到局部的水深及流速。