CN102353946B - 一种基于x波段雷达图像的海表面流反演方法 - Google Patents

Abstract

本发明公开的一种基于X波段雷达图像的海表面流反演方法，包括(1)采集时空域海杂波图像作为一个序列，得到子图像序列；(2)对子图像序列进行三维傅里叶变换；(3)根据色散关系构造带通滤波器；(4)对图像谱非线性影响进行校正；(5)计算隶属度；(6)加权计算；(7)初始估流；(8)迭代估流。本发明与现有流反演算法相比，提高了流反演精度，尤其是提高了低流速时流反演精度，且提高了反演结果稳定性。本发明使用依赖最大流速的色散关系带通滤波器对图像谱进行噪声滤除，去噪能力强。本发明中对图像谱进行校正得到海浪谱，能够反映真实的海态信息，然后将海浪谱作为最小二乘法的权值之一计算流，流反演结果接近真实海态。

Description

一种基于X波段雷达图像的海表面流反演方法

技术领域

本发明属于海浪参数反演技术领域，具体涉及一种基于X波段雷达图像的海表面流反演方法。

背景技术

海流是由于相邻海区间，海水长期存在温度、密度、气压的不同或长期受定向风的作用使海水产生水平方向的流动。它与海洋渔业、海上工程建设、海上航行安全和军事作战等有着十分密切关系。

X波段雷达工作在短脉冲模式时，发射的电磁波遇到海面上微尺度波发生Bragg散射，产生后向回波被雷达接收机收到，即形成海杂波图像。海杂波图像中含有丰富的海浪、海流信息，国内已有多家单位开展了相关研究，现有的传统的流反演算法包括以下几个步骤，如图1所示：

(1)雷达工作在短脉冲模式，采集N幅空间域海杂波连续图像，然后选取分析区域构成子图像序列I⁽³⁾(x，，t)。

(2)对子图像进行3维傅里叶变换得到三维图像谱F⁽³⁾(k，ω)。

(3)初始估流：

①海浪本身波长和频率符合色散关系，如式(1)所示。

${\mathrm{\ω}}_p=(p+1)\sqrt{\frac{g\left|k\right|}{p+1}\tanh \left(\frac{\left|k\right|d}{p+1}\right)}+k\·u---\left(1\right)$

式中，p为阶次，取0、1时分别代表0阶、1阶色散关系；ω_p为p阶次波频率；k为波数；u为雷达天线与海浪场间相对速度；d为海水深度。

②假定图像谱全部在0阶次色散关系上，取阈值C_fg，该阈值应高于非线性能量及背景噪声能量，能量高于C_fg的谱分量数目为n₀，一般选择C_fg的量值为0.2，n₀值在数十个。

③取极小值函数如下式所示：

$Q^2=\underset{i=1}{\overset{n_0}{\mathrm{\Σ}}}{({\mathrm{\ω}}_i-{\mathrm{\ω}}_0\left(k_2\right))}^2{F}^{\left(3\right)}(k_i,\mathrm{\ω})---\left(2\right)$

式中，ω_i为雷达观测的第i点的海浪频率其对应的能量为F⁽³⁾(k_i，ω)；LSM中第i个谱分量权值为它本身能量F⁽³⁾(k_i，ω)；n₀为LSM中谱分量点数目；ω₀(k_i)为理论上0阶色散关系给出的第i个谱分量频率。

④应用能量加权LSM，可得到较低精度的流。

(4)迭代估流：

①取一个比C_fg小的多的阈值C_it，要求非线性海浪能量大于C_it，能量高于C_it的谱分量数目为n₁。利用初始估值得到的较低精度流，由式(1)计算0阶次和1阶次波频率。

②判断实测的n₁个谱分量属于0阶次还是1阶次波，若|ω_i-ω₀(k_i)|＜|ω_i-ω₁(k_i)|，则属于0阶次波，若|ω_i-ω₀(k_i)|＞|ω_i-ω₁(k_i)|，则属于1阶次波。将判断好的数据根据不同极小值函数，应用能量加权LSM得到新表面流。

