基于SSOR預條件技術的快速相位解纏算法

劉志偉，張月園，張曉燕，何 姍，劉穎婷

?

基于SSOR預條件技術的快速相位解纏算法

*劉志偉1,2，張月園1，張曉燕1,2，何 姍1，劉穎婷1

(1. 華東交通大學信息工程學院，江西，南昌 330013；2. 毫米波國家重點實驗室，江蘇，南京 210096）

相位解纏是合成孔徑雷達干涉測量中的一個關鍵步驟和研究熱點。在眾多的解纏算法中，最小二乘相位解纏算法以其優(yōu)良的穩(wěn)定性受到人們的關注。該方法的核心思想是將相位解纏問題轉化為通過迭代方法求解大型線性方程組。然而，傳統(tǒng)的迭代方法存在收斂緩慢，耗時過長的缺點。針對這一問題，本文提出了一種利用對稱超松弛預條件技術加速相位解纏的新方法。數值仿真實驗表明，與傳統(tǒng)方法相比，該方法可以在精確恢復真實相位的前提下，大大提高相位解纏的效率。

相位解纏；干涉合成孔徑雷達；廣義最小余量法；對稱超松弛預條件

干涉合成孔徑雷達(Interferometric Synthetic Aperture Radar, InSAR)是利用合成孔徑雷達復圖像的相位數據來獲取地面目標三維空間信息和變化信息的一項技術[1]。InSAR通過兩副天線同時觀測(雙天線模式)或兩次近平行的觀測(重復軌道模式)獲取地面同一場景的信息，并利用合成孔徑雷達成像技術獲取相應的兩幅復圖像[2]。由于目標與兩天線位置的幾何關系不同，使得在兩幅復圖像之間存在相位差。于是，我們便可以從兩幅復圖像共軛相乘后獲得的干涉紋圖出發(fā)，反演出場景中的高程信息。

在利用InSAR測高的實現過程中，由于受到三角函數主值的影響，使得干涉圖中的相位不再是真實的相位差，而是被周期折疊后位于[-π, π]之間的相位主值[3]，它與真實的相位差之間存在著2kπ差別，為整數。為了能夠精確恢復場景的高程，必須由干涉圖中的相位主值恢復出真實相位差，這一過程便被稱為相位解纏[4]。相位解纏是InSAR處理中的關鍵步驟，其準確程度將直接決定數字高程圖(DEM)和地表形變探測的精度。

當前，相位解纏的方法主要分為兩大類: 路徑追蹤法和最小二乘法[5]。其中，基于最小二乘的相位解纏算法是一種具備良好穩(wěn)定性的全局最優(yōu)方法，它將相位解纏問題轉化為通過求解離散的泊松方程，尋找纏繞相位的離散偏導數與解纏相位的偏導數整體偏差最小的解。針對這類大型稀疏矩陣方程，通常采用迭代算法加以求解，如Gauss-Seidel，GMRES迭代法等[6]。然而，由于傳統(tǒng)的迭代方法存在不穩(wěn)定和收斂速度慢等缺點，使得相位解纏的求解時間往往過長，求解精度也難以保證。為了解決這一難題，本文將預條件思想[7-8]引入到基于最小二乘法的相位解纏過程中。通過理論分析及數值仿真實驗驗證，預條件方法可以在不增加額外計算量的基礎上，加速迭代算法的收斂，從而顯著提高相位解纏的效率，減少計算時間。

1 最小二乘相位解纏繞

基于理想情況下解纏相位梯度等于纏繞相位梯度的假設，相位解纏可以看成是一個優(yōu)化問題。最小二乘法是一種廣泛使用的優(yōu)化方法。它的主要思想是基于解纏相位和纏繞相位之差的平方和最小準則，其目標函數為：

其中，φ,j為解纏后相位的估值，且

代表纏繞操作，目的是對纏繞相位差分進行加或減2π, 從而確保Δ,j和Δ,j位于[-,]之間。

對每一個像素(,)，方程(1)可改為如下的離散泊松方程

(φ+1,j- 2φ,j+φ-1,j) + (φ,j+1- 2φ,j+φ,j-1) =ρ,j(2)

由上面的分析可以看出，最小二乘的相位解纏算法等效于求解紐曼邊界下的泊松方程，它可以表示為一個未知量為(=×) 的稀疏線性方程組：

Aφ = ρ (3)

其中，A是一個稀疏矩陣，φ和ρ分別代表φ,j及ρ,j的一維列向量。

2 基于預條件技術的相位解纏繞算法

最小二乘算法將相位解纏轉化為等價的泊松方程之后，由于生成的大型稀疏線性方程組階數通常很高，使得直接求解需要大量的計算資源和存儲空間。因此，通常需要采用迭代方法對(3)式進行求解。求解稀疏線性方程組的迭代法[6]通常采用超松弛(SOR)和Krylov子空間方法等。前者一般稱為靜態(tài)迭代方法，而后者稱為非靜態(tài)迭代方法，采用何種方法求解，視具體情況而定。

對于靜態(tài)迭代方法，以Gauss-Seidel迭代為例，首先將(3)式中的陣列φ,j初始化為0，然后采用(4)式的迭代形式進行計算，直至收斂。

φ,j= 0.25 × (φ+1,j+φ-1,j+φ,j+1+φ,j-1-ρ,j) (4)

實踐證明，靜態(tài)迭代方法盡管數學形式比較優(yōu)美，但對低頻誤差分量的消除能力卻很弱，從而導致了很低的收斂速度。對于求解性態(tài)較差的線性系統(tǒng)，往往采用Krylov子空間方法，其中最典型的代表是采用廣義最小余量法(GMRES)求解線性系統(tǒng)。該方法于1986年由Saad和Schultz提出[6]，只需要直接的矩陣矢量乘(Matrix Vector Product: MVP)的信息就可以順利進行迭代，獲得滿足收斂精度的解，同時還具有穩(wěn)定的收斂速度，因此可應用于一般的非Hermitian矩陣方程的求解。

假定待求解的矩陣方程為Aφ = ρ，對于任意給定的一個初始近似解φ0，其相應的余量為r0= ρ – Aφ0。GMRES的迭代過程通過構造Krylov子空間span{r0, Ar0, …, A-1r0}的一組正交規(guī)范基，并在此空間內獲得對精確解的逼近，代表Krylov子空間的維數。GMRES算法具有穩(wěn)定性好，操作簡單的特點，然而，對于條件數大，即性態(tài)很差的系數矩陣，GMRES迭代算法也會出現收斂十分緩慢的現象。這是因為從本質上來說，迭代算法的收斂特性是由待求解系數矩陣的性態(tài)，即系數矩陣的條件數決定的，條件數較小時矩陣方程收斂較快，反之亦然。由此，我們不難想到，如果在迭代開始之前，利用另一個非奇異矩陣M對系數矩陣A進行預處理，使得預處理后的矩陣M-1A的特征值分布在幾個相對集中的區(qū)域，那么便可使預條件后的矩陣方程相比原方程具有更加良好的性態(tài)，即其特征譜是聚集的，從而有利于加速迭代算法的求解。

因此，采用預條件技術后，對方程Aφ = ρ的求解轉化為求解方程

M-1Aφ = M-1ρ (5)

從這里可以看出，預處理后的線性系統(tǒng)是否比原系統(tǒng)更易于求解，取決于是否構造一個高效的預條件算子。預條件算子的構造一般需要遵循以下基本原則[6,8]

1) M應該是A或者A-1在某種意義下的一個良好的近似；

