藺祥楠， 禹衛(wèi)東

(1.中國(guó)科學(xué)院電子學(xué)研究所， 北京 100190； 2.中國(guó)科學(xué)院大學(xué)， 北京 100049)

0 引言

干涉合成孔徑雷達(dá)(Interferometric Synthetic Aperture Radar, InSAR)是通過(guò)計(jì)算兩幅配準(zhǔn)后的SAR復(fù)圖像的干涉相位差，結(jié)合地表高度和干涉相位差之間的關(guān)系，獲取地表三維信息的重要遙感技術(shù)[1]。傳統(tǒng)單頻率相位解纏前提是相鄰像素的絕對(duì)相位差小于π。實(shí)際就是要求待解纏的區(qū)域比較平緩，沒(méi)有陡變地形。此外，單頻率相位解纏算法計(jì)算方法復(fù)雜，且噪聲敏感性高，誤差傳遞現(xiàn)象嚴(yán)重。因而傳統(tǒng)單頻干涉SAR無(wú)法有效解決某些復(fù)雜地區(qū)的高程重建問(wèn)題。為了實(shí)現(xiàn)復(fù)雜地形的DEM(Digital Elevation Model)反演，提出了多通道InSAR測(cè)量技術(shù)。

多通道InSAR測(cè)量技術(shù)[2]包括多頻率InSAR測(cè)量技術(shù)[3]和多基線InSAR測(cè)量技術(shù)[4]。多頻率InSAR通常選取互質(zhì)的載頻波段來(lái)增大高度模糊數(shù)。多基線InSAR技術(shù)則是利用不同的基線來(lái)增大高度模糊數(shù)。模糊高度數(shù)的擴(kuò)大，有利于得到更為精確的解纏相位信息。

多通道干涉SAR技術(shù)發(fā)展迅速，高度重建算法種類(lèi)眾多。其中較為早期的一種方法是多通道最大似然估計(jì)(Maximum Likelihood Estimation, MLE)算法[5]。該算法利用干涉相位的概率統(tǒng)計(jì)特性直接進(jìn)行高程重建，重建速度快，但受干涉通道數(shù)目、基線長(zhǎng)度、工作頻率、相干系數(shù)的影響特別大，抗噪能力弱，在相干系數(shù)低的區(qū)域需要增加很多的通道數(shù)來(lái)保證精確度[6-7]。為了解決MLE中過(guò)高的通道數(shù)要求和重建精度的矛盾，基于Markov隨機(jī)模型的MRF-MAP(Markov Random Field Based Maximum a Posteriori Estimation)算法被提出來(lái)。該方法利用Markov模型作為先驗(yàn)分布，通過(guò)估計(jì)每個(gè)像素的超參數(shù)實(shí)現(xiàn)高程重建。至今，有許多人提出了改進(jìn)的MAP算法[8-9]。雖然MRF-MAP提高了反演精度，但該方法需估算每一點(diǎn)的超參數(shù)，運(yùn)算時(shí)間過(guò)長(zhǎng)，設(shè)備要求過(guò)高[10]。相比之下，ML算法在計(jì)算速度上擁有絕對(duì)的優(yōu)勢(shì)，適合做大數(shù)據(jù)量的數(shù)據(jù)處理。但ML估計(jì)結(jié)果噪聲點(diǎn)和噪聲聚集成斑塊的現(xiàn)象嚴(yán)重，本文提出了ML估計(jì)結(jié)果噪聲斑塊去除方法。具體方法是：首先利用聚類(lèi)分析[11](CA)思想去掉解纏結(jié)果中大量出錯(cuò)區(qū)域?qū)?yīng)的模糊類(lèi)；再利用變長(zhǎng)滑動(dòng)加權(quán)窗口判定并去除噪斑，得到正確的解纏相位。

1 基于MLE的雙頻干涉高程重建方法

雙頻系統(tǒng)可以通過(guò)兩個(gè)互質(zhì)的頻率增大高度模糊數(shù)，從而實(shí)現(xiàn)更復(fù)雜地形的高程反演。在基于MLE的高程重建算法中，高程重建問(wèn)題被轉(zhuǎn)化成了找到高程值h(p)使似然函數(shù)達(dá)到最大化的問(wèn)題。

在單通道InSAR中，對(duì)于圖像的每個(gè)像素點(diǎn)(i,j)，纏繞的干涉相位φ′與地表高程信息之間的關(guān)系為

(1)

式中，λ為載頻波長(zhǎng)，〈·〉2π為“模2π”處理操作，μ為一個(gè)比例系數(shù)，由InSAR系統(tǒng)收發(fā)模式以及飛行參數(shù)確定：

