基于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.