劉金昌，王慧妮，徐 吟，李康倫，張華睿，劉天鶴

1. 武漢工程大學(xué)土木工程與建筑學(xué)院，湖北 武漢 430074；

2. 中鐵十一局集團(tuán)有限公司，湖北 武漢 430061；

3. 武漢綜合交通研究院有限公司，湖北 武漢 430015

在大規(guī)模的地鐵施工建設(shè)過程中，由于受到地質(zhì)環(huán)境、周邊環(huán)境、施工方法和方案等因素的影響，由城市地鐵施工引起的地面塌陷、突涌水［1］等災(zāi)害屢見不鮮。地面塌陷是指在人為因素或者自然因素的作用下，地表巖、土體向下陷落，并形成塌陷坑（洞）的一種動(dòng)力地質(zhì)現(xiàn)象［2］。地面塌陷是突發(fā)性的，尤其是在城市建筑物密集的地方，一旦發(fā)生塌陷事故將會(huì)造成特別嚴(yán)重的經(jīng)濟(jì)損失和生命安全問題。以青島地鐵施工為例，2019 年5 月27 日，青島地鐵4 號(hào)線沙子口靜沙區(qū)間施工段發(fā)生坍塌，同年7 月4 日，青島地鐵1 號(hào)線勝利橋站施工圍擋處再次發(fā)生塌陷。

在城市地鐵施工區(qū)域開展地面塌陷的監(jiān)測(cè)和早期預(yù)警研究具有重要意義。在地面塌陷災(zāi)害形成過程中，地表早期會(huì)出現(xiàn)緩慢下沉、開裂等異?，F(xiàn)象［3］，可以利用塌陷災(zāi)變演化過程特性基于長(zhǎng)時(shí)序大范圍地表形變開展監(jiān)測(cè)，以實(shí)現(xiàn)潛在的地面塌陷災(zāi)害隱患點(diǎn)的早期識(shí)別及預(yù)警。合成孔徑雷達(dá)干涉測(cè)量（interferometric synthetic aperture radar，InSAR）技術(shù)是利用合成孔徑雷達(dá)（synthetic aperture radar，SAR）對(duì)同一地區(qū)觀測(cè)到的2 幅復(fù)數(shù)值影像（既有幅值又有相位的影像）數(shù)據(jù)進(jìn)行相干處理，以獲取地表高程信息的新技術(shù)，因其具有全天候、廣覆蓋、高精度、低成本、低風(fēng)險(xiǎn)獲取的特點(diǎn)［4-7］，現(xiàn)如今已經(jīng)廣泛應(yīng)用于冰川融化、地質(zhì)災(zāi)害監(jiān)測(cè)、城市地下水變化和地面塌陷等沉降的監(jiān)測(cè)中［8-11］。InSAR 技術(shù)的測(cè)量范圍大，形變測(cè)量精度可達(dá)到mm 級(jí)別，非常適合對(duì)城市地表的形變監(jiān)測(cè)研究。目前，國(guó)內(nèi)外學(xué)者利用InSAR 技術(shù)開展的地面沉降方面研究工作主要集中于地鐵沿線的沉降監(jiān)測(cè)［12-16］，但基于長(zhǎng)時(shí)序地表形變監(jiān)測(cè)數(shù)據(jù)結(jié)合歷史塌陷事故，開展地面塌陷預(yù)警的研究工作相對(duì)較少。

本文以青島地鐵工程施工過程中的2 次地鐵塌陷事故為研究對(duì)象，收集2018 年8 月至2020 年8 月期間青島市的60 景哨兵1 號(hào)（Sentinel-1A）衛(wèi)星數(shù)據(jù)，基于差分干涉測(cè)量短基線集時(shí)序分析（small baseline subset-InSAR，SBAS-InSAR）技術(shù)獲取塌陷事故前后期地表的時(shí)序形變過程數(shù)據(jù)，結(jié)合歷史塌陷事故開展地表時(shí)序形變與塌陷災(zāi)害演化過程相關(guān)性分析，建立利用SAR 衛(wèi)星遙感數(shù)據(jù)開展城市地鐵施工地面塌陷早期識(shí)別與監(jiān)測(cè)預(yù)警方法，為城市地鐵施工監(jiān)測(cè)預(yù)警提供技術(shù)支撐。

1 基本情況