(2)

式中，B⊥為垂直斜距平面的基線分量，r0為該點(diǎn)斜距，θ為該點(diǎn)斜視角。

根據(jù)InSAR干涉相位的統(tǒng)計(jì)特性，得到干涉相位的概率密度統(tǒng)計(jì)特性：

f(φ|φ′)=

(3)

將式(1)代入式(3)可以得到測(cè)量干涉相位與高程信息之間的概率密度分布函數(shù)：

f(φ|h)=

(4)

式中，γ稱(chēng)作干涉圖像的相關(guān)系數(shù)：

(5)

(6)

(7)

為了解決無(wú)數(shù)個(gè)最大似然解的問(wèn)題，增加獨(dú)立的觀測(cè)相位，采用雙頻干涉技術(shù)，式(4)中的各個(gè)變量加上角標(biāo)則得到了各個(gè)通道對(duì)應(yīng)的似然函數(shù)。此時(shí)，總的似然函數(shù)變成

(8)

利用聯(lián)合最大似然估計(jì)方法進(jìn)行高程重建，多頻率最大似然估計(jì)高程值如下：

(9)

(a) 通道1似然函數(shù)

(b) 通道2似然函數(shù)

(c) 通道3似然函數(shù)

(d) 多頻干涉系統(tǒng)似然函數(shù)圖1 單頻似然函數(shù)和多頻似然函數(shù)示意圖

2 MLE重建結(jié)果噪斑去除方法

MLE的雙頻高程重建算法魯棒性不強(qiáng)，得到的結(jié)果中有單點(diǎn)噪聲和噪聲斑。由文獻(xiàn)[6]可知，這些噪點(diǎn)高程都與周?chē)c(diǎn)的高程相差很大。此外，若高程范圍選擇得當(dāng)，MLE重建結(jié)果中噪聲像素所占比例較小?；谶@兩點(diǎn)，本文提出了一種噪聲斑塊去除方法。能在保證處理效率的同時(shí)，能夠提高重建精度。具體去噪方法分為兩個(gè)部分。

2.1 變鄰域錯(cuò)點(diǎn)標(biāo)記

根據(jù)文獻(xiàn)[11]，將像素的模糊矢量定義為該像素在各個(gè)干涉通道下的模糊數(shù)所組成的矢量。各個(gè)干涉通道下p像素的模糊數(shù)為

(10)

式中，?·」表示下取整運(yùn)算。這樣，p像素所對(duì)應(yīng)的模糊矢量即為

V(p)=[n1(p),n2(p),…,nk(p)]

(11)

根據(jù)聚類(lèi)分析的思想，將高程值按高度模糊數(shù)進(jìn)行聚類(lèi)分組。通過(guò)執(zhí)行以下幾個(gè)步驟來(lái)進(jìn)行錯(cuò)點(diǎn)標(biāo)記：

1) 構(gòu)建初始標(biāo)記矩陣，大小與圖像大小相同，初始值均為1。

2) 根據(jù)目標(biāo)地段高度范圍粗略計(jì)算出高度模糊區(qū)間，區(qū)間之外的點(diǎn)判定為錯(cuò)點(diǎn)，標(biāo)識(shí)符標(biāo)記為0。標(biāo)記出明顯出錯(cuò)的解纏壞斑。

3) 若某一個(gè)模糊類(lèi)所對(duì)應(yīng)的像素個(gè)數(shù)低于預(yù)期比例閾值，判定為噪聲并將該分組內(nèi)的像素標(biāo)識(shí)符標(biāo)記為0。

4) 通過(guò)目標(biāo)點(diǎn)高程與3×3鄰域內(nèi)重建正確點(diǎn)的高程均值的歸一化差值來(lái)判斷該點(diǎn)是否為噪聲點(diǎn)。將噪點(diǎn)像素的標(biāo)記置為0。

當(dāng)相鄰的一片像素點(diǎn)均為壞點(diǎn)時(shí)，3鄰域均值判定噪點(diǎn)的方法不再起作用。鑒于噪聲斑的大小不一，采用可變鄰域的判別器來(lái)判定突變?cè)肼暟撸瑥亩_(dá)到將連城小片的噪聲斑進(jìn)行標(biāo)記的目的。

2.2 非均勻加權(quán)濾波

前面的步驟已經(jīng)將誤差較大的錯(cuò)點(diǎn)標(biāo)記出來(lái)。傳統(tǒng)均值濾波會(huì)將錯(cuò)點(diǎn)帶入到加權(quán)中使得高程誤差變大。此外，MLE解纏結(jié)果中解纏錯(cuò)誤的部分區(qū)域也含有變化的高程信息，解纏出錯(cuò)的區(qū)域只是表現(xiàn)為該區(qū)域的整體抬高或者降低，不能將錯(cuò)點(diǎn)的像素權(quán)值全部設(shè)為0，因而采用比例加權(quán)的方式計(jì)算錯(cuò)點(diǎn)的加權(quán)系數(shù)。

