高雨欣，郭 睿，王 輝，張雙喜

（1.西北工業(yè)大學(xué) 電子信息學(xué)院，陜西 西安 710129；2.西北工業(yè)大學(xué) 自動(dòng)化學(xué)院，陜西 西安 710129；3.上海衛(wèi)星工程研究所，上海 201109）

0 引言

干涉合成孔徑雷達(dá)（Interferometric Synthetic Aperture Radar，InSAR）通過對(duì)同一地區(qū)獲取的兩幅或多幅相干圖像進(jìn)行干涉處理獲取高精度的地面高程信息，在地形測(cè)繪、地表變形監(jiān)測(cè)、土地覆蓋分類、目標(biāo)監(jiān)測(cè)、寬幅場(chǎng)景測(cè)速等遙感領(lǐng)域發(fā)揮著重要作用。在InSAR 數(shù)據(jù)處理過程中，干涉相位濾波是減少由熱噪聲、不同去相關(guān)和SAR圖像配準(zhǔn)引入噪聲的基本步驟。因此，先進(jìn)的相位濾波方法成為了干涉處理中獲得高質(zhì)量干涉圖的必須研究課題之一。

非局部濾波器可以有效地捕捉相位結(jié)構(gòu)，尤其在條紋濾波和細(xì)節(jié)保持方面表現(xiàn)優(yōu)異，其在InSAR相位去噪中的應(yīng)用引起了學(xué)者們極大的關(guān)注。DELEDALLE 等提出了一種基于幅度、相干系數(shù)和干涉相位聯(lián)合統(tǒng)計(jì)分布的非局部干涉相位估計(jì)方法。SICA 等將基于非局部方法的塊匹配3-D（The Block-Matching 3-D，BM3D）算法擴(kuò)展到InSAR 相位恢復(fù)領(lǐng)域。為進(jìn)一步改進(jìn)非局部均值濾波器的濾波結(jié)果，Stein 的無偏風(fēng)險(xiǎn)估計(jì)（Stein’s Unbiased Risk Estimate，SURE）準(zhǔn)則被引入非局部均值濾波。YANG 等提出了一種基于統(tǒng)計(jì)特性非局部均值和魯棒雙邊濾波器的混合算法，該算法使用SURE準(zhǔn)則估計(jì)最優(yōu)參數(shù)。LI 等提出了一種基于變分貝葉斯推理和SURE 準(zhǔn)則的聯(lián)合非局部去噪方法。目前關(guān)于SURE 準(zhǔn)則下非局部濾波在實(shí)測(cè)InSAR 干涉相位數(shù)據(jù)的噪聲抑制中的有效性研究較少。

基于此，提出了一種基于SURE 準(zhǔn)則的非局部均值濾波器自適應(yīng)參數(shù)選擇方法，并應(yīng)用于大場(chǎng)景下的InSAR 實(shí)測(cè)數(shù)據(jù)。首先，簡(jiǎn)要回顧了非局部均值算法和SURE 準(zhǔn)則在非局部均值濾波器上的應(yīng)用；然后，著重介紹了應(yīng)用于InSAR 實(shí)測(cè)數(shù)據(jù)相位濾波的SURE 準(zhǔn)則參數(shù)自適應(yīng)非局部均值濾波器；最后，通過仿真數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)評(píng)估了所提出方法的有效性，并將結(jié)果與其他2 種廣泛使用的濾波器進(jìn)行了比較。

1 理論背景

1.1 非局部均值算法

含噪圖像可以描述為

式中：（）、（）分別為實(shí)測(cè)圖像和原始圖像中的第個(gè)像素；為方差為σ的高斯噪聲。

非局部均值濾波器對(duì)（）的估計(jì)為

式中：（）為以為中心的搜索窗；（，）為由像素與周圍的鄰域均值計(jì)算得到的權(quán)值。

（，）按照如下定義：

式中：為非局部濾波參數(shù)，取決于噪聲水平，與噪聲標(biāo)準(zhǔn)差成正比；是鄰域的像素?cái)?shù)，例如鄰域大小為7×7 時(shí)，=［-3，3］×［-3，3］，=49。

1.2 SURE 在非局部均值濾波中的應(yīng)用

去噪圖像的均方誤差（Mean Square Error，MSE）定義式為

式中：‖·‖為歐式范數(shù)；為去噪圖像中的像素?cái)?shù)。

SURE 準(zhǔn)則提供了一種對(duì)真MSE 進(jìn)行無偏估計(jì)的方法，第個(gè)像素的SURE 值表示為

