CN113917490A - 激光测风雷达信号去噪方法及装置 - Google Patents

激光测风雷达信号去噪方法及装置 Download PDF

Info

Publication number
CN113917490A
CN113917490A CN202111067340.2A CN202111067340A CN113917490A CN 113917490 A CN113917490 A CN 113917490A CN 202111067340 A CN202111067340 A CN 202111067340A CN 113917490 A CN113917490 A CN 113917490A
Authority
CN
China
Prior art keywords
signal
modal
decomposition
component
laser radar
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
CN202111067340.2A
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.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
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 Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202111067340.2A priority Critical patent/CN113917490A/zh
Publication of CN113917490A publication Critical patent/CN113917490A/zh
Pending legal-status Critical Current

Links

Images

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
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88—Lidar systems specially adapted for specific applications
    • G01S17/95—Lidar systems specially adapted for specific applications for meteorological use
    • 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/48—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00
    • G01S7/4802—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • 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)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Electromagnetism (AREA)
  • Optical Radar Systems And Details Thereof (AREA)

Abstract

本发明公开了一种基于奇异值分解和变分模态分解的激光雷达回波信号去噪方法及装置，首先将带噪信号利用变分模态分解进行分解，通过各模态分量之间的互相关系数找出噪声所在的模态并去除，从而保留大部分的有效信息，再将剩余模态利用SVD进一步去噪，提高输出信号的信噪比，最后将经过SVD去噪后的模态进行重构，获得去噪信号。本发明能够有效提高远程相干测风激光雷达的探测距离，并保证其风速估计精度，以更低的脉冲累加数实现相同的探测距离从而提高相干测风激光雷达的时间分辨率。

Description

