CN108805904B - 一种基于卫星序列图像的运动舰船检测与跟踪方法 - Google Patents

Abstract

本发明涉及一种基于卫星序列图像的运动舰船检测与跟踪方法，属于视频卫星图像处理技术。本发明的步骤如下：不断获取包括运动舰船的卫星序列图像，输入相邻两帧图像；对两帧图像进行图像配准；计算两帧图像的差值图像并计算差值图像的多尺度显著图；基于显著图提取二值化前景图像之后，根据背景与舰船区域的灰度差异分别在两帧图像中提取运动舰船区域；将两帧图像中的运动舰船区域进行匹配，得到运动舰船匹配对；以每三帧图像的中间帧为关联帧，对运动舰船进行关联，从而实现舰船跟踪。

Description

一种基于卫星序列图像的运动舰船检测与跟踪方法

技术领域

本发明涉及一种基于卫星序列图像的运动舰船检测与跟踪方法，属于视频卫星图像处理技术。

背景技术

我国具有辽阔的海岸线，海洋面积广阔，舰船是海上重要的运输工具和军事目标，对相关海洋区域舰船目标进行监视与跟踪具有重要的现实意义。在民用领域，可以监控一些重要港口或海域，在海洋安全、海事管理、海洋交通管制、非法走私、非法捕鱼等应用前景广阔。在军事领域，舰船目标检测与跟踪技术在海洋侦查、监视领域具有重要的地位，可以监视重要港口或海域的舰船动态，分析舰船的位置、航向及航速，识别舰船类别等。

目前，常用的地基舰船跟踪手段，主要通过安装在海面浮漂上的数字摄像机来检测与跟踪海上舰船目标，但是这种方式覆盖范围小、隐蔽性极差等缺点使其在海洋舰船检测与跟踪方面的应用受限。而遥感技术实现了远距离和非接触的目标探测能力，具有侦查范围广、不受地面和空中条件限制、隐蔽性强、信息获取及时等优点。

随着遥感技术的发展，利用遥感图像进行舰船目标检测或跟踪成为可能。研究人员主要利用SAR和光学静态遥感图像开展舰船检测方法研究。对于SAR图像舰船目标检测与识别，主要是利用舰船目标与水体之间不同的电磁散射特性在SAR图像上所表现出的差异。相比于SAR图像数据，光学遥感图像舰船目标检测研究起步晚一些，但是，随着光学卫星分辨率的不断提高，利用光学静态图像舰船检测的研究不断增多。海陆背景下的舰船目标检测主要包括海陆分离、舰船目标检测、舰船目标分类等。但是，无论是SAR遥感图像还是光学静态遥感图像，只能实现基于单幅图像的舰船检测，对于舰船的运动信息无能为力。

近年来，随着高分辨率视频卫星系统的出现，特别是2013年1m分辨率的Skybox卫星的成功应用，我国也特别重视视频卫星的发展，分别于2015年和2017年发射了吉林一号和欧比特视频卫星，因此，利用卫星跟踪与监视大范围海域的运动舰船也将成为可能。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供一种基于卫星序列图像的运动舰船检测与跟踪方法，其优点在于实现卫星监视海域的多运动舰船目标的检测与跟踪，可有效获取多个运动舰船轨迹、航速、航向等动态信息。

本发明的技术解决方案是：

一种基于卫星序列图像的运动舰船检测与跟踪方法，该方法的步骤包括：

S1、通过卫星不断获取包括运动舰船的卫星序列图像，获取的卫星序列图像中，第t-1帧图像表示为I_t-1，第t帧图像表示为I_t，其中，t为大于等于2的正整数；

S2、以步骤S1中的第t-1帧图像I_t-1为参考图像，以第t帧图像I_t为待配准图像，采用基于金字塔的最小二乘匹配方法计算I_t-1和I_t间的仿射变换矩阵，依据得到的仿射变换矩阵对I_t进行变换得到配准后的第t帧图像

所述的采用基于金字塔的最小二乘匹配方法计算I_t-1和I_t间的仿射变换矩阵，依据得到的仿射变换矩阵对I_t进行变换得到配准后的第t帧图像

的方法如下：

S21，对第t-1帧图像和第t帧图像分别构建K层金字塔，得到第t-1帧K层金字塔和第t帧K层金字塔；所述金字塔优选高斯金字塔；所述K优选3～8；

S22，将第t-1帧K层金字塔中的最高层即第K层图像表示为I_t-1,K(x,y)，以I_t-1,K(x,y)为参考图像，将第t帧K层金字塔中的最高层即第K层图像表示为I_t,K(x',y')，以I_t,K(x',y')为待配准图像，采用最小二乘匹配方法计算I_t-1,K(x,y)和I_t,K(x',y')间的仿射变换矩阵；

对参考图像I_t-1,K(x,y)中的坐标(x,y)进行仿射变换：

其中，(x',y')为I_t,K(x',y')的坐标，M为仿射变换矩阵，

m₁～m₆为矩阵参数；设定仿射变换矩阵M的初始值为

根据最小二乘匹配原理，建立误差方程：

解上述误差方程得到m₁～m₆的值；

S23，保留仿射变换矩阵中的m₁～m₄的值不变，对m₅重新赋值为2m₅，对m₆重新赋值为2m₆，利用重新赋值后的仿射变换矩阵M'对第t帧K层金字塔中的第K-1层图像I_t,K-1进行仿射变换，得到变换后的图像I'_t,K-1；

