一种基于航海雷达图像的海面风速反演方法

摘要

一种基于航海雷达图像的海面风速反演方法，它属于海洋遥感技术领域。本发明解决了现有方法的风速反演精度低的问题。本发明的技术方案为：步骤1、离线读取雷达文件得到原始雷达图像集，分别对每个原始雷达图像进行处理获得每个原始雷达图像对应的逆风方向全距离积分面积；再根据采集每个原始雷达图像时的风速，建立风速与逆风方向全距离积分面积的函数关系；步骤2、在线读取风速待测的雷达图像，计算出风速待测的雷达图像的逆风方向全距离积分面积；步骤3、将步骤2中计算出的逆风方向全距离积分面积代入步骤1中建立的风速与逆风方向全距离积分面积的函数，得到采集步骤2读取的雷达图像时的海表面风速。本发明方法可以应用于海面风速反演。

说明书

技术领域

本发明属于海洋遥感技术领域，具体涉及一种基于航海雷达图像的海面风速反演方法。

背景技术

海面风场是研究海气之间互相作用的重要参数，海面风场通过调节热量、水汽、海气通量和颗粒物，调节大气和海洋之间的耦合作用，从而维持全球和区域气候。目前海面风速主要通过位于船的桅杆上的现场传感器(如风速计)获取，但风速计容易受大气湍流以及船体上层结构引起的气流紊乱的影响。即使是风速计被安装在无遮挡的位置，风参数测量的误差也可能很高(Thornhill E,Wall A,Mctavish S,et al.Ship anemometer biasmanagement[J].Ocean Engineering,2020,216:107843.)。近几十年来，X波段航海雷达以其高分辨率和及时反馈的优势，成为获取海表面风场信息的有效方法之一(王慧,卢志忠.基于波数能量谱的海面风向反演算法[J].华中科技大学学报(自然科学版),2014,42(12):96-100.)。并且，由于大部分船只都已经装备了X波段航海雷达，相对于其他遥感仪器，具有较小的额外成本，这些研究为航海雷达提取海面风场信息奠定了基础。

