CN102621531A - 一种基于x波段雷达图像的降雨干扰抑制方法 - Google Patents

Abstract

本发明公开了一种基于X波段雷达图像的降雨干扰抑制方法，属于遥感技术领域。所述的方法引入质量控制，确定存在问题的图像；接着引入三维表面参数来描述雷达回波图像的三维特性；然后结合信噪比及三维表面参数确定图像产生问题的原因，并确定降雨图像；最后利用频域滤波的方式处理被降雨影响较轻的图像。与现有降雨干扰抑制算法相比，本发明减少了需要进行识别判断的数据量，提高算法运行速度；比二维参数能提供更多的信息，描述更接近于真实表面，且能够给出所研究表面的直观图像以及与表面形态有关的足够信息；能更准确的确定雷达回波图像中受降雨影响的数据，以及影响的程度；能更准确地识别降雨资料；提高了海浪参数的反演精度。

Description

一种基于X波段雷达图像的降雨干扰抑制方法

技术领域

本发明属于遥感技术领域，具体涉及一种基于X波段雷达图像的降雨干扰抑制方法。

背景技术

在利用航海雷达进行海浪遥测时，经常存在一些干扰。这些干扰中最为常见的是同频干扰、固体目标物干扰和降雨干扰。在未降雨的资料中，图像信号在空间中的变化形式主要是近似于海面水面的波动，呈现周期性的变化。降雨的图像中，大片区域为降水影响，雨水在空间是随意分布的，所以导致图像中降雨信号的变动极为不规则，也因此造成了图像中各个像元的变化变大，使图像的变异数增大，降雨会影响雷达的测波能力，会改变海面的粗糙度，从而影响回波数据，最终影响海浪参数反演结果的准确度。

现有技术中在降雨数据的处理方面有一定的研究。利用海面回波强度的均值和差异系数作为判断标准，区分降雨和非降雨数据，并进行中值滤波处理；还有利用小波算法以及将中值滤波与小波算法相结合对降雨噪声进行抑制研究，都取得不错的效果。

已有的区分降雨与非降雨的方法中利用的是图像回波强度的均值及差异系数。

均值的定义如下所示：

$\mathrm{\μ}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i$

其中x_i表示雷达图像上第i个像素点的回波强度值，n表示像素点的个数。

差异系数CV的定义为：

$\mathrm{CV}=\frac{\mathrm{\σ}}{\mathrm{\μ}}\×100\%$

其中σ是图像回波强度的标准偏差，σ的定义为，

$\mathrm{\σ}=\sqrt{\frac{1}{n-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\mathrm{\μ})}^2}$

均值及差异系数等是评定表面粗糙度的二维参数，对于表面纹理的变化情况及排列规则，却不能反映出来；在降雨和非降雨资料之间均值及差异系数并没有非常明显的变化界限，在识别降雨资料时，准确性受限；中值滤波或者小波滤波的方法是对雷达图像各个像元之间非线性计算，不是具体针对降雨噪声的抑制方法。

发明内容

针对现有技术中的问题，本发明公开了一种基于X波段雷达图像的降雨干扰抑制方法。本发明提出的方法与现有技术的显著区别是：首先引入质量控制，确定存在问题的图像；接着引入三维表面粗糙度评定参数来描述雷达回波图像的三维特性；然后结合信噪比及三维表面参数确定图像产生问题的原因，并确定降雨图像；最后利用频域滤波的方式处理被降雨影响较轻的图像。

本发明提出的一种基于X波段雷达图像的降雨干扰抑制方法，包括下列步骤：

第一步，采集N幅空间海域杂波连续图像，然后选取分析区域构成图像序列。

第二步，对图像序列进行3维傅里叶变换得到三维波数-频率图像谱。

第三步，根据色散关系构造带通滤波器，滤除三维波数-频率图像谱中的非海浪信号，计算相对信号R_p，并根据相对信号R_p的值，确定存在问题的图像；具体为：相对信号R_p＞＝0.6时，雷达图像序列是可用的；相对信号R_p＜0.6时，雷达图像序列是存在问题的。

第四步，计算由上一步中确定出的存在问题的各组图像序列的三维表面参数。所述的三维表面参数包括回波强度的均值和表面纹理的视在特性比。

第五步，计算第三步中确定出的存在问题的各组图像序列的信噪比(SNR)：

带通滤波之后得到三维图像谱，经过积分得到二维图像谱，二维图像谱经过调制传递函数作用之后得到二维海浪谱，信噪比的计算公式定义如下：

$\mathrm{SNR}=\sqrt{\frac{{\mathrm{Eng}}_{\mathrm{wave}}}{{\mathrm{Eng}}_{\mathrm{total}}-{\mathrm{Eng}}_{\mathrm{image}}}}$

