用InSAR研究地面形變場DEM的提取

余景波，曹振坦，薛峰，趙玉江

(1.山東科技大學測繪科學與工程學院，山東青島 266510； 2.東方道邇濟南分公司，山東濟南 250101；3.山東省城鄉(xiāng)建設(shè)勘察院地信公司，山東濟南 250031)

用InSAR研究地面形變場DEM的提取

余景波1?，曹振坦1，薛峰2，趙玉江3

(1.山東科技大學測繪科學與工程學院，山東青島 266510； 2.東方道邇濟南分公司，山東濟南 250101；3.山東省城鄉(xiāng)建設(shè)勘察院地信公司，山東濟南 250031)

介紹了InSAR提取DEM的基本原理，以Earhview InSAR的兩種生成DEM數(shù)據(jù)處理模式處理ENVISAT衛(wèi)星獲取的巴姆地震影像數(shù)據(jù)，結(jié)果分別提取了巴姆地震同震形變場的DEM。對結(jié)果進行了對比分析，表明借助外部DEM可以提高提取的地面形變場DEM的精度。

Earthview InSAR；DEM生成模式；巴姆地震；ENVISAT衛(wèi)星；DEM的精度

1 引 言

合成孔徑雷達干涉測量(InSAR)在20世紀60年代末出現(xiàn)，是將雷達影像復數(shù)據(jù)的相位信息作為信息源提取地表三維信息的一項技術(shù)。1974年，Graham等人[1]第一次提出用InSAR技術(shù)進行制圖的設(shè)想，一直到1986年，通過Goldstein和Zebker[2]的改進，才真正實現(xiàn)數(shù)字化干涉雷達技術(shù)，并利用其生成數(shù)字地形圖。1988年，Goldstain等人[3]采用InSAR技術(shù)處理SEASAT星載SAR數(shù)據(jù)，獲取了死亡谷Cottonball盆地的地形圖，其地形數(shù)據(jù)與出版的USGS的地形圖很吻合。進入20世紀90年代，星載合成孔徑雷達的不斷升空，提供了豐富的SAR數(shù)據(jù)，使InSAR技術(shù)在地形測繪上的應(yīng)用取得了很大的進展。2000年2月美國的“奮進號”航天飛機，使用InSAR技術(shù)，在短時間內(nèi)，就成功獲取了覆蓋地球陸地絕大部分的In-SAR數(shù)據(jù)，用其生成了相當于1∶50 000比例尺精度的高程模型。通過十幾年的發(fā)展，InSAR已在地震、地質(zhì)、水文、地面沉降、火山觀測等[4]領(lǐng)域取得了初步應(yīng)用，同時取得了一些研究和應(yīng)用成果。

目前，世界多國的科研機構(gòu)進行了InSAR軟件的研制，如Earthview，Doris等。商用軟件Earthview包括了InSAR數(shù)據(jù)處理模塊，利用其處理雷達影像對可以生成DEM。

本文以Earthview InSAR模塊提供的兩種不同的DEM生成方式處理雷達影像對數(shù)據(jù)，提取了DEM，并對結(jié)果進行對比分析。

2 InSAR提取DEM的基本原理

InSAR通過兩副天線觀測，同時觀測或兩次平行觀測，獲取地面同一地區(qū)的復數(shù)影像對，因為目標與天線位置的幾何關(guān)系，在復圖像上產(chǎn)生了相位差，形成干涉條紋圖，此圖中包含了斜距向上的點與兩天線位置之差的精確信息[5]。因此，利用傳感器高度、雷達波長、波速視向和天線基線距之間的幾何關(guān)系，可以精確地測量出圖像中每一點的三維位置，從而生成該地區(qū)的DEM。

2.1 InSAR提取DEM的幾何原理

圖1是InSAR成像幾何示意圖[6]。圖中A1和A2分別表示兩個天線，R1和R2(R2＝R1＋△R)分別表示兩天線到地面目標的斜距，H為A1的軌道高度，B為基線距，水平方向和B的夾角為α，水平基線和垂直基線分別為Bp和Bv。

圖1 InSAR成像幾何示意圖

從天線A1發(fā)射波長為λ的信號，天線A1和A2同時接收從地面目標S返回的信號，天線A1和A2接收的從目標S返回的信號之間的測量相位可以分別表示為[7]。

則兩副天線所接收到目標S的信號之相位差為:

由圖1中所示的幾何關(guān)系模型，利用余弦定理可得:

由于通常R1>>B，R1>>△R，式(4)可以化簡為:

由式(3)、式(5)，可以解出θ:

由圖1中的幾何關(guān)系，可以推出下面的關(guān)系式:

由式(6)和式(7)就可確定目標點S的高程。也是就說，如果已知天線位置(H，B，α)和雷達系統(tǒng)的參數(shù)(θ)等，就可以從φ計算出地表的高程值h。

2.2 InSAR提取DEM的數(shù)據(jù)處理流程

InSAR提取DEM的數(shù)據(jù)處理，主要步驟[8]包括:復圖像的粗、精配準、生成干涉圖、相位解纏、提取數(shù)字高程模型。其具體處理數(shù)據(jù)流程圖如圖2所示。

圖2 InSAR提取DEM的數(shù)據(jù)處理流程圖

(1)選擇合適的InSAR成像對。選擇InSAR成像對，要根據(jù)不同的干涉要求，進行選取。要求影像對必須相干。在選取成像對時，要考慮傳感器類型、空間和時間基線、成像區(qū)域地形大氣狀況等的影響。

(2)圖像對的配準。當來自鄰近軌道上的兩幅SAR圖像配準時，它們的相位差圖像會顯出條文，條紋的變化包含著地表地形信息。如果兩幅圖像沒有精確配準，它們生成干涉條紋就會模糊不清，甚至生不成干涉條紋。通常，圖像配準誤差必須在1/8個像元以下才對干涉條紋的質(zhì)量沒有明顯的影響。

(3)圖像對做共軛相乘，并計算相干系數(shù)。圖像配準后，圖像上像元數(shù)據(jù)可以看成由實部數(shù)據(jù)和虛部數(shù)據(jù)組成，將它們共軛相乘，得到高質(zhì)量的干涉圖。采用最大似然估計器計算相干系數(shù)，根據(jù)相干系數(shù)生成相干圖，以便對干涉圖進行質(zhì)量評價。

(4)去“平地效應(yīng)”，干涉濾波和二次采樣。干涉紋圖中包含有形變和地形相位信息，用InSAR進行地形測繪，需要把其中地形相位信息消除，即去“平地效應(yīng)”，一般采用相位補償法消除平地效應(yīng)。干涉紋圖中的相位受到多種噪聲的影響，需要進行濾波處理，通過halfband濾波消除頂?shù)臀灰茖ο辔粩?shù)據(jù)的影響，方位向濾波提高了干涉紋圖的信噪比。二次采樣可以減少后期數(shù)據(jù)處理的數(shù)據(jù)量。

(5)相位解纏。主輔影像對在進行共軛相乘后，得到的相位差實際上是纏繞相位，其取值在(－π，π)之間，必須經(jīng)過相位解纏才可以得到真實的干涉相位，把恢復真實相位的過程稱為相位解纏。它是InSAR數(shù)據(jù)處理過程中關(guān)鍵一步，其結(jié)果好壞直接影響InSAR的最終數(shù)據(jù)產(chǎn)品(如DEM或地面形變圖)的質(zhì)量。

(6)相位到高度的轉(zhuǎn)換。由式(6)和式(7)得到解纏后的相位到高度的轉(zhuǎn)換公式:

利用式(8)可以實現(xiàn)解纏相位到高程數(shù)據(jù)的轉(zhuǎn)換。

(7)地理編碼。在完成相位到高度轉(zhuǎn)換后，還需要經(jīng)過從斜距到地距的轉(zhuǎn)換和地理編碼過程才成為與地圖匹配的數(shù)字地形圖，供給最終用戶使用。地理編碼是將影像數(shù)據(jù)和高程數(shù)據(jù)，從雷達成像時的坐標系統(tǒng)(高程、方位、距離)轉(zhuǎn)換到某一種較為通用的參考坐標系。

3 InSAR在提取地面形變場DEM的應(yīng)用實例

Earthview InSAR模塊提供了多種雷達數(shù)據(jù)處理模式，其中，在生成DEM上，有DEM和DEM with External DEM兩種方式。Earthview InSAR模塊結(jié)構(gòu)如圖3所示。

圖3 Earthview InSAR模塊結(jié)構(gòu)示意圖

在Earthview InSAR模塊中，生成DEM的這兩種，兩者所不同的是在雷達影像數(shù)據(jù)處理時，是否引入外部數(shù)字高程模型(DEM)。本文中引入SRTM數(shù)據(jù)作為外部DEM(External DEM)，SRTM是美國“奮進”號航天飛機于2000年2月測量得到。它覆蓋了地球80%以上的陸地表面。目前NASA開放了部分SRTM數(shù)據(jù)，并在網(wǎng)上公開發(fā)布。SRTM數(shù)據(jù)可通過JPL提供的網(wǎng)址下載。覆蓋伊朗巴姆地區(qū)的SRTM數(shù)據(jù)如圖4所示。

圖4 覆蓋巴姆SRTM數(shù)據(jù)示意圖

在Earthview InSAR模塊中，DEM和DEM with External DEM兩種雷達數(shù)據(jù)處理模式的基本流程圖如圖5所示。

圖5 Earthview InSAR中兩種生成DEM雷達數(shù)據(jù)處理模式流程圖(a～b)

3.1 實驗數(shù)據(jù)介紹

伊朗巴姆地震發(fā)生在2003年12月26日，震級為6.6級，震中位置為N29.010，E58.260，給當?shù)貛砹司薮蟮娜藛T傷亡和財產(chǎn)損傷。ENVISAT衛(wèi)星在巴姆地震后，監(jiān)測到了該次地震的全過程。在歐空局網(wǎng)站上可以免費下載這些數(shù)據(jù)。本文利用ENVISAT衛(wèi)星以降軌方式獲取的2幅雷達影像數(shù)據(jù)進行提取DEM處理，并對結(jié)果進行分析。實驗所用雷達影像數(shù)據(jù)基本特征如表1所示。

雷達影像干涉基本信息 表1

3.2 實驗結(jié)果及其分析

本文采用Earthview InSAR模塊，處理雷達數(shù)據(jù)提取DEM，其實驗主要工作如下:

(1)數(shù)據(jù)輸入:輸入雷達影像對主從圖像，確定研究范圍，完成雷達影像數(shù)據(jù)的讀取分析，并生成主從圖像干涉基線距信息。

(2)配準:首先計算出中心點偏差，然后將其作為像元級配準的輸入?yún)?shù)，采用基于窗口的自動配準技術(shù)進行配準，最后通過主從圖像過采樣處理以減少后期數(shù)據(jù)處理的數(shù)據(jù)量、最大相干估算、偏移量擬合進行精配準。

(3)干涉處理:在干涉處理前，先進行干涉圖像對的頻率域、時間域預濾波處理，接著進行去平地效應(yīng)的消除，再者必須去除地形相位的影響。在DEM with External DEM方式中中，通過外部DEM模擬生成處理地形相位。在干涉相位信息增強處理中，主要是生成干涉相位圖及掩膜、相位坡度圖，進行干涉相位校正、殘余平地效應(yīng)和地形相位消除等。

(4)相位解纏和形變量解算。

(5)地理編碼:采用WGS－84高程、UTM投影、空間分辨率20 m×20 m為標準將視線向形變量及其水平、垂直分量分別進行地理編碼和校正，生成最終的DEM。

圖6 實驗結(jié)果示意圖(a～h)

在完成以上實驗處理工作后，就可以得到兩種生成DEM模式處理雷達影像對的處理結(jié)果，如圖6所示。

在增強干涉圖中，圖6(a)和圖6(b)進行了去平地效應(yīng)和濾波處理，能夠看出清晰的干涉條紋。圖6 (a)中干涉條紋比圖6(b)中干涉條紋密集，這表明在圖6(a)中有大量殘余相位的存在，造成干涉條紋密集，掩蓋了地形變化情況。圖6(b)由于殘余相位被大量去除，表現(xiàn)出了清晰的地形變化情況，然而，兩圖都存在相位模糊現(xiàn)象，不能用來表示地面起伏情況，還需要進行相位解纏。這說明引入外部DEM，可以幫助消除殘余相位，提高干涉圖的質(zhì)量。

在生成增強干涉圖的同時，可以得到相干圖，相干圖是最常用、最直觀的干涉質(zhì)量圖。相干圖既可以表示研究區(qū)域的相干性，作為相位解纏等重要參考指標；又可以反映兩幅雷達影像獲取期間的地物變化信息，對地物和地貌進行分類。一般來說，相干圖上越明亮的區(qū)域表示干涉質(zhì)量越高。圖6(c)和圖6(d)明亮的區(qū)域較多，即兩相干圖的相干系數(shù)高，表明實驗所用兩幅雷達影像相干性很高，也就是說實驗所采用的兩幅雷達影像的干涉質(zhì)量高。

圖6(e)和圖6(f)反映相位解纏后的干涉相位情況，顯然圖6(e)比圖6(f)干涉質(zhì)量差些，這是由于生成它的增強干涉圖含有大量的相位殘余，影響了相位解纏質(zhì)量。這樣就造成解纏后的圖6(e)不能準確反映所研究區(qū)域的地面起伏狀況，也給最終的DEM精度帶來影響。

圖6(g)和圖6(h)是最終DEM，通過比較可以看出圖6(h)比圖6(g)明亮光滑，這說明圖6(h)地理編碼后的DEM精度比圖6(g)地理編碼后的DEM高，能夠較精確地反映所研究區(qū)域的地面起伏。這表明借助外部DEM有助于提高提取地面形變場DEM的精度。

4 結(jié) 論

本文分析InSAR提取DEM基本原理基礎(chǔ)上，以Earthview InSAR模塊為數(shù)據(jù)處理平臺，采用兩種DEM生成模式分別提取了伊朗巴姆地震的DEM，并對結(jié)果進行對比分析，可以得出借助外部DEM，可以提高提取的地面形變場DEM的精度。

可以預見，隨著高效SAR傳感系統(tǒng)、SAR數(shù)據(jù)處理軟件、新技術(shù)等的出現(xiàn)，將會有高自動化、更高精度、更具適用性和可靠性的DEM出現(xiàn)。

[1] GRAHAM L C.Synthetic Aperture Radar for Topographic Mapping[J].Proc IEEE，1974(62):762～763[2] Zebker H A，Goldstien R M.Topographic Mapping from Interferometry Synthetic Aperture Radar Observation[J].Journal of Geophysical Reseach，1986(91):4990～4993

[3] 王志勇，張繼賢，張永紅.從InSAR干涉測量提取DEM [J].測繪通報，2007(7):27～28

[4] 王超，張紅，劉智.星載合成孔徑雷達干涉測量[M].北京:科技出版社，2002:165～198

[5] 仇春平，王堅.從SAR影像提取DEM的方法研究[J].測繪通報，2006(6):7～8

[6] 胡波，朱建軍，張長書.InSAR提取DEM的原理與實踐[J].測繪工程，2008(17):57～58

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

[8] 苗放，梁軍，葉成名等.用InSAR技術(shù)提取數(shù)字高程模型的研究[J].物探化探計算技術(shù)，2007(29):158～160

The Studies of DEM Extraction of Ground Deformation Field by InSAR

Yu JingBo1，Cao ZhenTan1，Xue Feng2，Zhao YuJiang3
(1.Geomatics College，UUST，Qingdao 266510，China； 2.Jinan Branch Eastdawn，Jinan 250101，China；3.GIC of exploration of urban and rural construction，Jinan 250031，China)

The basic principles of InSAR extracting DEM were introduced，and then DEMS of the Bam earthquake coseismic deformation field were extracted by data processing mode of two generation of Earthview InSAR to DEM processing ENVISAT satellite image date from the Bam earthquake.That the external DEM with high accuracy that could increase the accuracy of the extracted DEM of the ground deformation field was shown by comparison and analysis to the result.

Earthview InSAR；DEM generation mode；Bam earthquake；EnVISAT satellite；DEM accuracy

1672－8262(2011)02－92－04

P236

A

2010—10—22

余景波(1983—)，男，碩士研究生，主要從事現(xiàn)代測量數(shù)據(jù)處理理論及應(yīng)用。

海島(礁)測繪技術(shù)國家測繪局重點實驗室資助項目(2010A01)