激光测风雷达信号去噪方法及装置
技术领域
本发明涉及光电技术领域，具体而言，涉及一种激光测风雷达信号去噪方法及装置。
背景技术
相干多普勒激光测风雷达利用激光的多普勒效应，将雷达回波信号和本征信号进行拍频从而计算其多普勒频移信息，进而反演出大气的风场信息。要获得远距离大气风场信息的关键是从若信噪比信号中提取有效信息，提高信号信噪比的方法有傅里叶变换、小波变换、奇异值分解、经验模态分解和变分模态分解等。
对激光雷达回波信号进行风速反演一般要使用多脉冲累加技术提高回波信号信噪比，但脉冲累加数会影响激光雷达时间分辨率。时间分辨率Δt和激光雷达光源的脉冲重复频率f以及脉冲累加数N有关，Δt＝N/f。
目前，科研人员已经利用一些改进的去噪算法和联合去噪算法实现了对雷达信号信噪比的提高，然而，当雷达信号信噪比很低，dB数为负时，这些算法不再适用。因此，需要开发一种应用于弱信噪比信号的，提高远程相干测风雷达探测距离的激光测风雷达信号去噪方法。
发明内容
本发明提供了一种可提高弱激光雷达信号信噪比的去噪方法，从而提高远程相干测风雷达探测距离。
为解决上述问题，本发明提供一种基于奇异值分解和变分模态分解的激光雷达回波信号去噪方法，包括以下步骤：
步骤1：通过信号采集卡采集原始激光雷达回波信号s(t)；
步骤2：利用预先设定的VMD分解参数，对原始激光雷达回波信号s(t)进行VMD分解，得到多个模态分量信号；
步骤3：计算步骤2得到的模态分量信号之间的互相关系数，筛选出信号模态分量和噪声模态分量；
步骤4：应用奇异值分解方法对筛选出的信号模态分量进一步降噪，之后将进一步降噪之后的信号模态分量叠加，进行回波信号的重构，得到去噪处理后的激光雷达回波信号。
可选地，对原始激光雷达回波信号s(t)进行VMD分解的具体步骤为：
VMD分解的核心是公式(1)所示的变分问题的构造和求解：
Figure BDA0003258901880000021
其中，K是需要分解的模态个数，{uk}、{ωk}分别对应分解后第k个模态分量和中心频率；δ(t)为狄拉克函数，
Figure BDA0003258901880000022
表示括号内的式子对时间t求导；j为虚数单位；引入拉格朗日算子λ，得到公式(2)所描述的增广拉格朗日表达式：
Figure BDA0003258901880000023
其中，α为二次惩罚因子，可降低高斯噪声的干扰；uk为第K个模态分量信号频率、ωk为第K个模态分量信号中心频率、λk为第K个模态分量信号对应的拉格朗日算子；
利用交替方向乘子法，通过搜寻增广拉格朗日函数的鞍点，交替寻优迭代
Figure BDA0003258901880000024
其具体步骤包括：
首先，预先指定最大迭代次数N，并初始化分量频率
Figure BDA0003258901880000025
及其对应的初始中心频率
Figure BDA0003258901880000026
以及初始拉格朗日乘数λ1；
然后根据公式(3)、(4)和(5)更新
Figure BDA0003258901880000031
和ωk；
Figure BDA0003258901880000032
Figure BDA0003258901880000033
Figure BDA0003258901880000034
其中
Figure BDA0003258901880000035
是第k个模态分量信号频率，ωk是第k个模态分量信号中心频率，
Figure BDA0003258901880000036
是原始激光雷达回波信号的频谱，γ为拉格朗日乘数的更新参数，取γ＝10-3；
最后判断更新迭代后的分量频率是否满足公式(6)所述收敛条件，若不满足，则继续更新，若满足，则迭代结束，进而得到VMD分解后的各分量信号；
Figure BDA0003258901880000037
式中，ε为收敛容差，取ε＝5×10-8；
得到的uk(k＝1,2,...,K)即为分解到的各个模态信号。
可选地，所述计算步骤2得到的模态分量信号之间的互相关系数，筛选出信号模态分量和噪声模态分量，具体步骤为：
根据VMD分解得到的模态分量信号，根据公式(7)计算各模态之间的互相关系数：
Figure BDA0003258901880000038
其中x(n)表示原始信号，blimfi(n)表示第i个模态分量，N为数据采样点数，
Figure BDA0003258901880000041
为原始信号均值，
Figure BDA0003258901880000042
为各模态分量均值；
根据计算得到的互相关系数，计算互相关系数的差分值，找到差分值中第一个局部极大值点对应的下标，该处即为高频噪声与信号的分解，将此前的分量滤除。
可选地，所述应用奇异值分解方法对筛选出的信号模态分量进一步降噪，之后将进一步降噪之后的信号模态分量叠加，进行回波信号的重构，得到去噪处理后的激光雷达回波信号，具体步骤为：
对于一长度为M的激光雷达信号分量，根据公式(8)构造其Hankel矩阵：
Figure BDA0003258901880000043
其中n＝(M+1)/2，xi(i＝1,2,...,M)为激光雷达信号值；
对根据(8)式构造的矩阵按公式(9)进行分解：
A＝UWVT (9)
令m＝M-n+1，则A的维度为m×n，则U、W、V的维度分别是U∈Rm×m，W∈Rm×n，V∈Rn×n，其中U、V分别是左右奇异矩阵，都是单位正交矩阵，W可以表示为公式(10)：
Figure BDA0003258901880000044
其中S是一个对角矩阵，S＝diag(s1,s2,…sn)；s1,s2…sn称为奇异值，且s1≥s2≥…sn≥0；
寻找奇异值序列的拐点，设拐点的序号为m，则前m个奇异值对应主要反映有信号，后n-m个奇异值主要反映噪声；将后n-m个奇异值置为零；
将置零后的奇异值矩阵还原为Hankel矩阵，从而得到经奇异值分解后的激光雷达信号分量；
将每一个经过奇异值分解去噪后的模态分量信号进行叠加，得到最终去噪激光雷达信号。
本发明提供一种基于奇异值分解和变分模态分解的激光雷达回波信号去噪装置，包括：
采集模块，用于通过信号采集卡采集原始激光雷达回波信号s(t)；
分解模块，用于利用预先设定的VMD分解参数，对原始激光雷达回波信号s(t)进行VMD分解，得到多个模态分量信号；
筛选模块，用于计算所述模态分量信号之间的互相关系数，筛选出信号模态分量和噪声模态分量；
去噪模块，用于应用奇异值分解方法对筛选出的信号模态分量进一步降噪，之后将进一步降噪之后的信号模态分量叠加，进行回波信号的重构，得到去噪处理后的激光雷达回波信号。
可选地，所述分解模块具体用于：
VMD分解的核心是公式(1)所示的变分问题的构造和求解：
Figure BDA0003258901880000051
其中，K是需要分解的模态个数，{uk}、{ωk}分别对应分解后第k个模态分量和中心频率；δ(t)为狄拉克函数，
Figure BDA0003258901880000052
表示括号内的式子对时间t求导；j为虚数单位；引入拉格朗日算子λ，得到公式(2)所描述的增广拉格朗日表达式：
Figure BDA0003258901880000061
其中，α为二次惩罚因子，可降低高斯噪声的干扰；uk为第K个模态分量信号频率、ωk为第K个模态分量信号中心频率、λk为第K个模态分量信号对应的拉格朗日算子；
利用交替方向乘子法，通过搜寻增广拉格朗日函数的鞍点，交替寻优迭代
Figure BDA0003258901880000062
其具体步骤包括：
首先，预先指定最大迭代次数N，并初始化分量频率
Figure BDA0003258901880000063
及其对应的初始中心频率
Figure BDA0003258901880000064
以及初始拉格朗日乘数λ1；
然后根据公式(3)、(4)和(5)更新
Figure BDA0003258901880000065
和ωk；
Figure BDA0003258901880000066
Figure BDA0003258901880000067
Figure BDA0003258901880000068
其中
Figure BDA0003258901880000069
是第k个模态分量信号频率，ωk是第k个模态分量信号中心频率，
Figure BDA00032589018800000610
是原始激光雷达回波信号的频谱，γ为拉格朗日乘数的更新参数，取γ＝10-3；
最后判断更新迭代后的分量频率是否满足公式(6)所述收敛条件，若不满足，则继续更新，若满足，则迭代结束，进而得到VMD分解后的各分量信号；
Figure BDA00032589018800000611
式中，ε为收敛容差，取ε＝5×10-8；
得到的uk(k＝1,2,...,K)即为分解到的各个模态信号。
可选地，所述筛选模块具体用于：
根据VMD分解得到的模态分量信号，根据公式(7)计算各模态之间的互相关系数：
Figure BDA0003258901880000071
其中x(n)表示原始信号，blimfi(n)表示第i个模态分量，N为数据采样点数，
Figure BDA0003258901880000072
为原始信号均值，
Figure BDA0003258901880000073
为各模态分量均值；
根据计算得到的互相关系数，计算互相关系数的差分值，找到差分值中第一个局部极大值点对应的下标，该处即为高频噪声与信号的分解，将此前的分量滤除。
可选地，所述去噪模块具体用于：
对于一长度为M的激光雷达信号分量，根据公式(8)构造其Hankel矩阵：
Figure BDA0003258901880000074
其中n＝(M+1)/2，xi(i＝1,2,...,M)为激光雷达信号值；
对根据(8)式构造的矩阵按公式(9)进行分解：
A＝UWVT (9)
令m＝M-n+1，则A的维度为m×n，则U、W、V的维度分别是U∈Rm×m，W∈Rm×n，V∈Rn×n，其中U、V分别是左右奇异矩阵，都是单位正交矩阵，W可以表示为公式(10)：
Figure BDA0003258901880000081
其中S是一个对角矩阵，S＝diag(s1,s2,…sn)；s1,s2…sn称为奇异值，且s1≥s2≥…sn≥0；
寻找奇异值序列的拐点，设拐点的序号为m，则前m个奇异值对应主要反映有信号，后n-m个奇异值主要反映噪声；将后n-m个奇异值置为零；
将置零后的奇异值矩阵还原为Hankel矩阵，从而得到经奇异值分解后的激光雷达信号分量；
将每一个经过奇异值分解去噪后的模态分量信号进行叠加，得到最终去噪激光雷达信号。
本发明与现有技术相比，具有以下优点和效果：
(1)采用互相关系数的差分值作为选择标准，对VMD分解后的模态进行筛选，提出了噪声相关模态，为奇异值分解提供了噪声含量较低的样本；
(2)应用奇异值分解方法进一步滤除信号相关模态的噪声，提高了去噪信号的信噪比；
(3)本发明较目前应用于激光雷达回波信号的去噪方法，对低信噪比信号有较强的去噪能力，能够在最大程度保留原始信号有效信息的同时提高信噪比，具有很好的价值和应用前景。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
图1为本发明实施例提供的基于奇异值分解和变分模态分解的激光雷达回波信号去噪方法的流程图；
图2为本发明实施例提供的带噪信号在50、100个300个脉冲累加数下的风速反演结果和探测距离的结果图；
图3为本发明实施例提供的去噪信号在50、100个300个脉冲累加数下的风速反演结果和探测距离的结果图；
图4为本发明实施例提供的去噪前后不同脉冲累加数反演的风速结果和标准风速的差值；
图5为本发明实施例提供的去噪后300个脉冲累加数反演的风速结果和标准风速的差异统计结果；
图6为本发明实施例提供的带噪信号在50和100个脉冲累加数下反演的风速结果和标准风速的差异统计结果；
图7为本发明实施例提供的去噪信号在50和100个脉冲累加数下反演的风速结果和标准风速的差异统计结果。
具体实施方式
为使本发明的上述目的、特征和优点能够更为明显易懂，下面结合附图对本发明的具体实施例做详细的说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
奇异值分解(Singular Value Decomposition,SVD)是一种矩阵分析方法，对信号进行奇异值分解，原始信号的大部分能量集中在奇异值较大处，噪声则对应较小的奇异值，通过对噪声对应的奇异值置零再将矩阵还原则可以实现信号的去噪，该方法可以对大部分噪声进行有效去除，但去噪信号和原始信号一致性较差；变分模态分解(Variational ModeDecomposition,VMD)是一种模态分解方法，该方法将非线性非稳态信号分解为几个模态分量，通过将噪声对应模态分量去除实现信号去噪，但保留的有效信号部分还存在部分噪声。
本发明提供了一种基于奇异值分解和变分模态分解的激光雷达回波信号去噪方法，可提高弱激光雷达信号信噪比，从而提高远程相干测风雷达探测距离。
参见图1所示的一种基于奇异值分解和变分模态分解的激光雷达回波信号去噪方法的流程示意图，包括以下步骤：
步骤1，通过信号采集卡采集原始激光雷达回波信号s(t)。
步骤2，设定VMD分解参数。
步骤3，利用预先设定的VMD分解参数，对原始激光雷达回波信号s(t)进行VMD分解，得到多个模态分量信号。
步骤4，计算步骤3得到的模态分量信号之间的互相关系数，筛选出信号模态分量和噪声模态分量。
步骤5，应用奇异值分解方法对筛选出的信号模态分量进一步降噪，之后将进一步降噪之后的信号模态分量叠加，进行回波信号的重构，得到去噪处理后的激光雷达回波信号。
本发明实施例能够有效提高远程相干测风激光雷达的探测距离，并保证其风速估计精度，以更低的脉冲累加数实现相同的探测距离从而提高相干测风激光雷达的时间分辨率。
可选地，VMD分解的具体步骤为：
步骤3.1：VMD分解的核心是公式(1)所示的变分问题的构造和求解：
Figure BDA0003258901880000101
其中，K是需要分解的模态个数，{uk}、{ωk}分别对应分解后第k个模态分量和中心频率；δ(t)为狄拉克函数，
Figure BDA0003258901880000111
表示括号内的式子对时间t求导；j为虚数单位；引入拉格朗日算子λ，得到公式(2)所描述的增广拉格朗日表达式：
Figure BDA0003258901880000112
其中，α为二次惩罚因子，可降低高斯噪声的干扰；uk为第K个模态分量信号频率、ωk为第K个模态分量信号中心频率、λk为第K个模态分量信号对应的拉格朗日算子；
步骤3.2：利用交替方向乘子法，通过搜寻增广拉格朗日函数的鞍点，交替寻优迭代
Figure BDA0003258901880000113
其具体步骤包括：
首先，预先指定最大迭代次数N，并初始化分量频率
Figure BDA0003258901880000114
及其对应的初始中心频率
Figure BDA0003258901880000115
以及初始拉格朗日乘数λ1；
然后根据公式(3)、(4)和(5)更新
Figure BDA0003258901880000116
和ωk；
Figure BDA0003258901880000117
Figure BDA0003258901880000118
Figure BDA0003258901880000119
其中
Figure BDA00032589018800001110
是第k个模态分量信号频率，ωk是第k个模态分量信号中心频率，
Figure BDA00032589018800001111
是原始激光雷达回波信号的频谱，γ为拉格朗日乘数的更新参数，取γ＝10-3；
最后判断更新迭代后的分量频率是否满足公式(6)所述收敛条件，若不满足，则继续更新，若满足，则迭代结束，进而得到VMD分解后的各分量信号；
Figure BDA0003258901880000121
式中，ε为收敛容差，取ε＝5×10-8；
得到的uk(k＝1,2,...,K)即为分解到的各个模态信号。
可选地，步骤4的具体步骤为：
步骤4.1：根据VMD分解得到的模态分量信号，根据公式(7)计算各模态之间的互相关系数：
Figure BDA0003258901880000122
其中x(n)表示原始信号，blimfi(n)表示第i个模态分量，N为数据采样点数，
Figure BDA0003258901880000123
为原始信号均值，
Figure BDA0003258901880000124
为各模态分量均值；
步骤4.2：根据计算得到的互相关系数，计算互相关系数的差分值，找到差分值中第一个局部极大值点对应的下标，该处即为高频噪声与信号的分解，将此前的分量滤除。
可选地，步骤5的具体步骤为：
步骤5.1：对于一长度为M的激光雷达信号分量，根据公式(8)构造其Hankel矩阵：
Figure BDA0003258901880000125
其中n＝(M+1)/2，xi(i＝1,2,...,M)为激光雷达信号值；
步骤5.2：对根据(8)式构造的矩阵按公式(9)进行分解：
A＝UWVT (9)
令m＝M-n+1，则A的维度为m×n，则U、W、V的维度分别是U∈Rm×m，W∈Rm×n，V∈Rn×n，其中U、V分别是左右奇异矩阵，都是单位正交矩阵，W可以表示为公式(10)：
Figure BDA0003258901880000131
其中S是一个对角矩阵，S＝diag(s1,s2,…sn)；s1,s2…sn称为奇异值，且s1≥s2≥…sn≥0；
步骤5.3：寻找奇异值序列的拐点，设拐点的序号为m，则前m个奇异值对应主要反映有信号，后n-m个奇异值主要反映噪声；将后n-m个奇异值置为零；
步骤5.4：将置零后的奇异值矩阵还原为Hankel矩阵，从而得到经奇异值分解后的激光雷达信号分量；
步骤5.5：将每一个经过奇异值分解去噪后的模态分量信号进行叠加，得到最终去噪激光雷达信号。
本发明实施例采用互相关系数的差分值作为选择标准，对VMD分解后的模态进行筛选，提出了噪声相关模态，为奇异值分解提供了噪声含量较低的样本；应用奇异值分解方法进一步滤除信号相关模态的噪声，提高了去噪信号的信噪比；较目前应用于激光雷达回波信号的去噪方法，对低信噪比信号有较强的去噪能力，能够在最大程度保留原始信号有效信息的同时提高信噪比，具有很好的价值和应用前景。
本发明实施例提供了一种基于奇异值分解和变分模态分解的激光雷达回波信号去噪方法，包括：
第一步：由信号采集卡采集激光雷达回波信号，并反演计算风俗结果。图2为本发明实施例提供的50、100个300个脉冲累加数下的风速反演结果和探测距离的结果图。由图可知，脉冲累加数的提高可以提高激光雷达的探测距离，但脉冲累加数过高会降低激光雷达的时间分辨率。
第二步：初始化VMD算法参数，进行VMD分解。
第三步：对VMD分解得到的模态计算其互相关系数的差分值，找到噪声对应模态并去除；
第四步：对剩余模态使用奇异值分解进一步去噪，对去噪模态进行重构，得到去噪雷达信号；
第五步：对去噪信号使用不同脉冲累加数进行风速反演计算，如图3所示的去噪信号在50、100个300个脉冲累加数下的风速反演结果和探测距离的结果图。可以看出，该去噪方法在相同的脉冲累加数下能够有效提高雷达的探测距离。对于50个脉冲累加数，如果忽略在22.5km处的错误值，其探测距离可以从18.45km提升至24km。
第六步：以去噪前300个脉冲累加下得到的风速结果作为标准风速，将去噪前后不同脉冲累加数得到的风速和标准风速对比，分析其风速差异的均值和标准差。
图4为本发明实施例提供的去噪前后不同脉冲累加数反演的风速结果和标准风速的差值，图5为本发明实施例提供的去噪信号在300个脉冲累加数反演的风速结果和标准风速的差异统计结果。可以看出，该去噪方法能够有效还原有效信息，从而保证去噪信号风速反演的精度。
图6为本发明实施例提供的带噪信号在50和100个脉冲累加下反演的风速结果和标准风速差值的差异统计结果，图7为本发明实施例提供的去噪信号在50和100个脉冲累加数下反演的风速结果和标准风速差值的差异统计结果。可以看出，相同的脉冲累加下，去噪信号的标准差远小于带噪信号的标准差，这说明该方法可以在提高风速有效探测距离的同时保证风速估计精度。
由于激光雷达的时间分辨率Δt和激光光源的脉冲重复频率f以及风速计算的脉冲累加数N有关，Δt＝N/f，因此该方法能够以低脉冲累加获得同样的探测距离，提高激光雷达的时间分辨率。
本发明还提供了一种基于奇异值分解和变分模态分解的激光雷达回波信号去噪装置，包括：
采集模块，用于通过信号采集卡采集原始激光雷达回波信号s(t)；
分解模块，用于利用预先设定的VMD分解参数，对原始激光雷达回波信号s(t)进行VMD分解，得到多个模态分量信号；
筛选模块，用于计算所述模态分量信号之间的互相关系数，筛选出信号模态分量和噪声模态分量；
去噪模块，用于应用奇异值分解方法对筛选出的信号模态分量进一步降噪，之后将进一步降噪之后的信号模态分量叠加，进行回波信号的重构，得到去噪处理后的激光雷达回波信号。
可选地，所述分解模块具体用于：
VMD分解的核心是公式(1)所示的变分问题的构造和求解：
Figure BDA0003258901880000151
其中，K是需要分解的模态个数，{uk}、{ωk}分别对应分解后第k个模态分量和中心频率；δ(t)为狄拉克函数，
Figure BDA0003258901880000152
表示括号内的式子对时间t求导；j为虚数单位；引入拉格朗日算子λ，得到公式(2)所描述的增广拉格朗日表达式：
Figure BDA0003258901880000153
其中，α为二次惩罚因子，可降低高斯噪声的干扰；uk为第K个模态分量信号频率、ωk为第K个模态分量信号中心频率、λk为第K个模态分量信号对应的拉格朗日算子；
利用交替方向乘子法，通过搜寻增广拉格朗日函数的鞍点，交替寻优迭代
Figure BDA0003258901880000161
其具体步骤包括：
首先，预先指定最大迭代次数N，并初始化分量频率
Figure BDA0003258901880000162
及其对应的初始中心频率
Figure BDA0003258901880000163
以及初始拉格朗日乘数λ1；
然后根据公式(3)、(4)和(5)更新
Figure BDA0003258901880000164
和ωk；
Figure BDA0003258901880000165
Figure BDA0003258901880000166
Figure BDA0003258901880000167
其中
Figure BDA0003258901880000168
是第k个模态分量信号频率，ωk是第k个模态分量信号中心频率，
Figure BDA0003258901880000169
是原始激光雷达回波信号的频谱，γ为拉格朗日乘数的更新参数，取γ＝10-3；
最后判断更新迭代后的分量频率是否满足公式(6)所述收敛条件，若不满足，则继续更新，若满足，则迭代结束，进而得到VMD分解后的各分量信号；
Figure BDA00032589018800001610
式中，ε为收敛容差，取ε＝5×10-8；
得到的uk(k＝1,2,...,K)即为分解到的各个模态信号。
可选地，所述筛选模块具体用于：
根据VMD分解得到的模态分量信号，根据公式(7)计算各模态之间的互相关系数：
Figure BDA0003258901880000171
其中x(n)表示原始信号，blimfi(n)表示第i个模态分量，N为数据采样点数，
Figure BDA0003258901880000172
为原始信号均值，
Figure BDA0003258901880000173
为各模态分量均值；
根据计算得到的互相关系数，计算互相关系数的差分值，找到差分值中第一个局部极大值点对应的下标，该处即为高频噪声与信号的分解，将此前的分量滤除。
可选地，所述去噪模块具体用于：
对于一长度为M的激光雷达信号分量，根据公式(8)构造其Hankel矩阵：
Figure BDA0003258901880000174
其中n＝(M+1)/2，xi(i＝1,2,...,M)为激光雷达信号值；
对根据(8)式构造的矩阵按公式(9)进行分解：
A＝UWVT (9)
令m＝M-n+1，则A的维度为m×n，则U、W、V的维度分别是U∈Rm×m，W∈Rm×n，V∈Rn×n，其中U、V分别是左右奇异矩阵，都是单位正交矩阵，W可以表示为公式(10)：
Figure BDA0003258901880000175
其中S是一个对角矩阵，S＝diag(s1,s2,…sn)；s1,s2…sn称为奇异值，且s1≥s2≥…sn≥0；
寻找奇异值序列的拐点，设拐点的序号为m，则前m个奇异值对应主要反映有信号，后n-m个奇异值主要反映噪声；将后n-m个奇异值置为零；
将置零后的奇异值矩阵还原为Hankel矩阵，从而得到经奇异值分解后的激光雷达信号分量；
将每一个经过奇异值分解去噪后的模态分量信号进行叠加，得到最终去噪激光雷达信号。
本领域技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程度来指令控制装置来完成，所述的程序可存储于一计算机可读取的存储介质中，所述程序在执行时可包括如上述各方法实施例的流程，其中所述的存储介质可为存储器、磁盘、光盘等。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。
本说明书中所描述的以上内容仅仅是对本发明所作的举例说明。本发明不受上述实施例的限制，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (8)

