一種基于多基線組合頻率估計的相位展開方法

謝先明，皮亦鳴

1.桂林電子科技大學信息與通信學院，廣西桂林541004；2.電子科技大學電子工程學院，四川成都611731

1 引 言

干涉相位展開是干涉合成孔徑雷達三維成像處理中的關鍵步驟，一直是干涉測量技術應用研究的熱點和難點［1－2］。利用多基線數(shù)據(jù)融合提高干涉相位精度是近些年來出現(xiàn)的一種新思想，它能克服單基線干涉SAR處理的缺陷，獲得更為精確的高程地圖（DEM），受到了越來越多的關注。在傳統(tǒng)的干涉SAR三維成像中，干涉基線長度與高度模糊數(shù)成反比，長基線對應的干涉圖條紋密集，能反映地形的精細結構，但增加了基線去相關，且條紋太密容易被噪聲污染，將增加相位展開誤差。另一方面，短基線對應的干涉相位圖條紋稀疏，有利于相位展開，但卻失去了地形的精細信息。多基線干涉SAR利用長短基線各自優(yōu)點，根據(jù)長短基線之間的關系更加準確地估計出長基線干涉相位，既有利于干涉相位展開，又能減小高度模糊數(shù)，保留地形的精細結構［3－7］。過去十幾年以來，有關多基線SAR干涉相位估計方法研究的文獻［8—14］已陸續(xù)發(fā)表。

目前，較為有效的多基線干涉相位估計方法主要有最小二乘（LS）和最大似然估計（ML）等方法［6，11］。其中LS算法能直接展開纏繞相位圖，但如果權值選擇不當或離散相位梯度估計不能反映真實的相位梯度，則會引入較大的誤差。盡管ML算法沒有LS算法中相鄰像素點的干涉相位差小于π的假設條件［6］，但它通常必須與其他方法結合才能獲得最終的展開相位。如文獻［6］中，ML算法只能獲得折疊的相位，還需用其他相位展開算法對折疊后的相位進行展開。另外，文獻［5］提出一種多基線迭代的相位展開方法，首先展開短基線干涉圖以獲得大致的DEM數(shù)據(jù)，再利用獲得的粗略DEM數(shù)據(jù)估計較長基線對應的干涉相位，以此類推，最終獲得最長基線對應的干涉相位估計值，迭代過程較為繁瑣。文獻［7］首先展開短基線對應的干涉圖，其次根據(jù)短基線和長基線的關系，估計出長基線對應的干涉相位被2π模糊的倍數(shù)，再根據(jù)長基線對應的干涉相位的主值和估計出來的模糊倍數(shù)恢復出長基線對應的干涉相位。但該方法會把前一級誤差擴大，如果擴大的誤差影響到后一級干涉相位模糊倍數(shù)估計時，則會導致誤差嚴重擴大而使得該相位展開后成為殘差點。文獻［8］把卡爾曼濾波（UKF）應用到多基線干涉相位估計之中，利用UKF強大的數(shù)據(jù)融合能力進行長基線干涉相位估計，具有較強的穩(wěn)健性和較高的估計精度。但該方法需對每一干涉圖進行濾波和展開處理，計算量較大，實時性不強。

本文提出一種多基線干涉相位估計方法，該方法選取適當?shù)幕€組合，用最大似然估計器［6］提取每一復像元隨基線變化的頻率，在展開短基線干涉相位基礎上，最終獲得長基線干涉相位估計值。該方法簡單有效，誤差傳遞較小，且不需要進行大量的迭代運算，甚至當短基線較短時，可不需要進行干涉相位展開。

2 多基線相位展開方法原理

設多基線干涉中垂直軌跡基線長度為Bi（i＝1，2，…，S），對應的基線傾角為ai（i＝1，2，…，S）。則歸一化的復干涉測量值可以表示為

式中，zi（k）為基線Bi所對應的干涉圖k像元復干涉測量值；φi（k）和φ′i（k）分別為基線Bi所對應的干涉圖k像元真實的以及含噪聲的干涉相位值；S指干涉基線的數(shù)目；?i（k）和δi（k）分別為雷達視線角以及幾何去相干等導致的干涉相位噪聲；λ為雷達波長。在多基線干涉中，通常有?i（k）＝?1（k），ai＝a1，因此，有

從而有

如果垂直軌跡基線長度Bi滿足則有是與i無關的；vk（i）＝exp［jδi（k）］可以看做為復高斯噪聲。于是，可以直接用最大似然估計器從序列y（k）＝［y1（k）y2（k）…yS－1（k）yS（k）］中獲得fk的估計值，于是φi（k）可以通過下式求得

