范雪婷李明巨張大騫趙朝賀

(1.江蘇省基礎(chǔ)地理信息中心,江蘇 南京210013;2.江蘇地質(zhì)調(diào)查研究院,江蘇 南京210018)

0 引 言

地表沉降是一項重要的地理國情信息。持續(xù)的地表沉降會導(dǎo)致建(構(gòu))筑物開裂、路面塌陷等災(zāi)害發(fā)生,破壞市政基礎(chǔ)設(shè)施和生產(chǎn)生活設(shè)施,嚴(yán)重時易引發(fā)城市內(nèi)澇,給城市建設(shè)和人民群眾生產(chǎn)生活造成較大影響。目前,傳統(tǒng)的地表沉降監(jiān)測技術(shù)手段主要是GNSS測量及水準(zhǔn)測量技術(shù),這些方法雖然監(jiān)測精度高,但是存在監(jiān)測點空間密度低、投入成本大,無法提供監(jiān)測全區(qū)域的細(xì)節(jié)沉降信息等問題。近幾年來發(fā)展起來的時間序列InSAR(Interferometry Synthetic Aperture Radar)技術(shù)可以獲取大范圍、高密度、連續(xù)時間序列的地表沉降數(shù)據(jù),逐步被應(yīng)用于各大城市的區(qū)域地面沉降監(jiān)測。本文主要采用多主影像相干目標(biāo)小基線InSAR方法(MCTSB-InSAR)[1],該方法的主體思路是將干涉相位分析與處理焦點完全集中在那些具備穩(wěn)定散射特性的像元點集上,基于這些穩(wěn)定點目標(biāo)進(jìn)行差分干涉相位建模,繼而從干涉相位中分離出單一信號分量并獲得所需的形變信息[2]。本文以無錫市為例,探索MCTSB-InSAR方法在大區(qū)域地表沉降監(jiān)測中的應(yīng)用效果,獲取無錫市2012年2月-2016年1月間的地表形變時空分布,并對監(jiān)測結(jié)果進(jìn)行精度評定。

1 監(jiān)測方法

MCTSB-InSAR技術(shù)是在永久散射體(Persistent Scatterer Interferometry,PS-InSAR)和小基線集(Small Baseline Subset,SBAS)方法基礎(chǔ)上進(jìn)行改進(jìn)的一種時間序列InSAR技術(shù),該方法綜合了PS-InSAR和SBAS兩種方法的優(yōu)勢,克服了相應(yīng)缺點。它不像PS-InSAR方法選取單一SAR影像為主影像,而是考慮滿足小基線要求的SAR影像任意組合干涉像對,從而降低了對影像數(shù)量的要求。MCTSB-InSAR也無需像SBAS方法要首先完成每個差分干涉圖的相位解纏,再進(jìn)行時序分析,而是通過考察空間鄰近的兩個高相干穩(wěn)定點的二次差分,同時完成形變參數(shù)解算和相位解纏[2]。另外MCTSB-InSAR方法可以實現(xiàn)多幅影像,特別是不同軌道沉降提取結(jié)果的無縫拼接,在跨軌道大區(qū)域沉降監(jiān)測中具有較大優(yōu)勢。MCTSB-InSAR技術(shù)監(jiān)測地面沉降算法原理主要包括組合干涉像對、穩(wěn)定點目標(biāo)提取、建立形變模型、非線性形變量估計、多幅InSAR地表沉降結(jié)果拼接5個方面。

1.1 基于特征區(qū)域相干性組合干涉像對

傳統(tǒng)小基線干涉像對組合方法主要為時空基線閾值法,該方法生成差分干涉圖后需要人工篩選,對干涉相位質(zhì)量較差的干涉對進(jìn)行人工剔除,導(dǎo)致工作量非常繁重。本文基于區(qū)域相干特性組合干涉像對,首先選取特征區(qū)域,通常選擇包含建筑物較多的高相干城區(qū),特征區(qū)域范圍應(yīng)較小,以利于快速處理。經(jīng)反復(fù)試驗,一般為500×500像元較為適宜。然后對所有干涉圖的相干性進(jìn)行定量評估,將平均相干系數(shù)由高到低排序,按影像數(shù)2-3倍自動篩選出排序在前的干涉像對,最后生成整個研究區(qū)域的干涉圖。與傳統(tǒng)通過設(shè)置時空基線硬閾值方法相比,該方法可實現(xiàn)優(yōu)質(zhì)干涉像對的自動、完備提取。

1.2 三級閾值穩(wěn)定點目標(biāo)提取方法

