星載InSAR技術(shù)支持下的昆明地表沉降監(jiān)測

麻源源，陳云波，左小清，麻衛(wèi)峰，吳文豪(. 昆明理工大學(xué)國土資源工程學(xué)院，云南 昆明 65009； . 昆明市規(guī)劃編制與信息中心，云南 昆明 650500； . 云南師范大學(xué)旅游與地理科學(xué)學(xué)院，云南 昆明 650000; . 湖南科技大學(xué)煤炭資源清潔利用與礦山環(huán)境保護湖南省重點實驗室，湖南 湘潭 0)

昆明地面沉降與高大建筑群的成片開發(fā)及地質(zhì)因素有關(guān)[1]。主要表現(xiàn)為在一定的范圍內(nèi)發(fā)生地面形變，形成多個沉降漏斗，影響了昆明城市規(guī)劃和發(fā)展。昆明市是目前為止中國大陸西部高原地區(qū)唯一一個發(fā)生沉降現(xiàn)象的城市[2]。合成孔徑雷達干涉測量(interferometry synthetic aperture rader，InSAR)技術(shù)，因大面積、非接觸、全天候、高精度等特點，與傳統(tǒng)的監(jiān)測測量(GPS和精密水準(zhǔn))技術(shù)相比有著一定的優(yōu)越性。在近幾年發(fā)展的時間序列差分干涉測量技術(shù)中，小基線集(small baseline subset，SBAS-InSAR)技術(shù)和永久散射體雷達干涉技術(shù)(persistent scatters InSAR，PS-InSAR)克服常規(guī)D-InSAR(differential InSAR)在地面形變監(jiān)測過程中受到失相干和大氣延遲的影響，可以進行長時間監(jiān)測及提高InSAR技術(shù)監(jiān)測精度。針對昆明沉降現(xiàn)象，國內(nèi)許多學(xué)者通過不同方法和技術(shù)進行了大量研究，文獻[3]通過二等水準(zhǔn)測量對昆明市區(qū)在1986—1987年期間的垂直沉降進行了分析。文獻[4]從礦物學(xué)和土質(zhì)結(jié)構(gòu)的角度對昆明在1987—1998年期間地面沉降進行了分析。文獻[5]利用歷史資料對1979—1998年期間昆明地面沉降發(fā)展過程及其特征進行了研究。文獻[6]利用ALOS PLASAR數(shù)據(jù)基于SBAS-InSAR技術(shù)對昆明市區(qū)在2008—2010年的沉降進行了分析。文獻[7]基于3S技術(shù)獲取高精度的昆明地面沉降信息。本文利用歐空局發(fā)布的Sentinel-1A數(shù)據(jù)，提出基于升降軌SBAS-InSAR技術(shù)和數(shù)據(jù)融合技術(shù)對昆明在2014年12月至2017年3月期間的地面沉降進行分析，試驗結(jié)果與上述學(xué)者所分析結(jié)果在形變趨勢上存在一致性。呈現(xiàn)出以小板橋、東菊新村、河尾村、曹家村、羅家村為中心的多個沉降漏斗，而且以河尾村、曹家村、東菊新村和羅家村為中心的沉降漏斗有連成一片的趨勢。其中沉降最嚴(yán)重的是以小板橋為中心的沉降漏斗，最大沉降速率達-54.2 mm/a。

1 小基線集(SBAS)技術(shù)原理

SBAS-InSAR技術(shù)先將時序SAR影像組成若干時空基線較短的干涉對集合，再利用奇異值分解方法將多個小基線集聯(lián)合起來求地表形變速率最小二乘解，有效地解決了不同InSAR數(shù)據(jù)集間空間基線過長造成時間不連續(xù)的問題。SBAS-InSAR技術(shù)處理SAR影像時，除影像相位中形變相位外，還有區(qū)域地形誤差帶來的相位，由獲取影像時大氣不均質(zhì)帶來的大氣相位異常、時空基線失相干現(xiàn)象和熱噪聲等引起相位誤差，在差分相位干涉圖中，兩個高相干點的觀測方程表示為[8]

(1)

2 SAR影像數(shù)據(jù)與處理

2.1 Sentinel-1數(shù)據(jù)來源

