改進的多基線相位解纏繞CANOPUS算法

劉會濤，張 歡，邢孟道，保 錚

（1.西安電子科技大學雷達信號處理國家重點實驗室，陜西西安710071；2.西安外事學院商學院，陜西西安710077）

改進的多基線相位解纏繞CANOPUS算法

劉會濤1，張 歡2，邢孟道1，保 錚1

（1.西安電子科技大學雷達信號處理國家重點實驗室，陜西西安710071；2.西安外事學院商學院，陜西西安710077）

中國余數定理（Chinese remainder theorem，CRT）方法是多基線相位解纏繞技術的一個基本方法，但是其較差的噪聲魯棒性問題限制了該方法在多基線相位解纏繞中的應用，然而基于聚類分析的多基線相位解纏繞技術能夠克服傳統的CRT算法噪聲魯棒性差的問題，本文根據CANOPUS算法中的聚類方法，提出用L∞-norm的距離測度定義兩點之間的距離，從而減少特小類的產生和降低噪點數，進而提高聚類分析的精度，并且改進CANOPUS算法的算法流程以提高聚類分析算法的執行效率，進而大幅度降低聚類分析的運算時間。通過用仿真數據和實測數據驗證可得，本文提出的改進聚類方法的聚類分析精度和算法執行效率更高，有效性通過實測數據實驗得到了驗證。

干涉合成孔徑雷達；多基線；聚類分析；相位解纏繞

0 引 言

干涉合成孔徑雷達（interferometric synthetic aperture radar，InSAR）通過利用兩幅SAR圖像的絕對相位差可以得到地表的高程信息。但是傳統的單基線InSAR技術存在無法有效解決復雜地形相位解纏繞問題和高程層疊效應問題［1－3］。然而通過增加不同頻率或者基線的雷達以不同的視角對成像場景進行觀測可以有效地克服單基線InSAR技術的技術瓶頸。多基線InSAR技術就是利用多個基線或者頻率雷達對同一場景進行觀測來實現復雜地形的三維測繪的。與單基線InSAR技術相比，多基線InSAR技術不僅可以有效地解決復雜地形的相位解纏繞問題和高程層疊效應問題，而且還可以提高模糊高度，進而有效地解決相位欠采樣問題［4］。然而多基線InSAR技術仍然存在相位解纏繞問題，而且與單基線相位解纏繞問題相比，多基線相位解纏繞問題還面臨著數據量大但可用信息較少的問題。

不同于單基線相位解纏繞方法，多基線相位解纏繞算法通過利用不同基線的干涉相位圖的多樣性實現相位解纏繞。理論上，只要基線長度滿足互質條件，纏繞相位可以根據中國余數定理（Chinese remainder theorem，CRT）的理論實現解纏繞［5－6］。但是由于CRT較差的噪聲魯棒性問題，通過傳統的CRT方法不能實現相位解纏繞。因此，要使傳統的CRT算法能夠應用于多基線相位解纏繞問題就必須增加CRT算法的噪聲魯棒性。文獻［7］通過深入研究CRT算法，提出了具有閉式解的魯棒的CRT算法。除此之外，還有一些算法通過研究多基線相位解纏繞所面臨的實際問題，將多基線相位解纏繞所獨有的特殊性質輔助傳統的CRT算法，以提高CRT算法的噪聲魯棒性。文獻［4］通過研究多基線數據所特有的聚類特性，將具有相同模糊數向量的像素聚為一個類，然后利用類內所有像素降低相位解纏繞所需數據的噪聲并將類內所有像素作為一個整體進行解纏繞聚類分析（cluster-analysis，CA）。再者，將CA算法和具有閉式解的魯棒的CRT算法相結合，從而在兩方面提高傳統CRT算法的噪聲魯棒性［8］。然而，由于CA算法中的聚類算法本身噪聲魯棒性不強，導致具有相同模糊數向量的像素未必能夠同時包含在一類內，為了改善CA算法的噪聲魯棒性，文獻［9］使用多基線數據的密度連通特性提高聚類的精確度，進而提高多基線相位解纏繞算法的噪聲魯棒性（cluster-analysis based noise robust phase-unwrapping algorithm，CANOPUS），同時該算法也是目前國際學術界提出的唯一一個能夠快速處理大規模多基線相位解纏繞的方法。除了基于傳統的CRT算法的多基線相位解纏繞算法之外，文獻［10］采用線性規劃的思想解決多基線相位解纏繞問題，提出了利用L1范數的多基線相位解纏繞算法。文獻［11］提出了用最大后驗估計（maximum a posterior estimation，MAP）的方法實現相位解纏繞，并且該算法在文獻［12］中得到了改進和發展。文獻［13－14］提出了利用最大似然估計（maximum likelihood，ML）的方法實現多基線相位解纏繞方法，文獻［15］利用CA算法中的聚類分析提出了一種改進的基于ML方法。文獻［16］發展和改進了MAP和ML估計方法，并給出了ML方法的克拉美羅界和MAP方法的誤差下界。文獻［17］提出用幅度信息輔助實現多基線相位解纏繞的方法。

