應(yīng)用升降軌SBAS-InSAR技術(shù)監(jiān)測昆明市地面沉降

范 軍,左小清,李 濤,陳乾福,麻源源

(1. 昆明理工大學(xué) 國土資源工程學(xué)院,云南 昆明 650093;2. 國家測繪地理信息局 衛(wèi)星測繪應(yīng)用中心,北京 100830)

0 引 言

地面沉降是一種十分嚴(yán)重的地災(zāi)現(xiàn)象,與人們的生活息息相關(guān).從2001~2017年,昆明市在西南地區(qū)無論是城市化建設(shè)還是軌道交通的飛速發(fā)展都已進(jìn)入先進(jìn)行列,但不可避免地導(dǎo)致地面發(fā)生嚴(yán)重的沉降.目前根據(jù)文獻(xiàn)顯示,2001年姜朝松等人利用1979~1998年水準(zhǔn)數(shù)據(jù)對昆明地面沉降的發(fā)展過程及其特征進(jìn)行研究[1],表明小板橋和河尾村出現(xiàn)兩個(gè)沉降漏斗;2003年孟國濤利用1987~1998年四期水準(zhǔn)數(shù)據(jù)對昆明南市區(qū)地面沉降進(jìn)行研究[2],表明廣衛(wèi)村-小板橋-河尾村-六甲和漁戶村-福海出現(xiàn)嚴(yán)重的沉降漏斗;2004年薛傳東在此基礎(chǔ)上對昆明市區(qū)地面沉降進(jìn)行研究[3],表明小板橋和河尾村沉降明顯加快,漁戶村-河尾村-大塘子-嚴(yán)家山等4處漏斗已連成一體.以上3位都是基于精密水準(zhǔn)測量方法監(jiān)測地面沉降,這些方法雖然精度滿足要求,但是存在工作力度大、周期長、成本高、效率低、空間分辨率低和覆蓋范圍小等缺點(diǎn),也很難滿足時(shí)間、空間尺度的需要.2007年高照中結(jié)合"3S"技術(shù)建立自動(dòng)監(jiān)測系統(tǒng)對昆明市地面進(jìn)行高精度沉降監(jiān)測[4],對比傳統(tǒng)水準(zhǔn)測量,GPS具有覆蓋范圍廣、精度高、速度快等優(yōu)點(diǎn),但作業(yè)周期長、測量結(jié)果不連續(xù)難以高效地監(jiān)測與分析整個(gè)城市的區(qū)域性沉降.合成孔徑雷達(dá)差分干涉測量(DInSAR)是近年來發(fā)展起來的一種監(jiān)測地表沉降的有效手段,能夠?qū)Φ乇磉M(jìn)行全天候、大范圍、亞厘米級精度的監(jiān)測,但是傳統(tǒng)DInSAR易受時(shí)空失相干與大氣效應(yīng)的影響,較難完成高精度長時(shí)間間隔的地表監(jiān)測.短基線干涉測量技術(shù)(SBAS)的提出,能夠充分有效地利用多景雷達(dá)影像,利用時(shí)序分析技術(shù)有效解決傳統(tǒng)DInSAR技術(shù)中一些問題,提高DInSAR地面監(jiān)測的精度[5-6].目前,SBAS技術(shù)在城市地面監(jiān)測已經(jīng)得到廣泛應(yīng)用.2016年尹振興選取2007~2010年SAR數(shù)據(jù)基于SBAS-InSAR技術(shù)對昆明地面沉降進(jìn)行監(jiān)測[7],得出小板橋、河尾村等沉降漏斗沉降速率趨于穩(wěn)定,并出現(xiàn)多處新的沉降漏斗.但是至今,國內(nèi)仍然很少有人對當(dāng)前時(shí)間的昆明市地面沉降情況進(jìn)行研究.因此,本文選取2014~2017年29景升軌和32景降軌的Sentinel-1A、1B數(shù)據(jù),基于SBAS技術(shù)對昆明市地面沉降進(jìn)行研究.

1 SBAS-InSAR技術(shù)原理

SBAS方法是由Berardino和Lanari等研究人員于2002年正式提出的一種InSAR時(shí)間序列分析方法[8].該方法通過選取數(shù)據(jù)自由組合所有短基線干涉對,得到一系列短基線差分干涉圖,這些差分干涉圖能夠有效地克服空間去相關(guān)現(xiàn)象,不僅連接所有差分中短基線集,同時(shí)連接較大空間基線中分開的孤立SAR數(shù)據(jù)集,來提高觀測數(shù)據(jù)時(shí)間采樣率.在求解形變速率時(shí),SBAS方法采用奇異值分解(SVD)法來獲得高相干點(diǎn)的沉降速率和時(shí)間序列.其過程主要原理如下:

