利用SBAS-InSAR技術(shù)監(jiān)測(cè)北京地面沉降

祝 杰 劉洋洋 邵銀星 陶志剛

（中國(guó)北京100045中國(guó)地震臺(tái)網(wǎng)中心）

0 引言

地面沉降是在自然或人為因素作用下，由于地下松散地層固結(jié)壓縮導(dǎo)致地殼表面標(biāo)高緩慢下降的一種局部地質(zhì)現(xiàn)象，具有形成時(shí)間長(zhǎng)、影響范圍廣、防治難度大和過(guò)程不可逆等特點(diǎn)（何慶成等，2006）。研究表明，我國(guó)已形成長(zhǎng)江三角洲、華北平原、汾渭盆地呈三足鼎立的地面嚴(yán)重沉降區(qū)，近百個(gè)城市正在遭受不同程度的地面沉降。北京地區(qū)的地面沉降經(jīng)歷了形成階段（20世紀(jì)50—70年代）、發(fā)展階段（20世紀(jì)70—80年代）、擴(kuò)展階段（20世紀(jì)80—90年代）和快速發(fā)展階段（2000年至今）（周朝棟等，2017），沉降演化過(guò)程呈現(xiàn)范圍不斷擴(kuò)大、變形逐漸加速的趨勢(shì)。

對(duì)地面沉降進(jìn)行大范圍、長(zhǎng)時(shí)間、高精度監(jiān)測(cè)是開(kāi)展沉降控制工作的基礎(chǔ)。基巖標(biāo)、分層標(biāo)、GNSS（Global Navigation Satellite System）、水準(zhǔn)點(diǎn)等監(jiān)測(cè)方式存在空間密度低、無(wú)法提供廣域形變信息等問(wèn)題，難以為防災(zāi)減災(zāi)提供可靠的決策依據(jù)（張永紅等，2016）。合成孔徑雷達(dá)干涉測(cè)量（Interferometric Synthetic Aperture Radar，InSAR）技術(shù)具有明顯的大范圍、高精度、高時(shí)空分辨率、全天候監(jiān)測(cè)等特點(diǎn)，彌補(bǔ)了現(xiàn)有常規(guī)沉陷監(jiān)測(cè)手段的不足，已在地震同震、震間和震后地殼形變監(jiān)測(cè)、地質(zhì)災(zāi)害隱患早期識(shí)別、區(qū)域地面沉降監(jiān)測(cè)等領(lǐng)域（王超等，2002；常占強(qiáng)等，2008；張景發(fā)等，2008；許強(qiáng)等，2019）成功應(yīng)用。在北京沉降研究中，李永生等（2013）利用2個(gè)相鄰軌道的ASAR數(shù)據(jù)提取北京地區(qū)沉降信息，分析了沉降空間分布特征，發(fā)現(xiàn)2003—2010年北京主要沉降區(qū)域快速下沉。張永紅等（2016）對(duì)京津冀地區(qū)1992—2014年3階段地面沉降進(jìn)行了InSAR時(shí)序監(jiān)測(cè)，也發(fā)現(xiàn)北京地區(qū)22年間地面沉降呈不斷加重的趨勢(shì)。Chen等（2016）分別利用ASAR和TerraSAR-X數(shù)據(jù)提取了北京2003—2011年的地表形變，得出最大形變超過(guò)100 mm/a的結(jié)論。上述研究均分析了北京地面沉降演化趨勢(shì)，為沉降防控提供了參考。2014年12月南水北調(diào)中線工程正式運(yùn)營(yíng)，每年向京輸水超1×108m3，改變了供水格局，也為沉降防控創(chuàng)造了有利條件（于海若等，2020；雷坤超，2023）。在北京水情發(fā)生顯著變化的背景下，地下水深度的增加和新的集中開(kāi)采區(qū)的出現(xiàn)，肯定會(huì)造成沉降分布出現(xiàn)新的特征，再加之北京近年地面沉降危險(xiǎn)性評(píng)估工作鮮有研究，因此對(duì)南水北調(diào)工程運(yùn)營(yíng)后北京地區(qū)近年的沉降開(kāi)展時(shí)序監(jiān)測(cè)尤為必要。鑒于此，本文利用SBAS（Small Baseline Subset）時(shí)間序列InSAR技術(shù)對(duì)覆蓋北京地區(qū)的2019年12月24日至2021年8月27日的52景Sentinel-1數(shù)據(jù)進(jìn)行差分干涉處理，提取了北京近期的地面沉降信息，分析了沉降時(shí)空分布和演化特征，并開(kāi)展了地面沉降危險(xiǎn)性評(píng)估，為北京地區(qū)地面沉降防治等工作提供了參考依據(jù)。

