1.   引言

中国余数定理(Chinese Remainder Theorem, CRT)由于其“化整为零”的特性，使得该定理在数字图像信号处理^[1]、编码^[2]、密码学^[3,4]及目标参数估计^[5]等领域获得了广泛应用。文献[6,7]将CRT应用到相位解模糊，提出的CRT算法不具备鲁棒性，对得到的余数中任一个略有误差就会导致巨大的重构错误，不能直接应用于有噪声的环境下。为了克服算法对余数噪声敏感的缺陷，一种基于搜索的鲁棒CRT重构算法在文献[8,9]中被提出，但是这种算法计算量非常大，虽然文献[10]指出将2维降低为1维搜索，当需要搜索的数量很大的时候，计算复杂度依然很高，在实际的工程应用中不是很理想。文献[11]提出的鲁棒闭式CRT方法克服运算量大的问题，但是选取第1个余数作为参考构造差分余数方程组，在重构多普勒频率时没有考虑余数误差对精度的影响，估计精度有一定的缺失。

本文基于脉冲多普勒(Pulsed Doppler, PD)雷达参差重频模式下，提出一种基于全相位离散傅里叶变换(all phase Discrete Fourier Transform, apDFT)相位差频谱校正^[12-17]的最优余数封闭式鲁棒中国余数定理(Closed-Form Robust Chinese Remainder Theorem, CFRCRT)算法^[18]，不仅克服了经典CRT算法对余数噪声敏感和搜索CRT算法计算量大的问题，并充分考虑余数误差的影响，选取最优余数作为参考构建差分余数方程组，而且通过apDFT相位差频谱校正准确提供频率余数，在很大程度上提高了多普勒频率测量精度。理论分析和仿真实验表明本文提出的基于相位差频谱校正的CFRCRT算法在测量精度和实时性能上可以满足工程上应用的需求。

2.   脉冲多普勒雷达频率估计模型

PD雷达在一个相干处理间隔(Coherent Processing Interval, CPI)内按脉冲重复频率(Pulse Repetition Frequency, PRF)发射脉冲，假设一个CPI为T，接收机接收到的回波信号是离散序列$x(n)$。

$$ \begin{aligned} x(n) \!=\! & s\left(n\right) \!+\! w\left(n\right) = a \cdot \exp \left({\rm{j}}2{\text{π}}F{t_{\rm s}}n\right) + w\left(n\right), \\ {\rm{}}& 0 \le n \le M - 1 \end{aligned} $ $

(1)

其中${t_{\rm{s}}} = {\rm PRT}\,$是采样时间间隔，即${f_{\rm{s}}}$为采样频率，在时间$T\;$内采样总数$M = {f_{\rm{s}}} \cdot T = {\rm{PRF}} \cdot T\,$, $s(n)$表示采样后的目标信号，$w(n)$表示噪声信号，$F\,$是待测的多普勒频率。$x(n)$的离散傅里叶变换(Discrete Fourier Transform, DFT)为

$$ \begin{aligned} X\left(k\right) = & \displaystyle\sum\limits_{n = 0}^{M - 1} {x\left(n\right)} W_M^{{\rm{nk}}} = \displaystyle\sum\limits_{n = 0}^{M - 1} {x\left(n\right)} \\ {\rm{}}& \cdot \exp {\left( - {\rm{j}}2{\text{π}}kn/M\right)}, 0 \le n \le M - 1 \end{aligned} $ $

(2)

如果${f_{\rm{s}}} < F\;$，估计值${\hat F}\;$在频谱上是模糊的，不能按照直接搜索峰值位置来估计多普勒频率。实际上，由于DFT的周期性，DFT频谱上的视在频率是真实多普勒频率$F\,$以${\rm{PRF}}\left( {{f_{\rm{s}}}} \right)$为模的余数，即

$$F \equiv \displaystyle \frac{{\hat r}}{T\;}od {f_{\rm{s}}}$ $

(3)

