基于InSAR的重慶東港地區(qū)地表形變監(jiān)測(cè)

李華蓉,戴雙璘

(重慶交通大學(xué) 智慧城市學(xué)院,重慶 400074)

大型機(jī)械長(zhǎng)期在高填方區(qū)域作業(yè)導(dǎo)致了重慶東港集裝箱碼頭區(qū)域地表形變較為嚴(yán)重,引發(fā)了東港區(qū)域支座墩柱出現(xiàn)裂縫以及軌道移位嚴(yán)重等現(xiàn)象,對(duì)工作人員的生命財(cái)產(chǎn)安全構(gòu)成了極大威脅。傳統(tǒng)地表形變監(jiān)測(cè)方法(如水準(zhǔn)測(cè)量監(jiān)測(cè)、GNSS監(jiān)測(cè)等)獲取形變信息精度高,但作業(yè)周期長(zhǎng),勞動(dòng)強(qiáng)度大,僅能獲取少量離散測(cè)點(diǎn)的形變信息,且測(cè)點(diǎn)分布稀疏、易被破壞、難以實(shí)時(shí)更新,具有一定局限性。

合成孔徑雷達(dá)干涉測(cè)量(interferometric synthetic aperture radar, InSAR)技術(shù)起源于20世紀(jì)60年代中后期,集成了合成孔徑雷達(dá)技術(shù)和微波干涉技術(shù)[1]。20世紀(jì)90年代,D-InSAR技術(shù)在InSAR技術(shù)基礎(chǔ)上發(fā)展而來(lái),但D-InSAR技術(shù)極易受大氣擾動(dòng)和時(shí)空失相干影響,無(wú)法應(yīng)用于高精度、變化緩慢的地表形變監(jiān)測(cè),只適用于地震等形變劇烈的地表形變監(jiān)測(cè)[2]。為克服上述缺點(diǎn),時(shí)序InSAR技術(shù)應(yīng)運(yùn)而生,時(shí)序InSAR技術(shù)主要以永久散射體干涉測(cè)量(PS-InSAR)技術(shù)和小基線集干涉測(cè)量(SBAS-InSAR)技術(shù)為代表[3-4],能實(shí)現(xiàn)低成本、高精度、高密度的大范圍地表形變監(jiān)測(cè)。

近年來(lái),學(xué)者們對(duì)時(shí)序InSAR技術(shù)進(jìn)行了深入研究。王義梅等[5]基于PS-InSAR和SBAS-InSAR技術(shù),處理了2012—2013年的15景TerraSAR-X影像,提取了鄭州市的沉降場(chǎng),并結(jié)合地質(zhì)條件及已有的研究成果進(jìn)行分析,表明時(shí)序InSAR技術(shù)在大范圍地表形變監(jiān)測(cè)應(yīng)用的有效性;楊國(guó)創(chuàng)等[6]聯(lián)合時(shí)序InSAR和GNSS技術(shù)對(duì)深圳和香港海岸線的地表形變進(jìn)行監(jiān)測(cè),將基于SBAS-InSAR獲得的矢量反演結(jié)果與基于PS-InSAR技術(shù)獲得的反演結(jié)果進(jìn)行比較分析,研究表明:這二者結(jié)果相近,形變結(jié)果均表明位于填海區(qū)的地鐵建筑及地基較淺的建筑物均存在明顯沉降現(xiàn)象;盧旺達(dá)等[7]利用PS-InSAR技術(shù)處理了2015—2018年的24景Sentinel-1A數(shù)據(jù),獲取了天津地區(qū)地面沉降結(jié)果,并利用SBAS-InSAR技術(shù)進(jìn)行驗(yàn)證,這兩種方法得到的累積形變量一致性高,表明該結(jié)果可為天津市災(zāi)害防治提供數(shù)據(jù)支持。上述研究均表明:時(shí)序InSAR技術(shù)在地表形變監(jiān)測(cè)領(lǐng)域均能取得較好效果。

基于此,筆者利用Sentinel-1數(shù)據(jù)和改進(jìn)SBAS-InSAR 技術(shù)對(duì)重慶東港集裝箱碼頭區(qū)域進(jìn)行地表形變監(jiān)測(cè),將基于PS-InSAR技術(shù)獲得的永久散射體(persistent scatter, PS)點(diǎn)作為地面控制點(diǎn),并引入SBAS-InSAR技術(shù)處理流程,與人工選取地面控制點(diǎn)的SBAS-InSAR技術(shù)進(jìn)行對(duì)比,并將兩種方法得到的形變結(jié)果與二等水準(zhǔn)測(cè)量數(shù)據(jù)進(jìn)行精度評(píng)定,探究改進(jìn)了SBAS-InSAR技術(shù)的可靠性。

