CN117169882A - 船载雷达海浪信息反演方法 - Google Patents
船载雷达海浪信息反演方法 Download PDFInfo
- Publication number
- CN117169882A CN117169882A CN202311129796.6A CN202311129796A CN117169882A CN 117169882 A CN117169882 A CN 117169882A CN 202311129796 A CN202311129796 A CN 202311129796A CN 117169882 A CN117169882 A CN 117169882A
- Authority
- CN
- China
- Prior art keywords
- radar
- wave
- sea
- image
- inversion
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 59
- 238000012937 correction Methods 0.000 claims abstract description 30
- 238000003384 imaging method Methods 0.000 claims abstract description 19
- 230000003595 spectral effect Effects 0.000 claims abstract description 16
- 238000012545 processing Methods 0.000 claims abstract description 10
- 230000008569 process Effects 0.000 claims description 17
- 238000001228 spectrum Methods 0.000 claims description 17
- 230000009466 transformation Effects 0.000 claims description 10
- 239000011159 matrix material Substances 0.000 claims description 3
- 239000013589 supplement Substances 0.000 claims description 3
- 238000013519 translation Methods 0.000 claims description 3
- 238000004590 computer program Methods 0.000 claims 4
- 230000008859 change Effects 0.000 abstract description 6
- 206010034719 Personality change Diseases 0.000 abstract description 5
- 238000013507 mapping Methods 0.000 abstract description 3
- 230000000694 effects Effects 0.000 description 9
- 238000004458 analytical method Methods 0.000 description 6
- 230000002829 reductive effect Effects 0.000 description 5
- 238000010183 spectrum analysis Methods 0.000 description 4
- 230000003068 static effect Effects 0.000 description 4
- 238000003702 image correction Methods 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 230000004044 response Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 230000001788 irregular Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000001629 suppression Effects 0.000 description 2
- 241001442234 Cosa Species 0.000 description 1
- 244000089409 Erythrina poeppigiana Species 0.000 description 1
- GIYXAJPCNFJEHY-UHFFFAOYSA-N N-methyl-3-phenyl-3-[4-(trifluoromethyl)phenoxy]-1-propanamine hydrochloride (1:1) Chemical compound Cl.C=1C=CC=CC=1C(CCNC)OC1=CC=C(C(F)(F)F)C=C1 GIYXAJPCNFJEHY-UHFFFAOYSA-N 0.000 description 1
- 235000009776 Rathbunia alamosensis Nutrition 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000002592 echocardiography Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000000670 limiting effect Effects 0.000 description 1
- 238000002156 mixing Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000001615 p wave Methods 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000005096 rolling process Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- 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
Abstract
船载雷达海浪信息反演方法，解决了雷达由于姿态变化使非快拍成像导致图像畸变，进而使反演精度低的问题，属于信号处理技术领域。本发明包括：建立非快拍成像校正模型，考虑雷达旋转期间船舶的姿态变化，根据载体艏向方位信息分割图像区域并采用区域滑窗的方式转换径向序数差，构建数据‑图像各个径向映射关系，实现雷达图像动态校正。针对校正后的雷达图像，确定雷达图像中的反演区域，重构海面波高程；对重构的海面波高程进行功率谱密度估计，通过波谱密度和峰值波数实现海浪参数的估计。
Description
技术领域
本发明涉及一种船载雷达海浪信息反演方法，属于涉及X波段导航雷达信号处理技术领域、随机海浪模型、海洋信息反演领域。
背景技术
船载X波段导航雷达从海面散射回波中获取海浪信息，实现了在船舶航行条件下的海洋环境信息探测。与岸基雷达海浪反演相比，船载雷达海浪反演通过导航雷达实时测量海浪，提供船舶航行中的海浪信息，能够实现对海域的灵活观测，适用性更为广泛。
现有的海浪反演算法大多需要利用现场波浪浮标等仪器进行校准，在相关研究(Qiu J,Zhang B,Chen Z,et al.A New Modulation Transfer Function With Range andAzimuth Dependence for Ocean Wave Spectra Retrieval From X-Band Marine RadarObservations[J].IEEE Geoscience and Remote Sensing Letters,2017,14(8):1373-1377.)(Wendy N,Juan V,Alejandro O.Estimation of sea state parameters using X-band marine radar technology in coastal areas,Proc.SPIE,2018vol.10773.)中，对海杂波图像的特性进行分析，应用三维傅里叶变换从图像中提取海浪信息，建立近似线性拟合关系，并通过浮标校准标定参数，实现了比较准确的估计效果。然而，傅里叶变换法需要建立在波浪场的空间均匀性和时间稳定性假设的基础上才能够成立，并且浮标等现场仪器不便于随船航行，反演算法仅能在岸基下针对固定海域进行研究和验证，需要设计在船载条件下实现免校准海浪反演的方法。
在船舶航行过程中，由于载体自身不规则运动以及雷达非快拍成像方式，引起雷达照射区域的重叠、缺失和非均匀成像等现象，导致雷达图像反演的海浪信息与真实值之间存在较大的误差。相关研究(Lund B,Collins C O,Tamura H,et al.Multi-directionalwave spectra from marine X-band radar[J].Ocean Dynamics,2016,66(8):973-988.)(吴凡.船载雷达海浪信息反演算法研究[D].哈尔滨工程大学,2022.)中指出，船舶航行过程中会由于姿态的变化会引起雷达反演区域发生偏移，并在单快拍成像假设下具体分析了反演区域偏移程度，但针对雷达非快拍成像导致的图像畸变问题却没有加以考虑。并且传统谱分析法等反演算法依赖浮标设备的标定信息，需要根据当前海域波浪特征校准反演参数，不适用于随船航行。
发明内容
针对雷达由于姿态变化使非快拍成像导致图像畸变，进而使反演精度低的问题，本发明提供一种有效提高反演精度的船载雷达海浪信息反演方法。
本发明的一种船载雷达海浪信息反演方法，包括：
S1、建立非快拍成像校正模型：
在当前帧雷达扫描成像期间，若船舶艏向按顺时针匀速转过的角度为θ1，确定k1，在当前帧雷达图像[1,M-k1]的雷达径向线范围内，通过均匀插值k1条相邻径向线展开图像中的海浪信息；
若船舶艏向按逆时针匀速转过的角度为确定k2，在当前帧雷达图像[1,M]的雷达径向线范围内，通过均匀抽取k2条径向线，同时取前一帧雷达图像中[M-k2+1,M]范围内的径向线，补充当前帧雷达图像缺失的径向线；
若船舶艏向非匀速转动，将当前帧雷达图像分割成若干扇形子图区域进行滑窗校正，将每个子区域视为匀速转动，根据顺时针或逆时针的情况，在各子区域内采用均匀插值或抽取径向线的方式降低径向累计偏移，通过惯导信息确定需要操作的径向线数；
雷达图像为M个圆心角/>相等、半径相同的扇形所构成的环形图像；
S2、针对校正后的雷达图像，确定雷达图像中的反演区域，重构海面波高程；
S3、对重构的海面波高程进行功率谱密度估计，通过波谱密度和峰值波数实现海浪参数的估计。
作为优选，S2中，通过二维连续小波变换，在反演区域的波数域中去除与海浪不相关的噪声分量，应用连续小波变换的逆过程重构海面波高程。
作为优选，通过二维连续小波变换将反演区域的海浪信号I转换至波数域W进行处理：
其中，r表示以雷达位置为原点的平面位置矢量，r＝(x,y)，x,y分别表示以东、北为正向的坐标，I(r)表示海洋遥感图像中r处的回波强度，表示母小波ψ(r)的缩放因子，/>为小波平移向量，/>是小波旋转角为/>的旋转矩阵，Cψ是与所用小波有关的容许常数，/>分别表示波数域下的二维雷达图像和母小波函数，波数向量k＝(kx,ky)，kx,ky分别表示k在x，y方向上的分量；
在波数域下加入高斯滤波器和Gabor滤波器/>去除与海浪不相关的噪声分量：
其中，是高斯滤波器x、y方向上的方差，/>是Gabor滤波器的方差；/> 分别表示峰值波数在x、y方向上的分量；
应用连续小波变换的逆过程重构海面波高程
其中，Cψ,δ表示以狄拉克函数为母小波实现逆CWT的容许常数，表示在r处波数域下的海浪信号；
根据确定雷达相距R0处的重构海面波高程η(R0,t)：
其中，Φ、/>分别表示反演区域掠射角度、雷达天线水平波束宽度、重构波高程图像的平均值。
作为优选，所述母小波为Morlet小波：
其中，k0表示中心波数，ε对应于各向异性参数，σ表示形状参数。
作为优选，所述母小波函数为函数为狄拉克函数(δ)。
作为优选，S3中，用Welch修正周期图法对重构的海面波高程进行功率谱密度估计，对每一段数据做加窗处理，通过波谱密度和峰值波数实现海浪参数的估计。
作为优选，海浪参数包括海浪波峰波长、波峰浪向、波高和波峰周期。
本发明的有益效果，本发明考虑雷达旋转期间船舶的姿态变化，根据载体艏向方位信息分割图像区域并采用区域滑窗的方式转换径向序数差，构建数据-图像各个径向映射关系，实现雷达图像动态校正。通过对图像校正前后的反演结果进行对比分析，验证了动态校正方案能够有效降低载体晃动所造成的反演误差；通过对不同算法的反演结果进行对比，验证了改进的海浪反演算法具有较高的反演精度。
附图说明
图1是船舶艏向、雷达及海域径向模型示意图。
图2是在艏向顺时针匀速转动情况下分析海域径向变化的示意图。
图3是在艏向逆时针匀速转动情况下分析海域径向变化的示意图。
图4是在艏向非匀速转动情况下分析海域径向变化的示意图：其中φ1,φ2,φ3,φ4分别表示4条惯导信息所提供的当前船舶艏向方位，a1,a2,a3,a4分别表示接收4条惯导信息的每个时间间隔内艏向变化角度θ1′,θ2′,θ3′,θ4′对应的径向线个数，即θ1′＝φ1,θ2′＝φ2-φ1,θ3′＝φ3-φ2,θ4′＝φ4-φ3，β1,β2,β3,β4分别是接收4条惯导信息的每个时间间隔内雷达旋转角度；
图5是选取雷达海杂波图像的反演区域。
图6是反演区域经过二维连续小波变换得到的小波频谱图像。
图7是滤波器频谱响应，其中(a)为高斯高通滤波器的频谱响应：能够去除与海浪不相关的低波数分量，并消除由于距离向平均雷达截面衰减而出现在f＝0Hz附近的峰值频谱强度；(b)为Gabor带通滤波器的频谱响应：能够去除波谱区域以外的分量，增强图像纹理特征，提高重构后海面波高程的波浪纹理识别度。
图8是重构的海面波高程图。
图9是Welch修正图周期法估计PSD。
图10是不同反演算法对海浪有效波高的反演结果对比：“o”表示谱分析法3-D FFT反演结果，“+”表示阴影抑制法反演结果，“△”表示本发明算法反演结果。其中第一列和第二列分别采用未经图像动态校正和校正后的图像反演结果。
图11是反演结果，其中(a)为波向的反演结果，(b)为波周期的反演结果。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
下面结合附图和具体实施例对本发明作进一步说明，但不作为本发明的限定。
本实施方式的船载雷达海浪信息反演方法，包括：
S1、建立非快拍成像校正模型：
首先，分析不规则载体姿态变化对导航雷达天线旋转期间成像结果的影响。雷达提供的首个径向数据方位方向通常与船舶艏向保持一致，当船舶平台完全静止时，下一帧雷达图像的首个径向与前一帧首个径向所代表的方位方向保持一致，即前后帧的首个径向夹角α＝0°。当船舶艏向发生偏移时，夹角α发生变化，单幅雷达图像存在畸变，其所代表的周围海域海浪信息存在重复、缺失等问题，若不加以校正，反演过程中所选的反演区域在地理坐标系下存在偏移，且随着反演区域距离雷达越远，其偏移程度越大，这将造成反演结果出现较大误差。
为了降低雷达图像畸变导致的反演误差，现有技术中考虑了船舶运动并非理想的线性运动、雷达扫描期间海面与船舶并非静止的这些实际影响因素，针对船舶水平运动导致的区域偏移程度进行具体分析，并提出了偏移校正补偿方案，得到了较为准确的反演结果。但对于雷达天线旋转期间的成像问题却没有详细地展开介绍。
本发明先对雷达和惯导数据进行处理，将接收到的导航雷达原始压缩数据解析为极坐标图像，与惯性测量单元提供的惯导信息通过队列形式混流，构建径向序与扫描方位之间的联系，使惯导低频数据绑定到雷达高频径向数据。将雷达极坐标数据记为Dmn＝[Am,Rn]m∈[1,M],n∈[1,N]，Am表示第m个径向，Rn表示第n个距离范围，则雷达数据以径向线矢量A的形式可表示为D＝[A1,A2,…,AM-1,AM]，Αm＝[R1,R2,…,RN-1,RN]。在雷达图像处理过程中，通常需要把极坐标图像转换到直角坐标系中，若单幅雷达图像的数据量是动态的，则导致坐标变换的运算复杂度极高，严重影响反演算法实时性。因此需要保证在雷达极坐标数据固定为M×N的条件下实现图像动态校正。
如图1所示，当船舶平台相对静止时，雷达所有径向线A均匀扫描360°的周围海域，得到的雷达图像可视为M个圆心角相等、半径相同的扇形所构成的环形图像，每个扇形角度为船舶所在海域的海浪信息径向线Xsea＝M条。当船舶平台存在艏摇运动或改变航向时，设正北方向为0°，艏向和首个径向方位均为正北方向，且将顺时针记为雷达旋转正方向，建立的非快拍成像校正模型包括如下三种图像校正情形：
情形1、如图2所示，在当前帧雷达扫描成像期间，若船舶艏向按顺时针匀速转过的角度为θ1，确定k1，则下一帧雷达数据的首个径向角度为θ1，当前帧雷达图像扫描过的角度为(360°+θ1)。此时雷达M条径向数据所包含的海域径向线Xsea＝M+k1条，在[0°,θ1]范围内的海浪信息发生重叠，导致算法选定的反演区域中心[Am,Rn]出现偏移，其偏移程度与径向线序m、转角θ和距离Rn的大小有关，m、θ1和Rn越大，偏移距离L越大：
与相对静止状态下的雷达图像相比，雷达数据有效径向线为M-k1条，图像中的海浪信息存在径向压缩现象。为了保证雷达数据满足M×N并去除重复内容，在当前帧雷达图像[1,M-k1]的雷达径向线范围内，通过均匀插值k1条相邻径向线展开图像中的海浪信息，降低径向旋转偏移，实现雷达图像重复情形下的动态校正。
情形2、如图3所示，若船舶艏向按逆时针匀速转过的角度为确定k2，则下一帧雷达数据的首个径向角度为(-θ2)，当前帧雷达扫描过的角度为(360°-θ2)。此时雷达数据所包含的海域径向线为Xsea＝M-k2条，在[-θ2,0°]范围内的海浪信息发生缺失，导致算法选定的反演区域中心[Am,Rn]出现偏移，其偏移距离L:
该运动状态下的雷达数据有效径向线为M条，但图像中的海浪信息存在径向扩充现象。为了保证雷达数据满足M×N并去除扩充的径向，在当前帧雷达图像[1,M]的雷达径向线范围内，通过均匀抽取k2条径向线，同时取前一帧雷达图像中[M-k2+1,M]范围内的径向线，补充当前帧雷达图像缺失的径向线，实现雷达图像缺失情形下的动态校正。
情形3、如图4所示，若船舶艏向非匀速转动，则需要借助惯导低频数据与雷达高频径向数据之间的联系，将当前帧雷达图像分割成若干扇形子图区域进行滑窗校正，将每个子区域视为匀速转动，根据顺时针或逆时针的情况，在各子区域内采用均匀插值或抽取径向线的方式降低径向累计偏移，通过惯导信息确定需要操作的径向线数；
设雷达第i帧绑定的惯导信息为其中Pi表示惯导信息个数，非固定值，Pi＜＜M，/>Ii与D之间的联系可表示为：
任取第Ip惯导信息，将序数为[mp-1+1,mp]范围内的径向线分割为子区域并展开分析，该子区域下的船舶艏向旋转可视为匀速转动(以顺时针为例)，假设的方位角度为φp-1且与海域径向线序数对应，期间艏向转过的角度为/>则径向/>的方位角度为φp+θ3，雷达扫描过的角度为φp-φp-1+θ3。此时该部分雷达数据所包含的海域径向线为Xsea＝mp-mp-1+k3条，与相对静止状态下的雷达图像相比，海域径向线序数大于雷达径向线序数，图像中的海浪信息存在顺时针偏移现象。然而，雷达图像各子区域均存在旋转偏移，的方位角度为与海域径向线序数并不一定是对应的，因此随着径向线序数的增大会出现偏移累计现象，且由于旋转速度不同导致偏移程度并不相同，影响相邻子区域的图像校正效果。为了保证整体雷达数据满足M×N，结合情形1和情形2的处理思路，在各子区域内采取均匀插值或抽取径向线的方式降低径向累计偏移，通过惯导信息确定需要操作的径向线数，实现雷达非均匀成像情形下的动态校正。
步骤2、针对校正后的雷达图像，确定雷达图像中的反演区域，重构海面波高程；
步骤3、对重构的海面波高程进行功率谱密度估计，通过波谱密度和峰值波数实现海浪参数的估计。
本实施方式中，当船舶存在横摇和纵摇运动时，不会改变艏向方向，因此不会造成雷达与海域径向线之间的差异。本实施方式针对横摇对雷达径向海浪信息的影响，并提出了各运动状态下的雷达图像动态校正方法，降低了所选反演区域的偏移程度，保证了雷达图像提供的海浪信息质量，再利用海浪反演技术实现海浪信息的精确估计，本实施方式的动态校正方案能够有效降低载体晃动所造成的反演误差，提高反演精度，解决了雷达由于姿态变化使非快拍成像导致的图像畸变的问题。
现有谱分析等反演方法对非平稳、非均匀的海浪信号处理难以取得精确的波高估计结果，且依赖浮标现场仪器校准线性拟合参数，不便于船载条件下的应用。为此，本实施方式提出了一种基于二维连续小波变换的海浪反演算法，能在不借助额外校准信息的情况下，实现海浪波高、浪向、波周期等参数的精确估计。本实施方式的步骤2中，通过二维连续小波变换，在反演区域的波数域中去除与海浪不相关的噪声分量，应用连续小波变换的逆过程重构海面波高程，具体包括：
首先，为了选取海杂波强度图像中合适的反演区域，减少回波数据的不确定性，需要粗估计主波波向，然后通过判断径向波束是否具有灰度最大方差值实现波向的检索。针对雷达回波强度随距离增大存在衰减的问题，需要对雷达图像进行校正，提高反演结果的可信度。雷达回波强度随距离在距离向上的关系为：
σ0∝R-n (4)
其中σ0为归一化雷达截面，R为海表面入射点到雷达天线的距离，n一般取2～4。考虑雷达图像中的阴影特征，采用低通滤波器实现对雷达图像的增强，以减小海面杂波图像获取过程中由于阴影效应引入的失真。滤波方法采用零相位巴特沃斯低通滤波器。根据最大方差波束检索过程中得到的波束k，能够确定反演区域的方位A0。由于雷达图像中的海杂波信号随距离的增大呈非线性衰减，当掠射角度在1～10°范围内时，主要受阴影调制作用的影响，且图像中包含了显著的海浪纹理特征的信号。考虑雷达量程、架设高度等因素，需要在合适的掠射角范围内选取R0，由此可以确定从雷达图像中选取的反演区域位置为(A0,R0)。
海浪的平均波周期大约在5～11s之间，在保证反演精度的前提下，分析雷达图像序列时在空间上通常要包含3～5个海浪波形。选取反演区域时，还需要将雷达极坐标系图像转换为笛卡尔坐标系图像，其大小一般选取笛卡尔坐标系下128*128的像素点，中心坐标为(R0 sinA0,R0cosA0)。
由于近岸条件下的海浪信号在空间和时间上显著变化，传统傅立叶分析方法的应用会受均匀性和平稳性限制。为了从非均匀波场图像中提取精确的谱，本文利用二维连续小波变换(2-D CWT)提取海洋遥感图像I(r)的频谱W，以揭示非均匀海浪特征：
其中，r表示以雷达位置为原点的平面位置矢量，r＝(x,y)，x,y分别表示以东、北为正向的坐标，I(r)表示海洋遥感图像中r处的回波强度，表示母小波ψ(r)的缩放因子，/>为小波平移向量，/>是小波旋转角为/>的旋转矩阵；在海浪线性水波理论应用中，通常选择Morlet小波作为母小波，它是空间域中的一种衰减正弦函数：
其中，k0表示中心波数，ε对应于各向异性参数，σ表示形状参数。一般在应用中默认为(k0,σ,ε)＝(6,1,1)。在计算过程中，可以利用具有傅里叶变换的二维连续小波变换(2-DCWTFT)实现2-D CWT，通过二维连续小波变换将反演区域的海浪信号I转换至波数域W进行处理：
Cψ是与所用小波有关的容许常数，分别表示波数域下的二维雷达图像和母小波函数，波数向量k＝(kx,ky)，kx,ky分别表示k在x，y方向上的分量；在波数域下加入高斯滤波器/>和Gabor滤波器/>去除与海浪不相关的噪声分量：
其中，是高斯滤波器x、y方向上的方差，/>是Gabor滤波器的方差；/> 分别表示峰值波数在x、y方向上的分量；
利用二维高斯高通滤波器能够去除与海浪不相关的低波数分量，并消除由于距离向平均雷达截面衰减而出现在f＝0Hz附近的峰值频谱强度。通过Gabor带通滤波器能够去除波谱区域以外的分量，增强图像纹理特征，提高重构后海面波高程的波浪纹理识别度。利用滤波后的方向频谱中出现的涌浪峰值，能够得到峰值波数kp，进而得到海浪波峰波长λp＝2π/kp和波峰浪向θp。
为了得到海面波高程图像η(r,t)，需要应用连续小波变换的逆过程实现重构。若根据上述过程中所用的Morlet小波函数直接进行2-D CWT逆变换，将存在相当大的复杂性，且信号重构误差大，但可以采用狄拉克函数(δ)作为2-D CWT逆变换的母小波函数，代替Morlet小波函数，应用连续小波变换的逆过程重构海面波高程
其中，Cψ,δ表示以狄拉克函数为母小波实现逆CWT的容许常数，表示在r处波数域下的海浪信号；
需要注意的是，此处重构的海面波高程图像对应的是真实海面波高程η(r,t)的灰度值，即，描述海杂波数据的电磁回波强度，而非海面的波浪高度，但它可以通过尺度变换得到真实海面波高程图。记/>表示与雷达相距R0处的重构海面波高程，根据/>确定雷达相距R0处的重构海面波高程η(R0,t)：
参数G与反演区域及雷达性能指标有关，定义为：
其中，Φ、分别表示反演区域掠射角度、雷达天线水平波束宽度、重构波高程图像的平均值。
为了从η(R0,t)图像中提取有效波高，本实施方式的步骤3用Welch修正周期图法对重构的海面波高程η(R0,t)的时间序列进行功率谱密度(PSD)估计，对每一段数据做加窗处理，降低PSD中不同谱峰的相互干扰，通过波谱密度和峰值波数实现海浪参数的估计。然而，较短的时间序列得到的功率谱估计结果可能并不理想，影响有效波高的反演性能，因此，需要选取相同反演区域下的连续图像序列，通常设置图像序列长度N＝32。将时间序列η(R0,t)每组128个样本分成m个相同大小的重叠Hamming窗来计算波谱密度S(f)，可以得到波峰周期Tp＝1/fp和有效波高Hs。工程上一般采用实现海浪有效波高估计，可通过角频域转换得到海浪频谱S(ω)，以及频谱能量E：
上述过程中，基于滤波和信号处理，利用航海雷达图像中的海杂波信号幅值反演海面波高程，实现了对有效波高的估计。由于没有用到表面波线性色散关系式，因此不需要浮标测量值作为外部参数校准。本实施方式的海浪参数包括海浪波峰波长、波峰浪向θp、波高和波峰周期：
海浪波峰波长λp＝2π/kp
波高
波峰周期Tp＝1/fp
本实施方式一方面，将惯性测量单元中获取的惯导信息绑定到雷达径向数据，分割图像区域并采用区域滑窗的方式转换径向序数差，构建数据-图像各个径向映射关系，实现雷达图像动态校正。另一方面，设计了基于二维连续小波变换的海浪反演算法，通过连续小波变换将其转换至波数域，利用高斯和Gabor滤波器去除与海浪不相关的低波数分量，经过小波变换逆过程重构海面波高程灰度图像。考虑雷达性能参数对其进行校正，得到真实海面波高程，最后利用Welch修正周期图法对波高程时间序列进行谱估计，提取波高、浪向、波周期等信息。所提方法提高了船载导航雷达海浪信息的反演精度。
具体实施例：船载雷达海洋信息反演方法流程如下：
本实施例的效果可通过以下实验验证说明。
①需要验证在非快拍成像条件下对雷达图像动态校正的合理性。
图4中呈现了对艏向非匀速、非定向旋转条件下的海域径向Xsea的变化分析。通过对单帧雷达图像采集过程中Xsea随艏向的变化可以看出，可将雷达图像分割为子区域校正，且当惯导信息频率较高时，每个子区域内的艏向转动更趋近于匀速变化。经实验验证，惯导信息频率为22条/帧时可取得较好的校正效果，因此针对雷达图像的动态校正方法是合理的。
②下述实验是为了验证本实施例所提的雷达图像动态校正方法及改进的海浪信息反演算法在实际应用过程中的有效性。
为进一步验证动态校正的应用效果，图10绘制了校正前(第一列)和校正后(第二列)的雷达图像反演波高对比结果。通过实验结果可以看出，分别采用谱分析、阴影抑制和本文算法对校正后的雷达图像进行处理，均获得了更好的反演效果，降低了反演误差。三种算法之间的对比可以看出，本实施例方法在相关系数方面有一定的提升，降低了波高反演的均方误差和平均偏差，能够实现海浪有效波高的准确估计。同时图11中呈现了本实施例方法在波向和波周期方面与浮标测量值之间的对比，可以看出本实施例方法具有较高的反演精度，最终实现了对有效波高、峰值波向和峰值波周期的精确估计，可以验证本实施例海浪信息反演方法的有效性。
虽然在本文中参照了特定的实施方式来描述本发明，但是应该理解的是，这些实施例仅仅是本发明的原理和应用的示例。因此应该理解的是，可以对示例性的实施例进行许多修改，并且可以设计出其他的布置，只要不偏离所附权利要求所限定的本发明的精神和范围。应该理解的是，可以通过不同于原始权利要求所描述的方式来结合不同的从属权利要求和本文中所述的特征。还可以理解的是，结合单独实施例所描述的特征可以使用在其他所述实施例中。
Claims (9)
1.船载雷达海浪信息反演方法，其特征在于，所述方法包括：
S1、建立非快拍成像校正模型：
在当前帧雷达扫描成像期间，若船舶艏向按顺时针匀速转过的角度为θ1，确定k1，在当前帧雷达图像[1,M-k1]的雷达径向线范围内，通过均匀插值k1条相邻径向线展开图像中的海浪信息；
若船舶艏向按逆时针匀速转过的角度为确定k2，在当前帧雷达图像[1,M]的雷达径向线范围内，通过均匀抽取k2条径向线，同时取前一帧雷达图像中[M-k2+1,M]范围内的径向线，补充当前帧雷达图像缺失的径向线；
若船舶艏向非匀速转动，将当前帧雷达图像分割成若干扇形子图区域进行滑窗校正，将每个子区域视为匀速转动，根据顺时针或逆时针的情况，在各子区域内采用均匀插值或抽取径向线的方式降低径向累计偏移，通过惯导信息确定需要操作的径向线数；
雷达图像为M个圆心角/>相等、半径相同的扇形所构成的环形图像；
S2、针对校正后的雷达图像，确定雷达图像中的反演区域，重构海面波高程；
S3、对重构的海面波高程进行功率谱密度估计，通过波谱密度和峰值波数实现海浪参数的估计。
2.根据权利要求1所述的船载雷达海浪信息反演方法，其特征在于，S2中，通过二维连续小波变换，在反演区域的波数域中去除与海浪不相关的噪声分量，应用连续小波变换的逆过程重构海面波高程。
3.根据权利要求2所述的船载雷达海浪信息反演方法，其特征在于，通过二维连续小波变换将反演区域的海浪信号I转换至波数域W进行处理：
其中，r表示以雷达位置为原点的平面位置矢量，r＝(x,y)，x,y分别表示以东、北为正向的坐标，I(r)表示海洋遥感图像中r处的回波强度，表示母小波ψ(r)的缩放因子，为小波平移向量，/>是小波旋转角为/>的旋转矩阵，Cψ是与所用小波有关的容许常数，/>分别表示波数域下的二维雷达图像和母小波函数，波数向量k＝(kx,ky)，kx,ky分别表示k在x，y方向上的分量；
在波数域下加入高斯滤波器和Gabor滤波器/>去除与海浪不相关的噪声分量：
其中，是高斯滤波器x、y方向上的方差，/>是Gabor滤波器的方差；/> 分别表示峰值波数在x、y方向上的分量；
应用连续小波变换的逆过程重构海面波高程
其中，Cψ,δ表示以狄拉克函数为母小波实现逆CWT的容许常数，表示在r处波数域下的海浪信号；
根据确定雷达相距R0处的重构海面波高程η(R0,t)：
其中，Φ、/>分别表示反演区域掠射角度、雷达天线水平波束宽度、重构波高程图像的平均值。
4.根据权利要求3所述的船载雷达海浪信息反演方法，其特征在于，所述母小波为Morlet小波：
其中，k0表示中心波数，ε对应于各向异性参数，σ表示形状参数。
5.根据权利要求3所述的船载雷达海浪信息反演方法，其特征在于，所述母小波函数为函数为狄拉克函数(δ)。
6.根据权利要求1所述的船载雷达海浪信息反演方法，其特征在于，S3中，用Welch修正周期图法对重构的海面波高程进行功率谱密度估计，对每一段数据做加窗处理，通过波谱密度和峰值波数实现海浪参数的估计。
7.根据权利要求6所述的船载雷达海浪信息反演方法，其特征在于，海浪参数包括海浪波峰波长、波峰浪向、波高和波峰周期。
8.一种计算机可读的存储设备，所述存储设备存储有计算机程序，其特征在于，所述计算机程序被执行时实现如权利要求1至7任一所述船载雷达海浪信息反演方法。
9.一种船载雷达海浪信息反演装置，包括存储设备、处理器以及存储在所述存储设备中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序实现如权利要求1至7任一所述船载雷达海浪信息反演方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311129796.6A CN117169882A (zh) | 2023-09-04 | 2023-09-04 | 船载雷达海浪信息反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311129796.6A CN117169882A (zh) | 2023-09-04 | 2023-09-04 | 船载雷达海浪信息反演方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117169882A true CN117169882A (zh) | 2023-12-05 |
Family
ID=88938981
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311129796.6A Pending CN117169882A (zh) | 2023-09-04 | 2023-09-04 | 船载雷达海浪信息反演方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117169882A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117554920A (zh) * | 2024-01-11 | 2024-02-13 | 之江实验室 | 一种水面探测方法、装置、存储介质及电子设备 |
-
2023
- 2023-09-04 CN CN202311129796.6A patent/CN117169882A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117554920A (zh) * | 2024-01-11 | 2024-02-13 | 之江实验室 | 一种水面探测方法、装置、存储介质及电子设备 |
| CN117554920B (zh) * | 2024-01-11 | 2024-04-02 | 之江实验室 | 一种水面探测方法、装置、存储介质及电子设备 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Dong et al. | A novel compressive sensing algorithm for SAR imaging | |
| CN111142105B (zh) | 复杂运动目标isar成像方法 | |
| CN109116352B (zh) | 一种圆扫isar模式船只超分辨率成像方法 | |
| CN103941257B (zh) | 一种基于波数能量谱的导航雷达图像反演海面风向的方法 | |
| CN111665489B (zh) | 一种基于目标特性的线谱提取方法 | |
| CN117169882A (zh) | 船载雷达海浪信息反演方法 | |
| CN111257886A (zh) | 一种利用单幅船载x波段雷达图像反演海浪参数的方法 | |
| CN112415515B (zh) | 一种机载圆迹sar对不同高度目标分离的方法 | |
| Marques et al. | Moving targets processing in SAR spatial domain | |
| CN113589287B (zh) | 合成孔径雷达稀疏成像方法、装置、电子设备及存储介质 | |
| CN113514833B (zh) | 一种基于海浪图像的海面任意点波向的反演方法 | |
| Huang et al. | Ocean remote sensing using X-band shipborne nautical radar—Applications in eastern Canada | |
| Kusk et al. | SAR focusing of P-band ice sounding data using back-projection | |
| CN112255607B (zh) | 一种海杂波的抑制方法 | |
| Lei et al. | A feature enhancement method based on the sub-aperture decomposition for rotating frame ship detection in SAR images | |
| CN107422320B (zh) | 一种消除降雨对x波段雷达观测海浪的影响的方法 | |
| CN106842201A (zh) | 一种基于序列图像的舰船目标isar交叉像判别方法 | |
| CN113406634B (zh) | 一种基于时域相位匹配的空间高速自旋目标isar三维成像方法 | |
| CN113204020B (zh) | 一种基于频谱分割的海浪波谱仪斑点噪声谱估计方法 | |
| CN115902791A (zh) | 基于s波段测波雷达时间多普勒谱的海浪反演方法及系统 | |
| CN115032601A (zh) | 一种基于空时联合滤波技术抑制图像序列中海杂波的航海雷达目标检测算法 | |
| CN114252878A (zh) | 一种基于逆合成孔径雷达对运动目标进行成像及横向定标的方法 | |
| Xu et al. | A novel SAR imaging method based on morphological component analysis | |
| CN107450058B (zh) | 基于FrFT和HT的雷达信号时频参数估计方法 | |
| Liu et al. | Sandglass transformation for synthetic aperture radar detection and imaging of ship at low signal-to-clutter-plus-noise ratio |
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 |