脉冲多普勒雷达出现速度模糊的原因在于运动目标真实多普勒频率超出了${\rm{PRF}}$下脉冲重频范围。为了有效解决速度模糊问题，脉冲多普勒雷达采用参差周期技术。在多重参差重频下，目标在各重频上测得存在模糊的多普勒频移，但是这些模糊值与真实值存在式(3)的函数关系，即有式(4)所示的同余方程组。

$$\left. {\begin{array}{l} {R = {n_1}{M_1} + {r_1}} \\ {R = {n_2}{M_2} + {r_2}} \\ \quad\quad\quad\; \vdots \\ {R = {n_L}{M_L} + {r_L}} \end{array}} \right\}$ $

(4)

其中，$R = F \!\cdot\! T,{M_i} \!=\! {\rm{PRF}}{_i} \!\cdot\! T \!=\! {f_{{\rm{s}}i}} \!\cdot\! T\;$, ${n_i}(1 \le i \le L)$是同余方程组的模糊度。

3.   基于频率校正的中国余数定理多普勒频率估计算法

3.1   封闭式鲁棒中国余数定理

由式(3)可以得到

$$F = \displaystyle \frac{1}{T\;}\left({n_i}{M_i} + {r_i}\right)$ $

(5)

为了解决经典CRT算法对噪声敏感的问题，封闭式鲁棒中国余数定理算法采用先确定式(5)中的模糊整数${n_i}$，所以PD雷达解多普勒频率模糊的问题回归到求解模糊度${n_i}(1 \le i \le L)$大小。不失一般性，假设${M_i}(1 \le i \le L)$之间存在最大公约数$G$，即$G{\rm{ = }}\gcd ({M_1},{M_2}, ·\!·\!· ,{M_L})$。令${{\text{Γ}}_{i}}={{{M}_{i}}}/{G}\;,1\le i\le L$，则所有${\text{Γ} _i}(1 \le i \le L)$都是互素的。式(4)可以等价为

$$\left. {\begin{array}{l} {{R_0} = {n_1}{\text{Γ} _1} + {q_1}} \\ {{R_0} = {n_2}{\text{Γ} _2} + {q_2}} \\ \quad\quad\quad \vdots \\ {{R_0} = {n_L}{\text{Γ} _L} + {q_L}} \end{array}} \right\}$ $

(6)

其中，${R_0}{\rm{ = }}\left\lfloor {\displaystyle\frac{R}{G}} \right\rfloor ,{q_i} = \left\lfloor {\displaystyle\frac{{{r_i}}}{G}} \right\rfloor (1 \le i \le L)$, $\left\lfloor \cdot \right\rfloor$表示向下取整运算，可以等价为

$$\left. {\begin{array}{l} {R = {R_0}G + {r^{\rm{c}}}} \\ {{r_i} = {q_i}G + {r^{\rm{c}}}} \end{array}} \right\}$ $

(7)

式(7)表明无误差的余数${r_i}$对模$G$有共同的余数${r^{\rm{c}}}$。

为了构建最优余数差分方程组，在余数存在误差的情况下需要选取最优余数作为参考。根据式(7)可知，在没有余数误差的情况下，余数${r_i}$对模$G$有共同的余数${r^{\rm{c}}}$，但是余数有误差时${r^{\rm{c}}}$不相同，所以首先需要从存在误差的余数${r_i}$中估计出最优${r^{\rm{c}}}$。由式(7)可知存在误差的余数对模$G$的余数为

$$\hat r_i^{\rm{c}} = {r_i}\;od \;G$ $

(8)

由于${r^{\rm{c}}}$是折叠余数，故定义特殊距离函数估计${r^{\rm{c}}}$

$${\hat r^{\rm{c}}} = \arg \;\mathop {\min }\limits_{\begin{array}{*{20}{c}} {0 \le m \le G - 1}\\ {\;\;\;m \in Z} \end{array}} \displaystyle \sum\limits_{i = 1}^L {{{\left( {\hat r_i^{\rm{c}} - m - {k_0}G} \right)}^2}} $ $

(9)

其中，${k_0} = \arg \!\!\mathop {\min }\limits_{\begin{array}{*{20}{c}} {k \in \left\{ { - 1,0,1} \right\}}\\ {0 \le n \le G - 1} \end{array}} \!\!\left| {m - n - kG} \right|$。则最优参考余数索引为