本文首先回顧了CANOPUS算法，然后提出了改進的CANOPUS聚類算法，以使聚類分析運行速度更快，聚類結果更有利于多基線相位解纏繞，最后通過仿真數據和實測數據實驗驗證本文提出的改進的聚類方法的有效性。

1 CANOPUS算法

根據多基線InSAR中相位和基線長度的約束關系，可以得到模糊數、纏繞相位和基線長度之間的關系為（以雙基線InSAR系統為例）

式中，k、φ和B分別表示模糊數、纏繞相位和有效基線長度，其下標用于區別不同的干涉相位圖。與傳統的CRT算法相比，CA算法［4］利用具有相同模糊數向量的像素其對應的式（1）中的截距應該相同這一特性，將具有相同模糊數向量的像素聚為一類，并將類內所有像素視為一個整體對截距進行降噪處理，然后再進行相位解纏繞。然而，由于現實因素的影響，具有相同模糊數向量像素的截距不再一成不變，該截距值的變化將導致具有相同模糊數的像素不能被聚為一類和具有不同模糊數向量的像素可能被聚為一類兩方面的劣勢。為了克服截距空變的影響，CANOPUS算法［9］提出用截距數據（下文稱之為待聚類數據）的密度信息來區分不同的類。為了更好地理解本文提出的改進策略，首先詳細介紹CANOPUS算法。

與CA算法相比，CANOPUS算法將更多的待聚類數據中的信息用于聚類分析，不僅利用了截距信息，而且還利用了截距在整個待聚類數據中的相對位置信息。以雙基線數據為例，CANOPUS算法將三維的待聚類數據中的一點視為聚類模式，包含三方面信息，即點在待聚類數據中的行和列信息以及其截距值。如此選擇聚類模式的好處是既能通過截距值識別具有相同模糊數的類又能克服截距空變帶來的影響。假設三維待聚類數據中的一點p，定義以點p為球心ε為半徑的球體為點p的ε鄰域，并定義點p的ε鄰域內包含的待聚類數據中的點的數量為點p的密度Nε（p）。根據待聚類數據中截距值變化的特點，即相鄰類之間截距值變化較大（多基線數據截距值固有的屬性）而同一類內截距值緩慢變化（受現實因素影響導致的截距值的空變現象），則處于類邊界的點的密度理應比類內點的密度小，因此只要選擇合適的判斷準則就能區分出哪些點是邊界點（密度小于最小點數的點即為邊界點）。圖1為重復航過的多基線InSAR系統得到的待聚類數據及其密度圖。從圖1（a）可以看出，由于現實因素影響導致的截距的空變特性相當嚴重，除此之外，不同類之間的某些像素存在截距值的相互干擾，從而導致CA算法失效。然而，僅僅通過肉眼，不同的類還是很容易區分開來的，這是因為相鄰類之間的邊界點存在著明顯的截距值的差。圖1（b）為圖1（a）對應的密度圖，從圖1（b）可以看出，相鄰類之間的邊界點處其密度明顯低于類內其他像素。

