空間網(wǎng)絡(luò)對(duì)時(shí)序InSAR相位解纏的影響
——以Delaunay與Dijkstra網(wǎng)絡(luò)為例

馬張烽，蔣 彌，李桂華，黃 騰

1. 河海大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 211100； 2. 中山大學(xué)測(cè)繪科學(xué)與技術(shù)學(xué)院，廣東 珠海 519082

相位解纏是時(shí)序InSAR(interferometric synthetic aperture radar)技術(shù)至關(guān)重要的一個(gè)數(shù)據(jù)處理步驟[1-3]。能否準(zhǔn)確估算出纏繞于[-π,π)的相位觀測(cè)的2π模糊度個(gè)數(shù)是InSAR監(jiān)測(cè)成敗的關(guān)鍵。該步驟主要面臨兩個(gè)挑戰(zhàn)：①判斷出存在模糊度的觀測(cè)點(diǎn)；②如何正確解算出對(duì)應(yīng)的模糊度。相位連續(xù)性假設(shè)是解決上述挑戰(zhàn)的基礎(chǔ)，即相鄰兩個(gè)待解纏點(diǎn)的真實(shí)相位差(相位梯度)絕對(duì)值應(yīng)小于π[4]。

基于此假設(shè)，最小二乘、最小費(fèi)用流和最小不連續(xù)圖等Lp范數(shù)方法[5]通過(guò)最小化纏繞相位梯度與真實(shí)相位梯度的范數(shù)距離來(lái)獲得可靠的解纏結(jié)果。雖然上述方法在很多案例里都有出色的表現(xiàn)[6]，但低相干區(qū)域的強(qiáng)相位噪聲、大氣湍流和高相位梯度(變形或局部地形)都在不同程度地影響著局部相位差，導(dǎo)致相位連續(xù)性假設(shè)在實(shí)際的干涉圖中并不能處處成立，導(dǎo)致解纏誤差的產(chǎn)生并進(jìn)而影響InSAR監(jiān)測(cè)結(jié)果。

為減少相位不連續(xù)性帶來(lái)的解纏誤差，學(xué)者們采用了不同的處理策略。文獻(xiàn)[7]提出的一維+二維方法將時(shí)域最小費(fèi)用流網(wǎng)絡(luò)規(guī)劃與空間解纏相結(jié)合，以緩解較大相位梯度對(duì)結(jié)果造成的影響，從而逼近相位連續(xù)性假設(shè)。文獻(xiàn)[8]提出了一種兩步三維解纏算法，該算法通過(guò)在時(shí)間維進(jìn)行低通濾波緩解高相位梯度，并根據(jù)噪聲程度估計(jì)準(zhǔn)確的權(quán)重，然后在空間維使用網(wǎng)絡(luò)將所有離散點(diǎn)進(jìn)行連接并進(jìn)行相位模糊度的計(jì)算，獲得了較好的解纏精度。文獻(xiàn)[9]提出了基于三角無(wú)旋性假設(shè)的三維解纏算法，該算法利用時(shí)間和空間維度上的三角相位閉合差來(lái)優(yōu)化無(wú)旋性約束，逼近相位連續(xù)性假設(shè)，從而獲得了較穩(wěn)健的解纏結(jié)果。

