基于SBAS-InSAR的北京地區地表沉降監測與分析

周　呂　郭際明　李　昕　胡紀元

1　武漢大學測繪學院，武漢市珞喻路129號，430079 2　北京市測繪設計研究院城市空間信息工程北京市重點實驗室，北京市羊坊店路15號，100038 3　桂林理工大學廣西空間信息與測繪重點實驗室，桂林市雁山街319號，541004 4　武漢大學精密工程與工業測量國家測繪地理信息局重點實驗室，武漢市珞喻路129號，430079

?

基于SBAS-InSAR的北京地區地表沉降監測與分析

周呂1,2,3郭際明1,4李昕1胡紀元1

1武漢大學測繪學院，武漢市珞喻路129號，430079 2北京市測繪設計研究院城市空間信息工程北京市重點實驗室，北京市羊坊店路15號，100038 3桂林理工大學廣西空間信息與測繪重點實驗室，桂林市雁山街319號，541004 4武漢大學精密工程與工業測量國家測繪地理信息局重點實驗室，武漢市珞喻路129號，430079

運用SBAS-InSAR獲取北京地區的地表沉降信息，采用18景ENVISAT ASAR影像完成北京地區2007～2010年地表沉降的時空分析。結果表明，北京地區沉降不均勻較為嚴重，在昌平區、順義區、通州區等區域出現多處沉降漏斗，且有連成一片并向東擴張的趨勢；大部分地區的平均沉降速率在-150 ～10 mm/a，沉降中心的最大沉降量超過400 mm；地表沉降受地下水開采與城市化影響明顯。

地表沉降；SBAS-InSAR；北京地區；時空分析

對于地表形變監測與分析，傳統的監測方法(如水準測量、GNSS測量、全站儀測量、分層標測量等[1])可以獲取監測點較高時間分辨率與測量精度的形變量，但監測結果空間分辨率低，難以有效地監測與分析整個城市的區域性形變。合成孔徑雷達差分干涉測量(differential interferometric synthetic aperture radar, D-InSAR)技術可以探測亞cm級的地表沉降，但其易受時間、空間失相干與大氣延遲的影響，較難完成高精度的長時間間隔地表監測[2]。永久散射體合成孔徑雷達干涉測量(persistent scatterer interferometric synthetic aperture radar, PS-InSAR)技術通過對高相干目標點進行相位分析，較好地克服了失相干與大氣延遲影響，但需要較多的SAR影像數據[3-4]。小基線集合成孔徑雷達干涉測量(small baseline subset interferometric synthetic aperture radar, SBAS-InSAR)技術將SAR影像數據組成若干個子集，采用最小二乘法求解子集的形變時間序列，同時利用奇異值分解法(singular value decomposition, SVD)將多個小基線集聯合求解，獲取整個時間跨度的形變序列，相對于PS-InSAR技術，SBAS-InSAR需要的SAR影像數目較少且獲取非線性形變信息的能力較強[5-12]。

本文選取覆蓋北京地區的18景C波段的ENVISAT ASAR 影像數據，采用SBAS-InSAR技術對地表沉降進行形變監測與分析，獲取了研究區域2007～2010年的沉降分布情況與沉降速率場，并驗證了SBAS-InSAR技術監測城市地表沉降形變的可行性。

1　SBAS-InSAR技術原理

假設獲取了N+1景SAR影像，且影像獲取時間t按時間順序(t0,…,tN)排列,依據干涉基線組合可以生成M幅干涉圖，且M滿足：

(1)

假設干涉圖j由tA時間獲取的影像與tB時間獲取的影像進行干涉生成(tB>tA),在去除平地效應與地形相位影響后，干涉圖j中距離向為r與方位向為x的某一像素的干涉相位可表示為：

(2)

