朱邦彥，李建成，儲征偉，唐　偉

(1. 武漢大學(xué)測繪學(xué)院，湖北 武漢 430079； 2. 南京市測繪勘察研究院有限公司，江蘇 南京 210019；3. 武漢大學(xué)測繪遙感信息工程國家重點實驗室，湖北 武漢 430079)

?

利用時序InSAR反演常州市地表沉降速率

朱邦彥1，李建成1，儲征偉2，唐偉3

(1. 武漢大學(xué)測繪學(xué)院，湖北 武漢 430079； 2. 南京市測繪勘察研究院有限公司，江蘇 南京 210019；3. 武漢大學(xué)測繪遙感信息工程國家重點實驗室，湖北 武漢 430079)

利用25景EnvisatASAR數(shù)據(jù)和29景高分辨率(3m分辨率)TerraSAR-X數(shù)據(jù)，采用永久散射體干涉測量技術(shù)(PSInSAR)研究了常州市2004—2014年的地表沉降速率。結(jié)果顯示，常州市的主要沉降區(qū)域發(fā)生在武進區(qū)，存在2個主要的沉降中心和1個大范圍的沉降條帶，2004—2010年間的地表沉降速率達到26mm/a，2010—2014年間沉降速率變緩，最大為24mm/a。兩組數(shù)據(jù)同時段的沉降量相關(guān)性達到0.78，并利用研究區(qū)域同期水準數(shù)據(jù)檢驗了本文的研究結(jié)果，兩者的平均速率差值均在5mm/a以內(nèi)，表明時序InSAR技術(shù)反演結(jié)果的可靠性和精度。X波段高分數(shù)據(jù)監(jiān)測到C波段無法監(jiān)測的高速路段存在5.3mm/a的沉降速率，與水準結(jié)果的RMS分別為2.5、4.2mm/a，表明TerraSAR-X比EnvisatASAR不僅具有更高密度的PS點，并且探測目標的位移具有更高的靈敏度和更高的精度。

PSInSAR；TerraSAR-X；高分辨率；失相干；永久散射體；地表沉降速率

常州市地處長江三角洲平原與太湖沖積平原交界處，由于地下水過度開采和城市的迅速發(fā)展，黏性土層壓縮固結(jié)、建筑容積率和密度增加，加劇了地面沉降，我國已有很多城市出現(xiàn)了地面塌陷、地裂縫等地質(zhì)災(zāi)害[1]。2003年底，常州市累計沉降量超過1000mm的區(qū)域面積達到102km2，累計沉降量超過800mm的區(qū)域面積達到218km2。利用時序InSAR技術(shù)監(jiān)測地表形變已經(jīng)取得了很多重要的研究成果[2]，為沉降監(jiān)測開辟了一條新道路。時序InSAR主要解決傳統(tǒng)InSAR技術(shù)受時間、空間去相干和大氣延遲影響的不足，目前主要方法有：一是永久散射體干涉測量方法(PSInSAR)[3-4]，最早由意大利Rocca教授研究小組提出[5]；二是短基線干涉測量方法(SBAS)[6]，最早由Berardino等于2002年提出[7]。Hooper等也提出了PS-InSAR和SBAS算法[8-9]。目前，時序InSAR技術(shù)已經(jīng)成功應(yīng)用于火山研究、滑坡探測、板塊運動、城市區(qū)域沉降等[10-13]。

PSInSAR主要思想是選取單一主圖像形成干涉圖，選取保持高相干性的點作為PS點(permanentscatterer)，采用一定方法構(gòu)成離散點網(wǎng)絡(luò)，通過分析PS點可靠的相位信息和時空特征，估算DEM誤差、大氣延遲及噪聲等，將其從差分干涉相位中逐個分離，最終獲得PS點的形變速率。Hooper等也對PSInSAR方法選取高相干目標點的方法進行了優(yōu)化[14]。傳統(tǒng)PSInSAR技術(shù)在纏繞相位中通過最優(yōu)搜索獲得形變模型的參數(shù)，而Hooper提出對時序相位進行三維時空解纏，從而獲取目標的形變速率。本文基于25景EnvisatASAR和29景TerraSAR-X數(shù)據(jù)，采用Hooper提出的PSInSAR算法，獲取常州地區(qū)2004—2014年期間地表沉降速率，與同期水準結(jié)果進行對比驗證，通過X波段和C波段的數(shù)據(jù)對比，分析高分辨率SAR數(shù)據(jù)對結(jié)果的影響。

一、PSInSAR算法

