朱邦彥，李建成，儲征偉，王 曉，孫靜雯，姚馮宇

(1. 南京市測繪勘察研究院股份有限公司，江蘇 南京 210019；2. 武漢大學(xué)測繪學(xué)院，湖北 武漢 430079； 3. 淮海工學(xué)院測繪與海洋信息學(xué)院，江蘇 連云港 222000)

20世紀(jì)90年代中后期以來，基于時序InSAR(MT-InSAR)分析方法的高精度差分干涉測量技術(shù)得到發(fā)展并成為當(dāng)前研究的熱點。MT-InSAR技術(shù)是通過多幅干涉圖按照一定的時間序列建模以去除傳統(tǒng)InSAR技術(shù)中相位觀測誤差、DEM誤差、軌道誤差等影響，通過統(tǒng)計分析時間序列上幅度和相位信息的穩(wěn)定性，探測不受時空去相關(guān)影響的穩(wěn)定點目標(biāo)(permanent scatterer，PS)。經(jīng)典的如Ferretti等在2001年提出的基于永久散射體的干涉測量技術(shù)(PSInSAR)[1]。Berardino等在2002年首次提出小基線集(small baseline subset，SBAS)技術(shù)，通過選取短時空基線SAR影像對進(jìn)行干涉組合，從而提高形變信息以提取密度和精度[2]。MT-InSAR技術(shù)已經(jīng)在城市及其基礎(chǔ)設(shè)施沉降監(jiān)測中得到越來越廣泛的應(yīng)用[3-9]。

隨著社會經(jīng)濟(jì)的快速發(fā)展，我國的地面沉降現(xiàn)象不斷加劇。截至2009年，全國超過50個城市發(fā)生了較嚴(yán)重的地面沉降。地面沉降已成為影響我國生態(tài)文明建設(shè)和可持續(xù)發(fā)展的一個重大問題。截止到2016年底，寧波市區(qū)沉降面積已達(dá)500 km2，快速發(fā)展建設(shè)中出現(xiàn)的大面積地面填土堆載和交通動負(fù)荷，成為現(xiàn)階段誘發(fā)寧波市區(qū)地面沉降的主因[10]。為滿足地面沉降高時空分辨率信息快速獲取的迫切需求，本文基于Sentinel-1影像對寧波地區(qū)地面沉降進(jìn)行研究。

Sentinel-1A是歐洲委員會哥白尼環(huán)境監(jiān)測計劃于2014年發(fā)射的第一顆衛(wèi)星[11]。寬幅影像幅度達(dá)250 km，12 d重復(fù)周期，可獲取全天候SAR圖像，非常適合對地面沉降的災(zāi)害影響進(jìn)行監(jiān)測和評估。本文基于21景Sentinel-1A衛(wèi)星SAR影像采用MT-InSAR方法，獲取寧波地區(qū)2016—2018年期間高精度的地表沉降信息，包括年平均沉降速率、時序累計沉降量等，并與同期水準(zhǔn)進(jìn)行比較，驗證結(jié)果的精度和可靠性，分析寧波地區(qū)的沉降時空分布特征，并以地鐵2號線為例，分析軌道交通沿線地面沉降情況。

1 MT-InSAR方法

1.1 差分干涉處理

為了提高干涉測量成功幾率，確保形變測量的可靠性，選取使得時間基線(T)、空間基線(B)及頻率(f)參數(shù)最佳的一幅公共主影像，以確保得到整體最優(yōu)的基線參數(shù)(γm)。

(1)

本文選擇21幅SAR數(shù)據(jù)中的2017年5月10日影像作為主影像，與其他SAR影像進(jìn)行干涉處理，得到干涉圖。使用精密軌道數(shù)據(jù)(precise orbit ephemerides，POD)減弱軌道誤差的影響。然后利用SRTM DEM模擬地形相位并去除，生成時序上的差分干涉圖。

1.2 PS點選擇

PS干涉相位模型僅適用于在時間序列上保持相位穩(wěn)定的像素。因此，為了進(jìn)行PS干涉處理，首先要在SAR圖像上選擇保持相位穩(wěn)定的散射點。Ferretti等提出了利用SAR振幅離差指數(shù)來選取PS點的方法，它依據(jù)PS點的穩(wěn)定性可以通過回波相位在時間序列上一定的統(tǒng)計特性來表達(dá)[1,12]。振幅離差指數(shù)的表達(dá)式為

(2)

式中，g為目標(biāo)反射能量(正實數(shù))；σA為時序振幅標(biāo)準(zhǔn)差；σnR、σnI分別為噪聲實部和虛部標(biāo)準(zhǔn)差；σφ、μA分別表示時序相位標(biāo)準(zhǔn)差和振幅均值；DA為振幅離差指數(shù)。

