張麗娜，伍吉倉(cāng)，李 濤，陳 杰

（同濟(jì)大學(xué) 測(cè)量與國(guó)土信息工程系，上海200092）

合成孔徑雷達(dá)干涉測(cè)量（InSAR）是近年來(lái)發(fā)展非常迅速的微波遙感對(duì)地觀測(cè)技術(shù).因其具有全天候、大范圍、高空間分辨率的特點(diǎn)，不僅可以快捷地提取出m級(jí)的地形信息，而且能夠觀測(cè)探測(cè)到亞厘米級(jí)的微小地表形變[1].因此，InSAR研究方法廣泛地應(yīng)用在地震、火山、山體滑坡等突發(fā)形變監(jiān)測(cè)領(lǐng)域[2－5].

對(duì)于地面沉降、斷層構(gòu)造運(yùn)動(dòng)等長(zhǎng)時(shí)間跨度的緩慢形變，受時(shí)間去相干和大氣延遲影響比較嚴(yán)重[6－7].意大利學(xué)者通過(guò)對(duì)大量SAR影像的研究發(fā)現(xiàn)，在建筑物相對(duì)集中的城區(qū)和巖石裸露的山區(qū)，存在一些相干點(diǎn)，即使在干涉效果很差的情況下，其相位和幅度信息仍能在長(zhǎng)時(shí)間范圍內(nèi)保持穩(wěn)定.利用這些高相干點(diǎn)上的相位信息，就可以有效地克服時(shí)間去相干和大氣延遲的干擾，提取長(zhǎng)時(shí)間序列的地表形變信息[8].由此發(fā)展出一系列基于相干點(diǎn)目標(biāo)的長(zhǎng)時(shí)間序列分析方法，最具代表性的有兩種：永久散射體干涉測(cè)量方法（permanent scatterer InSAR，PSInSAR）[8]和 短 基線 集 方 法[9].PSInSAR 依 賴于主影像的質(zhì)量，并且基線過(guò)長(zhǎng)容易導(dǎo)致幾何失相干.短基線方法則采用多主影像，并且選擇基線較短的干涉對(duì)，提高了影像的利用率，增加了多余觀測(cè)信息.

本文首先介紹了基于短基線差分干涉測(cè)量方法提取地面沉降的基本原理和過(guò)程，主要討論了短基線方法的關(guān)鍵技術(shù).然后利用ALOS PALSAR數(shù)據(jù)，采用短基線干涉測(cè)量的方法提取了上海城區(qū)的沉降信息.因?yàn)锳LOS數(shù)據(jù)受軌道誤差和電離層延遲影像比較大.本文采用多項(xiàng)式擬合的方法，將殘余的軌道誤差和大氣誤差等空間相關(guān)的部分?jǐn)M合成一個(gè)平面，從干涉相位中去除，從而提高影像的質(zhì)量.

1 短基線集方法

短基線集干涉測(cè)量方法的主要思想是犧牲信息的空間密度來(lái)獲取可靠的形變相位信息：將所有SAR數(shù)據(jù)任意兩兩干涉處理，選擇垂直基線和時(shí)間基線較短的干涉對(duì)進(jìn)行研究，提取在時(shí)間序列中保持穩(wěn)定的相干目標(biāo).根據(jù)相干目標(biāo)各組成部分的在時(shí)間維和空間維的特征不同，將地表形變相位分離出來(lái)，求得視線向的沉降速度.主要計(jì)算流程可用圖1表示[10－13]：

圖1 短基線算法流程Fig.1 Flow chart of small baseline approach process

圖1各步驟中，短基線方法的關(guān)鍵在于相干點(diǎn)的識(shí)別、三維相位解纏的實(shí)現(xiàn)以及線性速度的擬合.下面將詳細(xì)說(shuō)明這些算法和過(guò)程.

1.1 識(shí)別相干目標(biāo)