假設(shè)U(p,k)是以p為中心，可調(diào)控值k作為寬度的方形鄰域。p′是在U(p,k)內(nèi)除了p點(diǎn)之外的鄰域像素。

解纏正確的點(diǎn)的均值為

(12)

解纏錯(cuò)誤的點(diǎn)的均值為

(13)

從而計(jì)算出錯(cuò)點(diǎn)的歸一化加權(quán)系數(shù)為

(14)

得到錯(cuò)點(diǎn)修正DEM為

(15)

1)wp′為領(lǐng)域內(nèi)各點(diǎn)的第一個(gè)加權(quán)系數(shù)。在U(p,k)內(nèi)，若被標(biāo)記為解纏錯(cuò)點(diǎn)，則其值為0，若被標(biāo)記為正確的點(diǎn)，其值為1。

2)γN(p′)為鄰域各點(diǎn)的第二個(gè)加權(quán)系數(shù)，將其權(quán)值設(shè)定為U(p,k)鄰域內(nèi)解纏正確的像素的高程均值和解纏出錯(cuò)像素的高程均值之比。

綜合以上兩個(gè)部分，通過(guò)選取合適的模糊類(lèi)閾值，調(diào)節(jié)合適的標(biāo)記窗口和非均勻加權(quán)濾波窗口來(lái)去除噪聲。窗口選擇時(shí)可以利用小窗口方法去掉斑點(diǎn)噪聲，利用大窗口去除噪聲斑塊。而且由于本文方法不是將帶有高程信息的噪斑完全丟掉，如此便可以在保證解纏效率的同時(shí)提高重建精度。

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

實(shí)驗(yàn)分為兩個(gè)部分，分別是建模驗(yàn)證(3.1節(jié))和實(shí)測(cè)數(shù)據(jù)驗(yàn)證(3.2節(jié))。建立的仿真模型為一個(gè)城市地形和一個(gè)帶有懸崖的真實(shí)DEM模型，適合驗(yàn)證雙頻重建算法。仿真模型中的干涉圖均為根據(jù)舍選法[12]生成的干涉圖。實(shí)測(cè)數(shù)據(jù)為機(jī)載C/X雙頻干涉SAR數(shù)據(jù)，比較適合驗(yàn)證雙頻干涉算法。實(shí)驗(yàn)主機(jī)平臺(tái)內(nèi)存為8 GB，CPU主頻為 3.2 GHz。仿真代碼為Matlab代碼。

3.1 仿真數(shù)據(jù)處理

仿真數(shù)據(jù)的參數(shù)如表1所示。本節(jié)均采用該參數(shù)進(jìn)行仿真。第一個(gè)模型最大高程差為250 m，兩個(gè)通道的相鄰像素的相位差分別可達(dá)到34.79π和19.58π，適合用雙頻干涉算法進(jìn)行處理。第二個(gè)模型為美國(guó)Isolation Peak的真實(shí)地形仿真模型。峭壁的高度差可達(dá)137 m。對(duì)應(yīng)的兩個(gè)通道相鄰像素的相位差可達(dá)19.07π和10.73π。對(duì)模擬數(shù)據(jù)的處理結(jié)果用誤差圖和歸一化均方誤差σ2來(lái)進(jìn)行分析。σ2由如下公式給出：

(16)

表1 仿真模型主要參數(shù)

3.1.1 仿真城市地形

建立的仿真地形DEM投影圖如圖2(a)所示。其中三棱柱形高度為250 m。舍選法所得的兩個(gè)通道的干涉圖形如圖2(b)、圖2(c)所示。

(a) 仿真地形投影圖

(b) f=9.6 GHz干涉圖

(c) f=5.4 GHz干涉圖圖2 仿真城市地形數(shù)據(jù)

反演所得MLE結(jié)果如圖3(a)所示，可知MLE結(jié)果中含有大量的噪斑。利用本文算法去噪點(diǎn)所得到的結(jié)果如圖3(b)所示，二次去噪斑結(jié)果如圖3(c)所示。從結(jié)果看出，重建的高程基本將MLE結(jié)果中的噪斑去掉。

(a) MLE重建結(jié)果

(b) 噪點(diǎn)去除結(jié)果

(c) 噪斑去除結(jié)果圖3 仿真城市地形去噪斑

