孫宇超 魏長(zhǎng)壽 張明剛 李志進(jìn) 周立人

1 山東科技大學(xué)測(cè)繪與空間信息學(xué)院，青島市前灣港路579號(hào)，266590 2 內(nèi)蒙古科技大學(xué)礦業(yè)與煤炭學(xué)院，內(nèi)蒙古自治區(qū)包頭市阿爾丁大街7號(hào)，014010 3 國(guó)網(wǎng)山東省電力公司棗莊供電公司，山東省棗莊市黃河路999號(hào)，277800

延安新區(qū)位于濕陷性黃土溝壑地帶，土地結(jié)構(gòu)疏松，極易發(fā)生阻礙工程進(jìn)展的地質(zhì)災(zāi)害[1]，因此有必要對(duì)該地區(qū)地表形變時(shí)空變化進(jìn)行分析。目前，有學(xué)者對(duì)延安新區(qū)沉降機(jī)理進(jìn)行研究[2]，但是對(duì)沉降與抬升區(qū)域的時(shí)空分析還不夠深入。本文使用小基線集合成孔徑雷達(dá)干涉測(cè)量SBAS-InSAR技術(shù)[3]，利用覆蓋延安新區(qū)2016-07～2020-11的43景Sentinel-1A升軌雷達(dá)影像數(shù)據(jù)，獲取延安新區(qū)的地表形變信息，并運(yùn)用經(jīng)驗(yàn)正交函數(shù)(empirical orthogonal function，EOF)對(duì)SBAS點(diǎn)進(jìn)行分解，挖掘延安新區(qū)地面變化的主要時(shí)空演化特征，為該地區(qū)的城市建設(shè)、災(zāi)害防控等提供監(jiān)測(cè)依據(jù)。

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

延安新區(qū)(北區(qū))施工時(shí)間為2013-04～2018-08，最大挖方厚度與最大填方厚度分別為118 m和112 m，總施工面積22.3 km2。降雨是延安新區(qū)地下水的主要補(bǔ)給方式[4]。本文使用Sentinel-1A衛(wèi)星2016-07-03～2020-11-27共43幅升軌影像，數(shù)據(jù)類型為單視復(fù)數(shù)，模式為干涉寬幅，極化方式為VV極化。利用每幅影像對(duì)應(yīng)的精密軌道數(shù)據(jù)校正衛(wèi)星軌道，使用30 m分辨率的SRTM DEM作為地形數(shù)據(jù)反演形變速率。

2 研究方法與數(shù)據(jù)處理

2.1 SBAS-InSAR技術(shù)

SBAS-InSAR技術(shù)適用于分布式目標(biāo)的時(shí)序分析。相比于傳統(tǒng)的D-InSAR技術(shù)，SBAS-InSAR技術(shù)能夠在一定程度上克服空間失相干和時(shí)間失相干，在最小形變速率標(biāo)準(zhǔn)的基礎(chǔ)上對(duì)數(shù)個(gè)小基線數(shù)據(jù)集進(jìn)行組合，得到所需的干涉圖。相比于PS-InSAR技術(shù)，SBAS-InSAR技術(shù)能夠同時(shí)監(jiān)測(cè)出線性和非線性形變，具有更強(qiáng)的適用性。

2.2 經(jīng)驗(yàn)正交函數(shù)分析

EOF分析方法的分析對(duì)象為場(chǎng)的時(shí)間序列，該方法能夠在先驗(yàn)條件未知的前提下獲取研究區(qū)域時(shí)空數(shù)據(jù)中的空間分布特征和時(shí)間序列。隨時(shí)間變化的變量場(chǎng)可以被分解為只隨時(shí)間變化的時(shí)間函數(shù)以及不隨時(shí)間變化的空間函數(shù)，可以在沒有事先人為規(guī)定的前提下，由變量場(chǎng)序列本身特征來(lái)確定場(chǎng)的結(jié)構(gòu)，并且從分解出來(lái)的結(jié)構(gòu)中解讀出物理意義[5-6]。

2.3 數(shù)據(jù)處理