雖然時(shí)間序列InSAR相位解纏歷經(jīng)近20年發(fā)展，但是解纏算法的理論框架卻大相徑庭。即便如此，大多數(shù)方法在構(gòu)建空間網(wǎng)絡(luò)時(shí)均使用了Delaunay三角網(wǎng)[10]，并在此網(wǎng)絡(luò)的基礎(chǔ)上根據(jù)網(wǎng)絡(luò)邊所包含的兩點(diǎn)之間的相位差解算所有點(diǎn)的模糊度。盡管目前Delaunay三角網(wǎng)應(yīng)用廣泛，且在許多軟件中(如GAMMA軟件)，均被作為默認(rèn)的解纏網(wǎng)絡(luò)[11]，但也不乏一些負(fù)面的案例。例如在橋梁等三維結(jié)構(gòu)較復(fù)雜的建筑物的解纏過(guò)程中，Delaunay三角網(wǎng)并未獲得可靠的解纏結(jié)果[12-13]。考慮到Delaunay三角網(wǎng)的表現(xiàn)以及在解纏中的定量表現(xiàn)仍不明確，本文從理論和實(shí)測(cè)角度評(píng)估了Delaunay三角網(wǎng)對(duì)時(shí)序相位解纏的影響。結(jié)果表明，由于Delaunay三角剖分算法在構(gòu)網(wǎng)時(shí)僅考慮了拓?fù)潢P(guān)系，并沒(méi)有將相位連續(xù)性假設(shè)作為約束，因此容易觸發(fā)高相位梯度的邊，導(dǎo)致解纏錯(cuò)誤。為解決這一問(wèn)題，本文首先使用時(shí)間相干性度量干涉圖中指定邊的相位梯度，再引入圖論中的Dijkstra最短路徑算法[14]，引導(dǎo)選擇低相位梯度的邊，從而逼近相位連續(xù)性假設(shè)。為公平比較網(wǎng)絡(luò)性能，本文在對(duì)比過(guò)程中均采用了相同的最小費(fèi)用流方法[15]對(duì)網(wǎng)絡(luò)進(jìn)行解纏。模擬和真實(shí)數(shù)據(jù)的對(duì)比研究表明，本文構(gòu)網(wǎng)方法較現(xiàn)有Delaunay網(wǎng)絡(luò)能夠更為高效地執(zhí)行相位解纏，提高時(shí)序解纏精度。

1 面向Delaunay三角網(wǎng)的解纏方法

Delaunay三角網(wǎng)自Dirichlet 1850年首次提出在空間約束下根據(jù)離散點(diǎn)集合生成三角網(wǎng)的思想后[10]，經(jīng)過(guò)170年的發(fā)展已經(jīng)產(chǎn)生了諸多變種，本文僅對(duì)經(jīng)典Delaunay生成算法進(jìn)行簡(jiǎn)要回顧。

在二維平面中給定K個(gè)離散點(diǎn)集合Pi,i=1,2,…,K，對(duì)每一個(gè)離散點(diǎn)Pi都生成一個(gè)相對(duì)應(yīng)的凸包γi，γi的特征表現(xiàn)為：γi上所有包含的點(diǎn)離Pj的距離總是要遠(yuǎn)于與Pi之間的距離(凸包所包含的點(diǎn)并非待生成三角網(wǎng)的離散點(diǎn))。這些凸包構(gòu)成的圖形被稱為Voronoi圖。Voronoi圖中一個(gè)公共點(diǎn)包含3個(gè)凸包γm、γj與γk，且3個(gè)凸包相鄰且互相存在公共邊，將3個(gè)凸包所代表的離散點(diǎn)Pm、Pj與Pk相連，則構(gòu)成了基本的Delaunay三角形。重復(fù)三角形的生成則構(gòu)成了最后的Delaunay三角網(wǎng)。

Delaunay三角網(wǎng)具有以下兩個(gè)特征：①構(gòu)網(wǎng)唯一，即所有離散點(diǎn)集合只能生成唯一的三角網(wǎng)；②三角網(wǎng)內(nèi)所有三角形互不重疊。

在最小費(fèi)用流解纏(MCF)方法中，生成的三角網(wǎng)中每個(gè)三角閉合環(huán)Trii都代表著最小的獨(dú)立閉合圖形。在每個(gè)三角圖形單元中，3條邊的纏繞相位差之間存在著無(wú)旋性約束條件Ui[16]

Ui=round(-δφm,j/2π)+round(-δφj,k/2π)+

round(-δφk,m/2π)

(1)

式中，δφm,j、δφj,k及δφk,m分別為三角形中鄰近3個(gè)點(diǎn)m、j和k的纏繞相位梯度；round為取整算子。

假設(shè)生成的三角網(wǎng)中含有M個(gè)三角形，且由N個(gè)邊組成，則所有邊E的待求模糊度與所有無(wú)旋性約束所對(duì)應(yīng)的設(shè)計(jì)矩陣AM×N表示如下[16-17]

(2)