1 研究區(qū)概況及數(shù)據(jù)來(lái)源

1.1 研究區(qū)概況

北京地處華北平原西北緣，地理范圍為（115°25′—117°35′E，39°28′—41°05′N(xiāo)）（圖1），總體地勢(shì)西北高、東南低。西部山區(qū)屬太行山山脈，北部山區(qū)屬燕山山脈，東南部平原地區(qū)是由拒馬河、永定河、潮白河、北運(yùn)河、薊運(yùn)河五大水系共同作用形成的沖積扇群構(gòu)成。目前探明的北京地區(qū)發(fā)育的活動(dòng)斷裂帶主要有南口—孫河斷裂、黃莊—高麗營(yíng)斷裂、順義—良鄉(xiāng)斷裂、南苑—通縣斷裂、南口山前斷裂、紫荊關(guān)斷裂、永定河斷裂、大興凸起東源斷裂、夏墊斷裂、二十里長(zhǎng)山斷裂等。

圖1 北京地區(qū)地理位置及SAR影像覆蓋范圍Fig.1 Geographical location of Beijing and SAR image coverage

1.2 數(shù)據(jù)來(lái)源

選取覆蓋北京地區(qū)的52景142軌道的Sentinel-1A TOPS SAR數(shù)據(jù)進(jìn)行干涉處理，時(shí)間跨度為2019年12月24日至2021年8月27日。Sentinel-1A數(shù)據(jù)的空間分辨率是5 m×20 m，C波段降軌VV極化，主要參數(shù)見(jiàn)表1。

表1 SAR數(shù)據(jù)基本參數(shù)Table 1 Basic parameters of SAR data

還獲取了歐空局發(fā)布的POD精密定軌星歷數(shù)據(jù)，該數(shù)據(jù)通過(guò)GNSS下行后處理獲得，數(shù)據(jù)定位精度為5 cm，主要用于消除軌道誤差。選取的數(shù)字高程模型（Digital Elevation Model，DEM）是由美國(guó)航空航天局提供的30 m空間分辨率的SRTM3 DEM數(shù)據(jù)，用于進(jìn)行地理編碼和去除地形相位誤差。Sentinel-1A、POD精密定軌星歷數(shù)據(jù)下載網(wǎng)址：https://scihub.copernicus.eu/dhus/#/home；DEM下載網(wǎng)址：https://dwtkns.com/srtm30m/。

2 SBAS-InSAR處理

小基線集方法（SBAS-InSAR）是Berardino等（2002）提出的一種用于監(jiān)測(cè)地表形變時(shí)間演化的差分干涉測(cè)量算法，基本思想是將SAR數(shù)據(jù)按照一定的規(guī)則進(jìn)行組合，使同一子集內(nèi)干涉像對(duì)的時(shí)間和空間基線較小，而各子集之間的時(shí)間基線和空間基線都較大，得到一系列短空間基線差分干涉圖。通過(guò)集合內(nèi)的最小二乘和集合間的奇異值分解（SVD）法，得到整個(gè)時(shí)間序列地表形變信息的聯(lián)合求解結(jié)果。假設(shè)有N+1幅覆蓋同一地區(qū)的SAR影像，選擇一幅作為主影像，將其他影像配準(zhǔn)，設(shè)定合適的時(shí)間基線閾值和空間基線閾值，將時(shí)間與空間基線都小于該閾值的影像組成若干集合，其中每幅影像至少可與其他N幅影像中的一幅形成干涉像對(duì)。