1 數(shù)據(jù)集與研究方法

1.1 研究區(qū)概況與數(shù)據(jù)集

研究區(qū)域地處重慶市南岸區(qū),該區(qū)域地跨東經(jīng)106°44′—106°46′,北緯29°34′—29°36′,最低海拔為150 m,最高海拔為327 m,面積約為38.84 km2。研究區(qū)域的重慶東港集裝箱碼頭位于長(zhǎng)江江岸,地勢(shì)平坦,植被覆蓋較少。該區(qū)域?yàn)楦咛罘絽^(qū),長(zhǎng)期有大型機(jī)械在地面作業(yè),引發(fā)地面振動(dòng),導(dǎo)致填土不斷被壓實(shí),地表發(fā)生形變。在自然、人為因素的雙重作用下,該研究區(qū)在運(yùn)營(yíng)期間曾出現(xiàn)了纜車斜坡道蓋梁支座墩柱脫空現(xiàn)象,岸側(cè)蓋梁與支座墩柱中心線偏差較大,部分墩柱出現(xiàn)裂縫,靠近江側(cè)裂縫較多,軌道移位嚴(yán)重等。研究區(qū)地表形變較為嚴(yán)重,對(duì)纜車道運(yùn)行造成較大安全隱患,需定期對(duì)該研究區(qū)進(jìn)行地表形變監(jiān)測(cè)。研究區(qū)域及東港集裝箱碼頭區(qū)域如圖1。

圖1 研究區(qū)域Fig. 1 Research area

筆者選取31景Sentinel-1升軌、VV極化方式、單視復(fù)數(shù)、干涉寬幅模式的SAR衛(wèi)星影像,時(shí)間跨度為2018年1月4日—2019年12月25日,平均影像入射角為33°18′18″,每景影像時(shí)間間隔為24 d。此外,筆者還利用定位精度優(yōu)于5 cm的Sentine-1的衛(wèi)星精密軌道(precise orbit ephemerides, POD)數(shù)據(jù)修正軌道信息,去除因軌道誤差引起的系統(tǒng)性誤差。利用30 m分辨率的SRTM-1數(shù)據(jù)作為外部DEM數(shù)據(jù),去除平地相位和地形相位。

1.2 研究方法

在SBAS-InSAR處理流程中,須使用GCP文件修正干涉相位和解纏后的相位。GCP是被認(rèn)為形變?yōu)榱愕牡孛婵刂泣c(diǎn),GCP選取應(yīng)該遠(yuǎn)離形變區(qū)域、解纏錯(cuò)誤相位躍變區(qū)域,應(yīng)考慮平地區(qū)域及相干性較高區(qū)域。傳統(tǒng)SBAS-InSAR技術(shù)通常是操作人員根據(jù)GCP點(diǎn)的特性,手動(dòng)選取GCP點(diǎn),結(jié)果具有不穩(wěn)定性且依賴操作人員經(jīng)驗(yàn)判斷。

考慮PS-InSAR技術(shù)可從影像中提取出在各個(gè)時(shí)間都能保持高相干性的PS點(diǎn)。PS點(diǎn)常用的提取方法有振幅離差指數(shù)法和相干系數(shù)閾值法[8],振幅離差指數(shù)DA的計(jì)算如式(1)。

(1)

式中:μ為振幅標(biāo)準(zhǔn)差;σ為振幅均值。

通過(guò)振幅離差指數(shù)來(lái)判斷相位標(biāo)準(zhǔn)差大小,振幅離差指數(shù)越大,相位標(biāo)準(zhǔn)差越小[9]。相干系數(shù)閾值法是通過(guò)限制相干系數(shù)閾值對(duì)PS點(diǎn)進(jìn)行篩選,相干系數(shù)位于0～1,越接近1,相干性越好。

為提高基于SBAS-InSAR技術(shù)獲取的地表形變精度,減小人工選取GCP點(diǎn)帶來(lái)的誤差,筆者利用PS點(diǎn)選取特性與GCP選取特性一致,先對(duì)研究區(qū)基于PS-InSAR技術(shù)進(jìn)行處理,獲取穩(wěn)定PS點(diǎn),再通過(guò)限制振幅離差指數(shù)、相干性指數(shù)、年平均沉降速率這3個(gè)指標(biāo)閾值,從PS點(diǎn)中篩選出滿足條件的點(diǎn)作為GCP,引入SBAS-InSAR技術(shù)處理流程。改進(jìn)的SBAS-InSAR技術(shù)處理流程如圖2。

