SBAS InSAR與GM(1,1)模型在青島地鐵13號線地表形變監測中的應用

郭在潔 陶秋香 王鳳云 韓宇

山東科技大學測繪與空間信息學院，山東青島266590

地鐵的建設和運營中，由于人為因素和自然因素的影響，地質條件發生變化，可能引發地鐵沿線地表形變和坍塌［1-2］，對沿線地面建筑物、構筑物、地下管線、各種設施、城市道路等造成不同程度的損害。因此，對地鐵沿線的地表形變進行監測極為重要［3］。

隨著空間對地觀測技術的不斷發展，合成孔徑雷達干涉測量（Interferometric Synthetic Aperture Radar，InSAR）技術逐漸成為城市地面沉降、礦區地面沉降、山體滑坡等地表形變監測的重要技術手段［4］。賈煦等［5］利用歐空局37景ENVISAT ASAR數據，基于永久散射體合成孔徑雷達干涉測量（Persistent Scatterers Interferometric Synthetic Aperture Radar，PS-InSAR）技術獲取北京地鐵15號線的地表形變情況。劉琦等［6］利用85景Sentinel-1A SAR影像數據和PS-InSAR技術獲取了佛山市地鐵沿線的地表形變情況。劉凱斯等［7］基于PS-InSAR技術和熵值法研究了北京地鐵6號線沿線的地表形變分布情況。石曉宇等［8］利用合成孔徑雷達差分干涉測量（Differential InSAR，D-InSAR）技術和灰色理論模型GM（1，1）（Grey Model，GM）對礦區沉降進行了監測與預測。劉哲［9］基于DInSAR技術和GM（1，1）模型對礦山地表形變進行監測與預測。

本文利用小基線集（Small BAseline Subset，SBAS）InSAR技術對37景Sentinel-1A SAR影像數據進行處理，獲取了青島地鐵13號線積米崖站—世紀大道站沿線的地表形變信息，并對重點沉降區域進行分析；通過均勻提取地鐵中心線上高相干點的地表形變值，計算地鐵中心線的形變梯度，進而分析地鐵沿線的穩定性；通過提取地表形變較嚴重的高相干點的時序累積形變信息，利用GM（1，1）模型對其地表形變過程進行數值摸擬和預測，研究GM（1，1）模型在地鐵地表形變監測中的應用。

1 研究區概況和數據來源

1.1 研究區概況

青島地鐵13號線坐落于青島市黃島區，于2015年開建，2018年底正式開通運營。黃島區位于北緯35°35′～36°08′，東經119°30′～120°11′，南臨黃海，北靠膠州市，西鄰諸城市、五蓮縣和日照市。年平均降水量696.9 mm，年平均氣溫12.7℃，空氣濕度大，四季分明，具有明顯的海洋氣候特點。

1.2 數據來源

哨兵1號（Sentinel-1）衛星是全球環境與安全監測（Global Monitoring for Environment and Security，GMES）中的地球觀測衛星，由兩顆衛星組成，搭載有C波段的合成孔徑雷達。首顆衛星于2014年4月發射，其重訪周期為12 d。兩顆衛星協同工作，重訪周期縮短為6 d［6］。其中，干涉寬幅模式的幅寬為250 km，地面分辨率為5 m×20 m，數據可以從歐洲空間局網（https：//scihub.copernicus.eu/dhus/）下載［7］。選取覆蓋研究區的37景C波段、VV極化的Sentine-1A降軌數據，時間跨度為2018年1月16日至2020年8月9日。此外還用到了美國國家航空航天局（National Aeronautics and Space Administration，NASA）90 m的SRTM DEM高程數據消除地形相位引起的誤差。

2 原理與方法

2.1 SBAS InSAR原理

假設某一個區域內有t0，t1，…，tN（t0，t1，…，tN為影像對應的時刻）共N+1幅單視復數（Single Look Complex，SLC）SAR影像，通過設置時間基線和空間基線閾值可以得到M幅差分干涉圖，其中M的取值范圍［10］為(N+1)/2≤M≤N( )

N+1/2。利用時刻tA和tB獲取的2景影像生成第i幅差分干涉圖（i=1，2，…，M，且時刻tA早于時刻tB），那么生成的差分干涉圖上像元(x，r)（x為方位向坐標，r為距離向坐標）的解纏后的差分干涉相位Δφi(x，r)可表示為［10］

式中：φ(tA，x，r)和φ(tB，x，r)分別為像元(x，r)在時刻tA與tB相對于參考時刻t0的雷達視線向的累積形變相位；λ為雷達的波長；d(tA，x，r)和d(tB，x，r)分別為時刻tA與tB相對于參考時刻t0的雷達視線向的累積形變量。

對于參考時刻t0，有d(t0，x，r)≡0［10］。

