Goldstein枝切法對存在間斷相位缺陷的解纏研究

曾凡光，吳光敏，MAI John D，陳劍鳴*

(1.昆明理工大學(xué)理學(xué)院，昆明650500;2.香港城市大學(xué) 機械與生物醫(yī)學(xué)工程學(xué)系，香港)

引 言

作為干涉合成孔徑雷達(dá)、光學(xué)干涉、核磁共振圖像及新興的光學(xué)非接觸3維測量重建的關(guān)鍵步驟，2維相位解纏得到了迅速的發(fā)展和廣泛的應(yīng)用［1-2］。如果僅僅是將1維相位解纏的簡單推廣，2維相位解纏往往得不到滿意的結(jié)果，因為相位解纏的結(jié)果往往與1維分割的取向、信噪比和2維采樣率有著密切的關(guān)系［3］。在很多情況下，解纏誤差會沿著水平或垂直的方向傳遞。因此，相位解纏的難點就在于如何能很好地消除誤差對相位解纏的影響，來獲取精確的解纏相位［4］。

目前常用相位解纏的方法可以主要歸結(jié)為3類:(1)基于路徑積分的相位解纏［5］;(2)基于最小范數(shù)的相位解纏［6］;(3)基于網(wǎng)絡(luò)規(guī)劃的相位解纏算法［7］。噪聲或間斷相位在纏繞相位中較為常見，不同的解纏算法對于噪聲或間斷相位引起的殘差處理也不盡相同，主要包括:(1)找出優(yōu)化的積分路徑，避開殘差點造成的誤差傳遞［1］;(2)基于纏繞相位計算解纏相位的相位梯度估算值，將殘差點造成的誤差進(jìn)行平差處理［6］;(3)對殘差點不進(jìn)行任何處理，在解纏結(jié)果中予以保留［7］;Goldstein枝切法采用第1種方式處理，找出優(yōu)化的積分路徑［4］，通過建立合理的枝切線，隔絕噪聲區(qū)，阻止相位殘差點所造成的積分誤差在整幅圖像中的傳遞。

Goldstein枝切法是在殘差點電荷平衡的條件下用枝切線將附近的殘差點連接起來，或者把極性相反的殘差對連接在一起，或者將多個殘差點對組成的集合用枝切線連接起來［4］。通常殘差點與圖像邊界的連線也可以使殘差點達(dá)到極性平衡。Goldstein枝切線的連線策略是盡量使得枝切線最短［3］。Goldstein枝切法作為典型的路徑積分法，受間斷相位影響很大［3］。與此同時，殘差點的分布與枝切線的建立也有著緊密的內(nèi)在聯(lián)系。目前來看這類文獻(xiàn)較少見到。本文中對這一問題展開專門的探討和研究，通過模擬間斷相位代替真實數(shù)據(jù)的相位跳變或缺失等問題，通過3組實驗來研究不同的相位間斷形式(包括單間斷、雙間斷以及交叉雙間斷相位3種情況)對于解纏的影響。此外，還研究了枝切線搜索窗口半徑的對解纏的影響，驗證了存在“有效枝切線搜索窗口半徑”的結(jié)論。

1 Goldstein枝切法相位解纏原理

從數(shù)學(xué)角度而言，解纏的目的就是要正確地將相位重建起來。1維相位解纏以一個周期內(nèi)數(shù)據(jù)的連續(xù)性為前提條件的，沿著某一個方向上計算而進(jìn)行的。通常干涉雷達(dá)圖像數(shù)據(jù)是2維數(shù)據(jù)，通過兩幅復(fù)數(shù)圖像的配準(zhǔn)計算出每一像素點上的干涉相位差。2維相位是對1維相位的擴展［8］。

給定一個2維相位矩陣φi，j，為了對該相位矩陣進(jìn)行解纏，需要對每個點(i，j)加上或減去2π的整數(shù)倍值，從而得到一個連續(xù)的函數(shù) φi，j，k為整數(shù)矩陣:

圖1a中給出了1維纏繞相位，其中橫坐標(biāo)代表它的像素序列長度，縱坐標(biāo)代表的是纏繞相位值，圖1b中是圖1a的解纏結(jié)果值。

