CN105447867B - 基于isar图像的空间目标姿态估计方法 - Google Patents

Abstract

本发明公开了一种基于ISAR图像的空间目标姿态估计方法，包括以下步骤：(1)ISAR雷达系统接收设定时间段内的回波信号，采用距离‑多普勒算法对其进行处理，得到P幅目标区域的雷达图像，并依次进行滤波和对数变换，再进行目标与背景的分割，得到P幅目标‑背景的二值图像：(2)对每幅目标‑背景的二值图像进行形态学处理，得到提取边界后的图像；(3)提取每幅提取边界后的图像的N个线性结构；(4)将第p幅提取边界后的图像与第p‑1幅提取边界后的图像进行关联处理；(5)设定每个线性结构的姿态参数集，并构造每个线性结构的代价函数，在每个线性结构的姿态参数集中搜索使其代价函数达到最小的一组姿态参数，即完成对空间目标的姿态估计。

Description

基于ISAR图像的空间目标姿态估计方法

技术领域

本发明属于雷达技术领域，特别涉及一种基于ISAR图像的空间目标姿态估计方法，可用于对空间目标进行姿态估计。

背景技术

空间目标姿态估计对了解空间目标的动作意图以及判断空间目标的状态具有非常重要的意义。利用逆合成孔径雷达(Inverse Synthetic Aperture Radar，ISAR)对空间目标连续测量获得的序列图像来精确获得空间目标所配备太阳能翼、平板天线等重要载荷部件的绝对姿态的技术可实际应用于包括空间目标故障分析、威胁度评判等民用和军用领域，是目前能够实现空间目标姿态估计的实用技术。

目前，测定空间目标姿态常用主动和被动两种方式：主动测量是在飞行器中装载GPS接收机和惯性测量器对飞行器的位置和姿态进行测量，这种方法称为运动目标内测，运动目标内测要求飞行器配备专门的测量系统，系统复杂度高。被动测量是利用光电或雷达传感器对空间目标进行连续高分辨成像，根据成像序列信息来进行空间目标的姿态确定，被动测量方式具有与空间目标状态无关的优点，也可以直接应用到未知目标或者失控目标的姿态测量应用，但利用光电传感器进行空间目标姿态测定易受到天气、时间等因素的制约，在实际空间目标姿态测量中存在局限性。

发明内容

针对上述现有技术的不足，本发明的目的在于提出一种基于ISAR图像的空间目标姿态估计方法，该方法将获得的ISAR雷达图像依次进行滤波、对数变换、目标与背景的分割和形态学处理，得到提取边界后的图像，再通过Radon变换提取空间目标在三维空间中较为线束的线性结构，并将每一幅提取边界后的图像与前一幅提取边界后的图像进行关联处理后，构造代价函数，针对每个线性结构，搜索使代价函数为最小的姿态参数，完成对空间目标的姿态估计。

为了实现上述技术目的，本发明采用如下技术方案予以实现。

一种基于ISAR图像的空间目标姿态估计方法，其特征在于，包括以下步骤：

步骤1，由ISAR雷达系统接收设定时间段内的回波信号；对设定时间段内的回波信号进行距离-多普勒处理，得到目标区域的雷达图像序列，该序列包括P幅目标区域的雷达图像；对第p幅目标区域的雷达图像依次进行滤波和对数变换，得到第p幅预处理后的目标区域的雷达图像，p＝1，2，...，P；

步骤2，对第p幅预处理后的目标区域的雷达图像进行目标与背景的分割，得到第p幅目标-背景的二值图像；

步骤3，对第p幅目标-背景的二值图像进行形态学处理，提取目标的边界，得到第p幅提取边界后的图像；

步骤4，提取第p幅提取边界后的图像的N个线性结构；记录第p幅提取边界后的图像中每个线性结构的倾角、中心位置和长度，其中，N是正整数；

步骤5，将第一幅提取边界后的图像作为第一幅关联后的图像；再从第二幅提取边界后的图像开始，将第p幅提取边界后的图像与第p-1幅提取边界后的图像进行关联处理，即将第p幅提取边界后的图像的N个线性结构与第p-1幅提取边界后的图像的N个线性结构进行匹配标记，得到第二幅到第P幅关联后的图像；最终，得到P幅关联后的图像；