若 令φ=[φ(t1，x，r)，φ(t2，x，r)，…，φ(tN，x，r)]T，Δφ=[Δφ1，Δφ2，…，ΔφM]T，則主影像和輔影像的時間序列TIEi、TISi分別為

如果TIEi＞TISi，則任意的差分干涉對對應的相位值為

式中：φ(TIEi)為主影像的相位值，φ(TISi)為輔影像的相位值。

對任一像元，可以列出一個SBAS InSAR的線性形變模型［10］

式中：A是一個M×N的矩陣，矩陣行對應差分干涉圖，列對應SAR影像。

當所選數據均屬于一個基線集，即基線集數量L=1時，A為N階矩陣。若M=N，則方程只有唯一解；若M＞N，則方程的解不唯一。矩陣可以利用最小二乘法解得φ的估計值［10］，即

通常情況下，SAR影像的數據集分散在幾個子集中。當L＞1時ATA為奇異矩陣，矩陣A的秩為N-L+1，方程的解不唯一。為了得到唯一解，SBAS InSAR采用奇異值分解（Singular Value Decomposition，SVD）方法，將多個基線集聯合，求出最小范數意義上的最小二乘解，得到累積形變量［10］。

2.2 GM（1，1）模型

鄧聚龍教授在20世紀80年代提出了灰色系統理論，并用該理論解決信息不完備系統的復雜問題，根據已知數據建立模型來預測未來數據的發展趨勢［11］。GM（1，1）模型是灰色1階、1個變量的模型，具體建模過程如下。

式（11）為數列的預測公式，是對一次累加生成數列的預測值［11-12］。

利用GM（1，1）模型對SBAS InSAR監測的地表形變量進行數值摸擬和形變預測，具體步驟如下：

1）提取感興趣高相干像元（特征點）各成像時刻的累積形變量；

2）利用GM（1，1）模型對提取的特征點的前m個成像時刻的地表形變量進行數值摸擬，預測后（N+1-m）個成像時刻的形變量和后續SBAS InSAR未監測到的時刻的地表形變量；

3）對比分析SBAS InSAR地表形變監測結果和GM（1，1）模型模擬與預測的地表形變結果，對GM（1，1）模型進行精度驗證，為地鐵沿線后期沉降提供參考。

數據處理的具體流程見圖1。

圖1 數據處理流程

3 監測結果分析

3.1 地鐵沿線地表形變

青島地鐵13號線于2017年10月實現洞通，2018年2月實現軌通，2018年12月26日正式開通運營。利用圖1所示的數據處理流程和方法，對青島地鐵13號線積米崖站—世紀大道站37景Sentine-1A SAR數據進行處理，獲取2018年1月16日至2020年8月9日地鐵沿線的地表形變信息。對地鐵沿線600 m緩沖區內的地表形變進行分析得到地鐵沿線形變速率，見圖2?？芍罕O測期內地鐵沿線存在不同程度的地表形變，大部分區域的形變速率集中在-5～5 mm/a，形變比較明顯的區域為學院路站、辛屯站、兩河站、鳳凰山路站附近，最大沉降速率分別為-10、-18、-9、-10 mm/a。

從形變明顯的四個區域分別選取地表形變嚴重的點P1、P2、P3、P4（參見圖2），其地表形變過程見圖3。分別從建設期（運營前）和運營期進行分析。

圖2 地鐵沿線600 m緩沖區內的形變速率

圖3 4個特征點的時序形變結果

由圖3可知：①對于P1點，在建設期一直處于快速波動下沉狀態，在2018年月7月27日后趨于穩定，沉降值為-16 mm；地鐵運營后開始“沉降→抬升→……→沉降”的波動式緩慢沉降。②對于P2點，在建設期沉降值在2018年8月20日出現了較大的波動；在地鐵即將運營時，沉降值為-13 mm；地鐵運營后一直處于波動式下沉狀態；2019年8月27日后，沉降趨于穩定；2019年12月25日后沉降值直線式增大，其快速下沉的原因須實地勘察了解。③對于P3點，在建設期沉降值出現較大幅度的波動；地鐵即將運營時，沉降值僅為-8 mm；地鐵運營后呈現出波動式下沉，但是波動幅度并不大。④對于P4點，在建設期處于波動式下沉，但是整體波動幅度較?。坏罔F即將運營時，沉降值僅為-7 mm；運營后至2019年7月10日，該處的沉降處于直線式增大狀態；之后沉降趨于穩定。⑤至2020年8月9日，P1、P2、P3、P4點的累積沉降量分別為-25、-51、-22、-25 mm。

3.2 地鐵沿線地表穩定性

地表形變梯度描述地鐵沿線一定區間內形變速率或形變量變化快慢的過程，主要用來反映形變在空間上的均勻性［4］。不均勻地表形變所帶來的危害遠高于均勻地表形變。分析地鐵沿線的地表形變梯度對于地鐵線路的建設和后期的穩定運營具有重要意義。

