杜順季，邢 誠,2

（1.武漢大學(xué) 測繪學(xué)院，湖北 武漢 430079；2.精密工程與工業(yè)測量國家地理信息局重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079）

合成孔徑雷達(dá)作為主動式微波遙感技術(shù)，因不受時間和天氣條件的限制，成為對地觀測系統(tǒng)重要組成部分。1993年Massonnet等人在《Nature》上發(fā)表了利用ERS-1獲得與GPS一致的Landers地震形變場的文章[1]，由此DInSAR技術(shù)得到了廣泛研究和應(yīng)用。

目前能夠系統(tǒng)處理SAR數(shù)據(jù)軟件的主要是商用Gamma、ENVI的SARscape模 塊 和doris、ROI_PAC等開源軟件，商用軟件價格不菲而開源軟件一般基于Linux平臺的命令模式，不方便用戶操作和交互。本文利用以歐空局為主研發(fā)的跨平臺開源軟件NEST來進(jìn)行SAR圖像處理的研究，以伊朗Bam地震為例，利用ASAR影像基于兩軌法DInSAR技術(shù)獲取視線向的地表形變量，介紹了NEST軟件處理SAR數(shù)據(jù)的過程，得到了較好的結(jié)果。

1 二軌法DInSAR基本原理

根據(jù)去地形影響所采用的數(shù)據(jù)和方法，DInSAR可分為二軌法、三軌法和四軌法[2]。二軌法是利用外部高精度地面數(shù)字高程模型(DEM)和軌道數(shù)據(jù)來模擬地形相位，從原始干涉圖中減去該模擬的地形相位就得到新的差分干涉相位圖，去除了地形影響而僅包含了形變相位。

利用二軌法進(jìn)行形變相位的提取處理比較簡單，其優(yōu)點(diǎn)是無需對干涉圖相位解纏，避免了解纏困難[3]。在沒有外部DEM的情況下，可選擇三軌法和四軌法。但是利用三軌法和四軌法的地形像對獲取的DEM會進(jìn)一步引入大氣效應(yīng)的影響，從而引起形變測量的誤差。二軌法DEM高程誤差會直接影響形變相位，根據(jù)SAR測量幾何圖形和誤差傳播率可知：

式中，σΔr為形變測量誤差；σh為DEM高程誤差；B⊥為垂直基線分量；R為地面點(diǎn)至天線距離。由式(1)可知，當(dāng)垂直基線為0時，DEM高程誤差對形變測量結(jié)果無影響，這在地基SAR中是重要特性[4]，可以直接利用2幅圖像獲得地表形變，不需要考慮測量區(qū)域的地形相位，但是在星載SAR中基線不可能為零。目前美國航空局公開了全球的SRTM數(shù)字高程數(shù)據(jù)，高程精度可達(dá)16 m[5]，在平原地區(qū)絕對高程精度優(yōu)于5 m，基本上滿足了差分干涉測量的要求，為二軌法DInSAR的廣泛應(yīng)用提供了數(shù)據(jù)[6]。

2 NEST數(shù)據(jù)處理實(shí)踐

2.1 NEST軟件主要功能及特點(diǎn)

最新發(fā)布的NEST版本為4C-1.1，能夠?qū)AR數(shù)據(jù)進(jìn)行讀取、后處理、分析和可視化，支持ERS1/2、Envisat的ESA衛(wèi)星數(shù)據(jù)和諸如JERS、ALOS、TerreSAR-X、Radarsat1/2和Cosmo-Skymed等其他數(shù)據(jù)。主要功能包括有效的數(shù)據(jù)查看和管理、絕對定標(biāo)、軌道與濾波處理、自動配準(zhǔn)、干涉、地形改正、統(tǒng)計(jì)分析、海洋工具和全面集成doris。其特點(diǎn)為：良好的操作界面和靈活的參數(shù)設(shè)定、數(shù)據(jù)可視化和管理、自動下載軌道和DEM數(shù)據(jù)、批處理功能、數(shù)據(jù)統(tǒng)計(jì)和分析工具、多種數(shù)據(jù)格式輸出等。該軟件集成了doris軟件算法，既保證了良好的圖形化用戶界面，又保留了高效的批處理框架，實(shí)現(xiàn)了SAR處理的各方面功能并提供眾多分析工具，為SAR研究帶來極大的方便[7]。