圖2 改進(jìn)的SBAS-InSAR技術(shù)處理流程Fig. 2 Process flow chart of the improved SBAS-InSAR technology

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

1)數(shù)據(jù)預(yù)處理。使用精密軌道數(shù)據(jù)對(duì)SLC數(shù)據(jù)軌道信息進(jìn)行修正;對(duì)原始SLC進(jìn)行多視處理,得到強(qiáng)度信息圖;基于地理坐標(biāo)系對(duì)強(qiáng)度圖進(jìn)行裁剪,提取研究區(qū)域。

2)生成連接圖。根據(jù)空間基線為臨界基線的2%,時(shí)間基線為120 d的原則建立干涉對(duì)連接圖。共生成117對(duì)干涉對(duì),其中最大時(shí)間基線為120 d,最大的空間基線為95.95 m,干涉對(duì)基線的時(shí)間-空間分布如圖3。圖3中:空心點(diǎn)代表主影像,實(shí)心點(diǎn)代表輔影像,線條代表干涉對(duì)。

圖3 干涉對(duì)基線的時(shí)間-空間分布Fig. 3 Time-spatial distribution of interference versus baseline

3)差分干涉處理。按時(shí)間基線和最小原則選取2018年5月28日的影像作為超級(jí)主影像[10],輔影像基于超級(jí)主影像進(jìn)行配準(zhǔn),由每組像對(duì)的復(fù)數(shù)影像共軛相乘生成干涉相位圖[11]。利用外部DEM去除干涉圖中的平地相位和地形相位,采用最小費(fèi)流法(minimum cost flow, MCF)進(jìn)行相位解纏處理,解纏相干系數(shù)閾值設(shè)置為0.2,以去除相干性較低的像元。對(duì)該數(shù)據(jù)進(jìn)行距離向?yàn)?方位向?yàn)?的多視,獲得20 m制圖分辨率的InSAR產(chǎn)品,圖4為20180128、20180404像對(duì)的干涉結(jié)果。對(duì)生成的117對(duì)干涉結(jié)果進(jìn)行人工檢視,刪除相干性較低、解纏效果較差的干涉對(duì)。

圖4 20180128、 20180404像對(duì)的干涉結(jié)果Fig. 4 The result of interferometric graph between 20180128 and 20180404

4)進(jìn)行軌道精煉和重去平。軌道精煉與重去平的目的是估算和去除解纏后的相位圖中依然存在的殘余相位[12]。筆者采用改進(jìn)的SBAS-InSAR技術(shù),基于PS-InSAR技術(shù)獲取研究區(qū)內(nèi)的PS點(diǎn),通過(guò)三閾值判別法對(duì)PS點(diǎn)進(jìn)行篩選,選取振幅離差指數(shù)大于等于3.2、相干性指數(shù)大于0.75、年平均沉降速率位于[-1,1]之間的PS點(diǎn),共計(jì)28個(gè)。利用該GCP 點(diǎn)對(duì)基線參數(shù)進(jìn)行重新定義,計(jì)算出相位偏移,修改相位解纏圖像頭文件中的軌道參數(shù)。基于PS-InSAR技術(shù)獲取的GCP分布情況如圖5。

圖5 基于PS-InSAR技術(shù)獲取的GCP分布Fig. 5 GCP distribution map obtained by PS-InSAR technology

5)進(jìn)行SBAS反演估算與地理編碼。SBAS反演估算分為兩次。第一次反演是對(duì)軌道精煉和重去平后的干涉圖進(jìn)行殘余地形和形變相位的估算。第二次反演是通過(guò)在時(shí)間域上進(jìn)行高通濾波和在空間域上進(jìn)行低通濾波計(jì)算,計(jì)算和去除大氣相位,再使用上述GCP剔除殘余相位。此時(shí)獲取的形變結(jié)果在雷達(dá)坐標(biāo)系下,需將反演估算形變結(jié)果轉(zhuǎn)換到WGS-84坐標(biāo)系下。基于改進(jìn)SBAS-InSAR技術(shù)獲取的地表形變?nèi)鐖D6。