以第t-1帧K层金字塔中的第K-1层图像I_t-1,K-1为参考图像，变换后的图像I'_t,K-1为待配准图像，设定重新赋值后的仿射变换矩阵M'为I_t-1,K-1和I'_t,K-1配准的初始值，采用最小二乘匹配方法计算I_t-1,K-1和I'_t,K-1间的仿射变换矩阵；

S24，重复步骤S23，直到采用最小二乘匹配方法计算第t-1帧K层金字塔中的第1层图像I_t-1,1(x,y)(即I_t-1)和第t帧K层金字塔中的第1层图像I_t,1(x',y')(即I_t)间的仿射变换矩阵，利用仿射变换矩阵对I_t进行变换，得到配准后的第t帧图像

S3、差值图像的显著图计算：将步骤S2得到的

和I_t-1逐像素相减，取绝对值得到差值图像，利用视差注意模型计算差值图像的显著图；

所述利用视差注意模型计算差值图像的显著图，计算步骤如下：

S31，直接将两幅图像

和I_t-1逐像素相减取绝对值，得到差值图像D，并对得到的差值图像D构建L层的高斯金字塔，金字塔的第i层图像表示为D_i，i＝1,…,L；

S32，对于金字塔的第i层图像D_i，通过下式计算D_i的对数幅度谱和相位谱：

L_i(D_i)＝log(||F_i(D_i)||)

Φ_i(D_i)＝ph(F_i(D_i))

其中，

表示傅立叶变换，F_i(D_i)表示D_i的傅里叶变换值，||F_i(D_i)||表示幅度谱，L_i(D_i)表示对数幅度谱，ph(·)用于计算相位谱，Φ_i(D_i)表示相位谱；

S33，根据第i层图像的对数幅度谱L_i(D_i)，计算谱残差R_i(D_i)：

R_i(D_i)＝L_i(D_i)-h(D_i)*L_i(D_i)

其中，h(D_i)为均值滤波器，模板大小优选为3×3像素；

通过下式将谱残差R_i(D_i)和相位谱Φ_i(D_i)进行傅立叶反变换，并用高斯滤波器G₁进行平滑处理，得到第i层图像的显著图：

其中，高斯滤波器G₁的卷积核大小优选为3×3像素；

重复步骤S32与S33，直到得到所有层图像的显著图；

S34，显著图计算：对所有层图像的显著图分别进行插值运算，获得与差值图像D大小相同的插值显著图，然后将所有的插值显著图加权取平均得到平均显著图，最后用高斯滤波器G₂对平均显著图进行高斯平滑，得到差值图像D的显著图；所述高斯滤波器G₂的卷积核半径优选取5～15中的奇数；所述插值运算优选线性插值方法；

S4、运动舰船区域提取：计算步骤S3得到的差值图像D的显著图的均值μ和方差σ，设定阈值T＝μ+λ·σ，根据设定的阈值T对差值图像D的显著图进行阈值分割，获得差值图像D的显著图的二值化图像，其中，系数λ优选为1.0～2.0；进一步对二值化图像依次进行形态学膨胀、连通区域标记、对面积过大或过小的连通区域进行0值填充操作，得到二值化前景图像；根据得到的二值化前景图像分别检测I_t-1和

中的运动舰船区域；

所述进一步对二值化图像依次进行形态学膨胀、连通区域标记、对面积过大或过小的连通区域进行0值填充操作，得到二值化前景图像的方法，具体步骤包括：

S41，利用形态学的膨胀对二值化图像进行形态学膨胀，形态学参数选大小为5×5像素的模板；形态学膨胀的主要作用是把二值化图像中空洞、缺口、分离区域进行填充或连通，以确保运动目标的完整性；

S42，对膨胀处理后的二值化图像进行逐行扫描，找到一个像素值为1、未标记的第一像素点，并标记该点；检查该点的8邻域点，对像素值为1、未标记的邻域点再进行8邻域搜索，通过如此不断迭代的8邻域搜索，逐步标记为整个连通区域；接着再标记下一个未标记的区域，直到所获得的膨胀处理后的二值化图像的所有连通区域都被标记；

S43，由于膨胀处理的二值化图像中还存在一些离散的噪声像素点以及一些虚假舰船区域，因此，统计所有连通区域的面积，设定最大阈值和最小阈值，将连通区域面积大于最大阈值或小于最小阈值的像素值赋值为0，得到二值化前景图像，其中包括若干个像素值为1的连通区域。

所述根据得到的二值化前景图像分别检测I_t-1和

中的运动舰船区域的方法，具体步骤包括：

S44，根据步骤S43得到的二值化前景图像的若干个连通区域，在I_t-1中找到对应坐标位置的若干个图像区域；

S45，针对所述若干个图像区域中的每个区域，利用最大类间方差法(Otsu)对该图像区域进行分割，将分割结果进行第一次膨胀处理，将第一次膨胀处理结果与I_t-1中坐标位置对应的图像区域称为舰船区域；进一步将第一次膨胀处理结果继续第二次膨胀处理，将第二次膨胀处理结果与I_t-1中坐标位置对应的图像区域称为背景及舰船区域；将所述背景及舰船区域与舰船区域的对应像素相减，得到的区域称为背景区域；分别计算舰船区域、背景及舰船区域、背景区域的均值与方差，均值分别记为μ_FG、μ_FG+BG、μ_BG，方差分别记为σ_FG、σ_FG+BG、σ_BG，如果均值与方差满足如下条件：

μ_FG>μ_FG+BG>μ_BG

σ_FG>σ_FG+BG>σ_BG,σ_FG>γ·σ_BG

则将当前图像区域作为运动舰船区域；否则，则将当前图像区域去除；其中，系数γ优选1.5～2.0；

