基于SBAS-InSAR的北京地區(qū)地表沉降監(jiān)測與分析

周　呂　郭際明　李　昕　胡紀(jì)元

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

?

基于SBAS-InSAR的北京地區(qū)地表沉降監(jiān)測與分析

周呂1,2,3郭際明1,4李昕1胡紀(jì)元1

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

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

地表沉降；SBAS-InSAR；北京地區(qū)；時空分析

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

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

1　SBAS-InSAR技術(shù)原理

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

(1)

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

(2)

式中，φ(tB,x,r)和φ(tA,x,r)分別為tB與tA時刻SAR影像上的相位值，φdef,j(x,r)為tA時刻至tB時刻間衛(wèi)星視線向的形變相位，φ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á)傳播信號的波長，d(tB,x,r)和d(tA,x,r)分別是tB與tA時刻相對于參考時刻t0的雷達(dá)視線向的累積形變量，B⊥j為干涉圖j的垂直基線，Δz為DEM誤差，R為雷達(dá)與目標(biāo)點之間的距離，θ為入射角，φatm,j(tB,x,r)和φatm,j(tA,x,r)分別為tB與tA時刻SAR影像中的大氣相位分量。

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

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

(4)

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

為求解研究區(qū)域各高相干點的形變速率，可用兩景影像間的平均相位速率代替相位值，則式(4)變?yōu)椋?/p>

(5)

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

(6)

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

2　研究區(qū)概況

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

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

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

3.1數(shù)據(jù)處理

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

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

圖2　時空基線結(jié)構(gòu)Fig.2　The spatial-temporal structure of baselines

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

3.2結(jié)果分析

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

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

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

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

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

圖5　水準(zhǔn)測量與SBAS-InSAR結(jié)果對比Fig.5　Comparison between SBAS-InSAR results and leveling results

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

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

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

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

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

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

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

4　結(jié)　語

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

[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技術(shù)用于太原市地面沉降形變分析[J]. 大地測量與地球動力學(xué), 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數(shù)據(jù)處理研究，E-mail: zhoulv_whu@163.com。

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

10.14075/j.jgg.2016.09.009

1671-5942(2016)09-0793-05

P237

A

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