與傳統(tǒng)PSInSAR算法干涉配置方法相同，所有的影像均與選取的單一主影像(一般選擇時間中點即可)形成干涉圖。去除地形(DEM)相位后的相位值表示為

ψ=W(φd+φθ+φa+φo+φn)

(1)

式中，ψ為PSC(PS候選點)由主輔影像圖形成的干涉差分相位；W表示相位纏繞函數(shù)；φd為形變相位；φθ為側(cè)視角誤差；φa為大氣相位；φo為軌道誤差；φn為PSC散射信號、溫度噪聲和干涉配置誤差(如配準誤差)的貢獻總和。

首先，采用振幅離差閾值法，篩選出PSC，在距離向進行Goldstein濾波，并將地理編碼得到的PSC轉(zhuǎn)換到規(guī)則格網(wǎng)上(一般大小為100m)，確保格網(wǎng)間隔內(nèi)相位變化最小。根據(jù)相鄰相位之間的空間自相關(guān)性，利用低通濾波器過濾空間相關(guān)部分，得到

(2)

(3)

式中，λ為波長；B⊥,max為所有干涉圖中最大垂直基線距；ρ為傳感器到PSC所在像素的距離；θ為傳感器側(cè)視角。

(4)

最后，對PS點上的時序相位值進行三維時空解纏[15]，獲得絕對干涉相位值。相位解纏后，通過高通濾波器可以濾除部分大氣延遲相位，并獲得最終的形變相位。

二、數(shù)據(jù)處理與結(jié)果分析

1. 數(shù)據(jù)處理與形變分析

本文以常州市為研究對象，為驗證PSInSAR方法的可行性，試驗數(shù)據(jù)選取25景升軌EnvisatASAR影像和29景降軌TerraSAR-X影像，成像模式分別為IS2和Strip-Map，入射角范圍分別為16.8°—22.9°和21.3°—24.2°，時間范圍為2004年7月29日—2010年4月29日和2009年12月22日—2013年5月7日，時空基線參數(shù)見表1和表2，研究區(qū)域位置如圖1所示(由于TerraSAR-X覆蓋范圍較小(50km×30km)，因此基本未作切割)。其中，C波段EnvisatASAR利用Delft大學(xué)空間研究中心提供的精密軌道數(shù)據(jù)(orbitaldatarecords,ODR)。

首先采用振幅離差法，設(shè)閾值為0.5，篩選出PSC。然后，采用PSInSAR方法進行處理，EnvisatASAR共識別出223 920個PS點，每平方千米范圍內(nèi)包含約28個點，TerraSAR-X共識別出1 755 500個PS點，每平方千米包含約1170個PS點。最后，對PS點上的時序相位值進行三維時空解纏后，通過高通濾波器濾除主影像的大氣延遲相位，分別得到兩組時間段的研究區(qū)域地表沉降變化的空間分布特征，如圖2所示。

表1　Envisat ASAR影像參數(shù)

表2　TerraSAR-X影像參數(shù)

年沉降速率和累計沉降量在一定程度上反映了本地區(qū)總體的沉降速率和趨勢。本文將平均沉降速率和累計沉降量統(tǒng)一基準，由相位解纏后的PS點相位值分析，設(shè)以(120°E，31.8°N)為中心、半徑500 m的圓為參考穩(wěn)定區(qū)域，Envisat ASAR共44個參考PS點，TerraSAR-X共720個參考PS點，其他所有PS點都以該區(qū)域為參考基準(如圖2所示)。結(jié)果顯示，2004年7月—2013年5月間，常州市地表沉降呈不均勻分布，北部地區(qū)及周圍的山丘地帶未發(fā)生明顯沉降(地下水的抽取可以迅速得到長江水的補給有關(guān))。雖然2000年已經(jīng)禁止開采地下水，但是受到地下水抽取的延續(xù)性的影響，同時地下水補給微弱，且隨著工業(yè)化的發(fā)展和城市化的建設(shè)，建筑物容積率和荷載的增加，使沉降量和沉降區(qū)域面積也相應(yīng)增大。

圖2　研究區(qū)域的LOS方向平均沉降速率

