一种激光雷达回波信号散斑噪声的滤除方法与流程

本发明属于激光数据处理技术领域，特别涉及了一种激光雷达回波信号散斑噪声的滤除方法。

背景技术：

激光雷达以激光为探测工具，以感光元件作为接收器件，以光电转换器作为信号转换的核心的雷达。得益于激光这种特殊波段的辐射源，具有更高分辨率、更强抗电子干扰及反隐身能力、更丰富测图信息等优势，已日渐形成其无可替代的地位。随着激光雷达技术的进一步发展，激光雷达在民用和军用领域都占有一席之地。比如我国自主研发的嫦娥三号探测器，应用了激光测距仪和ccd相机，其中激光测距仪主要靠提取激光回波时间来判断目标距离，ccd相机则根据目标反射光的强度来提取光强信息，二者相互结合实现三维成像。

激光雷达系统一般安装在飞行器或汽车等运载平台上，利用激光测量目标到系统之间的距离。其工作原理是激光器持续向目标发射激光脉冲，脉冲光束穿过大气与目标接触后产生微弱后向散射回波，接收端利用望远镜将其接收，再通过光电探测器转换为电脉冲信号，用数字采集器将其采样后进行处理。由此可见，只需要从采集的回波数字信号中提取出激光脉冲从发射到接收之间的飞行时间差δt，通过r＝cδt/2即可计算出目标点的距离值。基于飞行时间探测原理，激光雷达系统借助扫描器对目标区域不同点进行持续地扫描探测，最终获得目标区域的多个探测点距离信息，与ccd获取的强度信息融合后，得到目标的三维图像。

然而，激光雷达系统实际获得的回波信号并非理想，可直接用以定位回波脉冲位置。在实际传播过程中，回波信号受到诸多噪声影响，如大气湍流、大气衰减、目标交互散斑、背景光噪声、散粒噪声、放大器噪声及暗电流噪声等。按噪声属性分类，可分为乘性噪声和加性噪声两种。在光学成像中，乘性噪声表现为散斑现象，由于光束通过不均匀介质或在粗糙表面反射而产生。例如，由于大气微粒的随机运动改变了大气的折射率，使得光束沿非直线传播，产生了闪烁光斑，即散斑。大气湍流的每一个散斑元概率密度分布函数：

上式中，s为探测信号能量；sav为探测信号平均能量；为目标平均后的有效湍流波动量。这里其中为大气湍流的结构常数，r为系统到目标的距离；γtarget＝(ρl/reff)^7/3中的reff为有效光斑半径、ρl为(一维)有效湍流半径。当激光光束照射在目标表面发射漫反射时，漫反射表面会产生一系列独立并且随机相位的光束，当这些光束相互干涉就会产生强度随机分布的散斑。目标反射的散斑元的概率分布分布函数：

这里m为探测器可以探测到的散斑元数，表示为m＝drec/dsp；drec是接收口径；dsp＝λ/φ，λ为激光波长，φ为光束发散角；γ(m)为伽玛函数，即：

考虑到激光雷达测距系统在探测过程中，尤其在恶劣的探测环境下，回波信号受到各种散斑噪声干扰，严重影响了系统的测距精度，限制了激光雷达的应用，需要滤除激光雷达回波信号的散斑噪声。

技术实现要素：

为了解决上述背景技术提出的技术问题，本发明旨在提供一种激光雷达回波信号散斑噪声的滤除方法，巧妙利用激光发射信号作为参考信号进行搜索匹配，通过独立成分分析算法分解寻找最优匹配位置，从而完成回波脉冲信号与散斑噪声的有效分离。

为了实现上述技术目的，本发明的技术方案为：

一种激光雷达回波信号散斑噪声的滤除方法，包括以下步骤：

(1)对雷达系统采集的激光回波信号进行加性噪声滤除，再对滤噪后的回波信号进行同态变换，将同态变换后的信号作为独立成分分析的第一路观测信号；

(2)对参考信号进行时域循环移位，n为移位步长，且初始移位步长为1，将移位后的参考信号与滤除加性噪声的回波信号相乘后进行同态变换，将同态变换后的信号作为独立成分分析的第二路观测信号；

(3)将第一路观测信号和第二路观测信号输入独立成分分析算法，得到两路分离结果；

(4)对两路分离结果做同态逆变换，得到一路分离信号和一路分离噪声，计算分离噪声的局部均方差，并判断是否大于设定阈值，若是，则判断本次循环得到的分离信号不符合要求，令n＝n+1；若不是，则判断本次循环得到的分离信号符合要求，保存本次分离结果，令n＝n+1；

(5)判断移位步长n是否大于回波信号波形的时域长度，若不是，返回步骤(2)，继续迭代，若是，则停止迭代；

