CN116540195A - 一种航海雷达图像风向检索方法 - Google Patents

一种航海雷达图像风向检索方法 Download PDF

Info

Publication number
CN116540195A
CN116540195A CN202310352687.4A CN202310352687A CN116540195A CN 116540195 A CN116540195 A CN 116540195A CN 202310352687 A CN202310352687 A CN 202310352687A CN 116540195 A CN116540195 A CN 116540195A
Authority
CN
China
Prior art keywords
radar image
wind direction
radar
intensity
azimuth
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202310352687.4A
Other languages
English (en)
Inventor
卢志忠
于环宇
文保天
张玉莹
孙雷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Engineering University
Original Assignee
Harbin Engineering University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN202310352687.4A priority Critical patent/CN116540195A/zh
Publication of CN116540195A publication Critical patent/CN116540195A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01P—MEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
    • G01P13/00—Indicating or recording presence, absence, or direction, of movement
    • G01P13/02—Indicating direction only, e.g. by weather vane
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/95—Radar or analogous systems specially adapted for specific applications for meteorological use
    • G01S13/956—Radar or analogous systems specially adapted for specific applications for meteorological use mounted on ship or other platform
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/28—Details of pulse systems
    • G01S7/285—Receivers
    • G01S7/292—Extracting wanted echo-signals
    • G01S7/2923—Extracting wanted echo-signals based on data belonging to a number of consecutive radar periods
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/28—Details of pulse systems
    • G01S7/285—Receivers
    • G01S7/295—Means for transforming co-ordinates or for evaluating data, e.g. using computers
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Electromagnetism (AREA)
  • Theoretical Computer Science (AREA)
  • Ocean & Marine Engineering (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种航海雷达图像风向检索方法，对原始雷达图像去除同频干扰，并进行归一化；在图像的极坐标系中，按照设定距离步长选取距离雷达图像中心距离相等的像素点作为一组，在每组像素点中剔除强度异常的像素点，然后在每组像素点中选取强度最大的像素点，并根据强度最大的像素点确定一维理想距离衰减数据然后应用最小二乘拟合将一维理想距离衰减数据拟合至理想距离衰减模型，并确定其回归参数；理想距离衰减模型与每个方位角的雷达图像像素点强度进行拟合比较，确定每个方位衰减水平分量；根据水平极化下航海雷达图像风向与方位角依赖性检索海面风向。本发明改善了风向检索精度，有效地减少雷达图像中异常值干扰以及低风速造成的误差。

Description

一种航海雷达图像风向检索方法
技术领域
本发明属于海洋遥感技术领域，涉及一种航海雷达图像风向检索方法，特别是一种基于回波强度距离衰减的航海雷达图像风向检索方法。
背景技术
海面风场是研究海气之间互相作用的重要参数，海面风场通过调节热量、水汽、海气通量和颗粒物,调节大气和海洋之间的耦合作用,从而维持全球和区域气候。目前海面风向主要通过位于船的桅杆上的现场传感器(如风速计)获取，但风速计容易受大气湍流以及船体上层结构引起的气流紊乱的影响。即使是风速计被安装在无遮挡的位置，风参数测量的误差也可能很高。近几十年来，X波段航海雷达以其高分辨率和及时反馈的优势，成为获取海表面风场信息的有效方法之一。并且，由于大部分船只都已经装备了X波段航海雷达，相对于其他遥感仪器，具有较小的额外成本，这些研究为航海雷达提取海面风场信息奠定了基础。
目前，针对航海雷达图像风向检索研究领域的主要成果如下。有两类方法可用于从X波段航海雷达图像序列中检索风向。第一类方法利用归一化雷达截面(NRCS)和雷达照射迎风方向之间的关系，但是应用此方法获取风向需要整个360°方位无遮挡的雷达图像。2012年，Lund在利用雷达后向散射强度对逆风方向的依赖，提出了一种使用最小二乘法将后向散射强度作为方位角的函数拟合到一个谐波函数上确定风向的方法(Lund B,GraberH C,Romeiser R.Wind Retrieval From Shipborne Nautical X-Band Radar Data[J].IEEE Transactions 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 vector retrieval from marineradar image sequences acquired at grazing angle[J].Journal of Atmospheric andOceanic Technology,2013,30(1):127-139.)。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.)。同年，Liu对Lund中的方法进行了改进，提出了双曲线拟合的方法，并增加了一条曲线拟合(Liu,Ying,Vicen-Bueno,et al.Comparison of Algorithms for Wind Parameters Extraction FromShipborne X-Band Marine Radar Images[J].IEEE Journal of Selected Topics inApplied Earth Observations&Remote Sensing,2015.)。但当风向受阻时，附加曲线拟合效果不佳。海洋是一种随机扰动的过程，并具有多种模型，上述方法在风浪的情况下能很好的工作，但是在涌浪条件下，可能无法获得良好的效果，由破碎的海浪引起的泡沫可能对电磁波的后向散射产生严重影响，因此可能影响风的测量。为了解决风向受遮挡以及海面风场模型不稳定的缺点，2022年，Yu等人提出了将离散小波变换和方位角数据扩展的方法，从单个受遮挡的雷达图像中提取风向(Yu H,Wang H,Lu Z.Wind-Direction Estimationfrom Single X-Band Marine Radar Image Improvement by Utilizing the DWT andAzimuth-Scale Expansion Method[J].Entropy,2022,24(6):747.)。第二类方法是基于雷达图像中风条纹，它属于低频信号，且其空间尺度在200～500m左右，并且与风向平行，通过低通滤波滤掉高频海浪信号后，可以得到风条纹，风条纹的方向即为风向。基于此理论，2021年王慧提出了能量谱理论(ESM)来检索雷达图像序列中的海面风向(Wang H,Qiu H,LuZ,et al.An Energy Spectrum Algorithm for Wind Direction Retrieval From X-BandMarine Radar Image Sequences[J].IEEE Journal of Selected Topics in AppliedEarth Observations and Remote Sensing,2021,14:4074-4088.)。
以上方法在不断发展的过程中，不仅考虑了雷达视场受到遮挡情况，也考虑了风向方位角受到阻挡时的情况，并获得了较好的结果。但是在同一个方位角中，不仅有遮挡的影响，还有固定目标遮挡后形成的目标尾迹，即使在排除目标的像素点也难以计算出此方位角的真实平均强度，即出现局部低值或者高值。因此在曲线拟合过程中，虽然排除了遮挡方向的数据，但可能保留了由于海况比较复杂而产生的多个方位角的数据，导致结果精度降低。此外，雷达回波强度会受到距离的调制，NRCS与海洋表面的距离依赖可以呈现在雷达图像中，雷达回波强度会随着距离的增大而衰减。较远的海面回波信号不仅会随着距离的增大不断减小，还会出现较多的零强度值像素点，这会造成这个方位角下求得平均回波强度较小，进而导致风向的检索精度降低。有研究表明，雷达图像中发生的任何径向距离依赖都是海洋表面NRCS的掠射角依赖的结果，NRCS与径向距离为三次方衰减关系。因此，当从航海雷达图像中检索风的信息时，不能忽略距离衰减的影响。
发明内容
针对上述现有技术，本发明要解决的技术问题是提供一种航海雷达图像风向检索方法，利用雷达图像中回波强度距离衰减特性从X波段航海雷达所获取的图像中检索海表面风向。
为解决上述技术问题，本发明的一种航海雷达图像风向检索方法，包括：
步骤1、对读取雷达文件得到的原始雷达图像去除同频干扰，并对雷达图像像素点强度进行归一化处理；
步骤2、在步骤1得到的雷达图像的极坐标系中，按照设定距离步长选取距离雷达图像中心距离相等的像素点作为一组，在每组像素点中剔除强度异常的像素点，然后在每组像素点中选取强度最大的像素点，并根据强度最大的像素点确定理想距离衰减数据Xra(r)，Xra(r)＝maxθ{X(r,θ)},θ∈[0,2π)，其中，r表示雷达天线到海表面的径向距离，θ表示雷达图像的方位角，X表示归一化后的雷达图像强度；然后应用最小二乘拟合将Xra(r)拟合至理想距离衰减模型D(r)，从而确定D(r)的回归参数；
步骤3、根据理想距离衰减模型D(r)与每个方位角的雷达图像像素点强度进行拟合比较，确定每个方位角的衰减水平分量；
步骤4、根据水平极化下航海雷达图像风向与方位角依赖性检索海面风向。
进一步的，步骤2所述在每组像素点中剔除强度异常的像素点包括：
将0-1的强度范围均匀划分为m个直方图窗格，然后将像素点分别存储至对应强度的直方图窗格中，像素点数量小于阈值T的窗格中的像素点被认为异常并剔除掉。
进一步的，所述阈值T满足：
T＝k×N
式中，N表示雷达图像中方位角线数，即同一距离上像素点的个数；k为0-1之间的常值。
进一步的，所述理想距离衰减模型D(r)为：
式中：r为雷达天线到海表面的径向距离，b0和b1为回归参数。
进一步的，步骤3所述根据理想距离衰减模型D(r)与每个方位角的雷达图像像素点强度进行拟合比较，确定每个方位角的衰减水平分量包括：
步骤3.1、初始化迭代次数N＝1，从方位角0°顺时针开始至2π，依次选取雷达图像中每一个方位角下的雷达图像像素点强度X(r,θ)，将X(r,θ)中回波强度为0的权值ω(r)置零，权值ω(r)为：
其中，r为雷达天线到海表面的径向距离，l为雷达的距离分辨率；n为径向距离的离散点位置，p为径向距离的最大离散点位置；
然后，通过非线性最小二乘法计算出每一个方位角的衰减水平分量Cθ，Cθ的计算公式为：
其中，δ为设定阈值，δ≤1。
步骤3.2、判断是否达到设定的迭代次数N，若达到，则输出衰减水平分量Cθ；否则，令N＝N+1，利用迭代约束对得到的Cθ的值进行优化，更新加权函数ω(r)以及阈值δ，具体为：阈值δ更新为δ/2；返回步骤3.1。
进一步的，步骤4所述根据水平极化下航海雷达图像风向与方位角依赖性检索海面风向包括：
步骤4.1、采用最小二乘法将衰减水平分量Cθ作为方位角的函数拟合为余弦函数：
Cθ＝a0+a1cos2(0.5(θ-w))
其中，θ为雷达图像的方位角，a0、a1为回归参数，w表示拟合函数峰值点的角度；
步骤4.2、将峰值点角度w转换为相对风向：
其中，θre表示相对风向；
步骤4.3、将相对风向转换为绝对风向：
θwind＝Rem(θre+θship,360)
其中，θwind表示绝对风向，θship表示船艏方向。
本发明的有益效果：本发明提出了一种基于回波强度距离衰减特性的航海雷达图像风向反演方法。利用雷达图像中回波强度与海面径向距离的依赖关系，推导出雷达图像的理想衰减模型。在低掠入射下工作的标准航海水平极化X波段雷达采集的图像在迎风方向只有一个回波强度峰值。通过理想衰减模型得出每个方位角度下的衰减函数并求出衰减水平分量，即平均回波强度，通过水平分量求出风向。排除遮挡区域后，通过灰度水平的方位角依赖性来检索风向。利用实测数据进行了实验验证，并通与将回波强度直接平均的单曲线拟合法比较，验证了该算法的有效性和优越性。与现有技术相比，利用本发明所提出的风向检索方法，其优点在于：
本发明利用航海雷达图像回波强度的距离衰减进行风向检索，克服了海表面上固定障碍物以及异常值的干扰，可以应用在有固定物影响的海域。
本发明在通过迭代优化以及最小二乘方法计算出每个方位角的衰减水平分量，避免了直接计算回波强度会受到异常值以及风速较小时过多的零点这一影响，提高了低风速时风向检索的精度。
附图说明
图1是极坐标系下的整幅雷达原始图像；
图2是极坐标系下的归一化后的整幅雷达原始图像；
图3是距离为7.5m处的回波强度直方统计图；
图4是雷达原始图像的一维理想衰减数据；
图5是根据理想衰减数据拟合出的理想衰减模型；
图6是每个方位角的距离衰减水平分量；
图7是根据衰减水平分量拟合出的余弦函数；
图8是本发明实施方式流程图。
具体实施方式
下面结合说明书附图和实施例对本发明做进一步说明。
结合图8，本发明分为以下几步：第一步为待检测雷达图像的读取和去除同频干扰以及归一化处理，第二步为确定航海雷达图像理想距离衰减模型，第三步为计算每个方位角的衰减水平分量，第四步为计算海表面风向。本发明包括以下步骤：
步骤1，读取雷达文件得到原始雷达图像，去除同频干扰并归一化处理，具体包括：
步骤1.1，利用雷达图像处理程序对雷达文件进行加载，记录雷达图像的采集时间、方位、径向距离、回波强度等信息。
步骤1.2，利用选定的滤波算法对雷达原始图像进行同频干扰抑制处理。
步骤1.3，将雷达图像数据归一化处理。
步骤2，确定航海雷达图像理想距离衰减模型。对于在雷达图像的极坐标中具有相同距离的像素，选取最高强度值的像素点推导理想距离衰减模型，具体包括：
步骤2.1，计算直方图统计中的阈值T。
阈值T的计算公式：
T＝k×N (1)
式中：
N-----雷达图像中方位角线数，即同一距离上像素点的个数
k-----经验参数
步骤2.2，从雷达图像中心开始选取具有相同距离的一维数据，将他们存储范围为0到1的m个直方图窗格中，其中数量小于阈值T的窗格被认为异常值并排除掉；
步骤2.3，对于步骤2.2中确定的具有相同距离r的像素点，选取最高强度值的像素点作为此距离下的值。之后，以图像中心为原点，选取更大的距离重复步骤2.2至步骤2.3直到最远的距离，从而确定一维理想距离衰减数据Xra(r)。
Xra(r)的计算公式如下：
Xra(r)＝maxθ{X(r,θ)}, θ∈[0,2π) (2)
式中：
r-----雷达天线到海表面的径向距离
θ-----雷达图像的方位角
X-----归一化后的雷达图像强度
步骤2.4，选取已有的理想距离衰减函数模型D(r)，应用最小二乘拟合将Xra(r)拟合至理想距离衰减模型，从而确定D(r)的回归参数。
理想距离衰减模型D(r)的公式如下：
式中：
r-----雷达天线到海表面的径向距离
b0-----回归参数
b1-----回归参数。
步骤3，根据理想距离衰减模型确定每个方位角的衰减水平分量，即相对平均回波强度：根据得到的理想衰减模型与每个方位角的一维数据进行拟合比较，确定每个方位角的衰减水平分量，具体包括：
步骤3.1，初始化加权函数ω以及阈值δ＝1。
ω(r)的初始化计算公式如下：
r-----雷达天线到海表面的径向距离
l-----雷达的距离分辨率
n-----径向距离的离散点位置
p-----径向距离的最大离散点位置
步骤3.2，从雷达图像中从方位角0°顺时针开始至2π结束，依次选取每一个方位角下的一维数据X(r,θ)，将X(r,θ)中回波强度为0的权值ω(r)置零，并通过非线性最小二乘方法最小化以下公式，计算出每一个方位角的衰减水平分量Cθ。
Cθ的计算公式如下：
l-----雷达的距离分辨率
n-----径向距离的离散点位置
p-----径向距离的最大离散点位置
步骤3.3，通过步骤3.2初步估计Cθ，然后使用迭代约束对Cθ的值进行优化，在迭代过程中更新加权函数ω以及阈值δ。通过公式(6)更新权重ω，阈值δ在每次迭代结束后更新为原来的一半。更新后的参数ω和δ被传递公式(5)中进行下一次迭代。整个过程需要迭代3次获得最后的衰减水平分量Cθ，每个方位角的衰减水平分量Cθ是与雷达照射方向存在余弦函数的依赖关系，并在迎风方向存在一个峰值。
加权函数ω的更新公式如下：
步骤4，检索海表面风向。利用水平极化下航海雷达图像风向与方位角依赖性检索海面风向，具体包括：
步骤4.1，应用最小二乘法将衰减水平分量Cθ作为方位角的函数拟合到一个余弦函数。
余弦函数如下：
Cθ＝a0+a1cos2(0.5(θ-w)) (7)
θ-----雷达图像的方位角
a0-----回归参数
a1-----回归参数
w-----拟合函数峰值点的角度
步骤4.2，步骤4.1公式(7)中w的为拟合函数峰值点的角度，将峰值点角度转换为相对风向。
相对风向计算公式如下：
θre-----相对风向
w-----拟合函数峰值点的角度
步骤4.3，将相对风向转换为绝对风向。
绝对风向计算公式如下：
θwind＝Rem(θre+θship,360) (9)
θwind-----绝对风向
θre-----相对风向
θship-----船艏方向。
下面结合具体参数给出实施例：
本发明实施用例所使用的是X波段航海雷达，该雷达安装在海监船主桅杆的前侧，在短脉冲下进行数据采集，监测径向距离范围为4.2km之内，角度分辨精度约为1度，每幅图像的采集时间约为2.7s，规定以32幅图像作为一个时间序列进行存储，单幅雷达图像的总线数为2048条，每根线上有560个点，距离分辨率为7.5m，方向分辨率约为0.18度。本发明所用的风参考数据来源于国家海洋局的海上船舶风速计数据，与航海雷达相同，该风速计也安装在船主桅杆上，风参数以分钟为单位进行记录，并用以验证雷达反演风向的精度。
结合附图1～8，本发明具体实施步骤为：
第一步，读取雷达文件得到原始雷达图像，去除同频干扰并归一化。
步骤1.1，利用雷达图像处理程序对雷达文件进行加载，记录雷达图像的采集时间、方位、径向距离、回波强度等信息。
步骤1.2，利用选定的滤波算法对雷达原始图像进行同频干扰抑制处理。本实施例选定中值滤波的方法对其进行同频干扰抑制，具体方法为将每一像素点的回波强度值以该点3×3邻域窗口内其余8个像元点回波强度的中值代替。
步骤1.3，将雷达图像数据归一化处理。本实施例雷达图像数据是将雷达的后向散射信息数字化并存储为14位的灰度深度图像序列，即数字化的后向散射强度范围为0至8192，因此将每个像素点的强度除以8192归一化至0～1的强度范围内，如图所示，图1为未归一化的雷达原始图像，图2为去除同频干扰并归一化后的雷达原始图像。
第二步，确定航海雷达图像理想距离衰减模型。理想衰减模型的确定包括以下步骤：
步骤2.1，利用公式(1)计算直方图统计中的阈值T。在本发明中N为2048，经验参数k选取为0.01，因此阈值T＝20.48。
阈值T的计算公式：
T＝k×N (1)
式中：
N-----雷达图像中方位角线数，即同一距离上像素点的个数
k-----经验参数
步骤2.2，从雷达图像中心开始选取具有相同距离的一维数据，将他们存储至强度范围为0～1的m个直方图窗格中，本实例将m选取为256，其中数量小于阈值T的窗格被认为异常值排除掉。如图3所示为距离为7.5m处的直方图强度值统计图，统计数量小于阈值20.48的窗格所在的像素点被认为为异常值并排除。
步骤2.3，对于步骤2.2中确定的具有相同距离r的像素点，选取最高强度值的像素点作为此距离下的值，在图3中我们选择了强度为0.5186作为7.5m处强度值。之后，以图像中心为原点，选取更大的距离重复步骤2.2至步骤2.3直到最远的距离，从而确定理想距离衰减数据Xra(r)，如图4所示为最终提取的Xra(r)。
Xra(r)的计算公式如下：
Xra(r)＝maxθ{X(r,θ)}, θ∈[0,2π) (2)
式中：
r-----雷达天线到海表面的径向距离
θ-----雷达图像的方位角
X-----归一化后的雷达图像强度
步骤2.4，根据相关文献选取已有的理想距离衰减函数模型D(r)，应用最小二乘拟合将Xra(r)拟合至理想距离衰减模型，从而确定D(r)的回归参数。将公式(3)的理想衰减模型拟合至图4中的Xra(r)，得到回归参数b0＝0.632，b1＝0.875，画出理想衰减模型的曲线图，如图5所示。
理想距离衰减模型D(r)的公式如下：
式中：
r-----雷达天线到海表面的径向距离
b0-----回归参数
b1-----回归参数
第三步，计算每个方位角的衰减水平分量。方位角的衰减水平分量表征此方位角的衰减情况与理想衰减模型的差异，用以衡量此方位角下的相对平均回波强度。包括以下步骤：
步骤3.1，初始化加权函数ω以及阈值δ＝1。
ω(r)的初始化计算公式如下：
r-----雷达天线到海表面的径向距离
l-----雷达的距离分辨率，本实例的距离分辨率为7.5m
n-----径向距离的离散点位置
p-----径向距离的最大离散点位置，本实例为560
步骤3.2，从雷达图像中从方位角0°顺时针开始至2π结束，依次选取每一个方位角下的一维数据X(r,θ)，将X(r,θ)中回波强度为0的权值ω(r)置零，并通过非线性最小二乘方法最小化以下公式，计算出每一个方位角的衰减水平分量Cθ。
Cθ的计算公式如下：
l-----雷达的距离分辨率，本实例的距离分辨率为7.5m
n-----径向距离的离散点位置
p-----径向距离的最大离散点位置，本实例为560
步骤3.3，通过步骤3.2初步估计Cθ，然后使用迭代约束对Cθ的值进行优化，在迭代过程中更新加权函数ω以及阈值δ。通过公式(6)更新权重ω，阈值δ在每次迭代结束后更新为原来的一半。更新后的参数ω和δ被传递公式(5)中进行下一次迭代。整个过程需要迭代3次获得最后的衰减水平分量Cθ，本发明中的方位角数据线一共有2048条，为了计算效率，选择每隔4条线选取一条线来计算，也就是说，最后会得到512个衰减水平分量Cθ值，每个方位角的衰减水平分量Cθ是与雷达照射方向存在余弦函数的依赖关系，并在迎风方向存在一个峰值。
如图6所示，虚线代表的为船舶桅杆的固定遮挡区域。
加权函数ω的更新公式如下：
第四步为计算海表面风向。利用水平极化下海洋雷达图像风向与方位角依赖性检索海面风向，即水平极化下的X波段航海雷达在逆风方向上只存在一个回波强度的峰值。包括以下步骤：
步骤4.1，应用最小二乘法将衰减水平分量Cθ作为方位角的函数拟合到一个余弦函数。首先去除固定遮挡区，本发明采用的数据在140°～200°会产生一个大约60°的遮挡区域，因此将这个区域的数据排除。并将去除遮挡区域的衰减水平分量Cθ拟合至公式(7)，得到参数为a0＝0.446，a1＝0.08，w＝-31.95。
余弦函数如下：
Cθ＝a0+a1cos2(0.5(θ-w)) (7)
θ-----雷达图像的方位角
a0-----回归参数
a1-----回归参数
w-----峰值点的位置
步骤4.2，步骤4.1公式(7)中w的为拟合曲线的峰值点位置，将峰值点位置转换为相对风向。根据公式(8)，本例相对风向θre＝360-31.95＝328.05°
相对风向计算公式如下：
θre-----相对风向
w-----峰值点的位置
步骤4.3，本例参考风向为24°，船首方向为40.6°，根据公式(9)，因此本例绝对风向θwind＝Rem(328.05+40.6，360)＝8.65°。与参考风向24°相差15.35°。
绝对风向计算公式如下：
θwind＝Rem(θre+θship,360) (9)
θwind-----绝对风向
θre-----相对风向
θship-----船艏方向
本发明的算法性能验证是基于国家海洋局海监船的测量数据实施的。实验测试时，雷达的监测范围为4.2km之内，每幅图像的采集时间约为2.7s，规定以32幅图像作为一个时间序列进行存储，雷达图像的总线数为2048条，每根线上有560个点，距离分辨率为7.5m，方向分辨率为0.18度。发明所用的风参考数据来源于国家海洋局的风速计数据，与航海雷达相同，该风速计也安装在船主桅杆上，风参数以分钟为单位进行记录，并用以验证雷达反演风场的精度。
为了更好的进行分析，选取了4个风速不同的时间段分别400幅雷达图像，一共1600幅雷达图像进行测试。并与直接叠加平均求每个方位角回波强度的单曲线拟合法进行了比较，其中第一阶段为仅存在固定遮挡区的中风速雷达数据，第二阶段为仅存在固定遮挡区的高风速雷达数据，第三阶段为仅存在固定遮挡区的低风速雷达数据，第四阶段也为高风速雷达数据，但是不仅存在桅杆造成的固定遮挡区域，而且雷达图像中也存在船舶等目标的干扰。总体结果如表一所示：
表一风向总体检索结果
由上表可知，基于距离衰减法的风向检索方法效果较好，远好于直接平均的单曲线拟合法，能很好的检索出雷达图像的风向，尤其是在低风速和存在固定物干扰的情况下。
本发明所提出的基于雷达图像距离衰减水平分量进行航海雷达图像风向检索的方法提高了风向检索的精度，该方法不仅克服了海况对雷达回波的影响因素，还充分利用了雷达图像的距离衰减特性，使得风向检索精度明显提高。

