李紅梅,陶秋香

(1.青島濱海學院,青島 266510;2.山東科技大學,山東省基礎地理信息與數(shù)字化技術重點實驗室,青島 266590)

?

基于InSAR時序分析的北京地區(qū)地表沉降監(jiān)測研究

李紅梅1,陶秋香2

(1.青島濱海學院,青島 266510;2.山東科技大學,山東省基礎地理信息與數(shù)字化技術重點實驗室,青島 266590)

針對常規(guī)大地測量技術在城市地表沉降監(jiān)測過程中時空分辨率較低的的局限性,本文利用PS和SBAS兩種InSAR時序分析技術對城市地表形變進行監(jiān)測。論文首先對SBAS-InSAR技術的原理進行詳細推導和分析,在此基礎上,對覆蓋北京地區(qū)的29景ENVISAT-ASAR數(shù)據(jù)采用PS和SBAS技術進行處理,分別提取了該區(qū)域地表年平均形變速率圖。通過與已發(fā)表成果對比發(fā)現(xiàn),PS和SBAS在大范圍、持續(xù)緩慢城市地表形變監(jiān)測中能夠探測一致的地表形變運動趨勢;并且在數(shù)據(jù)量有限的情況下SBAS技術能夠獲得更為可靠的監(jiān)測結果。

PS;SBAS;InSAR時序分析;城市地表形變監(jiān)測

0　引　言

隨著國民經(jīng)濟的迅速發(fā)展,城市地下水開采量也急劇增加,由此引發(fā)了一系列的城市地表沉降問題,嚴重制約著經(jīng)濟與環(huán)境的協(xié)調(diào)和可持續(xù)發(fā)展。北京地區(qū)擁有2000萬人口的特大型都市,城區(qū)人口密集,有關資料顯示該地區(qū)城市地表沉降越來越嚴重。因此,北京地區(qū)開展城市地面沉降監(jiān)測與防范相關的技術研究具有重要意義。

合成孔徑雷達差分干涉測量技術(D-InSAR)[1]以其監(jiān)測周期短、覆蓋范圍廣以及成本低的優(yōu)勢被廣泛應用于較大區(qū)域的城市地表沉降監(jiān)測中。常規(guī)D-InSAR技術側重研究短時間間隔的單次形變,無法獲取研究區(qū)域地面沉降隨時間的演化情況,而且其測量精度易受時空失相關和大氣延遲等因素的影響[2]。為了克服這些困難,米蘭工業(yè)大學的Ferretti與Prati等人基于干涉點目標分析的思想,提出了永久散射體(PS)法[3],并成功應用在Ancona地區(qū)的滑坡監(jiān)測中。國內(nèi)對PS-InSAR技術應用在城市地表形變監(jiān)測中的研宄工作起步較晚,最早由西南交通大學陳強、劉國祥等人于2006年在上海地區(qū)成功的利用PS技術提取了上海地面沉降場。隨著研究的不斷深入,學者們發(fā)現(xiàn)PS-InSAR技術監(jiān)測精度比較依賴于主影像的質(zhì)量,并且基線過長容易導致幾何失相干,因此該方法往往受到數(shù)據(jù)量較少的限制。為此Berardino等人于2002年提出了短基線集(SBAS)方法[4],該方法選擇基線較短的干涉對組成不同集合,集合之間基線較長,采用最小二乘法得到每個基線集合的地表形變時間序列,再利用奇異值分解法將多個集合聯(lián)合起來求解。不要求具有同一主影像,因此可以充分利用SAR數(shù)據(jù)生成較多的干涉圖,實現(xiàn)有限數(shù)據(jù)量條件下的高精度反演。

本文在對PS和SBAS技術進行分析的基礎上,以北京地區(qū)為試驗區(qū),對覆蓋該試驗區(qū)2007～2010年的29景ENVISAT-ASAR數(shù)據(jù)進行處理,獲取該區(qū)域地表年平均形變速率圖,分析了引起北京沉降的主要原因。通過對比PS和SBAS所得的結果進行表明, 兩種時序分析方法在大范圍、持續(xù)緩慢城市地表形變的監(jiān)測中可以獲取一致的地表形變運動趨勢;通過與已發(fā)表的成果對比可以發(fā)現(xiàn),在數(shù)據(jù)量有限的情況下SBAS技術監(jiān)測的形變結果更可靠。