2.2 像對數(shù)據(jù)選取與處理

歐空局在Bam地震后公開了震區(qū)的Envisat ASAR數(shù)據(jù)，為VV極化，地震時間為2003-12-26，這次使用的2景降軌的地震前后圖像，其相關(guān)參數(shù)見表1。

出該像對的垂直基線只有0.6 m，根據(jù)式(1)，用該數(shù)據(jù)計(jì)算形變場，地形高程引起的誤差影響值很小，能夠達(dá)到很高的形變監(jiān)測精度。DInSAR數(shù)據(jù)處理的主要步驟包括：影像配準(zhǔn)、干涉圖生成、DEM模擬相位、干涉差分、干涉圖相位濾波、相位解纏、形變量計(jì)算和地理編碼。NEST軟件實(shí)現(xiàn)了以上步驟并在每個操作時給出了詳細(xì)的參數(shù)控制選項(xiàng)，能根據(jù)要求進(jìn)行靈活的影像處理，相關(guān)的濾波和去噪處理可以選擇性進(jìn)行，其大致流程如圖1所示。

圖1 NEST處理流程圖

數(shù)據(jù)輸入階段，NEST軟件可通過Apply orbit file操作更新元數(shù)據(jù)中的軌道狀態(tài)矢量。對于ASAR數(shù)據(jù)提供doris或delft精密軌道文件的自動下載，這些文件提供每60 s的衛(wèi)星位置和速度。預(yù)處理中輻射校準(zhǔn)是為了得到反映地物后向散射系數(shù)的圖像，是關(guān)聯(lián)像素值和散射系數(shù)的重要工作[8]，對于SAR數(shù)據(jù)的定量分析有重要意義。

在InSAR處理中，圖像對的精確配準(zhǔn)是至關(guān)重要的，保證2幅復(fù)圖像的點(diǎn)對應(yīng)地面同一個點(diǎn)，一般達(dá)到1/8像元的亞像元級配準(zhǔn)精度才符合處理要求。NEST通過控制點(diǎn)選取操作進(jìn)行精配準(zhǔn)設(shè)置，可優(yōu)化插值函數(shù)、配準(zhǔn)窗口大小、多項(xiàng)式次數(shù)等參數(shù)，提供了最鄰近、線性、4點(diǎn)/6點(diǎn)立方卷積、sinc插值等眾多插值函數(shù)。本文采用精配準(zhǔn)，二次多項(xiàng)式和6點(diǎn)立方卷積插值得到的配準(zhǔn)殘差均值Mean=0.02，標(biāo)準(zhǔn)差Std=0.02，達(dá)到了很高的配準(zhǔn)精度，體現(xiàn)了該軟件算法的準(zhǔn)確性。

配準(zhǔn)之后可對圖像進(jìn)行斑點(diǎn)噪聲濾波，以便得到更好的干涉圖。NEST提供2種濾波方式：針對單幅圖像的經(jīng)典濾波(包括中值濾波、均值濾波、Frost濾波、Lee濾波和Gamma濾波等)和多時相濾波，分別采用5×5的Gamma濾波和多時相濾波進(jìn)行處理得到的圖像如圖2、圖3所示。可以看出，濾波總體上有平滑的效果，特別是單幅圖像的空間濾波使得圖像更平滑但降低了分辨率，多時相濾波保留了一定的邊緣信息，圖像分辨率沒有下降很多。多時相濾波引入時間維度，能夠綜合多幅圖像中的信息，將相同位置不同時相的像元以某種方式進(jìn)行組合，達(dá)到更好的斑點(diǎn)噪聲抑制效果[9]。

圖2 單幅圖像濾波

圖3 多時相濾波

