CN101957158A - 基于恒包络chirp信号模型的标量脱靶量快速测量方法 - Google Patents

Abstract

本发明涉及雷达领域，具体涉及一种基于无线电技术的靶标坐标系标量脱靶量快速测量方法。包括如下步骤：(1)根据系统要求，确定信号分段时长、步进拟合时长和总拟合时长；(2)发射单频信号，采样接收目标回波，并作预滤波；(3)按照恒包络chirp信号模型，利用分数阶Fourier变换来分段估计回波信号的多普勒频率，并步进拟合出足够时长的多普勒频率变化曲线；(4)根据匀速直线运动模型对得到的多普勒频率变化曲线作最小二乘拟合，估计出标量脱靶量参数。本发明在满足精度要求下缩短了测量耗时，可用于实时测量系统。此外，由于分段时间的延长，改善了积累效果，使得本发明具备较强的鲁棒性，在低信噪比下仍具有较好性能。

Description

基于恒包络chirp信号模型的标量脱靶量快速测量方法 

技术领域

本发明涉及雷达领域，具体涉及一种基于无线电技术的靶标坐标系标量脱靶量快速测量方法。 

背景技术

脱靶量测量是指对射击武器或制导武器遭遇段运动参数的测量，它是评估武器精度、分析武器系统误差因素必不可少的手段。脱靶量测量包括标量脱靶量测量和矢量脱靶量测量。标量脱靶量测量是指测量遭遇过程中弹与靶的相对最小距离，而矢量脱靶量测量不仅要获得标量脱靶量，而且还要测得弹与靶遭遇过程中的相对运动轨迹和相对速度矢量。 

从基本体制上来分，脱靶量测量可分为地面坐标系和靶标坐标系测量。前者以地面坐标为基准，同时测量出靶标与被测目标在地面坐标系中的位置与运动参数，而后换算成脱靶量。由于是以地面坐标系为准，测量误差与靶标和被测目标二者相对于地面测量设备的位置有关。当它们之间的距离较大时，一般测量误差较大，为提高精度，将使设备大为复杂化。而后者是将测量设备置于靶标上，因而测量精度只与被测目标和靶标的相对位置有关，而与二者至地面站的距离无关，故一般易于得到较高的测量精度，设备较为简单。 

从基本原理上来看，脱靶量测量可以借助无线电波、光波、声波等能量形式作为获取信息的媒介。目前，借助无线电波作为信息获取的媒介是脱靶量测量应用最广泛的技术手段。其中基于靶标坐标系的无线电测量方法因其技术难度小，设备简单，且可获得高的测量精度，而成为了国内外重点发展的脱靶量测量技术。 

基于无线电技术的靶标坐标系标量脱靶量测量方法主要有窄脉冲测量法和多普勒测量法。窄脉冲测量的原理是冲击雷达发射极窄的射频脉冲，然后确定被测目标反射回波所在的距离波门随时间的变化，并根据最近的距离波门直接读出标量脱靶量或根据距离波门时间序列经最小二乘拟合后得出标量脱靶量更准确的结果。这种测量方法的优点在于对被测目标与靶标的相对运动特性无限定，可用于高机动目标的测量。其缺点是在低空或地面、海面工作时，易受杂波等影响，且难以实现大脱靶量测量。多普勒测量方法的基本原理是依据弹靶交会过程中回波信号多普勒频率随时间变化的规律取决于被测目标与靶标的相对运动速度、相 对运动加速度、标量脱靶量和脱靶时刻。也就是说，如果能够测出交会过程中被测目标回波的多普勒频率随时间的实际变化曲线，然后基于选定的多普勒频率曲线模型通过最优化方法进行最优拟合，就能估计出被测目标相对于靶标的标量运动速度、脱靶距离等标量脱靶量参数。这种测量方法的优点在于易克服低空或地面、海面工作时杂波等的影响，以及可实现大脱靶量测量。其缺点是需要对弹靶相对运动规律作某种假设(例如：匀速直线运动)，不利于高机动目标测量。 

根据多普勒测量方法的原理，可以知道多普勒频率的测量是其中的关键环节。当前所使用的测频方法都是建立在恒包络恒频的假设前提下，然后采用基于快速Fourier变换(FFT)的谱分析或超分辨的现代谱分析等方法来进行频率测量。前者简写为FFT-MDM，其实时性较好。基于恒包络恒频的做法是对接收信号作分段处理，假定其在每小段内近似为恒频恒包络的谐波信号，对每小段信号作限带降噪后再来估计瞬时频率，得出相应的多普勒频率变化曲线。 