式中，φ(tB,x,r)和φ(tA,x,r)分別為tB與tA時刻SAR影像上的相位值，φdef,j(x,r)為tA時刻至tB時刻間衛星視線向的形變相位，φtopo,j(x,r)為參考DEM不精確引起的地形相位誤差，φatm,j(x,r)為大氣相位誤差，φnoise,j(x,r)為噪聲相位。其中，φdef,j(x,r)、φtopo,j(x,r)和φatm,j(x,r)可表示為：

(3)

式中，λ為雷達傳播信號的波長，d(tB,x,r)和d(tA,x,r)分別是tB與tA時刻相對于參考時刻t0的雷達視線向的累積形變量，B⊥j為干涉圖j的垂直基線，Δz為DEM誤差，R為雷達與目標點之間的距離，θ為入射角，φatm,j(tB,x,r)和φatm,j(tA,x,r)分別為tB與tA時刻SAR影像中的大氣相位分量。

為了獲取研究區域的形變序列，需要精確估計出地形相位誤差分量、大氣延遲相位誤差分量以及噪聲相位分量，并將這3個分量從干涉相位δφj(x,r)中去除。

對于整個集合中的所有干涉圖，在去除各項誤差分量之后，由式(2)可以得到一個方程組系統，其矩陣形式為：

(4)

式中，A為M×N的系數矩陣，且?j=1,…,M,M對應于干涉圖數量，N對應于SAR影像數量，φT=[φ(t1),…,φ(tN)]為每一景SAR影像中高相干點對應的相位值所組成的向量,δφT=[δφ1,…,δφM]為各差分干涉圖對應的解纏相位值所組成的向量。

為求解研究區域各高相干點的形變速率，可用兩景影像間的平均相位速率代替相位值，則式(4)變為：

(5)

式中，B為M×N的系數矩陣，vT可以表示為：

(6)

當系數矩陣B為滿秩(即M≥N)時，可用最小二乘法則求解出形變速率；而當M

2　研究區概況

北京市位于華北平原，地勢較為平坦，平均海拔為43.5 m。北京市的平原地區主要是由永定河、潮白河、溫榆河、泃河等河流聯合作用而形成的山前洪積、沖積平原。該地區多處沉降區位于平原地區。研究區域如圖1中的黑色方框所示，該區域西臨北京西山，東臨三河市燕郊鎮，南臨廊坊市，北臨順義區高麗營鎮，區域內多為平原，地形起伏較小。研究區域的中心經緯度為39°51′N、116°28′E，面積約為4 013 km2。

圖1　北京地區地理位置Fig.1　The geographical location of Beijing area

3　數據處理與結果分析

3.1數據處理

選取2007-08-01～2010-09-29的18景覆蓋研究區域的ENVISAT ASAR 0級影像數據作為實驗數據，其中影像數據采用波長為5.6 cm 的C波段，軌道方向為降軌，方位向與距離向的分辨率分別為4.570 m 與7.804 m，影像中心的入射角約為20.8°，極化方式為VV。采用美國宇航局提供的分辨率為90 m×90 m的SRTM DEM數據去除地形相位影響，同時利用歐洲空間局發布的DORIS軌道數據提高影像的軌道數據精度。

數據處理過程如下：1)將影像裁剪為本文的研究區域，對影像進行1×5(距離向×方位向)多視處理，選取2009-10-14獲取的影像為公共主影像，將所有輔影像進行配準并重采樣至公共主影像；2)選取時間基線閾值為400 d和空間基線閾值為900 m進行差分干涉圖處理，共生成64對小基線差分干涉圖集(圖2，圓點代表SAR影像，線段代表干涉圖，方框代表公共主影像)；3)采用DORIS軌道數據去除平地效應，同時利用SRTM3 DEM 數據消除地形相位；4)濾波處理，消除相關噪聲；5)采用最小費用流法對小基線干涉圖集進行相位解纏，結果見圖3；6)依據振幅與相位的穩定性選擇研究區域內的PS點，去除由DEM誤差引起的相位分量，通過奇異值分解(SVD)法求解高相干目標點的沉降速率。

圖2　時空基線結構Fig.2　The spatial-temporal structure of baselines