圖1 實測數據的待聚類數據及其密度圖

從圖1可以看出，不同的類被低密度像素區分開來，因此只要選擇合適的最小點數即可得到類的邊界點。然而如何確定哪些像素屬于同一類呢？CANOPUS算法提供了一種基于密度連通的聚類方法，該方法是由模式識別中經典的聚類算法（density-based spatial clustering of applications with noise，DBSCAN）改進而來［11］。在闡述該方法之前需首先定義密度大于最小點數MinPts的待聚類數據中的點為核心點，相反，稱其為非核心點。如果某點p到核心點q的距離小于ε，則稱點p是從點q直接密度可達的。如果存在一系列的點p1，p2，…，pn，p1＝p，pn＝q，并且點pi是從點pi＋1（i＝1，2，…，n－1）直接密度可達的，則稱點p是從點q密度可達的。如果存在點o使得點p和pi＋1（i＝1，2，…，n－1）點q是從點o密度可達的，則稱點p和點q是密度連通的。CANOPUS所提供的聚類方法是將被邊界點分開的區域內所有的密度連通的點視為一類。其算法流程如下［9，18］：

將待聚類數據定義為Xun；

①設置分類編號m＝1；

②WHILE Xun≠Φ

③ 選擇Xun中任意一點p；

④ IF p是非核心點THEN

⑤ 將p賦予噪點集合N；

⑥ Xun＝Xun－｛p｝；

⑦ ELSE IF p是核心點THEN

⑧ 找到所有與點p是密度可達的點p

⑨ 并將其賦予類Cm；

⑩ Xun＝Xun－Cm；

從上述CANOPUS中的聚類算法流程可以看出，聚類分析開始于類內的任何一個核心點，終止于非核心點，而非核心點是由最小點數MinPts唯一決定的。因此只要選擇的參數合適，即可實現對待聚類數據的聚類分析，并且該聚類結果僅與選擇的參數值ε和MinPts有關。除此之外，該聚類算法的時間復雜度僅為O（n），其中n為待聚類數據的大小。

2 改進的CANOPUS聚類方法

由于CANOPUS算法中的CA方法并未考慮多基線InSAR數據的特點，因此本文從兩個方面改進CANOPUS算法。首先，CANOPUS算法中的聚類算法在計算任意兩點之間的距離時采用的是L2-norm距離測度，如此選擇的好處是其幾何意義更直觀，但是存在算法執行時間長和聚類精度不理想的問題。因此本文提出用L∞-norm距離測度定義兩點之間的距離。選擇L∞-norm距離測度的好處主要表現在兩方面：一是易于判斷任何兩點之間是否直接密度可達；二是對距離越遠的點其容許的截距變化相對而言比L2-norm距離測度更大，因此聚類過程產生的噪點的數量會降低。其次，在不改變算法實現功能的條件下，改變算法運行的流程以提高算法的執行速度。改進后的算法流程如下：

① 將待聚類數據定義為Xun；

② FORi＝1：n

③ IF Xun（i）是核心點

④ Xun（i）賦值核心點標志位；

⑤ END IF

⑥ END FOR

⑦ 設置分類編號m＝1；

⑧ WHILE Xun≠Φ

⑨ 選擇Xun中任意一核心點p；

⑩ 找到所有與點p是密度可達的點并將其賦予類Cm；

