论文部分内容阅读
摘要： 为了从相干多普勒测风激光雷达的微弱气溶胶后向散射信号中高效、准确地提取多普勒频移信息，得到最接近CRB极限的无偏估计，利用气溶胶后向散射信号统计特性，建立了相干多普勒激光雷达仿真模型，并用最大概似法（ML）推导了在功率谱宽度已知情况下多普勒频率的估计算法。利用1.54 μm全光纤相干激光雷达系统获得风场测量数据，用ML算法进行多普勒频移估计，并将结果与传统的周期图最大值法（PM）所得结果进行比较。对比表明，ML算法与PM算法的数据处理结果有较好的一致性，并且ML算法的多普勒频率探测能力更优。
关键词： 最大概似法（ML）； 多普勒频移估计； 激光雷达； 周期图最大值法（PM）
中图分类号： TN 24文献标识码： Adoi： 10.3969/j.issn.10055630.2012.05.006
引言相干多普勒测风激光雷达能够在晴空下有效地获得高时间、空间分辨力和高精度的大气风速实时分布信息[1]，具有广泛的应用范围，包括三维风场测量、大气湍流探测、风切边预警和飞机尾迹涡流探测等[24]，是目前激光雷达的一个重要发展方向。相干多普勒测风激光雷达利用微弱的气溶胶后向散射信号，通过与本振光混频的方法提取风速引起的激光多普勒频移。在风速反演中，由于需要从含低信噪比的气溶胶后向散射信号中提取多普勒频移信息，高效、准确的算法对系统研制至关重要。目前较为普遍的多普勒频移反演算法主要有CRB[5]（cramerrao bound）、PP算法（pulse pair）、PPP算法（polypulsepair）、PM算法（periodgram maximum）和ML（maximum likelihood）等。CRB给出了方差最小的无偏估计，但是计算过于复杂；PP算法和PPP算法计算比较简单，但在低信噪比情况下，效率不高；ML（Maximum Likelihood）算法最接近CRB极限。文中主要对ML算法在实际系统中的应用进行研究，并对ML算法利用FFT适当降低计算量，同时对ML和PM算法的实际处理效果进行对比。ML算法通过确定最大概似函数，实现回波信号的特性参数。其中主要参数有功率谱宽度ω，信噪比SNR，以及多普勒频移fd（对应视向风速v）。如果要对这三个参数都进行估算，计算量非常大，因此根据实际系统情况，首先信号功率谱宽度ω，并通过仿真的方法，确定回波信号的信噪比SNR，通过数值仿真的方法，计算多普勒频移，从而大大降低了计算量，使得ML算法贴近实际应用。文中分别利用ML和PM算法，通过建立的仿真模型，以及实际系统采集到的回波信号进行数据反演处理，进行性能比较。1概率函数估计在大多数情况下，多普勒激光雷达信号可以描述为一个零均值圆型复数高斯随机过程z—[6]，采样时间为Ts，数组下标k的数组元素为zk，其表达式为zk=skexp（4πiυTsk/λ）+nk（1）式（1）中，skexp（4πiυTsk/λ）代表大气回波信号，i=-1，υ=λf/2为平均风速，λ为激光雷达波长，f为信号的平均频率，nk是非相关接收噪声。零均值圆型复数高斯随机过程可以由它的协方差矩阵R—（Θ—）完全描述[7]，Rkl （Θ—）=〈zkz′l〉（2）式（2）中，Θ—表示信号参量。R—（Θ—）的一般模型是Rn=SNRexp（4πiυTsn/λ-2π2ω2T2sn2）+δn（3）式（3）中，SNR=〈|sk|2〉/〈|nk|2〉，n=0时δn=1，其他情况下δn=0。光学仪器第34卷
第5期靳笑晗，等：最大概似法在多普勒激光测风雷达中的应用研究
计算中用到的大气回波信号参量Θ—包括υ，SNR和ω。令最大概似函数L（z—，Θ—）取最大值的Θ—即为ML估计值Θ⌒。对于零均值圆型复数高斯随机过程[8]，L（z—，Θ—）=-z—′R—（Θ—）z—-lnR—（Θ—）-Mlnπ（4）式（4）中，z—′代表z—的共轭转置，R—代表矩阵行列式，M是数组z—中的复数点个数。对于文中用到的信号模型，最大概似函数可写做L（z—，Θ—）=-z—′A—′Q-1（Φ）A—z—-lnQ—（Φ）-Mlnπ（5）式（5）中，A—是一个对角阵Akk=exp（4πiυTsk/λ），Q—（Φ）=R—（Θ—，υ=0），对于高斯模型矩阵中的参量Φ=（SNR，ω）。从估计理论可知，在相同条件下，弱信号的先验知识越少及所需估计的参数越多，则性能越差。当多普勒激光雷达发射具有高斯时域分布的信号，并且估计建立在对信号进行高频率采样的基础上时，ω可计算得到；对于地基激光雷达，通过计算多个连续的脉冲，可以得到SNR的有效估计，通过选择合适的大气模型，仿真得到SNR的有效值。ω和SNR确定后，通过计算概似函数取得最大值时对应的υ，即为ML估计得到的风速值，最大概似函数化简后[9]，L1（z—，（υ，SNR，ω））=-Re∑M-1m=0dmexp（-4πiTsυm/λ）（6）
dm=∑M-m-1k=0z′kzk+mDk，k+m（SNR，ω）（7）式（7）中，D—（SNR，ω）=Q—-1（SNR，ω）。估计的实际计算用到了快速傅里叶变换（FFT）算法。2ML和PM算法性能对比为了比较ML和PM算法，对两种算法进行了基于仿真模型和实际测量数据的对比反演。其中仿真模型采用了零平均高斯随机过程对气溶胶后向散射信号进行模拟，并结合激光雷达方程计算不同距离门的信噪比，用于替代ML算法在仿真数据和实验数据处理过程中一个未知参数SNR，从而减少ML算法的计算量。
2.1基于仿真模型的ML和PM算法对比气溶胶的浓度对系统的测量性能有很大影响，因此选择合适的大气模型，对仿真计算能否合理评估系统性能有着十分关键的作用。考虑到实际测量1.54 μm波长气溶胶后向散射系数比较困难，仿真过程中采用了1984年美国标准大气增强模式大气模型（1.64 μm），并进行相应的波长转换（1.64 μm到1.54 μm），实际测量环境是城市上空，为使计算与实验结果保持较好的一致性，将气溶胶后向散射系数乘2。图1为仿真采用的大气模型及仿真获得的气溶胶后向散射信号时域波形。 采用上述仿真模型获得在2.025 km距离门处的气溶胶后向散射信号，设定中频fd=5.614 MHz、风速υ=0（系统中频为56.15 MHz，对应多普勒频移为0，风速为0），脉冲功率为18 μJ。计算得到该距离门内SNR=6.77×10-4。为了有效实现弱信号的探测，分别在PM和ML计算过程中，对脉冲功率谱和最大概似函数进行累加，N代表脉冲累加次数。图2是在不同脉冲累加次数下ML算法得到的概似函数分布，图3是不同脉冲累加次数下PM算法得到的功率谱的分布。对比图2和图3可以看出，在脉冲累加次数为3 000次时，两种算法的频谱图都存在明显峰值，能够正确地提取多普勒频率；当脉冲累加次数减少到2 000、1 000次时，PM算法的功率谱图显示有用信号已经被噪声淹没，无法提取正确的多普勒频移，然而此时ML算法的频谱图中主峰值位置不变且明显，依然能够有效提取得到正确的结果。
2.2基于系统实测数据的对比除了在通过建立仿真模型的方法对两种算法进行对比，利用激光雷达对实际大气情况下测量的数据进行数据反演进行对比，对系统的实际应用有着很重要的评估作用。因此实验中，采用了实验室研制的一套1.54 μm全光纤相干激光雷达系统进行风速测量实验。该系统的参数如表1所示，分别进行垂直和水平探测，同时为了更好地对两种算法的探测概率进行对比，实验中对硬目标进行单脉冲探测，并对探测结果进行统计分析。根据2.1中对比结果，对回波信号进行3 000个脉冲的累加，然后提取多普勒频移值，反演视线风速。
名称参数波长（nm）1 539.7脉冲能量（μJ）50脉冲宽度（ns）500脉冲重复频率（kHz）10本振光功率（nW）260空间分辨力（m）75望远镜孔径（mm）50/100脉冲累加次数3 000系统聚焦范围（km）2
图4为垂直和水平探测时，用两种算法反演的视线风速分布图，其中实线代表ML估计结果，虚线代表PM估计结果。由图4（a）可知，垂直测量时，前26个距离门，即1.95 km范围内，ML估计与PM估计所得结果存在很好的一致性，观察频谱可知，3 000次脉冲累加时，26个距离门之后的有用信号已被噪声淹没，所以超过1.95 km反演得到的视线风速是错误的，可以忽略。由图4（b）可知，35个距离门以内ML估计与PM估计所得结果有很好的吻合。水平测量的有效范围拓展到了2.625 km，这是因为气溶胶的浓度稳定在一个较大的值，同样的测量距离时，水平探测的后向散射信号强度相比垂直探测的更大。
echo signal from hard target2.3单脉冲光硬目标回波数据的ML估计2.625 km处一栋建筑屋顶，图5为单脉冲在时域上的波形图，第一个脉冲是光纤端面的反射信号，第二个脉冲为硬目标的反射信号（16.8 μs）。分别用这两算法提取单脉冲硬目标反射信号的多普勒频移，并将1 500个脉冲得到的结果进行统计，比较两者性能。获得1 500组光脉冲回波数据。计算过程中，估算目标回波信号的信噪比为SNR=2.65，用ML算法和PM算法处理硬目标的回波数据（即第35个距离门内的采样数据），反演得到在[-1.5 m/s，1.3 m/s]范围内的视线风速概率分布如图6所示。
硬目标相对激光雷达是静止的，理论上多普勒频移为零。由图4可知，两种算法得到的视线风速值均为υ^=-0.378 2 m/s，在误差范围内。ML估计得到的视线风速在υ^处的概率为PML（υ^）=0.606，PM在υ^处的概率为PPM（υ^）=0.485，两种估计的方差分别为σML=0.052 6，σPM=0.067 5。可知ML估计的结果更集中，误差更小。3结论提出一种最接近CRB极限的ML算法，它能够克服其他算法在低信噪比情况下，效率不高的弊端，通过确定最大概似函数，实现回波信号的特性参数。利用该算法分别对仿真模型数据和系统实测数据进行多普勒频移反演，同时对ML和PM算法的实际处理效果进行对比。对仿真模型数据的处理结果显示，当系统信噪比SNR及相干信号功率谱宽ω先验已知，并设定正确时，ML估计具有更强的去噪能力，并不要求大量的脉冲累加。对系统实测数据的处理结果表明ML算法相比传统的PM算法，是一种更有效的多普勒频移反演工具。以上结果证明该算法在工程应用领域有很好的应用前景，特别是在一些不容易获得试验数据或者试验数据获得代价高的领域。参考文献：
[1]ASAKA K，YANAGISAWA T，HIRANO Y.1.5 μm eyesafe coherent lidar system for wind velocity measurement[J].SPIE，2001，4153：321-328.
[2]KOCH G J，BEYON J Y，PETZAR P，et al.Field testing of a highenergy 2 μm Doppler lidar[J].Journal of Applied Remote Sensing，2010，4：043512-043513.
[3]DOLFIBOUTEYRE A，CANAT G，VALLA M，et al.Pulsed 1.5 μm lidar for axial aircraft wake vortex detection based on highbrightness largecore fiber amplifier[J].IEEE，2009，15：441-450.
[4]INOKUCHI H，ENDO E，ANDO T，et al.Development of an airborne wind measurement system[J].SPIE，2009，7382：738205-738212.
[5]TIMOTHY J S，ALEXANDER B K.Variance boums on the estimation of reflectivity and polarizitio parameters in radar meterology[J].IEEE，1997，35：248-255.
关键词： 最大概似法（ML）； 多普勒频移估计； 激光雷达； 周期图最大值法（PM）
中图分类号： TN 24文献标识码： Adoi： 10.3969/j.issn.10055630.2012.05.006
引言相干多普勒测风激光雷达能够在晴空下有效地获得高时间、空间分辨力和高精度的大气风速实时分布信息[1]，具有广泛的应用范围，包括三维风场测量、大气湍流探测、风切边预警和飞机尾迹涡流探测等[24]，是目前激光雷达的一个重要发展方向。相干多普勒测风激光雷达利用微弱的气溶胶后向散射信号，通过与本振光混频的方法提取风速引起的激光多普勒频移。在风速反演中，由于需要从含低信噪比的气溶胶后向散射信号中提取多普勒频移信息，高效、准确的算法对系统研制至关重要。目前较为普遍的多普勒频移反演算法主要有CRB[5]（cramerrao bound）、PP算法（pulse pair）、PPP算法（polypulsepair）、PM算法（periodgram maximum）和ML（maximum likelihood）等。CRB给出了方差最小的无偏估计，但是计算过于复杂；PP算法和PPP算法计算比较简单，但在低信噪比情况下，效率不高；ML（Maximum Likelihood）算法最接近CRB极限。文中主要对ML算法在实际系统中的应用进行研究，并对ML算法利用FFT适当降低计算量，同时对ML和PM算法的实际处理效果进行对比。ML算法通过确定最大概似函数，实现回波信号的特性参数。其中主要参数有功率谱宽度ω，信噪比SNR，以及多普勒频移fd（对应视向风速v）。如果要对这三个参数都进行估算，计算量非常大，因此根据实际系统情况，首先信号功率谱宽度ω，并通过仿真的方法，确定回波信号的信噪比SNR，通过数值仿真的方法，计算多普勒频移，从而大大降低了计算量，使得ML算法贴近实际应用。文中分别利用ML和PM算法，通过建立的仿真模型，以及实际系统采集到的回波信号进行数据反演处理，进行性能比较。1概率函数估计在大多数情况下，多普勒激光雷达信号可以描述为一个零均值圆型复数高斯随机过程z—[6]，采样时间为Ts，数组下标k的数组元素为zk，其表达式为zk=skexp（4πiυTsk/λ）+nk（1）式（1）中，skexp（4πiυTsk/λ）代表大气回波信号，i=-1，υ=λf/2为平均风速，λ为激光雷达波长，f为信号的平均频率，nk是非相关接收噪声。零均值圆型复数高斯随机过程可以由它的协方差矩阵R—（Θ—）完全描述[7]，Rkl （Θ—）=〈zkz′l〉（2）式（2）中，Θ—表示信号参量。R—（Θ—）的一般模型是Rn=SNRexp（4πiυTsn/λ-2π2ω2T2sn2）+δn（3）式（3）中，SNR=〈|sk|2〉/〈|nk|2〉，n=0时δn=1，其他情况下δn=0。光学仪器第34卷
第5期靳笑晗，等：最大概似法在多普勒激光测风雷达中的应用研究
计算中用到的大气回波信号参量Θ—包括υ，SNR和ω。令最大概似函数L（z—，Θ—）取最大值的Θ—即为ML估计值Θ⌒。对于零均值圆型复数高斯随机过程[8]，L（z—，Θ—）=-z—′R—（Θ—）z—-lnR—（Θ—）-Mlnπ（4）式（4）中，z—′代表z—的共轭转置，R—代表矩阵行列式，M是数组z—中的复数点个数。对于文中用到的信号模型，最大概似函数可写做L（z—，Θ—）=-z—′A—′Q-1（Φ）A—z—-lnQ—（Φ）-Mlnπ（5）式（5）中，A—是一个对角阵Akk=exp（4πiυTsk/λ），Q—（Φ）=R—（Θ—，υ=0），对于高斯模型矩阵中的参量Φ=（SNR，ω）。从估计理论可知，在相同条件下，弱信号的先验知识越少及所需估计的参数越多，则性能越差。当多普勒激光雷达发射具有高斯时域分布的信号，并且估计建立在对信号进行高频率采样的基础上时，ω可计算得到；对于地基激光雷达，通过计算多个连续的脉冲，可以得到SNR的有效估计，通过选择合适的大气模型，仿真得到SNR的有效值。ω和SNR确定后，通过计算概似函数取得最大值时对应的υ，即为ML估计得到的风速值，最大概似函数化简后[9]，L1（z—，（υ，SNR，ω））=-Re∑M-1m=0dmexp（-4πiTsυm/λ）（6）
dm=∑M-m-1k=0z′kzk+mDk，k+m（SNR，ω）（7）式（7）中，D—（SNR，ω）=Q—-1（SNR，ω）。估计的实际计算用到了快速傅里叶变换（FFT）算法。2ML和PM算法性能对比为了比较ML和PM算法，对两种算法进行了基于仿真模型和实际测量数据的对比反演。其中仿真模型采用了零平均高斯随机过程对气溶胶后向散射信号进行模拟，并结合激光雷达方程计算不同距离门的信噪比，用于替代ML算法在仿真数据和实验数据处理过程中一个未知参数SNR，从而减少ML算法的计算量。
2.1基于仿真模型的ML和PM算法对比气溶胶的浓度对系统的测量性能有很大影响，因此选择合适的大气模型，对仿真计算能否合理评估系统性能有着十分关键的作用。考虑到实际测量1.54 μm波长气溶胶后向散射系数比较困难，仿真过程中采用了1984年美国标准大气增强模式大气模型（1.64 μm），并进行相应的波长转换（1.64 μm到1.54 μm），实际测量环境是城市上空，为使计算与实验结果保持较好的一致性，将气溶胶后向散射系数乘2。图1为仿真采用的大气模型及仿真获得的气溶胶后向散射信号时域波形。 采用上述仿真模型获得在2.025 km距离门处的气溶胶后向散射信号，设定中频fd=5.614 MHz、风速υ=0（系统中频为56.15 MHz，对应多普勒频移为0，风速为0），脉冲功率为18 μJ。计算得到该距离门内SNR=6.77×10-4。为了有效实现弱信号的探测，分别在PM和ML计算过程中，对脉冲功率谱和最大概似函数进行累加，N代表脉冲累加次数。图2是在不同脉冲累加次数下ML算法得到的概似函数分布，图3是不同脉冲累加次数下PM算法得到的功率谱的分布。对比图2和图3可以看出，在脉冲累加次数为3 000次时，两种算法的频谱图都存在明显峰值，能够正确地提取多普勒频率；当脉冲累加次数减少到2 000、1 000次时，PM算法的功率谱图显示有用信号已经被噪声淹没，无法提取正确的多普勒频移，然而此时ML算法的频谱图中主峰值位置不变且明显，依然能够有效提取得到正确的结果。
2.2基于系统实测数据的对比除了在通过建立仿真模型的方法对两种算法进行对比，利用激光雷达对实际大气情况下测量的数据进行数据反演进行对比，对系统的实际应用有着很重要的评估作用。因此实验中，采用了实验室研制的一套1.54 μm全光纤相干激光雷达系统进行风速测量实验。该系统的参数如表1所示，分别进行垂直和水平探测，同时为了更好地对两种算法的探测概率进行对比，实验中对硬目标进行单脉冲探测，并对探测结果进行统计分析。根据2.1中对比结果，对回波信号进行3 000个脉冲的累加，然后提取多普勒频移值，反演视线风速。
名称参数波长（nm）1 539.7脉冲能量（μJ）50脉冲宽度（ns）500脉冲重复频率（kHz）10本振光功率（nW）260空间分辨力（m）75望远镜孔径（mm）50/100脉冲累加次数3 000系统聚焦范围（km）2
图4为垂直和水平探测时，用两种算法反演的视线风速分布图，其中实线代表ML估计结果，虚线代表PM估计结果。由图4（a）可知，垂直测量时，前26个距离门，即1.95 km范围内，ML估计与PM估计所得结果存在很好的一致性，观察频谱可知，3 000次脉冲累加时，26个距离门之后的有用信号已被噪声淹没，所以超过1.95 km反演得到的视线风速是错误的，可以忽略。由图4（b）可知，35个距离门以内ML估计与PM估计所得结果有很好的吻合。水平测量的有效范围拓展到了2.625 km，这是因为气溶胶的浓度稳定在一个较大的值，同样的测量距离时，水平探测的后向散射信号强度相比垂直探测的更大。
echo signal from hard target2.3单脉冲光硬目标回波数据的ML估计2.625 km处一栋建筑屋顶，图5为单脉冲在时域上的波形图，第一个脉冲是光纤端面的反射信号，第二个脉冲为硬目标的反射信号（16.8 μs）。分别用这两算法提取单脉冲硬目标反射信号的多普勒频移，并将1 500个脉冲得到的结果进行统计，比较两者性能。获得1 500组光脉冲回波数据。计算过程中，估算目标回波信号的信噪比为SNR=2.65，用ML算法和PM算法处理硬目标的回波数据（即第35个距离门内的采样数据），反演得到在[-1.5 m/s，1.3 m/s]范围内的视线风速概率分布如图6所示。
硬目标相对激光雷达是静止的，理论上多普勒频移为零。由图4可知，两种算法得到的视线风速值均为υ^=-0.378 2 m/s，在误差范围内。ML估计得到的视线风速在υ^处的概率为PML（υ^）=0.606，PM在υ^处的概率为PPM（υ^）=0.485，两种估计的方差分别为σML=0.052 6，σPM=0.067 5。可知ML估计的结果更集中，误差更小。3结论提出一种最接近CRB极限的ML算法，它能够克服其他算法在低信噪比情况下，效率不高的弊端，通过确定最大概似函数，实现回波信号的特性参数。利用该算法分别对仿真模型数据和系统实测数据进行多普勒频移反演，同时对ML和PM算法的实际处理效果进行对比。对仿真模型数据的处理结果显示，当系统信噪比SNR及相干信号功率谱宽ω先验已知，并设定正确时，ML估计具有更强的去噪能力，并不要求大量的脉冲累加。对系统实测数据的处理结果表明ML算法相比传统的PM算法，是一种更有效的多普勒频移反演工具。以上结果证明该算法在工程应用领域有很好的应用前景，特别是在一些不容易获得试验数据或者试验数据获得代价高的领域。参考文献：
[1]ASAKA K，YANAGISAWA T，HIRANO Y.1.5 μm eyesafe coherent lidar system for wind velocity measurement[J].SPIE，2001，4153：321-328.
[2]KOCH G J，BEYON J Y，PETZAR P，et al.Field testing of a highenergy 2 μm Doppler lidar[J].Journal of Applied Remote Sensing，2010，4：043512-043513.
[3]DOLFIBOUTEYRE A，CANAT G，VALLA M，et al.Pulsed 1.5 μm lidar for axial aircraft wake vortex detection based on highbrightness largecore fiber amplifier[J].IEEE，2009，15：441-450.
[4]INOKUCHI H，ENDO E，ANDO T，et al.Development of an airborne wind measurement system[J].SPIE，2009，7382：738205-738212.
[5]TIMOTHY J S，ALEXANDER B K.Variance boums on the estimation of reflectivity and polarizitio parameters in radar meterology[J].IEEE，1997，35：248-255.