步骤6，设定第p幅关联后的图像的第i个线性结构的姿态参数包括仰角θ_p，i和方位角φ_p，i，则第i个线性结构的姿态参数集为C_i，C_i＝{(θ_1，i，φ_1，i)，...，(θ_p，i，φ_p，i)，...，(θ_P，i，φ_P，i)}，其中，-90°≤θ≤90°，0°≤φ≤180°；构造第i个线性结构的代价函数，并在第i个线性结构的姿态参数集C_i中搜索使其代价函数达到最小的一组姿态参数，将该姿态参数作为第i个线性结构的最终的姿态参数，i＝1，2，...，N，即完成对空间目标的姿态估计。

本发明与现有技术相比所具有的优点：

1)本发明利用空间目标的ISAR图像序列的相对投影变化关系，结合空间投影理论，通过二维ISAR图像序列对空间目标的三维绝对姿态实现精确反演，本发明具有很高的实用性。

2)本发明利用形态学、Radon变换等预处理技术有效提高了ISAR图像的线性结构提取的精确性，可针对具有复杂几何形态的空间目标进行姿态估计。

附图说明

图1是本发明的流程图；

图2是本发明仿真采用的空间目标的结构图，其中，空间三维坐标系各坐标轴的单位为米(m)；

图3a、图3b和图3c分别是提取的空间目标的第1幅、第6幅和第14幅提取边界后的图像的线性结构的示意图，其中，横坐标为像素点，纵坐标为像素点；

图4a、图4b和图4c分别是空间目标的第1幅、第6幅和第14幅提取边界后的图像的关联后的线性结构的示意图，其中，横坐标为像素点，纵坐标为像素点；

具体实施方式：

参照图1，本发明的基于ISAR图像的空间目标姿态估计方法，具体实现步骤如下：

步骤1，由ISAR雷达系统接收设定时间段内的回波信号；采用距离-多普勒算法对设定时间段内的回波信号进行处理，得到目标区域的雷达图像序列，该序列包括P幅目标区域的雷达图像；对第p幅目标区域的雷达图像依次进行滤波和对数变换，得到第p幅预处理后的目标区域的雷达图像，p＝1，2，...，P。

步骤1中，第p幅预处理后的目标区域的雷达图像中像素点(m，n)的像素值为F_p(m，n)：

F_p(m，n)＝c₁×ln[f_p(m，n)+c₂]+c₃

其中，f_p(m，n)为第p幅目标区域的雷达图像中像素点(m，n)的像素值，c₁、c₂和c₃对数变换常量。

步骤2，采用k-mean算法对第p幅预处理后的目标区域的雷达图像进行目标与背景的分割，得到第p幅目标-背景的二值图像。

步骤2的具体子步骤为：

2.1选取第p幅预处理后的目标区域的雷达图像中的最大像素值作为目标类ω₁的聚类中心μ_p，1，选取预处理后的目标区域的雷达图像中的最小像素值作为非目标类ω₂的聚类中心μ_p，2；

2.2计算分类阈值T，对第p幅预处理后的目标区域的雷达图像中的所有像素点进行分类，若像素点(m，n)的像素值F_p(m，n)＞T，则像素点(m，n)∈ω₁，反之，则像素点(m，n)∈ω₂；

2.3分别计算聚类后目标类ω₁的聚类中心μ′_p，1和聚类后非目标类ω₂的聚类中心μ′_p，2：

其中，N₁为属于目标类ω₁的像素点个数，N₂为属于非目标类ω₂的像素点个数；

2.4如果聚类后目标类ω₁的聚类中心μ′_p，1不等于目标类ω₁的聚类中心μ_p，1，或者聚类后非目标类ω₂的聚类中心μ′_p，2不等于非目标类ω₂的聚类中心μ_p，2，则令μ_p，1＝μ′_p，1，μ_p，2＝μ′_p，2，返回步骤2.2，

如果聚类后目标类ω₁的聚类中心μ′_p，1等于目标类ω₁的聚类中心μ_p，1，且聚类后非目标类ω₂的聚类中心μ′_p，2等于非目标类ω₂的聚类中心μ_p，2；并将N₁个属于目标类ω₁的像素点的像素值赋值为1，将N₂个属于非目标类ω₂的像素点的像素值赋值为0，得到第p幅目标-背景的二值图像。

步骤3，对第p幅目标-背景的二值图像进行形态学处理，提取目标的边界，得到第p幅提取边界后的图像。

步骤3的具体子步骤为：

3.1对第p幅目标-背景的二值图像进行膨胀处理，填补目标中的空洞，同时扩张目标-背景的二值图像的边界，得到第p幅膨胀图像BW1_p；