常州市沉降嚴重區(qū)域主要發(fā)生在武進區(qū)，Envisat ASAR數(shù)據(jù)處理結(jié)果顯示，2004—2010年間，江宜高速的東側(cè)沿湖濱路，在牛塘鎮(zhèn)和盧家巷的區(qū)域形成了一個沉降漏斗，漏斗中心的LOS(line of sight)方向平均沉降速率達到22 mm/a，沉降中心區(qū)的面積達到15 km2。沿著淹城路，從東寶路與滆湖西路交界處到港東村也形成了一個沉降漏斗，漏斗中心的LOS方向平均沉降速率達到23 mm/a，沉降中心區(qū)的面積達到7 km2。從武宜路與西湖路交界處至戴溪是沉降重災(zāi)區(qū)，形成了大范圍的沉降條帶區(qū)，與周邊若干個小漏斗基本已經(jīng)相互連接。沉降條帶區(qū)的LOS方向平均沉降速率最大達到26 mm/a，沉降條帶中心區(qū)的面積達到80 km2。TerraSAR-X數(shù)據(jù)處理結(jié)果顯示，2010—2014年間，由于地下水抽取的影響逐漸減弱，以及城市基礎(chǔ)建設(shè)的日趨完善，沉降區(qū)域的位置和面積未發(fā)生大的偏移和改變，但是沉降速率較2004—2010年變緩，兩個沉降漏斗中心的LOS方向年平均沉降速率均降到16 mm/a左右，沉降趨勢主要轉(zhuǎn)移至大范圍的沉降條帶，其中心平均沉降速率為24 mm/a(X波段數(shù)據(jù)范圍并未完全覆蓋沉降區(qū)，但沉降中心區(qū)域已經(jīng)覆蓋)。結(jié)合C波段和X波段的兩組數(shù)據(jù)，從2004年7月—2013年5月，兩個沉降漏斗中心的累計沉降量最大分別為180和189 mm，沉降條帶中心累計沉降量最大為229 mm。

2. 精度對比驗證

為了評價時序InSAR技術(shù)反演的可靠性和精度，本文采用兩種方法對比驗證。

第一，利用兩組不同波段相應(yīng)時段的沉降量分析相關(guān)性。利用Envisat ASAR和TerraSAR-X數(shù)據(jù)進行PSInSAR處理，在大沉降條帶區(qū)域內(nèi)，選擇兩者相同范圍內(nèi)的399個PS點，在前者PS點為中心、邊長為100 m的正方形內(nèi)，搜索后者PS點并取均值，利用兩者PS點上的時序解纏相位，計算出相應(yīng)時段2009年12月—2010年5月的沉降量并進行對比，得到相關(guān)性分析結(jié)果，兩者擬合得到的關(guān)系式為y=0.7812x-3.2882(如圖3所示)。圖3表明，Envisat ASAR與TerraSAR獲取的時段沉降量的趨勢一致，相關(guān)性系數(shù)達到了0.78，說明了時序InSAR技術(shù)反演結(jié)果的可靠性。

第二，利用研究區(qū)域的4個水準點的沉降速率對比驗證。設(shè)以水準點為中心、邊長為200 m的正方形內(nèi)形變梯度一致，根據(jù)反演結(jié)果，對范圍內(nèi)的PS點部分取中位數(shù)和作算術(shù)平均，考慮到LOS方向的側(cè)視角與入射角不同，每個PS點的入射角約為側(cè)視角加3°，利用該范圍內(nèi)的PS點時序分析得到的平均沉降速率進行對比，結(jié)果見表3。在正方形參考范圍內(nèi)，X波段TerraSAR-X獲取的PS點比C波段Envisat ASAR獲取的點密度高很多。L3點位于沉降漏斗中心區(qū)域，沉降速率較快，2004—2010年期間，Envisat ASAR獲取的結(jié)果與參考水準差值為4.9 mm/a，累計沉降量達到180 mm，2010—2014年間，TerraSAR-X獲取的結(jié)果與水準差值為0.4 mm/a，累計沉降量達到50 mm。L1點、L2點和L4點位于沉降漏斗邊緣，沉降速率較緩，Envisat ASAR獲取的結(jié)果與參考水準差值分別為2.0、5.0、4.2 mm/a，TerraSAR-X獲取的結(jié)果與水準差值分別為-1.2、4.7、0.6 mm/a。由此可見，時序InSAR反演地表沉降速率與水準結(jié)果非常吻合，驗證了試驗結(jié)果的可靠性和精度。

圖3　TerraSAR-X和Envisat ASAR同時段沉降量相關(guān)性

序號正方形內(nèi)C波段/X波段PS點數(shù)C波段/X波段InSAR結(jié)果/(mm/a)對應(yīng)水準結(jié)果/(mm/a)差值/(mm/a)L17/131-6.8/-2.9-8.8/-1.72.0/-1.2L26/315-11.8/-7.2-16.8/-11.95.0/4.7L34/127-26.1/-14.9-31/-15.34.9/0.4L47/268-15.0/-10.8-19.2/-11.44.2/0.6

