王朝霞 劉永信 張 杰 范陳清

①(內(nèi)蒙古大學(xué)計(jì)算機(jī)學(xué)院 呼和浩特 010021)

②(內(nèi)蒙古大學(xué)電子信息工程學(xué)院 呼和浩特 010021)

③(自然資源部第一海洋研究所 青島 266061)

1 引言

3維成像雷達(dá)高度計(jì)[1，2]綜合了傳統(tǒng)高度計(jì)高精度的測(cè)高能力和合成孔徑雷達(dá)高分辨率的成像能力，使海面高程的面測(cè)量成為可能，是遙感技術(shù)的最新發(fā)展前沿。2016年9月15日，搭載世界首例3維成像雷達(dá)高度計(jì)的“天宮二號(hào)”空間實(shí)驗(yàn)室成功發(fā)射，在國(guó)際上首次實(shí)驗(yàn)驗(yàn)證了采用小入射角和短基線干涉測(cè)量技術(shù)實(shí)現(xiàn)寬刈幅海平面高度厘米級(jí)精度測(cè)量的工作機(jī)理。因3維成像高度計(jì)測(cè)高采用了與干涉合成孔徑雷達(dá)(Interferometric Synthetic Aperture Radar， InSAR)基本相同的干涉測(cè)量技術(shù)[1]，故而在高程反演階段必須通過相位解纏這一環(huán)節(jié)將濾波后位于主值區(qū)間的纏繞相位轉(zhuǎn)換為絕對(duì)相位。相位解纏是干涉測(cè)量技術(shù)的關(guān)鍵環(huán)節(jié)，解纏的結(jié)果將直接影響最終的高程測(cè)量精度。在過去的數(shù)十年中，許多學(xué)者致力于相位解纏算法的研究和改進(jìn)，但如何在強(qiáng)噪聲干擾環(huán)境中進(jìn)一步提升干涉相位解纏的效率和速度，目前仍是一個(gè)富有挑戰(zhàn)的研究課題。

現(xiàn)有的微波干涉測(cè)量相位解纏方法主要有兩類，第1類是與展開路徑有關(guān)的局部法，典型的有枝切法、質(zhì)量圖導(dǎo)引法[3]、最小斷點(diǎn)法等；第2類是與展開路徑無(wú)關(guān)的全局法，典型的有最小二乘法[4]、最小費(fèi)用流法[5]、最小零范數(shù)法等。此外，也有綜合運(yùn)用以上兩類方法或提出了基于混合數(shù)學(xué)模型的方法[6]，以在解決實(shí)際問題中尋求更好的實(shí)驗(yàn)結(jié)果。其中，Goldstein等人[7]于1988年在研究In-SAR時(shí)提出的枝切法是一種非常經(jīng)典的算法。該方法通過識(shí)別殘差點(diǎn)，并根據(jù)殘差原理設(shè)置正確的枝切線，從而選擇合適的積分路徑，實(shí)現(xiàn)相位解纏。但殘差原理只是不連續(xù)相位存在的充分非必要條件，而且枝切線的放置因方法不同往往結(jié)果也不相同，不合理的枝切線放置可能導(dǎo)致相位跳躍和無(wú)法解纏的“孤島問題”[8]。針對(duì)枝切法存在的缺點(diǎn)和問題，近年來國(guó)內(nèi)外學(xué)者在理論研究和應(yīng)用中進(jìn)行了改進(jìn)。蔣銳等人[9]提出了等效殘差點(diǎn)的概念，并基于等效殘差點(diǎn)設(shè)置枝切線，解決了積分路徑穿過殘差點(diǎn)密集區(qū)所引起的展開相位跳變問題。張妍等人[10]通過虛擬組合殘差點(diǎn)改進(jìn)枝切線算法，在一定程度上解決了殘差點(diǎn)較密集區(qū)域在放置枝切線時(shí)易產(chǎn)生“孤島”的問題。De Souza等人[11]在獲取物體3維輪廓時(shí)，利用跳變相位解決了剩余殘差點(diǎn)的相位平衡問題，進(jìn)而通過尋求最小化處理時(shí)間實(shí)現(xiàn)了動(dòng)態(tài)全息技術(shù)中的干涉相位解纏。