3.2对第p幅膨胀图像BW1_p进行腐蚀处理，同时收缩第p幅膨胀图像BW1_p的边界，得到第p幅腐蚀图像BW2_p；

3.3根据第p幅膨胀图像BW1_p和第p幅腐蚀图像BW2_p，得到第p幅提取边界后的图像BW_p，其中，第p幅提取边界后的图像BW_p中像素点(m，n)的像素值为

其中，为第p幅膨胀图像BW1_p中像素点(m，n)的像素值，为第p幅膨胀图像BW2_p中像素点(m，n)的像素值。

步骤4，对第p幅提取边界后的图像进行Radon变换，提取第p幅提取边界后的图像的N个线性结构；记录第p幅提取边界后的图像中每个线性结构的倾角、中心位置和长度，其中，N是正整数。

步骤5，将第一幅提取边界后的图像作为第一幅关联后的图像；再从第二幅提取边界后的图像开始，将第p幅提取边界后的图像与第p-1幅提取边界后的图像进行关联处理，即将第p幅提取边界后的图像的N个线性结构与第p-1幅提取边界后的图像的N个线性结构进行匹配标记，得到第二幅到第P幅关联后的图像；最终，得到P幅关联后的图像。

步骤5中，将第p幅提取边界后的图像与第p-1幅提取边界后的图像进行关联处理，其具体子步骤为：

5.1计算第p幅提取边界后的图像的第i个线性结构相对于第p-1幅提取边界后的图像的第j个线性结构的变化量：

其中，Δψ_p(j，i)表示第p幅提取边界后的图像的第i个线性结构与第p-1幅提取边界后的图像的第j个线性结构的倾角变化量，表示第p幅提取边界后的图像的第i个线性结构的倾角，表示第p-1幅提取边界后的图像的第j个线性结构的倾角，

Δpos_p(j，i)表示第p幅提取边界后的图像的第i个线性结构与第p-1幅提取边界后的图像的第j个线性结构的中心位置变化量，表示第p幅提取边界后的图像的第i个线性结构的中心位置，表示第p-1幅提取边界后的图像的第j个线性结构的中心位置，

Δl_p(j，i)表示第p幅提取边界后的图像的第i个线性结构与第p-1幅提取边界后的图像的第j个线性结构的长度变化量，表示第p幅提取边界后的图像的第i个线性结构的长度，表示第p-1幅提取边界后的图像的第j个线性结构的长度，

i＝1，2，...，N，j＝1，2，...，N；

5.2计算第p幅提取边界后的图像的第i个线性结构与第p-1幅提取边界后的图像的第j个线性结构的相似值J_p(j，i)：

J_p(j，i)＝α₁×Δψ(j，i)+α₂×Δpos(j，i)+α₃×Δl(j，i)

其中，α₁、α₂、α₃为标准化因子；

定义第p幅提取边界后的图像的相似矩阵J_p为：

5.3设定l为搜索次数，并设置搜索次数l的初始值为1；

5.4对第p幅提取边界后的图像的相似矩阵J_p进行第l次搜索，搜索到其最小元素为J_p(j_l，i_l)，则将第p幅提取边界后的图像的第i_l个线性结构与第p-1幅提取边界后的图像的第j_l个线性结构进行匹配标记；

将第p幅提取边界后的图像的相似矩阵J_p的第j_l行的元素和第i_l列的元素去掉，保留其他元素，则将第p幅提取边界后的图像的相似矩阵J_p更新为(N-l)×(N-l)维的矩阵；

5.5如果搜索次数l＜N-1，则令搜索次数l曾加1，返回步骤5.4，

如果搜索次数l＝N-1，则停止第p幅提取边界后的图像的相似矩阵J_p进行搜索，第p幅提取边界后的图像的相似矩阵J_p中唯一的元素为J_p(j₀，i₀)，将第p幅提取边界后的图像的第i₀个线性结构与第p-1幅提取边界后的图像的第j₀个线性结构进行匹配标记，则将第p幅提取边界后的图像的N个线性结构与第p-1幅提取边界后的图像的N个线性结构均完成匹配标记，得到第二到第P幅关联后的图像。