重复步骤S45，当所有的图像区域都处理完，得到I_t-1中的运动舰船区域；

S46，针对图像

根据同样的原理重复步骤S44和S45，得到

中的运动舰船区域。

S5、两帧图像中的运动舰船匹配：将步骤S4得到的I_t-1和

中的运动舰船区域进行匹配，得到运动舰船匹配对；

所述的将步骤S4得到的I_t-1和

中的运动舰船区域进行匹配，得到运动舰船匹配对的方法，具体步骤包括：

S51，依次计算I_t-1中运动舰船区域与

中运动舰船区域的质心坐标，所述质心坐标为舰船区域内所有点的横、纵坐标各自的总和除以区域面积；

S52，对于I_t-1中的任一运动舰船区域，记为S₁；

S53，依次计算S₁与

中的任一运动舰船区域之间的质心距离，当得到的质心距离小于设定的阈值时，则将

中的当前运动舰船区域记为S₂；构造如下比值：

其中，C₁、C₂分别是S₁、S₂所包含的像素个数；

所述阈值一般取10～30像素；

计算S₁的归一化灰度直方图表示为H₁＝{h_1,1,h_1,2,…,h_1,n}，计算S₂的归一化灰度直方图表示为H₂＝{h_2,1,h_2,2,…,h_2,n}，其中，

n为直方图柱的个数，两个灰度直方图之间的巴氏距离为：

构造度量函数ψ：

ψ＝ω·(1-d_B)+(1-ω)·d_count

其中，ω为权重系数，优选ω＝0.5；

如果ψ小于设定的阈值，则将S₂作为S₁的匹配舰船区域；

重复步骤S53，

中的所有运动舰船区域都处理完毕，得到S₁的匹配舰船区域，如果

中存在多个S₁的匹配舰船区域，则选取度量函数ψ最大对应的区域作为S₁的匹配舰船区域，将S₁与其匹配舰船区域称为匹配舰船对。

S54，重复S52～S53，直到I_t-1中的所有运动舰船区域都处理完毕，得到若干匹配舰船对。

S6、计算步骤S2得到的I_t-1和I_t间的仿射变换矩阵的逆矩阵，利用得到的仿射变换矩阵的逆矩阵对运动舰船匹配对对应在

中的运动舰船区域进行变换，得到对应在I_t中的运动舰船区域，则得到第t-1帧图像I_t-1和第t帧图像I_t中的运动舰船匹配对；

S7、舰船关联与跟踪：将t+1赋值给t，重复步骤S1～S6，得到第t帧图像和第t+1帧图像中的运动舰船匹配对，通过第t帧图像实现第t-1帧、第t帧、第t+1帧图像之间的运动舰船关联，从而实现运动舰船的跟踪。

其中，通过第t帧图像实现第t-1帧、第t帧、第t+1帧图像之间的运动舰船关联的方法：

S71，设定

是通过I_t-1和I_t得到的运动舰船匹配对，p＝1,2,…,P，P为I_t-1和I_t中的匹配舰船对数，

是通过I_t和I_t+1得到的运动舰船匹配对，q＝1,2,…,Q，Q为I_t和I_t+1中的匹配舰船对数；

S72，对于通过I_t-1和I_t得到的在I_t中的任一运动舰船区域

通过I_t和I_t+1得到的I_t中的任一运动舰船匹配区域

定义重合系数为：

式中，

表示两个区域的并集，

表示两个区域的交集；如果重合系数R_p,q大于等于设定的阈值，则认为

与

对应同一个舰船，如果重合系数R_p,q小于设定的阈值，则认为

与

不是同一个舰船；

重复步骤S72，直到通过I_t-1和I_t得到的在I_t中的所有运动舰船区域和通过I_t和I_t+1得到的I_t中的所有运动舰船匹配区域都处理完，从而实现了第t-1帧I_t-1、第t帧I_t、第t+1帧I_t+1图像之间的运动舰船关联。所述阈值优选0.6～0.8。

本发明与现有技术相比有益效果为：

(1)本发明为了提取运动舰船区域，基于两帧配准图像的差分图像提取多尺度显著图，可以检测不同大小的运动船只，同时有效避免了差分图像存在的目标空洞、细小轮廓多等问题，确保了运动舰船的完整性。

(2)本发明为了有效地区分运动舰船和虚假区域，利用海面背景和舰船区域的周围灰度对比来确定运动舰船目标，对背景、缓慢光照变化以及背景图像噪声等引起的干扰具有较好的鲁棒性。

(3)本发明对每相邻两帧图像进行舰船匹配，以每三帧图像的中间帧为关联帧，提供一种运动舰船关联的跟踪算法，避免了所有序列图像都与某一固定参考图像配准再进行舰船检测与跟踪的问题。

附图说明

图1是本发明运动舰船检测与跟踪方法的流程图；

图2A为本发明实施例第t-1帧图像；

图2B为本发明实施例配准后的第t帧图像；

图2C为第t-1帧图像与配准后的第t帧的差值图像；

图2D为利用SR方法得到的差值图像的显著图；

图2E为本发明方法计算得到的显著图；

图3A是本发明实施例图2E的显著图的二值化前景图像；

图3B是本发明实施例确定运动舰船区域的示意图；

图4是本发明实施例运动舰船关联示意图。

具体实施方式

本发明实施例给出一种基于卫星序列图像的运动舰船检测与跟踪方法，具体步骤如图1所示，该方法包括以下步骤：

S1、通过卫星不断获取包括运动舰船的卫星序列图像，获取的卫星序列图像中，第t-1帧图像表示为I_t-1，第t帧图像表示为I_t，其中，t为大于等于2的正整数。