選擇合適的振幅離差指數(shù)閾值(本文設(shè)為0.3)，大于該值的像素被認(rèn)為具有較穩(wěn)定的散射體，可作為PS候選點目標(biāo)。

利用Delaunay不規(guī)則三角網(wǎng)建立PS候選點間的關(guān)系，其相位梯度可表示為

ΔφAPS+Δφnoise

(3)

(4)

通過相位解纏可從參考點解算出格網(wǎng)上每點的形變速率和DEM誤差。

1.3 APS估計

針對大氣相位屏APS(atmospheric phase screen)的估計，采用時空濾波的方法[1,14-15]。該方法是基于殘余相位各分量在時空上具有不同的統(tǒng)計特性。由于大氣在空間上的一定尺度范圍內(nèi)(如1 km)是相關(guān)的，在空間域上表現(xiàn)為平滑的低頻信號，而大氣條件隨時間變化很快，在時間域上呈現(xiàn)高頻信號。因此，對殘余相位進(jìn)行時空濾波就可以估計出ΔφAPS。

2 研究結(jié)果與精度驗證

2.1 研究區(qū)與處理結(jié)果

本文以寧波地區(qū)為研究對象。數(shù)據(jù)選取21景升軌Sentinel-1A影像，入射角范圍為36.17°～36.22°，VV極化，時間范圍為2016年7月26日至2017年12月12日，影像參數(shù)見表1，時空基線參數(shù)如圖1所示。POD精密軌道數(shù)據(jù)由歐空局提供。

圖1 Sentinel-1A影像時空基線

序號衛(wèi)星獲取日期軌道方向極化方式傳感器模式入射角/(°)1Sentinel-1A2016-07-26AscendingVVIW36.172Sentinel-1A2016-10-06AscendingVVIW36.173Sentinel-1A2016-10-30AscendingVVIW36.174Sentinel-1A2016-11-23AscendingVVIW36.175Sentinel-1A2016-12-05AscendingVVIW36.176Sentinel-1A2016-12-17AscendingVVIW36.177Sentinel-1A2017-01-22AscendingVVIW36.178Sentinel-1A2017-02-15AscendingVVIW36.179Sentinel-1A2017-03-11AscendingVVIW36.1710Sentinel-1A2017-04-04AscendingVVIW36.2211Sentinel-1A2017-04-16AscendingVVIW36.2212Sentinel-1A2017-05-10AscendingVVIW36.2213Sentinel-1A2017-06-27AscendingVVIW36.2214Sentinel-1A2017-07-21AscendingVVIW36.2215Sentinel-1A2017-08-26AscendingVVIW36.2216Sentinel-1A2017-09-07AscendingVVIW36.2217Sentinel-1A2017-09-19AscendingVVIW36.2218Sentinel-1A2017-10-13AscendingVVIW36.2219Sentinel-1A2017-11-06AscendingVVIW36.2220Sentinel-1A2017-11-30AscendingVVIW36.2221Sentinel-1A2017-12-12AscendingVVIW36.22

基于Sentinel-1A影像采用MT-InSAR方法進(jìn)行處理后，共識別出370 469個PS點，每平方千米范圍內(nèi)包含約376個PS點。這種密度已經(jīng)可以滿足MT-InSAR技術(shù)反演地面沉降和沉降分析的要求。年沉降速率和累計沉降量是地面沉降監(jiān)控和評價中常用的指標(biāo)，能夠正確定量地反映研究區(qū)的沉降情況。得到研究區(qū)地面沉降變化的空間分布特征如圖2 所示。

圖2 研究區(qū)PS點平均沉降速率和累計沉降量

2.2 精度驗證

本文利用水準(zhǔn)實測數(shù)據(jù)，對InSAR處理結(jié)果進(jìn)行精度評定。MT-InSAR結(jié)果以水準(zhǔn)測量為基準(zhǔn)，并轉(zhuǎn)化為高程方向，統(tǒng)一基準(zhǔn)后再進(jìn)行比較。本文選取了寧波市水準(zhǔn)觀測網(wǎng)中的7個水準(zhǔn)點與鄰近InSAR監(jiān)測點(<100 m)進(jìn)行比較分析。InSAR沉降速率和水準(zhǔn)沉降速率比較見表2。二者差值絕對值最大為2.7 mm/a，最小值為0.6 mm/a，均方根誤差RMSE為1.7 mm/a。由此可見，基于Sentinel-1A影像采用MT-InSAR技術(shù)獲取的寧波地區(qū)地面沉降結(jié)果與水準(zhǔn)觀測十分吻合，驗證了研究結(jié)果的精度和可靠性。

表2 InSAR反演和水準(zhǔn)觀測的平均沉降速率比較 mm/a

3 沉降時空特征分析

3.1 地面沉降分布特征

