CN105259553A - 基于距离-瞬时多普勒像的微动目标散射点航迹关联方法 - Google Patents

Abstract

本发明公开了一种基于距离-瞬时多普勒像的微动目标散射点航迹关联方法，其步骤包括：(1)接收逆合成孔径雷达回波；(2)逆合成孔径雷达回波预处理；(3)通过距离-瞬时多普勒像获得微动目标散射点的二维坐标；(4)二维航迹关联；(5)采用求根多重信号分类Root-MUSIC方法修正航迹矩阵；(6)获得微动目标散射点的航迹关联结果。本发明将散射点关联的域从距离-慢时间域转换到可分性更强的慢时间-距离-瞬时多普勒域，增强了微动目标散射点的可分性，对自旋、进动、章动等复杂运动形式具有普适性、航迹关联精度高，关联误差小等优点，能够为后续高分辨成像和目标特征提取提供有力保障。

Description

基于距离-瞬时多普勒像的微动目标散射点航迹关联方法

技术领域

本发明属于信号处理技术领域，更进一步涉及微动目标成像与特征提取领域中的基于距离-瞬时多普勒像的微动目标散射点航迹关联方法。本发明可以准确估计出弹道目标、空间碎片等微动目标散射点在距离-慢时间域的航迹，对自旋、进动、章动等复杂微动形式下的弹道目标、空间碎片散射点航迹关联精度高、关联误差小，能够为后续微动目标的高分辨成像和微动目标特征提取提供有力保障。

背景技术

当采用宽带雷达对微动目标进行高分辨观测时，微动目标的高分辨一维距离像(HRRP：high-resolutionrangeprofiles)反映出微动目标散射点分布沿雷达视线方向的投影。在距离-慢时间平面上，每个散射点的瞬时斜距随慢时间的变化曲线构成该散射点的航迹。

王昕等人在其发表的论文“基于二维ISAR图像序列的雷达目标三维重建方法”(《电子与信息学报》2013，35(10):2475-2480)中提出了一种基于卡尔曼滤波和最近邻准则的二维散射中心航迹关联方法。该方法的具体步骤是，首先从单部雷达获得的ISAR序列中提取目标散射中心二维位置，然后采用卡尔曼滤波和最近邻准则实现不同姿态下二维散射中心的关联。该方法可以解决多个不同姿态下的二维散射中心的关联问题，但是，该方法仍然存在不足之处是，当不同散射中心航迹存在较多交叉点时，卡尔曼滤波方法在交叉点附近关联效果差。

张颖康等人在其发表的论文“基于散射中心关联的三维成像算法”(《系统工程与电子技术》2011，33(9):1988-1994)中提出了一种基于随机抽样一致性的散射点关联方法，该方法的具体步骤是，判断各关联样本反投影内点数量，然后剔除错误的关联样本，完成散射点的航迹关联。该方法可以解决目标散射中心和姿态数目较多的情况下的散射点航迹关联问题，但是，该方法仍然存在不足之处是，在穷举所有姿态下一维散射中心的任意关联组合来提取真实三维散射中心位置时，算法复杂度高。

发明内容

本发明的目的在于克服上述现有技术的不足，提出一种基于距离-瞬时多普勒像的微动目标散射点航迹关联方法。该方法克服了基于卡尔曼滤波和最近邻准则的二维散射中心航迹关联方法，在航迹交叉点附近关联效果差的缺点；弥补了基于随机抽样一致性的散射点关联方法算法复杂度高的不足。

实现本发明的基本思路是：首先，通过距离-瞬时多普勒像的方法实现微动目标散射点航迹的粗关联，之后采用求根多重信号分类Root-MUSIC的方法修正粗关联航迹矩阵实现更加精确的细关联。

本发明的具体步骤如下：

(1)接收微动目标逆合成孔径雷达回波；

(2)逆合成孔径雷达回波预处理：

(2a)将逆合成孔径雷达到场景中心距离作为参考距离，将与逆合成孔径雷达发射信号载频、调频率相同、距离为参考距离的线性调频信号作为参考信号；

(2b)将参考信号取共轭后与接收的回波相乘，得到解线频调处理后的信号；

(3)获得微动目标散射点的二维坐标：

(3a)将解线频调处理后的信号在距离维做一维傅里叶变换，得到目标距离-慢时间域的逆合成孔径雷达回波；

(3b)将距离-慢时间域存在目标回波的每个距离单元做短时傅里叶变换，得到每个距离单元的二维时频图像；

(3c)将每个距离单元对应的二维时频图像沿距离维堆叠成距离-瞬时多普勒-慢时间三维矩阵；

(3d)将获得微动目标散射点二维坐标时的慢时间单元序号β初始化为1；

(3e)取三维矩阵在第β个慢时间单元对应的二维矩阵切片，得到第β个慢时间单元的微动目标距离-瞬时多普勒像；