設(shè)定空間基線閾值為150 m，時(shí)間基線閾值為36 d進(jìn)行短基線干涉影像對(duì)組合，生成的111個(gè)干涉像對(duì)連接如圖2所示，其中空間基線最小1.2 m，最大149.7 m，時(shí)間基線最小12 d，最大36 d。

圖2 Sentinel-1數(shù)據(jù)干涉對(duì)時(shí)空基線分布Fig.2 Spatiotemporal baseline distribution of Sentinel-1 data interferograms

對(duì)所有干涉像對(duì)進(jìn)行干涉處理，距離向方位向多視比設(shè)置為10:2，以保證數(shù)據(jù)解算中輸出的相位圖保持高分辨率。利用歐空局提供的POD精密定軌星歷數(shù)據(jù)進(jìn)行SLC（Single-Look Complex)數(shù)據(jù)軌道校正，去除軌道誤差。用SRTM3 DEM數(shù)據(jù)去除地形相位，得到111個(gè)差分干涉圖。采用Goldstein自適應(yīng)濾波方法對(duì)差分干涉圖進(jìn)行濾波，減弱大氣和各類噪聲誤差，采用最小費(fèi)用流算法（Eineder et al，1998）進(jìn)行相位解纏，相位解纏中相干性閾值設(shè)置為0.4，對(duì)于相干性閾值以下的區(qū)域不進(jìn)行解纏計(jì)算，采用高低通濾波方法進(jìn)一步消除大氣延遲相位影響，得到差分解纏相位（Gabriel et al，1989；Massonnet et al，1993；Wright et al，2004）。檢查所有相位解纏圖質(zhì)量，剔除相位解纏錯(cuò)誤和質(zhì)量較差的干涉像對(duì)。然后基于解纏的相位建立觀測(cè)方程，采用奇異值分解方法（李永生等，2013）求解形變速率的最小二乘范數(shù)解，對(duì)求取的形變速率在時(shí)域上進(jìn)行積分，得到形變時(shí)間序列結(jié)果。

3 北京地區(qū)InSAR時(shí)間序列結(jié)果及空間演化特征分析

3.1 北京地區(qū)InSAR地表形變監(jiān)測(cè)

運(yùn)用SBAS-InSAR技術(shù)提取北京地區(qū)2019年12月24日至2021年8月27日的地表形變信息。52景Sentinel-1A SAR數(shù)據(jù)提取的地表沿雷達(dá)視線方向的平均形變速率見(jiàn)圖3，可以看出，提取的高相干點(diǎn)位空間分布密度高，典型形變區(qū)域的高相干點(diǎn)位呈現(xiàn)清晰、連續(xù)的分布狀態(tài)，未發(fā)生明顯空洞現(xiàn)象。需說(shuō)明，InSAR技術(shù)提取到懷柔和密云部分地區(qū)存在微抬升現(xiàn)象，初步分析可能是數(shù)據(jù)處理過(guò)程中邊緣效應(yīng)導(dǎo)致的，實(shí)際并未發(fā)生。

3.2 空間分布特征分析

從2019年12月24日至2021年8月27日的地表平均形變速率可以看出，近垂向形變呈現(xiàn)明顯的不均勻、不規(guī)則、漏斗式運(yùn)動(dòng)特征，沉降現(xiàn)象在朝陽(yáng)、通州、大興和海淀4個(gè)區(qū)比較顯著，出現(xiàn)4個(gè)明顯的片狀沉降區(qū)域：通州沉降區(qū)[圖4（a）]、朝陽(yáng)金盞沉降漏區(qū)[圖4（b）]、海淀稻香湖公園沉降區(qū)[圖4（c）]和大興機(jī)場(chǎng)周邊沉降區(qū)[圖4（d）]。監(jiān)測(cè)到最大地表形變速率位于北京大興安定垃圾清埋場(chǎng)，約-156 mm/a，沉降面積約88×104m2，該區(qū)域沉降暫未展現(xiàn)進(jìn)一步向外侵蝕的趨勢(shì)。下面具體分析4個(gè)沉降區(qū)的地表形變。