與傳統CANOPUS算法流程相比，本文中提出的改進算法首先識別并標識核心點與非核心點，然后再進行聚類分析。其優勢表現在兩個方面，首先，通過優先確認核心點后避免了原CANOPUS算法第4行執行非核心點的可能，從而提高算法執行效率。更為關鍵的是在傳統的CANOPUS算法中執行第9行確定與核心點p密度可達的點時，需要首先確認與核心點p直接密度可達的點并確認哪些點為核心點，然后再以所確認的核心點為源查找所有與這些核心點直接密度可達的點，以此類推，以至于再也沒有其他點與此類內的核心點直接密度可達。再在傳統的CANOPUS算法流程中，沒有提前確定并表示核心點的前提下，執行第9行需要重復計算核心點鄰域內的點是否為核心點，存在重復計算的可能，因此其算法運算效率較低。然而在本文中提出的改進的CANOPUS算法中，核心點的確認與標示被提前完成，因此在執行本文中算法第10行時，不存在重復計算某些點是否為核心點的運算，因此本文提出的對傳統CANOPUS算法的改進能夠極大地提高算法的執行效率。

3 算法的性能分析

本文提出的改進的CANOPUS算法將通過大小為458像素×157像素的Isolation Peak國家公園仿真數據和大小為7 697像素×8 000像素的TerraSAR-X實測數據驗證。Isolation Peak國家公園仿真數據的干涉相位圖如圖2（a）和圖2（b）所示，其基線長度分別為345.27m和281.46m。本文中仿真數據的相干系數約為0.8，因此在圖2（a）和圖2（b）中的干涉相位圖的噪聲較大。圖2（c）是待聚類數據，圖2（d）是用傳統的CANOPUS算法中的聚類算法得到的聚類結果，其包含74個特小類（類內聚類點數不超過30個），包含41 021個噪點數。圖2（e）是用本文提出的改進的聚類算法得到的聚類結果，其包含62個特小類，包含37 430個噪點數。與傳統CANOPUS算法中的聚類算法相比，本文提出的改進措施在減少特小類和減少噪點數兩方面都有改進，而且需要指出的是噪聲越嚴重，本文提出的改進方法對聚類算法性能的改善越大。圖2（f）是用本文提出的改進的聚類算法實現的相位解纏繞結果。

圖2 仿真數據的多基線相位解纏繞實驗結果

TerraSAR-X實測數據的主要參數如表1所示。本文實驗采用的實測數據為單顆衛星多次對同一地區進行重復觀測獲得的，即重復航過多基線數據。通過其主要參數可以發現，本次實驗所用數據存在較嚴重的時間去相干，因此數據受噪聲影響較嚴重。相位解纏繞結果如圖3所示，圖3（a）為短基線的干涉相位圖，圖3（b）為長基線的干涉相位圖，圖3（c）為利用長短基線干涉圖和基線數據根據式（2）得到的截距圖。通過用本文提出的改進聚類方法對圖3（c）進行聚類分析得到的聚類結果如圖3（d）所示，與其對應的多基線相位解纏繞結果如圖3（e）所示，從圖3（e）中可以看出，解纏繞結果中仍然存在大量的噪點，而導致噪點產生的主要原因是河流湖泊等相干性差的區域的存在。為了能夠清晰地對比本文中改進的聚類算法與CANOPUS算法中的聚類算法的聚類結果的差別，本文將傳統的CANOPUS算法中的聚類分析結果呈現在圖3（f）中。為了減小圖像大小，圖2中所有的圖像均經過16倍降采樣。

表1 實驗數據主要參數

本文從3個方面比較了CANOPUS中的聚類算法以及改進的聚類算法的聚類效果。首先從聚類分析結果中包含的聚類數來看，如圖3（d）和圖3（f）所示，本文中提出的改進算法聚類數降低了約69.1%，說明采用改進的聚類方法后能夠有效降低特小類的產生，而特小類的產生實際上是因為聚類算法聚類的過程中將原本屬于一類的點分成多個類，如圖3（d）和圖3（f）中白色矩形框內所示聚類結果。特小類的產生不會導致多基線解纏繞CRT算法失敗，但是會降低基于聚類分析的CRT算法的噪聲魯棒性（這是因為同一類內點數越少，根據類內點實現的濾波效果越差）。其次，從聚類結果中的噪點數來看，改進后的聚類算法產生的噪點數降低了約13.7%，說明不可解纏繞點減少了，解纏繞后的相位圖像更加精細。從聚類結果中的聚類數和噪點數兩方面來看，本文中提出的采用L∞-norm距離測度代替L2-norm距離測度的方法是行之有效的。最后從聚類分析的執行時間來看，本文提出的改進后的聚類方法其運算時間比原聚類算法降低了24.8%（CPU：Intel（R）Core（TM）i5 3.20GHz，RAM：6GB），說明本文中對聚類算法的流程的改進能夠極大促進算法的執行效率。表2給出了上述3點聚類算法性能分析的對比結果。