據(jù)式（2）和式（3）、式（5）中偏導(dǎo)數(shù)表示為

通過計(jì)算得到與每個(gè)像素的濾波參數(shù)對(duì)應(yīng)的SURE 值，最小的SURE 值對(duì)應(yīng)的濾波參數(shù)最佳。

2 用于INSAR 的SURE 準(zhǔn)則非局部濾波

2.1 InSAR 相位噪聲模型

InSAR 相位噪聲為如下加性噪聲模型：

式中：?為實(shí)測(cè) 相位；?為無噪聲InSAR 相位；為零均值加性噪聲。

由于InSAR 相位在［-π，π］的區(qū)間內(nèi)存在相位跳變，無法直接應(yīng)用實(shí)域中的干涉相位進(jìn)行濾波，而為了正確提取目標(biāo)信息，發(fā)生在纏繞相位條紋中的相位跳變應(yīng)當(dāng)被保留。然而，由于復(fù)數(shù)域中的InSAR 相位是連續(xù)的，且使用復(fù)數(shù)據(jù)進(jìn)行干涉相位濾波的效果更佳，故采用復(fù)數(shù)域中的InSAR 相位模型：

式中：S、S分別為復(fù)數(shù)域中的實(shí)測(cè)相位和無噪聲相位；n為復(fù)數(shù)域中的相位噪聲。

S的實(shí)部和虛部可以分別表示為

式中：、分別為n的實(shí)部和虛部，可以將實(shí)部和虛部的噪聲分別視為零均值加性噪聲；是由相干性確定的相位質(zhì)量指標(biāo)。

2.2 提出的方法

所提出方法根據(jù)2.1 節(jié)中InSAR 相位噪聲模型，在實(shí)部和虛部中分別對(duì)InSAR 相位進(jìn)行SURE準(zhǔn)則快速非局部均值濾波，從而更好地提取目標(biāo)信息，流程如圖1 所示。

圖1 InSAR 相位濾波流程Fig.1 Flow chart of InSAR phase filtering

該方法側(cè)重于自適應(yīng)的濾波參數(shù)的選擇，使其適用于InSAR 實(shí)測(cè)復(fù)雜場(chǎng)景，尤其是城市區(qū)域的相位濾波。對(duì)于實(shí)部和虛部，選擇自適應(yīng)參數(shù)的步驟類似，該方法在InSAR 相位濾波中的應(yīng)用步驟如下：

噪聲方差估計(jì)：由于濾波參數(shù)很大程度上取決于實(shí)測(cè)數(shù)據(jù)噪聲方差σ，將實(shí)部和虛部分離后，首先采用零均值加性噪聲方差的估計(jì)器以獲得噪聲方差σ。

濾波參數(shù)設(shè)置：根據(jù)不同的數(shù)據(jù)集來設(shè)置濾波器的參數(shù)，如鄰域窗的半徑、搜索窗的半徑以及濾波參數(shù)。不同于傳統(tǒng)非局部均值濾波器采用固定的濾波參數(shù)，所提出方法通過比較不同數(shù)據(jù)集下不同范圍的濾波結(jié)果，給出最佳濾波參數(shù)范圍為=［0.55∶0.005∶0.6］×σ。

SURE 值計(jì)算：SURE 值的計(jì)算是該非局部均值濾波算法的核心。依據(jù)步驟2 中給出的濾波參數(shù)范圍執(zhí)行快速非局部均值濾波，同時(shí)計(jì)算每個(gè)像素的SURE 值。對(duì)于給定范圍中的每個(gè)濾波參數(shù)都可以獲得一個(gè)大小與原數(shù)據(jù)相同SURE值數(shù)組。

最佳濾波參數(shù)選擇：在完成對(duì)給定濾波參數(shù)范圍內(nèi)所有的非局部均值濾波之后，為輸入數(shù)據(jù)中的每個(gè)像素選擇最小SURE 值，并儲(chǔ)存對(duì)應(yīng)的最佳濾波參數(shù)。比起傳統(tǒng)非局部均值濾波器對(duì)所有像素統(tǒng)一選擇固定的濾波參數(shù)，所提出的方法依據(jù)SURE 準(zhǔn)則對(duì)每個(gè)像素選擇相應(yīng)的最佳濾波參數(shù)，對(duì)復(fù)雜場(chǎng)景中的不同區(qū)域具有更好的適應(yīng)能力。

采用步驟2 中的濾波器窗口半徑和以及獲取的自適應(yīng)濾波參數(shù)數(shù)組，對(duì)該部分（實(shí)部或虛部）進(jìn)行快速非局部均值濾波。最后，將實(shí)部和虛部的濾波結(jié)果組合，獲得濾波后的干涉相位，至此完成干涉相位濾波的所有步驟。

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

