加權最小二乘相位解纏的一種改進*

董春敏 劉國林 于勝文 劉偉科 周 偉

(1)山東科技大學測繪科學與工程學院，青島 266510 2)海島(礁)測繪技術國家測繪局重點實驗室，青島266510)

加權最小二乘相位解纏的一種改進*

董春敏1)劉國林1，2)于勝文1)劉偉科1)周 偉1)

(1)山東科技大學測繪科學與工程學院，青島 266510 2)海島(礁)測繪技術國家測繪局重點實驗室，青島266510)

基于最小二乘相位解纏原理，針對現(xiàn)有的加權和不加權最小二乘在殘差點過多時計算結果平滑問題，提出利用相位導數(shù)變化圖對其進行加權改進。結合真實數(shù)據(jù)，通過現(xiàn)有的加權與不加權最小二乘算法的比較和分析，證明該改進算法能有效處理殘差點過多的情況，解纏精度高，解纏結果更可靠。

InSAR;相位解纏;加權最小二乘;相位導數(shù)變化;殘差

1 前言

InSAR技術的應用已經涉及到很多領域，如地形測量、冰川研究、森林調查與制圖、海洋研究、獲取戰(zhàn)略、戰(zhàn)術機動和定期威脅軍事目標等。而相位解纏是利用InSAR技術獲取高精度目標高程信息的關鍵步驟，展開后的InSAR相位差直接關系到DEM信息提取的準確性和精確性。

近20年來，InSAR二維相位解纏算法研究發(fā)展十分迅速，它們大致可以分為3類。第一類是基于路徑跟蹤的算法［1］，主要有：Goldstein枝切法、質量圖指導算法、掩膜切割(Mask-cut)算法、Flynn最小不連續(xù)法等;第二類是基于最小范數(shù)的算法，主要有：不加權的最小二乘算法［2］、預解共軛梯度(PCG)算法、加權多網(wǎng)格算法、最小Lp范數(shù)算法［1］等;第三類是基于最優(yōu)估計的算法，主要有：基于網(wǎng)絡規(guī)劃的算法［3］、卡爾曼濾波法［4］、基于貝葉斯估計的方法［5］、遺傳算法［6］等。

Goldstein枝切法運算速度快，效率高，占用內存小，對噪聲小的區(qū)域能夠得到正確的解纏結果，但是對于含有大量不連續(xù)點的區(qū)域，由于不能設置正確的枝切線，解纏結果誤差大［7］;質量圖指導法對相位質量圖要求高，只有具備可靠的質量圖，解纏才能成功的進行;由于質量掩膜會造成一些孤立區(qū)域，Mask-cut算法適合于殘差點僅分布于低質量相位區(qū)的情況;Flynn最小不連續(xù)法和基于網(wǎng)絡規(guī)劃的算法適應性強，解纏結果較可靠，但是它們各自的算法復雜，所以耗費的時間較長;最小Lp范數(shù)算法解纏穩(wěn)定性較好，不需要識別殘差點，但是它們是“穿過”而不是“繞過”殘差點，解纏的同時會導致誤差向高質量區(qū)域傳遞。

本文在對最小二乘相位解纏算法研究基礎上，針對現(xiàn)有的最小二乘算法在殘差點過多時計算結果平滑問題，提出利用相位導數(shù)變化圖對其進行加權改進，利用MATLAB對新算法進行了實現(xiàn)，結合真實數(shù)據(jù)，借助原有加權最小二乘解纏算法，從目視和定量兩個方面對3種解纏算法進行了分析和比較。驗證了新算法的有效性。

2 最小二乘相位解纏

InSAR測量技術是利用從雷達復圖像數(shù)據(jù)中衍生出的相位信息提取地表三維信息。它的基本原理是通過兩部天線同時觀測，或者兩次平行的觀測，獲取地面上同一地物的兩幅復雷達圖像，利用兩幅復圖像相關運算得到的相位差反演地物的相對高程。但是實際上這樣得到的干涉相位差是真實相位差的卷疊，即InSAR干涉圖中與地面位置直接相關的相位是以2π為模的，所以為了計算每一點的高程必須給每一個相位測量值加上整數(shù)倍的相位周期，這種求解2π模糊性問題的技術稱為相位解纏。

對于纏繞相位函數(shù)Ψi，j(i=0，1，2，…，M-1;j =0，1，2，…，N-1)對應的解纏相位函數(shù)φi，j，則最小LP范數(shù)解必須滿足［1］

其中，

記

令δJ=0，將i，j互換，得到

對于邊界滿足ai，j=bi，j=0，那么有

假設擾動δφi，j是隨機的，則應有

求解式(10)，需要知道φ-1，j，φM，j，φi，-1，φi，N及對應的φ-1，j，φM，j，φi，-1，φi，N值，根據(jù)邊界條件有：

展開式(10)，并加上邊界條件，有：

當p=2，式(13)為離散形式的泊松方程

其中，