圖4 主要沉降漏斗區(qū)空間分布特征Fig.4 Spatial distribution characteristics of the main subsidence funnel zones

3.2.1 通州沉降區(qū)。通州沉降區(qū)主要由通州次渠鎮(zhèn)、亦莊火車(chē)站、北京環(huán)球影城、北馬頭村4個(gè)離散的沉降漏斗組成，呈現(xiàn)近NE向展布，沉降總面積達(dá)697×104m2，最大形變速率約-120 mm/a。其中，次渠鎮(zhèn)漏斗約186×104m2，平均形變速率約-25 mm/a。亦莊火車(chē)站站前街道的京津高速桂家墳橋—麥莊橋周邊區(qū)域沉降面積約60×104m2，最大形變速率達(dá)-120 mm/a。通州沉降區(qū)域中北京環(huán)球影城周邊地區(qū)的形變面積最大，約371×104m2，最大形變速率達(dá)-30 mm/a。北馬頭村周邊沉降區(qū)域面積約80×104m2，最大形變速率達(dá)-37 mm/a。

3.2.2 朝陽(yáng)金盞沉降區(qū)。該沉降區(qū)由朝陽(yáng)區(qū)金盞鄉(xiāng)、黎各莊鄉(xiāng)、馬各莊環(huán)繞的不連續(xù)沉降漏斗組成，呈扇形分布，總沉降面積約4500×104m2，最大地表形變速率達(dá)-110 mm/a。其中，金盞鄉(xiāng)沉降區(qū)主要由金盞鄉(xiāng)、雷莊村和黑橋村3個(gè)子沉降漏斗組成，總沉降面積約1578×104m2，金盞鄉(xiāng)、雷莊村和黑橋村3個(gè)子沉降區(qū)的最大形變速率分別約為-27 mm/a、-40 mm/a、-34 mm/a，雷莊村演播大樓附近地表形變速率最大，約-50 mm/a。黎各莊鄉(xiāng)附近沉降范圍涉及整個(gè)宋莊水井公園，面積約210×104m2，最大形變速率達(dá)-110 mm/a。馬各莊周邊地面沉降呈離散點(diǎn)狀分布，面積約144×104m2，位于馬各莊高安屯附近地表形變速率最大，達(dá)-82 mm/a，將形變速率結(jié)果疊加到谷歌地圖衛(wèi)星影像上進(jìn)行目視解譯推斷，結(jié)果顯示，馬各莊周邊地面沉降現(xiàn)象大概率由區(qū)域工程活動(dòng)引起。

3.2.3 海淀稻香湖公園沉降區(qū)。海淀稻香湖公園沉降區(qū)位于海淀西北部，由蘇家坨鎮(zhèn)溫陽(yáng)路、稻香湖公園和南沙河3個(gè)子沉降區(qū)組成，呈三角形分布。此區(qū)域沉降面積較小，約114×104m2，最大形變速率約-30 mm/a，根據(jù)形變分布可以看出，稻香湖公園和南沙河子沉降區(qū)有進(jìn)一步相向擴(kuò)張的運(yùn)動(dòng)趨勢(shì)。海淀區(qū)沉降速率最大的區(qū)域位于六里屯垃圾清埋場(chǎng)，達(dá)-152 mm/a，但該區(qū)域沉降現(xiàn)象僅發(fā)生在內(nèi)部，暫未呈現(xiàn)進(jìn)一步向外擴(kuò)張的趨勢(shì)。

