基于最大似然估计的相干多普勒雷达信号处理

邢晓晴，陆 威，庄子波

（1.中国民用航空飞行学院航空工程学院，四川 广汉 618307；2.哈尔滨工业大学可调谐（气体）激光技术重点实验室，哈尔滨 150001；3.中国民航大学飞行技术学院，天津 300300）

大气风场的快速精准测量一直是世界各国研究的热点，具有重大社会意义与较高的实用价值。一方面，大气风场可以表征大气中物质交换的动向，用于气候预测[1]。另一方面，大气风场的实时精准探测可有效保障飞机、火箭等飞行器的起降与发射安全。相比于微波雷达，激光雷达具有灵敏度较高的优势，并且是目前国际唯一公认的具备晴空测风能力的设备，具有广阔的发展前景。相干激光雷达利用回波信号与发射信号的多普勒频移差获取风场信息，在低空风场的探测领域具有明显优势。

相干激光雷达信号处理中，对回波信号的多普勒信息提取直接影响系统的测速精度，是系统设计的重点。目前常见的多普勒信息提取算法有脉冲对技术[2]、最大似然估计[3]、周期图最大似然估计和最大似然离散谱峰值估计（ML DSP）[4]等。上述算法中，脉冲对估计算法计算量较小，适用于高信噪比信号。Frehlich等[5]通过仿真对最大似然估计算法应用于微弱信号的情况进行了分析。ML DSP的核心是将信号频谱峰值对应的中心频率视为多普勒频率的最大似然解。最大似然离散谱峰值估计的是真实风速的概率，可选取概率密度函数对估计性能进行评价。概率密度函数由真实风速的高斯分布（“好”的估计）与均匀分布的噪声（“坏”的估计）组成。最大似然离散谱峰值估计是目前世界上广泛应用的风速反演算法。因此，通过Matlab编程与Simulink仿真工具箱的结合实现了最大似然离散谱峰值估计算法的应用，并将该方法成功应用于真实雷达信号的频率信息解算。

1 相干激光雷达结构

相干激光测风雷达系统由激光发射系统、信号探测系统、收发天线与信号处理系统组成，如图1所示。

图1 相干激光雷达系统结构Fig.1 Structure of coherent Doppler lidar

激光发射系统的主要功能是产生主振激光脉冲，为整个系统提供探测信号。激光器种子光采用线宽为3 kHz的单频窄线宽光纤激光器，由分束器分为两路。一路为本振光，另一路经过声光调制和放大器的作用，产生高能量的脉冲激光用于风场探测，即主振光。系统主振光与本振光频率差为80 MHz，可满足日常风场探测需求。主振光单脉冲能量为110 μJ，脉冲重复频率为10 kHz，波长为1 550 nm，处于人眼安全波段。

信号探测系统的主要功能是将系统收发信号隔离，对载有风场信息的回波信号进行探测。由于主振光单脉冲能量较高，会对光纤光环行器造成损害，影响隔离效率，因此该装置的偏振隔离系统采用空间光隔离结构。空间光隔离结构由1个偏振片与1个1/4波片组成。发射光依次经过偏振片与1/4波片。发射光的偏振方向与偏振片的透射方向一致，可全部通过偏振片。经过1/4波片后，发射光偏振状态为圆偏振，在探测目标处发生散射，形成回波信号。圆偏振光对各种偏振状态的探测目标均可实现有效探测。回波信号再依次经过1/4波片和偏振片。由于回波信号与发射信号偏振状态相同，因此，当回波信号经过1/4波片后，回波信号变为线偏振光，偏振方向与出射信号偏振方向垂直，故在偏振片处完全被反射进入合束器。合束器处实现回波信号与本振光的混频，由双平衡探测器探测出载有风场信息的中频信号。

信号处理系统的主要功能是对中频信号进行采集与处理，还原风场信息。中频信号经中频放大器进一步放大后，由示波器进行采集。系统采用安捷伦示波器进行数据采集，采样频率高达10 GHz。数据处理部分采用Matlab编程与Simulink仿真工具箱结合的方式，依据信号具体特征选取最佳信号处理方法，进一步提高系统的探测精度。

2 风速反演理论

2.1 估计风速的概率密度函数

相干激光测风雷达接收的信号为气溶胶粒子的后向散射信号。后向信号较为微弱，其理论极限由Cramer-Rao下界表征[6]。另一方面，估计的风速仅在一定概率上可代表真实风速，一般采用概率密度函数（PDF，probability density function）进行描述。概率密度函数由两部分组成，一部分是以风场引起的多普勒频移为中心的高斯分布，被称为“好”的估计；另一部分是由噪声引起的，通常为平均分布的白噪声，被称为“坏”的分布。当信号信噪比较低时，多普勒信号有可能淹没在噪声当中。可将风速的估计PDF表示为