4 結束語

與單基線相位解纏繞技術發展的程度相比，多基線相位解纏繞技術仍然面臨較多難題。為了克服傳統CRT算法噪聲魯棒性差的問題，相位解纏繞的科技工作者們提出了用聚類分析的方法來增強傳統CRT算法的噪聲魯棒性，本文提出的聚類分析方法是對CANOPUS算法中的聚類分析方法的改進，以使其聚類分析結果更加精細，運算時間更短。通過實測數據驗證可以發現，本文提出的改進的聚類算法能夠減少特小類的產生，降低噪點數和降低聚類分析時的運算時間。

［1］Yu H W，Xing M D，Bao Z.A fast phase unwrapping method for large-scale interferograms［J］.IEEE Trans.on Geoscience and Remote Sensing，2013，51（7）：4240－4248.

［2］Zhang Y，Feng D Z，Qu X N.Hybrid phase unwrapping algorithm combining branch-cut and surface-fitting for InSAR［J］.Journal of Xidian University，2012，39（5）：47－53.（張妍，馮大政，曲小寧.支切法與曲面擬合結合的InSAR相位展開算法［J］.西安電子科技大學學報，2012，39（5）：47－53.）

［3］Yu H W，Li Z F，Bao Z.Residues cluster-based segmentation and outlier-detection method for large-scale phase unwrapping［J］.IEEE Trans.on Image Processing，2011，20（10）：2865－2875.

［4］Yu H W，Li Z F，Bao Z.A cluster-analysis-based efficient multibaseline phase-unwrapping algorithm［J］.IEEE Trans.on Geoscience and Remote Sensing，2011，49（1）：478－487.

［5］Xu W，Chang E C，Kwoh L K，et al.Phase-unwrapping of SAR interferogram with multi-frequency or multi-baseline［C］∥Proc.of the IEEE International Geoscience and Remote Sensing Symposium，1994：730－732.

［6］Jin G W，Zhang H M，Xu Q，et al.Phase unwrapping algorithm with CRT for multi-band InSAR［J］.Journal of Xidian University，2011，38（6）：97－102.（靳國旺，張紅敏，徐青，等.多波段InSAR的CRT相位解纏繞方法［J］.西安電子科技大學學報，2011，38（6）：97－102.）

［7］Wang W J，Xia X G.A closed-form robust chinese remainder theorem and its performance analysis［J］.IEEE Trans.on Signal Processing，2010，58（11）：5655－5665.

［8］Yuan Z H，Deng Y K，Li F，et al.Multichannel InSAR DEM reconstruction through improved closed-form robust chinese remainder theorem［J］.IEEE Geoscience and Remote Sensing Letters，2013，10（6）：1314－1318.

［9］Liu H T，Xing M D，Bao Z.A cluster-analysis based noise robust phase-unwrapping algorithm for multi-baseline interferograms［J］.IEEE Trans.on Geoscience and Remote Sensing，2015，53（1）：494－504.

［10］Yu H W，Bao Z.L1-norm method for multi-baseline InSAR phase unwrapping［J］.Journal of Xidian University，2013，40（4）：37－41.（于瀚雯，保錚.利用L1范數的多基線InSAR相位解纏繞技術［J］.西安電子科技大學學報，2013，40（4）：37－41.）

［11］Ferraiuolo G，Pascazio V，Schirinzi G.Maximum a posteriori estimation of height profiles in InSAR imaging［J］.IEEE Geoscience and Remote Sensing Letters，2004，1（2）：66－70.