(a) 原始DEM投影

(b) MLE重建DEM投影

(c) MAP重建DEM投影

(d) 本文重建DEM投影圖4 仿真城市地形各種重建算法DEM與原始DEM對(duì)比

圖4分別列出原始DEM、MLE算法DEM、MAP估計(jì)DEM、本文算法DEM來(lái)進(jìn)行對(duì)比。3種重建算法的誤差圖分別如圖5(a)、圖5(b)、圖 5(c)所示。統(tǒng)計(jì)的各種重建方法歸一化均方誤差如表2所示。重建耗時(shí)上：MLE約13.53 s，本文方法約15.61 s；MAP估計(jì)達(dá)到1 493.14 s，是本文算法的將近100倍。從歸一化誤差來(lái)看，MLE算法最差，本文方法和MAP估計(jì)相近。

(a) MLE重建誤差圖

(b) MAP重建誤差圖

(c) 本文方法重建誤差圖圖5 仿真城市地形各種重建算法誤差圖

重建方法歸一化均方誤差MLE0.3268MAP0.0056本文方法0.0064

圖6給出各種算法重建結(jié)果與原始DEM三維對(duì)比。該圖證明本算法在精度上表現(xiàn)較為優(yōu)異。

(a) 原始DEM

(b) MLE重建DEM

(c) MAP重建DEM

(d) 本文算法重建DEM圖6 各重建結(jié)果三維視圖

3.1.2 仿真真實(shí)地形

美國(guó)Isolation Peak國(guó)家公園DEM仿真參數(shù)如表1所示。圖7(a)為該場(chǎng)景的DEM投影圖。兩個(gè)通道的干涉圖如圖7(b)、圖7(c)所示。

圖7 真實(shí)地形數(shù)據(jù)

進(jìn)行高程反演所得的MLE結(jié)果如圖8(a)所示。對(duì)MLE結(jié)果進(jìn)行噪點(diǎn)去除(結(jié)果如圖8(b)所示)和噪斑去除(結(jié)果如圖8(c)所示)。從結(jié)果中看出，重建的高程基本將MLE結(jié)果中的噪斑去掉。

圖9列出了各算法DEM與原始DEM。3種重建算法的誤差圖分別如圖10(a)、圖10(b)、圖10(c)所示。統(tǒng)計(jì)的誤差如表3所示。重建時(shí)間來(lái)看：MLE算法為2.14 s，本文方法為3.08 s，MAP估計(jì)法為258.14 s。說(shuō)明本文方法可以保證重建精度和重建效率。

圖8 真實(shí)地形去噪

圖9 真實(shí)地形各重建結(jié)果與原始DEM對(duì)比

圖10 仿真城市地形各種重建算法誤差圖

圖11為各重建結(jié)果三維視圖。算法優(yōu)劣此處不再贅述。

表3 真實(shí)地形重建算法歸一化均方誤差表

(a) 原始DEM

(b) MLE算法DEM

(c) MAP算法DEM

(d) 本文方法DEM圖11 各重建結(jié)果三維視圖

3.2 實(shí)測(cè)數(shù)據(jù)處理

實(shí)測(cè)數(shù)據(jù)采用了機(jī)載C/X雙頻干涉SAR測(cè)量數(shù)據(jù)，該機(jī)載平臺(tái)的具體參數(shù)如表4所示。第一塊實(shí)驗(yàn)數(shù)據(jù)為澄城縣城的干涉SAR數(shù)據(jù)。縣城中高樓高度有90～100 m。兩個(gè)通道相鄰像素之間的干涉相位差分別可達(dá)3.2π和5.8π。第二塊數(shù)據(jù)為白水縣東坡村附近數(shù)據(jù)，該塊數(shù)據(jù)中有溝壑和高架橋，地形復(fù)雜多變。兩塊數(shù)據(jù)均適合驗(yàn)證雙頻干涉處理算法。

表4 機(jī)載C/X雙頻系統(tǒng)主要參數(shù)

3.2.1 澄城縣實(shí)測(cè)數(shù)據(jù)處理

所選取的地域和對(duì)應(yīng)的光學(xué)圖像如圖12所示。圖中黃色旗幟所示的點(diǎn)為選取的控制點(diǎn)。對(duì)澄城縣的兩幅C/X波段圖像分別進(jìn)行配準(zhǔn)，得到兩幅干涉條紋圖分別如圖13(a)、圖13(b)所示，利用軌道法[13]去除平地效應(yīng)。再將C/X波段得到的干涉圖配準(zhǔn)。得到配準(zhǔn)的去平地干涉相位圖如圖13(c)、圖13(d)所示。MLE算法進(jìn)行高程重建結(jié)果如圖14(a)所示，利用本文方法去噪所得結(jié)果如圖14(b)所示。圖15(a)、圖15(b)所示為重建結(jié)果的三維視圖，可以更清楚地看出本文方法的有效性。

