祝 杰 劉洋洋 邵銀星 陶志剛

（中國北京100045中國地震臺網(wǎng)中心）

0 引言

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

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

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

1.1 研究區(qū)概況

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

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

1.2 數(shù)據(jù)來源

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

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

還獲取了歐空局發(fā)布的POD精密定軌星歷數(shù)據(jù)，該數(shù)據(jù)通過GNSS下行后處理獲得，數(shù)據(jù)定位精度為5 cm，主要用于消除軌道誤差。選取的數(shù)字高程模型（Digital Elevation Model，DEM）是由美國航空航天局提供的30 m空間分辨率的SRTM3 DEM數(shù)據(jù)，用于進行地理編碼和去除地形相位誤差。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)測地表形變時間演化的差分干涉測量算法，基本思想是將SAR數(shù)據(jù)按照一定的規(guī)則進行組合，使同一子集內(nèi)干涉像對的時間和空間基線較小，而各子集之間的時間基線和空間基線都較大，得到一系列短空間基線差分干涉圖。通過集合內(nèi)的最小二乘和集合間的奇異值分解（SVD）法，得到整個時間序列地表形變信息的聯(lián)合求解結果。假設有N+1幅覆蓋同一地區(qū)的SAR影像，選擇一幅作為主影像，將其他影像配準，設定合適的時間基線閾值和空間基線閾值，將時間與空間基線都小于該閾值的影像組成若干集合，其中每幅影像至少可與其他N幅影像中的一幅形成干涉像對。

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

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

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

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

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

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

3.2 空間分布特征分析

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

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

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

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

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

3.2.4 大興機場周邊沉降區(qū)。大興機場周邊沉降區(qū)由北京大興機場周邊的20余個小范圍沉降區(qū)組成，離散分布在大興機場北側、西側和南側，在機場東側未見明顯沉降現(xiàn)象?？偝两得娣e約420×104m2，最大形變速率約-80 mm/a。最近的義和場村沉降漏斗距機場東側飛機跑道僅約300 m，因此北京大興機場周邊區(qū)域的地面沉降態(tài)勢有必要重點關注。

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

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

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

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

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

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

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

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

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

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

在地表平均形變速率基礎上，根據(jù)北京地面沉降災害危險性分區(qū)量化指標（表2）（楊艷等，2010；楊艷，2015）對北京地區(qū)地面沉降危險性進行初步分析。北京地面沉降危險性評價結果如圖6所示。

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

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

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

6 結論

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

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

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

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

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