運城市地表沉降的時序InSAR監(jiān)測與分析

牛寶勝，楊德宏，成飛飛，高霞霞，石建楊

(昆明理工大學國土資源工程學院，云南 昆明 650093)

1 引 言

隨著地下水開采及大型建筑物的密集建造，城市地面沉降問題日益突出，尤其是在經濟發(fā)達的城市，可能導致嚴重的財產損失和基礎下部建筑的塌陷。在我國北方多數(shù)城市和沿海許多地區(qū)，都面臨著地面沉降的問題[1，2]。地表沉降是一種由人為或自然原因引起的一種具有緩變性、不均勻及不可逆性的地質災害[3]，由地表形變引發(fā)的地質災害已經成為我國一大熱點問題。傳統(tǒng)的地面沉降監(jiān)測手段主要依賴于精密水準測量和全球定位系統(tǒng)(GPS)等技術[4]，雖然測量精度高，但是水準測量存在作業(yè)周期長、耗費大、監(jiān)測區(qū)域小、復雜地形人力無法到達等缺點；全球定位系統(tǒng)測量在高程精度方面，不能滿足大范圍形變監(jiān)測的需要[5～7]，難以滿足城市地表形變精細化防控監(jiān)測的要求。

合成孔徑雷達干涉測量技術(D-InSAR)是近二十年發(fā)展起來的極具潛力的微波遙感新技術，與常規(guī)大地測量形變監(jiān)測技術相比，InSAR具有全天時、全天候、覆蓋范圍廣、成本低、安全和觀測連續(xù)等特點，能夠獲取高分辨率、高精度的地形及地表形變信息[8，9]。為克服傳統(tǒng)D-InSAR的時空失相關、大氣延遲和軌道誤差等影響[10]，由Ferretti[11]等2001年提出的PS-InSAR技術和Berandino[12]等2002年提出的SBAS-InSAR技術等新技術應運而生，被廣泛應用到有關形變監(jiān)測的多個領域，包括由地下水開采和建筑施工等引起的城市地面沉降等方面。

地面沉降嚴重威脅著城市的發(fā)展，但目前針對運城市的地面沉降研究較少，難以對地面沉降及其時空演化特征進行詳細表征。因此，為了完善和豐富運城市地面沉降的時空變化監(jiān)測，本文從歐空局獲得2019-1～2021-3時間段內的30景Sentel-1A升軌數(shù)據(jù)，綜合利用PS-InSAR與SBAS-InSAR兩種方法對運城市的沉降進行監(jiān)測并交叉驗證，分析運城地面沉降現(xiàn)狀和空間分布特點，為城市地面沉降防控及科學治理提供重要的依據(jù)。

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

2.1 研究區(qū)概況

運城市(圖1底圖使用空間分辨率為30 m的SRTM-1DEM數(shù)據(jù))位于山西省南端，與陜西、河南兩省隔黃河而相望，北與臨汾市毗連，東與晉城接壤?？偯娣e 14 182.78 km2，介于110°15′～112°04′E，34°35′～35°49′N之間。運城全年受季風活動影響，屬暖溫帶大陸性季風氣候，全市氣溫平均14.3℃，平均年降雨量 525.2 mm，區(qū)域內地形復雜，相對高差明顯，具有平原、山地、丘陵、盆地、臺地等多種地貌類型，地處山西斷陷盆地中地裂縫分布范圍最廣、發(fā)育規(guī)模最大的區(qū)域，地面沉降的形成和發(fā)展與區(qū)內數(shù)條隱伏斷裂有關[13]，且本地水資源匱乏，年際變化大，且分布不均勻，常年農田灌溉、居民用水等一直靠抽取地下水作為供水水源，已造成地下水位大幅下降，引起地面沉降等引發(fā)了一系列地質災害問題。

圖1 研究區(qū)域位置示意圖

2.2 數(shù)據(jù)源

為了研究該區(qū)域的地表沉降，主要選取了歐空局(ESA)于2014年發(fā)射的Sentine-1A衛(wèi)星的30景干涉寬幅(interferometric wide swath，IW)升軌SAR影像，極化方式為VV，重訪周期 12 d，影像入射角為39.02°，時間跨度為2019年1月～2021年3月，共計769天，影像具體參數(shù)如表1所示。同時采用定位精度優(yōu)于 5 m的AUX_POEORB精密軌道星歷參數(shù)用于主從影像間的空間基線估計，選取SRTM-1的 30 m分辨率的DEM用于地形相位去除與地理編碼。

研究區(qū)Sentine-1A影像參數(shù) 表1

3 數(shù)據(jù)處理方法

3.1 PS-InSAR技術