步骤6，设定第p幅关联后的图像的第i个线性结构的姿态参数包括仰角θ_p，i和方位角φ_p，i，则第i个线性结构的姿态参数集为C_i，C_i＝{(θ_1，i，φ_1，i)，...，(θ_p，i，φ_p，i)，...，(θ_P，i，φ_P，i)}，其中，-90°≤θ≤90°，0°≤φ≤180°；构造第i个线性结构的代价函数，并在第i个线性结构的姿态参数集C_i中搜索使其代价函数达到最小的一组姿态参数，将该姿态参数作为第i个线性结构的最终的姿态参数，i＝1，2，...，N，即完成对空间目标的姿态估计。

步骤6的具体子步骤为：

6.1提取所有P幅关联后的图像的第i个线性结构的倾角作为观测值，i＝1，2，...，N，N为线性结构的个数；

6.2构造第i个线性结构的代价函数J_i(θ，φ)：

其中，A_p为第p幅提取边界后的图像的投影矩阵，其表达式为：

其中，是第p幅提取边界后的图像对应的雷达视线方向与空间三维坐标系中水平面xoy的夹角，为第p幅提取边界后的图像对应的雷达视线方向在空间三维坐标系中水平面xoy的分量与y轴的夹角，ΔR＝c/(2×f_s)，c为光速，f_s为雷达采样频率，ΔR′＝-λ/(2×prf×N_d)，λ为雷达信号波长，prf为脉冲重复频率，N_d为距离多普勒像的方位脉冲数，“·”表示求导；

6.3在第i个线性结构的姿态参数集C_i中，搜索使第i个线性结构的代价函数J_i(θ，φ)达到最小的一组姿态参数将该姿态参数作为第i个线性结构的最终的姿态参数，即完成对空间目标的姿态估计。

本发明的效果可通多以下仿真实验作进一步说明：

1)仿真条件：

本发明仿真采用的空间目标的结构如图3所示，该目标的太阳能翼边界有明显的线性结构；本发明仿真的主要参数，如表1所示：

表1 ISAR雷达系统的主要参数

ISAR图像尺寸	512×512
		波长λ	0.018m
发射信号的带宽B	1GHz
		发射信号中心频率fc	16.7GHz
脉冲重复频率prf	100Hz

2.仿真内容及结果分析：

仿真1：用本发明方法对图3的空间目标的15幅提取边界后的图像进行线性结构的提取，其中，空间目标的第1幅、第6幅和第14幅提取边界后的图像的线性结构分别如图3a、图3b和图3c所示。

仿真2：用本发明方法对图3的空间目标的太阳能翼边界的三个线性结构进行关联，其中，空间目标的第1幅、第6幅和第14幅提取边界后的图像的关联后的线性结构分别如图4a、图4b和图4c所示，图中，1、2和3分别表示第1个、第2个和第3个线性结构。

仿真3，用本发明方法对关联后的三个线性结构，利用代价函数分别搜索其达到最小值的三个线性结构的姿态参数，结果如表2所示。

表2三个线性结构的姿态估计结果

线性结构	θ(度)	φ(度)	估计指向	估计指向与真实指向的夹角
					1	-9.6	171.3	(-0.14，0.98，0.17)	0.51度
2	-10	172.1	(-0.14，0.98，0.17)	0.37度
					3	75.5	121.5	(0.20，-0.14，0.97)	1.05度

3.仿真结果分析：

从图3a、图3b和图3c中可以看出，空间目标的太阳能翼边界的线性结构的提取基本正确，但是空间目标的非太阳能翼部分在一定程度上影响了太阳能翼边界的线性结构提取的精确度。

从图4a、图4b和图4c中可以看出，利用三个线性结构的变化信息可以很好地完成各幅提取边界后的图像之间线性结构的关联。

从表2中可以看出，估计的线性结构的指向与真实的线性结构的指向基本一致，平均误差在1度之内(为0.63度)，而且可以发现线性结构1与线性结构2基本平行(夹角为0.88度)，线性结构1与线性结构3基本垂直(夹角为90.12度)，线性结构2与线性结构3基本垂直(夹角为90.67度)，这与空间目标的结构一致，则空间目标的姿态可以确定。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围；这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (5)

1.一种基于ISAR图像的空间目标姿态估计方法，其特征在于，包括以下步骤：

步骤1，由ISAR雷达系统接收设定时间段内的回波信号；对设定时间段内的回波信号进行距离-多普勒处理，得到目标区域的雷达图像序列，该序列包括P幅目标区域的雷达图像；对第p幅目标区域的雷达图像依次进行滤波和对数变换，得到第p幅预处理后的目标区域的雷达图像，p＝1，2，...，P；