S2、以步骤S1中的第t-1帧图像I_t-1为参考图像，以第t帧图像I_t为待配准图像，采用基于金字塔的最小二乘匹配方法计算I_t-1和I_t间的仿射变换矩阵，依据得到的仿射变换矩阵对I_t进行变换得到配准后的第t帧图像

无论静止轨道卫星还是低轨卫星，卫星在成像过程中相机与运动目标都存在运动，因此，对序列图像进行配准是运动舰船检测与跟踪的前提。低轨视频卫星的帧频高(例如，SkySat卫星的帧频为30Hz)，序列图像之间的背景变化不大；静止轨道卫星的帧频虽然低(例如，MOIRE计划的帧频为1Hz)，但是其特有的静止轨道特性，序列图像间的背景变化更小。

具体的，本发明实施例中，考虑到相邻两帧图像之间差别不大，因此对相邻两帧图像的配准建立全局配准模型(参见S.Periaswamy,H.Farid,Elastic Registration in thePresence of Intensity Variations,IEEE Transactions on Medical Imaging,2003,22(7):865-874)，在全局模型基础上，采用最小二乘匹配方法，但是常规的最小二乘法计算量大、计算耗时。本发明为了提高匹配效率采用影像金字塔进行由粗到精的匹配，可大幅节省运算时间。因此，采用基于金字塔的最小二乘匹配方法计算I_t-1和I_t间的仿射变换矩阵，依据得到的仿射变换矩阵对I_t进行变换得到配准后的第t帧图像

的方法如下：

S21，对第t-1帧图像和第t帧图像分别构建K层高斯金字塔，得到第t-1帧K层金字塔和第t帧K层高斯金字塔。

S22，以第t-1帧K层金字塔中的最高层即第K层图像为参考图像，表示为I_t-1,K(x,y)，第t帧K层金字塔中的最高层即第K层图像为待配准图像，表示为I_t,K(x',y')，采用最小二乘匹配方法计算I_t-1,K(x,y)和I_t,K(x',y')间的仿射变换矩阵。

对参考图像I_t-1,K(x,y)中的坐标(x,y)进行仿射变换：

其中，M为仿射变换矩阵，

m₁～m₆为参数；设定仿射变换矩阵的初始值为

根据最小二乘匹配原理，建立误差方程：

解上述误差方程得到m₁～m₆。

S23，保留仿射变换矩阵中的m₁～m₄的值不变，对m₅、m₆重新赋值：m₅＝2m₅、m₆＝2m₆；利用赋值后的仿射变换矩阵对第t帧K层金字塔中的第K-1层图像I_t,K-1进行仿射变换，得到变换后的图像I'_t,K-1。

以第t-1帧K层金字塔中的第K-1层图像I_t-1,K-1为参考图像，变换后的图像I'_t,K-1为待配准图像，设定赋值后的仿射变换矩阵m₁～m₆为I_t-1,K-1和I'_t,K-1配准的初始值，采用最小二乘匹配方法计算I_t-1,K-1和I'_t,K-1间的仿射变换矩阵。

S24，重复步骤S23，直到采用最小二乘匹配方法计算第1层的两幅图像即I_t-1和I_t间的仿射变换矩阵，利用仿射变换矩阵对I_t进行变换，得到配准后的第t帧图像

例如，本发明实施例在基于最小二乘匹配过程中，金字塔层数K取值为5，图2A和图2B所示分别为配准后的第t-1帧图像和第t帧图像，图像大小为1024×768像素。

S3、差值图像的显著图计算：将步骤S2得到的

和I_t-1逐像素相减，取绝对值得到差值图像，利用视差注意模型计算差值图像的显著图。

由于差值方法是以各个像素点作为差分对象，容易存在很多细小的残留噪音。当运动目标在两帧图像中存在重叠时，导致提取出在运动方向上存在空洞、缺口、分离的目标甚至是目标的一部分；当运动目标运动速度过快时，此时就会检测出两个分离的运动目标。因此，本发明引入视觉注意机制，对差值图像进行处理，能够显著增强目标区域与背景区域的对比度，解决差值图像可能存在的运动目标空洞、缺口、分离等问题。利用频域残余谱(Spectral Residual，SR)方法(参见X.Hou,and L.Zhang,Saliency detection:Aspectral residual approach,IEEE Conference Computer Vision and PatternRecognition,2007,pp.1-8.)，将差值图像从时域变换到频域，而后在频域通过分析频谱等特性，最终得到显著性区域图像。

利用视差注意模型计算差值图像的显著图，计算步骤如下：

S31，直接将两幅图像

和I_t-1逐像素相减取绝对值，得到差值图像D，该方法具有算法简单、易于实现，对运动目标的检测灵敏度较高，并对得到的差值图像构建L层的高斯金字塔，金字塔的第i层图像表示为D_i，i＝1,…,L，本发明实施例中，金字塔层数为3。

S32，对于金字塔的第i层图像D_i，通过下式计算D_i的对数幅度谱和相位谱：

L_i(D_i)＝log(||F_i(D_i)||)

Φ_i(D_i)＝ph(F_i(D_i))

其中，

表示傅立叶变换，F_i(D_i)表示D_i的傅里叶变换值，||F_i(D_i)||表示幅度谱，L_i(D_i)表示对数幅度谱，ph(·)用于计算相位谱，Φ_i(D_i)表示相位谱。

S33，根据第i层图像的对数幅度谱L_i(D_i)，计算谱残差R_i(D_i)：

R_i(D_i)＝L_i(D_i)-h(D_i)*L_i(D_i)

其中，h(D_i)为均值滤波器，模板大小优选为取3×3像素。