其中：v代表实际风速；g为多普勒频移估计的标准差；b 为噪声估计的概率；b0为“坏”的估计，在（vmin，vmax）内均匀分布。

2.2 最大似然离散谱估计（ML DSP）

相干激光测风雷达中，风场反演的重要参数为回波信号频谱的平均频率。根据多普勒定律，该频率正比于散射物体的平均速度。大多数散射目标的谱线接近高斯分布，因此当离散谱分辨率较高时，可将频率分辨谱线中能量最大的谱线对应的频率作为中心频率的估计。当被估计频率的期望值等于频率分辨单元中心频率时，该估计为无偏估计。当被估计频率期望值不等于频率分辨单元中心频率时，该估计为有偏估计，偏差的最大值为频谱分辨宽度的50%。为更精确地估计频谱信息，可大致选定信号区间后，对其进行高斯拟合以提高频谱估计的准确度。

3 雷达信号处理方法

雷达信号处理的目的是提取回波信号中的多普勒信息[7-9]，通过多普勒信息计算对应风速。雷达回波信号处理的具体流程如图2所示。

图2 信号处理流程Fig.2 Signal processing procedure

根据多普勒信息反演风速，不同探测距离上的风场信息不同，为了测量在探测方向上的风场变化情况，首先要对回波进行多个距离门划分，距离门宽度通常由激光脉冲宽度决定。选取多组待测区域对应的距离门信号，对同一距离门内的信号进行时频转换，得到回波信号的离散频谱，通过非相干积累的方法对信号进行非相干积累后，通过频谱估计算法求解出对应风速。

功率谱的估计方法很多，其方法的选取决定了多普勒信息还原的准确性。采用最大似然功率谱估计方法，将多组观测值的联合概率密度分布视为待估函数，默认功率谱的概率密度函数为高斯分布，求出观测值中出现概率最大的参考值。最后通过多普勒频移公式对风速进行求解，获得风速信息。

4 信号获取与处理

频率测量精度直接影响系统的测速精度，是判定信号处理方法优劣的最重要参数之一。为验证最大似然估计法的功率谱估计性能，设计了一套检测系统测速精度的装置，并开展了一系列实验。

速度检测装置采用旋转电机搭载4片对称的叶片结构，如图3所示，该结构可为雷达提供固定运动速度的探测目标。通过设置电机不同的转向与转速，使探测点目标与系统分别进行相向、相离运动，探测点位置的速度范围为4.6～23.12 m/s，速度增减梯度为±1.85 m/s。

图3 速度检测装置实物图Fig.3 Velocity measurement device

雷达回波信号不仅包括探测目标的散射信号，还包括系统内部光学镜头的镜面反射信号以及白噪声。激光雷达探测目标的散射截面通常较小，散射信号一般较为微弱。内部光学镜头反射的信号距离发射端较近，反射的信号功率明显强于目标散射信号。适度增大目标与系统的距离[10]，使二者可在时间上进行有效区分。

雷达测速实验选取在室内进行。发射信号为波长1.55 μm的脉冲激光，激光脉冲宽度为300 ns，重复频率10 kHz。激光脉冲发射初期，系统内部的反射信号功率较高，将对回波信号的探测产生影响。因此实验将在光路中放置反射镜，增大电机反射信号光程，使其与镜头回波在时间上可区分。实验中，反射镜与雷达系统和速度检测装置的距离均为45 m。因此，电机与雷达系统的距离约为90 m。实验装置的相对位置关系如图4所示。

图4 实验装置结构示意图Fig.4 Experimental facilities

采用双平衡探测器对回波信号进行探测后，由安捷伦示波器进行高速采样，得到多组离散的雷达回波信号。选取1组雷达回波信号进行分析，如图5所示，回波信号可分为3段，a段为镜头回波信号，镜头回波信号较强，经中频放大器后探测器达到饱和，可观察到信号在峰值处有明显的削平现象。镜头回波的出现时刻可认定为激光发出的起始时刻。根据镜头前沿位置与激光测距方法，结合实验结构，判定b段主要为反射镜回波信号，c段为电机回波信号位置，其他位置为噪声。

图5 雷达回波信号Fig.5 CDL echo signal

对回波信号进行距离门划分，系统采用的激光脉冲宽度为300 ns，因此将距离门宽度选定为300 ns，对应系统的距离分辨率为45 m。选取电机回波位置（c段）的距离门进行快速傅里叶变换（FFT）得到回波信号的傅里叶频谱，如图6所示。由于进行时频转换的信号是有限长的，因此得到的傅里叶频谱为离散傅里叶频谱。信号相干时间为300 ns，对应频谱分辨率为3.3 MHz[11-12]。在此离散频谱下，系统对应的测速精度为2.54 m/s。采用ML DSP对信号进行频谱估计，可大幅提高系统的测速精度。

图6 信号频谱图Fig.6 Signal frequency spectrum