$${i_0} = \arg \;\mathop {\min }\limits_{1 \le i \le L} {\left( {r_i^{\rm{c}} - {{\hat r}^{\rm{c}}} - {k_0}G} \right)^2}$ $

(10)

将式(6)中$L$个等式分别减去第${i_0}$个等式，得

$$\left. {\begin{array}{l} {{n_{{i_0}}} \cdot {\text{Γ} _{{i_0}}} - {n_1} \cdot {\text{Γ} _1} = {q_{1,{i_0}}}} \\ { {n_{{i_0}}} \cdot {\text{Γ} _{{i_0}}} - {n_2} \cdot {\text{Γ} _2} = {q_{2,{i_0}}}} \\ \quad\quad\quad\quad\quad \vdots \\ {{n_{{i_0}}} \cdot {\text{Γ} _{{i_0}}} - {n_L} \cdot {\text{Γ} _L} = {q_{L,{i_0}}}} \end{array}} \right\}$ $

(11)

其中，${{q}_{i,{{i}_{0}}}}={{q}_{i}}-{{q}_{{{i}_{0}}}}=\left\lfloor {\left( {{r}_{i}}-{{r}_{{{i}_{0}}}} \right)}/{G}\; \right\rfloor$。

引理 1　如果${\text{Γ} _{{i_0}}}\left( {1 \le {i_0} \le L} \right)$和${\text{Γ} _i}(1 \le i \le L)$是互素的，即${\text{Γ} _i}(1 \le i \le L)$两两之间的公约数都是1，对于式(11)的方程组的解$({{\hat n}_{{i_0}}},{{\hat n}_i})$为

$$ \left. \begin{aligned} {\rm{}}& {{{{\hat n}}_{{i_0}}} = {q_{i,{i_0}}}{\text{Γ} _{i,{i_0}}} + k{\text{Γ} _i}} \\ {\rm{}} & {{{{\hat n}}_i} = \displaystyle \frac{{{q_{i,{i_0}}}{{{\bar {\text{Γ}}} }_{i,{i_0}}}{\text{Γ} _{{i_0}}} - {q_{i,{i_0}}}}}{{{\text{Γ} _i}}} + k{\text{Γ} _{{i_0}}}} \end{aligned} \right\} $ $

(12)

其中，${{\bar{ \text{Γ}}} _{i,{i_0}}}$是${\text{Γ} _{{i_0}}}$对于模${\text{Γ} _i}$的模逆元素。

证明：由于${\text{Γ} _{{i_0}}}$和${\text{Γ} _i}(1 \le i \le L,i \ne {i_0})$互素，根据贝祖定理可以得出${\text{Γ} _{{i_0}}}$对于${\text{Γ} _i}(1 \le i \le L)$的模逆元素存在，用${{\bar {\text{Γ}}} _{i,{i_0}}}$表示${\text{Γ} _{{i_0}}}$对于模${\text{Γ} _i}$的模逆元素。即

$$ {{\bar {\text{Γ}}} _{i,{i_0}}} \cdot {\text{Γ} _{{i_0}}} \equiv 1\;od \;{\text{Γ} _i} $ $

(13)

考虑式(11)中的第$i$个方程

$${n_{{i_0}}}{\text{Γ} _{{i_0}}} - {n_i}{\text{Γ} _i} = {q_{i,{i_0}}}$ $

(14)

将式(14)两端分别乘以${{\bar {\text{Γ}}} _{i,{i_0}}}$，然后取模${\text{Γ} _i}$运算

$$ {{\bar {\text{Γ}}} _{i,{i_0}}}({n_{{i_0}}}{\text{Γ} _{{i_0}}} - {n_i}{\text{Γ} _i})od {\text{Γ} _i} = {{\bar {\text{Γ}}} _{i,{i_0}}}{q_{i,{i_0}}}od {\text{Γ} _i} $ $

(15)

对式(15)整理可以得到

$$ {n_{{i_0}}} = {q_{i,{i_0}}} \cdot {{\bar {\text{Γ}}} _{i,{i_0}}}od {\text{Γ} _i} $ $

