基于X波段雷达图像反演近岸水深技术的创新

　0 引 言
　　对近岸浅水区域的海底地形的掌握对于近岸活动比如捕鱼，码头建设，铺设石油管道以及形态动力学的研究等都有着至关重要的作用。传统获取水深的方法包括船载声呐探测、机载激光测探、潜水器测量、超光谱图像测量等。但是它们普遍具有成本高、效率低的缺点，并且测量精度会受海水清澈度的影响。为了克服传统测量方法所存在的不足，基于X波段雷达图像序列反演水深的方法得到发展。
　　在浅水区域，局部的海底地形对表面波的传播有着重要的影响。当波移动至浅水区域，波的周期不发生变化，但是波的传播速度会发生变化，继而波长减小，波数增加。表面流的存在也影响表面波的传播，因此水深场以及表面流速场的反演方法都是基于这种传播变化——在物理学中被描述为表面波的色散关系。
　　1998年 Paul Bell运用连续的雷达图像序列之间的互相关性推导出了空间变化的表面波速[1]，同时利用从浮标中获取的频率信息，通过运用线性重力波的色散关系计算出了空间变化的水深，但是没有考虑海流的存在。Hessner等人运用一维FFT变换实现了对图像序列的频率分解[2]，某一固定频率的波所对应的波长通过确定局部空间的相位梯度计算得到。但是此方法的局限性在于它不能运用到包含同一频率但不同传播方向的波的波场，同时也没有考虑到海流的存在。
　　对时间序列的雷达图像进行3?D FFT变化，并取模的平方得到三维图像谱，由于波数和频率被色散关系联系在一起，因此线性表面波的信号应该很好地分布于其确定的三维形状上。色散关系的形状取决于水深和表面流速，因此通过拟合理论的色散关系和三维图像谱的坐标分布便可反演出大的空间范围内的平均水深以及流速[3?4]。但是此方法中的3?D FFT是针对全局范围的算子，因此假定了波场的均匀性以及稳定性。如果在深水区存在变化的流速或者浅水区存在变化的水深，波的折射将会产生，波场变成了非均匀场，以上方法不再适用，因此需要在局部空间范围内对波参数进行分析。
　　自1999年以来，Seemann等人针对非均匀波场做了一系列研究[5?10]，推导出了局部三维图像谱，同时反演出了局部范围内的水深以及流速。
　　本文将利用模拟的X波段雷达图像展开近岸浅水区域的水深的反演工作，该工作考虑到了表面波场的非均匀性，因此采用了局部反演算法，反演出了局部的水深值。
　　1 色散关系与水深以及流速的关系
　　色散关系描述了波数[k]和角频率[ω]之间的动力学关系，正常的色散关系适用于海表面重力波，线性色散关系可表示为：
　　[ωk，uc，d=±gktanhkd+] （1）
　式中：[g]表示为重力加速度；[d]为水深；[uc]为近表面流速。在式（1）中，第一部分称为固有频率[ζ=±gktanhkd，]第二部分称为多普勒频率[ωD=。]多普勒频率部分表明受表面流速的影响。在式（1）中，水深[d]和表面流速矢量[uc]在波数?频率域中影响色散关系的形状，因此色散关系的形状可以被用来反演这些参数值。图1显示了水深以及流速对色散关系的影响。
　　图1 三维波数?频率域中线性表面重力波的色散关系
　　2 水深及流速局部反演方法介绍
　　在浅水区域中，由于空间变化的水深，波的周期不变，既波场保持了稳定性，但是波长发生了变化，波场变成了非均匀场，因此需要在局部空间范围内对海态参数进行分析，得到空间分布的海态参数场。海洋表面波的特性由波长[λ、]波数[k、]角频率[ω、]振幅[ξ]和它们的传播方向[?]来描述。表面波场由一系列不同频率不同传播方向的单一成分的波（简称单波）叠加得到，因此其是多成分的，需要将其分解为单成分波。本文将按照以下步骤反演局部的水深及流速：
　　（1） 对时间序列的雷达图像进行3?D FFT变换，得到复数值的三维图像谱[G（k，ω）]；
　　（2） 对三维图像谱进行频率分解和方向分解得到单波成分的波谱 [Gk|ω，?]；
　　（3） 进行2?D 反FFT变化，到空间域，产生单波复数值的空间场[gx，y|ω，?]；
　　（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）
　　式中[ε]表示能量阈值。
　　加权方差函数是一个非线性函数，含有三个变量，求该函数最小值属于优化问题，本文采用拟牛顿法搜索最小值，并得到局部的水深及流速。
　　4 数据处理结果
　　模拟雷达图像时输入的非等值水深场如图5（a）所示，每8×8个像素点设置一个水深值，为减少模拟时的计算量，水深值只沿一维变化。反演得到的水深场如图5（b）所示，反演时选择的局部区域的大小为8×8个像素点。反演的水深值与输入的水深值吻合较好，平均误差约为2%，相比于过去的均匀场水深反演方法，该反演方法可将水深值的分辨率缩小到8×8个像素点。
　　5 结 语
　　利用X波段雷达图像可提取出重要的海态信息，比如水深、流速等等。均匀场的水深及流速的反演方法已相对成熟，本文的工作是针对非均匀场反演浅水水深值。由于实际的海况比较复杂，并且还没有得到可以用于比测的实际水深值，本文采用数值模拟的方法，通过输入非等值的水深仿真出非均匀波场及其雷达图像。利用仿真的雷达图像反演出局部水深值，并与输入的水深值进行对比，结果吻合较好，对利用实际的雷达图像反演非均匀场的水深具有重要的指导意义。本文的工作是基于岸基X波段雷达，对于船基X波段雷达来说，还要考虑运动补偿等因素，并且实际海况复杂多变，水深的反演过程有待进一步分析研究。
　　图5 输入的水深场与反演得到的水深场对比图
　　参考文献
　　. Coast. Eng.， 1999， 37（3）： 513?527.
　　， 1999： 16?20.
　　[3] SENET C， SEEMANN J， ZIEMER F. An iterative technique to determine the near surface current velocity from time series of sea surface images [C]// Proceedings of Oceans MTS/IEEE Conference?500 Years of Ocean Exploration. [S.l.]： IEEE，1997： 66?72.
　　 IEEE Transactions on Geoscience and Remote Sensing， 2002， 40（12）： 2606?2618.
　　. Procee?dings of SPIE 1999， 3808： 536?546.
　　： IEEE， 2000： 1898?1900.
　　[7] SEEMANN J， SENET C， ZIEMER F. Local analysis of inhomogeneous sea surfaces in coastal waters using nautical radar image sequences. Berlin， Germany： Springer?Verlag， 2000： 179?186.
　　. Proceedings of SPIE， 1999， 3808： 148?158.
　　： IEEE， 2000： 1769?1774.