1　InSAR時序分析技術原理

1.1PS-InSAR技術原理

2000年,Ferretti等[5]為克服常規(guī)D-InSAR技術受大氣延遲、時空去相關和DEM誤差等因素的限制,提出了PS-InSAR技術。該技術利用二維線性(或非線性) 相位回歸分析模型,對長時間間隔內(nèi)依然保持高相干的永久散射體(PS)的差分干涉相位進行回歸分析來獲取形變信息。

現(xiàn)有N+1景SAR影像。其中一景影像為公共主影像,另外N景影像為輔影像,將N景輔影像分別與公共主影像進行配準并做干涉處理。利用外部DEM做干涉處理后,得到N干涉圖,其干涉相位為

φdiff=φdef+φtopo+φatm+φflat+φnoise,

(1)

式中: φdef為兩次成像期間內(nèi)的地表形變產(chǎn)生的相位; φtopo表示由外部DEM誤差引起的地形殘差相位; φflat表示由衛(wèi)星軌道誤差導致的殘余平地相位; φatm為兩次成像時雷達信號穿過不均勻大氣層時引起的延遲相位; φnoise為噪聲相位。式中前兩項是主要相位貢獻值,后面三項為殘余相位,表示為δ φres,則:

φdiff=φdef+φtopo+φres

(2)

假設Δd=v·T,其中v為形變速率,則有

φdef=k1v·T+k2·Δherror+φres,

(3)

其中:T為兩次成像時間間隔;φres為殘余相位; Δherror為高程模糊數(shù);k1,k2為系數(shù)。

從上述公式可以看出,相位差分模型為二維線性相位模型,通過多次回歸分析可以得到地表形變的線性速率和殘余DEM誤差及其各類噪聲相位。

1.2SBAS技術原理

2001年,Berardino和Lanari等[6]提出了小基線集(SBAS)InSAR技術。該技術主要解決PS技術中時間上采樣過于稀疏的問題[7]。

設研究區(qū)有N+1幅SAR影像,獲取時間依次為t0,t1,…tn.根據(jù)時空基線條件選擇其中一幅影像作為公共主影像,其他N幅影像作為輔影像,得到M個干涉像對。假設在tA和tB時刻獲得的兩景SAR影像進行差分干涉處理后生成第j個差分干涉對,則差分干涉相位δφj可以表示為

δφj(r,x)=φ(tB)-φ(tA)

(4)

式中:d(tA)和d(tB)表示tA和tB時刻相對于參考時刻t0(假設d(t0)≡0)的累積形變信息;φ(tA)和φ(tB)表示相應相位值。

現(xiàn)用向量φT=[φ(t1),…,φ(tN)]表示N個相位圖,向量δφT=[δφ1,…,δφM]表示M個干涉相位圖,假設主圖像獲取時間先于輔圖像獲取時間,則:

δφj=φtISj-φtIMj,j=1,…,M,

(5)

式中:IM表示干涉像對中主圖像;IS表示干涉像對中輔圖像。

定義一個含有M個方程,N個未知參數(shù)的方程組,其表達式表示為

δφ=Aφ,

(6)

式中: 矩陣A[M×N]每一行對應一個干涉對,A[j,ISj]=1和A[j,IMj]=-1,其余元素為零。矩陣A是取決于干涉圖組合方式的一個近似關聯(lián)矩陣。當所有干涉對屬于同一個基線集時,A是一個列滿秩矩陣。當M=N時,方程組(9)有唯一解; 當M>N時,方程組是一個超定方程,此時在最小二乘約束下可以求得其唯一解

(7)

當所有干涉對屬于不同的基線集時(設有L(L>1)個基線集)，A的秩為N-L+1,秩虧,為此,采用奇異值分解的方法計算矩陣A的廣義逆矩陣A+,進而求解方程組(6)的最小二乘最小范數(shù)解,其過程為

對矩陣A進行奇異值分解

A+=US+VT,

(8)

式中,U為M×N的正交矩陣,其前N列是ATA特征向量,稱為矩陣A的左奇異矩陣,V為N×N的正交矩陣,它的所有列ATA是的特征向量,稱為矩陣A的右奇異矩陣,S為M×N矩陣。則:

A+=VS+UT,

(9)

則方程組(6)的最小二乘最小范數(shù)解為

(10)

式中: ui和vi分別是各自矩陣U和V的列向量。