單一閾值法用于穩(wěn)定點目標(biāo)的識別具有局限性。本文利用振幅、相干系數(shù)、振幅離差指數(shù)信息三級閾值探測法進(jìn)行穩(wěn)定點目標(biāo)提取[3]。首先設(shè)置絕對平均相干系數(shù)閾值,選取高相干點HCP1;設(shè)置相對平均相干系數(shù)閾值和平均幅度閾值T A,選取那些相干性相對較高,并且平均幅度又很大的高相干點HCP2;然后根據(jù)相干目標(biāo)的穩(wěn)定性,設(shè)置幅度離差指數(shù)閾值D T,選取穩(wěn)定點目標(biāo)HCP3;依據(jù)公式HCP=(HCP1∪HCP2)∩HCP3,將獲取的點目標(biāo)作為最終穩(wěn)定點目標(biāo)[4]。

1.3 基于區(qū)域Delaunay三角網(wǎng)的線性形變速率與高程誤差估計

1.3.1 Delaunay三角網(wǎng)建立

構(gòu)建合理的穩(wěn)定點目標(biāo)網(wǎng)絡(luò)是探測線性形變的關(guān)鍵步驟之一,MCTSB-InSAR技術(shù)利用局部Delaunay三角網(wǎng)連接方法[1],將整景影像分割成相互重疊的網(wǎng)格,在各網(wǎng)格中使用Delaunay三角網(wǎng)進(jìn)行連接。三角網(wǎng)中一條邊就代表相鄰點目標(biāo)構(gòu)成的一個點目標(biāo)對,由此構(gòu)建的三角網(wǎng)既簡潔,又增加了連接的冗余度,可有效避免整體網(wǎng)絡(luò)斷裂成不連通的子網(wǎng)絡(luò)。

1.3.2 線性形變速率與高程誤差參數(shù)計算

通過局部Delaunay三角網(wǎng)連接所有的點目標(biāo),考察第k個差分干涉圖上i、j兩頂點之間的差分干涉相位:

式(1)中,wrap{}表示相位纏繞算子;系數(shù)a1由該點的入射角、斜距以及雷達(dá)波長λ決定,a2=是垂直基線;Δh i,j為高程誤差之差;t k為干涉圖的時間基線;Δv i,j為兩點間的線性形變速率之差;殘余相位增量,為兩點間非線性相位、大氣相位和噪聲相位之差的和。

考慮三角網(wǎng)兩點間非線性形變相位差決定性較小,大氣影響為低頻信號,在空間上存在一個相關(guān)距離,因此設(shè)置三角網(wǎng)兩點間距在大氣影響一定距離內(nèi)(一般1～2km),大氣相位差可以忽略,以及各高相干點的單點噪聲很小,滿足在該條件下,利用整體相位的相關(guān)系數(shù)最大值進(jìn)行優(yōu)化處理,建立如下目標(biāo)函數(shù):

式(2)中,J表示單位復(fù)數(shù);γ為整體相位相干系數(shù),其最大值在[0,1]之間,γ值越大表示模型估計值對差分相位的擬合程度越好,Δv i,j和Δh i,j的估值越可靠。

待網(wǎng)絡(luò)上的邊全部實現(xiàn)最大化求解后,對模型相干系數(shù)進(jìn)行最大值檢驗計算,判斷哪些連接邊計算出來的結(jié)果是可靠的[5]。一般設(shè)定適當(dāng)閾值(本文設(shè)為0.7),將大于該閾值的模型相干系數(shù)最大值所對應(yīng)的邊作為可靠的邊。

1.3.3 三角網(wǎng)集成

以某一具有已知形變量和DEM誤差的高相干點為參考點,對篩選出的優(yōu)質(zhì)三角網(wǎng)連接邊,建立聯(lián)合大型方程組。對連接邊上解算出的速率和高程誤差增量進(jìn)行積分,獲得各相關(guān)點的形變速率和高程誤差的絕對值[6]。

1.4 非線性形變量估計

從原始差分干涉圖相位減去高程誤差和線性形變速率分量可獲得殘余相位,主要包含非線性形變,大氣影響及噪聲等相位。非線性形變速率由于空間上不相關(guān)性,因此為高頻信號,而在時間上為低頻信號;大氣相位在時間上具有不相關(guān)性,為高頻信號;噪聲相位體現(xiàn)為在時間和空間上都不相關(guān)的隨機(jī)高頻信號[7]。因此,可以利用殘余相位三個分量的不同頻率特征,采取適當(dāng)?shù)念l域濾波方法,將三者分離出來。將得到的非線性形變結(jié)果與線性形變結(jié)果疊加,從而得到時間序列的累計形變量[1]。