(16)

所以${n_{{i_0}}}$解的通式可以表示为

$$ {n_{{i_0}}} = {{\bar {\text{Γ}}} _{i,{i_0}}}{q_{i,{i_0}}} + k{\text{Γ} _i},\;\;\;k \in Z $ $

(17)

将式(17)代入式(14)可得

$$ {n_i} = \displaystyle \frac{{{q_{i,{i_0}}}{{{\bar {\text{Γ}}} }_{i,{i_0}}}{\text{Γ} _{{i_0}}} - {q_{i,{i_0}}}}}{{{\text{Γ} _i}}} + k{\text{Γ} _{{i_0}}},\;\;\;\;k \in Z,i \ne {i_0} $ $

(18)

即$\left( {{n_{{i_0}}},{n_i}} \right)$有具体形式

$$ \left. {\begin{array}{*{20}{c}} {{n_{{i_0}}} = {{{\bar {\text{Γ}}} }_{i,{i_0}}}{q_{i,{i_0}}} + k{\text{Γ} _i}{\rm{ }}}\qquad\qquad \\ {{n_i} = \displaystyle \frac{{{q_{i,{i_0}}}{{{\bar {\text{Γ}}} }_{i,{i_0}}}{\text{Γ} _{{i_0}}} - {q_{i,{i_0}}}}}{{{\text{Γ} _i}}} + k{\text{Γ} _{{i_0}}}} \end{array}} \right\} $ $

(19)

证毕

引理 2　假设所有的${\text{Γ} _i}(1 \le i \le L)$是两两互素的，且

$$0 \le F \le {\rm{l\;cm}}({M_1},{M_2}, ·\!·\!· ,{M_L}) = G{\text{Γ} _1}{\text{Γ} _2} ·\!·\!· {\text{Γ}_L}$ $

(20)

如果$\tau <{G}/{4}\;$(其中$\tau$是余数噪声${r_i}\left( {1 \le i \le L} \right)$的最大错误范围)，那么就有${{\hat n}_i} = {n_i},1 \le i \le L$。封闭式鲁棒CRT算法就可以准确地重构估计参数。多普勒频率的估计值${\hat F}$可以表示为

$$ {\hat F} = \frac{1}{T\;}\left[ {\frac{1}{L}\left( {\sum\limits_{i = 1}^L {\left( {{{{\hat n}}_i}{M_i} + {{{\bar r}}_i}} \right)} } \right)} \right] $ $

(21)

其中，$\left[ \cdot \right]$表示四舍五入运算，${{\bar r}_i}$表示重构的最优余数。当余数最大噪声误差$\tau <{G}/{4}\;$成立时，运用封闭式鲁棒CRT算法重构的多普勒频率的误差满足下列关系

$$\left| {F - {\hat F}} \right| \le \frac{\tau }{T\;}$ $

(22)

本文提出的封闭式鲁棒CRT具有封闭形式的解，不需要搜索操作，所以计算复杂度有明显的降低。封闭式鲁棒性的CRT算法(CFRCRT)的具体步骤。

步骤 1　依据式(8)计算$\hat r_i^{\rm{c}}$；

步骤 2　依据式(9)和式(10)计算最优余数的索引号${i_0}$；

步骤 3　由所给出的含有误差的余数${r_i}$计算出${q_{i,{i_0}}}\left( {1 \le i \le L,i \ne {i_0}} \right)$；

步骤 4　计算${{\hat \xi}_{i,{i_0}}} \equiv {\hat q_{i,{i_0}}}{{\bar {\text{Γ}}} _{i,{i_0}}}od {\text{Γ} _i}\left( {1 \!\le\! i \!\le\! L,}\right.$$\left.{i \ne {i_0}} \right)$；

步骤 5　计算${{\hat n}_{{i_0}}}\! =\! \!\left(\! {\displaystyle \sum\nolimits_{i = 1}^L \!\!{{{\hat \xi }_{i,{i_0}}}{b_{i,{i_0}}}\displaystyle \frac{{{\gamma _{{i_0}}}}}{{{\text{Γ} _i}}}} } \!\right)\!\!od {\gamma _{{i_0}}}$，其中${b_{i,{i_0}}} = \arg \;\mathop {\min }\limits_{b \in Z} \left( {b \cdot \displaystyle \frac{{{\gamma _{{i_0}}}}}{{{\text{Γ} _i}}} \equiv 1od {\text{Γ} _i}} \right)$；