步骤2，对第p幅预处理后的目标区域的雷达图像进行目标与背景的分割，得到第p幅目标-背景的二值图像；

步骤3，对第p幅目标-背景的二值图像提取目标的边界，得到第p幅提取边界后的图像；

步骤4，提取第p幅提取边界后的图像的N个线性结构；记录第p幅提取边界后的图像中每个线性结构的倾角、中心位置和长度，其中，N是正整数；

步骤5，将第一幅提取边界后的图像作为第一幅关联后的图像；再从第二幅提取边界后的图像开始，将第p幅提取边界后的图像与第p-1幅提取边界后的图像进行关联处理，即将第p幅提取边界后的图像的N个线性结构与第p-1幅提取边界后的图像的N个线性结构进行匹配标记，得到第二幅到第P幅关联后的图像；最终，得到P幅关联后的图像；

步骤6，设定第p幅关联后的图像的第i个线性结构的姿态参数包括仰角θ_p，i和方位角φ_p，i，则第i个线性结构的姿态参数集为C_i，C_i＝{(θ_1，i，φ_1，i)，...，(θ_p，i，φ_p，i)，...，(θ_P，i，φ_p，i)}，其中，-90°≤θ≤90°，0°≤φ≤180°；构造第i个线性结构的代价函数，并在第i个线性结构的姿态参数集C_i中搜索使其代价函数达到最小的一组姿态参数，将该姿态参数作为第i个线性结构的最终的姿态参数，i＝1，2，...，N，即完成对空间目标的姿态估计；

其中，步骤6的具体子步骤为：

6.1提取所有P幅关联后的图像的第i个线性结构的倾角作为观测值，i＝1，2，...，N，N为线性结构的个数；

6.2构造第i个线性结构的代价函数J_i(θ，φ)：

其中，A_p为第p幅提取边界后的图像的投影矩阵，其表达式为：

其中，是第p幅提取边界后的图像对应的雷达视线方向与空间三维坐标系中水平面xoy的夹角，为第p幅提取边界后的图像对应的雷达视线方向在空间三维坐标系中水平面xoy的分量与y轴的夹角，ΔR＝c/(2×f_s)，c为光速，f_s为雷达采样频率，ΔR′＝-λ/(2×prf×N_d)，λ为雷达信号波长，prf为脉冲重复频率，N_d为距离-多普勒向的方位脉冲数，“·”表示求导；

6.3在第i个线性结构的姿态参数集C_i中，搜索使第i个线性结构的代价函数J_i(θ，φ)达到最小的一组姿态参数将该姿态参数作为第i个线性结构的最终的姿态参数，即完成对空间目标的姿态估计。

2.如权利要求1所述的基于ISAR图像的空间目标姿态估计方法，其特征在于，步骤1中，第p幅预处理后的目标区域的雷达图像中像素点(m，n)的像素值为F_p(m，n)：

F_p(m，n)＝c₁×ln[f_p(m，n)+c₂]+c₃

其中，f_p(m，n)为第p幅目标区域的雷达图像中像素点(m，n)的像素值，c₁、c₂和c₃对数变换常量。

3.如权利要求1所述的基于ISAR图像的空间目标姿态估计方法，其特征在于，步骤2中，采用k-mean算法对第p幅预处理后的目标区域的雷达图像进行目标与背景的分割，步骤2的具体子步骤为：

2.1选取第p幅预处理后的目标区域的雷达图像中的最大像素值作为目标类ω₁的聚类中心μ_p，1，选取预处理后的目标区域的雷达图像中的最小像素值作为非目标类ω₂的聚类中心μ_p，2；

2.2计算分类阈值T，对第p幅预处理后的目标区域的雷达图像中的所有像素点进行分类，若像素点(m，n)的像素值F_p(m，n)＞T，则像素点(m，n)∈ω₁，反之，则像素点(m，n)∈ω₂；

2.3分别计算聚类后目标类ω₁的聚类中心μ′_p，1和聚类后非目标类ω₂的聚类中心μ′_p，2：

其中，N₁为属于目标类ω₁的像素点个数，N₂为属于非目标类ω₂的像素点个数；

2.4如果聚类后目标类ω₁的聚类中心μ′_p，1不等于目标类ω₁的聚类中心μ_p，1，或者聚类后非目标类ω₂的聚类中心μ′_p，2不等于非目标类ω₂的聚类中心μ_p，2，则令μ_p，1＝μ′_p，1，μ_p，2＝μ′_p，2，返回步骤2.2，