數(shù)據(jù)采用歐空局發(fā)布Sentinel-1A衛(wèi)星數(shù)據(jù)[10]，Sentinel-1A雷達衛(wèi)星于2014-04-03發(fā)射升空，是歐空局哥白尼計劃發(fā)射的首顆環(huán)境衛(wèi)星，該衛(wèi)星數(shù)據(jù)實施免費對外開放政策。Sentinel-1A數(shù)據(jù)在IW和EW模式下采用TOPS(terrain observation with progressive)模式，該模式控制天線旋轉(zhuǎn)方向與聚束模式正好相反，由后向前以固定速度旋轉(zhuǎn)，使得所有目標(biāo)在相同方位向天線模式下觀測，降低扇貝效應(yīng)和方位向變化模糊度[11-12]。Sentinel-1A數(shù)據(jù)具有精密軌道控制技術(shù)，采用交規(guī)盲區(qū)控制技術(shù)(across track dead-band control)確保任何重軌時衛(wèi)星每個軌道點與地固參考軌道的位置誤差均方根(RMS)在100 m范圍內(nèi)。試驗覆蓋區(qū)域如圖1所示。

圖1 試驗區(qū)域覆蓋示意圖

利用Sentinel-1A衛(wèi)星數(shù)據(jù)通過SBAS-InSAR技術(shù)可以對同一地區(qū)重訪周期內(nèi)地面沉降快速監(jiān)測，非常適合昆明城區(qū)的地表形變監(jiān)測，本文收集IW模式下Sentinel-1A數(shù)據(jù)，數(shù)據(jù)參數(shù)見表1。

表1 Sentinel-1數(shù)據(jù)參數(shù)

除此之外，衛(wèi)星軌道數(shù)據(jù)由歐空局提供精密軌道數(shù)據(jù)和美國宇航局提供的SRTM3 DEM，分辨率為90 m，精度為±30 m，用于去除地形相位。由于升降軌數(shù)據(jù)開始日期和結(jié)束日期均很接近，為升降軌模式下的觀測數(shù)據(jù)融合提供了可能性。

2.2 SBAS-InSAR主要流程

(1) 生成連接圖?；诳臻g基線閾值ΔB=200 m和時間基線閾值ΔT=550 d為基線條件，共選取升軌232對干涉數(shù)據(jù)和降軌227對干涉數(shù)據(jù)，其中降軌數(shù)據(jù)以獲取日期2015-10-16為主影像，升軌數(shù)據(jù)以獲取日期2015-10-02為主影像。Sentinel-1A數(shù)據(jù)升降軌干涉數(shù)據(jù)時空間基線分布如圖2所示。

圖2

(2) 影像配準(zhǔn)：將所有輔影像進行配準(zhǔn)并重采樣至主影像。Sentinel-1A衛(wèi)星采用TOPS模式，配準(zhǔn)精度要求達到千分之一像素級，相當(dāng)于方位向零多譜時間配準(zhǔn)真誤差小于2 μs，由于地形起伏會對配準(zhǔn)造成影響，多項式難以準(zhǔn)確地擬合出主副影像幾何偏移畸變的映射，因此利用幾何配準(zhǔn)誤差為系統(tǒng)誤差的特點采用哨兵衛(wèi)星精密軌道進行初配準(zhǔn)，再利用增強譜分集算法(enhanced spectral diversity，ESD)進行精配準(zhǔn)來消除片段間出現(xiàn)的相位跳變現(xiàn)象[13-14]。最后進行方位譜去斜，達到干涉需要的配準(zhǔn)精度獲取正確的干涉圖。

(3) 對干涉對進行集合D-InSAR處理：包括干涉對生成、去地平效應(yīng)、干涉圖濾波、相干性計算和相位解纏。干涉圖濾波采用Goldstein方法進行干涉圖濾波，相位解纏方法采用Delaunay三角網(wǎng)最小費用流量法(minimum cost flow，MCF)方法進行三維相位解纏時，升降軌γx均選取0.35作為閾值選取PS候選點。

(4) 軌道精煉和重去平：在參考DEM不準(zhǔn)確的情況下，需要用地面控制點(GCP)作為穩(wěn)定參考點，從這些點中計算出誤差相位從而去除并且修正SAR數(shù)據(jù)。使用同一組GCP數(shù)據(jù)修正升降軌干涉圖，目的是減少兩個圖像間相位誤差和使用不同GCP帶來的誤差。

(5) SBAS-InSAR數(shù)據(jù)反演。估算形變速率和殘余地形，基于線性模型計算所有像對形變。線性模型相比其他模型無需密集連接圖和高相干性即可得到可靠結(jié)果且穩(wěn)定性較好。計算時間序列上位移，在得到形變速率基礎(chǔ)上，本文通過時域高通濾波和頻域低通濾波去除殘余DEM誤差、大氣干擾與軌道誤差影響。