圖2a中給出了2維纏繞相位，其中橫和縱坐標(biāo)代表它的像素長度，豎直z坐標(biāo)代表的是纏繞相位值，圖2b中是圖2a的解纏結(jié)果值。

Fig.1 a—1-D wrapped phase b—1-D unwrapped phase

Fig.2 a—2-D wrapped phase b—2-D unwrapped phase

1.1 殘差點

實際的干涉相位數(shù)據(jù)由于噪聲和不連續(xù)點的影響，會出現(xiàn)不同路徑積分的結(jié)果［9-10］。為了描述有旋分量對相位解纏的影響，引入不連續(xù)點或殘差點的概念。在一幅原始纏繞數(shù)據(jù)中，定義一個2×2的窗口，以一個固定的方向，可以是順時針或者逆時針方向，利用下面的式子計算相鄰像素點之間的差值。其中q表示的4條邊的像素差值和，Δ表示2×2窗口內(nèi)的每條邊的像素差值，W表示纏繞操作，φW表示原始纏繞相位矩陣，i和j分別為行和列值。當(dāng)某一邊的差值超過范圍(－π，π)時，通過加減一個2π來處理，使得每條邊的差值范圍都保留在2π內(nèi)，最后得到的值再除以2π。

若q=0，則定義窗口內(nèi)最左上角的點(i，j)上的路徑為一致性路徑，否則就是不一致性路徑，且這個點(i，j)稱為殘差點。若q＞0，表示正殘差點;若q＜0，表示負(fù)殘差點。移動2×2窗口，直到找出整幅相位圖像的殘差點。

1.2 分支截斷理論

分支截斷理論的典型代表就是Goldstein枝切法，圖3是Goldstein枝切法的流程圖。通過計算相位圖像的殘差點，根據(jù)殘差點標(biāo)注的極性建立枝切線隔離區(qū)，最后通過建立的枝切線，在已知一個解纏相位點的情況下，標(biāo)注出它的4個領(lǐng)域像素點［11］。通過建立的枝切線尋找最優(yōu)積分路徑，避開殘差點所引起的誤差傳遞。如果是在枝切線外，則利用已知的解纏相位點進(jìn)行解纏繞，枝切線上的點留到最后進(jìn)行解纏，這種積分方法也俗稱洪水淹沒法［12］。

Fig.3 Flow chart of the Goldstein branch-cut algorithm

Goldstein枝切法的關(guān)鍵在于建立合理的枝切線，一條好的枝切線可以有效地阻斷相位殘差點所造成的解纏誤差傳遞［4-5］。根據(jù)殘差點的分布通常可以觀察判斷枝切線的放置是否有錯誤或存在被孤立的區(qū)域，而建立合理的枝切線又與枝切線內(nèi)的搜索半徑窗口有密切的關(guān)系。在枝切線搜索半徑窗口設(shè)置過小，它不能窮盡搜索到所有的殘差點分布;枝切線搜索半徑窗口設(shè)置過大，又有可能導(dǎo)致連接到圖像的邊界，造成封閉的“孤島”區(qū)域。

在如圖4a所示的6個殘差點中，黑色圓點表示負(fù)極性殘差點，白色圓圈表示正極性殘差點。圖4b所示枝切線搜索半徑窗口設(shè)置不合適，使得極性不平衡的殘差點連接到圖像的邊界處，導(dǎo)致錯誤的枝切線。當(dāng)建立圖4b中的枝切線時，在兩條枝切線半徑之間，其殘差點誤差并沒有完全隔離，當(dāng)解纏穿越兩條枝切線之間，會導(dǎo)致相位跳躍的情況。圖4c中給出了正確的枝切線連接，這條枝切連線由于涵蓋了圖像中所有殘差點，并且極性達(dá)到了平衡，有效地規(guī)避了殘差點對解纏的影響。

Fig.4 Branch-cut setting

2 間斷相位與枝切線搜索半徑的關(guān)系

為了對解纏結(jié)果進(jìn)行評價，采用人工模擬的高斯模型作為實驗數(shù)據(jù)，將解纏相位與真實人工相位的差值作為直觀的體現(xiàn)，并引入了反纏繞相位均方差進(jìn)行綜合評判。

