劉明泰,岳東杰

(河海大學(xué),江蘇 南京 210098)

0 引 言

近些年來(lái)相位圖像在各個(gè)領(lǐng)域得到廣泛應(yīng)用,如干涉合成孔徑雷達(dá)InSAR,磁共振成像 MRI以及激光雷達(dá)等相位成像技術(shù),然而它們?cè)趹?yīng)用過(guò)程中都面臨一個(gè)共同的問(wèn)題——相位解纏。傳統(tǒng)相位解纏算法包括路徑跟蹤法(枝切法[1]、質(zhì)量圖引導(dǎo)法[2]、最小不連續(xù)法[3])、最小范數(shù)法 (最小二乘法[4])、基于網(wǎng)絡(luò)流的優(yōu)化方法(最小費(fèi)用流法[5-6])等,都難以在噪聲和形變較大的區(qū)域恢復(fù)出真實(shí)的地形和形變相位。隨著深度學(xué)習(xí)模型的發(fā)展,在圖像處理領(lǐng)域提出了越來(lái)越多的模型如CNN、RNN、GAN、 Transformer等,這些模型也被用于InSAR的相位解纏研究。Spoorthietal等[7]認(rèn)為相位解纏是一個(gè)語(yǔ)義分割問(wèn)題,即預(yù)測(cè)每個(gè)像素的纏繞數(shù)量(整數(shù)為2π),然后將其添加到纏繞相位中以獲得最終解,通過(guò)卷積神經(jīng)網(wǎng)絡(luò)CNN構(gòu)建PhaseNet模型來(lái)求解纏繞整周數(shù),恢復(fù)展開(kāi)相位。Zhou Lifan等[8]將GAN 用于相位解纏,將解纏視作圖像到圖像的轉(zhuǎn)換問(wèn)題,通過(guò)對(duì)抗學(xué)習(xí)來(lái)訓(xùn)練相位解纏模型 PUGAN。Wu Zhipeng等[9]提出了適用于礦區(qū)快速挖掘引起快速變形的解纏網(wǎng)絡(luò) PUNet[9-10]。近幾年隨著擴(kuò)散模型DDPM[11]的提出,其在AI繪畫(huà)領(lǐng)域取得了許多令人驚嘆的研究成果,如stableDiffusion、NovelAI、Waifu 等。這些成果無(wú)不論證了擴(kuò)散模型的潛力?；诖吮疚膶U(kuò)散生成模型引入相位解纏領(lǐng)域,并與傳統(tǒng)相位解纏算法、基于深度學(xué)習(xí)的解纏模型進(jìn)行比較,從而評(píng)價(jià)模型的解纏能力。

1 擴(kuò)散概率模型DDPM

(1)

圖1 干涉相位圖擴(kuò)散過(guò)程

(2)

(3)

給定足夠的T和表現(xiàn)良好的噪聲表βt,最終時(shí)刻的潛在變量xt幾乎是各相同性的高斯分布。因此,如果知道確切的反向分布q(xt-1|xt),就可以對(duì)xt～N(0,1)進(jìn)行采樣,并反向運(yùn)行該過(guò)程,最終從q(x0)中獲得結(jié)果。上述反向運(yùn)算的過(guò)程稱為逆擴(kuò)散過(guò)程,由于逆向分布q(xt-1|xt)取決于整個(gè)數(shù)據(jù)分布,因此可使用神經(jīng)網(wǎng)絡(luò)來(lái)擬合每步的逆向分布,如公式(4)所示。

Pθ(xt-1|xt)=N[xt-1;uθ(xt,t),∑θ(xt,t)]

(4)

擴(kuò)散模型通過(guò)優(yōu)化目標(biāo)函數(shù)的下界(ELBO)來(lái)最小化負(fù)對(duì)數(shù)似然函數(shù),可以將變分下界寫(xiě)為:

LELBO=L0+L1+...+Lt

(5)

L0=-logpθ(x0|x1)

(6)

Lt-1=Dkl[q(xt-1|xt,x0)||pθ(xt-1|xt)]

(7)

Lt=Dkl[q(xt|x0)||p(xt)]

(8)

除了常數(shù)L0,ELBO的每一項(xiàng)都是兩個(gè)高斯分布之間的KL散度,其中Lt項(xiàng)并不依賴于神經(jīng)網(wǎng)絡(luò)中的參數(shù)θ,因此,可以只對(duì)Lt-1項(xiàng)進(jìn)行化簡(jiǎn),雖真實(shí)分布q(xt-1|xt)未知,但加上條件x0的后驗(yàn)分布q(xt-1|xt,x0)可以通過(guò)貝葉斯公式得到公式(8),最終優(yōu)化目標(biāo)函數(shù)的下界寫(xiě)作各個(gè)步驟t的后驗(yàn)高斯分布與神經(jīng)網(wǎng)絡(luò)擬合的高斯分布之間的KL散度,如下所示。