步骤 6　 计算${{\hat n}_i} \,=\, \displaystyle \frac{{{{{\hat n}}_{{i_0}}}{\text{Γ} _{{i_0}}} - {{\hat q}_{i,{i_0}}}}}{{{\text{Γ} _i}}}\,\left( {1 \,\le\, i \le L,}\right.$$\left. {{i \ne i_0}} \right)$；

步骤 7　计算最优余数值${{\bar r}_i} \!=\! G \cdot \left\lfloor {\displaystyle \frac{{{r_i}}}{G}} \right\rfloor + r_{{i_0}}^{\rm{c}}$；

步骤 8　计算${\hat F} = \displaystyle \frac{1}{T\;}\left[ {\frac{1}{L}\sum\nolimits_{i = 1}^L {\left( {{{{\hat n}}_i}{M_i} + {{{\bar r}}_i}} \right)} } \right]$。

3.2   相位差频谱校正算法

由于DFT的频谱泄漏和栅栏效应，在样本序列点数较少的情况下，从各路DFT得到的模糊频率估计值只能精确到整数，而丢弃小数，无法提供准确的CRT频率余数，这必然影响重构的多普勒频率精度。出于上述考虑，本文引入相位差频谱校正算法，提高多普勒频率的估计精度。

3.2.1   全相位谱分析

传统DFT谱分析出现频谱泄漏的主要原因是“截断”。然而对于一个离散的序列，取样的位置实际上代表了相位信息。如果考虑到样本所有可能的相位，则此频谱分析的结果是每个点的相位都是一样的。基于这样的考虑定义apDFT。

考虑一个有限长序列