(3f)用第β个慢时间单元的微动目标距离-瞬时多普勒像中每一个像素值除以第β个慢时间单元的微动目标距离-瞬时多普勒像中最大像素值，得到第β个慢时间单元归一化后的微动目标距离-瞬时多普勒像；

(3g)设定门限值th＝0.3；

(3h)将第β个慢时间单元归一化后的微动目标距离-瞬时多普勒像中像素值大于等于门限值的像素点赋值为1，将第β个慢时间单元归一化后的微动目标距离-瞬时多普勒像中像素值小于门限值的像素点赋值为0，得到二值图像；

(3i)将二值图像中像素值为0的像素点值赋值为1，像素值为1的像素点赋值为0，得到二值图像的补图像；

(3j)计算二值图像的补图像中每个像素点到其最近非零像素点的欧氏距离，得到二值图像补图像的距离变换矩阵；

(3k)将二值图像补图像的距离变换矩阵中的每个值取反，得到灰度图像；

(3l)采用分水岭图像分割方法，对灰度图像做图像分割，得到分割后的图像；

(3m)将分割后的图像中像素值为1和0的像素点均赋值为0，将其他像素点赋值为255，得到分割后的二值图像；

(3n)采用四邻域标记方法，得到分割后的二值图像的连通域，将分割后二值图像的每个连通域作为分割后二值图像中每个散射点对应的区域；

(3o)按照下式，计算分割后二值图像中每个散射点对应区域的区域重心横坐标值和纵坐标值；

$x_c=\frac{\underset{a}{\Σ}\underset{b}{\Σ}af(a,b)}{\underset{a}{\Σ}\underset{b}{\Σ}f(a,b)},y_c=\frac{\underset{a}{\Σ}\underset{b}{\Σ}bf(a,b)}{\underset{a}{\Σ}\underset{b}{\Σ}f(a,b)}$

其中，x_c表示分割后二值图像中每个散射点对应区域重心的横坐标值，y_c表示分割后二值图像中每个散射点对应区域重心的纵坐标值，a表示分割后二值图像中每个散射点对应区域中像素点的横坐标值，b表示分割后二值图像中每个散射点对应区域中像素点的纵坐标值，Σ表示求和操作，f(·)表示分割后二值图像中每个散射点对应的区域中像素点的像素值，f(a,b)表示分割后二值图像中每个散射点对应区域中横纵坐标分别为a、b的像素点的像素值；

(3p)将分割后二值图像中每个散射点对应区域的区域重心横坐标值和纵坐标值，作为分割后二值图像中每个散射点二维坐标的横坐标值和纵坐标值；

(3q)将分割后二值图像中每个散射点二维坐标的横坐标值和纵坐标值，作为微动目标距离-瞬时多普勒图像中每个散射点二维坐标的横坐标值和纵坐标值，并记录该横坐标值和纵坐标值；

(3r)判断是否获得所有时刻微动目标距离-瞬时多普勒图像中散射点的二维坐标，若是，执行步骤(4)，否则，令β＝β+1，执行步骤(3e)，其中，β表示获得微动目标散射点二维坐标时的慢时间单元序号；

(4)二维航迹关联：

(4a)生成列数为慢时间个数，行数为微动目标散射点个数的空矩阵，将该空矩阵作为航迹矩阵，将航迹矩阵中每个元素定义为二维向量；

(4b)将完成二维航迹关联时的慢时间单元序号l初始化为1，将航迹矩阵列编号i初始化为1，将散射点序号α初始化为1；

(4c)将航迹矩阵中的第一列元素值分别赋为微动目标各散射点第l个慢时间单元距离-瞬时多普勒图像中的二维坐标值；

(4d)令l＝l+1，其中l表示完成二维航迹关联时的慢时间单元序号，计算航迹矩阵第i列第α个散射点的二维坐标值与第l个慢时间单元距离-瞬时多普勒图像中所有散射点二维坐标值之间的欧氏距离；

(4f)将欧氏距离最小值所对应的第l个慢时间单元距离瞬时多普勒图像中散射点的二维坐标写入航迹矩阵第i+1列，第α行，完成第l个慢时间单元微动目标第α个散射点的航迹关联；

(4g)判断是否完成第l个慢时间单元微动目标所有散射点的航迹关联，若是，执行步骤(4h)，否则，令散射点序号α＝α+1，执行步骤(4f)；

(4h)判断是否完成所有时刻微动目标所有散射点的航迹关联，若是，得到航迹矩阵，执行步骤(5)；否则，航迹矩阵的列编号i＝i+1，执行步骤(4d)；

(4i)保留航迹矩阵每个元素的第一维坐标，删除每个元素的第二维坐标，得到更新后的航迹矩阵；

(5)修正航迹矩阵：

(5a)将修正航迹矩阵时的慢时间单元序号p初始化为1；

(5b)采用求根多重信号分类Root-MUSIC方法，对第p个慢时间单元解线频调处理后的信号做角频率估计，得到第p个慢时间单元解线频调处理后信号的n个角频率，其中，n表示微动目标散射点个数；