目前，有两类方法可用于从X波段航海雷达图像序列中反演风速。第一类方法是通过建立一个参考风速与雷达图像中后向散射强度信息相关的经验模型函数来反演风速。将雷达图像背景噪声的总方差与参考风速进行比较，建立了两者间的线性函数(Hatten H,Seemann J,Horstmann J,et al.Azimuthal dependence of the radar cross sectionand the spectral background noise of a nautical radar at grazing incidence[C]//IGARSS'98.Sensing and Managing the Environment.1998IEEE InternationalGeoscience and Remote Sensing.Symposium Proceedings.(Cat.No.98CH36174).IEEE,1998,5:2490-2492.以及Izquierdo P,Soares C G.Analysis of sea waves and windfrom X-band radar[J].Ocean Engineering,2005,32(11-12):1404-1419.)。2012年，Lund通过对海洋雷达图像的研究，发现回波强度与参考风速之间存在三次多项式关系，并通过最小二乘拟合建立了两者之间的三次多项式函数来反演风速(Lund B,Graber H C,Romeiser R.Wind Retrieval From Shipborne Nautical X-Band Radar Data[J].IEEETransactions on Geoscience&Remote Sensing,2012,50(10):3800-3811.)。随后，2013年Vicen-Bueno从三阶多项连续函数中反演风速，该函数取决于预先选定的强度值的最大范围距离(Vicen-Bueno R,Horstmann J,Terril E,et al.Real-time ocean wind vectorretrieval from marine radar image sequences acquired at grazing angle[J].Journal of Atmospheric and Oceanic Technology,2013,30(1):127-139.)。对于Decca和Furono雷达，2014年Liu又提出了利用双曲拟合法进行二次拟合，利用实测海洋雷达数据提取海面风速信息(Liu Y,Huang W,Gill E W,et al.Comparison of algorithms forwind parameters extraction from shipborne X-band marine radar images[J].IEEEJournal of Selected Topics in Applied Earth Observations and Remote Sensing,2014,8(2):896-906.)。2015年，Chen等人对于岸基雷达数据提出了一个风速经验函数，首先利用概率分布函数法建立雷达散射截面积法与风速的关系，然后在此基础上加入有效波高参数，使反演精度有所提高，但该模型只是适用于风浪情况(Chen Z,He Y,Zhang B,etal.Determination of nearshore sea surface wind vector from marine X-bandradar images[J].Ocean Engineering,2015,96(mar.1):79-85.)。另外，对于不同的雷达模型，需要对海面风和RCS进行建模。因此，模型函数法并不是普遍适用的。另一类从X波段海洋雷达图像序列中得出风速的方法是使用神经网络(Dankert H,Horstmann J,Rosenthal W.Ocean wind fields retrieved from radar-image sequences[J].Journalof Geophysical Research:Oceans,2003,108(C11).以及Dankert H,Horstmann J.Amarine radar wind sensor[J].Journal of Atmospheric and Oceanic Technology,2007,24(9):1629-1642.)，通过对风速和归一化雷达截面之间的关系进行参数化，并加入了海况信息和大气参数作为输入来获取模型。然而，神经网络在物理上是不明确的，因此它不能为雷达成像机制给出一个明确的模型。此外，神经网络的学习过程收敛缓慢，隐含层和隐含层节点数难以确定，这使得网络不适用于某些雷达数据，难以满足工程要求。

现有利用航海雷达图像进行风速反演方法普遍存在对雷达数据要求较高的缺点，需要360°无遮挡的雷达图像，这样才能获取到与风速具有相关性的雷达后向散射信息。另外，雷达后向散射信息还会受到固定目标的影响，固定物目标会影响到所选区域的后向散射信息，从而降低风速反演精度。

此外，雷达回波强度会受到距离的调制，回波强度与海洋表面的距离依赖可以呈现在雷达图像中，雷达回波强度会随着距离的增大而衰减。较远的海面回波信号不仅会随着距离的增大不断减小，还会出现较多的零强度值像素点，这会造成远距离的后向散射信息无法使用，但在高风速海况下，远距离的后向散射信息较为敏感，而低风速时对近距离的后向散射信息较为敏感，因此需要考虑所有距离的后向散射信息。有研究表明，雷达图像中发生的任何径向距离依赖都是海洋表面回波强度的掠射角依赖的结果，回波强度与径向距离为三次方衰减关系(Lund B,Graber H C,Romeiser R.Wind Retrieval From ShipborneNautical X-Band Radar Data[J].IEEE Transactions on Geoscience&Remote Sensing,2012,50(10):3800-3811.)。因此，当从航海雷达图像中反演风的信息时，不能忽略距离衰减的影响。

综上所述，由于实际雷达数据中存在遮挡、雷达后向散射信息受固定目标影响且雷达回波强度受距离的调制，因此，采用现有方法进行风速反演时仍然存在着风速反演精度低的问题。

发明内容

本发明的目的是为解决现有方法的风速反演精度低的问题，而提出的一种基于航海雷达图像的海面风速反演方法。

本发明为解决上述技术问题所采取的技术方案是：

一种基于航海雷达图像的海面风速反演方法，所述方法具体包括以下步骤：

步骤1、离线读取雷达文件得到原始雷达图像集，分别对每个原始雷达图像进行处理获得每个原始雷达图像对应的逆风方向全距离积分面积；

再根据采集每个原始雷达图像时的风速，建立风速与逆风方向全距离积分面积的函数关系；

步骤2、在线读取风速待测的雷达图像，计算出风速待测的雷达图像的逆风方向全距离积分面积；

步骤3、将步骤2中计算出的逆风方向全距离积分面积代入步骤1中建立的风速与逆风方向全距离积分面积的函数，得到采集步骤2读取的雷达图像时的海表面风速。

本发明的有益效果是：

本发明利用航海雷达图像回波强度的距离衰减进行风速反演，克服了距离衰减对风速反演的影响；利用通过直方图统计的方法克服了海表面上固定目标以及异常值的干扰，可以应用在有固定物影响的海域；通过风向与方位角的依赖关系计算出逆风方向的衰减函数值，然后通过衰减函数值计算计算逆风方向全距离积分面积，以反演风速，克服了雷达图像视场受到遮挡对风速反演的影响，可以从受遮挡的雷达图像中反演海表面风速。

通过实验表明，采用本发明方法可以显著提高风速反演的精度，风速反演的均方根误差仅为2.4m/s。

附图说明

图1为本发明的一种基于航海雷达图像的海面风速反演方法的流程图；

图2为极坐标系下的整幅雷达原始图像；

图3为极坐标系下的归一化后的整幅雷达原始图像；

图4为距离为7.5m处的回波强度直方统计图；

图5为雷达原始图像的一维理想衰减数据；

图6为根据理想衰减数据拟合出的理想衰减模型；

图7为根据衰减水平分量拟合出的余弦函数；

图8为逆风方向全距离积分面积曲线图；

图9为风速与逆风方向全距离积分面积的一次函数曲线图；

图10为两种方法测得的风速与参考风速的比对图。

具体实施方式

下面通过具体实施方式结合附图对本申请作进一步详细说明。显然，所描述的实施方式仅仅是本发明的一部分实施方式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通技术人员在没有做出创造性劳动前提下所获得的其他实施方式，都属于本发明保护的范围。

具体实施方式一、结合图1说明本实施方式。本实施方式所述的一种基于航海雷达图像的海面风速反演方法，所述方法具体包括以下步骤：

步骤1、离线读取雷达文件得到原始雷达图像集，分别对每个原始雷达图像进行处理获得每个原始雷达图像对应的逆风方向全距离积分面积；

再根据采集每个原始雷达图像时的风速，建立风速与逆风方向全距离积分面积的函数关系；

步骤2、在线读取风速待测的雷达图像，采用步骤1的方法计算出风速待测的雷达图像的逆风方向全距离积分面积；

步骤3、将步骤2中计算出的逆风方向全距离积分面积代入步骤1中建立的风速与逆风方向全距离积分面积的函数，得到采集步骤2读取的雷达图像时的海表面风速。

本发明利用雷达图像中回波强度与海面径向距离的依赖关系，推导出雷达图像的理想衰减模型。排除遮挡区域后，通过理想衰减模型得出每个方位角度的衰减函数并求出衰减水平分量，即相对平均回波强度。利用灰度水平的方位角依赖性来反演风向，通过逆风方向的衰减水平分量求出逆风方向的衰减模型，并建立衰减模型积分面积与参考风速的三次多项式关系来反演风速。

具体实施方式二：本实施方式与具体实施方式一不同的是，所述步骤1的具体过程为：

步骤1.1、离线读取雷达文件获得原始雷达图像，并记录原始雷达图像的方位、径向距离以及回波强度信息，对记录的回波强度信息进行同频干扰抑制处理后，再对处理后的回波强度信息进行归一化；

步骤1.2、基于归一化后的回波强度信息，计算逆风方向全距离积分面积；

步骤1.3、重复步骤1.1至步骤1.2的过程，获得H个原始雷达图像对应的逆风方向全距离积分面积；

应用最小二乘拟合建立采集原始雷达图像时的风速与逆风方向全距离积分面积的函数关系。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述对记录的回波强度信息进行同频干扰抑制处理采用的是滤波算法。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述步骤1.2的具体过程为：

步骤1.2.1、计算直方图统计中的阈值T：

T＝k×M

其中：M是原始雷达图像中方位角线数，即在原始雷达图像的极坐标系中，同一径向距离上像素点的个数，k是经验参数；

步骤1.2.2、在原始雷达图像的极坐标系中，选取出到原始雷达图像中心(即雷达天线在图像中的位置)径向距离相同的一组像素点，按照选取出像素点的归一化回波强度，将选取出的像素点存储到回波强度范围为0到1的m个直方图窗格中；

若某个直方图窗格中的像素点个数小于阈值T，则该直方图窗格中像素点的归一化回波强度被认为是异常值并排除掉；若某个直方图窗格中的像素点个数不小于阈值T，则该直方图窗格中像素点的归一化回波强度被认为是正常值并保留；

从保留下来的归一化回波强度值中选择出最大的回波强度值，将选择出的回波强度值作为当前径向距离(即当前情况下，选取出的像素点到原始雷达图像中心的径向距离)下的回波强度值；

步骤1.2.3、重复执行步骤1.2.2，直至各个像素点到原始雷达图像中心的径向距离值均被遍历到，分别得到各个径向距离下的回波强度值；

步骤1.2.4、应用最小二乘拟合将各个径向距离下的回波强度值拟合至理想距离衰减模型，确定出理想距离衰减模型的回归参数；

步骤1.2.5、计算逆风方向的衰减水平分量C

步骤1.2.6、基于逆风方向的衰减水平分量C

再基于逆风方向的衰减函数值计算逆风方向全距离积分面积S

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述从保留下来的回波强度值中选择出最大的回波强度值，其具体为：

其中，r是雷达天线到海表面的径向距离，即选取出的像素点到原始雷达图像中心的径向距离，{X(r，θ)}是在当前径向距离r下被保留下来的归一化回波强度值集合，θ是原始雷达图像的方位角，X

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是，所述步骤1.2.4的具体过程为：

理想距离衰减模型D(r)的公式如下：

其中：b

其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是，所述步骤1.2.5的具体过程为：

步骤1)、初始化权重以及阈值δ：