由分析可以看出，當p=2時，最小Lp范數(shù)的相位解纏問題等效于求解紐曼邊界的泊松方程，而此時的最小Lp范數(shù)算法也就是一般意義下的最小二乘相位解纏算法，它的目標是使纏繞相位Ψi，j與解纏相位φi，j的梯度差的平方和最小。最小二乘算法可以采用基于FFT或DCT的方法來加快算法的解纏速度。

基于FFT的最小二乘算法的具體實施步驟如下［8］：

1)計算0≤i≤M，0≤j≤N時的ρi，j值。

2)計算ρi，j的傅立葉變換。首先對每一行的ρi，j作鏡像對稱操作，得

然后進行傅立葉變換Fn，替換ρi，j(0≤n≤N)。行處理完后，對每一列進行同樣的操作，得到Pm，n。

3)計算Φm，n。

4)對Φm，n作反傅立葉變換，得到解纏函數(shù)φi，j的最小二乘估算值。

實踐證明，在殘差點過多時，基于FFT無權重最小二乘相位解纏算法的計算結果較平滑，但與實際數(shù)據(jù)存在很大的差異，因此有必要考慮引入權重，抑制誤差的傳播。

將最小二乘法得到的結果進行反纏繞，發(fā)現(xiàn)在原纏繞條紋密集的地方，經最小二乘法解纏后反纏繞時條紋明顯變稀。這是由于最小二乘法的平滑作用，使其在解纏過程中對真實相位的逼近過程中出現(xiàn)峰削尖、谷添底、陡坡變緩的趨勢。Bamler和Adam等人［10］指出，所有線性解纏方法都不可避免非零展開誤差。

3 改進加權最小二乘相位解纏

3.1 原有加權最小二乘

式中unit(.)為歸一化處理，filt［.］為低通濾波，求出權重wi，j后，用wi，j對ρi，j進行加權處理：

式中的K為控制加權強度，其值由實驗確定。權重wi，j取自二階差分的模，它對原相位中坡度不變的區(qū)域加權為0，僅對坡度有變化的區(qū)域有實際的加權，且坡度變化越快加權越重，因此能較好的抵消最小二乘的平滑作用。計算量上較最小二乘增加的不多［7］。

在干涉圖含噪聲多，相位不連續(xù)，殘差點過多的情況下，原有加權最小二乘法相位解纏的不足使其相比較不加權最小二乘解纏結果精度提高不多，甚至不如不加權最小二乘法，解纏效率低，耗時較長。

3.2 改進加權最小二乘法在采用加權最小二乘實現(xiàn)相位解纏時，權重的選取直接影響到算法的性能，甚至是決定算法成敗的關鍵。在無法獲得相干系數(shù)圖的情況下，相位導數(shù)變化圖是目前應用最廣，同時也是較可靠的一種權重選取方法。嚴格來講，相位導數(shù)變化描述了相位數(shù)據(jù)質量差的程度，而不是質量好的程度。改進的加權最小二乘相位解纏的權重依據(jù)相位導數(shù)變化圖來確定，即

定義梯度權重為：

加權相位拉普拉斯操作數(shù)為：

用Ci，j對ρi，j進行加權處理，即

4 實驗分析

試驗數(shù)據(jù)來源于濟寧某地區(qū)的InSAR圖像，通過GAMMA軟件對其進行形變提取，裁取感興趣區(qū)域(干涉條紋明顯)，得到了去平地和濾波處理后的干涉相位圖(160×146)如圖1(a)所示。由于干涉相位圖中仍然存在較大的噪聲，經過極大似然濾波處理后干涉相位圖如圖1(b)所示。可以看出，濾波后干涉條紋變得明顯清晰。殘差點分布圖如圖1 (c)所示(圖中黑點表示殘差)，3種解纏方法的解纏結果見圖23(圖像右邊的色度條表示相位值，單位為弧度)，解纏后重纏繞結果如圖3所示。

解纏方法目視評價：

1)從殘差點分布圖上可以看出，殘差點很多，主要分布在噪聲較多的區(qū)域。

2)3種解纏方法均能成功的進行解纏，解纏結果(圖2)相對比較平滑。由于各種方法對誤差的控制能力不同，所以中心區(qū)域的大小不同。

3)從圖3可以看出，由于最小二乘的平滑作用，重纏繞的條紋相比原纏繞條紋明顯變稀。改進的加權LS解纏方法相對來說能更好的恢復干涉條紋，其他兩種方法目視上區(qū)別不大。

圖1 去平地和濾波處理后的干涉相位Fig.1 Interfere phase after flat removing and filtering

圖2 3種解纏方法的解纏結果Fig.2 Results with three unwrarpping methods

圖3 解纏后重纏繞結果Fig.3 Rewrapping results after unwrapping

解纏方法定量評價：