其中Eng_total代表的是滤波前的三维波数-频率图像谱积分后得到的总能量，Eng_wave代表二维海浪谱能量，Eng_image代表二维图像谱能量。

第六步，根据回波强度的均值，表面纹理的视在特性比以及信噪比，确定引起图像序列产生问题的具体原因，并确定降雨图像。

第七步，针对引起图像序列产生问题的不同原因进行相应的处理。

如果本身采集的图像序列存在缺失或者海浪回波太小等，直接舍弃该组图像序列；对于降雨影响严重的图像序列，也直接舍弃；对于降雨影响较轻的图像序列，进行频域滤波处理。

对于降雨图像是指降雨影响较轻的图像序列，其波谱能量分布中，低频部分存在过多的能量，与正常未受降雨影响的图像有明显的差异，设置适当的截断频率，进行频域滤波处理。

本发明的优点在于：

(1)本发明提出一种基于X波段雷达图像的降雨干扰抑制方法，与现有降雨干扰抑制算

法相比，引入质量控制，减少了需要进行识别判断的数据量，提高算法运行速度。

(2)本发明提出一种基于X波段雷达图像的降雨干扰抑制方法，与现有降雨干扰抑制算

法相比，引入三维参数描述图像的回波强度变化，它比二维参数能提供更多的信息，描述

更接近于真实表面，且能够给出所研究表面的直观图像以及与表面形态有关的足够信息。

(3)本发明提出一种基于X波段雷达图像的降雨干扰抑制方法，与现有降雨干扰抑制算法相比，增加了信噪比做参考，信噪比SNR是波高的直接反应，在波谱分析上可以利用波高的改变来了解能量变化，这样能更准确的确定雷达回波图像中受降雨影响的数据，以及影响的程度。

(4)本发明提出一种基于X波段雷达图像的降雨干扰抑制方法，与现有降雨干扰抑制算法相比，能更准确地识别降雨资料。

(5)本发明提出一种基于X波段雷达图像的降雨干扰抑制方法，与现有降雨干扰抑制算法相比，提高了海浪参数的反演精度。

附图说明

图1a未受降雨影响的图像，图1b部分受降雨影响的图像，图1c降雨影响严重的图像；

图2a相对信号，图2b信噪比以及有效波高，图2c回波强度均值，图2d表面纹理视在特性比；

图3a未受降雨影响的能谱图，图3b受降雨影响的能谱图；

图4a有效波高对比结果，图4b周期对比结果，图4c波向对比结果；其中“-*-filter”代表的是本发明滤波后的结果，“-.-original”代表的是未经滤波的结果，另一个为相应的WAMOS或者WAVEX提供的参考值；

图5为本发明的降雨干扰抑制方法的流程图。

具体实施方式

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

本发明公开的一种基于X波段雷达图像的降雨干扰抑制方法，如图5所示，具体步骤如下：

第一步，采集N幅空间域海杂波连续图像，然后选取分析区域构成图像序列。所述的图像序列的个数根据需要的反演精度进行设定，本发明中优选图像序列3组以上。

第二步，对各组图像序列进行3维傅里叶变换得到三维波数-频率图像谱；

第三步，根据色散关系构造带通滤波器，滤除三维波数-频率图像谱中的非海浪信号，计算相对信号R_p，并根据相对信号R_p的值，确定存在问题的图像。

回波图像中的降雨、雪，陆地或者目标物等会对海浪波谱分析产生影响。带通滤波会使得不属于海浪的噪声被大量的滤除，但不能完全滤除。那么滤波前后的能量比例就是一个判定图像质量好坏的依据。

经验带通滤波器是根据色散关系构造的，它基于一个最大期望流速(雷达天线与海浪场的最大相对流速)U_max，所述的带通滤波器满足：