初始化权重为：

其中，r

将阈值δ的值初始化为1；

步骤2)、初始化方位角θ为0°；

步骤3)、在原始雷达图像中，选取出方位角θ上的全部离散点(每一个像素点即对应一个离散点)的归一化回波强度，并将归一化回波强度为0的离散点的权重置零；

通过非线性最小二乘方法最小化以下公式，计算出当前方位角θ下的衰减水平分量C

其中，X(r

步骤4)、基于步骤3)中计算出的衰减水平分量C

再基于更新后的衰减水平分量C

步骤5)、将顺时针方向定义为方位角增加方向(即从0°开始，顺时针方向上方位角的值是逐渐增加的)，令方位角θ＝Δθ+θ，Δθ为方位角增量，并返回步骤3)；

直至方位角θ达到2π时结束，获得每一个方位角下的最终衰减水平分量C

步骤6)、应用最小二乘法将最终衰减水平分量C

C

其中，a

步骤7)、根据回归参数a

C

其它步骤及参数与具体实施方式一至六之一相同。

具体实施方式八：本实施方式与具体实施方式一至七之一不同的是，所述步骤1.2.6的具体过程为：

其中，D(r)

其中，S

其它步骤及参数与具体实施方式一至七之一相同。

具体实施方式九：本实施方式与具体实施方式一至八之一不同的是，所述建立采集原始雷达图像时的风速与逆风方向全距离积分面积的函数关系，其具体过程为：