(6) SAR數(shù)據(jù)地理編碼和輻射定標(biāo)：將上述結(jié)果編碼到GCS-WGS-84坐標(biāo)系。

2.3 數(shù)據(jù)融合簡化模型

升降軌SBAS-InSAR形變速率可以區(qū)分垂直方向和水平方向的形變速率。由于升降軌兩組觀測數(shù)據(jù)幾乎從相反的兩個方向獲得。對于地表相同目標(biāo)而言，東西方向水平形變速率在兩組觀測數(shù)據(jù)中的表現(xiàn)是基本相反的變化，但是垂直方向上則表現(xiàn)為基本相同形變速率變化。視線向上的形變速率量(DLos)的構(gòu)成為

DLos=(UNsinφ-UEcosφ)sinθ+UVcosθ

(2)

式中，UN、UE和UV分別為南北、東西和垂直方向地表形變速率分量；φ為升降軌方位角；θ為入射角。哨兵-1雷達衛(wèi)星的干涉寬幅模式(IW)，升軌入射角θ約為40°，降軌入射角θ約為39.5°。雷達視線方向(line of sight，LOS)形變對UN方向形變不敏感，視線向上的形變速率量表達式可進一步簡化為

DLos≈±UNsinθ+UVcosθ

(3)

根據(jù)文獻[12]可知，地表變形以垂直方向上形變速率為主，水平方向上的形變速率極為微小，即UNsinθ?DLos。由于入射角在研究區(qū)的變化微弱，因此可用均值來估計地表垂直方向上的形變速率分量。由文獻[15]可知入射角變化的簡化對方位角的變化影響微弱，可以忽略不計。因而，可以推出升降軌觀測模式下垂直方向的形變速率量公式[15-16]為

(4)

(5)

利用式(4)可將升降軌高相干點上雷達視線方向(LOS)形變轉(zhuǎn)換為垂直方向地表形變速率分量UVA(D)，再根據(jù)式(5)求得升降軌數(shù)據(jù)融合后的垂直方向上速率分量UV，其前提是要求升降軌模式下同一觀測區(qū)域具有相同的觀測區(qū)域。升降軌下的昆明垂直沉降速率如圖3所示。

圖3

從圖3(a)和圖3(b)可知：兩種模式下的沉降漏斗中心位置、漏斗形狀、沉降區(qū)域大小及沉降速率大小幾乎一致，均在-55～10 mm/a之間。

3 升降軌數(shù)據(jù)相互驗證和數(shù)據(jù)融合

3.1 升降軌數(shù)據(jù)相互驗證

由于雷達波入射角和入射方向的影響,同一相干目標(biāo)在升降軌SAR觀測下對象的位置略有差異[17-18]，為了對升降軌觀測值進行相互對比驗證，本文采用地形校正方法統(tǒng)一升降軌坐標(biāo)系消除基準(zhǔn)偏差。依據(jù)單個點目標(biāo)進行直接比較難以實現(xiàn)，利用統(tǒng)計方法進行對比分析，從整體上驗證升降軌在地表形變垂直形變的相關(guān)性。從升降軌地表垂直沉降速率圖中提取11 072個同名點，其兩組數(shù)據(jù)計算得到的昆明升降軌垂直沉降相關(guān)分析如圖4所示。

圖4 昆明升降軌垂直沉降相關(guān)分析圖

由圖4得知，兩組數(shù)據(jù)線性擬合的方程系數(shù)K及兩組數(shù)據(jù)相干性計算得到判定系數(shù)R×R、相關(guān)系數(shù)coherence及均方差根誤差(RMES)[19]，見表2。

表2 升降軌垂直沉降速率關(guān)系

除此之外，圖4中擬合的線性方程Y=0.962 5X-0.799 5與理想狀態(tài)下的直線y=x方程對比可知兩直線非常接近。因此試驗結(jié)果表明兩組數(shù)據(jù)具有較高的一致性和相關(guān)性，當(dāng)SAR影像為30景以上，升降軌垂直沉降速率監(jiān)測互檢驗精度優(yōu)于2 mm/a。

3.2 升降軌數(shù)據(jù)數(shù)據(jù)融合

為了充分利用升降軌SBAS-InSAR技術(shù)監(jiān)測結(jié)果和避免因雷達疊掩、透視收縮、陰影等引起失真現(xiàn)象[20]，本文對垂直沉降速率通過式(5)進行數(shù)據(jù)均值融合[21-22]。目的是提高昆明地表沉降監(jiān)測精度和可靠性。圖5為融合后的2014年12月至2017年3月昆明垂直沉降速率圖。