谱残差能够描述一幅图像中的异常区域，因此可以用来进行显著目标检测。因此，通过下式将谱残差和相位谱进行傅立叶反变换，并用高斯滤波器G₁进行平滑处理，得到第i层图像的显著图：

其中，高斯滤波器G₁的卷积核大小优选为3×3像素。

重复步骤S32与S33，直到得到所有层图像的显著图。

S34，显著图计算：对所有层图像的显著图分别进行插值运算，获得与差值图像大小相同的插值显著图，然后将所有的插值显著图加权取平均得到平均显著图，最后用高斯滤波器G₂对平均显著图进行高斯平滑，得到差值图像的显著图；所述高斯滤波器G₂的卷积核半径优选取5～15中的奇数；所述插值运算优选线性插值方法。

图2C为所述两帧图像的差值图像，图2D为利用SR方法得到的差值图像的显著图，图2E为本发明方法计算得到的显著图，可以看出，本发明方法处理后变化区域得到明显加强，对后续舰船提取提供基础。

S4、运动舰船区域提取：计算步骤S3得到的差值图像的显著图的均值μ和方差σ，设定阈值T＝μ+λ·σ，根据设定的阈值T对差值图像的显著图进行阈值分割，获得差值图像的显著图的二值化图像，其中，系数λ优选为1.0～2.0；进一步对二值化图像依次进行形态学膨胀、连通区域标记、对面积过大或过小的连通区域进行0值填充操作，得到二值化前景图像；根据得到的二值化前景图像分别检测I_t-1和

中的运动舰船区域。

其中，进一步对二值化图像进行形态学膨胀、连通区域标记、对面积过大或过小的连通区域进行0值填充等操作，得到二值化前景图像，具体步骤包括：

S41，利用形态学的膨胀对二值化图像进行形态学膨胀，形态学参数选大小为5×5像素的模板；形态学膨胀的主要作用是把二值化图像中空洞、缺口、分离区域进行填充或连通，以确保运动目标的完整性；

S42，对膨胀处理后的二值化图像进行逐行扫描，找到一个像素值为1、未标记的第一像素点，并标记该点；检查该点的8邻域点，对像素值为1、未标记的邻域点再进行8邻域搜索，通过如此不断迭代的8邻域搜索，逐步标记为整个连通区域；接着再标记下一个未标记的区域，直到所获得的膨胀处理后的二值化图像的所有连通区域都被标记；

S43，由于膨胀处理的二值化图像中还存在一些离散的噪声像素点以及一些虚假舰船区域，因此，统计所有连通区域的面积，设定最大阈值和最小阈值，将连通区域面积大于最大阈值或小于最小阈值的像素值赋值为0，得到二值化前景图像，其中包括若干个像素值为1的连通区域。

例如，请参见图3A，其为图2E所示的显著图经过形态学膨胀、去除面积过大或过小的连通区域后得到的二值化前景图像。

具体的，本发明实施例根据海洋背景区域与舰船区域的灰度差别，进一步判定是否为运动舰船区域。根据得到的二值化前景图像分别检测I_t-1和

中的运动舰船区域的方法，具体步骤包括：

S44，根据步骤S43得到的二值化前景图像的若干个连通区域，在I_t-1中找到对应坐标位置的若干个图像区域；

S45，针对所述若干个图像区域中的每个区域(图3B[a])，利用最大类间方差法(Otsu)(参见N.Otsu,A threshold selection method from gray-level histograms,IEEE Transactions on Systems,Man,and Cybernetics,vol.9,no.1,pp.62–66,1979.)对该图像区域进行分割(图3B[b])，将分割结果进行第一次膨胀处理(模板大小3×3～5×5像素)(图3B[c])，进一步将第一次膨胀处理结果继续第二次膨胀处理(图3B[d])；将第一次膨胀处理结果与第t-1帧图像中坐标位置对应的图像区域称为舰船区域(图3B[e])，将第二次膨胀处理结果与第t-1帧图像中坐标位置对应的图像区域称为背景及舰船区域(图3B[f])；将所述背景及舰船区域与舰船区域的对应像素相减，得到的区域称为背景区域(图3B[g])；分别计算舰船区域、背景及舰船区域、背景区域的均值与方差，均值分别记为μ_FG、μ_FG+BG、μ_BG，方差分别记为σ_FG、σ_FG+BG、σ_BG，如果均值与方差满足如下条件：

μ_FG>μ_FG+BG>μ_BG

σ_FG>σ_FG+BG>σ_BG,σ_FG>γσ_BG

则将当前图像区域作为运动舰船区域；否则，则去除当前图像区域；其中，根据经验值，系数γ优选1.5～2.0。

重复步骤S45，当所有的图像区域都处理完，得到I_t-1中的运动舰船区域。

其中，Otsu被认为是阈值自动选取方法的最优方法之一，它主要依据图像的灰度特征将图像分割为背景和目标两部分。本发明实施例中的运动舰船区域与背景区域灰度特性差异明显，因此，可用Otsu方法将每个图像区域中的舰船与背景区分开来。设定图像区域的灰度级为L，灰度为i的像素数为n_i，图像区域的总像素数为total。最大类间方差的计算公式为：

其中，T是将图像区域分为舰船与背景两个区域的阈值，ω_ship、ω_bg分别是两个区域的像素数占图像的比例，μ_ship、μ_bg分别是两个区域的灰度均值，μ为图像区域的灰度均值。由于图像区域中的像素级数少，因此本发明实施例计算最大类间方差时只遍历图像区域存在的灰度级，这样可以进一步减少计算复杂性，提高图像区域分割速度。