(6)选择所保存的各次分离噪声的局部均方差最小的那一次分离结果作为信号与噪声的最佳分离结果，完成激光雷达回波信号的散斑噪声滤除。

进一步地，在步骤(1)中，设雷达系统采集的激光回波信号为e(t)＝m(t)s(t-t0)+n(t)，其中，n(t)为加性噪声，s(t)为发射源脉冲信号，t0为激光实际飞行时间，m(t)为总的散斑噪声；则加性噪声滤除后的回波信号为为e′(t)＝m(t)s(t-t0)，e′(t)经同态变换后得到第一路观测信号x1(t)＝log(m(t)s(t-t0))＝logs(t-t0)+logm(t)。

进一步地，在步骤(2)中，采样发射源脉冲信号s(t)作为参考信号，设采样率为δt，时移τ＝nδt，则s(t)的移位波形为s(t-τ)，将s(t-τ)与第一路观测信号x1(t)相乘后进行同态变换，得到第二路观测信号x2(t)＝log(s(t-τ)·x1(t))＝log(s(t-τ)·m(t)s(t-t0))＝logs(t-τ)+logs(t-t0)+logm(t)。

进一步地，步骤(3)的具体过程如下：

将两路观测信号作为独立成分分析算法的两路输入源信号，其输出在τ取不同值时有不同结果：

当τ≠t0时，s(t-τ)≠s(t-t0)，此时参考信号未匹配到实际回波信号位置，观测向量表示为

当τ＝t0时，s(t-τ)＝s(t-t0)，此时参考信号匹配到实际回波信号位置，观测向量表示为

其中，a为混合矩阵，s为要分离的源矩阵；独立成分分析算法利用观测向量x找到解混矩阵w，使得输出矩阵y＝wx＝[y1(t)y2(t)]^t最佳逼近源矩阵s，从而找到两路分离结果y1(t)和y2(t)。

进一步地，步骤(4)的具体步骤如下：

(a)对独立成分分析算法输出的两路分离结果做同态逆变换，即分别为分离信号和分离噪声；

(b)对分离噪声在时域范围内找出其阈值标准以上的所有数据点为的平均值，α为预设参数，为阈值标准；

(c)计算y2′(n)的均方差，即的局部均方差：为y2′(n)长度，为y2′(n)平均值；

(d)若σ＞ε，ε为预设的误差阈值，则判断当前迭代步参考信号未匹配实际回波信号，舍弃本迭代步的分离结果，令移位步长n＝n+1；反之，则判断当前迭代步参考信号在误差ε范围内匹配到实际回波信号，将当前迭代步的分离结果保存，令移位步长n＝n+1。

采用上述技术方案带来的有益效果：

本发明充分利用了激光器发射的源脉冲信号，在移位搜索的过程中，通过独立成分分析算法的分离结果寻找到匹配实际回波信号的位置，进而分离并滤除散斑信号。本发明可以应用在激光雷达探测环境恶劣的情况下，无需多次重复测量或其他先验条件，即可对信噪比低的弱回波信号完成较好的处理效果，实现了较高的测距精度。

附图说明

图1是本发明的方法流程图；

图2是ica算法数学模型示意图；

图3是fast-ica算法的流程图；

图4是分离结果自适应判决流程图；

图5包括(a)、(b)、(c)三幅图，分别是激光雷达原始回波信号、加性去噪信号及参考脉冲信号的波形图；

图6包括(a)、(b)、(c)、(d)四幅图，分别是参考脉冲于900ns、1100ns位置时从原始信号中ica分离出的信号及噪声结果图；

图7包括(a)、(b)两幅图，分别是参考脉冲于1000ns位置时从原始信号中ica分离出的信号及噪声结果图。

具体实施方式

以下将结合附图，对本发明的技术方案进行详细说明。

一种激光雷达回波信号散斑噪声的滤除方法，如图1所示，具体步骤如下。

步骤1：对原始激光回波信号进行加性噪声去除，经过同态变换得到的结果作为独立成分分析的第一路观测信号。

原始回波信号可表示为e(t)＝m(t)s(t-t0)+n(t)，其中s(t)为发射激光脉冲信号，t0为激光实际飞行时间，m(t)为总的散斑噪声，包括大气湍流噪声和目标反射噪声，n(t)为总的加性噪声，包括大气衰减、背景光噪声、电流噪声等。首先需对加性噪声成分进行消减，当m(t)s(t-t0)＞＞n′(t)时，e(t)≈m(t)s(t-t0)，那么所得到的信号可以认为主要受到散斑噪声影响，经同态变换变后作为ica第一路观测信号，表示为：

