基于PS-InSAR技術(shù)的太原地面沉降監(jiān)測(cè)研究

李 路,洪友堂

(中國(guó)地質(zhì)大學(xué)(北京)，土地科學(xué)技術(shù)學(xué)院，北京 100083)

地表形變是一種由人為或自然原因造成地表變形現(xiàn)象，世界上各個(gè)國(guó)家都受到不同程度的地表形變影響。據(jù)統(tǒng)計(jì)，自1949～2010年，我國(guó)由于地面沉降造成的經(jīng)濟(jì)損失累計(jì)高達(dá)5 000億元，同時(shí)也造成了重大人員傷亡事件發(fā)生。目前，傳統(tǒng)的地表形變測(cè)量?jī)x器主要有水準(zhǔn)儀、經(jīng)緯儀、全站儀、GPS等，它們雖然精度高但需要人員實(shí)地監(jiān)測(cè)，且耗費(fèi)大量人力物力。隨著科技發(fā)展，合成孔徑雷達(dá)差分干涉測(cè)量(DInSAR)技術(shù)憑借全天時(shí)不接觸、不易受環(huán)境影響、監(jiān)測(cè)效率高等優(yōu)勢(shì)，極其適用于對(duì)大面積地面進(jìn)行微小形變研究，精度可達(dá)毫米級(jí)，相比于傳統(tǒng)測(cè)量技術(shù)，具有極大優(yōu)勢(shì)。但是由于地球大氣層中的大氣噪聲、形變前后地面物體發(fā)生位置變化導(dǎo)致時(shí)空失相干等眾多因素存在，限制了該技術(shù)的發(fā)展。2002年，F(xiàn)erretti等人首次提出PS-InSAR技術(shù)并對(duì)研究區(qū)進(jìn)行滑坡監(jiān)測(cè)，處理結(jié)果與實(shí)測(cè)數(shù)據(jù)基本吻合。隨后，Berardino等人提出小基線集差分干涉測(cè)量(SBAS-InSAR)技術(shù)，并證明有效性。經(jīng)過(guò)十余年的發(fā)展，PS-InSAR技術(shù)與SBAS-InSAR技術(shù)得到了顯著發(fā)展與應(yīng)用。而PS-DInSAR技術(shù)主要依靠散射特性穩(wěn)定點(diǎn)目標(biāo)，而城鎮(zhèn)區(qū)域中的房角、電桿等穩(wěn)定地物較多，因此該技術(shù)極其適用于對(duì)城鎮(zhèn)區(qū)域進(jìn)行地表沉降監(jiān)測(cè)研究。

2015年，楊海瑞[1]等人研究了太原盆地地下水流和土體變形耦合的機(jī)理,并對(duì)未來(lái)20年太原盆地地下水水位及地面沉降情況進(jìn)行預(yù)測(cè)。2016年孫曉涵[2]等人全面分析了太原盆地地裂縫和地下水開(kāi)采、地面沉降在時(shí)間、空間、活動(dòng)上的相互關(guān)系。2017年趙強(qiáng)[3]等人結(jié)合太原盆地地面沉降的發(fā)育情況,從地下水開(kāi)采、粘性土分布不均、斷層構(gòu)造、地面載荷等方面,分析了影響地面沉降的因素。2018年周艷萍[4]等人根據(jù)五個(gè)沉降中心中30個(gè)典型的水準(zhǔn)觀測(cè)點(diǎn)的累積沉降量建立了灰色Verhulst預(yù)測(cè)模型，最后預(yù)測(cè)了2010年與2015年的地面沉降發(fā)展趨勢(shì)。本文利用分辨率為3 m的COSMO-SkyMed數(shù)據(jù)，提出基于PS-InSAR技術(shù)對(duì)2013年6月到2016年6月期間太原區(qū)域的地面沉降進(jìn)行分析，結(jié)果顯示太原地區(qū)共形成以萬(wàn)柏林、南堰、小店地區(qū)為中心的3個(gè)沉降漏斗，其中沉降最嚴(yán)重的是小店地區(qū)，最大累計(jì)沉降量達(dá)-152.7 mm。