2) 構造預條件矩陣M本身的計算復雜度和內存消耗要??；

3) 預條件操作的施加，即對任意一個矢量x，M-1x或Mx的計算量要小。

本文采用的對稱超松弛預條件(SSOR)來加速線性系統(tǒng)的求解，它完全滿足三條基本原則，相應預條件矩陣定義為[7]：

M = (D +E)D-1(D +F) (6)

其中，D、E、F分別表示系數矩陣A的對角陣、嚴格下三角陣以及嚴格上三角陣，為松弛因子，通常取值范圍為[0, 2]。我們注意到，由于SSOR預條件算子的分解因子是直接從系數矩陣A中抽取的，因此構造代價極小，同時還可以有效避免預條件算子構造過程中的不穩(wěn)定現象[6]。

3 數值算例

為了驗證結合預條件技術的迭代方法在相位解纏中的效果，針對典型的高斯山相位，我們進行了相應的數值仿真實驗。圖1和圖4分別為單個和兩個相鄰的高斯山相位，圖2和圖5為纏繞后的單個和兩個相鄰的高斯山相位，即纏繞相位。由圖2和圖5可以看出，纏繞相位為真實相位的主值，限于[-,]范圍之內。

圖1 真實的單個高斯山相位

圖2 單個高斯山纏繞相位

圖3 單個高斯山解纏后的相位

圖4 兩個相鄰的高斯山相位

圖5 兩個相鄰高斯山的纏繞相位

圖6 兩個相鄰高斯山解纏后的相位

圖3和圖6則分別反映了采用SSOR預條件技術之后，GMRES迭代方法進行相位解纏后的輸出結果。通過將圖3與圖1，圖6與圖4分別進行比較后不難看出，SSOR-GMRES算法可以準確的由纏繞相位恢復出真實相位，精確完成相位解纏的操作。

圖7 單個高斯山兩種解纏算法的迭代步數比較

圖8 兩個相鄰高斯山兩種解纏算法的迭代步數比較