如圖2所示，研究結(jié)果表明2016—2018年寧波地面沉降現(xiàn)象較為嚴(yán)重，沉降中心分布不集中且數(shù)量較多，廣泛分布在海曙區(qū)、鎮(zhèn)海區(qū)、江北區(qū)、鄞州區(qū)及北侖區(qū)。

位于海曙區(qū)的沉降漏斗連接成片，最大沉降速率為22.5 mm/a，沉降量達(dá)30.3 mm，是研究區(qū)域內(nèi)最嚴(yán)重的沉降中心。此外，海曙區(qū)的櫟社周邊等沉降漏斗也連接成片，最大沉降速率20.9 mm/a，最大沉降量28.2 mm。鎮(zhèn)海區(qū)存在3個沉降中心，分別位于鎮(zhèn)海煉化附近，最大沉降速率20.2 mm/a，最大沉降量27.3 mm；寧波幫博物館附近，最大沉降速率19.1 mm/a，最大沉降量25.8 mm；敬德小學(xué)附近，最大沉降速率21.1 mm/a，最大沉降量28.6 mm。江北區(qū)存在3個沉降中心，分別位于宏圖路與通寧路交叉口附近，最大沉降速率17.7 mm/a，最大沉降量24.0 mm；機(jī)場北路附近，最大沉降速率21.5 mm/a，最大沉降量29.0 mm；沈海高速寧波北樞紐附近，最大沉降速率16.3 mm/a，最大沉降量24.4 mm。鄞州區(qū)存在3個沉降中心，分別位于姜山鎮(zhèn)周邊，最大沉降速率20.5 mm/a，最大沉降量27.7 mm；寧波東站附近，最大沉降速率19.9 mm/a，最大沉降量27.0 mm；鄞州邱隘實驗小學(xué)以南，最大沉降速率17.3 mm/a，最大沉降量23.4 mm。北侖區(qū)僅存在一個沉降中心，位于北侖區(qū)中醫(yī)院附近，最大沉降速率20.3 mm/a，最大沉降量27.4 mm。

3.2 軌道交通沿線地面沉降

如圖3所示，主要交通公路、鐵路均有部分落入沉降區(qū)內(nèi)。如沈海高速公路穿過本研究范圍內(nèi)的一個沉降漏斗及兩個沉降帶，分別位于鄞州區(qū)的姜山鎮(zhèn)周邊，海曙區(qū)櫟社周邊等沉降漏斗連接成片的區(qū)域，古林鎮(zhèn)、高橋鎮(zhèn)等沉降漏斗連接成片的區(qū)域。北侖鐵路穿過本研究范圍內(nèi)的兩個沉降漏斗，分別位于鄞州區(qū)的寧波東站附近、鄞州邱隘實驗小學(xué)以南。以寧波地鐵2號線為例，該線路由櫟社國際機(jī)場站至清水浦站，已營運線路經(jīng)過海曙區(qū)、江北區(qū)及鎮(zhèn)海區(qū)。根據(jù)Sentinel-1A影像的InSAR反演結(jié)果，寧波地鐵2號線沿線地面沉降空間分布特征如圖3所示。2016—2018年間，地鐵2號線沿線的主要沉降中心集中在海曙區(qū)及鎮(zhèn)海區(qū)，其中最大沉降速率19.9 mm/a，最大沉降量26.9 mm，位于鎮(zhèn)海區(qū)。沉降較為嚴(yán)重的地鐵站點有路林站、三官堂站、櫟社國際機(jī)場站、櫟社站等。

圖3 寧波地鐵2號線沿線PS點平均沉降速率

4 結(jié) 語

本文基于21景Sentinel-1A影像采用MT-InSAR方法研究了寧波地區(qū)地面沉降，利用研究區(qū)同期水準(zhǔn)結(jié)果檢驗了本文的研究結(jié)果，兩者的差值均小于3.0 mm/a，RMSE為1.7 mm/a。說明MT-InSAR方法反演結(jié)果與水準(zhǔn)獲取結(jié)果吻合，驗證了MT-InSAR技術(shù)在寧波地區(qū)沉降監(jiān)測中的精度和可靠性，可以代替?zhèn)鹘y(tǒng)水準(zhǔn)作為一種沉降監(jiān)測的手段。

本文分析了2016—2018年期間的寧波地面沉降變化的時空分布特征。研究區(qū)域最大沉降速率為22.5 mm/a，沉降量達(dá)30.3 mm，位于海曙區(qū)。分析了軌道交通沿線地面沉降分布特征，主要交通公路、鐵路均有部分落入沉降區(qū)內(nèi)。本文研究結(jié)果可為城市公共安全及地質(zhì)災(zāi)害預(yù)防提供參考資料，為城市規(guī)劃發(fā)展、智慧城市建設(shè)及社會精細(xì)化治理提供科學(xué)依據(jù)。