FFT-MDM能够有效降低计算复杂度，满足脱靶量测量系统的实时性要求。但是，由于利用FFT的常规谱分析算法受到瑞利限的限制，因此希望分段的长度比较长，以提高频域分辨率。而脱靶量测量中由于多普勒频率的快速变化而对回波信号的分段通常较短。这就造成了利用恒包络恒频的假设来估计信号频率在脱靶点附近，尤其是在小脱靶量的情况下，会引起较大的测频误差。 

目前现代参数化谱估计算法已被应用于脱靶量的测量中，由于其具有超分辨能力，在信号分段长度较小时仍能较准确地估计出多普勒频率。但是由于参数化谱估计往往都要用到矩阵的奇异值分解，计算量较大，因此，使用该算法的脱靶量测量必须通过相当一段时间的后台处理。所以，只有在对实时性没有要求的应用场合才适合采用该算法。 

在匀速直线运动模型下，被测目标反射回波所产生的多普勒频率为 

$f_d\left(t\right)=\frac{2v}{\mathrm{\λ}}\·\frac{L-\mathrm{vt}}{\sqrt{r^2+{(L-\mathrm{vt})}^2}}$

式中，λ是发射信号波长，v是标量相对运动速度，r是标量脱靶量，L是起始时刻t_j的被测目标位置到脱靶点的距离，t_j一般设定为零时刻。当L-vt_e＝0的时刻t_e，即为脱靶点时刻。通过对多普勒频率变化曲线的分析，可以看出其大致可以分为三部分。前后两部分近似于水平直线，中间一部分在脱靶点附近，多普勒频率变化很大，可近似为一条斜线。三部分曲线连接处斜率变化较大，但段内曲线稳定性较好，因此，如果能估计出各小段曲线的斜率(调频率μ)和频移量(中心频率f_c)，就可以利用各小段斜线有效地拟合出多普勒频率变化曲线。 设第m小段信号的起始时刻为t_m，1，对上式以t_m，1为中心做泰勒级数展开，可得 

$f_{d,m}\left(t\right)=\frac{2v}{\mathrm{\λ}}\·\frac{L-{\mathrm{vt}}_{m,1}}{\sqrt{r^2+{(L-{\mathrm{vt}}_{m,1})}^2}}-\frac{{2v}^2}{\mathrm{\λ}}\·\frac{r^2}{\sqrt{{[r^2+{(L-{\mathrm{vt}}_{m,1})}^2]}^3}}\·(t-t_{m,1})+O(t-t_{m,1})$

式中，O(t-t_m，1)表示泰勒级数展开中(t-t_m，1)的二次方以上的项。从上式可以看出，忽略掉O(t-t_m，1)的影响，则可将该小段目标回波信号近似看成为chirp信号。其起始频率f_cm和调频率(chirp-rate)μ_m分别为： 