x1(t)＝log(m(t)s(t-t0))＝logs(t-t0)+logm(t)

步骤2：对参考脉冲进行时域循环移n位，n初始值为1，将移位后的参考信号与原始激光回波信号进行相乘，经过同态变换所得到的结果作为独立成分分析的第二路观测信号。

利用激光发射脉冲信号s(t)作为参考信号，对其采样波形s(n)移位n位，假设采样率为δt，那么时移为τ＝nδt，得到的s(t-τ)与第一路观测信号进行相乘，经同态变换后作为ica第二路观测信号，表示为：

x2(t)＝log(s(t-τ)·x1(t))＝log(s(t-τ)·m(t)s(t-t0))＝logs(t-τ)+logs(t-t0)+logm(t)

步骤3：将第一路观测信号和第二路观测信号输入独立成分分析算法输出两路分离结果。

当τ≠t0时，s(t-τ)≠s(t-t0)，此时参考信号未匹配到实际回波信号位置，观测向量表示为当τ＝t0时，s(t-τ)＝s(t-t0)，此时参考信号匹配到实际回波信号位置，观测向量表示为其中，a为混合矩阵，s为要分离的源矩阵。可以看出，当τ≠t0，混合矩阵a不是满秩矩阵，因此现有的样本不足够分离三路信号。这种情况下只能分离两路信号，一路为logs(t-τ)，另一路为logm(t)s(t-t0)，这意味着参考脉冲未匹配时，ica分离的两路信号并没有成功将回波信号中的有效信号s(t-t0)分离出来实现去除散斑的目的。当τ＝t0，混合矩阵a是满秩的，因此，ica算法最终可以分离出两路数据，一路为信号项logs(t-t0)，另一路为噪声项logm(t)。

以上的观测向量数学模型推演过程分析了通过ica算法分离的不同结果来判断参考脉冲是否匹配实际回波位置的理论可行性。而具体的ica算法模型及原理示意图如图2所示。独立成分分析是指从几个源信号的混合信号中分离出源信号的算法。而对于源信号，只需保证它们之间互不相关，而不知道关于源信号的任何具体信息，因此也称为盲分离。ica的输入数据模型表示为x＝as，其中s为源信号矩阵，a为混合矩阵，x为接收到的观测信号矩阵。混合矩阵a和信号源s都是未知的，只需要通过观测信号矩阵x即可估计出源信号s。ica算法需要找出一个解混矩阵w，使输出y＝wx尽可能地与源信号s一致，实现源信号的最佳估计。实施例中，x＝[x1(t)x2(t)]^t，找到解混矩阵w后输出y＝wx＝[y1(t)y2(t)]^t。这时估计的两路结果可能对应的是logs(t-τ)和logm(t)s(t-t0)，也可能是logs(t-t0)和logm(t)，需要通过步骤4来判断。

图3是fast-ica算法的流程图。fast-ica算法是一种独立成分分析的快速算法，以非高斯性最大化为基础，采用固定点迭代来使得非高斯性达到最大，该算法通过牛顿迭代方式对接收信号x的很多样值点进行大量数据处理，逐渐从样本值中分离出各个独立成分。考虑到回波信号最终只需分离有效信号和散斑噪声两路，这里以两路观测信号为例，结合图3阐释fast-ica算法实现过程：

(1)中心化：对观测信号x＝[x1,x2]的三个独立分量去均值；

(2)白化：对去均值后的观测信号进行白化：

a)求的协方差矩阵；

b)对该协方差矩阵进行特征分解，得到2个特征值组成的对角阵e与三个特征向量组成的特征矩阵d；

c)求得白化矩阵：

d)得到白化的最终数据

e)验证白化的有效性：cov(v)＝e(vv^t)＝i(计算机上应为近似单位对角阵)，说明v的各个分量vi互不相关；

(3)选择一个具有单位范数的初始化向量w，进入迭代环节：

w←e{vg(w^tv)}-e{g′(w^tv)}w

这里的g函数一般选取：g1(y)＝tanh(a1y)，g2(y)＝yexp(-a2y²/2)或g3(y)＝y³，其中参数a1,a2∈[1,2]；

(4)标准化w，w←w/||w||；如未收敛，返回(3)继续迭代；如收敛，由y＝wx得到分离结果。

步骤4：对两路结果作同态逆变换，得到一路为分离信号，另一路为分离噪声。计算分离噪声的局部均方差，判断是否大于初始设定阈值，如果是，则判断本次迭代中ica分离的信号不符合要求，令n＝n+1；如不是，则判断本次迭代得到的回波信号符合要求，保存本次结果，令n＝n+1。