聯(lián)合設(shè)計(jì)矩陣與無(wú)旋性約束可以求取所有邊的模糊度X，模糊度、設(shè)計(jì)矩陣與無(wú)旋性約束之間存在的線性關(guān)系及求解模糊度的目標(biāo)函數(shù)表示如下

(3)

式中，f代表相位梯度概率密度函數(shù)計(jì)算出的極大似然費(fèi)用[18]。

由于待求解模糊度的整數(shù)特征，且三角形個(gè)數(shù)M通常小于邊個(gè)數(shù)N，A存在秩虧現(xiàn)象，因此常規(guī)最小二乘無(wú)法求解?？紤]到部分邊不存在模糊度，即待求解K表示為稀疏特征，因此式(3)的解算需要由對(duì)于稀疏待求解有較好表現(xiàn)的整數(shù)線性規(guī)劃完成。由于采用的整數(shù)線性規(guī)劃求解式(3)時(shí)存在加權(quán)過(guò)程，這里的權(quán)值通常在網(wǎng)絡(luò)規(guī)劃中被認(rèn)為是通過(guò)路徑產(chǎn)生的費(fèi)用，因此加權(quán)的整數(shù)線性規(guī)劃也通常被稱為最小費(fèi)用流方法[19]，本文最小費(fèi)用流的求解過(guò)程均采用GUROBI軟件中的“intlinprog”整數(shù)線性規(guī)劃子函數(shù)完成。

需要注意的是，在上述目標(biāo)函數(shù)完成優(yōu)化后，僅計(jì)算出了每條邊的模糊度。因此還需要通過(guò)漫水填充(flood-fill)的方式，將邊的模糊度積分至每個(gè)點(diǎn)，最后將每個(gè)點(diǎn)的模糊度添加到初始纏繞相位上完成解纏。

2 面向Dijkstra網(wǎng)絡(luò)的時(shí)序解纏方法

2.1 Dijkstra最短路徑算法

最短路徑問(wèn)題是經(jīng)典圖論研究中的一個(gè)基本問(wèn)題，主要目標(biāo)是最小化指定起點(diǎn)i與終點(diǎn)j間的連通路徑的長(zhǎng)度或權(quán)重。Dijkstra方法是最短路徑搜索算法中最具代表性、計(jì)算效率最高的經(jīng)典算法之一，其核心思想在于[20]：給定一個(gè)無(wú)向雙連通圖G=(V,E)，其中V和E分別代表節(jié)點(diǎn)與邊；Dijkstra算法在初始搜索過(guò)程中，對(duì)起點(diǎn)i與終點(diǎn)j之間的路徑賦予0值，接著添加E中的一條連接邊至備選路徑并更新路徑長(zhǎng)度，在路徑的搜索中始終保持著松弛操作，即在搜索過(guò)程中可能駁回上次搜索的結(jié)果并重新開始搜索其他可能的連接邊，直至松弛操作結(jié)束后，獲得起始節(jié)點(diǎn)i到終點(diǎn)j由多條邊組成的最短路徑。

鑒于Dijkstra算法能夠?qū)θ我鈨牲c(diǎn)的路徑進(jìn)行篩選，根據(jù)權(quán)值或邊長(zhǎng)規(guī)劃最優(yōu)路徑(邊)集合，因而可用于提煉Delaunay三角網(wǎng)中任意兩個(gè)相干目標(biāo)的連接，即采用權(quán)值更小的邊集合來(lái)代替權(quán)值更大的邊。在之前的研究中，雖然文獻(xiàn)[21]已考慮用更新的網(wǎng)絡(luò)進(jìn)行形變參數(shù)(高程誤差和速率)估算，但是網(wǎng)絡(luò)形態(tài)根據(jù)具體規(guī)則進(jìn)行改變能夠提升解纏精度仍缺乏深入的研究。

2.2 Dijkstra網(wǎng)絡(luò)解纏

