地鐵施工塌陷的SBAS-InSAR長時序監測與早期識別

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

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

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

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

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

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

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

1 基本情況

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

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

2 實驗數據和方法

2.1 SAR 數據

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

圖1 時空基線圖Fig.1 Space-time baseline map

2.2 SBAS-InSAR 數據處理方法

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

圖2 SBAS-InSAR 技術流程圖Fig.2 SBAS-InSAR technical flowchart

（1）將監測區域的M+1 景SAR 影像數據集根據給定的干涉組合條件，生成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）以相干系數為基礎生成掩膜圖，如圖3（b）所示，在去地形干涉圖中選擇高相干點進行相位解纏處理，如圖3（c）所示，并對相干點進行定標處理。

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

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

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

當像素格網為1 時，

當像素格網為3 時，

當像素格網為4 時，

為了確保干涉測量網絡的一致性，在網絡感知中重新估計軌道參數。通過i反轉線性系統來估計每個SAR 采集的軌道參數。

結合重新估計的軌道參數，生成與干涉測量網絡一致的軌道校正圖，并校正每個干涉圖。

假設使用覆蓋同一研究區域的N+1 幅SAR 影像圖，則按時間順序對它們進行排序：

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

對生成的多視差分干涉圖進行濾波并解開相位，使用相干性在未包覆的差分干涉圖中選擇適當的參考點。對于任意2 個微分干涉圖（ta

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

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

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

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

變形速率的積分在干涉圖的主圖像和從圖像之間的時間間隔內。

直接在主圖像和從圖像間隔之間的每個時間間隔處進行速率積分。以矩陣形式書寫：

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

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

式（18）中：θ為入射角；Rd＝287.05 J·kg-1·K-1，Rν＝461.495 J·kg-1·K-1，分別是干燥和濕潤環境下的氣體常數；gm是z和zref之間重力加速度的加權平均值；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 結果與討論

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

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

以事故點為中心，選取方圓300 m 范圍內的時序形變點開展數理統計分析。分別在2 個事故點提取分析多個點位的時序形變信息，如圖5 所示。

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

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

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

（1）將時間信息和沉降信息一一對應：D0-S0，D1-S1，D2-S2，D3-S3，…，Dn-Sn；

（2）計算日均地表形變速率

式（19）中：Di（i=1，2，3，…，n）是日期，Si（i=1，2，3…，n）是沉降量。因為用的是前1 個周期的日期減去后1 個周期的日期，所以地面發生沉降時，沉降速率為正值。

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

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

實驗數據說明：日均沉降速率大于0.02 mm/d可作為監測預警依據，若在施工期間，連續在3 個監測周期（36 d）出現日均沉降速率大于0.02 mm/d時，應對施工單位發出塌陷事故預警，提醒施工單位此地存在塌陷災害隱患，及時采取更加精密的監測手段來監測地面沉降情況并作出事故預防措施。

4 結 論

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

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

（2）通過統計分析方法對沙子口和勝利橋兩處事故點長時序地表形變監測數據進行對比分析，將長時序累積形變量轉換為日均沉降速率，提出地表沉降監測早期識別指標為單周期日均沉陷速率大于0.02 mm/d，施工預警閾值為連續3 個觀測周期（36 d）大于0.02 mm/d。

因此，基于改進的SBAS-InSAR 獲取的高精度、長時序地表形變監測信息能夠為地鐵施工建設中的塌陷早期識別和監測預警提供輔助決策依據，本研究結果可為在建地鐵工程及未來規劃線路的建設提供大范圍監測的理論技術支撐。