（1）嚴(yán)格區(qū)分煤礦關(guān)閉和企業(yè)破產(chǎn)。煤礦關(guān)閉是行政行為，而企業(yè)破產(chǎn)是法律行為，二者政策大不相同。國內(nèi)煤礦關(guān)閉和企業(yè)破產(chǎn)之間沒有嚴(yán)格的區(qū)分，安置措施缺乏針對性。

反纏繞相位均方差是將相位進(jìn)行反纏繞［13］，把求得反纏繞相位和原始纏繞相位的均方差σ作為解纏精度評價的指標(biāo)之一，其中σ為:

式中，φW為原始纏繞相位矩陣，ψW為解纏相位經(jīng)過再次纏繞計算得到的相位矩陣;M，N為圖像的行列數(shù)。該指標(biāo)實質(zhì)是評價解纏前后的主相位值的保持性。

2.1 單間斷相位情況

在模擬高斯模型數(shù)據(jù)時，原始人工圖的大小為512pixel×512pixel。人工切割只有單間斷相位缺陷，并考慮在不同切割高度和切割深度的情況下，枝切線搜索窗口半徑的變化關(guān)系。

在切割高度為50pixel，切割深度為80pixel的情況時，計算它的有效枝切線搜索窗口半徑，如圖5所示。圖5a中給出了單間斷相位缺陷的2維顯示，圖5b中給出了單間斷相位缺陷的3維顯示，圖5c中給出了單纏繞相位缺陷的灰度2維顯示。

Fig.5 Gaussian distribution models for the single phase discontinuity

Fig.6 Goldstein branch-cut unwrapped result of the single phase discontinuity

通過計算發(fā)現(xiàn)，當(dāng)枝切線搜索窗口半徑處于65pixel～174pixel之間的任意一處值時，可以得到精確的解纏結(jié)果。圖6a是殘差點分布圖，圖6b是枝切線連線圖，圖6c是解纏相位的2維顯示圖，圖6d是解纏相位的3維灰度顯示圖，圖6e是反纏繞均方差曲線圖。其中橫坐標(biāo)表示的窗口半徑是枝切連線中尋找殘差點平衡時的窗口的半徑大小。從以上結(jié)果可以直觀看出，Goldstein枝切法能正確建立枝切線，獲得正確的解纏相位。從解纏相位與真實相位的差值圖也可以看出，在大部分區(qū)域中的解纏誤差為0。而在圖像的邊緣，由于只存在三領(lǐng)域，故在解纏時都被賦值為0，不參與解纏過程。在均方差與枝切線半徑的曲線圖(見圖6e)中，橫坐標(biāo)表示枝切線搜索窗口的半徑大小，縱坐標(biāo)表示反纏繞相位均方差的值。從圖6e中可以清楚地看出，在枝切線半徑范圍為65pixel～174pixel時，反纏繞相位均方差的值為0.1621。在枝切線搜索窗口半徑大于174pixel時，反纏繞相位均方差的值跳躍為0.5714。注意到在半徑小于65pixel時，反纏繞相位均方差的值仍然為0.1621，這主要與解纏結(jié)果有關(guān)。在枝切線半徑小于65pixel時，在間斷相位處，解纏結(jié)果有2π相位值的跳躍變化，而在計算反纏繞相位均方差值時，需要將解纏相位重新纏繞后再進(jìn)行計算，因此這里的2π相位值的跳躍變化不影響總體的反纏繞相位均方差值。

對于不同的切割高度和切割深度，有效枝切線搜索半徑大小范圍有所變化，在切割高度為40pixel，切割深度為40pixel時，枝切線搜索半徑窗口設(shè)置在73pixel～148pixel之間時，能夠得到正確的解纏相位。通過不同大小的切口計算得到的規(guī)律為:切口越大，能夠解纏的有效搜索半徑的下限值就越大，而上限值由枝切線是否與邊界接觸確定。當(dāng)枝切線連接到圖像的邊界，如果此時沒有搜索完所有的殘差點，遺留下未搜索的殘差點容易引起解纏誤差。如果已經(jīng)搜索完所有的殘差點，而此時連接到圖像邊界的枝切線如果沒有形成一個不可解纏的“孤島區(qū)域”，此時仍可以得到正確解纏。如果形成了一段封閉的“孤島”區(qū)域，此時解纏依然會有誤差。

