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

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

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

?

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

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

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

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

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

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

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

一、PSInSAR算法

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

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

(1)

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

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

(2)

(3)

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

(4)

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

二、數據處理與結果分析

1. 數據處理與形變分析

本文以常州市為研究對象，為驗證PSInSAR方法的可行性，試驗數據選取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日，時空基線參數見表1和表2，研究區域位置如圖1所示(由于TerraSAR-X覆蓋范圍較小(50km×30km)，因此基本未作切割)。其中，C波段EnvisatASAR利用Delft大學空間研究中心提供的精密軌道數據(orbitaldatarecords,ODR)。

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

表1　Envisat ASAR影像參數

表2　TerraSAR-X影像參數

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

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

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

2. 精度對比驗證

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

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

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

圖3　TerraSAR-X和Envisat ASAR同時段沉降量相關性

序號正方形內C波段/X波段PS點數C波段/X波段InSAR結果/(mm/a)對應水準結果/(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數據的影響

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

表4　TerraSAR-X和Envisat ASAR數據的傳感器參數對比

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

四、結束語

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

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

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

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

[2]王艷, 廖明生, 李德仁,等. 利用長時間序列相干目標獲取地面沉降場[J]. 地球物理學報, 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]李德仁，廖明生. 永久散射體雷達干涉測量技術[J]. 武漢大學學報(信息科學版), 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—)，男，博士生，研究方向為測量數據處理與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.