$E(k_x,k_y,\mathrm{\ω})=\left\{\begin{array}{cc}I(k_x,k_y,\mathrm{\ω})& k\∈[B_n,B_p]\\ 0& \mathrm{otherwise}\end{array}\right.$

其中

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

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

式中，I(k_x，k_y，ω)是滤波前的三维波数-频率图像谱，E(k_x，k_y，ω)为滤波后的三维波数-频率图像谱，Δω为频率分辨率，Δk为波数分辨率，ω为海浪频率，g为当地重力加速度，k为波数，B_n、B_p决定带通滤波器的带宽，并且假设B_n的定义式中，

非负。当

取值为负时，则令B_n＝0。

令波数-频率空间中滤波前的谱I(k_x，k_y，ω)能量为E_image，带通滤波后的谱E(k_x，k_y，ω)能量为E_wave，则有相对信号：

$R_p=\frac{E_{\mathrm{wave}}}{E_{\mathrm{image}}}$

相对信号R_p的阈值是一个质量控制的指标，一般认为相对信号R_p＞＝0.6时，雷达图像序列是可用的；相对信号R_p＜0.6时，雷达图像序列是存在问题的，这里面就包含了各种噪声的影响，不能直接进行波谱分析和海浪参数的反演计算等。下面步骤中将对存在问题的雷达图像序列进行处理。

第四步，计算存在问题的各组图像序列的三维表面参数。所述的三维表面参数包括回波强度的均值ave，表面纹理的视在特性比S_tr和表面的均方根偏差。

所述的回波强度的均值ave为：

$\mathrm{ave}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i$

其中x_i表示雷达图像上第i个像素点的回波强度值，n表示每幅图像上像素点的个数。

设M、N分别为在图像序列中每一幅图像内x向和y向的采样点数；z(x，y)为表面的偏离高度。

第一类参数基本上是二维相关参数的直接扩展，数学公式表示如下：

A.S_q——表面的均方根偏差：

$S_q=\sqrt{\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left[z(x_i,y_j)\right]}^2}$

其中，z(x_i，y_j)为图像内点(x_i，y_j)的表面偏离高度，i＝1，2，...M，j＝1，2，...N。

第二类参数描述的是表面水平方向上的粗糙度，定义如下：

定义标准化的区域自相关函数(AACF)：

$\mathrm{AACF}({\mathrm{\τ}}_x,{\mathrm{\τ}}_y)=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}z(x_i,y_j)z(x_i-{\mathrm{\τ}}_x,y_j-{\mathrm{\τ}}_y)}{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{z}^2(x_i,y_j)}$

其中，τ_x和τ_y是区域自相关函数AACF的每一点的位置指标，1-M≤τ_x≤M-1；1-M≤τ_y≤M-1，τ_x，τ_y∈Z(整数集)，z(x_i，y_j)和z(x_i-τ_x，y_j-τ_y)分别为表面偏离高度z(x，y)在点(x_i，y_j)和点(x_i-τ_x，y_j-τ_y)处的值。

B.最速下降自相关长度S_fd：是区域自相关函数AACF最快衰减到0.2的距离，

$S_{\mathrm{fd}}=\min \left(\sqrt{{\mathrm{\τ}}_x^2+{\mathrm{\τ}}_y^2}\right),\mathrm{AACF}({\mathrm{\τ}}_x,{\mathrm{\τ}}_y)\≤0.2$

相应的最慢下降自相关长度S_sd，是AACF最慢衰减到0.2的距离，

$S_{\mathrm{sd}}=\max \left(\sqrt{{\mathrm{\τ}}_x^2+{\mathrm{\τ}}_y^2}\right),\mathrm{AACF}({\mathrm{\τ}}_x,{\mathrm{\τ}}_y)\≤0.2$

C.表面纹理的视在特性比S_tr：

$S_{\mathrm{tr}}=\frac{S_{\mathrm{fd}}}{S_{\mathrm{sd}}}$

视在特性比S_tr是用于描述表面纹理特性的参数。显然，0＜S_tr＜1，对于表面有很强纹理方向的表面，S_tr＜＜1；然而对于纹理方向性不强或者各向异性的表面，S_tr的值较大并且趋于1。

第五步，计算第三步中确定出的存在问题的各组图像序列的信噪比(SNR)：

带通滤波之后得到滤除了非海浪信号的三维波数-频率图像谱，经过积分得到二维图像谱I(k_x，k_y)，二维图像谱经过调制传递函数作用之后得到二维海浪谱E(k_x，k_y)。

经验的调制传递函数由下面公式给出：

M(k)≈k^-β

其中β为经验参数，一般取值为1.2；k为波数。

相应的二维海浪谱E(k_x，k_y)定义如下：

E(k_x，k_y)＝M(k)·I(k_x，k_y)

信噪比的计算公式定义如下：

$\mathrm{SNR}=\sqrt{\frac{{\mathrm{Eng}}_{\mathrm{wave}}}{{\mathrm{Eng}}_{\mathrm{total}}-{\mathrm{Eng}}_{\mathrm{image}}}}$

其中Eng_total代表的是三维波数-频率图像谱积分后得到的总能量，Eng_wave代表二维海浪谱能量，Eng_inage代表二维图像谱能量。