PS-InSAR技術是利用多景(最少25景)同一地區(qū)的SAR影像，通過統(tǒng)計分析時間序列上幅度和相位信息的穩(wěn)定性，探測不受時間、空間基線去相關影響的穩(wěn)定點(PS)目標。將得到的PS點進行時序建模分析及線性參數(shù)求解，通過時間、空間濾波，提取PS點的形變量，進一步對各PS點在各段時間內的實際總變形量進行空間規(guī)則格網(wǎng)插值，即得到研究區(qū)內在各時間段內的形變場，具體流程如圖2所示。

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

考慮到時間基線、空間基線、多普勒中心的因素影響[14]，選取2020年3月8日成像的SAR影像為單一主影像，其余29景為輔助影像，共形成了29幅干涉像對，設置基線閾值500%，干涉對的時間、空間基線信息如圖3所示，其中黃點表示主影像位置。從圖3可以看出，空間基線主要分布在 -135 m～130 m之間，最長空間基線為 130.141 m，最短空間基線為 7.609 m。

圖3 PS時空基線分布

在PS點的選取上，采用振幅離差指數(shù)的方法對研究區(qū)內的穩(wěn)定散射目標點進行識別，相干系數(shù)閾值設置為0.75，振幅離差閾值為0.25，最后探測出PS點 647 067個。其表達式為：

(1)

將獲取的PS點利用Delaunay不規(guī)則三角網(wǎng)建立PS候選點之間的連接關系，使用帶權的最小二乘方法進行相位解纏得到每個點的線性形變速率和DEM的絕對誤差值，通過時間域和空間域濾波分離殘余相位的非線性形變和大氣相位，將線性形變和非線性形變疊加即可得到每個PS點上的形變時間序列。實驗獲得的形變結果如圖4(a)，其中不同的顏色代表不同的沉降速率，負值表示沉降，正值表示抬升。在ArcGIS軟件上采用克里金(Kriging)插值獲取運城市區(qū)的地表沉降分布，為了保證插值的準確性，選取了干涉點較多、沉降速率較快的城區(qū)進行插值[15]，插值結果如圖4(b)，不同的顏色代表不同的沉降區(qū)域，顏色越深，沉降速率越大。

圖4 PS監(jiān)測地面沉降結果

從圖4(a)、(b)可以看出，運城市南部區(qū)域有明顯的沉降，其他地區(qū)形變在空間上分布不均勻，且發(fā)生的沉降在東西方向上呈帶狀分布，形成了永濟-絳縣、稷山-新絳縣的帶狀沉降區(qū)域。2019-1到2021-2監(jiān)測時間段內的沉降速率在 -28.50 mm/y～20.47 mm/y之間，最大沉降速率為 28.50 mm/y。運城市中心城區(qū)沉降明顯，出現(xiàn)了多個沉降漏斗，運城關公機場周圍及空港開發(fā)區(qū)、運城學院等地沉降最為嚴重，沉降速率介于 -28.50 mm/y～ -7.96 mm/y之間，平均沉降速率為 23.15 mm/y。稷山縣城區(qū)地表穩(wěn)定，沒有沉降發(fā)生，但其西社鎮(zhèn)北部出現(xiàn)了以薛家莊村、麻古垛村、劉家莊村為中心的沉降漏斗，最大沉降速率為 27.50 mm/y；新絳縣城區(qū)域整體穩(wěn)定，轄區(qū)的北張鎮(zhèn)出現(xiàn)大幅度沉降，平均沉降速率為 19.10 mm/y；絳縣城區(qū)出現(xiàn)小幅沉降，轄區(qū)柳泉村-橫水鎮(zhèn)區(qū)域呈帶狀沉降，平均沉降速率為 6.44 mm/y；夏縣整體呈現(xiàn)不均勻沉降，沉降較快的區(qū)域包括尉郭鄉(xiāng)、禹王鄉(xiāng)、牛家凹村、井溝村、西陰村等。

3.2 SBAS-InSAR技術

SBAS技術是對相干目標(coherent target，CT)進行相位分析來獲取時序形變，通過選擇合適的空間基線和時間基線閾值組成差分干涉對，利用具有較短時-空基線的影像對產生干涉圖提高相干性，選取相干目標點利用線性相位變化模型進行建模和解算，并通過時空濾波去除大氣延遲，在減少DInSAR處理中的失相關影響及高程、大氣誤差的同時獲取地表的形變時間序列[12]，主要技術流程圖如圖5所示。

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

同樣以20200308影像為主影像，為了保證有足夠的干涉影像對，且限制時間和空間失相干的影響，選擇臨界基線閾值為45%和時間臨界基線閾值為365天進行短基線集干涉影像對組合。干涉處理時在方位向和距離向對影像進行了 4∶1的多視處理以削弱相位噪聲的影響。相干系數(shù)的閾值設為0.2，采用Goldstein方法進行濾波，最終生成的干涉基線分布如圖6所示，共生成284個干涉像對，最大臨界基線 257.362 m，最小臨界基線 0.030 m。