采用仿真數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)對(duì)提出方法進(jìn)行實(shí)驗(yàn)驗(yàn)證，與經(jīng)典Goldstein 濾波器和非自適應(yīng)的原始非局部濾波器的結(jié)果進(jìn)行對(duì)比分析。其中，Goldstein 濾波器的窗口大小均設(shè)置為13×13，濾波參數(shù)設(shè)置為0.5。非局部均值濾波器和所提出方法的鄰域窗均設(shè)置為10×10，搜索窗大小設(shè)置為23×23。所有實(shí)驗(yàn)均在配置為2.60 GHz Intel Core CPU和16 GB 內(nèi)存的PC 上進(jìn)行。

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

仿真數(shù)據(jù)為Matlab 函數(shù)隨機(jī)生成的常用仿真干涉條紋，并加入信噪比為1.406 dB 的高斯白噪聲，如圖2（a）所示。采用經(jīng)典Goldstein，原始非局部濾波和所提出的濾波方法進(jìn)行噪聲抑制后的干涉條紋圖如圖2（b）～圖2（d）所示。

圖2 仿真數(shù)據(jù)的濾波結(jié)果Fig.2 Filtering results of the simulation data

由上述可知，原始非局部均值和所提出方法的濾波結(jié)果比經(jīng)典Goldstein 濾波結(jié)果更平滑。在白框標(biāo)出的條紋密集區(qū)域，與圖2（b）與圖2（c）相比，圖2（d）中條紋連續(xù)性的保持效果最佳，顯示出所提出方法在條紋圖案細(xì)節(jié)保持方面更優(yōu)秀的性能。

此外，通過殘差數(shù)（Number of Residues，NOR）和峰值信噪比（Peak Signal to Noise Ratio，PSNR）對(duì)濾波結(jié)果進(jìn)行比較評(píng)估，同時(shí)給出每種方法的運(yùn)行時(shí)間。仿真數(shù)據(jù)濾波結(jié)果的定量指標(biāo)和運(yùn)行時(shí)間見表1，較小的NOR 和較高的PSNR 表明濾波器的降噪效果更好。在3 種濾波結(jié)果中，所提出的方法獲得了最大的PSNR。經(jīng)典Goldstein濾波減少了98.87%的含噪干涉圖的NOR，非局部均值濾波器減少了98.28%，提出的方法減少了98.50%，所有方法在噪聲抑制方面的效果差別不大。

表1 仿真數(shù)據(jù)濾波評(píng)估指標(biāo)對(duì)比Tab.1 Comparison of the estimated indexes for simulation data filtering

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

采用TanDEM-X 獲取的拉斯維加斯及其周邊區(qū)域的干涉實(shí)測(cè)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)驗(yàn)證。該城市場(chǎng)景包含各種類型的建筑，便于評(píng)估不同濾波方法在城市區(qū)域?yàn)V波的效果。含噪干涉相位如圖3（a）所示，經(jīng)典Goldstein 濾波器、原始非局部均值濾波器和所提出方法的濾波結(jié)果分別如圖3（b）～圖3（d）所示。

圖3 實(shí)測(cè)城市數(shù)據(jù)的濾波結(jié)果Fig.3 Filtering results of the measured urban data

圖3（c）中非局部均值濾波結(jié)果和圖3（d）中提出方法的濾波結(jié)果可以更好地抑制相位噪聲。然而，原始非局部均值濾波器的結(jié)果顯示出干涉相位被過度濾波，造成更多的細(xì)節(jié)損失。此外，經(jīng)典Goldstein 濾波器、原始非局部均值濾波器和所提出的方法的運(yùn)行時(shí)間分別為272.406 7、4 674.544 3、3 706.487 7 s。

將圖3（a）中標(biāo)出的均質(zhì)區(qū)域和高層建筑區(qū)域放大，以便更清楚地進(jìn)行細(xì)節(jié)比較，如圖4 和圖5所示。

圖4 實(shí)測(cè)城市數(shù)據(jù)均質(zhì)區(qū)域?yàn)V波結(jié)果Fig.4 Filtering results of the homogenous areas of the measured urban data

圖5 實(shí)測(cè)城市高層建筑區(qū)域Fig.5 High-rise building areas of the measured urban data