2.2 雙間斷相位情況

首先討論有兩處不相交的切割間斷相位的情況。在兩處切割高度分別為60pixel和50pixel，切割深度都為80pixel時，如圖7所示。圖7a是兩處間斷相位缺陷的2維顯示，圖7b是兩處間斷相位缺陷的3維顯示，圖7c是纏繞兩處相位缺陷的灰度2維顯示。由于兩處切割的距離為160pixel，故這里的枝切線半徑窗口上限設(shè)置為160pixel。當(dāng)超過160pixel時，枝切線搜索窗口容易連接到另一處間斷相位的殘差點，形成局部不可解纏的“孤島”區(qū)域。由前面計算得知，在切割高度為50pixel時，它的有效枝切線搜索半徑為65pixel～174pixel。在切割高度為60pixel，切割深度為80pixel時，有效枝切線搜索半徑不小于89pixel。由于搜索窗口的變化步長一般為2pixel，為了方便確定搜索窗口內(nèi)中心殘差點坐標(biāo)，一般的搜索半徑窗口都是奇數(shù)，因此搜索半徑的最大值應(yīng)該為157pixel，由于158pixel的搜索窗口中心和157pixel是一樣的，故這里設(shè)置158pixel。

Fig.7 Gaussian distribution models for the two disjointed phase discontinuities

Fig.8 Goldstein branch-cut unwrapped result for two disjoint phase discontinuites

經(jīng)過實際計算可知在枝切線半徑處于89pixel～158pixel之間，對于兩處不相交的間斷相位能夠正確解纏，在枝切線搜索半徑處于這個范圍之外，解纏結(jié)果會有誤差。圖8a是殘差點分布圖，圖8b是枝切線連線圖，圖8c是解纏相位的2維顯示圖，圖8d是解纏相位的3維灰度顯示圖，其中橫縱坐標(biāo)表示像素點個數(shù)，豎直軸表示解纏相位值，圖8e是反纏繞均方差曲線圖，橫坐標(biāo)代表的是枝切線搜索窗口的半徑。從均方差與枝切線半徑的曲線圖(見圖8e)可以看出，在枝切線半徑大于158pixel時，它的反纏繞均方差顯著增大，而在有效枝切線半徑范圍內(nèi)，它的值固定在0.1621。

當(dāng)兩處不相交的間斷相位切割高度同為50pixel，切割深度同為80pixel時，它的有效枝切線搜索半徑范圍為65pixel～158pixel。這說明在有兩處不相交的間斷相位時，如果事先知道每處間斷相位能夠正確解纏的枝切線范圍，并且已知兩處間斷相位的像素距離，便可計算在兩處單獨有效的枝切線搜索半徑的交集而獲得的枝切線的有效范圍。

如果兩處間斷相位發(fā)生了相交，解纏將會變得困難。圖9中給出了兩處間斷相位垂直相交的情況，圖9b中的橫縱坐標(biāo)表示相位像素點坐標(biāo)，豎直z軸表示的是它的真實相位值。解纏的結(jié)果如圖10所示。

Fig.9 Gaussian distribution models of the two intersecting phase discontinuities

Fig.10 Goldstein branch-cut unwrapped result for two intersecting phase discontinuities