V＝β

其中，V是风速，β

其它步骤及参数与具体实施方式一至八之一相同。

具体实施方式十：本实施方式与具体实施方式一至九之一不同的是，所述基于步骤3)中计算出的衰减水平分量C

其中，ω′(r

基于更新后的衰减水平分量C

阈值δ的更新方法为：阈值δ每一次更新后的值为前一次执行步骤3)时采用的阈值δ的一半，即初始化阈值δ＝1，第一次更新后阈值δ的值为0.5，最终阈值δ的值为0.25。

其它步骤及参数与具体实施方式一至九之一相同。

实施例

下面将结合附图对本发明提出的基于雷达图像距离衰减进行航海雷达图像风速反演的方法做进一步的详细说明。

实施流程如图1所示，具体可以分为以下几步：第一步、离线建立风速和逆风方向全距离积分面积的一次函数关系；第二步、在线读取雷达文件得到原始雷达图像，去除同频干扰并归一化处理，计算逆风方向全距离积分面积；第三步、计算海表面风速。

本发明实施例所使用的是X波段航海雷达，该雷达安装在海监船主桅杆的前侧，在短脉冲下进行数据采集，监测径向距离范围为4.2km之内，角度分辨精度约为1度，每幅图像的采集时间约为2.7s，规定以32幅图像作为一个时间序列进行存储，单幅雷达图像的总线数为3000条左右，每条线上有560个点，距离分辨率为7.5m。本发明所用的风参考数据来源于海试过程中的海上船舶风速计数据，与航海雷达相同，该风速计也安装在船主桅杆上，风参数以分钟为单位进行记录，并用以验证雷达反演风速的精度。