［12］Ferraioli G，Shabou A，Tupin F，et al.Multichannel phase unwrapping with graph cuts［J］.IEEE Geoscience and Remote Sensing Letters，2009，6（3）：66－70.

［13］Pascazio V，Schirinzi G.Multifrequency InSAR height reconstruction through maximum likelihood estimation of local planes parameters［J］.IEEE Trans.on Image Processing，2002，11（12）：1478－1489.

［14］Fornaro G，Pauciullo A，Sansosti E.Phase difference-based multi-channel phase unwrapping［J］.IEEE Trans.on Image Processing，2005，14（7）：960－972.

［15］Yuan Z H，Deng Y K，Li F，et al.Improved multichannel InSAR height reconstruction method based on maximum likelihood estimation［J］.Journal of Electronics ＆Information Technology，2013，35（9）：2161－2167.（袁志輝，鄧云凱，李飛，等.改進的基于最大似然估計的多通道InSAR高程重建方法［J］.電子與信息學報，2013，35（9）：2161－2167.）

［16］Ferraiuolo G，Meglio F，Pascazio V，et al.DEM reconstruction accuracy in multichannel SAR interferometry［J］.IEEE Trans.on Geoscience and Remote Sensing，2009，47（1）：191－201.

［17］Shabou A，Baselice F，Ferraioli G.Urban digital elevation model reconstruction using very high resolution multichannel InSAR data［J］.IEEE Trans.on Geoscience and Remote Sensing，2012，50（11）：4748－4758.

［18］Theodoridis S，Koutroumbas K.Pattern recognition［M］.4th ed.Burlington：Academic Press，2009：764－862.

Improved CANOPUS method for multi-baseline interferograms

LIU Hui-tao1，ZHANG Huan2，XING Meng-dao1，BAO Zheng1
（1.National Lab of Radar Signal Processing，Xidian University，Xi’an 710071，China；2.College of Bussiness，Xi’an International University，Xi’an 710077，China）

As a basic technique against the multi-baseline phase unwrapping problem，Chinese remainder theorem（CRT）restricts itself in the application of multi-baseline phase unwrapping due to its worse noise robustness.However，multi-baseline phase unwrapping algorithms based on the cluster-analysis are able to overcome the drawbacks of the traditional CRT method.Based on the cluster-analysis method in the CANOPUS algorithm，an L∞-norm is employed to define the distance between two elements to decrease the production of very small clusters and the number of noise points so as to improve the cluster-analysis（CA）performance.Besides，the procedure of the CA method in the CANOPUS algorithm is changed to improve the efficiency of the method and decrease the consuming time of the algorithm.According to the experiments on a simulated and repeat-pass real interferometric synthetic aperture radar（InSAR）dataset，the effectiveness and efficiency of the improved cluster-analysis method are tested.

interferometric synthetic aperture radar（InSAR）；multi-baseline；cluster-analysis（CA）；phase unwrapping

TN 957.52

A

10.3969／j.issn.1001-506X.2015.08.08

劉會濤（1986－），男，博士研究生，主要研究方向為InSAR相位解纏繞處理、多基線相位解纏繞。

E-mail：huitaoliuxd＠gmail.com

張 歡（1984－），女，助教，碩士，主要研究方向為應用統計學。

E-mail：huanzhangxd＠126.com

邢孟道（1975－），男，教授，博士，主要研究方向為雷達成像和目標識別。

E-mail：xmd＠xidian.edu.cn

保 錚（1927－），男，院士，主要研究方向為雷達信號處理。

E-mail：piaofei8＠gmail.com

1001－506X201508－1767－06

網址：www.sys－ele.com

2014－10－27；

2015－01－26；網絡優先出版日期：2015－03－30。

網絡優先出版地址：http：／／www.cnki.net／kcms／detail／11.2422.TN.20150330.1111.015.html

國家自然科學基金優秀青年基金項目（61222108）資助課題