將方程組(6)中對相位的求解轉化為對平均相位變化速率的求解問題,則待求參數(shù)向量為

(11)

則有:

δφ=Bv.

(12)

式中: B為M×N的矩陣,其中B[i,j]=tj+1-tj, (ISi+1≤j≤IMi,i=1,…,M),其它元素值為零。對B進行奇異值分解,解出各時間段平均速度v.

2　實驗分析

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

本文研究區(qū)如圖1所示,覆蓋整個北京地區(qū),中心地理坐標為北緯40.02°、東經(jīng)116.25°,覆蓋面積約為4 200km2.數(shù)據(jù)選用29景Track490的ENVISATASAR降軌數(shù)據(jù),時間跨度為2007年1月至2010年10月。實驗過程中還使用了荷蘭Delft大學發(fā)布的DORIS精軌數(shù)據(jù)和SRTM的 90m分辨率的外部DEM.

圖1　研究區(qū)域及周邊環(huán)境示意圖

2.2實驗處理與結果

1) 數(shù)據(jù)預處理

試驗收集的SAR影像是原始雷達數(shù)據(jù),進行時序分析之前需要對數(shù)據(jù)進行預處理,生成時序分析可以處理的干涉圖。主要包括原始數(shù)據(jù)成像、圖像配準、生成干涉圖主要步驟。干涉圖是數(shù)據(jù)處理的基礎,生成干涉圖的前提是兩景影像的精確配準,數(shù)據(jù)處理中采用相干系數(shù)法[8]逐一像素進行配準,保證配準精度。

2)PS時序分析

PS-處理過程中主要包括以下幾個關鍵環(huán)節(jié),分別是:PS點選擇、相位解纏、軌道誤差去除、地形相位去除,DEM誤差相位的估算與去除、大氣誤差相位的估算與去除等。

其中,永久散射特性的像素點(PS點)是PS時序分析法的處理對象。實驗采用Hooper博士提出的方法對PS點進行選取。研究區(qū)覆蓋城區(qū),建筑物較為密集,干涉效果比較好。試驗共選取了302 387個點,其密度大約為34個/km2.

確定PS點后,采用基于Delaunay不規(guī)則三角網(wǎng)法對纏繞的相位進行解纏[9]。解纏的方式是在選取的PS點上構建三角格網(wǎng),對三角格網(wǎng)按照時間域上的相位差和空間域上的相位差的順序逐步進行解纏。根據(jù)公式(1)顯示,每一個PS點的相位值包含地表形變相位、DEM誤差相位、大氣延遲和軌道誤差相位。為了提取地表形變相位值,需對其他相位成分加以估計并剔除。

在選取PS點時,計算出DEM誤差與垂直基線成比例關系,通過計算比例值大小判斷影像中對應的像素是否足夠穩(wěn)定,進而考慮其是否可以確定為PS點。在選取PS點的同時也估算出了DEM誤差相位。試驗中使用分辨率為90m、高程精度為16m的SRTM-DEM數(shù)據(jù)。為了消除大氣效應對數(shù)據(jù)處理的影響,2001年Ferretti假設大氣相位和軌道誤差相位是空間相關的,并且在時間域上和空間域上對解纏相位分別進行高通濾波和低通濾波來估算個影像的大氣效應誤差。圖2示出了PS處理過程中估計的各類誤差相位示例。

圖2　DEM誤差相位、大氣延遲和軌道誤差相位圖示例：(a) DEM誤差相位；(b) 主影像大氣延遲和軌道誤差相位

3)SBAS時序分析

與PS時序分析數(shù)據(jù)處理不同,SBAS時序分析的首要步驟是選取合適的差分干涉像對。實驗中選取了滿足空間基線小于800m、時間基線小于692天、相關系數(shù)大于0.65三個條件的110個干涉像對,組合方式如圖3所示。將110景干涉像對使用相干系數(shù)法進行精確配準,并做復共軛相乘得到原始干涉圖。

圖3　小基線集干涉圖組合方式

SBAS時序分析技術中選取的相位穩(wěn)定的點稱為緩慢失相關濾波相位像元點(SDFP),SDFP像元點是基于相位特征來定義的。試驗中共選取了290 364個SDFP點,其密度大約為33個/km2.