Claims (6)

1.一种航海雷达图像风向检索方法，其特征在于，包括：
步骤1、对读取雷达文件得到的原始雷达图像去除同频干扰，并对雷达图像像素点强度进行归一化处理；
步骤2、在步骤1得到的雷达图像的极坐标系中，按照设定距离步长选取距离雷达图像中心距离相等的像素点作为一组，在每组像素点中剔除强度异常的像素点，然后在每组像素点中选取强度最大的像素点，并根据强度最大的像素点确定理想距离衰减数据Xra(r)，Xra(r)＝maxθ{X(r,θ)},θ∈[0,2π)，其中，r表示雷达天线到海表面的径向距离，θ表示雷达图像的方位角，X表示归一化后的雷达图像强度；然后应用最小二乘拟合将Xra(r)拟合至理想距离衰减模型D(r)，从而确定D(r)的回归参数；
步骤3、根据理想距离衰减模型D(r)与每个方位角的雷达图像像素点强度进行拟合比较，确定每个方位角的衰减水平分量；
步骤4、根据水平极化下航海雷达图像风向与方位角依赖性检索海面风向。
2.根据权利要求1所述的一种航海雷达图像风向反演方法，其特征在于：步骤2所述在每组像素点中剔除强度异常的像素点包括：
将0-1的强度范围均匀划分为m个直方图窗格，然后将像素点分别存储至对应强度的直方图窗格中，像素点数量小于阈值T的窗格中的像素点被认为异常并剔除掉。
3.根据权利要求2所述的一种航海雷达图像风向反演方法，其特征在于：所述阈值T满足：
T＝k×N
式中，N表示雷达图像中方位角线数，即同一距离上像素点的个数；k为0-1之间的常值。
4.根据权利要求1所述的一种航海雷达图像风向反演方法，其特征在于：所述理想距离衰减模型D(r)为：
式中：r为雷达天线到海表面的径向距离，b0和b1为回归参数。
5.根据权利要求1所述的一种航海雷达图像风向反演方法，其特征在于：步骤3所述根据理想距离衰减模型D(r)与每个方位角的雷达图像像素点强度进行拟合比较，确定每个方位角的衰减水平分量包括：
步骤3.1、初始化迭代次数N＝1，从方位角0°顺时针开始至2π，依次选取雷达图像中每一个方位角下的雷达图像像素点强度X(r,θ)，将X(r,θ)中回波强度为0的权值ω(r)置零，权值ω(r)为：
其中，r为雷达天线到海表面的径向距离，l为雷达的距离分辨率；n为径向距离的离散点位置，p为径向距离的最大离散点位置；
然后，通过非线性最小二乘法计算出每一个方位角的衰减水平分量Cθ，Cθ的计算公式为：
其中，δ为设定阈值，δ≤1。
步骤3.2、判断是否达到设定的迭代次数N，若达到，则输出衰减水平分量Cθ；否则，令N＝N+1，利用迭代约束对得到的Cθ的值进行优化，更新加权函数ω(r)以及阈值δ，具体为：阈值δ更新为δ/2；返回步骤3.1。
6.根据权利要求1所述的一种航海雷达图像风向反演方法，其特征在于：步骤4所述根据水平极化下航海雷达图像风向与方位角依赖性检索海面风向包括：
步骤4.1、采用最小二乘法将衰减水平分量Cθ作为方位角的函数拟合为余弦函数：
Cθ＝a0+a1cos2(0.5(θ-w))
其中，θ为雷达图像的方位角，a0、a1为回归参数，w表示拟合函数峰值点的角度；
步骤4.2、将峰值点角度w转换为相对风向：
其中，θre表示相对风向；
步骤4.3、将相对风向转换为绝对风向：
θwind＝Rem(θre+θship,360)
其中，θwind表示绝对风向，θship表示船艏方向。
CN202310352687.4A 2023-04-04 2023-04-04 一种航海雷达图像风向检索方法 Pending CN116540195A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310352687.4A CN116540195A (zh) 2023-04-04 2023-04-04 一种航海雷达图像风向检索方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310352687.4A CN116540195A (zh) 2023-04-04 2023-04-04 一种航海雷达图像风向检索方法