$\left\{\begin{array}{c}{f}_{\mathrm{cm}}=\frac{2v}{\mathrm{\λ}}\·\frac{L-{\mathrm{vt}}_{m,1}}{\sqrt{r^2+{(L-{\mathrm{vt}}_{m,1})}^2}}\\ {\mathrm{\μ}}_m=\frac{{2v}^2}{\mathrm{\λ}}\·\frac{r^2}{\sqrt{{[r^2+{(L-{\mathrm{vt}}_{m,1})}^2]}^3}}\end{array}\right.$

也就是说，解决分段时长与测频精度之间的矛盾，且不带来过重的运算负担的可行措施就是：合理地延长分段时长，然后将每小段内信号按照恒包络chirp信号的假设来处理。目前，针对chirp信号的处理工具主要有：解线调技术、多项式相位变换、Wigner-Hough变换(Radon-Wigner变换)、分数阶Fourier变换等。分数阶Fourier变换的定义式如下： 

$X_p\left(u\right)=F_p\left[x\right]\left(u\right)={\∫}_{-\mathrm{\ω}}^{+\∞}x\left(t\right)K_p(t,u)\mathrm{dt}$

式中F_p(·)为p阶分数阶Fourier变换算子， 

$K_p(t,u)=\left\{\begin{array}{cc}\sqrt{(1-j\cot \mathrm{\α})}{e}^{\mathrm{j\π}(t^2\cot \mathrm{\α}-2\mathrm{ut}\csc \mathrm{\α}+u^2\cot \mathrm{\α})}& \mathrm{\α}\≠\mathrm{n\π}\\ \mathrm{\δ}(t-u)& \mathrm{\α}=2\mathrm{n\π}\\ \mathrm{\δ}(t+u)& \mathrm{\α}=(2n\±1)\mathrm{\π}\end{array}\right.$

其中α＝p·π/2。 

分数阶Fourier变换作为传统Fourier变换的广义形式，近年来其相关理论得到了迅速的发展。由于其实质是一种统一的时频变换，与常用二次型时频分布不同的是它不受交叉项困扰，且可以理解为chirp基分解，因此，它十分适于处理多分量chirp信号。此外，分数阶Fourier变换具有比较成熟的快速离散算法(与FFT计算量相当)。因此，本发明为解决已有多普勒测量技术的不足，以分数阶Fourier变换为处理工具，提供一种满足一定精度和实时性要求，同时对噪声不太敏感的标量脱靶量快速测量方法。 

发明内容

本发明是采用数字信号处理的方式来实现的，其具体步骤如下： 

1.确定信号分段时长a、步进拟合时长b和总拟合时长T。 

因为接收信号需要分段进行处理，为了保证各分段信号的调频率和起始频率估计精度，以及提高测量鲁棒性和对弱信号的估计能力(即：脱靶量估计的灵敏度)，要求分段时长在满足恒包络chirp信号的假设前提下尽可能的大；同时，为了保证多普勒频率变化曲线的拟合精度以提高脱靶量估计准确性，又希望拟合步长足够精细，因此，本发明采用分别确定分段时长a与步进拟合时长b的方式，并使得a＞b。其中，信号分段时长a需要满足系统要求的频率分辨率(S_f)和chirp-rate分辨率(S_c)，即 

max(·，·)表示两者中取大者；步进拟合时长b在满足实时性要求的前提下可取到下限，即C_a·T/b≤C，C为事先确定的系统完成一次测量所容许的最大计算量(一般用实数乘和实数加的次数来度量)，C_a为73次N点离散分数阶Fourier变换的计算量，N＝f_s·a+1，f_s为采样频率；而总拟合时长T＝M·b需要保证拟合完的多普勒频率曲线长度超过脱靶点时刻，即T＞t_e-t_j，t_j为起始时刻(一般设定为0)，t_e为脱靶点时刻，M为需要处理的分段信号个数。 

2.发射单频信号，接收被测目标的回波，并对回波信号采样后作预滤波处理，以提高信噪比。 

该预滤波器为常规带限数字滤波器(如：带限FIR滤波器或带限IIR滤波器)，其滤波器带宽可根据实际情况来确定。需要说明的是：如果还需要利用相位信息来进一步估计矢量脱靶量，则预滤波只能采用零相位数字滤波器。 

3.估计回波信号的多普勒频率以得到其变化曲线f_d(t)，t∈[t_j，t_j+t)。由于实际应用中是对离散时间信号进行处理，所以，需要对f_d(t)进行采样，得到f_d(t_k)，t_k∈{t_j，t_j+T_s，...，t_j+T-T_s}，k＝1，…，K，T_s＝1/f_s。 

(1)从起始时刻t_j开始对时长为a的第一个分段[t_j，t_j+a]的信号x₁(n)，n＝1，…，N按照恒包络chirp信号的假设来估计其调频率μ₁和起始频率f_c1，然后以f_c1+μ₁t，t ∈[t_j，t_j+b)拟合出[t_j，t_j+b)段多普勒频率； 

(2)然后对从t_j+b时刻开始的时长为a的第二个分段[t_j+b，t_j+b+a]信号x₂(n)作参数估计，即μ₂和f_c2，并依据μ₂和f_c2的估计值完成第二个时段[t_j+b，t_j+2b)的多普勒频率变化曲线的拟合，即f_c2+μ₂t，t ∈[t_j+b，t_j+2b)； 

(3)依次递推下去，对第m个分段[t_j+(m-1)b，t_j+(m-1)b+a]信号x_m(n)完成μ_m和f_cm的参数估计，并拟合出第m个时段[t_j+(m-1)b，t_j+mb)的多普勒频率变化曲线f_cm+μ_mt，t ∈[t_j+(m-1)b， t_j+mb)，直到多普勒频率的拟合长度达到总拟合时长T。 

(4)以采样频率f_s对拟合出的多普勒频率变化曲线进行采样以得到f_d(t_k)，t_k∈{t_j，t_j+T_s，...，t_j+T-T_s}。 

注：各分段信号调频率μ_m和起始频率f_cm的估计是通过在分数阶Fourier域内对各分段信号进行步进搜索来完成的。设第m分段信号为x_m(n)，其中，n＝1，…，N，表示该分段信号的样本序号，所对应的采样时刻为t_m，n，m＝1，…，M，表示分段序号。本发明中所采用的分数阶Fourier变换快速算法为H.M.Ozaktas等提出的算法，详见“Digital computation of the fractional Fouriertransform”，发表于IEEE Trans.Signal Processing第44卷第9期。 

第一步，首先在[-0.9，1]的分数阶次范围内以0.1为步长，对分段信号x_m(n)逐阶次作分数阶Fourier变换(不需要作0阶分数阶Fourier变换)，得到19个不同阶次的变换结果： 

$X_{{m,p}_l}\left(u_n\right)=F_{p_l}\left[x_m\right]\left(u_n\right),$ p_l-0.9，…，-0.1，0.1，…，1，l＝1，…，19 

然后，在分数阶次p_l和分数阶Fourier域变量u_n所构成的19×N网格上对 

作二维峰值搜索，得到峰值位置的网格坐标序号为(l_p1，n_u1)。则有 

第二步，接下来在[q₁-0.09，q₁+0.09]的分数阶次范围内以0.01为步长，对分段信号x_m(n)逐阶次作分数阶Fourier变换(不需要作q₁阶分数阶Fourier变换)，连同已有的q₁阶结果，则同样会有19个不同阶次的变换结果： 

$X_{{m,p}_l}\left(u_n\right)=F_{p_l}\left[x_m\right]\left(u_n\right),$ p_l＝q₁-0.09，…，q₁，…，q₁+0.09，l＝1，…，19 

同样，对这19个不同阶次的变换结果以分数阶次p_l和分数阶Fourier域变量u_n作二维峰值搜索，得到峰值位置的网格坐标序号为(l_p2，n_u2)。则有 

第三步，将步长缩小为0.001，然后以[q₂-0.009，q₂+0.009]为阶次搜索范围，令p_l＝q₂-0.009，…，q₂，…，q₂+0.009，重复第二步工作，得到峰值位置的网格坐标序号为(l_p3，n_u3)。则有q₃＝pl_p3， 

第四步，将步长缩小为0.0001，然后以[q₃-0.0009，q₃+0.0009]为阶次搜索范围，令p_l＝q₃-0.0009，…，q₃，…，q₃+0.0009，重复第二步工作，得到峰值位置的网格坐标序号为(l_p4， n_u4)。则有 

第五步，得到μ_m和f_cm的估计值： 

$\left\{\begin{array}{c}{\mathrm{\μ}}_m=-\cot \left(2\mathrm{\π}{q}_4\right)\·f_s/a\\ f_{\mathrm{cm}}={\mathrm{\ζ}}_4\·\csc \left(2\mathrm{\π}{q}_4\right)\·\sqrt{f_s/a}-{\mathrm{\μ}}_{m}a/2\end{array}\right.$

综上，最终得到精确到小数点后第4位的分数阶次峰值位置需要进行19+18+18+18＝73次N点离散分数阶Fourier变换，如果需要更高的估计精度，则可以进一步缩小步长来进行搜索。本发明所使用的离散分数阶Fourier变换快速算法已经能够做到与FFT算法的运算量相当，这就保证了本发明能够满足一定的实时性要求。 

4.根据设定的匀速直线运动模型对上一步所得到的多普勒频率变化曲线作最小二乘拟合，便可以估计出相应的标量脱靶量参数 

和 

其中， 是标量相对运动速度估计值， 

是起始时刻t_b的被测目标位置到脱靶点的距离估计值， 是标量脱靶量估计值。在连续时间情况下，多普勒周期数随时间变化情况为： 

$\mathrm{\φ}\left(t\right)={\∫}_0^t{f}_d\left(t\right)\mathrm{dt}$

经过推导后，可得： 

${\left[\frac{\mathrm{\λ\φ}\left(t\right)}{2}\right]}^2=\left[\begin{array}{ccc}\mathrm{\λ\φ}\left(t\right)& t& t^2\end{array}\right]\·\left[\begin{array}{c}R\\ g\\ h\end{array}\right]$

式中， 

g＝-2Lv，h＝v²。既然处理的是离散时间信号，则用求和代替积分，得到多普勒周期数 

$\mathrm{\φ}\left(t_k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{f}_d\left(t_i\right)\·T_s$

因为 

$\left[\begin{array}{c}{\left[\frac{\mathrm{\λ\φ}\left(t_1\right)}{2}\right]}^2\\ {\left[\frac{\mathrm{\λ\φ}\left(t_2\right)}{2}\right]}^2\\ .\\ .\\ .\\ {\left[\frac{\mathrm{\λ\φ}\left(t_K\right)}{2}\right]}^2\end{array}\right]=\left[\begin{array}{ccc}\mathrm{\λ\φ}\left(t_1\right)& t_1& {t_1}^2\\ \mathrm{\λ\φ}\left(t_2\right)& t_2& {t_2}^2\\ .& .& .\\ .& .& .\\ .& .& .\\ \mathrm{\λ\φ}\left(t_K\right)& t_K& {t_K}^2\end{array}\right]\·\left[\begin{array}{c}R\\ g\\ h\end{array}\right]$

利用最小二乘技术解上式可以获得估计值 进一步可以得到： 

$\left\{\begin{array}{c}\hat{v}=\sqrt{\hat{h}}\\ \hat{L}=-\hat{g}/2\hat{v}\\ \hat{r}=\sqrt{{\hat{R}}^2-{\hat{L}}^2}\end{array}\right.$

附图说明

图1是本发明所提供方法的流程图。 

图2是被测目标相对测量天线的运动模型图。 

图3是四种不同初始条件下的多普勒频率变化曲线图。 

图4是本发明中估计各分段信号调频率和起始频率的流程图。 

图5是两种不同初始条件下本发明与FFT-MDM方法的仿真比较。 

图6是0dB下本发明的标量脱靶量估计误差。 

具体实施方式

下面结合附图及实施例对本发明做详细说明： 

本发明的一种基于无线电技术的靶标坐标系标量脱靶量快速测量方法的总体流程如图1所示，包括以下步骤： 

设发射的雷达信号频率为2.9吉赫兹(即波长约为3/29米)，基带信号采样频率f_s＝78.125千赫兹。假定此次测量的标量脱靶量参数为：r＝36.58米，v＝1000米/秒，L＝300米。其匀速直线运动模型如图2所示，在该模型下，被测目标反射回波所产生的多普勒频率为 

$f_d\left(t\right)=\frac{2000}{3/29}\·\frac{300-1000t}{\sqrt{{36.58}^2+{(300-1000t)}^2}}$

1.确定信号分段时长a、步进拟合时长b和总拟合时长T。 

令t_j＝0，根据背景技术所述内容，可以知道脱靶点时刻t_e＝0.3秒，则时间t从t_j到t_e，分段信号的起始频率从19.191千赫兹下降到0赫兹，chirp-rate从937赫兹/秒上升到216.044兆赫兹/秒(如图3所示)。因此，不妨令频率分辨率S_f＝80赫兹，chirp-rate分辨率S_c＝50千赫兹/秒，即a必须不小于max(11.07毫秒，11.79毫秒)。既然采样间隔T_s＝1/f_s＝12.8纳秒，则意味着分段信号的样本数不能小于max(11.07/0.0128，11.79/0.0128)，即分段信号样本数必须大于922。为便于计算，不妨取分段信号样本数为1025，则a＝13.1072毫秒。 

步进拟合时长b的确定需要结合具体硬件情况。既然本发明所采用的分数阶Fourier变换离散算法的复杂度与FFT相当，C_a为73次1025点离散分数阶Fourier变换的计算量，那么不妨令b＝6.5536毫秒，即513个样本点的时间间隔。 

总拟合时长T需要保证拟合完的多普勒频率曲线长度超过脱靶点时刻，即T＞t_e-t_j＝0.3秒。既然0.3秒包含约23438个样本，那么不妨取T包含25600个样本，也就是，T＝0.3768秒。那么需要处理的分段信号个数M＝50。 

2.发射单频信号，接收被测目标的回波，并对回波信号采样后作预滤波处理，以提高信噪比。 

对接收到的回波信号降频到基带，然后以采样频率f_s＝78.125千赫兹对基带回波信号进行采样。采样后的信号按照分段时长(13.1072毫秒)和步进拟合时长(6.5536毫秒)分段，共得到50个分段信号，然后对每个分段信号都利用零相位FIR滤波器作降噪处理。该滤波器首先采用1025点汉明窗平滑截断误差，然后再用1025点FIR低通滤波器滤除掉20千赫兹以上的频率成分。 

3.估计回波信号的多普勒频率f_d(t_k)，t_k∈{0，0.0128毫秒，...，376.8-0.0128毫秒}，k＝1，…，25600。 

(1)从零时刻开始对时长为13.1072毫秒的第一个分段[0，13.1072毫秒]的信号x₁(n)，n＝1，…，1025按照恒包络chirp信号的假设来估计其调频率μ₁和起始频率f_c1(如图4所示)，然后以f_c1+μ₁t，t∈[0，6.5536毫秒)拟合出[0，6.5536毫秒)时段的多普勒频率； 

(2)然后对从6.5536毫秒时刻开始的时长为13.1072毫秒的第二个分段[6.5536毫秒，19.6608毫秒]信号x₂(n)作参数估计(如图4所示)，即μ₂和f_c2，并依据μ₂和f_c2的估计值完成第二个时段[6.5536毫秒，13.1072毫秒)的多普勒频率变化曲线的拟合，即f_c2+μ₂t，t ∈[6.5536毫秒，13.1072毫秒)； 

(3)依次递推下去，对第m个分段[(m-1)×6.5536毫秒，(m-1)×6.5536+13.1072毫秒]信号x_m(n)完成μ_m和f_cm的参数估计(如图4所示)，并拟合出第m个时段[(m-1)×6.5536毫秒，m×6.5536毫秒)的多普勒频率变化曲线f_cm+μ_mt，t∈[(m-1)×6.5536毫秒，m×6.5536毫秒)，直到多普勒频率的拟合长度达到总拟合时长0.3768秒。 

(4)以采样频率78.125千赫兹对拟合出的多普勒频率变化曲线进行采样以得到f_d(t_k)，t_k∈{0，0.0128毫秒，...，376.8-0.0128毫秒}。 

接下来以第45个分段[288.4毫秒，301.4656毫秒]的信号x₄₅(n)为例(该信号叠加有10分 贝的高斯白噪声)，说明参数估计过程(如图4所示)。 

第一步，首先在[-0.9，1]的分数阶次范围内以0.1为步长，对分段信号x₄₉(n)逐阶次作分数阶Fourier变换(不需要作0阶分数阶Fourier变换)，得到19个不同阶次的变换结果。然后，在分数阶次p_l和分数阶Fourier域变量u_n所构成的19×1025网格上对变换结果作二维峰值搜索，得到峰值位置的网格坐标序号为(18，575)。则有q₁＝0.9，ζ₁＝1.9375。 

第二步，接下来在[0.9-0.09，0.9+0.09]的分数阶次范围内以0.01为步长，对分段信号x₄₉(n)逐阶次作分数阶Fourier变换(不需要作0.9阶分数阶Fourier变换)，连同已有的0.9阶结果，则同样会有19个不同阶次的变换结果。同样，对这19个不同阶次的变换结果以分数阶次p_l和分数阶Fourier域变量u_n作二维峰值搜索，得到峰值位置的网格坐标序号为(15，549)。则有q₂＝0.95，ζ₂＝1.125。 

第三步，将步长缩小为0.001，然后以[0.95-0.009，0.95+0.009]为阶次搜索范围，令p_l＝0.95-0.009，…，0.95，…，0.95+0.009，重复第二步工作，得到峰值位置的网格坐标序号为(4，548)。则有q₃＝0.944，ζ₃＝1.0938。 

第四步，将步长缩小为0.0001，然后以[0.944-0.0009，0.944+0.0009]为阶次搜索范围，令p_l＝0.944-0.0009，…，0.944，…，0.944+0.0009，重复第二步工作，得到峰值位置的网格坐标序号为(11，548)。则有q₄＝0.9441，ζ₄＝1.0938。 

第五步，得到μ₄₅和f_c45的估计值： 

4.根据设定的匀速直线运动模型对上一步所得到的多普勒频率变化曲线f_d(t_k)作最小二乘拟合，便可以估计出相应的标量脱靶量参数 

和 

其中， 

是标量相对运动速度估计值， 

是起始时刻t_j的被测目标位置到脱靶点的距离估计值， 

是标量脱靶量估计值。 

首先利用下式求得多普勒周期数 

$\mathrm{\φ}\left(t_k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{f}_d\left(t_i\right)\·T_s,$ t_k∈{0，0.0128ms，...，376.8-0.0128ms}，k＝1，…，25600。 

因为 

$\left[\begin{array}{c}{\left[\frac{\mathrm{\λ\φ}\left(t_1\right)}{2}\right]}^2\\ {\left[\frac{\mathrm{\λ\φ}\left(t_2\right)}{2}\right]}^2\\ .\\ .\\ .\\ {\left[\frac{\mathrm{\λ\φ}\left(t_{25600}\right)}{2}\right]}^2\end{array}\right]=\left[\begin{array}{ccc}\mathrm{\λ\φ}\left(t_1\right)& t_1& {t_1}^2\\ \mathrm{\λ\φ}\left(t_2\right)& t_2& {t_2}^2\\ .& .& .\\ .& .& .\\ .& .& .\\ \mathrm{\λ\φ}\left(t_{25600}\right)& t_{25600}& {t_{25600}}^2\end{array}\right]\·\left[\begin{array}{c}R\\ g\\ h\end{array}\right]$

将上式写成矩阵的形式如下： 

Y＝U·V 

式中，Y为25600×1的列向量，U为25600×3的矩阵，V为3×1的列向量。利用最小二乘技术解上式可以获得估计值 如下： 

$\hat{V}={(U^T\·U)}^{-1}\·U^T\·Y$

式中， 

进一步可以得到： 

如图5所示，(1)本发明估计精度要明显优于FFT-MDM算法；(2)本发明具有较强的鲁棒性，即，随着信噪比的下降，脱靶量估计精度并没有明显的下降；(3)在多普勒频率变化较为剧烈的情况下，本发明仍具有较好的鲁棒性，而FFT-MDM方法则随着信噪比的下降其误差明显增大。 

如图6所示，本发明在低信噪比情况下对四组不同参数仍都具有较好的估计精度。 

Claims (5)

1.基于恒包络chirp信号模型的标量脱靶量快速测量方法，其特征在于：本发明包括以下四个步骤，其中：

(1)根据灵敏度、精度、单次测量耗时以及鲁棒性的要求，确定信号分段时长a、步进拟合时长b和总拟合时长T；

(2)发射单频信号，检测并跟踪被测目标的回波，在对回波信号采样后作预滤波处理，以提高信噪比；

(3)按照恒包络chirp信号的假设，利用分数阶Fourier变换来分段估计回波信号的多普勒频率，并步进拟合出足够时长的多普勒频率变化曲线；

(4)根据匀速直线运动模型对上一步所得到的多普勒频率变化曲线作最小二乘拟合，估计出相应的标量脱靶量参数。

2.根据权利要求1所述的基于恒包络chirp信号模型的标量脱靶量快速测量方法，其特征在于：步骤(1)中分别确定分段时长a与步进拟合时长b，并使得a＞b。其中，信号分段时长a需要满足系统要求的频率分辨率(S_f)和chirp-rate分辨率(S_c)，即

max(·，·)表示两者中取大者；步进拟合时长b在满足实时性要求的前提下可取到下限，即C_a·T/b≤C，C为事先确定的系统完成一次测量所容许的最大计算量(一般用实数乘和实数加的次数来度量)，C_a为73次N点离散分数阶Fourier变换的计算量，N＝f_s·a+1，f_s为采样频率。

3.根据权利要求1所述的基于恒包络chirp信号模型的标量脱靶量快速测量方法，其特征在于，步骤(1)中所确定的总拟合时长T＝M·b需要保证拟合完的多普勒频率曲线长度超过脱靶点时刻，即T＞t_e-t_j，t_j为起始时刻(一般设定为0)，t_e为脱靶点时刻，M为需要处理的分段信号个数。

4.根据权利要求1所述的基于恒包络chirp信号模型的标量脱靶量快速测量方法，其特征在于：步骤(3)中的多普勒频率变化曲线依如下步骤得到：

①从起始时刻t_j开始对时长为a的第一个分段[t_j，t_j+a]的信号x₁(n)，n＝1，…，N按照恒包络chirp信号的假设来估计其调频率μ₁和起始频率f_c1，然后以f_c1+μ₁t，t ∈[t_j，t_j+b)拟合出[t_j，t_j+b)段多普勒频率；

②然后对从t_j+b时刻开始的时长为a的第二个分段[t_j+b，t_j+b+a]信号x₂(n)作参数估计，即μ₂和f_c2，并依据μ₂和f_c2的估计值完成第二个时段[t_j+b，t_j+2b)的多普勒频率变化曲线的拟合，即f_c2+μ₂t，t ∈[t_j+b，t_j+2b)；

③依次递推下去，对第m个分段[t_j+(m-1)b，t_j+(m-1)b+a]信号x_m(n)完成μ_m和f_cm的参数估计，并拟合出第m个时段[t_j+(m-1)b，t_j+mb)的多普勒频率变化曲线f_cm+μ_mt，t ∈[t_j+(m-1)b，t_j+mb)，直到多普勒频率的拟合长度达到总拟合时长T。

④以采样频率f_s对拟合出的多普勒频率变化曲线进行采样以得到f_d(t_k)，t_k∈{t_j，t_j+T_s，...，t_j+T-T_s}。

5.根据权利要求1所述的基于恒包络chirp信号模型的标量脱靶量快速测量方法，其特征在于：步骤(3)中各分段信号调频率μ_m和起始频率f_cm依如下步骤得到：

设第m分段信号为x_m(n)，n＝1，…，N，表示该分段信号的样本序号，所对应的采样时刻为t_m，n。

第一步，首先在[-0.9，1]的分数阶次范围内以0.1为步长，对分段信号x_m(n)逐阶次作分数阶Fourier变换(不需要作0阶分数阶Fourier变换)，得到19个不同阶次的变换结果：

$X_{m,p_l}\left(u_n\right)=F_{p_l}\left[x_m\right]\left(u_n\right),$ p_l-0.9，…，-0.1，0.1，…，1，l＝1，…，19

所采用的分数阶Fourier变换离散算法为H.M.Ozaktas等提出的算法，详见“Digitalcomputation of the fractional Fourier transform”，发表于IEEE Transactions on SignalProcessing第44卷第9期。

然后，在分数阶次p_l和分数阶Fourier域变量u_n所构成的19×N网格上对

作二维峰值搜索，得到峰值位置的网格坐标序号为(l_p1，n_u1)。则有

第二步，接下来在[q₁-0.09，q₁+0.09]的分数阶次范围内以0.01为步长，对分段信号x_m(n)逐阶次作分数阶Fourier变换(不需要作q₁阶分数阶Fourier变换)，连同已有的q₁阶结果，则同样会有19个不同阶次的变换结果：

$X_{m,p_l}\left(u_n\right)=F_{p_l}\left[x_m\right]\left(u_n\right),$ p_l＝q₁-0.09，…，q₁，…，q₁+0.09，l＝1，…，19

同样，对这19个不同阶次的变换结果以分数阶次p_l和分数阶Fourier域变量u_n作二维峰值搜索，得到峰值位置的网格坐标序号为(l_p2，n_u2)。则有

第三步，将步长缩小为0.001，然后以[q₂-0.009，q₂+0.009]为阶次搜索范围，令p_l＝q₂-0.009，…，q₂，…，q₂+0.009，重复第二步工作，得到峰值位置的网格坐标序号为(l_p3，n_u3)。则有

第四步，将步长缩小为0.0001，然后以[q₃-0.0009，q₃+0.0009]为阶次搜索范围，令p_l＝q₃-0.0009，…，q₃，…，q₃+0.0009，重复第二步工作，得到峰值位置的网格坐标序号为(l_p4，n_u4)。则有

第五步，得到μ_m和f_cm的估计值：

$\left\{\begin{array}{c}{\mathrm{\μ}}_m=-\cot \left({2\mathrm{\πq}}_4\right)\·f_s/a\\ f_{\mathrm{cm}}={\mathrm{\ζ}}_4\·\csc \left({2\mathrm{\πq}}_4\right)\·\sqrt{f_s/a}-{\mathrm{\μ}}_{m}a/2\end{array}\right.$

Images