前已述及，Delaunay三角剖分準(zhǔn)則是相鄰3個(gè)凸包是否存在公共點(diǎn)以及兩兩之間是否存在公共邊，而未考慮生成的三角形中頂點(diǎn)與頂點(diǎn)間的相位連續(xù)性。由于變形信號(hào)、大氣湍流和相位噪聲的存在，Delaunay網(wǎng)絡(luò)中部分邊并不能滿足先驗(yàn)假設(shè)。這會(huì)引起相位差大于π的概率增大，導(dǎo)致模糊度的錯(cuò)誤估計(jì)。為了得到更好的解纏結(jié)果，需要消除或替換存在違背相位連續(xù)假設(shè)的邊。因此，這個(gè)優(yōu)化過(guò)程可以簡(jiǎn)化為一個(gè)網(wǎng)絡(luò)優(yōu)化問(wèn)題。本文提出利用Dijkstra算法去完成網(wǎng)絡(luò)的優(yōu)化，將點(diǎn)與點(diǎn)之間的相位變化的程度作為邊長(zhǎng)，變化程度越大的邊的邊長(zhǎng)越長(zhǎng)，變化程度越小的邊的邊長(zhǎng)越短。在路徑優(yōu)化中，相位變化程度越小的邊越容易被接受，變化越大的邊越容易被舍去，因此更容易滿足相位連續(xù)性假設(shè)。

對(duì)于時(shí)間序列數(shù)據(jù)，度量空間上兩個(gè)相干目標(biāo)i和j的相位變化程度的有效方法之一是時(shí)間相干性γi,j[21]

(4)

式中，δφi,j,p表示M幅干涉圖中第p幅干涉圖i，j像素間的相位差,在網(wǎng)絡(luò)搜索過(guò)程中，路徑的長(zhǎng)度為-10lg10(γ)。

由式(4)可見，時(shí)間相干性由干涉圖序列同一位置的邊的纏繞相位梯度求復(fù)數(shù)模疊加而來(lái)，相位梯度越小，時(shí)間相干性越接近1，相位梯度越大，時(shí)間相干性越接近0。因?yàn)闀r(shí)間相干性實(shí)際上反映了時(shí)序干涉圖相位梯度的平均水平，所以以時(shí)間相干性作為優(yōu)化準(zhǔn)則的網(wǎng)絡(luò)重構(gòu)方法有利于提高整個(gè)時(shí)序相位解纏的精度。

將Dijkstra算法與時(shí)間相干性結(jié)合實(shí)現(xiàn)網(wǎng)絡(luò)優(yōu)化的過(guò)程可以歸納為兩個(gè)基本步驟：①給定初始的無(wú)向雙連通圖G=(V,E)；②利用Dijkstra搜索給定起點(diǎn)集合S與終點(diǎn)集合D之間的最短路徑并記錄,直到遍歷所有待搜索路徑。在步驟②的搜索過(guò)程中本文將Delaunay三角網(wǎng)上的每條邊的端點(diǎn)作為Dijkstra搜索的起點(diǎn)與終點(diǎn)。由于計(jì)算所有邊的時(shí)間相干性十分耗時(shí)，本文采用KNN(K-nearest neighbor)算法搜索距離每個(gè)點(diǎn)鄰近的50個(gè)離散點(diǎn)來(lái)完成E的構(gòu)建，并將E與Delaunay三角網(wǎng)的邊集合合并生成的新的集合E。使用這種策略的依據(jù)是隨著兩點(diǎn)歐氏距離的減小，相位梯度可能更小。在步驟②中，根據(jù)相干性的高低，任何S與D直接相連的邊都有可能被集合E中的邊代替。

在完成Dijkstra網(wǎng)絡(luò)的搜索后，Delaunay三角網(wǎng)中時(shí)間相干性低的初始邊集合都被更新，替換成了搜索出的最短路徑的邊集合，由此生成了全新的網(wǎng)絡(luò)。值得注意的是，該算法并不能確保三角形完整，那些被取代了邊的三角形可能遭到破壞，也摧毀了MCF解纏三角無(wú)旋性約束，因此這些邊需要針對(duì)性地進(jìn)行三角形重構(gòu)。在重構(gòu)之前，需要根據(jù)簡(jiǎn)單的拓?fù)潢P(guān)系判斷出哪些邊仍不存在任何三角形中，這個(gè)拓?fù)潢P(guān)系可以描述為：如果一條邊所在的兩個(gè)端點(diǎn)，沒(méi)有另外兩條存在公共端點(diǎn)的邊與之相連，則判斷此邊不存在于任何三角形中。提取出這些邊后，根據(jù)每條邊的頂點(diǎn)生成待選邊集合，待選邊由兩個(gè)頂點(diǎn)的KNN生成的鄰近邊組成，選擇待選邊中存在公共端點(diǎn)且邊長(zhǎng)之和最短的兩條邊作為新三角的兩條邊。根據(jù)此思路可以生成的新的三角網(wǎng)絡(luò)并產(chǎn)生新的三角無(wú)旋約束。