因为在描述波浪特性的参数中，波高是表示波浪特性的重要参数之一，在波谱分析上可以利用波高的改变来了解能量变化，在工程设计和海洋相关问题的讨论上，波高更是重要的决策依据。而信噪比是波高的直接反应，所以后续的讨论将以信噪比作为讨论依据。

第六步，根据回波强度的均值ave，表面纹理的视在特性比S_tr以及信噪比SNR变化趋势确定引起图像序列产生问题的具体原因。

回波强度均值ave是回波强度的一个直观表现，它的变化反映了整幅图像回波强度的改变。表面纹理视在特性比反映了表面纹理的变化，通过它的大小可以判定表面纹理有无及纹理大小。

降雨、较小的海浪或者无浪时，回波图像都会比较平滑，此时表面纹理视在特性比S_tr的值就会比较大。如果有一定的海浪，存在降雨的影响，此时回波强度会比没有降雨的时候大，回波强度均值基本都偏大，回波强度中包含一定海浪的信息，经过滤波之后得到的信噪比会比浪小或者无浪的时候大；如果海浪小或者无浪时，回波强度均值基本都小，相应的信噪比也很小。具体的信噪比、回波强度均值和表面纹理视在特性比的大小的界定值要根据实际计算的结果设定。

根据各个参数的表征意义及具体数值，判定图像序列产生问题的原因。

在R_p＜0.6时，有：

SNR	ave	Str	原因
				大	大	大	降雨
小	小	大	海面波浪小或者基本没有海浪

在R_p＞＝0.6时，有：

另外，在R_p＞＝0.6的时候，S_tr和均值有同增(或者同减)的趋势，并且S_tr的值超越0.4，这种情况判定为降雨不是影响该时刻的图像序列中的所有图像而是影响其中的几幅图像。

第七步，针对引起图像序列产生问题的不同原因进行相应的处理：

如果本身采集的图像存在缺失或者海浪回波太小等，直接舍弃该组图像序列；对于降雨影响严重的图像序列，也直接舍弃；对于降雨影响较轻的图像序列，进行频域滤波处理。

降雨影响较轻的图像序列为降雨图像，其波谱能量分布中，低频部分存在过多的能量，与正常未受降雨影响的图像有明显的差异，设置适当的截断频率，对降雨图像进行频域滤波处理。

应用本发明提出的一种基于X波段雷达图像的降雨干扰抑制方法，采用海表面流反演方法进行实地实测，实验地点在水文观测站附近海域，实验数据选取时间为2010年10月24-25日，共计347组，其中24日的数据根据人工观测，有大部分降雨的影响。X波段雷达的架设高度为40米，平均旋转周期为2.39秒，工作在短脉冲模式，雷达探测半径约为2千米海域。每3分钟采集一组数据，每组包含32幅图像，使用现场WAMOS雷达提供的有效波高作波高的参考值，WAMOS每2分钟输出一次结果；以WAVEX雷达提供的波周期和波向作周期和波向的参考值，WAVEX每4分钟输出一次结果。实验图像区域中心位于船艏75度方向，近点离岸600米，面积960×960m²，此海域平均水深20米。

图1a～1c分别为未受降雨影响的图像，部分受降雨影响的图像以及被降雨严重模糊的图像，从这图中可以看出，降雨会模糊原来的海浪纹理，甚至会使得海浪纹理全部没有，这样必然后影响后续的海浪参数反演。

为了更直观的对比说明，本发明也计算了整段时间内数据的回波强度(取值范围为0-8191)的均值ave，相对信号R_p、表面纹理的视在特性比S_tr以及信噪比SNR，如图2a～2d所示。降雨、较小的海浪或者无浪时，回波图像都会比较平滑，此时表面纹理视在特性比S_tr的值就会比较大(基本都超过0.4)。如果有一定的海浪，存在降雨的影响，此时回波强度会比没有降雨的时候大，回波强度均值基本都大于2000，回波中包含一定海浪的信息，经过滤波之后得到的信噪比会比浪小或者无浪的时候大(基本都大于1)；如果海浪小或者无浪时，回波均值基本都小于1500，相应的信噪比也很小(基本在0.5附近)。本实施例中，设定信噪比、回波强度均值和表面纹理视在特性比的大小的界定值分别为：表面纹理视在特性比大于0.4认为是大，否则为小；回波强度均值大于2000认为大，否则为小；信噪比大于1认为大，否则为小。

在图3a、3b中给出了降雨和未降雨情况波谱能量的曲线，从中可以看出在低频率段，降雨能量曲线出现大能量分布。由于海浪的实际周期并不会出现在那么大的范围内，所以频域滤波过程中将频率小于0.06(对应周期15s)的能量清零了。

