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

郭在潔 陶秋香 王鳳云 韓宇

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

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

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

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

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

1.1 研究區(qū)概況

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

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

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

2 原理與方法

2.1 SBAS InSAR原理

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

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

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

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

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

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

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

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

式中：A是一個(gè)M×N的矩陣，矩陣行對應(yīng)差分干涉圖，列對應(yīng)SAR影像。

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

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

2.2 GM（1，1）模型

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

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

利用GM（1，1）模型對SBAS InSAR監(jiān)測的地表形變量進(jìn)行數(shù)值摸擬和形變預(yù)測，具體步驟如下：

1）提取感興趣高相干像元（特征點(diǎn)）各成像時(shí)刻的累積形變量；

2）利用GM（1，1）模型對提取的特征點(diǎn)的前m個(gè)成像時(shí)刻的地表形變量進(jìn)行數(shù)值摸擬，預(yù)測后（N+1-m）個(gè)成像時(shí)刻的形變量和后續(xù)SBAS InSAR未監(jiān)測到的時(shí)刻的地表形變量；

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

數(shù)據(jù)處理的具體流程見圖1。

圖1 數(shù)據(jù)處理流程

3 監(jiān)測結(jié)果分析

3.1 地鐵沿線地表形變

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

從形變明顯的四個(gè)區(qū)域分別選取地表形變嚴(yán)重的點(diǎn)P1、P2、P3、P4（參見圖2），其地表形變過程見圖3。分別從建設(shè)期（運(yùn)營前）和運(yùn)營期進(jìn)行分析。

圖2 地鐵沿線600 m緩沖區(qū)內(nèi)的形變速率

圖3 4個(gè)特征點(diǎn)的時(shí)序形變結(jié)果

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

3.2 地鐵沿線地表穩(wěn)定性

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

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

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

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

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

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

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

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

以地表形變梯度分析地鐵沿線地表形變的均勻性可以作為地鐵建設(shè)和運(yùn)營的重要參考依據(jù)，同時(shí)須要在運(yùn)營過程中加強(qiáng)監(jiān)測確認(rèn)。

3.3 GM（1，1）模型在地表形變監(jiān)測中的應(yīng)用

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

圖6 GM（1，1）模型模擬與預(yù)測的P1點(diǎn)地表形變曲線

圖7 GM（1，1）模型模擬與預(yù)測的P2點(diǎn)地表形變曲線

圖8 GM（1，1）模型模擬與預(yù)測的P3點(diǎn)地表形變曲線

圖9 GM（1，1）模型模擬與預(yù)測的P4點(diǎn)地表形變曲線

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

4 結(jié)論

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

2）通過地鐵中心線形變梯度分析，可以找出形變梯度偏大的區(qū)域。在日常地鐵運(yùn)營和維護(hù)中，須要重點(diǎn)監(jiān)控這些路段。

3）GM（1，1）模型模擬和預(yù)測的地表形變曲線呈現(xiàn)出指數(shù)函數(shù)曲線的持續(xù)沉降趨勢，更適合對地表持續(xù)沉降或持續(xù)抬升且形變過程波動較小的特征點(diǎn)進(jìn)行模擬和預(yù)測。

本文研究成果將有助于推進(jìn)InSAR技術(shù)和GM（1，1）模型在城市道路和地鐵沿線地表形變監(jiān)測中的應(yīng)用。