新生成的Dijkstra網(wǎng)絡(luò)與Delaunay三角網(wǎng)解纏的過(guò)程一致，均使用基于三角無(wú)旋約束的最小費(fèi)用流方法進(jìn)行解纏，兩者唯一的差異則為網(wǎng)絡(luò)形態(tài)的不一致。因此，對(duì)比兩者網(wǎng)絡(luò)解纏結(jié)果的優(yōu)劣可以判斷哪種網(wǎng)絡(luò)更適用于相位解纏。

3 試驗(yàn)與分析

本文使用模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)分別對(duì)基于Delaunay三角網(wǎng)與Dijkstra網(wǎng)絡(luò)的相位解纏進(jìn)行了評(píng)估。

3.1 模擬數(shù)據(jù)

本文模擬了一個(gè)邊長(zhǎng)為200 km的正方形區(qū)域，在正方形區(qū)域中隨機(jī)選擇了10 000個(gè)離散點(diǎn)，如圖1(a)所示，并模擬了50幅干涉圖，每幅時(shí)間間隔步長(zhǎng)為12 d，模擬的線性變形速率(mm/a)由Matlab中Peaks函數(shù)提供。第50幅干涉圖的形變相位如圖1(a)所示。在此基礎(chǔ)上，根據(jù)分形模型模擬了大氣噪聲(大氣絕對(duì)值相位不超過(guò)π)與相位隨機(jī)噪聲。此外，所有干涉圖的隨機(jī)噪聲的大小設(shè)置為與相干性0.3條件下的相位噪聲相等[22]。第50幅干涉圖大氣分量和噪聲分量的模擬相位分別如圖1(b)和圖1(c)所示。圖1(d)為模擬的包含變形、大氣和噪聲元素的第50幅干涉圖。

(1) 評(píng)估Dijkstra算法對(duì)Delaunay網(wǎng)絡(luò)進(jìn)行更新的能力。與圖1(e)中所展示的Delaunay網(wǎng)絡(luò)的時(shí)間相干性相比，圖1(i)中Dijkstra時(shí)間相干性有所提高。該結(jié)果驗(yàn)證了Dijkstra在尋找更可靠邊的能力。此外，不同網(wǎng)絡(luò)的解纏結(jié)果的準(zhǔn)確性也是評(píng)價(jià)不同網(wǎng)絡(luò)的重要指標(biāo)。圖1(g)和圖1(k)分別展示了第50幅干涉圖的解纏誤差。對(duì)照?qǐng)D1(f)和圖1(j)的解纏結(jié)果可以看出，用Delaunay網(wǎng)絡(luò)進(jìn)行解纏的結(jié)果誤差大于Dijkstra網(wǎng)絡(luò)的解纏結(jié)果，而誤差主要集中在相位梯度大的位置，這說(shuō)明通過(guò)引入時(shí)間相干性來(lái)考慮相位差的Dijkstra網(wǎng)絡(luò)更容易滿足相位連續(xù)性假設(shè)。在圖1(l)和1(h)中試驗(yàn)還比較了整個(gè)時(shí)序干涉圖(50景)中模糊度估計(jì)誤差的絕對(duì)值之和，一個(gè)模糊整數(shù)對(duì)應(yīng)2π解纏誤差。結(jié)果表明，基于Dijkstra網(wǎng)絡(luò)的模糊度估計(jì)誤差小于基于Delaunay網(wǎng)絡(luò)估算的模糊度誤差，尤其是在高相位梯度的區(qū)域。