(9)

為了便于計(jì)算,Ho等[11]將網(wǎng)絡(luò)擬合的分布P的方差設(shè)為常數(shù),根據(jù)高斯分布之間的KL散度公式,后驗(yàn)分布q(xt-1|xt,x0)與網(wǎng)絡(luò)擬合的分布Pθ(xt-1|xt)之間的KL散度可以寫(xiě)作公式(11),并且使用網(wǎng)絡(luò)去預(yù)測(cè)均值uθ(xt,t)的效果不如預(yù)測(cè)噪音ε的效果好,噪音ε與均值uθ的關(guān)系如下所示。

(10)

(11)

(12)

Lsimple=Et,x0{[ε-εθ(xt-t)]2}

(13)

2 實(shí)驗(yàn)數(shù)據(jù)以及訓(xùn)練、推理流程

在有噪聲且快速變形的干涉圖中,可靠的相位解纏是困難的,采集足夠的地面測(cè)量數(shù)據(jù)進(jìn)行訓(xùn)練是不可行的。因此,需要生成與真實(shí)干涉圖特征接近的合成干涉圖,并將其用于訓(xùn)練。干涉圖的相位可以分解為地形相位、形變相位、大氣湍流相位和噪聲。本文采用了兩種方法來(lái)獲取干涉圖數(shù)據(jù),一種是構(gòu)建干涉圖模擬器,通過(guò)干涉圖模擬器來(lái)模擬生成干涉相位圖數(shù)據(jù),另一種為使用LiCSAR系統(tǒng)自動(dòng)處理生成的Sentinel-1干涉圖數(shù)據(jù),訓(xùn)練樣本大小設(shè)置為256×256像素,覆蓋面積約為5 km×5 km,小樣本可以加快訓(xùn)練過(guò)程。

2.1 LiCSAR數(shù)據(jù)

LiCSAR的干涉圖產(chǎn)品包括干涉圖、相干圖、解纏圖,LiCSAR干涉圖的數(shù)據(jù)以Frame為單位(250 km×250 km),提供tif和png格式的數(shù)據(jù),空間分辨率110 m,為短時(shí)間基線干涉圖(每張影像與其后續(xù)3張影像進(jìn)行干涉),如圖2所示。由于訓(xùn)練網(wǎng)絡(luò)所需要數(shù)據(jù)的像素大小為256×256,因此要對(duì)下載得到的數(shù)據(jù)進(jìn)行預(yù)處理,對(duì)數(shù)據(jù)進(jìn)行旋轉(zhuǎn)和裁切,如圖3所示。

圖2 LiCSAR數(shù)據(jù)

圖3 預(yù)處理后的LiCSAR

2.2 干涉圖模擬器

干涉圖模擬器通過(guò)輸入隨機(jī)裁剪256×256像素的DEM地形數(shù)據(jù)用于生成地形相位圖4(a)。除此之外,使用二維高斯曲面和5×5的隨機(jī)向量生成二維扭曲高斯曲面來(lái)模擬形變相位圖4(b),而大氣相位圖4(c)可以使用分形柏林噪聲進(jìn)行模擬。最后,考慮到噪聲相位的影響,可以通過(guò)添加高斯噪聲進(jìn)行模擬。從而生成絕對(duì)相位圖4(d)并進(jìn)行相位纏繞運(yùn)算,得到干涉相位圖4(e),其流程如圖4所示。

圖4 模擬干涉相位圖

2.3 訓(xùn)練推理流程

基于去噪概率擴(kuò)散的相位解纏模型在訓(xùn)練過(guò)程中,根據(jù)真實(shí)相位x0以及t時(shí)刻隨機(jī)采樣的噪聲值εt來(lái)獲得擴(kuò)散過(guò)程中t時(shí)刻經(jīng)過(guò)加噪后的值xt,將xt和纏繞相位以及t時(shí)刻輸入網(wǎng)絡(luò)中得到預(yù)測(cè)的噪聲εθ,t,根據(jù)公式(13),計(jì)算εt和εθ,t的均方誤差MSE,來(lái)反向梯度傳播優(yōu)化網(wǎng)絡(luò)。在推理過(guò)程中,將纏繞相位和隨機(jī)采樣的噪聲值xt以及最終時(shí)刻T輸入網(wǎng)絡(luò),獲得T-1時(shí)刻的采樣值xT-1,之后將纏繞相位和xT-1以及T-1時(shí)刻輸入網(wǎng)絡(luò)獲得上一時(shí)刻的值。就這樣逐步采樣,最終獲得預(yù)測(cè)的真實(shí)相位,其訓(xùn)練和推理流程如圖5和表1所示。