1 PS-DInSAR原理及關(guān)鍵步驟

1.1 PS-DInSAR方法

該方法是意大利學(xué)者Ferretti等人于2002年提出，并對(duì)研究區(qū)進(jìn)行滑坡監(jiān)測(cè)，處理結(jié)果與實(shí)測(cè)數(shù)據(jù)基本吻合，證明了方法的有效性。該方法主要以獲得的在大時(shí)空基線情況下仍保持較高相干性(如房角、電桿等)的永久散射體(PS點(diǎn))為基礎(chǔ)，這些PS點(diǎn)受時(shí)空影響較小，更易于分離大氣、地形與軌道誤差相位，獲得精確的地表形變相位。該方法可以有效地消除或削弱大氣效應(yīng)、時(shí)空失相關(guān)的影響，超越傳統(tǒng)時(shí)間、空間基線距上的臨界基線條件限制，提高了數(shù)據(jù)利用率。經(jīng)過(guò)10多年的研究，利用該技術(shù)已經(jīng)成功在昆明、北京等地進(jìn)行地面沉降監(jiān)測(cè)研究，此外一些學(xué)者在主影像選擇、PS點(diǎn)選取、相位解纏等方面進(jìn)行了改進(jìn)，獲得了更為精確的地表結(jié)果信息。

PS-DInSAR方法的核心思想為：從一系列N景SAR影像中選取其中一景作為公共主影像，其他影像作為輔影像，將所有輔影像與主影像進(jìn)行配準(zhǔn)、重采樣、干涉處理與去除地形相位，獲得N-1景差分干涉圖，選擇PS點(diǎn)，建立函數(shù)模型，并利用回歸分析估計(jì)線性形變速率和高程殘差，利用時(shí)空濾波分離大氣相位成分、非線性形變成分與噪聲相位成分，最終獲得精確的地表形變相位信息[5]。具體流程圖如圖1所示。

1.2 方法處理關(guān)鍵步驟

1.2.1 PS點(diǎn)組合選取

在PS-DInSAR方法處理流程中，PS點(diǎn)提取是處理的第一步，PS點(diǎn)選擇的質(zhì)量和數(shù)量決定數(shù)據(jù)處理結(jié)果的優(yōu)劣性，這些點(diǎn)分布于裸露巖石、高大建筑物、大金屬支架等地物上，具有在大時(shí)空基線情況下仍保持較高相干性的特點(diǎn)。常用的點(diǎn)提取方法有相干系數(shù)法和振幅離差法。其中相干系數(shù)法為主要方法，依據(jù)公式(1)，計(jì)算目標(biāo)像元(a，b)在整個(gè)時(shí)序的相干系數(shù)，然后與給定閾值做比較，以確定其是否可選為PS點(diǎn)。通過(guò)幅度離差法可以濾除海面上的點(diǎn)。

(1)

式中，M、N為窗口的高度和寬度，φ(a,b)為地形相位項(xiàng)。

1.2.2 相位解纏優(yōu)化

相位解纏是 InSAR 數(shù)據(jù)處理流程中重要步驟之一，相位解纏的準(zhǔn)確度在很大程度上決定了獲取地表形變信息的準(zhǔn)確性，算法有路徑跟蹤算法、最小二乘算法、最小費(fèi)用流算法。本文為了更準(zhǔn)確計(jì)算殘余相位，首先通過(guò)積分求得平均DEM殘余和平均線性形變相位，再?gòu)牟罘指缮嫦辔恢袑⑵淙コ玫綒堄嘞辔?。為了提高殘余相位解纏的穩(wěn)健性，對(duì)其進(jìn)行空間域?yàn)V波處理后再解纏?？紤]到相位解纏中實(shí)際情況以提高解纏能力，引入權(quán)值思想，求得最小最優(yōu)解纏模糊數(shù)，見(jiàn)公式(2)。式中，K(a,b)為影像中(a,b)處的相位值，wm為權(quán)值，本文將像元振幅值作為權(quán)值，在M個(gè)時(shí)域內(nèi)進(jìn)行相位解纏。