S46，针对图像

根据同样的原理重复步骤S44和S45，得到

中的运动舰船区域。

S5、两帧图像中的运动舰船匹配：将步骤S4得到的I_t-1和

中的运动舰船区域进行匹配，得到运动舰船匹配对。

S51，依次计算I_t-1中运动舰船区域与

中运动舰船区域的质心坐标，所述质心坐标为舰船区域内所有点的横、纵坐标各自的总和除以区域面积。

S52，对于I_t-1中的任一运动舰船区域，记为S₁。

S53，依次计算S₁与

中的任一运动舰船区域之间的质心距离，当得到的质心距离小于设定的阈值时，则将

中的当前运动舰船区域记为S₂；构造如下比值：

其中，C₁、C₂分别是S₁、S₂所包含的像素个数；

所述阈值一般取10～30像素；

计算S₁的归一化灰度直方图表示为H₁＝{h_1,1,h_1,2,…,h_1,n}，计算S₂的归一化灰度直方图表示为H₂＝{h_2,1,h_2,2,…,h_2,n}，其中，

n为直方图柱的个数，两个灰度直方图之间的巴氏距离(参见D.Comaniciu,V.Ramesh,and P.Meer,Kernel-Based Object Tracking,IEEE Trans.Pattern Anal.Mach.Intell.,25(5):564-577,2003.)为：

巴氏距离越小，两个区域的相似度就越大，反之就越小。

构造度量函数ψ：

ψ＝ω·(1-d_B)+(1-ω)·d_count

其中，ω为权重系数，本发明实施例选取ω＝0.5。

如果ψ小于设定的阈值，则将S₂作为S₁的匹配舰船区域；

重复步骤S53，

中的所有运动舰船区域都处理完毕，得到S₁的匹配舰船区域，如果

中存在多个S₁的匹配舰船区域，则选取度量函数ψ最大对应的区域作为S₁的匹配舰船区域，将S₁与其匹配舰船区域称为匹配舰船对。

S54，重复S52～S53，直到I_t-1中的所有运动舰船区域都处理完毕，得到若干匹配舰船对。

S6、计算步骤S2得到的I_t-1和I_t间的仿射变换矩阵的逆矩阵，利用得到的仿射变换矩阵的逆矩阵对运动舰船匹配对对应在

中的运动舰船区域进行变换，得到对应在I_t中的运动舰船区域，则得到第t-1帧图像I_t-1和第t帧图像I_t中的运动舰船匹配对；

S7、舰船关联与跟踪：将t+1赋值给t，重复步骤S1～S6，得到第t帧图像和第t+1帧图像中的运动舰船匹配对，通过第t帧图像实现第t-1帧、第t帧、第t+1帧图像之间的运动舰船关联，从而实现运动舰船的跟踪。

其中，通过第t帧图像实现第t-1帧、第t帧、第t+1帧图像之间的运动舰船关联的方法：

S71，设定

是通过I_t-1和I_t得到的运动舰船匹配对，p＝1,2,…,P，P为I_t-1和I_t中的匹配舰船对数，

是通过I_t和I_t+1得到的运动舰船匹配对，q＝1,2,…,Q，Q为I_t和I_t+1中的匹配舰船对数；

S72，对于通过I_t-1和I_t得到的在I_t中的任一运动舰船区域

通过I_t和I_t+1得到的I_t中的任一运动舰船匹配区域

定义重合系数为：

式中，

表示两个区域的并集，

表示两个区域的交集。

一般来说，如果所述第i个舰船区域与第j个舰船区域是同一个舰船，则重合系数理论上应该接近1，但是考虑到实际情况，很难达到1，因此，如果重合系数R_p,q大于设定的阈值，则认为

与

对应同一个舰船，所述阈值取值0.6～0.8，从而实现了第t-1帧I_t-1、第t帧I_t、第t+1帧I_t+1图像之间的运动舰船关联。

如图4所示，第t-1帧图像与第t帧图像中得到两对匹配的运动舰船(参见图4[a]与4[b1]中的S1与S2)，第t帧图像与第t+1帧图像中得到三对匹配的运动舰船(参见图4[b2]与4[c]的S3、S4与S5)。因此，可以通过第t帧建立关联，图4[b1]中的S1与图4[b2]中的S3属于同一舰船，图4[b1]中的S2与图4[b2]中的S4属于同一舰船，图4[b2]中的S3是第t帧图像中新出现的舰船。

Claims (9)

1.一种基于卫星序列图像的运动舰船检测与跟踪方法，其特征在于该方法的步骤包括：

S1、通过卫星不断获取包括运动舰船的卫星序列图像；

S2、对步骤S1中获取的包括运动舰船的卫星序列图像中相邻的两幅原始图像进行配准，得到配准后的两幅图像；

S3、计算步骤S2配准后的两幅图像的差值图像，并计算差值图像的显著图；

S4、根据步骤S3得到的差值图像的显著图提取两幅图像中的运动舰船区域；

S5、将步骤S4得到的两幅图像中的运动舰船区域进行匹配，得到运动舰船匹配对；

S6、将步骤S5得到的运动舰船匹配对变换到步骤S2中原始两幅图像中，得到原始两幅图像中的运动舰船匹配对；

S7、对步骤S1中获取的包括运动舰船的卫星序列图像中相邻的两幅原始图像进行配准，得到配准后的两幅图像；

所述步骤S2中的相邻的两幅原始图像中有且仅有一幅为所述步骤S7中的相邻的两幅原始图像中的一幅；

S8、根据步骤S3-S6，对S7中原始两幅图像中的运动舰船匹配对；