表1 訓(xùn)練和推理流程

圖5 訓(xùn)練和推理流程

3 實(shí)驗(yàn)結(jié)果及分析

3.1 評(píng)測(cè)指標(biāo)

由于相位解纏的結(jié)果是絕對(duì)相位中各像素之間的相對(duì)值,其結(jié)果并非唯一解(即解纏相位與真實(shí)相位相差一個(gè)常數(shù),其解纏結(jié)果依然能準(zhǔn)確反映出地形的起伏變化),采用原始相位和解纏相位之間的均方誤差并不能很好反映解纏結(jié)果的精度,因此本文采用原始相位和解纏相位之間差值的方差作為評(píng)測(cè)的L1指標(biāo),采用原始相位的相位梯度和解纏相位的相位梯度作為評(píng)測(cè)的L2指標(biāo):

L1=E[(Ψ-Ψθ)2]-E2(Ψ-Ψθ)

(14)

(15)

Gver=Ψi,j+1-Ψi,j

(16)

Ghor=Ψi+1,j-Ψi,j

(17)

3.2 不同相位解纏方法的結(jié)果比較

在相位解纏領(lǐng)域,有許多不同的方法可用于解纏干涉相位。在本文的實(shí)驗(yàn)中,將本文提出的基于去噪擴(kuò)散概率模型的相位解纏模型與4種傳統(tǒng)的相位解纏方法,包括枝切法、質(zhì)量圖引導(dǎo)法、最小二乘法、最小費(fèi)用流相位解纏方法以及3種深度學(xué)習(xí)相位解纏模型PhaseNet2.0、PUNet、PU-GAN進(jìn)行比較。本實(shí)驗(yàn)采樣用3 000個(gè)樣本對(duì)這些方法進(jìn)行了測(cè)試,并對(duì)解纏結(jié)果進(jìn)行比較和分析,其中3 000幅干涉相位圖樣本組成的測(cè)試數(shù)據(jù)集在L1指標(biāo)下的誤差分布如圖6所示,在L2指標(biāo)下的誤差分布如圖7所示。

圖6 數(shù)據(jù)集上在L1指標(biāo)的相位展開(kāi)誤差分布

圖7 數(shù)據(jù)集上在L2指標(biāo)的相位展開(kāi)誤差分布

結(jié)果表明,DDPU的解纏結(jié)果在L1指標(biāo)下,測(cè)試樣本的誤差σ多分布于0～2的區(qū)間內(nèi),|σ|<2的樣本占86.53%,相比之下傳統(tǒng)相位解纏算法中的枝切法和質(zhì)量圖引導(dǎo)法在不同信噪比的樣本測(cè)試中,噪聲的影響導(dǎo)致效果十分不理想,解纏結(jié)果的誤差大多分布于4～16的區(qū)間,最小二乘法的效果表現(xiàn)一般,解纏效果優(yōu)于枝切法和質(zhì)量圖引導(dǎo)法,較差于最小費(fèi)用流方法,其解纏結(jié)果的誤差多分布于2～8的區(qū)間內(nèi)。最小費(fèi)用流的相位解纏方法是4種傳統(tǒng)相位解纏算法中表現(xiàn)最好的算法,在一些情況下甚至優(yōu)于個(gè)別基于深度學(xué)習(xí)的相位解纏算法。在深度學(xué)習(xí)相位解纏算法中,DDPU的解纏效果要明顯優(yōu)于PUNet、PhaseNet和PU-GAN,DDPU通過(guò)抑制大偏差的發(fā)生,顯著提高了展開(kāi)干涉圖的可靠性。

在L2指標(biāo)下,測(cè)試樣本的解纏相位梯度和原始相位梯度之間的均方誤差多分布于0～0.1和0.1～0.2的區(qū)間內(nèi),占比83.2%。根據(jù)上述結(jié)果,可以看出使用對(duì)抗生成模型的PU-GAN解纏所得的數(shù)據(jù)較為平滑,不會(huì)出現(xiàn)較大的誤差,因此其解纏結(jié)果的相位梯度多分布于中間值0.1～0.4的區(qū)間內(nèi),基于擴(kuò)散概率的相位解纏模型DDPU不僅在整體上能更好符合原始相位,在梯度變化處也能更好地恢復(fù)出原始相位即大多數(shù)樣本處于0～0.1這個(gè)低誤差區(qū)間內(nèi)。