为验证本发明提出的基于X波段雷达图像的降雨干扰抑制方法的性能，将不作处理得到的结果和本发明得到的结果分别与现场传雷达提供的真实值进行比对，本发明实验中参数设置如表1所示。

表1：实验参数设置

参数	数值
		频率分辨率Δω	0.082s-1
波数分辨率Δk	0.0065m-1

空间x向分辨率Δx	7.5m
		空间y向分辨率Δy	7.5m
实验海区平均水深d	20m
		非线性能量校正指数β	1.2
图像谱中正频率数目Nω	15
		最大期望流速Umax	3m/s

为清晰的看出改进算法本发明的性能，图4a～4c中分别给出了对于降雨影响的数据不经滤波得到的反演结果、本发明得到反演结果与参考值进行了对比(包括有效波高(Hs)，周期(Tp)，波向(Pdir))，误差统计结果如表2所示。

设X_i为算法反演输出，即测量值，其中i＝1，2，...N；Y_i为真值(WAMOS或者WAVEX提供)，其中i＝1，2，...N，则定义：

测量值均值为： $\stackrel{\‾}{X}=\frac{{\mathrm{\Σ}}_{i=1}^N{X}_i}{N};$

真值均值为： $\stackrel{\‾}{Y}=\frac{{\mathrm{\Σ}}_{i=1}^N{Y}_i}{N};$

序列偏差为：A_i＝Y_i-X；

平均差为：

均方差为： $\mathrm{VAR}=\frac{{\mathrm{\Σ}}_{i=1}^N{A}_i^2}{N-1};$

标准均方根误差为：

测量误差为：

则相对误差为：

对波向和有效波高而言，不需要计算相对误差。

表2误差统计结果

从表2可以看出，针对降雨图像做的滤波处理，周期和波向都得到了很大程度的改善；有效波高的变化不大，误差都在允许范围之内。

Claims (5)

1.一种基于X波段雷达图像的降雨干扰抑制方法，其特征在于：

第一步，采集N幅空间域海杂波连续图像，然后选取分析区域构成图像序列；

第二步，对各组图像序列进行三维傅里叶变换得到三维波数-频率图像谱；

第三步，根据色散关系构造带通滤波器，滤除三维波数-频率图像谱中的非海浪信号，计算相对信号R_p，并根据相对信号R_p的值，确定存在问题的图像序列，具体为：相对信号R_p≥0.6时，雷达图像序列是可用的；相对信号R_p＜0.6时，雷达图像序列是存在问题的；

第四步，计算存在问题的各组图像序列的三维表面参数；所述的三维表面参数包括回波强度的均值和表面纹理的视在特性比；

第五步，计算第三步中确定出的存在问题的各组图像序列的信噪比；

第六步，根据回波强度的均值，表面纹理的视在特性比以及信噪比，确定引起图像序列产生问题的具体原因，并确定降雨图像；

第七步，针对引起图像序列产生问题的不同原因进行相应的处理。

2.根据权利要求1所述的一种基于X波段雷达图像的降雨干扰抑制方法，其特征在于：第三步中所述的带通滤波器满足：

其中，

式中，I(k_x，k_y，ω)是滤波前的三维波数-频率图像谱，E(k_x，k_y，ω)为滤波后的三维波数-频率图像谱，Δω为频率分辨率，Δk为波数分辨率，ω为海浪频率，g为当地重力加速度，k为波数，B_n、B_p为带通滤波器的带宽，并且假设B_n的定义式中， 

非负；当 取值为负时，则令B_n＝0；

令波数-频率空间中滤波前的谱I(k_x，k_y，ω)能量为E_image，带通滤波后的谱E(k_x，k_y，ω)能量为E_wave，则有相对信号： 

。

3.根据权利要求1所述的一种基于X波段雷达图像的降雨干扰抑制方法，其特征在于：第四步所述的回波强度的均值为：

其中x_i表示雷达图像上第i个像素点的回波强度值，n表示每幅图像上像素点的个数。

4.根据权利要求1所述的一种基于X波段雷达图像的降雨干扰抑制方法，其特征在于：第四步所述的表面纹理的视在特性比S_tr为：

其中，S_fd为最速下降自相关长度，是区域自相关函数AACF最快衰减到0.2的距离，S_sd是最慢下降自相关长度，是AACF最慢衰减到0.2的距离，具体为：

。

5.根据权利要求1所述的一种基于X波段雷达图像的降雨干扰抑制方法，其特征在于：第四步所述的信噪比的计算公式定义如下：

其中Eng_total代表的是三维波数-频率图像谱积分后得到的总能量，Eng_wave代表二维海浪谱能量，Eng_image代表二维图像谱能量。 

Images