以積米崖站為起點，世紀大道站為終點，分別提取地鐵中心線在2018年2月9日、2018年12月18日、2020年8月9日的累積形變量，見圖4。

圖4 地鐵中心線地表累積形變量

從圖4知：①在監測初期，地鐵沿線地表處于穩定狀態，地表形變非常小，幾乎無沉降發生。②相比監測初期，在地鐵即將開通時期（2018年12月18日），地表形變量和波動幅度都發生了較大變化，在距離起點7.0～7.8 km內，地表由抬升變為下沉，形變由2 mm減小到-18 mm。③在監測末期（2020年8月9日），地鐵沿線的地表形變波動趨勢基本和地鐵即將開通時期的趨勢一致，但地表形變量增大，在距離起點5.0～5.8 km處，沉降量由4 mm增加到14 mm；在18.5 km處，沉降量由5 mm增加到28 mm，變化幅度較大。

為了分析地鐵沿線地表形變的穩定性，沿著地鐵中心線的走向，提取高相干像元2020年8月9日的地表形變值，計算兩個高相干像元間的形變梯度。計算公式為

式中：kz，z+1為第z和第z+1像元間的地表形變梯度；xz為第z像元的地表形變值；xz+1為第z+1像元的地表形變值；s為第z和第z+1像元間的距離。

根據GB 50157—2013《地鐵設計規范》，正線的最大坡度不宜大于0.3，困難地段不高于0.35。算得的地鐵中心線地表形變梯度見圖5?？芍罔F沿線地表形變梯度整體處于穩定狀態，形變梯度大于0.3的位置主要在距離起點1.0、7.0～9.0、16.4、15.5、18.5 km處。其中，18.5 km處的地表形變梯度達到最大值1.3，在日常地鐵運營和維護中，須重點監控該路段。

圖5 地鐵中心線地表形變梯度

以地表形變梯度分析地鐵沿線地表形變的均勻性可以作為地鐵建設和運營的重要參考依據，同時須要在運營過程中加強監測確認。

3.3 GM（1，1）模型在地表形變監測中的應用

以地表形變較嚴重的P1、P2、P3、P4點為例，利用GM（1，1）模型對4個點SBAS InSAR監測的前31期、34期、37期地表形變數據進行模擬，對后6期、后3期以及未來3期的地表形變進行預測，結果見圖6—圖9。可知：①根據GM（1，1）模型對P1、P2、P3、P4點SBAS InSAR監測的地表形變數據進行模擬與預測的結果，得到P1、P2、P3、P4點在2018年1月16日至2020年8月9日地表一直在緩慢沉降，其沉降趨勢相似，都是以指數函數曲線的趨勢緩慢沉降，無法真實反映“沉降→抬升→……→沉降”的地表形變過程。②GM（1，1）模型模擬與預測的P2、P3點地表形變曲線要比P1、P4點地表形變曲線精度要高，也更符合實際情況。這是因為P2、P3點在37景SAR影像監測期間總體呈現出波動式下沉的趨勢；P1點在最后3個成像時刻，地表趨于穩定，甚至略有抬升；P4點在后15個成像時刻，地表趨于穩定，在最后4個成像時刻，地表逐漸抬升。③比較GM（1，1）模型模擬與預測的P2、P3點地表形變曲線發現，地表形變過程波動越小，形變數據的波動越小，GM（1，1）模型的模擬與預測精度就越高。圖6—圖9中第1個成像時刻和第2成像時刻的坡度較大也是由于地表形變過程的波動造成的。

圖6 GM（1，1）模型模擬與預測的P1點地表形變曲線

圖7 GM（1，1）模型模擬與預測的P2點地表形變曲線

圖8 GM（1，1）模型模擬與預測的P3點地表形變曲線

圖9 GM（1，1）模型模擬與預測的P4點地表形變曲線

可見，GM（1，1）模型更適合對地表持續沉降或持續抬升且形變過程波動較小的特征點進行模擬和預測。

4 結論

1）利用SBAS InSAR技術能夠監測和分析地鐵沿線、周邊區域及主要特征點大范圍、長時間序列的地表形變時空分布，直觀形象地展示其地表形變過程。

2）通過地鐵中心線形變梯度分析，可以找出形變梯度偏大的區域。在日常地鐵運營和維護中，須要重點監控這些路段。

3）GM（1，1）模型模擬和預測的地表形變曲線呈現出指數函數曲線的持續沉降趨勢，更適合對地表持續沉降或持續抬升且形變過程波動較小的特征點進行模擬和預測。

本文研究成果將有助于推進InSAR技術和GM（1，1）模型在城市道路和地鐵沿線地表形變監測中的應用。