电机回波信号较强，对应的回波信号信噪比较高，认为噪声并未淹没信号。首先，以频谱幅值最高处为中心，结合信号的频谱宽度截取待拟合区域。通过Matlab工具箱对区域内谱线进行高斯拟合，如图7所示。图7中横轴代表信号频率，纵轴代表该频率下对应的频谱幅值。离散点代表对信号进行时频转换后的离散频谱，曲线代表对离散频谱进行高斯拟合后的结果。由频谱泄露的相关理论计算，离散频谱呈以回波信号频率为中心的高斯分布，当信号信噪比较高时，噪声引入的有偏估计较小。故通过Matlab实现对频谱的高斯拟合，拟合后的中心位置即为该估计方法对频谱的还原结果。最后，根据多普勒效应公式，计算得出风速估计值。

图7 信号的最大似然离散谱估计示意图Fig.7 Maximum likelihood discretization estimation of echo signal

通过16组不同速度运动目标的测速实验，采用上述最大似然估计法分别对信号进行频谱估计，计算出系统对运动目标的测量值。通过对频谱进行高斯拟合得到回波信号频率测量值，并根据频率计算速度，将雷达测速结果与电机叶片实际转速对比，如表1所示。

由表1对测量速度与目标实际速度的误差分析发现，本次开展的16组实验中系统实际测速误差均小于0.5 m/s。结合拟合相关系数对测速误差进行分析，结果显示，当拟合的相关系数较高时，系统测速误差相对较小。

表1 电机运动参数与测速结果对照表Tab.1 Motor parameters and measured velocities

对测速误差进行均方根计算得到系统测速误差（均方根）为0.26 m/s。对实际速度与测量速度结果进行直线拟合，结果如图8所示，所得直线斜率为0.998 6，可决系数为0.999 7，具有良好的线性和较高的精度。

图8 真实速度与实测速度拟合曲线Fig.8 Fitting curve of measured and real velocities

5 结语

最大似然离散谱估计作为一种重要的风速估计方法，目前在世界相干激光雷达领域得到广泛应用。主要通过编程与仿真工具箱的结合，实现了雷达信号频谱的最大似然估计。开展了一系列激光雷达速度检测实验，通过最大似然离散谱估计对信号进行处理，测速误差（均方根）小于0.3 m/s。对测量结果与实际速度进行直线拟合，斜率为0.998 6，可决系数为0.999 7，结果具有良好的线性。证明了最大似然离散谱估计算法的准确性，对相干激光雷达信号处理方法的研究具有借鉴意义。

参考文献：

[1]DIAO WEIFENG,ZHANG XIN,LIU JIQIAO,et al.All fiber pulsed coherent lidar development for wind profiles measurements in boundary layers[J].Chinese Optics Letters,2014,12(7)：71-74.

[2]WANG HENG,ZENG HAIYING,LI MIN,et al.Maximum likelihood estimationofclockskewinwirelesssensornetworkswithperiodicalclock correctionunderexponentialdelays[J].IEEETransactionsonSignal Processing,2017,65(10)：2714-2724.

[3]靳笑晗,汪岳峰,竹孝鹏,等.最大概似法在多普勒激光测风雷达中的应用研究[J].光学仪器,2012,34(5)：28-33.

[4]JIN XIAOHAN,WANG YEFENG,ZHU XIAOPENG,et al.Performance of maximum likelihood estimators of Doppler frequency for coherent Doppler lidar[J].Optical Instruments,2012,34(5)：28-33.

[5]FREHLICH R G,YADLOWSKY M J.Performance of mean frequency estimators for Doppler radar and lidar[J].Journal of Atmospheric&Oceanic Technology,1994,11(5)：1217.

[6]贾晓东.1.55 μm相干测风激光雷达样机的研制[D].合肥：中国科学技术大学,2015.

[7]COSTA E,PAULA E R D,REZENDE L F C,et al.Equatorial scintillation calculations based on coherent scatter radar and C/NOFS data[J].Radio Science,2016,46(2)：1-21.

[8]竹孝鹏,刘继桥,刁伟峰,等.相干多普勒测风激光雷达研究[J].红外,2012,33(2)：8-12.

[9]MATHEWS J D.Incoherent scatter radar probing of the 60-100-km atmosphere and ionosphere[J].IEEE Transactions on Geoscience&Remote Sensing,1986,GE-24(5)：765-776.

[10]HYSELL D L,WOODMAN R F.Imaging coherent backscatter radar observations of topside equatorial spread F[J].Radio Science,2016,32(6)：2309-2320.

[11]MEYER M G,SAHR J D.Passive coherent scatter radar interferometer implementation,observations,and analysis[J].Radio Science,2016,39(3)：1-10.

[12]邢晓晴.相干激光测风雷达大气探测方法实验研究[D].哈尔滨：哈尔滨工业大学,2016.