本文考慮利用幅度和相位兩者的穩(wěn)定性來(lái)挑選相干目標(biāo).先用幅度穩(wěn)定性來(lái)初步篩選出候選相干點(diǎn).幅度穩(wěn)定性判斷采用Hooper等提出的振幅差離散度指標(biāo)[13].幅度離散度越小，說(shuō)明相干像元越穩(wěn)定.然后對(duì)候選點(diǎn)的相位進(jìn)行時(shí)間序列分析，利用相位的時(shí)間相關(guān)性測(cè)度來(lái)進(jìn)行計(jì)算，確定其相位穩(wěn)定性，剔除不合格點(diǎn)，最終得到穩(wěn)定的相干目標(biāo).

1.2 三維解纏

因?yàn)榻饫p相位和纏繞相位之間相差2π的整數(shù)倍，所以相位解纏實(shí)際上是一個(gè)求取最小化整數(shù)變量問(wèn)題.本文假設(shè)相位差在時(shí)間維和二維空間上都小于半個(gè)周期，實(shí)現(xiàn)近似三維解纏算法：先在時(shí)間維上解纏，然后把這個(gè)結(jié)果作為初始值，在二維空間上進(jìn)行優(yōu)化處理[14].

時(shí)間維的解纏.假設(shè)所有相鄰采樣點(diǎn)之間的相位差由大部分的平滑的變形信號(hào)和很少的隨機(jī)噪聲組成，而且其在時(shí)間維上小于半個(gè)周期，這樣經(jīng)過(guò)時(shí)間上的低通濾波之后，就可以先在時(shí)間維上對(duì)相位差進(jìn)行解纏，估計(jì)出相位差在時(shí)間維上的變化量.

空間域的解纏采用網(wǎng)絡(luò)最小費(fèi)用流法.利用時(shí)間域的解纏結(jié)果構(gòu)建先驗(yàn)概率函數(shù)，由先驗(yàn)函數(shù)的負(fù)對(duì)數(shù)函數(shù)提取費(fèi)用函數(shù)，然后用最小費(fèi)用流法（SNAPHU）來(lái)計(jì)算總體最小費(fèi)用解決方案.

1.3 提取LOS平均沉降速度

對(duì)解纏后的相位進(jìn)行時(shí)間域和空間域的處理，將大氣延遲、軌道誤差、地形誤差等干擾信號(hào)分離出去之后，這時(shí)的信號(hào)主要包含地形形變相位以及很少量的噪聲，可表示為

式中：φdef為干涉對(duì)兩次成像期間的形變相位；φnoi為殘余誤差相位和噪聲相位.

假設(shè)在主輔影像成像時(shí)間段Δt之間，地面點(diǎn)x發(fā)生了地形位移，一般認(rèn)為位移在時(shí)間上的變化為線性的，vx為點(diǎn)x的平均形變速度，λ為波長(zhǎng)，式（1）變?yōu)?/p>

令α＝4π·Δt／λ，則可得到：

若研究區(qū)中總共檢測(cè)出H個(gè)相干點(diǎn)，則可以組成H個(gè)方程.用最小二乘方法對(duì)方程組求解，即可得到視線向（LOS）的平均線性速度.

2 上海市區(qū)地面沉降實(shí)驗(yàn)及分析

上海是受地面沉降影響最嚴(yán)重的城市之一.從2000年起，地面沉降得到控制，目前上海市中心城區(qū)的年均沉降量在10 mm以內(nèi)[15]，線性趨勢(shì)較為明顯.近年來(lái)隨著城市大規(guī)模改建和郊區(qū)新城的發(fā)展，地面沉降開始出現(xiàn)在郊區(qū).本文選擇上海外環(huán)線以內(nèi)的城區(qū)作為實(shí)驗(yàn)區(qū)，以斯坦福大學(xué)開源軟件Sta MPS中的短基線模塊為基礎(chǔ)，利用日本ALOS衛(wèi)星的PALSAR數(shù)據(jù)開展上海地面沉降的研究.

2.1 實(shí)驗(yàn)數(shù)據(jù)

本文的DEM數(shù)據(jù)采用ASTER GDEM.SAR影像是日本的L波段的ALOS PALSAR數(shù)據(jù).PALSAR數(shù)據(jù)具有兩種：FBS和FBD.在處理時(shí)可對(duì)FBS數(shù)據(jù)進(jìn)行重采樣，使其與FBD數(shù)據(jù)具有相同的距離向分辨率.文中所用數(shù)據(jù)為2007—2010年間上海地區(qū)的17景PALSAR影像，見表1.

