基于MSA 特征和模拟退火优化的遥感图像多目标关联算法

李晖晖，滑 立，杨 宁，刘 坤

（1.西北工业大学 自动化学院，西安710072；2.上海海事大学 信息工程学院，上海200135）

0 引 言

多源遥感图像融合［1］能够综合利用多源互补的图像信息以提高信息的准确性与可靠性，因此在诸多领域得到广泛应用，包括军事、测绘、地质、农业及灾情监测等。而目标关联是融合的先决条件，概括来说就是在两幅（或多幅）不同时间（或空间）获取的遥感图像中，将来自同一目标的信息进行匹配。然后才能进行融合检测、识别、跟踪等一系列处理。

以往的目标关联方法主要指状态滤波类方法，将目标视为点对象，利用雷达型数据提供的目标位置、速度、方位等运动特征进行关联，适合密集采样的序列图像。如文献［2］提出了基于目标质心及质心偏移量的联合概率数据关联（Joint probabilistic data association，JPDA）方法，首先将目标与背景分割，然后计算目标区域的质心，以估计目标的运动信息，最后利用JPDA 方法［3］实现关联与跟踪。文献［4］提出了在红外系统中将目标检测和跟踪联合设计的方法，利用多假设跟踪（Multiple hypothesis tracking，MHT）算法进行目标图像特征点的匹配跟踪。但由于当前遥感成像技术一般只能获取采样稀疏的遥感图像，很难预测目标的状态量，因此信息融合领域中传统的利用状态特征进行关联的方法并不适合遥感图像的目标关联。需要利用目标图像自身特征进行匹配，来建立新的目标关联准则。文献［5］结合Gabor小波提取特征以及迭代最小平方和的方法，实现目标尺度变化、旋转及平移等几何变形时的跟踪。文献［6］首先提取目标图像的8个不变矩特征，然后通过训练神经网络实现目标的匹配和定位。

首先需要解决图像特征的提取问题。不变矩方法是解决图像特征不变性的常用方法，它能够克服视点变化对特征量的干扰，最具代表性的就是Hu 矩 特 征［7］。2005 年Rahtu 等［8］结 合 多 尺度几何分析方法，在Hu矩特征的基础上，研究了可以同时捕获空间特性和图像强度的多维不变描述子，构造了多尺度自卷积矩（Multi scale autoconvolution，MSA）特征，实践证明这是目前最稳健的不变矩特征之一［8］。然后采用距离度量法进行特征匹配，但由于遥感信息的不确定性和特征提取算法的不精确性导致特征匹配结果存在误差，因此需要进行关联修正来消除多目标对应关系的模糊性。文献［9］对上述问题提供了一种解决思路，首先利用目标特征匹配结果构造一个多目标关联代价矩阵（Association cost matrix，ACM），然后采用模拟退火算法求解使整体关联代价最小的关联代价矩阵，即为最终多目标关联结果。

本文借鉴文献［9］的思想，将基于MSA 特征匹配和关联代价矩阵最优化的多目标关联算法应用于海地机场飞机多目标关联问题上，在组合优化过程中，着重对模拟退火算法进行了改进：在原始算法的搜索过程中设定内、外循环迭代次数，并设计了一种新的温度更新函数，它具有一定的自适应性，改进了温度控制方式，提高了搜索效率。通过海地机场6类飞机灰度图像进行多组实验，结果证明了本文算法具有很好的多目标关联效果。

1 多尺度自卷积矩（MSA）

多尺度自卷积矩（MSA）是迄今为止最为有效的不变矩。MSA 变换的基本原理：图像中任意一点的坐标可以用其他随机不共线三点的坐标线性表示，对图像进行仿射变换时，这四个点线性表示的系数不会发生变化，因此可利用系数不变性质构造仿射不变量，具有较好的稳定性。

定义图像的仿射变换为A＝A（T，t），设图像仿射变换前后的坐标分别为x 和x′，即x′＝A（x）＝Tx＋t。设f 为图像的灰度函数，仿射变换后的图像为f′（x′）＝f（Tx＋t）。

设x0，x1，x2∈R2是f（x）定义域中的3个点，则点μα，β 可线性表示为：

μ′α，β可表示为：

式中：x′0，x′1，x′2分别 对应于x0，x1，x2的仿射变换。可见随机变量f（μα，β）和f′（μ′α，β）有着相同的分布，因此其数学期望也相等，这种相等关系独立于仿射变换。MSA 方法即为变量f（μα，β）的数学期望值：