圖12 澄城縣Google Earth光學(xué)圖像與干涉SAR數(shù)據(jù)圖像

(a) C波段干涉相位圖

(b) X波段干涉相位圖

(c) C波段去平地干涉相位圖

(d) X波段去平地干涉相位圖圖13 各個(gè)通道去平地前后干涉相位圖

(a) MLE重建結(jié)果

(b) 本文方法結(jié)果圖14 MLE算法重建結(jié)果與本文方法去噪結(jié)果

對(duì)圖12中所選取的控制點(diǎn)進(jìn)行高程統(tǒng)計(jì)，統(tǒng)計(jì)結(jié)果如表5所示。從結(jié)果中可以看出，MLE算法噪聲過(guò)多，本文方法反演高度接近于真實(shí)高度。由于使用MAP耗時(shí)過(guò)長(zhǎng)，這里沒(méi)有采用。

(a) MLE結(jié)果三維DEM

(b) 本文結(jié)果三維DEM圖15 MLE算法重建結(jié)果與本文方法去噪結(jié)果三維視圖

m

3.2.2 白水縣東坡村實(shí)測(cè)數(shù)據(jù)處理

首先給出的是光學(xué)圖像和SAR圖像如圖16所示。黃色旗幟為誤差分析控制點(diǎn)。其中1號(hào)點(diǎn)地區(qū)為溝壑，2號(hào)點(diǎn)白色條狀地形為高架橋的影像，3號(hào)點(diǎn)地區(qū)為平地上的河流。

圖16 白水縣東坡村附近高架Google Earth光學(xué)圖像與干涉SAR數(shù)據(jù)圖像

兩幅干涉條紋圖分別如圖17(a)、圖17(b)所示，配準(zhǔn)的去平地干涉相位圖如圖17(c)、圖17(d)所示。

(a) C波段干涉條紋

(b) X波段干涉條紋

(c) C波段去平地干涉相位

(d) X波段去平地干涉相位圖17 各個(gè)通道去平地前后干涉相位圖

MLE算法重建結(jié)果如圖18(a)所示，利用本文方法進(jìn)行重建結(jié)果如圖18(b)所示。可見(jiàn)噪斑去除效果明顯，具有良好的重建性能。

(a) MLE重建結(jié)果

(b) 本文方法重建結(jié)果圖18 MLE算法重建結(jié)果與本文方法去噪結(jié)果

控制點(diǎn)高程統(tǒng)計(jì)結(jié)果如表6所示。在某些地帶MLE算法噪聲過(guò)大，例如本例中的3號(hào)控制點(diǎn)河流部分，高架邊緣部分。統(tǒng)計(jì)結(jié)果表明，本文方法有效提高了重建精度。

表6 控制點(diǎn)高度統(tǒng)計(jì)表 m

圖19(a)、圖19(b)列出了MLE和本文方法的重建結(jié)果三維視圖。為了減少內(nèi)存消耗，在生成三維圖形時(shí)，在方位向作了降采樣處理，因而三維視圖中方位向點(diǎn)數(shù)變少。

(a) MLE重建結(jié)果

(b) 本文方法重建結(jié)果圖19 MLE算法重建結(jié)果與本文方法去噪結(jié)果三維視圖

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

本文對(duì)傳統(tǒng)的基于MLE的雙頻重建算法進(jìn)行了改進(jìn)，可以去除MLE重建結(jié)果中的噪點(diǎn)和噪聲斑塊，在保證重建效率的同時(shí)提高重建精度?；谏衔哪M模型的論證分析，本文的方法在精度上與MAP估計(jì)相近，但是在重建效率上要遠(yuǎn)遠(yuǎn)優(yōu)于MAP估計(jì)。上文中的實(shí)測(cè)數(shù)據(jù)處理結(jié)果表明，本文提出的方法可以用在大規(guī)模數(shù)據(jù)處理過(guò)程中，擁有較為良好的處理效果。但是本文方法也存在不足，例如閾值的選取和加權(quán)窗口的選取一般都是選取經(jīng)驗(yàn)值，后期需要添加自適應(yīng)的閾值控制方法來(lái)降低后處理的難度。此外，更加合理的去噪權(quán)值應(yīng)該為越臨近所求像素點(diǎn)的權(quán)值越大，越遠(yuǎn)的像素點(diǎn)權(quán)值越小，這也是以后工作中一個(gè)要調(diào)整的方面。