延安新區地表垂直形變時空演化的經驗正交函數分析

孫宇超 魏長壽 張明剛 李志進 周立人

1 山東科技大學測繪與空間信息學院，青島市前灣港路579號，266590 2 內蒙古科技大學礦業與煤炭學院，內蒙古自治區包頭市阿爾丁大街7號，014010 3 國網山東省電力公司棗莊供電公司，山東省棗莊市黃河路999號，277800

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

1 研究區概況及數據來源

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

2 研究方法與數據處理

2.1 SBAS-InSAR技術

SBAS-InSAR技術適用于分布式目標的時序分析。相比于傳統的D-InSAR技術，SBAS-InSAR技術能夠在一定程度上克服空間失相干和時間失相干，在最小形變速率標準的基礎上對數個小基線數據集進行組合，得到所需的干涉圖。相比于PS-InSAR技術，SBAS-InSAR技術能夠同時監測出線性和非線性形變，具有更強的適用性。

2.2 經驗正交函數分析

EOF分析方法的分析對象為場的時間序列，該方法能夠在先驗條件未知的前提下獲取研究區域時空數據中的空間分布特征和時間序列。隨時間變化的變量場可以被分解為只隨時間變化的時間函數以及不隨時間變化的空間函數，可以在沒有事先人為規定的前提下，由變量場序列本身特征來確定場的結構，并且從分解出來的結構中解讀出物理意義[5-6]。

2.3 數據處理

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

圖1 時空基線

將獲得的矢量結果構建成i個觀測點、j個觀測時間的矩陣X：

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

(1)

方程求解前，需要先對矩陣內的數值進行正則化，將原始數據轉化成無量綱的形式。計算得到矩陣X的協方差矩陣A為：

(2)

式中，N為SAR影像的觀測期數。根據實對稱矩陣分解定理得：

A=VΛVT

(3)

式中，對角矩陣Λ(λ1,λ2,…,λi)為矩陣X的特征值，V為特征向量，每個非零特征根對應的特征向量值為EOF模態。時間系數T可表示為：

T=VTX

(4)

特征根越大表示模態越重要，方差貢獻率也越高。每個模態的貢獻率P可表示為：

(5)

一般情況下，通過分析前幾個占比較高的EOF模態得分和對應時間系數，就可以解釋原始數據的主要時空變化規律。

3 結果與討論

3.1 SBAS-InSAR形變結果

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

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

3.2 EOF分析地面沉降時空分布特征

對2016-07～2020-11共43景Sentinel-1A影像得到的SBAS點進行EOF分析，結果如表1所示。由表可見，第1模態的方差貢獻率高達85.24%，第2模態為12.52%，第3模態為0.71%。前2個模態的累積方差貢獻率為97.76%，超過總體的95%，因此前2個模態就可以解釋延安新區的地表形變時空演化特征。

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

3.2.1 第1模態特征

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

圖3 第1模態與年平均形變速率對比

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

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

圖5 高程與年平均形變速率對比

由圖6可見，第1模態的時間系數均為負值。2016-07～2017-01時間系數絕對值快速上升時段正是延安新區快速建設時期。2017-01～2020-11時間系數絕對值維持在0.9～1的范圍內，絕對值在2018-10達到最大，說明此時空間分布最為典型，而且這個時間點正是延安新區(北區)建設完成的節點，隨后整體的地表形變逐漸趨于穩定狀態。通過分析延安新區抬升與沉降的形成機理和發生變化的時間節點可知，挖山造地正是影響延安新區(北區)地表形變的最主要原因。

圖6 第1模態時間系數

3.2.2 第2模態特征

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

圖7 第2模態與年平均形變速率對比

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

圖8 時間系數與特征點

4 結 語

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

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