一种基于X波段导航雷达图像的有效波高反演方法

一种基于x波段导航雷达图像的有效波高反演方法
技术领域
1.本发明涉及图像反演技术领域，具体涉及一种基于x波段导航雷达图像的有效波高反演方法。


背景技术：

2.利用导航雷达可以反演海浪的有效波高、波峰峰向和主波周期等参数。目前，通过导航雷达反演有效波高的算法有谱分析法、阴影法和正交分解法等，但在工程实践中应用的主要是基于三维傅里叶变换的谱分析法。
3.谱分析法是选取雷达图像序列中的特定区域进行三维傅里叶变换得到图像谱，然后基于色散带通滤波器提取出海浪谱，进而获取海浪信号的信噪比，最后利用有效波高与信噪比均方根之间的线性经验关系来估测有效波高。
4.在实际应用中，有效波高与信噪比的均方根之间并不是完全呈线性关系，而且这种基于线性关系表示的方法需要外部专用设备如浮标来进行标定。此外，需要大量有效观测数据进行标定，才能使反演的有效波高满足精度要求，从而大大增加研发成本和周期。更重要的是，由于雷达系统的差异、海域环境的变化、信噪比计算方法不同等因素的影响，特别是在近岸海域海底地形复杂，以及海岸线、岛屿等对海浪的反射、衍射等影响，对于每部雷达以及不同海域都需要重新标定系数。因此，在实际应用中严重影响了谱分析法的广泛使用。
5.同时，阴影统计法由于具有无需外部参考设备进行标定的优势，但是，在利用阴影统计法反演波高过程中需使用波周期，而波周期通常由外部参考设备如浮标提供。此外，在利用波高、波陡和波周期之间的关系反演有效波高的过程中，通常假定海浪满足一阶色散关系且忽略了水深、海表层流对波高反演的影响，与实际应用的情况不符，不能使反演的有效波高满足精度的要求。


技术实现要素：