圖5 2014年12月至2017年3月期間昆明垂直沉降速率序列圖

4 沉降分析

由圖5分析可知昆明城區(qū)2014—2017年期間地面沉降分布的基本情況，昆明五華區(qū)沉降較慢，相對穩(wěn)定，而昆明西山區(qū)和官渡區(qū)內(nèi)的滇池北岸附近區(qū)域沉降較快，形成多個沉降漏斗。由于缺少同期的水準(zhǔn)測量和GPS測量數(shù)據(jù)，無法準(zhǔn)確地評價本文監(jiān)測數(shù)據(jù)的精度，本文將對比不同歷史時期[2-5]昆明城區(qū)的沉降速率及對昆明城區(qū)地表沉降的演變規(guī)律進行分析，不同時期昆明城區(qū)地表垂直沉降見表3。

表3 不同時期昆明城區(qū)地表垂直沉降 mm/a

試驗結(jié)果分析：SBAS-InSAR技術(shù)監(jiān)測昆明城區(qū)地面沉降分布與歷史沉降監(jiān)測數(shù)據(jù)得到的沉降漏斗分布是基本一致的，但也有新的情況變化。2014—2017年間沉降基本情況是五華區(qū)內(nèi)的嚴(yán)家山一帶和吳家營-大塘子一帶，該地區(qū)較為穩(wěn)定，與歷史時期的沉降速率相比變化不大，沉降速率保持在-13 mm/a以內(nèi)。其主要原因是該地區(qū)沉積物主要以巖石為主，其壓縮性很小，地面形變量不大。而昆明地面沉降嚴(yán)重區(qū)域主要分布在官渡區(qū)和西山區(qū)，其中靠近滇池的區(qū)域出現(xiàn)大面積沉降，形成多個漏斗中心，形成以河尾村、曹家村、小板橋、羅家村和東菊新村為中心的5個沉降漏斗，5個沉降中心垂直沉降速率均達到-20 mm/a以上且有連成一片的趨勢。并且羅家村-曹家村一帶在歷史時期并沒有出現(xiàn)沉降現(xiàn)象，而在2014—2017年間卻形成了兩個沉降漏斗，沉降速率達到-25 mm/a以上。在幾個沉降區(qū)域中，小板橋一帶的沉降最明顯,其最大平均地表沉降速率達-54.2 mm/a。說明目前昆明滇池附近和小板橋沉降情況較為嚴(yán)重，并且嚴(yán)重性逐漸增大，嚴(yán)重區(qū)域正在不斷擴大。主要原因與該地區(qū)的地質(zhì)因素有很大的關(guān)系，昆明地區(qū)的地基沉降主要以湖濱地帶軟土沉降為主，近年來滇池周圍大量施工，高層建筑加壓會導(dǎo)致軟土的水被擠出，土層被壓密實。加上工業(yè)用水劇增，導(dǎo)致承壓水位降低，地層向更為致密狀態(tài)變化而引起地基下沉。對滇池附近進行實地調(diào)查發(fā)現(xiàn)，滇池北岸周圍存在道路地基錯位、房屋傾斜，房體與地面分離、屋內(nèi)地面下沉等現(xiàn)象，如圖6所示。

圖6 外業(yè)實地核查照片

外業(yè)實地調(diào)查結(jié)果表明，昆明地面沉降已經(jīng)危害到公共安全和人民的居住安全。該調(diào)查也驗證了InSAR技術(shù)可以快速、直觀、有效地發(fā)現(xiàn)地面沉降嚴(yán)重區(qū)域及沉降中心，說明InSAR技術(shù)用于城市地面形變監(jiān)測是可行的。

5 結(jié) 語

本文利用Sentinel-1A數(shù)據(jù)，基于升降軌SBAS-InSAR技術(shù)監(jiān)測昆明城區(qū)地表沉降情況，獲取升降軌兩種不同模式下的昆明城區(qū)2014—2017年間地表沉降速率及城區(qū)的沉降漏斗分布情況，通過對兩組觀測數(shù)據(jù)進行相互驗證表明衛(wèi)星兩種飛行模式在地面監(jiān)測具有一致性。由于獲取升降軌SAR數(shù)據(jù)的起始和結(jié)束時間十分接近，具有相近的觀測時間段，因此對兩種數(shù)據(jù)均值融合可以避免衛(wèi)星單一軌道下觀測數(shù)據(jù)失真。總體而言，SBAS-InSAR技術(shù)監(jiān)測地表信息效果理想，與傳統(tǒng)形變監(jiān)測技術(shù)相比可以在短時間內(nèi)得到大面積、高空間密度的地表形變場，具有成本較低無需布設(shè)地面控制點等優(yōu)點。