S9、对步骤S6得到的运动舰船匹配对和步骤S8得到的运动舰船匹配对进行关联，完成运动舰船的检测与跟踪；

所述的步骤S9中，对步骤S6得到的运动舰船匹配对和步骤S8得到的运动舰船匹配对进行关联的方法为：

定义步骤S6得到的运动舰船匹配对为第t-1帧和第t帧得到的运动舰船匹配对；步骤S8得到的运动舰船匹配对为第t帧和第t+1帧得到的运动舰船匹配对，定义第t+1帧的图像表示为I_t+1；

即通过第t帧图像实现第t-1帧、第t帧、第t+1帧图像之间的运动舰船关联的方法：

S91，设定

是通过I_t-1和I_t得到的运动舰船匹配对，p＝1,2,…,P，P为I_t-1和I_t中的匹配舰船对数，

是通过I_t和I_t+1得到的运动舰船匹配对，q＝1,2,…,Q，Q为I_t和I_t+1中的匹配舰船对数；

S92，对于通过I_t-1和I_t得到的在I_t中的任一运动舰船区域

通过I_t和I_t+1得到的I_t中的任一运动舰船匹配区域

定义重合系数为：

式中，

表示两个区域的并集，

表示两个区域的交集；如果重合系数R_p,q大于等于设定的阈值，则认为

与

对应同一个舰船，如果重合系数R_p,q小于设定的阈值，则认为

与

不是同一个舰船；

重复步骤S92，直到通过I_t-1和I_t得到的在I_t中的所有运动舰船区域和通过I_t和I_t+1得到的I_t中的所有运动舰船匹配区域都处理完，从而实现了第t-1帧I_t-1、第t帧I_t、第t+1帧I_t+1图像之间的运动舰船关联。

2.根据权利要求1所述的一种基于卫星序列图像的运动舰船检测与跟踪方法，其特征在于：所述的步骤S2中，相邻的两幅原始图像进行配准的方法为：

定义两幅原始图像分别为I_t-1和I_t，其中I_t-1为第t-1帧图像，I_t为第t帧图像；t为大于等于2的正整数；

以第t-1帧图像I_t-1为参考图像，以第t帧图像I_t为待配准图像，采用基于金字塔的最小二乘匹配方法计算I_t-1和I_t间的仿射变换矩阵，依据得到的仿射变换矩阵对I_t进行变换得到配准后的第t帧图像

3.根据权利要求2所述的一种基于卫星序列图像的运动舰船检测与跟踪方法，其特征在于：所述的采用基于金字塔的最小二乘匹配方法计算I_t-1和I_t间的仿射变换矩阵，依据得到的仿射变换矩阵对I_t进行变换得到配准后的第t帧图像

的方法如下：

S21，对第t-1帧图像和第t帧图像分别构建K层金字塔，得到第t-1帧K层金字塔和第t帧K层金字塔；

S22，将第t-1帧K层金字塔中的最高层即第K层图像表示为I_t-1,K(x,y)，以I_t-1,K(x,y)为参考图像，将第t帧K层金字塔中的最高层即第K层图像表示为I_t,K(x',y')，以I_t,K(x',y')为待配准图像，采用最小二乘匹配方法计算I_t-1,K(x,y)和I_t,K(x',y')间的仿射变换矩阵；

对参考图像I_t-1,K(x,y)中的坐标(x,y)进行仿射变换：

其中，(x',y')为I_t,K(x',y')的坐标，M为仿射变换矩阵，

m₁～m₆为矩阵参数；设定仿射变换矩阵M的初始值为

根据最小二乘匹配原理，建立误差方程：

解上述误差方程得到m₁～m₆的值；

S23，保留仿射变换矩阵中的m₁～m₄的值不变，对m₅重新赋值为2m₅，对m₆重新赋值为2m₆，利用重新赋值后的仿射变换矩阵M'对第t帧K层金字塔中的第K-1层图像I_t,K-1(x,y)进行仿射变换，得到变换后的图像I'_t,K-1(x,y)；

以第t-1帧K层金字塔中的第K-1层图像I_t-1,K-1(x,y)为参考图像，变换后的图像I'_t,K-1(x,y)为待配准图像，设定重新赋值后的仿射变换矩阵M'为I_t-1,K-1 (x,y)和I'_t,K-1(x,y)配准的初始值，采用最小二乘匹配方法计算I_t-1,K-1(x,y)和I'_t,K-1(x,y)间的仿射变换矩阵；

S24，重复步骤S23，直到采用最小二乘匹配方法计算第t-1帧K层金字塔中的第1层图像I_t-1,1(x,y)即I_t-1和第t帧K层金字塔中的第1层图像I_t,1(x',y')即I_t间的仿射变换矩阵，利用仿射变换矩阵对I_t进行变换，得到配准后的第t帧图像

4.根据权利要求1所述的一种基于卫星序列图像的运动舰船检测与跟踪方法，其特征在于：所述的步骤S3中，采用视差注意模型计算差值图像的显著图。

5.根据权利要求4所述的一种基于卫星序列图像的运动舰船检测与跟踪方法，其特征在于：所述采 用视差注意模型计算差值图像的显著图，计算步骤如下：

S31，直接将两幅图像

和I_t-1逐像素相减取绝对值，得到差值图像D，并对得到的差值图像D构建L层的高斯金字塔，金字塔的第i层图像表示为D_i，i＝1,…,L；

S32，对于金字塔的第i层图像D_i，通过下式计算D_i的对数幅度谱和相位谱：

L_i(D_i)＝log(||F_i(D_i)||)

Φ_i(D_i)＝ph(F_i(D_i))