3.2.4 大興機(jī)場(chǎng)周邊沉降區(qū)。大興機(jī)場(chǎng)周邊沉降區(qū)由北京大興機(jī)場(chǎng)周邊的20余個(gè)小范圍沉降區(qū)組成，離散分布在大興機(jī)場(chǎng)北側(cè)、西側(cè)和南側(cè)，在機(jī)場(chǎng)東側(cè)未見(jiàn)明顯沉降現(xiàn)象。總沉降面積約420×104m2，最大形變速率約-80 mm/a。最近的義和場(chǎng)村沉降漏斗距機(jī)場(chǎng)東側(cè)飛機(jī)跑道僅約300 m，因此北京大興機(jī)場(chǎng)周邊區(qū)域的地面沉降態(tài)勢(shì)有必要重點(diǎn)關(guān)注。

4 北京地區(qū)地面沉降時(shí)間演化特征分析

選取6個(gè)特征點(diǎn)（圖5）的形變量時(shí)間序列來(lái)表征北京典型區(qū)域地面沉降隨時(shí)間的演化過(guò)程。特征點(diǎn)選取依據(jù)為：優(yōu)先選擇位于局部沉降速率較大的漏斗區(qū)且演化過(guò)程中未出現(xiàn)明顯突跳的點(diǎn)。這種依據(jù)能使選點(diǎn)基本符合局部區(qū)域沉降的共性特征。下面具體分析6個(gè)特征點(diǎn)形變量時(shí)間序列。

圖5 6個(gè)特征點(diǎn)形變量時(shí)間序列（a）首都機(jī)場(chǎng)跑道特征點(diǎn)1；（b）首都機(jī)場(chǎng)跑道特征點(diǎn)2；（c）金盞地區(qū)特征點(diǎn)；（d）故宮特征點(diǎn)；（e）高安屯特征點(diǎn)；（f）大興機(jī)場(chǎng)特征點(diǎn)Fig.5 Time-series of deformation quantity for six feature points

（1）從首都機(jī)場(chǎng)跑道上2個(gè)特征點(diǎn)的形變量變化過(guò)程可以看出，首都機(jī)場(chǎng)部分飛機(jī)跑道正在遭受著地面沉降現(xiàn)象，2點(diǎn)位置在研究期內(nèi)的平均形變量約-31 mm，形變量基本呈線性發(fā)展趨勢(shì)，沉降速率分別是-19 mm/a和-22 mm/a，因此可知首都機(jī)場(chǎng)雖正在遭受地面沉降，但整體形變速率較小且較穩(wěn)定，影響不大。

（2）金盞地區(qū)的地面沉降也基本呈線性發(fā)展趨勢(shì)，累計(jì)形變量約為-166 mm，形變速率約-127 mm/a，形變速率較大，地面沉降形勢(shì)較嚴(yán)重。

（3）故宮在研究期內(nèi)垂向形變基本保持為零，因此故宮未發(fā)生地面沉降現(xiàn)象。

（4）高安屯地區(qū)也是一個(gè)明顯的沉降漏斗，累計(jì)形變量持續(xù)加大，約-167 mm，平均形變速率達(dá)-121 mm/a。

（5）在研究期內(nèi)大興機(jī)場(chǎng)周邊出現(xiàn)20余個(gè)離散且面積很小的沉降漏斗，選擇其中一個(gè)漏斗分析其沉降量的變化，可以看出該特征點(diǎn)形變也呈線性趨勢(shì)，累計(jì)形變量約-122 mm，形變速率約為-71 mm/a。

從6個(gè)特征點(diǎn)的形變量時(shí)間序列及地表形變空間分布可以看出，北京朝陽(yáng)金盞地區(qū)和高安屯地區(qū)的地表形變速率較大，均超過(guò)-120 mm/a，值得進(jìn)一步關(guān)注。

5 北京地區(qū)地面沉降危險(xiǎn)性淺析

在地表平均形變速率基礎(chǔ)上，根據(jù)北京地面沉降災(zāi)害危險(xiǎn)性分區(qū)量化指標(biāo)（表2）（楊艷等，2010；楊艷，2015）對(duì)北京地區(qū)地面沉降危險(xiǎn)性進(jìn)行初步分析。北京地面沉降危險(xiǎn)性評(píng)價(jià)結(jié)果如圖6所示。