本文在研究3維成像高度計(jì)高程反演時(shí)，提出了一種改進(jìn)的枝切線干涉相位解纏方法，利用JVC(Jonker-Volgenant-Castanon)全局最優(yōu)分配算法在正負(fù)殘差點(diǎn)之間放置枝切線，以進(jìn)一步縮短枝切線總長(zhǎng)度，避免“孤島”產(chǎn)生，提升干涉相位圖中像素的解纏率，降低高程反演誤差。

2 相位解纏與Goldstein枝切線算法

由于復(fù)數(shù)相位的周期性，對(duì)配準(zhǔn)后兩幅相干復(fù)圖像作復(fù)共軛相乘得到的干涉相位圖中，相位值是被周期折疊后的相位主值，位于(-π，π]之間，與真實(shí)相位值之間相差2 π的整數(shù)倍。

其中，Φ(i，j)為真實(shí)相位值，Ψ(i，j)為纏繞相位值，k(i，j)為整數(shù)。相位解纏的目標(biāo)就是恢復(fù)被模糊掉的周期分量2k(i，j)π。

一般情況下，干涉相位圖的采樣率均滿足Nyquist原理，因此干涉相位圖中相鄰像素點(diǎn)之差不超過半個(gè)周期，也就是差的絕對(duì)值不大于 π。在理想的無(wú)噪聲環(huán)境中，相位解纏只需根據(jù)干涉相位圖中像素值在距離向和方位向的偏導(dǎo)數(shù)進(jìn)行簡(jiǎn)單的積分即可。

然而在實(shí)際環(huán)境中，由于受噪聲等諸多因素的影響，干涉相位中存在不連續(xù)及不同方向上積分結(jié)果不一致的殘差點(diǎn)[12](相鄰4個(gè)像素點(diǎn)組成的正方形回路的相位差積分不為0，而是-2π(負(fù)殘差)或2π(正殘差))。根據(jù)相位解纏的基本原則，采用任何一種包含單個(gè)或多個(gè)正負(fù)殘差點(diǎn)數(shù)不等的積分路徑進(jìn)行解纏，都會(huì)產(chǎn)生不一致的結(jié)果。

為避免相位解纏后出現(xiàn)不一致的現(xiàn)象，Goldstein枝切線算法首先查找并確定干涉相位圖中所有的正負(fù)殘差點(diǎn)，然后在數(shù)量相等的相鄰正負(fù)殘差點(diǎn)之間放置“平衡”枝切線，如果正負(fù)殘差點(diǎn)數(shù)量不等，則把剩余的殘差點(diǎn)連接到圖像最近的邊界上，這樣也認(rèn)為是“平衡”的。進(jìn)而控制積分路徑，使其繞過這些枝切線，采用式(3)進(jìn)行相位解纏，即可得到一致解纏結(jié)果。因所有積分路徑包含的正負(fù)殘差點(diǎn)值總和都為0，所以不會(huì)形成全局誤差。但同時(shí)在枝切線另一側(cè)的像素點(diǎn)將會(huì)有相位不連續(xù)情況出現(xiàn)。為使不連續(xù)像素點(diǎn)數(shù)量最小，Goldstein算法在生成枝切線時(shí)要求滿足總長(zhǎng)度最短的原則。此外，在生成枝切線時(shí)，在殘差點(diǎn)相對(duì)密集部分有時(shí)形成多條枝切線環(huán)繞的像素區(qū)域(即“孤島現(xiàn)象”)，造成被環(huán)繞區(qū)域像素?zé)o法解纏的情況，這些都是改進(jìn)枝切線算法時(shí)需要重點(diǎn)考慮的問題。

3 基于JVC算法生成枝切線的解纏方法

JVC算法[13]是由Jonker， Volgenant和Castanon 3人聯(lián)合提出的一種全局最優(yōu)線性分配算法，常用于兩個(gè)數(shù)據(jù)點(diǎn)集之間的一對(duì)一分配。JVC算法將線性分配問題描述為：在約束條件

成立的前提下，求最小值