图4给出步骤(4)中对ica分离结果进行自适应判决的方法流程图。首先，对ica分离的数字信号结果做同态逆变换，即和分别称为分离信号和分离噪声。接着，对分离噪声找出其在时域范围内幅值在阈值标准以上的所有数据点y2′(n)，即满足为的平均值，α为预设参数，作为阈值标准。然后，计算y2′(n)的均方差，即的局部均方差：这里为y2′(n)长度，为y2′(n)平均值。与预设误差ε比较，如果σ＞ε，则判断当前迭代步参考脉冲未匹配实际回波脉冲，则舍弃本迭代步的结果，直接增加移位步长n＝n+1，进入步骤5；反之，则判断当前迭代步参考脉冲在误差ε范围内匹配到实际回波脉冲，将当前迭代步的ica分离结果保存，再增加移位步长n＝n+1，进入步骤5。

步骤5：判断移位步长n是否大于回波信号波形的时域长度，若不是，返回步骤2，继续迭代，若是，则停止迭代。

步骤6：所保存的各次分离噪声信号的局部均方差均符合阈值标准，选择局部均方差最小的那一次循环结果作为信号与散斑噪声分离的最佳结果。至此，完成了激光雷达回波信号的散斑噪声消减。

在图5中，(a)、(b)、(c)分别是原始回波信号、加性去噪信号及参考脉冲信号的波形图。在本实施例中，激光雷达的采样率为1ghz，因此波形图中相邻点距为1ns。采用小波阈值去噪算法对加性噪声进行消减，小波参数设定为层数8，小波基为sym8。去噪后利用五点三次平滑算法对毛刺进行消除，得到的结果如图5中的(b)所示。图5中的(c)中用于ica处理的参考脉冲信号即为激光雷达的发射信号，波长为532nm，脉宽为5ns。

图6展示了迭代算法过程中当n＝900和n＝1100迭代得到的结果，均包括分离信号和分离噪声两路结果。当n＝900时，采样间隔为1ns，参考脉冲移位τ＝900ns，经过ica分离得到的信号和噪声分别如图6中的(a)和(b)所示。从图6中的(b)的分离噪声波形上看，对应参考脉冲并未匹配到正确的回波脉冲信号位置。类似地，当n＝1100时，参考脉冲移位1100ns，经过ica分离得到的信号和噪声分别如图6中的(c)和(d)所示，依然可定性地判断此情况亦未匹配到正确的回波脉冲信号位置。注意这里图6中的(a)-(d)为ica算法输出后经过同态逆变换的结果。

图7展示了迭代算法过程中当n＝1000迭代得到的分离信号和分离噪声两路结果。由于采样间隔为1ns，参考脉冲移位τ＝1000ns，经过ica分离得到的信号和噪声分别如图7中的(a)和(b)所示。从图7中的(b)的分离噪声波形上看，第二路的噪声只保留了散斑噪声，而有效的回波脉冲信号均被分离到第一路中，印证了参考脉冲准确匹配的情况。

在迭代算法中，需要借助步骤4提出的局部均方差判决方法，通过数值分析而非人工目测来自动判断参考脉冲是否匹配，进而决定当前迭代步的处理结果是否需要保留。如图4所示，局部均方差为的平均值，α为预设参数，为y2′(n)长度，为y2′(n)平均值。实施例中，设定α＝0.1，计算得到900ns、1100ns和1000ns的σ分别为0.2534，0.2531和0.1049。预设误差ε＝0.1，那么只有1000ns符合标准，而900ns和1100ns因为超出标准证明分离噪声在时域范围内有尖峰，参考脉冲未匹配，从而舍弃这两组结果。需要说明的是，之所以要找到分离噪声在某个阈值标准以上的数据均方差，而非整个时域数据均方差的原因是，整个时域数据量很大，脉冲尖峰的振动偏差会被周边数据点振动偏差拉低，产生比无尖峰的情况更低的均方差，从而使得均方差越大来判决时域内产生尖峰的方法出现错误。实施例中，如果针对整个时域，900ns、1100ns和1000ns的均方差分别为0.0183，0.0182和0.0722。有尖峰的分离噪声均方差反而很低，以此判断整个时域内数据分布较为均匀并无尖峰将产生错误。因此，要缩小范围，取尖峰附近的数据点计算其均方差。范围则是通过前述的阈值标准来确定的。

在本实施例中，经过迭代后保存的结果均在误差ε范围内，寻找到最佳的数据位置为1000ns。因此，散斑去噪后的有效回波脉冲信号如图7中的图(a)所示，从而可以准确地确定激光脉冲的飞行时间为1000ns，对应的目标点距离为3×10⁸m/s×1000×10^-9s/2＝150m。

实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。