1.一种基于奇异值分解和变分模态分解的激光雷达回波信号去噪方法，其特征在于，包括以下步骤：
步骤1：通过信号采集卡采集原始激光雷达回波信号s(t)；
步骤2：利用预先设定的VMD分解参数，对原始激光雷达回波信号s(t)进行VMD分解，得到多个模态分量信号；
步骤3：计算步骤2得到的模态分量信号之间的互相关系数，筛选出信号模态分量和噪声模态分量；
步骤4：应用奇异值分解方法对筛选出的信号模态分量进一步降噪，之后将进一步降噪之后的信号模态分量叠加，进行回波信号的重构，得到去噪处理后的激光雷达回波信号。
2.根据权利要求1所述的方法，其特征在于，对原始激光雷达回波信号s(t)进行VMD分解的具体步骤为：
VMD分解的核心是公式(1)所示的变分问题的构造和求解：
Figure FDA0003258901870000011
其中，K是需要分解的模态个数，{uk}、{ωk}分别对应分解后第k个模态分量和中心频率；δ(t)为狄拉克函数，
Figure FDA0003258901870000012
表示括号内的式子对时间t求导；j为虚数单位；引入拉格朗日算子λ，得到公式(2)所描述的增广拉格朗日表达式：
Figure FDA0003258901870000013
其中，α为二次惩罚因子，可降低高斯噪声的干扰；uk为第K个模态分量信号频率、ωk为第K个模态分量信号中心频率、λk为第K个模态分量信号对应的拉格朗日算子；
利用交替方向乘子法，通过搜寻增广拉格朗日函数的鞍点，交替寻优迭代
Figure FDA0003258901870000021
其具体步骤包括：
首先，预先指定最大迭代次数N，并初始化分量频率
Figure FDA0003258901870000022
及其对应的初始中心频率
Figure FDA0003258901870000023
以及初始拉格朗日乘数λ1；
然后根据公式(3)、(4)和(5)更新
Figure FDA0003258901870000024
和ωk；
Figure FDA0003258901870000025
Figure FDA0003258901870000026
Figure FDA0003258901870000027
其中
Figure FDA0003258901870000028
是第k个模态分量信号频率，ωk是第k个模态分量信号中心频率，
Figure FDA0003258901870000029
是原始激光雷达回波信号的频谱，γ为拉格朗日乘数的更新参数，取γ＝10-3；
最后判断更新迭代后的分量频率是否满足公式(6)所述收敛条件，若不满足，则继续更新，若满足，则迭代结束，进而得到VMD分解后的各分量信号；
Figure FDA00032589018700000210
式中，ε为收敛容差，取ε＝5×10-8；
得到的uk(k＝1,2,...,K)即为分解到的各个模态信号。
3.根据权利要求1或2所述的方法，其特征在于，所述计算步骤2得到的模态分量信号之间的互相关系数，筛选出信号模态分量和噪声模态分量，具体步骤为：
根据VMD分解得到的模态分量信号，根据公式(7)计算各模态之间的互相关系数：
Figure FDA0003258901870000031
其中x(n)表示原始信号，blimfi(n)表示第i个模态分量，N为数据采样点数，
Figure FDA0003258901870000032
为原始信号均值，
Figure FDA0003258901870000033
为各模态分量均值；
根据计算得到的互相关系数，计算互相关系数的差分值，找到差分值中第一个局部极大值点对应的下标，该处即为高频噪声与信号的分解，将此前的分量滤除。
4.根据权利要求3所述的方法，其特征在于，所述应用奇异值分解方法对筛选出的信号模态分量进一步降噪，之后将进一步降噪之后的信号模态分量叠加，进行回波信号的重构，得到去噪处理后的激光雷达回波信号，具体步骤为：
对于一长度为M的激光雷达信号分量，根据公式(8)构造其Hankel矩阵：
Figure FDA0003258901870000034
其中n＝(M+1)/2，xi(i＝1,2,...,M)为激光雷达信号值；
对根据(8)式构造的矩阵按公式(9)进行分解：
A＝UWVT (9)
令m＝M-n+1，则A的维度为m×n，则U、W、V的维度分别是U∈Rm×m，W∈Rm×n，V∈Rn×n，其中U、V分别是左右奇异矩阵，都是单位正交矩阵，W可以表示为公式(10)：
Figure FDA0003258901870000041
其中S是一个对角矩阵，S＝diag(s1,s2,…sn)；s1,s2...sn称为奇异值，且s1≥s2≥…sn≥0；
寻找奇异值序列的拐点，设拐点的序号为m，则前m个奇异值对应主要反映有信号，后n-m个奇异值主要反映噪声；将后n-m个奇异值置为零；
将置零后的奇异值矩阵还原为Hankel矩阵，从而得到经奇异值分解后的激光雷达信号分量；
将每一个经过奇异值分解去噪后的模态分量信号进行叠加，得到最终去噪激光雷达信号。
5.一种基于奇异值分解和变分模态分解的激光雷达回波信号去噪装置，其特征在于，包括：
采集模块，用于通过信号采集卡采集原始激光雷达回波信号s(t)；
分解模块，用于利用预先设定的VMD分解参数，对原始激光雷达回波信号s(t)进行VMD分解，得到多个模态分量信号；
筛选模块，用于计算所述模态分量信号之间的互相关系数，筛选出信号模态分量和噪声模态分量；
去噪模块，用于应用奇异值分解方法对筛选出的信号模态分量进一步降噪，之后将进一步降噪之后的信号模态分量叠加，进行回波信号的重构，得到去噪处理后的激光雷达回波信号。
6.根据权利要求5所述的装置，其特征在于，所述分解模块具体用于：
VMD分解的核心是公式(1)所示的变分问题的构造和求解：
Figure FDA0003258901870000042
其中，K是需要分解的模态个数，{uk}、{ωk}分别对应分解后第k个模态分量和中心频率；δ(t)为狄拉克函数，
Figure FDA0003258901870000051
表示括号内的式子对时间t求导；j为虚数单位；引入拉格朗日算子λ，得到公式(2)所描述的增广拉格朗日表达式：
Figure FDA0003258901870000052
其中，α为二次惩罚因子，可降低高斯噪声的干扰；uk为第K个模态分量信号频率、ωk为第K个模态分量信号中心频率、λk为第K个模态分量信号对应的拉格朗日算子；
利用交替方向乘子法，通过搜寻增广拉格朗日函数的鞍点，交替寻优迭代
Figure FDA0003258901870000053
其具体步骤包括：
首先，预先指定最大迭代次数N，并初始化分量频率
Figure FDA0003258901870000054
及其对应的初始中心频率
Figure FDA0003258901870000055
以及初始拉格朗日乘数λ1；
然后根据公式(3)、(4)和(5)更新
Figure FDA0003258901870000056
和ωk；
Figure FDA0003258901870000057
Figure FDA0003258901870000058
Figure FDA0003258901870000059
其中
Figure FDA00032589018700000510
是第k个模态分量信号频率，ωk是第k个模态分量信号中心频率，
Figure FDA00032589018700000511
是原始激光雷达回波信号的频谱，γ为拉格朗日乘数的更新参数，取γ＝10-3；
最后判断更新迭代后的分量频率是否满足公式(6)所述收敛条件，若不满足，则继续更新，若满足，则迭代结束，进而得到VMD分解后的各分量信号；
Figure FDA0003258901870000061
式中，ε为收敛容差，取ε＝5×10-8；
得到的uk(k＝1,2,...,K)即为分解到的各个模态信号。
7.根据权利要求5或6所述的装置，其特征在于，所述筛选模块具体用于：
根据VMD分解得到的模态分量信号，根据公式(7)计算各模态之间的互相关系数：
Figure FDA0003258901870000062
其中x(n)表示原始信号，blimfi(n)表示第i个模态分量，N为数据采样点数，
Figure FDA0003258901870000063
为原始信号均值，
Figure FDA0003258901870000064
为各模态分量均值；
根据计算得到的互相关系数，计算互相关系数的差分值，找到差分值中第一个局部极大值点对应的下标，该处即为高频噪声与信号的分解，将此前的分量滤除。
8.根据权利要求7所述的装置，其特征在于，所述去噪模块具体用于：
对于一长度为M的激光雷达信号分量，根据公式(8)构造其Hankel矩阵：
Figure FDA0003258901870000065
其中n＝(M+1)/2，xi(i＝1,2,...,M)为激光雷达信号值；
对根据(8)式构造的矩阵按公式(9)进行分解：
A＝UWVT (9)
令m＝M-n+1，则A的维度为m×n，则U、W、V的维度分别是U∈Rm×m，W∈Rm×n，V∈Rn×n，其中U、V分别是左右奇异矩阵，都是单位正交矩阵，W可以表示为公式(10)：
Figure FDA0003258901870000071
其中S是一个对角矩阵，S＝diag(s1,s2,…sn)；s1,s2...sn称为奇异值，且s1≥s2≥…sn≥0；
寻找奇异值序列的拐点，设拐点的序号为m，则前m个奇异值对应主要反映有信号，后n-m个奇异值主要反映噪声；将后n-m个奇异值置为零；
将置零后的奇异值矩阵还原为Hankel矩阵，从而得到经奇异值分解后的激光雷达信号分量；
将每一个经过奇异值分解去噪后的模态分量信号进行叠加，得到最终去噪激光雷达信号。
CN202111067340.2A 2021-09-13 2021-09-13 激光测风雷达信号去噪方法及装置 Pending CN113917490A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111067340.2A CN113917490A (zh) 2021-09-13 2021-09-13 激光测风雷达信号去噪方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111067340.2A CN113917490A (zh) 2021-09-13 2021-09-13 激光测风雷达信号去噪方法及装置