參考文獻：

[1] 郭樂萍，岳建平，岳順.SBAS技術(shù)在南京河西地表沉降監(jiān)測中的應(yīng)用[J].測繪通報，2017(3):26-28.

[2] 姜朝松，邵德晟，樊友心,等.昆明市地面沉降發(fā)展過程及其特征[J].地震研究，2001，24(1):55-60.

[3] 黃崐，萬建平，吳彬生.昆明地區(qū)近期地面垂直升降運動及分析[J].地震研究，1989，12(3):249-253.

[4] 薛傳東，劉星，李保珠,等.昆明市區(qū)地面沉降的機理分析[J].中國地質(zhì)災(zāi)害與防治學(xué)報，2004，15(3):47-54.

[5] 孟國濤.昆明南市區(qū)地面沉降研究[D].昆明：昆明理工大學(xué)，2003.

[6] 尹振興，鐘麗云，許兵,等.基于SBAS-InSAR的昆明地面沉降監(jiān)測研究[J].地礦測繪，2016，32(4):1-5.

[7] 高照忠，羅志清，魏海霞.3S技術(shù)在昆明市地面沉降監(jiān)測中的應(yīng)用[J].地礦測繪，2007，23(1):27-29.

[8] BERARDINO P，F(xiàn)ORNARO G，LANARI R，et al.A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J].IEEE Transactions on Geoscience and Remote Sensing，2003，40(11):2375-2383.

[9] 周志偉，鄢子平，劉蘇,等.永久散射體與短基線雷達干涉測量在城市地表形變中的應(yīng)用[J].武漢大學(xué)學(xué)報(信息科學(xué)版)，2011，36(8):928-931.

[10] 張艷梅，王萍，羅想,等.利用Sentinel-1數(shù)據(jù)和SBAS-InSAR技術(shù)監(jiān)測西安地表沉降[J].測繪通報，2017(4):93-97.

[11] 吳文豪，周志偉，李陶，等.精密軌道支持下的哨兵衛(wèi)星TOPS模式干涉處理[J].測繪學(xué)報，2017，46(9)：1156-1164.

[12] SOWTER A,AMAT MBC,CIGNA F,et al.Mexico City Land Subsidence in 2014—2015 with Sentinel-1 IW TOPS：Results Using the Intermittent SBAS (ISBAS) Technique[J].International Journal of Applied Earth Observations and Geoinformation,2016,52：230-242.

[13] 吳文豪.哨兵雷達衛(wèi)星TOPS模式干涉處理研究[D].武漢：武漢大學(xué),2016.

[14] 白澤朝，靳國旺，張紅敏，等.天津地區(qū)Sentinel-1A雷達影像PSInSAR地面沉降監(jiān)測[J].測繪科學(xué)技術(shù)學(xué)報,2017,34(3):283-288.

[15] 王濤，顧麗娟，詹華明,等.基于D-InSAR技術(shù)的天津地區(qū)地面沉降監(jiān)測[J].測繪科學(xué)，2013，38(6)：23-28.

[16] 王艷，張玲，葛大慶,等.升降軌PSInSAR觀測反演沉降與水平向位移試驗[J].國土資源遙感，2014，26(4):97-102.

[17] 劉波，晏王波，范雪婷,等.基于SBAS技術(shù)的武進區(qū)地面沉降監(jiān)測應(yīng)用研究[J].現(xiàn)代測繪，2016，39(4):9-12.

[18] 王艷，葛大慶，張玲,等.升降軌PSInSAR地面沉降監(jiān)測結(jié)果的互檢驗與時序融合[J].國土資源遙感，2014，26(4):125-130.

[19] 劉一霖.黃河三角洲地面沉降時序InSAR技術(shù)監(jiān)測與地下流體開采相關(guān)性分析[D].青島：中國科學(xué)院研究生院(海洋研究所)，2016.

[20] 丁瓊，劉國祥，蔡國林,等.InSAR DEM精度與地形特征的關(guān)系分析[J].測繪科學(xué)，2009，34(1):147-148.

[21] 王琴，陳蜜，劉書軍,等.利用升降軌道SAR數(shù)據(jù)獲取DEM的試驗研究[J].測繪通報，2015(6):39-43.

[22] 王霞迎，張菊清，張勤,等.升降軌InSAR與GPS數(shù)據(jù)集成反演西安形變場[J].測繪學(xué)報，2016，45(7):810-817.