圖6 改進(jìn)SBAS-InSAR技術(shù)獲取的形變Fig. 6 Deformation map obtained by the improved SBAS-InSAR technology

3 精度評(píng)定

筆者對(duì)研究區(qū)進(jìn)行了同時(shí)段的二等水準(zhǔn)監(jiān)測(cè),將獲得的該區(qū)域形變數(shù)據(jù)作為真值,對(duì)上述基于SBAS-InSAR技術(shù)反演的形變數(shù)據(jù)進(jìn)行精度評(píng)定。同時(shí),與手動(dòng)選取地面控制點(diǎn)反演出的變形數(shù)據(jù)進(jìn)行對(duì)比研究,驗(yàn)證上述方法的可靠性和先進(jìn)性。

3.1 二等水準(zhǔn)監(jiān)測(cè)結(jié)果

本研究對(duì)A區(qū)域(東港集裝箱碼頭)進(jìn)行二等水準(zhǔn)測(cè)量,獲取了2018年1月—2019年12月間的形變數(shù)據(jù)。共布設(shè)48個(gè)水準(zhǔn)點(diǎn),具體數(shù)據(jù)如表1。需注意的是:水準(zhǔn)測(cè)量數(shù)據(jù)是基于重慶獨(dú)立坐標(biāo)系獲得的,因此需要在ArcGIS中利用特征點(diǎn)對(duì)水準(zhǔn)點(diǎn)進(jìn)行空間校正將其轉(zhuǎn)換到WGS 84坐標(biāo)系下,隨后才能與SBAS-InSAR技術(shù)反演的形變數(shù)據(jù)進(jìn)行精度評(píng)定,空間校正后的水準(zhǔn)點(diǎn)如圖7。

表1 水準(zhǔn)測(cè)量獲取的形變數(shù)據(jù)Table 1 Deformation data obtained by leveling measurement

圖7 空間校正后的水準(zhǔn)點(diǎn)Fig. 7 The levelling point after spatial correction

3.2 手動(dòng)選取GCP方法的反演結(jié)果

為探究基于PS-InSAR技術(shù)選取GCP方法的可靠性,筆者同時(shí)結(jié)合相位解纏圖、相干系數(shù)圖和Google earth歷史影像,根據(jù)GCP選點(diǎn)原則手動(dòng)選取28個(gè)GCP,分布情況如圖8。與圖5相比,手動(dòng)選取的GCP位置更加集中,均位于公路和建筑物范圍,而基于PS-InSAR技術(shù)獲取的GCP分布大部分位于建筑物,少部分位于裸地,總體分布更加均勻。

圖8 手動(dòng)選取的GCP分布Fig. 8 Manually selected GCP distribution map

基于手動(dòng)選取的GCP反演研究區(qū)地表形變的結(jié)果如圖9。與圖6相比,未改進(jìn)SBAS-InSAR技術(shù)獲取的形變結(jié)果最大形變量和最小形變量相差無(wú)幾,但4個(gè)形變嚴(yán)重區(qū)域的平均沉降量較小。

圖9 未改進(jìn)SBAS-InSAR技術(shù)獲取的形變Fig. 9 Deformation map obtained by the unimproved SBAS-InSAR technology

由于InSAR技術(shù)獲取形變量為雷達(dá)視線方向(light of sight, LOS)的一維形變量,即地表各方向形變量在雷達(dá)視線方向上的投影。需先進(jìn)行形變場(chǎng)解算,獲取垂直方向形變量,再與水準(zhǔn)測(cè)量獲取的形變量進(jìn)行精度評(píng)定[13]。根據(jù)文獻(xiàn)[14],列出視線向形變量垂直向轉(zhuǎn)換計(jì)算,如式(2)。

d=Δr/cosθ

(2)

式中:d為垂直向形變;Δr為視線向形變;θ為雷達(dá)入射角。

SBAS-InSAR技術(shù)獲取的垂直形變量如表2。

表2 改進(jìn)前后獲取的垂直形變數(shù)據(jù)Table 2 Vertical deformation data obtained before and after improvement

3.3 對(duì)比分析

將基于改進(jìn)SBAS-InSAR技術(shù)獲取的垂直形變和基于未改進(jìn)SBAS-InSAR技術(shù)獲取的垂直形變分別與水準(zhǔn)測(cè)量數(shù)據(jù)進(jìn)行對(duì)比。采用平均絕對(duì)誤差EMA和均方根誤差ERMS作為精度評(píng)定指標(biāo)[15],EMA和ERMS的計(jì)算如式(3)、式(4):