根據(jù)3 000個(gè)樣本在指標(biāo)L1和指標(biāo)L2下8種方法的解纏結(jié)果誤差以及最大值和最小值,繪制箱型圖來(lái)分析各個(gè)方法的魯棒性,如圖8所示。

圖8 數(shù)據(jù)集上在L1指標(biāo)L2指標(biāo)的誤差箱型圖

由圖8可知,DDPU方法所得到的解纏誤差分布相對(duì)于其他方法來(lái)說(shuō)更加緊湊,即四分位距(上下四分位之差)更小。這表明DDPU方法在不同信噪比條件下具有更強(qiáng)的穩(wěn)定性和魯棒性,能夠有效地抑制大誤差的發(fā)生。相比之下,其他方法的解纏誤差分布比較散亂,四分位距較大,說(shuō)明這些方法在不同程度的噪聲情況下容易出現(xiàn)大幅度的偏差,穩(wěn)定性相對(duì)較差。因此,DDPU方法在實(shí)際應(yīng)用中具有更高的可靠性和準(zhǔn)確性,能夠更好地滿足解纏需求。而且,DDPU解纏誤差的均值和中位數(shù)也小于其他解纏方法,這表明DDPU的解纏結(jié)果更加準(zhǔn)確。

為了更加直觀展示DDPU的解纏效果,采用三幅不同噪聲水平的模擬干涉圖,并進(jìn)行了相位展開(kāi)處理,其結(jié)果如圖9所示,相位展開(kāi)誤差如圖10所示。其中,樣本2的噪聲水平低,條紋清晰,使得以上各種方法都能取得較好的解纏結(jié)果,然而隨著信噪比的下降即干涉相位中噪聲水平的提高如樣本1和樣本3,傳統(tǒng)方法中的枝切法、質(zhì)量圖引導(dǎo)法便很難獲得正確的解纏結(jié)果,最小二乘法的解纏誤差覆蓋了整個(gè)干涉圖,最小費(fèi)用流和其他深度學(xué)習(xí)相位解纏算法在地形起伏大的區(qū)域均有一定幅度的解纏誤差,而本文提出的DDPU模型相對(duì)于其他相位解纏方法可以克服噪聲的影響,精確地從干涉圖中恢復(fù)出解纏相位。

圖9 不同噪聲水平的干涉圖實(shí)驗(yàn)結(jié)果

圖10 不同方法解纏結(jié)果的誤差

3.3 DDPU模型評(píng)估

為了評(píng)估DDPU相位解纏模型對(duì)于不同噪聲水平和不同相位梯度干涉相位圖的解纏效果,本文設(shè)計(jì)了一系列實(shí)驗(yàn)來(lái)對(duì)其進(jìn)行全面的評(píng)測(cè),旨在深入了解DDPU模型在不同噪聲水平和相位梯度情況下的魯棒性和性能表現(xiàn),以驗(yàn)證其在實(shí)際應(yīng)用中的可靠性和適用性。將不同噪聲水平的測(cè)試樣本傳入DDPU模型中逐步獲得解纏相位,使用評(píng)測(cè)指標(biāo)L1對(duì)DDPU模型在不同噪聲和相位梯度條件下的解纏效果進(jìn)行定量評(píng)估和比較,其結(jié)果如圖11所示。

圖11 DDPU對(duì)不同噪聲水平樣本的解纏誤差

結(jié)果表明,四組不同噪聲水平的干涉圖樣本隨著最大相位梯度的增大,起初,相位展開(kāi)誤差的標(biāo)準(zhǔn)差分布較為密集且都集中在較小的區(qū)域,之后隨著最大相位梯度的增大,DDPU的解纏誤差開(kāi)始增大并且分布的更加離散。此外,DDPU表現(xiàn)出優(yōu)秀的抗噪聲能力,考慮到可接收的誤差范圍為小于1 rad,當(dāng)噪聲水平小于0.8π時(shí),對(duì)于最大相位梯度小于4rad的干涉相位圖,DDPU均能夠獲得良好的解纏結(jié)果。

4 結(jié) 語(yǔ)

將計(jì)算機(jī)視覺(jué)領(lǐng)域中的擴(kuò)散生成模型引入相位解纏,通過(guò)與現(xiàn)有的幾種相位解纏方法進(jìn)行比較,發(fā)現(xiàn)擴(kuò)散生成模型即使在嚴(yán)重的噪聲條件下也能達(dá)到優(yōu)秀的解纏性能,其對(duì)噪聲的抵抗力非常適合于低信噪比的相位數(shù)據(jù),且需要精確相位解纏結(jié)果的應(yīng)用。展望未來(lái),我們希望擴(kuò)展該網(wǎng)絡(luò),以解決特定應(yīng)用的挑戰(zhàn)。