三、高分辨率SAR數(shù)據(jù)的影響

近幾年發(fā)射的高分辨率X波段的SAR傳感器在地理編碼的精度和對目標微小位移的探測方面具有明顯優(yōu)勢。相比之前歐洲空間局發(fā)射的C波段的Envisat，德國航空航天中心發(fā)射X波段的TerraSAR-X SAR傳感器具有更短的重復(fù)周期和更高的空間分辨率。ESA的ERS數(shù)據(jù)是最先用來進行PSInSAR分析的，隨后C波段Envisat的空間分辨率達到30 m，重復(fù)周期為35 d，這使得針對相位信息提取分析的算法和技術(shù)得以發(fā)展、實現(xiàn)和運用。為了分析兩種波段數(shù)據(jù)的差異，表4為本文用到的兩組數(shù)據(jù)傳感器參數(shù)的對比，分別為常州市相同區(qū)域的Envisat ASAR和TerraSAR-X數(shù)據(jù)。TerraSAR-X重復(fù)周期為11 d，更短的重復(fù)周期減弱了時間失相干的影響，而且分辨率小于3 m，更高的空間分辨率能更好體現(xiàn)散射體細節(jié)的變化。X波段的波長為λ=3.1 cm，大致是C波段的一半，這使得探測目標的位移具有更高的靈敏度(大約是2倍)，并且使位移的探測精度得到提升。

表4　TerraSAR-X和Envisat ASAR數(shù)據(jù)的傳感器參數(shù)對比

由于一個PS點取決于分辨單元內(nèi)占主導(dǎo)地位的后向散射信號，也就意味著SAR圖像空間分辨率與PS點密度有著密切的關(guān)系。分辨率越高，散射中心主導(dǎo)周圍散射體的可能性越大，從而PS點越密集。從另一個角度講，X波段的SAR干涉圖在人工地物較少區(qū)域卻更容易受到時間失相干的影響，這是因為X波段對散射體分布的差異導(dǎo)致的相位變化更加敏感。但是，在具有人工地物和裸露巖石的區(qū)域，高分辨率數(shù)據(jù)獲得PS點的空間密度比C波段的高得多。如表3所示，4個水準點相同范圍內(nèi)的PS點個數(shù)，C波段和X波段識別的PS點數(shù)差兩個數(shù)量級，這主要由兩種因素引起：一是更短的重復(fù)周期，減弱了時間去相干的影響；二是X波段的SAR傳感器具有更高的軌道穩(wěn)定性，容許相對更低的空間基線。圖4為研究區(qū)局部放大圖，將獲取的PS點疊加到衛(wèi)星圖像上，可以看出東西方向橫跨滆湖的常合高速路段，C波段數(shù)據(jù)基本沒有識別的PS點，而X波段數(shù)據(jù)識別出2382個PS點，2010—2013年間，跨滆湖路段的LOS方向平均沉降速率在5.3 mm/a以內(nèi)?？梢姼叻直媛蔛AR數(shù)據(jù)可以用來監(jiān)測公路、橋梁等目標微小的位移和反映其細節(jié)的變化。并且，利用Envisat ASAR和TerraSAR-X數(shù)據(jù)通過PSInSAR技術(shù)識別出來的PS點，與研究區(qū)域內(nèi)的4個參考水準點作比較之后，RMS分別為4.2和2.5 mm/a，說明高分辨率SAR數(shù)據(jù)探測目標的位移具有更高的精度。

四、結(jié)束語

本文利用PSInSAR時序分析方法，基于25景Envisat ASAR和29景TerraSAR-X數(shù)據(jù)，得到了2004—2014年間的地表沉降變化的空間分布特征，利用兩組不同波段相應(yīng)時段的數(shù)據(jù)分析了相關(guān)性，并與水準點進行了對比驗證，最后分析了X波段高分數(shù)據(jù)對結(jié)果的影響。結(jié)果顯示，常州市地面沉降主要發(fā)生在武進區(qū)，存在2個主要的沉降中心和1個大范圍的沉降條帶，2010—2014年間，2個沉降漏斗中心的沉降速率降至16 mm/a，大范圍沉降條帶中心的沉降速率仍然有24 mm/a，成為常州市的主要沉降區(qū)域。兩組不同波段的數(shù)據(jù)相關(guān)性達到0.78，與水準點的差值均在5 mm/a以內(nèi)，驗證了時序InSAR技術(shù)在常州市地表沉降研究中的可靠性和精度。X波段高分辨率的數(shù)據(jù)識別PS點的密度比C波段的高出2個數(shù)量級，能夠監(jiān)測到跨滆湖路段存在平均沉降速率5.3 mm/a的變化，前者與水準點的RMS比后者低將近1/2。