表2 北京地面沉降危險(xiǎn)性評(píng)價(jià)形變速率分區(qū)Table 2 Risk assessment of land subsidence in Beijing based on deformation rate zoning

圖6 北京地表形變速率危險(xiǎn)性分區(qū)Fig.6 Risk zoning of surface deformation rate in Beijing

從北京地區(qū)地面沉降危險(xiǎn)性分區(qū)結(jié)果（表3）可知，在2019年12月24日至2021年8月27日期間，北京地區(qū)嚴(yán)重沉降區(qū)占比較小，僅為0.04%，主要位于海淀區(qū)六里屯垃圾清埋場(chǎng)、大興區(qū)安定垃圾清理場(chǎng)和宋莊水井公園；較嚴(yán)重沉降區(qū)和一般沉降區(qū)占據(jù)的沉降面積比例約為6%，主要位于前述的通州沉降區(qū)、朝陽(yáng)金盞沉降區(qū)、海淀稻香湖公園沉降區(qū)和大興機(jī)場(chǎng)周邊沉降區(qū)。這些區(qū)域的沉降速率均超過(guò)-30 mm/a，屬于地面沉降的高危險(xiǎn)區(qū)，進(jìn)一步發(fā)育可能會(huì)對(duì)城市基礎(chǔ)設(shè)施、重要建筑物造成不可挽回的損害。輕微沉降區(qū)和較穩(wěn)定區(qū)占比分別為46.65%和47.15%，二者之和達(dá)93.80%，可認(rèn)為在研究期內(nèi)北京絕大部分地區(qū)未發(fā)生明顯的地面沉降現(xiàn)象。

6 結(jié)論

在南水北調(diào)新水情下，北京地區(qū)地面沉降發(fā)生明顯變化，為揭示地面沉降分布及其演化特征，利用SBAS-InSAR技術(shù)成功提取了北京地區(qū)2019年12月24日至2021年8月27日地面沉降信息，進(jìn)而分析地面沉降時(shí)空分布特征和演化規(guī)律，開(kāi)展地面沉降危險(xiǎn)性評(píng)估分析。主要結(jié)論如下：

（1）2019年12月至2021年8月北京地區(qū)出現(xiàn)地面沉降現(xiàn)象，沉降分布在空間上呈現(xiàn)明顯的不均勻、不規(guī)則、漏斗式運(yùn)動(dòng)特征。但總體上沉降面積占比較小，最大形變速率約-156 mm/a。

（2）地面沉降現(xiàn)象在朝陽(yáng)、通州、大興、海淀4個(gè)區(qū)比較顯著，出現(xiàn)了4個(gè)明顯的片狀沉降區(qū)域，分別是通州沉降區(qū)、朝陽(yáng)金盞沉降區(qū)、海淀稻香湖公園沉降區(qū)和大興機(jī)場(chǎng)周邊沉降區(qū)。通州沉降區(qū)4個(gè)沉降漏斗的地面沉降面積達(dá)697×104m2，最大平均沉降速率約-120 mm/a。金盞沉降區(qū)總沉降面積約4500×104m2，最大地表形變速率達(dá)-110 mm/a。海淀稻香湖公園沉降區(qū)沉降面積較小為114×104m2，最大形變速率約-30 mm/a。大興機(jī)場(chǎng)周邊沉降區(qū)沉降面積約420×104m2，最大形變速率約-80 mm/a。

（3）從典型形變點(diǎn)的形變量時(shí)間序列可以看出，北京朝陽(yáng)金盞和高安屯地區(qū)形變量基本呈線性發(fā)展趨勢(shì)，地表形變速率較大，均超過(guò)-120 mm/a，值得進(jìn)一步關(guān)注。

（4）地面沉降速率分類結(jié)果表明，嚴(yán)重沉降區(qū)、較嚴(yán)重和一般沉降區(qū)占北京總面積的6.2%，主要位于上述4個(gè)地面沉降區(qū)，危險(xiǎn)程度較高，值得重點(diǎn)關(guān)注，其他區(qū)域地面沉降危險(xiǎn)性較低。