其中，c[i，j]為代價(jià)，x[i，j]為權(quán)值。

(2) 采用上述基于JVC算法生成枝切線的方法，將殘差點(diǎn)分為兩類，分別在殘差點(diǎn)與其最近邊界點(diǎn)之間或正負(fù)殘差點(diǎn)之間放置枝切線；最后遍歷所有的正負(fù)殘差點(diǎn)，若還有未連接枝切線的，則在其與最近邊界點(diǎn)之間放置枝切線，以此確保枝切線的極性平衡。

(3) 在枝切線相對(duì)較為稀疏的區(qū)域，選取一個(gè)非枝切線上的像素點(diǎn)作為起點(diǎn)，繞過所有枝切線采用式(3)的方法進(jìn)行積分，即可得到解纏相位。

(4) 對(duì)于枝切線上未能解纏的像素點(diǎn)，根據(jù)解纏后相位圖應(yīng)具有連續(xù)性的原理，將枝切線上點(diǎn)的像素值用7×7鄰域內(nèi)已解纏像素點(diǎn)的均值替換，完成全部像素點(diǎn)的相位解纏。

4 方法驗(yàn)證

4.1 3維成像高度計(jì)海面仿真干涉相位圖驗(yàn)證

為驗(yàn)證算法的可行性和有效性，本文根據(jù)“一發(fā)雙收”天線3維成像高度計(jì)的成像原理[1]仿真了Ku波段的海面干涉相位圖。首先采用P-M海浪譜和雙尺度模型模擬海洋表面，求得海面上每一點(diǎn)的高程數(shù)據(jù)。然后利用Delaunay三角剖分法將海面劃分為三角形小面元，采用物理光學(xué)(Physical Optics，PO)模型及其Kirchhoff近似解計(jì)算仿真海面區(qū)域的后向散射系數(shù)。最后根據(jù)3維成像雷達(dá)高度計(jì)的成像原理使用后向投影(Back Projection， BP)算法進(jìn)行仿真成像，得到兩幅仿真3維成像高度計(jì)海面相干復(fù)圖像。經(jīng)圖像配準(zhǔn)、去平地效應(yīng)及濾波后得到如圖1(a)所示的干涉相位。

在3維成像高度計(jì)海面仿真干涉相位圖中提取殘差點(diǎn)，共得到正殘差點(diǎn)578個(gè)，負(fù)殘差點(diǎn)577個(gè)。分別采用本文算法和Goldstein枝切線算法生成枝切線，結(jié)果如圖1(b)和圖1(c)所示。Goldstein枝切線算法生成的枝切線總長(zhǎng)度為3313像素，本文算法生成的枝切線總長(zhǎng)度為1723像素。

為驗(yàn)證本文算法的可行性和有效性，分別采用Goldstein枝切線算法、四向加權(quán)最小二乘法[14]、快速傅里葉變換(Fast Fourier Transform， FFT)算法[15]和本文算法對(duì)3維成像高度計(jì)海面仿真干涉相位圖進(jìn)行解纏，結(jié)果如圖2所示。Goldstein枝切線算法的像素解纏率為76.99%，本文算法的像素解纏率為88.03%，采用7×7鄰域內(nèi)的解纏像素均值替換枝切線上點(diǎn)的像素值后，解纏率可達(dá)100%。觀察解纏后的4幅干涉相位圖，可見四向加權(quán)最小二乘法解纏后相鄰相素之間的連續(xù)性在4種算法之中最好，階躍式跳變最少，本文算法與FFT算法次之，傳統(tǒng)Goldstein枝切線算法相對(duì)最差。

為定量比較4種算法的解纏結(jié)果，利用式(9)分別將解纏相位轉(zhuǎn)換為圖3(a)—圖3(d)所示的相對(duì)數(shù)字高程。

其中，h 為相對(duì)數(shù)字高程， H為3維成像高度計(jì)雷達(dá)平臺(tái)高度(本文取393 km)， Rm為雷達(dá)主天線到海面高程點(diǎn)的斜距， λ為電磁波波長(zhǎng)，φ 為解纏后的相位值， B為基線長(zhǎng)度(本文中取10 m)。