青島地處山東半島東南部，位于東經(jīng)119°30′~121°00′、北緯35°35′~37°09′。青島所處大地構(gòu)造位置為新華夏隆起帶次級(jí)構(gòu)造單元——膠南隆起區(qū)東北緣和膠萊凹陷區(qū)中南部。區(qū)內(nèi)缺失整個(gè)古生界地層及部分中生界地層，但白堊系青山組火山巖層發(fā)育充分，在青島市出露十分廣泛。巖漿巖以元古代膠南期月季山式片麻狀花崗巖及中生代燕山晚期的艾山式花崗閃長(zhǎng)巖和嶗山式花崗巖為主。市區(qū)全部坐落于該類花崗巖之上，建筑地基條件優(yōu)良。構(gòu)造以斷裂構(gòu)造為主。

沙子口事故發(fā)生于2019 年5 月27 日，青島地鐵4 號(hào)線沙子口靜沙區(qū)間施工段發(fā)生坍塌。勝利橋事故發(fā)生于2019 年7 月4 日，青島地鐵1 號(hào)線勝利橋站施工圍擋處發(fā)生直徑約10 m 的塌陷。

2 實(shí)驗(yàn)數(shù)據(jù)和方法

2.1 SAR 數(shù)據(jù)

本研究選擇了地球觀測(cè)衛(wèi)星Sentinel-1 的遙感數(shù)據(jù)。該衛(wèi)星配備了C 波段SAR，由2 顆衛(wèi)星A，B 組成，衛(wèi)星重訪期為12 d。本實(shí)驗(yàn)搜集了覆蓋研究區(qū)2018-08－2020-08 期間的60 景Sentinel-1A 數(shù)據(jù)，所選數(shù)據(jù)極化方式均為VV（垂直極化發(fā)射，垂直極化接收）極化，所處波段均為C 波段，軌道方向?yàn)樯?，方位向和距離向分辨率分別為20 m 和5 m。研究區(qū)域使用的外部數(shù)字高程模型（digital elevation model，DEM）數(shù)據(jù)為德國(guó)航空航天中心使用TerraSAR-X/TanDEM-X 衛(wèi)星獲取的全球陸地高程數(shù)據(jù)，空間分辨率為30 m。為了提高Sentinel-1 影像的軌道精度，采用了歐空局發(fā)布的精確軌道數(shù)據(jù)，將研究區(qū)域數(shù)據(jù)按照小基線原則（只保留平均相干性大于0.4 的干涉對(duì)，時(shí)間基線設(shè)置為12 d/24 d，空間基線由SBAS-InSAR 方法進(jìn)行計(jì)算）組合而成干涉對(duì)的時(shí)空基線圖，如圖1 所示。

圖1 時(shí)空基線圖Fig.1 Space-time baseline map

2.2 SBAS-InSAR 數(shù)據(jù)處理方法

SBAS-InSAR 是一種新開發(fā)的星載合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)（differential intereferometric synthetic aperture radar，D-InSAR）時(shí)間序列分析方法，它將采集到的所有SAR 數(shù)據(jù)組合成干涉對(duì)，然后設(shè)置空間基線和時(shí)間基線閾值，并獲得幾個(gè)時(shí)空基線小于閾值的小基線集。每個(gè)集合的表面變形時(shí)間序列可以通過奇異值分解（singular value decomposition，SVD）求解集合來解決。這有效地解決了由于基線時(shí)間長(zhǎng)不連續(xù)性概率引起的基線組間差異小基線，提高了監(jiān)測(cè)時(shí)間分辨率，有效解決了一般方程秩缺陷問題，最終得到覆蓋整個(gè)地面沉降觀測(cè)時(shí)間的序列，SBAS-InSAR 技術(shù)的數(shù)據(jù)處理流程如圖2 所示。

圖2 SBAS-InSAR 技術(shù)流程圖Fig.2 SBAS-InSAR technical flowchart

（1）將監(jiān)測(cè)區(qū)域的M+1 景SAR 影像數(shù)據(jù)集根據(jù)給定的干涉組合條件，生成N景小基線距的干涉圖，如圖3（a）所示，將干涉圖中的地形相位去除生成去地形干涉圖。