經(jīng)過配準(zhǔn)的SAR圖像共軛相乘得到干涉圖，再消除平地相位，該相位的存在使得條紋密集，影響后續(xù)處理和分析，必須消除。NEST軟件在同一個操作中實(shí)現(xiàn)，基于軌道參數(shù)和成像區(qū)域信息，通過分布在整幅影像的控制點(diǎn)進(jìn)行最小二乘擬合計(jì)算平地相位，在干涉圖中減去該部分相位。同時SAR成像的斜視幾何效應(yīng)，產(chǎn)生透視收縮、陰影、頂?shù)椎怪玫葓D像畸變，二軌差分中NEST自動下載相應(yīng)地區(qū)的DEM數(shù)據(jù)，依據(jù)軌道和影像數(shù)據(jù)進(jìn)行地形改正，同時可完成地理編碼，支持ENVI、HDF5、GeoTIFF等格式輸出，特別是可輸出為Google Earth支持的KMZ格式，方便用戶直觀分析結(jié)果。以上所有操作也可以在Graph Builder的命令行可視化工具中完成批處理，該窗口分為了2部分，上半部分是所執(zhí)行命令過程圖，下半部分是每條命令的詳細(xì)參數(shù)設(shè)置，也是該軟件強(qiáng)大的特色工具之一。

經(jīng)過這一系列處理可得到二軌差分結(jié)果見圖4，每條干涉條紋表示28 mm形變量，南面10個條紋表示衛(wèi)星視線向形變約28 cm，北面7個干涉條紋，形變約19.6 cm，最大的相對位移變化大約為48 cm，與德國地學(xué)研究中心(GFZ)等機(jī)構(gòu)公布的結(jié)果一致[10]，深入分析見文獻(xiàn)[11]，輸出到Google Earth疊加結(jié)果如圖5所示。

圖4 Bam形變干涉相位圖

圖5 地理編碼輸出到Google Earth結(jié)果圖

3 結(jié) 語

本文介紹了二軌法差分干涉的主要步驟，并對Bam震區(qū)的SAR數(shù)據(jù)進(jìn)行了處理。結(jié)果表明，DInSAR技術(shù)能有效獲取地表變形，且NEST軟件可以有效地進(jìn)行SAR數(shù)據(jù)處理，良好的可視化操作界面使其使用更加方便，自動下載軌道和DEM數(shù)據(jù)并可對每個步驟詳細(xì)設(shè)置，支持批處理，更加便捷和靈活，讓SAR數(shù)據(jù)處理更加簡單，有利于相關(guān)問題的學(xué)習(xí)和研究。

[1]Massonnet D, Rossi M, Carmona C, et al. The Displacement Field of the Landers Earthquake Mapped by Radar Interferometry[J]. Nature,1993, 364(6433): 138-142

[2]廖明生，林琿. 雷達(dá)干涉測量: 原理與信號處理基礎(chǔ)[M]. 北京: 測繪出版社, 2003

[3]張紅，王超，吳濤，等. 基于相干目標(biāo)的DInSAR方法研究[M]. 北京: 科學(xué)出版社, 2009

[4]Casagli N, Catani F, Del Ventisette C, et al. Monitoring,Prediction, and Early Warning Using Ground-based Radar Interferometry[J]. Landslides, 2010, 7(3): 291-301

[5]陳俊勇. 對SRTM3和GTOPO30地形數(shù)據(jù)質(zhì)量的評估[J].武漢大學(xué)學(xué)報:信息科學(xué)版,2005, 30(11): 941-944

[6]何敏，陸曉燕，何秀鳳. 利用DInSAR二軌法監(jiān)測徐州大屯中心區(qū)地表形變[J]. 地理空間信息,2011,9(5): 3-5

[7]ESA. NEST User Manual .[EB/OL].http://nest.array.ca/web/nest/documentation,2012-02-21

[8]袁禮海，葛家龍，江凱,等. SAR輻射定標(biāo)精度設(shè)計(jì)與分析[J]. 雷達(dá)科學(xué)與技術(shù),2009, 7(1): 35-39

[9]龔麗霞，吳樊，曾琪明,等. SAR 圖像多時相組合濾波方法[J].計(jì)算機(jī)工程,2011, 37(21):3-5

[10]Xia Ye. Bam Earthquake: D-INSAR Preliminary Result.[EB/OL]. http://www.gfz-potsdam.de/portal/gfz/Struktur/Departments/Department+1/sec14/topics/radar/sar/earthquake_bam,2012-11-23

[11]Fialko Y, Sandwell D, Simons M, et al. Three-dimensional Deformation Caused by the Bam, Iran, Earthquake and the Origin of Shallow Slip Deficit[J]. Nature,2005, 435(7040): 295-299