单纯从数学上考虑积分元素在积分域内的任意性，可知F（α，β）就是仿射变换的不变特征。设γ＝1－α－β，则μα，β ＝αx1＋βx2＋γx0。概率密度函数pUα，β（u）＝（pα×pβ×pγ）（u），可将式（3）改写为：

2 关联代价矩阵（ACM）

2.1 构造ACM

研究两幅遥感图像中的目标关联问题，可将其看作2维分配问题［10］。假设前一图像中有M个目标，后一图像中有N 个目标。现在定义一个二值分配变量：

于是可用此形成整体关联矩阵：a ＝｛a（m，n）；m＝0，1，…，M；n ＝0，1，…，N｝，关联的目的是找到最优关联矩阵使得下面的全局关联代价最小：

提取出各目标的MSA 特征后，通过计算目标特征之间的欧氏距离可获得任意两目标的MSA 特征匹配代价m ＝0，1，…，M；n＝0，1，…，N，其中为前一幅图像中目标m 的k 维MSA 特征，为后一幅图像中的k 维MSA 特征。C ＝｛a（m，n）c（m，n）；（m ＝0，1，…，M）；（n ＝0，1，…，N）｝为关联代价矩阵。当m ＝n＝0时，表示两幅图像之间没有发生匹配，c（0，0）可设为较大的正数。当m或n 有一个为0时，表示存在漏检或虚警等问题，此时令

2.2 ACM 最优化

假设同一幅图像中目标种类各不相同，根据多目标关联的具体应用，优化建模需要满足如下条件：两幅图像最多只能分配给对方一个目标。数学表述即为“关联矩阵a 中每一行和每一列只有一个值为1”，以此构造可行的关联矩阵。设C＝｛C1，C2，…，Ck｝为所有可行的关联代价矩阵的集合，E（Ci）为对应于Ci的目标函数值。构造关联代价矩阵最优化的数学模型为：E（C＊）＝minE（Ci），C＊∈C，∀Ci∈C。目标函数的最小值对应关联代价矩阵的最优解，即关联代价矩阵的能量最小［11］。为方便表述，记a （m ，n） 为amn，考虑到实际的约束条件，进一步将目标函数定义如下：

式中：前两项为惩罚项，对应具体的约束条件，只有当前两项均为零时，第三项才是多目标关联的实际代价。不考虑关联为空的情况，设定A、B 为较大的正数，只有当关联矩阵a 中每一行和每一列有且只有一个元素为1、其余元素均为零时，求得的目标函数值才会最小，对应的关联矩阵即为最优关联结果。

3 模拟退火算法

3.1 传统模拟退火算法

在实际优化问题中，有些目标函数具有非凸性，存在局部最优解，尤其是当优化问题规模增大时，迅速增加的局部最优解数目将会导致全局优化结果不够准确，同时又增加了计算量。模拟退火 算 法［12］（Simulated annealing，SA）是 基 于Monte－Carlo迭代求解策略的一种随机寻优算法，它来源于对固体退火过程的模拟，挖掘物理工程中固体物质的退火降温过程与组合优化问题的相似性，从理论上来说，它是一种近似全局最优算法。1983年，Krickpatric等人将模拟退火算法成功应用于大规模组合优化问题的求解［12］。

模拟退火过程大致分为以下三个阶段：

（1）升温阶段：加热时，固体内部粒子随着温度的升高，运动增强，当温度足够高，固体变为无序状，内能增大。升温阶段对应算法的设定初始温度（充分大）。

（2）平衡阶段：退火过程中，根据热力学定律可知在每一温度下，系统自由能逐步减到最小时，即达到平衡态。等温下热平衡过程可用Monte Carlo方法模拟，为使结果精确，需要大量采样，导致计算量很大。为提高搜索效率，1953 年，Metropolis等人提出重要性采样法，粒子在当前温度T 趋于平衡态的概率为χ ＝exp（－ΔE／kT），其 中 ΔE 为 内 能 改 变 量，k 为Boltzmann常数。依据概率χ ＞r＝random［0，1］接受该状态为重要状态。上述接受新状态的准则称Metropolis准则［13］，相对Monte Carlo方法来说，计算量明显减少。

（3）冷却阶段：随着温度缓慢降低，固体内部粒子运动渐趋有序，同时能量下降。本阶段，采用合适的温度衰减函数即退火策略来控制算法进程。最简单也最常用的控制参数衰减函数为：Tk＝aTk－1＝akT0（k＝1，2，…），其中，a 是常值，取0.5～0.9。重复执行新状态产生及判断是否接受新状态的过程，逐渐降温直到趋于零时，达到能量最小状态，即求得全局最优解。