圖3 SBAS-InSAR 干涉及相位解纏過程：（a）干涉，（b）掩膜，（c）相位解纏，（d）去除軌道傾斜Fig.3 SBAS-InSAR interferogram and phase unwrapping process：（a）interference，（b）mask，（c）phase unwrapping，（d）removal of orbital tilt

（2）以相干系數(shù)為基礎(chǔ)生成掩膜圖，如圖3（b）所示，在去地形干涉圖中選擇高相干點(diǎn)進(jìn)行相位解纏處理，如圖3（c）所示，并對(duì)相干點(diǎn)進(jìn)行定標(biāo)處理。

數(shù)據(jù)處理流程中生成的干涉圖、去除軌道傾斜改正、相位解纏、掩膜如圖3 所示。

（3）根據(jù)不同時(shí)段的平均形變速率、高程誤差、差分相位進(jìn)行模型方程組的建立，利用SVD 算法獲得最小范數(shù)條件下的最小二乘解。在開展SAR 數(shù)據(jù)時(shí)序分析前，使用斜率格網(wǎng)法估計(jì)軌道效應(yīng)，進(jìn)行軌道誤差校正，如圖3（d）去除軌道傾斜所示。使用最小二乘法在每個(gè)干涉圖上獨(dú)立估計(jì)軌道誤差。默認(rèn)情況下，對(duì)于由指數(shù)為i和j的2 個(gè)SAR 采集組成的干涉圖，采用最小范數(shù)反演算法解決中dij(x,y)的值是方位x，y處的像素值。

其中aij、bij、cij和eij稱為干涉圖ij的軌道參數(shù)。

當(dāng)像素格網(wǎng)為1 時(shí)，

當(dāng)像素格網(wǎng)為3 時(shí)，

當(dāng)像素格網(wǎng)為4 時(shí)，

為了確保干涉測(cè)量網(wǎng)絡(luò)的一致性，在網(wǎng)絡(luò)感知中重新估計(jì)軌道參數(shù)。通過i反轉(zhuǎn)線性系統(tǒng)來估計(jì)每個(gè)SAR 采集的軌道參數(shù)。

結(jié)合重新估計(jì)的軌道參數(shù)，生成與干涉測(cè)量網(wǎng)絡(luò)一致的軌道校正圖，并校正每個(gè)干涉圖。

假設(shè)使用覆蓋同一研究區(qū)域的N+1 幅SAR 影像圖，則按時(shí)間順序?qū)λ鼈冞M(jìn)行排序：

N個(gè)干涉對(duì)是通過最大空間基線閾值獲得的M對(duì)多視差分干涉對(duì)生成的，M滿足：

對(duì)生成的多視差分干涉圖進(jìn)行濾波并解開相位，使用相干性在未包覆的差分干涉圖中選擇適當(dāng)?shù)膮⒖键c(diǎn)。對(duì)于任意2 個(gè)微分干涉圖（ta

式（11）中：?x,i是i圖像和x像素的干涉相位，?x,i(tb)和?x,i(ta)是x像素2 個(gè)不同時(shí)間的干涉相位。?flat,x,i是平坦相位的相位差，?ele,x,i是地形相位的相位差，?disp,x,i是地面變形的相位差，?atm,x,i是大氣相位的誤差，?noise是噪聲的誤差。

公式中的平坦相位、地形相位、變形相位、大氣相位可以進(jìn)一步表示為：

式（12）~式（13）中：λ是雷達(dá)波長(zhǎng)，R是斜率距離，Bn是垂直基線，Δq是DEM 高程差，θ是入射角。

假設(shè)2 個(gè)不同干涉圖之間的變形速率為Vi(i= 1,2,3,…,n)，則ta和tb之間的累積變形可以表示為：

變形速率的積分在干涉圖的主圖像和從圖像之間的時(shí)間間隔內(nèi)。

直接在主圖像和從圖像間隔之間的每個(gè)時(shí)間間隔處進(jìn)行速率積分。以矩陣形式書寫：

通過使用SVD方法，可以獲得系數(shù)矩陣B的廣義逆矩陣，從而得到數(shù)據(jù)向量的最小范數(shù)解。最后，每個(gè)時(shí)間段的速度決定了每個(gè)時(shí)間段的變形量。