(2)

2 研究區(qū)概況與數(shù)據(jù)處理

2.1 研究區(qū)概況

研究區(qū)域?yàn)樘谐菂^(qū)及其部分周邊地區(qū)，主要位于山西省中北部，太原盆地北部。地理坐標(biāo)范圍為東經(jīng)112°19′～112°48′，北緯37°40′～38°1′，覆蓋面積為36×32 km2。以研究區(qū)地形圖為底圖，疊加影像的強(qiáng)度均值圖，如圖2所示。據(jù)太原市114個(gè)二等水準(zhǔn)測(cè)量點(diǎn)實(shí)測(cè)資料，1956～2003年地面沉降范圍南北長(zhǎng)約39 km，東西寬約15 km，至2003年已形成西張、萬(wàn)柏林、下元和吳家堡4個(gè)沉降中心[6-7]。最大沉降中心為吳家堡地區(qū)， 年均沉降速率為63.0 mm/a。近年來(lái)，由于超采地下水引發(fā)的地面沉降也愈演愈烈，沉降范圍向盆地邊緣擴(kuò)展，沉降漏斗面積擴(kuò)大，對(duì)太原地區(qū)進(jìn)行地表變形監(jiān)測(cè)刻不容緩[8]。

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

圖2 研究區(qū)域

本文使用瑞士的GAMMA雷達(dá)數(shù)據(jù)處理軟件對(duì)37景高分辨率COSMO-SkyMed SAR數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，時(shí)間跨度從2013年6月～2016年6月，X波段(3.1 cm)，地面分辨率為3 m×3 m，極化方式為HH，入射角為24.94°～28.36°，右視，升軌,幅寬為40 km×47 km。采用美國(guó)NASA提供的分辨率為30 m的SRTM1 DEM數(shù)據(jù)。

表1 研究區(qū)域數(shù)據(jù)詳細(xì)情況

選擇時(shí)空基線居中的2014年8月9日影像作為單一主影像，對(duì)配準(zhǔn)好的輔影像進(jìn)行干涉處理，形成36對(duì)干涉對(duì)，如圖3所示。平均垂直基線385 d，平均時(shí)間基線179 d，可以看出基線分布較穩(wěn)定，相干性較高，DEM誤差對(duì)形變相位影響較小，結(jié)果更精確。

本文采用相干系數(shù)法和振幅離差法組合選點(diǎn)方法，相干系數(shù)閾值設(shè)置為1.5，PS點(diǎn)數(shù)為33 735個(gè)，振幅離差閾值設(shè)置為0.3，PS點(diǎn)數(shù)為3 772 458個(gè)，共選擇初始PS點(diǎn)3 787 263個(gè)。本次研究區(qū)域?yàn)樘菂^(qū)，因此PS點(diǎn)密度較高，平均每平方公里PS點(diǎn)達(dá)兩千余個(gè)，且廣泛分布于研究區(qū)域中的房頂、橋梁、火車站、機(jī)場(chǎng)等地，更有利于進(jìn)行區(qū)域地表變形分析?；谶x取的初始候選點(diǎn)，提取其高程、基線、相位等信息，并利用DEM高程數(shù)據(jù)去除地形相位，形成基于干涉點(diǎn)的差分干涉相位。對(duì)各點(diǎn)的差分干涉相位先進(jìn)行空間濾波降噪處理，采用的是50×50的均值濾波窗口，之后，對(duì)獲得的穩(wěn)健的點(diǎn)差分干涉相位進(jìn)行多路徑回歸分析處理，來(lái)估計(jì)高程改正、線性形變、殘余相位等。對(duì)于殘余相位再次進(jìn)行纏繞、濾波、解纏處理，以精確估計(jì)大氣相位。將獲得的高程改正、大氣相位代入原式進(jìn)行修正，以進(jìn)行迭代處理。通過(guò)迭代，PS點(diǎn)數(shù)量逐步降低，殘差標(biāo)準(zhǔn)差降低，點(diǎn)質(zhì)量得到提升，模型也越來(lái)越穩(wěn)健，以達(dá)到優(yōu)化目的。采用基于加權(quán)最小費(fèi)用流相位解纏算法，并使用最小二乘算法來(lái)校正基線信息，以獲得可靠的結(jié)果。另外，利用SVD算法對(duì)每個(gè)干涉點(diǎn)目標(biāo)的時(shí)間序列相位進(jìn)行處理，求算2013年6月～2016年6月期間的年形變速率與累計(jì)形變量。

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