6.针对现有技术存在的不足，本发明提出一种基于x波段导航雷达图像的有效波高反演方法，以解决现有波高反演方法需要外部参考设备进行标定、水深和海表层流等因素对波高反演的影响，导致反演波高的精度低的问题。
7.为达到上述目的，本发明采用如下技术方案：一种基于x波段导航雷达图像的有效波高反演方法，其特征在于，包括以下步骤：
8.采集原始导航雷达图像，根据原始导航雷达图像分别选取基于谱分析技术的分析区域和阴影统计法的分析区域；
9.对选取的所述基于阴影统计法的分析区域进行边缘检测，得出雷达图像的阴影分割阈值；
10.利用雷达图像的阴影分割阈值获取阴影图像，进而计算分析区域中的阴影比例，得出每个分区的阴影比例函数；
11.根据史密斯函数和所述阴影比例函数进行均方根波陡的计算，得到均方根波陡；
12.对选取的所述基于谱分析技术的分析区域图像序列进行三维傅立叶变换后得到波数频率谱，利用色散带通滤波器对其进行滤波处理，再通过调制传递函数修正后得到海浪谱，提取主波波数；
13.根据所述均方根波陡和所述主波波数，得出有效波高。
14.优选的，选取所述基于谱分析技术的分析区域时，选取的分析区域至少包含3个完整的海浪波形。
15.优选的，对选取所述基于阴影统计法的雷达图像分析区域进行边缘检测的具体步骤为：
16.对雷达图像i(r,θ)在相邻的8个方向上分别基于差分算子进行卷积运算，得到不同方向上的梯度图像ii；
17.选取前n％对应的所述雷达图像的像元值作为阈值ti，对得到的梯度图像ii进行边缘提取；
18.边缘提取的计算公式为：
[0019][0020]
通过将八个不同方向上提取的边缘图像进行叠加，得到完整的边缘图像；
[0021]
完整的边缘图像的计算公式为：
[0022][0023]
优选的，所述雷达图像阴影分割阈值的具体计算步骤为：
[0024]
将所述边缘图像i
t
的每一个孤立的1，找到在该位置处原始导航雷达图像的像元值η，通过边缘图像与原始导航雷达图像之间的对应关系建立一个图像直方图函数fh(η)，对fh(η)取众数，得到阴影分割阈值τs；
[0025]
阴影分割阈值τs计算公式为：
[0026]
τs＝mod e(fh(η))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)。
[0027]
优选的，得出所述每个分区的阴影比例函数的具体步骤为：
[0028]
利用所述阴影阈值τs将原始导航雷达图像分割成阴影区和非阴影区，阴影图像is的计算公式为：
[0029][0030]
将所述阴影图像在方位向上分割成区，并将阴影图像按35～50m的径向步长平均分割为固定的块数；
[0031]
通过计算每个分割块中的阴影比例，得到每个分区的阴影比例函数s(θ)。
[0032]
优选的，得到所述均方根波陡的具体步骤为：
[0033]
将所述边缘图像分割区的每一块阴影比例函数s(θ)拟合估计，得到每个分区的均方根波陡σ
rms
；
[0034]
将方位向上每个分区中的波陡进行平均，得到平均均方根波陡，其计算公式如下：
[0035][0036]
其中，m为阴影图像分割成的区的总数。
[0037]
优选的，所述雷达图像序列的波数频率谱的具体计算步骤为：
[0038]
将所述选取的分析区域雷达图像序列η(x,y,t)进行三维傅里叶变换得到
[0039]
三维波数频率图像谱f(k
x
,ky,ω)，其计算公式为：
[0040][0041]
其中，l
x
*ly为选取的分析区域，雷达图像序列总时间为t；为波数，ω为角频率；三维波数频率能量谱的表达式为：
[0042][0043]
优选的，利用所述色散带通滤波器对其进行滤波处理的具体步骤为：
[0044]
对所述雷达图像序列进行三维傅里叶变换后，得到三维波数频率图像谱，假定在选定的区域内，海表面的浪场在空间上匀质和在时间上平稳；
[0045]
在线性波理论基础下，海浪满足的色散关系为：
[0046][0047]
其中，g为重力加速度，为波数，λ为海浪波长，d为水深；为相对于雷达天线平台运动的海表层流，
[0048]
所述通过调制传递函数修正后得到海浪谱的具体步骤为：
[0049]
对三维波数频率能量谱i(k
x
,ky,ω)进行滤波处理后，对频率ω进行积分得到二维图像谱i(k
x
,ky)，利用经验调制传递函数对二维图像谱i(k
x
,ky)进行修正：
[0050]
e(k
x
,ky)＝|m(k
x
,ky)|2·
i(k
x
,ky)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)；
[0051]
其中，|m(kx,ky)|2为经验调制传递函数，e(k
x
,ky)为二维海浪谱；调制传递函数为：
[0052]
|m(k
x
,ky)|2≈k-β
[0053]
其中，β＝1.2为经验系数。
[0054]
优选的，提取所述主波波数的具体步骤为：
[0055]
笛卡尔坐标下的海浪谱e(k
x
,ky)与极坐标下的海浪谱e(k,θ)的关系为：
[0056]
e(k,θ)＝e(k
x
,ky)k
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)；
[0057]
其中：
[0058][0059]
arctan2(
·
)为反正切函数；寻找海浪谱e(k,θ)中能量最大值对应的位置，即为主波波数km。
[0060]
优选的，得出所述有效波高的具体计算步骤为：
[0061]
有效波高hs与均方根波陡和主波波数km的关系为：
[0062][0063]
其中，c为经验系数。
[0064]
与现有技术相比，本方案产生的有益效果是：
[0065]
1、本方案与现有的方法相比具有无需外部参考设备进行标定的优势，能够完成有效波高反演的任务，且提高了阴影统计法反演波高的精度。
[0066]
2、本方案不需要大量有效的观测数据进行标定，能够降低研发成本和缩短研发周期。
[0067]
3、本方案在反演波高过程中无需外部参考设备提供波周期，降低了系统复杂性，对于每部雷达以及不同海域无需标定系数。
[0068]
4、本方案基于谱分析技术对雷达图像序列进行三维傅里叶变换得到波数频率图像谱后，利用色散关系带通滤波器提取海浪成份，在滤波过程中同时考虑了水深、海表层流等因素对波高反演的影响，能够提高有效波高的反演精度。
附图说明
[0069]
为了更清楚地说明本发明具体实施方式，下面将对具体实施方式中所需要使用的附图作简单地介绍。在所有附图中，各元件或部分并不一定按照实际的比例绘制。
[0070]
图1为本发明一种基于x波段导航雷达图像的有效波高反演方法的流程图；
[0071]
图2为本发明一种基于x波段导航雷达图像的有效波高反演方法中采集的原始导航雷达图像；
[0072]
图3为本发明一种基于x波段导航雷达图像的有效波高反演方法中选取的基于谱分析技术的分析区域和阴影统计法的的分析区域；
[0073]
图4为本发明一种基于x波段导航雷达图像的有效波高反演方法中对雷达图像中分析区域进行卷积运算后得到的梯度图像；
[0074]
图5为本发明一种基于x波段导航雷达图像的有效波高反演方法中利用阴影阈值分割雷达图像后得到的阴影图像；
[0075]
图6为本发明一种基于x波段导航雷达图像的有效波高反演方法中本发明中基于谱分析技术得到的二维海浪谱；
[0076]
图7为本发明一种基于x波段导航雷达图像的有效波高反演方法中为浮标记录的有效波高和基于雷达图像反演的有效波高的结果图；
[0077]
图8为本发明一种基于x波段导航雷达图像的有效波高反演方法中为现有的阴影统计法反演的有效波高的散点图；
[0078]
图9为本发明一种基于x波段导航雷达图像的有效波高反演方法中为本发明中提出方法反演的有效波高的散点图。
具体实施方式
[0079]
下面将结合附图对本发明技术方案的实施例进行详细的描述。以下实施例仅用于更加清楚地说明本发明的技术方案，因此只作为示例，而不能以此来限制本发明的保护范
围。
[0080]
请参阅图1-图9，一种基于x波段导航雷达图像的有效波高反演方法，包括以下步骤：
[0081]
步骤一、请一并参阅图2，采集原始导航雷达图像。图2中图像右半部分为对海有效观测区域，期间海浪的有效波高为1.8m。
[0082]
根据原始导航雷达图像选取基于谱分析技术的分析区域和阴影统计法的分析区域。
[0083]
请一并参阅图3，针对阴影统计法从雷达图像中选取分析区域，在距离方向上，选取400～2500m范围内的雷达回波信号的用于有效波高反演。
[0084]
针对谱分析法从雷达图像序列中选取分析区域，根据掠射角区间的最佳范围，选取的分析区域距离雷达平台实际距离约为600m，此时导航雷达回波图像分辨率比较好，图像中的海浪纹理特征明显。
[0085]
本步骤中，雷达图像中选取基于谱分析技术的分析区域时，从雷达图像的序列中选取的分析区域包含3个完整的海浪波形。同时，也可选用从雷达图像的序列中选取的分析区域包含4个或5个完整的海浪波形。
[0086]
本发明中采集的雷达图像距离分辨率为7.5m，为便于后续的傅里叶变换分析在空间上选取128*128的像素点，即分析区域为960m*960m。实验中导航雷达的旋转周期约为2.3s，一般选取40～100s的雷达图像时间序列，为了方便在时间域上进行傅里叶变换，选取连续的32幅雷达海杂波图像进行分析，时间总长度约为80s。
[0087]
可知，基于三维傅里叶变换法进行海浪信息反演，从32幅连续的航海雷达图像选取距离船艏角度固定，离雷达平台的距离约为600m，大小为960m*960m的分析区域。图3中的黑色矩形框为从雷达图像中选取的分析区域用于提取主波波数。
[0088]
步骤二、请一并参阅图4，对选取的阴影统计法的雷达图像分析区域进行边缘检测，得出雷达图像的阴影分割阈值。
[0089]
对选取的阴影统计法的雷达图像分析区域进行边缘检测的具体步骤为：
[0090]
对雷达图像i(r,θ)在相邻的8个方向上分别基于差分算子进行卷积运算，得到8个不同方向上的梯度图像ii。
[0091]
选取前n％对应的雷达图像的像元值作为阈值ti，对得到的梯度图像ii进行边缘提取。本步骤中n取10，即取前10％的像素值作为梯度阈值对边缘梯度进行图像边缘提取，通过将滤波处理后的8个方向的边缘梯度图像进行叠加得到一幅梯度图像，如图4所示，从图中能够看到海杂波条纹的边界。
[0092]
边缘提取的计算公式为：
[0093][0094]
通过将八个不同方向上提取的边缘图像进行叠加，得到完整的边缘图像。
[0095]
完整的边缘图像的计算公式为：
[0096][0097]
同时，请一并参阅图5，雷达图像阴影分割阈值的具体计算步骤为：
[0098]
将边缘图像i
t
的每一个孤立的1，找到在该位置处原始导航雷达图像的像元值η，通过边缘图像与原始导航雷达图像之间的对应关系建立一个图像直方图函数fh(η)，对fh(η)取众数，得到阴影分割阈值τs。
[0099]
阴影分割阈值τs计算公式为：
[0100]
τs＝mod e(fh(η))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)。
[0101]
对获得的阴影边缘梯度图像进行直方图统计，并对直方图统计结果取众数得到阴影分割阈值τs。
[0102]
具体的，图5中白色的部分表示被海浪遮挡住的阴影区域，黑色表示能够被雷达观测到的海浪，从图中可以清晰的观察到海浪的纹理特征，同时随着与雷达天线距离的增加，图像中阴影也在逐渐增多。
[0103]
步骤三、利用雷达图像的阴影分割阈值获取阴影图像，进而计算分析区域中的阴影比例，得出每个分区的阴影比例函数。
[0104]
得出每个分区的阴影比例函数的具体步骤为：
[0105]
利用阴影阈值τs将原始导航雷达图像分割成阴影区和非阴影区，阴影图像is的计算公式为：
[0106][0107]
将阴影图像在方位向上分割成区，并将阴影图像按35～50m的径向步长平均分割为固定的块数。
[0108]
具体的，即本步骤中先将雷达阴影图像在方位角度向上每隔10
°
左右划分为一个分区，再将阴影图像按35～50m的径向步长平均分割为固定的50块左右。
[0109]
通过计算每个分割块中的阴影比例，得到每个分区的阴影比例函数s(θ)。即在每10
°
一个的分区可以得到阴影比例函数，由阴影比例函数能够得到照度曲线。
[0110]
步骤四、根据史密斯函数和所述阴影比例函数进行均方根波陡的计算，得到均方根波陡。
[0111]
得到均方根波陡的具体步骤为：
[0112]
将边缘图像分割区的每一块阴影比例函数s(θ)拟合估计，得到每个分区的均方根波陡σ
rms
。通过利用史密斯拟合函数对进行照度曲线拟合求取每个分区上的均方根波陡σ
rms
。通过对方位向上每个分区中的波陡进行平均能够得到平均均方根波陡
[0113]
将方位向上每个分区中的波陡进行平均，得到平均均方根波陡，其计算公式如下：
[0114][0115]
其中，m为阴影图像分割成的区的总数。
[0116]
步骤五、请一并参阅图6，采用三维傅立叶方法计算雷达图像序列的波数频率谱，利用色散带通滤波器对其进行滤波处理，再通过调制传递函数修正后得到海浪谱，提取主波波数。
[0117]
雷达图像序列的波数频率谱的具体计算步骤为：
[0118]
将选取的分析区域雷达图像序列η(x,y,t)进行三维傅里叶变换得到三维
[0119]
波数频率图像谱f(k
x
,ky,ω)，其计算公式为：
[0120][0121]
其中，l
x
*ly为选取的分析区域，雷达图像序列总时间为t；为波数，ω为角频率；三维波数频率能量谱的表达式为：
[0122][0123]
利用色散带通滤波器对其进行滤波处理的具体步骤为：
[0124]
对雷达图像序列η(x,y,t)进行三维傅里叶变换后，得到三维波数频率图像谱。该能量谱中不仅包含海浪能量，还包含了大量不属于海浪的噪声，利用色散关系带通滤波器滤除图像谱中的非海浪成份。假定在选定的区域内，海表面的浪场在空间上匀质和在时间上平稳；
[0125]
在线性波理论基础下，海浪满足的色散关系公式为：
[0126][0127]
其中，g为重力加速度，为波数，λ为海浪波长，d为水深。为相对于雷达天线平台运动的海表层流，
[0128]
通过调制传递函数修正后得到海浪谱的具体步骤为：
[0129]
对三维波数频率能量谱i(k
x
,ky,ω)进行滤波处理后，对频率ω进行积分得到二维图像谱i(k
x
,ky)；
[0130]
由于雷达成像机制的非线性特性，二维图像谱i(k
x
,ky)与真实的二维海浪谱之间存在一定的差异。利用经验调制传递函数对图像谱i(k
x
,ky)进行修正：
[0131]
e(k
x
,ky)＝|m(k
x
,ky)|2·
i(k
x
,ky)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)；
[0132]
其中，|m(k
x
,ky)|2为经验调制传递函数，e(k
x
,ky)为二维海浪谱；调制传递函数为：
[0133]
|m(k
x
,ky)|2≈k-β
[0134]
其中，β＝1.2为经验系数。
[0135]
提取主波波数的具体步骤为：
[0136]
笛卡尔坐标下的海浪谱e(k
x
,ky)与极坐标下的海浪谱e(k,θ)的关系为：
[0137]
e(k,θ)＝e(k
x
,ky)k
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)。
[0138]
其中：
[0139][0140]
arctan2(
·
)为反正切函数；寻找海浪谱e(k,θ)中能量最大值对应的位置，即为主波波数km。
[0141]
步骤六、根据均方根波陡和主波波数，得出有效波高。
[0142]
具体的，有效波高的具体计算步骤为：
[0143]
有效波高hs与均方根波陡和主波波数km的关系为：
[0144][0145]
其中，c为经验系数。
[0146]
实验分析
[0147]
实验期间布放的波浪浮标作为参考设备。
[0148]
为验证本发明的有效性，请一并参阅图7、图8和图9，将本发明的反演结果与现有的阴影统计法的反演结果进行比较，现有的阴影统计法的反演结果和性能分析。
[0149]
其中，图7为现有的阴影统计法反演的有效波高和本发明中的方法反演的有效波高、浮标记录的有效波高的对比图。叉号表示现有的阴影统计法的反演结果，三角形表示本发明提出方法的反演结果，圆圈表示浮标记录。从图7可以看出，与现有的阴影统计法求得的有效波高相比，利用本发明提出方法计算得到的有效波高更接近浮标测量值，而基于现有的阴影统计法反演得到的有效波高与真值曲线偏差较大。本发明提出的方法能够完成有效波高的反演任务。
[0150]
图8为现有的阴影统计法反演有效波高与浮标测量值的散点图。图9为基于本发明提出的方法反演的有效波高分别与浮标测量值的散点图。从图8可以看出现有的阴影统计法的散点分布相对较散，而本发明提出的方法的散点分布相对集中。本发明提出的方法测得的波高与外部参考波高的相关系数更大，且本求取的有效波高绝对误差更小，验证了本发明提出的方法比现有的阴影统计法更有优势，证明了本发明能够提高有效波高的反演精度。
[0151]
本发明提出的基于阴影统计法和谱分析法相结合从x波段导航雷达图像反演有效波高的方法。首先，从单幅雷达图像选取的分析区域，利用现有的阴影统计法从分析区域中提取均方根波陡；其次，从雷达图像序列选取的分析区域，利用谱分析法进从分析区域中提取二维海浪谱，进而获取主波波数；最后，根据均方根波陡、主波波数与有效波高间的关系，求取有效波高。本发明提出的方法的相关性系数和均方根方差分别为0.65和0.35，相较于现有的阴影统计法相关性系数提高了0.12，均方根方差降低了0.05。本发明提出的方法与现有的阴影统计法相比，反演精度的更高，提高了方法的适用性，使得反演精度完全达到工程应用需求。
[0152]
以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。