3.2 改进模拟退火算法寻优

模拟退火算法是一种启发式搜索方法，它具有原理简单、适合求解组合优化问题中的全局（或近似）最优解［14］、适用范围广等优点，但也存在一些缺点：当变量过多、构造的目标函数较复杂时，为获得较精确的全局最优解，初始温度需要足够高且降温幅度要小，这将导致迭代搜索过程缓慢，求解时间太长，工业应用较困难；在Metropolis准则依概率接受新状态环节中，有可能丢失当前最优解而陷入局部最优解邻域。

针对上述问题，为确保算法的优化质量，同时又提高算法的搜索效率（时间性能），本文研究了以下改进算法：选择合适的初始温度和终止温度，并设计高效的退火策略，改进对温度的控制方式，设计一个新的温度更新函数，根据某一温度下状态被接受的次数来决定降温幅度，保证温度更新有一定的自适应性。在搜索过程中设置双阈值，即内、外循环的迭代次数，若目标函数值在当前温度下保持s＿max步不变，记为当前最优解；判断在T （i ＋1） 温度下求得的关联矩阵与上一温度 T （）i 下求得的关联矩阵是否一致，若外循环连续iter＿max步降温过程中搜索到的最优关联矩阵均保持不变，可认为得到了较高质量的全局最优解。具体步骤如下：

4 实验结果及分析

为了检验所提方法的可行性，选用如图1所示的海地机场6 类飞机目标灰度图像进行实验（这些图像来自于某些海地机场的IKONOS卫星图像，如图2 所示，它的分辨率为1 m，大小为3000×3000像素），每幅图像大小为128×128像素。考虑到实际图像会受到噪声干扰、遮挡、亮度变化等情况的影响，我们除了对图像进行尺度缩放、旋转、仿射变换外，还对图像加入了不同等级的高斯白噪声、遮挡以及亮度变化。构成了多组6类待关联的飞机目标，下面给出了原始飞机目标和其中6组仿真图像，如图3所示。

图1 海地机场6类飞机目标Fig.1 Six class aircraft targets of Haiti airport

图2 IKONOS卫星图像（海地机场）Fig.2 IKONOS satellite image（Haiti airport）

图3 6组仿真图像Fig.3 Six group simulation image

以图2中的飞机目标和图3中其仿真的待关联目标进行关联为例：首先在集合｛－1，－0.75，－0.5，－0.25，0，0.25，0.5，0.75，1｝中取任意两元素的组合构成29对（α，β）值，进而计算各飞机目标的29维MSA 特征向量，通过距离度量法获得两两匹配概率；随机初始化关联矩阵a，并求得相应的关联代价矩阵C 如下：

然后构造整体关联代价矩阵目标函数，利用改进的模拟退火算法寻优，得到使得ACM 能量最小的关联矩阵。此外，提取出MSA 特征后，利用最近邻（NN）算法进行关联实验，与得到的最终关联矩阵作比较。由于仿真图像目标组与原始目标组的排列顺序一样，只对每个目标图像进行变换、加入干扰等，理论上求得的最终关联矩阵应该是单位矩阵。实验结果与理论一致，可见利用关联代价矩阵最优化能够克服多目标特征匹配引起的关联模糊性，关联结果更加准确。

引入关联正确率指标进行统计实验分析，图4给出了结合MSA 特征和ACM 最优化的多目标关联结果，图5 给出了结合MSA 特征和NN方法实现多目标关联的结果。可以看出，本文基于MSA 和ACM 最优化的多目标关联算法结果明显优于利用最近邻算法进行关联的结果。当图像发生尺度缩放、旋转、噪声干扰、仿射变换、遮挡以及亮度变化等一系列变化时，最近邻算法无法解决关联模糊性，而ACM 考虑了整体关联代价，具有一定的抗模糊性，可获得更为有效的关联结果。

图4 MSA＋ACMFig.4 Association accurate ratio of MSA＋ACM

图5 MSA＋NNFig.5 Association accurate ratio of MSA＋NN