③将新表面流代入式(1)，得到新的0阶次和1阶次波色散关系，重复上述步骤(一般重复步数为10次可满足工程上精度需要)，将得到不断准确的表面流。

目前测表面流算法中存在着如下不足：①流速较低时，测量结果误差偏大；②测量结果不稳定，经常出现奇异值，在流速低时表现的尤为突出。

发明内容

针对现有技术中存在的缺点，本发明公开了一种基于X波段雷达图像的海表面流反演方法。本发明提出的方法与现有技术的显著区别是：①根据海浪特性对图像谱进行了噪声滤除；②对图像谱进行非线性校正；③根据图像谱结构风险最小化原则，引入隶属度函数加权最小二乘法解算流。

本发明提出一种基于X波段雷达图像的海表面流反演方法，包括以下步骤：

(1)雷达工作在短脉冲模式，连续采集时间长度为t的N幅空间域海杂波连续图像作为一个序列，然后选取分析区域构成子图像序列I⁽³⁾(x，y，t)；

(2)对子图像序列进行三维傅里叶变换得到三维波数频率图像谱F⁽³⁾(k，ω)；

(3)根据色散关系构造带通滤波器，滤除三维波数频率图像谱中的非海浪信号，所述的滤波器满足

$B_n=\frac{(\mathrm{\ω}-\mathrm{\Δ\ω}/2-U_{\max }{\mathrm{\ω}}^2/g)}{g}-\frac{\sqrt{2}\mathrm{\Δk}}{2}$

$B_n=\frac{(\mathrm{\ω}+\mathrm{\Δ\ω}/2+U_{\max }{\mathrm{\ω}}^2/g)}{g}+\frac{\sqrt{2}\mathrm{\Δk}}{2}$

式中，B_n和B_p均为滤波器带宽；ω为海浪频率，Δω为频率分辨率；k为波数，Δk为波数分辨率；U_max为雷达天线与海浪场最大相对流速；E⁽³⁾(k，ω)为滤波后海浪图像谱；

(4)引入调制传递函数MTF对图像谱E⁽³⁾(k，ω)的非线性影响进行校正，得到海浪谱，调制传递函数MTF满足：

MTF＝|k|^β

$E_M^{\left(3\right)}(k,\mathrm{\ω})={\left|k\right|}^{\mathrm{\β}}\·E^{\left(3\right)}(k,\mathrm{\ω})$

其中，β为指数，

为海浪谱；

(5)对于固定频率ω_i上有n_i个海浪谱分量，海浪谱分量波数的模组成集合为{k_j|k₁ k₂...k_ni}，中心为k_i0，集合最大半径ri定义为：

$k_{i0}=\frac{1}{n_i}\underset{j=1}{\overset{n_i}{\mathrm{\Σ}}}{k}_j$

$r_i=\underset{k_j}{\max }\left|\right|k_j-k_{i0}\left|\right|$

式中，k_j为海浪谱分量，海浪谱分量波数的模组成的集合中谱分量隶属度μ(k_ij)为：

$\mathrm{\μ}\left(k_{\mathrm{ij}}\right)=1-\frac{\left|\right|k_j-k_{i0}\left|\right|}{r_i+\mathrm{\δ}}$

其中，δ为表征k_j距离中心k_i0距离的参数；

(6)采用步骤(4)校正后海浪谱和步骤(5)得到隶属度函数的的乘积加权LSM，得到极小值函数Q²为：