3 沉降結(jié)果分析與驗(yàn)證

3.1 沉降結(jié)果

圖4 研究區(qū)域年沉降速率圖

以光學(xué)遙感影像作為底圖，數(shù)據(jù)處理獲得太原研究區(qū)域的年沉降速率圖，如圖4所示，將區(qū)域A、B的年沉降速率情況放大至右側(cè)。A區(qū)域?yàn)槿f(wàn)柏林地區(qū)；B地區(qū)為南堰地區(qū)；C地區(qū)為小店地區(qū)；D地區(qū)為吳家堡地區(qū)；E地區(qū)為西張地區(qū)。可以很明顯得看出，萬(wàn)柏林、南堰、小店地區(qū)有明顯沉降出現(xiàn)，而吳家堡、西張地區(qū)沉降平穩(wěn)。在太原南部的小店地區(qū)沉降現(xiàn)象最為嚴(yán)重，最大值累計(jì)沉降量為-152.7 mm，平均年沉降速率達(dá)50.901 mm/a。

3.2 地表變形分析

查閱太原的沉降歷史資料發(fā)現(xiàn)，在2000年前市區(qū)主要有四個(gè)沉降中心，分別為西張、萬(wàn)柏林、下元和吳家堡，而小店地區(qū)沉降不明顯，在2003年后，這四個(gè)地區(qū)的沉降現(xiàn)象趨于緩和，局部區(qū)域出現(xiàn)反彈，南部區(qū)域沉降現(xiàn)象加劇。因此本文主要基于此四個(gè)沉降中心與小店地區(qū)進(jìn)行地表變形分析。

A區(qū)域-萬(wàn)柏林地區(qū)，地質(zhì)構(gòu)造為太原西山向斜煤盆地，汾河流域陷落區(qū)，屬于舊有沉降區(qū)，1989年至2000年該地區(qū)年沉降速率為46.73 mm/a，隨后數(shù)年緩和很多，在2010～2013年僅為-14.6 mm/a。根據(jù)《太原市萬(wàn)柏林和平老工業(yè)區(qū)搬遷改造實(shí)施方案 (2013年～2020年)》，重點(diǎn)推進(jìn)下元、南寒兩個(gè)村的改造，從數(shù)據(jù)處理結(jié)果看，在2013年6月至2016年6月期間，該地區(qū)最大累計(jì)沉降量達(dá)-93.61 mm，年平均沉降速率達(dá)-31.20 mm/a，符合區(qū)域?qū)嶋H情況。該地區(qū)在2016年初有回升趨勢(shì)，可以初步估計(jì)該區(qū)沉降速率有所緩和，查閱資料發(fā)現(xiàn)該地區(qū)基層社區(qū)較多，共有7個(gè)社區(qū)居委會(huì)和2個(gè)村委會(huì)，該時(shí)間段內(nèi)新建高層建筑諸多，正常沉降，暫時(shí)沒(méi)有安全隱患。萬(wàn)柏林地區(qū)2013年6月～2016年6月的時(shí)間序列沉降量圖(如圖5所示)。對(duì)萬(wàn)柏林地區(qū)的37°52′5″緯度線上、東西向約2 km，做時(shí)序分析，沉降漏斗剖面圖表現(xiàn)年沉降速率變化情況，如圖6所示，從圖中可以看出該地區(qū)很明顯為沉降漏斗，最低點(diǎn)的年沉降速率最大可達(dá)31 mm/a。