(3)

(4)

式中:d為SBAS-InSAR技術(shù)反演的垂直形變量;d1為水準(zhǔn)測(cè)量獲得的形變數(shù)據(jù)。

利用S-W test對(duì)絕對(duì)誤差進(jìn)行正態(tài)性檢驗(yàn)。精度評(píng)定表、S-W檢驗(yàn)結(jié)果和頻率分布直方圖如表3、表4和圖10。

表3 精度評(píng)定Table 3 Accuracy evaluation

表4 S-W檢驗(yàn)結(jié)果Table 4 S-W test results

圖10 頻率分布直方圖Fig. 10 Frequency distribution histogram

4 地表形變結(jié)果與分析

筆者分別基于PS-InSAR技術(shù)篩選GCP和手動(dòng)選取GCP的SBAS-InSAR技術(shù)獲取研究區(qū)地表形變場(chǎng)。由結(jié)果可知:研究區(qū)形變量約為-60～30 mm,研究區(qū)共存在4個(gè)形變嚴(yán)重區(qū)域,其中:A區(qū)域?yàn)橹貞c東港集裝箱碼頭,B區(qū)域?yàn)橹貞c鐵路樞紐東環(huán)線拌和站,C區(qū)域?yàn)闁|港家園,D區(qū)域?yàn)橹貞c永翔現(xiàn)代物流產(chǎn)業(yè)園。A區(qū)域地表形變情況復(fù)雜,最大形變量達(dá)-52 mm,最大形變速率達(dá)-28 mm/a-1;C區(qū)域?yàn)榄h(huán)狀形變,外環(huán)形變約為-15 mm,內(nèi)部較為穩(wěn)定,幾乎沒(méi)有形變;B、 D區(qū)域呈漏斗狀形變場(chǎng),最大形變值為-30 mm。結(jié)合Google Earth歷史影像(圖11)和相關(guān)資料可知:A區(qū)域?yàn)楦咛罘絽^(qū)域,大型機(jī)械作業(yè)及貨物運(yùn)輸使填土受壓不均,導(dǎo)致地表沉降現(xiàn)象嚴(yán)重。2018—2019年,B、D區(qū)域正處于施工期,實(shí)際形變情況與漏斗狀形變結(jié)果一致;C區(qū)域的東港家園周邊處于施工期,與環(huán)狀形變結(jié)果相符。

圖11 Google Earth歷史影像Fig. 11 Google Earth historical imagery

由精度評(píng)定結(jié)果可知:改進(jìn)后SBAS-InSAR技術(shù)的EMA和ERMS分別為6.27、7.85 mm,未改進(jìn)SBAS-InSAR技術(shù)的EMA和ERMS分別為6.67、 8.57 mm。這兩種方法的ERMS都滿足文獻(xiàn)[14]要求,改進(jìn)后SBAS-InSAR技術(shù)的EMA和ERM更小。利用S-W檢驗(yàn)法對(duì)兩種方法的EMA進(jìn)行正態(tài)性檢驗(yàn),由檢驗(yàn)結(jié)果可知:這兩種方法的EMA均滿足正態(tài)分布,其中,改進(jìn)后SBAS-InSAR技術(shù)的EMA的S-W檢驗(yàn)P值更大,更符合正態(tài)分布。

5 結(jié) 語(yǔ)

筆者分別將基于PS-InSAR技術(shù)篩選的GCP和手動(dòng)選取的GCP引入SBAS-InSAR技術(shù)處理流程,對(duì)兩種方法反演地表形變進(jìn)行結(jié)果分析和精度評(píng)定。研究結(jié)果表明:SBAS-InSAR技術(shù)可作為一種有力的新型監(jiān)測(cè)手段進(jìn)行地表形變監(jiān)測(cè),利用PS-InSAR技術(shù)篩選GCP方法比手動(dòng)篩選GCP方法更加自動(dòng)化,獲取的形變結(jié)果精度更高。

但本研究仍存在不足之處:如缺乏坐標(biāo)轉(zhuǎn)換參數(shù),無(wú)法將水準(zhǔn)點(diǎn)準(zhǔn)確地轉(zhuǎn)換至WGS-84坐標(biāo)系,在空間上造成誤差;SAR影像覆蓋時(shí)間段與水準(zhǔn)測(cè)量時(shí)間段不能完全一致,在時(shí)間上造成誤差。