1.5 多幅大區(qū)域InSAR地表沉降結(jié)果拼接

大范圍監(jiān)測區(qū)域需要多幅影像才能實現(xiàn)全覆蓋,但是各景影像獲取時間、影像數(shù)目都不一致,得到的地表沉降結(jié)果具有差異性,因此實現(xiàn)多景影像沉降結(jié)果的無縫拼接是一大挑戰(zhàn)。MCTSB-InSAR方法從測量平差角度出發(fā),通過提取影像間重疊區(qū)域的同名點結(jié)果,采用最小二乘平差方法算出各影像結(jié)果的改正數(shù),實現(xiàn)多幅點目標(biāo)結(jié)果的一次性無縫拼接。

2 實驗及結(jié)果分析

2.1 實驗數(shù)據(jù)介紹

以無錫市作為實驗區(qū)。選取加拿大C波段的Radarsat-2寬幅雷達(dá)影像作為沉降監(jiān)測的數(shù)據(jù)源,其單景影像覆蓋面積為150 km×150 km,影像分辨率約為30 m,共需兩幅影像實現(xiàn)無錫市全覆蓋(圖1)。每幅各有25期影像,數(shù)據(jù)獲取時段為2012年2月—2016年1月。采用30 m分辨率SRTM DEM數(shù)據(jù)去除地形貢獻(xiàn)項。

圖1 無錫市Radarsat-2影像覆蓋范圍

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

借助中國測繪科學(xué)研究院張永紅課題組研發(fā)的InSAR地表形變監(jiān)測系統(tǒng)(GDEMSI)進(jìn)行數(shù)據(jù)預(yù)處理和時序SAR影像的地表形變反演工作。利用特征區(qū)域相干性法組合干涉像對,選取包含建(構(gòu))筑物較多的500×500像元大小的城鎮(zhèn)地區(qū)作為特征區(qū),共生成63個最佳干涉像對(圖2)。

圖2 小基線對時空分布圖

設(shè)置絕對平均相干系數(shù)閾值0.72,相對平均相干系數(shù)閾值0.65和平均幅度閾值1.2,幅度離差指數(shù)閾值0.3,根據(jù)三級閾值法探測高精度穩(wěn)定點目標(biāo)。在此基礎(chǔ)上,創(chuàng)建Delaunay三角網(wǎng)差分相位函數(shù)模型,求解獲得線性形變信息和DEM誤差。最后依據(jù)殘余相位分量的不同頻率特征分離出非線性形變相位,與線性形變疊加獲取該幅影像完整形變值。在獲取兩幅監(jiān)測點結(jié)果后,通過提取影像間重疊區(qū)域的同名點結(jié)果,采用最小二乘平差方法對兩幅點目標(biāo)結(jié)果進(jìn)行拼接,最終形成區(qū)域協(xié)調(diào)一致的沉降結(jié)果。

2.3 結(jié)果分析

經(jīng)處理,最終獲取無錫市489 781個點目標(biāo),平均每平方千米監(jiān)測點數(shù)達(dá)106個,可見監(jiān)測點分布密度之大。從無錫市2012年2月—2016年1月沉降速率分布,可以看出,在監(jiān)測時段內(nèi),無錫市的沉降相對較為集中,主要分布在江陰南部和錫山北部地區(qū),涉及江陰市云亭街道、周莊鎮(zhèn)、徐霞客鎮(zhèn)、長涇鎮(zhèn)、華士鎮(zhèn)以及錫山區(qū)東港鎮(zhèn)等。全市最大沉降速率為28 mm/a,四年最大累積沉降量達(dá)94.6 mm,位于江陰市周莊鎮(zhèn);沉降速率超過10 mm/a的面積達(dá)97.3 km2。筆者收集到的無錫市2007—2012年的InSAR監(jiān)測結(jié)果,對比發(fā)現(xiàn)無錫市沉降范圍略有擴(kuò)大,沉降速率有所趨緩。

無錫市是中國重要的紡織品制造基地和出口基地,至2016年末,全市共有規(guī)模以上紡織服裝企業(yè)772家,每年生成大量的紡織品。通過現(xiàn)場調(diào)研發(fā)現(xiàn),無錫沉降區(qū)存在大量的紡織企業(yè),由于紡織行業(yè)是用水較多的行業(yè)之一,據(jù)估算,平均工業(yè)產(chǎn)值每100元就消耗3 m3的水,需水壓力較大[9]。因此,長期過量開采地下水是無錫市地面沉降的主要原因。2000年,江蘇省通了《關(guān)于在蘇錫常地區(qū)限期禁止開采地下水的決定》,2000-2005年無錫市共封井1 100眼,區(qū)域地質(zhì)環(huán)境發(fā)生了顯著改善。但根據(jù)InSAR監(jiān)測結(jié)果,江陰南部和錫山北部地區(qū)地面沉降仍相對較為嚴(yán)重,推測深部地下水禁采以后,江陰的工廠紛紛轉(zhuǎn)為開采淺部地下水,導(dǎo)致沉降范圍擴(kuò)大。

