基于短基線差分干涉法的上海地面沉降分析

張麗娜，伍吉倉，李 濤，陳 杰

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

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

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

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

1 短基線集方法

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

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

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

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

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

1.2 三維解纏

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

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

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

1.3 提取LOS平均沉降速度

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

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

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

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

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

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

上海是受地面沉降影響最嚴(yán)重的城市之一.從2000年起，地面沉降得到控制，目前上海市中心城區(qū)的年均沉降量在10 mm以內(nèi)[15]，線性趨勢較為明顯.近年來隨著城市大規(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.在處理時可對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é)果及分析

地形誤差相位與垂直基線成正比，因此在干涉圖組合時對垂直基線的要求比較嚴(yán)格.考慮到數(shù)據(jù)比較少，本次試驗(yàn)限定垂直基線800 m，時間基線1 000 d，干涉數(shù)據(jù)基線組合情況如圖2所示.圖2中每條連線代表一個干涉對組合，本次實(shí)驗(yàn)共產(chǎn)生了40個干涉對.滿足上面基線組合生成干涉圖如圖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ù)號表示沉降.在研究區(qū)約750 km的范圍內(nèi)，總共提取出了65 819個相干目標(biāo)，相干目標(biāo)點(diǎn)的密度達(dá)到了84個·km－2.從圖中可以看出，在2007—2010年間，圖像的上半部分，在虹口區(qū)與楊浦區(qū)均有沉降分布，楊浦區(qū)的四平路附近存在一個明顯的漏斗狀沉降；圖像中部也就是中心城區(qū)，總體比較平穩(wěn)，沒有明顯的沉降；圖像的下半部分，在閔行區(qū)東部有較大的沉降發(fā)生，而浦東新區(qū)的南部出現(xiàn)一個大范圍的不規(guī)則的沉降條帶.總體來看，中心城區(qū)沉降比較緩慢，較大的沉降主要發(fā)生在新開發(fā)區(qū)域，說明上海的沉降開始明顯地由中心城區(qū)向新城區(qū)和郊區(qū)轉(zhuǎn)移，這與實(shí)際情況是相符的.同時，也應(yīng)當(dāng)注意到，圖像中存在一些突變點(diǎn)，比如大片沉降區(qū)域中間夾雜著一些抬升速度較大的藍(lán)點(diǎn)，這可能是由于所選的相干目標(biāo)并不是穩(wěn)定點(diǎn)造成的.

從上述結(jié)果中，利用短基線方法監(jiān)測地面沉降的主要優(yōu)點(diǎn)有以下幾個方面：① 利用短基線方法，在沒有地面控制點(diǎn)的情況下，仍能夠提取城區(qū)地面沉降范圍和沉降速度信息；② 利用短基線方法，即使是在計(jì)算機(jī)的內(nèi)存有限的情況下，可以通過分塊分析再合并的方法，提取大范圍區(qū)域的沉降信息；③圖4中相干目標(biāo)點(diǎn)的密度遠(yuǎn)遠(yuǎn)大于傳統(tǒng)的水準(zhǔn)測量點(diǎn)密度.利用InSAR提取的沉降結(jié)果作為先驗(yàn)信息，對發(fā)生沉降的區(qū)域增加水準(zhǔn)觀測的頻率和密度.因此相對于水準(zhǔn)測量來說，利用短基線InSAR方法是一個非常有益的補(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ù)相對比.因只有上海部分地區(qū)的水準(zhǔn)數(shù)據(jù)，而且水準(zhǔn)點(diǎn)與ALOS的觀測時間不一致，故比較的結(jié)果是初步的.圖4中三角形表示水準(zhǔn)點(diǎn)位置.水準(zhǔn)測量得到的沉降速率與本文得到的沉降速率及其差值列于圖5中.

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

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

4 結(jié)論

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

從觀測結(jié)果可以看出，利用短基線方法監(jiān)測的優(yōu)點(diǎn)主要體現(xiàn)在兩個方面：（1）利用短基線方法，在沒有地面控制點(diǎn)的情況下，仍能夠提取大時間跨度、大區(qū)域、高精度的地面沉降信息.（2）短基線方法提取的相干目標(biāo)的密度遠(yuǎn)遠(yuǎn)大于傳統(tǒng)的水準(zhǔn)測量，能夠體現(xiàn)出沉降區(qū)域的范圍和趨勢.利用InSAR提取的大范圍沉降結(jié)果作為參考，對那些沉降觀測量比較大的區(qū)域增加水準(zhǔn)觀測的頻率和密度.因此，在地面沉降監(jiān)測應(yīng)用中，基于InSAR技術(shù)的短基線集方法是傳統(tǒng)水準(zhǔn)測量的有效補(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] 王艷，廖明生，李德仁，等.利用長時間序列相干目標(biāo)獲取地面沉降場[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ī)劃和國土資源管理局.上海市地質(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.