其中，

表示傅立叶变换，F_i(D_i)表示D_i的傅里叶变换值，||F_i(D_i)||表示幅度谱，L_i(D_i)表示对数幅度谱，ph(·)用于计算相位谱，Φ_i(D_i)表示相位谱；

S33，根据第i层图像的对数幅度谱L_i(D_i)，计算谱残差R_i(D_i)：

R_i(D_i)＝L_i(D_i)-h(D_i)*L_i(D_i)

其中，h(D_i)为均值滤波器；

通过下式将谱残差R_i(D_i)和相位谱Φ_i(D_i)进行傅立叶反变换，并用高斯滤波器G₁进行平滑处理，得到第i层图像的显著图：

重复步骤S32与S33，直到得到所有层图像的显著图；

S34，显著图计算：对所有层图像的显著图分别进行插值运算，获得与差值图像D大小相同的插值显著图，然后将所有的插值显著图加权取平均得到平均显著图，最后用高斯滤波器G₂对平均显著图进行高斯平滑，得到差值图像D的显著图。

6.根据权利要求1所述的一种基于卫星序列图像的运动舰船检测与跟踪方法，其特征在于：所述的步骤S4中，提取两幅图像中的运动舰船区域的方法为：计算差值图像的显著图的均值μ和方差σ，设定阈值T＝μ+λ·σ，系数λ为1.0～2.0；根据设定的阈值T对差值图像D的显著图进行阈值分割，获得差值图像D的显著图的二值化图像，进一步对二值化图像依次进行形态学膨胀、连通区域标记、将连通区域面积大于最大阈值或小于最小阈值的像素值赋值为0，得到二值化前景图像；根据得到的二值化前景图像分别检测I_t-1和

中的运动舰船区域；所述的最大阈值和最小阈值均为设定值。

7.根据权利要求6所述的一种基于卫星序列图像的运动舰船检测与跟踪方法，其特征在于：所述根据得到的二值化前景图像分别检测I_t-1和

中的运动舰船区域的方法，具体步骤包括：

S44，根据二值化前景图像的若干个连通区域，在I_t-1中找到对应坐标位置的若干个图像区域；

S45，针对所述若干个图像区域中的每个区域，利用最大类间方差法对该图像区域进行分割，将分割结果进行第一次膨胀处理，将第一次膨胀处理结果与I_t-1中坐标位置对应的图像区域称为舰船区域；进一步将第一次膨胀处理结果继续第二次膨胀处理，将第二次膨胀处理结果与I_t-1中坐标位置对应的图像区域称为背景及舰船区域；将所述背景及舰船区域与舰船区域的对应像素相减，得到的区域称为背景区域；分别计算舰船区域、背景及舰船区域、背景区域的均值与方差，均值分别记为μ_FG、μ_FG+BG、μ_BG，方差分别记为σ_FG、σ_FG+BG、σ_BG，如果均值与方差满足如下条件：

μ_FG>μ_FG+BG>μ_BG

σ_FG>σ_FG+BG>σ_BG,σ_FG>γ·σ_BG

系数γ为1.5～2.0；

则将当前图像区域作为运动舰船区域；否则，则将当前图像区域去除；

重复步骤S45，当所有的图像区域都处理完，得到I_t-1中的运动舰船区域；

S46，针对图像

根据同样的原理重复步骤S44和S45，得到

中的运动舰船区域。

8.根据权利要求1所述的一种基于卫星序列图像的运动舰船检测与跟踪方法，其特征在于：所述的步骤S5中，得到运动舰船匹配对的方法，具体步骤包括：

S51，依次计算I_t-1中运动舰船区域与

中运动舰船区域的质心坐标，所述质心坐标为舰船区域内所有点的横、纵坐标各自的总和除以区域面积；

S52，对于I_t-1中的任一运动舰船区域，记为S₁；

S53，依次计算S₁与

中的任一运动舰船区域之间的质心距离，当得到的质心距离小于设定的阈值时，则将

中的当前运动舰船区域记为S₂；构造如下比值：

其中，C₁、C₂分别是S₁、S₂所包含的像素个数；

计算S₁的归一化灰度直方图表示为H₁＝{h_1,1,h_1,2,…,h_1,n}，计算S₂的归一化灰度直方图表示为H₂＝{h_2,1,h_2,2,…,h_2,n}，其中，

n为直方图柱的个数，两个灰度直方图之间的巴氏距离为：

构造度量函数ψ：

ψ＝ω·(1-d_B)+(1-ω)·d_count

其中，ω为权重系数；

如果ψ小于设定的阈值，则将S₂作为S₁的匹配舰船区域；

重复步骤S53，

中的所有运动舰船区域都处理完毕，得到S₁的匹配舰船区域，如果

中存在多个S₁的匹配舰船区域，则选取度量函数ψ最大对应的区域作为S₁的匹配舰船区域，将S₁与其匹配舰船区域称为匹配舰船对；

S54，重复S52～S53，直到I_t-1中的所有运动舰船区域都处理完毕，得到若干匹配舰船对。

9.根据权利要求1所述的一种基于卫星序列图像的运动舰船检测与跟踪方法，其特征在于：所述的步骤S6中，将步骤S5得到的运动舰船匹配对变换到步骤S2中原始两幅图像中的方法为：

计算I_t-1和I_t间的仿射变换矩阵的逆矩阵，利用得到的仿射变换矩阵的逆矩阵对运动舰船匹配对对应在

中的运动舰船区域进行变换，得到对应在I_t中的运动舰船区域，则得到第t-1帧图像I_t-1和第t帧图像I_t中的运动舰船匹配对。

Images