獲取同一研究區(qū)N+1幅SAR影像,按照影像獲取的時(shí)間順序(t0,t1,…,tN)排列,選取其中一幅影像為主影像,其他SAR影像配準(zhǔn)到這幅主影像上,在短基線干涉組合條件限制下,生成M幅差分干涉圖,且M滿足:

假設(shè)對于任意干涉圖j在tA與tB時(shí)間內(nèi)獲取的影像進(jìn)行干涉生成(tB>tA),方位向坐標(biāo)x和距離向坐標(biāo)r的像素的干涉相位可以表示為:

式中,j∈(1,…,M);λ為雷達(dá)信號的波長;φB(x,r)和φA(x,r)分別為tB與tA時(shí)刻SAR影像的相位值;d(tB,x,r)和d(tA,x,r)分別為tB與tA時(shí)刻相對于d(t0,x,r)=0的雷達(dá)視線方向的累積形變量值;φjtopo(x,r)為差分干涉圖中殘余的地形相位誤差,主要是由于參考DEM精度較低而導(dǎo)致的,如果殘余地形相位較小,可在解算過程中忽略;φjnoise(x,r)為噪聲相位;φjatm(tB,tA,x,r)為大氣相位誤差.如果不考慮殘余地形相位、噪聲相位和大氣相位誤差,那么式(2)可表示為:

為了更好地表示地面沉降時(shí)間序列,可將式(3)相位值表示為兩個(gè)獲取時(shí)間段的平均相位速率Vj與該時(shí)間段之間的乘積.

則第j幅干涉圖的相位值可表示為:

即第j幅干涉圖的相位值表示為各時(shí)間段平均速率在主、從影像時(shí)間間隔上的積分,可表示為新的矩陣方程:

式(6)中,B是一個(gè)MXN的系數(shù)矩陣,由于SBAS采用多個(gè)主影像時(shí)空基線獲取干涉像對的方法,因此矩陣B產(chǎn)生秩虧,運(yùn)用奇異值分解法(SVD)對矩陣B進(jìn)行解算,得到矩陣B的廣義逆矩陣,進(jìn)而得到速率的最小范數(shù)解,最后在各個(gè)時(shí)間段對速率積分就可得到各個(gè)時(shí)間段的累積沉降量[9].

2 實(shí)驗(yàn)區(qū)概況及實(shí)驗(yàn)數(shù)據(jù)

昆明市區(qū)坐落于昆明晚新生代斷陷盆地內(nèi),廣泛發(fā)育第四紀(jì)松散沉積層,以湖沼相粉砂和軟粘土為主,夾多層淤泥和草煤層等松軟沉積物[3].近年來,由于昆明市城市與軌道交通建設(shè)的飛速發(fā)展,市區(qū)發(fā)生大規(guī)模地面沉降現(xiàn)象,沉降范圍日益擴(kuò)大,新的沉降中心不斷產(chǎn)生.為了進(jìn)一步研究昆明市地面沉降,本文研究區(qū)域選擇為升降軌模式下Sentinel-1A、1B(簡稱S1A、S1B)數(shù)據(jù)的昆明市區(qū),對升降軌實(shí)驗(yàn)數(shù)據(jù)進(jìn)行裁剪得到本文研究區(qū)域,如圖1所示紅色表示升軌,黑色表示降軌,中心經(jīng)緯度為25.00°N,102.71°E,覆蓋面積約為400 km2,該研究區(qū)主要覆蓋官渡區(qū),五華區(qū),盤龍區(qū),西山區(qū)等部分區(qū)域.

圖1 研究區(qū)地理位置Fig.1 Geographical location of the study area

本文從歐空局共獲取2014-12-06~2017-02-05間的29景升軌(Track-26)和2014-10-09~2017-01-20間的32景降軌(Track-62)的S1A、S1B的Level-0級原始影像數(shù)據(jù),S1A在2014年發(fā)射,S1B在2016年發(fā)射.在升降軌模式下入射角分別為39.5°和39.6°,C波段(波長5.6 cm),分辨率5mX20m,觀測模式為干涉寬幅(IW),極化方式為單極化(VV)[10-12].具體影像數(shù)據(jù)及相關(guān)參數(shù)見表1.