圖1 模擬干涉圖及Delaunay網(wǎng)絡(luò)與Dijkstra網(wǎng)絡(luò)解纏性能比較Fig.1 Simulated interferograms and the phase unwrapping performance comparison between Delaunay network and Dijkstra network

(2) 進(jìn)行蒙特卡羅模擬試驗(yàn)，測(cè)試Dijkstra相比Delaunay具體可以減少的解纏誤差百分比。在模擬試驗(yàn)中，相干性被設(shè)置為不同的值，范圍在0.05～0.95之間。在每次仿真中都記錄了包含解纏誤差的點(diǎn)的百分比。給定的不同的相干性數(shù)值下，模擬試驗(yàn)重復(fù)1000次。1000次模擬試驗(yàn)的平均結(jié)果如圖2所示。由圖2可以看出，在大多數(shù)相干情況下，Dijkstra網(wǎng)絡(luò)的解纏誤差均比Delaunay網(wǎng)絡(luò)小約3%～8%，這驗(yàn)證了網(wǎng)絡(luò)質(zhì)量對(duì)相位解纏的重要性。當(dāng)相干性接近1時(shí)，表明Delaunay網(wǎng)絡(luò)上的邊都具有很高的相干性，因此Dijkstra算法更新前后網(wǎng)絡(luò)幾乎無(wú)差異，即當(dāng)前的邊集合就是最優(yōu)的路徑。

圖2 各相干性條件下，兩種網(wǎng)絡(luò)的解纏誤差比例Fig.2 Percentages of points with unwrapping errors

3.2 真實(shí)數(shù)據(jù)

本文選取COMET LICS[23]平臺(tái)處理的Sentinel-1A Track128上93幅未經(jīng)濾波的小基線集干涉圖(子集最長(zhǎng)時(shí)間基線36 d，如圖3中時(shí)空基線圖所示)作為真實(shí)數(shù)據(jù)，影像觀測(cè)時(shí)間覆蓋2019-01-09—2020-02-09，影像觀測(cè)時(shí)間間隔為12 d。Track128影像覆蓋區(qū)域見圖3，圖3(b)為切割后試驗(yàn)區(qū)的Sentinel-2光學(xué)影像，該地區(qū)覆蓋著大面積植被，相干性較低且地形較為復(fù)雜，適合作為判斷解纏方法性能的區(qū)域。在選取93幅干涉圖后，利用LICS平臺(tái)提供的所有干涉圖的相干性計(jì)算出平均相干性，選取平均相干性大于0.25的點(diǎn)作為相干目標(biāo)點(diǎn)，進(jìn)行構(gòu)網(wǎng)和解纏。

圖3 研究區(qū)域概況及干涉圖Fig.3 Overview of the study region and the general information of the used interferograms

圖4展示了Delaunay三角網(wǎng)與優(yōu)化后的Dijkstra網(wǎng)絡(luò)基于93幅干涉圖計(jì)算出的時(shí)間相干性?？梢钥闯鯠elaunay三角網(wǎng)中時(shí)間相干性較低的區(qū)域在經(jīng)過(guò)Dijkstra優(yōu)化后總體上得到了提升，由圖4(a)與圖4(b)中局部區(qū)域的縮放圖可以看出，時(shí)間相干性較低的邊被周圍相干性較高的邊所取代，時(shí)間相干性較高的邊則被進(jìn)一步保留，因此可以進(jìn)一步說(shuō)明Dijkstra網(wǎng)絡(luò)在生成過(guò)程中，對(duì)時(shí)間相干性較低的邊進(jìn)行了更新，時(shí)間相干性較低的邊會(huì)被兩條或兩條以上的邊所替代，而原本時(shí)間相干性較高的邊則不會(huì)被替換從而被保留了下來(lái)。因此，最后生成的Dijkstra網(wǎng)絡(luò)較Delaunay網(wǎng)絡(luò)整體時(shí)間相干性有所上升，這也直接反映出在Delaunay三角網(wǎng)中相位差較大的邊已經(jīng)被減少，并被相位差較小的邊所替代。