(5c)按照下式，计算解线调频处理后信号的瞬时斜距：

$R_p=\frac{\mathrm{\ω}\·c\·F_s}{-4\mathrm{\π}\mathrm{\γ}}$

其中，R_p表示第p个慢时间单元解线调频处理后信号的瞬时斜距，ω表示第p个慢时间单元解线调频处理后信号的角频率，c表示光速，γ表示调频率，F_s表示逆合成孔径雷达发射信号的距离采样频率，其取值范围为F_s≥B，B表示逆合成孔径雷达发射信号的带宽；

(5d)分别计算第p个慢时间单元每个瞬时斜距与航迹矩阵第p列中每个散射点对应瞬时斜距的欧氏距离，用欧氏距离最小值对应的瞬时斜距值代替航迹矩阵第p列中对应散射点的瞬时斜距，得到第p列修正后的航迹矩阵；

(5e)判断是否完成所有慢时间单元航迹矩阵的修正，若是，执行步骤(6)，否则，令p＝p+1，执行步骤(5b)；

(6)获得微动目标散射点的航迹关联结果。

与现有技术相比，本发明具有以下优点。

第一，本发明通过距离-瞬时多普勒像获得微动目标散射点的二维坐标的方法，将散射点关联的域从传统的距离-慢时间域转换到可分性更强的慢时间-距离-瞬时多普勒域，该方法克服了基于卡尔曼滤波和最近邻准则的二维散射中心航迹关联方法，在航迹交叉点附近关联效果差的缺点，使得本发明具有了航迹关联误差小，航迹关联精度高的优点。

第二，本发明通过采用求根多重信号分类Root-MUSIC方法，修正航迹矩阵，弥补了基于随机抽样一致性的散射点关联方法算法复杂度高的不足，使得本发明具有运算复杂度低，算法简单易行的优点。

附图说明

图1为本发明的流程图；

图2为本发明的仿真图。

具体实施方式

下面结合附图，对本发明作进一步的详细描述。

下面结合附图1，对本发明具体实施方式作进一步的详细描述。

步骤1，接收微动目标逆合成孔径雷达回波。

步骤2，逆合成孔径雷达回波预处理。

将逆合成孔径雷达到场景中心距离作为参考距离，将与逆合成孔径雷达发射信号载频、调频率相同、距离为参考距离的线性调频信号作为参考信号，其中参考信号如下：

$S(\hat{t},t_m)=rect\left(\frac{\hat{t}-2R_{ref}/c}{T_p}\right)e^{j2\mathrm{\π}\left(f_c\right(\hat{t}-2R_{ref}/c)+\frac{1}{2}\mathrm{\γ}{\left(\hat{t}-2R_{ref}/c\right)}^2)}$