（4）開展估算大氣相位與非線性形變。對(duì)于給定數(shù)據(jù)集，選擇與SAR 場(chǎng)景的空間覆蓋范圍重疊的格網(wǎng)點(diǎn)。大氣變量以精確的壓力水平提供。將這些值垂直插值到表面和參考高程之間的規(guī)則網(wǎng)格中，計(jì)算每個(gè)選定網(wǎng)格點(diǎn)上的延遲函數(shù)作為高程的函數(shù)，高程處的視線方向（line of sight，LOS）單路徑延遲

式（18）中：θ為入射角；Rd＝287.05 J·kg-1·K-1，Rν＝461.495 J·kg-1·K-1，分別是干燥和濕潤(rùn)環(huán)境下的氣體常數(shù)；gm是z和zref之間重力加速度的加權(quán)平均值；P是干燥空氣分壓（Pa）；e是水蒸氣分壓（Pa）；T是溫度（K）。式（18）中：k1＝0.776 K·Pa-1，k2＝0.716 K·Pa-1，k3＝3.75×103K2·Pa-1。

3 結(jié)果與討論

根據(jù)上述數(shù)據(jù)處理流程，使用改進(jìn)的SBASInSAR 時(shí)序分析方法獲取了青島市范圍內(nèi)2018 年8 月9 日至2020 年8 月22 日期間的地表累計(jì)形變計(jì)算結(jié)果，將結(jié)果疊加在天地圖影像上，如圖4所示。

圖4 地表累計(jì)形變圖：（a）沙子口事故區(qū)域，（b）勝利橋事故區(qū)域Fig.4 Cumulative surface deformation maps：（a）Shazikou accident area，（b）Shengliqiao accident area

以事故點(diǎn)為中心，選取方圓300 m 范圍內(nèi)的時(shí)序形變點(diǎn)開展數(shù)理統(tǒng)計(jì)分析。分別在2 個(gè)事故點(diǎn)提取分析多個(gè)點(diǎn)位的時(shí)序形變信息，如圖5 所示。

圖5 事故區(qū)域部分點(diǎn)位累計(jì)沉降：（a）沙子口，（b）勝利橋Fig.5 Cumulative subsidence of some points in accident areas：（a）Shazikou，（b）Shengliqiao

將InSAR 數(shù)據(jù)獲取的地表變形量數(shù)據(jù)轉(zhuǎn)化為日均形變速率，處理后得到該點(diǎn)的日均形變速率，單位為mm/d，如圖6 所示。形變速率計(jì)算與分析過程如下：

圖6 沙子口和勝利橋研究區(qū)域日均沉降速率Fig.6 Average daily subsidence rates of Shazikou and Shengliqiao study areas

（1）將時(shí)間信息和沉降信息一一對(duì)應(yīng)：D0-S0，D1-S1，D2-S2，D3-S3，…，Dn-Sn；

（2）計(jì)算日均地表形變速率

式（19）中：Di（i=1，2，3，…，n）是日期，Si（i=1，2，3…，n）是沉降量。因?yàn)橛玫氖乔? 個(gè)周期的日期減去后1 個(gè)周期的日期，所以地面發(fā)生沉降時(shí)，沉降速率為正值。

基于SBAS-InSAR 時(shí)序形變觀測(cè)信息對(duì)沙子口研究區(qū)開展時(shí)序分析，如圖6 沙子口事故所示，2018 年8 月21 日，該區(qū)域便開始呈現(xiàn)下沉趨勢(shì)，下沉速率在0.011 4 mm/d 左右，有地表沉降形變但速率相對(duì)較為緩慢。2018 年8 月21 日至2019 年1 月12 日該研究區(qū)沉降速率呈現(xiàn)加速增長(zhǎng)趨勢(shì)，于2019 年1 月12 日沉降速率突破0.02 mm/d，在這一段時(shí)間內(nèi)該研究區(qū)累計(jì)沉降量為2.37 mm。從2019 年1 月12 日至2019 年5 月27 日發(fā)生地面塌陷事故，在連續(xù)132 d 內(nèi)的11 次觀測(cè)（Sentinel-1 衛(wèi)星的重訪周期為12 d）中，研究區(qū)日均沉降速率均大于0.02 mm/d，期間地面累計(jì)沉降量為3.29 mm。在事故發(fā)生后，地表形變沉降速率顯著下降，地表形變逐漸趨于穩(wěn)定。