Publications (1)

Publication Number Publication Date
CN116540195A true CN116540195A (zh) 2023-08-04

Family

ID=87451378

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310352687.4A Pending CN116540195A (zh) 2023-04-04 2023-04-04 一种航海雷达图像风向检索方法

Country Status (1)

Country Link
CN (1) CN116540195A (zh)

Similar Documents

Publication Publication Date Title
Gao A parzen-window-kernel-based CFAR algorithm for ship detection in SAR images
CN105319537B (zh) 基于空间相关性的航海雷达同频干扰抑制方法
CN106569193B (zh) 基于前-后向收益参考粒子滤波的海面小目标检测方法
Liu et al. Comparison of algorithms for wind parameters extraction from shipborne X-band marine radar images
CN109324315B (zh) 基于双层次块稀疏性的空时自适应处理雷达杂波抑制方法
CN113808174B (zh) 基于全卷积网络和卡尔曼滤波的雷达小目标跟踪方法
CN111796288B (zh) 一种基于杂波频谱补偿技术的三坐标雷达动目标处理方法
CN112255607B (zh) 一种海杂波的抑制方法
CN117169882A (zh) 船载雷达海浪信息反演方法
CN115236664A (zh) 一种航海雷达图像反演有效波高的方法
CN113514833B (zh) 一种基于海浪图像的海面任意点波向的反演方法
CN108318881A (zh) 基于k参数的航海雷达图像降雨识别方法
Wang et al. An energy spectrum algorithm for wind direction retrieval from X-band marine radar image sequences
CN108957430A (zh) 基于距离-多普勒图的高频雷达射频干扰区域提取方法
CN117436329A (zh) 一种基于机器学习的合成孔径雷达图像风向反演方法
CN116540195A (zh) 一种航海雷达图像风向检索方法
CN115407282B (zh) 一种短基线下基于干涉相位的sar有源欺骗干扰检测方法
Sun et al. A Wave Texture Difference Method for Rainfall Detection Using X‐Band Marine Radar
Lu et al. Research on rainfall identification based on the echo differential value from X-band navigation radar image
CN113406634B (zh) 一种基于时域相位匹配的空间高速自旋目标isar三维成像方法
CN115205216A (zh) 一种基于显著性和加权引导滤波的红外小目标检测方法
CN116664974A (zh) 一种基于航海雷达图像的海面风速反演方法
Fan et al. An automatic correction method of marine radar rainfall image based on continuous wavelet transform
CN113552563A (zh) 垂测信息和高频地波雷达杂波信息的对应性分析方法
CN112731386B (zh) 基于杂波图匹配的机场跑道异物检测方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination