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

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

用InSAR研究地面形變場(chǎng)DEM的提取

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

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

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

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

1 引 言

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

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

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

2 InSAR提取DEM的基本原理

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

2.1 InSAR提取DEM的幾何原理

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

圖1 InSAR成像幾何示意圖

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

則兩副天線所接收到目標(biāo)S的信號(hào)之相位差為:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

在Earthview InSAR模塊中，生成DEM的這兩種，兩者所不同的是在雷達(dá)影像數(shù)據(jù)處理時(shí)，是否引入外部數(shù)字高程模型(DEM)。本文中引入SRTM數(shù)據(jù)作為外部DEM(External DEM)，SRTM是美國“奮進(jìn)”號(hào)航天飛機(jī)于2000年2月測(cè)量得到。它覆蓋了地球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兩種雷達(dá)數(shù)據(jù)處理模式的基本流程圖如圖5所示。

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

3.1 實(shí)驗(yàn)數(shù)據(jù)介紹

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

雷達(dá)影像干涉基本信息 表1

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

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

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

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

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

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

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

圖6 實(shí)驗(yàn)結(jié)果示意圖(a～h)

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

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

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

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

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

4 結(jié) 論

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

可以預(yù)見，隨著高效SAR傳感系統(tǒng)、SAR數(shù)據(jù)處理軟件、新技術(shù)等的出現(xiàn)，將會(huì)有高自動(dòng)化、更高精度、更具適用性和可靠性的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干涉測(cè)量提取DEM [J].測(cè)繪通報(bào)，2007(7):27～28

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

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

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

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

[8] 苗放，梁軍，葉成名等.用InSAR技術(shù)提取數(shù)字高程模型的研究[J].物探化探計(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)代測(cè)量數(shù)據(jù)處理理論及應(yīng)用。

海島(礁)測(cè)繪技術(shù)國家測(cè)繪局重點(diǎn)實(shí)驗(yàn)室資助項(xiàng)目(2010A01)