$$ x(n) = \left\{ {\begin{array}{*{20}{c}} {x(n),\;\;\; - N + 1 \le n \le N - 1}\\ {{\rm{0}}\;,\;\;\;\;\;\;\;{\text{其它}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\quad\,} \end{array}} \right. $ $

(23)

对式(23)的有限长序列进行加窗处理

$${x_m}(n) = x(n){w_1}(n + m),\;\;0 \le m \le N - 1$ $

(24)

对式(24)的$m$个序列进行$N$点的周期延拓

$${\tilde x_m}(n) = \sum\limits_{r = - \infty }^\infty {{x_m}} (n + rN),\;\;0 \le m \le N - 1$ $

(25)

对式(25)的$m$个周期序列分别进行加窗处理

$${y_m}(n) \!=\! \displaystyle \sum\limits_{r \!=\! - \infty }^\infty {{w_2}(n){x_m}} (n \!+\! rN),\;\;0 \le m \le N - 1$ $

(26)

定义$N$点的全相位信号序列为

$$ {x_{{\rm{ap}}}}(n) = \left\{ {\begin{array}{*{20}{c}} {\sum\limits_{m = 0}^{N - 1} {{y_m}(n),\;\;0 \le n \le N - 1} }\\ {0,\;\;\;\;\;\;\;\;\;\;\;\;{\text{其它}}\;\;\;\;\;\;\;\;\;\;\;\;} \end{array}} \right. $ $

(27)

对${x_{{\rm{ap}}}}(n)$做DFT变换

$$ \begin{aligned} {\rm{}}&{Y_{\rm{a}}}(k) = \displaystyle \frac{{{\rm a} \cdot {{\rm e}^{{\rm{j}}{\theta _0}}}}}{{{N^2}}}{\left( {\displaystyle \frac{{\sin \left[ {(\beta - k){\text{π}}} \right]}}{{\sin \left[ {(\beta - k){\text{π}}/N} \right]}}} \right)^2},\\ {\rm{}}& k = 0,1, ·\!·\!· ,N - 1 \end{aligned} $ $

(28)

可以得到相位谱为

$$ {\phi _{\rm a}}(k) = {\theta _0},\;\;\;k = 0,1, ·\!·\!· ,N - 1 $ $

(29)

令$N \!=\! 8,A \!=\! 1,{w_0} \!=\! \beta \!\cdot\! \displaystyle \frac{{2{\text{π}}}}{8},{\theta _0} = {30^ \circ }$，传统DFT和全相位DFT仿真如图1和图2。

图 1  传统DFT幅频响应和相频响应

下载: 全尺寸图片 幻灯片

图 2  全相位DFT幅频响应和相频响应

下载: 全尺寸图片 幻灯片

综上所述，全相位DFT在分析单频复正弦信号的频谱特性时有以下优点：

(1)相位谱${\phi _{\rm a}}(k)$与频率指数$k$是相互独立的，相位信息可以在整个频率轴中抽取。

(2)全相位DFT旁瓣衰减的速度是传统DFT旁瓣的衰减速度的2次方，全相位DFT频谱泄漏的程度要远远小于传统DFT。

3.2.2   相位差频率校正原理

信号的频率为$\beta \Delta w$，其中估计值$\beta$可以分成两部分，整数部分${k^ {*} }$和小数部分$\Delta {k^ {*} }$，即

$$\beta = {k^ {*} } + \Delta {k^ {*} }$ $

(30)

利用传统DFT和全相位DFT频谱分析，可以搜索到峰值位置，然后在峰值位置求取相位差

$$ \begin{aligned} \Delta \phi \left({k^ {*} }\right) = & {\phi _X}\left({k^ {*} }\right) - {\phi_{\rm a}}\left({k^ {*} }\right) = \left(1 - 1/N\right) \\ {\rm{}}& \cdot\left(\beta - {k^ {*} }\right){\text{π}} = \left(1 - 1/N\right)\Delta {k^ {*} }{\text{π}} \end{aligned} $ $

(31)

所以$\beta$的估计值是

$$\hat \beta = {k^*} + \hat \delta = {k^ {*} } + \frac{{\Delta \phi ({k^ {*} })N}}{{(N - 1){\text{π}}}}$ $

(32)

其中第2项$\hat \delta = \displaystyle \frac{{\Delta \phi \left( {{k^*}} \right)N}}{{(N - 1){\text{π}}}}$称为频率误差估计，该项是由传统DFT和apDFT的相位差得到，这将明显提高最后的频率估计精度。相位差频谱校正是基于全相位频谱分析的一种算法，由于apDFT具有恒定不变的相位的特性，所以其估计的频率均方根误差和AM算法基本相同，即得到的估计均方根误差在最好的情况下大约为克拉美罗下界的1.01倍。而且相位差频偏估计中采用两次快速傅里叶变换，所以计算复杂度为$O\left( {3N\lg N} \right)$，实时性可以满足工程的需求。

3.3   基于频谱校正的CFRCRT多普勒频率估计方法

CFRCRT多普勒频率估计方法的流程如图3所示。

图 3  多普勒频率估计流程图

下载: 全尺寸图片 幻灯片

在脉冲多普勒雷达系统中，发射机采用参差周期工作方式，即采用多个重复频率周期工作。为了得到$L$路频率余数信息，在一个CPI内，即相等时间$T\,$内，分别采用不同的${\rm{PR}}{{\rm{F}}_i}(1 \le i \le L)$发射脉冲信号。具体步骤如下所述：

步骤 1　发射机以脉冲重复频率${\rm{PR}}{{\rm{F}}_i}(1 \le i \le L)$发射脉冲信号，接收机接收长度为${M_i}(1 \le i \le L)$的回波信号序列${s_i}(n)$；

步骤 2　对第$i(1 \le i \le L)$路回波信号序列${s_i}(n)$做${M_i} = {\rm{PR}}{{\rm{F}}_i} \cdot T\left( {1 \le i \le L} \right)$点DFT并记录峰值位置${\hat k_{{\rm{p}}i}}$；

步骤 3　利用频率余数校正算法计算出频偏${\hat \delta _i}(1 \le i \le L)$，得到新的频率余数估计值

$${\hat r_i} = {\hat k_{{\rm{p}}i}} + {\hat \delta _i}$ $

(33)

步骤 4　把${M_i} = {f_{{\rm{s}}i}} \cdot T = {{\rm PRF}_i} \cdot T\,(1 \le i \le L)$(对应模值组)，${r_i}(1 \le i \le L)$(对应余数组)代入CFRCRT算法，重构多普勒频率的估计值${\hat F}\;$。

4.   仿真实验及分析

实验 1　为了验证本文提出的封闭式鲁棒CRT(CFRCRT)能够精确测量多普勒频率，并且在实时性能能够满足工程要求，做如下实验：假设脉冲多普勒雷达采用3个重复频率参差工作，脉冲重复频率分别是：${f_{{\rm{s}}1}} = {\rm{PR}}{{\rm{F}}_1} = 5\;{\rm{kHz}},{f_{{\rm{s}}2}} = {\rm{PR}}{{\rm{F}}_2} =$$7\;{\rm{kHz}},{f_{{\rm{s}}3}} = {\rm{PR}}{{\rm{F}}_3} = 9\;{\rm{kHz}}$。

根据中国余数定理，最大可测多普勒频率为${f_{{\rm{s}}1}},$${f_{{\rm{s}}2}},{f_{{\rm{s}}3}}$的最小公倍数，即$f_{\max}={\rm{LCM}}({f_{{\rm{s}}1}},{f_{{\rm{s}}2}},{f_{{\rm{s}}3}}) =$$3.15 \times {10^5}\;{\rm{Hz}}$。本次仿真实验不采用频谱校正算法，对序列进行DFT变换，通过搜索峰值位置直接找到频率余数。在$(0,\left. {{f_{\max }}} \right]$范围内任意选取一个频率值，每次实验重复100次的蒙特卡洛频率估计实验，统计其均方根误差(如图4所示)。在$\left( {\left. {0,{f_{\max }}} \right]} \right.$范围内任意选取100个频率，统计其消耗的时间(如图5所示)。

图 4  中国余数定理算法测量频率的精度

下载: 全尺寸图片 幻灯片

图 5  中国余数定理算法的时间复杂度

下载: 全尺寸图片 幻灯片

通过图4，图5可以看出，3种基于CRT算法的测量精度可以比拟。对比于文献[11]封闭式鲁棒CRT算法，本文通过估计共同余数选取最优参考余数构建差分方程组，并通过最优余数重构多普勒频率，大大减小了余数误差的影响，所以提出的封闭式鲁棒CRT算法相比于文献[11]提出的算法测量精度更高一些，本文的封闭式鲁棒CRT和文献[11]所提出的封闭式鲁棒CRT算法估计每个多普勒频率都需要时间${10^{ - 3}}\;{\rm{s}}$，计算复杂度可以比拟，但相比于搜索鲁棒CRT算法计算复杂度大幅度减少，能够满足工程实时性的要求。

实验 2　在信噪比SNR=15 dB的情况下，从$(0,\left. {{f_{\max }}} \right]$范围内任意选取3个频率，分别利用基于apDFT相位差谱校正算法、AM频谱校正算法和Candan频谱校正的CFRCRT方案进行测量。其中$F$是待测的多普勒频率，${{\hat{f}}_{ri}}={{{{\hat{k}}}_{i}}+{{\delta }_{i}}}/{T}\;(1\le i\le 3)$是经过频谱校正算法校正过的最后频率余数，$\Delta F$是多普勒频率估计误差。说明：表1中有2组实验数据，每一组数据分别是按照apDFT相位差频谱校正、AM频谱校正和Candan频谱校正算法的顺序进行测量。

表 1  基于谱校正的中国余数定理的方案测量结果(Hz)

$F$	${\hat f_{{\rm{r}}1}}$	余数理论值	${\hat f_{{\rm{r}}2}}$	余数理论值	${\hat f_{{\rm{r}}3}}$	余数理论值	$\Delta F\;$
${\rm{5}}{\rm{.5122}} \times {\rm{1}}{{\rm{0}}^{\rm{3}}}$	${\rm{512}}{\rm{.60}}$	$512$	${\rm{5512}}{\rm{.65}}$	$5512$	${\rm{5512}}{\rm{.90}}$	$5512$	${\rm{1}}.57 \times {\rm{1}}{{\rm{0}}^{{\rm{ - 1}}}}$
	$512.58$		$5512.57$		$5512.58$		$1.66 \times {10^{ - 1}}$
	$515.04$		$5516.01$		$5516.97$		$3.41 \times {10^0}$
${\rm{3}}{\rm{.1157}} \times {\rm{1}}{{\rm{0}}^{\rm{5}}}$	$1569.27$	$1570$	$3569.69$	$3570$	$5569.80$	$5570$	$2.03 \times {10^{ - 2}}$
	$1568.72$		$3568.70$		$5568.95$		$4.99 \times {10^{ - 2}}$
	$1574.25$		$3566.42$		$5574.39$		$8.42 \times {10^{ - 1}}$

下载: 导出CSV 

| 显示表格

从表1可以看出，信噪比为15 dB的实验环境下，基于Candan频谱校正的CFRCRT方法对多普勒频率测量的均方根误差基本在${10^0} \sim {10^{ - 1}}$数量级，而基于apDFT相位差和AM频偏校正的CFRCRT方法对多普勒频率测量的均方根误差基本都在${10^{ - 1}} \sim {10^{ - 2}}$数量级。相比之下，基于相位差频谱校正的CFRCRT算法精度性能更高些。

实验 3　本次仿真实验将信噪比环境变化的范围设定为${\rm{SNR}} = 0 \sim 20\;{\rm{dB}}$，待测多普勒频率在$\left( {0,\left. {{f_{\max }}} \right]} \right.$范围内随机取一个。每次重复100次蒙特卡洛多普勒频率估计实验，并统计基于apDFT相位差频谱校正、AM频谱校正和Candan频谱校正的CFRCRT算法测量出的多普勒频率与真实值的均方根误差(如图6所示)。从$\left( {0,\left. {{f_{\max }}} \right]} \right.$范围任取100个待测多普勒频率，在信噪比SNB=15 dB情况下，分别利用3种算法进行100次的蒙特卡洛运算时间测量实验(如图7所示)。

图 6  基于频谱校正的封闭式鲁棒CRT测量精度

下载: 全尺寸图片 幻灯片

图 7  基于频谱校正的封闭式鲁棒CRT计算复杂度

下载: 全尺寸图片 幻灯片

从图6可以看出基于apDFT相位差频谱校正的CFRCRT算法在测量精度和抗噪声性能上比基于AM频谱校正的CFRCRT算法稍微好一些，在信噪比15 dB时，精度可以达到${10^{ - 2}}\;{\rm{s}}$数量级。并且从图7可以看出基于apDFT相位差频谱校正的CFRCRT算法在实时性能上明显好于基于AM频谱校正的CFRCRT算法，其计算时间可以达到${10^{ - 3}}\;{\rm{s}}$数量级。

5.   结束语

基于脉冲多普勒雷达体制，利用中国余数定理算法解多普勒频率模糊，引入频谱校正算法有效地改善多普勒频率估计精度，本文分别在仿真实验中分析验证了3种基于频谱校正算法的性能，总结归纳了3种算法的优缺点。在精确度性能上：相位差频谱校正优于AM频谱校正且优于Candan频谱校正；在计算复杂度方面：Candan频谱校正可比拟相位差频谱校正且优于AM频谱校正。但是相位差频谱校正和Candan频谱校正在计算复杂度上基本一致，都可以达到${10^{ - 3}}$ s数量级。从而可以得出一个有意义的结论：基于apDFT相位差频谱校正的封闭式鲁棒中国余数定理多普勒频率估计算法具有比较强的抗干扰性能；算法的计算复杂度不高，具有很强的实时性能，适合工程上的应用。

中国余数定理在解模糊方面是个非常有效的工具，文献[19,20]提出了采用中国余数定理估计多目标参数的算法，但是文献[19,20]提出的算法对余数噪声敏感，不具有鲁棒性。由于多目标参数的余数与目标对应关系事先未知，故基于本文的鲁棒封闭式中国余数定理估计多目标的方法需要进一步详细研究。