其中，表示参考信号，表示快时间，t_m表示慢时间，rect(·)表示矩形窗函数， $rect\left(u\right)=\left\{\begin{array}{cc}1& \left|u\right|\≤\frac{1}{2}\\ 0& \left|u\right|>\frac{1}{2}\end{array}\right.,$ R_ref表示参考距离，T_p表示参考信号的脉宽，j表示虚数单位，f_c表示载频，c表示光速，γ表示调频率。

将参考信号取共轭后与接收的回波相乘，得到解线频调处理后的信号：

$s_{P2}(f^{\′},t_m)={\mathrm{\σ}}_P{a}_B\left(f^{\′}\right)\exp (-j\frac{4\mathrm{\π}}{c}(f_c+f^{\′}\left){\mathrm{\ΔR}}_P\right(t_m\left)\right)$

其中，exp(·)表示取以e为底的指数操作，j表示虚数单位，σ_p表示散射点P的回波幅度，c表示光速，f_c表示载频，表示多普勒频率，γ表示调频率，表示快时间，R_ref表示参考距离，t_m表示慢时间，a_B(f′)表示矩形窗，当|f′|≤B/2时，a_B(f′)＝1，否则，a_B(f′)＝0，B表示逆合成孔径雷达发射信号带宽， ${\mathrm{\ΔR}}_P(\hat{t},t_m)=R_P(\hat{t},t_m)-R_{ref}\left(t_m\right),$ 表示散射点P到场景中心的距离。

步骤3，获得微动目标散射点的二维坐标。

(3a)将解线频调处理后的信号在距离维做一维傅里叶变换，得到目标距离-慢时间域的逆合成孔径雷达回波；

(3b)将距离-慢时间域存在目标回波的每个距离单元做短时傅里叶变换，得到每个距离单元的二维时频图像，其中短时傅里叶变换的具体步骤如下：

STFT_s(t′,f)＝∫s·w(t-t′)exp{-j2πft}dt

其中，STFT(·)表示短时傅里叶变换操作，STFT_s(·)表示对存在目标回波的距离单元信号做短时傅里叶变换操作，t′表示窗函数中心所在的时间点，f表示瞬时多普勒，∫表示积分操作，s表示解线频调处理后信号中存在目标回波的距离单元信号，w(·)表示窗函数，通常选高斯窗或者汉明窗，窗长为N_a表示慢时间单元的总数，w(t-t′)表示沿时间轴平移后的窗函数，exp(·)表示取以e为底的指数操作，j表示虚数单位，t表示慢时间。

(3c)将每个距离单元对应的二维时频图像沿距离维堆叠成距离-瞬时多普勒-慢时间三维矩阵；

(3d)将获得微动目标散射点二维坐标时的慢时间单元序号β初始化为1；

(3e)取三维矩阵在第β个慢时间单元对应的二维矩阵切片，得到第β个慢时间单元的微动目标距离-瞬时多普勒像；

(3f)用第β个慢时间单元的微动目标距离-瞬时多普勒像中每一个像素值除以第β个慢时间单元的微动目标距离-瞬时多普勒像中最大像素值，得到第β个慢时间单元归一化后的微动目标距离-瞬时多普勒像；

(3g)设定门限值th＝0.3；

(3h)将第β个慢时间单元归一化后的微动目标距离-瞬时多普勒像中像素值大于等于门限值的像素点赋值为1，将第β个慢时间单元归一化后的微动目标距离-瞬时多普勒像中像素值小于门限值的像素点赋值为0，得到二值图像；

(3i)将二值图像中像素值为0的像素点值赋值为1，像素值为1的像素点赋值为0，得到二值图像的补图像；

(3j)计算二值图像的补图像中每个像素点到其最近非零像素点的欧氏距离，得到二值图像补图像的距离变换矩阵；

(3k)将二值图像补图像的距离变换矩阵中的每个值取反，得到灰度图像；

(3l)采用分水岭图像分割方法，对灰度图像做图像分割，得到分割后的图像。

分水岭图像分割方法的具体步骤如下：

第一步：对灰度图像中的各像素点按灰度值大小从小到大排序，将所有像素点中相同灰度值的像素点作为同一个层级；

第二步：初始化层级编号k，令k＝1；

第三步：处理第k个层级所有的像素点，如果该像素点的邻域已经被标识属于某一个区域，则将该像素点加入一个先进先出队列；

第四步：判断先进先出队列是否为空，若非空，则将先进先出队列的元素编号n初始化为1，执行本步骤的第五步；否则，执行本步骤的第七步；

第五步：取出先进先出队列中的第n个元素，扫描第n个元素的邻域像素，如果第n个元素及其邻域像素的灰度值属于同一层级，则用邻域像素的标识刷新先进先出队列的第n个元素像素的标识；

第六步：令先进先出队列的元素编号n＝n+1，重复本步骤的第五步，直到先进先出队列为空；

第七步：扫描第k个层级所有的像素点，判断第k个层级像素点是否都被标识，若有像素点未被标识，则将该像素点的灰度值加1，若第k个层级像素点都被标识，执行本步骤的第八步；

第八步：判断是否处理完所有层级的像素点，若是，得到分割后的图像，否则，执行本步骤的第二步。

(3m)将分割后的图像中像素值为1和0的像素点均赋值为0，将其它像素点赋值为255，得到分割后的二值图像；

(3n)采用四邻域标记方法，得到分割后的二值图像的连通域，将分割后二值图像的每个连通域作为分割后二值图像中每个散射点对应的区域。

四邻域标记法的具体步骤如下：

第一步：将分割后的二值图像的像素点编号m初始化为K+2，其中K表示距离-瞬时多普勒像的列数；

第二步：判断m是否为K的整数倍，若是，则m＝m+1，执行第三步，否则直接执行第三步，其中K表示距离-瞬时多普勒像的列数；

第三步：若第m个像素值的左边、上边的像素值与第m个像素值都不相同，则将该像素值所在的区域标记与为第m个像素值的左边、上边的像素值所在的区域都不相同的一个区域，执行第七步；

第四步：若第m个像素值的左边的像素值与第m个像素值相同，第m个像素点的上边的像素值与第m个像素值不相同，则将第m个像素点标记为与第m个像素点的左边的像素值所在的区域相同的一个区域，执行第七步；

第五步：若第m个像素值的左边的像素值与第m个像素值不相同，第m个像素值的上边的像素值与第m个像素值相同，则将第m个像素点标记为与第m个像素点的上边的像素值所在的区域相同的一个区域，执行第七步；

第六步：若第m个像素点的左边、上边的像素值与第m个像素值都相同，则将第m个像素点标记为与第m个像素点的左边、上边的像素值中像素值较小的像素值所在区域的相同的一个区域，执行第七步；

第七步：判断分割后二值图像中的全部像素点是否标记完，若是，则得到分割后的二值图像的连通域；否则，将分割后的二值图像的像素点编号m＝m+1，执行第二步。

(3o)按照下式，计算分割后二值图像中每个散射点对应区域的区域重心横坐标值和纵坐标值；

$x_c=\frac{\underset{a}{\Σ}\underset{b}{\Σ}af(a,b)}{\underset{a}{\Σ}\underset{b}{\Σ}f(a,b)},y_c=\frac{\underset{a}{\Σ}\underset{b}{\Σ}bf(a,b)}{\underset{a}{\Σ}\underset{b}{\Σ}f(a,b)}$

其中，x_c表示分割后二值图像中每个散射点对应区域重心的横坐标值，y_c表示分割后二值图像中每个散射点对应区域重心的纵坐标值，a表示分割后二值图像中每个散射点对应区域中像素点的横坐标值，b表示分割后二值图像中每个散射点对应区域中像素点的纵坐标值，Σ表示求和操作，f(·)表示分割后二值图像中每个散射点对应的区域中像素点的像素值，f(a,b)表示分割后二值图像中每个散射点对应区域中横纵坐标分别为a、b的像素点的像素值。

(3p)将分割后二值图像中每个散射点对应区域的区域重心横坐标值和纵坐标值，作为分割后二值图像中每个散射点二维坐标的横坐标值和纵坐标值；

(3q)将分割后二值图像中每个散射点二维坐标的横坐标值和纵坐标值，作为微动目标距离-瞬时多普勒图像中每个散射点二维坐标的横坐标值和纵坐标值，并记录该横坐标值和纵坐标值；

(3r)判断是否获得所有时刻微动目标距离-瞬时多普勒图像中散射点的二维坐标，若是，执行步骤(4)，否则，令β＝β+1，执行步骤(3e)，其中，β表示获得微动目标散射点二维坐标时的慢时间单元序号；

步骤4，二维航迹关联。

(4a)生成列数为慢时间个数，行数为微动目标散射点个数的空矩阵，将该空矩阵作为航迹矩阵，将航迹矩阵中每个元素定义为二维向量；

(4b)将完成二维航迹关联时的慢时间单元序号l初始化为1，将航迹矩阵列编号i初始化为1，将散射点序号α初始化为1；

(4c)将航迹矩阵中的第一列元素值分别赋为微动目标各散射点第l个慢时间单元距离-瞬时多普勒图像中的二维坐标值；

(4d)令l＝l+1，其中l表示完成二维航迹关联时的慢时间单元序号，计算航迹矩阵第i列第α个散射点的二维坐标值与第l个慢时间单元距离-瞬时多普勒图像中所有散射点二维坐标值之间的欧氏距离；

(4f)将欧氏距离最小值所对应的第l个慢时间单元距离瞬时多普勒图像中散射点的二维坐标写入航迹矩阵第i+1列，第α行，完成第l个慢时间单元微动目标第α个散射点的航迹关联；

(4g)判断是否完成第l个慢时间单元微动目标所有散射点的航迹关联，若是，执行步骤(4h)，否则，令散射点序号α＝α+1，执行步骤(4f)；

(4h)判断是否完成所有时刻微动目标所有散射点的航迹关联，若是，得到航迹矩阵，执行步骤(5)，否则，航迹矩阵的列编号i＝i+1，执行步骤(4d)；

(4i)保留航迹矩阵每个元素的第一维坐标，删除每个元素的第二维坐标，得到更新后的航迹矩阵。

步骤5，修正航迹矩阵。

(5a)将修正航迹矩阵时的慢时间单元序号p初始化为1。

(5b)采用求根多重信号分类Root-MUSIC方法，对第p个慢时间单元解线频调处理后的信号做角频率估计，得到第p个慢时间单元解线频调处理后信号的n个角频率，其中，n表示微动目标散射点个数。

(5c)按照下式，计算解线频调处理后信号的瞬时斜距：

$R_p=\frac{\mathrm{\ω}\·c\·F_s}{-4\mathrm{\π}\mathrm{\γ}}$

其中，R_p表示第p个慢时间单元解线频调处理后信号的瞬时斜距，ω表示第p个慢时间单元解线频调处理后信号的角频率，c表示光速，γ表示调频率，F_s表示逆合成孔径雷达发射信号的距离采样频率，其取值范围为F_s≥B，B表示逆合成孔径雷达发射信号的带宽；

(5d)分别计算第p个慢时间单元每个瞬时斜距与航迹矩阵第p列中每个散射点对应瞬时斜距的欧氏距离，用欧氏距离最小值对应的瞬时斜距值代替航迹矩阵第p列中对应散射点的瞬时斜距，得到第p列修正后的航迹矩阵；

(5e)判断是否完成所有慢时间单元航迹矩阵的修正，若是，执行步骤(6)，否则，令p＝p+1，执行步骤(5b)。

步骤6，获得微动目标散射点的航迹关联结果。

下面结合附图2对本发明的效果做进一步说明。

1、仿真实验条件：

附图2所示的仿真图在MATLAB7.0软件下进行的，仿真数据的参数如下：雷达带宽为2GHz，载频为10GHz，PRF为2000Hz，观测时间为1s，方位向回波为2000次，距离向采样点数为128，雷达视线的方位角为50.6°，俯仰角为45.3°，仿真数据中目标微动形式为进动，目标包含9个强散射点，自旋角频率为2πrad/s，锥旋角频率为0.8πrad/s，进动角为6.2°。

2、仿真实验内容：

图2(a)是本发明仿真实验过程中所假设的微动目标散射点的分布图，其中，横轴表示散射点的横坐标，单位为米，纵轴表示散射点的纵坐标，单位为米，竖轴表示散射点的竖坐标，单位为米，圆点表示微动目标的散射点所在的位置。图2(b)是本发明解线频调处理后的信号在距离维做一维傅里叶变换，得到目标的距离-慢时间域逆合成孔径雷达回波图。其中，水平坐标表示慢时间，单位为秒，垂直坐标表示瞬时斜距，单位为米。图2(c)是本发明采用距离-瞬时多普勒像获取微动目标散射点二维坐标，进行二维航迹关联后的结果图。其中，水平坐标表示慢时间，单位为秒，垂直坐标表示瞬时斜距，单位均为米，分别用圆圈、加号、五角星、黑点、乘号、星号、方形、菱形、向上的三角形以及向下的三角形对九条曲线进行了标识，以代表不同散射点的航迹。图2(d)是本发明采用求根多重信号分类Root-MUSIC方法，修正航迹矩阵后航迹关联的结果图。其中，水平坐标表示慢时间，单位为秒，垂直坐标表示瞬时斜距，单位为米，分别用圆圈、加号、五角星、黑点、乘号、星号、方形、菱形、向上的三角形以及向下的三角形对九条曲线进行了标识，以代表不同的航迹。

3、仿真效果分析：

从图2(a)仿真实验过程中假设的微动目标散射点分布图中可以看出，微动目标共九个散射点，九个散射点坐标依次是(0,0,1)m，(0.3,03,-0.8)m，(0.3,03,-0.8)m，(-0.3,-03,-0.8)m，(0,-0.3,-0.8)m，(0.5,0,-2)m，(0,0.5,-2)m，(-0.5,0,-2)m，(0,-0.5,-2)m。从图2(b)仿真实验过程中解线频调处理后的信号在距离维做一维傅里叶变换，得到的目标的距离-慢时间域逆合成孔径雷达回波图中可以看出瞬时斜距随慢时间的变化情况，但是，由该图像不能得到每一时刻微动目标散射点的具体二维坐标。由图2(c)本发明仿真实验过程中采于距离-瞬时多普勒像获取微动目标散射点二维坐标，进行二维航迹关联后的结果图中可以看出共有九条散射点航迹，并能够获得每一时刻微动目标散射点的具体二维坐标，但是仍存在不足之处是，在微动目标散射点航迹交叉点处存在关联误差。因此，本发明又采用了求根多重信号分类Root-MUSIC方法，修正航迹矩阵，得到修正后航迹关联的结果图，如图2(d)所示。从图2(d)的结果图中可以看出，采用了求根多重信号分类Root-MUSIC方法修正航迹矩阵，能够弥补图2(c)中微动目标散射点航迹交叉点处存在关联误差的不足，得到了更好的关联效果。

Claims (4)

1.基于距离-瞬时多普勒像微动目标散射点航迹关联方法，包括如下步骤：

(1)接收微动目标逆合成孔径雷达回波；

(2)逆合成孔径雷达回波预处理：

(2a)将逆合成孔径雷达到场景中心距离作为参考距离，将与逆合成孔径雷达发射信号载频、调频率相同、距离为参考距离的线性调频信号作为参考信号；

(2b)将参考信号取共轭后与接收的回波相乘，得到解线频调处理后的信号；

(3)获得微动目标散射点的二维坐标：

(3a)将解线频调处理后的信号在距离维做一维傅里叶变换，得到目标距离-慢时间域的逆合成孔径雷达回波；

(3b)将距离-慢时间域存在目标回波的每个距离单元做短时傅里叶变换，得到每个距离单元的二维时频图像；

(3c)将每个距离单元对应的二维时频图像沿距离维堆叠成距离-瞬时多普勒-慢时间三维矩阵；

(3d)将获得微动目标散射点二维坐标时的慢时间单元序号β初始化为1；

(3e)取三维矩阵在第β个慢时间单元对应的二维矩阵切片，得到第β个慢时间单元的微动目标距离-瞬时多普勒像；

(3f)用第β个慢时间单元的微动目标距离-瞬时多普勒像中每一个像素值除以第β个慢时间单元的微动目标距离-瞬时多普勒像中最大像素值，得到第β个慢时间单元归一化后的微动目标距离-瞬时多普勒像；

(3g)设定门限值th＝0.3；

(3h)将第β个慢时间单元归一化后的微动目标距离-瞬时多普勒像中像素值大于等于门限值的像素点赋值为1，将第β个慢时间单元归一化后的微动目标距离-瞬时多普勒像中像素值小于门限值的像素点赋值为0，得到二值图像；

(3i)将二值图像中像素值为0的像素点值赋值为1，像素值为1的像素点赋值为0，得到二值图像的补图像；

(3j)计算二值图像的补图像中每个像素点到其最近非零像素点的欧氏距离，得到二值图像补图像的距离变换矩阵；

(3k)将二值图像补图像的距离变换矩阵中的每个值取反，得到灰度图像；

(3l)采用分水岭图像分割方法，对灰度图像做图像分割，得到分割后的图像；

(3m)将分割后的图像中像素值为1和0的像素点均赋值为0，将其他像素点赋值为255，得到分割后的二值图像；

(3n)采用四邻域标记方法，得到分割后的二值图像的连通域，将分割后二值图像的每个连通域作为分割后二值图像中每个散射点对应的区域；

(3o)按照下式，计算分割后二值图像中每个散射点对应区域的区域重心横坐标值和纵坐标值；

$x_c=\frac{\underset{a}{\mathrm{\Σ}}\underset{b}{\mathrm{\Σ}}af(a,b)}{\underset{a}{\mathrm{\Σ}}\underset{b}{\mathrm{\Σ}}f(a,b)},y_c=\frac{\underset{a}{\mathrm{\Σ}}\underset{b}{\mathrm{\Σ}}bf(a,b)}{\underset{a}{\mathrm{\Σ}}\underset{b}{\mathrm{\Σ}}f(a,b)}$

其中，x_c表示分割后二值图像中每个散射点对应区域重心的横坐标值，y_c表示分割后二值图像中每个散射点对应区域重心的纵坐标值，a表示分割后二值图像中每个散射点对应区域中像素点的横坐标值，b表示分割后二值图像中每个散射点对应区域中像素点的纵坐标值，Σ表示求和操作，f(·)表示分割后二值图像中每个散射点对应的区域中像素点的像素值，f(a,b)表示分割后二值图像中每个散射点对应区域中横纵坐标分别为a、b的像素点的像素值；

(3p)将分割后二值图像中每个散射点对应区域的区域重心横坐标值和纵坐标值，作为分割后二值图像中每个散射点二维坐标的横坐标值和纵坐标值；

(3q)将分割后二值图像中每个散射点二维坐标的横坐标值和纵坐标值，作为微动目标距离-瞬时多普勒图像中每个散射点二维坐标的横坐标值和纵坐标值，并记录该横坐标值和纵坐标值；

(3r)判断是否获得所有时刻微动目标距离-瞬时多普勒图像中散射点的二维坐标，若是，执行步骤(4)，否则，令β＝β+1，执行步骤(3e)，其中，β表示获得微动目标散射点二维坐标时的慢时间单元序号；

(4)二维航迹关联：

(4a)生成列数为慢时间个数，行数为微动目标散射点个数的空矩阵，将该空矩阵作为航迹矩阵，将航迹矩阵中每个元素定义为二维向量；

(4b)将完成二维航迹关联时的慢时间单元序号l初始化为1，将航迹矩阵列编号i初始化为1，将散射点序号α初始化为1；

(4c)将航迹矩阵中的第一列元素值分别赋为微动目标各散射点第l个慢时间单元距离-瞬时多普勒图像中的二维坐标值；

(4d)令l＝l+1，其中，l表示完成二维航迹关联时的慢时间单元序号，计算航迹矩阵第i列第α个散射点的二维坐标值与第l个慢时间单元距离-瞬时多普勒图像中所有散射点二维坐标值之间的欧氏距离；

(4e)将欧氏距离最小值所对应的第l个慢时间单元距离瞬时多普勒图像中散射点的二维坐标写入航迹矩阵第i+1列，第α行，完成第l个慢时间单元微动目标第α个散射点的航迹关联；

(4f)判断是否完成第l个慢时间单元微动目标所有散射点的航迹关联，若是，执行步骤(4g)，否则，令散射点序号α＝α+1，执行步骤(4e)；

(4g)判断是否完成所有时刻微动目标所有散射点的航迹关联，若是，得到航迹矩阵，执行步骤(5)；否则，航迹矩阵的列编号i＝i+1，执行步骤(4d)；

(4h)保留航迹矩阵每个元素的第一维坐标，删除每个元素的第二维坐标，得到更新后的航迹矩阵；

(5)修正航迹矩阵：

(5a)将修正航迹矩阵时的慢时间单元序号p初始化为1；

(5b)采用求根多重信号分类Root-MUSIC方法，对第p个慢时间单元解线频调处理后的信号做角频率估计，得到第p个慢时间单元解线频调处理后信号的n个角频率，其中，n表示微动目标散射点个数；

(5c)按照下式，计算解线调频处理后信号的瞬时斜距：

$R_p=\frac{\mathrm{\ω}\·c\·F_s}{-4\mathrm{\π}\mathrm{\γ}}$

其中，R_p表示第p个慢时间单元解线调频处理后信号的瞬时斜距，ω表示第p个慢时间单元解线调频处理后信号的角频率，c表示光速，γ表示调频率，F_s表示逆合成孔径雷达发射信号的距离采样频率，其取值范围为F_s≥B，B表示逆合成孔径雷达发射信号的带宽；

(5d)分别计算第p个慢时间单元每个瞬时斜距与航迹矩阵第p列中每个散射点对应瞬时斜距的欧氏距离，用欧氏距离最小值对应的瞬时斜距值代替航迹矩阵第p列中对应散射点的瞬时斜距，得到第p列修正后的航迹矩阵；

(5e)判断是否完成所有慢时间单元航迹矩阵的修正，若是，执行步骤(6)，否则，令p＝p+1，执行步骤(5b)；

(6)获得微动目标散射点的航迹关联结果。

2.根据权利要求1所述的基于距离-瞬时多普勒像的微动目标散射点航关联方法，其特征在于，步骤(3b)中所述的短时傅里叶变换的方法的具体步骤如下：

STFT_s(t′,f)＝∫s·w(t-t′)exp{-j2πft}dt

其中，STFT(·)表示短时傅里叶变换操作，STFT_s(·)表示对存在目标回波的距离单元信号做短时傅里叶变换操作，t′表示窗函数中心所在的时间点，f表示瞬时多普勒，∫表示积分操作，s表示解线频调处理后信号中存在目标回波的距离单元信号，w(·)表示窗函数，通常选高斯窗或者汉明窗，窗长为N_a表示慢时间单元的总数，w(t-t′)表示沿时间轴平移后的窗函数，exp(·)表示取以e为底的指数操作，j表示虚数单位，t表示慢时间。

3.根据权利要求1所述的基于距离-瞬时多普勒像的微动目标散射点航迹关联方法，其特征在于，步骤(3l)中所述分水岭图像分割方法的具体步骤如下：

第一步：对灰度图像中的各像素点按灰度值大小从小到大排序，将所有像素点中相同灰度值的像素点作为同一个层级；

第二步：初始化层级编号k，令k＝1；

第三步：处理第k个层级所有的像素点，如果该像素点的邻域已经被标识属于某一个区域，则将该像素点加入一个先进先出队列；

第四步：判断先进先出队列是否为空，若非空，则将先进先出队列的元素编号n初始化为1，执行第五步；否则，执行第七步；

第五步：取出先进先出队列中的第n个元素，扫描第n个元素的邻域像素，如果第n个元素及其邻域像素的灰度值属于同一层级，则用邻域像素的标识刷新先进先出的队列的第n个元素像素的标识；

第六步：令先进先出队列的元素编号n＝n+1，重复第五步，直到先进先出的队列为空；

第七步：扫描的像素点，判断第k个层级像素点是否都被标识，若有像素点未被标识，则将该像素点的灰度值加1，若第k个层级像素点都被标识，执行第八步；

第八步：判断是否处理完所有层级的像素点，若是，得到分割后的图像，否则，令层级编号k＝k+1，执行第二步。

4.根据权利要求1所述的基于距离-瞬时多普勒像的微动目标散射点航迹关联方法，其特征在于，步骤(3n)中所述的四邻域标记方法的具体步骤如下：

第一步：将分割后的二值图像的像素点编号m初始化为K+2，其中，K表示距离-瞬时多普勒像的列数；

第二步：判断m是否为K的整数倍，若是，则m＝m+1，执行第三步，否则直接执行第三步，其中，K表示距离-瞬时多普勒像的列数；

第三步：若第m个像素值的左边、上边的像素值与第m个像素值都不相同，则将该像素值所在的区域标记与为第m个像素值的左边、上边的像素值所在的区域都不相同的一个区域，执行第七步；

第四步：若第m个像素值的左边的像素值与第m个像素值相同，第m个像素点的上边的像素值与第m个像素值不相同，则将第m个像素点标记为与第m个像素点的左边的像素值所在的区域相同的一个区域，执行第七步；

第五步：若第m个像素值的左边的像素值与第m个像素值不相同，第m个像素点的上边的像素值与第m个像素值相同，则将第m个像素点标记为与第m个像素点的上边的像素值所在的区域相同的一个区域，执行第七步；

第六步：若第m个像素点的左边、上边的像素值与第m个像素值都相同，则将第m个像素点标记为与第m个像素点的左边、上边的像素值中像素值较小的像素值所在区域的相同的一个区域，执行第七步；

第七步：判断分割后二值图像中的全部像素点是否标记完，若是，则得到分割后的二值图像的连通域；否则，将分割后的二值图像的像素点编号m＝m+1，执行第二步。