求解研究區的形變速率場時，選取區域內一個穩定點作為參考點R(依據已有的水準測量數據選取)。假定該點的形變速率為0，則各PS點的沉降速率均是相對于該參考點而言。獲取研究區各PS點的沉降速率后，采用克里金插值算法計算出整個研究區域內的沉降形變速率場。最后，依據監測時段內速度的積分便可得到該時段的形變量。

3.2結果分析

圖3為對各小基線差分干涉圖進行解纏后的相位形變圖，正值(藍色)表示在視線方向上輔影像相對于主影像沉降，負值(紅色)表示輔影像相對于主影像上升。從各時間段內的相位形變圖可以看出，研究區內形變較明顯，存在沉降漏斗。研究區內共識別出202 742個高相干目標點，平均每1 km2包含約51個高相干目標點。

圖3　解纏后的小基線差分干涉圖Fig.3　Small baseline interferograms stack after phase unwrapping

圖4為2007-08-01～2010-09-29研究區的沉降形變平均速率圖。可以看出，北京地區從北至南的主要沉降區分別為昌平區、順義區、朝陽區、通州區和大興區，與Ng等[13]采用PS-InSAR技術利用ALOS PALSAR 數據獲取的形變結果一致，驗證了本文采用SBAS-InSAR技術監測北京地區地表沉降形變的可行性與可靠性。

圖4　研究區LOS向平均速率Fig.4　Mean LOS velocity of the studied area

為進一步驗證實驗結果的可靠性，收集研究區域內2010～2011年的12個水準點測量數據(水準點L1至L12的分布見圖4)，將水準測量與SBAS-InSAR重疊時段的處理結果進行對比，見圖5。可以看出，水準測量結果與SBAS-InSAR獲取的結果趨勢符合較好。對比分析兩種方法的數據結果可以發現，兩種方法的最大相對誤差為 7.1 mm/a, 最小相對誤差為 -0.4 mm/a，說明兩種方法獲取的結果具有較好的一致性。

圖5　水準測量與SBAS-InSAR結果對比Fig.5　Comparison between SBAS-InSAR results and leveling results

沉降區主要分布于潮白河、溫榆河和泃河流域的沖積、洪積扇平原上，昌平沉降區、順義沉降區、朝陽沉降區與通州沉降區逐漸連成一片并有向東擴張的趨勢，這與近幾年北京地區的城市化擴張緊密相關。同時，在東部出現了新的沉降區即燕郊鎮沉降區。

由圖4可知，研究區內大部分地區的平均沉降速率在-150～10 mm/a之間，不均勻沉降情況明顯。北京中心城區、海淀區、房山區、豐臺區與石景山區的地表沉降量小且相對穩定，大部分地區的沉降速率小于10 mm/a，這與市區嚴格控制地下水的開采、北京地區西部的地表沉積物多為壓縮性較小的粗顆粒沙卵礫石有關。昌平區出現多處沉降漏斗，其中沙河鎮與上莊鎮沉降中心較為嚴重，年沉降速率超過70 mm/a，且兩處沉降中心的最大累計沉降量均超過180 mm。朝陽區的來廣營鄉、孫河鄉和金盞鄉一帶沉降較為明顯，其中金盞鄉沉降區最為嚴重，該區域最大平均沉降速率超過130 mm/a，累計最大沉降量達到300 mm。

通州區的管莊鄉、三間房鄉和豆各莊鄉一帶(圖4中的區域S)沉降較為嚴重。圖6為通州區的沉降帶區域S的沉降形變速率情況。圖7為區域S淺地表空間的利用情況，該區域內地鐵八通線、地鐵6號線、通惠河沿東西向通過，并有多條高速公路與鐵路穿過，同時其地下水開采較為嚴重，因此其沉降可能同時受動載荷、地下空間應用、地質構造以及地下水開采等影響。對比圖6與圖7可以發現，沉降中心主要位于同時有地鐵、公路、鐵路以及河流穿過的區域，且最大年平均沉降速率超過150 mm/a，最嚴重沉降中心的累計沉降量超過400 mm，不均勻沉降較為明顯，導致該地區出現部分建筑物墻面開裂現象。說明復雜的淺地表空間利用對地表沉降具有一定的貢獻率，同時地表沉降嚴重時會威脅建筑物的安全。