為了進一步說明SSOR-GMRES算法在處理相位解纏問題上的有效性，我們分別采用SSOR-GMRES算法和傳統(tǒng)的Gauss-Seidel迭代法對圖2和圖4所示的纏繞高斯山相位進行相位解纏，并在統(tǒng)一的收斂精度(10-3)的前提下，將不同算法的迭代求解步數進行了比較，如圖7、圖8和表1所示。由圖表可以看出，SSOR-GMRES算法所需的收斂步數遠遠小于傳統(tǒng)的Gauss-Seidel迭代法，從而證明了預條件技術可用來加速干涉SAR中的相位解纏，具有較大的實用價值和現實意義。

表1 兩種解纏算法的迭代步數比較(迭代精度10-3)

4 結論

相位解纏是InSAR數據處理流程中的一個關鍵步驟，同時也是一個難點，因此，對解纏算法的研究具有重要意義。本文提出了一種將預條件技術與迭代法相結合，快速求解相位解纏問題的新方法，并就該方法的數學原理及實現手段進行了詳細闡述。理論分析及數值算例結果表明，預條件技術是一種加速迭代收斂的有效方法, 在求解相位解纏問題時具有收斂速度快、計算效率高、解纏精度高等優(yōu)點，具有較大的實用價值和廣泛的應用前景。

[1] 保錚. 雷達成像技術[M]. 北京:電子工業(yè)出版社,2005.

[2] 王磊. 干涉合成孔徑雷達信號處理的研究[D]. 北京: 中國科學院電子學研究所, 2001.

[3] Loffeld O, Nies H, Knedlik S, et al. Phase Unwrapping for SAR Interferometry - A data fusion approach by Kalman Filtering[J]. IEEE Trans. Geosicence and Remote Sensing, 2008, 46(1): 47-58.

[4] Goldstein R M, Zebker H A , Werner C. Satellite radar interferometry: Two-dimension phase unwrapping[J]. Radio Science, 1998, 23(4): 713-720.

[5] Bao M, Kwoh L K, Singh K. A improved least-squares method for INSAR phase unwrapping[C]. IEEE International symposium on Geoscience and remote sensing, 1998: 62-64.

[6] Yousef Saad. Iterative methods for sparse linear systems[M]. PWS Publishing Company, 1996

[7] Chen R S, Edward K N, Yung C H, et al. Application of SSOR preconditioned conjugate gradient algorithm to edge-FEM for 3-dimensional full wave Electromagnetic boundary value problems[J]. IEEE Trans. Microwave Theory and Techniques, 2002, 50(4):1165 -1172.

[8] 芮平亮. 電磁散射分析中的快速迭代求解技術[D]. 南京：南京理工大學, 2007.

PHASE UNWRAPPING BY MEANS OF PRECONDITIONING TECHNIQUE FOR INTERFEROMETRIC SAR BASED ON SSOR PRECONDITIONER Preconditioning technique for Interferometric SAR

*LIU Zhi-wei, ZHANG Yue-yuan, ZHANG Xiao-yan, HE Shan, LIU Ying-ting

(1. School of Information Engineering, East China Jiaotong University, Nanchang, Jiangxi 330013, China;2. State Key Laboratory of Millimeter Wave, Southeast University, Nanjing, Jiangsu 210096, China)

To solve the two-dimensional phase unwrapping problem, least square phase unwrapping is a robust approach, which is equivalent to the solution of a large sparse linear equation based on the iterative method. However, the convergence rate of the iterative method is usually very slow. To efficiently solve these linear systems arising from phase unwrapping problems, we propose the preconditioning technique. Numerical results demonstrate that the present method can reduce the number of iterations significantly, when no additional computing time is required to construct the preconditioning matrix.

phase unwrapping; interferometric synthetic aperture radar (InSAR); general minimal residual method (GMRES); SSOR preconditioner

1674-8085(2014)01-0046-05

O175.2

A

10.3969/j.issn.1674-8085.2014.01.010

2013-09-10；

2013-10-12

國家自然科學基金項目(61261005)；毫米波國家重點實驗室開放課題(K201326)；江西省科技廳青年基金項目 (20122BAB211018)；華東交通大學校立科研課題 (11XX01)

*劉志偉(1982-)，男，江西南昌人，講師，博士，主要從事計算電磁學研究(E-mail: zwliu1982@hotmail.com);

張月園(1981-)，女，浙江金華人，講師，碩士，主要從事計算機應用技術研究(E-mail: aney0360_cn@sina.com);

張曉燕(1979-)，女，云南楚雄人，副教授，博士，主要從事計算電磁學研究(E-mail: xy_zhang3129@sina.com);

何 姍(1989-)，女，安徽安慶人，碩士生，主要從事計算電磁學研究(E-mail: lucky-shan@163.com);

劉穎婷(1988-)，女，湖南冷水江人，碩士生，主要從事計算電磁學研究(E-mail: yingtingmm@163.com).