劉志偉，張?jiān)聢@，張曉燕，何 姍，劉穎婷

?

基于SSOR預(yù)條件技術(shù)的快速相位解纏算法

*劉志偉1,2，張?jiān)聢@1，張曉燕1,2，何 姍1，劉穎婷1

(1. 華東交通大學(xué)信息工程學(xué)院，江西，南昌 330013；2. 毫米波國家重點(diǎn)實(shí)驗(yàn)室，江蘇，南京 210096）

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

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

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

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

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

1 最小二乘相位解纏繞

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

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

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

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

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

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

Aφ = ρ (3)

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

2 基于預(yù)條件技術(shù)的相位解纏繞算法

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

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

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

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

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

因此，采用預(yù)條件技術(shù)后，對方程Aφ = ρ的求解轉(zhuǎn)化為求解方程

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

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

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

2) 構(gòu)造預(yù)條件矩陣M本身的計(jì)算復(fù)雜度和內(nèi)存消耗要??；

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

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

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

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

3 數(shù)值算例

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

圖1 真實(shí)的單個高斯山相位

圖2 單個高斯山纏繞相位

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

[2] 王磊. 干涉合成孔徑雷達(dá)信號處理的研究[D]. 北京: 中國科學(xué)院電子學(xué)研究所, 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] 芮平亮. 電磁散射分析中的快速迭代求解技術(shù)[D]. 南京：南京理工大學(xué), 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

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

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

張?jiān)聢@(1981-)，女，浙江金華人，講師，碩士，主要從事計(jì)算機(jī)應(yīng)用技術(shù)研究(E-mail: aney0360_cn@sina.com);

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

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

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