结合图1至图10，本发明具体实施步骤为：

第一步，离线建立风速(风速参考数据来源于海试过程中的海上船舶风速计数据)和逆风方向全距离积分面积的一次函数关系。

步骤1.1、离线读取雷达文件得到原始雷达图像，去除同频干扰并归一化处理。其中步骤1.1包括以下步骤：

步骤1.1.1、利用雷达图像处理程序对雷达文件进行加载，记录雷达图像的采集时间、方位、径向距离、回波强度等信息。

步骤1.1.2、利用选定的滤波算法对雷达原始图像进行同频干扰抑制处理。本实施例选定中值滤波的方法对其进行同频干扰抑制，具体方法为将每一像素点的回波强度值以该点3×3邻域窗口内其余8个像元点回波强度的中值代替。

步骤1.1.3、将雷达图像数据归一化处理。本实施例雷达图像数据是将雷达的后向散射信息数字化并存储为14位的灰度深度图像序列，即数字化的后向散射强度范围为0至8192，因此将每个像素点的强度除以8192归一化至0～1的强度范围内。图2为未归一化的雷达原始图像，图3为去除同频干扰并归一化后的雷达原始图像。

步骤1.2、计算逆风方向全距离积分面积。在步骤1.1得到的雷达图像的极坐标系中，通过直方图统计获得理想的衰减数据，并应用最小二乘拟合将衰减数据拟合至理想距离衰减模型，从而确定回归参数；然后根据理想距离衰减模型与每个方位角的雷达图像像素点强度进行拟合比较，确定每个方位角的衰减水平分量。利用水平极化下航海雷达图像风向与方位角依赖性得到拟合曲线并确定逆风方向的衰减水平分量。最后确定逆风方向的衰减函数并计算逆风方向全距离的积分面积。

其中步骤1.2具体包括以下步骤：

步骤1.2.1、利用公式(1)计算直方图统计中的阈值T。在本发明中均匀取了3000条角度数据中的360个，M为360，经验参数k选取为0.01，因此阈值T＝3.6。

阈值T的计算公式：

T＝k×M    (1)

式中：

M代表雷达图像中方位角线数，即同一距离上像素点的个数；

k代表经验参数；

步骤1.2.2、从雷达图像中心开始选取具有相同距离的一维数据，将他们存储至强度范围为0～1的m个直方图窗格中，本实施例将m选取为256，其中数量小于阈值T的窗格被认为异常值排除掉。如图4所示，为距离为7.5m处的直方图强度值统计图，统计数量小于阈值3.6的窗格所在的像素点被认为为异常值并排除。

步骤1.2.3、对于步骤1.2.2中确定的具有相同距离r的像素点，选取最高强度值的像素点作为此距离下的值，在图4中我们选择了强度为0.1521作为7.5m处强度值。之后，以图像中心为原点，选取更大的距离重复步骤1.2.2至步骤1.2.3直到最远的距离，从而确定理想距离衰减数据，如图5所示为最终提取的理想距离衰减数据。

距离r下的最高强度值X

X

式中：