表1 升降軌Sentinel-1影像數(shù)據(jù)Tab.1 Sentinel-1 image data for ascending and descending

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

對升降軌Level-1級單視復(fù)數(shù)影像(SLC)數(shù)據(jù)進(jìn)行裁剪得到本文研究區(qū)域,綜合考慮時(shí)間基線、空間基線、多普勒基線3種的臨界值,如圖2a升軌選取2015年8月15日為超級主影像,并對所有影像進(jìn)行4X1(距離向X方位向)的多視處理.試驗(yàn)中設(shè)置的時(shí)間基線和空間基線閾值分別為550 d和±250 m,將所有的配對的干涉像對都配準(zhǔn)到主影像上,進(jìn)行多主影像干涉對組合,剔除質(zhì)量較差的干涉對,最終選取361個(gè)具有較高質(zhì)量的差分干涉對如圖2a,黃色代表超級主影像,綠色代表SAR影像,線段代表干涉對.其中最大時(shí)間基線540d,最小時(shí)間基線12 d,平均時(shí)間基線為235d,最大空間基線為227.718m,最小空間基線為3.028m,平均空間基線為55.561m.如圖2b降軌選取2015年8月5日為超級主影像,最終選取426個(gè)差分干涉對,其中最大時(shí)間基線540d,最小時(shí)間基線24d,平均時(shí)間基線為240d,最大空間基線為212.352 m,最小空間基線為1.748m,平均空間基線為54.370 m.由于升軌和降軌的空間基線分布較小、相干性較穩(wěn)定,所構(gòu)集合網(wǎng)絡(luò)更強(qiáng),保證后續(xù)計(jì)算結(jié)果的可靠性[13].

圖2 時(shí)空基線分布Fig.2 Spatial and temporal baseline distribution

對于本次試驗(yàn),在數(shù)據(jù)處理過程中,采用美國宇航局NASA提供的分辨率為90m的SRTM數(shù)據(jù)去除地形相位和平地效應(yīng)的影響,升降軌相干系數(shù)閾值設(shè)置為0.35(相干系數(shù)小于該閾值的像元不進(jìn)行解纏)采用Delaunay MCF方法進(jìn)行3D相位解纏得到差分干涉相位[14];采用Goldstein方法對干涉圖進(jìn)行濾波處理,提高干涉圖信號質(zhì)量,減少由空間基線和時(shí)間基線引起的失相干噪聲[15];利用歐空局(ESA)發(fā)布的精密軌道數(shù)據(jù)提高影像的軌道精度,但是在生成干涉圖中仍然存在軌道的殘余相位誤差.本文采用自動(dòng)精煉的方法,利用GCP對所有數(shù)據(jù)進(jìn)行重去平;為了對生成干涉圖進(jìn)行優(yōu)化,采用同樣方法進(jìn)行二次解纏;通過大氣高通或低通濾波去除大氣相位;最后將升降軌下SBAS結(jié)果分別進(jìn)行地理編碼,使不同軌道獲取的形變觀測值處于相同地理坐標(biāo)系下.采用奇異值分解法(SVD)法求解高相干點(diǎn)時(shí)間序列的沉降速率,最后通過積分獲得各個(gè)時(shí)段的累積形變量[10-11].

4 結(jié)果分析

4.1 地面沉降結(jié)果分析

經(jīng)過上述SBAS處理,在升降軌模式下得到2014年10月~2017年2月昆明市區(qū)地面沉降平均速率圖(底圖為昆明市Google地圖)如圖3所示,本次沉降量選取3個(gè)相對較大的矩形區(qū)域?yàn)锳、B、C,并計(jì)算各個(gè)區(qū)域內(nèi)地面沉降平均速率(沿雷達(dá)視線方向),紅色代表遠(yuǎn)離衛(wèi)星,表示地面下沉,藍(lán)色代表接近衛(wèi)星,表示地表上升[16].

圖3 升(左)降(右)軌模式下研究區(qū)年平均沉降速率Fig.3 Annual average subsidence rate of the study area in ascending(left:Track-26) and descending (right:Track-62)