圖10a是殘差點分布圖，圖10b是枝切線連線圖，圖10c是解纏相位的2維顯示圖，圖10d是解纏相位的3維灰度顯示圖，其橫縱坐標(biāo)表示像素點坐標(biāo)值，豎直z軸表示解纏相位值。從解纏相位與真實相位的誤差圖(見圖10e)中可以看出，解纏結(jié)果中存在一個明顯的三角形的誤差區(qū)域。在均方差曲線圖(見圖10f)中，橫坐標(biāo)代表枝切線搜索殘差點窗口的半徑大小，縱坐標(biāo)代表計算得到的反纏繞均方差值。由于誤差引起的均方差存在著起伏跳躍式變化，無論它們的切割高度和切割深度是否相同，兩處間斷處都含有分布較為緊密的殘差點。尤其在間斷相位相交處，非常容易形成不可解纏的“孤島”區(qū)域，使得總體解纏結(jié)果有誤差。枝切線搜索半徑設(shè)置過大，容易連接間斷相位處的所有殘差點，造成由多條枝切線形成的封閉區(qū)域。如果枝切線搜索半徑過小，容易造成殘差點的未完全隔離，導(dǎo)致解纏誤差的全局傳遞。從所討論的情況來看，很難找到一個合適的枝切線搜索半徑能夠建立合理的枝切線，隔離所有殘差點的誤差，避免誤差傳遞。

3 結(jié)論

(1)在單相位間斷缺陷和雙不交叉的相位間斷缺陷的情況下，Goldstein枝切法仍然具有較好的解纏效果。

(2)對于雙交叉相位間斷缺陷，Goldstein枝切法在這一局部區(qū)域無法得到正確的解纏結(jié)果。

(3)不同的間斷相位缺陷對應(yīng)著不同的有效枝切線半徑范圍。在雙不相交切割間斷相位的情況時，可以根據(jù)有效枝切線半徑的交集得到正確的解纏結(jié)果;而在雙相交切割間斷相位的情況時，Goldstein枝切法無法得到有效的枝切線半徑范圍，也就無法得到正確的解纏相位。由此可以得出，Goldstein枝切法解纏存在著有效枝切線半徑。當(dāng)枝切線搜索半徑處于有效范圍之外，解纏結(jié)果一般都存在誤差。

本文中的實驗結(jié)果對于采用或聯(lián)合采用Goldstein枝切法進(jìn)行的相位解纏理論研究和應(yīng)用具有參考價值。

［1］ LI B，QIAN X F，LI X H，et al.Phase-unwrapping algorithm based on radial shearing principle［J］.Laser Technology，2013，37(1):44-47(in Chinese).

［2］ YANG F T，LUO J L，LIU Z Q，et al.Comparison of six phase unwrapping algorithms［J］.Laser Technology，2008，32(3):323-326(in Chinese).

［3］ GHIGLIA D C，PRITT M D.Two-dimensional phase unwrapping:theory，algorithms，and software［M］.New York，USA:Wiley，1998:103-121.

［4］ GOLDSTEIN R M，ZEBKER H A，WERNER C L.Satellite radar interferometry:two-dimensional phase unwrapping［J］.Radio Science，1988，23(4):713-720.

［5］ ITOH K.Analysis of the phase unwrapping algorithm［J］.Applied Optics，1982，21(14):2470-2470.

［6］ BONE D J.Fourier fringe analysis:the two-dimensional phase unwrapping problem［J］.Applied Optics，1991，30(25):3627-3632.

［7］ COSTANTINI M.A novel phase unwrapping method based on network programming［J］.IEEE Transactions on Geoscience and Remote Sensing，1998，36(3):813-821.

［8］ FLYNN T J.Two-dimensional phase unwrapping with minimum weighted discontinuity［J］.Journal of the Optical Society of A-merica，1997，14(10):2692-2701.

［9］ OESCH D W，SANCHEZ D J，TEWKSBURY-CHRISTLE C M.Aggregate behavior of branch points-persistent pairs［J］.Optics Express，2012，20(2):1046-1059.

［10］ HANSSEN R F.Radar interferometry:data interpretation and error analysis［M］.Dordrecht，Netherlands:Kluwer Academic Publishers，2001:22-23.

［11］ BAO Zh，XIN M D，WANG T.Radar imaging technology［M］.Beijing:Publishing House of Electronics Industry，2005:310-311(in Chinese).

［12］ ASUNDI A，WENSEN Z. Fast phase-unwrapping algorithm based on a gray-scale mask and flood fill［J］.Applied Optics，1998，37(23):5416-5420.

［13］ LEI S，ZHANG S.Digital sinusoidal fringe pattern generation:defocusing binary patterns vs focusing sinusoidal patterns［J］.Optics and Lasers in Engineering，2010，48(5):561-569.