圖5給出了2019-08-01—2019-09-06與2019-10-12—2019-11-05兩期干涉圖解纏的結(jié)果。從目視解譯可知，Delaunay解纏結(jié)果(圖5(c)—(d))中存在較多相位不連續(xù)的區(qū)域，這些區(qū)域表現(xiàn)為與鄰近區(qū)域相差2π，說(shuō)明這些不連續(xù)的區(qū)域存在解纏誤差，且這些區(qū)域大多為圖4中時(shí)間相干性較低的區(qū)域。而Dijkstra解纏后相位(圖5(e)—(f))在空間上比較連續(xù)，相位不連續(xù)的區(qū)域被減少。值得注意的是，在Delaunay解纏結(jié)果中存在相位不連續(xù)的區(qū)域，在Dijkstra的解纏結(jié)果中存在較少與鄰近區(qū)域有2π跳變的情況，所以可以判斷Dijkstra網(wǎng)絡(luò)時(shí)間相干性的優(yōu)化對(duì)相位解纏起到了積極的作用。

圖5 纏繞相位圖與兩種網(wǎng)絡(luò)解纏相位圖Fig.5 Wrapped and unwrapped interferograms

定量方面，本文沿用當(dāng)前解纏誤差檢驗(yàn)的方法，通過(guò)統(tǒng)計(jì)時(shí)間維的三角閉合差作為解纏誤差的評(píng)價(jià)標(biāo)準(zhǔn)[24-25]。相位三角閉合差Δψ表示時(shí)空基線平面(圖3(c))形成三角閉合環(huán)的3個(gè)干涉圖解纏后的相位差，閉合差對(duì)應(yīng)的模糊度估計(jì)錯(cuò)誤個(gè)數(shù)μ[25-26]

(5)

為方便表述，這里將處于相位三角中的干涉圖分別用數(shù)字1、2、3表示。另外，本文還使用解纏后干涉序列反算出的時(shí)間序列相位與解纏序列的殘差來(lái)衡量解纏精度，其值越小，說(shuō)明解纏越可靠。

圖6展示了在兩種量化標(biāo)準(zhǔn)下各方法的評(píng)估結(jié)果。相比圖6(a)中Delaunay三角網(wǎng)的模糊度誤差個(gè)數(shù)，圖6(b)中Dijkstra網(wǎng)絡(luò)的模糊度誤差個(gè)數(shù)減少，且減少區(qū)域?qū)?yīng)時(shí)間相干性提升的區(qū)域，因此可以判斷考慮相位梯度，并以此改善網(wǎng)絡(luò)形態(tài)的Dijkstra網(wǎng)絡(luò)的確能夠更好地滿足相位連續(xù)性假設(shè)，從而減少解纏誤差。對(duì)比模糊度誤差個(gè)數(shù)直方圖，Delaunay峰值為45而Dijkstra為30。相比而言，Dijkstra網(wǎng)絡(luò)減少了約33%的模糊度解算錯(cuò)誤造成的不閉合三角環(huán)數(shù)。此外，對(duì)比兩者解纏結(jié)果的時(shí)間相干性及其統(tǒng)計(jì)直方圖可以發(fā)現(xiàn)，Dijkstra網(wǎng)絡(luò)相比Delaunay三角網(wǎng)整體時(shí)間相干性更高，整體相干性提升約0.05，這也代表提升時(shí)間相干性的Dijkstra網(wǎng)絡(luò)更有益于InSAR時(shí)間序列的解算。

圖6 兩種網(wǎng)絡(luò)解纏結(jié)果中存在的模糊度錯(cuò)誤個(gè)數(shù)以及對(duì)應(yīng)解纏結(jié)果估計(jì)出的時(shí)間相干性Fig.6 The sum of absolute value of ambiguity error and temporal coherence respectively for Delaunay network and Dijkstra network