r代表雷达天线到海表面的径向距离，选取范围为0-4.2km；

θ代表雷达图像的方位角；

X代表归一化后的雷达图像强度；

当r遍历到全部距离值后，得到理想距离衰减数据。

步骤1.2.4、选取已有的理想距离衰减函数模型D(r)，应用最小二乘拟合将理想距离衰减数据拟合至理想距离衰减模型，从而确定D(r)的回归参数b

理想距离衰减模型D(r)的公式如下：

式中：

r代表雷达天线到海表面的径向距离；

b

b

步骤1.2.5、初始化加权函数ω以及阈值δ＝1。

ω(r

式中：

r

n＝1，2，3，...，p，p是每个方位角对应的离散点个数，本实施例中为560；

步骤1.2.6、在雷达图像中从方位角0°顺时针开始至2π结束，依次选取每一个方位角下的一维数据X(r，θ)，将X(r，θ)中回波强度为0的权值ω(r

C

其中，n代表径向距离的离散点位置；

步骤1.2.7、通过步骤1.2.6初步估计C

加权函数ω的更新公式如下：

步骤1.2.8、应用最小二乘法将衰减水平分量C

余弦函数如下：

C

θ代表雷达图像的方位角；

a

a

w代表拟合函数峰值点的角度；

步骤1.2.9、根据步骤1.2.8中的公式(7)得到的拟合函数计算逆风方向的衰减水平分量C

计算公式如下：

C

式中：

C

a

a

步骤1.2.10、将逆风方向的衰减水平分量C

计算公式如下：

式中：

D(r)

C

b

b

步骤1.2.11、计算逆风方向的全距离积分面积S

计算公式如下：

D(r)

S

r

r代表雷达天线到海表面的径向距离；

步骤1.3、建立风速和逆风方向的全距离积分面积的一次函数关系。

步骤1.3.1、重复步骤1.1及步骤1.2，获取H个雷达图像数据的积分面积，本发明H＝2000。

步骤1.3.2、如图9所示，应用最小二乘拟合建立参考风速V作为积分面积S

计算公式如下：

V＝β

V代表风速；

S

β

β

第二步、在线读取雷达文件得到原始雷达图像，去除同频干扰并归一化处理，与步骤1.2离线读取雷达文件得到原始雷达图像的方法相同。计算逆风方向全距离积分面积，与步骤1.2离线计算逆风方向全距离积分面积的方法相同。

第三步、计算海表面风速。包括以下步骤：

将所求雷达图像的S

本实施例得到的海表面风速为：V＝76.35×0.284-10.23＝11.45m/s，与参考风速12.1m/s相差0.65m/s。

本发明的算法性能验证是基于海试的测量数据实施的。实验测试时，雷达的监测范围为4.2km之内，每幅图像的采集时间约为2.7s，规定以32幅图像作为一个时间序列进行存储，雷达图像的总线数为3000条左右，每条线上有560个点，距离分辨率为7.5m。本发明所用的风参考数据来源于海试的风速计数据，与航海雷达相同，该风速计也安装在船主桅杆上，风参数以分钟为单位进行记录，并用以验证雷达反演风场的精度。本发明随机选用了2012年3月7日到2012年3月8日的2000组数据来建立风速与雷达图像的关系，并选用了2012年3月7日到2012年3月8日的15914组数据来验证本发明的精度。

为了更好的进行分析，本发明的方法与直接叠加平均求取每个方位角回波强度的单曲线拟合法进行了比较，如图10所示，总体结果如表1所示：

表1风速总体计算结果

由表1可知，本发明基于距离衰减法的风速反演方法效果较好，远好于直接平均的单曲线拟合法，能很好的反演出雷达图像的风速。

本发明所提出的基于雷达图像距离衰减水平分量进行航海雷达图像风速反演的方法提高了风速反演的精度，不仅克服了海况对雷达回波的影响因素，还充分利用了雷达图像的距离衰减特性，使得风速反演精度明显提高。

本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。