統(tǒng)計(jì)圖3(a)—圖3(d)所示數(shù)字高程與圖3(e)所示的原始模擬數(shù)字高程的均方根誤差(Root Mean Square Error， RMSE)，并在同樣的軟硬件環(huán)境(Intel Core i5-4590 4核，16 GB DDR3 1666 MHz內(nèi)存，Win7 x64旗艦版操作系統(tǒng))下分析4種算法的時(shí)間復(fù)雜度及實(shí)測(cè)運(yùn)行時(shí)間，詳細(xì)對(duì)比結(jié)果如表1所示。表中，分析時(shí)間復(fù)雜度時(shí)， M， N分別為干涉相位圖的行、列數(shù)，m ， n分別為正、負(fù)殘差點(diǎn)的個(gè)數(shù)，s=max(m，n)， W為四向加權(quán)最小二乘法控制加權(quán)窗口的尺寸， I為FFT算法的迭代次數(shù)。

圖1 3維成像高度計(jì)海面仿真干涉相位及其枝切線示意圖

圖2 3維成像高度計(jì)海面仿真干涉相位解纏結(jié)果

由圖3及表1可知，雖然四向加權(quán)最小二乘法解纏后相位值的連續(xù)性和算法運(yùn)行時(shí)間優(yōu)于其他3種算法，但解纏相位對(duì)應(yīng)的數(shù)字高程誤差大于本文算法和FFT算法，這是因?yàn)樗南蚣訖?quán)最小二乘法在進(jìn)行相位值擬合時(shí)導(dǎo)致了誤差的全局?jǐn)U散。采用4次迭代FFT算法時(shí)，雖然運(yùn)行時(shí)間優(yōu)于本文算法，但轉(zhuǎn)換后的高程誤差約為本文算法的2倍。本文算法的執(zhí)行時(shí)間優(yōu)于Goldstein枝切線算法、略次于FFT算法，但解纏精度均高于其他3種算法，能夠滿足3維成像高度計(jì)厘米級(jí)測(cè)高誤差的要求。

就時(shí)間復(fù)雜度而言，本文算法主要取決于殘差點(diǎn)的稠密程度，越稀疏時(shí)間復(fù)雜度越低。當(dāng)m ＜M， n ＜N ， 1+log2s ＜W2時(shí)，本文算法的時(shí)間復(fù)雜度將會(huì)低于其他3種算法；當(dāng)nlog2s ＞MN時(shí)，本文算法的時(shí)間復(fù)雜度將高于其他3種算法。介于上述兩種情況中間時(shí)，本文算法的時(shí)間復(fù)雜度一般低于Goldstein枝切線法，高于四向加權(quán)最小二乘法，接近FFT算法。

圖3 海面仿真干涉相位圖解纏后對(duì)應(yīng)的相對(duì)數(shù)字高程

表1 3維成像高度計(jì)海面仿真干涉相位圖4種算法解纏結(jié)果對(duì)比

4.2 實(shí)際Etna火山地區(qū)干涉相位圖驗(yàn)證

為說明算法的通用性和有效性，利用本文算法和Goldstein枝切線算法、四向加權(quán)最小二乘法、FFT算法分別對(duì)ERS1/2兩次經(jīng)過Etna火山地區(qū)時(shí)獲取的相干SLC(Single Look Complex)數(shù)據(jù)進(jìn)行解纏。本文從經(jīng)過多視處理的兩幅復(fù)圖像中分別選取500像素×500像素的區(qū)域進(jìn)行實(shí)驗(yàn)，經(jīng)配準(zhǔn)、去平地效應(yīng)、濾波后的干涉相位圖如圖4(a)所示。提取殘差點(diǎn)，得到正殘差點(diǎn)3383個(gè)、負(fù)殘差點(diǎn)3389個(gè)。本文算法和Goldstein枝切線算法生成的枝切線如圖4(b)、圖4(c)所示，本文算法和Goldstein枝切線算法生成的枝切線長(zhǎng)度分別為13546像素和15664像素。解纏后的結(jié)果如圖5所示。將解纏相位轉(zhuǎn)換為圖6(a)—圖6(d)所示的相對(duì)數(shù)字高程，并與圖6(e)所示對(duì)應(yīng)經(jīng)緯度的ASTER GDEMV2 30 m分辨率數(shù)字高程數(shù)據(jù)[16]進(jìn)行比較，計(jì)算RMSE，4種算法的詳細(xì)對(duì)比結(jié)果如表2所示(統(tǒng)計(jì)算法執(zhí)行時(shí)間的硬件環(huán)境與4.1節(jié)相同)。