對(duì)于均質(zhì)區(qū)域，Goldstein 濾波器的結(jié)果明顯存在更多的斑點(diǎn)噪聲，且表2 中雖然Goldstein 濾波器獲得最高的PSNR，但其只能降低62.34%的NOR。在圖4（c）和圖4（d）中，原始非局部均值濾波器和所提出的方法在均質(zhì)區(qū)域中均呈現(xiàn)出較好的噪聲抑制效果。原始非局部均值濾波器減少了99.88%的NOR，而提出的方法減少了97.80%的NOR，且取得了更高的PSNR。這2 種方法在均質(zhì)區(qū)域的噪聲抑制方面具有相似的性能。

表2 均質(zhì)區(qū)域?yàn)V波評(píng)估指標(biāo)對(duì)比Tab.2 Comparison of the estimated indexes for the homogeneous areas

由于高層建筑的干涉條紋、邊緣在城市建筑重建和高度信息的提取中起關(guān)鍵作用，在相位濾波中應(yīng)著重保留這些基本細(xì)節(jié)。所提出方法的另一個(gè)優(yōu)點(diǎn)在于能夠更好地保持建筑區(qū)域的邊緣信息。如圖5 所示，與經(jīng)典Goldstein 濾波器和非局部均值濾波器相比，提出方法在圖5（d）中更好地保留了建筑邊緣形狀和條紋細(xì)節(jié)。

邊緣保持度（Edge Preservation Degree of Ratio of Average，EPD-ROA）可用于驗(yàn)證干涉圖中建筑邊緣保持效果。指標(biāo)更接近1 時(shí)，相應(yīng)的方法在邊緣保持方面的效果更好。3 種濾波方法在建筑區(qū)域的邊緣保持度見表3，所提出方法的EPD-ROA比2 種傳統(tǒng)方法更接近于1，表明該方法在水平和垂直方向上均取得了最好的邊緣保持效果。

表3 建筑區(qū)域EPD-ROA 對(duì)比Tab.3 Comparison of the EPD-ROAs in the building areas

為進(jìn)一步評(píng)價(jià)改進(jìn)方法的整體濾波效果，采用最小二乘法對(duì)上述3 種濾波器濾波后的干涉相位解纏，以比較相位濾波結(jié)果對(duì)后續(xù)相位解纏處理的影響。解纏前的噪聲干涉相位如圖6（a）所示，采用3 種濾波方法進(jìn)行濾波后的干涉相位解纏結(jié)果分別如圖6（b）～圖6（d）所示。

總體而言，Goldstein 濾波后解纏結(jié)果和所提方法濾波后的結(jié)果相似。而原始非局部均值濾波后的解纏結(jié)果對(duì)干涉圖的相位保持效果較差。比較圖6（b）和圖6（d）中標(biāo)出的建筑區(qū)域，提出的方法比Goldstein 濾波器更有效地保持建筑結(jié)構(gòu)的形狀，建筑邊緣更為清晰。此外，所提出的方法在均質(zhì)區(qū)的解纏結(jié)果比Goldstein 濾波器更平滑。

圖6 實(shí)測(cè)城市數(shù)據(jù)集的解纏結(jié)果Fig.6 Results of the unwrapped phase of the measured urban data

綜合考慮濾波后的纏繞相位結(jié)果和解纏相位的結(jié)果，所提出方法在均質(zhì)區(qū)域與城市建筑區(qū)域均能進(jìn)行有效的噪聲抑制和細(xì)節(jié)保持，適用于包含多種類型建筑的復(fù)雜城市場(chǎng)景。

4 結(jié)束語

針對(duì)傳統(tǒng)InSAR 相位濾波方法不能實(shí)現(xiàn)自適應(yīng)，且在噪聲抑制和細(xì)節(jié)保持能力上不能兼得的問題，介紹了一種基于SURE 準(zhǔn)則獲取自適應(yīng)濾波參數(shù)的改進(jìn)非局部均值濾波器，并將其應(yīng)用于InSAR實(shí)測(cè)數(shù)據(jù)的干涉相位濾波。提出的方法利用InSAR 相位復(fù)數(shù)域模型，將其分為實(shí)部和虛部分別濾波，并依據(jù)SURE 準(zhǔn)則估計(jì)最優(yōu)濾波參數(shù)實(shí)現(xiàn)自適應(yīng)，最終合并實(shí)部和虛部的結(jié)果得到濾波后的InSAR 相位。仿真數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)的實(shí)驗(yàn)結(jié)果均表明，該方法在時(shí)間消耗和細(xì)節(jié)保持方面優(yōu)于非自適應(yīng)的InSAR 非局部均值濾波器，尤其是在條紋圖案的連續(xù)性和邊緣信息的保留方面更優(yōu)。