A區(qū)域位于官渡區(qū),由升降軌模式下沉測量結(jié)果圖4所示,沉降漏斗主要位于滇池北岸以及東北岸,分布范圍廣泛,出現(xiàn)多處沉降漏斗,其中沉降最嚴(yán)重為彼岸小區(qū)-楓林盛景小區(qū)-星宇園小區(qū)-子君欣景小區(qū)-羊甫村-義路村-廣衛(wèi)村-螺螄灣北等多處沉降漏斗連成一體,構(gòu)成研究區(qū)域內(nèi)最大的漏斗沉降區(qū)域,范圍有明顯加劇趨勢,沉降漏斗形成面積已達(dá)到5.63km2,升降軌模式最大年沉降速率分別為-37.405mm/a和-38.975mm/a,最大沉降量達(dá)到89 mm.沉降原因主要由于居民地和交通網(wǎng)絡(luò)密切集中.由于昆明地區(qū)的地面沉降主要以軟土沉降為主,因此地面受到建筑用地施壓造成軟粘土中的水被擠出,土層被壓密實(shí),導(dǎo)致地基下沉;其次,這些沉降區(qū)域位于昆磨高速公路和東環(huán)城高速公路交界的廣衛(wèi)立交橋的西北處附近,此處有正在建設(shè)的昆明地鐵4號線,其中螺螄灣北為4號地鐵站,年沉降速率分別為-31.119 mm/a和-32.483 mm/a,這些區(qū)域受到地鐵隧道開挖和相應(yīng)的工程降水措施會直接擾動(dòng)軟土,同時(shí)地面上下車荷載交通流量引起的土體變形和地層移動(dòng)可能影響地表結(jié)構(gòu),共同觸發(fā)該路段以及其周邊區(qū)域的沉降[17].小板橋位于云大西路和昆磨高速公路與廣昆高速公路交界的鳴泉村立交橋西南處附近,年沉降速率為-20.174 mm/a和-21.431 mm/a,沉降原因不僅受到地質(zhì)因素的影響,同時(shí)還受到相對集中的居民地和交通流量荷載的影響,造成該區(qū)域沉降.曉東村,雨龍村,小板橋,竹園村和陳旗營等沉降中心靠近地鐵1號線已連成一片,曉東村為1號地鐵站,在沉降區(qū)域表現(xiàn)最為嚴(yán)重,年沉降速率分別為-23.686mm/a和-23.626mm/a.

圖4 升(左)降(右)軌模式下A區(qū)域年平均沉降速率Fig.4 Annual average subsidence rate of A area in ascending (left)and descending (right)

B區(qū)域位于官渡區(qū),為昆明老城區(qū),靠近滇池,由監(jiān)測結(jié)果圖5可知,福保路附近一帶且出現(xiàn)多處沉降漏斗,沉降區(qū)域有明顯向市中心擴(kuò)張趨勢,從南向北分別為丁姚村,韓家村,龔家村,金家村,六甲村,葉家村,羅家村等沉降區(qū)域,其中六甲村出現(xiàn)多處沉降漏斗,韓家村-金家村-龔家村等沉降區(qū)域連成一片產(chǎn)生多處新的沉降中心,構(gòu)成研究區(qū)較為嚴(yán)重的漏斗沉降區(qū).在升降軌模式下,最大年沉降速率分別達(dá)到-25.381 mm/a與-25.809 mm/a,最大累積沉降量達(dá)到-59mm,該地區(qū)沉降主要因?yàn)轫n家村,金家村,龔家村等區(qū)域靠近滇池,滇池湖水面是盆地最低的徑流基準(zhǔn)面,大量地下水向滇池方向徑流和排泄;其次滇池附近地區(qū)主要存在大量軟土、飽水粉砂層,并夾有多層淤泥、草煤層等松軟沉積物,這些地層形成時(shí)間短、固結(jié)度低、孔隙率高、含水量大,導(dǎo)致滇池周邊地區(qū)地基太軟,隨著軟土下沉而出現(xiàn)明顯的沉降現(xiàn)象[2].此外,由于B區(qū)域位于老城區(qū),從理論上分析,一般可能是鄰近地區(qū)有人類工程活動(dòng)影響,如深基坑開挖和地下水超采、城市建筑加載等.

圖5 升(左)降(右)軌模式下B區(qū)域年平均沉降速率Fig.5 Annual average subsidence rate of B area in ascending (left)and descending (right)

C區(qū)域位于西山區(qū),由監(jiān)測結(jié)果圖6可知,沉降漏斗出現(xiàn)在春苑小區(qū)和棕樹營兩個(gè)區(qū)域,在升降軌模式下,最大年沉降速率分別為-20.944mm/a和-19.669 mm/a,而其他大部分區(qū)域無明顯的升降變化,則認(rèn)為地面相對穩(wěn)定.