只需要展開短基線B1對應的干涉圖φ1，就可以通過式（6）獲得最長基線BS所對應的干涉相位φS（k）。

而基線間隔ΔB必須滿足采樣定理，即

需要注意的是滿足式（7）要求的基線間隔非常小，因此很難滿足實用要求。實際的干涉圖由于平地效應的影響，干涉條紋密集不利于相位展開，通常需要先去除平地效應，再進行展開，且去平地效應之后的干涉圖能大致反映場景的相對高度。InSAR幾何關系如圖1所示，B為基線的長度，θ1為主天線入射角，雷達高度為H，α為基線傾角，R1和R2為主副天線到目標點P的斜距。

圖1 InSAR成像幾何關系圖Fig.1 InSAR geometry diagram

根據(jù)InSAR成像的幾何關系，垂直雷達視線方向基線長度為

去平地效應后的干涉相位為

因此，基線Bi（i＝1，2，…，S）去平地效應后的干涉相位可表示為

式中，h（k）和R1（k）分別是干涉圖k像元所對應的場景目標點高度以及到主天線的斜距，它們大致范圍通常是已知的同樣可以用最大似然估計器獲取fk的估計值。這時基線間隔ΔB應該滿足

3 基線間隔均勻時的頻率估計

以去平地效應后干涉圖為例來估計每一像素隨基線變化的頻率。去平地效應的復干涉圖k像元可以表示為

式中

對復干涉圖k像元信號作關于變量i的傅里葉變換，并作取模運算，即可得

式中，F(xiàn)［x］表示取x的傅里葉變換。圖2即為圖3所示的干涉圖對應的復干涉圖像在（23，25）像素的歸一化功率密度譜。本文使用最大似然頻率估計器來估計每一像素隨基線變化的頻率fk。具體實現(xiàn)時，最大似然頻率估計值通過搜索窗口內一維信號頻譜最大值處的頻率而得到。為了進一步提高估精度，可以通過給窗口內的信號末尾補零來增加頻率域的采樣數(shù)目。

圖2 復干涉像元功率密度譜Fig.2 Power spectrum of complex pixel

4 基線間隔非均勻時的頻率估計

前述方法實質上相當于對干涉基線B以滿足采樣定理的方式進行均勻抽樣，要求基線以相同的長度ΔB增加，條件比較嚴格，實現(xiàn)起來可能比較困難。因此，筆者也研究了基線間隔非均勻時的相位估計方法。

由前面的推導知道，去平地效應后復干涉圖k像元可以表示為

一般的，一維連續(xù)單頻信號x（t）的頻譜X（f）被定義為

由于不可能獲得時間無限長的信號測量值，因此，在實際的處理中通常把時間限定在一定范圍內，這里假設為信號時間寬度為T0。為了便于實際的計算，需要對連續(xù)的信號x（t）進行離散化處理，用離散的N0點采樣值代替連續(xù)信號x（t）。對于非均勻N0點采樣，用集合｛t1，t2，…，tN－1，tN｝代表采樣時刻，則連續(xù)信號x（t）的頻譜X（f）近似表示為［15］

式中，δ（t－tn）為脈沖函數(shù)，頻譜X（f）是連續(xù)的，需要進一步離散化為

通常采樣點數(shù)目是非常有限的，這導致頻率域采樣間隔較大，不利于高精度的估計其主頻。因此，為了進一步提高估精度，可以通過給窗口內的信號末尾補零來減小頻率域采樣間隔，假設補零后的信號采樣點數(shù)目為Nd，則此離散化頻率譜可進一步表示為

同樣的，通過索窗口內一維信號頻譜最大值的頻率，從而得到其頻率估計值。

5 仿真結果與分析

5.1 基線間隔均勻條件下仿真結果

為了驗證本文方法的有效性，對仿真多基線干涉SAR數(shù)據(jù)進行了處理。仿真基本參數(shù)如下：軌道高度為590km，下視角為45°，波長為0.136 6m，基線傾向角為10°，仿真場景為錐形場景（場景高度為390m），地面分辨率為8m×8m。本文仿真中，最短基線為200m，最長基線為550m，基線間隔為50m，相應的基線條數(shù)為8。

圖3（a）到圖3（h）為基線間隔為50m時的纏繞相位圖（圖3（e）、（f）省略），信噪比（SNR）為1.41dB。圖4（a）為干涉基線長度為550m時的真實干涉相位曲面。圖4（b）是枝切法［16］直接展開圖3（h）（即長基線干涉圖）的結果，其展開相位誤差及誤差統(tǒng)計直方圖見圖4（c）和圖4（d）。可以看出盡管枝切法展開相位誤差范圍較大，且主要分布在（0.2，0.5）區(qū)間，其平均展開相位誤差為0.35rad左右。