圖4 Etna火山地區(qū)干涉相位及枝切線示意圖

圖5 Etna火山地區(qū)干涉相位圖解纏結(jié)果

由圖5、圖6及表2可知，雖然FFT算法和四向加權(quán)最小二乘法兩種全局性算法解纏后相位的連續(xù)性和執(zhí)行時(shí)間優(yōu)于本文算法，但解纏后相位所代表的高程走勢(shì)與原干涉相位圖不一致，且與GDEMV2之間的RMSE相對(duì)較大，這是由于兩種算法均引起了誤差的全局性擴(kuò)散。傳統(tǒng)Goldstein枝切線算法解纏雖未改變高程走勢(shì)，但未解纏相位和相鄰像素值之間的階躍式跳變較多，解纏結(jié)果連續(xù)性較差，且與GDEMV2的RMSE較大。本文算法兼顧了解纏的精確性和相位連續(xù)性，消除了大部分相鄰相素之間階躍式跳變，避免了誤差的全局?jǐn)U散。

分別對(duì)比圖2(b)、圖2(c)與圖4(b)、圖4(c)中兩種算法生成的枝切線，可以看出采用本文算法生成枝切線可以有效防止其相互連接形成“孤島”，使得干涉相位圖中除枝切線以外的像素點(diǎn)都可以通過積分進(jìn)行解纏。這是因?yàn)?，不論是與邊界點(diǎn)直接相連生成枝切線，還是通過JVC線性分配后在正負(fù)殘差點(diǎn)對(duì)之間放置枝切線，每個(gè)殘差點(diǎn)只存在于一條枝切線上，任何一條枝切線都是兩點(diǎn)之間的連線(忽略像素?cái)?shù)字化時(shí)取整的因素后可近似為直線段)，是無(wú)法連成閉環(huán)的。

圖6 Etna火山地區(qū)干涉相位解纏后對(duì)應(yīng)的相對(duì)數(shù)字高程

表2 Etna火山地區(qū)干涉相位圖4種算法解纏結(jié)果對(duì)比

5 結(jié)束語(yǔ)

本文為縮短枝切線相位解纏算法中生成枝切線的總長(zhǎng)度，提升干涉相位圖的解纏率和準(zhǔn)確性，提出了采用JVC全局最優(yōu)線性分配算法在殘差點(diǎn)之間放置枝切線，以此對(duì)Goldstein枝切線算法進(jìn)行改進(jìn)。通過對(duì)3維成像高度計(jì)海面仿真干涉相位圖和Etna火山地區(qū)干涉相位圖進(jìn)行解纏實(shí)驗(yàn)，結(jié)果表明在同樣的軟硬件運(yùn)行環(huán)境中，本文算法在生成枝切線總長(zhǎng)度、像素解纏率和運(yùn)行時(shí)間上均優(yōu)于Goldstein枝切線法，解纏精度優(yōu)于Goldstein枝切線法、四向加權(quán)最小二乘法和FFT算法，同時(shí)能夠較為有效地防止“孤島現(xiàn)象”產(chǎn)生，為干涉相位解纏提供了一種可行的方法。然而，在殘差點(diǎn)數(shù)量少且正負(fù)殘差點(diǎn)距離較近的干涉相位圖中，對(duì)Goldstein枝切線算法的改進(jìn)并不明顯；在信噪比低、殘差點(diǎn)密集的干涉相位圖(或局部區(qū)域)中，算法運(yùn)行時(shí)間相對(duì)其他算法較長(zhǎng)，且無(wú)法完全避免“孤島現(xiàn)象”的產(chǎn)生，這是下一步需要重點(diǎn)考慮的優(yōu)化方向。