首先進(jìn)行SBAS-InSAR數(shù)據(jù)處理，將2016-10-19的影像設(shè)置為超級(jí)主影像，設(shè)置最大空間基線為臨界基線的37%，最大時(shí)間基線為175 d，生成的時(shí)空基線如圖1所示。對(duì)所有干涉像對(duì)進(jìn)行干涉處理前首先需要進(jìn)行多視處理，根據(jù)以往經(jīng)驗(yàn)，設(shè)置多視視數(shù)為4∶1。然后使用Goldstein方法對(duì)干涉結(jié)果進(jìn)行濾波處理，使用Delaunay 最小費(fèi)用流法(minimum cost flow，MCF)進(jìn)行3D相位解纏。在選擇地面控制點(diǎn)時(shí)，應(yīng)該避開形變區(qū)域，使用奇異值分解法(singular value decomposition，SVD)反演獲得研究區(qū)域可靠的形變序列。最后，將SBAS-InSAR處理后得到的柵格結(jié)果轉(zhuǎn)為矢量格式，方便接下來(lái)進(jìn)行EOF分析。

圖1 時(shí)空基線

將獲得的矢量結(jié)果構(gòu)建成i個(gè)觀測(cè)點(diǎn)、j個(gè)觀測(cè)時(shí)間的矩陣X：

i=1,2,…,m,j=1,2,…,n

(1)

方程求解前，需要先對(duì)矩陣內(nèi)的數(shù)值進(jìn)行正則化，將原始數(shù)據(jù)轉(zhuǎn)化成無(wú)量綱的形式。計(jì)算得到矩陣X的協(xié)方差矩陣A為：

(2)

式中，N為SAR影像的觀測(cè)期數(shù)。根據(jù)實(shí)對(duì)稱矩陣分解定理得：

A=VΛVT

(3)

式中，對(duì)角矩陣Λ(λ1,λ2,…,λi)為矩陣X的特征值，V為特征向量，每個(gè)非零特征根對(duì)應(yīng)的特征向量值為EOF模態(tài)。時(shí)間系數(shù)T可表示為：

T=VTX

(4)

特征根越大表示模態(tài)越重要，方差貢獻(xiàn)率也越高。每個(gè)模態(tài)的貢獻(xiàn)率P可表示為：

(5)

一般情況下，通過(guò)分析前幾個(gè)占比較高的EOF模態(tài)得分和對(duì)應(yīng)時(shí)間系數(shù)，就可以解釋原始數(shù)據(jù)的主要時(shí)空變化規(guī)律。

3 結(jié)果與討論

3.1 SBAS-InSAR形變結(jié)果

從Sentinel-1A雷達(dá)數(shù)據(jù)集中可以得到延安新區(qū)視線向年平均地面形變速率。已知雷達(dá)波入射角，就可以把視線向的形變轉(zhuǎn)換為垂直向的形變(圖2)，其形變結(jié)果和空間分布與前人所得結(jié)果基本一致[7]。由圖可見，研究區(qū)域整體形變速率為-56～32 mm/a，絕大部分區(qū)域形變比較穩(wěn)定，形變量為-5～5 mm/a。西北至東南方向有一條帶狀沉降區(qū)域，東南方向沉降最嚴(yán)重，沉降速率達(dá)-56～-45 mm/a。東北部可以看出有明顯沉降，沉降速率為-45～-20 mm/a。該區(qū)域同時(shí)存在明顯的地表抬升，抬升速率為15～32 mm/a，但由于該區(qū)域失相干嚴(yán)重，SBAS點(diǎn)并沒有將其完全覆蓋。

圖2 2016～2020年延安新區(qū)(北區(qū))年平均形變速率

3.2 EOF分析地面沉降時(shí)空分布特征

對(duì)2016-07～2020-11共43景Sentinel-1A影像得到的SBAS點(diǎn)進(jìn)行EOF分析，結(jié)果如表1所示。由表可見，第1模態(tài)的方差貢獻(xiàn)率高達(dá)85.24%，第2模態(tài)為12.52%，第3模態(tài)為0.71%。前2個(gè)模態(tài)的累積方差貢獻(xiàn)率為97.76%，超過(guò)總體的95%，因此前2個(gè)模態(tài)就可以解釋延安新區(qū)的地表形變時(shí)空演化特征。