Publications (1)

Publication Number Publication Date
CN113917490A true CN113917490A (zh) 2022-01-11

Family

ID=79234716

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111067340.2A Pending CN113917490A (zh) 2021-09-13 2021-09-13 激光测风雷达信号去噪方法及装置

Country Status (1)

Country Link
CN (1) CN113917490A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114755654A (zh) * 2022-06-14 2022-07-15 中达天昇(江苏)电子科技有限公司 一种基于图像拟态技术的残损雷达信号修复方法
CN116541696A (zh) * 2023-07-07 2023-08-04 北京理工大学 一种脉冲体制引信回波信噪比估计方法
CN116631429A (zh) * 2023-07-25 2023-08-22 临沂金诺视讯数码科技有限公司 基于voip呼叫volte通话的语音视频处理方法及系统

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114755654A (zh) * 2022-06-14 2022-07-15 中达天昇(江苏)电子科技有限公司 一种基于图像拟态技术的残损雷达信号修复方法
CN116541696A (zh) * 2023-07-07 2023-08-04 北京理工大学 一种脉冲体制引信回波信噪比估计方法
CN116541696B (zh) * 2023-07-07 2023-09-19 北京理工大学 一种脉冲体制引信回波信噪比估计方法
CN116631429A (zh) * 2023-07-25 2023-08-22 临沂金诺视讯数码科技有限公司 基于voip呼叫volte通话的语音视频处理方法及系统
CN116631429B (zh) * 2023-07-25 2023-10-10 临沂金诺视讯数码科技有限公司 基于voip呼叫volte通话的语音视频处理方法及系统