根據(jù)不連續(xù)點數(shù)目、ε值［1］和解纏重纏繞結果與纏繞相位的差值3個指標評價改進的加權LS解纏方法的性能。不連續(xù)點數(shù)目越小，抗相位畸變的性能越好;ε值越小，則解纏質量越高;解纏重纏繞結果與纏繞相位的差值越小，可靠性越好。表1給出了3種方法的不連續(xù)點數(shù)目、ε值和解纏時間。表2給出了3種方法解纏重纏繞結果與纏繞相位的差值指標。

表1 3種解纏方法不連續(xù)點數(shù)目、ε值和解纏時間表Tab.1 Comparison among discontinuous point numbers，ε values and time of phase unwrapping with three unwrapping methods

表2 3種解纏方法解纏重纏繞結果與纏繞相位的差值表Tab.2 Differences between rewrapped results with three unwrapping methods and original wrap phase

從表1～2可以看出，改進的加權LS解纏方法不連續(xù)點數(shù)目和值、誤差絕對值的均值和中誤差均小于其他兩種解纏方法。解纏時間較不加權LS解纏方法增加的不多。說明本文改進方法在殘差點過多的情況下，抗相位畸變能力更好，解纏質量更高，可靠性更好。而原有加權最小二乘方法在本實驗中對于噪聲多，殘差點豐富的干涉圖解纏結果反而在解纏質量和效率上都不如不加權LS解纏方法，這也證明了傳統(tǒng)加權方法的不足。

4 結論

針對基于FFT變換的不加權LS解纏方法和原加權LS解纏方法在殘差點過多情況下解纏結果較平滑這一情況，提出了一種利用相位導數(shù)變化圖進行加權改進的算法。利用真實數(shù)據(jù)，在目視和定量兩個方面分析和評價了3種方法的性能。試驗證明，改進的算法在殘差點過多的情況下，具有更好的抗相位畸變能力以及更高的解纏質量和可靠性。本文提出的算法不僅對小像元圖像適用，對于較大型程序的編制和干涉圖的解纏也有一定的適用性和參考價值。

1 Ghiglia D C and Pritt M D.Two-dimensional phase unwrapping：Theory，algorithms and software［M］.John Wiley＆Sons，1998.

2 Ghiglia D C and Romero L A.Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods［J］.J Opt Soc Am A.，1994，11(1)：107-117.

3 Costantini M.A phase unwrapping method based on network programming［A］.Proc.Fringe’96 Workshop-ERS SAR interferometry［C］.Zurich：［s.n.］，1996.

4 劉國林，等.卡爾曼濾波在InSAR噪聲消除與相位解纏中的應用［J］.大地測量與地球動力學，2006，(2)：66-70.

5 Dias J M B and Leitao J M N.InSAR phase unwrapping：a Bayesian approach［J］.Geoscience and Remote Sensing Symposium，2001，IGARSS'01：396-400.

6 趙爭，張繼賢，張過.遺傳算法在InSAR相位解纏中的應用［J］.測繪科學，2002，27(3)：37-40.

7 許才軍，王華.InSAR相位解纏算法比較及誤差分析［J］.武漢大學學報(信息科學版)，2004，29(1)：67-71.

8 Mark D，Pritt and Jerome S Shipman.Least-squares two-dimensional phase unwrapping using FFT’s［J］.IEEE Trans，1994，32(3)：706-708.

9 李平湘，樣杰.雷達干涉測量原理與應用［M］.北京：測繪出版社，2006.

10 Bamler R，et al.Nosie-induced slope disotrion in 2-D phase unwrapping by linear estimators with application to SAR interferometry［J］.IEEE Trans，1998，36(3)：913-921.

AN IMPROVED ALGORITHM OF WEIGHTED LEAST-SQUARES FOR PHASE UNWRAPPING

Dong Chunmin1)，Liu Guolin1，2)，Yu Shengwen1)，Liu Weike1)and Zhou Wei1)

(1)Geomatics College，Shandong University of Science and Technology，Qingdao 266510 2)Key Laboratory of Surveying and Mapping Technology on Island and Reef，State Bureau of Surveying and Mapping，Qingdao266510)

Phase unwrapping is the key step in topographic elevation mapping and reconstruction of three-dimensional terrain model by Interferometric Synthetic Aperture Radar(InSAR).Considering the problem that the existing weighted and unweighted least-squares phase unwrapping results need smoothing，an improved algorithm using phase derivative variance map based on the principle of Least-squares phase unwrapping is proposed.With the real InSAR data，by comparison of the improved algorithm with the existing weighted and unweighted least-squares algorithm，it is verified that the proposed algorithm can effectively deal with the situation that there are too many residues and can get a higher accurate and reliable phase unwrapping results.

InSAR;phase unwrapping;weighted least-squares;phase derivative variance;residue

1671-5942(2011)Supp.-0090-05

2011-01-16

國家自然科學基金(40874001);海島(礁)測繪技術國家測繪局重點實驗室資助項目(2010A01)

董春敏，女，1985年生，山東濱州人，碩士，主要從事現(xiàn)代測量數(shù)據(jù)處理理論及應用研究.E-mail：chunmin0919@126.com

P227

A