表1 實(shí)驗(yàn)所用的PALSAR數(shù)據(jù)Tab.1 PALSAR data

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

地形誤差相位與垂直基線成正比，因此在干涉圖組合時(shí)對(duì)垂直基線的要求比較嚴(yán)格.考慮到數(shù)據(jù)比較少，本次試驗(yàn)限定垂直基線800 m，時(shí)間基線1 000 d，干涉數(shù)據(jù)基線組合情況如圖2所示.圖2中每條連線代表一個(gè)干涉對(duì)組合，本次實(shí)驗(yàn)共產(chǎn)生了40個(gè)干涉對(duì).滿足上面基線組合生成干涉圖如圖3所示.

圖2 數(shù)據(jù)基線組合圖Fig.2 Baseline combination chart

圖3 全部干涉圖Fig.3 All interferograms

圖4 LOS沉降速度Fig.4 Subsidence rate in LOS

從圖3中，可以看出干涉圖整體顏色比較平滑.利用短基線方法得到所有相干目標(biāo)點(diǎn)上視線向的平均速度，如圖4所示.

由圖4可以看出，相干目標(biāo)的分布均勻，LOS向速度在－27.4—4.9 mm·年－1之間，其中負(fù)號(hào)表示沉降.在研究區(qū)約750 km的范圍內(nèi)，總共提取出了65 819個(gè)相干目標(biāo)，相干目標(biāo)點(diǎn)的密度達(dá)到了84個(gè)·km－2.從圖中可以看出，在2007—2010年間，圖像的上半部分，在虹口區(qū)與楊浦區(qū)均有沉降分布，楊浦區(qū)的四平路附近存在一個(gè)明顯的漏斗狀沉降；圖像中部也就是中心城區(qū)，總體比較平穩(wěn)，沒有明顯的沉降；圖像的下半部分，在閔行區(qū)東部有較大的沉降發(fā)生，而浦東新區(qū)的南部出現(xiàn)一個(gè)大范圍的不規(guī)則的沉降條帶.總體來(lái)看，中心城區(qū)沉降比較緩慢，較大的沉降主要發(fā)生在新開發(fā)區(qū)域，說(shuō)明上海的沉降開始明顯地由中心城區(qū)向新城區(qū)和郊區(qū)轉(zhuǎn)移，這與實(shí)際情況是相符的.同時(shí)，也應(yīng)當(dāng)注意到，圖像中存在一些突變點(diǎn)，比如大片沉降區(qū)域中間夾雜著一些抬升速度較大的藍(lán)點(diǎn)，這可能是由于所選的相干目標(biāo)并不是穩(wěn)定點(diǎn)造成的.

從上述結(jié)果中，利用短基線方法監(jiān)測(cè)地面沉降的主要優(yōu)點(diǎn)有以下幾個(gè)方面：① 利用短基線方法，在沒有地面控制點(diǎn)的情況下，仍能夠提取城區(qū)地面沉降范圍和沉降速度信息；② 利用短基線方法，即使是在計(jì)算機(jī)的內(nèi)存有限的情況下，可以通過(guò)分塊分析再合并的方法，提取大范圍區(qū)域的沉降信息；③圖4中相干目標(biāo)點(diǎn)的密度遠(yuǎn)遠(yuǎn)大于傳統(tǒng)的水準(zhǔn)測(cè)量點(diǎn)密度.利用InSAR提取的沉降結(jié)果作為先驗(yàn)信息，對(duì)發(fā)生沉降的區(qū)域增加水準(zhǔn)觀測(cè)的頻率和密度.因此相對(duì)于水準(zhǔn)測(cè)量來(lái)說(shuō)，利用短基線InSAR方法是一個(gè)非常有益的補(bǔ)充.

3 與水準(zhǔn)數(shù)據(jù)的比較

將LOS方向的沉降速度轉(zhuǎn)換到垂直方向，在距離水準(zhǔn)點(diǎn)100m的范圍內(nèi)搜索相干目標(biāo)點(diǎn)，選取距離最近的相干目標(biāo)與水準(zhǔn)數(shù)據(jù)相對(duì)比.因只有上海部分地區(qū)的水準(zhǔn)數(shù)據(jù)，而且水準(zhǔn)點(diǎn)與ALOS的觀測(cè)時(shí)間不一致，故比較的結(jié)果是初步的.圖4中三角形表示水準(zhǔn)點(diǎn)位置.水準(zhǔn)測(cè)量得到的沉降速率與本文得到的沉降速率及其差值列于圖5中.