Similar Documents

Publication Publication Date Title
CN113917490A (zh) 激光测风雷达信号去噪方法及装置
CN108845306B (zh) 基于变分模态分解的激光雷达回波信号去噪方法
CN111239697B (zh) 低秩矩阵分解的多维域联合sar宽带干扰抑制方法
CN103049892B (zh) 基于相似块矩阵秩最小化的非局部图像去噪方法
US20220020123A1 (en) Bayesian image denoising method based on distribution constraint of noisy images
CN110865357A (zh) 一种基于参数优化vmd的激光雷达回波信号降噪方法
CN103473755B (zh) 基于变化检测的sar图像稀疏去噪方法
CN103454622B (zh) 基于稀疏约束的宽带雷达目标复回波去噪方法
CN101685158B (zh) 基于隐马尔科夫树模型的sar图像去噪方法
CN103454621B (zh) 基于匹配追踪的宽带雷达目标复回波去噪方法
CN112327259B (zh) 一种sar图像中干扰信号的消除方法和装置
CN104515984A (zh) 基于贝叶斯压缩感知的宽带雷达目标复回波去噪方法
CN103077507B (zh) 基于Beta算法的多尺度SAR图像降噪方法
CN116699526A (zh) 一种基于稀疏与低秩模型的车载毫米波雷达干扰抑制方法
CN112578471A (zh) 一种探地雷达杂波噪声去除方法
CN111768349A (zh) 一种基于深度学习的espi图像降噪方法及系统
CN113238193B (zh) 一种多分量联合重构的sar回波宽带干扰抑制方法
CN110118958B (zh) 基于变分编码-解码网络的宽带雷达复回波去噪方法
Devapal et al. Discontinuity adaptive SAR image despeckling using curvelet-based BM3D technique
CN112014813B (zh) 一种天波雷达电离层污染校正方法及系统
CN107907542B (zh) 一种ivmd及能量估计相结合的dspi相位滤波方法
CN102509268B (zh) 基于免疫克隆选择的非下采样轮廓波域图像去噪方法
CN114966687A (zh) 基于低秩和非局部自相似的稀疏isar成像方法及系统
CN112230200A (zh) 一种基于激光雷达回波信号的改进型组合降噪方法
CN112363217A (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