在勝利橋區(qū)域，如圖6 勝利橋事故所示，從2018 年8 月21 日開始，該區(qū)域沉降速率處于-0.010 5 mm/d 左右，這是由于青島市所處板塊的地殼運(yùn)動(dòng)，呈現(xiàn)緩慢的上升趨勢(shì)，但是到2018 年10 月20 日開始出現(xiàn)沉降，沉降速率開始快速增大。出現(xiàn)沉降后，該研究點(diǎn)的沉降速率越來越快，于2019 年5 月12 日沉降速率突破0.02 mm/d，此時(shí)累計(jì)沉降量為1.41 mm。從2019 年5 月12 日到2019 年7 月4 日發(fā)生地面塌陷事故，連續(xù)4 次衛(wèi)星觀測(cè)（53 d）的日均沉降速率均大于0.02 mm/d，期間地面累計(jì)沉降量為1.44 mm。勝利橋區(qū)域在日均沉降速率最大（0.024 5 mm/d）處發(fā)生塌陷事故，事故發(fā)生后，該區(qū)域沉降速率逐漸降低，于2020 年1 月7 日沉降速率開始小于0 mm/d，期間累計(jì)沉降量為4.88 mm。最后于2020 年3 月開始趨于穩(wěn)定的上升狀態(tài)。

實(shí)驗(yàn)數(shù)據(jù)說明：日均沉降速率大于0.02 mm/d可作為監(jiān)測(cè)預(yù)警依據(jù)，若在施工期間，連續(xù)在3 個(gè)監(jiān)測(cè)周期（36 d）出現(xiàn)日均沉降速率大于0.02 mm/d時(shí)，應(yīng)對(duì)施工單位發(fā)出塌陷事故預(yù)警，提醒施工單位此地存在塌陷災(zāi)害隱患，及時(shí)采取更加精密的監(jiān)測(cè)手段來監(jiān)測(cè)地面沉降情況并作出事故預(yù)防措施。

4 結(jié) 論

通過獲取2018 年8 月—2020 年8 月期間青島市的60 景Sentinel-1A 數(shù)據(jù)，基于改進(jìn)的SBASInSAR 方法，計(jì)算長(zhǎng)時(shí)序、高精度地表形變數(shù)據(jù)，通過與青島市2 次事故對(duì)比分析，提出了地鐵施工塌陷監(jiān)測(cè)早期識(shí)別方法和監(jiān)測(cè)預(yù)警指標(biāo)。主要結(jié)論如下：

（1）本文中通過改進(jìn)的SBAS-InSAR 方法，利用SVD 算法獲得最小范數(shù)條件下的最小二乘解，使用斜率格網(wǎng)法估計(jì)軌道效應(yīng)，并進(jìn)行軌道誤差校正，同時(shí)結(jié)合高程、氣溫、氣壓及濕度估計(jì)值構(gòu)建估算函數(shù)，使用大氣變量的3D 分布來確定每個(gè)干涉圖的每個(gè)像素上的大氣相位延遲，以提高SBAS-InSAR 方法處理的地表形變數(shù)據(jù)精度。實(shí)驗(yàn)表明該方法獲取的地表監(jiān)測(cè)數(shù)據(jù)能夠滿足對(duì)城市地鐵施工監(jiān)測(cè)的需求。

（2）通過統(tǒng)計(jì)分析方法對(duì)沙子口和勝利橋兩處事故點(diǎn)長(zhǎng)時(shí)序地表形變監(jiān)測(cè)數(shù)據(jù)進(jìn)行對(duì)比分析，將長(zhǎng)時(shí)序累積形變量轉(zhuǎn)換為日均沉降速率，提出地表沉降監(jiān)測(cè)早期識(shí)別指標(biāo)為單周期日均沉陷速率大于0.02 mm/d，施工預(yù)警閾值為連續(xù)3 個(gè)觀測(cè)周期（36 d）大于0.02 mm/d。

因此，基于改進(jìn)的SBAS-InSAR 獲取的高精度、長(zhǎng)時(shí)序地表形變監(jiān)測(cè)信息能夠?yàn)榈罔F施工建設(shè)中的塌陷早期識(shí)別和監(jiān)測(cè)預(yù)警提供輔助決策依據(jù)，本研究結(jié)果可為在建地鐵工程及未來規(guī)劃線路的建設(shè)提供大范圍監(jiān)測(cè)的理論技術(shù)支撐。