圖5 相干點(diǎn)與水準(zhǔn)點(diǎn)的沉降速度對(duì)比Fig.5 Comparison of subsidence rate of coherent target and leveling

圖5中沉降速率觀測(cè)值為負(fù)號(hào)代表下降.對(duì)比三期觀測(cè)結(jié)果可以看出，隨著時(shí)間的推移，17水準(zhǔn)點(diǎn)里有14個(gè)點(diǎn)的沉降觀測(cè)量在變小，說(shuō)明在這些點(diǎn)附近沉降趨勢(shì)是減弱的，這與實(shí)際情況相吻合.本文重點(diǎn)對(duì)比分析2 0 0 1—2 0 0 6年的水準(zhǔn)觀測(cè)結(jié)果與2007—2010年的短基線觀測(cè)結(jié)果，有1，2，3，4，6，8，9，14等8個(gè)點(diǎn)的兩種觀測(cè)量差異很小，在8 mm以內(nèi)，說(shuō)明這些點(diǎn)上本文采用的短基線結(jié)果是比較可靠的；5，7，11，12，13，15這六個(gè)點(diǎn)，通過(guò)對(duì)比1995—2001和2001—2006年的沉降速率可知，兩期水準(zhǔn)數(shù)據(jù)也都相差較大，而且兩期水準(zhǔn)的差異與兩種觀測(cè)方法的差異很接近，說(shuō)明這些點(diǎn)的沉降速率變小了，本文的結(jié)果是比較符合事實(shí)的.10，16，17這三個(gè)點(diǎn)，其兩種觀測(cè)結(jié)果差異比較大但是兩期水準(zhǔn)數(shù)據(jù)的差異比較小，這說(shuō)明這些點(diǎn)附近的相干目標(biāo)（PS）點(diǎn)觀測(cè)量有誤差.經(jīng)分析，主要有兩個(gè)原因：（1）相干目標(biāo)點(diǎn)與水準(zhǔn)點(diǎn)的不匹配.本文在水準(zhǔn)點(diǎn)周圍100 m以內(nèi)搜索相干目標(biāo)點(diǎn)，在這個(gè)距離內(nèi)，即使地表地物相同，水準(zhǔn)點(diǎn)與相干目標(biāo)點(diǎn)的沉降情況也可能不同.（2）相干目標(biāo)點(diǎn)的選取和沉降計(jì)算過(guò)程中存在誤差.由于地形誤差、軌道誤差、大氣延遲等因素的影響，選擇相干目標(biāo)的時(shí)候可能引進(jìn)了部分不穩(wěn)定點(diǎn)，導(dǎo)致由InSAR所得到的相干目標(biāo)點(diǎn)的沉降速度不可靠.

4 結(jié)論

本文討論了基于短基線方法提取地面沉降的基本過(guò)程和原理，并以上海城區(qū)為例，利用PALSAR數(shù)據(jù)提取了上海城區(qū)的平均地面沉降速率.楊浦區(qū)存在一個(gè)明顯的漏斗狀沉降；在閔行區(qū)和浦東新區(qū)出現(xiàn)了一個(gè)較大范圍的不規(guī)則的沉降條帶.總體來(lái)看，中心城區(qū)的沉降比較緩慢，較大的沉降主要都發(fā)生在新開發(fā)區(qū)域，說(shuō)明上海的沉降開始明顯地由中心城區(qū)向新城區(qū)和郊區(qū)轉(zhuǎn)移.將短基線提取的沉降結(jié)果與水準(zhǔn)數(shù)據(jù)相比較，可以推斷本文提取的沉降速率是能夠達(dá)到厘米級(jí)或者更高的精度.