圖3 仿真干涉圖Fig.3 Synthetic interferogram

圖4 枝切法直接展開長基線干涉圖結果Fig.4 Branch－cut algorithm resolution

本文方法首先利用枝切法展開了圖3（a）所示的干涉圖（即最短基線干涉圖），其展開結果見圖5（a）。在展開最短基線干涉圖的基礎之上，利用最大似然估計器提取每一復像元隨基線變化的頻率，再利用公式（10）即可獲得圖3（h）所示的長基線干涉圖的展開結果，見圖5（b）。本文方法展開相位誤差及誤差統(tǒng)計直方圖見圖5（c）和圖5（d），可以看出本文方法展開相位誤差范圍很小，且誤差分布是以0為中心分布，其平均展開相位誤差相當小。

另外，也用LS算法進行多基線相位展開，其展開相位圖見圖6（a）。LS算法展開相位誤差及誤差統(tǒng)計直方圖見圖6（b）和圖6（c），可以看出多基線LS算法展開相位誤差范圍比本文方法展開誤差范圍大，且相位誤差主要分布在（0.1，0.55）區(qū)間，其平均展開相位誤差為0.31rad左右。為了定量分析相位展開算法精度，表1列出了在不同信噪比條件下多基線LS算法、枝切法（直接展開長基線干涉圖）以及本文方法的均方根誤差。另外，表1也列出了多基線ML算法展開誤差，需要注意的是ML算法通常通過2基線或者3基線聯(lián)合處理來估計較長基線的干涉相位，更多基線的聯(lián)合處理是非常困難的（且暫時尚未發(fā)現(xiàn)有文獻進行3基線以上的聯(lián)合處理）。表1列出的ML算法誤差是ML算法聯(lián)合處理B1基線干涉圖、B2基線干涉圖和B3干涉圖的結果。

圖5 本文方法展開長基線干涉圖結果Fig.5 Proposed algorithm resolution

圖6 LS方法展開結果Fig.6 LS algorithm resolution

表1 相位展開精度比較Tab.1 Comparison of phase unwrapping accuracy

5.2 基線間隔非均勻條件下仿真結果與分析

基線間隔均勻條件比較嚴格，現(xiàn)實起來可能比較困難。于是，對基線間隔非均勻時的相位估計也進行了研究，用本文方法實現(xiàn)了非均勻基線間隔條件下長基線干涉相位估計。仿真場景與主要參數(shù)與上一節(jié)相同，仿真基線長度分別為：B1＝200m、B2＝260m、B3＝290m、B4＝350m、B5＝400m、B6＝440m、B7＝490m、B8＝550m。圖7（a）為本文方法展開信噪比為1.41dB的干涉圖的結果，圖7（b）是展開誤差圖，圖7（c）為展開誤差直方圖。表2列出了在基線間隔非均勻與均勻條件下本文方法展開精度，可以看出非均勻間隔條件下本文方法展開精度要略微高一些。這是因為非均勻頻率估計具有更好的頻率檢測能力。

圖7 本文方法展開結果Fig.7 Proposed algorithm resolution

表2 基線間隔非均勻條件本文方法相位展開誤差Tab.2 Phase error of proposed algorithm under the condition of non－uniform baseline sampling

6 結 論

本文提出一種多基線干涉相位估計方法，并在多基線仿真數(shù)據(jù)處理中獲得了較好的效果。與最小二乘多基線展開方法等相比，本文方法具有簡單有效、計算量較小和誤差傳遞較小的特點。但由于該方法需要用最大似然估計器提取每一像元隨基線變化的頻率，基線間隔不能太大，否則提取出來的頻率是折疊的，這將引入嚴重的誤差。因此用本文方法來估計較長基線干涉相位時通常需要的干涉基線條數(shù)較多，代價較高，但這對于多航跡重復軌道高精度三維測高系統(tǒng)來說是可以接受的。另外，筆者也正在研究增大基線間隔（亦即基線間隔不必滿足采樣定理）的方法，這也是下一步工作的研究重點。

［1］ XIANG Maosheng，LI Shukai.InSAR Phase Unwrapping Using Poisson’s Equation with Dirichlet Boundary Conditions［J］.Acta Geodaetica et Cartographica Sinica，1998，27（3）：204－210.（向茂生，李樹楷.用狄氏邊界的泊松方程進行InSAR相位解纏的研究［J］.測繪學報，1998，27（3）：204－210.）