圖7統(tǒng)計(jì)了兩種網(wǎng)絡(luò)在解纏過(guò)程中的費(fèi)用和解纏效率。對(duì)于93幅干涉圖的解纏，Delaunay三角網(wǎng)消耗的費(fèi)用(圖7(a))均高于Dijkstra網(wǎng)絡(luò)，相應(yīng)的解纏消耗時(shí)間(圖7(b))也高于Dijkstra。這兩者的增益均來(lái)自Dijkstra網(wǎng)絡(luò)中高相位梯度的邊被減少，因此需要解算的存在模糊度的邊也被減少，相對(duì)應(yīng)解纏消耗時(shí)間則減少。從這些試驗(yàn)中可以得出如下結(jié)論：就相位解纏而言，Delaunay三角網(wǎng)從精度與效率方面均非最優(yōu)選網(wǎng)絡(luò)，Dijkstra網(wǎng)絡(luò)在這兩個(gè)方面均優(yōu)于Delaunay三角網(wǎng)。

圖7 Delaunay三角網(wǎng)與Dijkstra網(wǎng)絡(luò)的解纏費(fèi)用與消耗的時(shí)間對(duì)比Fig.7 The total cost and computational time for unwrapping based on Delaunay network and Dijkstra network

由于選取區(qū)域缺乏實(shí)測(cè)數(shù)據(jù)進(jìn)行InSAR結(jié)果的驗(yàn)證，且解纏整體的提升效果難以通過(guò)點(diǎn)的對(duì)比進(jìn)行體現(xiàn)，因此，本文進(jìn)行了進(jìn)一步試驗(yàn)以說(shuō)明解纏對(duì)形變解算所產(chǎn)生的影響，本文計(jì)算了從小基線集相位到單主影像轉(zhuǎn)換過(guò)程中的均方根誤差，對(duì)兩種網(wǎng)絡(luò)解纏相位的均方根誤差進(jìn)行了對(duì)比。因?yàn)閮烧呦辔恢挥薪饫p網(wǎng)絡(luò)的不同，因此均方根誤差的大小變化可以認(rèn)定為是由解纏誤差所帶來(lái)的。均方根結(jié)果如圖8所示為表征解纏結(jié)果對(duì)形變解算的影響。

圖8(a)、(b)分別為Delaunay與Dijkstra網(wǎng)絡(luò)解纏結(jié)果的均方根誤差，圖8(c)為兩者結(jié)果的直方圖對(duì)比，根據(jù)兩者的對(duì)比可以看出，Dijkstra網(wǎng)絡(luò)小基線集相位到單主影像轉(zhuǎn)換的模型均方根誤差明顯更小，因此可以確定Dijkstra網(wǎng)絡(luò)對(duì)于相位解纏精度的提升有益于形變量的解算。

圖8 小基線集與單主影像相位轉(zhuǎn)換過(guò)程中的均方根誤差Fig.8 RMSE of single-reference phase converted from SBAS phase

4 結(jié) 論

本文針對(duì)InSAR時(shí)間序列相位解纏中常用的Delaunay三角網(wǎng)性能進(jìn)行了評(píng)估，并提出了Dijkstra網(wǎng)絡(luò)與之進(jìn)行了解纏性能的對(duì)比。Dijkstra網(wǎng)絡(luò)利用相位梯度計(jì)算時(shí)間相干性，并將時(shí)間相干性作為衡量路徑長(zhǎng)度的指標(biāo)，在此基礎(chǔ)上，利用Dijkstra算法對(duì)Delaunay三角網(wǎng)進(jìn)行優(yōu)化，以低相位梯度的邊取代了高相位梯度的邊，最大化了網(wǎng)絡(luò)的時(shí)間相干性，從而逼近了相位連續(xù)性假設(shè)。模擬和真實(shí)試驗(yàn)均表明Delaunay三角網(wǎng)并不是相位解纏的最優(yōu)選擇，Dijkstra網(wǎng)絡(luò)在精度與效率方面具有更好的表現(xiàn)。尤其在精度方面，Dijkstra網(wǎng)絡(luò)的解纏模糊度誤差個(gè)數(shù)較Delaunay三角網(wǎng)減少接近約33%。此外本文研究結(jié)果證明空間參考網(wǎng)絡(luò)的質(zhì)量?jī)?yōu)化也是提升解纏性能的關(guān)鍵步驟之一。

致謝：感謝ESA/Copernicus提供的Sentinel-1數(shù)據(jù)。