圖4　相同研究區(qū)域的LOS方向PS點平均沉降速率(底圖為衛(wèi)星圖像)

由本文的結(jié)果看出，常州市武進區(qū)地面沉降嚴重，需要繼續(xù)加強地下水和施工管理，嚴格控制建筑物的密度、容積率，采取相應(yīng)的管理措施和預(yù)防手段。時序InSAR技術(shù)在城市區(qū)域沉降監(jiān)測方面具有獨特的優(yōu)勢，并且高分辨率SAR數(shù)據(jù)在監(jiān)測公路、橋梁等目標微小的位移和細節(jié)方面具有巨大的應(yīng)用前景。

[1]李永生, 張景發(fā), 李振洪, 等. 利用短基線集干涉測量時序分析方法監(jiān)測北京市地面沉降[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版)，2013, 38(11): 1374-1377.

[2]王艷, 廖明生, 李德仁,等. 利用長時間序列相干目標獲取地面沉降場[J]. 地球物理學(xué)報, 2007, 50(2): 598-604.

[3]FERRETTI A, PRATI C, ROCCA F. Permanent Scatterers in SAR Interferometry [J]. IEEE Transaction on Geoscience and Remote Sensing, 2001, 39(1) :8-20.

[4]李德仁，廖明生. 永久散射體雷達干涉測量技術(shù)[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版), 2004, 29(8): 664- 667.

[5]FERRETTI A, PRATI C, ROCCA F. Nonlinear Subsidence Rate Estimation Using Permanent Scatterers in Differential SAR Interferometry [J]. IEEE Transactions on Geoscience and Remote Sensing. 2000, 38(5): 202-212.

[6]LANARI R, MORA O. A Small-baseline Approach for Investigating Deformations on Full-resolution Differential SAR Interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(7): 1377-1386.

[7]BERARDINO P, FORNARO 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-2383.

[8]HOOPER A, ZEBKER H, 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): 1-5.

[9]HOOPER A. Persistent Scatterer Radar Interferometry for Crustal Deformation Studies and Modeling of Volcanic Deformation [D]. Stanford: Stanford University, 2006.

[10]HOOPER A, SEGALL P, ZEBKER H. Persistent Scatterer Interferometric Synthetic Aperture Radar for Crustal Deformation Analysis, with Application to Volcan Alcedo, Galapagos [J].Journal of Geophysical Research, 2007，112(B07407): 1-21.

[11]HILLEY G E, BURGMANN R, FERRETTI A，et al. Dynamics of Slow-moving Landslides from Permanent Scatterer Analysis[J]. Science, 2004, 304(5679): 1952-1955.

[12]YANG G H, JIANG Z S, LIU G Y, et al. Possible Relation of Horizontal Movement Field in North China to Kunlun Mountain Ms 8.1 Earthquake [J].Journal of Geodesy and Geodynamics, 2007, 27(2): 10-15.

[13]PERISSIN D, FERRETTI A. Urban Target Recognition by Means of Repeated Spaceborne SAR Images[J]. IEEE Transaction on Geoscience and Remote Sensing, 2007, 45(12): 4043-4058.

[14]HOOPER A. A Multi-temporal InSAR Method Incorporating both Persistent Scatterer and Small Baselines Approaches [J]. Geophysical Research Letters, 2008, 35(L16302): 1-5.

[15]HOOPER A, ZEBKER H. Phase Unwrapping in Three Dimensions with Application to InSAR Time Series[J]. Journal of the Optical Society of America A: Opeics, Image Science, and Vision, 2007, 24(9): 2737-2747.

InversionofEarthSurfaceDisplacementRatesBasedonTemporalInSARMethodinChangzhouCity

ZHUBangyan，LIJiancheng，CHUZhengwei，TANGWei

2015-09-23

國家973計劃(2013CB733300)

朱邦彥(1986—)，男，博士生，研究方向為測量數(shù)據(jù)處理與InSAR時間序列分析。E-mail：byz@whu.edu.cn

P237

B

0494-0911(2016)05-0026-06

引文格式： 朱邦彥，李建成，儲征偉，等. 利用時序InSAR反演常州市地表沉降速率[J].測繪通報，2016(5)：26-31.DOI:10.13474/j.cnki.11-2246.2016.0148.