圖3 無錫市年沉降速率分布圖(2012年2月—2016年1月)

3 精度評定

為了驗證InSAR監(jiān)測結(jié)果的精度,搜集了無錫市30個水準(zhǔn)點數(shù)據(jù),2012—2015年每年12月份施測一次。因InSAR監(jiān)測點和水準(zhǔn)點位置的不完全一致性,根據(jù)約定臨近點原則對初始獲取的水準(zhǔn)點進(jìn)行了篩選:即以實測點位為中心,如果在一定距離范圍內(nèi)(本文設(shè)為80 m)存在至少一個InSAR點目標(biāo),則選擇該實測點參與精度評定。實驗中利用19個水準(zhǔn)點數(shù)據(jù)對InSAR監(jiān)測得到的無錫市地表平均沉降速率進(jìn)行比對驗證分析[8](表1)。

表1 無錫市水準(zhǔn)與InSAR沉降速率對比/mm·a-1

從表1可以看出,二者差值絕對值最大值為4.7 mm/a,誤差均方差為1.9 mm/a。該誤差實質(zhì)上還包括兩種數(shù)據(jù)時間不一致性和空間不一致性導(dǎo)致的誤差,可見MCTSB-InSAR方法的理論測量精度應(yīng)高于該值。

為了進(jìn)一步分析InSAR時序沉降監(jiān)測結(jié)果的可靠性,選取4個典型水準(zhǔn)監(jiān)測點(CJDL33、CJ-123、CJ078、I-145),利用2012-2015年的時序?qū)崪y成果與InSAR監(jiān)測點的時序沉降量進(jìn)行對比分析(圖4)。水準(zhǔn)點CJDL33位于東港鎮(zhèn)沉降中心,CJ-123位于祝塘鎮(zhèn)沉降中心,兩個水準(zhǔn)點表現(xiàn)為不斷下沉;CJ078位于惠山區(qū)堰橋街道的京滬高速鐵路沿線附近,呈輕微抬升趨勢;I-145位于江陰青陽鎮(zhèn),水準(zhǔn)點和InSAR監(jiān)測點表現(xiàn)都較為穩(wěn)定,形變量在+5～-5 mm范圍內(nèi)波動。通過對比發(fā)現(xiàn),無論是位于沉降中心還是抬升區(qū)亦或穩(wěn)定區(qū),二者的監(jiān)測結(jié)果具有較高的一致性。

綜上分析,利用時間序列InSAR技術(shù)監(jiān)測得到的年平均形變速率和時序累積形變量具都有較高的監(jiān)測精度,實驗結(jié)果能夠反映了區(qū)域的基本沉降特征。

圖4 無錫市時序InSAR沉降量與水準(zhǔn)沉降量曲線對比圖

4 結(jié) 語

(1)基于區(qū)域相干特性組合干涉像對,實現(xiàn)了優(yōu)質(zhì)干涉像對的自動、完備提取,極大提高了干涉圖質(zhì)量與生成效率。

(2)本文既顧及相干點回波信號的高信噪特征,又兼顧相干點的穩(wěn)定性,采用振幅、相干系數(shù)、振幅離差指數(shù)信息三級閾值探測法實現(xiàn)了穩(wěn)定點目標(biāo)的高精度、大密度提取。

(3)本文以 Radarsat-2為數(shù)據(jù)源,利用MCTSB-InSAR方法實現(xiàn)了無錫市2012年2月-2016年1月的地表形變的沉降信息反演。結(jié)果表明,在該階段內(nèi),無錫市沉降主要發(fā)生在江陰南部和錫山北部地區(qū),涉及江陰市云亭街道、周莊鎮(zhèn)、徐霞客鎮(zhèn)、長涇鎮(zhèn)、華士鎮(zhèn)以及錫山區(qū)東港鎮(zhèn)等,最大沉降速率為28 mm/a,發(fā)生在江陰市周莊鎮(zhèn)。利用水準(zhǔn)實測數(shù)據(jù)開展精度評定,二者具有較高一致性,驗證了MCTSB-InSAR方法利用Radarsat-2寬幅雷達(dá)影像獲取大區(qū)域地表形變信息的可靠性。