［9］是用MSA 特征和ACM 最优化算法解决舰船目标的关联问题，本文借鉴其思路，将该算法应用于飞机图像的多目标关联问题上，并着重对传统模拟退火算法进行了改进。利用本文改进的模拟退火算法和传统算法对每组目标群图像分别独立运行50次，记录每组运行时间（s）和目标关联结果数据，计算关联正确率（%），进行统计分析，其对比结果如表1所示。可以看出，改进的模拟退火算法既保持了多目标关联结果的准确性，又减少了计算量，在时间性能（搜索效率）上有了较大的提高。

表1 传统SA和改进的SA实验结果对比Table 1 Experiment comparison results of the traditional SA and improved SA

5 结束语

针对无法准确估计采样稀疏的遥感图像中目标的状态信息，而利用图像特征匹配的目标关联算法又无法处理大场景中多个目标关联引起的模糊性的问题。本文研究了基于MSA 特征和关联代价矩阵最优化的遥感图像多目标关联算法，并着重改进了模拟退火算法寻优过程。通过仿真实验表明：MSA 特征可克服遥感成像中视点变化或目标姿态变化等因素的影响，具有良好的仿射不变性和抗干扰性，是有效的关联量；在关联准则的约束下构造目标函数对ACM 进行最优化，得到的全局最优解有效消除了多目标之间特征关联的模糊性；改进的模拟退火算法能够改善执行多次迭代搜索过程的时效性，可快速得到全局最优解。

参考文献：

［1］雷琳.多源遥感图像舰船目标特征提取和融合技术研究［D］.长沙：国防科技大学电子科学与工程学院，2008.Lei Lin.Ship feature extraction and fusion in multiple remote sensing images［D］.Changsha：School of Electronic Science and Engineering，National University of Defense Technology，2008.

［2］Bar－shalom Y，Sherlukde H M，Pattipati K R.Use of measurements from an image sensor for precision target tracking［J］.IEEE Trans AES，1989，25（6）：863－871.

［3］Bar Shalom Y.Extension of the probabilistic data association filter in multi－target tracking［C］∥Proceedings of the 5th symposium on nonlinear estimation，1974：16－21.

［4］Blackman S，Dempster R J，Broida T J.Multiple hypothesis track confirmation for infrared surveillance systems［J］.IEEE Trans AES，1993，29（3）：810－823.

［5］董学志，宋建中，韩广良.一种利用Gabor小波特征的目标跟踪方法［J］.光学技术，2003，29（4）：484－486.Dong Xue－zhi，Song Jian－zhong，Han Guang－liang.Target tracking method using wavelet feature［J］.Optical Technique，2003，29（4）：484－486.

［6］姚剑，刘其真，张斌.模糊技术与神经网络的混合算法在运动目标识别与跟踪中的应用［J］.计算机工程与应用，2000，1：62－64.Yao Jian，Liu Qi－zhen，Zhang Bin.The application of combination of fuzzy algorithms and neural network in tracking moving object［J］.Computer Engineering and Applications，2000，1：62－64.

［7］Hu M K.Pattern recognition by moment invariant［C］∥Proc IRE，1961，49：1428－1436.

［8］Rahtu E，Salo M，Heikkila J.Affine invariant pattern recognition using multi－scale auto convolution［J］.IEEE Transactions on Pattern analysis and Machine Intelligence，2005，27（6）：908－918.

［9］雷琳，蔡红苹，唐涛，等.基于MSA 特征的遥感图像多目标关联算法［J］.遥感学报，2008，12（4）：586－592.Lei Lin，Cai Hong－ping，Tang Tao，et al.A MSA feature－based multiple targets association algorithm in remote sensing images［J］.Journal of Remote Sensing，2008，12（4）：586－592.

［10］Bar－Shalom Y，Kirubarajan T，Gokberk C.Tracking with classification aided multiframe data assoeiation［J］.IEEE Transactions on Aerospace and Electronic Systems，2005，41（3）：868－878.

［11］Mandal A K，Pal S，De A K，et al.Novel approach to identify good tracer clouds from a sequence of satellite images［J］.IEEE Transactions on Geoscience and Remote Sensing，2005，43（4）：813－818.

［12］Krickpatric S，Gelett J C D，Vecchi M P.Optimization by simulated annealing［J］.Science，1983，220（4598）：671－680.

［13］Winkler G.Image Analysis，Random Fields and Dynamic Monte Carlo Methods［M］.Berlin：Springer－Verlag，1999.

［14］高尚.模拟退火算法中的退火策略研究［J］.航空计算技术，2002，32（4）：20－26.Gao Shang.Research on annealing strategy in simulated annealing algorithm［J］.Aeronautical Computer Technique，2002，32（4）：20－26.