圖6　S區域的沉降形變速率Fig.6　Subsidence deformation velocity in S area

圖7　S區域的淺地表空間利用情況Fig.7　Shallow surface space utilization in S area

相對于通州沉降區、朝陽沉降區與昌平沉降區，順義沉降區與大興沉降區的沉降量較小。

2007～2010年，北京地區地下水開采量逐年增長，使得地下水位基本處于持續下降狀態。同時，受經濟發展、人口增加以及政策影響，北京的城市化面積逐漸增長，使得北京地區的地表沉降日趨嚴重，通州沉降區、朝陽沉降區、昌平沉降區等較為嚴重的沉降漏斗的地表沉降有逐漸連成一片并向東擴張的趨勢。

4　結　語

本文采用SBAS-InSAR技術，利用18景ENVISAT ASAR 數據對北京地區的地表沉降形變進行時空分析，獲取了該地區2007-08-01～2010-09-29的平均沉降速率以及各時段內的相位形變，并與前人的研究結果和實測水準數據進行對比，證明了本文方法的有效性。通過對研究區域的分析得出，2007～2010年北京地區沉降不均勻較為明顯，并出現多處沉降漏斗，且各沉降漏斗逐漸連成一片并有向東發展的趨勢，多處沉降中心的年平均沉降速率超過100 mm/a；較嚴重的沉降區主要分布在潮白河、溫榆河和泃河等流域的沖積、洪積扇平原上，且部分地區的累計沉降量達到400 mm；地下水的過度開采與北京城市化的擴張嚴重影響地表形變的穩定性，使得地表沉降的幅度與范圍逐漸增大。

[1]Poland M, Bürgmann R, Dzurisin D, et al. Constraints on the Mechanism of Long-Term, Steady Subsidence at Medicine Lake Volcano, Northern California, from GPS, Leveling, and InSAR[J]. Journal of Volcanology and Geothermal Research, 2006, 150(1):55-78

[2]Gabriel A K, Goldstein R M,Zebker H A. Mapping Small Elevation Changes over Large Areas: Differential Radar Interferometry [J].Journal of Geophysical Research, 1989, 94(B7):9 183-9 191

[3]Hooper A. A Multi-Temporal InSAR Method Incorporating Both Persistent Scatterer and Small Baseline Approaches[J]. Geophysical Research Letters, 2008, 35(16):96-106[4]劉媛媛, 張勤, 趙超英,等. PS-InSAR技術用于太原市地面沉降形變分析[J]. 大地測量與地球動力學, 2014, 34(2):6-9(Liu Yuanyuan, Zhang Qin, Zhao Chaoying, et al. Analysis of Ground Subsidence Deformation in Taiyuan City with PS-InSAR Technique[J].Journal of Geodesy and Geodynamics, 2014, 34(2):6-9)

[5]Beradino 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 & Remote Sensing, 2002, 40(11):2 375-2 383

[6]Casu F, Manzo M, Lanari R. A Quantitative Assessment of the SBAS Algorithm Performance for Surface Deformation Retrieval from DInSAR Data[J]. Remote Sensing of Environment, 2006, 102(3-4):195-210

[7]Chaussard E, Wdowinski S, Cabral-Cano E, et al. Land Subsidence in Central Mexico Detected by ALOS InSAR Time-Series[J]. Remote Sensing of Environment, 2014, 140(1):94-106

[8]Dong S C, Samsonov S, Yin H W, et al. Time-Series Analysis of Subsidence Associated with Rapid Urbanization in Shanghai, China Measured with SBAS InSAR Method[J]. Environmental Earth Sciences, 2014, 72(3):677-691