圖5 A區(qū)域地表監(jiān)測(cè)累計(jì)沉降量

圖6 A區(qū)域沉降漏斗剖面

B地區(qū)-南堰地區(qū)，與吳家堡地區(qū)相近，同屬舊有沉降區(qū)，在2000年以前為太原市區(qū)最大的沉降區(qū)，本文數(shù)據(jù)處理結(jié)果發(fā)現(xiàn)研究期間該地區(qū)累計(jì)沉降量達(dá)-105.78 mm，平均沉降速率為35.26 mm/a。從遙感影像上發(fā)現(xiàn)該地區(qū)存在大面積的建筑用地，說(shuō)明正在進(jìn)行施工建設(shè)或拆遷工程，從提取的PS點(diǎn)數(shù)量較少可以說(shuō)明此地時(shí)空失相干較為嚴(yán)重，也可以證明這點(diǎn)。

C地區(qū)-小店地區(qū)，位于太原市區(qū)的東南部，晉中盆地的北端，有明顯沉降漏斗出現(xiàn)。在2000年前該地區(qū)為非沉降中心，在2003～2010年的年沉降速率為27.9 mm/a，在2010～2013年的年沉降速率為36.6 mm/a，隨著太原地區(qū)經(jīng)濟(jì)南移，使得小店地區(qū)經(jīng)濟(jì)發(fā)展迅猛，小店地區(qū)地面沉降現(xiàn)象日趨嚴(yán)重、沉降速率日趨增加、沉降中心逐漸形成。從結(jié)果看，在2013年6月～2016年6月期間，小店高新區(qū)的沉降中心范圍逐漸擴(kuò)大，沉降量逐漸變大，最大沉降量達(dá)152.7 mm，年沉降速率達(dá)50.901 mm/a。整個(gè)區(qū)域沉降趨勢(shì)較大且有加劇跡象，整個(gè)沉降漏斗面積達(dá)40 km2。在小店區(qū)域的中部、東部部分地區(qū)有較嚴(yán)重沉降現(xiàn)象，而東南部較穩(wěn)定。對(duì)緯度為37°44′40″的東西向5 km區(qū)域做時(shí)序分析，形成沉降漏斗剖面圖，如圖8所示。從圖中可以發(fā)現(xiàn)小店地區(qū)主要有2個(gè)沉降漏斗出現(xiàn)：在a處沉降漏斗區(qū)域，東西向約500 m，面積較小，然而僅僅200 m差異卻有較大沉降起伏，西向高點(diǎn)年沉降速率為32.825 mm/a，東向高點(diǎn)年沉降速率為32.034 mm/a，而兩者間的沉降最低點(diǎn)為47.71 mm/a，有每年15 mm的沉降差異?？疾彀l(fā)現(xiàn)該處有數(shù)棟高聳建筑，且地處在b處沉降漏斗區(qū)域，呈現(xiàn)自距起點(diǎn)2.7 km處(也為a沉降漏斗東向沉降最高點(diǎn))起，向東側(cè)沉降愈發(fā)嚴(yán)重，在4.13 km處達(dá)到了小店地區(qū)沉降漏斗最大年沉降速率點(diǎn)，年沉降速率達(dá)50.901 mm/a。在4.7 km處有明顯沉降浮動(dòng)出現(xiàn)，根據(jù)距離、坐標(biāo)、位置判斷為田莊斷層，與2014年劉媛媛等人進(jìn)行太原市地面沉降監(jiān)測(cè)后得出的田莊斷層結(jié)論一致，之后剖面向東進(jìn)入農(nóng)用地范圍，沉降穩(wěn)定。