表1 延安新區(qū)(北區(qū))地面沉降EOF分解方差貢獻(xiàn)率

3.2.1 第1模態(tài)特征

第1模態(tài)的方差貢獻(xiàn)率為85.24%，特征向量值有正有負(fù)，正值和負(fù)值覆蓋區(qū)域與年平均形變速率圖大致相同。作為占比最高的一個(gè)模態(tài)，第1模態(tài)反映的是延安新區(qū)(北區(qū))地表形變的整體空間分布特征。使用反距離權(quán)重法(inverse distance weighted，IDW)對(duì)SBAS點(diǎn)和第1模態(tài)的得分進(jìn)行插值，結(jié)果如圖3所示。由圖可見，兩者的空間分布高度相似，相關(guān)性為90.43%，因此能準(zhǔn)確反映出延安新區(qū)的沉降和抬升區(qū)域。

圖3 第1模態(tài)與年平均形變速率對(duì)比

由圖3(a)可見，地表形變區(qū)域與填挖方區(qū)域有著高度一致的對(duì)應(yīng)關(guān)系。為進(jìn)一步證實(shí)這一結(jié)論，選擇覆蓋延安新區(qū)(北區(qū))的2018-08-22～2020-11-27共22景Sentinel-1A數(shù)據(jù)進(jìn)行SBAS-InSAR分析，得到的年平均形變速率如圖4所示。選擇2個(gè)橫穿填挖方區(qū)域的剖面CC′和DD′，對(duì)比剖面上各點(diǎn)的形變速率和2012年施工前的高程數(shù)據(jù)(圖5)。由圖可見，CC′剖面點(diǎn)號(hào)0～30、56～93、120～130，DD′剖面點(diǎn)號(hào)48～73、88～103為填方區(qū)域，對(duì)應(yīng)的形變速率均為負(fù)值。根據(jù)以往的研究可知，造成地表沉降的主要原因有3個(gè)：1)填方體在填充后自身發(fā)生壓縮下沉；2)原地基在填方后被施加大量載荷，破壞其內(nèi)在應(yīng)力結(jié)構(gòu)，發(fā)生壓縮下沉；3)原地基和填方體的土質(zhì)均屬于濕陷性黃土[8]，容易受到土體濕度增加的影響從而發(fā)生壓縮下沉。研究資料顯示，延安新區(qū)土體的含水率與填方深度成正比[9]，CC′與DD′剖面各填方區(qū)域中填方厚度隨高程的降低而增加，沉降速率的絕對(duì)值相比其他剖面點(diǎn)較大。CC′剖面點(diǎn)號(hào)36～50、100～114、132～142，DD′剖面點(diǎn)號(hào)25～40、114～123為挖方區(qū)域，對(duì)應(yīng)的形變速率為正值。在挖山造地的過(guò)程中，地表荷載量的減少打破原有的應(yīng)力平衡，為達(dá)到新的平衡，挖方區(qū)域的土體會(huì)出現(xiàn)一定程度的回彈[10]。由圖5可見，挖方區(qū)域中高程值越大的點(diǎn)號(hào)挖方越深，卸荷量越大，抬升速率相對(duì)較高。在挖方厚度較小的區(qū)域土體回彈量小，且后期施工建設(shè)必然會(huì)對(duì)該區(qū)域重新加載，所以沒有出現(xiàn)明顯的抬升。

圖4 2018-08～2020-11年平均形變速率

圖5 高程與年平均形變速率對(duì)比

由圖6可見，第1模態(tài)的時(shí)間系數(shù)均為負(fù)值。2016-07～2017-01時(shí)間系數(shù)絕對(duì)值快速上升時(shí)段正是延安新區(qū)快速建設(shè)時(shí)期。2017-01～2020-11時(shí)間系數(shù)絕對(duì)值維持在0.9～1的范圍內(nèi)，絕對(duì)值在2018-10達(dá)到最大，說(shuō)明此時(shí)空間分布最為典型，而且這個(gè)時(shí)間點(diǎn)正是延安新區(qū)(北區(qū))建設(shè)完成的節(jié)點(diǎn)，隨后整體的地表形變逐漸趨于穩(wěn)定狀態(tài)。通過(guò)分析延安新區(qū)抬升與沉降的形成機(jī)理和發(fā)生變化的時(shí)間節(jié)點(diǎn)可知，挖山造地正是影響延安新區(qū)(北區(qū))地表形變的最主要原因。