從觀測(cè)結(jié)果可以看出，利用短基線方法監(jiān)測(cè)的優(yōu)點(diǎn)主要體現(xiàn)在兩個(gè)方面：（1）利用短基線方法，在沒有地面控制點(diǎn)的情況下，仍能夠提取大時(shí)間跨度、大區(qū)域、高精度的地面沉降信息.（2）短基線方法提取的相干目標(biāo)的密度遠(yuǎn)遠(yuǎn)大于傳統(tǒng)的水準(zhǔn)測(cè)量，能夠體現(xiàn)出沉降區(qū)域的范圍和趨勢(shì).利用InSAR提取的大范圍沉降結(jié)果作為參考，對(duì)那些沉降觀測(cè)量比較大的區(qū)域增加水準(zhǔn)觀測(cè)的頻率和密度.因此，在地面沉降監(jiān)測(cè)應(yīng)用中，基于InSAR技術(shù)的短基線集方法是傳統(tǒng)水準(zhǔn)測(cè)量的有效補(bǔ)充.

致謝：本文采用的ALOS衛(wèi)星PALSAR數(shù)據(jù)由日本航天局提供.

[1] Prati C，Rocca F，Monti G A.SAR interferometry experiments with ERS－1[C]／／Proceedings First ERS－1 Symposium.Tokyo：ESA，1993：211－218.

[2] Prati C，F(xiàn)erretti A，Perssin D.Recent advance on surface ground deformation measurement by means of repeated spaceborne SAR observations[J].Journal of Geodynamics，2010，49：161.

[3] Massonnet D，F(xiàn)eigl K L.Radar interferometry and its application to changes in the earth’s surface[J].Reviews of Geophysics，1998，36（4）：441.

[4] Rosen P A，Hensley S，Joughin I R，et al.Synthetic aperture radar interferometry[J].Proceedings of the IEEE，2000，88（3）：333.

[5] Bamler R，Hartl P.Topical review：synthetic aperture radar interferometry[J].Inverse problems，1998，14（4）：R1.

[6] Joaquin M S，Ramon H，Bert M K，et al.Physical analysis of atmospheric delay signal observed in stacked radar interferometric data[C]／／Geoscience and Remote Sensing Symposium.[S.l.]：IEEE，2003：2112－2115.

[7] 王艷，廖明生，李德仁，等.利用長(zhǎng)時(shí)間序列相干目標(biāo)獲取地面沉降場(chǎng)[J].地球物理學(xué)報(bào)，2007，50（2）：598.WANG Yan，LIAO Mingsheng，LI Deren，et al..Subsidence velocity retrieval from long－term coherent targets in radar interferometric stacks[J].Chinese Journal Geophysics，2007，50（2）：598.

[8] Ferretti A，Prati C，Rocca F.Permanent scatterers in SAR interferometry[J].IEEE TRANS.Geoscience and Remote Sensing，2001，39（1）：8.

[9] Berardino P，F(xiàn)ornaro G，Lanari R，et al.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].IEEE Transactions on Geoscience and Remote Sensing，2002，40（11）：2375.

[10] 張紅，王超，吳濤，等.基于相干目標(biāo)的DInSAR方法研究[M].北京：科學(xué)出版社，2009.ZHANG Hong，WANG Chao，WU Tao，et al.Study of DInSAR approach based on coherent targets[M].Beijing：Science Press，2009.

[11] Hooper A.Persistent scatterer radar interferometry for crustal deformation studies and modeling of volcanic deformation[D].Delft：Department of Remote Sensing and Geoscience of Stanford University，2006.

[12] Hooper A，Zebker H A，Segall P，et al.A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers[J].Geophysical Research Letters，2004，31：L23611.

[13] Hooper A.StaMPS／MTI Manual[Z／OL]. [2009－11－24].http：／／radar.tudelft.nl／～ahooper／Stamps／StaMPS.

[14] Hooper A，Zebker H A.Phase unwrapping in three dimensions with application to InSAR time series[J].Optical Society of America，2007，24：2737.

[15] 上海市規(guī)劃和國(guó)土資源管理局.上海市地質(zhì)環(huán)境公報(bào)（2008年）[R／OL].[2009－05－25].http：／／www.shgtj.gov.cn.Shanghai Municipal Planning and Land Resources Administration.Shanghai geological environmental bulletin（2008）[R／OL].[2009－05－25].http：／／www.shgtj.gov.cn.