圖6 SBAS時空基線分布

相位解纏過程中，由于用于探測的緩慢失相關濾波相位像素點(slowly-decorrelating filtered phase，SDFP)點[16]是離散的，在相位解纏前，將滿足相干系數(shù)閾值條件的像元作為節(jié)點生成不規(guī)則三角網(wǎng)(Delaunay)，鑒于研究區(qū)域存在植被和水體，基于這些離散點的格網(wǎng)利用最小費用網(wǎng)絡流法(MCF)進行相位解纏。相位解纏完成后，將PS-InSAR獲取的GCP文件轉換到SAR坐標下，用于精確估算軌道參數(shù)，去除殘余相位和軌道誤差，即軌道精煉和重去平。將前面得到的相干點建立線性模型，采用奇異值分解法(SVD)估算形變速率和高程系數(shù)。

圖7 從殘余相位中分出大氣相位

圖8 剩余相位分離出非線性相位

圖9 基于SBAS-InSAR方法地面沉降監(jiān)測結果

4 結果對比及分析

4.1 交叉驗證

由于缺少實測同期水準監(jiān)測數(shù)據(jù)，本文通過比較PS(圖7)和SBAS(圖10)得到形變速率結果進行精度驗證。不難看出，兩種結果所呈現(xiàn)的沉降范圍及量級具有較好的一致性，主要沉降區(qū)分布在運城市區(qū)、稷山縣、夏縣、絳縣區(qū)域。區(qū)別在于SBAS獲得點的密度大于PS的結果，個別區(qū)域的沉降速率存在較大偏差。

兩種方法反演得到的地面沉降速率統(tǒng)計指標如表2所示，其中，由PS與SBAS兩種方法獲得的地面沉降速率平均值分別為 -3.02 mm每年(mm/y)和 -0.54 mm/y，相差 2.08 mm，兩種方法的反演結果很接近，且SBAS方法得到的相干點數(shù)是PS的15倍，結合圖7與圖9，可以看出SBAS一定程度上提高了相干點在非城區(qū)的密度。PS與SBAS方法反演地面沉降的速率分別是 -28.5 mm/y、-88.99 mm/y，分別位于稷山縣的劉家莊村和馬首官莊村。SBAS方法的標準差為5.93，稍小于PS方法的標準差6.70，表明SBAS方法反演的沉降速率分布更集中。

兩種方法沉降速率監(jiān)測結果對比 表2

上面的分析可以看出，兩種方法得到的部分點的沉降速率有較大偏差，可能受個別因素的影響(如房屋倒塌等)，因此對其沉降速率分布進行統(tǒng)計，如圖10所示，PS和SBAS方法得到的相干點速率在 -10 mm/y～0 mm/y之間的點各占相干點總數(shù)的43.12%、40.18%，說明運城市地表沉降相對穩(wěn)定，這與運城市對地下水的壓采有關，表明治理成果顯著。沉降速率大于 10 mm/y分別占16.49%、5.10%，表明小部分區(qū)域發(fā)生較大沉降，且已經形成沉降漏斗。整體來看，出現(xiàn)沉降的區(qū)域分別占整個運城市的59.61%、49.09%，結合圖7與圖10，沉降區(qū)域主要分布在運城市區(qū)、夏縣、稷山縣等地。

圖10 地面沉降速率數(shù)值分布

在沉降形成的漏斗區(qū)域，本文選取6對典型點進行分析，如圖11所示，在運城市區(qū)選取了運城學院、軍屯村、圣惠北路附近三個點，在絳縣選取了橫水鎮(zhèn)的典型點，在夏縣選取了馬村的典型點，在新絳縣選取了西行莊村的典型點，繪制了在2019年～2021年期間的時序沉降折線圖。從圖中可以看出，在研究期間，兩種數(shù)據(jù)處理方法得到的6個點的沉降趨勢具有較高的一致性，圖11(a)、(e)中PS和SDFP(SBAS)點的沉降量基本保持一致，運城學院在2020年6月18號以后SDFP點沉降加劇，在2021年后PS點的沉降超過SDFP點。圖11(b)、(c)、(d)、(f)中PS和SDFP點的沉降量保持一致，受到了共模系統(tǒng)誤差的影響，SDFP點沉降量整個大于PS點?？傮w來看，它們的沉降趨勢一直呈持續(xù)加劇趨勢。

圖11 特征點沉降序列