［2］ LIU Zilong，CAI Bin，DONG Zhen.A New InSAR Phase Unwrapping Method Based on Local Frequency Estimation and Multigrid Technique［J］.Acta Geodaetica et Cartographica Sinica，2010，39（1）：82－87.（劉子龍，蔡斌，董臻.基于局部頻率和多網(wǎng)格技術的InSAR相位解纏算法［J］.測繪學報，2010，39（1）：82－87.）

［3］ LOMBARDINI F，LOMBARDO P.Maximum Likelihood Array SAR Interferometry［C］∥IEEE Digital Signal Processing Workshop Proceedings.Loen：IEEE，1996：358－361.

［4］ LOMBARDO P，LOMBARDINI F.Multi－baseline SAR Interferometry for Terrain Slope Adaptivity［C］∥Proceedings of the 1997IEEE National Radar Conference.Syracuse：IEEE，1997：196－201.

［5］ THOMPSON D G，ROBERTSON A E，ARNOLD D V，et al.Multi－baseline Interferometric SAR for Iterative Height Estimation［C］∥IGARSS’99Proceedings：IEEE 1999International Geoscience and Remote Sensing Symposium：1.Hamburg：IEEE，1999：251－253.

［6］ ZHANG Qiuling，WANG Yanfei.Improving the Interferometric Phase Accuracy of Distributed Satellites InSAR System with Multibaseline Data Fusion［J］.Journal of Electronics and Information Technology，28（11）：2011－2014.（張秋玲，王巖飛.利用多基線數(shù)據(jù)融合提高分布式衛(wèi)星InSAR系統(tǒng)的干涉相位精度［J］.電子與信息學報，2006，28（11）：2011－2013.）

［7］ LIU Hui，ZHOU Yinqing，XU Huaping.An Estimate Method for Multi－baseline INSAR Phase［J］.Journal of Astronautics，2008，29（6）：1992－1994.（劉慧，周蔭清，徐華平.多基線干涉SAR的相位估計［J］.宇航學報，2008，29（6）：1992－1994.）

［8］ XIE Xianming，PI Yiming.Multi－baseline Phase Unwrapping Algorithm Based on the Unscented Kalman Filter［J］.IET Radar，Sonar ＆Navigation，2011，5（3）：296－304.

［9］ SHI X J，ZHAN Y H.A Multi－baseline Data Fusion Algorithm for Distributed Satellites SAR Interferometry by Combining Iterative and Maximum－likelihood Methods［C］∥Proceedings of 2009Asia－Pacific Conference on Synthetic Aperture Radar（APSAR’09）.Xi’an：［s.n.］，2009：165－168.

［10］ XU W，CHANG E C，KWOH L K，et al.Phase Unwrapping of SAR Interferogram with Multi－frequency or Multi－baseline［C］∥1994International Geoscience and Remote Sensing Symposium（IGARSS’94）.Piscataway：［s.n.］，1994：730－732.

［11］ GHIGLIA D C，WAHL D E.Interferometric Synthetic Aperture Radar Terrain Elevation Mapping from Multiple Observations［C］∥Proceedings of 6th IEEE Digital Signal Processing Workshop.Yosemite National Park：［s.n.］，1994：33－36.

［12］ KIM M G，GRIFFITHS H D.Phase Unwrapping of Multibaseline Interferometry Using Kalman Filtering［J］.Seventh International Conference on Image Processing and Its Applications.Manchester：IEEE，1999：813－817.

［13］ LI Z F，BAO Z.A Joint Image Coregistration，Phase Noise Suppression，and Phase Unwrapping Method Based on Subspace Projection for Multi－baseline InSAR Systems［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45（3）：584－591.

［14］ LACHAISE M，EINEDER M，F(xiàn)RITZ T.Multi Baseline SAR Acquisition Concepts and Phase Unwrapping Algorithms for the TanDEM－X Mission［C］∥2007International Geoscience and Remote Sensing Symposium（IGARSS’07），Barcelona：［s.n.］，2007：5272－5276.

［15］ GAO Yukai，ZHANG Wei.Spectrum Estimation from Non－uniformly Sampled Signals［J］.Electrical Measurement ＆Instrumentation，2009，46（2）：8－11.（高玉凱，張維.非均勻采樣信號的頻譜分析方法［J］.電測與儀表，2009，46（2）：8－11.）

［16］ GOLDSTEIN R M，ZERBER H A，WERNER C L.Satellite Radar Interferometry：Two－dimensional Phase Unwrapping［J］.Radio Science，1988，23（4）：713－720.