確定出SDFP點后,實驗采用與PS技術相位解纏相同的方法進行相位解纏,獲取解纏后差分相位。建立如公式(9)的矩陣方程,并求解得到地表形變量和形變速率。圖4示出了SBAS-InSAR處理過程中估計的各類誤差相位的示例。

經(jīng)過以上一系列的時間序列分析處理,最終得到的研究區(qū)年平均雷達視線向(LOS) 形變速率場,如圖 5所示。其中得到的形變量都是相對于主圖像20 080 908里面的PS/SDFP點的形變量進行分析得到的。形變速率圖中,紅色點位(負值)表示背離衛(wèi)星的方向運動,即該點的運動趨勢是沿雷達視線向方向下降;藍色點位(正值)表示朝著衛(wèi)星的方向運動,即該點的運動趨勢是沿雷達視線向方向上升。為了更好的顯示研究區(qū)的沉降情況,分別將PS法與SBAS法得到的高相干點的年平均沉降速率文件導入ArcGIS,利用Kriging插值方法得到整個研究區(qū)的沉降圖如圖5所示。

圖4　SBAS-InSAR技術估算的DEM誤差相位、大氣延遲和軌道誤差相位圖示例(a) DEM誤差相位；(b) 主影像大氣延遲和軌道誤差相位

圖5　研究區(qū)2007～2010年期間的沉降圖(a) PS法; (b) SBAS法

2.3結果分析

如圖5(a)表示PS技術獲取的研究區(qū)地表年平均沉降速率圖,該區(qū)域年平均沉降速率為21.2～-51.6mm/yr.北京城區(qū)范圍內(nèi)呈現(xiàn)紅色,表現(xiàn)為地表下沉的運動趨勢;周邊區(qū)為藍色,表現(xiàn)為地表抬升的運動趨勢。圖5(b)表示SBAS技術獲取的研究區(qū)地表年平均沉降速率圖,該區(qū)域年平均沉降速率為31.7～-109.9mm/yr.同樣的該監(jiān)測結果顯示城區(qū)范圍內(nèi)呈現(xiàn)紅色,表現(xiàn)為地表下沉的運動趨勢;周邊區(qū)為藍色,表現(xiàn)為地表抬升的運動趨勢。兩種技術監(jiān)測出的該區(qū)域地表沉降保持一致。

仔細分析后發(fā)現(xiàn),該區(qū)域沉降由中部到東部逐步加重的原因有兩點:1) 該區(qū)域主要是沖積扇的中下部地區(qū),土質(zhì)較細,過度抽取地下水會導致地面沉降; 2) 隨著一些重工業(yè)的東移,地表構筑物的建設以及地下工程建設加深了該區(qū)域的地表沉降。

西部以及周邊區(qū)域地表略有抬升的原因是該區(qū)域相干性低,受噪聲污染嚴重。仔細分析發(fā)現(xiàn):1) 該區(qū)域植被覆蓋較密集,這一點造成了數(shù)據(jù)處理過程中選取PS點和SDFP點過程產(chǎn)生誤差,增加了后續(xù)處理過程中噪聲的產(chǎn)生。2) 該區(qū)域地處山區(qū),地形起伏比較大,地形相位估計需要高分辨率的DEM數(shù)據(jù),但本次試驗中由于各種原因僅收集到了90m分辨率的數(shù)據(jù),精度較低,這一點也增加了數(shù)據(jù)處理的噪聲。

2.4精度分析

由于種種原因,未能獲取研究區(qū)外部的水準數(shù)據(jù)或GPS數(shù)據(jù),因此本文將PS和SBAS技術所得結果與已發(fā)表的文獻所述結果進行對比,表1給出了兩種方法所得結果與文獻[10]的對比情況。

從表1可知,本文利用SBAS技術獲取的北京地區(qū)地表形變監(jiān)測結果相對于PS技術的結果在數(shù)值上更接近已發(fā)表文獻的監(jiān)測結果,主要有以下兩個原因:1)SAR數(shù)據(jù)量較少,PS-InSAR技術在識別PS點時對數(shù)據(jù)量有一定的要求,實驗中由于資料有限僅收集了29景數(shù)據(jù),影響了PS-InSAR的監(jiān)測結果;2)SBAS-InSAR技術將29景SAR數(shù)據(jù)組合成若干個集合,形成了109景時間和空間基線均較小的差分干涉圖,彌補了PS技術依賴于單一公共主影像質(zhì)量的缺陷,從而獲取了相對更可靠的監(jiān)測結果。

