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

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

1 引 言

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

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

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

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

2.1 研究區(qū)概況

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

圖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°，時(shí)間跨度為2019年1月～2021年3月，共計(jì)769天，影像具體參數(shù)如表1所示。同時(shí)采用定位精度優(yōu)于 5 m的AUX_POEORB精密軌道星歷參數(shù)用于主從影像間的空間基線估計(jì)，選取SRTM-1的 30 m分辨率的DEM用于地形相位去除與地理編碼。

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

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

3.1 PS-InSAR技術(shù)

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

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

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

圖3 PS時(shí)空基線分布

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

(1)

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

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

從圖4(a)、(b)可以看出，運(yùn)城市南部區(qū)域有明顯的沉降，其他地區(qū)形變?cè)诳臻g上分布不均勻，且發(fā)生的沉降在東西方向上呈帶狀分布，形成了永濟(jì)-絳縣、稷山-新絳縣的帶狀沉降區(qū)域。2019-1到2021-2監(jiān)測(cè)時(shí)間段內(nèi)的沉降速率在 -28.50 mm/y～20.47 mm/y之間，最大沉降速率為 28.50 mm/y。運(yùn)城市中心城區(qū)沉降明顯，出現(xiàn)了多個(gè)沉降漏斗，運(yùn)城關(guān)公機(jī)場(chǎng)周圍及空港開(kāi)發(fā)區(qū)、運(yùn)城學(xué)院等地沉降最為嚴(yán)重，沉降速率介于 -28.50 mm/y～ -7.96 mm/y之間，平均沉降速率為 23.15 mm/y。稷山縣城區(qū)地表穩(wěn)定，沒(méi)有沉降發(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技術(shù)

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

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

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

圖6 SBAS時(shí)空基線分布

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

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

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

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

4 結(jié)果對(duì)比及分析

4.1 交叉驗(yàn)證

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

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

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

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

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

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

圖11 特征點(diǎn)沉降序列

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

圖12 PS與SBAS年均沉降速率關(guān)系圖

4.2 原因分析

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

圖13 運(yùn)城市地裂縫分布

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

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

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

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

5 結(jié) 論

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