$Q^2=\underset{i=1}{\overset{N_{\mathrm{\ω}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_i}{\mathrm{\Σ}}}{({\mathrm{\ω}}_i-{\mathrm{\ω}}_p\left(k_{\mathrm{ij}}\right))}^2{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)$

其中，ω_i为实测海浪谱数据

中第i个频率分量，ω_p为根据色散关系给出的频率分量；N_ω为正频率数目，N_i为ω_i频率上谱分量数目；μ(k_ij)分别为第ω_i频率上第j个谱大小和隶属度大小，k_ij分别为第i个频率中第j个谱分量的波数分量；

(7)初始估流：

当阶次p满足p＝0的0阶色散关系时，设定阈值C_fg满，将极小值函数Q²分别对表面流在X方向的分量u_x、表面流在Y方向的分量u_y求偏导，并使其偏导数均为零，根据谱数量n₂的个数计算表面流公式为：

u＝A·B

$A={\left[\begin{array}{cc}\mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x^2& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y\\ \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y^2\end{array}\right]}^{-1}$

$B=\left[\begin{array}{c}\mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x({\mathrm{\ω}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k|d})\\ \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y({\mathrm{\ω}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k|d})\end{array}\right]$

其中，k_x、k_y为波数k_ij在X方向和Y方向的分量，ω₀为0阶波频率，A和B表示矩阵；

(8)迭代估流：

①先选取小于阈值C_fg的阈值C_it，使得阈值C_it高于背景噪声，然后根据步骤(7)中的初始估流过程得到流计算为0和1阶次波频率；

②判断大于阈值C_it的n₃个海浪谱分量属于0阶还是1阶，若|ω_i-ω₀(k_i)|＜|ω_i-ω₁(k_i)|，则属于0阶次波，若|ω_i-ω₀(k_i)|＞|ω_i-ω₁(k_i)|，则属于1阶次波，其中ω_i为实测海浪谱数据中第i个频率分量，ω₀(k_i)为0阶波频率，ω₁(k_i)为1阶波频率；判断完成后，根据步骤(6)得到两个极小值函数，应用加权LSM得到新的表面流，计算公式为：

u＝A′·B′              (15)

$A^{\′}={\left[\begin{array}{cc}\mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x^2& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y\\ \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y^2\end{array}\right]}^{-1}---\left(16\right)$

$B^{\′}=\left[\begin{array}{c}\underset{n_{30}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x(\mathrm{\ω}-\sqrt{g\left|k\right|\tanh \left(\right|k|d})+\underset{n_{31}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x(\mathrm{\ω}-2\sqrt{\frac{g\left|k\right|}{2}\tanh (\frac{\left|k\right|d}{2}})\\ \underset{n_{30}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)k_y(\mathrm{\ω}-\sqrt{g\left|k\right|\tanh \left(\right|k|d})+\underset{n_{31}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y(\mathrm{\ω}-2\sqrt{\frac{g\left|k\right|}{2}\tanh (\frac{\left|k\right|d}{2}})\end{array}\right]---\left(17\right)$

其中，n₃₀、n₃₁分别为n₃个海浪谱分量中0、1阶次波数量，u为表层流，A′和B′表示矩阵；

③将新的表面流代入色散关系式

中，得到新的0阶次和1阶次波频率，返回步骤②，进行迭代计算，得到不断精确的表层流，得到流速与流向。

本发明的优点在于：

(1)；本发明提出的一种基于X波段雷达图像的海表面流反演方法，与现有流反演算法相比，提高了流反演精度，尤其是提高了低流速时流反演精度；

(2)；本发明提出的一种基于X波段雷达图像的海表面流反演方法，与现有流反演算法相比，提高了反演结果稳定性。

(3)；本发明提出的一种基于X波段雷达图像的海表面流反演方法，使用依赖最大流速的色散关系带通滤波器对图像谱进行噪声滤除，与现有流反演算法使用的滤波方法相比，去噪能力更强。

(4)本发明提出的一种基于X波段雷达图像的海表面流反演方法，对图像谱进行校正得到海浪谱，海浪谱能够反映真实的海态信息，然后把海浪谱作为最小二乘法的权值之一计算流，与现有流反演算法相比，流反演结果更接近真实海态。

(5)本发明提出的一种基于X波段雷达图像的海表面流反演方法，根据结构风险最小化原则，把谱分量隶属度作为最小二乘的权值之一计算流，与现有流反演算法相比，有效的控制了强噪声、强边缘海浪信号的影响。

附图说明

图1：现有技术中传统的测流方法的流程图；

图2：本发明提出的一种基于X波段雷达图像的海表面流反演方法的流程图；

图3：雷达实测区域示意图；

图4：流速结果比较示意图；

图5：流向结果比较示意图。

具体实施方式

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

本发明公开的一种基于X波段雷达图像的海表面流反演方法，如图2所示，具体包括以下步骤：

(1)雷达工作在短脉冲模式，连续采集时间长度为t的N幅空间域海杂波连续图像作为一个序列(N根据连续采样时间长度t和雷达转速计算得到，其数值为连续采样时间长度t除以雷达转速后取整数)，然后选取分析区域构成子图像序列I⁽³⁾(x，y，t)。在文献，周蓓硕士论文，2008年一文X波段雷达海表面流信息提取技术研究中给出连续采集时间长度t应取40～100s。在本发明中依据此文献取连续采集时间长度t为40～100s。

(2)对子图像序列进行三维傅里叶变换得到三维波数频率图像谱F⁽³⁾(k，ω)。

(3)根据色散关系构造带通滤波器，滤除三维波数频率图像谱中的非海浪信号，所述的滤波器为不依赖实测流速的带通滤波器，满足

$B_n=\frac{(\mathrm{\ω}-\mathrm{\Δ\ω}/2-U_{\max }{\mathrm{\ω}}^2/g)}{g}-\frac{\sqrt{2}\mathrm{\Δk}}{2}---\left(4\right)$

$B_n=\frac{(\mathrm{\ω}+\mathrm{\Δ\ω}/2+U_{\max }{\mathrm{\ω}}^2/g)}{g}+\frac{\sqrt{2}\mathrm{\Δk}}{2}---\left(5\right)$

式中，B_n和B_p均为滤波器带宽；ω为海浪频率，Δω为频率分辨率；k为波数，Δk为波数分辨率；U_max为粗略估计的雷达天线与海浪场最大相对流速；E⁽³⁾(k，ω)为滤波后海浪图像谱。

经滤波后非海浪信号被滤除，得到滤波后海浪图像谱E⁽³⁾(k，ω)，仅含有海浪信号，所述的滤波器适用于流速较小(＜1.5m/s)的背景下，如岸基或海上石油平台背景。

(4)由于图像谱和海浪谱，仍有区别，引入调制传递函数(modulation transfer function，MTF)对图像谱E⁽³⁾(k，ω)的非线性影响进行校正得到海浪谱，大量实验事实表明，MTF按指数衰减，发明中只关心谱的相对值，因此MTF由下式给出。

MTF＝|k|^β                     (6)

$E_M^{\left(3\right)}(k,\mathrm{\ω})={\left|k\right|}^{\mathrm{\β}}\·E^{\left(3\right)}(k,\mathrm{\ω})---\left(7\right)$

式中，β为指数，β值根据长时间大量的实验标定，它取决于雷达成像的非线性，对于同一部雷达只需要标定一次即可，在文献：崔利民、何宜军，IEEE International Geoscience andRemote Sensing Symposium，2009年一文MEASURMENTS OF OCEAN WAVE SPECTRAWITH VERTICAL POLARIZATION X-BAND RADAR IMAGE SEQUENCES中，作者标定出X波段雷达水平极化下β取值1.17；在文献：周蓓硕士论文，2008年一文X波段雷达海表面流信息提取技术研究中，标定出X波段雷达水平极化下β取值1.21；在文献：P.Izquierdo、J.C.Nieto Borge、C.Guedes Soares等，JOURNAL OF WATERWAY，PORT，COASTAL，ANDOCEAN ENGINEERING期刊，2005年一文Comparison of Wave Spectra from NauticalRadar Images and Scalar Buoy Data中使用β值1.2。本专利中，考虑到β值取决于雷达成像机制的非线性，要求β值范围为1.17≤β≤1.21，雷达E⁽³⁾(k，ω)为图像谱，

为海浪谱。

(5)基于海浪谱分量到谱中心的距离来度量隶属度函数大小，海浪谱分量距离谱中心越近，隶属度越大，反之，则越小。对于固定频率ω_i上有n_i个海浪谱分量，它们波数的模组成集合为{k_j|k₁ k₂...k_ni}，中心为k_i0，集合最大半径ri定义如下：

$k_{i0}=\frac{1}{n_i}\underset{j=1}{\overset{n_i}{\mathrm{\Σ}}}{k}_j---\left(8\right)$

$r_i=\underset{k_j}{\max }\left|\right|k_j-k_{i0}\left|\right|---\left(9\right)$

式中，k_j为海浪谱分量。根据距离(最大半径r_i)确定隶属度时，海浪谱分量波数模集合中谱分量隶属度μ(k_ij)为：

$\mathrm{\μ}\left(k_{\mathrm{ij}}\right)=1-\frac{\left|\right|k_j-k_{i0}\left|\right|}{r_i+\mathrm{\δ}}---\left(10\right)$

其中δ为很小的正数，表征k_j距离中心k_i0距离的参数，本发明中δ值数量级范围为10^-4～10^-6，避免出现μ(k_ij)＝0的情况。

(6)为根据色散关系对实测海浪谱拟合得到表面流，采用步骤(4)校正后海浪谱和步骤(5)得到隶属度函数的的乘积加权LSM(最小二乘法加权)，得到极小值函数Q²为：

$Q^2=\underset{i=1}{\overset{N_{\mathrm{\ω}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_i}{\mathrm{\Σ}}}{({\mathrm{\ω}}_i-{\mathrm{\ω}}_p\left(k_{\mathrm{ij}}\right))}^2{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)---\left(11\right)$

式中，ω_i为实测海浪谱数据

中第i个频率分量，ω_p为根据理论上色散关系(背景技术中的公式(1))给出的频率分量；N_ω为正频率数目，N_i为ω_i频率上谱分量数目；

μ(k_ij)分别为第ω_i频率上第j个谱大小和隶属度大小，k_ij分别为第i个频率中第j个谱分量的波数分量。

(7)初始估流：仅考虑背景技术公式(1)中阶次p＝0时色散关系，取阈值C_fg，使得谱能量大于阈值C_fg的谱数量n₂。在文献：周蓓硕士论文，2008年一文X波段雷达海表面流信息提取技术研究，给出C_fg量值为0.7；在文献：Senet C M，Seemann J，Ziemer F，IEEE Transactionson Geoscience and Remote Sensing期刊，2001年一文The Near-Surface Current VelocityDetermined from Image Sequences of the Sea Surface中给出C_fg量值为0.2。本专利中考虑阈值C_fg主要靠经验给出，选择有很大的余地，但要求0.2≤C_fg≤0.7。将极小值函数Q²分别对表面流在X方向的分量u_x、表面流在Y方向的分量u_y求偏导(X向是指水平面坐标系中东向，Y向是指水平面坐标系中北向)，并使其偏导数均为零，根据谱数量n₂的个数计算表面流公式如下：

u＝A·B                     (12)

$A={\left[\begin{array}{cc}\mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x^2& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y\\ \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y^2\end{array}\right]}^{-1}---\left(13\right)$

$B=\left[\begin{array}{c}\mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x({\mathrm{\ω}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k|d})\\ \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y({\mathrm{\ω}}_0-\sqrt{g\left|k\right|\tanh \left(\right|k|d})\end{array}\right]---\left(14\right)$

式中，k_x、k_y为波数k_ij在X方向和Y方向的分量，ω₀为理论上0阶波频率，可由色散关系公式(1)p＝0时求出。A和B表示矩阵。

(8)迭代估流：由于大于1阶的海浪谱分量较少，因此仅考虑0、1阶次波影响(大于1阶谱分量较少)。

①先选取比阈值C_fg小的多的阈值C_it，使得C_it高于背景噪声，这样图像谱中大于阈值C_it部分包含了非线性能量。在文献：Senet C M，Seemann J，Ziemer F，IEEE Transactions onGeoscience and Remote Sensing期刊，2001年一文The Near-Surface Current Velocity Determinedfrom Image Sequences of the Sea Surface中给出C_it量值为0.02；大量的实验发现阈值C_it在0.015～0.03内取值时具有较好效果，在本专利中要求阈值C_it满足0.015≤C_it≤0.03，然后根据步骤(7)所述的初始估流得到流计算为0和1阶次波频率。

②判断大于阈值C_it的n₃个海浪谱分量属于0阶还是1阶，若|ω_i-ω₀(k_i)|＜|ω_i-ω₁(k_i)|，则属于0阶次波，若|ω_i-ω₀(k_i)|＞|ω_i-ω₁(k_i)|，则属于1阶次波。其中ω_i为实测海浪谱数据

中第i个频率分量，ω₀(k_i)为理论上0阶波频率，可由公式(1)p＝0时求出；ω₁(k_i)理论上1阶波频率，可由式(1)p＝1时求出。将判断好的数据根据步骤(6)得到不同的极小值函数，这样就有了两个极小值函数(当为0阶次波时，步骤(6)中ω_p＝ω₀；当为1阶次波时ω_p＝ω₁)，应用加权LSM得到新的表面流，计算公式如下：

u＝A′·B′                       (15)

$A^{\′}={\left[\begin{array}{cc}\mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x^2& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y\\ \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y^2\end{array}\right]}^{-1}---\left(16\right)$

$B^{\′}=\left[\begin{array}{c}\underset{n_{30}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x(\mathrm{\ω}-\sqrt{g\left|k\right|\tanh \left(\right|k|d})+\underset{n_{31}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x(\mathrm{\ω}-2\sqrt{\frac{g\left|k\right|}{2}\tanh (\frac{\left|k\right|d}{2}})\\ \underset{n_{30}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)k_y(\mathrm{\ω}-\sqrt{g\left|k\right|\tanh \left(\right|k|d})+\underset{n_{31}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y(\mathrm{\ω}-2\sqrt{\frac{g\left|k\right|}{2}\tanh (\frac{\left|k\right|d}{2}})\end{array}\right]---\left(17\right)$

式中，n₃₀、n₃₁分别为n₃个海浪谱分量中0、1阶次波数量，u为表层流，A′和B′表示矩阵。

③将新表面流代入色散关系式(1)

中，得到新的0阶次频率ω₀和1阶次波频率ω₁，返回步骤②，进行迭代计算，每经过一次迭代，表面流流速和流向是变化的，即公式(1)中u是变化的，每次迭代过程都会重新确定色散关系方程，返回步骤②后ω₀和ω₁是更新的。参考文献：崔利民博士论文“X波段雷达海浪和海流遥感机理及信息提取算法”给出迭代次数范围为10～15次，经大量实验发现迭代次数大于10次时流反演结果精度提高很小，在本发明中考虑到节省算法时间，取迭代次数选取为10次，将得到不断精确流速与流向。

应用本发明提出的一种基于X波段雷达图像的海表面流反演方法进行实地实测，实验地点在国家海洋局平潭水文观测站附近海域，受地形影响此海域有较强的潮流，遇大潮时潮差超过6米。X波段雷达的架设高度为40米，平均旋转周期为2.5秒，工作在短脉冲模式，雷达探测半径约为2千米海域。每4分钟采集一组数据，每组包含32幅图像，2009年10月6日至8日，共采集不连续数据780组。使用已标定现场传感器测量的表面流作为真值，现场传感器每4分钟输出一次结果。现场传感器与本次实验测量区域相同，测量区域中心位于船艏75度方向，近点离岸600米，面积960*960m²，此海域平均水深25米，如图3所示。为验证本发明提出的基于X波段雷达图像的海表面流反演方法的性能，将应用传统算法(流反演算法)和本发明得到的结果分别与现场传感器的真实值进行比对，然后给出误差统计，本发明实验中参数设置如表1所示。

表1：实验参数设置

参数	数值
		频率分辨率Δω	0.079s-1
波数分辨率Δk	0.0093m-1
		空间x向分辨率Δx	7.5m
空间y向分辨率Δy	7.5m
		实验海区平均水深d	25m
实验海区最大流速Umax	1.5m/s
		非线性能量校正指数β	1.2
图像谱中正频率数目Nω	16
		初始估流阈值Cit	0.2
迭代估流阈值Cfg	0.02

为清晰的看出改进算法本发明的性能，图4中分别给出了现场传感器流速与本发明得到的改进流速的对比、现场传感器流速与原流速对比、改进流速与原流速对比；图5中分别给出了现场传感器流向与本发明得到的改进流向对比、现场传感器流向与原流向对比、改进流向与原流向对比；本发明结果和现有原算法(流反演算法)结果的误差统计特性如表2所示，本发明流速较原流速精度提高66.7％，流向精度提高61.3％。

表2本发明与现有算法的误差统计特征表

Claims (6)

1.一种基于X波段雷达图像的海表面流反演方法，其特征在于：包括以下步骤：

（1）雷达工作在短脉冲模式，连续采集时间长度为t的N幅空间域海杂波连续图像作为一个序列，然后选取分析区域构成子图像序列I⁽³⁾(x,y,t)；

（2）对子图像序列进行三维傅里叶变换得到三维波数频率图像谱F⁽³⁾(k,ω)；

（3）根据色散关系构造带通滤波器，滤除三维波数频率图像谱中的非海浪信号，所述的滤波器满足

$B_n=\frac{(\mathrm{\ω}-\mathrm{\Δ\ω}/2-U_{\max }{\mathrm{\ω}}^2/g)}{g}-\frac{\sqrt{2}\mathrm{\Δk}}{2}$

$B_n=\frac{(\mathrm{\ω}+\mathrm{\Δ\ω}/2+U_{\max }{\mathrm{\ω}}^2/g)}{g}+\frac{\sqrt{2}\mathrm{\Δk}}{2}$

式中，B_n和B_p均为滤波器带宽；ω为海浪频率，Δω为频率分辨率；k为波数，Δk为波数分辨率；U_max为雷达天线与海浪场最大相对流速；E⁽³⁾(k,ω)为滤波后海浪图像谱；

（4）引入调制传递函数MTF对图像谱E⁽³⁾(k,ω)的非线性影响进行校正，得到海浪谱，调制传递函数MTF满足：

MTF＝|k|^β

$E_M^{\left(3\right)}(k,\mathrm{\ω})={\left|k\right|}^{\mathrm{\β}}\·E^{\left(3\right)}(k,\mathrm{\ω})$

其中，β为指数，

为海浪谱；

（5）对于固定频率ω_i上有n_i个海浪谱分量，海浪谱分量波数的模组成集合为{k_j|k₁ k₂...k_ni}，中心为k_i0，集合最大半径r_i定义为：

$k_{i0}=\frac{1}{n_i}\underset{j=1}{\overset{n_i}{\mathrm{\Σ}}}{k}_j$

$r_i=\underset{k_j}{\max }\left|\right|k_j-k_{\mathrm{io}}\left|\right|$

式中，k_j为海浪谱分量，海浪谱分量波数的模组成的集合中谱分量隶属度μ(k_ij)为：

$\mathrm{\μ}\left(k_{\mathrm{ij}}\right)=1-\frac{\left|\right|k_j-k_{\mathrm{io}}\left|\right|}{r_i+\mathrm{\δ}}$

其中，δ为表征k_j距离中心k_i0距离的参数；

（6）采用步骤（4）校正后海浪谱和步骤（5）得到隶属度函数的的乘积加权LSM，得到极小值函数Q²为：

$Q^2=\underset{i=1}{\overset{N_{\mathrm{\ω}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_i}{\mathrm{\Σ}}}{({\mathrm{\ω}}_i-{\mathrm{\ω}}_p\left(k_{\mathrm{ij}}\right))}^2{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)$

其中，ω_i为实测海浪谱数据

中第i个频率分量，ω_p为根据色散关系给出的频率分量；N_ω为正频率数目，N_i为ω_i频率上谱分量数目；

μ(k_ij)分别为第ω_i频率上第j个谱大小和隶属度大小,k_ij分别为第i个频率中第j个谱分量的波数分量；

（7）初始估流：

当阶次p满足p=0的0阶色散关系时，设定阈值C_fg，将极小值函数Q²分别对表面流在X方向的分量u_x、表面流在Y方向的分量u_y求偏导，并使其偏导数均为零，根据谱数量n₂的个数计算表面流公式为：

u=A·B

$A={\left[\begin{array}{cc}\mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x^2& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y\\ \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y^2\end{array}\right]}^{-1}$

$B=\left[\begin{array}{c}\mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x({\mathrm{\ω}}_o-\sqrt{g\left|k\right|\tanh \left(\right|k\left|d\right)}\\ \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y({\mathrm{\ω}}_o-\sqrt{g\left|k\right|\tanh \left(\right|k\left|d\right)}\end{array}\right]$

其中，k_x、k_y为波数k_ij在X方向和Y方向的分量，ω₀为0阶波频率，A和B表示矩阵；

（8）迭代估流：

①先选取小于阈值C_fg的阈值C_it，使得阈值C_it高于背景噪声，然后根据步骤（7）中的初始估流过程得到流计算为0和1阶次波频率；

②判断大于阈值C_it的n₃个海浪谱分量属于0阶还是1阶，若|ω_i-ω₀(k_i)|＜|ω_i-ω₁(k_i)|，则属于0阶次波，若|ω_i-ω₀(k_i)|＞|ω_i-ω₁(k_i)|，则属于1阶次波，其中ω_i为实测海浪谱数据

中第i个频率分量，ω₀(k_i)为0阶波频率，ω₁(k_i)为1阶波频率；判断完成后，根据步骤(6)得到两个极小值函数，应用加权LSM得到新的表面流，计算公式为：

u=A′·B′   (15)

$A\′={\left[\begin{array}{cc}\mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x^2& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y\\ \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x{k}_y& \mathrm{\Σ}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y^2\end{array}\right]}^{-1}---\left(16\right)$

$B^{\′}=\left[\begin{array}{c}\underset{n_{30}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x(\mathrm{\ω}-\sqrt{g\left|k\right|\tanh \left(\right|k\left|d\right)}+\underset{n_{31}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_x(\mathrm{\ω}-2\sqrt{\frac{g\left|k\right|}{2}\tanh (\frac{\left|k\right|d}{2}})\\ \underset{n_{30}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)k_y(\mathrm{\ω}-\sqrt{g\left|k\right|\tanh \left(\right|k\left|d\right)}+\underset{n_{31}}{\mathrm{\Σ}}{E}_M^{\left(3\right)}(k_{\mathrm{ij}},{\mathrm{\ω}}_i)\mathrm{\μ}\left(k_{\mathrm{ij}}\right)k_y(\mathrm{\ω}-2\sqrt{\frac{g\left|k\right|}{2}\tanh (\frac{\left|k\right|d}{2}})\end{array}\right]---\left(17\right)$

其中，n₃₀、n₃₁分别为n₃个海浪谱分量中0、1阶次波数量，u为表层流，A′和B′表示矩阵，d为海水深度；

③将新的表面流代入色散关系式

中，得到新的0阶次和1阶次波频率，返回步骤②，进行迭代计算，得到不断精确的表层流，得到流速与流向。

2.根据权利要求1所述的一种基于X波段雷达图像的海表面流反演方法，其特征在于：步骤（1）中时间长度t的取值满足40s≤t≤100s。

3.根据权利要求1所述的一种基于X波段雷达图像的海表面流反演方法，其特征在于：步骤（4）中指数β值的取值范围满足1.17≤β≤1.21。

4.根据权利要求1所述的一种基于X波段雷达图像的海表面流反演方法，其特征在于：步骤（5）中参数δ的数量级范围为10^-4～10^-6。

5.根据权利要求1所述的一种基于X波段雷达图像的海表面流反演方法，其特征在于：步骤（7）中取阈值C_fg满足0.2≤C_fg≤0.7。

6.根据权利要求1所述的一种基于X波段雷达图像的海表面流反演方法，其特征在于：步骤（8）中阈值C_it满足0.015≤C_it≤0.03。

Images