-
近年来，激光雷达作为重要的感知技术得到了不断的发展，被广泛应用于焊工安全、大气科学、自动驾驶、气象预报等重要领域[1-2]。然而，随着技术的发展，商用最为广泛的脉冲激光雷达的局限性也逐步展现出来，其在远距离探测、人眼安全、抗干扰能力、同时测距测速等方面的问题显得尤为突出[3-4]。研究人员开始将目光转向调频连续波 (Frequency Modulated Continuous Wave, FMCW) 激光雷达。FMCW激光雷达采用相干探测原理，增强了抗干扰能力并大大提升了测距范围，同时激光发射为连续波形式，不仅平均功率小、人眼安全性高，而且具备同时测距测速能力[5-6]。
根据FMCW激光雷达原理，雷达发射信号与回波信号混频后输出中频信号的频谱是测量目标距离和速度的关键，且直接影响到激光雷达测距测速精度和实时性等重要技术指标[7]。中频信号频谱获取采用FFT实现，但FFT存在固有的栅栏效应会导致频率计算的精度受到限制[8]。
目前提升FFT精度通常有两种方式，一种为增加FFT运算点数，这种方式比较简单，但是随着点数增加会导致FFT运算延迟迅速增大，在实际应用中容易引起实时性降低的问题；另一种为频率估计方式，频率估计不需要增加FFT运算点数，在提升精度的同时保证了实时性要求。因此，不少学者针对频率估计做了大量研究，如能量重心法、Rife算法、相位差法等。其中，Rife算法由于具有原理简单、易于工程实现等优点，获得更多关注。但是，在大噪声环境下，Rife算法的估计频率在靠近FFT量化频率点时，会出现修正方向判断错误的现象，实际测量时精度会有明显下降，无法满足高精度测量需求。为解决该问题，参考文献[9]提出I-Rife算法，其通过计算两个量化频率点的中间频谱模值确定修正方向，估计精度较高，但需要运算多个N点DFT，计算量较大。参考文献[10]采用基于相角和频谱模值相结合的方式，通过求解量化频率点与两侧相角差，并结合其模值大小进行修正方向判断，虽然计算量有所下降，但是相角计算时需要求解反正切函数，不利于硬件实现。参考文献[11]提出一种双门限判决修正Rife算法，其结合了I-Rife、Rife算法的优点，虽然精度有所提升，但需要进行频谱搬移，比I-Rife算法的计算量有所降低，但是计算方式仍较为复杂，同样不宜于硬件实现。
为同时满足FMCW激光雷达中频信号频率高精度和实时获取等方面的需求，文中提出一种修正Rife算法，该算法从易于硬件实现角度出发，基于相角判据完成频率估计，在保证精度的同时，计算量较小，具有较好的工程应用价值。
-
根据实际需求，FMCW激光雷达的光源需要对光载波频率进行不同形式的调制，常见方式有锯齿波、三角波和梯形波。为实现对目标同时测距测速，文中采用梯形波调制，其基本原理如图1(a) 所示。激光发射信号频率从f0开始随时间线性增加，线性调制时间为T，扫频带宽为B，其中线性上扫频、中间扫频和下扫频调制时间相等，对应梯形波调制周期为3T。由于激光器与目标之间存在距离，激光回波信号相对于发射信号会存在时间τ延迟。而当激光器与目标存在相对运动时，激光回波信号还会产生多普勒频移，使得回波激光信号与发射信号存在频域方向的平移
dfD。由于激光信号随时间线性变化使激光回波信号与激光发射信号产生了如图1(b) 所示的频率差 df，频率差由距离延迟τ和多普勒频率dfD构成，结合三角形相似和相干解调原理，目标的距离L和相对速度V可由df1和df3得出[12]： 
图 1 (a) 梯形波形式调制下激光发射信号和激光回波信号时频域关系图；(b) 经相干检测所得差拍中频信号时频域关系图；(c) 差拍中频电信号的时域波形
Figure 1. (a) Time frequency domain diagram of transmitted signal and received echo signal under the trapezoidal modulation; (b) Time frequency domain diagram of the beat signal after coherent detection; (c) Time domain waveform of the received beat signal
$$ L{\text{ = }}\dfrac{{c \cdot \tau }}{2} = \dfrac{{c \cdot ({\rm d}{f_1} + {\rm d}{f_3}) \cdot T}}{{4B}} $$ (1) $$ V{\text{ = }}\dfrac{{\left| {{\rm d}{f_1} - {\rm d}{f_3}} \right|}}{4} \cdot \lambda $$ (2) 式中：τ为激光雷达发射波和回波之间的时间间隔；c为光速；B为调制信号扫频带宽；T为调制信号线性扫频时间；λ为激光波长。
激光雷达中频信号经光电探测器和模拟前端电路转换后产生如图1(c) 所示的中频电信号，可以看出在延迟τ之后中频信号呈现较低频率正弦波形式，通过对每一段电信号进行傅里叶变换，可得到用于目标信息解算的频率值。
-
由FMCW原理可知，在一个线性扫频时间T内，激光中频信号为单一频率的正弦形式，信号经采样后为：
$$ x(n) = a{{\rm e} ^{j(2{\text{π }}{f_0}n/{f_{\text{s}}} + {\theta _0})}}{\text{ }}n = 1,{\text{ }}2,{\text{ }} \cdot \cdot \cdot ,{\text{ }}N - 1 $$ (3) 式中：a、f0、θ0分别为中频信号的幅度、频率和初相；fs为采样频率；N为采样点数。
序列x(n)的FFT记为：
$$ X(k) = \dfrac{{a\sin [{\text{π }}(k - {k_0})]}}{{2\sin [{\text{π }}(k - {k_0})/N]}}{{\rm e}^{j[{\theta _0} - \frac{{N - 1}}{N}(k - {k_0}){\text{π }}]}} $$ (4) 由Rife算法可得中频信号的频率估计值为：
$$ f = \dfrac{{{f_{\text{s}}}}}{N}({k_{\text{m}}} + r\dfrac{{\left| {X({k_{\text{m}}} + r)} \right|}}{{\left| {X({k_{\text{m}}} + r)} \right| + \left| {X({k_{\text{m}}})} \right|}}) $$ (5) 式中：km=int[f0T]，int[]表示取最接近的整数，km为X(k)幅度最大值对应的谱线索引值，对应频率为fm；r为修正方向；|X(km)|为频谱对应模值，|X(km+1)|> |X(km−1)|时r=1，否则r=−1。根据公式(5)可知归一化频率误差δ=(f0−km×Δf)/Δf，δ表示FFT频谱模值最大值对应的频率与实际信号频率的相对偏差，δ的取值范围为−0.5~0.5，其中Δf=fs/N为FFT的频谱分辨率。
在理想情况下，Rife算法能够很好地估计出信号频率，但是存在大噪声的情况下，|X(km+1)|> |X(km−1)|可能出现判断错误的情况，导致修正方向r出现错误，使得Rife算法产生较大误差。
-
为降低Rife算法修正方向发生错误的概率，文中结合参考文献[10]，采用相角方法求解修正方向r。记公式(4)为：
$$ \begin{split} X(k)&{\text{ = }}A(k){{\rm e}^{jB(k)}}{\text{ = }} \\ &A(k)\cos [B(k)] + jA(k)\sin [B(k)] \\ \end{split} $$ (6) 其中：
$$ A(k) = \dfrac{{a\sin [{\text{π }}(k - {k_0})]}}{{2\sin [{\text{π }}(k - {k_0})/N]}} $$ (7) $$ B(k) = {\theta _0} - \dfrac{{N - 1}}{N}(k - {k_0}){\text{π }} \approx {\theta _0} + ({k_0} - k){\text{π }} $$ (8) 由公式(6)可得X(k)的相角由A(k)和B(k)确定。不失一般性，设定δ取值范围为0~0.5，在km处有：
$$ \left\{ \begin{gathered} A({k_{\text{m}}}) > 0 \\ A({k_{\text{m}}} - 1) < 0 \\ A({k_{\text{m}}} + 1) > 0 \\ \end{gathered} \right. $$ (9) 假设B(km)位于第一象限，由公式(8)可得，B(km−1)位于第三象限，B(km+1)也位于第三象限。结合A(k)，并将B(k)代入公式(6)可得X(km)和X(km−1)的相角位于第一象限，X(km+1)的相角位于第三象限。B(km)位于其他象限时，可得出同样结论。因而可得出对应修正方向为：
$$ r{\text{ = }}\left\{ \begin{gathered} {\text{ }}1,{\text{ |}}{\theta _{{k_{\text{m}}}}} - {\theta _{{k_{{{\rm{m}}}+ 1}} }}{\text{| > |}}{\theta _{{k_{\text{m}}}}} - {\theta _{{k_{{{\rm{m}}- 1}}} }}{\text{| }} \\ - 1,{\text{ |}}{\theta _{{k_{\text{m}}}}} - {\theta _{{k_{{{\rm{m}}+ 1}}} }}{\text{| < |}}{\theta _{{k_{\text{m}}}}} - {\theta _{{k_{{{\rm{m}}- 1}}} }}{\text{|}} \\ \end{gathered} \right. $$ (10) 式中：θ为频谱对应的相角。在FPGA硬件中求解θ时，为避免涉及到计算量较大的反正切函数，将公式(4)中km处频谱值记为：
$$ X({k_{\text{m}}}){\text{ = |}}X({k_{\text{m}}})| \cdot {{\rm e}^{j{\theta _{{k_{\text{m}}}}}}} $$ (11) 同理可得X(km+1)、X(km−1)的频谱值，其中X(km)的转置X(km)T为：
$$ X({k_{\text{m}}})^{\rm T}= {\text{|}}X({k_{\text{m}}})| \cdot {{\rm e}^{ - j{\theta _{{k_{\text{m}}}}}}} $$ (12) 以X(km)为参考，进行复数相乘，即可得到θkm与左右两侧相角差：
$$ \begin{split} & X({k_{{\text{m}} - 1}}) \cdot X({k_{\text{m}}})^{\rm T}{\text{ = |}}X({k_{\text{m}}})|{\text{|}}X({k_{\text{m}}} - 1)| \cdot {{\rm e}^{j({\theta _{{k_{{{\rm{m}}- 1}}} }} - {\theta _{{k_{\text{m}}}}})}} \\ & X({k_{{\text{m}} + 1}}) \cdot X({k_{\text{m}}})^{\rm T}{\text{ = |}}X({k_{\text{m}}})|{\text{|}}X({k_{\text{m}}} + 1)| \cdot {{\rm e}^{j({\theta _{{k_{{{\rm{m}}+ 1}}} }} - {\theta _{{k_{\text{m}}}}})}} \\ \end{split} $$ (13) 此时虽然已经取得相角差，但是在实际算法实现中仍需求解反正切函数。为降低计算复杂度，根据上述X(k)相角分析可知，当0<δ<0.5时，X(km)和X(km−1)的相角位于同一象限，与X(km+1)的相角所处象限成原点对称。由此可得，θkm−θkm+1相角差位于第二或第三象限，θkm−θkm−1相角差位于第一或第四象限。由三角函数性质可得：
$$ r{\text{ = }}\left\{ \begin{gathered} {\text{ }}1,{\text{ cos(}}{\theta _{{k_{\text{m}}}}} - {\theta _{{k_{{{\rm{m}}+ 1}}} }}{\text{)}} < \cos {\text{(}}{\theta _{{k_{\text{m}}}}} - {\theta _{{k_{{{\rm{m}- 1}}}} }}{\text{) }} \\ - 1,{\text{ cos(}}{\theta _{{k_{\text{m}}}}} - {\theta _{{k_{{{\rm{m}}+ 1}}} }}{\text{)}} > {\text{cos(}}{\theta _{{k_{\text{m}}}}} - {\theta _{{k_{{{\rm{m}}- 1}}} }}{\text{) }} \\ \end{gathered} \right. $$ (14) 当−0.5<δ<0时情况类似。因此记：
$$ \begin{split} Ang\left[ k \right] =& {\rm Re}\left( {X\left( k \right) \cdot X{{\left( {{k_{\text{m}}}} \right)}^{\rm{T}}}} \right)= \\ &{\text{ |}}X\left( k \right){\text{||}}X\left( {{k_{\text{m}}}} \right){\text{|}} \cdot \cos ({\theta _k} - {\theta _{{k_{\text{m}}}}}) \\ \end{split} $$ (15) 最终可通过计算km两侧复数乘法实部确定插值修正方向r：
$$ r{\text{ = }}\left\{ \begin{gathered} {\text{ }}1,{\text{ }}Ang{\text{[}}{k_{\text{m}}} - {\text{1] > }}Ang{\text{[}}{k_{\text{m}}}{\text{ + 1] }} \\ - 1,{\text{ }}Ang{\text{[}}{k_{\text{m}}} - {\text{1] < }}Ang{\text{[}}{k_{\text{m}}}{\text{ + 1] }} \\ \end{gathered} \right. $$ (16) 将上述修正方向r代入公式(5)即可得最终频率估计值。
为研究该算法在大噪声情况下的可行性，文中基于工程中常见的高斯白噪声对其进行分析。设定加性高斯白噪声序列为z(n)，含噪声的信号序列为s(n)，对s(n)添加长度为N的对称窗wN(n)后可得加窗后的采样序列sw(n)为：
$$ {s_{\text{w}}}(n) = x(n){w_N}(n) + z(n){w_N}(n) $$ (17) 对公式(17)做FFT后可得：
$$ {S_{\text{w}}}(k) = {X_{\text{w}}}(k) + {Z_{\text{w}}}(k) $$ (18) 式中：Sw(k)、Xw(k)、Zw(k)分别为s(n)、x(n)、z(n)加窗后的FFT频谱，其中Xw(k)为无噪声信号频谱，Zw(k)为噪声频谱。
修正Rife算法通过寻找频谱幅值最大值，以确定相角并判断修正方向。不失一般性，设定δ取值范围为0~0.5，认为存在噪声情况下，Sw(k)总能在谱线km或km+1处取最大值，并记为kr。最大值寻找正确时，kr=km，寻找错误时kr=km+1，设在谱线kr处，由FFT频谱可得：
$$ {S_{\text{w}}}({k_{\text{r}}}) = \sqrt {S_{\text{w}}^{\text{R}}{{({k_{\text{r}}})}^2} + S_{\text{w}}^{\text{I}}{{({k_{\text{r}}})}^2}} $$ (19) $$ S_{\text{w}}^{\text{P}}({k_{\text{r}}}) = \arctan \left(\frac{{S_{\text{w}}^{\text{I}}({k_{\text{r}}})}}{{S_{\text{w}}^{\text{R}}({k_{\text{r}}})}}\right) $$ (20) 式中：
$S_{\rm w}^{\rm R}({k_{\rm r}}) = X_{\rm w}^{\rm R}({k_{\rm r}}) + Z_{\rm w}^{\rm R}({k_{\rm r}})$ 为kr处含噪声的离散频谱序列实部，由无噪声信号频谱实部$X_{\rm w}^{\rm R}({k_{\rm r}})$ 和噪声频谱实部序列$Z_{\rm w}^{\rm R}({k_{\rm r}})$ 组成，虚部同理可得。$S_{\rm w}^{\rm R}({k_{\rm r}})$ 为kr处含噪声的离散频谱相角序列。对公式(19)、(20)利用二元函数泰勒展开，并忽略高阶项可得：
$$ {S_{\text{w}}}({k_{\text{r}}}) \approx {X_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{) + }}Z{A_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{)}} $$ (21) $$ S_{\text{w}}^{\text{P}}({k_{\text{r}}}) \approx X_{\text{w}}^{\text{P}}{\text{(}}{k_{\text{r}}}{\text{) + }}Z{P_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{)}} $$ (22) 式中：
$X_{\rm w}^{\rm P}({k_{\rm r}})$ 为无噪声时信号的相角，当kr=km时，$X_{\rm w}^{\rm P}({k_{\rm m}}) = {\theta _{k_{\rm m}}}$ 。而ZAw(kr)和ZPw(kr)可表示为：$$ \begin{split} Z{A_{\text{w}}}({k_{\text{r}}}){\text{ = }}& \frac{{X_{\text{w}}^{\text{R}}{\text{(}}{k_{\text{r}}}{\text{)}}}}{{X_{\text{w}}^{}{\text{(}}{k_{\text{r}}}{\text{)}}}}Z_{\text{w}}^{\text{R}}{\text{(}}{k_{\text{r}}}{\text{) + }}\frac{{X_{\text{w}}^{\text{I}}{\text{(}}{k_{\text{r}}}{\text{)}}}}{{X_{\text{w}}^{}{\text{(}}{k_{\text{r}}}{\text{)}}}}Z_{\text{w}}^{\text{I}}{\text{(}}{k_{\text{r}}}{\text{)}}= \\ &{\text{ cos(}}X_{\text{w}}^{\text{P}}{\text{(}}{k_{\text{r}}}{\text{))}}Z_{\text{w}}^{\text{R}}{\text{(}}{k_{\text{r}}}{\text{)}} + {\text{sin(}}X_{\text{w}}^{\text{P}}{\text{(}}{k_{\text{r}}}{\text{))}}Z_{\text{w}}^{\text{I}}{\text{(}}{k_{\text{r}}}{\text{)}} \\ \end{split} $$ (23) $$ \begin{split} Z{P_{\text{w}}}({k_{\text{r}}}) {\text{ = }}& \frac{{X_{\text{w}}^{\text{R}}{\text{(}}{k_{\text{r}}}{\text{)}}}}{{X_{\text{w}}^2{\text{(}}{k_{\text{r}}}{\text{)}}}}Z_{\text{w}}^{\text{I}}{\text{(}}{k_{\text{r}}}{\text{)}} - \frac{{X_{\text{w}}^{\text{I}}{\text{(}}{k_{\text{r}}}{\text{)}}}}{{X_{\text{w}}^2{\text{(}}{k_{\text{r}}}{\text{)}}}}Z_{\text{w}}^{\text{R}}{\text{(}}{k_{\text{r}}}{\text{)}}{\text{ = }} \\ & \frac{{\cos (X_{\text{w}}^{\text{P}}{\text{(}}{k_{\text{r}}}{\text{)}})Z_{\text{w}}^{\text{I}}{\text{(}}{k_{\text{r}}}{\text{)}} - \sin (X_{\text{w}}^{\text{P}}{\text{(}}{k_{\text{r}}}{\text{)}})Z_{\text{w}}^{\text{R}}{\text{(}}{k_{\text{r}}}{\text{)}}}}{{X_{\text{w}}^{}{\text{(}}{k_{\text{r}}}{\text{)}}}} \\ \end{split} $$ (24) 由于噪声频谱序列Zw(k)的实部和虚部不相关，结合参考文献[13]可得幅值的相对误差百分比εabs和相角对应的相对误差百分比εph为：
$$ \begin{gathered} {\varepsilon _{{\text{abs}}}}{ = _{}} \\ \frac{{{S_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{)}} - {X_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{)}}}}{{{X_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{)}}}}\sim N(0,\frac{{{P_{\text{w}}}}}{{SNR \cdot N \cdot {W^2}(\nabla {f^1})}}) \\ \end{gathered} $$ (25) $$ \begin{gathered} {\varepsilon _{{\text{ph}}}} = \\ \frac{{S_{\text{w}}^{\text{P}}{\text{(}}{k_{\text{r}}}) - X_{\text{w}}^{\text{P}}{\text{(}}{k_{\text{r}}}{\text{)}}}}{{2{\text{π }}}}\sim N(0,\frac{{{P_{\text{w}}}}}{{4{{\text{π }}^2}N \cdot SNR \cdot {W^2}(\nabla {f^1})}}) \\ \end{gathered} $$ (26) 式中：SNR为信噪比，定义为SNR=a2/2σ2；W2(
$\nabla {f^1}$ )为归一化窗谱函数；$\nabla {f^1}$ 为窗长归一化频率误差；Pw为窗功率。由上述公式(25)、(26)可知，存在噪声环境下，基于相角的分析精度要优于幅值的分析精度。设定含噪声情况下，相角误差率门限
$\varepsilon _{\rm ph}^v$ ，根据数理统计和公式(26)可得，要使相角误差率以1−α的概率落在可接受范围，需满足条件：$$ {u_{1 - \alpha /2}}{\left[ {\frac{{{P_{\text{w}}}}}{{N \cdot SNR}}} \right]^{\tfrac{1}{2}}} \cdot \frac{1}{{\left| {W(\nabla {f^1})} \right|}} \leqslant \varepsilon _{{\text{ph}}}^{\text{v}} $$ (27) 式中：u1−α/2为正态分布1−α/2分位数。在实际FPGA频谱计算时，往往采用矩形窗方式，对应Pw=1。要想确定上述不等式关系，还需确定
$\nabla {f^1}$ 的取值范围，根据上述相角判据原理，要想取得可信结果，需满足kr在km处取得最大值，但当幅值最大值取错时，有：$$ {S_{\text{w}}}({k_{\text{r}}}) < {S_{\text{w}}}({k_{\text{r}}} + 1) $$ (28) 设幅值取错的概率为20%以内，由公式(21)可得：
$$ - 0.2{X_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{ + 1)}} \leqslant Z{A_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{ + 1)}} \leqslant 0.2{X_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{ + 1)}} $$ (29) $$ - 0.2{X_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{ + 1)}} \leqslant Z{A_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{)}} \leqslant 0.2{X_{\text{w}}}{\text{(}}{k_{\text{r}}}{\text{ + 1)}} $$ (30) 取ZAw(kr)和ZAw(kr+1)极限，代入公式(21)可得：
$$ {X_{\text{w}}}({k_{\text{r}}}) < 1.4{X_{\text{w}}}({k_{\text{r}}} + 1) $$ (31) 根据矩形窗函数原理，结合公式(31)可以得到
$\nabla {f^1}$ ≥0.4167。考虑最大值取错的情况，实际正确取值范围为[0,0.5833]，对应W(0.5833)=0.5271。取α=0.01，
$\varepsilon^{\mathrm{v}}_ \mathrm{ph}$ =20%，N=1024，当公式(27)成立时，在$\nabla {f^1}$ 的可能区间内，能够保证采用相角分析时用含噪声相角$S_{\rm w}^{\rm P}({k_{\rm r}})$ 来估计$X_{\rm w}^{\rm P}({k_{\rm r}})$ 的误差99%的可能性落在20%以内，且此时对应的信噪比小于−10 dB。综上所述，理论分析了在大噪声环境下基于相角为判据具有可行性。 -
为更直观表现该算法的特性，对文中修正算法的性能进行验证。在信噪比为−10~10 dB的范围内，设f0为15625000 Hz，采样率fs为250 MHz，FFT点数N为1024，量化频率Δf为244140 Hz。从f0−Δf/2~f0+Δf/2取1500个点，在Matlab中进行1000次蒙特卡罗仿真实验。图2为文中算法(Z-Rife)和Rife、A-P-Rife、I-Rife算法频率估计标准差在不同信噪比情况下和克拉美-罗下界 (Cramer-Rao Lower Bound,CRLB) 的比较。
图 2 4种算法频率估计标准差与CRLB比较
Figure 2. Comparison of four algorithms for standard deviation of frequency estimation with CRLB
从图2中可以看出，其他3种算法的估计精度均优于Rife算法，A-P-Rife算法优于Rife算法，低于I-Rife算法和文中算法，而文中算法在估计精度上逼近I-Rife算法，接近CRLB。
上述几种算法均基于FFT频谱进行插值计算，为直观反映不同算法的计算复杂度，表1给出了不同算法基于N点FFT所增加的运算量。
表 1 4种算法基于N点FFT增加计算量对比
Table 1. Calculation comparison of four algorithms based on N-point FFT
文中算法基于相角实现修正方向的确定，计算公式(15)时，仅需1次复数乘法和加法。由于确定修正方向涉及到计算谱线km两侧的相角，所以完成最终频率估计共需要2次乘法、1次除法和4次加法。计算量远小于I-Rife算法，略小于A-P-Rife算法，相对于FFT运算所增加计算量可忽略不计。
-
文中采用FPGA为平台对修正Rife算法进行硬件实现，主要由以下步骤完成：
(1) 在空闲状态下检测到FFT输出有效信号时，算法开始以流水方式检测FFT输出频谱的最大值。为缓存最大值及两侧数据，采用三级寄存器级联，并利用峰值检测电路检测三级寄存器的中间寄存器值是否为最大值，如果是，则将三级寄存器数据更新，否则，保持原来数据。待有效信号结束时，即可获取X(km)、X(km−1)、X(km+1)和对应模值。
(2) FFT有效信号结束后，跳转至该步骤，通过式(15)计算左右两侧谱线对应的相角差值Ang(km−1)、Ang(km+1)，由于Ang的最终值为复数的实部，因此只需求解对应的实部即可，避免资源浪费。
(3) 根据Ang(km−1)和Ang(km+1)的大小确定修正方向r，并将步骤(1)缓存的频谱模值代入公式(5)即可求得最终的频率估计值。
-
为验证该算法整体性能，基于实验室现有条件，搭建了如图3所示光学测试平台进行测试。测试平台由以下模块组成，其中①~⑧分别是：稳压源、任意函数发生器、DFB激光源、电光强度调制器、准直器、光电探测器、算法硬件处理平台和上位机。
图 3 系统测试平台
Figure 3. Test platform of the system
-
为验证文中修正Rife算法精度，以传统Rife算法为参考，采用函数发生器为信号源，采样点数N=1024，在SNR为−10 dB情况下以FFT量化频率点左右栅栏中心为界，等间隔选取15个频率点，对Rife算法和文中算法(Z-Rife)进行精度测试[14]。
根据图4计算出2种算法的均方根误差和平均误差，如表2所示。从表中可知文中算法的均方根误差(RMSE)和平均误差(ME)相较于Rife算法分别降低了50.7%和69.6%，相对Rife算法文中算法在频率求解精度方面有明显提升。
图 4 2种算法绝对频率误差分布图
Figure 4. Distribution diagram of two algorithms absolute frequency error
表 2 不同算法的测频误差
Table 2. Frequency error of different algorithms
Algorithm RMSE/Hz ME/Hz Rife 37258 35598 Z-Rife 18371 10821 -
为验证该算法在系统中的处理精度，利用如图3所示实验测试平台，模拟生成FMCW激光雷达光信号测试该算法。其中DFB激光器激光波长λ为1550 nm，系统调频时间T设定为10 μs，B为1 GHz。结合激光雷达传输方程[15]，考虑不同距离情况下激光雷达衰减，设置电光调制器工作在正交偏置点，并根据FMCW激光雷达原理生成如图5所示激光回波信号，将信号输入任意函数发生器后输出至电光强度调制器，实现FMCW激光雷达回波光信号输出，光信号经光电探测模块接收后输入至算法硬件处理平台进行处理。
图 5 FMCW激光雷达中频回波信号
Figure 5. FMCW IF echo signal
根据实际城市道路速度要求和城市中行人、非机动车辆和机动车辆等不同场景，给出不同距离和相对速度，并与算法实际测试数据进行对比，实测结果如表3所示。
表 3 实验结果明细表
Table 3. Experimental results
L
/mV
/km·h−1L’
/mV’
/km·h−1L_error
/mV_error
/km·h−11 10 0.98 9.98 0.02 0.02 5 20 4.98 19.97 0.02 0.03 20 30 20.01 29.84 0.01 0.16 20 60 19.99 60.04 0.01 0.04 20 90 20.04 90.00 0.04 0.00 50 30 50.04 29.99 0.04 0.01 50 90 50.00 90.01 0.00 0.01 50 120 50.02 120.01 0.02 0.01 80 30 80.03 29.87 0.03 0.13 80 90 80.05 89.98 0.05 0.02 80 140 80.04 140.01 0.04 0.01 112 30 112.04 29.99 0.04 0.01 112 90 112.02 89.95 0.02 0.05 112 140 120.05 139.85 0.05 0.15 表中，L和V为实际目标的距离和速度，L’和V’为算法实测结果，L_error和V_error为对应距离和速度测试误差。
通过实验可知所提算法能够实时解算对应目标的距离和速度，同时经过数据分析，该算法由于ADC采样率的限制，最大测距范围为112 m，测速范围满足城市道路70 km/h的需求，其中该算法延迟约为100 μs，测距误差不大于5 cm，测速误差不大于0.16 km/h。由于该修正算法在进行频率估计时随着频率和噪声不同，估计误差会有不同表现，且由图4中修正Rife算法最大频率误差可知，该算法对应最大测距误差不大于5 cm，与表3测试结果相符。
-
为满足FMCW激光雷达系统高精度测距测速的需求，文中通过对正弦波的幅值和相角特性进行分析，提出了一种易于硬件实现的修正Rife算法。经Matlab仿真和计算量分析，该算法精度优于传统Rife算法和A-P-Rife算法，略小于I-Rife算法，且更接近CRLB，在计算量上小于I-Rife算法和A-P-Rife算法，略大于Rife算法，整体结果表明该算法性能最佳。
搭建测试平台，将该算法进行硬件移植，比较该算法和传统Rife算法的精度，结果显示：该算法相较于传统Rife算法在计算量仅增加两个乘法和加法的情况下平均误差降低了69.6%，均方根误差降低了50.7%。同时，模拟FMCW激光雷达回波光信号完成该算法测试，测试结果显示该算法能够处理中频信号频谱信息，在112 m范围内能够满足实时测距测速需求，测距误差不大于5 cm，测速误差不大于0.16 km/h。综上，该算法具有较高的工程应用价值。
Design and hardware implementation of modified Rife algorithm for FMCW LiDAR
-
摘要: FMCW激光雷达以其高精度、抗干扰能力强、同时测距测速等特点得到了广泛研究。针对FFT固有栅栏效应引入测距、测速误差的问题，通过分析频谱幅值和相角的规律，并结合正弦函数原理提出了一种易于硬件实现的修正Rife算法，有效地降低了传统Rife算法在估计频率接近FFT量化频率点时的误差。通过仿真和FPGA验证，修正Rife算法在信噪比为−10 dB时相较于传统Rife算法平均误差降低了69.6%，均方根误差降低了50.7%，而计算量仅增加了两个乘法和加法，与N点FFT计算量相比可忽略不计。最后，通过搭建光学测试平台，模拟激光雷达中频回波信号验证了该算法的有效性。测试结果显示，该算法可在112 m范围内实现同时测距测速，测距误差不大于5 cm，测速误差不大于0.16 km/h，满足实时性要求。Abstract: FMCW LiDAR is being widely studied for its high accuracy, strong anti-interference capability and simultaneous ranging and speed measurement. The inherent fence effect of FFT will introduce errors in ranging and speed measurement. To solve this problem, firstly, this paper analyzes the law of spectrum amplitude and phase angle, and then combines with the principle of sine function, proposes a modified Rife algorithm that is easy to implement in hardware. When the estimated frequency is close to the FFT quantization frequency point, this method can effectively reduce the error of the Rife algorithm. The simulation and FPGA verification results show that when the SNR is −10 dB, the mean error (ME) and root mean square error (RMSE) of the algorithm are reduced by 69.6% and 50.7%, respectively, compared with the traditional Rife algorithm. The calculation amount is only increased by two multiplications and additions, which is negligible compared to N-point FFT computation. Finally, in order to verify the effectiveness of the algorithm, this paper build an optical test platform to simulate the intermediate frequency echo signal of LiDAR. The results show that the algorithm can achieve simultaneous ranging and speed measurement within 112 m. The ranging error is no more than 5 cm, and the speed measurement error is no more than 0.16 km/h. The algorithm meets the demand of real-time ranging and speed measurement.
-
Key words:
- FMCW LiDAR /
- frequency estimation /
- modified Rife algorithm /
- FPGA
-
图 1 (a) 梯形波形式调制下激光发射信号和激光回波信号时频域关系图；(b) 经相干检测所得差拍中频信号时频域关系图；(c) 差拍中频电信号的时域波形
Figure 1. (a) Time frequency domain diagram of transmitted signal and received echo signal under the trapezoidal modulation; (b) Time frequency domain diagram of the beat signal after coherent detection; (c) Time domain waveform of the received beat signal
图 2 4种算法频率估计标准差与CRLB比较
Figure 2. Comparison of four algorithms for standard deviation of frequency estimation with CRLB
图 4 2种算法绝对频率误差分布图
Figure 4. Distribution diagram of two algorithms absolute frequency error
表 1 4种算法基于N点FFT增加计算量对比
Table 1. Calculation comparison of four algorithms based on N-point FFT
下载: 导出CSV
表 2 不同算法的测频误差
Table 2. Frequency error of different algorithms
Algorithm RMSE/Hz ME/Hz Rife 37258 35598 Z-Rife 18371 10821
下载: 导出CSV
表 3 实验结果明细表
Table 3. Experimental results
L
/mV
/km·h−1L’
/mV’
/km·h−1L_error
/mV_error
/km·h−11 10 0.98 9.98 0.02 0.02 5 20 4.98 19.97 0.02 0.03 20 30 20.01 29.84 0.01 0.16 20 60 19.99 60.04 0.01 0.04 20 90 20.04 90.00 0.04 0.00 50 30 50.04 29.99 0.04 0.01 50 90 50.00 90.01 0.00 0.01 50 120 50.02 120.01 0.02 0.01 80 30 80.03 29.87 0.03 0.13 80 90 80.05 89.98 0.05 0.02 80 140 80.04 140.01 0.04 0.01 112 30 112.04 29.99 0.04 0.01 112 90 112.02 89.95 0.02 0.05 112 140 120.05 139.85 0.05 0.15
下载: 导出CSV
-
[1] Wu Xiru, Xue Qiwei. 3D vehicle detection for unmanned driving system based on LiDAR [J]. Optics and Precision Engineering, 2022, 30(4): 489-497. (in Chinese) doi: 10.37188/OPE.20223004.0489 [2] Wei Yu, Jiang Shilei, Sun Guobin, et al. Design of solid-state array laser radar receiving optical system [J]. Chinese Optics, 2020, 13(3): 517-526. (in Chinese) [3] Ahmod Z, Liao Y M, Kuo S I, et al. High-power and high-responsivity avalanche photodiodes for self-heterodyne FMCW LiDAR system applications [J]. IEEE Access, 2021, 9: 85661-85671. doi: 10.1109/ACCESS.2021.3089082 [4] Wang C S, Zhang Y S, Zheng J L, et al. Frequency modulated continuous wave dual frequency LiDAR based on a monolithic integrated two-section DFB laser [J]. Chinese Optics Letters, 2021, 19(11): 111402. doi: 10.3788/COL202119.111402 [5] Gao Li, Zhang Xiaoli, Ma Jingting, et al. Quantum enhanced doppler LiDAR based on integrated quantum squeezed light source (Invited) [J]. Infrared and Laser Engineering, 2021, 50(3): 20210031. (in Chinese) doi: 10.3788/IRLA20210031 [6] Chen Peng, Zhao Jiguang, Song Yishuo, et al. Comparison on detection performance of FMCW and pulsed LiDAR in aerosol environment [J]. Infrared and Laser Engineering, 2020, 49(6): 20190399. (in Chinese) doi: 10.3788/irla.22_2019-0399 [7] Qu Xinghua, Zhi Guangtao, Zhang Fumin, et al. Improvement of resolution of frequency modulated continuous wave laser ranging system by signal splicing [J]. Optics and Precision Engineering, 2015, 23(1): 40-47. (in Chinese) doi: 10.3788/OPE.20152301.0040 [8] Zhou Yongsheng, Ma Xunpeng, Zhao Yiming, et al. Frequency estimation of the weak signal of the coherent wind LiDAR [J]. Infrared and Laser Engineering, 2018, 47(3): 0306002. (in Chinese) doi: 10.3788/IRLA201847.0306002 [9] Zou Lirong, Zhu Li, Zhang Sheng, et al. LFMCW radar moving target ranging method based on improved rife algorithm [J]. Journal of Microwaves, 2019, 35(S1): 233-238. (in Chinese) [10] Sun Hongjun, Wang Xiaowei. Modified rife algorithm for frequency estimation of sinusoid wave based on amplitude and phase criterion [J]. Journal of Tianjin University (Science and Technology), 2018, 51(8): 810-816. (in Chinese) [11] Cui Songqi, Wang Aihua, Tang Le, et al. Frequency estimation of sinusoid based on the dual-threshold decision modified Rife algorithm [J]. Transactions of Beijing Institute of Technology, 2017, 37(3): 292-296. (in Chinese) doi: 10.15918/j.tbit1001-0645.2017.03.013 [12] Lim H S, Park H M, Lee J E, et al. Lane-by-lane traffic monitoring using 24.1 GHz FMCW radar system [J]. IEEE Access, 2021, 9: 14677-14687. doi: 10.1109/ACCESS.2021.3052876 [13] 丁康, 谢明, 杨志坚, 离散频谱分析校正理论与技术[M]. 北京: 科学出版社, 2008: 170-192. Ding Kang, Xie Ming, Yang Zhijian. The Theory and Technology of Discrete Spectrum Correction[M]. Beijing: Science Press, 2008: 170-192. (in Chinese) [14] Zhan Qidong, Tu Yaqing. Analysis and implementation of Rife-based ranging algorithm for linear frequency modulated continuous wave radar [J]. Acta Armamentarii, 2014, 35(5): 748-752. (in Chinese) [15] Sun Junling, Ma Pengge, Sun Guangmin, et al. Multi-pulse laser radar target signal simulation based on target echo waveform model [J]. Infrared and Laser Engineering, 2016, 45(7): 0726006. (in Chinese) doi: 10.3788/irla201645.0726006 -
点击查看大图
计量
- 文章访问数: 203
- HTML全文浏览量: 67
- PDF下载量: 95
- 被引次数: 0