為進一步確定監(jiān)測結果的精度，對PS和SBAS進行內符合精度驗證，在運城市城區(qū)選取兩種方法中經緯度相同的點735個同名點，比較兩種方法的年均沉降速率。以PS點沉降速率為橫軸，SBAS點沉降速率為縱軸，得到的散點圖如圖12所示，由圖可知，兩者的相關系數(shù)R2達到了 0.887 2，進一步驗證了PS和SBAS的沉降速率具有較高的一致性。

圖12 PS與SBAS年均沉降速率關系圖

4.2 原因分析

(1)構造活動：研究區(qū)域地處山西地塹系西南端的一個半封閉斷陷盆地(運城盆地)，盆地南、東、北三面環(huán)山，是一個強烈的沉降盆地，內部形成多個斷裂帶(如圖13所示)，結合PS與SBAS的沉降速率圖，可以發(fā)現(xiàn)稷山縣與新絳縣的交界帶附、夏縣的尉郭鄉(xiāng)至鹽湖區(qū)陶村鎮(zhèn)等地的沉降是由斷裂活動引起的，說明斷裂帶的出現(xiàn)往往伴隨著地面沉降，兩者相輔相成。

圖13 運城市地裂縫分布

(2)地下水開采：在運城盆地內，多數(shù)地面沉降與地下水的超采有關[17]，地下水的超采導致地下水水位下降，從而引起可壓縮土層的壓縮，進而引起地面沉降。運城屬于嚴重缺水地區(qū)，中深層地下水是運城市工農業(yè)生產和城鎮(zhèn)生活用水的主要水源，開采量占總用水量的75%～80%。從圖14可以看出，自20世紀80年代以來，運城盆地地下水開始大量開采，持續(xù)的開采地下水致使水位不斷下降[18]，產生了沉降漏斗，2002年漏斗中心水位降到了最低點，2007年完成關井壓采121眼后，2008年水位有所回升。截至2015年底，地下水超采仍占全市總面積的28%[20]。受區(qū)域超采的影響，漏斗區(qū)域水位仍呈下降趨勢，史雙雙[19]監(jiān)測2012年運城市最大沉降速率為 50 mm/y～60 mm/y，位于夏縣禹王鄉(xiāng)、尉郭鄉(xiāng)及裴介鎮(zhèn)中間地帶，從2019年～2021年的監(jiān)測結果來看，禹王鄉(xiāng)、尉郭鄉(xiāng)及裴介鎮(zhèn)中間地帶的最大沉降速率為 42.25 mm/y，相較2012年速率減緩，這與運城市實施地下水壓采措施有關。但農業(yè)灌溉用水仍以大水漫灌為主，根據(jù)2020年的《運城市統(tǒng)計年鑒》可知，較上一年用水量增加7.3%，其中2019年底農業(yè)用水13億m3，部分農村地區(qū)仍采用地下水澆灌農田，對地面沉降起到了加速作用。可見超量開采是運城城區(qū)地面沉降形成的最主要原因。

圖14 地下水開采量與降落漏斗中心水位變化曲線

(3)地震與密集高程建筑的影響：運城地震活動比較頻繁，自公元1485年到現(xiàn)在，有記錄的Ms5級以上的破壞地震中，5.0級～5.9級10次、6.0級～6.9級2次，對運城地面沉降造成一定程度的影響。同時伴隨著人口的快速增長和經濟的發(fā)展，大規(guī)模城市工程建設的增多，也對地面沉降產生一定促進作用，截至2019年底，運城市的房屋建筑面積為 1 980.8萬m2，比上年增長18.1%，密集的高層建筑在大規(guī)模城市改造建設中地面沉降效應明顯，在研究區(qū)域中，圖15中姚孟街道周邊的沉降等就是新建樓房引起的，但其沉降量達到一定程度就會趨于穩(wěn)定。

圖15 姚孟街道周邊(光學影像)

5 結 論

本文利用PS-InSAR和SBAS-InSAR技術處理了運城市30景哨兵數(shù)據(jù)，獲取了運城市2019-1～2021-2的地面形變信息，由于沒有獲取到監(jiān)測時間段內的水準數(shù)據(jù)，通過對比分析與已有資料來驗證成果的可靠性。結果表明，運城市城區(qū)、夏縣、稷山、絳縣、新絳縣的部分地區(qū)沉降較為嚴重，PS的最大沉降速率達到了 28.50 mm/y，位于運城市市區(qū)，SBAS的最大沉降速率達到了 84.44 mm/y，其位于稷山縣。結合已有資料，地面沉降的主要原因在于地下水的超采，部分是城市工程建設增加等因素的影響，對運城市沉降防治提供了一定的參考。試驗區(qū)周圍有大量的農田，除冬季外，可能會產生嚴重的失相干，需要后續(xù)實驗進一步討論與分析。