如果聚类后目标类ω₁的聚类中心μ′_p，1等于目标类ω₁的聚类中心μ_p，1，且聚类后非目标类ω₂的聚类中心μ′_p，2等于非目标类ω₂的聚类中心μ_p，2，将N₁个属于目标类ω₁的像素点的像素值赋值为1，将N₂个属于非目标类ω₂的像素点的像素值赋值为0，得到第p幅目标-背景的二值图像。

4.如权利要求1所述的基于ISAR图像的空间目标姿态估计方法，其特征在于，步骤3的具体子步骤为：

3.1对第p幅目标-背景的二值图像进行膨胀处理，填补目标中的空洞，同时扩张目标-背景的二值图像的边界，得到第p幅膨胀图像BW1_p；

3.2对第p幅膨胀图像BW1_p进行腐蚀处理，同时收缩第p幅膨胀图像BW1_p的边界，得到第p幅腐蚀图像BW2_p；

3.3根据第p幅膨胀图像BW1_p和第p幅腐蚀图像BW2_p，得到第p幅提取边界后的图像BW_p，其中，第p幅提取边界后的图像BW_p中像素点(m，n)的像素值为

其中，为第p幅膨胀图像BW1_p中像素点(m，n)的像素值，为第p幅膨胀图像BW2_p中像素点(m，n)的像素值。

5.如权利要求1所述的基于ISAR图像的空间目标姿态估计方法，其特征在于，步骤5中，将第p幅提取边界后的图像与第p-1幅提取边界后的图像进行关联处理，其具体子步骤为：

5.1计算第p幅提取边界后的图像的第i个线性结构相对于第p-1幅提取边界后的图像的第j个线性结构的变化量：

其中，Δψ_p(j，i)表示第p幅提取边界后的图像的第i个线性结构与第p-1幅提取边界后的图像的第j个线性结构的倾角变化量，表示第p幅提取边界后的图像的第i个线性结构的倾角，表示第p-1幅提取边界后的图像的第j个线性结构的倾角，

Δpos_p(j，i)表示第p幅提取边界后的图像的第i个线性结构与第p-1幅提取边界后的图像的第j个线性结构的中心位置变化量，表示第p幅提取边界后的图像的第i个线性结构的中心位置，表示第p-1幅提取边界后的图像的第j个线性结构的中心位置，

Δl_p(j，i)表示第p幅提取边界后的图像的第i个线性结构与第p-1幅提取边界后的图像的第j个线性结构的长度变化量，表示第p幅提取边界后的图像的第i个线性结构的长度，表示第p-1幅提取边界后的图像的第j个线性结构的长度，

i＝1，2，...，N，j＝1，2，...，N；

5.2计算第p幅提取边界后的图像的第i个线性结构与第p-1幅提取边界后的图像的第j个线性结构的相似值J_p(j，i)：

J_p(j，i)＝α₁×Δψ_p(j，i)+α₂×Δpos_p(j，i)+α₃×Δl_p(j，i)

其中，α₁、α₂、α₃为标准化因子；

定义第p幅提取边界后的图像的相似矩阵J_p为：

5.3设定l为搜索次数，并设置搜索次数l的初始值为1；

5.4对第p幅提取边界后的图像的相似矩阵J_p进行第l次搜索，搜索到其最小元素为J_p(j_l，i_l)，则将第p幅提取边界后的图像的第i_l个线性结构与第p-1幅提取边界后的图像的第j_l个线性结构进行匹配标记；

将第p幅提取边界后的图像的相似矩阵J_p的第j_l行的元素和第i_l列的元素去掉，保留其他元素，则将第p幅提取边界后的图像的相似矩阵J_p更新为(N-l)×(N-l)维的矩阵；

5.5如果搜索次数l＜N-1，则令搜索次数l增加1，返回步骤5.4，

如果搜索次数l＝N-1，则停止对第p幅提取边界后的图像的相似矩阵J_p的搜索，第p幅提取边界后的图像的相似矩阵J_p中唯一的元素为J_p(j₀，i₀)，将第p幅提取边界后的图像的第i₀个线性结构与第p-1幅提取边界后的图像的第j₀个线性结构进行匹配标记，则将第p幅提取边界后的图像的N个线性结构与第p-1幅提取边界后的图像的N个线性结构均完成匹配标记，得到第二到第P幅关联后的图像。