圖6 第1模態(tài)時(shí)間系數(shù)

3.2.2 第2模態(tài)特征

第2模態(tài)的方差貢獻(xiàn)率為12.52%，特征向量值有正有負(fù)。由圖7(b)可見，第2模態(tài)得分為正的區(qū)域分布在東北方向，得分為負(fù)的區(qū)域分布在西南方向，二者形成較為明顯的分界線。B區(qū)域正值較大，這與B區(qū)域較晚的開發(fā)建設(shè)時(shí)間有一定的關(guān)系。據(jù)資料顯示，特征點(diǎn)P1、P2、P3所在區(qū)域的工程任務(wù)完成時(shí)間為2014年初到2016年底，而特征點(diǎn)P4所在的區(qū)域直到2018-08才徹底開發(fā)完成。

圖7 第2模態(tài)與年平均形變速率對(duì)比

第2模態(tài)的時(shí)間系數(shù)和特征點(diǎn)累積形變量如圖8所示。由圖可見，第2模態(tài)的時(shí)間系數(shù)從2016-07開始一直處于增長(zhǎng)狀態(tài)，對(duì)應(yīng)選取的P1、P2、P3、P4特征點(diǎn)一直處于不斷下沉的狀態(tài)，直到2018-08后時(shí)間系數(shù)變?yōu)檎?，此時(shí)所選特征點(diǎn)的累積形變量也由正值變?yōu)樨?fù)值。隨后時(shí)間系數(shù)曲線逐漸趨于穩(wěn)定，特征點(diǎn)持續(xù)下沉的趨勢(shì)也得到緩解，說(shuō)明此時(shí)第2模態(tài)的分布特征基本成型。第2模態(tài)的分布特征與填挖方區(qū)域和時(shí)間對(duì)應(yīng)較好，延安新區(qū)(北區(qū))西南方向從2012年開始施工，2016年初整體沉降速率已經(jīng)趨于穩(wěn)定；2016-02～2018-08東北方向開始施工，包括該區(qū)域內(nèi)的挖方填方以及地表建筑建設(shè)，上述工程使得該區(qū)域在這一時(shí)段內(nèi)的形變速率加快，項(xiàng)目整體完工后，地表形變速率不斷下降直至達(dá)到穩(wěn)定狀態(tài)。

圖8 時(shí)間系數(shù)與特征點(diǎn)

4 結(jié) 語(yǔ)

1)根據(jù)獲得的地表形變信息可知，發(fā)生沉降的區(qū)域大多集中在A、B區(qū)域，最大沉降速率為-56 mm/a。發(fā)生抬升的區(qū)域覆蓋范圍較廣，最大抬升速率為32 mm/a。工程前期整體形變速率較快，工程結(jié)束時(shí)形變速率逐漸趨于穩(wěn)定。

2)通過(guò)EOF分解SBAS點(diǎn)得到空間分布特征和時(shí)間系數(shù)，第1模態(tài)的空間分布與地表形變的整體空間分布相似，相關(guān)性高達(dá)90.43%。對(duì)比施工前的高程與施工結(jié)束后的形變速率，證實(shí)工程填挖方與形變速率有關(guān)。第1模態(tài)時(shí)間系數(shù)趨于穩(wěn)定的時(shí)間點(diǎn)與工程結(jié)束的時(shí)間點(diǎn)相吻合，證明延安新區(qū)(北區(qū))的工程建設(shè)是影響其地表形變的主要原因；第2模態(tài)在空間分布上具有明顯的分界線，分界線的西南方向開發(fā)時(shí)間為2012～2016年，后期地表形變速率已經(jīng)趨于穩(wěn)定，分界線東北方向從2016年開始施工，直到2018年才建設(shè)完成，這一點(diǎn)在第2模態(tài)時(shí)間系數(shù)的變化中也有所體現(xiàn)。