[9]Lanari R, Mora O, Manunta M, et al. A Small-Baseline Approach for Investigating Deformations on Full-Resolution Differential SAR Interferograms[J]. IEEE Transactions on Geoscience & Remote Sensing, 2004, 42(7):1 377-1 386

[10] Lanari R, Casu F, Manzo M, et al. An Overview of the Small Baseline Subset Algorithm: A DInSAR Technique for Surface Deformation Analysis[J].Pure & Applied Geophysics, 2007,164(4):637-661

[11] Usai S. A Least Squares Database Approach for SAR Interferometric Data[J]. IEEE Transactions on Geoscience & Remote Sensing, 2003, 41(4):753-760

[12] Usai S, Klees R. SAR Interferometry on a Very Long Time Scale: A Study of the Interferometric Characteristics of Man-Made Features[J]. IEEE Transactions on Geoscience & Remote Sensing, 1999, 37(4):2 118-2 123

[13] Ng A H M, Ge L, Li X, et al. Monitoring Ground Deformation in Beijing, China with Persistent Scatterer SAR Interferometry[J]. Journal of Geodesy, 2012, 86(6):375-392

Foundation support:National Natural Science Foundation of China,No.41474004,41461089; Fund of Beijing Key Laboratory of Urban Spatial Information Engineering,No. 2016204;Fund of Guangxi Key Laboratory of Spatial Information and Geomatics,No. 15-140-07-32; Open Fund of Key Laboratory of Precise Engineering and Industry Surveying, NASMG,No. PF2013-10.

About the first author:ZHOU Lü, PhD candidate, majors in InSAR data processing, E-mail: zhoulv_whu@163.com.

Monitoring and Analyzing on Ground Settlement in Beijing Area Based on SBAS-InSAR

ZHOULü1,2,3GUOJiming1,4LIXin1HUJiyuan1

1School of Geodesy and Geomatics, Wuhan University, 129 Luoyu Road, Wuhan 430079,China 2Beijing Key Laboratory of Urban Spatial Information Engineering, Beijing Institute of Surveying and Mapping,15 Yangfangdian Road,Beijing 100038,China 3Guangxi Key Laboratory for Spatial Information and Geomatics, Guilin University of Technology,319 Yanshan Street,Guilin 541004,China 4Key Laboratory of Precise Engineering and Industry Surveying, NASMG, Wuhan University,129 Luoyu Road, Wuhan 430079,China

In this paper, SBAS-InSAR is used to obtain high resolution ground subsidence information for the Beijing region. A spatial-temporal analysis of the ground subsidence in the region during the years of 2007 to 2010 is performed utilizing eighteen ENVISAT ASAR images. The results show that subsidence in the Beijing region is severely uneven; that multiple subsidence funnels formed in Changping district, Shunyi district, Tongzhou district, etc. are interconnected and have an eastward expansion trend; that the subsidence velocities in most areas are in the range of -150 mm/a to 10 mm/a and the maximum subsidence is over 400 mm; and that ground subsidence is influenced by groundwater exploitation and urbanization significantly.

ground subsidence; SBAS-InSAR; Beijing area; spatial-temporal analysis

GUO Jiming, professor, PhD supervisor, majors in precise engineering surveying and deformation monitoring, E-mail: jmguo@sgg.whu.edu.cn.

2015-09-18

周呂，博士生，主要從事InSAR數據處理研究，E-mail: zhoulv_whu@163.com。

郭際明，教授，博士生導師，主要從事精密工程測量與形變監測研究，E-mail: jmguo@sgg.whu.edu.cn。

10.14075/j.jgg.2016.09.009

1671-5942(2016)09-0793-05

P237

A

項目來源：國家自然科學基金(41474004，41461089)；城市空間信息工程北京市重點實驗室基金(2016204)；廣西空間信息與測繪重點實驗室基金(15-140-07-32)；精密工程與工業測量國家測繪地理信息局重點實驗室開放基金(PF2013-10)。