表1　兩種方法所得的北京地區(qū)年平均沉降速

3　結束語

本文利用PS和SBAS技術分別對覆蓋北京地區(qū)的29景ENVISAT ASAR數(shù)據(jù)進行了處理,獲取了該地區(qū)2007年至2010年期間的沉降時間序列圖,發(fā)現(xiàn)研究區(qū)的東部沉降比較嚴重,部分地區(qū)年平均沉降速率達到90 mm/a以上。同時分析了整個研究區(qū)的沉降情況及形成原因,得到造成嚴重沉降的主要原因為地下水的過量開采。PS和SBAS法所得年平均速率標準差以及與已發(fā)表文獻的對比均表明, SBAS技術相對于PS技術監(jiān)測地面沉降能夠得到更為可靠的結果,與實際情況比較相符。

[1] GABRIEL K, GOLDSTEIN M, ZEBKER H A. Mapping small elevation changes over large areas:differential radar interferometer[J].Journal of Geophysical Research, 1989, 94(B7): 9183-9191.

[2] 許文斌,李志偉,丁曉麗,等. 利用MERIS 水汽數(shù)據(jù)改正ASAR干涉圖中的大氣影響[J]. 地球物理學報,2010, 53(5): 1073-1084.

[3] FERRETTI A, ROCCA F, PRATI C. Permanent scatterers in SAR interferometry [C]//Proceedings of IGARSS 1999, Hamburg, Germany, 28 June-2 July, 1999.

[4] BERARDINO P, FORNARO G, LANARI R,etal. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J]. IEEE Transactions on Geosciences and Remote Sensing, 2002, 40(11), 2375-2383.

[5] FERRETTI A, PRATI C, ROCCA F. Permanent scatterers in SAR interferometry[J].Geoscience and Remote Sensing, 2001, 39(1):8-20.

[6] BERARDINO P, FORNARO G, LANARI R,etal. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].Geoscience and Remote Sensing, 2002, 40(11):2375-2383.

[7] LANARI R, MORA O, MANUNTA M,etal. A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms[J].Geoscience and Remote Sensing, 2004, 42(7):1377-1386

[8] 陳富龍,張紅,王超.高分辨率SAR影像同名點自動匹配技術[J].中國圖象圖形學報,2006,11(9):1276-1281.

[9] 吳濤,張紅,王超,等.多基線距DInSAR技術反演城市地表緩慢形變[J].科學通報,2008,53(15):1849-1857.

[10]李永生,張景發(fā),李振洪等.利用短基線集干涉測量時序分析方法監(jiān)測北京市地面沉降[J].武漢大學學報(信息科學版),2013,38(11):1374-1377.

Research of the land Deformation Monitoring in Beijing Based on the InSAR Tme Series Analysis Method

LI Hongmei1,TAO Qiuxiang2

(1.QingdaoBinhaiUniversity,Qingdao266510,China; 2.GeomaticsCollege,ShandongUniversityofScienceandTechnology,Qingdao266590,China)

Aiming at the limitations of low spatial and temporal resolution of the traditional geodetic measurements in the process of urban surface subsidence monitoring. Two InSAR time series analysis methods were introduced into the surface subsidence monitoring. The paper firstly derived and analysed the basic principles of SBAS-InSAR technique detailedly. On this basis, we processed and analysed 29 scenes of ENVISAT-ASAR images covering the area using the PS-InSAR and SBAS-InSAR technique, getting the annual average rate of deformation of the area. Compared with published resultsshows that PS and SBAS are of consistency in deformation tendency of the wide range and slowly urban surface subsidence and in the case of a limited quantity of the data SBAS technique can obtain a more reliable monitoring result.

PS; SBAS; InSAR time series analysis; urban surface subsidence monitoring

2015-12-11

山東省優(yōu)秀中青年科學家科研獎勵基金項目(編號:BS2013SF013)

P642.26

A

1008-9268(2016)03-0078-06

李紅梅(1978-),女,碩士,講師,研究方向為短時交通流預測及深基礎變形監(jiān)測預測。

陶秋香(1977-),女,博士,副教授,研究方向為微波遙感技術與應用、InSAR關鍵技術、數(shù)據(jù)處理等。

doi:10.13442/j.gnss.1008-9268.2016.03.016

聯(lián)系人： 陶秋香 E-mail: qiuxiangtao@163.com