圖6 升(左)降(右)軌模式下C區(qū)域年平均沉降速率Fig.6 Annual average subsidence rate of C area in ascending (left)and descending (right)

4.2 升降軌結(jié)果分析對比

由于文中缺少同期的測量數(shù)據(jù),通過采用點(diǎn)-點(diǎn)的驗(yàn)證方法,對升降軌模式下沉降漏斗中同名點(diǎn)的沉降速率進(jìn)行精度驗(yàn)證[18].首先分別在升降軌沉降速率圖上提取同名點(diǎn),同名點(diǎn)的獲取原則是采用點(diǎn)位追蹤算法在具有明顯的沉降漏斗處選取均勻分布的高相干點(diǎn).這些同名點(diǎn)具有相同的地理坐標(biāo),其中每一個(gè)同名點(diǎn)的沉降速率代表該沉降漏斗區(qū)域的沉降速率.兩者在相同時(shí)間段內(nèi)獲取的地面沉降速率對比見表2.

表2 升降軌模式下同名點(diǎn)沉降速率對比表Tab.2 Comparison table of subsidence rate under the same point in ascending and descending

從上表可以看出,升降軌模式下同名點(diǎn)的沉降速率整體都保持了較好的一致性,大部分同名點(diǎn)位獲取的沉降速率差值都在2mm/a以內(nèi),中誤差為±0.896mm/a,圖7為升降軌同名點(diǎn)沉降速率線性回歸分析圖[19].

圖7 升降軌同名點(diǎn)沉降速率線性回歸分析Fig.7 Linear regression analysis of average subsidence rate on the same point of ascending and descending

由圖7可知,升降軌模式下研究區(qū)同名點(diǎn)沉降速率的整體分布符合線性函數(shù)特征.兩者的相關(guān)系數(shù)達(dá)到了0.988,由于升降軌結(jié)果相互獨(dú)立,兩者顯示地面同名點(diǎn)沉降速率不僅一致,而且幅度也基本吻合,表明升降軌監(jiān)測結(jié)果精度的可靠性[20].圖8表示升軌與降軌模式下A區(qū)域和B區(qū)域時(shí)間序列的最大累積沉降量.

圖8 升降軌模式下A、B區(qū)域時(shí)間序列最大累積沉降量Fig.8 The maximum cumulative subsidence of A and B regional time series in ascending and descending

5 結(jié)束語

由于本文沒有收集到同期的實(shí)測水準(zhǔn)數(shù)據(jù),只是通過單獨(dú)采用29景升軌和32景降軌的S1A、S1B數(shù)據(jù)利用SBAS-InSAR技術(shù)分別對昆明市區(qū)的地面沉降進(jìn)行監(jiān)測,通過對比升降軌監(jiān)測結(jié)果和精度驗(yàn)證,可知在兩種不同類型觀測數(shù)據(jù)模式下,沉降結(jié)果基本吻合保持一致,表明了本文方法的可靠性.試驗(yàn)證明該研究區(qū)共出現(xiàn)3處較大的形變區(qū),其中A區(qū)域螺螄灣北,小板橋南一帶地面沉降基本呈線性沉降規(guī)律,最大年平均沉降速率為-38.975 mm/a,最大累積沉降量達(dá)到89 mm.B區(qū)域靠近滇池一帶,福保路周邊多處區(qū)域出現(xiàn)沉降漏斗連成一體,最大年平均沉降速率為-25.809mm/a,最大累積沉降量達(dá)到59mm,這些嚴(yán)重的沉降現(xiàn)象不僅與工程地質(zhì)環(huán)境密切相關(guān),近幾年來昆明主城區(qū)城市化和軌道交通建設(shè)的飛快發(fā)展,居民區(qū)的興建與密集的道路網(wǎng)增大地面荷載力,地鐵施工造成地下降水量增大,導(dǎo)致廣泛分布的粘松散型軟土層容易受到外力作用產(chǎn)生地面沉降,成為昆明市地面沉降主要原因.由于本文升降軌SAR數(shù)據(jù)均為S1A、S1B數(shù)據(jù),對于試驗(yàn)結(jié)果可能會存在偶然性,下一步工作是要采用其他SAR數(shù)據(jù)參與SBAS-InSAR技術(shù)監(jiān)測地面沉降,并開展野外水準(zhǔn)測量工作與InSAR處理結(jié)果進(jìn)行對比.