D地區(qū)-吳家堡地區(qū)，自1956年沉降觀測(cè)以來(lái)，吳家堡地區(qū)一直為沉降中心，在1981～1989年年沉降速率達(dá)到了114 mm/a，隨后逐步緩和，根據(jù)最新的2010～2013年的沉降觀測(cè)數(shù)據(jù)顯示為6.6 mm/a，幾近穩(wěn)定。在本次監(jiān)測(cè)時(shí)間段內(nèi)，吳家堡地區(qū)最大沉降速率為7.929 mm/a，最大沉降量為-23.787 mm。其余多數(shù)點(diǎn)年沉降速率在2～6 mm/a之間，整個(gè)地區(qū)穩(wěn)定。

圖7 C區(qū)域地表監(jiān)測(cè)累計(jì)沉降量

圖8 C區(qū)域沉降漏斗剖面

E地區(qū)-西張地區(qū)，資料查詢發(fā)現(xiàn)該地區(qū)沉降不明顯，且在2003～2010年的年沉降速率為2.3 mm/a，而本文結(jié)果顯示該地區(qū)無(wú)沉降。

總之，根據(jù)太原市區(qū)年平均沉降速率圖與分析結(jié)果顯示，在2013～2016年期間，太原地表總計(jì)形成萬(wàn)柏林、南堰、小店三個(gè)沉降中心，吳家堡、西張地區(qū)地表沉降現(xiàn)象較為穩(wěn)定，而小店地區(qū)沉降現(xiàn)象最為嚴(yán)重，應(yīng)引起注意。

3.3 可靠性分析

2014年，劉媛媛[7]等人使用20景ENVISAT ASAR數(shù)據(jù)采用PS-InSAR技術(shù)對(duì)太原市區(qū)進(jìn)行地面沉降監(jiān)測(cè)研究，結(jié)果顯示2010～2013年太原形成吳家堡、小店兩個(gè)沉降中心，并且在小店地區(qū)的剖線上發(fā)現(xiàn)田莊斷層，本文也在小店地區(qū)緯度為37°44′40″東西向的5 km剖線上發(fā)現(xiàn)田莊斷層，且剖面趨勢(shì)與劉媛媛等人研究結(jié)果基本一致。2015年，劉瑾[9]等人對(duì)太原地表各階段的各沉降漏斗中心的地下水水位變化與地面沉降速率進(jìn)行相關(guān)聯(lián)性分析，北部區(qū)域基本停止沉降， 中部區(qū)域趨緩，南部區(qū)域呈惡化趨勢(shì)。2018年，周艷萍[4]基于灰色Verhulst模型對(duì)山西太原地面沉降進(jìn)行趨勢(shì)分析，得出西張穩(wěn)定，萬(wàn)柏林和下元沉降減緩，吳家堡變化不大，小店年均沉降速率為45 mm/a的結(jié)果。2018年，許慧鵬[10]等人基于水準(zhǔn)測(cè)量對(duì)太原地面沉降進(jìn)行時(shí)空特征分析與模擬，小店地區(qū)沉降加劇。2016年，李軍的《太原市現(xiàn)今地面沉降特征分析》[11]中，GPS監(jiān)測(cè)點(diǎn)TY07為觀象臺(tái)，2009年7月～2014年6月的累計(jì)沉降量為-115.5 mm，年平均速率為23.1 mm/a。眾多論文數(shù)據(jù)或分析結(jié)果與本文分析結(jié)論一致，證明本文數(shù)據(jù)處理成果的可靠性。

4 結(jié) 論

采用PS-InSAR技術(shù)對(duì)覆蓋太原市區(qū)的37景2013年6月～2016年6月的COSMO-SkyMed數(shù)據(jù)進(jìn)行處理，得到該區(qū)域的形變速率圖和時(shí)間序列累計(jì)沉降量圖。結(jié)果顯示在該研究期內(nèi)太原共計(jì)形成萬(wàn)柏林、南堰、小店三個(gè)沉降中心，尤其注意的是小店地區(qū)沉降范圍逐漸擴(kuò)大，沉降量逐漸變大，且最大沉降量達(dá)-152.7 mm，年沉降速率達(dá)50.901 mm/a，而吳家堡、西張地區(qū)沉降現(xiàn)象較為穩(wěn)定。